This project has moved and is read-only. For the latest updates, please go here.

Bug in DiGamma?

Sep 5, 2010 at 12:25 PM

double x = 3;

MathNet.Numerics.SpecialFunctions.DiGamma(x/2 + 1/2);
Should be: 1 - gamma = 0.42278...; Is: 0.0364899...

MathNet.Numerics.SpecialFunctions.DiGamma(x/2);
Should be: 2 - gamma - 2 ln(2) = 0.0364899... Is: 0.0364899...

Regards

Sep 5, 2010 at 2:01 PM

Hi Peter,

Octave and Mathematica give digamma(2) = 0.4227843 as does Math.NET Numerics.

Perhaps I'm not following the problem.  

Regards,

Marcus

Sep 5, 2010 at 2:39 PM

 If I call
 |  double x = 3;
 |  Console.WriteLine(MathNet.Numerics.SpecialFunctions.DiGamma(x / 2 + 1 / 2)
 |  + " " + MathNet.Numerics.SpecialFunctions.DiGamma(x / 2));
 
 with the code below, which I have downloaded yesterday from your site, I get:
 
 .036489974  .036489974
 
 Wolfram Alpha gives
 digamma(2)    = 0.422784335..
 digamma(3/2) = 0.036489973..

This is the same what I said in my first message.

/// <summary>
/// Computes the Digamma function which is mathematically defined as the derivative of the logarithm of the gamma function.
/// This implementation is based on
///     Jose Bernardo
///     Algorithm AS 103:
///     Psi ( Digamma ) Function,
///     Applied Statistics,
///     Volume 25, Number 3, 1976, pages 315-317.
/// Using the modifications as in Tom Minka's lightspeed toolbox.
/// </summary>
/// <param name="x">The argument of the digamma function.</param>
/// <returns>The value of the DiGamma function at <paramref name="x"/>.</returns>
public static double DiGamma(double x)
{
    const double C = 12.0;
    const double D1 = -0.57721566490153286;
    const double D2 = 1.6449340668482264365;
    const double S = 1e-6;
    const double S3 = 1.0 / 12.0;
    const double S4 = 1.0 / 120.0;
    const double S5 = 1.0 / 252.0;
    const double S6 = 1.0 / 240.0;
    const double S7 = 1.0 / 132.0;

    if (Double.IsNegativeInfinity(x) || Double.IsNaN(x))
    {
        return Double.NaN;
    }

    // Handle special cases.
    if (x <= 0 && Math.Floor(x) == x)
    {
        return Double.NegativeInfinity;
    }

    // Use inversion formula for negative numbers.
    if (x < 0)
    {
        return DiGamma(1.0 - x) + (Math.PI / Math.Tan(-Math.PI * x));
    }

    if (x <= S)
    {
        return D1 - (1 / x) + (D2 * x);
    }

    double result = 0;
    while (x < C)
    {
        result -= 1 / x;
        x++;
    }

    if (x >= C)
    {
        var r = 1 / x;
        result += Math.Log(x) - (0.5 * r);
        r *= r;

        result -= r * (S3 - (r * (S4 - (r * (S5 - (r * (S6 - (r * S7))))))));
    }

    return result;
}

Sep 5, 2010 at 2:54 PM

OK. The problem is that  1/2 evaluates to zero since they are both integers - so you are computing digamma(1.5). Change it to 1.0/2.0.

Sep 5, 2010 at 4:35 PM

Thanks cuda!