近年、機械学習は様々な分野で重要な役割を果たしており、その応用範囲はますます広がっています。特に、回帰問題において、従来の機械学習モデルは、与えられた特徴量に基づいて単一の「最も可能性の高い」予測値、つまり点推定を返すことが一般的でした。しかし、実際の多くの問題では、予測値だけでなく、その不確実性を評価することが不可欠です。

例えば、医療分野での患者の生存予測や、気象予測などの分野では、予測結果に対する信頼度や、ある範囲内に結果が収まる確率を知ることが非常に重要になります。従来の点推定では、これらの要求に応えることは困難でした。そこで、確率分布推定というアプローチが重要になります。これは、単一の予測値ではなく、結果の確率分布全体を推定することで、予測の不確実性を定量的に評価することを可能にするものです。

今回は、この確率分布推定を効率的に行うためのブースティングアルゴリズムであるNGBoost（Natural Gradient Boosting）について解説します。NGBoostは、柔軟性、拡張性、使いやすさを兼ね備えており、様々な分野での応用が期待されています。具体的には、NGBoostがどのようにして予測の不確実性を捉え、実用的な問題解決に貢献するのか、そのメカニズムと利点について掘り下げていきます。

従来の機械学習手法の欠点

従来の機械学習モデルは、回帰問題において、点推定に焦点を当ててきました。これは、モデルが与えられた入力（特徴量）に対して、単一の「最も可能性の高い」予測値を返すというものです。例えば、明日の気温を予測する場合、モデルは「明日の気温は16℃になる」といった単一の数値を予測します。

しかし、このアプローチにはいくつかの重要な限界があります。

• 不確実性の無視: 点推定は、予測の不確実性を考慮していません。実際には、予測には必ず誤差が伴い、その誤差の範囲や分布を知ることが重要となる場面が多くあります。例えば、医療分野での生存予測では、単に「生存期間はX日」という予測だけでなく、「X日以上生存する確率はY%」といった、確率的な情報が求められます。また、気象予測においても、単一の気温予測だけでなく、「18℃から20℃になる確率はZ%」といった、区間推定が役立つことがあります。
• 限定的な情報: 点推定は、予測に関する情報を単一の値に限定してしまいます。これにより、予測の背後にある様々な可能性や、予測に対する信頼度を把握することが困難になります。
• 実用上の課題: 多くの実務的な問題では、単一の予測値だけでは不十分です。例えば、意思決定を行う際には、起こりうる様々なシナリオを考慮し、それぞれのシナリオに対する確率を評価する必要があります。点推定では、このような多角的な視点を提供することができません。

上記の問題を解決するためには、確率分布推定というアプローチが重要になります。確率分布推定では、単一の予測値ではなく、結果の確率分布全体を推定します。これにより、予測の不確実性を定量的に評価することが可能になり、より信頼性の高い予測や、より実用的な意思決定を支援することができます。

従来の機械学習モデルにおける点推定は、特に以下のような状況において問題となります。

• 医療分野: 患者の生存時間予測、疾患リスク評価など、結果の不確実性を考慮する必要がある場面
• 気象予測: 気温、降水量などの予測、異常気象の発生確率など、予測の幅を考慮する必要がある場面
• 金融分野: 株価予測、リスク評価など、市場の変動性を考慮する必要がある場面

これらの分野では、単に「最も可能性の高い」予測値だけでなく、「どの程度の確率でその予測が起こりうるのか」を知ることが非常に重要です。点推定は、このような確率的な情報を提供することができないため、これらの分野では不十分であると言えます。

NGBoostとは？

NGBoost（Natural Gradient Boosting）は、確率分布推定のための新しいブースティングアルゴリズムです。従来の勾配ブースティングを拡張し、予測値だけでなく、予測の不確実性を表現する確率分布全体を効率的に推定することを目的としています。NGBoostは、柔軟性、拡張性、使いやすさを兼ね備え、さまざまな分野への応用が可能です。

NGBoostの基本的な考え方

NGBoostは、以下の3つの主要な要素を組み合わせることで、確率分布推定を実現します。

