Changeset 294

Timestamp:

03/12/09 22:18:43 (15 years ago)

Author:

smidl

Message:

tr2245

Files:

• CMake_modules/FindITPP.cmake (modified) (1 diff)
• bdm/estim/ekf_templ.h (modified) (2 diffs)
• bdm/math/chmat.h (modified) (1 diff)
• bdm/stat/libBM.h (modified) (1 diff)
• bdm/stat/libEF.h (modified) (4 diffs)
• pmsm/TR2245/CMakeLists.txt (modified) (1 diff)
• pmsm/TR2245/unitsteps.cfg (modified) (2 diffs)
• pmsm/TR2245/unitsteps.cpp (modified) (1 diff)
• tests/chmat_test.cpp (modified) (2 diffs)
• tests/testSmp.cpp (modified) (1 diff)

CMake_modules/FindITPP.cmake

@@ -23 +23 @@
   /usr/lib64
   /usr/lib
-  NO_DEFAULT_PATH
+  #NO_DEFAULT_PATH
 )
 

bdm/estim/ekf_templ.h

@@ -27 +27 @@
 };
 
-//!Extended Kalman filter in Choleski form with unknown \c Q
-class EKFCh_unQ : public EKFCh {
+//!Extended Kalman filter in Choleski form with unknown diagonal \c Q
+class EKFCh_dQ : public EKFCh {
 public:
         void condition ( const vec &Q0 ) {
@@ -35 +35 @@
                 preA.set_submatrix ( dimy+dimx,dimy,Q._Ch() );
         };
+};
+
+//!Extended Kalman filter in Choleski form with unknown \c Q
+class EKFCh_chQ : public EKFCh {
+        public:
+                void condition ( const vec &chQ0 ) {
+                        Q.setCh ( chQ0);
+                //from EKF
+                        preA.set_submatrix ( dimy+dimx,dimy,Q._Ch() );
+                };
 };
 

bdm/math/chmat.h

@@ -66 +66 @@
         //! Access function
         mat & _Ch() {return Ch;}
+        //! Access function
+        const mat & _Ch() const {return Ch;}
         //! Access functions
         void setD ( const vec &nD ) {Ch=diag ( sqrt(nD) );}
+        //! Access functions
+        void setCh ( const vec &chQ ) {
+                it_assert_debug(chQ.length()==dim*dim,""); 
+                copy_vector(dim*dim, chQ._data(), Ch._data()); 
+        }
         //! Access functions
         void setD ( const vec &nD, int i ) {for ( int j=i;j<nD.length();j++ ) {Ch( j,j ) =sqrt(nD ( j-i ));}} //Fixme can be more general

bdm/stat/libBM.h

@@ -319 +319 @@
                 return temp;
         };
