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

@@ -245,6 +245,9 @@
     <contributor>
       <name>J&#246;rg Weimar</name>
     </contributor>
+    <contributor>
+      <name>Christian Winter</name>
+    </contributor>
     <contributor>
       <name>Xiaogang Zhang</name>
     </contributor>

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);
     }
 
     /**

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:<br/>
+     * {@code erf(X_CRIT) < 0.5},<br/>
+     * {@code erf(Math.nextUp(X_CRIT) > 0.5},<br/>
+     * {@code erfc(X_CRIT) = 0.5}, and<br/>
+     * {@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);
+        final double ret = Gamma.regularizedGammaP(0.5, x * x, 1.0e-15, 10000);
-        if (x < 0) {
+        return x < 0 ? -ret : ret;
-            ret = -ret;
-        }
-        return 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);
+    }
 }
 

src/site/xdoc/changes.xml

@@ -52,6 +52,10 @@ The <action> type attribute can be add,update,fix,remove.
     If the output is not quite correct, check for invisible trailing spaces!
      -->
     <release version="3.0" date="TBD" description="TBD">
+      <action dev="psteitz" type="update" issue="MATH-364" due-to="Christian Winter">
+        Added erf(double,double) to Erf and used this to improve tail probability
+        accuracy in NormalDistributionImpl. 
+      </action>
       <action dev="psteitz" type="fix" issue="MATH-654">
         Enabled reseeding of the random generators used by EmpiricalDistributionImpl
         and ValueServer.  Modified ValueServer to pass its RandomData instance to

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