exp_family.h
Go to the documentation of this file.00001 00013 #ifndef EF_H 00014 #define EF_H 00015 00016 00017 #include "../shared_ptr.h" 00018 #include "../base/bdmbase.h" 00019 #include "../math/chmat.h" 00020 00021 namespace bdm { 00022 00023 00025 extern Uniform_RNG UniRNG; 00027 extern Normal_RNG NorRNG; 00029 extern Gamma_RNG GamRNG; 00030 00031 00037 class eEF : public epdf { 00038 public: 00039 // eEF() :epdf() {}; 00041 eEF () : epdf () {}; 00043 virtual double lognc() const = 0; 00044 00046 virtual double evallog_nn ( const vec &val ) const NOT_IMPLEMENTED(0); 00047 00049 virtual double evallog ( const vec &val ) const { 00050 double tmp; 00051 tmp = evallog_nn ( val ) - lognc(); 00052 return tmp; 00053 } 00055 virtual vec evallog_mat ( const mat &Val ) const { 00056 vec x ( Val.cols() ); 00057 for ( int i = 0; i < Val.cols(); i++ ) { 00058 x ( i ) = evallog_nn ( Val.get_col ( i ) ) ; 00059 } 00060 return x - lognc(); 00061 } 00063 virtual vec evallog_mat ( const Array<vec> &Val ) const { 00064 vec x ( Val.length() ); 00065 for ( int i = 0; i < Val.length(); i++ ) { 00066 x ( i ) = evallog_nn ( Val ( i ) ) ; 00067 } 00068 return x - lognc(); 00069 } 00070 00072 virtual void pow ( double p ) NOT_IMPLEMENTED_VOID; 00073 }; 00074 00075 00077 class BMEF : public BM { 00078 public: 00080 double frg; 00081 protected: 00083 double last_lognc; 00085 double frg_sched_factor; 00086 public: 00088 BMEF ( double frg0 = 1.0 ) : BM (), frg ( frg0 ), last_lognc(0.0),frg_sched_factor(0.0) {} 00090 BMEF ( const BMEF &B ) : BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ),frg_sched_factor(B.frg_sched_factor) {} 00092 virtual void set_statistics ( const BMEF* BM0 ) NOT_IMPLEMENTED_VOID; 00093 00095 virtual void bayes_weighted ( const vec &data, const vec &cond = empty_vec, const double w = 1.0 ) { 00096 if (frg_sched_factor>0) { 00097 frg = frg*(1-frg_sched_factor)+frg_sched_factor; 00098 } 00099 }; 00100 //original Bayes 00101 void bayes ( const vec &yt, const vec &cond = empty_vec ); 00102 00104 virtual void flatten ( const BMEF * B, double weight=1.0 ) NOT_IMPLEMENTED_VOID;; 00105 00106 00107 void to_setting ( Setting &set ) const 00108 { 00109 BM::to_setting( set ); 00110 UI::save(frg, set, "frg"); 00111 UI::save( frg_sched_factor, set, "frg_sched_factor" ); 00112 } 00113 00130 void from_setting( const Setting &set) { 00131 BM::from_setting(set); 00132 if ( !UI::get ( frg, set, "frg" ) ) 00133 frg = 1.0; 00134 if ( !UI::get ( frg_sched_factor, set, "frg_sched_factor" ) ) 00135 frg_sched_factor = 0.0; 00136 } 00137 00138 void validate() { 00139 BM::validate(); 00140 } 00141 00142 }; 00143 00149 class mgdirac: public pdf { 00150 protected: 00151 shared_ptr<fnc> g; 00152 public: 00153 vec samplecond(const vec &cond) { 00154 bdm_assert_debug(cond.length()==g->dimensionc(),"given cond in not compatible with g"); 00155 vec tmp = g->eval(cond); 00156 return tmp; 00157 } 00158 double evallogcond ( const vec &yt, const vec &cond ) { 00159 return std::numeric_limits< double >::max(); 00160 } 00161 00171 void from_setting(const Setting& set); 00172 void to_setting(Setting &set) const; 00173 void validate(); 00174 }; 00175 UIREGISTER(mgdirac); 00176 00177 00178 template<class sq_T, template <typename> class TEpdf> 00179 class mlnorm; 00180 00186 template<class sq_T> 00187 class enorm : public eEF { 00188 protected: 00190 vec mu; 00192 sq_T R; 00193 public: 00196 00197 enorm () : eEF (), mu (), R () {}; 00198 enorm ( const vec &mu, const sq_T &R ) { 00199 set_parameters ( mu, R ); 00200 } 00201 void set_parameters ( const vec &mu, const sq_T &R ); 00211 void from_setting ( const Setting &root ); 00212 void to_setting ( Setting &root ) const ; 00213 00214 void validate(); 00216 00219 00221 void dupdate ( mat &v, double nu = 1.0 ); 00222 00224 double bhattacharyya(const enorm<sq_T> &e2) { 00225 bdm_assert(dim == e2.dimension(), "enorms of differnt dimensions"); 00226 sq_T P=R; 00227 P.add(e2._R()); 00228 00229 double tmp = 0.125*P.invqform(mu - e2._mu()) + 0.5*(P.logdet() - 0.5*(R.logdet() + e2._R().logdet())); 00230 return tmp; 00231 } 00232 00233 vec sample() const; 00234 00235 double evallog_nn ( const vec &val ) const; 00236 double lognc () const; 00237 vec mean() const { 00238 return mu; 00239 } 00240 vec variance() const { 00241 return diag ( R.to_mat() ); 00242 } 00243 mat covariance() const { 00244 return R.to_mat(); 00245 } 00246 // mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why? 00247 shared_ptr<pdf> condition ( const RV &rvn ) const; 00248 00249 // target not typed to mlnorm<sq_T, enorm<sq_T> > & 00250 // because that doesn't compile (perhaps because we 00251 // haven't finished defining enorm yet), but the type 00252 // is required 00253 void condition ( const RV &rvn, pdf &target ) const; 00254 00255 shared_ptr<epdf> marginal ( const RV &rvn ) const; 00256 void marginal ( const RV &rvn, enorm<sq_T> &target ) const; 00258 00261 00262 vec& _mu() { 00263 return mu; 00264 } 00265 const vec& _mu() const { 00266 return mu; 00267 } 00268 void set_mu ( const vec mu0 ) { 00269 mu = mu0; 00270 } 00271 sq_T& _R() { 00272 return R; 00273 } 00274 const sq_T& _R() const { 00275 return R; 00276 } 00278 00279 }; 00280 UIREGISTER2 ( enorm, chmat ); 00281 SHAREDPTR2 ( enorm, chmat ); 00282 UIREGISTER2 ( enorm, ldmat ); 00283 SHAREDPTR2 ( enorm, ldmat ); 00284 UIREGISTER2 ( enorm, fsqmat ); 00285 SHAREDPTR2 ( enorm, fsqmat ); 00286 00289 typedef enorm<ldmat> egauss; 00290 UIREGISTER(egauss); 00291 00292 00293 //forward declaration 00294 class mstudent; 00295 00300 template<class sq_T> 00301 class estudent : public eEF { 00302 protected: 00304 vec mu; 00306 sq_T H; 00308 double delta; 00309 public: 00310 double evallog_nn(const vec &val) const { 00311 double tmp = -0.5*H.logdet() - 0.5*(delta + dim) * log(1+ H.invqform(val - mu)/delta); 00312 return tmp; 00313 } 00314 double lognc() const { 00315 //log(pi) = 1.14472988584940 00316 double tmp = -lgamma(0.5*(delta+dim))+lgamma(0.5*delta) + 0.5*dim*(log(delta) + 1.14472988584940); 00317 return tmp; 00318 } 00319 void marginal (const RV &rvm, estudent<sq_T> &marg) const { 00320 ivec ind = rvm.findself_ids(rv); // indices of rvm in rv 00321 marg._mu() = mu(ind); 00322 marg._H() = sq_T(H,ind); 00323 marg._delta() = delta; 00324 marg.validate(); 00325 } 00326 shared_ptr<epdf> marginal(const RV &rvm) const { 00327 shared_ptr<estudent<sq_T> > tmp = new estudent<sq_T>; 00328 marginal(rvm, *tmp); 00329 return tmp; 00330 } 00331 vec sample() const { 00332 enorm<sq_T> en; 00333 en._mu()=mu; 00334 en._R()=H; 00335 en._R()*=delta/(delta-2); 00336 en.validate(); 00337 return en.sample(); 00338 } 00339 00340 vec mean() const { 00341 return mu; 00342 } 00343 mat covariance() const { 00344 return delta/(delta-2)*H.to_mat(); 00345 } 00346 vec variance() const { 00347 return diag(covariance()); 00348 } 00352 vec& _mu() { 00353 return mu; 00354 } 00356 sq_T& _H() { 00357 return H; 00358 } 00360 double& _delta() { 00361 return delta; 00362 } 00365 void from_setting(const Setting &set) { 00366 epdf::from_setting(set); 00367 mat H0; 00368 UI::get(H0,set, "H"); 00369 H= H0; // conversion!! 00370 UI::get(delta,set,"delta"); 00371 UI::get(mu,set,"mu"); 00372 } 00373 void to_setting(Setting &set) const { 00374 epdf::to_setting(set); 00375 UI::save(H.to_mat(), set, "H"); 00376 UI::save(delta, set, "delta"); 00377 UI::save(mu, set, "mu"); 00378 } 00379 void validate() { 00380 eEF::validate(); 00381 dim = H.rows(); 00382 } 00383 }; 00384 UIREGISTER2(estudent,fsqmat); 00385 UIREGISTER2(estudent,ldmat); 00386 UIREGISTER2(estudent,chmat); 00387 00388 // template < class sq_T, template <typename> class TEpdf = estudent > 00389 // class mstudent : public pdf_internal< TEpdf<sq_T> > { 00390 // protected: 00391 // //! Internal epdf that arise by conditioning on \c rvc 00392 // mat A; 00393 // //! Constant additive term 00394 // vec mu_const; 00395 // 00396 // public: 00397 // void condition( const vec &cond ) { 00398 // iepdf._mu()=A*val*mu_const; 00399 // } 00400 // 00401 // }; 00402 00416 class egiw : public eEF { 00419 00420 LOG_LEVEL(egiw,logvartheta); 00421 00422 protected: 00424 ldmat V; 00426 double nu; 00428 int dimx; 00430 int nPsi; 00431 public: 00434 egiw() : eEF(),dimx(0) {}; 00435 egiw ( int dimx0, ldmat V0, double nu0 = -1.0 ) : eEF(),dimx(0) { 00436 set_parameters ( dimx0, V0, nu0 ); 00437 validate(); 00438 }; 00439 00440 void set_parameters ( int dimx0, ldmat V0, double nu0 = -1.0 ); 00442 00443 vec sample() const; 00444 mat sample_mat ( int n ) const; 00445 vec mean() const; 00446 vec variance() const; 00447 //mat covariance() const; 00448 void sample_mat ( mat &Mi, chmat &Ri ) const; 00449 00450 void factorize ( mat &M, ldmat &Vz, ldmat &Lam ) const; 00452 vec est_theta() const; 00453 00455 ldmat est_theta_cov() const; 00456 00458 void mean_mat ( mat &M, mat&R ) const; 00460 double evallog_nn ( const vec &val ) const; 00461 double lognc () const; 00462 void pow ( double p ) { 00463 V *= p; 00464 nu *= p; 00465 }; 00466 00468 shared_ptr<epdf> marginal(const RV &rvm) const { 00469 bdm_assert(dimx==1, "Not supported"); 00470 //TODO - this is too trivial!!! 00471 ivec ind = rvm.findself_ids(rv); 00472 if (min(ind)==0) { //assume it si 00473 shared_ptr<estudent<ldmat> > tmp = new estudent<ldmat>; 00474 mat M; 00475 ldmat Vz; 00476 ldmat Lam; 00477 factorize(M,Vz,Lam); 00478 00479 tmp->_mu() = M.get_col(0); 00480 ldmat H; 00481 Vz.inv(H); 00482 H *=Lam._D()(0)/nu; 00483 tmp->_H() = H; 00484 tmp->_delta() = nu; 00485 tmp->validate(); 00486 return tmp; 00487 } 00488 return NULL; 00489 } 00492 00493 ldmat& _V() { 00494 return V; 00495 } 00496 const ldmat& _V() const { 00497 return V; 00498 } 00499 double& _nu() { 00500 return nu; 00501 } 00502 const double& _nu() const { 00503 return nu; 00504 } 00505 const int & _dimx() const { 00506 return dimx; 00507 } 00508 00534 void from_setting ( const Setting &set ); 00536 void to_setting ( Setting& set ) const; 00537 void validate(); 00538 void log_register ( bdm::logger& L, const string& prefix ); 00539 00540 void log_write() const; 00542 }; 00543 UIREGISTER ( egiw ); 00544 SHAREDPTR ( egiw ); 00545 00549 template<class sq_T> 00550 class egw_ls: public eEF{ 00551 public: 00552 vec theta; 00553 sq_T P; 00554 double omega; 00555 double nu; 00556 00557 vec mean() const{ 00558 return concat(theta, omega); 00559 } 00560 mat covariance() const { 00561 sq_T tmp=P; 00562 tmp*=nu/((nu-2)*omega); 00563 return tmp.to_mat();//<======= error - missing omega 00564 } 00565 vec variance() const { 00566 return diag(covariance());//<======= error - missing omega 00567 } 00568 vec sample() const NOT_IMPLEMENTED(vec(0)); 00569 double lognc() const {return 0.0;} //TODO 00570 }; 00571 00580 class eDirich: public eEF { 00581 protected: 00583 vec beta; 00584 public: 00587 00588 eDirich () : eEF () {}; 00589 eDirich ( const eDirich &D0 ) : eEF () { 00590 set_parameters ( D0.beta ); 00591 validate(); 00592 }; 00593 eDirich ( const vec &beta0 ) { 00594 set_parameters ( beta0 ); 00595 validate(); 00596 }; 00597 void set_parameters ( const vec &beta0 ) { 00598 beta = beta0; 00599 dim = beta.length(); 00600 } 00602 00604 vec sample() const { 00605 vec y ( beta.length() ); 00606 for ( int i = 0; i < beta.length(); i++ ) { 00607 GamRNG.setup ( beta ( i ), 1 ); 00608 #pragma omp critical 00609 y ( i ) = GamRNG(); 00610 } 00611 return y / sum ( y ); 00612 } 00613 00614 vec mean() const { 00615 return beta / sum ( beta ); 00616 }; 00617 vec variance() const { 00618 double gamma = sum ( beta ); 00619 return elem_mult ( beta, ( gamma - beta ) ) / ( gamma*gamma* ( gamma + 1 ) ); 00620 } 00622 double evallog_nn ( const vec &val ) const { 00623 double tmp; 00624 tmp = ( beta - 1 ) * log ( val ); 00625 return tmp; 00626 } 00627 00628 double lognc () const { 00629 double tmp; 00630 double gam = sum ( beta ); 00631 double lgb = 0.0; 00632 for ( int i = 0; i < beta.length(); i++ ) { 00633 lgb += lgamma ( beta ( i ) ); 00634 } 00635 tmp = lgb - lgamma ( gam ); 00636 return tmp; 00637 } 00638 00640 vec& _beta() { 00641 return beta; 00642 } 00643 00652 void from_setting ( const Setting &set ); 00653 00654 void validate(); 00655 00656 void to_setting ( Setting &set ) const; 00657 }; 00658 UIREGISTER ( eDirich ); 00659 00668 class eBeta: public eEF { 00669 public: 00671 vec alpha; 00673 vec beta; 00674 public: 00677 00678 eBeta () : eEF () {}; 00679 eBeta ( const eBeta &B0 ) : eEF (), alpha(B0.alpha),beta(B0.beta) { 00680 validate(); 00681 }; 00683 00685 vec sample() const { 00686 vec y ( beta.length() ); // all vectors 00687 for ( int i = 0; i < beta.length(); i++ ) { 00688 GamRNG.setup ( alpha ( i ), 1 ); 00689 #pragma omp critical 00690 double Ga = GamRNG(); 00691 00692 GamRNG.setup ( beta ( i ), 1 ); 00693 #pragma omp critical 00694 double Gb = GamRNG(); 00695 00696 y ( i ) = Ga/(Ga+Gb); 00697 } 00698 return y; 00699 } 00700 00701 vec mean() const { 00702 return elem_div(alpha, alpha + beta); // dot-division 00703 }; 00704 vec variance() const { 00705 vec apb=alpha+beta; 00706 return elem_div (elem_mult ( alpha, beta) , 00707 elem_mult ( elem_mult(apb,apb), apb+1 ) ); 00708 } 00710 double evallog_nn ( const vec &val ) const { 00711 double tmp; 00712 tmp = ( alpha - 1 ) * log ( val ) + (beta-1)*log(1-val); 00713 return tmp; 00714 } 00715 00716 double lognc () const { 00717 double lgb = 0.0; 00718 for ( int i = 0; i < beta.length(); i++ ) { 00719 lgb += -lgamma ( alpha(i)+beta(i) ) + lgamma(alpha(i)) + lgamma(beta(i)); 00720 } 00721 return lgb; 00722 } 00723 00734 void from_setting ( const Setting &set ) { 00735 UI::get(alpha, set, "alpha", UI::compulsory); 00736 UI::get(beta, set, "beta", UI::compulsory); 00737 } 00738 00739 void validate() { 00740 bdm_assert(alpha.length()==beta.length(), "eBeta:: alpha and beta length do not match"); 00741 dim = alpha.length(); 00742 } 00743 00744 void to_setting ( Setting &set ) const { 00745 UI::save(alpha, set, "alpha"); 00746 UI::save(beta, set, "beta"); 00747 } 00748 }; 00749 UIREGISTER ( eBeta ); 00750 00762 class mDirich: public pdf_internal<eDirich> { 00763 protected: 00765 double k; 00767 vec &_beta; 00769 vec betac; 00770 public: 00771 mDirich() : pdf_internal<eDirich>(), _beta ( iepdf._beta() ) {}; 00772 void condition ( const vec &val ) { 00773 _beta = val / k + betac; 00774 }; 00775 00794 void from_setting ( const Setting &set ); 00795 void to_setting (Setting &set) const; 00796 void validate(); 00797 }; 00798 UIREGISTER ( mDirich ); 00799 00814 class mBeta: public pdf_internal<eBeta> { 00816 vec k; 00818 vec betac; 00819 00820 public: 00821 void condition ( const vec &val ) { 00822 this->iepdf.alpha = elem_div(val , k) + betac; 00823 this->iepdf.beta = elem_div (1-val , k) + betac; 00824 }; 00825 00845 void from_setting ( const Setting &set ); 00846 00847 void to_setting (Setting &set) const; 00848 00849 void validate() { 00850 pdf_internal<eBeta>::validate(); 00851 bdm_assert(betac.length()==dimension(),"Incomaptible betac"); 00852 bdm_assert(k.length()==dimension(),"Incomaptible k"); 00853 dimc = iepdf.dimension(); 00854 } 00856 }; 00857 UIREGISTER(mBeta); 00858 00860 class multiBM : public BMEF { 00861 protected: 00863 eDirich est; 00865 vec β 00866 public: 00868 multiBM () : BMEF (), est (), beta ( est._beta() ) { 00869 if ( beta.length() > 0 ) { 00870 last_lognc = est.lognc(); 00871 } else { 00872 last_lognc = 0.0; 00873 } 00874 } 00876 multiBM ( const multiBM &B ) : BMEF ( B ), est ( B.est ), beta ( est._beta() ) {} 00878 void set_statistics ( const BM* mB0 ) { 00879 const multiBM* mB = dynamic_cast<const multiBM*> ( mB0 ); 00880 beta = mB->beta; 00881 } 00882 void bayes ( const vec &yt, const vec &cond = empty_vec ); 00883 00884 double logpred ( const vec &yt ) const; 00885 00886 void flatten ( const BMEF* B , double weight); 00887 00889 const eDirich& posterior() const { 00890 return est; 00891 }; 00893 void set_parameters ( const vec &beta0 ) { 00894 est.set_parameters ( beta0 ); 00895 est.validate(); 00896 if ( evalll ) { 00897 last_lognc = est.lognc(); 00898 } 00899 } 00900 00901 void to_setting ( Setting &set ) const { 00902 BMEF::to_setting ( set ); 00903 UI::save( &est, set, "prior" ); 00904 } 00905 00915 void from_setting (const Setting &set ) { 00916 BMEF::from_setting ( set ); 00917 UI::get( est, set, "prior" ); 00918 } 00919 }; 00920 UIREGISTER( multiBM ); 00921 00931 class egamma : public eEF { 00932 protected: 00934 vec alpha; 00936 vec beta; 00937 public : 00940 egamma () : eEF (), alpha ( 0 ), beta ( 0 ) {}; 00941 egamma ( const vec &a, const vec &b ) { 00942 set_parameters ( a, b ); 00943 validate(); 00944 }; 00945 void set_parameters ( const vec &a, const vec &b ) { 00946 alpha = a, beta = b; 00947 }; 00949 00950 vec sample() const; 00951 double evallog_nn ( const vec &val ) const; 00952 double lognc () const; 00954 vec& _alpha() { 00955 return alpha; 00956 } 00958 vec& _beta() { 00959 return beta; 00960 } 00961 vec mean() const { 00962 return elem_div ( alpha, beta ); 00963 } 00964 vec variance() const { 00965 return elem_div ( alpha, elem_mult ( beta, beta ) ); 00966 } 00967 00979 void from_setting ( const Setting &set ); 00980 00981 void to_setting ( Setting &set ) const; 00982 void validate(); 00983 }; 00984 UIREGISTER ( egamma ); 00985 SHAREDPTR ( egamma ); 00986 01003 class eigamma : public egamma { 01004 protected: 01005 public : 01010 01011 vec sample() const { 01012 return 1.0 / egamma::sample(); 01013 }; 01015 vec mean() const { 01016 return elem_div ( beta, alpha - 1 ); 01017 } 01018 vec variance() const { 01019 vec mea = mean(); 01020 return elem_div ( elem_mult ( mea, mea ), alpha - 2 ); 01021 } 01022 }; 01023 /* 01025 class emix : public epdf { 01026 protected: 01027 int n; 01028 vec &w; 01029 Array<epdf*> Coms; 01030 public: 01032 emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {}; 01033 void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;} 01034 vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;}; 01035 }; 01036 */ 01037 01039 class euni: public epdf { 01040 protected: 01042 vec low; 01044 vec high; 01046 vec distance; 01048 double nk; 01050 double lnk; 01051 public: 01054 euni () : epdf () {} 01055 euni ( const vec &low0, const vec &high0 ) { 01056 set_parameters ( low0, high0 ); 01057 } 01058 void set_parameters ( const vec &low0, const vec &high0 ) { 01059 distance = high0 - low0; 01060 low = low0; 01061 high = high0; 01062 nk = prod ( 1.0 / distance ); 01063 lnk = log ( nk ); 01064 } 01066 01067 double evallog ( const vec &val ) const { 01068 if ( any ( val < low ) && any ( val > high ) ) { 01069 return -inf; 01070 } else return lnk; 01071 } 01072 vec sample() const { 01073 vec smp ( dim ); 01074 #pragma omp critical 01075 UniRNG.sample_vector ( dim , smp ); 01076 return low + elem_mult ( distance, smp ); 01077 } 01079 vec mean() const { 01080 return ( high - low ) / 2.0; 01081 } 01082 vec variance() const { 01083 return ( pow ( high, 2 ) + pow ( low, 2 ) + elem_mult ( high, low ) ) / 3.0; 01084 } 01085 01086 01100 void from_setting ( const Setting &set ); 01101 void to_setting (Setting &set) const; 01102 void validate(); 01103 }; 01104 UIREGISTER ( euni ); 01105 01107 class mguni : public pdf_internal<euni> { 01109 shared_ptr<fnc> mean; 01111 vec delta; 01112 public: 01113 void condition ( const vec &cond ) { 01114 vec mea = mean->eval ( cond ); 01115 iepdf.set_parameters ( mea - delta, mea + delta ); 01116 } 01117 01127 void from_setting ( const Setting &set ) { 01128 pdf::from_setting ( set ); //reads rv and rvc 01129 UI::get ( delta, set, "delta", UI::compulsory ); 01130 mean = UI::build<fnc> ( set, "mean", UI::compulsory ); 01131 iepdf.set_parameters ( -delta, delta ); 01132 } 01133 void to_setting (Setting &set) const { 01134 pdf::to_setting ( set ); 01135 UI::save( iepdf.mean(), set, "delta"); 01136 UI::save(mean, set, "mean"); 01137 } 01138 void validate() { 01139 pdf_internal<euni>::validate(); 01140 dimc = mean->dimensionc(); 01141 01142 } 01143 01144 }; 01145 UIREGISTER ( mguni ); 01151 template < class sq_T, template <typename> class TEpdf = enorm > 01152 class mlnorm : public pdf_internal< TEpdf<sq_T> > { 01153 protected: 01155 mat A; 01157 vec mu_const; 01158 // vec& _mu; //cached epdf.mu; !!!!!! WHY NOT? 01159 public: 01162 mlnorm() : pdf_internal< TEpdf<sq_T> >() {}; 01163 mlnorm ( const mat &A, const vec &mu0, const sq_T &R ) : pdf_internal< TEpdf<sq_T> >() { 01164 set_parameters ( A, mu0, R ); 01165 validate(); 01166 } 01167 01169 void set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 ) { 01170 this->iepdf.set_parameters ( zeros ( A0.rows() ), R0 ); 01171 A = A0; 01172 mu_const = mu0; 01173 } 01174 01177 void condition ( const vec &cond ) { 01178 this->iepdf._mu() = A * cond + mu_const; 01179 //R is already assigned; 01180 } 01181 01183 const vec& _mu_const() const { 01184 return mu_const; 01185 } 01187 const mat& _A() const { 01188 return A; 01189 } 01191 mat _R() const { 01192 return this->iepdf._R().to_mat(); 01193 } 01195 sq_T __R() const { 01196 return this->iepdf._R(); 01197 } 01198 01200 template<typename sq_M> 01201 friend std::ostream &operator<< ( std::ostream &os, mlnorm<sq_M, enorm> &ml ); 01202 01214 void from_setting ( const Setting &set ) { 01215 pdf::from_setting ( set ); 01216 01217 UI::get ( A, set, "A", UI::compulsory ); 01218 UI::get ( mu_const, set, "const", UI::optional); 01219 mat R0; 01220 UI::get ( R0, set, "R", UI::compulsory ); 01221 set_parameters ( A, mu_const, R0 ); 01222 } 01223 01224 void to_setting (Setting &set) const { 01225 pdf::to_setting(set); 01226 UI::save ( A, set, "A"); 01227 UI::save ( mu_const, set, "const"); 01228 UI::save ( _R(), set, "R"); 01229 } 01230 01231 void validate() { 01232 pdf_internal<TEpdf<sq_T> >::validate(); 01233 if (mu_const.length()==0) { // default in from_setting 01234 mu_const=zeros(A.rows()); 01235 } 01236 bdm_assert ( A.rows() == mu_const.length(), "mlnorm: A vs. mu mismatch" ); 01237 bdm_assert ( A.rows() == _R().rows(), "mlnorm: A vs. R mismatch" ); 01238 this->dimc = A.cols(); 01239 01240 } 01241 }; 01242 UIREGISTER2 ( mlnorm, ldmat ); 01243 SHAREDPTR2 ( mlnorm, ldmat ); 01244 UIREGISTER2 ( mlnorm, fsqmat ); 01245 SHAREDPTR2 ( mlnorm, fsqmat ); 01246 UIREGISTER2 ( mlnorm, chmat ); 01247 SHAREDPTR2 ( mlnorm, chmat ); 01248 01251 typedef mlnorm<fsqmat> mlgauss; 01252 UIREGISTER(mlgauss); 01253 01255 template<class sq_T> 01256 class mgnorm : public pdf_internal< enorm< sq_T > > { 01257 private: 01258 // vec μ WHY NOT? 01259 shared_ptr<fnc> g; 01260 01261 public: 01263 mgnorm() : pdf_internal<enorm<sq_T> >() { } 01265 inline void set_parameters ( const shared_ptr<fnc> &g0, const sq_T &R0 ); 01266 inline void condition ( const vec &cond ); 01267 01268 01288 void from_setting ( const Setting &set ) { 01289 pdf::from_setting ( set ); 01290 shared_ptr<fnc> g = UI::build<fnc> ( set, "g", UI::compulsory ); 01291 01292 mat R; 01293 vec dR; 01294 if ( UI::get ( dR, set, "dR" ) ) 01295 R = diag ( dR ); 01296 else 01297 UI::get ( R, set, "R", UI::compulsory ); 01298 01299 set_parameters ( g, R ); 01300 //validate(); 01301 } 01302 01303 01304 void to_setting (Setting &set) const { 01305 UI::save( g,set, "g"); 01306 UI::save(this->iepdf._R().to_mat(),set, "R"); 01307 01308 } 01309 01310 01311 01312 void validate() { 01313 this->iepdf.validate(); 01314 bdm_assert ( g->dimension() == this->iepdf.dimension(), "incompatible function" ); 01315 this->dim = g->dimension(); 01316 this->dimc = g->dimensionc(); 01317 this->iepdf.validate(); 01318 } 01319 01320 }; 01321 01322 UIREGISTER2 ( mgnorm, chmat ); 01323 UIREGISTER2 ( mgnorm, ldmat ); 01324 SHAREDPTR2 ( mgnorm, chmat ); 01325 01326 01334 class mlstudent : public mlnorm<ldmat, enorm> { 01335 protected: 01337 ldmat Lambda; 01339 ldmat &_R; 01341 ldmat Re; 01342 public: 01343 mlstudent () : mlnorm<ldmat, enorm> (), 01344 Lambda (), _R ( iepdf._R() ) {} 01346 void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 ) { 01347 iepdf.set_parameters ( mu0, R0 );// was Lambda, why? 01348 A = A0; 01349 mu_const = mu0; 01350 Re = R0; 01351 Lambda = Lambda0; 01352 } 01353 01354 void condition ( const vec &cond ); 01355 01356 void validate() { 01357 mlnorm<ldmat, enorm>::validate(); 01358 bdm_assert ( A.rows() == mu_const.length(), "mlstudent: A vs. mu mismatch" ); 01359 bdm_assert ( _R.rows() == A.rows(), "mlstudent: A vs. R mismatch" ); 01360 01361 } 01362 }; 01363 01373 class mgamma : public pdf_internal<egamma> { 01374 protected: 01375 01377 double k; 01378 01380 vec &_beta; 01381 01382 public: 01384 mgamma() : pdf_internal<egamma>(), k ( 0 ), 01385 _beta ( iepdf._beta() ) { 01386 } 01387 01389 void set_parameters ( double k, const vec &beta0 ); 01390 01391 void condition ( const vec &val ) { 01392 _beta = k / val; 01393 }; 01394 01405 void from_setting ( const Setting &set ); 01406 void to_setting (Setting &set) const; 01407 void validate(); 01408 }; 01409 UIREGISTER ( mgamma ); 01410 SHAREDPTR ( mgamma ); 01411 01421 class migamma : public pdf_internal<eigamma> { 01422 protected: 01424 double k; 01425 01427 vec &_alpha; 01428 01430 vec &_beta; 01431 01432 public: 01435 migamma() : pdf_internal<eigamma>(), 01436 k ( 0 ), 01437 _alpha ( iepdf._alpha() ), 01438 _beta ( iepdf._beta() ) { 01439 } 01440 01441 migamma ( const migamma &m ) : pdf_internal<eigamma>(), 01442 k ( 0 ), 01443 _alpha ( iepdf._alpha() ), 01444 _beta ( iepdf._beta() ) { 01445 } 01447 01449 void set_parameters ( int len, double k0 ) { 01450 k = k0; 01451 iepdf.set_parameters ( ( 1.0 / ( k*k ) + 2.0 ) *ones ( len ) /*alpha*/, ones ( len ) /*beta*/ ); 01452 }; 01453 01454 void validate () { 01455 pdf_internal<eigamma>::validate(); 01456 dimc = dimension(); 01457 }; 01458 01459 void condition ( const vec &val ) { 01460 _beta = elem_mult ( val, ( _alpha - 1.0 ) ); 01461 }; 01462 }; 01463 01475 class mgamma_fix : public mgamma { 01476 protected: 01478 double l; 01480 vec refl; 01481 public: 01483 mgamma_fix () : mgamma (), refl () {}; 01485 void set_parameters ( double k0 , vec ref0, double l0 ) { 01486 mgamma::set_parameters ( k0, ref0 ); 01487 refl = pow ( ref0, 1.0 - l0 ); 01488 l = l0; 01489 }; 01490 01491 void validate () { 01492 mgamma::validate(); 01493 dimc = dimension(); 01494 }; 01495 01496 void condition ( const vec &val ) { 01497 vec mean = elem_mult ( refl, pow ( val, l ) ); 01498 _beta = k / mean; 01499 }; 01500 }; 01501 01502 01515 class migamma_ref : public migamma { 01516 protected: 01518 double l; 01520 vec refl; 01521 public: 01523 migamma_ref () : migamma (), refl () {}; 01524 01526 void set_parameters ( double k0 , vec ref0, double l0 ) { 01527 migamma::set_parameters ( ref0.length(), k0 ); 01528 refl = pow ( ref0, 1.0 - l0 ); 01529 l = l0; 01530 }; 01531 01532 void validate() { 01533 migamma::validate(); 01534 dimc = dimension(); 01535 }; 01536 01537 void condition ( const vec &val ) { 01538 vec mean = elem_mult ( refl, pow ( val, l ) ); 01539 migamma::condition ( mean ); 01540 }; 01541 01553 void from_setting ( const Setting &set ); 01554 01555 void to_setting (Setting &set) const; 01556 }; 01557 01558 01559 UIREGISTER ( migamma_ref ); 01560 SHAREDPTR ( migamma_ref ); 01561 01571 class elognorm: public enorm<ldmat> { 01572 public: 01573 vec sample() const { 01574 return exp ( enorm<ldmat>::sample() ); 01575 }; 01576 vec mean() const { 01577 vec var = enorm<ldmat>::variance(); 01578 return exp ( mu - 0.5*var ); 01579 }; 01580 01581 }; 01582 01589 class mlognorm : public pdf_internal<elognorm> { 01590 protected: 01592 double sig2; 01593 01595 vec μ 01596 public: 01598 mlognorm() : pdf_internal<elognorm>(), 01599 sig2 ( 0 ), 01600 mu ( iepdf._mu() ) { 01601 } 01602 01604 void set_parameters ( int size, double k ) { 01605 sig2 = 0.5 * log ( k * k + 1 ); 01606 iepdf.set_parameters ( zeros ( size ), 2*sig2*eye ( size ) ); 01607 }; 01608 01609 void validate() { 01610 pdf_internal<elognorm>::validate(); 01611 dimc = iepdf.dimension(); 01612 } 01613 01614 void condition ( const vec &val ) { 01615 mu = log ( val ) - sig2;//elem_mult ( refl,pow ( val,l ) ); 01616 }; 01617 01628 void from_setting ( const Setting &set ); 01629 01630 void to_setting (Setting &set) const; 01631 }; 01632 01633 UIREGISTER ( mlognorm ); 01634 SHAREDPTR ( mlognorm ); 01635 01640 class eWishartCh : public epdf { 01641 protected: 01643 chmat Sigma; 01645 int p; 01647 double delta; 01648 public: 01650 void set_parameters ( const mat &Y0, const double delta0 ) { 01651 Sigma = chmat ( Y0 ); 01652 delta = delta0; 01653 p = Sigma.rows(); 01654 } 01656 void set_parameters ( const chmat &Y0, const double delta0 ) { 01657 Sigma = Y0; 01658 delta = delta0; 01659 p = Sigma.rows(); 01660 } 01661 01662 virtual void validate () { 01663 epdf::validate(); 01664 dim = p * p; 01665 } 01666 01669 chmat sample_mat() const { 01670 mat A_T = zeros ( p, p ); // A transpose 01671 01672 //sample diagonal 01673 for ( int i = 0; i < p; i++ ) { 01674 GamRNG.setup ( 0.5* ( delta - i ) , 0.5 ); // no +1 !! index if from 0 01675 #pragma omp critical 01676 A_T ( i, i ) = sqrt ( GamRNG() ); 01677 } 01678 //do the rest 01679 for ( int i = 0; i < p; i++ ) { 01680 for ( int j = i + 1; j < p; j++ ) { 01681 #pragma omp critical 01682 A_T ( i, j ) = NorRNG.sample(); 01683 } 01684 } 01685 chmat tmp; 01686 tmp._Ch()=A_T*Sigma._Ch(); 01687 return tmp;// return upper triangular part of the decomposition 01688 } 01689 01690 vec sample () const { 01691 return vec ( sample_mat().to_mat()._data(), p*p ); 01692 } 01693 01694 virtual vec mean() const NOT_IMPLEMENTED(0); 01695 01697 virtual vec variance() const NOT_IMPLEMENTED(0); 01698 01699 virtual double evallog ( const vec &val ) const NOT_IMPLEMENTED(0); 01700 01702 void setY ( const mat &Ch0 ) { 01703 copy_vector ( dim, Ch0._data(), Sigma._Ch()._data() ); 01704 } 01705 01707 void _setY ( const vec &ch0 ) { 01708 copy_vector ( dim, ch0._data(), Sigma._Ch()._data() ); 01709 } 01710 01712 const chmat& getY() const { 01713 return Sigma; 01714 } 01715 }; 01716 01718 01720 class eiWishartCh: public epdf { 01721 protected: 01723 eWishartCh W; 01725 int p; 01727 double delta; 01728 public: 01730 void set_parameters ( const mat &Y0, const double delta0 ) { 01731 delta = delta0; 01732 W.set_parameters ( inv ( Y0 ), delta0 ); 01733 p = Y0.rows(); 01734 } 01735 01736 virtual void validate () { 01737 epdf::validate(); 01738 W.validate(); 01739 dim = W.dimension(); 01740 } 01741 01742 01743 vec sample() const { 01744 mat iCh; 01745 iCh = inv ( W.sample_mat()._Ch() ); 01746 return vec ( iCh._data(), dim ); 01747 } 01749 void _setY ( const vec &y0 ) { 01750 mat Ch ( p, p ); 01751 mat iCh ( p, p ); 01752 copy_vector ( dim, y0._data(), Ch._data() ); 01753 01754 iCh = inv ( Ch ); 01755 W.setY ( iCh ); 01756 } 01757 01758 virtual double evallog ( const vec &val ) const { 01759 chmat X ( p ); 01760 const chmat& Y = W.getY(); 01761 01762 copy_vector ( p*p, val._data(), X._Ch()._data() ); 01763 chmat iX ( p ); 01764 X.inv ( iX ); 01765 // compute 01766 // \frac{ |\Psi|^{m/2}|X|^{-(m+p+1)/2}e^{-tr(\Psi X^{-1})/2} }{ 2^{mp/2}\Gamma_p(m/2)}, 01767 mat M = Y.to_mat() * iX.to_mat(); 01768 01769 double log1 = 0.5 * p * ( 2 * Y.logdet() ) - 0.5 * ( delta + p + 1 ) * ( 2 * X.logdet() ) - 0.5 * trace ( M ); 01770 //Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!! 01771 01772 /* if (0) { 01773 mat XX=X.to_mat(); 01774 mat YY=Y.to_mat(); 01775 01776 double log2 = 0.5*p*log(det(YY))-0.5*(delta+p+1)*log(det(XX))-0.5*trace(YY*inv(XX)); 01777 cout << log1 << "," << log2 << endl; 01778 }*/ 01779 return log1; 01780 }; 01781 01782 virtual vec mean() const NOT_IMPLEMENTED(0); 01783 01785 virtual vec variance() const NOT_IMPLEMENTED(0); 01786 }; 01787 01789 class rwiWishartCh : public pdf_internal<eiWishartCh> { 01790 protected: 01792 double sqd; 01794 vec refl; 01796 double l; 01798 int p; 01799 01800 public: 01801 rwiWishartCh() : sqd ( 0 ), l ( 0 ), p ( 0 ) {} 01803 void set_parameters ( int p0, double k, vec ref0, double l0 ) { 01804 p = p0; 01805 double delta = 2 / ( k * k ) + p + 3; 01806 sqd = sqrt ( delta - p - 1 ); 01807 l = l0; 01808 refl = pow ( ref0, 1 - l ); 01809 iepdf.set_parameters ( eye ( p ), delta ); 01810 }; 01811 01812 void validate() { 01813 pdf_internal<eiWishartCh>::validate(); 01814 dimc = iepdf.dimension(); 01815 } 01816 01817 void condition ( const vec &c ) { 01818 vec z = c; 01819 int ri = 0; 01820 for ( int i = 0; i < p*p; i += ( p + 1 ) ) {//trace diagonal element 01821 z ( i ) = pow ( z ( i ), l ) * refl ( ri ); 01822 ri++; 01823 } 01824 01825 iepdf._setY ( sqd*z ); 01826 } 01827 }; 01828 01830 enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 }; 01831 01833 void resample(const vec &w, ivec &ind, RESAMPLING_METHOD=SYSTEMATIC); 01834 01839 class eEmp: public epdf { 01840 protected : 01842 int n; 01844 vec w; 01846 Array<vec> samples; 01847 public: 01850 eEmp () : epdf (), w (), samples () {}; 01852 eEmp ( const eEmp &e ) : epdf ( e ), w ( e.w ), samples ( e.samples ) {}; 01854 01856 void set_statistics ( const vec &w0, const epdf &pdf0 ); 01858 void set_statistics ( const epdf &pdf0 , int n ) { 01859 set_statistics ( ones ( n ) / n, pdf0 ); 01860 }; 01862 void set_samples ( const epdf* pdf0 ); 01864 void set_parameters ( int n0, bool copy = true ) { 01865 n = n0; 01866 w.set_size ( n0, copy ); 01867 samples.set_size ( n0, copy ); 01868 }; 01870 void set_parameters ( const Array<vec> &Av ) { 01871 n = Av.size(); 01872 w = 1 / n * ones ( n ); 01873 samples = Av; 01874 }; 01875 virtual void validate (); 01877 vec& _w() { 01878 return w; 01879 }; 01881 const vec& _w() const { 01882 return w; 01883 }; 01885 Array<vec>& _samples() { 01886 return samples; 01887 }; 01889 const vec& _sample ( int i ) const { 01890 return samples ( i ); 01891 }; 01893 const Array<vec>& _samples() const { 01894 return samples; 01895 }; 01897 void resample ( RESAMPLING_METHOD method = SYSTEMATIC ); 01898 01900 vec sample() const NOT_IMPLEMENTED(0); 01901 01903 double evallog ( const vec &val ) const NOT_IMPLEMENTED(0); 01904 01905 vec mean() const { 01906 vec pom = zeros ( dim ); 01907 for ( int i = 0; i < n; i++ ) { 01908 pom += samples ( i ) * w ( i ); 01909 } 01910 return pom; 01911 } 01912 vec variance() const { 01913 vec pom = zeros ( dim ); 01914 for ( int i = 0; i < n; i++ ) { 01915 pom += pow ( samples ( i ), 2 ) * w ( i ); 01916 } 01917 return pom - pow ( mean(), 2 ); 01918 } 01920 void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const; 01921 01922 void to_setting ( Setting &set ) const; 01923 01934 void from_setting ( const Setting &set ); 01935 }; 01936 UIREGISTER(eEmp); 01937 01938 01940 01941 template<class sq_T> 01942 void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) { 01943 //Fixme test dimensions of mu0 and R0; 01944 mu = mu0; 01945 R = R0; 01946 validate(); 01947 }; 01948 01949 template<class sq_T> 01950 void enorm<sq_T>::from_setting ( const Setting &set ) { 01951 epdf::from_setting ( set ); //reads rv 01952 01953 UI::get ( mu, set, "mu", UI::compulsory ); 01954 mat Rtmp;// necessary for conversion 01955 UI::get ( Rtmp, set, "R", UI::compulsory ); 01956 R = Rtmp; // conversion 01957 } 01958 01959 template<class sq_T> 01960 void enorm<sq_T>::validate() { 01961 eEF::validate(); 01962 bdm_assert ( mu.length() == R.rows(), "mu and R parameters do not match" ); 01963 dim = mu.length(); 01964 } 01965 01966 template<class sq_T> 01967 void enorm<sq_T>::to_setting ( Setting &set ) const { 01968 epdf::to_setting ( set ); //reads rv 01969 UI::save ( mu, set, "mu"); 01970 UI::save ( R.to_mat(), set, "R"); 01971 } 01972 01973 01974 01975 template<class sq_T> 01976 void enorm<sq_T>::dupdate ( mat &v, double nu ) { 01977 // 01978 }; 01979 01980 // template<class sq_T> 01981 // void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) { 01982 // // 01983 // }; 01984 01985 template<class sq_T> 01986 vec enorm<sq_T>::sample() const { 01987 vec x ( dim ); 01988 #pragma omp critical 01989 NorRNG.sample_vector ( dim, x ); 01990 vec smp = R.sqrt_mult ( x ); 01991 01992 smp += mu; 01993 return smp; 01994 }; 01995 01996 // template<class sq_T> 01997 // double enorm<sq_T>::eval ( const vec &val ) const { 01998 // double pdfl,e; 01999 // pdfl = evallog ( val ); 02000 // e = exp ( pdfl ); 02001 // return e; 02002 // }; 02003 02004 template<class sq_T> 02005 double enorm<sq_T>::evallog_nn ( const vec &val ) const { 02006 // 1.83787706640935 = log(2pi) 02007 double tmp = -0.5 * ( R.invqform ( mu - val ) );// - lognc(); 02008 return tmp; 02009 }; 02010 02011 template<class sq_T> 02012 inline double enorm<sq_T>::lognc () const { 02013 // 1.83787706640935 = log(2pi) 02014 double tmp = 0.5 * ( R.cols() * 1.83787706640935 + R.logdet() ); 02015 return tmp; 02016 }; 02017 02018 02019 // template<class sq_T> 02020 // vec mlnorm<sq_T>::samplecond (const vec &cond, double &lik ) { 02021 // this->condition ( cond ); 02022 // vec smp = epdf.sample(); 02023 // lik = epdf.eval ( smp ); 02024 // return smp; 02025 // } 02026 02027 // template<class sq_T> 02028 // mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) { 02029 // int i; 02030 // int dim = rv.count(); 02031 // mat Smp ( dim,n ); 02032 // vec smp ( dim ); 02033 // this->condition ( cond ); 02034 // 02035 // for ( i=0; i<n; i++ ) { 02036 // smp = epdf.sample(); 02037 // lik ( i ) = epdf.eval ( smp ); 02038 // Smp.set_col ( i ,smp ); 02039 // } 02040 // 02041 // return Smp; 02042 // } 02043 02044 02045 template<class sq_T> 02046 shared_ptr<epdf> enorm<sq_T>::marginal ( const RV &rvn ) const { 02047 enorm<sq_T> *tmp = new enorm<sq_T> (); 02048 shared_ptr<epdf> narrow ( tmp ); 02049 marginal ( rvn, *tmp ); 02050 return narrow; 02051 } 02052 02053 template<class sq_T> 02054 void enorm<sq_T>::marginal ( const RV &rvn, enorm<sq_T> &target ) const { 02055 bdm_assert ( isnamed(), "rv description is not assigned" ); 02056 ivec irvn = rvn.dataind ( rv ); 02057 02058 sq_T Rn ( R, irvn ); // select rows and columns of R 02059 02060 target.set_rv ( rvn ); 02061 target.set_parameters ( mu ( irvn ), Rn ); 02062 } 02063 02064 template<class sq_T> 02065 shared_ptr<pdf> enorm<sq_T>::condition ( const RV &rvn ) const { 02066 mlnorm<sq_T> *tmp = new mlnorm<sq_T> (); 02067 shared_ptr<pdf> narrow ( tmp ); 02068 condition ( rvn, *tmp ); 02069 return narrow; 02070 } 02071 02072 template<class sq_T> 02073 void enorm<sq_T>::condition ( const RV &rvn, pdf &target ) const { 02074 typedef mlnorm<sq_T> TMlnorm; 02075 02076 bdm_assert ( isnamed(), "rvs are not assigned" ); 02077 TMlnorm &uptarget = dynamic_cast<TMlnorm &> ( target ); 02078 02079 RV rvc = rv.subt ( rvn ); 02080 bdm_assert ( ( rvc._dsize() + rvn._dsize() == rv._dsize() ), "wrong rvn" ); 02081 //Permutation vector of the new R 02082 ivec irvn = rvn.dataind ( rv ); 02083 ivec irvc = rvc.dataind ( rv ); 02084 ivec perm = concat ( irvn , irvc ); 02085 sq_T Rn ( R, perm ); 02086 02087 //fixme - could this be done in general for all sq_T? 02088 mat S = Rn.to_mat(); 02089 //fixme 02090 int n = rvn._dsize() - 1; 02091 int end = R.rows() - 1; 02092 mat S11 = S.get ( 0, n, 0, n ); 02093 mat S12 = S.get ( 0, n , rvn._dsize(), end ); 02094 mat S22 = S.get ( rvn._dsize(), end, rvn._dsize(), end ); 02095 02096 vec mu1 = mu ( irvn ); 02097 vec mu2 = mu ( irvc ); 02098 mat A = S12 * inv ( S22 ); 02099 sq_T R_n ( S11 - A *S12.T() ); 02100 02101 uptarget.set_rv ( rvn ); 02102 uptarget.set_rvc ( rvc ); 02103 uptarget.set_parameters ( A, mu1 - A*mu2, R_n ); 02104 uptarget.validate(); 02105 } 02106 02108 class dirac: public epdf { 02109 public: 02110 vec point; 02111 public: 02112 double evallog (const vec &dt) const { 02113 return -inf; 02114 } 02115 vec mean () const { 02116 return point; 02117 } 02118 vec variance () const { 02119 return zeros(point.length()); 02120 } 02121 void qbounds ( vec &lb, vec &ub, double percentage = 0.95 ) const { 02122 lb = point; 02123 ub = point; 02124 } 02126 const vec& _point() { 02127 return point; 02128 } 02129 void set_point(const vec& p) { 02130 point =p; 02131 dim=p.length(); 02132 } 02133 vec sample() const { 02134 return point; 02135 } 02136 }; 02137 02138 02140 02141 template<class sq_T> 02142 void mgnorm<sq_T >::set_parameters ( const shared_ptr<fnc> &g0, const sq_T &R0 ) { 02143 g = g0; 02144 this->iepdf.set_parameters ( zeros ( g->dimension() ), R0 ); 02145 } 02146 02147 template<class sq_T> 02148 void mgnorm<sq_T >::condition ( const vec &cond ) { 02149 this->iepdf._mu() = g->eval ( cond ); 02150 }; 02151 02153 template<class sq_T> 02154 std::ostream &operator<< ( std::ostream &os, mlnorm<sq_T> &ml ) { 02155 os << "A:" << ml.A << endl; 02156 os << "mu:" << ml.mu_const << endl; 02157 os << "R:" << ml._R() << endl; 02158 return os; 02159 }; 02160 02161 } 02162 #endif //EF_H
Generated on 2 Dec 2013 for mixpp by
1.4.7
