Changeset 32 for bdm/stat

Timestamp:

03/03/08 13:00:32 (16 years ago)

Author:

smidl

Message:

test KF : estimation of R in KF is not possible! Likelihood of y_t is growing when R -> 0

Location:

bdm/stat

Files:

• libBM.cpp (modified) (6 diffs)
• libBM.h (modified) (4 diffs)
• libEF.cpp (modified) (1 diff)
• libEF.h (modified) (14 diffs)

bdm/stat/libBM.cpp

@@ -7 +7 @@
 using std::cout;
 
-void RV::init( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times, ivec in_obs ) {
+void RV::init( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times ) {
         //
         int i;
@@ -16 +16 @@
                 ( len != in_names.length() ) || \
                 ( len != in_sizes.length() ) || \
-                ( len != in_times.length() ) || \
-                ( len != in_obs.length() ) ) {
+                ( len != in_times.length() )) {
                 it_error( "RV::RV inconsistent length of input vectors." );
         }
@@ -25 +24 @@
         sizes = in_sizes;
         times = in_times;
-        obs = in_obs;
         size = 0;
         for(i=0;i<len;i++){size+=sizes(i);}
 };
 
-RV::RV ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times, ivec in_obs ) {
-        init ( in_ids, in_names, in_sizes, in_times, in_obs );
+RV::RV ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times ) {
+        init ( in_ids, in_names, in_sizes, in_times );
 }
 
-RV::RV () {};
+RV::RV () : size(0),len(0){};
+
+void RV::add (const RV &rv2) {
+        // TODO
+        size+=rv2.size;
+        len +=rv2.len;
+        ids=concat(ids,rv2.ids);
+        sizes=concat(sizes,rv2.sizes);
+        times=concat(times,rv2.times);
+        names=concat(names,rv2.names);
+//      return *this;
+};
 
 RV::RV ( ivec in_ids ) {
@@ -46 +55 @@
         }
 
-        init ( in_ids, A, ones_i( len ), zeros_i( len ), zeros_i( len ) );
+        init ( in_ids, A, ones_i( len ), zeros_i( len ) );
 }
 
 RV RV::subselect(ivec ind){
-        return RV(ids(ind), names(to_Arr(ind)), sizes(ind), times(ind), obs(ind));
+        return RV(ids(ind), names(to_Arr(ind)), sizes(ind), times(ind));
 }
 
@@ -56 +65 @@
 
 RV RV::operator()(ivec ind){
-        return RV(ids(ind), names(to_Arr(ind)), sizes(ind), times(ind), obs(ind));
+        return RV(ids(ind), names(to_Arr(ind)), sizes(ind), times(ind));
 }
 
@@ -77 +86 @@
         return indlist;
 }
+
+RV concat(const RV &rv1, const RV &rv2 ){
+        RV pom = rv1;
+        pom.add(rv2); 
+        return pom;
+}

bdm/stat/libBM.h

@@ -24 +24 @@
 * More?...
 */
+
 class RV {
+protected:
         //! size = sum of sizes
         int size;
@@ -36 +38 @@
 
 private:
-        void init ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times, ivec in_obs );
+        void init ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times);
 public:
         //! Full constructor which is called by the others
-        RV ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times, ivec in_obs );
+        RV ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times);
         //! default constructor
         RV ( ivec ids );
         //! Empty constructor will be set later
         RV ();
-        
+
         //! Printing output e.g. for debugging.
         friend std::ostream &operator<< ( std::ostream &os, const RV &rv );
@@ -50 +52 @@
         //! Return length (number of scalars) of the RV.
         int count() const {return size;} ;
+
         //TODO why not inline and later??
-        
+
         //! Find indexes of another rv in self
-        ivec find(RV rv2);
+        ivec find ( RV rv2 );
         //! Add (concat) another variable to the current one
-        RV add(RV rv2);
+        void add (const RV &rv2 );
+        //! Add (concat) another variable to the current one
+        friend RV concat (const RV &rv1, const RV &rv2 );
         //! Subtract  another variable from the current one
-        RV subt(RV rv2);
+        RV subt ( RV rv2 );
         //! Select only variables at indeces ind
-        RV subselect(ivec ind);
+        RV subselect ( ivec ind );
         //! Select only variables at indeces ind
-        RV operator()(ivec ind);
+        RV operator() ( ivec ind );
         //! Generate new \c RV with \c time shifted by delta.
