Index: bdm/estim/merger.cpp
===================================================================
--- bdm/estim/merger.cpp (revision 176)
+++ bdm/estim/merger.cpp (revision 176)
@@ -0,0 +1,62 @@
+#include <itpp/itbase.h>
+#include "merger.h"
+#include "arx.h"
+
+vec merger::lognorm_merge ( mat &lW ) {
+	int nu=lW.rows();
+	vec mu = sum ( lW ) /nu; //mean of logs
+	vec lam = sum ( pow ( lW,2 ) )-nu*pow ( mu,2 );
+	double coef=0.0;
+	switch ( nu ) {
+		case 2:
+			coef=sqrt ( beta*2 ) * ( 1-0.5*sqrt ( ( 4*beta-3 ) /beta ) );
+			return exp ( coef*sqrt ( lam ) + mu );
+			break;
+		case 3://Ration of Bessel
+			break;
+		case 4:
+			break;
+	}
+	return vec(0);
+}
+
+void merger::merge ( const epdf* g0 ) {
+	mat Smp ( rv.count(),Ns );
+	mat lW=zeros ( n,Ns );
+
+	vec lw_src ( Ns );
+	vec lw_mix ( Ns );
+	vec vec0 ( 0 );
+
+	//Make sure g0 is compatible
+	it_assert_debug ( rv.equal(g0->_rv()),"Incompatible g0" );
+
+	// Initial component in the mixture model from g0
+	mat V0=1e-8*eye ( rv.count() +1 );
+	ARX A0 ( rv, V0, 3.0 ); //initial guess of Mix: zero mean, large variance
+
+	//generate samples from the proposal g0
+	for ( int i =0; i<Ns; i++ ) {Smp.set_col ( i,g0->sample() );}
+
+	//Initialize Mixture model
+	Mix.init ( &A0, Smp, Nc );
+
+	bool converged=false;
+	while ( !converged ) {
+		Mix.bayesB ( Smp );
+		//Generate weight for each particle
+		for ( int i=0;i<n;i++ ) {
+			//Split according to dependency in rvs
+			//======== Same RVs ===========
+			if ( rvsinrv ( i ).length() ==rv.count() ) {
+				// no need for conditioning or marginalization
+				lw_src=mpdfs ( i )->_epdf().evalpdflog_m ( Smp );
+				lw_mix=Mix.logpred_m ( Smp );
+			}
+			lW.set_row ( i, lw_src-lw_mix );
+		}
+		converged = true;
+	}
+
+}
+
Index: bdm/estim/merger.h
===================================================================
--- bdm/estim/merger.h (revision 176)
+++ bdm/estim/merger.h (revision 176)
@@ -0,0 +1,70 @@
+/*!
+  \file
+  \brief Mergers for combination of pdfs
+  \author Vaclav Smidl.
+
+  -----------------------------------
+  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
+
+  Using IT++ for numerical operations
+  -----------------------------------
+*/
+
+#ifndef MER_H
+#define MER_H
+
+#include <itpp/itbase.h>
+#include "mixef.h"
+
+using namespace itpp;
+using std::string;
+
+/*!
+@brief Function for general combination of pdfs
+
+Mixtures of Gaussian densities are used internally. Switching to other densities should be trivial.
+*/
+
+class merger : public compositepdf, public epdf {
+protected:
+	//!Internal mixture of EF models
+	MixEF Mix;
+	//!Number of samples used in approximation
+	int Ns;
+	//!Number of components in a mixture
+	int Nc;
+	//!Prior on the log-normal merging model
+	double beta;
+public:
+//!Default constructor
+	merger ( const Array<mpdf*> &S ) :
+			compositepdf ( S ), epdf ( getrv ( false ) ),
+			Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) )
+	{ beta=2.0; Ns=100; Nc=10;}
+	//! Set internal parameters used in approximation
+	void setparameters ( double beta0, int Ns0, int Nc0 ) {	beta=beta0;Ns=Ns0;Nc=Nc0;}
+	//!Initialize the proposal density. This function must be called before merge()!
+	void init() {
+		Array<vec> Smps ( n );
+		//Gibbs sampling
+		for ( int i=0;i<n;i++ ) {Smps ( i ) =zeros ( 0 );}
+	}
+	//!Create a mixture density using known proposal
+	void merge ( const epdf* g0 );
+	//!Create a mixture density, make sure to call init() before the first call
+	void merge () {merge ( & ( Mix._epdf() ) );};
+
+	//! Merge log-likelihood values
+	vec lognorm_merge ( mat &lW );
+	//! sample from merged density
+	//! weight w is a
+	vec sample ( )const { return Mix._epdf().sample();}
+	double evalpdflog ( const vec &dt ) const{ return Mix._epdf().evalpdflog ( dt );}
+	vec mean()const {return Mix._epdf().mean();}
+	//! for future use
+	virtual ~merger() {};
+};
+
+
+
+#endif // MER_H
Index: bdm/estim/mixef.cpp
===================================================================
--- bdm/estim/mixef.cpp (revision 176)
+++ bdm/estim/mixef.cpp (revision 176)
@@ -0,0 +1,111 @@
+#include "mixef.h"
+#include <vector>
+
+using namespace itpp;
+
+
+void MixEF::init ( BMEF* Com0, const mat &Data, int c ) {
+	//prepare sizes
+	Coms.set_size ( c );
+	n=c;
+	weights.set_parameters ( 1/ ( double ) c*ones ( c ) ); //assume at least one observation in each comp.
+	//est will be done at the end
+	//
+	int i;
+	int ndat = Data.cols();
+	//Estimate  Com0 from all data
+	Coms ( 0 ) = ( BMEF* ) Com0->_copy_();
+//	Coms(0)->set_evalll(false);
+	Coms ( 0 )->bayesB ( Data );
+	// Flatten it to its original shape
+	Coms ( 0 )->flatten ( Com0 );
+
+	//Copy it to the rest 
+	for ( i=1;i<n;i++ ) {
+		//copy Com0 and create new rvs for them
+		Coms ( i ) = ( BMEF* ) Coms ( 0 )->_copy_ ( true );
+	}
+	//Pick some data for each component and update it
+	for ( i=0;i<n;i++ ) {
+		//pick one datum
+		int ind=ndat*UniRNG.sample();
+		Coms ( i )->bayes ( Data.get_col ( ind ),ndat/n );
+	}
+
+	//est already exists - must be deleted before build_est() can be used
+	delete est;
+	build_est();
+
+}
+void MixEF::bayesB ( const mat &Data ) {
+	this->bayes ( Data );
+}
+
+void MixEF::bayes ( const vec &data ) {
+
+};
+
+void MixEF::bayes ( const mat &data ) {
+	int ndat=data.cols();
+	int t,i,niter;
+	bool converged;
+
+	multiBM weights0 ( weights );
+
+	Array<BMEF*> Coms0 ( n );
+	for ( i=0;i<n;i++ ) {Coms0 ( i ) = ( BMEF* ) Coms ( i )->_copy_();}
+
+	niter=0;
+	mat W=ones ( n,ndat ) / n;
+	mat Wlast=ones ( n,ndat ) / n;
+	vec w ( n );
+	vec ll ( n );
+	// tmp for weights
+	vec wtmp = zeros ( n );
+	//Estim
+	while ( !converged ) {
+		// Copy components back to their initial values
+		// All necessary information is now in w and Coms0.
+		Wlast = W;
+		//
+		for ( t=0;t<ndat;t++ ) {
+			for ( i=0;i<n;i++ ) {
+				ll ( i ) =Coms ( i )->logpred ( data.get_col ( t ) );
+				wtmp =0.0; wtmp ( i ) =1.0;
+				ll ( i ) += weights.logpred ( wtmp );
+			}
+			w = exp ( ll-max ( ll ) );
+			W.set_col ( t, w/sum ( w ) );
+		}
+
+		for ( i=0;i<n;i++ ) {
+			Coms ( i )-> set_statistics ( Coms0 ( i ) );
+		}
+		weights.set_statistics ( &weights0 );
+
+		for ( t=0;t<ndat;t++ ) {
+			for ( i=0;i<n;i++ ) {
+				Coms ( i )-> bayes ( data.get_col ( t ),W ( i,t ) );
+			}
+			weights.bayes ( W.get_col ( t ) );
+		}
+
+		niter++;
+		//TODO better convergence rule.
+		converged = ( sumsum ( abs ( W-Wlast ) ) /n<0.001 );
+	}
+
+	//Clean Coms0
+	for ( i=0;i<n;i++ ) {delete Coms0 ( i );}
+};
+
+
+double MixEF::logpred ( const vec &dt ) const {
+
+	vec w=weights._epdf().mean();
+	double exLL=0.0;
+	for ( int i=0;i<n;i++ ) {
+		exLL+=w ( i ) *exp ( Coms ( i )->logpred ( dt ) );
+	}
+	return log ( exLL );
+}
Index: bdm/estim/mixef.h
===================================================================
--- bdm/estim/mixef.h (revision 176)
+++ bdm/estim/mixef.h (revision 176)
@@ -0,0 +1,98 @@
+/*!
+  \file
+  \brief Bayesian Filtering for mixtures of exponential family (EF) members
+  \author Vaclav Smidl.
+
+  -----------------------------------
+  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
+
+  Using IT++ for numerical operations
+  -----------------------------------
+*/
+
+#ifndef MEF_H
+#define MEF_H
+
+#include <itpp/itbase.h>
+#include "../stat/libFN.h"
+#include "../stat/libEF.h"
+#include "../stat/emix.h"
+
+using namespace itpp;
+
+/*!
+* \brief Mixture of Exponential Family Densities
+
+An approximate estimation method for models with latent discrete variable, such as
+mixture models of the following kind:
+\f[
+f(y_t|\psi_t, \Theta) = \sum_{i=1}^{n} w_i f(y_t|\psi_t, \theta_i)
+\f]
+where  \f$\psi\f$ is a known function of past outputs, \f$w=[w_1,\ldots,w_n]\f$ are component weights, and component parameters \f$\theta_i\f$ are assumed to be mutually independent. \f$\Theta\f$ is an aggregation af all component parameters and weights, i.e. \f$\Theta = [\theta_1,\ldots,\theta_n,w]\f$.
+
+The characteristic feature of this model is that if the exact values of the latent variable were known, estimation of the parameters can be handled by a single model. For example, for the case of mixture models, posterior density for each component parameters would be a BayesianModel from Exponential Family.
+
+This class uses EM-style type algorithms for estimation of its parameters. Under this simplification, the posterior density is a product of exponential family members, hence approximate estimation project this class itself belongs to the exponential family.
+
+TODO: Extend BM to use rvc.
+*/
+class MixEF: public BM {
+protected:
+	//!Number of components
+	int n;
+	//! Models for Components of \f$\theta_i\f$
+	Array<BMEF*> Coms;
+	//! Statistics for weights
+	multiBM weights;
+	//!Posterior on component parameters
+	eprod* est;
+	////!Indeces of component rvc in common rvc
+
+	//! Auxiliary function for use in constructors
+	void build_est() {
+		Array<const epdf*> epdfs ( n+1 );
+		for ( int i=0;i<Coms.length();i++ ) {
+//			it_assert_debug(!x,"MixEF::MixEF : Incompatible components");
+			epdfs ( i ) =& ( Coms ( i )->_epdf() );
+		}
+		// last in the product is the weight
+		epdfs ( n ) =& ( weights._epdf() );
+		est = new eprod ( epdfs );
+	}
+
+public:
+	//! Full constructor
+	MixEF ( const Array<BMEF*> &Coms0, const vec &alpha0 ) :
+			BM ( RV() ), n ( Coms0.length() ), Coms ( n ),
+			weights ( RV ( "{w }", vec_1 ( n ) ),alpha0 ) {
+		it_assert_debug ( n>0,"MixEF::MixEF : Empty Component list" );
+
+		for ( int i=0;i<n;i++ ) {Coms ( i ) = ( BMEF* ) Coms0 ( i )->_copy_();}
+		build_est();
+	};
+	MixEF () :
+			BM ( RV() ), n ( 0 ), Coms ( 0 ),
+			weights ( RV ( "{w }", vec_1 ( 0 ) ),vec ( 0 ) ) {build_est();}
+	//! Initializing the mixture by a random pick of centroids from data
+	//! \param Com0 Initial component - necessary to determine its type.
+	//! \param Data Data on which the initialization will be done
+	//! \param c Initial number of components, default=5
+	void init ( BMEF* Com0, const mat &Data, int c=5 );
+	//Destructor
+	~MixEF() {
+		delete est;
+		for ( int i=0;i<n;i++ ) {delete Coms ( i );}
+	}
+	//! Recursive EM-like algorithm (QB-variant), see Karny et. al, 2006
+	void bayes ( const vec &dt );
+	//! EM algorithm
+	void bayes ( const mat &dt );
+	void bayesB ( const mat &dt );
+	double logpred ( const vec &dt ) const;
+	const epdf& _epdf() const {return *est;}
+};
+
+
+#endif // MEF_H
+
+
Index: bdm/stat/libBM.cpp
===================================================================
--- bdm/stat/libBM.cpp (revision 175)
+++ bdm/stat/libBM.cpp (revision 176)
@@ -91,5 +91,5 @@
 // }
 
