Changeset 810

Timestamp:

02/21/10 20:58:41 (15 years ago)

Author:

smidl

Message:

egiw with empty normal part

Location:

library

Files:

• bdm/estim/arx.h (modified) (1 diff)
• bdm/stat/exp_family.cpp (modified) (7 diffs)
• tests/testsuite/egiw.cfg (modified) (1 diff)
• tests/testsuite/enorm.cfg (modified) (2 diffs)

library/bdm/estim/arx.h

@@ -152 +152 @@
                 vec dV0 ( prior._V().rows() );
                 dV0.set_subvector ( 0, prior._dimx() - 1, 1.0 );
-                dV0.set_subvector ( prior._dimx(), dV0.length() - 1, 1e-5 );
+                if (dV0.length()>prior._dimx())
+                        dV0.set_subvector ( prior._dimx(), dV0.length() - 1, 1e-5 );
 
                 prior.set_parameters ( prior._dimx(), ldmat ( dV0 ) );

library/bdm/stat/exp_family.cpp

@@ -90 +90 @@
         Omi.setCh ( Omega.sample_mat() );
 
-        mat Z = randn ( M.rows(), M.cols() );
-        Mi = M + Omi._Ch() * Z * inv ( Vz._L() * diag ( sqrt ( Vz._D() ) ) );
+        if (M._datasize()>0){
+                mat Z = randn ( M.rows(), M.cols() );
+                Mi = M + Omi._Ch() * Z * inv ( Vz._L() * diag ( sqrt ( Vz._D() ) ) );
+        }
         Omi.inv ( Ri );
 }
 
 double egiw::evallog_nn ( const vec &val ) const {
-        bdm_assert_debug(val.length()==nPsi+dimx,"Incorrect cond in egiw::evallog_nn" );
+        bdm_assert_debug(val.length()==dimx*(nPsi+dimx),"Incorrect cond in egiw::evallog_nn" );
         
         int vend = val.length() - 1;
@@ -110 +112 @@
                 return -0.5* ( nu*log ( r ) + Vq / r );
         } else {
-                mat Th = reshape ( val ( 0, nPsi * dimx - 1 ), nPsi, dimx );
+                mat Tmp;
+                if (nPsi>0){
+                        mat Th = reshape ( val ( 0, nPsi * dimx - 1 ), nPsi, dimx );
+                        Tmp = concat_vertical ( -eye ( dimx ), Th );
+                } else {
+                        Tmp = -eye(dimx);
+                }
                 fsqmat R ( reshape ( val ( nPsi*dimx, vend ), dimx, dimx ) );
                 double ldetR = R.logdet();
-                if ( ldetR ) return -inf;
-                mat Tmp = concat_vertical ( -eye ( dimx ), Th );
+                if ( !isfinite(ldetR) ) return -inf;
                 fsqmat iR ( dimx );
                 R.inv ( iR );
@@ -131 +138 @@
 #define Inf std::numeric_limits<double>::infinity()
 
-        double nkG = 0.5 * dimx * ( -nPsi * log2pi + sum ( log ( D ( dimx, D.length() - 1 ) ) ) );
+        double nkG = 0.5 * dimx * ( -nPsi * log2pi + sum ( log ( D.mid ( dimx, nPsi ) ) ) );
         // temporary for lgamma in Wishart
         double lg = 0;
@@ -169 +176 @@
         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
+
+        if (dimx<=end){
+                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;
+        }
+/*      if ( 0 ) { // test with Peterka
                 mat VF = V.to_mat();
                 mat Vf = VF ( 0, dimx - 1, 0, dimx - 1 );
@@ -183 +192 @@
 
                 mat Lam2 = Vf - Vzf.T() * inv ( VZ ) * Vzf;
-        }
+        }*/
 }
 
@@ -227 +236 @@
         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 );
-
+        if (dimx<l){
+                const ldmat tmp ( V, linspace ( dimx, l - 1 ) );
+                tmp.inv ( itmp );
+        }
         // following Wikipedia notation
         // m=nu-nPsi-dimx-1, p=dimx
@@ -275 +285 @@
 
         ldmat ldR ( L ( 0, dimx - 1, 0, dimx - 1 ), D ( 0, dimx - 1 ) / ( nu - nPsi - 2*dimx - 2 ) ); //exp val of R
-        mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) );
 
         // set mean value
-        mat Lpsi = L ( dimx, end, 0, dimx - 1 );
-        M = iLsub * Lpsi;
+        if (dimx<=end){
+                mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) );
+                mat Lpsi = L ( dimx, end, 0, dimx - 1 );
+                M = iLsub * Lpsi;
+        }
         R = ldR.to_mat()  ;
 }

library/tests/testsuite/egiw.cfg

@@ -45 +45 @@
         };
   tolerance = 0.01;
-} );
+},
+{
+  class = "epdf_harness";
+  epdf = {
+    class = "egiw";
+    V = ( "matrix", 2, 2, [ 10.0, 5.0, 5.0, 15.0 ] );
+    nu = 10.0;
+    dimx = 2;
+    rv : 
+    {
+      class = "RV";
+      names = ( "a", "b" );
+    };
+  };
+  mean = [1.0, 0.5, 0.5, 1.5];
+  variance = [1.0, 1.0];
+  R = ( "matrix", 2, 2, [ 1.0, 0.0, 0.0, 1.0 ] );
+  support = {
+        class = "rectangular_support";
+        ranges = ( [ 0.0, 4.0 ], [-2.0, 2.0], [-2.0, 2.0], [0.0, 4.0 ] );
+        gridsizes = [ 10, 10, 10, 10];
+        };
+  tolerance = 0.01;
+});
 
 

library/tests/testsuite/enorm.cfg

@@ -103 +103 @@
     class = "enorm<chmat>";
     mu = ( "matrix", 1, 2, [ 0.0, 2.0 ] );
-    R = ( "matrix", 2, 2, [ 2.0, 0.0, 0.0, 0.5 ] );
+    R = ( "matrix", 2, 2, [ 1.0, 0.3, 0.3, 0.5 ] );
     rv : 
     {
@@ -111 +111 @@
   };
   mean = [ 0.0, 2.0 ];
-  variance = [ 2.0, 0.5 ];
-  R = ( "matrix", 2, 2, [ 2.0, 0.0, 0.0, 0.5 ] );
-  tolerance = 0.5;
+  variance = [ 1.0, 0.5 ];
+  R = ( "matrix", 2, 2, [ 1.0, 0.3, 0.3, 0.5 ] );
 } );