Changeset 225

Timestamp:

01/12/09 13:28:18 (15 years ago)

Author:

smidl

Message:

prestavba mpf_test do inv-gamma

Files:

• bdm/estim/libKF.cpp (modified) (2 diffs)
• bdm/estim/libPF.h (modified) (2 diffs)
• bdm/stat/libBM.h (modified) (2 diffs)
• bdm/stat/libEF.cpp (modified) (6 diffs)
• bdm/stat/libEF.h (modified) (4 diffs)
• matlab/pmsm/mpf_test_disp.m (modified) (2 diffs)

bdm/estim/libKF.cpp

@@ -202 +202 @@
 
 //      mat Sttm = _P->to_mat();
+//      cout << preA <<endl;
+//      cout << "_mu:" << _mu <<endl;
 
 //      if ( !qr ( preA,Q,postA ) ) it_warning ( "QR in kalman unstable!" );
         if ( !qr ( preA,postA ) ) {it_warning ( "QR in kalman unstable!" );}
+
 
         ( _Ry._Ch() ) =postA ( 0,dimy-1, 0,dimy-1 );
@@ -229 +232 @@
 //              cout << "P:" <<_P.to_mat() <<endl;
 //              cout << "Ry:" <<_Ry.to_mat() <<endl;
-//              cout << "_K:" <<_K <<endl;
+//      cout << "_mu:" <<_mu <<endl;
+//      cout << "dt:" <<dt <<endl;
 
         if ( evalll==true ) { //likelihood of observation y

bdm/estim/libPF.h

@@ -52 +52 @@
         void set_est ( const epdf &epdf0 );
         void bayes ( const vec &dt );
+        //!access function
+        vec* __w(){return &_w;}
 };
 
@@ -131 +133 @@
         }
 
+        //!Access function
+        BM* _BM(int i){return Bms[i];}
 //SimStr:
         double SSAT;

bdm/stat/libBM.h

@@ -220 +220 @@
 public:
 
-        //! Returns the required moment of the epdf
-//      virtual fnc moment ( const int order = 1 );
         //! Returns a sample from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv \param ll is a return value of log-likelihood of the sample.
         virtual vec samplecond ( const vec &cond, double &ll ) {
@@ -256 +254 @@
         //!access function
         epdf& _epdf() {return *ep;}
+        //!access function
+        epdf* _e() {return ep;}
 };
 

bdm/stat/libEF.cpp

@@ -59 +59 @@
                      - 0.5*xdim* ( m*log2 + 0.5* ( xdim-1 ) *log2pi )  - lg;
 
-        it_assert_debug(((-nkG-nkW)>-Inf) && ((-nkG-nkW)<Inf), "ARX improper");
+        it_assert_debug ( ( ( -nkG-nkW ) >-Inf ) && ( ( -nkG-nkW ) <Inf ), "ARX improper" );
         return -nkG-nkW;
 }
@@ -78 +78 @@
                 m ( end ) = D ( 0 ) / ( nu -nPsi -2*xdim -2 );
                 return m;
-        } else {
+        }
+        else {
                 mat M;
                 mat R;
-                mean_mat(M,R);
-                return cvectorize (concat_vertical(M,R));
-        }
-
-}
-void egiw::mean_mat(mat &M, mat&R) const {
+                mean_mat ( M,R );
+                return cvectorize ( concat_vertical ( M,R ) );
+        }
+
+}
+void egiw::mean_mat ( mat &M, mat&R ) const {
         const mat &L= V._L();
         const vec &D= V._D();
         int end = L.rows()-1;
-                
+
         ldmat ldR ( L ( 0,xdim-1,0,xdim-1 ), D ( 0,xdim-1 ) / ( nu -nPsi -2*xdim -2 ) ); //exp val of R
         mat iLsub=ltuinv ( L ( xdim,end,xdim,end ) );
-        
+
         // set mean value
         mat Lpsi = L ( xdim,end,0,xdim-1 );
@@ -105 +106 @@
 
         for ( i=0; i<rv.count(); i++ ) {
-                GamRNG.setup ( alpha ( i ),beta ( i ) );
+                if ( beta ( i ) >std::numeric_limits<double>::epsilon() ) {
+                        GamRNG.setup ( alpha ( i ),beta ( i ) );
+                }
+                else {
+                        GamRNG.setup ( alpha ( i ),std::numeric_limits<double>::epsilon() );
+                }
 #pragma omp critical
                 smp ( i ) = GamRNG();
@@ -136 +142 @@
         }
         double tmp=res-lognc();;
-        it_assert_debug(std::isfinite(tmp),"Infinite value");
+        it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
         return tmp;
 }
@@ -152 +158 @@
 
 //MGamma
-mgamma::mgamma ( const RV &rv,const RV &rvc ) : mEF ( rv,rvc ), epdf ( rv ) {vec* tmp; epdf._param ( tmp,_beta );};
-
 void mgamma::set_parameters ( double k0 ) {
         k=k0;
@@ -255 +259 @@
 }
 
