Comment tracer le cdf empirique (ecdf) ?

Question

Comment puis-je tracer la FCD empirique d'un tableau de nombres dans matplotlib en Python ? Je cherche l'analogue cdf de la fonction "hist" de pylab.

Une chose à laquelle je peux penser est :

from scipy.stats import cumfreq
a = array([...]) # my array of numbers
num_bins =  20
b = cumfreq(a, num_bins)
plt.plot(b)

Answer 1

Si vous aimez linspace et préférez les one-liners, vous pouvez le faire :

plt.plot(np.sort(a), np.linspace(0, 1, len(a), endpoint=False))

Étant donné mes goûts, je le fais presque toujours :

# a is the data array
x = np.sort(a)
y = np.arange(len(x))/float(len(x))
plt.plot(x, y)

Ce qui fonctionne pour moi même s'il y a >O(1e6) les valeurs des données. Si vous avez vraiment besoin de réduire l'échantillonage, je mettrais

x = np.sort(a)[::down_sampling_step]

Modifier pour répondre au commentaire/rédaction sur la raison pour laquelle j'utilise endpoint=False ou le y comme défini ci-dessus. Voici quelques détails techniques.

La CDF empirique est généralement définie formellement comme suit

CDF(x) = "number of samples <= x"/"number of samples"

afin de correspondre exactement à cette définition formelle, il faudrait utiliser y = np.arange(1,len(x)+1)/float(len(x)) de sorte que nous obtenions y = [1/N, 2/N ... 1] . Cet estimateur est un estimateur sans biais qui convergera vers la vraie CDF dans la limite d'échantillons infinis. Wikipédia réf. .

J'ai tendance à utiliser y = [0, 1/N, 2/N ... (N-1)/N] car (a) il est plus facile à coder/plus idiomatique, (b) mais est toujours formellement justifié puisqu'on peut toujours échanger CDF(x) con 1-CDF(x) dans la preuve de convergence, et (c) fonctionne avec la méthode de sous-échantillonnage (facile) décrite ci-dessus.

Dans certains cas particuliers, il est utile de définir

y = (arange(len(x))+0.5)/len(x)

qui est intermédiaire entre ces deux conventions. Ce qui, en fait, dit "il y a une 1/(2N) chance d'avoir une valeur inférieure à la plus basse que j'ai vue dans mon échantillon, et une 1/(2N) chance d'une valeur supérieure à la plus grande que j'ai vue jusqu'à présent.

Notez que la sélection de cette convention interagit avec l'option where utilisé dans le plt.step s'il semble plus utile d'afficher la FCD comme une fonction constante par morceaux. Afin de correspondre exactement à la définition formelle mentionnée ci-dessus, il faudrait utiliser where=pre la proposition y=[0,1/N..., 1-1/N] convention, ou where=post avec le y=[1/N, 2/N ... 1] convention, mais pas l'inverse.

Cependant, pour de grands échantillons et des distributions raisonnables, la convention donnée dans le corps de la réponse est facile à écrire, est un estimateur sans biais de la vraie CDF et fonctionne avec la méthodologie de sous-échantillonnage.

Answer 2

Vous pouvez utiliser le ECDF de la fonction scikits.statsmodels bibliothèque :

import numpy as np
import scikits.statsmodels as sm
import matplotlib.pyplot as plt

sample = np.random.uniform(0, 1, 50)
ecdf = sm.tools.ECDF(sample)

x = np.linspace(min(sample), max(sample))
y = ecdf(x)
plt.step(x, y)

Avec la version 0.4 scicits.statsmodels a été renommé en statsmodels . ECDF est maintenant situé dans le distributions (alors que statsmodels.tools.tools.ECDF est amortie).

import numpy as np
import statsmodels.api as sm # recommended import according to the docs
import matplotlib.pyplot as plt

sample = np.random.uniform(0, 1, 50)
ecdf = sm.distributions.ECDF(sample)

x = np.linspace(min(sample), max(sample))
y = ecdf(x)
plt.step(x, y)
plt.show()

Answer 3

Cela semble être (presque) exactement ce que vous voulez. Deux choses :

Tout d'abord, les résultats sont un tuple de quatre éléments. Le troisième est la taille des bacs. Le deuxième est le point de départ du plus petit bac. Le premier est le nombre de points dans les limites ou en dessous de chaque bac. (Le dernier est le nombre de points en dehors des limites, mais comme vous n'en avez pas défini, tous les points seront mis en bin).

Deuxièmement, vous devrez modifier l'échelle des résultats pour que la valeur finale soit égale à 1, afin de respecter les conventions habituelles d'une FCD, mais sinon, c'est correct.

Voici ce qu'il fait sous le capot :

def cumfreq(a, numbins=10, defaultreallimits=None):
    # docstring omitted
    h,l,b,e = histogram(a,numbins,defaultreallimits)
    cumhist = np.cumsum(h*1, axis=0)
    return cumhist,l,b,e

Il fait l'histogramme, puis produit une somme cumulative des comptes dans chaque case. La ième valeur du résultat est donc le nombre de valeurs du tableau inférieures ou égales au maximum de la ième case. Ainsi, la valeur finale est juste la taille du tableau initial.

Enfin, pour le tracer, vous devrez utiliser la valeur initiale de l'emplacement et la taille de l'emplacement pour déterminer les valeurs de l'axe des x dont vous aurez besoin.

Une autre option consiste à utiliser numpy.histogram qui peut effectuer la normalisation et renvoie les bords des bacs. Vous devrez faire la somme cumulative des comptes résultants vous-même.

a = array([...]) # your array of numbers
num_bins = 20
counts, bin_edges = numpy.histogram(a, bins=num_bins, normed=True)
cdf = numpy.cumsum(counts)
pylab.plot(bin_edges[1:], cdf)

( bin_edges[1:] est le bord supérieur de chaque bac).

Answer 4

Avez-vous essayé l'argument cumulative=True de pyplot.hist ?

Answer 5

Que voulez-vous faire avec le CDF ? La tracer, c'est un début. Vous pouvez essayer plusieurs valeurs différentes, comme ceci :

from __future__ import division
import numpy as np
from scipy.stats import cumfreq
import pylab as plt

hi = 100.
a = np.arange(hi) ** 2
for nbins in ( 2, 20, 100 ):
    cf = cumfreq(a, nbins)  # bin values, lowerlimit, binsize, extrapoints
    w = hi / nbins
    x = np.linspace( w/2, hi - w/2, nbins )  # care
    # print x, cf
    plt.plot( x, cf[0], label=str(nbins) )

plt.legend()
plt.show()

Histogramme énumère diverses règles pour le nombre de bacs, par ex. num_bins ~ sqrt( len(a) ) .

(En petits caractères : deux choses bien différentes se passent ici,

• binning / histogramme des données brutes
• plot interpole une courbe lisse à travers les 20 valeurs binées.

L'un ou l'autre peut s'écarter considérablement des données qui sont "touffues". ou qui ont de longues queues, même pour des données 1d - les données 2d, 3d deviennent de plus en plus difficiles.
Voir aussi Estimation de la densité et utilisation de l'estimation de la densité du noyau gaussien de scipy ).