Priors on Constrained Covariance Matrix

Modeling
1

I am new to Stan but I am in awe of it. Thanks so much for this great piece of software.

When it comes to programming my first proper Stan program, my problem is that I have a covariance matrix, Omega:

parameters {
cov_matrix[3] Omega;
}

and I would like to impose constraints on it plus priors that observe these constraints. Specifically, I would like to normalize Omega[2:3, 2:3] to be an identity matrix and Omega[1,2:3] (and thus Omega[2:3,1]) to be positive. I don’t know how to do while ensuring that Omega remains positive definite and does not mess up the HMC sampling. I would appreciate any help or pointers you could give me.

For context, I am trying to implement a network econometric model in Stan (reference below) for my research, and this is the last bit missing, I think. In coding the model, I observed the advice to start with the most simple model, incrementally adding layers of complexity, making sure I am able to recover parameters from simulated data.

–Christian

Reference
Hsieh, C. S., & Lee, L. F. (2016). A social interactions model with endogenous friendship formation and selectivity. Journal of Applied Econometrics, 31(2), 301-319.

2

Broadly, you can do this by defining a subset of the parameters (I believe in this case it would just be a 2x2 matrix) in the parameter block, and then transforming it into a new variable (i.e. a 3x3 matrix with 2:3,2:3 being the identity) in the transformed parameter block.

3

Thanks, Aaron. I tried creating Omega in the transformed parameters block, but the initial values are rejected because Omega is not positive definite. Here is a minimal example to illustrate my problem:

data{
  // no data
}
transformed data {
 // identity matrix to restrict Omega
 matrix[2, 2] Z_cov;
 Z_cov = diag_matrix(rep_vector(1,2)); 
}
parameters{
  // bounded vector to restrict elements of Omega
  vector<lower=0>[2] Z_e_cov;
}
transformed parameters{
  // restricting Omega
  cov_matrix[3] Omega;
  Omega[2:3,2:3] = Z_cov;
  Omega[2:3,1] = Z_e_cov;
  Omega[1,2:3] = Z_e_cov';
}
model
{
  // sampling won't work: Omega is not positive definite
}
4

You also need to set Omega[1,1]. This will also let you set the range of the other parameters.

You can define Z_e_cov to be on the range [0,1]. Then you can renoromlize by multiplying by sqrt(Omega[1,1]) to ensure the matrix is positive definite. You may need a Jacobian correction depending on where you place your priors.

5

Thanks a lot again. Like this you mean? I added a normal prior on Omega[1,1] and the parser does warn me about possibly needing a Jacobian correction. However, the initial draws for Omega are still not positive definite? Apologies if I am missing something really obvious.

Possibly a minor note, but now the parser also warns me that letting Omega[1,1] appear on the right-hand side of an assignment causes an inefficient deep copy.

data{
  // no data
}
transformed data {
 // identity matrix to restrict Omega
 matrix[2, 2] Z_cov;
 Z_cov = diag_matrix(rep_vector(1,2)); 
}
parameters{
  // bounded vector to restrict elements of Omega
  vector<lower=0, upper=1>[2] Z_e_cov;
  real<lower=0> sigma_sq;
}
transformed parameters{
  // restricting Omega
  cov_matrix[3] Omega;
  Omega[1,1] = sigma_sq;
  Omega[2:3,2:3] = Z_cov;
  Omega[2:3,1] = sqrt(Omega[1,1])*Z_e_cov;
  Omega[1,2:3] = sqrt(Omega[1,1])*Z_e_cov';
}
model
{
  Omega[1,1] ~ normal(1,0.5) T[0,];
}
6

That may just be a numeric stability issue, where the eigenvalue becomes very close to 0. You can probably remove the cov_matrix constraint and just make it a normal matrix because the positive definiteness should be guaranteed by construction. Use sqrt(sigma_sq) rather than sqrt(Omega[1,1]) to avoid the deep copy.

7

All worked out. Thanks so much for your help!

8

No problem. I’m glad it worked out.