Specifying Correlation in Stan

Modeling
1 
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
2

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
3

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

2 Likes
4

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
5

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]); 
}

Duplicate declaration of variable, name=error; attempt to redeclare as matrix in transformed parameter; previously declared as matrix in parameter
error in ‘model14ae1914545_2ef64b1018424bfaa8c4b0cbcd92eefb’ at line 18, column 72

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.

6

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

7

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

8

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.)

9

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.

10

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)?

11

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.