Changeset 1033 for library

Show
Ignore:
Timestamp:
06/02/10 23:51:04 (14 years ago)
Author:
smidl
Message:

New ebeta distribution + tests

Location:
library
Files:
1 added
4 modified

Legend:

Unmodified
Added
Removed
  • library/bdm/estim/particles.h

    r1000 r1033  
    439439 
    440440                if ( particles(0)->_rv()._dsize() > 0 ) { 
    441                         bdm_assert (  particles(0)->_rv()._dsize() == est.dimension(), "Mismatch of RV " +particles(0)->_rv().to_string() +  
    442                         " of size (" +num2str(particles(0)->_rv()._dsize())+"and dimension of posterior ("+num2str(est.dimension()) + ")" ); 
     441                        bdm_assert (  particles(0)->_rv()._dsize() == est.dimension(), "MPF:: Mismatch of RV " +particles(0)->_rv().to_string() +  
     442                        " of size (" +num2str(particles(0)->_rv()._dsize())+") and dimension of posterior ("+num2str(est.dimension()) + ")" ); 
    443443                } 
    444444        } 
  • library/bdm/stat/exp_family.cpp

    r1015 r1033  
    775775                dimc = _beta.length(); 
    776776        } 
    777 } 
     777 
     778 
     779void mBeta::from_setting ( const Setting &set ) { 
     780        pdf::from_setting ( set ); // reads rv and rvc 
     781        if ( _rv()._dsize() > 0 ) { 
     782                rvc = _rv().copy_t ( -1 ); 
     783        } 
     784        if ( !UI::get ( iepdf.beta, set, "beta", UI::optional ) ) { 
     785                iepdf.beta = ones ( _rv()._dsize() ); 
     786        } 
     787        if ( !UI::get ( iepdf.alpha, set, "alpha", UI::optional ) ) { 
     788                iepdf.alpha = ones ( _rv()._dsize() ); 
     789        } 
     790        if ( !UI::get ( betac, set, "betac", UI::optional ) ) { 
     791                betac = 0.1 * ones ( _rv()._dsize() ); 
     792        } 
     793         
     794        UI::get ( k, set, "k", UI::compulsory ); 
     795} 
     796 
     797void mBeta::to_setting  (Setting  &set) const { 
     798        pdf::to_setting(set); 
     799        UI::save( iepdf.beta, set, "beta"); 
     800        UI::save( betac, set, "betac"); 
     801        UI::save ( k, set, "k" ); 
     802} 
     803 
     804} 
  • library/bdm/stat/exp_family.h

    r1020 r1033  
    552552UIREGISTER ( eDirich ); 
    553553 
     554/*! Product of Beta distributions 
     555 
     556Continuous Dirichlet density of \f$n\f$-dimensional variable \f$x\f$ 
     557\f[ 
     558f(x|\alpha,\beta)  = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\, x^{\alpha-1}(1-x)^{\beta-1}  
     559\f] 
     560is a simplification of Dirichlet to univariate case. 
     561*/ 
     562class eBeta: public eEF { 
     563        public: 
     564                //!sufficient statistics 
     565                vec alpha; 
     566                //!sufficient statistics 
     567                vec beta; 
     568        public: 
     569                //!\name Constructors 
     570                //!@{ 
     571                         
     572                        eBeta () : eEF () {}; 
     573                        eBeta ( const eBeta &B0 ) : eEF (), alpha(B0.alpha),beta(B0.beta) { 
     574                                validate(); 
     575                        }; 
     576                        //!@} 
     577                         
     578                        //! using sampling procedure from wikipedia 
     579                        vec sample() const { 
     580                                vec y ( beta.length() ); // all vectors 
     581                                for ( int i = 0; i < beta.length(); i++ ) { 
     582                                        GamRNG.setup ( alpha ( i ), 1 ); 
     583                                        #pragma omp critical 
     584                                        double Ga = GamRNG(); 
     585                                         
     586                                        GamRNG.setup ( beta ( i ), 1 ); 
     587                                        #pragma omp critical 
     588                                        double Gb = GamRNG(); 
     589                                         
     590                                        y ( i ) = Ga/(Ga+Gb); 
     591                                } 
     592                                return y; 
     593                        } 
     594                         
     595                        vec mean() const { 
     596                                return elem_div(alpha, alpha + beta); // dot-division 
     597                        }; 
     598                        vec variance() const { 
     599                                vec apb=alpha+beta; 
     600                                return elem_div (elem_mult ( alpha, beta) , 
     601                                                                 elem_mult ( elem_mult(apb,apb), apb+1 ) ); 
     602                        } 
     603                        //! In this instance, val is ... 
     604                        double evallog_nn ( const vec &val ) const { 
     605                                double tmp; 
     606                                tmp = ( alpha - 1 ) * log ( val ) + (beta-1)*log(1-val); 
     607                                return tmp; 
     608                        } 
     609                         
     610                        double lognc () const { 
     611                                double lgb = 0.0; 
     612                                for ( int i = 0; i < beta.length(); i++ ) { 
     613                                        lgb += -lgamma ( alpha(i)+beta(i) ) + lgamma(alpha(i)) + lgamma(beta(i)); 
     614                                } 
     615                                return lgb; 
     616                        } 
     617                         
     618                        /*! configuration structure 
     619                        \code 
     620                        class = 'eBeta'; 
     621                        beta  = [];           //parametr beta 
     622                        alpha = [];           //parametr alpha 
     623                        \endcode 
     624                        */ 
     625                        void from_setting ( const Setting &set ){ 
     626                                UI::get(alpha, set, "alpha", UI::compulsory); 
     627                                UI::get(beta, set, "beta", UI::compulsory); 
     628                        } 
     629                        void validate(){ 
     630                                bdm_assert(alpha.length()==beta.length(), "eBeta:: alpha and beta length do not match"); 
     631                                dim = alpha.length(); 
     632                        }; 
     633                        void to_setting ( Setting &set ) const{ 
     634                                UI::save(alpha, set, "alpha"); 
     635                                UI::save(beta, set, "beta"); 
     636                        } 
     637}; 
     638UIREGISTER ( eBeta ); 
     639 
    554640/*! Random Walk on Dirichlet 
    555641Using simple assignment 
     
    594680UIREGISTER ( mDirich ); 
    595681 
     682/*! \brief Random Walk with vector Beta distribution 
     683Using simple assignment 
     684\f{align}  
     685\alpha &= rvc / k + \beta_c \\ 
     686\beta &= (1-rvc) / k + \beta_c \\ 
     687\f{align} 
     688for each element of alpha and beta, mean value = rvc, variance = (k+1)*mean*mean; 
     689 
     690The greater \f$ k \f$ is, the greater is the variance of the random walk; 
     691 
     692\f$ \beta_c \f$ is used as regularizing element to avoid corner cases, i.e. when one element of rvc is zero. 
     693By default is it set to 0.1; 
     694*/ 
     695 
     696class mBeta: public pdf_internal<eBeta>{ 
     697        //! vector of constants \f$ k \f$ of the random walk  
     698        vec k; 
     699        //! stabilizing coefficient \f$ \beta_c \f$ 
     700        vec betac; 
     701 
     702        public: 
     703                void condition ( const vec &val ) { 
     704                        this->iepdf.alpha =  elem_div(val , k) + betac; 
     705                        this->iepdf.beta =  elem_div (1-val , k) + betac; 
     706                }; 
     707                /*! Create Beta random walk 
     708                \f[ f(rv|rvc) = \prod Beta(rvc,k) \f] 
     709                from structure 
     710                \code 
     711                class = 'mBeta'; 
     712                k = 1;                      // multiplicative constant k 
     713                --- optional --- 
     714                rv = RV({'name'},size)      // description of RV 
     715                beta  = [];                 // initial value of beta 
     716                betac = [];                 // initial value of beta stabilizing constant 
     717                \endcode 
     718                */ 
     719                void from_setting ( const Setting &set ); 
     720                void    to_setting  (Setting  &set) const; 
     721                void validate(){ 
     722                         
     723                        pdf_internal<eBeta>::validate(); 
     724                        bdm_assert(betac.length()==dimension(),"Incomaptible betac"); 
     725                        bdm_assert(k.length()==dimension(),"Incomaptible k"); 
     726                        dimc = iepdf.dimension(); 
     727                } 
     728                //!  
     729}; 
     730UIREGISTER(mBeta); 
    596731 
    597732//! \brief Estimator for Multinomial density 
  • library/tests/testsuite/epdf_test.cpp

    r886 r1033  
    9090        epdf_harness::test_config ( "edirich.cfg" ); 
    9191} 
     92TEST ( ebeta_test ) { 
     93        epdf_harness::test_config ( "ebeta.cfg" ); 
     94}