Index: library/bdm/estim/particles.cpp
===================================================================
--- library/bdm/estim/particles.cpp (revision 850)
+++ library/bdm/estim/particles.cpp (revision 887)
@@ -5,22 +5,4 @@
 using std::endl;
 
-void PF::bayes_gensmp ( const vec &ut ) {
-	if ( ut.length() > 0 ) {
-		vec cond ( par->dimensionc() );
-		cond.set_subvector ( par->dimension(), ut );
-#pragma parallel for
-		for ( int i = 0; i < n; i++ ) {
-			cond.set_subvector ( 0, _samples ( i ) );
-			_samples ( i ) = par->samplecond ( cond );
-			lls ( i ) = 0;
-		}
-	} else {
-#pragma parallel for
-		for ( int i = 0; i < n; i++ ) {
-			_samples ( i ) = par->samplecond ( _samples ( i ) );
-			lls ( i ) = 0;
-		}
-	}
-}
 
 void PF::bayes_weights() {
@@ -29,39 +11,37 @@
 	// compute weights
 	for ( int  i = 0; i < n; i++ ) {
-		_w ( i ) *= exp ( lls ( i ) - mlls ); // multiply w by likelihood
+		w ( i ) *= exp ( lls ( i ) - mlls ); // multiply w by likelihood
 	}
 
 	//renormalize
-	double sw = sum ( _w );
+	double sw = sum ( w );
 	if ( !std::isfinite ( sw ) ) {
 		for ( int i = 0; i < n; i++ ) {
-			if ( !std::isfinite ( _w ( i ) ) ) {
-				_w ( i ) = 0;
+			if ( !std::isfinite ( w ( i ) ) ) {
+				w ( i ) = 0;
 			}
 		}
-		sw = sum ( _w );
+		sw = sum ( w );
 		if ( !std::isfinite ( sw ) || sw == 0.0 ) {
 			bdm_error ( "Particle filter is lost; no particle is good enough." );
 		}
 	}
-	_w /= sw;
+	w /= sw;
 }
 
 void PF::bayes ( const vec &yt, const vec &cond ) {
-	const vec &ut = cond; //todo
-
-	int i;
-	// generate samples - time step
-	bayes_gensmp ( ut );
+
+/*	if (auxiliary){
+		...
+	}*/
 	// weight them - data step
-	for ( i = 0; i < n; i++ ) {
-		lls ( i ) += obs->evallogcond ( yt, _samples ( i ) ); //+= because lls may have something from gensmp!
+	for (int i = 0; i < n; i++ ) {
+		particles(i)->bayes(yt,cond);
+		lls ( i ) = particles(i)->_ll(); //+= because lls may have something from gensmp!
 	}
 
 	bayes_weights();
 
-	if ( do_resampling() ) {
-		est.resample ( resmethod );
-	}
+	if ( do_resampling() ) { resample ( ); }
 
 }
@@ -75,186 +55,186 @@
 // }
 
-vec MPF::mpfepdf::mean() const {
-	const vec &w = pf->posterior()._w();
-	vec pom = zeros ( BMs ( 0 )->posterior ().dimension() );
-	//compute mean of BMs
-	for ( int i = 0; i < w.length(); i++ ) {
-		pom += BMs ( i )->posterior().mean() * w ( i );
-	}
-	return concat ( pf->posterior().mean(), pom );
-}
-
-vec MPF::mpfepdf::variance() const {
-	const vec &w = pf->posterior()._w();
-
-	vec pom = zeros ( BMs ( 0 )->posterior ().dimension() );
-	vec pom2 = zeros ( BMs ( 0 )->posterior ().dimension() );
-	vec mea;
-
-	for ( int i = 0; i < w.length(); i++ ) {
-		// save current mean
-		mea = BMs ( i )->posterior().mean();
-		pom += mea * w ( i );
-		//compute variance
-		pom2 += ( BMs ( i )->posterior().variance() + pow ( mea, 2 ) ) * w ( i );
-	}
-	return concat ( pf->posterior().variance(), pom2 - pow ( pom, 2 ) );
-}
-
-void MPF::mpfepdf::qbounds ( vec &lb, vec &ub, double perc ) const {
-	//bounds on particles
-	vec lbp;
-	vec ubp;
-	pf->posterior().qbounds ( lbp, ubp );
-
-	//bounds on Components
-	int dimC = BMs ( 0 )->posterior().dimension();
-	int j;
-	// temporary
-	vec lbc ( dimC );
-	vec ubc ( dimC );
-	// minima and maxima
-	vec Lbc ( dimC );
-	vec Ubc ( dimC );
-	Lbc = std::numeric_limits<double>::infinity();
-	Ubc = -std::numeric_limits<double>::infinity();
-
-	for ( int i = 0; i < BMs.length(); i++ ) {
-		// check Coms
-		BMs ( i )->posterior().qbounds ( lbc, ubc );
-		//save either minima or maxima
-		for ( j = 0; j < dimC; j++ ) {
-			if ( lbc ( j ) < Lbc ( j ) ) {
-				Lbc ( j ) = lbc ( j );
-			}
-			if ( ubc ( j ) > Ubc ( j ) ) {
-				Ubc ( j ) = ubc ( j );
-			}
-		}
-	}
-	lb = concat ( lbp, Lbc );
-	ub = concat ( ubp, Ubc );
-}
-
-
-
-void MPF::bayes ( const vec &yt, const vec &cond ) {
-	// follows PF::bayes in most places!!!
-	int i;
-	int n = pf->__w().length();
-	vec &lls = pf->_lls();
-	Array<vec> &samples = pf->__samples();
-
-	// generate samples - time step
-	pf->bayes_gensmp ( vec ( 0 ) );
-	// weight them - data step
-#pragma parallel for
-	for ( i = 0; i < n; i++ ) {
-		vec bm_cond ( BMs ( i )->dimensionc() );
-		this2bm.fill_cond ( yt, cond, bm_cond );
-		pf2bm.filldown ( samples ( i ), bm_cond );
-		BMs ( i ) -> bayes ( this2bm.pushdown ( yt ), bm_cond );
-		lls ( i ) += BMs ( i )->_ll();
-	}
-
-	pf->bayes_weights();
-
-	ivec ind;
-	if ( pf->do_resampling() ) {
-		pf->resample ( ind );
-
-#pragma omp parallel for
-		for ( i = 0; i < n; i++ ) {
-			if ( ind ( i ) != i ) {//replace the current Bm by a new one
-				delete BMs ( i );
-				BMs ( i ) = (BM*) BMs ( ind ( i ) )->_copy(); //copy constructor
-			}
-		};
-	}
-};
-
-
-void MPF_ARXg::bayes ( const vec &yt, const vec &cond ) {
-	// follows PF::bayes in most places!!!
-	int i;
-	int n = pf->__w().length();
-	vec &lls = pf->_lls();
-	Array<vec> &samples = pf->__samples();
-	
-	// generate samples - time step
-	for (int i =0;i<n; i++){
-		mat M; chmat R;
-		BMsp(i)->posterior().sample_mat(M,R);
-		vec diff=R._Ch().T()*randn(samples(i).length());
-		samples(i) = g->eval(samples(i)) + diff;
-		////////////
-//		samples(i) = cond;
-		/////////
-		BMsp(i)->bayes(diff);
-	}
-	// weight them - data step
-	enorm<ldmat> ooo;
-	ooo.set_parameters(zeros(2),0.1*eye(2));
-	ooo.validate();
-	
-	#pragma parallel for
-	for ( i = 0; i < n; i++ ) {
-		vec tmp=yt - h->eval(samples(i));
-		BMso ( i ) -> bayes ( tmp );
-		lls ( i ) += BMso ( i )->_ll();
-		/////
-		ooo._mu()=h->eval(samples(i));
-		lls(i) = ooo.evallog_nn(yt);
-	}
-	
-	pf->bayes_weights();
-	
-	ivec ind;
-	if ( pf->do_resampling() ) {
-		pf->resample ( ind );
-		
-		#pragma omp parallel for
-		for ( i = 0; i < n; i++ ) {
-			if ( ind ( i ) != i ) {//replace the current Bm by a new one
-				delete BMso ( i );
-				BMso ( i ) = (ARX*) BMso ( ind ( i ) )->_copy(); //copy constructor
-				delete BMsp ( i );
-				BMsp ( i ) = (ARX*) BMsp ( ind ( i ) )->_copy(); //copy constructor
-			}
-		};
-	}
-};
-
-void MPF_ARXg::from_setting(const libconfig::Setting& set) {
-	BM::from_setting(set);
-	
-	pf = new PF;
-	// prior must be set before BM
-	pf->prior_from_set ( set );
-	pf->resmethod_from_set ( set );
-		
-	// read functions g and h
-	g=UI::build<fnc>(set,"g");
-	h=UI::build<fnc>(set,"h");
-	
-	set_dim( g->dimension());
-	dimy = h->dimension();
-	
-	shared_ptr<ARX> arxo=UI::build<ARX>(set,"arxo");
-	shared_ptr<ARX> arxp=UI::build<ARX>(set,"arxp");
-	int n = pf->__samples().length();
-	BMso.set_length(n);
-	BMsp.set_length(n);
-	for(int i=0; i<n; i++){
-		BMso(i)=arxo->_copy();
-		BMsp(i)=arxp->_copy();
-	}
-	
-	yrv = arxo->_yrv();
-	//rvc = arxo->_rvc();
-	
-	validate();
-	
-}
+// vec MPF::mpfepdf::mean() const {
+// 	const vec &w = pf->posterior()._w();
+// 	vec pom = zeros ( BMs ( 0 )->posterior ().dimension() );
+// 	//compute mean of BMs
+// 	for ( int i = 0; i < w.length(); i++ ) {
+// 		pom += BMs ( i )->posterior().mean() * w ( i );
+// 	}
+// 	return concat ( pf->posterior().mean(), pom );
+// }
+// 
+// vec MPF::mpfepdf::variance() const {
+// 	const vec &w = pf->posterior()._w();
+// 
+// 	vec pom = zeros ( BMs ( 0 )->posterior ().dimension() );
+// 	vec pom2 = zeros ( BMs ( 0 )->posterior ().dimension() );
+// 	vec mea;
+// 
+// 	for ( int i = 0; i < w.length(); i++ ) {
+// 		// save current mean
+// 		mea = BMs ( i )->posterior().mean();
+// 		pom += mea * w ( i );
+// 		//compute variance
+// 		pom2 += ( BMs ( i )->posterior().variance() + pow ( mea, 2 ) ) * w ( i );
+// 	}
+// 	return concat ( pf->posterior().variance(), pom2 - pow ( pom, 2 ) );
+// }
+// 
+// void MPF::mpfepdf::qbounds ( vec &lb, vec &ub, double perc ) const {
+// 	//bounds on particles
+// 	vec lbp;
+// 	vec ubp;
+// 	pf->posterior().qbounds ( lbp, ubp );
+// 
+// 	//bounds on Components
+// 	int dimC = BMs ( 0 )->posterior().dimension();
+// 	int j;
+// 	// temporary
+// 	vec lbc ( dimC );
+// 	vec ubc ( dimC );
+// 	// minima and maxima
+// 	vec Lbc ( dimC );
+// 	vec Ubc ( dimC );
+// 	Lbc = std::numeric_limits<double>::infinity();
+// 	Ubc = -std::numeric_limits<double>::infinity();
+// 
+// 	for ( int i = 0; i < BMs.length(); i++ ) {
+// 		// check Coms
+// 		BMs ( i )->posterior().qbounds ( lbc, ubc );
+// 		//save either minima or maxima
+// 		for ( j = 0; j < dimC; j++ ) {
+// 			if ( lbc ( j ) < Lbc ( j ) ) {
+// 				Lbc ( j ) = lbc ( j );
+// 			}
+// 			if ( ubc ( j ) > Ubc ( j ) ) {
+// 				Ubc ( j ) = ubc ( j );
+// 			}
+// 		}
+// 	}
+// 	lb = concat ( lbp, Lbc );
+// 	ub = concat ( ubp, Ubc );
+// }
+// 
+// 
+// 
+// void MPF::bayes ( const vec &yt, const vec &cond ) {
+// 	// follows PF::bayes in most places!!!
+// 	int i;
+// 	int n = pf->__w().length();
+// 	vec &lls = pf->_lls();
+// 	Array<vec> &samples = pf->__samples();
+// 
+// 	// generate samples - time step
+// 	pf->bayes_gensmp ( vec ( 0 ) );
+// 	// weight them - data step
+// #pragma parallel for
+// 	for ( i = 0; i < n; i++ ) {
+// 		vec bm_cond ( BMs ( i )->dimensionc() );
+// 		this2bm.fill_cond ( yt, cond, bm_cond );
+// 		pf2bm.filldown ( samples ( i ), bm_cond );
+// 		BMs ( i ) -> bayes ( this2bm.pushdown ( yt ), bm_cond );
+// 		lls ( i ) += BMs ( i )->_ll();
+// 	}
+// 
+// 	pf->bayes_weights();
+// 
+// 	ivec ind;
+// 	if ( pf->do_resampling() ) {
+// 		pf->resample ( ind );
+// 
+// #pragma omp parallel for
+// 		for ( i = 0; i < n; i++ ) {
+// 			if ( ind ( i ) != i ) {//replace the current Bm by a new one
+// 				delete BMs ( i );
+// 				BMs ( i ) = (BM*) BMs ( ind ( i ) )->_copy(); //copy constructor
+// 			}
+// 		};
+// 	}
+// };
+// 
+// 
+// void MPF_ARXg::bayes ( const vec &yt, const vec &cond ) {
+// 	// follows PF::bayes in most places!!!
+// 	int i;
+// 	int n = pf->__w().length();
+// 	vec &lls = pf->_lls();
+// 	Array<vec> &samples = pf->__samples();
+// 	
+// 	// generate samples - time step
+// 	for (int i =0;i<n; i++){
+// 		mat M; chmat R;
+// 		BMsp(i)->posterior().sample_mat(M,R);
+// 		vec diff=R._Ch().T()*randn(samples(i).length());
+// 		samples(i) = g->eval(samples(i)) + diff;
+// 		////////////
+// //		samples(i) = cond;
+// 		/////////
+// 		BMsp(i)->bayes(diff);
+// 	}
+// 	// weight them - data step
+// 	enorm<ldmat> ooo;
+// 	ooo.set_parameters(zeros(2),0.1*eye(2));
+// 	ooo.validate();
+// 	
+// 	#pragma parallel for
+// 	for ( i = 0; i < n; i++ ) {
+// 		vec tmp=yt - h->eval(samples(i));
+// 		BMso ( i ) -> bayes ( tmp );
+// 		lls ( i ) += BMso ( i )->_ll();
+// 		/////
+// 		ooo._mu()=h->eval(samples(i));
+// 		lls(i) = ooo.evallog_nn(yt);
+// 	}
+// 	
+// 	pf->bayes_weights();
+// 	
+// 	ivec ind;
+// 	if ( pf->do_resampling() ) {
+// 		pf->resample ( ind );
+// 		
+// 		#pragma omp parallel for
+// 		for ( i = 0; i < n; i++ ) {
+// 			if ( ind ( i ) != i ) {//replace the current Bm by a new one
+// 				delete BMso ( i );
+// 				BMso ( i ) = (ARX*) BMso ( ind ( i ) )->_copy(); //copy constructor
+// 				delete BMsp ( i );
+// 				BMsp ( i ) = (ARX*) BMsp ( ind ( i ) )->_copy(); //copy constructor
+// 			}
+// 		};
+// 	}
+// };
+// 
+// void MPF_ARXg::from_setting(const libconfig::Setting& set) {
+// 	BM::from_setting(set);
+// 	
+// 	pf = new PF;
+// 	// prior must be set before BM
+// 	pf->prior_from_set ( set );
+// 	pf->resmethod_from_set ( set );
+// 		
+// 	// read functions g and h
+// 	g=UI::build<fnc>(set,"g");
+// 	h=UI::build<fnc>(set,"h");
+// 	
+// 	set_dim( g->dimension());
+// 	dimy = h->dimension();
+// 	
+// 	shared_ptr<ARX> arxo=UI::build<ARX>(set,"arxo");
+// 	shared_ptr<ARX> arxp=UI::build<ARX>(set,"arxp");
+// 	int n = pf->__samples().length();
+// 	BMso.set_length(n);
+// 	BMsp.set_length(n);
+// 	for(int i=0; i<n; i++){
+// 		BMso(i)=arxo->_copy();
+// 		BMsp(i)=arxp->_copy();
+// 	}
+// 	
+// 	yrv = arxo->_yrv();
+// 	//rvc = arxo->_rvc();
+// 	
+// 	validate();
+// 	
+// }
 
 
Index: library/bdm/estim/particles.h
===================================================================
--- library/bdm/estim/particles.h (revision 870)
+++ library/bdm/estim/particles.h (revision 887)
@@ -16,6 +16,167 @@
 
 #include "../estim/arx_ext.h"
+#include "../stat/emix.h"
 
 namespace bdm {
+	
+//! class used in PF
+class MarginalizedParticle : public BM{
+	protected:
+	//! discrte particle
+	dirac est_emp;
+	//! internal Bayes Model
+	shared_ptr<BM> bm; 
+	//! pdf with for transitional par
+	shared_ptr<pdf> par; // pdf for non-linear part
+	//! link from this to bm
+	shared_ptr<datalink_part> cond2bm;
+	//! link from cond to par
+	shared_ptr<datalink_part> cond2par;
+	//! link from emp 2 par
+	shared_ptr<datalink_part> emp2bm;
+	//! link from emp 2 par
+	shared_ptr<datalink_part> emp2par;
+	
+	//! custom posterior - product
+	class eprod_2:public eprod_base {
+		protected:
+		MarginalizedParticle &mp;
+		public:
+		eprod_2(MarginalizedParticle &m):mp(m){}
+		const epdf* factor(int i) const {return (i==0) ? &mp.bm->posterior() : &mp.est_emp;}
+		const int no_factors() const {return 2;}
+	} est;
+	
+	public:
+		MarginalizedParticle():est(*this){};
+		BM* _copy() const{return new MarginalizedParticle(*this);};
+		void bayes(const vec &dt, const vec &cond){
+			vec par_cond(par->dimensionc());
+			cond2par->filldown(cond,par_cond); // copy ut
+			emp2par->filldown(est_emp._point(),par_cond); // copy xt-1
+			
+			//sample new particle
+			est_emp.set_point(par->samplecond(par_cond));
+			//if (evalll)
+			vec bm_cond(bm->dimensionc());
+			cond2bm->filldown(cond, bm_cond);// set e.g. ut
+			emp2bm->filldown(est_emp._point(), bm_cond);// set e.g. ut
+			bm->bayes(dt, bm_cond);
+			ll=bm->_ll();
+		}
+		const eprod_2& posterior() const {return est;}
+	
+		void set_prior(const epdf *pdf0){
+			const eprod *ep=dynamic_cast<const eprod*>(pdf0);
+			if (ep){ // full prior
+				bdm_assert(ep->no_factors()==2,"Incompatible prod");
+				bm->set_prior(ep->factor(0));
+				est_emp.set_point(ep->factor(1)->sample());
+			} else {
+				// assume prior is only for emp;
+				est_emp.set_point(pdf0->sample());
+			}
+		}
+
+		/*! parse structure
+		\code
+		class = "BootstrapParticle";
+		parameter_pdf = {class = 'epdf_offspring', ...};
+		observation_pdf = {class = 'epdf_offspring',...};
+		\endcode
+		If rvs are set, then it checks for compatibility.
+		*/
+		void from_setting(const Setting &set){
+			BM::from_setting ( set );
+			par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory );
+			bm = UI::build<BM> ( set, "bm", UI::compulsory );
+		}
+		void validate(){
+			yrv = bm->_rv();
+			dimy = bm->dimension();
+			set_rv( par->_rv());
+			set_dim( par->dimension());
+
+			rvc = par->_rvc().subt(par->_rv().copy_t(-1));
+			rvc.add(bm->_rvc()); //
+			
+			cond2bm=new datalink_part;
+			cond2par=new datalink_part;
+			emp2bm  =new datalink_part;
+			emp2par =new datalink_part;
+			cond2bm->set_connection(bm->_rvc(), rvc);
+			cond2par->set_connection(par->_rvc(), rvc);
+			emp2bm->set_connection(bm->_rvc(), _rv());
+			emp2par->set_connection(par->_rvc(), _rv().copy_t(-1));
+			
+			dimc = rvc._dsize();
+		};
+};
+UIREGISTER(MarginalizedParticle);
+
+//! class used in PF
+class BootstrapParticle : public BM{
+	dirac est;
+	shared_ptr<pdf> par;
+	shared_ptr<pdf> obs;
+	shared_ptr<datalink_part> cond2par;
+	shared_ptr<datalink_part> cond2obs;
+	shared_ptr<datalink_part> xt2obs;
+	shared_ptr<datalink_part> xtm2par;
+	public:
+		BM* _copy() const{return new BootstrapParticle(*this);};
+		void bayes(const vec &dt, const vec &cond){
+			vec par_cond(par->dimensionc());
+			cond2par->filldown(cond,par_cond); // copy ut
+			xtm2par->filldown(est._point(),par_cond); // copy xt-1
+			
+			//sample new particle
+			est.set_point(par->samplecond(par_cond));
+			//if (evalll)
+			vec obs_cond(obs->dimensionc());
+			cond2obs->filldown(cond, obs_cond);// set e.g. ut
+			xt2obs->filldown(est._point(), obs_cond);// set e.g. ut
+			ll=obs->evallogcond(dt,obs_cond);
+		}
+		const dirac& posterior() const {return est;}
+		
+		void set_prior(const epdf *pdf0){est.set_point(pdf0->sample());}
+		
+		/*! parse structure
+		\code
+		class = "BootstrapParticle";
+		parameter_pdf = {class = 'epdf_offspring', ...};
+		observation_pdf = {class = 'epdf_offspring',...};
+		\endcode
+		If rvs are set, then it checks for compatibility.
+		*/
+		void from_setting(const Setting &set){
+			BM::from_setting ( set );
+			par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory );
+			obs = UI::build<pdf> ( set, "observation_pdf", UI::compulsory );
+		}
+		void validate(){
+			yrv = obs->_rv();
+			dimy = obs->dimension();
+			set_rv( par->_rv());
+			set_dim( par->dimension());
+			
+			rvc = par->_rvc().subt(par->_rv().copy_t(-1));
+			rvc.add(obs->_rvc()); //
+			
+			cond2obs=new datalink_part;
+			cond2par=new datalink_part;
+			xt2obs  =new datalink_part;
+			xtm2par =new datalink_part;
+			cond2obs->set_connection(obs->_rvc(), rvc);
+			cond2par->set_connection(par->_rvc(), rvc);
+			xt2obs->set_connection(obs->_rvc(), _rv());
+			xtm2par->set_connection(par->_rvc(), _rv().copy_t(-1));
+			
+			dimc = rvc._dsize();
+		};
+};
+UIREGISTER(BootstrapParticle);
+
 
 /*!
@@ -33,17 +194,20 @@
 	LOG_LEVEL(PF,weights,samples);
 	
+	class pf_mix: public emix_base{
+		Array<BM*> &bms;
+		public:
+			pf_mix(vec &w0, Array<BM*> &bms0):emix_base(w0),bms(bms0){}
+			const epdf* component(const int &i)const{return &(bms(i)->posterior());}
+			int no_coms() const {return bms.length();}
+	};
 protected:
 	//!number of particles;
 	int n;
 	//!posterior density
-	eEmp est;
-	//! pointer into \c eEmp
-	vec &_w;
-	//! pointer into \c eEmp
-	Array<vec> &_samples;
-	//! Parameter evolution model
-	shared_ptr<pdf> par;
-	//! Observation model
-	shared_ptr<pdf> obs;
+	pf_mix est;
+	//! weights;
+	vec w;
+	//! particles
+	Array<BM*> particles;
 	//! internal structure storing loglikelihood of predictions
 	vec lls;
@@ -62,6 +226,5 @@
 	//! \name Constructors
 	//!@{
-	PF ( ) : est(), _w ( est._w() ), _samples ( est._samples() ) {
-	};
+	PF ( ) : est(w,particles) { };
 
 	void set_parameters ( int n0, double res_th0 = 0.5, RESAMPLING_METHOD rm = SYSTEMATIC ) {
@@ -70,34 +233,32 @@
 		resmethod = rm;
 	};
-	void set_model ( shared_ptr<pdf> par0, shared_ptr<pdf> obs0 ) {
-		par = par0;
-		obs = obs0;
+	void set_model ( const BM *particle0, const epdf *prior) {
+		if (n>0){
+			particles.set_length(n);
+			for (int i=0; i<n;i++){
+				particles(i) = particle0->_copy();
+				particles(i)->set_prior(prior);
+			}
+		}
 		// set values for posterior
-		est.set_rv ( par->_rv() );
+		est.set_rv ( particle0->posterior()._rv() );
 	};
 	void set_statistics ( const vec w0, const epdf &epdf0 ) {
-		est.set_statistics ( w0, epdf0 );
+		//est.set_statistics ( w0, epdf0 );
 	};
-	void set_statistics ( const eEmp &epdf0 ) {
+/*	void set_statistics ( const eEmp &epdf0 ) {
 		bdm_assert_debug ( epdf0._rv().equal ( par->_rv() ), "Incompatible input" );
 		est = epdf0;
-	};
+	};*/
 	//!@}
 
-	//! Set posterior density by sampling from epdf0
-	//! bayes I - generate samples and add their weights to lls
-	virtual void bayes_gensmp ( const vec &ut );
-	//! bayes II - compute weights of the
+	//! bayes compute weights of the
 	virtual void bayes_weights();
 	//! important part of particle filtering - decide if it is time to perform resampling
 	virtual bool do_resampling() {
-		double eff = 1.0 / ( _w * _w );
+		double eff = 1.0 / ( w * w );
 		return eff < ( res_threshold*n );
 	}
 	void bayes ( const vec &yt, const vec &cond );
-	//!access function
-	vec& __w() {
-		return _w;
-	}
 	//!access function
 	vec& _lls() {
@@ -109,5 +270,5 @@
 	}
 	//! return correctly typed posterior (covariant return)
-	const eEmp& posterior() const {
+	const pf_mix& posterior() const {
 		return est;
 	}
@@ -127,20 +288,48 @@
 	void from_setting ( const Setting &set ) {
 		BM::from_setting ( set );
-		par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory );
-		obs = UI::build<pdf> ( set, "observation_pdf", UI::compulsory );
-
-		prior_from_set ( set );
+		
+		shared_ptr<BM> bm0 = UI::build<BM>(set, "particle",UI::compulsory);
+		
+		shared_ptr<epdf> pri = UI::build<epdf> ( set, "prior", UI::compulsory );
+		n =0;
+		UI::get(n,set,"n",UI::optional);;
+		if (n>0){
+			particles.set_length(n);
+			for(int i=0;i<n;i++){particles(i)=bm0->_copy();}
+			w = ones(n)/n;
+		}
+		set_prior(pri.get());
+		// set resampling method
 		resmethod_from_set ( set );
-		// set resampling method
 		//set drv
-		//find potential input - what remains in rvc when we subtract rv
-		RV u = par->_rvc().remove_time().subt ( par->_rv() );
-		//find potential input - what remains in rvc when we subtract x_t
-		RV obs_u = obs->_rvc().remove_time().subt ( par->_rv() );
-
-		u.add ( obs_u ); // join both u, and check if they do not overlap
-
-		set_yrv ( obs->_rv() );
-		rvc = u;
+
+		set_yrv ( bm0->_rv() );
+		rvc = bm0->_rvc();
+		BM::set_rv(bm0->_rv());
+		yrv=bm0->_yrv();
+	}
+	
+	void set_prior(const epdf *pri){
+		const emix_base *emi=dynamic_cast<const emix_base*>(pri);
+		if (emi) {
+			bdm_assert(particles.length()>0, "initial particle is not assigned");
+			n = emi->_w().length();
+			int old_n = particles.length();
+			if (n!=old_n){
+				particles.set_length(n,true);
+			} 
+			for(int i=old_n;i<n;i++){particles(i)=particles(0)->_copy();}
+			
+			for (int i =0; i<n; i++){
+				particles(i)->set_prior(emi->_com(i));
+			}
+		} else {
+			// try to find "n"
+			bdm_assert(n>0, "Field 'n' must be filled when prior is not of type emix");
+			for (int i =0; i<n; i++){
+				particles(i)->set_prior(pri);
+			}
+			
+		}
 	}
 	//! auxiliary function reading parameter 'resmethod' from configuration file
@@ -165,51 +354,37 @@
 		};
 		if ( !UI::get ( res_threshold, set, "res_threshold", UI::optional ) ) {
-			res_threshold = 0.5;
+			res_threshold = 0.9;
 		}
 		//validate();
