root/bdm/estim/merger.cpp @ 193

Revision 193, 4.7 kB (checked in by smidl, 16 years ago)

oprava merger_iter a jeho casti

Line 
1#include <itpp/itbase.h>
2#include "merger.h"
3#include "arx.h"
4
5vec merger::lognorm_merge ( mat &lW ) {
6        int nu=lW.rows();
7        vec mu = sum ( lW ) /nu; //mean of logs
8        vec lam = sum ( pow ( lW,2 ) )-nu*pow ( mu,2 );
9        double coef=0.0;
10        switch ( nu ) {
11                case 2:
12                        coef=sqrt ( beta*2 ) * ( 1-0.5*sqrt ( ( 4*beta-3 ) /beta ) );
13                        return exp ( coef*sqrt ( lam ) + mu );
14                        break;
15                case 3://Ration of Bessel
16                        break;
17                case 4:
18                        break;
19                default: // Approximate conditional density
20                        break;
21        }
22        return vec ( 0 );
23}
24
25void merger::merge ( const epdf* g0 ) {
26        it_file dbg ( "merger_debug.it" );
27
28        it_assert_debug ( rv.equal ( g0->_rv() ),"Incompatible g0" );
29        //Empirical density - samples
30        eEmp eSmp ( rv,Ns );
31        eSmp.set_parameters ( ones ( Ns ), g0 );
32        Array<vec> &Smp = eSmp._samples(); //aux
33        vec &w = eSmp._w(); //aux
34
35        mat Smp_ex =ones ( rv.count() +1,Ns ); // Extended samples for the ARX model - the last row is ones
36        for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );}
37
38        dbg << Name ( "Smp_0" ) << Smp_ex;
39
40        // Stuff for merging
41        vec lw_src ( Ns );
42        vec lw_mix ( Ns );
43        mat lW=zeros ( n,Ns );
44        vec vec0 ( 0 );
45
46        // Initial component in the mixture model
47        mat V0=1e-8*eye ( rv.count() +1 );
48        ARX A0 ( RV ( "{th_r  }", vec_1 ( rv.count() * ( rv.count() +1 ) ) ),\
49                 V0, rv.count() *rv.count() +3.0 ); //initial guess of Mix: zero mean, large variance
50
51        // ============= MAIN LOOP ==================
52        bool converged=false;
53        int niter = 0;
54        char str[100];
55
56        epdf* Mpred;
57        vec Mix_pdf ( Ns );
58        while ( !converged ) {
59                //Re-estimate Mix
60                //Re-Initialize Mixture model
61                Mix.init ( &A0, Smp_ex, Nc );
62                Mix.bayesB ( Smp_ex, w*Ns );
63                Mpred = Mix.predictor ( rv ); // Allocation => must be deleted at the end!!
64
65                if ( 1 ) {
66                        // Generate new samples
67                        eSmp.set_samples ( Mpred );
68                        for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );}
69                }
70                else {
71                        for ( int ii=0;ii<10;ii++ ) {
72                                for ( int jj=0; jj<10; jj++ ) {
73                                        Smp ( ii+jj*10 ) =vec_2 ( -1.0+6*ii/10.0, -1.0+6*jj/10.0 );
74                                }
75                        }
76                        for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );}
77                }
78
79                sprintf ( str,"Mpdf%d",niter );
80                for ( int i=0;i<Ns;i++ ) {Mix_pdf ( i ) = Mix.logpred ( Smp_ex.get_col ( i ) );}
81                dbg << Name ( str ) << Mix_pdf;
82
83                sprintf ( str,"Smp%d",niter );
84                dbg << Name ( str ) << Smp_ex;
85
86                //Importace weighting
87                for ( int i=0;i<n;i++ ) {
88                        lw_src=0.0;
89                        //======== Same RVs ===========
90                        //Split according to dependency in rvs
91                        if ( mpdfs ( i )->_rv().count() ==rv.count() ) {
92                                // no need for conditioning or marginalization
93                                for ( int j=0;j<Ns; j++ ) { // Smp is Array<> => for cycle
94                                        lw_src ( j ) =mpdfs ( i )->_epdf().evalpdflog ( Smp ( j ) );
95                                }
96                        }
97                        else {
98                                // compute likelihood of marginal on the conditional variable
99                                if ( mpdfs ( i )->_rvc().count() >0 ) {
100                                        // Make marginal on rvc_i
101                                        epdf* tmp_marg = Mpred->marginal ( mpdfs ( i )->_rvc() );
102                                        //compute vector of lw_src
103                                        for ( int k=0;k<Ns;k++ ) {
104                                                lw_src ( k ) += tmp_marg->evalpdflog ( dls ( i )->get_val ( Smp ( i ) ) );
105                                        }
106                                        delete tmp_marg;
107                                }
108                                // Compute likelihood of the missing variable
109                                if ( rv.count() > ( mpdfs ( i )->_rv().count() + mpdfs ( i )->_rvc().count() ) ) {
110                                        ///////////////
111                                        cout << Mpred->mean() <<endl;
112                                        // There are variales unknown to mpdfs(i) : rvzs
113                                        mpdf* tmp_cond = Mpred->condition ( rvzs ( i ) );
114                                        // Compute likelihood
115                                        vec lw_dbg=lw_src;
116                                        for ( int k= 0; k<Ns; k++ ) {
117                                                lw_src ( k ) += log (
118                                                                    tmp_cond->evalcond (
119                                                                        zdls ( i )->get_val ( Smp ( k ) ),
120                                                                        zdls ( i )->get_cond ( Smp ( k ) ) ) );
121                                        }
122                                        delete tmp_cond;
123                                }
124                                // Compute likelihood of the partial source
125                                for ( int k= 0; k<Ns; k++ ) {
126                                        mpdfs ( i )->condition ( dls ( i )->get_cond ( Smp ( k ) ) );
127                                        lw_src ( k ) += mpdfs ( i )->_epdf().evalpdflog ( dls ( i )->get_val ( Smp ( k ) ) );
128                                }
129                               
130                        }
131                        lW.set_row ( i, lw_src ); // do not divide by mix
132                }
133                //Importance of the mixture
134                for ( int j=0;j<Ns;j++ ) {
135                        lw_mix ( j ) =Mix.logpred ( Smp_ex.get_col ( j ) );
136                }
137                sprintf ( str,"lW%d",niter );
138                dbg << Name ( str ) << lW;
139
140                w = lognorm_merge ( lW ); //merge
141
142                sprintf ( str,"w%d",niter );
143                dbg << Name ( str ) << w;
144                sprintf ( str,"lw_m%d",niter );
145                dbg << Name ( str ) << lw_mix;
146
147                //Importance weighting
148                w /=exp ( lw_mix ); // hoping that it is not numerically sensitive...
149                //renormalize
150                w /=sum ( w );
151
152                sprintf ( str,"w_is_%d",niter );
153                dbg << Name ( str ) << w;
154
155//              eSmp.resample(); // So that it can be used in bayes
156//              for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );}
157
158                sprintf ( str,"Smp_res%d",niter );
159                dbg << Name ( str ) << Smp;
160
161                // ==== stopping rule ===
162                niter++;
163                converged = ( niter>4 );
164        }
165
166}
Note: See TracBrowser for help on using the browser.