[Rd] Bad values obtained from digamma implementation (PR#13714)
ted.dunning at gmail.com
ted.dunning at gmail.com
Sun May 24 06:35:12 CEST 2009
I was writing and testing an implementation of digamma and discovered that
the values produced by R are very poor for some arguments.
For example:
> print(digamma(1e-15),20)
[1] -562949953421312.5
whereas the correct value to 20 places as computed by Mathematica (and my
own independent implementation) is
-1.0000000000000005772e15
This is about a factor of 2 different from the correct value.
Here is a reference implementation in Java of the algorithm by Bernardo as
published as AS103 that does not suffer this problem. I hope it helps.
/**
* Computes gamma related functions.
* <p/>
* This is an independently written implementation of the algorithm
described in
* <p/>
* Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied
Statistics, 1976
* <p/>
*
* I have changed some of the constants 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.
*
* 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. */
public class Gamma {
public static final double GAMMA = 0.577215664901532860606512090082;
private static final double C = 49;
private static final double S = 1e-5;
public static double digamma(double x) {
if (x > 0 && x <= S) {
// use method 5 from Bernardo AS103
// accurate to O(x)
return -GAMMA - 1 / x;
}
if (x >= C) {
// 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 Math.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 /
120 - inv / 252));
}
return digamma(x + 1) - 1 / x;
}
}
And here are some test cases:
public class TestGamma {
@Test
public void digammaLargeArgs() {
double eps = 1e-8;
Assert.assertEquals(4.6001618527380874002, Gamma.digamma(100), eps);
Assert.assertEquals(3.9019896734278921970, Gamma.digamma(50), eps);
Assert.assertEquals(2.9705239922421490509, Gamma.digamma(20), eps);
Assert.assertEquals(2.9958363947076465821, Gamma.digamma(20.5),
eps);
Assert.assertEquals(2.2622143570941481605, Gamma.digamma(10.1),
eps);
Assert.assertEquals(2.1168588189004379233, Gamma.digamma(8.8), eps);
Assert.assertEquals(1.8727843350984671394, Gamma.digamma(7), eps);
Assert.assertEquals(0.42278433509846713939, Gamma.digamma(2), eps);
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);
}
@Test
public void digammaSmallArgs() {
// values for negative powers of 10 from 1 to 30 as computed by
webMathematica with 20 digits
// see functions.wolfram.com
double[] expected = {-10.423754940411076795, -100.56088545786867450,
-1000.5755719318103005,
-10000.577051183514335, -100000.57719921568107,
-1.0000005772140199687e6, -1.0000000577215500408e7,
-1.0000000057721564845e8, -1.0000000005772156633e9,
-1.0000000000577215665e10, -1.0000000000057721566e11,
-1.0000000000005772157e12, -1.0000000000000577216e13,
-1.0000000000000057722e14, -1.0000000000000005772e15, -1e+16,
-1e+17, -1e+18, -1e+19, -1e+20, -1e+21, -1e+22, -1e+23,
-1e+24, -1e+25, -1e+26,
-1e+27, -1e+28, -1e+29, -1e+30};
for (double n = 1; n < 30; n++) {
System.out.printf("10^-n = %.10f\n", Math.pow(10, -n));
checkRelativeError(String.format("Test %.0f: ", n),
expected[(int) (n - 1)], Gamma.digamma(Math.pow(10.0, -n)), 1e-8);
}
private void checkRelativeError(String msg, double expected, double
actual, double tolerance) {
System.out.printf("%s = %.2f (%s)\n", msg, Math.abs(expected -
actual) / actual, Math.abs(expected - actual) > Math.abs(tolerance *
actual));
Assert.assertEquals(msg, expected, actual, Math.abs(tolerance *
actual));
}
}
Feel free to use this code in any way you like. I wrote it and granted it
to the Apache Software Foundation; I will grant it
to the R community under the same terms if useful.
My contact information is Ted Dunning, ted.dunning at gmail.com
--please do not edit the information below--
Version:
platform = i386-apple-darwin8.11.1
arch = i386
os = darwin8.11.1
system = i386, darwin8.11.1
status =
major = 2
minor = 9.0
year = 2009
month = 04
day = 17
svn rev = 48333
language = R
version.string = R version 2.9.0 (2009-04-17)
Locale:
en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8
Search Path:
.GlobalEnv, package:stats, package:graphics, package:grDevices,
package:utils, package:datasets, package:methods, Autoloads, package:base
[[alternative HTML version deleted]]
More information about the R-devel
mailing list