-        //! Returns \param N samples from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv \param ll is a return value of log-likelihood of the sample.
+        //! Returns \param N samples from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv 
         virtual mat samplecond_m ( const vec &cond, int N ) {
                 this->condition ( cond );

bdm/stat/libEF.h

@@ -19 +19 @@
 //#include <std>
 
-namespace bdm {
+namespace bdm
+{
 
 
 //! Global Uniform_RNG
-extern Uniform_RNG UniRNG;
+        extern Uniform_RNG UniRNG;
 //! Global Normal_RNG
-extern Normal_RNG NorRNG;
+        extern Normal_RNG NorRNG;
 //! Global Gamma_RNG
-extern Gamma_RNG GamRNG;
-
-/*!
-* \brief General conjugate exponential family posterior density.
-
-* More?...
-*/
-
-class eEF : public epdf {
-public:
+        extern Gamma_RNG GamRNG;
+
+        /*!
+        * \brief General conjugate exponential family posterior density.
+
+        * More?...
+        */
+
+        class eEF : public epdf
+        {
+                public:
 //      eEF() :epdf() {};
-        //! default constructor
-        eEF ( ) :epdf ( ) {};
-        //! logarithm of the normalizing constant, \f$\mathcal{I}\f$
-        virtual double lognc() const =0;
-        //!TODO decide if it is really needed
-        virtual void dupdate ( mat &v ) {it_error ( "Not implemented" );};
-        //!Evaluate normalized log-probability
-        virtual double evallog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
-        //!Evaluate normalized log-probability
-        virtual double evallog ( const vec &val ) const {double tmp;tmp= evallog_nn ( val )-lognc();it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); return tmp;}
-        //!Evaluate normalized log-probability for many samples
-        virtual vec evallog ( const mat &Val ) const {
-                vec x ( Val.cols() );
-                for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog_nn ( Val.get_col ( i ) ) ;}
-                return x-lognc();
-        }
-        //!Power of the density, used e.g. to flatten the density
-        virtual void pow ( double p ) {it_error ( "Not implemented" );};
-};
-
-/*!
-* \brief Exponential family model.
-
-* More?...
-*/
-
-class mEF : public mpdf {
-
-public:
+                        //! default constructor
+                        eEF ( ) :epdf ( ) {};
+                        //! logarithm of the normalizing constant, \f$\mathcal{I}\f$
+                        virtual double lognc() const =0;
+                        //!TODO decide if it is really needed
+                        virtual void dupdate ( mat &v ) {it_error ( "Not implemented" );};
+                        //!Evaluate normalized log-probability
+                        virtual double evallog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
+                        //!Evaluate normalized log-probability
+                        virtual double evallog ( const vec &val ) const {double tmp;tmp= evallog_nn ( val )-lognc();it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); return tmp;}
+                        //!Evaluate normalized log-probability for many samples
+                        virtual vec evallog ( const mat &Val ) const
+                        {
+                                vec x ( Val.cols() );
+                                for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog_nn ( Val.get_col ( i ) ) ;}
+                                return x-lognc();
+                        }
+                        //!Power of the density, used e.g. to flatten the density
+                        virtual void pow ( double p ) {it_error ( "Not implemented" );};
+        };
+
+        /*!
+        * \brief Exponential family model.
+
+        * More?...
+        */
+
+        class mEF : public mpdf
+        {
+
+                public:
+                        //! Default constructor
+                        mEF ( ) :mpdf ( ) {};
+        };
+
+//! Estimator for Exponential family
+        class BMEF : public BM
+        {
+                protected:
+                        //! forgetting factor
+                        double frg;
+                        //! cached value of lognc() in the previous step (used in evaluation of \c ll )
+                        double last_lognc;
+                public:
+                        //! Default constructor (=empty constructor)
+                        BMEF ( double frg0=1.0 ) :BM ( ), frg ( frg0 ) {}
+                        //! Copy constructor
+                        BMEF ( const BMEF &B ) :BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
+                        //!get statistics from another model
+                        virtual void set_statistics ( const BMEF* BM0 ) {it_error ( "Not implemented" );};
+                        //! Weighted update of sufficient statistics (Bayes rule)
+                        virtual void bayes ( const vec &data, const double w ) {};
+                        //original Bayes
+                        void bayes ( const vec &dt );
+                        //!Flatten the posterior according to the given BMEF (of the same type!)
+                        virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
+                        //!Flatten the posterior as if to keep nu0 data
+//      virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
+
+                        BMEF* _copy_ () const {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
+        };
+
+        template<class sq_T>
+        class mlnorm;
+
+        /*!
+        * \brief Gaussian density with positive definite (decomposed) covariance matrix.
+
+        * More?...
+        */
+        template<class sq_T>
+        class enorm : public eEF
+        {
+                protected:
+                        //! mean value
+                        vec mu;
+                        //! Covariance matrix in decomposed form
+                        sq_T R;
+                public:
+                        //!\name Constructors
+                        //!@{
+
+                        enorm ( ) :eEF ( ), mu ( ),R ( ) {};
+                        enorm ( const vec &mu,const sq_T &R ) {set_parameters ( mu,R );}
+                        void set_parameters ( const vec &mu,const sq_T &R );
+                        //!@}
+
+                        //! \name Mathematical operations
+                        //!@{
+
+                        //! dupdate in exponential form (not really handy)
+                        void dupdate ( mat &v,double nu=1.0 );
+
+                        vec sample() const;
+                        mat sample ( int N ) const;
+                        double evallog_nn ( const vec &val ) const;
+                        double lognc () const;
+                        vec mean() const {return mu;}
+                        vec variance() const {return diag ( R.to_mat() );}
+//      mlnorm<sq_T>* condition ( const RV &rvn ) const ;
+                        mpdf* condition ( const RV &rvn ) const ;
+//      enorm<sq_T>* marginal ( const RV &rv ) const;
+                        epdf* marginal ( const RV &rv ) const;
+                        //!@}
+
+                        //! \name Access to attributes
+                        //!@{
+
+                        vec& _mu() {return mu;}
+                        void set_mu ( const vec mu0 ) { mu=mu0;}
+                        sq_T& _R() {return R;}
+                        const sq_T& _R() const {return R;}
+                        //!@}
+
+        };
+
+        /*!
+        * \brief Gauss-inverse-Wishart density stored in LD form
+
+        * For \f$p\f$-variate densities, given rv.count() should be \f$p\times\f$ V.rows().
+        *
+        */
+        class egiw : public eEF
+        {
+                protected:
+                        //! Extended information matrix of sufficient statistics
+                        ldmat V;
+                        //! Number of data records (degrees of freedom) of sufficient statistics
+                        double nu;
+                        //! Dimension of the output
+                        int dimx;
+                        //! Dimension of the regressor
+                        int nPsi;
+                public:
+                        //!\name Constructors
+                        //!@{
+                        egiw() :eEF() {};
+                        egiw ( int dimx0, ldmat V0, double nu0=-1.0 ) :eEF() {set_parameters ( dimx0,V0, nu0 );};
+
+                        void set_parameters ( int dimx0, ldmat V0, double nu0=-1.0 )
+                        {
+                                dimx=dimx0;
+                                nPsi = V0.rows()-dimx;
+                                dim = dimx* ( dimx+nPsi ); // size(R) + size(Theta)
+
+                                V=V0;
+                                if ( nu0<0 )
+                                {
+                                        nu = 0.1 +nPsi +2*dimx +2; // +2 assures finite expected value of R
+                                        // terms before that are sufficient for finite normalization
+                                }
+                                else
+                                {
+                                        nu=nu0;
+                                }
+                        }
+                        //!@}
+
+                        vec sample() const;
+                        vec mean() const;
+                        vec variance() const;
+                        void mean_mat ( mat &M, mat&R ) const;
+                        //! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
+                        double evallog_nn ( const vec &val ) const;
+                        double lognc () const;
+                        void pow ( double p ) {V*=p;nu*=p;};
+
+                        //! \name Access attributes
+                        //!@{
+
+                        ldmat& _V() {return V;}
+                        const ldmat& _V() const {return V;}
+                        double& _nu()  {return nu;}
+                        const double& _nu() const {return nu;}
+                        //!@}
+        };
+
+        /*! \brief Dirichlet posterior density
+
+        Continuous Dirichlet density of \f$n\f$-dimensional variable \f$x\f$
+        \f[
+        f(x|\beta) = \frac{\Gamma[\gamma]}{\prod_{i=1}^{n}\Gamma(\beta_i)} \prod_{i=1}^{n}x_i^{\beta_i-1}
+        \f]
+        where \f$\gamma=\sum_i \beta_i\f$.
+        */
+        class eDirich: public eEF
+        {
+                protected:
+                        //!sufficient statistics
+                        vec beta;
+                        //!speedup variable
+                        double gamma;
+                public:
+                        //!\name Constructors
+                        //!@{
+
+                        eDirich () : eEF ( ) {};
+                        eDirich ( const eDirich &D0 ) : eEF () {set_parameters ( D0.beta );};
+                        eDirich ( const vec &beta0 ) {set_parameters ( beta0 );};
+                        void set_parameters ( const vec &beta0 )
+                        {
+                                beta= beta0;
+                                dim = beta.length();
+                                gamma = sum ( beta );
+                        }
+                        //!@}
+
+                        vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
+                        vec mean() const {return beta/gamma;};
+                        vec variance() const {return elem_mult ( beta, ( beta+1 ) ) / ( gamma* ( gamma+1 ) );}
+                        //! In this instance, val is ...
+                        double evallog_nn ( const vec &val ) const
+                        {
+                                double tmp; tmp= ( beta-1 ) *log ( val );               it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
+                                return tmp;
+                        };
+                        double lognc () const
+                        {
+                                double tmp;
+                                double gam=sum ( beta );
+                                double lgb=0.0;
+                                for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
+                                tmp= lgb-lgamma ( gam );
+                                it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
+                                return tmp;
+                        };
+                        //!access function
+                        vec& _beta()  {return beta;}
+                        //!Set internal parameters
+        };
+
+//! \brief Estimator for Multinomial density
+        class multiBM : public BMEF
+        {
+                protected:
+                        //! Conjugate prior and posterior
+                        eDirich est;
+                        //! Pointer inside est to sufficient statistics
+                        vec &beta;
+                public:
+                        //!Default constructor
+                        multiBM ( ) : BMEF ( ),est ( ),beta ( est._beta() )
+                        {
+                                if ( beta.length() >0 ) {last_lognc=est.lognc();}
+                                else{last_lognc=0.0;}
+                        }
+                        //!Copy constructor
+                        multiBM ( const multiBM &B ) : BMEF ( B ),est ( B.est ),beta ( est._beta() ) {}
+                        //! Sets sufficient statistics to match that of givefrom mB0
+                        void set_statistics ( const BM* mB0 ) {const multiBM* mB=dynamic_cast<const multiBM*> ( mB0 ); beta=mB->beta;}
+                        void bayes ( const vec &dt )
+                        {
+                                if ( frg<1.0 ) {beta*=frg;last_lognc=est.lognc();}
+                                beta+=dt;
+                                if ( evalll ) {ll=est.lognc()-last_lognc;}
+                        }
+                        double logpred ( const vec &dt ) const
+                        {
+                                eDirich pred ( est );
+                                vec &beta = pred._beta();
+
+                                double lll;
+                                if ( frg<1.0 )
+                                        {beta*=frg;lll=pred.lognc();}
+                                else
+                                        if ( evalll ) {lll=last_lognc;}
+                                        else{lll=pred.lognc();}
+
+                                beta+=dt;
+                                return pred.lognc()-lll;
+                        }
+                        void flatten ( const BMEF* B )
+                        {
+                                const multiBM* E=dynamic_cast<const multiBM*> ( B );
+                                // sum(beta) should be equal to sum(B.beta)
+                                const vec &Eb=E->beta;//const_cast<multiBM*> ( E )->_beta();
+                                beta*= ( sum ( Eb ) /sum ( beta ) );
+                                if ( evalll ) {last_lognc=est.lognc();}
+                        }
+                        const epdf& posterior() const {return est;};
+                        const eDirich* _e() const {return &est;};
+                        void set_parameters ( const vec &beta0 )
+                        {
+                                est.set_parameters ( beta0 );
+                                if ( evalll ) {last_lognc=est.lognc();}
+                        }
+        };
+
+        /*!
+         \brief Gamma posterior density
+
+         Multivariate Gamma density as product of independent univariate densities.
+         \f[
+         f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i)
+         \f]
+        */
+
+        class egamma : public eEF
+        {
+                protected:
+                        //! Vector \f$\alpha\f$
+                        vec alpha;
+                        //! Vector \f$\beta\f$
+                        vec beta;
+                public :
+                        //! \name Constructors
+                        //!@{
+                        egamma ( ) :eEF ( ), alpha ( 0 ), beta ( 0 ) {};
+                        egamma ( const vec &a, const vec &b ) {set_parameters ( a, b );};
+                        void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;dim = alpha.length();};
+                        //!@}
+
+                        vec sample() const;
+                        //! TODO: is it used anywhere?
+//      mat sample ( int N ) const;
+                        double evallog ( const vec &val ) const;
+                        double lognc () const;
+                        //! Returns poiter to alpha and beta. Potentially dengerous: use with care!
+                        vec& _alpha() {return alpha;}
+                        vec& _beta() {return beta;}
+                        vec mean() const {return elem_div ( alpha,beta );}
+                        vec variance() const {return elem_div ( alpha,elem_mult ( beta,beta ) ); }
+        };
+
+        /*!
+         \brief Inverse-Gamma posterior density
+
+         Multivariate inverse-Gamma density as product of independent univariate densities.
+         \f[
+         f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i)
+         \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)
+        \f]
+        This relation is used in sampling.
+         */
+
+        class eigamma : public egamma
+        {
+                protected:
+                public :
+                        //! \name Constructors
+                        //! 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 );}
+        };
+        /*
+        //! Weighted mixture of epdfs with external owned components.
+        class emix : public epdf {
+        protected:
+                int n;
+                vec &w;
+                Array<epdf*> Coms;
+        public:
         //! Default constructor
-        mEF ( ) :mpdf ( ) {};
-};
-
-//! Estimator for Exponential family
-class BMEF : public BM {
-protected:
-        //! forgetting factor
-        double frg;
-        //! cached value of lognc() in the previous step (used in evaluation of \c ll )
-        double last_lognc;
-public:
-        //! Default constructor (=empty constructor)
-        BMEF ( double frg0=1.0 ) :BM ( ), frg ( frg0 ) {}
-        //! Copy constructor
-        BMEF ( const BMEF &B ) :BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
-        //!get statistics from another model
-        virtual void set_statistics ( const BMEF* BM0 ) {it_error ( "Not implemented" );};
-        //! Weighted update of sufficient statistics (Bayes rule)
-        virtual void bayes ( const vec &data, const double w ) {};
-        //original Bayes
-        void bayes ( const vec &dt );
-        //!Flatten the posterior according to the given BMEF (of the same type!)
-        virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
-        //!Flatten the posterior as if to keep nu0 data
-//      virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
-
-        BMEF* _copy_ () const {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
-};
-
-template<class sq_T>
-class mlnorm;
-
-/*!
-* \brief Gaussian density with positive definite (decomposed) covariance matrix.
-
-* More?...
-*/
-template<class sq_T>
-class enorm : public eEF {
-protected:
-        //! mean value
-        vec mu;
-        //! Covariance matrix in decomposed form
-        sq_T R;
-public:
-        //!\name Constructors
-        //!@{
-
-        enorm ( ) :eEF ( ), mu ( ),R ( ) {};
-        enorm ( const vec &mu,const sq_T &R ) {set_parameters ( mu,R );}
-        void set_parameters ( const vec &mu,const sq_T &R );
-        //!@}
-
-        //! \name Mathematical operations
-        //!@{
-
-        //! dupdate in exponential form (not really handy)
-        void dupdate ( mat &v,double nu=1.0 );
-
-        vec sample() const;
-        mat sample ( int N ) const;
-        double evallog_nn ( const vec &val ) const;
-        double lognc () const;
-        vec mean() const {return mu;}
-        vec variance() const {return diag ( R.to_mat() );}
-//      mlnorm<sq_T>* condition ( const RV &rvn ) const ;
-        mpdf* condition ( const RV &rvn ) const ;
-//      enorm<sq_T>* marginal ( const RV &rv ) const;
-        epdf* marginal ( const RV &rv ) const;
-        //!@}
-
-        //! \name Access to attributes
-        //!@{
-
-        vec& _mu() {return mu;}
-        void set_mu ( const vec mu0 ) { mu=mu0;}
-        sq_T& _R() {return R;}
-        const sq_T& _R() const {return R;}
-        //!@}
-
-};
-
-/*!
-* \brief Gauss-inverse-Wishart density stored in LD form
-
-* For \f$p\f$-variate densities, given rv.count() should be \f$p\times\f$ V.rows().
-*
-*/
-class egiw : public eEF {
-protected:
-        //! Extended information matrix of sufficient statistics
-        ldmat V;
-        //! Number of data records (degrees of freedom) of sufficient statistics
-        double nu;
-        //! Dimension of the output
-        int dimx;
-        //! Dimension of the regressor
-        int nPsi;
-public:
-        //!\name Constructors
-        //!@{
-        egiw() :eEF() {};
-        egiw ( int dimx0, ldmat V0, double nu0=-1.0 ) :eEF() {set_parameters ( dimx0,V0, nu0 );};
-
-        void set_parameters ( int dimx0, ldmat V0, double nu0=-1.0 ) {
-                dimx=dimx0;
-                nPsi = V0.rows()-dimx;
-                dim = dimx* ( dimx+nPsi ); // size(R) + size(Theta)
-
-                V=V0;
-                if ( nu0<0 ) {
-                        nu = 0.1 +nPsi +2*dimx +2; // +2 assures finite expected value of R
-                        // terms before that are sufficient for finite normalization
-                }
-                else {
-                        nu=nu0;
-                }
-        }
-        //!@}
-
-        vec sample() const;
-        vec mean() const;
-        vec variance() const;
-        void mean_mat ( mat &M, mat&R ) const;
-        //! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
-        double evallog_nn ( const vec &val ) const;
-        double lognc () const;
-        void pow ( double p ) {V*=p;nu*=p;};
-
-        //! \name Access attributes
-        //!@{
-
-        ldmat& _V() {return V;}
-        const ldmat& _V() const {return V;}
-        double& _nu()  {return nu;}
-        const double& _nu() const {return nu;}
-        //!@}
-};
-
-/*! \brief Dirichlet posterior density
-
-Continuous Dirichlet density of \f$n\f$-dimensional variable \f$x\f$
-\f[
-f(x|\beta) = \frac{\Gamma[\gamma]}{\prod_{i=1}^{n}\Gamma(\beta_i)} \prod_{i=1}^{n}x_i^{\beta_i-1}
-\f]
-where \f$\gamma=\sum_i \beta_i\f$.
-*/
-class eDirich: public eEF {
-protected:
-        //!sufficient statistics
-        vec beta;
-        //!speedup variable
-        double gamma;
-public:
-        //!\name Constructors
-        //!@{
-
-        eDirich () : eEF ( ) {};
-        eDirich ( const eDirich &D0 ) : eEF () {set_parameters ( D0.beta );};
-        eDirich ( const vec &beta0 ) {set_parameters ( beta0 );};
-        void set_parameters ( const vec &beta0 ) {
-                beta= beta0;
-                dim = beta.length();
-                gamma = sum ( beta );
-        }
-        //!@}
-
-        vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
-        vec mean() const {return beta/gamma;};
-        vec variance() const {return elem_mult ( beta, ( beta+1 ) ) / ( gamma* ( gamma+1 ) );}
-        //! In this instance, val is ...
-        double evallog_nn ( const vec &val ) const {
-                double tmp; tmp= ( beta-1 ) *log ( val );               it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
-                return tmp;
-        };
-        double lognc () const {
-                double tmp;
-                double gam=sum ( beta );
-                double lgb=0.0;
-                for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
-                tmp= lgb-lgamma ( gam );
-                it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
-                return tmp;
-        };
-        //!access function
-        vec& _beta()  {return beta;}
-        //!Set internal parameters
-};
-
-//! \brief Estimator for Multinomial density
-class multiBM : public BMEF {
-protected:
-        //! Conjugate prior and posterior
-        eDirich est;
-        //! Pointer inside est to sufficient statistics
-        vec &beta;
-public:
-        //!Default constructor
-        multiBM ( ) : BMEF ( ),est ( ),beta ( est._beta() ) {if ( beta.length() >0 ) {last_lognc=est.lognc();}else{last_lognc=0.0;}}
-        //!Copy constructor
-        multiBM ( const multiBM &B ) : BMEF ( B ),est ( B.est ),beta ( est._beta() ) {}
-        //! Sets sufficient statistics to match that of givefrom mB0
-        void set_statistics ( const BM* mB0 ) {const multiBM* mB=dynamic_cast<const multiBM*> ( mB0 ); beta=mB->beta;}
-        void bayes ( const vec &dt ) {
-                if ( frg<1.0 ) {beta*=frg;last_lognc=est.lognc();}
-                beta+=dt;
-                if ( evalll ) {ll=est.lognc()-last_lognc;}
-        }
-        double logpred ( const vec &dt ) const {
-                eDirich pred ( est );
-                vec &beta = pred._beta();
-
-                double lll;
-                if ( frg<1.0 )
-                        {beta*=frg;lll=pred.lognc();}
-                else
-                        if ( evalll ) {lll=last_lognc;}
-                        else{lll=pred.lognc();}
-
-                beta+=dt;
-                return pred.lognc()-lll;
-        }
-        void flatten ( const BMEF* B ) {
-                const multiBM* E=dynamic_cast<const multiBM*> ( B );
-                // sum(beta) should be equal to sum(B.beta)
-                const vec &Eb=E->beta;//const_cast<multiBM*> ( E )->_beta();
-                beta*= ( sum ( Eb ) /sum ( beta ) );
-                if ( evalll ) {last_lognc=est.lognc();}
-        }
-        const epdf& posterior() const {return est;};
-        const eDirich* _e() const {return &est;};
-        void set_parameters ( const vec &beta0 ) {
-                est.set_parameters ( beta0 );
-                if ( evalll ) {last_lognc=est.lognc();}
-        }
-};
-
-/*!
- \brief Gamma posterior density
-
- Multivariate Gamma density as product of independent univariate densities.
- \f[
- f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i)
- \f]
-*/
-
-class egamma : public eEF {
-protected:
-        //! Vector \f$\alpha\f$
-        vec alpha;
-        //! Vector \f$\beta\f$
-        vec beta;
-public :
-        //! \name Constructors
-        //!@{
-        egamma ( ) :eEF ( ), alpha(0), beta(0) {};
-        egamma ( const vec &a, const vec &b ) {set_parameters ( a, b );};
-        void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;dim = alpha.length();};
-        //!@}
-
-        vec sample() const;
-        //! TODO: is it used anywhere?
-//      mat sample ( int N ) const;
-        double evallog ( const vec &val ) const;
-        double lognc () const;
-        //! Returns poiter to alpha and beta. Potentially dengerous: use with care!
-        vec& _alpha() {return alpha;}
-        vec& _beta() {return beta;}
-        vec mean() const {return elem_div ( alpha,beta );}
-        vec variance() const {return elem_div ( alpha,elem_mult ( beta,beta ) ); }
-};
-
-/*!
- \brief Inverse-Gamma posterior density
-
- Multivariate inverse-Gamma density as product of independent univariate densities.
- \f[
- f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i)
- \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)
-\f]
-This relation is used in sampling.
- */
-
-class eigamma : public egamma {
-protected:
-public :
-        //! \name Constructors
-        //! 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 );}
-};
-/*
-//! 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;}
-};
-*/
+                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:
+        class euni: public epdf
+        {
+                protected:
 //! lower bound on support
-        vec low;
+                        vec low;
 //! upper bound on support
-        vec high;
+                        vec high;
 //! internal
-        vec distance;
+                        vec distance;
 //! normalizing coefficients
-        double nk;
+                        double nk;
 //! cache of log( \c nk )
-        double lnk;
-public:
-        //! \name Constructors
-        //!@{
-        euni ( ) :epdf ( ) {}
-        euni ( const vec &low0, const vec &high0 ) {set_parameters ( low0,high0 );}
-        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 );
-                dim = low.length();
-        }
-        //!@}
-
-        double eval ( const vec &val ) const  {return nk;}
-        double evallog ( const vec &val ) const  {return lnk;}
-        vec sample() const {
-                vec smp ( dim );
+                        double lnk;
+                public:
+                        //! \name Constructors
+                        //!@{
+                        euni ( ) :epdf ( ) {}
+                        euni ( const vec &low0, const vec &high0 ) {set_parameters ( low0,high0 );}
+                        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 );
+                                dim = low.length();
+                        }
+                        //!@}
+
+                        double eval ( const vec &val ) const  {return nk;}
+                        double evallog ( const vec &val ) const  {return lnk;}
+                        vec sample() const
+                        {
+                                vec smp ( dim );
 #pragma omp critical
-                UniRNG.sample_vector ( dim ,smp );
-                return low+elem_mult ( distance,smp );
-        }
-        //! set values of \c low and \c high
-        vec mean() const {return ( high-low ) /2.0;}
-        vec variance() const {return ( pow ( high,2 ) +pow ( low,2 ) +elem_mult ( high,low ) ) /3.0;}
-};
-
-
-/*!
- \brief Normal distributed linear function with linear function of mean value;
-
- Mean value \f$mu=A*rvc+mu_0\f$.
-*/
-template<class sq_T>
-class mlnorm : public mEF {
-protected:
-        //! Internal epdf that arise by conditioning on \c rvc
-        enorm<sq_T> epdf;
-        mat A;
-        vec mu_const;
-        vec& _mu; //cached epdf.mu;
-public:
-        //! \name Constructors
-        //!@{
-        mlnorm ( ) :mEF (),epdf ( ),A ( ),_mu ( epdf._mu() ) {ep =&epdf; };
-        mlnorm ( const  mat &A, const vec &mu0, const sq_T &R ) :epdf ( ),_mu ( epdf._mu() ) {
-                ep =&epdf; set_parameters ( A,mu0,R );
-        };
-        //! Set \c A and \c R
-        void set_parameters ( const  mat &A, const vec &mu0, const sq_T &R );
-        //!@}
-        //! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
-        void condition ( const vec &cond );
-
-        //!access function
-        vec& _mu_const() {return mu_const;}
-        //!access function
-        mat& _A() {return A;}
-        //!access function
-        mat _R() {return epdf._R().to_mat();}
-
-        template<class sq_M>
-        friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M> &ml );
-};
+                                UniRNG.sample_vector ( dim ,smp );
+                                return low+elem_mult ( distance,smp );
+                        }
+                        //! set values of \c low and \c high
+                        vec mean() const {return ( high-low ) /2.0;}
+                        vec variance() const {return ( pow ( high,2 ) +pow ( low,2 ) +elem_mult ( high,low ) ) /3.0;}
+        };
+
+
+        /*!
+         \brief Normal distributed linear function with linear function of mean value;
+
+         Mean value \f$mu=A*rvc+mu_0\f$.
+        */
+        template<class sq_T>
+        class mlnorm : public mEF
+        {
+                protected:
+                        //! Internal epdf that arise by conditioning on \c rvc
+                        enorm<sq_T> epdf;
+                        mat A;
+                        vec mu_const;
+                        vec& _mu; //cached epdf.mu;
+                public:
+                        //! \name Constructors
+                        //!@{
+                        mlnorm ( ) :mEF (),epdf ( ),A ( ),_mu ( epdf._mu() ) {ep =&epdf; };
+                        mlnorm ( const  mat &A, const vec &mu0, const sq_T &R ) :epdf ( ),_mu ( epdf._mu() )
+                        {
+                                ep =&epdf; set_parameters ( A,mu0,R );
+                        };
+                        //! Set \c A and \c R
+                        void set_parameters ( const  mat &A, const vec &mu0, const sq_T &R );
+                        //!@}
+                        //! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
+                        void condition ( const vec &cond );
+
+                        //!access function
+                        vec& _mu_const() {return mu_const;}
+                        //!access function
+                        mat& _A() {return A;}
+                        //!access function
+                        mat _R() {return epdf._R().to_mat();}
+
+                        template<class sq_M>
+                        friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M> &ml );
+        };
 
 //! Mpdf with general function for mean value
