[176] | 1 | /*! |
---|
| 2 | \file |
---|
| 3 | \brief Bayesian Filtering for mixtures of exponential family (EF) members |
---|
| 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 | |
---|
[394] | 13 | #ifndef MIXTURES_H |
---|
| 14 | #define MIXTURES_H |
---|
[176] | 15 | |
---|
[262] | 16 | |
---|
[384] | 17 | #include "../math/functions.h" |
---|
| 18 | #include "../stat/exp_family.h" |
---|
[394] | 19 | #include "../stat/emix.h" |
---|
[997] | 20 | #include "arx.h" |
---|
[176] | 21 | |
---|
[286] | 22 | namespace bdm { |
---|
[176] | 23 | |
---|
[536] | 24 | //! enum switch for internal approximation used in MixEF |
---|
[189] | 25 | enum MixEF_METHOD { EM = 0, QB = 1}; |
---|
| 26 | |
---|
[176] | 27 | /*! |
---|
| 28 | * \brief Mixture of Exponential Family Densities |
---|
| 29 | |
---|
| 30 | An approximate estimation method for models with latent discrete variable, such as |
---|
| 31 | mixture models of the following kind: |
---|
| 32 | \f[ |
---|
[1077] | 33 | f(y_t|\psi_t, heta) = \sum_{i=1}^{n} w_i f(y_t|\psi_t, heta_i) |
---|
[176] | 34 | \f] |
---|
[1077] | 35 | where \f$\psi\f$ is a known function of past outputs, \f$w=[w_1,\ldots,w_n]\f$ are component weights, and component parameters \f$ heta_i\f$ are assumed to be mutually independent. \f$ heta\f$ is an aggregation af all component parameters and weights, i.e. \f$ heta = [ heta_1,\ldots, heta_n,w]\f$. |
---|
[176] | 36 | |
---|
| 37 | The characteristic feature of this model is that if the exact values of the latent variable were known, estimation of the parameters can be handled by a single model. For example, for the case of mixture models, posterior density for each component parameters would be a BayesianModel from Exponential Family. |
---|
| 38 | |
---|
[180] | 39 | This class uses EM-style type algorithms for estimation of its parameters. Under this simplification, the posterior density is a product of exponential family members, hence under EM-style approximate estimation this class itself belongs to the exponential family. |
---|
[176] | 40 | |
---|
[1014] | 41 | Two methods are provided: |
---|
| 42 | - QB where the probability of being from one component is computed using predictors, |
---|
| 43 | - EM where the data is assigned to component with the highest likelihood (winner takes all) |
---|
| 44 | |
---|
| 45 | These methods are stored in attribute options: |
---|
| 46 | - method: ["EM","QB"] estimation method as mentioned above, QB is default |
---|
| 47 | - max_niter: maximum of iterations in bayes_batch() |
---|
| 48 | |
---|
[176] | 49 | */ |
---|
[197] | 50 | class MixEF: public BMEF { |
---|
[176] | 51 | protected: |
---|
[1077] | 52 | //! Models for Components of \f$ heta_i\f$ |
---|
[1064] | 53 | Array<BMEF*> Coms; |
---|
| 54 | //! Statistics for weights |
---|
| 55 | multiBM weights; |
---|
| 56 | //aux |
---|
| 57 | friend class eprod_mix; |
---|
[1066] | 58 | |
---|
| 59 | //! \brief Posterior on component parameters |
---|
[1064] | 60 | class eprod_mix: public eprod_base { |
---|
| 61 | protected: |
---|
| 62 | const MixEF &mix; // pointer to parent.n |
---|
| 63 | public: |
---|
| 64 | eprod_mix(const MixEF &m):mix(m) {} |
---|
| 65 | const epdf* factor(int i) const { |
---|
| 66 | return (i==(mix.Coms.length()-1)) ? &mix.weights.posterior() : &mix.Coms(i)->posterior(); |
---|
| 67 | } |
---|
| 68 | const int no_factors()const { |
---|
| 69 | return mix.Coms.length()+1; |
---|
| 70 | } |
---|
| 71 | } est; |
---|
| 72 | ////!indices of component rvc in common rvc |
---|
[286] | 73 | |
---|
[1077] | 74 | class Options: public root { |
---|
[1064] | 75 | public: |
---|
| 76 | //! Flag for a method that is used in the inference |
---|
| 77 | MixEF_METHOD method; |
---|
[286] | 78 | |
---|
[1064] | 79 | //! maximum number of iterations |
---|
| 80 | int max_niter; |
---|
| 81 | |
---|
[1077] | 82 | Options():method(QB),max_niter(10) {}; |
---|
[1064] | 83 | |
---|
[1077] | 84 | /*! Create object from the following structure |
---|
| 85 | |
---|
[1064] | 86 | \code |
---|
[1077] | 87 | --- optional fields --- |
---|
| 88 | method = '...'; % 'EM' or 'QB', methods of computing bayes |
---|
| 89 | max_iter = []; % maximum number of iterations in bayes_batch |
---|
| 90 | --- inherited fields --- |
---|
| 91 | bdm::root::from_setting |
---|
[1064] | 92 | \endcode |
---|
[1077] | 93 | If the optional fields are not given, they will be filled as follows: |
---|
| 94 | \code |
---|
| 95 | max_niter = 10; |
---|
| 96 | method = 'QB'; |
---|
| 97 | \endcode |
---|
[1064] | 98 | */ |
---|
| 99 | void from_setting(const Setting &set) { |
---|
| 100 | string meth; |
---|
| 101 | UI::get(meth,set,"method",UI::optional); |
---|
| 102 | if (meth=="EM") { |
---|
| 103 | method=EM; |
---|
| 104 | } |
---|
| 105 | max_niter =10; |
---|
| 106 | UI::get(max_niter,set,"max_niter",UI::optional); |
---|
| 107 | }; |
---|
[1077] | 108 | |
---|
[1064] | 109 | void to_setting(Setting &set)const { |
---|
| 110 | string meth=(method==EM ? "EM" : "QB"); |
---|
| 111 | UI::save(meth,set,"method"); |
---|
| 112 | UI::save(max_niter,set,"max_niter"); |
---|
| 113 | }; |
---|
| 114 | }; |
---|
| 115 | |
---|
[1077] | 116 | Options options; |
---|
[176] | 117 | public: |
---|
[1064] | 118 | //! Full constructor |
---|
| 119 | MixEF ( const Array<BMEF*> &Coms0, const vec &alpha0 ) : |
---|
| 120 | BMEF ( ), Coms ( Coms0.length() ), |
---|
| 121 | weights (), est(*this), options() { |
---|
| 122 | for ( int i = 0; i < Coms0.length(); i++ ) { |
---|
| 123 | Coms ( i ) = ( BMEF* ) Coms0 ( i )->_copy(); |
---|
| 124 | } |
---|
| 125 | weights.set_parameters(alpha0); |
---|
| 126 | weights.validate(); |
---|
| 127 | } |
---|
[565] | 128 | |
---|
[1064] | 129 | //! Constructor of empty mixture |
---|
| 130 | MixEF () : |
---|
| 131 | BMEF ( ), Coms ( 0 ), |
---|
| 132 | weights (), est(*this), options() { |
---|
| 133 | } |
---|
| 134 | //! Copy constructor |
---|
| 135 | MixEF ( const MixEF &M2 ) : BMEF ( ), Coms ( M2.Coms.length() ), |
---|
| 136 | weights ( M2.weights ), est(*this), options(M2.options) { |
---|
| 137 | for ( int i = 0; i < M2.Coms.length(); i++ ) { |
---|
| 138 | Coms ( i ) = (BMEF*) M2.Coms ( i )->_copy(); |
---|
| 139 | } |
---|
| 140 | } |
---|
[565] | 141 | |
---|
[1064] | 142 | //! Initializing the mixture by a random pick of centroids from data |
---|
| 143 | //! \param Com0 Initial component - necessary to determine its type. |
---|
| 144 | //! \param Data Data on which the initialization will be done |
---|
| 145 | //! \param c Initial number of components, default=5 |
---|
| 146 | //! when the number of data records (columns of Data) is equal to the number of requested components, each data is used, otherwise, they are picked randomly. |
---|
| 147 | void init ( BMEF* Com0, const mat &Data, const int c = 5 ); |
---|
| 148 | //Destructor |
---|
| 149 | //! Recursive EM-like algorithm (QB-variant), see Karny et. al, 2006 |
---|
| 150 | void bayes ( const vec &yt, const vec &cond ); |
---|
| 151 | //! batch weighted Bayes rule |
---|
| 152 | double bayes_batch_weighted ( const mat &yt, const mat &cond, const vec &wData ); |
---|
| 153 | double bayes_batch ( const mat &yt, const mat &cond) { |
---|
| 154 | return bayes_batch_weighted(yt,cond,ones(yt.cols())); |
---|
| 155 | }; |
---|
| 156 | double logpred ( const vec &yt, const vec &cond ) const; |
---|
| 157 | //! return correctly typed posterior (covariant return) |
---|
| 158 | const eprod_mix& posterior() const { |
---|
| 159 | return est; |
---|
| 160 | } |
---|
[746] | 161 | |
---|
[1064] | 162 | emix* epredictor(const vec &cond=vec()) const; |
---|
| 163 | //! Flatten the density as if it was not estimated from the data |
---|
| 164 | void flatten ( const BMEF* M2, double weight ); |
---|
| 165 | //! Access function |
---|
| 166 | BMEF* _Coms ( int i ) { |
---|
| 167 | return Coms ( i ); |
---|
| 168 | } |
---|
[286] | 169 | |
---|
[1064] | 170 | //!Set which method is to be used |
---|
| 171 | void set_method ( MixEF_METHOD M ) { |
---|
| 172 | options.method = M; |
---|
| 173 | } |
---|
[737] | 174 | |
---|
[1064] | 175 | void to_setting ( Setting &set ) const { |
---|
| 176 | BMEF::to_setting( set ); |
---|
| 177 | UI::save ( Coms, set, "Coms" ); |
---|
| 178 | UI::save ( &weights, set, "weights" ); |
---|
| 179 | UI::save (options, set, "options"); |
---|
| 180 | } |
---|
[1077] | 181 | |
---|
| 182 | /*! Create object from the following structure |
---|
| 183 | |
---|
[1064] | 184 | \code |
---|
[1077] | 185 | class = 'MixEF'; |
---|
| 186 | Coms = { list of bdm::BMEF }; % list of components containing any offsprings of Bayesian models BMEF, bdm::BMEF::from_setting |
---|
| 187 | weights = [...]; % vector containing weights of components |
---|
| 188 | --- optional fields --- |
---|
| 189 | options = configuration of bdm::MixEF::Options; % see MixEF::Options, bdm::MixEF::Options::from_setting |
---|
| 190 | --- inherited fields --- |
---|
| 191 | bdm::BMEF::from_setting |
---|
[1064] | 192 | \endcode |
---|
| 193 | */ |
---|
| 194 | void from_setting (const Setting &set ) { |
---|
| 195 | BMEF::from_setting( set ); |
---|
| 196 | UI::get ( Coms, set, "Coms" ); |
---|
| 197 | UI::get ( weights, set, "weights" ); |
---|
| 198 | UI::get (options, set, "options",UI::optional); |
---|
| 199 | } |
---|
[176] | 200 | }; |
---|
[746] | 201 | UIREGISTER ( MixEF ); |
---|
[176] | 202 | |
---|
[1085] | 203 | //! \brief Product of ARX models forming an ARX model with potentially reduced structure |
---|
[997] | 204 | class ARXprod: public ProdBMBase { |
---|
[1064] | 205 | Array<shared_ptr<ARX> > arxs; |
---|
| 206 | public: |
---|
| 207 | ARX* bm(int i) { |
---|
| 208 | return arxs(i).get(); |
---|
| 209 | } |
---|
| 210 | int no_bms() { |
---|
| 211 | return arxs.length(); |
---|
| 212 | } |
---|
[997] | 213 | }; |
---|
| 214 | UIREGISTER(ARXprod); |
---|
| 215 | |
---|
[254] | 216 | } |
---|
[384] | 217 | #endif // MIXTURES_H |
---|