work/git/mixpp/bdm/stat/libEF.h

Go to the documentation of this file.

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


extern Uniform_RNG UniRNG;
extern Normal_RNG NorRNG;
extern Gamma_RNG GamRNG;

class eEF : public epdf {
public:
//      eEF() :epdf() {};
        eEF ( const RV &rv ) :epdf ( rv ) {};
        virtual double lognc() const =0;
        virtual void dupdate ( mat &v ) {it_error ( "Not implemented" );};
        virtual double evallog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
        virtual double evallog ( const vec &val ) const {double tmp;tmp= evallog_nn ( val )-lognc();it_assert_debug(std::isfinite(tmp),"Infinite value"); return tmp;}
        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();
        }
        virtual void pow ( double p ) {it_error ( "Not implemented" );};
};

class mEF : public mpdf {

public:
        mEF ( const RV &rv0, const RV &rvc0 ) :mpdf ( rv0,rvc0 ) {};
};

class BMEF : public BM {
protected:
        double frg;
        double last_lognc;
public:
        BMEF ( const RV &rv, double frg0=1.0 ) :BM ( rv ), frg ( frg0 ) {}
        BMEF ( const BMEF &B ) :BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
        virtual void set_statistics ( const BMEF* BM0 ) {it_error ( "Not implemented" );};
        virtual void bayes ( const vec &data, const double w ) {};
        //original Bayes
        void bayes ( const vec &dt );
        virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
//      virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}

        BMEF* _copy_ ( bool changerv=false ) {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
};

template<class sq_T>
class mlnorm;

template<class sq_T>
class enorm : public eEF {
protected:
        vec mu;
        sq_T R;
        int dim;
public:
        enorm ( const RV &rv );
        void set_parameters ( const vec &mu,const sq_T &R );
        //void tupdate ( double phi, mat &vbar, double nubar );
        void dupdate ( mat &v,double nu=1.0 );

        vec sample() const;
        mat sample ( int N ) const;
//      double eval ( const vec &val ) 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;
//Access methods
        vec& _mu() {return mu;}

        void set_mu ( const vec mu0 ) { mu=mu0;}

        sq_T& _R() {return R;}
        const sq_T& _R() const {return R;}

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

class egiw : public eEF {
protected:
        ldmat V;
        double nu;
        int xdim;
        int nPsi;
public:
        egiw ( RV rv, mat V0, double nu0=-1.0 ) : eEF ( rv ), V ( V0 ), nu ( nu0 ) {
                xdim = rv.count() /V.rows();
                it_assert_debug ( rv.count() ==xdim*V.rows(),"Incompatible V0." );
                nPsi = V.rows()-xdim;
                //set mu to have proper normalization and 
                if (nu0<0){
                        nu = 0.1 +nPsi +2*xdim +2; // +2 assures finite expected value of R
                        // terms before that are sufficient for finite normalization
                }
        }
        egiw ( RV rv, ldmat V0, double nu0=-1.0 ) : eEF ( rv ), V ( V0 ), nu ( nu0 ) {
                xdim = rv.count() /V.rows();
                it_assert_debug ( rv.count() ==xdim*V.rows(),"Incompatible V0." );
                nPsi = V.rows()-xdim;
                if (nu0<0){
                        nu = 0.1 +nPsi +2*xdim +2; // +2 assures finite expected value of R
                        // terms before that are sufficient for finite normalization
                }
        }

        vec sample() const;
        vec mean() const;
        vec variance() const{it_error("Not implemented"); return vec(0);};
        void mean_mat ( mat &M, mat&R ) const;
        double evallog_nn ( const vec &val ) const;
        double lognc () const;

        //Access
        ldmat& _V() {return V;}
        const ldmat& _V() const {return V;}
        double& _nu()  {return nu;}
        const double& _nu() const {return nu;}
        void pow ( double p ) {V*=p;nu*=p;};
};

class eDirich: public eEF {
protected:
        vec beta;
        double gamma;
public:
        eDirich ( const RV &rv, const vec &beta0 ) : eEF ( rv ),beta ( beta0 ) {it_assert_debug ( rv.count() ==beta.length(),"Incompatible statistics" ); };
        eDirich ( const eDirich &D0 ) : eEF ( D0.rv ),beta ( D0.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));}
        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;
        };
        vec& _beta()  {return beta;}
        void set_parameters ( const vec &beta0 ) {
                if ( beta0.length() !=beta.length() ) {
                        it_assert_debug ( rv.length() ==1,"Undefined" );
                        rv.set_size ( 0,beta0.length() );
                }
                beta= beta0;
                gamma = sum(beta);
        }
};