-template<class sq_T>
-class mgnorm : public mEF {
-protected:
-        //! Internal epdf that arise by conditioning on \c rvc
-        enorm<sq_T> epdf;
-        vec &mu;
-        fnc* g;
-public:
-        //!default constructor
-        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->dimension() ), R0 );}
-        void condition ( const vec &cond ) {mu=g->eval ( cond );};
-};
-
-/*! (Approximate) Student t density with linear function of mean value
-
-The internal epdf of this class is of the type of a Gaussian (enorm).
-However, each conditioning is trying to assure the best possible approximation by taking into account the zeta function. See [] for reference.
-
-Perhaps a moment-matching technique?
-*/
-class mlstudent : public mlnorm<ldmat> {
-protected:
-        ldmat Lambda;
-        ldmat &_R;
-        ldmat Re;
-public:
-        mlstudent ( ) :mlnorm<ldmat> (),
-                        Lambda (),      _R ( epdf._R() ) {}
-        void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 ) {
-                it_assert_debug ( A0.rows() ==mu0.length(),"" );
-                it_assert_debug ( R0.rows() ==A0.rows(),"" );
-
-                epdf.set_parameters ( mu0,Lambda ); //
-                A = A0;
-                mu_const = mu0;
-                Re=R0;
-                Lambda = Lambda0;
-        }
-        void condition ( const vec &cond ) {
-                _mu = A*cond + mu_const;
-                double zeta;
-                //ugly hack!
-                if ( ( cond.length() +1 ) ==Lambda.rows() ) {
-                        zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) );
-                }
-                else {
-                        zeta = Lambda.invqform ( cond );
-                }
-                _R = Re;
-                _R*= ( 1+zeta );// / ( nu ); << nu is in Re!!!!!!
-        };
-
-};
-/*!
-\brief  Gamma random walk
-
-Mean value, \f$\mu\f$, of this density is given by \c rvc .
-Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
-This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
-
-The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
-*/
-class mgamma : public mEF {
-protected:
-        //! Internal epdf that arise by conditioning on \c rvc
-        egamma epdf;
-        //! Constant \f$k\f$
-        double k;
-        //! cache of epdf.beta
-        vec &_beta;
-
-public:
-        //! Constructor
-        mgamma ( ) : mEF ( ), epdf (), _beta ( epdf._beta() ) {ep=&epdf;};
-        //! Set value of \c k
-        void set_parameters ( double k, const vec &beta0 );
-        void condition ( const vec &val ) {_beta=k/val;};
-};
-
-/*!
-\brief  Inverse-Gamma random walk
-
-Mean value, \f$ \mu \f$, of this density is given by \c rvc .
-Standard deviation of the random walk is proportional to one \f$ k \f$-th the mean.
-This is achieved by setting \f$ \alpha=\mu/k^2+2 \f$ and \f$ \beta=\mu(\alpha-1)\f$.
-
-The standard deviation of the walk is then: \f$ \mu/\sqrt(k)\f$.
- */
-class migamma : public mEF {
-protected:
-        //! Internal epdf that arise by conditioning on \c rvc
-        eigamma epdf;
-        //! Constant \f$k\f$
-        double k;
-        //! cache of epdf.alpha
-        vec &_alpha;
-        //! cache of epdf.beta
-        vec &_beta;
-
-public:
-        //! \name Constructors
-        //!@{
-        migamma ( ) : mEF (), epdf ( ), _alpha ( epdf._alpha() ), _beta ( epdf._beta() ) {ep=&epdf;};
-        migamma ( const migamma &m ) : mEF (), epdf ( m.epdf ), _alpha ( epdf._alpha() ), _beta ( epdf._beta() ) {ep=&epdf;};
-        //!@}
-
-        //! Set value of \c k
-        void set_parameters ( int len, double k0 ) {
-                k=k0;
-                epdf.set_parameters ( ( 1.0/ ( k*k ) +2.0 ) *ones ( len ) /*alpha*/, ones ( len ) /*beta*/ );
-                dimc = dimension();
-        };
-        void condition ( const vec &val ) {
-                _beta=elem_mult ( val, ( _alpha-1.0 ) );
-        };
-};
-
-/*!
-\brief  Gamma random walk around a fixed point
-
-Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
-\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
-
-Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
-This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
-
-The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
-*/
-class mgamma_fix : public mgamma {
-protected:
-        //! parameter l
-        double l;
-        //! reference vector
-        vec refl;
-public:
-        //! Constructor
-        mgamma_fix ( ) : mgamma ( ),refl () {};
-        //! Set value of \c k
-        void set_parameters ( double k0 , vec ref0, double l0 ) {
-                mgamma::set_parameters ( k0, ref0 );
-                refl=pow ( ref0,1.0-l0 );l=l0;
-                dimc=dimension();
-        };
-
-        void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); _beta=k/mean;};
-};
-
-
-/*!
-\brief  Inverse-Gamma random walk around a fixed point
-
-Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
-\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
-
-==== Check == vv =
-Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
-This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
-
-The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
- */
-class migamma_ref : public migamma {
-protected:
-        //! parameter l
-        double l;
-        //! reference vector
-        vec refl;
-public:
-        //! Constructor
-        migamma_ref ( ) : migamma (),refl ( ) {};
-        //! Set value of \c k
-        void set_parameters ( double k0 , vec ref0, double l0 ) {
-                migamma::set_parameters ( ref0.length(), k0 );
-                refl=pow ( ref0,1.0-l0 );
-                l=l0;
-                dimc = dimension();
-        };
-
-        void condition ( const vec &val ) {
-                vec mean=elem_mult ( refl,pow ( val,l ) );
-                migamma::condition ( mean );
-        };
-};
-
-/*! Log-Normal probability density
- only allow diagonal covariances!
-
-Density of the form \f$ \log(x)\sim \mathcal{N}(\mu,\sigma^2), i.e.
-\f[
-x \sim \frac{1}{x\sigma\sqrt{2\pi}}\exp{-\frac{1}{2\sigma^2}(\log(x)-\mu)}
-\f]
-
-*/
-class elognorm: public enorm<ldmat>{
-        public:
-                vec sample() const {return exp(enorm<ldmat>::sample());};
-                vec mean() const {vec var=enorm<ldmat>::variance();return exp(mu - 0.5*var);};
-                
-};
-
-/*!
-\brief  Log-Normal random walk
-
-Mean value, \f$\mu\f$, is...
-
-==== Check == vv =
-Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
-This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
-
-The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
- */
-class mlognorm : public mpdf {
-        protected:
-                elognorm eno;
-        //! parameter 1/2*sigma^2
-                double sig2;
-                //! access
-                vec &mu;
-        public:
-        //! Constructor
-                mlognorm ( ) : eno (), mu(eno._mu()) {ep=&eno;};
-        //! Set value of \c k
-                void set_parameters ( int size, double k) {
-                        sig2 = 0.5*log(k*k+1);
-                        eno.set_parameters(zeros(size),2*sig2*eye(size));
+        template<class sq_T>
+        class mgnorm : public mEF
+        {
+                protected:
+                        //! Internal epdf that arise by conditioning on \c rvc
+                        enorm<sq_T> epdf;
+                        vec &mu;
+                        fnc* g;
+                public:
+                        //!default constructor
+                        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->dimension() ), R0 );}
+                        void condition ( const vec &cond ) {mu=g->eval ( cond );};
+        };
+
+        /*! (Approximate) Student t density with linear function of mean value
+
+        The internal epdf of this class is of the type of a Gaussian (enorm).
+        However, each conditioning is trying to assure the best possible approximation by taking into account the zeta function. See [] for reference.
+
+        Perhaps a moment-matching technique?
+        */
+        class mlstudent : public mlnorm<ldmat>
+        {
+                protected:
+                        ldmat Lambda;
+                        ldmat &_R;
+                        ldmat Re;
+                public:
+                        mlstudent ( ) :mlnorm<ldmat> (),
+                                        Lambda (),      _R ( epdf._R() ) {}
+                        void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 )
+                        {
+                                it_assert_debug ( A0.rows() ==mu0.length(),"" );
+                                it_assert_debug ( R0.rows() ==A0.rows(),"" );
+
+                                epdf.set_parameters ( mu0,Lambda ); //
+                                A = A0;
+                                mu_const = mu0;
+                                Re=R0;
+                                Lambda = Lambda0;
+                        }
+                        void condition ( const vec &cond )
+                        {
+                                _mu = A*cond + mu_const;
+                                double zeta;
+                                //ugly hack!
+                                if ( ( cond.length() +1 ) ==Lambda.rows() )
+                                {
+                                        zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) );
+                                }
+                                else
+                                {
+                                        zeta = Lambda.invqform ( cond );
+                                }
+                                _R = Re;
+                                _R*= ( 1+zeta );// / ( nu ); << nu is in Re!!!!!!
+                        };
+
+        };
+        /*!
+        \brief  Gamma random walk
+
+        Mean value, \f$\mu\f$, of this density is given by \c rvc .
+        Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
+        This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
+
+        The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
+        */
+        class mgamma : public mEF
+        {
+                protected:
+                        //! Internal epdf that arise by conditioning on \c rvc
+                        egamma epdf;
+                        //! Constant \f$k\f$
+                        double k;
+                        //! cache of epdf.beta
+                        vec &_beta;
+
+                public:
+                        //! Constructor
+                        mgamma ( ) : mEF ( ), epdf (), _beta ( epdf._beta() ) {ep=&epdf;};
+                        //! Set value of \c k
+                        void set_parameters ( double k, const vec &beta0 );
+                        void condition ( const vec &val ) {_beta=k/val;};
+        };
+
+        /*!
+        \brief  Inverse-Gamma random walk
+
+        Mean value, \f$ \mu \f$, of this density is given by \c rvc .
+        Standard deviation of the random walk is proportional to one \f$ k \f$-th the mean.
+        This is achieved by setting \f$ \alpha=\mu/k^2+2 \f$ and \f$ \beta=\mu(\alpha-1)\f$.
+
+        The standard deviation of the walk is then: \f$ \mu/\sqrt(k)\f$.
+         */
+        class migamma : public mEF
+        {
+                protected:
+                        //! Internal epdf that arise by conditioning on \c rvc
+                        eigamma epdf;
+                        //! Constant \f$k\f$
+                        double k;
+                        //! cache of epdf.alpha
+                        vec &_alpha;
+                        //! cache of epdf.beta
+                        vec &_beta;
+
+                public:
+                        //! \name Constructors
+                        //!@{
+                        migamma ( ) : mEF (), epdf ( ), _alpha ( epdf._alpha() ), _beta ( epdf._beta() ) {ep=&epdf;};
+                        migamma ( const migamma &m ) : mEF (), epdf ( m.epdf ), _alpha ( epdf._alpha() ), _beta ( epdf._beta() ) {ep=&epdf;};
+                        //!@}
+
+                        //! Set value of \c k
+                        void set_parameters ( int len, double k0 )
+                        {
+                                k=k0;
+                                epdf.set_parameters ( ( 1.0/ ( k*k ) +2.0 ) *ones ( len ) /*alpha*/, ones ( len ) /*beta*/ );
+                                dimc = dimension();
+                        };
+                        void condition ( const vec &val )
+                        {
+                                _beta=elem_mult ( val, ( _alpha-1.0 ) );
+                        };
+        };
+
+        /*!
+        \brief  Gamma random walk around a fixed point
+
+        Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
+        \f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
+
+        Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
+        This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
+
+        The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
+        */
+        class mgamma_fix : public mgamma
+        {
+                protected:
+                        //! parameter l
+                        double l;
+                        //! reference vector
+                        vec refl;
+                public:
+                        //! Constructor
+                        mgamma_fix ( ) : mgamma ( ),refl () {};
+                        //! Set value of \c k
+                        void set_parameters ( double k0 , vec ref0, double l0 )
+                        {
+                                mgamma::set_parameters ( k0, ref0 );
+                                refl=pow ( ref0,1.0-l0 );l=l0;
+                                dimc=dimension();
+                        };
+
+                        void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); _beta=k/mean;};
+        };
+
+
+        /*!
+        \brief  Inverse-Gamma random walk around a fixed point
+
+        Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
+        \f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
+
+        ==== Check == vv =
+        Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
+        This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
+
+        The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
+         */
+        class migamma_ref : public migamma
+        {
+                protected:
+                        //! parameter l
+                        double l;
+                        //! reference vector
+                        vec refl;
+                public:
+                        //! Constructor
+                        migamma_ref ( ) : migamma (),refl ( ) {};
+                        //! Set value of \c k
+                        void set_parameters ( double k0 , vec ref0, double l0 )
+                        {
+                                migamma::set_parameters ( ref0.length(), k0 );
+                                refl=pow ( ref0,1.0-l0 );
+                                l=l0;
+                                dimc = dimension();
+                        };
+
+                        void condition ( const vec &val )
+                        {
+                                vec mean=elem_mult ( refl,pow ( val,l ) );
+                                migamma::condition ( mean );
+                        };
+        };
+
+        /*! Log-Normal probability density
+         only allow diagonal covariances!
+
+        Density of the form \f$ \log(x)\sim \mathcal{N}(\mu,\sigma^2), i.e.
+        \f[
+        x \sim \frac{1}{x\sigma\sqrt{2\pi}}\exp{-\frac{1}{2\sigma^2}(\log(x)-\mu)}
+        \f]
+
+        */
+        class elognorm: public enorm<ldmat>
+        {
+                public:
+                        vec sample() const {return exp ( enorm<ldmat>::sample() );};
+                        vec mean() const {vec var=enorm<ldmat>::variance();return exp ( mu - 0.5*var );};
+
+        };
+
+        /*!
+        \brief  Log-Normal random walk
+
+        Mean value, \f$\mu\f$, is...
+
+        ==== Check == vv =
+        Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
+        This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
+
+        The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
+         */
+        class mlognorm : public mpdf
+        {
+                protected:
+                        elognorm eno;
+                        //! parameter 1/2*sigma^2
+                        double sig2;
+                        //! access
+                        vec &mu;
+                public:
+                        //! Constructor
+                        mlognorm ( ) : eno (), mu ( eno._mu() ) {ep=&eno;};
+                        //! Set value of \c k
+                        void set_parameters ( int size, double k )
+                        {
+                                sig2 = 0.5*log ( k*k+1 );
+                                eno.set_parameters ( zeros ( size ),2*sig2*eye ( size ) );
+
+                                dimc = size;
+                        };
+
+                        void condition ( const vec &val )
+                        {
+                                mu=log ( val )-sig2;//elem_mult ( refl,pow ( val,l ) );
+                        };
+        };
+
+        /*! inverse Wishart density defined on Choleski decomposition
+
+        */
+        class eWishartCh : public epdf
+        {
+                protected:
+                        //! Upper-Triagle of Choleski decomposition of \f$ \Psi \f$
+                        chmat Y;
+                        //! dimension of matrix  \f$ \Psi \f$
+                        int p;
+                        //! degrees of freedom  \f$ \nu \f$
+                        double delta;
+                public:
+                        void set_parameters ( const mat &Y0, const double delta0 ) {Y=chmat ( Y0 );delta=delta0; p=Y.rows(); dim = p*p; }
+                        mat sample_mat() const
+                        {
+                                mat X=zeros ( p,p );
+
+                                //sample diagonal
+                                for ( int i=0;i<p;i++ )
+                                {
+                                        GamRNG.setup ( 0.5* ( delta-i ) , 0.5 ); // no +1 !! index if from 0
+#pragma omp critical
+                                        X ( i,i ) =sqrt ( GamRNG() );
+                                }
+                                //do the rest
+                                for ( int i=0;i<p;i++ )
+                                {
+                                        for ( int j=i+1;j<p;j++ )
+                                        {
+#pragma omp critical
+                                                X ( i,j ) =NorRNG.sample();
+                                        }
+                                }
+                                return X*Y._Ch();// return upper triangular part of the decomposition
+                        }
+                        vec sample () const
+                        {
+                                return vec ( sample_mat()._data(),p*p );
+                        }
+                        //! fast access function y0 will be copied into Y.Ch.
+                        void setY ( const mat &Ch0 ) {copy_vector ( dim,Ch0._data(), Y._Ch()._data() );}
+                        //! fast access function y0 will be copied into Y.Ch.
+                        void _setY ( const vec &ch0 ) {copy_vector ( dim, ch0._data(), Y._Ch()._data() ); }
+                        //! access function
+                        const chmat& getY()const {return Y;}
+        };
+
+        class eiWishartCh: public epdf
+        {
+                protected:
+                        eWishartCh W;
+                        int p;
+                        double delta;
+                public:
+                        void set_parameters ( const mat &Y0, const double delta0) {
+                                delta = delta0;
+                                W.set_parameters ( inv ( Y0 ),delta0 ); 
+                                dim = W.dimension(); p=Y0.rows();
+                        }
+                        vec sample() const {mat iCh; iCh=inv ( W.sample_mat() ); return vec ( iCh._data(),dim );}
+                        void _setY ( const vec &y0 )
+                        {
+                                mat Ch ( p,p );
+                                mat iCh ( p,p );
+                                copy_vector ( dim, y0._data(), Ch._data() );
+                                
+                                iCh=inv ( Ch );
+                                W.setY ( iCh );
+                        }
+                        virtual double evallog ( const vec &val ) const {
+                                chmat X(p);
+                                const chmat& Y=W.getY();
+                                 
+                                copy_vector(p*p,val._data(),X._Ch()._data());
+                                chmat iX(p);X.inv(iX);
+                                // compute  
+//                              \frac{ |\Psi|^{m/2}|X|^{-(m+p+1)/2}e^{-tr(\Psi X^{-1})/2} }{ 2^{mp/2}\Gamma_p(m/2)},
+                                mat M=Y.to_mat()*iX.to_mat();
+                                
+                                double log1 = 0.5*p*(2*Y.logdet())-0.5*(delta+p+1)*(2*X.logdet())-0.5*trace(M); 
+                                //Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!
+                                
+/*                              if (0) {
+                                        mat XX=X.to_mat();
+                                        mat YY=Y.to_mat();
+                                        
+                                        double log2 = 0.5*p*log(det(YY))-0.5*(delta+p+1)*log(det(XX))-0.5*trace(YY*inv(XX)); 
+                                        cout << log1 << "," << log2 << endl;
+                                }*/
+                                return log1;                            
+                        };
                         
-                        dimc = size;
-                };
-
-                void condition ( const vec &val ) {
-                        mu=log(val)-sig2;//elem_mult ( refl,pow ( val,l ) );
-                };
-};
-
-/*! inverse Wishart density
-
-*/
-class iW : public epdf{
-        public:
-                vec sample(){return vec();}
-};
+        };
+
+        class rwiWishartCh : public mpdf
+        {
+                protected:
+                        eiWishartCh eiW;
+                        //!square root of \f$ \nu-p-1 \f$ - needed for computation of \f$ \Psi \f$ from conditions
+                        double sqd;
+                        //reference point for diagonal
+                        vec refl;
+                        double l;
+                        int p;
+                public:
+                        void set_parameters ( int p0, double k, vec ref0, double l0  )
+                        {
+                                p=p0;
+                                double delta = 2/(k*k)+p+3;
+                                sqd=sqrt ( delta-p-1 );
+                                l=l0;
+                                refl=pow(ref0,1-l);
+                                
+                                eiW.set_parameters ( eye ( p ),delta );
+                                ep=&eiW;
+                                dimc=eiW.dimension();
+                        }
+                        void condition ( const vec &c ) {
+                                vec z=c;
+                                int ri=0;
+                                for(int i=0;i<p*p;i+=(p+1)){//trace diagonal element
+                                        z(i) = pow(z(i),l)*refl(ri);
+                                        ri++;
+                                }
+
+                                eiW._setY ( sqd*z );
+                        }
+        };
 
 //! 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;
