Covariance Matrix not symmetric

Modeling
,
1

Hi,

I am receiving an Informational message that the scale parameter is not symmetric and that the covariance matrix is not symmetric. It looks like this the warning ( the code works fine) as below.
I am using symmetric matrices as you will see below in the initial values setting.

Any idea?

Warning Message

Exception: inv_wishart_lpdf: scale parameter is not symmetric. scale parameter[2,3] = -1.33734e-09, but scale parameter[3,2] = -1.13908e-08 (in
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: multi_normal_lpdf: Covariance matrix is not symmetric. Covariance matrix[1,2] = 1.27815e-08, but Covariance matrix[2,1] = 2.29241e-08

Stan Code

data {
int<lower=0> N;
vector[16] x; 
matrix[N,16] y;
}
parameters {
vector[N] D; 
vector[N] S0 ; 
vector[N] f; 
vector[N] Dstar; 
real <lower=0> sigma;
vector[3] mu;
matrix[3,3] Sigma;
vector[3] mu_hyper;
matrix[3,3] Sigma_hyper;
real<lower=2> df_Sigma_hyper;
matrix[3,3] Psi_Sigma_hyper;
}
model {
for (n in 1:N){
S0[n] ~ uniform(50, 250);
matrix[n,3] theta;
theta[n,1] = D[n];
theta[n,2] = Dstar[n];
theta[n,3] = f[n];
theta[n] ~ multi_normal(mu, Sigma);
}
sigma ~ uniform(0,10); 
mu ~ multi_normal(mu_hyper, Sigma_hyper);
Sigma ~ inv_wishart(df_Sigma_hyper, Psi_Sigma_hyper);
for (n in 1:N){
for (k in 1:16){
y[n,k] ~ normal(S0[n]*((1-f[n])*exp(-D[n]*x[k])+f[n]*exp(-Dstar[n]*x[k])),sigma); 
}
}
}

My initial values are as below :

mu, zeros(3, 1)
Sigma, (diag([1e-3, 1e-3, 1e-3])+diag([1e-3, 1e-3, 1e-3])‘) / 2
mu_hyper, zeros(3, 1)
Sigma_hyper, (diag([1e3, 1e3, 1e3])+diag([1e3, 1e3, 1e3])’) / 2
Psi_Sigma_hyper’,(diag([1e3, 1e3, 1e3])+diag([1e3, 1e3, 1e3])‘) / 2
df_Sigma_hyper’, 3

2

I believe the covariance matrix parameters need to be declared as cov_matrix instead of matrix. For example,

matrix[3,3] Sigma;

becomes

cov_matrix[3] Sigma;
4

I wonder how you are sending the initial values in to Stan. In the past, I think I have seen these sorts of messages when there are problems with the initial values.

6

I am not quite sure there, partly because I have never used Stan with Matlab. If you can post a dataset that leads to the warning, I might be able to see whether the problem also happens in R. But if the message only appears at the start of sampling, I suspect it is not a major problem.

Also, I have not spent much time on the model, but I suspect you can avoid the loop in the model block for increased efficiency. Maybe declare theta in the parameter block as

vector[3] theta[N];

and then, in the model block, remove the loop and do

theta ~ multi_normal(mu, Sigma);

If you want to keep them, you could define D, Dstar, and f in a transformed parameters block via entries of theta.

7

I had a bit more of a look at it via rstan (with data sent privately) and got a different error message (double free or corruption (out). I was able to fit a simplified model without hyperpriors and with theta declared as a parameter, but I don’t have time to figure out what else is going wrong. I also wonder whether f should be bounded between 0 and 1, because it seems like you are taking a weighted average of two components. I also think more vectorization is possible to avoid the double loop but ignored it here.

The model I got to work (some lines are commented out using //):

data {
  int<lower=0> N;
  vector[16] x; 
  matrix[N,16] y;
}
parameters {
  vector[3] theta[N];
  real S0[N];
  real <lower=0> sigma;
  //vector[3] mu;
  //cov_matrix[3] Sigma;
  //vector[3] mu_hyper;
  //cov_matrix[3] Sigma_hyper;
  //real<lower=2> df_Sigma_hyper;
  //cov_matrix[3] Psi_Sigma_hyper;
} transformed parameters {
  real D[N] = theta[,1];
  real Dstar[N] = theta[,2];
  real f[N] = theta[,3];
} model {
  //theta ~ multi_normal(mu, Sigma);
  //sigma ~ uniform(0,300); 
  //mu ~ normal(0, 1); //multi_normal(mu_hyper, Sigma_hyper);
  //Sigma ~ inv_wishart(3, diag_matrix(rep_vector(1,3)));
  for (n in 1:N){
    for (k in 1:16){
      y[n,k] ~ normal(S0[n] * ((1 - f[n])*exp(-D[n]*x[k]) + f[n]*exp(-Dstar[n]*x[k])), sigma);
    }
  }
}