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