• ベース学習器（Base Learner）：予測モデルの基本となる学習器です。例えば、決定木（Regression Tree）などが用いられます。NGBoostは、任意のベース学習器を使用できるモジュール性を備えています。
• パラメトリック確率分布（Parametric Probability Distribution）：予測対象となる確率分布の形式を定めます。例えば、正規分布（Normal）、ラプラス分布（Laplace）など、さまざまな分布を選択できます。NGBoostは、連続パラメータを持つ任意の分布を扱うことができます。
• スコアリングルール（Scoring Rule）：推定された確率分布と実際の観測値とのずれを評価するための指標です。例えば、最尤推定（MLE）やCRPS（Continuous Ranked Probability Score）などが用いられます。NGBoostは、任意のスコアリングルールを使用できます。

NGBoostは、これらの要素を組み合わせ、勾配ブースティングの枠組みの中で、確率分布のパラメータを推定します。具体的には、以下のステップで学習が進みます。

自然勾配（Natural Gradients）の重要性

NGBoostの重要な特徴の一つは、自然勾配（Natural Gradients）を使用している点です。従来の勾配ブースティングでは、通常の勾配が用いられますが、これはパラメータの再パラメータ化に対して不変ではありません。つまり、パラメータの表現方法を変えると、学習のダイナミクスが大きく変わってしまう可能性があります。

一方、自然勾配は、パラメータ空間における距離が、確率分布空間における距離に対応するように調整された勾配であるため、再パラメータ化に対して不変です。これにより、パラメータの表現方法に依存せず、より安定した学習が可能になります。また、自然勾配を用いることで、各パラメータが「最適に事前スケーリング」された状態で学習されるため、学習の収束も早くなります。

マルチ・パラメータ・ブースティング

NGBoostは、マルチ・パラメータ・ブースティングというアプローチを採用しています。これは、確率分布の複数のパラメータを同時に学習する手法です。例えば、正規分布の場合、平均（\(\mu\)）と標準偏差（\(\sigma\)）の2つのパラメータを同時に学習します。従来の勾配ブースティングでは、通常、平均のみを学習しますが、NGBoostでは、分布の形状を決定する複数のパラメータを同時に学習することで、より柔軟な確率分布推定を可能にします。

スコアリングルール

NGBoostでは、任意の適切なスコアリングルールを使用できます。スコアリングルールとは、予測された確率分布と実際の観測値とのずれを評価するための指標です。NGBoostでは、最尤推定（MLE）やCRPS（Continuous Ranked Probability Score）などのさまざまなスコアリングルールを使用できます。

• 最尤推定（MLE）：観測されたデータに対して、最も尤もらしい確率分布のパラメータを推定する手法です。
• CRPS：予測分布と実際の観測値とのずれを測るスコアリングルールです。MLEよりも外れ値にロバストであるとされています。

NGBoostの利点

NGBoostは、以下の点で既存の手法よりも優れています。

• 柔軟性：任意のベース学習器、確率分布、スコアリングルールを使用できる。
• 拡張性：さまざまな問題（回帰、分類、生存分析など）に適用できる。
• 使いやすさ：特別な知識や調整を必要とせず、すぐに使用できる。
• 高い性能：既存の確率分布推定手法と同等またはそれ以上の性能を発揮する。
• 自然勾配による安定した学習：パラメータの再パラメータ化の影響を受けにくく、学習が安定する。
• マルチ・パラメータ・ブースティングによる柔軟な分布推定：確率分布の複数のパラメータを同時に学習できるため、予測の不確実性をより正確に表現できる。

NGBoostの計算コスト

NGBoostの計算コストは、以下の2つの要因によって増加します。

1. 各パラメータに対する学習器の系列：NGBoostでは、確率分布のパラメータごとに学習器の系列を適合させる必要があるため、計算コストはパラメータ数（\(p\)）に比例して増加します。
2. 自然勾配の計算：各データ点に対して、自然勾配を計算する必要があり、\(p \times p\)の行列の逆行列を計算する必要があるため、計算コストは\(p^3\)に比例して増加します。

