root/bdm/stat/libEF.h @ 215

Revision 215, 20.1 kB (checked in by smidl, 15 years ago)

PMSM new sim.cpp

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