Changeset 198 for bdm

Timestamp:

11/04/08 14:54:34 (16 years ago)

Author:

smidl

Message:

opravy + zavedeni studenta + zakomentovani debug v mergeru

Location:

bdm

Files:

• estim/arx.cpp (modified) (3 diffs)
• estim/arx.h (modified) (1 diff)
• estim/merger.cpp (modified) (9 diffs)
• estim/merger.h (modified) (3 diffs)
• estim/mixef.cpp (modified) (2 diffs)
• estim/mixef.h (modified) (1 diff)
• stat/libBM.h (modified) (1 diff)
• stat/libEF.cpp (modified) (2 diffs)
• stat/libEF.h (modified) (14 diffs)
• stat/loggers.h (modified) (2 diffs)

bdm/estim/arx.cpp

@@ -33 +33 @@
                 lll = pred.lognc();
         }
-        else 
-                if(evalll){lll=last_lognc;}else{lll=pred.lognc();}
+        else
+                if ( evalll ) {lll=last_lognc;}
+                else{lll=pred.lognc();}
 
         V.opupdt ( dt,1.0 );
@@ -42 +43 @@
 }
 
-ARX* ARX::_copy_ ( bool changerv) {
+ARX* ARX::_copy_ ( bool changerv ) {
         ARX* Tmp=new ARX ( *this );
-        if ( changerv ) {Tmp->rv.newids(); Tmp->est._renewrv(Tmp->rv);}
+        if ( changerv ) {Tmp->rv.newids(); Tmp->est._renewrv ( Tmp->rv );}
         return Tmp;
 }
@@ -55 +56 @@
 }
 
-enorm<ldmat>* ARX::predictor(const RV &rv){
-        mat mu(rv.count(), 1);
-        mat R(rv.count(),rv.count());
-        enorm<ldmat>* tmp; 
-        tmp=new enorm<ldmat>(rv);
+enorm<ldmat>* ARX::predictor ( const RV &rv, const vec &rgr ) const {
+        mat mu ( rv.count(), V.rows()-rv.count() );
+        mat R ( rv.count(),rv.count() );
+        enorm<ldmat>* tmp;
+        tmp=new enorm<ldmat> ( rv );
+
+        est.mean_mat ( mu,R ); //mu =
+        //correction for student-t  -- TODO check if correct!!
+        //R*=nu/(nu-2);
+        mat p_mu=mu.T() *rgr;   //the result is one column
+        tmp->set_parameters ( p_mu.get_col ( 0 ),ldmat ( R ) );
+        return tmp;
+}
+
+mlnorm<ldmat>* ARX::predictor ( const RV &rv, const RV &rvc ) const {
+        int dif=V.rows() - rv.count() - rvc.count();
+        it_assert_debug ( ( dif==0 ) || ( dif==1 ), "Give RVs do not match" );
+
+        mat mu ( rv.count(), V.rows()-rv.count() );
+        mat R ( rv.count(),rv.count() );
+        mlnorm<ldmat>* tmp;
+        tmp=new mlnorm<ldmat> ( rv,rvc );
+
+        est.mean_mat ( mu,R ); //mu =
+        mu = mu.T();
+        //correction for student-t  -- TODO check if correct!!
+
+        if ( dif==0 ) { // no constant term
+                tmp->set_parameters ( mu, zeros ( rv.count() ), ldmat ( R ) );
+        }
+        else {
+                //Assume the constant term is the last one:
+                tmp->set_parameters ( mu.get_cols (0,mu.cols()-2 ), mu.get_col ( mu.cols()-1 ), ldmat ( R ) );
+        }
+        return tmp;
+}
+
+mlstudent* ARX::predictor_student ( const RV &rv, const RV &rvc ) const {
+        int dif=V.rows() - rv.count() - rvc.count();
+        it_assert_debug ( ( dif==0 ) || ( dif==1 ), "Give RVs do not match" );
+
+        mat mu ( rv.count(), V.rows()-rv.count() );
+        mat R ( rv.count(),rv.count() );
+        mlstudent* tmp;
+        tmp=new mlstudent ( rv,rvc );
+
+        est.mean_mat ( mu,R ); //
+        mu = mu.T();
         
-        est.mean_mat(mu,R);
-        //here I know that mu can have only one column
-        tmp->set_parameters(mu.get_col(0),ldmat(R));
+        int xdim = rv.count();
+        int end = V._L().rows()-1;
+        ldmat Lam ( V._L() ( xdim,end,xdim,end ), V._D() ( xdim,end ) );  //exp val of R
+
+
+        if ( dif==0 ) { // no constant term
+                tmp->set_parameters ( mu, zeros ( rv.count() ), ldmat ( R ), Lam);
+        }
+        else {
+                //Assume the constant term is the last one:
+                tmp->set_parameters ( mu.get_cols (0,mu.cols()-2 ), mu.get_col ( mu.cols()-1 ), ldmat ( R ), Lam);
+        }
         return tmp;
 }

bdm/estim/arx.h

@@ -72 +72 @@
                 if(evalll){last_lognc=est.lognc();}
         }
