Zenn
🚜

NumPyで正確な移動平均を取る

2021/06/19に公開
Python
numpy
tech

NumPy.convolveのmode

NumPyには移動平均を取るnumpy.convolveという関数があります。これは正確には畳み込みをする関数なのですが、[1/n, 1/n, ..., 1/n]という窓を使うことで移動平均を取ることができます。

# xx: 1次元のnumpy.array
# n_conv: 何個で移動平均を取るか

b = np.ones(n_conv)/n_conv 
# >> [1/n_size, 1/n_size, ..., 1/n_size]

xx_mean = np.convolve(xx, b, mode="full")

numpy.convolveのmode

convolveのmodeにはfull, same, validの3つが用意されています。

full

フィルタが少しでもかかるパターンは全部計算するオプションです。

出力の要素数はlen(xx)+len(b)-1になります。

same

出力の要素数が入力の要素数と一致するように計算するオプションです。

計算はフィルタが半分だけはみ出した状態からスタートします。終端は要素数を合わせるために半分はみ出たところで停止します。

出力の要素数は長いほうに合わせるためmax(len(xx),len(b))になります。

valid

完璧に計算できる部分だけを出力するオプションです。端ではみ出る部分は無視します。

出力の要素数はabs(len(xx)-len(b))+1になります。

端点付近の挙動とその補正

modeがfullかsameの場合、フィルタがはみ出る区間が存在します。移動平均を考える場合これらの区間は正しい計算がなされません。

シンプルにxx=[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] （長さ10）、窓が6の場合で計算してみます。移動平均も[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] と出力されるのが理想的です。

ですが、残念なことにfullでもsameでもvalidでも理想的な出力は得られません。

xx_full = np.convolve(xx,b,mode="full")
# はみ出ている区間は正しい値でないうえに、長さも変わっている
# array([0.16666667, 0.33333333, 0.5       , 0.66666667, 0.83333333,
#      1.        , 1.        , 1.        , 1.        , 1.        ,
#      0.83333333, 0.66666667, 0.5       , 0.33333333, 0.16666667])

xx_mean = np.convolve(xx,b,mode="mean")
# 長さは同じだがはみ出ている区間はやはり値が狂う
# array([0.5       , 0.66666667, 0.83333333, 1.        , 1.        ,
#      1.        , 1.        , 1.        , 0.83333333, 0.66666667])

xx_valid = np.convolve(xx,b,mode="valid")
# はみ出ている区間はカットされている
# array([1., 1., 1., 1., 1.])

原理は下記のとおりです。まず正しく計算される区間では、出力は

xx_conv[i] = xx[i]*1/6 + xx[i+1]*1/6 + ... + xx[i+6]*1/6

と計算されます。これは正しい平均値の計算になっています。

ところが、はみ出た区間でははみ出た部分は0ですが係数は変わりません。例えば2個分だけはみ出ているところでは有効な4個分の平均値を取るのが理想ですが、分母はあいかわらず6のままです。すなわちその要素は4/6倍と小さな値で出力されます。

補正して理想的な移動平均を出力する

先述の誤差を補正して理想的な移動平均をを得る処理を作ってみます。

まず要素を出し入れすると厄介なので、正しい要素数を返すsameオプションをベースにします。

sameの場合、両端からlen(b)//2ずつはみ出ます。両端から要素i個分離れたところでは(i+1+size//2)/len(b)倍だけ小さな値になっているので、その逆数をかけて補正します。

下記に補正を合わせて理想的な移動平均を返す関数を書いておきます。

# xxに対してsize個での移動平均を取る
def valid_convolve(xx, size):
    b = np.ones(size)/size
    xx_mean = np.convolve(xx, b, mode="same")

    n_conv = math.ceil(size/2)

    # 補正部分
    xx_mean[0] *= size/n_conv
    for i in range(1, n_conv):
        xx_mean[i] *= size/(i+n_conv)
        xx_mean[-i] *= size/(i + n_conv - (size % 2)) 
	# size%2は奇数偶数での違いに対応するため

    return xx_mean

参考

Discussion