5

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?

Srivatsan
  • 9,225
  • 13
  • 58
  • 83
Ankur Ankan
  • 2,953
  • 2
  • 23
  • 38
  • numerical integrals are approximated by a weighted sum of function values. If I had to guess I would say that the function values are all 0, because almost all points of the range have a density of (almost) 0. – cel Jun 18 '15 at 11:21
  • it works though if the ranges are from `-inf` to `inf` – Srivatsan Jun 18 '15 at 11:23
  • 3
    @ThePredator This is because in that case scipy will transform the integral into one over the range ``(0,1]`` with the transformation ``x = (1-t)/t``. This then well samples points near the origin. – Simon Gibbons Jun 18 '15 at 11:27

1 Answers1

5

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)
Simon Gibbons
  • 6,969
  • 1
  • 21
  • 34