PRML-3.5 证据近似

3.5证据近似

解决两个超参数$\alpha,\beta$

---

如果我们引入$\alpha, \beta$上的超先验，那么预测分布可以通过边缘化$w,\alpha,\beta$来获得：

$p(t|\textbf{t})=\int\int\int p(t|w,\beta)p(w|\textbf{t},\alpha,\beta)p(\alpha,\beta|\textbf{t})dwd\alpha d\beta \tag{3.74}$

其中$p(t|w,\beta)$由式（3.8）给出，$p(w|\textbf{t},\alpha,\beta)$由式（3.49）给出，其中的$m_N, S_N$分别由（3.53）（3.54）给出。

$$p(t|x,w,\beta) = \mathcal{N}(t|y(x,w), \beta^{-1}) \tag{3.8}$$

$p(w|t) = \mathcal{N}(w|m_N,S_N) \tag{3.49}$

$$\begin{eqnarray} m_N &=& \beta S_N\Phi^Tt \tag{3.53} \\ S_N^{-1} &=& \alpha I + \beta\Phi^T\Phi \tag{3.54} \end{eqnarray}$$

---

第一项和第二项都有了，解决第三项

为了让记号简洁，我们省略了对于输入变量$x$的依赖。如果后验分布$p(\alpha,\beta|\textbf{t})$在$\hat{\alpha},\hat{\beta}$周围是尖峰，那么把$\alpha,\beta$固定为$\hat{\alpha},\hat{\beta}$，然后对$w$求积分：

$因为\alpha,\beta被固定了，所以认为是常数，所以\alpha,\beta两个积分可以拿掉$

$p(t|\textbf{t}) \simeq p(t|\textbf{t},\hat{\alpha},\hat{\beta}) = \int p(t|w,\hat{\beta})p(w|\textbf{t}, \hat{\alpha},\hat{\beta})dw \tag{3.75}$

就得到了预测分布。

---

那么如何才能将$\alpha,\beta$固定住呢?
就要用贝叶斯定理

根据贝叶斯定理，$\alpha,\beta$的后验分布由

$p(\alpha,\beta|\textbf{t}) \propto p(\textbf{t}|\alpha,\beta)p(\alpha,\beta) \tag{3.76}$

给出。如果先验相对较平，那么在证据框架中$\hat{\alpha},\hat{\beta}$可以通过最大化边缘似然函数$p(\textbf{t}|\alpha,\beta)$来获得。我们先得到线性基函数模型的边缘似然函数，然后再找出最大值。这将使我们能够不通过交叉验证，而直接从训练数据确定这些超参数的值。回忆一下比值$\alpha / \beta$类似于正则化参数。

解读一下
$先验相对较平,意思就是近似均匀分布,那么此时的后验最大就是变换似然函数p(\textbf{t}|\alpha,\beta)最大$

此外，值得注意的是，如果在$\alpha,\beta$上定义共轭（Gamma）先验分布，那么对式（3.74）中对这些超参数的边缘化可以通过$w$上的学生t分布（见第2.3.7节）解析的计算出来。虽然得到的结果在$w$上的积分不再有解析解，但是我们可以对这个积分求近似，例如，可以使用基于以后验概率分布的众数为中心的局部高斯近似的拉普拉斯近似方法(见第4.4节)，给证据框架提供了另一种实用的方法（Buntine and Weigend, 1991）。然而，作为关于$w$的被积函数的众数通常是强倾斜的(不准确)，所以拉普拉斯近似方法不能描述概率质量中的大部分信息。这就导致最终的结果要比最大化证据的方法给出的结果差（MacKay, 1999）。

回到证据框架中，我们注意到有两种方法可以用来最大化对数证据。我们可以解析的得到证据函数，然后令它的导数等于零，来得到了$\alpha,\beta$的再估计方程（将在3.5.2节讨论）。另一种方法是，我们使用一种被称为期望最大化（EM）算法的技术，这将在9.3.4节讨论，并在那证明这两种方法会收敛到同一个解。

