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 2022/01/26 20:04:20 UTC

[commons-statistics] 05/06: Add special cases for binomial distribution PMF

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 5fb858a5d1723faa48211378aef8639101975a20
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Wed Jan 26 19:41:07 2022 +0000

    Add special cases for binomial distribution PMF
    
    When x=0 or x=n the power function can be used to evaluate the PMF to 1
    ULP. This is more accurate than exp(logPMF) when the expected result is
    very small.
    
    Increased precision of reference test values.
    
    Added specific bounds test using exact pmf computation.
    
    Updated test tolerance from default 1e-12 to 1e-15.
    
    Added special case tests for p=0, p=1, n=0, n=1.
    
    Updated code to handle p=-0.0.
---
 .../distribution/BinomialDistribution.java         |  37 +++-
 .../distribution/SaddlePointExpansionUtils.java    |   6 +-
 .../distribution/BinomialDistributionTest.java     | 214 +++++++++++++++++++++
 .../distribution/test.binomial.1.properties        |  25 ++-
 .../distribution/test.binomial.2.properties        |   3 +-
 .../distribution/test.binomial.3.properties        |   3 +-
 .../distribution/test.binomial.4.properties        |  13 +-
 7 files changed, 289 insertions(+), 12 deletions(-)

diff --git a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/BinomialDistribution.java b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/BinomialDistribution.java
index 5b5a1c0..e12cb68 100644
--- a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/BinomialDistribution.java
+++ b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/BinomialDistribution.java
@@ -37,10 +37,17 @@ import org.apache.commons.numbers.gamma.RegularizedBeta;
  * @see <a href="https://mathworld.wolfram.com/BinomialDistribution.html">Binomial distribution (MathWorld)</a>
  */
 public final class BinomialDistribution extends AbstractDiscreteDistribution {
+    /** 1/2. */
+    private static final float HALF = 0.5f;
+
     /** The number of trials. */
     private final int numberOfTrials;
     /** The probability of success. */
     private final double probabilityOfSuccess;
+    /** Cached value for pmf(x=0). */
+    private final double pmf0;
+    /** Cached value for pmf(x=n). */
+    private final double pmfn;
 
     /**
      * @param trials Number of trials.
@@ -50,6 +57,16 @@ public final class BinomialDistribution extends AbstractDiscreteDistribution {
                                  double p) {
         probabilityOfSuccess = p;
         numberOfTrials = trials;
+        // Special pmf cases where the power function is more accurate:
+        //   (n choose k) == 1 for k=0, k=n
+        //   pmf = p^k (1-p)^(n-k)
+        // Note: This handles the edge case of n=0: pmf(k=0) = 1, else 0
+        if (probabilityOfSuccess >= HALF) {
+            pmf0 = Math.pow(1 - probabilityOfSuccess, numberOfTrials);
+        } else {
+            pmf0 = Math.exp(numberOfTrials * Math.log1p(-probabilityOfSuccess));
+        }
+        pmfn = Math.pow(probabilityOfSuccess, numberOfTrials);
     }
 
     /**
@@ -68,7 +85,8 @@ public final class BinomialDistribution extends AbstractDiscreteDistribution {
                                             trials);
         }
         ArgumentUtils.checkProbability(p);
-        return new BinomialDistribution(trials, p);
+        // Avoid p = -0.0 to avoid returning -0.0 for some probability computations.
+        return new BinomialDistribution(trials, Math.abs(p));
     }
 
     /**
@@ -92,7 +110,16 @@ public final class BinomialDistribution extends AbstractDiscreteDistribution {
     /** {@inheritDoc} */
     @Override
     public double probability(int x) {
-        return Math.exp(logProbability(x));
+        if (x < 0 || x > numberOfTrials) {
+            return 0;
+        } else if (x == 0) {
+            return pmf0;
+        } else if (x == numberOfTrials) {
+            return pmfn;
+        }
+        return Math.exp(SaddlePointExpansionUtils.logBinomialProbability(x,
+                        numberOfTrials, probabilityOfSuccess,
+                        1.0 - probabilityOfSuccess));
     }
 
     /** {@inheritDoc} **/
@@ -103,6 +130,8 @@ public final class BinomialDistribution extends AbstractDiscreteDistribution {
         } else if (x < 0 || x > numberOfTrials) {
             return Double.NEGATIVE_INFINITY;
         }
+        // Special cases for x=0, x=n
+        // are handled in the saddle point expansion
         return SaddlePointExpansionUtils.logBinomialProbability(x,
                 numberOfTrials, probabilityOfSuccess,
                 1.0 - probabilityOfSuccess);
@@ -115,6 +144,8 @@ public final class BinomialDistribution extends AbstractDiscreteDistribution {
             return 0.0;
         } else if (x >= numberOfTrials) {
             return 1.0;
+        } else if (x == 0) {
+            return pmf0;
         }
         return RegularizedBeta.complement(probabilityOfSuccess,
                                           x + 1.0, (double) numberOfTrials - x);
@@ -127,6 +158,8 @@ public final class BinomialDistribution extends AbstractDiscreteDistribution {
             return 1.0;
         } else if (x >= numberOfTrials) {
             return 0.0;
+        } else if (x == numberOfTrials - 1) {
+            return pmfn;
         }
         return RegularizedBeta.value(probabilityOfSuccess,
                                      x + 1.0, (double) numberOfTrials - x);
diff --git a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/SaddlePointExpansionUtils.java b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/SaddlePointExpansionUtils.java
index c465f42..31ff62d 100644
--- a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/SaddlePointExpansionUtils.java
+++ b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/SaddlePointExpansionUtils.java
@@ -155,14 +155,16 @@ final class SaddlePointExpansionUtils {
     static double logBinomialProbability(int x, int n, double p, double q) {
         if (x == 0) {
             if (p < ONE_TENTH) {
-                return -getDeviancePart(n, n * q) - n * p;
+                // Subtract from 0 avoids returning -0.0 for p=0.0
+                return 0.0 - getDeviancePart(n, n * q) - n * p;
             } else if (n == 0) {
                 return 0;
             }
             return n * Math.log(q);
         } else if (x == n) {
             if (q < ONE_TENTH) {
-                return -getDeviancePart(n, n * p) - n * q;
+                // Subtract from 0 avoids returning -0.0 for p=1.0
+                return 0.0 - getDeviancePart(n, n * p) - n * q;
             }
             return n * Math.log(p);
         }
diff --git a/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/BinomialDistributionTest.java b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/BinomialDistributionTest.java
index 261042b..8a46506 100644
--- a/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/BinomialDistributionTest.java
+++ b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/BinomialDistributionTest.java
@@ -16,8 +16,13 @@
  */
 package org.apache.commons.statistics.distribution;
 
+import java.math.BigDecimal;
+import java.math.BigInteger;
 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;
+import org.junit.jupiter.params.provider.ValueSource;
 
 /**
  * Test cases for {@link BinomialDistribution}.
@@ -46,6 +51,11 @@ class BinomialDistributionTest extends BaseDiscreteDistributionTest {
         return new String[] {"NumberOfTrials", "ProbabilityOfSuccess"};
     }
 
+    @Override
+    protected double getRelativeTolerance() {
+        return 1e-15;
+    }
+
     //-------------------- Additional test cases -------------------------------
 
     @Test
@@ -60,4 +70,208 @@ class BinomialDistributionTest extends BaseDiscreteDistributionTest {
             Assertions.assertEquals(trials / 2, p);
         }
     }
+
+    /**
+     * Test special case of probability of success 0.0.
+     */
+    @ParameterizedTest
+    @ValueSource(ints = {0, 1, 2, 3, 10})
+    void testProbabilityOfSuccess0(int n) {
+        // The sign of p should not matter.
+        // Exact equality checks no -0.0 values are generated.
+        for (final double p : new double[] {-0.0, 0.0}) {
+            final BinomialDistribution dist = BinomialDistribution.of(n, p);
+            for (int k = -1; k <= n + 1; k++) {
+                Assertions.assertEquals(k == 0 ? 1.0 : 0.0, dist.probability(k));
+                Assertions.assertEquals(k == 0 ? 0.0 : Double.NEGATIVE_INFINITY, dist.logProbability(k));
+                Assertions.assertEquals(k >= 0 ? 1.0 : 0.0, dist.cumulativeProbability(k));
+                Assertions.assertEquals(k >= 0 ? 0.0 : 1.0, dist.survivalProbability(k));
+            }
+        }
+    }
+
+    /**
+     * Test special case of probability of success 1.0.
+     */
+    @ParameterizedTest
+    @ValueSource(ints = {0, 1, 2, 3, 10})
+    void testProbabilityOfSuccess1(int n) {
+        final BinomialDistribution dist = BinomialDistribution.of(n, 1);
+        // Exact equality checks no -0.0 values are generated.
+        for (int k = -1; k <= n + 1; k++) {
+            Assertions.assertEquals(k == n ? 1.0 : 0.0, dist.probability(k));
+            Assertions.assertEquals(k == n ? 0.0 : Double.NEGATIVE_INFINITY, dist.logProbability(k));
+            Assertions.assertEquals(k >= n ? 1.0 : 0.0, dist.cumulativeProbability(k));
+            Assertions.assertEquals(k >= n ? 0.0 : 1.0, dist.survivalProbability(k));
+        }
+    }
+
+    @ParameterizedTest
+    @ValueSource(doubles = {0, 1, 0.01, 0.99, 1e-17, 0.3645257e-8, 0.123415276368128, 0.67834532657232434})
+    void testNumberOfTrials0(double p) {
+        final BinomialDistribution dist = BinomialDistribution.of(0, p);
+        // Edge case where the probability is ignored when computing the result
+        for (int k = -1; k <= 2; k++) {
+            Assertions.assertEquals(k == 0 ? 1.0 : 0.0, dist.probability(k));
+            Assertions.assertEquals(k == 0 ? 0.0 : Double.NEGATIVE_INFINITY, dist.logProbability(k));
+            Assertions.assertEquals(k >= 0 ? 1.0 : 0.0, dist.cumulativeProbability(k));
+            Assertions.assertEquals(k >= 0 ? 0.0 : 1.0, dist.survivalProbability(k));
+        }
+    }
+
+    @ParameterizedTest
+    @ValueSource(doubles = {0, 1, 0.01, 0.99, 1e-17, 0.3645257e-8, 0.123415276368128, 0.67834532657232434})
+    void testNumberOfTrials1(double p) {
+        final BinomialDistribution dist = BinomialDistribution.of(1, p);
+        // Edge case where the probability should be exact
+        Assertions.assertEquals(0.0, dist.probability(-1));
+        Assertions.assertEquals(1 - p, dist.probability(0));
+        Assertions.assertEquals(p, dist.probability(1));
+        Assertions.assertEquals(0.0, dist.probability(2));
+        Assertions.assertEquals(Double.NEGATIVE_INFINITY, dist.logProbability(-1));
+        // Current implementation does not use log1p so allow an error tolerance
+        TestUtils.assertEquals(Math.log1p(-p), dist.logProbability(0), DoubleTolerances.ulps(1));
+        Assertions.assertEquals(Math.log(p), dist.logProbability(1));
+        Assertions.assertEquals(Double.NEGATIVE_INFINITY, dist.logProbability(2));
+        Assertions.assertEquals(0.0, dist.cumulativeProbability(-1));
+        Assertions.assertEquals(1 - p, dist.cumulativeProbability(0));
+        Assertions.assertEquals(1.0, dist.cumulativeProbability(1));
+        Assertions.assertEquals(1.0, dist.cumulativeProbability(2));
+        Assertions.assertEquals(1.0, dist.survivalProbability(-1));
+        Assertions.assertEquals(p, dist.survivalProbability(0));
+        Assertions.assertEquals(0.0, dist.survivalProbability(1));
+        Assertions.assertEquals(0.0, dist.survivalProbability(2));
+    }
+
+    /**
+     * Special case for x=0.
+     * This hits cases where the SaddlePointExpansionUtils are not used for
+     * probability functions. It ensures the edge case handling in BinomialDistribution
+     * matches the original logic in the saddle point expansion. This x=0 logic is used
+     * by the related hypergeometric distribution and covered by test cases for that
+     * distribution ensuring it is correct.
+     */
+    @ParameterizedTest
+    @ValueSource(doubles = {0, 1, 0.01, 0.99, 1e-17, 0.3645257e-8, 0.123415276368128, 0.67834532657232434})
+    void testX0(double p) {
+        for (final int n : new int[] {0, 1, 10}) {
+            final BinomialDistribution dist = BinomialDistribution.of(n, p);
+            final double expected = SaddlePointExpansionUtils.logBinomialProbability(0, n, p, 1 - p);
+            Assertions.assertEquals(expected, dist.logProbability(0));
+        }
+    }
+
+    /**
+     * Test the probability functions at the lower and upper bounds when the p-values
+     * are very small. The expected results should be within 1 ULP of an exact
+     * computation for pmf(x=0) and pmf(x=n). These values can be computed using
+     * java.lang.Math functions.
+     *
+     * <p>The next value, e.g. pmf(x=1) and cdf(x=1), is asserted to the specified
+     * relative error tolerance. These values require computations using p and 1-p
+     * and are less exact.
+     *
+     * @param n Number of trials
+     * @param p Probability of success
+     * @param eps1 Relative error tolerance for pmf(x=1)
+     * @param epsn1 Relative error tolerance for pmf(x=n-1)
+     */
+    @ParameterizedTest
+    @CsvSource({
+        // Min p-value is shown for reference.
+        "100, 0.50, 8e-15, 8e-15", // 7.888609052210118E-31
+        "100, 0.01, 1e-15, 3e-14", // 1.000000000000002E-200
+        "100, 0.99, 4e-15, 1e-15", // 1.0000000000000887E-200
+        "140, 0.01, 1e-15, 2e-13", // 1.0000000000000029E-280
+        "140, 0.99, 2e-13, 1e-15", // 1.0000000000001244E-280
+        "50, 0.001, 1e-15, 5e-14", // 1.0000000000000011E-150
+        "50, 0.999, 3e-14, 1e-15", // 1.0000000000000444E-150
+    })
+    void testBounds(int n, double p, double eps1, double epsn1) {
+        final BinomialDistribution dist = BinomialDistribution.of(n, p);
+        final BigDecimal prob0 = binomialProbability(n, p, 0);
+        final BigDecimal probn = binomialProbability(n, p, n);
+        final double p0 = prob0.doubleValue();
+        final double pn = probn.doubleValue();
+
+        // Require very small non-zero probabilities to make the test difficult.
+        // Check using 2^-53 so that at least one p-value satisfies 1 - p == 1.
+        final double minp = Math.min(p0, pn);
+        Assertions.assertTrue(minp < 0x1.0p-53, () -> "Test should target small p-values: " + minp);
+        Assertions.assertTrue(minp > Double.MIN_NORMAL, () -> "Minimum P-value should not be sub-normal: " + minp);
+
+        // Almost exact at the bounds
+        final DoubleTolerance tol1 = DoubleTolerances.ulps(1);
+        TestUtils.assertEquals(p0, dist.probability(0), tol1, "pmf(0)");
+        TestUtils.assertEquals(pn, dist.probability(n), tol1, "pmf(n)");
+        // Consistent at the bounds
+        Assertions.assertEquals(dist.probability(0), dist.cumulativeProbability(0), "pmf(0) != cdf(0)");
+        Assertions.assertEquals(dist.probability(n), dist.survivalProbability(n - 1), "pmf(n) != sf(n-1)");
+
+        // Test probability and log probability are consistent.
+        // Avoid log when p-value is close to 1.
+        if (p0 < 0.9) {
+            TestUtils.assertEquals(Math.log(p0), dist.logProbability(0), tol1, "log(pmf(0)) != logpmf(0)");
+        } else {
+            TestUtils.assertEquals(p0, Math.exp(dist.logProbability(0)), tol1, "pmf(0) != exp(logpmf(0))");
+        }
+        if (pn < 0.9) {
+            TestUtils.assertEquals(Math.log(pn), dist.logProbability(n), tol1, "log(pmf(n)) != logpmf(n)");
+        } else {
+            TestUtils.assertEquals(pn, Math.exp(dist.logProbability(n)), tol1, "pmf(n) != exp(logpmf(n))");
+        }
+
+        // The next probability is accurate to the specified tolerance.
+        final BigDecimal prob1 = binomialProbability(n, p, 1);
+        final BigDecimal probn1 = binomialProbability(n, p, n - 1);
+        TestUtils.assertEquals(prob1.doubleValue(), dist.probability(1), DoubleTolerances.relative(eps1), "pmf(1)");
+        TestUtils.assertEquals(probn1.doubleValue(), dist.probability(n - 1), DoubleTolerances.relative(epsn1), "pmf(n-1)");
+
+        // Check the cumulative functions
+        final double cdf1 = prob0.add(prob1).doubleValue();
+        final double sfn2 = probn.add(probn1).doubleValue();
+        TestUtils.assertEquals(cdf1, dist.cumulativeProbability(1), DoubleTolerances.relative(eps1), "cmf(1)");
+        TestUtils.assertEquals(sfn2, dist.survivalProbability(n - 2), DoubleTolerances.relative(epsn1), "sf(n-2)");
+    }
+
+    /**
+     * Compute the binomial distribution probability mass function using exact
+     * arithmetic.
+     *
+     * <p>This has no error handling for invalid arguments.
+     *
+     * <p>Warning: BigDecimal has a limit on the size of the exponent for the power
+     * function. This method has not been extensively tested with very small
+     * p-values, large n or large k. Use of a MathContext to round intermediates may be
+     * required to reduce memory consumption. The binomial coefficient may not
+     * compute for large n and k ~ n/2.
+     *
+     * @param n Number of trials (must be positive)
+     * @param p Probability of success (in [0, 1])
+     * @param k Number of successes (must be positive)
+     * @return pmf(X=k)
+     */
+    private static BigDecimal binomialProbability(int n, double p, int k) {
+        final int nmk = n - k;
+        final int m = Math.min(k, nmk);
+        // Probability component: p^k * (1-p)^(n-k)
+        final BigDecimal bp = new BigDecimal(p);
+        final BigDecimal result = bp.pow(k).multiply(
+                 BigDecimal.ONE.subtract(bp).pow(nmk));
+        // Compute the binomial coefficient
+        // Simple edge cases first.
+        if (m == 0) {
+            return result;
+        } else if (m == 1) {
+            return result.multiply(new BigDecimal(n));
+        }
+        // See org.apache.commons.numbers.combinatorics.BinomialCoefficient
+        BigInteger nCk = BigInteger.ONE;
+        int i = n - m + 1;
+        for (int j = 1; j <= m; j++) {
+            nCk = nCk.multiply(BigInteger.valueOf(i)).divide(BigInteger.valueOf(j));
+            i++;
+        }
+        return new BigDecimal(nCk).multiply(result);
+    }
 }
diff --git a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.binomial.1.properties b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.binomial.1.properties
index 5064881..81c639d 100644
--- a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.binomial.1.properties
+++ b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.binomial.1.properties
@@ -23,12 +23,27 @@ upper = 10
 cdf.points = -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
 # Reference values are from R, version 2.15.3.
 cdf.values = \
-  0d, 5.9049e-06, 0.0001436859, 0.0015903864, 0.0105920784,  0.0473489874,\
-  0.1502683326, 0.3503892816, 0.6172172136, 0.8506916541, 0.9717524751, 1d, 1d
+  0.0000000000000000000e+00 5.9049000000000118537e-06 \
+  1.4368590000000025682e-04 1.5903864000000027951e-03 \
+  1.0592078400000015312e-02 4.7348987400000042136e-02 \
+  1.5026833260000005410e-01 3.5038928160000010203e-01 \
+  6.1721721360000014744e-01 8.5069165410000002758e-01 \
+  9.7175247509999995721e-01 1.0000000000000000000e+00 \
+  1.0000000000000000000e+00
 pmf.values = \
-  0d, 0.0000059049d, 0.000137781d, 0.0014467005,\
-  0.009001692, 0.036756909, 0.1029193452, 0.200120949, 0.266827932,\
-  0.2334744405, 0.121060821, 0.0282475249, 0d
+  0.0000000000000000000e+00 5.9049000000000059245e-06 \
+  1.3778100000000039404e-04 1.4467005000000031888e-03 \
+  9.0016920000000021085e-03 3.6756909000000038967e-02 \
+  1.0291934520000002584e-01 2.0012094900000007569e-01 \
+  2.6682793199999993439e-01 2.3347444049999999116e-01 \
+  1.2106082099999986024e-01 2.8247524899999980341e-02 \
+  0.0000000000000000000e+00
+sf.values = \
+  1.00000000000000000000 0.99999409509999992451 0.99985631409999997654 \
+  0.99840961360000002323 0.98940792160000001765 0.95265101259999995786 \
+  0.84973166739999994590 0.64961071839999995348 0.38278278639999985256 \
+  0.14930834590000002793 0.02824752489999998728 0.00000000000000000000 \
+  0.00000000000000000000
 icdf.points = \
   0, 0.001, 0.010, 0.025, 0.050, 0.100,\
   0.999, 0.990, 0.975, 0.950, 0.900, 1
diff --git a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.binomial.2.properties b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.binomial.2.properties
index dfa3c15..92bac80 100644
--- a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.binomial.2.properties
+++ b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.binomial.2.properties
@@ -25,7 +25,8 @@ upper = 100
 cdf.points = 100
 cdf.values = 1.0
 pmf.points = 0
-pmf.values = 1e-200
+# BigDecimal.ONE.subtract(new BigDecimal(0.99)).pow(100, MathContext.DECIMAL128)
+pmf.values = 1.000000000000088817841970016428095E-200
 # computed using R version 3.4.4
 cdf.hp.points = 82, 81
 cdf.hp.values = 1.4061271955993513664e-17, 6.1128083336354843707e-19
diff --git a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.binomial.3.properties b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.binomial.3.properties
index 46b8d97..09cb32c 100644
--- a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.binomial.3.properties
+++ b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.binomial.3.properties
@@ -23,7 +23,8 @@ upper = 100
 # Set values for the required fields.
 cdf.points = 100
 cdf.values = 1.0
-pmf.values = 1e-200
+# new BigDecimal(0.01).pow(100, MathContext.DECIMAL128)
+pmf.values = 1.000000000000002081668171172170658E-200
 # computed using R version 3.4.4
 sf.hp.points = 18, 19
 sf.hp.values = 6.1128083336353977038e-19, 2.4944165604029235392e-20
diff --git a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.binomial.4.properties b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.binomial.4.properties
index 1c4c3d6..e459342 100644
--- a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.binomial.4.properties
+++ b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.binomial.4.properties
@@ -14,7 +14,6 @@
 # limitations under the License.
 
 parameters = 10 0.3
-tolerance.relative = 1e-10
 # 10 * 0.3
 mean = 3
 # 10 * 0.3 * 0.7
@@ -33,6 +32,18 @@ pmf.values = \
   1.2106082099999991575e-01, 2.3347444049999999116e-01, 2.6682793199999993439e-01, 2.0012094900000007569e-01,\
   1.0291934520000002584e-01, 3.6756909000000004273e-02, 9.0016919999999864960e-03, 1.4467005000000008035e-03,\
   1.3778099999999990615e-04, 5.9048999999999949131e-06, 0.0000000000000000000e+00
+# sf(x=9) of 5.9048999999999915249e-06 replaced with exact computation
+# using BinomialDistributionTest.binomialProbability.
+# The values differ by 7 ulp.
+sf.values = \
+  1.0000000000000000000e+00 9.7175247509999995721e-01,\
+  8.5069165410000002758e-01 6.1721721360000025847e-01,\
+  3.5038928159999976897e-01 1.5026833259999988757e-01,\
+  4.7348987399999951931e-02 1.0592078399999996230e-02,\
+  1.5903864000000001930e-03 1.4368589999999998577e-04,\
+  0.000005904899999999997814748020630304745,\
+  0.0000000000000000000e+00,\
+  0.0000000000000000000e+00
 icdf.points = \
   0, 0.001, 0.010, 0.025, 0.050, 0.100,\
   0.999, 0.990, 0.975, 0.950, 0.900, 1