A quick look in the source code shows that scipy.stats.norm.pdf
simply returns the value for x
of the pdf using NumPy:
def _norm_pdf(x):
return np.exp(-x**2/2.0) / _norm_pdf_C
where _norm_pdf_C = np.sqrt(2*np.pi)
.
For the cdf, since we talk of a normal distribution, special functions are used (for the relation between them and the normal distribution, see here).
SciPy implements special functions directly in C. In particular, the cumulative distribution function is computed from ndtr.c
. So, even if NumPy is really fast, C is still faster in this case, I assume.
EDIT
Sorry, I just realized my answer does not completely answer your question.
First of all, NumPy also implements mathematical operations in C.
Therefore, to understand why the difference in the times, one should understand what happens in C.
- If you look at this question, it seems that numerical values and hardware architecture can affect the time.
So I checked again the C implementation for the cdf, and I saw that constants and coefficients of the polynomials that evaluate the special functions are not computed but stored in arrays and variables! For example, 1/sqrt(2)
is contained in NPY_SQRT1_2
. This could be the reason why cdf is faster than pdf!
Thus I tried to calculate the pdf having the constants already initialized:
import scipy.stats as st
from datetime import datetime
import numpy as np
num_iter = 100000
x_lower = 0.25
x_upper = 0.75
const = np.sqrt(2*np.pi)
time_start = datetime.now()
for x in np.arange(x_lower, x_upper, (x_upper - x_lower) / (num_iter - 1)):
# y = st.norm.pdf(x)
y = np.exp((x**2 / 2)) / const
time_end = datetime.now()
print(time_end - time_start)
time_start = datetime.now()
for x in np.arange(x_lower, x_upper, (x_upper - x_lower) / (num_iter - 1)):
y = st.norm.cdf(x)
time_end = datetime.now()
This code gives me:
0:00:00.202531
0:00:07.703083
Notice that also norm.pdf
pre-initializes the denominator of the pdf, but in the for loop you are calling the method every time, slowing things down.
P.S.: if you try to get rid of the loops in your original code and have simply x = np.arange(x_lower, x_upper, (x_upper - x_lower) / (num_iter - 1))
, cdf is again faster. The reason can be that cdf is computed with polynomials approximations. But I did not find information on how exactly C handles the exponential to get a comparison.