root/library/bdm/stat/exp_family.h @ 634

Revision 634, 36.9 kB (checked in by smidl, 15 years ago)

Dirichlet sampling with unit test

  • Property svn:eol-style set to native
RevLine 
[8]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
[262]16
[461]17#include "../shared_ptr.h"
[384]18#include "../base/bdmbase.h"
[262]19#include "../math/chmat.h"
[8]20
[294]21namespace bdm
22{
[8]23
[32]24
25//! Global Uniform_RNG
[488]26extern Uniform_RNG UniRNG;
[33]27//! Global Normal_RNG
[488]28extern Normal_RNG NorRNG;
[33]29//! Global Gamma_RNG
[488]30extern Gamma_RNG GamRNG;
[32]31
[488]32/*!
33* \brief General conjugate exponential family posterior density.
[8]34
[488]35* More?...
36*/
[28]37
[488]38class eEF : public epdf
39{
40        public:
[32]41//      eEF() :epdf() {};
[488]42                //! default constructor
43                eEF () : epdf () {};
44                //! logarithm of the normalizing constant, \f$\mathcal{I}\f$
45                virtual double lognc() const = 0;
[565]46
[488]47                //!Evaluate normalized log-probability
[565]48                virtual double evallog_nn (const vec &val) const {
49                        bdm_error ("Not implemented");
50                        return 0.0;
51                }
52
[488]53                //!Evaluate normalized log-probability
54                virtual double evallog (const vec &val) const {
55                        double tmp;
56                        tmp = evallog_nn (val) - lognc();
57                        return tmp;
58                }
59                //!Evaluate normalized log-probability for many samples
60                virtual vec evallog_m (const mat &Val) const {
61                        vec x (Val.cols());
62                        for (int i = 0;i < Val.cols();i++) {x (i) = evallog_nn (Val.get_col (i)) ;}
63                        return x -lognc();
64                }
65                //!Evaluate normalized log-probability for many samples
66                virtual vec evallog_m (const Array<vec> &Val) const {
67                        vec x (Val.length());
68                        for (int i = 0;i < Val.length();i++) {x (i) = evallog_nn (Val (i)) ;}
69                        return x -lognc();
70                }
[565]71
[488]72                //!Power of the density, used e.g. to flatten the density
[565]73                virtual void pow (double p) {
74                        bdm_error ("Not implemented");
75                }
[488]76};
[8]77
[33]78
[170]79//! Estimator for Exponential family
[488]80class BMEF : public BM
81{
82        protected:
83                //! forgetting factor
84                double frg;
85                //! cached value of lognc() in the previous step (used in evaluation of \c ll )
86                double last_lognc;
87        public:
88                //! Default constructor (=empty constructor)
89                BMEF (double frg0 = 1.0) : BM (), frg (frg0) {}
90                //! Copy constructor
91                BMEF (const BMEF &B) : BM (B), frg (B.frg), last_lognc (B.last_lognc) {}
92                //!get statistics from another model
[565]93                virtual void set_statistics (const BMEF* BM0) {
94                        bdm_error ("Not implemented");
95                }
96
[488]97                //! Weighted update of sufficient statistics (Bayes rule)
98                virtual void bayes (const vec &data, const double w) {};
99                //original Bayes
100                void bayes (const vec &dt);
[565]101
[488]102                //!Flatten the posterior according to the given BMEF (of the same type!)
[565]103                virtual void flatten (const BMEF * B) {
104                        bdm_error ("Not implemented");
105                }
[198]106
[565]107                BMEF* _copy_ () const {
108                        bdm_error ("function _copy_ not implemented for this BM");
109                        return NULL;
110                }
[488]111};
[170]112
[504]113template<class sq_T, template <typename> class TEpdf>
[488]114class mlnorm;
[178]115
[488]116/*!
117* \brief Gaussian density with positive definite (decomposed) covariance matrix.
[8]118
[488]119* More?...
120*/
121template<class sq_T>
122class enorm : public eEF
123{
124        protected:
125                //! mean value
126                vec mu;
127                //! Covariance matrix in decomposed form
128                sq_T R;
129        public:
130                //!\name Constructors
131                //!@{
[270]132
[488]133                enorm () : eEF (), mu (), R () {};
134                enorm (const vec &mu, const sq_T &R) {set_parameters (mu, R);}
135                void set_parameters (const vec &mu, const sq_T &R);
[621]136                /*! Create Normal density
137                \f[ f(rv) = N(\mu, R) \f]
138                from structure
139                \code
140                class = 'enorm<ldmat>', (OR) 'enorm<chmat>', (OR) 'enorm<fsqmat>';
141                mu    = [];                  // mean value
142                R     = [];                  // variance, square matrix of appropriate dimension
143                \endcode
144                */
[488]145                void from_setting (const Setting &root);
146                void validate() {
[620]147                        bdm_assert (mu.length() == R.rows(), "mu and R parameters do not match");
[488]148                        dim = mu.length();
149                }
150                //!@}
[270]151
[488]152                //! \name Mathematical operations
153                //!@{
[270]154
[488]155                //! dupdate in exponential form (not really handy)
156                void dupdate (mat &v, double nu = 1.0);
[28]157
[488]158                vec sample() const;
[450]159
[488]160                double evallog_nn (const vec &val) const;
161                double lognc () const;
162                vec mean() const {return mu;}
163                vec variance() const {return diag (R.to_mat());}
[299]164//      mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
[504]165                shared_ptr<mpdf> condition ( const RV &rvn ) const;
166
167                // target not typed to mlnorm<sq_T, enorm<sq_T> > &
168                // because that doesn't compile (perhaps because we
169                // haven't finished defining enorm yet), but the type
170                // is required
171                void condition ( const RV &rvn, mpdf &target ) const;
172
173                shared_ptr<epdf> marginal (const RV &rvn ) const;
174                void marginal ( const RV &rvn, enorm<sq_T> &target ) const;
[488]175                //!@}
[270]176
[488]177                //! \name Access to attributes
178                //!@{
[270]179
[488]180                vec& _mu() {return mu;}
[583]181                const vec& _mu() const {return mu;}
[488]182                void set_mu (const vec mu0) { mu = mu0;}
183                sq_T& _R() {return R;}
184                const sq_T& _R() const {return R;}
185                //!@}
[28]186
[488]187};
[535]188UIREGISTER2 (enorm, chmat);
[529]189SHAREDPTR2 ( enorm, chmat );
[535]190UIREGISTER2 (enorm, ldmat);
[529]191SHAREDPTR2 ( enorm, ldmat );
[535]192UIREGISTER2 (enorm, fsqmat);
[529]193SHAREDPTR2 ( enorm, fsqmat );
[8]194
[388]195
[488]196/*!
197* \brief Gauss-inverse-Wishart density stored in LD form
[96]198
[488]199* For \f$p\f$-variate densities, given rv.count() should be \f$p\times\f$ V.rows().
200*
201*/
202class egiw : public eEF
203{
204        protected:
205                //! Extended information matrix of sufficient statistics
206                ldmat V;
207                //! Number of data records (degrees of freedom) of sufficient statistics
208                double nu;
209                //! Dimension of the output
210                int dimx;
211                //! Dimension of the regressor
212                int nPsi;
213        public:
214                //!\name Constructors
215                //!@{
216                egiw() : eEF() {};
217                egiw (int dimx0, ldmat V0, double nu0 = -1.0) : eEF() {set_parameters (dimx0, V0, nu0);};
[270]218
[629]219                void set_parameters (int dimx0, ldmat V0, double nu0 = -1.0);
[488]220                //!@}
[96]221
[488]222                vec sample() const;
223                vec mean() const;
224                vec variance() const;
[330]225
[488]226                //! LS estimate of \f$\theta\f$
227                vec est_theta() const;
[330]228
[488]229                //! Covariance of the LS estimate
230                ldmat est_theta_cov() const;
[330]231
[536]232                //! expected values of the linear coefficient and the covariance matrix are written to \c M and \c R , respectively
[488]233                void mean_mat (mat &M, mat&R) const;
234                //! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
235                double evallog_nn (const vec &val) const;
236                double lognc () const;
237                void pow (double p) {V *= p;nu *= p;};
[96]238
[488]239                //! \name Access attributes
240                //!@{
[270]241
[488]242                ldmat& _V() {return V;}
243                const ldmat& _V() const {return V;}
244                double& _nu()  {return nu;}
245                const double& _nu() const {return nu;}
[621]246                /*! Create Gauss-inverse-Wishart density
247                \f[ f(rv) = GiW(V,\nu) \f]
248                from structure
249                \code
250                class = 'egiw';
251                V     = [];               // square matrix
252                nu    = [];               // scalar \nu ((almost) degrees of freedom)
253                                                                  // when missing, it will be computed to obtain proper pdf
254                dimx  = [];               // dimension of the wishart part
255                rv = RV({'name'})         // description of RV
256                rvc = RV({'name'})        // description of RV in condition
257                \endcode
258                */
259               
[488]260                void from_setting (const Setting &set) {
[621]261                        epdf::from_setting(set);
262                        if (!UI::get (nu, set, "nu", UI::compulsory)) {nu=-1;}
[488]263                        UI::get (dimx, set, "dimx", UI::compulsory);
264                        mat V;
265                        UI::get (V, set, "V", UI::compulsory);
266                        set_parameters (dimx, V, nu);
267                }
268                //!@}
269};
[529]270UIREGISTER ( egiw );
271SHAREDPTR ( egiw );
[96]272
[488]273/*! \brief Dirichlet posterior density
[173]274
[488]275Continuous Dirichlet density of \f$n\f$-dimensional variable \f$x\f$
276\f[
277f(x|\beta) = \frac{\Gamma[\gamma]}{\prod_{i=1}^{n}\Gamma(\beta_i)} \prod_{i=1}^{n}x_i^{\beta_i-1}
278\f]
279where \f$\gamma=\sum_i \beta_i\f$.
280*/
281class eDirich: public eEF
282{
283        protected:
284                //!sufficient statistics
285                vec beta;
286        public:
287                //!\name Constructors
288                //!@{
[270]289
[488]290                eDirich () : eEF () {};
291                eDirich (const eDirich &D0) : eEF () {set_parameters (D0.beta);};
292                eDirich (const vec &beta0) {set_parameters (beta0);};
293                void set_parameters (const vec &beta0) {
294                        beta = beta0;
295                        dim = beta.length();
296                }
297                //!@}
[270]298
[634]299                //! using sampling procedure from wikipedia
[565]300                vec sample() const {
[634]301                        vec y(beta.length());
302                        for (int i=0; i<beta.length(); i++){
303                                GamRNG.setup(beta(i),1);
304                                y(i)=GamRNG.sample();
305                        }
306                        return y/sum(y);
[565]307                }
308
[488]309                vec mean() const {return beta / sum (beta);};
[634]310                vec variance() const {double gamma = sum (beta); return elem_mult (beta, (gamma-beta)) / (gamma*gamma* (gamma + 1));}
[488]311                //! In this instance, val is ...
312                double evallog_nn (const vec &val) const {
313                        double tmp; tmp = (beta - 1) * log (val);
314                        return tmp;
[565]315                }
316
[488]317                double lognc () const {
318                        double tmp;
319                        double gam = sum (beta);
320                        double lgb = 0.0;
321                        for (int i = 0;i < beta.length();i++) {lgb += lgamma (beta (i));}
322                        tmp = lgb - lgamma (gam);
323                        return tmp;
[565]324                }
325
[488]326                //!access function
327                vec& _beta()  {return beta;}
[634]328                /*! configuration structure
329                \code
330                class = 'eDirich';   
331                beta  = [];           //parametr beta
332                \endcode
333                */
334                void from_setting(const Setting &set){
335                        epdf::from_setting(set);
336                        UI::get(beta,set, "beta", UI::compulsory);
337                        validate();
338                }
339                void validate() {
340                        //check rv
341                        dim = beta.length();
342                }
[488]343};
[634]344UIREGISTER(eDirich);
[96]345
[181]346//! \brief Estimator for Multinomial density
[488]347class multiBM : public BMEF
348{
349        protected:
350                //! Conjugate prior and posterior
351                eDirich est;
352                //! Pointer inside est to sufficient statistics
353                vec &beta;
354        public:
355                //!Default constructor
356                multiBM () : BMEF (), est (), beta (est._beta()) {
357                        if (beta.length() > 0) {last_lognc = est.lognc();}
358                        else{last_lognc = 0.0;}
359                }
360                //!Copy constructor
361                multiBM (const multiBM &B) : BMEF (B), est (B.est), beta (est._beta()) {}
362                //! Sets sufficient statistics to match that of givefrom mB0
363                void set_statistics (const BM* mB0) {const multiBM* mB = dynamic_cast<const multiBM*> (mB0); beta = mB->beta;}
364                void bayes (const vec &dt) {
365                        if (frg < 1.0) {beta *= frg;last_lognc = est.lognc();}
366                        beta += dt;
367                        if (evalll) {ll = est.lognc() - last_lognc;}
368                }
369                double logpred (const vec &dt) const {
370                        eDirich pred (est);
371                        vec &beta = pred._beta();
[176]372
[488]373                        double lll;
374                        if (frg < 1.0)
375                                {beta *= frg;lll = pred.lognc();}
376                        else
377                                if (evalll) {lll = last_lognc;}
378                                else{lll = pred.lognc();}
[170]379
[488]380                        beta += dt;
381                        return pred.lognc() - lll;
382                }
383                void flatten (const BMEF* B) {
384                        const multiBM* E = dynamic_cast<const multiBM*> (B);
385                        // sum(beta) should be equal to sum(B.beta)
386                        const vec &Eb = E->beta;//const_cast<multiBM*> ( E )->_beta();
387                        beta *= (sum (Eb) / sum (beta));
388                        if (evalll) {last_lognc = est.lognc();}
389                }
[536]390                //! reimplemnetation of BM::posterior()
391                const eDirich& posterior() const {return est;};
392                //! constructor function               
[488]393                void set_parameters (const vec &beta0) {
394                        est.set_parameters (beta0);
395                        if (evalll) {last_lognc = est.lognc();}
396                }
397};
[170]398
[488]399/*!
400 \brief Gamma posterior density
[32]401
[488]402 Multivariate Gamma density as product of independent univariate densities.
403 \f[
404 f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i)
405 \f]
406*/
[32]407
[488]408class egamma : public eEF
409{
410        protected:
411                //! Vector \f$\alpha\f$
412                vec alpha;
413                //! Vector \f$\beta\f$
414                vec beta;
415        public :
416                //! \name Constructors
417                //!@{
418                egamma () : eEF (), alpha (0), beta (0) {};
419                egamma (const vec &a, const vec &b) {set_parameters (a, b);};
420                void set_parameters (const vec &a, const vec &b) {alpha = a, beta = b;dim = alpha.length();};
421                //!@}
[270]422
[488]423                vec sample() const;
424                double evallog (const vec &val) const;
425                double lognc () const;
[536]426                //! Returns pointer to internal alpha. Potentially dengerous: use with care!
[488]427                vec& _alpha() {return alpha;}
[536]428                //! Returns pointer to internal beta. Potentially dengerous: use with care!
[488]429                vec& _beta() {return beta;}
430                vec mean() const {return elem_div (alpha, beta);}
431                vec variance() const {return elem_div (alpha, elem_mult (beta, beta)); }
[225]432
[621]433                /*! Create Gamma density
434                \f[ f(rv|rvc) = \Gamma(\alpha, \beta) \f]
435                from structure
436                \code
437                class = 'egamma';
438                alpha = [...];         // vector of alpha
439                beta = [...];          // vector of beta
440                rv = RV({'name'})      // description of RV
441                \endcode
442                */
[488]443                void from_setting (const Setting &set) {
444                        epdf::from_setting (set); // reads rv
445                        UI::get (alpha, set, "alpha", UI::compulsory);
446                        UI::get (beta, set, "beta", UI::compulsory);
447                        validate();
448                }
449                void validate() {
[620]450                        bdm_assert (alpha.length() == beta.length(), "parameters do not match");
[488]451                        dim = alpha.length();
452                }
453};
454UIREGISTER (egamma);
[529]455SHAREDPTR ( egamma );
456
[488]457/*!
458 \brief Inverse-Gamma posterior density
[225]459
[488]460 Multivariate inverse-Gamma density as product of independent univariate densities.
461 \f[
462 f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i)
463 \f]
[283]464
[488]465Vector \f$\beta\f$ has different meaning (in fact it is 1/beta as used in definition of iG)
[225]466
[488]467 Inverse Gamma can be converted to Gamma using \f[
468 x\sim iG(a,b) => 1/x\sim G(a,1/b)
469\f]
470This relation is used in sampling.
471 */
[270]472
[488]473class eigamma : public egamma
474{
[294]475        protected:
[488]476        public :
477                //! \name Constructors
478                //! All constructors are inherited
479                //!@{
480                //!@}
[32]481
[488]482                vec sample() const {return 1.0 / egamma::sample();};
483                //! Returns poiter to alpha and beta. Potentially dangerous: use with care!
484                vec mean() const {return elem_div (beta, alpha - 1);}
485                vec variance() const {vec mea = mean(); return elem_div (elem_mult (mea, mea), alpha - 2);}
486};
487/*
488//! Weighted mixture of epdfs with external owned components.
489class emix : public epdf {
490protected:
491        int n;
492        vec &w;
493        Array<epdf*> Coms;
494public:
495//! Default constructor
496        emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
497        void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
498        vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
499};
500*/
501
[32]502//!  Uniform distributed density on a rectangular support
503
[488]504class euni: public epdf
505{
506        protected:
[32]507//! lower bound on support
[488]508                vec low;
[32]509//! upper bound on support
[488]510                vec high;
[32]511//! internal
[488]512                vec distance;
[32]513//! normalizing coefficients
[488]514                double nk;
[33]515//! cache of log( \c nk )
[488]516                double lnk;
517        public:
518                //! \name Constructors
519                //!@{
520                euni () : epdf () {}
521                euni (const vec &low0, const vec &high0) {set_parameters (low0, high0);}
522                void set_parameters (const vec &low0, const vec &high0) {
523                        distance = high0 - low0;
524                        low = low0;
525                        high = high0;
526                        nk = prod (1.0 / distance);
527                        lnk = log (nk);
528                        dim = low.length();
529                }
530                //!@}
[270]531
[488]532                double evallog (const vec &val) const  {
533                        if (any (val < low) && any (val > high)) {return inf;}
534                        else return lnk;
535                }
536                vec sample() const {
537                        vec smp (dim);
[270]538#pragma omp critical
[488]539                        UniRNG.sample_vector (dim , smp);
540                        return low + elem_mult (distance, smp);
541                }
542                //! set values of \c low and \c high
543                vec mean() const {return (high -low) / 2.0;}
544                vec variance() const {return (pow (high, 2) + pow (low, 2) + elem_mult (high, low)) / 3.0;}
[621]545                /*! Create Uniform density
546                \f[ f(rv) = U(low,high) \f]
547                from structure
548                 \code
549                 class = 'euni'
550                 high = [...];          // vector of upper bounds
551                 low = [...];           // vector of lower bounds
552                 rv = RV({'name'});     // description of RV
553                 \endcode
554                 */
[488]555                void from_setting (const Setting &set) {
556                        epdf::from_setting (set); // reads rv and rvc
[471]557
[488]558                        UI::get (high, set, "high", UI::compulsory);
559                        UI::get (low, set, "low", UI::compulsory);
[612]560                        set_parameters(low,high);
561                        validate();
[488]562                }
[612]563                void validate() {
564                        bdm_assert(high.length()==low.length(), "Incompatible high and low vectors");
565                        dim = high.length();
[620]566                        bdm_assert (min (distance) > 0.0, "bad support");
[612]567                }
[488]568};
[612]569UIREGISTER(euni);
[32]570
[488]571/*!
572 \brief Normal distributed linear function with linear function of mean value;
[32]573
[536]574 Mean value \f$ \mu=A*\mbox{rvc}+\mu_0 \f$.
[488]575*/
576template < class sq_T, template <typename> class TEpdf = enorm >
577class mlnorm : public mpdf_internal< TEpdf<sq_T> >
578{
579        protected:
580                //! Internal epdf that arise by conditioning on \c rvc
581                mat A;
[536]582                //! Constant additive term
[488]583                vec mu_const;
[487]584//                      vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
[488]585        public:
586                //! \name Constructors
587                //!@{
588                mlnorm() : mpdf_internal< TEpdf<sq_T> >() {};
589                mlnorm (const mat &A, const vec &mu0, const sq_T &R) : mpdf_internal< TEpdf<sq_T> >() {
590                        set_parameters (A, mu0, R);
591                }
[461]592
[488]593                //! Set \c A and \c R
[620]594                void set_parameters (const  mat &A0, const vec &mu0, const sq_T &R0) { 
[488]595                        this->iepdf.set_parameters (zeros (A0.rows()), R0);
596                        A = A0;
597                        mu_const = mu0;
598                        this->dimc = A0.cols();
599                }
600                //!@}
601                //! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
602                void condition (const vec &cond) {
603                        this->iepdf._mu() = A * cond + mu_const;
[487]604//R is already assigned;
[488]605                }
[198]606
[488]607                //!access function
[604]608                const vec& _mu_const() const {return mu_const;}
[488]609                //!access function
[604]610                const mat& _A() const {return A;}
[488]611                //!access function
[604]612                mat _R() const { return this->iepdf._R().to_mat(); }
[471]613
[536]614                //! Debug stream
[488]615                template<typename sq_M>
616                friend std::ostream &operator<< (std::ostream &os,  mlnorm<sq_M, enorm> &ml);
[8]617
[621]618                /*! Create Normal density with linear function of mean value
619                 \f[ f(rv|rvc) = N(A*rvc+const, R) \f]
620                from structure
621                \code
622                class = 'mlnorm<ldmat>', (OR) 'mlnorm<chmat>', (OR) 'mlnorm<fsqmat>';
623                A     = [];                  // matrix or vector of appropriate dimension
624                const = [];                  // vector of constant term
625                R     = [];                  // square matrix of appropriate dimension
626                \endcode
627                */
[488]628                void from_setting (const Setting &set) {
629                        mpdf::from_setting (set);
630
631                        UI::get (A, set, "A", UI::compulsory);
632                        UI::get (mu_const, set, "const", UI::compulsory);
633                        mat R0;
634                        UI::get (R0, set, "R", UI::compulsory);
635                        set_parameters (A, mu_const, R0);
[620]636                        validate();
[488]637                };
[620]638                void validate() {
639                        bdm_assert (A.rows() == mu_const.length(), "mlnorm: A vs. mu mismatch");
640                        bdm_assert (A.rows() == _R().rows(), "mlnorm: A vs. R mismatch");
641                       
642                }
[488]643};
[535]644UIREGISTER2 (mlnorm,ldmat);
[529]645SHAREDPTR2 ( mlnorm, ldmat );
[535]646UIREGISTER2 (mlnorm,fsqmat);
[529]647SHAREDPTR2 ( mlnorm, fsqmat );
[535]648UIREGISTER2 (mlnorm, chmat);
[529]649SHAREDPTR2 ( mlnorm, chmat );
[488]650
[280]651//! Mpdf with general function for mean value
[488]652template<class sq_T>
653class mgnorm : public mpdf_internal< enorm< sq_T > >
654{
[527]655        private:
[487]656//                      vec &mu; WHY NOT?
[527]657                shared_ptr<fnc> g;
658
[488]659        public:
660                //!default constructor
661                mgnorm() : mpdf_internal<enorm<sq_T> >() { }
662                //!set mean function
[527]663                inline void set_parameters (const shared_ptr<fnc> &g0, const sq_T &R0);
[488]664                inline void condition (const vec &cond);
[357]665
666
[621]667                /*! Create Normal density with given function of mean value
668                \f[ f(rv|rvc) = N( g(rvc), R) \f]
669                from structure
[488]670                \code
[621]671                class = 'mgnorm';
672                g.class =  'fnc';      // function for mean value evolution
673                g._fields_of_fnc = ...;
[357]674
[621]675                R = [1, 0;             // covariance matrix
676                        0, 1];
677                        --OR --
678                dR = [1, 1];           // diagonal of cavariance matrix
[357]679
[621]680                rv = RV({'name'})      // description of RV
681                rvc = RV({'name'})     // description of RV in condition
[488]682                \endcode
683                */
[357]684
[488]685                void from_setting (const Setting &set) {
[621]686                        mpdf::from_setting(set);
[527]687                        shared_ptr<fnc> g = UI::build<fnc> (set, "g", UI::compulsory);
[357]688
[488]689                        mat R;
690                        vec dR;
691                        if (UI::get (dR, set, "dR"))
692                                R = diag (dR);
693                        else
694                                UI::get (R, set, "R", UI::compulsory);
[280]695
[488]696                        set_parameters (g, R);
[621]697                        validate();
[488]698                }
[621]699                void validate() {
700                        bdm_assert(g->dimension()==this->dimension(),"incompatible function");
701                }
[488]702};
[357]703
[535]704UIREGISTER2 (mgnorm, chmat);
[529]705SHAREDPTR2 ( mgnorm, chmat );
[357]706
[262]707
[488]708/*! (Approximate) Student t density with linear function of mean value
[262]709
[488]710The internal epdf of this class is of the type of a Gaussian (enorm).
711However, each conditioning is trying to assure the best possible approximation by taking into account the zeta function. See [] for reference.
[270]712
[488]713Perhaps a moment-matching technique?
714*/
715class mlstudent : public mlnorm<ldmat, enorm>
716{
717        protected:
[536]718                //! Variable \f$ \Lambda \f$ from theory
[488]719                ldmat Lambda;
[536]720                //! Reference to variable \f$ R \f$
[488]721                ldmat &_R;
[536]722                //! Variable \f$ R_e \f$
[488]723                ldmat Re;
724        public:
725                mlstudent () : mlnorm<ldmat, enorm> (),
726                                Lambda (),      _R (iepdf._R()) {}
[536]727                //! constructor function
[488]728                void set_parameters (const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0) {
[576]729                        iepdf.set_parameters (mu0, R0);// was Lambda, why?
[488]730                        A = A0;
731                        mu_const = mu0;
732                        Re = R0;
733                        Lambda = Lambda0;
734                }
735                void condition (const vec &cond) {
736                        iepdf._mu() = A * cond + mu_const;
737                        double zeta;
738                        //ugly hack!
739                        if ( (cond.length() + 1) == Lambda.rows()) {
740                                zeta = Lambda.invqform (concat (cond, vec_1 (1.0)));
741                        } else {
742                                zeta = Lambda.invqform (cond);
[294]743                        }
[488]744                        _R = Re;
745                        _R *= (1 + zeta);// / ( nu ); << nu is in Re!!!!!!
746                };
[294]747
[613]748                void validate() {
[620]749                        bdm_assert (A.rows() == mu_const.length(), "mlstudent: A vs. mu mismatch");
750                        bdm_assert (_R.rows() == A.rows(), "mlstudent: A vs. R mismatch");
[613]751                       
752                }
[488]753};
754/*!
755\brief  Gamma random walk
[198]756
[488]757Mean value, \f$\mu\f$, of this density is given by \c rvc .
758Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
759This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
[32]760
[488]761The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
762*/
763class mgamma : public mpdf_internal<egamma>
764{
765        protected:
[461]766
[488]767                //! Constant \f$k\f$
768                double k;
[461]769
[488]770                //! cache of iepdf.beta
771                vec &_beta;
[32]772
[488]773        public:
774                //! Constructor
775                mgamma() : mpdf_internal<egamma>(), k (0),
776                                _beta (iepdf._beta()) {
777                }
[461]778
[488]779                //! Set value of \c k
780                void set_parameters (double k, const vec &beta0);
[461]781
[488]782                void condition (const vec &val) {_beta = k / val;};
[621]783                /*! Create Gamma density with conditional mean value
784                \f[ f(rv|rvc) = \Gamma(k, k/rvc) \f]
785                from structure
786                \code
787                  class = 'mgamma';
788                  beta = [...];          // vector of initial alpha
789                  k = 1.1;               // multiplicative constant k
790                  rv = RV({'name'})      // description of RV
791                  rvc = RV({'name'})     // description of RV in condition
792                 \endcode
793                */
[488]794                void from_setting (const Setting &set) {
795                        mpdf::from_setting (set); // reads rv and rvc
796                        vec betatmp; // ugly but necessary
797                        UI::get (betatmp, set, "beta", UI::compulsory);
798                        UI::get (k, set, "k", UI::compulsory);
799                        set_parameters (k, betatmp);
800                }
801};
802UIREGISTER (mgamma);
[529]803SHAREDPTR (mgamma);
[32]804
[488]805/*!
806\brief  Inverse-Gamma random walk
[225]807
[488]808Mean value, \f$ \mu \f$, of this density is given by \c rvc .
809Standard deviation of the random walk is proportional to one \f$ k \f$-th the mean.
810This is achieved by setting \f$ \alpha=\mu/k^2+2 \f$ and \f$ \beta=\mu(\alpha-1)\f$.
[461]811
[488]812The standard deviation of the walk is then: \f$ \mu/\sqrt(k)\f$.
813 */
814class migamma : public mpdf_internal<eigamma>
815{
816        protected:
817                //! Constant \f$k\f$
818                double k;
[461]819
[488]820                //! cache of iepdf.alpha
821                vec &_alpha;
[225]822
[488]823                //! cache of iepdf.beta
824                vec &_beta;
[461]825
[488]826        public:
827                //! \name Constructors
828                //!@{
829                migamma() : mpdf_internal<eigamma>(),
830                                k (0),
831                                _alpha (iepdf._alpha()),
832                                _beta (iepdf._beta()) {
833                }
[225]834
[488]835                migamma (const migamma &m) : mpdf_internal<eigamma>(),
836                                k (0),
837                                _alpha (iepdf._alpha()),
838                                _beta (iepdf._beta()) {
839                }
840                //!@}
[225]841
[488]842                //! Set value of \c k
843                void set_parameters (int len, double k0) {
844                        k = k0;
845                        iepdf.set_parameters ( (1.0 / (k*k) + 2.0) *ones (len) /*alpha*/, ones (len) /*beta*/);
846                        dimc = dimension();
847                };
848                void condition (const vec &val) {
849                        _beta = elem_mult (val, (_alpha - 1.0));
850                };
851};
[357]852
[60]853
[488]854/*!
855\brief  Gamma random walk around a fixed point
[60]856
[488]857Mean 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
858\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
[60]859
[488]860Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
861This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
[294]862
[488]863The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
864*/
865class mgamma_fix : public mgamma
866{
867        protected:
868                //! parameter l
869                double l;
870                //! reference vector
871                vec refl;
872        public:
873                //! Constructor
874                mgamma_fix () : mgamma (), refl () {};
875                //! Set value of \c k
876                void set_parameters (double k0 , vec ref0, double l0) {
877                        mgamma::set_parameters (k0, ref0);
878                        refl = pow (ref0, 1.0 - l0);l = l0;
879                        dimc = dimension();
880                };
[60]881
[488]882                void condition (const vec &val) {vec mean = elem_mult (refl, pow (val, l)); _beta = k / mean;};
883};
[60]884
[225]885
[488]886/*!
887\brief  Inverse-Gamma random walk around a fixed point
[225]888
[488]889Mean 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
890\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
[225]891
[488]892==== Check == vv =
893Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
894This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
[225]895
[488]896The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
897 */
898class migamma_ref : public migamma
899{
900        protected:
901                //! parameter l
902                double l;
903                //! reference vector
904                vec refl;
905        public:
906                //! Constructor
907                migamma_ref () : migamma (), refl () {};
908                //! Set value of \c k
909                void set_parameters (double k0 , vec ref0, double l0) {
910                        migamma::set_parameters (ref0.length(), k0);
911                        refl = pow (ref0, 1.0 - l0);
912                        l = l0;
913                        dimc = dimension();
914                };
[357]915
[488]916                void condition (const vec &val) {
917                        vec mean = elem_mult (refl, pow (val, l));
918                        migamma::condition (mean);
919                };
[357]920
921
[621]922                /*! Create inverse-Gamma density with conditional mean value
923                \f[ f(rv|rvc) = i\Gamma(k, k/(rvc^l \circ ref^{(1-l)}) \f]
924                from structure
[488]925                \code
[621]926                class = 'migamma_ref';
927                ref = [1e-5; 1e-5; 1e-2 1e-3];            // reference vector
928                l = 0.999;                                // constant l
929                k = 0.1;                                  // constant k
930                rv = RV({'name'})                         // description of RV
931                rvc = RV({'name'})                        // description of RV in condition
[488]932                \endcode
933                 */
934                void from_setting (const Setting &set);
[225]935
[488]936                // TODO dodelat void to_setting( Setting &set ) const;
937};
[357]938
939
[488]940UIREGISTER (migamma_ref);
[529]941SHAREDPTR (migamma_ref);
[294]942
[488]943/*! Log-Normal probability density
944 only allow diagonal covariances!
[294]945
[488]946Density of the form \f$ \log(x)\sim \mathcal{N}(\mu,\sigma^2) \f$ , i.e.
947\f[
948x \sim \frac{1}{x\sigma\sqrt{2\pi}}\exp{-\frac{1}{2\sigma^2}(\log(x)-\mu)}
949\f]
[294]950
[621]951Function from_setting loads mu and R in the same way as it does for enorm<>!
[488]952*/
953class elognorm: public enorm<ldmat>
954{
955        public:
956                vec sample() const {return exp (enorm<ldmat>::sample());};
957                vec mean() const {vec var = enorm<ldmat>::variance();return exp (mu - 0.5*var);};
[285]958
[488]959};
[285]960
[488]961/*!
962\brief  Log-Normal random walk
[285]963
[488]964Mean value, \f$\mu\f$, is...
[285]965
[488]966 */
967class mlognorm : public mpdf_internal<elognorm>
968{
969        protected:
970                //! parameter 1/2*sigma^2
971                double sig2;
[461]972
[488]973                //! access
974                vec &mu;
975        public:
976                //! Constructor
977                mlognorm() : mpdf_internal<elognorm>(),
978                                sig2 (0),
979                                mu (iepdf._mu()) {
980                }
[285]981
[488]982                //! Set value of \c k
983                void set_parameters (int size, double k) {
984                        sig2 = 0.5 * log (k * k + 1);
985                        iepdf.set_parameters (zeros (size), 2*sig2*eye (size));
[285]986
[488]987                        dimc = size;
988                };
[357]989
[488]990                void condition (const vec &val) {
991                        mu = log (val) - sig2;//elem_mult ( refl,pow ( val,l ) );
992                };
[357]993
[621]994                /*! Create logNormal random Walk
995                \f[ f(rv|rvc) = log\mathcal{N}( \log(rvc)-0.5\log(k^2+1), k I) \f]
996                from structure
[488]997                \code
[621]998                class = 'mlognorm';
999                k   = 0.1;               // "variance" k
1000                mu0 = 0.1;               // Initial value of mean
1001                rv  = RV({'name'})       // description of RV
1002                rvc = RV({'name'})       // description of RV in condition
[488]1003                \endcode
[621]1004                */
[488]1005                void from_setting (const Setting &set);
[357]1006
[488]1007                // TODO dodelat void to_setting( Setting &set ) const;
[285]1008
[488]1009};
[294]1010
[488]1011UIREGISTER (mlognorm);
[529]1012SHAREDPTR (mlognorm);
[294]1013
[488]1014/*! inverse Wishart density defined on Choleski decomposition
1015
1016*/
1017class eWishartCh : public epdf
1018{
1019        protected:
1020                //! Upper-Triagle of Choleski decomposition of \f$ \Psi \f$
1021                chmat Y;
1022                //! dimension of matrix  \f$ \Psi \f$
1023                int p;
1024                //! degrees of freedom  \f$ \nu \f$
1025                double delta;
1026        public:
[536]1027                //! Set internal structures
[488]1028                void set_parameters (const mat &Y0, const double delta0) {Y = chmat (Y0);delta = delta0; p = Y.rows(); dim = p * p; }
[536]1029                //! Sample matrix argument
[488]1030                mat sample_mat() const {
1031                        mat X = zeros (p, p);
1032
1033                        //sample diagonal
1034                        for (int i = 0;i < p;i++) {
1035                                GamRNG.setup (0.5* (delta - i) , 0.5);   // no +1 !! index if from 0
[294]1036#pragma omp critical
[488]1037                                X (i, i) = sqrt (GamRNG());
1038                        }
1039                        //do the rest
1040                        for (int i = 0;i < p;i++) {
1041                                for (int j = i + 1;j < p;j++) {
[294]1042#pragma omp critical
[488]1043                                        X (i, j) = NorRNG.sample();
[294]1044                                }
1045                        }
[488]1046                        return X*Y._Ch();// return upper triangular part of the decomposition
1047                }
1048                vec sample () const {
1049                        return vec (sample_mat()._data(), p*p);
1050                }
1051                //! fast access function y0 will be copied into Y.Ch.
1052                void setY (const mat &Ch0) {copy_vector (dim, Ch0._data(), Y._Ch()._data());}
1053                //! fast access function y0 will be copied into Y.Ch.
1054                void _setY (const vec &ch0) {copy_vector (dim, ch0._data(), Y._Ch()._data()); }
1055                //! access function
1056                const chmat& getY() const {return Y;}
1057};
[294]1058
[536]1059//! Inverse Wishart on Choleski decomposition
1060/*! Being computed by conversion from `standard' Wishart
1061*/
[488]1062class eiWishartCh: public epdf
1063{
1064        protected:
[536]1065                //! Internal instance of Wishart density
[488]1066                eWishartCh W;
[536]1067                //! size of Ch
[488]1068                int p;
[536]1069                //! parameter delta
[488]1070                double delta;
1071        public:
[536]1072                //! constructor function
[488]1073                void set_parameters (const mat &Y0, const double delta0) {
1074                        delta = delta0;
1075                        W.set_parameters (inv (Y0), delta0);
1076                        dim = W.dimension(); p = Y0.rows();
1077                }
1078                vec sample() const {mat iCh; iCh = inv (W.sample_mat()); return vec (iCh._data(), dim);}
[536]1079                //! access function
[488]1080                void _setY (const vec &y0) {
1081                        mat Ch (p, p);
1082                        mat iCh (p, p);
1083                        copy_vector (dim, y0._data(), Ch._data());
1084
1085                        iCh = inv (Ch);
1086                        W.setY (iCh);
1087                }
1088                virtual double evallog (const vec &val) const {
1089                        chmat X (p);
1090                        const chmat& Y = W.getY();
1091
1092                        copy_vector (p*p, val._data(), X._Ch()._data());
1093                        chmat iX (p);X.inv (iX);
1094                        // compute
[294]1095//                              \frac{ |\Psi|^{m/2}|X|^{-(m+p+1)/2}e^{-tr(\Psi X^{-1})/2} }{ 2^{mp/2}\Gamma_p(m/2)},
[488]1096                        mat M = Y.to_mat() * iX.to_mat();
[285]1097
[488]1098                        double log1 = 0.5 * p * (2 * Y.logdet()) - 0.5 * (delta + p + 1) * (2 * X.logdet()) - 0.5 * trace (M);
1099                        //Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!
[461]1100
[488]1101                        /*                              if (0) {
1102                                                                mat XX=X.to_mat();
1103                                                                mat YY=Y.to_mat();
[461]1104
[488]1105                                                                double log2 = 0.5*p*log(det(YY))-0.5*(delta+p+1)*log(det(XX))-0.5*trace(YY*inv(XX));
1106                                                                cout << log1 << "," << log2 << endl;
1107                                                        }*/
1108                        return log1;
1109                };
[285]1110
[488]1111};
1112
[536]1113//! Random Walk on inverse Wishart
[510]1114class rwiWishartCh : public mpdf_internal<eiWishartCh>
[488]1115{
1116        protected:
1117                //!square root of \f$ \nu-p-1 \f$ - needed for computation of \f$ \Psi \f$ from conditions
1118                double sqd;
[536]1119                //!reference point for diagonal
[488]1120                vec refl;
[536]1121                //! power of the reference
[488]1122                double l;
[536]1123                //! dimension
[488]1124                int p;
1125
1126        public:
[510]1127                rwiWishartCh() : sqd (0), l (0), p (0) {}
[536]1128                //! constructor function
[488]1129                void set_parameters (int p0, double k, vec ref0, double l0) {
1130                        p = p0;
1131                        double delta = 2 / (k * k) + p + 3;
1132                        sqd = sqrt (delta - p - 1);
1133                        l = l0;
1134                        refl = pow (ref0, 1 - l);
1135
[510]1136                        iepdf.set_parameters (eye (p), delta);
1137                        dimc = iepdf.dimension();
[488]1138                }
1139                void condition (const vec &c) {
1140                        vec z = c;
1141                        int ri = 0;
1142                        for (int i = 0;i < p*p;i += (p + 1)) {//trace diagonal element
1143                                z (i) = pow (z (i), l) * refl (ri);
1144                                ri++;
[294]1145                        }
[285]1146
[510]1147                        iepdf._setY (sqd*z);
[488]1148                }
1149};
1150
[32]1151//! Switch between various resampling methods.
[488]1152enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
1153/*!
1154\brief Weighted empirical density
[32]1155
[488]1156Used e.g. in particle filters.
1157*/
1158class eEmp: public epdf
1159{
1160        protected :
1161                //! Number of particles
1162                int n;
1163                //! Sample weights \f$w\f$
1164                vec w;
1165                //! Samples \f$x^{(i)}, i=1..n\f$
1166                Array<vec> samples;
1167        public:
1168                //! \name Constructors
1169                //!@{
1170                eEmp () : epdf (), w (), samples () {};
1171                //! copy constructor
1172                eEmp (const eEmp &e) : epdf (e), w (e.w), samples (e.samples) {};
1173                //!@}
[280]1174
[488]1175                //! Set samples and weights
1176                void set_statistics (const vec &w0, const epdf &pdf0);
1177                //! Set samples and weights
1178                void set_statistics (const epdf &pdf0 , int n) {set_statistics (ones (n) / n, pdf0);};
1179                //! Set sample
1180                void set_samples (const epdf* pdf0);
1181                //! Set sample
1182                void set_parameters (int n0, bool copy = true) {n = n0; w.set_size (n0, copy);samples.set_size (n0, copy);};
[569]1183                //! Set samples
1184                void set_parameters (const Array<vec> &Av) {
[620]1185                        bdm_assert(Av.size()>0,"Empty samples"); 
[569]1186                        n = Av.size(); 
1187                        epdf::set_parameters(Av(0).length());
1188                        w=1/n*ones(n);
1189                        samples=Av;
1190                };
[488]1191                //! Potentially dangerous, use with care.
1192                vec& _w()  {return w;};
1193                //! Potentially dangerous, use with care.
1194                const vec& _w() const {return w;};
1195                //! access function
1196                Array<vec>& _samples() {return samples;};
1197                //! access function
1198                const Array<vec>& _samples() const {return samples;};
1199                //! 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.
1200                ivec resample (RESAMPLING_METHOD method = SYSTEMATIC);
[565]1201
1202                //! inherited operation : NOT implemented
1203                vec sample() const {
1204                        bdm_error ("Not implemented");
1205                        return vec();
1206                }
1207
1208                //! inherited operation : NOT implemented
1209                double evallog (const vec &val) const {
1210                        bdm_error ("Not implemented");
1211                        return 0.0;
1212                }
1213
[488]1214                vec mean() const {
1215                        vec pom = zeros (dim);
1216                        for (int i = 0;i < n;i++) {pom += samples (i) * w (i);}
1217                        return pom;
1218                }
1219                vec variance() const {
1220                        vec pom = zeros (dim);
1221                        for (int i = 0;i < n;i++) {pom += pow (samples (i), 2) * w (i);}
1222                        return pom -pow (mean(), 2);
1223                }
1224                //! For this class, qbounds are minimum and maximum value of the population!
1225                void qbounds (vec &lb, vec &ub, double perc = 0.95) const {
1226                        // lb in inf so than it will be pushed below;
1227                        lb.set_size (dim);
1228                        ub.set_size (dim);
1229                        lb = std::numeric_limits<double>::infinity();
1230                        ub = -std::numeric_limits<double>::infinity();
1231                        int j;
1232                        for (int i = 0;i < n;i++) {
1233                                for (j = 0;j < dim; j++) {
1234                                        if (samples (i) (j) < lb (j)) {lb (j) = samples (i) (j);}
1235                                        if (samples (i) (j) > ub (j)) {ub (j) = samples (i) (j);}
[294]1236                                }
1237                        }
[488]1238                }
1239};
[32]1240
1241
[8]1242////////////////////////
1243
[488]1244template<class sq_T>
1245void enorm<sq_T>::set_parameters (const vec &mu0, const sq_T &R0)
1246{
[28]1247//Fixme test dimensions of mu0 and R0;
[488]1248        mu = mu0;
1249        R = R0;
1250        validate();
1251};
[8]1252
[488]1253template<class sq_T>
1254void enorm<sq_T>::from_setting (const Setting &set)
1255{
1256        epdf::from_setting (set); //reads rv
[384]1257
[488]1258        UI::get (mu, set, "mu", UI::compulsory);
1259        mat Rtmp;// necessary for conversion
1260        UI::get (Rtmp, set, "R", UI::compulsory);
1261        R = Rtmp; // conversion
1262        validate();
1263}
[8]1264
[488]1265template<class sq_T>
1266void enorm<sq_T>::dupdate (mat &v, double nu)
1267{
1268        //
1269};
1270
[178]1271// template<class sq_T>
1272// void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
1273//      //
1274// };
[8]1275
[488]1276template<class sq_T>
1277vec enorm<sq_T>::sample() const
1278{
1279        vec x (dim);
[270]1280#pragma omp critical
[488]1281        NorRNG.sample_vector (dim, x);
1282        vec smp = R.sqrt_mult (x);
[12]1283
[488]1284        smp += mu;
1285        return smp;
1286};
[8]1287
[214]1288// template<class sq_T>
1289// double enorm<sq_T>::eval ( const vec &val ) const {
1290//      double pdfl,e;
1291//      pdfl = evallog ( val );
1292//      e = exp ( pdfl );
1293//      return e;
1294// };
[8]1295
[488]1296template<class sq_T>
1297double enorm<sq_T>::evallog_nn (const vec &val) const
1298{
1299        // 1.83787706640935 = log(2pi)
1300        double tmp = -0.5 * (R.invqform (mu - val));// - lognc();
1301        return  tmp;
1302};
[28]1303
[488]1304template<class sq_T>
1305inline double enorm<sq_T>::lognc () const
1306{
1307        // 1.83787706640935 = log(2pi)
1308        double tmp = 0.5 * (R.cols() * 1.83787706640935 + R.logdet());
1309        return tmp;
1310};
[28]1311
[8]1312
[192]1313// template<class sq_T>
1314// vec mlnorm<sq_T>::samplecond (const  vec &cond, double &lik ) {
1315//      this->condition ( cond );
1316//      vec smp = epdf.sample();
1317//      lik = epdf.eval ( smp );
1318//      return smp;
1319// }
[8]1320
[192]1321// template<class sq_T>
1322// mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) {
1323//      int i;
1324//      int dim = rv.count();
1325//      mat Smp ( dim,n );
1326//      vec smp ( dim );
1327//      this->condition ( cond );
[198]1328//
[192]1329//      for ( i=0; i<n; i++ ) {
1330//              smp = epdf.sample();
1331//              lik ( i ) = epdf.eval ( smp );
1332//              Smp.set_col ( i ,smp );
1333//      }
[198]1334//
[192]1335//      return Smp;
1336// }
[28]1337
[8]1338
[488]1339template<class sq_T>
[504]1340shared_ptr<epdf> enorm<sq_T>::marginal ( const RV &rvn ) const
[488]1341{
[504]1342        enorm<sq_T> *tmp = new enorm<sq_T> ();
1343        shared_ptr<epdf> narrow(tmp);
1344        marginal ( rvn, *tmp );
1345        return narrow;
1346}
1347
1348template<class sq_T>
1349void enorm<sq_T>::marginal ( const RV &rvn, enorm<sq_T> &target ) const
1350{
[620]1351        bdm_assert (isnamed(), "rv description is not assigned");
[488]1352        ivec irvn = rvn.dataind (rv);
[178]1353
[504]1354        sq_T Rn (R, irvn);  // select rows and columns of R
[280]1355
[504]1356        target.set_rv ( rvn );
1357        target.set_parameters (mu (irvn), Rn);
[488]1358}
[178]1359
[488]1360template<class sq_T>
[504]1361shared_ptr<mpdf> enorm<sq_T>::condition ( const RV &rvn ) const
[488]1362{
[504]1363        mlnorm<sq_T> *tmp = new mlnorm<sq_T> ();
1364        shared_ptr<mpdf> narrow(tmp);
1365        condition ( rvn, *tmp );
1366        return narrow;
1367}
[178]1368
[504]1369template<class sq_T>
1370void enorm<sq_T>::condition ( const RV &rvn, mpdf &target ) const
1371{
1372        typedef mlnorm<sq_T> TMlnorm;
1373
[620]1374        bdm_assert (isnamed(), "rvs are not assigned");
[504]1375        TMlnorm &uptarget = dynamic_cast<TMlnorm &>(target);
[270]1376
[488]1377        RV rvc = rv.subt (rvn);
[620]1378        bdm_assert ( (rvc._dsize() + rvn._dsize() == rv._dsize()), "wrong rvn");
[488]1379        //Permutation vector of the new R
1380        ivec irvn = rvn.dataind (rv);
1381        ivec irvc = rvc.dataind (rv);
1382        ivec perm = concat (irvn , irvc);
1383        sq_T Rn (R, perm);
[178]1384
[488]1385        //fixme - could this be done in general for all sq_T?
1386        mat S = Rn.to_mat();
1387        //fixme
1388        int n = rvn._dsize() - 1;
1389        int end = R.rows() - 1;
1390        mat S11 = S.get (0, n, 0, n);
1391        mat S12 = S.get (0, n , rvn._dsize(), end);
1392        mat S22 = S.get (rvn._dsize(), end, rvn._dsize(), end);
[178]1393
[488]1394        vec mu1 = mu (irvn);
1395        vec mu2 = mu (irvc);
1396        mat A = S12 * inv (S22);
1397        sq_T R_n (S11 - A *S12.T());
[178]1398
[504]1399        uptarget.set_rv (rvn);
1400        uptarget.set_rvc (rvc);
1401        uptarget.set_parameters (A, mu1 - A*mu2, R_n);
[488]1402}
[178]1403
[488]1404////
1405///////
1406template<class sq_T>
[527]1407void mgnorm<sq_T >::set_parameters (const shared_ptr<fnc> &g0, const sq_T &R0) {
1408        g = g0;
1409        this->iepdf.set_parameters (zeros (g->dimension()), R0);
1410}
1411
[488]1412template<class sq_T>
1413void mgnorm<sq_T >::condition (const vec &cond) {this->iepdf._mu() = g->eval (cond);};
[28]1414
[536]1415//! \todo unify this stuff with to_string()
[488]1416template<class sq_T>
1417std::ostream &operator<< (std::ostream &os,  mlnorm<sq_T> &ml)
1418{
1419        os << "A:" << ml.A << endl;
1420        os << "mu:" << ml.mu_const << endl;
1421        os << "R:" << ml._R() << endl;
1422        return os;
1423};
[28]1424
[254]1425}
[8]1426#endif //EF_H
Note: See TracBrowser for help on using the browser.