CODEWITHSUNDEEP
pythonrnumpyscipynumerical-integration
davidhalpern
davidhalpern

Reputation: 31

Numerical integration over mean of normal in scipy

I'm getting some weird output from the scipy integrate.quad function when integrating a normal pdf. Here is the code I'm trying to use:

inpdf = lambda c: norm.pdf(50, loc=c, scale = 1)
result = integrate.quad(inpdf, -np.inf, np.inf)
print result

Which returns (3.281718223506531e-99, 0.0) which is obviously wrong. When I change the x value to 0, in the snippet below, I get appropriate output:

inpdf = lambda c: norm.pdf(0, loc=c, scale = 1)
result = integrate.quad(inpdf, -np.inf, np.inf)
print result

returns (0.9999999999999998, 1.0178191320905743e-08) So pretty close to 1 which is what it should be. The exact same thing happens using R's integrate function so maybe this is just a known property of the quadrature algorithm? Does anyone know why this is the case?

Upvotes: 2

Views: 301

Answers (1)

Jonathan Lisic
Jonathan Lisic

Reputation: 1840

At least in R, and probably in scipy this this issue is fairly straightforward.

What you are asking the program to do is to integrate a function over the entire real line. The quadrature algorithm tries to find where the the function is non-zero which in the case of a normal pdf is everywhere, but the significant portions are hard to find.

In the case of R the integrate function evaluates at a number of points, 100 by default, and more than likely they won't occur in the small interval that your function is reasonably large. If instead you tried something like

> integrate(dnorm,lower=45,upper=55,mean=50)
0.9999994 with absolute error < 8.7e-10

you get the right result or close enough. If you move further away the quadrature algorithm starts failing again

> integrate(dnorm,lower=-1000,upper=1000,mean=50)
0 with absolute error < 0

Upvotes: 5

Related Questions