root/bdm/stat/libEF.h @ 193

Revision 193, 17.8 kB (checked in by smidl, 16 years ago)

oprava merger_iter a jeho casti

  • Property svn:eol-style set to native
Line 
1/*!
2  \file
3  \brief Probability distributions for Exponential Family models.
4  \author Vaclav Smidl.
5
6  -----------------------------------
7  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
8
9  Using IT++ for numerical operations
10  -----------------------------------
11*/
12
13#ifndef EF_H
14#define EF_H
15
16#include <itpp/itbase.h>
17#include "../math/libDC.h"
18#include "libBM.h"
19#include "../itpp_ext.h"
20//#include <std>
21
22using namespace itpp;
23
24
25//! Global Uniform_RNG
26extern Uniform_RNG UniRNG;
27//! Global Normal_RNG
28extern Normal_RNG NorRNG;
29//! Global Gamma_RNG
30extern Gamma_RNG GamRNG;
31
32/*!
33* \brief General conjugate exponential family posterior density.
34
35* More?...
36*/
37
38class eEF : public epdf {
39public:
40//      eEF() :epdf() {};
41        //! default constructor
42        eEF ( const RV &rv ) :epdf ( rv ) {};
43        //! logarithm of the normalizing constant, \f$\mathcal{I}\f$
44        virtual double lognc() const =0;
45        //!TODO decide if it is really needed
46        virtual void dupdate ( mat &v ) {it_error ( "Not implemented" );};
47        //!Evaluate normalized log-probability
48        virtual double evalpdflog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
49        //!Evaluate normalized log-probability
50        virtual double evalpdflog ( const vec &val ) const {return evalpdflog_nn ( val )-lognc();}
51        //!Evaluate normalized log-probability for many samples
52        virtual vec evalpdflog ( const mat &Val ) const {
53                vec x ( Val.cols() );
54                for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evalpdflog_nn ( Val.get_col ( i ) ) ;}
55                return x-lognc();
56        }
57        //!Power of the density, used e.g. to flatten the density
58        virtual void pow ( double p ) {it_error ( "Not implemented" );};
59};
60
61/*!
62* \brief Exponential family model.
63
64* More?...
65*/
66
67class mEF : public mpdf {
68
69public:
70        //! Default constructor
71        mEF ( const RV &rv0, const RV &rvc0 ) :mpdf ( rv0,rvc0 ) {};
72};
73
74//! Estimator for Exponential family
75class BMEF : public BM {
76protected:
77        //! forgetting factor
78        double frg;
79        //! cached value of lognc() in the previous step (used in evaluation of \c ll )
80        double last_lognc;
81public:
82        //! Default constructor
83        BMEF ( const RV &rv, double frg0=1.0 ) :BM ( rv ), frg ( frg0 ) {}
84        //! Copy constructor
85        BMEF ( const BMEF &B ) :BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
86        //!get statistics from another model
87        virtual void set_statistics ( const BMEF* BM0 ) {it_error ( "Not implemented" );};
88        //! Weighted update of sufficient statistics (Bayes rule)
89        virtual void bayes ( const vec &data, const double w ) {};
90        //original Bayes
91        void bayes ( const vec &dt );
92        //!Flatten the posterior according to the given BMEF (of the same type!)
93        virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
94        //!Flatten the posterior as if to keep nu0 data
95        virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
96};
97
98template<class sq_T>
99class mlnorm;
100
101/*!
102* \brief Gaussian density with positive definite (decomposed) covariance matrix.
103
104* More?...
105*/
106template<class sq_T>
107class enorm : public eEF {
108protected:
109        //! mean value
110        vec mu;
111        //! Covariance matrix in decomposed form
112        sq_T R;
113        //! dimension (redundant from rv.count() for easier coding )
114        int dim;
115public:
116        //!Default constructor
117        enorm ( const RV &rv );
118        //! Set mean value \c mu and covariance \c R
119        void set_parameters ( const vec &mu,const sq_T &R );
120        ////! tupdate in exponential form (not really handy)
121        //void tupdate ( double phi, mat &vbar, double nubar );
122        //! dupdate in exponential form (not really handy)
123        void dupdate ( mat &v,double nu=1.0 );
124
125        vec sample() const;
126        //! TODO is it used?
127        mat sample ( int N ) const;
128        double eval ( const vec &val ) const ;
129        double evalpdflog_nn ( const vec &val ) const;
130        double lognc () const;
131        vec mean() const {return mu;}
132//      mlnorm<sq_T>* condition ( const RV &rvn ) const ;
133        mpdf* condition ( const RV &rvn ) const ;
134//      enorm<sq_T>* marginal ( const RV &rv ) const;
135        epdf* marginal ( const RV &rv ) const;
136//Access methods
137        //! returns a pointer to the internal mean value. Use with Care!
138        vec& _mu() {return mu;}
139
140        //! access function
141        void set_mu ( const vec mu0 ) { mu=mu0;}
142
143        //! returns pointers to the internal variance and its inverse. Use with Care!
144        sq_T& _R() {return R;}
145
146        //! access method
147//      mat getR () {return R.to_mat();}
148};
149
150/*!
151* \brief Gauss-inverse-Wishart density stored in LD form
152
153* For \f$p\f$-variate densities, given rv.count() should be \f$p\times\f$ V.rows().
154*
155*/
156class egiw : public eEF {
157protected:
158        //! Extended information matrix of sufficient statistics
159        ldmat V;
160        //! Number of data records (degrees of freedom) of sufficient statistics
161        double nu;
162        //! Dimension of the output
163        int xdim;
164        //! Dimension of the regressor
165        int nPsi;
166public:
167        //!Default constructor, assuming
168        egiw ( RV rv, mat V0, double nu0 ) : eEF ( rv ), V ( V0 ), nu ( nu0 ) {
169                xdim = rv.count() /V.rows();
170                it_assert_debug ( rv.count() ==xdim*V.rows(),"Incompatible V0." );
171                nPsi = V.rows()-xdim;
172        }
173        //!Full constructor for V in ldmat form
174        egiw ( RV rv, ldmat V0, double nu0 ) : eEF ( rv ), V ( V0 ), nu ( nu0 ) {
175                xdim = rv.count() /V.rows();
176                it_assert_debug ( rv.count() ==xdim*V.rows(),"Incompatible V0." );
177                nPsi = V.rows()-xdim;
178        }
179
180        vec sample() const;
181        vec mean() const;
182        void mean_mat ( mat &M, mat&R ) const;
183        //! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
184        double evalpdflog_nn ( const vec &val ) const;
185        double lognc () const;
186
187        //Access
188        //! returns a pointer to the internal statistics. Use with Care!
189        ldmat& _V() {return V;}
190        //! returns a pointer to the internal statistics. Use with Care!
191        double& _nu() {return nu;}
192        void pow ( double p );
193};
194
195/*! \brief Dirichlet posterior density
196
197Continuous Dirichlet density of \f$n\f$-dimensional variable \f$x\f$
198\f[
199f(x|\beta) = \frac{\Gamma[\gamma]}{\prod_{i=1}^{n}\Gamma(\beta_i)} \prod_{i=1}^{n}x_i^{\beta_i-1}
200\f]
201where \f$\gamma=\sum_i \beta_i\f$.
202*/
203class eDirich: public eEF {
204protected:
205        //!sufficient statistics
206        vec beta;
207public:
208        //!Default constructor
209        eDirich ( const RV &rv, const vec &beta0 ) : eEF ( rv ),beta ( beta0 ) {it_assert_debug ( rv.count() ==beta.length(),"Incompatible statistics" ); };
210        //! Copy constructor
211        eDirich ( const eDirich &D0 ) : eEF ( D0.rv ),beta ( D0.beta ) {};
212        vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
213        vec mean() const {return beta/sum ( beta );};
214        //! In this instance, val is ...
215        double evalpdflog_nn ( const vec &val ) const {return ( beta-1 ) *log ( val );};
216        double lognc () const {
217                double gam=sum ( beta );
218                double lgb=0.0;
219                for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
220                return lgb-lgamma ( gam );
221        };
222        //!access function
223        vec& _beta()  {return beta;}
224        //!Set internal parameters
225        void set_parameters ( const vec &beta0 ) {
226                if ( beta0.length() !=beta.length() ) {
227                        it_assert_debug ( rv.length() ==1,"Undefined" );
228                        rv.set_size ( 0,beta0.length() );
229                }
230                beta= beta0;
231        }
232};
233
234//! \brief Estimator for Multinomial density
235class multiBM : public BMEF {
236protected:
237        //! Conjugate prior and posterior
238        eDirich est;
239        //! Pointer inside est to sufficient statistics
240        vec &beta;
241public:
242        //!Default constructor
243        multiBM ( const RV &rv, const vec beta0 ) : BMEF ( rv ),est ( rv,beta0 ),beta ( est._beta() ) {last_lognc=est.lognc();}
244        //!Copy constructor
245        multiBM ( const multiBM &B ) : BMEF ( B ),est ( rv,B.beta ),beta ( est._beta() ) {}
246        //! Sets sufficient statistics to match that of givefrom mB0
247        void set_statistics ( const BM* mB0 ) {const multiBM* mB=dynamic_cast<const multiBM*> ( mB0 ); beta=mB->beta;}
248        void bayes ( const vec &dt ) {
249                if ( frg<1.0 ) {beta*=frg;last_lognc=est.lognc();}
250                beta+=dt;
251                if ( evalll ) {ll=est.lognc()-last_lognc;}
252        }
253        double logpred ( const vec &dt ) const {
254                eDirich pred ( est );
255                vec &beta = pred._beta();
256
257                double lll;
258                if ( frg<1.0 )
259                        {beta*=frg;lll=pred.lognc();}
260                else
261                        if ( evalll ) {lll=last_lognc;}
262                        else{lll=pred.lognc();}
263
264                beta+=dt;
265                return pred.lognc()-lll;
266        }
267        void flatten ( const BMEF* B ) {
268                const eDirich* E=dynamic_cast<const eDirich*> ( B );
269                // sum(beta) should be equal to sum(B.beta)
270                const vec &Eb=const_cast<eDirich*> ( E )->_beta();
271                est.pow ( sum ( beta ) /sum ( Eb ) );
272                if ( evalll ) {last_lognc=est.lognc();}
273        }
274        const epdf& _epdf() const {return est;};
275        void set_parameters ( const vec &beta0 ) {
276                est.set_parameters ( beta0 );
277                rv = est._rv();
278                if ( evalll ) {last_lognc=est.lognc();}
279        }
280};
281
282/*!
283 \brief Gamma posterior density
284
285 Multivariate Gamma density as product of independent univariate densities.
286 \f[
287 f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i)
288 \f]
289*/
290
291class egamma : public eEF {
292protected:
293        //! Vector \f$\alpha\f$
294        vec alpha;
295        //! Vector \f$\beta\f$
296        vec beta;
297public :
298        //! Default constructor
299        egamma ( const RV &rv ) :eEF ( rv ) {};
300        //! Sets parameters
301        void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;};
302        vec sample() const;
303        //! TODO: is it used anywhere?
304//      mat sample ( int N ) const;
305        double evalpdflog ( const vec &val ) const;
306        double lognc () const;
307        //! Returns poiter to alpha and beta. Potentially dengerous: use with care!
308        void _param ( vec* &a, vec* &b ) {a=&alpha;b=&beta;};
309        vec mean() const {vec pom ( alpha ); pom/=beta; return pom;}
310};
311/*
312//! Weighted mixture of epdfs with external owned components.
313class emix : public epdf {
314protected:
315        int n;
316        vec &w;
317        Array<epdf*> Coms;
318public:
319//! Default constructor
320        emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
321        void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
322        vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
323        vec sample() {it_error ( "Not implemented" );return 0;}
324};
325*/
326
327//!  Uniform distributed density on a rectangular support
328
329class euni: public epdf {
330protected:
331//! lower bound on support
332        vec low;
333//! upper bound on support
334        vec high;
335//! internal
336        vec distance;
337//! normalizing coefficients
338        double nk;
339//! cache of log( \c nk )
340        double lnk;
341public:
342        //! Defualt constructor
343        euni ( const RV rv ) :epdf ( rv ) {}
344        double eval ( const vec &val ) const  {return nk;}
345        double evalpdflog ( const vec &val ) const  {return lnk;}
346        vec sample() const {
347                vec smp ( rv.count() );
348#pragma omp critical
349                UniRNG.sample_vector ( rv.count(),smp );
350                return low+elem_mult ( distance,smp );
351        }
352        //! set values of \c low and \c high
353        void set_parameters ( const vec &low0, const vec &high0 ) {
354                distance = high0-low0;
355                it_assert_debug ( min ( distance ) >0.0,"bad support" );
356                low = low0;
357                high = high0;
358                nk = prod ( 1.0/distance );
359                lnk = log ( nk );
360        }
361        vec mean() const {vec pom=high; pom-=low; pom/=2.0; return pom;}
362};
363
364
365/*!
366 \brief Normal distributed linear function with linear function of mean value;
367
368 Mean value \f$mu=A*rvc+mu_0\f$.
369*/
370template<class sq_T>
371class mlnorm : public mEF {
372        //! Internal epdf that arise by conditioning on \c rvc
373        enorm<sq_T> epdf;
374        mat A;
375        vec mu_const;
376        vec& _mu; //cached epdf.mu;
377public:
378        //! Constructor
379        mlnorm (const RV &rv, const RV &rvc );
380        //! Set \c A and \c R
381        void set_parameters ( const  mat &A, const vec &mu0, const sq_T &R );
382//      //!Generate one sample of the posterior
383//      vec samplecond (const vec &cond, double &lik );
384//      //!Generate matrix of samples of the posterior
385//      mat samplecond (const vec &cond, vec &lik, int n );
386        //! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
387        void condition (const vec &cond );
388       
389        template<class sq_M>
390        friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M> &ml );
391};
392
393/*!
394\brief  Gamma random walk
395
396Mean value, \f$\mu\f$, of this density is given by \c rvc .
397Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
398This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
399
400The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
401*/
402class mgamma : public mEF {
403protected:
404        //! Internal epdf that arise by conditioning on \c rvc
405        egamma epdf;
406        //! Constant \f$k\f$
407        double k;
408        //! cache of epdf.beta
409        vec* _beta;
410
411public:
412        //! Constructor
413        mgamma ( const RV &rv,const RV &rvc );
414        //! Set value of \c k
415        void set_parameters ( double k );
416        void condition ( const vec &val ) {*_beta=k/val;};
417};
418
419/*!
420\brief  Gamma random walk around a fixed point
421
422Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
423\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
424
425Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
426This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
427
428The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
429*/
430class mgamma_fix : public mgamma {
431protected:
432        //! parameter l
433        double l;
434        //! reference vector
435        vec refl;
436public:
437        //! Constructor
438        mgamma_fix ( const RV &rv,const RV &rvc ) : mgamma ( rv,rvc ),refl ( rv.count() ) {};
439        //! Set value of \c k
440        void set_parameters ( double k0 , vec ref0, double l0 ) {
441                mgamma::set_parameters ( k0 );
442                refl=pow ( ref0,1.0-l0 );l=l0;
443        };
444
445        void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); *_beta=k/mean;};
446};
447
448//! Switch between various resampling methods.
449enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
450/*!
451\brief Weighted empirical density
452
453Used e.g. in particle filters.
454*/
455class eEmp: public epdf {
456protected :
457        //! Number of particles
458        int n;
459        //! Sample weights \f$w\f$
460        vec w;
461        //! Samples \f$x^{(i)}, i=1..n\f$
462        Array<vec> samples;
463public:
464        //! Default constructor
465        eEmp ( const RV &rv0 ,int n0 ) :epdf ( rv0 ),n ( n0 ),w ( n ),samples ( n ) {};
466        //! Set samples and weights
467        void set_parameters ( const vec &w0, const epdf* pdf0 );
468        //! Set sample
469        void set_samples ( const epdf* pdf0 );
470        //! Potentially dangerous, use with care.
471        vec& _w()  {return w;};
472        //! access function
473        Array<vec>& _samples() {return samples;};
474        //! Function performs resampling, i.e. removal of low-weight samples and duplication of high-weight samples such that the new samples represent the same density.
475        ivec resample ( RESAMPLING_METHOD method = SYSTEMATIC );
476        //! inherited operation : NOT implemneted
477        vec sample() const {it_error ( "Not implemented" );return 0;}
478        //! inherited operation : NOT implemneted
479        double evalpdflog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
480        vec mean() const {
481                vec pom=zeros ( rv.count() );
482                for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
483                return pom;
484        }
485};
486
487
488////////////////////////
489
490template<class sq_T>
491enorm<sq_T>::enorm ( const RV &rv ) :eEF ( rv ), mu ( rv.count() ),R ( rv.count() ),dim ( rv.count() ) {};
492
493template<class sq_T>
494void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) {
495//Fixme test dimensions of mu0 and R0;
496        mu = mu0;
497        R = R0;
498};
499
500template<class sq_T>
501void enorm<sq_T>::dupdate ( mat &v, double nu ) {
502        //
503};
504
505// template<class sq_T>
506// void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
507//      //
508// };
509
510template<class sq_T>
511vec enorm<sq_T>::sample() const {
512        vec x ( dim );
513        NorRNG.sample_vector ( dim,x );
514        vec smp = R.sqrt_mult ( x );
515
516        smp += mu;
517        return smp;
518};
519
520template<class sq_T>
521mat enorm<sq_T>::sample ( int N ) const {
522        mat X ( dim,N );
523        vec x ( dim );
524        vec pom;
525        int i;
526
527        for ( i=0;i<N;i++ ) {
528                NorRNG.sample_vector ( dim,x );
529                pom = R.sqrt_mult ( x );
530                pom +=mu;
531                X.set_col ( i, pom );
532        }
533
534        return X;
535};
536
537template<class sq_T>
538double enorm<sq_T>::eval ( const vec &val ) const {
539        double pdfl,e;
540        pdfl = evalpdflog ( val );
541        e = exp ( pdfl );
542        return e;
543};
544
545template<class sq_T>
546double enorm<sq_T>::evalpdflog_nn ( const vec &val ) const {
547        // 1.83787706640935 = log(2pi)
548        return  -0.5* ( R.invqform ( mu-val ) );// - lognc();
549};
550
551template<class sq_T>
552inline double enorm<sq_T>::lognc () const {
553        // 1.83787706640935 = log(2pi)
554        return 0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
555};
556
557template<class sq_T>
558mlnorm<sq_T>::mlnorm ( const RV &rv0, const RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu ( epdf._mu() ) {
559        ep =&epdf;
560}
561
562template<class sq_T>
563void mlnorm<sq_T>::set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 ) {
564        epdf.set_parameters ( zeros ( rv.count() ),R0 );
565        A = A0;
566        mu_const = mu0;
567}
568
569// template<class sq_T>
570// vec mlnorm<sq_T>::samplecond (const  vec &cond, double &lik ) {
571//      this->condition ( cond );
572//      vec smp = epdf.sample();
573//      lik = epdf.eval ( smp );
574//      return smp;
575// }
576
577// template<class sq_T>
578// mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) {
579//      int i;
580//      int dim = rv.count();
581//      mat Smp ( dim,n );
582//      vec smp ( dim );
583//      this->condition ( cond );
584//
585//      for ( i=0; i<n; i++ ) {
586//              smp = epdf.sample();
587//              lik ( i ) = epdf.eval ( smp );
588//              Smp.set_col ( i ,smp );
589//      }
590//
591//      return Smp;
592// }
593
594template<class sq_T>
595void mlnorm<sq_T>::condition (const vec &cond ) {
596        _mu = A*cond + mu_const;
597//R is already assigned;
598}
599
600template<class sq_T>
601epdf* enorm<sq_T>::marginal ( const RV &rvn ) const {
602        ivec irvn = rvn.dataind ( rv );
603
604        sq_T Rn ( R,irvn );
605        enorm<sq_T>* tmp = new enorm<sq_T>( rvn );
606        tmp->set_parameters ( mu ( irvn ), Rn );
607        return tmp;
608}
609
610template<class sq_T>
611mpdf* enorm<sq_T>::condition ( const RV &rvn ) const {
612
613        RV rvc = rv.subt ( rvn );
614        it_assert_debug ( ( rvc.count() +rvn.count() ==rv.count() ),"wrong rvn" );
615        //Permutation vector of the new R
616        ivec irvn = rvn.dataind ( rv );
617        ivec irvc = rvc.dataind ( rv );
618        ivec perm=concat ( irvn , irvc );
619        sq_T Rn ( R,perm );
620
621        //fixme - could this be done in general for all sq_T?
622        mat S=Rn.to_mat();
623        //fixme
624        int n=rvn.count()-1;
625        int end=R.rows()-1;
626        mat S11 = S.get ( 0,n, 0, n );
627        mat S12 = S.get (0, n , rvn.count(), end );
628        mat S22 = S.get ( rvn.count(), end, rvn.count(), end );
629
630        vec mu1 = mu ( irvn );
631        vec mu2 = mu ( irvc );
632        mat A=S12*inv ( S22 );
633        sq_T R_n ( S11 - A *S12.T() );
634
635        mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( rvn,rvc );
636
637        tmp->set_parameters ( A,mu1-A*mu2,R_n );
638        return tmp;
639}
640
641///////////
642
643template<class sq_T>
644std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_T> &ml ){
645        os << "A:"<< ml.A<<endl;
646        os << "mu:"<< ml.mu_const<<endl;
647        os << "R:" << ml.epdf._R().to_mat() <<endl;
648        return os;
649};
650
651#endif //EF_H
Note: See TracBrowser for help on using the browser.