Skip to content

Commit

Permalink
Improved MarginalDistribution
Browse files Browse the repository at this point in the history
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.
  • Loading branch information
regislebrun committed Apr 4, 2024
1 parent f565121 commit 745f6ab
Show file tree
Hide file tree
Showing 4 changed files with 127 additions and 5 deletions.
6 changes: 6 additions & 0 deletions lib/etc/openturns.conf.in
Original file line number Diff line number Diff line change
Expand Up @@ -542,6 +542,12 @@
<LogNormalFactory-MaximumIteration value_int="50" />
<LogNormalFactory-EstimationMethod value_int="0" />

<!-- MarginalDistribution parameters -->
<MarginalDistribution-UsePDF value_bool="true" />
<MarginalDistribution-MaximumError value_float="1.0e-7" />
<MarginalDistribution-MaximumSubIntervals value_int="128" />
<MarginalDistribution-Rule value_str="G15K31" />

<!-- Meixner parameters -->
<MeixnerDistribution-MaximumAbsoluteError value_float="1.0e-12" />
<MeixnerDistribution-MaximumConstraintError value_float="1.0e-12" />
Expand Down
6 changes: 6 additions & 0 deletions lib/src/Base/Common/ResourceMap.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -1169,6 +1169,12 @@ void ResourceMap::loadDefaultConfiguration()
addAsUnsignedInteger("LogNormalFactory-EstimationMethod", 0);
addAsUnsignedInteger("LogNormalFactory-MaximumIteration", 50);

// MarginalDistribution parameters //
addAsBool("MarginalDistribution-UsePDF", true);
addAsScalar("MarginalDistribution-MaximumError", 1.0e-7);
addAsString("MarginalDistribution-Rule", "G15K31");
addAsUnsignedInteger("MarginalDistribution-MaximumSubIntervals", 128);

// Meixner parameters //
addAsScalar("MeixnerDistribution-MaximumAbsoluteError", 1.0e-12);
addAsScalar("MeixnerDistribution-MaximumConstraintError", 1.0e-12);
Expand Down
83 changes: 79 additions & 4 deletions lib/src/Uncertainty/Distribution/MarginalDistribution.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@

#include "openturns/MarginalDistribution.hxx"
#include "openturns/Uniform.hxx"
#include "openturns/ParametricFunction.hxx"
#include "openturns/PersistentObjectFactory.hxx"

BEGIN_NAMESPACE_OPENTURNS
Expand Down Expand Up @@ -141,10 +142,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
Expand Down Expand Up @@ -192,6 +193,36 @@ 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())));
return integrationAlgorithm_.integrate(kernel, marginalInterval)[0];
}
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
{
Expand Down Expand Up @@ -270,7 +301,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();
Expand Down Expand Up @@ -306,7 +337,9 @@ Distribution MarginalDistribution::getMarginal(const Indices & indices) const
Indices marginalIndices(outputDimension);
for (UnsignedInteger i = 0; i < outputDimension; ++i)
marginalIndices[i] = indices_[indices[i]];
return new MarginalDistribution(distribution_, marginalIndices);
MarginalDistribution marginal(distribution_, marginalIndices);
marginal.setIntegrationAlgorithm(getIntegrationAlgorithm());
return marginal.clone();
}

/* Get the isoprobabilistic transformation */
Expand Down Expand Up @@ -381,6 +414,48 @@ Bool MarginalDistribution::isIntegral() const
return distribution_.isIntegral();
}

/* UsePDF accessor */
void MarginalDistribution::setUsePDF(const Bool usePDF)
{
usePDF_ = usePDF;
}

Bool MarginalDistribution::getUsePDF() const
{
return usePDF_;
}

/* AlgoIntegration accessor */
void MarginalDistribution::setIntegrationAlgorithm(const IntegrationAlgorithm & algo)
{
integrationAlgorithm_ = algo;
}

IntegrationAlgorithm MarginalDistribution::getIntegrationAlgorithm() const
{
return integrationAlgorithm_;
}

/* 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
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,9 @@

#include "openturns/Distribution.hxx"
#include "openturns/DistributionImplementation.hxx"
#include "openturns/IntegrationAlgorithm.hxx"
#include "openturns/IteratedQuadrature.hxx"
#include "openturns/GaussKronrod.hxx"

BEGIN_NAMESPACE_OPENTURNS

Expand Down Expand Up @@ -67,7 +70,24 @@ public:
void setIndices(const Indices & indices);
Indices getIndices() const;

private:
/** UsePDF accessor */
void setUsePDF(const Bool usePDF);
Bool getUsePDF() const;

/** Integration algorithm accessor */
void setIntegrationAlgorithm(const IntegrationAlgorithm & algo);
IntegrationAlgorithm getIntegrationAlgorithm() 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);
Expand All @@ -82,6 +102,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;
Expand Down Expand Up @@ -177,6 +201,17 @@ 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");

/** Integration algorithm used to compute the PDF. By default it is an IteratedQuadrature. */
IntegrationAlgorithm integrationAlgorithm_ = IteratedQuadrature(
GaussKronrod(ResourceMap::GetAsUnsignedInteger("MarginalDistribution-MaximumSubIntervals"),
ResourceMap::GetAsScalar("MarginalDistribution-MaximumError"),
GaussKronrod::GetRuleFromName(ResourceMap::GetAsString("MarginalDistribution-Rule"))
));

}; /* class MarginalDistribution */


Expand Down

0 comments on commit 745f6ab

Please sign in to comment.