NaN交じりの配列処理はnanxxx()関数よりもmasked arrayがおすすめ

# はじめに

NaN交じりのnumpy配列に関する処理として，`np.nanmean()`や`np.nansum()`などの関数がしばしば記事で紹介されるが，より複雑な処理も考慮すれば，`np.ma.masked_array`を用いる方が応用幅が広い．本記事では，平均の計算を題材として説明を行う．

# 通常の配列（NaNなし）

3×3の二次元行列を扱う．Aは平均を求める配列，Wは重みづけ平均の計算に用いる重み行列である．

```py:
A = np.random.randint(1, 10, (3, 3))
W = np.arange(9).reshape(3, 3)
print(A, W, sep='\n')
'''
[[3 8 2]
 [1 9 8]
 [2 3 9]]
[[0 1 2]
 [3 4 5]
 [6 7 8]]
'''
```

## np.mean()

numpy配列の平均をとる代表的なメソッドは`mean()`である．関数および配列自体のメソッドとして使用が可能である．

```py:
print(np.mean(A))
print(np.mean(A, axis=0))
print(np.mean(A, axis=1))
'''
5.0
[2.         6.66666667 6.33333333]
[4.33333333 6.         4.66666667]
'''
# ↑と全く同様の計算
print(A.mean())
print(A.mean(axis=0))
print(A.mean(axis=1))
```

## np.average()

`mean()`としばしば対比されるのが`average()`．`mean()`と同様，軸を指定した平均を算出することができる．

```py:
print(np.average(A))
print(np.average(A, axis=0))
print(np.average(A, axis=1))
```

それに加え，`average()`では，重みづけ平均の計算が可能．
重みづけ平均とは，各元と重みとの積の総和を重みの総和で割った値である．

```py:
print(np.sum(A*W)/np.sum(W))
print(np.sum(A*W, axis=0)/np.sum(W, axis=0))
print(np.sum(A*W, axis=1)/np.sum(W, axis=1))
'''
5.444444444444445
[1.66666667 5.41666667 7.73333333]
[4.         6.58333333 5.        ]
'''
# ↑と全く同様の計算
print(np.average(A, weights=W))
print(np.average(A, axis=0, weights=W))
print(np.average(A, axis=1, weights=W))
```

# NaN交じりの配列

上記の行列Aの対角成分を`np.nan`に置換する．
intで定義された配列にはfloatの`np.nan`を定義できないため，`astype()`によって型変換を行う．

```py:
NA = A.copy().astype(float)
NA[np.eye(3, dtype=bool)] = np.nan
print(NA)
'''
array([[nan,  8.,  2.],
       [ 1., nan,  8.],
       [ 2.,  3., nan]])
'''
```

`mean()`による演算はnp.nanを返す．

```py:
print(np.mean(NA)) # np.nan
print(NA.mean()) # np.nan
```

## np.nanmean()

NaN交じりの配列の計算には，`nanmean()`が使用できる．

```py:
print(np.nanmean(NA))
print(np.nanmean(NA, axis=0))
print(np.nanmean(NA, axis=1))
'''
4.0
[1.5 5.5 5. ]
[5.  4.5 2.5]
'''
```

しかし，NaN交じりの配列に対する重みづけ平均に対応する関数は存在しない．

# masked array

masked arrayは下記のように定義され，通常の`data`配列に`mask`配列を重ねた形となっている．

```py:
MA = np.ma.masked_array(A, mask=np.eye(3, dtype=bool))
print(MA, MA.data, MA.mask, sep='\n')
'''
[[-- 8 2]
 [1 -- 8]
 [2 3 --]]
[[3 8 2]
 [1 9 8]
 [2 3 9]]
[[ True False False]
 [False  True False]
 [False False  True]]
'''
```

## np.ma.mean()

`np.mean()`に該当する関数は`np.ma.mean()`である．メソッドも存在する．
実はmasked arrayに対して`np.mean()`を使用しても結果は全く変わらない．

```py:
print(np.ma.mean(MA))
print(np.ma.mean(MA, axis=0))
print(np.ma.mean(MA, axis=1))
'''
4.0
[1.5 5.5 5.0]
[5.0 4.5 2.5]
'''
# ↑と全く同様の計算
print(MA.mean())
print(MA.mean(axis=0))
print(MA.mean(axis=1))
```

## np.ma.average()

普通に平均をとる場合には，`np.average()`を用いることも可能．
しかしながら，重みづけ平均の計算には，`np.ma.average()`を用いなければならない．

```py:
print(np.ma.average(MA, weights=W)) # 88/24
print(np.ma.average(MA, axis=0, weights=W))
print(np.ma.average(MA, axis=1, weights=W))
'''
3.6666666666666665
[1.6666666666666667 3.625 6.285714285714286]
[4.0 5.375 2.5384615384615383]
'''
```

上記では，重みづけ平均の分母の総和にもNaNが適用される．
しかし，`np.average()`を用いてしまうと重みにNaNが適用されず，誤った計算となってしまう．

```py:NG!!!
print(np.average(MA, weights=W)) # 88/36
print(np.average(MA, axis=0, weights=W))
print(np.average(MA, axis=1, weights=W))
'''
2.4444444444444446
[1.6666666666666667 2.4166666666666665 2.933333333333333]
[4.0 3.5833333333333335 1.5714285714285714]
'''
```

# おわりに

というわけで，NaN交じりの重みづけ平均にはmasked arrayを用いるのが一つの解であり，その他のNaN交じりの処理もmasked arrayであればほぼ実現可能である．
課題としては，np.arrayを想定したコードに対する導入障壁がある．`data`配列と`mask`配列を同時に保持するために，それぞれに対する挙動を考慮しながらのコーディングにはやや慣れを要する．状況に応じてメリットを生かすのがよい．