-	}
-	//! load prior information from set and set internal structures accordingly
-	void prior_from_set ( const Setting & set ) {
-		shared_ptr<epdf> pri = UI::build<epdf> ( set, "prior", UI::compulsory );
-
-		eEmp *test_emp = dynamic_cast<eEmp*> ( & ( *pri ) );
-		if ( test_emp ) { // given pdf is sampled
-			est = *test_emp;
-		} else {
-			//int n;
-			if ( !UI::get ( n, set, "n", UI::optional ) ) {
-				n = 10;
-			}
-			// sample from prior
-			set_statistics ( ones ( n ) / n, *pri );
-		}
-		n = est._w().length();
-		lls=zeros(n);
 	}
 
 	void validate() {
 		BM::validate();
-		bdm_assert ( par, "PF::parameter_pdf not specified." );
-		n = _w.length();
+		est.validate();
+		bdm_assert ( n>0, "empty particle pool" );
+		n = w.length();
 		lls = zeros ( n );
 
-		if ( par->_rv()._dsize() > 0 ) {
-			bdm_assert ( par->_rv()._dsize() == est.dimension(), "Mismatch of RV and dimension of posterior" );
+		if ( particles(0)->_rv()._dsize() > 0 ) {
+			bdm_assert (  particles(0)->_rv()._dsize() == est.dimension(), "Mismatch of RV and dimension of posterior" );
 		}
 	}
 	//! resample posterior density (from outside - see MPF)
-	void resample ( ivec &ind ) {
-		est.resample ( ind, resmethod );
+	void resample ( ) {
+		ivec ind = zeros_i ( n );
+		bdm::resample(w,ind,resmethod);
+		// copy the internals according to ind
+		for (int i = 0; i < n; i++ ) {
+			if ( ind ( i ) != i ) {
+				particles( i ) = particles( ind ( i ) )->_copy();
+			}
+			w ( i ) = 1.0 / n;
+		}
 	}
 	//! access function
-	Array<vec>& __samples() {
-		return _samples;
-	}
-
-	virtual double logpred ( const vec &yt ) const NOT_IMPLEMENTED(0);
-
-	virtual epdf* epredictor() const NOT_IMPLEMENTED(NULL);
-	
-	virtual pdf* predictor() const NOT_IMPLEMENTED(NULL);
+	Array<BM*>& _particles() {
+		return particles;
+	}
+
 };
 UIREGISTER ( PF );
@@ -222,201 +397,200 @@
 */
 
-class MPF : public BM  {
-	//! \var log_level_enums means 
-	//! TODO DOPLNIT
-	LOG_LEVEL(MPF,means);
-
-protected:
-	//! particle filter on non-linear variable
-	shared_ptr<PF> pf;
-	//! Array of Bayesian models
-	Array<BM*> BMs;
-
-	//! internal class for pdf providing composition of eEmp with external components
-
-	class mpfepdf : public epdf  {
-		//! pointer to particle filter
-		shared_ptr<PF> &pf;
-		//! pointer to Array of BMs
-		Array<BM*> &BMs;
-	public:
-		//! constructor
-		mpfepdf ( shared_ptr<PF> &pf0, Array<BM*> &BMs0 ) : epdf(), pf ( pf0 ), BMs ( BMs0 ) { };
-		//! a variant of set parameters - this time, parameters are read from BMs and pf
-		void read_parameters() {
-			rv = concat ( pf->posterior()._rv(), BMs ( 0 )->posterior()._rv() );
-			dim = pf->posterior().dimension() + BMs ( 0 )->posterior().dimension();
-			bdm_assert_debug ( dim == rv._dsize(), "Wrong name " );
-		}
-		vec mean() const;
-
-		vec variance() const;
-
-		void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const;
-
-		vec sample() const NOT_IMPLEMENTED(0);
-
-		double evallog ( const vec &val ) const NOT_IMPLEMENTED(0);		
-	};
-
-	//! Density joining PF.est with conditional parts
-	mpfepdf jest;
-
-	//! datalink from global yt and cond (Up) to BMs yt and cond (Down)
-	datalink_m2m this2bm;
-	//! datalink from global yt and cond (Up) to PFs yt and cond (Down)
-	datalink_m2m this2pf;
-	//!datalink from PF part to BM
-	datalink_part pf2bm;
-
-public:
-	//! Default constructor.
-	MPF () :  jest ( pf, BMs ) {};
-	//! set all parameters at once
-	void set_pf ( shared_ptr<pdf> par0, int n0, RESAMPLING_METHOD rm = SYSTEMATIC ) {
-		if (!pf) pf=new PF;
-		pf->set_model ( par0, par0 ); // <=== nasty!!!
-		pf->set_parameters ( n0, rm );
-		pf->set_rv(par0->_rv());
-		BMs.set_length ( n0 );
-	}
-	//! set a prototype of BM, copy it to as many times as there is particles in pf
-	void set_BM ( const BM &BMcond0 ) {
-
-		int n = pf->__w().length();
-		BMs.set_length ( n );
-		// copy
-		//BMcond0 .condition ( pf->posterior()._sample ( 0 ) );
-		for ( int i = 0; i < n; i++ ) {
-			BMs ( i ) = (BM*) BMcond0._copy();
-		}
-	};
-
-	void bayes ( const vec &yt, const vec &cond );
-
-	const epdf& posterior() const {
-		return jest;
-	}
-
-	//!Access function
-	const BM* _BM ( int i ) {
-		return BMs ( i );
-	}
-	PF& _pf() {return *pf;}
-
-
-	virtual double logpred ( const vec &yt ) const NOT_IMPLEMENTED(0); 
-		
-	virtual epdf* epredictor() const NOT_IMPLEMENTED(NULL); 
-	
-	virtual pdf* predictor() const NOT_IMPLEMENTED(NULL); 
-
-
-	/*! configuration structure for basic PF
-	\code
-	BM              = BM_class;           // Bayesian filtr for analytical part of the model
-	parameter_pdf   = pdf_class;         // transitional pdf for non-parametric part of the model
-	prior           = epdf_class;         // prior probability density
-	--- optional ---
-	n               = 10;                 // number of particles
-	resmethod       = 'systematic', or 'multinomial', or 'stratified'
-										  // resampling method
-	\endcode
-	*/
-	void from_setting ( const Setting &set ) {
-		BM::from_setting( set );
-
-		shared_ptr<pdf> par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory );
-
-		pf = new PF;
-		// prior must be set before BM
-		pf->prior_from_set ( set );
-		pf->resmethod_from_set ( set );
-		pf->set_model ( par, par ); // too hackish!
-
-		shared_ptr<BM> BM0 = UI::build<BM> ( set, "BM", UI::compulsory );
-		set_BM ( *BM0 );
-
-		//set drv
-		//??set_yrv(concat(BM0->_yrv(),u) );
-		set_yrv ( BM0->_yrv() );
-		rvc = BM0->_rvc().subt ( par->_rv() );
-		//find potential input - what remains in rvc when we subtract rv
-		RV u = par->_rvc().subt ( par->_rv().copy_t ( -1 ) );
-		rvc.add ( u );
-		dimc = rvc._dsize();
-		validate();
-	}
-
-	void validate() {
-		BM::validate();
-		try {
-			pf->validate();
-		} catch ( std::exception ) {
-			throw UIException ( "Error in PF part of MPF:" );
-		}
-		jest.read_parameters();
-		this2bm.set_connection ( BMs ( 0 )->_yrv(), BMs ( 0 )->_rvc(), yrv, rvc );
-		this2pf.set_connection ( pf->_yrv(), pf->_rvc(), yrv, rvc );
-		pf2bm.set_connection ( BMs ( 0 )->_rvc(), pf->posterior()._rv() );
-	}
-};
-UIREGISTER ( MPF );
+// class MPF : public BM  {
+// 	//! Introduces new option 
+// 	//! \li means - meaning TODO
+// 	LOG_LEVEL(MPF,means);
+// protected:
+// 	//! particle filter on non-linear variable
+// 	shared_ptr<PF> pf;
+// 	//! Array of Bayesian models
+// 	Array<BM*> BMs;
+// 
+// 	//! internal class for pdf providing composition of eEmp with external components
+// 
+// 	class mpfepdf : public epdf  {
+// 		//! pointer to particle filter
+// 		shared_ptr<PF> &pf;
+// 		//! pointer to Array of BMs
+// 		Array<BM*> &BMs;
+// 	public:
+// 		//! constructor
+// 		mpfepdf ( shared_ptr<PF> &pf0, Array<BM*> &BMs0 ) : epdf(), pf ( pf0 ), BMs ( BMs0 ) { };
+// 		//! a variant of set parameters - this time, parameters are read from BMs and pf
+// 		void read_parameters() {
+// 			rv = concat ( pf->posterior()._rv(), BMs ( 0 )->posterior()._rv() );
+// 			dim = pf->posterior().dimension() + BMs ( 0 )->posterior().dimension();
+// 			bdm_assert_debug ( dim == rv._dsize(), "Wrong name " );
+// 		}
+// 		vec mean() const;
+// 
+// 		vec variance() const;
+// 
+// 		void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const;
+// 
+// 		vec sample() const NOT_IMPLEMENTED(0);
+// 
+// 		double evallog ( const vec &val ) const NOT_IMPLEMENTED(0);		
+// 	};
+// 
+// 	//! Density joining PF.est with conditional parts
+// 	mpfepdf jest;
+// 
+// 	//! datalink from global yt and cond (Up) to BMs yt and cond (Down)
+// 	datalink_m2m this2bm;
+// 	//! datalink from global yt and cond (Up) to PFs yt and cond (Down)
+// 	datalink_m2m this2pf;
+// 	//!datalink from PF part to BM
+// 	datalink_part pf2bm;
+// 
+// public:
+// 	//! Default constructor.
+// 	MPF () :  jest ( pf, BMs ) {};
+// 	//! set all parameters at once
+// 	void set_pf ( shared_ptr<pdf> par0, int n0, RESAMPLING_METHOD rm = SYSTEMATIC ) {
+// 		if (!pf) pf=new PF;
+// 		pf->set_model ( par0, par0 ); // <=== nasty!!!
+// 		pf->set_parameters ( n0, rm );
+// 		pf->set_rv(par0->_rv());
+// 		BMs.set_length ( n0 );
+// 	}
+// 	//! set a prototype of BM, copy it to as many times as there is particles in pf
+// 	void set_BM ( const BM &BMcond0 ) {
+// 
+// 		int n = pf->__w().length();
+// 		BMs.set_length ( n );
+// 		// copy
+// 		//BMcond0 .condition ( pf->posterior()._sample ( 0 ) );
+// 		for ( int i = 0; i < n; i++ ) {
+// 			BMs ( i ) = (BM*) BMcond0._copy();
+// 		}
+// 	};
+// 
+// 	void bayes ( const vec &yt, const vec &cond );
+// 
+// 	const epdf& posterior() const {
+// 		return jest;
+// 	}
+// 
+// 	//!Access function
+// 	const BM* _BM ( int i ) {
+// 		return BMs ( i );
+// 	}
+// 	PF& _pf() {return *pf;}
+// 
+// 
+// 	virtual double logpred ( const vec &yt ) const NOT_IMPLEMENTED(0); 
+// 		
+// 	virtual epdf* epredictor() const NOT_IMPLEMENTED(NULL); 
+// 	
+// 	virtual pdf* predictor() const NOT_IMPLEMENTED(NULL); 
+// 
+// 
+// 	/*! configuration structure for basic PF
+// 	\code
+// 	BM              = BM_class;           // Bayesian filtr for analytical part of the model
+// 	parameter_pdf   = pdf_class;         // transitional pdf for non-parametric part of the model
+// 	prior           = epdf_class;         // prior probability density
+// 	--- optional ---
+// 	n               = 10;                 // number of particles
+// 	resmethod       = 'systematic', or 'multinomial', or 'stratified'
+// 										  // resampling method
+// 	\endcode
+// 	*/
+// 	void from_setting ( const Setting &set ) {
+// 		BM::from_setting( set );
+// 
+// 		shared_ptr<pdf> par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory );
+// 
+// 		pf = new PF;
+// 		// prior must be set before BM
+// 		pf->prior_from_set ( set );
+// 		pf->resmethod_from_set ( set );
+// 		pf->set_model ( par, par ); // too hackish!
+// 
+// 		shared_ptr<BM> BM0 = UI::build<BM> ( set, "BM", UI::compulsory );
+// 		set_BM ( *BM0 );
+// 
+// 		//set drv
+// 		//??set_yrv(concat(BM0->_yrv(),u) );
+// 		set_yrv ( BM0->_yrv() );
+// 		rvc = BM0->_rvc().subt ( par->_rv() );
+// 		//find potential input - what remains in rvc when we subtract rv
+// 		RV u = par->_rvc().subt ( par->_rv().copy_t ( -1 ) );
+// 		rvc.add ( u );
+// 		dimc = rvc._dsize();
+// 		validate();
+// 	}
+// 
+// 	void validate() {
+// 		BM::validate();
+// 		try {
+// 			pf->validate();
+// 		} catch ( std::exception ) {
+// 			throw UIException ( "Error in PF part of MPF:" );
+// 		}
+// 		jest.read_parameters();
+// 		this2bm.set_connection ( BMs ( 0 )->_yrv(), BMs ( 0 )->_rvc(), yrv, rvc );
+// 		this2pf.set_connection ( pf->_yrv(), pf->_rvc(), yrv, rvc );
+// 		pf2bm.set_connection ( BMs ( 0 )->_rvc(), pf->posterior()._rv() );
+// 	}
+// };
+// UIREGISTER ( MPF );
 
 /*! ARXg for estimation of state-space variances
 */
