有没有更好的方法来实现直方图？答案

【问题标题】：Is there a better way of implementing a histogram?有没有更好的方法来实现直方图？
【发布时间】：2020-05-04 12:33:58
【问题描述】：

我有一个 2D uint16 numpy 数组，我想计算这个数组的直方图。我使用的功能是：

def calc_hist(source):
    hist = np.zeros(2**16, dtype='uint16')
    for i in range(source.shape[0]):
       for j in range(source.shape[1]):
           hist[source[i, j]] = hist[source[i, j]] + 1
    return hist

此函数执行时间过长。据我了解，numpy 模块中有一个直方图函数，但我不知道如何使用它。我试过了：

hist,_ = np.histogram(source.flatten(), bins=range(2**16))

但我得到的结果与我自己的函数不同。如何调用 numpy.histogram 来获得相同的结果？还是有其他选择？

【问题讨论】：

您能否添加数据样本！
您能否修复您发布的代码以使其实际工作？
它有效，你可以随机化一个 uint16 二维数组..
不，它没有。剪下问题中的代码，将其粘贴到文件或 IDE 中并尝试运行它。即使您修复了zeroes，代码中仍然存在不平衡的括号
您的代码使用准确的数字，因此它始终是范围 (0, 2^16-1)。 numpy 直方图默认为范围 (a.min(), a.max())，然后拆分为 num 个 bin。您可能需要明确地给它（浮点）范围（0, 2^16）以获得相似的结果......

标签： python performance numpy histogram

【解决方案1】：

对于数据类型为uint16 的输入，numpy.bincount 应该可以正常工作：

hist = np.bincount(source.ravel(), minlength=2**16)

您的函数所做的几乎与 bincount 完全相同，但 bincount 是用 C 实现的。

例如，以下检查bincount 的使用是否与calc_hist 函数的结果相同：

In [159]: rng = np.random.default_rng()

In [160]: x = rng.integers(0, 2**16, size=(1000, 1000))                                      

In [161]: h1 = calc_hist(x)

In [162]: h2 = np.bincount(x.ravel(), minlength=2**16)

In [163]: (h1 == h2).all()  # Verify that they are the same.
Out[163]: True

使用 ipython %timeit 命令检查性能。您可以看到使用bincount 的速度要快得多。

In [164]: %timeit calc_hist(x)
2.66 s ± 21.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [165]: %timeit np.bincount(x.ravel(), minlength=2**16)
3.13 ms ± 100 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

【讨论】：

【解决方案2】：

正如 Corley Brigman 指出的那样，传递 bins=range(x) 确定了 bin 边缘 [1]。因此，您最终会得到 x-1 个带有对应边 [0, 1)、[1, 2)、...、[x-1, x] 的 bin。

在您的情况下，您将有 2^16-1 个垃圾箱。要修复它，只需使用range(2**16+1)。

[1]https://numpy.org/doc/stable/reference/generated/numpy.histogram.html?highlight=histogram#numpy.histogram

【讨论】：

只有在最后一个 bin 是唯一不正确的 bin 时，这才是正确的答案