注意
转到末尾 下载完整的示例代码。或通过 JupyterLite 或 Binder 在您的浏览器中运行此示例
比较线性贝叶斯回归器#
此示例比较了两种不同的贝叶斯回归器
a 贝叶斯岭回归
在第一部分，我们使用普通最小二乘法 (OLS) 模型作为基线，用于比较模型系数与真实系数。然后，我们展示了这些模型的估计是通过迭代地最大化观测值的边际对数似然函数来完成的。
在最后一节中，我们使用多项式特征扩展来拟合X和y之间的非线性关系，并绘制ARD和贝叶斯岭回归的预测值和不确定性。
# Authors: The scikit-learn developers
# SPDX-License-Identifier: BSD-3-Clause
模型恢复真实权重的稳健性#
生成合成数据集#
我们生成一个数据集，其中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,
)
拟合回归器#
我们现在拟合两个贝叶斯模型和OLS模型，以便稍后比较模型系数。
import pandas as pd
from sklearn.linear_model import ARDRegression, BayesianRidge, LinearRegression
olr = LinearRegression().fit(X, y)
brr = BayesianRidge(compute_score=True, max_iter=30).fit(X, y)
ard = ARDRegression(compute_score=True, max_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")
由于添加了噪声，没有一个模型能够恢复真实权重。实际上，所有模型的非零系数都超过10个。与OLS估计量相比，使用贝叶斯岭回归的系数略微向零偏移，这使它们更加稳定。ARD回归提供更稀疏的解：一些非信息性系数被精确地设置为零，而其他系数则向零偏移。一些非信息性系数仍然存在并保持较大的值。
绘制边际对数似然函数#
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")
实际上，这两个模型都最小化了对数似然函数，直到由max_iter参数定义的任意截止点。
具有多项式特征扩展的贝叶斯回归#
生成合成数据集#
我们创建一个目标，它是输入特征的非线性函数。添加服从标准均匀分布的噪声。
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
rng = np.random.RandomState(0)
n_samples = 110
# sort the data to make plotting easier later
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))
# extrapolation
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)
绘制多项式回归及其分数的标准误差#
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")
误差线表示查询点预测高斯分布的一个标准差。请注意，在使用两个模型的默认参数时，ARD回归能够最好地捕捉真实值，但是进一步减小贝叶斯岭的lambda_init超参数可以减少其偏差（参见示例使用贝叶斯岭回归进行曲线拟合）。最后，由于多项式回归的内在局限性，这两个模型在进行外推时都会失败。
脚本总运行时间：（0分钟0.714秒）
相关示例
