Skip to content

Instantly share code, notes, and snippets.

@gaberoo
Last active June 25, 2024 18:05
Show Gist options
  • Save gaberoo/4619102 to your computer and use it in GitHub Desktop.
Save gaberoo/4619102 to your computer and use it in GitHub Desktop.
R code to run an **MCMC** chain using a **Metropolis-Hastings** algorithm with a Gaussian proposal distribution. Although there are hundreds of these in various packages, none that I could find returned the likelihood values along with the samples from the posterior distribution. However, if you have these likelihood values, it's very easy to ca…
##############################################################################
# Metropolis-Hastings MCMC
#
# Runs a Metropolis-Hasting MCMC chain for a given likelihood function.
# Proposal steps are sampled from a Gaussian distribution, either in a single
# step or sequentially over the parameter space.
#
# Input:
# theta : starting value of the chain
# lik.fun : likelihood function
# prior.fun : prior probability function
# ... : parameters passed to lik.fun and prior.fun
# V : variance-covariance matrix for the proposal
# maxit : length of the chain
# sequential : TRUE or FALSE whether proposals are chosen for all parameters
# in a single step or sequentially.
#
# Output:
# list()
# chain : matrix with the posterior samples
# lik : likelihood values of those posterior samples
##############################################################################
metrop <- function(theta,lik.fun,prior.fun,
V=diag(theta),maxit=1000,thin=1,
sequential=FALSE,pb=FALSE,...) {
if (! sequential) {
if (! require(mvtnorm)) stop("You need to install the mvtnorm library.")
}
if (pb) { .pb <- txtProgressBar(0,maxit,style=3) }
ndim <- length(theta)
chain <- matrix(NA,nrow=maxit,ncol=ndim)
lik <- vector(length=maxit)
last <- theta
chain[1,] <- theta
last.lik <- lik.fun(theta,...)
lik[1] <- last.lik
last.prior <- prior.fun(theta,...)
it <- 1
naccept <- 0
if (is.finite(last.lik)) {
for (it in seq(2,maxit)) {
if (pb) setTxtProgressBar(.pb,it)
if (sequential) {
for (j in 1:thin) {
for (i in 1:length(theta)) {
if (V[i,i] > 0) {
accept <- FALSE
proposal <- last
proposal[i] <- rnorm(1,mean=last[i],sd=sqrt(V[i,i]))
proposal.prior <- prior.fun(proposal,...)
if (is.finite(proposal.prior)) {
proposal.lik <- lik.fun(proposal,...)
alpha <- exp(proposal.lik+proposal.prior-last.lik-last.prior)
if (alpha > runif(1)) accept <- TRUE
}
if (accept) {
last <- proposal
last.lik <- proposal.lik
last.prior <- proposal.prior
naccept <- naccept + 1/sum(diag(V) > 0)
}
}
}
}
chain[it,] <- last
lik[it] <- last.lik
} else {
for (j in 1:thin) {
accept <- FALSE
proposal <- rmvnorm(1,mean=last,sigma=V)
proposal.prior <- prior.fun(proposal,...)
if (is.finite(proposal.prior)) {
proposal.lik <- lik.fun(proposal,...)
alpha <- exp(proposal.lik+proposal.prior-last.lik-last.prior)
if (alpha > runif(1)) accept <- TRUE
}
if (accept) {
last <- proposal
last.lik <- proposal.lik
last.prior <- proposal.prior
naccept <- naccept+1
}
}
chain[it,] <- last
lik[it] <- last.lik
if (! pb) {
message(sprintf("Acceptance ratio = %.4f",naccept/(it*thin)))
}
}
}
}
if (pb) close(.pb)
message(sprintf("Acceptance ratio = %.4f",naccept/(maxit*thin)))
return(list(chain=chain,lik=lik))
}
##############################################################################
# Estimate Deviance Information Criterion (DIC)
#
# References:
# Bayesian Data Analysis.
# Gelman, A., Carlin, J., Stern, H., and Rubin D.
# Second Edition, 2003
#
# Bayesian predictive information criterion for the evaluation of
# hierarchical Bayesian and empirical Bayes models.
# Ando, T.
# Biometrika, 2007
#
# Input:
# x : matrix of posterior samples
# lik : vector of the likelihood of the posterior samples
# lik.fun : function that calculates the likelihood
# ... : other parameters that are passed to 'lik.fun'
#
# Output:
# list()
# DIC : Deviance Information Criterion
# IC : Bayesian Predictive Information Criterion
# pD : Effective number of parameters (pD = Dbar - Dhat)
# pV : Effective number of parameters (pV = var(D)/2)
# Dbar : Expected value of the deviance over the posterior
# Dhat : Deviance at the mean posterior estimate
##############################################################################
calc.dic <- function(x,lik,lik.fun,...) {
D.bar <- -2*mean(lik)
theta.bar <- summary(x)$statistics[,"Mean"]
D.hat <- -2*lik.fun(theta.bar,...)
pD <- D.bar - D.hat
pV <- var(-2*lik)/2
list(DIC=pD+D.bar,IC=2*pD+D.bar,pD=pD,pV=pV,Dbar=D.bar,Dhat=D.hat)
}
##############################################################################
# Estimate Marginal Likelihood
#
# Reference:
# Marginal likelihood from the Metropolis--Hastings output.
# Chib, S. and Jeliazkov, I.,
# Journal of the American Statistical Association, 2001
#
# Input:
# x : mcmc object or mcmc.list object,
# the samples from the posterior
# lik : vector, likelihood values for the samples in x
# lik.fun : function, calculates the likelihood
# prior.fun : function, calculates the prior probability
# num.samples : integer, number of samples used in the estimate
# log : boolean, if log-likelihood or likelihood is used
#
# Output:
# list()
# ln.m  : marginal liklihood (log if log-liklihood is used)
# ln.lik.star : reference likelihood (log if log-liklihood is used)
# The reference sample from the posterior is the
# sample that has the largest likelihood.
# ln.pi.star : refernce prior probability
# (log if log-likelihood is used)
# ln.pi.hat : posterior ordinate at reference sample
# (log if log-likelihood is used)
##############################################################################
marginal.likelihood <- function(x,lik,V,lik.fun,prior.fun,...,
num.samples=1000,log=TRUE) {
if (class(x) != "mcmc" & class(x) != "mcmc.list") {
stop("x must be an mcmc or mcmcList object.")
}
# get mean parameters
y <- summary(x)
theta.star <- y$statistics[,"Mean"]
lik.star <- lik.fun(theta.star,...)
# get samples from posterior
g <- sample(1:nrow(x),num.samples,replace=TRUE)
theta.g <- x[g,]
q.g <- dmvnorm(theta.g,mean=theta.star,sigma=V,log=FALSE)
#lik.g <- apply(theta.g,1,lik.fun,...)
lik.g <- lik[g]
alpha.g <- sapply(lik.g,function(l) min(1,exp(lik.star-l)))
# get samples from proposal
theta.j <- rmvnorm(num.samples,mean=theta.star,sigma=V)
lik.j <- apply(theta.j,1,lik.fun,...)
alpha.j <- sapply(lik.j,function(l) min(1,exp(l-lik.star)))
pi.hat <- mean(alpha.g*q.g)/mean(alpha.j)
pi.star <- 0
if (!is.null(prior.fun)) pi.star <- prior.fun(theta.star,...)
ln.m <- lik.star + pi.star - log(pi.hat)
if (! log) ln.m <- exp(ln.m)
list(ln.m=ln.m,ln.lik.star=lik.star,ln.pi.star=pi.star,ln.pi.hat=log(pi.hat))
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment