だから僕はpandasを辞めた【NumPyだけでgroupby.mean()する３つの方法 篇】
な、なんだこのふざけた記事は
低スペPCにpandasはもうこりごりだよ～（トホホ）という人はデータをNumPyでごちゃごちゃ扱っていることであろう。しかしNumPyにはpd.DataFrame.groupby的な関数がないので、泣く泣くpandasでごちゃごちゃしている人もいるかもしれない。しかししかし馬鹿な、pandasでできることがNumPyでできないはずがないはずだ、という人は頑張ってNumPyのドキュメントをひっくり返して、いろんな関数を眺めて考えてみてほしい（ただしforループはNG）。そうです、できますできるんです。
いや考えたってわからないし、という人のために、幾星霜過去にある賢者が『Numpyのみでグループ別に集計し平均を出す方法（for文、Pandasなし）』という記事を書いてくれた（神）。しかし、科学技術の発展した現代、よりよい方法がありそうな気がしたので、Numpyの関数の勉強がてらいろいろ考えてみた、その戦いの記録をつたない日本語で綴ったメモをpandasの引退者に捧げよう、というのがこの記事である。
やること
今回は、グループごとに平均値と標準偏差を求める、つまりpd.DataFrame.groupby.mean()・pd.DataFrame.groupby.std()をNumPyだけで計算したい。
用いるデータ
おなじみのアヤメのデータ。行はシャッフルしております。
import numpy as np
import pandas as pd
df = (pd.read_csv('https://raw.githubusercontent.com/pandas-dev/'
'pandas/master/pandas/tests/io/data/csv/iris.csv')
.sample(frac=1, random_state=1))
df.head()
| SepalLength | SepalWidth | PetalLength | PetalWidth | Name | |
|---|---|---|---|---|---|
| 14 | 5.8 | 4 | 1.2 | 0.2 | Iris-setosa |
| 98 | 5.1 | 2.5 | 3 | 1.1 | Iris-versicolor |
| 75 | 6.6 | 3 | 4.4 | 1.4 | Iris-versicolor |
| 16 | 5.4 | 3.9 | 1.3 | 0.4 | Iris-setosa |
| 131 | 7.9 | 3.8 | 6.4 | 2 | Iris-virginica |
Numpyのオブジェクトに変換する前に注意点。NumPyで効率よくデータを扱う場合には、数値型以外のデータ（文字列・オブジェクト型など）は数値型に変換するのが必須である。文字列操作をしたい場合は標準のリストを使うのが良い。
今回の「Name」列はカテゴリーを表しているだけで、データ上文字列に意味はないので、いわゆるラベルエンコードを行って数字に変換する。これはnp.unique(array, return_inverse=True)で可能。
unique_names, df['Name'] = np.unique(df['Name'].to_numpy(), return_inverse=True)
# (array(['Iris-setosa', 'Iris-versicolor', 'Iris-virginica'], dtype=object),
# array([0, 1, 1, 0, 2, 1, 2, 0, 0, 2, 1, 0, 2, 1, 1, 0, 1, 1, 0, 0, 1, 1,
# ...,
# 1, 0, 1, 1, 1, 1, 2, 0, 0, 2, 1, 2, 1, 2, 2, 1, 2, 0], dtype=int64))
# これでも可
df['Name'], unique_names = pd.factorize(df['Name'], sort=True)
なお、このとき、それなりにデータ数がある場合は どうやら1000行を超える余談
np.unique()よりもpd.factorize()の方が高速であることが知られている。import benchit
def Numpy_unique(df):
return np.unique(df.iloc[:, 0].to_numpy(), return_inverse=True)
def Pandas_factorize(df):
return pd.factorize(df.iloc[:, 0], sort=True)
funcs = [Numpy_unique, Pandas_factorize]
inputs = [pd.DataFrame(list('abcde')).sample(n=n, replace=True, random_state=1)
for n in 2**np.arange(21)]
t = benchit.timings(funcs, inputs)
t.plot(logy=True, logx=True)
Series（のみならずndarrayであっても）を変換する場合はpd.factorize()を使用すべきなのかもしれない。
NumPy配列に変換（ここでは説明の簡便のためレコード配列を使う）。pandasと決別する瞬間。
arr = df.to_records(index=False)
groupbyする方法
その１（onehot）
まず、上掲の記事の手法を紹介。
アヤメの種類をOne-hotベクトル化します。これと集計したい花弁の長さのベクトルをかけ合わせて、列単位の平均を取ります。
one-hotベクトルをそのまま平均を取ってしまうと0の影響を受けてしまうので事前にNaNに置き換えます。
onehot = np.eye(arr['Name'].max()+1)[arr['Name']]
onehot[onehot == 0] = np.nan
onehot_value = onehot * arr['PetalWidth'][:, None]
result_mean = np.nanmean(onehot_value, 0)
# [0.244 1.326 2.026]
result_std = np.nanstd(onehot_value, 0)
# [0.10613199 0.19576517 0.27188968]
この手法の問題点としてnp.nanmean()やnp.nanstd()は遅いということが知られている。平均の式は、
mean = \frac{1}{n} \sum_{i=1}^{N} x_i
であるから、np.nansum()を使った式に置き換えると以下のようになる。
result_mean = np.nansum(onehot_value, 0)/(~np.isnan(onehot_value)).sum(0)
result_std = np.sqrt(np.nansum(onehot_value**2, 0)
/ (~np.isnan(onehot_value)).sum(0) - result_mean**2)
さらにonehot_value中のnp.nanは0を置き換えたものであるが、そもそも0を何個足しても合計はかわらないのだから、np.nanに置換せずかつそのままnp.sum()を使うことができる。
onehotにnp.nanを代入する必要がないということは、floatデータ型にする必要がないので、以下の方法でより高速にonehotを得ることができる。
%timeit np.eye(arr['Name'].max()+1)[arr['Name']]
# 21.2 µs ± 919 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
# array([[1., 0., 0.],
# [0., 1., 0.],
# [0., 1., 0.],
# ...,
# [1., 0., 0.]])
%timeit np.arange(arr['Name'].max()+1) == arr['Name'][:, None]
# 20.1 µs ± 205 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
# array([[ True, False, False],
# [False, True, False],
# [False, True, False],
# ...,
# [ True, False, False]])
よって、
onehot = np.arange(arr['Name'].max()+1) == arr['Name'][:, None]
onehot_value = onehot * arr['PetalWidth'][:, None]
sums = onehot_value.sum(0)
counts = onehot.sum(0)
result_mean = sums/counts
result_std = np.sqrt((onehot_value**2).sum(0)/counts - result_mean**2)
その２（reduceat）
groupby のポイントは、どうやってグループごとに計算するかということにある。
その１では、値をグループごとの異なる列に移動して平均することでgroupby.mean()と同様の結果を得ることができた。
g | A B B C A B A C C
v
A | 1 0 0 0 1 0 1 0 0 | mean_A
B | 0 1 1 0 0 1 0 0 0 | -> mean_B
C | 0 0 0 1 0 0 0 1 1 | mean_C
しかし、元の配列において、要素がグループごとにまとまって並んでいれば、以下のような考え方ができる。
g | A A A B B B C C C
----A--- ----B--- ----C---
v
mean_A mean_B mean_C
そこで、まずグループごとに並べるために、Name列に基づいてソートを行う。
sorter_index = np.argsort(arr['Name'])
name_s = arr['Name'][sorter_index]
# array([0, 0, 0, ..., 1, 1, 1, ..., 2, 2, 2], dtype=int64)
pwidth_s = arr['PetalWidth'][sorter_index]
続いて、グループの境界を得る。これはソート済Name列の各隣接する値を見ていけばいい。例えば、
ex_array = np.array([0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3])
# 0------ 1------ 2------ 3------ ←この境界の位置がほしい
ex_array[1:] != ex_array[1:]
"""考え方
0 [0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3] # ex_array[1:]
[0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3] 3 # ex_array[:-1]
[F, F, T, F, F, T, F, F, T, F, F] # ex_array[1:] != ex_array[1:]
（上２つが一致すればF, 一致しなければT）
[0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3] # ex_array
[T, F, F, T, F, F, T, F, F, T, F, F, T] # 本当にほしかった配列
0------ 1------ 2------ 3------
# ex_array[1:] != ex_array[1:] の前後にそれぞれ1個ずつTrueを加えれば、
「T～次のTの前」がグループになる。
"""
# 前後を補完したバージョン
np.concatenate(([True], ex_array[1:] != ex_array[:-1], [True]))
# array([ True, False, False, False, False, False, False, False, False, False, False, False, True])
np.ndarray.nonzero()で、この配列のTrueの位置を取得すれば、それが各グループの境界の位置になる。
cut_index = np.concatenate(([True], ex_array[1:] != ex_array[:-1], [True])).nonzero()[0]
# array([ 0, 3, 6, 9, 12], dtype=int64)
"""考え方
0 1 2 3 4 5 6 7 8 9 10 11 12 # ex_arrayのindex
[0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3] # ex_array
[0, 3, 6, 9, 12] # cut_index: 各グループの境界（始点）
"""
繰り返しになるが、平均は「合計÷個数」なのでそれぞれを求める。グループ内の要素の個数は、cut_indexの隣接する値の差分だから、np.diff()やnp.ediff1d()を使えば一発で求められる。
counts = np.ediff1d(cut_index)
# array([3, 3, 3, 3], dtype=int64)
続いて、グループ内の要素の合計を求める。ところでところで、np.add(x1, x2)という、２つの数字や配列を足すことができるNumpyユニバーサル関数がある。
np.add(1, 2)
# 3
np.add([0, 1], [0, 2])
# array([0, 3])
基本的にnp.add()はただの+でしかないので、この関数自体を使うことはまずない。重要なのはユニバーサル関数のメソッドである。例えばnp.ufunc.reduce(array)というのがあって、これはarrayから要素を取り出していって関数を適用させることができる。どういうことかいうと例えば、
ex_array = np.arange(10)
result = ex_array[0]
for value in ex_array[1:]:
result = np.add(result, value)
print(result)
# 45
np.add.reduce(ex_array)
# 45
さらに、np.ufunc.reduceat(array, indices)を使うとnp.ufunc.reduce(a[indices[i]:indices[i+1]])の配列が得られる。
ex_array = np.arange(10)
indices = [0, 1, 6, len(ex_array)]
results = []
for i in range(len(indices) - 1):
result = np.add.reduce(ex_array[indices[i]:indices[i+1]])
results.append(result)
print(results)
# [0, 15, 30]
indices = [0, 1, 6]
np.add.reduceat(ex_array, indices)
# array([0, 15, 30], dtype=int32)
ここにおけるindicesはまさに先程求めた「グループの境界」である。つまりまとめると、
sorter_index = np.argsort(arr['Name'])
name_s = arr['Name'][sorter_index]
pwidth_s = arr['PetalWidth'][sorter_index]
cut_index, = np.concatenate(([True], name_s[1:] != name_s[:-1], [True])).nonzero()
sums = np.add.reduceat(pwidth_s, cut_index[:-1])
counts = np.ediff1d(cut_index)
result_mean = sums/counts
# [0.244 1.326 2.026]
result_std = np.array([np.std(pwidth_s[s:e], ddof=0)
for s, e in zip(cut_index[:-1], cut_index[1:])])
# [0.10613199 0.19576517 0.27188968]
この手法の悲しいところは、標準偏差を求める際はうまくnp.ufunc.reduceat()が使えないために、forループに頼ることになることだ（その代わり直接np.std()を使うことができる）。
その３（bincount）
その１とかその２が馬鹿らしくなるくらい、最も簡単な方法。
Name列はラベルエンコードにより得た配列のため、{0, 1, 2}が並んでいる。np.bincount()を使うと、これを一気に数えられる。どういうことかというと例えば、
ex_array = np.array([3, 4, 3, 3, 2, 2, 4, 4, 1, 3])
results = []
for value in range(ex_array.max()+1):
result = sum(ex_array == value)
results.append(result)
print(results)
# [0, 1, 2, 4, 3]
# ex_arrayには、0が0個、1が1個、2が2個、3が4個、4が3個
np.bincount(ex_array)
# array([0, 1, 2, 4, 3], dtype=int64)
np.bincount()にはweightsオプションがあり、重み付けができる。つまり例えば、
ex_array = np.array([1, 2, 2, 3, 3, 3, 3, 4, 4, 4])
weights = np.array([2, 1, 1, 1, 1, 1, 1, 2, 3, 4])
results = []
for value in range(ex_array.max()+1):
idx = np.nonzero(ex_array == value)[0]
result = sum(weights[idx])
results.append(result)
print(results)
# [0, 2, 2, 4, 9]
np.bincount(ex_array, weights)
# array([0., 2., 2., 4., 9.])
この重み付けを利用すれば、要素の合計を求めることができる。したがって、
counts = np.bincount(arr['Name'])
sums = np.bincount(arr['Name'], arr['PetalWidth'])
result_mean = sums/counts
# [0.244 1.326 2.026]
deviation_array = result_mean[arr['Name']] - arr['PetalWidth']
result_std = np.sqrt(np.bincount(arr['Name'], deviation_array**2) / counts)
# [0.10613199 0.19576517 0.27188968]
まとめ（タイトルに釣られてきた人はここだけ読めばいい）
import benchit
import numpy as np
import pandas as pd
def type1(value, group):
onehot = np.eye(group.max()+1)[group]
onehot[onehot == 0] = np.nan
onehot_value = onehot * value[:, None]
result_mean = np.nanmean(onehot_value, 0)
result_std = np.nanstd(onehot_value, 0)
return result_mean, result_std
def type1rev(value, group):
onehot = np.arange(group.max()+1) == group[:, None]
onehot_value = onehot * value[:, None]
sums = onehot_value.sum(0)
counts = onehot.sum(0)
result_mean = sums/counts
result_std = np.sqrt((onehot_value**2).sum(0)/counts - result_mean**2)
return result_mean, result_std
def type2(value, group):
sorter_index = np.argsort(group)
group_s = group[sorter_index]
value_s = value[sorter_index]
cut_index, = np.concatenate(([True], group_s[1:] != group_s[:-1], [True])).nonzero()
sums = np.add.reduceat(value_s, cut_index[:-1])
counts = np.ediff1d(cut_index)
result_mean = sums/counts
result_std = np.array([np.std(value_s[s:e], ddof=0)
for s, e in zip(cut_index[:-1], cut_index[1:])])
return result_mean, result_std
def type3(value, group):
counts = np.bincount(group)
sums = np.bincount(group, value)
result_mean = sums/counts
result_std = np.sqrt(
np.bincount(group, (result_mean[group] - value)**2) / counts)
return result_mean, result_std
df = (pd.read_csv('https://raw.githubusercontent.com/pandas-dev/'
'pandas/master/pandas/tests/io/data/csv/iris.csv')
.sample(frac=1, random_state=1))
arr_value = df['PetalWidth'].to_numpy()
arr_group = np.unique(df['Name'].to_numpy(), return_inverse=True)[1]
# 比較
funcs = [type1, type1rev, type2, type3]
inputs = {n: [arr_value[np.random.randint(0, 150, n)],
arr_group[np.random.randint(0, 150, n)]]
for n in 2 ** np.arange(1, 21)}
t = benchit.timings(funcs, inputs, multivar=True, input_name='Array-length')
t.plot(logy=True, logx=True).savefig('a.png')
その３の勝利（グループの数にもよるだろうけど）。
間違ってたらコメントお願いします。
おまけ：今回登場した関数たち
覚えておこう！（クリックすると公式ドキュメントに飛びます）
np.nansum() - 非数（NaN）をゼロとして扱って、指定された軸上の配列要素の合計を返します。
np.nanmean() - NaNを無視して、指定された軸に沿って算術平均を計算します。
np.nanstd() - NaNを無視して、指定された軸に沿って標準偏差を計算します。
np.add() - 引数を要素ごとに足し算します。
np.ufunc.reduce() - １つの軸に沿ってufuncを適用することにより、配列の次元を1つ減らします。
np.ufunc.reduceat() - 単一の軸で指定されたスライスを使用して（ローカル）リデュースを実行します。
np.argsort() - 配列をソートするインデックスを返します。
np.nonzero() - ゼロ以外の要素のインデックスを返します。
np.bincount() - 自然数のみから構成される配列において、各値の出現回数をカウントします。
np.unique() - 配列の一意の要素を見つけます。
np.ediff1d() - 配列の連続した要素間の差。
np.concatenate() - 既存の軸に沿って配列のシーケンスを結合します。