Index: library/bdm/stat/merger.cpp
===================================================================
--- library/bdm/stat/merger.cpp (revision 737)
+++ library/bdm/stat/merger.cpp (revision 739)
@@ -8,4 +8,85 @@
 		Npoints ( 0 ), DBG ( false ), dbg_file ( 0 ) {
 	set_sources ( S );
+}
+
+void merger_base::set_sources ( const Array<shared_ptr<pdf> > &Sources ) {
+	pdfs = Sources;
+	Nsources = pdfs.length();
+	//set sizes
+	dls.set_size ( Sources.length() );
+	rvzs.set_size ( Sources.length() );
+	zdls.set_size ( Sources.length() );
+
+	rv = get_composite_rv ( pdfs, /* checkoverlap = */ false );
+
+	RV rvc;
+	// Extend rv by rvc!
+	for ( int i = 0; i < pdfs.length(); i++ ) {
+		RV rvx = pdfs ( i )->_rvc().subt ( rv );
+		rvc.add ( rvx ); // add rv to common rvc
+	}
+
+	// join rv and rvc - see descriprion
+	rv.add ( rvc );
+	// get dimension
+	dim = rv._dsize();
+
+	// create links between sources and common rv
+	RV xytmp;
+	for ( int i = 0; i < pdfs.length(); i++ ) {
+		//Establich connection between pdfs and merger
+		dls ( i ) = new datalink_m2e;
+		dls ( i )->set_connection ( pdfs ( i )->_rv(), pdfs ( i )->_rvc(), rv );
+
+		// find out what is missing in each pdf
+		xytmp = pdfs ( i )->_rv();
+		xytmp.add ( pdfs ( i )->_rvc() );
+		// z_i = common_rv-xy
+		rvzs ( i ) = rv.subt ( xytmp );
+		//establish connection between extension (z_i|x,y)s and common rv
+		zdls ( i ) = new datalink_m2e;
+		zdls ( i )->set_connection ( rvzs ( i ), xytmp, rv ) ;
+	};
+}
+
+void merger_base::set_support ( rectangular_support &Sup ) {
+	Npoints = Sup.points();
+	eSmp.set_parameters ( Npoints, false );
+	Array<vec> &samples = eSmp._samples();
+	eSmp._w() = ones ( Npoints ) / Npoints; //unifrom size of bins
+	//set samples
+	samples ( 0 ) = Sup.first_vec();
+	for ( int j = 1; j < Npoints; j++ ) {
+		samples ( j ) = Sup.next_vec();
+	}
+}
+
+void merger_base::merge () {
+	validate();
+
+	//check if sources overlap:
+	bool OK = true;
+	for ( int i = 0; i < pdfs.length(); i++ ) {
+		OK &= ( rvzs ( i )._dsize() == 0 ); // z_i is empty
+		OK &= ( pdfs ( i )->_rvc()._dsize() == 0 ); // y_i is empty
+	}
+
+	if ( OK ) {
+		mat lW = zeros ( pdfs.length(), eSmp._w().length() );
+
+		vec emptyvec ( 0 );
+		for ( int i = 0; i < pdfs.length(); i++ ) {
+			for ( int j = 0; j < eSmp._w().length(); j++ ) {
+				lW ( i, j ) = pdfs ( i )->evallogcond ( eSmp._samples() ( j ), emptyvec );
+			}
+		}
+
+		vec w_nn = merge_points ( lW );
+		vec wtmp = exp ( w_nn - max ( w_nn ) );
+		//renormalize
+		eSmp._w() = wtmp / sum ( wtmp );
+	} else {
+		bdm_error ( "Sources are not compatible - use merger_mix" );
+	}
 }
 
@@ -62,4 +143,40 @@
 }
 
