PRML-3.3 贝叶斯线性回归
记号说明
x_train
y_train
(basis function)
bias parameter
摩尔彭罗斯伪逆(Moore-Penrose pseudo-inverse)
3.3.1
引入模型参数的先验概率分布,吧噪声精度参数当做已知参数
于是对应的共轭先验是形式为
均值和方差分别为$ m_0, S_0 S_0$是已知参数)
由于共轭的性质,后验也是高斯分布
其中
证明
贝叶斯定理
代入3.10 ,3.48
有
其中
,
,同色的相结合
3.50式
3.51式
用配方法,不确定可以反向乘开来验证
,因为后面一项和w没有关系
方法2 利用前面2.3.3 章节的结论
对于$ x y x $的条件高斯分布:
那么$ y x y $的条件高斯分布为:
其中
0均值的先验
后验分布的对数是由对数似然与先验的对数的和给出的一个关于的函数,形式为:
关于来最大化这个后验分布等价于最小化加上一个二次正则项的平方和误差函数,正则项对应的式(3.27)。
贝叶斯线性回归的顺序学习
看图说话
0.模型选用简单的
1.如果数据点时顺序到达的,那么任何一个阶段的后验概率分布都可以看成后续数据点的先验
2.假设先验的和后验的都是已知的
3.当新数据点被观测到的时候,当前的后验分布变成了先验分布
4.第一行对应于未观测到任何数据点之前的情况,因为w的先验假设是正态的,所以第一行中间一列的图像是正态的,第一行右边一列是6个y_test的值,用公式,其中都是先验采样出来的,也就是这一格都是正态分布采样出来的,看起来杂乱无章
5.左中两列,都有一个白色的+号,表示的是的真实数据
6.左列是似然函数,一行二列的先验 二行一列的似然 = 二行二列的后验,下面的都是以此类推
7.第二行是观测到第一个观测数据后的更新后验的情况
8.第三行是观测到第二个观测后的情况
9.观测到20个数据的情况
证明-习题3.8
Solution.
记 ,则已知 个样本的后验估计可以表示为
记 ,则
故
另外
即得到了增加第 个样本时的更新公式。可以看到,递推式和通项公式形式上十分接近。
其他形式的参数先验也可以被考虑。例如,我们可以推广高斯先验分布,得到:
当时对应着高斯分布,并且只有在这种情形下的先验才与式(3.10)给出的似然函数共轭。找到的后验概率分布的最大值等价于找到正则化误差函数(3.29)的最小值。在高斯先验的情况下,后验概率分布的众数等于均值,但是当时这个性质就不成立了。
python实现
点击查看代码
import numpy as np
from scipy.stats import multivariate_normal
import matplotlib.pyplot as plt
%matplotlib inline
from prml.preprocess import GaussianFeature, PolynomialFeature, SigmoidalFeature
from prml.linear import (
BayesianRegression,
EmpiricalBayesRegression,
LinearRegression,
RidgeRegression
)
np.random.seed(1234)
def create_toy_data(func, sample_size, std, domain=[0, 1]):
x = np.linspace(domain[0], domain[1], sample_size)
np.random.shuffle(x)
t = func(x) + np.random.normal(scale=std, size=x.shape)
return x, t
# 自定义图3.7
def linear(x):
return -0.3 + 0.5 * x
x_train, y_train = create_toy_data(linear, 20, 0.1, [-1, 1])
x = np.linspace(-1, 1, 100)
w0, w1 = np.meshgrid(
np.linspace(-1, 1, 100),
np.linspace(-1, 1, 100))
w = np.array([w0, w1]).transpose(1, 2, 0)
feature = PolynomialFeature(degree=1)
X_train = feature.transform(x_train)
X = feature.transform(x)
model = BayesianRegression(alpha=1., beta=100.)
for begin, end in [[0, 0], [0, 1], [1, 2], [2, 3], [3, 20]]:
model.fit(X_train[begin: end], y_train[begin: end])
plt.subplot(1, 2, 1)
plt.scatter(-0.3, 0.5, s=200, marker="x")
plt.contour(w0, w1, multivariate_normal.pdf(w, mean=model.w_mean, cov=model.w_cov))
plt.gca().set_aspect('equal')
plt.xlabel("$w_0$")
plt.ylabel("$w_1$")
plt.title("prior/posterior")
plt.subplot(1, 2, 2)
plt.scatter(x_train[:end], y_train[:end], s=100, facecolor="none", edgecolor="steelblue", lw=1)
plt.plot(x, model.predict(X, sample_size=6), c="orange")
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
点击查看代码
import numpy as np
from prml.linear.regression import Regression
class BayesianRegression(Regression):
"""
贝叶斯回归模型
w ~ N(w|0, alpha^(-1)I)
y = X @ w
t ~ N(t|X @ w, beta^(-1))
"""
def __init__(self, alpha: float = 1., beta: float = 1.):
self.alpha = alpha
self.beta = beta
self.w_mean = None
self.w_precision = None
def _is_prior_defined(self) -> bool:
return self.w_mean is not None and self.w_precision is not None
def _get_prior(self, ndim: int) -> tuple: # -> tuple 表示返回一个tuple
'''
:param ndim:
:return: w_mean 先验的均值
w_precision 先验的精度
'''
if self._is_prior_defined():
return self.w_mean, self.w_precision
else:
return np.zeros(ndim), self.alpha * np.eye(ndim)
def fit(self, X: np.ndarray, t: np.ndarray):
"""
给定训练集,更新贝叶斯参数
Parameters
----------
X : (N, n_features) np.ndarray
训练集输入
t : (N,) np.ndarray
训练集输出
"""
# 先验的均值,先验的精度
mean_prev, precision_prev = self._get_prior(
np.size(X, 1)) # mean_prev =m_0 precision_prev =S_0= \alpha I 见PRML 3.48
w_precision = precision_prev + self.beta * X.T @ X # 这里是计算 w_precision=S_N^{-1}=S_0^{-1}+\beta X^T X,PRML 3.54
# 计算m(x) 也就是后验的均值
w_mean = np.linalg.solve(
w_precision,
precision_prev @ mean_prev + self.beta * X.T @ t
# w_mean=m_N ,解方程 S_N^{-1} w_mean =S_0^{-1} * m_0 +\beta X^T t ,PRML 3.53
)
self.w_mean = w_mean
self.w_precision = w_precision
self.w_cov = np.linalg.inv(self.w_precision) # w_cov = S_N , 协方差矩阵
def predict(self, X: np.ndarray, return_std: bool = False, sample_size: int = None):
"""
return mean (and standard deviation) of predictive distribution
Parameters
----------
X : (N, n_features) np.ndarray
independent variable
return_std : bool, optional
flag to return standard deviation (the default is False)
sample_size : int, optional
number of samples to draw from the predictive distribution
(the default is None, no sampling from the distribution)
Returns
-------
y : (N,) np.ndarray
mean of the predictive distribution
y_std : (N,) np.ndarray
standard deviation of the predictive distribution
y_sample : (N, sample_size) np.ndarray
samples from the predictive distribution
"""
# 得到拟合的正态分布后,对这个正态分布采样
if sample_size is not None:
w_sample = np.random.multivariate_normal(
self.w_mean, self.w_cov, size=sample_size
)
y_sample = X @ w_sample.T
return y_sample
y = X @ self.w_mean
if return_std:
y_var = 1 / self.beta + np.sum(X @ self.w_cov * X, axis=1)
y_std = np.sqrt(y_var)
return y, y_std
return y
3.3.2 预测分布
预测分布如下
公式3.57
目标变量的条件分布由式(3.8)给出,后验分布由式(3.49)给出。
目标变量的条件分布:
后验分布
得到的预测分布的形式为:
其中预测分布的方差是
当额外的数据点被观测到的时候,后验概率分布会变窄,从而可以证明
证明-习题3.11
Hint.
记
根据提示(3.110恒等式),,令,则有,则有。
3.3.3 贝叶斯模型的等价核解释
预测均值E(y_test)与训练集目标值y_train之间的线性关系3.60 3.61
其中函数
被称为平滑矩阵(smoother matrix)或等价核(equivalent kernel)。像这样,通过训练集的目标值的线性组合来进行预测的回归函数被称为线性平滑(linear smoother)
用核函数表示线性回归是解决回归问题的另一种方法。在给定观测数据集的条件下,直接定义一个局部的核函数,而不是引入一组基函数来隐式地定义了一个等价的核,然后使用这个核函数对新的输入变量$ x x $值做预测。可以证明权值的和为1,即:
对于所有成立。
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· C# 集成 DeepSeek 模型实现 AI 私有化(本地部署与 API 调用教程)
· spring官宣接入deepseek,真的太香了~