Index: bdm/stat/libEF.cpp
===================================================================
--- bdm/stat/libEF.cpp (revision 170)
+++ bdm/stat/libEF.cpp (revision 178)
@@ -49,4 +49,5 @@
 #define	logpi 1.144729885849400163877476189
 #define log2pi 1.83787706640935
+#define Inf std::numeric_limits<double>::infinity()
 
 	double nkG = 0.5* xdim* ( -nPsi *log2pi + sum ( log ( D ( xdim,D.length()-1 ) ) ) );
@@ -58,4 +59,5 @@
 	             - 0.5*xdim* ( m*log2 + 0.5* ( xdim-1 ) *log2pi )  - lg;
 
+	it_assert_debug(((-nkG-nkW)>-Inf) && ((-nkG-nkW)<Inf), "ARX improper");
 	return -nkG-nkW;
 }
@@ -244,5 +246,5 @@
 }
 
-void eEmp::set_parameters ( const vec &w0, epdf* epdf0 ) {
+void eEmp::set_parameters ( const vec &w0, const epdf* epdf0 ) {
 	//it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
 	w=w0;
@@ -253,2 +255,11 @@
 	for ( int i=0;i<n;i++ ) {samples ( i ) =epdf0->sample();}
 }
+
+void eEmp::set_samples (  const epdf* epdf0 ) {
+	//it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
+	w=1;
+	w/=sum ( w );//renormalize
+
+	for ( int i=0;i<n;i++ ) {samples ( i ) =epdf0->sample();}
+}
+
Index: bdm/stat/emix.cpp
===================================================================
--- bdm/stat/emix.cpp (revision 175)
+++ bdm/stat/emix.cpp (revision 178)
@@ -3,11 +3,19 @@
 using namespace itpp;
 
-void emix::set_parameters ( const vec &w0, const Array<epdf*> &Coms0 ) {
-	w = w0/sum(w0);
+void emix::set_parameters ( const vec &w0, const Array<epdf*> &Coms0, bool copy ) {
+	w = w0/sum ( w0 );
 	int i;
 	for ( i=0;i<w.length();i++ ) {
 		it_assert_debug ( rv.equal ( Coms0 ( i )->_rv() ),"RVs do not match!" );
 	}
-	Coms = Coms0;
+	if ( copy ) {
+		Coms.set_length(Coms0.length());
+		for ( i=0;i<w.length();i++ ) {*Coms ( i ) =*Coms0 ( i );}
+		destroyComs=true;
+	}
+	else {
+		Coms = Coms0;
+		destroyComs=false;
+	}
 }
 
@@ -37,5 +45,5 @@
 // 			rvc.add ( mpdfs ( i )->_rvc().subt ( rv ) ); //add rv to common rvs.
 // 		};
-// 		
+//
 // //		independent = true;
 // 		//test rvc of mpdfs and fill rvinds
Index: bdm/stat/libBM.cpp
===================================================================
--- bdm/stat/libBM.cpp (revision 176)
+++ bdm/stat/libBM.cpp (revision 178)
@@ -157,5 +157,6 @@
 RV RV::subt ( const RV &rv2 ) const {
 	ivec res = this->findself ( rv2 ); // nonzeros
-	ivec valid = itpp::find ( res == -1 ); //-1 => value not found => it remains
+	ivec valid;
+	if (tsize>0) {valid= itpp::find ( res == -1 );} //-1 => value not found => it remains
 	return ( *this ) ( valid ); //keep those that were not found in rv2
 }
@@ -196,5 +197,6 @@
 void compositepdf::setrvc(const RV &rv, RV &rvc){
 	for (int i = 0;i < n;i++ ) {
-		rvc.add ( mpdfs ( i )->_rvc().subt ( rv ) ); //add rv to common rvc
+		RV rvx = mpdfs ( i )->_rvc().subt ( rv );
+		rvc.add ( rvx ); //add rv to common rvc
 	};
 }
@@ -202,7 +204,10 @@
 void compositepdf::setindices(const RV &rv){
 	for (int i = 0;i < n;i++ ) {
-			// find ith rv in common rv
+		// find ith rv in common rv
 		rvsinrv ( i ) = mpdfs ( i )->_rv().dataind ( rv );
 		rvcsinrv ( i ) = mpdfs ( i )->_rvc().dataind ( rv );
+		rvinrvcs ( i ) = rv.dataind(mpdfs ( i )->_rvc());
+		//rvcsinrv and rvinrvcs must be 1-to-1
+		it_assert_debug(rvcsinrv ( i ).length()==rvinrvcs(i).length(),"" );
 	}
 }
Index: bdm/stat/libEF.h
===================================================================
--- bdm/stat/libEF.h (revision 176)
+++ bdm/stat/libEF.h (revision 178)
@@ -44,7 +44,7 @@
 	virtual double lognc() const =0;
 	//!TODO decide if it is really needed
-	virtual void dupdate ( mat &v ) {it_error ( "Not implemneted" );};
+	virtual void dupdate ( mat &v ) {it_error ( "Not implemented" );};
 	//!Evaluate normalized log-probability
-	virtual double evalpdflog_nn ( const vec &val ) const{it_error ( "Not implemneted" );return 0.0;};
+	virtual double evalpdflog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
 	//!Evaluate normalized log-probability
 	virtual double evalpdflog ( const vec &val ) const {return evalpdflog_nn ( val )-lognc();}
@@ -90,7 +90,12 @@
 	//original Bayes
 	void bayes ( const vec &dt );
-	//!Flatten the posterior
-	virtual void flatten ( BMEF * B ) {it_error ( "Not implemented" );}
-};
+	//!Flatten the posterior according to the given BMEF (of the same type!)
+	virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
+	//!Flatten the posterior as if to keep nu0 data
+	virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
+};
+
+template<class sq_T>
+class mlnorm;
 
 /*!
@@ -100,5 +105,4 @@
 */
 template<class sq_T>
