Changeset 810 for library/bdm/stat

Timestamp:

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

Author:

smidl

Message:

egiw with empty normal part

Files:

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

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