Changeset 379 for bdm/estim/merger.cpp
- Timestamp:
- 06/17/09 23:53:11 (15 years ago)
- Files:
-
- 1 modified
-
bdm/estim/merger.cpp (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
bdm/estim/merger.cpp
r311 r379 5 5 namespace bdm 6 6 { 7 vec merger ::lognorm_merge( mat &lW )7 vec merger_base::merge_points ( mat &lW ) 8 8 { 9 9 int nu=lW.rows(); 10 vec mu = sum ( lW ) /nu; //mean of logs 11 // vec lam = sum ( pow ( lW,2 ) )-nu*pow ( mu,2 ); ======= numerically unsafe! 12 vec lam = sum ( pow ( lW-outer_product ( ones ( lW.rows() ),mu ),2 ) ); 13 double coef=0.0; 14 vec sq2bl=sqrt ( 2*beta*lam ); //this term is everywhere 15 switch ( nu ) 10 switch ( METHOD ) 16 11 { 17 case 2: 18 coef= ( 1-0.5*sqrt ( ( 4.0*beta-3.0 ) /beta ) ); 19 return coef*sq2bl + mu ; 12 case ARITHMETIC: 13 return log ( sum ( exp ( lW ) ) ); //ugly! 20 14 break; 21 case 3://Ratio of Bessel 22 coef = sqrt ( ( 3*beta-2 ) /3*beta ); 23 return log ( besselk ( 0,sq2bl*coef ) ) - log ( besselk ( 0,sq2bl ) ) + mu; 15 case GEOMETRIC: 16 return sum ( lW ) /nu; 24 17 break; 25 case 4: 26 break; 27 default: // Approximate conditional density 18 case LOGNORMAL: 19 vec mu = sum ( lW ) /nu; //mean of logs 20 // vec lam = sum ( pow ( lW,2 ) )-nu*pow ( mu,2 ); ======= numerically unsafe! 21 vec lam = sum ( pow ( lW-outer_product ( ones ( lW.rows() ),mu ),2 ) ); 22 double coef=0.0; 23 vec sq2bl=sqrt ( 2*beta*lam ); //this term is everywhere 24 switch ( nu ) 25 { 26 case 2: 27 coef= ( 1-0.5*sqrt ( ( 4.0*beta-3.0 ) /beta ) ); 28 return coef*sq2bl + mu ; 29 break; 30 case 3://Ratio of Bessel 31 coef = sqrt ( ( 3*beta-2 ) /3*beta ); 32 return log ( besselk ( 0,sq2bl*coef ) ) - log ( besselk ( 0,sq2bl ) ) + mu; 33 break; 34 case 4: 35 break; 36 default: // Approximate conditional density 37 break; 38 } 28 39 break; 29 40 } … … 31 42 } 32 43 33 void merger ::merge ( const epdf* g0)44 void merger_mix::merge ( ) 34 45 { 35 36 it_assert_debug ( rv.equal ( g0->_rv() ),"Incompatible g0" );37 //Empirical density - samples38 if ( !fix_smp )39 {40 eSmp.set_statistics ( ones ( Ns ), g0 );41 }42 43 46 Array<vec> &Smp = eSmp._samples(); //aux 44 47 vec &w = eSmp._w(); //aux 45 48 46 mat Smp_ex =ones ( dim +1,N s ); // Extended samples for the ARX model - the last row is ones47 for ( int i=0;i<N s;i++ ) { set_col_part ( Smp_ex,i,Smp ( i ) );}48 49 if ( DBG ) *dbg << Name ( "Smp_0" ) << Smp_ex;49 mat Smp_ex =ones ( dim +1,Npoints ); // Extended samples for the ARX model - the last row is ones 50 for ( int i=0;i<Npoints;i++ ) { set_col_part ( Smp_ex,i,Smp ( i ) );} 51 52 if ( DBG ) *dbg_file << Name ( "Smp_0" ) << Smp_ex; 50 53 51 54 // Stuff for merging 52 vec lw_src ( N s );53 vec lw_mix ( N s );54 vec lw ( N s );55 mat lW=zeros ( n,Ns );55 vec lw_src ( Npoints ); // weights of the ith source 56 vec lw_mix ( Npoints ); // weights of the approximating mixture 57 vec lw ( Npoints ); // tmp 58 mat lW=zeros ( Nsources,Npoints ); // array of weights of all sources 56 59 vec vec0 ( 0 ); 57 60 58 61 //initialize importance weights 59 if ( !fix_smp ) 60 for ( int i=0;i<Ns;i++ ) 61 { 62 lw_mix ( i ) =g0->evallog ( Smp ( i ) ); 63 } 62 lw_mix = 1.0; // assuming uniform grid density -- otherwise 64 63 65 64 // Initial component in the mixture model … … 67 66 ARX A0; A0.set_statistics ( dim, V0 ); //initial guess of Mix: 68 67 69 Mix.init ( &A0, Smp_ex, Nc );68 Mix.init ( &A0, Smp_ex, Ncoms ); 70 69 //Preserve initial mixture for repetitive estimation via flattening 71 70 MixEF Mix_init ( Mix ); … … 74 73 bool converged=false; 75 74 int niter = 0; 76 char str[100];75 char dbg_str[100]; 77 76 78 77 emix* Mpred=Mix.epredictor ( ); 79 vec Mix_pdf ( N s );78 vec Mix_pdf ( Npoints ); 80 79 while ( !converged ) 81 80 { … … 83 82 //Re-Initialize Mixture model 84 83 Mix.flatten ( &Mix_init ); 85 Mix.bayesB ( Smp_ex, w*N s );84 Mix.bayesB ( Smp_ex, w*Npoints ); 86 85 delete Mpred; 87 86 Mpred = Mix.epredictor ( ); // Allocation => must be deleted at the end!! … … 89 88 90 89 // This will be active only later in iterations!!! 91 if ( ( !fix_smp ) & ( 1./sum_sqr ( w ) <effss_coef*Ns ) )90 if ( 1./sum_sqr ( w ) <effss_coef*Npoints ) 92 91 { 93 92 // Generate new samples 94 93 eSmp.set_samples ( Mpred ); 95 for ( int i=0;i<N s;i++ )94 for ( int i=0;i<Npoints;i++ ) 96 95 { 97 96 //////////// !!!!!!!!!!!!! 98 if ( Smp ( i ) ( 2 ) <0 ) {Smp ( i ) ( 2 ) = 0.01; }97 //if ( Smp ( i ) ( 2 ) <0 ) {Smp ( i ) ( 2 ) = 0.01; } 99 98 set_col_part ( Smp_ex,i,Smp ( i ) ); 100 99 //Importance of the mixture … … 102 101 lw_mix ( i ) = Mpred->evallog ( Smp ( i ) ); 103 102 } 104 if ( 0)103 if ( DBG ) 105 104 { 106 105 cout<<"Resampling =" << 1./sum_sqr ( w ) << endl; 107 106 cout << Mix._e()->mean() <<endl; 108 cout << sum ( Smp_ex,2 ) /N s <<endl;109 cout << Smp_ex*Smp_ex.T() /N s << endl;107 cout << sum ( Smp_ex,2 ) /Npoints <<endl; 108 cout << Smp_ex*Smp_ex.T() /Npoints << endl; 110 109 } 111 110 } 112 111 if ( DBG ) 113 112 { 114 sprintf ( str,"Mpred_mean%d",niter );115 *dbg << Name (str ) << Mpred->mean();116 sprintf ( str,"Mpred_var%d",niter );117 *dbg << Name (str ) << Mpred->variance();118 119 120 sprintf ( str,"Mpdf%d",niter );121 for ( int i=0;i<N s;i++ ) {Mix_pdf ( i ) = Mix.logpred ( Smp_ex.get_col ( i ) );}122 *dbg << Name (str ) << Mix_pdf;123 124 sprintf ( str,"Smp%d",niter );125 *dbg << Name (str ) << Smp_ex;113 sprintf ( dbg_str,"Mpred_mean%d",niter ); 114 *dbg_file << Name ( dbg_str ) << Mpred->mean(); 115 sprintf ( dbg_str,"Mpred_var%d",niter ); 116 *dbg_file << Name ( dbg_str ) << Mpred->variance(); 117 118 119 sprintf ( dbg_str,"Mpdf%d",niter ); 120 for ( int i=0;i<Npoints;i++ ) {Mix_pdf ( i ) = Mix.logpred ( Smp_ex.get_col ( i ) );} 121 *dbg_file << Name ( dbg_str ) << Mix_pdf; 122 123 sprintf ( dbg_str,"Smp%d",niter ); 124 *dbg_file << Name ( dbg_str ) << Smp_ex; 126 125 127 126 } 128 127 //Importace weighting 129 for ( int i=0;i< n;i++ )128 for ( int i=0;i<mpdfs.length();i++ ) 130 129 { 131 130 lw_src=0.0; … … 135 134 { 136 135 // no need for conditioning or marginalization 137 for ( int j=0;j<N s; j++ ) // Smp is Array<> => for cycle136 for ( int j=0;j<Npoints; j++ ) // Smp is Array<> => for cycle 138 137 { 139 138 lw_src ( j ) =mpdfs ( i )->_epdf().evallog ( Smp ( j ) ); … … 148 147 epdf* tmp_marg = Mpred->marginal ( mpdfs ( i )->_rvc() ); 149 148 //compute vector of lw_src 150 for ( int k=0;k<N s;k++ )149 for ( int k=0;k<Npoints;k++ ) 151 150 { 152 151 // Here val of tmp_marg = cond of mpdfs(i) ==> calling dls->get_cond … … 167 166 // Compute likelihood 168 167 vec lw_dbg=lw_src; 169 for ( int k= 0; k<N s; k++ )168 for ( int k= 0; k<Npoints; k++ ) 170 169 { 171 170 lw_src ( k ) += log ( … … 181 180 } 182 181 // Compute likelihood of the partial source 183 for ( int k= 0; k<N s; k++ )182 for ( int k= 0; k<Npoints; k++ ) 184 183 { 185 184 mpdfs ( i )->condition ( dls ( i )->get_cond ( Smp ( k ) ) ); … … 188 187 189 188 } 190 // it_assert_debug(std::isfinite(sum(lw_src)),"bad");189 // it_assert_debug(std::isfinite(sum(lw_src)),"bad"); 191 190 lW.set_row ( i, lw_src ); // do not divide by mix 192 191 } 193 lw = lognorm_merge( lW ); //merge192 lw = merger_base::merge_points ( lW ); //merge 194 193 195 194 //Importance weighting 196 if ( !fix_smp ) 197 lw -= lw_mix; // hoping that it is not numerically sensitive... 195 lw -= lw_mix; // hoping that it is not numerically sensitive... 198 196 w = exp ( lw-max ( lw ) ); 199 197 … … 212 210 if ( DBG ) 213 211 { 214 sprintf ( str,"lW%d",niter );215 *dbg << Name (str ) << lW;216 sprintf ( str,"w%d",niter );217 *dbg << Name (str ) << w;218 sprintf ( str,"lw_m%d",niter );219 *dbg << Name (str ) << lw_mix;212 sprintf ( dbg_str,"lW%d",niter ); 213 *dbg_file << Name ( dbg_str ) << lW; 214 sprintf ( dbg_str,"w%d",niter ); 215 *dbg_file << Name ( dbg_str ) << w; 216 sprintf ( dbg_str,"lw_m%d",niter ); 217 *dbg_file << Name ( dbg_str ) << lw_mix; 220 218 } 221 219 // ==== stopping rule ===
