[R] Metropolis-Hastings in R

[Patrick Breheny]

This sounds like homework (in which case you should ask your instructor 
for help), but among the problems I see just glancing at it:

* runif(-0.5,0.5) should be runif(1,-0.5,0.5)
* there are a host of problems involved with the fact that you have set 
theta1 up as a two-column matrix, but repeatedly manipulate it as though 
its rows are scalars
* you're missing a negative in the normal likelihood.

-- 
Patrick Breheny
Assistant Professor
Department of Biostatistics
Department of Statistics
University of Kentucky


On 03/14/2012 01:21 PM, Michael Williams wrote:
> Hi all,
> I'm trying to write a MH algorithm in R for a standard normal distribution,
> I've been trying for a good week or so now with multiple attempts and have
> finally given up trying to do it on my own as I'm beginning to run out of
> time for this, would somebody please tell me what is wrong with my latest
> attempt:
>
>
> n=100
> mu=0
> sigma=1
> lik<-function(theta) exp(((theta-mu)^2)/2*sigma)
> alpha<-function(theta,phi) min(lik(phi)/lik(theta),1)
> theta1<-matrix(0,n,2)
> theta1[1,]<-mu
> for(i in 2:n){
> theta<-theta1[i-1,]
> phi<-theta+runif(-0.5,0.5)
> k<-rbinom(1,1,alpha(theta,phi))
> k1<-k1+k
> theta1[i,]<-theta+k*(phi-theta)
> }
> plot(theta1)
>
>