-void eEmp::set_samples (  const epdf* epdf0 ) {
+void eEmp::set_samples ( const epdf* epdf0 ) {
         //it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
         w=1;

bdm/stat/libEF.h

@@ -317 +317 @@
 public :
         //! Default constructor
-        egamma ( const RV &rv ) :eEF ( rv ) {};
+        egamma ( const RV &rv ) :eEF ( rv ), alpha(rv.count()), beta(rv.count()) {};
         //! Sets parameters
         void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;};
@@ -328 +328 @@
         void _param ( vec* &a, vec* &b ) {a=αb=β};
         vec mean() const {vec pom ( alpha ); pom/=beta; return pom;}
+};
+
+/*!
+ \brief Inverse-Gamma posterior density
+
+ Multivariate inverse-Gamma density as product of independent univariate densities.
+ \f[
+ f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i)
+ \f]
+
+ Inverse Gamma can be converted to Gamma using \[
+ x\sim iG(a,b) => 1/x\sim G(a,1/b)
+\]
+This relation is used in sampling.
+ */
+
+class eigamma : public eEF {
+        protected:
+        //! Vector \f$\alpha\f$
+                vec* alpha;
+        //! Vector \f$\beta\f$ (in fact it is 1/beta as used in definition of iG)
+                vec* beta;
+                //!internal egamma
+                egamma eg;
+        public :
+        //! Default constructor
+                eigamma ( const RV &rv ) :eEF ( rv ), eg(rv) {eg._param(alpha,beta);};
+        //! Sets parameters
+                void set_parameters ( const vec &a, const vec &b ) {*alpha=a,*beta=b;};
+                vec sample() const {return 1.0/eg.sample();};
+        //! TODO: is it used anywhere?
+//      mat sample ( int N ) const;
+                double evallog ( const vec &val ) const {return eg.evallog(val);};
+                double lognc () const {return eg.lognc();};
+        //! Returns poiter to alpha and beta. Potentially dangerous: use with care!
+                void _param ( vec* &a, vec* &b ) {a=alpha;b=beta;};
+                vec mean() const {return elem_div(*beta,*alpha-1);}
 };
 /*
@@ -471 +508 @@
 public:
         //! Constructor
-        mgamma ( const RV &rv,const RV &rvc );
+        mgamma ( const RV &rv,const RV &rvc ): mEF ( rv,rvc ), epdf ( rv ) {vec* tmp; epdf._param ( tmp,_beta );ep=&epdf;};
         //! Set value of \c k
         void set_parameters ( double k );
         void condition ( const vec &val ) {*_beta=k/val;};
+};
+
+/*!
+\brief  Inverse-Gamma random walk
+
+Mean value, \f$\mu\f$, of this density is given by \c rvc .
+Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
+This is achieved by setting \f$\alpha=\mu/k+2\f$ and \f$\beta=\mu(\alpha-1)\f$.
+
+The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
+ */
+class migamma : public mEF {
+        protected:
+        //! Internal epdf that arise by conditioning on \c rvc
+                eigamma epdf;
+        //! Constant \f$k\f$
+                double k;
+        //! cache of epdf.beta
+                vec* _beta;
+                //! chaceh of epdf.alpha
+                vec* _alpha;
+
+        public:
+        //! Constructor
+                migamma ( const RV &rv,const RV &rvc ): mEF ( rv,rvc ), epdf ( rv ) {epdf._param ( _alpha,_beta );ep=&epdf;};
+        //! Set value of \c k
+                void set_parameters ( double k0 ){k=k0;*_alpha=1.0/(k*k)+2;};
+                void condition ( const vec &val ) {
+                        *_beta=elem_mult(val,(*_alpha-1));
+                };
 };
 
@@ -506 +573 @@
 };
 
+
+/*!
+\brief  Inverse-Gamma random walk around a fixed point
+
+Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
+\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
+
+==== Check == vv =
+Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
+This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
+
+The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
+ */
+class migamma_fix : public migamma {
+        protected:
+        //! parameter l
+                double l;
+        //! reference vector
+                vec refl;
+        public:
+        //! Constructor
+                migamma_fix ( const RV &rv,const RV &rvc ) : migamma ( rv,rvc ),refl ( rv.count() ) {};
+        //! Set value of \c k
+                void set_parameters ( double k0 , vec ref0, double l0 ) {
+                        migamma::set_parameters ( k0 );
+                        refl=pow ( ref0,1.0-l0 );l=l0;
+                };
+
+                void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); migamma::condition(mean);};
+};
 //! Switch between various resampling methods.
 enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };

matlab/pmsm/mpf_test_disp.m

@@ -1 +1 @@
 close all
-itload('../mpf_test.it')
+itload('../../mpf_test.it')
 figure(1);
 for i =1:4
@@ -22 +22 @@
         plot(xthM(:,i)')
     if i<4
-        set(gca,'XTick',[]);
+        %set(gca,'XTick',[]);
     else
         xlabel('sample [t]');