-        //! Sample weights \f$w\f$
-        vec w;
-        //! Samples \f$x^{(i)}, i=1..n\f$
-        Array<vec> samples;
-public:
-        //! \name Constructors
-        //!@{
-        eEmp ( ) :epdf ( ),w ( ),samples ( ) {};
-        eEmp ( const eEmp &e ) : epdf ( e ), w ( e.w ), samples ( e.samples ) {};
-        //!@}
-
-        //! Set samples and weights
-        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_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;};
-        //! Potentially dangerous, use with care.
-        const vec& _w() const {return w;};
-        //! access function
-        Array<vec>& _samples() {return samples;};
-        //! access function
-        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 );
-        //! inherited operation : NOT implemneted
-        vec sample() const {it_error ( "Not implemented" );return 0;}
-        //! inherited operation : NOT implemneted
-        double evallog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
-        vec mean() const {
-                vec pom=zeros ( dim );
-                for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
-                return pom;
-        }
-        vec variance() const {
-                vec pom=zeros ( dim );
-                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 );}
-                        }
-                }
-        }
-};
+        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;
+                        //! Sample weights \f$w\f$
+                        vec w;
+                        //! Samples \f$x^{(i)}, i=1..n\f$
+                        Array<vec> samples;
+                public:
+                        //! \name Constructors
+                        //!@{
+                        eEmp ( ) :epdf ( ),w ( ),samples ( ) {};
+                        eEmp ( const eEmp &e ) : epdf ( e ), w ( e.w ), samples ( e.samples ) {};
+                        //!@}
+
+                        //! Set samples and weights
+                        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_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;};
+                        //! Potentially dangerous, use with care.
+                        const vec& _w() const {return w;};
+                        //! access function
+                        Array<vec>& _samples() {return samples;};
+                        //! access function
+                        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 );
+                        //! inherited operation : NOT implemneted
+                        vec sample() const {it_error ( "Not implemented" );return 0;}
+                        //! inherited operation : NOT implemneted
+                        double evallog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
+                        vec mean() const
+                        {
+                                vec pom=zeros ( dim );
+                                for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
+                                return pom;
+                        }
+                        vec variance() const
+                        {
+                                vec pom=zeros ( dim );
+                                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 );}
+                                        }
+                                }
+                        }
+        };
 
 
 ////////////////////////
 
