Suppose the random variable X ∼ N(μ,σ2) distribution
and I am given the observations x =(0.7669, 1.6709, 0.8944, 1.0321, 0.0793, 0.1033, 1.2709, 0.7798, 0.6483, 0.3256)
Given prior distribution μ ∼ N(0, (100)^2) and σ2 ∼ Inverse − Gamma(5/2, 10/2).
I am trying to give a MCMC estimation of the parameter.
Here is my sample code, I am trying to use Random Walk MCMC to do the estimation, but it seems that it does not converge:
x = c(0.7669,1.6709,0.8944, 1.0321, 0.0793, 0.1033, 1.2709, 0.7798, 0.6483, 0.3256)
mu.old = rnorm(1,0,100); #initial guess for mu
alpha = 5/2;beta = 10/2
sigma.old = sqrt(rinvgamma(1, shape = alpha, scale = 1/beta)) #initial guess for sigma
burn.in = 50000
rep = 100000
acc = 0
lik.old = sum(dnorm(x,mu.old,sd = sigma.old),log = T)
res = c(mu.old,sigma.old,lik.old)
for (i in (1:rep)){
mu.new = mu.old + rnorm(1,0,10)
sigma.new = sigma.old + rnorm(1,0,10)
if ((mu.new <=0)||sigma.new <= 0){
res = rbind(res,c(mu.old,sigma.old,lik.old))
}
else{
lik.new = sum(dnorm(x,mu.new,sigma.new,log = T))
r = exp(lik.new - lik.old)
u = runif(1,0,1)
if (u<r){
res = rbind(res,c(mu.new,sigma.new,lik.new))
mu.old = mu.new
sigma.old = sigma.new
lik.old = lik.new
acc = acc + 1
}
else{
res = rbind(res,c(mu.old,sigma.old,lik.old))
}
}
if (i %% 10000 == 0){
plot(res[,1],type = 'l')
}
}
Can anybody gives advice on how to find a better MCMC methods?
There are a few problems with the code.
You calculate the log-likelihood, but ignore the log-prior. So effectively you're using a uniform prior, not the one you specified.
The initial calculation of the log likelihood has parentheses in the wrong place:
sum(dnorm(x,mu.old,sd = sigma.old),log = T)
The last parameter should be log = TRUE
, inside the dnorm()
call, not outside it. (Avoid using T
for TRUE
. It's not worth the saving of 3 chars when typing.)
Using rbind()
to grow your res
array is very slow. It's better to allocate the whole thing in advance, and assign results to specific rows, e.g.
# initialize
res <- matrix(NA, nrow = reps, ncol = 3)
# in the loop
res[i,] <- c(mu.old,sigma.old,lik.old)