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/08 09:23:36 UTC
[commons-numbers] branch master updated: NUMBERS-182: Avoid overflow in LogBeta for tiny arguments
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
The following commit(s) were added to refs/heads/master by this push:
new 23897b0 NUMBERS-182: Avoid overflow in LogBeta for tiny arguments
23897b0 is described below
commit 23897b03807b609fb757ba675c751238cb4a4391
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Sat Jan 8 09:23:32 2022 +0000
NUMBERS-182: Avoid overflow in LogBeta for tiny arguments
---
.../org/apache/commons/numbers/gamma/LogBeta.java | 8 ++++++--
.../apache/commons/numbers/gamma/LogBetaTest.java | 23 ++++++++++++++++++++++
src/changes/changes.xml | 3 +++
3 files changed, 32 insertions(+), 2 deletions(-)
diff --git a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/LogBeta.java b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/LogBeta.java
index 6decedb..6f08a9f 100644
--- a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/LogBeta.java
+++ b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/LogBeta.java
@@ -242,8 +242,12 @@ public final class LogBeta {
// The original NSWC implementation was
// LogGamma.value(a) + (LogGamma.value(b) - LogGamma.value(a + b));
// but the following command turned out to be more accurate.
- return Math.log(Gamma.value(a) * Gamma.value(b) /
- Gamma.value(a + b));
+ // Note: Check for overflow that occurs if a and/or b are tiny.
+ final double beta = Gamma.value(a) * Gamma.value(b) / Gamma.value(a + b);
+ if (Double.isFinite(beta)) {
+ return Math.log(beta);
+ }
+ return LogGamma.value(a) + (LogGamma.value(b) - LogGamma.value(a + b));
}
}
diff --git a/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/LogBetaTest.java b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/LogBetaTest.java
index 114d055..d92fba7 100644
--- a/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/LogBetaTest.java
+++ b/commons-numbers-gamma/src/test/java/org/apache/commons/numbers/gamma/LogBetaTest.java
@@ -21,6 +21,8 @@ import java.lang.reflect.Method;
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;
/**
* Tests for {@link LogBeta}.
@@ -699,4 +701,25 @@ class LogBetaTest {
() -> sumDeltaMinusDeltaSum(10, 9)
);
}
+
+ /**
+ * Hit edge case for log(beta(a,b)) = log(gamma(a)*gamma(b)/gamma(a+b))
+ * where gamma(a)*gamma(b) will overflow.
+ *
+ * <p>See NUMBERS-182.
+ */
+ @ParameterizedTest
+ @CsvSource({
+ // Computed using matlab betaln(a, b)
+ "1e-155, 1.5e-155, 357.41151503784306",
+ // MIN NORMAL, MIN_NORMAL * 3/4
+ "2.2250738585072014E-308, 1.668805393880401E-308, 709.2437163926513",
+ // MIN_VALUE, MIN_VALUE * 2
+ "4.9E-324, 1.0E-323, 744.84553702948938",
+ // MIN_VALUE, MIN_VALUE
+ "4.9E-324, 4.9E-324, 745.13321910194111",
+ })
+ void testTinyAB(double a, double b, double expected) {
+ Assertions.assertEquals(expected, LogBeta.value(a, b));
+ }
}
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index fac324f..065e42a 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -72,6 +72,9 @@ N.B. the Performance testing module requires Java 9+.
(The unit tests require Java 8+)
">
+ <action dev="aherbert" type="fix" issue="NUMBERS-182">
+ "LogBeta": Avoid overflow for tiny arguments.
+ </action>
<action dev="aherbert" type="add" issue="NUMBERS-180">
"GammaRatio": Compute the ratio of two gamma functions.
</action>