import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
%matplotlib inline
之前的讨论中,我们使用统计上的极大似然估计的方法来找到模型的参数,我们看到,为了得到较好的结果,我们需要控制模型的复杂度,即选择一个合适的基函数个数;或者通过加上正则项来解决。
通常模型的复杂性不能通过极大似然函数来单独决定,因为极大似然估计通常都伴随着过拟合的现象,我们可以考虑使用独立的验证集来控制模型复杂度,但这通常是非常昂贵的。
我们考虑用 Bayes 的角度来看待线性回归的问题。
考虑单目标变量 $t$ 的情况,我们先考虑给我们的参数 $\bf w$ 加上一个先验分布,这里我们先认为噪音的参数 $\beta$ 是已知的。在 $N$ 个数据点的情况下,之前我们得到 $p(\mathsf t|\mathbf w)$ 的分布:
$$ p(\mathsf t|\mathbf{X,w},\beta) = \prod_{n=1}^N \mathcal N(t_n|\mathbf w^\top\mathbf \phi(\mathbf x_n),\beta^{-1}) $$是一个 $\bf w$ 二次型的指数形式,因此,我们考虑使用高斯分布的先验:
$$ p(\mathbf w)=\mathcal N(\mathbf w|\mathbf m_0,\mathbf S_0) $$利用${\rm posterior} \propto {\rm prior} \times {\rm likelihood}$ 以及共轭性,我们知道后验分布也是一个高斯分布:
$$ p(\mathbf w|\mathsf t)=\mathcal N(\mathbf w|\mathbf m_N,\mathbf S_N) $$不难计算得到:
$$ \begin{aligned} \mathbf m_N & = \mathbf S_N(\mathbf S_0^{-1}\mathbf m_0+\beta\mathbf\Phi^\top\mathsf t) \\ \mathbf S_N^{-1} & = \mathbf S_0^{-1} + \beta\mathbf{\Phi^\top\Phi} \end{aligned} $$其中
$$ \mathbf\Phi = \begin{pmatrix} \mathbf\phi(\mathbf x_1)^\top \\ \mathbf\phi(\mathbf x_2)^\top \\ \vdots \\ \mathbf\phi(\mathbf x_N)^\top \\ \end{pmatrix} = \begin{pmatrix} \phi_0(\mathbf x_1) & \phi_1(\mathbf x_1) & \cdots & \phi_{M-1}(\mathbf x_1) \\ \phi_0(\mathbf x_2) & \phi_1(\mathbf x_2) & \cdots & \phi_{M-1}(\mathbf x_2) \\ \vdots & \vdots & \ddots & \vdots\\ \phi_0(\mathbf x_N) & \phi_1(\mathbf x_N) & \cdots & \phi_{M-1}(\mathbf x_n)\\ \end{pmatrix} $$在高斯分布的情况下,均值就是众数,因此,最大后验给出:
$$ \mathbf w_{MAP} = \mathbf m_N $$如果我们考虑一个无限宽的先验:$\mathbf S_0 = \alpha\mathbf I, \alpha\to 0$,那么我们的最大后验结果就是极大似然的结果;类似的,如果 $N=0$,那么后验分布就是先验分布,与之前的讨论类似,我们可以将这个模型改成在线的。
如果我们考虑一个零均值的各向同性的先验:
$$ p(\mathbf w|\alpha) = \mathcal N(\mathbf 0,\alpha\mathbf I) $$我们有:
$$ \begin{aligned} \mathbf m_N & = \beta\mathbf S_N\mathbf\Phi^\top\mathsf t \\ \mathbf S_N^{-1} & = \alpha\mathbf I + \beta\mathbf{\Phi^\top\Phi} \end{aligned} $$此时,后验分布的对数为:
$$ \ln p(w|t) = \underbrace{-\frac{\beta}{2}\sum_{n=1}^N \{t_n-\mathbf{w^\top\phi(x}_n)\}^2}_{\rm likelihood} - \underbrace{\frac{\alpha}{2} \mathbf{w^\top w}}_{\rm prior} + {\rm const} $$因此,在这种情况下,最大后验分布,相当于加了一个二范数正则的最小二乘,正则系数 $\lambda = \alpha/\beta$。
from scipy.stats import uniform, norm
from matplotlib.mlab import bivariate_normal
from scipy.stats import multivariate_normal
from numpy import linalg
a0, a1 = -0.3, 0.5
bt = 25.0
s = np.sqrt(1.0 / bt)
ap = 2.0
xx = uniform.rvs(loc=-1, scale=2, size=(20,))
tt = a0 + a1 * xx + norm.rvs(loc=0, scale=s, size=(20,))
_, axes = plt.subplots(4, 3, figsize=(8, 11))
x, y = np.mgrid[-1:1:.02, -1:1:.02]
nums = [0, 1, 2, 20]
phi = np.vstack((np.ones(20), xx)).T
for idx in range(4):
n = nums[idx]
if idx == 0:
axes[idx][0].axis("off")
axes[idx][0].set_title("likelihood", fontsize="x-large")
axes[idx][1].set_title("prior/posterior", fontsize="x-large")
axes[idx][2].set_title("dataspace", fontsize="x-large")
else:
z1 = 1 / np.sqrt(2*np.pi) / s *np.exp(-(tt[n-1] - x - y * xx[n-1]) ** 2 / 2 / s ** 2)
axes[idx][0].contourf(x, y, z1)
axes[idx][0].set_xticks([-1, 0, 1])
axes[idx][0].set_ylim(-1, 1)
axes[idx][0].set_yticks([-1, 0, 1])
axes[idx][0].text(0.5, -1.2, '$w_0$', fontsize='xx-large')
axes[idx][0].text(-1.4, 0.5, '$w_1$', fontsize='xx-large')
S = linalg.inv(ap * np.eye(2) + bt * np.dot(phi[:n].T, phi[:n]))
mu = bt * np.dot(np.dot(S, phi[:n].T), tt[:n])
z = bivariate_normal(x, y,
mux=mu[0], muy=mu[1],
sigmax=np.sqrt(S[0][0]),
sigmay=np.sqrt(S[1][1]),
sigmaxy=S[0][1])
ww = multivariate_normal.rvs(mean=mu,
cov=S,
size=6)
ww0, ww1 = ww[:, 0], ww[:, 1]
for idy in range(6):
axes[idx][2].plot(x, ww0[idy] + ww1[idy] * x, 'r')
axes[idx][1].contourf(x, y, z)
axes[idx][1].set_xticks([-1, 0, 1])
axes[idx][1].set_ylim(-1, 1)
axes[idx][1].set_yticks([-1, 0, 1])
axes[idx][1].text(0.5, -1.2, '$w_0$', fontsize='xx-large')
axes[idx][1].text(-1.4, 0.5, '$w_1$', fontsize='xx-large')
axes[idx][2].plot(xx[:n], tt[:n], 'o')
axes[idx][2].set_xticks([-1, 0, 1])
axes[idx][2].set_ylim(-1, 1)
axes[idx][2].set_yticks([-1, 0, 1])
axes[idx][2].text(0.5, -1.2, '$x$', fontsize='xx-large')
axes[idx][2].text(-1.4, 0.5, '$y$', fontsize='xx-large')
plt.show()
我们考虑一个单变量 $x$,单目标 $t$ 的情况,并使用一个线性模型:$y(x,\mathbf w)=w_0+w_1x$。我们使用函数 $f(x,\mathbf a)=a_0+a_1x$ 生成一组 $x\in U[-1,1]$ 之间的数据,并加上零均值 $\beta=25$ 的高斯噪声,这里 $a_0 = -0.3, a_1 = 0.5$,令 $\alpha=2$,考虑数据点不断增加时,后验分布的变化情况,如上图所示。
当然我们可以一般化高斯分布先验为(高斯分布对应于 $q=2$):
$$ p(\mathbf w|\alpha)=\left[\frac{q}{2}\left(\frac{\alpha}{2}\right)^{1/q} \frac{1}{\Gamma(1/q)}\right]^M \exp\left(-\frac{\alpha}{2}\sum_{j=1}^M |w_j|^q\right) $$此时 MAP 对应于更一般的正则项形式:
$$ \frac{1}{2}\sum_{n=1}^N \{t_n-\mathbf w^\top\mathbf \phi(\mathbf x_n)\}^2 + \frac \lambda 2 \sum_{j=1}^M |w_j|^q $$在这种情况下,其 MAP 的值不一定为均值。
通常我们更关心对一个新输入 $\bf x$,其预测值 $t$ 的分布:
$$ p(t|{\bf x}, {\sf t}, \alpha, \beta) = \int p(t|\mathbf{x,w}, \beta) p({\bf w}|{\sf t}, \alpha, \beta)d{\bf w} $$因为 $$ \begin{aligned} p(t|\mathbf{x,w},\beta)&=\mathcal N(t|y(\mathbf{x,w}),\beta^{-1}) \\ p(\mathbf w|\mathsf t, \alpha, \beta)&=\mathcal N(\mathbf w|\mathbf m_N,\mathbf S_N) \end{aligned} $$
都是高斯分布,考虑:
$$ y(\mathbf{x,w}) = \mathbf{w^\top \phi(x)} $$我们有:
$$ p(t|{\bf x}, {\sf t}, \alpha, \beta) = \int \mathcal N(t|\mathbf{w^\top \phi(x)},\beta^{-1}) \mathcal N(\mathbf w|\mathbf m_N,\mathbf S_N)d{\bf w} $$对比之前2.3.3节的的结论(给定 $p(x)$ 和 $p(y|x)$ 求 $p(y)$,这里 $\bf w$ 相当于 $x$,$t$ 相当于 $y$),我们有:
$$ p(t|{\bf x}, {\sf t}, \alpha, \beta) =\mathcal N(t|\mathbf m_N^\top \mathbf{\phi(x)}, \sigma^2_N({\bf x})) $$其中
$$ \sigma^2_N({\bf x}) = \frac{1}{\beta} + \mathbf{\phi(x)}^\top \mathbf S_N \mathbf{\phi(x)} $$方差的第一部分表示数据中的噪声,第二部分反映了参数 $\bf w$ 的不确定性。
一半来说,当有新的数据点进来了,后验分布 $p(\mathbf w|\mathsf t, \alpha, \beta)$ 变得更窄,参数 $\bf w$ 的不确定性减小,从而有
$$ \sigma^2_{N+1}({\bf x}) \leq \sigma^2_N({\bf x}) $$当 $N\to\infty$ 是,第二项趋于 0,从而方差趋于 $\frac{1}{\beta}$。
def phi(x, M):
return x[:,None] ** np.arange(M + 1)
fig, axes = plt.subplots(2, 2, figsize=(10, 6))
axes = axes.flatten()
NN = [1, 2, 4, 25]
for N, ax in zip(NN, axes):
# 生成 0,1 之间等距的 N 个 数
x_tr = np.linspace(0, 1, N+2)
x_tr = x_tr[1:-1]
# 计算 t
t_tr = np.sin(2 * np.pi * x_tr) + 0.25 * np.random.randn(N)
# 加正则项的解
M = 8
alpha = 5e-3
beta = 11.1
lam = alpha / beta
phi_x_tr = phi(x_tr, M)
A_0 = phi_x_tr.T.dot(phi_x_tr) + lam * np.eye(M+1)
y_0 = t_tr.dot(phi_x_tr)
# 求解 Aw=y
coeff = np.linalg.solve(A_0, y_0)[::-1]
f = np.poly1d(coeff)
# 绘图
xx = np.linspace(0, 1, 500)
# Bayes估计的均值和标准差
S = np.linalg.inv(A_0 * beta)
m_xx = beta * phi(xx, M).dot(S).dot(y_0)
s_xx = np.sqrt(1 / beta + phi(xx, M).dot(S).dot(phi(xx, M).T).diagonal())
ax.plot(x_tr, t_tr, 'co')
ax.plot(xx, np.sin(2 * np.pi * xx), 'g')
ax.plot(xx, f(xx), 'r')
ax.fill_between(xx, m_xx-s_xx, m_xx+s_xx, color="pink")
ax.set_xlim(-0.1, 1.1)
ax.set_ylim(-1.5, 1.5)
ax.set_xticks([0, 1])
ax.set_yticks([-1, 0, 1])
ax.set_xlabel("$x$", fontsize="x-large")
ax.set_ylabel("$t$", fontsize="x-large")
plt.show()
上图是 $N = 1,2,4,25$ 时,使用线性基,方差的变化情况,红色部分的宽度表示两倍标准差。
上面讨论的是 $\beta$ 已知的情况,如果两者都不知道,那么我们与 2.3.6 的情况一样,引入一个 Gaussian-gamma
分布,然后得到预测值的分布为学生 t 分布。
将线性基函数模型参数的后验均值带入模型:
$$ \mathbf m_N = \beta\mathbf S_N\mathbf\Phi^\top\mathsf t $$我们有:
$$ y(\mathbf{x, m}_N) = \mathbf m_N^\top \mathbf{\phi(x)} = \beta \mathbf{\phi(x)}^\top \mathbf S_N\mathbf\Phi^\top\mathsf t = \sum_{n=1}^N \beta \mathbf{\phi(x)}^\top \mathbf S_N\mathbf \phi(\mathbf x_n) t_n $$我们将$y$ 看成是一个 $t_n$ 的线性组合,我们有:
$$ y(\mathbf{x, m}_N) = \sum_{n=1}^N k(\mathbf x, \mathbf x_n) t_n $$其中:
$$ k(\mathbf x, \mathbf x') = \beta \mathbf{\phi(x)}^\top \mathbf S_N\mathbf \phi(\mathbf x') $$这就是所谓的平滑矩阵(smoother matrix
)或者等价核(equivalent kernel
)。上面的这种使用训练集预测值 $t_n$ 的线性组合的情况叫做线性平滑(linear smoothers
)。
注意每个等价核函数的计算要使用到训练集中的每个 $\mathbf x_n$(因为 $\mathbf S_N$)
此外,我们有:
$$ {\rm cov}[y(\mathbf x), y(\mathbf x')] = {\rm cov}[\mathbf{\phi(x)}^\top \mathbf w, \mathbf w^\top \mathbf{\phi(x')} ] = \mathbf{\phi(x)}^\top \mathbf S_N \mathbf{\phi(x')} = \beta^{-1} k(\mathbf x, \mathbf x') $$核函数表示两点间的距离远近,两个点如果核函数值大,则相关性高。
我们可以将基函数的计算完全变成核函数的形式,因此,我们可以通过使用等价核,而不显式地引入基函数。
在我们的模型中,可以证明,我们的等价核满足:
$$ \sum_{n=1}^N k(\mathbf x, \mathbf x_n) = 1 $$这里忽略证明的细节,从直觉上,我们如果让所有的 $t_n$ 等于 1,基函数中有一个常数项,然后训练集数目 $N$ 大于独立的基函数数目,那么预测值应该为 1,从而上面的式子成立。
从一般意义上,我们的等价核可以改写为:
$$ k(\mathbf x, \mathbf z) = \mathbf{\psi(x)}^\top \mathbf \psi(\mathbf z) $$其中
$$ \mathbf{\psi(x)} = \beta^{1/2} \mathbf S_N^{1/2} \mathbf{\phi(x)} $$