1. 前回までのまとめ
少しシリーズが長くなってきたので、過去記事をどのようなモチベーションで書いてきたかを整理します。
syleir.hatenablog.com
分割時系列解析とはなにか？どのようなモチベーションでこの分析がやりたいのか？モデリングを提示し、基本的な考え方や求めたいものについてまとめました。
syleir.hatenablog.com
今までITSはRやStataでばかり実装されてきており、Pythonでの実装は確認できなかったので、Statmodelsというライブラリを用いて、折れ線回帰、セグメント回帰をPythonで実装し、グラフの描画をしました。
syleir.hatenablog.com
介入前後での傾きを求めるモデルについて、式変形を交えながら解説し、モデリングと結果の解釈について解説しました。
2. 今回の目標
時系列データには少なからず周期性が反映されます。季節性や周期性をもつ時系列データに対し、今までのような直線での折れ線回帰を用いた場合、当てはまりが悪いことがあります。そのような場合にどのような処理をしてITSを実装していくのかをPythonを使いながら説明していきます。
3. 目次
5. 周期回帰分析について
別名をFourier analysisとか、Harmonic analysisとか言います。今回の記事の肝です。
フーリエ級数展開を利用して周期関数を回帰分析しようというモチベーションです。
周期性とフーリエ級数、フーリエ変換の相性が良いのは過去のこの辺の記事でも触れています。
フーリエ級数展開というのは
周期の関数
がフーリエ級数展開できるとき、（実務上大体の関数はできます。できるための条件は専門書を参照ください。）
と無限級数展開することを言います。
周期回帰分析では、この無限級数展開を途中で打ち切り、残りを誤差項と考え、実際のデータからを推定するものになります。
例えば、n=1で切るとすると、を推測することができれば、１次のフーリエ級数で近似した予測関数が得られることになります。何次の項まで考えるかはハイパーパラメータで、分析の時にデータを見ながら設計することになります。この係数も通常の回帰分析と同じように解析的に計算することができるので、周期関数を回帰分析で推測できるようになっています。
6. 今回のモデリング
ただしh(T)はフーリエ級数展開できる周期関数：つまり上で定義した関数とします。
今までのモデリングに加えて、h(T)が追加されています。この周期関数でのモデリングを追加することで周期性に対してrobustなモデルになります。
7. Rでのharmonic関数
Rは比較的簡単にharmonic関数という関数を用いることで簡単にフーリエ級数を用意できます。
harmonic(x, nfreq, period, intercept)というような引数を取りますが、
x:データ間隔を表す引数、月なのか日なのか
nfreq: １期あたりに含まれる周期の数（上の式でいう）
period: １期あたりのデータ点（上の式でいう）
intercept: 切片の有無(True: 上の式で、False:
)
とすれば良いです。例えば年単位の周期性を考えたく、データが1月に1つの今回のデータ（12データで1年）において、2次までのフーリエ級数展開の式を用意したければ、harmonic(month,2, 12,False) とすれば良いです。
8. PythonにおけるRでのharmonic関数に相当する実装
Pythonでは関数一つでこれを解決する実装は（少なくとも自分の実力・調査範囲では）難しいです。より直接的にのデータ列を用意していくしかありません。ただ、DeterministicProcessというライブラリを使うことで少し楽ができます。これはデータ列を生成するためのライブラリで、季節ごとのダミー変数とか、月毎のダミー変数とか、今回利用するフーリエ級数列を用意することができる便利なライブラリです。
これを2次までのフーリエ級数を求める方針で実装してみると下のようになります。
#ライブラリの読み込み from statsmodels.tsa.deterministic import DeterministicProcess, Fourier #周期12, 2次までのフーリエ級数を用意する fourier = Fourier(period=12, order = 2) #データをDeterministicProcessにより生成 dp = DeterministicProcess( index=data.index, constant=False, order=0, drop=True, additional_terms=[fourier]) H = dp.in_sample()
pythonでのFourier関数はFourier(period, order)という引数を取ります。
period： 1周期の長さ()に相当します。Rのharmonic()におけるperiodと全く同じです。
order: フーリエ成分の数() Rのnfreqと全く同じです。
xはありませんが、前処理で調節するか、DeterministicProcessにおいてindexを目的のdatetime型に設定することで実現可能です。おそらくコードを読んでいても難しいと思うので結果を載せます。Hを最初の方だけ出力するとこうなります。
試しにH['sin(1,12)']を描画してみると次のようになります。
plt.figure(figsize=(16,8)) plt.plot(H['sin(1,12)']) plt.xlabel('time',fontsize=18) plt.xticks(np.arange(0, 48, 12)) plt.ylabel('y', fontsize=18)
ちゃんととなっていることがわかります。
これらのデータ列の線形結合がharmonic関数に相当するので、これらを変数に追加して今までのモデリングに適応してみれば良いです。
9. 変数の抽出
## 研究開始時点から経過した時間T t = np.arange(0, (len(data))) data["t"] = t ## 処置前後を表す二値変数Xt x = np.where(data['date']>="2022-01-01", 1, 0) data["xt"]=x ### 処置タイミングTiを取得、介入からの経過時間Xt(T-Ti) ti = 24 data['xt(t-ti)'] = data['xt']*(data['t']-ti) #変数の抽出 X1 = data[['t','xt','xt(t-ti)']] H = dp.in_sample()[["sin(1,12)", "cos(1,12)","sin(2,12)","cos(2,12)"]] y = data[['y']] #X1とHを統合 X = X1.join(H)
これでデータの準備ができました。現在Xはこのようになっています。
モデルの作成と学習
ここからは前回記事とほぼ同じです。適宜省略します。
# 回帰モデル作成 mod = sm.OLS(y, sm.add_constant(X)) # 訓練 result = mod.fit() # 統計サマリを表示 print(result.summary())
このように分析結果を出力することができます。原著論文では2次まで採用していたので2次で実装していますが、結果を見てみると、sin(2, 12) , cos(2, 12) は有意な変数ではないので1次までのモデリングでよかったかもしれません。
10. 推測を行う
#Xの0列目に'cep'というXと同じ長さの[1, 1, 1, ...]という列を挿入 X.insert(0, 'cep', [1]*len(X)) #学習したモデルから推測を行う y_predict = result.predict(X)
ここで、推測してみたデータを描画してみます。
#グラフを描画 plt.figure(figsize=(16,8)) plt.title('ITS_test') plt.plot(y, 'o') plt.plot(y_predict) plt.xlabel('time',fontsize=18) plt.ylabel('event', fontsize=18)
赤実線が予測した線、青点が実際のデータですが、いい感じに学習できていることがわかります。
11. 反事実データの作成、推測
#反事実の作成 cf_X = X.copy(deep=True) cf_X['xt'] = 0 cf_X['xt(t-ti)'] = 0 #反事実を予測する X_cf = cf_X[['t','xt','xt(t-ti)','sin(1,12)','cos(1,12)','sin(2,12)','cos(2,12)']] X_cf.insert(0, 'cep', [1]*len(X_cf)) y_cf_predict = result.predict(X_cf) #事実データに対する信頼区間を返す predict = result.get_prediction(X) frame = predict.summary_frame(alpha=0.05) y_predict_l = frame['mean_ci_lower'] y_predict_u = frame['mean_ci_upper'] #反事実データに対する信頼区間を返す predict_cf = result.get_prediction(X_cf) frame_cf = predict_cf.summary_frame(alpha=0.05) y_predict_cf_l = frame_cf['mean_ci_lower'] y_predict_cf_u = frame_cf['mean_ci_upper']
12. 描画
###描画 plt.style.use('fivethirtyeight') #itv（intervention time）介入したところが何番目のデータか #tiは上で定義済み（介入した時間） itv = 24 plt.figure(figsize=(16,8)) plt.title('ITS_test') plt.plot(y_predict) plt.plot(y_cf_predict[(itv-1):], 'r--') plt.plot(y, 'o') plt.fill_between(data["date"], y_predict_l, y_predict_u, where= data["date"] >= "2022-01-01", color='b', alpha=0.2) plt.fill_between(data["date"], y_predict_cf_l, y_predict_cf_u, where= data["date"] >= "2022-01-01", color='r', alpha=0.2) plt.xlabel('date',fontsize=18) plt.xticks(np.arange(0, 48, 12)) plt.ylabel('event', fontsize=18)
出力がこれです。今までのモデルと比較してみましょう。
今までの折れ線回帰、セグメント回帰での出力はこれですから、相当当てはまりが改善していることがわかります。
pre-postの傾きの変化のモデルも前回の式変形の両辺にharmonic項を足すだけなので前回と同様に変数を用意して作ることができます。やってみてください。
おわりに
Rでのharmonic関数に相当する実装をPythonで行い、折れ線回帰分析に周期回帰分析を組み合わせたモデリングを行うことで周期性・季節性を考慮したITSがPythonで実装できました。ある程度は実務に耐えられるレベルまで昇華できた気がしますがいかがでしょうか。
少し難しかったでしょうか。過去の記事も見ながらじっくり取り組んでみてください。過去にはあまりフーリエ解析の記事は書けていないので基礎があやふやな方はこの辺の本を読んでみてください。統計検定準1級レベルでもたまに出るのでこの辺は押さえておくと良いです。
フーリエ解析の初歩の初歩を教えてくれる良書です。漫画だからといってナメてたんですけど、初学者にはわかりやすかったのでおすすめです。次にフーリエ解析について学ぶならこの本がおすすめです。どう実務に生かしていくかを含めてわかります。
全実装
#ライブラリのインポート import pandas as pd import numpy as np import matplotlib.pyplot as plt import statsmodels.api as sm from statsmodels.tsa.deterministic import DeterministicProcess, Fourier #ファイルをDataFrameとして読み込む data = pd.read_csv("dummy_data.csv") #周期12, 2次までのフーリエ級数を用意する fourier = Fourier(period=12, order = 2) #データをDeterministicProcessにより生成 dp = DeterministicProcess( index=data.index, constant=False, order=0, drop=True, additional_terms=[fourier]) H = dp.in_sample() ## 研究開始時点から経過した時間T t = np.arange(0, (len(data))) data["t"] = t ## 処置前後を表す二値変数Xt x = np.where(data['date']>="2022-01-01", 1, 0) data["xt"]=x ### 処置タイミングTiを取得、介入からの経過時間Xt(T-Ti) ti = 24 data['xt(t-ti)'] = data['xt']*(data['t']-ti) #変数の抽出 X1 = data[['t','xt','xt(t-ti)']] H = dp.in_sample()[["sin(1,12)", "cos(1,12)","sin(2,12)","cos(2,12)"]] y = data[['y']] #X1とHを統合 X = X1.join(H) # 回帰モデル作成 mod = sm.OLS(y, sm.add_constant(X)) # 訓練 result = mod.fit() # 統計サマリを表示 print(result.summary()) #Xの0列目に'cep'というXと同じ長さの[1, 1, 1, ...]という列を挿入 X.insert(0, 'cep', [1]*len(X)) #学習したモデルから推測を行う y_predict = result.predict(X) #反事実の作成 cf_X = X.copy(deep=True) cf_X['xt'] = 0 cf_X['xt(t-ti)'] = 0 #反事実を予測する X_cf = cf_X[['t','xt','xt(t-ti)','sin(1,12)','cos(1,12)','sin(2,12)','cos(2,12)']] X_cf.insert(0, 'cep', [1]*len(X_cf)) y_cf_predict = result.predict(X_cf) #事実データに対する信頼区間を返す predict = result.get_prediction(X) frame = predict.summary_frame(alpha=0.05) y_predict_l = frame['mean_ci_lower'] y_predict_u = frame['mean_ci_upper'] #反事実データに対する信頼区間を返す predict_cf = result.get_prediction(X_cf) frame_cf = predict_cf.summary_frame(alpha=0.05) y_predict_cf_l = frame_cf['mean_ci_lower'] y_predict_cf_u = frame_cf['mean_ci_upper'] ###描画 plt.style.use('fivethirtyeight') #itv（intervention time）介入したところが何番目のデータか #tiは上で定義済み（介入した時間） itv = 24 plt.figure(figsize=(16,8)) plt.title('ITS_test') plt.plot(y_predict) plt.plot(y_cf_predict[(itv-1):], 'r--') plt.plot(y, 'o') plt.fill_between(data["date"], y_predict_l, y_predict_u, where= data["date"] >= "2022-01-01", color='b', alpha=0.2) plt.fill_between(data["date"], y_predict_cf_l, y_predict_cf_u, where= data["date"] >= "2022-01-01", color='r', alpha=0.2) plt.xlabel('date',fontsize=18) plt.xticks(np.arange(0, 48, 12)) plt.ylabel('event', fontsize=18)
