Bike Sharing Data Setを使って時系列データの予測をやってみた

はじめに

時系列データの回帰分析をする機会があったので、その練習がてら試したことをまとめました。

実装

利用するデータ

今回は以下のBike Sharing Data Setを利用します。
archive.ics.uci.edu

カラム名がわかりやすいので、あまり解説はいらないと思いますが、以下のカラムを持っています。
(hour.csv, day.csvの2つがありますが、hour.csvの方を利用します)

#	カラム名	大まかな意味
1	instant	index
2	dteday	そのデータを観測した年月日時
3	season	季節。1が冬、2が春、3が夏、4が秋
4	yr	観測した年。2011年が0, 2012年が1
5	mnth	観測した月
6	hr	観測した時刻
7	holiday	休みなら1, 平日なら0
8	weekday	月曜〜日曜までを0~7で表現
9	workingday	holidayでも土日でもなければ1
10	weathersit	天気。1は概ね晴れ、2は曇り系、3は雪、4は強い雨
11	temp	標準化した気温
12	atemp	標準化した体感気温
13	hum	湿度
14	windspeed	風速
15	casual	レンタルサイクルの観光者利用数
16	registered	レンタルサイクルの登録利用者数
17	cnt	トーナルのレンタルサイクルの利用者数

24時間 * 2年分のデータになります。

ライブラリのインポート

今回はこの辺使います。ここについては特に説明はしません。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datetime import date
from datetime import datetime as dt
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
import seaborn as sns
import statsmodels.api as sm

欠損値の処理

理由はわからないですが、よく見ると時間単位でデータの欠損がありますので、まずはそれの処理をします。
どこが抜けているのかよくわからないので、2011/1/1 ~ 2012/12/31までの1時間単位のdatetimeを生成して、outer join相当のことをします。

# データの読み込み
data = pd.read_csv("hour.csv")

# dteday(文字列)をdatetimeオブジェクトに変換
def to_datetime(x):
    sp = x.split("-")
    return datetime.datetime(year=int(sp[0]), month=int(sp[1]), day=int(sp[2]))
data["date"] = data[["dteday"]].applymap(to_datetime)

# 1時間単位で時刻データの作成
delta = datetime.timedelta(hours=1)
curr = datetime.datetime(2011, 1, 1, hour=0)
time_series = []
for i in range(24 * 365 * 2 + 24):
    time_series.append([curr, datetime.datetime(curr.year, curr.month, curr.day), curr.hour])
    curr += delta
time_df = pd.DataFrame(time_series, columns=["datetime", "date", "hr"])

# outer join
df_with_nan = pd.merge(time_df, data, on=["date", "hr"],how="outer")

これでできました。あと、今後不要なカラムはここでドロップしてしまいます。

df_with_nan = df_with_nan.drop(columns=["instant", "dteday"])

欠損値の処理

上記に伴って欠損しているフィールドがいくつもあります。
今回は適当に前後の値の中央値を取ることでfillします。
特にweathersitはこの方法では良くないのですが、一旦簡単に進めるために全部これで進めます。

rows = df_with_nan.values
for i, row in enumerate(rows):
    if np.isnan(row[3]):
        n = i+1
        while np.isnan(rows[n][3]):
            n += 1
        for j in range(2, df_with_nan.shape[1]):
            row[j] = (rows[i-1][j] + rows[n][j]) / 2
            if j == 9: # weathersit
                row[j] = np.round(row[j])
data = pd.DataFrame(rows, columns=df_with_nan.columns)
# 上記の方法だと、すべてのカラムがobject型になるので適切な型に戻す
for c in data.columns:
    data[c] = data[[c]].applymap(lambda x: x)

特徴量エンジニアリング①

ここから先の処理のために、週番号を追加します。

data["week_num"] = data[["datetime"]].applymap(lambda x:x.isocalendar()[1])

また、よく見ると、hr, weekdayはきれいな数字ではなく小数点以下の値を持っているため、丸めてしまいます。

