# Priors on Constrained Covariance Matrix

**URL:** <https://discourse.mc-stan.org/t/priors-on-constrained-covariance-matrix/3858>\
**Category:** Modeling\
**Created:** [April 13, 2018, 12:34am UTC](https://discourse.mc-stan.org/t/priors-on-constrained-covariance-matrix/3858 "2018-04-13T00:34:37Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Christian1](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/christian1/32/2167_2.png) [@Christian1](https://discourse.mc-stan.org/u/Christian1)\
**Post date:** [April 13, 2018, 12:34am UTC](https://discourse.mc-stan.org/t/priors-on-constrained-covariance-matrix/3858/1 "2018-04-13T00:34:37Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![aaronjg](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/aaronjg/32/3768_2.png) [@aaronjg](https://discourse.mc-stan.org/u/aaronjg)\
**Post date:** [April 13, 2018, 12:41am UTC](https://discourse.mc-stan.org/t/priors-on-constrained-covariance-matrix/3858/2 "2018-04-13T00:41:37Z")

</div>

> [@Christian1](#):
>
> 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])

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.

---

<div class="post-metadata">

**Author:** ![Christian1](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/christian1/32/2167_2.png) [@Christian1](https://discourse.mc-stan.org/u/Christian1)\
**Post date:** [April 13, 2018, 1:17am UTC](https://discourse.mc-stan.org/t/priors-on-constrained-covariance-matrix/3858/3 "2018-04-13T01:17:47Z")

</div>

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
}
```

---

<div class="post-metadata">

**Author:** ![aaronjg](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/aaronjg/32/3768_2.png) [@aaronjg](https://discourse.mc-stan.org/u/aaronjg)\
**Post date:** [April 13, 2018, 1:33am UTC](https://discourse.mc-stan.org/t/priors-on-constrained-covariance-matrix/3858/4 "2018-04-13T01:33:34Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Christian1](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/christian1/32/2167_2.png) [@Christian1](https://discourse.mc-stan.org/u/Christian1)\
**Post date:** [April 13, 2018, 2:16am UTC](https://discourse.mc-stan.org/t/priors-on-constrained-covariance-matrix/3858/5 "2018-04-13T02:16:31Z")

</div>

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,];
}
```

---

<div class="post-metadata">

**Author:** ![aaronjg](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/aaronjg/32/3768_2.png) [@aaronjg](https://discourse.mc-stan.org/u/aaronjg)\
**Post date:** [April 13, 2018, 2:59am UTC](https://discourse.mc-stan.org/t/priors-on-constrained-covariance-matrix/3858/6 "2018-04-13T02:59:08Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Christian1](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/christian1/32/2167_2.png) [@Christian1](https://discourse.mc-stan.org/u/Christian1)\
**Post date:** [April 13, 2018, 3:33am UTC](https://discourse.mc-stan.org/t/priors-on-constrained-covariance-matrix/3858/7 "2018-04-13T03:33:31Z")

</div>

All worked out. Thanks so much for your help!

---

<div class="post-metadata">

**Author:** ![aaronjg](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/aaronjg/32/3768_2.png) [@aaronjg](https://discourse.mc-stan.org/u/aaronjg)\
**Post date:** [April 13, 2018, 6:02pm UTC](https://discourse.mc-stan.org/t/priors-on-constrained-covariance-matrix/3858/8 "2018-04-13T18:02:39Z")

</div>

No problem. I’m glad it worked out.