トリムして加重移動平均をとる

はじめに

トリムして移動平均をとりたい気分の時がある。さらに直近の数値の影響を大きくしたいので加重平均もとりたい。
基本的に同じ傾向で推移しているけど、たまにグンッてピーキーに数値が上がっちゃう時系列データの時とか、グンッて上がった時の数値は参考外にして移動平均とりたいなぁ。そう思うこともあるはず。
え、無いって？それでも書いちゃう。rolling().applyのところでちょっと苦労したから書き残したい。

トリム平均

平均をとりたいデータ群の数値を大きさ順に並べて、片側、もしくは両側のN%を排除して平均をとる手法。
以下10個のデータがあったとする。
$[10,24,31,34,65,86,87,88,99,101]$

普通に平均をとると、
$(10+24+31+34+65+86+87+88+99+101)\div10=62.50‬$
片側10％・片側10％のデータを除去したとすると、
$(24+31+34+65+86+87+88+99)\div8=64.25$

こんな感じで除去して平均をとるのがトリム(刈り込み)平均。
移動平均だとウィンドウサイズ内の数値に対してトリムして平均をとるイメージ。
メリットは、異常値を排除できる点。平均が異常値の方に引っ張られるのを防ぐことができる。

加重平均

いろんなところで説明が書かれているので詳細は省略。
数値に重みをつけて平均をとるということ。
加重移動平均だと直近の数値の影響が大きくなるように重みをつけることが多い。

今回はトリム加重移動平均をとりたいので、移動平均のウィンドウサイズ内の数値でトリムした後に加重移動平均をとる。

Pythonで実践してみる

パッケージインポート

import.py

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import datetime as dt
from dateutil.relativedelta import relativedelta
sns.set()

データ作成

スマホアプリのアクティブユーザー数の推移とかをイメージにデータを作る

make_data.py

# 減衰曲線
def exp_func(x, a, b):
    return b*(x**a)
x=np.linspace(1,36,36)
data=exp_func(x, -0.5, 100000)

# データフレーム化
df=pd.DataFrame({'x':x,'y':data})
df=df.astype(int)

# 時系列データの想定なので月のカラムを作成
init_mon=dt.datetime(2017,df['x'].values[0],1)
months=[init_mon]
for mon in range(1,len(df)):
    months.append(init_mon + relativedelta(months=mon))
df['month']=months
df.index=months
display(df.head())

# plot
fig=plt.figure(figsize=(12,6))
ax=plt.subplot(1,1,1)
df.plot.bar('month','y',ax=ax)
plt.show()

データの加工（異常値を入れる）

スマホアプリってキャンペーンとかイベント次第でアクティブユーザー数が一時的にかなり上がったりしそうなのでそのような状況を想定してデータの数値を変えてみる。

change_data.py

df2=df.copy()
df2.loc[(df2.index.month==1)&(df2.index.year>=2018), 'y']=df2['y']*1.6
df2.loc[(df2.index.month==2)&(df2.index.year>=2018), 'y']=df2['y']*1.4
df2.loc[(df2.index.month==3)&(df2.index.year>=2018), 'y']=df2['y']*1.2

fig=plt.figure(figsize=(12,6))
ax=plt.subplot(1,1,1)
df2.plot.bar('month','y',ax=ax)
plt.show()

1周年とかのキャンペーンでアクティブユーザー数が上がった感じのデータが完成した。

トリム加重移動平均をとってみる

今回はトリム加重移動平均を予測値として使うという方針でコードを組んでみる。
3か月後の数値の予測値をトリム加重移動平均の値とする。
（例. 2018-08の数値の予測は2018-05時点でトリム加重移動平均をとった値とする）

まずはトリム単純移動平均の関数

sma.py

