From 5786c5c944d3d8c5a18742c71d1cba6dc1a5eeab Mon Sep 17 00:00:00 2001 From: Phil Steitz Date: Thu, 1 Sep 2011 01:24:37 +0000 Subject: [PATCH] Added erf(double,double) to Erf and used this to improve tail probability accuracy in NormalDistributionImpl. JIRA: MATH-364. Reported and patched by Christian Winter. git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1163888 13f79535-47bb-0310-9956-ffa450edef68 --- pom.xml | 3 ++ .../distribution/NormalDistributionImpl.java | 22 ++++++++- .../org/apache/commons/math/special/Erf.java | 45 ++++++++++++++++--- src/site/xdoc/changes.xml | 4 ++ .../apache/commons/math/special/ErfTest.java | 18 ++++++++ 5 files changed, 85 insertions(+), 7 deletions(-) diff --git a/pom.xml b/pom.xml index c092d885b..7e7627359 100644 --- a/pom.xml +++ b/pom.xml @@ -245,6 +245,9 @@ Jörg Weimar + + Christian Winter + Xiaogang Zhang diff --git a/src/main/java/org/apache/commons/math/distribution/NormalDistributionImpl.java b/src/main/java/org/apache/commons/math/distribution/NormalDistributionImpl.java index eab5978a6..8df7369ae 100644 --- a/src/main/java/org/apache/commons/math/distribution/NormalDistributionImpl.java +++ b/src/main/java/org/apache/commons/math/distribution/NormalDistributionImpl.java @@ -21,6 +21,7 @@ import java.io.Serializable; import org.apache.commons.math.MathException; import org.apache.commons.math.exception.NotStrictlyPositiveException; +import org.apache.commons.math.exception.NumberIsTooLargeException; import org.apache.commons.math.exception.util.LocalizedFormats; import org.apache.commons.math.special.Erf; import org.apache.commons.math.util.FastMath; @@ -40,8 +41,10 @@ public class NormalDistributionImpl extends AbstractContinuousDistribution public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; /** Serializable version identifier. */ private static final long serialVersionUID = 8589540077390120676L; - /** &sqrt;(2 π) */ + /** √(2 π) */ private static final double SQRT2PI = FastMath.sqrt(2 * FastMath.PI); + /** √(2) */ + private static final double SQRT2 = FastMath.sqrt(2.0); /** Mean of this distribution. */ private final double mean; /** Standard deviation of this distribution. */ @@ -125,7 +128,22 @@ public class NormalDistributionImpl extends AbstractContinuousDistribution if (FastMath.abs(dev) > 40 * standardDeviation) { return dev < 0 ? 0.0d : 1.0d; } - return 0.5 * (1 + Erf.erf(dev / (standardDeviation * FastMath.sqrt(2)))); + return 0.5 * (1 + Erf.erf(dev / (standardDeviation * SQRT2))); + } + + /** + * {@inheritDoc} + */ + @Override + public double cumulativeProbability(double x0, double x1) throws MathException { + if (x0 > x1) { + throw new NumberIsTooLargeException(LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT, + x0, x1, true); + } + final double denom = standardDeviation * SQRT2; + final double v0 = (x0 - mean) / denom; + final double v1 = (x1 - mean) / denom; + return 0.5 * Erf.erf(v0, v1); } /** diff --git a/src/main/java/org/apache/commons/math/special/Erf.java b/src/main/java/org/apache/commons/math/special/Erf.java index 8cbb8b5b1..f936be2f8 100644 --- a/src/main/java/org/apache/commons/math/special/Erf.java +++ b/src/main/java/org/apache/commons/math/special/Erf.java @@ -25,6 +25,19 @@ import org.apache.commons.math.util.FastMath; * @version $Id$ */ public class Erf { + + /** + * The number {@code X_CRIT} is used by {@link #erf(double, double)} internally. + * This number solves {@code erf(x)=0.5} within 1ulp. + * More precisely, the current implementations of + * {@link #erf(double)} and {@link #erfc(double)} satisfy:
+ * {@code erf(X_CRIT) < 0.5},
+ * {@code erf(Math.nextUp(X_CRIT) > 0.5},
+ * {@code erfc(X_CRIT) = 0.5}, and
+ * {@code erfc(Math.nextUp(X_CRIT) < 0.5} + */ + private static final double X_CRIT = 0.4769362762044697; + /** * Default constructor. Prohibit instantiation. */ @@ -54,11 +67,8 @@ public class Erf { if (FastMath.abs(x) > 40) { return x > 0 ? 1 : -1; } - double ret = Gamma.regularizedGammaP(0.5, x * x, 1.0e-15, 10000); - if (x < 0) { - ret = -ret; - } - return ret; + final double ret = Gamma.regularizedGammaP(0.5, x * x, 1.0e-15, 10000); + return x < 0 ? -ret : ret; } /** @@ -91,5 +101,30 @@ public class Erf { final double ret = Gamma.regularizedGammaQ(0.5, x * x, 1.0e-15, 10000); return x < 0 ? 2 - ret : ret; } + + /** + * Returns the difference between erf(x1) and erf(x2). + * + * The implementation uses either erf(double) or erfc(double) + * depending on which provides the most precise result. + * + * @param x1 the first value + * @param x2 the second value + * @return erf(x2) - erf(x1) + */ + public static double erf(double x1, double x2) { + if(x1 > x2) { + return -erf(x2, x1); + } + + return + x1 < -X_CRIT ? + x2 < 0.0 ? + erfc(-x2) - erfc(-x1) : + erf(x2) - erf(x1) : + x2 > X_CRIT && x1 > 0.0 ? + erfc(x1) - erfc(x2) : + erf(x2) - erf(x1); + } } diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml index 7e92a1465..9b1dcb59a 100644 --- a/src/site/xdoc/changes.xml +++ b/src/site/xdoc/changes.xml @@ -52,6 +52,10 @@ The type attribute can be add,update,fix,remove. If the output is not quite correct, check for invisible trailing spaces! --> + + Added erf(double,double) to Erf and used this to improve tail probability + accuracy in NormalDistributionImpl. + Enabled reseeding of the random generators used by EmpiricalDistributionImpl and ValueServer. Modified ValueServer to pass its RandomData instance to diff --git a/src/test/java/org/apache/commons/math/special/ErfTest.java b/src/test/java/org/apache/commons/math/special/ErfTest.java index 765bccc4f..f0abcafd0 100644 --- a/src/test/java/org/apache/commons/math/special/ErfTest.java +++ b/src/test/java/org/apache/commons/math/special/ErfTest.java @@ -195,4 +195,22 @@ public class ErfTest { TestUtils.assertRelativelyEquals(ref[i][1], result, 1E-13); } } + + /** + * Test the implementation of Erf.erf(double, double) for consistency with results + * obtained from Erf.erf(double) and Erf.erfc(double). + */ + @Test + public void testTwoArgumentErf() throws Exception { + double[] xi = new double[]{-2.0, -1.0, -0.9, -0.1, 0.0, 0.1, 0.9, 1.0, 2.0}; + for(double x1 : xi) { + for(double x2 : xi) { + double a = Erf.erf(x1, x2); + double b = Erf.erf(x2) - Erf.erf(x1); + double c = Erf.erfc(x1) - Erf.erfc(x2); + Assert.assertEquals(a, b, 1E-15); + Assert.assertEquals(a, c, 1E-15); + } + } + } }