-class MPF_ARXg :public BM{
-	protected:
-	shared_ptr<PF> pf;
-	//! pointer to Array of BMs
-	Array<ARX*> BMso;
-	//! pointer to Array of BMs
-	Array<ARX*> BMsp;
-	//!parameter evolution
-	shared_ptr<fnc> g;
-	//!observation function
-	shared_ptr<fnc> h;
-	
-	public:
-		void bayes(const vec &yt, const vec &cond );
-		void from_setting(const Setting &set) ;
-		void validate() {
-			bdm_assert(g->dimensionc()==g->dimension(),"not supported yet");
-			bdm_assert(h->dimensionc()==g->dimension(),"not supported yet");			
-		}
-
-		double logpred(const vec &cond) const NOT_IMPLEMENTED(0.0);
-		epdf* epredictor() const NOT_IMPLEMENTED(NULL);
-		pdf* predictor() const NOT_IMPLEMENTED(NULL);
-		const epdf& posterior() const {return pf->posterior();};
-		
-		void log_register( logger &L, const string &prefix ){
-			BM::log_register(L,prefix);
-			logrec->ids.set_size ( 3 );
-			logrec->ids(1)= L.add_vector(RV("Q",dimension()*dimension()), prefix+L.prefix_sep()+"Q");
-			logrec->ids(2)= L.add_vector(RV("R",dimensiony()*dimensiony()), prefix+L.prefix_sep()+"R");
-			
-		};
-		void log_write() const {
-			BM::log_write();
-			mat mQ=zeros(dimension(),dimension());
-			mat pom=zeros(dimension(),dimension());
-			mat mR=zeros(dimensiony(),dimensiony());
-			mat pom2=zeros(dimensiony(),dimensiony());
-			mat dum;
-			const vec w=pf->posterior()._w();
-			for (int i=0; i<w.length(); i++){
-				BMsp(i)->posterior().mean_mat(dum,pom);
-				mQ += w(i) * pom;
-				BMso(i)->posterior().mean_mat(dum,pom2);
-				mR += w(i) * pom2;
-				
-			}
-			logrec->L.log_vector ( logrec->ids ( 1 ), cvectorize(mQ) );
-			logrec->L.log_vector ( logrec->ids ( 2 ), cvectorize(mR) );
-			
-		}
-};
-UIREGISTER(MPF_ARXg);
+// class MPF_ARXg :public BM{
+// 	protected:
+// 	shared_ptr<PF> pf;
+// 	//! pointer to Array of BMs
+// 	Array<ARX*> BMso;
+// 	//! pointer to Array of BMs
+// 	Array<ARX*> BMsp;
+// 	//!parameter evolution
+// 	shared_ptr<fnc> g;
+// 	//!observation function
+// 	shared_ptr<fnc> h;
+// 	
+// 	public:
+// 		void bayes(const vec &yt, const vec &cond );
+// 		void from_setting(const Setting &set) ;
+// 		void validate() {
+// 			bdm_assert(g->dimensionc()==g->dimension(),"not supported yet");
+// 			bdm_assert(h->dimensionc()==g->dimension(),"not supported yet");			
+// 		}
+// 
+// 		double logpred(const vec &cond) const NOT_IMPLEMENTED(0.0);
+// 		epdf* epredictor() const NOT_IMPLEMENTED(NULL);
+// 		pdf* predictor() const NOT_IMPLEMENTED(NULL);
+// 		const epdf& posterior() const {return pf->posterior();};
+// 		
+// 		void log_register( logger &L, const string &prefix ){
+// 			BM::log_register(L,prefix);
+// 			logrec->ids.set_size ( 3 );
+// 			logrec->ids(1)= L.add_vector(RV("Q",dimension()*dimension()), prefix+L.prefix_sep()+"Q");
+// 			logrec->ids(2)= L.add_vector(RV("R",dimensiony()*dimensiony()), prefix+L.prefix_sep()+"R");
+// 			
+// 		};
+// 		void log_write() const {
+// 			BM::log_write();
+// 			mat mQ=zeros(dimension(),dimension());
+// 			mat pom=zeros(dimension(),dimension());
+// 			mat mR=zeros(dimensiony(),dimensiony());
+// 			mat pom2=zeros(dimensiony(),dimensiony());
+// 			mat dum;
+// 			const vec w=pf->posterior()._w();
+// 			for (int i=0; i<w.length(); i++){
+// 				BMsp(i)->posterior().mean_mat(dum,pom);
+// 				mQ += w(i) * pom;
+// 				BMso(i)->posterior().mean_mat(dum,pom2);
+// 				mR += w(i) * pom2;
+// 				
+// 			}
+// 			logrec->L.log_vector ( logrec->ids ( 1 ), cvectorize(mQ) );
+// 			logrec->L.log_vector ( logrec->ids ( 2 ), cvectorize(mR) );
+// 			
+// 		}
+// };
+// UIREGISTER(MPF_ARXg);
 
 
Index: library/bdm/stat/exp_family.h
===================================================================
--- library/bdm/stat/exp_family.h (revision 878)
+++ library/bdm/stat/exp_family.h (revision 887)
@@ -97,9 +97,4 @@
 	virtual void flatten ( const BMEF * B ) NOT_IMPLEMENTED_VOID;
 
