勾配ブースティングで「幅を持たせた予測」を行う
背景と目的
ビジネスにおいて機械学習による回帰モデルを構築する際に、単に「予測値（点予測）」が知りたいというよりも、「予測値が（どれくらいの確率で）どれくらいの範囲に収まりそうか（区間予測）」を知りたい場面があります。例えば天気予報を眺めてみると台風の進路予測には予報円が出ており、ある程度幅を持たせて進路を予測するといったことが行われています。あるいは需要予測を行う場合に「95%の確率で40個~60個の間に入りそうだ」という幅を持たせた予測をすることができれば、「売れ残りがでると損失が大きいので、明日の仕入れ量は40個にしよう」といった意思決定に繋げることができるかもしれません。
そこで、本記事ではテーブルデータの回帰タスクで広く利用されており皆様にとって馴染み深い勾配ブースティングアルゴリズムを用いて、「幅を持たせた予測」を実装する方法を紹介します。
各手法の概要
この区間予測を実装する方法として筆者に思いつくのは以下の２つの手法です。
- NGBoost
- LightGBM
- Quantile Regression (分位点回帰)
- 事後的な予測誤差の推定
NGBoost
NGBoostとは勾配ブースティングによる区間予測を実施することを目的に開発されたライブラリです。公式ページは以下です。
実際、公開されている公式のユーザーガイドに従って分析を行うことで簡単に区間予測を実現することができます。
一方著者が使った感じでは大規模なデータになると LightGBM と比べて推論速度が大幅に悪化したり、そもそも特徴量の欠損値に対応していなかったりすることがあるため実用上の利用は少し怖いかなという印象でした。
LightGBM
言わずもがな非常に広く利用されているライブラリです。LightGBM を用いた幅を持たせた予測の実現方法として筆者に思いつくのは以下の２パターンです。
Quantile Regression (分位点回帰)
まず一つ目の手法として Quantile Regression(分位点回帰) というやり方が存在します。Quantile Regression では「予測値よりも目的変数が低い値をとる確率が
例えば
事後的な予測誤差のあてはめ
二つ目の手法として、「事後的に誤差の従う分布」を推定するというアプローチが考えられます。この手法は以下のページが参考になります。
これは「ひとまず LightGBM の通常の予測」を行うことで目的変数の平均的な値を近似できたと見なし、この近似的な目的変数の平均と実際の目的変数の値の誤差を確率分布でモデリングするアプローチです。
例えば学習データにおいて、「あるＡ店のおにぎりの販売個数の予測値と実際の値が大きくずれている一方で、Ｂ店のおにぎりの販売個数の予測値と実際の値があまりずれていない」という傾向が確認されたとしましょう。つまり「Ｂ店の販売個数は割と自信をもって予測できる一方でＡ店の販売個数の予測には自信がない」状態です。そこで、この誤差を正規分布を使って近似すると、Ａ店の正規分布の分散は高く推測されるがＢ店の正規分布の分散は低く推定されるといった形になります。
理論と実装
以下では、LightGBM を用いた２つの手法について簡単な理論と実装を確認します。実装は Google Colaboratory を用いて行いますので適宜そちらをご参照ください。それぞれに簡単な理論説明をつけていますが、必要に応じて読み飛ばしていただいて構いません。
なお、実装はこちらの記事内でも全て網羅していますが、github にも上げておきますので適宜ご参照ください。
まず２つの実装において共通して必要となるコードの準備をします。
必要なパッケージのインストール
!pip install japanize-matplotlib
必要なライブラリのインポート
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import lightgbm as lgbm
import seaborn as sns
ダミーデータの作成
ground truth (真の値) が
np.random.seed(1)
num_of_sample = 3000
x = np.random.uniform(low=-5, high=5, size=num_of_sample)
# 後のmatplotlibを用いた描画に備えて小さい順に並び替えておく。そうしないと描画がぐちゃぐちゃになる。
x = np.sort(x)
# ground truth (真の値) は sin カーブであるとする。
ground_truth = np.sin(x)
# 目的変数の実際の値は ground truth に誤差が乗った状態で現れる。
# 今回は正規分布に従って誤差が発生すると考えるが、その正規分布の分散は cos カーブの2乗で表現されるものとする。
# 2乗するのは分散にマイナスの値が表れるのを防ぐため。
true_var = np.cos(x)**2
error = np.random.randn(num_of_sample) * np.sqrt(true_var)
# 目的変数の値 = 真の値 + 誤差
y = ground_truth + error
# lightgbm に学習データとして渡すとき用に整形
x_for_lgbm = x.reshape((-1,1))
y_for_lgbm = y.reshape((-1,))
目的変数の値と ground truth を図示します。
fig, ax = plt.subplots(figsize=(12,6))
sns.scatterplot(
x=x, y=y, ax=ax, label='y'
)
sns.lineplot(
x=x, y=ground_truth, ax=ax, label='ground_truth', color='red'
)
ax.grid()
ax.legend()
LightGBMを用いた分位点回帰
理論
まず分位点回帰について確認します。分位点回帰は損失関数を Pinball loss 関数と呼ばれるものにすることにより実装します。
alphas = [0.1, 0.5, 0.9]
colors = ['blue', 'red', 'green']
labels = ['10%分位点', '50%分位点(中央値)', "90%分位点"]
def get_pinball_loss(y_diff: float, alpha: float):
if y_diff > 0:
return y_diff * alpha
else:
return -y_diff * (1 - alpha)
y_diffs = np.linspace(-5, 5, 100)
fig, ax = plt.subplots(figsize=(12,6))
for i, alpha in enumerate(alphas):
color = colors[i]
label = labels[i]
pinball_loss = [
get_pinball_loss(y_diff, alpha) for y_diff in y_diffs
]
sns.lineplot(
x=y_diffs,
y=pinball_loss,
label=label,
color=color,
ax=ax
)
ax.grid()
ax.set_xlabel('予測誤差（正の値であれば予測値>実際の値）')
fig.show()
図より、分位点の値によって上振れ、下振れのペナルティの重みが変わることが理解できます。例えば10%分位点（つまり小さめの予測）の時は、目的変数の値が予測値を上回った場合の損失はそれほど大きくありませんが、下回った場合の損失は大きくなっています。これにより、なるべく予測値を目的変数の値より小さくならないように定めるようになります。
以上の説明はあくまで直感的な説明に過ぎませんが、この損失関数を利用して予測値を導出することにより、ちょうど
累積型分布関数は「確率変数がある値以下になる確率」を示すものなので、
一応証明チックなものをつけておくので興味のある方はご参照ください。
証明チックなもの
pinball loss 関数の期待値は以下のように記述できます。
ここで、
最初の変形が難しい場合は 以下の解説が分かりやすいです。
したがって、
後は、
よって二階微分を計算すると
実装
さて以下では分位点予測を実装していきましょう。ここでは objective に "quantile" を指定し、alpha に分位点を与えるだけで実装できます。
lower_quantile = 0.025
upper_quantile = 0.975
predict_dict = {}
for quantile in [lower_quantile, upper_quantile]:
quantile_params = {
"objective": "quantile", # quantile を指定することで分位点予測となる。
"alpha": quantile # alpha には分位点の値を指定する。
}
quantile_regressor = lgbm.train(
quantile_params,
lgbm.Dataset(x_for_lgbm, y_for_lgbm),
)
quantile_predict = quantile_regressor.predict(x_for_lgbm)
predict_dict[quantile] = quantile_predict
可視化してみましょう。
fig, ax = plt.subplots(figsize=(12,6))
# ground truth
ax.plot(x, ground_truth, c='red', label='ground truth')
# 2.5% 分位予測
ax.plot(x, predict_dict[lower_quantile], c='green', label=f'{lower_quantile*100:.1f}%分位点予測')
# 97.5% 分位予測
ax.plot(x, predict_dict[upper_quantile], c='purple', label=f'{upper_quantile*100:.1f}%分位点予測')
# 描画で指定している alpha は透明度です。
ax.fill_between(
x,
predict_dict[lower_quantile],
predict_dict[upper_quantile],
color='orange',
alpha = 0.5,
label='90%予測区間'
)
sns.scatterplot(x=x,y=y,ax=ax, alpha=0.5)
ax.grid()
ax.legend()
fig.show()
とてもいい感じに推測できていることがうかがえます。それでは実際にこの予測区間内に収まっている正解データの割合を計算してみましょう。
upper = predict_dict[upper_quantile]
lower = predict_dict[lower_quantile]
total = 0
for i in range(num_of_sample):
upper_i = upper[i]
lower_i = lower[i]
y_i = y[i]
# 区間内に含まれたデータの数をカウント
total += (lower_i <= y_i <= upper_i)
# 割合を計算する。
print(total/num_of_sample)
LightGBM を用いた事後的な予測誤差の推定
理論
事後的な予測誤差の当てはめでは、まずは通常通り LightGBM による学習を行い目的変数
さてそのために、
これをパラメータ
ここで
さらに、
最後に、これを対数尤度を損失関数に置き換えることを考えます。対数尤度は「当てはまりの良さ」を表す指標なので、負の値を取ることで(推測値
さて、ここでは勾配の計算の都合上この
まず一階微分値は以下のように計算できます。
次に二階微分値は以下のように計算できます。
ここまでくれば準備完了です！実装に入りましょう！
実装
まず初めにすることは目的変数 y_pred) の導出です。これは以下のように実装できます。
# まず目的変数に対する予測値を導出する。
y_reg_params = {
"objective": 'regression'
}
y_regressor = lgbm.train(
y_reg_params,
lgbm.Dataset(x_for_lgbm, y_for_lgbm)
)
y_pred = y_regressor.predict(x_for_lgbm)
次に、この y_pred を用いて損失関数の対数分散についての一階微分と二階微分を導出する関数を定義し、これに基づいて標準偏差を推測しましょう。これは以下のように実装できます。
# 次に、分散推測モデルを構築し、分散の推測を行う。
def grad_hess_for_logvar(log_var: np.ndarray, data: lgbm.Dataset):
'''
推測された分散対数値に対する損失の一階微分及び二階微分を導出する
y_pred については引数としてではなくグローバル変数として与えていることに注意。
'''
y_true = data.get_label()
grad = 1/2 - 1/2 * (y_true - y_pred)**2 / np.exp(log_var)
hess = 1/2 * (y_true - y_pred)**2 / np.exp(log_var)
return grad, hess
log_var_reg_params = {
"objective": "regression"
}
log_var_regressor = lgbm.train(
log_var_reg_params,
lgbm.Dataset(x_for_lgbm, y_for_lgbm),
fobj=grad_hess_for_logvar # fobj 引数に自作の grad, hess 導出関数を定義する。
)
# 推測された標準偏差を導出する。
log_var_pred = log_var_regressor.predict(x_for_lgbm)
std_pred = np.sqrt(np.exp(log_var_pred))
この標準偏差を用いることで、95%信頼区間の上限と下限を計算することができます。
# 推測された分散から95%区間予測を導出する。
# 謎な場合には https://bellcurve.jp/statistics/course/8888.html などが分かりやすいです。
upper_pred = y_pred + 1.96 * std_pred
lower_pred = y_pred - 1.96 * std_pred
最後にこの結果を可視化してみましょう。
fig, ax = plt.subplots(figsize=(12,6))
ax.scatter(x, y, color='blue', s=1, alpha=0.5)
ax.plot(x, ground_truth, color='red', label='ground truth', linestyle='--')
ax.plot(x, y_pred, color='orange', label='pred')
ax.plot(x, upper_pred, color='purple', linestyle='--', label='97.5%分位点予測')
ax.plot(x, lower_pred, color='green', linestyle='--', label='2.5%分位点予測')
ax.grid()
ax.legend()
plt.show()
以上で実装は終了です。お疲れさまでした！
まとめ
本記事では勾配ブースティングを用いてどのように予測区間を導出できるかを概観しました。さらに LightGBM を用いた分位点予測と事後的な確率分布の当てはめについては実装まで含めて紹介いたしました。
分位点予測のメリットは、「目的変数の従う確率的な分布」に前提を置くことなく
一方で誤差に対する事後的な確率分布のあてはめについては先ほど紹介した分位点回帰とは異なり、どういった確率分布での誤差分布の近似を行うかといった前提を分析者が仮定する（今回の例でいえば正規分布）必要があり、それ自体が適切でないと分析が妥当なものになりません。また、予測値を「実際の平均値」と見なして誤差をモデリングする都合上、導出された予測値が平均値のよい近似でなければ誤った分析となってしまいます。例えば過学習をしていた場合誤差を低く見積もってしまうことになります。しかしながら、分位点回帰とは異なり誤差の従う確率分布を推定するため、いちいち必要な分位点についてモデルをたくさん作る必要がない点がメリットになります。
以上になります。何か皆様の参考になるところがあれば幸いです。長文最後まで読んでいただきありがとうございました！
Discussion
途中で紹介頂いた「LightGBMで予測区間を推定する」記事の筆者です。
分位点回帰の発想は無かったので非常に勉強になりました！ありがとうございます！
Yuhsak Inoue さま
コメントいただきありがとうございます！まさか記事を執筆くださった方にコメントいただけると思っていてなくて、本当に嬉しいです。
こちらこそ、当時「予測区間の推定ってどうやるんだろう...？」と色々調べていた時にとても勉強になって記事として執筆させていただきました。有益な情報発信を本当に有難うございます！