Index: bdm/estim/arx.h
===================================================================
--- bdm/estim/arx.h (revision 271)
+++ bdm/estim/arx.h (revision 283)
@@ -55,5 +55,5 @@
 	ARX ( const double frg0=1.0 ) : BMEF ( frg0 ),est (), V ( est._V() ), nu ( est._nu() ) {};
 	ARX ( const ARX &A0 ) : BMEF (),est ( A0.est ), V ( est._V() ), nu ( est._nu() ) {};
-	ARX* _copy_();
+	ARX* _copy_() const;
 	void set_parameters ( double frg0 ) {frg=frg0;}
 	void set_statistics ( int dimx0, const ldmat V0, double nu0=-1.0 ) {est.set_parameters ( dimx0,V0,nu0 );last_lognc=est.lognc();dimx=dimx0;}
Index: bdm/estim/libPF.h
===================================================================
--- bdm/estim/libPF.h (revision 281)
+++ bdm/estim/libPF.h (revision 283)
@@ -40,11 +40,14 @@
 	mpdf *obs;
 
+	//! which resampling method will be used
+	RESAMPLING_METHOD resmethod;
+
 	//! \name Options
 	//!@{
 
-	//! Log resampling times
-	bool opt_L_res;
 	//! Log all samples
 	bool opt_L_smp;
+	//! Log all samples
+	bool opt_L_wei;
 	//!@}
 
@@ -52,17 +55,18 @@
 	//! \name Constructors
 	//!@{
-	PF ( ) :est(), _w ( est._w() ),_samples ( est._samples() ) {};
-	PF ( mpdf *par0, mpdf *obs0, epdf *epdf0, int n0 ) :
-			est ( ),_w ( est._w() ),_samples ( est._samples() )
-	{ set_parameters ( par0,obs0,n0 ); set_statistics ( ones ( n0 ),epdf0 ); };
-	void set_parameters ( mpdf *par0, mpdf *obs0, int n0 )
-	{ par = par0; obs=obs0; n=n0; est.set_n ( n );};
-	void set_statistics ( const vec w0, epdf *epdf0 ) {est.set_parameters ( w0,epdf0 );};
+	PF ( ) :est(), _w ( est._w() ),_samples ( est._samples() ), opt_L_smp ( false ), opt_L_wei ( false ) {LIDs.set_size ( 5 );};
+	/*	PF ( mpdf *par0, mpdf *obs0, epdf *epdf0, int n0 ) :
+				est ( ),_w ( est._w() ),_samples ( est._samples() ),opt_L_smp(false), opt_L_wei(false)
+		{ set_parameters ( par0,obs0,n0 ); set_statistics ( ones ( n0 ),epdf0 ); };*/
+	void set_parameters ( mpdf *par0, mpdf *obs0, int n0, RESAMPLING_METHOD rm=SYSTEMATIC )
+	{ par = par0; obs=obs0; n=n0; resmethod= rm;};
+	void set_statistics ( const vec w0, epdf *epdf0 ) {est.set_statistics ( w0,epdf0 );};
 	//!@}
 	//! Set posterior density by sampling from epdf0
