Changeset 730 for library/bdm/stat

Timestamp:

11/17/09 20:55:03 (15 years ago)

Author:

smidl

Message:

Reworked epdf_harness - testing is now more general.

Corrections of existing tests.

Location:

library/bdm/stat

Files:

• discrete.cpp (modified) (1 diff)
• discrete.h (modified) (3 diffs)
• exp_family.cpp (modified) (1 diff)
• exp_family.h (modified) (1 diff)

library/bdm/stat/discrete.cpp

@@ -122 +122 @@
 }
 
+void grid_fnc::set_values (const epdf &ep ) {
+        if ( sup.points() > 0 ) {
+                values ( 0 ) =  exp(ep.evallog( sup.first_vec() ));
+                for ( int j = 1; j < sup.points(); j++ ) {
+                        values ( j ) = exp(ep.evallog( sup.next_vec() ));
+                }
+        }
 }
+
+}

library/bdm/stat/discrete.h

@@ -62 +62 @@
         //! Get next active vector, call ONLY after first_vector()!
         const vec& next_vec();
-
+        
+        //! return active vector, call ONLY after first_vector()!
+        const vec& act_vec(){return actvec;};
+        
         //! \todo to je asi navic .. v predkovi!
         ivec nearest_point ( const vec &val );
@@ -140 +143 @@
         //! constructor function fills values by calling function \c f , double f(vec&), given by a pointer
         void set_values ( double ( *evalptr ) ( const vec& ) );
-
+        
+        //! constructor function fills values by calling function \c f , double f(vec&), given by a pointer
+        void set_values ( const epdf &ep);
+        
         //! get value nearest to the given point
         double nearest_val ( const vec &val ) {
@@ -149 +155 @@
                 return vec_1 ( nearest_val ( val ) );
         }
+        const vec & _values() const { return values; }
 };
 UIREGISTER ( grid_fnc );

library/bdm/stat/exp_family.cpp

@@ -38 +38 @@
         
         return concat (cvectorize(M),cvectorize(R.to_mat()));
+}
+
+mat egiw::sample_mat(int n) const {
+        // TODO - correct approach - convert to product of norm * Wishart
+        mat M;
+        ldmat Vz;
+        ldmat Lam;
+        factorize(M,Vz,Lam);
+        
+        chmat ChLam(Lam.to_mat());
+        chmat iChLam;
+        ChLam.inv(iChLam);
+        
+        eWishartCh Omega; //inverse Wishart, result is R, 
+        Omega.set_parameters(iChLam,nu-2*nPsi-dimx); // 2*nPsi is there to match numercial simulations - check if analytically correct
+        
+        mat OmChi;
+        mat Z(M.rows(),M.cols());
+
+        mat Mi;
+        mat RChiT;
+        mat tmp(dimension(), n);
+        for (int i=0; i<n;i++){
+                OmChi=Omega.sample_mat();
+                RChiT=inv(OmChi);
+                Z=randn(M.rows(), M.cols());
+                Mi = M + RChiT * Z * inv(Vz._L().T() *diag(sqrt(Vz._D())));
+        
+                tmp.set_col(i,concat (cvectorize(Mi),cvectorize(RChiT*RChiT.T())));
+        }
+        return tmp;
 }
 

library/bdm/stat/exp_family.h

@@ -221 +221 @@
 
                 vec sample() const;
+                mat sample_mat(int n) const;
                 vec mean() const;
                 vec variance() const;