-RV RV::subselect ( ivec ind ) const {
+RV RV::subselect (const ivec &ind ) const {
 	RV ret;
 	ret.init ( ids ( ind ), names ( to_Arr ( ind ) ), sizes ( ind ), times ( ind ) );
@@ -99,5 +99,5 @@
 void RV::t ( int delta ) { times += delta;}
 
-RV RV::operator() ( ivec ind ) const {
+RV RV::operator() (const ivec &ind ) const {
 	RV ret;
 	if ( ind.length() >0 ) {
@@ -140,5 +140,5 @@
 }
 
-ivec RV::dataind ( RV rv2 ) const {
+ivec RV::dataind (const RV &rv2 ) const {
 	ivec res ( 0 );
 	if ( rv2.count()>0 ) {
@@ -155,5 +155,5 @@
 }
 
-RV RV::subt ( const RV rv2 ) const {
+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
Index: bdm/stat/libEF.h
===================================================================
--- bdm/stat/libEF.h (revision 173)
+++ bdm/stat/libEF.h (revision 176)
@@ -91,5 +91,5 @@
 	void bayes ( const vec &dt );
 	//!Flatten the posterior
-	virtual void flatten ( BMEF * B) {it_error ( "Not implemented" );}
+	virtual void flatten ( BMEF * B ) {it_error ( "Not implemented" );}
 };
 
@@ -206,5 +206,5 @@
 	vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
 	vec mean() const {return beta/sum ( beta );};
-	//! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
+	//! In this instance, val is ...
 	double evalpdflog_nn ( const vec &val ) const {return ( beta-1 ) *log ( val );};
 	double lognc () const {
@@ -216,4 +216,12 @@
 	//!access function
 	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());
+		}
+		beta= beta0;
+	}
 };
 
@@ -239,5 +247,5 @@
 		eDirich pred ( est );
 		vec &beta = pred._beta();
-		
+
 		double lll;
 		if ( frg<1.0 )
@@ -250,13 +258,17 @@
 		return pred.lognc()-lll;
 	}
-	void flatten (BMEF* B ) {
-		eDirich* E=dynamic_cast<eDirich*>(B);
+	void flatten ( BMEF* B ) {
+		eDirich* E=dynamic_cast<eDirich*> ( B );
 		// sum(beta) should be equal to sum(B.beta)
 		const vec &Eb=E->_beta();
-		est.pow ( sum(beta)/sum(Eb) );
+		est.pow ( sum ( beta ) /sum ( Eb ) );
+		if ( evalll ) {last_lognc=est.lognc();}
+	}
+	const epdf& _epdf() const {return est;};
+	void set_parameters ( const vec &beta0 ) {
+		est.set_parameters(beta0);
+		rv = est._rv();
 		if(evalll){last_lognc=est.lognc();}
 	}
-	const epdf& _epdf() const {return est;};
-	//!access funct
 };
 
Index: bdm/stat/emix.h
===================================================================
--- bdm/stat/emix.h (revision 175)
+++ bdm/stat/emix.h (revision 176)
@@ -78,5 +78,5 @@
 	Additional parameter overlap is left for future use. Do not set to true for mprod.
 	*/
-	mprod ( Array<mpdf*> mFacs): compositepdf(mFacs), mpdf(getrv(true),RV()), rvcinds(n), epdfs(n)
+	mprod ( Array<mpdf*> mFacs): compositepdf(mFacs), mpdf(getrv(true),RV()), epdfs(n), rvcinds(n)
 	{ 	setrvc(rv,rvc);
 		setrvcinrv(rvc,rvcinds);
Index: bdm/stat/libBM.h
===================================================================
--- bdm/stat/libBM.h (revision 175)
+++ bdm/stat/libBM.h (revision 176)
@@ -80,9 +80,9 @@
 	bool add ( const RV &rv2 );
 	//! Subtract  another variable from the current one
-	RV subt ( const RV rv2 ) const;
+	RV subt ( const RV &rv2 ) const;
 	//! Select only variables at indeces ind
-	RV subselect ( ivec ind ) const;
+	RV subselect ( const ivec &ind ) const;
 	//! Select only variables at indeces ind
-	RV operator() ( ivec ind ) const;
+	RV operator() ( const ivec &ind ) const;
 	//! Shift \c time shifted by delta.
 	void t ( int delta );
@@ -90,5 +90,5 @@
 	str tostr() const;
 	//! generate indeces into \param crv data vector that form data vector of self.
-	ivec dataind(RV crv) const;
+	ivec dataind(const RV &crv) const;
 
 	//!access function
@@ -103,4 +103,12 @@
 	//!access function
 	std::string name ( int at ) {return names ( at );};
+	
+	//!access function
+	void set_id ( int at, int id0 ) {ids ( at )=id0;};
+	//!access function
+	void set_size ( int at, int size0 ) {sizes ( at )=size0; tsize=sum(sizes);};
+	//!access function
+	void set_time ( int at, int time0 ) {times ( at )=time0;};
+
 	//!Assign unused ids to this rv 
 	void newids();
@@ -160,5 +168,5 @@
 
 	//! Compute log-probability of multiple values argument \c val
-	virtual vec evalpdflog ( const mat &Val ) const {
+	virtual vec evalpdflog_m ( const mat &Val ) const {
 		vec x ( Val.cols() );
 		for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evalpdflog( Val.get_col(i) ) ;}
@@ -318,4 +326,6 @@
 	//! I.e. marginal likelihood of the data with the posterior integrated out.
 	virtual double logpred(const vec &dt)const{it_error("Not implemented");return 0.0;}
+	//! Matrix version of logpred
+	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;}
 	
 	//! Destructor for future use;
Index: bdm/stat/merger.cpp
===================================================================
--- bdm/stat/merger.cpp (revision 163)
+++  (revision )
@@ -1,3 +1,0 @@
-#include <itpp/itbase.h>
-#include "merger.h"
-
Index: bdm/stat/merger.h
===================================================================
--- bdm/stat/merger.h (revision 168)
+++  (revision )
@@ -1,108 +1,0 @@
-/*!
-  \file
-  \brief Mergers for combination of pdfs
-  \author Vaclav Smidl.
-
-  -----------------------------------
-  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
-
-  Using IT++ for numerical operations
-  -----------------------------------
-*/
-
-#ifndef MER_H
-#define MER_H
-
-#include <itpp/itbase.h>
-#include "libBM.h"
-#include "libEF.h"
-#include "emix.h"
-
-using namespace itpp;
-using std::string;
-
-/*!
-@brief Function for general combination of pdfs 
-
-Parent of simple arithmentic or geometric combinations
-*/
-
-class merger : public mprod{
-protected:	
-	//! Additional pdf on the part in condition (if undefined);
-	enorm<fsqmat> condpdf;
-	//! Find potential overlaps in rv
-	Array<ivec> overlaps;
-public:
-//!Default constructor
-	merger (const Array<mpdf*> &in_sources) : mprod(in_sources,true), condpdf(rvc), overlaps(rv.count()) {
-		if (rvc.count()>0) {
-			vec mu = zeros(rvc.count());
-			mat R = 100*eye(rvc.count());
-			condpdf.set_parameters(mu,R);
-		}
-		for (int i=0;i<rv.count();i++){ // cycle over rv
-			ivec tmp(0);
-			for (int j=0;j<n;j++){//cycle over rvinds
-				ivec id=itpp::find(rvinds(j)==i);
-				tmp=concat(tmp,id); // add i to tmp if found
-			}
-			overlaps(i)=tmp;
-		}
-	}
-	//! sample from merged density
-	//! weight w is a 
-/*	vec sample(double &w, const vec &smp0){
-		// result
-		vec smp=smp0;
-		vec cond=condpdf.sample(); // Just like in samplecond, here it is not given!
-		// temporary
-		mat smpi=zeros(rv.count()+rvc.count(), n);
-		vec condi;
-		vec ll(n);
-		
-		// Consider arithmetic mean as a proposal density
-		ll = 0;
-		for ( int i = ( n - 1 );i >= 0;i-- ) {
-			if ( rvcinds ( i ).length() > 0 ) {
-				condi = zeros ( rvcinrv.length() + rvcinds.length() );
-				// copy data in condition
-				set_subvector ( condi,rvcinds ( i ), cond );
-				// copy data in already generated sample
-				set_subvector ( condi,rvcinrv ( i ), smp );
-
-				mpdfs ( i )->condition ( condi );
-			}
-			smpi.set_col(i, epdfs ( i )->sample());
-			// add ith likelihood contribution
-			ll+=epdfs(i)->evalpdflog(smpi);
-		}
-		// Now lets analyze samples smpi
-		for (int i=0;i<rv.count();i++){
-			// number of components
-			int noc = overlaps(i).length();
-			if (noc<2) { // only one sample in this dimension, take it
-				smp(i) = smpi(i,overlaps(i)(0));
-			}
-			else {
-				//pick dimension
-				double unis = UniRNG.sample();
-				dim = 0;
-				while(unis<(dim+1)/noc) {dim++;} // 
-				//in dim we now have randomly picked dimension
-				smp(i) = smpi(i,overlaps(i)(dim));
-			}
-		}
-		
-		// now we have sample how to evaluate weight?
-		return smp;
-		
-		// copied from mprod.sample
-	};*/
-//	project 
-	//! for future use
-	virtual ~merger() {};
-};
-
-
-#endif // MER_H
Index: bdm/CMakeLists.txt
===================================================================
--- bdm/CMakeLists.txt (revision 168)
+++ bdm/CMakeLists.txt (revision 176)
@@ -2,6 +2,6 @@
 
 SET(BdmMath math/libDC.cpp math/libDC.h math/chmat.cpp math/chmat.h)
-SET(BdmStat stat/libDS.cpp stat/libDS.h stat/libFN.cpp stat/libFN.h stat/libBM.cpp stat/libBM.h stat/libEF.cpp stat/libEF.h stat/loggers.cpp stat/loggers.h stat/emix.cpp stat/emix.h stat/merger.cpp)
-SET(BdmEstim estim/libKF.cpp estim/libKF.h estim/libPF.cpp estim/libPF.h estim/arx.cpp estim/arx.h estim/mixef.cpp estim/mixef.h)
+SET(BdmStat stat/libDS.cpp stat/libDS.h stat/libFN.cpp stat/libFN.h stat/libBM.cpp stat/libBM.h stat/libEF.cpp stat/libEF.h stat/loggers.cpp stat/loggers.h stat/emix.cpp stat/emix.h)
+SET(BdmEstim estim/libKF.cpp estim/libKF.h estim/libPF.cpp estim/libPF.h estim/arx.cpp estim/arx.h estim/mixef.cpp estim/mixef.h estim/merger.cpp)
 SET(BdmUI userinfo.cpp userinfo.h)
 
