[Rd] Bad values obtained from digamma implementation (PR#13714)

[ted.dunning at gmail.com]

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