root/library/bdm/stat/exp_family.h @ 1383

/*!
  \file
  \brief Probability distributions for Exponential Family models.
  \author Vaclav Smidl.

  -----------------------------------
  BDM++ - C++ library for Bayesian Decision Making under Uncertainty

  Using IT++ for numerical operations
  -----------------------------------
*/

#ifndef EF_H
#define EF_H


#include "../shared_ptr.h"
#include "../base/bdmbase.h"
#include "../math/chmat.h"

namespace bdm {


//! Global Uniform_RNG
extern Uniform_RNG UniRNG;
//! Global Normal_RNG
extern Normal_RNG NorRNG;
//! Global Gamma_RNG
extern Gamma_RNG GamRNG;


/*!
* \brief Abstract class of 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;

    //!Evaluate normalized log-probability
    virtual double evallog_nn ( const vec &val ) const NOT_IMPLEMENTED(0);

    //!Evaluate normalized log-probability
    virtual double evallog ( const vec &val ) const {
        double tmp;
        tmp = evallog_nn ( val ) - lognc();
        return tmp;
    }
    //!Evaluate normalized log-probability for many samples
    virtual vec evallog_mat ( 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();
    }
    //!Evaluate normalized log-probability for many samples
    virtual vec evallog_mat ( const Array<vec> &Val ) const {
        vec x ( Val.length() );
        for ( int i = 0; i < Val.length(); i++ ) {
            x ( i ) = evallog_nn ( Val ( i ) ) ;
        }
        return x - lognc();
    }

    //!Power of the density, used e.g. to flatten the density
    virtual void pow ( double p ) NOT_IMPLEMENTED_VOID;
};


//! Estimator for Exponential family
class BMEF : public BM {
public:
    //! forgetting factor
    double frg;
protected:
    //! cached value of lognc() in the previous step (used in evaluation of \c ll )
    double last_lognc;
    //! factor k = [0..1] for scheduling of forgetting factor: \f$ frg_t = (1-k) * frg_{t-1} + k \f$, default 0
    double frg_sched_factor;
public:
    //! Default constructor (=empty constructor)
    BMEF ( double frg0 = 1.0 ) : BM (), frg ( frg0 ), last_lognc(0.0),frg_sched_factor(0.0) {}
    //! Copy constructor
    BMEF ( const BMEF &B ) : BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ),frg_sched_factor(B.frg_sched_factor) {}
    //!get statistics from another model
    virtual void set_statistics ( const BMEF* BM0 ) NOT_IMPLEMENTED_VOID;

    //! Weighted update of sufficient statistics (Bayes rule)
    virtual void bayes_weighted ( const vec &data, const vec &cond = empty_vec, const double w = 1.0 ) {
        if (frg_sched_factor>0) {
            frg = frg*(1-frg_sched_factor)+frg_sched_factor;
        }
    };
    //original Bayes
    void bayes ( const vec &yt, const vec &cond = empty_vec );

    //!Flatten the posterior according to the given BMEF (of the same type!)
    virtual void flatten ( const BMEF * B, double weight=1.0 ) NOT_IMPLEMENTED_VOID;;


    void to_setting ( Setting &set ) const
    {
        BM::to_setting( set );
        UI::save(frg, set, "frg");
        UI::save( frg_sched_factor, set, "frg_sched_factor" );
    }

   /*! Create object from the following structure

    \code
    class = 'BMEF';
     --- optional fields ---
    frg = [];                   % forgetting factor
    frg_sched_factor = [];      % factor for scheduling of forgetting factor: a number from [0..1]
    --- inherited fields ---
    bdm::BM::from_setting
    \endcode
    If the optional fields are not given, they will be filled as follows:
    \code
    frg = 1;                    % default forgetting factor
    frg_sched_factor = 0;
    \endcode
    */
    void from_setting( const Setting &set) {
        BM::from_setting(set);
        if ( !UI::get ( frg, set, "frg" ) )
            frg = 1.0;
        if ( !UI::get ( frg_sched_factor, set, "frg_sched_factor" ) )
            frg_sched_factor = 0.0;
    }

    void validate() {
        BM::validate();
    }

};

/*! \brief Dirac delta density with predefined transformation

Density of the type:\f[ f(x_t | y_t) = \delta (x_t - g(y_t)) \f]
where \f$ x_t \f$ is the \c rv, \f$ y_t \f$ is the \c rvc and g is a deterministic transformation of class fn.
*/
class mgdirac: public pdf {
protected:
    shared_ptr<fnc> g;
public:
    vec samplecond(const vec &cond) {
        bdm_assert_debug(cond.length()==g->dimensionc(),"given cond in not compatible with g");
        vec tmp = g->eval(cond);
        return tmp;
    }
    double evallogcond ( const vec &yt, const vec &cond ) {
        return std::numeric_limits< double >::max();
    }

    /*! Create object from the following structure

    \code
    class = 'mgdirac';
    g = function bdm::fnc;      % any offspring of fnc, bdm::fnc::from_setting
    --- inherited fields ---
    bdm::pdf::from_setting
    \endcode
    */
    void from_setting(const Setting& set);
    void to_setting(Setting &set) const;
    void validate();
};
UIREGISTER(mgdirac);


