Гистограмма Numpy на многомерном массиве

учитывая np.array формы (n_days, n_lat, n_lon), я хотел бы вычислить гистограмму с фиксированными ячейками для каждой ячейки lat-lon (т.е. распределение ежедневных значений).

Простое решение проблемы — пройтись по ячейкам и вызвать 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, а затем использовал причудливую индексацию для обновления 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

но это не работает, потому что оператор добавления на месте (+=), похоже, не поддерживает несколько обновлений одной и той же ячейки.

спасибо, Питер


person Peter Prettenhofer    schedule 17.09.2013    source источник
comment
похоже, что github.com/numpy/numpy/pull/2821 обратился к причудливой индексации и проблема на месте. Причина, по которой numpy не разрешает несколько обновлений, заключается в том, что a[idx] += 1 не будет таким же, как a[idx] = a[idx] + 1 .   -  person Peter Prettenhofer    schedule 17.09.2013
comment
Используйте np.histogram2d с аргументом ключевого слова weights.   -  person Jaime    schedule 17.09.2013
comment
@ Джейме, как бы я использовал weights? Я не хочу делать двумерную гистограмму.   -  person Peter Prettenhofer    schedule 19.09.2013
comment
Также есть функция np.histogramdd.   -  person Jaime    schedule 19.09.2013


Ответы (2)


Вы можете использовать numpy.apply_along_axis, чтобы устранить цикл.

hist, bin_edges = apply_along_axis(lambda x: histogram(x, bins=bins), 0, B)
person Greg Whittier    schedule 18.09.2013
comment
@PeterPrettenhofer только что исправил опечатку. у лямбды буквы были переставлены. Надеюсь, что это работает для вас. - person Greg Whittier; 19.09.2013

Может это работает?:

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))
person Raphael Roth    schedule 17.09.2013
comment
это дает мне в основном то же, что и bins.searchsorted(B), массив формы (n_days, n_lat, n_lon), но сложная часть заключается в том, как преобразовать его в (n_bins, n_lat, n_lon). - person Peter Prettenhofer; 17.09.2013