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