- Main Page
- Related Pages
- Classes
- Files

- File List
- File Members

# exp_family.h

Go to the documentation of this file.

    00001
    00013 #ifndef EF_H
    00014 #define EF_H
    00015
    00016
    00017 #include "../shared_ptr.h"
    00018 #include "../base/bdmbase.h"
    00019 #include "../math/chmat.h"
    00020
    00021 namespace bdm
    00022 {
    00023
    00024
    00026 extern Uniform_RNG UniRNG;
    00028 extern Normal_RNG NorRNG;
    00030 extern Gamma_RNG GamRNG;
    00031
    00038 class eEF : public epdf
    00039 {
    00040         public:
    00041 //      eEF() :epdf() {};
    00043                 eEF () : epdf () {};
    00045                 virtual double lognc() const = 0;
    00046
    00048                 virtual double evallog_nn (const vec &val) const {
    00049                         bdm_error ("Not implemented");
    00050                         return 0.0;
    00051                 }
    00052
    00054                 virtual double evallog (const vec &val) const {
    00055                         double tmp;
    00056                         tmp = evallog_nn (val) - lognc();
    00057                         return tmp;
    00058                 }
    00060                 virtual vec evallog_m (const mat &Val) const {
    00061                         vec x (Val.cols());
    00062                         for (int i = 0;i < Val.cols();i++) {x (i) = evallog_nn (Val.get_col (i)) ;}
    00063                         return x -lognc();
    00064                 }
    00066                 virtual vec evallog_m (const Array<vec> &Val) const {
    00067                         vec x (Val.length());
    00068                         for (int i = 0;i < Val.length();i++) {x (i) = evallog_nn (Val (i)) ;}
    00069                         return x -lognc();
    00070                 }
    00071
    00073                 virtual void pow (double p) {
    00074                         bdm_error ("Not implemented");
    00075                 }
    00076 };
    00077
    00078
    00080 class BMEF : public BM
    00081 {
    00082         protected:
    00084                 double frg;
    00086                 double last_lognc;
    00087         public:
    00089                 BMEF (double frg0 = 1.0) : BM (), frg (frg0) {}
    00091                 BMEF (const BMEF &B) : BM (B), frg (B.frg), last_lognc (B.last_lognc) {}
    00093                 virtual void set_statistics (const BMEF* BM0) {
    00094                         bdm_error ("Not implemented");
    00095                 }
    00096
    00098                 virtual void bayes (const vec &data, const double w) {};
    00099                 //original Bayes
    00100                 void bayes (const vec &dt);
    00101
    00103                 virtual void flatten (const BMEF * B) {
    00104                         bdm_error ("Not implemented");
    00105                 }
    00106
    00107                 BMEF* _copy_ () const {
    00108                         bdm_error ("function _copy_ not implemented for this BM");
    00109                         return NULL;
    00110                 }
    00111 };
    00112
    00113 template<class sq_T, template <typename> class TEpdf>
    00114 class mlnorm;
    00115
    00121 template<class sq_T>
    00122 class enorm : public eEF
    00123 {
    00124         protected:
    00126                 vec mu;
    00128                 sq_T R;
    00129         public:
    00132
    00133                 enorm () : eEF (), mu (), R () {};
    00134                 enorm (const vec &mu, const sq_T &R) {set_parameters (mu, R);}
    00135                 void set_parameters (const vec &mu, const sq_T &R);
    00145                 void from_setting (const Setting &root);
    00146                 void validate() {
    00147                         bdm_assert (mu.length() == R.rows(), "mu and R parameters do not match");
    00148                         dim = mu.length();
    00149                 }
    00151
    00154
    00156                 void dupdate (mat &v, double nu = 1.0);
    00157
    00158                 vec sample() const;
    00159
    00160                 double evallog_nn (const vec &val) const;
    00161                 double lognc () const;
    00162                 vec mean() const {return mu;}
    00163                 vec variance() const {return diag (R.to_mat());}
    00164 //      mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
    00165                 shared_ptr<mpdf> condition ( const RV &rvn ) const;
    00166
    00167                 // target not typed to mlnorm<sq_T, enorm<sq_T> > &
    00168                 // because that doesn't compile (perhaps because we
    00169                 // haven't finished defining enorm yet), but the type
    00170                 // is required
    00171                 void condition ( const RV &rvn, mpdf &target ) const;
    00172
    00173                 shared_ptr<epdf> marginal (const RV &rvn ) const;
    00174                 void marginal ( const RV &rvn, enorm<sq_T> &target ) const;
    00176
    00179
    00180                 vec& _mu() {return mu;}
    00181                 const vec& _mu() const {return mu;}
    00182                 void set_mu (const vec mu0) { mu = mu0;}
    00183                 sq_T& _R() {return R;}
    00184                 const sq_T& _R() const {return R;}
    00186
    00187 };
    00188 UIREGISTER2 (enorm, chmat);
    00189 SHAREDPTR2 ( enorm, chmat );
    00190 UIREGISTER2 (enorm, ldmat);
    00191 SHAREDPTR2 ( enorm, ldmat );
    00192 UIREGISTER2 (enorm, fsqmat);
    00193 SHAREDPTR2 ( enorm, fsqmat );
    00194
    00195
    00202 class egiw : public eEF
    00203 {
    00204         protected:
    00206                 ldmat V;
    00208                 double nu;
    00210                 int dimx;
    00212                 int nPsi;
    00213         public:
    00216                 egiw() : eEF() {};
    00217                 egiw (int dimx0, ldmat V0, double nu0 = -1.0) : eEF() {set_parameters (dimx0, V0, nu0);};
    00218
    00219                 void set_parameters (int dimx0, ldmat V0, double nu0 = -1.0);
    00221
    00222                 vec sample() const;
    00223                 vec mean() const;
    00224                 vec variance() const;
    00225
    00227                 vec est_theta() const;
    00228
    00230                 ldmat est_theta_cov() const;
    00231
    00233                 void mean_mat (mat &M, mat&R) const;
    00235                 double evallog_nn (const vec &val) const;
    00236                 double lognc () const;
    00237                 void pow (double p) {V *= p;nu *= p;};
    00238
    00241
    00242                 ldmat& _V() {return V;}
    00243                 const ldmat& _V() const {return V;}
    00244                 double& _nu()  {return nu;}
    00245                 const double& _nu() const {return nu;}
    00260                 void from_setting (const Setting &set) {
    00261                         epdf::from_setting(set);
    00262                         if (!UI::get (nu, set, "nu", UI::compulsory)) {nu=-1;}
    00263                         UI::get (dimx, set, "dimx", UI::compulsory);
    00264                         mat V;
    00265                         UI::get (V, set, "V", UI::compulsory);
    00266                         set_parameters (dimx, V, nu);
    00267                 }
    00269 };
    00270 UIREGISTER ( egiw );
    00271 SHAREDPTR ( egiw );
    00272
    00281 class eDirich: public eEF
    00282 {
    00283         protected:
    00285                 vec beta;
    00286         public:
    00289
    00290                 eDirich () : eEF () {};
    00291                 eDirich (const eDirich &D0) : eEF () {set_parameters (D0.beta);};
    00292                 eDirich (const vec &beta0) {set_parameters (beta0);};
    00293                 void set_parameters (const vec &beta0) {
    00294                         beta = beta0;
    00295                         dim = beta.length();
    00296                 }
    00298
    00300                 vec sample() const {
    00301                         vec y(beta.length());
    00302                         for (int i=0; i<beta.length(); i++){
    00303                                 GamRNG.setup(beta(i),1);
    00304                                 #pragma omp critical
    00305                                 y(i)=GamRNG();
    00306                         }
    00307                         return y/sum(y);
    00308                 }
    00309
    00310                 vec mean() const {return beta / sum (beta);};
    00311                 vec variance() const {double gamma = sum (beta); return elem_mult (beta, (gamma-beta)) / (gamma*gamma* (gamma + 1));}
    00313                 double evallog_nn (const vec &val) const {
    00314                         double tmp; tmp = (beta - 1) * log (val);
    00315                         return tmp;
    00316                 }
    00317
    00318                 double lognc () const {
    00319                         double tmp;
    00320                         double gam = sum (beta);
    00321                         double lgb = 0.0;
    00322                         for (int i = 0;i < beta.length();i++) {lgb += lgamma (beta (i));}
    00323                         tmp = lgb - lgamma (gam);
    00324                         return tmp;
    00325                 }
    00326
    00328                 vec& _beta()  {return beta;}
    00335                 void from_setting(const Setting &set){
    00336                         epdf::from_setting(set);
    00337                         UI::get(beta,set, "beta", UI::compulsory);
    00338                         validate();
    00339                 }
    00340                 void validate() {
    00341                         //check rv
    00342                         dim = beta.length();
    00343                 }
    00344 };
    00345 UIREGISTER(eDirich);
    00346
    00358 class mDirich: public mpdf_internal<eDirich> {
    00359         protected:
    00361                 double k;
    00363                 vec &_beta;
    00365                 vec betac;
    00366         public:
    00367                 mDirich(): mpdf_internal<eDirich>(), _beta(iepdf._beta()){};
    00368                 void condition (const vec &val) {_beta =  val/k+betac; };
    00381                 void from_setting (const Setting &set) {
    00382                         mpdf::from_setting (set); // reads rv and rvc
    00383                         if (_rv()._dsize()>0){
    00384                                 rvc = _rv().copy_t(-1);
    00385                         }
    00386                         vec beta0;
    00387                         if (!UI::get (beta0, set, "beta0", UI::optional)){
    00388                                 beta0 = ones(_rv()._dsize());
    00389                         }
    00390                         if (!UI::get (betac, set, "betac", UI::optional)){
    00391                                 betac = 0.1*ones(_rv()._dsize());
    00392                         }
    00393                         _beta = beta0;
    00394
    00395                         UI::get (k, set, "k", UI::compulsory);
    00396                         validate();
    00397                 }
    00398                 void validate() {
    00399                         iepdf.validate();
    00400                         bdm_assert(_beta.length()==betac.length(),"beta0 and betac are not compatible");
    00401                         if (_rv()._dsize()>0){
    00402                                 bdm_assert( (_rv()._dsize()==dimension()) , "Size of rv does not match with beta");
    00403                         }
    00404                         dimc = _beta.length();
    00405                 };
    00406 };
    00407 UIREGISTER(mDirich);
    00408
    00410 class multiBM : public BMEF
    00411 {
    00412         protected:
    00414                 eDirich est;
    00416                 vec &beta;
    00417         public:
    00419                 multiBM () : BMEF (), est (), beta (est._beta()) {
    00420                         if (beta.length() > 0) {last_lognc = est.lognc();}
    00421                         else{last_lognc = 0.0;}
    00422                 }
    00424                 multiBM (const multiBM &B) : BMEF (B), est (B.est), beta (est._beta()) {}
    00426                 void set_statistics (const BM* mB0) {const multiBM* mB = dynamic_cast<const multiBM*> (mB0); beta = mB->beta;}
    00427                 void bayes (const vec &dt) {
    00428                         if (frg < 1.0) {beta *= frg;last_lognc = est.lognc();}
    00429                         beta += dt;
    00430                         if (evalll) {ll = est.lognc() - last_lognc;}
    00431                 }
    00432                 double logpred (const vec &dt) const {
    00433                         eDirich pred (est);
    00434                         vec &beta = pred._beta();
    00435
    00436                         double lll;
    00437                         if (frg < 1.0)
    00438                                 {beta *= frg;lll = pred.lognc();}
    00439                         else
    00440                                 if (evalll) {lll = last_lognc;}
    00441                                 else{lll = pred.lognc();}
    00442
    00443                         beta += dt;
    00444                         return pred.lognc() - lll;
    00445                 }
    00446                 void flatten (const BMEF* B) {
    00447                         const multiBM* E = dynamic_cast<const multiBM*> (B);
    00448                         // sum(beta) should be equal to sum(B.beta)
    00449                         const vec &Eb = E->beta;//const_cast<multiBM*> ( E )->_beta();
    00450                         beta *= (sum (Eb) / sum (beta));
    00451                         if (evalll) {last_lognc = est.lognc();}
    00452                 }
    00454                 const eDirich& posterior() const {return est;};
    00456                 void set_parameters (const vec &beta0) {
    00457                         est.set_parameters (beta0);
    00458                         if (evalll) {last_lognc = est.lognc();}
    00459                 }
    00460 };
    00461
    00471 class egamma : public eEF
    00472 {
    00473         protected:
    00475                 vec alpha;
    00477                 vec beta;
    00478         public :
    00481                 egamma () : eEF (), alpha (0), beta (0) {};
    00482                 egamma (const vec &a, const vec &b) {set_parameters (a, b);};
    00483                 void set_parameters (const vec &a, const vec &b) {alpha = a, beta = b;dim = alpha.length();};
    00485
    00486                 vec sample() const;
    00487                 double evallog (const vec &val) const;
    00488                 double lognc () const;
    00490                 vec& _alpha() {return alpha;}
    00492                 vec& _beta() {return beta;}
    00493                 vec mean() const {return elem_div (alpha, beta);}
    00494                 vec variance() const {return elem_div (alpha, elem_mult (beta, beta)); }
    00495
    00506                 void from_setting (const Setting &set) {
    00507                         epdf::from_setting (set); // reads rv
    00508                         UI::get (alpha, set, "alpha", UI::compulsory);
    00509                         UI::get (beta, set, "beta", UI::compulsory);
    00510                         validate();
    00511                 }
    00512                 void validate() {
    00513                         bdm_assert (alpha.length() == beta.length(), "parameters do not match");
    00514                         dim = alpha.length();
    00515                 }
    00516 };
    00517 UIREGISTER (egamma);
    00518 SHAREDPTR ( egamma );
    00519
    00536 class eigamma : public egamma
    00537 {
    00538         protected:
    00539         public :
    00544
    00545                 vec sample() const {return 1.0 / egamma::sample();};
    00547                 vec mean() const {return elem_div (beta, alpha - 1);}
    00548                 vec variance() const {vec mea = mean(); return elem_div (elem_mult (mea, mea), alpha - 2);}
    00549 };
    00550 /*
    00552 class emix : public epdf {
    00553 protected:
    00554         int n;
    00555         vec &w;
    00556         Array<epdf*> Coms;
    00557 public:
    00559         emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
    00560         void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
    00561         vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
    00562 };
    00563 */
    00564
    00566
    00567 class euni: public epdf
    00568 {
    00569         protected:
    00571                 vec low;
    00573                 vec high;
    00575                 vec distance;
    00577                 double nk;
    00579                 double lnk;
    00580         public:
    00583                 euni () : epdf () {}
    00584                 euni (const vec &low0, const vec &high0) {set_parameters (low0, high0);}
    00585                 void set_parameters (const vec &low0, const vec &high0) {
    00586                         distance = high0 - low0;
    00587                         low = low0;
    00588                         high = high0;
    00589                         nk = prod (1.0 / distance);
    00590                         lnk = log (nk);
    00591                         dim = low.length();
    00592                 }
    00594
    00595                 double evallog (const vec &val) const  {
    00596                         if (any (val < low) && any (val > high)) {return inf;}
    00597                         else return lnk;
    00598                 }
    00599                 vec sample() const {
    00600                         vec smp (dim);
    00601 #pragma omp critical
    00602                         UniRNG.sample_vector (dim , smp);
    00603                         return low + elem_mult (distance, smp);
    00604                 }
    00606                 vec mean() const {return (high -low) / 2.0;}
    00607                 vec variance() const {return (pow (high, 2) + pow (low, 2) + elem_mult (high, low)) / 3.0;}
    00618                 void from_setting (const Setting &set) {
    00619                         epdf::from_setting (set); // reads rv and rvc
    00620
    00621                         UI::get (high, set, "high", UI::compulsory);
    00622                         UI::get (low, set, "low", UI::compulsory);
    00623                         set_parameters(low,high);
    00624                         validate();
    00625                 }
    00626                 void validate() {
    00627                         bdm_assert(high.length()==low.length(), "Incompatible high and low vectors");
    00628                         dim = high.length();
    00629                         bdm_assert (min (distance) > 0.0, "bad support");
    00630                 }
    00631 };
    00632 UIREGISTER(euni);
    00633
    00639 template < class sq_T, template <typename> class TEpdf = enorm >
    00640 class mlnorm : public mpdf_internal< TEpdf<sq_T> >
    00641 {
    00642         protected:
    00644                 mat A;
    00646                 vec mu_const;
    00647 //                      vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
    00648         public:
    00651                 mlnorm() : mpdf_internal< TEpdf<sq_T> >() {};
    00652                 mlnorm (const mat &A, const vec &mu0, const sq_T &R) : mpdf_internal< TEpdf<sq_T> >() {
    00653                         set_parameters (A, mu0, R);
    00654                 }
    00655
    00657                 void set_parameters (const  mat &A0, const vec &mu0, const sq_T &R0) {
    00658                         this->iepdf.set_parameters (zeros (A0.rows()), R0);
    00659                         A = A0;
    00660                         mu_const = mu0;
    00661                         this->dimc = A0.cols();
    00662                 }
    00665                 void condition (const vec &cond) {
    00666                         this->iepdf._mu() = A * cond + mu_const;
    00667 //R is already assigned;
    00668                 }
    00669
    00671                 const vec& _mu_const() const {return mu_const;}
    00673                 const mat& _A() const {return A;}
    00675                 mat _R() const { return this->iepdf._R().to_mat(); }
    00676
    00678                 template<typename sq_M>
    00679                 friend std::ostream &operator<< (std::ostream &os,  mlnorm<sq_M, enorm> &ml);
    00680
    00691                 void from_setting (const Setting &set) {
    00692                         mpdf::from_setting (set);
    00693
    00694                         UI::get (A, set, "A", UI::compulsory);
    00695                         UI::get (mu_const, set, "const", UI::compulsory);
    00696                         mat R0;
    00697                         UI::get (R0, set, "R", UI::compulsory);
    00698                         set_parameters (A, mu_const, R0);
    00699                         validate();
    00700                 };
    00701                 void validate() {
    00702                         bdm_assert (A.rows() == mu_const.length(), "mlnorm: A vs. mu mismatch");
    00703                         bdm_assert (A.rows() == _R().rows(), "mlnorm: A vs. R mismatch");
    00704
    00705                 }
    00706 };
    00707 UIREGISTER2 (mlnorm,ldmat);
    00708 SHAREDPTR2 ( mlnorm, ldmat );
    00709 UIREGISTER2 (mlnorm,fsqmat);
    00710 SHAREDPTR2 ( mlnorm, fsqmat );
    00711 UIREGISTER2 (mlnorm, chmat);
    00712 SHAREDPTR2 ( mlnorm, chmat );
    00713
    00715 template<class sq_T>
    00716 class mgnorm : public mpdf_internal< enorm< sq_T > >
    00717 {
    00718         private:
    00719 //                      vec &mu; WHY NOT?
    00720                 shared_ptr<fnc> g;
    00721
    00722         public:
    00724                 mgnorm() : mpdf_internal<enorm<sq_T> >() { }
    00726                 inline void set_parameters (const shared_ptr<fnc> &g0, const sq_T &R0);
    00727                 inline void condition (const vec &cond);
    00728
    00729
    00748                 void from_setting (const Setting &set) {
    00749                         mpdf::from_setting(set);
    00750                         shared_ptr<fnc> g = UI::build<fnc> (set, "g", UI::compulsory);
    00751
    00752                         mat R;
    00753                         vec dR;
    00754                         if (UI::get (dR, set, "dR"))
    00755                                 R = diag (dR);
    00756                         else
    00757                                 UI::get (R, set, "R", UI::compulsory);
    00758
    00759                         set_parameters (g, R);
    00760                         validate();
    00761                 }
    00762                 void validate() {
    00763                         bdm_assert(g->dimension()==this->dimension(),"incompatible function");
    00764                 }
    00765 };
    00766
    00767 UIREGISTER2 (mgnorm, chmat);
    00768 SHAREDPTR2 ( mgnorm, chmat );
    00769
    00770
    00778 class mlstudent : public mlnorm<ldmat, enorm>
    00779 {
    00780         protected:
    00782                 ldmat Lambda;
    00784                 ldmat &_R;
    00786                 ldmat Re;
    00787         public:
    00788                 mlstudent () : mlnorm<ldmat, enorm> (),
    00789                                 Lambda (),      _R (iepdf._R()) {}
    00791                 void set_parameters (const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0) {
    00792                         iepdf.set_parameters (mu0, R0);// was Lambda, why?
    00793                         A = A0;
    00794                         mu_const = mu0;
    00795                         Re = R0;
    00796                         Lambda = Lambda0;
    00797                 }
    00798                 void condition (const vec &cond) {
    00799                         iepdf._mu() = A * cond + mu_const;
    00800                         double zeta;
    00801                         //ugly hack!
    00802                         if ( (cond.length() + 1) == Lambda.rows()) {
    00803                                 zeta = Lambda.invqform (concat (cond, vec_1 (1.0)));
    00804                         } else {
    00805                                 zeta = Lambda.invqform (cond);
    00806                         }
    00807                         _R = Re;
    00808                         _R *= (1 + zeta);// / ( nu ); << nu is in Re!!!!!!
    00809                 };
    00810
    00811                 void validate() {
    00812                         bdm_assert (A.rows() == mu_const.length(), "mlstudent: A vs. mu mismatch");
    00813                         bdm_assert (_R.rows() == A.rows(), "mlstudent: A vs. R mismatch");
    00814
    00815                 }
    00816 };
    00826 class mgamma : public mpdf_internal<egamma>
    00827 {
    00828         protected:
    00829
    00831                 double k;
    00832
    00834                 vec &_beta;
    00835
    00836         public:
    00838                 mgamma() : mpdf_internal<egamma>(), k (0),
    00839                                 _beta (iepdf._beta()) {
    00840                 }
    00841
    00843                 void set_parameters (double k, const vec &beta0);
    00844
    00845                 void condition (const vec &val) {_beta = k / val;};
    00857                 void from_setting (const Setting &set) {
    00858                         mpdf::from_setting (set); // reads rv and rvc
    00859                         vec betatmp; // ugly but necessary
    00860                         UI::get (betatmp, set, "beta", UI::compulsory);
    00861                         UI::get (k, set, "k", UI::compulsory);
    00862                         set_parameters (k, betatmp);
    00863                 }
    00864 };
    00865 UIREGISTER (mgamma);
    00866 SHAREDPTR (mgamma);
    00867
    00877 class migamma : public mpdf_internal<eigamma>
    00878 {
    00879         protected:
    00881                 double k;
    00882
    00884                 vec &_alpha;
    00885
    00887                 vec &_beta;
    00888
    00889         public:
    00892                 migamma() : mpdf_internal<eigamma>(),
    00893                                 k (0),
    00894                                 _alpha (iepdf._alpha()),
    00895                                 _beta (iepdf._beta()) {
    00896                 }
    00897
    00898                 migamma (const migamma &m) : mpdf_internal<eigamma>(),
    00899                                 k (0),
    00900                                 _alpha (iepdf._alpha()),
    00901                                 _beta (iepdf._beta()) {
    00902                 }
    00904
    00906                 void set_parameters (int len, double k0) {
    00907                         k = k0;
    00908                         iepdf.set_parameters ( (1.0 / (k*k) + 2.0) *ones (len) /*alpha*/, ones (len) /*beta*/);
    00909                         dimc = dimension();
    00910                 };
    00911                 void condition (const vec &val) {
    00912                         _beta = elem_mult (val, (_alpha - 1.0));
    00913                 };
    00914 };
    00915
    00916
    00928 class mgamma_fix : public mgamma
    00929 {
    00930         protected:
    00932                 double l;
    00934                 vec refl;
    00935         public:
    00937                 mgamma_fix () : mgamma (), refl () {};
    00939                 void set_parameters (double k0 , vec ref0, double l0) {
    00940                         mgamma::set_parameters (k0, ref0);
    00941                         refl = pow (ref0, 1.0 - l0);l = l0;
    00942                         dimc = dimension();
    00943                 };
    00944
    00945                 void condition (const vec &val) {vec mean = elem_mult (refl, pow (val, l)); _beta = k / mean;};
    00946 };
    00947
    00948
    00961 class migamma_ref : public migamma
    00962 {
    00963         protected:
    00965                 double l;
    00967                 vec refl;
    00968         public:
    00970                 migamma_ref () : migamma (), refl () {};
    00972                 void set_parameters (double k0 , vec ref0, double l0) {
    00973                         migamma::set_parameters (ref0.length(), k0);
    00974                         refl = pow (ref0, 1.0 - l0);
    00975                         l = l0;
    00976                         dimc = dimension();
    00977                 };
    00978
    00979                 void condition (const vec &val) {
    00980                         vec mean = elem_mult (refl, pow (val, l));
    00981                         migamma::condition (mean);
    00982                 };
    00983
    00984
    00997                 void from_setting (const Setting &set);
    00998
    00999                 // TODO dodelat void to_setting( Setting &set ) const;
    01000 };
    01001
    01002
    01003 UIREGISTER (migamma_ref);
    01004 SHAREDPTR (migamma_ref);
    01005
    01016 class elognorm: public enorm<ldmat>
    01017 {
    01018         public:
    01019                 vec sample() const {return exp (enorm<ldmat>::sample());};
    01020                 vec mean() const {vec var = enorm<ldmat>::variance();return exp (mu - 0.5*var);};
    01021
    01022 };
    01023
    01030 class mlognorm : public mpdf_internal<elognorm>
    01031 {
    01032         protected:
    01034                 double sig2;
    01035
    01037                 vec &mu;
    01038         public:
    01040                 mlognorm() : mpdf_internal<elognorm>(),
    01041                                 sig2 (0),
    01042                                 mu (iepdf._mu()) {
    01043                 }
    01044
    01046                 void set_parameters (int size, double k) {
    01047                         sig2 = 0.5 * log (k * k + 1);
    01048                         iepdf.set_parameters (zeros (size), 2*sig2*eye (size));
    01049
    01050                         dimc = size;
    01051                 };
    01052
    01053                 void condition (const vec &val) {
    01054                         mu = log (val) - sig2;//elem_mult ( refl,pow ( val,l ) );
    01055                 };
    01056
    01068                 void from_setting (const Setting &set);
    01069
    01070                 // TODO dodelat void to_setting( Setting &set ) const;
    01071
    01072 };
    01073
    01074 UIREGISTER (mlognorm);
    01075 SHAREDPTR (mlognorm);
    01076
    01080 class eWishartCh : public epdf
    01081 {
    01082         protected:
    01084                 chmat Y;
    01086                 int p;
    01088                 double delta;
    01089         public:
    01091                 void set_parameters (const mat &Y0, const double delta0) {Y = chmat (Y0);delta = delta0; p = Y.rows(); dim = p * p; }
    01093                 mat sample_mat() const {
    01094                         mat X = zeros (p, p);
    01095
    01096                         //sample diagonal
    01097                         for (int i = 0;i < p;i++) {
    01098                                 GamRNG.setup (0.5* (delta - i) , 0.5);   // no +1 !! index if from 0
    01099 #pragma omp critical
    01100                                 X (i, i) = sqrt (GamRNG());
    01101                         }
    01102                         //do the rest
    01103                         for (int i = 0;i < p;i++) {
    01104                                 for (int j = i + 1;j < p;j++) {
    01105 #pragma omp critical
    01106                                         X (i, j) = NorRNG.sample();
    01107                                 }
    01108                         }
    01109                         return X*Y._Ch();// return upper triangular part of the decomposition
    01110                 }
    01111                 vec sample () const {
    01112                         return vec (sample_mat()._data(), p*p);
    01113                 }
    01115                 void setY (const mat &Ch0) {copy_vector (dim, Ch0._data(), Y._Ch()._data());}
    01117                 void _setY (const vec &ch0) {copy_vector (dim, ch0._data(), Y._Ch()._data()); }
    01119                 const chmat& getY() const {return Y;}
    01120 };
    01121
    01123
    01125 class eiWishartCh: public epdf
    01126 {
    01127         protected:
    01129                 eWishartCh W;
    01131                 int p;
    01133                 double delta;
    01134         public:
    01136                 void set_parameters (const mat &Y0, const double delta0) {
    01137                         delta = delta0;
    01138                         W.set_parameters (inv (Y0), delta0);
    01139                         dim = W.dimension(); p = Y0.rows();
    01140                 }
    01141                 vec sample() const {mat iCh; iCh = inv (W.sample_mat()); return vec (iCh._data(), dim);}
    01143                 void _setY (const vec &y0) {
    01144                         mat Ch (p, p);
    01145                         mat iCh (p, p);
    01146                         copy_vector (dim, y0._data(), Ch._data());
    01147
    01148                         iCh = inv (Ch);
    01149                         W.setY (iCh);
    01150                 }
    01151                 virtual double evallog (const vec &val) const {
    01152                         chmat X (p);
    01153                         const chmat& Y = W.getY();
    01154
    01155                         copy_vector (p*p, val._data(), X._Ch()._data());
    01156                         chmat iX (p);X.inv (iX);
    01157                         // compute
    01158 //                              \frac{ |\Psi|^{m/2}|X|^{-(m+p+1)/2}e^{-tr(\Psi X^{-1})/2} }{ 2^{mp/2}\Gamma_p(m/2)},
    01159                         mat M = Y.to_mat() * iX.to_mat();
    01160
    01161                         double log1 = 0.5 * p * (2 * Y.logdet()) - 0.5 * (delta + p + 1) * (2 * X.logdet()) - 0.5 * trace (M);
    01162                         //Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!
    01163
    01164                         /*                              if (0) {
    01165                                                                 mat XX=X.to_mat();
    01166                                                                 mat YY=Y.to_mat();
    01167
    01168                                                                 double log2 = 0.5*p*log(det(YY))-0.5*(delta+p+1)*log(det(XX))-0.5*trace(YY*inv(XX));
    01169                                                                 cout << log1 << "," << log2 << endl;
    01170                                                         }*/
    01171                         return log1;
    01172                 };
    01173
    01174 };
    01175
    01177 class rwiWishartCh : public mpdf_internal<eiWishartCh>
    01178 {
    01179         protected:
    01181                 double sqd;
    01183                 vec refl;
    01185                 double l;
    01187                 int p;
    01188
    01189         public:
    01190                 rwiWishartCh() : sqd (0), l (0), p (0) {}
    01192                 void set_parameters (int p0, double k, vec ref0, double l0) {
    01193                         p = p0;
    01194                         double delta = 2 / (k * k) + p + 3;
    01195                         sqd = sqrt (delta - p - 1);
    01196                         l = l0;
    01197                         refl = pow (ref0, 1 - l);
    01198
    01199                         iepdf.set_parameters (eye (p), delta);
    01200                         dimc = iepdf.dimension();
    01201                 }
    01202                 void condition (const vec &c) {
    01203                         vec z = c;
    01204                         int ri = 0;
    01205                         for (int i = 0;i < p*p;i += (p + 1)) {//trace diagonal element
    01206                                 z (i) = pow (z (i), l) * refl (ri);
    01207                                 ri++;
    01208                         }
    01209
    01210                         iepdf._setY (sqd*z);
    01211                 }
    01212 };
    01213
    01215 enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
    01221 class eEmp: public epdf
    01222 {
    01223         protected :
    01225                 int n;
    01227                 vec w;
    01229                 Array<vec> samples;
    01230         public:
    01233                 eEmp () : epdf (), w (), samples () {};
    01235                 eEmp (const eEmp &e) : epdf (e), w (e.w), samples (e.samples) {};
    01237
    01239                 void set_statistics (const vec &w0, const epdf &pdf0);
    01241                 void set_statistics (const epdf &pdf0 , int n) {set_statistics (ones (n) / n, pdf0);};
    01243                 void set_samples (const epdf* pdf0);
    01245                 void set_parameters (int n0, bool copy = true) {n = n0; w.set_size (n0, copy);samples.set_size (n0, copy);};
    01247                 void set_parameters (const Array<vec> &Av) {
    01248                         bdm_assert(Av.size()>0,"Empty samples");
    01249                         n = Av.size();
    01250                         epdf::set_parameters(Av(0).length());
    01251                         w=1/n*ones(n);
    01252                         samples=Av;
    01253                 };
    01255                 vec& _w()  {return w;};
    01257                 const vec& _w() const {return w;};
    01259                 Array<vec>& _samples() {return samples;};
    01261                 const vec& _sample(int i) const {return samples(i);};
    01263                 const Array<vec>& _samples() const {return samples;};
    01266                 void resample ( ivec &index, RESAMPLING_METHOD method = SYSTEMATIC);
    01267
    01269                 void resample (RESAMPLING_METHOD method = SYSTEMATIC){ivec ind; resample(ind,method);};
    01270
    01272                 vec sample() const {
    01273                         bdm_error ("Not implemented");
    01274                         return vec();
    01275                 }
    01276
    01278                 double evallog (const vec &val) const {
    01279                         bdm_error ("Not implemented");
    01280                         return 0.0;
    01281                 }
    01282
    01283                 vec mean() const {
    01284                         vec pom = zeros (dim);
    01285                         for (int i = 0;i < n;i++) {pom += samples (i) * w (i);}
    01286                         return pom;
    01287                 }
    01288                 vec variance() const {
    01289                         vec pom = zeros (dim);
    01290                         for (int i = 0;i < n;i++) {pom += pow (samples (i), 2) * w (i);}
    01291                         return pom -pow (mean(), 2);
    01292                 }
    01294                 void qbounds (vec &lb, vec &ub, double perc = 0.95) const {
    01295                         // lb in inf so than it will be pushed below;
    01296                         lb.set_size (dim);
    01297                         ub.set_size (dim);
    01298                         lb = std::numeric_limits<double>::infinity();
    01299                         ub = -std::numeric_limits<double>::infinity();
    01300                         int j;
    01301                         for (int i = 0;i < n;i++) {
    01302                                 for (j = 0;j < dim; j++) {
    01303                                         if (samples (i) (j) < lb (j)) {lb (j) = samples (i) (j);}
    01304                                         if (samples (i) (j) > ub (j)) {ub (j) = samples (i) (j);}
    01305                                 }
    01306                         }
    01307                 }
    01308 };
    01309
    01310
    01312
    01313 template<class sq_T>
    01314 void enorm<sq_T>::set_parameters (const vec &mu0, const sq_T &R0)
    01315 {
    01316 //Fixme test dimensions of mu0 and R0;
    01317         mu = mu0;
    01318         R = R0;
    01319         validate();
    01320 };
    01321
    01322 template<class sq_T>
    01323 void enorm<sq_T>::from_setting (const Setting &set)
    01324 {
    01325         epdf::from_setting (set); //reads rv
    01326
    01327         UI::get (mu, set, "mu", UI::compulsory);
    01328         mat Rtmp;// necessary for conversion
    01329         UI::get (Rtmp, set, "R", UI::compulsory);
    01330         R = Rtmp; // conversion
    01331         validate();
    01332 }
    01333
    01334 template<class sq_T>
    01335 void enorm<sq_T>::dupdate (mat &v, double nu)
    01336 {
    01337         //
    01338 };
    01339
    01340 // template<class sq_T>
    01341 // void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
    01342 //      //
    01343 // };
    01344
    01345 template<class sq_T>
    01346 vec enorm<sq_T>::sample() const
    01347 {
    01348         vec x (dim);
    01349 #pragma omp critical
    01350         NorRNG.sample_vector (dim, x);
    01351         vec smp = R.sqrt_mult (x);
    01352
    01353         smp += mu;
    01354         return smp;
    01355 };
    01356
    01357 // template<class sq_T>
    01358 // double enorm<sq_T>::eval ( const vec &val ) const {
    01359 //      double pdfl,e;
    01360 //      pdfl = evallog ( val );
    01361 //      e = exp ( pdfl );
    01362 //      return e;
    01363 // };
    01364
    01365 template<class sq_T>
    01366 double enorm<sq_T>::evallog_nn (const vec &val) const
    01367 {
    01368         // 1.83787706640935 = log(2pi)
    01369         double tmp = -0.5 * (R.invqform (mu - val));// - lognc();
    01370         return  tmp;
    01371 };
    01372
    01373 template<class sq_T>
    01374 inline double enorm<sq_T>::lognc () const
    01375 {
    01376         // 1.83787706640935 = log(2pi)
    01377         double tmp = 0.5 * (R.cols() * 1.83787706640935 + R.logdet());
    01378         return tmp;
    01379 };
    01380
    01381
    01382 // template<class sq_T>
    01383 // vec mlnorm<sq_T>::samplecond (const  vec &cond, double &lik ) {
    01384 //      this->condition ( cond );
    01385 //      vec smp = epdf.sample();
    01386 //      lik = epdf.eval ( smp );
    01387 //      return smp;
    01388 // }
    01389
    01390 // template<class sq_T>
    01391 // mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) {
    01392 //      int i;
    01393 //      int dim = rv.count();
    01394 //      mat Smp ( dim,n );
    01395 //      vec smp ( dim );
    01396 //      this->condition ( cond );
    01397 //
    01398 //      for ( i=0; i<n; i++ ) {
    01399 //              smp = epdf.sample();
    01400 //              lik ( i ) = epdf.eval ( smp );
    01401 //              Smp.set_col ( i ,smp );
    01402 //      }
    01403 //
    01404 //      return Smp;
    01405 // }
    01406
    01407
    01408 template<class sq_T>
    01409 shared_ptr<epdf> enorm<sq_T>::marginal ( const RV &rvn ) const
    01410 {
    01411         enorm<sq_T> *tmp = new enorm<sq_T> ();
    01412         shared_ptr<epdf> narrow(tmp);
    01413         marginal ( rvn, *tmp );
    01414         return narrow;
    01415 }
    01416
    01417 template<class sq_T>
    01418 void enorm<sq_T>::marginal ( const RV &rvn, enorm<sq_T> &target ) const
    01419 {
    01420         bdm_assert (isnamed(), "rv description is not assigned");
    01421         ivec irvn = rvn.dataind (rv);
    01422
    01423         sq_T Rn (R, irvn);  // select rows and columns of R
    01424
    01425         target.set_rv ( rvn );
    01426         target.set_parameters (mu (irvn), Rn);
    01427 }
    01428
    01429 template<class sq_T>
    01430 shared_ptr<mpdf> enorm<sq_T>::condition ( const RV &rvn ) const
    01431 {
    01432         mlnorm<sq_T> *tmp = new mlnorm<sq_T> ();
    01433         shared_ptr<mpdf> narrow(tmp);
    01434         condition ( rvn, *tmp );
    01435         return narrow;
    01436 }
    01437
    01438 template<class sq_T>
    01439 void enorm<sq_T>::condition ( const RV &rvn, mpdf &target ) const
    01440 {
    01441         typedef mlnorm<sq_T> TMlnorm;
    01442
    01443         bdm_assert (isnamed(), "rvs are not assigned");
    01444         TMlnorm &uptarget = dynamic_cast<TMlnorm &>(target);
    01445
    01446         RV rvc = rv.subt (rvn);
    01447         bdm_assert ( (rvc._dsize() + rvn._dsize() == rv._dsize()), "wrong rvn");
    01448         //Permutation vector of the new R
    01449         ivec irvn = rvn.dataind (rv);
    01450         ivec irvc = rvc.dataind (rv);
    01451         ivec perm = concat (irvn , irvc);
    01452         sq_T Rn (R, perm);
    01453
    01454         //fixme - could this be done in general for all sq_T?
    01455         mat S = Rn.to_mat();
    01456         //fixme
    01457         int n = rvn._dsize() - 1;
    01458         int end = R.rows() - 1;
    01459         mat S11 = S.get (0, n, 0, n);
    01460         mat S12 = S.get (0, n , rvn._dsize(), end);
    01461         mat S22 = S.get (rvn._dsize(), end, rvn._dsize(), end);
    01462
    01463         vec mu1 = mu (irvn);
    01464         vec mu2 = mu (irvc);
    01465         mat A = S12 * inv (S22);
    01466         sq_T R_n (S11 - A *S12.T());
    01467
    01468         uptarget.set_rv (rvn);
    01469         uptarget.set_rvc (rvc);
    01470         uptarget.set_parameters (A, mu1 - A*mu2, R_n);
    01471 }
    01472
    01475 template<class sq_T>
    01476 void mgnorm<sq_T >::set_parameters (const shared_ptr<fnc> &g0, const sq_T &R0) {
    01477         g = g0;
    01478         this->iepdf.set_parameters (zeros (g->dimension()), R0);
    01479 }
    01480
    01481 template<class sq_T>
    01482 void mgnorm<sq_T >::condition (const vec &cond) {this->iepdf._mu() = g->eval (cond);};
    01483
    01485 template<class sq_T>
    01486 std::ostream &operator<< (std::ostream &os,  mlnorm<sq_T> &ml)
    01487 {
    01488         os << "A:" << ml.A << endl;
    01489         os << "mu:" << ml.mu_const << endl;
    01490         os << "R:" << ml._R() << endl;
    01491         return os;
    01492 };
    01493
    01494 }
    01495 #endif //EF_H