# Quadrature from x=0 to x=Inf?
**URL:** https://discourse.julialang.org/t/quadrature-from-x-0-to-x-inf/44751
**Category:** General Usage
**Created:** [August 11, 2020, 4:23pm UTC](https://discourse.julialang.org/t/quadrature-from-x-0-to-x-inf/44751 "2020-08-11T16:23:24Z")
**Posts on this page:** 12
**Page:** 1
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [August 11, 2020, 4:23pm UTC](https://discourse.julialang.org/t/quadrature-from-x-0-to-x-inf/44751/1 "2020-08-11T16:23:24Z")
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?
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [August 11, 2020, 4:26pm UTC](https://discourse.julialang.org/t/quadrature-from-x-0-to-x-inf/44751/2 "2020-08-11T16:26:01Z")
QuadGK already supports this:
```julia
julia> using QuadGK
julia> quadgk(x -> exp(-x), 0, Inf, rtol=1e-5)
(0.9999999997018256, 2.779309769038174e-6)
```
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [August 11, 2020, 8:18pm UTC](https://discourse.julialang.org/t/quadrature-from-x-0-to-x-inf/44751/3 "2020-08-11T20:18:30Z")
> [@stevengj](#):
>
> `quadgk(x -> exp(-x), 0, Inf, rtol=1e-5)`
Thanks for tip. I have problems making it work. Consider:
```julia
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:
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
julia> quadgk(x->poisson(x),0,Inf)
(0.8338114480884391, 5.717008672714917e-9)
```
However, for \lambda = 2, `quadgk` crashes:
```julia
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?
---
### Author:  [@giordano](https://discourse.julialang.org/u/giordano)
#### Post date: [August 11, 2020, 8:25pm UTC](https://discourse.julialang.org/t/quadrature-from-x-0-to-x-inf/44751/4 "2020-08-11T20:25:20Z")
Can you modify you poisson function to show the input when the output is nan? Is it that 0.9375?
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [August 11, 2020, 8:28pm UTC](https://discourse.julialang.org/t/quadrature-from-x-0-to-x-inf/44751/5 "2020-08-11T20:28:23Z")
Not exactly sure how to do that. Some kind of `try - catch`?
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [August 11, 2020, 8:29pm UTC](https://discourse.julialang.org/t/quadrature-from-x-0-to-x-inf/44751/6 "2020-08-11T20:29:39Z")
I should say that Wolfram Alpha computes the result without problem. [But is inconvenient to work with.]
---
### Author:  [@giordano](https://discourse.julialang.org/u/giordano)
#### Post date: [August 11, 2020, 8:32pm UTC](https://discourse.julialang.org/t/quadrature-from-x-0-to-x-inf/44751/7 "2020-08-11T20:32:08Z")
I’m not at the computer, but I’d just do something like
```julia
if isnan(....)
println(...)
end
```
before returning
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [August 11, 2020, 8:38pm UTC](https://discourse.julialang.org/t/quadrature-from-x-0-to-x-inf/44751/8 "2020-08-11T20:38:00Z")
OK – I’m trying this:
```julia
# "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
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...
```
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [August 11, 2020, 8:47pm UTC](https://discourse.julialang.org/t/quadrature-from-x-0-to-x-inf/44751/9 "2020-08-11T20:47:56Z")
Hm… could it be that my function is numerically poor? I returned to the original `poisson` function:
```julia
julia> poisson(1871.5213495195621,2)
NaN
julia> poisson(BigFloat(1871.5213495195621),2)
1.946532642091951927364768066088547511008985129095238634528284797438114570307282e-4751
```
Would there be a way around this?
---
### Author:  [@giordano](https://discourse.julialang.org/u/giordano)
#### Post date: [August 11, 2020, 8:56pm UTC](https://discourse.julialang.org/t/quadrature-from-x-0-to-x-inf/44751/10 "2020-08-11T20:56:01Z")
> [@BLI](#):
>
> could it be that my function is numerically poor?
It is:
```julia
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`
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [August 11, 2020, 8:56pm UTC](https://discourse.julialang.org/t/quadrature-from-x-0-to-x-inf/44751/11 "2020-08-11T20:56:27Z")
OK – I solved the problem by rephrasing the function:
```julia
# "Continuous" Poisson distribution
function poisson_ln(x,λ=1)
fX_ln = x*log(λ) - λ - log(gamma(x+1))
return exp(fX_ln)
end
```
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [August 11, 2020, 9:38pm UTC](https://discourse.julialang.org/t/quadrature-from-x-0-to-x-inf/44751/12 "2020-08-11T21:38:03Z")
> [@BLI](#):
>
> `quadgk` crashes:
(That’s an exception, not a “crash”.) The basic issue is that your `poisson` function returns `NaN` for large arguments:
```julia
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
julia> function poisson(x,λ=1)
return exp(-λ + x * log(λ) - loggamma(x+1))
end
```
which then gives:
```julia
julia> poisson(1e4,2)
0.0
```
at which point `quadgk` has no problem:
```julia
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.)
> [@BLI](#):
>
> OK – I solved the problem by rephrasing the function:
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.)