Hierarchical multi_normal

Modeling
1

I would like to use a model containing input-dependent noise, simply put where y ~ N(0,C) and C = D + diag(v) where v ~ N(0, E). I have coded this in Stan (in a manner that is similar to my model) below:

data {
  int<lower=1> N;
  vector[N] y;
  matrix[N,N] Covar1;
}
parameters {
  real<lower=0> var1;
  real<lower=0> var2;
  vector[N] var3;
  real<lower=0> var4;
}
model {
  matrix[N, N] Covar2;
  matrix[N, N] Covar3;
  
  Covar2 = var1 * Covar1 + diag_matrix(var3);
  Covar3 = var2 * Covar1 + diag_matrix(rep_vector(var4,N));

  var1 ~ normal(0,1000);
  var2 ~ normal(0,1000);
  var4 ~ normal(0,1000);

  var3 ~ multi_normal(rep_vector(1000,N), Covar3);
  y ~ multi_normal(rep_vector(0, N), Covar2);
}

I get an error saying that the Covar2 matrix is not positive definite and “Initialization between (-2, 2) failed after 100 attempts.”. This would be the case if var3 was negative, but it’s not as its mean is set to 1000 and the Covar3 matrix does not contain values to make var3 negative (I have tested this through simulation).

If I replace the 4th model line as such:

Covar2 = var1 * Covar1 + diag_matrix(var3+rep_vector(100,N));

the model runs fine.

How do I go about modelling without resorting to this addition “trick”? What am I missing? Is there an initialization issue here? I am new to Stan.

2

The matrix needs to be positive definite. One way to ensure positive definiteness is add a large enough positive constant to the diagonal.

Unless var1 can take on values in the thousands, this prior’s unlikely to provide any useful information and you might as well just use an improper uniform if that’s what you really want.

Using diag_matrix is almost always a bad idea because it adds a bunch of zeroes. Much better to add in a loop (there’s a new add-diagonal function coming out in Stan 2.18 I think to encapsulate this). Loops aren’t slow in Stan like they are in R or Python.

Stan’s much more efficient when parameterized with Cholesky factors, but it looks like you’re scaling constants and trying to encourage them to be positive definite.

You certainly don’t want to be adding negative values to the diagonal, as that leads you away from positive definiteness. so maybe var3 just needs a lower bound of zero if you know Covar1 is already positive definite.

3

I didn’t realize you could lower bound vectors:

vector<lower=0>[N] var3

That fixed my issue. Thanks for the optimization tips too :)