Below is an example for coding up the Metropolis algorithm in R. Here,

1) the parameter of interest is "b",

2) "rq()" is generating a random value from q, my candidate density, which I'm saying is just a N(b,1), and

3) "h()" is the unnormalized posterior, which for me is just N(0,4)

Also, you don't have to make a function like my "met_b", but it cleans things up a bit.

Obviously, when you're doing the homework, h() may depend on other parameters and will certainly depend on the data, y.

#################functions...##########################

met_b=function(b){

#given the current position of the chain, b, this will generate the next

bs = rq(1, b) #sample 1 value from q( |b)

r = h(bs)/h(b) #you may need to use some algebra to simplify this expression,

#especially in multivariate settings where h(x) can be super small

if(r >= 1){

return(bs)

}else if(r >= runif(1)){

return(bs)

}else{

return(b)

}

}

rq=function(num,x){

return(rnorm(num,mean=x,sd=1))

}

h=function(x){

return(dnorm(x,sd=2))

}

###################main part...########################

T=1000 #the number of iterations you want... start small, and go from there

b=array(dim=c(T)) #makes b a vector of length T.

# If b is a k-vector, you can replace "c(T)" with "c(T,k)"

b[1]=0 #initialize b

for(it in 2:T){

b[it] = met_b(b[it-1])

}

par(mfrow=c(1,2))

plot(y=b,x=1:T,type='l',main='History')

hist(b,main='Density',freq=FALSE)

## Leave a comment