« Spatial Talk April 30 | Main | MT2 is Posted »

Hints for sampling algorithms

Some students have asked for help getting started with the sampling algorithm code, so here is a little "pseudo code" to get you started. You'll have to fill in the problem-specific computations where indicated. I haven't debugged this, it's just a starting point! :) Pay particular attention to the dimension of the parameter you're drawing and the function h(), which you can either write out directly in the met() function or write as a separate function that you use within met().

## Sample Metropolis-Hastings algorithm
met <- function(data,initials,NREP){
# Constants
n <- dim(data)[1]
p <- length(initials)
# Empty matrix to hold samples
chain <- matrix(NA,NREP+1,p)
# Initialize the chain
chain[1,] <- initials
for (i in 2:(NREP+1)){
# Parameters of the proposal
meen <-
var <-
# Draw a proposal value
thetastar <- rnorm(1,meen, var)
# Compute the rejection ratio
logr <- h(thetastar) + dnorm(chain[i-1,],meen,var) - h(chain[i-1,]) - dnorm(thetastar,meen,var)
# Decide whether to accept or reject
if (logr > 0){
chain[i,] <- thetastar
} elseif (logr < log(runif(1))) {
chain[i,] <- thetastar
} else {
chain[i,] <- chain[i-1,]
}
}
return(chain[-1,])
}