Changeset 299 for bdm

Show
Ignore:
Timestamp:
03/19/09 15:38:55 (15 years ago)
Author:
smidl
Message:

merger - cull de sac :-(

Location:
bdm
Files:
3 modified

Legend:

Unmodified
Added
Removed

  • bdm/estim/merger.cpp

    r286 r299  
    33#include "arx.h" 
    44 
    5 namespace bdm{ 
    6 vec merger::lognorm_merge ( mat &lW ) { 
    7         int nu=lW.rows(); 
    8         vec mu = sum ( lW ) /nu; //mean of logs 
    9         // vec lam = sum ( pow ( lW,2 ) )-nu*pow ( mu,2 ); ======= numerically unsafe! 
    10         vec lam = sum ( pow ( lW-outer_product ( ones ( lW.rows() ),mu ),2 ) ); 
    11         double coef=0.0; 
    12         vec sq2bl=sqrt ( 2*beta*lam ); //this term is everywhere 
    13         switch ( nu ) { 
    14                 case 2: 
    15                         coef= ( 1-0.5*sqrt ( ( 4.0*beta-3.0 ) /beta ) ); 
    16                         return  coef*sq2bl + mu ; 
    17                         break; 
    18                 case 3://Ratio of Bessel 
    19                         coef = sqrt ( ( 3*beta-2 ) /3*beta ); 
    20                         return log ( besselk ( 0,sq2bl*coef ) ) - log ( besselk ( 0,sq2bl ) ) +  mu; 
    21                         break; 
    22                 case 4: 
    23                         break; 
    24                 default: // Approximate conditional density 
    25                         break; 
     5namespace bdm 
     6{ 
     7        vec merger::lognorm_merge ( mat &lW ) 
     8        { 
     9                int nu=lW.rows(); 
     10                vec mu = sum ( lW ) /nu; //mean of logs 
     11                // vec lam = sum ( pow ( lW,2 ) )-nu*pow ( mu,2 ); ======= numerically unsafe! 
     12                vec lam = sum ( pow ( lW-outer_product ( ones ( lW.rows() ),mu ),2 ) ); 
     13                double coef=0.0; 
     14                vec sq2bl=sqrt ( 2*beta*lam ); //this term is everywhere 
     15                switch ( nu ) 
     16                { 
     17                        case 2: 
     18                                coef= ( 1-0.5*sqrt ( ( 4.0*beta-3.0 ) /beta ) ); 
     19                                return  coef*sq2bl + mu ; 
     20                                break; 
     21                        case 3://Ratio of Bessel 
     22                                coef = sqrt ( ( 3*beta-2 ) /3*beta ); 
     23                                return log ( besselk ( 0,sq2bl*coef ) ) - log ( besselk ( 0,sq2bl ) ) +  mu; 
     24                                break; 
     25                        case 4: 
     26                                break; 
     27                        default: // Approximate conditional density 
     28                                break; 
     29                } 
     30                return vec ( 0 ); 
    2631        } 
    27         return vec ( 0 ); 
     32 
     33        void merger::merge ( const epdf* g0 ) 
     34        { 
     35 
     36                it_assert_debug ( rv.equal ( g0->_rv() ),"Incompatible g0" ); 
     37                //Empirical density - samples 
     38                if (!fix_smp) 
     39                        eSmp.set_statistics ( ones ( Ns ), g0 ); 
     40                 
     41                Array<vec> &Smp = eSmp._samples(); //aux 
     42                vec &w = eSmp._w(); //aux 
     43 
     44                mat Smp_ex =ones ( dim +1,Ns ); // Extended samples for the ARX model - the last row is ones 
     45                for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );} 
     46 
     47                if ( DBG )      *dbg << Name ( "Smp_0" ) << Smp_ex; 
     48 
     49                // Stuff for merging 
     50                vec lw_src ( Ns ); 
     51                vec lw_mix ( Ns ); 
     52                vec lw ( Ns ); 
     53                mat lW=zeros ( n,Ns ); 
     54                vec vec0 ( 0 ); 
     55 
     56                // Initial component in the mixture model 
     57                mat V0=1e-8*eye ( dim +1 ); 
     58                ARX A0; A0.set_statistics ( dim, V0 ); //initial guess of Mix: 
     59 
     60                Mix.init ( &A0, Smp_ex, Nc ); 
     61                //Preserve initial mixture for repetitive estimation via flattening 
     62                MixEF Mix_init ( Mix ); 
     63 
     64                // ============= MAIN LOOP ================== 
     65                bool converged=false; 
     66                int niter = 0; 
     67                char str[100]; 
     68 
     69                emix* Mpred=Mix.epredictor ( ); 
     70                vec Mix_pdf ( Ns ); 
     71                while ( !converged ) 
     72                { 
     73                        //Re-estimate Mix 
     74                        //Re-Initialize Mixture model 
     75                        Mix.flatten ( &Mix_init ); 
     76                        Mix.bayesB ( Smp_ex, w*Ns ); 
     77                        delete Mpred; 
     78                        Mpred = Mix.epredictor ( ); // Allocation => must be deleted at the end!! 
     79                        Mpred->set_rv ( rv ); //the predictor predicts rv of this merger 
     80 
     81                        // This will be active only later in iterations!!! 
     82                        if ((!fix_smp)&( 1./sum_sqr ( w ) <0.5*Ns )) 
     83                        { 
     84                                // Generate new samples 
     85                                eSmp.set_samples ( Mpred ); 
     86                                for ( int i=0;i<Ns;i++ ) 
     87                                { 
     88                                        //////////// !!!!!!!!!!!!! 
     89                                        if ( Smp ( i ) ( 2 ) <0 ) {Smp ( i ) ( 2 ) = 0.01; } 
     90                                        set_col_part ( Smp_ex,i,Smp ( i ) ); 
     91                                } 
     92                                if ( 0 ) 
     93                                { 
     94                                        cout<<"Resampling =" << 1./sum_sqr ( w ) << endl; 
     95                                        cout << sum ( Smp_ex,2 ) /Ns <<endl; 
     96                                        cout << Smp_ex*Smp_ex.T() /Ns << endl; 
     97                                } 
     98                        } 
     99                        if ( DBG ) 
     100                        { 
     101                                sprintf ( str,"Mpred_mean%d",niter ); 
     102                                *dbg << Name ( str ) << Mpred->mean(); 
     103 
     104 
     105                                sprintf ( str,"Mpdf%d",niter ); 
     106                                for ( int i=0;i<Ns;i++ ) {Mix_pdf ( i ) = Mix.logpred ( Smp_ex.get_col ( i ) );} 
     107                                *dbg << Name ( str ) << Mix_pdf; 
     108 
     109                                sprintf ( str,"Smp%d",niter ); 
     110                                *dbg << Name ( str ) << Smp_ex; 
     111 
     112                        } 
     113                        //Importace weighting 
     114                        for ( int i=0;i<n;i++ ) 
     115                        { 
     116                                lw_src=0.0; 
     117                                //======== Same RVs =========== 
     118                                //Split according to dependency in rvs 
     119                                if ( mpdfs ( i )->dimension() ==dim ) 
     120                                { 
     121                                        // no need for conditioning or marginalization 
     122                                        for ( int j=0;j<Ns; j++ )   // Smp is Array<> => for cycle 
     123                                        { 
     124                                                lw_src ( j ) =mpdfs ( i )->_epdf().evallog ( Smp ( j ) ); 
     125                                        } 
     126                                } 
     127                                else 
     128                                { 
     129                                        // compute likelihood of marginal on the conditional variable 
     130                                        if ( mpdfs ( i )->dimensionc() >0 ) 
     131                                        { 
     132                                                // Make marginal on rvc_i 
     133                                                epdf* tmp_marg = Mpred->marginal ( mpdfs ( i )->_rvc() ); 
     134                                                //compute vector of lw_src 
     135                                                for ( int k=0;k<Ns;k++ ) 
     136                                                { 
     137                                                        // Here val of tmp_marg = cond of mpdfs(i) ==> calling dls->get_cond 
     138                                                        lw_src ( k ) += tmp_marg->evallog ( dls ( i )->get_cond ( Smp ( k ) ) ); 
     139                                                } 
     140                                                delete tmp_marg; 
     141 
     142//                                      sprintf ( str,"marg%d",niter ); 
     143//                                      *dbg << Name ( str ) << lw_src; 
     144 
     145                                        } 
     146                                        // Compute likelihood of the missing variable 
     147                                        if ( dim > ( mpdfs ( i )->dimension() + mpdfs ( i )->dimensionc() ) ) 
     148                                        { 
     149                                                /////////////// 
     150                                                // There are variales unknown to mpdfs(i) : rvzs 
     151                                                mpdf* tmp_cond = Mpred->condition ( rvzs ( i ) ); 
     152                                                // Compute likelihood 
     153                                                vec lw_dbg=lw_src; 
     154                                                for ( int k= 0; k<Ns; k++ ) 
     155                                                { 
     156                                                        lw_src ( k ) += log ( 
     157                                                                            tmp_cond->evallogcond ( 
     158                                                                                zdls ( i )->pushdown ( Smp ( k ) ), 
     159                                                                                zdls ( i )->get_cond ( Smp ( k ) ) ) ); 
     160                                                        if ( !std::isfinite ( lw_src ( k ) ) ) 
     161                                                        { 
     162                                                                lw_src ( k ) = -1e16; cout << "!"; 
     163                                                        } 
     164                                                } 
     165                                                delete tmp_cond; 
     166                                        } 
     167                                        // Compute likelihood of the partial source 
     168                                        for ( int k= 0; k<Ns; k++ ) 
     169                                        { 
     170                                                mpdfs ( i )->condition ( dls ( i )->get_cond ( Smp ( k ) ) ); 
     171                                                lw_src ( k ) += mpdfs ( i )->_epdf().evallog ( dls ( i )->pushdown ( Smp ( k ) ) ); 
     172                                        } 
     173 
     174                                } 
     175//                      it_assert_debug(std::isfinite(sum(lw_src)),"bad"); 
     176                                lW.set_row ( i, lw_src ); // do not divide by mix 
     177                        } 
     178                        //Importance of the mixture 
     179                        for ( int j=0;j<Ns;j++ ) 
     180                        { 
     181                                lw_mix ( j ) =Mix.logpred ( Smp_ex.get_col ( j ) ); 
     182                        } 
     183                        lw = lognorm_merge ( lW ); //merge 
     184 
     185                        if ( DBG ) 
     186                        { 
     187                                sprintf ( str,"lW%d",niter ); 
     188                                *dbg << Name ( str ) << lW; 
     189                                sprintf ( str,"w%d",niter ); 
     190                                *dbg << Name ( str ) << w; 
     191                                sprintf ( str,"lw_m%d",niter ); 
     192                                *dbg << Name ( str ) << lw_mix; 
     193                        } 
     194                        //Importance weighting 
     195                        if (!fix_smp) 
     196                                lw -=  lw_mix; // hoping that it is not numerically sensitive... 
     197                        w = exp ( lw-max ( lw ) ); 
     198                        //renormalize 
     199                        double sumw =sum ( w ); 
     200                        if ( std::isfinite ( sumw ) ) 
     201                        { 
     202                                w = w/sumw; 
     203                        } 
     204                        else 
     205                        { 
     206                                it_file itf ( "merg_err.it" ); 
     207                                itf << Name ( "w" ) << w; 
     208                        } 
     209 
     210                        if ( DBG ) 
     211                        { 
     212                                sprintf ( str,"w_is_%d",niter ); 
     213                                *dbg << Name ( str ) << w; 
     214                        } 
     215                        // ==== stopping rule === 
     216                        niter++; 
     217                        converged = ( niter>20 ); 
     218                } 
     219                delete Mpred; 
     220//              cout << endl; 
     221 
     222        } 
     223 
    28224} 
    29  
    30 void merger::merge ( const epdf* g0 ) { 
    31 //      it_file dbg ( "merger_debug.it" ); 
    32  
    33         it_assert_debug ( rv.equal ( g0->_rv() ),"Incompatible g0" ); 
    34         //Empirical density - samples 
    35         eSmp.set_statistics ( ones ( Ns ), g0 ); 
    36         Array<vec> &Smp = eSmp._samples(); //aux 
    37         vec &w = eSmp._w(); //aux 
    38  
    39         mat Smp_ex =ones ( dim +1,Ns ); // Extended samples for the ARX model - the last row is ones 
    40         for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );} 
    41  
    42 //      dbg << Name ( "Smp_0" ) << Smp_ex; 
    43  
    44         // Stuff for merging 
    45         vec lw_src ( Ns ); 
    46         vec lw_mix ( Ns ); 
    47         vec lw ( Ns ); 
    48         mat lW=zeros ( n,Ns ); 
    49         vec vec0 ( 0 ); 
    50  
    51         // Initial component in the mixture model 
    52         mat V0=1e-8*eye ( dim +1 ); 
    53         ARX A0; A0.set_statistics(dim, V0); //initial guess of Mix:  
    54  
    55         Mix.init ( &A0, Smp_ex, Nc ); 
    56         //Preserve initial mixture for repetitive estimation via flattening 
    57         MixEF Mix_init ( Mix ); 
    58  
    59         // ============= MAIN LOOP ================== 
    60         bool converged=false; 
    61         int niter = 0; 
    62         char str[100]; 
    63  
    64         emix* Mpred=Mix.epredictor (  ); 
    65         vec Mix_pdf ( Ns ); 
    66         while ( !converged ) { 
    67                 //Re-estimate Mix 
    68                 //Re-Initialize Mixture model 
    69                 Mix.flatten ( &Mix_init ); 
    70                 Mix.bayesB ( Smp_ex, w*Ns ); 
    71                 delete Mpred; 
    72                 Mpred = Mix.epredictor ( ); // Allocation => must be deleted at the end!! 
    73                 Mpred->set_rv(rv); //the predictor predicts rv of this merger 
    74  
    75                 // This will be active only later in iterations!!! 
    76                 if ( 1./sum_sqr ( w ) <0.5*Ns ) { 
    77                         // Generate new samples 
    78                         eSmp.set_samples ( Mpred ); 
    79                         for ( int i=0;i<Ns;i++ ) { 
    80                                 //////////// !!!!!!!!!!!!! 
    81                                 if ( Smp ( i ) ( 2 ) <0 ) {Smp ( i ) ( 2 ) = 0.01; } 
    82                                 set_col_part ( Smp_ex,i,Smp ( i ) ); 
    83                         } 
    84                         if(0){cout<<"Resampling =" << 1./sum_sqr ( w ) << endl; 
    85                         cout << sum ( Smp_ex,2 ) /Ns <<endl; 
    86                         cout << Smp_ex*Smp_ex.T() /Ns << endl;} 
    87                 } 
    88 //              sprintf ( str,"Mpred_mean%d",niter ); 
    89 //              dbg << Name ( str ) << Mpred->mean(); 
    90  
    91  
    92 //              sprintf ( str,"Mpdf%d",niter ); 
    93 //              for ( int i=0;i<Ns;i++ ) {Mix_pdf ( i ) = Mix.logpred ( Smp_ex.get_col ( i ) );} 
    94 //              dbg << Name ( str ) << Mix_pdf; 
    95  
    96 //              sprintf ( str,"Smp%d",niter ); 
    97 //              dbg << Name ( str ) << Smp_ex; 
    98  
    99                 //Importace weighting 
    100                 for ( int i=0;i<n;i++ ) { 
    101                         lw_src=0.0; 
    102                         //======== Same RVs =========== 
    103                         //Split according to dependency in rvs 
    104                         if ( mpdfs ( i )->dimension() ==dim ) { 
    105                                 // no need for conditioning or marginalization 
    106                                 for ( int j=0;j<Ns; j++ ) { // Smp is Array<> => for cycle 
    107                                         lw_src ( j ) =mpdfs ( i )->_epdf().evallog ( Smp ( j ) ); 
    108                                 } 
    109                         } 
    110                         else { 
    111                                 // compute likelihood of marginal on the conditional variable 
    112                                 if ( mpdfs ( i )->dimensionc() >0 ) { 
    113                                         // Make marginal on rvc_i 
    114                                         epdf* tmp_marg = Mpred->marginal ( mpdfs ( i )->_rvc() ); 
    115                                         //compute vector of lw_src 
    116                                         for ( int k=0;k<Ns;k++ ) { 
    117                                                 // Here val of tmp_marg = cond of mpdfs(i) ==> calling dls->get_cond 
    118                                                 lw_src ( k ) += tmp_marg->evallog ( dls ( i )->get_cond ( Smp ( k ) ) ); 
    119                                         } 
    120                                         delete tmp_marg; 
    121  
    122 //                                      sprintf ( str,"marg%d",niter ); 
    123 //                                      dbg << Name ( str ) << lw_src; 
    124  
    125                                 } 
    126                                 // Compute likelihood of the missing variable 
    127                                 if ( dim > ( mpdfs ( i )->dimension() + mpdfs ( i )->dimensionc() ) ) { 
    128                                         /////////////// 
    129                                         // There are variales unknown to mpdfs(i) : rvzs 
    130                                         mpdf* tmp_cond = Mpred->condition ( rvzs ( i ) ); 
    131                                         // Compute likelihood 
    132                                         vec lw_dbg=lw_src; 
    133                                         for ( int k= 0; k<Ns; k++ ) { 
    134                                                 lw_src ( k ) += log ( 
    135                                                                     tmp_cond->evallogcond ( 
    136                                                                         zdls ( i )->pushdown ( Smp ( k ) ), 
    137                                                                         zdls ( i )->get_cond ( Smp ( k ) ) ) ); 
    138                                                 if ( !std::isfinite ( lw_src ( k ) ) ) { 
    139                                                         lw_src ( k ) = -1e16; cout << "!"; 
    140                                                 } 
    141                                         } 
    142                                         delete tmp_cond; 
    143                                 } 
    144                                 // Compute likelihood of the partial source 
    145                                 for ( int k= 0; k<Ns; k++ ) { 
    146                                         mpdfs ( i )->condition ( dls ( i )->get_cond ( Smp ( k ) ) ); 
    147                                         lw_src ( k ) += mpdfs ( i )->_epdf().evallog ( dls ( i )->pushdown ( Smp ( k ) ) ); 
    148                                 } 
    149  
    150                         } 
    151 //                      it_assert_debug(std::isfinite(sum(lw_src)),"bad"); 
    152                         lW.set_row ( i, lw_src ); // do not divide by mix 
    153                 } 
    154                 //Importance of the mixture 
    155                 for ( int j=0;j<Ns;j++ ) { 
    156                         lw_mix ( j ) =Mix.logpred ( Smp_ex.get_col ( j ) ); 
    157                 } 
    158 //              sprintf ( str,"lW%d",niter ); 
    159 //              dbg << Name ( str ) << lW; 
    160  
    161                 lw = lognorm_merge ( lW ); //merge 
    162  
    163 //              sprintf ( str,"w%d",niter ); 
    164 //              dbg << Name ( str ) << w; 
    165 //              sprintf ( str,"lw_m%d",niter ); 
    166 //              dbg << Name ( str ) << lw_mix; 
    167  
    168                 //Importance weighting 
    169                 lw -=  lw_mix; // hoping that it is not numerically sensitive... 
    170                 w = exp ( lw-max ( lw ) ); 
    171                 //renormalize 
    172                 double sumw =sum ( w ); 
    173                 if ( std::isfinite ( sumw ) ) { 
    174                         w = w/sumw; 
    175                 } 
    176                 else { 
    177                         it_file itf ( "merg_err.it" ); 
    178                         itf << Name ( "w" ) << w; 
    179                 } 
    180  
    181 //              sprintf ( str,"w_is_%d",niter ); 
    182 //              dbg << Name ( str ) << w; 
    183  
    184 //              eSmp.resample(); // So that it can be used in bayes 
    185 //              for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );} 
    186  
    187 //              sprintf ( str,"Smp_res%d",niter ); 
    188 //              dbg << Name ( str ) << Smp; 
    189  
    190                 // ==== stopping rule === 
    191                 niter++; 
    192                 converged = ( niter>20 ); 
    193         } 
    194         delete Mpred; 
    195         cout << endl; 
    196  
    197 } 
    198  
    199 } 

  • bdm/estim/merger.h

    r286 r299  
    1717#include "mixef.h" 
    1818 
    19 namespace bdm{ 
    20 using std::string; 
    21  
    22 /*! 
    23 @brief Function for general combination of pdfs 
    24  
    25 Mixtures of Gaussian densities are used internally. Switching to other densities should be trivial. 
    26 */ 
    27  
    28 class merger : public compositepdf, public epdf { 
    29 protected: 
    30         //!Internal mixture of EF models 
    31         MixEF Mix; 
    32         //! Data link for each mpdf in mpdfs 
    33         Array<datalink_m2e*> dls; 
    34         //! Array of rvs that are not modelled by mpdfs at all (aux) 
    35         Array<RV> rvzs; 
    36         //! Data Links of rv0 mpdfs - these will be conditioned the (rv,rvc) of mpdfs 
    37         Array<datalink_m2e*> zdls; 
    38  
    39         //!Number of samples used in approximation 
    40         int Ns; 
    41         //!Number of components in a mixture 
    42         int Nc; 
    43         //!Prior on the log-normal merging model 
    44         double beta; 
    45         //! Projection to empirical density 
    46         eEmp eSmp; 
    47  
    48 public: 
     19namespace bdm 
     20{ 
     21        using std::string; 
     22 
     23        /*! 
     24        @brief Function for general combination of pdfs 
     25 
     26        Mixtures of Gaussian densities are used internally. Switching to other densities should be trivial. 
     27        */ 
     28 
     29        class merger : public compositepdf, public epdf 
     30        { 
     31                protected: 
     32                        //!Internal mixture of EF models 
     33                        MixEF Mix; 
     34                        //! Data link for each mpdf in mpdfs 
     35                        Array<datalink_m2e*> dls; 
     36                        //! Array of rvs that are not modelled by mpdfs at all (aux) 
     37                        Array<RV> rvzs; 
     38                        //! Data Links of rv0 mpdfs - these will be conditioned the (rv,rvc) of mpdfs 
     39                        Array<datalink_m2e*> zdls; 
     40 
     41                        //!Number of samples used in approximation 
     42                        int Ns; 
     43                        //!Number of components in a mixture 
     44                        int Nc; 
     45                        //!Prior on the log-normal merging model 
     46                        double beta; 
     47                        //! Projection to empirical density 
     48                        eEmp eSmp; 
     49 
     50                        //! debug or not debug 
     51                        bool DBG; 
     52                        //! debugging file 
     53                        it_file* dbg; 
     54                        //! Flag if the samples are fixed or not 
     55                        bool fix_smp; 
     56                public: 
    4957//!Default constructor 
    50         merger ( const Array<mpdf*> &S ) : 
    51                         compositepdf ( S ), epdf ( ), 
    52                         Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ), dls ( n ), rvzs ( n ), zdls ( n ), eSmp() { 
    53                 RV ztmp; 
    54                 rv = getrv(false); 
    55                 dim = rv._dsize(); 
    56                 // Extend rv by rvc! 
    57                 RV rvc; setrvc ( rv,rvc ); 
    58                 rv.add ( rvc ); 
    59                 for ( int i=0;i<n;i++ ) { 
    60                         //Establich connection between mpdfs and merger 
    61                         dls ( i ) = new datalink_m2e;dls(i)->set_connection( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), rv ); 
    62                         // find out what is missing in each mpdf 
    63                         ztmp= mpdfs ( i )->_rv(); 
    64                         ztmp.add ( mpdfs ( i )->_rvc() ); 
    65                         rvzs ( i ) =rv.subt ( ztmp ); 
    66                         zdls ( i ) = new datalink_m2e; zdls(i)->set_connection ( rvzs ( i ), ztmp, rv ) ; 
    67                 }; 
    68                 //Set Default values of parameters 
    69                 beta=2.0; 
    70                 Ns=100; 
    71                 Nc=10; 
    72                 Mix.set_method ( EM ); 
    73         } 
     58                        merger ( const Array<mpdf*> &S ) : 
     59                                        compositepdf ( S ), epdf ( ), 
     60                                        Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ), dls ( n ), rvzs ( n ), zdls ( n ), eSmp() 
     61                        { 
     62                                RV ztmp; 
     63                                rv = getrv ( false ); 
     64                                RV rvc; setrvc ( rv,rvc ); // Extend rv by rvc! 
     65                                rv.add ( rvc ); 
     66                                // get dimension 
     67                                dim = rv._dsize(); 
     68 
     69                                for ( int i=0;i<n;i++ ) 
     70                                { 
     71                                        //Establich connection between mpdfs and merger 
     72                                        dls ( i ) = new datalink_m2e;dls ( i )->set_connection ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), rv ); 
     73                                        // find out what is missing in each mpdf 
     74                                        ztmp= mpdfs ( i )->_rv(); 
     75                                        ztmp.add ( mpdfs ( i )->_rvc() ); 
     76                                        rvzs ( i ) =rv.subt ( ztmp ); 
     77                                        zdls ( i ) = new datalink_m2e; zdls ( i )->set_connection ( rvzs ( i ), ztmp, rv ) ; 
     78                                }; 
     79                                //Set Default values of parameters 
     80                                beta=2.0; 
     81                                Ns=100; 
     82                                Nc=10; 
     83                                Mix.set_method ( EM ); 
     84                                DBG = false; 
     85                                fix_smp = false; 
     86                        } 
     87                        //! set debug file 
     88                        void debug_file ( const string fname ) { if ( DBG ) delete dbg; dbg = new it_file ( fname ); if ( dbg ) DBG=true;} 
    7489//! Set internal parameters used in approximation 
    75         void set_parameters ( double beta0, int Ns0, int Nc0 ) {beta=beta0;Ns=Ns0;Nc=Nc0;eSmp.set_parameters(Ns0,false);} 
     90                        void set_parameters ( double beta0, int Ns0, int Nc0 ) {beta=beta0;Ns=Ns0;Nc=Nc0;eSmp.set_parameters ( Ns0,false );} 
     91                        void set_grid ( Array<vec> &XYZ ) 
     92                        { 
     93                                int dim=XYZ.length(); ivec szs ( dim ); 
     94                                for(int i=0; i<dim;i++){szs=XYZ(i).length();} 
     95                                Ns=prod(szs); 
     96                                eSmp.set_parameters(Ns,false); 
     97                                Array<vec> &samples=eSmp._samples(); 
     98                                eSmp._w()=ones(Ns)/Ns; 
     99                                                 
     100                                //set samples 
     101                                ivec is=zeros_i(dim);//indeces of dimensions in for cycle; 
     102                                vec smpi(dim); 
     103                                for(int i=0; i<Ns; i++){ 
     104                                        for(int j=0; j<dim; j++){smpi(j)=XYZ(j)(is(j)); /* jty vector*/ } 
     105                                        samples(i)=smpi; 
     106                                        // shift indeces 
     107                                        for (int j=0;j<dim;j++){ 
     108                                                if (is(j)==szs(j)-1) { //j-th index is full 
     109                                                        is(j)=0; //shift back 
     110                                                        is(j+1)++; //increase th next dimension; 
     111                                                        if (is(j+1)<szs(j+1)-1) break; 
     112                                                } else { 
     113                                                        is(j)++; break; 
     114                                                } 
     115                                        } 
     116                                } 
     117                                 
     118                                fix_smp = true; 
     119                        } 
    76120//!Initialize the proposal density. This function must be called before merge()! 
    77         void init() { ////////////// NOT FINISHED 
    78                 Array<vec> Smps ( n ); 
    79                 //Gibbs sampling 
    80                 for ( int i=0;i<n;i++ ) {Smps ( i ) =zeros ( 0 );} 
    81         } 
     121                        void init()   ////////////// NOT FINISHED 
     122                        { 
     123                                Array<vec> Smps ( n ); 
     124                                //Gibbs sampling 
     125                                for ( int i=0;i<n;i++ ) {Smps ( i ) =zeros ( 0 );} 
     126                        } 
    82127//!Create a mixture density using known proposal 
    83         void merge ( const epdf* g0 ); 
     128                        void merge ( const epdf* g0 ); 
    84129//!Create a mixture density, make sure to call init() before the first call 
    85         void merge () {merge ( & ( Mix.posterior() ) );}; 
     130                        void merge () {merge ( & ( Mix.posterior() ) );}; 
    86131 
    87132//! Merge log-likelihood values 
    88         vec lognorm_merge ( mat &lW ); 
     133                        vec lognorm_merge ( mat &lW ); 
    89134//! sample from merged density 
    90135//! weight w is a 
    91         vec sample ( ) const { return Mix.posterior().sample();} 
    92         double evallog ( const vec &dt ) const { 
    93                 vec dtf=ones ( dt.length() +1 ); 
    94                 dtf.set_subvector ( 0,dt ); 
    95                 return Mix.logpred ( dtf ); 
    96         } 
    97         vec mean() const { 
    98                 const Vec<double> &w = eSmp._w(); 
    99                 const Array<vec> &S = eSmp._samples(); 
    100                 vec tmp=zeros ( dim);  
    101                 for ( int i=0; i<Ns; i++ ) { 
    102                         tmp+=w ( i ) *S ( i ); 
    103                 } 
    104                 return tmp; 
    105         } 
    106         mat covariance() const { 
    107                 const vec &w = eSmp._w(); 
    108                 const Array<vec> &S = eSmp._samples(); 
    109                  
    110                 vec mea = mean();  
    111                  
    112                 cout << sum(w) << "," << w*w <<endl; 
    113                  
    114                 mat Tmp=zeros(dim, dim); 
    115                 for ( int i=0; i<Ns; i++ ) { 
    116                         Tmp+=w ( i ) *outer_product(S ( i ), S(i)); 
    117                 } 
    118                 return Tmp-outer_product(mea,mea); 
    119         } 
    120         vec variance() const { 
    121                 const vec &w = eSmp._w(); 
    122                 const Array<vec> &S = eSmp._samples(); 
    123                  
    124                 vec tmp=zeros(dim); 
    125                 for ( int i=0; i<Ns; i++ ) { 
    126                         tmp+=w ( i ) *pow(S ( i ),2); 
    127                 } 
    128                 return tmp-pow(mean(),2); 
    129         } 
     136                        vec sample ( ) const { return Mix.posterior().sample();} 
     137                        double evallog ( const vec &dt ) const 
     138                        { 
     139                                vec dtf=ones ( dt.length() +1 ); 
     140                                dtf.set_subvector ( 0,dt ); 
     141                                return Mix.logpred ( dtf ); 
     142                        } 
     143                        vec mean() const 
     144                        { 
     145                                const Vec<double> &w = eSmp._w(); 
     146                                const Array<vec> &S = eSmp._samples(); 
     147                                vec tmp=zeros ( dim ); 
     148                                for ( int i=0; i<Ns; i++ ) 
     149                                { 
     150                                        tmp+=w ( i ) *S ( i ); 
     151                                } 
     152                                return tmp; 
     153                        } 
     154                        mat covariance() const 
     155                        { 
     156                                const vec &w = eSmp._w(); 
     157                                const Array<vec> &S = eSmp._samples(); 
     158 
     159                                vec mea = mean(); 
     160 
     161                                cout << sum ( w ) << "," << w*w <<endl; 
     162 
     163                                mat Tmp=zeros ( dim, dim ); 
     164                                for ( int i=0; i<Ns; i++ ) 
     165                                { 
     166                                        Tmp+=w ( i ) *outer_product ( S ( i ), S ( i ) ); 
     167                                } 
     168                                return Tmp-outer_product ( mea,mea ); 
     169                        } 
     170                        vec variance() const 
     171                        { 
     172                                const vec &w = eSmp._w(); 
     173                                const Array<vec> &S = eSmp._samples(); 
     174 
     175                                vec tmp=zeros ( dim ); 
     176                                for ( int i=0; i<Ns; i++ ) 
     177                                { 
     178                                        tmp+=w ( i ) *pow ( S ( i ),2 ); 
     179                                } 
     180                                return tmp-pow ( mean(),2 ); 
     181                        } 
    130182//! for future use 
    131         virtual ~merger() { 
    132                 for ( int i=0; i<n; i++ ) { 
    133                         delete dls ( i ); 
    134                         delete zdls ( i ); 
    135                 } 
     183                        virtual ~merger() 
     184                        { 
     185                                for ( int i=0; i<n; i++ ) 
     186                                { 
     187                                        delete dls ( i ); 
     188                                        delete zdls ( i ); 
     189                                } 
     190                                if ( DBG ) delete dbg; 
     191                        }; 
     192 
     193//! Access function 
     194                        MixEF& _Mix() {return Mix;} 
     195//! Access function 
     196                        eEmp& _Smp() {return eSmp;} 
    136197        }; 
    137198 
    138 //! Access function 
    139         MixEF& _Mix() {return Mix;} 
    140 //! Access function 
    141         eEmp& _Smp() {return eSmp;} 
    142 }; 
    143  
    144199} 
    145200 

  • bdm/stat/libEF.h

    r294 r299  
    139139                        vec mean() const {return mu;} 
    140140                        vec variance() const {return diag ( R.to_mat() );} 
    141 //      mlnorm<sq_T>* condition ( const RV &rvn ) const ; 
     141//      mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why? 
    142142                        mpdf* condition ( const RV &rvn ) const ; 
    143 //      enorm<sq_T>* marginal ( const RV &rv ) const; 
    144                         epdf* marginal ( const RV &rv ) const; 
     143        enorm<sq_T>* marginal ( const RV &rv ) const; 
     144//                      epdf* marginal ( const RV &rv ) const; 
    145145                        //!@} 
    146146 
    … …  
    10821082 
    10831083        template<class sq_T> 
    1084         epdf* enorm<sq_T>::marginal ( const RV &rvn ) const 
     1084        enorm<sq_T>* enorm<sq_T>::marginal ( const RV &rvn ) const 
    10851085        { 
    10861086                it_assert_debug ( isnamed(), "rv description is not assigned" );