Outputting covariance matrices from CmdStan

Developers
1

I can’t seem to round trip covariance matrices in CmdStanPy, though I think this may be an underlying CmdStan issue that is just easier to see with CmdStanPy. I’m generating some covariance matrices with the wishart_rng and then trying to read them in as data into another model. When I try that, I get errors about some small percentage of these matrices not being positive definite.

It looks like the CSVs that CmdStan outputs are not super precise; is there a way to change that (and see that change in CmdStanPy)? Am I doing something else wrong?

Generate:

data {
  int K;
  int N;
}
generated quantities {
  cov_matrix[K] covs[N];
  for (n in 1:N) 
    covs[n] = wishart_rng(K, diag_matrix(rep_vector(1.0, K)));
}

check:

data {
  int K;
  int N;
  cov_matrix[K] covs[N];
}

This works in Rstan though:

R code 
library("rstan")
options(mc.cores = parallel::detectCores()).
rstan_options(auto_write = TRUE)
draws =stan(model_code="
data {
  int K;
  int N;
}
generated quantities {
  cov_matrix[K] covs[N];
  for (n in 1:N) 
    covs[n] = wishart_rng(K, diag_matrix(rep_vector(1.0, K)));
}",
data=list(N=2, K=10),
algorithm="Fixed_param")

draws = extract(draws)
draws$covs[1,]

checker = stan_model(model_code="
data {
  int K;
  int N;
  cov_matrix[K] covs[N];
}")
bleh = lapply(1:(dim(draws$covs)[1]), 
       function (i) {
         sampling(checker,
                  data=list(N=2, K=10, covs=draws$covs[1,,,]),
                  algorithm="Fixed_param",
                  iter=1, chains=1)
       })
sum(is.na(bleh))

Thanks,
Sean

1 Like
2

The precision in cmdstan CSVs is definitely the problem here as you need a lot of precision here to pass the strict positive definitnes check in Stan. Especially with large matrices.

Rstan reads the doubles directly, as far as I understand.

We would have to make the precision in the CSV files an argument at some point probably.
So not sure you can do anything about that at this point.

2 Likes
3

I thought there was an option to increase the output precision in CmdStan. If not, we could add it.

1 Like
4

I agree. Issue created: https://github.com/stan-dev/cmdstan/issues/922

2 Likes