-
 class enorm : public eEF {
 protected:
@@ -110,11 +114,10 @@
 	int dim;
 public:
-//	enorm() :eEF() {};
 	//!Default constructor
-	enorm ( RV &rv );
+	enorm ( const RV &rv );
 	//! Set mean value \c mu and covariance \c R
 	void set_parameters ( const vec &mu,const sq_T &R );
-	//! tupdate in exponential form (not really handy)
-	void tupdate ( double phi, mat &vbar, double nubar );
+	////! tupdate in exponential form (not really handy)
+	//void tupdate ( double phi, mat &vbar, double nubar );
 	//! dupdate in exponential form (not really handy)
 	void dupdate ( mat &v,double nu=1.0 );
@@ -124,8 +127,9 @@
 	mat sample ( int N ) const;
 	double eval ( const vec &val ) const ;
-	double evalpdflog ( const vec &val ) const;
+	double evalpdflog_nn ( const vec &val ) const;
 	double lognc () const;
 	vec mean() const {return mu;}
-
+	mlnorm<sq_T>* condition ( const RV &rvn );
+	enorm<sq_T>* marginal ( const RV &rv );
 //Access methods
 	//! returns a pointer to the internal mean value. Use with Care!
@@ -139,5 +143,5 @@
 
 	//! access method
-	mat getR () {return R.to_mat();}
+//	mat getR () {return R.to_mat();}
 };
 
@@ -215,10 +219,10 @@
 	};
 	//!access function
-	vec& _beta() {return beta;}
+	vec& _beta()  {return beta;}
 	//!Set internal parameters