template<class sq_T, template <typename> class TEpdf>
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 );
    /*! Create Normal density
    \f[ f(rv) = N(\mu, R) \f]
    from structure
    \code
    class = 'enorm<ldmat>', (OR) 'enorm<chmat>', (OR) 'enorm<fsqmat>';
    mu    = [];                  // mean value
    R     = [];                  // variance, square matrix of appropriate dimension
    \endcode
    */
    void from_setting ( const Setting &root );
    void to_setting ( Setting &root ) const ;

    void validate();
    //!@}

    //! \name Mathematical operations
    //!@{

    //! dupdate in exponential form (not really handy)
    void dupdate ( mat &v, double nu = 1.0 );

    //! evaluate bhattacharya distance
    double bhattacharyya(const enorm<sq_T> &e2) {
        bdm_assert(dim == e2.dimension(), "enorms of differnt dimensions");
        sq_T P=R;
        P.add(e2._R());

        double tmp = 0.125*P.invqform(mu - e2._mu()) + 0.5*(P.logdet() - 0.5*(R.logdet() + e2._R().logdet()));
        return tmp;
    }

    vec sample() const;

    double evallog_nn ( const vec &val ) const;
    double lognc () const;
    vec mean() const {
        return mu;
    }
    vec variance() const {
        return diag ( R.to_mat() );
    }
    mat covariance() const {
        return R.to_mat();
    }
    //    mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
    shared_ptr<pdf> condition ( const RV &rvn ) const;

    // target not typed to mlnorm<sq_T, enorm<sq_T> > &
    // because that doesn't compile (perhaps because we
    // haven't finished defining enorm yet), but the type
    // is required
    void condition ( const RV &rvn, pdf &target ) const;

    shared_ptr<epdf> marginal ( const RV &rvn ) const;
    void marginal ( const RV &rvn, enorm<sq_T> &target ) const;
    //!@}

    //! \name Access to attributes
    //!@{

    vec& _mu() {
        return mu;
    }
    const vec& _mu() const {
        return mu;
    }
    void set_mu ( const vec mu0 ) {
        mu = mu0;
    }
    sq_T& _R() {
        return R;
    }
    const sq_T& _R() const {
        return R;
    }
    //!@}

};
UIREGISTER2 ( enorm, chmat );
SHAREDPTR2 ( enorm, chmat );
UIREGISTER2 ( enorm, ldmat );
SHAREDPTR2 ( enorm, ldmat );
UIREGISTER2 ( enorm, fsqmat );
SHAREDPTR2 ( enorm, fsqmat );

//! \class bdm::egauss
//!\brief Gaussian (Normal) distribution. Same as enorm<fsqmat>.
typedef enorm<ldmat> egauss;
UIREGISTER(egauss);


//forward declaration
class mstudent;

/*! distribution of multivariate Student t density

Based on article by Genest and Zidek,
*/
template<class sq_T>
class estudent : public eEF {
protected:
    //! mena value
    vec mu;
    //! matrix H
    sq_T H;
    //! degrees of freedom
    double delta;
public:
    double evallog_nn(const vec &val) const {
        double tmp = -0.5*H.logdet() - 0.5*(delta + dim) * log(1+ H.invqform(val - mu)/delta);
        return tmp;
    }
    double lognc() const {
        //log(pi) = 1.14472988584940
        double tmp = -lgamma(0.5*(delta+dim))+lgamma(0.5*delta) + 0.5*dim*(log(delta) + 1.14472988584940);
        return tmp;
    }
    void marginal (const RV &rvm, estudent<sq_T> &marg) const {
        ivec ind = rvm.findself_ids(rv); // indices of rvm in rv
        marg._mu() = mu(ind);
        marg._H() = sq_T(H,ind);
        marg._delta() = delta;
        marg.validate();
    }
    shared_ptr<epdf> marginal(const RV &rvm) const {
        shared_ptr<estudent<sq_T> > tmp = new estudent<sq_T>;
        marginal(rvm, *tmp);
        return tmp;
    }
    vec sample() const {
                enorm<sq_T> en;
                en._mu()=mu;
                en._R()=H;
                en._R()*=delta/(delta-2);
                en.validate();
                return en.sample(); 
        }

    vec mean() const {
        return mu;
    }
    mat covariance() const {
        return delta/(delta-2)*H.to_mat();
    }
    vec variance() const {
        return diag(covariance());
    }
    //! \name access
    //! @{
    //! access function
    vec& _mu() {
        return mu;
    }
    //! access function
    sq_T& _H() {
        return H;
    }
    //! access function
    double& _delta() {
        return delta;
    }
    //!@}
    //! todo
    void from_setting(const Setting &set) {
        epdf::from_setting(set);
        mat H0;
        UI::get(H0,set, "H");
        H= H0; // conversion!!
        UI::get(delta,set,"delta");
        UI::get(mu,set,"mu");
    }
    void to_setting(Setting &set) const {
        epdf::to_setting(set);
        UI::save(H.to_mat(), set, "H");
        UI::save(delta, set, "delta");
        UI::save(mu, set, "mu");
    }
    void validate() {
        eEF::validate();
        dim = H.rows();
    }
};
UIREGISTER2(estudent,fsqmat);
UIREGISTER2(estudent,ldmat);
UIREGISTER2(estudent,chmat);

// template < class sq_T, template <typename> class TEpdf = estudent >
// class mstudent : public pdf_internal< TEpdf<sq_T> > {
// protected:
//      //! Internal epdf that arise by conditioning on \c rvc
//      mat A;
//      //! Constant additive term
//      vec mu_const;
//      
// public:
//      void condition( const vec &cond ) {
//              iepdf._mu()=A*val*mu_const;
//      }
//              
// };

/*!
* \brief Gauss-inverse-Wishart density stored in LD form

* For \f$p\f$-variate densities, 
* 
* \f$ M,R \sim GiW( V_t, \nu_t) \propto |R|^{0.5\nu}\exp(-1/2 tr(R^{-1}[I,M] V_t [I;M'])) \f$
*
* Factorizes as:
* \f$ M|R \sim N( \hat{M}, R \otimes Vz^{-1}) \propto |R|^{-0.5dim(psi)}\exp(-1/2 tr((M-\hat{M})R^{-1}(M-\hat{M})Vz) \f$
* \f$ R \sim iW( \Lambda,\delta) \propto |R|^{-0.5(\nu - dim(psi))} |\Lambda|^{}\exp(-1/2 tr(R^{-1}\Lambda) \f$
* where in standard notation \f$ |R|^{-0.5(\delta + p +1)}\f$, i.e. 
* \f$ \delta = \nu-dim(psi) -p-1 \f$
*/
class egiw : public eEF {
    //! \var log_level_enums logvartheta
    //! Log variance of the theta part

    LOG_LEVEL(egiw,logvartheta);

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(),dimx(0) {};
    egiw ( int dimx0, ldmat V0, double nu0 = -1.0 ) : eEF(),dimx(0) {
        set_parameters ( dimx0, V0, nu0 );
        validate();
    };

