# Specifying Correlation in Stan

``````model {
y1 ~ normal(mu1, sigma1);
y2 ~ normal(mu2, sigma2);
}
``````

Hi fellas,

I am relatively new to Stan and have a simple question which I have not been able to figure out how to code in stan for quite sometime now…

I simply wanted to specify a correlation between sigma1 and simga2, the correlation is unknown. Based on my (extensive searching), it seems the best way to do it is by specifying it in the parameters block, using a covariance matrix. However I have no idea how to do this and I haven’t even used the parameters block before.

Would really appreciate it if anyone can help with the code.

1 Like

Do you have any prior information on how `sigma1` and `sigma2` are correlated?

(I think) the most general case is:

``````parameters {
cholesky_factor_corr[2] Lcorr;
vector[2] prior_mu;
vector<lower=0>[2] prior_scale;
vector[2] Z;
}

transformed parameters {
vector[2] sigmas = prior_mu + diag_pre_multiply(prior_scale, Lcorr)*Z;
}

model {
prior_mu ~ multi_normal(0, 1); // you can put prior information about the sigmas here
Z ~ std_normal();
prior_scale ~ std_normal();
Lcorr ~ lkj_corr_cholesky(2.0);

y1 ~ normal(mu1, sigmas[1]);
y2 ~ normal(mu2, sigmas[2]);
}``````
3 Likes

This is also a good read to understand better what is going on.

2 Likes

I would recommend imposing the correlation on the log scale instead. The non-centered approach here doesn’t guarantee that the sigmas remain positive, so you get into rejection sampling and slow mixing. Instead, can parameterize in terms of log(sigma), then just exponentiate to get the positive values out.

2 Likes

Thanks a lot, mathDR! It worked well and the reading was also very useful.

There are actually contemporaneous covariances between the 2. I tried to modify the code to make it work but failed and the link you gave didn’t talk about this. Would you mind taking another look?:)

Here’s my modified code:

``````parameters {
cholesky_factor_corr[2] Lcorr;
matrix[N,2] prior_mu;
vector<lower=0>[2] prior_scale;
matrix[N,2] Z;
}

transformed parameters {
matrix[N,2] sigma = prior_mu + diag_pre_multiply(prior_scale, Lcorr)*Z;
}

model {
Z ~ std_normal();
prior_scale ~ std_normal();
Lcorr ~ lkj_corr_cholesky(2.0);

y1[n] ~ normal(mu1[n], sigma[n,1]);
y2[n] ~ normal(mu2[n], sigma[n,2]);
}
``````

## 16: 17: transformed parameters { 18: matrix[N,2] error = prior_mu + diag_pre_multiply(prior_scale, Lcorr)*Z; ^ 19: }

Error in stanc(file = file, model_code = model_code, model_name = model_name, :
failed to parse Stan model ‘2ef64b1018424bfaa8c4b0cbcd92eefb’ due to the above error.

Thanks for the great suggestion! I will update it in my final model!

Oh! This is a really nice idea! Is this a standard idea for positive parameters?

Hmmm, this seems to think you have `Lcorr` declared in both the parameters block and the transformed parameters block.

Can you post a minimal example (for example, the psuedocode above has `n` but it isn’t declared anywhere, etc.)

Sure, my code is rather long (200+ rows) and the only thing I am trying to modify is to use correlated contemporaneous error. This code should represent the setup…

``````  "data
{
int<lower=0> N;
vector[N] y1;
vector[N] y2;
}

parameters {
matrix[N,2] sigma;

cholesky_factor_corr[2] Lcorr;
matrix[N,2] prior_mu;
matrix<lower=0>[2,N] prior_scale;
matrix[N,2] Z;
}

transformed parameters {
matrix[N,2] sigma = prior_mu + diag_pre_multiply(prior_scale, Lcorr)*Z;
}

model {
Z ~ std_normal();
prior_scale ~ std_normal();
Lcorr ~ lkj_corr_cholesky(2.0);

y1[n] ~ normal(mu1[n], sigma[n,1]);
y2[n] ~ normal(mu2[n], sigma[n,2]);
}
``````

I think my goal is close to this thread(Covariance Matrix Not Symmetric for LKJ prior), however I am too unfamiliar with Stan’s syntax(never used BUGS, etc, before) to figure out the correct way to setup the model.

Okay I don’t see an issue with the snippet you posted. Is that from your code? Also, can you post the stack trace you get (the error messages)?

The issue with the posted snippet is that `sigma` is declared twice. If we delete the declaration of `sigma` from the `parameters` block, then we progress to a different error because the arguments to `diag_pre_multiply` are ill-typed. I can’t tell what the intent is behind modifying @mathDR’s code so that `prior_scale` is a matrix rather than a vector. If we change it back to vector, then next we get an issue with the `Z ~ std_normal()` statement. We can fix with `to_vector(Z) ~ std_normal()`. Then we get some errors with out-of-scope variables `n`, `mu1`, and `mu2`. But if we comment these lines out, then the program compiles at least.