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