Python 3.4 Scipy Integrate.quad Dropoff

I'm trying to compute the integral of a Gaussian in python like so: from math import exp from scipy import stats, integrate import scipy.interpolate as interpolate from numpy impor

Solution 1:

I think the key to understand this problem is to recall that numerical integration methods calculate a weighted sum of function values at specific knots.

The gaussian quickly goes to zero as you deviate from the mean, so basically on the interval between (-50, 50) most of the function values are zero. If the integration method fails to sample points from your small area where the function is non-zero, it will see the whole function as completely flat and thus gives you the integral 0.

So what can you do?

Instead of choosing the fixed interval of (-50,50), choose an interval based on smaller sigma values, to avoid integrating over an overly large interval of zeros.

If you go only 5, 10 or 20 standard deviations to the left and to the right, you will not see this issue, and you still have a very accurate integration result.

This is the result if you integrate from 10 standard deviations to the left and to the right.