我正在使用以下 R 代码,取自已发表的论文(以下引文)。这是代码:
int2=function(x,r,n,p) {
(1+x)^((n-1-p)/2)*(1+(1-r^2)*x)^(-(n-1)/2)*x^(-3/2)*exp(-n/(2*x))}
integrate(f=int2,lower=0,upper=Inf,n=530,r=sqrt(.245),p=3, stop.on.error=FALSE)
当我运行它时,我收到错误“非有限函数值”。然而 Maple 能够将其计算为 4.046018765*10^27。
我尝试在包 pracma 中使用“integral”,这给了我一个不同的错误:
if (delta < tol) break 中的错误:缺少需要 TRUE/FALSE 的值
总体目标是计算两个积分的比率,如 Wetzels & Wagenmakers (2012)“相关性的默认贝叶斯假设检验”( http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3505519/ ) 中所述。整个函数如下:
jzs.pcorbf = function(r0, r1, p0, p1, n) {
int = function(r,n,p,g) {
(1+g)^((n-1-p)/2)*(1+(1-r^2)*g)^(-(n-1)/2)*g^(-3/2)*exp(-n/(2*g))};
bf10=integrate(int, lower=0,upper=Inf,r=r1,p=p1,n=n)$value/
integrate(int,lower=0,upper=Inf,r=r0,p=p0,n=n)$value;
return(bf10)
}
谢谢!
最佳答案
问题是您的积分函数在其域中使用 NaN
值调用时生成 x
值。您正在从 0 集成到 Infinity
,所以让我们检查一个有效的 x
值 1000:
int2(1000, sqrt(0.245), 530, 3)
# [1] NaN
您的目标乘以四部分:
x <- 1000
r <- sqrt(0.245)
n <- 530
p <- 3
(1+x)^((n-1-p)/2)
# [1] Inf
(1+(1-r^2)*x)^(-(n-1)/2)
# [1] 0
x^(-3/2)
# [1] 3.162278e-05
exp(-n/(2*x))
# [1] 0.7672059
我们现在可以看到问题在于您将无穷大乘以 0(或者说,数值上等于无穷大乘以数值上等于 0 的数值),这导致了数值问题。不是计算
a*b*c*d
,而是计算 exp(log(a) + log(b) + log(c) + log(d))
(使用 log(a*b*c*d) = log(a)+log(b)+log(c)+log(d)
的身份)会更稳定。另一个快速说明——值 x=0 需要一个特殊情况。int3 = function(x, r, n, p) {
loga <- ((n-1-p)/2) * log(1+x)
logb <- (-(n-1)/2) * log(1+(1-r^2)*x)
logc <- -3/2 * log(x)
logd <- -n/(2*x)
return(ifelse(x == 0, 0, exp(loga + logb + logc + logd)))
}
integrate(f=int3,lower=0,upper=Inf,n=530,r=sqrt(.245),p=3, stop.on.error=FALSE)
# 1.553185e+27 with absolute error < 2.6e+18
关于r - 在 R 中使用集成()时的 "non-finite function value",我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/28653672/