Changeset 725 for library/bdm/stat

Timestamp:

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

Author:

smidl

Message:

Sampling from egiw

Location:

library/bdm/stat

Files:

• exp_family.cpp (modified) (3 diffs)
• exp_family.h (modified) (5 diffs)

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 {