data["hr"] = data[["hr"]].applymap(lambda x: np.round(x))
data["weekday"] = data[["weekday"]].applymap(lambda x: np.round(x))

可視化

準備が整ったので、可視化していきます。

まずは週単位で利用者数を合計して、利用者数の推移を見ていきます。

fig,ax = plt.subplots(figsize=(15,10))
counts = data.groupby(["yr", "week_num"])[["cnt", "casual", "registered"]].sum()
counts.plot(ax=ax)
plt.grid()

なんとなく、周期性がある（夏多くて冬少ない）のと、全体的には上昇傾向があることが見えると思います。

ちょっと深堀りして、季節ごとの利用者数を箱ひげで書いてみました。

sns.boxplot(x="season",y="cnt",data=data)

冬が少ないのは上で見たとおりでした。

ということは、温度と利用者数も関係ありそうな気がしてきますので、温度(atemp)と利用者数を可視化してみます。

fig,ax = plt.subplots(figsize=(15,10))
sns.scatterplot(x="atemp", y="cnt", data=data, ax=ax)

主観的ですが、基本的に温度が上がるほど利用数は増加しますが、0.7~0.8を境に利用者数は減少傾向にあるように見えます。

関係ありそうといえば、時間帯によっても差がありそうですね。
1日の中の使用率という尺度で可視化しました（そうしないと絶対値の影響を受けてしまうため）

fig,ax = plt.subplots(figsize=(15,10))
ratio_by_hour = data.loc[:]
tmp = pd.merge(ratio_by_hour, data.groupby(["date"])[["cnt"]].sum().reset_index(), on=["date"])
tmp["cnt_ratio"] = tmp["cnt_x"] / tmp["cnt_y"]
tmp["cnt_ratio"] = tmp["cnt_ratio"].astype(float)
sns.boxplot(x="hr", y="cnt_ratio", data=tmp, ax=ax)

通勤時間帯と日中帯、その他で分別できるように思います。

関係ありそうといえば、曜日も関係ありそうです。
可視化してみましょう。

sns.boxplot(x="weekday",y="cnt",data=data)

あれ、ほとんど変わらないですね。案外曜日というのは影響しないようです。

あとは天気ですかね。見てみましょう。

やはり天気が悪くなるほど利用者数は減少傾向ですね。

特徴量エンジニアリング②

先程可視化で見つけた、

• 日中(9~16時)は比較的利用率が高く、通勤時間(8, 17, 18時)はとても利用率が高い。
• atempが0.7 ~ 0.8を超える辺りから利用率は減少傾向
• 天気、季節の影響があるが、それぞれ連続値で扱うのは不適に見えるため、ダミー変数化する

という3点を反映するような特徴量を生成します。

data["daytime"] = ((9 <= data["hr"]) & (data["hr"] <= 16)) * 1
data["transittime"] = ((8 == data["hr"]) | (data["hr"] == 17) | (data["hr"] == 18)) * 1
data["overheat"] = (0.8 < data["atemp"]) * 1
data[["spring", "summer","autom"]] = pd.get_dummies(data["season"],drop_first=True)
data[["weather2", "weather3","weather4"]] = pd.get_dummies(data["weathersit"],drop_first=True)

そして、学習データは以下のようにまとめます。

X_train = data.drop(columns=["datetime", "date","cnt", "casual", "registered", "yr", "hr", "week_num", "weekday", "temp", "season"])
Y_train = data[["cnt", "casual", "registered"]]

今回はcntを目的変数にしますが、casualやregisteredも試せるように3つ持たせておきます。
結果的に、説明変数は以下のような形になります。

線形回帰でモデリング

とりあえず何も考えずに線形回帰器にいれて、実際の値と比べてみます。

from sklearn.linear_model import LinearRegression
lr = LinearRegression()
lr.fit(X_train, Y_train["cnt"])
preds = lr.predict(X_train)
preds_df = pd.concat([data, pd.DataFrame(preds[:, np.newaxis], columns=["preds"])], axis=1)

