お弁当の需要予測２

今回はSIGNATEの「お弁当需要予測」の続きをやっていきます！
ds-blog.tbtech.co.jp

チュートリアル

このお弁当需要予測のデータセットには、チュートリアルが投稿されていました。
SIGNATEのデータセットでチュートリアルがあるのは珍しいと思いますが
活用しない手は無いので！チュートリアルを抜粋して見ていきます。
signate.jp

データの内容を再掲します。

データの可視化

まずチュートリアルでは、データを読み込んだ後
インデックス（表の左端にある0,1,2,,の部分）を日付データである「datetime」に変更しています。

train.index = pd.to_datetime(train["datetime"])

欠損値処理

次に欠損値処理を行います。
[1]「payday」の欠けた部分を全て0で補完します。
[2] （apply：行や列ごとにまとめて関数を適用してくれる ／ lambda：関数を簡易的に定義できる）
　「precipitation」で、もし値ｘが「--」なら x=-1 、それ以外はそのままで型をfloatにします。
[3][4] 「event」と「remarks」の欠けた部分全てを「なし」で補完します。
[5] 「datatime」をx.splitで「-」で分解した[1]は月になるので、[”month”]に入れています。

train["payday"] = train["payday"].fillna(0)
train["precipitation"] = train["precipitation"].apply(lambda x : -1 if x == "--" else float(x))
train["event"] = train["event"].fillna("なし")
train["remarks"] = train["remarks"].fillna("なし")
train["month"] = train["datetime"].apply(lambda x : int(x.split("-")[1]))

販売数ｙの推移

続いて販売数の推移に着目します。
「ｙ」の推移をプロットしてみます。

train["y"].plot(figsize=(15,4))

日数が経つにつれ販売数が減っていることがわかります。
また、チュートリアルでは「落ちて行っているのにも関わらず、スパイクしている為、
売り上げ数に寄与している何かしらのファクターがある？」とあります。
確かに、落ちていっている中にも販売数が多くなる日がいくつも見られます。

箱ひげ図

説明変数に対する販売数「ｙ」の分布を、箱ひげ図で可視化します。
[1] 複数グラフを表示するので予めfigとaxを定義します。２×２個のグラフを表示します。
[2]~[4],[6] sns.boxplotで箱ひげ図を表示します。axではグラフの場所を指定します。
[5] 「remarks」の座標値が重ならないよう、斜めに表示させています。
[7] 座標値が重ならないよう調整しています。

fig, ax = plt.subplots(2,2,figsize=(12,7))
sns.boxplot(x="week",y="y",data=train,ax=ax[0][0])
sns.boxplot(x="weather",y="y",data=train,ax=ax[0][1])
sns.boxplot(x="remarks",y="y",data=train,ax=ax[1][0])
ax[1][0].set_xticklabels(ax[1][0].get_xticklabels(),rotation=30)
sns.boxplot(x="event",y="y",data=train,ax=ax[1][1])
plt.tight_layout()

顕著に変化が現れたのは、「remarks」が「お楽しみメニュー」と「なし」の値でした。
「お楽しみメニュー」の日はそうでない日に比べ、販売数が多いようです。

この結果から、「お楽しみメニュー」の日を除いた日の販売数の推移を確認します。

train[train["remarks"]!="お楽しみメニュー"]["y"].plot(figsize=(15,4))

先程のグラフより、なだらかになりました！「お楽しみメニュー」の影響がありそうです。

お楽しみメニューの分析

更に「お楽しみメニュー」の中身にも着目してみます。
「お楽しみメニュー」の日のみの販売数の推移を表示します。

train[train["remarks"]=="お楽しみメニュー"]["y"].plot(figsize=(15,4))

同じ「お楽しみメニュー」でも販売数に差があることがわかります。

そこで「お楽しみメニュー」で販売数の多い順に並べてみました。
＊このコードはチュートリアルにありません
[1] 「お楽しみメニュー」の日を抽出し、funというDataframeを作ります。
販売数が多い順に並べ替え、表示させます。

fun = train[train["remarks"]=="お楽しみメニュー"]
fun.sort_values('y', ascending=False)

なんと！上位８位までがカレーであることがわかりました。

箱ひげ図でも確認してみます。
[1] 「name」にカレーが含まれればx=1、それ以外はX=0を「curry」という新しい変数に入れます。

train["curry"] = train["name"].apply(lambda x : 1 if x.find("カレー")>=0 else 0)
sns.boxplot(x="curry",y="y",data=train)

やっぱりカレーは人気なんですね^^*

実装