def sma(roll_list):
    # roll_listにNanがある場合除去
    roll_list=roll_list[~np.isnan(roll_list)]
    # roll_listを昇順に並べる
    sorted_roll_list=sorted(roll_list)
    # 昇順に並べたroll_listの半分の長さを定義
    harf_span=round(len(sorted_roll_list)/2)
    if harf_span > 0:
        # roll_listの中央値以下の数値を取得し、平均をとる
        harf_index=np.where(roll_list < sorted_roll_list[harf_span])
        roll_list_harf=roll_list[harf_index]
        sma = np.sum(roll_list_harf) / len(roll_list_harf)
    else:
        # roll_listの長さが1以下なので中央値が取れない
        # roll_listの値をそのまま使用する
        roll_list_harf=roll_list[0]
        sma = roll_list_harf
    return sma

その次にトリム加重移動平均の関数

sma.py

def wma(roll_list):
    # roll_listにNanがある場合除去
    roll_list=roll_list[~np.isnan(roll_list)]
    # roll_listを昇順に並べる
    sorted_roll_list=sorted(roll_list)
    # 昇順に並べたroll_listの半分の長さを定義
    harf_span=round(len(sorted_roll_list)/2)
    # roll_listの中央値以下の数値を取得
    harf_index=np.where(roll_list < sorted_roll_list[harf_span])
    roll_list_harf=roll_list[harf_index]
    # roll_listの中央値以下の数値の重みを計算し、加重移動平均をとる
    weight = np.arange(len(roll_list_harf)) + 1
    wma = np.sum(weight * roll_list_harf) / weight.sum()
    return wma

次に、上記の関数を使いウィンドウサイズは6か月でそのうち数値が低い3か月分以外をトリムし、加重移動平均をとる。
pandasのrollingを使用すると簡単に移動平均が取れるし、applyを使うと自分で作った関数を適用できる。
ちなみに上記の関数smaとwmaに出てくるroll_listというのはpandasのrollingで指定したウィンドウサイズ(period)で取得されるデータの配列を指す。データの配列は何もしないとSeries型として関数に入れられるっぽい。smaやwmaの関数はndarrayを想定して組んだので、Series型だとエラーとなる。エラーを防ぐためにapplyの引数にraw = Trueを入れてndarray型にしている。

moving_mean.py

# 3か月前の数値のカラムを作る
df2['y_shift'] = df2['y'].shift(3)
# 3か月前の数値のSMA
df2['y_sma'] = df2['y_shift'].rolling(6,min_periods=1).apply(sma, raw = True)
# 3か月前の数値のWMA
df2['y_wma'] = df2['y_shift'].rolling(6,min_periods=1).apply(wma, raw = True)
# WMAの計算ができないNULLをSMA値にする
df2.loc[pd.isna(df2['y_wma']), 'y_wma']=df2['y_sma']

display(df2.head(10))

# plot
fig=plt.figure(figsize=(12,6))
ax=plt.subplot(1,1,1)
df2.plot.bar('month','y',ax=ax,color='b',alpha=0.9)
df2.plot.bar('month','y_wma',ax=ax,color='r',alpha=0.6)
plt.show()

データフレーム黄色線が実際の値、水色線が直近の3か月の数値、緑色線が直近6か月のうち数値が低い3か月分の数値
グラフの青が実際の数値で、赤が3か月前までのデータでトリム加重移動平均をとった数値（予測値）

例えば、2018年6月に注目してみる。
・実際の値：
　$23,570$
・トリムせず3か月分の加重移動平均をとる場合：2018年01月－2018年03月の加重平均
　$(3 \times 30,983 + 2 \times 37,416 + 1 \times 44,376) \div (3+2+1) = 35,359$
・トリムして3か月分の加重移動平均をとる場合：2017年10月－2018年03月の数値が低い3か月分で加重平均
　$(3 \times 30,982 + 2 \times 28,867 + 1 \times 30,151) \div (3+2+1) = 30,139$

トリムして加重移動平均をとった方が2018年01月の急な上昇の影響をその後の計算に与えなくても済むようになる。

まとめ

ちょくちょく異常値が入ってくるような時系列データで、異常値の影響を無くした移動平均をとりたいというときには便利かなと思う。