在本文中，首先使用普通最小二乘（OLS）模型作为基线，以比较不同模型系数与真实系数的关系。随后，展示了如何通过迭代最大化观测值的边际对数似然来估计这类模型。在最后部分，使用多项式特征扩展来拟合X和y之间的非线性关系，并为ARD和贝叶斯岭回归绘制了预测和不确定性。
生成了一个合成数据集，其中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({
"真实生成过程的权重": true_weights,
"ARD回归": ard.coef_,
"贝叶斯岭回归": brr.coef_,
"线性回归": olr.coef_,
})
现在将比较每个模型的系数与真实生成模型的权重。由于添加了噪声，没有模型能够恢复真实的权重。实际上，所有模型都有超过10个非零系数。与OLS估计器相比，使用贝叶斯岭回归的系数稍微向零偏移，这稳定了它们。ARD回归提供了一个更稀疏的解决方案：一些非信息系数被精确地设置为零，而其他一些则更接近零。一些非信息系数仍然存在并保留大值。
绘制了ARD和贝叶斯岭回归的边际对数似然，以展示模型如何随着迭代次数的增加而最小化对数似然。
import numpy as np
ard_scores = -np.array(ard.scores_)
brr_scores = -np.array(brr.scores_)
# 绘制代码省略...
实际上，两种模型都最小化了对数似然，直到由max_iter参数定义的任意截止值。
创建了一个目标，它是输入特征的非线性函数。添加了遵循标准均匀分布的噪声。
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({
"输入特征": X,
"目标": 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)
# 绘图代码省略...
误差条表示查询点的预测高斯分布的一个标准差。注意，当使用两个模型中的默认参数时，ARD回归最好地捕获了真实情况，但是进一步减少贝叶斯岭回归的lambda_init超参数可以减少其偏差（参见示例“使用贝叶斯岭回归进行曲线拟合”）。最后，由于多项式回归的固有局限性，两种模型在进行外推时都失败了。