    void set_parameters ( int dimx0, ldmat V0, double nu0 = -1.0 );
    //!@}

    vec sample() const;
    mat sample_mat ( int n ) const;
    vec mean() const;
    vec variance() const;
    //mat covariance() const;
    void sample_mat ( mat &Mi, chmat &Ri ) const;

    void factorize ( mat &M, ldmat &Vz, ldmat &Lam ) const;
    //! LS estimate of \f$\theta\f$
    vec est_theta() const;

    //! Covariance of the LS estimate
    ldmat est_theta_cov() const;

    //! expected values of the linear coefficient and the covariance matrix are written to \c M and \c R , respectively
    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;
    };

    //! marginal density (only student for now)
    shared_ptr<epdf> marginal(const RV &rvm) const {
        bdm_assert(dimx==1, "Not supported");
        //TODO - this is too trivial!!!
        ivec ind = rvm.findself_ids(rv);
        if (min(ind)==0) { //assume it si
            shared_ptr<estudent<ldmat> > tmp = new estudent<ldmat>;
            mat M;
            ldmat Vz;
            ldmat Lam;
            factorize(M,Vz,Lam);

            tmp->_mu() = M.get_col(0);
            ldmat H;
            Vz.inv(H);
            H *=Lam._D()(0)/nu;
            tmp->_H() = H;
            tmp->_delta() = nu;
            tmp->validate();
            return tmp;
        }
        return NULL;
    }
    //! \name Access attributes
    //!@{

    ldmat& _V() {
        return V;
    }
    const ldmat& _V() const {
        return V;
    }
    double& _nu()  {
        return nu;
    }
    const double& _nu() const {
        return nu;
    }
    const int & _dimx() const {
        return dimx;
    }

    /*! Create object from the following structure
    \code

    class = 'egiw';
    dimx    = [...];       % dimension of the wishart part
    V.L     = [...];       % L part of matrix V
    V.D     = [...];       % D part of matrix V
    -or- fV = [...];       % full matrix V
    -or- dV = [...];       % vector of diagonal of V (when V not given)

    rv = RV({'names',...},[sizes,...],[times,...]);   % description of RV
    rvc = RV({'names',...},[sizes,...],[times,...]);  % description of RV in condition

    --- optional fields ---
    nu      = [];             % scalar \nu ((almost) degrees of freedom)
    --- inherited fields ---
    bdm::eEF::from_setting
    \endcode

    fulfilling formula \f[ f(rv) = GiW(V,\nu) \f]

    If \nu is not given, it will be computed to obtain proper pdf.

    \sa log_level_enums
    */
    void from_setting ( const Setting &set );
    //! see egiw::from_setting
    void to_setting ( Setting& set ) const;
    void validate();
    void log_register ( bdm::logger& L, const string& prefix );

    void log_write() const;
    //!@}
};
UIREGISTER ( egiw );
SHAREDPTR ( egiw );

//! \brief Gauss-Wishart with recursion on moments
//! Using precision as parameter
//! following notation of [Karny Andrysek 2009], precision
template<class sq_T>
class egw_ls: public eEF{
        public:
                vec theta;
                sq_T P;
                double omega;
                double nu;
                
                vec mean() const{
                        return concat(theta, omega);
                }
                mat covariance() const {
                        sq_T tmp=P;
                        tmp*=nu/((nu-2)*omega);
                        return tmp.to_mat();//<======= error - missing omega
                }
                vec variance() const {
                        return diag(covariance());//<======= error - missing omega
                }
                vec sample() const NOT_IMPLEMENTED(vec(0));
                double lognc() const {return 0.0;} //TODO
};

/*! \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;
public:
    //!\name Constructors
    //!@{

    eDirich () : eEF () {};
    eDirich ( const eDirich &D0 ) : eEF () {
        set_parameters ( D0.beta );
        validate();
    };
    eDirich ( const vec &beta0 ) {
        set_parameters ( beta0 );
        validate();
    };
    void set_parameters ( const vec &beta0 ) {
        beta = beta0;
        dim = beta.length();
    }
    //!@}

    //! using sampling procedure from wikipedia
    vec sample() const {
        vec y ( beta.length() );
        for ( int i = 0; i < beta.length(); i++ ) {
            GamRNG.setup ( beta ( i ), 1 );
#pragma omp critical
            y ( i ) = GamRNG();
        }
        return y / sum ( y );
    }

    vec mean() const {
        return beta / sum ( beta );
    };
    vec variance() const {
        double gamma = sum ( beta );
        return elem_mult ( beta, ( gamma - beta ) ) / ( gamma*gamma* ( gamma + 1 ) );
    }
    //! In this instance, val is ...
    double evallog_nn ( const vec &val ) const {
        double tmp;
        tmp = ( beta - 1 ) * log ( val );
        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 );
        return tmp;
    }

    //!access function
    vec& _beta()  {
        return beta;
    }

    /*! Create object from the following structure
    \code
    class = 'eDirich';
    beta  = [...];           % vector parameter beta
    --- inherited fields ---
    bdm::eEF::from_setting
    \endcode
    */
    void from_setting ( const Setting &set );

    void validate();

    void to_setting ( Setting &set ) const;
};
UIREGISTER ( eDirich );

/*! \brief Product of Beta distributions

Continuous Dirichlet density of \f$n\f$-dimensional variable \f$x\f$
\f[
f(x|\alpha,\beta)  = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\, x^{\alpha-1}(1-x)^{\beta-1}
\f]
is a simplification of Dirichlet to univariate case.
*/
class eBeta: public eEF {
public:
    //!sufficient statistics
    vec alpha;
    //!sufficient statistics
    vec beta;
public:
    //!\name Constructors
    //!@{

    eBeta () : eEF () {};
    eBeta ( const eBeta &B0 ) : eEF (), alpha(B0.alpha),beta(B0.beta) {
        validate();
    };
    //!@}