-        void t(int delta);
-        //! generate a list of indeces, i.e. which 
+        void t ( int delta );
+        //! generate a list of indeces, i.e. which
         ivec indexlist();
 };
 
 
-
-
 //! Class representing function $f(x)$ of variable $x$ represented by \c rv
+
 class fnc {
 protected:
         int dimy;
-public: 
+public:
         //! function evaluates numerical value of $f(x)$ at $x=cond$
-        virtual vec eval(const vec &cond)
-        {
-                return vec(0);
+        virtual vec eval ( const vec &cond ) {
+                return vec ( 0 );
         }; //Fixme: virtual?
 
         //! access function
-        int _dimy()const{return dimy;}
-
-                //! Destructor for future use;
-                virtual ~fnc(){};
+        int _dimy() const{return dimy;}
+
+        //! Destructor for future use;
+        virtual ~fnc() {};
 };
 
 
 //! Probability density function with numerical statistics, e.g. posterior density.
+
 class epdf {
+protected:
         RV rv;
 public:
         //!default constructor
-        epdf():rv(ivec(0)){};
+        epdf() :rv ( ivec ( 0 ) ) {};
+
+        //!default constructor
+        epdf ( const RV &rv0 ) :rv ( rv0 ) {};
+
         //! Returns the required moment of the epdf
 //      virtual vec moment ( const int order = 1 );
         //! Returns a sample from the density, \f$x \sim epdf(rv)\f$
-        virtual vec sample ()=0;
+        virtual vec sample () const =0;
         //! Compute probability of argument \c val
-        virtual double eval(const vec &val){return 0.0;};
+        virtual double eval ( const vec &val ) const {return exp(this->evalpdflog(val));};
+
         //! Compute log-probability of argument \c val
-        virtual double evalpdflog(const vec &val){return 0.0;};
-
-                //! Destructor for future use;
-                virtual ~epdf(){};
-};
+        virtual double evalpdflog ( const vec &val ) const =0;
+
+        //! return expected value 
+        virtual vec mean() const =0;
+        
+        //! Destructor for future use;
+        virtual ~epdf() {};
+        //! access function
+        RV _rv() const {return rv;}
+};
+
 
 //! Conditional probability density, e.g. modeling some dependencies.
+//TODO Samplecond can be generalized
+
 class mpdf {
+protected:
         //! modeled random variable
         RV rv;
         //! random variable in condition
         RV rvc;
+        //! pointer to internal epdf
+        epdf* ep;
 public:
 
         //! Returns the required moment of the epdf
 //      virtual fnc moment ( const int order = 1 );
-        //! Returns a sample from the density conditioned on \c cond, $x \sim epdf(rv|cond)$
-        virtual vec samplecond (vec &cond, double lik){return vec(0);};
-        virtual void condition (vec &cond){};
-
-                //! Destructor for future use;
-                virtual ~mpdf(){};
+        //! Returns a sample from the density conditioned on \c cond, $x \sim epdf(rv|cond)$. \param ll is a return value of log-likelihood of the sample.
+        virtual vec samplecond ( vec &cond, double &ll ) {this->condition(cond);vec temp= ep->sample();ll=ep->evalpdflog(temp);return temp;};
+
+        virtual void condition ( const vec &cond ) {};
+        
+        virtual double evalcond (const vec &dt, const vec &cond ) {this->condition(cond);return ep->eval(dt);};
+
+        //! Destructor for future use;
+        virtual ~mpdf() {};
+
+        //! Default constructor
+        mpdf ( const RV &rv0, const RV &rvc0 ) :rv ( rv0 ),rvc ( rvc0 ) {};
+        //! access function
+        RV _rvc(){return rvc;}
+        //!access function
+        epdf& _epdf(){return *ep;}
 };
 
@@ -133 +164 @@
 
 */
