如何使用NumPy获取累积分布函数？

Question

我想使用NumPy创建一个累积分布函数(CDF)，我的代码如下：

histo = np.zeros(4096, dtype = np.int32)
for x in range(0, width):
   for y in range(0, height):
      histo[data[x][y]] += 1
      q = 0 
   cdf = list()
   for i in histo:
      q = q + i
      cdf.append(q)

我正在遍历数组，但程序执行时间很长。有没有一个内置函数可以实现这个功能？

- omar

---

请参阅经验分布函数。 - djvg

Answer 1

使用直方图是一种解决方法，但它涉及对数据进行分组。 对于绘制经验数据的CDF来说，这并非必要。 假设F（x）表示小于x的条目计数，则当我们看到一个测量时，它会增加一个。 因此，如果我们对样本进行排序，然后在每个点上将计数增加1（或者将分数增加1/N），并将一个绘制到另一个上，我们将看到“精确”的（即未经分组的）经验CDF。

以下代码示例演示了该方法

import numpy as np
import matplotlib.pyplot as plt

N = 100
Z = np.random.normal(size = N)
# method 1
H,X1 = np.histogram( Z, bins = 10, normed = True )
dx = X1[1] - X1[0]
F1 = np.cumsum(H)*dx
#method 2
X2 = np.sort(Z)
F2 = np.array(range(N))/float(N)

plt.plot(X1[1:], F1)
plt.plot(X2, F2)
plt.show()

它输出以下内容

- Dan

---

根据numpy.histogram文档：__normed__等同于__density__参数，但对于不等宽度的箱子会产生错误的结果。从1.15.0版本开始更改：实际发出DeprecationWarnings。 - Oliver Prislan

你会如何处理在 Z 中的重复数值？ - djvg

Answer 2

我不太确定你的代码在做什么，但如果你有由numpy.histogram返回的hist和bin_edges数组，你可以使用numpy.cumsum来生成直方图内容的累积和。

>>> import numpy as np
>>> hist, bin_edges = np.histogram(np.random.randint(0,10,100), normed=True)
>>> bin_edges
array([ 0. ,  0.9,  1.8,  2.7,  3.6,  4.5,  5.4,  6.3,  7.2,  8.1,  9. ])
>>> hist
array([ 0.14444444,  0.11111111,  0.11111111,  0.1       ,  0.1       ,
        0.14444444,  0.14444444,  0.08888889,  0.03333333,  0.13333333])
>>> np.cumsum(hist)
array([ 0.14444444,  0.25555556,  0.36666667,  0.46666667,  0.56666667,
        0.71111111,  0.85555556,  0.94444444,  0.97777778,  1.11111111])

- user545424

---

然而，这引入了一个分箱步骤，对于累积分布来说是不必要的。 - hans_meine

“normed”这个关键词因为其混乱/错误的行为在Numpy 1.6中已经被弃用，将会在Numpy 2.0中被移除。如果bin不在[0,1]之间，则代码存在错误。请添加x=np.cumsum(hist); x=(x - x.min()) / x.ptp()。 - ArtificiallyIntelligence

@hans_meine 没错。这个问题有更好的解决方案吗？ - a06e

@becko 丹的回复中包含了基于直方图和“精确”解决方案（“方法2”）两种方法。 - hans_meine

啊，是的，我错过了那个。谢谢 @hans_meine - a06e

Answer 3

更新到numpy版本1.9.0。user545424的答案在1.9.0中无法正常工作。以下方法可行：

>>> import numpy as np
>>> arr = np.random.randint(0,10,100)
>>> hist, bin_edges = np.histogram(arr, density=True)
>>> hist = array([ 0.16666667,  0.15555556,  0.15555556,  0.05555556,  0.08888889,
    0.08888889,  0.07777778,  0.04444444,  0.18888889,  0.08888889])
>>> hist
array([ 0.1       ,  0.11111111,  0.11111111,  0.08888889,  0.08888889,
    0.15555556,  0.11111111,  0.13333333,  0.1       ,  0.11111111])
>>> bin_edges
array([ 0. ,  0.9,  1.8,  2.7,  3.6,  4.5,  5.4,  6.3,  7.2,  8.1,  9. ])
>>> np.diff(bin_edges)
array([ 0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9])
>>> np.diff(bin_edges)*hist
array([ 0.09,  0.1 ,  0.1 ,  0.08,  0.08,  0.14,  0.1 ,  0.12,  0.09,  0.1 ])
>>> cdf = np.cumsum(hist*np.diff(bin_edges))
>>> cdf
array([ 0.15,  0.29,  0.43,  0.48,  0.56,  0.64,  0.71,  0.75,  0.92,  1.  ])
>>>

- offwhitelotus

---

user12287，我感觉修改别人的答案有点奇怪。此外，不同版本可能会有不同的答案。 - offwhitelotus

Answer 4

为了补充Dan的解决方案。 如果样本中有多个相同的值，您可以使用numpy.unique：

Z = np.array([1,1,1,2,2,4,5,6,6,6,7,8,8])
X, F = np.unique(Z, return_index=True)
F=F/X.size

plt.plot(X, F)

- Alex

---

这将给你大于1的F值。也许你想使用F = F / float(F.max())（还要注意整数除法会给使用Python 2x的人带来问题）。 - ali_m

这个答案已经有些年头了，感谢你的评论和回答。我在每一个回答中都看到了我三年前的初级方法。 - omar

@Alex 这不太正确，因为对于出现多次的条目，它应该上升超过1/N。你是正确的，我的解决方案只对最后一个出现的情况正确，但它会正确绘制。 - Dan

原则上，您正在使用计数，但Python在F中使用从零开始的索引，因此也许您的意思是 (F + 1) / (F[-1] + 1)。 - Jthorpe

Answer 5

现有的答案要么使用柱状图，要么不能很好地处理重复值（无论是忽略重复值还是产生包含相同x值的多个y值的CDF）。我建议采用以下方法：

x, CDF_counts = np.unique(data, return_counts = True)
y = np.cumsum(CDF_counts)/np.sum(CDF_counts)

- Andy

Answer 6

我不确定是否有现成的答案，确切的做法是定义一个函数，例如：

def _cdf(x,data):
    return(sum(x>data))

这将非常快。

- user1505725