# Bayesian inference on the correlation parameter of a bivariate normal

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?

• First, your lif will be easier if you take (WLOG) $\mu_1 = \mu_2 = 0$, $\sigma_1 = \sigma_2 = 1$ and bother about the shift and the scaling once you are done. Did you try taking an inverse Wishart prior and then conditioning on $\sigma_1, \sigma_2$? I imagine the integration might not be too horrible, especially if you take the simplification I suggested. I can't do the entire thing now but I think this is the way to go, given that inverse wishart is the conjugate prior for the covariance. I can try that later but tell me first if this is something you tried. Oct 2, 2017 at 3:26
• Thanks for the suggestion. I did not try that but I did think about that. It does not make sense to me that the prior for $\rho$ depends on $\sigma_{11}$ and $\sigma_{22}$. My intuition tells me that the posterior under that prior would the that the posterior under the inverse Wishart prior assuming that $\sigma_{11}$ and $\sigma_{22}$ unknown and then conditioning on $\sigma_{11}$ and $\sigma_{22}$. Oct 2, 2017 at 3:31
• Also, this seems relevant: stats.stackexchange.com/questions/58420/… Oct 2, 2017 at 3:31
• I think we are suggesting the same thing. Alternatively, you can take the inverse wishart, condition that for you $\sigma_2, \sigma_1$ and use the result as a conjugate prior. Oct 2, 2017 at 3:35
• But that requires that the inverse wishart distribution conditioning on $\sigma_{11}$ and $\sigma_{22}$ having a tractable form, which I am not sure of. Will try to work on the math to see if that is true. Oct 2, 2017 at 3:43

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])} • Thanks for the answer and the bound is really clever. Looking at $f(\rho)$, seems that we can use a inverse gamma prior truncated between $(0,1)$, this motivates me to transform the variable $\alpha = 1-\rho^2$. But then I realize that under the transformation you will need to decide the sign of $\rho$. So maybe we can transform $\rho=h\times \sqrt{1-\alpha}$ where $\alpha = 1-\rho^2$ and $h=\pm 1$. I feel we might be able to directly sample $(\alpha, h)$ with the help of the rejection sampling strategy you proposed. I will work on the math details and post a follow-up. Thanks again! Oct 2, 2017 at 16:54
• I tried this change of variable (on the positive real line as sign does not matter) but the Jacobian gets in the way... Oct 3, 2017 at 4:56
• +1 Note that the "$\beta$" in $(1)$ is not the same as the $\beta$ in the enveloping distribution.
– whuber
Oct 3, 2017 at 14:03

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

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. • For many values of$\beta$and$\alpha$(especially small$\alpha\$, which one might use for relatively uninformative priors) it looks like a Laplace approximation could be extremely poor, so it would be prudent to characterize the situations where it could work and to quantify how much error it might introduce.
– whuber
Oct 3, 2017 at 22:44