From 5c107ab245c8bfac5446ebd764b688abd8c0a73f Mon Sep 17 00:00:00 2001 From: regislebrun Date: Wed, 9 Oct 2024 22:49:40 +0200 Subject: [PATCH] Improved JointByConditioningDistribution Now, the compute[Sequential]Conditional[PDF|CDF[Quantile] methods have a dedicated implementation, resulting into a huge performance boost. --- .../JointByConditioningDistribution.cxx | 121 ++++++++++++++++++ .../JointByConditioningDistribution.hxx | 15 +++ ...JointByConditioningDistribution_std.expout | 6 + .../t_JointByConditioningDistribution_std.py | 20 +++ 4 files changed, 162 insertions(+) diff --git a/lib/src/Uncertainty/Distribution/JointByConditioningDistribution.cxx b/lib/src/Uncertainty/Distribution/JointByConditioningDistribution.cxx index f72a559df9..3e62260708 100644 --- a/lib/src/Uncertainty/Distribution/JointByConditioningDistribution.cxx +++ b/lib/src/Uncertainty/Distribution/JointByConditioningDistribution.cxx @@ -503,6 +503,127 @@ void JointByConditioningDistribution::computeCovariance() const isAlreadyComputedCovariance_ = true; } +/* Compute the PDF of Xi | X1, ..., Xi-1. x = Xi, y = (X1,...,Xi-1) */ +Scalar JointByConditioningDistribution::computeConditionalPDF(const Scalar x, + const Point & y) const +{ + const UnsignedInteger conditioningDimension = y.getDimension(); + if (conditioningDimension >= getDimension()) throw InvalidArgumentException(HERE) << "Error: cannot compute a conditional PDF with a conditioning point of dimension greater or equal to the distribution dimension."; + // Special case for a conditioning only in the conditioning part + const UnsignedInteger conditioningDistributionDimension = conditioningDistribution_.getDimension(); + if ((conditioningDimension < conditioningDistributionDimension)) + return conditioningDistribution_.computeConditionalPDF(x, y); + // The conditioning part is fully conditioned, let's evaluate the link function + Point fixedConditioningPart(conditioningDistributionDimension); + std::copy(y.begin(), y.begin() + conditioningDistributionDimension, fixedConditioningPart.begin()); + Distribution conditioned(conditionedDistribution_); + conditioned.setParameter(linkFunction_(fixedConditioningPart)); + Point fixedConditionedPart(conditioningDistributionDimension - conditioningDistributionDimension); + if (fixedConditionedPart.getDimension() > 0) + std::copy(y.begin() + conditioningDistributionDimension, y.end(), fixedConditionedPart.begin()); + return conditioned.computeConditionalPDF(x, fixedConditionedPart); +} + +/* Compute the PDF of Xi | X1, ..., Xi-1. x = Xi, y = (X1,...,Xi-1) */ +Point JointByConditioningDistribution::computeSequentialConditionalPDF(const Point & x) const +{ + if (x.getDimension() != dimension_) throw InvalidArgumentException(HERE) << "Error: expected a point of dimension=" << dimension_ << ", got dimension=" << x.getDimension(); + const UnsignedInteger conditioningDistributionDimension = conditioningDistribution_.getDimension(); + Point result(dimension_); + Point conditioningArgument(conditioningDistributionDimension); + std::copy(x.begin(), x.begin() + conditioningDistributionDimension, conditioningArgument.begin()); + const Point conditioningConditionalPDF(conditioningDistribution_.computeSequentialConditionalPDF(conditioningArgument)); + std::copy(conditioningConditionalPDF.begin(), conditioningConditionalPDF.begin() + conditioningDistributionDimension, result.begin()); + Distribution conditioned(conditionedDistribution_); + conditioned.setParameter(linkFunction_(conditioningArgument)); + Point conditionedArgument(conditioned.getDimension()); + std::copy(x.begin() + conditioningDistributionDimension, x.end(), conditionedArgument.begin()); + const Point conditionedConditionalPDF(conditioned.computeSequentialConditionalPDF(conditionedArgument)); + std::copy(conditionedConditionalPDF.begin(), conditionedConditionalPDF.end(), result.begin() + conditioningDistributionDimension); + return result; +} + +/* Compute the CDF of Xi | X1, ..., Xi-1. x = Xi, y = (X1,...,Xi-1) */ +Scalar JointByConditioningDistribution::computeConditionalCDF(const Scalar x, + const Point & y) const +{ + const UnsignedInteger conditioningDimension = y.getDimension(); + if (conditioningDimension >= getDimension()) throw InvalidArgumentException(HERE) << "Error: cannot compute a conditional CDF with a conditioning point of dimension greater or equal to the distribution dimension."; + // Special case for a conditioning only in the conditioning part + const UnsignedInteger conditioningDistributionDimension = conditioningDistribution_.getDimension(); + if ((conditioningDimension < conditioningDistributionDimension)) + return conditioningDistribution_.computeConditionalCDF(x, y); + // The conditioning part is fully conditioned, let's evaluate the link function + Point fixedConditioningPart(conditioningDistributionDimension); + std::copy(y.begin(), y.begin() + conditioningDistributionDimension, fixedConditioningPart.begin()); + Distribution conditioned(conditionedDistribution_); + conditioned.setParameter(linkFunction_(fixedConditioningPart)); + Point fixedConditionedPart(conditioningDimension - conditioningDistributionDimension); + if (fixedConditionedPart.getDimension() > 0) + std::copy(y.begin() + conditioningDistributionDimension, y.end(), fixedConditionedPart.begin()); + return conditioned.computeConditionalCDF(x, fixedConditionedPart); +} + +/* Compute the CDF of Xi | X1, ..., Xi-1. x = Xi, y = (X1,...,Xi-1) */ +Point JointByConditioningDistribution::computeSequentialConditionalCDF(const Point & x) const +{ + if (x.getDimension() != dimension_) throw InvalidArgumentException(HERE) << "Error: expected a point of dimension=" << dimension_ << ", got dimension=" << x.getDimension(); + const UnsignedInteger conditioningDistributionDimension = conditioningDistribution_.getDimension(); + Point result(dimension_); + Point conditioningArgument(conditioningDistributionDimension); + std::copy(x.begin(), x.begin() + conditioningDistributionDimension, conditioningArgument.begin()); + const Point conditioningConditionalCDF(conditioningDistribution_.computeSequentialConditionalCDF(conditioningArgument)); + std::copy(conditioningConditionalCDF.begin(), conditioningConditionalCDF.begin() + conditioningDistributionDimension, result.begin()); + Distribution conditioned(conditionedDistribution_); + conditioned.setParameter(linkFunction_(conditioningArgument)); + Point conditionedArgument(conditioned.getDimension()); + std::copy(x.begin() + conditioningDistributionDimension, x.end(), conditionedArgument.begin()); + const Point conditionedConditionalCDF(conditioned.computeSequentialConditionalCDF(conditionedArgument)); + std::copy(conditionedConditionalCDF.begin(), conditionedConditionalCDF.end(), result.begin() + conditioningDistributionDimension); + return result; +} + +/* Compute the quantile of Xi | X1, ..., Xi-1, i.e. x such that CDF(x|y) = q with x = Xi, y = (X1,...,Xi-1) */ +Scalar JointByConditioningDistribution::computeConditionalQuantile(const Scalar q, + const Point & y) const +{ + const UnsignedInteger conditioningDimension = y.getDimension(); + if (conditioningDimension >= getDimension()) throw InvalidArgumentException(HERE) << "Error: cannot compute a conditional quantile with a conditioning point of dimension greater or equal to the distribution dimension."; + if ((q < 0.0) || (q > 1.0)) throw InvalidArgumentException(HERE) << "Error: cannot compute a conditional quantile for a probability level outside of [0, 1]"; + // Special case for a conditioning only in the conditioning part + const UnsignedInteger conditioningDistributionDimension = conditioningDistribution_.getDimension(); + if ((conditioningDimension < conditioningDistributionDimension)) + return conditioningDistribution_.computeConditionalQuantile(q, y); + // The conditioning part is fully conditioned, let's evaluate the link function + Point fixedConditioningPart(conditioningDistributionDimension); + std::copy(y.begin(), y.begin() + conditioningDistributionDimension, fixedConditioningPart.begin()); + Distribution conditioned(conditionedDistribution_); + conditioned.setParameter(linkFunction_(fixedConditioningPart)); + Point fixedConditionedPart(conditioningDimension - conditioningDistributionDimension); + if (fixedConditionedPart.getDimension() > 0) + std::copy(y.begin() + conditioningDistributionDimension, y.end(), fixedConditionedPart.begin()); + return conditioned.computeConditionalQuantile(q, fixedConditionedPart); +} + +/* Compute the quantile of Xi | X1, ..., Xi-1. x = Xi, y = (X1,...,Xi-1) */ +Point JointByConditioningDistribution::computeSequentialConditionalQuantile(const Point & q) const +{ + if (q.getDimension() != dimension_) throw InvalidArgumentException(HERE) << "Error: expected a point of dimension=" << dimension_ << ", got dimension=" << q.getDimension(); + const UnsignedInteger conditioningDistributionDimension = conditioningDistribution_.getDimension(); + Point result(dimension_); + Point conditioningArgument(conditioningDistributionDimension); + std::copy(q.begin(), q.begin() + conditioningDistributionDimension, conditioningArgument.begin()); + const Point conditioningConditionalQuantile(conditioningDistribution_.computeSequentialConditionalQuantile(conditioningArgument)); + std::copy(conditioningConditionalQuantile.begin(), conditioningConditionalQuantile.begin() + conditioningDistributionDimension, result.begin()); + Distribution conditioned(conditionedDistribution_); + conditioned.setParameter(linkFunction_(conditioningConditionalQuantile)); + Point conditionedArgument(conditioned.getDimension()); + std::copy(q.begin() + conditioningDistributionDimension, q.end(), conditionedArgument.begin()); + const Point conditionedConditionalQuantile(conditioned.computeSequentialConditionalQuantile(conditionedArgument)); + std::copy(conditionedConditionalQuantile.begin(), conditionedConditionalQuantile.end(), result.begin() + conditioningDistributionDimension); + return result; +} + /* Method save() stores the object through the StorageManager */ void JointByConditioningDistribution::save(Advocate & adv) const { diff --git a/lib/src/Uncertainty/Distribution/openturns/JointByConditioningDistribution.hxx b/lib/src/Uncertainty/Distribution/openturns/JointByConditioningDistribution.hxx index 18b1f66790..30eed91089 100644 --- a/lib/src/Uncertainty/Distribution/openturns/JointByConditioningDistribution.hxx +++ b/lib/src/Uncertainty/Distribution/openturns/JointByConditioningDistribution.hxx @@ -102,6 +102,21 @@ public: /** Parameters description accessor */ Description getParameterDescription() const override; + /** Compute the PDF of Xi | X1, ..., Xi-1. x = Xi, y = (X1,...,Xi-1) */ + using DistributionImplementation::computeConditionalPDF; + Scalar computeConditionalPDF(const Scalar x, const Point & y) const override; + Point computeSequentialConditionalPDF(const Point & x) const override; + + /** Compute the CDF of Xi | X1, ..., Xi-1. x = Xi, y = (X1,...,Xi-1) */ + using DistributionImplementation::computeConditionalCDF; + Scalar computeConditionalCDF(const Scalar x, const Point & y) const override; + Point computeSequentialConditionalCDF(const Point & x) const override; + + /** Compute the quantile of Xi | X1, ..., Xi-1, i.e. x such that CDF(x|y) = q with x = Xi, y = (X1,...,Xi-1) */ + using DistributionImplementation::computeConditionalQuantile; + Scalar computeConditionalQuantile(const Scalar q, const Point & y) const override; + Point computeSequentialConditionalQuantile(const Point & q) const override; + /** Method save() stores the object through the StorageManager */ void save(Advocate & adv) const override; diff --git a/python/test/t_JointByConditioningDistribution_std.expout b/python/test/t_JointByConditioningDistribution_std.expout index 9bc990119c..e4e13ece4d 100644 --- a/python/test/t_JointByConditioningDistribution_std.expout +++ b/python/test/t_JointByConditioningDistribution_std.expout @@ -23,3 +23,9 @@ anotherSample covariance= [[ 0.0830803 0.000313403 0.0767514 ] Point point= [0.5, 1.5, 1.0] pdf=0.251589 cdf=0.181637 Quantile= [0.982796,1.9828,3.87457] CDF(quantile)= 0.95 +conditional PDF=1.09340e-01 +conditional CDF=9.08789e-01 +conditional quantile=2.50000e+00 +sequential conditional PDF= [1,1,0.10934] +sequential conditional CDF( [0.5, 1.5, 2.5] )= [0.5,0.5,0.908789] +sequential conditional quantile( [0.5,0.5,0.908789] )= [0.5,1.5,2.5] diff --git a/python/test/t_JointByConditioningDistribution_std.py b/python/test/t_JointByConditioningDistribution_std.py index b143c7fe93..99d968e88d 100755 --- a/python/test/t_JointByConditioningDistribution_std.py +++ b/python/test/t_JointByConditioningDistribution_std.py @@ -56,6 +56,26 @@ print("Quantile=", quantile) print("CDF(quantile)= %.12g" % distribution.computeCDF(quantile)) +x = 2.5 +y = [0.5, 1.5] +print("conditional PDF=%.5e" % distribution.computeConditionalPDF(x, y)) +condCDF = distribution.computeConditionalCDF(x, y) +print("conditional CDF=%.5e" % condCDF) +q = condCDF +print("conditional quantile=%.5e" % distribution.computeConditionalQuantile(q, y)) +pt = [i + 0.5 for i in range(dim)] +print( + "sequential conditional PDF=", distribution.computeSequentialConditionalPDF(pt) +) +resCDF = distribution.computeSequentialConditionalCDF(pt) +print("sequential conditional CDF(", pt, ")=", resCDF) +print( + "sequential conditional quantile(", + resCDF, + ")=", + distribution.computeSequentialConditionalQuantile(resCDF), +) + ot.Log.Show(ot.Log.TRACE) checker = ott.DistributionChecker(distribution) checker.skipMoments() # slow