root/bdm/stat/libEF.h @ 163

/*!
  \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 <itpp/itbase.h>
#include "../math/libDC.h"
#include "libBM.h"
#include "../itpp_ext.h"
//#include <std>

using namespace itpp;


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

/*!
* \brief General conjugate exponential family posterior density.

* More?...
*/

class eEF : public epdf {
public:
//      eEF() :epdf() {};
        //! default constructor
        eEF ( const RV &rv ) :epdf ( rv ) {};
        //! logarithm of the normalizing constant, \f$\mathcal{I}\f$ 
        virtual double lognc()const =0;
        //!TODO decide if it is really needed
        virtual void tupdate ( double phi, mat &vbar, double nubar ) {};
        //!TODO decide if it is really needed
        virtual void dupdate ( mat &v,double nu=1.0 ) {};
};

/*!
* \brief Exponential family model.

* More?...
*/

class mEF : public mpdf {

public:
        //! Default constructor
        mEF ( const RV &rv0, const RV &rvc0 ) :mpdf ( rv0,rvc0 ) {};
};

/*!
* \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;
        //! dimension (redundant from rv.count() for easier coding )
        int dim;
public:
//      enorm() :eEF() {};
        //!Default constructor
        enorm ( RV &rv );
        //! Set mean value \c mu and covariance \c R
        void set_parameters ( const vec &mu,const sq_T &R );
        //! tupdate in exponential form (not really handy)
        void tupdate ( double phi, mat &vbar, double nubar );
        //! dupdate in exponential form (not really handy)
        void dupdate ( mat &v,double nu=1.0 );

        vec sample() const;
        //! TODO is it used?
        mat sample ( int N ) const;
        double eval ( const vec &val ) const ;
        double evalpdflog ( const vec &val ) const;
        double lognc () const;
        vec mean() const {return mu;}

//Access methods
        //! returns a pointer to the internal mean value. Use with Care!
        vec& _mu() {return mu;}
        
        //! access function
        void set_mu(const vec mu0) { mu=mu0;}

        //! returns pointers to the internal variance and its inverse. Use with Care!
        sq_T& _R() {return R;}

        //! access method
        mat getR () {return R.to_mat();}
};

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

* More?...
*/
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;
public:
        //!Default constructor
        egiw(RV rv, mat V0, double nu0): eEF(rv), V(V0), nu(nu0) {
                it_assert_debug(rv.count()==V.rows(),"Incompatible V0.");
        }

        vec sample() const;
        vec mean() const;
        double evalpdflog ( const vec &val ) const;
        double lognc () const;

        //Access
        //! returns a pointer to the internal statistics. Use with Care!
        ldmat& _V() {return V;}
        //! returns a pointer to the internal statistics. Use with Care!
        double& _nu() {return nu;}

};

/*!
 \brief Gamma posterior density

 Multvariate 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 :
        //! Default constructor
        egamma ( const RV &rv ) :eEF ( rv ) {};
        //! Sets parameters
        void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;};
        vec sample() const;
        //! TODO: is it used anywhere?
//      mat sample ( int N ) const;
        double evalpdflog ( const vec &val ) const;
        double lognc () const;
        //! Returns poiter to alpha and beta. Potentially dengerous: use with care!
        void _param ( vec* &a, vec* &b ) {a=&alpha;b=&beta;};
        vec mean() const {vec pom ( alpha ); pom/=beta; return pom;}
};
/*
//! Weighted mixture of epdfs with external owned components.
class emix : public epdf {
protected:
        int n;
        vec &w;
        Array<epdf*> Coms;
public:
//! Default constructor
        emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
        void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
        vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
        vec sample() {it_error ( "Not implemented" );return 0;}
};
*/

//!  Uniform distributed density on a rectangular support

class euni: public epdf {
protected:
//! lower bound on support
        vec low;
//! upper bound on support
        vec high;
//! internal
        vec distance;
//! normalizing coefficients
        double nk;
//! cache of log( \c nk )
        double lnk;
public:
        //! Defualt constructor
        euni ( const RV rv ) :epdf ( rv ) {}
        double eval ( const vec &val ) const  {return nk;}
        double evalpdflog ( const vec &val ) const  {return lnk;}
        vec sample() const {
                vec smp ( rv.count() ); 
                #pragma omp critical
                UniRNG.sample_vector ( rv.count(),smp );
                return low+elem_mult(distance,smp);
        }
        //! set values of \c low and \c high
        void set_parameters ( const vec &low0, const vec &high0 ) {
                distance = high0-low0;
                it_assert_debug ( min ( distance ) >0.0,"bad support" );
                low = low0;
                high = high0;
                nk = prod ( 1.0/distance );
                lnk = log ( nk );
        }
        vec mean() const {vec pom=high; pom-=low; pom/=2.0; return pom;}
};