以上の結果を踏まえて、改めてデータ取り込みから処理していきます。
チュートリアルでは以下のような方針が書かれています。
・日数が経過するにつれて減衰している為、売上数と日数の単回帰モデルを軸に検討する
　ただし、2014-05以前はやや傾向が異なる為、学習データから除く
・お楽しみメニューやカレー等は大きく寄与はしていそうだが、非線形な関係であることも考慮し
　Random Forestを用いて単回帰モデルの結果を修正するモデルも作成する。

前処理

それではデータ取り込み後、データ処理をしていきます。
[1] [2] trainとtestデータそれぞれに「t」という変数にラベルを付けます。
[3] trainとtestのデータを結合させます。sort=Trueはカラムの並べ替えです。

train["t"] = 1
test["t"] = 0
dat = pd.concat([train,test],sort=True).reset_index(drop=True)

trainとtestデータを結合するのは、「weather」にtestには無いものがtrainに含まれているため
別でダミー変換すると、変数の数に違いが生じるため、一緒にダミー変換をするためです。
後からtrainとtestを分割しやすいようにラベルを付けました。

[1] 「datetime」をdatetime型に変更し、それをindexに指定します。
[2] datの['2014-05-01']以降のデータに絞り、datへ再度入れなおします。
[3] indexをリセットします。

dat.index = pd.to_datetime(dat["datetime"])
dat = dat["2014-05-01":]
dat = dat.reset_index(drop=True)

[1] dat.indexは先程リセットされたので、1,2,3,,,という数です。
　これを「days」に入れます。これにより、先頭の日からの経過日数とすることができます。
[2] 「precipitation」の欠損値処理で、先程と同じです。
[3] 「remarks」のうち「お楽しみメニュー」の日は1、それ以外は0を「fun」という変数に入れます。
[4] 「curry」という変数の定義で、先程と同じです。
[5] colsという変数に該当のカラム名を入れます。

dat["days"] = dat.index
dat["precipitation"] = dat["precipitation"].apply(lambda x : -1 if x=="--" else x).astype(np.float)
dat["fun"] = dat["remarks"].apply(lambda x: 1 if x=="お楽しみメニュー" else 0)
dat["curry"] = dat["name"].apply(lambda x : 1 if x.find("カレー")>=0 else 0)
cols = ["precipitation","weather","days","fun","curry","y"]

学習

学習に必要なライブラリをインポートします。
[1] k-分割交差検証（trainとtestデータの分割）のライブラリ
[2] 評価関数である平均二乗誤差のライブラリ
[3] 線形回帰のライブラリ
[4] 回帰用のランダムフォレストのライブラリ

from sklearn.model_selection import KFold
from sklearn.metrics import mean_squared_error as MSE
from sklearn.linear_model import LinearRegression as LR
from sklearn.ensemble import RandomForestRegressor as RF

学習用関数の定義

学習を実行するための関数を定義します。
[1] learningという関数で、trainXとy_trainという引数を受けます。
　 trainXとy_trainは後で定義されます。
[2] 線形回帰をmodel1にセットします。
[3] ランダムフォレストをmodel2にセットします。
　 n_estimatorsは木の数、max_depthは決定木の深さの最大値です。
[4] model1を実行します。
　 trainXの「days」をreshapeしたものをxとして、y_trainをyとして単回帰分析します。
[5] [4]で学習したmodel1を使ってtrainXの予測値を算出し、predに入れます。
[6] 実測値のy_trainと予測値のpredの差をpred_subに入れます。
[7] model2を実行します。trainXのカラム名が「ｙ」でない列がxに、y_trainがyに指定されています。
[8] 最後にこの関数の返り値はmodel1とmodel2になります。

def learning(trainX,y_train):
    model1 = LR()
    model2 = RF(n_estimators=100,max_depth=4,random_state=777)
    model1.fit(trainX["days"].values.reshape(-1,1),y_train)
    pred = model1.predict(trainX["days"].values.reshape(-1,1))
    pred_sub = y_train - pred
    model2.fit(trainX.iloc[:, ~trainX.columns.str.match("y")],pred_sub)
    return model1, model2

trainデータで検証

trainデータのみでデータを分割し、その中でtrainとtestに分け、学習から評価までを行います。
[1] k-分割交差検証をkfにセットします。n_splitsは分割数を示します。
[2] ラベルが1(=train)であり、先程colsという変数に入れたカラム名のデータを抽出しtrに入れます。

kf = KFold(n_splits=5,random_state=777)
tr = dat[dat["t"]==1][cols]