    //! using sampling procedure from wikipedia
    vec sample() const {
        vec y ( beta.length() ); // all vectors
        for ( int i = 0; i < beta.length(); i++ ) {
            GamRNG.setup ( alpha ( i ), 1 );
#pragma omp critical
            double Ga = GamRNG();

            GamRNG.setup ( beta ( i ), 1 );
#pragma omp critical
            double Gb = GamRNG();

            y ( i ) = Ga/(Ga+Gb);
        }
        return y;
    }

    vec mean() const {
        return elem_div(alpha, alpha + beta); // dot-division
    };
    vec variance() const {
        vec apb=alpha+beta;
        return elem_div (elem_mult ( alpha, beta) ,
                         elem_mult ( elem_mult(apb,apb), apb+1 ) );
    }
    //! In this instance, val is ...
    double evallog_nn ( const vec &val ) const {
        double tmp;
        tmp = ( alpha - 1 ) * log ( val ) + (beta-1)*log(1-val);
        return tmp;
    }

    double lognc () const {
        double lgb = 0.0;
        for ( int i = 0; i < beta.length(); i++ ) {
            lgb += -lgamma ( alpha(i)+beta(i) ) + lgamma(alpha(i)) + lgamma(beta(i));
        }
        return lgb;
    }

    /*! Create object from the following structure

    \code
    class = 'eBeta';
    alpha = [...];           % vector parameter alpha
    beta  = [...];           % vector parameter beta of the same length as alpha
    \endcode

    Class does not call bdm::eEF::from_setting
    */
    void from_setting ( const Setting &set ) {
        UI::get(alpha, set, "alpha", UI::compulsory);
        UI::get(beta, set, "beta", UI::compulsory);
    }

    void validate() {
        bdm_assert(alpha.length()==beta.length(), "eBeta:: alpha and beta length do not match");
        dim = alpha.length();
    }

    void to_setting ( Setting &set ) const {
        UI::save(alpha, set, "alpha");
        UI::save(beta, set, "beta");
    }
};
UIREGISTER ( eBeta );

/*! \brief Random Walk on Dirichlet

Using simple assignment
 \f[ \beta = rvc / k + \beta_c \f]
 hence, mean value = rvc, variance = (k+1)*mean*mean;

 The greater k is, the greater is the variance of the random walk;

 \f$ \beta_c \f$ is used as regularizing element to avoid corner cases, i.e. when one element of rvc is zero.
 By default is it set to 0.1;
*/
class mDirich: public pdf_internal<eDirich> {
protected:
    //! constant \f$ k \f$ of the random walk
    double k;
    //! cache of beta_i
    vec &_beta;
    //! stabilizing coefficient \f$ \beta_c \f$
    vec betac;
public:
    mDirich() : pdf_internal<eDirich>(), _beta ( iepdf._beta() ) {};
    void condition ( const vec &val ) {
        _beta =  val / k + betac;
    };

    /*! Create object from the following structure
    \code
    class = 'mDirich';
    k = 1;                      % multiplicative constant k
    --- optional ---
    beta0 = [...];              % initial values of beta
    betac = [...];              % initial values of beta stabilizing coefficients
    --- inherited fields ---
    bdm::pdf::from_setting
    \endcode
    fulfilling form  \f[ f(rv|rvc) = Di(rvc*k) \f]

    If the optional fields are not given, they will be filled as follows:
    \code
    beta0 = [1,1,1,...];                
    betac = 0.1 * [1,1,1,...];
    \endcode
    */
    void from_setting ( const Setting &set );
    void     to_setting  (Setting  &set) const;
    void validate();
};
UIREGISTER ( mDirich );

/*! \brief Random Walk with vector Beta distribution

Using simple assignment
\f{eqnarray*}
\alpha & = & rvc / k + \beta_c \\
\beta & = &(1-rvc) / k + \beta_c \\
\f}
for each element of alpha and beta, mean value = rvc, variance = (k+1)*mean*mean;

The greater \f$ k \f$ is, the greater is the variance of the random walk;

\f$ \beta_c \f$ is used as regularizing element to avoid corner cases, i.e. when one element of rvc is zero.
By default is it set to 0.1;
*/
class mBeta: public pdf_internal<eBeta> {
    //! vector of constants \f$ k \f$ of the random walk
    vec k;
    //! stabilizing coefficient \f$ \beta_c \f$
    vec betac;

public:
    void condition ( const vec &val ) {
        this->iepdf.alpha =  elem_div(val , k) + betac;
        this->iepdf.beta =  elem_div (1-val , k) + betac;
    };

    /*! Create object from the following structure
    \code
    class = 'mBeta';
    k = [...];          % vector of multiplicative constants k
    --- optional fields ---
    beta  = [...];      % initial values of beta
    betac = [...];      % initial values of beta stabilizing constants
    --- inherited fields ---
    bdm::pdf::from_setting
    \endcode
    fulfilling form \f[ f(rv|rvc) = \prod Beta(rvc,k) \f]

    If the optional fields are not given, they will be filled as follows:
    \code
    beta  = [1,1,1,...];
    betac = 0.1 * [1,1,1,...];
    \endcode
    
    */
    void from_setting ( const Setting &set );

    void to_setting  (Setting  &set) const;

    void validate() {
        pdf_internal<eBeta>::validate();
        bdm_assert(betac.length()==dimension(),"Incomaptible betac");
        bdm_assert(k.length()==dimension(),"Incomaptible k");
        dimc = iepdf.dimension();
    }
    //!
};
UIREGISTER(mBeta);

//! \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 &yt, const vec &cond = empty_vec );

    double logpred ( const vec &yt ) const;

    void flatten ( const BMEF* B , double weight);

    //! return correctly typed posterior (covariant return)
    const eDirich& posterior() const {
        return est;
    };
    //! constructor function
    void set_parameters ( const vec &beta0 ) {
        est.set_parameters ( beta0 );
        est.validate();
        if ( evalll ) {
            last_lognc = est.lognc();
        }
    }

    void to_setting ( Setting &set ) const {
        BMEF::to_setting ( set );
        UI::save( &est, set, "prior" );
    }

    /*! Create object from the following structure

    \code
    class = 'MultiBM';
    prior = configuration of bdm::eDirich;       % any offspring of eDirich, bdm::eDirich::from_setting
    --- inherited fields ---
    bdm::BMEF::from_setting
    \endcode
    */
    void from_setting (const Setting &set )  {
        BMEF::from_setting ( set );
        UI::get( est, set, "prior" );
    }
};
UIREGISTER( multiBM );

