- Timestamp:
- 11/10/08 15:40:30 (16 years ago)
- Location:
- bdm
- Files:
-
- 3 modified
Legend:
- Unmodified
- Added
- Removed
-
bdm/estim/merger.cpp
r204 r205 16 16 case 3://Ratio of Bessel 17 17 coef = sqrt ( ( 3*beta-2 ) /3*beta ); 18 return log ( besselk ( 0,sq2bl*coef )) - log( besselk ( 0,sq2bl ) ) + mu;18 return log ( besselk ( 0,sq2bl*coef ) ) - log ( besselk ( 0,sq2bl ) ) + mu; 19 19 break; 20 20 case 4: … … 31 31 it_assert_debug ( rv.equal ( g0->_rv() ),"Incompatible g0" ); 32 32 //Empirical density - samples 33 eEmp eSmp ( rv,Ns );34 33 eSmp.set_parameters ( ones ( Ns ), g0 ); 35 34 Array<vec> &Smp = eSmp._samples(); //aux … … 62 61 char str[100]; 63 62 64 epdf* Mpred ;63 epdf* Mpred=Mix.predictor(rv); 65 64 vec Mix_pdf ( Ns ); 66 65 while ( !converged ) { … … 69 68 Mix.flatten ( &Mix_init ); 70 69 Mix.bayesB ( Smp_ex, w*Ns ); 70 delete Mpred; 71 71 Mpred = Mix.predictor ( rv ); // Allocation => must be deleted at the end!! 72 73 // sprintf ( str,"Mpred_mean%d",niter ); 74 // dbg << Name ( str ) << Mpred->mean(); 75 76 if ( 1 ) { 72 73 if ( 1./sum_sqr ( w ) <0.5*Ns ) { 77 74 // Generate new samples 78 75 eSmp.set_samples ( Mpred ); 79 76 for ( int i=0;i<Ns;i++ ) { 80 77 //////////// !!!!!!!!!!!!! 81 if ( Smp ( i ) ( 1 ) <0 ) {Smp ( i ) ( 1 ) = GamRNG(); }78 if ( Smp ( i ) ( 1 ) <0 ) {Smp ( i ) ( 1 ) = 0.01; } 82 79 set_col_part ( Smp_ex,i,Smp ( i ) ); 83 80 } 84 {cout<<" Eff. sample size=" << 1./sum_sqr ( w ) << endl;}81 {cout<<"Resampling =" << 1./sum_sqr ( w ) << endl;} 85 82 cout << sum ( Smp_ex,2 ) /Ns <<endl; 83 cout << Smp_ex*Smp_ex.T()/Ns << endl; 86 84 } 87 else 88 {cout<<"Eff. sample size=" << 1./sum_sqr ( w ) << endl;} 85 // sprintf ( str,"Mpred_mean%d",niter ); 86 // dbg << Name ( str ) << Mpred->mean(); 87 89 88 90 89 // sprintf ( str,"Mpdf%d",niter ); … … 135 134 zdls ( i )->get_cond ( Smp ( k ) ) ) ); 136 135 if ( !std::isfinite ( lw_src ( k ) ) ) { 137 cout << endl;136 lw_src ( k ) = -1e16; cout << "!"; 138 137 } 139 138 } … … 166 165 //Importance weighting 167 166 lw -= lw_mix; // hoping that it is not numerically sensitive... 168 w = exp (lw-max(lw));167 w = exp ( lw-max ( lw ) ); 169 168 //renormalize 170 169 double sumw =sum ( w ); … … 173 172 } 174 173 else { 175 174 176 175 } 177 { 178 double eff = 1./sum_sqr ( w ); 179 if ( eff<2 ) { 180 int mi= max_index ( w ); 181 cout << w(mi) <<endl; 182 cout << lW.get_col ( mi ) <<endl; 183 mat mm(2,1);mm.set_col(0,lW.get_col(mi)); 184 cout << lognorm_merge(mm) <<endl; 185 cout << lw_mix ( mi ) <<endl; 186 cout << lw ( mi ) <<endl; 187 cout << Smp ( mi ) <<endl; 188 cout << Mix._Coms(0)->_e()->mean() <<endl; 189 } 190 cout << "Eff: " << eff <<endl; 191 } 176 192 177 // sprintf ( str,"w_is_%d",niter ); 193 178 // dbg << Name ( str ) << w; -
bdm/estim/merger.h
r198 r205 43 43 //!Prior on the log-normal merging model 44 44 double beta; 45 //! Projection to empirical density 46 eEmp eSmp; 47 45 48 public: 46 49 //!Default constructor 47 50 merger ( const Array<mpdf*> &S ) : 48 51 compositepdf ( S ), epdf ( getrv ( false ) ), 49 Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ), dls ( n ), rvzs ( n ), zdls ( n ) {52 Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ), dls ( n ), rvzs ( n ), zdls ( n ), eSmp(rv,0) { 50 53 RV ztmp; 51 54 // Extend rv by rvc! … … 68 71 } 69 72 //! Set internal parameters used in approximation 70 void set_parameters ( double beta0, int Ns0, int Nc0 ) { beta=beta0;Ns=Ns0;Nc=Nc0;}73 void set_parameters ( double beta0, int Ns0, int Nc0 ) {beta=beta0;Ns=Ns0;Nc=Nc0;eSmp.set_n(Ns0,false);} 71 74 //!Initialize the proposal density. This function must be called before merge()! 72 75 void init() { ////////////// NOT FINISHED … … 90 93 return Mix.logpred ( dtf ); 91 94 } 92 vec mean() const {return Mix._epdf().mean();} 95 vec mean() const { 96 const Vec<double> &w = eSmp._w(); 97 const Array<vec> &S = eSmp._samples(); 98 vec tmp=zeros ( rv.count() ); 99 for ( int i=0; i<Ns; i++ ) { 100 tmp+=w ( i ) *S ( i ); 101 } 102 return tmp; 103 } 104 mat variance() const { 105 const vec &w = eSmp._w(); 106 const Array<vec> &S = eSmp._samples(); 107 108 vec mea = mean(); 109 110 cout << sum(w) << "," << w*w <<endl; 111 112 mat Tmp=zeros(rv.count(), rv.count()); 113 for ( int i=0; i<Ns; i++ ) { 114 Tmp+=w ( i ) *outer_product(S ( i ), S(i)); 115 } 116 return Tmp-outer_product(mea,mea); 117 } 93 118 //! for future use 94 119 virtual ~merger() { … … 101 126 //! Access function 102 127 MixEF& _Mix() {return Mix;} 128 //! Access function 129 eEmp& _Smp() {return eSmp;} 103 130 }; 104 131 -
bdm/stat/libEF.h
r204 r205 200 200 ldmat& _V() {return V;} 201 201 //! returns a pointer to the internal statistics. Use with Care! 202 double& _nu() {return nu;} 202 const ldmat& _V() const {return V;} 203 //! returns a pointer to the internal statistics. Use with Care! 204 double& _nu() {return nu;} 205 const double& _nu() const {return nu;} 203 206 void pow ( double p ) {V*=p;nu*=p;}; 204 207 }; … … 520 523 //! Set sample 521 524 void set_samples ( const epdf* pdf0 ); 525 //! Set sample 526 void set_n ( int n0, bool copy=true ){w.set_size(n0,copy);samples.set_size(n0,copy);}; 522 527 //! Potentially dangerous, use with care. 523 528 vec& _w() {return w;}; 529 //! Potentially dangerous, use with care. 530 const vec& _w() const {return w;}; 524 531 //! access function 525 532 Array<vec>& _samples() {return samples;}; 533 //! access function 534 const Array<vec>& _samples() const {return samples;}; 526 535 //! Function performs resampling, i.e. removal of low-weight samples and duplication of high-weight samples such that the new samples represent the same density. 527 536 ivec resample ( RESAMPLING_METHOD method = SYSTEMATIC );