Improved MarginalDistribution

Now it is possible to switch between a CDF-based approach to compute the marginal PDF or a PDF-based approach (by default). It results into a huge speed-up of the PDF computations, see https://openturns.discourse.group/t/skew-distributions-in-ot/348/2.
It fixes openturns#2067.

lib/etc/openturns.conf.in

@@ -542,6 +542,9 @@
   <LogNormalFactory-MaximumIteration  value_int="50" />
   <LogNormalFactory-EstimationMethod  value_int="0"  />
 
+  <!-- MarginalDistribution parameters -->
+  <MarginalDistribution-UsePDF      value_bool="true"         />
+
   <!-- Meixner parameters -->
   <MeixnerDistribution-MaximumAbsoluteError      value_float="1.0e-12"         />
   <MeixnerDistribution-MaximumConstraintError    value_float="1.0e-12"         />

lib/src/Base/Common/ResourceMap.cxx

@@ -1169,6 +1169,9 @@ void ResourceMap::loadDefaultConfiguration()
   addAsUnsignedInteger("LogNormalFactory-EstimationMethod", 0);
   addAsUnsignedInteger("LogNormalFactory-MaximumIteration", 50);
 
+  // MarginalDistribution parameters //
+  addAsBool("MarginalDistribution-UsePDF", true);
+
   // Meixner parameters //
   addAsScalar("MeixnerDistribution-MaximumAbsoluteError", 1.0e-12);
   addAsScalar("MeixnerDistribution-MaximumConstraintError", 1.0e-12);

lib/src/Uncertainty/Distribution/MarginalDistribution.cxx

@@ -22,6 +22,9 @@
 
 #include "openturns/MarginalDistribution.hxx"
 #include "openturns/Uniform.hxx"
+#include "openturns/ParametricFunction.hxx"
+#include "openturns/IteratedQuadrature.hxx"
+#include "openturns/GaussKronrod.hxx"
 #include "openturns/PersistentObjectFactory.hxx"
 
 BEGIN_NAMESPACE_OPENTURNS
@@ -141,10 +144,10 @@ void MarginalDistribution::setDistributionAndIndices(const Distribution & distri
   setDimension(dimension);
   // Compute the range
   // From the underlying distribution
-  Point distributionLowerBound(distribution.getRange().getLowerBound());
+  const Point distributionLowerBound(distribution.getRange().getLowerBound());
   lowerBound_ = distributionLowerBound;
   Interval::BoolCollection distributionFiniteLowerBound(distribution.getRange().getFiniteLowerBound());
-  Point distributionUpperBound(distribution.getRange().getUpperBound());
+  const Point distributionUpperBound(distribution.getRange().getUpperBound());
   upperBound_ = distributionUpperBound;
   Interval::BoolCollection distributionFiniteUpperBound(distribution.getRange().getFiniteUpperBound());
   // For the marginal distribution
@@ -192,6 +195,37 @@ Sample MarginalDistribution::getSample(const UnsignedInteger size) const
   return distribution_.getSample(size).getMarginal(indices_);
 }
 
+/* Get the PDF of the MarginalDistribution */
+Scalar MarginalDistribution::computePDF(const Point & point) const
+{
+  if (point.getDimension() != getDimension()) throw InvalidArgumentException(HERE) << "Error: expected a point of dimension=" << getDimension() << ", got dimension=" << point.getDimension();
+  // First, check if the marginal distribution is just a eordering of the underlying distribution
+  if (distribution_.getDimension() == getDimension())
+    return distribution_.computePDF(expandPoint(point));
+  if (usePDF_)
+    {
+      if (isContinuous())
+        {
+          // Build the relevant parametric function to be integrated over the remaining parameters
+          const ParametricFunction kernel(PDFWrapper(distribution_.getImplementation()->clone()), indices_, point);
+          const Interval marginalInterval(distribution_.getRange().getMarginal(indices_.complement(distribution_.getDimension())));
+          const Scalar pdf = IteratedQuadrature(GaussKronrod()).integrate(kernel, marginalInterval)[0];
+          return pdf;
+        }
+      if (isDiscrete())
+        {
+          const Point probabilities(distribution_.getProbabilities());
+          // We modify the support in-place to speed-up the computation
+          Sample support(distribution_.getImplementation()->getSupport());
+          for (UnsignedInteger i = 0; i < support.getSize(); ++i)
+              for (UnsignedInteger j = 0; j < indices_.getSize(); ++j)
+                support(i, indices_[j]) = point[j];
+          return distribution_.computePDF(support).asPoint().dot(probabilities);
+        }
+    } // use PDF
+  return DistributionImplementation::computePDF(point);
+}
+
 /* Get the CDF of the MarginalDistribution */
 Scalar MarginalDistribution::computeCDF(const Point & point) const
 {
@@ -270,7 +304,7 @@ CorrelationMatrix MarginalDistribution::getSpearmanCorrelation() const
   return spearmanCorrelation;
 }
 
