root/bdm/stat/libEF.cpp @ 32

#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 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; //will be added
        int i;

        for ( i=0; i<rv.count(); i++ ) {
                res += ( alpha ( i ) - 1 ) *std::log ( val ( i ) ) + alpha ( i ) *std::log ( beta ( i ) ) - beta ( i ) *val ( i ) - lgamma ( alpha ( 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 ) {
        *_beta=k / 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();}
}