/*!
 \brief Normal distributed linear function with linear function of mean value;

 Mean value \f$mu=A*rvc\f$.
*/
template<class sq_T>
class mlnorm : public mEF {
        //! Internal epdf that arise by conditioning on \c rvc
        enorm<sq_T> epdf;
        mat A;
        vec& _mu; //cached epdf.mu;
public:
        //! Constructor
        mlnorm ( RV &rv,RV &rvc );
        //! Set \c A and \c R
        void set_parameters ( const  mat &A, const sq_T &R );
        //!Generate one sample of the posterior
        vec samplecond ( vec &cond, double &lik );
        //!Generate matrix of samples of the posterior
        mat samplecond ( vec &cond, vec &lik, int n );
        //! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
        void condition ( vec &cond );
};

/*!
\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 ( const RV &rv,const RV &rvc );
        //! Set value of \c k
        void set_parameters ( double k );
        void condition ( const vec &val ) {*_beta=k/val;};
};

/*!
\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 ( const RV &rv,const RV &rvc ) : mgamma ( rv,rvc ),refl ( rv.count() ) {};
        //! Set value of \c k
        void set_parameters ( double k0 , vec ref0, double l0 ) {
                mgamma::set_parameters ( k0 );
                refl=pow ( ref0,1.0-l0 );l=l0;
        };

        void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); *_beta=k/mean;};
};

//! 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:
        //! Default constructor
        eEmp ( const RV &rv0 ,int n0 ) :epdf ( rv0 ),n ( n0 ),w ( n ),samples ( n ) {};
        //! Set sample
        void set_parameters ( const vec &w0, epdf* pdf0 );
        //! Potentially dangerous, use with care.
        vec& _w()  {return w;};
        //! access function
        Array<vec>& _samples() {return samples;};
        //! Function performs resampling, i.e. removal of low-weight samples and duplication of high-weight samples such that the new samples represent the same density.
        ivec resample ( RESAMPLING_METHOD method = SYSTEMATIC );
        //! inherited operation : NOT implemneted
        vec sample() const {it_error ( "Not implemented" );return 0;}
        //! inherited operation : NOT implemneted
        double evalpdflog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
        vec mean() const {
                vec pom=zeros ( rv.count() );
                for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
                return pom;
        }
};


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

template<class sq_T>
enorm<sq_T>::enorm ( RV &rv ) :eEF ( rv ), mu ( rv.count() ),R ( rv.count() ),dim ( rv.count() ) {};

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;
};

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 );
        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++ ) {
                NorRNG.sample_vector ( dim,x );
                pom = R.sqrt_mult ( x );
                pom +=mu;
                X.set_col ( i, pom );
        }

        return X;
};

template<class sq_T>
double enorm<sq_T>::eval ( const vec &val ) const {
        double pdfl,e;
        pdfl = evalpdflog ( val );
        e = exp ( pdfl );
        return e;
};

template<class sq_T>
double enorm<sq_T>::evalpdflog ( const vec &val ) const {
        // 1.83787706640935 = log(2pi)
        return  -0.5* (  +R.invqform ( mu-val ) ) - lognc();
};

template<class sq_T>
inline double enorm<sq_T>::lognc () const {
        // 1.83787706640935 = log(2pi)
        return -0.5* ( R.cols() * 1.83787706640935 +R.logdet());
};

template<class sq_T>
mlnorm<sq_T>::mlnorm ( RV &rv0,RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu(epdf._mu()) { ep =&epdf;
}

template<class sq_T>
void mlnorm<sq_T>::set_parameters ( const mat &A0, const sq_T &R0 ) {
        epdf.set_parameters ( zeros ( rv.count() ),R0 );
        A = A0;
}

template<class sq_T>
vec mlnorm<sq_T>::samplecond ( 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 ( 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>
void mlnorm<sq_T>::condition ( vec &cond ) {
        _mu = A*cond;
//R is already assigned;
}

///////////


#endif //EF_H