Changeset 205

Timestamp:

11/10/08 15:40:30 (15 years ago)

Author:

smidl

Message:

merger posledni verze

Files:

• bdm/estim/merger.cpp (modified) (7 diffs)
• bdm/estim/merger.h (modified) (4 diffs)
• bdm/stat/libEF.h (modified) (2 diffs)
• mpdm/merg_giw.cpp (modified) (5 diffs)
• mpdm/merger_iter_cond.cpp (added)
• tests/CMakeLists.txt (modified) (1 diff)
• tests/merger_2d_test.cpp (added)

bdm/estim/merger.cpp

@@ -16 +16 @@
                 case 3://Ratio of Bessel
                         coef = sqrt ( ( 3*beta-2 ) /3*beta );
-                        return log( besselk ( 0,sq2bl*coef )) - log( besselk ( 0,sq2bl ) ) +  mu;
+                        return log ( besselk ( 0,sq2bl*coef ) ) - log ( besselk ( 0,sq2bl ) ) +  mu;
                         break;
                 case 4:
@@ -31 +31 @@
         it_assert_debug ( rv.equal ( g0->_rv() ),"Incompatible g0" );
         //Empirical density - samples
-        eEmp eSmp ( rv,Ns );
         eSmp.set_parameters ( ones ( Ns ), g0 );
         Array<vec> &Smp = eSmp._samples(); //aux
@@ -62 +61 @@
         char str[100];
 
-        epdf* Mpred;
+        epdf* Mpred=Mix.predictor(rv);
         vec Mix_pdf ( Ns );
         while ( !converged ) {
@@ -69 +68 @@
                 Mix.flatten ( &Mix_init );
                 Mix.bayesB ( Smp_ex, w*Ns );
+                delete Mpred;
                 Mpred = Mix.predictor ( rv ); // Allocation => must be deleted at the end!!
-
-//              sprintf ( str,"Mpred_mean%d",niter );
-//              dbg << Name ( str ) << Mpred->mean();
-
-                if ( 1 ) {
+                
+                if ( 1./sum_sqr ( w ) <0.5*Ns ) {
                         // Generate new samples
                         eSmp.set_samples ( Mpred );
                         for ( int i=0;i<Ns;i++ ) {
                                 //////////// !!!!!!!!!!!!!
-                                if ( Smp ( i ) ( 1 ) <0 ) {Smp ( i ) ( 1 ) = GamRNG(); }
+                                if ( Smp ( i ) ( 1 ) <0 ) {Smp ( i ) ( 1 ) = 0.01; }
                                 set_col_part ( Smp_ex,i,Smp ( i ) );
                         }
-                        {cout<<"Eff. sample size=" << 1./sum_sqr ( w ) << endl;}
+                        {cout<<"Resampling =" << 1./sum_sqr ( w ) << endl;}
                         cout << sum ( Smp_ex,2 ) /Ns <<endl;
+                        cout << Smp_ex*Smp_ex.T()/Ns << endl; 
                 }
-                else
-                        {cout<<"Eff. sample size=" << 1./sum_sqr ( w ) << endl;}
+//              sprintf ( str,"Mpred_mean%d",niter );
+//              dbg << Name ( str ) << Mpred->mean();
+
 
 //              sprintf ( str,"Mpdf%d",niter );
@@ -135 +134 @@
                                                                         zdls ( i )->get_cond ( Smp ( k ) ) ) );
                                                 if ( !std::isfinite ( lw_src ( k ) ) ) {
-                                                        cout << endl;
+                                                        lw_src ( k ) = -1e16; cout << "!";
                                                 }
                                         }
@@ -166 +165 @@
                 //Importance weighting
                 lw -=  lw_mix; // hoping that it is not numerically sensitive...
-                w = exp(lw-max(lw));
+                w = exp ( lw-max ( lw ) );
                 //renormalize
                 double sumw =sum ( w );
@@ -173 +172 @@
                 }
                 else {
-                        
+
                 }
-                {
-                        double eff = 1./sum_sqr ( w );
-                        if ( eff<2 ) {
-                                int mi= max_index ( w );
-                                cout << w(mi) <<endl;
-                                cout << lW.get_col ( mi ) <<endl;
-                                mat mm(2,1);mm.set_col(0,lW.get_col(mi));
-                                cout << lognorm_merge(mm) <<endl;
-                                cout << lw_mix ( mi ) <<endl;
-                                cout << lw ( mi ) <<endl;
-                                cout << Smp ( mi ) <<endl;
-                                cout << Mix._Coms(0)->_e()->mean() <<endl;
-                        }
-                        cout << "Eff: " << eff <<endl;
-                }
+
 //              sprintf ( str,"w_is_%d",niter );
 //              dbg << Name ( str ) << w;

bdm/estim/merger.h

@@ -43 +43 @@
         //!Prior on the log-normal merging model
         double beta;
