Computing the marginal likelihood of a Gaussian model in R with integrate()

Question

I am trying to compute the marginal likelihood of a Gaussian model in R. More precisely, I am trying to integrate the likelihood over both a Gaussian prior on mu and a Gaussian prior on sigma, with some observations yi.

In other words, I am trying to compute:

I tried to write this in R using the following function (following a similar SA question here: 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

}

Using this function, I can compute the marginal likelihood of the above model for a set of observations:

set.seed(666)

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

[1] 9.704133e-24

However, the results I obtain with this procedure are quite different from results I obtain with grid approximation, or using other packages/softwares.

My question is then: is it possible to do this double integration with integrate ? If it is, how would you do that ?

score 1 · Answer 1 · answered Mar 12 '18 at 12:54

integrate() only takes in univariate functions. That is, the function you put in must be one-dimensional.

In general, such a problem is better tackled using specialised tools, either using something bridgesampling, ie. through the bridgesampling package if you have MCMC output or the cubature package if you have more general multivariate integration problems.

However, if we absolutely must do this using integrate() twice, we can make this work, but some errors need to be taken out of the code, and . Something like the following would work, although numerically the result seems to be zero most of the time, which is why you would generally try to obtain the log-marginal likelihood.

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)

Computing the marginal likelihood of a Gaussian model in R with integrate()

1 Answers1