Weird quantile failure for a Gamma distribution

[the-noble-argon]

Had a random failure here, normally this sort of things happens when your quantile is greater than 1 but in this case, it's only 0.98989898989899

dG = Gamma{Float64}(0.010316813105574363, 0.43719886895594257)
q    = 0.010101010101010102
quantile(dG, 1-q)

ERROR: DomainError with -5.9123520191840475:
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
Stacktrace:
 [1] throw_complex_domainerror(f::Symbol, x::Float64)
   @ Base.Math .\math.jl:33
 [2] sqrt
   @ .\math.jl:567 [inlined]
 [3] gamma_inc_inv_qsmall(a::Float64, q::Float64)
   @ SpecialFunctions C:\Users\usr\.julia\packages\SpecialFunctions\9pXme\src\gamma_inc.jl:707
 [4] __gamma_inc_inv(a::Float64, minpq::Float64, pcase::Bool)
   @ SpecialFunctions C:\Users\usr\.julia\packages\SpecialFunctions\9pXme\src\gamma_inc.jl:945
 [5] _gamma_inc_inv
   @ C:\Users\usr\.julia\packages\SpecialFunctions\9pXme\src\gamma_inc.jl:933 [inlined]
 [6] gamma_inc_inv
   @ C:\Users\usr.julia\packages\SpecialFunctions\9pXme\src\gamma_inc.jl:915 [inlined]
 [7] gammainvcdf
   @ C:\Users\usr\.julia\packages\StatsFuns\dTYga\src\distrs\gamma.jl:98 [inlined]
 [8] quantile(d::Gamma{Float64}, q::Float64)
   @ Distributions C:\Users\usr\.julia\packages\Distributions\T2SAc\src\univariates.jl:627
 [9] top-level scope
   @ REPL[30]:1

[devmotion]

I checked and the example works with StatsFuns@0.9.15 but errors in the way you describe with StatsFuns@0.9.16. So it seems this bug was introduced by JuliaStats/StatsFuns.jl#113. Looks like a bug in SpecialFunctions at first glance.

[devmotion]

Yes, this is a bug in SpecialFunctions:

julia> using SpecialFunctions

julia> gamma_inc_inv(0.010316813105574363, 1-0.010101010101010102, 0.010101010101010102)
ERROR: DomainError with -5.912352019184924:
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
Stacktrace:
 [1] throw_complex_domainerror(f::Symbol, x::Float64)
   @ Base.Math ./math.jl:33
 [2] sqrt
   @ ./math.jl:567 [inlined]
 [3] gamma_inc_inv_qsmall(a::Float64, q::Float64)
   @ SpecialFunctions ~/.julia/packages/SpecialFunctions/9pXme/src/gamma_inc.jl:707
 [4] __gamma_inc_inv(a::Float64, minpq::Float64, pcase::Bool)
   @ SpecialFunctions ~/.julia/packages/SpecialFunctions/9pXme/src/gamma_inc.jl:945
 [5] _gamma_inc_inv
   @ ~/.julia/packages/SpecialFunctions/9pXme/src/gamma_inc.jl:933 [inlined]
 [6] gamma_inc_inv(a::Float64, p::Float64, q::Float64)
   @ SpecialFunctions ~/.julia/packages/SpecialFunctions/9pXme/src/gamma_inc.jl:915
 [7] top-level scope
   @ REPL[12]:1

(jl_j536nH) pkg> st
      Status `/tmp/jl_j536nH/Project.toml`
  [276daf66] SpecialFunctions v2.1.2

[devmotion]

I opened an issue in the SpecialFunctions repo: JuliaMath/SpecialFunctions.jl#385

[the-noble-argon]

Thanks

[the-noble-argon]

So I saw the pull request being made for "SpecialFunctions.jl" and updated all my packages. I'm seeing fewer failures (the previous example now works) but I still ran into a failure here

dG = Gamma{Float64}(0.0016546512046778552, 0.21745036229505915)
q = 1-0.7070707070707071
quantile(dG,q)

