1

我正在尝试为由集成定义的分布编写似然函数。我正在使用 integration() 函数，但是当我尝试在函数的其余部分中使用它时，我得到了错误：

“B(alpha + i, beta + 6 - i)/B(alpha, beta) 中的错误：二元运算符的非数字参数”

例如，积分值是“9.501501，绝对误差 < 0.00078”。我试过使用 trunc() 但这也无济于事。我对 R 比较陌生，所以有一个已知的解决方案吗？任何帮助，将不胜感激！

B <- function(a,b){ 
   integrand <- function(t){(t^(a-1))*((1-t)^(b-1))} 
   integrate(integrand,lower=0,upper=1) 
} 
betalik <- function(alpha,beta){ 
    likelihood <- 0 Z <- c(37,22,25,29,34,49) 
    for(i in 1:6) 
       likelihood <- likelihood + 
         Z[i]*log((B(alpha+i,beta+6-i))/B(alpha,beta)) 
    return(likelihood)

}

多里安，

4

2 回答 2

4

?integrate使用来自....的集成示例

integrate(dnorm, -1.96, 1.96)
# 0.9500042 with absolute error < 1e-11

这里到底发生了什么？好吧，集成函数创建了一个“集成”类 S3 对象，它基本上是一个包含多个字段的列表，然后应用print()到它，依次执行print.integrate()。该函数的输出实际上不是字符串“0.9500042 with absolute error < 1e-11”，它只是以这种方式显示。

要了解integrate()实际创建的 R 对象是什么，请执行以下操作

obj = integrate(dnorm, -1.96, 1.96)
str(obj)
# List of 5
# $ value       : num 0.95
# $ abs.error   : num 1.05e-11
# $ subdivisions: int 1
# $ message     : chr "OK"
# $ call        : language integrate(f = dnorm, lower = -1.96, upper = 1.96)
# - attr(*, "class")= chr "integrate"

因此，如果您只想要积分的值，那么您必须提取value由该函数创建的列表的字段，以便进行计算。例如

10*integrate(dnorm, -1.96, 1.96)$value
# [1] 9.500042
于 2012-05-24T17:53:11.743 回答
4

跟进上面的评论和答案......

这是您的原始函数，格式更好，集成了@Fhnuzoag 关于$value从integrate()结果中提取组件的评论：

B <- function(a,b){
    integrand <- function(t){(t^(a-1))*((1-t)^(b-1))}
    integrate(integrand,lower=0,upper=1)$value
}
betalik <- function(alpha,beta){
    likelihood <- 0
    Z <- c(37,22,25,29,34,49)
    for(i in 1:6) likelihood <- likelihood +
        Z[i]*log((B(alpha+i,beta+6-i))/B(alpha,beta))
    return(likelihood)
}

在这里，我们检查@dason 的评论，即您的B函数等同于 R 的内置beta函数（但 R 的函数肯定更快，而且很可能更准确）：

all.equal(B(1.1,2.7),beta(1.1,2.7))  ## TRUE

我更喜欢单独指定“数据”：

Z <- c(37,22,25,29,34,49)

一个新版本的似然函数，它使用内置lbeta（log-beta）函数，并被矢量化：

blik2 <- function(alpha,beta,Z) {
    index <- 1:6
    sum(Z*(lbeta(alpha+index,beta+6-index)-lbeta(alpha,beta)))
}

all.equal(blik2(1.1,2.7,Z),betalik(1.1,2.7))
## small difference, blik2 is *probably* more
##   accurate ... "Mean relative difference: 6.406495e-08"

（进一步概括这一点并将6代码中的值替换为length(Z)... 可能会更好）

于 2012-05-24T18:02:03.823 回答