-
-        enorm<ldmat>* predictor(const RV &rv); 
+        //! Conditional version of the predictor
+        enorm<ldmat>* predictor(const RV &rv0, const vec &rgr) const; 
+        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));}
+        mlnorm<ldmat>* predictor(const RV &rv0, const RV &rvc0) const;
+        mlstudent* predictor_student(const RV &rv0, const RV &rvc0) const;
         //! Brute force structure estimation.\return indeces of accepted regressors.
         ivec structure_est ( egiw Eg0 );

bdm/estim/merger.cpp

@@ -8 +8 @@
         vec lam = sum ( pow ( lW,2 ) )-nu*pow ( mu,2 );
         double coef=0.0;
+        vec sq2bl=sqrt ( 2*beta*lam ); //this term is everywhere
         switch ( nu ) {
                 case 2:
-                        coef=sqrt ( beta*2 ) * ( 1-0.5*sqrt ( ( 4*beta-3 ) /beta ) );
-                        return exp ( coef*sqrt ( lam ) + mu );
+                        coef= ( 1-0.5*sqrt ( ( 4*beta-3 ) /beta ) );
+                        return exp ( coef*sq2bl + mu );
                         break;
-                case 3://Ration of Bessel
+                case 3://Ratio of Bessel
+                        coef = sqrt ( ( 3*beta-2 ) /3*beta );
+                        return elem_mult (
+                                   elem_div ( besselk ( 0,sq2bl*coef ), besselk ( 0,sq2bl ) ),
+                                   exp ( mu ) );
                         break;
                 case 4:
@@ -24 +29 @@
 
 void merger::merge ( const epdf* g0 ) {
-        it_file dbg ( "merger_debug.it" );
+//      it_file dbg ( "merger_debug.it" );
 
         it_assert_debug ( rv.equal ( g0->_rv() ),"Incompatible g0" );
@@ -36 +41 @@
         for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );}
 
-        dbg << Name ( "Smp_0" ) << Smp_ex;
+//      dbg << Name ( "Smp_0" ) << Smp_ex;
 
         // Stuff for merging
@@ -51 +56 @@
         Mix.init ( &A0, Smp_ex, Nc );
         //Preserve initial mixture for repetitive estimation via flattening
-        MixEF Mix_init(Mix);
+        MixEF Mix_init ( Mix );
 
         // ============= MAIN LOOP ==================
@@ -63 +68 @@
                 //Re-estimate Mix
                 //Re-Initialize Mixture model
-                Mix.flatten(&Mix_init);
+                Mix.flatten ( &Mix_init );
                 Mix.bayesB ( Smp_ex, w*Ns );
                 Mpred = Mix.predictor ( rv ); // Allocation => must be deleted at the end!!
 
-                sprintf ( str,"Mpred_mean%d",niter );
-                dbg << Name ( str ) << Mpred->mean();
+//              sprintf ( str,"Mpred_mean%d",niter );
+//              dbg << Name ( str ) << Mpred->mean();
 
-                // Generate new samples
-                eSmp.set_samples ( Mpred );
-                for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );}
+                if ( niter<10 ) {
+                        // Generate new samples
+                        eSmp.set_samples ( Mpred );
+                        for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );}
 
-                sprintf ( str,"Mpdf%d",niter );
-                for ( int i=0;i<Ns;i++ ) {Mix_pdf ( i ) = Mix.logpred ( Smp_ex.get_col ( i ) );}
-                dbg << Name ( str ) << Mix_pdf;
+                }
+//              sprintf ( str,"Mpdf%d",niter );
+//              for ( int i=0;i<Ns;i++ ) {Mix_pdf ( i ) = Mix.logpred ( Smp_ex.get_col ( i ) );}
+//              dbg << Name ( str ) << Mix_pdf;
 
