Index: library/bdm/estim/arx.h
===================================================================
--- library/bdm/estim/arx.h (revision 1176)
+++ library/bdm/estim/arx.h (revision 1189)
@@ -221,6 +221,6 @@
 	
 	void bayes(const vec &dt, const vec &psi){
-	  ldmat &Pbeta = est.P;
-	  ldmat &Palpha = alternative.P;
+//	  ldmat &Pbeta = est.P;
+//	  ldmat &Palpha = alternative.P;
 	  
 	  
Index: library/bdm/stat/exp_family.h
===================================================================
--- library/bdm/stat/exp_family.h (revision 1175)
+++ library/bdm/stat/exp_family.h (revision 1189)
@@ -388,6 +388,13 @@
 * \brief Gauss-inverse-Wishart density stored in LD form
 
-* For \f$p\f$-variate densities, given rv.count() should be \f$p    imes\f$ V.rows().
+* For \f$p\f$-variate densities, 
+* 
+* \f$ M,R \sim GiW( V_t, \nu_t) \propto |R|^{0.5\nu}\exp(-1/2 tr(R^{-1}[I,M] V_t [I;M'])) \f$
 *
+* Factorizes as:
+* \f$ M|R \sim N( \hat{M}, R \otimes Vz^{-1}) \propto |R|^{-0.5dim(psi)}\exp(-1/2 tr((M-\hat{M})R^{-1}(M-\hat{M})Vz) \f$
+* \f$ R \sim iW( \Lambda,\delta) \propto |R|^{-0.5(\nu - dim(psi))} |\Lambda|^{}\exp(-1/2 tr(R^{-1}\Lambda) \f$
+* where in standard notation \f$ |R|^{-0.5(\delta + p +1)}\f$, i.e. 
+* \f$ \delta = \nu-dim(psi) -p-1 \f$
 */
 class egiw : public eEF {
@@ -426,5 +433,5 @@
 
     void factorize ( mat &M, ldmat &Vz, ldmat &Lam ) const;
-    //! LS estimate of \f$    heta\f$
+    //! LS estimate of \f$\theta\f$
     vec est_theta() const;
 
@@ -1612,9 +1619,11 @@
 
 /*! \brief Inverse Wishart density defined on Choleski decomposition
+ * 
+ * here \f$\Omega \sim |\Sigma|^{-0.5(\delta}|\Omega|^{0.5(\delta -p -1)} \exp(-1/2 tr(\Omega\Sigma^{-1}))\$f
 */
 class eWishartCh : public epdf {
 protected:
     //! Upper-Triagle of Choleski decomposition of \f$ \Psi \f$
-    chmat Y;
+    chmat Sigma;
     //! dimension of matrix  \f$ \Psi \f$
     int p;
@@ -1624,13 +1633,13 @@
     //! Set internal structures
     void set_parameters ( const mat &Y0, const double delta0 ) {
-        Y = chmat ( Y0 );
+        Sigma = chmat ( Y0 );
         delta = delta0;
-        p = Y.rows();
+        p = Sigma.rows();
     }
     //! Set internal structures
     void set_parameters ( const chmat &Y0, const double delta0 ) {
-        Y = Y0;
+        Sigma = Y0;
         delta = delta0;
-        p = Y.rows();
+        p = Sigma.rows();
     }
 
@@ -1640,7 +1649,8 @@
     }
 
-    //! Sample matrix argument
-    mat sample_mat() const {
-        mat X = zeros ( p, p );
+    //! Sample matrix argument - lower 
+    //! Using Bartlet decomposition: W=L A A^T L^T where A is lower triag and L is choleski factor of Sigma.
+    chmat sample_mat() const {
+        mat A_T = zeros ( p, p ); // A transpose
 
         //sample diagonal
@@ -1648,5 +1658,5 @@
             GamRNG.setup ( 0.5* ( delta - i ) , 0.5 );   // no +1 !! index if from 0
 #pragma omp critical
-            X ( i, i ) = sqrt ( GamRNG() );
+            A_T ( i, i ) = sqrt ( GamRNG() );
         }
         //do the rest
@@ -1654,12 +1664,14 @@
             for ( int j = i + 1; j < p; j++ ) {
 #pragma omp critical
-                X ( i, j ) = NorRNG.sample();
+                A_T ( i, j ) = NorRNG.sample();
             }
         }
-        return X*Y._Ch();// return upper triangular part of the decomposition
+        chmat tmp;
+		tmp._Ch()=A_T*Sigma._Ch();
+        return tmp;// return upper triangular part of the decomposition
     }
 
     vec sample () const {
-        return vec ( sample_mat()._data(), p*p );
+        return vec ( sample_mat().to_mat()._data(), p*p );
     }
 
@@ -1673,15 +1685,15 @@
     //! fast access function y0 will be copied into Y.Ch.
     void setY ( const mat &Ch0 ) {
-        copy_vector ( dim, Ch0._data(), Y._Ch()._data() );
+        copy_vector ( dim, Ch0._data(), Sigma._Ch()._data() );
     }
 
     //! fast access function y0 will be copied into Y.Ch.
     void _setY ( const vec &ch0 ) {
-        copy_vector ( dim, ch0._data(), Y._Ch()._data() );
+        copy_vector ( dim, ch0._data(), Sigma._Ch()._data() );
     }
 
     //! access function
     const chmat& getY() const {
-        return Y;
+        return Sigma;
     }
 };
@@ -1715,5 +1727,5 @@
     vec sample() const {
         mat iCh;
-        iCh = inv ( W.sample_mat() );
+        iCh = inv ( W.sample_mat()._Ch() );
         return vec ( iCh._data(), dim );
     }
Index: library/bdm/stat/exp_family.cpp
===================================================================
--- library/bdm/stat/exp_family.cpp (revision 1102)
+++ library/bdm/stat/exp_family.cpp (revision 1189)
@@ -49,22 +49,24 @@
     ChLam.inv ( iChLam );
 
-    eWishartCh Omega; //inverse Wishart, result is R,
-    Omega.set_parameters ( iChLam, nu - 2*nPsi - dimx ); // 2*nPsi is there to match numercial simulations - check if analytically correct
+    eWishartCh Omega; //iWishart, same degrees of freedom as in Wishard in terms of delta
+	// delta = nu - npsi- p -1
+    Omega.set_parameters ( iChLam, nu - nPsi -dimx -1);
     Omega.validate();
 
-    mat OmChi;
+    chmat OmChi;
     mat Z ( M.rows(), M.cols() );
 
     mat Mi;
-    mat RChiT;
+    chmat Ri;
     mat tmp ( dimension(), n );
     M=M.T();// ugly hack == decide what to do with M.
     for ( int i = 0; i < n; i++ ) {
         OmChi = Omega.sample_mat();
-        RChiT = inv ( OmChi );
-        Z = randn ( M.rows(), M.cols() );
-        Mi = M + RChiT * Z * inv ( Vz._L().T() * diag ( sqrt ( Vz._D() ) ) );
-
-        tmp.set_col ( i, concat ( cvectorize ( Mi ), cvectorize ( RChiT*RChiT.T() ) ) );
+        OmChi.inv(Ri );
+		if (M._datasize()>0){
+			Z = randn ( M.rows(), M.cols() );
+			Mi = M + Ri._Ch().T() * Z * inv ( Vz._L().T() * diag ( sqrt ( Vz._D() ) ) );
+		}
+        tmp.set_col ( i, concat ( cvectorize ( Mi ), cvectorize ( Ri.to_mat() ) ) );
     }
     return tmp;
@@ -79,21 +81,21 @@
     factorize ( M, Vz, Lam );
 
-    chmat Ch;
-    Ch.setCh ( Lam._L() *diag ( sqrt ( Lam._D() ) ) );
+    chmat Ch(Lam.to_mat());
     chmat iCh;
     Ch.inv ( iCh );
 
     eWishartCh Omega; //inverse Wishart, result is R,
-    Omega.set_parameters ( iCh, nu - 2*nPsi - dimx ); // 2*nPsi is there to match numercial simulations - check if analytically correct
+    Omega.set_parameters ( iCh, nu - nPsi - dimx -1 ); // 2*nPsi is there to match numercial simulations - check if analytically correct
     Omega.validate();
 
     chmat Omi;
-    Omi.setCh ( Omega.sample_mat() );
-
+    Omi= Omega.sample_mat();
+	Omi.inv ( Ri );
+	
     if (M._datasize()>0) {
+		M = M.T();
         mat Z = randn ( M.rows(), M.cols() );
-        Mi = M + Omi._Ch() * Z * inv ( Vz._L() * diag ( sqrt ( Vz._D() ) ) );
-    }
-    Omi.inv ( Ri );
+        Mi = M + Ri._Ch().T() * Z * inv ( Vz._L() * diag ( sqrt ( Vz._D() ) ) );
+    }
 }
 
