3
$\begingroup$

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?

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

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)
      )
  )
$\endgroup$
1
  • $\begingroup$ thank you! This worked $\endgroup$ – rubikscube09 Jun 3 '20 at 3:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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