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.