+        //! Projection to empirical density
+        eEmp eSmp;
+
 public:
 //!Default constructor
         merger ( const Array<mpdf*> &S ) :
                         compositepdf ( S ), epdf ( getrv ( false ) ),
-                        Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ), dls ( n ), rvzs ( n ), zdls ( n ) {
+                        Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ), dls ( n ), rvzs ( n ), zdls ( n ), eSmp(rv,0) {
                 RV ztmp;
                 // Extend rv by rvc!
@@ -68 +71 @@
         }
 //! Set internal parameters used in approximation
-        void set_parameters ( double beta0, int Ns0, int Nc0 ) {        beta=beta0;Ns=Ns0;Nc=Nc0;}
+        void set_parameters ( double beta0, int Ns0, int Nc0 ) {beta=beta0;Ns=Ns0;Nc=Nc0;eSmp.set_n(Ns0,false);}
 //!Initialize the proposal density. This function must be called before merge()!
         void init() { ////////////// NOT FINISHED
@@ -90 +93 @@
                 return Mix.logpred ( dtf );
         }
-        vec mean() const {return Mix._epdf().mean();}
+        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 variance() 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);
+        }
 //! for future use
         virtual ~merger() {
@@ -101 +126 @@
 //! Access function
         MixEF& _Mix() {return Mix;}
+//! Access function
+        eEmp& _Smp() {return eSmp;}
 };
 

bdm/stat/libEF.h

@@ -200 +200 @@
         ldmat& _V() {return V;}
         //! returns a pointer to the internal statistics. Use with Care!
-        double& _nu() {return nu;}
+        const ldmat& _V() const {return V;}
+        //! returns a pointer to the internal statistics. Use with Care!
+        double& _nu()  {return nu;}
+        const double& _nu() const {return nu;}
         void pow ( double p ) {V*=p;nu*=p;};
 };
@@ -520 +523 @@
         //! Set sample
         void set_samples ( const epdf* pdf0 );
+        //! Set sample
+        void set_n ( int n0, bool copy=true ){w.set_size(n0,copy);samples.set_size(n0,copy);};
         //! Potentially dangerous, use with care.
         vec& _w()  {return w;};
+        //! Potentially dangerous, use with care.
+        const vec& _w() const {return w;};
         //! access function
         Array<vec>& _samples() {return samples;};
+        //! access function
+        const Array<vec>& _samples() const {return samples;};
         //! Function performs resampling, i.e. removal of low-weight samples and duplication of high-weight samples such that the new samples represent the same density.
         ivec resample ( RESAMPLING_METHOD method = SYSTEMATIC );

mpdm/merg_giw.cpp

@@ -17 +17 @@
         RV uu=u1; uu.add ( u2 );
 
+        double a1t = 1.5;
+        double a2t = 0.8;
+        double sqr=0.10;
         // Full system
-        vec thg ( "1 1" ); //Simulated system - zero for constant term
-        double sqr=1.0;
+        vec thg =vec_2 ( a1t,a2t ); //Simulated system - zero for constant term
         vec Th = concat ( thg, sqr ); //Full parameter
 
@@ -28 +30 @@
         RV all =a1; all.add ( a2 ); all.add ( r );
         RV allj =a1; allj.add ( r ); allj.add ( a2 );
-        vec Thj("1 1 1");
+        vec Thj=vec_3 ( a1t,sqr,a2t );
         // Setup values
 
@@ -37 +39 @@
         ARX P1 ( concat ( a1,r ), V0, -1 );
         ARX P2 ( concat ( a2,r ), V0, -1 );
-        ARX PG ( all, V0g, -1 );
+        ARX PG ( all, V0g, -1 ); //or -1?
+//      ARX PGk ( all, V0g, -1 );
 
         //Test estimation
@@ -44 +47 @@
 
         // Logging
-        dirfilelog L ( "exp/merg_giw",1 );
+        dirfilelog L ( "exp/merg_giw",ndat );
         int Li_Data = L.add ( RV ( "{Y U1 U2 }" ), "" );
-        int Li_LL   = L.add ( RV ( "{1 2 G M }" ), "LL" );
+        int Li_LL   = L.add ( RV ( "{G M }" ), "LL" );
+        int Li_P1m   = L.add ( RV ( "{a1 r }" ), "P1" );
+        int Li_P2m   = L.add ( RV ( "{a2 r }" ), "P2" );
         int Li_Gm   = L.add ( RV ( "{a1 a2 r }" ), "G" );
         int Li_Mm   = L.add ( RV ( "{a1 r a2 }" ), "M" );
@@ -54 +59 @@
         vec yt ( 1 );
 
-        vec LLs ( 4 );
+        vec LLs ( 2 );
         vec rgrg ( 2 );
-        
+
         //Proposal
-        enorm<ldmat> g0 ( a1 ); g0.set_parameters ( "1  ",mat("1") );
+        enorm<ldmat> g0 ( a1 ); g0.set_parameters ( "1  ",mat ( "1" ) );
         egamma g1 ( r ); g1.set_parameters ( "2  ", "2" );
