计算 R 中伯努利似然函数的积分

Question

我正在尝试在 R 中集成一个非常简单的似然函数（三个伯努利试验），但似乎我还没有理解集成函数是如何工作的。我的代码如下：

integrate(
    function(theta) prod(
        dbinom(
            x = c(1,0,1), #x is the data: "success, fail, success"
            size = 1, #Just one Bernoulli trial times the length of x.
            prob = theta, #Since, this is a likelihood function, the parameter value is the variable
            log = F #We use a likelihood instead of log-likelihood
        )
    ),
    lower = 0, #The parameter theta cannot be smaller than zero...
    upper = 1 #...nor greater than one.
)

但是，我收到此错误消息：

Error in integrate(function(theta) prod(dbinom(x = c(1, 0, 1), size = 1,  : 
  evaluation of function gave a result of wrong length

现在，为什么结果长度不对呢？结果的长度始终为 1，因为匿名函数使用 prod 函数，该函数又创建函数 dbinom 返回的所有向量元素的乘积（该向量的长度为 3，因为其第一个参数 x 的长度为 3） .

那么结果应该是什么才具有正确的长度？

Answer 1

问题代码有两点错误：

• 你必须优化，而不是集成；
• (0, 1) 中许多因素的乘积将接近于零，这就是为什么您需要对对数似然求和。

fn <- function(theta, x) sum(dbinom(x, size = 1, prob = theta, log = TRUE))

MLE <- function(x) {
  f <- function(theta, x) fn(theta, x)
  optimise(f, c(0, 1), x = x, maximum = TRUE)$maximum
}

# the question's sample
MLE(c(1, 0, 1))
#> [1] 0.666675

# other samples, increasingly larger
set.seed(2023)
x <- rbinom(10, 1, p = 2/3)
MLE(x)
#> [1] 0.9999339

x2 <- rbinom(100, 1, p = 2/3)
MLE(x2)
#> [1] 0.5600005

x3 <- rbinom(1000, 1, p = 2/3)
MLE(x3)
#> [1] 0.691006

x4 <- rbinom(10000, 1, p = 2/3)
MLE(x4)
#> [1] 0.6746081

^{创建于 2023-07-31，使用 reprex v2.0.2}

---

编辑

我误解了这个问题，我以为它是在给定样本的情况下要求 MLE

c(1, 0, 1)

。
这是一个更正的

sum/log

版本，无非是Avraham的答案中的解决方案。结果和绝对误差是一样的。

fn <- function(theta) sum(dbinom(c(1, 0, 1), size = 1, prob = theta, log = TRUE)) |> exp()
integrate(Vectorize(fn), 0, 1)
#> 0.08333333 with absolute error < 9.3e-16

^{创建于 2023-08-01，使用 reprex v2.0.2}

Answer 2

这里的问题不是

dbinom

，而是

prod

。

dbinom

是 一个向量化函数，但是

prod

根据其定义，返回“一个数字（类型为“double”）或长度为 1 的复数向量。”

如评论中所述，最简单的方法是将函数简单地包装在

Vectorize

调用中的

integrate

中。例如：

fn <- function(theta) prod(dbinom(c(1, 0, 1), size = 1, prob = theta, log = FALSE))
integrate(Vectorize(fn), 0, 1)
0.08333333 with absolute error < 9.3e-16