+vec merger_base::mean() const {
+	const Vec<double> &w = eSmp._w();
+	const Array<vec> &S = eSmp._samples();
+	vec tmp = zeros ( dim );
+	for ( int i = 0; i < Npoints; i++ ) {
+		tmp += w ( i ) * S ( i );
+	}
+	return tmp;
+}
+
+mat merger_base::covariance() const {
+	const vec &w = eSmp._w();
+	const Array<vec> &S = eSmp._samples();
+
+	vec mea = mean();
+
+// 			cout << sum (w) << "," << w*w << endl;
+
+	mat Tmp = zeros ( dim, dim );
+	for ( int i = 0; i < Npoints; i++ ) {
+		Tmp += w ( i ) * outer_product ( S ( i ), S ( i ) );
+	}
+	return Tmp - outer_product ( mea, mea );
+}
+
+vec merger_base::variance() const {
+	const vec &w = eSmp._w();
+	const Array<vec> &S = eSmp._samples();
+
+	vec tmp = zeros ( dim );
+	for ( int i = 0; i < Nsources; i++ ) {
+		tmp += w ( i ) * pow ( S ( i ), 2 );
+	}
+	return tmp - pow ( mean(), 2 );
+}
+
 void merger_mix::merge ( ) {
 	Array<vec> &Smp = eSmp._samples(); //aux
Index: library/bdm/stat/emix.cpp
===================================================================
--- library/bdm/stat/emix.cpp (revision 737)
+++ library/bdm/stat/emix.cpp (revision 739)
@@ -34,4 +34,53 @@
 
 	return Coms ( i )->sample();
+}
+
+vec emix::mean() const {
+	int i;
+	vec mu = zeros ( dim );
+	for ( i = 0; i < w.length(); i++ ) {
+		mu += w ( i ) * Coms ( i )->mean();
+	}
+	return mu;
+}
+
+vec emix::variance() const {
+	//non-central moment
+	vec mom2 = zeros ( dim );
+	for ( int i = 0; i < w.length(); i++ ) {
+		mom2 += w ( i ) * ( Coms ( i )->variance() + pow ( Coms ( i )->mean(), 2 ) );
+	}
+	//central moment
+	return mom2 - pow ( mean(), 2 );
+}
+
+double emix::evallog ( const vec &val ) const {
+	int i;
+	double sum = 0.0;
+	for ( i = 0; i < w.length(); i++ ) {
+		sum += w ( i ) * exp ( Coms ( i )->evallog ( val ) );
+	}
+	if ( sum == 0.0 ) {
+		sum = std::numeric_limits<double>::epsilon();
+	}
+	double tmp = log ( sum );
+	bdm_assert_debug ( std::isfinite ( tmp ), "Infinite" );
+	return tmp;
+}
+
+vec emix::evallog_mat ( const mat &Val ) const {
+	vec x = zeros ( Val.cols() );
+	for ( int i = 0; i < w.length(); i++ ) {
+		x += w ( i ) * exp ( Coms ( i )->evallog_mat ( Val ) );
+	}
+	return log ( x );
+};
+
+mat emix::evallog_coms ( const mat &Val ) const {
+	mat X ( w.length(), Val.cols() );
+	for ( int i = 0; i < w.length(); i++ ) {
+		X.set_row ( i, w ( i ) *exp ( Coms ( i )->evallog_mat ( Val ) ) );
+	}
+	return X;
 }
 
@@ -203,4 +252,37 @@
 };
 
