Index: bdm/stat/libEF.cpp
===================================================================
--- bdm/stat/libEF.cpp (revision 214)
+++ bdm/stat/libEF.cpp (revision 225)
@@ -59,5 +59,5 @@
 	             - 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,20 +78,21 @@
 		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,5 +106,10 @@
 
 	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,5 +142,5 @@
 	}
 	double tmp=res-lognc();;
-	it_assert_debug(std::isfinite(tmp),"Infinite value");
+	it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
 	return tmp;
 }
@@ -152,6 +158,4 @@
 
 //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,5 +259,5 @@
 }
 
-void eEmp::set_samples (  const epdf* epdf0 ) {
+void eEmp::set_samples ( const epdf* epdf0 ) {
 	//it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
 	w=1;
Index: bdm/stat/libEF.h
===================================================================
--- bdm/stat/libEF.h (revision 215)
+++ bdm/stat/libEF.h (revision 225)
@@ -317,5 +317,5 @@
 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,4 +328,41 @@
 	void _param ( vec* &a, vec* &b ) {a=&alpha;b=&beta;};
 	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,8 +508,38 @@
 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,4 +573,34 @@
 };
 
+
+/*!
+\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 };
Index: bdm/stat/libBM.h
===================================================================
--- bdm/stat/libBM.h (revision 223)
+++ bdm/stat/libBM.h (revision 225)
@@ -220,6 +220,4 @@
 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,4 +254,6 @@
 	//!access function
 	epdf& _epdf() {return *ep;}
+	//!access function
+	epdf* _e() {return ep;}
 };
 
