## 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)

## Code used in Tuesday's R/WinBUGS lecture

oh, and:
HW 2 is due Thursday and
HW 3 is "officially" assigned (due 2/21)

## R Programming: Files from Prof. Koopmeiners' Lecture

In case you're interested, here are all of the necessary files from Prof. Koopmeiners' lecture on R from Jan. 24th

If you have any questions with this program, or R programming in general, don't hesitate to email me

## First Blog Entry for 2012!

I'll occasionally write things on here like course reminders (like midterm deadlines), updates (regarding changes in office hours), and some homework help.

I'm tentatively planning on having the following office hours:
Tuesday 11-1
Thursday 12-2.

I also check my email pretty frequently, so feel free to ask me any questions you may have.

Harrison

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