numpyのhistogram2dでbin内の中央値を利用する方法

Question

前提・実現したいこと

(x, y)座標と特定の値valueのカラムで構成されたデータが手元にあり、これをもとに二次元のカラーメッシュを作ろうとしています。具体的には、xy座標をもとにグリッドを作成し、それぞれのグリッドに含まれるデータポイントのvalueの合計値、平均値、中央値に基づいたカラーメッシュが作りたいです。
合計値、平均値は以下の方法で計算できました。

import numpy as np
import pandas as pd

step = 5
xx = np.arange(0, 200+step, step)
yy = np.arange(0, 200+step, step)

df = pd.read_csv(data_path)

# 合計値
sum_hist, _x, _y = np.histogram2d(df['x'], df['y'], bins=(xx, yy), weights=df['value'])
# 平均値
count, _x, _y = np.histogram2d(df['x'], df['y'], bins=(xx, yy))
mean_hist = np.divide(sum_hist, count)

sum_hist、mean_histにはそれぞれ、xxとyyに基づいて分割されたデータポイントの持つvalueの合計値、平均値がそれぞれのグリッドに格納された二次元配列が生成されています。histogram2dはweightを入力するとグリッド内のデータの持つ合計値、入力しないとグリッド内の単純なデータ数が計算されるため、np.divideを使ってやればグリッド内のvalueの平均値が計算できます。

発生している問題

これと似たようなnpの関数を使って（つまり、Pythonネイティブのfor-loopを利用せず）各グリッド内のvalueの中央値を計算する方法はあるでしょうか？中央値を計算するにはどうしてもグリッド内のデータを整列する必要があるため、histogram2dをそのまま使う方法は思いつきませんでした。データが巨大なため、for-loopを使うのはとても遅くなってしまいます。numpyに詳しい方、あるいはアルゴリズムに詳しい方、ご教授お願いします。

環境

バグではないので関係はないと思いますが、環境はpython3.7, numpy1.15です。

Answer 1

scipy.stats.binned_statistic_2d をお使いください。
この関数を使用すると、ビンごとの平均、中央値、和の計算が行えます。

コード

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binned_statistic_2d


np.random.seed(1)

x, y, z = np.random.uniform(0, 10, (3, 1000))

# ビン
xbins = np.linspace(0, 10, 11)
ybins = np.linspace(0, 10, 11)

# 合計値
sum_, xedges, yedges, binnumber = binned_statistic_2d(
    x, y, z, statistic="sum", bins=(xbins, ybins)
)

# 中央地
median, xedges, yedges, binnumber = binned_statistic_2d(
    x, y, z, statistic="median", bins=(xbins, ybins)
)

# 平均値
mean, xedges, yedges, binnumber = binned_statistic_2d(
    x, y, z, statistic="mean", bins=(xbins, ybins)
)

fig, [ax1, ax2, ax3] = plt.subplots(1, 3, figsize=(10, 3))
ax1.pcolormesh(sum_)
ax1.set_title("sum")

ax2.pcolormesh(median)
ax2.set_title("median")

ax3.pcolormesh(mean)
ax3.set_title("mean")

plt.show()