gamma_inc_inv throws an DomainError for valid arguments

[devmotion]

First observed in JuliaStats/Distributions.jl#1510:

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

The same error occurs also in SpecialFunctions 2.0.0, so it wasn't introduced by #377.

[devmotion]

The problem is that the expression inside the sqrt in

SpecialFunctions.jl/src/gamma_inc.jl

Line 707 in 8b4aba0

eta = sqrt(-2/a*log(q*gammax(a)*sqrt(twoπ/a)))

is negative. Since a >= 0 and 0 <= q <= 1, AFAICT it is non-negative iff q < exp(a * (log(a) - 1)) / gamma(a). It seems in general this is not implied by the condition in

SpecialFunctions.jl/src/gamma_inc.jl

Line 944 in 8b4aba0

elseif !pcase && minpq < min(0.02, exp(-1.5*a)/gamma(a)) && a < 10 #small q

.

[andreasnoack]

I've just tried the original Fortran version accompanying the paper on which our implementation is based and it also fails for this input. From a first read of the paper, it wasn't obvious where eta comes from but I'll try to look more carefully.

[devmotion]

I had a quick look at the Fortran code and the paper as well yesterday. I assumed it would error since I couldn't spot an obvious difference from our implementation but didn't check it. It wasn't immediately clear to me either from skimming through the paper why it fails and what's the explanation for the if condition. There was a note in the zip file with the Fortran code about a more recent paper with an improved (?) code - maybe the problem is fixed there?

[Deduction42]

Just to help track down when things started to go off the rails, the input that triggered all this was from Distributions.jl. I recently updated and some previously smooth code started failing. Before updating yesterday, things ran quite smoothly. In fact, I was able to verify this by running the same code on a different machine that I haven't updated yet.

Distributions v0.25.35 (an older version that works):

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

However, on the same machine, I have a version of SpecialFunctions which fails for the inputs

SpecialFunctions v1.8.1:

gamma_inc_inv(0.010316813105574363, 1-0.010101010101010102, 0.010101010101010102)
DomainError with -5.912352019184924:
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).

[devmotion]

Thanks, this will be fixed soon. There's already a PR to SpecialFunctions with a fix.

It's not an issue with Distributions but with SpecialFunctions (therefore I opened the issue here after I was notified about it in an issue in Distributions). You don't see it in the computer where you haven't updated your packages because it means you haven't updated StatsFuns: In StatsFuns 0.9.16 many (inv)(log)cdf functions that used Rmath up until StatsFuns 0.9.15 were replaced with native Julia implementations. And this uncovered the bug here in SpecialFunctions. Currently the example will work with any version of Distributions if you use StatsFuns 0.9.15 but fail if you use 0.9.16.

[Deduction42]

Ah okay, well when that's up and running, the solution I'm building will put it to a very thorough test. I make thousands of calls to Normal and Gamma quantiles to assess scale for quadrature. I could probably do some guesswork to capture feasible ranges but quantiles make for a very clean and stable solution.