Bayesian inference on the correlation parameter of a bivariate normal

Question

Suppose that the data $y_i$'s are from the following bivariate normal

$$y_i\sim \mathcal{N}\bigg(\mu,\left[ {\begin{array}{cc} \sigma_{11} & \sqrt{\sigma_{11}\sigma_{22}}\rho \\ \sqrt{\sigma_{11}\sigma_{22}}\rho & \sigma_{22} \\ \end{array} } \right]\bigg).$$

Suppose that $\mu$, $\sigma_{11}$ and $\sigma_{22}$ are all known and one wants to learn the posterior distribution of $\rho$ under some prior distribution, for instance,

$$\dfrac{\rho+1}{2}\sim beta(2,2).$$

My question is, can the posterior be directly sampled from? Is there any conjugate prior that can result in some tractable posterior?

I worked through the tedious math and have the following

$$L(y_1,\ldots,y_n|\rho)\propto(1-\rho^2)^{-\frac{n}{2}}\exp\bigg\{-\dfrac{\sum_{i=1}^{n}\tilde{y}_{i1}^2 - 2\rho\tilde{y}_{i1}\tilde{y}_{i2}+\tilde{y}_{i2}^2}{2(1-\rho^2)}\bigg \},$$ where $\tilde{y}_{i1} = (y_{i1}-\mu_1)/\sqrt{\sigma_{11}}$ and $\tilde{y}_{i2} = (y_{i2}-\mu_2)/\sqrt{\sigma_{22}}$. However this does not remind me of any possible conjugate prior.

Or, if there is no conjugate prior available, could any one suggest a good rejection sampler strategy? What could an efficient rejection proposal distribution?

Any suggestions?

Answer 1

Since $$L(y_1,\ldots,y_n|\rho)\propto(1-\rho^2)^{-\frac{n}{2}}\exp\bigg\{-\dfrac{\sum_{i=1}^{n}\tilde{y}_{i1}^2 - 2\rho\tilde{y}_{i1}\tilde{y}_{i2}+\tilde{y}_{i2}^2}{2(1-\rho^2)}\bigg \}$$ is a function of $\rho$ of the form $$(1-\rho^2)^{-\alpha}\exp\bigg\{-\dfrac{\beta}{1-\rho^2}-\dfrac{\gamma\rho}{1-\rho^2}\bigg \}\qquad (1)$$this leads to an exponential family choice of a conjugate prior (with $\alpha,\beta>0$ and $|\gamma|<\beta$). While this is not a standard distribution, as far as I know, an accept-reject solution may be available, using the bound$$\dfrac{\sum_{i=1}^{n}\tilde{y}_{i1}^2 - 2\rho\tilde{y}_{i1}\tilde{y}_{i2}+\tilde{y}_{i2}^2}{2(1-\rho^2)}\ge\dfrac{\sum_{i=1}^{n}\min[(\tilde{y}_{i1}-\tilde{y}_{i2})^2,(\tilde{y}_{i1}+\tilde{y}_{i2})^2]}{2(1-\rho^2)}$$ may help, even though I have not been able to find a simple way to simulate from $$f(\rho)\propto(1-\rho^2)^{-\alpha}\exp\bigg\{-\dfrac{\beta}{1-\rho^2}\bigg\}$$ Obviously, since the above density (1) is bounded, a brute-force accept-reject method based on a Uniform works, if slowly [in the picture below the acceptance rate is 2.35%!]:

targ=function(x,a,b,c){
  (1-x*x)^{-a}*exp(-b/(1-x*x)-c*x/(1-x*x))}

upb=function(a,b,c){
  return(optimise(targ,maximum=TRUE,a=a,b=b,c=c,inte=c(-1,1))$obj)}

simz=function(n,a=1,b=1,c=0){
  bon=upb(a,b,c)
  rejcz=integrate(targ,low=-1,upp=1,a=a,b=b,c=c)$val/2/bon
  uniz=runif(ceiling(2*n/rejcz),min=-1,max=1)
  vuniz=runif(ceiling(2*n/rejcz))
  samplz=uniz[vuniz<targ(uniz,a,b,c)/bon]
  return(samplz[1:n])}

Answer 2

It seems that a Laplace approximation works quite well. Below I define the log-likelihood and its gradient. Note that I change the variable so that the support is on the real line for better Laplace approximation performance. I use a logit transformation, i.e., $\rho = \dfrac{2}{e^{-x}+1}-1$.

likfcn <- function(x, a, b, c, log = FALSE) {
  rho = 2 / (exp(-x) + 1) - 1
  aux = 1 - rho^2
  llik = -a * log(aux) - b / aux - c * rho / aux
  if (log) {
    return(llik)
  } else {
    return(exp(llik))
  }
}

llikGradient <- function(x, a, b, c) {
  rho = 2 / (exp(-x) + 1) - 1
  aux = 1 - rho^2
  llik_gradient_rho = (2 * a * rho * aux - 2 * b * rho - c - c * rho^2) / aux^2
  return(llik_gradient_rho * 2 * exp(x) / (1 + exp(x))^2)
}

Then I simulate some data from some simple bivariate gaussian and calculate the corresponding a, b and c.

# simulation data
set.seed(2017)
suppressMessages(require(mvtnorm))
Sigma = matrix(c(1, .5, .5, 1), 2, 2)
n = 20
y = rmvnorm(n, c(0, 0), Sigma)
a = n / 2
b = sum(y[, 1]^2 + y[, 2]^2) / 2
c = -sum(y[, 1] * y[, 2])

I then maximize the log-likelihood and obtain the hessian at the maximum. With the information I can create a Laplace approximation.

fn <- function(x) {
  -likfcn(x, a, b, c, log = TRUE)
}
gr <- function(x) {
  -llikGradient(x, a, b, c)
}

optim_res = optim(par = 0, fn = fn, gr = gr, method = "BFGS", lower = -Inf, upper = Inf,hessian = TRUE)

# laplace approximation
laplace_mean = optim_res$par
laplace_var = 1 / optim_res$hessian

Then I compare the Laplace approximation to the true posterior.

# compare the laplace approximation to the true posterior
x_rho = seq(-.99, .99, .01)
targ=function(x,a,b,c){
  (1-x*x)^{-a}*exp(-b/(1-x*x)-c*x/(1-x*x))}
normalizing_const = integrate(targ,a,b,c,low=-1,upp=1)$val
y_dens_true = targ(x_rho,a,b,c) / normalizing_const
plot(x_rho, y_dens_true, 'l')
x = log((x_rho+1)/2/(1-(x_rho+1)/2))
y_dens_laplace = dnorm(x, laplace_mean, sqrt(laplace_var)) * (1 / (1+x_rho) + 1 / (1-x_rho))
lines(x_rho, y_dens_laplace, col = "red")

where the black line is the density plot of the truth and the red is the Laplace approximation. As can be seen, Laplace approximation performs well, especially noting that the data size $n=20$ is not big.