Geek Answers Handbook

Scipy.integrate.quad gives the wrong result on large ranges

I am trying to integrate over the sum of two "half" normal distributions. scipy.integrate.quadworks great when I try to integrate into a small range, but returns 0 when I do this for large ranges. Here is the code:

mu1 = 0
mu2 = 0
std1 = 1
std2 = 1

def integral_fun(x):
    nor1 = 0.5 * ((1 / (np.sqrt(2 * np.pi) * std1)) * (np.e ** ((-(x-mu1) ** 2) / (2 * std1 **2))))
    nor2 = 0.5 * ((1 / (np.sqrt(2 * np.pi) * std2)) * (np.e ** ((-(x-mu2) ** 2) / (2 * std2 **2))))
    return nor1 + nor2


integrate.quad(integral_fun, -5, 5)
Out[54]: (0.9999994266968564, 8.668320228277793e-10)

integrate.quad(integral_fun, -10, 10)
Out[55]: (1.0000000000000002, 8.671029607900576e-10)

integrate.quad(integral_fun, -100000, 100000)
Out[56]: (0.0, 0.0)

Why is this happening?

+4
source share
1 answer

The reason is that your function peaked very strongly in a very small area of ​​your integration area and is virtually zero everywhere, quadnever finds this peak and, therefore, only sees that the integrand is null.

, , , .

, points , quad .

In [3]: integrate.quad(integral_fun, -100000, 100000, points=[-10,10])
Out[3]: (1.0000000000000002, 8.671029607900576e-10)
+3

More articles:


All Articles