用integrate()计算R中高斯模型的边际似然答案

【问题标题】：Computing the marginal likelihood of a Gaussian model in R with integrate()用integrate()计算R中高斯模型的边际似然
【发布时间】：2018-08-19 13:04:54
【问题描述】：

我正在尝试计算 R 中高斯模型的边际似然。更准确地说，我正在尝试将 mu 上的高斯先验和 sigma 上的高斯先验的似然性与一些观察值 yi 相结合。

换句话说，我正在尝试计算：

我尝试使用以下函数在 R 中编写此代码（在此处遵循类似的 SA 问题：Quadrature to approximate a transformed beta distribution in R）：

marglik <- function(data) {

    integrand <-
        Vectorize(function(data, mu, sigma) {
            prod(dnorm(data, mu, sigma) ) * dnorm(mu, 110, 1) * dnorm(sigma, 10, 1)
            } )

    integrate(integrand, lower = 0, upper = Inf, mu = 100, sigma = 10)$value

}

使用这个函数，我可以计算上述模型对于一组观察值的边际似然：

set.seed(666)

d <- rnorm(100, mean = 107.5, sd = 2.5)
marglik(data = d)

[1] 9.704133e-24

但是，我通过此过程获得的结果与我通过网格近似或使用其他软件包/软件获得的结果大不相同。

那么我的问题是：是否可以使用集成进行这种双重集成？如果是，你会怎么做？

【问题讨论】：

标签： r bayesian numerical-integration integrate log-likelihood

【解决方案1】：

integrate() 只接受单变量函数。也就是说，你输入的函数必须是一维的。

一般来说，这样的问题最好使用专门的工具来解决，或者使用桥采样，即。通过bridgesampling package（如果您有 MCMC 输出）或cubature package（如果您有更一般的多元集成问题）。

但是，如果我们绝对必须使用integrate() 两次，我们可以完成这项工作，但需要从代码中删除一些错误，并且 .像下面这样的东西会起作用，尽管大多数时候结果似乎为零，这就是为什么你通常会尝试获得对数边际似然。

marglik <- function(data) {

  # Function that integrates over mu for given sigma.
  mu_integrand <- Vectorize(function(sigma) {
    mu_given_sigma_fun <- Vectorize(function(mu) {
      prod(dnorm(data, mu, sigma) ) * dnorm(mu, 110, 1) * dnorm(sigma, 10, 1)
      })
    integrate(mu_given_sigma_fun, lower = -Inf, upper = Inf)$value
  })

  integrate(mu_integrand, lower = 0, upper = Inf)$value

}

set.seed(666)

d <- rnorm(100, mean = 110, sd = 10)
marglik(data = d)

【讨论】：