-/* Get the Spearman correlation of the distribution */
+/* Get the Kendall tau of the distribution */
 CorrelationMatrix MarginalDistribution::getKendallTau() const
 {
   const UnsignedInteger dimension = getDimension();
@@ -298,6 +332,7 @@ Distribution MarginalDistribution::getMarginal(const UnsignedInteger i) const
 /* Get the distribution of the marginal distribution corresponding to indices dimensions */
 Distribution MarginalDistribution::getMarginal(const Indices & indices) const
 {
+", dimension=" << getDimension());
   const UnsignedInteger dimension = getDimension();
   if (!indices.check(dimension)) throw InvalidArgumentException(HERE) << "The indices of a marginal distribution must be in the range [0, dim-1] and must be different";
   if (dimension == 1) return clone();
@@ -381,6 +416,37 @@ Bool MarginalDistribution::isIntegral() const
   return distribution_.isIntegral();
 }
 
+/* UsePDF accessor */
+void MarginalDistribution::setUsePDF(const Bool usePDF)
+{
+  usePDF_ = usePDF;
+}
+
+Bool MarginalDistribution::getUsePDF() const
+{
+  return usePDF_;
+}
+
+/* Get the support of a discrete distribution that intersect a given interval */
+Sample MarginalDistribution::getSupport(const Interval & interval) const
+{
+  if (interval.getDimension() != getDimension()) throw InvalidArgumentException(HERE) << "Error: the given interval has a dimension that does not match the distribution dimension.";
+  const Interval extendedInterval(expandPoint(interval.getLowerBound(), false), expandPoint(interval.getUpperBound()));
+  return distribution_.getSupport(extendedInterval).getMarginal(indices_);
+}
+
+/* Get the support on the whole range */
+Sample MarginalDistribution::getSupport() const
+{
+  return distribution_.getSupport().getMarginal(indices_);
+}
+
+/* Get the discrete probability levels */
+Point MarginalDistribution::getProbabilities() const
+{
+  return distribution_.getProbabilities();
+}
+
 /* Method to expand a given point in the marginal space to a point in the underlying distribution space */
 Point MarginalDistribution::expandPoint(const Point & point,
                                         const Bool upper) const

lib/src/Uncertainty/Distribution/openturns/MarginalDistribution.hxx

@@ -67,7 +67,20 @@ public:
   void setIndices(const Indices & indices);
   Indices getIndices() const;
 
-private:
+  /** UsePDF accessor */
+  void setUsePDF(const Bool usePDF);
+  Bool getUsePDF() const;
+
+  /** Get the support of a discrete distribution that intersect a given interval */
+  using DistributionImplementation::getSupport;
+  Sample getSupport() const override;
+  Sample getSupport(const Interval & interval) const override;
+
+  /** Get the discrete probability levels */
+  using DistributionImplementation::getProbabilities;
+  Point getProbabilities() const override;
+
+ private:
   /** Set the distribution and the indices in one shot */
   void setDistributionAndIndices(const Distribution & distribution,
                                  const Indices & indices);
@@ -82,6 +95,10 @@ public:
   Point getRealization() const override;
   Sample getSample(const UnsignedInteger size) const override;
 
+  /** Get the PDF of the MarginalDistribution */
+  using DistributionImplementation::computePDF;
+  Scalar computePDF(const Point & point) const override;
+
   /** Get the CDF of the MarginalDistribution */
   using DistributionImplementation::computeCDF;
   Scalar computeCDF(const Point & point) const override;
@@ -177,6 +194,10 @@ private:
   /** The upper bound of the underlying distribution */
   Point upperBound_;
 
+  /** Flag to tell if the computations have to be done using the CDF or the PDF
+      of the underlying distribution */
+  Bool usePDF_ = ResourceMap::GetAsBool("MarginalDistribution-UsePDF");
+
 }; /* class MarginalDistribution */