+
 class DS {
 protected:
         //!Observed variables, returned by \c getdata().
-        RV Drv; 
+        RV Drv;
         //!Action variables, accepted by \c write().
         RV Urv; //
 public:
-        //! Returns full vector of observed data 
-        void getdata(vec &dt);
+        //! Returns full vector of observed data
+        void getdata ( vec &dt );
         //! Returns data records at indeces.
-        void getdata(vec &dt, ivec &indeces);
+        void getdata ( vec &dt, ivec &indeces );
         //! Accepts action variable and schedule it for application.
-        void write(vec &ut);
-        //! Accepts action variables at specific indeces        
-        void write(vec &ut, ivec &indeces);
-        /*! \brief Method that assigns random variables to the datasource. 
+        void write ( vec &ut );
+        //! Accepts action variables at specific indeces
+        void write ( vec &ut, ivec &indeces );
+        /*! \brief Method that assigns random variables to the datasource.
         Typically, the datasource will be constructed without knowledge of random variables. This method will associate existing variables with RVs.
-        
+
         (Inherited from m3k, may be deprecated soon).
         */
-        void linkrvs(RV &drv, RV &urv);
-        
-        //! Moves from $t$ to $t+1$, i.e. perfroms the actions and reads response of the system. 
+        void linkrvs ( RV &drv, RV &urv );
+
+        //! Moves from $t$ to $t+1$, i.e. perfroms the actions and reads response of the system.
         void step();
-                
+
 };
 
 /*! \brief Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities.
-        
-*/
+
+*/
+
 class BM {
-public:
+protected:
+        //!Random variable of the posterior
+        RV rv;
         //!Logarithm of marginalized data likelihood.
-        double ll;
+ double ll;
         //!  If true, the filter will compute likelihood of the data record and store it in \c ll . Set to false if you want to save time.
         bool evalll;
+public:
 
         //!Default constructor
-        BM():ll(0),evalll(true){};
-        
+        BM(const RV &rv0) :rv(rv0), ll ( 0 ),evalll ( true ) {//Fixme: test rv 
+        };
+
         /*! \brief Incremental Bayes rule
         @param dt vector of input data
         */
-        virtual void bayes ( const vec &dt) = 0;
+        virtual void bayes ( const vec &dt ) = 0;
         //! Batch Bayes rule (columns of Dt are observations)
         void bayes ( mat Dt );
         //! Returns a pointer to the epdf representing posterior density on parameters. Use with care!
-        virtual epdf* _epdf()=0;
-
-                //! Destructor for future use;
-                virtual ~BM(){};
+        virtual epdf& _epdf()=0;
+
+        //! Destructor for future use;
+        virtual ~BM() {};
+        //!access function
+        const RV& _rv() const {return rv;}
+        //!access function
+        double _ll() const {return ll;}
+};
+
+/*!
+\brief Conditional Bayesian Filter
+
+Evaluates conditional filtering density $f(rv|rvc,data)$ 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 {
+protected:
+        RV rvc;
+public:
+        //! Substitute \c val for \c rvc.
+        virtual void condition ( const vec &val ) =0;
+        BMcond(RV &rv0):rvc(rv0){};
+        virtual ~BMcond(){};
+        const RV& _rvc() const {return rvc;}
 };
 

bdm/stat/libEF.cpp

@@ -1 +1 @@
 #include <itpp/itbase.h>
+#include <itpp/base/bessel.h>
 #include "libEF.h"
+#include <math.h>
 
 using namespace itpp;
 
+Uniform_RNG UniRNG;
+Normal_RNG NorRNG;
+Gamma_RNG GamRNG;
+
 using std::cout;
 
+vec egamma::sample() const {
+        vec smp ( rv.count() );
+        int i;
+
+        for ( i=0; i<rv.count(); i++ ) {
+                GamRNG.setup ( alpha ( i ),beta ( i ) );
+                smp ( i ) = GamRNG();
+        }
+
+        return smp;
+}
+
+mat egamma::sample ( int N ) const {
+        mat Smp ( rv.count(),N );
+        int i,j;
+
+        for ( i=0; i<rv.count(); i++ ) {
+                GamRNG.setup ( alpha ( i ),beta ( i ) );
+
+                for ( j=0; j<N; j++ ) {
+                        Smp ( i,j ) = GamRNG();
+                }
+        }
+
+        return Smp;
+}
+
+double egamma::evalpdflog ( const vec &val ) const {
+        double res = 0.0; //will be added
+        int i;
+
+        for ( i=0; i<rv.count(); i++ ) {
+                res += ( alpha ( i ) - 1 ) *std::log ( val ( i ) ) + alpha ( i ) *std::log ( beta ( i ) ) - beta ( i ) *val ( i ) - lgamma ( alpha ( i ) );
+        }
+
+        return res;
+}
+
+//MGamma
+mgamma::mgamma ( const RV &rv,const RV &rvc ) : mEF ( rv,rvc ), epdf ( rv ) {vec* tmp; epdf._param ( tmp,_beta );};
+
+void mgamma::set_parameters ( double k0 ) {
+        k=k0;
+        ep = &epdf;
+        epdf.set_parameters ( k*ones ( rv.count() ),*_beta );
+};
+
+vec mgamma::samplecond ( vec &cond, double &ll ) {
+        *_beta=k / cond ;
+        vec smp = epdf.sample();
+        ll = epdf.evalpdflog ( smp );
+        return smp;
+};
+
+//Fixme repetition of mlnorm.samplecond.
+mat mgamma::samplecond ( vec &cond, vec &lik, int n ) {
+        int i;
+        int dim = rv.count();
+        mat Smp ( dim,n );
+        vec smp ( dim );
+        this->condition ( cond );
+
+        for ( i=0; i<n; i++ ) {
+                smp = epdf.sample();
+                lik ( i ) = epdf.eval ( smp );
+                Smp.set_col ( i ,smp );
+        }
+
+        return Smp;
+};
+
+ivec eEmp::resample ( RESAMPLING_METHOD method ) {
+        ivec ind=zeros_i ( n );
+        ivec N_babies = zeros_i ( n );
+        vec cumDist = cumsum ( w );
+        vec u ( n );
+        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:
+                it_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 ) ++;
+        }
+
+        // We have assigned new babies for each Particle
+        // Now, we fill the resulting index such that:
+        // * 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++ ) {
+                if ( N_babies ( i ) > 0 ) {
+                        ind ( i ) = i;
+                        N_babies ( i ) --; //this index was now replicated;
+                } else {
+                        // test if the parent has been fully replicated
+                        // 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;
+        }
+
+        return ind;
+}
+
+void eEmp::set_parameters ( const vec &w0, epdf* epdf0 ) {
+//it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
+        w=w0;
+        w/=sum ( w0 );//renormalize
+        n=w.length();
+        samples.set_size ( n );
+
+        for ( int i=0;i<n;i++ ) {samples ( i ) =epdf0->sample();}
+}

bdm/stat/libEF.h

@@ -17 +17 @@
 #include "../math/libDC.h"
 #include "libBM.h"
+#include "../itpp_ext.h"
 //#include <std>
 
 using namespace itpp;
 
+
+//! Global Uniform_RNG
+extern Uniform_RNG UniRNG;
+extern Normal_RNG NorRNG;
+extern Gamma_RNG GamRNG;
+
+
 /*!
 * \brief General conjugate exponential family posterior density.
@@ -30 +38 @@
 
 public:
+//      eEF() :epdf() {};
         //! default constructor
-        eEF() :epdf() {};
+        eEF ( const RV &rv ) :epdf ( rv ) {};
 
         virtual void tupdate ( double phi, mat &vbar, double nubar ) {};
@@ -41 +50 @@
 
 public:
-
+        mEF ( const RV &rv0, const RV &rvc0 ) :mpdf ( rv0,rvc0 ) {};
 };
 
@@ -61 +70 @@
         //! indicator if \c _iR is chached
         bool cached;
-        //! Random number generator for sampling; Fixme: what about *static* correct ?
-        Normal_RNG RNG;
         //! dimension (redundant from rv.count() for easier coding )
         int dim;
 public:
-        enorm() :eEF() {};
+//      enorm() :eEF() {};
 
         enorm ( RV &rv );
@@ -74 +81 @@
         void dupdate ( mat &v,double nu=1.0 );
 
-        vec sample();
-        mat sample ( int N );
-        double eval ( const vec &val );
-        double evalpdflog ( const vec &val );
+        vec sample() const;
+        mat sample ( int N ) const;
+        double eval ( const vec &val ) const ;
+        double evalpdflog ( const vec &val ) const;
+        vec mean()const {return mu;}
 
 //Access methods
@@ -94 +102 @@
 
 /*!
- \brief
-*/
-template<class sq_T>
-
+ \brief Gamma posterior density
+
+ Multvariate Gamma density as product of independent univariate densities.
+ \f[
+ f(x|a,b) = \prod f(x_i|a_i,b_i)
+ \f]
+*/
+
+class egamma : public eEF {
+protected:
+        vec alpha;
+        vec beta;
+public :
+        //! Default constructor
+        egamma ( const RV &rv ) :eEF ( rv ) {};
+        //! Sets parameters
+        void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;};
+        vec sample() const;
+        mat sample ( int N ) const;
+        double evalpdflog ( const vec &val ) const;
+        //! Returns poiter to alpha and beta. Potentially dengerous: use with care!
+        void _param ( vec* &a, vec* &b ) {a=&alpha;b=&beta;};
+        vec mean()const {vec pom(alpha); pom/=beta; return pom;}
+};
+
+//! Weighted mixture of epdfs with external owned components.
+class emix : public epdf {
+protected:
+        int n;
+        vec &w;
+        Array<epdf*> Coms;
+public:
+//! Default constructor
+        emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
+        void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
+        vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
+        vec sample() {it_error ( "Not implemented" );return 0;}
+};
+
+//!  Uniform distributed density on a rectangular support
+
+class euni: public epdf {
+protected:
+//! lower bound on support
+        vec low;
+//! upper bound on support
+        vec high;
+//! internal
+        vec distance;
+//! normalizing coefficients
+        double nk,lnk;
+public:
+        euni ( const RV rv ) :epdf ( rv ) {}
+        double eval ( const vec &val ) const  {return nk;}
+        double evalpdflog ( const vec &val ) const  {return lnk;}
+        vec sample() const {
+                vec smp ( rv.count() ); UniRNG.sample_vector ( rv.count(),smp );
+                return low+distance*smp;
+        }
+        void set_parameters ( const vec &low0, const vec &high0 ) {
+                distance = high0-low0;
+                it_assert_debug ( min ( distance ) >0.0,"bad support" );
+                low = low0;
+                high = high0;
+                nk = prod ( 1.0/distance );
+                lnk = log ( nk );
+        }
+        vec mean() const {vec pom=high; pom-=low; pom/=2.0; return pom;}
+};
+
+
+/*!
+ \brief Normal distributed linear function with linear function of mean value;
+
+ Mean value $mu=A*rvc$.
+*/
+template<class sq_T>
 class mlnorm : public mEF {
         enorm<sq_T> epdf;
+        vec* _mu; //cached epdf.mu;
         mat A;
 public:
         //! Constructor
-        mlnorm ( RV &rv,RV &rvc, mat &A, sq_T &R );
+        mlnorm ( RV &rv,RV &rvc );
+        void set_parameters ( const  mat &A, const sq_T &R );
         //!Generate one sample of the posterior
         vec samplecond ( vec &cond, double &lik );
+        //!Generate matrix of samples of the posterior
         mat samplecond ( vec &cond, vec &lik, int n );
         void condition ( vec &cond );
 };
 
