Re-implement the Clopper-Pearson interval using the beta distribution

This requires less conversion of arguments before and after the
underlying call to the regularized beta function.

Author: aherbert

commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/interval/ClopperPearsonInterval.java

@@ -16,7 +16,7 @@
  */
 package org.apache.commons.math4.legacy.stat.interval;
 
-import org.apache.commons.statistics.distribution.FDistribution;
+import org.apache.commons.statistics.distribution.BetaDistribution;
 
 /**
  * Implements the Clopper-Pearson method for creating a binomial proportion confidence interval.
@@ -37,22 +37,17 @@ public ConfidenceInterval createInterval(int numberOfTrials,
         double lowerBound = 0;
         double upperBound = 1;
 
-        final double alpha = 0.5 * (1 - confidenceLevel);
+        // alpha = 1 - confidence level
+        final double halfAlpha = 0.5 * (1 - confidenceLevel);
+        final int n = numberOfTrials;
+        final int x = numberOfSuccesses;
 
         if (numberOfSuccesses > 0) {
-            final FDistribution distributionLowerBound = FDistribution.of(2.0 * (numberOfTrials - numberOfSuccesses + 1),
-                                                                          2.0 * numberOfSuccesses);
-            final double fValueLowerBound = distributionLowerBound.inverseSurvivalProbability(alpha);
-            lowerBound = numberOfSuccesses /
-                (numberOfSuccesses + (numberOfTrials - numberOfSuccesses + 1) * fValueLowerBound);
+            lowerBound = BetaDistribution.of(x, n - x + 1).inverseCumulativeProbability(halfAlpha);
         }
 
         if (numberOfSuccesses < numberOfTrials) {
-            final FDistribution distributionUpperBound = FDistribution.of(2.0 * (numberOfSuccesses + 1),
-                                                                          2.0 * (numberOfTrials - numberOfSuccesses));
-            final double fValueUpperBound = distributionUpperBound.inverseSurvivalProbability(alpha);
-            upperBound = (numberOfSuccesses + 1) * fValueUpperBound /
-                (numberOfTrials - numberOfSuccesses + (numberOfSuccesses + 1) * fValueUpperBound);
+            upperBound = BetaDistribution.of(x + 1, n - x).inverseSurvivalProbability(halfAlpha);
         }
 
         return new ConfidenceInterval(lowerBound, upperBound, confidenceLevel);

commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/stat/interval/ClopperPearsonIntervalTest.java

@@ -92,16 +92,32 @@ public void testCase3() {
         check(expected, createBinomialConfidenceInterval().createInterval(trials, successes, confidenceLevel));
     }
 
+    @Test
+    public void testCase4() {
+        final int successes = 20;
+        final int trials = 400;
+        final double confidenceLevel = 0.95;
+
+        // Check using scipy.stats.beta
+        // k, n, alpha = 20, 400, 0.05
+        // p_l, p_u = beta.ppf([alpha / 2, 1 - alpha / 2], [k, k + 1], [n - k + 1, n - k])
+        final ConfidenceInterval expected = new ConfidenceInterval(0.030805241143265938,
+                                                                   0.07616697275514255,
+                                                                   confidenceLevel);
+
+        check(expected, createBinomialConfidenceInterval().createInterval(trials, successes, confidenceLevel));
+    }
+
     private void check(ConfidenceInterval expected,
                        ConfidenceInterval actual) {
         final double relTol = 1.0e-6; // Reasonable relative tolerance for floating point comparison
         // Compare bounds using a relative tolerance
         Assert.assertEquals(expected.getLowerBound(),
                             actual.getLowerBound(),
-                            relTol * (1.0 + Math.abs(expected.getLowerBound())));
+                            relTol * expected.getLowerBound());
         Assert.assertEquals(expected.getUpperBound(),
                             actual.getUpperBound(),
-                            relTol * (1.0 + Math.abs(expected.getUpperBound())));
+                            relTol * expected.getUpperBound());
         // The confidence level must be exact
         Assert.assertEquals(expected.getConfidenceLevel(),
                             actual.getConfidenceLevel(),