Changeset 629

Timestamp:

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

Author:

smidl

Message:

egiw.variance works for multidimensional + cleanup in tests

Location:

library

Files:

• bdm/stat/exp_family.cpp (modified) (3 diffs)
• bdm/stat/exp_family.h (modified) (1 diff)
• tests/CMakeLists.txt (modified) (1 diff)
• tests/arx_straux_test.cpp (modified) (1 diff)
• tests/egiw.cfg (modified) (4 diffs)
• tests/egiw_harness.cpp (deleted)
• tests/egiw_harness.h (deleted)
• tests/egiw_test.cpp (modified) (1 diff)

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

library/bdm/stat/exp_family.h

@@ -217 +217 @@
                 egiw (int dimx0, ldmat V0, double nu0 = -1.0) : eEF() {set_parameters (dimx0, V0, nu0);};
 
-                void set_parameters (int dimx0, ldmat V0, double nu0 = -1.0) {
-                        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;
-                        }
-                }
+                void set_parameters (int dimx0, ldmat V0, double nu0 = -1.0);
                 //!@}
 

library/tests/CMakeLists.txt

@@ -7 +7 @@
 link_directories (./unittest-cpp)
 
-SET(testutil_src egiw_harness.cpp egiw_harness.h epdf_harness.cpp epdf_harness.h mat_checks.cpp mat_checks.h 
+SET(testutil_src epdf_harness.cpp epdf_harness.h mat_checks.cpp mat_checks.h 
         mpdf_harness.cpp mpdf_harness.h square_mat_point.cpp square_mat_point.h test_util.cpp test_util.h)
 

library/tests/arx_straux_test.cpp

@@ -5 +5 @@
 
 TEST(test_arx_straux){
-mat A="[   0.8147    0.9134    0.2785;"
+mat A="0.8147    0.9134    0.2785;"
     "0.9058    0.6324    0.5469;"
-    "0.1270    0.0975    0.9575]";
+    "0.1270    0.0975    0.9575";
 
-mat B="[0.9649    0.9572    0.1419;"
+mat B="0.9649    0.9572    0.1419;"
     "0.1576    0.4854    0.4218;"
-    "0.9706    0.8003    0.9157]";
+    "0.9706    0.8003    0.9157";
 //when updateing matrices do not forget to update CHECK_EQUAL below!!!
         

library/tests/egiw.cfg

@@ -1 +1 @@
 data = (
 {
-  class = "egiw_harness";
+  class = "epdf_harness";
   epdf = {
     class = "egiw";
@@ -13 +13 @@
   };
   mean = [ 1.1, 0.1 ];
-  lognc = 3.39463;
   variance = [ 0.01, 8e-05 ];
   support = ( "matrix", 2, 2, [ -2.0, 4.0, 0.01, 2.0 ] );
@@ -20 +19 @@
 },
 {
-  class = "egiw_harness";
+  class = "epdf_harness";
   epdf = {
     class = "egiw";
-    V = ( "matrix", 2, 2, [ 20.0, 8.0, 8.0, 4.0 ] );
-    nu = 4.0;
+    V = ( "matrix", 2, 2, [ 200.0, 80.0, 80.0, 40.0 ] );
+    nu = 40.0;
     dimx = 1;
     rv : 
@@ -32 +31 @@
     };
   };
-  mean = [ 2, -4 ];
-  lognc = 0.451583;
-  variance = [ -1, -16 ];
+  mean = [2.0, 1.14286];
+  variance = [0.0285714, 0.0395795];
 } );
 

library/tests/egiw_test.cpp

@@ -5 +5 @@
 #include "itpp_ext.h"
 #include "epdf_harness.h"
-#include "egiw_harness.h"
 #include "mat_checks.h"
 #include "UnitTest++.h"