pythonで正確な予測を実現したい

カルマンフィルタを用いての予測

一年間のバスの利用者数のデータを用いて翌年1年間の利用者数を予測したいと考えています。
そこでカルマンフィルタを用いて予測を行ったのですが、季節成分のみでは正確な予測が出来ませんでした。
そこで、利用者数が圧倒的に少ない日（お盆休み、、など）を抽出してダミー変数化してモデルに入れたいと考えていますが、どのように組み込めばいいでしょうか。このサイトを参考にしています。→https://logics-of-blue.com/python-state-space-models/
初心者で、分かりにくいと思いますがよろしくお願いいたします。

発生している問題・エラーメッセージ

該当のソースコード

#---ライブラリのインポート
import pandas as pd
import matplotlib.pyplot as plt

# グラフを横長にする
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6

# 統計モデル
import statsmodels.api as sm

df = pd.read_csv('バス利用者.csv', encoding = "shift-jis")

df['ymd'] = pd.to_datetime(df['ymd'])
df['B_Year'] = df['ymd'].dt.year
df['B_Month']= df['ymd'].dt.month.fillna(0.0).astype(int)
df['B_Day'] = df['ymd'].dt.day.fillna(0.0).astype(int)
df['B_Hour'] = df['ymd'].dt.hour
df['B_Minute'] = df['ymd'].dt.minute
df['B_Date'] = df['B_Year'].astype(str).str.zfill(4).astype(str) + df['B_Month'].astype(str).str.zfill(2).astype(str) + df['B_Day'].astype(str).str.zfill(2).astype(str)
df['B_DayofWeek'] = df['ymd'].dt.dayofweek.fillna(-1.0).astype(int)#曜日

df['KosyaTm'] = pd.to_datetime(df['KosyaTm'])
df['A_Month']= df['KosyaTm'].dt.month
df['A_Day'] = df['KosyaTm'].dt.day
df['A_Hour'] = df['KosyaTm'].dt.hour
df['A_Minute'] = df['KosyaTm'].dt.minute

B_BSList = BSList.rename('JosyaTeiNm')
A_BSList = BSList.rename('KosyaTeiNm')

df = pd.merge(df, B_BSList,on='JosyaTeiNm',how='left', indicator=True)

#---ODペアの付与
df['ODPair'] = df['JosyaTeiNm'] + '-' + df['KosyaTeiNm']

df_1 = df.groupby(['B_Date']).count()['ODPair'].reset_index()
df_1['B_Date'] = pd.to_datetime(df_1['B_Date'], format='%Y%m%d')
df_1 = df_1.set_index('B_Date')#B_Dateをインデックスに


# 日付形式にする
ts = df_1['ODPair'] 
ts.head()

# プロット
plt.plot(ts)



################################
# ローカルレベルモデルの推定
################################
mod_local_level = sm.tsa.UnobservedComponents(ts, 'local level')

# 最尤法によるパラメタの推定
res_local_level = mod_local_level.fit()

# 推定されたパラメタ一覧
print(res_local_level.summary())

# 推定された状態・トレンドの描画
rcParams['figure.figsize'] = 15, 15
fig = res_local_level.plot_components()
fig2 = res_local_level.plot_diagnostics()


################################
# ローカル線形トレンドモデル
################################
mod_trend = sm.tsa.UnobservedComponents(ts,'local linear trend')

# 最尤法によるパラメタの推定
# ワーニングが出たのでBFGS法で最適化する
res_trend = mod_trend.fit(method='bfgs')

# 推定されたパラメタ一覧
print(res_trend.summary())

# 推定された状態・トレンドの描画
rcParams['figure.figsize'] = 15, 20
fig = res_trend.plot_components()
fig2 = res_trend.plot_diagnostics()

################################
# 季節変動ありのローカルレベルモデル
################################
mod_season_local_level = sm.tsa.UnobservedComponents(
    ts,
    'local level',
    seasonal = 6
    
)

# まずはNelder-Meadでパラメタを推定し、その結果を初期値としてまた最適化する。2回目はBFGSを使用。
res_season_local_level = mod_season_local_level.fit(
    method='bfgs', 
    maxiter=500, 
    start_params=mod_season_local_level.fit(method='nm', maxiter=500).params,
)

# 推定されたパラメタ一覧
print(res_season_local_level.summary())

# 推定された状態・トレンド・季節の影響の描画
rcParams['figure.figsize'] = 15, 20
fig = res_season_local_level.plot_components()
fig2 = res_season_local_level.plot_diagnostics()

################################
# 季節変動ありのローカル線形トレンドモデル
################################
mod_season_trend = sm.tsa.UnobservedComponents(
    ts,
    'local linear trend',
    seasonal = 6
 
)

# まずはNelder-Meadでパラメタを推定し、その結果を初期値としてまた最適化する。2回目はBFGSを使用。
res_season_trend = mod_season_trend.fit(
    method='bfgs', 
    maxiter=500, 
    start_params=mod_season_trend.fit(method='nm', maxiter=500).params,
)

# 推定されたパラメタ一覧
print(res_season_trend.summary())

