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

Revision 735, 41.7 kB (checked in by smidl, 14 years ago)

Initial implementation of mixef.init() -- does not work due to missing to_setting

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