しかし、実際には、一般的に使用される分布のパラメータ数は1つか2つであり、5つ以上のパラメータを持つ分布は非常にまれであるため、\(p\)に関する計算コストは大きな問題にはなりません。また、データ数が大きい場合は、ミニバッチをサブサンプリングすることで計算コストを削減することも可能です。

PythonによるNGBoostの実装

ここからは、Pythonを用いてNGBoostを実装する方法について解説していきます。NGBoostのPython実装は、オープンソースとしてPyPIにパッケージが公開されています。まずは、pipコマンドを使用して必要なパッケージをインストールします。

$ pip install ngboost

# conda を利用する場合はこちら
$ conda install -c conda-forge ngboost

では、ここからコードを実装していきます。まずは、必要なパッケージを読み込みます。

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from ngboost import NGBRegressor
from ngboost.distns import Normal
from ngboost.scores import LogScore

from scipy.stats import norm

from sklearn.datasets import fetch_california_housing
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.model_selection import train_test_split
from sklearn.tree import DecisionTreeRegressor

次に使用するデータセットですが、今回はscikit-learnに同梱されているCalifornia Housing Dataset（カリフォルニア住宅価格データセット）を使用します。データを読み込んだら、次のようにして説明変数を確認します。説明変数は8つで、すべて数値データであることが分かります。

dataset = fetch_california_housing(as_frame=True)

dataset["data"].head()

   MedInc  HouseAge  AveRooms  AveBedrms  Population  AveOccup  Latitude  Longitude
0  8.3252      41.0  6.984127   1.023810       322.0  2.555556     37.88    -122.23
1  8.3014      21.0  6.238137   0.971880      2401.0  2.109842     37.86    -122.22
2  7.2574      52.0  8.288136   1.073446       496.0  2.802260     37.85    -122.24
3  5.6431      52.0  5.817352   1.073059       558.0  2.547945     37.85    -122.25
4  3.8462      52.0  6.281853   1.081081       565.0  2.181467     37.85    -122.25

合わせて、ヒストグラムを表示して目的変数についても確認しておきます。グラフから分かるように、目的変数は正規分布から外れていることが分かります。

skew = dataset['target'].skew()  # 歪度
kurt = dataset['target'].kurt()  # 尖度

dataset["target"].plot(kind="hist", bins=30, alpha=0.5, color="b")
dataset["target"].plot(kind="kde", secondary_y=True, color="skyblue")
plt.title(f"Probability distribution (Skewness: {skew:.03f}, Kurtosis: {kurt:.03f})")
plt.xlim([0.0, 5.2])
plt.show()

では、データセットを訓練データとテストデータに分割しましょう。

今回は、決定木を弱識別器としてNGBoostを学習させます。また、説明変数がすべて数値データであるため、データの前処理は省略します。

X_train, X_test, Y_train, Y_test = train_test_split(dataset["data"], dataset["target"], test_size=0.1, shuffle=True)

print(f"訓練データ: {len(X_train)}サンプル, テストデータ: {len(X_test)}サンプル")

訓練データ: 18576サンプル, テストデータ: 2064サンプル

次に、NGBoostの学習を行います。弱識別器には決定木を使用します。また、目的変数の分布は正規分布から外れていますが、ここでは正規分布であると仮定して学習を行います。具体的な実装は次の通りで、scikit-learnのような使い勝手でモデルを学習させることができます。なお、NGBoostのハイパーパラメータについては、コード中のコメントを参照してください。

base_learner = DecisionTreeRegressor(
    criterion="squared_error",
    min_samples_split=2,
    min_samples_leaf=1,
    min_weight_fraction_leaf=0.0,
    max_depth=4,
    splitter="best",
    random_state=None,
)

