QuadGK.jl
QuadGK.jl copied to clipboard
DomainError in transformed space
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)
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.