[SciPy-User] integrate.quad(stats.norm(loc=100.0).pdf, -np.inf, np.inf) = 0 ?

[Daniele Nicolodi]

Hello,

sorry for the cryptic Subject, but I didn't know how else to title this.

I don't have much experience with numerical integration, so it may very
well be that I'm doing something stupid, please advise :-)

I need to compute some numerical integrals of some functions that are
very "localized" in their domain (trying to compute some posterior
probability density).  A Gaussian (normal) function is a good example.
I find that if the region in which the function is significantly non
zero is away from the center of the integration interval, the
integration routines "step over" that domain interval and do not "see"
the function having values different from zero.

For example, integrating a Gaussian centered around zero, works as
expected:

integrate.quad(stats.norm(loc=0.0).pdf, -np.inf, np.inf)
(0.9999999999999998, 1.0178191320905743e-08)

but if I shift the Gaussian such that the function value at 0 is
numerically zero, the integration fails:

integrate.quad(stats.norm(loc=100.0).pdf, -np.inf, np.inf)
(0.0, 0.0)

same if I have an integration interval not centered:

integrate.quad(stats.norm(loc=0.0).pdf, -100, 1000)
(4.471689909323402e-30, 8.454234364421014e-30)

vs:

integrate.quad(stats.norm(loc=0.0).pdf, -100, 100)
(1.0000000000000002, 1.0346447361664605e-12)

I definitely see why those examples are problematic, but I don't know
what I can do about it other than evaluating the functions on a grid to
restrict somehow the integration limits.

Does anyone have any suggestion?

Thanks! Cheers,
Daniele