# 可視化
fig,ax = plt.subplots(figsize=(15,10))
counts = preds_df.groupby(["yr", "week_num"])[["cnt", "preds"]].sum()
counts.plot(ax=ax)
plt.grid()

なんだか微妙な感じですね。
個人的には全体にある上昇傾向が掴めてない点が得に気にかかりました。

また、実際の予測結果（時間単位）を見てみると、

sorted(preds)[:15]

なんかマイナスの予測値を出してしまっています。
今回利用台数の予測をするのに、マイナスの値が出てしまうのはなんとも微妙というか、ありえない予測モデルだなぁと思うわけです。

上記2つの課題を1つずつ潰していきましょう。

特徴量エンジニアリング3: 上昇傾向への対応

上昇傾向があるときはラグ特徴量が効くと聞いたので実際に試してみます。
今回は2011/1/1からの経過日数を表すフィールドを追加します。

base = datetime.datetime(2011, 1, 1, hour=0)
data["lag"] = data[["date"]].applymap(lambda x: (x-base).total_seconds() / (24*60*60))

また、lagはスケールが大きいので標準化しておきます。
(実は他のパラメータは0~1の範囲に既に収まっています)

ss = StandardScaler()
data["lag_ss"] = ss.fit_transform(data[["lag"]])
X_train = data.drop(columns=["datetime", "date","cnt", "casual", "registered", "yr", "hr", "week_num", "weekday", "temp", "season", "lag"])
Y_train = data[["cnt", "casual", "registered"]]

では、改めてモデル化してみましょう。結果がマイナスになることは一旦無視して線形回帰で試します。

lr_lag = LinearRegression()
lr_lag.fit(X_train, Y_train["cnt"])

preds = lr_lag.predict(X_train)
preds_df = pd.concat([data, pd.DataFrame(preds[:, np.newaxis], columns=["preds"])], axis=1)

fig,ax = plt.subplots(figsize=(15,10))
counts = preds_df.groupby(["yr", "week_num"])[["cnt", "preds"]].sum()
counts.plot(ax=ax)
plt.grid()

おっ！今度は上昇傾向を捉えられていますね！！

でも相変わらず、結果にはマイナスの値が含まれています。

ポアソン回帰でモデリング

今回予測したい利用される自転車の台数は正の離散値です。
こういう分布はきっとポアソン分布なんだろうと思って、ポアソン回帰を利用してみることにします。
調べた限り、scikit-learnでポアソン回帰は使えなかったので、statsmodelというパッケージを利用します。

# statsmodelで回帰をするときは定数項を足す
X_train_po = sm.add_constant(X_train)
glm = sm.GLM(Y_train["cnt"], X_train_po, family=sm.families.Poisson())
model = glm.fit()
preds_po = model.predict(X_train_po)
preds_df = pd.concat([data, pd.DataFrame(preds_po[:, np.newaxis], columns=["preds"])], axis=1)
fig,ax = plt.subplots(figsize=(15,10))
counts = preds_df.groupby(["yr", "week_num"])[["cnt", "preds"]].sum()
counts.plot(ax=ax)
plt.grid()

おお！いい感じに予測できてますね。
では、実際の値はというと・・・

おお！ちゃんと正の値になってますね！これなら使えそうです。

まとめ

今回はBike Sharing Data Setを使って時系列データの回帰予測をしてみました。
なんとなく、一般化線形モデルについて少しは理解が進んだようにも思うし、時系列データの取り扱いが少しわかった気もします。

ただ、今回は以下については触れていません

1. 予測誤差の評価
2. テストデータセットに対する検証

1.については今回はどういったモデルを作っていけばよいのか、おおまかなところをやったのみで、特定のモデルを徹底して精度向上をさせていくというものではないので、面倒になるのであえてやっていませんが、本来はもちろんやるべきことです。
2.については今回はデータがで2周期分しかなく、どこで分割すればいいんじゃい、となったのでやりませんでした。もしやり方を知っていたり、アイディアのある方はぜひ教えてほしいです。

というわけで、今回はこれでおしまいです。