-template<class sq_T>
-void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) {
+        template<class sq_T>
+        void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 )
+        {
 //Fixme test dimensions of mu0 and R0;
-        mu = mu0;
-        R = R0;
-        dim = mu0.length();
-};
-
-template<class sq_T>
-void enorm<sq_T>::dupdate ( mat &v, double nu ) {
-        //
-};
+                mu = mu0;
+                R = R0;
+                dim = mu0.length();
+        };
+
+        template<class sq_T>
+        void enorm<sq_T>::dupdate ( mat &v, double nu )
+        {
+                //
+        };
 
 // template<class sq_T>
@@ -802 +981 @@
 // };
 
-template<class sq_T>
-vec enorm<sq_T>::sample() const {
-        vec x ( dim );
-#pragma omp critical
-        NorRNG.sample_vector ( dim,x );
-        vec smp = R.sqrt_mult ( x );
-
-        smp += mu;
-        return smp;
-};
-
-template<class sq_T>
-mat enorm<sq_T>::sample ( int N ) const {
-        mat X ( dim,N );
-        vec x ( dim );
-        vec pom;
-        int i;
-
-        for ( i=0;i<N;i++ ) {
+        template<class sq_T>
+        vec enorm<sq_T>::sample() const
+        {
+                vec x ( dim );
 #pragma omp critical
                 NorRNG.sample_vector ( dim,x );
-                pom = R.sqrt_mult ( x );
-                pom +=mu;
-                X.set_col ( i, pom );
-        }
-
-        return X;
-};
+                vec smp = R.sqrt_mult ( x );
+
+                smp += mu;
+                return smp;
+        };
+
+        template<class sq_T>
+        mat enorm<sq_T>::sample ( int N ) const
+        {
+                mat X ( dim,N );
+                vec x ( dim );
+                vec pom;
+                int i;
+
+                for ( i=0;i<N;i++ )
+                {
+#pragma omp critical
+                        NorRNG.sample_vector ( dim,x );
+                        pom = R.sqrt_mult ( x );
+                        pom +=mu;
+                        X.set_col ( i, pom );
+                }
+
+                return X;
+        };
 
 // template<class sq_T>