-	void set_est ( const epdf &epdf0 );
-	void set_options(const string &opt){
-		opt_L_res= ( s.find ( "logres" ) !=string::npos );
-		opt_L_smp= ( s.find ( "logsmp" ) !=string::npos );		
+//	void set_est ( const epdf &epdf0 );
+	void set_options ( const string &opt ) {
+		BM::set_options(opt);
+		opt_L_wei= ( opt.find ( "logweights" ) !=string::npos );
+		opt_L_smp= ( opt.find ( "logsamples" ) !=string::npos );
 	}
 	void bayes ( const vec &dt );
@@ -78,7 +82,6 @@
 
 template<class BM_T>
-
 class MPF : public PF {
-	BM_T* Bms[10000];
+	Array<BM_T*> BMs;
 
 	//! internal class for MPDF providing composition of eEmp with external components
@@ -94,8 +97,17 @@
 				Coms ( _w.length() ) {
 		};
+		//! read statistics from MPF
+		void read_statistics ( Array<BM_T*> &A ) {
+			dim = E.dimension() +A ( 0 )->posterior().dimension();
+			for ( int i=0; i<_w.length() ;i++ ) {Coms ( i ) = A ( i )->_e();}
+		}
+		//! needed in resampling
 		void set_elements ( int &i, double wi, const epdf* ep )
 		{_w ( i ) =wi; Coms ( i ) =ep;};
 
-		void set_n ( int n ) {E.set_n ( n ); Coms.set_length ( n );}
+		void set_parameters ( int n ) {
+			E.set_parameters ( n, false );
+			Coms.set_length ( n );
+		}
 		vec mean() const {
 			// ugly
@@ -114,4 +126,33 @@
 			return concat ( E.variance(),pom2-pow ( pom,2 ) );
 		}
+		void qbounds ( vec &lb, vec &ub, double perc=0.95 ) const {
+			//bounds on particles
+			vec lbp;
+			vec ubp;
+			E.qbounds ( lbp,ubp );
+
+			//bounds on Components
+			int dimC=Coms ( 0 )->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<_w.length();i++ ) {
+				// check Coms
+				Coms ( i )->qbounds ( lbc,ubc );
+				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);
+		}
 
 		vec sample() const {it_error ( "Not implemented" );return 0;}
@@ -125,44 +166,40 @@
 	//! Log means of BMs
 	bool opt_L_mea;
-	
+
 public:
 	//! Default constructor.
-	MPF ( mpdf *par0, mpdf *obs0, int n, const BM_T &BMcond0 ) : PF (), jest ( est ) {
-		PF::set_parameters ( par0,obs0,n );
-		jest.set_n ( n );
-		//
-		//TODO test if rv and BMcond.rv are compatible.
-//		rv.add ( rvlin );
-		//
-
-		if ( n>10000 ) {it_error ( "increase 10000 here!" );}
-
-		for ( int i=0;i<n;i++ ) {
-			Bms[i] = new BM_T ( BMcond0 ); //copy constructor
-			const epdf& pom=Bms[i]->posterior();
-			jest.set_elements ( i,1.0/n,&pom );
-		}
+	MPF () : PF (), jest ( est ) {};
+	void set_parameters ( mpdf *par0, mpdf *obs0, int n0, RESAMPLING_METHOD rm=SYSTEMATIC ) {
+		PF::set_parameters ( par0, obs0, n0, rm );
+		jest.set_parameters ( n0 );//duplication of rm
+		BMs.set_length ( n0 );
+	}
+	void set_statistics ( epdf *epdf0, const BM_T* BMcond0 ) {
+
+		PF::set_statistics ( ones ( n ) /n, epdf0 );
+		// copy
+		for ( int i=0;i<n;i++ ) { BMs ( i ) = new BM_T ( *BMcond0 ); BMs ( i )->condition ( _samples ( i ) );}
+
+		jest.read_statistics ( BMs );
+		//options
 	};
-
-	~MPF() {
-	}
 
 	void bayes ( const vec &dt );
 	const epdf& posterior() const {return jest;}
 	const epdf* _e() const {return &jest;} //Fixme: is it useful?
-	//! Set postrior of \c rvc to samples from epdf0. Statistics of Bms are not re-computed! Use only for initialization!
-	void set_est ( const epdf& epdf0 ) {
-		PF::set_est ( epdf0 );  // sample params in condition
-		// copy conditions to BMs
-
-		for ( int i=0;i<n;i++ ) {Bms[i]->condition ( _samples ( i ) );}
-	}
-	void set_options(const string &opt){
-		PF:set_options(opt);
-		opt_L_mea = ( s.find ( "logmeans" ) !=string::npos );
-	}
-	
+	//! Set postrior of \c rvc to samples from epdf0. Statistics of BMs are not re-computed! Use only for initialization!
+	/*	void set_est ( const epdf& epdf0 ) {
+			PF::set_est ( epdf0 );  // sample params in condition
+			// copy conditions to BMs
+
+			for ( int i=0;i<n;i++ ) {BMs(i)->condition ( _samples ( i ) );}
+		}*/
+	void set_options ( const string &opt ) {
+		PF::set_options ( opt );
+		opt_L_mea = ( opt.find ( "logmeans" ) !=string::npos );
+	}
+
 	//!Access function
-	BM* _BM ( int i ) {return Bms[i];}
+	BM* _BM ( int i ) {return BMs ( i );}
 };
 
@@ -180,7 +217,7 @@
 		_samples ( i ) = par->samplecond ( _samples ( i ) );
 		llsP ( i ) = par->_e()->evallog ( _samples ( i ) );
-		Bms[i]->condition ( _samples ( i ) );
-		Bms[i]->bayes ( dt );
-		lls ( i ) = Bms[i]->_ll(); // lls above is also in proposal her must be lls(i) =, not +=!!
+		BMs ( i )->condition ( _samples ( i ) );
+		BMs ( i )->bayes ( dt );
+		lls ( i ) = BMs ( i )->_ll(); // lls above is also in proposal her must be lls(i) =, not +=!!
 		if ( lls ( i ) >mlls ) mlls=lls ( i ); //find maximum likelihood (for numerical stability)
 	}
@@ -204,5 +241,5 @@
 	double eff = 1.0/ ( _w*_w );
 	if ( eff < ( 0.3*n ) ) {
-		ind = est.resample();
+		ind = est.resample ( resmethod );
 		// Resample Bms!
 
@@ -215,7 +252,7 @@
 				// poor-man's solution: replicate constructor here
 				// copied from MPF::MPF
-				delete Bms[i];
-				Bms[i] = new BM_T ( *Bms[ind ( i ) ] ); //copy constructor
-				const epdf& pom=Bms[i]->posterior();
+				delete BMs ( i );
+				BMs ( i ) = new BM_T ( *BMs ( ind ( i ) ) ); //copy constructor
+				const epdf& pom=BMs ( i )->posterior();
 				jest.set_elements ( i,1.0/n,&pom );
 			}
Index: bdm/estim/merger.cpp
===================================================================
--- bdm/estim/merger.cpp (revision 278)
+++ bdm/estim/merger.cpp (revision 283)
@@ -33,5 +33,5 @@
 	it_assert_debug ( rv.equal ( g0->_rv() ),"Incompatible g0" );
 	//Empirical density - samples
