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

# work/git/mixpp/bdm/stat/libEF.h

Go to the documentation of this file.

    00001
    00013 #ifndef EF_H
    00014 #define EF_H
    00015
    00016 #include <itpp/itbase.h>
    00017 #include "../math/libDC.h"
    00018 #include "libBM.h"
    00019 #include "../itpp_ext.h"
    00020 //#include <std>
    00021
    00022 using namespace itpp;
    00023
    00024
    00026 extern Uniform_RNG UniRNG;
    00028 extern Normal_RNG NorRNG;
    00030 extern Gamma_RNG GamRNG;
    00031
    00038 class eEF : public epdf {
    00039 public:
    00040 //      eEF() :epdf() {};
    00042         eEF ( const RV &rv ) :epdf ( rv ) {};
    00044         virtual double lognc() const =0;
    00046         virtual void dupdate ( mat &v ) {it_error ( "Not implemented" );};
    00048         virtual double evallog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
    00050         virtual double evallog ( const vec &val ) const {double tmp;tmp= evallog_nn ( val )-lognc();it_assert_debug(std::isfinite(tmp),"Infinite value"); return tmp;}
    00052         virtual vec evallog ( const mat &Val ) const {
    00053                 vec x ( Val.cols() );
    00054                 for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog_nn ( Val.get_col ( i ) ) ;}
    00055                 return x-lognc();
    00056         }
    00058         virtual void pow ( double p ) {it_error ( "Not implemented" );};
    00059 };
    00060
    00067 class mEF : public mpdf {
    00068
    00069 public:
    00071         mEF ( const RV &rv0, const RV &rvc0 ) :mpdf ( rv0,rvc0 ) {};
    00072 };
    00073
    00075 class BMEF : public BM {
    00076 protected:
    00078         double frg;
    00080         double last_lognc;
    00081 public:
    00083         BMEF ( const RV &rv, double frg0=1.0 ) :BM ( rv ), frg ( frg0 ) {}
    00085         BMEF ( const BMEF &B ) :BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
    00087         virtual void set_statistics ( const BMEF* BM0 ) {it_error ( "Not implemented" );};
    00089         virtual void bayes ( const vec &data, const double w ) {};
    00090         //original Bayes
    00091         void bayes ( const vec &dt );
    00093         virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
    00095 //      virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
    00096
    00097         BMEF* _copy_ ( bool changerv=false ) {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
    00098 };
    00099
    00100 template<class sq_T>
    00101 class mlnorm;
    00102
    00108 template<class sq_T>
    00109 class enorm : public eEF {
    00110 protected:
    00112         vec mu;
    00114         sq_T R;
    00116         int dim;
    00117 public:
    00119         enorm ( const RV &rv );
    00121         void set_parameters ( const vec &mu,const sq_T &R );
    00123         //void tupdate ( double phi, mat &vbar, double nubar );
    00125         void dupdate ( mat &v,double nu=1.0 );
    00126
    00127         vec sample() const;
    00129         mat sample ( int N ) const;
    00130 //      double eval ( const vec &val ) 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;
    00139 //Access methods
    00141         vec& _mu() {return mu;}
    00142
    00144         void set_mu ( const vec mu0 ) { mu=mu0;}
    00145
    00147         sq_T& _R() {return R;}
    00148         const sq_T& _R() const {return R;}
    00149
    00151 //      mat getR () {return R.to_mat();}
    00152 };
    00153
    00160 class egiw : public eEF {
    00161 protected:
    00163         ldmat V;
    00165         double nu;
    00167         int xdim;
    00169         int nPsi;
    00170 public:
    00172         egiw ( RV rv, mat V0, double nu0=-1.0 ) : eEF ( rv ), V ( V0 ), nu ( nu0 ) {
    00173                 xdim = rv.count() /V.rows();
    00174                 it_assert_debug ( rv.count() ==xdim*V.rows(),"Incompatible V0." );
    00175                 nPsi = V.rows()-xdim;
    00176                 //set mu to have proper normalization and
    00177                 if (nu0<0){
    00178                         nu = 0.1 +nPsi +2*xdim +2; // +2 assures finite expected value of R
    00179                         // terms before that are sufficient for finite normalization
    00180                 }
    00181         }
    00183         egiw ( RV rv, ldmat V0, double nu0=-1.0 ) : eEF ( rv ), V ( V0 ), nu ( nu0 ) {
    00184                 xdim = rv.count() /V.rows();
    00185                 it_assert_debug ( rv.count() ==xdim*V.rows(),"Incompatible V0." );
    00186                 nPsi = V.rows()-xdim;
    00187                 if (nu0<0){
    00188                         nu = 0.1 +nPsi +2*xdim +2; // +2 assures finite expected value of R
    00189                         // terms before that are sufficient for finite normalization
    00190                 }
    00191         }
    00192
    00193         vec sample() const;
    00194         vec mean() const;
    00195         vec variance() const{it_error("Not implemented"); return vec(0);};
    00196         void mean_mat ( mat &M, mat&R ) const;
    00198         double evallog_nn ( const vec &val ) const;
    00199         double lognc () const;
    00200
    00201         //Access
    00203         ldmat& _V() {return V;}
    00205         const ldmat& _V() const {return V;}
    00207         double& _nu()  {return nu;}
    00208         const double& _nu() const {return nu;}
    00209         void pow ( double p ) {V*=p;nu*=p;};
    00210 };
    00211
    00220 class eDirich: public eEF {
    00221 protected:
    00223         vec beta;
    00225         double gamma;
    00226 public:
    00228         eDirich ( const RV &rv, const vec &beta0 ) : eEF ( rv ),beta ( beta0 ) {it_assert_debug ( rv.count() ==beta.length(),"Incompatible statistics" ); };
    00230         eDirich ( const eDirich &D0 ) : eEF ( D0.rv ),beta ( D0.beta ) {};
    00231         vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
    00232         vec mean() const {return beta/gamma;};
    00233         vec variance() const {return elem_mult(beta,(beta+1))/ (gamma*(gamma+1));}
    00235         double evallog_nn ( const vec &val ) const {double tmp; tmp=( beta-1 ) *log ( val );            it_assert_debug(std::isfinite(tmp),"Infinite value");
    00236         return tmp;};
    00237         double lognc () const {
    00238                 double tmp;
    00239                 double gam=sum ( beta );
    00240                 double lgb=0.0;
    00241                 for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
    00242                 tmp= lgb-lgamma ( gam );
    00243                 it_assert_debug(std::isfinite(tmp),"Infinite value");
    00244                 return tmp;
    00245         };
    00247         vec& _beta()  {return beta;}
    00249         void set_parameters ( const vec &beta0 ) {
    00250                 if ( beta0.length() !=beta.length() ) {
    00251                         it_assert_debug ( rv.length() ==1,"Undefined" );
    00252                         rv.set_size ( 0,beta0.length() );
    00253                 }
    00254                 beta= beta0;
    00255                 gamma = sum(beta);
    00256         }
    00257 };
    00258
    00260 class multiBM : public BMEF {
    00261 protected:
    00263         eDirich est;
    00265         vec &beta;
    00266 public:
    00268         multiBM ( const RV &rv, const vec beta0 ) : BMEF ( rv ),est ( rv,beta0 ),beta ( est._beta() ) {if(beta.length()>0){last_lognc=est.lognc();}else{last_lognc=0.0;}}
    00270         multiBM ( const multiBM &B ) : BMEF ( B ),est ( rv,B.beta ),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& _epdf() const {return est;};
    00300         const eDirich* _e() const {return &est;};
    00301         void set_parameters ( const vec &beta0 ) {
    00302                 est.set_parameters ( beta0 );
    00303                 rv = est._rv();
    00304                 if ( evalll ) {last_lognc=est.lognc();}
    00305         }
    00306 };
    00307
    00317 class egamma : public eEF {
    00318 protected:
    00320         vec alpha;
    00322         vec beta;
    00323 public :
    00325         egamma ( const RV &rv ) :eEF ( rv ), alpha(rv.count()), beta(rv.count()) {};
    00327         void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;};
    00328         vec sample() const;
    00330 //      mat sample ( int N ) const;
    00331         double evallog ( const vec &val ) const;
    00332         double lognc () const;
    00334         void _param ( vec* &a, vec* &b ) {a=&alpha;b=&beta;};
    00335         vec mean() const {return elem_div(alpha,beta);}
    00336         vec variance() const {return elem_div(alpha,elem_mult(beta,beta)); }
    00337 };
    00338
    00353 class eigamma : public eEF {
    00354         protected:
    00356                 vec* alpha;
    00358                 vec* beta;
    00360                 egamma eg;
    00361         public :
    00363                 eigamma ( const RV &rv ) :eEF ( rv ), eg(rv) {eg._param(alpha,beta);};
    00365                 void set_parameters ( const vec &a, const vec &b ) {*alpha=a,*beta=b;};
    00366                 vec sample() const {return 1.0/eg.sample();};
    00368 //      mat sample ( int N ) const;
    00369                 double evallog ( const vec &val ) const {return eg.evallog(val);};
    00370                 double lognc () const {return eg.lognc();};
    00372                 void _param ( vec* &a, vec* &b ) {a=alpha;b=beta;};
    00373                 vec mean() const {return elem_div(*beta,*alpha-1);}
    00374                 vec variance() const {vec mea=mean(); return elem_div(elem_mult(mea,mea),*alpha-2);}
    00375 };
    00376 /*
    00378 class emix : public epdf {
    00379 protected:
    00380         int n;
    00381         vec &w;
    00382         Array<epdf*> Coms;
    00383 public:
    00385         emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
    00386         void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
    00387         vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
    00388         vec sample() {it_error ( "Not implemented" );return 0;}
    00389 };
    00390 */
    00391
    00393
    00394 class euni: public epdf {
    00395 protected:
    00397         vec low;
    00399         vec high;
    00401         vec distance;
    00403         double nk;
    00405         double lnk;
    00406 public:
    00408         euni ( const RV rv ) :epdf ( rv ) {}
    00409         double eval ( const vec &val ) const  {return nk;}
    00410         double evallog ( const vec &val ) const  {return lnk;}
    00411         vec sample() const {
    00412                 vec smp ( rv.count() );
    00413 #pragma omp critical
    00414                 UniRNG.sample_vector ( rv.count(),smp );
    00415                 return low+elem_mult ( distance,smp );
    00416         }
    00418         void set_parameters ( const vec &low0, const vec &high0 ) {
    00419                 distance = high0-low0;
    00420                 it_assert_debug ( min ( distance ) >0.0,"bad support" );
    00421                 low = low0;
    00422                 high = high0;
    00423                 nk = prod ( 1.0/distance );
    00424                 lnk = log ( nk );
    00425         }
    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:
    00446         mlnorm ( const RV &rv, const RV &rvc );
    00448         void set_parameters ( const  mat &A, const vec &mu0, const sq_T &R );
    00449 //      //!Generate one sample of the posterior
    00450 //      vec samplecond (const vec &cond, double &lik );
    00451 //      //!Generate matrix of samples of the posterior
    00452 //      mat samplecond (const vec &cond, vec &lik, int n );
    00454         void condition ( const vec &cond );
    00455
    00457         vec& _mu_const() {return mu_const;}
    00459         mat& _A() {return A;}
    00461         mat _R() {return epdf._R().to_mat();}
    00462
    00463         template<class sq_M>
    00464         friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M> &ml );
    00465 };
    00466
    00469 class mlstudent : public mlnorm<ldmat> {
    00470 protected:
    00471         ldmat Lambda;
    00472         ldmat &_R;
    00473         ldmat Re;
    00474 public:
    00475         mlstudent ( const RV &rv0, const RV &rvc0 ) :mlnorm<ldmat> ( rv0,rvc0 ),
    00476                         Lambda ( rv0.count() ),
    00477                         _R ( epdf._R() ) {}
    00478         void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0) {
    00479                 epdf.set_parameters ( zeros ( rv.count() ),Lambda );
    00480                 A = A0;
    00481                 mu_const = mu0;
    00482                 Re=R0;
    00483                 Lambda = Lambda0;
    00484         }
    00485         void condition ( const vec &cond ) {
    00486                 _mu = A*cond + mu_const;
    00487                 double zeta;
    00488                 //ugly hack!
    00489                 if ((cond.length()+1)==Lambda.rows()){
    00490                         zeta = Lambda.invqform ( concat(cond, vec_1(1.0)) );
    00491                 } else {
    00492                         zeta = Lambda.invqform ( cond );
    00493                 }
    00494                 _R = Re;
    00495                 _R*=( 1+zeta );// / ( nu ); << nu is in Re!!!!!!
    00496         };
    00497
    00498 };
    00508 class mgamma : public mEF {
    00509 protected:
    00511         egamma epdf;
    00513         double k;
    00515         vec* _beta;
    00516
    00517 public:
    00519         mgamma ( const RV &rv,const RV &rvc ): mEF ( rv,rvc ), epdf ( rv ) {vec* tmp; epdf._param ( tmp,_beta );ep=&epdf;};
    00521         void set_parameters ( double k );
    00522         void condition ( const vec &val ) {*_beta=k/val;};
    00523 };
    00524
    00534 class migamma : public mEF {
    00535         protected:
    00537                 eigamma epdf;
    00539                 double k;
    00541                 vec* _beta;
    00543                 vec* _alpha;
    00544
    00545         public:
    00547                 migamma ( const RV &rv,const RV &rvc ): mEF ( rv,rvc ), epdf ( rv ) {epdf._param ( _alpha,_beta );ep=&epdf;};
    00549                 void set_parameters ( double k0 ){k=k0;*_alpha=1.0/(k*k)+2;};
    00550                 void condition ( const vec &val ) {
    00551                         *_beta=elem_mult(val,(*_alpha-1));
    00552                 };
    00553 };
    00554
    00566 class mgamma_fix : public mgamma {
    00567 protected:
    00569         double l;
    00571         vec refl;
    00572 public:
    00574         mgamma_fix ( const RV &rv,const RV &rvc ) : mgamma ( rv,rvc ),refl ( rv.count() ) {};
    00576         void set_parameters ( double k0 , vec ref0, double l0 ) {
    00577                 mgamma::set_parameters ( k0 );
    00578                 refl=pow ( ref0,1.0-l0 );l=l0;
    00579         };
    00580
    00581         void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); *_beta=k/mean;};
    00582 };
    00583
    00584
    00597 class migamma_fix : public migamma {
    00598         protected:
    00600                 double l;
    00602                 vec refl;
    00603         public:
    00605                 migamma_fix ( const RV &rv,const RV &rvc ) : migamma ( rv,rvc ),refl ( rv.count() ) {};
    00607                 void set_parameters ( double k0 , vec ref0, double l0 ) {
    00608                         migamma::set_parameters ( k0 );
    00609                         refl=pow ( ref0,1.0-l0 );l=l0;
    00610                 };
    00611
    00612                 void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); migamma::condition(mean);};
    00613 };
    00615 enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
    00621 class eEmp: public epdf {
    00622 protected :
    00624         int n;
    00626         vec w;
    00628         Array<vec> samples;
    00629 public:
    00631         eEmp ( const RV &rv0 ,int n0 ) :epdf ( rv0 ),n ( n0 ),w ( n ),samples ( n ) {};
    00633         void set_parameters ( const vec &w0, const epdf* pdf0 );
    00635         void set_samples ( const epdf* pdf0 );
    00637         void set_n ( int n0, bool copy=true ){w.set_size(n0,copy);samples.set_size(n0,copy);};
    00639         vec& _w()  {return w;};
    00641         const vec& _w() const {return w;};
    00643         Array<vec>& _samples() {return samples;};
    00645         const Array<vec>& _samples() const {return samples;};
    00647         ivec resample ( RESAMPLING_METHOD method = SYSTEMATIC );
    00649         vec sample() const {it_error ( "Not implemented" );return 0;}
    00651         double evallog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
    00652         vec mean() const {
    00653                 vec pom=zeros ( rv.count() );
    00654                 for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
    00655                 return pom;
    00656         }
    00657         vec variance() const {
    00658                 vec pom=zeros ( rv.count() );
    00659                 for ( int i=0;i<n;i++ ) {pom+=pow(samples ( i ),2) *w ( i );}
    00660                 return pom-pow(mean(),2);
    00661         }
    00662 };
    00663
    00664
    00666
    00667 template<class sq_T>
    00668 enorm<sq_T>::enorm ( const RV &rv ) :eEF ( rv ), mu ( rv.count() ),R ( rv.count() ),dim ( rv.count() ) {};
    00669
    00670 template<class sq_T>
    00671 void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) {
    00672 //Fixme test dimensions of mu0 and R0;
    00673         mu = mu0;
    00674         R = R0;
    00675 };
    00676
    00677 template<class sq_T>
    00678 void enorm<sq_T>::dupdate ( mat &v, double nu ) {
    00679         //
    00680 };
    00681
    00682 // template<class sq_T>
    00683 // void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
    00684 //      //
    00685 // };
    00686
    00687 template<class sq_T>
    00688 vec enorm<sq_T>::sample() const {
    00689         vec x ( dim );
    00690         NorRNG.sample_vector ( dim,x );
    00691         vec smp = R.sqrt_mult ( x );
    00692
    00693         smp += mu;
    00694         return smp;
    00695 };
    00696
    00697 template<class sq_T>
    00698 mat enorm<sq_T>::sample ( int N ) const {
    00699         mat X ( dim,N );
    00700         vec x ( dim );
    00701         vec pom;
    00702         int i;
    00703
    00704         for ( i=0;i<N;i++ ) {
    00705                 NorRNG.sample_vector ( dim,x );
    00706                 pom = R.sqrt_mult ( x );
    00707                 pom +=mu;
    00708                 X.set_col ( i, pom );
    00709         }
    00710
    00711         return X;
    00712 };
    00713
    00714 // template<class sq_T>
    00715 // double enorm<sq_T>::eval ( const vec &val ) const {
    00716 //      double pdfl,e;
    00717 //      pdfl = evallog ( val );
    00718 //      e = exp ( pdfl );
    00719 //      return e;
    00720 // };
    00721
    00722 template<class sq_T>
    00723 double enorm<sq_T>::evallog_nn ( const vec &val ) const {
    00724         // 1.83787706640935 = log(2pi)
    00725         double tmp=-0.5* ( R.invqform ( mu-val ) );// - lognc();
    00726         return  tmp;
    00727 };
    00728
    00729 template<class sq_T>
    00730 inline double enorm<sq_T>::lognc () const {
    00731         // 1.83787706640935 = log(2pi)
    00732         double tmp=0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
    00733         return tmp;
    00734 };
    00735
    00736 template<class sq_T>
    00737 mlnorm<sq_T>::mlnorm ( const RV &rv0, const RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu ( epdf._mu() ) {
    00738         ep =&epdf;
    00739 }
    00740
    00741 template<class sq_T>
    00742 void mlnorm<sq_T>::set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 ) {
    00743         epdf.set_parameters ( zeros ( rv.count() ),R0 );
    00744         A = A0;
    00745         mu_const = mu0;
    00746 }
    00747
    00748 // template<class sq_T>
    00749 // vec mlnorm<sq_T>::samplecond (const  vec &cond, double &lik ) {
    00750 //      this->condition ( cond );
    00751 //      vec smp = epdf.sample();
    00752 //      lik = epdf.eval ( smp );
    00753 //      return smp;
    00754 // }
    00755
    00756 // template<class sq_T>
    00757 // mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) {
    00758 //      int i;
    00759 //      int dim = rv.count();
    00760 //      mat Smp ( dim,n );
    00761 //      vec smp ( dim );
    00762 //      this->condition ( cond );
    00763 //
    00764 //      for ( i=0; i<n; i++ ) {
    00765 //              smp = epdf.sample();
    00766 //              lik ( i ) = epdf.eval ( smp );
    00767 //              Smp.set_col ( i ,smp );
    00768 //      }
    00769 //
    00770 //      return Smp;
    00771 // }
    00772
    00773 template<class sq_T>
    00774 void mlnorm<sq_T>::condition ( const vec &cond ) {
    00775         _mu = A*cond + mu_const;
    00776 //R is already assigned;
    00777 }
    00778
    00779 template<class sq_T>
    00780 epdf* enorm<sq_T>::marginal ( const RV &rvn ) const {
    00781         ivec irvn = rvn.dataind ( rv );
    00782
    00783         sq_T Rn ( R,irvn );
    00784         enorm<sq_T>* tmp = new enorm<sq_T> ( rvn );
    00785         tmp->set_parameters ( mu ( irvn ), Rn );
    00786         return tmp;
    00787 }
    00788
    00789 template<class sq_T>
    00790 mpdf* enorm<sq_T>::condition ( const RV &rvn ) const {
    00791
    00792         RV rvc = rv.subt ( rvn );
    00793         it_assert_debug ( ( rvc.count() +rvn.count() ==rv.count() ),"wrong rvn" );
    00794         //Permutation vector of the new R
    00795         ivec irvn = rvn.dataind ( rv );
    00796         ivec irvc = rvc.dataind ( rv );
    00797         ivec perm=concat ( irvn , irvc );
    00798         sq_T Rn ( R,perm );
    00799
    00800         //fixme - could this be done in general for all sq_T?
    00801         mat S=Rn.to_mat();
    00802         //fixme
    00803         int n=rvn.count()-1;
    00804         int end=R.rows()-1;
    00805         mat S11 = S.get ( 0,n, 0, n );
    00806         mat S12 = S.get ( 0, n , rvn.count(), end );
    00807         mat S22 = S.get ( rvn.count(), end, rvn.count(), end );
    00808
    00809         vec mu1 = mu ( irvn );
    00810         vec mu2 = mu ( irvc );
    00811         mat A=S12*inv ( S22 );
    00812         sq_T R_n ( S11 - A *S12.T() );
    00813
    00814         mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( rvn,rvc );
    00815
    00816         tmp->set_parameters ( A,mu1-A*mu2,R_n );
    00817         return tmp;
    00818 }
    00819
    00821
    00822 template<class sq_T>
    00823 std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_T> &ml ) {
    00824         os << "A:"<< ml.A<<endl;
    00825         os << "mu:"<< ml.mu_const<<endl;
    00826         os << "R:" << ml.epdf._R().to_mat() <<endl;
    00827         return os;
    00828 };
    00829
    00830 #endif //EF_H