Negative Binomial CDF and Quantile do not always roundtrip

[mborland]

Failure can be found here: https://github.com/scipy/scipy/pull/17432/checks?check_run_id=10961417302
Found on PR: scipy/scipy#17432

Minimal reproducer:

import numpy as np
import numpy.testing as npt
from scipy import stats

dist = stats.nbinom(5, 0.5)
supp = 6
cdf_supp = dist.cdf(supp)  # ideally 0.7255859375
cdf_supp0 = cdf_supp - 10*np.spacing(cdf_supp)  # 0.7255859374999989
npt.assert_array_equal(dist.ppf(cdf_supp0), supp)  # both should be 6

Note in their code says:

In very rare cases, the finite precision calculation of ppf(cdf(supp))
can produce an array in which an element is off by one. We nudge the
CDF values down by 10 ULPs help to avoid this.

At 15 ULPs it seems to be correct.

[mdhaber]

Note that this failure only occured on our macos_arm64_test. This is the same platform in which we see the overflow warnings.

Actually, they could be related. I didn't realize it when I wrote

Maybe more accurately, I don't recall a case in which there is an overflow warning and the code produces an incorrect solution : )

but we are silencing overflow warnings from the nbinom PPF calculation, too.

[jzmaddock]

I wonder if this is truly a bug, but more particularly this comes down to what we actually want the quantile to do here as different libraries take slightly different approaches.

Our goal for quantile is to find integer k such that negative_binomial(k, n, p) >= probability, and we approach that by bracketing the root, and then refining the brackets until both ends generate a probability >= the target. Now this is conceptually the right thing to do, and for random probability values it homes in on the result very rapidly. However, any attempt to round trip like this will automatically generate an edge case where the two ends of the bracketed root are either side of but very close to an integer - so the algorithm keeps refining the brackets until one end randomly pops over the integer boundary, and both ends then agree which way to round (by rounding UP). An alternative is to treat it as a continuous function, find the root and round to nearest, but that can generate absurdities too. Or we could introduce some sort of fudge factor whereby if it looks like it's dead on an integer we pick that one rather than risking an arbitrary round up?

So I guess my question would be: is the original bug a result of a needed feature, or simply someone doing a naive sanity check?

[mdhaber]

Our goal for quantile is to find integer k such that negative_binomial(k, n, p) >= probability

In general, the PPF is the smallest integer k such that the CDF >= probability, right?

is the original bug a result of a needed feature, or simply someone doing a naive sanity check?

I can't speak for whether users will actually need the edge behavior to be perfect for science/engineering use. Nonetheless, users check cases like this to verify their understanding of the function and to check that it's working as they expect, and they will find and report edge cases like this as a bug. This can cause confusion and erode confidence in SciPy, so we try to get the edge cases right, too.

an edge case where the two ends of the bracketed root are either side of but very close to an integer - so the algorithm keeps refining the brackets

It sounds like you are using a typical bracketing rootfinder, say bisection, and the ends of the bracket are a and b. When solving for the PPF, do you evaluate the CDF at a and b or floor(a) and floor(b)? (It would not matter in SciPy because our cdf(x) is equivalent to cdf(floor(x)) for discrete distributions.)

Either way, can't you terminate as soon as the bracket width drops below 1? At that point (or as soon as floor(b) == floor(a +1), really) you can round up from the left or down from the right. From what you described, it sounded like the bracket width was getting very tiny.

Please forgive me if I'm confusing myself. I've certainly done that before when thinking about this problem!

[jzmaddock]

Me too!

I think there may be a solution though:

1. Terminate the root finder if the bracket on the root drops below 1 then:
2. Return the ceil of either end if the two ends do not straddle an integer, or
3. Double check the ceil of the left end otherwise to see if it satisfies the condition, before falling back to the ceil of the right end.

We might even save one or two CDF evaluations per quantile call that way.

[mdhaber]

I think the situation is improved by terminating early, yup. Less likely to run into the trouble you described.

This is making me wonder about how

one end randomly pops over the integer boundary

In the round-trip case. If the CDF were numerically monotonic, this wouldn't happen, right? This would mean that the bracket was no longer valid because the root is at the boundary, but now both ends of the bracket are on one side of the boundary. Not a bracket.

But since there are more floating point values within the support of a discrete distribution than there are between 0 and 1 (the range of the CDF), it seems that the CDF will always be a bit step-wise. So when "one end randomly pops over the integer boundary", does this happen such that cdf(a) == cdf(b)==desired_prob?

I think if this is the case and you round up, I think you run the risk of returning one integer too high.

This would be very rare if you terminate as soon as possible, but I can see how it would happen if you continue to refine the bracket.

This brings another question to mind - how is the CDF defined if cdf(x) is not equivalent to cdf(floor(x))? Is the result at non-integral x have a statistical interesting or is it like an interpolant?

[jzmaddock]

Just looking through our code at this.... looks like this case should already be handled (albeit it may be possible to be more efficient), so I pushed this test case to investigate more what's actually happening rather than just speculating: https://github.com/boostorg/math/pull/938/files#diff-f08ad2489004236ba31ecf0b1d938d4433525c9bfecd328bba33398f65642e18

But it's passing the tests, what am I missing? https://drone.cpp.al/boostorg/math/1163/1/2

[mdhaber]

But it's passing the tests, what am I missing?

Maybe the architecture?
#935 (comment)

[jzmaddock]

Ah, arm64, missed that, thought it was M1, I think we can emulate arm64/Linux at least, will try and see.

[mborland]

Ah, arm64, missed that, thought it was M1, I think we can emulate arm64/Linux at least, will try and see.

macOS ARM should be limited to M1

[mborland]

@jzmaddock I don't know if it matters but these are the policies being used in the build:

cpp_args = [
   '-DBOOST_MATH_DOMAIN_ERROR_POLICY=ignore_error',
   '-DBOOST_MATH_EVALUATION_ERROR_POLICY=user_error',
   '-DBOOST_MATH_OVERFLOW_ERROR_POLICY=user_error'
   '-DBOOST_MATH_OVERFLOW_ERROR_POLICY=user_error',
   '-DBOOST_MATH_PROMOTE_DOUBLE_POLICY=false'
 ]

On M1 long doubles are just doubles.