Given known bivariate normal means and variances, update correlation estimate, $P(\rho)$, with new data?

Question

I'm dealing with two correlated random variables which are modeled via a bivariate normal distribution. I have values for the means ($\mu_x, \mu_y$) and individual variances ($\sigma_x, \sigma_y$) of the variables (effectively the marginalize normal distributions). Given new observations of the random variables, could I update my estimate via a Bayesian inference for the correlation coefficient, $P(\rho)$?

Thoughts so far:

Starting with the sampling distribution:

$P(x,y|\rho) = \frac{1}{2 \pi \sigma_x \sigma_y \sqrt{1-\rho^2}} \exp\left( -\frac{1}{2(1-\rho^2)}\left[ \frac{(x-\mu_x)^2}{\sigma_x^2} + \frac{(y-\mu_y)^2}{\sigma_y^2} - \frac{2\rho(x-\mu_x)(y-\mu_y)}{\sigma_x \sigma_y} \right] \right)$

We want to update the prior, $P(\rho)$ (not sure of the prior form) with some observation $(x_i,y_i)$:

$P(\rho|x_i,y_i)=P(x_i,y_i)^{-1}*P(x_i,y_i|\rho)*P(\rho)$

What is throwing me off is the form of the prior distribution. I know for my process that the $\rho$ is bounded due to positive correlation. How can I factor this information into the definition of the prior? Additionally, could I implement this procedure in a recursive manner?

Other thoughts:

Perhaps I am overthinking this problem entirely. I could just compute the Pearson correlation coefficient directly given all measured data. However, my expected sample size will be small, and I also know bounds for the correlation ($0<\rho<1$). Given the small sample and a priori known constraints, I still am leaning Bayesian.

Answer 1

A quick Google search led me to this article: Estimating the Correlation in Bivariate Normal Data With Known Variances and Small Sample Sizes. In this article, they discussed several possible priors: "uniform", "Jeffreys", and "arc-sine". Specifically, take the "uniform" prior for example, it assumes that $\rho$ follows a uniform distribution on $[-1,1]$.

This leads to a density of the form: $$ f(\rho | X, Y) \propto \prod\limits_i{P(x_i, y_i | \rho)} $$

If a full Bayesian estimator is desired, then equation $\hat{\rho}^{(6)}$ on page 35 of that article needs to be computed. This can be done using a Monte Carlo integration.

Or, if one only needs the MAP estimate of $\rho$, it can be done using Gibbs sampling. I followed one citation of the first article: (Barnard 2000) Modeling Covariance Matrices in terms of Standard Deviations and Correlations, with Application to Shrinkage. In this citation, it is discussed in detail how a uniform prior for $\rho$ is derived out of the Inverse-Wishart distribution for the covariance matrix $\Sigma$ and how Gibbs sampling was carried out. (Equation (8)).

I also tried searching for other references, but was unable to find possible conjugate priors for $\rho$ and approaches to get the Bayesian estimation of $\rho$ without doing numerical integration or sampling. It is an interesting problem to look at by the way.