PRML-3.3 贝叶斯线性回归
记号说明
\(1.输入集\textbf{X}=\{x_1,...,x_N\}是N个观测值,某一个观测\{x_n\},其中n=1,2,...,N,通俗讲就是\)x_train\(,或者文中称为\mathcal{D}\)
\(2.观测对应的目标值\textbf{t}=\{t_1,...,t_n\},通俗讲就是\)y_train
\(3.模型函数 t=y(x),输入变量x，输出对应的t的预测\)
\(4.预测分布p(t|x)\)
\(5.线性模型y(x,w)=w_0+w_1x_1+...+w_Dx_D=w_0+\sum\limits_{j=1}^{M-1}w_j\phi(x)=w^T\phi(x)\)
\(6.M 多项式拟合的最高阶数\)
\(7.\phi_j(x)被称为基函数\)（basis function）\(j=0,...,M-1\)
\(8.w_0被称为偏置参数\) bias parameter
\(9.w = (w_0,...,w_{M-1})^T\)
$10.\phi(x) = (\phi_0,...,\phi_{M-1})^T $
\(11.\phi_j(x)=x^j,多项式基函数\)
\(12.\beta,先验为高斯分布的精度,p(t|x,w,\beta) = \mathcal{N}(t|y(x,w), \beta^{-1})\)
\(13.\boldsymbol{\Phi}=\left(\begin{array}{cccc}
\phi_{0}\left(\boldsymbol{x}_{1}\right) & \phi_{1}\left(\boldsymbol{x}_{1}\right) & \cdots & \phi_{M-1}\left(\boldsymbol{x}_{1}\right) \\
\phi_{0}\left(\boldsymbol{x}_{2}\right) & \phi_{1}\left(\boldsymbol{x}_{2}\right) & \cdots & \phi_{M-1}\left(\boldsymbol{x}_{2}\right) \\
\vdots & \vdots & \ddots & \vdots \\
\phi_{0}\left(\boldsymbol{x}_{N}\right) & \phi_{1}\left(\boldsymbol{x}_{N}\right) & \cdots & \phi_{M-1}\left(\boldsymbol{x}_{N}\right)
\end{array}\right)_{N\times M}\)
\(14. \Phi^+ \equiv (\Phi^T\Phi)^{-1}\Phi^T,\)摩尔彭罗斯伪逆（Moore-Penrose pseudo-inverse）
\(15.E_{D}(\boldsymbol{w})=\frac{1}{2} \sum_{n=1}^{N}\left\{t_{n}-w_{0}-\sum_{j=1}^{M-1} w_{j} \phi_{j}\left(\boldsymbol{x}_{n}\right)\right\}^{2}= \frac{1}{2}\sum\limits_{n=1}^N\{t_n-w^T\phi(x_n)\}^2,误差函数\)
\(15.E_W(w) = \frac{1}{2}w^Tw ,正则化项\)
3.3.1
引入模型参数\(w\)的先验概率分布,吧噪声精度参数\(\beta\)当做已知参数
于是对应的共轭先验是形式为
均值和方差分别为$ m_0, S_0 \(的高斯分布。(\)S_0$是已知参数)
由于共轭的性质，后验也是高斯分布
其中
证明
贝叶斯定理\(p(w|t)\propto p(t|w)p(w)\)
代入3.10 ,3.48
\(p(\textbf{t}|\textbf{X},w,\beta) = \prod\limits_{n=1}^N\mathcal{N}(t_n|w^T\phi(x_n),\beta^{-1})\)
\(p(w) = \mathcal{N}(w|m_0,S_0)\)
有
\(p(w|t)\propto exp[-\frac{\beta}{2}(t_1-w^T\phi_1)^2] exp[-\frac{\beta}{2}(t_2-w^T\phi_2)^2]...exp[-\frac{\beta}{2}(t_n-w^T\phi_n)^2]exp[-\frac{1}{2}(w-m_0)^TS_0^{-1}(w-m_0)]\)
\(=exp[-\frac{\beta}{2}\{(t_1-w^T\phi_1)^2+...+(t_n-w^T\phi_n)^2\}]exp[-\frac{1}{2}(w-m_0)^TS_0^{-1}(w-m_0)]\)
\(=exp[-\frac{\beta}{2}\begin{pmatrix}
t_1-w^T\phi_1\\
...\\
t_n-w^T\phi_n
\end{pmatrix}^T\begin{pmatrix}
t_1-w^T\phi_1\\
...\\
t_n-w^T\phi_n
\end{pmatrix}]exp[-\frac{1}{2}(w-m_0)^TS_0^{-1}(w-m_0)]\)
\(=exp[-\frac{\beta}{2}(T-\Phi w)^T(T-\Phi w)]exp[-\frac{1}{2}(w-m_0)^TS_0^{-1}(w-m_0)]\)
其中
\(T=\begin{pmatrix}
t_1\\
...\\
t_n
\end{pmatrix}^T\)，\(\Phi w^T=\begin{pmatrix}
\phi_1^Tw\\
...\\
\phi_n^Tw
\end{pmatrix}=\begin{pmatrix}
w^T\phi_1\\
...\\
w^T\phi_n
\end{pmatrix}^T\)
\(=exp[-\frac{1}{2}(\beta T^TT-\color{red}{\beta T^T\Phi w}-\color{green}{\beta(\Phi w)^TT }+ \color{blue}{\beta w^T\Phi^T\Phi^T w}+ \color{blue}{w^TS_0^{-1}w} -\color{red}{m_0^TS_0^{-1}w}-\color{green}{w^TS_0^{-1}m_0 }+m_0^TS_0^{-1}m_0)]\)
\(=exp[-\frac{1}{2}(\color{blue}{w^T(S_0^{-1}+\beta \Phi^T\Phi)w}-\color{green}{w^T(S_0^{-1}m_0+\beta\Phi^TT)}-\color{red}{(\beta T^T\Phi + m_0^TS_0^{-1})w}+m_0^TS_0^{-1}m_0+\beta T^T T)]\),同色的相结合
3.50式\(m_N=S_N(S_0^{-1}m_0+\beta \Phi^TT)\)
3.51式\(S_N^{-1}=S_0^{-1}+\beta\Phi^T\Phi\)
用配方法，不确定可以反向乘开来验证
\(=exp[-\frac{1}{2}(w-m_N)^TS_N^{-1}(w-m_N)]exp[-\frac{1}{2}(m_0^TS_0^{-1}m_0+\beta ^TT+m_N^TS_N^Tm_N)]\)
\(=C\times exp[-\frac{1}{2}(w-m_N)^TS_N^{-1}(w-m_N)]\),因为后面一项和w没有关系
方法2 利用前面2.3.3 章节的结论
对于$ x \(的边缘高斯分布和\) y \(关于\) x $的条件高斯分布：
那么$ y \(的边缘分布和\) x \(关于\) y $的条件高斯分布为：
其中
0均值的先验
后验分布的对数是由对数似然与先验的对数的和给出的一个关于\(w\)的函数，形式为：
关于\(w\)来最大化这个后验分布等价于最小化加上一个二次正则项的平方和误差函数，正则项对应$$ \lambda = \alpha/\beta $$的式（3.27）。
贝叶斯线性回归的顺序学习
看图说话
0.模型选用简单的\(y(x,w)=w_0+w_1x\)
1.如果数据点时顺序到达的，那么任何一个阶段的后验概率分布都可以看成后续数据点的先验
2.假设先验的\(\beta\)和后验的\(\alpha\)都是已知的
3.当新数据点被观测到的时候，当前的后验分布变成了先验分布
4.第一行对应于未观测到任何数据点之前的情况，因为w的先验假设是正态的，所以第一行中间一列的图像是正态的，第一行右边一列是6个y_test的值,用公式\(y(x,w)=w_0+w_1x\),其中\(w_0,w_1\)都是先验采样出来的，也就是这一格都是正态分布采样出来的，看起来杂乱无章
5.左中两列，都有一个白色的+号，表示的是\(w_0,w_1\)的真实数据
6.左列是似然函数，一行二列的先验 \(\times\) 二行一列的似然 = 二行二列的后验，下面的都是以此类推
7.第二行是观测到第一个观测数据后的更新后验的情况
8.第三行是观测到第二个观测后的情况
9.观测到20个数据的情况
证明-习题3.8
Solution.
记 \(\Phi_N=[\phi_1^T;\dots;\phi_N^T],\mathbf{t}_N=[t_1,\dots,t_N]^T\)，则已知 \(N\) 个样本的后验估计可以表示为
记 \(\Phi_{N+1}=[\Phi_N;\phi_{N+1}],\mathbf{t}_{N+1}=[\mathbf{t}_N^T,t_{N+1}]^T\)，则
故
另外
即得到了增加第 \(N+1\) 个样本时的更新公式。可以看到，递推式和通项公式形式上十分接近。
其他形式的参数先验也可以被考虑。例如，我们可以推广高斯先验分布，得到：
当\(q = 2\)时对应着高斯分布，并且只有在这种情形下的先验才与式（3.10）给出的似然函数共轭。找到\(w\)的后验概率分布的最大值等价于找到正则化误差函数（3.29）的最小值。在高斯先验的情况下，后验概率分布的众数等于均值，但是当\(q \neq 2\)时这个性质就不成立了。
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
\(目标\)
$ p(t|T, \alpha, \beta) = \int p(t|w,\beta)p(w|T,\alpha,\beta)dw \tag{3.57} $
\(1.p(t,T, \alpha, \beta)=\int p(t,T, \alpha, \beta,w)dw\)
\(=\int p(t|T, \alpha, \beta,w)p(T, \alpha, \beta,w)dw\)
\(=\int p(t|T, \alpha, \beta,w)p(w|T, \alpha, \beta)p(T, \alpha, \beta)dw\)
\(2.p(t,T, \alpha, \beta)=p(t|T, \alpha, \beta)p(T, \alpha, \beta)\)
\(1=2，p(x,X,T)和w无关，可以约去，所以有\)
\(p(t|T, \alpha, \beta)=\int p(t|T, \alpha, \beta,w)p(w|T, \alpha, \beta)dw\)
\(其中p(t|T, \alpha, \beta,w)，t是验证数据集参数，和\alpha,T 这两个训练数据集参数无关,所以 p(t|T, \alpha, \beta,w)=p(t| \beta,w)\)
\(最后得到\)
\(p(t|T, \alpha, \beta) = \int p(t|w,\beta)p(w|T,\alpha,\beta)dw\)
目标变量的条件分布由式（3.8）给出，后验分布由式（3.49）给出。
目标变量的条件分布：
后验分布
$ p(w|t) = \mathcal{N}(w|m_N,S_N) \tag{3.49} $
得到的预测分布的形式为：
其中预测分布的方差\(\sigma_N^2(x)\)是
$ \sigma_N^2(x) = \frac{1}{\beta} + \phi(x)^TS_N\phi(x) \tag{3.59} $
当额外的数据点被观测到的时候，后验概率分布会变窄，从而可以证明\(\sigma_{N+1}^2(x) \leq \sigma_N^2(x)\)
证明-习题3.11
Hint.
记\(\Phi_N = [\phi^T_0;\dots;\phi^T_N]\)
根据提示(3.110恒等式)\((M+v v^T)^{-1}=M^{-1}-\frac{(M^{-1}v)(v^T M^{-1})}{1+v^T M^{-1}v}\)，\((M+v v^T)^{-1}\preceq M^{-1}\)，令\(M=S_N^{-1}\)，则有\(S_{N+1}\preceq S_N\)，则有\(\sigma^2_{N+1}\leq \sigma^2_{N}\)。
3.3.3 贝叶斯模型的等价核解释
预测均值E(y_test)与训练集目标值y_train之间的线性关系3.60 3.61
其中函数
$ k(x,x') = \beta\phi(x)^TS_N\phi(x') \tag{3.62} $
被称为平滑矩阵（smoother matrix）或等价核（equivalent kernel）。像这样，通过训练集的目标值的线性组合来进行预测的回归函数被称为线性平滑（linear smoother）
用核函数表示线性回归是解决回归问题的另一种方法。在给定观测数据集的条件下，直接定义一个局部的核函数，而不是引入一组基函数来隐式地定义了一个等价的核，然后使用这个核函数对新的输入变量$ x \(做预测。这就是将在6.4节详细讨论的用于回归问题（以及分类问题）的一个称为高斯过程（Gaussian process）的实用框架。 我们已经看到，有效核定义了训练数据集里的目标值组合的权值，然后对新的\) x $值做预测。可以证明权值的和为1，即：
$ \sum\limits_{n=1}^Nk(x,x_n) = 1 \tag{3.64} $
对于所有$ x $成立。
