Empirical cdf in python similiar to matlab's one
Question: I have some code in matlab, that I would like to rewrite into python. It's
simple program, that computes some distribution and plot it in double-log
scale.

The problem I occured is with computing cdf. Here is matlab code:

    
    
    for D = 1:10
        delta = D / 10;
        for k = 1:n
            N_delta = poissrnd(delta^-alpha,1);
            Y_k_delta = ( (1 - randn(N_delta)) / (delta.^alpha) ).^(-1/alpha);
            Y_k_delta = Y_k_delta(Y_k_delta > delta);
            X(k) = sum(Y_k_delta);
            %disp(X(k))
    
        end
        [f,x] = ecdf(X);
    
        plot(log(x), log(1-f))
        hold on
    end
    

In matlab one I can simply use:

    
    
    [f,x] = ecdf(X);
    

to get cdf (f) at points x.
[Here](http://www.mathworks.com/help/stats/ecdf.html) is documentation for it.  
In python it is more complicated:

    
    
    import numpy as np
    from scipy.stats import norm
    import matplotlib.pyplot as plt
    from statsmodels.distributions.empirical_distribution import ECDF
    
    alpha = 1.5
    n = 1000
    X = []
    for delta in range(1,5):
        delta = delta/10.0
        for k in range(1,n + 1):
            N_delta = np.random.poisson(delta**(-alpha), 1)
            Y_k_delta = ( (1 - np.random.random(N_delta)) / (delta**alpha) )**(-1/alpha)
            Y_k_delta = [i for i in Y_k_delta if i > delta]
            X.append(np.sum(Y_k_delta))
    
        ecdf = ECDF(X)
    
        x = np.linspace(min(X), max(X))
        f = ecdf(x)
        plt.plot(np.log(f), np.log(1-f))
    
    plt.show()
    

It makes my plot look very strange, definetly not smooth as matlab's one.  
I think the problem is that I do not understand `ECDF` function or it works
differently than in matlab.  
I implemented [this](http://stackoverflow.com/questions/3209362/how-to-plot-
empirical-cdf-in-matplotlib-in-python) solution (the most points one) for my
python code, but it looks like it doesn't work correctly.


Answer: Once you have your sample, you can easily compute the ECDF using a combination
of
[`np.unique`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.unique.html)*
and
[`np.cumsum`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.cumsum.html):

    
    
    import numpy as np
    
    def ecdf(sample):
    
        # convert sample to a numpy array, if it isn't already
        sample = np.atleast_1d(sample)
    
        # find the unique values and their corresponding counts
        quantiles, counts = np.unique(sample, return_counts=True)
    
        # take the cumulative sum of the counts and divide by the sample size to
        # get the cumulative probabilities between 0 and 1
        cumprob = np.cumsum(counts).astype(np.double) / sample.size
    
        return quantiles, cumprob
    

For example:

    
    
    from scipy import stats
    from matplotlib import pyplot as plt
    
    # a normal distribution with a mean of 0 and standard deviation of 1
    n = stats.norm(loc=0, scale=1)
    
    # draw some random samples from it
    sample = n.rvs(100)
    
    # compute the ECDF of the samples
    qe, pe = ecdf(sample)
    
    # evaluate the theoretical CDF over the same range
    q = np.linspace(qe[0], qe[-1], 1000)
    p = n.cdf(q)
    
    # plot
    fig, ax = plt.subplots(1, 1)
    ax.hold(True)
    ax.plot(q, p, '-k', lw=2, label='Theoretical CDF')
    ax.plot(qe, pe, '-r', lw=2, label='Empirical CDF')
    ax.set_xlabel('Quantile')
    ax.set_ylabel('Cumulative probability')
    ax.legend(fancybox=True, loc='right')
    
    plt.show()
    

[![enter image description
here](http://i.stack.imgur.com/92VBM.png)](http://i.stack.imgur.com/92VBM.png)

* * *

* If you're using a version of numpy older than 1.9.0 then `np.unique` won't accept the `return_counts` keyword argument, and you'll get a `TypeError`:
    
    
    TypeError: unique() got an unexpected keyword argument 'return_counts'
    

In that case, a workaround would be to get the set of "inverse" indices and
use `np.bincount` to count the occurrences:

    
    
    quantiles, idx = np.unique(sample, return_inverse=True)
    counts = np.bincount(idx)
    

