Metropolis Hastings for Posterior of Bivariate Normal

Question

As an exercise, I am trying to implement metropolis hastings to draw samples from the posterior distribution of a bivariate normal: $$ (X,Y) \sim N \left( (0,0)\begin{bmatrix}1 & \rho \\ \rho & 1 \end{bmatrix}\right) $$ to estimate the parameter $\rho$. The setup is as follows:

1. We have a Jeffreys prior for $\rho$, that is, the distribution of $\rho$ is proportional to $1/(1-\rho^2)^{1/2}$.
2. We compute the posterior distribution for $\rho$, and arrive at the fact that: $$ f(\rho | \{(x_i,y_i)\}_{i=1}^{n}) =\propto \frac{1}{2\pi^2}\frac{1}{{1-\rho^2}} \prod_{i=1}^n\exp\left(\frac{-1}{2(1-\rho^2)}[x_i^2 -2\rho x_iy_i+y_i^2]\right) $$

We draw samples from a uniform random walk kernel. That is, given an estimate $\rho_m$, we draw an estimate:

$$ \rho^* \sim \mathrm{Unif} (\rho_n-0.1, \rho_n + 0.1) $$

The acceptance function is thus given by: $$ \alpha = \min \left(1, \frac{f(\rho^*|\{(x_i,y_i)\}_{i=1}^{n})}{f(\rho_m|\{(x_i,y_i)\}_{i=1}^{n})}\right) $$ Where $(x_i,y_i)_{i=1}^n$ are samples that have been drawn before running the chain. We start with $\rho_0 = 0.1$.

I have implemented this using the following R code:

gensamples <- function (rho, N){ #Draw correlated normals
      X1 = rnorm(N)
      X2 = rnorm(N)
      X3 = rho*X1 + sqrt(1-rho^2)*X2
      Y1 = X1
      Y2 = X3
      samples = matrix(c(Y1,Y2),nrow = N, ncol=2)
      return (samples)
}

l_ratio <- function(samples,rho,rho_) #Likelihood ratio
  return (
    exp(
    sum(
    -1/(2*(1-rho**2))*(samples[,1]**2-2*rho*samples[,1]*samples[,2]+samples[,2]**2)  
  + 
    1/(2*(1-rho_**2))*(samples[,1]**2-2*rho_*samples[,1]*samples[,2]+samples[,2]**2)
  )
  )
  )


prior_ratio <- function(rho,rho_)
  return (
    (1/(1-rho**2)**(1/2))
    /
    (1/(1-rho_**2)**(1/2))
          )

posterior_ratio<- function(samples,rho,rho_){ #Use Bayes Formula
  return(l_ratio(samples,rho,rho_)*prior_ratio(rho,rho_))

}

samples = gensamples(rho = 0.2,1000)
burn_in = 10000
iterations = burn_in + 1000
rho_0 = 0.1
rho = rho_0
s = c(0)
for (i in 1:iterations){
  rho_ = runif(1, min = rho -0.1, max = rho+0.1)
  alpha = min(1, 1/posterior_ratio(samples,rho,rho_))
  if (runif(1)<alpha){
    rho = rho_
  }
  if (i >burn_in)
    s = c(s,rho)
}
n = seq_along(s)
m = cumsum(s)/n
m2 = cumsum(s*s)/n
v = (m2 -m*m)*(n/(n-1))
plot(m,type = 'l')
plot(v,type = 'l')

However, it is giving me issues. A quick look at plots tells me the chain converges, but it seems to be very biased. If I use $0.2$, like in the sample above, the usual estimate comes out to about $0.1$-$0.15$. Could anyone let me know if I'm doing something wrong in the calculation?

Answer 1

The posterior should be $$f(\rho | \{(x_i,y_i)\}_{i=1}^{n}) \propto \frac{1}{{(1-\rho^2})^{n+1/2}} \exp\left(\frac{-1}{2(1-\rho^2)}\sum_{i=1}^n[x_i^2 -2\rho x_iy_i+y_i^2]\right)$$ and the part $(1-\rho^2)^{n}$ is missing from the likelihood ratio in the R code.

The proposal being $\mathrm{Unif} (\rho_n-0.1, \rho_n + 0.1)$, there is a positive probability that the simulated value stands outside $(-1,1)$ for values of $\rho_{n}$ close enough to $\pm 1$. The target density in the R code should therefore be set to zero outside $(-1,1)$ to accommodate such entries.

There is thus an issue with the likelihood function as coded since, if I use instead

library(mvtnorm)
l_ratio <- function(samples,rho,rho_)
  return(
   exp(
    sum( dmvnorm(samples,sigma=matrix(c(1,rho,rho,1),2),log=TRUE)) -
      sum( dmvnorm(samples,sigma=matrix(c(1,rho_,rho_,1),2),log=TRUE))   
      )
  )

I recover an MCMC sample converging to the generating value of $\rho$. In the R code provided in the question, it should be

l_ratio <- function(samples,rho,rho_) #Likelihood ratio
  return (
          sum(.5*log(abs(1-rho_**2))+
        1/(2*(1-rho**2))*(samples[,1]**2-2*rho*samples[,1]*
                            samples[,2]+samples[,2]**2)  
        - .5*log(abs(1-rho**2))-
          1/(2*(1-rho**2))*(samples[,1]**2-2*rho*samples[,1]*
                            samples[,2]+samples[,2]**2)
      )
  )