Hi all,

I’m encountering a bug (?) when trying to integrate out two parameters in a double integral via nested `integrate_1d`

calls. Specifically, doing the integration causes other generated quantities to be zero. I am not sure if it is even a bug and I would like to come up with a MWE, but maybe before doing that and heading over to github, somebody has an idea.

The problem is quite straightforward. This is my generated quantities block:

```
generated quantities {
vector<lower=0>[3] var_props[C, K];
{
matrix[N_R, C] garbage_matrix;
for (c in 1:C) {
////////////////////// Calculate variance proportions
for (k in 1:K) {
var_props[c, k] = [
square(loading_T[k] * sd_ability_T),
square(loading_R[k] * sd_ability_R[c]),
square(loading_I[k] * sd_ability_I)
]';
var_props[c, k] /= sum(var_props[c, k]);
}
}
for (c in 1:C) {
real distr_item_pars[JKp2Kp3];
////////////////////// Log likelihood (marginal)
distr_item_pars[1:3] = { mean_ability[c], sd_ability_R[c], sd_ability_I };
distr_item_pars[4:Kp3] = to_array_1d(loading_R);
distr_item_pars[Kp4:K2p3] = to_array_1d(loading_I);
for (k in 1:K) {
distr_item_pars[(K2p3p1 + (k - 1) * J) : (K2p3 + k * J)] = to_array_1d(threshold[k]);
}
for (r in 1:N_R) {
// garbage_matrix[r, c] = integrate_1d(
// cc_gpcm_likelihood_lvl2, // Function invoking integrate_1d again
// negative_infinity(),
// positive_infinity(),
// append_array(
// distr_item_pars,
// to_array_1d(loading_T * ability_T[targets[rater_starts[r] : rater_ends[r]]]')
// ),
// { tol },
// integration_data[r],
// tol
// );
}
}
}
}
```

When I run it like that with the `garbage matrix`

, i.e., the double `integrate_1d`

part commented out, I get ordinary draws for the `var_props`

:

```
> double_integral_without_garbage_cmdstanfit$draws("var_props")
# A draws_array: 100 iterations, 4 chains, and 27 variables
, , variable = var_props[1,1,1]
chain
iteration 1 2 3 4
1 0.28 0.29 0.23 0.32
2 0.35 0.30 0.25 0.36
3 0.28 0.23 0.25 0.31
4 0.29 0.32 0.27 0.27
5 0.35 0.24 0.27 0.29
...
```

However, when I de-comment the aforementioned part, then, strangely, the draws for `var_props`

are zero even though their calculation precedes the double integration and both calculations are independent.

```
> double_integral_with_garbage_cmdstanfit$draws("var_props")
# A draws_array: 100 iterations, 4 chains, and 27 variables
, , variable = var_props[1,1,1]
chain
iteration 1 2 3 4
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 0 0 0 0
...
```

The draws for the other (primary or transformed, not generated) parameters, however, are identical and sensible in both cases (provided I’m setting a seed), so this issue seems to concern solely the generated quantities block. I thought it might be a memory issue, but `cmdstan`

consumes almost no memory while running. I’m using `cmdstanr`

because with `rstan`

there is another (or maybe it is the same) issue in that `rstan`

somehow fails to create a fit object. Moreover, if I replace the `integrate_1d`

call with just `garbage_matrix[r, c] = 1`

, the issue is not caused, so it should be specifically the `integrate_1d`

call.

In the output I’m getting a lot of “error exceeds tolerance” exceptions, but I’m sure this is due to the model as well as the double integration part, but in any case, it should have no bearing on other generated quantities.

Maybe I’m overlooking something obvious or this is a known issue but I couldn’t find something that sounded similar. In the meantime, I’ll try to make a MWE and see if it is a unique problem of my model (integrand?) or if this issue arises in toy cases, too.

Thank you!