-	eSmp.set_parameters ( ones ( Ns ), g0 );
+	eSmp.set_statistics ( ones ( Ns ), g0 );
 	Array<vec> &Smp = eSmp._samples(); //aux
 	vec &w = eSmp._w(); //aux
Index: bdm/estim/merger.h
===================================================================
--- bdm/estim/merger.h (revision 278)
+++ bdm/estim/merger.h (revision 283)
@@ -71,5 +71,5 @@
 	}
 //! Set internal parameters used in approximation
-	void set_parameters ( double beta0, int Ns0, int Nc0 ) {beta=beta0;Ns=Ns0;Nc=Nc0;eSmp.set_n(Ns0,false);}
+	void set_parameters ( double beta0, int Ns0, int Nc0 ) {beta=beta0;Ns=Ns0;Nc=Nc0;eSmp.set_parameters(Ns0,false);}
 //!Initialize the proposal density. This function must be called before merge()!
 	void init() { ////////////// NOT FINISHED
Index: bdm/estim/libKF.cpp
===================================================================
--- bdm/estim/libKF.cpp (revision 279)
+++ bdm/estim/libKF.cpp (revision 283)
@@ -69,5 +69,5 @@
 
 	dimx = pfxu0->_dimx();
-	dimy = phxu0->_dimy();
+	dimy = phxu0->dimension();
 	dimu = pfxu0->_dimu();
 
@@ -169,6 +169,6 @@
 	phxu = phxu0;
 
-	dimx = pfxu0->_dimy();
-	dimy = phxu0->_dimy();
+	dimx = pfxu0->dimension();
+	dimy = phxu0->dimension();
 	dimu = pfxu0->_dimu();
 	// set size of mu - just for constant terms of A and C 
@@ -249,17 +249,3 @@
 }
 
-void KFcondQR::condition ( const vec &QR ) {
-	it_assert_debug ( QR.length() == ( dimx+dimy ),"KFcondRQ: conditioning by incompatible vector" );
-
-	Q.setD ( QR ( 0, dimx-1 ) );
-	R.setD ( QR ( dimx, -1 ) );
-};
-
-void KFcondR::condition ( const vec &R0 ) {
-	it_assert_debug ( R0.length() == ( dimy ),"KFcondR: conditioning by incompatible vector" );
-
-	R.setD ( R0 );
-};
-
-
-}
+}
Index: bdm/estim/arx.cpp
===================================================================
--- bdm/estim/arx.cpp (revision 270)
+++ bdm/estim/arx.cpp (revision 283)
@@ -44,5 +44,5 @@
 }
 
