Как построить эмпирический Cdf в Matplotlib в Python? - PullRequest
55 голосов
/ 09 июля 2010

Как я могу построить эмпирический CDF массива чисел в matplotlib в Python? Я ищу в формате cdf аналог функции "Hist" в pylab.

Одна вещь, о которой я могу думать, это:

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

Это правильно, хотя? Есть ли более легкий / лучший способ?

спасибо.

Ответы [ 15 ]

81 голосов
/ 27 июля 2012

Если вам нравится linspace и вы предпочитаете однострочники, вы можете сделать:

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

Учитывая мои вкусы, я почти всегда делаю:

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

Что работает для менядаже если есть >O(1e6) значений данных.Если вам действительно нужно уменьшить выборку, я бы установил

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

Редактировать , чтобы ответить на комментарий / редактировать, почему я использую endpoint=False или y, как определено выше.Ниже приведены некоторые технические детали.

Эмпирический CDF обычно формально определяется как

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

, чтобы точно соответствовать этому формальному определению, вам нужно будет использовать y = np.arange(1,len(x)+1)/float(len(x)), чтобы мы получилиy = [1/N, 2/N ... 1].Эта оценка является объективной оценкой, которая будет сходиться к истинному CDF в пределе бесконечных выборок Ссылка на Википедию .

Я склонен использовать y = [0, 1/N, 2/N ... (N-1)/N], поскольку (а) легчеcode / more idomatic, (b), но все еще формально оправдано, поскольку можно всегда поменять CDF(x) на 1-CDF(x) в доказательстве сходимости, и (c) работает с (простым) методом понижающей дискретизации, описанным выше.

В некоторых частных случаях полезно определить

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

, который является промежуточным между этими двумя соглашениями.Который, по сути, говорит, что «есть вероятность 1/(2N) значения меньше, чем наименьшее значение, которое я видел в моей выборке, и 1/(2N) шанс значения больше, чем наибольшее значение, которое я видел до сих пор»..

Однако для больших выборок и разумных распределений соглашение, приведенное в основной части ответа, легко написать, является объективной оценкой истинного CDF и работает с методологией понижающей дискретизации.

69 голосов
/ 11 июля 2010

Вы можете использовать функцию ECDF из библиотеки scikits.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)

С версией 0.4 scicits.statsmodels был переименован в statsmodels. ECDF теперь находится в модуле distributions (тогда как statsmodels.tools.tools.ECDF не рекомендуется).

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()
16 голосов
/ 09 июля 2010

Это выглядит (почти) именно то, что вы хотите. Две вещи:

Во-первых, результаты представляют собой набор из четырех элементов. Третий размер бункеров. Второй является отправной точкой самого маленького мусорного ведра. Первый - это количество точек в или ниже каждой ячейки. (Последнее - это количество точек, выходящих за пределы, но поскольку вы не установили ни одного, все точки будут помечены.)

Во-вторых, вам нужно изменить масштаб результатов, чтобы конечное значение равнялось 1, чтобы следовать обычным соглашениям CDF, но в остальном это правильно.

Вот что он делает под капотом:

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

Он выполняет гистограмму, а затем производит кумулятивную сумму отсчетов в каждом бине. Таким образом, i-е значение результата - это число значений массива, меньшее или равное максимуму i-го бина. Итак, окончательное значение - это просто размер исходного массива.

Наконец, для построения графика вам нужно использовать начальное значение корзины и размер корзины, чтобы определить, какие значения по оси X вам понадобятся.

Другой вариант - использовать numpy.histogram, который может выполнить нормализацию и вернуть ребра корзины. Вам нужно будет самостоятельно накапливать итоговую сумму.

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:] - верхний край каждой ячейки.)

15 голосов
/ 28 апреля 2011

Вы пробовали кумулятивный = истинный аргумент pyplot.hist?

7 голосов
/ 17 апреля 2016

Однострочник на основе ответа Дэйва:

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

Редактировать: hans_meine также предложил это в комментариях.

3 голосов
/ 24 декабря 2014

Мы можем просто использовать функцию step из matplotlib, которая создает пошаговый график, который является определением эмпирического CDF:

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)

Последняя вертикальная линия на max(data) была добавлена ​​вручную. В противном случае сюжет просто останавливается на уровне 1 - 1/len(data).

В качестве альтернативы мы можем использовать опцию where='post' для step()

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

В этом случае начальная вертикальная линия от нуля не отображается.

3 голосов
/ 29 мая 2013

Если вы хотите отобразить фактический истинный ECDF (который, как заметил Дэвид Б, является пошаговой функцией, которая увеличивает 1 / n в каждой из n точек данных), я предлагаю написать код для генерации двух точек «заговора» для каждой точки данных :

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)

Таким образом, вы получите график с n шагами, характерными для ECDF, что особенно удобно для наборов данных, которые достаточно малы, чтобы шаги были видимыми. Кроме того, нет необходимости выполнять какое-либо объединение с гистограммами (что может привести к смещению в нарисованном ECDF).

3 голосов
/ 14 июля 2010

У меня есть тривиальное дополнение к методу AFoglia, чтобы нормализовать 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

Нормализация гисто делает его целым единицей, что означает, что cdf не будет нормализованВы должны сами масштабировать его.

3 голосов
/ 09 июля 2010

Что вы хотите сделать с CDF? Для начала это начало. Вы можете попробовать несколько разных значений, например:

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

Гистограмма перечисляет различные правила для количества бинов, например, num_bins ~ sqrt( len(a) ).

(мелкий шрифт: здесь происходят две совершенно разные вещи,

  • биннинг / гистограмма необработанных данных
  • plot интерполирует плавную кривую через, скажем, 20 значений.

Любой из этих способов может оказаться слишком сложным для данных, которые являются «клочковыми» или имеет длинные хвосты, даже для 1d данных - 2d, 3d данные становятся все труднее.
Смотрите также Density_estimation а также с использованием оценки плотности ядра scipy gaussian ).

2 голосов
/ 17 февраля 2017

Это однострочник в морском заливе с использованием параметра cumulative = True. Вот, пожалуйста,

import seaborn as sns
sns.kdeplot(a, cumulative=True)
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...