Bug in quantile() for some special cases of a MixtureModel

[dfok]

The problem

For rather exotic instances of a MixtureModel the quantile() method fails. The underlying reason is that the initial bisection interval (in quantile_bisect()) does not contain the solution. However, both ends of the interval are very close to the solution. The current implementation of quantile_bisect() checks for equality of the start and end of the interval. Changing this to approximately equal would fix this bug.

Example

using Distributions
d = MixtureModel([Normal(0, 1), Normal(eps(), 1)], [0.999, 0.001])
quantile(d, 0.001)

gives as output

ERROR: ArgumentError: [-3.090232306167818, -3.0902323061678176] is not a valid bracketing interval for `quantile(d, 0.001)`
Stacktrace:
 [1] quantile_bisect(d::MixtureModel{Univariate, Continuous, Normal{Float64}, Categorical{Float64, Vector{Float64}}}, p::Float64, lx::Float64, rx::Float64)
   @ Distributions ~/.julia/packages/Distributions/ji8PW/src/quantilealgs.jl:21
 [2] quantile(d::MixtureModel{Univariate, Continuous, Normal{Float64}, Categorical{Float64, Vector{Float64}}}, p::Float64)
   @ Distributions ~/.julia/packages/Distributions/ji8PW/src/mixtures/mixturemodel.jl:447
 [3] top-level scope
   @ REPL[6]:1

I'm currently using

julia> versioninfo()
Julia Version 1.10.3
Commit 0b4590a5507 (2024-04-30 10:59 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × 12th Gen Intel(R) Core(TM) i7-1260P
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, alderlake)
Threads: 1 default, 0 interactive, 1 GC (on 16 virtual cores)

and Distributions v0.25.109.

The solution(?)

Changing line 8 of quantilealgs.jl

Distributions.jl/src/quantilealgs.jl

Line 8 in b356da0

if rx == lx

to if rx ≈ lx seems to do the trick. I'm not sure whether this has bad side effects. However, after this change all tests in Distributions.jl still succeed.

An appropriate new test case to catch this bug would be

d = MixtureModel([Normal(0, 1), Normal(eps(), 1)], [0.999, 0.001])
@test cdf(d, quantile(d, 0.001) ) ≈ 0.001

I can submit this as a pull request if desired (would need to figure out how to do that).