Hello,

I wanted to generate a correlated bivariate normal distribution based on two data points (one for each of the 2 variables), and one value of a standard deviation that describes the uncertainty in those values. These are the 4 parameters in the data block.

My data are therefore just used to describe marginal normal distributions, providing no information on the correlation, but I have specified prior models for the parameters in the covariance matrix, specifically, the standard deviations, and the correlation coefficient (based on a separate analysis).

I expected that the posterior bivariate distribution would show correlated samples, because of the inclusion of a correlation coefficient. The posterior model of the correlation coefficient looks good (unchanged from my prior description), but the posterior model of the bivariate normal distribution looks like a plot of 2 independent variables.

Below is the model I am describing. If you can give me any advice on what I am overlooking, that would be much appreciated, but please bear in mind that I am not a statistician.

```
data {
real <lower = 0> x1Meas; // Observed value of first parameter
real <lower = 0> x1Error; // Std dev of uncertainty of measurement of first parameter
real <lower = 0> x2Meas; // Observed value of second parameter
real <lower = 0> x2Error; // Std dev of uncertainty of measurement of secondparameter
}
parameters {
vector <lower = 0> [2] mu; // Means of multi-normal distribution
real <lower = 0> sigma[2]; // Std devs of multi-normal distribution
real <lower = -1, upper = 1> rho; // Correlation coefficient
vector <lower = 0> [2] y; // Parameeter assigned to multi-normal distribution
}
transformed parameters {
cov_matrix[2] cov = [[ sigma[1] ^ 2 , sigma[1] * sigma[2] * rho],
[sigma[1] * sigma[2] * rho, sigma[2] ^ 2 ]];
}
model {
// Likelihood
y ~ multi_normal(mu, cov);
// Priors
rho ~ normal (0.7, 0.2) T[-1, 1];
sigma ~ exponential(1);
mu[1] ~ normal(x1Meas, x1Error) T[0,];
mu[2] ~ normal(x2Meas, x2Error) T[0,];
x1Meas ~ normal (4, 1) T[0,];
x1Error ~ exponential(1);
x2Meas ~ normal (20, 5) T[0,];
x2Error ~ exponential(0.5);
}
```

In summary, I expected correlated samples (from parameter y), but they appear to be independent.

Many thanks.