@@ -839 +1021 @@
 // };
 
-template<class sq_T>
-double enorm<sq_T>::evallog_nn ( const vec &val ) const {
-        // 1.83787706640935 = log(2pi)
-        double tmp=-0.5* ( R.invqform ( mu-val ) );// - lognc();
-        return  tmp;
-};
-
-template<class sq_T>
-inline double enorm<sq_T>::lognc () const {
-        // 1.83787706640935 = log(2pi)
-        double tmp=0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
-        return tmp;
-};
-
-template<class sq_T>
-void mlnorm<sq_T>::set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 ) {
-        it_assert_debug ( A0.rows() ==mu0.length(),"" );
-        it_assert_debug ( A0.rows() ==R0.rows(),"" );
-
-        epdf.set_parameters ( zeros ( A0.rows() ),R0 );
-        A = A0;
-        mu_const = mu0;
-        dimc=A0.cols();
-}
+        template<class sq_T>
+        double enorm<sq_T>::evallog_nn ( const vec &val ) const
+        {
+                // 1.83787706640935 = log(2pi)
+                double tmp=-0.5* ( R.invqform ( mu-val ) );// - lognc();
+                return  tmp;
+        };
+
+        template<class sq_T>
+        inline double enorm<sq_T>::lognc () const
+        {
+                // 1.83787706640935 = log(2pi)
+                double tmp=0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
+                return tmp;
+        };
+
+        template<class sq_T>
+        void mlnorm<sq_T>::set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 )
+        {
+                it_assert_debug ( A0.rows() ==mu0.length(),"" );
+                it_assert_debug ( A0.rows() ==R0.rows(),"" );
+
+                epdf.set_parameters ( zeros ( A0.rows() ),R0 );
+                A = A0;
+                mu_const = mu0;
+                dimc=A0.cols();
+        }
 
 // template<class sq_T>