-	double logpred ( const vec &yt ) const NOT_IMPLEMENTED(0);
-	
-	virtual epdf* epredictor() const NOT_IMPLEMENTED(NULL);
-
-	virtual pdf* predictor() const NOT_IMPLEMENTED(NULL);
 
 	void to_setting ( Setting &set ) const
@@ -267,4 +262,8 @@
 UIREGISTER2 ( enorm, fsqmat );
 SHAREDPTR2 ( enorm, fsqmat );
+
+typedef enorm<ldmat> egauss;
+UIREGISTER(egauss);
+
 
 //forward declaration
@@ -1596,4 +1595,8 @@
 //! Switch between various resampling methods.
 enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
+
+//! Shortcut for multinomial.sample(int n). Various simplifications are allowed see RESAMPLING_METHOD 
+void resample(const vec &w, ivec &ind, RESAMPLING_METHOD=SYSTEMATIC);
+
 /*!
 \brief Weighted empirical density
@@ -1664,12 +1667,5 @@
 	};
 	//! Function performs resampling, i.e. removal of low-weight samples and duplication of high-weight samples such that the new samples represent the same density.
-	//! The vector with indeces of new samples is returned in variable \c index.
-	void resample ( ivec &index, RESAMPLING_METHOD method = SYSTEMATIC );
-
-	//! Resampling without returning index of new particles.
-	void resample ( RESAMPLING_METHOD method = SYSTEMATIC ) {
-		ivec ind;
-		resample ( ind, method );
-	};
+	void resample ( RESAMPLING_METHOD method = SYSTEMATIC );
 
 	//! inherited operation : NOT implemented
@@ -1874,4 +1870,19 @@
 }
 
+/*! Dirac delta function distribution */
+class dirac: public epdf{
+	protected: 
+		vec point;
+	public:
+		double evallog (const vec &dt) const {return -inf;}
+		vec mean () const {return point;}
+		vec variance () const {return pow(point,2);}
+		void qbounds ( vec &lb, vec &ub, double percentage = 0.95 ) const { lb = point; ub = point;}
+		//! access
+		const vec& _point() {return point;}
+		void set_point(const vec& p){point =p; dim=p.length();}
+		vec sample() const {return point;}
+	};
+
 ////
 ///////
