root/bdm/stat/libEF.h @ 178

Revision 178, 17.2 kB (checked in by smidl, 16 years ago)

Changes in basic structures! new methods

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