-        enorm<ldmat> g2 ( a2 ); g2.set_parameters ( "1  ",mat("1") );
+        enorm<ldmat> g2 ( a2 ); g2.set_parameters ( "1  ",mat ( "1" ) );
 
-        Array<const epdf*> A(3); A(0) = &g0; A(1)=&g1; A(2) = &g2;
-        eprod G0(A);
-        
+        Array<const epdf*> A ( 3 ); A ( 0 ) = &g0; A ( 1 ) =&g1; A ( 2 ) = &g2;
+        eprod G0 ( A );
+
         for ( t=0; t<ndat; t++ ) {
                 // True system
                 rgrg ( 0 ) = pow ( sin ( ( t/40.0 ) *pi ),3 );
-                rgrg ( 1 ) = pow ( cos ( ( t/40.0 ) *pi ),3 );
+                rgrg ( 1 ) = pow ( cos ( ( t/40.0 +0.1 ) *pi ),3 );
 
                 Yt ( t ) = thg*rgrg + sqr * NorRNG();
 
                 // Bayes for all
-                P1.bayes ( concat ( Yt ( t ),vec_1(rgrg ( 0 ) )) );
-                P2.bayes ( concat ( Yt ( t ),vec_1(rgrg ( 1 ) )) );
+                P1.bayes ( concat ( Yt ( t ),vec_1 ( rgrg ( 0 ) ) ) );
+                P2.bayes ( concat ( Yt ( t ),vec_1 ( rgrg ( 1 ) ) ) );
                 PG.bayes ( concat ( Yt ( t ),rgrg ) );
 
+
+                // crippling PGk by substituting zeros.
+                /*      ldmat &Vk=const_cast<egiw*>(PGk._e())->_V();  //PG ldmat does not like 0!
+                        mat fVk=PG._e()->_V().to_mat();
+                        fVk(1,2) = 0.0;
+                        fVk(2,1) = 0.0;
+                        Vk = ldmat(fVk);
+                */      //PGk is now krippled
+
                 // Merge estimates
-                mepdf eG1(P1._e());
-                mepdf eG2(P2._e());
+                mepdf eG1 ( P1._e() );
+                mepdf eG2 ( P2._e() );
                 Array<mpdf*> A ( 2 ); A ( 0 ) =&eG1;A ( 1 ) =&eG2;
                 merger M ( A );
-                M.set_parameters ( 10.0, 100,2 );
+                M.set_parameters ( 1.5, 100,3 ); //M._Mix().set_method(QB);
+                //M2.set_parameters ( 100.0, 1000,3 ); //M2._Mix().set_method(QB);
                 M.merge ( &G0 );
+                //M2.merge ( &G0 );
 
                 //Likelihood
                 yt ( 0 ) = Yt ( t );
 
-                LLs ( 0 ) = P1._e()->evalpdflog ( get_vec(Th, "1 2") );
-                LLs ( 1 ) = P2._e()->evalpdflog ( get_vec(Th, "3 2")  );
-                LLs ( 2 ) = PG._e()->evalpdflog (Th );
-                LLs ( 3 ) = M._Mix().logpred ( concat(Thj, vec_1(1.0)) );
+//              LLs ( 0 ) = P1._e()->evalpdflog ( get_vec(Th, "1 2") );
+//              LLs ( 1 ) = P2._e()->evalpdflog ( get_vec(Th, "3 2")  );
+                LLs ( 0 ) = PG._e()->evalpdflog ( Th );
+                LLs ( 1 ) = M._Mix().logpred ( concat ( Thj, vec_1 ( 1.0 ) ) );
+//              LLs ( 3 ) = M2._Mix().logpred ( concat(Thj, vec_1(1.0)) );
                 L.logit ( Li_LL, LLs ); //log-normal
 
                 //Logger
                 L.logit ( Li_Data, vec_3 ( Yt ( t ), rgrg ( 0 ), rgrg ( 1 ) ) );
+                L.logit ( Li_P1m, P1._e()->mean() );
+                L.logit ( Li_P2m, P2._e()->mean() );
                 L.logit ( Li_Gm, PG._e()->mean() );
-                emix *tm =M._Mix().predictor(allj);
-                L.logit ( Li_Mm, tm->mean() );
-                delete tm;
+                L.logit ( Li_Mm, M.mean() );
+
                 L.step ( );
+
+                cout << "Vg: " << PG._e()->_V().to_mat() <<endl;
+                vec mea = M.mean();
+                cout << "Ve: " << M.variance() <<endl;
         }
         L.finalize( );
         L.itsave ( "merg_egiw.it" );
+        cout << endl;
 }

tests/CMakeLists.txt

@@ -42 +42 @@
 TEST(arx_elem_test)
 TEST(merger_test)
+TEST(merger_2d_test)
 TEST(merger_iter_test)
 TEST(mixef_test)