0

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

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

在此处输入图像描述

我尝试使用以下函数在 R 中编写此代码（在此处遵循类似的 SA 问题：Quadrature toapproximate a changed 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

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

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

4

1 回答 1

1

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

一般来说，使用专门的工具可以更好地解决这样的问题，或者使用桥接采样，即。如果您有 MCMC 输出，则通过bridgesampling 包；如果您有更一般的多变量集成问题，则通过 cubature 包。

但是，如果我们绝对必须使用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)
于 2018-03-12T12:54:23.980 回答