Reputation: 7338
I have an array of random samples from a normal distribution where I want to evaluate the CDF of each element in place
import numpy as np
arr = np.random.normal(0, 1, 10000)
arr
array([-0.03960733, -0.58329607, -1.55133923, ..., -0.94473672,
1.24757701, -0.66197476])
I know I could do this using scipy.stats.norm().cdf
, but I am restricted to only using numpy.
I have found this SO post that outlines how to do something similar using numpy.histogram
and numpy.cumsum
. How can I extend this (using only numpy) to evaluate the CDF of each element in place, so the output array is as below
from scipy import stats
stats.norm().cdf(arr)
array([0.48420309, 0.279847 , 0.06041021, ..., 0.17239665, 0.893907 ,
0.2539937 ])
Upvotes: 1
Views: 273
Reputation: 7338
It seems this can be achieved using numpy.argsort()
twice to get the ranks of each random sample in arr
. There is some round-off error however
import numpy as np
arr = np.random.normal(0, 1, 10000)
arr
array([-0.24822623, -0.49071664, -0.75405418, ..., -0.59249804,
-0.9140224 , 0.18904534])
x = arr.argsort().argsort() # ranks of each entry in `arr`
y = np.arange(len(arr)) / len(arr)
numpy_cdfs = y[x] # sort `y` by ranks
numpy_cdfs
array([0.3973, 0.307 , 0.2204, ..., 0.2713, 0.1745, 0.5696])
If we compare to scipy, we need to set the absolute tolerance to 1e-2 (quite high).
from scipy import stats
scipy_cdfs = stats.norm().cdf(arr)
scipy_cdfs
array([0.40197969, 0.31181344, 0.22540834, ..., 0.27675857, 0.18035254,
0.57497136])
np.allclose(numpy_cdfs, scipy_cdfs, atol=1e-2)
True
This error will decrease the more samples we have.
Upvotes: 1