# 推定された状態・トレンド・季節の影響の描画
rcParams['figure.figsize'] = 15, 20
fig = res_season_trend.plot_components()
fig2 = res_season_trend.plot_diagnostics()

################################
# 季節変動ありのローカル線形トレンドモデル
# ただし、トレンドの分散は無し
################################
mod_season_trend_d = sm.tsa.UnobservedComponents(
    ts,
    'local linear deterministic trend',
    seasonal = 6
    
)

# まずはNelder-Meadでパラメタを推定し、その結果を初期値としてまた最適化する。2回目はBFGSを使用。
res_season_trend_d = mod_season_trend_d.fit(
    method='bfgs', 
    maxiter=500, 
    start_params=mod_season_trend_d.fit(method='nm', maxiter=500).params,
)

# 推定されたパラメタ一覧
print(res_season_trend_d.summary())

# 推定された状態・トレンド・季節の影響の描画
rcParams['figure.figsize'] = 15, 20
fig = res_season_trend_d.plot_components()
fig2 = res_season_trend_d.plot_diagnostics()

################################
# 季節変動ありのローカル線形トレンドモデル
# ただし、トレンドと観測誤差の分散は無し
################################
mod_season_rw = sm.tsa.UnobservedComponents(
    ts,
    'random walk with drift',
    seasonal = 6
    

)

# まずはNelder-Meadでパラメタを推定し、その結果を初期値としてまた最適化する。2回目はBFGSを使用。
res_season_rw = mod_season_rw.fit(
    method='bfgs', 
    maxiter=500, 
    start_params=mod_season_rw.fit(method='nm', maxiter=500).params,
)

# 推定されたパラメタ一覧
print(res_season_rw.summary())

# 推定された状態・トレンド・季節の影響の描画
rcParams['figure.figsize'] = 15, 20
fig = res_season_rw.plot_components()
fig2 = res_season_rw.plot_diagnostics()

################################
# 今まで計算してきたモデルのAICを格納する
################################
aic_list = pd.DataFrame(index=range(6), columns=["model", "aic"])

aic_list.loc[0]["model"] = "res_local_level"
aic_list.loc[0]["aic"] = res_local_level.aic

aic_list.loc[1]["model"] = "res_trend"
aic_list.loc[1]["aic"] = res_trend.aic

aic_list.loc[2]["model"] = "res_season_local_level"
aic_list.loc[2]["aic"] = res_season_local_level.aic

aic_list.loc[3]["model"] = "res_season_trend"
aic_list.loc[3]["aic"] = res_season_trend.aic

aic_list.loc[4]["model"] = "res_season_trend_d"
aic_list.loc[4]["aic"] = res_season_trend_d.aic

aic_list.loc[5]["model"] = "res_season_rw"
aic_list.loc[5]["aic"] = res_season_rw.aic

# 結果の表示
aic_list

試したこと

補足情報（FW/ツールのバージョンなど）

ここにより詳細な情報を記載してください。

coffeebar

2021/12/21 00:24

ryotaro18782さんがやりたいのは『「休日」カラムを作って休日は1、それ以外は0として、説明変数に加える』ということでしょうか？　それなら、それで問題ないと思いますが。それとも「休日」の抽出方法をお聞きになっている？

ryotaro18782

2021/12/21 05:05

休日カラムを作って説明変数を加えるということです！一年分のデータしかないので、季節変動は予測に使えないと思い、休日をダミー変数としていれたいと考えています。分かりにくくて申し訳ないです。一応休日の抽出はモジュールとしてこんな感じで作っているのですがこのコードは変数としてモデルに入れることは可能でしょうか？？ def CheckHoliday(daytime): #--- #---祝日のリスト作成 #--- holidays = {} holidays[2020,1,1] = '元日' holidays[2020,1,2] = '正月休み' holidays[2020,1,3] = '正月休み' year = daytime.year month = daytime.month day = daytime.day if (year,month,day) in holidays: result = 1 else: result = 0 return result def CheckBizDay(daytime):#平日であれば1，土日祝日であれば0を返す関数 if daytime.dayofweek >=5 | CheckHoliday(daytime) ==1: return 0 else: return 1

coffeebar

2021/12/21 07:58

質問の意図がよく分かりません。「モデルに入れることは可能でしょうか？？」と言われても、入れればいいんじゃないですか？としか。「入れてみたらエラーが出て、その対処法が分かりません。」「入れてみたら精度が下がったのですが、入れ方に間違いがありますか？」とかいう質問なら、答えようもありますが。

ryotaro18782

2021/12/22 04:58

質問が抽象的で申し訳ありません。入れ方が分からなかったので質問した所存です。このような形で入れたいと思っていて、エラーは出ずにコードは回るのですが、cond_holiday = trueの部分（ダミー変数化した部分）が反映されていないです。この対処法がわからないので教えていただきたいです。度々申し訳ありません。 cond_holiday = (df['B_HolidayBizDay'] == 1) df.loc[cond_holiday ,'B_HolidayBizDay'].mean() mod_season_rw = sm.tsa.UnobservedComponents( ts, 'random walk with drift', cond_holiday = true )