- 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;
    00047                 virtual double evallog_nn (const vec &val) const{it_error ("Not implemented");return 0.0;};
    00049                 virtual double evallog (const vec &val) const {
    00050                         double tmp;
    00051                         tmp = evallog_nn (val) - lognc();
    00052 //                              it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
    00053                         return tmp;
    00054                 }
    00056                 virtual vec evallog_m (const mat &Val) const {
    00057                         vec x (Val.cols());
    00058                         for (int i = 0;i < Val.cols();i++) {x (i) = evallog_nn (Val.get_col (i)) ;}
    00059                         return x -lognc();
    00060                 }
    00062                 virtual vec evallog_m (const Array<vec> &Val) const {
    00063                         vec x (Val.length());
    00064                         for (int i = 0;i < Val.length();i++) {x (i) = evallog_nn (Val (i)) ;}
    00065                         return x -lognc();
    00066                 }
    00068                 virtual void pow (double p) {it_error ("Not implemented");};
    00069 };
    00070
    00071
    00073 class BMEF : public BM
    00074 {
    00075         protected:
    00077                 double frg;
    00079                 double last_lognc;
    00080         public:
    00082                 BMEF (double frg0 = 1.0) : BM (), frg (frg0) {}
    00084                 BMEF (const BMEF &B) : BM (B), frg (B.frg), last_lognc (B.last_lognc) {}
    00086                 virtual void set_statistics (const BMEF* BM0) {it_error ("Not implemented");};
    00088                 virtual void bayes (const vec &data, const double w) {};
    00089                 //original Bayes
    00090                 void bayes (const vec &dt);
    00092                 virtual void flatten (const BMEF * B) {it_error ("Not implemented");}
    00094 //      virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
    00095
    00096                 BMEF* _copy_ () const {it_error ("function _copy_ not implemented for this BM"); return NULL;};
    00097 };
    00098
    00099 template<class sq_T, template <typename> class TEpdf>
    00100 class mlnorm;
    00101
    00107 template<class sq_T>
    00108 class enorm : public eEF
    00109 {
    00110         protected:
    00112                 vec mu;
    00114                 sq_T R;
    00115         public:
    00118
    00119                 enorm () : eEF (), mu (), R () {};
    00120                 enorm (const vec &mu, const sq_T &R) {set_parameters (mu, R);}
    00121                 void set_parameters (const vec &mu, const sq_T &R);
    00122                 void from_setting (const Setting &root);
    00123                 void validate() {
    00124                         it_assert (mu.length() == R.rows(), "parameters mismatch");
    00125                         dim = mu.length();
    00126                 }
    00128
    00131
    00133                 void dupdate (mat &v, double nu = 1.0);
    00134
    00135                 vec sample() const;
    00136
    00137                 double evallog_nn (const vec &val) const;
    00138                 double lognc () const;
    00139                 vec mean() const {return mu;}
    00140                 vec variance() const {return diag (R.to_mat());}
    00141 //      mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
    00142                 shared_ptr<mpdf> condition ( const RV &rvn ) const;
    00143
    00144                 // target not typed to mlnorm<sq_T, enorm<sq_T> > &
    00145                 // because that doesn't compile (perhaps because we
    00146                 // haven't finished defining enorm yet), but the type
    00147                 // is required
    00148                 void condition ( const RV &rvn, mpdf &target ) const;
    00149
    00150                 shared_ptr<epdf> marginal (const RV &rvn ) const;
    00151                 void marginal ( const RV &rvn, enorm<sq_T> &target ) const;
    00153
    00156
    00157                 vec& _mu() {return mu;}
    00158                 void set_mu (const vec mu0) { mu = mu0;}
    00159                 sq_T& _R() {return R;}
    00160                 const sq_T& _R() const {return R;}
    00162
    00163 };
    00164 UIREGISTER2 (enorm, chmat);
    00165 SHAREDPTR2 ( enorm, chmat );
    00166 UIREGISTER2 (enorm, ldmat);
    00167 SHAREDPTR2 ( enorm, ldmat );
    00168 UIREGISTER2 (enorm, fsqmat);
    00169 SHAREDPTR2 ( enorm, fsqmat );
    00170
    00171
    00178 class egiw : public eEF
    00179 {
    00180         protected:
    00182                 ldmat V;
    00184                 double nu;
    00186                 int dimx;
    00188                 int nPsi;
    00189         public:
    00192                 egiw() : eEF() {};
    00193                 egiw (int dimx0, ldmat V0, double nu0 = -1.0) : eEF() {set_parameters (dimx0, V0, nu0);};
    00194
    00195                 void set_parameters (int dimx0, ldmat V0, double nu0 = -1.0) {
    00196                         dimx = dimx0;
    00197                         nPsi = V0.rows() - dimx;
    00198                         dim = dimx * (dimx + nPsi); // size(R) + size(Theta)
    00199
    00200                         V = V0;
    00201                         if (nu0 < 0) {
    00202                                 nu = 0.1 + nPsi + 2 * dimx + 2; // +2 assures finite expected value of R
    00203                                 // terms before that are sufficient for finite normalization
    00204                         } else {
    00205                                 nu = nu0;
    00206                         }
    00207                 }
    00209
    00210                 vec sample() const;
    00211                 vec mean() const;
    00212                 vec variance() const;
    00213
    00215                 vec est_theta() const;
    00216
    00218                 ldmat est_theta_cov() const;
    00219
    00221                 void mean_mat (mat &M, mat&R) const;
    00223                 double evallog_nn (const vec &val) const;
    00224                 double lognc () const;
    00225                 void pow (double p) {V *= p;nu *= p;};
    00226
    00229
    00230                 ldmat& _V() {return V;}
    00231                 const ldmat& _V() const {return V;}
    00232                 double& _nu()  {return nu;}
    00233                 const double& _nu() const {return nu;}
    00234                 void from_setting (const Setting &set) {
    00235                         UI::get (nu, set, "nu", UI::compulsory);
    00236                         UI::get (dimx, set, "dimx", UI::compulsory);
    00237                         mat V;
    00238                         UI::get (V, set, "V", UI::compulsory);
    00239                         set_parameters (dimx, V, nu);
    00240                         shared_ptr<RV> rv = UI::build<RV> (set, "rv", UI::compulsory);
    00241                         set_rv (*rv);
    00242                 }
    00244 };
    00245 UIREGISTER ( egiw );
    00246 SHAREDPTR ( egiw );
    00247
    00256 class eDirich: public eEF
    00257 {
    00258         protected:
    00260                 vec beta;
    00261         public:
    00264
    00265                 eDirich () : eEF () {};
    00266                 eDirich (const eDirich &D0) : eEF () {set_parameters (D0.beta);};
    00267                 eDirich (const vec &beta0) {set_parameters (beta0);};
    00268                 void set_parameters (const vec &beta0) {
    00269                         beta = beta0;
    00270                         dim = beta.length();
    00271                 }
    00273
    00274                 vec sample() const {it_error ("Not implemented");return vec_1 (0.0);};
    00275                 vec mean() const {return beta / sum (beta);};
    00276                 vec variance() const {double gamma = sum (beta); return elem_mult (beta, (beta + 1)) / (gamma* (gamma + 1));}
    00278                 double evallog_nn (const vec &val) const {
    00279                         double tmp; tmp = (beta - 1) * log (val);
    00280 //                              it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
    00281                         return tmp;
    00282                 };
    00283                 double lognc () const {
    00284                         double tmp;
    00285                         double gam = sum (beta);
    00286                         double lgb = 0.0;
    00287                         for (int i = 0;i < beta.length();i++) {lgb += lgamma (beta (i));}
    00288                         tmp = lgb - lgamma (gam);
    00289 //                              it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
    00290                         return tmp;
    00291                 };
    00293                 vec& _beta()  {return beta;}
    00295 };
    00296
    00298 class multiBM : public BMEF
    00299 {
    00300         protected:
    00302                 eDirich est;
    00304                 vec &beta;
    00305         public:
    00307                 multiBM () : BMEF (), est (), beta (est._beta()) {
    00308                         if (beta.length() > 0) {last_lognc = est.lognc();}
    00309                         else{last_lognc = 0.0;}
    00310                 }
    00312                 multiBM (const multiBM &B) : BMEF (B), est (B.est), beta (est._beta()) {}
    00314                 void set_statistics (const BM* mB0) {const multiBM* mB = dynamic_cast<const multiBM*> (mB0); beta = mB->beta;}
    00315                 void bayes (const vec &dt) {
    00316                         if (frg < 1.0) {beta *= frg;last_lognc = est.lognc();}
    00317                         beta += dt;
    00318                         if (evalll) {ll = est.lognc() - last_lognc;}
    00319                 }
    00320                 double logpred (const vec &dt) const {
    00321                         eDirich pred (est);
    00322                         vec &beta = pred._beta();
    00323
    00324                         double lll;
    00325                         if (frg < 1.0)
    00326                                 {beta *= frg;lll = pred.lognc();}
    00327                         else
    00328                                 if (evalll) {lll = last_lognc;}
    00329                                 else{lll = pred.lognc();}
    00330
    00331                         beta += dt;
    00332                         return pred.lognc() - lll;
    00333                 }
    00334                 void flatten (const BMEF* B) {
    00335                         const multiBM* E = dynamic_cast<const multiBM*> (B);
    00336                         // sum(beta) should be equal to sum(B.beta)
    00337                         const vec &Eb = E->beta;//const_cast<multiBM*> ( E )->_beta();
    00338                         beta *= (sum (Eb) / sum (beta));
    00339                         if (evalll) {last_lognc = est.lognc();}
    00340                 }
    00342                 const eDirich& posterior() const {return est;};
    00344                 void set_parameters (const vec &beta0) {
    00345                         est.set_parameters (beta0);
    00346                         if (evalll) {last_lognc = est.lognc();}
    00347                 }
    00348 };
    00349
    00359 class egamma : public eEF
    00360 {
    00361         protected:
    00363                 vec alpha;
    00365                 vec beta;
    00366         public :
    00369                 egamma () : eEF (), alpha (0), beta (0) {};
    00370                 egamma (const vec &a, const vec &b) {set_parameters (a, b);};
    00371                 void set_parameters (const vec &a, const vec &b) {alpha = a, beta = b;dim = alpha.length();};
    00373
    00374                 vec sample() const;
    00375                 double evallog (const vec &val) const;
    00376                 double lognc () const;
    00378                 vec& _alpha() {return alpha;}
    00380                 vec& _beta() {return beta;}
    00381                 vec mean() const {return elem_div (alpha, beta);}
    00382                 vec variance() const {return elem_div (alpha, elem_mult (beta, beta)); }
    00383
    00392                 void from_setting (const Setting &set) {
    00393                         epdf::from_setting (set); // reads rv
    00394                         UI::get (alpha, set, "alpha", UI::compulsory);
    00395                         UI::get (beta, set, "beta", UI::compulsory);
    00396                         validate();
    00397                 }
    00398                 void validate() {
    00399                         it_assert (alpha.length() == beta.length(), "parameters do not match");
    00400                         dim = alpha.length();
    00401                 }
    00402 };
    00403 UIREGISTER (egamma);
    00404 SHAREDPTR ( egamma );
    00405
    00422 class eigamma : public egamma
    00423 {
    00424         protected:
    00425         public :
    00430
    00431                 vec sample() const {return 1.0 / egamma::sample();};
    00433                 vec mean() const {return elem_div (beta, alpha - 1);}
    00434                 vec variance() const {vec mea = mean(); return elem_div (elem_mult (mea, mea), alpha - 2);}
    00435 };
    00436 /*
    00438 class emix : public epdf {
    00439 protected:
    00440         int n;
    00441         vec &w;
    00442         Array<epdf*> Coms;
    00443 public:
    00445         emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
    00446         void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
    00447         vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
    00448         vec sample() {it_error ( "Not implemented" );return 0;}
    00449 };
    00450 */
    00451
    00453
    00454 class euni: public epdf
    00455 {
    00456         protected:
    00458                 vec low;
    00460                 vec high;
    00462                 vec distance;
    00464                 double nk;
    00466                 double lnk;
    00467         public:
    00470                 euni () : epdf () {}
    00471                 euni (const vec &low0, const vec &high0) {set_parameters (low0, high0);}
    00472                 void set_parameters (const vec &low0, const vec &high0) {
    00473                         distance = high0 - low0;
    00474                         it_assert_debug (min (distance) > 0.0, "bad support");
    00475                         low = low0;
    00476                         high = high0;
    00477                         nk = prod (1.0 / distance);
    00478                         lnk = log (nk);
    00479                         dim = low.length();
    00480                 }
    00482
    00483                 double evallog (const vec &val) const  {
    00484                         if (any (val < low) && any (val > high)) {return inf;}
    00485                         else return lnk;
    00486                 }
    00487                 vec sample() const {
    00488                         vec smp (dim);
    00489 #pragma omp critical
    00490                         UniRNG.sample_vector (dim , smp);
    00491                         return low + elem_mult (distance, smp);
    00492                 }
    00494                 vec mean() const {return (high -low) / 2.0;}
    00495                 vec variance() const {return (pow (high, 2) + pow (low, 2) + elem_mult (high, low)) / 3.0;}
    00504                 void from_setting (const Setting &set) {
    00505                         epdf::from_setting (set); // reads rv and rvc
    00506
    00507                         UI::get (high, set, "high", UI::compulsory);
    00508                         UI::get (low, set, "low", UI::compulsory);
    00509                 }
    00510 };
    00511
    00512
    00518 template < class sq_T, template <typename> class TEpdf = enorm >
    00519 class mlnorm : public mpdf_internal< TEpdf<sq_T> >
    00520 {
    00521         protected:
    00523                 mat A;
    00525                 vec mu_const;
    00526 //                      vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
    00527         public:
    00530                 mlnorm() : mpdf_internal< TEpdf<sq_T> >() {};
    00531                 mlnorm (const mat &A, const vec &mu0, const sq_T &R) : mpdf_internal< TEpdf<sq_T> >() {
    00532                         set_parameters (A, mu0, R);
    00533                 }
    00534
    00536                 void set_parameters (const  mat &A0, const vec &mu0, const sq_T &R0) {
    00537                         it_assert_debug (A0.rows() == mu0.length(), "");
    00538                         it_assert_debug (A0.rows() == R0.rows(), "");
    00539
    00540                         this->iepdf.set_parameters (zeros (A0.rows()), R0);
    00541                         A = A0;
    00542                         mu_const = mu0;
    00543                         this->dimc = A0.cols();
    00544                 }
    00547                 void condition (const vec &cond) {
    00548                         this->iepdf._mu() = A * cond + mu_const;
    00549 //R is already assigned;
    00550                 }
    00551
    00553                 vec& _mu_const() {return mu_const;}
    00555                 mat& _A() {return A;}
    00557                 mat _R() { return this->iepdf._R().to_mat(); }
    00558
    00560                 template<typename sq_M>
    00561                 friend std::ostream &operator<< (std::ostream &os,  mlnorm<sq_M, enorm> &ml);
    00562
    00563                 void from_setting (const Setting &set) {
    00564                         mpdf::from_setting (set);
    00565
    00566                         UI::get (A, set, "A", UI::compulsory);
    00567                         UI::get (mu_const, set, "const", UI::compulsory);
    00568                         mat R0;
    00569                         UI::get (R0, set, "R", UI::compulsory);
    00570                         set_parameters (A, mu_const, R0);
    00571                 };
    00572 };
    00573 UIREGISTER2 (mlnorm,ldmat);
    00574 SHAREDPTR2 ( mlnorm, ldmat );
    00575 UIREGISTER2 (mlnorm,fsqmat);
    00576 SHAREDPTR2 ( mlnorm, fsqmat );
    00577 UIREGISTER2 (mlnorm, chmat);
    00578 SHAREDPTR2 ( mlnorm, chmat );
    00579
    00581 template<class sq_T>
    00582 class mgnorm : public mpdf_internal< enorm< sq_T > >
    00583 {
    00584         private:
    00585 //                      vec &mu; WHY NOT?
    00586                 shared_ptr<fnc> g;
    00587
    00588         public:
    00590                 mgnorm() : mpdf_internal<enorm<sq_T> >() { }
    00592                 inline void set_parameters (const shared_ptr<fnc> &g0, const sq_T &R0);
    00593                 inline void condition (const vec &cond);
    00594
    00595
    00623                 void from_setting (const Setting &set) {
    00624                         shared_ptr<fnc> g = UI::build<fnc> (set, "g", UI::compulsory);
    00625
    00626                         mat R;
    00627                         vec dR;
    00628                         if (UI::get (dR, set, "dR"))
    00629                                 R = diag (dR);
    00630                         else
    00631                                 UI::get (R, set, "R", UI::compulsory);
    00632
    00633                         set_parameters (g, R);
    00634                 }
    00635 };
    00636
    00637 UIREGISTER2 (mgnorm, chmat);
    00638 SHAREDPTR2 ( mgnorm, chmat );
    00639
    00640
    00648 class mlstudent : public mlnorm<ldmat, enorm>
    00649 {
    00650         protected:
    00652                 ldmat Lambda;
    00654                 ldmat &_R;
    00656                 ldmat Re;
    00657         public:
    00658                 mlstudent () : mlnorm<ldmat, enorm> (),
    00659                                 Lambda (),      _R (iepdf._R()) {}
    00661                 void set_parameters (const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0) {
    00662                         it_assert_debug (A0.rows() == mu0.length(), "");
    00663                         it_assert_debug (R0.rows() == A0.rows(), "");
    00664
    00665                         iepdf.set_parameters (mu0, Lambda);  //
    00666                         A = A0;
    00667                         mu_const = mu0;
    00668                         Re = R0;
    00669                         Lambda = Lambda0;
    00670                 }
    00671                 void condition (const vec &cond) {
    00672                         iepdf._mu() = A * cond + mu_const;
    00673                         double zeta;
    00674                         //ugly hack!
    00675                         if ( (cond.length() + 1) == Lambda.rows()) {
    00676                                 zeta = Lambda.invqform (concat (cond, vec_1 (1.0)));
    00677                         } else {
    00678                                 zeta = Lambda.invqform (cond);
    00679                         }
    00680                         _R = Re;
    00681                         _R *= (1 + zeta);// / ( nu ); << nu is in Re!!!!!!
    00682                 };
    00683
    00684 };
    00694 class mgamma : public mpdf_internal<egamma>
    00695 {
    00696         protected:
    00697
    00699                 double k;
    00700
    00702                 vec &_beta;
    00703
    00704         public:
    00706                 mgamma() : mpdf_internal<egamma>(), k (0),
    00707                                 _beta (iepdf._beta()) {
    00708                 }
    00709
    00711                 void set_parameters (double k, const vec &beta0);
    00712
    00713                 void condition (const vec &val) {_beta = k / val;};
    00723                 void from_setting (const Setting &set) {
    00724                         mpdf::from_setting (set); // reads rv and rvc
    00725                         vec betatmp; // ugly but necessary
    00726                         UI::get (betatmp, set, "beta", UI::compulsory);
    00727                         UI::get (k, set, "k", UI::compulsory);
    00728                         set_parameters (k, betatmp);
    00729                 }
    00730 };
    00731 UIREGISTER (mgamma);
    00732 SHAREDPTR (mgamma);
    00733
    00743 class migamma : public mpdf_internal<eigamma>
    00744 {
    00745         protected:
    00747                 double k;
    00748
    00750                 vec &_alpha;
    00751
    00753                 vec &_beta;
    00754
    00755         public:
    00758                 migamma() : mpdf_internal<eigamma>(),
    00759                                 k (0),
    00760                                 _alpha (iepdf._alpha()),
    00761                                 _beta (iepdf._beta()) {
    00762                 }
    00763
    00764                 migamma (const migamma &m) : mpdf_internal<eigamma>(),
    00765                                 k (0),
    00766                                 _alpha (iepdf._alpha()),
    00767                                 _beta (iepdf._beta()) {
    00768                 }
    00770
    00772                 void set_parameters (int len, double k0) {
    00773                         k = k0;
    00774                         iepdf.set_parameters ( (1.0 / (k*k) + 2.0) *ones (len) /*alpha*/, ones (len) /*beta*/);
    00775                         dimc = dimension();
    00776                 };
    00777                 void condition (const vec &val) {
    00778                         _beta = elem_mult (val, (_alpha - 1.0));
    00779                 };
    00780 };
    00781
    00782
    00794 class mgamma_fix : public mgamma
    00795 {
    00796         protected:
    00798                 double l;
    00800                 vec refl;
    00801         public:
    00803                 mgamma_fix () : mgamma (), refl () {};
    00805                 void set_parameters (double k0 , vec ref0, double l0) {
    00806                         mgamma::set_parameters (k0, ref0);
    00807                         refl = pow (ref0, 1.0 - l0);l = l0;
    00808                         dimc = dimension();
    00809                 };
    00810
    00811                 void condition (const vec &val) {vec mean = elem_mult (refl, pow (val, l)); _beta = k / mean;};
    00812 };
    00813
    00814
    00827 class migamma_ref : public migamma
    00828 {
    00829         protected:
    00831                 double l;
    00833                 vec refl;
    00834         public:
    00836                 migamma_ref () : migamma (), refl () {};
    00838                 void set_parameters (double k0 , vec ref0, double l0) {
    00839                         migamma::set_parameters (ref0.length(), k0);
    00840                         refl = pow (ref0, 1.0 - l0);
    00841                         l = l0;
    00842                         dimc = dimension();
    00843                 };
    00844
    00845                 void condition (const vec &val) {
    00846                         vec mean = elem_mult (refl, pow (val, l));
    00847                         migamma::condition (mean);
    00848                 };
    00849
    00870                 void from_setting (const Setting &set);
    00871
    00872                 // TODO dodelat void to_setting( Setting &set ) const;
    00873 };
    00874
    00875
    00876 UIREGISTER (migamma_ref);
    00877 SHAREDPTR (migamma_ref);
    00878
    00888 class elognorm: public enorm<ldmat>
    00889 {
    00890         public:
    00891                 vec sample() const {return exp (enorm<ldmat>::sample());};
    00892                 vec mean() const {vec var = enorm<ldmat>::variance();return exp (mu - 0.5*var);};
    00893
    00894 };
    00895
    00907 class mlognorm : public mpdf_internal<elognorm>
    00908 {
    00909         protected:
    00911                 double sig2;
    00912
    00914                 vec &mu;
    00915         public:
    00917                 mlognorm() : mpdf_internal<elognorm>(),
    00918                                 sig2 (0),
    00919                                 mu (iepdf._mu()) {
    00920                 }
    00921
    00923                 void set_parameters (int size, double k) {
    00924                         sig2 = 0.5 * log (k * k + 1);
    00925                         iepdf.set_parameters (zeros (size), 2*sig2*eye (size));
    00926
    00927                         dimc = size;
    00928                 };
    00929
    00930                 void condition (const vec &val) {
    00931                         mu = log (val) - sig2;//elem_mult ( refl,pow ( val,l ) );
    00932                 };
    00933
    00952                 void from_setting (const Setting &set);
    00953
    00954                 // TODO dodelat void to_setting( Setting &set ) const;
    00955
    00956 };
    00957
    00958 UIREGISTER (mlognorm);
    00959 SHAREDPTR (mlognorm);
    00960
    00964 class eWishartCh : public epdf
    00965 {
    00966         protected:
    00968                 chmat Y;
    00970                 int p;
    00972                 double delta;
    00973         public:
    00975                 void set_parameters (const mat &Y0, const double delta0) {Y = chmat (Y0);delta = delta0; p = Y.rows(); dim = p * p; }
    00977                 mat sample_mat() const {
    00978                         mat X = zeros (p, p);
    00979
    00980                         //sample diagonal
    00981                         for (int i = 0;i < p;i++) {
    00982                                 GamRNG.setup (0.5* (delta - i) , 0.5);   // no +1 !! index if from 0
    00983 #pragma omp critical
    00984                                 X (i, i) = sqrt (GamRNG());
    00985                         }
    00986                         //do the rest
    00987                         for (int i = 0;i < p;i++) {
    00988                                 for (int j = i + 1;j < p;j++) {
    00989 #pragma omp critical
    00990                                         X (i, j) = NorRNG.sample();
    00991                                 }
    00992                         }
    00993                         return X*Y._Ch();// return upper triangular part of the decomposition
    00994                 }
    00995                 vec sample () const {
    00996                         return vec (sample_mat()._data(), p*p);
    00997                 }
    00999                 void setY (const mat &Ch0) {copy_vector (dim, Ch0._data(), Y._Ch()._data());}
    01001                 void _setY (const vec &ch0) {copy_vector (dim, ch0._data(), Y._Ch()._data()); }
    01003                 const chmat& getY() const {return Y;}
    01004 };
    01005
    01007
    01009 class eiWishartCh: public epdf
    01010 {
    01011         protected:
    01013                 eWishartCh W;
    01015                 int p;
    01017                 double delta;
    01018         public:
    01020                 void set_parameters (const mat &Y0, const double delta0) {
    01021                         delta = delta0;
    01022                         W.set_parameters (inv (Y0), delta0);
    01023                         dim = W.dimension(); p = Y0.rows();
    01024                 }
    01025                 vec sample() const {mat iCh; iCh = inv (W.sample_mat()); return vec (iCh._data(), dim);}
    01027                 void _setY (const vec &y0) {
    01028                         mat Ch (p, p);
    01029                         mat iCh (p, p);
    01030                         copy_vector (dim, y0._data(), Ch._data());
    01031
    01032                         iCh = inv (Ch);
    01033                         W.setY (iCh);
    01034                 }
    01035                 virtual double evallog (const vec &val) const {
    01036                         chmat X (p);
    01037                         const chmat& Y = W.getY();
    01038
    01039                         copy_vector (p*p, val._data(), X._Ch()._data());
    01040                         chmat iX (p);X.inv (iX);
    01041                         // compute
    01042 //                              \frac{ |\Psi|^{m/2}|X|^{-(m+p+1)/2}e^{-tr(\Psi X^{-1})/2} }{ 2^{mp/2}\Gamma_p(m/2)},
    01043                         mat M = Y.to_mat() * iX.to_mat();
    01044
    01045                         double log1 = 0.5 * p * (2 * Y.logdet()) - 0.5 * (delta + p + 1) * (2 * X.logdet()) - 0.5 * trace (M);
    01046                         //Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!
    01047
    01048                         /*                              if (0) {
    01049                                                                 mat XX=X.to_mat();
    01050                                                                 mat YY=Y.to_mat();
    01051
    01052                                                                 double log2 = 0.5*p*log(det(YY))-0.5*(delta+p+1)*log(det(XX))-0.5*trace(YY*inv(XX));
    01053                                                                 cout << log1 << "," << log2 << endl;
    01054                                                         }*/
    01055                         return log1;
    01056                 };
    01057
    01058 };
    01059
    01061 class rwiWishartCh : public mpdf_internal<eiWishartCh>
    01062 {
    01063         protected:
    01065                 double sqd;
    01067                 vec refl;
    01069                 double l;
    01071                 int p;
    01072
    01073         public:
    01074                 rwiWishartCh() : sqd (0), l (0), p (0) {}
    01076                 void set_parameters (int p0, double k, vec ref0, double l0) {
    01077                         p = p0;
    01078                         double delta = 2 / (k * k) + p + 3;
    01079                         sqd = sqrt (delta - p - 1);
    01080                         l = l0;
    01081                         refl = pow (ref0, 1 - l);
    01082
    01083                         iepdf.set_parameters (eye (p), delta);
    01084                         dimc = iepdf.dimension();
    01085                 }
    01086                 void condition (const vec &c) {
    01087                         vec z = c;
    01088                         int ri = 0;
    01089                         for (int i = 0;i < p*p;i += (p + 1)) {//trace diagonal element
    01090                                 z (i) = pow (z (i), l) * refl (ri);
    01091                                 ri++;
    01092                         }
    01093
    01094                         iepdf._setY (sqd*z);
    01095                 }
    01096 };
    01097
    01099 enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
    01105 class eEmp: public epdf
    01106 {
    01107         protected :
    01109                 int n;
    01111                 vec w;
    01113                 Array<vec> samples;
    01114         public:
    01117                 eEmp () : epdf (), w (), samples () {};
    01119                 eEmp (const eEmp &e) : epdf (e), w (e.w), samples (e.samples) {};
    01121
    01123                 void set_statistics (const vec &w0, const epdf &pdf0);
    01125                 void set_statistics (const epdf &pdf0 , int n) {set_statistics (ones (n) / n, pdf0);};
    01127                 void set_samples (const epdf* pdf0);
    01129                 void set_parameters (int n0, bool copy = true) {n = n0; w.set_size (n0, copy);samples.set_size (n0, copy);};
    01131                 vec& _w()  {return w;};
    01133                 const vec& _w() const {return w;};
    01135                 Array<vec>& _samples() {return samples;};
    01137                 const Array<vec>& _samples() const {return samples;};
    01139                 ivec resample (RESAMPLING_METHOD method = SYSTEMATIC);
    01141                 vec sample() const {it_error ("Not implemented");return 0;}
    01143                 double evallog (const vec &val) const {it_error ("Not implemented");return 0.0;}
    01144                 vec mean() const {
    01145                         vec pom = zeros (dim);
    01146                         for (int i = 0;i < n;i++) {pom += samples (i) * w (i);}
    01147                         return pom;
    01148                 }
    01149                 vec variance() const {
    01150                         vec pom = zeros (dim);
    01151                         for (int i = 0;i < n;i++) {pom += pow (samples (i), 2) * w (i);}
    01152                         return pom -pow (mean(), 2);
    01153                 }
    01155                 void qbounds (vec &lb, vec &ub, double perc = 0.95) const {
    01156                         // lb in inf so than it will be pushed below;
    01157                         lb.set_size (dim);
    01158                         ub.set_size (dim);
    01159                         lb = std::numeric_limits<double>::infinity();
    01160                         ub = -std::numeric_limits<double>::infinity();
    01161                         int j;
    01162                         for (int i = 0;i < n;i++) {
    01163                                 for (j = 0;j < dim; j++) {
    01164                                         if (samples (i) (j) < lb (j)) {lb (j) = samples (i) (j);}
    01165                                         if (samples (i) (j) > ub (j)) {ub (j) = samples (i) (j);}
    01166                                 }
    01167                         }
    01168                 }
    01169 };
    01170
    01171
    01173
    01174 template<class sq_T>
    01175 void enorm<sq_T>::set_parameters (const vec &mu0, const sq_T &R0)
    01176 {
    01177 //Fixme test dimensions of mu0 and R0;
    01178         mu = mu0;
    01179         R = R0;
    01180         validate();
    01181 };
    01182
    01183 template<class sq_T>
    01184 void enorm<sq_T>::from_setting (const Setting &set)
    01185 {
    01186         epdf::from_setting (set); //reads rv
    01187
    01188         UI::get (mu, set, "mu", UI::compulsory);
    01189         mat Rtmp;// necessary for conversion
    01190         UI::get (Rtmp, set, "R", UI::compulsory);
    01191         R = Rtmp; // conversion
    01192         validate();
    01193 }
    01194
    01195 template<class sq_T>
    01196 void enorm<sq_T>::dupdate (mat &v, double nu)
    01197 {
    01198         //
    01199 };
    01200
    01201 // template<class sq_T>
    01202 // void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
    01203 //      //
    01204 // };
    01205
    01206 template<class sq_T>
    01207 vec enorm<sq_T>::sample() const
    01208 {
    01209         vec x (dim);
    01210 #pragma omp critical
    01211         NorRNG.sample_vector (dim, x);
    01212         vec smp = R.sqrt_mult (x);
    01213
    01214         smp += mu;
    01215         return smp;
    01216 };
    01217
    01218 // template<class sq_T>
    01219 // double enorm<sq_T>::eval ( const vec &val ) const {
    01220 //      double pdfl,e;
    01221 //      pdfl = evallog ( val );
    01222 //      e = exp ( pdfl );
    01223 //      return e;
    01224 // };
    01225
    01226 template<class sq_T>
    01227 double enorm<sq_T>::evallog_nn (const vec &val) const
    01228 {
    01229         // 1.83787706640935 = log(2pi)
    01230         double tmp = -0.5 * (R.invqform (mu - val));// - lognc();
    01231         return  tmp;
    01232 };
    01233
    01234 template<class sq_T>
    01235 inline double enorm<sq_T>::lognc () const
    01236 {
    01237         // 1.83787706640935 = log(2pi)
    01238         double tmp = 0.5 * (R.cols() * 1.83787706640935 + R.logdet());
    01239         return tmp;
    01240 };
    01241
    01242
    01243 // template<class sq_T>
    01244 // vec mlnorm<sq_T>::samplecond (const  vec &cond, double &lik ) {
    01245 //      this->condition ( cond );
    01246 //      vec smp = epdf.sample();
    01247 //      lik = epdf.eval ( smp );
    01248 //      return smp;
    01249 // }
    01250
    01251 // template<class sq_T>
    01252 // mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) {
    01253 //      int i;
    01254 //      int dim = rv.count();
    01255 //      mat Smp ( dim,n );
    01256 //      vec smp ( dim );
    01257 //      this->condition ( cond );
    01258 //
    01259 //      for ( i=0; i<n; i++ ) {
    01260 //              smp = epdf.sample();
    01261 //              lik ( i ) = epdf.eval ( smp );
    01262 //              Smp.set_col ( i ,smp );
    01263 //      }
    01264 //
    01265 //      return Smp;
    01266 // }
    01267
    01268
    01269 template<class sq_T>
    01270 shared_ptr<epdf> enorm<sq_T>::marginal ( const RV &rvn ) const
    01271 {
    01272         enorm<sq_T> *tmp = new enorm<sq_T> ();
    01273         shared_ptr<epdf> narrow(tmp);
    01274         marginal ( rvn, *tmp );
    01275         return narrow;
    01276 }
    01277
    01278 template<class sq_T>
    01279 void enorm<sq_T>::marginal ( const RV &rvn, enorm<sq_T> &target ) const
    01280 {
    01281         it_assert_debug (isnamed(), "rv description is not assigned");
    01282         ivec irvn = rvn.dataind (rv);
    01283
    01284         sq_T Rn (R, irvn);  // select rows and columns of R
    01285
    01286         target.set_rv ( rvn );
    01287         target.set_parameters (mu (irvn), Rn);
    01288 }
    01289
    01290 template<class sq_T>
    01291 shared_ptr<mpdf> enorm<sq_T>::condition ( const RV &rvn ) const
    01292 {
    01293         mlnorm<sq_T> *tmp = new mlnorm<sq_T> ();
    01294         shared_ptr<mpdf> narrow(tmp);
    01295         condition ( rvn, *tmp );
    01296         return narrow;
    01297 }
    01298
    01299 template<class sq_T>
    01300 void enorm<sq_T>::condition ( const RV &rvn, mpdf &target ) const
    01301 {
    01302         typedef mlnorm<sq_T> TMlnorm;
    01303
    01304         it_assert_debug (isnamed(), "rvs are not assigned");
    01305         TMlnorm &uptarget = dynamic_cast<TMlnorm &>(target);
    01306
    01307         RV rvc = rv.subt (rvn);
    01308         it_assert_debug ( (rvc._dsize() + rvn._dsize() == rv._dsize()), "wrong rvn");
    01309         //Permutation vector of the new R
    01310         ivec irvn = rvn.dataind (rv);
    01311         ivec irvc = rvc.dataind (rv);
    01312         ivec perm = concat (irvn , irvc);
    01313         sq_T Rn (R, perm);
    01314
    01315         //fixme - could this be done in general for all sq_T?
    01316         mat S = Rn.to_mat();
    01317         //fixme
    01318         int n = rvn._dsize() - 1;
    01319         int end = R.rows() - 1;
    01320         mat S11 = S.get (0, n, 0, n);
    01321         mat S12 = S.get (0, n , rvn._dsize(), end);
    01322         mat S22 = S.get (rvn._dsize(), end, rvn._dsize(), end);
    01323
    01324         vec mu1 = mu (irvn);
    01325         vec mu2 = mu (irvc);
    01326         mat A = S12 * inv (S22);
    01327         sq_T R_n (S11 - A *S12.T());
    01328
    01329         uptarget.set_rv (rvn);
    01330         uptarget.set_rvc (rvc);
    01331         uptarget.set_parameters (A, mu1 - A*mu2, R_n);
    01332 }
    01333
    01336 template<class sq_T>
    01337 void mgnorm<sq_T >::set_parameters (const shared_ptr<fnc> &g0, const sq_T &R0) {
    01338         g = g0;
    01339         this->iepdf.set_parameters (zeros (g->dimension()), R0);
    01340 }
    01341
    01342 template<class sq_T>
    01343 void mgnorm<sq_T >::condition (const vec &cond) {this->iepdf._mu() = g->eval (cond);};
    01344
    01346 template<class sq_T>
    01347 std::ostream &operator<< (std::ostream &os,  mlnorm<sq_T> &ml)
    01348 {
    01349         os << "A:" << ml.A << endl;
    01350         os << "mu:" << ml.mu_const << endl;
    01351         os << "R:" << ml._R() << endl;
    01352         return os;
    01353 };
    01354
    01355 }
    01356 #endif //EF_H