root/library/bdm/stat/exp_family.cpp @ 826

Revision 826, 13.5 kB (checked in by smidl, 14 years ago)

win32 fixes

  • Property svn:eol-style set to native
Line 
1#include <math.h>
2
3#include <itpp/base/bessel.h>
4#include "exp_family.h"
5
6namespace bdm {
7
8Uniform_RNG UniRNG;
9Normal_RNG NorRNG;
10Gamma_RNG GamRNG;
11
12using std::cout;
13
14///////////
15
16void BMEF::bayes ( const vec &yt, const vec &cond ) {
17        this->bayes_weighted ( yt, cond, 1.0 );
18};
19
20void 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
34vec egiw::sample() const {
35        mat M;
36        chmat R;
37        sample_mat ( M, R );
38
39        return concat ( cvectorize ( M ), cvectorize ( R.to_mat() ) );
40}
41
42mat egiw::sample_mat ( int n ) const {
43        // TODO - correct approach - convert to product of norm * Wishart
44        mat M;
45        ldmat Vz;
46        ldmat Lam;
47        factorize ( M, Vz, Lam );
48
49        chmat ChLam ( Lam.to_mat() );
50        chmat iChLam;
51        ChLam.inv ( iChLam );
52
53        eWishartCh Omega; //inverse Wishart, result is R,
54        Omega.set_parameters ( iChLam, nu - 2*nPsi - dimx ); // 2*nPsi is there to match numercial simulations - check if analytically correct
55
56        mat OmChi;
57        mat Z ( M.rows(), M.cols() );
58
59        mat Mi;
60        mat RChiT;
61        mat tmp ( dimension(), n );
62        for ( int i = 0; i < n; i++ ) {
63                OmChi = Omega.sample_mat();
64                RChiT = inv ( OmChi );
65                Z = randn ( M.rows(), M.cols() );
66                Mi = M + RChiT * Z * inv ( Vz._L().T() * diag ( sqrt ( Vz._D() ) ) );
67
68                tmp.set_col ( i, concat ( cvectorize ( Mi ), cvectorize ( RChiT*RChiT.T() ) ) );
69        }
70        return tmp;
71}
72
73void egiw::sample_mat ( mat &Mi, chmat &Ri ) const {
74
75        // TODO - correct approach - convert to product of norm * Wishart
76        mat M;
77        ldmat Vz;
78        ldmat Lam;
79        factorize ( M, Vz, Lam );
80
81        chmat Ch;
82        Ch.setCh ( Lam._L() *diag ( sqrt ( Lam._D() ) ) );
83        chmat iCh;
84        Ch.inv ( iCh );
85
86        eWishartCh Omega; //inverse Wishart, result is R,
87        Omega.set_parameters ( iCh, nu - 2*nPsi - dimx ); // 2*nPsi is there to match numercial simulations - check if analytically correct
88
89        chmat Omi;
90        Omi.setCh ( Omega.sample_mat() );
91
92        if (M._datasize()>0){
93                mat Z = randn ( M.rows(), M.cols() );
94                Mi = M + Omi._Ch() * Z * inv ( Vz._L() * diag ( sqrt ( Vz._D() ) ) );
95        }
96        Omi.inv ( Ri );
97}
98
99double egiw::evallog_nn ( const vec &val ) const {
100        bdm_assert_debug(val.length()==dimx*(nPsi+dimx),"Incorrect cond in egiw::evallog_nn" );
101       
102        int vend = val.length() - 1;
103
104        if ( dimx == 1 ) { //same as the following, just quicker.
105                double r = val ( vend ); //last entry!
106                if ( r < 0 ) return -inf;
107                vec Psi ( nPsi + dimx );
108                Psi ( 0 ) = -1.0;
109                Psi.set_subvector ( 1, val ( 0, vend - 1 ) ); // fill the rest
110
111                double Vq = V.qform ( Psi );
112                return -0.5* ( nu*log ( r ) + Vq / r );
113        } else {
114                mat Tmp;
115                if (nPsi>0){
116                        mat Th = reshape ( val ( 0, nPsi * dimx - 1 ), nPsi, dimx );
117                        Tmp = concat_vertical ( -eye ( dimx ), Th );
118                } else {
119                        Tmp = -eye(dimx);
120                }
121                fsqmat R ( reshape ( val ( nPsi*dimx, vend ), dimx, dimx ) );
122                double ldetR = R.logdet();
123                if ( !std::isfinite(ldetR) ) return -inf;
124                fsqmat iR ( dimx );
125                R.inv ( iR );
126
127                return -0.5* ( nu*ldetR + trace ( iR.to_mat() *Tmp.T() *V.to_mat() *Tmp ) );
128        }
129}
130
131double egiw::lognc() const {
132        const vec& D = V._D();
133
134        double m = nu - nPsi - dimx - 1;
135#define log2  0.693147180559945286226763983
136#define logpi 1.144729885849400163877476189
137#define log2pi 1.83787706640935
138#define Inf std::numeric_limits<double>::infinity()
139
140        double nkG = 0.5 * dimx * ( -nPsi * log2pi + sum ( log ( D.mid ( dimx, nPsi ) ) ) );
141        // temporary for lgamma in Wishart
142        double lg = 0;
143        for ( int i = 0; i < dimx; i++ ) {
144                lg += lgamma ( 0.5 * ( m - i ) );
145        }
146
147        double nkW = 0.5 * ( m * sum ( log ( D ( 0, dimx - 1 ) ) ) ) \
148                     - 0.5 * dimx * ( m * log2 + 0.5 * ( dimx - 1 ) * log2pi )  - lg;
149
150//      bdm_assert_debug ( ( ( -nkG - nkW ) > -Inf ) && ( ( -nkG - nkW ) < Inf ), "ARX improper" );
151        if ( -nkG - nkW == Inf ) {
152                cout << "??" << endl;
153        }
154        return -nkG - nkW;
155}
156
157vec egiw::est_theta() const {
158        if ( dimx == 1 ) {
159                const mat &L = V._L();
160                int end = L.rows() - 1;
161
162                mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) );
163
164                vec L0 = L.get_col ( 0 );
165
166                return iLsub * L0 ( 1, end );
167        } else {
168                bdm_error ( "ERROR: est_theta() not implemented for dimx>1" );
169                return vec();
170        }
171}
172
173void egiw::factorize ( mat &M, ldmat &Vz, ldmat &Lam ) const {
174        const mat &L = V._L();
175        const vec &D = V._D();
176        int end = L.rows() - 1;
177
178        Lam = ldmat ( L ( 0, dimx - 1, 0, dimx - 1 ), D ( 0, dimx - 1 ) );  //exp val of R
179
180        if (dimx<=end){
181                Vz = ldmat ( L ( dimx, end, dimx, end ), D ( dimx, end ) );
182                mat iLsub = ltuinv ( Vz._L() );
183                // set mean value
184                mat Lpsi = L ( dimx, end, 0, dimx - 1 );
185                M = iLsub * Lpsi;
186        }
187/*      if ( 0 ) { // test with Peterka
188                mat VF = V.to_mat();
189                mat Vf = VF ( 0, dimx - 1, 0, dimx - 1 );
190                mat Vzf = VF ( dimx, end, 0, dimx - 1 );
191                mat VZ = VF ( dimx, end, dimx, end );
192
193                mat Lam2 = Vf - Vzf.T() * inv ( VZ ) * Vzf;
194        }*/
195}
196
197ldmat egiw::est_theta_cov() const {
198        if ( dimx == 1 ) {
199                const mat &L = V._L();
200                const vec &D = V._D();
201                int end = D.length() - 1;
202
203                mat Lsub = L ( 1, end, 1, end );
204//              mat Dsub = diag ( D ( 1, end ) );
205
206                ldmat LD ( inv ( Lsub ).T(), 1.0 / D ( 1, end ) );
207                return LD;
208
209        } else {
210                bdm_error ( "ERROR: est_theta_cov() not implemented for dimx>1" );
211                return ldmat();
212        }
213
214}
215
216vec egiw::mean() const {
217
218        if ( dimx == 1 ) {
219                const vec &D = V._D();
220                int end = D.length() - 1;
221
222                vec m ( dim );
223                m.set_subvector ( 0, est_theta() );
224                m ( end ) = D ( 0 ) / ( nu - nPsi - 2 * dimx - 2 );
225                return m;
226        } else {
227                mat M;
228                mat R;
229                mean_mat ( M, R );
230                return concat ( cvectorize ( M ), cvectorize ( R ) );
231        }
232
233}
234
235vec egiw::variance() const {
236        int l = V.rows();
237        // cut out rest of lower-right part of V
238        // invert it
239        ldmat itmp ( l );
240        if (dimx<l){
241                const ldmat tmp ( V, linspace ( dimx, l - 1 ) );
242                tmp.inv ( itmp );
243        }
244        // following Wikipedia notation
245        // m=nu-nPsi-dimx-1, p=dimx
246        double mp1p = nu - nPsi - 2 * dimx; // m-p+1
247        double mp1m = mp1p - 2;     // m-p-1
248
249        if ( dimx == 1 ) {
250                double cove = V._D() ( 0 ) / mp1m ;
251
252                vec var ( l );
253                var.set_subvector ( 0, diag ( itmp.to_mat() ) *cove );
254                var ( l - 1 ) = cove * cove / ( mp1m - 2 );
255                return var;
256        } else {
257                ldmat Vll ( V, linspace ( 0, dimx - 1 ) ); // top-left part of V
258                mat Y = Vll.to_mat();
259                mat varY ( Y.rows(), Y.cols() );
260
261                double denom = ( mp1p - 1 ) * mp1m * mp1m * ( mp1m - 2 );         // (m-p)(m-p-1)^2(m-p-3)
262
263                int i, j;
264                for ( i = 0; i < Y.rows(); i++ ) {
265                        for ( j = 0; j < Y.cols(); j++ ) {
266                                varY ( i, j ) = ( mp1p * Y ( i, j ) * Y ( i, j ) + mp1m * Y ( i, i ) * Y ( j, j ) ) / denom;
267                        }
268                }
269                vec mean_dR = diag ( Y ) / mp1m; // corresponds to cove
270                vec var_th = diag ( itmp.to_mat() );
271                vec var_Th ( mean_dR.length() *var_th.length() );
272                // diagonal of diag(mean_dR) \kron diag(var_th)
273                for ( int i = 0; i < mean_dR.length(); i++ ) {
274                        var_Th.set_subvector ( i*var_th.length(), var_th*mean_dR ( i ) );
275                }
276
277                return concat ( var_Th, cvectorize ( varY ) );
278        }
279}
280
281void egiw::mean_mat ( mat &M, mat&R ) const {
282        const mat &L = V._L();
283        const vec &D = V._D();
284        int end = L.rows() - 1;
285
286        ldmat ldR ( L ( 0, dimx - 1, 0, dimx - 1 ), D ( 0, dimx - 1 ) / ( nu - nPsi - 2*dimx - 2 ) ); //exp val of R
287
288        // set mean value
289        if (dimx<=end){
290                mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) );
291                mat Lpsi = L ( dimx, end, 0, dimx - 1 );
292                M = iLsub * Lpsi;
293        }
294        R = ldR.to_mat()  ;
295}
296
297void egiw::log_register ( bdm::logger& L, const string& prefix ) {
298        if ( log_level == 3 ) {
299                root::log_register ( L, prefix );
300                logrec->ids.set_length ( 2 );
301                int th_dim = dimension() - dimx * ( dimx + 1 ) / 2;
302                logrec->ids ( 0 ) = L.add_vector ( RV ( "", th_dim ), prefix + logrec->L.prefix_sep() + "mean" );
303                logrec->ids ( 1 ) = L.add_vector ( RV ( "", th_dim * th_dim ), prefix + logrec->L.prefix_sep() + "variance" );
304        } else {
305                epdf::log_register ( L, prefix );
306        }
307}
308
309void egiw::log_write() const {
310        if ( log_level == 3 ) {
311                mat M;
312                ldmat Lam;
313                ldmat Vz;
314                factorize ( M, Vz, Lam );
315                logrec->L.log_vector ( logrec->ids ( 0 ), est_theta() );
316                logrec->L.log_vector ( logrec->ids ( 1 ), cvectorize ( est_theta_cov().to_mat() ) );
317        } else {
318                epdf::log_write();
319        }
320
321}
322
323void multiBM::bayes ( const vec &yt, const vec &cond ) {
324        if ( frg < 1.0 ) {
325                beta *= frg;
326                last_lognc = est.lognc();
327        }
328        beta += yt;
329        if ( evalll ) {
330                ll = est.lognc() - last_lognc;
331        }
332}
333
334double multiBM::logpred ( const vec &yt ) const {
335        eDirich pred ( est );
336        vec &beta = pred._beta();
337
338        double lll;
339        if ( frg < 1.0 ) {
340                beta *= frg;
341                lll = pred.lognc();
342        } else if ( evalll ) {
343                lll = last_lognc;
344        } else {
345                lll = pred.lognc();
346        }
347
348        beta += yt;
349        return pred.lognc() - lll;
350}
351void multiBM::flatten ( const BMEF* B ) {
352        const multiBM* E = dynamic_cast<const multiBM*> ( B );
353        // sum(beta) should be equal to sum(B.beta)
354        const vec &Eb = E->beta;//const_cast<multiBM*> ( E )->_beta();
355        beta *= ( sum ( Eb ) / sum ( beta ) );
356        if ( evalll ) {
357                last_lognc = est.lognc();
358        }
359}
360
361vec egamma::sample() const {
362        vec smp ( dim );
363        int i;
364
365        for ( i = 0; i < dim; i++ ) {
366                if ( beta ( i ) > std::numeric_limits<double>::epsilon() ) {
367                        GamRNG.setup ( alpha ( i ), beta ( i ) );
368                } else {
369                        GamRNG.setup ( alpha ( i ), std::numeric_limits<double>::epsilon() );
370                }
371#pragma omp critical
372                smp ( i ) = GamRNG();
373        }
374
375        return smp;
376}
377
378// mat egamma::sample ( int N ) const {
379//      mat Smp ( rv.count(),N );
380//      int i,j;
381//
382//      for ( i=0; i<rv.count(); i++ ) {
383//              GamRNG.setup ( alpha ( i ),beta ( i ) );
384//
385//              for ( j=0; j<N; j++ ) {
386//                      Smp ( i,j ) = GamRNG();
387//              }
388//      }
389//
390//      return Smp;
391// }
392
393double egamma::evallog ( const vec &val ) const {
394        double res = 0.0; //the rest will be added
395        int i;
396
397        if ( any ( val <= 0. ) ) return -inf;
398        if ( any ( beta <= 0. ) ) return -inf;
399        for ( i = 0; i < dim; i++ ) {
400                res += ( alpha ( i ) - 1 ) * std::log ( val ( i ) ) - beta ( i ) * val ( i );
401        }
402        double tmp = res - lognc();;
403        bdm_assert_debug ( std::isfinite ( tmp ), "Infinite value" );
404        return tmp;
405}
406
407double egamma::lognc() const {
408        double res = 0.0; //will be added
409        int i;
410
411        for ( i = 0; i < dim; i++ ) {
412                res += lgamma ( alpha ( i ) ) - alpha ( i ) * std::log ( beta ( i ) ) ;
413        }
414
415        return res;
416}
417
418void mgamma::set_parameters ( double k0, const vec &beta0 ) {
419        k = k0;
420        iepdf.set_parameters ( k * ones ( beta0.length() ), beta0 );
421        dimc = iepdf.dimension();
422        dim = iepdf.dimension();
423}
424
425void eEmp::resample ( ivec &ind, RESAMPLING_METHOD method ) {
426        ind = zeros_i ( n );
427        ivec N_babies = zeros_i ( n );
428        vec cumDist = cumsum ( w );
429        vec u ( n );
430        int i, j, parent;
431        double u0;
432
433        switch ( method ) {
434        case MULTINOMIAL:
435                u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n );
436
437                for ( i = n - 2; i >= 0; i-- ) {
438                        u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) );
439                }
440
441                break;
442
443        case STRATIFIED:
444
445                for ( i = 0; i < n; i++ ) {
446                        u ( i ) = ( i + UniRNG.sample() ) / n;
447                }
448
449                break;
450
451        case SYSTEMATIC:
452                u0 = UniRNG.sample();
453
454                for ( i = 0; i < n; i++ ) {
455                        u ( i ) = ( i + u0 ) / n;
456                }
457
458                break;
459
460        default:
461                bdm_error ( "PF::resample(): Unknown resampling method" );
462        }
463
464        // U is now full
465        j = 0;
466
467        for ( i = 0; i < n; i++ ) {
468                while ( u ( i ) > cumDist ( j ) ) j++;
469
470                N_babies ( j ) ++;
471        }
472        // We have assigned new babies for each Particle
473        // Now, we fill the resulting index such that:
474        // * particles with at least one baby should not move *
475        // This assures that reassignment can be done inplace;
476
477        // find the first parent;
478        parent = 0;
479        while ( N_babies ( parent ) == 0 ) parent++;
480
481        // Build index
482        for ( i = 0; i < n; i++ ) {
483                if ( N_babies ( i ) > 0 ) {
484                        ind ( i ) = i;
485                        N_babies ( i ) --; //this index was now replicated;
486                } else {
487                        // test if the parent has been fully replicated
488                        // if yes, find the next one
489                        while ( ( N_babies ( parent ) == 0 ) || ( N_babies ( parent ) == 1 && parent > i ) ) parent++;
490
491                        // Replicate parent
492                        ind ( i ) = parent;
493
494                        N_babies ( parent ) --; //this index was now replicated;
495                }
496
497        }
498
499// copy the internals according to ind
500        for ( i = 0; i < n; i++ ) {
501                if ( ind ( i ) != i ) {
502                        samples ( i ) = samples ( ind ( i ) );
503                }
504                w ( i ) = 1.0 / n;
505        }
506}
507
508void eEmp::set_statistics ( const vec &w0, const epdf &epdf0 ) {
509        dim = epdf0.dimension();
510        w = w0;
511        w /= sum ( w0 );//renormalize
512        n = w.length();
513        samples.set_size ( n );
514
515        for ( int i = 0; i < n; i++ ) {
516                samples ( i ) = epdf0.sample();
517        }
518}
519
520void eEmp::set_samples ( const epdf* epdf0 ) {
521        w = 1;
522        w /= sum ( w );//renormalize
523
524        for ( int i = 0; i < n; i++ ) {
525                samples ( i ) = epdf0->sample();
526        }
527}
528
529void migamma_ref::from_setting ( const Setting &set ) {
530        vec ref;
531        UI::get ( ref, set, "ref" , UI::compulsory );
532        set_parameters ( set["k"], ref, set["l"] );
533}
534
535void mlognorm::from_setting ( const Setting &set ) {
536        vec mu0;
537        UI::get ( mu0, set, "mu0", UI::compulsory );
538        set_parameters ( mu0.length(), set["k"] );
539        condition ( mu0 );
540}
541
542void mlstudent::condition ( const vec &cond ) {
543        if ( cond.length() > 0 ) {
544                iepdf._mu() = A * cond + mu_const;
545        } else {
546                iepdf._mu() =  mu_const;
547        }
548        double zeta;
549        //ugly hack!
550        if ( ( cond.length() + 1 ) == Lambda.rows() ) {
551                zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) );
552        } else {
553                zeta = Lambda.invqform ( cond );
554        }
555        _R = Re;
556        _R *= ( 1 + zeta );// / ( nu ); << nu is in Re!!!!!!
557}
558
559void eEmp::qbounds ( vec &lb, vec &ub, double perc ) const {
560        // lb in inf so than it will be pushed below;
561        lb.set_size ( dim );
562        ub.set_size ( dim );
563        lb = std::numeric_limits<double>::infinity();
564        ub = -std::numeric_limits<double>::infinity();
565        int j;
566        for ( int i = 0; i < n; i++ ) {
567                for ( j = 0; j < dim; j++ ) {
568                        if ( samples ( i ) ( j ) < lb ( j ) ) {
569                                lb ( j ) = samples ( i ) ( j );
570                        }
571                        if ( samples ( i ) ( j ) > ub ( j ) ) {
572                                ub ( j ) = samples ( i ) ( j );
573                        }
574                }
575        }
576}
577
578
579};
Note: See TracBrowser for help on using the browser.