+/*!
+\brief  Gamma random walk
+
+Mean value, $\mu$, of this density is given by \c rvc .
+Standard deviation of the random walk is proportional to one $k$-th the mean.
+This is achieved by setting $\alpha=k$ and $\beta=k/\mu$.
+
+The standard deviation of the walk is then: $\mu/\sqrt(k)$.
+*/
+class mgamma : public mEF {
+        egamma epdf;
+        double k;
+        vec* _beta;
+
+public:
+        //! Constructor
+        mgamma ( const RV &rv,const RV &rvc );
+        void set_parameters ( double k );
+        //!Generate one sample of the posterior
+        vec samplecond ( vec &cond, double &lik );
+        //!Generate matrix of samples of the posterior
+        mat samplecond ( vec &cond, vec &lik, int n );
+        void condition ( const vec &val ) {*_beta=k/val;};
+};
+
+//! Switch between various resampling methods.
+enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
+/*!
+\brief Weighted empirical density
+
+Used e.g. in particle filters.
+*/
+class eEmp: public epdf {
+protected :
+        //! Number of particles
+        int n;
+        vec w;
+        Array<vec> samples;
+public:
+        eEmp ( const RV &rv0 ,int n0) :epdf ( rv0 ),n(n0),w(n),samples(n) {};
+        void set_parameters ( const vec &w0, epdf* pdf0 );
+        //! Potentially dangerous, use with care.
+        vec& _w()  {return w;};
+        Array<vec>& _samples() {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 );
+        vec sample() const {it_error ( "Not implemented" );return 0;}
+        double evalpdflog(const vec &val) const {it_error ( "Not implemented" );return 0.0;}
+        vec mean()const {vec pom=zeros(rv.count()); 
+                for (int i=0;i<n;i++){pom+=samples(i)*w(i);}
+                return pom;
+        }
+};
+
+
 ////////////////////////
 
 template<class sq_T>
