mMCMC: Run multiple MCMC chains.

In dereksonderegger/STA578: Functions for a statistical computing course

Description

This function generalizes the MCMC function to allow for running multiple chains and alows for multi-diminsionality in the distribution.

Usage

mMCMC(df, start, rprop, dprop = NULL, N = 1000, num.chains = 4)

Arguments

df	A function that is the target distribution
start	The starting point for the MCMC
rprop	A function that takes the current chain step and proposes a new value
dprop	The density function of the proposal distribution. If NULL, it is assumed to be symetric.
N	How many chain steps to calculate.
num.chains	How many chains to create.

Value

An object of type mcmc.list. This type of object can be used by functions in the coda package.

Examples

# Target Distribution
df <- function(x){
   return(.7*dnorm(x, mean=2, sd=1) + .3*dnorm(x, mean=5, sd=1))
}

# Proposal distribution
rprop <- function(x){ x + runif(1, -2, 2)}
dprop <- function(x.star, x){ dunif(x.star-x, -2, 2) }

# Create 4 chains
chains <- mMCMC(df, 4, rprop, dprop, N=1000, num.chains=4)
traceplot(chains)

##############################
##  A multi-variate example
##############################
dtarget <- function(x){
  dmvnorm(x, mean=c(3,10), sigma=matrix(c(3,3,3,7), nrow=2))
}

x1 <- seq(-6,12,length=101)
x2 <- seq(-11,31, length=101)
contour.data <- expand.grid(x1=x1, x2=x2)
contour.data$Z <- apply(contour.data, MARGIN=1, dtarget)

target.map <- ggplot(contour.data, aes(x=x1, y=x2)) +
  stat_contour(aes(z=Z))
target.map

rprop <- function(x){
  rmvnorm(1, mean=x, sigma=diag(c(1,1)))
}
start <- c(x1=0, x2=0)

chain <- mMCMC(df=dtarget, start, rprop=rprop,
               N=20, num.chains=1)
plot_2D_chains(target.map, chains)

chain <- mMCMC(df=dtarget, start, rprop=rprop,
               N=100, num.chains=1)
plot_2D_chains(target.map, chains)

chain <- mMCMC(df=dtarget, start, rprop=rprop,
               N=1000, num.chains=1)
plot_2D_chains(target.map, chains)
traceplot( chain )
plot(chain)

chains <- mMCMC(df=dtarget, start, rprop=rprop,
               N=1000, num.chains=4)
plot_2D_chains(target.map, chains)
plot(chains)

dereksonderegger/STA578 documentation built on May 15, 2019, 3:58 a.m.