Index: library/bdm/math/chmat.h
===================================================================
--- library/bdm/math/chmat.h (revision 723)
+++ library/bdm/math/chmat.h (revision 725)
@@ -57,4 +57,5 @@
 	void inv ( chmat &Inv ) const	{
 		( Inv.Ch ) = itpp::inv ( Ch ).T();
+		Inv.dim = dim;
 	}; //Fixme: can be more efficient
 	;
Index: library/bdm/math/square_mat.cpp
===================================================================
--- library/bdm/math/square_mat.cpp (revision 565)
+++ library/bdm/math/square_mat.cpp (revision 725)
@@ -115,6 +115,6 @@
 		}
 	}
-	mat V2 = L.transpose() * diag ( D ) * L;
-	return V2;
+	//mat V2 = L.transpose() * diag ( D ) * L;
+	return V;
 }
 
Index: library/bdm/stat/exp_family.h
===================================================================
--- library/bdm/stat/exp_family.h (revision 723)
+++ library/bdm/stat/exp_family.h (revision 725)
@@ -223,5 +223,7 @@
 		vec mean() const;
 		vec variance() const;
-
+		void sample_mat(mat &Mi, chmat &Ri)const;
+			
+		void factorize(mat &M, ldmat &Vz, ldmat &Lam) const;
 		//! LS estimate of \f$\theta\f$
 		vec est_theta() const;
@@ -245,4 +247,5 @@
 		const double& _nu() const {return nu;}
 		const int & _dimx() const {return dimx;}
+			
 		/*! Create Gauss-inverse-Wishart density 
 		\f[ f(rv) = GiW(V,\nu) \f]
@@ -259,5 +262,5 @@
 		\endcode
 		*/
-		
+				
 		void from_setting (const Setting &set) {
 			epdf::from_setting(set);
@@ -274,6 +277,38 @@
 			}
 		}
+		
+		void to_setting ( Setting& set ) const{
+			epdf::to_setting(set);
+			UI::save(dimx,set,"dimx");
+			UI::save(V.to_mat(),set,"V");
+			UI::save(nu,set,"nu");			
+		};
+		
 		void validate(){
 			// check sizes, rvs etc.
+		}
+		void log_register( bdm::logger& L, const string& prefix ){
+			if (log_level==3){
+				root::log_register(L,prefix);
+				logrec->ids.set_length(2);
+				int th_dim=dimension()-dimx*(dimx+1)/2;
+				logrec->ids(0)=L.add(RV("",th_dim), prefix + logrec->L.prefix_sep() +"mean");
+				logrec->ids(1)=L.add(RV("",th_dim*th_dim),prefix + logrec->L.prefix_sep() + "variance"); 
+			} else {
+				epdf::log_register(L,prefix);
+			}
+		}
+		void log_write() const {
+			if (log_level==3){
+				mat M;
+				ldmat Lam;
+				ldmat Vz;
+				factorize(M,Vz,Lam);
+				logrec->L.logit(logrec->ids(0), est_theta() );
+				logrec->L.logit(logrec->ids(1), cvectorize(est_theta_cov().to_mat()));
+			} else {
+				epdf::log_write();
+			}
+			
 		}
 		//!@}
@@ -1138,4 +1173,6 @@
 		//! Set internal structures
 		void set_parameters (const mat &Y0, const double delta0) {Y = chmat (Y0);delta = delta0; p = Y.rows(); dim = p * p; }
+		//! Set internal structures
+		void set_parameters (const chmat &Y0, const double delta0) {Y = Y0;delta = delta0; p = Y.rows(); dim = p * p; }
 		//! Sample matrix argument
 		mat sample_mat() const {
Index: library/bdm/stat/exp_family.cpp
===================================================================
--- library/bdm/stat/exp_family.cpp (revision 679)
+++ library/bdm/stat/exp_family.cpp (revision 725)
@@ -33,6 +33,33 @@
 
 vec egiw::sample() const {
-	bdm_warning ( "Function not implemented" );
-	return vec_1 ( 0.0 );
+	mat M;
+	chmat R;
+	sample_mat(M,R);
+	
+	return concat (cvectorize(M),cvectorize(R.to_mat()));
+}
+
+void egiw::sample_mat(mat &Mi, chmat &Ri)const{
+	
+	// TODO - correct approach - convert to product of norm * Wishart
+	mat M;
+	ldmat Vz;
+	ldmat Lam;
+	factorize(M,Vz,Lam);
+	
+	chmat Ch;
+	Ch.setCh(Lam._L()*diag(sqrt(Lam._D())));
+	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
+	
+	chmat Omi;
+	Omi.setCh(Omega.sample_mat());
+	
+	mat Z=randn(M.rows(), M.cols());
+	Mi = M+ Omi._Ch() * Z * inv(Vz._L()*diag(sqrt(Vz._D())));
+	Omi.inv(Ri);
 }
 
@@ -104,4 +131,26 @@
 }
 
+void egiw::factorize(mat &M, ldmat &Vz, ldmat &Lam) const{	
+	const mat &L = V._L();
+	const vec &D = V._D();
+	int end = L.rows() - 1;
+	
+	Vz=ldmat(L ( dimx, end, dimx, end ), D(dimx,end));
+	mat iLsub = ltuinv ( Vz._L());
+	// set mean value
+	mat Lpsi = L ( dimx, end, 0, dimx - 1 );
+	M = iLsub * Lpsi;
+	
+	Lam =ldmat ( L (0, dimx-1, 0, dimx-1 ), D (0, dimx-1 ) );  //exp val of R
+	 if (1){ // test with Peterka
+		 mat VF=V.to_mat();
+		 mat Vf=VF(0,dimx-1,0, dimx-1);
+		 mat Vzf = VF(dimx,end,0,dimx-1);
+		 mat VZ = VF(dimx,end,dimx,end);
+		 
+		 mat Lam2 = Vf-Vzf.T()*inv(VZ)*Vzf;
+	 }
+}
+
 ldmat egiw::est_theta_cov() const {
 	if ( dimx == 1 ) {
@@ -111,7 +160,8 @@
 
 		mat Lsub = L ( 1, end, 1, end );
-		mat Dsub = diag ( D ( 1, end ) );
-
-		return inv ( transpose ( Lsub ) * Dsub * Lsub );
+//		mat Dsub = diag ( D ( 1, end ) );
+
+		ldmat LD(inv(Lsub).T(), 1.0/D(1,end));
+		return LD;
 
 	} else {
