R:如何让gamma()返回实际数字而不是Inf

Question

运行gamma(200)返回R中的Inf，是否有可能让R以某种方式返回实际数字?看起来gamma(171.6)以上的任何内容都将在r中返回Inf

Answer 1

问题是你不能用双精度表示它:

gamma(200) # too large value
#R> [1] Inf
lgamma(200) # but log is not
#R> [1] 857.9337
exp(857) # the issue!
#R> [1] Inf
.Machine$double.xmax # maximum double value
#R> [1] 1.797693e+308
gamma(171) # almost there!
#R> [1] 7.257416e+306

您可以使用lgamma来代替gamma函数的对数。否则你将需要使用比R的浮点精度更高的第三方库。

谷歌搜索表明，Rmpfr::igamma函数可能是你想要的，如果你不能与伽马函数的对数工作:

Rmpfr::igamma(171, 0)
#R> 1 'mpfr' number of precision  53   bits 
#R> [1] 7.257415615307999e+306
Rmpfr::igamma(200, 0)
#R> 1 'mpfr' number of precision  53   bits 
#R> [1] 3.9432893368239526e+372

Answer 2

使用Benjamin Cristoffersen提出的lgamma，您可以将有效数和指数(以10为基数)作为单独变量计算:

(res <- gamma(100))
9.332622e+155
# Natural logarithm of result
(ln_res <- lgamma(100))
359.1342
# log base 10 of result
(log10_res <- ln_res/log(10))
155.97
# decimal part of the number above, raised to the 10th power
(significand_res <- 10 ^ (log10_res %% 1))
9.332622
# non-decimal part
(exp_res <- log10_res %/% 1)
155

对于gamma(200)，返回:3.9432 * 10 ^ 372