[SciPy-User] stats.norm and integrate.quad

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

[SciPy-User] stats.norm and integrate.quad

Andrew Nelson
Hi all,
I'm trying to integrate stats.norm.pdf in the range [-np.inf, x] using integrate.quad. I'm having a few issues where the integral dives to zero when it should be one - https://gist.github.com/andyfaff/c8a94ab87333c42791e58e7d6e12a05e

It does this when norm.pdf(x) reaches the smallest number expressible by a float, and for some unknown reason when x ~ 20.

Is there anyway I can mitigate against this behaviour using `integrate.quad` options? I'm hoping this is expected behaviour.

cheers,
Andrew.


p.s.  I'm using integration instead of norm.cdf for a specific reason.

--
_____________________________________
Dr. Andrew Nelson


_____________________________________

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: stats.norm and integrate.quad

josef.pktd


On Tue, Feb 6, 2018 at 6:01 PM, Andrew Nelson <[hidden email]> wrote:
Hi all,
I'm trying to integrate stats.norm.pdf in the range [-np.inf, x] using integrate.quad. I'm having a few issues where the integral dives to zero when it should be one - https://gist.github.com/andyfaff/c8a94ab87333c42791e58e7d6e12a05e

It does this when norm.pdf(x) reaches the smallest number expressible by a float, and for some unknown reason when x ~ 20.

Is there anyway I can mitigate against this behaviour using `integrate.quad` options? I'm hoping this is expected behaviour.

AFAIR, I never figured out the quad options for this (there are no gaussian weights)

Things that I tried in the past:

use integration limit using information about where tails are close to zero 

>>> from scipy import stats
>>> stats.norm.ppf(1e-10)
-6.3613409024040557
>>> stats.norm.ppf(1e-15)
-7.9413453261709979
>>> stats.norm.isf(1e-15)
7.9413453261709979

break up in sub intervals with known central region (I only tried the points option iin quad for a few bad cases, AFAIR)

>>> from scipy.integrate import quad
>>> quad(stats.norm.pdf, -np.inf, 40)
(5.505747798875221e-28, 1.0434548229074602e-27)

>>> quad(stats.norm.pdf, -np.inf, 0)
(0.4999999999999999, 5.089095674629994e-09)
>>> quad(stats.norm.pdf, 0, 40)
(0.5000000000000001, 4.335515230512748e-10)


use more symmetric bounds so integrate quad looks at where the spike or action is

>>> quad(stats.norm.pdf, -40, 40)
(1.0000000000000002, 8.671029619552417e-10)


use an indicator function and integrate over full range -inf, inf

>>> quad(lambda x: stats.norm.pdf(x) * (x <= 40), -np.inf, np.inf)
(0.9999999999999998, 1.0178191349259989e-08)
>>> quad(lambda x: stats.norm.pdf(x) * (x <= 0), -np.inf, np.inf)
(0.4999999999999999, 5.089095674629994e-09)


Josef

 

cheers,
Andrew.


p.s.  I'm using integration instead of norm.cdf for a specific reason.

--
_____________________________________
Dr. Andrew Nelson


_____________________________________

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user



_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user