+double mprod::evallogcond ( const vec &val, const vec &cond ) {
+	int i;
+	double res = 0.0;
+	for ( i = pdfs.length() - 1; i >= 0; i-- ) {
+		/*			if ( pdfs(i)->_rvc().count() >0) {
+						pdfs ( i )->condition ( dls ( i )->get_cond ( val,cond ) );
+					}
+					// add logarithms
+					res += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );*/
+		res += pdfs ( i )->evallogcond (
+		           dls ( i )->pushdown ( val ),
+		           dls ( i )->get_cond ( val, cond )
+		       );
+	}
+	return res;
+}
+
+vec mprod::evallogcond_mat ( const mat &Dt, const vec &cond ) {
+	vec tmp ( Dt.cols() );
+	for ( int i = 0; i < Dt.cols(); i++ ) {
+		tmp ( i ) = evallogcond ( Dt.get_col ( i ), cond );
+	}
+	return tmp;
+}
+
+vec mprod::evallogcond_mat ( const Array<vec> &Dt, const vec &cond ) {
+	vec tmp ( Dt.length() );
+	for ( int i = 0; i < Dt.length(); i++ ) {
+		tmp ( i ) = evallogcond ( Dt ( i ), cond );
+	}
+	return tmp;
+}
+
 void mprod::set_elements ( const Array<shared_ptr<pdf> > &mFacs ) {
 	pdfs = mFacs;
@@ -236,4 +318,38 @@
 
 	return Coms ( i )->samplecond ( cond );
+}
+
+vec eprod::mean() const {
+	vec tmp ( dim );
+	for ( int i = 0; i < epdfs.length(); i++ ) {
+		vec pom = epdfs ( i )->mean();
+		dls ( i )->pushup ( tmp, pom );
+	}
+	return tmp;
+}
+
+vec eprod::variance() const {
+	vec tmp ( dim ); //second moment
+	for ( int i = 0; i < epdfs.length(); i++ ) {
+		vec pom = epdfs ( i )->mean();
+		dls ( i )->pushup ( tmp, pow ( pom, 2 ) );
+	}
+	return tmp - pow ( mean(), 2 );
+}
+vec eprod::sample() const {
+	vec tmp ( dim );
+	for ( int i = 0; i < epdfs.length(); i++ ) {
+		vec pom = epdfs ( i )->sample();
+		dls ( i )->pushup ( tmp, pom );
+	}
+	return tmp;
+}
+double eprod::evallog ( const vec &val ) const {
+	double tmp = 0;
+	for ( int i = 0; i < epdfs.length(); i++ ) {
+		tmp += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );
+	}
+	bdm_assert_debug ( std::isfinite ( tmp ), "Infinite" );
+	return tmp;
 }
 
@@ -267,2 +383,3 @@
 // 		}
 // 	};
+
Index: library/bdm/stat/exp_family.h
===================================================================
--- library/bdm/stat/exp_family.h (revision 737)
+++ library/bdm/stat/exp_family.h (revision 739)
@@ -323,28 +323,7 @@
 		// 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_vector ( RV ( "", th_dim ), prefix + logrec->L.prefix_sep() + "mean" );
-			logrec->ids ( 1 ) = L.add_vector ( 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.log_vector ( logrec->ids ( 0 ), est_theta() );
-			logrec->L.log_vector ( logrec->ids ( 1 ), cvectorize ( est_theta_cov().to_mat() ) );
-		} else {
-			epdf::log_write();
-		}
-
-	}
+	void log_register ( bdm::logger& L, const string& prefix );
+
+	void log_write() const;
 	//!@}
 };
@@ -526,40 +505,10 @@
 		beta = mB->beta;
 	}
-	void bayes ( const vec &yt, const vec &cond = empty_vec ) {
-		if ( frg < 1.0 ) {
-			beta *= frg;
-			last_lognc = est.lognc();
-		}
-		beta += yt;
-		if ( evalll ) {
-			ll = est.lognc() - last_lognc;
-		}
-	}
-	double logpred ( const vec &yt ) const {
-		eDirich pred ( est );
-		vec &beta = pred._beta();
-
-		double lll;
-		if ( frg < 1.0 ) {
-			beta *= frg;
-			lll = pred.lognc();
-		} else if ( evalll ) {
-			lll = last_lognc;
-		} else {
-			lll = pred.lognc();
-		}
-
-		beta += yt;
-		return pred.lognc() - lll;
-	}
-	void flatten ( const BMEF* B ) {
-		const multiBM* E = dynamic_cast<const multiBM*> ( B );
-		// sum(beta) should be equal to sum(B.beta)
-		const vec &Eb = E->beta;//const_cast<multiBM*> ( E )->_beta();
-		beta *= ( sum ( Eb ) / sum ( beta ) );
-		if ( evalll ) {
-			last_lognc = est.lognc();
-		}
-	}
+	void bayes ( const vec &yt, const vec &cond = empty_vec );
+
+	double logpred ( const vec &yt ) const;
+
+	void flatten ( const BMEF* B );
+
 	//! return correctly typed posterior (covariant return)
 	const eDirich& posterior() const {
@@ -971,20 +920,6 @@
 		Lambda = Lambda0;
 	}
