Quantile bisect widen interval

[arnaud-ma]

Fixes #1869

The issue

Specifically for mixture models, in this code

Distributions.jl/src/mixtures/mixturemodel.jl

Lines 446 to 450 in da420cd

	function quantile(d::UnivariateMixture{Continuous}, p::Real)
	ps = probs(d)
	min_q, max_q = extrema(quantile(component(d, i), p) for (i, pi) in enumerate(ps) if pi > 0)
	quantile_bisect(d, p, min_q, max_q)
	end

Let's call $F(x) = \text{cdf}(d, x)$ and $F_i(x) = \text{cdf}(\text{component}(d, i), x)$, same $F^{-1}$ and $F_i^{-1}$ for quantile. quantile_bisect requires that

$$ F(q_\text{min}) \leq p \leq F(q_\text{max}) $$

otherwise an error is thrown. Mathematically, in the code above, the condition is guaranteed. But if we take the definition of $q_\text{min}$ and $q_\text{max}$, the condition becomes

$$ F(\min_i F_i^{-1}(p)) \leq p \leq F(\max_i F_i^{-1}(p)) . $$

And numerically, for some specific edge cases dealing with very small numbers, this condition is not satisfied. For example,

julia> d = MixtureModel([Normal(0, 1), Normal(eps(), 1)], [0.5, 0.5])
MixtureModel{Normal{Float64}}(K = 2)
components[1] (prior = 0.5000): Normal{Float64}(μ=0.0, σ=1.0)
components[2] (prior = 0.5000): Normal{Float64}(μ=2.220446049250313e-16, σ=1.0)


julia> d1, d2 = components(d)
2-element Vector{Normal{Float64}}:
 Normal{Float64}(μ=0.0, σ=1.0)
 Normal{Float64}(μ=2.220446049250313e-16, σ=1.0)

julia> p = 0.001
0.001

julia> cdf(d, quantile(d1, p)), cdf(d, quantile(d2, p))
(0.000999999999999983, 0.0009999999999999853)

julia> quantile(d, p)
ERROR: ArgumentError: [-3.090232306167818, -3.0902323061678176] is not a valid bracketing interval for `quantile(d, 0.001)`

I found another example different than "almost identical" components, but still with the same issue:

julia> d = MixtureModel(
           [Exponential(2.1656391815606704e-5), Exponential(79.77140291062331)],
           [8.83763312432469e-25, 1.0],
       )
MixtureModel{Exponential{Float64}}(K = 2)
components[1] (prior = 0.0000): Exponential{Float64}(θ=2.1656391815606704e-5)
components[2] (prior = 1.0000): Exponential{Float64}(θ=79.77140291062331)


julia> d1, d2 = components(d)
2-element Vector{Exponential{Float64}}:
 Exponential{Float64}(θ=2.1656391815606704e-5)
 Exponential{Float64}(θ=79.77140291062331)

julia> p = 0.1066425045830068
0.1066425045830068

julia> cdf(d1, quantile(d1, p)) <= p, cdf(d2, quantile(d2, p)) <= p
(true, true)

julia> quantile(d, p)
ERROR: ArgumentError: [2.442157681386963e-6, 8.995697253353159] is not a valid bracketing interval for `quantile(d, 0.1066425045830068)`

I did not found any issues such that $p < F(q_\text{min})$ though. It's always $p > F(q_\text{max})$. Maybe it is impossible with mixtures?

The fix

Because it seems specific to mixture models, we could simply make the interval wider in quantile(d::UnivariateMixture{Continuous}, p::Real) by adding some small number to max_q and subtracting some small number from min_q. But I think it would be better to make the quantile_bisect function more robust, by allowing it to widen the interval if the condition is not satisfied. It does not impact the old case at all, and it would fix the issue for mixture models, and potentially for other / future distributions with similar issues.

[codecov-commenter]

Codecov Report

❌ Patch coverage is 66.66667% with 10 lines in your changes missing coverage. Please review.
✅ Project coverage is 86.54%. Comparing base (da420cd) to head (dc1f6ef).
⚠️ Report is 13 commits behind head on master.

Files with missing lines	Patch %	Lines
src/quantilealgs.jl	66.66%	10 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2055      +/-   ##
==========================================
- Coverage   86.62%   86.54%   -0.08%     
==========================================
  Files         148      148              
  Lines        8903     8928      +25     
==========================================
+ Hits         7712     7727      +15     
- Misses       1191     1201      +10     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[arnaud-ma]

For the missing coverage, as I said above, I tried but did not find any distribution that triggers $p < F(q_\text{min})$

[arnaud-ma]

I guess now that it is using Roots.jl the bug is "upstream". One day it will be solved!