Failure of an integral from 0 to infinity
- 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.