-                sprintf ( str,"Smp%d",niter );
-                dbg << Name ( str ) << Smp_ex;
+//              sprintf ( str,"Smp%d",niter );
+//              dbg << Name ( str ) << Smp_ex;
 
                 //Importace weighting
@@ -99 +106 @@
                                         //compute vector of lw_src
                                         for ( int k=0;k<Ns;k++ ) {
-                                                lw_src ( k ) += tmp_marg->evalpdflog ( dls ( i )->get_val ( Smp ( i ) ) );
+                                                // Here val of tmp_marg = cond of mpdfs(i) ==> calling dls->get_cond
+                                                lw_src ( k ) += tmp_marg->evalpdflog ( dls ( i )->get_cond ( Smp ( k ) ) );
                                         }
                                         delete tmp_marg;
+
+//                                      sprintf ( str,"marg%d",niter );
+//                                      dbg << Name ( str ) << lw_src;
+
                                 }
                                 // Compute likelihood of the missing variable
@@ -125 +137 @@
 
                         }
+//                      it_assert_debug(std::isfinite(sum(lw_src)),"bad");
                         lW.set_row ( i, lw_src ); // do not divide by mix
                 }
@@ -131 +144 @@
                         lw_mix ( j ) =Mix.logpred ( Smp_ex.get_col ( j ) );
                 }
-                sprintf ( str,"lW%d",niter );
-                dbg << Name ( str ) << lW;
+//              sprintf ( str,"lW%d",niter );
+//              dbg << Name ( str ) << lW;
 
                 w = lognorm_merge ( lW ); //merge
 
-                sprintf ( str,"w%d",niter );
-                dbg << Name ( str ) << w;
-                sprintf ( str,"lw_m%d",niter );
-                dbg << Name ( str ) << lw_mix;
+//              sprintf ( str,"w%d",niter );
+//              dbg << Name ( str ) << w;
+//              sprintf ( str,"lw_m%d",niter );
+//              dbg << Name ( str ) << lw_mix;
 
                 //Importance weighting
@@ -146 +159 @@
                 w /=sum ( w );
 
-                sprintf ( str,"w_is_%d",niter );
-                dbg << Name ( str ) << w;
+//              sprintf ( str,"w_is_%d",niter );
+//              dbg << Name ( str ) << w;
 
 //              eSmp.resample(); // So that it can be used in bayes
 //              for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );}
 
-                sprintf ( str,"Smp_res%d",niter );
-                dbg << Name ( str ) << Smp;
+//              sprintf ( str,"Smp_res%d",niter );
+//              dbg << Name ( str ) << Smp;
 
                 // ==== stopping rule ===
                 niter++;
-                converged = ( niter>9 );
+                converged = ( niter>5 );
         }
 

bdm/estim/merger.h

@@ -35 +35 @@
         Array<RV> rvzs;
         //! Data Links of rv0 mpdfs - these will be conditioned the (rv,rvc) of mpdfs
-        Array<datalink_m2e*> zdls; 
-        
+        Array<datalink_m2e*> zdls;
+
         //!Number of samples used in approximation
         int Ns;
@@ -47 +47 @@
         merger ( const Array<mpdf*> &S ) :
                         compositepdf ( S ), epdf ( getrv ( false ) ),
-                        Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ), dls(n), rvzs(n), zdls(n) {
-                                RV ztmp;
+                        Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ), dls ( n ), rvzs ( n ), zdls ( n ) {
+                RV ztmp;
+                // Extend rv by rvc!
+                RV rvc; setrvc ( rv,rvc );
+                rv.add ( rvc );
                 for ( int i=0;i<n;i++ ) {
                         //Establich connection between mpdfs and merger
-                        dls ( i ) = new datalink_m2e ( mpdfs ( i )->_rv(), mpdfs(i)->_rvc(), rv );
+                        dls ( i ) = new datalink_m2e ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), rv );
                         // find out what is missing in each mpdf
