ekmett / ekmett/integration

Failure of an integral from 0 to infinity

Open
#1 8 comments 0 reactions 0 assignees View on GitHub
Dominant language
Haskell
Stars
12
Forks
2
PR merge metrics
No merged PRs in 30d

Description

Hello,
I've found a failure of `nonNegative` for a certain function.

```
module Bug
where
import Data.Number.Erf (normcdf)
import Numeric.Integration.TanhSinh

pnorm :: Double -> Double
pnorm = normcdf

integrand :: Double -> Double -> Double -> Double -> Double
integrand nu delta t x = pnorm(t*x /sqrt nu - delta) * x**(nu-1) * exp(-x*x/2)

integral_0_to_b :: Double -> Double -> Double -> Double -> [Result]
integral_0_to_b nu delta t = trap (integrand nu delta t) 0

integral_0_to_Inf :: Double -> Double -> Double -> [Result]
integral_0_to_Inf nu t delta = nonNegative trap (integrand nu delta t)

-- works:
test1 :: [Result]
test1 = integral_0_to_Inf 23 0 0
-- fails:
test2 :: [Result]
test2 = integral_0_to_Inf 24 0 0
-- works:
test3 :: [Result]
test3 = integral_0_to_b 24 0 0 30
```

Results:
```
*Bug> test1
[Result {result = 1.8259133852602253e10, errorEstimate = 1.825319374957239e10, evaluations = 13},Result {result = 1.0014205519217815e10, errorEstimate = 8.244928333384438e9, evaluations = 25},Result {result = 8.626254877110512e9, errorEstimate = 1.3879506421073036e9, evaluations = 49},Result {result = 8.616102649530426e9, errorEstimate = 1.01522275800848e7, evaluations = 97},Result {result = 8.616102660994505e9, errorEstimate = 11.464078903198242, evaluations = 193},Result {result = 8.616102660994507e9, errorEstimate = 1.9073486328125e-6, evaluations = 385},Result {result = 8.616102660994507e9, errorEstimate = 1.9073486328125e-6, evaluations = 769}]
*Bug> test2
[Result {result = NaN, errorEstimate = NaN, evaluations = 13},Result {result = NaN, errorEstimate = NaN, evaluations = 25},Result {result = NaN, errorEstimate = NaN, evaluations = 49},Result {result = NaN, errorEstimate = NaN, evaluations = 97},Result {result = NaN, errorEstimate = NaN, evaluations = 193},Result {result = NaN, errorEstimate = NaN, evaluations = 385},Result {result = NaN, errorEstimate = NaN, evaluations = 769}]
*Bug> test3
[Result {result = 8.316940226856241e10, errorEstimate = 8.31634003061991e10, evaluations = 13},Result {result = 4.498200156324874e10, errorEstimate = 3.818740070531367e10, evaluations = 25},Result {result = 4.087451716834612e10, errorEstimate = 4.1074843949026146e9, evaluations = 49},Result {result = 4.0874803200000015e10, errorEstimate = 286031.65388703346, evaluations = 97},Result {result = 4.0874803200000015e10, errorEstimate = 3.5762786865234375e-6, evaluations = 193},Result {result = 4.08748032e10, errorEstimate = 1.0728836059570313e-5, evaluations = 385},Result {result = 4.08748032e10, errorEstimate = 1.0728836059570313e-5, evaluations = 769}]
```

The result should be `gamma(nu/2)*2^((nu-2)/2-1)`:
```
*Bug> import Math.Gamma
*Bug Math.Gamma> gamma(23/2) * 2**((23-2)/2-1)
8.616102660994509e9
*Bug Math.Gamma> gamma(24/2) * 2**((24-2)/2-1)
4.08748032e10
```

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.