root/bdm/stat/libEF.h @ 14

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

using namespace itpp;

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

* More?...
*/
class eEF : public epdf {

public:
        virtual void tupdate( double phi, mat &vbar, double nubar ) {};
        virtual void dupdate( mat &v,double nu=1.0 ) {};
};

class mEF : public mpdf {

public:

};

/*!
* \brief General exponential family density

* More?...
*/
template<class sq_T>
class enorm : public eEF {
        int dim;
        vec mu;
        sq_T R;
public:
        enorm( RV &rv, vec &mu, sq_T &R );
        enorm();
        void tupdate( double phi, mat &vbar, double nubar );
        void dupdate( mat &v,double nu=1.0 );
        vec sample();
        mat sample(int N);
        double eval( const vec &val );
        Normal_RNG RNG;
};

/*!
 \brief
*/
template<class sq_T>
class mlnorm : public mEF {
        enorm<sq_T> epdf;
        mat A;
public:
        //! Constructor
        mlnorm( RV &rv,RV &rvc, mat &A, sq_T &R );
        //!Generate one sample of the posterior
        vec samplecond( vec &cond, double &lik );
        mat samplecond( vec &cond, vec &lik, int n );
        void condition( vec &cond );
};

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

template<class sq_T>
enorm<sq_T>::enorm( RV &rv, vec &mu0, sq_T &R0 ) {
        dim = rv.count();
        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() {
        vec x( dim );
        RNG.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 ) {
        mat X( dim,N );
        vec x( dim );
        vec pom;
        int i;
        for ( i=0;i<N;i++ ) {
                RNG.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 ) {
        //
};


template<class sq_T>
enorm<sq_T>::enorm() {};

template<class sq_T>
mlnorm<sq_T>::mlnorm( RV &rv,RV &rvc, mat &A, sq_T &R ) {
        int dim = rv.count();
        vec mu( dim );

        epdf = enorm<sq_T>( rv,mu,R );
}

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<dim; 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 ) {
        epdf.mu = A*cond;
//R is already assigned;
}

#endif //EF_H