Changeset 310

Timestamp:

04/09/09 14:33:08 (15 years ago)

Author:

smidl

Message:

Published example 2d

Files:

• bdm/estim/merger.cpp (modified) (4 diffs)
• bdm/estim/merger.h (modified) (2 diffs)
• bdm/stat/emix.h (modified) (1 diff)
• bdm/stat/libEF.h (modified) (2 diffs)
• mpdm/TR2244/merger_iter_cond.cpp (modified) (4 diffs)
• mpdm/TR2244/merger_iter_debug.m (added)

bdm/estim/merger.cpp

@@ -88 +88 @@
 
                         // This will be active only later in iterations!!!
-                        if ( ( !fix_smp ) & ( 1./sum_sqr ( w ) <0.9*Ns ) )
+                        if ( ( !fix_smp ) & ( 1./sum_sqr ( w ) <effss_coef*Ns ) )
                         {
                                 // Generate new samples
@@ -98 +98 @@
                                         set_col_part ( Smp_ex,i,Smp ( i ) );
                                         //Importance of the mixture
-                                        lw_mix ( i ) =Mix.logpred (Smp_ex.get_col(i) );
+                                        //lw_mix ( i ) =Mix.logpred (Smp_ex.get_col(i) );
+                                        lw_mix ( i ) = Mpred->evallog(Smp(i));
                                 }
                                 if ( 1 )
                                 {
                                         cout<<"Resampling =" << 1./sum_sqr ( w ) << endl;
+                                        cout << Mix._e()->mean() <<endl;
                                         cout << sum ( Smp_ex,2 ) /Ns <<endl;
                                         cout << Smp_ex*Smp_ex.T() /Ns << endl;
@@ -111 +113 @@
                                 sprintf ( str,"Mpred_mean%d",niter );
                                 *dbg << Name ( str ) << Mpred->mean();
+                                sprintf ( str,"Mpred_var%d",niter );
+                                *dbg << Name ( str ) << Mpred->variance();
 
 
@@ -216 +220 @@
                         // ==== stopping rule ===
                         niter++;
-                        converged = ( niter>20 );
+                        converged = ( niter>40 );
                 }
                 delete Mpred;

bdm/estim/merger.h

@@ -47 +47 @@
                         //! Projection to empirical density
                         eEmp eSmp;
+                        //! coefficient of resampling
+                        double effss_coef;
 
                         //! debug or not debug
@@ -88 +90 @@
                         void debug_file ( const string fname ) { if ( DBG ) delete dbg; dbg = new it_file ( fname ); if ( dbg ) DBG=true;}
 //! Set internal parameters used in approximation
-                        void set_parameters ( double beta0, int Ns0, int Nc0 ) {beta=beta0;Ns=Ns0;Nc=Nc0;eSmp.set_parameters ( Ns0,false );}
+                        void set_parameters ( double beta0, int Ns0, int Nc0, double effss_coef0=0.5 ) {beta=beta0;
+                                Ns=Ns0;
+                                Nc=Nc0;
+                                effss_coef=effss_coef0;
+                                eSmp.set_parameters ( Ns0,false );
+                        }
                         void set_grid ( Array<vec> &XYZ )
                         {

bdm/stat/emix.h

@@ -112 +112 @@
                 //non-central moment
                 vec mom2 = zeros ( dim );
-                for ( int i = 0;i < w.length();i++ ) {mom2 += w ( i ) * pow ( Coms ( i )->mean(),2 ); }
+                for ( int i = 0;i < w.length();i++ ) {mom2 += w ( i ) * (Coms(i)->variance() + pow ( Coms ( i )->mean(),2 )); }
                 //central moment
                 return mom2-pow ( mean(),2 );

bdm/stat/libEF.h

@@ -230 +230 @@
                         //!sufficient statistics
                         vec beta;
-                        //!speedup variable
-                        double gamma;
                 public:
                         //!\name Constructors
@@ -243 +241 @@
                                 beta= beta0;
                                 dim = beta.length();
-                                gamma = sum ( beta );
                         }
                         //!@}
 
                         vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
-                        vec mean() const {return beta/gamma;};
-                        vec variance() const {return elem_mult ( beta, ( beta+1 ) ) / ( gamma* ( gamma+1 ) );}
+                        vec mean() const {return beta/sum(beta);};
+                        vec variance() const {double gamma =sum(beta); return elem_mult ( beta, ( beta+1 ) ) / ( gamma* ( gamma+1 ) );}
                         //! In this instance, val is ...
                         double evallog_nn ( const vec &val ) const

mpdm/TR2244/merger_iter_cond.cpp

@@ -25 +25 @@
 
         //Differneces in constant term are essential
-        f1.set_parameters ( "0.4","1",mat ( "0.04" ) );
-        f2.set_parameters ( "0.2", mat ( "0.08" ) );
+        f1.set_parameters ( "1","1",mat ( "0.01" ) );
+        f2.set_parameters ( "2", mat ( "0.02" ) );
 
         mepdf mf2(&f2);
@@ -40 +40 @@
 /*      Cov(1,2)=0.29;
         Cov(2,1)=0.29;*/
-        g0.set_parameters ( vec ( "0 0" ),Cov);// +1*ones ( 3,3 ) );
+        g0.set_parameters ( vec ( "2 1" ),Cov);// +1*ones ( 3,3 ) );
 
 //      enorm<ldmat>* teste=g0.marginal(concat(u1,u2));
 //      mlnorm<ldmat>* testm=(mlnorm<ldmat>*)teste->condition(u2);
         
-        M.set_parameters ( 1.2,1000,1);
+        M.set_parameters ( 1.2,1000,1,0.5);
         
         Array<vec> YUU(2); 
-        YUU(0)=linspace(-2.9,3.1,20);
-        YUU(1)=linspace(-5,5,20);
-        //YUU(2)=linspace(-3,3,20);
+        YUU(0)=linspace(1.5,2.5,20);
+        YUU(1)=linspace(0.5,1.5,20);
+        
+        YUU(0)=linspace(0,4,20);
+        YUU(1)=linspace(-1,3,20);
+//YUU(2)=linspace(-3,3,20);
 
-        M.set_grid(YUU);
+        //M.set_grid(YUU);
         
         int Ntrials=1;
@@ -68 +71 @@
                 epdf* MP = MM._Coms ( 0 )->epredictor ( );
                 MP->set_rv(all);
+                
+                //Missing piece
+                enorm<ldmat>* miss=     (enorm<ldmat>*) MP->marginal(u1);
+                cout << "Missing: " << miss->_mu() << " +- " << miss->_R().to_mat() << endl;
 
                 RV yu1 = y; yu1.add ( u1 );
@@ -102 +109 @@
         cout << "R1: " << R1mean << " +- " << 2*sqrt(sum_sqr(R1s)/Ntrials-R1mean*R1mean) <<endl;
         cout << "R2: " << R2mean << " +- " << 2*sqrt(sum_sqr(R2s)/Ntrials-R2mean*R2mean) <<endl;
+        
+        
 }