だから僕はpandasを辞めた【NumPyだけでgroupby.mean()する３つの方法 篇】

# だから僕はpandasを辞めた【NumPyだけでgroupby.mean()する３つの方法 篇】

## な、なんだこのふざけた記事は

低スペPCにpandasはもうこりごりだよ～（トホホ）という人はデータをNumPyでごちゃごちゃ扱っていることであろう。しかしNumPyには`pd.DataFrame.groupby`的な関数がないので、泣く泣くpandasでごちゃごちゃしている人もいるかもしれない。しかししかし馬鹿な、pandasでできることがNumPyでできないはずがないはずだ、という人は頑張ってNumPyのドキュメントをひっくり返して、いろんな関数を眺めて考えてみてほしい（ただしforループはNG）。そうです、できますできるんです。

いや考えたってわからないし、という人のために、幾星霜過去にある賢者が『[Numpyのみでグループ別に集計し平均を出す方法（for文、Pandasなし）](https://blog.shikoan.com/numpy-grouping-no-pandas/)』という記事を書いてくれた（神）。しかし、科学技術の発展した現代、よりよい方法がありそうな気がしたので、Numpyの関数の勉強がてらいろいろ考えてみた、その戦いの記録をつたない日本語で綴ったメモをpandasの引退者に捧げよう、というのがこの記事である。

## やること

今回は、グループごとに平均値と標準偏差を求める、つまり`pd.DataFrame.groupby.mean()`・`pd.DataFrame.groupby.std()`をNumPyだけで計算したい。

## 用いるデータ

おなじみのアヤメのデータ。行はシャッフルしております。

```py
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)`で可能。

```py
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)
```

<details><summary>余談</summary><div>

なお、このとき、それなりにデータ数がある場合は`np.unique()`よりも`pd.factorize()`の方が高速であることが知られている。

```py
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)
```

> ![pd_to_np_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/394333/9bae2c06-2539-d30a-d467-14f9660e7fdb.png)


どうやら1000行を超える`Series`（のみならず`ndarray`であっても）を変換する場合は`pd.factorize()`を使用すべきなのかもしれない。</div></details>

NumPy配列に変換（ここでは説明の簡便のためレコード配列を使う）。pandasと決別する瞬間。

```py
arr = df.to_records(index=False)
```

## groupbyする方法

### その１（onehot）

まず、上掲の記事の手法を紹介。

> アヤメの種類をOne-hotベクトル化します。これと集計したい花弁の長さのベクトルをかけ合わせて、列単位の平均を取ります。
> one-hotベクトルをそのまま平均を取ってしまうと0の影響を受けてしまうので事前にNaNに置き換えます。

```py:その１
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()`は遅いということが知られている。平均の式は、

```math
mean = \frac{1}{n} \sum_{i=1}^{N} x_i
```

であるから、`np.nansum()`を使った式に置き換えると以下のようになる。

```py:
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`を得ることができる。

```py
%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]])
```

よって、

```py:その１改
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()`と同様の結果を得ることができた。

```text:その１の考え方
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
```

しかし、元の配列において、要素がグループごとにまとまって並んでいれば、以下のような考え方ができる。

```text:その２の考え方
g |  A  A  A  B  B  B  C  C  C
    ----A--- ----B--- ----C---
                 v
     mean_A   mean_B   mean_C
```

そこで、まずグループごとに並べるために、Name列に基づいてソートを行う。

```py
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列の各隣接する値を見ていけばいい。例えば、

```py
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の位置を取得すれば、それが各グループの境界の位置になる。

```py
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()`を使えば一発で求められる。

```py
counts = np.ediff1d(cut_index)
# array([3, 3, 3, 3], dtype=int64)
```

続いて、グループ内の要素の合計を求める。ところでところで、`np.add(x1, x2)`という、２つの数字や配列を足すことができる**Numpyユニバーサル関数**がある。

```py
np.add(1, 2)
# 3
np.add([0, 1], [0, 2])
# array([0, 3])
```

基本的に`np.add()`はただの`+`でしかないので、この関数自体を使うことはまずない。重要なのはユニバーサル関数のメソッドである。例えば`np.ufunc.reduce(array)`というのがあって、これは`array`から要素を取り出していって関数を適用させることができる。どういうことかいうと例えば、

```py
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]])`の配列が得られる。

```py
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`はまさに先程求めた「グループの境界」である。つまりまとめると、

```py:その２
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()`を使うと、これを一気に数えられる。どういうことかというと例えば、

```py
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`オプションがあり、重み付けができる。つまり例えば、

```py
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.])
```

この重み付けを利用すれば、要素の合計を求めることができる。したがって、

```py:その３
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]
```

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

```py: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')
```

![pd_to_np_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/394333/f3a5d73f-8822-2ddb-30ed-a0bffb3ddff8.png)


その３の勝利（グループの数にもよるだろうけど）。

間違ってたらコメントお願いします。

---

### おまけ：今回登場した関数たち

覚えておこう！（クリックすると公式ドキュメントに飛びます）

[`np.nansum()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.nansum.html#numpy.nansum) - 非数（NaN）をゼロとして扱って、指定された軸上の配列要素の合計を返します。
[`np.nanmean()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.nanmean.html#numpy.nanmean) - NaNを無視して、指定された軸に沿って算術平均を計算します。
[`np.nanstd()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.nanstd.html#numpy.nanstd) - NaNを無視して、指定された軸に沿って標準偏差を計算します。

[`np.add()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.add.html#numpy.add) - 引数を要素ごとに足し算します。
[`np.ufunc.reduce()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.ufunc.reduce.html#numpy.ufunc.reduce) - １つの軸に沿ってufuncを適用することにより、配列の次元を1つ減らします。
[`np.ufunc.reduceat()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.ufunc.reduceat.html#numpy.ufunc.reduceat) - 単一の軸で指定されたスライスを使用して（ローカル）リデュースを実行します。

[`np.argsort()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.argsort.html#numpy.argsort) - 配列をソートするインデックスを返します。
[`np.nonzero()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.nonzero.html#numpy.nonzero) - ゼロ以外の要素のインデックスを返します。

[`np.bincount()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.bincount.html#numpy.bincount) - 自然数のみから構成される配列において、各値の出現回数をカウントします。
[`np.unique()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.unique.html#numpy.unique) - 配列の一意の要素を見つけます。

[`np.ediff1d()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.ediff1d.html#numpy.ediff1d) - 配列の連続した要素間の差。

[`np.concatenate()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.concatenate.html#numpy.concatenate) - 既存の軸に沿って配列のシーケンスを結合します。