最大化边缘似然函数的方法
1.假设alpha.beta的先验是Gamma分布，但是这样w就没有解析解了
2.拉普拉斯近似 4.4章节
3.解析计算证据函数，求导=0 本章3.5.2
4.EM算法

3.5.1 计算证据函数

边缘似然函数$p(\textbf{t}|\alpha,\beta)$是通过对权重参数$w$的积分来获得的，即：

$p(\textbf{t}|\alpha,\beta) = \int p(\textbf{t}|w,\beta)p(w|\alpha)dw \tag{3.77}$

计算这个积分的一种方法是再次使用式（2.115）给出的线性-高斯模型的条件分布的结果(之前也常用的对应字符后直接写结果)。这里，我们使用另一种通过对指数项配平方，然后使用高斯分布的标准化系数的基本形式，来计算这个积分。

$$\begin{eqnarray} \ln p(\textbf{t}|w, \beta) &=& \sum\limits_{n=1}^N \ln\mathcal{N}(t_n|w^T\phi(x_n),\beta^{-1}) \ &=& \frac{N}{2}\ln\beta-\frac{N}{2}\ln(2\pi) - \beta E_D(w) \tag{3.11} \end{eqnarray}$$

$$E_D(w) = \frac{1}{2}\sum\limits_{n=1}^N\{t_n - w^T\phi(w_n)\}^2 \tag{3.12}$$

$p(w|\alpha) = \mathcal{N}(w|0,\alpha^{-1}I) \tag{3.52}$

从（3.11）（3.12）和（3.52）中，得到我们可以把证据函数写成

$$p(\textbf{t}|\alpha,\beta) = \left(\frac{\beta}{2\pi}\right)^{N/2}\left(\frac{\alpha}{2\pi}\right)^{M/2}\int exp\{-E(w)\}dw \tag{3.78}$$

的形式。其中$M$是$w$的维数，且定义

$\begin{eqnarray} E(w) &=& \beta E_D(w) + \alpha E_W(w) \ &=& \frac{\beta}{2}\Vert \textbf{t}-\Phi w \Vert^2 + \frac{\alpha}{2}w^Tw \tag{3.79} \end{eqnarray}$

---

3.78 证明-习题3.17
3.11 改形式
$p(\textbf{t}|w, \beta) = \prod\limits_{n=1}{N}N(t_n|w^T\phi(x_n),\beta^{-1})=(\frac{\beta}{2\pi})^{\frac{N}{2}}exp(-\beta E_D(w))$
$其中E_D(w)=\frac{1}{2}||t-\Phi w||^2$
3.52
$p(w|\alpha)=N(w|0,\alpha^{-1}I)=(\frac{\alpha}{2\pi})^{\frac{M}{2}}exp(-\frac{\alpha}{2}w^Tw)$
相乘即可得3.78,3.79

---

(如果忽略一些比例常数，式（3.79）可以看成与正则化的平方和误差函数(3.27)相等。)

$$\frac{1}{2}\sum\limits_{n=1}^N\{t_n-w^T\phi(x_n)\}^2 + \frac{\lambda}{2}w^Tw \tag{3.27}$$

现在，对$w$配平方，可得：

$E(w) = E(m_N) + \frac{1}{2}(w-m_N)^TA(w-m_N) \tag{3.80}$

其中，我们引入了

$A = \alpha I + \beta \Phi^T\Phi \tag{3.81}$

和

$E(m_N) = \frac{\beta}{2}\Vert \textbf{t}-\Phi m_N \Vert^2 + \frac{\alpha}{2}m_N^Tm_N \tag{3.82}$

注意$A$对应的是误差函数的二阶导数矩阵

$A = \nabla\nabla E(w) \tag{3.83}$

这被称为Hessian矩阵(之前证过，可以类比二阶函数)。这里我们定义$m_N$为：

