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/22 13:13:59 UTC
[commons-numbers] 04/05: NUMBERS-181: New public API for beta functions
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-numbers.git
commit e7b1d6c7a9a5d654189650a1dd30180ae27552f7
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Thu Jan 20 23:17:31 2022 +0000
NUMBERS-181: New public API for beta functions
---
.../org/apache/commons/numbers/gamma/Beta.java | 58 +++++++
.../commons/numbers/gamma/IncompleteBeta.java | 127 +++++++++++++++
.../commons/numbers/gamma/RegularizedBeta.java | 176 +++++++++++----------
.../org/apache/commons/numbers/gamma/BetaTest.java | 35 ++++
.../commons/numbers/gamma/IncompleteBetaTest.java | 112 +++++++++++++
.../commons/numbers/gamma/RegularizedBetaTest.java | 175 +++++++++++++++++---
6 files changed, 580 insertions(+), 103 deletions(-)
diff --git a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/Beta.java b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/Beta.java
new file mode 100644
index 0000000..7f5a9d9
--- /dev/null
+++ b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/Beta.java
@@ -0,0 +1,58 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.numbers.gamma;
+
+/**
+ * <a href="https://mathworld.wolfram.com/BetaFunction.html">Beta function</a>.
+ *
+ * <p>\[ B(a, b) = \frac{\Gamma(a)\ \Gamma(b)}{\Gamma(a+b)} = \frac{(a-1)!\ (b-1)!}{(a+b-1)!} \]
+ *
+ * <p>where \( \Gamma(z) \) is the gamma function.
+ *
+ * <p>This code has been adapted from the <a href="https://www.boost.org/">Boost</a>
+ * {@code c++} implementation {@code <boost/math/special_functions/beta.hpp>}.
+ *
+ * @see
+ * <a href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_beta/beta_function.html">
+ * Boost C++ Beta function</a>
+ * @since 1.1
+ */
+public final class Beta {
+
+ /** Private constructor. */
+ private Beta() {
+ // intentionally empty.
+ }
+
+ /**
+ * Computes the value of the
+ * <a href="https://mathworld.wolfram.com/BetaFunction.html">
+ * beta function</a> B(a, b).
+ *
+ * <p>\[ B(a, b) = \frac{\Gamma(a)\ \Gamma(b)}{\Gamma(a+b)} \]
+ *
+ * <p>where \( \Gamma(z) \) is the gamma function.
+ *
+ * @param a Parameter {@code a}.
+ * @param b Parameter {@code b}.
+ * @return the beta function \( B(a, b) \).
+ */
+ public static double value(double a,
+ double b) {
+ return BoostBeta.beta(a, b);
+ }
+}
diff --git a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/IncompleteBeta.java b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/IncompleteBeta.java
new file mode 100644
index 0000000..dc4911f
--- /dev/null
+++ b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/IncompleteBeta.java
@@ -0,0 +1,127 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.numbers.gamma;
+
+/**
+ * <a href="https://mathworld.wolfram.com/IncompleteBetaFunction.html">
+ * Incomplete Beta function</a>.
+ *
+ * <p>\[ B_x(a,b) = \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt \]
+ *
+ * <p>This code has been adapted from the <a href="https://www.boost.org/">Boost</a>
+ * {@code c++} implementation {@code <boost/math/special_functions/beta.hpp>}.
+ *
+ * @see
+ * <a href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_beta/ibeta_function.html">
+ * Boost C++ Incomplete Beta functions</a>
+ * @since 1.1
+ */
+public final class IncompleteBeta {
+
+ /** Private constructor. */
+ private IncompleteBeta() {
+ // intentionally empty.
+ }
+
+ /**
+ * Computes the value of the
+ * <a href="https://mathworld.wolfram.com/IncompleteBetaFunction.html">
+ * incomplete beta function</a> B(x, a, b).
+ *
+ * <p>\[ B_x(a,b) = \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt \]
+ *
+ * @param x Value.
+ * @param a Parameter {@code a}.
+ * @param b Parameter {@code b}.
+ * @return the incomplete beta function \( B_x(a, b) \).
+ * @throws ArithmeticException if the series evaluation fails to converge.
+ */
+ public static double value(double x,
+ double a,
+ double b) {
+ return BoostBeta.beta(a, b, x);
+ }
+
+ /**
+ * Computes the value of the
+ * <a href="https://mathworld.wolfram.com/IncompleteBetaFunction.html">
+ * incomplete beta function</a> B(x, a, b).
+ *
+ * <p>\[ B_x(a,b) = \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt \]
+ *
+ * @param x the value.
+ * @param a Parameter {@code a}.
+ * @param b Parameter {@code b}.
+ * @param epsilon Tolerance in series evaluation.
+ * @param maxIterations Maximum number of iterations in series evaluation.
+ * @return the incomplete beta function \( B_x(a, b) \).
+ * @throws ArithmeticException if the series evaluation fails to converge.
+ */
+ public static double value(double x,
+ final double a,
+ final double b,
+ double epsilon,
+ int maxIterations) {
+ return BoostBeta.beta(a, b, x, new Policy(epsilon, maxIterations));
+ }
+
+ /**
+ * Computes the complement of the
+ * <a href="https://mathworld.wolfram.com/IncompleteBetaFunction.html">
+ * incomplete beta function</a> B(x, a, b).
+ *
+ * <p>\[ B(a, b) - B_x(a,b) = B_{1-x}(b, a) \]
+ *
+ * <p>where \( B(a, b) \) is the beta function.
+ *
+ * @param x Value.
+ * @param a Parameter {@code a}.
+ * @param b Parameter {@code b}.
+ * @return the complement of the incomplete beta function \( B(a, b) - B_x(a, b) \).
+ * @throws ArithmeticException if the series evaluation fails to converge.
+ */
+ public static double complement(double x,
+ double a,
+ double b) {
+ return BoostBeta.betac(a, b, x);
+ }
+
+ /**
+ * Computes the complement of the
+ * <a href="https://mathworld.wolfram.com/IncompleteBetaFunction.html">
+ * incomplete beta function</a> B(x, a, b).
+ *
+ * <p>\[ B(a, b) - B_x(a,b) = B_{1-x}(b, a) \]
+ *
+ * <p>where \( B(a, b) \) is the beta function.
+ *
+ * @param x the value.
+ * @param a Parameter {@code a}.
+ * @param b Parameter {@code b}.
+ * @param epsilon Tolerance in series evaluation.
+ * @param maxIterations Maximum number of iterations in series evaluation.
+ * @return the complement of the incomplete beta function \( B(a, b) - B_x(a, b) \).
+ * @throws ArithmeticException if the series evaluation fails to converge.
+ */
+ public static double complement(double x,
+ final double a,
+ final double b,
+ double epsilon,
+ int maxIterations) {
+ return BoostBeta.betac(a, b, x, new Policy(epsilon, maxIterations));
+ }
+}
diff --git a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/RegularizedBeta.java b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/RegularizedBeta.java
index 04cd683..ec3bbbe 100644
--- a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/RegularizedBeta.java
+++ b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/RegularizedBeta.java
@@ -16,18 +16,22 @@
*/
package org.apache.commons.numbers.gamma;
-import org.apache.commons.numbers.fraction.ContinuedFraction;
-
/**
- * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
+ * <a href="https://mathworld.wolfram.com/RegularizedBetaFunction.html">
* Regularized Beta function</a>.
- * <p>
- * This class is immutable.
- * </p>
+ *
+ * <p>\[ I_x(a,b) = \frac{1}{B(a, b)} \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt \]
+ *
+ * <p>where \( B(a, b) \) is the beta function.
+ *
+ * <p>This code has been adapted from the <a href="https://www.boost.org/">Boost</a>
+ * {@code c++} implementation {@code <boost/math/special_functions/beta.hpp>}.
+ *
+ * @see
+ * <a href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_beta/ibeta_function.html">
+ * Boost C++ Incomplete Beta functions</a>
*/
public final class RegularizedBeta {
- /** Maximum allowed numerical error. */
- private static final double DEFAULT_EPSILON = 1e-14;
/** Private constructor. */
private RegularizedBeta() {
@@ -36,105 +40,115 @@ public final class RegularizedBeta {
/**
* Computes the value of the
- * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
+ * <a href="https://mathworld.wolfram.com/RegularizedBetaFunction.html">
* regularized beta function</a> I(x, a, b).
*
+ * <p>\[ I_x(a,b) = \frac{1}{B(a, b)} \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt \]
+ *
+ * <p>where \( B(a, b) \) is the beta function.
+ *
* @param x Value.
* @param a Parameter {@code a}.
* @param b Parameter {@code b}.
- * @return the regularized beta function I(x, a, b).
- * @throws ArithmeticException if the algorithm fails to converge.
+ * @return the regularized beta function \( I_x(a, b) \).
+ * @throws ArithmeticException if the series evaluation fails to converge.
*/
public static double value(double x,
double a,
double b) {
- return value(x, a, b, DEFAULT_EPSILON, Integer.MAX_VALUE);
+ return BoostBeta.ibeta(a, b, x);
}
-
/**
* Computes the value of the
- * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
+ * <a href="https://mathworld.wolfram.com/RegularizedBetaFunction.html">
* regularized beta function</a> I(x, a, b).
*
- * The implementation of this method is based on:
- * <ul>
- * <li>
- * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
- * Regularized Beta Function</a>.
- * </li>
- * <li>
- * <a href="http://functions.wolfram.com/06.21.10.0001.01">
- * Regularized Beta Function</a>.
- * </li>
- * </ul>
+ * <p>\[ I_x(a,b) = \frac{1}{B(a, b)} \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt \]
+ *
+ * <p>where \( B(a, b) \) is the beta function.
*
* @param x the value.
* @param a Parameter {@code a}.
* @param b Parameter {@code b}.
- * @param epsilon When the absolute value of the nth item in the
- * series is less than epsilon the approximation ceases to calculate
- * further elements in the series.
- * @param maxIterations Maximum number of "iterations" to complete.
- * @return the regularized beta function I(x, a, b).
- * @throws ArithmeticException if the algorithm fails to converge.
+ * @param epsilon Tolerance in series evaluation.
+ * @param maxIterations Maximum number of iterations in series evaluation.
+ * @return the regularized beta function \( I_x(a, b) \).
+ * @throws ArithmeticException if the series evaluation fails to converge.
*/
public static double value(double x,
final double a,
final double b,
double epsilon,
int maxIterations) {
- if (Double.isNaN(x) ||
- Double.isNaN(a) ||
- Double.isNaN(b) ||
- x < 0 ||
- x > 1 ||
- a <= 0 ||
- b <= 0) {
- return Double.NaN;
- } else if (x == 0) {
- return 0.0;
- } else if (x == 1) {
- return 1.0;
- } else if (b == 1) {
- return Math.pow(x, a);
- } else if (a == 1) {
- // 1 - (1-x)^b
- // When x >= 0.5 then 1-x is exact
- if (x >= 0.5) {
- return 1.0 - Math.pow(1 - x, b);
- }
- // Otherwise support very small x with high-precision exp/log functions
- // 1 - exp(b * log(1-x))
- return -Math.expm1(b * Math.log1p(-x));
- } else if (x > (a + 1) / (2 + b + a) &&
- 1 - x <= (b + 1) / (2 + b + a)) {
- return 1 - value(1 - x, b, a, epsilon, maxIterations);
- } else {
- final ContinuedFraction fraction = new ContinuedFraction() {
- /** {@inheritDoc} */
- @Override
- protected double getA(int n, double x) {
- if ((n & 0x1) == 0) { // even
- final double m = n / 2d;
- return (m * (b - m) * x) /
- ((a + (2 * m) - 1) * (a + (2 * m)));
- }
- final double m = (n - 1d) / 2d;
- return -((a + m) * (a + b + m) * x) /
- ((a + (2 * m)) * (a + (2 * m) + 1));
- }
+ return BoostBeta.ibeta(a, b, x, new Policy(epsilon, maxIterations));
+ }
- /** {@inheritDoc} */
- @Override
- protected double getB(int n, double x) {
- return 1;
- }
- };
+ /**
+ * Computes the complement of the
+ * <a href="https://mathworld.wolfram.com/RegularizedBetaFunction.html">
+ * regularized beta function</a> I(x, a, b).
+ *
+ * <p>\[ 1 - I_x(a,b) = I_{1-x}(b, a) \]
+ *
+ * @param x Value.
+ * @param a Parameter {@code a}.
+ * @param b Parameter {@code b}.
+ * @return the complement of the regularized beta function \( 1 - I_x(a, b) \).
+ * @throws ArithmeticException if the series evaluation fails to converge.
+ * @since 1.1
+ */
+ public static double complement(double x,
+ double a,
+ double b) {
+ return BoostBeta.ibetac(a, b, x);
+ }
- return Math.exp((a * Math.log(x)) + (b * Math.log1p(-x)) -
- Math.log(a) - LogBeta.value(a, b)) /
- fraction.evaluate(x, epsilon, maxIterations);
- }
+ /**
+ * Computes the complement of the
+ * <a href="https://mathworld.wolfram.com/RegularizedBetaFunction.html">
+ * regularized beta function</a> I(x, a, b).
+ *
+ * <p>\[ 1 - I_x(a,b) = I_{1-x}(b, a) \]
+ *
+ * @param x the value.
+ * @param a Parameter {@code a}.
+ * @param b Parameter {@code b}.
+ * @param epsilon Tolerance in series evaluation.
+ * @param maxIterations Maximum number of iterations in series evaluation.
+ * @return the complement of the regularized beta function \( 1 - I_x(a, b) \).
+ * @throws ArithmeticException if the series evaluation fails to converge.
+ * @since 1.1
+ */
+ public static double complement(double x,
+ final double a,
+ final double b,
+ double epsilon,
+ int maxIterations) {
+ return BoostBeta.ibetac(a, b, x, new Policy(epsilon, maxIterations));
+ }
+
+ /**
+ * Computes the derivative of the
+ * <a href="https://mathworld.wolfram.com/RegularizedBetaFunction.html">
+ * regularized beta function</a> I(x, a, b).
+ *
+ * <p>\[ \frac{\delta}{\delta x} I_x(a,b) = \frac{(1-x)^{b-1} x^{a-1}}{B(a, b)} \]
+ *
+ * <p>where \( B(a, b) \) is the beta function.
+ *
+ * <p>This function has uses in some statistical distributions.
+ *
+ * @param x Value.
+ * @param a Parameter {@code a}.
+ * @param b Parameter {@code b}.
+ * @return the derivative of the regularized beta function \( I_x(a, b) \).
+ * @throws ArithmeticException if the series evaluation fails to converge.
+ * @since 1.1
+ */
+ public static double derivative(double x,
+ double a,
+ double b) {
+ return BoostBeta.ibetaDerivative(a, b, x);
}
}
diff --git a/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/BetaTest.java b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/BetaTest.java
new file mode 100644
index 0000000..bc17514
--- /dev/null
+++ b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/BetaTest.java
@@ -0,0 +1,35 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.numbers.gamma;
+
+import org.junit.jupiter.params.ParameterizedTest;
+import org.junit.jupiter.params.provider.CsvFileSource;
+
+/**
+ * Tests for {@link Beta}.
+ *
+ * <p>The class directly calls the methods in {@link BoostBeta}. This test ensures
+ * the arguments are passed through correctly. Accuracy of the function is tested
+ * in {@link BoostBetaTest}.
+ */
+class BetaTest {
+ @ParameterizedTest
+ @CsvFileSource(resources = "beta_med_data.csv")
+ void testBeta(double a, double b, double beta) {
+ TestUtils.assertEquals(beta, Beta.value(a, b), 200);
+ }
+}
diff --git a/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/IncompleteBetaTest.java b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/IncompleteBetaTest.java
new file mode 100644
index 0000000..2ab1358
--- /dev/null
+++ b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/IncompleteBetaTest.java
@@ -0,0 +1,112 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.numbers.gamma;
+
+import org.junit.jupiter.api.Assertions;
+import org.junit.jupiter.params.ParameterizedTest;
+import org.junit.jupiter.params.provider.CsvFileSource;
+import org.junit.jupiter.params.provider.CsvSource;
+
+/**
+ * Tests for {@link IncompleteBeta}.
+ *
+ * <p>The class directly calls the methods in {@link BoostBeta}. This test ensures
+ * the arguments are passed through correctly. Accuracy of the function is tested
+ * in {@link BoostBetaTest}.
+ */
+class IncompleteBetaTest {
+ private static final double EPS = Policy.getDefault().getEps();
+ private static final int MAX_ITER = Policy.getDefault().getMaxIterations();
+
+ @ParameterizedTest
+ @CsvFileSource(resources = "ibeta_med_data.csv")
+ void testIBeta(double a, double b, double x, double beta) {
+ TestUtils.assertEquals(beta, IncompleteBeta.value(x, a, b), 190);
+ }
+
+ @ParameterizedTest
+ @CsvFileSource(resources = "ibeta_med_data.csv")
+ void testIBetaComplement(double a, double b, double x, double beta, double betac) {
+ TestUtils.assertEquals(betac, IncompleteBeta.complement(x, a, b), 130);
+ }
+
+ /**
+ * Test the incomplete beta function uses the policy containing the epsilon and
+ * maximum iterations for series evaluations. The data targets each method computed
+ * using a series component to check the policy is not ignored.
+ *
+ * @see BoostBetaTest#testIBetaPolicy(double, double, double, double, double, double, double)
+ */
+ @ParameterizedTest
+ @CsvSource(value = {
+ // Target ibetaSeries
+ "4.201121919322759E-5,2.1881177963223308E-4,0.6323960423469543,23803.707603529016,4569.595256369948,0.838947362634037,0.16105263736596298",
+ "0.22512593865394592,0.4898320138454437,0.7996731996536255,4.764013990849158,0.9820281524501243,0.8290948573018541,0.17090514269814597",
+ "4.623167842510156E-5,4.340034502092749E-5,0.135563462972641,21628.337679706918,23043.143905699715,0.48416432390665337,0.5158356760933467",
+ "2.9415024982881732E-5,4.1924233664758503E-4,0.3082362115383148,33995.42298781772,2386.0630988783787,0.9344154580933705,0.0655845419066295",
+ "1.184685606858693E-5,0.015964560210704803,0.3082362115383148,84409.7658651171,63.42758275377908,0.9992491395179357,7.508604820642986E-4",
+ "3.529437162796967E-5,2.2326681573758833E-5,0.6323960423469543,28333.671885777032,44788.91695876111,0.3874817937042091,0.6125182062957909",
+ "0.06715317070484161,2.306319236755371,0.9084427952766418,13.787272071732604,0.001859217475218142,0.999865167903893,1.3483209610697333E-4",
+ "0.3183284401893616,3.165504217147827,0.07764927297830582,1.3374998709679642,0.6794195418585712,0.6631399660602084,0.3368600339397915",
+ "0.15403440594673157,4.049813747406006,0.34629878401756287,4.872033861103044,0.08561968850485947,0.9827297959310547,0.01727020406894529",
+ "1.3317101001739502,0.7650398015975952,0.6445860862731934,0.47144799136487586,0.5594135526519237,0.4573339592510717,0.5426660407489283",
+ "0.11902070045471191,7.269547462463379,0.19963125884532928,6.225047194692518,0.08420413075357451,0.9866538632858122,0.013346136714187858",
+ "2.664715051651001,0.6914005279541016,0.8443243503570557,0.3338388990912521,0.3587830340198169,0.4819929648946269,0.518007035105373",
+ "1.0562920570373535,6.812713623046875,0.8258343935012817,0.12732498812245932,9.807107058749213E-7,0.9999922976379849,7.702362015088557E-6",
+ "1.7118667364120483,3.0191311836242676,0.07594671100378036,0.0064151684204504875,0.10850933283432233,0.05582072012850319,0.9441792798714969",
+ // Target betaSmallBLargeASeries
+ "0.04634225741028786,0.34317314624786377,0.24176712334156036,20.363670894714268,3.6307737402387885,0.8486827348798152,0.15131726512018484",
+ "2.113992923113983E-5,1.7535277947899885E-5,0.8350250720977783,47305.46963235176,57026.27394012006,0.4534139659948463,0.5465860340051537",
+ "4.005068331025541E-5,42.84983825683594,0.12707412242889404,24964.03974453176,4.764518849491958E-4,0.9999999809144722,1.9085527852527327E-8",
+ "67.90167999267578,0.8324270844459534,0.9676981568336487,0.002669338211636283,0.031098353139101195,0.0790500654578503,0.9209499345421497",
+ "0.395370751619339,0.004023698624223471,0.9058013558387756,4.307485901473997,246.2348442100959,0.017192647244702382,0.9828073527552976",
+ "3.444607973098755,66.36054992675781,0.09654488414525986,1.4741027361579568E-6,8.307589573110104E-8,0.9466497330301025,0.05335026696989754",
+ "1.0665277242660522,4.985442161560059,0.2400285303592682,0.12523918055824373,0.0476465037156954,0.7244045745268357,0.2755954254731643",
+ // Target ibetaFraction2
+ "1.319732904434204,4.903014659881592,0.33251503109931946,0.0837704419861451,0.021604794441302123,0.7949727547593459,0.20502724524065405",
+ "485.7690734863281,190.16734313964844,0.6323960423469543,7.8023253024461E-182,7.885435919806278E-176,9.894592590329194E-7,0.9999990105407409",
+ })
+ void testIBetaPolicy(double a, double b, double x, double beta, double betac, double ibeta, double ibetac) {
+ // Low iterations should fail to converge
+ Assertions.assertThrows(ArithmeticException.class, () -> IncompleteBeta.value(x, a, b, EPS, 1), "ibeta");
+ Assertions.assertThrows(ArithmeticException.class, () -> IncompleteBeta.complement(x, a, b, EPS, 1), "ibetac");
+
+ // Low epsilon should not be as accurate
+
+ // Ignore infinite
+ if (Double.isFinite(beta)) {
+ final double u1 = IncompleteBeta.value(x, a, b);
+ final double u2 = IncompleteBeta.value(x, a, b, 1e-3, MAX_ITER);
+ assertCloser("beta", beta, u1, u2);
+ }
+ if (Double.isFinite(betac)) {
+ final double l1 = IncompleteBeta.complement(x, a, b);
+ final double l2 = IncompleteBeta.complement(x, a, b, 1e-3, MAX_ITER);
+ assertCloser("betac", betac, l1, l2);
+ }
+ }
+
+ /**
+ * Assert x is closer to the expected result than y.
+ */
+ private static void assertCloser(String msg, double expected, double x, double y) {
+ final double dx = Math.abs(expected - x);
+ final double dy = Math.abs(expected - y);
+ Assertions.assertTrue(dx < dy,
+ () -> String.format("%s %s : %s (%s) : %s (%s)", msg, expected, x, dx, y, dy));
+ }
+}
diff --git a/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/RegularizedBetaTest.java b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/RegularizedBetaTest.java
index 509f6af..5262587 100644
--- a/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/RegularizedBetaTest.java
+++ b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/RegularizedBetaTest.java
@@ -19,35 +19,39 @@ package org.apache.commons.numbers.gamma;
import org.junit.jupiter.api.Assertions;
import org.junit.jupiter.api.Test;
import org.junit.jupiter.params.ParameterizedTest;
+import org.junit.jupiter.params.provider.CsvFileSource;
import org.junit.jupiter.params.provider.CsvSource;
/**
* Tests for {@link RegularizedBeta}.
*/
class RegularizedBetaTest {
+ private static final double EPS = Policy.getDefault().getEps();
+ private static final int MAX_ITER = Policy.getDefault().getMaxIterations();
@ParameterizedTest
@CsvSource({
- "0.0, 0, 1, 2",
- "1.0, 1, 1, 2",
- "0.75, 0.5, 1, 2",
+ "0, 1, 2, 0",
+ "1, 1, 2, 1",
+ "0.5, 1, 2, 0.75",
// Invalid variants
- "NaN, NaN, 1, 2",
- "NaN, 0.5, NaN, 2",
- "NaN, 0.5, 1, NaN",
- "NaN, -0.5, 1, 2",
- "NaN, 0.5, -1, 2",
- "NaN, 0.5, 1, -2",
- "NaN, 0.5, 0, 2",
- "NaN, 0.5, 1, 0",
- "NaN, 1.5, 1, 2",
+ "1, 2, NaN, NaN",
+ "0.5, 2, NaN, NaN",
+ "0.5, 1, NaN, NaN",
+ "-0.5, 1, 2, NaN",
+ "0.5, -1, 2, NaN",
+ "0.5, 1, -2, NaN",
+ "0.5, 0, 0, NaN",
+ "1.5, 1, 2, NaN",
+ // Special case where a xor b is zero
+ "0.5, 0, 2, 1",
+ "0.5, 1, 0, 0",
})
- void testRegularizedBeta(double expected,
- double x,
- double a,
- double b) {
- final double actual = RegularizedBeta.value(x, a, b);
- Assertions.assertEquals(expected, actual, 1e-15);
+ void testRegularizedBetaArguments(double x,
+ double a,
+ double b,
+ double expected) {
+ assertRegularizedBeta(x, a, b, expected);
}
@Test
@@ -84,15 +88,15 @@ class RegularizedBetaTest {
@Test
void testZeroAndOne() {
// NUMBERS: 170
- Assertions.assertEquals(1.0, RegularizedBeta.value(1.0, 1e17, 0.5));
- Assertions.assertEquals(0.0, RegularizedBeta.value(0.0, 1e17, 0.5));
+ assertRegularizedBeta(0.0, 1e17, 0.5, 0);
+ assertRegularizedBeta(1.0, 1e17, 0.5, 1);
// a and b do not matter
final double[] v = {0.1, 0.5, 1, 2, 10};
for (final double a : v) {
for (final double b : v) {
- Assertions.assertEquals(0.0, RegularizedBeta.value(0.0, b, a));
- Assertions.assertEquals(1.0, RegularizedBeta.value(1.0, b, a));
+ assertRegularizedBeta(0, a, b, 0);
+ assertRegularizedBeta(1, a, b, 1);
}
}
}
@@ -130,4 +134,131 @@ class RegularizedBetaTest {
}
}
}
+
+ @ParameterizedTest
+ @CsvSource({
+ // Generated using matlab's betainc function with optional 'upper'
+ // argument for the complement
+ "0.125, 0.75, 0.25, 0.065817016037708564, 0.93418298396229149",
+ "0.125, 0.75, 8.75, 0.78161774393594907, 0.21838225606405093",
+ "0.125, 0.75, 13, 0.88429832262475594, 0.11570167737524407",
+ "0.125, 5.5, 0.25, 8.8583485327920987e-07, 0.99999911416514675",
+ "0.125, 5.5, 8.75, 0.0080298818621658066, 0.9919701181378342",
+ "0.125, 5.5, 13, 0.031403972369053644, 0.96859602763094632",
+ "0.125, 15.75, 0.25, 2.275822511739568e-16, 0.99999999999999978",
+ "0.125, 15.75, 8.75, 1.1608725295814249e-09, 0.99999999883912749",
+ "0.125, 15.75, 13, 3.5489227844330215e-08, 0.99999996451077211",
+ "0.375, 0.75, 0.25, 0.16604052887639467, 0.83395947112360536",
+ "0.375, 0.75, 8.75, 0.9905382314971497, 0.0094617685028503297",
+ "0.375, 0.75, 13, 0.99882118884718074, 0.0011788111528192643",
+ "0.375, 5.5, 0.25, 0.00045786386155093555, 0.99954213613844911",
+ "0.375, 5.5, 8.75, 0.48323155833692188, 0.51676844166307812",
+ "0.375, 5.5, 13, 0.77664105547610296, 0.22335894452389699",
+ "0.375, 15.75, 0.25, 9.3965393190826217e-09, 0.99999999060346068",
+ "0.375, 15.75, 8.75, 0.0035432647390343341, 0.99645673526096568",
+ "0.375, 15.75, 13, 0.030538591713165242, 0.9694614082868348",
+ "0.8125, 0.75, 0.25, 0.40160289105220459, 0.59839710894779541",
+ "0.8125, 0.75, 8.75, 0.99999978600605499, 2.1399394497066033e-07",
+ "0.8125, 0.75, 13, 0.99999999984153987, 1.584600819904228e-10",
+ "0.8125, 5.5, 0.25, 0.06149698416108057, 0.93850301583891937",
+ "0.8125, 5.5, 8.75, 0.99978985623633465, 0.00021014376366534313",
+ "0.8125, 5.5, 13, 0.99999931052227398, 6.8947772598025737e-07",
+ "0.8125, 15.75, 0.25, 0.003971077054507376, 0.99602892294549261",
+ "0.8125, 15.75, 8.75, 0.97219596511873307, 0.027804034881266888",
+ "0.8125, 15.75, 13, 0.9993055839090208, 0.00069441609097920377",
+ })
+ void testRegularizedBeta(double x, double a, double b, double ibeta, double ibetac) {
+ final double eps = 1e-13;
+ Assertions.assertEquals(ibeta, RegularizedBeta.value(x, a, b), ibeta * eps);
+ Assertions.assertEquals(ibetac, RegularizedBeta.complement(x, a, b), ibetac * eps);
+ }
+
+ /**
+ * Test the incomplete beta function uses the policy containing the epsilon and
+ * maximum iterations for series evaluations. The data targets each method computed
+ * using a series component to check the policy is not ignored.
+ *
+ * @see BoostBetaTest#testIBetaPolicy(double, double, double, double, double, double, double)
+ */
+ @ParameterizedTest
+ @CsvSource(value = {
+ // Target ibetaSeries
+ "4.201121919322759E-5,2.1881177963223308E-4,0.6323960423469543,23803.707603529016,4569.595256369948,0.838947362634037,0.16105263736596298",
+ "0.22512593865394592,0.4898320138454437,0.7996731996536255,4.764013990849158,0.9820281524501243,0.8290948573018541,0.17090514269814597",
+ "4.623167842510156E-5,4.340034502092749E-5,0.135563462972641,21628.337679706918,23043.143905699715,0.48416432390665337,0.5158356760933467",
+ "2.9415024982881732E-5,4.1924233664758503E-4,0.3082362115383148,33995.42298781772,2386.0630988783787,0.9344154580933705,0.0655845419066295",
+ "1.184685606858693E-5,0.015964560210704803,0.3082362115383148,84409.7658651171,63.42758275377908,0.9992491395179357,7.508604820642986E-4",
+ "3.529437162796967E-5,2.2326681573758833E-5,0.6323960423469543,28333.671885777032,44788.91695876111,0.3874817937042091,0.6125182062957909",
+ "0.06715317070484161,2.306319236755371,0.9084427952766418,13.787272071732604,0.001859217475218142,0.999865167903893,1.3483209610697333E-4",
+ "0.3183284401893616,3.165504217147827,0.07764927297830582,1.3374998709679642,0.6794195418585712,0.6631399660602084,0.3368600339397915",
+ "0.15403440594673157,4.049813747406006,0.34629878401756287,4.872033861103044,0.08561968850485947,0.9827297959310547,0.01727020406894529",
+ "1.3317101001739502,0.7650398015975952,0.6445860862731934,0.47144799136487586,0.5594135526519237,0.4573339592510717,0.5426660407489283",
+ "0.11902070045471191,7.269547462463379,0.19963125884532928,6.225047194692518,0.08420413075357451,0.9866538632858122,0.013346136714187858",
+ "2.664715051651001,0.6914005279541016,0.8443243503570557,0.3338388990912521,0.3587830340198169,0.4819929648946269,0.518007035105373",
+ "1.0562920570373535,6.812713623046875,0.8258343935012817,0.12732498812245932,9.807107058749213E-7,0.9999922976379849,7.702362015088557E-6",
+ "1.7118667364120483,3.0191311836242676,0.07594671100378036,0.0064151684204504875,0.10850933283432233,0.05582072012850319,0.9441792798714969",
+ // Target betaSmallBLargeASeries
+ "0.04634225741028786,0.34317314624786377,0.24176712334156036,20.363670894714268,3.6307737402387885,0.8486827348798152,0.15131726512018484",
+ "2.113992923113983E-5,1.7535277947899885E-5,0.8350250720977783,47305.46963235176,57026.27394012006,0.4534139659948463,0.5465860340051537",
+ "4.005068331025541E-5,42.84983825683594,0.12707412242889404,24964.03974453176,4.764518849491958E-4,0.9999999809144722,1.9085527852527327E-8",
+ "67.90167999267578,0.8324270844459534,0.9676981568336487,0.002669338211636283,0.031098353139101195,0.0790500654578503,0.9209499345421497",
+ "0.395370751619339,0.004023698624223471,0.9058013558387756,4.307485901473997,246.2348442100959,0.017192647244702382,0.9828073527552976",
+ "3.444607973098755,66.36054992675781,0.09654488414525986,1.4741027361579568E-6,8.307589573110104E-8,0.9466497330301025,0.05335026696989754",
+ "1.0665277242660522,4.985442161560059,0.2400285303592682,0.12523918055824373,0.0476465037156954,0.7244045745268357,0.2755954254731643",
+ // Target ibetaFraction2
+ "1.319732904434204,4.903014659881592,0.33251503109931946,0.0837704419861451,0.021604794441302123,0.7949727547593459,0.20502724524065405",
+ "485.7690734863281,190.16734313964844,0.6323960423469543,7.8023253024461E-182,7.885435919806278E-176,9.894592590329194E-7,0.9999990105407409",
+ })
+ void testIBetaPolicy(double a, double b, double x, double beta, double betac, double ibeta, double ibetac) {
+ // Low iterations should fail to converge
+ Assertions.assertThrows(ArithmeticException.class, () -> RegularizedBeta.value(x, a, b, EPS, 1), "ibeta");
+ Assertions.assertThrows(ArithmeticException.class, () -> RegularizedBeta.complement(x, a, b, EPS, 1), "ibetac");
+
+ // Low epsilon should not be as accurate
+
+ // Ignore 0 or 1
+ if ((int) ibeta != ibeta) {
+ final double p1 = RegularizedBeta.value(x, a, b);
+ final double p2 = RegularizedBeta.value(x, a, b, 1e-3, MAX_ITER);
+ assertCloser("ibeta", ibeta, p1, p2);
+ }
+ if ((int) ibetac != ibetac) {
+ final double q1 = RegularizedBeta.complement(x, a, b);
+ final double q2 = RegularizedBeta.complement(x, a, b, 1e-3, MAX_ITER);
+ assertCloser("ibetac", ibetac, q1, q2);
+ }
+ }
+
+ @ParameterizedTest
+ @CsvFileSource(resources = "ibeta_derivative_med_data.csv")
+ void testIBetaDerivative(double a, double b, double x, double expected) {
+ final double actual = RegularizedBeta.derivative(x, a, b);
+ TestUtils.assertEquals(expected, actual, 150);
+ }
+
+ /**
+ * Assert the regularized beta exactly matches the expected result.
+ * The complement is tested as 1-expected.
+ *
+ * @param x Argument x
+ * @param a Argument a
+ * @param b Argument b
+ * @param expected Expected result
+ */
+ private static void assertRegularizedBeta(double x, double a, double b, double expected) {
+ Assertions.assertEquals(expected, RegularizedBeta.value(x, a, b), 1e-15);
+ Assertions.assertEquals(1 - expected, RegularizedBeta.complement(x, a, b), 1e-15);
+ Assertions.assertEquals(expected, RegularizedBeta.value(x, a, b, EPS, MAX_ITER), 1e-15);
+ Assertions.assertEquals(1 - expected, RegularizedBeta.complement(x, a, b, EPS, MAX_ITER), 1e-15);
+ }
+
+ /**
+ * Assert x is closer to the expected result than y.
+ */
+ private static void assertCloser(String msg, double expected, double x, double y) {
+ final double dx = Math.abs(expected - x);
+ final double dy = Math.abs(expected - y);
+ Assertions.assertTrue(dx < dy,
+ () -> String.format("%s %s : %s (%s) : %s (%s)", msg, expected, x, dx, y, dy));
+ }
}