Incorrect quantiles for beta-binomial distribution #2090

It looks like quantile and cdf calculations for beta-binomial distributions can sometimes be off from their true value (here off by >20%):

julia> using Distributions

julia> d = BetaBinomial(1000, 0.1, 20)
BetaBinomial{Float64}(n=1000, α=0.1, β=20.0)

julia> maximum(d)
1000

julia> quantile(d, 1.0)
777

julia> cdf(d, 777)          
0.9999999999999984

This bug occurs specifically at quantile(d, 1.0).

Possible cause (AI-diagnosed)

In univariate/discrete/betabinomial.jl we have the following:

cdf(d::BetaBinomial, k::Int) = integerunitrange_cdf(d, k)     # accurate recurrence

_categorical(d::BetaBinomial) = Categorical(map(Base.Fix1(pdf, d), support(d)))
quantile(d::BetaBinomial, p::Float64) = quantile(_categorical(d), p) - 1

quantile re-derives the whole probability vector from pointwise pdf calls (each an independent logbeta-based formula evaluated in isolation) and hands it to Categorical, whose quantile does a plain left-to-right running sum against p. For p = 1.0 that running sum only needs to accumulate any tiny excess (here +1.3e-15 total, distributed unevenly across 1001 terms) to cross the >= 1.0 threshold hundreds of indices early, since most of the mass in this skewed case sits at very small k. cdf, defined independently via integerunitrange_cdf, does not share this error and correctly stays under 1.0 until k = 1000. The two functions on the same object were never checked against each other. [copied out of a longer response, may lack context]