Index: library/bdm/stat/emix.cpp
===================================================================
--- library/bdm/stat/emix.cpp (revision 878)
+++ library/bdm/stat/emix.cpp (revision 886)
@@ -3,8 +3,8 @@
 namespace bdm {
 
-void emix::validate (){
-	bdm_assert ( Coms.length() > 0, "There has to be at least one component." );
-
-	bdm_assert ( Coms.length() == w.length(), "It is obligatory to define weights of all the components." );
+void emix_base::validate (){
+	bdm_assert ( no_coms() > 0, "There has to be at least one component." );
+
+	bdm_assert ( no_coms() == w.length(), "It is obligatory to define weights of all the components." );
 
 	double sum_w = sum ( w );
@@ -12,10 +12,10 @@
 	w = w / sum_w;
 
-	dim = Coms ( 0 )->dimension();
-	RV rv_tmp = Coms ( 0 )->_rv() ;
-	bool isnamed = Coms ( 0 )->isnamed();
-	for ( int i = 1; i < Coms.length(); i++ ) {
-		bdm_assert ( dim == ( Coms ( i )->dimension() ), "Component sizes do not match!" );
-		isnamed &= Coms(i)->isnamed() & Coms(i)->_rv().equal(rv_tmp);
+	dim = component ( 0 )->dimension();
+	RV rv_tmp = component ( 0 )->_rv() ;
+	bool isnamed = component( 0 )->isnamed();
+	for ( int i = 1; i < no_coms(); i++ ) {
+		bdm_assert ( dim == ( component ( i )->dimension() ), "Component sizes do not match!" );
+		isnamed &= component(i)->isnamed() & component(i)->_rv().equal(rv_tmp);
 	}
 	if (isnamed)
@@ -23,16 +23,7 @@
 }
 
-void emix::from_setting ( const Setting &set ) {
-	UI::get ( Coms, set, "pdfs", UI::compulsory );
-
-	if ( !UI::get ( w, set, "weights", UI::optional ) ) {
-		int len = Coms.length();
-		w.set_length ( len );
-		w = 1.0 / len;
-	}
-}
-
-
-vec emix::sample() const {
+
+
+vec emix_base::sample() const {
 	//Sample which component
 	vec cumDist = cumsum ( w );
@@ -46,21 +37,21 @@
 	}
 
-	return Coms ( i )->sample();
-}
-
-vec emix::mean() const {
+	return component ( i )->sample();
+}
+
+vec emix_base::mean() const {
 	int i;
 	vec mu = zeros ( dim );
 	for ( i = 0; i < w.length(); i++ ) {
-		mu += w ( i ) * Coms ( i )->mean();
+		mu += w ( i ) * component ( i )->mean();
 	}
 	return mu;
 }
 
-vec emix::variance() const {
+vec emix_base::variance() const {
 	//non-central moment
 	vec mom2 = zeros ( dim );
 	for ( int i = 0; i < w.length(); i++ ) {
-		mom2 += w ( i ) * ( Coms ( i )->variance() + pow ( Coms ( i )->mean(), 2 ) );
+		mom2 += w ( i ) * ( component( i )->variance() + pow ( component ( i )->mean(), 2 ) );
 	}
 	//central moment
@@ -68,9 +59,9 @@
 }
 
-double emix::evallog ( const vec &val ) const {
+double emix_base::evallog ( const vec &val ) const {
 	int i;
 	double sum = 0.0;
 	for ( i = 0; i < w.length(); i++ ) {
-		sum += w ( i ) * exp ( Coms ( i )->evallog ( val ) );
+		sum += w ( i ) * exp ( component ( i )->evallog ( val ) );
 	}
 	if ( sum == 0.0 ) {
@@ -82,21 +73,21 @@
 }
 
-vec emix::evallog_mat ( const mat &Val ) const {
+vec emix_base::evallog_mat ( const mat &Val ) const {
 	vec x = zeros ( Val.cols() );
 	for ( int i = 0; i < w.length(); i++ ) {
-		x += w ( i ) * exp ( Coms ( i )->evallog_mat ( Val ) );
+		x += w ( i ) * exp ( component( i )->evallog_mat ( Val ) );
 	}
 	return log ( x );
 };
 
-mat emix::evallog_coms ( const mat &Val ) const {
+mat emix_base::evallog_coms ( const mat &Val ) const {
 	mat X ( w.length(), Val.cols() );
 	for ( int i = 0; i < w.length(); i++ ) {
-		X.set_row ( i, w ( i ) *exp ( Coms ( i )->evallog_mat ( Val ) ) );
+		X.set_row ( i, w ( i ) *exp ( component( i )->evallog_mat ( Val ) ) );
 	}
 	return X;
 }
 
-shared_ptr<epdf> emix::marginal ( const RV &rv ) const {
+shared_ptr<epdf> emix_base::marginal ( const RV &rv ) const {
 	emix *tmp = new emix();
 	shared_ptr<epdf> narrow ( tmp );
@@ -105,10 +96,10 @@
 }
 
-void emix::marginal ( const RV &rv, emix &target ) const {
+void emix_base::marginal ( const RV &rv, emix &target ) const {
 	bdm_assert ( isnamed(), "rvs are not assigned" );
 
-	Array<shared_ptr<epdf> > Cn ( Coms.length() );
-	for ( int i = 0; i < Coms.length(); i++ ) {
-		Cn ( i ) = Coms ( i )->marginal ( rv );
+	Array<shared_ptr<epdf> > Cn ( no_coms() );
+	for ( int i = 0; i < no_coms(); i++ ) {
+		Cn ( i ) = component ( i )->marginal ( rv );
 	}
 
@@ -118,5 +109,5 @@
 }
 
-shared_ptr<pdf> emix::condition ( const RV &rv ) const {
+shared_ptr<pdf> emix_base::condition ( const RV &rv ) const {
 	bdm_assert ( isnamed(), "rvs are not assigned" );
 	mratio *tmp = new mratio ( this, rv );
@@ -124,151 +115,14 @@
 }
 
-void egiwmix::set_parameters ( const vec &w0, const Array<egiw*> &Coms0, bool copy ) {
-	w = w0 / sum ( w0 );
-	int i;
-	for ( i = 0; i < w.length(); i++ ) {
-		bdm_assert_debug ( dim == ( Coms0 ( i )->dimension() ), "Component sizes do not match!" );
-	}
-	if ( copy ) {
-		Coms.set_length ( Coms0.length() );
-		for ( i = 0; i < w.length(); i++ ) {
-			bdm_error ( "Not implemented" );
-			// *Coms ( i ) = *Coms0 ( i );
-		}
-		destroyComs = true;
-	} else {
-		Coms = Coms0;
-		destroyComs = false;
-	}
-}
-
-void 	egiwmix::validate (){
-   dim = Coms ( 0 )->dimension();
-}
-
-vec egiwmix::sample() const {
-	//Sample which component
-	vec cumDist = cumsum ( w );
-	double u0;
-#pragma omp critical
-	u0 = UniRNG.sample();
-
-	int i = 0;
-	while ( ( cumDist ( i ) < u0 ) && ( i < ( w.length() - 1 ) ) ) {
-		i++;
-	}
-
-	return Coms ( i )->sample();
-}
-
-vec egiwmix::mean() const {
-	int i;
-	vec mu = zeros ( dim );
-	for ( i = 0; i < w.length(); i++ ) {
-		mu += w ( i ) * Coms ( i )->mean();
-	}
-	return mu;
-}
-
-vec egiwmix::variance() const {
-	// non-central moment
-	vec mom2 = zeros ( dim );
-	for ( int i = 0; i < w.length(); i++ ) {
-		// pow is overloaded, we have to use another approach
-		mom2 += w ( i ) * ( Coms ( i )->variance() + elem_mult ( Coms ( i )->mean(), Coms ( i )->mean() ) );
-	}
-	// central moment
-	// pow is overloaded, we have to use another approach
-	return mom2 - elem_mult ( mean(), mean() );
-}
-
-shared_ptr<epdf> egiwmix::marginal ( const RV &rv ) const {
-	emix *tmp = new emix();
-	shared_ptr<epdf> narrow ( tmp );
-	marginal ( rv, *tmp );
-	return narrow;
-}
-
-void egiwmix::marginal ( const RV &rv, emix &target ) const {
-	bdm_assert_debug ( isnamed(), "rvs are not assigned" );
-
-	Array<shared_ptr<epdf> > Cn ( Coms.length() );
-	for ( int i = 0; i < Coms.length(); i++ ) {
-		Cn ( i ) = Coms ( i )->marginal ( rv );
-	}
-
-	target._w() = w;
-	target._Coms() = Cn;
-	target.validate();
-}
-
-egiw* 	egiwmix::approx() {
-	// NB: dimx == 1 !!!
-	// The following code might look a bit spaghetti-like,
-	// consult Dedecius, K. et al.: Partial forgetting in AR models.
-
-	double sumVecCommon;			      	// common part for many terms in eq.
-	int len = w.length();				// no. of mix components
-	int dimLS = Coms ( 1 )->_V()._D().length() - 1; 	// dim of LS
-	vec vecNu ( len );					// vector of dfms of components
-	vec vecD ( len );					// vector of LS reminders of comps.
-	vec vecCommon ( len );				// vector of common parts
-	mat matVecsTheta;				// matrix which rows are theta vects.
-
-	// fill in the vectors vecNu, vecD and matVecsTheta
-	for ( int i = 0; i < len; i++ ) {
-		vecNu.shift_left ( Coms ( i )->_nu() );
-		vecD.shift_left ( Coms ( i )->_V()._D() ( 0 ) );
-		matVecsTheta.append_row ( Coms ( i )->est_theta() );
-	}
-
-	// calculate the common parts and their sum
-	vecCommon = elem_mult ( w, elem_div ( vecNu, vecD ) );
-	sumVecCommon = sum ( vecCommon );
-
-	// LS estimator of theta
-	vec aprEstTheta ( dimLS );
-	aprEstTheta.zeros();
-	for ( int i = 0; i < len; i++ ) {
-		aprEstTheta +=  matVecsTheta.get_row ( i ) * vecCommon ( i );
-	}
-	aprEstTheta /= sumVecCommon;
-
-
-	// LS estimator of dfm
-	double aprNu;
-	double A = log ( sumVecCommon );		// Term 'A' in equation
-
-	for ( int i = 0; i < len; i++ ) {
-		A += w ( i ) * ( log ( vecD ( i ) ) - psi ( 0.5 * vecNu ( i ) ) );
-	}
-
-	aprNu = ( 1 + sqrt ( 1 + 2 * ( A - LOG2 ) / 3 ) ) / ( 2 * ( A - LOG2 ) );
-
-
-	// LS reminder (term D(0,0) in C-syntax)
-	double aprD = aprNu / sumVecCommon;
-
-	// Aproximation of cov
-	// the following code is very numerically sensitive, thus
-	// we have to eliminate decompositions etc. as much as possible
-	mat aprC = zeros ( dimLS, dimLS );
-	for ( int i = 0; i < len; i++ ) {
-		aprC += Coms ( i )->est_theta_cov().to_mat() * w ( i );
-		vec tmp = ( matVecsTheta.get_row ( i ) - aprEstTheta );
-		aprC += vecCommon ( i ) * outer_product ( tmp, tmp );
-	}
-
-	// Construct GiW pdf :: BEGIN
-	ldmat aprCinv ( inv ( aprC ) );
-	vec D = concat ( aprD, aprCinv._D() );
-	mat L = eye ( dimLS + 1 );
-	L.set_submatrix ( 1, 0, aprCinv._L() * aprEstTheta );
-	L.set_submatrix ( 1, 1, aprCinv._L() );
-	ldmat aprLD ( L, D );
-
-	egiw* aprgiw = new egiw ( 1, aprLD, aprNu );
-	return aprgiw;
-};
+void emix::from_setting ( const Setting &set ) {
+	UI::get ( Coms, set, "pdfs", UI::compulsory );
+	UI::get ( w, set, "weights", UI::compulsory );
+}
+
+void 	emix::validate (){
+	emix_base::validate();
+	dim = Coms ( 0 )->dimension();
+}
+
 
 double mprod::evallogcond ( const vec &val, const vec &cond ) {
@@ -378,8 +232,8 @@
 }
 
-vec eprod::mean() const {
+vec eprod_base::mean() const {
 	vec tmp ( dim );
-	for ( int i = 0; i < epdfs.length(); i++ ) {
-		vec pom = epdfs ( i )->mean();
+	for ( int i = 0; i < no_factors(); i++ ) {
+		vec pom = factor( i )->mean();
 		dls ( i )->pushup ( tmp, pom );
 	}
@@ -387,57 +241,29 @@
 }
 
-vec eprod::variance() const {
+vec eprod_base::variance() const {
 	vec tmp ( dim ); //second moment
-	for ( int i = 0; i < epdfs.length(); i++ ) {
-		vec pom = epdfs ( i )->mean();
+	for ( int i = 0; i < no_factors(); i++ ) {
+		vec pom = factor ( i )->mean();
 		dls ( i )->pushup ( tmp, pow ( pom, 2 ) );
 	}
 	return tmp - pow ( mean(), 2 );
 }
-vec eprod::sample() const {
+vec eprod_base::sample() const {
 	vec tmp ( dim );
-	for ( int i = 0; i < epdfs.length(); i++ ) {
-		vec pom = epdfs ( i )->sample();
+	for ( int i = 0; i < no_factors(); i++ ) {
+		vec pom = factor ( i )->sample();
 		dls ( i )->pushup ( tmp, pom );
 	}
 	return tmp;
 }
-double eprod::evallog ( const vec &val ) const {
+double eprod_base::evallog ( const vec &val ) const {
 	double tmp = 0;
-	for ( int i = 0; i < epdfs.length(); i++ ) {
-		tmp += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );
-	}
-	bdm_assert_debug ( std::isfinite ( tmp ), "Infinite" );
-	return tmp;
-}
-
-}
-// mprod::mprod ( Array<pdf*> mFacs, bool overlap) : pdf ( RV(), RV() ), n ( mFacs.length() ), epdfs ( n ), pdfs ( mFacs ), rvinds ( n ), rvcinrv ( n ), irvcs_rvc ( n ) {
-// 		int i;
-// 		bool rvaddok;
-// 		// Create rv
-// 		for ( i = 0;i < n;i++ ) {
-// 			rvaddok=rv.add ( pdfs ( i )->_rv() ); //add rv to common rvs.
-// 			// If rvaddok==false, pdfs overlap => assert error.
-// 			epdfs ( i ) = & ( pdfs ( i )->posterior() ); // add pointer to epdf
-// 		};
-// 		// Create rvc
-// 		for ( i = 0;i < n;i++ ) {
-// 			rvc.add ( pdfs ( i )->_rvc().subt ( rv ) ); //add rv to common rvs.
-// 		};
-//
-// //		independent = true;
-// 		//test rvc of pdfs and fill rvinds
-// 		for ( i = 0;i < n;i++ ) {
-// 			// find ith rv in common rv
-// 			rvsinrv ( i ) = pdfs ( i )->_rv().dataind ( rv );
-// 			// find ith rvc in common rv
-// 			rvcinrv ( i ) = pdfs ( i )->_rvc().dataind ( rv );
-// 			// find ith rvc in common rv
-// 			irvcs_rvc ( i ) = pdfs ( i )->_rvc().dataind ( rvc );
-// 			//
-// /*			if ( rvcinrv ( i ).length() >0 ) {independent = false;}
-// 			if ( irvcs_rvc ( i ).length() >0 ) {independent = false;}*/
-// 		}
-// 	};
-
+	for ( int i = 0; i < no_factors(); i++ ) {
+		tmp += factor ( i )->evallog ( dls ( i )->pushdown ( val ) );
+	}
+	//bdm_assert_debug ( std::isfinite ( tmp ), "Infinite" );
+	return tmp;
+}
+
+}
+
Index: library/bdm/stat/emix.h
===================================================================
--- library/bdm/stat/emix.h (revision 878)
+++ library/bdm/stat/emix.h (revision 886)
@@ -97,10 +97,59 @@
 
 
-
-
 private:
 	// not implemented
 	mratio ( const mratio & );
 	mratio &operator= ( const mratio & );
+};
+
+class emix; //forward
+
+/*! Base class (Interface) for mixtures 
+*/
+class emix_base : public epdf {
+	protected:
+	//! reference to vector of weights
+	vec &w;
+	//! function returning ith component
+	virtual const epdf * component(const int &i) const=0;
+	
+	virtual int no_coms() const=0;
+	
+	public:
+		
+		emix_base(vec &w0): w(w0){}
+		
+		void validate ();
+		
+		vec sample() const;
+		
+		vec mean() const;
+		
+		vec variance() const;
+		
+		double evallog ( const vec &val ) const;
+		
+		vec evallog_mat ( const mat &Val ) const;
+		
+		//! Auxiliary function that returns pdflog for each component
+		mat evallog_coms ( const mat &Val ) const;
+		
+		shared_ptr<epdf> marginal ( const RV &rv ) const;
+		//! Update already existing marginal density  \c target
+		void marginal ( const RV &rv, emix &target ) const;
+		shared_ptr<pdf> condition ( const RV &rv ) const;
+		
+		//Access methods	
+		//! returns a reference to the internal weights. Use with Care!
+		vec& _w() {
+			return w;
+		}
+		
+		const vec& _w() const {
+			return w;
+		}
+		//!access
+		const epdf* _com(int i) const {return component(i);}
+		
 };
 
@@ -115,8 +164,8 @@
 
 */
-class emix : public epdf {
+class emix : public emix_base {
 protected:
 	//! weights of the components
-	vec w;
+	vec weights;
 
 	//! Component (epdfs)
@@ -125,53 +174,11 @@
 public:
 	//! Default constructor
-	emix ( ) : epdf ( ) { }
-
-	virtual void validate ();
-
-	vec sample() const;
-
-	vec mean() const;
-
-	vec variance() const;
-
-	double evallog ( const vec &val ) const;
-
-	vec evallog_mat ( const mat &Val ) const;
-
-	//! Auxiliary function that returns pdflog for each component
-	mat evallog_coms ( const mat &Val ) const;
-
-	shared_ptr<epdf> marginal ( const RV &rv ) const;
-	//! Update already existing marginal density  \c target
-	void marginal ( const RV &rv, emix &target ) const;
-	shared_ptr<pdf> condition ( const RV &rv ) const;
-
-	//Access methods	
-	//! returns a reference to the internal weights. Use with Care!
-	vec& _w() {
-		return w;
-	}
-
-	/*!
-	  \brief returns a reference to the internal array of components. Use with Care! Set components \c Coms
-
-	  Shared pointers in Coms are kept inside this instance and
-	  shouldn't be modified after being passed to this method.
-	*/
-	Array<shared_ptr<epdf> >& _Coms ( ) {
-		return Coms;
-	}
-
-	//! returns a reference to the internal components specified by index i. Use with Care!
-	shared_ptr<epdf> _Coms ( int i ) {
-		return Coms ( i );
-	}
-
-	void set_rv ( const RV &rv ) {
-		epdf::set_rv ( rv );
-		for ( int i = 0; i < Coms.length(); i++ ) {
-			Coms ( i )->set_rv ( rv );
-		}
-	}
+	emix ( ) : emix_base ( weights) { }
+	
+	const epdf* component(const int &i) const {return Coms(i).get();}
+	void validate();
+	
+
+	int no_coms() const {return Coms.length(); }
 
 	//! Load from structure with elements:
@@ -184,84 +191,19 @@
 	//!@}
 	void from_setting ( const Setting &set );
+
+	void set_rv ( const RV &rv ) {
+		epdf::set_rv ( rv );
+		for ( int i = 0; i < no_coms(); i++ ) {
+			Coms( i )->set_rv ( rv );
+		}
+	}
+	
+	Array<shared_ptr<epdf> >& _Coms ( ) {
+			return Coms;
+		}
 };
 SHAREDPTR ( emix );
 UIREGISTER ( emix );
 
-/*!
-* \brief Mixture of egiws
-
-*/
-class egiwmix : public egiw {
-protected:
-	//! weights of the components
-	vec w;
-	//! Component (epdfs)
-	Array<egiw*> Coms;
-	//!Flag if owning Coms
-	bool destroyComs;
-public:
-	//!Default constructor
-	egiwmix ( ) : egiw ( ) {};
-
-	//! Set weights \c w and components \c Coms
-	//!By default Coms are copied inside. Parameter \c copy can be set to false if Coms live externally. Use method ownComs() if Coms should be destroyed by the destructor.
-	void set_parameters ( const vec &w, const Array<egiw*> &Coms, bool copy = false );
-	
-	//!return expected value
-	void validate();
-	
-	vec mean() const;
-
-	//!return a sample from the density
-	vec sample() const;
-
-	//!return the expected variance
-	vec variance() const;
-
-	// TODO!!! Defined to follow ANSI and/or for future development
-	void mean_mat ( mat &M, mat&R ) const {};
-	double evallog_nn ( const vec &val ) const {
-		return 0;
-	};
-	double lognc () const {
-		return 0;
-	}
-
-	shared_ptr<epdf> marginal ( const RV &rv ) const;
-	//! marginal density update
-	void marginal ( const RV &rv, emix &target ) const;
-
-//Access methods
-	//! returns a pointer to the internal mean value. Use with Care!
-	vec& _w() {
-		return w;
-	}
-	virtual ~egiwmix() {
-		if ( destroyComs ) {
-			for ( int i = 0; i < Coms.length(); i++ ) {
-				delete Coms ( i );
-			}
-		}
-	}
-	//! Auxiliary function for taking ownership of the Coms()
-	void ownComs() {
-		destroyComs = true;
-	}
-
-	//!access function
-	egiw* _Coms ( int i ) {
-		return Coms ( i );
-	}
-
-	void set_rv ( const RV &rv ) {
-		egiw::set_rv ( rv );
-		for ( int i = 0; i < Coms.length(); i++ ) {
-			Coms ( i )->set_rv ( rv );
-		}
-	}
-
-	//! Approximation of a GiW mix by a single GiW pdf
-	egiw* approx();
-};
 
 /*! \brief Chain rule decomposition of epdf
@@ -329,42 +271,53 @@
 SHAREDPTR ( mprod );
 
+
 //! Product of independent epdfs. For dependent pdfs, use mprod.
-class eprod: public epdf {
+class eprod_base: public epdf {
 protected:
-	//! Components (epdfs)
-	Array<const epdf*> epdfs;
 	//! Array of indeces
 	Array<datalink*> dls;
+	//! interface for a factor
+	virtual const epdf* factor(int i) const NOT_IMPLEMENTED(NULL); 
+	//!number of factors
+	virtual const int no_factors() const =0;
 public:
 	//! Default constructor
-	eprod () : epdfs ( 0 ), dls ( 0 ) {};
+	eprod_base () :  dls (0) {};
 	//! Set internal
-	void set_parameters ( const Array<const epdf*> &epdfs0, bool named = true ) {
-		epdfs = epdfs0;//.set_length ( epdfs0.length() );
-		dls.set_length ( epdfs.length() );
+	vec mean() const;
+
+	vec variance() const;
+
+	vec sample() const;
+
+	double evallog ( const vec &val ) const;
+
+	//!Destructor
+	~eprod_base() {
+		for ( int i = 0; i < dls.length(); i++ ) {
+			delete dls ( i );
+		}
+	}
+	void validate() {
+		dls.set_length ( no_factors() );
+		
 		bool independent = true;
-		if ( named ) {
-			for ( int i = 0; i < epdfs.length(); i++ ) {
-				independent = rv.add ( epdfs ( i )->_rv() );
-				bdm_assert_debug ( independent, "eprod:: given components are not independent." );
-			}
-			dim = rv._dsize();
-			
-		} else {
-			dim = 0;
-			for ( int i = 0; i < epdfs.length(); i++ ) {
-				dim += epdfs ( i )->dimension();
-			}
-		}
+		dim = 0;
+		for ( int i = 0; i < no_factors(); i++ ) {
+			independent = rv.add ( factor ( i )->_rv() );
+			dim += factor ( i )->dimension();
+			bdm_assert_debug ( independent, "eprod:: given components are not independent." );
+		};
+		
 		//
 		int cumdim = 0;
 		int dimi = 0;
 		int i;
-		for ( i = 0; i < epdfs.length(); i++ ) {
+		for ( i = 0; i < no_factors(); i++ ) {
 			dls ( i ) = new datalink;
-			if ( named ) {
-				dls ( i )->set_connection ( epdfs ( i )->_rv() , rv );
-			} else {
-				dimi = epdfs ( i )->dimension();
+			if ( isnamed() ) { // rvs are complete
+				dls ( i )->set_connection ( factor ( i )->_rv() , rv );
+			} else { //rvs are not reliable
+				dimi = factor ( i )->dimension();
 				dls ( i )->set_connection ( dimi, dim, linspace ( cumdim, cumdim + dimi - 1 ) );
 				cumdim += dimi;
@@ -373,28 +326,28 @@
 	
 	}
-
-	
-	vec mean() const;
-
-	vec variance() const;
-
-	vec sample() const;
-
-	double evallog ( const vec &val ) const;
-
-	//!access function
-	const epdf* operator () ( int i ) const {
-		bdm_assert_debug ( i < epdfs.length(), "wrong index" );
-		return epdfs ( i );
-	}
-
-	//!Destructor
-	~eprod() {
-		for ( int i = 0; i < epdfs.length(); i++ ) {
-			delete dls ( i );
-		}
-	}
-};
-
+};
+
+class eprod: public eprod_base{
+	protected:
+		Array<shared_ptr<epdf> > factors;
+	public:
+		const epdf* factor(int i) const {return factors(i).get();}
+		const int no_factors() const {return factors.length();}
+	void set_parameters ( const Array<shared_ptr<epdf> > &epdfs0) {
+		factors = epdfs0;
+	}
+};
+
+//! similar to eprod but used only internally -- factors are external pointers
+class eprod_internal: public eprod_base{
+	protected:
+		Array<epdf* > factors;
+		const epdf* factor(int i) const {return factors(i);}
+		const int no_factors() const {return factors.length();}
+	public:
+		void set_parameters ( const Array<epdf *> &epdfs0) {
+			factors = epdfs0;
+		}
+};
 
 /*! \brief Mixture of pdfs with constant weights, all pdfs are of equal RV and RVC
