Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Improved digamma function, especially for negative values.

This commit is contained in:
Eric Prescott-Gagnon 2016-01-08 14:18:30 -05:00
parent cbae75b900
commit 5366158572
2 changed files with 26 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
src/main/java/org/apache/commons/math4/special/Gamma.java
View File

@ -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&apos;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;
}
/**

1
src/test/java/org/apache/commons/math4/special/GammaTest.java
View File

@ -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