Metropolis–Hastings algorithm variance isn't converging in R?

Question

I'm trying to simulate a sample from a t distribution with 4 degrees of freedom. The candidate density I'm using is a normal(0,1) distribution. Although the mean does converge to 0, the variance keeps converging to 1.5 instead of 2 as a t distribution with 4 degrees of freedom should have. It's especially confusing since the histogram matches the desired t distribution density almost perfectly. The code I am using is this:

normMH <- function(n=100000) {
  t <- rep(NA,n)
  t[1] <- 0 #initial value
  for (i in 2:n) {
    y <- rnorm(1,0,1)
    r <- dt(y,df=4)/dt(t[i-1],df=4)*dnorm(t[i-1],0,1)/dnorm(y,0,1) # acceptance ratio
    accprob <- min(1,r)
    if (runif(1) < accprob)
      t[i] <- y
    else
      t[i] <- t[i-1]
  }
  return(t)
}

When I plot the variance of the chain at each iteration it converges to 1.5, even though a t distribution with 4 degrees of freedom has variance 2.

Answer 1

You're using a proposal which hardly ever moves the chain out into the heavy tails of the target $t$-distribution. Once the chain finally finds its way into the tails it typically stays there for a large number of iterations (as it should, see figure below). The chain will eventually converge but only extremely slowly.

If you increase the variance of your proposal or make it dependent on the current value of the chain, say by using y <- rnorm(t[i-1],1), your algorithm should converge faster.