Index: library/bdm/stat/exp_family.cpp
===================================================================
--- library/bdm/stat/exp_family.cpp (revision 878)
+++ library/bdm/stat/exp_family.cpp (revision 887)
@@ -423,5 +423,18 @@
 }
 
-void eEmp::resample ( ivec &ind, RESAMPLING_METHOD method ) {
+void eEmp::resample ( RESAMPLING_METHOD method ) {
+	ivec ind = zeros_i ( n );
+	bdm::resample(w,ind,method);
+	// copy the internals according to ind
+	for (int i = 0; i < n; i++ ) {
+		if ( ind ( i ) != i ) {
+			samples ( i ) = samples ( ind ( i ) );
+		}
+		w ( i ) = 1.0 / n;
+	}
+}
+
+void resample ( const vec &w, ivec &ind, RESAMPLING_METHOD method ) {
+	int n = w.length();
 	ind = zeros_i ( n );
 	ivec N_babies = zeros_i ( n );
@@ -430,42 +443,42 @@
 	int i, j, parent;
 	double u0;
-
+	
 	switch ( method ) {
-	case MULTINOMIAL:
-		u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n );
-
-		for ( i = n - 2; i >= 0; i-- ) {
-			u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) );
-		}
-
-		break;
-
-	case STRATIFIED:
-
-		for ( i = 0; i < n; i++ ) {
-			u ( i ) = ( i + UniRNG.sample() ) / n;
-		}
-
-		break;
-
-	case SYSTEMATIC:
-		u0 = UniRNG.sample();
-
-		for ( i = 0; i < n; i++ ) {
-			u ( i ) = ( i + u0 ) / n;
-		}
-
-		break;
-
-	default:
-		bdm_error ( "PF::resample(): Unknown resampling method" );
-	}
-
+		case MULTINOMIAL:
+			u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n );
+			
+			for ( i = n - 2; i >= 0; i-- ) {
+				u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) );
+			}
+			
+			break;
+			
+		case STRATIFIED:
+			
+			for ( i = 0; i < n; i++ ) {
+				u ( i ) = ( i + UniRNG.sample() ) / n;
+			}
+			
+			break;
+			
+		case SYSTEMATIC:
+			u0 = UniRNG.sample();
+			
+			for ( i = 0; i < n; i++ ) {
+				u ( i ) = ( i + u0 ) / n;
+			}
+			
+			break;
+			
+		default:
+			bdm_error ( "PF::resample(): Unknown resampling method" );
+	}
+	
 	// U is now full
 	j = 0;
