Markov Chain Monte Carlo

mcmc.R

## The observed data. 
X <- rnorm(1000, 2, 5)

loglik <- function(mu, sigma, X){
  ## The likelihood function
  sum( dnorm(X, mean=mu, sd=sigma, log=TRUE) )
}

## initial starting point 
pars <- c(mu=5, sigma=15)

## An important parameter
step_size <- 1
step_fn <- function(pars){
  ## A function to propose some new parameters
  if( runif(1) > 0.5){
    pars["mu"] <- rnorm(1, pars["mu"], step_size)
  } else {
    pars["sigma"] <- rnorm(1, pars["sigma"], 1)
  }
  pars
}

Q <- function(pars,pars_prime){ 
  ## The proposal distribution used in step_fn (Normal in our case)
  if(pars["mu"] == pars_prime["mu"]){ #if sigma changed
    ## Note we'll always work with logs, to avoid rounding errors
    dnorm(pars["sigma"], pars_prime["sigma"], step_size, log=TRUE)
  } else { #mu changed, so proposal density reflects mu
    dnorm(pars["mu"], pars_prime["mu"], step_size, log=TRUE)
  }
}

## Time over which to MCMC
MaxTime <- 10000
## Record this stuff 
history <- list(mu = numeric(MaxTime), sigma=numeric(MaxTime), lik=numeric(MaxTime))

## Ready to rock and roll
for(i in 1:MaxTime){
  ## Record the current position
  history$mu[i] = pars["mu"]
  history$sigma[i] = pars["sigma"]

  ## Calculate its likelihood
  L <- loglik(pars["mu"], pars["sigma"], X)
  ## oh, we'll store that too
  history$lik[i] <- L # log likelihood 

  ## Propose a new position
  proposed_pars <- step_fn(pars)
  ## Calculate its likelihood
  proposed_L <- loglik(proposed_pars["mu"], proposed_pars["sigma"], X)

  ## Calculate alpha, the acceptance probability
  ## Note that with logs, we use sums not products.  then exponentiate to get alpha
  alpha <- exp(proposed_L + Q(pars, proposed_pars) - L - Q(proposed_pars, pars) )

  ## If Q is symmetric, we just need likelihood ratio
#  alpha <- exp(proposed_L - L)

  ## Accept always if proposed improves likelihood
  if (alpha > 1){
    pars <- proposed_pars
#   print(paste("pars =", round(pars["mu"],3), "alpha=", round(alpha,3)))

  ## Accept with probability alpha if proposal decreases likelihood
  } else if(runif(1) <= alpha){
    pars <- proposed_pars
#   print(paste("accepted anyway, pars =", round(pars["mu"],9), 
#         round(proposed_pars["mu"],9), "alpha=", round(alpha,3)))
 } else {
    # nothing happens
 }

}

# Parameter posterior distribution reflects time spent in each area
# disgards burn-in
burnin <- 100

par(mfrow=c(1,3))
hist(history$mu[-burnin], main="Posterior for mu")
hist(history$sigma[-burnin], main="Posterior for sigma")

# look at the likelihood over time -- fuzzy catipillar?  
plot(history$lik, main="likelihood", ylab="loglik", xlab="iteration")