Example of Metropolis Algorithm

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)