search

JuliaStats / Distributions.jl

A Julia package for probability distributions and associated functions.
Other
1.08k stars 410 forks source link

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

Open dfok opened 3 weeks ago

dfok commented 3 weeks ago

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 https://github.com/JuliaStats/Distributions.jl/blob/b356da03a189d023cdb8467c61806a8a11dcb262/src/quantilealgs.jl#L8 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 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).