2010-07-09 14 views
51

¿Cómo puedo trazar la CDF empírica de una matriz de números en matplotlib en Python? Estoy buscando el análogo cdf de la función "hist" de pylab.Cómo trazar pldf empírico en matplotlib en Python?

Una cosa que puedo pensar es:

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

¿Eso es correcto, aunque? ¿Hay una manera más fácil/mejor?

gracias.

Respuesta

15

Eso parece ser (casi) exactamente lo que quiere. Dos cosas:

En primer lugar, los resultados son una tupla de cuatro elementos. El tercero es el tamaño de los contenedores. El segundo es el punto de inicio del bin más pequeño. El primero es la cantidad de puntos en el interior o debajo de cada contenedor. (El último es el número de puntos fuera de los límites, pero como no ha configurado ninguno, todos los puntos se agruparán).

En segundo lugar, querrá volver a escalar los resultados para que el valor final sea 1, a siga las convenciones usuales de un CDF, pero de lo contrario es correcto.

Esto es lo que lo hace bajo el capó:

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 

Se hace el histografía, produce entonces una suma acumulada de los recuentos en cada bandeja. Por lo tanto, el valor ith del resultado es el número de valores de matriz menores o iguales al máximo de la i-ésima categoría. Entonces, el valor final es solo el tamaño de la matriz inicial.

Finalmente, para trazarlo, necesitará usar el valor inicial de la bandeja y el tamaño de la bandeja para determinar qué valores de eje x necesitará.

Otra opción es usar numpy.histogram que puede hacer la normalización y devuelve los bordes de la bandeja. Tendrá que hacer la suma acumulativa de los recuentos resultantes usted mismo.

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:] es el borde superior de cada bandeja.)

+17

: este código en realidad no le dan la CDF empírica (una función escalonada aumentando en 1/n en cada uno de los n puntos de datos). En cambio, este código proporciona una estimación de la CDF basada en una estimación del PDF basada en un histograma. Esta estimación basada en el histograma puede ser manipulada/sesgada mediante la selección cuidadosa/inadecuada de los contenedores, por lo que no es una caracterización tan buena de la verdadera CDF como la ECDF real. –

+2

También me desagrada el punto de que esto impone el binning; vea la respuesta corta de Dave, que simplemente usa 'numpy.sort' para trazar el CDF sin binning. –

3

¿Qué desea hacer con el CDF? Para trazarlo, eso es un comienzo. Usted podría tratar de unos valores diferentes, así:

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() 

Histogram enumera diversas reglas para el número de contenedores, por ejemplo, num_bins ~ sqrt(len(a)).

(Letra pequeña: dos cosas muy diferentes están pasando aquí,

  • de agrupación/histografía los datos en bruto
  • plot interpola una curva suave a través de los 20 valores agrupados decir

. Cualquiera de estos puede ir demasiado lejos en datos que son "grumosos" o tiene colas largas, incluso para datos 1d - 2d, los datos en 3D se vuelven cada vez más difíciles.
Véase también Density_estimation y using scipy gaussian kernel density estimation ).

65

Usted puede utilizar la función de la biblioteca ECDFscikits.statsmodels:

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) 

Con la versión 0.4 scicits.statsmodels se renombró a statsmodels. ECDF ahora se encuentra en el módulo distributions (mientras que statsmodels.tools.tools.ECDF está depreciado).

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() 
+2

@bmu (y @Luca): impresionante; ¡Gracias por amablemente hacer que el código sea actual con el modelo actual de stats! – ars

+0

Para scikits.statsmodels v0.3.1 tuve que 'importar scikits.statsmodels.tools como smtools' y' ecdf = smtools.tools.EDCF (...) ' – alexei

3

I tienen una adición trivial con el método de AFoglia, para normalizar la CDF

n_counts,bin_edges = np.histogram(myarray,bins=11,normed=True) 
cdf = np.cumsum(n_counts) # cdf not normalized, despite above 
scale = 1.0/cdf[-1] 
ncdf = scale * cdf 

La normalización de la histo hace su integral unidad, lo que significa la cdf no será normalizada. Tienes que escalarlo tú mismo.

13

¿Has probado el argumento cumulative = True para pyplot.hist?

+1

Muy buen comentario. Aún así, eso impone binning; ver la respuesta de Dave usando np.sort. –

+0

Buena y sencilla opción, pero la desventaja es una personalización limitada de la gráfica de líneas resultante, p. no pude encontrar la forma de agregar marcadores. Fui a la respuesta 'scikits.statsmodels'. – alexei

62

Si te gusta linspace y prefiere una sola línea, que puede hacer:

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

Teniendo en cuenta mis gustos, casi siempre hago:

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

que funciona para mí, incluso si no son >O(1e6) valores de datos. Si realmente necesita la muestra abajo fijaría

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

Editar para responder a comentar/edición el por qué uso o la endpoint=Falseyvals como se definió anteriormente. Los siguientes son algunos detalles técnicos.

La CDF empírica es usualmente definido formalmente como

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

con el fin de coincidir exactamente con esta definición formal de lo que tendría que utilizar yvals = np.arange(1,len(sorted_)+1)/float(len(sorted_)) de modo que consigamos yvals = [1/N, 2/N ... 1]. Este estimador es un estimador insesgado que convergerá al verdadero CDF en el límite de muestras infinitas Wikipedia ref..

que tienden a utilizar yvals = [0, 1/N, 2/N ... (N-1)/N] ya que (a) es más fácil de código/más idiomáticas, (b), pero aún se justifica formalmente desde siempre se puede intercambiar CDF(x) con 1-CDF(x) en la prueba de convergencia, y (c) trabaja con el (fácil) método de reducción de muestreo descrito anteriormente.

