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