Solved – Problem with integrating exponential function in R

numerical integrationr

I'm trying to calculate the finite integral for the CDF of the normal distribution, when I plug the equation into Wolfram Alpha and integrate

exp((-t^2)/2) dt from -inf to 1000000

I get: 2.5 (for just the integral), when I run the same equation through R I get 0. On further digging I think the discrepency lies with the 'exp' function. I'm guessing that means I'm doing something wrong or I missed a step somewhere.

Any help would be greatly appreciated.

Best Answer

It actually is because numerical integration has difficulty when only a relatively small part of the space between the ranges is substantially different from zero; it has trouble finding that portion of the range to concentrate on. You can see this by increasing the right side bound; on my system it works fine up to 19.32 but then has trouble.

> integrate(f, -Inf, 19.32)
2.506628 with absolute error < 3.7e-05
> integrate(f, -Inf, 19.33)
6.43579e-05 with absolute error < 0.00012

But why not use pnorm?

> pnorm(1e6)*sqrt(2*pi)
[1] 2.506628
Related Question