root/library/tests/stat_mixtures_test.cpp @ 384

#include "stat/exp_family.h"
#include "stat/mixtures.h"
using namespace bdm;

//These lines are needed for use of cout and endl
using std::cout;
using std::endl;

double normcoef ( const epdf* ep,const vec &xb, const vec &yb, int Ngr=100 ) {
        mat PPdf ( Ngr+1,Ngr+1 );
        vec rgr ( 2 );

        int i=0,j=0;
        double xstep= ( xb ( 1 )-xb ( 0 ) ) /Ngr;
        double ystep= ( yb ( 1 )-yb ( 0 ) ) /Ngr;

        for ( double x=xb ( 0 );x<=xb ( 1 );x+= xstep,i++ ) {
                rgr ( 0 ) =x;j=0;
                for ( double y=yb ( 0 );y<=yb ( 1 );y+=ystep,j++ ) {
                        rgr ( 1 ) =y;
                        PPdf ( i,j ) =exp ( ep->evallog ( rgr ) );
                }
        }
        return sumsum ( PPdf ) *xstep*ystep;
}

int main() {
        RNG_randomize();

        RV x ( "{x }" );
        RV y ( "{y }" );
        RV xy=concat(x,y);
        
        enorm<ldmat> E1; E1.set_rv ( xy );
                
        it_file iti("xxx");
        mat R; vec mu;
        iti >> Name("R") >>R;
        iti >> Name("mu") >>mu;

        cout << mu << R <<endl;
        E1.set_parameters( mu, R);//"1.00054 1.0455" , mat ( "0.740142 -0.259015; -0.259015 1.0302" ));
        enorm<ldmat> E2;E2.set_rv ( xy );
        E2.set_parameters( "-1.2 -0.1" , mat ( "1 0.4; 0.4 0.5" ));

        Array<epdf*> A1(1);
        A1(0) = &E1;
        
        emix M1; M1.set_rv(xy);
        M1.set_parameters(vec("1"), A1, false);
        cout << "======== test if ARX and emix with one ARX are the same ==="<<endl;
        
        epdf* Mm = M1.marginal(y);
        epdf* Am = E1.marginal(y);
        mpdf* Mc = M1.condition(y);
        mpdf* Ac = E1.condition(y);
        
        cout << *((mlnorm<ldmat>*)Ac) <<endl;
        
        cout << "Mix marg at 0: " << Mm->evallog(vec_1(0.0)) <<endl;
        cout << "ARX marg at 0: " << Am->evallog(vec_1(0.0)) <<endl;
        cout << "Mix cond at 0,0: " << Mc->evallogcond(vec_1(0.0),vec_1(0.0)) <<endl;
        cout << "ARX cond at 0,0: " << Ac->evallogcond(vec_1(0.0),vec_1(0.0)) <<endl;
        
        cout << "======== Mixture with two components ==="<<endl;
        Array<epdf*> A2(2);
        A2(0) = &E1;
        A2(1) = &E2;
        
        emix M2; M2.set_rv(xy);
        M2.set_parameters(vec("1"), A2, false);
        
        
        cout << "Mixture normalization: " << normcoef(&M2, vec("-3 3 "), vec("-3 3 ")) <<endl;
        
        int N=3;
        vec ll ( N );
        vec ll2 ( N );
        mat Smp=M2.sample_m ( N );
        ll = M2.evallog_m(Smp);
        
        vec Emu = sum ( Smp,2 ) /N;
        mat Er = ( Smp*Smp.transpose() ) /N - outer_product ( Emu,Emu );

        cout << "original empirical mean: " <<Emu <<endl;
        cout << "original empirical variance: " <<Er <<endl;

        epdf* Mg = (epdf*)M2.marginal ( y );
        mpdf* Cn = (mpdf*)M2.condition ( x );

        it_file it("emix_test.it");
        it << Name("Smp") << Smp;

        cout<< "marginal mean: " << Mg->mean() <<endl;

        cout << "========== putting them back together ======= "<<endl;
        mepdf mMg ( Mg );
        Array<mpdf*> AA ( 2 );
        AA ( 0 ) = Cn;
        AA ( 1 ) = &mMg;
        mprod mEp ( AA );

        for (int j=0; j<N; j++){
                ll2(j)=mEp.evallogcond (Smp.get_col(j), vec ( 0 ));
        }
        it << Name("ll") << ll;
        it << Name("ll2") << ll2;
        
        
}