1

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 ])
PyRsquared
  • 6,970
  • 11
  • 50
  • 86

1 Answers1

1

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.

PyRsquared
  • 6,970
  • 11
  • 50
  • 86