Search code examples
restimationmcmc

Needing advice on MCMC estimation in R


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?


Solution

  • There are a few problems with the code.

    1. You calculate the log-likelihood, but ignore the log-prior. So effectively you're using a uniform prior, not the one you specified.

    2. 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.)

    3. 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)