Repository: commons-math Updated Branches: refs/heads/develop f695c9ce3 -> 1325484c3
Improved digamma function, especially for negative values. Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/53661585 Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/53661585 Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/53661585 Branch: refs/heads/develop Commit: 5366158572bd88f18df7883ab823ea62b7701d77 Parents: cbae75b Author: Eric Prescott-Gagnon <eric.prescottgag...@jda.com> Authored: Fri Jan 8 14:18:30 2016 -0500 Committer: Eric Prescott-Gagnon <eric.prescottgag...@jda.com> Committed: Thu May 5 13:48:13 2016 -0400 ---------------------------------------------------------------------- .../org/apache/commons/math4/special/Gamma.java | 42 ++++++++++++-------- .../apache/commons/math4/special/GammaTest.java | 1 + 2 files changed, 26 insertions(+), 17 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-math/blob/53661585/src/main/java/org/apache/commons/math4/special/Gamma.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/special/Gamma.java b/src/main/java/org/apache/commons/math4/special/Gamma.java index 70aaec3..c59049c 100644 --- a/src/main/java/org/apache/commons/math4/special/Gamma.java +++ b/src/main/java/org/apache/commons/math4/special/Gamma.java @@ -427,18 +427,16 @@ public class Gamma { * <p>Computes the digamma function of x.</p> * * <p>This is an independently written implementation of the algorithm described in - * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.</p> + * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976. + * A reflection formula (https://en.wikipedia.org/wiki/Digamma_function#Reflection_formula) + * is incorporated to improve performance on negative values.</p> * * <p>Some of the constants have been changed to increase accuracy at the moderate expense - * of run-time. The result should be accurate to within 10^-8 absolute tolerance for - * x >= 10^-5 and within 10^-8 relative tolerance for x > 0.</p> - * - * <p>Performance for large negative values of x will be quite expensive (proportional to - * |x|). Accuracy for negative values of x should be about 10^-8 absolute for results - * less than 10^5 and 10^-8 relative for results larger than that.</p> + * of run-time. The result should be accurate to within 10^-8 relative tolerance for + * 0 < x < 10^-5 and within 10^-8 absolute tolerance otherwise.</p> * * @param x Argument. - * @return digamma(x) to within 10-8 relative or absolute error whichever is smaller. + * @return digamma(x) to within 10-8 relative or absolute error whichever is larger. * @see <a href="http://en.wikipedia.org/wiki/Digamma_function";>Digamma</a> * @see <a href="http://www.uv.es/~bernardo/1976AppStatist.pdf";>Bernardo's original article </a> * @since 2.0 @@ -448,22 +446,32 @@ public class Gamma { return x; } + double digamma = 0.0; + if (x < 0) { + // use reflection formula to fall back into positive values + digamma -= FastMath.PI / FastMath.tan(FastMath.PI * x); + x = 1 - x; + } + if (x > 0 && x <= S_LIMIT) { // use method 5 from Bernardo AS103 // accurate to O(x) - return -GAMMA - 1 / x; + return digamma -GAMMA - 1 / x; } - if (x >= C_LIMIT) { - // use method 4 (accurate to O(1/x^8) - double inv = 1 / (x * x); - // 1 1 1 1 - // log(x) - --- - ------ + ------- - ------- - // 2 x 12 x^2 120 x^4 252 x^6 - return FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv / 252)); + while (x < C_LIMIT) { + digamma -= 1.0 / x; + x += 1.0; } - return digamma(x + 1) - 1 / x; + // use method 4 (accurate to O(1/x^8) + double inv = 1 / (x * x); + // 1 1 1 1 + // log(x) - --- - ------ + ------- - ------- + // 2 x 12 x^2 120 x^4 252 x^6 + digamma += FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv / 252)); + + return digamma; } /** http://git-wip-us.apache.org/repos/asf/commons-math/blob/53661585/src/test/java/org/apache/commons/math4/special/GammaTest.java ---------------------------------------------------------------------- diff --git a/src/test/java/org/apache/commons/math4/special/GammaTest.java b/src/test/java/org/apache/commons/math4/special/GammaTest.java index 51982a3..c320526 100644 --- a/src/test/java/org/apache/commons/math4/special/GammaTest.java +++ b/src/test/java/org/apache/commons/math4/special/GammaTest.java @@ -107,6 +107,7 @@ public class GammaTest { Assert.assertEquals(-100.56088545786867450, Gamma.digamma(0.01), eps); Assert.assertEquals(-4.0390398965921882955, Gamma.digamma(-0.8), eps); Assert.assertEquals(4.2003210041401844726, Gamma.digamma(-6.3), eps); + Assert.assertEquals(-3.110625123035E-5, Gamma.digamma(1.4616), eps); } @Test
