Changeset 725

Timestamp:

11/17/09 00:53:54 (15 years ago)

Author:

smidl

Message:

Sampling from egiw

Location:

library

Files:

• bdm/math/chmat.h (modified) (1 diff)
• bdm/math/square_mat.cpp (modified) (1 diff)
• bdm/stat/exp_family.cpp (modified) (3 diffs)
• bdm/stat/exp_family.h (modified) (5 diffs)
• tests/epdf_harness.cpp (modified) (4 diffs)
• tests/testsuite/egiw.cfg (modified) (2 diffs)
• tests/testsuite/epdf_test.cpp (modified) (2 diffs)

library/bdm/math/chmat.h

@@ -57 +57 @@
         void inv ( chmat &Inv ) const   {
                 ( Inv.Ch ) = itpp::inv ( Ch ).T();
+                Inv.dim = dim;
         }; //Fixme: can be more efficient
         ;

library/bdm/math/square_mat.cpp

@@ -115 +115 @@
                 }
         }
-        mat V2 = L.transpose() * diag ( D ) * L;
-        return V2;
+        //mat V2 = L.transpose() * diag ( D ) * L;
+        return V;
 }
 

library/bdm/stat/exp_family.cpp

@@ -33 +33 @@
 
 vec egiw::sample() const {
-        bdm_warning ( "Function not implemented" );
-        return vec_1 ( 0.0 );
+        mat M;
+        chmat R;
+        sample_mat(M,R);
+        
+        return concat (cvectorize(M),cvectorize(R.to_mat()));
+}
+
+void egiw::sample_mat(mat &Mi, chmat &Ri)const{
+        
+        // TODO - correct approach - convert to product of norm * Wishart
+        mat M;
+        ldmat Vz;
+        ldmat Lam;
+        factorize(M,Vz,Lam);
+        
+        chmat Ch;
+        Ch.setCh(Lam._L()*diag(sqrt(Lam._D())));
+        chmat iCh;
+        Ch.inv(iCh);
+        
+        eWishartCh Omega; //inverse Wishart, result is R, 
+        Omega.set_parameters(iCh,nu-2*nPsi-dimx); // 2*nPsi is there to match numercial simulations - check if analytically correct
+        
+        chmat Omi;
+        Omi.setCh(Omega.sample_mat());
+        
+        mat Z=randn(M.rows(), M.cols());
+        Mi = M+ Omi._Ch() * Z * inv(Vz._L()*diag(sqrt(Vz._D())));
+        Omi.inv(Ri);
 }
 
@@ -104 +131 @@
 }
 
