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

Go to the documentation of this file.

#ifndef MX_H
#define MX_H

#include "libBM.h"
#include "libEF.h"
//#include <std>

using namespace itpp;

class emix : public epdf {
protected:
        vec w;
        Array<epdf*> Coms;
        bool destroyComs;
public:
        emix (const RV &rv ) : epdf ( rv ) {};
        void set_parameters ( const vec &w, const Array<epdf*> &Coms, bool copy=true );

        vec sample() const;
        vec mean() const {
                int i; vec mu = zeros ( rv.count() );
                for ( i = 0;i < w.length();i++ ) {mu += w ( i ) * Coms ( i )->mean(); }
                return mu;
        }
        double evalpdflog ( const vec &val ) const {
                int i;
                double sum = 0.0;
                for ( i = 0;i < w.length();i++ ) {sum += w ( i ) * Coms ( i )->evalpdflog ( val );}
                return log ( sum );
        };

//Access methods
        vec& _w() {return w;}
        virtual ~emix(){if (destroyComs){for(int i=0;i<Coms.length();i++){delete Coms(i);}}}
        void ownComs(){destroyComs=true;}
};

class mprod: public compositepdf, public mpdf {
protected:
        Array<epdf*> epdfs;
        Array<ivec> rvcinds;
public:
        mprod ( Array<mpdf*> mFacs): compositepdf(mFacs), mpdf(getrv(true),RV()), epdfs(n), rvcinds(n)
        {       setrvc(rv,rvc);
                setrvcinrv(rvc,rvcinds);
                setindices(rv); 
        for(int i=0;i<n;i++){epdfs(i)=&(mpdfs(i)->_epdf());}
        };

        double evalpdflog ( const vec &val ) const {
                int i;
                double res = 0.0;
                for ( i = n - 1;i > 0;i++ ) {
                        if ( rvcinds ( i ).length() > 0 )
                                {mpdfs ( i )->condition ( val ( rvcinds ( i ) ) );}
                        // add logarithms
                        res += epdfs ( i )->evalpdflog ( val ( rvsinrv ( i ) ) );
                }
                return res;
        }
        vec samplecond ( const vec &cond, double &ll ) {
                vec smp=zeros ( rv.count() );
                vec condi;
                vec smpi;
                ll = 0;
                for ( int i = ( n - 1 );i >= 0;i-- ) {
                        if (( rvcinds ( i ).length() > 0 )||( rvcsinrv ( i ).length() > 0 )) {
                                condi = zeros ( rvcsinrv(i).length() + rvcinds(i).length() );
                                // copy data in condition
                                set_subvector ( condi,rvcinds ( i ), cond );
                                // copy data in already generated sample
                                set_subvector ( condi,rvinrvcs ( i ), get_vec(smp,rvcsinrv(i)) );

                                mpdfs ( i )->condition ( condi );
                        }
                        smpi = epdfs ( i )->sample();
                        // copy contribution of this pdf into smp
                        set_subvector ( smp,rvsinrv ( i ), smpi );
                        // add ith likelihood contribution
                        ll+=epdfs ( i )->evalpdflog ( smpi );
                }
                return smp;
        }
        mat samplecond ( const vec &cond, vec &ll, int N ) {
                mat Smp ( rv.count(),N );
                for ( int i=0;i<N;i++ ) {Smp.set_col ( i,samplecond ( cond,ll ( i ) ) );}
                return Smp;
        }

        ~mprod() {};
};

class eprod: public epdf {
protected:
        Array<const epdf*> epdfs;
        Array<ivec> rvinds;
public:
        eprod ( const Array<const epdf*> epdfs0 ) : epdf ( RV() ),epdfs ( epdfs0 ),rvinds ( epdfs.length() ) {
                bool independent=true;
                for ( int i=0;i<epdfs.length();i++ ) {
                        independent=rv.add ( epdfs ( i )->_rv() );
                        it_assert_debug ( independent==true, "eprod:: given components are not independent ." );
                        rvinds ( i ) = ( epdfs ( i )->_rv() ).dataind ( rv );
                }
        }

        vec mean() const {
                vec tmp ( rv.count() );
                for ( int i=0;i<epdfs.length();i++ ) {
                        vec pom = epdfs ( i )->mean();
                        set_subvector ( tmp,rvinds ( i ), pom );
                }
                return tmp;
        }
        vec sample() const {
                vec tmp ( rv.count() );
                for ( int i=0;i<epdfs.length();i++ ) {
                        vec pom = epdfs ( i )->sample();
                        set_subvector ( tmp,rvinds ( i ), pom );
                }
                return tmp;
        }
        double evalpdflog ( const vec &val ) const {
                double tmp=0;
                for ( int i=0;i<epdfs.length();i++ ) {
                        tmp+=epdfs(i)->evalpdflog(val(rvinds(i)));
                }
                return tmp;
        }
        const epdf* operator () (int i) const {it_assert_debug(i<epdfs.length(),"wrong index");return epdfs(i);}
};


class mmix : public mpdf {
protected:
        Array<mpdf*> Coms;
        emix Epdf;
public:
        mmix ( RV &rv, RV &rvc ) : mpdf ( rv, rvc ), Epdf ( rv ) {ep = &Epdf;};
        void set_parameters ( const vec &w, const Array<mpdf*> &Coms ) {
                Array<epdf*> Eps ( Coms.length() );

                for ( int i = 0;i < Coms.length();i++ ) {
                        Eps ( i ) = & ( Coms ( i )->_epdf() );
                }
                Epdf.set_parameters ( w, Eps );
        };

        void condition ( const vec &cond ) {
                for ( int i = 0;i < Coms.length();i++ ) {Coms ( i )->condition ( cond );}
        };
};
#endif //MX_H

---