random walk Metropolis for a equal mixture of two normals

convergenceToADistribution.R

mn = 3
f = function(x) {
	n = length(x)
	apply(.5*cbind(dnorm(x,-mn),dnorm(x,mn)),1,sum)
}

# Random-walk Metropolis
set.seed(1)
n.reps = 1e5
x.reps = rep(NA,n.reps)
x.reps[1] = 0

sigma = 1
for (i in 2:n.reps) {
	proposal = rnorm(1,x.reps[i-1],sigma)
	log.mh = log(f(proposal))-log(f(x.reps[i-1]))
	accept = log(runif(1))<log.mh 
	
	x.reps[i] = ifelse(accept, proposal, x.reps[i-1])
	
	# Auto-tune MH proposal
	if (i<n.reps/2)	sigma = ifelse(accept, sigma*1.001, sigma/1.001)
	stopifnot(!is.na(x.reps[i]))
}
length(unique(x.reps))/n.reps # empirical acceptance rate

png("convergenceToDistribution.png")
par(mfrow=c(2,2))
curve(f,-2*mn, 2*mn, main="Posterior distribution")
hist(x.reps[seq(n.reps/2,n.reps)], 100, freq=F, main="Histogram approximation")
curve(f,col="red", add=T)
plot(x.reps,type="l", main="Traceplot")
plot(x.reps[1:4000], type="l", main="Partial traceplot")
dev.off()