9、贝叶斯岭回归

import numpy as np

import matplotlib.pyplot as plt

from scipy import stats

from sklearn.linear_model import BayesianRidge, LinearRegression

plt.rcParams['font.sans-serif'] = ['SimHei']

plt.rcParams['axes.unicode_minus'] = False

# 用高斯加权生成模拟数据

np.random.seed(0)

n_samples, n_features = 100, 100

X = np.random.randn(n_samples, n_features)  # Create Gaussian data

# 创建精度λ为4的权重。

lambda_ = 4.

w = np.zeros(n_features)

# 只保留10个权重

relevant_features = np.random.randint(0, n_features, 10)

for i in relevant_features:

    w[i] = stats.norm.rvs(loc=0, scale=1. / np.sqrt(lambda_))

# 创建噪声的精度α为50。

alpha_ = 50.

noise = stats.norm.rvs(loc=0, scale=1. / np.sqrt(alpha_), size=n_samples)

# 创建目标

y = np.dot(X, w) + noise

# 贝叶斯岭回归与OLS的拟合比较

clf = BayesianRidge(compute_score=True)

clf.fit(X, y)

ols = LinearRegression()

ols.fit(X, y)

# 绘制真实权重、估计权重、权重直方图，以及带有标准差的预测

lw = 2

plt.figure(figsize=(6, 5))

plt.title("模型权重")

plt.plot(clf.coef_, color='lightgreen', linewidth=lw,

        label="Bayesian Ridge estimate")

plt.plot(w, color='gold', linewidth=lw, label="Ground truth")

plt.plot(ols.coef_, color='navy', linestyle='--', label="OLS estimate")

plt.xlabel("Features")

plt.ylabel("Values of the weights")

plt.legend(loc="best", prop=dict(size=12))

plt.figure(figsize=(6, 5))

plt.title("权重直方图")

plt.hist(clf.coef_, bins=n_features, color='gold', log=True,

        edgecolor='black')

plt.scatter(clf.coef_[relevant_features], np.full(len(relevant_features), 5.),

            color='navy', label="Relevant features")

plt.ylabel("Features")

plt.xlabel("Values of the weights")

plt.legend(loc="upper left")

plt.figure(figsize=(6, 5))

plt.title("边际对数似然")

plt.plot(clf.scores_, color='navy', linewidth=lw)

plt.ylabel("Score")

plt.xlabel("Iterations")

# 绘制多项式回归的预测

def f(x, noise_amount):

    y = np.sqrt(x) * np.sin(x)

    noise = np.random.normal(0, 1, len(x))

    return y + noise_amount * noise

degree = 10

X = np.linspace(0, 10, 100)

y = f(X, noise_amount=0.1)

clf_poly = BayesianRidge()

clf_poly.fit(np.vander(X, degree), y)

X_plot = np.linspace(0, 11, 25)

y_plot = f(X_plot, noise_amount=0)

y_mean, y_std = clf_poly.predict(np.vander(X_plot, degree), return_std=True)

plt.figure(figsize=(6, 5))

plt.errorbar(X_plot, y_mean, y_std, color='navy',

            label="Polynomial Bayesian Ridge Regression", linewidth=lw)

plt.plot(X_plot, y_plot, color='gold', linewidth=lw,

        label="Ground Truth")

plt.ylabel("Output y")

plt.xlabel("Feature X")

plt.legend(loc="lower left")

plt.title("贝叶斯岭回归")

plt.show()