merger.h

Go to the documentation of this file.

#ifndef MER_H
#define MER_H


#include "mixef.h"

namespace bdm{
using std::string;

class merger : public compositepdf, public epdf {
protected:
        MixEF Mix;
        Array<datalink_m2e*> dls;
        Array<RV> rvzs;
        Array<datalink_m2e*> zdls;

        int Ns;
        int Nc;
        double beta;
        eEmp eSmp;

public:
        merger ( const Array<mpdf*> &S ) :
                        compositepdf ( S ), epdf ( getrv ( false ) ),
                        Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ), dls ( n ), rvzs ( n ), zdls ( n ), eSmp(rv,0) {
                RV ztmp;
                // Extend rv by rvc!
                RV rvc; setrvc ( rv,rvc );
                rv.add ( rvc );
                for ( int i=0;i<n;i++ ) {
                        //Establich connection between mpdfs and merger
                        dls ( i ) = new datalink_m2e ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), rv );
                        // find out what is missing in each mpdf
                        ztmp= mpdfs ( i )->_rv();
                        ztmp.add ( mpdfs ( i )->_rvc() );
                        rvzs ( i ) =rv.subt ( ztmp );
                        zdls ( i ) = new datalink_m2e ( rvzs ( i ), ztmp, rv ) ;
                };
                //Set Default values of parameters
                beta=2.0;
                Ns=100;
                Nc=10;
                Mix.set_method ( EM );
        }
        void set_parameters ( double beta0, int Ns0, int Nc0 ) {beta=beta0;Ns=Ns0;Nc=Nc0;eSmp.set_n(Ns0,false);}
        void init() { 
                Array<vec> Smps ( n );
                //Gibbs sampling
                for ( int i=0;i<n;i++ ) {Smps ( i ) =zeros ( 0 );}
        }
        void merge ( const epdf* g0 );
        void merge () {merge ( & ( Mix._epdf() ) );};

        vec lognorm_merge ( mat &lW );
        vec sample ( ) const { return Mix._epdf().sample();}
        double evallog ( const vec &dt ) const {
                vec dtf=ones ( dt.length() +1 );
                dtf.set_subvector ( 0,dt );
                return Mix.logpred ( dtf );
        }
        vec mean() const {
                const Vec<double> &w = eSmp._w();
                const Array<vec> &S = eSmp._samples();
                vec tmp=zeros ( rv.count() ); 
                for ( int i=0; i<Ns; i++ ) {
                        tmp+=w ( i ) *S ( i );
                }
                return tmp;
        }
        mat covariance() const {
                const vec &w = eSmp._w();
                const Array<vec> &S = eSmp._samples();
                
                vec mea = mean(); 
                
                cout << sum(w) << "," << w*w <<endl;
                
                mat Tmp=zeros(rv.count(), rv.count());
                for ( int i=0; i<Ns; i++ ) {
                        Tmp+=w ( i ) *outer_product(S ( i ), S(i));
                }
                return Tmp-outer_product(mea,mea);
        }
        vec variance() const {
                const vec &w = eSmp._w();
                const Array<vec> &S = eSmp._samples();
                
                vec tmp=zeros(rv.count());
                for ( int i=0; i<Ns; i++ ) {
                        tmp+=w ( i ) *pow(S ( i ),2);
                }
                return tmp-pow(mean(),2);
        }
        virtual ~merger() {
                for ( int i=0; i<n; i++ ) {
                        delete dls ( i );
                        delete zdls ( i );
                }
        };

        MixEF& _Mix() {return Mix;}
        eEmp& _Smp() {return eSmp;}
};

}

#endif // MER_H

---