Gamma function returns unstable value?

Question

Gamma function should not take any negative value as an argument. Look at the code below where strange thing happens. Is this some problem with R?

I was using function optim to optimize some function containing:

gamma(sum(alpha))

with respect to alpha. R returns negative alpha.

> gamma(sum(alpha))
[1] 3.753+14
>sum(alpha)
[1] -3
gamma(-3)
[1] NaN
Warning message:
In gamma(-3) NaN's produced.

Can somebody explain? Or any suggestion for the optimization?

Thanks!

Answer 1

Gamma function is "not defined" at negative integer argument values so R returns Not a Number (NaN). The reason of the "strange" behaviour is decimal representation of numbers in R. In case the number differs from the nearest integer not very much, R rounds it during printing (in fact when you type alpha, R is calling for print(alpha). Please see the examples of such a behaviour below.

gamma(-3)
# [1] NaN
# Warning message:
#  In gamma(-3) : NaNs produced
x <- -c(1, 2, 3) / 2 - 1e-15
x
# [1] -0.5 -1.0 -1.5
sum(x)
# [1] -3
gamma(sum(x))
# [1] 5.361428e+13
curve(gamma, xlim = c(-3.5, -2.5)) 

Please see a graph below which explains the behaviour of gamma-function near negative integers: