DomainError in transformed space

[tbeason]

This is possibly related to #30, but not sure. I am integrating a function that is basically the Gamma PDF (very spiked near zero, but tapers off as x -> infinity). The function does not return NaN for any value (I fix gpdf(0)=0, or else it would give NaN here).

quadgk(x->gpdf(x,l1,d1,s1),0,Inf)
ERROR: DomainError with 0.9999999999999964:
integrand produced NaN in the interval (0.9999999999999929, 1.0)
Stacktrace:
 [1] evalrule(::getfield(QuadGK, Symbol("##12#18")){getfield(Main, Symbol("##41#42")),Float64}, ::Float64, ::Float64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::typeof(LinearAlgebra.norm)) at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\evalrule.jl:41
 [2] adapt(::Function, ::Array{QuadGK.Segment{Float64,Float64,Float64},1}, ::Float64, ::Float64, ::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Int64, ::Float64, ::Float64, ::Int64, ::typeof(LinearAlgebra.norm)) at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:39
 [3] do_quadgk(::getfield(QuadGK, Symbol("##12#18")){getfield(Main, Symbol("##41#42")),Float64}, ::Tuple{Float64,Float64}, ::Int64, ::Nothing, ::Nothing, ::Int64, ::typeof(LinearAlgebra.norm)) at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:28
 [4] handle_infinities(::getfield(Main, Symbol("##41#42")), ::Tuple{Float64,Float64}, ::Int64, ::Nothing, ::Nothing, ::Int64, ::Function) at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:83
 [5] #quadgk#21(::Nothing, ::Nothing, ::Int64, ::Int64, ::Function, ::typeof(quadgk), ::Function, ::Float64, ::Float64) at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:155
 [6] #quadgk#20 at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:155 [inlined]
 [7] quadgk(::Function, ::Int64, ::Float64) at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:151
 [8] top-level scope at REPL[109]:1

The error indicates that the NaN arises in the handle_infinities function, which I guess it where you transform the integration domain from (0,Inf) to something usable. If I use a large number instead of Inf, I do not get the error.

MWE

using SpecialFunctions, QuadGK
l1 = exp(-1.2604534253)
d1 = -0.916695708
s1 = 8.428294026

function gpdf(t,lambda,delta,sigma)
    t == 0 && return 0.0
    d2 = delta^(-2)
    constantpiece = lambda*abs(delta/sigma)/gamma(d2)
    LT = lambda*t
    varpiece = LT^(1/(delta*sigma)-1) * exp(-LT^(delta/sigma))
    val = constantpiece*varpiece
    return val
end

quadgk(x->gpdf(x,l1,d1,s1),0,Inf)
quadgk(x->gpdf(x,l1,d1,s1),0,1e30)

[stevengj]

For integrating from t = 0 to Inf, it uses the transformation f(t)dt goes to f(u/(1-u)) / (1-u)^2 du for u=0..1

My guess is that you are asking for more accuracy than it is possible for QuadGK to attain for your integrand, due to roundoff error limitations, so it keeps refining and refining near u=1 until finally 1-u == 0 due to roundoff errors, and you get f(Inf) / 0 == 0 / 0 which gives NaN.

(Can you print out the t values where your gpdf function is being evaluated? I'm guessing it is evaluating a point where t=Inf.)

Maybe try lowering the rtol you are requesting.