[1] trainsとtestsという空のリストを準備します。
[2] kf.spit(tr)でtrの分割を実行し、train_indexとtest_indexにインデックスが入ります。
[3] [4] trainとtestのデータそれぞれに'tt'という変数を作り、ラベルを付けます。
[5] 'tt'をint型に変換します。
[6]trを一括でダミー変換しtmpに入れます。
[7] [9] trainXとtestXにそれぞれのラベルが付いたデータを入れます。
[8] [10] trainXとtestX内にある’tt’は不要なので削除します。
[11] [12] y_trainとy_testにそれぞれのラベルが付いている「y」の値を入れます。
[13] 先程定義した学習用関数にtrainXとy_tainを渡し、実行します。
　　実行した結果はそれぞれmodel1とmodel2に入ります。
[14] model1とmodel2にtrainXの値を入れ予測値を出します。
[15] [14]と同様にtestXを入れて推論を行います。
[16] 評価関数RMSEを求めるためにMSEに与えた推論値の0.5乗してルートを取ります。
[17] [18] 先程作った空のリストに計算したRMSEの値を追加していきます。
[19] AVGとは平均値を意味します。
[20] RMSEを入れたリストの平均を出し、それぞれ表示します。

trains ,tests = [] , []
for train_index, test_index in kf.split(tr):
    tr.loc[train_index,"tt"] = 1
    tr.loc[test_index,"tt"] = 0
    tr["tt"] = tr["tt"].astype(np.int)
    tmp = pd.get_dummies(tr)
    
    trainX = tmp[tmp["tt"]==1]
    del trainX["tt"]
    testX = tmp[tmp["tt"]==0]
    del testX["tt"]
    y_train = tmp[tmp["tt"]==1]["y"]
    y_test = tmp[tmp["tt"]==0]["y"]
    
    model1, model2 = learning(trainX, y_train)
    
    pred_train = model1.predict(trainX["days"].values.reshape(-1,1)) + model2.predict(trainX.iloc[:, ~trainX.columns.str.match("y")])
    pred_test = model1.predict(testX["days"].values.reshape(-1,1)) + model2.predict(testX.iloc[:, ~testX.columns.str.match("y")])
    
    print("TRAIN:",MSE(y_train,pred_train)**0.5, "VARIDATE",MSE(y_test, pred_test)**0.5)
    trains.append(MSE(y_train,pred_train)**0.5)
    tests.append(MSE(y_test, pred_test)**0.5)
print("AVG")
print(np.array(trains).mean(), np.array(tests).mean())

これでtrainデータに対する評価値を出すことが出来ました。

testデータで検証

[1] colsという変数にカラム名を入れます。
[2] ダミー変換をしますtmpに入れます。
[3] [5]ラベルが1はtrainXに、ラベルが0はtestXに入れます。
[4] [6] 不要になった't'を削除します。
[7] [8] それぞれのラベルが付いた「y」の値をy_trainとy_testに入れます。

cols = ["precipitation","weather","days","fun","curry","y","t"]
tmp = pd.get_dummies(dat[cols])
trainX = tmp[tmp["t"]==1]
del trainX["t"]
testX = tmp[tmp["t"]==0]
del testX["t"]
y_train = tmp[tmp["t"]==1]["y"]
y_test = tmp[tmp["t"]==0]["y"]

[1] データの準備ができたので、学習用関数で実行します。
[2] 学習した結果をモデルで推論を行います。
[3] 新たにpというDataFrameを作り、「actual」にy_train(実測値)を、「pred」にpred(予測値)を入れます。
[4] 上記のDataFrameをグラフに表示します。
[5] 評価関数RMSEの値を表示します。

model1, model2 = learning(trainX,y_train)
pred = model1.predict(trainX["days"].values.reshape(-1,1)) + model2.predict(trainX.iloc[:,~trainX.columns.str.match("y")])

p = pd.DataFrame({"actual":y_train,"pred":pred})
p.plot(figsize=(15,4))
print("RMSE",MSE(y_train,pred)**0.5)

山の部分がかなりフィットしています。
下に凸の部分は少し乖離が見られます。

同様にtrainデータで検証を行い、結果を下に表示します。

model1, model2 = learning(trainX,y_train)
pred = model1.predict(testX["days"].values.reshape(-1,1)) + model2.predict(testX.iloc[:,~testX.columns.str.match("y")])
plt.figure(figsize=(15,4))
plt.plot(pred)

以上でチュートリアルが終わります。

まとめ

チュートリアルでは上に凸するデータの学習が良くできていました。
課題として小さな波形状の予測が不十分な点が挙げられていたので、
ここを自分でも改善してみたいと思います。

他の人が作ったコードを読み解くのは、自分で作るのとは
また違った形で勉強になりなりました。
これからも、このような勉強方法も取り入れたいと思います。