root/bdm/stat/libEF.h @ 173

Revision 173, 15.4 kB (checked in by smidl, 16 years ago)

drobnosti

  • 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 implemneted" );};
47        //!Evaluate normalized log-probability
48        virtual double evalpdflog_nn ( const vec &val ) const{it_error ( "Not implemneted" );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
93        virtual void flatten ( BMEF * B) {it_error ( "Not implemented" );}
94};
95
96/*!
97* \brief Gaussian density with positive definite (decomposed) covariance matrix.
98
99* More?...
100*/
101template<class sq_T>
102
103class enorm : public eEF {
104protected:
105        //! mean value
106        vec mu;
107        //! Covariance matrix in decomposed form
108        sq_T R;
109        //! dimension (redundant from rv.count() for easier coding )
110        int dim;
111public:
112//      enorm() :eEF() {};
113        //!Default constructor
114        enorm ( RV &rv );
115        //! Set mean value \c mu and covariance \c R
116        void set_parameters ( const vec &mu,const sq_T &R );
117        //! tupdate in exponential form (not really handy)
118        void tupdate ( double phi, mat &vbar, double nubar );
119        //! dupdate in exponential form (not really handy)
120        void dupdate ( mat &v,double nu=1.0 );
121
122        vec sample() const;
123        //! TODO is it used?
124        mat sample ( int N ) const;
125        double eval ( const vec &val ) const ;
126        double evalpdflog ( const vec &val ) const;
127        double lognc () const;
128        vec mean() const {return mu;}
129
130//Access methods
131        //! returns a pointer to the internal mean value. Use with Care!
132        vec& _mu() {return mu;}
133
134        //! access function
135        void set_mu ( const vec mu0 ) { mu=mu0;}
136
137        //! returns pointers to the internal variance and its inverse. Use with Care!
138        sq_T& _R() {return R;}
139
140        //! access method
141        mat getR () {return R.to_mat();}
142};
143
144/*!
145* \brief Gauss-inverse-Wishart density stored in LD form
146
147* For \f$p\f$-variate densities, given rv.count() should be \f$p\times\f$ V.rows().
148*
149*/
150class egiw : public eEF {
151protected:
152        //! Extended information matrix of sufficient statistics
153        ldmat V;
154        //! Number of data records (degrees of freedom) of sufficient statistics
155        double nu;
156        //! Dimension of the output
157        int xdim;
158        //! Dimension of the regressor
159        int nPsi;
160public:
161        //!Default constructor, assuming
162        egiw ( RV rv, mat V0, double nu0 ) : eEF ( rv ), V ( V0 ), nu ( nu0 ) {
163                xdim = rv.count() /V.rows();
164                it_assert_debug ( rv.count() ==xdim*V.rows(),"Incompatible V0." );
165                nPsi = V.rows()-xdim;
166        }
167        //!Full constructor for V in ldmat form
168        egiw ( RV rv, ldmat 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
174        vec sample() const;
175        vec mean() const;
176        void mean_mat ( mat &M, mat&R ) const;
177        //! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
178        double evalpdflog_nn ( const vec &val ) const;
179        double lognc () const;
180
181        //Access
182        //! returns a pointer to the internal statistics. Use with Care!
183        ldmat& _V() {return V;}
184        //! returns a pointer to the internal statistics. Use with Care!
185        double& _nu() {return nu;}
186        void pow ( double p );
187};
188
189/*! \brief Dirichlet posterior density
190
191Continuous Dirichlet density of \f$n\f$-dimensional variable \f$x\f$
192\f[
193f(x|\beta) = \frac{\Gamma[\gamma]}{\prod_{i=1}^{n}\Gamma(\beta_i)} \prod_{i=1}^{n}x_i^{\beta_i-1}
194\f]
195where \f$\gamma=\sum_i \beta_i\f$.
196*/
197class eDirich: public eEF {
198protected:
199        //!sufficient statistics
200        vec beta;
201public:
202        //!Default constructor
203        eDirich ( const RV &rv, const vec &beta0 ) : eEF ( rv ),beta ( beta0 ) {it_assert_debug ( rv.count() ==beta.length(),"Incompatible statistics" ); };
204        //! Copy constructor
205        eDirich ( const eDirich &D0 ) : eEF ( D0.rv ),beta ( D0.beta ) {};
206        vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
207        vec mean() const {return beta/sum ( beta );};
208        //! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
209        double evalpdflog_nn ( const vec &val ) const {return ( beta-1 ) *log ( val );};
210        double lognc () const {
211                double gam=sum ( beta );
212                double lgb=0.0;
213                for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
214                return lgb-lgamma ( gam );
215        };
216        //!access function
217        vec& _beta() {return beta;}
218};
219
220//! Estimator for Multinomial density
221class multiBM : public BMEF {
222protected:
223        //! Conjugate prior and posterior
224        eDirich est;
225        vec &beta;
226public:
227        //!Default constructor
228        multiBM ( const RV &rv, const vec beta0 ) : BMEF ( rv ),est ( rv,beta0 ),beta ( est._beta() ) {last_lognc=est.lognc();}
229        //!Copy constructor
230        multiBM ( const multiBM &B ) : BMEF ( B ),est ( rv,B.beta ),beta ( est._beta() ) {}
231
232        void set_statistics ( const BM* mB0 ) {const multiBM* mB=dynamic_cast<const multiBM*> ( mB0 ); beta=mB->beta;}
233        void bayes ( const vec &dt ) {
234                if ( frg<1.0 ) {beta*=frg;last_lognc=est.lognc();}
235                beta+=dt;
236                if ( evalll ) {ll=est.lognc()-last_lognc;}
237        }
238        double logpred ( const vec &dt ) const {
239                eDirich pred ( est );
240                vec &beta = pred._beta();
241               
242                double lll;
243                if ( frg<1.0 )
244                        {beta*=frg;lll=pred.lognc();}
245                else
246                        if ( evalll ) {lll=last_lognc;}
247                        else{lll=pred.lognc();}
248
249                beta+=dt;
250                return pred.lognc()-lll;
251        }
252        void flatten (BMEF* B ) {
253                eDirich* E=dynamic_cast<eDirich*>(B);
254                // sum(beta) should be equal to sum(B.beta)
255                const vec &Eb=E->_beta();
256                est.pow ( sum(beta)/sum(Eb) );
257                if(evalll){last_lognc=est.lognc();}
258        }
259        const epdf& _epdf() const {return est;};
260        //!access funct
261};
262
263/*!
264 \brief Gamma posterior density
265
266 Multivariate Gamma density as product of independent univariate densities.
267 \f[
268 f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i)
269 \f]
270*/
271
272class egamma : public eEF {
273protected:
274        //! Vector \f$\alpha\f$
275        vec alpha;
276        //! Vector \f$\beta\f$
277        vec beta;
278public :
279        //! Default constructor
280        egamma ( const RV &rv ) :eEF ( rv ) {};
281        //! Sets parameters
282        void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;};
283        vec sample() const;
284        //! TODO: is it used anywhere?
285//      mat sample ( int N ) const;
286        double evalpdflog ( const vec &val ) const;
287        double lognc () const;
288        //! Returns poiter to alpha and beta. Potentially dengerous: use with care!
289        void _param ( vec* &a, vec* &b ) {a=&alpha;b=&beta;};
290        vec mean() const {vec pom ( alpha ); pom/=beta; return pom;}
291};
292/*
293//! Weighted mixture of epdfs with external owned components.
294class emix : public epdf {
295protected:
296        int n;
297        vec &w;
298        Array<epdf*> Coms;
299public:
300//! Default constructor
301        emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
302        void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
303        vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
304        vec sample() {it_error ( "Not implemented" );return 0;}
305};
306*/
307
308//!  Uniform distributed density on a rectangular support
309
310class euni: public epdf {
311protected:
312//! lower bound on support
313        vec low;
314//! upper bound on support
315        vec high;
316//! internal
317        vec distance;
318//! normalizing coefficients
319        double nk;
320//! cache of log( \c nk )
321        double lnk;
322public:
323        //! Defualt constructor
324        euni ( const RV rv ) :epdf ( rv ) {}
325        double eval ( const vec &val ) const  {return nk;}
326        double evalpdflog ( const vec &val ) const  {return lnk;}
327        vec sample() const {
328                vec smp ( rv.count() );
329#pragma omp critical
330                UniRNG.sample_vector ( rv.count(),smp );
331                return low+elem_mult ( distance,smp );
332        }
333        //! set values of \c low and \c high
334        void set_parameters ( const vec &low0, const vec &high0 ) {
335                distance = high0-low0;
336                it_assert_debug ( min ( distance ) >0.0,"bad support" );
337                low = low0;
338                high = high0;
339                nk = prod ( 1.0/distance );
340                lnk = log ( nk );
341        }
342        vec mean() const {vec pom=high; pom-=low; pom/=2.0; return pom;}
343};
344
345
346/*!
347 \brief Normal distributed linear function with linear function of mean value;
348
349 Mean value \f$mu=A*rvc\f$.
350*/
351template<class sq_T>
352class mlnorm : public mEF {
353        //! Internal epdf that arise by conditioning on \c rvc
354        enorm<sq_T> epdf;
355        mat A;
356        vec& _mu; //cached epdf.mu;
357public:
358        //! Constructor
359        mlnorm ( RV &rv,RV &rvc );
360        //! Set \c A and \c R
361        void set_parameters ( const  mat &A, const sq_T &R );
362        //!Generate one sample of the posterior
363        vec samplecond ( vec &cond, double &lik );
364        //!Generate matrix of samples of the posterior
365        mat samplecond ( vec &cond, vec &lik, int n );
366        //! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
367        void condition ( vec &cond );
368};
369
370/*!
371\brief  Gamma random walk
372
373Mean value, \f$\mu\f$, of this density is given by \c rvc .
374Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
375This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
376
377The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
378*/
379class mgamma : public mEF {
380protected:
381        //! Internal epdf that arise by conditioning on \c rvc
382        egamma epdf;
383        //! Constant \f$k\f$
384        double k;
385        //! cache of epdf.beta
386        vec* _beta;
387
388public:
389        //! Constructor
390        mgamma ( const RV &rv,const RV &rvc );
391        //! Set value of \c k
392        void set_parameters ( double k );
393        void condition ( const vec &val ) {*_beta=k/val;};
394};
395
396/*!
397\brief  Gamma random walk around a fixed point
398
399Mean 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
400\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
401
402Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
403This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
404
405The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
406*/
407class mgamma_fix : public mgamma {
408protected:
409        //! parameter l
410        double l;
411        //! reference vector
412        vec refl;
413public:
414        //! Constructor
415        mgamma_fix ( const RV &rv,const RV &rvc ) : mgamma ( rv,rvc ),refl ( rv.count() ) {};
416        //! Set value of \c k
417        void set_parameters ( double k0 , vec ref0, double l0 ) {
418                mgamma::set_parameters ( k0 );
419                refl=pow ( ref0,1.0-l0 );l=l0;
420        };
421
422        void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); *_beta=k/mean;};
423};
424
425//! Switch between various resampling methods.
426enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
427/*!
428\brief Weighted empirical density
429
430Used e.g. in particle filters.
431*/
432class eEmp: public epdf {
433protected :
434        //! Number of particles
435        int n;
436        //! Sample weights \f$w\f$
437        vec w;
438        //! Samples \f$x^{(i)}, i=1..n\f$
439        Array<vec> samples;
440public:
441        //! Default constructor
442        eEmp ( const RV &rv0 ,int n0 ) :epdf ( rv0 ),n ( n0 ),w ( n ),samples ( n ) {};
443        //! Set sample
444        void set_parameters ( const vec &w0, epdf* pdf0 );
445        //! Potentially dangerous, use with care.
446        vec& _w()  {return w;};
447        //! access function
448        Array<vec>& _samples() {return samples;};
449        //! 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.
450        ivec resample ( RESAMPLING_METHOD method = SYSTEMATIC );
451        //! inherited operation : NOT implemneted
452        vec sample() const {it_error ( "Not implemented" );return 0;}
453        //! inherited operation : NOT implemneted
454        double evalpdflog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
455        vec mean() const {
456                vec pom=zeros ( rv.count() );
457                for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
458                return pom;
459        }
460};
461
462
463////////////////////////
464
465template<class sq_T>
466enorm<sq_T>::enorm ( RV &rv ) :eEF ( rv ), mu ( rv.count() ),R ( rv.count() ),dim ( rv.count() ) {};
467
468template<class sq_T>
469void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) {
470//Fixme test dimensions of mu0 and R0;
471        mu = mu0;
472        R = R0;
473};
474
475template<class sq_T>
476void enorm<sq_T>::dupdate ( mat &v, double nu ) {
477        //
478};
479
480template<class sq_T>
481void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
482        //
483};
484
485template<class sq_T>
486vec enorm<sq_T>::sample() const {
487        vec x ( dim );
488        NorRNG.sample_vector ( dim,x );
489        vec smp = R.sqrt_mult ( x );
490
491        smp += mu;
492        return smp;
493};
494
495template<class sq_T>
496mat enorm<sq_T>::sample ( int N ) const {
497        mat X ( dim,N );
498        vec x ( dim );
499        vec pom;
500        int i;
501
502        for ( i=0;i<N;i++ ) {
503                NorRNG.sample_vector ( dim,x );
504                pom = R.sqrt_mult ( x );
505                pom +=mu;
506                X.set_col ( i, pom );
507        }
508
509        return X;
510};
511
512template<class sq_T>
513double enorm<sq_T>::eval ( const vec &val ) const {
514        double pdfl,e;
515        pdfl = evalpdflog ( val );
516        e = exp ( pdfl );
517        return e;
518};
519
520template<class sq_T>
521double enorm<sq_T>::evalpdflog ( const vec &val ) const {
522        // 1.83787706640935 = log(2pi)
523        return  -0.5* ( +R.invqform ( mu-val ) ) - lognc();
524};
525
526template<class sq_T>
527inline double enorm<sq_T>::lognc () const {
528        // 1.83787706640935 = log(2pi)
529        return -0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
530};
531
532template<class sq_T>
533mlnorm<sq_T>::mlnorm ( RV &rv0,RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu ( epdf._mu() ) {
534        ep =&epdf;
535}
536
537template<class sq_T>
538void mlnorm<sq_T>::set_parameters ( const mat &A0, const sq_T &R0 ) {
539        epdf.set_parameters ( zeros ( rv.count() ),R0 );
540        A = A0;
541}
542
543template<class sq_T>
544vec mlnorm<sq_T>::samplecond ( vec &cond, double &lik ) {
545        this->condition ( cond );
546        vec smp = epdf.sample();
547        lik = epdf.eval ( smp );
548        return smp;
549}
550
551template<class sq_T>
552mat mlnorm<sq_T>::samplecond ( vec &cond, vec &lik, int n ) {
553        int i;
554        int dim = rv.count();
555        mat Smp ( dim,n );
556        vec smp ( dim );
557        this->condition ( cond );
558
559        for ( i=0; i<n; i++ ) {
560                smp = epdf.sample();
561                lik ( i ) = epdf.eval ( smp );
562                Smp.set_col ( i ,smp );
563        }
564
565        return Smp;
566}
567
568template<class sq_T>
569void mlnorm<sq_T>::condition ( vec &cond ) {
570        _mu = A*cond;
571//R is already assigned;
572}
573
574///////////
575
576
577#endif //EF_H
Note: See TracBrowser for help on using the browser.