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