This is an automated email from the ASF dual-hosted git repository. aherbert pushed a commit to branch master in repository https://gitbox.apache.org/repos/asf/commons-statistics.git
commit 57315f2c7083ead3e18be7980d06f3ee98297323 Author: Alex Herbert <[email protected]> AuthorDate: Fri Oct 22 14:16:16 2021 +0100 STATISTICS-37: Update Normal distribution with high precision erfc Commons Numbers 1.1 improved the implementation of erfc and its inverse. Update the truncated normal distribution and unit tests to verify the increased precision. Added better use of the probability in a range for the truncated normal. Add notes about possible cancellation is computation of the moments for the truncated normal. --- .../distribution/TruncatedNormalDistribution.java | 143 +++++++++++++-------- .../TruncatedNormalDistributionTest.java | 42 ++++++ .../distribution/test.truncatednormal.1.properties | 11 +- .../distribution/test.truncatednormal.7.properties | 6 +- 4 files changed, 146 insertions(+), 56 deletions(-) diff --git a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/TruncatedNormalDistribution.java b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/TruncatedNormalDistribution.java index 8cb94b0..97a781b 100644 --- a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/TruncatedNormalDistribution.java +++ b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/TruncatedNormalDistribution.java @@ -20,7 +20,7 @@ package org.apache.commons.statistics.distribution; /** * Implementation of the truncated normal distribution. * - * @see <a href="http://en.wikipedia.org/wiki/Truncated_normal_distribution";> + * @see <a href="https://en.wikipedia.org/wiki/Truncated_normal_distribution";> * Truncated normal distribution (Wikipedia)</a> */ public final class TruncatedNormalDistribution extends AbstractContinuousDistribution { @@ -28,10 +28,8 @@ public final class TruncatedNormalDistribution extends AbstractContinuousDistrib * This is immutable and thread-safe and can be used across instances. */ private static final NormalDistribution STANDARD_NORMAL = NormalDistribution.of(0, 1); - /** Mean of parent normal distribution. */ - private final double parentMean; - /** Standard deviation of parent normal distribution. */ - private final double parentSd; + /** Parent normal distribution. */ + private final NormalDistribution parentNormal; /** Mean of this distribution. */ private final double mean; /** Variance of this distribution. */ @@ -41,17 +39,14 @@ public final class TruncatedNormalDistribution extends AbstractContinuousDistrib /** Upper bound of this distribution. */ private final double upper; - /** Stored value of @{code standardNormal.cumulativeProbability((lower - mean) / sd)} for faster computations. */ - private final double cdfAlpha; - /** - * Stored value of @{code standardNormal.cumulativeProbability((upper - mean) / sd) - cdfAlpha} - * for faster computations. - */ + /** Stored value of {@code parentNormal.probability(lower, upper)}. This is used to + * normalise the probability computations. */ private final double cdfDelta; - /** parentSd * cdfDelta. */ - private final double parentSdByCdfDelta; - /** log(parentSd * cdfDelta). */ - private final double logParentSdByCdfDelta; + /** log(cdfDelta). */ + private final double logCdfDelta; + /** Stored value of {@code parentNormal.cumulativeProbability(lower)}. Used to map + * a probability into the range of the parent normal distribution. */ + private final double cdfAlpha; /** * @param mean Mean for the parent distribution. @@ -63,48 +58,66 @@ public final class TruncatedNormalDistribution extends AbstractContinuousDistrib this.lower = lower; this.upper = upper; - parentMean = mean; - parentSd = sd; + // Use an instance for the parent normal distribution to maximise accuracy + // in range computations using the error function + parentNormal = NormalDistribution.of(mean, sd); - final double alpha = (lower - mean) / sd; - final double beta = (upper - mean) / sd; - - final double cdfBeta = STANDARD_NORMAL.cumulativeProbability(beta); - cdfAlpha = STANDARD_NORMAL.cumulativeProbability(alpha); - cdfDelta = cdfBeta - cdfAlpha; - - parentSdByCdfDelta = parentSd * cdfDelta; - logParentSdByCdfDelta = Math.log(parentSdByCdfDelta); + cdfDelta = parentNormal.probability(lower, upper); + logCdfDelta = Math.log(cdfDelta); + // Used to map the inverseCumulativeProbability + cdfAlpha = parentNormal.cumulativeProbability(lower); // Calculation of variance and mean. + // + // Use the equations provided on Wikipedia: + // https://en.wikipedia.org/wiki/Truncated_normal_distribution#Moments + + final double alpha = (lower - mean) / sd; + final double beta = (upper - mean) / sd; final double pdfAlpha = STANDARD_NORMAL.density(alpha); final double pdfBeta = STANDARD_NORMAL.density(beta); - final double pdfCdfDelta = (pdfAlpha - pdfBeta) / cdfDelta; - final double alphaBetaDelta = (alpha * pdfAlpha - beta * pdfBeta) / cdfDelta; - if (lower == Double.NEGATIVE_INFINITY) { - if (upper == Double.POSITIVE_INFINITY) { + // lower or upper may be infinite or the density is zero. + + double mu; + double var; + + if (lower == Double.NEGATIVE_INFINITY || pdfAlpha == 0) { + if (upper == Double.POSITIVE_INFINITY || pdfBeta == 0) { // No truncation - this.mean = mean; - variance = sd * sd; + mu = mean; + var = sd * sd; } else { - // One-sided lower tail truncation - final double betaRatio = pdfBeta / cdfBeta; - this.mean = mean - sd * betaRatio; - variance = sd * sd * (1 - beta * betaRatio - betaRatio * betaRatio); + // One sided truncation (of upper tail) + final double betaRatio = pdfBeta / cdfDelta; + mu = mean - sd * betaRatio; + var = sd * sd * (1 - beta * betaRatio - betaRatio * betaRatio); } } else { - if (upper == Double.POSITIVE_INFINITY) { - // One-sided upper tail truncation + if (upper == Double.POSITIVE_INFINITY || pdfBeta == 0) { + // One sided truncation (of lower tail) final double alphaRatio = pdfAlpha / cdfDelta; - this.mean = mean + sd * alphaRatio; - variance = sd * sd * (1 + alpha * alphaRatio - alphaRatio * alphaRatio); + mu = mean + sd * alphaRatio; + var = sd * sd * (1 + alpha * alphaRatio - alphaRatio * alphaRatio); } else { // Two-sided truncation - this.mean = mean + pdfCdfDelta * parentSd; - variance = sd * sd * (1 + alphaBetaDelta - pdfCdfDelta * pdfCdfDelta); + // Note: + // This computation is numerically unstable and requires improvement. + + // Do not use z = cdfDelta which can create cancellation. + final double cdfBeta = parentNormal.cumulativeProbability(upper); + final double z = cdfBeta - cdfAlpha; + final double pdfCdfDelta = (pdfAlpha - pdfBeta) / z; + final double alphaBetaDelta = (alpha * pdfAlpha - beta * pdfBeta) / z; + mu = mean + pdfCdfDelta * sd; + var = sd * sd * (1 + alphaBetaDelta - pdfCdfDelta * pdfCdfDelta); } } + + // The mean should be clipped to the range [lower, upper]. + // The variance should be less than the variance of the parent normal distribution. + this.mean = clipToRange(mu); + variance = Math.min(var, sd * sd); } /** @@ -137,7 +150,13 @@ public final class TruncatedNormalDistribution extends AbstractContinuousDistrib if (x < lower || x > upper) { return 0; } - return STANDARD_NORMAL.density((x - parentMean) / parentSd) / parentSdByCdfDelta; + return parentNormal.density(x) / cdfDelta; + } + + /** {@inheritDoc} */ + @Override + public double probability(double x0, double x1) { + return parentNormal.probability(clipToRange(x0), clipToRange(x1)) / cdfDelta; } /** {@inheritDoc} */ @@ -146,7 +165,7 @@ public final class TruncatedNormalDistribution extends AbstractContinuousDistrib if (x < lower || x > upper) { return Double.NEGATIVE_INFINITY; } - return STANDARD_NORMAL.logDensity((x - parentMean) / parentSd) - logParentSdByCdfDelta; + return parentNormal.logDensity(x) - logCdfDelta; } /** {@inheritDoc} */ @@ -157,7 +176,18 @@ public final class TruncatedNormalDistribution extends AbstractContinuousDistrib } else if (x >= upper) { return 1; } - return (STANDARD_NORMAL.cumulativeProbability((x - parentMean) / parentSd) - cdfAlpha) / cdfDelta; + return parentNormal.probability(lower, x) / cdfDelta; + } + + /** {@inheritDoc} */ + @Override + public double survivalProbability(double x) { + if (x <= lower) { + return 1; + } else if (x >= upper) { + return 0; + } + return parentNormal.probability(x, upper) / cdfDelta; } /** {@inheritDoc} */ @@ -170,12 +200,9 @@ public final class TruncatedNormalDistribution extends AbstractContinuousDistrib } else if (p == 1) { return upper; } - final double x = STANDARD_NORMAL.inverseCumulativeProbability(cdfAlpha + p * cdfDelta) * parentSd + parentMean; - // Clip to support to handle floating-point error at the support bound - if (x <= lower) { - return lower; - } - return x < upper ? x : upper; + // Linearly map p to the range [lower, upper] + final double x = parentNormal.inverseCumulativeProbability(cdfAlpha + p * cdfDelta); + return clipToRange(x); } /** @@ -223,4 +250,18 @@ public final class TruncatedNormalDistribution extends AbstractContinuousDistrib public boolean isSupportConnected() { return true; } + + /** + * Clip to the value to the range [lower, upper]. + * This is used to handle floating-point error at the support bound. + * + * @param x the x + * @return x clipped to the range + */ + private double clipToRange(double x) { + if (x <= lower) { + return lower; + } + return x < upper ? x : upper; + } } diff --git a/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/TruncatedNormalDistributionTest.java b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/TruncatedNormalDistributionTest.java index f4c5fe1..160a68c 100644 --- a/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/TruncatedNormalDistributionTest.java +++ b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/TruncatedNormalDistributionTest.java @@ -17,6 +17,10 @@ package org.apache.commons.statistics.distribution; +import org.junit.jupiter.api.Assertions; +import org.junit.jupiter.params.ParameterizedTest; +import org.junit.jupiter.params.provider.CsvSource; + /** * Test class for {@link TruncatedNormalDistribution}. * Extends {@link BaseContinuousDistributionTest}. See javadoc of that class for details. @@ -47,4 +51,42 @@ class TruncatedNormalDistributionTest extends BaseContinuousDistributionTest { // The constructor arguments do not match the mean and SD of the truncated distribution. return new String[] {null, null, "SupportLowerBound", "SupportUpperBound"}; } + + /** + * Hit the edge cases where the lower and upper bound are not infinite but the + * CDF of the parent distribution is either 0 or 1. This is effectively no truncation. + * Big finite bounds should be handled as if infinite when computing the moments. + * + * @param mean Mean for the parent distribution. + * @param sd Standard deviation for the parent distribution. + * @param lower Lower bound (inclusive) of the distribution, can be {@link Double#NEGATIVE_INFINITY}. + * @param upper Upper bound (inclusive) of the distribution, can be {@link Double#POSITIVE_INFINITY}. + */ + @ParameterizedTest + @CsvSource({ + "0.0, 1.0, -4, 6", + "1.0, 2.0, -4, 6", + "3.45, 6.78, -8, 10", + }) + void testEffectivelyNoTruncation(double mean, double sd, double lower, double upper) { + double inf = Double.POSITIVE_INFINITY; + double max = Double.MAX_VALUE; + TruncatedNormalDistribution dist1; + TruncatedNormalDistribution dist2; + // truncation of upper tail + dist1 = TruncatedNormalDistribution.of(mean, sd, -inf, upper); + dist2 = TruncatedNormalDistribution.of(mean, sd, -max, upper); + Assertions.assertEquals(dist1.getMean(), dist2.getMean(), "Mean"); + Assertions.assertEquals(dist1.getVariance(), dist2.getVariance(), "Variance"); + // truncation of lower tail + dist1 = TruncatedNormalDistribution.of(mean, sd, lower, inf); + dist2 = TruncatedNormalDistribution.of(mean, sd, lower, max); + Assertions.assertEquals(dist1.getMean(), dist2.getMean(), "Mean"); + Assertions.assertEquals(dist1.getVariance(), dist2.getVariance(), "Variance"); + // no truncation + dist1 = TruncatedNormalDistribution.of(mean, sd, -inf, inf); + dist2 = TruncatedNormalDistribution.of(mean, sd, -max, max); + Assertions.assertEquals(dist1.getMean(), dist2.getMean(), "Mean"); + Assertions.assertEquals(dist1.getVariance(), dist2.getVariance(), "Variance"); + } } diff --git a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.1.properties b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.1.properties index 9b5fe50..620b8b5 100644 --- a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.1.properties +++ b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.1.properties @@ -14,14 +14,19 @@ # limitations under the License. parameters = 1.9, 1.3, -1.1, 3.4 -mean = 1.63375792365723 -variance = 1.03158703914439 +# Limited by survival function values. +# SciPy computes (sf(x) - sf(b)) / (sf(a) - sf(b). +# TruncatedNormal computes using the range: probability(x, b) / probability(a, b) +tolerance.relative = 2e-12 +mean = 1.6337579236572282 +variance = 1.0315870391443911 lower = -1.1 upper = 3.4 # Computed using Python with SciPy v1.6.0. # mean, std, clip_a, clip_b = 1.9, 1.3, -1.1, 3.4 # a, b = (clip_a - mean) / std, (clip_b - mean) / std -# truncnorm.var(a, b, loc=mean, scale=std) +# t = truncnorm(a, b, loc=mean, scale=std) +# t.mean(), t.var(), t.cdf(0.5), etc. cdf.points = \ -1.1, -1.09597275767544, -1.0609616183922, -0.79283350106842,\ -0.505331829887808, -0.192170173599874, 0.21173317261645,\ diff --git a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.7.properties b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.7.properties index dc85eae..5084cb1 100644 --- a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.7.properties +++ b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.7.properties @@ -15,8 +15,10 @@ # Test a narrow truncation range. parameters = 7.1, 9.9, 7.0999999, 7.1000001 -# Limited by the mean computed by scipy -# TODO: This test case should be investigated to verify the low tolerance required. +# Limited by the computation of mean and variance +# TODO: The mean/variance computation uses direct formulas and is numerically +# unstable for cases such as this (small range around the mean) +# or where the [lower, upper] interval is far from the mean. tolerance.relative = 1e-7 mean = 7.1 variance = 1.13584123966337e-07