-	void condition ( const vec &cond ) {
-		if ( cond.length() > 0 ) {
-			iepdf._mu() = A * cond + mu_const;
-		} else {
-			iepdf._mu() =  mu_const;
-		}
-		double zeta;
-		//ugly hack!
-		if ( ( cond.length() + 1 ) == Lambda.rows() ) {
-			zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) );
-		} else {
-			zeta = Lambda.invqform ( cond );
-		}
-		_R = Re;
-		_R *= ( 1 + zeta );// / ( nu ); << nu is in Re!!!!!!
-	};
+
+	void condition ( const vec &cond ); 
 
 	void validate() {
@@ -1526,22 +1461,5 @@
 	}
 	//! For this class, qbounds are minimum and maximum value of the population!
-	void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const {
-		// lb in inf so than it will be pushed below;
-		lb.set_size ( dim );
-		ub.set_size ( dim );
-		lb = std::numeric_limits<double>::infinity();
-		ub = -std::numeric_limits<double>::infinity();
-		int j;
-		for ( int i = 0; i < n; i++ ) {
-			for ( j = 0; j < dim; j++ ) {
-				if ( samples ( i ) ( j ) < lb ( j ) ) {
-					lb ( j ) = samples ( i ) ( j );
-				}
-				if ( samples ( i ) ( j ) > ub ( j ) ) {
-					ub ( j ) = samples ( i ) ( j );
-				}
-			}
-		}
-	}
+	void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const;
 };
 
Index: library/bdm/stat/merger.h
===================================================================
--- library/bdm/stat/merger.h (revision 737)
+++ library/bdm/stat/merger.h (revision 739)
@@ -93,55 +93,9 @@
 
 	//! Function setting the main internal structures
-	void set_sources ( const Array<shared_ptr<pdf> > &Sources ) {
-		pdfs = Sources;
-		Nsources = pdfs.length();
-		//set sizes
-		dls.set_size ( Sources.length() );
-		rvzs.set_size ( Sources.length() );
-		zdls.set_size ( Sources.length() );
-
-		rv = get_composite_rv ( pdfs, /* checkoverlap = */ false );
-
-		RV rvc;
-		// Extend rv by rvc!
-		for ( int i = 0; i < pdfs.length(); i++ ) {
-			RV rvx = pdfs ( i )->_rvc().subt ( rv );
-			rvc.add ( rvx ); // add rv to common rvc
-		}
-
-		// join rv and rvc - see descriprion
-		rv.add ( rvc );
-		// get dimension
-		dim = rv._dsize();
-
-		// create links between sources and common rv
-		RV xytmp;
-		for ( int i = 0; i < pdfs.length(); i++ ) {
-			//Establich connection between pdfs and merger
-			dls ( i ) = new datalink_m2e;
-			dls ( i )->set_connection ( pdfs ( i )->_rv(), pdfs ( i )->_rvc(), rv );
-
-			// find out what is missing in each pdf
-			xytmp = pdfs ( i )->_rv();
-			xytmp.add ( pdfs ( i )->_rvc() );
-			// z_i = common_rv-xy
-			rvzs ( i ) = rv.subt ( xytmp );
-			//establish connection between extension (z_i|x,y)s and common rv
-			zdls ( i ) = new datalink_m2e;
-			zdls ( i )->set_connection ( rvzs ( i ), xytmp, rv ) ;
-		};
-	}
+	void set_sources ( const Array<shared_ptr<pdf> > &Sources );
+
 	//! Set support points from rectangular grid
