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

Timestamp:

09/18/09 00:17:11 (15 years ago)

Author:

smidl

Message:

egiw.variance works for multidimensional + cleanup in tests

Files:

• library/bdm/stat/exp_family.cpp (modified) (3 diffs)

library/bdm/stat/exp_family.cpp

@@ -17 +17 @@
         this->bayes ( dt, 1.0 );
 };
+
+void egiw::set_parameters (int dimx0, ldmat V0, double nu0) {
+        dimx = dimx0;
+        nPsi = V0.rows() - dimx;
+        dim = dimx * (dimx + nPsi); // size(R) + size(Theta)
+        
+        V = V0;
+        if (nu0 < 0) {
+                nu = 0.1 + nPsi + 2 * dimx + 2; // +2 assures finite expected value of R
+                // terms before that are sufficient for finite normalization
+        } else {
+                nu = nu0;
+        }
+}
 
 vec egiw::sample() const {
@@ -119 +133 @@
                 mat R;
                 mean_mat ( M, R );
-                return cvectorize ( concat_vertical ( M, R ) );
+                return concat( cvectorize (  M),cvectorize( R ) );
         }
 
@@ -125 +139 @@
 
 vec egiw::variance() const {
-
+        int l = V.rows();
+        // cut out rest of lower-right part of V
+        const ldmat tmp ( V, linspace ( dimx, l - 1 ) );
+        // invert it
+        ldmat itmp ( l );
+        tmp.inv ( itmp );
+        
+        // following Wikipedia notation
+        // m=nu-nPsi-dimx-1, p=dimx 
+        double mp1p=nu-nPsi-2*dimx; // m-p+1
+        double mp1m=mp1p-2;         // m-p-1
+        
         if ( dimx == 1 ) {
-                int l = V.rows();
-                const ldmat tmp ( V, linspace ( 1, l - 1 ) );
-                ldmat itmp ( l );
-                tmp.inv ( itmp );
-                double cove = V._D() ( 0 ) / ( nu - nPsi - 2 * dimx - 2 );
+                double cove = V._D() ( 0 ) / mp1m ;
 
                 vec var ( l );
                 var.set_subvector ( 0, diag ( itmp.to_mat() ) *cove );
-                var ( l - 1 ) = cove * cove / ( nu - nPsi - 2 * dimx - 2 );
+                var ( l - 1 ) = cove * cove / (mp1m-2);
                 return var;
         } else {
-                bdm_error ( "not implemented" );
-                return vec();
+                ldmat Vll( V, linspace(0,dimx-1)); // top-left part of V
+                mat Y=Vll.to_mat();
+                mat varY(Y.rows(), Y.cols()); 
+                
+                double denom = (mp1p-1)*mp1m*mp1m*(mp1m-2);         // (m-p)(m-p-1)^2(m-p-3)
+                
+                int i,j;
+                for ( i=0; i<Y.rows(); i++){
+                        for ( j=0; j<Y.cols(); j++){
+                                varY(i,j) = (mp1p*Y(i,j)*Y(i,j) + mp1m * Y(i,i)* Y(j,j)) /denom;
+                        }
+                }
+                vec mean_dR = diag(Y)/mp1m; // corresponds to cove
+                vec var_th=diag ( itmp.to_mat() );
+                vec var_Th ( mean_dR.length()*var_th.length() );
+                // diagonal of diag(mean_dR) \kron diag(var_th)
+                for (int i=0; i<mean_dR.length(); i++){
+                        var_Th.set_subvector(i*var_th.length(), var_th*mean_dR(i));  
+                }
+                
+                return concat(var_Th, cvectorize(varY));
         }
 }