ngb = NGBRegressor(
    Dist=Normal,                   # 目的変数の確率分布
    Score=LogScore,                # 予測と観測データの比較指標
    Base=base_learner,             # 弱識別器

    n_estimators=2000,             # ブースティング反復回数
    learning_rate=1.0E-2,          # 学習率
    minibatch_frac=1.0,            # 各反復で使用するデータの割合（行）
    col_sample=1.0,                # 各反復で使用するデータの割合（列）
    verbose=True,                  # 学習中のログ出力
    verbose_eval=100,              # ログ出力頻度
    tol=1.0E-4,                    # 最適化の許容誤差
    random_state=None,             # 再現性維持のための乱数シード
    validation_fraction=0.1,       # 早期停止用検証データの割合
    early_stopping_rounds=None,    # 早期停止までの最大反復回数（損失増加時）
)

ngb.fit(X_train, Y_train)

[iter 0] loss=1.5640 val_loss=0.0000 scale=1.0000 norm=1.1068
[iter 100] loss=1.0726 val_loss=0.0000 scale=2.0000 norm=1.4747
[iter 200] loss=0.7715 val_loss=0.0000 scale=2.0000 norm=1.3191
[iter 300] loss=0.6241 val_loss=0.0000 scale=2.0000 norm=1.2962
[iter 400] loss=0.5587 val_loss=0.0000 scale=1.0000 norm=0.6522
[iter 500] loss=0.5168 val_loss=0.0000 scale=1.0000 norm=0.6533
[iter 600] loss=0.4837 val_loss=0.0000 scale=1.0000 norm=0.6505
[iter 700] loss=0.4521 val_loss=0.0000 scale=1.0000 norm=0.6435
[iter 800] loss=0.4281 val_loss=0.0000 scale=1.0000 norm=0.6398
[iter 900] loss=0.4077 val_loss=0.0000 scale=1.0000 norm=0.6360
[iter 1000] loss=0.3874 val_loss=0.0000 scale=1.0000 norm=0.6316
[iter 1100] loss=0.3672 val_loss=0.0000 scale=1.0000 norm=0.6265
[iter 1200] loss=0.3504 val_loss=0.0000 scale=1.0000 norm=0.6231
[iter 1300] loss=0.3335 val_loss=0.0000 scale=1.0000 norm=0.6194
[iter 1400] loss=0.3218 val_loss=0.0000 scale=1.0000 norm=0.6182
[iter 1500] loss=0.3094 val_loss=0.0000 scale=0.5000 norm=0.3077
[iter 1600] loss=0.2978 val_loss=0.0000 scale=1.0000 norm=0.6130
[iter 1700] loss=0.2873 val_loss=0.0000 scale=0.5000 norm=0.3057
[iter 1800] loss=0.2778 val_loss=0.0000 scale=1.0000 norm=0.6095
[iter 1900] loss=0.2684 val_loss=0.0000 scale=1.0000 norm=0.6077

学習が完了したら、テストデータを用いてモデルの性能を評価します。以下では、平均二乗誤差とR²スコアを用いて評価を行います。

# 平均二乗誤差
mse_train = mean_squared_error(ngb.predict(X_train), Y_train)
mse_test = mean_squared_error(ngb.predict(X_test), Y_test)
print(f"MSE: TrainScore={mse_train}, TestScore={mse_test}")

# R2スコア
r2_train = r2_score(ngb.predict(X_train), Y_train)
r2_test = r2_score(ngb.predict(X_test), Y_test)
print(f"R2 Score: TrainScore={r2_train}, TestScore={r2_test}")

MSE: TrainScore=0.1725062719521442, TestScore=0.2161245167553553
R2 Score: TrainScore=0.8428209847323425, TestScore=0.7881002381084635

NGBoostの利点は、確率分布を予測できることにあります。その方法について紹介します。pred_distメソッドを利用することで、次のように目的変数の確率分布を予測できます。今回は、目的変数に正規分布を仮定しているため、予測される確率分布のパラメータは平均（loc）と標準偏差（scale）になります。

# 確率分布の予測
Y_dists = ngb.pred_dist(X_test)
print(Y_dists[0].param)

{'loc': np.float64(1.6180486741477984), 'scale': np.float64(0.24827396011543082)}

また、次のようにして、予測した確率分布から負の対数尤度（Negative Log Likelihood）を計算することもできます。

# 負の対数尤度 (Negative Log Likelihood)
nll_test = -Y_dists.logpdf(Y_test).mean()

