[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 | |
---|
[477] | 16 | void BMEF::bayes ( const vec &dt ) { |
---|
| 17 | this->bayes ( dt, 1.0 ); |
---|
| 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 | |
---|
[565] | 83 | bdm_assert_debug ( ( ( -nkG - nkW ) > -Inf ) && ( ( -nkG - nkW ) < Inf ), "ARX improper" ); |
---|
[477] | 84 | return -nkG - nkW; |
---|
[96] | 85 | } |
---|
| 86 | |
---|
[330] | 87 | vec egiw::est_theta() const { |
---|
[477] | 88 | if ( dimx == 1 ) { |
---|
[330] | 89 | const mat &L = V._L(); |
---|
| 90 | int end = L.rows() - 1; |
---|
| 91 | |
---|
[477] | 92 | mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) ); |
---|
[330] | 93 | |
---|
[477] | 94 | vec L0 = L.get_col ( 0 ); |
---|
[330] | 95 | |
---|
[477] | 96 | return iLsub * L0 ( 1, end ); |
---|
| 97 | } else { |
---|
[565] | 98 | bdm_error ( "ERROR: est_theta() not implemented for dimx>1" ); |
---|
| 99 | return vec(); |
---|
[330] | 100 | } |
---|
| 101 | } |
---|
| 102 | |
---|
| 103 | ldmat egiw::est_theta_cov() const { |
---|
| 104 | if ( dimx == 1 ) { |
---|
| 105 | const mat &L = V._L(); |
---|
| 106 | const vec &D = V._D(); |
---|
| 107 | int end = D.length() - 1; |
---|
| 108 | |
---|
[477] | 109 | mat Lsub = L ( 1, end, 1, end ); |
---|
| 110 | mat Dsub = diag ( D ( 1, end ) ); |
---|
[330] | 111 | |
---|
[477] | 112 | return inv ( transpose ( Lsub ) * Dsub * Lsub ); |
---|
[330] | 113 | |
---|
[477] | 114 | } else { |
---|
[565] | 115 | bdm_error ( "ERROR: est_theta_cov() not implemented for dimx>1" ); |
---|
| 116 | return ldmat(); |
---|
[330] | 117 | } |
---|
| 118 | |
---|
| 119 | } |
---|
| 120 | |
---|
[96] | 121 | vec egiw::mean() const { |
---|
[168] | 122 | |
---|
[477] | 123 | if ( dimx == 1 ) { |
---|
| 124 | const vec &D = V._D(); |
---|
| 125 | int end = D.length() - 1; |
---|
[170] | 126 | |
---|
[270] | 127 | vec m ( dim ); |
---|
[330] | 128 | m.set_subvector ( 0, est_theta() ); |
---|
[477] | 129 | m ( end ) = D ( 0 ) / ( nu - nPsi - 2 * dimx - 2 ); |
---|
[170] | 130 | return m; |
---|
[477] | 131 | } else { |
---|
[170] | 132 | mat M; |
---|
| 133 | mat R; |
---|
[477] | 134 | mean_mat ( M, R ); |
---|
[629] | 135 | return concat( cvectorize ( M),cvectorize( R ) ); |
---|
[168] | 136 | } |
---|
[170] | 137 | |
---|
| 138 | } |
---|
[262] | 139 | |
---|
| 140 | vec egiw::variance() const { |
---|
[629] | 141 | int l = V.rows(); |
---|
| 142 | // cut out rest of lower-right part of V |
---|
| 143 | const ldmat tmp ( V, linspace ( dimx, l - 1 ) ); |
---|
| 144 | // invert it |
---|
| 145 | ldmat itmp ( l ); |
---|
| 146 | tmp.inv ( itmp ); |
---|
| 147 | |
---|
| 148 | // following Wikipedia notation |
---|
| 149 | // m=nu-nPsi-dimx-1, p=dimx |
---|
| 150 | double mp1p=nu-nPsi-2*dimx; // m-p+1 |
---|
| 151 | double mp1m=mp1p-2; // m-p-1 |
---|
| 152 | |
---|
[477] | 153 | if ( dimx == 1 ) { |
---|
[629] | 154 | double cove = V._D() ( 0 ) / mp1m ; |
---|
[477] | 155 | |
---|
| 156 | vec var ( l ); |
---|
| 157 | var.set_subvector ( 0, diag ( itmp.to_mat() ) *cove ); |
---|
[629] | 158 | var ( l - 1 ) = cove * cove / (mp1m-2); |
---|
[262] | 159 | return var; |
---|
[477] | 160 | } else { |
---|
[629] | 161 | ldmat Vll( V, linspace(0,dimx-1)); // top-left part of V |
---|
| 162 | mat Y=Vll.to_mat(); |
---|
| 163 | mat varY(Y.rows(), Y.cols()); |
---|
| 164 | |
---|
| 165 | double denom = (mp1p-1)*mp1m*mp1m*(mp1m-2); // (m-p)(m-p-1)^2(m-p-3) |
---|
| 166 | |
---|
| 167 | int i,j; |
---|
| 168 | for ( i=0; i<Y.rows(); i++){ |
---|
| 169 | for ( j=0; j<Y.cols(); j++){ |
---|
| 170 | varY(i,j) = (mp1p*Y(i,j)*Y(i,j) + mp1m * Y(i,i)* Y(j,j)) /denom; |
---|
| 171 | } |
---|
| 172 | } |
---|
| 173 | vec mean_dR = diag(Y)/mp1m; // corresponds to cove |
---|
| 174 | vec var_th=diag ( itmp.to_mat() ); |
---|
| 175 | vec var_Th ( mean_dR.length()*var_th.length() ); |
---|
| 176 | // diagonal of diag(mean_dR) \kron diag(var_th) |
---|
| 177 | for (int i=0; i<mean_dR.length(); i++){ |
---|
| 178 | var_Th.set_subvector(i*var_th.length(), var_th*mean_dR(i)); |
---|
| 179 | } |
---|
| 180 | |
---|
| 181 | return concat(var_Th, cvectorize(varY)); |
---|
[262] | 182 | } |
---|
| 183 | } |
---|
| 184 | |
---|
[225] | 185 | void egiw::mean_mat ( mat &M, mat&R ) const { |
---|
[477] | 186 | const mat &L = V._L(); |
---|
| 187 | const vec &D = V._D(); |
---|
| 188 | int end = L.rows() - 1; |
---|
[225] | 189 | |
---|
[477] | 190 | ldmat ldR ( L ( 0, dimx - 1, 0, dimx - 1 ), D ( 0, dimx - 1 ) / ( nu - nPsi - 2*dimx - 2 ) ); //exp val of R |
---|
| 191 | mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) ); |
---|
[225] | 192 | |
---|
[170] | 193 | // set mean value |
---|
[477] | 194 | mat Lpsi = L ( dimx, end, 0, dimx - 1 ); |
---|
| 195 | M = iLsub * Lpsi; |
---|
| 196 | R = ldR.to_mat() ; |
---|
[96] | 197 | } |
---|
| 198 | |
---|
[32] | 199 | vec egamma::sample() const { |
---|
[477] | 200 | vec smp ( dim ); |
---|
[32] | 201 | int i; |
---|
| 202 | |
---|
[477] | 203 | for ( i = 0; i < dim; i++ ) { |
---|
| 204 | if ( beta ( i ) > std::numeric_limits<double>::epsilon() ) { |
---|
| 205 | GamRNG.setup ( alpha ( i ), beta ( i ) ); |
---|
| 206 | } else { |
---|
| 207 | GamRNG.setup ( alpha ( i ), std::numeric_limits<double>::epsilon() ); |
---|
[225] | 208 | } |
---|
[168] | 209 | #pragma omp critical |
---|
[32] | 210 | smp ( i ) = GamRNG(); |
---|
| 211 | } |
---|
| 212 | |
---|
| 213 | return smp; |
---|
| 214 | } |
---|
| 215 | |
---|
[102] | 216 | // mat egamma::sample ( int N ) const { |
---|
| 217 | // mat Smp ( rv.count(),N ); |
---|
| 218 | // int i,j; |
---|
[168] | 219 | // |
---|
[102] | 220 | // for ( i=0; i<rv.count(); i++ ) { |
---|
| 221 | // GamRNG.setup ( alpha ( i ),beta ( i ) ); |
---|
[168] | 222 | // |
---|
[102] | 223 | // for ( j=0; j<N; j++ ) { |
---|
| 224 | // Smp ( i,j ) = GamRNG(); |
---|
| 225 | // } |
---|
| 226 | // } |
---|
[168] | 227 | // |
---|
[102] | 228 | // return Smp; |
---|
| 229 | // } |
---|
[32] | 230 | |
---|
[211] | 231 | double egamma::evallog ( const vec &val ) const { |
---|
[96] | 232 | double res = 0.0; //the rest will be added |
---|
| 233 | int i; |
---|
| 234 | |
---|
[477] | 235 | if ( any ( val <= 0. ) ) return -inf; |
---|
| 236 | if ( any ( beta <= 0. ) ) return -inf; |
---|
| 237 | for ( i = 0; i < dim; i++ ) { |
---|
| 238 | res += ( alpha ( i ) - 1 ) * std::log ( val ( i ) ) - beta ( i ) * val ( i ); |
---|
[96] | 239 | } |
---|
[477] | 240 | double tmp = res - lognc();; |
---|
[565] | 241 | bdm_assert_debug ( std::isfinite ( tmp ), "Infinite value" ); |
---|
[214] | 242 | return tmp; |
---|
[96] | 243 | } |
---|
| 244 | |
---|
| 245 | double egamma::lognc() const { |
---|
[32] | 246 | double res = 0.0; //will be added |
---|
| 247 | int i; |
---|
| 248 | |
---|
[477] | 249 | for ( i = 0; i < dim; i++ ) { |
---|
| 250 | res += lgamma ( alpha ( i ) ) - alpha ( i ) * std::log ( beta ( i ) ) ; |
---|
[32] | 251 | } |
---|
| 252 | |
---|
| 253 | return res; |
---|
| 254 | } |
---|
| 255 | |
---|
[477] | 256 | void mgamma::set_parameters ( double k0, const vec &beta0 ) { |
---|
[461] | 257 | k = k0; |
---|
[487] | 258 | iepdf.set_parameters ( k * ones ( beta0.length() ), beta0 ); |
---|
| 259 | dimc = iepdf.dimension(); |
---|
[461] | 260 | } |
---|
[32] | 261 | |
---|
[477] | 262 | ivec eEmp::resample ( RESAMPLING_METHOD method ) { |
---|
| 263 | ivec ind = zeros_i ( n ); |
---|
[32] | 264 | ivec N_babies = zeros_i ( n ); |
---|
| 265 | vec cumDist = cumsum ( w ); |
---|
| 266 | vec u ( n ); |
---|
[477] | 267 | int i, j, parent; |
---|
[32] | 268 | double u0; |
---|
| 269 | |
---|
| 270 | switch ( method ) { |
---|
[477] | 271 | case MULTINOMIAL: |
---|
| 272 | u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n ); |
---|
[32] | 273 | |
---|
[477] | 274 | for ( i = n - 2; i >= 0; i-- ) { |
---|
| 275 | u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) ); |
---|
| 276 | } |
---|
[32] | 277 | |
---|
[477] | 278 | break; |
---|
[32] | 279 | |
---|
[477] | 280 | case STRATIFIED: |
---|
[32] | 281 | |
---|
[477] | 282 | for ( i = 0; i < n; i++ ) { |
---|
| 283 | u ( i ) = ( i + UniRNG.sample() ) / n; |
---|
| 284 | } |
---|
[32] | 285 | |
---|
[477] | 286 | break; |
---|
[32] | 287 | |
---|
[477] | 288 | case SYSTEMATIC: |
---|
| 289 | u0 = UniRNG.sample(); |
---|
[32] | 290 | |
---|
[477] | 291 | for ( i = 0; i < n; i++ ) { |
---|
| 292 | u ( i ) = ( i + u0 ) / n; |
---|
| 293 | } |
---|
[32] | 294 | |
---|
[477] | 295 | break; |
---|
[32] | 296 | |
---|
[477] | 297 | default: |
---|
[565] | 298 | bdm_error ( "PF::resample(): Unknown resampling method" ); |
---|
[32] | 299 | } |
---|
| 300 | |
---|
| 301 | // U is now full |
---|
| 302 | j = 0; |
---|
| 303 | |
---|
[477] | 304 | for ( i = 0; i < n; i++ ) { |
---|
[32] | 305 | while ( u ( i ) > cumDist ( j ) ) j++; |
---|
| 306 | |
---|
| 307 | N_babies ( j ) ++; |
---|
| 308 | } |
---|
| 309 | // We have assigned new babies for each Particle |
---|
| 310 | // Now, we fill the resulting index such that: |
---|
| 311 | // * particles with at least one baby should not move * |
---|
| 312 | // This assures that reassignment can be done inplace; |
---|
| 313 | |
---|
| 314 | // find the first parent; |
---|
[477] | 315 | parent = 0; |
---|
| 316 | while ( N_babies ( parent ) == 0 ) parent++; |
---|
[32] | 317 | |
---|
| 318 | // Build index |
---|
[477] | 319 | for ( i = 0; i < n; i++ ) { |
---|
[32] | 320 | if ( N_babies ( i ) > 0 ) { |
---|
| 321 | ind ( i ) = i; |
---|
| 322 | N_babies ( i ) --; //this index was now replicated; |
---|
[477] | 323 | } else { |
---|
[32] | 324 | // test if the parent has been fully replicated |
---|
| 325 | // if yes, find the next one |
---|
[477] | 326 | while ( ( N_babies ( parent ) == 0 ) || ( N_babies ( parent ) == 1 && parent > i ) ) parent++; |
---|
[32] | 327 | |
---|
| 328 | // Replicate parent |
---|
| 329 | ind ( i ) = parent; |
---|
| 330 | |
---|
| 331 | N_babies ( parent ) --; //this index was now replicated; |
---|
| 332 | } |
---|
| 333 | |
---|
| 334 | } |
---|
| 335 | |
---|
| 336 | // copy the internals according to ind |
---|
[477] | 337 | for ( i = 0; i < n; i++ ) { |
---|
| 338 | if ( ind ( i ) != i ) { |
---|
| 339 | samples ( i ) = samples ( ind ( i ) ); |
---|
[32] | 340 | } |
---|
[477] | 341 | w ( i ) = 1.0 / n; |
---|
[32] | 342 | } |
---|
| 343 | |
---|
| 344 | return ind; |
---|
| 345 | } |
---|
| 346 | |
---|
[488] | 347 | void eEmp::set_statistics ( const vec &w0, const epdf &epdf0 ) { |
---|
| 348 | dim = epdf0.dimension(); |
---|
[477] | 349 | w = w0; |
---|
| 350 | w /= sum ( w0 );//renormalize |
---|
| 351 | n = w.length(); |
---|
[32] | 352 | samples.set_size ( n ); |
---|
| 353 | |
---|
[477] | 354 | for ( int i = 0; i < n; i++ ) { |
---|
[488] | 355 | samples ( i ) = epdf0.sample(); |
---|
[477] | 356 | } |
---|
[32] | 357 | } |
---|
[178] | 358 | |
---|
[225] | 359 | void eEmp::set_samples ( const epdf* epdf0 ) { |
---|
[477] | 360 | w = 1; |
---|
| 361 | w /= sum ( w );//renormalize |
---|
[178] | 362 | |
---|
[477] | 363 | for ( int i = 0; i < n; i++ ) { |
---|
| 364 | samples ( i ) = epdf0->sample(); |
---|
| 365 | } |
---|
[178] | 366 | } |
---|
| 367 | |
---|
[477] | 368 | void migamma_ref::from_setting ( const Setting &set ) { |
---|
[357] | 369 | vec ref; |
---|
[477] | 370 | UI::get ( ref, set, "ref" , UI::compulsory ); |
---|
| 371 | set_parameters ( set["k"], ref, set["l"] ); |
---|
[357] | 372 | } |
---|
| 373 | |
---|
[477] | 374 | void mlognorm::from_setting ( const Setting &set ) { |
---|
| 375 | vec mu0; |
---|
| 376 | UI::get ( mu0, set, "mu0", UI::compulsory ); |
---|
| 377 | set_parameters ( mu0.length(), set["k"] ); |
---|
| 378 | condition ( mu0 ); |
---|
[357] | 379 | } |
---|
| 380 | |
---|
[262] | 381 | }; |
---|