Changeset 810 for library/bdm/stat/exp_family.cpp
- Timestamp:
- 02/21/10 20:58:41 (14 years ago)
- Files:
-
- 1 modified
Legend:
- Unmodified
- Added
- Removed
-
library/bdm/stat/exp_family.cpp
r761 r810 90 90 Omi.setCh ( Omega.sample_mat() ); 91 91 92 mat Z = randn ( M.rows(), M.cols() ); 93 Mi = M + Omi._Ch() * Z * inv ( Vz._L() * diag ( sqrt ( Vz._D() ) ) ); 92 if (M._datasize()>0){ 93 mat Z = randn ( M.rows(), M.cols() ); 94 Mi = M + Omi._Ch() * Z * inv ( Vz._L() * diag ( sqrt ( Vz._D() ) ) ); 95 } 94 96 Omi.inv ( Ri ); 95 97 } 96 98 97 99 double egiw::evallog_nn ( const vec &val ) const { 98 bdm_assert_debug(val.length()== nPsi+dimx,"Incorrect cond in egiw::evallog_nn" );100 bdm_assert_debug(val.length()==dimx*(nPsi+dimx),"Incorrect cond in egiw::evallog_nn" ); 99 101 100 102 int vend = val.length() - 1; … … 110 112 return -0.5* ( nu*log ( r ) + Vq / r ); 111 113 } else { 112 mat Th = reshape ( val ( 0, nPsi * dimx - 1 ), nPsi, dimx ); 114 mat Tmp; 115 if (nPsi>0){ 116 mat Th = reshape ( val ( 0, nPsi * dimx - 1 ), nPsi, dimx ); 117 Tmp = concat_vertical ( -eye ( dimx ), Th ); 118 } else { 119 Tmp = -eye(dimx); 120 } 113 121 fsqmat R ( reshape ( val ( nPsi*dimx, vend ), dimx, dimx ) ); 114 122 double ldetR = R.logdet(); 115 if ( ldetR ) return -inf; 116 mat Tmp = concat_vertical ( -eye ( dimx ), Th ); 123 if ( !isfinite(ldetR) ) return -inf; 117 124 fsqmat iR ( dimx ); 118 125 R.inv ( iR ); … … 131 138 #define Inf std::numeric_limits<double>::infinity() 132 139 133 double nkG = 0.5 * dimx * ( -nPsi * log2pi + sum ( log ( D ( dimx, D.length() - 1) ) ) );140 double nkG = 0.5 * dimx * ( -nPsi * log2pi + sum ( log ( D.mid ( dimx, nPsi ) ) ) ); 134 141 // temporary for lgamma in Wishart 135 142 double lg = 0; … … 169 176 int end = L.rows() - 1; 170 177 171 Vz = ldmat ( L ( dimx, end, dimx, end ), D ( dimx, end ) );172 mat iLsub = ltuinv ( Vz._L() );173 // set mean value174 mat Lpsi = L ( dimx, end, 0, dimx - 1 );175 M = iLsub * Lpsi;176 177 178 Lam = ldmat ( L ( 0, dimx - 1, 0, dimx - 1 ), D ( 0, dimx - 1 ) ); //exp val of R 178 if ( 1 ) { // test with Peterka 179 180 if (dimx<=end){ 181 Vz = ldmat ( L ( dimx, end, dimx, end ), D ( dimx, end ) ); 182 mat iLsub = ltuinv ( Vz._L() ); 183 // set mean value 184 mat Lpsi = L ( dimx, end, 0, dimx - 1 ); 185 M = iLsub * Lpsi; 186 } 187 /* if ( 0 ) { // test with Peterka 179 188 mat VF = V.to_mat(); 180 189 mat Vf = VF ( 0, dimx - 1, 0, dimx - 1 ); … … 183 192 184 193 mat Lam2 = Vf - Vzf.T() * inv ( VZ ) * Vzf; 185 } 194 }*/ 186 195 } 187 196 … … 227 236 int l = V.rows(); 228 237 // cut out rest of lower-right part of V 229 const ldmat tmp ( V, linspace ( dimx, l - 1 ) );230 238 // invert it 231 239 ldmat itmp ( l ); 232 tmp.inv ( itmp ); 233 240 if (dimx<l){ 241 const ldmat tmp ( V, linspace ( dimx, l - 1 ) ); 242 tmp.inv ( itmp ); 243 } 234 244 // following Wikipedia notation 235 245 // m=nu-nPsi-dimx-1, p=dimx … … 275 285 276 286 ldmat ldR ( L ( 0, dimx - 1, 0, dimx - 1 ), D ( 0, dimx - 1 ) / ( nu - nPsi - 2*dimx - 2 ) ); //exp val of R 277 mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) );278 287 279 288 // set mean value 280 mat Lpsi = L ( dimx, end, 0, dimx - 1 ); 281 M = iLsub * Lpsi; 289 if (dimx<=end){ 290 mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) ); 291 mat Lpsi = L ( dimx, end, 0, dimx - 1 ); 292 M = iLsub * Lpsi; 293 } 282 294 R = ldR.to_mat() ; 283 295 }