Changeset 637

Timestamp:

09/27/09 00:58:02 (15 years ago)

Author:

smidl

Message:

new Dirichlet random walk + numerics in egiw

Location:

library/bdm/stat

Files:

• exp_family.cpp (modified) (4 diffs)
• exp_family.h (modified) (2 diffs)

library/bdm/stat/exp_family.cpp

@@ -81 +81 @@
                      - 0.5 * dimx * ( m * log2 + 0.5 * ( dimx - 1 ) * log2pi )  - lg;
 
-        bdm_assert_debug ( ( ( -nkG - nkW ) > -Inf ) && ( ( -nkG - nkW ) < Inf ), "ARX improper" );
+//      bdm_assert_debug ( ( ( -nkG - nkW ) > -Inf ) && ( ( -nkG - nkW ) < Inf ), "ARX improper" );
+if (-nkG - nkW==Inf){
+        cout << "??" <<endl;
+}
         return -nkG - nkW;
 }
@@ -260 +263 @@
 }
 
-ivec eEmp::resample ( RESAMPLING_METHOD method ) {
-        ivec ind = zeros_i ( n );
+void eEmp::resample ( ivec &ind, RESAMPLING_METHOD method ) {
+        ind = zeros_i ( n );
         ivec N_babies = zeros_i ( n );
         vec cumDist = cumsum ( w );
@@ -334 +337 @@
         }
 
-        // copy the internals according to ind
+// copy the internals according to ind
         for ( i = 0; i < n; i++ ) {
                 if ( ind ( i ) != i ) {
@@ -341 +344 @@
                 w ( i ) = 1.0 / n;
         }
-
-        return ind;
 }
 

library/bdm/stat/exp_family.h

@@ -344 +344 @@
 UIREGISTER(eDirich);
 
+/*! Random Walk on Dirichlet
+Using simple assignment 
+ \f[ \beta = rvc / k + \beta_c \f]
+ hence, mean value = rvc, variance = (k+1)*mean*mean;
+ 
+ The greater k is, the greater is the variance of the random walk;
+ 
+ \f$ \beta_c \f$ is used as regularizing element to avoid corner cases, i.e. when one element of rvc is zero.
+ By default is it set to 0.1;
+*/
+
+class mDirich: public mpdf_internal<eDirich> {
+        protected:
+                //! constant \f$ k \f$ of the random walk
+                double k;
+                //! cache of beta_i
+                vec &_beta;
+                //! stabilizing coefficient \f$ \beta_c \f$
+                vec betac;
+        public:
+                mDirich(): mpdf_internal<eDirich>(), _beta(iepdf._beta()){};
+                void condition (const vec &val) {_beta =  val/k+betac; };
+                /*! Create Dirichlet random walk 
+                \f[ f(rv|rvc) = Di(rvc*k) \f]
+                from structure
+                \code
+                class = 'mDirich';
+                k = 1;                      // multiplicative constant k
+                --- optional ---
+                rv = RV({'name'},size)      // description of RV
+                beta0 = [];                 // initial value of beta
+                betac = [];                 // initial value of beta
+                \endcode
+                */
+                void from_setting (const Setting &set) {
+                        mpdf::from_setting (set); // reads rv and rvc
+                        if (_rv()._dsize()>0){
+                                rvc = _rv().copy_t(-1);
+                        }
+                        vec beta0; 
+                        if (!UI::get (beta0, set, "beta0", UI::optional)){
+                                beta0 = ones(_rv()._dsize());
+                        }
+                        if (!UI::get (betac, set, "betac", UI::optional)){
+                                betac = 0.1*ones(_rv()._dsize());
+                        }
+                        _beta = beta0;
+                        
+                        UI::get (k, set, "k", UI::compulsory);
+                        validate();
+                }
+                void validate() { 
+                        iepdf.validate();
+                        bdm_assert(_beta.length()==betac.length(),"beta0 and betac are not compatible");
+                        if (_rv()._dsize()>0){
+                                bdm_assert( (_rv()._dsize()==dimension()) , "Size of rv does not match with beta");
+                        }
+                        dimc = _beta.length();
+                };
+};
+UIREGISTER(mDirich);
+
 //! \brief Estimator for Multinomial density
 class multiBM : public BMEF
@@ -1196 +1258 @@
                 Array<vec>& _samples() {return samples;};
                 //! access function
+                const vec& _sample(int i) const {return samples(i);};
+                //! access function
                 const Array<vec>& _samples() const {return samples;};
                 //! Function performs resampling, i.e. removal of low-weight samples and duplication of high-weight samples such that the new samples represent the same density.
-                ivec resample (RESAMPLING_METHOD method = SYSTEMATIC);
-
+                //! The vector with indeces of new samples is returned in variable \c index. 
+                void resample ( ivec &index, RESAMPLING_METHOD method = SYSTEMATIC);
+
+                //! Resampling without returning index of new particles. 
+                void resample (RESAMPLING_METHOD method = SYSTEMATIC){ivec ind; resample(ind,method);};
+                
                 //! inherited operation : NOT implemented
                 vec sample() const {