Prophet?ARIMA?外れ値検出?

安倍政権による北朝鮮ミサイル打ち上げ関与の陰謀論をデータから検証してみた-Qiita

トンデモ記事かとタイトルに釣られて読んだけど,面白かった.趣旨自体はまあ「ネタ」だからどうでも良いけど(陰謀論は所詮陰謀論),重要なのは「Prophet」モデルが妥当かどうかじゃないだろうか.内閣支持率の様な時系列データに於いて,時間の流れというのはかなり重要だと思うけど,「Prophet」の様な平滑化モデルはそれを無視して,空間的な変化だけを捉える.それ自体が悪い事では無いし,必ずしも計量時系列的な分析(例えばARIMAモデル)の方が良いかどうかは分からないけど,例えば,周期要因を除外して,ミサイル打ち上げという突発的な外生的要因がどれだけ影響しているかをみる場合は,平滑化モデルでも良いのかもしれないけど,少なくともこの緑の線はぐにゃぐにゃしすぎの様な.「Prophet」は使った事が無いので,もしかしたら盛大な勘違いをしているかもしれないけど,Seasonal Decompositionは腐る程利用してきたけども,これって,要はSeasonal Decompositionと同じ話だと思う.で,そうしてみると,周期性変動をちゃんと捉えられていないというか,周期が合っていない気がする(追記:トレンドをみているのではなくて,予測応答をみているっぽいので,周期性変動が含まれているのは当然か).

平滑化モデルって何にでも当て嵌められる反面,本当にそれで尤もらしいと云えるのかが難しい様な.
(「Prophet」は特段難しい事を考えなくても割と尤もらしい結果が得られるので有名らしいけど)

この場合,全系列でみて変動がみられないって事が何より重要で,強いて「外れ値」としてみるなら,単にone-class SVMなりで外れ値検出をして,ミサイル発射が有意に影響を及ぼすなら外れ値としてみられると考えてその関係性をみれば良いんじゃないかな.

30分で出来る所まで.ミサイルに関するデータは,登録しないといけないっぽいので記事をみながら合わせてみる(後で「Prophet」も試してみたい).

%matplotlib inline
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import changefinder
from sklearn.ensemble import IsolationForest


df = pd.read_csv('http://www.sekai-kabuka.com/shiji/data.csv')
df.head()
df.index = pd.to_datetime(df.Date)
df.drop(['Date', 'Unnamed: 3'], axis=1, inplace=True)
df['Ratio'] = df.Yes / (df.Yes + df.No) * 100
df.Ratio['2016':].plot()

ダウンロード

・ChangeFinderARIMA(ARIMAモデルを当て嵌めて変化点検出)

cf = changefinder.ChangeFinderARIMA(term=30, smooth=7, order=(1, 0, 0))

ret = []
for i in df.Ratio['2016':].values:
    score = cf.update(i)
    ret.append(score)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(ret)
ax2 = ax.twinx()
ax2.plot(df.Ratio['2016':].values,'r')
plt.show()

ダウンロード (1)

・IsolationFrestによる外れ値検出

from scipy import stats
from sklearn.preprocessing import MinMaxScaler

mms = MinMaxScaler()
norm = mms.fit_transform(df.Ratio['2016':].values)
clf = IsolationForest()
clf.fit(norm.reshape(norm.size, -1), norm.ravel())
scores_pred = clf.decision_function(norm.reshape(norm.size, -1))
mms_scores_pred = mms.fit_transform(scores_pred)
thresholdover = [stats.scoreatpercentile(mms_scores_pred, 20) for _ in range(len(scores_pred))]
thresholdunder = [stats.scoreatpercentile(mms_scores_pred, 50) - stats.scoreatpercentile(mms_scores_pred, 20) for _ in range(len(scores_pred))]

scoredf = pd.DataFrame({'under': thresholdunder, 'over': thresholdover}, index = pd.date_range(start='2016', periods=len(scores_pred), freq='d'))
normdf = pd.DataFrame({'norm Ratio': norm}, index = pd.date_range(start='2016', periods=len(scores_pred), freq='d')) 
fig = plt.figure()
ax = fig.add_subplot(111)
scoredf.plot(ax=ax, c='r')
normdf.plot(ax = ax)

ダウンロード (2)

急ぎ足なので色々乱雑だけど,2016年4月頃の支持率が急落している箇所以外は,概ね平常運転と云えるのではないか.ChangeFinderARIMAで変化点として検出されている箇所と,ミサイル発射のタイミングは,ちょっと時間をかけないと何とも言えないけど,まあミサイルが原因というよりは様々な内生的,外生的要因による変動はもちろんあるよね,程度にみえる.

ほんの30分でも(コードは汚いけど)色々な事はできる.でも,本当に結果が尤もらしいのかどうか,その検証にこそ時間が掛かるのであって,手法の当て嵌め,解析自体はPythonでも一瞬だし,でも例えばARIMAモデルの次数は妥当かとか(本当は(p, d, q)の推定をするべきだろうし),外れ値検出も明確な仮定とそれに基づくモデルがあってこそだろうし(trainデータとして何を正常値とするか,どこに閾値をおくか),平滑化モデルは時間を無視する分,空間に対してどういう意図を持つのか,モデリング時のデザインが重要で,記事のProphetって余りにもぐにゃぐにゃしていてきちんとトレンド取れているのかな?っていう疑問.こういうのをみて思うのは,数理的妥当性と説得力の2つが大事で,数理的妥当性は数学的合理性や統計学的合理性の話でより学術的な話だけど,実務的には説得力が結局全てで,じゃあこれって説得力あるんだろうか?

