Changeset 739

Show
Ignore:
Timestamp:
11/25/09 18:02:21 (14 years ago)
Author:
mido
Message:

the rest of h to cpp movements, with exception of from_setting and validate to avoid conflicts with Sarka

Location:
library/bdm
Files:
10 modified

Legend:

Unmodified
Added
Removed

  • library/bdm/estim/kalman.cpp

    r737 r739  
    150150} 
    151151 
     152void StateCanonical::connect_mlnorm ( const mlnorm<fsqmat> &ml ) { 
     153        //get ids of yrv 
     154        const RV &yrv = ml._rv(); 
     155        //need to determine u_t - it is all in _rvc that is not in ml._rv() 
     156        RV rgr0 = ml._rvc().remove_time(); 
     157        RV urv = rgr0.subt ( yrv ); 
     158 
     159        //We can do only 1d now... :( 
     160        bdm_assert ( yrv._dsize() == 1, "Only for SISO so far..." ); 
     161 
     162        // create names for 
     163        RV xrv; //empty 
     164        RV Crv; //empty 
     165        int td = ml._rvc().mint(); 
     166        // assuming strictly proper function!!! 
     167        for ( int t = -1; t >= td; t-- ) { 
     168                xrv.add ( yrv.copy_t ( t ) ); 
     169                Crv.add ( urv.copy_t ( t ) ); 
     170        } 
     171 
     172        // get mapp 
     173        th2A.set_connection ( xrv, ml._rvc() ); 
     174        th2C.set_connection ( Crv, ml._rvc() ); 
     175        th2D.set_connection ( urv, ml._rvc() ); 
     176 
     177        //set matrix sizes 
     178        this->A = zeros ( xrv._dsize(), xrv._dsize() ); 
     179        for ( int j = 1; j < xrv._dsize(); j++ ) { 
     180                A ( j, j - 1 ) = 1.0;    // off diagonal 
     181        } 
     182        this->B = zeros ( xrv._dsize(), 1 ); 
     183        this->B ( 0 ) = 1.0; 
     184        this->C = zeros ( 1, xrv._dsize() ); 
     185        this->D = zeros ( 1, urv._dsize() ); 
     186        this->Q = zeros ( xrv._dsize(), xrv._dsize() ); 
     187        // R is set by update 
     188 
     189        //set cache 
     190        this->A1row = zeros ( xrv._dsize() ); 
     191        this->C1row = zeros ( xrv._dsize() ); 
     192        this->D1row = zeros ( urv._dsize() ); 
     193 
     194        update_from ( ml ); 
     195        validate(); 
     196}; 
     197 
     198void StateCanonical::update_from ( const mlnorm<fsqmat> &ml ) { 
     199 
     200        vec theta = ml._A().get_row ( 0 ); // this 
     201 
     202        th2A.filldown ( theta, A1row ); 
     203        th2C.filldown ( theta, C1row ); 
     204        th2D.filldown ( theta, D1row ); 
     205 
     206        R = ml._R(); 
     207 
     208        A.set_row ( 0, A1row ); 
     209        C.set_row ( 0, C1row + D1row ( 0 ) *A1row ); 
     210        D.set_row ( 0, D1row ); 
     211} 
     212 
     213void StateFromARX::connect_mlnorm ( const mlnorm<chmat> &ml, RV &xrv, RV &urv ) { 
     214        //get ids of yrv 
     215        const RV &yrv = ml._rv(); 
     216        //need to determine u_t - it is all in _rvc that is not in ml._rv() 
     217        const RV &rgr = ml._rvc(); 
     218        RV rgr0 = rgr.remove_time(); 
     219        urv = rgr0.subt ( yrv ); 
     220 
     221        // create names for state variables 
     222        xrv = yrv; 
     223 
     224        int y_multiplicity = -rgr.mint ( yrv ); 
     225        int y_block_size = yrv.length() * ( y_multiplicity ); // current yt + all delayed yts 
     226        for ( int m = 0; m < y_multiplicity - 1; m++ ) { // ========= -1 is important see arx2statespace_notes 
     227                xrv.add ( yrv.copy_t ( -m - 1 ) ); //add delayed yt 
     228        } 
     229        //! temporary RV for connection to ml.rvc, since notation of xrv and ml.rvc does not match 
     230        RV xrv_ml = xrv.copy_t ( -1 ); 
     231 
     232        // add regressors 
     233        ivec u_block_sizes ( urv.length() ); // no_blocks = yt + unique rgr 
     234        for ( int r = 0; r < urv.length(); r++ ) { 
     235                RV R = urv.subselect ( vec_1 ( r ) ); //non-delayed element of rgr 
     236                int r_size = urv.size ( r ); 
     237                int r_multiplicity = -rgr.mint ( R ); 
     238                u_block_sizes ( r ) = r_size * r_multiplicity ; 
     239                for ( int m = 0; m < r_multiplicity; m++ ) { 
     240                        xrv.add ( R.copy_t ( -m - 1 ) ); //add delayed yt 
     241                        xrv_ml.add ( R.copy_t ( -m - 1 ) ); //add delayed yt 
     242                } 
     243        } 
     244        // add constant 
     245        if ( any ( ml._mu_const() != 0.0 ) ) { 
     246                have_constant = true; 
     247                xrv.add ( RV ( "bdm_reserved_constant_one", 1 ) ); 
     248        } else { 
     249                have_constant = false; 
     250        } 
     251 
     252 
     253        // get mapp 
     254        th2A.set_connection ( xrv_ml, ml._rvc() ); 
     255        th2B.set_connection ( urv, ml._rvc() ); 
     256 
     257        //set matrix sizes 
     258        this->A = zeros ( xrv._dsize(), xrv._dsize() ); 
     259        //create y block 
     260        diagonal_part ( this->A, yrv._dsize(), 0, y_block_size - yrv._dsize() ); 
     261 
     262        this->B = zeros ( xrv._dsize(), urv._dsize() ); 
     263        //add diagonals for rgr 
     264        int active_x = y_block_size; 
     265        for ( int r = 0; r < urv.length(); r++ ) { 
     266                diagonal_part ( this->A, active_x + urv.size ( r ), active_x, u_block_sizes ( r ) - urv.size ( r ) ); 
     267                this->B.set_submatrix ( active_x, 0, eye ( urv.size ( r ) ) ); 
     268                active_x += u_block_sizes ( r ); 
     269        } 
     270        this->C = zeros ( yrv._dsize(), xrv._dsize() ); 
     271        this->C.set_submatrix ( 0, 0, eye ( yrv._dsize() ) ); 
     272        this->D = zeros ( yrv._dsize(), urv._dsize() ); 
     273        this->R.setCh ( zeros ( yrv._dsize(), yrv._dsize() ) ); 
     274        this->Q.setCh ( zeros ( xrv._dsize(), xrv._dsize() ) ); 
     275        // Q is set by update 
     276 
     277        update_from ( ml ); 
     278        validate(); 
     279} 
     280         
     281void StateFromARX::update_from ( const mlnorm<chmat> &ml ) { 
     282        vec Arow = zeros ( A.cols() ); 
     283        vec Brow = zeros ( B.cols() ); 
     284        //  ROW- WISE EVALUATION ===== 
     285        for ( int i = 0; i < ml._rv()._dsize(); i++ ) { 
     286 
     287                vec theta = ml._A().get_row ( i ); 
     288 
     289                th2A.filldown ( theta, Arow ); 
     290                if ( have_constant ) { 
     291                        // constant is always at the end no need for datalink 
     292                        Arow ( A.cols() - 1 ) = ml._mu_const() ( i ); 
     293                } 
     294                this->A.set_row ( i, Arow ); 
     295 
     296                th2B.filldown ( theta, Brow ); 
     297                this->B.set_row ( i, Brow ); 
     298        } 
     299        this->Q._Ch().set_submatrix ( 0, 0, ml.__R()._Ch() ); 
     300 
     301} 
    152302 
    153303 

  • library/bdm/estim/kalman.h

    r737 r739  
    457457public: 
    458458        //! set up this object to match given mlnorm 
    459         void connect_mlnorm ( const mlnorm<fsqmat> &ml ) { 
    460                 //get ids of yrv 
    461                 const RV &yrv = ml._rv(); 
    462                 //need to determine u_t - it is all in _rvc that is not in ml._rv() 
    463                 RV rgr0 = ml._rvc().remove_time(); 
    464                 RV urv = rgr0.subt ( yrv ); 
    465  
    466                 //We can do only 1d now... :( 
    467                 bdm_assert ( yrv._dsize() == 1, "Only for SISO so far..." ); 
    468  
    469                 // create names for 
    470                 RV xrv; //empty 
    471                 RV Crv; //empty 
    472                 int td = ml._rvc().mint(); 
    473                 // assuming strictly proper function!!! 
    474                 for ( int t = -1; t >= td; t-- ) { 
    475                         xrv.add ( yrv.copy_t ( t ) ); 
    476                         Crv.add ( urv.copy_t ( t ) ); 
    477                 } 
    478  
    479                 // get mapp 
    480                 th2A.set_connection ( xrv, ml._rvc() ); 
    481                 th2C.set_connection ( Crv, ml._rvc() ); 
    482                 th2D.set_connection ( urv, ml._rvc() ); 
    483  
    484                 //set matrix sizes 
    485                 this->A = zeros ( xrv._dsize(), xrv._dsize() ); 
    486                 for ( int j = 1; j < xrv._dsize(); j++ ) { 
    487                         A ( j, j - 1 ) = 1.0;    // off diagonal 
    488                 } 
    489                 this->B = zeros ( xrv._dsize(), 1 ); 
    490                 this->B ( 0 ) = 1.0; 
    491                 this->C = zeros ( 1, xrv._dsize() ); 
    492                 this->D = zeros ( 1, urv._dsize() ); 
    493                 this->Q = zeros ( xrv._dsize(), xrv._dsize() ); 
    494                 // R is set by update 
    495  
    496                 //set cache 
    497                 this->A1row = zeros ( xrv._dsize() ); 
    498                 this->C1row = zeros ( xrv._dsize() ); 
    499                 this->D1row = zeros ( urv._dsize() ); 
    500  
    501                 update_from ( ml ); 
    502                 validate(); 
    503         }; 
     459        void connect_mlnorm ( const mlnorm<fsqmat> &ml ); 
     460 
    504461        //! fast function to update parameters from ml - not checked for compatibility!! 
    505         void update_from ( const mlnorm<fsqmat> &ml ) { 
    506  
    507                 vec theta = ml._A().get_row ( 0 ); // this 
    508  
    509                 th2A.filldown ( theta, A1row ); 
    510                 th2C.filldown ( theta, C1row ); 
    511                 th2D.filldown ( theta, D1row ); 
    512  
    513                 R = ml._R(); 
    514  
    515                 A.set_row ( 0, A1row ); 
    516                 C.set_row ( 0, C1row + D1row ( 0 ) *A1row ); 
    517                 D.set_row ( 0, D1row ); 
    518  
    519         } 
     462        void update_from ( const mlnorm<fsqmat> &ml ); 
    520463}; 
    521464/*! 
    … …  
    553496        //!While mlnorm typically assumes that \f$ u_t \rightarrow y_t \f$ in state space it is \f$ u_{t-1} \rightarrow y_t \f$ 
    554497        //! For consequences in notation of internal variable xt see arx2statespace_notes.lyx. 
    555         void connect_mlnorm ( const mlnorm<chmat> &ml, RV &xrv, RV &urv ) { 
    556  
    557                 //get ids of yrv 
    558                 const RV &yrv = ml._rv(); 
    559                 //need to determine u_t - it is all in _rvc that is not in ml._rv() 
    560                 const RV &rgr = ml._rvc(); 
    561                 RV rgr0 = rgr.remove_time(); 
    562                 urv = rgr0.subt ( yrv ); 
    563  
    564                 // create names for state variables 
    565                 xrv = yrv; 
    566  
    567                 int y_multiplicity = -rgr.mint ( yrv ); 
    568                 int y_block_size = yrv.length() * ( y_multiplicity ); // current yt + all delayed yts 
    569                 for ( int m = 0; m < y_multiplicity - 1; m++ ) { // ========= -1 is important see arx2statespace_notes 
    570                         xrv.add ( yrv.copy_t ( -m - 1 ) ); //add delayed yt 
    571                 } 
    572                 //! temporary RV for connection to ml.rvc, since notation of xrv and ml.rvc does not match 
    573                 RV xrv_ml = xrv.copy_t ( -1 ); 
    574  
    575                 // add regressors 
    576                 ivec u_block_sizes ( urv.length() ); // no_blocks = yt + unique rgr 
    577                 for ( int r = 0; r < urv.length(); r++ ) { 
    578                         RV R = urv.subselect ( vec_1 ( r ) ); //non-delayed element of rgr 
    579                         int r_size = urv.size ( r ); 
    580                         int r_multiplicity = -rgr.mint ( R ); 
    581                         u_block_sizes ( r ) = r_size * r_multiplicity ; 
    582                         for ( int m = 0; m < r_multiplicity; m++ ) { 
    583                                 xrv.add ( R.copy_t ( -m - 1 ) ); //add delayed yt 
    584                                 xrv_ml.add ( R.copy_t ( -m - 1 ) ); //add delayed yt 
    585                         } 
    586                 } 
    587                 // add constant 
    588                 if ( any ( ml._mu_const() != 0.0 ) ) { 
    589                         have_constant = true; 
    590                         xrv.add ( RV ( "bdm_reserved_constant_one", 1 ) ); 
    591                 } else { 
    592                         have_constant = false; 
    593                 } 
    594  
    595  
    596                 // get mapp 
    597                 th2A.set_connection ( xrv_ml, ml._rvc() ); 
    598                 th2B.set_connection ( urv, ml._rvc() ); 
    599  
    600                 //set matrix sizes 
    601                 this->A = zeros ( xrv._dsize(), xrv._dsize() ); 
    602                 //create y block 
    603                 diagonal_part ( this->A, yrv._dsize(), 0, y_block_size - yrv._dsize() ); 
    604  
    605                 this->B = zeros ( xrv._dsize(), urv._dsize() ); 
    606                 //add diagonals for rgr 
    607                 int active_x = y_block_size; 
    608                 for ( int r = 0; r < urv.length(); r++ ) { 
    609                         diagonal_part ( this->A, active_x + urv.size ( r ), active_x, u_block_sizes ( r ) - urv.size ( r ) ); 
    610                         this->B.set_submatrix ( active_x, 0, eye ( urv.size ( r ) ) ); 
    611                         active_x += u_block_sizes ( r ); 
    612                 } 
    613                 this->C = zeros ( yrv._dsize(), xrv._dsize() ); 
    614                 this->C.set_submatrix ( 0, 0, eye ( yrv._dsize() ) ); 
    615                 this->D = zeros ( yrv._dsize(), urv._dsize() ); 
    616                 this->R.setCh ( zeros ( yrv._dsize(), yrv._dsize() ) ); 
    617                 this->Q.setCh ( zeros ( xrv._dsize(), xrv._dsize() ) ); 
    618                 // Q is set by update 
    619  
    620                 update_from ( ml ); 
    621                 validate(); 
    622         }; 
     498        void connect_mlnorm ( const mlnorm<chmat> &ml, RV &xrv, RV &urv ); 
     499 
    623500        //! fast function to update parameters from ml - not checked for compatibility!! 
    624         void update_from ( const mlnorm<chmat> &ml ) { 
    625  
    626                 vec Arow = zeros ( A.cols() ); 
    627                 vec Brow = zeros ( B.cols() ); 
    628                 //  ROW- WISE EVALUATION ===== 
    629                 for ( int i = 0; i < ml._rv()._dsize(); i++ ) { 
    630  
    631                         vec theta = ml._A().get_row ( i ); 
    632  
    633                         th2A.filldown ( theta, Arow ); 
    634                         if ( have_constant ) { 
    635                                 // constant is always at the end no need for datalink 
    636                                 Arow ( A.cols() - 1 ) = ml._mu_const() ( i ); 
    637                         } 
    638                         this->A.set_row ( i, Arow ); 
    639  
    640                         th2B.filldown ( theta, Brow ); 
    641                         this->B.set_row ( i, Brow ); 
    642                 } 
    643                 this->Q._Ch().set_submatrix ( 0, 0, ml.__R()._Ch() ); 
    644  
    645         }; 
     501        void update_from ( const mlnorm<chmat> &ml ); 
     502 
    646503        //! access function 
    647504        bool _have_constant() const { 

  • library/bdm/estim/particles.cpp

    r737 r739  
    7575// } 
    7676 
     77vec MPF::mpfepdf::mean() const { 
     78        const vec &w = pf->posterior()._w(); 
     79        vec pom = zeros ( BMs ( 0 )->posterior ().dimension() ); 
     80        //compute mean of BMs 
     81        for ( int i = 0; i < w.length(); i++ ) { 
     82                pom += BMs ( i )->posterior().mean() * w ( i ); 
     83        } 
     84        return concat ( pf->posterior().mean(), pom ); 
     85} 
     86 
     87vec MPF::mpfepdf::variance() const { 
     88        const vec &w = pf->posterior()._w(); 
     89 
     90        vec pom = zeros ( BMs ( 0 )->posterior ().dimension() ); 
     91        vec pom2 = zeros ( BMs ( 0 )->posterior ().dimension() ); 
     92        vec mea; 
     93 
     94        for ( int i = 0; i < w.length(); i++ ) { 
     95                // save current mean 
     96                mea = BMs ( i )->posterior().mean(); 
     97                pom += mea * w ( i ); 
     98                //compute variance 
     99                pom2 += ( BMs ( i )->posterior().variance() + pow ( mea, 2 ) ) * w ( i ); 
     100        } 
     101        return concat ( pf->posterior().variance(), pom2 - pow ( pom, 2 ) ); 
     102} 
     103 
     104void MPF::mpfepdf::qbounds ( vec &lb, vec &ub, double perc ) const { 
     105        //bounds on particles 
     106        vec lbp; 
     107        vec ubp; 
     108        pf->posterior().qbounds ( lbp, ubp ); 
     109 
     110        //bounds on Components 
     111        int dimC = BMs ( 0 )->posterior().dimension(); 
     112        int j; 
     113        // temporary 
     114        vec lbc ( dimC ); 
     115        vec ubc ( dimC ); 
     116        // minima and maxima 
     117        vec Lbc ( dimC ); 
     118        vec Ubc ( dimC ); 
     119        Lbc = std::numeric_limits<double>::infinity(); 
     120        Ubc = -std::numeric_limits<double>::infinity(); 
     121 
     122        for ( int i = 0; i < BMs.length(); i++ ) { 
     123                // check Coms 
     124                BMs ( i )->posterior().qbounds ( lbc, ubc ); 
     125                //save either minima or maxima 
     126                for ( j = 0; j < dimC; j++ ) { 
     127                        if ( lbc ( j ) < Lbc ( j ) ) { 
     128                                Lbc ( j ) = lbc ( j ); 
     129                        } 
     130                        if ( ubc ( j ) > Ubc ( j ) ) { 
     131                                Ubc ( j ) = ubc ( j ); 
     132                        } 
     133                } 
     134        } 
     135        lb = concat ( lbp, Lbc ); 
     136        ub = concat ( ubp, Ubc ); 
     137} 
     138 
     139 
     140 
    77141void MPF::bayes ( const vec &yt, const vec &cond ) { 
    78142        // follows PF::bayes in most places!!! 

  • library/bdm/estim/particles.h

    r737 r739  
    236236                        bdm_assert_debug ( dim == rv._dsize(), "Wrong name " ); 
    237237                } 
    238                 vec mean() const { 
    239                         const vec &w = pf->posterior()._w(); 
    240                         vec pom = zeros ( BMs ( 0 )->posterior ().dimension() ); 
    241                         //compute mean of BMs 
    242                         for ( int i = 0; i < w.length(); i++ ) { 
    243                                 pom += BMs ( i )->posterior().mean() * w ( i ); 
    244                         } 
    245                         return concat ( pf->posterior().mean(), pom ); 
    246                 } 
    247                 vec variance() const { 
    248                         const vec &w = pf->posterior()._w(); 
    249  
    250                         vec pom = zeros ( BMs ( 0 )->posterior ().dimension() ); 
    251                         vec pom2 = zeros ( BMs ( 0 )->posterior ().dimension() ); 
    252                         vec mea; 
    253  
    254                         for ( int i = 0; i < w.length(); i++ ) { 
    255                                 // save current mean 
    256                                 mea = BMs ( i )->posterior().mean(); 
    257                                 pom += mea * w ( i ); 
    258                                 //compute variance 
    259                                 pom2 += ( BMs ( i )->posterior().variance() + pow ( mea, 2 ) ) * w ( i ); 
    260                         } 
    261                         return concat ( pf->posterior().variance(), pom2 - pow ( pom, 2 ) ); 
    262                 } 
    263  
    264                 void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const { 
    265                         //bounds on particles 
    266                         vec lbp; 
    267                         vec ubp; 
    268                         pf->posterior().qbounds ( lbp, ubp ); 
    269  
    270                         //bounds on Components 
    271                         int dimC = BMs ( 0 )->posterior().dimension(); 
    272                         int j; 
    273                         // temporary 
    274                         vec lbc ( dimC ); 
    275                         vec ubc ( dimC ); 
    276                         // minima and maxima 
    277                         vec Lbc ( dimC ); 
    278                         vec Ubc ( dimC ); 
    279                         Lbc = std::numeric_limits<double>::infinity(); 
    280                         Ubc = -std::numeric_limits<double>::infinity(); 
    281  
    282                         for ( int i = 0; i < BMs.length(); i++ ) { 
    283                                 // check Coms 
    284                                 BMs ( i )->posterior().qbounds ( lbc, ubc ); 
    285                                 //save either minima or maxima 
    286                                 for ( j = 0; j < dimC; j++ ) { 
    287                                         if ( lbc ( j ) < Lbc ( j ) ) { 
    288                                                 Lbc ( j ) = lbc ( j ); 
    289                                         } 
    290                                         if ( ubc ( j ) > Ubc ( j ) ) { 
    291                                                 Ubc ( j ) = ubc ( j ); 
    292                                         } 
    293                                 } 
    294                         } 
    295                         lb = concat ( lbp, Lbc ); 
    296                         ub = concat ( ubp, Ubc ); 
    297                 } 
     238                vec mean() const; 
     239 
     240                vec variance() const; 
     241 
     242                void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const; 
    298243 
    299244                vec sample() const { 

  • library/bdm/stat/emix.cpp

    r737 r739  
    3434 
    3535        return Coms ( i )->sample(); 
     36} 
     37 
     38vec emix::mean() const { 
     39        int i; 
     40        vec mu = zeros ( dim ); 
     41        for ( i = 0; i < w.length(); i++ ) { 
     42                mu += w ( i ) * Coms ( i )->mean(); 
     43        } 
     44        return mu; 
     45} 
     46 
     47vec emix::variance() const { 
     48        //non-central moment 
     49        vec mom2 = zeros ( dim ); 
     50        for ( int i = 0; i < w.length(); i++ ) { 
     51                mom2 += w ( i ) * ( Coms ( i )->variance() + pow ( Coms ( i )->mean(), 2 ) ); 
     52        } 
     53        //central moment 
     54        return mom2 - pow ( mean(), 2 ); 
     55} 
     56 
     57double emix::evallog ( const vec &val ) const { 
     58        int i; 
     59        double sum = 0.0; 
     60        for ( i = 0; i < w.length(); i++ ) { 
     61                sum += w ( i ) * exp ( Coms ( i )->evallog ( val ) ); 
     62        } 
     63        if ( sum == 0.0 ) { 
     64                sum = std::numeric_limits<double>::epsilon(); 
     65        } 
     66        double tmp = log ( sum ); 
     67        bdm_assert_debug ( std::isfinite ( tmp ), "Infinite" ); 
     68        return tmp; 
     69} 
     70 
     71vec emix::evallog_mat ( const mat &Val ) const { 
     72        vec x = zeros ( Val.cols() ); 
     73        for ( int i = 0; i < w.length(); i++ ) { 
     74                x += w ( i ) * exp ( Coms ( i )->evallog_mat ( Val ) ); 
     75        } 
     76        return log ( x ); 
     77}; 
     78 
     79mat emix::evallog_coms ( const mat &Val ) const { 
     80        mat X ( w.length(), Val.cols() ); 
     81        for ( int i = 0; i < w.length(); i++ ) { 
     82                X.set_row ( i, w ( i ) *exp ( Coms ( i )->evallog_mat ( Val ) ) ); 
     83        } 
     84        return X; 
    3685} 
    3786 
    … …  
    203252}; 
    204253 
     254double mprod::evallogcond ( const vec &val, const vec &cond ) { 
     255        int i; 
     256        double res = 0.0; 
     257        for ( i = pdfs.length() - 1; i >= 0; i-- ) { 
     258                /*                      if ( pdfs(i)->_rvc().count() >0) { 
     259                                                pdfs ( i )->condition ( dls ( i )->get_cond ( val,cond ) ); 
     260                                        } 
     261                                        // add logarithms 
     262                                        res += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );*/ 
     263                res += pdfs ( i )->evallogcond ( 
     264                           dls ( i )->pushdown ( val ), 
     265                           dls ( i )->get_cond ( val, cond ) 
     266                       ); 
     267        } 
     268        return res; 
     269} 
     270 
     271vec mprod::evallogcond_mat ( const mat &Dt, const vec &cond ) { 
     272        vec tmp ( Dt.cols() ); 
     273        for ( int i = 0; i < Dt.cols(); i++ ) { 
     274                tmp ( i ) = evallogcond ( Dt.get_col ( i ), cond ); 
     275        } 
     276        return tmp; 
     277} 
     278 
     279vec mprod::evallogcond_mat ( const Array<vec> &Dt, const vec &cond ) { 
     280        vec tmp ( Dt.length() ); 
     281        for ( int i = 0; i < Dt.length(); i++ ) { 
     282                tmp ( i ) = evallogcond ( Dt ( i ), cond ); 
     283        } 
     284        return tmp; 
     285} 
     286 
    205287void mprod::set_elements ( const Array<shared_ptr<pdf> > &mFacs ) { 
    206288        pdfs = mFacs; 
    … …  
    236318 
    237319        return Coms ( i )->samplecond ( cond ); 
     320} 
     321 
     322vec eprod::mean() const { 
     323        vec tmp ( dim ); 
     324        for ( int i = 0; i < epdfs.length(); i++ ) { 
     325                vec pom = epdfs ( i )->mean(); 
     326                dls ( i )->pushup ( tmp, pom ); 
     327        } 
     328        return tmp; 
     329} 
     330 
     331vec eprod::variance() const { 
     332        vec tmp ( dim ); //second moment 
     333        for ( int i = 0; i < epdfs.length(); i++ ) { 
     334                vec pom = epdfs ( i )->mean(); 
     335                dls ( i )->pushup ( tmp, pow ( pom, 2 ) ); 
     336        } 
     337        return tmp - pow ( mean(), 2 ); 
     338} 
     339vec eprod::sample() const { 
     340        vec tmp ( dim ); 
     341        for ( int i = 0; i < epdfs.length(); i++ ) { 
     342                vec pom = epdfs ( i )->sample(); 
     343                dls ( i )->pushup ( tmp, pom ); 
     344        } 
     345        return tmp; 
     346} 
     347double eprod::evallog ( const vec &val ) const { 
     348        double tmp = 0; 
     349        for ( int i = 0; i < epdfs.length(); i++ ) { 
     350                tmp += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) ); 
     351        } 
     352        bdm_assert_debug ( std::isfinite ( tmp ), "Infinite" ); 
     353        return tmp; 
    238354} 
    239355 
    … …  
    267383//              } 
    268384//      }; 
     385 

  • library/bdm/stat/emix.h

    r737 r739  
    129129 
    130130        vec sample() const; 
    131         vec mean() const { 
    132                 int i; 
    133                 vec mu = zeros ( dim ); 
    134                 for ( i = 0; i < w.length(); i++ ) { 
    135                         mu += w ( i ) * Coms ( i )->mean(); 
    136                 } 
    137                 return mu; 
    138         } 
    139         vec variance() const { 
    140                 //non-central moment 
    141                 vec mom2 = zeros ( dim ); 
    142                 for ( int i = 0; i < w.length(); i++ ) { 
    143                         mom2 += w ( i ) * ( Coms ( i )->variance() + pow ( Coms ( i )->mean(), 2 ) ); 
    144                 } 
    145                 //central moment 
    146                 return mom2 - pow ( mean(), 2 ); 
    147         } 
    148         double evallog ( const vec &val ) const { 
    149                 int i; 
    150                 double sum = 0.0; 
    151                 for ( i = 0; i < w.length(); i++ ) { 
    152                         sum += w ( i ) * exp ( Coms ( i )->evallog ( val ) ); 
    153                 } 
    154                 if ( sum == 0.0 ) { 
    155                         sum = std::numeric_limits<double>::epsilon(); 
    156                 } 
    157                 double tmp = log ( sum ); 
    158                 bdm_assert_debug ( std::isfinite ( tmp ), "Infinite" ); 
    159                 return tmp; 
    160         }; 
    161  
    162         vec evallog_mat ( const mat &Val ) const { 
    163                 vec x = zeros ( Val.cols() ); 
    164                 for ( int i = 0; i < w.length(); i++ ) { 
    165                         x += w ( i ) * exp ( Coms ( i )->evallog_mat ( Val ) ); 
    166                 } 
    167                 return log ( x ); 
    168         }; 
     131 
     132        vec mean() const; 
     133 
     134        vec variance() const; 
     135 
     136        double evallog ( const vec &val ) const; 
     137 
     138        vec evallog_mat ( const mat &Val ) const; 
    169139 
    170140        //! Auxiliary function that returns pdflog for each component 
    171         mat evallog_coms ( const mat &Val ) const { 
    172                 mat X ( w.length(), Val.cols() ); 
    173                 for ( int i = 0; i < w.length(); i++ ) { 
    174                         X.set_row ( i, w ( i ) *exp ( Coms ( i )->evallog_mat ( Val ) ) ); 
    175                 } 
    176                 return X; 
    177         }; 
     141        mat evallog_coms ( const mat &Val ) const; 
    178142 
    179143        shared_ptr<epdf> marginal ( const RV &rv ) const; 
    … …  
    334298        void set_elements ( const Array<shared_ptr<pdf> > &mFacs ); 
    335299 
    336         double evallogcond ( const vec &val, const vec &cond ) { 
    337                 int i; 
    338                 double res = 0.0; 
    339                 for ( i = pdfs.length() - 1; i >= 0; i-- ) { 
    340                         /*                      if ( pdfs(i)->_rvc().count() >0) { 
    341                                                         pdfs ( i )->condition ( dls ( i )->get_cond ( val,cond ) ); 
    342                                                 } 
    343                                                 // add logarithms 
    344                                                 res += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );*/ 
    345                         res += pdfs ( i )->evallogcond ( 
    346                                    dls ( i )->pushdown ( val ), 
    347                                    dls ( i )->get_cond ( val, cond ) 
    348                                ); 
    349                 } 
    350                 return res; 
    351         } 
    352         vec evallogcond_mat ( const mat &Dt, const vec &cond ) { 
    353                 vec tmp ( Dt.cols() ); 
    354                 for ( int i = 0; i < Dt.cols(); i++ ) { 
    355                         tmp ( i ) = evallogcond ( Dt.get_col ( i ), cond ); 
    356                 } 
    357                 return tmp; 
    358         }; 
    359         vec evallogcond_mat ( const Array<vec> &Dt, const vec &cond ) { 
    360                 vec tmp ( Dt.length() ); 
    361                 for ( int i = 0; i < Dt.length(); i++ ) { 
    362                         tmp ( i ) = evallogcond ( Dt ( i ), cond ); 
    363                 } 
    364                 return tmp; 
    365         }; 
    366  
     300        double evallogcond ( const vec &val, const vec &cond ); 
     301 
     302        vec evallogcond_mat ( const mat &Dt, const vec &cond ); 
     303 
     304        vec evallogcond_mat ( const Array<vec> &Dt, const vec &cond ); 
    367305 
    368306        //TODO smarter... 
    … …  
    441379        } 
    442380 
    443         vec mean() const { 
    444                 vec tmp ( dim ); 
    445                 for ( int i = 0; i < epdfs.length(); i++ ) { 
    446                         vec pom = epdfs ( i )->mean(); 
    447                         dls ( i )->pushup ( tmp, pom ); 
    448                 } 
    449                 return tmp; 
    450         } 
    451         vec variance() const { 
    452                 vec tmp ( dim ); //second moment 
    453                 for ( int i = 0; i < epdfs.length(); i++ ) { 
    454                         vec pom = epdfs ( i )->mean(); 
    455                         dls ( i )->pushup ( tmp, pow ( pom, 2 ) ); 
    456                 } 
    457                 return tmp - pow ( mean(), 2 ); 
    458         } 
    459         vec sample() const { 
    460                 vec tmp ( dim ); 
    461                 for ( int i = 0; i < epdfs.length(); i++ ) { 
    462                         vec pom = epdfs ( i )->sample(); 
    463                         dls ( i )->pushup ( tmp, pom ); 
    464                 } 
    465                 return tmp; 
    466         } 
    467         double evallog ( const vec &val ) const { 
    468                 double tmp = 0; 
    469                 for ( int i = 0; i < epdfs.length(); i++ ) { 
    470                         tmp += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) ); 
    471                 } 
    472                 bdm_assert_debug ( std::isfinite ( tmp ), "Infinite" ); 
    473                 return tmp; 
    474         } 
     381        vec mean() const; 
     382 
     383        vec variance() const; 
     384 
     385        vec sample() const; 
     386 
     387        double evallog ( const vec &val ) const; 
     388 
    475389        //!access function 
    476390        const epdf* operator () ( int i ) const { 

  • library/bdm/stat/exp_family.cpp

    r737 r739  
    279279        M = iLsub * Lpsi; 
    280280        R = ldR.to_mat()  ; 
     281} 
     282 
     283void egiw::log_register ( bdm::logger& L, const string& prefix ) { 
     284        if ( log_level == 3 ) { 
     285                root::log_register ( L, prefix ); 
     286                logrec->ids.set_length ( 2 ); 
     287                int th_dim = dimension() - dimx * ( dimx + 1 ) / 2; 
     288                logrec->ids ( 0 ) = L.add_vector ( RV ( "", th_dim ), prefix + logrec->L.prefix_sep() + "mean" ); 
     289                logrec->ids ( 1 ) = L.add_vector ( RV ( "", th_dim * th_dim ), prefix + logrec->L.prefix_sep() + "variance" ); 
     290        } else { 
     291                epdf::log_register ( L, prefix ); 
     292        } 
     293} 
     294 
     295void egiw::log_write() const { 
     296        if ( log_level == 3 ) { 
     297                mat M; 
     298                ldmat Lam; 
     299                ldmat Vz; 
     300                factorize ( M, Vz, Lam ); 
     301                logrec->L.log_vector ( logrec->ids ( 0 ), est_theta() ); 
     302                logrec->L.log_vector ( logrec->ids ( 1 ), cvectorize ( est_theta_cov().to_mat() ) ); 
     303        } else { 
     304                epdf::log_write(); 
     305        } 
     306 
     307} 
     308 
     309void multiBM::bayes ( const vec &yt, const vec &cond ) { 
     310        if ( frg < 1.0 ) { 
     311                beta *= frg; 
     312                last_lognc = est.lognc(); 
     313        } 
     314        beta += yt; 
     315        if ( evalll ) { 
     316                ll = est.lognc() - last_lognc; 
     317        } 
     318} 
     319 
     320double multiBM::logpred ( const vec &yt ) const { 
     321        eDirich pred ( est ); 
     322        vec &beta = pred._beta(); 
     323 
     324        double lll; 
     325        if ( frg < 1.0 ) { 
     326                beta *= frg; 
     327                lll = pred.lognc(); 
     328        } else if ( evalll ) { 
     329                lll = last_lognc; 
     330        } else { 
     331                lll = pred.lognc(); 
     332        } 
     333 
     334        beta += yt; 
     335        return pred.lognc() - lll; 
     336} 
     337void multiBM::flatten ( const BMEF* B ) { 
     338        const multiBM* E = dynamic_cast<const multiBM*> ( B ); 
     339        // sum(beta) should be equal to sum(B.beta) 
     340        const vec &Eb = E->beta;//const_cast<multiBM*> ( E )->_beta(); 
     341        beta *= ( sum ( Eb ) / sum ( beta ) ); 
     342        if ( evalll ) { 
     343                last_lognc = est.lognc(); 
     344        } 
    281345} 
    282346 
    … …  
    462526} 
    463527 
     528void mlstudent::condition ( const vec &cond ) { 
     529        if ( cond.length() > 0 ) { 
     530                iepdf._mu() = A * cond + mu_const; 
     531        } else { 
     532                iepdf._mu() =  mu_const; 
     533        } 
     534        double zeta; 
     535        //ugly hack! 
     536        if ( ( cond.length() + 1 ) == Lambda.rows() ) { 
     537                zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) ); 
     538        } else { 
     539                zeta = Lambda.invqform ( cond ); 
     540        } 
     541        _R = Re; 
     542        _R *= ( 1 + zeta );// / ( nu ); << nu is in Re!!!!!! 
     543} 
     544 
     545void eEmp::qbounds ( vec &lb, vec &ub, double perc ) const { 
     546        // lb in inf so than it will be pushed below; 
     547        lb.set_size ( dim ); 
     548        ub.set_size ( dim ); 
     549        lb = std::numeric_limits<double>::infinity(); 
     550        ub = -std::numeric_limits<double>::infinity(); 
     551        int j; 
     552        for ( int i = 0; i < n; i++ ) { 
     553                for ( j = 0; j < dim; j++ ) { 
     554                        if ( samples ( i ) ( j ) < lb ( j ) ) { 
     555                                lb ( j ) = samples ( i ) ( j ); 
     556                        } 
     557                        if ( samples ( i ) ( j ) > ub ( j ) ) { 
     558                                ub ( j ) = samples ( i ) ( j ); 
     559                        } 
     560                } 
     561        } 
     562} 
     563 
     564 
    464565}; 

  • library/bdm/stat/exp_family.h

    r737 r739  
    323323                // check sizes, rvs etc. 
    324324        } 
    325         void log_register ( bdm::logger& L, const string& prefix ) { 
    326                 if ( log_level == 3 ) { 
    327                         root::log_register ( L, prefix ); 
    328                         logrec->ids.set_length ( 2 ); 
    329                         int th_dim = dimension() - dimx * ( dimx + 1 ) / 2; 
    330                         logrec->ids ( 0 ) = L.add_vector ( RV ( "", th_dim ), prefix + logrec->L.prefix_sep() + "mean" ); 
    331                         logrec->ids ( 1 ) = L.add_vector ( RV ( "", th_dim * th_dim ), prefix + logrec->L.prefix_sep() + "variance" ); 
    332                 } else { 
    333                         epdf::log_register ( L, prefix ); 
    334                 } 
    335         } 
    336         void log_write() const { 
    337                 if ( log_level == 3 ) { 
    338                         mat M; 
    339                         ldmat Lam; 
    340                         ldmat Vz; 
    341                         factorize ( M, Vz, Lam ); 
    342                         logrec->L.log_vector ( logrec->ids ( 0 ), est_theta() ); 
    343                         logrec->L.log_vector ( logrec->ids ( 1 ), cvectorize ( est_theta_cov().to_mat() ) ); 
    344                 } else { 
    345                         epdf::log_write(); 
    346                 } 
    347  
    348         } 
     325        void log_register ( bdm::logger& L, const string& prefix ); 
     326 
     327        void log_write() const; 
    349328        //!@} 
    350329}; 
    … …  
    526505                beta = mB->beta; 
    527506        } 
    528         void bayes ( const vec &yt, const vec &cond = empty_vec ) { 
    529                 if ( frg < 1.0 ) { 
    530                         beta *= frg; 
    531                         last_lognc = est.lognc(); 
    532                 } 
    533                 beta += yt; 
    534                 if ( evalll ) { 
    535                         ll = est.lognc() - last_lognc; 
    536                 } 
    537         } 
    538         double logpred ( const vec &yt ) const { 
    539                 eDirich pred ( est ); 
    540                 vec &beta = pred._beta(); 
    541  
    542                 double lll; 
    543                 if ( frg < 1.0 ) { 
    544                         beta *= frg; 
    545                         lll = pred.lognc(); 
    546                 } else if ( evalll ) { 
    547                         lll = last_lognc; 
    548                 } else { 
    549                         lll = pred.lognc(); 
    550                 } 
    551  
    552                 beta += yt; 
    553                 return pred.lognc() - lll; 
    554         } 
    555         void flatten ( const BMEF* B ) { 
    556                 const multiBM* E = dynamic_cast<const multiBM*> ( B ); 
    557                 // sum(beta) should be equal to sum(B.beta) 
    558                 const vec &Eb = E->beta;//const_cast<multiBM*> ( E )->_beta(); 
    559                 beta *= ( sum ( Eb ) / sum ( beta ) ); 
    560                 if ( evalll ) { 
    561                         last_lognc = est.lognc(); 
    562                 } 
    563         } 
     507        void bayes ( const vec &yt, const vec &cond = empty_vec ); 
     508 
     509        double logpred ( const vec &yt ) const; 
     510 
     511        void flatten ( const BMEF* B ); 
     512 
    564513        //! return correctly typed posterior (covariant return) 
    565514        const eDirich& posterior() const { 
    … …  
    971920                Lambda = Lambda0; 
    972921        } 
    973         void condition ( const vec &cond ) { 
    974                 if ( cond.length() > 0 ) { 
    975                         iepdf._mu() = A * cond + mu_const; 
    976                 } else { 
    977                         iepdf._mu() =  mu_const; 
    978                 } 
    979                 double zeta; 
    980                 //ugly hack! 
    981                 if ( ( cond.length() + 1 ) == Lambda.rows() ) { 
    982                         zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) ); 
    983                 } else { 
    984                         zeta = Lambda.invqform ( cond ); 
    985                 } 
    986                 _R = Re; 
    987                 _R *= ( 1 + zeta );// / ( nu ); << nu is in Re!!!!!! 
    988         }; 
     922 
     923        void condition ( const vec &cond );  
    989924 
    990925        void validate() { 
    … …  
    15261461        } 
    15271462        //! For this class, qbounds are minimum and maximum value of the population! 
    1528         void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const { 
    1529                 // lb in inf so than it will be pushed below; 
    1530                 lb.set_size ( dim ); 
    1531                 ub.set_size ( dim ); 
    1532                 lb = std::numeric_limits<double>::infinity(); 
    1533                 ub = -std::numeric_limits<double>::infinity(); 
    1534                 int j; 
    1535                 for ( int i = 0; i < n; i++ ) { 
    1536                         for ( j = 0; j < dim; j++ ) { 
    1537                                 if ( samples ( i ) ( j ) < lb ( j ) ) { 
    1538                                         lb ( j ) = samples ( i ) ( j ); 
    1539                                 } 
    1540                                 if ( samples ( i ) ( j ) > ub ( j ) ) { 
    1541                                         ub ( j ) = samples ( i ) ( j ); 
    1542                                 } 
    1543                         } 
    1544                 } 
    1545         } 
     1463        void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const; 
    15461464}; 
    15471465 

  • library/bdm/stat/merger.cpp

    r737 r739  
    88                Npoints ( 0 ), DBG ( false ), dbg_file ( 0 ) { 
    99        set_sources ( S ); 
     10} 
     11 
     12void merger_base::set_sources ( const Array<shared_ptr<pdf> > &Sources ) { 
     13        pdfs = Sources; 
     14        Nsources = pdfs.length(); 
     15        //set sizes 
     16        dls.set_size ( Sources.length() ); 
     17        rvzs.set_size ( Sources.length() ); 
     18        zdls.set_size ( Sources.length() ); 
     19 
     20        rv = get_composite_rv ( pdfs, /* checkoverlap = */ false ); 
     21 
     22        RV rvc; 
     23        // Extend rv by rvc! 
     24        for ( int i = 0; i < pdfs.length(); i++ ) { 
     25                RV rvx = pdfs ( i )->_rvc().subt ( rv ); 
     26                rvc.add ( rvx ); // add rv to common rvc 
     27        } 
     28 
     29        // join rv and rvc - see descriprion 
     30        rv.add ( rvc ); 
     31        // get dimension 
     32        dim = rv._dsize(); 
     33 
     34        // create links between sources and common rv 
     35        RV xytmp; 
     36        for ( int i = 0; i < pdfs.length(); i++ ) { 
     37                //Establich connection between pdfs and merger 
     38                dls ( i ) = new datalink_m2e; 
     39                dls ( i )->set_connection ( pdfs ( i )->_rv(), pdfs ( i )->_rvc(), rv ); 
     40 
     41                // find out what is missing in each pdf 
     42                xytmp = pdfs ( i )->_rv(); 
     43                xytmp.add ( pdfs ( i )->_rvc() ); 
     44                // z_i = common_rv-xy 
     45                rvzs ( i ) = rv.subt ( xytmp ); 
     46                //establish connection between extension (z_i|x,y)s and common rv 
     47                zdls ( i ) = new datalink_m2e; 
     48                zdls ( i )->set_connection ( rvzs ( i ), xytmp, rv ) ; 
     49        }; 
     50} 
     51 
     52void merger_base::set_support ( rectangular_support &Sup ) { 
     53        Npoints = Sup.points(); 
     54        eSmp.set_parameters ( Npoints, false ); 
     55        Array<vec> &samples = eSmp._samples(); 
     56        eSmp._w() = ones ( Npoints ) / Npoints; //unifrom size of bins 
     57        //set samples 
     58        samples ( 0 ) = Sup.first_vec(); 
     59        for ( int j = 1; j < Npoints; j++ ) { 
     60                samples ( j ) = Sup.next_vec(); 
     61        } 
     62} 
     63 
     64void merger_base::merge () { 
     65        validate(); 
     66 
     67        //check if sources overlap: 
     68        bool OK = true; 
     69        for ( int i = 0; i < pdfs.length(); i++ ) { 
     70                OK &= ( rvzs ( i )._dsize() == 0 ); // z_i is empty 
     71                OK &= ( pdfs ( i )->_rvc()._dsize() == 0 ); // y_i is empty 
     72        } 
     73 
     74        if ( OK ) { 
     75                mat lW = zeros ( pdfs.length(), eSmp._w().length() ); 
     76 
     77                vec emptyvec ( 0 ); 
     78                for ( int i = 0; i < pdfs.length(); i++ ) { 
     79                        for ( int j = 0; j < eSmp._w().length(); j++ ) { 
     80                                lW ( i, j ) = pdfs ( i )->evallogcond ( eSmp._samples() ( j ), emptyvec ); 
     81                        } 
     82                } 
     83 
     84                vec w_nn = merge_points ( lW ); 
     85                vec wtmp = exp ( w_nn - max ( w_nn ) ); 
     86                //renormalize 
     87                eSmp._w() = wtmp / sum ( wtmp ); 
     88        } else { 
     89                bdm_error ( "Sources are not compatible - use merger_mix" ); 
     90        } 
    1091} 
    1192 
    … …  
    62143} 
    63144 
     145vec merger_base::mean() const { 
     146        const Vec<double> &w = eSmp._w(); 
     147        const Array<vec> &S = eSmp._samples(); 
     148        vec tmp = zeros ( dim ); 
     149        for ( int i = 0; i < Npoints; i++ ) { 
     150                tmp += w ( i ) * S ( i ); 
     151        } 
     152        return tmp; 
     153} 
     154 
     155mat merger_base::covariance() const { 
     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 < Npoints; i++ ) { 
     165                Tmp += w ( i ) * outer_product ( S ( i ), S ( i ) ); 
     166        } 
     167        return Tmp - outer_product ( mea, mea ); 
     168} 
     169 
     170vec merger_base::variance() const { 
     171        const vec &w = eSmp._w(); 
     172        const Array<vec> &S = eSmp._samples(); 
     173 
     174        vec tmp = zeros ( dim ); 
     175        for ( int i = 0; i < Nsources; i++ ) { 
     176                tmp += w ( i ) * pow ( S ( i ), 2 ); 
     177        } 
     178        return tmp - pow ( mean(), 2 ); 
     179} 
     180 
    64181void merger_mix::merge ( ) { 
    65182        Array<vec> &Smp = eSmp._samples(); //aux 

  • library/bdm/stat/merger.h

    r737 r739  
    9393 
    9494        //! Function setting the main internal structures 
    95         void set_sources ( const Array<shared_ptr<pdf> > &Sources ) { 
    96                 pdfs = Sources; 
    97                 Nsources = pdfs.length(); 
    98                 //set sizes 
    99                 dls.set_size ( Sources.length() ); 
    100                 rvzs.set_size ( Sources.length() ); 
    101                 zdls.set_size ( Sources.length() ); 
    102  
    103                 rv = get_composite_rv ( pdfs, /* checkoverlap = */ false ); 
    104  
    105                 RV rvc; 
    106                 // Extend rv by rvc! 
    107                 for ( int i = 0; i < pdfs.length(); i++ ) { 
    108                         RV rvx = pdfs ( i )->_rvc().subt ( rv ); 
    109                         rvc.add ( rvx ); // add rv to common rvc 
    110                 } 
    111  
    112                 // join rv and rvc - see descriprion 
    113                 rv.add ( rvc ); 
    114                 // get dimension 
    115                 dim = rv._dsize(); 
    116  
    117                 // create links between sources and common rv 
    118                 RV xytmp; 
    119                 for ( int i = 0; i < pdfs.length(); i++ ) { 
    120                         //Establich connection between pdfs and merger 
    121                         dls ( i ) = new datalink_m2e; 
    122                         dls ( i )->set_connection ( pdfs ( i )->_rv(), pdfs ( i )->_rvc(), rv ); 
    123  
    124                         // find out what is missing in each pdf 
    125                         xytmp = pdfs ( i )->_rv(); 
    126                         xytmp.add ( pdfs ( i )->_rvc() ); 
    127                         // z_i = common_rv-xy 
    128                         rvzs ( i ) = rv.subt ( xytmp ); 
    129                         //establish connection between extension (z_i|x,y)s and common rv 
    130                         zdls ( i ) = new datalink_m2e; 
    131                         zdls ( i )->set_connection ( rvzs ( i ), xytmp, rv ) ; 
    132                 }; 
    133         } 
     95        void set_sources ( const Array<shared_ptr<pdf> > &Sources ); 
     96 
    13497        //! Set support points from rectangular grid 
    135         void set_support ( rectangular_support &Sup ) { 
    136                 Npoints = Sup.points(); 
    137                 eSmp.set_parameters ( Npoints, false ); 
    138                 Array<vec> &samples = eSmp._samples(); 
    139                 eSmp._w() = ones ( Npoints ) / Npoints; //unifrom size of bins 
    140                 //set samples 
    141                 samples ( 0 ) = Sup.first_vec(); 
    142                 for ( int j = 1; j < Npoints; j++ ) { 
    143                         samples ( j ) = Sup.next_vec(); 
    144                 } 
    145         } 
     98        void set_support ( rectangular_support &Sup ); 
     99 
    146100        //! Set support points from dicrete grid 
    147101        void set_support ( discrete_support &Sup ) { 
    … …  
    180134 
    181135        //!Merge given sources in given points 
    182         virtual void merge () { 
    183                 validate(); 
    184  
    185                 //check if sources overlap: 
    186                 bool OK = true; 
    187                 for ( int i = 0; i < pdfs.length(); i++ ) { 
    188                         OK &= ( rvzs ( i )._dsize() == 0 ); // z_i is empty 
    189                         OK &= ( pdfs ( i )->_rvc()._dsize() == 0 ); // y_i is empty 
    190                 } 
    191  
    192                 if ( OK ) { 
    193                         mat lW = zeros ( pdfs.length(), eSmp._w().length() ); 
    194  
    195                         vec emptyvec ( 0 ); 
    196                         for ( int i = 0; i < pdfs.length(); i++ ) { 
    197                                 for ( int j = 0; j < eSmp._w().length(); j++ ) { 
    198                                         lW ( i, j ) = pdfs ( i )->evallogcond ( eSmp._samples() ( j ), emptyvec ); 
    199                                 } 
    200                         } 
    201  
    202                         vec w_nn = merge_points ( lW ); 
    203                         vec wtmp = exp ( w_nn - max ( w_nn ) ); 
    204                         //renormalize 
    205                         eSmp._w() = wtmp / sum ( wtmp ); 
    206                 } else { 
    207                         bdm_error ( "Sources are not compatible - use merger_mix" ); 
    208                 } 
    209         }; 
    210  
     136        virtual void merge (); 
    211137 
    212138        //! Merge log-likelihood values in points using method specified by parameter METHOD 
    … …  
    216142        //! sample from merged density 
    217143//! weight w is a 
    218         vec mean() const { 
    219                 const Vec<double> &w = eSmp._w(); 
    220                 const Array<vec> &S = eSmp._samples(); 
    221                 vec tmp = zeros ( dim ); 
    222                 for ( int i = 0; i < Npoints; i++ ) { 
    223                         tmp += w ( i ) * S ( i ); 
    224                 } 
    225                 return tmp; 
    226         } 
    227         mat covariance() const { 
    228                 const vec &w = eSmp._w(); 
    229                 const Array<vec> &S = eSmp._samples(); 
    230  
    231                 vec mea = mean(); 
    232  
    233 //                      cout << sum (w) << "," << w*w << endl; 
    234  
    235                 mat Tmp = zeros ( dim, dim ); 
    236                 for ( int i = 0; i < Npoints; i++ ) { 
    237                         Tmp += w ( i ) * outer_product ( S ( i ), S ( i ) ); 
    238                 } 
    239                 return Tmp - outer_product ( mea, mea ); 
    240         } 
    241         vec variance() const { 
    242                 const vec &w = eSmp._w(); 
    243                 const Array<vec> &S = eSmp._samples(); 
    244  
    245                 vec tmp = zeros ( dim ); 
    246                 for ( int i = 0; i < Nsources; i++ ) { 
    247                         tmp += w ( i ) * pow ( S ( i ), 2 ); 
    248                 } 
    249                 return tmp - pow ( mean(), 2 ); 
    250         } 
     144        vec mean() const; 
     145 
     146        mat covariance() const; 
     147 
     148        vec variance() const; 
     149 
    251150        //!@} 
    252151