/*!
 \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 );
        validate();
    };
    void set_parameters ( const vec &a, const vec &b ) {
        alpha = a, beta = b;
    };
    //!@}

    vec sample() const;
        double evallog_nn ( const vec &val ) const;
        double lognc () const;
    //! Returns pointer to internal alpha. Potentially dengerous: use with care!
    vec& _alpha() {
        return alpha;
    }
    //! Returns pointer to internal beta. Potentially dengerous: use with care!
    vec& _beta() {
        return beta;
    }
    vec mean() const {
        return elem_div ( alpha, beta );
    }
    vec variance() const {
        return elem_div ( alpha, elem_mult ( beta, beta ) );
    }

    /*! Create object from the following structure

    \code
    class = 'egamma';
    alpha = [...];         % vector alpha
    beta = [...];          % vector beta
    --- inherited fields ---
    bdm::eEF::from_setting
    \endcode
     fulfilling formula \f[ f(rv|rvc) = \Gamma(\alpha, \beta) \f]
    */
    void from_setting ( const Setting &set );

    void to_setting ( Setting &set ) const;
    void validate();
};
UIREGISTER ( egamma );
SHAREDPTR ( egamma );

/*!
 \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;};
};
*/

//! \brief 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;
//! 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;
        low = low0;
        high = high0;
        nk = prod ( 1.0 / distance );
        lnk = log ( nk );
    }
    //!@}

    double evallog ( const vec &val ) const  {
        if ( any ( val < low ) && any ( val > high ) ) {
            return -inf;
        } else 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;
    }


     /*! Create object from the following structure
     \code
    
     class = 'euni'
     high = [...];          % vector of upper bounds
     low = [...];           % vector of lower bounds     
      rv = RV({'names',...},[sizes,...],[times,...]);   % description of RV
     --- inherited fields ---
     bdm::epdf::from_setting
     \endcode

     fulfilling form \f[ f(rv) = U(low,high) \f]
     */
    void from_setting ( const Setting &set );
    void     to_setting  (Setting  &set) const;
    void validate();
};
UIREGISTER ( euni );

//! Uniform density with conditional mean value
class mguni : public pdf_internal<euni> {
    //! function of the mean value
    shared_ptr<fnc> mean;
    //! distance from mean to both sides
    vec delta;
public:
    void condition ( const vec &cond ) {
        vec mea = mean->eval ( cond );
        iepdf.set_parameters ( mea - delta, mea + delta );
    }

    /*! Create object from the following structure
    \code
    class = 'mguni';
    mean = function bdm::fnc;       % any offspring of fnc, bdm::fnc::from_setting
    delta = [...];                  % distance from mean to both sides
    --- inherited fields ---
    bdm::pdf::from_setting
    \endcode
    */
    void from_setting ( const Setting &set ) {
        pdf::from_setting ( set ); //reads rv and rvc
        UI::get ( delta, set, "delta", UI::compulsory );
        mean = UI::build<fnc> ( set, "mean", UI::compulsory );
        iepdf.set_parameters ( -delta, delta );
    }
    void     to_setting  (Setting  &set) const {
        pdf::to_setting ( set );
        UI::save( iepdf.mean(), set, "delta");
        UI::save(mean, set, "mean");
    }
    void validate() {
        pdf_internal<euni>::validate();
        dimc = mean->dimensionc();

    }

};
UIREGISTER ( mguni );
/*!
 \brief Normal distributed linear function with linear function of mean value;

 Mean value \f$ \mu=A*\mbox{rvc}+\mu_0 \f$.
*/
template < class sq_T, template <typename> class TEpdf = enorm >
class mlnorm : public pdf_internal< TEpdf<sq_T> > {
protected:
        //! Internal epdf that arise by conditioning on \c rvc
        mat A;
        //! Constant additive term
        vec mu_const;
        //            vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
public:
    //! \name Constructors
    //!@{
    mlnorm() : pdf_internal< TEpdf<sq_T> >() {};
    mlnorm ( const mat &A, const vec &mu0, const sq_T &R ) : pdf_internal< TEpdf<sq_T> >() {
        set_parameters ( A, mu0, R );
        validate();
    }

    //! Set \c A and \c R
    void set_parameters ( const  mat &A0, const vec &mu0, const sq_T &R0 ) {
        this->iepdf.set_parameters ( zeros ( A0.rows() ), R0 );
        A = A0;
        mu_const = mu0;
    }

    //!@}
    //! 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 ) {
        this->iepdf._mu() = A * cond + mu_const;
//R is already assigned;
    }

    //!access function
    const vec& _mu_const() const {
        return mu_const;
    }
    //!access function
    const mat& _A() const {
        return A;
    }
    //!access function
    mat _R() const {
        return this->iepdf._R().to_mat();
    }
    //!access function
    sq_T __R() const {
        return this->iepdf._R();
    }

    //! Debug stream
    template<typename sq_M>
    friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M, enorm> &ml );

    /*! Create Normal density with linear function of mean value
     \f[ f(rv|rvc) = N(A*rvc+const, R) \f]
    from structure
    \code
    class = 'mlnorm<ldmat>', (OR) 'mlnorm<chmat>', (OR) 'mlnorm<fsqmat>';
    A     = [];                  // matrix or vector of appropriate dimension
    R     = [];                  // square matrix of appropriate dimension
    --- optional ---
    const = zeros(A.rows);       // vector of constant term
    \endcode
    */
    void from_setting ( const Setting &set ) {
        pdf::from_setting ( set );

        UI::get ( A, set, "A", UI::compulsory );
        UI::get ( mu_const, set, "const", UI::optional);
        mat R0;
        UI::get ( R0, set, "R", UI::compulsory );
        set_parameters ( A, mu_const, R0 );
    }

    void to_setting (Setting &set) const {
        pdf::to_setting(set);
        UI::save ( A, set, "A");
        UI::save ( mu_const, set, "const");
        UI::save ( _R(), set, "R");
    }

    void validate() {
        pdf_internal<TEpdf<sq_T> >::validate();
        if (mu_const.length()==0) { // default in from_setting
            mu_const=zeros(A.rows());
        }
        bdm_assert ( A.rows() == mu_const.length(), "mlnorm: A vs. mu mismatch" );
        bdm_assert ( A.rows() == _R().rows(), "mlnorm: A vs. R mismatch" );
        this->dimc = A.cols();

    }
};
UIREGISTER2 ( mlnorm, ldmat );
SHAREDPTR2 ( mlnorm, ldmat );
UIREGISTER2 ( mlnorm, fsqmat );
SHAREDPTR2 ( mlnorm, fsqmat );
UIREGISTER2 ( mlnorm, chmat );
SHAREDPTR2 ( mlnorm, chmat );