En algunos casos particular, es útil definir

yvals = (arange(len(sorted_))+0.5)/len(sorted_) 

que es intermedia entre estas dos convenciones. Que, en efecto, dice "hay una probabilidad de 1/(2N) de un valor inferior al más bajo que he visto en mi muestra, y una probabilidad de 1/(2N) de un valor mayor que el más grande que he visto hasta ahora.

Sin embargo, para muestras grandes y distribuciones razonables, la convención dada en el cuerpo principal de la respuesta es fácil de escribir, es un estimador insesgado del verdadero CDF y funciona con la metodología de reducción de muestreo.

+3

Esta respuesta debería recibir más votos positivos, ya que es la única hasta el momento que no impone binning. Solo simplifiqué el código un poco, usando linspace. –

+1

@hans_meine su edición, es decir, 'yvals = linspace (0,1, len (ordenado))', produce 'yvals' que no son un estimador insesgado de la verdadera CDF. – Dave

+0

Entonces, deberíamos haber usado linspace con 'endpoint = False', ¿verdad? –

3

Si desea visualizar el ECDF verdadero real (que como notó David B es una función de paso que aumenta 1/n en cada uno de los n puntos de datos), mi sugerencia es escribir código para generar dos puntos de "trazado" para cada punto de datos:

a = array([...]) # your array of numbers 
sorted=np.sort(a) 
x2 = [] 
y2 = [] 
y = 0 
for x in sorted: 
    x2.extend([x,x]) 
    y2.append(y) 
    y += 1.0/len(a) 
    y2.append(y) 
plt.plot(x2,y2) 

de esta manera se obtendrá una parcela con los n pasos que son característicos de un ECDF, lo cual es bueno, especialmente para los conjuntos de datos que son lo suficientemente pequeños para los pasos que sean visibles. Además, no hay necesidad de hacer ningún binning con histogramas (lo que corre el riesgo de introducir un sesgo en el ECDF extraído).

2

Sólo puede utilizar la función de stepmatplotlib, lo que hace que una parcela paso a paso, que es la definición de la CDF empírica:

import numpy as np 
from matplotlib import pyplot as plt 

data = np.random.randn(11) 

levels = np.linspace(0, 1, len(data) + 1) # endpoint 1 is included by default 
plt.step(sorted(list(data) + [max(data)]), levels) 

La línea vertical final en max(data) se añadió manualmente. De lo contrario, la trama simplemente se detiene en el nivel 1 - 1/len(data).

alternativa que puede utilizar la opción where='post' a step()

levels = np.linspace(1./len(data), 1, len(data)) 
plt.step(sorted(data), levels, where='post') 

en cuyo caso no se representa la línea vertical inicial de cero.

1

(Esto es una copia de mi respuesta a la pregunta: Plotting CDF of a pandas series in python)

Un CDF o parcela función de distribución acumulada es básicamente un gráfico con el eje X los valores ordenados y en el eje Y el acumulado distribución. Entonces, crearía una nueva serie con los valores ordenados como índice y la distribución acumulativa como valores.

En primer lugar crear una serie ejemplo:

import pandas as pd 
import numpy as np 
ser = pd.Series(np.random.normal(size=100)) 

Ordenar la serie:

ser = ser.order() 

Ahora, antes de continuar, añadir de nuevo el último (y más grande) de valor. Este paso es importante, especialmente para los pequeños tamaños de muestra con el fin de obtener una CDF imparcial:

ser[len(ser)] = ser.iloc[-1] 

crear una nueva serie con los valores ordenados como el índice y la distribución acumulada como valores

cum_dist = np.linspace(0.,1.,len(ser)) 
ser_cdf = pd.Series(cum_dist, index=ser) 

Por último, la trama la función como pasos:

ser_cdf.plot(drawstyle='steps') 
5

de una sola línea sobre la base de la respuesta de Dave:

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

Editar: Esto también fue sugerido por hans_meine en los comentarios.

+1

Esta es la respuesta más directa, que resuelve el problema con elegancia. ¡Esta debería ser la respuesta aceptada! – Alex

1

Se trata de utilizar el bokeh

`` `

from bokeh.plotting import figure, show 
from statsmodels.distributions.empirical_distribution import ECDF 
ecdf = ECDF(pd_series) 
p = figure(title="tests", tools="save", background_fill_color="#E8DDCB") 
p.line(ecdf.x,ecdf.y) 
show(p) 

` ``

1

Suponiendo que Vals sostiene sus valores, a continuación, sólo tiene que trazar la CDF de la siguiente manera:

y = numpy.arange(0, 101) 
x = numpy.percentile(vals, y) 
plot(x, y) 

Para escalarlo entre 0 y 1, simplemente divida y por 100.

0

Es un trazador de líneas en nadado usando el parámetro cumulativo = True. Aquí tiene,

import seaborn as sns 
sns.kdeplot(a, cumulative=True) 
0

Ninguna de las respuestas cubre hasta el momento lo que quería cuando aterricé aquí, que es:

def empirical_cdf(x, data): 
    "evaluate ecdf of data at points x" 
    return np.mean(data[None, :] <= x[:, None], axis=1) 

Se evalúa la CDF empírica de un determinado conjunto de datos en una matriz de puntos x, que no tiene que ser ordenado. No hay binning intermedio y no hay bibliotecas externas.

un método equivalente que se adapta mejor para x grande es para ordenar los datos y utilizar np.searchsorted: Sólo una nota rápida

def empirical_cdf(x, data): 
    "evaluate ecdf of data at points x" 
    data = np.sort(data) 
    return np.searchsorted(data, x)/float(data.size) 
Cuestiones relacionadas