-	void set_support ( rectangular_support &Sup ) {
-		Npoints = Sup.points();
-		eSmp.set_parameters ( Npoints, false );
-		Array<vec> &samples = eSmp._samples();
-		eSmp._w() = ones ( Npoints ) / Npoints; //unifrom size of bins
-		//set samples
-		samples ( 0 ) = Sup.first_vec();
-		for ( int j = 1; j < Npoints; j++ ) {
-			samples ( j ) = Sup.next_vec();
-		}
-	}
+	void set_support ( rectangular_support &Sup );
+
 	//! Set support points from dicrete grid
 	void set_support ( discrete_support &Sup ) {
@@ -180,33 +134,5 @@
 
 	//!Merge given sources in given points
-	virtual void merge () {
-		validate();
-
-		//check if sources overlap:
-		bool OK = true;
-		for ( int i = 0; i < pdfs.length(); i++ ) {
-			OK &= ( rvzs ( i )._dsize() == 0 ); // z_i is empty
-			OK &= ( pdfs ( i )->_rvc()._dsize() == 0 ); // y_i is empty
-		}
-
-		if ( OK ) {
-			mat lW = zeros ( pdfs.length(), eSmp._w().length() );
-
-			vec emptyvec ( 0 );
-			for ( int i = 0; i < pdfs.length(); i++ ) {
-				for ( int j = 0; j < eSmp._w().length(); j++ ) {
-					lW ( i, j ) = pdfs ( i )->evallogcond ( eSmp._samples() ( j ), emptyvec );
-				}
-			}
-
-			vec w_nn = merge_points ( lW );
-			vec wtmp = exp ( w_nn - max ( w_nn ) );
-			//renormalize
-			eSmp._w() = wtmp / sum ( wtmp );
-		} else {
-			bdm_error ( "Sources are not compatible - use merger_mix" );
-		}
-	};
-
+	virtual void merge ();
 
 	//! Merge log-likelihood values in points using method specified by parameter METHOD
@@ -216,37 +142,10 @@
 	//! sample from merged density
 //! weight w is a
-	vec mean() const {
-		const Vec<double> &w = eSmp._w();
-		const Array<vec> &S = eSmp._samples();
-		vec tmp = zeros ( dim );
-		for ( int i = 0; i < Npoints; i++ ) {
-			tmp += w ( i ) * S ( i );
-		}
-		return tmp;
-	}
-	mat covariance() const {
-		const vec &w = eSmp._w();
-		const Array<vec> &S = eSmp._samples();
-
-		vec mea = mean();
-
-// 			cout << sum (w) << "," << w*w << endl;
-
-		mat Tmp = zeros ( dim, dim );
-		for ( int i = 0; i < Npoints; i++ ) {
-			Tmp += w ( i ) * outer_product ( S ( i ), S ( i ) );
-		}
-		return Tmp - outer_product ( mea, mea );
-	}
-	vec variance() const {
-		const vec &w = eSmp._w();
-		const Array<vec> &S = eSmp._samples();
-
-		vec tmp = zeros ( dim );
-		for ( int i = 0; i < Nsources; i++ ) {
-			tmp += w ( i ) * pow ( S ( i ), 2 );
-		}
-		return tmp - pow ( mean(), 2 );
-	}
+	vec mean() const;
+
+	mat covariance() const;
+
+	vec variance() const;
+
 	//!@}
 
Index: library/bdm/stat/emix.h
===================================================================
--- library/bdm/stat/emix.h (revision 737)
+++ library/bdm/stat/emix.h (revision 739)
@@ -129,51 +129,15 @@
 
 	vec sample() const;
-	vec mean() const {
-		int i;
-		vec mu = zeros ( dim );
-		for ( i = 0; i < w.length(); i++ ) {
-			mu += w ( i ) * Coms ( i )->mean();
-		}
-		return mu;
-	}
-	vec variance() const {
-		//non-central moment
-		vec mom2 = zeros ( dim );
-		for ( int i = 0; i < w.length(); i++ ) {
-			mom2 += w ( i ) * ( Coms ( i )->variance() + pow ( Coms ( i )->mean(), 2 ) );
-		}
-		//central moment
-		return mom2 - pow ( mean(), 2 );
-	}
-	double evallog ( const vec &val ) const {
-		int i;
-		double sum = 0.0;
-		for ( i = 0; i < w.length(); i++ ) {
-			sum += w ( i ) * exp ( Coms ( i )->evallog ( val ) );
-		}
-		if ( sum == 0.0 ) {
-			sum = std::numeric_limits<double>::epsilon();
-		}
-		double tmp = log ( sum );
-		bdm_assert_debug ( std::isfinite ( tmp ), "Infinite" );
-		return tmp;
-	};
-
-	vec evallog_mat ( const mat &Val ) const {
-		vec x = zeros ( Val.cols() );
-		for ( int i = 0; i < w.length(); i++ ) {
-			x += w ( i ) * exp ( Coms ( i )->evallog_mat ( Val ) );
-		}
-		return log ( x );
-	};
+
+	vec mean() const;
+
+	vec variance() const;
+
+	double evallog ( const vec &val ) const;
+
+	vec evallog_mat ( const mat &Val ) const;
 
 	//! Auxiliary function that returns pdflog for each component
-	mat evallog_coms ( const mat &Val ) const {
-		mat X ( w.length(), Val.cols() );
-		for ( int i = 0; i < w.length(); i++ ) {
-			X.set_row ( i, w ( i ) *exp ( Coms ( i )->evallog_mat ( Val ) ) );
-		}
-		return X;
-	};
+	mat evallog_coms ( const mat &Val ) const;
 
 	shared_ptr<epdf> marginal ( const RV &rv ) const;
@@ -334,35 +298,9 @@
 	void set_elements ( const Array<shared_ptr<pdf> > &mFacs );
 
-	double evallogcond ( const vec &val, const vec &cond ) {
-		int i;
-		double res = 0.0;
-		for ( i = pdfs.length() - 1; i >= 0; i-- ) {
-			/*			if ( pdfs(i)->_rvc().count() >0) {
-							pdfs ( i )->condition ( dls ( i )->get_cond ( val,cond ) );
-						}
-						// add logarithms
-						res += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );*/
-			res += pdfs ( i )->evallogcond (
-			           dls ( i )->pushdown ( val ),
-			           dls ( i )->get_cond ( val, cond )
-			       );
-		}
-		return res;
-	}
-	vec evallogcond_mat ( const mat &Dt, const vec &cond ) {
-		vec tmp ( Dt.cols() );
-		for ( int i = 0; i < Dt.cols(); i++ ) {
-			tmp ( i ) = evallogcond ( Dt.get_col ( i ), cond );
-		}
-		return tmp;
-	};
-	vec evallogcond_mat ( const Array<vec> &Dt, const vec &cond ) {
-		vec tmp ( Dt.length() );
-		for ( int i = 0; i < Dt.length(); i++ ) {
-			tmp ( i ) = evallogcond ( Dt ( i ), cond );
-		}
-		return tmp;
-	};
-
+	double evallogcond ( const vec &val, const vec &cond );
+
+	vec evallogcond_mat ( const mat &Dt, const vec &cond );
+
+	vec evallogcond_mat ( const Array<vec> &Dt, const vec &cond );
 
 	//TODO smarter...
@@ -441,36 +379,12 @@
 	}
 