$$m_N = \beta A^{-1}\Phi^T\textbf{t} \tag{3.84}$$

使用式（3.54）(0均值w后验)，可得$A = S_N^{-1}$，因此式（3.84）等价于之前的定义的式（3.53），所以它表示后验分布的均值。

$$\begin{eqnarray} m_N &=& \beta S_N\Phi^Tt \tag{3.53} \\ S_N^{-1} &=& \alpha I + \beta\Phi^T\Phi \tag{3.54} \end{eqnarray}$$

---

3.80证明 - 习题 3.18

---

通过标准化多元高斯分布的标准化系数，就可以简单的得到关于$w$的积分，即：
3.78后面的积分项

$$\begin{array}{l} \int \exp \{-E(\boldsymbol{w})\} \mathrm{d} \boldsymbol{w} \\ =\exp \left\{-E\left(\boldsymbol{m}_{N}\right)\right\} \int \exp \left\{-\frac{1}{2}\left(\boldsymbol{w}-\boldsymbol{m}_{N}\right)^{T} \boldsymbol{A}\left(\boldsymbol{w}-\boldsymbol{m}_{N}\right)\right\} \mathrm{d} \boldsymbol{w} \\ =\exp \left\{-E\left(\boldsymbol{m}_{N}\right)\right\}(2 \pi)^{\frac{M}{2}}|\boldsymbol{A}|^{-\frac{1}{2}} \end{array} \tag{3.85}$$

使用式（3.78）我们可以把对数边缘似然写成
$\ln p(\textbf{t}|\alpha,\beta) = \frac{M}{2}\ln \alpha + \frac{N}{2}\ln \beta - E(m_N) - \frac{1}{2}\ln\vert A \vert - \frac{N}{2}\ln(2\pi) \tag{3.86}$
这就是证据函数所需要的表达式。

---

3.85证明 - 习题 3.19

---

模型证据与多项式阶数之间的关系

代码中的经验贝叶斯公式要和3.5.2结合起来看

点击查看代码

# 3.5
def cubic(x):
    return x * (x - 5) * (x + 5)


x_train, y_train = create_toy_data(cubic, 30, 10, [-5, 5])
x_test = np.linspace(-5, 5, 100)
evidences = []
models = []
for i in range(8):
    feature = PolynomialFeature(degree=i)
    X_train = feature.transform(x_train)
    model = EmpiricalBayesRegression(alpha=100., beta=100.)  # Empirical 经验的意思
    model.fit(X_train, y_train, max_iter=100)
    evidences.append(model.log_evidence(X_train, y_train))
    models.append(model)

degree = np.nanargmax(evidences)
regression = models[degree]  # 9个模型，拿出 evidences 最好的那个

X_test = PolynomialFeature(degree=int(degree)).transform(x_test)
y, y_std = regression.predict(X_test, return_std=True)

plt.scatter(x_train, y_train, s=50, facecolor="none", edgecolor="steelblue", label="observation")
plt.plot(x_test, cubic(x_test), label="x(x-5)(x+5)")
plt.plot(x_test, y, label="prediction")
plt.fill_between(x_test, y - y_std, y + y_std, alpha=0.5, label="std", color="orange")
plt.legend()
plt.show()

plt.plot(evidences)
plt.title("Model evidence")
plt.xlabel("degree")
plt.ylabel("log evidence")
plt.show()