-	void set_parameters(const vec &beta0){
-		if(beta0.length()!=beta.length()){
-			it_assert_debug(rv.length()==1,"Undefined");
-			rv.set_size(0,beta0.length());
+	void set_parameters ( const vec &beta0 ) {
+		if ( beta0.length() !=beta.length() ) {
+			it_assert_debug ( rv.length() ==1,"Undefined" );
+			rv.set_size ( 0,beta0.length() );
 		}
 		beta= beta0;
@@ -258,8 +262,8 @@
 		return pred.lognc()-lll;
 	}
-	void flatten ( BMEF* B ) {
-		eDirich* E=dynamic_cast<eDirich*> ( B );
+	void flatten ( const BMEF* B ) {
+		const eDirich* E=dynamic_cast<const eDirich*> ( B );
 		// sum(beta) should be equal to sum(B.beta)
-		const vec &Eb=E->_beta();
+		const vec &Eb=const_cast<eDirich*> ( E )->_beta();
 		est.pow ( sum ( beta ) /sum ( Eb ) );
 		if ( evalll ) {last_lognc=est.lognc();}
@@ -267,7 +271,7 @@
 	const epdf& _epdf() const {return est;};
 	void set_parameters ( const vec &beta0 ) {
-		est.set_parameters(beta0);
+		est.set_parameters ( beta0 );
 		rv = est._rv();
-		if(evalll){last_lognc=est.lognc();}
+		if ( evalll ) {last_lognc=est.lognc();}
 	}
 };
@@ -359,5 +363,5 @@
  \brief Normal distributed linear function with linear function of mean value;
 
- Mean value \f$mu=A*rvc\f$.
+ Mean value \f$mu=A*rvc+mu_0\f$.
 */
 template<class sq_T>
@@ -366,16 +370,17 @@
 	enorm<sq_T> epdf;
 	mat A;
+	vec mu_const;
 	vec& _mu; //cached epdf.mu;
 public:
 	//! Constructor
-	mlnorm ( RV &rv,RV &rvc );
+	mlnorm (const RV &rv, const RV &rvc );
 	//! Set \c A and \c R
-	void set_parameters ( const  mat &A, const sq_T &R );
+	void set_parameters ( const  mat &A, const vec &mu0, const sq_T &R );
 	//!Generate one sample of the posterior
-	vec samplecond ( vec &cond, double &lik );
+	vec samplecond (const vec &cond, double &lik );
 	//!Generate matrix of samples of the posterior
-	mat samplecond ( vec &cond, vec &lik, int n );
+	mat samplecond (const vec &cond, vec &lik, int n );
 	//! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
-	void condition ( vec &cond );
+	void condition (const vec &cond );
 };
 
@@ -453,6 +458,8 @@
 	//! Default constructor
 	eEmp ( const RV &rv0 ,int n0 ) :epdf ( rv0 ),n ( n0 ),w ( n ),samples ( n ) {};
+	//! Set samples and weights
+	void set_parameters ( const vec &w0, const epdf* pdf0 );
 	//! Set sample
-	void set_parameters ( const vec &w0, epdf* pdf0 );
+	void set_samples ( const epdf* pdf0 );
 	//! Potentially dangerous, use with care.
 	vec& _w()  {return w;};
@@ -476,5 +483,5 @@
 
 template<class sq_T>
-enorm<sq_T>::enorm ( RV &rv ) :eEF ( rv ), mu ( rv.count() ),R ( rv.count() ),dim ( rv.count() ) {};
+enorm<sq_T>::enorm ( const RV &rv ) :eEF ( rv ), mu ( rv.count() ),R ( rv.count() ),dim ( rv.count() ) {};
 
 template<class sq_T>
@@ -490,8 +497,8 @@
 };
 
-template<class sq_T>
-void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
-	//
-};
+// template<class sq_T>
+// void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
+// 	//
+// };
 
 template<class sq_T>
@@ -531,7 +538,7 @@
 
 template<class sq_T>
