root/bdm/stat/libEF.cpp @ 96

#include <itpp/itbase.h>
#include <itpp/base/bessel.h>
#include "libEF.h"
#include <math.h>

using namespace itpp;

Uniform_RNG UniRNG;
Normal_RNG NorRNG;
Gamma_RNG GamRNG;

using std::cout;

vec egiw::sample() const {
        it_warning("Function not implemented");
        return vec_1(0.0);
}

double egiw::evalpdflog( const vec &val ) const {
        int nPsi = rv.count()-1; // assuming 1dim y
        double k = nu + nPsi + 2;
        
        double r = val(nPsi); //last entry!
        vec Psi(nPsi+1);
        Psi(0) = -1.0;
        Psi.set_subvector(1,val); // fill the rest
        
        return -0.5*( k*log(r) + V.qform(Psi)) - lognc();
}

double egiw::lognc() const{
        const vec& D = V._D();
        int nPsi = D.length()-1; // assuming 1dim y

// log(2) = 0.693147180559945286226763983
// log(pi) = 1.144729885849400163877476189
return lgamma(0.5*nu) + 0.5*((1.0-nu)*log(D(0)) - V.logdet() + (nu+nPsi)*0.693147180559945286226763983 + nPsi*1.144729885849400163877476189);
}

vec egiw::mean() const {
        const mat &L= V._L();
        const vec &D= V._D();
        
        int end = L.rows()-1;
        vec L0 = L.get_col(0);
        
        vec m(D.length());
        mat iLsub = ltuinv(L(1,end,1,end));
        m.set_subvector(0,iLsub*L0(1,end));
        m(end)= D(0)/(nu-2.0);
        
        return m;
}

vec egamma::sample() const {
        vec smp ( rv.count() );
        int i;

        for ( i=0; i<rv.count(); i++ ) {
                GamRNG.setup ( alpha ( i ),beta ( i ) );
                smp ( i ) = GamRNG();
        }

        return smp;
}

mat egamma::sample ( int N ) const {
        mat Smp ( rv.count(),N );
        int i,j;

        for ( i=0; i<rv.count(); i++ ) {
                GamRNG.setup ( alpha ( i ),beta ( i ) );

                for ( j=0; j<N; j++ ) {
                        Smp ( i,j ) = GamRNG();
                }
        }

        return Smp;
}

double egamma::evalpdflog ( const vec &val ) const {
        double res = 0.0; //the rest will be added
        int i;

        for ( i=0; i<rv.count(); i++ ) {
                res += ( alpha ( i ) - 1 ) *std::log ( val ( i ) ) - beta ( i ) *val ( i );
        }

        return res-lognc();
}

double egamma::lognc() const {
        double res = 0.0; //will be added
        int i;

        for ( i=0; i<rv.count(); i++ ) {
                res += lgamma ( alpha ( i ) ) - alpha ( i ) *std::log ( beta ( i ) ) ;
        }

        return res;
}

//MGamma
mgamma::mgamma ( const RV &rv,const RV &rvc ) : mEF ( rv,rvc ), epdf ( rv ) {vec* tmp; epdf._param ( tmp,_beta );};

void mgamma::set_parameters ( double k0 ) {
        k=k0;
        ep = &epdf;
        epdf.set_parameters ( k*ones ( rv.count() ),*_beta );
};

vec mgamma::samplecond ( vec &cond, double &ll ) {
        this->condition(cond );
        vec smp = epdf.sample();
        ll = epdf.evalpdflog ( smp );
        return smp;
};

//Fixme repetition of mlnorm.samplecond.
mat mgamma::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;
};

ivec eEmp::resample ( RESAMPLING_METHOD method ) {
        ivec ind=zeros_i ( n );
        ivec N_babies = zeros_i ( n );
        vec cumDist = cumsum ( w );
        vec u ( n );
        int i,j,parent;
        double u0;

        switch ( method ) {
        case MULTINOMIAL:
                u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n );

                for ( i = n - 2;i >= 0;i-- ) {
                        u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) );
                }

                break;

        case STRATIFIED:

                for ( i = 0;i < n;i++ ) {
                        u ( i ) = ( i + UniRNG.sample() ) / n;
                }

                break;

        case SYSTEMATIC:
                u0 = UniRNG.sample();

                for ( i = 0;i < n;i++ ) {
                        u ( i ) = ( i + u0 ) / n;
                }

                break;

        default:
                it_error ( "PF::resample(): Unknown resampling method" );
        }

        // U is now full
        j = 0;

        for ( i = 0;i < n;i++ ) {
                while ( u ( i ) > cumDist ( j ) ) j++;

                N_babies ( j ) ++;
        }

        // We have assigned new babies for each Particle
        // Now, we fill the resulting index such that:
        // * particles with at least one baby should not move *
        // This assures that reassignment can be done inplace;

        // find the first parent;
        parent=0; while ( N_babies ( parent ) ==0 ) parent++;

        // Build index
        for ( i = 0;i < n;i++ ) {
                if ( N_babies ( i ) > 0 ) {
                        ind ( i ) = i;
                        N_babies ( i ) --; //this index was now replicated;
                } else {
                        // test if the parent has been fully replicated
                        // if yes, find the next one
                        while ( ( N_babies ( parent ) ==0 ) || ( N_babies ( parent ) ==1 && parent>i ) ) parent++;

                        // Replicate parent
                        ind ( i ) = parent;

                        N_babies ( parent ) --; //this index was now replicated;
                }

        }

        // copy the internals according to ind
        for ( i=0;i<n;i++ ) {
                if ( ind ( i ) !=i ) {
                        samples ( i ) =samples ( ind ( i ) );
                }
                w ( i ) = 1.0/n;
        }

        return ind;
}

void eEmp::set_parameters ( const vec &w0, epdf* epdf0 ) {
//it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
        w=w0;
        w/=sum ( w0 );//renormalize
        n=w.length();
        samples.set_size ( n );

        for ( int i=0;i<n;i++ ) {samples ( i ) =epdf0->sample();}
}