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)
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
