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);
+ }
+ }
+ }
}