这里，我们已经假定先验分布的形式为式（1.65），参数$\alpha$的值固定为$alpha = 5 × 10^{−3}$。这个图像的形式非常有指导意义。回头看图1.4，得到$M = 0$的多项式对数据的拟合效果非常差，结果模型证据的值也相对较小。$M = 1$的多项式对于数据的拟合效果有了显著的提升，因此模型证据变大了。但是，对于$M = 2$的多项式，因为产生数据的正弦函数是奇函数，所以在多项式展开中没有偶次项，所以拟合效果又变差了。事实上，图1.5给出的数据残差从$M = 1$到$M = 2$只有微小的减小。由于复杂的模型有着更大的复杂度惩罚项，因此从$M = 1$到$M = 2$，模型证据实际上减小了。当$M = 3$时，我们对于数据的拟合效果有了很大的提升，如图1.4所示，因此模型证据再次增大，给出了多项式拟合的最高的模型证据。进一步增加$M$的值，只能少量地提升拟合的效果，但是模型的复杂度却越来越复杂，这导致整体的模型证据下降。再次看图1.5，我们看到泛化错误在$M = 3$到$M = 8$之间几乎为常数，因此只基于这幅图很难对模型做出选择。然而，模型证据的值明显地倾向于选择$M = 3$的模型，因为这是能很好地解释观测数据的最简单的模型。

3.5.2 最大化证据函数

让我们先考虑关于$\alpha$来最大化$p(\textbf{t}|\alpha,\beta)$。首先定义特征方程：

$(\beta\Phi^T\Phi)u_i = \lambda_iu_i \tag{3.87}$

根据式（3.81）得到$A$的特征值$\alpha + \lambda_i$。

$A = \alpha I + \beta \Phi^T\Phi \tag{3.81}$

$\ln p(\textbf{t}|\alpha,\beta) = \frac{M}{2}\ln \alpha + \frac{N}{2}\ln \beta - E(m_N) - \frac{1}{2}\ln\vert A \vert - \frac{N}{2}\ln(2\pi) \tag{3.86}$

$E(m_N) = \frac{\beta}{2}\Vert \textbf{t}-\Phi m_N \Vert^2 + \frac{\alpha}{2}m_N^Tm_N \tag{3.82}$

现在考虑式（3.86）中涉及$\ln \vert A \vert$的项关于$\alpha$的导数。得到：
这里用到了一个线性代数的定理,矩阵的行列式=矩阵的特征值之积

---

$A=P\Lambda P^{-1},相似矩阵$

$det(A) = det(P\Lambda P^{-1})$

$det(A) = det(\Lambda P P^{-1})$

$det(A) = \lambda_1\lambda_2...\lambda_n$

---

$\frac{d}{d\alpha}\ln\vert A \vert = \frac{d}{d\alpha}\ln\prod\limits_i(\lambda_i + \alpha) = \frac{d}{d\alpha}\sum\limits_i\ln(\lambda_i + \alpha) = \sum\limits_i\frac{1}{\lambda_i + \alpha} \tag{3.88}$

因此式（3.86）关于$\alpha$的驻点满足

$0 = \frac{M}{2\alpha} - \frac{1}{2}m_N^Tm_N - \frac{1}{2}\sum\limits_i\frac{1}{\lambda_i + \alpha} \tag{3.89}$

乘以$2\alpha$，整理，可得

$\alpha m_N^Tm_N = M - \alpha\sum\limits_i\frac{1}{\lambda_i+\alpha} = \gamma \tag{3.90}$

由于关于$i$的求和项有$M$个，因此$\gamma$可以写成：

$\gamma = \sum\limits_i\frac{\lambda_i}{\alpha + \lambda_i} \tag{3.91}$

$0 \le \gamma \le M$

稍后讨论$\gamma$的解释。根据方程（3.90），得到使得边缘似然函数最大化的$\alpha$满足：

$\alpha = \frac{\gamma}{m_N^Tm_N} \tag{3.92}$

这里对$\alpha求导有两种方式，书上一种，后面习题3.21还有一种$

---

$$\begin{eqnarray} m_N &=& \beta S_N\Phi^Tt \tag{3.53} \\ S_N^{-1} &=& \alpha I + \beta\Phi^T\Phi \tag{3.54} \end{eqnarray}$$

