marginalLikelihood: Calcluated the marginal likelihood from a set of MCMC samples

Description

Calcluated the marginal likelihood from a set of MCMC samples

Usage

marginalLikelihood(sampler, numSamples = 1000, method = "Chib", ...)

Arguments

sampler

an object that implements the getSample function, i.e. a mcmc / smc Sampler (list)

numSamples

number of samples to use. How this works, and if it requires recalculating the likelihood, depends on the method

method

method to choose. Currently available are "Chib" (default), the harmonic mean "HM", and sampling from the prior "prior". See details

...

further arguments passed to getSample

Details

The function currently implements three ways to calculate the marginal likelihood. The recommended way is the method "Chib" (Chib and Jeliazkov, 2001). which is based on MCMC samples, but performs additional calculations. Despite being the current recommendation, note there are some numeric issues with this algorithm that may limit reliability for larger dimensions.

The harmonic mean approximation, is implemented only for comparison. Note that the method is numerically unrealiable and usually should not be used.

The third method is simply sampling from the prior. While in principle unbiased, it will only converge for a large number of samples, and is therefore numerically inefficient.

References

Chib, Siddhartha, and Ivan Jeliazkov. "Marginal likelihood from the Metropolis-Hastings output." Journal of the American Statistical Association 96.453 (2001): 270-281.

Examples

Run this code

# Harmonic mean works OK for a low-dim case with 

likelihood <- function(x) sum(dnorm(x, log = TRUE))
prior = createUniformPrior(lower = rep(-1,2), upper = rep(1,2))
bayesianSetup <- createBayesianSetup(likelihood = likelihood, prior = prior)
out = runMCMC(bayesianSetup = bayesianSetup, settings = list(iterations = 5000))

plot(out)


marginalLikelihood(out, numSamples = 500)[[1]]
marginalLikelihood(out, method = "HM", numSamples = 500)[[1]]
marginalLikelihood(out, method = "Prior", numSamples =  500)[[1]]

# True marginal likelihood (brute force approximation)

marginalLikelihood(out, method = "Prior", numSamples =  10000)[[1]]


# Harmonic mean goes totally wrong for higher dimendsions - wide prior. 
# Could also be a problem of numeric stability of the implementation

likelihood <- function(x) sum(dnorm(x, log = TRUE))
prior = createUniformPrior(lower = rep(-10,3), upper = rep(10,3))
bayesianSetup <- createBayesianSetup(likelihood = likelihood, prior = prior)
out = runMCMC(bayesianSetup = bayesianSetup, settings = list(iterations = 5000))

plot(out)

marginalLikelihood(out, numSamples = 500)[[1]]
marginalLikelihood(out, method = "HM", numSamples = 500)[[1]]
marginalLikelihood(out, method = "Prior", numSamples =  500)[[1]]

# True marginal likelihood (brute force approximation)

marginalLikelihood(out, method = "Prior", numSamples =  10000)[[1]]

Run the code above in your browser using DataLab

Last chance! 50% off unlimited learning