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)

(in reverse alphabetical order...)

stacks_BUGS.txt

MACfrailty_BUGS.txt

MAC.txt

dugongsNL_BUGS.txt

dugongsBin_BUGS.txt

dugongs_BUGS.txt

crprot_BUGS.txt

Ex 1 - hv.r

Ex 2 - Pumps!.txt

Ex 3 - InterStim.txt

oh, and:

HW 2 is due Thursday and

HW 3 is "officially" assigned (due 2/21)

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

R Slides

grades_example.R

grades_example.txt

grades_example.csv

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

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