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

# emix.h

Go to the documentation of this file.

    00001
    00013 #ifndef MX_H
    00014 #define MX_H
    00015
    00016 #include "libBM.h"
    00017 #include "libEF.h"
    00018 //#include <std>
    00019
    00020 namespace bdm{
    00021
    00022 //this comes first because it is used inside emix!
    00023
    00036 class mratio: public mpdf {
    00037 protected:
    00039         const epdf* nom;
    00041         epdf* den;
    00043         bool destroynom;
    00045         datalink_m2e dl;
    00046 public:
    00049         mratio ( const epdf* nom0, const RV &rv, bool copy=false ) :mpdf ( ), dl ( rv,rvc,nom0->_rv() ) {
    00050                 ep->set_rv(_rv());
    00051                 set_rvc(nom0->_rv().subt ( rv ) );
    00052                 if ( copy ) {
    00053 //                      nom = nom0->_copy_();
    00054                         it_error ( "todo" );
    00055                         destroynom=true;
    00056                 }
    00057                 else {
    00058                         nom = nom0;
    00059                         destroynom = false;
    00060                 }
    00061                 it_assert_debug ( rvc.length() >0,"Makes no sense to use this object!" );
    00062                 den = nom->marginal ( rvc );
    00063         };
    00064         double evallogcond ( const vec &val, const vec &cond ) {
    00065                 double tmp;
    00066                 vec nom_val ( ep->dimension() + dimc );
    00067                 dl.pushup_cond ( nom_val,val,cond );
    00068                 tmp = exp ( nom->evallog ( nom_val ) - den->evallog ( cond ) );
    00069                 it_assert_debug(std::isfinite(tmp),"Infinite value");
    00070                 return tmp;
    00071         }
    00073         void ownnom() {destroynom=true;}
    00075         ~mratio() {delete den; if ( destroynom ) {delete nom;}}
    00076 };
    00077
    00088 class emix : public epdf {
    00089 protected:
    00091         vec w;
    00093         Array<epdf*> Coms;
    00095         bool destroyComs;
    00096 public:
    00098         emix (  ) : epdf ( ) {};
    00101         void set_parameters ( const vec &w, const Array<epdf*> &Coms, bool copy=false );
    00102
    00103         vec sample() const;
    00104         vec mean() const {
    00105                 int i; vec mu = zeros ( dim);
    00106                 for ( i = 0;i < w.length();i++ ) {mu += w ( i ) * Coms ( i )->mean(); }
    00107                 return mu;
    00108         }
    00109         vec variance() const {
    00110                 //non-central moment
    00111                 vec mom2 = zeros(dim);
    00112                 for ( int i = 0;i < w.length();i++ ) {mom2 += w ( i ) * pow(Coms ( i )->mean(),2); }
    00113                 //central moment
    00114                 return mom2-pow(mean(),2);
    00115         }
    00116         double evallog ( const vec &val ) const {
    00117                 int i;
    00118                 double sum = 0.0;
    00119                 for ( i = 0;i < w.length();i++ ) {sum += w ( i ) * exp ( Coms ( i )->evallog ( val ) );}
    00120                 if (sum==0.0){sum=std::numeric_limits<double>::epsilon();}
    00121                 double tmp=log ( sum );
    00122                 it_assert_debug(std::isfinite(tmp),"Infinite");
    00123                 return tmp;
    00124         };
    00125         vec evallog_m ( const mat &Val ) const {
    00126                 vec x=zeros ( Val.cols() );
    00127                 for ( int i = 0; i < w.length(); i++ ) {
    00128                         x+= w ( i ) *exp ( Coms ( i )->evallog_m ( Val ) );
    00129                 }
    00130                 return log ( x );
    00131         };
    00133         mat evallog_M ( const mat &Val ) const {
    00134                 mat X ( w.length(), Val.cols() );
    00135                 for ( int i = 0; i < w.length(); i++ ) {
    00136                         X.set_row ( i, w ( i ) *exp ( Coms ( i )->evallog_m ( Val ) ) );
    00137                 }
    00138                 return X;
    00139         };
    00140
    00141         emix* marginal ( const RV &rv ) const;
    00142         mratio* condition ( const RV &rv ) const; //why not mratio!!
    00143
    00144 //Access methods
    00146         vec& _w() {return w;}
    00147         virtual ~emix() {if ( destroyComs ) {for ( int i=0;i<Coms.length();i++ ) {delete Coms ( i );}}}
    00149         void ownComs() {destroyComs=true;}
    00150
    00152         epdf* _Coms ( int i ) {return Coms ( i );}
    00153 };
    00154
    00163 class mprod: public compositepdf, public mpdf {
    00164 protected:
    00166         Array<epdf*> epdfs;
    00168         Array<datalink_m2m*> dls;
    00170         epdf dummy;
    00171 public:
    00174         mprod ( Array<mpdf*> mFacs ) : compositepdf ( mFacs ), mpdf (), epdfs ( n ), dls ( n ) {
    00175                 ep=&dummy;
    00176                 RV rv=getrv(true);
    00177                 set_rv(rv);dummy.set_parameters(rv._dsize());
    00178                 setrvc ( ep->_rv(),rvc );
    00179                 // rv and rvc established = > we can link them with mpdfs
    00180                 for ( int i = 0;i < n;i++ ) {
    00181                         dls ( i ) = new datalink_m2m ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), _rv(), _rvc() );
    00182                 }
    00183
    00184                 for ( int i=0;i<n;i++ ) {
    00185                         epdfs ( i ) =& ( mpdfs ( i )->_epdf() );
    00186                 }
    00187         };
    00188
    00189         double evallogcond ( const vec &val, const vec &cond ) {
    00190                 int i;
    00191                 double res = 0.0;
    00192                 for ( i = n - 1;i >= 0;i-- ) {
    00193                         /*                      if ( mpdfs(i)->_rvc().count() >0) {
    00194                                                         mpdfs ( i )->condition ( dls ( i )->get_cond ( val,cond ) );
    00195                                                 }
    00196                                                 // add logarithms
    00197                                                 res += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );*/
    00198                         res += mpdfs ( i )->evallogcond (
    00199                                    dls ( i )->pushdown ( val ),
    00200                                    dls ( i )->get_cond ( val, cond )
    00201                                );
    00202                 }
    00203                 return res;
    00204         }
    00205         //TODO smarter...
    00206         vec samplecond ( const vec &cond ) {
    00208                 vec smp= std::numeric_limits<double>::infinity() * ones ( ep->dimension() );
    00209                 vec smpi;
    00210                 // Hard assumption here!!! We are going backwards, to assure that samples that are needed from smp are already generated!
    00211                 for ( int i = ( n - 1 );i >= 0;i-- ) {
    00212                         if ( mpdfs ( i )->dimensionc() ) {
    00213                                 mpdfs ( i )->condition ( dls ( i )->get_cond ( smp ,cond ) ); // smp is val here!!
    00214                         }
    00215                         smpi = epdfs ( i )->sample();
    00216                         // copy contribution of this pdf into smp
    00217                         dls ( i )->pushup ( smp, smpi );
    00218                         // add ith likelihood contribution
    00219                 }
    00220                 return smp;
    00221         }
    00222         mat samplecond ( const vec &cond,  int N ) {
    00223                 mat Smp ( dimension(),N );
    00224                 for ( int i=0;i<N;i++ ) {Smp.set_col ( i,samplecond ( cond ) );}
    00225                 return Smp;
    00226         }
    00227
    00228         ~mprod() {};
    00229 };
    00230
    00232 class eprod: public epdf {
    00233 protected:
    00235         Array<const epdf*> epdfs;
    00237         Array<datalink*> dls;
    00238 public:
    00239         eprod ( const Array<const epdf*> epdfs0 ) : epdf (  ),epdfs ( epdfs0 ),dls ( epdfs.length() ) {
    00240                 bool independent=true;
    00241                 for ( int i=0;i<epdfs.length();i++ ) {
    00242                         independent=rv.add ( epdfs ( i )->_rv() );
    00243                         it_assert_debug ( independent==true, "eprod:: given components are not independent ." );
    00244                 }
    00245                 for ( int i=0;i<epdfs.length();i++ ) {
    00246                         dls ( i ) = new datalink ( epdfs ( i )->_rv() , rv );
    00247                 }
    00248         }
    00249
    00250         vec mean() const {
    00251                 vec tmp ( dim );
    00252                 for ( int i=0;i<epdfs.length();i++ ) {
    00253                         vec pom = epdfs ( i )->mean();
    00254                         dls ( i )->pushup ( tmp, pom );
    00255                 }
    00256                 return tmp;
    00257         }
    00258         vec variance() const {
    00259                 vec tmp ( dim ); //second moment
    00260                 for ( int i=0;i<epdfs.length();i++ ) {
    00261                         vec pom = epdfs ( i )->mean();
    00262                         dls ( i )->pushup ( tmp, pow(pom,2) );
    00263                 }
    00264                 return tmp-pow(mean(),2);
    00265         }
    00266         vec sample() const {
    00267                 vec tmp ( dim );
    00268                 for ( int i=0;i<epdfs.length();i++ ) {
    00269                         vec pom = epdfs ( i )->sample();
    00270                         dls ( i )->pushup ( tmp, pom );
    00271                 }
    00272                 return tmp;
    00273         }
    00274         double evallog ( const vec &val ) const {
    00275                 double tmp=0;
    00276                 for ( int i=0;i<epdfs.length();i++ ) {
    00277                         tmp+=epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );
    00278                 }
    00279                 it_assert_debug(std::isfinite(tmp),"Infinite");
    00280                 return tmp;
    00281         }
    00283         const epdf* operator () ( int i ) const {it_assert_debug ( i<epdfs.length(),"wrong index" );return epdfs ( i );}
    00284
    00286         ~eprod() {for ( int i=0;i<epdfs.length();i++ ) {delete dls ( i );}}
    00287 };
    00288
    00289
    00293 class mmix : public mpdf {
    00294 protected:
    00296         Array<mpdf*> Coms;
    00298         emix Epdf;
    00299 public:
    00301         mmix ( ) : mpdf ( ), Epdf () {ep = &Epdf;};
    00303         void set_parameters ( const vec &w, const Array<mpdf*> &Coms ) {
    00304                 Array<epdf*> Eps ( Coms.length() );
    00305
    00306                 for ( int i = 0;i < Coms.length();i++ ) {
    00307                         Eps ( i ) = & ( Coms ( i )->_epdf() );
    00308                 }
    00309                 Epdf.set_parameters ( w, Eps );
    00310         };
    00311
    00312         void condition ( const vec &cond ) {
    00313                 for ( int i = 0;i < Coms.length();i++ ) {Coms ( i )->condition ( cond );}
    00314         };
    00315 };
    00316
    00317 }
    00318 #endif //MX_H