Here’s the test:

```
TEST(ProbBetaBinomial, lcdf_matches_mathematica) {
int n = 8;
int N = 10;
double alpha = 3.0;
double beta = 1.0;
EXPECT_FLOAT_EQ(-0.5500463, (stan::math::beta_binomial_lcdf(n, N, alpha, beta)));
}
```

The message is:

```
[ RUN ] ProbBetaBinomial.lcdf_matches_mathematica
unknown file: Failure
C++ exception with description "called from function 'F32', hypergeometric
function 3F2 does not meet convergence conditions with given arguments.
a1: 1, a2: 12, a3: -1, b1: 10, b2: -1, z: 1" thrown in the test body.
[ FAILED ] ProbBetaBinomial.lcdf_matches_mathematica (0 ms)
[----------] 3 tests from ProbBetaBinomial (0 ms total)
```

With the arguments given to F32 (a1: 1, a2: 12, a3: -1, b1: 10, b2: -1, z: 1) this makes sense—b2 is a negative integer and the denominator is not allowed to be a negative integer so it’s not supposed to converge. Also a1+a2+a3 > b1+b2 so it’s not supposed to converge. The current behavior is to throw and that’s fine, but the arguments I’m feeding here seem like perfectly normal arguments to the beta-binomial.

I think we need a continuation to have this working or an alternative call to the F32. This is a test I added so maybe I don’t understand the beta-binomial arguments? Any suggestions?

The beta-binomial CDF should be well-defined whenever `alpha > 0`

, `beta > 0`

, `N >= 0`

, and `0 <= n <= N`

. I’m not sure why we need to get mathematica involved. We can test a CDF with a finite range using the log PMF (`beta_binomial_lpmf`

), with `log_sum_exp`

to sum on the log scale.

The derivative calculation blows up when the shape parameters happen to be integers.

Actually even the LCDF calculation itself is undefined (currently) for one set of parameters, but it is only for a *single* set of parameters (well, whenever Beta == 1.0 and some other conditions). Currently that causes a rejection at that point in parameter space which I think is in line with what happens in other related situations (e.g.-when we have a laplace density the gradient is not defined at the peak). I’m looking at the 3F2 identities to see if there’s a way to modify the calculation for that point without resorting to summing the pmf.

Mathematica is involved for a sanity check, unless you happen to have a better source for 3F2 values and gradients.

Mathematica’s fine for generating true answers. We don’t need it for beta-binomial cdfs, which we can define with the beta-binomial pmf, but I have no idea what 3F2 is.

We should be able to handle integer-valued parameters for the beta binomial cdfs—they should be well defined and they’re definitely likely to be used by users. Is there a way to fall back on what we used to be doing if the inputs are integers or has this always been broken? If it’s always been broken, then no need to fix it along with everything else, but if it used to work, it’d be a shame to break it.

It used to be broken. I’m evaluating a solution that avoids using the PMF

so if that works I’ll put it in, otherwise I’ll write up a separate issue

and skip it. 3F2 is a hypergeometric function, look at the definition of

the beta-binomial CDF (say in Wikipedia).

Sounds good and thanks for the pointer.