max_max_mir
max_max_mir

Reputation: 1736

What is probability density function in the context of scipy.stats.norm?

This is a very basic question, but I can't seem to find a good answer. What exactly does scipy calculate for

scipy.stats.norm(50,10).pdf(45)

I understand that the probability of a particular value like 45 in a gaussian with mean 50 and std dev 10 is 0. So what exactly is pdf calculating? Is it the area under the gaussian curve, and if so, what is the range of values on the x axis?

Upvotes: 12

Views: 17514

Answers (2)

Stefan Zobel
Stefan Zobel

Reputation: 3212

The probability density function of the normal distribution expressed in Python is

from math import pi
from math import exp
from scipy import stats


def normal_pdf(x, mu, sigma):
    return 1.0 / (sigma * (2.0 * pi)**(1/2)) * exp(-1.0 * (x - mu)**2 / (2.0 * (sigma**2)))

(compare that to the wikipedia definition). And this is exactly what scipy.stats.norm().pdf() computes: the value of the pdf at point x for a given mu, sigma.

Note that this is not a probability (= area under the pdf) but rather the value of the pdf at the point x you pass to pdf(x) (and that value can very well be greater than 1.0!). You can see that, for example, for N(0, 0.1) at x = 0:

val = stats.norm(0, 0.1).pdf(0)

print(val)

val = normal_pdf(0, 0, 0.1)

print(val)

which gives the output

3.98942280401

3.989422804014327

Not at all a probability = area under the curve!

Note that this doesn't contradict the statement that the probability of particular value like x = 0 is 0 because, formally, the area under the pdf for a point (i.e., an interval of length 0) is zero (if f is a continuous function on [a, b] and F is its antiderivative on [a, b], then the definite integral of f over [a, b] = F(a) - F(b). Here, a = b = x hence the value of the integral is F(x) - F(x) = 0).

Upvotes: 16

suvy
suvy

Reputation: 693

what you are getting is pdf at value x for a normal pdf function with mean 50 and standard deviation 10. check the function here)

easy to visualize using

npdf=norm(50,10)
plt.plot(range(0,100), npdf.pdf(range(0,100)), 'k-', lw=2)`

you could also generate random variables from the normal pdf you created using

npdf.rvs(1000) #1000 numbers 
hist=plt.hist(n.rvs(10000),bins=100,normed=True)

theoretical pdf and normalized histogram from random variables

Upvotes: 1

Related Questions