@@ -889 +1074 @@
 // }
 
-template<class sq_T>
-void mlnorm<sq_T>::condition ( const vec &cond ) {
-        _mu = A*cond + mu_const;
+        template<class sq_T>
+        void mlnorm<sq_T>::condition ( const vec &cond )
+        {
+                _mu = A*cond + mu_const;
 //R is already assigned;
-}
-
-template<class sq_T>
-epdf* enorm<sq_T>::marginal ( const RV &rvn ) const {
-        it_assert_debug ( isnamed(), "rv description is not assigned" );
-        ivec irvn = rvn.dataind ( rv );
-
-        sq_T Rn ( R,irvn ); //select rows and columns of R
-
-        enorm<sq_T>* tmp = new enorm<sq_T>;
-        tmp->set_rv ( rvn );
-        tmp->set_parameters ( mu ( irvn ), Rn );
-        return tmp;
-}
-
-template<class sq_T>
-mpdf* enorm<sq_T>::condition ( const RV &rvn ) const {
-
-        it_assert_debug ( isnamed(),"rvs are not assigned" );
-
-        RV rvc = rv.subt ( rvn );
-        it_assert_debug ( ( rvc._dsize() +rvn._dsize() ==rv._dsize() ),"wrong rvn" );
-        //Permutation vector of the new R
-        ivec irvn = rvn.dataind ( rv );
-        ivec irvc = rvc.dataind ( rv );
-        ivec perm=concat ( irvn , irvc );
-        sq_T Rn ( R,perm );
-
-        //fixme - could this be done in general for all sq_T?
-        mat S=Rn.to_mat();
-        //fixme
-        int n=rvn._dsize()-1;
-        int end=R.rows()-1;
-        mat S11 = S.get ( 0,n, 0, n );
-        mat S12 = S.get ( 0, n , rvn._dsize(), end );
-        mat S22 = S.get ( rvn._dsize(), end, rvn._dsize(), end );
-
-        vec mu1 = mu ( irvn );
-        vec mu2 = mu ( irvc );
-        mat A=S12*inv ( S22 );
-        sq_T R_n ( S11 - A *S12.T() );
-
-        mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( );
-        tmp->set_rv ( rvn ); tmp->set_rvc ( rvc );
-        tmp->set_parameters ( A,mu1-A*mu2,R_n );
-        return tmp;
-}
+        }
+
+        template<class sq_T>
+        epdf* enorm<sq_T>::marginal ( const RV &rvn ) const
+        {
+                it_assert_debug ( isnamed(), "rv description is not assigned" );
+                ivec irvn = rvn.dataind ( rv );
+
+                sq_T Rn ( R,irvn ); //select rows and columns of R
+
+                enorm<sq_T>* tmp = new enorm<sq_T>;
+                tmp->set_rv ( rvn );
+                tmp->set_parameters ( mu ( irvn ), Rn );
+                return tmp;
+        }
+
+        template<class sq_T>
+        mpdf* enorm<sq_T>::condition ( const RV &rvn ) const
+        {
+
+                it_assert_debug ( isnamed(),"rvs are not assigned" );
+
+                RV rvc = rv.subt ( rvn );
+                it_assert_debug ( ( rvc._dsize() +rvn._dsize() ==rv._dsize() ),"wrong rvn" );
+                //Permutation vector of the new R
+                ivec irvn = rvn.dataind ( rv );
+                ivec irvc = rvc.dataind ( rv );
+                ivec perm=concat ( irvn , irvc );
+                sq_T Rn ( R,perm );
+
+                //fixme - could this be done in general for all sq_T?
+                mat S=Rn.to_mat();
+                //fixme
+                int n=rvn._dsize()-1;
+                int end=R.rows()-1;
+                mat S11 = S.get ( 0,n, 0, n );
+                mat S12 = S.get ( 0, n , rvn._dsize(), end );
+                mat S22 = S.get ( rvn._dsize(), end, rvn._dsize(), end );
+
+                vec mu1 = mu ( irvn );
+                vec mu2 = mu ( irvc );
+                mat A=S12*inv ( S22 );
+                sq_T R_n ( S11 - A *S12.T() );
+
+                mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( );
+                tmp->set_rv ( rvn ); tmp->set_rvc ( rvc );
+                tmp->set_parameters ( A,mu1-A*mu2,R_n );
+                return tmp;
+        }
 
 ///////////
 