-	vec mean() const {
-		vec tmp ( dim );
-		for ( int i = 0; i < epdfs.length(); i++ ) {
-			vec pom = epdfs ( i )->mean();
-			dls ( i )->pushup ( tmp, pom );
-		}
-		return tmp;
-	}
-	vec variance() const {
-		vec tmp ( dim ); //second moment
-		for ( int i = 0; i < epdfs.length(); i++ ) {
-			vec pom = epdfs ( i )->mean();
-			dls ( i )->pushup ( tmp, pow ( pom, 2 ) );
-		}
-		return tmp - pow ( mean(), 2 );
-	}
-	vec sample() const {
-		vec tmp ( dim );
-		for ( int i = 0; i < epdfs.length(); i++ ) {
-			vec pom = epdfs ( i )->sample();
-			dls ( i )->pushup ( tmp, pom );
-		}
-		return tmp;
-	}
-	double evallog ( const vec &val ) const {
-		double tmp = 0;
-		for ( int i = 0; i < epdfs.length(); i++ ) {
-			tmp += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );
-		}
-		bdm_assert_debug ( std::isfinite ( tmp ), "Infinite" );
-		return tmp;
-	}
+	vec mean() const;
+
+	vec variance() const;
+
+	vec sample() const;
+
+	double evallog ( const vec &val ) const;
+
 	//!access function
 	const epdf* operator () ( int i ) const {
Index: library/bdm/stat/exp_family.cpp
===================================================================
--- library/bdm/stat/exp_family.cpp (revision 737)
+++ library/bdm/stat/exp_family.cpp (revision 739)
@@ -279,4 +279,68 @@
 	M = iLsub * Lpsi;
 	R = ldR.to_mat()  ;
+}
+
+void egiw::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_vector ( RV ( "", th_dim ), prefix + logrec->L.prefix_sep() + "mean" );
+		logrec->ids ( 1 ) = L.add_vector ( RV ( "", th_dim * th_dim ), prefix + logrec->L.prefix_sep() + "variance" );
+	} else {
+		epdf::log_register ( L, prefix );
+	}
+}
+
+void egiw::log_write() const {
+	if ( log_level == 3 ) {
+		mat M;
+		ldmat Lam;
+		ldmat Vz;
+		factorize ( M, Vz, Lam );
+		logrec->L.log_vector ( logrec->ids ( 0 ), est_theta() );
+		logrec->L.log_vector ( logrec->ids ( 1 ), cvectorize ( est_theta_cov().to_mat() ) );
+	} else {
+		epdf::log_write();
+	}
+
+}
+
+void multiBM::bayes ( const vec &yt, const vec &cond ) {
+	if ( frg < 1.0 ) {
+		beta *= frg;
+		last_lognc = est.lognc();
+	}
+	beta += yt;
+	if ( evalll ) {
+		ll = est.lognc() - last_lognc;
+	}
+}
+
+double multiBM::logpred ( const vec &yt ) const {
+	eDirich pred ( est );
+	vec &beta = pred._beta();
+
+	double lll;
+	if ( frg < 1.0 ) {
+		beta *= frg;
+		lll = pred.lognc();
+	} else if ( evalll ) {
+		lll = last_lognc;
+	} else {
+		lll = pred.lognc();
+	}
+
+	beta += yt;
+	return pred.lognc() - lll;
+}
+void multiBM::flatten ( const BMEF* B ) {
+	const multiBM* E = dynamic_cast<const multiBM*> ( B );
+	// sum(beta) should be equal to sum(B.beta)
+	const vec &Eb = E->beta;//const_cast<multiBM*> ( E )->_beta();
+	beta *= ( sum ( Eb ) / sum ( beta ) );
+	if ( evalll ) {
+		last_lognc = est.lognc();
+	}
 }
 
