Changeset 198 for bdm/stat/libEF.h

Timestamp:

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

Author:

smidl

Message:

opravy + zavedeni studenta + zakomentovani debug v mergeru

Files:

• bdm/stat/libEF.h (modified) (14 diffs)

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;