1 /* _______________________________________________________________________
2
3 PECOS: Parallel Environment for Creation Of Stochastics
4 Copyright (c) 2011, Sandia National Laboratories.
5 This software is distributed under the GNU Lesser General Public License.
6 For more information, see the README file in the top Pecos directory.
7 _______________________________________________________________________ */
8
9 //- Class: IntervalRandomVariable
10 //- Description: A random variable described by discrete values + probabilities
11 //- Owner: Mike Eldred
12 //- Revised by:
13 //- Version:
14
15 #ifndef INTERVAL_RANDOM_VARIABLE_HPP
16 #define INTERVAL_RANDOM_VARIABLE_HPP
17
18 #include "RandomVariable.hpp"
19 #include "HistogramBinRandomVariable.hpp"
20 #include "DiscreteSetRandomVariable.hpp"
21 #include "pecos_stat_util.hpp"
22
23 namespace Pecos {
24
25
26 /// Derived random variable class for interval random variables
27
28 /** Manages basic probability assignments (BPAs: pairings of
29 intervals with probabilities) for types int and real. */
30
31 template <typename T>
32 class IntervalRandomVariable: public RandomVariable
33 {
34 public:
35
36 //
37 //- Heading: Constructors and destructor
38 //
39
40 /// default constructor
41 IntervalRandomVariable();
42 /// alternate constructor
43 IntervalRandomVariable(const std::map<std::pair<T, T>, Real>& bpa);
44 /// destructor
45 ~IntervalRandomVariable();
46
47 //
48 //- Heading: Virtual function redefinitions
49 //
50
51 Real cdf(Real x) const;
52 Real ccdf(Real x) const;
53 Real inverse_cdf(Real p_cdf) const;
54 Real inverse_ccdf(Real p_ccdf) const;
55
56 Real pdf(Real x) const;
57 Real pdf_gradient(Real x) const;
58 Real pdf_hessian(Real x) const;
59
60 Real mean() const;
61 //Real median() const;
62 Real mode() const;
63 Real standard_deviation() const;
64 Real variance() const;
65
66 RealRealPair moments() const;
67 RealRealPair distribution_bounds() const;
68
69 //Real coefficient_of_variation() const;
70
71 void pull_parameter(short dist_param,
72 std::map<std::pair<T, T>, Real>& bpa) const;
73 void push_parameter(short dist_param,
74 const std::map<std::pair<T, T>, Real>& bpa);
75
76 void copy_parameters(const RandomVariable& rv);
77
78 //
79 //- Heading: Member functions
80 //
81
82 /// update intervalBPA
83 void update(const std::map<std::pair<T, T>, Real>& bpa);
84
85 /// verify that valueProbPairs is defined
86 void check_vpp() const;
87 /// activate valueProbPairs tracking by copying intervalBPA
88 void activate_vpp();
89 /// update valueProbPairs if currently in use
90 void update_vpp_if_active();
91
92 //
93 //- Heading: Static member functions (global utilities)
94 //
95
96 static void moments_from_params(const std::map<std::pair<T, T>, Real>& bpa,
97 Real& mean, Real& std_dev);
98
99 //static Real pdf(Real x, const RealVector& dist_params);
100
101 protected:
102
103 //
104 //- Heading: Data
105 //
106
107 /// interval-probability pairs for basic probability assignment (BPA)
108 std::map<std::pair<T, T>, Real> intervalBPA;
109
110 /// collapse (overlapping, disjoint) intervals into a histogram-like format
111 /// (defined if needed for use in moments/PDF/CDF)
112 std::map<T, Real> valueProbPairs;
113 };
114
115
116 //// GENERIC ////
117
118
119 template <typename T>
IntervalRandomVariable()120 IntervalRandomVariable<T>::IntervalRandomVariable():
121 RandomVariable(BaseConstructor())
122 { }
123
124
125 template <typename T>
126 IntervalRandomVariable<T>::
IntervalRandomVariable(const std::map<std::pair<T,T>,Real> & bpa)127 IntervalRandomVariable(const std::map<std::pair<T, T>, Real>& bpa):
128 RandomVariable(BaseConstructor()), intervalBPA(bpa)
129 { }
130
131
132 template <typename T>
~IntervalRandomVariable()133 IntervalRandomVariable<T>::~IntervalRandomVariable()
134 { }
135
136
137 template <typename T>
check_vpp() const138 void IntervalRandomVariable<T>::check_vpp() const
139 {
140 if (valueProbPairs.empty()) {
141 PCerr << "Error: valueProbPairs not activated in IntervalRandomVariable<T>."
142 << std::endl;
143 abort_handler(-1);
144 }
145 }
146
147
148 template <typename T>
activate_vpp()149 void IntervalRandomVariable<T>::activate_vpp()
150 {
151 if (valueProbPairs.empty()) // activate once to avoid copying repeatedly
152 intervals_to_xy_pdf(intervalBPA, valueProbPairs);
153 }
154
155
156 template <typename T>
update_vpp_if_active()157 void IntervalRandomVariable<T>::update_vpp_if_active()
158 {
159 if (!valueProbPairs.empty()) // update valueProbPairs if already activated
160 intervals_to_xy_pdf(intervalBPA, valueProbPairs);
161 }
162
163
164 template <typename T>
165 void IntervalRandomVariable<T>::
update(const std::map<std::pair<T,T>,Real> & bpa)166 update(const std::map<std::pair<T, T>, Real>& bpa)
167 {
168 intervalBPA = bpa;
169 update_vpp_if_active(); // update valueProbPairs if in use, else leave empty
170 }
171 // specializations could be used for assigning ranVarType, or could employ
172 // std::is_same for type identification. Simplest: ranVarType assigned at
173 // bottom of RandomVariable::get_random_variable().
174
175
176 template <typename T>
177 void IntervalRandomVariable<T>::
pull_parameter(short dist_param,std::map<std::pair<T,T>,Real> & bpa) const178 pull_parameter(short dist_param, std::map<std::pair<T, T>, Real>& bpa) const
179 {
180 // could specialize template, but case aggregation seems adequate
181
182 switch (dist_param) {
183 case CIU_BPA: case DIU_BPA: bpa = intervalBPA; break;
184 default:
185 PCerr << "Error: update failure for distribution parameter " << dist_param
186 << " in IntervalRandomVariable::pull_parameter(T)." << std::endl;
187 abort_handler(-1); break;
188 }
189 }
190
191
192 template <typename T>
193 void IntervalRandomVariable<T>::
push_parameter(short dist_param,const std::map<std::pair<T,T>,Real> & bpa)194 push_parameter(short dist_param, const std::map<std::pair<T, T>, Real>& bpa)
195 {
196 // could specialize template, but case aggregation seems adequate
197
198 switch (dist_param) {
199 case CIU_BPA: case DIU_BPA:
200 intervalBPA = bpa; break;
201 default:
202 PCerr << "Error: update failure for distribution parameter " << dist_param
203 << " in IntervalRandomVariable::push_parameter(T)." << std::endl;
204 abort_handler(-1); break;
205 }
206 update_vpp_if_active(); // update valueProbPairs if in use, else leave empty
207 }
208
209
210 template <typename T>
copy_parameters(const RandomVariable & rv)211 void IntervalRandomVariable<T>::copy_parameters(const RandomVariable& rv)
212 {
213 switch (ranVarType) {
214 case CONTINUOUS_INTERVAL_UNCERTAIN:
215 rv.pull_parameter(CIU_BPA, intervalBPA); break;
216 case DISCRETE_INTERVAL_UNCERTAIN:
217 rv.pull_parameter(DIU_BPA, intervalBPA); break;
218 default:
219 PCerr << "Error: update failure for RandomVariable type " << rv.type()
220 << " in IntervalRandomVariable::copy_parameters(T)." << std::endl;
221 abort_handler(-1); break;
222 }
223 update_vpp_if_active(); // update valueProbPairs if in use, else leave empty
224 }
225
226
227 template <typename T>
cdf(Real x) const228 Real IntervalRandomVariable<T>::cdf(Real x) const
229 {
230 //activate_vpp(); // can't do this since non-const
231
232 if (valueProbPairs.empty()) {
233 std::map<T, Real> xy_map; intervals_to_xy_pdf(intervalBPA, xy_map);
234 return DiscreteSetRandomVariable<T>::cdf(x, xy_map);
235 }
236 else
237 return DiscreteSetRandomVariable<T>::cdf(x, valueProbPairs);
238 }
239
240
241 template <typename T>
ccdf(Real x) const242 Real IntervalRandomVariable<T>::ccdf(Real x) const
243 {
244 //activate_vpp(); // can't do this since non-const
245
246 if (valueProbPairs.empty()) {
247 std::map<T, Real> xy_map; intervals_to_xy_pdf(intervalBPA, xy_map);
248 return DiscreteSetRandomVariable<T>::ccdf(x, xy_map);
249 }
250 else
251 return DiscreteSetRandomVariable<T>::ccdf(x, valueProbPairs);
252 }
253
254
255 template <typename T>
inverse_cdf(Real p_cdf) const256 Real IntervalRandomVariable<T>::inverse_cdf(Real p_cdf) const
257 {
258 //activate_vpp(); // can't do this since non-const
259
260 if (valueProbPairs.empty()) {
261 std::map<T, Real> xy_map; intervals_to_xy_pdf(intervalBPA, xy_map);
262 return DiscreteSetRandomVariable<T>::inverse_cdf(p_cdf, xy_map);
263 }
264 else
265 return DiscreteSetRandomVariable<T>::inverse_cdf(p_cdf, valueProbPairs);
266 }
267
268
269 template <typename T>
inverse_ccdf(Real p_ccdf) const270 Real IntervalRandomVariable<T>::inverse_ccdf(Real p_ccdf) const
271 {
272 //activate_vpp(); // can't do this since non-const
273
274 if (valueProbPairs.empty()) {
275 std::map<T, Real> xy_map; intervals_to_xy_pdf(intervalBPA, xy_map);
276 return DiscreteSetRandomVariable<T>::inverse_ccdf(p_ccdf, xy_map);
277 }
278 else
279 return DiscreteSetRandomVariable<T>::inverse_ccdf(p_ccdf, valueProbPairs);
280 }
281
282
283 template <typename T>
pdf(Real x) const284 Real IntervalRandomVariable<T>::pdf(Real x) const
285 {
286 //activate_vpp(); // can't do this since non-const
287
288 if (valueProbPairs.empty()) {
289 std::map<T, Real> xy_map; intervals_to_xy_pdf(intervalBPA, xy_map);
290 return DiscreteSetRandomVariable<T>::pdf(x, xy_map);
291 }
292 else
293 return DiscreteSetRandomVariable<T>::pdf(x, valueProbPairs);
294 }
295
296
297 template <typename T>
pdf_gradient(Real x) const298 Real IntervalRandomVariable<T>::pdf_gradient(Real x) const
299 { return 0.; }
300
301
302 template <typename T>
pdf_hessian(Real x) const303 Real IntervalRandomVariable<T>::pdf_hessian(Real x) const
304 { return 0.; }
305
306
307 template <typename T>
moments() const308 RealRealPair IntervalRandomVariable<T>::moments() const
309 {
310 //activate_vpp(); // can't do this since non-const
311
312 RealRealPair moms;
313 if (valueProbPairs.empty()) {
314 std::map<T, Real> xy_map; intervals_to_xy_pdf(intervalBPA, xy_map);
315 DiscreteSetRandomVariable<T>::
316 moments_from_params(xy_map, moms.first, moms.second);
317 }
318 else
319 DiscreteSetRandomVariable<T>::
320 moments_from_params(valueProbPairs, moms.first, moms.second);
321 return moms;
322 }
323
324
325 template <typename T>
mean() const326 Real IntervalRandomVariable<T>::mean() const
327 { return moments().first; }
328
329
330 template <typename T>
standard_deviation() const331 Real IntervalRandomVariable<T>::standard_deviation() const
332 { return moments().second; }
333
334
335 template <typename T>
variance() const336 Real IntervalRandomVariable<T>::variance() const
337 { Real std_dev = moments().second; return std_dev * std_dev; }
338
339
340 template <typename T>
mode() const341 Real IntervalRandomVariable<T>::mode() const
342 {
343 //activate_vpp(); // can't do this since non-const
344
345 if (valueProbPairs.empty()) {
346 std::map<T, Real> xy_map; intervals_to_xy_pdf(intervalBPA, xy_map);
347 return DiscreteSetRandomVariable<T>::mode(xy_map);
348 }
349 else
350 return DiscreteSetRandomVariable<T>::mode(valueProbPairs);
351 }
352
353
354 template <typename T>
distribution_bounds() const355 RealRealPair IntervalRandomVariable<T>::distribution_bounds() const
356 {
357 T l_bnd, u_bnd;
358 if (valueProbPairs.empty()) {
359 typename std::map<std::pair<T, T>, Real>::const_iterator
360 cit = intervalBPA.begin();
361 const std::pair<T, T>& bnds0 = cit->first;
362 l_bnd = bnds0.first; u_bnd = bnds0.second; ++cit;
363 for (; cit!=intervalBPA.end(); ++cit) {
364 const std::pair<T, T>& bnds = cit->first;
365 if (bnds.first < l_bnd) l_bnd = bnds.first;
366 if (bnds.second > u_bnd) u_bnd = bnds.second;
367 }
368 }
369 else {
370 l_bnd = valueProbPairs.begin()->first;
371 u_bnd = (--valueProbPairs.end())->first;
372 }
373
374 return RealRealPair((Real)l_bnd, (Real)u_bnd);
375 }
376
377
378 /// for T-valued histogram, return a real-valued mean and std dev
379 template <typename T>
380 void IntervalRandomVariable<T>::
moments_from_params(const std::map<std::pair<T,T>,Real> & bpa,Real & mean,Real & std_dev)381 moments_from_params(const std::map<std::pair<T, T>, Real>& bpa,
382 Real& mean, Real& std_dev)
383 {
384 typename std::map<T, Real> xy_map;
385 intervals_to_xy_pdf(bpa, xy_map);
386 DiscreteSetRandomVariable<T>::moments_from_params(xy_map, mean, std_dev);
387 }
388
389
390 //// SPECIALIZATIONS ////
391
392
393 template <>
cdf(Real x) const394 inline Real IntervalRandomVariable<Real>::cdf(Real x) const
395 {
396 //activate_vpp(); // can't do this since non-const
397
398 if (valueProbPairs.empty()) {
399 RealRealMap xy_map; intervals_to_xy_pdf(intervalBPA, xy_map);
400 return HistogramBinRandomVariable::cdf(x, xy_map);
401 }
402 else
403 return HistogramBinRandomVariable::cdf(x, valueProbPairs);
404 }
405
406
407 template <>
ccdf(Real x) const408 inline Real IntervalRandomVariable<Real>::ccdf(Real x) const
409 {
410 //activate_vpp(); // can't do this since non-const
411
412 if (valueProbPairs.empty()) {
413 RealRealMap xy_map; intervals_to_xy_pdf(intervalBPA, xy_map);
414 return HistogramBinRandomVariable::ccdf(x, xy_map);
415 }
416 else
417 return HistogramBinRandomVariable::ccdf(x, valueProbPairs);
418 }
419
420
421 template <>
inverse_cdf(Real x) const422 inline Real IntervalRandomVariable<Real>::inverse_cdf(Real x) const
423 {
424 //activate_vpp(); // can't do this since non-const
425
426 if (valueProbPairs.empty()) {
427 RealRealMap xy_map; intervals_to_xy_pdf(intervalBPA, xy_map);
428 return HistogramBinRandomVariable::inverse_cdf(x, xy_map);
429 }
430 else
431 return HistogramBinRandomVariable::inverse_cdf(x, valueProbPairs);
432 }
433
434
435 template <>
inverse_ccdf(Real x) const436 inline Real IntervalRandomVariable<Real>::inverse_ccdf(Real x) const
437 {
438 //activate_vpp(); // can't do this since non-const
439
440 if (valueProbPairs.empty()) {
441 RealRealMap xy_map; intervals_to_xy_pdf(intervalBPA, xy_map);
442 return HistogramBinRandomVariable::inverse_ccdf(x, xy_map);
443 }
444 else
445 return HistogramBinRandomVariable::inverse_ccdf(x, valueProbPairs);
446 }
447
448
449 template <>
pdf(Real x) const450 inline Real IntervalRandomVariable<Real>::pdf(Real x) const
451 {
452 //activate_vpp(); // can't do this since non-const
453
454 if (valueProbPairs.empty()) {
455 RealRealMap xy_map; intervals_to_xy_pdf(intervalBPA, xy_map);
456 return HistogramBinRandomVariable::pdf(x, xy_map);
457 }
458 else
459 return HistogramBinRandomVariable::pdf(x, valueProbPairs);
460 }
461
462
463 template <>
moments() const464 inline RealRealPair IntervalRandomVariable<Real>::moments() const
465 {
466 //activate_vpp(); // can't do this since non-const
467
468 RealRealPair moms;
469 if (valueProbPairs.empty()) {
470 RealRealMap xy_map; intervals_to_xy_pdf(intervalBPA, xy_map);
471 HistogramBinRandomVariable::
472 moments_from_params(xy_map, moms.first, moms.second);
473 }
474 else
475 HistogramBinRandomVariable::
476 moments_from_params(valueProbPairs, moms.first, moms.second);
477 return moms;
478 }
479
480
481 template <>
mode() const482 inline Real IntervalRandomVariable<Real>::mode() const
483 {
484 //activate_vpp(); // can't do this since non-const
485
486 if (valueProbPairs.empty()) {
487 RealRealMap xy_map; intervals_to_xy_pdf(intervalBPA, xy_map);
488 return HistogramBinRandomVariable::mode(xy_map);
489 }
490 else
491 return HistogramBinRandomVariable::mode(valueProbPairs);
492 }
493
494
495 template <>
496 inline void IntervalRandomVariable<Real>::
moments_from_params(const RealRealPairRealMap & bpa,Real & mean,Real & std_dev)497 moments_from_params(const RealRealPairRealMap& bpa, Real& mean, Real& std_dev)
498 {
499 RealRealMap xy_map; intervals_to_xy_pdf(bpa, xy_map);
500 HistogramBinRandomVariable::moments_from_params(xy_map, mean, std_dev);
501 }
502
503 } // namespace Pecos
504
505 #endif
506