-                        ztmp= mpdfs(i)->_rv();
-                        ztmp.add(mpdfs(i)->_rvc());
-                        rvzs(i)=rv.subt(ztmp);
-                        zdls(i) = new datalink_m2e (rvzs(i), ztmp, rv ) ;
+                        ztmp= mpdfs ( i )->_rv();
+                        ztmp.add ( mpdfs ( i )->_rvc() );
+                        rvzs ( i ) =rv.subt ( ztmp );
+                        zdls ( i ) = new datalink_m2e ( rvzs ( i ), ztmp, rv ) ;
                 };
                 //Set Default values of parameters
@@ -90 +93 @@
 //! for future use
         virtual ~merger() {
-                for (int i=0; i<n; i++){
-                        delete dls(i);
-                        delete zdls(i);                 
+                for ( int i=0; i<n; i++ ) {
+                        delete dls ( i );
+                        delete zdls ( i );
                 }
         };

bdm/estim/mixef.cpp

@@ -132 +132 @@
 }
 
-emix* MixEF::predictor ( const RV &rv ) {
+emix* MixEF::predictor ( const RV &rv ) const {
         Array<epdf*> pC ( n );
         for ( int i=0;i<n;i++ ) {pC ( i ) =Coms ( i )->predictor ( rv );}
@@ -149 +149 @@
                 Coms(i)->flatten(Mix2->Coms(i));
         }
-        //Faltten weights
-        weights.flatten(&(Mix2->weights));
+        //Flatten weights = make them equal!!
+        weights.set_statistics(&(Mix2->weights));
 }

bdm/estim/mixef.h

@@ -105 +105 @@
         double logpred ( const vec &dt ) const;
         const epdf& _epdf() const {return *est;}
-        emix* predictor(const RV &rv);
+        emix* predictor(const RV &rv) const;
         //! Flatten the density as if it was not estimated from the data
         void flatten(const BMEF* M2);

bdm/stat/libBM.h

@@ -426 +426 @@
 
         //!Constructs a predictive density (marginal density on data)
-        virtual epdf* predictor ( const RV &rv ) {it_error ( "Not implemented" );return NULL;};
+        virtual epdf* predictor ( const RV &rv ) const {it_error ( "Not implemented" );return NULL;};
 
         //! Destructor for future use;

bdm/stat/libEF.cpp

@@ -82 +82 @@
                 mat R;
                 mean_mat(M,R);
-                return cvectorize (concat_vertical(M.T(),R));
+                return cvectorize (concat_vertical(M,R));
         }
 
@@ -95 +95 @@
         
         // set mean value
-        M= L ( xdim,end,0,xdim-1 ).T()*iLsub;
+        mat Lpsi = L ( xdim,end,0,xdim-1 );
+        M= iLsub*Lpsi;
         R= ldR.to_mat()  ;
 }

bdm/stat/libEF.h

@@ -94 +94 @@
         //!Flatten the posterior as if to keep nu0 data
 //      virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
-        
+
         BMEF* _copy_ ( bool changerv=false ) {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
 };
@@ -192 +192 @@
         //! returns a pointer to the internal statistics. Use with Care!
         double& _nu() {return nu;}
-        void pow ( double p ){V*=p;nu*=p;};
+        void pow ( double p ) {V*=p;nu*=p;};
 };
 
@@ -239 +239 @@
         //! Conjugate prior and posterior
         eDirich est;
-        //! Pointer inside est to sufficient statistics 
+        //! Pointer inside est to sufficient statistics
         vec β
 public:
@@ -271 +271 @@
                 // sum(beta) should be equal to sum(B.beta)
                 const vec &Eb=E->beta;//const_cast<multiBM*> ( E )->_beta();
-                beta*= ( sum ( Eb) /sum ( beta ) );
+                beta*= ( sum ( Eb ) /sum ( beta ) );
                 if ( evalll ) {last_lognc=est.lognc();}
         }
@@ -372 +372 @@
 template<class sq_T>
 class mlnorm : public mEF {
+protected:
         //! Internal epdf that arise by conditioning on \c rvc
         enorm<sq_T> epdf;
@@ -379 +380 @@
 public:
         //! Constructor
-        mlnorm (const RV &rv, const RV &rvc );
+        mlnorm ( const RV &rv, const RV &rvc );
         //! Set \c A and \c R
         void set_parameters ( const  mat &A, const vec &mu0, const sq_T &R );
@@ -387 +388 @@
 //      mat samplecond (const vec &cond, vec &lik, int n );
         //! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
-        void condition (const vec &cond );
-        
+        void condition ( const vec &cond );
+
+        //!access function
+        vec& _mu_const() {return mu_const;}
+        //!access function
+        mat& _A() {return A;}
+        //!access function
+        mat _R() {return epdf._R().to_mat();}
+
         template<class sq_M>
         friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M> &ml );
 };
 
