invictus
invictus

Reputation: 1971

how to integrate normal distribution with numpy and scipy?

trying to integrate a normal distribution

import numpy as np
import scipy

this does not work, throws a type error:

def f_np(x):
  return np.random.normal(size=x)

integrate.quad(f_manual, -np.inf, np.inf)

mtrand.pyx in numpy.random.mtrand.RandomState.normal()

_common.pyx in numpy.random._common.cont()

TypeError: 'float' object cannot be interpreted as an integer

but it works if I manually enter the pdf:

def f_manual(x):
   return 1/(1 * np.sqrt(2 * np.pi)) * np.exp(-(x - 0)**2 / (2 * 1**2) )

integrate.quad(f_manual, -np.inf, np.inf)

(0.9999999999999997, 1.017819145094224e-08)

any idea why?

Upvotes: 2

Views: 1693

Answers (1)

user650654
user650654

Reputation: 6128

Several things going on here.

np.random.normal draws samples from the normal distribution. The size parameter specifies the number of samples you want. If you specify 10 you'll get an array with 10 samples. If you specify a tuple, like (4, 5) you'll get a 4x5 array. Also, np.inf is a float and np.random.normal is expecting an integer or a tuple of integers for the size parameter.

What you have in f_manual is a deterministic function (i.e. the PDF) which returns the value of the PDF at the value x.

These are two different things.

scipy has a function to return the PDF of a Gaussian: scipy.stats.norm.pdf

import scipy.stats
scipy.integrate.quad(scipy.stats.norm.pdf, -np.inf, np.inf)
# (0.9999999999999998, 1.0178191320905743e-08)

scipy also has a CDF function that returns the integral from -inf to x:

scipy.stats.norm.cdf(np.inf)
# 1.0

scipy.stats.norm.cdf(0)
# 0.5

Upvotes: 6

Related Questions