//! \class mlgauss
//!\brief Normal distribution with linear function of mean value. Same as mlnorm<fsqmat>.
typedef mlnorm<fsqmat> mlgauss;
UIREGISTER(mlgauss);

//! pdf with general function for mean value
template<class sq_T>
class mgnorm : public pdf_internal< enorm< sq_T > > {
private:
//            vec &mu; WHY NOT?
    shared_ptr<fnc> g;

public:
    //!default constructor
    mgnorm() : pdf_internal<enorm<sq_T> >() { }
    //!set mean function
    inline void set_parameters ( const shared_ptr<fnc> &g0, const sq_T &R0 );
    inline void condition ( const vec &cond );


    /*! Create Normal density with given function of mean value
    \f[ f(rv|rvc) = N( g(rvc), R) \f]
    from structure
    \code
    class = 'mgnorm';
    g.class =  'fnc';      // function for mean value evolution
    g._fields_of_fnc = ...;

    R = [1, 0;             // covariance matrix
        0, 1];
        --OR --
    dR = [1, 1];           // diagonal of cavariance matrix

    rv = RV({'name'})      // description of RV
    rvc = RV({'name'})     // description of RV in condition
    \endcode
    */


    void from_setting ( const Setting &set ) {
        pdf::from_setting ( set );
        shared_ptr<fnc> g = UI::build<fnc> ( set, "g", UI::compulsory );

        mat R;
        vec dR;
        if ( UI::get ( dR, set, "dR" ) )
            R = diag ( dR );
        else
            UI::get ( R, set, "R", UI::compulsory );

        set_parameters ( g, R );
        //validate();
    }


    void to_setting  (Setting  &set) const {
        UI::save( g,set, "g");
        UI::save(this->iepdf._R().to_mat(),set, "R");

    }



    void validate() {
        this->iepdf.validate();
        bdm_assert ( g->dimension() == this->iepdf.dimension(), "incompatible function" );
        this->dim = g->dimension();
        this->dimc = g->dimensionc();
        this->iepdf.validate();
    }

};

UIREGISTER2 ( mgnorm, chmat );
UIREGISTER2 ( mgnorm, ldmat );
SHAREDPTR2 ( mgnorm, chmat );


/*! \brief (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, enorm> {
protected:
    //! Variable \f$ \Lambda \f$ from theory
    ldmat Lambda;
    //! Reference to variable \f$ R \f$
    ldmat &_R;
    //! Variable \f$ R_e \f$
    ldmat Re;
public:
    mlstudent () : mlnorm<ldmat, enorm> (),
        Lambda (),    _R ( iepdf._R() ) {}
    //! constructor function
    void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 ) {
        iepdf.set_parameters ( mu0, R0 );// was Lambda, why?
        A = A0;
        mu_const = mu0;
        Re = R0;
        Lambda = Lambda0;
    }

    void condition ( const vec &cond );

    void validate() {
        mlnorm<ldmat, enorm>::validate();
        bdm_assert ( A.rows() == mu_const.length(), "mlstudent: A vs. mu mismatch" );
        bdm_assert ( _R.rows() == A.rows(), "mlstudent: A vs. R mismatch" );

    }
};

/*!
\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 pdf_internal<egamma> {
protected:

    //! Constant \f$k\f$
    double k;

    //! cache of iepdf.beta
    vec &_beta;

public:
    //! Constructor
    mgamma() : pdf_internal<egamma>(), k ( 0 ),
        _beta ( iepdf._beta() ) {
    }

    //! Set value of \c k
    void set_parameters ( double k, const vec &beta0 );

    void condition ( const vec &val ) {
        _beta = k / val;
    };

    /*! Create object from the following structure
    \code
    class = 'mgamma';    
    beta = [...];      % vector of initial beta
    k = x;             %  multiplicative scalar constant k
    --- inherited fields ---
    bdm::pdf::from_setting
    \endcode
    fulfilling form \f[ f(rv|rvc) = \Gamma(k, k/rvc) \f]
    */
    void from_setting ( const Setting &set );
    void     to_setting  (Setting  &set) const;
    void validate();
};
UIREGISTER ( mgamma );
SHAREDPTR ( mgamma );

/*!
\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 pdf_internal<eigamma> {
protected:
    //! Constant \f$k\f$
    double k;

    //! cache of iepdf.alpha
    vec &_alpha;

    //! cache of iepdf.beta
    vec &_beta;

public:
    //! \name Constructors
    //!@{
    migamma() : pdf_internal<eigamma>(),
        k ( 0 ),
        _alpha ( iepdf._alpha() ),
        _beta ( iepdf._beta() ) {
    }

    migamma ( const migamma &m ) : pdf_internal<eigamma>(),
        k ( 0 ),
        _alpha ( iepdf._alpha() ),
        _beta ( iepdf._beta() ) {
    }
    //!@}

    //! Set value of \c k
    void set_parameters ( int len, double k0 ) {
        k = k0;
        iepdf.set_parameters ( ( 1.0 / ( k*k ) + 2.0 ) *ones ( len ) /*alpha*/, ones ( len ) /*beta*/ );
    };

    void     validate () {
        pdf_internal<eigamma>::validate();
        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;
    };

    void validate () {
        mgamma::validate();
        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;
    };

    void validate() {
        migamma::validate();
        dimc = dimension();
    };

    void condition ( const vec &val ) {
        vec mean = elem_mult ( refl, pow ( val, l ) );
        migamma::condition ( mean );
    };

    /*! Create object from the following structure
    \code
    class = 'migamma_ref';
    ref = [...];           % reference vector
    l = [];                % constant scalar l
    k = [];                % constant scalar k    
    --- inherited fields ---
    bdm::migamma::from_setting
    \endcode
    fulfilling form \f[ f(rv|rvc) = i\Gamma(k, k/(rvc^l \circ ref^{(1-l)}) \f]
     */
    void from_setting ( const Setting &set );

    void to_setting  (Setting  &set) const;
};