-enorm<sq_T>::enorm ( RV &rv ) :eEF(), mu ( rv.count() ),R(rv.count()),_iR(rv.count()),cached ( false ),dim ( rv.count() ) {};
+enorm<sq_T>::enorm ( RV &rv ) :eEF(rv), mu ( rv.count() ),R ( rv.count() ),_iR ( rv.count() ),cached ( false ),dim ( rv.count() ) {};
 
 template<class sq_T>
@@ -120 +259 @@
         mu = mu0;
         R = R0;
-        if(_iR.rows()!=R.rows()) _iR=R; // memory allocation!
+        if ( _iR.rows() !=R.rows() ) _iR=R; // memory allocation!
         R.inv ( _iR ); //update cache
         cached=true;
@@ -136 +275 @@
 
 template<class sq_T>
-vec enorm<sq_T>::sample() {
+vec enorm<sq_T>::sample() const {
         vec x ( dim );
-        RNG.sample_vector ( dim,x );
+        NorRNG.sample_vector ( dim,x );
         vec smp = R.sqrt_mult ( x );
 
@@ -146 +285 @@
 
 template<class sq_T>
-mat enorm<sq_T>::sample ( int N ) {
+mat enorm<sq_T>::sample ( int N )const {
         mat X ( dim,N );
         vec x ( dim );
@@ -153 +292 @@
 
         for ( i=0;i<N;i++ ) {
-                RNG.sample_vector ( dim,x );
+                NorRNG.sample_vector ( dim,x );
                 pom = R.sqrt_mult ( x );
                 pom +=mu;
@@ -163 +302 @@
 
 template<class sq_T>
-double enorm<sq_T>::eval ( const vec &val ) {
-        return exp ( evalpdflog ( val ) );
-};
-
-template<class sq_T>
-double enorm<sq_T>::evalpdflog ( const vec &val ) {
-        if ( !cached ) {R.inv ( _iR );cached=true;}
-
-        return -0.5* ( R.cols() *0.79817986835811504957 +R.logdet() +_iR.qform ( mu-val ) );
-};
-
-
-template<class sq_T>
-mlnorm<sq_T>::mlnorm ( RV &rv,RV &rvc, mat &A, sq_T &R ) {
-        int dim = rv.count();
-        vec mu ( dim );
-
-        epdf = enorm<sq_T> ( rv,mu,R );
+double enorm<sq_T>::eval ( const vec &val ) const {
+        double pdfl,e;
+        pdfl = evalpdflog ( val );
+        e = exp ( pdfl );
+        return e;
+};
+
+template<class sq_T>
+double enorm<sq_T>::evalpdflog ( const vec &val ) const {
+        if ( !cached ) {it_error("this should not happen, see cached");}
+
+        // 1.83787706640935 = log(2pi)
+        return -0.5* ( R.cols() * 1.83787706640935 +R.logdet() +_iR.qform ( mu-val ) );
+};
+
+
+template<class sq_T>
+mlnorm<sq_T>::mlnorm ( RV &rv0,RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv ),A ( rv0.count(),rv0.count() ) {
+        _mu = epdf._mu();
+}
+
+template<class sq_T>
+void mlnorm<sq_T>::set_parameters ( const mat &A0, const sq_T &R0 ) {
+        epdf.set_parameters ( zeros ( rv.count() ),R0 );
+        A = A0;
 }
 
@@ -199 +345 @@
         this->condition ( cond );
 
-        for ( i=0; i<dim; i++ ) {
+        for ( i=0; i<n; i++ ) {
                 smp = epdf.sample();
                 lik ( i ) = epdf.eval ( smp );
@@ -210 +356 @@
 template<class sq_T>
 void mlnorm<sq_T>::condition ( vec &cond ) {
-        epdf.mu = A*cond;
+        *_mu = A*cond;
 //R is already assigned;
 }
@@ -216 +362 @@
 ///////////
 
-//#ifdef HAVE_EXTERN_TEMPLATE
-
-//  extern template class enorm<fsqmat>;
-//  extern template class enorm<ldmat>;
-//#endif
 
 #endif //EF_H