public static double digamma(double x){
  assert x > 0 : x;
  double r, f, t;
  r=0;
  while (x <= 5) {
    r-=1 / x;
    x+=1;
  }
  f=1 / (x * x);
  t=f * (-1 / 12.0 + f * (1 / 120.0 + f * (-1 / 252.0 + f * (1 / 240.0 + f * (-1 / 132.0 + f * (691 / 32760.0 + f * (-1 / 12.0 + f * 3617 / 8160.0)))))));
  return r + Math.log(x) - 0.5 / x + t;
}

public static double digamma ( double x ) { assert x > 0 : x ; double r , f , t ; r = 0 ; while ( x <= 5 ) { r -= 1 / x ; x += 1 ; } f = 1 / ( x * x ) ; t = f * ( - 1 / 12.0 + f * ( 1 / 120.0 + f * ( - 1 / 252.0 + f * ( 1 / 240.0 + f * ( - 1 / 132.0 + f * ( 691 / 32760.0 + f * ( - 1 / 12.0 + f * 3617 / 8160.0 ) ) ) ) ) ) ) ; return r + Math . log ( x ) - 0.5 / x + t ; }
stolen from radford neal ' s fbm package . digamma ( x ) is defined as ( d / dx ) log gamma ( x ) . it is computed here using an asymptotic expansion when x > 5 . for x < = 5 , the recurrence relation digamma ( x ) = digamma ( x + 1 ) - 1 / x is used repeatedly . see venables & ripley , modern applied statistics with s - plus , pp . 151 - 152 . compute the digamma function . returns - inf if the argument is an integer less than or equal to zero .