@@ -462,3 +526,40 @@
 }
 
+void mlstudent::condition ( const vec &cond ) {
+	if ( cond.length() > 0 ) {
+		iepdf._mu() = A * cond + mu_const;
+	} else {
+		iepdf._mu() =  mu_const;
+	}
+	double zeta;
+	//ugly hack!
+	if ( ( cond.length() + 1 ) == Lambda.rows() ) {
+		zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) );
+	} else {
+		zeta = Lambda.invqform ( cond );
+	}
+	_R = Re;
+	_R *= ( 1 + zeta );// / ( nu ); << nu is in Re!!!!!!
+}
+
+void eEmp::qbounds ( vec &lb, vec &ub, double perc ) const {
+	// lb in inf so than it will be pushed below;
+	lb.set_size ( dim );
+	ub.set_size ( dim );
+	lb = std::numeric_limits<double>::infinity();
+	ub = -std::numeric_limits<double>::infinity();
+	int j;
+	for ( int i = 0; i < n; i++ ) {
+		for ( j = 0; j < dim; j++ ) {
+			if ( samples ( i ) ( j ) < lb ( j ) ) {
+				lb ( j ) = samples ( i ) ( j );
+			}
+			if ( samples ( i ) ( j ) > ub ( j ) ) {
+				ub ( j ) = samples ( i ) ( j );
+			}
+		}
+	}
+}
+
+
 };
