多维数组上的Numpy直方图

Question

给定一个形状为(n_days, n_lat, n_lon)的 np.array 数组，我想为每个经纬度单元格计算具有固定区间的直方图（即每日值的分布）。

问题的简单解决方法是循环遍历单元格并对每个单元格调用np.histogram函数：

bins = np.linspace(0, 1.0, 10)
B = np.rand(n_days, n_lat, n_lon)
H = np.zeros((n_bins, n_lat, n_lon), dtype=np.int32)
for lat in range(n_lat):
    for lon in range(n_lon):
        H[:, lat, lon] = np.histogram(A[:, lat, lon], bins=bins)[0]
# note: code not tested

但这种方法速度较慢。有没有更有效率的解决方案，不需要使用循环？

我研究了np.searchsorted以获取B中每个值的bin索引，然后使用fancy indexing来更新H：

bin_indices = bins.searchsorted(B)
H[bin_indices.ravel(), idx[0], idx[1]] += 1  # where idx is a index grid given by np.indices
# note: code not tested

但是这样做行不通，因为就地加法运算符（+=）似乎不支持同一单元格的多次更新。

谢谢， 彼得

- Peter Prettenhofer

---

似乎 https://github.com/numpy/numpy/pull/2821 解决了花式索引和原地操作的问题。Numpy 不允许多次更新的原因是 "a[idx] += 1" 不等同于 "a[idx] = a[idx] + 1"。 - Peter Prettenhofer

使用np.histogram2d函数并带上weights关键字参数。 - Jaime

@Jaime，我该如何使用“weights”？我不想做一个二维直方图。 - Peter Prettenhofer

还有一个 np.histogramdd 函数。 - Jaime

Answer 1

你可以使用 numpy.apply_along_axis() 来消除循环。

import numpy as np

hist, bin_edges = np.apply_along_axis(lambda x: np.histogram(x, bins=bins), 0, B)

- Greg Whittier

---

@PeterPrettenhofer 刚刚修正了一个打字错误。lambda函数的字母位置错了。希望这样能对你有用。 - Greg Whittier

Answer 2

也许这个可以工作：

import numpy as np
n_days=31
n_lat=10
n_lon=10
n_bins=10
bins = np.linspace(0, 1.0, n_bins)
B = np.random.rand(n_days, n_lat, n_lon)


# flatten to 1D
C=np.reshape(B,n_days*n_lat*n_lon)
# use digitize to get the index of the bin to which the numbers belong
D=np.digitize(C,bins)-1
# reshape the results back to the original shape
result=np.reshape(D,(n_days, n_lat, n_lon))

- Raphael Roth

---

这基本上给了我与“bins.searchsorted(B)”相同形状为“(n_days, n_lat, n_lon)”的数组，但棘手的部分是如何将其转换为“(n_bins, n_lat, n_lon)”。 - Peter Prettenhofer