print(f"NLL={nll_test}")

NLL=0.5540062261492721

それでは、実際に予測した確率分布と正解値をグラフに表示し、予測された確率分布がどの程度妥当であるかを確認しましょう。次のように、グラフ中の曲線が予測された確率分布（正規分布）で、点線が正解値（ラベルデータ）に対応します。確認すると分かりますが、予測された確率分布の\(\pm 1 \sigma\)範囲内に正解値があることから、予測された確率分布はおおむね妥当であると考えられます。

target_index = 0

loc = Y_dists[target_index].params["loc"]
scale = Y_dists[target_index].params["scale"]

x = np.arange(loc - 5.0 * scale, loc + 5.0 * scale, 0.02)
y = norm.pdf(x, loc, scale)

plt.plot(x, y, color="m", label="Predicted distribution")
plt.vlines(Y_test.iloc[target_index], 0.0, np.max(y), color="b", linestyles="dotted", label="Correct")
plt.title(f"Probability distribution (loc: {loc:.03f}, scale: {scale:.03f})")
plt.xlabel(dataset["target_names"][0])
plt.ylabel("Probability Density")
plt.legend(loc="best")
plt.show()

最後に、予測した確率分布のパラメータ（今回は平均と標準偏差）について、説明変数の重要度（Feature Importance）を確認します。

feature_importance_loc = ngb.feature_importances_[0]
feature_importance_scale = ngb.feature_importances_[1]

df_loc = pd.DataFrame({"feature": dataset["feature_names"], "importance": feature_importance_loc})
df_loc = df_loc.sort_values("importance", ascending=True)

df_scale = pd.DataFrame({"feature": dataset["feature_names"], "importance": feature_importance_scale})
df_scale = df_scale.sort_values('importance', ascending=True)

fig, ax = plt.subplots(1, 2, figsize = (10, 5), tight_layout=True)
fig.suptitle("Feature importance plot for distribution parameters", fontsize=17)

ax[0].barh(df_loc["feature"], df_loc["importance"], color="b", height=0.5)
ax[0].set_title("loc")
ax[0].set_xlabel("Importance")
ax[0].set_ylabel("Feature Names")

ax[1].barh(df_scale["feature"], df_scale["importance"], color="b", height=0.5)
ax[1].set_title("scale")
ax[1].set_xlabel("Importance")
ax[1].set_ylabel("Feature Names")

plt.tight_layout()
plt.show()

以上のように、NGBoostはscikit-learnのように簡単に利用できます。今回は回帰問題を扱いましたが、NGBoostには分類問題に対応したモデルも用意されています。また、目的変数の確率分布についても正規分布を指定しましたが、t分布や指数分布など一般的な確率分布であれば利用可能です。今回ご紹介しきれなかった内容については、GitHubのリポジトリやユーザーガイドなどをご参照ください。

おわりに

NGBoostは、最先端の確率的予測を実現するアルゴリズムで、マルチ・パラメータ・ブースティングと自然勾配を組み合わせることで、予測分布のパラメータを効率的に推定します。 この手法は、柔軟性（分類、回帰、生存分析に対応）、拡張性（大規模データに対応）、使いやすさ（専門知識不要）を兼ね備えています。

元論文では、実験結果をもとに以下の点を主張しています。

• マルチ・パラメータ・ブースティングと自然勾配の組み合わせが重要
• 均一分散の仮定は必ずしも最適ではない
• 自然勾配が学習ダイナミクスを改善
• 点推定でも優れた性能を発揮

また、今後の展望として、以下の点が挙げられます。

• 分類問題や生存分析への応用
• 共同予測（複数の結果を同時予測）
• 高次の不変性を持つ自然勾配
• より優れたベース学習器の利用
• 推論問題への応用

NGBoostは、確率的予測において優れた性能を発揮しますが、今後の研究により、その理論的側面や応用範囲がさらに拡大することが期待されます。

More Information

• arXiv:1910.03225, Tony Duan et al., 「NGBoost: Natural Gradient Boosting for Probabilistic Prediction」, https://arxiv.org/abs/1910.03225