ERROR: DomainError with -9.5893e-320:
log will only return a complex result if called with a complex argument. Try log(Complex(x)).
Stacktrace:
 [1] throw_complex_domainerror(f::Symbol, x::Float64)
   @ Base.Math .\math.jl:33
 [2] _log(x::Float64, base::Val{:ℯ}, func::Symbol)
   @ Base.Math .\special\log.jl:304
 [3] log
   @ .\special\log.jl:269 [inlined]
 [4] __gamma_inc_inv(a::Float64, minpq::Float64, pcase::Bool)
   @ SpecialFunctions C:\Users\usr\.julia\packages\SpecialFunctions\rNt4H\src\gamma_inc.jl:971
 [5] _gamma_inc_inv
   @ C:\Users\usr\.julia\packages\SpecialFunctions\rNt4H\src\gamma_inc.jl:935 [inlined]
 [6] gamma_inc_inv
   @ C:\Users\usr\.julia\packages\SpecialFunctions\rNt4H\src\gamma_inc.jl:917 [inlined]
 [7] gammainvcdf
   @ C:\Users\usr\.julia\packages\StatsFuns\dTYga\src\distrs\gamma.jl:98 [inlined]
 [8] quantile(d::Gamma{Float64}, q::Float64)
   @ Distributions C:\Users\usr\.julia\packages\Distributions\T2SAc\src\univariates.jl:627
 [9] top-level scope
   @ REPL[4]:1

It looks like we're getting closer here, because the log-domain error is for -9.5893e-320 which may as well be zero for all intents and purposes.

[the-noble-argon]

I was also able to trigger a different kind of error with this:

dG = Gamma{Float64}(1.0309015068677239, 0.03724188228734454)
q = 1-0.020202020202020204

quantile(dG, q)
ERROR: DomainError with (1.0309015068677239, -380.7187661907244, 0):
`a` and `x` must be greater than 0 ---- Domain : (0, Inf)
Stacktrace:
 [1] _gamma_inc(a::Float64, x::Float64, ind::Int64)
   @ SpecialFunctions C:\Users\usr\.julia\packages\SpecialFunctions\rNt4H\src\gamma_inc.jl:787
 [2] gamma_inc
   @ C:\Users\usr\.julia\packages\SpecialFunctions\rNt4H\src\gamma_inc.jl:781 [inlined]
 [3] __gamma_inc_inv(a::Float64, minpq::Float64, pcase::Bool)
   @ SpecialFunctions C:\Users\usr\.julia\packages\SpecialFunctions\rNt4H\src\gamma_inc.jl:981
 [4] _gamma_inc_inv
   @ C:\Users\usr\.julia\packages\SpecialFunctions\rNt4H\src\gamma_inc.jl:935 [inlined]
 [5] gamma_inc_inv
   @ C:\Users\usr\.julia\packages\SpecialFunctions\rNt4H\src\gamma_inc.jl:917 [inlined]
 [6] gammainvcdf
   @ C:\Users\usr\.julia\packages\StatsFuns\dTYga\src\distrs\gamma.jl:98 [inlined]
 [7] quantile(d::Gamma{Float64}, q::Float64)
   @ Distributions C:\Users\usr\.julia\packages\Distributions\T2SAc\src\univariates.jl:627
 [8] top-level scope
   @ REPL[10]:1

[Deduction42]

So the last pull request in SpecialFunctions.jl fixed a lot of the earlier issues, but there's still a couple cases where it fails. You mentioned I should pin StatsFuns.jl to version 0.9.15 and not use 0.9.16? I did it on my computer and it worked, but how would I do this inside a Packages.toml file?

I'm saying this because I thought Distributions.jl -> Project.toml specified that its version of StatsFuns had to be "0.9.15" (I figured that the convention would be that Distributions.jl wouldn't be compatible with 0.9.16 because StatsFuns.jl is a pre-release). I clearly still have more to learn about Julia's package manager. How would I fix this in my own package so that only StatsFuns "0.9.15" or earlier is used for its Distributions.jl dependency?

Use something like this?
StatsFuns = "0.9.14 - 0.9.15"

[devmotion]

I clearly still have more to learn about Julia's package manager.

Maybe the Pkg docs are helpful: https://pkgdocs.julialang.org/v1/

How would I fix this in my own package so that only StatsFuns "0.9.15" or earlier is used for its Distributions.jl dependency?

Use something like this? StatsFuns = "0.9.14 - 0.9.15"

You can use

StatsFuns = "= 0.9.15"

if it should be exactly StatsFuns 0.9.15 (https://pkgdocs.julialang.org/v1/compatibility/#Equality-specifier), and e.g.

StatsFuns = "0.9.13 - 0.9.15"

if it should be 0.9.13, 0.9.14, or 0.9.15.

[Deduction42]

Oh, shoot. I missed Equality specifiers when I was skimming the documentation. Thanks. That makes sense.

[the-noble-argon]

I verified this behavior no longer happens with Distributions.jl v0.25.62 using StatsFuns.jl v1.01. Instead I get a reasonable result

dG = Gamma{Float64}(0.010316813105574363, 0.43719886895594257)
q    = 0.010101010101010102
quantile(dG, 1-q)
0.11912515....

so it's fair to call this issue closed. Thanks for all your work!