UIREGISTER ( migamma_ref );
SHAREDPTR ( migamma_ref );

/*! \brief Log-Normal probability density - it allows only diagonal covariances!

Density of the form \f$ \log(x)\sim \mathcal{N}(\mu,\sigma^2) \f$ , i.e.
\f[
x \sim \frac{1}{x\sigma\sqrt{2\pi}}\exp{-\frac{1}{2\sigma^2}(\log(x)-\mu)}
\f]

Function from_setting loads mu and R in the same way as it does for enorm<>!
*/
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...

 */
class mlognorm : public pdf_internal<elognorm> {
protected:
    //! parameter 1/2*sigma^2
    double sig2;

    //! access
    vec &mu;
public:
    //! Constructor
    mlognorm() : pdf_internal<elognorm>(),
        sig2 ( 0 ),
        mu ( iepdf._mu() ) {
    }

    //! Set value of \c k
    void set_parameters ( int size, double k ) {
        sig2 = 0.5 * log ( k * k + 1 );
        iepdf.set_parameters ( zeros ( size ), 2*sig2*eye ( size ) );
    };

    void validate() {
        pdf_internal<elognorm>::validate();
        dimc = iepdf.dimension();
    }

    void condition ( const vec &val ) {
        mu = log ( val ) - sig2;//elem_mult ( refl,pow ( val,l ) );
    };

    /*! Create object from the following structure
    \code
    class = 'mlognorm';
    k   = [];               % "variance" k
    mu0 = [];               % initial value of mean
    --- inherited fields ---
    bdm::pdf_internal<elognorm>::from_setting
    \endcode
    fulfilling form \f[ f(rv|rvc) = log\mathcal{N}( \log(rvc)-0.5\log(k^2+1), k I) \f]
    */
    void from_setting ( const Setting &set );

    void to_setting  (Setting  &set) const;
};

UIREGISTER ( mlognorm );
SHAREDPTR ( mlognorm );

/*! \brief Inverse Wishart density defined on Choleski decomposition
 * 
 * here \f$\Omega \sim |\Sigma|^{-0.5(\delta}|\Omega|^{0.5(\delta -p -1)} \exp(-1/2 tr(\Omega\Sigma^{-1}))\$f
*/
class eWishartCh : public epdf {
protected:
    //! Upper-Triagle of Choleski decomposition of \f$ \Psi \f$
    chmat Sigma;
    //! dimension of matrix  \f$ \Psi \f$
    int p;
    //! degrees of freedom  \f$ \nu \f$
    double delta;
public:
    //! Set internal structures
    void set_parameters ( const mat &Y0, const double delta0 ) {
        Sigma = chmat ( Y0 );
        delta = delta0;
        p = Sigma.rows();
    }
    //! Set internal structures
    void set_parameters ( const chmat &Y0, const double delta0 ) {
        Sigma = Y0;
        delta = delta0;
        p = Sigma.rows();
    }

    virtual void validate () {
        epdf::validate();
        dim = p * p;
    }

    //! Sample matrix argument - lower 
    //! Using Bartlet decomposition: W=L A A^T L^T where A is lower triag and L is choleski factor of Sigma.
    chmat sample_mat() const {
        mat A_T = zeros ( p, p ); // A transpose

        //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
            A_T ( 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
                A_T ( i, j ) = NorRNG.sample();
            }
        }
        chmat tmp;
                tmp._Ch()=A_T*Sigma._Ch();
        return tmp;// return upper triangular part of the decomposition
    }

    vec sample () const {
        return vec ( sample_mat().to_mat()._data(), p*p );
    }

    virtual vec mean() const NOT_IMPLEMENTED(0);

    //! return expected variance (not covariance!)
    virtual vec variance() const NOT_IMPLEMENTED(0);

    virtual double evallog ( const vec &val ) const NOT_IMPLEMENTED(0);

    //! fast access function y0 will be copied into Y.Ch.
    void setY ( const mat &Ch0 ) {
        copy_vector ( dim, Ch0._data(), Sigma._Ch()._data() );
    }

    //! fast access function y0 will be copied into Y.Ch.
    void _setY ( const vec &ch0 ) {
        copy_vector ( dim, ch0._data(), Sigma._Ch()._data() );
    }

    //! access function
    const chmat& getY() const {
        return Sigma;
    }
};

//! Inverse Wishart on Choleski decomposition
/*! Being computed by conversion from `standard' Wishart
*/
class eiWishartCh: public epdf {
protected:
    //! Internal instance of Wishart density
    eWishartCh W;
    //! size of Ch
    int p;
    //! parameter delta
    double delta;
public:
    //! constructor function
    void set_parameters ( const mat &Y0, const double delta0 ) {
        delta = delta0;
        W.set_parameters ( inv ( Y0 ), delta0 );
        p = Y0.rows();
    }

    virtual void     validate () {
        epdf::validate();
        W.validate();
        dim = W.dimension();
    }


    vec sample() const {
        mat iCh;
        iCh = inv ( W.sample_mat()._Ch() );
        return vec ( iCh._data(), dim );
    }
    //! access function
    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;
    };

    virtual vec mean() const NOT_IMPLEMENTED(0);

    //! return expected variance (not covariance!)
    virtual vec variance() const NOT_IMPLEMENTED(0);
};

