Changeset 192 for bdm/estim

Show
Ignore:
Timestamp:
10/22/08 10:46:40 (16 years ago)
Author:
smidl
Message:

modification of datalinks and switch mprod and merger to use them

Location:
bdm/estim
Files:
2 modified

Legend:

Unmodified
Added
Removed
  • bdm/estim/merger.cpp

    r190 r192  
    3333        vec &w = eSmp._w(); //aux 
    3434 
    35         mat Smp_ex =ones ( rv.count() +1,Ns ); // Extended samples for the ARX model - the last row is ones  
     35        mat Smp_ex =ones ( rv.count() +1,Ns ); // Extended samples for the ARX model - the last row is ones 
    3636        for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );} 
    3737 
     
    6060                //Re-Initialize Mixture model 
    6161                Mix.init ( &A0, Smp_ex, Nc ); 
    62                 Mix.bayesB ( Smp_ex, ones(Ns));//w*Ns ); 
    63                 Mpred = Mix.predictor(rv); // Allocation => must be deleted at the end!! 
    64          
     62                Mix.bayesB ( Smp_ex, ones ( Ns ) );//w*Ns ); 
     63                Mpred = Mix.predictor ( rv ); // Allocation => must be deleted at the end!! 
     64 
    6565                // Generate new samples 
    6666                eSmp.set_samples ( Mpred ); 
     
    7979                        //======== Same RVs =========== 
    8080                        //Split according to dependency in rvs 
    81                         if ( rvsinrv ( i ).length() ==rv.count() ) { 
     81                        if ( mpdfs ( i )->_rv().count() ==rv.count() ) { 
    8282                                // no need for conditioning or marginalization 
    8383                                for ( int j=0;j<Ns; j++ ) { // Smp is Array<> => for cycle 
     
    8686                        } 
    8787                        else { 
    88                                 vec smpk; 
    8988                                // compute likelihood of marginal on the conditional variable 
    9089                                if ( mpdfs ( i )->_rvc().count() >0 ) { 
    9190                                        // Make marginal on rvc_i 
    9291                                        epdf* tmp_marg = Mpred->marginal ( mpdfs ( i )->_rvc() ); 
    93                                         for(int k=0;k<Ns;k++){ 
    94                                         lw_src(k) += tmp_marg->evalpdflog ( get_vec(Smp(i), irv_rvcs(i) ) );} 
     92                                        //compute vector of lw_src 
     93                                        for ( int k=0;k<Ns;k++ ) { 
     94                                                lw_src ( k ) += tmp_marg->evalpdflog ( dls ( i )->get_val ( Smp ( i ) ) ); 
     95                                        } 
    9596                                        delete tmp_marg; 
    9697                                } 
    9798                                // Compute likelihood of the missing variable 
    98                                 if ( rv.count() > mpdfs ( i )->_rv().count() + mpdfs ( i )->_rvc().count() ) { 
    99                                         // There are variales unknown to mpdfs 
    100                                         RV z = ( rv.subt ( mpdfs ( i )->_rv() ) ).subt ( mpdfs ( i )->_rvc() ); 
    101                                         mpdf* tmp_cond = Mpred->condition ( z );  
    102                                         // Indeces of rest rv in Smp 
    103                                         ivec zinrv=z.dataind ( rv ) ; 
    104                                         // Indeces of rest rvc in Smp 
    105                                         ivec zinrvc=tmp_cond->_rvc().dataind ( rv ); 
     99                                if ( rv.count() > ( mpdfs ( i )->_rv().count() + mpdfs ( i )->_rvc().count() ) ) { 
     100                                        // There are variales unknown to mpdfs(i) : rvzs 
     101                                        mpdf* tmp_cond = Mpred->condition ( rvzs ( i ) ); 
    106102                                        // Compute likelihood 
    107103                                        for ( int k= 0; k<Ns; k++ ) { 
    108                                                 smpk=Smp( k ); 
    109                                                 lw_src ( k ) += log( tmp_cond->evalcond ( get_vec ( smpk,zinrv ), get_vec ( smpk,zinrvc ) )); 
     104                                                lw_src ( k ) += log ( 
     105                                                                    tmp_cond->evalcond ( 
     106                                                                        zdls ( i )->get_val ( Smp ( k ) ), 
     107                                                                        zdls ( i )->get_cond ( Smp ( k ) ) ) ); 
    110108                                        } 
    111109                                        delete tmp_cond; 
     
    113111                                // Compute likelihood of the partial source 
    114112                                for ( int k= 0; k<Ns; k++ ) { 
    115                                         smpk=Smp( k ); 
    116                                         mpdfs ( i )->condition ( get_vec ( smpk,irv_rvcs(i) ) ); 
    117                                         lw_src ( k ) += mpdfs ( i )->_epdf().evalpdflog ( get_vec ( smpk, rvsinrv ( i ) ) ); 
     113                                        mpdfs ( i )->condition ( dls ( i )->get_cond ( Smp(k) ) ); 
     114                                        lw_src ( k ) += mpdfs ( i )->_epdf().evalpdflog ( dls ( i )->get_val ( Smp(k) ) ); 
    118115                                } 
    119116                        } 
     
    128125 
    129126                w = lognorm_merge ( lW ); //merge 
    130                  
     127 
    131128                sprintf ( str,"w%d",niter ); 
    132129                dbg << Name ( str ) << w; 
     
    150147                // ==== stopping rule === 
    151148                niter++; 
    152                 converged = ( niter>20); 
     149                converged = ( niter>20 ); 
    153150        } 
    154151 
  • bdm/estim/merger.h

    r190 r192  
    3030        //!Internal mixture of EF models 
    3131        MixEF Mix; 
     32        //! Data link for each mpdf in mpdfs 
     33        Array<datalink_m2e*> dls; 
     34        //! Array of rvs that are not modelled by mpdfs at all (aux) 
     35        Array<RV> rvzs; 
     36        //! Data Links of rv0 mpdfs - these will be conditioned the (rv,rvc) of mpdfs 
     37        Array<datalink_m2e*> zdls;  
     38         
    3239        //!Number of samples used in approximation 
    3340        int Ns; 
     
    4047        merger ( const Array<mpdf*> &S ) : 
    4148                        compositepdf ( S ), epdf ( getrv ( false ) ), 
    42                         Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ) 
    43         {setindices(rv); beta=2.0; Ns=100; Nc=10; Mix.set_method(EM);} 
    44         //! Set internal parameters used in approximation 
     49                        Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ), dls(n), rvzs(n), zdls(n) { 
     50                                RV ztmp; 
     51                for ( int i=0;i<n;i++ ) { 
     52                        //Establich connection between mpdfs and merger 
     53                        dls ( i ) = new datalink_m2e ( mpdfs ( i )->_rv(), mpdfs(i)->_rvc(), rv ); 
     54                        // find out what is missing in each mpdf 
     55                        ztmp= mpdfs(i)->_rv(); 
     56                        ztmp.add(mpdfs(i)->_rvc()); 
     57                        rvzs(i)=rv.subt(ztmp); 
     58                        zdls(i) = new datalink_m2e (rvzs(i), ztmp, rv ) ; 
     59                }; 
     60                //Set Default values of parameters 
     61                beta=2.0; 
     62                Ns=100; 
     63                Nc=10; 
     64                Mix.set_method ( EM ); 
     65        } 
     66//! Set internal parameters used in approximation 
    4567        void set_parameters ( double beta0, int Ns0, int Nc0 ) {        beta=beta0;Ns=Ns0;Nc=Nc0;} 
    46         //!Initialize the proposal density. This function must be called before merge()! 
    47         void init() { 
     68//!Initialize the proposal density. This function must be called before merge()! 
     69        void init() { ////////////// NOT FINISHED 
    4870                Array<vec> Smps ( n ); 
    4971                //Gibbs sampling 
    5072                for ( int i=0;i<n;i++ ) {Smps ( i ) =zeros ( 0 );} 
    5173        } 
    52         //!Create a mixture density using known proposal 
     74//!Create a mixture density using known proposal 
    5375        void merge ( const epdf* g0 ); 
    54         //!Create a mixture density, make sure to call init() before the first call 
     76//!Create a mixture density, make sure to call init() before the first call 
    5577        void merge () {merge ( & ( Mix._epdf() ) );}; 
    5678 
    57         //! Merge log-likelihood values 
     79//! Merge log-likelihood values 
    5880        vec lognorm_merge ( mat &lW ); 
    59         //! sample from merged density 
    60         //! weight w is a 
    61         vec sample ( )const { return Mix._epdf().sample();} 
    62         double evalpdflog ( const vec &dt ) const{  
    63                 vec dtf=ones(dt.length()+1);     
    64                 dtf.set_subvector(0,dt);  
    65                 return Mix.logpred ( dtf );} 
    66         vec mean()const {return Mix._epdf().mean();} 
    67         //! for future use 
    68         virtual ~merger() {}; 
    69          
    70         //! Access function 
     81//! sample from merged density 
     82//! weight w is a 
     83        vec sample ( ) const { return Mix._epdf().sample();} 
     84        double evalpdflog ( const vec &dt ) const { 
     85                vec dtf=ones ( dt.length() +1 ); 
     86                dtf.set_subvector ( 0,dt ); 
     87                return Mix.logpred ( dtf ); 
     88        } 
     89        vec mean() const {return Mix._epdf().mean();} 
     90//! for future use 
     91        virtual ~merger() { 
     92                for (int i=0; i<n; i++){ 
     93                        delete dls(i); 
     94                        delete zdls(i);                  
     95                } 
     96        }; 
     97 
     98//! Access function 
    7199        MixEF& _Mix() {return Mix;} 
    72100};