LabEx
ログイン登録

線形ベイズ回帰の比較

Beginner

This tutorial is from open-source community. Access the source code

はじめに

この実験では、合成データセットを使って、2 つの異なるベイズ回帰モデル：自動相関決定（Automatic Relevance Determination：ARD）とベイズリッジ回帰（Bayesian Ridge Regression）を比較します。実験の最初の部分では、通常最小二乗法（Ordinary Least Squares：OLS）モデルをベースラインとして、モデルの係数を真の係数と比較します。最後のセクションでは、多項式特徴量拡張を使ってXとyの間の非線形関係をフィットさせ、ARD とベイズリッジ回帰の予測と不確定性をプロットします。

VM のヒント

VM の起動が完了したら、左上隅をクリックしてノートブックタブに切り替え、Jupyter Notebook を使って練習しましょう。

時々、Jupyter Notebook が読み込み終了するまで数秒待つ必要があります。Jupyter Notebook の制限により、操作の検証は自動化できません。

学習中に問題に遭遇した場合は、Labby にお問い合わせください。セッション後にフィードバックを提供してください。すぐに問題を解決いたします。

合成データセットの生成

Xとyが線形に関連付けられた合成データセットを生成します。Xの 10 個の特徴量を使ってyを生成します。他の特徴量はyの予測には役立ちません。また、n_samples == n_featuresのデータセットも生成します。このような設定は OLS モデルにとってチャレンジングであり、潜在的に任意の大きな重みをもたらす可能性があります。重みに事前分布を持ち、ペナルティを課すことで問題を軽減します。最後に、ガウスノイズを追加します。

from sklearn.datasets import make_regression

X, y, true_weights = make_regression(
    n_samples=100,
    n_features=100,
    n_informative=10,
    noise=8,
    coef=True,
    random_state=42,
)

回帰モデルのフィッティング

後でモデルの係数を比較するために、2 つのベイズモデルと OLS をフィッティングさせます。

import pandas as pd
from sklearn.linear_model import ARDRegression, LinearRegression, BayesianRidge

olr = LinearRegression().fit(X, y)
brr = BayesianRidge(compute_score=True, n_iter=30).fit(X, y)
ard = ARDRegression(compute_score=True, n_iter=30).fit(X, y)
df = pd.DataFrame(
    {
        "Weights of true generative process": true_weights,
        "ARDRegression": ard.coef_,
        "BayesianRidge": brr.coef_,
        "LinearRegression": olr.coef_,
    }
)

真の係数と推定係数をプロットする

各モデルの係数を真の生成モデルの重みと比較します。

import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib.colors import SymLogNorm

plt.figure(figsize=(10, 6))
ax = sns.heatmap(
    df.T,
    norm=SymLogNorm(linthresh=10e-4, vmin=-80, vmax=80),
    cbar_kws={"label": "coefficients' values"},
    cmap="seismic_r",
)
plt.ylabel("linear model")
plt.xlabel("coefficients")
plt.tight_layout(rect=(0, 0, 1, 0.95))
_ = plt.title("Models' coefficients")

限界対数尤度をプロットする

両モデルの限界対数尤度をプロットします。

import numpy as np

ard_scores = -np.array(ard.scores_)
brr_scores = -np.array(brr.scores_)
plt.plot(ard_scores, color="navy", label="ARD")
plt.plot(brr_scores, color="red", label="BayesianRidge")
plt.ylabel("Log-likelihood")
plt.xlabel("Iterations")
plt.xlim(1, 30)
plt.legend()
_ = plt.title("Models log-likelihood")

合成データセットの生成

入力特徴量の非線形関数であるターゲットを作成します。標準一様分布に従うノイズを追加します。

from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler

rng = np.random.RandomState(0)
n_samples = 110

## データをソートして、後でプロットを簡単にする
X = np.sort(-10 * rng.rand(n_samples) + 10)
noise = rng.normal(0, 1, n_samples) * 1.35
y = np.sqrt(X) * np.sin(X) + noise
full_data = pd.DataFrame({"input_feature": X, "target": y})
X = X.reshape((-1, 1))

## 外挿
X_plot = np.linspace(10, 10.4, 10)
y_plot = np.sqrt(X_plot) * np.sin(X_plot)
X_plot = np.concatenate((X, X_plot.reshape((-1, 1))))
y_plot = np.concatenate((y - noise, y_plot))

回帰モデルのフィッティング

潜在的に過適合するために、10 次の多項式を試します。ただし、ベイズ線形モデルは多項式係数のサイズを正則化します。ARDRegression と BayesianRidge のデフォルト設定はfit_intercept=Trueなので、PolynomialFeatures は追加のバイアス特徴を導入しないはずです。return_std=Trueを設定することで、ベイズ回帰モデルはモデルパラメータの事後分布の標準偏差を返します。

ard_poly = make_pipeline(
    PolynomialFeatures(degree=10, include_bias=False),
    StandardScaler(),
    ARDRegression(),
).fit(X, y)
brr_poly = make_pipeline(
    PolynomialFeatures(degree=10, include_bias=False),
    StandardScaler(),
    BayesianRidge(),
).fit(X, y)

y_ard, y_ard_std = ard_poly.predict(X_plot, return_std=True)
y_brr, y_brr_std = brr_poly.predict(X_plot, return_std=True)

スコアの標準誤差付きの多項式回帰のプロット

誤差バーは、照会ポイントの予測ガウス分布の 1 標準偏差を表します。両モデルのデフォルトパラメータを使用した場合、ARD 回帰が真の値を最も良く捉えていることに注意してください。ただし、ベイジアンリッジのlambda_initハイパーパラメータをさらに減らすことで、そのバイアスを低減できます。最後に、多項式回帰の固有の制限のため、両モデルは外挿時に失敗します。

ax = sns.scatterplot(
    data=full_data, x="input_feature", y="target", color="black", alpha=0.75
)
ax.plot(X_plot, y_plot, color="black", label="Ground Truth")
ax.plot(X_plot, y_brr, color="red", label="BayesianRidge with polynomial features")
ax.plot(X_plot, y_ard, color="navy", label="ARD with polynomial features")
ax.fill_between(
    X_plot.ravel(),
    y_ard - y_ard_std,
    y_ard + y_ard_std,
    color="navy",
    alpha=0.3,
)
ax.fill_between(
    X_plot.ravel(),
    y_brr - y_brr_std,
    y_brr + y_brr_std,
    color="red",
    alpha=0.3,
)
ax.legend()
_ = ax.set_title("Polynomial fit of a non-linear feature")

まとめ

この実験では、合成データセットを使って 2 つの異なるベイズ回帰モデルを比較します。実験の最初の部分では、真の係数に対するモデルの係数を比較するために、最小二乗法（OLS）モデルをベースラインとして使用します。最後のセクションでは、多項式特徴展開を使ってXとyの間の非線形関係をフィットさせるため、ARD とベイジアンリッジ回帰の予測と不確定性をプロットします。