Changeset 1033

Timestamp:

06/02/10 23:51:04 (14 years ago)

Author:

smidl

Message:

New ebeta distribution + tests

Location:

library

Files:

• bdm/estim/particles.h (modified) (1 diff)
• bdm/stat/exp_family.cpp (modified) (1 diff)
• bdm/stat/exp_family.h (modified) (2 diffs)
• tests/testsuite/ebeta.cfg (added)
• tests/testsuite/epdf_test.cpp (modified) (1 diff)

library/bdm/estim/particles.h

@@ -439 +439 @@
 
                 if ( particles(0)->_rv()._dsize() > 0 ) {
-                        bdm_assert (  particles(0)->_rv()._dsize() == est.dimension(), "Mismatch of RV " +particles(0)->_rv().to_string() + 
-                        " of size (" +num2str(particles(0)->_rv()._dsize())+"and dimension of posterior ("+num2str(est.dimension()) + ")" );
+                        bdm_assert (  particles(0)->_rv()._dsize() == est.dimension(), "MPF:: Mismatch of RV " +particles(0)->_rv().to_string() + 
+                        " of size (" +num2str(particles(0)->_rv()._dsize())+") and dimension of posterior ("+num2str(est.dimension()) + ")" );
                 }
         }

library/bdm/stat/exp_family.cpp

@@ -775 +775 @@
                 dimc = _beta.length();
         }
-}
+
+
+void mBeta::from_setting ( const Setting &set ) {
+        pdf::from_setting ( set ); // reads rv and rvc
+        if ( _rv()._dsize() > 0 ) {
+                rvc = _rv().copy_t ( -1 );
+        }
+        if ( !UI::get ( iepdf.beta, set, "beta", UI::optional ) ) {
+                iepdf.beta = ones ( _rv()._dsize() );
+        }
+        if ( !UI::get ( iepdf.alpha, set, "alpha", UI::optional ) ) {
+                iepdf.alpha = ones ( _rv()._dsize() );
+        }
+        if ( !UI::get ( betac, set, "betac", UI::optional ) ) {
+                betac = 0.1 * ones ( _rv()._dsize() );
+        }
+        
+        UI::get ( k, set, "k", UI::compulsory );
+}
+
+void mBeta::to_setting  (Setting  &set) const {
+        pdf::to_setting(set);
+        UI::save( iepdf.beta, set, "beta");
+        UI::save( betac, set, "betac");
+        UI::save ( k, set, "k" );
+}
+
+}

library/bdm/stat/exp_family.h

@@ -552 +552 @@
 UIREGISTER ( eDirich );
 
