- Timestamp:
- 10/15/08 19:08:06 (16 years ago)
- Location:
- bdm
- Files:
-
- 9 modified
Legend:
- Unmodified
- Added
- Removed
-
bdm/estim/arx.cpp
r174 r180 54 54 set_parameters ( A0->V,A0->nu ); 55 55 } 56 57 enorm<ldmat>* ARX::predictor(const RV &rv){ 58 mat mu(rv.count(), 1); 59 mat R(rv.count(),rv.count()); 60 enorm<ldmat>* tmp; 61 tmp=new enorm<ldmat>(rv); 62 63 est.mean_mat(mu,R); 64 tmp->set_parameters(mu.get_col(0),ldmat(R)); 65 return tmp; 66 } 67 56 68 /*! \brief Return the best structure 57 69 @param Eg a copy of GiW density that is being examined -
bdm/estim/arx.h
r174 r180 66 66 const epdf& _epdf() const {return est;} 67 67 double logpred ( const vec &dt ) const; 68 void flatten ( BMEF* B ) {69 ARX* A=dynamic_cast<ARX*>(B);68 void flatten (const BMEF* B ) { 69 const ARX* A=dynamic_cast<const ARX*>(B); 70 70 // nu should be equal to B.nu 71 71 est.pow ( A->nu/nu); … … 73 73 } 74 74 75 enorm<ldmat>* predictor(const RV &rv); 75 76 //! Brute force structure estimation.\return indeces of accepted regressors. 76 77 ivec structure_est ( egiw Eg0 ); -
bdm/estim/merger.cpp
r176 r180 18 18 break; 19 19 } 20 return vec (0);20 return vec ( 0 ); 21 21 } 22 22 23 23 void merger::merge ( const epdf* g0 ) { 24 mat Smp ( rv.count(),Ns ); 25 mat lW=zeros ( n,Ns ); 24 it_file dbg ( "merger_debug.it" ); 26 25 26 it_assert_debug ( rv.equal ( g0->_rv() ),"Incompatible g0" ); 27 //Empirical density - samples 28 eEmp eSmp ( rv,Ns ); 29 eSmp.set_parameters ( ones ( Ns ), g0 ); 30 Array<vec> &Smp = eSmp._samples(); //aux 31 vec &w = eSmp._w(); //aux 32 33 mat Smp_ex =ones ( rv.count() +1,Ns ); // Extended samples The last row is ones for the ARX model 34 for ( int i=0;i<Ns;i++ ) { set_col_part ( Smp_ex,i,Smp ( i ) );} 35 36 dbg << Name ( "Smp0" ) << Smp_ex; 37 38 // Stuff for merging 27 39 vec lw_src ( Ns ); 28 40 vec lw_mix ( Ns ); 41 mat lW=zeros ( n,Ns ); 29 42 vec vec0 ( 0 ); 30 43 31 //Make sure g0 is compatible 32 it_assert_debug ( rv.equal(g0->_rv()),"Incompatible g0" ); 44 // Initial component in the mixture model 45 mat V0=1e-8*eye ( rv.count() +1 ); 46 ARX A0 ( RV ( "{th_r }", vec_1 ( rv.count() * ( rv.count() +1 ) ) ),\ 47 V0, 4.0 ); //initial guess of Mix: zero mean, large variance 33 48 34 // Initial component in the mixture model from g035 mat V0=1e-8*eye ( rv.count() +1 );36 ARX A0 ( rv, V0, 3.0 ); //initial guess of Mix: zero mean, large variance37 49 38 //generate samples from the proposal g039 for ( int i =0; i<Ns; i++ ) {Smp.set_col ( i,g0->sample() );}40 50 41 //Initialize Mixture model 42 Mix.init ( &A0, Smp, Nc ); 51 // ============= MAIN LOOP ================== 52 bool converged=false; 53 int niter = 0; 54 char str[100]; 55 56 vec Mix_pdf(Ns); 57 while ( !converged ) { 58 //Re-estimate Mix 59 //Re-Initialize Mixture model 60 Mix.init ( &A0, Smp_ex, Nc ); 61 Mix.bayesB ( Smp_ex ); 43 62 44 bool converged=false; 45 while ( !converged ) { 46 Mix.bayesB ( Smp ); 47 //Generate weight for each particle 63 // Generate new samples 64 eSmp.set_samples ( Mix.predictor ( rv ) ); 65 for ( int i=0;i<Ns;i++ ) { set_col_part ( Smp_ex,i,Smp ( i ) );} 66 67 sprintf ( str,"Mpdf%d",niter ); 68 for (int i=0;i<Ns;i++){Mix_pdf(i) = Mix.logpred(Smp_ex.get_col(i));} 69 dbg << Name ( str ) << Mix_pdf; 70 71 sprintf ( str,"Smp%d",niter ); 72 dbg << Name ( str ) << Smp_ex; 73 74 //Importace weighting 48 75 for ( int i=0;i<n;i++ ) { 49 76 //Split according to dependency in rvs … … 51 78 if ( rvsinrv ( i ).length() ==rv.count() ) { 52 79 // no need for conditioning or marginalization 53 lw_src=mpdfs ( i )->_epdf().evalpdflog_m ( Smp ); 54 lw_mix=Mix.logpred_m ( Smp ); 80 for ( int j=0;j<Ns;j++ ) { 81 lw_src ( j ) =mpdfs ( i )->_epdf().evalpdflog ( Smp ( j ) ); 82 } 55 83 } 56 lW.set_row ( i, lw_src -lw_mix );84 lW.set_row ( i, lw_src ); // do not divide by mix 57 85 } 58 converged = true; 86 //Importance of the mixture 87 for ( int j=0;j<Ns;j++ ) { 88 lw_mix ( j ) =Mix.logpred ( Smp_ex.get_col ( j ) ); 89 } 90 sprintf ( str,"lW%d",niter ); 91 dbg << Name ( str ) << lW; 92 93 w = lognorm_merge ( lW ); //merge 94 95 sprintf ( str,"w%d",niter ); 96 dbg << Name ( str ) << w; 97 98 //Importance weighting 99 w /=exp ( lw_mix ); // hoping that it is not numerically sensitive... 100 //renormalize 101 w /=sum ( w ); 102 103 eSmp.resample(); // So that it can be used in bayes 104 for ( int i=0;i<Ns;i++ ) { set_col_part ( Smp_ex,i,Smp ( i ) );} 105 106 sprintf ( str,"Smp_res%d",niter ); 107 dbg << Name ( str ) << Smp; 108 109 // ==== stopping rule === 110 niter++; 111 converged = ( niter>3 ); 59 112 } 60 113 -
bdm/estim/merger.h
r177 r180 41 41 compositepdf ( S ), epdf ( getrv ( false ) ), 42 42 Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ) 43 { beta=2.0; Ns=100; Nc=10;}43 {setindices(rv); beta=2.0; Ns=100; Nc=10;} 44 44 //! Set internal parameters used in approximation 45 45 void set_parameters ( double beta0, int Ns0, int Nc0 ) { beta=beta0;Ns=Ns0;Nc=Nc0;} … … 60 60 //! weight w is a 61 61 vec sample ( )const { return Mix._epdf().sample();} 62 double evalpdflog ( const vec &dt ) const{ return Mix._epdf().evalpdflog ( dt);}62 double evalpdflog ( const vec &dt ) const{ vec dtf=zeros(dt.length()+1); dtf.set_subvector(0,dt); return Mix.logpred ( dtf );} 63 63 vec mean()const {return Mix._epdf().mean();} 64 64 //! for future use -
bdm/estim/mixef.cpp
r176 r180 9 9 Coms.set_size ( c ); 10 10 n=c; 11 weights.set_parameters ( 1/ ( double ) c*ones ( c ) ); //assume at least one observation in each comp.11 weights.set_parameters ( ones ( c ) ); //assume at least one observation in each comp. 12 12 //est will be done at the end 13 13 // … … 30 30 //pick one datum 31 31 int ind=ndat*UniRNG.sample(); 32 Coms ( i )->bayes ( Data.get_col ( ind ), ndat/n);32 Coms ( i )->bayes ( Data.get_col ( ind ),1.0 ); 33 33 } 34 34 … … 110 110 return log ( exLL ); 111 111 } 112 113 emix* MixEF::predictor(const RV &rv){ 114 Array<epdf*> pC(n); 115 for(int i=0;i<n;i++){pC(i)=Coms(i)->predictor(rv);} 116 emix* tmp; 117 tmp = new emix(rv); 118 tmp->set_parameters(weights._epdf().mean(), pC, false); 119 tmp->ownComs(); 120 return tmp; 121 } -
bdm/estim/mixef.h
r177 r180 33 33 The characteristic feature of this model is that if the exact values of the latent variable were known, estimation of the parameters can be handled by a single model. For example, for the case of mixture models, posterior density for each component parameters would be a BayesianModel from Exponential Family. 34 34 35 This class uses EM-style type algorithms for estimation of its parameters. Under this simplification, the posterior density is a product of exponential family members, hence approximate estimation projectthis class itself belongs to the exponential family.35 This class uses EM-style type algorithms for estimation of its parameters. Under this simplification, the posterior density is a product of exponential family members, hence under EM-style approximate estimation this class itself belongs to the exponential family. 36 36 37 37 TODO: Extend BM to use rvc. … … 48 48 eprod* est; 49 49 ////!Indeces of component rvc in common rvc 50 50 51 51 //! Auxiliary function for use in constructors 52 52 void build_est() { … … 91 91 double logpred ( const vec &dt ) const; 92 92 const epdf& _epdf() const {return *est;} 93 emix* predictor(const RV &rv); 94 //! Flatten the density as if it was not estimated from the data 95 void flatten(double sumw=1.0); 93 96 }; 94 97 -
bdm/itpp_ext.cpp
r145 r180 48 48 ov ( iv ( i ) ) = v ( i ); 49 49 } 50 } 51 52 vec get_vec(const vec &v, const ivec &indexlist){ 53 int size = indexlist.size(); 54 vec temp(size); 55 for (int i = 0; i < size; ++i) { 56 temp(i) = v._data()[indexlist(i)]; 57 } 58 return temp; 50 59 } 51 60 -
bdm/itpp_ext.h
r145 r180 22 22 ivec linspace ( int from, int to ); 23 23 24 vec get_vec(const vec &v, const ivec &indexlist); 25 24 26 bvec operator& ( const bvec &a, const bvec &b ); 25 27 bvec operator| ( const bvec &a, const bvec &b ); … … 29 31 30 32 void set_subvector ( vec &ov, const ivec &iv, const vec &v ); 33 34 template<class Num_T> inline 35 void set_col_part(mat &M, int c, const Vec<Num_T> &v) 36 { 37 copy_vector(v.size(), v._data(), M._data() + c*M.rows()); 38 } 31 39 32 40 /*! -
bdm/math/libDC.h
r168 r180 15 15 16 16 #include <itpp/itbase.h> 17 #include "../itpp_ext.h" 17 18 18 19 using namespace itpp; … … 172 173 }; 173 174 174 /*! \brief Matrix stored in LD form, ( typically known as UD)175 /*! \brief Matrix stored in LD form, (commonly known as UD) 175 176 176 177 Matrix is decomposed as follows: \f[M = L'DL\f] where only \f$L\f$ and \f$D\f$ matrices are stored. … … 180 181 { 181 182 public: 182 183 183 //! Construct by copy of L and D. 184 184 ldmat ( const mat &L, const vec &D ); 185 185 //! Construct by decomposition of full matrix V. 186 186 ldmat (const mat &V ); 187 //! Construct by restructuring of V0 accordint to permutation vector perm. 188 ldmat (const ldmat &V0, const ivec &perm):sqmat(V0.rows()){ ldform(V0.L.get_cols(perm), V0.D);}; 187 189 //! Construct diagonal matrix with diagonal D0 188 190 ldmat ( vec D0 ); … … 191 193 //! Default initialization with proper size 192 194 ldmat(const int dim0); 193 195 194 196 //! Destructor for future use; 195 197 virtual ~ldmat(){}; … … 206 208 double qform (const vec &v ) const; 207 209 double invqform (const vec &v ) const; 208 // sqmat& operator -= ( const sqmat & ld2 );209 210 void clear(); 210 211 int cols() const; … … 212 213 vec sqrt_mult ( const vec &v ) const; 213 214 215 214 216 /*! \brief Matrix inversion preserving the chosen form. 215 217 @param Inv a space where the inverse is stored. … … 267 269 268 270 }; 269 270 271 271 272 //////// Operations: