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

はじめに

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

通常の配列（NaNなし）

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

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()である．関数および配列自体のメソッドとして使用が可能である．

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()と同様，軸を指定した平均を算出することができる．

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

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

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()によって型変換を行う．

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を返す．

print(np.mean(NA)) # np.nan
print(NA.mean()) # np.nan

np.nanmean()

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

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配列を重ねた形となっている．

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()を使用しても結果は全く変わらない．

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()を用いなければならない．

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が適用されず，誤った計算となってしまう．

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配列を同時に保持するために，それぞれに対する挙動を考慮しながらのコーディングにはやや慣れを要する．状況に応じてメリットを生かすのがよい．