この場合,変に解析するよりも,生のデータをみせた方が説得力ってあるんじゃないのかな.これってあまり良いプロモーションだと思えない.

追記(Prophetについて):
Windowsだと「pip install fbprophet」ではエラーが出て失敗する(pystan側の問題?色々な所で質問や回答が上がっているけど,長すぎて読み切れない).「conda install -c conda-forge fbprophet」ならすんなり通るので,大人しくそれで.

%matplotlib inline
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from fbprophet import Prophet


df = pd.read_csv('http://www.sekai-kabuka.com/shiji/data.csv')
df.index = pd.to_datetime(df.Date)
df.drop(['Date', 'Unnamed: 3'], axis=1, inplace=True)
df.Ratio = df.Yes / (df.Yes + df.No) * 100
df2 = pd.DataFrame({'y': df.Ratio['2016':], 'ds': df.Ratio['2016':].index})

prpht = Prophet()
prpht.fit(df2)
future = prpht.make_future_dataframe(periods=365)
prpht_pred = prpht.predict(future)
prpht_pred.index = prpht_pred.ds
prpht_pred.drop('ds', axis=1, inplace=True)

fig = plt.figure()
ax = fig.add_subplot(111)
prpht_pred.trend.plot(ax=ax, c='r')
ax.fill_between(prpht_pred.index, prpht_pred.trend_lower, prpht_pred.trend_upper, color='r', alpha=0.2)
df.Ratio['2016':].plot(ax = ax)

ダウンロード (3)

fig = plt.figure()
ax = fig.add_subplot(111)
prpht_pred.yhat.plot(ax=ax, c='r')
ax.fill_between(prpht_pred.index, prpht_pred.yhat_lower, prpht_pred.yhat_upper, color='r', alpha=0.2)
df.Ratio['2016':].plot(ax = ax)

ダウンロード (5)

ちょっと勘違いしていた.Seasonal Decompositonと基本的な考え方は一緒だけど,プロセスは全然違って,あれは周期性を仮定して,トレンドと周期性変動を求めるが(なので前提がすごく大事),ProphetはMCMCに基づいて平均誤差(幾つものランダムなフィッティング区間と予測区間を於いて,それぞれの誤差の平均が最小になる様な感じ)を最小化する様な多項式回帰モデルを求めて,非常にシンプルな方法で尤もらしい予想をしようっていう話の様だ.だから,何も考えずにフィッティングできる.凄く手軽に使えるので,(厳密な話はともかく何かしらビジュアライゼーションしたい時に)確かにこれは便利そう.

future2 = prpht.make_future_dataframe(periods=0)
prpht_pred = prpht.predict(future2)
prpht.plot(prpht_pred)
prpht_pred.index = prpht_pred.ds
prpht_pred.drop('ds', axis=1, inplace=True)

fig = plt.figure()
ax = fig.add_subplot(111)
prpht_pred.yhat.plot(ax=ax, c='r')
ax.fill_between(prpht_pred.index, prpht_pred.yhat_lower, prpht_pred.yhat_upper, color='r', alpha=0.2)
df.Ratio['2016':].plot(ax = ax)

ダウンロード (7)

ダウンロード (4)

記事と同様のグラフ.てっきり,トレンドをみて検討しているのかと思っていたけど,モデルの予測応答をみていたからギザギザだったのか.更に云えば,元々のデータがかなりノイジーなので,本当はリサンプリングしたり,スムージングしたり,プリパレーションした方が良い.

ただ,これだとやっぱり,外れ値の検討,或いは外生的要因の影響について,あまり尤もらしい議論ではないだろう.明らかに影響がみられる部分について,「影響が無いとは云えない」位には云えるかも知れないけど.

Prophetに基づいて内閣支持率のトレンドを求め,それが真値(母比率)の推定(マクロにみた世論の動向)であると仮定する.また,このデータ(内閣支持率)は,N=1,000程度だと考えられるので,母比率の信頼区間95%を考えると,

fig = plt.figure()
ax = fig.add_subplot(111)
prpht_pred.trend.plot(ax=ax, c='r')
ax.fill_between(prpht_pred.index, prpht_pred.trend-3.1, prpht_pred.trend+3.1, color='r', alpha=0.2)
df.Ratio['2016':].plot(ax = ax)

ダウンロード (8)

この方がまだ尤もらしいんじゃないだろうか.

・その他
ダウンロード (6)

ダウンロード (10)

ダウンロード (9)

この手の手法は,パラダイムシフトとか外生的要因が多く存在するとか,周期性変動が容易に変わり得るデータ(要は複雑系)には向いていない.

ダウンロード (1)

カテゴリー: 未分類 パーマリンク

Prophet?ARIMA?外れ値検出? への3件のフィードバック

  1. ピンバック: 時系列データの変化点検出 | 粉末@それは風のように (日記)

  2. ピンバック: TensorFlowで線形分類器による分類 | 粉末@それは風のように (日記)

  3. ピンバック: お山を探して間隔をみる | 粉末@それは風のように (日記)

コメントは受け付けていません。