在R中计算伯努利似然函数积分时遇长度错误的问题

错误原因 #

integrate做数值积分时，会一次性传入**多个theta值（向量形式）**来计算，你的匿名函数没处理这种场景：

• 当theta是向量时，dbinom(x=c(1,0,1), prob=theta)会返回一个3行N列的矩阵（N是theta的长度），每列对应一个theta对应的三个伯努利概率。
• 而prod()默认会把整个矩阵的所有元素扁平化后求积，最终只返回1个值，这和输入的theta向量长度N不匹配，所以触发“结果长度错误”。

正确的返回长度 #

当传入长度为N的theta向量时，函数必须返回长度为N的向量，每个元素对应单个theta值的似然计算结果。

修正代码 #

方法1：直接写似然表达式（最简单高效） #

三次试验是2次成功1次失败，似然函数可直接写成theta^2 * (1-theta)：

integrate(
  function(theta) theta^2 * (1 - theta),
  lower = 0,
  upper = 1
)

方法2：修改原代码支持向量化输入 #

用apply按列求积，每列对应一个theta的似然值：

integrate(
  function(theta) {
    probs <- dbinom(x = c(1,0,1), size = 1, prob = theta, log = F)
    apply(probs, 2, prod)
  },
  lower = 0,
  upper = 1
)

或者用matrixStats包的rowProds（需转置矩阵，因为dbinom返回的行对应x、列对应theta）：

library(matrixStats)
integrate(
  function(theta) {
    probs <- dbinom(x = c(1,0,1), size = 1, prob = theta, log = F)
    rowProds(t(probs))
  },
  lower = 0,
  upper = 1
)

运行修正后的代码会得到正确结果：积分值为0.08333333（即1/12），和解析解一致。

内容的提问来源于stack exchange，提问作者Aku-Ville Lehtimäki