Quadrature from x=0 to x=Inf?

General Usage
1

I have previously used QuadGK over a finite interval. Now I need to do numeric integration from x=0 to x=Inf, and wonder if there is a routine that supports this? [I need to check scaling, etc. of the “continuous” Poisson distribution.]

Any suggestions?

2

QuadGK already supports this:

julia> using QuadGK

julia> quadgk(x -> exp(-x), 0, Inf, rtol=1e-5)
(0.9999999997018256, 2.779309769038174e-6)
5 Likes
3

Thanks for tip. I have problems making it work. Consider:

using SpecialFunctions, QuadGK, Plots
# "Continuous" Poisson distribution
function poisson(x,λ=1)
    return λ^x*exp(-λ)/gamma(x+1)
end

I can plot the “continuous” Poisson distribution, see grey lines below:

image
image1704×1110 214 KB

The markers have been generated using the Distribution.jl package.
The “continuous” Poisson function is not a proper pdf because the integral from zero to infinity differs from unity. So I’d like to find scaling constants depending on \lambda.

With \lambda = 1, it works:

julia> quadgk(x->poisson(x),0,Inf)
(0.8338114480884391, 5.717008672714917e-9)

However, for \lambda = 2, quadgk crashes:

julia> quadgk(x->poisson(x,2),0,Inf)
DomainError with 0.9375:
integrand produced NaN in the interval (0.875, 1.0)

Stacktrace:
 [1] evalrule(::QuadGK.var"#14#23"{var"#35#36",Float64}, ::Float64, ::Float64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::typeof(LinearAlgebra.norm))...

Can you replicate this crash, or is it related to a previously reported LinearAlgebra problem on my i9-9900 home computer?

4

Can you modify you poisson function to show the input when the output is nan? Is it that 0.9375?

5

Not exactly sure how to do that. Some kind of try - catch?

6

I should say that Wolfram Alpha computes the result without problem. [But is inconvenient to work with.]

7

I’m not at the computer, but I’d just do something like

if isnan(....)
    println(...)
end

before returning

8

OK – I’m trying this:

# "Continuous" Poisson distribution
function poisson(x,λ=1)
    if isnan(λ^x*exp(-λ)/gamma(x+1))
        println(x)
    end
    return λ^x*exp(-λ)/gamma(x+1)
end

leading to

julia> quadgk(x->poisson(x,2),0,Inf)
1871.5213495195621
DomainError with 0.9375:
integrand produced NaN in the interval (0.875, 1.0)

Stacktrace:
 [1] evalrule(::QuadGK.var"#14#23"{var"#5#6",Float64}, ::Float64, ::Float64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::typeof(LinearAlgebra.norm)) at C:\Users...
9

Hm… could it be that my function is numerically poor? I returned to the original poisson function:

julia> poisson(1871.5213495195621,2)
NaN

julia> poisson(BigFloat(1871.5213495195621),2)
1.946532642091951927364768066088547511008985129095238634528284797438114570307282e-4751

Would there be a way around this?

10

It is:

julia> function poisson(x,λ=1)
           if isnan(λ^x*exp(-λ)/gamma(x+1))
               @show λ^x; @show exp(-λ); @show gamma(x+1)
           end
           return λ^x*exp(-λ)/gamma(x+1)
       end
poisson (generic function with 2 methods)

julia> poisson(1871.5213495195621,2)
λ ^ x = Inf
exp(-λ) = 0.1353352832366127
gamma(x + 1) = Inf
NaN

You have Inf/Inf

2 Likes
11

OK – I solved the problem by rephrasing the function:

# "Continuous" Poisson distribution
function poisson_ln(x,λ=1)
    fX_ln = x*log(λ) - λ - log(gamma(x+1))
    return exp(fX_ln)
end
12

(That’s an exception, not a “crash”.) The basic issue is that your poisson function returns NaN for large arguments:

julia> poisson(1e4,2)
NaN

This is due to a spurious overflow that arises in how you defined it: λ^x == 2^(1e4) == Inf and gamma(1e4+1) == Inf, so you are computing Inf / Inf which gives NaN. In fact, the final result is representable if you rearrange your computation of the poisson function to combine all of the exponents into a single exp call

julia> function poisson(x,λ=1)
           return exp(-λ + x * log(λ) - loggamma(x+1))
       end

which then gives:

julia> poisson(1e4,2)
0.0

at which point quadgk has no problem:

julia> quadgk(x->poisson(x,2),0,Inf)
(0.9470194210852405, 3.1332612554072925e-9)

(On a separate note, realize that integrating over an infinite interval with a generic routine like quadgk is internally transformed into a singular integrand on a finite interval. Because the integrand is effectively singular, it can be more challenging to integrate accurately, so sometimes you have to increase the requested tolerance rtol from the default ~ 1e-8.)

(Alternatively, you could use Gauss–Legendre quadrature for this sort of integrand, but that requires more manual intervention in terms of properly rescaling the integrand and selecting the number of quadrature points.)

Ah, sounds like you mostly figured it out for yourself. However, be sure to use loggamma(x+1) and not log(gamma(x+1)) — the latter can spuriously overflow. e.g. for x=1e4, the log(gamma(x+1)) call gives Inf whereas loggamma(x+1) gives 82108.92783681436. (In this particular case, the overflow may be harmless because the final exp call underflows to 0.0 in both cases.)

4 Likes