Changeset 225 for bdm/stat/libEF.h

Timestamp:

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

Author:

smidl

Message:

prestavba mpf_test do inv-gamma

Files:

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

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 };