Book Notes on Pattern Recognition and Machine Learning (PRML)
在之前的内容中，我们提出了一系列的分类回归模型。在最后一节我们考虑结合若干模型以提供更好的预测结果。在具体介绍模型结合的不同方式前，我们首先要区分模型结合与贝叶斯模型平均的区别。
以高斯混合模型为例，我们就结合了若干个高斯分布模型。我们通过二进制隐变量 z 来指示具体模型的选择，所以该模型输出的是 \(p(x,z)\) 联合分布。我们可以通过边缘化隐变量来得到观测数据的分布，假设我们用了 K 个高斯模型
\[p(x) = \sum_z p(x,z) \Rightarrow p(x) = \sum_{k=1}^K \pi_k N(x\vert u_k,\Sigma_k)\]对于独立同分布的数据，每个数据点 \(x_i\) 都有唯一对应的隐变量 \(z_i\) 控制，对于整个数据集 \(X=\{x_1,...,x_n\}\) ，每个数据是由不同的子模型成分构成的，我们可以写出其概率分布
\[p(X) = \prod_{i=1}^n p(x_i) = \prod_{i=1}^n(\sum_{z_i}p(x_i,z_i))\]而在贝叶斯模型平均中，我们有若干个不同下标 \(h=1,...,H\) 对应的模型，比如一个是混合高斯，另一个是神经网络等。我们用 \(p(h)\) 表示具体模型选择的先验概率，整个数据集都是同一种模型生成得到的，即
\[p(X) = \sum_{h=1}^H p(X|h)p(h)\]随着数据量的增大，后验概率 \(p(h\vert X)\) 会越来越趋近于某一个具体模型，这对应整个数据集都是同一种模型生成得到的
最简单的一种组合模型的方法是将组合内的模型输出值进行平均，尤其是在回归任务中。接下来我们来分析为什么这种平均的方式可以提升效果。
我们假设组合模型内总共有 M 个模型，回归的真实函数我们用 \(h(x)\) 表示，第 m 个模型的预测输出我们用 \(y_m(x)\) 表示，那么我们可以得到
\[y_m(x) = h(x) + \epsilon_m(x)\]这里的 \(\epsilon_m(x)\) 表示第 m 个模型下的预测误差。我们采用最小二乘作为损失函数，首先考虑单个模型的 Loss
\[L_m(X) = E_x[\epsilon_m^2(x)]\]如果我们只采用第 m 个模型进行预测，那么这就是对应的 Loss，我们对单个模型预测的 Loss 平均作为统计意义上的单模型 Loss，可以得到
\[L_1(X) = \frac{1}{M}\sum_{m=1}^M E_x[\epsilon_m^2(x)]\]接下来我们考虑组合模型的输出，将所有模型的预测进行平均作为最终预测结果，此时的最小二乘损失函数写为
\[L_2(X) = E_x[\{\frac{1}{M}\sum_{m=1}^My_m(x) -h(x)\}^2] = E_x [\{\frac{1}{M}\sum_{m=1}^M \epsilon_m(x)\}^2]\]如果我们假设误差的均值为零，\(E_x[\epsilon_m(x)] = 0\) 且不同模型之间的误差是不相关的 \(E_x[\epsilon_i(x)\epsilon_j(x)] = 0\) ，那么可以得到
\[L_2(X) = \frac{1}{M} L_1(X)\]也就是说，组合模型的误差只有单个模型误差的 1/M 。但实际中，不同模型间的误差是高度相关的，我们之前的假设太严格了。即使是这样，我们依然有 \(L_2\leq L_1\) ，保证了组合模型的误差不大于其中的单个模型
Boosting 是一种结合多个分类器以实现更好效果的方法，也可以扩展到回归任务。这里我们介绍最常用的 adaptive boosting，简称为 AdaBoost 。以二分类问题为例，假设数据集数据为 \(\{(x_i,y_i^\star)\vert i=1,...,n\}\) 并且 \(y_i^\star \in \{+1,-1\}\)
在 AdaBoost 中，我们为每一个训练数据分配权重，初始时所有数据的权重相同。在训练过程中，我们串行地训练组合模型中的子模型。根据上一个子模型的分类预测结果，修改数据权重。对于正确预测的数据，权重减小；分类错误的数据，其权重或增加。这种权重会传递至下一个子模型，并继续被修改。
而我们最终的预测，是基于所有的子模型预测结果，结合每个模型之前对应的数据权重，通过加权平均的方式输出最终的分类预测。接下来我们给出具体的算法流程以更好地说明 AdaBoost。和之前假设一样，我们假设一共有 M 个子模型，模型参数用 \(w\) 表示
第一步，初始化。将第一个模型的权重全部初始化为 1/n ，即
\[w_i^{(1)} = 1/n \quad for\ i = 1,2,...,n\]第二步，串行训练组合模型中的子模型，从第 1 个一直训练到第 M 个模型。我们假设当前正在训练第 m 个模型。
训练第 m 个模型，我们用 \(y_m(x)\) 表示其预测输出。采用带权重的损失函数，即
\[L_m = \sum_{i=1}^n \frac{w^{(m)}_i I(y_m(x_i)\neq y_i^*)}{\sum_{j=1}^n w^{(m)}_j}\]统计该模型的错误分类率，计算方式为
\[\epsilon_m = \frac{\sum_{i=1}^n w_i^{(m)} I(y_m(x_i)\neq y_i^*)}{\sum_{i=1}^n w_i^{(m)}}\]如果错误分类率 \(\epsilon_m\) 大于0.5，说明此时模型表现甚至不如随机猜测，可以直接结束算法流程。如果小于0.5，则用这个错误分类率 \(\epsilon_m\) 来计算权重系数 \(\alpha\) ，即
\[\alpha_m = \frac{1}{2}\ln \frac{1-\epsilon_m}{\epsilon_m}\]更新数据的权重 w ，并传至下一个(第 m+1 个子模型)模型
\[w_i^{(m+1)} = w_i^{(m)} e^{-y_i^*\alpha_iy_m(x_i)}\]最后当所有子模型都训练完毕后，结合权重系数加权平均，给出 AdaBoost 的最终预测结果，即
\[y_{pred}(x) = sign(\sum_{m=1}^M\alpha_m y_m(x))\]对于最终的预测结果 \(y_{pred}\) ，我们同样需要损失函数来计算更新。这里我们考虑使用指数损失函数，我们假设选用 m 个子模型组合作为最终的预测，记为 \(f_m(x)\) ，可以写出损失函数
\[L = \sum_{i=1}^n\frac{ w_i^{(m)}exp\{-y_i^*f_m(x_i)\}}{\sum_{j=1}^n w_j^{(m)}}, \quad where\ f_m(x) = \sum_{i=1}^m \alpha_i y_i(x)\]我们首先来看第一个问题，为什么这里选用指数损失函数，用这种损失函数是否是正确的？
如果最终的预测模型 \(f_m(x)\) 可以让指数损失函数最小，我们令其偏导为零
\[\frac{\partial L}{\partial f_m} = -e^{-f_m(x)}P(y^*=1|x) + e^{f_m(x)}P(y^*=-1|x) = 0\]可以得到
\[f_m(x) = \frac{1}{2}\ln \frac{P(y^*=1|x)}{P(y^*=-1|x)} \Rightarrow sign(f_m(x)) = argmax_y P(y^*=y|x)\]可以看出，\(sign(f_m(x))\) 达到了贝叶斯最优错误率。也就是说，当指数损失函数最小时，分类错误率也将达到最小。我们选用指数损失函数并不会对该任务本身造成影响。
接下来我们分析 AdaBoost 算法中，模型 \(y_m\) 更新公式及损失函数的由来。同样在这里，我们只考虑第 m 个子模型，假设之前的模型已经训练好了，将指数损失函数中的 \(f_m(x)\) 进行拆分
\[L(y_m) = \sum_{i=1}^n \frac{w^{(m)}_i e^{-y_i^*f_{m-1}(x_i)}e^{-y_i^*y_m(x_i)}}{\sum_{j=1}^n w^{(m)}_j}\]因为是二分类问题，\((y_i^\star)^2 = y_m(x_i)^2 = 1\) ，我们可以对上式中的后一项 \(e^{-y_i^\star y_m(x_i)}\) 进行泰勒展开得到
\[e^{-y_i^*y_m(x_i)} \approx 1 - y_i^*y_m(x_i) + 1/2\]代入上式可以发现，对于理想的子模型 \(y_m(x)\) ，应该满足
\[y_m(x) = argmin\ L(y_m) = argmax\ \sum_{i=1}^n \frac{w^{(m)}_i y_i^*y_m(x_i)}{\sum_{j=1}^n w^{(m)}_j}\]我们利用正负一的二分类条件，可以得到如下关系
\[y_i^*y_m(x_i) = 1-2I(y_m(x_i)\neq y_i^*)\]所以我们就可以得到之前算法中的损失函数形式，即
\[y_m(x) = argmin\ \sum_{i=1}^n \frac{w^{(m)}_i I(y_m(x_i)\neq y_i^*)}{\sum_{j=1}^n w^{(m)}_j} = argmin\ L_m\]接下来我们分析 AdaBoost 算法中，模型权重 \(\alpha_m\) 更新公式的由来，此时我们已经训练好了对应的模型 \(y_m(x)\)。当我们训练好第 m 个模型 \(y_m(x)\) 时，该模型权重 \(\alpha_m\) 的选取应当使指数损失函数最小。假设在最后的优化中，前 m-1 个子模型及其权重参数 \(\alpha\) 是固定的。只有第 m 个模型 \(y_m(x)\) 和对应的权重 \(\alpha_m\) 可以优化。把最后的第 m 个子模型从损失函数中单独拆分出来，可以写成
\[L = \sum_{i=1}^n \frac{w^{(m)}_iexp\{-y_i^*f_{m-1}(x_i) -y_i^*\alpha_my_m(x_i)\}}{\sum_{j=1}^n w^{(m)}_j} = L(\alpha_my_m) \times constant\]进一步改写 \(L(\alpha_my_m)\) 可以得到
\[L(\alpha_my_m) = \sum_{i=1}^n\frac{ w_i^{(m)}exp\{-y_i^* \alpha_m y_m(x_i)\}}{\sum_{j=1}^n w_j^{(m)}} = e^{-\alpha_m}(1-\epsilon_m) + e^{\alpha_m}\epsilon_m\]将其对 \(\alpha_m\) 求导为零就可以得到之前算法中的模型权重更新公式，即
\[\frac{\partial L(\alpha_my_m)}{\partial \alpha_m} = 0\Rightarrow \alpha_m = \frac{1}{2}\ln \frac{1-\epsilon_m}{\epsilon_m}\]最后，我们来分析 AdaBoost 算法中，数据权重 w 更新公式的由来。假设我们已经完成了前 m 个子模型的训练，需要以此修改传递到第 (m+1) 个模型的数据权重。根据归一化，我们写出针对数据 i 的权重更新公式
\[w_i^{(m+1)} = \frac{e^{-y_i^*f_m(x_i)}}{\sum_{j=1}^n w_j^{(m)} e^{-y_j^*f_m(x_j)}} = \frac{e^{-y_i^*f_{m-1}(x_i)} e^{-y_i^*\alpha_iy_m(x_i)}}{\sum_{j=1}^n w_j^{(m)} e^{-y_j^*f_m(x_j)}}\]我们在分子分母同乘针对前 m-1 项的和，以得到 \(w_i^{(m)}\) 的形式，代入后得到之前算法中的数据权重更新公式，即
\[w_i^{(m+1)} = w_i^{(m)} e^{-y_i^*\alpha_iy_m(x_i)} \times constant\]至此我们就推导出了算法中所有更新公式的来源。最后我们来比较一下这里的指数误差函数和之前的交叉熵函数的优劣。
相对于 Boosting 的串行集成学习，Bagging 是并行式集成学习的一种。这种方法主要是基于自主采样法 (bootstrap sampling)。因此，在介绍 Bagging 之前，我们先介绍这种采样方法。
bootstrap sampling 的采样流程是这样的，假如我们已知包含 n 个样本的数据集，想要获得包含 n 个样本的采样集。我们每次从数据集中随机选择一个样本，将其复制放入采样集中，如此重复 n 次即可得到采样集。
从这个流程中也可以看出，我们的采样集中很可能包含相同的样本，有些数据集中的数据也可能不会出现在采样集中。每次采样，每个数据被采中的概率为 1/n。我们可以通过极限逼近，对于某一个数据，其不会出现在采样集中的概率为
\[\lim_{n\rightarrow 0}\ (1-\frac{1}{n})^n = \frac{1}{e} \approx 0.368\]根据这个概率，结合大量数据的统计意义。我们可以知道，初始数据集中约有 \(63.2\%\) 的数据出现在采样集中。
接下来对 Bagging 进行介绍。我们假设总共有 m 个子模型。Bagging 的基本流程就是通过 m 次 bootstrap sampling 得到 m 个样本数为 n 的采样集，分别用来训练这 m 个子模型，最后将这 m 个子模型的输出进行结合作为最终的预测输出。
对于不同的任务，可以选取不同的模型输出结合策略。
针对回归问题，我们可以简单采用平均法，将不同子模型的输出平均作为最终输出，即
\[y_{pred}(x) = \frac{1}{m}\sum_{i=1}^m y_i(x)\]也可以采用加权平均法，即
\[y_{pred}(x) = \frac{1}{m}\sum_{i=1}^m w_i y_i(x)\]这里的权重一般是从训练数据中学习得到。一般而言，在子模型性能相差大时宜使用加权平 均，子模型性能接近时宜使用简单平均
针对分类问题，主要采用的是投票法。即针对每个数据，选取不同子模型分类预测中，出现次数最高的那个类别作为最终的预测分类。
最后，我们介绍一种通过一个单独的子模型结合其他子模型的模型混合方法，其中的典型代表是 Stacking 。这种方法的原理很简单，将不同子模型的输出作为这个单独子模型（也称为次级模型，meta-learner）的输入，次级模型的输出即为整个混合模型的输出。
这里存在一个问题，次级模型的训练集是利用子模型产生的，若直接用子模型的训练集来产生次级训练集，则过拟合风险会比较大。因此，一般是通过使用交叉验证或留一法这样的方式，用训练子模型时未使用的样本来产生次级模型的训练样本。
有研究表明，将子模型的输出类概率作为次级模型的输入属性，用多响应线性回归(Multi-response Linear Regression，简称 MLR) 作为次级模型效果较好。
接下来我们考虑条件分布的混合模型。最简单的情况下，混合系数与输入变量相互独立，我们将介绍线性回归混合模型和 logistic 回归混合模型。如果混合系数依赖于输入变量，我们可以得到 mixture of experts 模型。更进一步，如果我们让混合模型中的每一个子模型同样是 mixture of experts 模型，那么可以得到层次 mixture of experts 模型
我们假设混合 K 个回归模型，每一个回归模型用参数 \(w_k\) 表示。混合系数我们用 \(\pi\) 表示。同时，假设这 K 个模型下的噪声方差由相同的参数 \(\beta\) 控制。现在我们可以写出混合分布
\[p(y|\theta) = \sum_{k=1}^K \pi_k N(y\vert w_k^Tx,\beta^{-1})\]这里的 \(\theta\) 代表所有可调整参数 \(\{w,\pi,\beta\}\) 。在给定数据集 \(\{(x_i,y_i^*)\vert i=1,...,n\}\) 的情况下，可以写出 log 似然方程
\[\ln p(Y^*\vert \theta) = \sum_{i=1}^n \ln (\sum_{k=1}^K \pi_kN(y_i^*\vert w_k^Tx_i,\beta^{-1}))\]为了最大化这个似然方程，我们可以引入二进制隐变量 z 结合 EM 算法进行求解。混合系数 \(\pi\) 的求法和之前介绍的高斯混合模型完全相同，不同的是回归模型 \(w\) 和噪声方差 \(\beta\) 的更新式。
我们首先来看回归模型的求解，我们针对第 k 个模型 \(w_k\) ，将 EM 的目标函数 \(Q(\theta,\theta^{t})\) 中关于回归模型的部分拆分出来可以得到
\[Q(\theta,\theta^{t}) = \sum_{i=1}^n \gamma_{ik}\{-\frac{\beta}{2}(y_i^*-w_k^Tx_i)^2\} + constant\]接着对回归模型参数求导并令其为零，消去常数系数 \(\beta\) 可得
\[\sum_{i=1}^n \gamma_{ik}(y_i^*-w_k^Tx_i)x_i = 0\]和单纯的回归模型求解不同的是，这里多了一个数据权重 \(\gamma_{ik}\) 。这个权重表示第 i 个数据在第 k 个回归模型中的权重，会随着 EM 算法的迭代而更新。在带权参数的影响下，模型的解也会发生变化，此时的解析解可以写为
\[w_k = (X^TR_kX)^{-1}X^TR_kY^* ,\quad where\ \ R_k = diag(\gamma_{nk})\]和没有权参数的单纯线性回归模型的解比较一下，\(w_k = (X^TX)^{-1}X^TY^*\) ，也就是多了对权重的考虑。
最后看下对噪声方差 \(\beta\) 的优化。同样的方法，我们将 EM 的目标函数 \(Q(\theta,\theta^{t})\) 中关于 \(\beta\) 的部分拆分出来
\[Q(\theta,\theta^{t}) = \sum_{i=1}^n\sum_{k=1}^K \gamma_{ik}\{\frac{\ln \beta}{2}-\frac{\beta}{2}(y_i^*-w_k^Tx_i)^2\}\]因为 \(\beta\) 是一个标量，我们可以直接对其求导令导数为零得到 EM 中的更新公式
\[\frac{1}{\beta} = \frac{1}{n}\sum_{i=1}^n\sum_{k=1}^K \gamma_{ik}(y_i^*-w_k^Tx_i)^2\]我们求出了 EM 算法中 \(\{\pi,w,\beta\}\) 的更新公式，最后按照算法流程更新即可。
我们假设混合 K 个回归模型，每一个回归模型用参数 \(w_k\) 表示，\(y_k(x) = \sigma(w_k^Tx)\)。混合系数我们用 \(\pi\) 表示。同时，假设是 0-1 二分类问题。现在我们可以写出混合分布
\[p(y^*\vert x,\theta) = \sum_{k=1}^K \pi_k y_k^{y^*}[1-y_k]^{1-y^*}\]这里的 \(\theta\) 代表所有可调整参数 \(\{w,\pi\}\) 。在给定数据集 \(\{(x_i,y_i^*)\vert i=1,...,n\}\) 的情况下，可以写出似然方程
\[p(Y^*\vert \theta) = \prod_{i=1}^n(\sum_{k=1}^K \pi_k y_{ik}^{y^*_i}[1-y_{ik}]^{1-y^*_i})\]和上面一样的方法，我们引入隐变量 Z 结合 EM 算法最大化似然函数，此时目标方程为
\[Q(\theta,\theta^t) = E_z[\ln p(Y^*Z\vert\theta)] = \sum_{i=1}^n\sum_{k=1}^K\gamma_{ik}\{\ln \pi_k + y_i^*(\ln y_{ik}) + (1-y_i^*)\ln (1-y_{ik})\}\]这里混合系数 \(\pi\) 的更新和之前一致，但 logistic 回归模型 \(w\) 没有解析解，必须通过其他的迭代算法去逼近近似。这种混合模型可以扩展到多分类任务，只需将二分类的 sigmoid 函数替换为 softmax 函数
最后简要介绍一下 Mixtures of experts 模型，在这个模型中，我们认为混合系数也是和输入变量相关的。和之前的两个混合模型相比，变化为 \(\pi_k \rightarrow \pi_k(x)\)
此时的条件概率分布可以写为
\[p(Y^*\vert x) = \sum_{k=1}^K \pi_k(x) p_k(Y^*\vert x)\]这里的混合系数 \(\pi_k(x)\) 被称为 gating 函数，它需要满足和之前一样的约束条件，即
\[0\leq \pi_k(x) \leq 1,\quad \sum_{k=1}^K \pi_k(x) =1\]而这里的子概率密度 \(p_k(y^*\vert x)\) 被称为 experts。
Mixtures of experts 和我们之前在神经网络中提到的 mixture density network 关联很密切。前者的优点是可以通过 EM 算法结合凸优化求解，后者的优点是子模型与混合系数共享神经网络的隐藏单元，对输入空间划分约束更少