-double enorm<sq_T>::evalpdflog ( const vec &val ) const {
+double enorm<sq_T>::evalpdflog_nn ( const vec &val ) const {
 	// 1.83787706640935 = log(2pi)
-	return  -0.5* ( +R.invqform ( mu-val ) ) - lognc();
+	return  -0.5* ( R.invqform ( mu-val ) );// - lognc();
 };
 
@@ -539,20 +546,21 @@
 inline double enorm<sq_T>::lognc () const {
 	// 1.83787706640935 = log(2pi)
-	return -0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
-};
-
-template<class sq_T>
-mlnorm<sq_T>::mlnorm ( RV &rv0,RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu ( epdf._mu() ) {
+	return 0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
+};
+
+template<class sq_T>
+mlnorm<sq_T>::mlnorm ( const RV &rv0, const RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu ( epdf._mu() ) {
 	ep =&epdf;
 }
 
 template<class sq_T>
-void mlnorm<sq_T>::set_parameters ( const mat &A0, const sq_T &R0 ) {
+void mlnorm<sq_T>::set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 ) {
 	epdf.set_parameters ( zeros ( rv.count() ),R0 );
 	A = A0;
+	mu_const = mu0;
 }
 
 template<class sq_T>
-vec mlnorm<sq_T>::samplecond ( vec &cond, double &lik ) {
+vec mlnorm<sq_T>::samplecond (const  vec &cond, double &lik ) {
 	this->condition ( cond );
 	vec smp = epdf.sample();
@@ -562,5 +570,5 @@
 
 template<class sq_T>
-mat mlnorm<sq_T>::samplecond ( vec &cond, vec &lik, int n ) {
+mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) {
 	int i;
 	int dim = rv.count();
@@ -579,9 +587,50 @@
 
 template<class sq_T>
-void mlnorm<sq_T>::condition ( vec &cond ) {
-	_mu = A*cond;
+void mlnorm<sq_T>::condition (const vec &cond ) {
+	_mu = A*cond + mu_const;
 //R is already assigned;
 }
 
+template<class sq_T>
+enorm<sq_T>* enorm<sq_T>::marginal ( const RV &rvn ) {
+	ivec irvn = rvn.dataind ( rv );
+
+	sq_T Rn ( R,irvn );
+	enorm<sq_T>* tmp = new enorm<sq_T>( rvn );
+	tmp->set_parameters ( mu ( irvn ), Rn );
+	return tmp;
+}
+
+template<class sq_T>
+mlnorm<sq_T>* enorm<sq_T>::condition ( const RV &rvn ) {
+
+	RV rvc = rv.subt ( rvn );
+	it_assert_debug ( ( rvc.count() +rvn.count() ==rv.count() ),"wrong rvn" );
+	//Permutation vector of the new R
+	ivec irvn = rvn.dataind ( rv );
+	ivec irvc = rvc.dataind ( rv );
+	ivec perm=concat ( irvn , irvc );
+	sq_T Rn ( R,perm );
+
+	//fixme - could this be done in general for all sq_T?
+	mat S=R.to_mat();
+	//fixme
+	int n=rvn.count()-1;
+	int end=R.rows()-1;
+	mat S11 = S.get ( 0,n, 0, n );
+	mat S12 = S.get ( rvn.count(), end, 0, n );
+	mat S22 = S.get ( rvn.count(), end, rvn.count(), end );
+
+	vec mu1 = mu ( irvn );
+	vec mu2 = mu ( irvc );
+	mat A=S12*inv ( S22 );
+	sq_T R_n ( S11 - A *S12.T() );
+
+	mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( rvn,rvc );
+
+	tmp->set_parameters ( A,mu1-A*mu2,R_n );
+	return tmp;
+}
+
 ///////////
 
Index: bdm/stat/emix.h
===================================================================
--- bdm/stat/emix.h (revision 176)
+++ bdm/stat/emix.h (revision 178)
@@ -36,9 +36,11 @@
 	//! Component (epdfs)
 	Array<epdf*> Coms;
+	//!Flag if owning Coms
+	bool destroyComs;
 public:
 	//!Default constructor
-	emix ( RV &rv ) : epdf ( rv ) {};
-	//! Set weights \c w and components \c R
-	void set_parameters ( const vec &w, const Array<epdf*> &Coms );
+	emix (const RV &rv ) : epdf ( rv ) {};
+	//! Set weights \c w and components \c Coms , Coms are not copied!
+	void set_parameters ( const vec &w, const Array<epdf*> &Coms, bool copy=true );
 
 	vec sample() const;
@@ -58,4 +60,7 @@
 	//! returns a pointer to the internal mean value. Use with Care!
 	vec& _w() {return w;}
+	virtual ~emix(){if (destroyComs){for(int i=0;i<Coms.length();i++){delete Coms(i);}}}
+	//! Auxiliary function for taking ownership of the Coms()
+	void ownComs(){destroyComs=true;}
 };
 
@@ -102,10 +107,10 @@
 		ll = 0;
 		for ( int i = ( n - 1 );i >= 0;i-- ) {
-			if ( rvcinds ( i ).length() > 0 ) {
-				condi = zeros ( rvcsinrv.length() + rvcinds.length() );
+			if (( rvcinds ( i ).length() > 0 )||( rvcsinrv ( i ).length() > 0 )) {
+				condi = zeros ( rvcsinrv(i).length() + rvcinds(i).length() );
 				// copy data in condition
 				set_subvector ( condi,rvcinds ( i ), cond );
 				// copy data in already generated sample
-				set_subvector ( condi,rvcsinrv ( i ), smp );
+				set_subvector ( condi,rvinrvcs ( i ), get_vec(smp,rvcsinrv(i)) );
 
 				mpdfs ( i )->condition ( condi );
Index: bdm/stat/libBM.h
===================================================================
--- bdm/stat/libBM.h (revision 176)
+++ bdm/stat/libBM.h (revision 178)
@@ -140,4 +140,5 @@
 };
 
+class mpdf;
 
 //! Probability density function with numerical statistics, e.g. posterior density.
@@ -173,4 +174,8 @@
 		return x;
 	}
+	//! Return conditional density on the given RV, the remaining rvs will be in conditioning
+	mpdf* condition(const RV &rv){it_warning("Not implemented"); return NULL;}
+	//! Return marginal density on the given RV, the remainig rvs are intergrated out
+	epdf* marginal(const RV &rv){it_warning("Not implemented"); return NULL;}
 
 	//! return expected value
@@ -183,4 +188,5 @@
 	//! modifier function - useful when copying epdfs
 	void _renewrv(const RV &in_rv){rv=in_rv;}
+	//!
 };
 
@@ -213,5 +219,5 @@
 	};
 	//! Update \c ep so that it represents this mpdf conditioned on \c rvc = cond
-	virtual void condition ( const vec &cond ) {};
+	virtual void condition ( const vec &cond ) {it_error("Not implemented");};
 
 	//! Shortcut for conditioning and evaluation of the internal epdf. In some cases,  this operation can be implemented efficiently.
@@ -252,6 +258,8 @@
 		//! Indeces of rvc in common rv
 		Array<ivec> rvcsinrv;
+		//! Indeces of common rv in rvc
+		Array<ivec> rvinrvcs;
 	public:
-		compositepdf(Array<mpdf*> A0): n(A0.length()), mpdfs(A0), rvsinrv(n), rvcsinrv(n){};
+		compositepdf(Array<mpdf*> A0): n(A0.length()), mpdfs(A0), rvsinrv(n), rvcsinrv(n),rvinrvcs(n){};
 		RV getrv(bool checkoverlap=false);
 		void setrvc(const RV &rv, RV &rvc);
@@ -329,4 +337,7 @@
 	vec logpred_m(const mat &dt)const{vec tmp(dt.cols());for(int i=0;i<dt.cols();i++){tmp(i)=logpred(dt.get_col(i));}return tmp;}
 	
+	//!Constructs a predictive density (marginal density on data)
+	virtual epdf* predictor(const RV &rv){it_error("Not implemented");return NULL;};
+	
 	//! Destructor for future use;
 	virtual ~BM() {};