-
+	
 	for ( i = 0; i < n; i++ ) {
 		while ( u ( i ) > cumDist ( j ) ) j++;
-
+		
 		N_babies ( j ) ++;
 	}
@@ -474,9 +487,9 @@
 	// * particles with at least one baby should not move *
 	// This assures that reassignment can be done inplace;
-
+	
 	// find the first parent;
 	parent = 0;
 	while ( N_babies ( parent ) == 0 ) parent++;
-
+	
 	// Build index
 	for ( i = 0; i < n; i++ ) {
@@ -488,19 +501,11 @@
 			// if yes, find the next one
 			while ( ( N_babies ( parent ) == 0 ) || ( N_babies ( parent ) == 1 && parent > i ) ) parent++;
-
+			
 			// Replicate parent
 			ind ( i ) = parent;
-
+			
 			N_babies ( parent ) --; //this index was now replicated;
 		}
-
-	}
-
-// copy the internals according to ind
-	for ( i = 0; i < n; i++ ) {
-		if ( ind ( i ) != i ) {
-			samples ( i ) = samples ( ind ( i ) );
-		}
-		w ( i ) = 1.0 / n;
+		
 	}
 }
Index: library/tests/stresssuite/particle_stress.cpp
===================================================================
--- library/tests/stresssuite/particle_stress.cpp (revision 722)
+++ library/tests/stresssuite/particle_stress.cpp (revision 887)
@@ -37,5 +37,5 @@
 
 	ivec ind;
-	emp.resample ( ind );
+	resample ( v, ind );
 
 	cout << ind << endl;
Index: library/tests/stresssuite/resample_stress.cpp
===================================================================
--- library/tests/stresssuite/resample_stress.cpp (revision 722)
+++ library/tests/stresssuite/resample_stress.cpp (revision 887)
@@ -38,5 +38,5 @@
 
 	ivec ind;
-	emp.resample ( ind );
+	resample ( v,ind );
 
 	cout << ind << endl;
