- Main Page
- Related Pages
- Classes
- Files

- File List
- File Members

# emix.h

Go to the documentation of this file.

    00001
    00013 #ifndef EMIX_H
    00014 #define EMIX_H
    00015
    00016 #define LOG2  0.69314718055995
    00017
    00018 #include "../shared_ptr.h"
    00019 #include "exp_family.h"
    00020
    00021 namespace bdm {
    00022
    00023 //this comes first because it is used inside emix!
    00024
    00037 class mratio: public mpdf {
    00038 protected:
    00040         const epdf* nom;
    00041
    00043         shared_ptr<epdf> den;
    00044
    00046         bool destroynom;
    00048         datalink_m2e dl;
    00050         epdf iepdf;
    00051 public:
    00054         mratio ( const epdf* nom0, const RV &rv, bool copy = false ) : mpdf ( ), dl ( ),iepdf() {
    00055                 // adjust rv and rvc
    00056                 rvc = nom0->_rv().subt ( rv );
    00057                 dimc = rvc._dsize();
    00058                 set_ep ( iepdf );
    00059                 iepdf.set_parameters ( rv._dsize() );
    00060                 iepdf.set_rv ( rv );
    00061
    00062                 //prepare data structures
    00063                 if ( copy ) {
    00064                         it_error ( "todo" );
    00065                         destroynom = true;
    00066                 } else {
    00067                         nom = nom0;
    00068                         destroynom = false;
    00069                 }
    00070                 it_assert_debug ( rvc.length() > 0, "Makes no sense to use this object!" );
    00071
    00072                 // build denominator
    00073                 den = nom->marginal ( rvc );
    00074                 dl.set_connection ( rv, rvc, nom0->_rv() );
    00075         };
    00076         double evallogcond ( const vec &val, const vec &cond ) {
    00077                 double tmp;
    00078                 vec nom_val ( dimension() + dimc );
    00079                 dl.pushup_cond ( nom_val, val, cond );
    00080                 tmp = exp ( nom->evallog ( nom_val ) - den->evallog ( cond ) );
    00081                 return tmp;
    00082         }
    00084         void ownnom() {
    00085                 destroynom = true;
    00086         }
    00088         ~mratio() {
    00089                 if ( destroynom ) {
    00090                         delete nom;
    00091                 }
    00092         }
    00093 };
    00094
    00105 class emix : public epdf {
    00106 protected:
    00108         vec w;
    00110         Array<shared_ptr<epdf> > Coms;
    00111
    00112 public:
    00114         emix ( ) : epdf ( ) {};
    00117         void set_parameters ( const vec &w, const Array<shared_ptr<epdf> > &Coms );
    00118
    00119         vec sample() const;
    00120         vec mean() const {
    00121                 int i;
    00122                 vec mu = zeros ( dim );
    00123                 for ( i = 0; i < w.length(); i++ ) {
    00124                         mu += w ( i ) * Coms ( i )->mean();
    00125                 }
    00126                 return mu;
    00127         }
    00128         vec variance() const {
    00129                 //non-central moment
    00130                 vec mom2 = zeros ( dim );
    00131                 for ( int i = 0; i < w.length(); i++ ) {
    00132                         mom2 += w ( i ) * ( Coms ( i )->variance() + pow ( Coms ( i )->mean(), 2 ) );
    00133                 }
    00134                 //central moment
    00135                 return mom2 - pow ( mean(), 2 );
    00136         }
    00137         double evallog ( const vec &val ) const {
    00138                 int i;
    00139                 double sum = 0.0;
    00140                 for ( i = 0; i < w.length(); i++ ) {
    00141                         sum += w ( i ) * exp ( Coms ( i )->evallog ( val ) );
    00142                 }
    00143                 if ( sum == 0.0 ) {
    00144                         sum = std::numeric_limits<double>::epsilon();
    00145                 }
    00146                 double tmp = log ( sum );
    00147                 it_assert_debug ( std::isfinite ( tmp ), "Infinite" );
    00148                 return tmp;
    00149         };
    00150         vec evallog_m ( const mat &Val ) const {
    00151                 vec x = zeros ( Val.cols() );
    00152                 for ( int i = 0; i < w.length(); i++ ) {
    00153                         x += w ( i ) * exp ( Coms ( i )->evallog_m ( Val ) );
    00154                 }
    00155                 return log ( x );
    00156         };
    00158         mat evallog_M ( const mat &Val ) const {
    00159                 mat X ( w.length(), Val.cols() );
    00160                 for ( int i = 0; i < w.length(); i++ ) {
    00161                         X.set_row ( i, w ( i ) *exp ( Coms ( i )->evallog_m ( Val ) ) );
    00162                 }
    00163                 return X;
    00164         };
    00165
    00166         shared_ptr<epdf> marginal ( const RV &rv ) const;
    00168         void marginal ( const RV &rv, emix &target ) const;
    00169         shared_ptr<mpdf> condition ( const RV &rv ) const;
    00170
    00171 //Access methods
    00173         vec& _w() {
    00174                 return w;
    00175         }
    00176
    00178         shared_ptr<epdf> _Coms ( int i ) {
    00179                 return Coms ( i );
    00180         }
    00181
    00182         void set_rv ( const RV &rv ) {
    00183                 epdf::set_rv ( rv );
    00184                 for ( int i = 0; i < Coms.length(); i++ ) {
    00185                         Coms ( i )->set_rv ( rv );
    00186                 }
    00187         }
    00188 };
    00189 SHAREDPTR( emix );
    00190
    00195 class egiwmix : public egiw {
    00196 protected:
    00198         vec w;
    00200         Array<egiw*> Coms;
    00202         bool destroyComs;
    00203 public:
    00205         egiwmix ( ) : egiw ( ) {};
    00206
    00209         void set_parameters ( const vec &w, const Array<egiw*> &Coms, bool copy = false );
    00210
    00212         vec mean() const;
    00213
    00215         vec sample() const;
    00216
    00218         vec variance() const;
    00219
    00220         // TODO!!! Defined to follow ANSI and/or for future development
    00221         void mean_mat ( mat &M, mat&R ) const {};
    00222         double evallog_nn ( const vec &val ) const {
    00223                 return 0;
    00224         };
    00225         double lognc () const {
    00226                 return 0;
    00227         }
    00228
    00229         shared_ptr<epdf> marginal ( const RV &rv ) const;
    00230         void marginal ( const RV &rv, emix &target ) const;
    00231
    00232 //Access methods
    00234         vec& _w() {
    00235                 return w;
    00236         }
    00237         virtual ~egiwmix() {
    00238                 if ( destroyComs ) {
    00239                         for ( int i = 0; i < Coms.length(); i++ ) {
    00240                                 delete Coms ( i );
    00241                         }
    00242                 }
    00243         }
    00245         void ownComs() {
    00246                 destroyComs = true;
    00247         }
    00248
    00250         egiw* _Coms ( int i ) {
    00251                 return Coms ( i );
    00252         }
    00253
    00254         void set_rv ( const RV &rv ) {
    00255                 egiw::set_rv ( rv );
    00256                 for ( int i = 0; i < Coms.length(); i++ ) {
    00257                         Coms ( i )->set_rv ( rv );
    00258                 }
    00259         }
    00260
    00262         egiw* approx();
    00263 };
    00264
    00273 class mprod: public mpdf {
    00274 private:
    00275         Array<shared_ptr<mpdf> > mpdfs;
    00276
    00277 protected:
    00279         Array<datalink_m2m*> dls;
    00280
    00282         epdf iepdf;
    00283
    00284 public:
    00286         mprod() { }
    00287
    00290         mprod ( const Array<shared_ptr<mpdf> > &mFacs ) {
    00291                 set_elements ( mFacs );
    00292         }
    00294         void set_elements (const Array<shared_ptr<mpdf> > &mFacs );
    00295
    00296         double evallogcond ( const vec &val, const vec &cond ) {
    00297                 int i;
    00298                 double res = 0.0;
    00299                 for ( i = mpdfs.length() - 1; i >= 0; i-- ) {
    00300                         /*                      if ( mpdfs(i)->_rvc().count() >0) {
    00301                                                         mpdfs ( i )->condition ( dls ( i )->get_cond ( val,cond ) );
    00302                                                 }
    00303                                                 // add logarithms
    00304                                                 res += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );*/
    00305                         res += mpdfs ( i )->evallogcond (
    00306                                    dls ( i )->pushdown ( val ),
    00307                                    dls ( i )->get_cond ( val, cond )
    00308                                );
    00309                 }
    00310                 return res;
    00311         }
    00312         vec evallogcond_m ( const mat &Dt, const vec &cond ) {
    00313                 vec tmp ( Dt.cols() );
    00314                 for ( int i = 0; i < Dt.cols(); i++ ) {
    00315                         tmp ( i ) = evallogcond ( Dt.get_col ( i ), cond );
    00316                 }
    00317                 return tmp;
    00318         };
    00319         vec evallogcond_m ( const Array<vec> &Dt, const vec &cond ) {
    00320                 vec tmp ( Dt.length() );
    00321                 for ( int i = 0; i < Dt.length(); i++ ) {
    00322                         tmp ( i ) = evallogcond ( Dt ( i ), cond );
    00323                 }
    00324                 return tmp;
    00325         };
    00326
    00327
    00328         //TODO smarter...
    00329         vec samplecond ( const vec &cond ) {
    00331                 vec smp = std::numeric_limits<double>::infinity() * ones ( dimension() );
    00332                 vec smpi;
    00333                 // Hard assumption here!!! We are going backwards, to assure that samples that are needed from smp are already generated!
    00334                 for ( int i = ( mpdfs.length() - 1 ); i >= 0; i-- ) {
    00335                         // generate contribution of this mpdf
    00336                         smpi = mpdfs(i)->samplecond(dls ( i )->get_cond ( smp , cond ));
    00337                         // copy contribution of this pdf into smp
    00338                         dls ( i )->pushup ( smp, smpi );
    00339                 }
    00340                 return smp;
    00341         }
    00342         mat samplecond ( const vec &cond,  int N ) {
    00343                 mat Smp ( dimension(), N );
    00344                 for ( int i = 0; i < N; i++ ) {
    00345                         Smp.set_col ( i, samplecond ( cond ) );
    00346                 }
    00347                 return Smp;
    00348         }
    00349
    00357         void from_setting ( const Setting &set ) {
    00358                 Array<shared_ptr<mpdf> > atmp; //temporary Array
    00359                 UI::get ( atmp, set, "mpdfs", UI::compulsory );
    00360                 set_elements ( atmp );
    00361         }
    00362
    00363 };
    00364 UIREGISTER ( mprod );
    00365 SHAREDPTR ( mprod );
    00366
    00368 class eprod: public epdf {
    00369 protected:
    00371         Array<const epdf*> epdfs;
    00373         Array<datalink*> dls;
    00374 public:
    00376         eprod () : epdfs ( 0 ), dls ( 0 ) {};
    00378         void set_parameters ( const Array<const epdf*> &epdfs0, bool named = true ) {
    00379                 epdfs = epdfs0;//.set_length ( epdfs0.length() );
    00380                 dls.set_length ( epdfs.length() );
    00381
    00382                 bool independent = true;
    00383                 if ( named ) {
    00384                         for ( int i = 0; i < epdfs.length(); i++ ) {
    00385                                 independent = rv.add ( epdfs ( i )->_rv() );
    00386                                 it_assert_debug ( independent == true, "eprod:: given components are not independent." );
    00387                         }
    00388                         dim = rv._dsize();
    00389                 } else {
    00390                         dim = 0;
    00391                         for ( int i = 0; i < epdfs.length(); i++ ) {
    00392                                 dim += epdfs ( i )->dimension();
    00393                         }
    00394                 }
    00395                 //
    00396                 int cumdim = 0;
    00397                 int dimi = 0;
    00398                 int i;
    00399                 for ( i = 0; i < epdfs.length(); i++ ) {
    00400                         dls ( i ) = new datalink;
    00401                         if ( named ) {
    00402                                 dls ( i )->set_connection ( epdfs ( i )->_rv() , rv );
    00403                         } else {
    00404                                 dimi = epdfs ( i )->dimension();
    00405                                 dls ( i )->set_connection ( dimi, dim, linspace ( cumdim, cumdim + dimi - 1 ) );
    00406                                 cumdim += dimi;
    00407                         }
    00408                 }
    00409         }
    00410
    00411         vec mean() const {
    00412                 vec tmp ( dim );
    00413                 for ( int i = 0; i < epdfs.length(); i++ ) {
    00414                         vec pom = epdfs ( i )->mean();
    00415                         dls ( i )->pushup ( tmp, pom );
    00416                 }
    00417                 return tmp;
    00418         }
    00419         vec variance() const {
    00420                 vec tmp ( dim ); //second moment
    00421                 for ( int i = 0; i < epdfs.length(); i++ ) {
    00422                         vec pom = epdfs ( i )->mean();
    00423                         dls ( i )->pushup ( tmp, pow ( pom, 2 ) );
    00424                 }
    00425                 return tmp - pow ( mean(), 2 );
    00426         }
    00427         vec sample() const {
    00428                 vec tmp ( dim );
    00429                 for ( int i = 0; i < epdfs.length(); i++ ) {
    00430                         vec pom = epdfs ( i )->sample();
    00431                         dls ( i )->pushup ( tmp, pom );
    00432                 }
    00433                 return tmp;
    00434         }
    00435         double evallog ( const vec &val ) const {
    00436                 double tmp = 0;
    00437                 for ( int i = 0; i < epdfs.length(); i++ ) {
    00438                         tmp += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );
    00439                 }
    00440                 it_assert_debug ( std::isfinite ( tmp ), "Infinite" );
    00441                 return tmp;
    00442         }
    00444         const epdf* operator () ( int i ) const {
    00445                 it_assert_debug ( i < epdfs.length(), "wrong index" );
    00446                 return epdfs ( i );
    00447         }
    00448
    00450         ~eprod() {
    00451                 for ( int i = 0; i < epdfs.length(); i++ ) {
    00452                         delete dls ( i );
    00453                 }
    00454         }
    00455 };
    00456
    00457
    00461 class mmix : public mpdf {
    00462 protected:
    00464         Array<shared_ptr<mpdf> > Coms;
    00466         vec w;
    00468         epdf dummy_epdf;
    00469 public:
    00471         mmix() : Coms(0), dummy_epdf() { set_ep(dummy_epdf);    }
    00472
    00474         void set_parameters ( const vec &w0, const Array<shared_ptr<mpdf> > &Coms0 ) {
    00476                 Coms = Coms0;
    00477                 w=w0;
    00478
    00479                 if (Coms0.length()>0){
    00480                         set_rv(Coms(0)->_rv());
    00481                         dummy_epdf.set_parameters(Coms(0)->_rv()._dsize());
    00482                         set_rvc(Coms(0)->_rvc());
    00483                         dimc = rvc._dsize();
    00484                 }
    00485         }
    00486         double evallogcond (const vec &dt, const vec &cond) {
    00487                 double ll=0.0;
    00488                 for (int i=0;i<Coms.length();i++){
    00489                         ll+=Coms(i)->evallogcond(dt,cond);
    00490                 }
    00491                 return ll;
    00492         }
    00493
    00494         vec samplecond (const vec &cond);
    00495
    00496 };
    00497
    00498 }
    00499 #endif //MX_H