JuliaMath / JuliaMath/QuadGK.jl

DomainError in transformed space

Open
#38 1 comment 0 reactions 0 assignees View on GitHub
Dominant language
Julia
Stars
314
Forks
41
PR merge metrics
No merged PRs in 30d

Description

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)
```

Contributor guide

No contributing guide indexed for this repository

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.