だから僕は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)

余談

なお、このとき、それなりにデータ数がある場合は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)

どうやら1000行を超える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]

まとめ（タイトルに釣られてきた人はここだけ読めばいい）

numpy_groupby.py

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() - 既存の軸に沿って配列のシーケンスを結合します。