JuliaStats / StatsFuns.jl

Mathematical functions related to statistics.
Other
233 stars 41 forks source link

gammalogcdf(k, θ, x) returns error for valid arguments #150

Closed the-noble-argon closed 1 year ago

the-noble-argon commented 1 year ago

I'm getting some errors in my code for some edge cases of the gammalogcdf function:

k = 42648.50647826826
θ = 2.2498007956420723e-5
x = 0.6991377135675367
gammalogcdf(k, θ, x)

which returns

ERROR: DomainError with -1.1090323347913147:
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:301
 [3] log
   @ .\special\log.jl:267 [inlined]
 [4] _gammalogcdf(k::Float64, θ::Float64, x::Float64)
   @ StatsFuns C:\Users\user\.julia\packages\StatsFuns\0IeNQ\src\distrs\gamma.jl:60
 [5] gammalogcdf(k::Float64, θ::Float64, x::Float64)
   @ StatsFuns C:\Users\user\.julia\packages\StatsFuns\0IeNQ\src\distrs\gamma.jl:47
 [6] top-level scope
   @ c:\Users\user\Documents\julia_debug.jl:7

Internally this is because of the log(drummond1F1(1.0, 1 + k, xdθ)) expression on line 60 of gamma.jl. For this input we get

k = 42648.50647826826
xdθ = 31075.538550870202
drummond1F1(1.0, 1 + k, xdθ)
-1.1090323347913147

which is a negative number that you can't take the log of.

the-noble-argon commented 1 year ago

I know that for this case, k is absolutely huge. Maybe you can case out the situation where k is really big and so the gamma distribution is essentially Normal?