1 // -*- C++ -*-
2 /**
3 * @brief The ProductDistribution distribution
4 *
5 * Copyright 2005-2021 Airbus-EDF-IMACS-ONERA-Phimeca
6 *
7 * This library is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU Lesser General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public License
18 * along with this library. If not, see <http://www.gnu.org/licenses/>.
19 *
20 */
21 #include <cstdlib>
22 #include <cmath>
23
24 #include "openturns/ProductDistribution.hxx"
25 #include "openturns/GaussKronrod.hxx"
26 #include "openturns/PersistentObjectFactory.hxx"
27 #include "openturns/UniVariateFunctionImplementation.hxx"
28 #include "openturns/Uniform.hxx"
29
30 BEGIN_NAMESPACE_OPENTURNS
31
32 CLASSNAMEINIT(ProductDistribution)
33
34 static const Factory<ProductDistribution> Factory_ProductDistribution;
35
36 /* Default constructor */
ProductDistribution()37 ProductDistribution::ProductDistribution()
38 : ContinuousDistribution()
39 , p_left_(new Uniform(0.0, 1.0))
40 , p_right_(new Uniform(0.0, 1.0))
41 , algo_()
42 {
43 setName("ProductDistribution");
44 setDimension(1);
45 // Adjust the truncation interval and the distribution range
46 computeRange();
47 }
48
49 /* Parameters constructor to use when the two bounds are finite */
ProductDistribution(const Distribution & left,const Distribution & right)50 ProductDistribution::ProductDistribution(const Distribution & left,
51 const Distribution & right)
52 : ContinuousDistribution()
53 , p_left_(left.getImplementation())
54 , p_right_(right.getImplementation())
55 , algo_()
56 {
57 setName("ProductDistribution");
58 setLeft(left);
59 setRight(right);
60 computeRange();
61 }
62
63 /* Comparison operator */
operator ==(const ProductDistribution & other) const64 Bool ProductDistribution::operator ==(const ProductDistribution & other) const
65 {
66 if (this == &other) return true;
67 return (p_left_ == other.getLeft().getImplementation()) && (p_right_ == other.getRight().getImplementation());
68 }
69
equals(const DistributionImplementation & other) const70 Bool ProductDistribution::equals(const DistributionImplementation & other) const
71 {
72 const ProductDistribution* p_other = dynamic_cast<const ProductDistribution*>(&other);
73 return p_other && (*this == *p_other);
74 }
75
76 /* String converter */
__repr__() const77 String ProductDistribution::__repr__() const
78 {
79 OSS oss;
80 oss << "class=" << ProductDistribution::GetClassName()
81 << " name=" << getName()
82 << " left=" << p_left_->__repr__()
83 << " right=" << p_right_->__repr__();
84 return oss;
85 }
86
__str__(const String &) const87 String ProductDistribution::__str__(const String & ) const
88 {
89 OSS oss;
90 oss << getClassName() << "(" << p_left_->__str__() << " * " << p_right_->__str__() << ")";
91 return oss;
92 }
93
94 /* Virtual constructor */
clone() const95 ProductDistribution * ProductDistribution::clone() const
96 {
97 return new ProductDistribution(*this);
98 }
99
100 /* Compute the numerical range of the distribution given the parameters values */
computeRange()101 void ProductDistribution::computeRange()
102 {
103 const Scalar a = p_left_->getRange().getLowerBound()[0];
104 const Scalar b = p_left_->getRange().getUpperBound()[0];
105 const Scalar c = p_right_->getRange().getLowerBound()[0];
106 const Scalar d = p_right_->getRange().getUpperBound()[0];
107 const Scalar ac = a * c;
108 const Scalar ad = a * d;
109 const Scalar bc = b * c;
110 const Scalar bd = b * d;
111 setRange(Interval(std::min(std::min(ac, ad), std::min(bc, bd)), std::max(std::max(ac, ad), std::max(bc, bd))));
112 }
113
114 /* Get one realization of the distribution */
getRealization() const115 Point ProductDistribution::getRealization() const
116 {
117 return Point(1, p_left_->getRealization()[0] * p_right_->getRealization()[0]);
118 }
119
120 namespace
121 {
122 // Class used to wrap the kernel of the integral defining the PDF of the product
123 class PDFKernelProductDistribution: public UniVariateFunctionImplementation
124 {
125 public:
PDFKernelProductDistribution(const Pointer<DistributionImplementation> & p_left,const Pointer<DistributionImplementation> & p_right,const Scalar x)126 PDFKernelProductDistribution(const Pointer<DistributionImplementation> & p_left,
127 const Pointer<DistributionImplementation> & p_right,
128 const Scalar x)
129 : UniVariateFunctionImplementation()
130 , p_left_(p_left)
131 , p_right_(p_right)
132 , x_(x)
133 , isZero_(std::abs(x) < ResourceMap::GetAsScalar("Distribution-DefaultQuantileEpsilon")), pdf0_(isZero_ ? p_right->computePDF(0.0) : 0.0)
134 {
135 // Nothing to do
136 };
137
clone() const138 PDFKernelProductDistribution * clone() const
139 {
140 return new PDFKernelProductDistribution(*this);
141 }
142
operator ()(const Scalar u) const143 Scalar operator() (const Scalar u) const
144 {
145 const Scalar value = p_left_->computePDF(u);
146 if (value == 0.0) return 0.0;
147 const Scalar absU = std::abs(u);
148 // x_ == 0
149 if (isZero_)
150 {
151 if (pdf0_ == 0.0) return 0.0;
152 if (absU == 0.0) return SpecFunc::MaxScalar;
153 return value * pdf0_ / absU;
154 }
155 // x_ != 0
156 if (absU == 0.0)
157 {
158 const Scalar epsilon = 1e-7;
159 return value * 0.5 * (p_right_->computePDF(x_ / epsilon) + p_right_->computePDF(-x_ / epsilon)) / epsilon;
160 }
161 return value * p_right_->computePDF(x_ / u) / absU;
162 };
163
164 private:
165 const Pointer<DistributionImplementation> p_left_;
166 const Pointer<DistributionImplementation> p_right_;
167 const Scalar x_;
168 const Bool isZero_;
169 const Scalar pdf0_;
170
171 }; // class PDFKernelProductDistribution
172
173 // Class used to wrap the kernel of the integral defining the CDF of the product
174 class CDFKernelProductDistribution: public UniVariateFunctionImplementation
175 {
176 public:
CDFKernelProductDistribution(const Pointer<DistributionImplementation> & p_left,const Pointer<DistributionImplementation> & p_right,const Scalar x)177 CDFKernelProductDistribution(const Pointer<DistributionImplementation> & p_left,
178 const Pointer<DistributionImplementation> & p_right,
179 const Scalar x)
180 : UniVariateFunctionImplementation()
181 , p_left_(p_left)
182 , p_right_(p_right)
183 , x_(x)
184 , isZero_(std::abs(x) < ResourceMap::GetAsScalar("Distribution-DefaultQuantileEpsilon"))
185 , cdf0_(isZero_ ? p_right->computeCDF(0.0) : 0.0)
186 , ccdf0_(isZero_ ? p_right->computeComplementaryCDF(0.0) : 0.0)
187 {
188 // Nothing to do
189 };
190
clone() const191 CDFKernelProductDistribution * clone() const
192 {
193 return new CDFKernelProductDistribution(*this);
194 }
195
operator ()(const Scalar u) const196 Scalar operator() (const Scalar u) const
197 {
198 const Scalar value = p_left_->computePDF(u);
199 if (value == 0.0) return 0.0;
200 // x_ == 0
201 if (isZero_) return value * cdf0_;
202 if (u == 0.0) return (x_ < 0.0 ? 0.0 : value);
203 return value * p_right_->computeCDF(x_ / u);
204 };
205
206 private:
207 const Pointer<DistributionImplementation> p_left_;
208 const Pointer<DistributionImplementation> p_right_;
209 const Scalar x_;
210 const Bool isZero_;
211 const Scalar cdf0_;
212 const Scalar ccdf0_;
213 }; // struct CDFKernelProductDistribution
214
215 // Class used to wrap the kernel of the integral defining the CDF of the product
216 class ComplementaryCDFKernelProductDistributionProductDistribution: public UniVariateFunctionImplementation
217 {
218 public:
ComplementaryCDFKernelProductDistributionProductDistribution(const Pointer<DistributionImplementation> & p_left,const Pointer<DistributionImplementation> & p_right,const Scalar x)219 ComplementaryCDFKernelProductDistributionProductDistribution(const Pointer<DistributionImplementation> & p_left,
220 const Pointer<DistributionImplementation> & p_right,
221 const Scalar x)
222 : UniVariateFunctionImplementation()
223 , p_left_(p_left)
224 , p_right_(p_right)
225 , x_(x)
226 , isZero_(std::abs(x) < ResourceMap::GetAsScalar("Distribution-DefaultQuantileEpsilon"))
227 , cdf0_(isZero_ ? p_right->computeCDF(0.0) : 0.0)
228 , ccdf0_(isZero_ ? p_right->computeComplementaryCDF(0.0) : 0.0)
229 {
230 // Nothing to do
231 };
232
clone() const233 ComplementaryCDFKernelProductDistributionProductDistribution * clone() const
234 {
235 return new ComplementaryCDFKernelProductDistributionProductDistribution(*this);
236 }
237
operator ()(const Scalar u) const238 Scalar operator() (const Scalar u) const
239 {
240 const Scalar value = p_left_->computePDF(u);
241 if (value == 0.0) return 0.0;
242 // x_ == 0
243 if (isZero_) return value * ccdf0_;
244 if (u == 0.0) return (x_ < 0.0 ? 0.0 : value);
245 return value * p_right_->computeComplementaryCDF(x_ / u);
246 };
247
248 private:
249 const Pointer<DistributionImplementation> p_left_;
250 const Pointer<DistributionImplementation> p_right_;
251 const Scalar x_;
252 const Bool isZero_;
253 const Scalar cdf0_;
254 const Scalar ccdf0_;
255 }; // struct ComplementaryCDFKernelProductDistributionProductDistribution
256
257 // Class used to wrap the kernel of the integral defining the CDF of the product
258 class CFKernelProductDistribution: public EvaluationImplementation
259 {
260 public:
CFKernelProductDistribution(const Pointer<DistributionImplementation> & p_left,const Pointer<DistributionImplementation> & p_right,const Scalar x)261 CFKernelProductDistribution(const Pointer<DistributionImplementation> & p_left,
262 const Pointer<DistributionImplementation> & p_right,
263 const Scalar x)
264 : EvaluationImplementation()
265 , p_left_(p_left)
266 , p_right_(p_right)
267 , x_(x)
268 {
269 // Nothing to do
270 };
271
clone() const272 CFKernelProductDistribution * clone() const
273 {
274 return new CFKernelProductDistribution(*this);
275 }
276
operator ()(const Point & point) const277 Point operator() (const Point & point) const
278 {
279 Point value(2);
280 const Scalar u = point[0];
281 const Complex phi(p_right_->computeCharacteristicFunction(u * x_));
282 const Scalar pdf = p_left_->computePDF(point);
283 value[0] = pdf * phi.real();
284 value[1] = pdf * phi.imag();
285 return value;
286 };
287
getInputDimension() const288 UnsignedInteger getInputDimension() const
289 {
290 return 1;
291 }
292
getOutputDimension() const293 UnsignedInteger getOutputDimension() const
294 {
295 return 2;
296 }
297
298 private:
299 const Pointer<DistributionImplementation> p_left_;
300 const Pointer<DistributionImplementation> p_right_;
301 const Scalar x_;
302 }; // struct CFKernelProductDistributionWrapper
303 } // anonymous namespace
304
305 /* Get the PDF of the distribution: PDF(x) = \int_R PDF_left(u) * PDF_right(x / u) * du / |u| */
computePDF(const Point & point) const306 Scalar ProductDistribution::computePDF(const Point & point) const
307 {
308 if (point.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=1, here dimension=" << point.getDimension();
309 return computePDF(point[0]);
310 }
311
computePDF(const Scalar x) const312 Scalar ProductDistribution::computePDF(const Scalar x) const
313 {
314 const Scalar a = getRange().getLowerBound()[0];
315 const Scalar b = getRange().getUpperBound()[0];
316 if ((x < a) || (x > b)) return 0.0;
317 const Scalar aLeft = p_left_->getRange().getLowerBound()[0];
318 const Scalar bLeft = p_left_->getRange().getUpperBound()[0];
319 const Scalar aRight = p_right_->getRange().getLowerBound()[0];
320 const Scalar bRight = p_right_->getRange().getUpperBound()[0];
321 // First, the case where the joint support of left and right is included in a unique quadrant
322 if ((aLeft >= 0.0) && (aRight >= 0.0))
323 {
324 LOGDEBUG("In ProductDistribution::computePDF, Q1");
325 const Scalar value = computePDFQ1(x, aLeft, bLeft, aRight, bRight);
326 LOGDEBUG(OSS() << "pdf=" << value);
327 return value;
328 }
329 if ((bLeft <= 0.0) && (aRight >= 0.0))
330 {
331 LOGDEBUG("In ProductDistribution::computePDF, Q2");
332 const Scalar value = computePDFQ2(x, aLeft, bLeft, aRight, bRight);
333 LOGDEBUG(OSS() << "pdf=" << value);
334 return value;
335 }
336 if ((bLeft <= 0.0) && (bRight <= 0.0))
337 {
338 LOGDEBUG("In ProductDistribution::computePDF, Q3");
339 const Scalar value = computePDFQ3(x, aLeft, bLeft, aRight, bRight);
340 LOGDEBUG(OSS() << "pdf=" << value);
341 return value;
342 }
343 if ((aLeft >= 0.0) && (bRight <= 0.0))
344 {
345 LOGDEBUG("In ProductDistribution::computePDF, Q3");
346 const Scalar value = computePDFQ4(x, aLeft, bLeft, aRight, bRight);
347 LOGDEBUG(OSS() << "pdf=" << value);
348 return value;
349 }
350 // Second, the case where the support is in Q1 U Q2
351 if (aRight > 0.0)
352 {
353 LOGDEBUG("In ProductDistribution::computePDF, Q1 U Q2");
354 const Scalar q1 = computePDFQ1(x, 0.0, bLeft, aRight, bRight);
355 const Scalar q2 = computePDFQ2(x, aLeft, 0.0, aRight, bRight);
356 LOGDEBUG(OSS() << "value Q1=" << q1 << ", value Q2=" << q2 << ", pdf=" << q1 + q2);
357 return q1 + q2;
358 }
359 // Third, the case where the support is in Q3 U Q4
360 if (bRight <= 0.0)
361 {
362 LOGDEBUG("In ProductDistribution::computePDF, Q3 U Q4");
363 const Scalar q3 = computePDFQ3(x, aLeft, 0.0, aRight, bRight);
364 const Scalar q4 = computePDFQ4(x, 0.0, bLeft, aRight, bRight);
365 LOGDEBUG(OSS() << "value Q3=" << q3 << ", value Q4=" << q4 << ", pdf=" << q3 + q4);
366 return q3 + q4;
367 }
368 // Fourth, the case where the support is in Q1 U Q4
369 if (aLeft >= 0.0)
370 {
371 LOGDEBUG("In ProductDistribution::computePDF, Q1 U Q4");
372 const Scalar q1 = computePDFQ1(x, aLeft, bLeft, 0.0, bRight);
373 const Scalar q4 = computePDFQ4(x, aLeft, bLeft, aRight, 0.0);
374 LOGDEBUG(OSS() << "value Q1=" << q1 << ", value Q4=" << q4 << ", pdf=" << q1 + q4);
375 return q1 + q4;
376 }
377 // Fifth, the case where the support is in Q2 U Q3
378 if (bLeft <= 0.0)
379 {
380 LOGDEBUG("In ProductDistribution::computePDF, Q2 U Q3");
381 const Scalar q2 = computePDFQ2(x, aLeft, bLeft, 0.0, bRight);
382 const Scalar q3 = computePDFQ3(x, aLeft, bLeft, aRight, 0.0);
383 LOGDEBUG(OSS() << "value Q2=" << q2 << ", value Q3=" << q3 << ", pdf=" << q2 + q3);
384 return q2 + q3;
385 }
386 // Sixth, the case where the support is in Q1 U Q2 U Q3 U Q4
387 LOGDEBUG("In ProductDistribution::computePDF, Q1 U Q2 U Q3 U Q4");
388 const Scalar q1 = computePDFQ1(x, 0.0, bLeft, 0.0, bRight);
389 const Scalar q2 = computePDFQ2(x, aLeft, 0.0, 0.0, bRight);
390 const Scalar q3 = computePDFQ3(x, aLeft, 0.0, aRight, 0.0);
391 const Scalar q4 = computePDFQ4(x, 0.0, bLeft, aRight, 0.0);
392 LOGDEBUG(OSS() << "value Q1=" << q1 << "value Q2=" << q2 << ", value Q3=" << q3 << "value Q4=" << q4 << ", pdf=" << q1 + q2 + q3 + q4);
393 return q1 + q2 + q3 + q4;
394 }
395
396 /* Get the PDF of the distribution: PDF(x) = \int_Q1 PDF_left(u) * PDF_right(x / u) * du / |u| when left >= 0, right >= 0 */
computePDFQ1(const Scalar x,const Scalar a,const Scalar b,const Scalar c,const Scalar d) const397 Scalar ProductDistribution::computePDFQ1(const Scalar x,
398 const Scalar a,
399 const Scalar b,
400 const Scalar c,
401 const Scalar d) const
402 {
403 LOGDEBUG(OSS() << "In ProductDistribution::computePDFQ1, x=" << x << ", a=" << a << ", b=" << b << ", c=" << c << ", d=" << d);
404 const Scalar ac = a * c;
405 const Scalar ad = a * d;
406 const Scalar bc = b * c;
407 const Scalar bd = b * d;
408 LOGDEBUG(OSS() << "ac=" << ac << ", ad=" << ad << ", bc=" << bc << ", bd=" << bd);
409 // Here the support is included into [ac, bd]
410 if ((x < ac) || (x >= bd)) return 0.0;
411 const PDFKernelProductDistribution pdfKernel(p_left_, p_right_, x);
412 if (c == 0.0)
413 {
414 LOGDEBUG("c == 0.0");
415 if (x < ad)
416 {
417 LOGDEBUG(OSS() << x << " < " << ad);
418 const Scalar value = algo_.integrate(pdfKernel, a, b);
419 LOGDEBUG(OSS() << "value=" << value);
420 return value;
421 }
422 LOGDEBUG(OSS() << ad << " <= " << x);
423 const Scalar value = algo_.integrate(pdfKernel, x / d, b);
424 LOGDEBUG(OSS() << "value=" << value);
425 return value;
426 }
427 if (ad <= bc)
428 {
429 LOGDEBUG("ad <= bc");
430 if (x < ad)
431 {
432 LOGDEBUG(OSS() << x << " < " << ad);
433 const Scalar value = algo_.integrate(pdfKernel, a, x / c);
434 LOGDEBUG(OSS() << "value=" << value);
435 return value;
436 }
437 if (x < bc)
438 {
439 LOGDEBUG(OSS() << x << " < " << bc);
440 const Scalar value = algo_.integrate(pdfKernel, x / d, x / c);
441 LOGDEBUG(OSS() << "value=" << value);
442 return value;
443 }
444 LOGDEBUG(OSS() << x << " < " << bd);
445 const Scalar value = algo_.integrate(pdfKernel, x / d, b);
446 LOGDEBUG(OSS() << "value=" << value);
447 return value;
448 }
449 LOGDEBUG("ad > bc");
450 if (x < bc)
451 {
452 LOGDEBUG(OSS() << x << " < " << bc);
453 const Scalar value = algo_.integrate(pdfKernel, a, x / c);
454 LOGDEBUG(OSS() << "value=" << value);
455 return value;
456 }
457 if (x < ad)
458 {
459 LOGDEBUG(OSS() << x << " < " << ad);
460 const Scalar value = algo_.integrate(pdfKernel, a, b);
461 LOGDEBUG(OSS() << "value=" << value);
462 return value;
463 }
464 LOGDEBUG(OSS() << x << " < " << bd);
465 const Scalar value = algo_.integrate(pdfKernel, x / d, b);
466 LOGDEBUG(OSS() << "value=" << value);
467 return value;
468 }
469
470 /* Get the PDF of the distribution: PDF(x) = \int_Q2 PDF_left(u) * PDF_right(x / u) * du / |u| when left <= 0, right >= 0 */
computePDFQ2(const Scalar x,const Scalar a,const Scalar b,const Scalar c,const Scalar d) const471 Scalar ProductDistribution::computePDFQ2(const Scalar x,
472 const Scalar a,
473 const Scalar b,
474 const Scalar c,
475 const Scalar d) const
476 {
477 LOGDEBUG(OSS() << "In ProductDistribution::computePDFQ2, x=" << x << ", a=" << a << ", b=" << b << ", c=" << c << ", d=" << d);
478 const Scalar ac = a * c;
479 const Scalar ad = a * d;
480 const Scalar bc = b * c;
481 const Scalar bd = b * d;
482 LOGDEBUG(OSS() << "ac=" << ac << ", ad=" << ad << ", bc=" << bc << ", bd=" << bd);
483 // Here the support is included into [ad, bc]
484 if ((x < ad) || (x >= bc)) return 0.0;
485 const PDFKernelProductDistribution pdfKernel(p_left_, p_right_, x);
486 if (c == 0.0)
487 {
488 LOGDEBUG("c == 0.0");
489 if (x < bd)
490 {
491 LOGDEBUG(OSS() << x << " < " << bd);
492 const Scalar value = algo_.integrate(pdfKernel, a, x / d);
493 LOGDEBUG(OSS() << "value=" << value);
494 return value;
495 }
496 LOGDEBUG(OSS() << x << " < " << 0.0);
497 const Scalar value = algo_.integrate(pdfKernel, a, b);
498 LOGDEBUG(OSS() << "value=" << value);
499 return value;
500 }
501 if (ac <= bd)
502 {
503 LOGDEBUG("ac <= bd");
504 if (x < ac)
505 {
506 LOGDEBUG(OSS() << x << " < " << ac);
507 const Scalar value = algo_.integrate(pdfKernel, a, x / d);
508 LOGDEBUG(OSS() << "value=" << value);
509 return value;
510 }
511 if (x < bd)
512 {
513 LOGDEBUG(OSS() << x << " < " << bd);
514 const Scalar value = algo_.integrate(pdfKernel, x / c, x / d);
515 LOGDEBUG(OSS() << "value=" << value);
516 return value;
517 }
518 LOGDEBUG(OSS() << x << " < " << bc);
519 const Scalar value = algo_.integrate(pdfKernel, x / c, b);
520 LOGDEBUG(OSS() << "value=" << value);
521 return value;
522 }
523 LOGDEBUG("ac > bd");
524 if (x < bd)
525 {
526 LOGDEBUG(OSS() << x << " < " << bd);
527 const Scalar value = algo_.integrate(pdfKernel, a, x / d);
528 LOGDEBUG(OSS() << "value=" << value);
529 return value;
530 }
531 if (x < ac)
532 {
533 LOGDEBUG(OSS() << x << " < " << ac);
534 const Scalar value = algo_.integrate(pdfKernel, a, b);
535 LOGDEBUG(OSS() << "value=" << value);
536 return value;
537 }
538 LOGDEBUG(OSS() << x << " < " << bc);
539 const Scalar value = algo_.integrate(pdfKernel, x / c, b);
540 LOGDEBUG(OSS() << "value=" << value);
541 return value;
542 }
543
544 /* Get the PDF of the distribution: PDF(x) = \int_Q3 PDF_left(u) * PDF_right(x / u) * du / |u| when left <= 0, right <= 0 */
computePDFQ3(const Scalar x,const Scalar a,const Scalar b,const Scalar c,const Scalar d) const545 Scalar ProductDistribution::computePDFQ3(const Scalar x,
546 const Scalar a,
547 const Scalar b,
548 const Scalar c,
549 const Scalar d) const
550 {
551 LOGDEBUG(OSS() << "In ProductDistribution::computePDFQ3, x=" << x << ", a=" << a << ", b=" << b << ", c=" << c << ", d=" << d);
552 const Scalar ac = a * c;
553 const Scalar ad = a * d;
554 const Scalar bc = b * c;
555 const Scalar bd = b * d;
556 LOGDEBUG(OSS() << "ac=" << ac << ", ad=" << ad << ", bc=" << bc << ", bd=" << bd);
557 // Here the support is included into [bd, ac]
558 if ((x < bd) || (x >= ac)) return 0.0;
559 const PDFKernelProductDistribution pdfKernel(p_left_, p_right_, x);
560 if (d == 0.0)
561 {
562 LOGDEBUG("d == 0.0");
563 if (x < bc)
564 {
565 LOGDEBUG(OSS() << x << " < " << bc);
566 const Scalar value = algo_.integrate(pdfKernel, a, b);
567 LOGDEBUG(OSS() << "value=" << value);
568 return value;
569 }
570 LOGDEBUG(OSS() << x << " < " << ac);
571 const Scalar value = algo_.integrate(pdfKernel, a, x / c);
572 LOGDEBUG(OSS() << "value=" << value);
573 return value;
574 }
575 if (ad <= bc)
576 {
577 LOGDEBUG("ad <= bc");
578 if (x < ad)
579 {
580 LOGDEBUG(OSS() << x << " < " << ad);
581 const Scalar value = algo_.integrate(pdfKernel, x / d, b);
582 LOGDEBUG(OSS() << "value=" << value);
583 return value;
584 }
585 if (x < bc)
586 {
587 LOGDEBUG(OSS() << x << " < " << bc);
588 const Scalar value = algo_.integrate(pdfKernel, a, b);
589 LOGDEBUG(OSS() << "value=" << value);
590 return value;
591 }
592 LOGDEBUG(OSS() << x << " < " << ac);
593 const Scalar value = algo_.integrate(pdfKernel, a, x / c);
594 LOGDEBUG(OSS() << "value=" << value);
595 return value;
596 }
597 LOGDEBUG("ad > bc");
598 if (x < bc)
599 {
600 LOGDEBUG(OSS() << x << " < " << bc);
601 const Scalar value = algo_.integrate(pdfKernel, x / d, b);
602 LOGDEBUG(OSS() << "value=" << value);
603 return value;
604 }
605 if (x < ad)
606 {
607 LOGDEBUG(OSS() << x << " < " << ad);
608 const Scalar value = algo_.integrate(pdfKernel, x / d, x / c);
609 LOGDEBUG(OSS() << "value=" << value);
610 return value;
611 }
612 LOGDEBUG(OSS() << x << " < " << ac);
613 const Scalar value = algo_.integrate(pdfKernel, a, x / c);
614 LOGDEBUG(OSS() << "value=" << value);
615 return value;
616 }
617
618 /* Get the PDF of the distribution: PDF(x) = \int_Q4 PDF_left(u) * PDF_right(x / u) * du / |u| when left >= 0, right <= 0 */
computePDFQ4(const Scalar x,const Scalar a,const Scalar b,const Scalar c,const Scalar d) const619 Scalar ProductDistribution::computePDFQ4(const Scalar x,
620 const Scalar a,
621 const Scalar b,
622 const Scalar c,
623 const Scalar d) const
624 {
625 LOGDEBUG(OSS() << "In ProductDistribution::computePDFQ4, x=" << x << ", a=" << a << ", b=" << b << ", c=" << c << ", d=" << d);
626 const Scalar ac = a * c;
627 const Scalar ad = a * d;
628 const Scalar bc = b * c;
629 const Scalar bd = b * d;
630 LOGDEBUG(OSS() << "ac=" << ac << ", ad=" << ad << ", bc=" << bc << ", bd=" << bd);
631 // Here the support is included into [bc, ad]
632 if ((x < bc) || (x >= ad)) return 0.0;
633 const PDFKernelProductDistribution pdfKernel(p_left_, p_right_, x);
634 if (d == 0.0)
635 {
636 LOGDEBUG("d == 0.0");
637 if (x < ac)
638 {
639 LOGDEBUG(OSS() << x << " < " << ac);
640 const Scalar value = algo_.integrate(pdfKernel, x / c, b);
641 LOGDEBUG(OSS() << "value=" << value);
642 return value;
643 }
644 LOGDEBUG(OSS() << x << " < " << 0.0);
645 const Scalar value = algo_.integrate(pdfKernel, a, b);
646 LOGDEBUG(OSS() << "value=" << value);
647 return value;
648 }
649 if (bd <= ac)
650 {
651 LOGDEBUG("bd <= ac");
652 if (x < bd)
653 {
654 LOGDEBUG(OSS() << x << " < " << bd);
655 const Scalar value = algo_.integrate(pdfKernel, x / c, b);
656 LOGDEBUG(OSS() << "value=" << value);
657 return value;
658 }
659 if (x < ac)
660 {
661 LOGDEBUG(OSS() << x << " < " << ac);
662 const Scalar value = algo_.integrate(pdfKernel, x / c, x / d);
663 LOGDEBUG(OSS() << "value=" << value);
664 return value;
665 }
666 LOGDEBUG(OSS() << x << " < " << ad);
667 const Scalar value = algo_.integrate(pdfKernel, a, x / d);
668 LOGDEBUG(OSS() << "value=" << value);
669 return value;
670 }
671 LOGDEBUG("In ProductDistribution::computePDFQ4, bd > ac");
672 LOGDEBUG("bd > ac");
673 if (x < ac)
674 {
675 LOGDEBUG(OSS() << x << " < " << ac);
676 const Scalar value = algo_.integrate(pdfKernel, x / c, b);
677 LOGDEBUG(OSS() << "value=" << value);
678 return value;
679 }
680 if (x < bd)
681 {
682 LOGDEBUG(OSS() << x << " < " << bd);
683 const Scalar value = algo_.integrate(pdfKernel, a, b);
684 LOGDEBUG(OSS() << "value=" << value);
685 return value;
686 }
687 LOGDEBUG(OSS() << bd << " <= " << x << " < " << ad);
688 const Scalar value = algo_.integrate(pdfKernel, a, x / d);
689 LOGDEBUG(OSS() << "value=" << value);
690 return value;
691 }
692
693 /* Get the CDF of the distribution */
computeCDF(const Point & point) const694 Scalar ProductDistribution::computeCDF(const Point & point) const
695 {
696 if (point.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=1, here dimension=" << point.getDimension();
697
698 return computeCDF(point[0]);
699 }
700
computeCDF(const Scalar x) const701 Scalar ProductDistribution::computeCDF(const Scalar x) const
702 {
703 const Scalar a = getRange().getLowerBound()[0];
704 const Scalar b = getRange().getUpperBound()[0];
705 if (x <= a) return 0.0;
706 if (x >= b) return 1.0;
707 const Scalar aLeft = p_left_->getRange().getLowerBound()[0];
708 const Scalar bLeft = p_left_->getRange().getUpperBound()[0];
709 Scalar value = 0.0;
710 // First, compute the negative part
711 if (aLeft < 0)
712 {
713 const ComplementaryCDFKernelProductDistributionProductDistribution cdfKernel(p_left_, p_right_, x);
714 value += algo_.integrate(cdfKernel, aLeft, std::min(bLeft, 0.0));
715 }
716 if (bLeft >= 0)
717 {
718 const CDFKernelProductDistribution cdfKernel(p_left_, p_right_, x);
719 value += algo_.integrate(cdfKernel, std::max(0.0, aLeft), bLeft);
720 }
721 return value;
722 }
723
724 /* Compute the probability content of an interval */
computeProbability(const Interval & interval) const725 Scalar ProductDistribution::computeProbability(const Interval & interval) const
726 {
727 return computeProbabilityContinuous(interval);
728 }
729
730 /* Get the characteristic function of the distribution, i.e. phi(u) = E(exp(I*u*X)) */
computeCharacteristicFunction(const Scalar x) const731 Complex ProductDistribution::computeCharacteristicFunction(const Scalar x) const
732 {
733 const Scalar muLeft = p_left_->getMean()[0];
734 const Scalar muRight = p_right_->getMean()[0];
735 const Scalar varLeft = p_left_->getCovariance()(0, 0);
736 const Scalar varRight = p_right_->getCovariance()(0, 0);
737 if (x * x * (varLeft + muLeft * muLeft + varRight + muRight * muRight) < 2.0 * SpecFunc::ScalarEpsilon) return Complex(1.0, -x * muLeft * muRight);
738 if (std::abs(x) > ResourceMap::GetAsScalar("ProductDistribution-LargeCharacteristicFunctionArgument")) return ContinuousDistribution::computeCharacteristicFunction(x);
739 const Scalar aLeft = p_left_->getRange().getLowerBound()[0];
740 const Scalar bLeft = p_left_->getRange().getUpperBound()[0];
741 const CFKernelProductDistribution cfKernel(p_left_, p_right_, x);
742 Scalar negativeError = 0.0;
743 const Point negativePart(algo_.integrate(cfKernel, Interval(aLeft, muLeft), negativeError));
744 Scalar positiveError = 0.0;
745 const Point positivePart(algo_.integrate(cfKernel, Interval(muLeft, bLeft), positiveError));
746 Complex value(negativePart[0] + positivePart[0], negativePart[1] + positivePart[1]);
747 return value;
748 }
749
750 /* Compute the mean of the distribution */
computeMean() const751 void ProductDistribution::computeMean() const
752 {
753 mean_ = Point(1, p_left_->getMean()[0] * p_right_->getMean()[0]);
754 isAlreadyComputedMean_ = true;
755 }
756
757 /* Compute the covariance of the distribution */
computeCovariance() const758 void ProductDistribution::computeCovariance() const
759 {
760 covariance_ = CovarianceMatrix(1);
761 const Scalar meanLeft = p_left_->getMean()[0];
762 const Scalar meanRight = p_right_->getMean()[0];
763 const Scalar varLeft = p_left_->getCovariance()(0, 0);
764 const Scalar varRight = p_right_->getCovariance()(0, 0);
765 covariance_(0, 0) = meanLeft * meanLeft * varRight + meanRight * meanRight * varLeft + varLeft * varRight;
766 isAlreadyComputedCovariance_ = true;
767 }
768
769 /* Get the skewness of the distribution */
getSkewness() const770 Point ProductDistribution::getSkewness() const
771 {
772 const Scalar meanLeft = p_left_->getMean()[0];
773 const Scalar meanRight = p_right_->getMean()[0];
774 const Scalar varLeft = p_left_->getCovariance()(0, 0);
775 const Scalar varRight = p_right_->getCovariance()(0, 0);
776 const Scalar mu3Left = p_left_->getSkewness()[0] * std::pow(varLeft, 1.5);
777 const Scalar mu3Right = p_right_->getSkewness()[0] * std::pow(varRight, 1.5);
778 const Scalar variance = meanLeft * meanLeft * varRight + meanRight * meanRight * varLeft + varLeft * varRight;
779 return Point(1, (mu3Left * mu3Right + mu3Left * std::pow(meanRight, 3.0) + mu3Right * std::pow(meanLeft, 3.0) + 3.0 * (mu3Left * varRight * meanRight + mu3Right * varLeft * meanLeft) + 6.0 * varLeft * varRight * meanLeft * meanRight) / std::pow(variance, 1.5));
780 }
781
782 /* Get the kurtosis of the distribution */
getKurtosis() const783 Point ProductDistribution::getKurtosis() const
784 {
785 const Scalar meanLeft = p_left_->getMean()[0];
786 const Scalar meanLeft2 = meanLeft * meanLeft;
787 const Scalar meanLeft4 = meanLeft2 * meanLeft2;
788 const Scalar meanRight = p_right_->getMean()[0];
789 const Scalar meanRight2 = meanRight * meanRight;
790 const Scalar meanRight4 = meanRight2 * meanRight2;
791 const Scalar varLeft = p_left_->getCovariance()(0, 0);
792 const Scalar varRight = p_right_->getCovariance()(0, 0);
793 const Scalar mu3Left = p_left_->getSkewness()[0] * std::pow(varLeft, 1.5);
794 const Scalar mu3Right = p_right_->getSkewness()[0] * std::pow(varRight, 1.5);
795 const Scalar mu4Left = p_left_->getKurtosis()[0] * varLeft * varLeft;
796 const Scalar mu4Right = p_right_->getKurtosis()[0] * varRight * varRight;
797 const Scalar variance = meanLeft * meanLeft * varRight + meanRight * meanRight * varLeft + varLeft * varRight;
798 return Point(1, (mu4Left * mu4Right + mu4Left * meanRight4 + mu4Right * meanLeft4 + 4.0 * (mu4Left * mu3Right * meanRight + mu4Right * mu3Left * meanLeft) + 6.0 * (varLeft * meanLeft2 * varRight * meanRight2 + mu4Left * varRight * meanRight2 + mu4Right * varLeft * meanLeft2) + 12.0 * (mu3Left * meanLeft * mu3Right * meanRight + mu3Left * meanLeft * varRight * meanRight2 + mu3Right * meanRight * varLeft * meanLeft2)) / (variance * variance));
799 }
800
801 /* Get the raw moments of the distribution */
getMoment(const UnsignedInteger n) const802 Point ProductDistribution::getMoment(const UnsignedInteger n) const
803 {
804 return Point(1, p_left_->getMoment(n)[0] * p_right_->getMoment(n)[0]);
805 }
806
807
808 /* Parameters value accessor */
getParameter() const809 Point ProductDistribution::getParameter() const
810 {
811 Point point(p_left_->getParameter());
812 point.add(p_right_->getParameter());
813 return point;
814 }
815
setParameter(const Point & parameter)816 void ProductDistribution::setParameter(const Point & parameter)
817 {
818 const UnsignedInteger leftSize = p_left_->getParameterDimension();
819 const UnsignedInteger rightSize = p_right_->getParameterDimension();
820 if (parameter.getSize() != leftSize + rightSize)
821 throw InvalidArgumentException(HERE) << "Error: expected " << leftSize + rightSize << " values, got " << parameter.getSize();
822 Point newLeftParameters(leftSize);
823 Point newRightParameters(rightSize);
824 std::copy(parameter.begin(), parameter.begin() + leftSize, newLeftParameters.begin());
825 std::copy(parameter.begin() + leftSize, parameter.end(), newRightParameters.begin());
826 Distribution newLeft(p_left_);
827 Distribution newRight(p_right_);
828 newLeft.setParameter(newLeftParameters);
829 newRight.setParameter(newRightParameters);
830 const Scalar w = getWeight();
831 *this = ProductDistribution(newLeft, newRight);
832 setWeight(w);
833 }
834
835 /* Parameters description accessor */
getParameterDescription() const836 Description ProductDistribution::getParameterDescription() const
837 {
838 Description description(p_left_->getParameterDescription());
839 description.add(p_right_->getParameterDescription());
840 return description;
841 }
842
843 /* Check if the distribution is elliptical */
isElliptical() const844 Bool ProductDistribution::isElliptical() const
845 {
846 return (p_left_->isElliptical() && (std::abs(p_left_->getRange().getLowerBound()[0] + p_left_->getRange().getUpperBound()[0]) < ResourceMap::GetAsScalar("Distribution-DefaultQuantileEpsilon"))) || (p_right_->isElliptical() && (std::abs(p_right_->getRange().getLowerBound()[0] + p_right_->getRange().getUpperBound()[0]) < ResourceMap::GetAsScalar("Distribution-DefaultQuantileEpsilon")));
847 }
848
849 /* Left accessor */
setLeft(const Distribution & left)850 void ProductDistribution::setLeft(const Distribution & left)
851 {
852 if (left.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: can multiply only distribution with dimension=1, here dimension=" << left.getDimension();
853 if (!left.isContinuous()) throw InvalidArgumentException(HERE) << "Error: can multiply only continuous distributions";
854 p_left_ = left.getImplementation();
855 isAlreadyComputedMean_ = false;
856 isAlreadyComputedCovariance_ = false;
857 isAlreadyCreatedGeneratingFunction_ = false;
858 isParallel_ = p_left_->isParallel();
859 computeRange();
860 }
861
getLeft() const862 Distribution ProductDistribution::getLeft() const
863 {
864 return p_left_;
865 }
866
867 /* Right accessor */
setRight(const Distribution & right)868 void ProductDistribution::setRight(const Distribution & right)
869 {
870 if (right.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: can multiply only distribution with dimension=1, here dimension=" << right.getDimension();
871 if (!right.isContinuous()) throw InvalidArgumentException(HERE) << "Error: can multiply only continuous distributions";
872 p_right_ = right.getImplementation();
873 isAlreadyComputedMean_ = false;
874 isAlreadyComputedCovariance_ = false;
875 isAlreadyCreatedGeneratingFunction_ = false;
876 isParallel_ = p_right_->isParallel();
877 computeRange();
878 }
879
getRight() const880 Distribution ProductDistribution::getRight() const
881 {
882 return p_right_;
883 }
884
isContinuous() const885 Bool ProductDistribution::isContinuous() const
886 {
887 return p_left_->isContinuous() && p_right_->isContinuous();
888 }
889
890 /* Tell if the distribution is discrete */
isDiscrete() const891 Bool ProductDistribution::isDiscrete() const
892 {
893 return p_left_->isDiscrete() && p_right_->isDiscrete();
894 }
895
896 /* Tell if the distribution is integer valued */
isIntegral() const897 Bool ProductDistribution::isIntegral() const
898 {
899 return p_left_->isIntegral() && p_right_->isIntegral();
900 }
901
902 /* Get the PDF singularities inside of the range - 1D only */
getSingularities() const903 Point ProductDistribution::getSingularities() const
904 {
905 if (getRange().getLowerBound()[0] >= 0.0) return Point(0);
906 if (getRange().getUpperBound()[0] <= 0.0) return Point(0);
907 return Point(1);
908 }
909
910 /* Method save() stores the object through the StorageManager */
save(Advocate & adv) const911 void ProductDistribution::save(Advocate & adv) const
912 {
913 ContinuousDistribution::save(adv);
914 adv.saveAttribute( "left_", *p_left_ );
915 adv.saveAttribute( "right_", *p_right_ );
916 }
917
918 /* Method load() reloads the object from the StorageManager */
load(Advocate & adv)919 void ProductDistribution::load(Advocate & adv)
920 {
921 ContinuousDistribution::load(adv);
922 Distribution left;
923 adv.loadAttribute( "left_", left );
924 p_left_ = left.getImplementation();
925 Distribution right;
926 adv.loadAttribute( "right_", right );
927 p_right_ = right.getImplementation();
928 computeRange();
929 }
930
931 END_NAMESPACE_OPENTURNS
932