Double integrate_1d leads to other generated quantities being zero

Interfaces CmdStan
,
1

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!

2

I think what you are experiencing is that an exception occurs in the generated quantities block and that invalidates everything else in the GQ block. So you should fix the “error exceeds tolerance” exceptions.

FYI: in the next release this would produces NaNs instead of zeros.

1 Like
3

Thank you for the quick answer, that is very helpful!

Is there a way to disable this kind of error checking? Because I know that in this case, it doesn’t matter (because for a related problem I did sinh-sinh-integration with a fixed instead of relative error-dependent number of quadrature points and it works just fine) and the errors only arise in the extremes of the integral (e.g., with a tolerance of \sqrt\epsilon: error estimate of integral 4.1223e-294 exceeds the given relative tolerance times norm of integral).

4

Without digging into C++ I dont think there is anything you can do. See generated quantities block returns zeroes for ALL quantities if there's a problem with one · Issue #3114 · stan-dev/stan · GitHub for some more details if you require them. If you do want to tinker with the C++ let me know I can point you in the right direction.

1 Like