You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2021/12/11 08:43:52 UTC
[commons-statistics] 01/03: STATISTICS-46: Update the truncated normal computation of moments
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 62829a20efe3dfd2123ab4227c6c4ec109cd85eb
Author: aherbert <ah...@apache.org>
AuthorDate: Fri Nov 26 18:10:06 2021 +0000
STATISTICS-46: Update the truncated normal computation of moments
The moments cannot be computed when the truncation is very close to zero
as the PDF is flat. Use a uniform distribution in this case.
Extract the logic for the computation of the mean and variance to
methods called on demand.
---
LICENSE | 35 ++
NOTICE | 5 +
commons-statistics-distribution/LICENSE | 35 ++
commons-statistics-distribution/NOTICE | 5 +
.../distribution/TruncatedNormalDistribution.java | 353 ++++++++++++++++-----
.../TruncatedNormalDistributionTest.java | 300 +++++++++++++++++
.../distribution/test.truncatednormal.5.properties | 8 +-
.../distribution/test.truncatednormal.6.properties | 9 +-
.../distribution/test.truncatednormal.7.properties | 35 +-
.../distribution/test.truncatednormal.8.properties | 51 +++
10 files changed, 747 insertions(+), 89 deletions(-)
diff --git a/LICENSE b/LICENSE
index 261eeb9..2256462 100644
--- a/LICENSE
+++ b/LICENSE
@@ -199,3 +199,38 @@
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
+
+================================================================================
+
+APACHE COMMONS STATISTICS SUBCOMPONENTS:
+
+Apache Commons Statistics includes a number of subcomponents with separate
+copyright notices and license terms. Your use of these subcomponents is
+subject to the terms and conditions of the following licenses.
+
+
+For the Commons Statistics Distributions (commons-statistics-distributions-x.x.x.jar)
+component:
+
+ * Part of the distributions package is licensed under the
+ MIT "Expat" License
+
+ Copyright (c) 2017: Jorge Fernandez-de-Cossio-Diaz.
+
+ Permission is hereby granted, free of charge, to any person obtaining a
+ copy of this software and associated documentation files (the "Software"),
+ to deal in the Software without restriction, including without limitation
+ the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ and/or sell copies of the Software, and to permit persons to whom the
+ Software is furnished to do so, subject to the following conditions:
+
+ The above copyright notice and this permission notice shall be included in
+ all copies or substantial portions of the Software.
+
+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ DEALINGS IN THE SOFTWARE.
diff --git a/NOTICE b/NOTICE
index d00085d..49dd9a1 100644
--- a/NOTICE
+++ b/NOTICE
@@ -3,3 +3,8 @@ Copyright 2018-2021 The Apache Software Foundation
This product includes software developed at
The Apache Software Foundation (http://www.apache.org/).
+
+The org.apache.commons.statistics.distribution package contains derivative
+work originating from the Julia implementation of the moments of the
+truncated normal distribution.
+Copyright 2017 Jorge Fernandez-de-Cossio-Diaz.
diff --git a/commons-statistics-distribution/LICENSE b/commons-statistics-distribution/LICENSE
index 261eeb9..0de3dec 100644
--- a/commons-statistics-distribution/LICENSE
+++ b/commons-statistics-distribution/LICENSE
@@ -199,3 +199,38 @@
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
+
+================================================================================
+
+APACHE COMMONS STATISTICS SUBCOMPONENTS:
+
+Apache Commons Statistics includes a number of subcomponents with separate
+copyright notices and license terms. Your use of these subcomponents is
+subject to the terms and conditions of the following licenses.
+
+
+For the Commons Statistics Distributions (commons-statistics-distributions-x.x.x.jar)
+component:
+
+ * Part of the distributions package is licensed under the
+ MIT Expat License
+
+ Copyright (c) 2017: Jorge Fernandez-de-Cossio-Diaz.
+
+ Permission is hereby granted, free of charge, to any person obtaining a
+ copy of this software and associated documentation files (the "Software"),
+ to deal in the Software without restriction, including without limitation
+ the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ and/or sell copies of the Software, and to permit persons to whom the
+ Software is furnished to do so, subject to the following conditions:
+
+ The above copyright notice and this permission notice shall be included in
+ all copies or substantial portions of the Software.
+
+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ DEALINGS IN THE SOFTWARE.
diff --git a/commons-statistics-distribution/NOTICE b/commons-statistics-distribution/NOTICE
index d00085d..49dd9a1 100644
--- a/commons-statistics-distribution/NOTICE
+++ b/commons-statistics-distribution/NOTICE
@@ -3,3 +3,8 @@ Copyright 2018-2021 The Apache Software Foundation
This product includes software developed at
The Apache Software Foundation (http://www.apache.org/).
+
+The org.apache.commons.statistics.distribution package contains derivative
+work originating from the Julia implementation of the moments of the
+truncated normal distribution.
+Copyright 2017 Jorge Fernandez-de-Cossio-Diaz.
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 55e6368..1dec468 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
@@ -17,6 +17,10 @@
package org.apache.commons.statistics.distribution;
+import org.apache.commons.numbers.gamma.Erf;
+import org.apache.commons.numbers.gamma.ErfDifference;
+import org.apache.commons.numbers.gamma.Erfcx;
+
/**
* Implementation of the truncated normal distribution.
*
@@ -24,16 +28,26 @@ package org.apache.commons.statistics.distribution;
* Truncated normal distribution (Wikipedia)</a>
*/
public final class TruncatedNormalDistribution extends AbstractContinuousDistribution {
- /** A standard normal distribution used for calculations.
- * This is immutable and thread-safe and can be used across instances. */
- private static final NormalDistribution STANDARD_NORMAL = NormalDistribution.of(0, 1);
+
+ /** The max allowed value for x where (x*x) will not overflow.
+ * This is a limit on computation of the moments of the truncated normal
+ * as some calculations assume x*x is finite. Value is sqrt(MAX_VALUE). */
+ private static final double MAX_X = 0x1.fffffffffffffp511;
+
+ /** The min allowed probability range of the parent normal distribution.
+ * Set to 0.0. This may be too low for accurate usage. It is a signal that
+ * the truncation is invalid. */
+ private static final double MIN_P = 0.0;
+
+ /** sqrt(2). */
+ private static final double ROOT2 = 1.414213562373095048801688724209698078;
+ /** Normalisation constant 2 / sqrt(2 pi) = sqrt(2 / pi). */
+ private static final double ROOT_2_PI = 0.797884560802865405726436165423365309;
+ /** Normalisation constant sqrt(2 pi) / 2 = sqrt(pi / 2). */
+ private static final double ROOT_PI_2 = 1.253314137315500251207882642405522626;
/** Parent normal distribution. */
private final NormalDistribution parentNormal;
- /** Mean of this distribution. */
- private final double mean;
- /** Variance of this distribution. */
- private final double variance;
/** Lower bound of this distribution. */
private final double lower;
/** Upper bound of this distribution. */
@@ -52,76 +66,21 @@ public final class TruncatedNormalDistribution extends AbstractContinuousDistrib
private final double sfBeta;
/**
- * @param mean Mean for the parent distribution.
- * @param sd Standard deviation for the parent distribution.
+ * @param parent Parent distribution.
+ * @param z Probability of the parent distribution for {@code [lower, upper]}.
* @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}.
*/
- private TruncatedNormalDistribution(double mean, double sd, double lower, double upper) {
+ private TruncatedNormalDistribution(NormalDistribution parent, double z, double lower, double upper) {
+ this.parentNormal = parent;
this.lower = lower;
this.upper = upper;
- // Use an instance for the parent normal distribution to maximise accuracy
- // in range computations using the error function
- parentNormal = NormalDistribution.of(mean, sd);
-
- cdfDelta = parentNormal.probability(lower, upper);
+ cdfDelta = z;
logCdfDelta = Math.log(cdfDelta);
// Used to map the inverse probability.
cdfAlpha = parentNormal.cumulativeProbability(lower);
sfBeta = parentNormal.survivalProbability(upper);
-
- // 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);
-
- // lower or upper may be infinite or the density is zero.
-
- double mu;
- double v;
-
- if (lower == Double.NEGATIVE_INFINITY || pdfAlpha == 0) {
- if (upper == Double.POSITIVE_INFINITY || pdfBeta == 0) {
- // No truncation
- mu = mean;
- v = sd * sd;
- } else {
- // One sided truncation (of upper tail)
- final double betaRatio = pdfBeta / cdfDelta;
- mu = mean - sd * betaRatio;
- v = sd * sd * (1 - beta * betaRatio - betaRatio * betaRatio);
- }
- } else {
- if (upper == Double.POSITIVE_INFINITY || pdfBeta == 0) {
- // One sided truncation (of lower tail)
- final double alphaRatio = pdfAlpha / cdfDelta;
- mu = mean + sd * alphaRatio;
- v = sd * sd * (1 + alpha * alphaRatio - alphaRatio * alphaRatio);
- } else {
- // Two-sided truncation
- // 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;
- v = 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(v, sd * sd);
}
/**
@@ -135,7 +94,8 @@ public final class TruncatedNormalDistribution extends AbstractContinuousDistrib
* @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}.
* @return the distribution
- * @throws IllegalArgumentException if {@code sd <= 0} or if {@code upper <= lower}.
+ * @throws IllegalArgumentException if {@code sd <= 0}; if {@code lower >= upper}; or if
+ * the truncation covers no probability range in the parent distribution.
*/
public static TruncatedNormalDistribution of(double mean, double sd, double lower, double upper) {
if (sd <= 0) {
@@ -144,8 +104,26 @@ public final class TruncatedNormalDistribution extends AbstractContinuousDistrib
if (lower >= upper) {
throw new DistributionException(DistributionException.INVALID_RANGE_LOW_GTE_HIGH, lower, upper);
}
- return new TruncatedNormalDistribution(mean, sd, lower, upper);
+ // Use an instance for the parent normal distribution to maximise accuracy
+ // in range computations using the error function
+ final NormalDistribution parent = NormalDistribution.of(mean, sd);
+
+ // If there is no computable range then raise an exception.
+ final double z = parent.probability(lower, upper);
+ if (z <= MIN_P) {
+ // Map the bounds to a standard normal distribution for the message
+ final double a = (lower - mean) / sd;
+ final double b = (upper - mean) / sd;
+ throw new DistributionException(
+ "Excess truncation of standard normal : CDF(%s, %s) = %s", a, b, z);
+ }
+
+ // Here we have a meaningful truncation. Note that excess truncation may not be optimal.
+ // For example truncation close to zero where the PDF is constant can be approximated
+ // using a uniform distribution.
+
+ return new TruncatedNormalDistribution(parent, z, lower, upper);
}
/** {@inheritDoc} */
@@ -160,6 +138,10 @@ public final class TruncatedNormalDistribution extends AbstractContinuousDistrib
/** {@inheritDoc} */
@Override
public double probability(double x0, double x1) {
+ if (x0 > x1) {
+ throw new DistributionException(DistributionException.INVALID_RANGE_LOW_GT_HIGH,
+ x0, x1);
+ }
return parentNormal.probability(clipToRange(x0), clipToRange(x1)) / cdfDelta;
}
@@ -232,7 +214,11 @@ public final class TruncatedNormalDistribution extends AbstractContinuousDistrib
*/
@Override
public double getMean() {
- return mean;
+ final double u = parentNormal.getMean();
+ final double s = parentNormal.getStandardDeviation();
+ final double a = (lower - u) / s;
+ final double b = (upper - u) / s;
+ return u + moment1(a, b) * s;
}
/**
@@ -243,7 +229,11 @@ public final class TruncatedNormalDistribution extends AbstractContinuousDistrib
*/
@Override
public double getVariance() {
- return variance;
+ final double u = parentNormal.getMean();
+ final double s = parentNormal.getStandardDeviation();
+ final double a = (lower - u) / s;
+ final double b = (upper - u) / s;
+ return variance(a, b) * s * s;
}
/** {@inheritDoc} */
@@ -259,16 +249,233 @@ public final class TruncatedNormalDistribution extends AbstractContinuousDistrib
}
/**
- * Clip to the value to the range [lower, upper].
+ * Clip the value to the range [lower, upper].
* This is used to handle floating-point error at the support bound.
*
- * @param x the x
+ * @param x Value x
* @return x clipped to the range
*/
private double clipToRange(double x) {
+ return clip(x, lower, upper);
+ }
+
+ /**
+ * Clip the value to the range [lower, upper].
+ *
+ * @param x Value x
+ * @param lower Lower bound (inclusive)
+ * @param upper Upper bound (inclusive)
+ * @return x clipped to the range
+ */
+ private static double clip(double x, double lower, double upper) {
if (x <= lower) {
return lower;
}
return x < upper ? x : upper;
}
+
+ // Calculation of variance and mean can suffer from cancellation.
+ //
+ // Use formulas from Jorge Fernandez-de-Cossio-Diaz adapted under the
+ // terms of the MIT "Expat" License (see NOTICE and LICENSE).
+ //
+ // These formulas use the complementary error function
+ // erfcx(z) = erfc(z) * exp(z^2)
+ // This avoids computation of exp terms for the Gaussian PDF and then
+ // dividing by the error functions erf or erfc:
+ // exp(-0.5*x*x) / erfc(x / sqrt(2)) == 1 / erfcx(x / sqrt(2))
+ // At large z the erfcx function is computable but exp(-0.5*z*z) and
+ // erfc(z) are zero. Use of these formulas allows computation of the
+ // mean and variance for the usable range of the truncated distribution
+ // (cdf(a, b) != 0). The variance is not accurate when it approaches
+ // machine epsilon (2^-52) at extremely narrow truncations and the
+ // computation -> 0.
+ //
+ // See: https://github.com/cossio/TruncatedNormal.jl
+
+ /**
+ * Compute the first moment (mean) of the truncated standard normal distribution.
+ *
+ * <p>Assumes {@code a <= b}.
+ *
+ * @param a Lower bound
+ * @param b Upper bound
+ * @return the first moment
+ */
+ static double moment1(double a, double b) {
+ // Assume a <= b
+ if (a == b) {
+ return a;
+ }
+ if (Math.abs(a) > Math.abs(b)) {
+ // Subtract from zero to avoid generating -0.0
+ return 0 - moment1(-b, -a);
+ }
+
+ // Here:
+ // |a| <= |b|
+ // a < b
+ // 0 < b
+
+ if (a <= -MAX_X) {
+ // No truncation
+ return 0;
+ }
+ if (b >= MAX_X) {
+ // One-sided truncation
+ return ROOT_2_PI / Erfcx.value(a / ROOT2);
+ }
+
+ // pdf = exp(-0.5*x*x) / sqrt(2*pi)
+ // cdf = erfc(-x/sqrt(2)) / 2
+ // Compute:
+ // -(pdf(b) - pdf(a)) / cdf(b, a)
+ // Note:
+ // exp(-0.5*b*b) - exp(-0.5*a*a)
+ // Use cancellation of powers:
+ // exp(-0.5*(b*b-a*a)) * exp(-0.5*a*a) - exp(-0.5*a*a)
+ // expm1(-0.5*(b*b-a*a)) * exp(-0.5*a*a)
+
+ // dx = -0.5*(b*b-a*a)
+ final double dx = 0.5 * (b + a) * (b - a);
+ double m;
+ if (a <= 0) {
+ // Opposite signs
+ m = ROOT_2_PI * -Math.expm1(-dx) * Math.exp(-0.5 * a * a) / ErfDifference.value(a / ROOT2, b / ROOT2);
+ } else {
+ final double z = Math.exp(-dx) * Erfcx.value(b / ROOT2) - Erfcx.value(a / ROOT2);
+ if (z == 0) {
+ // Occurs when a and b have large magnitudes and are very close
+ return (a + b) * 0.5;
+ }
+ m = ROOT_2_PI * Math.expm1(-dx) / z;
+ }
+
+ // Clip to the range
+ return clip(m, a, b);
+ }
+
+ /**
+ * Compute the second moment of the truncated standard normal distribution.
+ *
+ * <p>Assumes {@code a <= b}.
+ *
+ * @param a Lower bound
+ * @param b Upper bound
+ * @return the first moment
+ */
+ private static double moment2(double a, double b) {
+ // Assume a < b.
+ // a == b is handled in the variance method
+ if (Math.abs(a) > Math.abs(b)) {
+ return moment2(-b, -a);
+ }
+
+ // Here:
+ // |a| <= |b|
+ // a < b
+ // 0 < b
+
+ if (a <= -MAX_X) {
+ // No truncation
+ return 1;
+ }
+ if (b >= MAX_X) {
+ // One-sided truncation.
+ // For a -> inf : moment2 -> a*a
+ // This occurs when erfcx(z) is approximated by (1/sqrt(pi)) / z and terms
+ // cancel. z > 6.71e7, a > 9.49e7
+ return 1 + ROOT_2_PI * a / Erfcx.value(a / ROOT2);
+ }
+
+ // pdf = exp(-0.5*x*x) / sqrt(2*pi)
+ // cdf = erfc(-x/sqrt(2)) / 2
+ // Compute:
+ // 1 - (b*pdf(b) - a*pdf(a)) / cdf(b, a)
+ // = (cdf(b, a) - b*pdf(b) -a*pdf(a)) / cdf(b, a)
+
+ // Note:
+ // For z -> 0:
+ // sqrt(pi / 2) * erf(z / sqrt(2)) -> z
+ // z * Math.exp(-0.5 * z * z) -> z
+ // Both computations below have cancellation as b -> 0 and the
+ // second moment is not computable as the fraction P/Q
+ // since P < ulp(Q). This always occurs when b < MIN_X
+ // if MIN_X is set at the point where
+ // exp(-0.5 * z * z) / sqrt(2 pi) == 1 / sqrt(2 pi).
+ // This is JDK dependent due to variations in Math.exp.
+ // For b < MIN_X the second moment can be approximated using
+ // a uniform distribution: (b^3 - a^3) / (3b - 3a).
+ // In practice it also occurs when b > MIN_X since any a < MIN_X
+ // is effectively zero for part of the computation. A
+ // threshold to transition to a uniform distribution
+ // approximation is a compromise. Also note it will not
+ // correct computation when (b-a) is small and is far from 0.
+ // Thus the second moment is left to be inaccurate for
+ // small ranges (b-a) and the variance -> 0 when the true
+ // variance is close to or below machine epsilon.
+
+ double m;
+
+ if (a <= 0) {
+ // Opposite signs
+ final double ea = ROOT_PI_2 * Erf.value(a / ROOT2);
+ final double eb = ROOT_PI_2 * Erf.value(b / ROOT2);
+ final double fa = ea - a * Math.exp(-0.5 * a * a);
+ final double fb = eb - b * Math.exp(-0.5 * b * b);
+ // Assume fb >= fa && eb >= ea
+ // If fb <= fa this is a tiny range around 0
+ m = (fb - fa) / (eb - ea);
+ // Clip to the range
+ m = clip(m, 0, 1);
+ } else {
+ final double dx = 0.5 * (b + a) * (b - a);
+ final double ex = Math.exp(-dx);
+ final double ea = ROOT_PI_2 * Erfcx.value(a / ROOT2);
+ final double eb = ROOT_PI_2 * Erfcx.value(b / ROOT2);
+ final double fa = ea + a;
+ final double fb = eb + b;
+ m = (fa - fb * ex) / (ea - eb * ex);
+ // Clip to the range
+ m = clip(m, a * a, b * b);
+ }
+ return m;
+ }
+
+ /**
+ * Compute the variance of the truncated standard normal distribution.
+ *
+ * <p>Assumes {@code a <= b}.
+ *
+ * @param a Lower bound
+ * @param b Upper bound
+ * @return the first moment
+ */
+ static double variance(double a, double b) {
+ if (a == b) {
+ return 0;
+ }
+
+ final double m1 = moment1(a, b);
+ double m2 = moment2(a, b);
+ // variance = m2 - m1*m1
+ // rearrange x^2 - y^2 as (x-y)(x+y)
+ m2 = Math.sqrt(m2);
+ final double var = (m2 - m1) * (m2 + m1);
+
+ // Detect floating-point error.
+ if (var >= 1) {
+ // Note:
+ // Extreme truncations in the tails can compute a variance above 1,
+ // for example if m2 is infinite: m2 - m1*m1 > 1
+ // Detect no truncation as the terms a and b lie far either side of zero;
+ // otherwise return 0 to indicate very small unknown variance.
+ return a < -1 && b > 1 ? 1 : 0;
+ } else if (var <= 0) {
+ // Floating-point error can create negative variance so return 0.
+ return 0;
+ }
+
+ return var;
+ }
}
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 160a68c..2d8290e 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,7 +17,10 @@
package org.apache.commons.statistics.distribution;
+import org.apache.commons.numbers.gamma.Erf;
+import org.apache.commons.numbers.gamma.Erfcx;
import org.junit.jupiter.api.Assertions;
+import org.junit.jupiter.api.Test;
import org.junit.jupiter.params.ParameterizedTest;
import org.junit.jupiter.params.provider.CsvSource;
@@ -42,6 +45,8 @@ class TruncatedNormalDistributionTest extends BaseContinuousDistributionTest {
{0.0, 0.0, -1.0, 1.0},
{0.0, -0.1, -1.0, 1.0},
{0.0, 1.0, 1.0, -1.0},
+ // No usable probability range
+ {0.0, 1.0, 100.0, 101.0},
};
}
@@ -89,4 +94,299 @@ class TruncatedNormalDistributionTest extends BaseContinuousDistributionTest {
Assertions.assertEquals(dist1.getMean(), dist2.getMean(), "Mean");
Assertions.assertEquals(dist1.getVariance(), dist2.getVariance(), "Variance");
}
+
+ /**
+ * Test mean cases adapted from the source implementation for the truncated
+ * normal moments.
+ *
+ * @see <a href="https://github.com/cossio/TruncatedNormal.jl/blob/master/test/tnmom1.jl">
+ * cossio TruncatedNormal moment1 tests</a>
+ */
+ @Test
+ void testMean() {
+ assertMean(Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, 0, 0);
+ assertMean(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, 0);
+ assertMean(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY, 0);
+ assertMean(0, Double.POSITIVE_INFINITY, Math.sqrt(2 / Math.PI), 1e-15);
+ assertMean(Double.NEGATIVE_INFINITY, 0, -Math.sqrt(2 / Math.PI), 1e-15);
+
+ for (int x = -10; x <= 10; x++) {
+ final double expected = Math.sqrt(2 / Math.PI) / Erfcx.value(x / Math.sqrt(2));
+ assertMean(x, Double.POSITIVE_INFINITY, expected, 1e-15);
+ }
+
+ for (int i = -100; i <= 100; i++) {
+ final double x = Math.exp(i);
+ assertMean(-x, x, 0, 0);
+ final double expected = -Math.sqrt(2 / Math.PI) * Math.expm1(-x * x / 2) / Erf.value(x / Math.sqrt(2));
+ assertMean(0, x, expected, 1e-15);
+ }
+
+ assertMean(1e-44, 1e-43, 5.4999999999999999999999999999999999999999e-44, 1e-15);
+
+ assertMean(100, 115, 100.00999800099926070518490239457545847490332879043, 1e-15);
+ assertMean(-1e6, -999000, -999000.00000100100100099899498898098, 1e-15);
+ assertMean(+1e6, Double.POSITIVE_INFINITY, +1.00000000000099999999999800000e6, 1e-15);
+ assertMean(Double.NEGATIVE_INFINITY, -1e6, -1.00000000000099999999999800000e6, 1e-15);
+
+ assertMean(-1e200, 1e200, 0, 1e-15);
+ assertMean(0, +1e200, +0.797884560802865355879892119869, 1e-15);
+ assertMean(-1e200, 0, -0.797884560802865355879892119869, 1e-15);
+
+ assertMean(50, 70, -2, 3, 50.171943499898757645751683644632860837133138152489, 1e-15);
+ assertMean(-100.0, 0.0, 0.0, 2.0986317998643735, -1.6744659119217125058885983754999713622460154892645, 1e-15);
+ assertMean(0.0, 0.9, 0.0, 0.07132755843183151, 0.056911157632522598806524588414964004271754161737065, 1e-15);
+ assertMean(-100.0, 100.0, 0.0, 17.185261847875548, 0, 1e-15);
+ assertMean(-100.0, 0.5, 0.0, 0.47383322897860064, -0.1267981330521791493635176736743283314399, 1e-15);
+ assertMean(-100.0, 100.0, 0.0, 17.185261847875548, 0, 1e-15);
+
+ for (int i = -10; i <= 10; i++) {
+ final double a = Math.exp(i);
+ for (int j = -10; j <= 10; j++) {
+ final double b = Math.exp(j);
+ if (a <= b) {
+ final double mean = TruncatedNormalDistribution.moment1(a, b);
+ Assertions.assertTrue(a <= mean && mean <= b);
+ }
+ }
+ }
+
+ // https://github.com/JuliaStats/Distributions.jl/issues/827, 1e-15);
+ assertMean(0, 1000, 1000000, 1, 999.99999899899899900100501101901899090472046236710608108591983, 6e-14);
+ }
+
+ /**
+ * Test variance cases adapted from the source implementation for the truncated
+ * normal moments.
+ *
+ * @see <a href="https://github.com/cossio/TruncatedNormal.jl/blob/master/test/tnvar.jl">
+ * cossio TruncatedNormal variance tests</a>
+ */
+ @Test
+ void testVariance() {
+ assertVariance(Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, 1, 0);
+ assertVariance(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, 0, 0);
+ assertVariance(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY, 0, 0);
+ assertVariance(0, Double.POSITIVE_INFINITY, 1 - 2 / Math.PI, 1e-15);
+ assertVariance(Double.NEGATIVE_INFINITY, 0, 1 - 2 / Math.PI, 1e-15);
+
+ for (int x = -10; x <= 10; x++) {
+ final double expected = 1 + Math.sqrt(2 / Math.PI) * x / Erfcx.value(x / Math.sqrt(2)) -
+ (2 / Math.PI) / Math.pow(Erfcx.value(x / Math.sqrt(2)), 2);
+ assertVariance(x, Double.POSITIVE_INFINITY, expected, 1e-11);
+ }
+
+ assertVariance(50, 70, 0.0003990431868038995479099272265360593305365, 1e-9);
+
+ assertVariance(50, 70, -2, 3, 0.029373438107168350377591231295634273607812172191712, 1e-11);
+ assertVariance(-100.0, 0.0, 0.0, 2.0986317998643735, 1.6004193412141677189841357987638847137391508803335, 1e-15);
+ assertVariance(0.0, 0.9, 0.0, 0.07132755843183151, 0.0018487407287725028827020557707636415445504260892486, 1e-15);
+ assertVariance(-100.0, 100.0, 0.0, 17.185261847875548, 295.333163899557735486302841237124507431445, 1e-15);
+ assertVariance(-100.0, 0.5, 0.0, 0.47383322897860064, 0.145041095812679283837328561547251019229612, 1e-15);
+ assertVariance(-100.0, 100.0, 0.0, 17.185261847875548, 295.333163899557735486302841237124507431445, 1e-15);
+ assertVariance(-10000, 10000, 0, 1, 1, 1e-15);
+
+ // https://github.com/JuliaStats/Distributions.jl/issues/827
+ Assertions.assertTrue(TruncatedNormalDistribution.variance(999000, 1e6) >= 0);
+ Assertions.assertTrue(TruncatedNormalDistribution.variance(-1000000, 1000 - 1000000) >= 0);
+
+ // These tests are marked as broken in the reference implementation.
+ // They present extreme deviations of the truncation bounds from the mean.
+ //assertVariance(1e6, Double.POSITIVE_INFINITY, 9.99999999994000000000050000000e-13, 1e-15);
+ //assertVariance(999000, 1e6, 1.00200300399898194688784897455e-12, 1e-15);
+ //assertVariance(-1e6, -999000, 1.00200300399898194688784897455e-12, 1e-15);
+ }
+
+ /**
+ * Test cases for computation of the moments. This hits edge cases including truncations
+ * too extreme to have a probability range for the distribution.
+ * The test ensures that the moments are computable for parameterisations
+ * where the bounds fall within +/- 40 standard deviations from the mean.
+ *
+ * <p>Test data generated using a 128-bit implementation of the method using GCC lib quadmath
+ * and Boost C++ Error function routines adapted to compute erfcx. Data verified using
+ * the Julia implementation:
+ * <pre>
+ * import Pkg
+ * Pkg.add(url="https://github.com/cossio/TruncatedNormal.jl")
+ * using TruncatedNormal
+ *
+ * tnmean(1.23, 4.56) # 1.7122093853640246
+ * tnvar(1.23, 4.56) # 0.1739856461219162
+ *
+ * # Using BigFloat does not work on hard cases of the variance
+ * tnvar(BigFloat(1.0), BigFloat(1.0000000000000002))
+ * </pre>
+ *
+ * <p>Computation of the mean is stable. Computation of the variance is not accurate as it
+ * approaches machine epsilon (2^-52). Using Julia's BigFloat support does not allow computation
+ * of the difficult cases listed below for the variance.
+ *
+ * @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}.
+ * @param mean Expected mean
+ * @param variance Expected variance
+ * @param meanRelativeError Relative error tolerance for the mean
+ * @param varianceRelativeError Relative error tolerance for the variance
+ * (if set to negative the variance is allowed to be within 1.5 * epsilon of zero)
+ */
+ @ParameterizedTest
+ @CsvSource({
+ // Equal bounds
+ "1.23, 1.23, 1.23, 0, 0, 0",
+ "1.23, 4.56, 1.7122093853640246, 0.1739856461219162, 1e-15, 5e-15",
+
+ // Effectively no truncation
+ "-55, 60, 0, 1, 0, 0",
+
+ // Long tail
+ "-100, 101, 1.3443134677817230433408433600205167e-2172, 1, 1e-15, 1e-15",
+ "-40, 101, 1.46327025083830317873709720033828097e-348, 1, 1e-15, 1e-15",
+ "-30, 101, 1.47364613487854751904949326604507453e-196, 1, 1e-15, 1e-15",
+ "-20, 101, 5.52094836215976318958273568278700042e-88, 1, 1e-15, 1e-15",
+ "-10, 101, 7.69459862670641934633909221175249367e-23, 0.999999999999999999999230540137329438, 1e-15, 1e-15",
+ "-5, 101, 1.48671994090490571244174411946057083e-06, 0.999992566398085139288753504945569711, 1e-15, 1e-15",
+ "-1, 101, 0.287599970939178361228670127385217202, 0.629686285776605400861244494862843017, 1e-15, 1e-15",
+ "0, 101, 0.797884560802865355879892119868763748, 0.363380227632418656924464946509942526, 1e-15, 1e-15",
+ "1, 101, 1.52513527616098120908909053639057876, 0.199097665570348791553367979096726767, 1e-15, 1e-14",
+ "5, 101, 5.18650396712584211561650896200523673, 0.032696434617112225345315807700917674, 1e-15, 1e-13",
+ "10, 101, 10.0980932339625119628436416537120371, 0.00944537782565626116413681765035684208, 1e-15, 1e-11",
+ "20, 101, 20.0497530685278505422140233087209891, 0.00246326161505216359968528619980015911, 1e-15, 1e-11",
+ "30, 101, 30.033259667433677037071124100012257, 0.00110377151189009100113674138540728116, 1e-15, 1e-10",
+ "40, 101, 40.0249688472072637232448709953697417, 0.000622668378591388773498879400697584317, 1e-15, 2e-9",
+ "100, 101, 100.009998000999260705184902394575471, 9.99400499482634503612772420030347819e-05, 1e-15, 2e-8",
+
+ // One-sided truncation
+ "-5, Infinity, 1.4867199409049057124417441194605712e-06, 0.999992566398085139288753504945569711, 1e-14, 1e-14",
+ "-3, Infinity, 0.00443783904212566379330210431090259846, 0.98666678845825919379095350748267984, 1e-15, 1e-15",
+ "-1, Infinity, 0.287599970939178361228670127385217154, 0.629686285776605400861244494862843306, 1e-15, 1e-15",
+ "0, Infinity, 0.797884560802865355879892119868763748, 0.363380227632418656924464946509942526, 1e-15, 1e-15",
+ "1, Infinity, 1.52513527616098120908909053639057876, 0.199097665570348791553367979096726767, 1e-15, 1e-15",
+ "3, Infinity, 3.28309865493043650692809222681220005, 0.0705591867852681168624020577420568271, 1e-15, 2e-14",
+ "20, Infinity, 20.0497530685278505422140233087209891, 0.00246326161505216359968528619980015911, 1e-15, 1e-11",
+ "100, Infinity, 100.009998000999260705184902394575471, 9.99400499482634503612772420030347819e-05, 1e-15, 4e-8",
+ // The variance method is inaccurate at this extreme
+ "1e4, Infinity, 10000.0000999999980000000999999925986, 9.99999940000005002391967510312099493e-09, 1e-15, 0.8",
+ "1e6, Infinity, 1000000.00000099999999999800000000016, 9.99999999770471649802883928921316157e-13, 1e-15, 1.0",
+ // XXX: The expected variance here is incorrect. It will be small but may be non zero.
+ // The computation will return 0. This hits an edge case in the code that detects when the
+ // variance computation fails.
+ "1e100, Infinity, 1.00000000000000001590289110975991788e+100, 0, 1e-15, -1",
+
+ // XXX: The expected variance here is incorrect. It will be small but may be non zero.
+ // This hits an edge case where the computed variance (infinity) is above 1
+ "1e290, 1e300, 1.00000000000000006172783352786715689e+290, 0, 1e-15, -1",
+
+ // Small ranges.
+ "1, 1.1000000000000001, 1.04912545221799091312759556239135752, 0.000832596851563726615564931035799390151, 1e-15, 2e-12",
+ "5, 5.0999999999999996, 5.04581083165668427678725919870992629, 0.000822546087919772895415146023240560636, 1e-15, 2e-11",
+ "35, 35.100000000000001, 35.025438801080858717764612789648226, 0.000494605845872597846399929727938197022, 1e-15, 2e-9",
+
+ // (b-a) = 1 ULP
+ // XXX: The expected variance here is incorrect.
+ // It is upper limited to the variance of a uniform distribution.
+ // The computation will return 0. This hits an edge case in the code that detects when the
+ // variance computation fails.
+ // Spans p=8.327e-17 of the parent normal distribution
+ "1, 1.0000000000000002, 1.00000000000000011091535982917837267, 0, 1e-15, -1",
+ // Spans p=1.626e-19 of the parent normal distribution
+ "4, 4.0000000000000009, 4.00000000000000044406536771487238653, 0, 1e-15, -1",
+ // Spans p=1.925e-37 of the parent normal distribution
+ "10, 10.000000000000002, 10.0000000000000008883225369216741152, 0, 1e-15, -1",
+
+ // Test for truncation close to zero.
+ // At z <= ~1.5e-8, exp(-0.5 * z * z) / sqrt(2 pi) == 1 / sqrt(2 pi)
+ // and the PDF is constant. It can be approximated as a uniform distribution.
+ // Here the mean is computable but the variance computation -> 0.
+ // The epsilons for the variance allow the test to pass if the second moment
+ // uses a uniform distribution approximation: (b^3 - a^3) / (3b - 3a).
+ // This is not done at present and the variance computes incorrectly and close to 0.
+ // The largest span covers only 5.8242e-8 of the probability range of the parent normal
+ // and these are not practical truncations.
+ "-7.299454196351098e-8, 7.299454196351098e-8, 0, 1.77606771882092042827020676955306864e-15, 1e-15, -1e-15",
+ "-7.299454196351098e-8, 3.649727098175549e-8, -1.82486354908777262111748030604612676e-08, 9.99038091836768051420202283759953002e-16, 1e-15, -1e-15",
+ "-7.299454196351098e-8, 1.8248635490877744e-8, -2.7372953236316597672674778496667655e-08, 6.93776452664422342699175710737901419e-16, 1e-15, -2e-15",
+ "-7.299454196351098e-8, 0, -3.64972709817554726791073610445021429e-08, 4.44016929705230343732389204118195096e-16, 1e-15, -2e-15",
+ "-7.299454196351098e-8, -1.8248635490877744e-8, -4.56215887271943497112157190855901547e-08, 2.49759522957641055973442997155578316e-16, 3e-10, -5e-9",
+ "-7.299454196351098e-8, -3.649727098175549e-8, -5.47459064726332272497430210977513379e-08, 1.11004232421306844799494326433537718e-16, 3e-10, -2e-8",
+ "-3.649727098175549e-8, 3.649727098175549e-8, 0, 4.44016929705230343602092590994317462e-16, 1e-15, -1e-15",
+ "-3.649727098175549e-8, 1.8248635490877744e-8, -9.12431774543886994224314381693928319e-09, 2.49759522959192087703300220816741702e-16, 1e-15, -1e-15",
+ "-3.649727098175549e-8, 0, -1.82486354908777424165810069993271136e-08, 1.11004232426307600672725101668733634e-16, 1e-15, -2e-15",
+ "-3.649727098175549e-8, -1.8248635490877744e-8, -2.73729532363166159037567578213044937e-08, 2.77510581119222125912321725734620803e-17, 3e-10, -2e-8",
+ "-1.8248635490877744e-8, 1.8248635490877744e-8, 0, 1.11004232426307600757649604002272128e-16, 1e-15, -1e-15",
+ "-1.8248635490877744e-8, 9.124317745438872e-9, -4.5621588727194358257035396943085424e-09, 6.24398807397980267185125584627689296e-17, 1e-15, -1e-15",
+ "-1.8248635490877744e-8, 0, -9.12431774543887196791891930929818729e-09, 2.77510581065769011631256630479419225e-17, 1e-15, -1e-15",
+ "-9.124317745438872e-9, 9.124317745438872e-9, 0, 2.77510581065769013145632047586655539e-17, 1e-15, -1e-15",
+ "-9.124317745438872e-9, 4.562158872719436e-9, -2.28107943635971801967451582038414367e-09, 1.56099701849495071020338587207547036e-17, 1e-15, -1e-15",
+ "-9.124317745438872e-9, 0, -4.5621588727194360789130116308534264e-09, 6.93776452664422554954584023114952882e-18, 1e-15, -1e-15",
+
+ // The variance method is inaccurate at this extreme.
+ // Spans p=8.858e-17 of the parent normal distribution
+ "0, 2.220446049250313e-16, 1.11022302462515654042363166809081572e-16, 4.14074938043255708407035257655783112e-33, 1e-15, -1e-2",
+ })
+ void testAdditionalMoments(double lower, double upper,
+ double mean, double variance,
+ double meanRelativeError, double varianceRelativeError) {
+ assertMean(lower, upper, mean, meanRelativeError);
+ if (varianceRelativeError < 0) {
+ // Known problem case.
+ // Allow small absolute variances using an absolute threshold of
+ // machine epsilon (2^-52) * 1.5. Any true variance approaching machine epsilon
+ // is allowed to be computed as small or zero but cannot be too large.
+ final double var = TruncatedNormalDistribution.variance(lower, upper);
+ Assertions.assertTrue(var >= 0, () -> "Variance is not positive: " + var);
+ Assertions.assertEquals(var, TruncatedNormalDistribution.variance(-upper, -lower));
+ TestUtils.assertEquals(variance, var,
+ DoubleTolerances.relative(-varianceRelativeError).or(DoubleTolerances.absolute(1.5 * 0x1.0p-52)),
+ () -> String.format("variance(%s, %s)", lower, upper));
+ } else {
+ assertVariance(lower, upper, variance, varianceRelativeError);
+ }
+ }
+
+ /**
+ * Assert the mean of the truncated normal distribution is within the provided relative error.
+ */
+ private static void assertMean(double lower, double upper, double expected, double eps) {
+ final double mean = TruncatedNormalDistribution.moment1(lower, upper);
+ Assertions.assertEquals(0 - mean, TruncatedNormalDistribution.moment1(-upper, -lower));
+ TestUtils.assertEquals(expected, mean, DoubleTolerances.relative(eps),
+ () -> String.format("mean(%s, %s)", lower, upper));
+ }
+
+ /**
+ * Assert the mean of the truncated normal distribution is within the provided relative error.
+ * Helper method using range [lower, upper] of the parent normal distribution with the specified
+ * mean and standard deviation.
+ */
+ private static void assertMean(double lower, double upper, double u, double s, double expected, double eps) {
+ final double a = (lower - u) / s;
+ final double b = (upper - u) / s;
+ final double mean = u + TruncatedNormalDistribution.moment1(a, b) * s;
+ TestUtils.assertEquals(expected, mean, DoubleTolerances.relative(eps),
+ () -> String.format("mean(%s, %s, %s, %s)", lower, upper, u, s));
+ }
+
+ /**
+ * Assert the variance of the truncated normal distribution is within the provided relative error.
+ */
+ private static void assertVariance(double lower, double upper, double expected, double eps) {
+ final double var = TruncatedNormalDistribution.variance(lower, upper);
+ Assertions.assertEquals(var, TruncatedNormalDistribution.variance(-upper, -lower));
+ TestUtils.assertEquals(expected, var, DoubleTolerances.relative(eps),
+ () -> String.format("variance(%s, %s)", lower, upper));
+ }
+
+ /**
+ * Assert the variance of the truncated normal distribution is within the provided relative error.
+ * Helper method using range [lower, upper] of the parent normal distribution with the specified
+ * mean and standard deviation.
+ */
+ private static void assertVariance(double lower, double upper, double u, double s, double expected, double eps) {
+ final double a = (lower - u) / s;
+ final double b = (upper - u) / s;
+ final double var = TruncatedNormalDistribution.variance(a, b) * s * s;
+ TestUtils.assertEquals(expected, var, DoubleTolerances.relative(eps),
+ () -> String.format("variance(%s, %s, %s, %s)", lower, upper, u, s));
+ }
}
diff --git a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.5.properties b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.5.properties
index 82ba9e2..dbd39a1 100644
--- a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.5.properties
+++ b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.5.properties
@@ -15,10 +15,10 @@
# Test a truncation range that is completely below the mean.
parameters = 0.0, 1.0, -Infinity, -5.0
-# Limited by the CDF inverse mapping test
-tolerance.relative = 1e-7
-mean = -5.186503967125854
-variance = 0.03269643461704752
+# Limited by the SF inverse mapping test
+tolerance.relative = 2e-8
+mean = -5.18650396712584211561650896200523673
+variance = 0.0326964346171122253453158077009336653
lower = -Infinity
upper = -5
# Computed using Python with SciPy v1.6.0.
diff --git a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.6.properties b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.6.properties
index ad35a08..964fcfc 100644
--- a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.6.properties
+++ b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.6.properties
@@ -15,11 +15,10 @@
# Test a truncation range that is completely above the mean.
parameters = 0.0, 1.0, 5.0, Infinity
-# Limited by the mean computed by scipy
-# TODO: This test case should be investigated to verify the low tolerance required.
-tolerance.relative = 1e-6
-mean = 5.186503967125854
-variance = 0.03269643461704752
+# Limited by the CDF inverse mapping test
+tolerance.relative = 6e-9
+mean = 5.18650396712584211561650896200523673
+variance = 0.0326964346171122253453158077009336653
lower = 5
upper = Infinity
# Computed using Python with SciPy v1.6.0: truncnorm(5, inf)
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 5084cb1..33269ea 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
@@ -14,20 +14,41 @@
# limitations under the License.
# Test a narrow truncation range.
+# Note:
+# This is a very narrow truncation.
+# The normalised value for the bound [(lower-u)/sd, (upper-u)/sd] is:
+# [-1.0101010039621505E-8, 1.0101010129336497E-8]
+# Here exp(-0.5*x*x) == 1 and the computation loses precision and variance -> 0
+
parameters = 7.1, 9.9, 7.0999999, 7.1000001
-# 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
+
+# Limited by the computation of variance (which -> 0). Solved with absolute tolerance.
+# Also note that the relative error is lowered to pass the tests for logpdf, cdf and sf.
+
+tolerance.relative = 5e-9
+tolerance.absolute = 1e-14
+
+# scipy computes mean = 7.0999999999999996, var = 1.1358413484763297e-07 (too high)
+
+# Computed using 128-bit quad precision test c++ program using Boost erf for the
+# error function with variations on the computation creates values:
+mean = 7.10000000000000008881784197001250783
+# Variance has different values.
+# Wikipedia formula
+# 3.33333332242280748659034271839979628e-15
+# Alternate formula to avoid cancellation
+# 3.33333332242280748952284174437850924e-15
+variance = 3.33333332242280748952284174437850924e-15
+
lower = 7.0999999
upper = 7.1000001
+
# Computed using Python with SciPy v1.6.0:
# mean, std, clip_a, clip_b = 7.1, 9.9, 7.0999999, 7.1000001
# a, b = (clip_a - mean) / std, (clip_b - mean) / std
# truncnorm.var(a, b, loc=mean, scale=std)
+
cdf.points = 7.0999999, 7.1, 7.1000001
cdf.values = 0, 0.5, 1
pdf.values = 5000000.00238838, 5000000.00238838, 5000000.00238838
+logpdf.values = 15.42494847087604981084, 15.42494847087604981084, 15.42494847087604981084
diff --git a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.8.properties b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.8.properties
new file mode 100644
index 0000000..2258149
--- /dev/null
+++ b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.truncatednormal.8.properties
@@ -0,0 +1,51 @@
+# Licensed to the Apache Software Foundation (ASF) under one or more
+# contributor license agreements. See the NOTICE file distributed with
+# this work for additional information regarding copyright ownership.
+# The ASF licenses this file to You under the Apache License, Version 2.0
+# (the "License"); you may not use this file except in compliance with
+# the License. You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+# Test a truncation range that is completely above the mean but not at an infinite bound
+parameters = 0.0, 1.0, 1.0, 3.0
+# Computed using Python with SciPy v1.7.1: truncnorm(1, 3)
+mean = 1.510049513243984
+variance = 0.17345290492412202
+lower = 1
+upper = 3
+cdf.points = \
+ 0,\
+ 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9,\
+ 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9,\
+ 3.0, 4
+cdf.values = \
+ 0. , 0. ,\
+ 0.14614373969386054708, 0.27707628554938984466,\
+ 0.39321464289643304202, 0.49520624553490927289,\
+ 0.58388382351767098566, 0.66021885674404623412,\
+ 0.72527594829998764947, 0.78017009730223874087,\
+ 0.82602841697393314391, 0.86395737262641647547,\
+ 0.89501614591112199637, 0.92019629967422933436,\
+ 0.94040754724196873315, 0.95646913702463720597,\
+ 0.96910615492812091087, 0.97894992212345721683,\
+ 0.98654161672181917009, 0.99233826278893277895,\
+ 0.99672029432308406616, 1. , 1.
+pdf.values = \
+ 0. , 1.53822305117971835919,\
+ 1.38489993418456891483, 1.23445291403092372029,\
+ 1.089400873018817828 , 0.95182687696337731076,\
+ 0.82335146775529599594, 0.70513069338904532657,\
+ 0.59787587548371834423, 0.501891101223343572 ,\
+ 0.41712384425398174592, 0.34322395575373698673,\
+ 0.27960647448289188688, 0.22551420861233653636,\
+ 0.18007675313763649161, 0.14236343172640528176,\
+ 0.1114285040919437697 , 0.08634778616397101314,\
+ 0.06624653531224142244, 0.05031902148341711734,\
+ 0.03784062141893289322, 0.02817353793573459941, 0.