[357] | 1 | #include <math.h> |
---|
[262] | 2 | |
---|
[32] | 3 | #include <itpp/base/bessel.h> |
---|
[384] | 4 | #include "exp_family.h" |
---|
[13] | 5 | |
---|
[477] | 6 | namespace bdm { |
---|
[13] | 7 | |
---|
[32] | 8 | Uniform_RNG UniRNG; |
---|
| 9 | Normal_RNG NorRNG; |
---|
| 10 | Gamma_RNG GamRNG; |
---|
| 11 | |
---|
[13] | 12 | using std::cout; |
---|
| 13 | |
---|
[487] | 14 | /////////// |
---|
| 15 | |
---|
[679] | 16 | void BMEF::bayes( const vec &yt, const vec &cond ) { |
---|
| 17 | this->bayes_weighted ( yt, cond, 1.0 ); |
---|
[477] | 18 | }; |
---|
[170] | 19 | |
---|
[629] | 20 | void egiw::set_parameters (int dimx0, ldmat V0, double nu0) { |
---|
| 21 | dimx = dimx0; |
---|
| 22 | nPsi = V0.rows() - dimx; |
---|
| 23 | dim = dimx * (dimx + nPsi); // size(R) + size(Theta) |
---|
| 24 | |
---|
| 25 | V = V0; |
---|
| 26 | if (nu0 < 0) { |
---|
| 27 | nu = 0.1 + nPsi + 2 * dimx + 2; // +2 assures finite expected value of R |
---|
| 28 | // terms before that are sufficient for finite normalization |
---|
| 29 | } else { |
---|
| 30 | nu = nu0; |
---|
| 31 | } |
---|
| 32 | } |
---|
| 33 | |
---|
[96] | 34 | vec egiw::sample() const { |
---|
[565] | 35 | bdm_warning ( "Function not implemented" ); |
---|
[168] | 36 | return vec_1 ( 0.0 ); |
---|
[96] | 37 | } |
---|
| 38 | |
---|
[211] | 39 | double egiw::evallog_nn ( const vec &val ) const { |
---|
[477] | 40 | int vend = val.length() - 1; |
---|
[168] | 41 | |
---|
[477] | 42 | if ( dimx == 1 ) { //same as the following, just quicker. |
---|
[170] | 43 | double r = val ( vend ); //last entry! |
---|
[477] | 44 | if ( r < 0 ) return -inf; |
---|
| 45 | vec Psi ( nPsi + dimx ); |
---|
[168] | 46 | Psi ( 0 ) = -1.0; |
---|
[477] | 47 | Psi.set_subvector ( 1, val ( 0, vend - 1 ) ); // fill the rest |
---|
[168] | 48 | |
---|
[477] | 49 | double Vq = V.qform ( Psi ); |
---|
| 50 | return -0.5* ( nu*log ( r ) + Vq / r ); |
---|
| 51 | } else { |
---|
| 52 | mat Th = reshape ( val ( 0, nPsi * dimx - 1 ), nPsi, dimx ); |
---|
| 53 | fsqmat R ( reshape ( val ( nPsi*dimx, vend ), dimx, dimx ) ); |
---|
| 54 | double ldetR = R.logdet(); |
---|
| 55 | if ( ldetR ) return -inf; |
---|
| 56 | mat Tmp = concat_vertical ( -eye ( dimx ), Th ); |
---|
[270] | 57 | fsqmat iR ( dimx ); |
---|
[168] | 58 | R.inv ( iR ); |
---|
[170] | 59 | |
---|
[395] | 60 | return -0.5* ( nu*ldetR + trace ( iR.to_mat() *Tmp.T() *V.to_mat() *Tmp ) ); |
---|
[168] | 61 | } |
---|
[96] | 62 | } |
---|
| 63 | |
---|
[168] | 64 | double egiw::lognc() const { |
---|
[96] | 65 | const vec& D = V._D(); |
---|
| 66 | |
---|
[477] | 67 | double m = nu - nPsi - dimx - 1; |
---|
[168] | 68 | #define log2 0.693147180559945286226763983 |
---|
| 69 | #define logpi 1.144729885849400163877476189 |
---|
| 70 | #define log2pi 1.83787706640935 |
---|
[178] | 71 | #define Inf std::numeric_limits<double>::infinity() |
---|
[168] | 72 | |
---|
[477] | 73 | double nkG = 0.5 * dimx * ( -nPsi * log2pi + sum ( log ( D ( dimx, D.length() - 1 ) ) ) ); |
---|
[168] | 74 | // temporary for lgamma in Wishart |
---|
[477] | 75 | double lg = 0; |
---|
| 76 | for ( int i = 0; i < dimx; i++ ) { |
---|
| 77 | lg += lgamma ( 0.5 * ( m - i ) ); |
---|
| 78 | } |
---|
[168] | 79 | |
---|
[477] | 80 | double nkW = 0.5 * ( m * sum ( log ( D ( 0, dimx - 1 ) ) ) ) \ |
---|
| 81 | - 0.5 * dimx * ( m * log2 + 0.5 * ( dimx - 1 ) * log2pi ) - lg; |
---|
[168] | 82 | |
---|
[637] | 83 | // bdm_assert_debug ( ( ( -nkG - nkW ) > -Inf ) && ( ( -nkG - nkW ) < Inf ), "ARX improper" ); |
---|
[665] | 84 | if (-nkG - nkW==Inf){ |
---|
| 85 | cout << "??" <<endl; |
---|
| 86 | } |
---|
[477] | 87 | return -nkG - nkW; |
---|
[96] | 88 | } |
---|
| 89 | |
---|
[330] | 90 | vec egiw::est_theta() const { |
---|
[477] | 91 | if ( dimx == 1 ) { |
---|
[330] | 92 | const mat &L = V._L(); |
---|
| 93 | int end = L.rows() - 1; |
---|
| 94 | |
---|
[477] | 95 | mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) ); |
---|
[330] | 96 | |
---|
[477] | 97 | vec L0 = L.get_col ( 0 ); |
---|
[330] | 98 | |
---|
[477] | 99 | return iLsub * L0 ( 1, end ); |
---|
| 100 | } else { |
---|
[565] | 101 | bdm_error ( "ERROR: est_theta() not implemented for dimx>1" ); |
---|
| 102 | return vec(); |
---|
[330] | 103 | } |
---|
| 104 | } |
---|
| 105 | |
---|
| 106 | ldmat egiw::est_theta_cov() const { |
---|
| 107 | if ( dimx == 1 ) { |
---|
| 108 | const mat &L = V._L(); |
---|
| 109 | const vec &D = V._D(); |
---|
| 110 | int end = D.length() - 1; |
---|
| 111 | |
---|
[477] | 112 | mat Lsub = L ( 1, end, 1, end ); |
---|
| 113 | mat Dsub = diag ( D ( 1, end ) ); |
---|
[330] | 114 | |
---|
[477] | 115 | return inv ( transpose ( Lsub ) * Dsub * Lsub ); |
---|
[330] | 116 | |
---|
[477] | 117 | } else { |
---|
[565] | 118 | bdm_error ( "ERROR: est_theta_cov() not implemented for dimx>1" ); |
---|
| 119 | return ldmat(); |
---|
[330] | 120 | } |
---|
| 121 | |
---|
| 122 | } |
---|
| 123 | |
---|
[96] | 124 | vec egiw::mean() const { |
---|
[168] | 125 | |
---|
[477] | 126 | if ( dimx == 1 ) { |
---|
| 127 | const vec &D = V._D(); |
---|
| 128 | int end = D.length() - 1; |
---|
[170] | 129 | |
---|
[270] | 130 | vec m ( dim ); |
---|
[330] | 131 | m.set_subvector ( 0, est_theta() ); |
---|
[477] | 132 | m ( end ) = D ( 0 ) / ( nu - nPsi - 2 * dimx - 2 ); |
---|
[170] | 133 | return m; |
---|
[477] | 134 | } else { |
---|
[170] | 135 | mat M; |
---|
| 136 | mat R; |
---|
[477] | 137 | mean_mat ( M, R ); |
---|
[629] | 138 | return concat( cvectorize ( M),cvectorize( R ) ); |
---|
[168] | 139 | } |
---|
[170] | 140 | |
---|
| 141 | } |
---|
[262] | 142 | |
---|
| 143 | vec egiw::variance() const { |
---|
[629] | 144 | int l = V.rows(); |
---|
| 145 | // cut out rest of lower-right part of V |
---|
| 146 | const ldmat tmp ( V, linspace ( dimx, l - 1 ) ); |
---|
| 147 | // invert it |
---|
| 148 | ldmat itmp ( l ); |
---|
| 149 | tmp.inv ( itmp ); |
---|
| 150 | |
---|
| 151 | // following Wikipedia notation |
---|
| 152 | // m=nu-nPsi-dimx-1, p=dimx |
---|
| 153 | double mp1p=nu-nPsi-2*dimx; // m-p+1 |
---|
| 154 | double mp1m=mp1p-2; // m-p-1 |
---|
| 155 | |
---|
[477] | 156 | if ( dimx == 1 ) { |
---|
[629] | 157 | double cove = V._D() ( 0 ) / mp1m ; |
---|
[477] | 158 | |
---|
| 159 | vec var ( l ); |
---|
| 160 | var.set_subvector ( 0, diag ( itmp.to_mat() ) *cove ); |
---|
[629] | 161 | var ( l - 1 ) = cove * cove / (mp1m-2); |
---|
[262] | 162 | return var; |
---|
[477] | 163 | } else { |
---|
[629] | 164 | ldmat Vll( V, linspace(0,dimx-1)); // top-left part of V |
---|
| 165 | mat Y=Vll.to_mat(); |
---|
| 166 | mat varY(Y.rows(), Y.cols()); |
---|
| 167 | |
---|
| 168 | double denom = (mp1p-1)*mp1m*mp1m*(mp1m-2); // (m-p)(m-p-1)^2(m-p-3) |
---|
| 169 | |
---|
| 170 | int i,j; |
---|
| 171 | for ( i=0; i<Y.rows(); i++){ |
---|
| 172 | for ( j=0; j<Y.cols(); j++){ |
---|
| 173 | varY(i,j) = (mp1p*Y(i,j)*Y(i,j) + mp1m * Y(i,i)* Y(j,j)) /denom; |
---|
| 174 | } |
---|
| 175 | } |
---|
| 176 | vec mean_dR = diag(Y)/mp1m; // corresponds to cove |
---|
| 177 | vec var_th=diag ( itmp.to_mat() ); |
---|
| 178 | vec var_Th ( mean_dR.length()*var_th.length() ); |
---|
| 179 | // diagonal of diag(mean_dR) \kron diag(var_th) |
---|
| 180 | for (int i=0; i<mean_dR.length(); i++){ |
---|
| 181 | var_Th.set_subvector(i*var_th.length(), var_th*mean_dR(i)); |
---|
| 182 | } |
---|
| 183 | |
---|
| 184 | return concat(var_Th, cvectorize(varY)); |
---|
[262] | 185 | } |
---|
| 186 | } |
---|
| 187 | |
---|
[225] | 188 | void egiw::mean_mat ( mat &M, mat&R ) const { |
---|
[477] | 189 | const mat &L = V._L(); |
---|
| 190 | const vec &D = V._D(); |
---|
| 191 | int end = L.rows() - 1; |
---|
[225] | 192 | |
---|
[477] | 193 | ldmat ldR ( L ( 0, dimx - 1, 0, dimx - 1 ), D ( 0, dimx - 1 ) / ( nu - nPsi - 2*dimx - 2 ) ); //exp val of R |
---|
| 194 | mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) ); |
---|
[225] | 195 | |
---|
[170] | 196 | // set mean value |
---|
[477] | 197 | mat Lpsi = L ( dimx, end, 0, dimx - 1 ); |
---|
| 198 | M = iLsub * Lpsi; |
---|
| 199 | R = ldR.to_mat() ; |
---|
[96] | 200 | } |
---|
| 201 | |
---|
[32] | 202 | vec egamma::sample() const { |
---|
[477] | 203 | vec smp ( dim ); |
---|
[32] | 204 | int i; |
---|
| 205 | |
---|
[477] | 206 | for ( i = 0; i < dim; i++ ) { |
---|
| 207 | if ( beta ( i ) > std::numeric_limits<double>::epsilon() ) { |
---|
| 208 | GamRNG.setup ( alpha ( i ), beta ( i ) ); |
---|
| 209 | } else { |
---|
| 210 | GamRNG.setup ( alpha ( i ), std::numeric_limits<double>::epsilon() ); |
---|
[225] | 211 | } |
---|
[168] | 212 | #pragma omp critical |
---|
[32] | 213 | smp ( i ) = GamRNG(); |
---|
| 214 | } |
---|
| 215 | |
---|
| 216 | return smp; |
---|
| 217 | } |
---|
| 218 | |
---|
[102] | 219 | // mat egamma::sample ( int N ) const { |
---|
| 220 | // mat Smp ( rv.count(),N ); |
---|
| 221 | // int i,j; |
---|
[168] | 222 | // |
---|
[102] | 223 | // for ( i=0; i<rv.count(); i++ ) { |
---|
| 224 | // GamRNG.setup ( alpha ( i ),beta ( i ) ); |
---|
[168] | 225 | // |
---|
[102] | 226 | // for ( j=0; j<N; j++ ) { |
---|
| 227 | // Smp ( i,j ) = GamRNG(); |
---|
| 228 | // } |
---|
| 229 | // } |
---|
[168] | 230 | // |
---|
[102] | 231 | // return Smp; |
---|
| 232 | // } |
---|
[32] | 233 | |
---|
[211] | 234 | double egamma::evallog ( const vec &val ) const { |
---|
[96] | 235 | double res = 0.0; //the rest will be added |
---|
| 236 | int i; |
---|
| 237 | |
---|
[477] | 238 | if ( any ( val <= 0. ) ) return -inf; |
---|
| 239 | if ( any ( beta <= 0. ) ) return -inf; |
---|
| 240 | for ( i = 0; i < dim; i++ ) { |
---|
| 241 | res += ( alpha ( i ) - 1 ) * std::log ( val ( i ) ) - beta ( i ) * val ( i ); |
---|
[96] | 242 | } |
---|
[477] | 243 | double tmp = res - lognc();; |
---|
[565] | 244 | bdm_assert_debug ( std::isfinite ( tmp ), "Infinite value" ); |
---|
[214] | 245 | return tmp; |
---|
[96] | 246 | } |
---|
| 247 | |
---|
| 248 | double egamma::lognc() const { |
---|
[32] | 249 | double res = 0.0; //will be added |
---|
| 250 | int i; |
---|
| 251 | |
---|
[477] | 252 | for ( i = 0; i < dim; i++ ) { |
---|
| 253 | res += lgamma ( alpha ( i ) ) - alpha ( i ) * std::log ( beta ( i ) ) ; |
---|
[32] | 254 | } |
---|
| 255 | |
---|
| 256 | return res; |
---|
| 257 | } |
---|
| 258 | |
---|
[477] | 259 | void mgamma::set_parameters ( double k0, const vec &beta0 ) { |
---|
[461] | 260 | k = k0; |
---|
[487] | 261 | iepdf.set_parameters ( k * ones ( beta0.length() ), beta0 ); |
---|
| 262 | dimc = iepdf.dimension(); |
---|
[678] | 263 | dim = iepdf.dimension(); |
---|
[461] | 264 | } |
---|
[32] | 265 | |
---|
[637] | 266 | void eEmp::resample ( ivec &ind, RESAMPLING_METHOD method ) { |
---|
| 267 | ind = zeros_i ( n ); |
---|
[32] | 268 | ivec N_babies = zeros_i ( n ); |
---|
| 269 | vec cumDist = cumsum ( w ); |
---|
| 270 | vec u ( n ); |
---|
[477] | 271 | int i, j, parent; |
---|
[32] | 272 | double u0; |
---|
| 273 | |
---|
| 274 | switch ( method ) { |
---|
[477] | 275 | case MULTINOMIAL: |
---|
| 276 | u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n ); |
---|
[32] | 277 | |
---|
[477] | 278 | for ( i = n - 2; i >= 0; i-- ) { |
---|
| 279 | u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) ); |
---|
| 280 | } |
---|
[32] | 281 | |
---|
[477] | 282 | break; |
---|
[32] | 283 | |
---|
[477] | 284 | case STRATIFIED: |
---|
[32] | 285 | |
---|
[477] | 286 | for ( i = 0; i < n; i++ ) { |
---|
| 287 | u ( i ) = ( i + UniRNG.sample() ) / n; |
---|
| 288 | } |
---|
[32] | 289 | |
---|
[477] | 290 | break; |
---|
[32] | 291 | |
---|
[477] | 292 | case SYSTEMATIC: |
---|
| 293 | u0 = UniRNG.sample(); |
---|
[32] | 294 | |
---|
[477] | 295 | for ( i = 0; i < n; i++ ) { |
---|
| 296 | u ( i ) = ( i + u0 ) / n; |
---|
| 297 | } |
---|
[32] | 298 | |
---|
[477] | 299 | break; |
---|
[32] | 300 | |
---|
[477] | 301 | default: |
---|
[565] | 302 | bdm_error ( "PF::resample(): Unknown resampling method" ); |
---|
[32] | 303 | } |
---|
| 304 | |
---|
| 305 | // U is now full |
---|
| 306 | j = 0; |
---|
| 307 | |
---|
[477] | 308 | for ( i = 0; i < n; i++ ) { |
---|
[32] | 309 | while ( u ( i ) > cumDist ( j ) ) j++; |
---|
| 310 | |
---|
| 311 | N_babies ( j ) ++; |
---|
| 312 | } |
---|
| 313 | // We have assigned new babies for each Particle |
---|
| 314 | // Now, we fill the resulting index such that: |
---|
| 315 | // * particles with at least one baby should not move * |
---|
| 316 | // This assures that reassignment can be done inplace; |
---|
| 317 | |
---|
| 318 | // find the first parent; |
---|
[477] | 319 | parent = 0; |
---|
| 320 | while ( N_babies ( parent ) == 0 ) parent++; |
---|
[32] | 321 | |
---|
| 322 | // Build index |
---|
[477] | 323 | for ( i = 0; i < n; i++ ) { |
---|
[32] | 324 | if ( N_babies ( i ) > 0 ) { |
---|
| 325 | ind ( i ) = i; |
---|
| 326 | N_babies ( i ) --; //this index was now replicated; |
---|
[477] | 327 | } else { |
---|
[32] | 328 | // test if the parent has been fully replicated |
---|
| 329 | // if yes, find the next one |
---|
[477] | 330 | while ( ( N_babies ( parent ) == 0 ) || ( N_babies ( parent ) == 1 && parent > i ) ) parent++; |
---|
[32] | 331 | |
---|
| 332 | // Replicate parent |
---|
| 333 | ind ( i ) = parent; |
---|
| 334 | |
---|
| 335 | N_babies ( parent ) --; //this index was now replicated; |
---|
| 336 | } |
---|
| 337 | |
---|
| 338 | } |
---|
| 339 | |
---|
[637] | 340 | // copy the internals according to ind |
---|
[477] | 341 | for ( i = 0; i < n; i++ ) { |
---|
| 342 | if ( ind ( i ) != i ) { |
---|
| 343 | samples ( i ) = samples ( ind ( i ) ); |
---|
[32] | 344 | } |
---|
[477] | 345 | w ( i ) = 1.0 / n; |
---|
[32] | 346 | } |
---|
| 347 | } |
---|
| 348 | |
---|
[488] | 349 | void eEmp::set_statistics ( const vec &w0, const epdf &epdf0 ) { |
---|
| 350 | dim = epdf0.dimension(); |
---|
[477] | 351 | w = w0; |
---|
| 352 | w /= sum ( w0 );//renormalize |
---|
| 353 | n = w.length(); |
---|
[32] | 354 | samples.set_size ( n ); |
---|
| 355 | |
---|
[477] | 356 | for ( int i = 0; i < n; i++ ) { |
---|
[488] | 357 | samples ( i ) = epdf0.sample(); |
---|
[477] | 358 | } |
---|
[32] | 359 | } |
---|
[178] | 360 | |
---|
[225] | 361 | void eEmp::set_samples ( const epdf* epdf0 ) { |
---|
[477] | 362 | w = 1; |
---|
| 363 | w /= sum ( w );//renormalize |
---|
[178] | 364 | |
---|
[477] | 365 | for ( int i = 0; i < n; i++ ) { |
---|
| 366 | samples ( i ) = epdf0->sample(); |
---|
| 367 | } |
---|
[178] | 368 | } |
---|
| 369 | |
---|
[477] | 370 | void migamma_ref::from_setting ( const Setting &set ) { |
---|
[357] | 371 | vec ref; |
---|
[477] | 372 | UI::get ( ref, set, "ref" , UI::compulsory ); |
---|
| 373 | set_parameters ( set["k"], ref, set["l"] ); |
---|
[357] | 374 | } |
---|
| 375 | |
---|
[477] | 376 | void mlognorm::from_setting ( const Setting &set ) { |
---|
| 377 | vec mu0; |
---|
| 378 | UI::get ( mu0, set, "mu0", UI::compulsory ); |
---|
| 379 | set_parameters ( mu0.length(), set["k"] ); |
---|
| 380 | condition ( mu0 ); |
---|
[357] | 381 | } |
---|
| 382 | |
---|
[262] | 383 | }; |
---|