+/*! Product of Beta distributions
+
+Continuous Dirichlet density of \f$n\f$-dimensional variable \f$x\f$
+\f[
+f(x|\alpha,\beta)  = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\, x^{\alpha-1}(1-x)^{\beta-1} 
+\f]
+is a simplification of Dirichlet to univariate case.
+*/
+class eBeta: public eEF {
+        public:
+                //!sufficient statistics
+                vec alpha;
+                //!sufficient statistics
+                vec beta;
+        public:
+                //!\name Constructors
+                //!@{
+                        
+                        eBeta () : eEF () {};
+                        eBeta ( const eBeta &B0 ) : eEF (), alpha(B0.alpha),beta(B0.beta) {
+                                validate();
+                        };
+                        //!@}
+                        
+                        //! using sampling procedure from wikipedia
+                        vec sample() const {
+                                vec y ( beta.length() ); // all vectors
+                                for ( int i = 0; i < beta.length(); i++ ) {
+                                        GamRNG.setup ( alpha ( i ), 1 );
+                                        #pragma omp critical
+                                        double Ga = GamRNG();
+                                        
+                                        GamRNG.setup ( beta ( i ), 1 );
+                                        #pragma omp critical
+                                        double Gb = GamRNG();
+                                        
+                                        y ( i ) = Ga/(Ga+Gb);
+                                }
+                                return y;
+                        }
+                        
+                        vec mean() const {
+                                return elem_div(alpha, alpha + beta); // dot-division
+                        };
+                        vec variance() const {
+                                vec apb=alpha+beta;
+                                return elem_div (elem_mult ( alpha, beta) ,
+                                                                 elem_mult ( elem_mult(apb,apb), apb+1 ) );
+                        }
+                        //! In this instance, val is ...
+                        double evallog_nn ( const vec &val ) const {
+                                double tmp;
+                                tmp = ( alpha - 1 ) * log ( val ) + (beta-1)*log(1-val);
+                                return tmp;
+                        }
+                        
+                        double lognc () const {
+                                double lgb = 0.0;
+                                for ( int i = 0; i < beta.length(); i++ ) {
+                                        lgb += -lgamma ( alpha(i)+beta(i) ) + lgamma(alpha(i)) + lgamma(beta(i));
+                                }
+                                return lgb;
+                        }
+                        
+                        /*! configuration structure
+                        \code
+                        class = 'eBeta';
+                        beta  = [];           //parametr beta
+                        alpha = [];           //parametr alpha
+                        \endcode
+                        */
+                        void from_setting ( const Setting &set ){
+                                UI::get(alpha, set, "alpha", UI::compulsory);
+                                UI::get(beta, set, "beta", UI::compulsory);
+                        }
+                        void validate(){
+                                bdm_assert(alpha.length()==beta.length(), "eBeta:: alpha and beta length do not match");
+                                dim = alpha.length();
+                        };
+                        void to_setting ( Setting &set ) const{
+                                UI::save(alpha, set, "alpha");
+                                UI::save(beta, set, "beta");
+                        }
+};
+UIREGISTER ( eBeta );
+
 /*! Random Walk on Dirichlet
 Using simple assignment
@@ -594 +680 @@
 UIREGISTER ( mDirich );
 
+/*! \brief Random Walk with vector Beta distribution
+Using simple assignment
+\f{align} 
+\alpha &= rvc / k + \beta_c \\
+\beta &= (1-rvc) / k + \beta_c \\
+\f{align}
+for each element of alpha and beta, mean value = rvc, variance = (k+1)*mean*mean;
+
+The greater \f$ k \f$ 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 mBeta: public pdf_internal<eBeta>{
+        //! vector of constants \f$ k \f$ of the random walk 
+        vec k;
+        //! stabilizing coefficient \f$ \beta_c \f$
+        vec betac;
+
+        public:
+                void condition ( const vec &val ) {
+                        this->iepdf.alpha =  elem_div(val , k) + betac;
+                        this->iepdf.beta =  elem_div (1-val , k) + betac;
+                };
+                /*! Create Beta random walk
+                \f[ f(rv|rvc) = \prod Beta(rvc,k) \f]
+                from structure
+                \code
+                class = 'mBeta';
+                k = 1;                      // multiplicative constant k
+                --- optional ---
+                rv = RV({'name'},size)      // description of RV
+                beta  = [];                 // initial value of beta
+                betac = [];                 // initial value of beta stabilizing constant
+                \endcode
+                */
+                void from_setting ( const Setting &set );
+                void    to_setting  (Setting  &set) const;
+                void validate(){
+                        
+                        pdf_internal<eBeta>::validate();
+                        bdm_assert(betac.length()==dimension(),"Incomaptible betac");
+                        bdm_assert(k.length()==dimension(),"Incomaptible k");
+                        dimc = iepdf.dimension();
+                }
+                //! 
+};
+UIREGISTER(mBeta);
 
 //! \brief Estimator for Multinomial density

library/tests/testsuite/epdf_test.cpp

@@ -90 +90 @@
         epdf_harness::test_config ( "edirich.cfg" );
 }
+TEST ( ebeta_test ) {
+        epdf_harness::test_config ( "ebeta.cfg" );
+}