# Metropolis Hastings for Posterior of Bivariate Normal

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?

• By uncorrelated do you mean just drawing from $[0,1]$, without taking into account the position of the previous sample $\rho_n$? – rubikscube09 May 31 '20 at 18:35
• No problem, I understand it's truly awful code - one of my first attempts using R – rubikscube09 May 31 '20 at 18:46

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

• thank you! This worked – rubikscube09 Jun 3 '20 at 3:24