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