Changeset 730 for library/bdm/stat/exp_family.cpp

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.

Files:

• library/bdm/stat/exp_family.cpp (modified) (1 diff)

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;
 }