注意，这是$\alpha$的一个隐式解，不仅仅因为$\gamma$与$\alpha$相关，还因为后验分布众数$m_N$本身也与$\alpha$的选择有关。因此我们使用迭代的方法求解。首先我们选择一个$\alpha$的初始值，再把这个初始值代入式（3.53）求得$m_N$，利用式（3.91）计算得到$\gamma$，然后这些值再被代入式（3.92）来重新估计$\alpha$的值。 不断执行这个过程直至收敛。注意，由于矩阵$\Phi^T\Phi$是固定的，因此我们只需要在最开始的时候计算一次特征值，然后接下来只需乘以$\beta$就可以得到$\lambda_i$的值(3.87特征值的定义)。

需要强调的是，$\alpha$的值是完全通过训练集的数据来确定的。最大似然方法不同，模型的最优的复杂度不需要独立的数据集。

类似地，我们可以关于$\beta$最大化对数边缘似然函数（3.86）。为了达到这个目的，我们注意到由式（3.87）定义的特征值$\lambda_i$正比于$\beta$，即$d\lambda_i/d\beta = \lambda_i/\beta$，得到：

$\frac{d}{d\beta}\ln\vert A\vert = \frac{d}{d\beta}\sum\limits_i\ln(\lambda_i + \alpha) = \frac{1}{\beta}\sum\limits_i\frac{\lambda_i}{\lambda_i+\alpha} = \frac{\gamma}{\beta} \tag{3.93}$

因此，边缘似然的驻点满足：

$$0 = \frac{N}{2\beta}-\frac{1}{2}\sum\limits_{n=1}^N\{t_n - m_N^T\phi(x_n)\}^2 - \frac{\gamma}{2\beta} \tag{3.94}$$

整理可得：

$$\frac{1}{\beta} = \frac{1}{N - \gamma}\sum\limits_{n=1}^N\{t_n - m_N^T\phi(x_n)\}^2 \tag{3.95}$$

与之前一样，这是$\beta$的一个隐式解，可以通过迭代的方法解出。首先选择$\beta$的一个初始值，然后使用这个初始值计算$m_N, \gamma$，然后使用式（3.95）重新估计$\beta$的值，重复直到收敛。如果$\alpha, \beta$的值都是由数据确定的，那么他们的值可以在每次更新$\gamma$之后一起重新估计。

3.5.3 参数的有效数量

对图片做一个解释
在$w_2$方向上,边缘最大似然对应的特征值$\lambda_2$更大(注意这里的椭圆的一长一短两个方向长度对应的是特征值的倒数)，也就是会出现书上的,对于$\lambda_i \gg \alpha$的方向，对应的参数$w_i$会很接近最大似然值，且比值$\lambda_i/(\lambda_i + \alpha)$会很接近1，表示这次的似然学习是有效的，这个参数是有效参数，反之对于$w_1$这个方向，$\lambda_1$更小，表示似然函数学习了老半天，没什么用，$w_{map}$还是更加接近于先验(绿圈)，这个时候就出现了$\gamma$ 接近于0的情况，这个参数是无效参数
所以有一句结论
因此式（3.91）定义的$\gamma$度量了好确定的参数的有效总数。

如果考虑数据点的数量大于参数数量$ N \gg M $的极限情况，那么根据式（3.87），因为$ \Phi^T\Phi $涉及到数据点的隐式求和，因此特征值$ \lambda_i $随着数据集规模的增加而增大，所以所有参数都可以很好的确定。在这种情况下$ \gamma = M $，且$ \alpha,\beta $的再估计方程变成了

$$\begin{eqnarray} \alpha &=& \frac{M}{2E_W(m_N)} \tag{3.98} \\ \beta &=& \frac{N}{2E_D(m_N)} \tag{3.99} \end{eqnarray}$$

其中，$E_W, E_D$分别由（3.25）（3.26）定义。因为这些结果不需要计算Hessian矩阵的一系列特征值，所以可以当成完整证据再估计公式的简化计算的近似。