How does this code find the CDF?

Question

How does the below code give CDF? Can someone please explain what
np.arange(len(sorted_data))/float(len(sorted_data)-1) does?

import numpy as np

import matplotlib.pyplot as plt

data = np.loadtxt('Filename.txt')

sorted_data = np.sort(data)

yvals=np.arange(len(sorted_data))/float(len(sorted_data)-1)

plt.plot(sorted_data,yvals)

plt.show()

Answer 1

The code you posted implements the empirical cumulative distribution function. We want to estimate unknown cumulative distribution function, $P(X \le x)$, from the data. The function applies the most basic definition of probability, based on counts

$$ P(\text{something}) = \frac{\text{how often something happen}}{\text{sample size}} $$

In this case, we are are calculating cumulative probabilities, so it is

$$ P(X \le x) = \frac{\#(X \le x)}{N} $$

where $\#$ is the count operator, that counts the values in $X$ meeting the condition.

The function however is incorrect, as it should be rather

yvals = np.arange(1, len(sorted_data)+1)/float(len(sorted_data))

because for the first element of sorted_data we don't want it to return probability equal to zero (obviously we observed this value, so it cannot happen with zero probability!).

But let's discuss it in more details. We are dealing with sorted_data, the data that is already sorted for us in ascending order. For each element in sorted_data, we assign some value from yvals. To calculate the values, we use the np.arange(1, len(sorted_data)+1) function, that returns series $1,2,\dots,N$. This series corresponds to empirical cumulative counts, since first value in sorted_data appeared only once, there was two values smaller or equal to the second value from sorted_data, all the values were smaller or equal then the last value from sorted_data etc. Next, the values are divided by float(len(sorted_data)), so by the sample size $N$ (the float function is needed only for python 2, for it to calculate the floating-point division rather then integer). So we end up with fractions $\tfrac{1}{N},\tfrac{2}{N},\dots,\tfrac{N}{N}$, that correspond to probabilities.

Notice that this assumes that we are dealing with continuous variable in here, as the function ignores the fact that there could be duplicated elements in sorted_data. If you are dealing with discrete variable, you should rather calculate this on counts of sorted values.

If you wonder if we could use similar function to calculate the non-cumulative probabilities, $P(X=x)$, then the answer is: yes, but only for discrete variable (and large sample size to get reasonable approximation). For continuous variable it gets more complicated, because there is infinitely many possible values of $X$, so $P(X=x)=0$, and for this reason we have probability densities, i.e. "probabilities per foot". What follows, for continuous variables we have histograms, kernel density estimators and a few other, less popular, estimators.

Answer 2

When we learn a programming language, one of the skill that I think is crucial is to print your intermediate step for small cases.

Suppose your sorted_data is of length $5$

np.arange(len(sorted_data))

created the numpy array

array([0, 1, 2, 3, 4])

We want to normalize these to uniformly separated points from $0$ to $1$. So we want to divide it by $4.0$ (notice the $.0$, it matters in programming language which data type are you using).

Hence for the example of an array of length $5$, you obtain the following.

array([0.  , 0.25, 0.5 , 0.75, 1.  ])

I will leave it as an exploratary task for you to figure out if it is really a cdf if your data array is something like $[1,2,3,4,5]$, $[1,4,5,8,9]$, and also $[1,2,2,2,5]$.