- Main Page
- Related Pages
- Modules
- Classes
- Files

# libEF.h

Go to the documentation of this file.

    00001
    00013 #ifndef EF_H
    00014 #define EF_H
    00015
    00016
    00017 #include "libBM.h"
    00018 #include "../math/chmat.h"
    00019 //#include <std>
    00020
    00021 namespace bdm {
    00022
    00023
    00025 extern Uniform_RNG UniRNG;
    00027 extern Normal_RNG NorRNG;
    00029 extern Gamma_RNG GamRNG;
    00030
    00037 class eEF : public epdf {
    00038 public:
    00039 //      eEF() :epdf() {};
    00041         eEF ( ) :epdf ( ) {};
    00043         virtual double lognc() const =0;
    00045         virtual void dupdate ( mat &v ) {it_error ( "Not implemented" );};
    00047         virtual double evallog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
    00049         virtual double evallog ( const vec &val ) const {double tmp;tmp= evallog_nn ( val )-lognc();it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); return tmp;}
    00051         virtual vec evallog ( const mat &Val ) const {
    00052                 vec x ( Val.cols() );
    00053                 for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog_nn ( Val.get_col ( i ) ) ;}
    00054                 return x-lognc();
    00055         }
    00057         virtual void pow ( double p ) {it_error ( "Not implemented" );};
    00058 };
    00059
    00066 class mEF : public mpdf {
    00067
    00068 public:
    00070         mEF ( ) :mpdf ( ) {};
    00071 };
    00072
    00074 class BMEF : public BM {
    00075 protected:
    00077         double frg;
    00079         double last_lognc;
    00080 public:
    00082         BMEF ( double frg0=1.0 ) :BM ( ), frg ( frg0 ) {}
    00084         BMEF ( const BMEF &B ) :BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
    00086         virtual void set_statistics ( const BMEF* BM0 ) {it_error ( "Not implemented" );};
    00088         virtual void bayes ( const vec &data, const double w ) {};
    00089         //original Bayes
    00090         void bayes ( const vec &dt );
    00092         virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
    00094 //      virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
    00095
    00096         BMEF* _copy_ ( bool changerv=false ) const {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
    00097 };
    00098
    00099 template<class sq_T>
    00100 class mlnorm;
    00101
    00107 template<class sq_T>
    00108 class enorm : public eEF {
    00109 protected:
    00111         vec mu;
    00113         sq_T R;
    00114 public:
    00117
    00118         enorm ( ) :eEF ( ), mu ( ),R ( ) {};
    00119         enorm ( const vec &mu,const sq_T &R ) {set_parameters ( mu,R );}
    00120         void set_parameters ( const vec &mu,const sq_T &R );
    00122
    00125
    00127         void dupdate ( mat &v,double nu=1.0 );
    00128
    00129         vec sample() const;
    00130         mat sample ( int N ) const;
    00131         double evallog_nn ( const vec &val ) const;
    00132         double lognc () const;
    00133         vec mean() const {return mu;}
    00134         vec variance() const {return diag ( R.to_mat() );}
    00135 //      mlnorm<sq_T>* condition ( const RV &rvn ) const ;
    00136         mpdf* condition ( const RV &rvn ) const ;
    00137 //      enorm<sq_T>* marginal ( const RV &rv ) const;
    00138         epdf* marginal ( const RV &rv ) const;
    00140
    00143
    00144         vec& _mu() {return mu;}
    00145         void set_mu ( const vec mu0 ) { mu=mu0;}
    00146         sq_T& _R() {return R;}
    00147         const sq_T& _R() const {return R;}
    00149
    00150 };
    00151
    00158 class egiw : public eEF {
    00159 protected:
    00161         ldmat V;
    00163         double nu;
    00165         int dimx;
    00167         int nPsi;
    00168 public:
    00171         egiw() :eEF() {};
    00172         egiw ( int dimx0, ldmat V0, double nu0=-1.0 ) :eEF() {set_parameters ( dimx0,V0, nu0 );};
    00173
    00174         void set_parameters ( int dimx0, ldmat V0, double nu0=-1.0 ) {
    00175                 dimx=dimx0;
    00176                 nPsi = V0.rows()-dimx;
    00177                 dim = dimx* ( dimx+nPsi ); // size(R) + size(Theta)
    00178
    00179                 V=V0;
    00180                 if ( nu0<0 ) {
    00181                         nu = 0.1 +nPsi +2*dimx +2; // +2 assures finite expected value of R
    00182                         // terms before that are sufficient for finite normalization
    00183                 }
    00184                 else {
    00185                         nu=nu0;
    00186                 }
    00187         }
    00189
    00190         vec sample() const;
    00191         vec mean() const;
    00192         vec variance() const;
    00193         void mean_mat ( mat &M, mat&R ) const;
    00195         double evallog_nn ( const vec &val ) const;
    00196         double lognc () const;
    00197         void pow ( double p ) {V*=p;nu*=p;};
    00198
    00201
    00202         ldmat& _V() {return V;}
    00203         const ldmat& _V() const {return V;}
    00204         double& _nu()  {return nu;}
    00205         const double& _nu() const {return nu;}
    00207 };
    00208
    00217 class eDirich: public eEF {
    00218 protected:
    00220         vec beta;
    00222         double gamma;
    00223 public:
    00226
    00227         eDirich () : eEF ( ) {};
    00228         eDirich ( const eDirich &D0 ) : eEF () {set_parameters ( D0.beta );};
    00229         eDirich ( const vec &beta0 ) {set_parameters ( beta0 );};
    00230         void set_parameters ( const vec &beta0 ) {
    00231                 beta= beta0;
    00232                 dim = beta.length();
    00233                 gamma = sum ( beta );
    00234         }
    00236
    00237         vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
    00238         vec mean() const {return beta/gamma;};
    00239         vec variance() const {return elem_mult ( beta, ( beta+1 ) ) / ( gamma* ( gamma+1 ) );}
    00241         double evallog_nn ( const vec &val ) const {
    00242                 double tmp; tmp= ( beta-1 ) *log ( val );               it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
    00243                 return tmp;
    00244         };
    00245         double lognc () const {
    00246                 double tmp;
    00247                 double gam=sum ( beta );
    00248                 double lgb=0.0;
    00249                 for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
    00250                 tmp= lgb-lgamma ( gam );
    00251                 it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
    00252                 return tmp;
    00253         };
    00255         vec& _beta()  {return beta;}
    00257 };
    00258
    00260 class multiBM : public BMEF {
    00261 protected:
    00263         eDirich est;
    00265         vec &beta;
    00266 public:
    00268         multiBM ( ) : BMEF ( ),est ( ),beta ( est._beta() ) {if ( beta.length() >0 ) {last_lognc=est.lognc();}else{last_lognc=0.0;}}
    00270         multiBM ( const multiBM &B ) : BMEF ( B ),est ( B.est ),beta ( est._beta() ) {}
    00272         void set_statistics ( const BM* mB0 ) {const multiBM* mB=dynamic_cast<const multiBM*> ( mB0 ); beta=mB->beta;}
    00273         void bayes ( const vec &dt ) {
    00274                 if ( frg<1.0 ) {beta*=frg;last_lognc=est.lognc();}
    00275                 beta+=dt;
    00276                 if ( evalll ) {ll=est.lognc()-last_lognc;}
    00277         }
    00278         double logpred ( const vec &dt ) const {
    00279                 eDirich pred ( est );
    00280                 vec &beta = pred._beta();
    00281
    00282                 double lll;
    00283                 if ( frg<1.0 )
    00284                         {beta*=frg;lll=pred.lognc();}
    00285                 else
    00286                         if ( evalll ) {lll=last_lognc;}
    00287                         else{lll=pred.lognc();}
    00288
    00289                 beta+=dt;
    00290                 return pred.lognc()-lll;
    00291         }
    00292         void flatten ( const BMEF* B ) {
    00293                 const multiBM* E=dynamic_cast<const multiBM*> ( B );
    00294                 // sum(beta) should be equal to sum(B.beta)
    00295                 const vec &Eb=E->beta;//const_cast<multiBM*> ( E )->_beta();
    00296                 beta*= ( sum ( Eb ) /sum ( beta ) );
    00297                 if ( evalll ) {last_lognc=est.lognc();}
    00298         }
    00299         const epdf& posterior() const {return est;};
    00300         const eDirich* _e() const {return &est;};
    00301         void set_parameters ( const vec &beta0 ) {
    00302                 est.set_parameters ( beta0 );
    00303                 if ( evalll ) {last_lognc=est.lognc();}
    00304         }
    00305 };
    00306
    00316 class egamma : public eEF {
    00317 protected:
    00319         vec alpha;
    00321         vec beta;
    00322 public :
    00325         egamma ( ) :eEF ( ), alpha(0), beta(0) {};
    00326         egamma ( const vec &a, const vec &b ) {set_parameters ( a, b );};
    00327         void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;dim = alpha.length();};
    00329
    00330         vec sample() const;
    00332 //      mat sample ( int N ) const;
    00333         double evallog ( const vec &val ) const;
    00334         double lognc () const;
    00336         vec& _alpha() {return alpha;}
    00337         vec& _beta() {return beta;}
    00338         vec mean() const {return elem_div ( alpha,beta );}
    00339         vec variance() const {return elem_div ( alpha,elem_mult ( beta,beta ) ); }
    00340 };
    00341
    00358 class eigamma : public egamma {
    00359 protected:
    00360 public :
    00365
    00366         vec sample() const {return 1.0/egamma::sample();};
    00368         vec mean() const {return elem_div ( beta,alpha-1 );}
    00369         vec variance() const {vec mea=mean(); return elem_div ( elem_mult ( mea,mea ),alpha-2 );}
    00370 };
    00371 /*
    00373 class emix : public epdf {
    00374 protected:
    00375         int n;
    00376         vec &w;
    00377         Array<epdf*> Coms;
    00378 public:
    00380         emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
    00381         void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
    00382         vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
    00383         vec sample() {it_error ( "Not implemented" );return 0;}
    00384 };
    00385 */
    00386
    00388
    00389 class euni: public epdf {
    00390 protected:
    00392         vec low;
    00394         vec high;
    00396         vec distance;
    00398         double nk;
    00400         double lnk;
    00401 public:
    00404         euni ( ) :epdf ( ) {}
    00405         euni ( const vec &low0, const vec &high0 ) {set_parameters ( low0,high0 );}
    00406         void set_parameters ( const vec &low0, const vec &high0 ) {
    00407                 distance = high0-low0;
    00408                 it_assert_debug ( min ( distance ) >0.0,"bad support" );
    00409                 low = low0;
    00410                 high = high0;
    00411                 nk = prod ( 1.0/distance );
    00412                 lnk = log ( nk );
    00413                 dim = low.length();
    00414         }
    00416
    00417         double eval ( const vec &val ) const  {return nk;}
    00418         double evallog ( const vec &val ) const  {return lnk;}
    00419         vec sample() const {
    00420                 vec smp ( dim );
    00421 #pragma omp critical
    00422                 UniRNG.sample_vector ( dim ,smp );
    00423                 return low+elem_mult ( distance,smp );
    00424         }
    00426         vec mean() const {return ( high-low ) /2.0;}
    00427         vec variance() const {return ( pow ( high,2 ) +pow ( low,2 ) +elem_mult ( high,low ) ) /3.0;}
    00428 };
    00429
    00430
    00436 template<class sq_T>
    00437 class mlnorm : public mEF {
    00438 protected:
    00440         enorm<sq_T> epdf;
    00441         mat A;
    00442         vec mu_const;
    00443         vec& _mu; //cached epdf.mu;
    00444 public:
    00447         mlnorm ( ) :mEF (),epdf ( ),A ( ),_mu ( epdf._mu() ) {ep =&epdf; };
    00448         mlnorm ( const  mat &A, const vec &mu0, const sq_T &R ) :epdf ( ),_mu ( epdf._mu() ) {
    00449                 ep =&epdf; set_parameters ( A,mu0,R );
    00450         };
    00452         void set_parameters ( const  mat &A, const vec &mu0, const sq_T &R );
    00455         void condition ( const vec &cond );
    00456
    00458         vec& _mu_const() {return mu_const;}
    00460         mat& _A() {return A;}
    00462         mat _R() {return epdf._R().to_mat();}
    00463
    00464         template<class sq_M>
    00465         friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M> &ml );
    00466 };
    00467
    00469 template<class sq_T>
    00470 class mgnorm : public mEF {
    00471 protected:
    00473         enorm<sq_T> epdf;
    00474         vec &mu;
    00475         fnc* g;
    00476 public:
    00478         mgnorm() :mu ( epdf._mu() ) {ep=&epdf;}
    00480         void set_parameters ( fnc* g0, const sq_T &R0 ) {g=g0; epdf.set_parameters ( zeros ( g->dimension() ), R0 );}
    00481         void condition ( const vec &cond ) {mu=g->eval ( cond );};
    00482 };
    00483
    00491 class mlstudent : public mlnorm<ldmat> {
    00492 protected:
    00493         ldmat Lambda;
    00494         ldmat &_R;
    00495         ldmat Re;
    00496 public:
    00497         mlstudent ( ) :mlnorm<ldmat> (),
    00498                         Lambda (),      _R ( epdf._R() ) {}
    00499         void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 ) {
    00500                 it_assert_debug ( A0.rows() ==mu0.length(),"" );
    00501                 it_assert_debug ( R0.rows() ==A0.rows(),"" );
    00502
    00503                 epdf.set_parameters ( mu0,Lambda ); //
    00504                 A = A0;
    00505                 mu_const = mu0;
    00506                 Re=R0;
    00507                 Lambda = Lambda0;
    00508         }
    00509         void condition ( const vec &cond ) {
    00510                 _mu = A*cond + mu_const;
    00511                 double zeta;
    00512                 //ugly hack!
    00513                 if ( ( cond.length() +1 ) ==Lambda.rows() ) {
    00514                         zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) );
    00515                 }
    00516                 else {
    00517                         zeta = Lambda.invqform ( cond );
    00518                 }
    00519                 _R = Re;
    00520                 _R*= ( 1+zeta );// / ( nu ); << nu is in Re!!!!!!
    00521         };
    00522
    00523 };
    00533 class mgamma : public mEF {
    00534 protected:
    00536         egamma epdf;
    00538         double k;
    00540         vec &_beta;
    00541
    00542 public:
    00544         mgamma ( ) : mEF ( ), epdf (), _beta ( epdf._beta() ) {ep=&epdf;};
    00546         void set_parameters ( double k, const vec &beta0 );
    00547         void condition ( const vec &val ) {_beta=k/val;};
    00548 };
    00549
    00559 class migamma : public mEF {
    00560 protected:
    00562         eigamma epdf;
    00564         double k;
    00566         vec &_alpha;
    00568         vec &_beta;
    00569
    00570 public:
    00573         migamma ( ) : mEF (), epdf ( ), _alpha ( epdf._alpha() ), _beta ( epdf._beta() ) {ep=&epdf;};
    00574         migamma ( const migamma &m ) : mEF (), epdf ( m.epdf ), _alpha ( epdf._alpha() ), _beta ( epdf._beta() ) {ep=&epdf;};
    00576
    00578         void set_parameters ( int len, double k0 ) {
    00579                 k=k0;
    00580                 epdf.set_parameters ( ( 1.0/ ( k*k ) +2.0 ) *ones ( len ) /*alpha*/, ones ( len ) /*beta*/ );
    00581                 dimc = dimension();
    00582         };
    00583         void condition ( const vec &val ) {
    00584                 _beta=elem_mult ( val, ( _alpha-1.0 ) );
    00585         };
    00586 };
    00587
    00599 class mgamma_fix : public mgamma {
    00600 protected:
    00602         double l;
    00604         vec refl;
    00605 public:
    00607         mgamma_fix ( ) : mgamma ( ),refl () {};
    00609         void set_parameters ( double k0 , vec ref0, double l0 ) {
    00610                 mgamma::set_parameters ( k0, ref0 );
    00611                 refl=pow ( ref0,1.0-l0 );l=l0;
    00612                 dimc=dimension();
    00613         };
    00614
    00615         void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); _beta=k/mean;};
    00616 };
    00617
    00618
    00631 class migamma_ref : public migamma {
    00632 protected:
    00634         double l;
    00636         vec refl;
    00637 public:
    00639         migamma_ref ( ) : migamma (),refl ( ) {};
    00641         void set_parameters ( double k0 , vec ref0, double l0 ) {
    00642                 migamma::set_parameters ( ref0.length(), k0 );
    00643                 refl=pow ( ref0,1.0-l0 );
    00644                 l=l0;
    00645                 dimc = dimension();
    00646         };
    00647
    00648         void condition ( const vec &val ) {
    00649                 vec mean=elem_mult ( refl,pow ( val,l ) );
    00650                 migamma::condition ( mean );
    00651         };
    00652 };
    00653
    00663 class elognorm: public enorm<ldmat>{
    00664         public:
    00665                 vec sample() const {return exp(enorm<ldmat>::sample());};
    00666                 vec mean() const {vec var=enorm<ldmat>::variance();return exp(mu - 0.5*var);};
    00667
    00668 };
    00669
    00681 class mlognorm : public mpdf {
    00682         protected:
    00683                 elognorm eno;
    00685                 double sig2;
    00687                 vec &mu;
    00688         public:
    00690                 mlognorm ( ) : eno (), mu(eno._mu()) {ep=&eno;};
    00692                 void set_parameters ( int size, double k) {
    00693                         sig2 = 0.5*log(k*k+1);
    00694                         eno.set_parameters(zeros(size),2*sig2*eye(size));
    00695
    00696                         dimc = size;
    00697                 };
    00698
    00699                 void condition ( const vec &val ) {
    00700                         mu=log(val)-sig2;//elem_mult ( refl,pow ( val,l ) );
    00701                 };
    00702 };
    00703
    00707 class iW : public epdf{
    00708         public:
    00709                 vec sample(){return vec();}
    00710 };
    00711
    00713 enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
    00719 class eEmp: public epdf {
    00720 protected :
    00722         int n;
    00724         vec w;
    00726         Array<vec> samples;
    00727 public:
    00730         eEmp ( ) :epdf ( ),w ( ),samples ( ) {};
    00731         eEmp ( const eEmp &e ) : epdf ( e ), w ( e.w ), samples ( e.samples ) {};
    00733
    00735         void set_statistics ( const vec &w0, const epdf* pdf0 );
    00737         void set_statistics ( const epdf* pdf0 , int n ) {set_statistics ( ones ( n ) /n,pdf0 );};
    00739         void set_samples ( const epdf* pdf0 );
    00741         void set_parameters ( int n0, bool copy=true ) {n=n0; w.set_size ( n0,copy );samples.set_size ( n0,copy );};
    00743         vec& _w()  {return w;};
    00745         const vec& _w() const {return w;};
    00747         Array<vec>& _samples() {return samples;};
    00749         const Array<vec>& _samples() const {return samples;};
    00751         ivec resample ( RESAMPLING_METHOD method=SYSTEMATIC );
    00753         vec sample() const {it_error ( "Not implemented" );return 0;}
    00755         double evallog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
    00756         vec mean() const {
    00757                 vec pom=zeros ( dim );
    00758                 for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
    00759                 return pom;
    00760         }
    00761         vec variance() const {
    00762                 vec pom=zeros ( dim );
    00763                 for ( int i=0;i<n;i++ ) {pom+=pow ( samples ( i ),2 ) *w ( i );}
    00764                 return pom-pow ( mean(),2 );
    00765         }
    00767         void qbounds ( vec &lb, vec &ub, double perc=0.95 ) const {
    00768                 // lb in inf so than it will be pushed below;
    00769                 lb.set_size(dim);
    00770                 ub.set_size(dim);
    00771                 lb = std::numeric_limits<double>::infinity();
    00772                 ub = -std::numeric_limits<double>::infinity();
    00773                 int j;
    00774                 for ( int i=0;i<n;i++ ) {
    00775                         for ( j=0;j<dim; j++ ) {
    00776                                 if ( samples ( i ) ( j ) <lb ( j ) ) {lb ( j ) =samples ( i ) ( j );}
    00777                                 if ( samples ( i ) ( j ) >ub ( j ) ) {ub ( j ) =samples ( i ) ( j );}
    00778                         }
    00779                 }
    00780         }
    00781 };
    00782
    00783
    00785
    00786 template<class sq_T>
    00787 void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) {
    00788 //Fixme test dimensions of mu0 and R0;
    00789         mu = mu0;
    00790         R = R0;
    00791         dim = mu0.length();
    00792 };
    00793
    00794 template<class sq_T>
    00795 void enorm<sq_T>::dupdate ( mat &v, double nu ) {
    00796         //
    00797 };
    00798
    00799 // template<class sq_T>
    00800 // void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
    00801 //      //
    00802 // };
    00803
    00804 template<class sq_T>
    00805 vec enorm<sq_T>::sample() const {
    00806         vec x ( dim );
    00807 #pragma omp critical
    00808         NorRNG.sample_vector ( dim,x );
    00809         vec smp = R.sqrt_mult ( x );
    00810
    00811         smp += mu;
    00812         return smp;
    00813 };
    00814
    00815 template<class sq_T>
    00816 mat enorm<sq_T>::sample ( int N ) const {
    00817         mat X ( dim,N );
    00818         vec x ( dim );
    00819         vec pom;
    00820         int i;
    00821
    00822         for ( i=0;i<N;i++ ) {
    00823 #pragma omp critical
    00824                 NorRNG.sample_vector ( dim,x );
    00825                 pom = R.sqrt_mult ( x );
    00826                 pom +=mu;
    00827                 X.set_col ( i, pom );
    00828         }
    00829
    00830         return X;
    00831 };
    00832
    00833 // template<class sq_T>
    00834 // double enorm<sq_T>::eval ( const vec &val ) const {
    00835 //      double pdfl,e;
    00836 //      pdfl = evallog ( val );
    00837 //      e = exp ( pdfl );
    00838 //      return e;
    00839 // };
    00840
    00841 template<class sq_T>
    00842 double enorm<sq_T>::evallog_nn ( const vec &val ) const {
    00843         // 1.83787706640935 = log(2pi)
    00844         double tmp=-0.5* ( R.invqform ( mu-val ) );// - lognc();
    00845         return  tmp;
    00846 };
    00847
    00848 template<class sq_T>
    00849 inline double enorm<sq_T>::lognc () const {
    00850         // 1.83787706640935 = log(2pi)
    00851         double tmp=0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
    00852         return tmp;
    00853 };
    00854
    00855 template<class sq_T>
    00856 void mlnorm<sq_T>::set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 ) {
    00857         it_assert_debug ( A0.rows() ==mu0.length(),"" );
    00858         it_assert_debug ( A0.rows() ==R0.rows(),"" );
    00859
    00860         epdf.set_parameters ( zeros ( A0.rows() ),R0 );
    00861         A = A0;
    00862         mu_const = mu0;
    00863         dimc=A0.cols();
    00864 }
    00865
    00866 // template<class sq_T>
    00867 // vec mlnorm<sq_T>::samplecond (const  vec &cond, double &lik ) {
    00868 //      this->condition ( cond );
    00869 //      vec smp = epdf.sample();
    00870 //      lik = epdf.eval ( smp );
    00871 //      return smp;
    00872 // }
    00873
    00874 // template<class sq_T>
    00875 // mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) {
    00876 //      int i;
    00877 //      int dim = rv.count();
    00878 //      mat Smp ( dim,n );
    00879 //      vec smp ( dim );
    00880 //      this->condition ( cond );
    00881 //
    00882 //      for ( i=0; i<n; i++ ) {
    00883 //              smp = epdf.sample();
    00884 //              lik ( i ) = epdf.eval ( smp );
    00885 //              Smp.set_col ( i ,smp );
    00886 //      }
    00887 //
    00888 //      return Smp;
    00889 // }
    00890
    00891 template<class sq_T>
    00892 void mlnorm<sq_T>::condition ( const vec &cond ) {
    00893         _mu = A*cond + mu_const;
    00894 //R is already assigned;
    00895 }
    00896
    00897 template<class sq_T>
    00898 epdf* enorm<sq_T>::marginal ( const RV &rvn ) const {
    00899         it_assert_debug ( isnamed(), "rv description is not assigned" );
    00900         ivec irvn = rvn.dataind ( rv );
    00901
    00902         sq_T Rn ( R,irvn ); //select rows and columns of R
    00903
    00904         enorm<sq_T>* tmp = new enorm<sq_T>;
    00905         tmp->set_rv ( rvn );
    00906         tmp->set_parameters ( mu ( irvn ), Rn );
    00907         return tmp;
    00908 }
    00909
    00910 template<class sq_T>
    00911 mpdf* enorm<sq_T>::condition ( const RV &rvn ) const {
    00912
    00913         it_assert_debug ( isnamed(),"rvs are not assigned" );
    00914
    00915         RV rvc = rv.subt ( rvn );
    00916         it_assert_debug ( ( rvc._dsize() +rvn._dsize() ==rv._dsize() ),"wrong rvn" );
    00917         //Permutation vector of the new R
    00918         ivec irvn = rvn.dataind ( rv );
    00919         ivec irvc = rvc.dataind ( rv );
    00920         ivec perm=concat ( irvn , irvc );
    00921         sq_T Rn ( R,perm );
    00922
    00923         //fixme - could this be done in general for all sq_T?
    00924         mat S=Rn.to_mat();
    00925         //fixme
    00926         int n=rvn._dsize()-1;
    00927         int end=R.rows()-1;
    00928         mat S11 = S.get ( 0,n, 0, n );
    00929         mat S12 = S.get ( 0, n , rvn._dsize(), end );
    00930         mat S22 = S.get ( rvn._dsize(), end, rvn._dsize(), end );
    00931
    00932         vec mu1 = mu ( irvn );
    00933         vec mu2 = mu ( irvc );
    00934         mat A=S12*inv ( S22 );
    00935         sq_T R_n ( S11 - A *S12.T() );
    00936
    00937         mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( );
    00938         tmp->set_rv ( rvn ); tmp->set_rvc ( rvc );
    00939         tmp->set_parameters ( A,mu1-A*mu2,R_n );
    00940         return tmp;
    00941 }
    00942
    00944
    00945 template<class sq_T>
    00946 std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_T> &ml ) {
    00947         os << "A:"<< ml.A<<endl;
    00948         os << "mu:"<< ml.mu_const<<endl;
    00949         os << "R:" << ml.epdf._R().to_mat() <<endl;
    00950         return os;
    00951 };
    00952
    00953 }
    00954 #endif //EF_H