+/*! (Approximate) Student t density with linear function of mean value
+*/
+class mlstudent : public mlnorm<ldmat> {
+protected:
+        ldmat Lambda;
+        ldmat &_R;
+        ldmat Re;
+public:
+        mlstudent ( const RV &rv0, const RV &rvc0 ) :mlnorm<ldmat> ( rv0,rvc0 ),
+                        Lambda ( rv0.count() ),
+                        _R ( epdf._R() ) {}
+        void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0) {
+                epdf.set_parameters ( zeros ( rv.count() ),Lambda );
+                A = A0;
+                mu_const = mu0;
+                Re=R0;
+                Lambda = Lambda0;
+        }
+        void condition ( const vec &cond ) {
+                _mu = A*cond + mu_const;
+                double zeta;
+                //ugly hack!
+                if ((cond.length()+1)==Lambda.rows()){
+                        zeta = Lambda.invqform ( concat(cond, vec_1(1.0)) );
+                } else {
+                        zeta = Lambda.invqform ( cond );
+                }
+                _R = Re;
+                _R*=( 1+zeta );// / ( nu ); << nu is in Re!!!!!!
+                
+                cout << _mu <<" +- " << sqrt(_R.to_mat()) << endl;
+        };
+
+};
 /*!
 \brief  Gamma random walk
@@ -548 +590 @@
 double enorm<sq_T>::evalpdflog_nn ( const vec &val ) const {
         // 1.83787706640935 = log(2pi)
-        return  -0.5* ( R.invqform ( mu-val ) );// - lognc();
+        double tmp=-0.5* ( R.invqform ( mu-val ) );// - lognc();
+        return  tmp;
 };
 
@@ -554 +597 @@
 inline double enorm<sq_T>::lognc () const {
         // 1.83787706640935 = log(2pi)
-        return 0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
+        double tmp=0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
+        return tmp;
 };
 
@@ -584 +628 @@
 //      vec smp ( dim );
 //      this->condition ( cond );
-// 
+//
 //      for ( i=0; i<n; i++ ) {
 //              smp = epdf.sample();
@@ -590 +634 @@
 //              Smp.set_col ( i ,smp );
 //      }
-// 
+//
 //      return Smp;
 // }
 
 template<class sq_T>
-void mlnorm<sq_T>::condition (const vec &cond ) {
+void mlnorm<sq_T>::condition ( const vec &cond ) {
         _mu = A*cond + mu_const;
 //R is already assigned;
@@ -605 +649 @@
 
         sq_T Rn ( R,irvn );
-        enorm<sq_T>* tmp = new enorm<sq_T>( rvn );
+        enorm<sq_T>* tmp = new enorm<sq_T> ( rvn );
         tmp->set_parameters ( mu ( irvn ), Rn );
         return tmp;
@@ -627 +671 @@
         int end=R.rows()-1;
         mat S11 = S.get ( 0,n, 0, n );
-        mat S12 = S.get (0, n , rvn.count(), end );
+        mat S12 = S.get ( 0, n , rvn.count(), end );
         mat S22 = S.get ( rvn.count(), end, rvn.count(), end );
 
@@ -644 +688 @@
 
 template<class sq_T>
-std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_T> &ml ){
+std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_T> &ml ) {
         os << "A:"<< ml.A<<endl;
         os << "mu:"<< ml.mu_const<<endl;

bdm/stat/loggers.h

@@ -46 +46 @@
 
 //! log this vector
-        virtual void logit ( int id, vec v ) =0;
+        virtual void logit ( int id, const vec &v ) =0;
 
 //! Shifts storage position for another time step.
@@ -86 +86 @@
         }
         void step() {if ( ind<maxlen ) ind++; else it_error ( "memlog::ind is too high;" );}
-        void logit ( int id, vec v ) {
+        void logit ( int id, const vec &v ) {
                 it_assert_debug(id<vectors.length(),"Logger was not initialized, run init().");
                 vectors ( id ).set_row ( ind,v );}