-ARX* ARX::_copy_ ( ) {
+ARX* ARX::_copy_ ( ) const {
 	ARX* Tmp=new ARX ( *this );
 	return Tmp;
Index: bdm/estim/ekf_templ.h
===================================================================
--- bdm/estim/ekf_templ.h (revision 279)
+++ bdm/estim/ekf_templ.h (revision 283)
@@ -19,8 +19,6 @@
 
 //!Extended Kalman filter with unknown \c Q and \c R
-class EKFful_unQR : public EKFfull , public BMcond {
+class EKFful_unQR : public EKFfull  {
 public:
-	//! Default constructor
-	EKFful_unQR ( ) :EKFfull ( ),BMcond ( ) {};
 	void condition ( const vec &QR0 ) {
 		Q=diag(QR0(0,dimx-1));
@@ -30,8 +28,6 @@
 
 //!Extended Kalman filter in Choleski form with unknown \c Q
-class EKFCh_unQ : public EKFCh , public BMcond {
+class EKFCh_unQ : public EKFCh {
 public:
-	//! Default constructor
-	EKFCh_unQ (  ) :EKFCh ( ),BMcond ( ) {};
 	void condition ( const vec &Q0 ) {
 		Q.setD ( Q0,0 );
@@ -42,8 +38,6 @@
 
 //!Extended Kalman filter with unknown parameters in \c IM
-class EKFCh_cond : public EKFCh , public BMcond {
+class EKFCh_cond : public EKFCh  {
 	public:
-	//! Default constructor
-		EKFCh_cond ( ) :EKFCh ( ),BMcond ( ) {};
 		void condition ( const vec &val ) {
 			pfxu->condition ( val );
Index: bdm/estim/libKF.h
===================================================================
--- bdm/estim/libKF.h (revision 279)
+++ bdm/estim/libKF.h (revision 283)
@@ -19,5 +19,5 @@
 #include "../math/chmat.h"
 
-namespace bdm{
+namespace bdm {
 
 /*!
@@ -49,5 +49,5 @@
 	friend std::ostream &operator<< ( std::ostream &os, const KalmanFull &kf );
 	//! For EKFfull;
-	KalmanFull(){};
+	KalmanFull() {};
 };
 
@@ -56,7 +56,7 @@
 * \brief Kalman filter with covariance matrices in square root form.
 
-Parameter evolution model:\f[ x_t = A x_{t-1} + B u_t + Q^{1/2} e_t \f] 
+Parameter evolution model:\f[ x_t = A x_{t-1} + B u_t + Q^{1/2} e_t \f]
 Observation model: \f[ y_t = C x_{t-1} + C u_t + Q^{1/2} w_t. \f]
-Where $e_t$ and $w_t$ are independent vectors Normal(0,1)-distributed disturbances. 
+Where $e_t$ and $w_t$ are independent vectors Normal(0,1)-distributed disturbances.
 */
 template<class sq_T>
@@ -77,5 +77,5 @@
 	mat A;
 	//! Matrix B
-	mat B; 
+	mat B;
 	//! Matrix C
 	mat C;
@@ -112,7 +112,7 @@
 	//! Set estimate values, used e.g. in initialization.
 	void set_est ( const vec &mu0, const sq_T &P0 ) {
-		sq_T pom(dimy);
+		sq_T pom ( dimy );
 		est.set_parameters ( mu0,P0 );
-		P0.mult_sym(C,pom);
+		P0.mult_sym ( C,pom );
 		fy.set_parameters ( C*mu0, pom );
 	};
@@ -133,16 +133,21 @@
 Trivial example:
 \include kalman_simple.cpp
- 
-*/
-class KalmanCh : public Kalman<chmat>{
+
+*/
+class KalmanCh : public Kalman<chmat> {
 protected:
 //! pre array (triangular matrix)
-mat preA;
+	mat preA;
 //! post array (triangular matrix)
-mat postA;
-
-public:
-	//! Default constructor
-	KalmanCh ():Kalman<chmat>(),preA(),postA(){};
+	mat postA;
+
+public:
+	//! copy constructor
+	BM* _copy_() const {
+		KalmanCh* K=new KalmanCh;
+		K->set_parameters ( A,B,C,D,Q,R );
+		K->set_statistics ( _mu,_P );
+		return K;
+	}
 	//! Set parameters with check of relevance
 	void set_parameters ( const mat &A0,const mat &B0,const mat &C0,const mat &D0,const chmat &Q0,const chmat &R0 );
@@ -150,18 +155,18 @@
 		est.set_parameters ( mu0,P0 );
 	};
-	
-	
+
+
 	/*!\brief  Here dt = [yt;ut] of appropriate dimensions
-	
+
 	The following equality hold::\f[
-\left[\begin{array}{cc}
-R^{0.5}\\
-P_{t|t-1}^{0.5}C' & P_{t|t-1}^{0.5}CA'\\
- & Q^{0.5}\end{array}\right]<\mathrm{orth.oper.}>=\left[\begin{array}{cc}
-R_{y}^{0.5} & KA'\\
- & P_{t+1|t}^{0.5}\\
-\\\end{array}\right]\f]
-
-Thus this object evaluates only predictors! Not filtering densities.
+	\left[\begin{array}{cc}
+	R^{0.5}\\
+	P_{t|t-1}^{0.5}C' & P_{t|t-1}^{0.5}CA'\\
+	& Q^{0.5}\end{array}\right]<\mathrm{orth.oper.}>=\left[\begin{array}{cc}
+	R_{y}^{0.5} & KA'\\
+	& P_{t+1|t}^{0.5}\\
+	\\\end{array}\right]\f]
+
+	Thus this object evaluates only predictors! Not filtering densities.
 	*/
 	void bayes ( const vec &dt );
@@ -179,6 +184,6 @@
 	//! Observation Model h(x,u)
 	diffbifn* phxu;
-	
-	enorm<fsqmat> E; 
+
+	enorm<fsqmat> E;
 public:
 	//! Default constructor
@@ -189,9 +194,9 @@
 	void bayes ( const vec &dt );
 	//! set estimates
-	void set_est (vec mu0, mat P0){mu=mu0;P=P0;};
+	void set_est ( vec mu0, mat P0 ) {mu=mu0;P=P0;};
 	//!dummy!
-	const epdf& posterior()const{return E;};
-	const enorm<fsqmat>* _e()const{return &E;};
-	const mat _R(){return P;}
+	const epdf& posterior() const{return E;};
+	const enorm<fsqmat>* _e() const{return &E;};
+	const mat _R() {return P;}
 };
 
@@ -210,4 +215,6 @@
 	//! Default constructor
 	EKF ( RV rvx, RV rvy, RV rvu );
+	//! copy constructor
+	EKF<sq_T>* _copy_() const { return new EKF<sq_T>(this); }
 	//! Set nonlinear functions for mean values and covariance matrices.
 	void set_parameters ( diffbifn* pfxu, diffbifn* phxu, const sq_T Q0, const sq_T R0 );
@@ -223,5 +230,5 @@
 
 class EKFCh : public KalmanCh {
-	protected:
+protected:
 	//! Internal Model f(x,u)
 	diffbifn* pfxu;
@@ -229,4 +236,11 @@
 	diffbifn* phxu;
 public:
+	//! copy constructor duplicated - calls different set_parameters
+	BM* _copy_() const {
+		EKFCh* E=new EKFCh;
+		E->set_parameters ( pfxu,phxu,Q,R );
+		E->set_statistics ( _mu,_P );
+		return E;
+	}
 	//! Set nonlinear functions for mean values and covariance matrices.
 	void set_parameters ( diffbifn* pfxu, diffbifn* phxu, const chmat Q0, const chmat R0 );
@@ -239,24 +253,31 @@
 */
 
-class KFcondQR : public Kalman<ldmat>, public BMcond {
+class KFcondQR : public Kalman<ldmat> {
 //protected:
 public:
+	void condition ( const vec &QR ) {
+		it_assert_debug ( QR.length() == ( dimx+dimy ),"KFcondRQ: conditioning by incompatible vector" );
+
+		Q.setD ( QR ( 0, dimx-1 ) );
+		R.setD ( QR ( dimx, -1 ) );
+	};
+};
+
+/*!
+\brief Kalman Filter with conditional diagonal matrices R and Q.
+*/
+
+class KFcondR : public Kalman<ldmat> {
+//protected:
+public:
 	//!Default constructor
-	KFcondQR ( ) : Kalman<ldmat> ( ),BMcond ( ) {};
-
-	void condition ( const vec &RQ );
-};
-
-/*!
-\brief Kalman Filter with conditional diagonal matrices R and Q.
-*/
-
-class KFcondR : public Kalman<ldmat>, public BMcond {
-//protected:
-public:
-	//!Default constructor
-	KFcondR ( ) : Kalman<ldmat> ( ),BMcond ( ) {};
-
-	void condition ( const vec &R );
+	KFcondR ( ) : Kalman<ldmat> ( ) {};
+
+	void condition ( const vec &R0 ) {
+		it_assert_debug ( R0.length() == ( dimy ),"KFcondR: conditioning by incompatible vector" );
+
+		R.setD ( R0 );
+	};
+
 };
 
@@ -267,6 +288,6 @@
 		dimx ( K0.dimx ), dimy ( K0.dimy ),dimu ( K0.dimu ),
 		A ( K0.A ), B ( K0.B ), C ( K0.C ), D ( K0.D ),
-		Q(K0.Q), R(K0.R),
-		est ( K0.est ), fy ( K0.fy ), _yp(fy._mu()),_Ry(fy._R()), _mu(est._mu()), _P(est._R()) {
+		Q ( K0.Q ), R ( K0.R ),
+		est ( K0.est ), fy ( K0.fy ), _yp ( fy._mu() ),_Ry ( fy._R() ), _mu ( est._mu() ), _P ( est._R() ) {
 
 // copy values in pointers
@@ -279,5 +300,5 @@
 
 template<class sq_T>
-Kalman<sq_T>::Kalman ( ) : BM (), est ( ), fy (),  _yp(fy._mu()), _Ry(fy._R()), _mu(est._mu()), _P(est._R()) {
+Kalman<sq_T>::Kalman ( ) : BM (), est ( ), fy (),  _yp ( fy._mu() ), _Ry ( fy._R() ), _mu ( est._mu() ), _P ( est._R() ) {
 };
 
@@ -287,5 +308,5 @@
 	dimy = C0.rows();
 	dimu = B0.cols();
-	
+
 	it_assert_debug ( A0.cols() ==dimx, "Kalman: A is not square" );
 	it_assert_debug ( B0.rows() ==dimx, "Kalman: B is not compatible" );
@@ -307,5 +328,5 @@
 	it_assert_debug ( dt.length() == ( dimy+dimu ),"KalmanFull::bayes wrong size of dt" );
 
-	sq_T iRy(dimy);
+	sq_T iRy ( dimy );
 	vec u = dt.get ( dimy,dimy+dimu-1 );
 	vec y = dt.get ( 0,dimy-1 );
@@ -328,7 +349,7 @@
 	sq_T pom ( ( int ) Pfull.rows() );
 	iRy.mult_sym_t ( C*Pfull,pom );
-	(_P ) -= pom; // P = P -PC'iRy*CP;
-	(_yp ) = C* _mu  +D*u; //y prediction
-	(_mu ) += _K* ( y- _yp  );
+	( _P ) -= pom; // P = P -PC'iRy*CP;
+	( _yp ) = C* _mu  +D*u; //y prediction
+	( _mu ) += _K* ( y- _yp );
 
 
@@ -340,5 +361,5 @@
 
 };
- 
+
 
 
@@ -369,5 +390,5 @@
 	it_assert_debug ( dt.length() == ( dimy+dimu ),"KalmanFull::bayes wrong size of dt" );
 
-	sq_T iRy(dimy,dimy);
+	sq_T iRy ( dimy,dimy );
 	vec u = dt.get ( dimy,dimy+dimu-1 );
 	vec y = dt.get ( 0,dimy-1 );
@@ -393,5 +414,5 @@
 	sq_T pom ( ( int ) Pfull.rows() );
 	iRy.mult_sym_t ( C*Pfull,pom );
-	(_P ) -= pom; // P = P -PC'iRy*CP;
+	( _P ) -= pom; // P = P -PC'iRy*CP;
 	_yp = phxu->eval ( _mu,u ); //y prediction
 	( _mu ) += _K* ( y-_yp );
Index: bdm/estim/libPF.cpp
===================================================================
--- bdm/estim/libPF.cpp (revision 278)
+++ bdm/estim/libPF.cpp (revision 283)
@@ -1,5 +1,5 @@
 #include "libPF.h"
 
-namespace bdm{
+namespace bdm {
 
 using std::endl;
@@ -13,6 +13,6 @@
 	for ( i=0;i<n;i++ ) {
 		//generate new samples from paramater evolution model;
-		_samples ( i ) = par->samplecond ( _samples ( i ));
-		lls(i )= par->_e()->evallog(_samples(i));
+		_samples ( i ) = par->samplecond ( _samples ( i ) );
+		lls ( i ) = par->_e()->evallog ( _samples ( i ) );
 		lls ( i ) *= obs->evallogcond ( dt,_samples ( i ) );
 
@@ -30,15 +30,16 @@
 	_w ( i ) /=sum; //?
 
-	ind = est.resample();
+	ind = est.resample(resmethod);
 
 }
 
-void PF::set_est ( const epdf &epdf0 ) {
-	int i;
+// void PF::set_est ( const epdf &epdf0 ) {
+// 	int i;
+// 
+// 	for ( i=0;i<n;i++ ) {
+// 		_samples ( i ) = epdf0.sample();
+// 	}
+// }
 
-	for ( i=0;i<n;i++ ) {
-		_samples ( i ) = epdf0.sample();
-	}
-}
 
 }
Index: bdm/stat/libEF.cpp
===================================================================
--- bdm/stat/libEF.cpp (revision 281)
+++ bdm/stat/libEF.cpp (revision 283)
@@ -183,5 +183,5 @@
 };
 
-ivec eEmp::resample ( RESAMPLING_METHOD method ) {
+ivec eEmp::resample (RESAMPLING_METHOD method) {
 	ivec ind=zeros_i ( n );
 	ivec N_babies = zeros_i ( n );
@@ -268,5 +268,5 @@
 }
 
-void eEmp::set_parameters ( const vec &w0, const epdf* epdf0 ) {
+void eEmp::set_statistics ( const vec &w0, const epdf* epdf0 ) {
 	//it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
 	dim = epdf0->dimension();
Index: bdm/stat/libDS.h
===================================================================
--- bdm/stat/libDS.h (revision 271)
+++ bdm/stat/libDS.h (revision 283)
@@ -28,4 +28,5 @@
 */
 class MemDS : public DS {
+	protected:
 	//! internal matrix of data
 	mat Data;
@@ -45,5 +46,33 @@
 	void step();
 	//!Default constructor
+	MemDS () {};
 	MemDS ( mat &Dat, ivec &rowid, ivec &delays );
+};
+
+/*! Read Data Matrix from an IT file
+
+*/
+class FileDS: public MemDS {
+
+public:
+	FileDS ( const string &fname, const string &varname ) :MemDS() {
+		it_file it ( fname );
+		it << Name ( varname ); 
+		it >> Data;
+		time =0;
+		//rowid and delays are ignored
+	}
+	void getdata ( vec &dt ) {
+		it_assert_debug ( dt.length() ==Data.rows(),"" );
+		dt = Data.get_col(time);
+	};
+	void getdata ( vec &dt, const ivec &indeces ){
+		it_assert_debug ( dt.length() ==indeces.length(),"" );
+		vec tmp(indeces.length());
+		tmp = Data.get_col(time);
+		dt = tmp(indeces);
+	};
+	//! returns number of data in the file;
+	int ndat(){return Data.cols();}
 };
 
@@ -96,20 +125,20 @@
 	{ model.set_parameters ( Th0, mu0, sqR0 );};
 	//! Set
-	void set_drv(RV &yrv, RV &urv, RV &rrv){
+	void set_drv ( RV &yrv, RV &urv, RV &rrv ) {
 		Rrv = rrv;
 		Urv = urv;
-		dt_size = yrv._dsize()+urv._dsize();
-		
-		RV drv = concat(yrv,urv);
+		dt_size = yrv._dsize() +urv._dsize();
+
+		RV drv = concat ( yrv,urv );
 		Drv = drv;
 		int td = rrv.mint();
-		H.set_size(drv._dsize()*(-td+1));
-		U.set_size(Urv._dsize());
-		for (int i=-1;i>=td;i--){
-			drv.t(-1);
-			Drv.add(drv); //shift u1
+		H.set_size ( drv._dsize() * ( -td+1 ) );
+		U.set_size ( Urv._dsize() );
+		for ( int i=-1;i>=td;i-- ) {
+			drv.t ( -1 );
+			Drv.add ( drv ); //shift u1
 		}
-		rgrlnk.set_connection(rrv,Drv);
-		
+		rgrlnk.set_connection ( rrv,Drv );
+
 		dtsize = Drv._dsize();
 		utsize = Urv._dsize();
@@ -121,5 +150,5 @@
 	virtual void log_add ( logger &L ) {
 		//DS::log_add ( L ); too long!!
-		L_dt=L.add ( Drv(0,dt_size),"" );
+		L_dt=L.add ( Drv ( 0,dt_size ),"" );
 		L_ut=L.add ( Urv,"" );
 
@@ -131,7 +160,7 @@
 	virtual void logit ( logger &L ) {
 		//DS::logit ( L );
-		L.logit( L_dt, H.left(dt_size));
-		L.logit(L_ut, U);
-		
+		L.logit ( L_dt, H.left ( dt_size ) );
+		L.logit ( L_ut, U );
+
 		mat &A =model._A();
 		mat R =model._R();
Index: bdm/stat/libEF.h
===================================================================
--- bdm/stat/libEF.h (revision 281)
+++ bdm/stat/libEF.h (revision 283)
@@ -94,5 +94,5 @@
 //	virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
 
-	BMEF* _copy_ ( bool changerv=false ) {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
+	BMEF* _copy_ ( bool changerv=false ) const {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
 };
 
@@ -348,4 +348,6 @@
  \f]
 
+Vector \f$\beta\f$ has different meaning (in fact it is 1/beta as used in definition of iG)
+
  Inverse Gamma can be converted to Gamma using \f[
  x\sim iG(a,b) => 1/x\sim G(a,1/b)
@@ -354,30 +356,16 @@
  */
 
-class eigamma : public eEF {
-protected:
-	//!internal egamma
-	egamma eg;
-	//! Vector \f$\alpha\f$
-	vec &alpha;
-	//! Vector \f$\beta\f$ (in fact it is 1/beta as used in definition of iG)
-	vec &beta;
+class eigamma : public egamma {
+protected:
 public :
 	//! \name Constructors
-	//!@{
-	eigamma ( ) :eEF ( ), eg(),alpha ( eg._alpha() ), beta ( eg._beta() ) {};
-	eigamma ( const vec &a, const vec &b ) :eEF ( ), eg(),alpha ( eg._alpha() ), beta ( eg._beta() ) {eg.set_parameters ( a,b );};
-	void set_parameters ( const vec &a, const vec &b ) {eg.set_parameters ( a,b );};
-	//!@}
-
-	vec sample() const {return 1.0/eg.sample();};
-	//! TODO: is it used anywhere?
-//	mat sample ( int N ) const;
-	double evallog ( const vec &val ) const {return eg.evallog ( val );};
-	double lognc () const {return eg.lognc();};
+	//! All constructors are inherited
+	//!@{
+	//!@}
+
+	vec sample() const {return 1.0/egamma::sample();};
 	//! Returns poiter to alpha and beta. Potentially dangerous: use with care!
 	vec mean() const {return elem_div ( beta,alpha-1 );}
 	vec variance() const {vec mea=mean(); return elem_div ( elem_mult ( mea,mea ),alpha-2 );}
-	vec& _alpha() {return alpha;}
-	vec& _beta() {return beta;}
 };
 /*
@@ -490,5 +478,5 @@
 	mgnorm() :mu ( epdf._mu() ) {ep=&epdf;}
 	//!set mean function
-	void set_parameters ( fnc* g0, const sq_T &R0 ) {g=g0; epdf.set_parameters ( zeros ( g->_dimy() ), R0 );}
+	void set_parameters ( fnc* g0, const sq_T &R0 ) {g=g0; epdf.set_parameters ( zeros ( g->dimension() ), R0 );}
 	void condition ( const vec &cond ) {mu=g->eval ( cond );};
 };
@@ -686,9 +674,11 @@
 
 	//! Set samples and weights
-	void set_parameters ( const vec &w0, const epdf* pdf0 );
+	void set_statistics ( const vec &w0, const epdf* pdf0 );
+	//! Set samples and weights
+	void set_statistics ( const epdf* pdf0 , int n ) {set_statistics ( ones ( n ) /n,pdf0 );};
 	//! Set sample
 	void set_samples ( const epdf* pdf0 );
 	//! Set sample
-	void set_n ( int n0, bool copy=true ) {n=n0; w.set_size ( n0,copy );samples.set_size ( n0,copy );};
+	void set_parameters ( int n0, bool copy=true ) {n=n0; w.set_size ( n0,copy );samples.set_size ( n0,copy );};
 	//! Potentially dangerous, use with care.
 	vec& _w()  {return w;};
@@ -700,5 +690,5 @@
 	const Array<vec>& _samples() const {return samples;};
 	//! 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.
-	ivec resample ( RESAMPLING_METHOD method = SYSTEMATIC );
+	ivec resample ( RESAMPLING_METHOD method=SYSTEMATIC );
 	//! inherited operation : NOT implemneted
 	vec sample() const {it_error ( "Not implemented" );return 0;}
@@ -714,4 +704,19 @@
 		for ( int i=0;i<n;i++ ) {pom+=pow ( samples ( i ),2 ) *w ( i );}
 		return pom-pow ( mean(),2 );
+	}
+	//! 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 );}
+			}
+		}
 	}
 };
Index: bdm/stat/libBM.h
===================================================================
--- bdm/stat/libBM.h (revision 281)
+++ bdm/stat/libBM.h (revision 283)
@@ -126,7 +126,7 @@
 	int length() const {return len;} ;
 	int id ( int at ) const{return ids ( at );};
-	int size ( int at ) const {return RV_SIZES ( ids(at) );};
+	int size ( int at ) const {return RV_SIZES ( ids ( at ) );};
 	int time ( int at ) const{return times ( at );};
-	std::string name ( int at ) const {return RV_NAMES ( ids(at) );};
+	std::string name ( int at ) const {return RV_NAMES ( ids ( at ) );};
 	void set_time ( int at, int time0 ) {times ( at ) =time0;};
 	//!@}
@@ -204,5 +204,5 @@
 
 	//! access function
-	int _dimy() const{return dimy;}
+	int dimension() const{return dimy;}
 };
 
@@ -259,4 +259,9 @@
 	//! return expected variance (not covariance!)
 	virtual vec variance() const {it_error ( "not implemneted" );return vec ( 0 );};
+	//! Lower and upper bounds of \c percentage % quantile, returns mean-2*sigma as default
+	virtual void qbounds ( vec &lb, vec &ub, double percentage=0.95 ) const {
+		vec mea=mean(); vec std=sqrt(variance()); 
+		lb = mea-2*std; ub=mea+2*std;
+	};
 	//!@}
 
@@ -592,4 +597,21 @@
 /*! \brief Bayesian Model of a system, i.e. all uncertainty is modeled by probabilities.
 
+This object represents exact or approximate evaluation of the Bayes rule:
+\f[
+f(\theta_t | d_1,\ldots,d_t) = \frac{f(y_t|\theta_t,\cdot) f(\theta_t|d_1,\ldots,d_{t-1})}{f(y_t|d_1,\ldots,d_{t-1})}
+\f]
+
+Access to the resulting posterior density is via function \c posterior().
+
+As a "side-effect" it also evaluates log-likelihood of the data, which can be accessed via function _ll().
+It can also evaluate predictors of future values of \f$y_t\f$, see functions epredictor() and predictor().
+
+Alternatively, it can evaluate posterior density conditioned by a known constant, \f$ c_t \f$:
+\f[
+f(\theta_t | c_t, d_1,\ldots,d_t) \propto  f(y_t,\theta_t|c_t,\cdot, d_1,\ldots,d_{t-1})
+\f]
+
+The value of \f$ c_t \f$ is set by function condition().
+
 */
 
@@ -606,9 +628,9 @@
 	//! @{
 
-	BM () :ll ( 0 ),evalll ( true) {};
+	BM () :ll ( 0 ),evalll ( true ), LIDs ( 3 ), opt_L_bounds ( false ) {};
 	BM ( const BM &B ) :  drv ( B.drv ), ll ( B.ll ), evalll ( B.evalll ) {}
 	//! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually!
-	//! Prototype: \code BM* _copy_(){return new BM(*this);} \endcode
-	virtual BM* _copy_ () {return NULL;};
+	//! Prototype: \code BM* _copy_() const {return new BM(*this);} \endcode
+	virtual BM* _copy_ () const {return NULL;};
 	//!@}
 
@@ -634,4 +656,20 @@
 	//!@}
 
+	//! \name Extension to conditional BM
+	//! This extension is useful e.g. in Marginalized Particle Filter (\ref bdm::MPF).
+	//! Alternatively, it can be used for automated connection to DS when the condition is observed
+	//!@{
+
+	//! Name of extension variable
+	RV rvc;
+	//! access function
+	const RV& _rvc() const {return rvc;}
+
+	//! Substitute \c val for \c rvc.
+	virtual void condition ( const vec &val ) {it_error ( "Not implemented!" );};
+
+	//!@}
+
+
 	//! \name Access to attributes
 	//!@{
@@ -646,30 +684,40 @@
 	//!@}
 
-};
-
-/*!
-\brief Conditional Bayesian Filter
-
-Evaluates conditional filtering density \f$f(rv|rvc,data)\f$ for a given \c rvc which is specified in each step by calling function \c condition.
-
-This is an interface class used to assure that certain BM has operation \c condition .
-
-*/
-
-class BMcond :public bdmroot {
-protected:
-	//!dimension of the conditioning variable
-	int dimc;
-	//! Identificator of the conditioning variable
-	RV rvc;
-public:
-	//! Substitute \c val for \c rvc.
-	virtual void condition ( const vec &val ) =0;
-	//! Default constructor
-	BMcond ( ) :rvc ( ) {};
-	//! Destructor for future use
-	virtual ~BMcond() {};
-	//! access function
-	const RV& _rvc() const {return rvc;}
+	//! \name Logging of results
+	//!@{
+
+	//! Set boolean options from a string
+	void set_options ( const string &opt ) {
+		opt_L_bounds= ( opt.find ( "logbounds" ) !=string::npos );
+	}
+	//! IDs of storages in loggers
+	ivec LIDs;
+
+	//! Option for logging bounds
+	bool opt_L_bounds;
+	//! Add all logged variables to a logger
+	void log_add ( logger *L, const string &name="" ) {
+		// internal
+		RV r;
+		if ( posterior().isnamed() ) {r=posterior()._rv();}
+		else{r=RV ( "est", posterior().dimension() );};
+
+		// Add mean value
+		LIDs ( 0 ) =L->add ( r,name );
+		if ( opt_L_bounds ) {
+			LIDs ( 1 ) =L->add ( r,name+"_lb" );
+			LIDs ( 2 ) =L->add ( r,name+"_ub" );
+		}
+	}
+	void logit ( logger *L ) {
+		L->logit ( LIDs ( 0 ), posterior().mean() );
+		if ( opt_L_bounds ) {
+			vec ub,lb;
+			posterior().qbounds(lb,ub);
+			L->logit ( LIDs ( 1 ), lb );
+			L->logit ( LIDs ( 2 ), ub );
+		}
+	}
+	//!@}
 };
 
