Numerical precision problems in R?

Question

I have a problem with the following function in R:

test <- function(alpha, beta, n){
  result <- exp(lgamma(alpha) + lgamma(n + beta) - lgamma(alpha + beta + n) - (lgamma(alpha) + lgamma(beta) - lgamma(alpha + beta)))
  return(result)
}

Now if you insert the following values:

betabinom(-0.03292708, -0.3336882, 10)

It should fail and result in a NaN. That is because if we implement the exact function in Excel, we would get a result that is not a number. The implementation in Excel is simple, for J32 is a cell for alpha, K32 beta and L32 for N. The implementation of the resulting cell is given below:

=EXP(GAMMALN(J32)+GAMMALN(L32+K32)-GAMMALN(J32+K32+L32)-(GAMMALN(J32)+GAMMALN(K32)-GAMMALN(J32+K32)))

So this seems to give the correct answer, because the function is only defined for alpha and beta greater than zero and n greater or equal to zero. Therefore I am wondering what is happening here? I have also tried the package Rmpf to increase the numerical accuracy, but that does not seem to do anything.

Thanks

Answer 1

tl;dr log(gamma(x)) is defined more generally than you think, or than Excel thinks. If you want your function not to accept negative values of alpha and beta, or to return NaN, just test manually and return the appropriate values (if (alpha<0 || beta<0) return(NaN)).

It's not a numerical accuracy problem, it's a definition issue. The Gamma function is defined for negative real values: ?lgamma says:

The gamma function is defined by (Abramowitz and Stegun section 6.1.1, page 255)

Gamma(x) = integral_0^Inf t^(x-1) exp(-t) dt

for all real ‘x’ except zero and negative integers (when ‘NaN’ is returned).

Furthermore, referring to lgamma ...

... and the natural logarithm of the absolute value of the gamma function ...

(emphasis in original)

curve(lgamma(x),-1,1)

gamma(-0.1)          ## -10.68629
log(gamma(-0.1)+0i)  ## 2.368961+3.141593i
log(abs(gamma(-0.1)) ## 2.368961
lgamma(-0.1)         ## 2.368961

Wolfram Alpha agrees with second calculation.