//! Random Walk on inverse Wishart
class rwiWishartCh : public pdf_internal<eiWishartCh> {
protected:
    //!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;
    //! power of the reference
    double l;
    //! dimension
    int p;

public:
    rwiWishartCh() : sqd ( 0 ), l ( 0 ), p ( 0 ) {}
    //! constructor function
    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 );
        iepdf.set_parameters ( eye ( p ), delta );
    };

    void validate() {
        pdf_internal<eiWishartCh>::validate();
        dimc = iepdf.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++;
        }

        iepdf._setY ( sqd*z );
    }
};

//! Switch between various resampling methods.
enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };

//! Shortcut for multinomial.sample(int n). Various simplifications are allowed see RESAMPLING_METHOD
void resample(const vec &w, ivec &ind, RESAMPLING_METHOD=SYSTEMATIC);

/*! \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 () {};
    //! copy constructor
    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 );
    };
    //! Set samples
    void set_parameters ( const Array<vec> &Av ) {
        n = Av.size();
        w = 1 / n * ones ( n );
        samples = Av;
    };
    virtual void     validate ();
    //! 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 vec& _sample ( int i ) const {
        return samples ( i );
    };
    //! 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.
    void resample ( RESAMPLING_METHOD method = SYSTEMATIC );

    //! inherited operation : NOT implemented
    vec sample() const NOT_IMPLEMENTED(0);

    //! inherited operation : NOT implemented
    double evallog ( const vec &val ) const NOT_IMPLEMENTED(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;

    void to_setting ( Setting &set ) const;

    /*! Create object from the following structure

    \code
    class = 'eEmp';
    samples = [...];  % array of samples
    w = [...];        % weights of samples stored in vector
    --- inherited fields ---
    bdm::epdf::from_setting
    \endcode
    */
    void from_setting ( const Setting &set );
};
UIREGISTER(eEmp);


////////////////////////

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;
    validate();
};

template<class sq_T>
void enorm<sq_T>::from_setting ( const Setting &set ) {
    epdf::from_setting ( set ); //reads rv

    UI::get ( mu, set, "mu", UI::compulsory );
    mat Rtmp;// necessary for conversion
    UI::get ( Rtmp, set, "R", UI::compulsory );
    R = Rtmp; // conversion
}

template<class sq_T>
void enorm<sq_T>::validate() {
    eEF::validate();
    bdm_assert ( mu.length() == R.rows(), "mu and R parameters do not match" );
    dim = mu.length();
}

template<class sq_T>
void enorm<sq_T>::to_setting ( Setting &set ) const {
    epdf::to_setting ( set ); //reads rv
    UI::save ( mu, set, "mu");
    UI::save ( R.to_mat(), set, "R");
}



template<class sq_T>
void enorm<sq_T>::dupdate ( mat &v, double nu ) {
    //
};

// template<class sq_T>
// void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
//     //
// };

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>
// double enorm<sq_T>::eval ( const vec &val ) const {
//     double pdfl,e;
//     pdfl = evallog ( val );
//     e = exp ( pdfl );
//     return e;
// };

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>
// vec mlnorm<sq_T>::samplecond (const  vec &cond, double &lik ) {
//     this->condition ( cond );
//     vec smp = epdf.sample();
//     lik = epdf.eval ( smp );
//     return smp;
// }

// template<class sq_T>
// mat mlnorm<sq_T>::samplecond (const 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;
// }


template<class sq_T>
shared_ptr<epdf> enorm<sq_T>::marginal ( const RV &rvn ) const {
    enorm<sq_T> *tmp = new enorm<sq_T> ();
    shared_ptr<epdf> narrow ( tmp );
    marginal ( rvn, *tmp );
    return narrow;
}

template<class sq_T>
void enorm<sq_T>::marginal ( const RV &rvn, enorm<sq_T> &target ) const {
    bdm_assert ( isnamed(), "rv description is not assigned" );
    ivec irvn = rvn.dataind ( rv );

    sq_T Rn ( R, irvn );  // select rows and columns of R

    target.set_rv ( rvn );
    target.set_parameters ( mu ( irvn ), Rn );
}

template<class sq_T>
shared_ptr<pdf> enorm<sq_T>::condition ( const RV &rvn ) const {
    mlnorm<sq_T> *tmp = new mlnorm<sq_T> ();
    shared_ptr<pdf> narrow ( tmp );
    condition ( rvn, *tmp );
    return narrow;
}

template<class sq_T>
void enorm<sq_T>::condition ( const RV &rvn, pdf &target ) const {
    typedef mlnorm<sq_T> TMlnorm;

    bdm_assert ( isnamed(), "rvs are not assigned" );
    TMlnorm &uptarget = dynamic_cast<TMlnorm &> ( target );

    RV rvc = rv.subt ( rvn );
    bdm_assert ( ( 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() );

    uptarget.set_rv ( rvn );
    uptarget.set_rvc ( rvc );
    uptarget.set_parameters ( A, mu1 - A*mu2, R_n );
    uptarget.validate();
}

/*! \brief Dirac delta function distribution */
class dirac: public epdf {
public:
    vec point;
public:
    double evallog (const vec &dt) const {
        return -inf;
    }
    vec mean () const {
        return point;
    }
    vec variance () const {
        return zeros(point.length());
    }
    void qbounds ( vec &lb, vec &ub, double percentage = 0.95 ) const {
        lb = point;
        ub = point;
    }
    //! access
    const vec& _point() {
        return point;
    }
    void set_point(const vec& p) {
        point =p;
        dim=p.length();
    }
    vec sample() const {
        return point;
    }
};


///////////

template<class sq_T>
void mgnorm<sq_T >::set_parameters ( const shared_ptr<fnc> &g0, const sq_T &R0 ) {
    g = g0;
    this->iepdf.set_parameters ( zeros ( g->dimension() ), R0 );
}

template<class sq_T>
void mgnorm<sq_T >::condition ( const vec &cond ) {
    this->iepdf._mu() = g->eval ( cond );
};

//!     odo unify this stuff with to_string()
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._R() << endl;
    return os;
};

}
#endif //EF_H