Changeset 198

Show
Ignore:
Timestamp:
11/04/08 14:54:34 (16 years ago)
Author:
smidl
Message:

opravy + zavedeni studenta + zakomentovani debug v mergeru

Files:
17 modified

Legend:

Unmodified
Added
Removed

  • bdm/estim/arx.cpp

    r182 r198  
    3333                lll = pred.lognc(); 
    3434        } 
    35         else  
    36                 if(evalll){lll=last_lognc;}else{lll=pred.lognc();} 
     35        else 
     36                if ( evalll ) {lll=last_lognc;} 
     37                else{lll=pred.lognc();} 
    3738 
    3839        V.opupdt ( dt,1.0 ); 
    … …  
    4243} 
    4344 
    44 ARX* ARX::_copy_ ( bool changerv) { 
     45ARX* ARX::_copy_ ( bool changerv ) { 
    4546        ARX* Tmp=new ARX ( *this ); 
    46         if ( changerv ) {Tmp->rv.newids(); Tmp->est._renewrv(Tmp->rv);} 
     47        if ( changerv ) {Tmp->rv.newids(); Tmp->est._renewrv ( Tmp->rv );} 
    4748        return Tmp; 
    4849} 
    … …  
    5556} 
    5657 
    57 enorm<ldmat>* ARX::predictor(const RV &rv){ 
    58         mat mu(rv.count(), 1); 
    59         mat R(rv.count(),rv.count()); 
    60         enorm<ldmat>* tmp;  
    61         tmp=new enorm<ldmat>(rv); 
     58enorm<ldmat>* ARX::predictor ( const RV &rv, const vec &rgr ) const { 
     59        mat mu ( rv.count(), V.rows()-rv.count() ); 
     60        mat R ( rv.count(),rv.count() ); 
     61        enorm<ldmat>* tmp; 
     62        tmp=new enorm<ldmat> ( rv ); 
     63 
     64        est.mean_mat ( mu,R ); //mu = 
     65        //correction for student-t  -- TODO check if correct!! 
     66        //R*=nu/(nu-2); 
     67        mat p_mu=mu.T() *rgr;   //the result is one column 
     68        tmp->set_parameters ( p_mu.get_col ( 0 ),ldmat ( R ) ); 
     69        return tmp; 
     70} 
     71 
     72mlnorm<ldmat>* ARX::predictor ( const RV &rv, const RV &rvc ) const { 
     73        int dif=V.rows() - rv.count() - rvc.count(); 
     74        it_assert_debug ( ( dif==0 ) || ( dif==1 ), "Give RVs do not match" ); 
     75 
     76        mat mu ( rv.count(), V.rows()-rv.count() ); 
     77        mat R ( rv.count(),rv.count() ); 
     78        mlnorm<ldmat>* tmp; 
     79        tmp=new mlnorm<ldmat> ( rv,rvc ); 
     80 
     81        est.mean_mat ( mu,R ); //mu = 
     82        mu = mu.T(); 
     83        //correction for student-t  -- TODO check if correct!! 
     84 
     85        if ( dif==0 ) { // no constant term 
     86                tmp->set_parameters ( mu, zeros ( rv.count() ), ldmat ( R ) ); 
     87        } 
     88        else { 
     89                //Assume the constant term is the last one: 
     90                tmp->set_parameters ( mu.get_cols (0,mu.cols()-2 ), mu.get_col ( mu.cols()-1 ), ldmat ( R ) ); 
     91        } 
     92        return tmp; 
     93} 
     94 
     95mlstudent* ARX::predictor_student ( const RV &rv, const RV &rvc ) const { 
     96        int dif=V.rows() - rv.count() - rvc.count(); 
     97        it_assert_debug ( ( dif==0 ) || ( dif==1 ), "Give RVs do not match" ); 
     98 
     99        mat mu ( rv.count(), V.rows()-rv.count() ); 
     100        mat R ( rv.count(),rv.count() ); 
     101        mlstudent* tmp; 
     102        tmp=new mlstudent ( rv,rvc ); 
     103 
     104        est.mean_mat ( mu,R ); // 
     105        mu = mu.T(); 
    62106         
    63         est.mean_mat(mu,R); 
    64         //here I know that mu can have only one column 
    65         tmp->set_parameters(mu.get_col(0),ldmat(R)); 
     107        int xdim = rv.count(); 
     108        int end = V._L().rows()-1; 
     109        ldmat Lam ( V._L() ( xdim,end,xdim,end ), V._D() ( xdim,end ) );  //exp val of R 
     110 
     111 
     112        if ( dif==0 ) { // no constant term 
     113                tmp->set_parameters ( mu, zeros ( rv.count() ), ldmat ( R ), Lam); 
     114        } 
     115        else { 
     116                //Assume the constant term is the last one: 
     117                tmp->set_parameters ( mu.get_cols (0,mu.cols()-2 ), mu.get_col ( mu.cols()-1 ), ldmat ( R ), Lam); 
     118        } 
    66119        return tmp; 
    67120} 

  • bdm/estim/arx.h

    r180 r198  
    7272                if(evalll){last_lognc=est.lognc();} 
    7373        } 
    74  
    75         enorm<ldmat>* predictor(const RV &rv);  
     74        //! Conditional version of the predictor 
     75        enorm<ldmat>* predictor(const RV &rv0, const vec &rgr) const;  
     76        enorm<ldmat>* predictor(const RV &rv0) const {it_assert_debug(rv0.count()==V.rows()-1,"Regressor is not only 1");return predictor(rv0,vec_1(1.0));} 
     77        mlnorm<ldmat>* predictor(const RV &rv0, const RV &rvc0) const; 
     78        mlstudent* predictor_student(const RV &rv0, const RV &rvc0) const; 
    7679        //! Brute force structure estimation.\return indeces of accepted regressors. 
    7780        ivec structure_est ( egiw Eg0 ); 

  • bdm/estim/merger.cpp

    r197 r198  
    88        vec lam = sum ( pow ( lW,2 ) )-nu*pow ( mu,2 ); 
    99        double coef=0.0; 
     10        vec sq2bl=sqrt ( 2*beta*lam ); //this term is everywhere 
    1011        switch ( nu ) { 
    1112                case 2: 
    12                         coef=sqrt ( beta*2 ) * ( 1-0.5*sqrt ( ( 4*beta-3 ) /beta ) ); 
    13                         return exp ( coef*sqrt ( lam ) + mu ); 
     13                        coef= ( 1-0.5*sqrt ( ( 4*beta-3 ) /beta ) ); 
     14                        return exp ( coef*sq2bl + mu ); 
    1415                        break; 
    15                 case 3://Ration of Bessel 
     16                case 3://Ratio of Bessel 
     17                        coef = sqrt ( ( 3*beta-2 ) /3*beta ); 
     18                        return elem_mult ( 
     19                                   elem_div ( besselk ( 0,sq2bl*coef ), besselk ( 0,sq2bl ) ), 
     20                                   exp ( mu ) ); 
    1621                        break; 
    1722                case 4: 
    … …  
    2429 
    2530void merger::merge ( const epdf* g0 ) { 
    26         it_file dbg ( "merger_debug.it" ); 
     31//      it_file dbg ( "merger_debug.it" ); 
    2732 
    2833        it_assert_debug ( rv.equal ( g0->_rv() ),"Incompatible g0" ); 
    … …  
    3641        for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );} 
    3742 
    38         dbg << Name ( "Smp_0" ) << Smp_ex; 
     43//      dbg << Name ( "Smp_0" ) << Smp_ex; 
    3944 
    4045        // Stuff for merging 
    … …  
    5156        Mix.init ( &A0, Smp_ex, Nc ); 
    5257        //Preserve initial mixture for repetitive estimation via flattening 
    53         MixEF Mix_init(Mix); 
     58        MixEF Mix_init ( Mix ); 
    5459 
    5560        // ============= MAIN LOOP ================== 
    … …  
    6368                //Re-estimate Mix 
    6469                //Re-Initialize Mixture model 
    65                 Mix.flatten(&Mix_init); 
     70                Mix.flatten ( &Mix_init ); 
    6671                Mix.bayesB ( Smp_ex, w*Ns ); 
    6772                Mpred = Mix.predictor ( rv ); // Allocation => must be deleted at the end!! 
    6873 
    69                 sprintf ( str,"Mpred_mean%d",niter ); 
    70                 dbg << Name ( str ) << Mpred->mean(); 
     74//              sprintf ( str,"Mpred_mean%d",niter ); 
     75//              dbg << Name ( str ) << Mpred->mean(); 
    7176 
    72                 // Generate new samples 
    73                 eSmp.set_samples ( Mpred ); 
    74                 for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );} 
     77                if ( niter<10 ) { 
     78                        // Generate new samples 
     79                        eSmp.set_samples ( Mpred ); 
     80                        for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );} 
    7581 
    76                 sprintf ( str,"Mpdf%d",niter ); 
    77                 for ( int i=0;i<Ns;i++ ) {Mix_pdf ( i ) = Mix.logpred ( Smp_ex.get_col ( i ) );} 
    78                 dbg << Name ( str ) << Mix_pdf; 
     82                } 
     83//              sprintf ( str,"Mpdf%d",niter ); 
     84//              for ( int i=0;i<Ns;i++ ) {Mix_pdf ( i ) = Mix.logpred ( Smp_ex.get_col ( i ) );} 
     85//              dbg << Name ( str ) << Mix_pdf; 
    7986 
    80                 sprintf ( str,"Smp%d",niter ); 
    81                 dbg << Name ( str ) << Smp_ex; 
     87//              sprintf ( str,"Smp%d",niter ); 
     88//              dbg << Name ( str ) << Smp_ex; 
    8289 
    8390                //Importace weighting 
    … …  
    99106                                        //compute vector of lw_src 
    100107                                        for ( int k=0;k<Ns;k++ ) { 
    101                                                 lw_src ( k ) += tmp_marg->evalpdflog ( dls ( i )->get_val ( Smp ( i ) ) ); 
     108                                                // Here val of tmp_marg = cond of mpdfs(i) ==> calling dls->get_cond 
     109                                                lw_src ( k ) += tmp_marg->evalpdflog ( dls ( i )->get_cond ( Smp ( k ) ) ); 
    102110                                        } 
    103111                                        delete tmp_marg; 
     112 
     113//                                      sprintf ( str,"marg%d",niter ); 
     114//                                      dbg << Name ( str ) << lw_src; 
     115 
    104116                                } 
    105117                                // Compute likelihood of the missing variable 
    … …  
    125137 
    126138                        } 
     139//                      it_assert_debug(std::isfinite(sum(lw_src)),"bad"); 
    127140                        lW.set_row ( i, lw_src ); // do not divide by mix 
    128141                } 
    … …  
    131144                        lw_mix ( j ) =Mix.logpred ( Smp_ex.get_col ( j ) ); 
    132145                } 
    133                 sprintf ( str,"lW%d",niter ); 
    134                 dbg << Name ( str ) << lW; 
     146//              sprintf ( str,"lW%d",niter ); 
     147//              dbg << Name ( str ) << lW; 
    135148 
    136149                w = lognorm_merge ( lW ); //merge 
    137150 
    138                 sprintf ( str,"w%d",niter ); 
    139                 dbg << Name ( str ) << w; 
    140                 sprintf ( str,"lw_m%d",niter ); 
    141                 dbg << Name ( str ) << lw_mix; 
     151//              sprintf ( str,"w%d",niter ); 
     152//              dbg << Name ( str ) << w; 
     153//              sprintf ( str,"lw_m%d",niter ); 
     154//              dbg << Name ( str ) << lw_mix; 
    142155 
    143156                //Importance weighting 
    … …  
    146159                w /=sum ( w ); 
    147160 
    148                 sprintf ( str,"w_is_%d",niter ); 
    149                 dbg << Name ( str ) << w; 
     161//              sprintf ( str,"w_is_%d",niter ); 
     162//              dbg << Name ( str ) << w; 
    150163 
    151164//              eSmp.resample(); // So that it can be used in bayes 
    152165//              for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );} 
    153166 
    154                 sprintf ( str,"Smp_res%d",niter ); 
    155                 dbg << Name ( str ) << Smp; 
     167//              sprintf ( str,"Smp_res%d",niter ); 
     168//              dbg << Name ( str ) << Smp; 
    156169 
    157170                // ==== stopping rule === 
    158171                niter++; 
    159                 converged = ( niter>9 ); 
     172                converged = ( niter>5 ); 
    160173        } 
    161174 

  • bdm/estim/merger.h

    r192 r198  
    3535        Array<RV> rvzs; 
    3636        //! Data Links of rv0 mpdfs - these will be conditioned the (rv,rvc) of mpdfs 
    37         Array<datalink_m2e*> zdls;  
    38          
     37        Array<datalink_m2e*> zdls; 
     38 
    3939        //!Number of samples used in approximation 
    4040        int Ns; 
    … …  
    4747        merger ( const Array<mpdf*> &S ) : 
    4848                        compositepdf ( S ), epdf ( getrv ( false ) ), 
    49                         Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ), dls(n), rvzs(n), zdls(n) { 
    50                                 RV ztmp; 
     49                        Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ), dls ( n ), rvzs ( n ), zdls ( n ) { 
     50                RV ztmp; 
     51                // Extend rv by rvc! 
     52                RV rvc; setrvc ( rv,rvc ); 
     53                rv.add ( rvc ); 
    5154                for ( int i=0;i<n;i++ ) { 
    5255                        //Establich connection between mpdfs and merger 
    53                         dls ( i ) = new datalink_m2e ( mpdfs ( i )->_rv(), mpdfs(i)->_rvc(), rv ); 
     56                        dls ( i ) = new datalink_m2e ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), rv ); 
    5457                        // find out what is missing in each mpdf 
    55                         ztmp= mpdfs(i)->_rv(); 
    56                         ztmp.add(mpdfs(i)->_rvc()); 
    57                         rvzs(i)=rv.subt(ztmp); 
    58                         zdls(i) = new datalink_m2e (rvzs(i), ztmp, rv ) ; 
     58                        ztmp= mpdfs ( i )->_rv(); 
     59                        ztmp.add ( mpdfs ( i )->_rvc() ); 
     60                        rvzs ( i ) =rv.subt ( ztmp ); 
     61                        zdls ( i ) = new datalink_m2e ( rvzs ( i ), ztmp, rv ) ; 
    5962                }; 
    6063                //Set Default values of parameters 
    … …  
    9093//! for future use 
    9194        virtual ~merger() { 
    92                 for (int i=0; i<n; i++){ 
    93                         delete dls(i); 
    94                         delete zdls(i);                  
     95                for ( int i=0; i<n; i++ ) { 
     96                        delete dls ( i ); 
     97                        delete zdls ( i ); 
    9598                } 
    9699        }; 

  • bdm/estim/mixef.cpp

    r197 r198  
    132132} 
    133133 
    134 emix* MixEF::predictor ( const RV &rv ) { 
     134emix* MixEF::predictor ( const RV &rv ) const { 
    135135        Array<epdf*> pC ( n ); 
    136136        for ( int i=0;i<n;i++ ) {pC ( i ) =Coms ( i )->predictor ( rv );} 
    … …  
    149149                Coms(i)->flatten(Mix2->Coms(i)); 
    150150        } 
    151         //Faltten weights 
    152         weights.flatten(&(Mix2->weights)); 
     151        //Flatten weights = make them equal!! 
     152        weights.set_statistics(&(Mix2->weights)); 
    153153} 

  • bdm/estim/mixef.h

    r197 r198  
    105105        double logpred ( const vec &dt ) const; 
    106106        const epdf& _epdf() const {return *est;} 
    107         emix* predictor(const RV &rv); 
     107        emix* predictor(const RV &rv) const; 
    108108        //! Flatten the density as if it was not estimated from the data 
    109109        void flatten(const BMEF* M2); 

  • bdm/stat/libBM.h

    r193 r198  
    426426 
    427427        //!Constructs a predictive density (marginal density on data) 
    428         virtual epdf* predictor ( const RV &rv ) {it_error ( "Not implemented" );return NULL;}; 
     428        virtual epdf* predictor ( const RV &rv ) const {it_error ( "Not implemented" );return NULL;}; 
    429429 
    430430        //! Destructor for future use; 

  • bdm/stat/libEF.cpp

    r197 r198  
    8282                mat R; 
    8383                mean_mat(M,R); 
    84                 return cvectorize (concat_vertical(M.T(),R)); 
     84                return cvectorize (concat_vertical(M,R)); 
    8585        } 
    8686 
    … …  
    9595         
    9696        // set mean value 
    97         M= L ( xdim,end,0,xdim-1 ).T()*iLsub; 
     97        mat Lpsi = L ( xdim,end,0,xdim-1 ); 
     98        M= iLsub*Lpsi; 
    9899        R= ldR.to_mat()  ; 
    99100} 

  • bdm/stat/libEF.h

    r197 r198  
    9494        //!Flatten the posterior as if to keep nu0 data 
    9595//      virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );} 
    96          
     96 
    9797        BMEF* _copy_ ( bool changerv=false ) {it_error ( "function _copy_ not implemented for this BM" ); return NULL;}; 
    9898}; 
    … …  
    192192        //! returns a pointer to the internal statistics. Use with Care! 
    193193        double& _nu() {return nu;} 
    194         void pow ( double p ){V*=p;nu*=p;}; 
     194        void pow ( double p ) {V*=p;nu*=p;}; 
    195195}; 
    196196 
    … …  
    239239        //! Conjugate prior and posterior 
    240240        eDirich est; 
    241         //! Pointer inside est to sufficient statistics  
     241        //! Pointer inside est to sufficient statistics 
    242242        vec &beta; 
    243243public: 
    … …  
    271271                // sum(beta) should be equal to sum(B.beta) 
    272272                const vec &Eb=E->beta;//const_cast<multiBM*> ( E )->_beta(); 
    273                 beta*= ( sum ( Eb) /sum ( beta ) ); 
     273                beta*= ( sum ( Eb ) /sum ( beta ) ); 
    274274                if ( evalll ) {last_lognc=est.lognc();} 
    275275        } 
    … …  
    372372template<class sq_T> 
    373373class mlnorm : public mEF { 
     374protected: 
    374375        //! Internal epdf that arise by conditioning on \c rvc 
    375376        enorm<sq_T> epdf; 
    … …  
    379380public: 
    380381        //! Constructor 
    381         mlnorm (const RV &rv, const RV &rvc ); 
     382        mlnorm ( const RV &rv, const RV &rvc ); 
    382383        //! Set \c A and \c R 
    383384        void set_parameters ( const  mat &A, const vec &mu0, const sq_T &R ); 
    … …  
    387388//      mat samplecond (const vec &cond, vec &lik, int n ); 
    388389        //! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it. 
    389         void condition (const vec &cond ); 
    390          
     390        void condition ( const vec &cond ); 
     391 
     392        //!access function 
     393        vec& _mu_const() {return mu_const;} 
     394        //!access function 
     395        mat& _A() {return A;} 
     396        //!access function 
     397        mat _R() {return epdf._R().to_mat();} 
     398 
    391399        template<class sq_M> 
    392400        friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M> &ml ); 
    393401}; 
    394402 
     403/*! (Approximate) Student t density with linear function of mean value 
     404*/ 
     405class mlstudent : public mlnorm<ldmat> { 
     406protected: 
     407        ldmat Lambda; 
     408        ldmat &_R; 
     409        ldmat Re; 
     410public: 
     411        mlstudent ( const RV &rv0, const RV &rvc0 ) :mlnorm<ldmat> ( rv0,rvc0 ), 
     412                        Lambda ( rv0.count() ), 
     413                        _R ( epdf._R() ) {} 
     414        void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0) { 
     415                epdf.set_parameters ( zeros ( rv.count() ),Lambda ); 
     416                A = A0; 
     417                mu_const = mu0; 
     418                Re=R0; 
     419                Lambda = Lambda0; 
     420        } 
     421        void condition ( const vec &cond ) { 
     422                _mu = A*cond + mu_const; 
     423                double zeta; 
     424                //ugly hack! 
     425                if ((cond.length()+1)==Lambda.rows()){ 
     426                        zeta = Lambda.invqform ( concat(cond, vec_1(1.0)) ); 
     427                } else { 
     428                        zeta = Lambda.invqform ( cond ); 
     429                } 
     430                _R = Re; 
     431                _R*=( 1+zeta );// / ( nu ); << nu is in Re!!!!!! 
     432                 
     433                cout << _mu <<" +- " << sqrt(_R.to_mat()) << endl; 
     434        }; 
     435 
     436}; 
    395437/*! 
    396438\brief  Gamma random walk 
    … …  
    548590double enorm<sq_T>::evalpdflog_nn ( const vec &val ) const { 
    549591        // 1.83787706640935 = log(2pi) 
    550         return  -0.5* ( R.invqform ( mu-val ) );// - lognc(); 
     592        double tmp=-0.5* ( R.invqform ( mu-val ) );// - lognc(); 
     593        return  tmp; 
    551594}; 
    552595 
    … …  
    554597inline double enorm<sq_T>::lognc () const { 
    555598        // 1.83787706640935 = log(2pi) 
    556         return 0.5* ( R.cols() * 1.83787706640935 +R.logdet() ); 
     599        double tmp=0.5* ( R.cols() * 1.83787706640935 +R.logdet() ); 
     600        return tmp; 
    557601}; 
    558602 
    … …  
    584628//      vec smp ( dim ); 
    585629//      this->condition ( cond ); 
    586 //  
     630// 
    587631//      for ( i=0; i<n; i++ ) { 
    588632//              smp = epdf.sample(); 
    … …  
    590634//              Smp.set_col ( i ,smp ); 
    591635//      } 
    592 //  
     636// 
    593637//      return Smp; 
    594638// } 
    595639 
    596640template<class sq_T> 
    597 void mlnorm<sq_T>::condition (const vec &cond ) { 
     641void mlnorm<sq_T>::condition ( const vec &cond ) { 
    598642        _mu = A*cond + mu_const; 
    599643//R is already assigned; 
    … …  
    605649 
    606650        sq_T Rn ( R,irvn ); 
    607         enorm<sq_T>* tmp = new enorm<sq_T>( rvn ); 
     651        enorm<sq_T>* tmp = new enorm<sq_T> ( rvn ); 
    608652        tmp->set_parameters ( mu ( irvn ), Rn ); 
    609653        return tmp; 
    … …  
    627671        int end=R.rows()-1; 
    628672        mat S11 = S.get ( 0,n, 0, n ); 
    629         mat S12 = S.get (0, n , rvn.count(), end ); 
     673        mat S12 = S.get ( 0, n , rvn.count(), end ); 
    630674        mat S22 = S.get ( rvn.count(), end, rvn.count(), end ); 
    631675 
    … …  
    644688 
    645689template<class sq_T> 
    646 std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_T> &ml ){ 
     690std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_T> &ml ) { 
    647691        os << "A:"<< ml.A<<endl; 
    648692        os << "mu:"<< ml.mu_const<<endl; 

  • bdm/stat/loggers.h

    r162 r198  
    4646 
    4747//! log this vector 
    48         virtual void logit ( int id, vec v ) =0; 
     48        virtual void logit ( int id, const vec &v ) =0; 
    4949 
    5050//! Shifts storage position for another time step. 
    … …  
    8686        } 
    8787        void step() {if ( ind<maxlen ) ind++; else it_error ( "memlog::ind is too high;" );} 
    88         void logit ( int id, vec v ) { 
     88        void logit ( int id, const vec &v ) { 
    8989                it_assert_debug(id<vectors.length(),"Logger was not initialized, run init()."); 
    9090                vectors ( id ).set_row ( ind,v );} 

  • matlab/contour_2.m

    r194 r198  
    77 % 
    88 if nargin<4 
    9      contour(X,Y,Z) 
     9     contour(X,Y,Z,7) 
    1010 else 
    1111     contour(X,Y,Z,str) 

  • matlab/merger_iter_debug.m

    r193 r198  
    66YL= XL; 
    77 
     8iters = 0:9 
     9 
    810figure(1); 
    9 for it=0:4 
    10     figure(it+1) 
    11     subplot(2,2,1); 
     11niters = length(iters); 
     12noi = 1; 
     13for it=iters 
     14%    figure(it+1) 
     15    subplot(5,niters,noi); 
    1216    si = num2str(it); 
    1317    eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),exp(Mpdf' si '))']); 
    14     title('Proposal density '); 
     18    if it==iters(1), 
     19    ylabel('Proposal density '); 
     20    end 
    1521    set(gca,'XLim',XL); 
    1622    set(gca,'YLim',YL); 
    1723         
    18     subplot(2,2,2); 
     24    subplot(5,niters,noi+niters); 
    1925    eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),exp(lW' si '(1,:)))']); 
    20     title('First source');  
     26    if it==iters(1), 
     27    ylabel('First source');  
     28    end 
    2129    set(gca,'XLim',XL); 
    2230    set(gca,'YLim',YL); 
    2331     
    24     subplot(2,2,3); 
     32    subplot(5,niters,noi+2*niters); 
    2533    eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),exp(lW' si '(2,:)))']); 
    26     %eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),exp(lw_cond' si '))']); 
    27     title('Second source'); 
     34    if it==iters(1), 
     35    ylabel('Second source'); 
     36    end 
    2837    set(gca,'XLim',XL); 
    2938    set(gca,'YLim',YL); 
    3039     
    3140     
    32     subplot(2,2,4); 
     41    subplot(5,niters,noi+3*niters); 
    3342    hold off 
    34     eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),w' si ','':'')']); 
    35     title('Merged density');  
     43    eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),w' si ')']); 
     44    if it==iters(1), 
     45    ylabel('Merged density');  
     46    end 
    3647    set(gca,'XLim',XL); 
    3748    set(gca,'YLim',YL); 
    38     hold on 
     49     
     50    subplot(5,niters,noi+4*niters); 
     51    hold off 
    3952    eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),w_is_' si ')']); 
    40      
     53    if it==iters(1), 
     54    ylabel('Importance density');  
     55    end 
     56    set(gca,'XLim',XL); 
     57    set(gca,'YLim',YL); 
     58    noi = noi+1; 
    4159end 
    4260 
    43     
    44 itload('../merger_iter_test.it'); 
    45 XG = reshape(Grid(1,:),Npoints,Npoints); 
    46 YG = reshape(Grid(2,:),Npoints,Npoints); 
     61for it=0:25 
     62    si = num2str(it); 
     63    eval(['std_w(it+1) = std(w_is_' si '*200);']); 
     64    eval(['prop_mean(1:2,it+1) = Mpred_mean' si ';']); 
     65end 
    4766 
    48 figure(it+2); 
    49 M1 = reshape(exp(Res1),Npoints,Npoints); 
    50 contour(XG,YG,M1,7); 
     67figure(2) 
     68subplot(1,3,1) 
     69hold off 
     70plot(std_w) 
     71xlabel('iteration'); 
     72ylabel('standard deviation'); 
     73title('Non-normalized importance weights') 
    5174 
    52 figure(it+3); 
     75subplot(1,3,2); 
     76hist (w_is_15*200); 
     77title('Histogram of importance weights'); 
     78 
     79subplot(1,3,3); 
    5380hold off 
    54 mm = max(max(Res2)); 
    55 for i=1:size(Res2,1) 
    56     M2 = reshape(Res2(i,:),Npoints,Npoints); 
    57     contour(XG,YG,M2,[0:mm/7:mm]); 
    58     hold on 
    59 end 
     81plot(prop_mean'); 
     82hold on 
     83plot(ones(size(prop_mean,2),1)*[3 2],'--') 
     84xlabel('iteration'); 
     85title('Mean value of the merged density'); 
     86 
     87% itload('../merger_iter_test.it'); 
     88% XG = reshape(Grid(1,:),Npoints,Npoints); 
     89% YG = reshape(Grid(2,:),Npoints,Npoints); 
     90%  
     91% figure(2); 
     92% M1 = reshape(exp(Res1),Npoints,Npoints); 
     93% contour(XG,YG,M1,7); 
     94%  
     95% figure(3); 
     96% hold off 
     97% mm = max(max(Res2)); 
     98% for i=1:size(Res2,1) 
     99%     M2 = reshape(Res2(i,:),Npoints,Npoints); 
     100%     contour(XG,YG,M2,[0:mm/7:mm]); 
     101%     hold on 
     102% end 
     103 

  • mpdm/CMakeLists.txt

    r161 r198  
    2020target_link_libraries (merg_pred ${BdmLibs}) 
    2121 
     22add_executable (merger_iter_cond merger_iter_cond.cpp) 
     23target_link_libraries (merger_iter_cond ${BdmLibs}) 

  • mpdm/merg_pred.cpp

    r176 r198  
    11#include <estim/arx.h> 
     2#include <estim/merger.h> 
    23#include <stat/libEF.h> 
    34#include <stat/loggers.h> 
    … …  
    1415        RV u1 ( "{u1 }" ); 
    1516        RV u2 ( "{u2 }" ); 
     17        RV ym=y; ym.t(-1); 
     18        RV yy = y; yy.add(ym); 
     19        RV uu=u1; uu.add(u2); 
    1620 
    1721        // Full system 
    18         vec thg ( "0.7 1 1" ); //Simulated system 
     22        vec thg ( "0.7 1 1 0" ); //Simulated system - zero for constant term 
    1923        //y=a y_t-1 + u1 + u2 
    2024        double sqr=0.1; 
    … …  
    2226 
    2327        // Estimated systems ARX(2) 
    24         RV thri ( "{thr_i }",vec_1 ( 2+1 ) ); 
    25         RV thrg ( "{thr_g }",vec_1 ( 3+1 ) ); 
     28        RV thri ( "{thr_i }",vec_1 ( 3+1 ) ); 
     29        RV thrg ( "{thr_g }",vec_1 ( 4+1 ) ); 
    2630        // Setup values 
    2731 
    … …  
    2933        mat V0 = 0.001*eye ( thri.count() ); V0 ( 0,0 ) *= 10; // 
    3034        mat V0g = 0.001*eye ( thrg.count() ); V0g ( 0,0 ) *= 10; // 
    31         double nu0 = ord+1; 
    32         double frg = 0.99; 
     35        double nu0 = ord+6.0; 
     36        double frg = 1.0;//0.99; 
    3337 
    3438        ARX P1 ( thri, V0, nu0, frg ); 
    … …  
    3640        ARX PG ( thrg, V0g, nu0, frg ); 
    3741        //Test estimation 
    38         int ndat = 10000; 
     42        int ndat = 500; 
    3943        int t; 
    4044 
    … …  
    4246        dirfilelog L ( "exp/merg",ndat ); 
    4347 
    44         int Eth1_log = L.add ( thri,"P1" ); 
    45         int Eth2_log = L.add ( thri,"P2" ); 
    46         int Ethg_log = L.add ( thrg,"PG" ); 
    47         int Data_log = L.add ( RV ( "{Y U1 U2 }" ), "" ); 
    48         int LL_log   = L.add ( RV ( "{1 2 G }" ), "LL" ); 
     48        int Li_Eth1 = L.add ( thri,"P1" ); 
     49        int Li_Eth2 = L.add ( thri,"P2" ); 
     50        int Li_Ethg = L.add ( thrg,"PG" ); 
     51        int Li_Data = L.add ( RV ( "{Y U1 U2 }" ), "" ); 
     52        int Li_LL   = L.add ( RV ( "{1 2 G }" ), "LL" ); 
     53        int Li_Pred = L.add ( RV ( "{1 2 G ar ge ln}" ), "Pred" ); 
    4954 
    5055        L.init(); 
    5156 
    5257        vec Yt ( ndat ); 
     58        vec yt(1); 
    5359 
    5460        Yt.set_subvector ( 0,randn ( ord ) ); //initial values 
    55         vec rgrg ( thg.length() ); 
    56         vec rgri ( thri.count() ); 
     61        vec rgrg ( thrg.count() -1); // constant terms are in! 
     62        vec rgr1 ( thri.count() -1); 
     63        vec rgr2 ( thri.count() -1); 
    5764 
     65        vec PredLLs(5); 
     66        vec PostLLs(3); 
     67        vec PredLLs_m=zeros(5); 
     68        vec PostLLs_m=zeros(3); 
     69        ivec ind_r1 = "0 1 3"; 
     70        ivec ind_r2 = "0 2 3"; 
    5871        for ( t=0; t<ndat; t++ ) { 
    5972                // True system 
    … …  
    6275                        rgrg ( 1 ) = pow(sin ( ( t/40.0 ) *pi ),3); 
    6376                        rgrg ( 2 ) = pow(cos ( ( t/40.0 ) *pi ),3); 
     77                        rgrg (3) = 1.0; // constant term 
     78                         
     79                        rgr1(0) = rgrg(0); rgr1(1) = rgrg(1); rgr1(2) = rgrg(3); // no u2 
     80                        rgr2(0) = rgrg(0); rgr2(1) = rgrg(2); rgr2(2) = rgrg(3); // no u1 
    6481 
    6582                        Yt ( t ) = thg*rgrg + sqr * NorRNG(); 
    6683 
     84                        // Test predictors 
     85                        if (t>2){ 
     86                                mlnorm<ldmat>* P1p = P1.predictor_student(y,concat(ym,u1)); 
     87                                mlnorm<ldmat>* P2p = P2.predictor_student(y,concat(ym,u2)); 
     88                                mlnorm<ldmat>* Pgp = PG.predictor_student(y,concat(ym,uu)); 
     89                                 
     90                                Array<mpdf*> A(2); A(0)=P1p;A(1)=P2p; 
     91                                merger M(A); 
     92                                enorm<ldmat> g0(concat(yy,uu)); g0.set_parameters("0 0 0 0 ",3*eye(4)); 
     93                                M.set_parameters(1.2, 100,1); 
     94                                M.merge(&g0); 
     95                                 
     96                                yt(0) = Yt(t); 
     97                                double P1pl = P1p->evalcond(yt,rgr1);    
     98                                double P2pl = P2p->evalcond(yt,rgr2);    
     99                                double PGpl = Pgp->evalcond(yt,rgrg);    
     100                                PredLLs(0) = log(P1pl); PredLLs(1) = log(P2pl); PredLLs(2) = log(PGpl); 
     101                                { 
     102                                        ARX* Apred = (ARX*)M._Mix()._Coms(0); 
     103                                        enorm<ldmat>* MP= Apred->predictor(concat(yy,uu)); 
     104                                        enorm<ldmat>* mP1p = (enorm<ldmat>*)MP->marginal(concat(yy,u1)); 
     105                                        enorm<ldmat>* mP2p = (enorm<ldmat>*)MP->marginal(concat(yy,u2)); 
     106                                        mlnorm<ldmat>* cP1p = (mlnorm<ldmat>*)mP1p->condition(y); 
     107                                        mlnorm<ldmat>* cP2p = (mlnorm<ldmat>*)mP2p->condition(y); 
     108 
     109                                        double cP1pl = cP1p->evalcond(yt,rgr1);  
     110                                        double cP2pl = cP2p->evalcond(yt,rgr2);  
     111                                        PredLLs(3) = log(cP1pl); 
     112                                        PredLLs(4) = log(cP2pl); 
     113                                } 
     114 
     115                                L.logit(Li_Pred, PredLLs_m*frg+ PredLLs); //log-normal 
     116                                PredLLs_m *=frg; 
     117                                PredLLs_m += PredLLs; 
     118                                 
     119                                delete P1p; 
     120                                delete P2p; 
     121                                delete Pgp; 
     122                        } 
     123                         
    67124                        // 1st 
    68                         rgri ( 0 ) =Yt ( t ); 
    69                         rgri ( 1 ) =Yt ( t-1 ); 
    70                         rgri ( 2 ) =rgrg ( 1 ); 
    71                         P1.bayes ( rgri ); 
     125                        P1.bayes ( concat(Yt(t),rgr1) ); 
    72126                        // 2nd 
    73                         rgri ( 2 ) =rgrg ( 2 ); 
    74                         P2.bayes ( rgri ); 
     127                        P2.bayes ( concat(Yt(t),rgr2) ); 
    75128 
    76129                        //Global 
    … …  
    79132                        //Merger 
    80133                } 
    81                 L.logit ( Eth1_log, P1._epdf().mean() ); 
    82                 L.logit ( Eth2_log, P2._epdf().mean() ); 
    83                 L.logit ( Ethg_log, PG._epdf().mean() ); 
    84                 L.logit ( Data_log, vec_3 ( Yt ( t ), rgrg ( 1 ), rgrg ( 2 ) ) ); 
    85                 L.logit ( LL_log, vec_3 ( P1._ll(), P2._ll(), PG._ll() ) ); 
     134                L.logit ( Li_Eth1, P1._epdf().mean() ); 
     135                L.logit ( Li_Eth2, P2._epdf().mean() ); 
     136                L.logit ( Li_Ethg, PG._epdf().mean() ); 
     137                L.logit ( Li_Data, vec_3 ( Yt ( t ), rgrg ( 1 ), rgrg ( 2 ) ) ); 
     138                PostLLs *= frg; 
     139                PostLLs += vec_3 ( P1._ll(), P2._ll(), PG._ll() ); 
     140                L.logit ( Li_LL, PostLLs ); 
    86141                L.step (  ); 
    87142        } 

  • tests/arx_test.cpp

    r170 r198  
    1414         
    1515        //Test constructor 
    16         mat V0 = 0.00001*eye(ord+1); V0(0.0)*= 10; // 
    17         double nu0 = ord+1; 
     16        mat V0 = 0.00001*eye(ord+1); V0(0.0)*= 100; // 
     17        double nu0 = ord+4; 
    1818         
    1919        RV thr("{theta_and_r }",vec_1(ord+1)); 
    … …  
    2222                                 
    2323        //Test estimation 
    24         int ndat = 10000; 
     24        int ndat = 100; 
    2525        int t,j; 
    2626        vec Yt(ndat); 
    … …  
    2828        Yt.set_subvector(0,randn(ord)); //initial values 
    2929        vec rgr(ord); 
     30         
    3031         
    3132        cout << Ar_ep.mean()<<endl; 
    … …  
    3738                Ar.bayes(Psi); 
    3839                LL(t) = Ar._ll(); 
     40                 
     41                cout << "y: " << Yt(t) << endl; 
     42                mlstudent*      Pr = Ar.predictor_student(RV("{y }"),RV("{y1 y2 y3 y4 }")); 
     43                cout << Ar._ll() <<" , " << log(Pr->evalcond(vec_1(Yt(t)),rgr)) <<endl; 
     44                delete Pr; 
    3945        } 
    4046        cout << Ar_ep.mean()<<endl; 

  • tests/merger_iter_test.cpp

    r193 r198  
    1111int main() { 
    1212 
    13         //RNG_randomize(); 
     13        RNG_randomize(); 
    1414         
    1515        RV x ( "{x }","1" ); 
    … …  
    1919         
    2020        enorm<fsqmat> f1 ( xy ); 
    21         enorm<fsqmat> f2 ( x ); 
     21        enorm<fsqmat> f2 ( xy ); 
     22        enorm<fsqmat> f3(y); 
    2223                 
    23         f1.set_parameters ( "3 2",mat ( "0.5 0; 0 0.5" ) ); 
    24         f2.set_parameters ( "3",mat ( "0.5" ) ); 
     24        f1.set_parameters ( "4 3",mat ( "0.4 0.3; 0.3 0.4" ) ); 
     25        f2.set_parameters ( "1 3",mat ( "0.3 -0.2; -0.2 0.3" ) ); 
     26        f3.set_parameters ( "2",mat("0.4") ); 
    2527         
    26         Array<mpdf* > A ( 2 ); 
     28        Array<mpdf* > A ( 3 ); 
    2729        mepdf A1(f1); 
    2830        mepdf A2(f2); 
     31        mepdf A3(f3); 
    2932        A ( 0 ) =&A1; 
    3033        A ( 1 ) =&A2; 
     34        A ( 2 ) =&A3; 
    3135         
    3236        int Npoints=100; 
    … …  
    4246        merger M ( A ); 
    4347        enorm<fsqmat> g0(xy); 
    44         g0.set_parameters(vec("1 1"),mat("1 0; 0 1")); 
     48        g0.set_parameters(vec("4 4"),mat("1 0; 0 1")); 
    4549         
    46         M.set_parameters(1.2,100,1); 
     50        M.set_parameters(1.2,400,5); 
    4751        M.merge(&g0); 
    4852         

  • tests/merger_test.cpp

    r182 r198  
    1111int main() { 
    1212 
     13        RNG_randomize(); 
     14         
    1315        RV x ( "{x }","1" ); 
    1416 
    … …  
    4244        g0.set_parameters(vec("0.0"),mat("100.0")); 
    4345         
    44         M.set_parameters(1.2,1000,3); 
     46        M.set_parameters(1.2,200,3); 
    4547        M.merge(&g0); 
    4648         
    … …  
    4951         
    5052        it_file it("merger_test.it"); 
     53        it << Name("x_grid") << x_grid; 
    5154        it << Name("lf1") << l_f1; 
    5255        it << Name("lf2") << l_f2;