+void egiw::factorize(mat &M, ldmat &Vz, ldmat &Lam) const{      
+        const mat &L = V._L();
+        const vec &D = V._D();
+        int end = L.rows() - 1;
+        
+        Vz=ldmat(L ( dimx, end, dimx, end ), D(dimx,end));
+        mat iLsub = ltuinv ( Vz._L());
+        // set mean value
+        mat Lpsi = L ( dimx, end, 0, dimx - 1 );
+        M = iLsub * Lpsi;
+        
+        Lam =ldmat ( L (0, dimx-1, 0, dimx-1 ), D (0, dimx-1 ) );  //exp val of R
+         if (1){ // test with Peterka
+                 mat VF=V.to_mat();
+                 mat Vf=VF(0,dimx-1,0, dimx-1);
+                 mat Vzf = VF(dimx,end,0,dimx-1);
+                 mat VZ = VF(dimx,end,dimx,end);
+                 
+                 mat Lam2 = Vf-Vzf.T()*inv(VZ)*Vzf;
+         }
+}
+
 ldmat egiw::est_theta_cov() const {
         if ( dimx == 1 ) {
@@ -111 +160 @@
 
                 mat Lsub = L ( 1, end, 1, end );
-                mat Dsub = diag ( D ( 1, end ) );
-
-                return inv ( transpose ( Lsub ) * Dsub * Lsub );
+//              mat Dsub = diag ( D ( 1, end ) );
+
+                ldmat LD(inv(Lsub).T(), 1.0/D(1,end));
+                return LD;
 
         } else {

library/bdm/stat/exp_family.h

@@ -223 +223 @@
                 vec mean() const;
                 vec variance() const;
-
+                void sample_mat(mat &Mi, chmat &Ri)const;
+                        
+                void factorize(mat &M, ldmat &Vz, ldmat &Lam) const;
                 //! LS estimate of \f$\theta\f$
                 vec est_theta() const;
@@ -245 +247 @@
                 const double& _nu() const {return nu;}
                 const int & _dimx() const {return dimx;}
+                        
                 /*! Create Gauss-inverse-Wishart density 
                 \f[ f(rv) = GiW(V,\nu) \f]
@@ -259 +262 @@
                 \endcode
                 */
-                
+                                
                 void from_setting (const Setting &set) {
                         epdf::from_setting(set);
@@ -274 +277 @@
                         }
                 }
+                
+                void to_setting ( Setting& set ) const{
+                        epdf::to_setting(set);
+                        UI::save(dimx,set,"dimx");
+                        UI::save(V.to_mat(),set,"V");
+                        UI::save(nu,set,"nu");                  
+                };
+                
                 void validate(){
                         // check sizes, rvs etc.
+                }
+                void log_register( bdm::logger& L, const string& prefix ){
+                        if (log_level==3){
+                                root::log_register(L,prefix);
+                                logrec->ids.set_length(2);
+                                int th_dim=dimension()-dimx*(dimx+1)/2;
+                                logrec->ids(0)=L.add(RV("",th_dim), prefix + logrec->L.prefix_sep() +"mean");
+                                logrec->ids(1)=L.add(RV("",th_dim*th_dim),prefix + logrec->L.prefix_sep() + "variance"); 
+                        } else {
+                                epdf::log_register(L,prefix);
+                        }
+                }
+                void log_write() const {
+                        if (log_level==3){
+                                mat M;
+                                ldmat Lam;
+                                ldmat Vz;
+                                factorize(M,Vz,Lam);
+                                logrec->L.logit(logrec->ids(0), est_theta() );
+                                logrec->L.logit(logrec->ids(1), cvectorize(est_theta_cov().to_mat()));
+                        } else {
+                                epdf::log_write();
+                        }
+                        
                 }
                 //!@}
@@ -1138 +1173 @@
                 //! Set internal structures
                 void set_parameters (const mat &Y0, const double delta0) {Y = chmat (Y0);delta = delta0; p = Y.rows(); dim = p * p; }
+                //! Set internal structures
+                void set_parameters (const chmat &Y0, const double delta0) {Y = Y0;delta = delta0; p = Y.rows(); dim = p * p; }
                 //! Sample matrix argument
                 mat sample_mat() const {

library/tests/epdf_harness.cpp

@@ -59 +59 @@
         }
 
-        if ( R.rows() > 0 ) {
+        if ( variance.length() > 0 ) {
                 check_sample_mean();
+        }
+        if (R.rows()>0){
                 check_covariance();
         }
@@ -114 +116 @@
         vec yb = support.get_row ( 1 );
 
-        int tc = 0;
-        Array<vec> actual ( CurrentContext::max_trial_count );
-        do {
-                vec emu = num_mean2 ( hepdf.get(), xb, yb, nbins ( 0 ), nbins ( 1 ) );
-                actual ( tc ) = emu;
-                ++tc;
-        } while ( ( tc < CurrentContext::max_trial_count ) &&
-                  !UnitTest::AreClose ( mean, actual ( tc - 1 ), tolerance ) );
-        if ( ( tc == CurrentContext::max_trial_count ) &&
-                ( !UnitTest::AreClose ( mean, actual ( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
-                UnitTest::MemoryOutStream stream;
-                stream << CurrentContext::format_context ( __LINE__ ) << "expected " << mean << " +/- " << tolerance << " but was " << actual;
-
-                UnitTest::TestDetails details ( *UnitTest::CurrentTest::Details(), 0, false );
-
-                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
-        }
+        vec actual;
+        actual = num_mean2 ( hepdf.get(), xb, yb, nbins ( 0 ), nbins ( 1 ) );
+
+        CHECK_CLOSE(mean, actual, tolerance);
 }
 
@@ -137 +126 @@
         vec yb = support.get_row ( 1 );
 
-        int tc = 0;
-        Array<double> actual ( CurrentContext::max_trial_count );
-        do {
-                double nc = normcoef ( hepdf.get(), xb, yb, nbins ( 0 ), nbins ( 1 ) );
-                actual ( tc ) = nc;
-                ++tc;
-        } while ( ( tc < CurrentContext::max_trial_count ) &&
-                  !UnitTest::AreClose ( 1.0, actual ( tc - 1 ), tolerance ) );
-        if ( ( tc == CurrentContext::max_trial_count ) &&
-                ( !UnitTest::AreClose ( 1.0, actual ( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
-                UnitTest::MemoryOutStream stream;
-                stream << CurrentContext::format_context ( __LINE__ ) << "expected " << mean << " +/- " << tolerance << " but was " << actual;
-
-                UnitTest::TestDetails details ( *UnitTest::CurrentTest::Details(), 0, false );
-
-                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
-        }
+        double nc = normcoef ( hepdf.get(), xb, yb, nbins ( 0 ), nbins ( 1 ) );
+        CHECK_CLOSE(1.0,nc,0.01);
 }
 
@@ -168 +142 @@
         } while ( ( tc < CurrentContext::max_trial_count ) &&
                   !UnitTest::AreClose ( mean, actual ( tc - 1 ), delta ) );
+                          
         if ( ( tc == CurrentContext::max_trial_count ) &&
                 ( !UnitTest::AreClose ( mean, actual ( CurrentContext::max_trial_count - 1 ), delta ) ) ) {

library/tests/testsuite/egiw.cfg

@@ -12 +12 @@
     };
   };
-  mean = [ 1.1, 0.1 ];
-  variance = [ 0.01, 8e-05 ];
+  mean = [ 1.1, 0.2 ];
+  variance = [ 0.02, 0.01333 ];
   support = ( "matrix", 2, 2, [ -2.0, 4.0, 0.01, 2.0 ] );
   nbins = [ 100, 200 ];
-  tolerance = 0.2;
+  tolerance = 0.01;
 },
 {
@@ -33 +33 @@
   mean = [2.0, 1.14286];
   variance = [0.0285714, 0.0395795];
+  support = ( "matrix", 2, 2, [ 0.0, 4.0, 0.01, 3.0 ] );
+  nbins = [ 100, 200 ];
+  tolerance = 0.01;
 } );
 

library/tests/testsuite/epdf_test.cpp

@@ -43 +43 @@
 
 TEST ( ewishart_test ) {
-        mat wM = "1.0 0.9; 0.9 1.0";
+        mat wM = "10.0 0.9; 0.9 1.0";
         eWishartCh eW;
         eW.set_parameters ( wM / 100, 100 );
@@ -53 +53 @@
         }
 
-        mat observed ( "0.978486 0.88637; 0.88637 0.992141" );
         mat actual = mea / 100;
-        CHECK_CLOSE ( observed, actual, 0.1 );
+        CHECK_CLOSE ( wM, actual, 0.1 );
 }