scipy.integrate.quad gives wrong result on large r

Question

I am trying to integrate over the sum of two 'half' normal distributions. scipy.integrate.quad works fine when I try to integrate over a small range but returns 0 when I do it for large ranges. Here's 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?

Answer 1

The reason here is that your function is only very strongly peaked in a very small region of your integration region and is effectively zero everywhere else, quad never finds this peak and thus only see's the integrand being zero.

Since in this case you know where the peaks are, it would be reasonable to split the limits of the integration so that you consider the regions around the peaks separately.

To do this you can use the points argument in a slightly bastardized way to force quad to consider the peaks separately.

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