Example of Metropolis Algorithm

| No Comments

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

About this Entry

This page contains a single entry by quic0038 published on March 20, 2012 12:08 AM.

Code for Thursday's Lecture was the previous entry in this blog.

My Spatial Lecture! is the next entry in this blog.

Find recent content on the main index or look in the archives to find all content.

Categories

Pages

Powered by Movable Type 4.31-en