-template<class sq_T>
-std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_T> &ml ) {
-        os << "A:"<< ml.A<<endl;
-        os << "mu:"<< ml.mu_const<<endl;
-        os << "R:" << ml.epdf._R().to_mat() <<endl;
-        return os;
-};
+        template<class sq_T>
+        std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_T> &ml )
+        {
+                os << "A:"<< ml.A<<endl;
+                os << "mu:"<< ml.mu_const<<endl;
+                os << "R:" << ml.epdf._R().to_mat() <<endl;
+                return os;
+        };
 
 }

pmsm/TR2245/CMakeLists.txt

@@ -1 +1 @@
 # Make sure the compiler can find include files from our Bdm library.
 EXEC (unitsteps pmsmsim)
+EXEC (wishart pmsmsim)
+EXEC (pmsm_wishart pmsmsim)

pmsm/TR2245/unitsteps.cfg

@@ -1 +1 @@
 ndat = 9000;
-Npart = 200;
+Npart = 1000;
 
 //simulator 
@@ -13 +13 @@
 };
 // For EKF
-dQ=[1e-6, 1e-6, 0.001, 0.0001];
-dR=[1e-6, 1e-6];
+dQ=[1e-5, 1e-5, 1e-4, 1e-5];
+dR=[1e-8, 1e-8];
 

pmsm/TR2245/unitsteps.cpp

@@ -68 +68 @@
 
         RV rQ ( "{Q }","4" );
-        EKFCh_unQ KFEp ;
+        EKFCh_dQ KFEp ;
         KFEp.set_parameters ( &fxu,&hxu,Q,R );
         KFEp.set_est ( mu0, chmat ( zeros ( 4 ) ) );
 
-        MPF<EKFCh_unQ> M;
+        MPF<EKFCh_dQ> M;
         M.set_parameters ( evolQ,evolQ,Npart );
         // initialize

tests/chmat_test.cpp

@@ -17 +17 @@
    cout << "Testing constructors:" << endl
        << "A = " << A << endl
+       << "Ch = " << Ch._Ch() << endl
        << "Ch.to_mat() = " << Ch.to_mat() << endl << endl;
        
@@ -49 +50 @@
    cout << "A+vv' = " << A+outer_product(v,v) << endl
                         <<      "opupdt(Ch,v) = " << Ch2.to_mat() << endl <<endl;
-                                
+
+   vec vCh=vec(Ch._Ch()._data(),9);
+   cout << "vectorized Ch: " << vCh <<endl;
+   chmat nCh(3);
+   nCh.setCh(vCh);
+   cout << "Ch from vec: " << nCh._Ch() <<endl;
+
 }

tests/testSmp.cpp

@@ -116 +116 @@
         disp(concat(eN.mean(),eG.mean()), epV,Smp);
         
+        cout << "======= eWishart ======== " << endl;
+        mat wM="1.0 0.9; 0.9 1.0";
+        eWishartCh eW; eW.set_parameters(wM/100,100);
+        mat mea=zeros(2,2);
+        mat Ch(2,2);
+        for (int i=0;i<100;i++){Ch=eW.sample_mat(); mea+=Ch.T()*Ch;}
+        cout << mea /100 <<endl;
+        
+        cout << "======= rwiWishart ======== " << endl;
+        rwiWishartCh rwW; rwW.set_parameters(2,0.1,"1 1",0.9);
+        mea=zeros(2,2);
+        mat wMch=chol(wM);
+        for (int i=0;i<100;i++){ 
+                vec tmp=rwW.samplecond(vec(wMch._data(),4));
+                 copy_vector(4,tmp._data(), Ch._data()); 
+                 mea+=Ch.T()*Ch;
+        }
+        cout << mea /100 <<endl;
         //Exit program:
         return 0;