[176] | 1 | /*! |
---|
| 2 | \file |
---|
| 3 | \brief Mergers for combination of pdfs |
---|
| 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 | |
---|
[384] | 13 | #ifndef MERGER_H |
---|
| 14 | #define MERGER_H |
---|
[176] | 15 | |
---|
[262] | 16 | |
---|
[384] | 17 | #include "../estim/mixtures.h" |
---|
[556] | 18 | #include "discrete.h" |
---|
[176] | 19 | |
---|
[477] | 20 | namespace bdm { |
---|
[384] | 21 | using std::string; |
---|
[176] | 22 | |
---|
[384] | 23 | //!Merging methods |
---|
| 24 | enum MERGER_METHOD {ARITHMETIC = 1, GEOMETRIC = 2, LOGNORMAL = 3}; |
---|
[176] | 25 | |
---|
[1085] | 26 | /*! \brief Abstract common class for mergers |
---|
[1072] | 27 | |
---|
[1095] | 28 | Merging is defined as a combination of information from \c source pdfs into a single pdf called the \c merger. |
---|
| 29 | |
---|
| 30 | \f[ |
---|
| 31 | f_{merged}(x) \leftarrow {f_1(x),f_2(x),f_3(x)} |
---|
| 32 | \f] |
---|
| 33 | |
---|
| 34 | In the basic form, the random variables of all source pdfs must be identical, e.g. \f$ x \f$ from the example above. |
---|
| 35 | |
---|
| 36 | Extension of the merger to a more demanding scenarios with fragmental sources is available in offspring bdm::merger_mix. |
---|
| 37 | |
---|
| 38 | This abstract class defines only the interface to using merger: operation merge(), and operation merger() which returns the result. |
---|
[1072] | 39 | */ |
---|
| 40 | class MergerBase: public root{ |
---|
| 41 | public: |
---|
| 42 | //! weights of the sources -- no restrictions |
---|
| 43 | vec w; |
---|
[1079] | 44 | //! source pdfs -- conditioning allowed |
---|
[1072] | 45 | Array<shared_ptr<epdf> > sources; |
---|
| 46 | //! merge |
---|
[1079] | 47 | virtual void merge() NOT_IMPLEMENTED_VOID; |
---|
[1072] | 48 | //! |
---|
| 49 | //! check if all epdfs are on the same support |
---|
[1083] | 50 | virtual void set_sources(const Array<shared_ptr<pdf> > &A) { |
---|
| 51 | int n=A.length(); |
---|
[1072] | 52 | |
---|
| 53 | bdm_assert(n>0,"merger has no sources to merge"); |
---|
| 54 | // check if weights are ok |
---|
| 55 | if (w.length()!=n) { |
---|
| 56 | w = ones(n)/n; |
---|
| 57 | } |
---|
| 58 | |
---|
| 59 | // check compatibility of sources -- no types are needed |
---|
[1083] | 60 | int dim0 = A(0)->dimension(); |
---|
[1072] | 61 | for (int i=0; i<n; i++){ |
---|
[1083] | 62 | const epdf *si=dynamic_cast<const epdf*>(A(i).get()); |
---|
| 63 | if (si){ |
---|
[1114] | 64 | sources(i) = shared_ptr<epdf>(const_cast<epdf*>(si)); |
---|
[1083] | 65 | } |
---|
[1072] | 66 | bdm_assert(sources(i)->dimension()==dim0, "Merger: Incompatible dimensions of sources"); |
---|
| 67 | if (sources(i)->isnamed()){ |
---|
| 68 | //check if rv match |
---|
| 69 | } |
---|
| 70 | } |
---|
| 71 | } |
---|
| 72 | //! return the result of merging |
---|
| 73 | virtual const epdf& merger()=0; |
---|
| 74 | }; |
---|
| 75 | |
---|
[1085] | 76 | //! \brief Merger into normal density, max entropy approximator for 2 moments (mean+variance) |
---|
[1072] | 77 | template<class sq_T> |
---|
| 78 | class ENormMerger: public MergerBase{ |
---|
| 79 | protected: |
---|
| 80 | //!internal epdf |
---|
| 81 | enorm<sq_T> iepdf; |
---|
| 82 | public: |
---|
| 83 | ENormMerger():method(GEOMETRIC){}; |
---|
| 84 | MERGER_METHOD method; |
---|
| 85 | void merge(){ |
---|
| 86 | int n = sources.length(); |
---|
| 87 | int dim = sources(0)->dimension(); |
---|
| 88 | sq_T Var(zeros(dim)); |
---|
| 89 | vec mea=zeros(dim); |
---|
| 90 | // go through all sources |
---|
[1079] | 91 | epdf * si; |
---|
[1072] | 92 | for (int i=0; i<n; i++){ |
---|
[1079] | 93 | si = dynamic_cast<epdf*>(sources(i).get()); // transform pdf into epdf |
---|
| 94 | |
---|
| 95 | if (!si) {bdm_error("Can't merger this type of pdf: " + sources(i)->to_string());} |
---|
| 96 | |
---|
| 97 | sq_T Ci = si->covariance(); |
---|
[1072] | 98 | sq_T wCi = Ci; wCi*=w(i); |
---|
[1079] | 99 | vec mi = si->mean(); |
---|
[1072] | 100 | switch (method) { |
---|
| 101 | case ARITHMETIC:{ |
---|
| 102 | Var += wCi; |
---|
| 103 | Var.opupdt(mi,w(i)); // + mean * mean' |
---|
| 104 | mea += w(i)*mi; |
---|
| 105 | break; |
---|
| 106 | } |
---|
| 107 | case GEOMETRIC:{ |
---|
| 108 | Var += wCi; |
---|
| 109 | sq_T iCi; |
---|
| 110 | Ci.inv(iCi); |
---|
| 111 | mea += iCi.to_mat()*w(i)*mi; |
---|
| 112 | break; |
---|
| 113 | } |
---|
| 114 | default: |
---|
| 115 | bdm_error("Method not implemneted"); |
---|
| 116 | } |
---|
| 117 | } |
---|
| 118 | // post pocessing |
---|
| 119 | switch (method) { |
---|
| 120 | case ARITHMETIC: |
---|
| 121 | iepdf._R()=Var; |
---|
| 122 | iepdf._R().opupdt(mea,-1.0); // Var - mean * mean' |
---|
| 123 | iepdf._mu()=mea; |
---|
| 124 | break; |
---|
| 125 | case GEOMETRIC: |
---|
| 126 | iepdf._R()=Var; |
---|
| 127 | iepdf._mu() = Var.to_mat()*mea; // mean=sqrt(Var)*mea |
---|
| 128 | break; |
---|
| 129 | default: |
---|
| 130 | bdm_error("Method not implemneted"); |
---|
| 131 | } |
---|
| 132 | } |
---|
| 133 | //! access function |
---|
| 134 | const epdf& merger(){return iepdf;} |
---|
| 135 | }; |
---|
| 136 | |
---|
[384] | 137 | /*! |
---|
[1068] | 138 | @brief Base class (interface) for for general combination of pdfs on discrete support |
---|
[176] | 139 | |
---|
[384] | 140 | Mixtures of Gaussian densities are used internally. Switching to other densities should be trivial. |
---|
[198] | 141 | |
---|
[384] | 142 | The merged pdfs are expected to be of the form: |
---|
| 143 | \f[ f(x_i|y_i), i=1..n \f] |
---|
| 144 | where the resulting merger is a density on \f$ \cup [x_i,y_i] \f$ . |
---|
| 145 | Note that all variables will be joined. |
---|
[205] | 146 | |
---|
[384] | 147 | As a result of this feature, each source must be extended to common support |
---|
| 148 | \f[ f(z_i|y_i,x_i) f(x_i|y_i) f(y_i) i=1..n \f] |
---|
| 149 | where \f$ z_i \f$ accumulate variables that were not in the original source. |
---|
| 150 | These extensions are calculated on-the-fly. |
---|
[299] | 151 | |
---|
[384] | 152 | However, these operations can not be performed in general. Hence, this class merges only sources on common support, \f$ y_i={}, z_i={}, \forall i \f$. |
---|
| 153 | For merging of more general cases, use offsprings merger_mix and merger_grid. |
---|
[1079] | 154 | |
---|
| 155 | \TODO use sources for extensions |
---|
[384] | 156 | */ |
---|
[1079] | 157 | class MergerDiscrete : public MergerBase { |
---|
[477] | 158 | protected: |
---|
[1079] | 159 | //! partial sources |
---|
| 160 | Array<shared_ptr< pdf > > part_sources; |
---|
| 161 | |
---|
[1064] | 162 | //! Data link for each pdf in pdfs |
---|
| 163 | Array<datalink_m2e*> dls; |
---|
[507] | 164 | |
---|
[1064] | 165 | //! Array of rvs that are not modelled by pdfs at all, \f$ z_i \f$ |
---|
| 166 | Array<RV> rvzs; |
---|
[507] | 167 | |
---|
[1064] | 168 | //! Data Links for extension \f$ f(z_i|x_i,y_i) \f$ |
---|
| 169 | Array<datalink_m2e*> zdls; |
---|
[507] | 170 | |
---|
[1064] | 171 | //! number of support points |
---|
| 172 | int Npoints; |
---|
[507] | 173 | |
---|
[1064] | 174 | //! number of sources |
---|
| 175 | int Nsources; |
---|
[384] | 176 | |
---|
[1064] | 177 | //! switch of the methoh used for merging |
---|
| 178 | MERGER_METHOD METHOD; |
---|
| 179 | //! Default for METHOD |
---|
| 180 | static const MERGER_METHOD DFLT_METHOD; |
---|
[384] | 181 | |
---|
[1064] | 182 | //!Prior on the log-normal merging model |
---|
| 183 | double beta; |
---|
| 184 | //! default for beta |
---|
| 185 | static const double DFLT_beta; |
---|
[384] | 186 | |
---|
[1064] | 187 | //! Projection to empirical density (could also be piece-wise linear) |
---|
| 188 | eEmp eSmp; |
---|
[384] | 189 | |
---|
[1064] | 190 | //! debug or not debug |
---|
| 191 | bool DBG; |
---|
[423] | 192 | |
---|
[1064] | 193 | //! debugging file |
---|
| 194 | it_file* dbg_file; |
---|
[477] | 195 | public: |
---|
[1064] | 196 | //! \name Constructors |
---|
| 197 | //! @{ |
---|
[423] | 198 | |
---|
[1064] | 199 | //! Default constructor |
---|
[1079] | 200 | MergerDiscrete () : Npoints ( 0 ), Nsources ( 0 ), DBG ( false ), dbg_file ( 0 ) { |
---|
[1064] | 201 | } |
---|
[384] | 202 | |
---|
[1064] | 203 | //!Constructor from sources |
---|
[1079] | 204 | MergerDiscrete ( const Array<shared_ptr<pdf> > &S ); |
---|
[384] | 205 | |
---|
[1064] | 206 | //! Function setting the main internal structures |
---|
[1079] | 207 | void set_sources ( const Array<shared_ptr<pdf> > &PartSources ); |
---|
[384] | 208 | |
---|
[1064] | 209 | //! Set support points from rectangular grid |
---|
| 210 | void set_support ( rectangular_support &Sup ); |
---|
[739] | 211 | |
---|
[1064] | 212 | //! Set support points from dicrete grid |
---|
| 213 | void set_support ( discrete_support &Sup ) { |
---|
| 214 | Npoints = Sup.points(); |
---|
| 215 | eSmp.set_parameters ( Sup._Spoints() ); |
---|
| 216 | eSmp.validate(); |
---|
| 217 | } |
---|
[957] | 218 | |
---|
| 219 | |
---|
[1064] | 220 | //! set debug file |
---|
| 221 | void set_debug_file ( const string fname ) { |
---|
| 222 | if ( DBG ) delete dbg_file; |
---|
| 223 | dbg_file = new it_file ( fname ); |
---|
| 224 | DBG = ( dbg_file != 0 ); |
---|
| 225 | } |
---|
[957] | 226 | |
---|
[1064] | 227 | //! Set internal parameters used in approximation |
---|
[1079] | 228 | void set_method ( MERGER_METHOD MTH = GEOMETRIC, double beta0 = 1.2 ) { |
---|
[1064] | 229 | METHOD = MTH; |
---|
| 230 | beta = beta0; |
---|
| 231 | } |
---|
| 232 | //! Set support points from a pdf by drawing N samples |
---|
| 233 | void set_support ( const epdf &overall, int N ) { |
---|
| 234 | eSmp.set_statistics ( overall, N ); |
---|
| 235 | Npoints = N; |
---|
| 236 | } |
---|
[477] | 237 | |
---|
[1064] | 238 | //! Destructor |
---|
[1079] | 239 | virtual ~MergerDiscrete() { |
---|
[1064] | 240 | for ( int i = 0; i < Nsources; i++ ) { |
---|
| 241 | delete dls ( i ); |
---|
| 242 | delete zdls ( i ); |
---|
| 243 | } |
---|
| 244 | if ( DBG ) delete dbg_file; |
---|
| 245 | }; |
---|
| 246 | //!@} |
---|
[388] | 247 | |
---|
[1064] | 248 | //! \name Mathematical operations |
---|
| 249 | //!@{ |
---|
[388] | 250 | |
---|
[1064] | 251 | //!Merge given sources in given points |
---|
[1079] | 252 | void merge (); |
---|
[388] | 253 | |
---|
[1064] | 254 | //! Merge log-likelihood values in points using method specified by parameter METHOD |
---|
| 255 | vec merge_points ( mat &lW ); |
---|
[388] | 256 | |
---|
[1079] | 257 | //!@} |
---|
[477] | 258 | |
---|
[1064] | 259 | //! \name Access to attributes |
---|
| 260 | //! @{ |
---|
[384] | 261 | |
---|
[1064] | 262 | //! Access function |
---|
| 263 | eEmp& _Smp() { |
---|
| 264 | return eSmp; |
---|
| 265 | } |
---|
[388] | 266 | |
---|
[1079] | 267 | eEmp & merger() {return eSmp;} |
---|
| 268 | |
---|
[1068] | 269 | /*! Create object from the following structure |
---|
| 270 | |
---|
| 271 | \code |
---|
[1086] | 272 | class = 'MergerDiscrete'; |
---|
[1068] | 273 | method = 'arithmetic' or 'geometric' or 'lognormal'; % name of the model used for merging |
---|
| 274 | --- compulsory only for lognormal merging model --- |
---|
| 275 | beta = x; % scalar prior on the log-normal merging model |
---|
| 276 | --- optional fields --- |
---|
| 277 | dbg_file = '...'; % name of the file to store debug informations |
---|
| 278 | --- inherited fields --- |
---|
[1086] | 279 | bdm::MergerBase::from_setting |
---|
[1068] | 280 | \endcode |
---|
| 281 | */ |
---|
[1064] | 282 | void from_setting ( const Setting& set ); |
---|
[388] | 283 | |
---|
[1064] | 284 | void to_setting (Setting &set) const ; |
---|
| 285 | |
---|
| 286 | void validate() ; |
---|
| 287 | //!@} |
---|
[384] | 288 | }; |
---|
[1079] | 289 | UIREGISTER ( MergerDiscrete ); |
---|
| 290 | SHAREDPTR ( MergerDiscrete ); |
---|
[384] | 291 | |
---|
[1068] | 292 | //! \brief Merger using importance sampling with mixture proposal density |
---|
[1079] | 293 | class merger_mix : public MergerDiscrete { |
---|
[477] | 294 | protected: |
---|
[1064] | 295 | //!Internal mixture of EF models |
---|
| 296 | MixEF Mix; |
---|
| 297 | //!Number of components in a mixture |
---|
| 298 | int Ncoms; |
---|
| 299 | //! coefficient of resampling [0,1] |
---|
| 300 | double effss_coef; |
---|
| 301 | //! stop after niter iterations |
---|
| 302 | int stop_niter; |
---|
[384] | 303 | |
---|
[1064] | 304 | //! default value for Ncoms |
---|
| 305 | static const int DFLT_Ncoms; |
---|
| 306 | //! default value for efss_coef; |
---|
| 307 | static const double DFLT_effss_coef; |
---|
[388] | 308 | |
---|
[477] | 309 | public: |
---|
[1064] | 310 | //!\name Constructors |
---|
| 311 | //!@{ |
---|
| 312 | merger_mix () : Ncoms ( 0 ), effss_coef ( 0 ), stop_niter ( 0 ) { } |
---|
[507] | 313 | |
---|
[1064] | 314 | merger_mix ( const Array<shared_ptr<pdf> > &S ) : |
---|
| 315 | Ncoms ( 0 ), effss_coef ( 0 ), stop_niter ( 0 ) { |
---|
| 316 | set_sources ( S ); |
---|
| 317 | } |
---|
[507] | 318 | |
---|
[1064] | 319 | //! Set sources and prepare all internal structures |
---|
| 320 | void set_sources ( const Array<shared_ptr<pdf> > &S ) { |
---|
[1079] | 321 | MergerDiscrete::set_sources ( S ); |
---|
[1064] | 322 | //Nsources = S.length(); |
---|
| 323 | } |
---|
[507] | 324 | |
---|
[1064] | 325 | //! Set internal parameters used in approximation |
---|
| 326 | void set_parameters ( int Ncoms0 = DFLT_Ncoms, double effss_coef0 = DFLT_effss_coef ) { |
---|
| 327 | Ncoms = Ncoms0; |
---|
| 328 | effss_coef = effss_coef0; |
---|
| 329 | } |
---|
| 330 | //!@} |
---|
[388] | 331 | |
---|
[1064] | 332 | //! \name Mathematical operations |
---|
| 333 | //!@{ |
---|
[384] | 334 | |
---|
[1064] | 335 | //!Merge values using mixture approximation |
---|
| 336 | void merge (); |
---|
[384] | 337 | |
---|
[1064] | 338 | //! sample from the approximating mixture |
---|
| 339 | vec sample () const { |
---|
| 340 | return Mix.posterior().sample(); |
---|
| 341 | } |
---|
| 342 | //! loglikelihood computed on mixture models |
---|
| 343 | double evallog ( const vec &yt ) const { |
---|
| 344 | vec dtf = ones ( yt.length() + 1 ); |
---|
| 345 | dtf.set_subvector ( 0, yt ); |
---|
| 346 | return Mix.logpred ( dtf, vec(0) ); |
---|
| 347 | } |
---|
| 348 | //!@} |
---|
[477] | 349 | |
---|
[1064] | 350 | //!\name Access functions |
---|
| 351 | //!@{ |
---|
[192] | 352 | //! Access function |
---|
[1064] | 353 | MixEF& _Mix() { |
---|
| 354 | return Mix; |
---|
| 355 | } |
---|
| 356 | //! Access function |
---|
| 357 | emix* proposal() { |
---|
| 358 | emix* tmp = Mix.epredictor(); |
---|
[1079] | 359 | tmp->set_rv ( merger()._rv() ); |
---|
[1064] | 360 | return tmp; |
---|
| 361 | } |
---|
[1068] | 362 | /*! Create object from the following structure |
---|
| 363 | |
---|
| 364 | \code |
---|
| 365 | class = 'merger_mix'; |
---|
| 366 | --- optional fields --- |
---|
| 367 | ncoms = []; % number of components in a mixture |
---|
| 368 | effss_coef = []; % coefficient of resampling from interval [0,1] |
---|
| 369 | stop_niter = []; % number of interations |
---|
| 370 | --- inherited fields --- |
---|
| 371 | bdm::merger_base::from_setting |
---|
| 372 | \endcode |
---|
| 373 | If the optional fields are not given, they will be filled as follows: |
---|
| 374 | \code |
---|
| 375 | ncoms = 10; |
---|
| 376 | effss_coef = 0.9; |
---|
| 377 | stop_niter = 10; |
---|
| 378 | \endcode |
---|
| 379 | |
---|
| 380 | */ |
---|
[1064] | 381 | void from_setting ( const Setting& set ); |
---|
[1068] | 382 | |
---|
| 383 | void to_setting (Setting &set) const; |
---|
[1064] | 384 | void validate(); |
---|
[388] | 385 | |
---|
[1079] | 386 | // const emix &merger(){return *Mpred;} |
---|
[1064] | 387 | //! @} |
---|
[384] | 388 | }; |
---|
[477] | 389 | UIREGISTER ( merger_mix ); |
---|
[529] | 390 | SHAREDPTR ( merger_mix ); |
---|
[176] | 391 | |
---|
[254] | 392 | } |
---|
[176] | 393 | |
---|
| 394 | #endif // MER_H |
---|