class multiBM : public BMEF {
protected:
        eDirich est;
        vec &beta;
public:
        multiBM ( const RV &rv, const vec beta0 ) : BMEF ( rv ),est ( rv,beta0 ),beta ( est._beta() ) {if(beta.length()>0){last_lognc=est.lognc();}else{last_lognc=0.0;}}
        multiBM ( const multiBM &B ) : BMEF ( B ),est ( rv,B.beta ),beta ( est._beta() ) {}
        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& _epdf() const {return est;};
        const eDirich* _e() const {return &est;};
        void set_parameters ( const vec &beta0 ) {
                est.set_parameters ( beta0 );
                rv = est._rv();
                if ( evalll ) {last_lognc=est.lognc();}
        }
};

class egamma : public eEF {
protected:
        vec alpha;
        vec beta;
public :
        egamma ( const RV &rv ) :eEF ( rv ), alpha(rv.count()), beta(rv.count()) {};
        void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;};
        vec sample() const;
//      mat sample ( int N ) const;
        double evallog ( const vec &val ) const;
        double lognc () const;
        void _param ( vec* &a, vec* &b ) {a=&alpha;b=&beta;};
        vec mean() const {return elem_div(alpha,beta);}
        vec variance() const {return elem_div(alpha,elem_mult(beta,beta)); }
};

class eigamma : public eEF {
        protected:
                vec* alpha;
                vec* beta;
                egamma eg;
        public :
                eigamma ( const RV &rv ) :eEF ( rv ), eg(rv) {eg._param(alpha,beta);};
                void set_parameters ( const vec &a, const vec &b ) {*alpha=a,*beta=b;};
                vec sample() const {return 1.0/eg.sample();};
//      mat sample ( int N ) const;
                double evallog ( const vec &val ) const {return eg.evallog(val);};
                double lognc () const {return eg.lognc();};
                void _param ( vec* &a, vec* &b ) {a=alpha;b=beta;};
                vec mean() const {return elem_div(*beta,*alpha-1);}
                vec variance() const {vec mea=mean(); return elem_div(elem_mult(mea,mea),*alpha-2);}
};
/*
class emix : public epdf {
protected:
        int n;
        vec &w;
        Array<epdf*> Coms;
public:
        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;}
};
*/


class euni: public epdf {
protected:
        vec low;
        vec high;
        vec distance;
        double nk;
        double lnk;
public:
        euni ( const RV rv ) :epdf ( rv ) {}
        double eval ( const vec &val ) const  {return nk;}
        double evallog ( 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 );
        }
        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 {return (high-low)/2.0;}
        vec variance() const {return (pow(high,2)+pow(low,2)+elem_mult(high,low))/3.0;}
};


template<class sq_T>
class mlnorm : public mEF {
protected:
        enorm<sq_T> epdf;
        mat A;
        vec mu_const;
        vec& _mu; //cached epdf.mu;
public:
        mlnorm ( const RV &rv, const RV &rvc );
        void set_parameters ( const  mat &A, const vec &mu0, const sq_T &R );
//      //!Generate one sample of the posterior
//      vec samplecond (const vec &cond, double &lik );
//      //!Generate matrix of samples of the posterior
//      mat samplecond (const vec &cond, vec &lik, int n );
        void condition ( const vec &cond );

        vec& _mu_const() {return mu_const;}
        mat& _A() {return A;}
        mat _R() {return epdf._R().to_mat();}

        template<class sq_M>
        friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M> &ml );
};

