- Timestamp:
- 10/22/08 10:46:40 (16 years ago)
- Location:
- bdm/estim
- Files:
-
- 2 modified
Legend:
- Unmodified
- Added
- Removed
-
bdm/estim/merger.cpp
r190 r192 33 33 vec &w = eSmp._w(); //aux 34 34 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 36 36 for ( int i=0;i<Ns;i++ ) { set_col_part ( Smp_ex,i,Smp ( i ) );} 37 37 … … 60 60 //Re-Initialize Mixture model 61 61 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 65 65 // Generate new samples 66 66 eSmp.set_samples ( Mpred ); … … 79 79 //======== Same RVs =========== 80 80 //Split according to dependency in rvs 81 if ( rvsinrv ( i ).length() ==rv.count() ) {81 if ( mpdfs ( i )->_rv().count() ==rv.count() ) { 82 82 // no need for conditioning or marginalization 83 83 for ( int j=0;j<Ns; j++ ) { // Smp is Array<> => for cycle … … 86 86 } 87 87 else { 88 vec smpk;89 88 // compute likelihood of marginal on the conditional variable 90 89 if ( mpdfs ( i )->_rvc().count() >0 ) { 91 90 // Make marginal on rvc_i 92 91 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 } 95 96 delete tmp_marg; 96 97 } 97 98 // 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 ) ); 106 102 // Compute likelihood 107 103 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 ) ) ) ); 110 108 } 111 109 delete tmp_cond; … … 113 111 // Compute likelihood of the partial source 114 112 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) ) ); 118 115 } 119 116 } … … 128 125 129 126 w = lognorm_merge ( lW ); //merge 130 127 131 128 sprintf ( str,"w%d",niter ); 132 129 dbg << Name ( str ) << w; … … 150 147 // ==== stopping rule === 151 148 niter++; 152 converged = ( niter>20 );149 converged = ( niter>20 ); 153 150 } 154 151 -
bdm/estim/merger.h
r190 r192 30 30 //!Internal mixture of EF models 31 31 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 32 39 //!Number of samples used in approximation 33 40 int Ns; … … 40 47 merger ( const Array<mpdf*> &S ) : 41 48 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 45 67 void set_parameters ( double beta0, int Ns0, int Nc0 ) { beta=beta0;Ns=Ns0;Nc=Nc0;} 46 47 void init() { 68 //!Initialize the proposal density. This function must be called before merge()! 69 void init() { ////////////// NOT FINISHED 48 70 Array<vec> Smps ( n ); 49 71 //Gibbs sampling 50 72 for ( int i=0;i<n;i++ ) {Smps ( i ) =zeros ( 0 );} 51 73 } 52 74 //!Create a mixture density using known proposal 53 75 void merge ( const epdf* g0 ); 54 76 //!Create a mixture density, make sure to call init() before the first call 55 77 void merge () {merge ( & ( Mix._epdf() ) );}; 56 78 57 79 //! Merge log-likelihood values 58 80 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 71 99 MixEF& _Mix() {return Mix;} 72 100 };