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 (

prior_ratio <- function(rho,rho_)
  return (

posterior_ratio<- function(samples,rho,rho_){ #Use Bayes Formula


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?

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

1 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

l_ratio <- function(samples,rho,rho_)
    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 (
        - .5*log(abs(1-rho**2))-
  • $\begingroup$ thank you! This worked $\endgroup$ Jun 3, 2020 at 3:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.