class mlstudent : public mlnorm<ldmat> {
protected:
        ldmat Lambda;
        ldmat &_R;
        ldmat Re;
public:
        mlstudent ( const RV &rv0, const RV &rvc0 ) :mlnorm<ldmat> ( rv0,rvc0 ),
                        Lambda ( rv0.count() ),
                        _R ( epdf._R() ) {}
        void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0) {
                epdf.set_parameters ( zeros ( rv.count() ),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!!!!!!
        };

};
class mgamma : public mEF {
protected:
        egamma epdf;
        double k;
        vec* _beta;

public:
        mgamma ( const RV &rv,const RV &rvc ): mEF ( rv,rvc ), epdf ( rv ) {vec* tmp; epdf._param ( tmp,_beta );ep=&epdf;};
        void set_parameters ( double k );
        void condition ( const vec &val ) {*_beta=k/val;};
};

class migamma : public mEF {
        protected:
                eigamma epdf;
                double k;
                vec* _beta;
                vec* _alpha;

        public:
                migamma ( const RV &rv,const RV &rvc ): mEF ( rv,rvc ), epdf ( rv ) {epdf._param ( _alpha,_beta );ep=&epdf;};
                void set_parameters ( double k0 ){k=k0;*_alpha=1.0/(k*k)+2;};
                void condition ( const vec &val ) {
                        *_beta=elem_mult(val,(*_alpha-1));
                };
};

class mgamma_fix : public mgamma {
protected:
        double l;
        vec refl;
public:
        mgamma_fix ( const RV &rv,const RV &rvc ) : mgamma ( rv,rvc ),refl ( rv.count() ) {};
        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;};
};


class migamma_fix : public migamma {
        protected:
                double l;
                vec refl;
        public:
                migamma_fix ( const RV &rv,const RV &rvc ) : migamma ( rv,rvc ),refl ( rv.count() ) {};
                void set_parameters ( double k0 , vec ref0, double l0 ) {
                        migamma::set_parameters ( k0 );
                        refl=pow ( ref0,1.0-l0 );l=l0;
                };

                void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); migamma::condition(mean);};
};
enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
class eEmp: public epdf {
protected :
        int n;
        vec w;
        Array<vec> samples;
public:
        eEmp ( const RV &rv0 ,int n0 ) :epdf ( rv0 ),n ( n0 ),w ( n ),samples ( n ) {};
        void set_parameters ( const vec &w0, const epdf* pdf0 );
        void set_samples ( const epdf* pdf0 );
        void set_n ( int n0, bool copy=true ){w.set_size(n0,copy);samples.set_size(n0,copy);};
        vec& _w()  {return w;};
        const vec& _w() const {return w;};
        Array<vec>& _samples() {return samples;};
        const Array<vec>& _samples() const {return samples;};
        ivec resample ( RESAMPLING_METHOD method = SYSTEMATIC );
        vec sample() const {it_error ( "Not implemented" );return 0;}
        double evallog ( 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;
        }
        vec variance() const {
                vec pom=zeros ( rv.count() );
                for ( int i=0;i<n;i++ ) {pom+=pow(samples ( i ),2) *w ( i );}
                return pom-pow(mean(),2);
        }
};



template<class sq_T>
enorm<sq_T>::enorm ( const 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 = 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>
mlnorm<sq_T>::mlnorm ( const RV &rv0, const 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 vec &mu0, const sq_T &R0 ) {
        epdf.set_parameters ( zeros ( rv.count() ),R0 );
        A = A0;
        mu_const = mu0;
}

// 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>
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 {
        ivec irvn = rvn.dataind ( rv );

        sq_T Rn ( R,irvn );
        enorm<sq_T>* tmp = new enorm<sq_T> ( rvn );
        tmp->set_parameters ( mu ( irvn ), Rn );
        return tmp;
}

template<class sq_T>
mpdf* enorm<sq_T>::condition ( const RV &rvn ) const {

        RV rvc = rv.subt ( rvn );
        it_assert_debug ( ( rvc.count() +rvn.count() ==rv.count() ),"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.count()-1;
        int end=R.rows()-1;
        mat S11 = S.get ( 0,n, 0, n );
        mat S12 = S.get ( 0, n , rvn.count(), end );
        mat S22 = S.get ( rvn.count(), end, rvn.count(), 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> ( rvn,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;
};

#endif //EF_H

---