《机器学习》

简介

苏老师的《机器学习》。一开始的内容和吴恩达的有些相似，就记一些重要的推导；后面的主要来源于《统计学习方法》和PRML，就记得比较多。

曲线拟合（线性回归）

$$\Phi(x)=\left [\begin{aligned}&\phi^T(x_1)\\&\phi^T(x_2)\\ &\dots\\&\phi^T(x_n)\end{aligned}\right]=\left [\begin{aligned}&1 &x_1 &\dots &x_1^{d-1}\\&1 &x_2 &\dots &x_2^{d-1}\\&&&\dots\\&1 &x_n &\dots &x_n^{d-1}\end{aligned}\right],y=\left [\begin{aligned}y_1\\y_2\\ \dots\\y_n\end{aligned}\right]$$

经验损失：$D(\mathcal{L})=\frac{1}{n}\sum_{i=1}^n\mathcal{L}(y_i,\hat y_i)$

期望损失：$E(\mathcal{L})=\int\mathcal{L}(x,y)p(x,y)dxdy$

设：$h=\Phi w-y$，其中 $x\in\mathbb R^{n\times 1},\phi(x_1)\in\mathbb R^{d\times 1},\Phi\in\mathbb R^{n\times d}, w \in \mathbb R^{d\times 1}, y\in \mathbb R^{n\times 1}$，最小二乘法为，

$$\begin{align*} L_D(w)&=\frac{1}{2}h^Th\\ &=\frac{1}{2}(\Phi w-y)^T(\Phi w-y)\text{ ———展开}\\ &=\frac{1}{2}(w^T\Phi^T\Phi w-w^T\Phi^Ty-y^T\Phi w+y^Ty)\text{ ———}w^T\Phi^Ty=y\Phi w\\ &=\frac{1}{2}(w^T\Phi^T\Phi w-2w^T\Phi^Ty+y^Ty) \end{align*}$$

根据 $\nabla_x(x^TAx)=(A+A^T)x$，对 $w$ 求偏导：

$$\nabla_w L_D(w)=\Phi^T\Phi w-\Phi^Ty=0\\ w=(\Phi^T\Phi)^{-1}\Phi^Ty$$

不同角度看解析解

核

$$\Phi^Ty=\left(\Phi_1,\Phi_2,...,\Phi_n \right)\left(\begin{aligned}y_1\\y_2\\...\\y_n\end{aligned}\right)=\sum_{i=1}^ny_i\phi(x_i)$$

新来一个样本 $\forall x\in \mathbb{R}^d$，则

$$\begin{align*} \hat y&=\phi^T(x)\sum_{i=1}^n(\Phi^T\Phi)^{-1}y_i\phi(x_i)\\ &=\sum_{i=1}^n y_i\phi^T(x)(\Phi^T\Phi)^{-1}\phi(x_i)\\ &=\sum_{i=1}^ny_ik(x,x_i) \end{align*}$$

$k(x,x_i)$ 即为一种核函数。

空间

真实值：$y \in \mathbb{R}^n$，预测值：$\hat y =\Phi w\in \mathbb{R}^n$，则虚线部分为真实值与预测值的误差 $\Phi w-y$。且虚线部分与 $\Phi$ 所构成的空间的每个分量 $\phi_i$ 都垂直，即

$$\left \{\begin{aligned}\phi_1^T(\Phi w-y)=0\\\phi_2^T(\Phi w-y)=0\\\dots\\\phi_d^T(\Phi w-y)=0 \end{aligned}\right.$$

即，$\Phi^T(\Phi w-y)=0$。

最大似然估计

可以将每个样本和标签看成计算出的参数 $w$，将特征空间投影到标签空间后，再附上一个随机变量，

$$y_i=w^Tx_i+\varepsilon,\varepsilon \sim N(0,1)$$

则 $D=\{\langle x_i,y_i \rangle\}_{i=1}^n$。因为 $\varepsilon$ 属于高斯分布，则 $y_i$ 也属于高斯分布，$y_i \sim N(w^Tx_i,1)$。将其展开：

$$p(y_i|x_i;w)=\frac{1}{\sqrt {2\pi}}\text{exp}(-\frac{(y_i-w^Tx_i)^2}{2})$$

则似然函数为：

$$\begin{align*} \mathcal{l}(w)&=P(Y|X;w)\\ &=\prod_{i=1}^np(y_i;x_i,w)\\ &=\prod_{i=1}^n\frac{1}{\sqrt {2\pi}}\text{exp}(-\frac{(y_i-w^Tx_i)^2}{2}) \end{align*}$$

对数似然函数为：

$$\begin{align*} L(w)&=\log l(w)\\ &=\text{log}\prod_{i=1}^n\frac{1}{\sqrt {2\pi}}\text{exp}(-\frac{(y_i-w^Tx_i)^2}{2})\\ &=-n\text{log}\sqrt{2\pi}-\frac{1}{2}\sum_{i=1}^n(w^Tx_i-y_i)^2 \end{align*}$$

其中第一项为定值，第二项为最小二乘法的损失函数。

最大后验概率，引入正则化

$D$ 为数据集，则 $p(w|D)\propto p(w,D)\propto p(w)p(D|w)$。

共轭分布定义：$p(w)$ 与 $p(w|D)$ 同分布。

假设 $w \sim N(\vec 0,\lambda I)$，即协方差 $\Sigma=\lambda I$。

$$\begin{align*} p(w)&=\frac{1}{\sqrt{(2\pi)^d|\Sigma|}}\text{exp}(-\frac{1}{2}w^T\Sigma^{-1}w)\\ &=\frac{1}{\sqrt{(2\pi)^d\lambda}}\text{exp}(-\frac{1}{2\lambda}w^Tw) \end{align*}$$

似然对数为

$$\begin{align*} \text{log }p(w)=\text{log }\frac{1}{\sqrt{2\pi\lambda}}-\frac{1}{2\lambda}w^Tw \end{align*}$$

同样地，第一项为定值，第二项为正则化。

感知机（线性分类）

这里讲的感知机是单层的，用来作分类。

不像多层神经网络主要靠的是多层的非线性作用进行回归和分类。

所以，这里单纯讲了两种用于分类的激活函数。

问题背景：训练感知机进行二分类。

[-1,1]的阶跃函数

激活函数：

$$y=\text{sgn}(w^Tx+b)$$

损失函数：

$$L(w,b)=-\sum_{x^{(i)}\in M}y^{(i)}(w^Tx^{(i)}+b)=\sum_{i=1}^n\text{max}(0,-y^{(i)}(w^Tx^{(i)}+b))$$

其中， $M$ 是指分错的类别。

补充：如果改为SVM（支持向量机）进行分类的话，则损失函数变为：

$$L(w,b)=-\sum_{x^{(i)}\in M}y^{(i)}(w^Tx^{(i)}+b)=\sum_{i=1}^n\text{max}(0,1-y^{(i)}(w^Tx^{(i)}+b))$$

这是因为SVM多限定了一个最大距离。

对损失函数求导：

$$\frac{\partial L}{\partial w}=-y^{(i)}x^{(i)}\\ \frac{\partial L}{\partial b}=-y^{(i)}$$

所以随机下降的更新规则为：

$$w\leftarrow w-(-\eta y^{(i)}x^{(i)})=w+\eta y^{(i)}x^{(i)}\\ b\leftarrow b-(-\eta y^{(i)})=b+\eta y^{(i)}$$

Logistic函数

$$h(x)=\frac{1}{1+e^{-(w^Tx+b)}}$$

假设：

$$\begin{align*} &P(y=1|x;w)=h(x)\\ &P(y=0|x;w)=1-h(x) \end{align*}$$

以上为样本 $x$ 被分为某一类分别的概率，则 $x$ 被分到真实标签的概率为：

$$p(y|x;w)=(h(x))^y(1-h(x))^{1-y}\\$$

则似然函数为：

$$\begin{align*} L(w)&=p(\vec y|X;w) \\ &=\prod_{i=1}^np(y^{(i)}|x^{(i)};w) \\ &=\prod_{i=1}^n(h(x^{(i)}))^{y^{(i)}}(1-h(x^{(i)}))^{1-y^{(i)}} \end{align*}$$

似然函数越大越好，而损失函数越小越好，所以我们将似然函数取负再取对数，

$$\begin{align*} -l(w)&=-\log \space L(w) \\ &=-\sum_{i=1}^n\left[ y^{(i)}\log \space h(x^{(i)})+(1-y^{(i)})\log(1-h(x^{(i)}))\right] \end{align*}$$

求导：

$$\begin{align*} -\frac{\partial l(w)}{\partial w_j} &= -(y\frac{1}{h(x)}-(1-y)\frac{1}{1-h(x)})\frac{\partial}{\partial w_j}h(x) \\ &= -(y\frac{1}{h(x)}-(1-y)\frac{1}{1-h(x)})h(x)(1-h(x))\frac{\partial w^Tx}{\partial w_j}\\ &=-(y(1-h(x))-(1-y)h(x))x_j \\ &=-(y-h(x))x_j \end{align*}$$

所以，对于单个样本来说，

$$-\frac{\partial l(w)}{\partial w}=-(y-h(x))x$$

所以随机下降的更新规则为：

$$w\leftarrow w-\eta[-(y^{(i)}-h(x^{(i)}))x^{(i)}]=w+\eta y^{(i)}x^{(i)}-\eta h(x^{(i)})x^{(i)}$$

支持向量机SVM

基本概念

支持向量机（support vector machine，SVM）是一种二类分类模型。它的基本模型是定义在特征空间上的间隔最大的线性分类器，间隔最大使它有别于感知机。支持向量机的学习策略就是间隔最大化，可形式化为一个求解凸二次规划的问题，也等价于正则化的合页损失函数的最小化问题。支持向量机的学习算法是求解凸二次规划的最优化算法。

线性可分

如果一个线性函数能够将样本分开，称这些数据样本是线性可分的。那么什么是线性函数呢？其实很简单，在二维空间中就是一条直线，在三维空间中就是一个平面，以此类推，如果不考虑空间维数，这样的线性函数统称为超平面。我们看一个简单的二维空间的例子，O代表正类，X代表负类，样本是线性可分的，但是很显然不只有这一条直线可以将样本分开，而是有无数条，我们所说的线性可分支持向量机就对应着能将数据正确划分并且间隔最大的直线。

间隔最大

一般来说，一个点距离分离超平面的远近可以表示分类预测的确信度，如图中的 $A, B$ 两个样本点，$B$ 点被预测为正类的确信度要大于 $A$ 点，所以SVM的目标是寻找一个超平面，使得离超平面较近的异类点之间能有更大的间隔，即不必考虑所有样本点，只需让求得的超平面使得离它近的点间隔最大。

函数间隔与几何间隔

给定一个训练样本 $(x_i,y_i)$，$x$ 是特征，$y$ 是标签，$i$ 表示第 $i$ 个样本。定义函数间隔为：

$$\hat \gamma_i=y_i(w·x_i+b)$$

函数间隔可以表示分类预测的正确性及确信度，但是选择分离超平面时，只有函数间隔还不够。因为只要成比例的改变 $w$ 和 $b$ ，例如将他们改为 $2w$ 和 $2b$ ，超平面并没有改变，但函数间隔却成为原来的2倍。

这一事实告诉我们，可以对分离超平面的法向量 $w$ 加某些约束，如规范化，$||w||=1$，使得间隔是确定的。这时候函数间隔就成为了几何间隔。

$$\gamma_i=y_i\left(\frac{w}{||w||}·x_i+\frac{b}{||w||}\right)$$

可以把几何间隔理解为函数间隔的归一化形式。

间隔最大化

间隔最大化的直观解释是：对训练数据集找到几何间隔最大的超平面意味着以充分大的确信度对训练数据进行分类。也就说，不仅将正负实例分开，而且对最难分的实例点（离超平面最近的点）也有足够大的确信度将它们分开。这样的超平面应该对未知的新实例又很大的分类预测能力。

这个问题可以表示为下面的约束最优化问题：

$$\max_{w,b}\left(\min_{i=1,2,\dots,n}\gamma_i\right)\Leftrightarrow\left\{\begin{aligned}&\max_{w,b}\gamma\\&\text{s.t. }y_i\left(\frac{w^Tx_i+b}{||w||}\right)\ge \gamma \\\end{aligned}\right.$$

即我们希望最大化超平面 $(w,b)$ 关于训练数据集的几何间隔 $\gamma$。约束条件表示的是超平面 $(w,b)$ 关于每个训练样本的几何间隔至少是 $\gamma$。

根据几何间隔和函数间隔的关系是，可以改写为，

$$\left\{\begin{aligned}&\max_{w,b}\gamma\\&\text{s.t. }y_i\left(\frac{w^Tx_i+b}{||w||}\right)\ge \gamma \\\end{aligned}\right.\Leftrightarrow\left\{\begin{aligned}&\max_{w,b}\frac{\hat \gamma}{||w||}\\&\text{s.t. }y_i\left(\frac{w^Tx_i+b}{||w||}\right)\ge \frac{\hat \gamma}{||w||}\\\end{aligned}\right.\Leftrightarrow\left\{\begin{aligned}&\max_{w,b}\frac{\hat \gamma}{||w||}\\&\text{s.t. }y_i\left({w^Tx_i+b}\right)\ge \hat \gamma \\\end{aligned}\right.$$

函数间隔 $\hat \gamma$ 的取值并不影响最优化问题的解。这样，就可以取 $\hat \gamma=1$。将 $\hat \gamma=1$ 代入上面的最优化问题。注意到最大化 $\frac{1}{||w||}$ 和最小化 $\frac{1}{2}||w||^2$ 是等价的，所以我们可以改写为

$$\left\{\begin{aligned}&\max_{w,b}\frac{\hat \gamma}{||w||}\\&\text{s.t. }y_i\left({w^Tx_i+b}\right)\ge \hat \gamma \\\end{aligned}\right.\Leftrightarrow\left\{\begin{aligned}&\max_{w,b}\frac{1}{||w||}\\&\text{s.t. }y_i\left({w^Tx_i+b}\right)\ge 1 \\\end{aligned}\right.\Leftrightarrow\left\{\begin{aligned}&\min_{w,b}\frac{1}{2}||w||^2\\&\text{s.t. }y_i\left({w^Tx_i+b}\right)-1\ge 0 \\\end{aligned}\right.$$

这就是一个凸二次规划问题。

硬间隔最大化

对偶问题

原始问题为：

$$\min_{x\in \mathbb R^n}f(x)\\ \begin{aligned} \text{s.t. }&c_i\le 0,i=1,\dots,k\\ &h_j=0,j=1,\dots,l \end{aligned}$$

拉格朗日函数为：

$$L(x,\alpha,\beta)=f(x)+\sum_{i=1}^k\alpha_ic_i(x)+\sum_{j=1}^l\beta_jh_j(x),\ \alpha_i\ge0$$

考虑 $x$ 的函数：

$$\theta_p(x)=\max_{\alpha_i,\beta_j:\alpha_i\ge0}L(x,\alpha,\beta)$$

其中，$p$ 表示原问题。

1. 如果不满足条件，即 $c_i(x)>0$ 或 $h_j(x)\neq 0$，则：

   $$\theta_p(x)=\max_{\alpha_i,\beta_j:\alpha_i\ge0}\left[f(x)+\sum_{i=1}^k\alpha_ic_i(x)+\sum_{j=1}^l\beta_jh_j(x)\right]=+\infty$$

2. 如果 $x$ 满足约束条件，则：

   $$\theta_p(x)=f(x)$$

所以，

$$\theta_p(x)=\left\{\begin{aligned}&f(x)&&,x满足原始问题约束\\&+\infty&&,其他\end{aligned}\right.$$

考虑极小化的问题：

$$\min_x\theta_p(x)=\min_x\max_{\alpha_i,\beta_j:\alpha_i\ge0}L(x,\alpha,\beta)$$

它和原问题是等价的，即它们有相同的解，称为广义拉格朗日的极小极大问题。即在满足约束的条件下，通过 $x$ 使问题最小化。

定义原始问题的最优值：

$$p^\ast=\min_x\theta_p(x)$$

对偶问题定义为：

$$\theta_d(\alpha,\beta)=\min_xL(x,\alpha,\beta)$$

极大化上述问题：

$$\max_{\alpha_i,\beta_j:\alpha_i\ge 0}\theta_d(\alpha,\beta)=\max_{\alpha_i,\beta_j:\alpha\ge0}\min_xL(x,\alpha,\beta)$$

上述问题称为广义拉格朗日函数的极大极小问题。

定义对偶问题的最优值：

$$d^\ast=\max_{\alpha_i,\beta_j:\alpha_i\ge0}\theta_d(\alpha,\beta)$$

原始问题与对偶问题的关系：对任意的 $\alpha$，$\beta$ 和 $x$，有：

$$\theta_d(\alpha,\beta)=\min_xL(x,\alpha,\beta)\le L(x,\alpha,\beta)\le \max_{\alpha_i,\beta_j:\alpha_i\ge0}L(x,\alpha,\beta)=\theta_p(x)$$

即，

$$\theta_d(\alpha,\beta)\le \theta_p(x)$$

即，

$$\max_{\alpha_i,\beta_j:\alpha\ge0}\theta_d(\alpha,\beta)\le\min_x\theta_p(x)$$

即，

$$d^\ast=\max_{\alpha_i,\beta_j:\alpha_i\ge 0}\min_x L(x,\alpha,\beta)\le\min_x\max_{\alpha_i,\beta_j:\alpha_i\ge0}L(x,\alpha,\beta)=p^\ast$$

KKT：假设 $L(x,\alpha,\beta)$ 和 $c_i(x)$ 是凸函数，$h_j(x)$ 是仿射函数，且约束 $c_i(x)$ 严格可执行，则满足以下条件后，$x^\ast$，$\alpha^\ast$ 和 $\beta^\ast$ 为极值解。

$$\left\{\begin{aligned} &\nabla_xL(x^\ast,\alpha^\ast,\beta^\ast)=0\\ &\nabla_\alpha L(x^\ast,\alpha^\ast,\beta^\ast)=0\\ &\nabla_\beta L(x^\ast,\alpha^\ast,\beta^\ast)=0\\ &\alpha_i^\ast c_i(x^\ast)=0,i=1,\dots,k\\ &c_i(x^\ast)\le0,i=1,\dots,k\\ &\alpha_i^\ast\ge 0,i=1,\dots,k\\ &h_j(x^\ast)=0,i=1,\dots,l \end{aligned}\right.$$

其中，$\alpha_i^\ast c_i(x^\ast)=0$ 称为KKT的对偶互补条件。若 $\alpha_i^\ast>0$，则 $c_i(x^\ast)=0$。

回到原问题

注意：支持向量在不同的地方存在着不同的定义。第一种定义：在间隔超平面上，即使等式 $y(x_1w_1+x_2w_2+b)=1$ 成立的数据点为支持向量；第二种定义，当 $\alpha_i>0$ 时，样本 $x_i$ 为支持向量，第二种的支持向量与第一种的差别关键在于 $\alpha_i(y_i(w^Tx_i+b)-1)=0$ 中，$\alpha_i$ 和 $y_i(w^Tx_i+b)-1$ 有可能同时为0。

$$\left\{\begin{aligned}&\min_{w,b}\frac{1}{2}||w||^2\\&\text{s.t. }y_i\left({w^Tx_i+b}\right)-1\ge 0 \ ,i=1,\dots,N \\\end{aligned}\right.\Leftrightarrow\left\{\begin{aligned}&\min_{w,b}\frac{1}{2}||w||^2\\&\text{s.t. }1-y_i\left({w^Tx_i+b}\right)\le 0,\ i=1,\dots,N \\\end{aligned}\right.$$

拉格朗日函数为：

$$\begin{align*} L(w,b,\alpha)&=\frac{1}{2}||w||^2+\sum_{i=1}^N\alpha_i-\sum_{i=1}^N\alpha_iy_iw^Tx_i-\sum_{i=1}^N\alpha_iy_ib\\ &=\frac{1}{2}w^Tw+\sum_{i=1}^N\alpha_i-w^T(\sum_{i=1}^N\alpha_iy_ix_i)-b(\sum_{i=1}^N\alpha_iy_i) \end{align*}$$

求偏导取极值，

$$\left\{\begin{aligned}&\frac{\partial L(w,b,\alpha)}{\partial w}=w-\sum_{i=1}^N\alpha_iy_ix_i=0,\ w=\sum_{i=1}^N\alpha_iy_ix_i\\&\frac{\partial L(w,b,\alpha)}{\partial b}=-\sum_{i=1}^N\alpha_iy_i=0\end{aligned}\right.$$

将上式结果代入，

$$\begin{align*} L(w,b,\alpha)&=\frac{1}{2}w^Tw+\sum_{i=1}^N\alpha_i-w^T(\sum_{i=1}^N\alpha_iy_ix_i)-b(\sum_{i=1}^N\alpha_iy_i)\\ &=\frac{1}{2}w^Tw+\sum_{i=1}^N\alpha_i-w^Tw\\ &=\sum_{i=1}^N\alpha_i-\frac{1}{2}w^Tw\\ &=-\frac{1}{2}\langle\sum_{i=1}^N\alpha_iy_ix_i,\sum_{i=1}^N\alpha_iy_ix_i\rangle+\sum_{i=1}^N\alpha_i \end{align*}$$

即问题变成了，

$$\min_{w,b}L(w,b,\alpha)=-\frac{1}{2}\langle\sum_{i=1}^N\alpha_iy_ix_i,\sum_{i=1}^N\alpha_iy_ix_i\rangle+\sum_{i=1}^N\alpha_i$$

它的对偶问题为：

$$\left\{\begin{aligned}&\max_{\alpha_i:\alpha_i\ge0}-\frac{1}{2}\langle\sum_{i=1}^N\alpha_iy_ix_i,\sum_{i=1}^N\alpha_iy_ix_i\rangle+\sum_{i=1}^N\alpha_i\\&\text{s.t. }\sum_{i=1}^N\alpha_iy_i=0\end{aligned}\right.$$

符号取反，

$$\left\{\begin{aligned}&\min_{\alpha_i:\alpha_i\ge0}\frac{1}{2}\langle\sum_{i=1}^N\alpha_iy_ix_i,\sum_{i=1}^N\alpha_iy_ix_i\rangle-\sum_{i=1}^N\alpha_i\\&\text{s.t. }\sum_{i=1}^N\alpha_iy_i=0\end{aligned}\right.$$

把目标函数展开，

$$\left\{\begin{aligned}&\min_{\alpha_i:\alpha_i\ge0}\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_j\langle x_i,x_j\rangle-\sum_{i=1}^N\alpha_i\\&\text{s.t. }\sum_{i=1}^N\alpha_iy_i=0\end{aligned}\right.$$

设 $\alpha^\ast=(\alpha_1^\ast,\alpha_2^\ast,\dots,\alpha_N^\ast)$ 是上述问题的最优解，设原始问题最优解 $w^\ast$ 和 $b^\ast$。

KKT条件为：

$$\left\{\begin{aligned} &\nabla_wL(w^\ast,b^\ast,\alpha^\ast)=w^\ast-\sum_{i=1}^N\alpha_i^\ast y_ix_i=0,w^\ast=\sum_{i=1}^N\alpha_i^\ast y_ix_i\\ &\nabla_bL(w^\ast,b^\ast,\alpha^\ast)=-\sum_{i=1}^N\alpha_i^\ast y_i=0\\ &\alpha_i^\ast(y_i(w^\ast x_i+b^\ast)-1)=0,i=1,\dots,N\\ &y_i(w^\ast x_i+b^\ast)-1\ge0,i=1,\dots,N\\ &\alpha_i\ge 0,i=1,\dots,N \end{aligned}\right.$$

由上式可得：$\exists \alpha_j^\ast\neq 0\Leftrightarrow y_j(w^{\ast T}x_j+b^\ast)=1$，因为 $y_j=\pm 1$，所以 $w^{\ast T} x_j+b^\ast=\pm 1$，得：

$$y_j(w^{\ast T} x_j+b^\ast)=1=y_j^2\ \Rightarrow\ b^\ast=y_j-w^{\ast T} x_j=y_j-\sum_{i=1}^N\alpha_i^\ast y_i x_i ·x_j$$

这表示只有部分点能起到作用，即支持向量。

最终，分类决策函数为：

$$f(x)=\text{sign}(\sum_{i=1}^N\alpha_i^\ast y_i(x_i·x)+b^\ast)$$

其中，$b^\ast=y_j-w^{\ast T} x_j$。

注意：有两个地方需要求偏导：1. KKT；2. 取极值。KKT条件中的偏导为0，是最优解 $w^\ast, b^\ast$；取极值，是在原问题得出后，偏导为0，代入，得对偶问题。

软间隔最大化

引入松弛变量，$\xi_i$ 表示样本逾越的程度。

• $\xi_i<0,y_i(w^Tx_i+b)=1-\xi_i>1$，未逾越且很远；
• $\xi_i=0,y_i(w^Tx_i+b)=1-\xi_i=1$，在最大间隔边界上；
• $0<\xi_i<1$，跨越了最大间隔边界，依然线性可分不过间隔变窄；
• $\xi_i=1,y_i(w^Tx_i+b)=1-\xi_i=0$，到达“中线”，间隔称为“线”；
• $\xi_i>1,y_i(w^Tx_i+b)=1-\xi_i<0$，跨越了“中线”，线性不可分。

引入松弛变量，约束条件变为：

$$y_i(w^Tx_i+b)\ge 1-\xi_i$$

即允许部分犯错。

问题变成：

$$\left\{\begin{aligned}&\min_{w,b,\xi}\frac{1}{2}||w||^2+C\sum_{i=1}^N\xi_i\\&\text{s.t. }y_i\left({w^Tx_i+b}\right)\ge1-\xi_i \ ,i=1,\dots,N \\&\xi_i\ge0,i=1,\dots,N\end{aligned}\right.$$

拉格朗日函数：

$$L(w,b,\xi,\alpha,\mu)=\frac{1}{2}||w||^2+C\sum_{i=1}^N\xi_i-\sum_{i=1}^N\alpha_i(y_i(w^Tx_i+b)-1+\xi_i)-\sum_{i=1}^N\mu_i\xi_i$$

求偏导为：

$$\left\{\begin{aligned}&\frac{\partial L(w,b,\xi,\alpha,\mu)}{\partial w}=w-\sum_{i=1}^N\alpha_iy_ix_i=0, w=\sum_{i=1}^N\alpha_iy_ix_i\\&\frac{\partial L(w,b,\xi,\alpha,\mu)}{\partial b}=-\sum_{i=1}^N\alpha_iy_i=0\\&\frac{\partial L(w,b,\xi,\alpha,\mu)}{\partial \xi_i}=C-\alpha_i-\mu_i=0\end{aligned}\right.$$

代入拉格朗日函数，

$$\begin{align*} \min_{w,b,\xi}L(w,b,\xi,\alpha,\mu)&=-\frac{1}{2}w^Tw+C\sum_{i=1}^N\xi_i+\sum_{i=1}^N\alpha_i-\sum_{i=1}^N\alpha_i\xi_i-\sum_{i=1}^N\mu_i\xi_i\\ &=-\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_j\langle x_i,x_j\rangle+\sum_{i=1}^N\alpha_i \end{align*}$$

对偶问题为：

$$\left\{\begin{aligned}&\max_\alpha&&-\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^Ny_iy_j\langle x_i,x_j\rangle+\sum_{i=1}^N\alpha_i\\&\text{s.t. }&&\sum_{i=1}^N\alpha_iy_i=0\\&&&c-\alpha_i-y_i=0\\ &&&\alpha_i\ge 0\\ &&&\mu_i\ge 0,i=1,2,\dots,N\end{aligned}\right.$$

符号取反为：

$$\left\{\begin{aligned}&\min_\alpha&&\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^Ny_iy_j\langle x_i,x_j\rangle-\sum_{i=1}^N\alpha_i\\&\text{s.t. }&&\sum_{i=1}^N\alpha_iy_i=0\\&&&0\le \alpha_i\le C,i=1,\dots,N\end{aligned}\right.$$

同样地，设 $\alpha^\ast$ 为上述问题最优解，原始最优解为 $w^\ast$ 和 $b^\ast$。

KKT条件为：

$$\left\{\begin{aligned} &\nabla_wL(w^\ast,b^\ast,\xi^\ast,\alpha^\ast,\mu^\ast)=w^\ast-\sum_{i=1}^N\alpha_i^\ast y_ix_i=0, w^\ast=\sum_{i=1}^N\alpha_i^\ast y_ix_i\\ &\nabla_b L(w^\ast,b^\ast,\xi^\ast,\alpha^\ast,\mu^\ast)=-\sum_{i=1}^N\alpha_i^\ast y_i=0\\ &\nabla_{\xi_i} L(w^\ast,b^\ast,\xi^\ast,\alpha^\ast,\mu^\ast)=C-\alpha_i^\ast-\mu_i^\ast=0\\ &\alpha_i^\ast(y_i(w^{\ast T} x_i+b^\ast)-1+\xi_i^\ast)=0,i=1,\dots,N\\ &\mu_i^\ast \xi_i^\ast =0,i=1,\dots,N\\ &y_i(w^{\ast T} x_i+b^\ast)-1+\xi_i^\ast\ge0,i=1,\dots,N\\ &\xi_i\ge 0,\alpha_i\ge 0,\mu_i\ge 0,i=1,\dots,N \end{aligned}\right.$$

若存在 $\alpha^\ast$ 的一个分量 $\alpha_j^\ast$，

1. $0<\alpha_j^\ast<C$，则： $$y_i(w^{\ast T} x_i+b^\ast)-1=0$$ 则： $$b^\ast=y_j-\sum_{i=1}^Ny_i\alpha_i^\ast\langle x_i,x_j\rangle$$ 决策函数为： $$f(x)=\text{sign}(\sum_{i=1}^N\alpha_i^\ast y_i\langle x,x_i\rangle+b^\ast)$$ 满足该条件的点位于间隔平面上，即都满足于 $y_i(w^{\ast T} x_i+b^\ast)-1=0$ ，且 $w^\ast, b^\ast$ 都唯一。

1. $a_j^\ast=0$，全部为零，不可能，不考虑。

1. $\alpha_j^\ast=C$，即去除$\alpha_j^\ast=0$ 的数据点后，没有 $0<\alpha_j^\ast<C$ 的数据点，只有 $\alpha_j^\ast=C$ 中，正负样本的数量相等。这种情况下，根据 $$\alpha_i^\ast(y_i(w^{\ast T} x_i+b^\ast)-1+\xi_i^\ast)=0$$ 则： $$y_i(w^{\ast T} x_i+b^\ast)-1+\xi_i^\ast=0$$ 之后，取 $\xi_i^\ast$ 的值，平均多个点求出 $\bar b^\ast$。

序列最小最优化算法

当训练样本容量很大时，我们需要一些快速实现算法来提高支持向量机的效率。

SMO算法是一种启发式算法，其基本思路是：如果所有变量的解都满足此最优化问题的KKT条件，那么这个最优化问题的解就得到了。整个SMO算法包括两个部分：求解两个变量二次规划的解析方法和选择变量的启发式方法。下面给出SMO算法的具体步骤：

1. 取初值 $\alpha^{(0)}=0$，令 $k=0$；
2. 选取优化变量 $\alpha_1^{(k)},\alpha_2^{(k)}$，解析求解两个变量的最优化问题，求得最优解 $\alpha_1^{(k+1)},\alpha_2^{(k+1)}$，更新为 $\alpha$ 为 $\alpha^{(k+1)}$；
3. 若在精度 $\varepsilon$ 范围内满足停机条件，则转(4)；否则令 $k=k+1$，转(2)；
4. 取 $\hat \alpha = \alpha^{(k+1)}$。

SVM的损失函数

SVM的损失函数为：

$$\min_{w,b}\sum_{i=1}^N\max(0,1-y_i(w^Tx_i+b))+\lambda||w||^2$$

原始最优化问题是：

$$\left\{\begin{aligned}&\min_{w,b,\xi}\frac{1}{2}||w||^2+C\sum_{i=1}^N\xi_i\\&\text{s.t. }y_i\left({w^Tx_i+b}\right)\ge1-\xi_i \ ,i=1,\dots,N \\&\xi_i\ge0,i=1,\dots,N\end{aligned}\right.$$

要让两者等价，首先有：

$$\xi_i\ge1-y_i(w^Tx_i+b)$$

也就是说，当样本点被正确分类且函数间隔 $y_i(w^Tx_i+b)$ 大于1时，损失为0，否则损失是 $1-y_i(w^Tx_i+b)$ ，图表示为：

称为合页损失（或铰链损失，hinge loss）。

对比一下：

• 感知机：$L(w,b)=-\sum_{x_i\in M}y_i(w^Tx_i+b)=\sum_{i=1}^N\max(0,-y_i(w^Tx_i+b))$
• SVM：$L(w,b)=-\sum_{x_i\in M}y_i(w^Tx_i+b)=\sum_{i=1}^N\max(0,1-y_i(w^Tx_i+b))$

其中绿线表示SVM的损失函数，紫线表示感知机的损失函数，当 $f(x)$ 大于0的时候，感知机的loss直接就为0了，也就是说，它只需要分对就行了，“做事”比较粗糙，而对于SVM来说，$f(x)$ 大于0，损失依然存在，直到 $f(x)$ 大于等于1，它才对分类结果比较满意，loss才降为0，换句话来说，SVM“做事”比较严谨。

而当 $\lambda=\frac{1}{2C}$ 时，可以看出SVM是建立在结构风险最小化的基础上的。

核函数

其中 $\langle x_i,x_j\rangle$ 计算的是数据矩阵的第 $i$ 维与第 $j$ 维的内积，这其实是比较低维的输入空间，如果将它通过核技巧，将它映射到高维，变得线性可分，再依靠SVM在线性可分问题上强大的分类能力，问题似乎变得简单起来了。

常见的核函数有：

• Linear：$k(x_i,x_j)=x_i^Tx_j$
• Polynomial：$k(x_i,x_j)=(ax_i^Tx_j+c)^d$
• Gaussian：$k(x_i,x_j)=\exp(-\frac{||x_i-x_j||^2}{2\sigma^2})$
• Sigmoid：$k(x_i,x_j)=\tanh (ax_i^Tx_j+c)$

Logistic回归分类模型

（注：因为上课所讲的和PPT里所用的符号有冲突，这里个人统一了一下，可能会有遗漏。）

接上节课所说的，Logistic回归分类模型的对数损失函数为：

$$L(w)=\sum_{i=1}^n[-(y_i\text{ln}\sigma_i+(1-y_i)\text{ln}(1-\sigma_i))],\tag 1$$

其中，$\sigma_i=\frac{1}{1+\text{exp}(-w^Tx_i)}$。

梯度下降法

求解损失函数的梯度，即对 $w$ 的偏导为：

$$\frac{\partial L}{\partial w}=\sum_{i=1}^n(\sigma_i-y_i)x_i,$$

其矩阵形式为：

$$\frac{\partial L}{\partial w}=X^T(\vec \sigma-\vec y).\tag 2$$

利用梯度下降法的算法流程为：

牛顿法

简单地说，牛顿法是二阶逼近，梯度下降法是一阶逼近。

牛顿法是利用泰勒级数来进行逼近。首先，对 $f(x)$ 的级数展开至两次，

$$f(x_0+\Delta x)=f(x_0)+f'(x_0)\Delta x+\frac{f''(x_0)}{2}\Delta x^2+R(\Delta x)$$

现假设：

$$g(\Delta x)=f(x_0)+f'(x_0)\Delta x+\frac{f''(x_0)}{2}\Delta x^2$$

$g(\Delta x)$ 是一个关于 $\Delta x$ 的一元二次方程，一元二次方程的最小值求法：

$$f(x)=ax^2+bx+c,$$

$f(x)$ 的极值点为 $-\frac{b}{2a}$。所以，$g(\Delta x)$ 的极值为：

$$-\frac{f'(x_0)}{f''(x_0)},$$

此时 $f(x_0+\Delta x)$ 的极值为：

$$x_1 = x_0-\frac{f'(x_0)}{f''(x_0)},$$

以此类推，

$$x_n=x_{n-1}-\frac{f'(x_{n-1})}{f''(x_{n-1})}$$

当 $f(x)$ 是一个多元函数时，则为：

$$x_n=x_{n-1}-\frac{\nabla f(x_{n-1})}{\mathbb{H}(f(x_{n-1}))}$$

---

将牛顿法应用到Logistic回归分类模型，则总的算法流程为：

其中，$g$ 为一阶导 $\frac{\partial L}{\partial w}$，所以我们需要求出 $L(w)$ 的二阶偏导，

$$\frac{\partial^2L(w)}{\partial w_j \partial w_k}=\sum_{i=1}^nx_i^{(j)}\sigma_i(1-\sigma_i)x_i^{(k)},$$

其中，$j$ 和 $k$ 表示的是第几维，$i$ 表示的是第几个样本。写成矩阵形式为：

$$\mathbb{H}=X^TBX,\tag 3$$

其中，

$$B=\left(\begin{align}\sigma_1(1-\sigma_1)\\&\ddots\\&& \sigma_n(1-\sigma_n) \end{align}\right).$$

$w^{t+1}$ 的更新方式如上图所示，

$$w^{t+1}\leftarrow w^t-\mathbb{H}^{-1}g,$$

将 $L(w)$ 的一阶偏导和二阶偏导，即(2)式和(3)式代入得，

$$w^{t+1}\leftarrow w^t-(X^TBX)^{-1}X^T(\vec \sigma-\vec y).$$

对上式进行变换，

$$\begin{align*} w^{t+1}&=w^t-(X^TBX)^{-1}X^T(\vec \sigma - \vec y)\\ &=(X^TBX)^{-1}X^TBXw^t-(X^TBX)^{-1}X^TBB^{-1}(\vec \sigma - \vec y)\\ &=(X^TBX)^{-1}X^TB(Xw^t-B^{-1}(\vec \sigma-\vec y))\\ &=(X^TBX)^{-1}X^TBS^t \end{align*}\tag 4$$

其中，$S^t=Xw^t-B^{-1}(\vec \sigma-\vec y)$。

---

这里先不管上面的，介绍一种算法：迭代重加权最小二乘法（Iterative Reweighted Least Squares，IRLS），即是说，每次迭代更新参数时，它会有一个随着迭代次数 $t$ 变化的权重 $\widetilde W^{t}$ 在影响。

它被应用于p范数损失函数的优化，本质上是用2范数来近似替代p范数。

算法表示为：

$$\beta^{(t+1)}=\text{arg min}_\beta \sum_{i=1}^n\tilde w_i^t|\tilde y_i-\widetilde X_i\beta|^2=(\widetilde X^T\widetilde W^t\widetilde X)^{-1}\widetilde X^T\widetilde W^t\widetilde y.$$

与上面的(4)式比对，

$$\begin{align*} \widetilde X&=X\\ \widetilde W^t&=B\\ \tilde y&=S^t \end{align*}$$

---

所以，(4)式在某种程度上，可以看作是以下问题的梯度下降法，

$$\text{arg min}_w \sum_{i=1}^n\sigma_i(1-\sigma_i)(w^Tx_i-S_i^t).$$

引用

1. 牛顿法和梯度下降法的学习
2. Iteratively reweighted least squares

产生式模型

之前我们学习的几种算法，包括Logistic回归、SVM和感知机，都是属于判别式模型。它们是直接对 $p(y|x)$ 进行建模或者直接学习输入空间到输出空间的映射关系，其中，$x$ 是某类样例的特征，$y$ 是某类样例的分类标记。

今天要讲的产生式模型，又称“生成式模型”。其中经典算法有：朴素贝叶斯、隐式马尔可夫模型和高斯判别模型。产生式模型实对 $p(x|y)$ （条件概率）和 $p(y)$ （先验概率）进行建模，然后按贝叶斯法则求出后验概率 $p(y|x)$：

$$p(y|x)=\frac{p(x|y)p(y)}{p(x)},$$

使得后验概率最大的类别 $y$ 即是新样例的预测值：

$$\text{arg max}_y p(y|x)=\text{arg max}_y \frac{p(x|y)p(y)}{p(x)} = \text{arg max}_y p(x|y)p(y).$$

朴素贝叶斯

朴素贝叶斯主要解决的是特征为离散值的问题。

它主要基于以下两个假设：

1. 各个维度独立，即 $P(X_1,\cdots,X_d|Y)=\prod_{j=1}^dP(X_j|Y)$，这也是它被称为”朴素“的原因；
2. 各个维度的取值为0或1，即 $X_i\in\{0, 1\}$。

而与朴素贝叶斯的假设相反的，称为贝叶斯网络，它假设各个维度都不独立。贝叶斯网络还有很多名字：有向无环概率图、信念网络、置信网络等，有些听起来感觉很叼，但它们都指的是贝叶斯网络。

---

假设现在有一些数据 $D$，标签为0或1，

$$D=\{\langle x_i,y_i\rangle\}_{i=1}^n,x_i\in\mathbb{R}^d,y_i\in\{0,1\}$$

二分类问题，可用伯努利分布来刻画：

1. 标签的分布为：

	$y_i=1$	$y_i=0$
$p(y_i)$	$\varphi_1$	$1-\varphi_1$

$$p(y_i)=(\varphi_1)^{y_i}(1-\varphi_1)^{1-y_i}.\tag 5$$

1. 当样本的标签为0时，其第 $j$ 维特征的分布为：

	$x_{ij}=1$	$x_{ij}=0$
$p(x_{ij}	y_i=0)$	$\varphi_{j0}$	$1-\varphi_{j0}$

$$p(x_{ij}|y_i=0)=(\varphi_{j0})^{x_{ij}}(1-\varphi_{j0})^{1-x_{ij}},$$

其中，$x_{ij}$ 表示的是第 $i$ 样本的第 $j$ 维特征。

1. 当样本的标签为1时，其第 $j$ 维特征的分布为：

	$x_{ij}=1$	$x_{ij}=0$
$p(x_{ij}	y_i=1)$	$\varphi_{j1}$	$1-\varphi_{j1}$

$$p(x_{ij}|y_i=1)=(\varphi_{j1})^{x_{ij}}(1-\varphi_{j1})^{1-x_{ij}},$$

似然函数为：

$$\begin{align*} L(\phi)&=\prod_{i=1}^n p(\langle x_i, y_i\rangle)\\ &=\prod_{i=1}^n \left[p(y_i)p(x_i|y_i)\right]\\ &=\prod_{i=1}^n \left[p(y_i)\prod_{j=1}^d p(x_{ij}|y_i)\right]\\ &=\prod_{i=1}^n \left[p(y_i)\times\left(\prod_{j=1}^dp(x_{ij}|y_i=1)\right)^{1(y_i=1)}\times \left(\prod_{j=1}^dp(x_{ij}|y_i=0)\right)^{1(y_i=0)}\right] \end{align*}$$

对数似然函数为：

$$\begin{align*} \text{ln}L(\phi)&=\text{ln}\prod_{i=1}^n \left[p(y_i)\times\left(\prod_{j=1}^dp(x_{ij}|y_i=1)\right)^{1(y_i=1)}\times \left(\prod_{j=1}^dp(x_{ij}|y_i=0)\right)^{1(y_i=0)}\right]\\ &=\sum_{i=1}^n\left[\text{ln}p(y_i)+1(y_i=1)\text{ln}\left(\prod_{j=1}^dP(x_{ij}|y_i=1)\right)+1(y_i=0)\text{ln}\left(\prod_{j=1}^dP(x_{ij}|y_i=0)\right)\right]\\ &=\sum_{i=1}^n\left[\text{ln}\left((\varphi_1)^{y_i}(1-\varphi_1)^{1-y_i}\right)+1(y_i=0)\text{ln}\left(\prod_{j=1}^d(\varphi_{j0})^{x_{ij}}(1-\varphi_{j0})^{1-x_{ij}}\right)\\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ +1(y_i=1)\text{ln}\left(\prod_{j=1}^d(\varphi_{j1})^{x_{ij}}(1-\varphi_{j1})^{1-x_{ij}}\right)\right]\\ &=\sum_{i=1}^n\left(y_i\text{ln}\varphi_1+(1-y_i)\text{ln}(1-\varphi_1)\right)+\sum_{i=1}^n\left[1(y_i=0)\sum_{j=1}^d\left(x_{ij}\text{ln}\varphi_{j0}+(1-x_{ij})\text{ln}(1-\varphi_{j0})\right)\right]\\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ +\sum_{i=1}^n\left[1(y_i=1)\sum_{j=1}^d\left(x_{ij}\text{ln}\varphi_{j1}+(1-x_{ij})\text{ln}(1-\varphi_{j1})\right)\right] \end{align*}$$

为了最大化对数似然函数，

$$\varphi_{MLE}=\text{arg max}_\varphi L(\phi)=\text{arg max}_\varphi \text{ln}L(\phi),$$

对其各个参数求偏导，并令偏导数为零，即：

$$\begin{align*} &\frac{\partial \text{ln}L}{\partial \varphi_1}=\sum_{i=1}^n\left(y_i·\frac{1}{\varphi_1}-(1-y_i)·\frac{1}{(1-\varphi_1)}\right)=0\tag 6\\ &\frac{\partial \text{ln}L}{\partial \varphi_{j0}}=\sum_{i=1}^n\left(1(y_i=0)\left(x_{ij}·\frac{1}{\varphi_{j0}}-(1-x_{ij})·\frac{1}{(1-\varphi_{j0})}\right)\right)=0 \tag 7\\ &\frac{\partial \text{ln}L}{\partial \varphi_{j1}}=\sum_{i=1}^n\left(1(y_i=1)\left(x_{ij}·\frac{1}{\varphi_{j1}}-(1-x_{ij})·\frac{1}{(1-\varphi_{j1})}\right)\right)=0 \tag 8\\ \end{align*}$$

接下来，我们从(6)式开始推，

$$\begin{align*} \sum_{i=1}^n\left(y_i·\frac{1}{\varphi_1}-(1-y_i)·\frac{1}{(1-\varphi_1)}\right)&=0\\ \sum_{i=1}^n\left(\frac{y_i(1-\varphi_1)-(1-y_i)\varphi_1}{\varphi_1(1-\varphi_1)}\right)&=0\\ \sum_{i=1}^n\left(\frac{y_i-\varphi_1}{\varphi_1(1-\varphi_1)}\right)&=0\\ \frac{\sum_{i=1}^n(y_i-\varphi_1)}{\varphi_1(1-\varphi_1)}&=0\\ \sum_{i=1}^n(y_i-\varphi_1)&=0\\ \sum_{i=1}^ny_i-n\varphi_1&=0\\ \varphi_1 &=\frac{\sum_{i=1}^n1\{y_i=1\}}{n}\tag 9 \end{align*}$$

然后是(7)式，类似上式，步骤稍微跳跃一下

$$\begin{align*} \sum_{i=1}^n\left(1(y_i=0)\left(x_{ij}·\frac{1}{\varphi_{j0}}-(1-x_{ij})·\frac{1}{(1-\varphi_{j0})}\right)\right)&=0\\ \sum_{i=1}^n\left(1(y_i=0)\left(\frac{x_{ij}(1-\varphi_{j0})-(1-x_{ij})\varphi_{j0}}{\varphi_{j0}(1-\varphi_{j0})}\right)\right)&=0\\ \sum_{i=1}^n\left(1(y_i=0)\left(x_{ij}-\varphi_{j0}\right)\right)&=0\\ \sum_{i=1}^n1\{y_i=1\wedge x_{ij}=1\}-\sum_{i=1}^n1\{y_i=0\}\varphi_{j0}&=0\\ \varphi_{j0}&=\frac{\sum_{i=1}^n1\{y_i=0\wedge x_{ij}=1\}}{\sum_{i=1}^n1\{y_i=0\}}\tag{10} \end{align*}$$

同理可得，

$$\varphi_{j1}=\frac{\sum_{i=1}^n1\{y_i=1\wedge x_{ij}=1\}}{\sum_{i=1}^n1\{y_i=1\}}.\tag{11}$$

推导完成。

因为朴素贝叶斯模型简单，所以只需要少量样本就足够训练。

零概率问题：零概率问题就是在计算新实例的概率时，如果某个分量在训练集中从没有出现过，会导致整个实例的概率计算结果为0，这显然是不合理的。为了解决这个问题，引出拉普拉斯平滑：

$$P(w|c)=\frac{\text{num}(w,c)+\varepsilon}{\text{num}(c)+2\varepsilon}$$

属性缺失问题：如果数据中某个可能维度的值缺失，可以不用管，不会受到影响。原因是不用管意味着计算概率时忽略该维度的影响，即没有对该维度加入先验。

而在最开始提到的贝叶斯网络，即概率图，该模型中每个结点之间都是相关的，所以复杂度最大能达到 $C_d^2$ ，其中 $d$ 为样本 $x_i$ 的特征维数。

策略

朴素贝叶斯法将实例分到后验概率最大的类中，这等价于期望风险最小化。假设选择0-1损失函数：

$$L(Y,f(X))=\left\{\begin{align}&1,&&Y\neq f(X)\\&0,&&Y=f(X)\end{align}\right.$$

式中 $f(X)$ 是分类决策函数。这时，期望风险函数为，

$$\begin{align*} R_{\text{exp}}(f)&=E[L(Y,f(X))]\\ &=E_X \left[\sum_{k=1}^KL(c_k,f(X))P(c_k|X)\right] \end{align*}$$

为了使期望风险最小化，只需对 $X=x$ 逐个极小化，

$$Y\in\{c_1,c_2,\dots,c_K\}\\ \begin{align*} R_{\text{exp}}(f)&=E[L(Y,f(x))]\\ &=\iint L(Y=y,f(X=x))P(X=x,Y=y)dxdy\\ &=\iint L(Y=y,f(X=x))P(Y=y|X=x)p(X=x)dxdy\\ &=\int \left\{\int[L(Y=y,f(X=x))P(Y=y|X=x)]dy\right\}p(X=x)dx\\ &=\int \left\{\sum_{k=1}^KL(c_k,f(X=x))P(c_k|X=x)\right\}p(X=x)dx\\ &=E_X \left[\sum_{k=1}^KL(c_k,f(X))P(c_k|X)\right] \end{align*}$$

由此得到：

$$\begin{align*} f(x)&=\text{arg min}_{y\in Y}\sum_{k=1}^KL(c_k,y)P(c_k|X=x)\\ &=\text{arg min}_{y\in Y}\sum_{k=1}^KP(y\neq c_k|X=x)\\ &=\text{arg min}_{y\in Y}(1-P(y=c_k|X=x))\\ &=\text{arg max}_{y\in Y}P(y=c_k|X=x) \end{align*}$$

这样一来，根据期望风险最小化准则就得到了后验概率最大化准则：

$$f(x)=\text{arg max}_{c_k}P(c_k|X=x)$$

高斯判别分析

在朴素贝叶斯中，我们假设 $x$ 是离散的，服从多项分布（包括伯努利分布）。在高斯判别分析（Gaussian discriminant analysis，GDA）中，我们假设 $x$ 是连续的，服从高斯分布。根据以下公式，判别样本 $x$ 的标签。

$$\begin{align*} y^\ast&=\text{arg max}_{y\in\{0,1\}}p(y|x)\\ &=\text{arg max}_{y\in\{0,1\}}p(y)p(x|y)\\ &\Leftrightarrow\frac{p(y=1)p(x|y=1)}{p(y=0)p(x|y=0)}?>1\\ &\Leftrightarrow\text{ln}\left[\frac{p(y=1)p(x|y=1)}{p(y=0)p(x|y=0)}\right]?>0 \end{align*}$$

该模型基于两点假设，$y$ 服从伯努利分布，特征 $x$ 是多元高斯分布。写出分布概率密度为：

$$\begin{align*} &p(y=1)=\phi,&&p(y=0)=1-\phi\\ &p(x|y=1)\sim N(\mu_1,\Sigma),&&p(x|y=1)=(2\pi)^{-\frac{d}{2}}|\Sigma|^{-\frac{1}{2}}\text{exp}\left\{-\frac{1}{2}(x-\mu_1)^T\Sigma(x-\mu_1)\right\}\\ &p(x|y=0)\sim N(\mu_0,\Sigma),&&p(x|y=0)=(2\pi)^{-\frac{d}{2}}|\Sigma|^{-\frac{1}{2}}\text{exp}\left\{-\frac{1}{2}(x-\mu_0)^T\Sigma(x-\mu_0)\right\} \end{align*}$$

所以，比较 $p(y=1)p(x|y=1)$ 和 $p(y=0)p(x|y=0)$ 的大小：

$$\begin{align*} &\text{ln}\left[\frac{p(y=1)p(x|y=1)}{p(y=0)p(x|y=0)}\right]\\ =&\text{ln}\frac{\phi \times(2\pi)^{-\frac{d}{2}}|\Sigma|^{-\frac{1}{2}}\text{exp}\left\{-\frac{1}{2}(x-\mu_1)^T\Sigma(x-\mu_1)\right\}}{(1-\phi) \times(2\pi)^{-\frac{d}{2}}|\Sigma|^{-\frac{1}{2}}\text{exp}\left\{-\frac{1}{2}(x-\mu_0)^T\Sigma(x-\mu_0)\right\}}\\ =&(\mu_1^T\Sigma-\mu_0^T\Sigma)x+\frac{1}{2}\mu_0^T\Sigma\mu_0-\frac{1}{2}\mu_1\Sigma\mu_1+\text{ln}\frac{\phi}{1-\phi}\\ =&w^Tx+b \end{align*}$$

可以看到GDA实际上是一个线性模型，这是一个判别式模型。产生式模型有其对应的判别式模型，这样的组合叫做“产生式-判别式”对；判别式模型不一定有其对应的产生式模型。

• 训练样本多的情况下，采用判别式；
• 训练样本少的情况下，采用产生式。

隐马尔可夫模型

隐马尔可夫模型是关于时序的概率模型，描述由一个隐藏的马尔可夫链随机生成不可观测的状态随机序列，再由各个状态生成一个观测而产生观测随机序列的过程。

隐马尔可夫链随机生成的状态的序列，成为状态序列；每个状态生成一个观测，而由此产生的观测的随机序列，称为观测序列。序列的每一个位置又可以看作是一个时刻。

隐马尔可夫模型由初始概率分布、状态转移概率和观测概率（又称“发射概率”）确定。隐马尔可夫模型的形式定义如下，设 $Q$ 是所有可能的状态的集合，$V$ 是所有可能的观测的集合，

• 状态序列：$O_t\in Q=\{q_1,q_2,\dots,q_N\}$
• 观测序列：$I_t\in V=\{v_1,v_2,\dots,v_M\}$

其中，$N$ 是可能的状态数，$M$ 是可能的观测数。

隐马尔可夫模型可用三元符号表示，即 $\lambda=(A, B, \Pi)$，

• $A$ —状态转移概率矩阵：$A=[a_{ij}]_{N\times N}, a_{ij}=P(I_{t+1}=q_j|I_t=q_i)$
• $B$—观测概率矩阵：$B=[b_j(k)]_{N\times M},b_j(k)=P(O_t=v_k|I_t=q_j)$
• $\Pi$—初始状态概率矩阵：$\Pi=(\pi_i)_{1\times N},\pi_i=P(I_1=q_i)$

隐马尔可夫模型有三个基本问题：

1. 概率计算问题：给定模型 $\lambda=(A,B,\Pi)$ 和观测序列 $O=(o_1,o_2,\dots,o_T)$，计算在模型 $\lambda$ 下观测序列 $O$ 出现的概率 $P(O|\lambda)$。
2. 学习问题：已知观测序列 $O=(o_1,o_2,\dots,o_T)$，估计模型 $\lambda=(A,B,\Pi)$ 参数，使得在该模型下观测序列概率 $P(O|\lambda)$ 最大。即用极大似然估计的方法估计参数。
3. 预测问题：也成为解码问题。已知模型 $\lambda=(A,B,\Pi)$ 和观测序列 $O=(o_1,o_2,\dots,o_T)$，求对给定观测序列条件概率 $P(I|O)$ 最大的状态序列 $I=(i_1,i_2,\dots,i_T)$。

第一个问题对应的算法是前向—后向算法（Forward-Backward algorithm），第二个则是模型训练问题，第三个问题可用维特比算法解决。

推理（概率计算）问题

已知：模型 $\lambda=(A,B,\Pi)$ 和观测序列 $O=(o_1,o_2,\dots,o_T)$，求：$P(O|\lambda)$。

1.直接计算法

通过列举所有可能的长度为 $T$ 的状态序列 $I=(i_1,i_2,\dots,i_T)$，求各个状态序列 $I$ 与观测序列 $O=(o_1,o_2,\dots,o_T)$ 的联合概率 $P(O,I|\lambda)$，然后对所有可能的状态序列求和，得到 $P(O|\lambda)$。

$$\begin{align*} P(O)&=P(O_1,O_2,\dots,O_T) \\ &=\sum_I P(O,I)=\sum_I \left[P(O|I)P(I)\right]\\ &=\sum_I \left[P(O_1,O_2,\dots,O_T|I_1,I_2,\dots,I_T)P(I_1,I_2,\dots,I_T)\right]\\ &=\sum_I \left[P(O_1|I_1)\times P(O_2|I_2)\times \dots \times P(O_T|I_T)\times P(I_1)\times P(I_2|I_1)\times \dots\times P(I_T|I_{T-1})\right]\\ &= \sum_{I_1,\dots,I_T}\left[\left(\prod_{t=1}^TP(O_t|I_t)\right)\times\left(\prod_{t=1}^{T-1}P(I_{t+1}|I_t)\right)\times P(I_1)\right]\\ &=\sum_{I_1}\sum_{I_2}\sum_{I_3}\sum_{I_4}\dots\sum_{I_T}\left[\left(\prod_{t=1}^TP(O_t|I_t)\right)\times\left(\prod_{t=1}^{T-1}P(I_{t+1}|I_t)\right)\times P(I_1)\right] \end{align*}$$

稍微分析一下就可以发现，对于某个长度为 $T$ 的状态序列，求和需要花费 $O(T)$ 时间，由于状态序列共有 $N^T$ 种，故直接计算的算法时间复杂度为 $O(TN^T)$。对于稍长的马尔可夫链来说，其联合分布的计算时间是难以忍受的。

2.前向算法

直接计算法通过枚举所有可能的状态序列来求解概率，由于每次计算新状态序列对应的概率时，会重复计算每个求和号中的结果，所以直接计算法存在大量重复的运算。

考虑一般的多级求和问题的求解：

$$\sum_{i=1}^N\sum_{j=1}^N\sum_{k=1}^Ni\times j\times k=\sum_{i=1}^N\left(i\sum_{j=1}^N\left(j\sum_{k=1}^Nk\right)\right)$$ $$\begin{align*} &\sum_{i=1}^N\sum_{j=1}^N\sum_{k=1}^Nf_0(i)\times f_1(i,j)\times f_2(j)\times f_3(j,k)\times f_4(k)\\ =&\sum_{k=1}^Nf_4(k)\times \left(\sum_{j=1}^Nf_3(j,k)\times f_2(j)\times\left(\sum_{i=1}^Nf_0(i)\times f_1(i,j)\right)\right) \end{align*}$$

每个靠外的求和问题都依赖于内部求和问题的求解。故每次求解时，可以保存子问题的求解结果，避免重复计算。

假设现在只有一个观测：

$$\begin{align*} P(O)&=P(o_1)\\ &=\sum_{I_1}P(O_1,I_1)\\ &=\sum_{I_1}P(O_1|I_1)P(I_1)\\ &=\sum_{i=1}^N P(O_1=o_1|I_1=i)P(I_1=i)\\ &=\sum_{i=1}^N \alpha_1(i) \end{align*}$$

两个观测：

$$\begin{align*} P(O)&=P(o_1,o_2)=\sum_{I_2}\sum_{I_1}P(O_1,O_2,I_1,I_2)\\ &=\sum_{I_2}\sum_{I_1}P(O_2|I_2)P(I_2|I_1)P(O_1|I_1)P(I_1)\\ &=\sum_{I_2}\left(P(O_2|I_2)\sum_{I_1}P(I_2|I_1)P(O_1|I_1)P(I_1)\right)\\ &=\sum_{I_2}\left(P(O_2|I_2)\sum_{i=1}^NP(I_2|I_1=i)\alpha_1(i)\right)\\ &=\sum_{j=1}^N\left(P(O_2=o_2|I_2=j)\sum_{i=1}^NP(I_2=j|I_1=i)\alpha_1(i)\right)\\ &=\sum_{j=1}^N\alpha_2(j) \end{align*}$$

三个观测：

$$\begin{align*} P(O)&=P(o_1,o_2,o_3)=\sum_{I_3}\sum_{I_2}\sum_{I_1}P(O_1,O_2,O_3,I_1,I_2,I_3)\\ &=\sum_{I_3}\sum_{I_2}\sum_{I_1}P(O_3|I_3)P(I_3|I_2)P(O_2|I_2)P(I_2|I_1)P(O_1|I_1)P(I_1)\\ &=\sum_{I_3}P(O_3|I_3)\left(\sum_{I_2}P(I_3|I_2)P(O_2|I_2)\left(\sum_{I_1}P(I_2|I_1)P(O_1|I_1)P(I_1)\right)\right)\\ &=\sum_{I_3}P(O_3|I_3)\left(\sum_{i=1}^NP(I_3|I_2=i)\alpha_2(i)\right)\\ &=\sum_{j=1}^NP(O_3=o_3|I_3=j)\left(\sum_{i=1}^NP(I_3=j|I_2=i)\alpha_2(i)\right)\\ &=\sum_{j=1}^N\alpha_3(j) \end{align*}$$

总结规律，得：

$$\begin{align*} \alpha_1(j)&=P(O_1=o_1|I_1=j)P(I_1=j)\\ \alpha_2(j)&=P(O_2=O_2|I_2=j)\left(\sum_{i=1}^NP(I_2=j|I_1=i)\alpha_1(i)\right)\\ \alpha_3(j)&=P(O_3=o_3|I_3=j)\left(\sum_{i=1}^NP(I_3=j|I_2=i)\alpha_2(i)\right)\\ \dots\\ \alpha_{t+1}(j)&=P(O_{t+1}=o_{t+1}|I_{t+1}=j)\sum_{i=1}^NP(I_{t+1}=j|I_t=i)\alpha_t(i)\\ \alpha_{t+1}(j)&=\left(\sum_{i=1}^N\alpha_t(i)a_{ij}\right)b_j(o_{t+1}) \end{align*}$$

有了上述推导之后，再来看概率计算问题：

对于原本时间复杂度很高的问题，可通过变量消除法来解决，首选计算出 $\alpha_1(i_1)$，消去 $I_1$，然后据此可得到 $\alpha_2(i_2)$，依次类推，从内到外，依次消除变量。

本质：sum-product过程

在 $t+1$ 时刻，状态取值为 $j$ 的概率 $\alpha_{t+1}(j)$ 等于在 $t$ 时刻状态取值为 $i$ 的概率 $\alpha_t(i)$ 乘以从状态 $i$ 转移到状态 $j$ 的概率 $a_{ij}$，因为在 $t$ 时刻状态有 $N$ 种可能，因此对所有状态求和（sum），然后乘以在 $t+1$ 时刻，状态为 $j$ 时产生观测 $o_{t+1}$ 的概率 $b_j(o_{t+1})$（product）。

从概率的角度看待前向算法

$$\begin{align*} \alpha_t(i)&=P(o_1,o_2,\dots,o_t,I_t=i)\\ \alpha_{t+1}(j)&=P(o_1,o_2,\dots,o_t,o_{t+1},I_{t+1}=j)\\ &=\sum_iP(o_1,o_2,\dots,o_t,I_t=i,o_{t+1},I_{t+1}=j)\\ &=\sum_i\left(P(o_1,o_2,\dots,o_t,I_t=i)\times P(I_{t+1}=j|o_1,o_2,\dots,o_t,I_t=i)\\\times P(o_{t+1}|o_1,o_2,\dots,o_t,I_t=i,I_{t+1}=j)\right)\\ &=\sum_{i}\left(\alpha_t(i)P(I_{t+1}=j|I_t=i)P(o_{t+1}|I_{t+1}=j)\right)\\ &=\left(\sum_i\alpha_t(i)a_{ij}\right)\times b_j(o_{t+1}) \end{align*}$$

公式第一行通过增加一个 $t$ 时刻的状态 $I_t=i$，在对 $i$ 求和与原式等价，并在第二行使用贝叶斯公式展开第一项变为前向概率，后续推导继续使用贝叶斯公式化简。值得注意的是，由于 $t+1$ 时刻的观测 $o_{t+1}$ 与前面的观测 $o_1\sim o_t$ 以及 $t$ 时刻的状态 $I_t=i$ 并没有关系，故后一项条件概率可化简为第三行的形式。

有了前面的分析，观测序列概率的前向算法如下：

一道例题

考虑盒子和球模型 $\lambda=(A,B,\Pi)$，状态集合 $Q=\{1,2,3\}$，观测集合 $V=\{$红, 白$\}$，

$$A=\left[\begin{aligned}&0.5&0.2&&0.3\\&0.3&0.5&&0.2\\&0.2 &0.3&&0.5\end{aligned}\right],B=\left[\begin{aligned}&0.5&0.5\\&0.4&0.6\\&0.7 &0.3\end{aligned}\right], \pi=(0.2\ \ 0.4\ \ 0.4)$$

设 $T=3$，$X=($红, 白, 红$)$，使用前向算法计算 $P(X|\lambda)$。

1. 初始值

   第一观测 $x_1=$红，则计算出这种观测来自的盒子的概率分布为：

   $$\begin{align*} \alpha_1(Y_1)&=\pi b_1(x_1=红)\\ &=(0.2\ 0.4\ 0.4)·(0.5\ 0.4\ 0.7)\\ &=(0.1\ 0.16\ 0.28) \end{align*}$$

2. 递推

   由前向公式向前递推计算到 $T=3$ 时的前向概率

   $$\begin{align*} \alpha_2(Y_2)&=b_2(x_2=白)\sum_{y_1}\alpha_1(Y_1)a_{y_1y_2}\\ &=(0.5\ 0.6\ 0.3)·\left(\begin{aligned}0.1\times 0.5+0.16\times0.3+0.28\times0.2\\0.1\times 0.2+0.16\times0.5+0.28\times0.3\\0.1\times 0.3+0.16\times0.2+0.28\times0.5\end{aligned}\right)^T\\ &=(0.5\ 0.6\ 0.3)·\left(\begin{aligned}0.154\\0.184\\0.202\end{aligned}\right)^T\\ &=(0.077\ 0.1104\ 0.0606) \end{align*}$$

​ 进一步计算第三时刻的前向概率：

$$\begin{align*} \alpha_3(Y_3)&=b_3(x_3=红)\sum_{y_2}\alpha_2(Y_2)a_{y_2y_3}\\ &=(0.5\ 0.4\ 0.7)·\left(\begin{aligned}0.077\times 0.5+0.1104\times0.3+0.0606\times0.2\\0.077\times 0.2+0.1104\times0.5+0.0606\times0.3\\0.077\times 0.3+0.1104\times0.2+0.0606\times0.5\end{aligned}\right)^T\\ &=(0.5\ 0.4\ 0.7)·\left(\begin{aligned}0.08374\\0.08878\\0.07548\end{aligned}\right)^T\\ &=(0.04187\ 0.035512\ 0.052836) \end{align*}$$

1. 停止

   计算得到边缘分布为：

   $$P(x_1,x_2,x_3;\lambda)=\sum_{y_3}\alpha_3(x_3)=0.130218$$

3.后向算法

有前向算法就有后向算法。后向算法通过从后往前逆向计算“状态路径的最优结构”实现边缘分布的求解如下图所示。

$$\begin{align*} P(O)&=P(O_1,O_2,\dots,O_T)\\ &=\sum_IP(O_1,O_2,\dots,O_T,I_1,I_2,\dots,I_T)\\ &=\sum_{I_1,\dots,I_T}P(I_1)P(O_1|I_1)P(I_2|I_1)P(O_2|I_2)\dots P(I_T|I_{T-1})P(O_T|I_T)\\ &=\sum_{I_1}P(I_1)P(O_1|I_1)\left(\sum_{I_2}P(I_2|I_1)P(O_2|I_2)\dots\left(\sum_{I_T}P(I_T|I_{T-1})P(O_T|I_T)1\right)\right)\\ &=\sum_{I_1}P(I_1)P(O_1|I_1)\beta_1(I_1)\\ &=\sum_{i=1}^N\pi_ib_i(o_1)\beta_1(i) \end{align*}$$

递推公式为：

$$\begin{align*} \beta_t(i)&=\sum_{j=1}^NP(I_{t+1}+j|I_t=i)P(O_{t+1}=o_{t+1}|I_{t+1}=j)\beta_{t+1}(I_{t+1}=j)\\ &=\sum_{j=1}^Na_{ij}b_j(o_{t+1})\beta_{t+1}(j) \end{align*}$$

其中，$i=1,2,\dots,N, t=T-1,T-2,\dots,1$ 。

4.前向-后向

$$\begin{align*} P(O|\lambda)&=\sum_{i=1}^N\sum_{j=1}^NP(o_1,\dots,o_t,I_t=i,I_{t+1}=j,o_{t+1},\dots,o_T)\\ &=\sum_{i=1}^N\sum_{j=1}^NP(o_1,\dots,o_t,I_t=i)\times P(I_{t+1}=j,o_{t+1},\dots,o_T|o_1,\dots,o_t,I_t=i)\\ &=\sum_{i=1}^N\sum_{j=1}^N\alpha_t(i)\times P(I_{t+1}=j,o_{t+1},\dots,o_T|I_t=i)\\ &=\sum_{i=1}^N\sum_{j=1}^N\alpha_t(i)\times P(I_{t+1}=j|I_t=i)\times P(o_{t+1},\dots,o_T|I_t=i,I_{t+1}=j)\\ &=\sum_{i=1}^N\sum_{j=1}^N\alpha_t(i)\times a_{ij}\times P(o_{t+1},o_{t+2},\dots,o_T|I_{t+1}=j)\\ &=\sum_{i=1}^N\sum_{j=1}^N\alpha_t(i)\times a_{ij}\times P(o_{t+1}|I_{t+1}=j)\times P(o_{t+2},\dots,o_T|I_{t+1}=j,o_{t+1})\\ &=\sum_{i=1}^N\sum_{j=1}^N\alpha_t(i)\times a_{ij}\times b_j(o_{t+1})\times\beta_{t+1}(j),\ t=1,2,\dots,T-1\\ &=\sum_{i=1}^N\sum_{j=1}^N\alpha_t(i)\times a_{ij}\times b_j(o_{t+1})\times\beta_{t+1}(j)\\ &=\sum_{j=1}^N\left(\sum_{i=1}^N\alpha_t(i)\times a_{ij} \times b_j(o_{t+1})\right)\times \beta_{t+1}(j)\\ &=\sum_{j=1}^N \alpha_{t+1}(j)\times\beta_{t+1}(j)=\sum_{i=1}^N\alpha_t(i)\times \beta_t(i) \end{align*}$$

• 当 $t=T-1$ 时，

  $$P(O|\lambda)=\sum_{j=1}^N\alpha_{t+1}(j)\times\beta_{t+1}(j)=\sum_{j=1}^N \alpha_T(j)\times \beta_T(j)=\sum_{i=1}^N \alpha_T(i)$$

• 当 $t=1$ 时，

  $$P(O|\lambda)=\sum_{i=1}^N\alpha_t(i)\times\beta_t(i)=\sum_{i=1}^N\alpha_1(i)\times\beta_1(i)=\sum_{i=1}^N\pi_ib_1(o_1)\times\beta_1(i)$$

另外几个重要的量

• $\gamma_t(i)=P(i_t=i|O,\lambda)=\frac{P(i_t=i,O|\lambda)}{P(O|\lambda)}=\frac{\alpha_t(i)\beta_t(i)}{\sum_{j=1}^N\alpha_t(j)\times\beta_t(j)}$
• $\xi_t(i,j)=P(i_t=i,i_{t+1}=j|O,\lambda)=\frac{P(i_t=i,i_{t+1}=j,O|\lambda)}{P(O|\lambda)}=\frac{\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)}{\sum_{i=1}^N\sum_{j=1}^N\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)}$
• $\gamma_t(i)=P(i_t=i|O,\lambda)=\sum_jP(i_t=i,i_{t+1}=j|O,\lambda)=\sum_j\xi_t(i,j)$

其他

• 在观测 $O$ 下，状态 $i$ 出现的期望值

  $$\sum_{t=1}^T\gamma_t(i)$$

• 在观测 $O$ 下，由状态 $i$ 转移的期望值

  $$\sum_{t=1}^{T-1}\gamma_t(i)$$

• 在观测 $O$ 下，由状态 $i$ 转移到状态 $j$ 的期望值

  $$\sum_{t=1}^{T-1}\xi_t(i,j)$$

预测（解码）问题

已知：模型 $\lambda=(A,B,\Pi)$ 和观测序列 $O=(o_1,o_2,\dots,o_T)$，求：$I^\ast=\text{arg max}_IP(I|O,\lambda)$。

$$I^\ast=\text{arg max}_I P(I|O,\lambda)=\text{arg max}_I \frac{P(I|\lambda)P(O|I,\lambda)}{P(O|\lambda)}=\text{arg max}_IP(I|\lambda)P(O|I,\lambda)$$

近似算法

$$\gamma_t(i)=P(i_t=i|O,\lambda)=\frac{\alpha_t(i)\beta_t(i)}{\sum_{j=1}^N\alpha_t(j)\times\beta_t(j)}\\ i_t^\ast=\text{arg max}_{1\le i\le N}[\gamma_t(i)],\ t=1,2,\dots,T$$

这存在两个问题：

1. 不能保证整体最优；

2. 存在 $a_{ij}=0$ 的状态。

Viterbi算法

• 假设现在只有一个观测：

$$\begin{align*} &\text{max}_{i\in\{1,2,\dots,N\}}P(I_1=i)\times P(O_1|I_1=i)\\ =&\text{max}_{i\in\{1,2,\dots,N\}}\delta_1(i)\\ =&\text{max}_{i\in\{1,2,\dots,N\}}\pi_ib_i(o_1) \end{align*}$$

• 两个观测：

  

  $$\begin{align*} &\text{max}_{I_2}\text{max}_{I_1}P(I_1)\times P(O_1|I_1)\times P(I_2|I_1)\times P(O_2|I_2)\\ =&\text{max}_{I_2}\left\{\text{max}_{I_1}\{P(I_1)\times P(O_1|I_1)\}\times P(I_2|I_1)\times P(O_2|I_2)\right\}\\ =&\text{max}_{I_2}\left\{\left[\text{max}_{I_1}\delta_1(i)\times P(I_2|I_1=i)\right]\times P(O_2|I_2)\right\}\\ =&\text{max}_{j\in\{1,2,\dots,N\}}\left\{\left[\text{max}_{I_1}\delta_1(i)\times P(I_2=j|I_1=i)\right]\times P(O_2=o_2|I_2=j)\right\}\\ =&\text{max}_{j\in\{1,2,\dots,N\}}\delta_2(j) \end{align*}$$ $$\delta_2(j)=\text{max}_{i\in\{1,2,\dots,N\}}[\delta_1(i)a_{ij}]\times b_j(o_2)$$

• 三个观测：

  

$$\delta_3(j)=\text{max}_{i\in\{1,2,\dots,N\}}[\delta_2(i)a_{ij}]\times b_j(o_3)$$

总的算法流程如下：

参数学习问题

已知：观测序列 $O=(o_1,o_2,\dots,o_T)$，求：模型 $\lambda=(A,B,\Pi)$。

根据一系列状态变量和观测变量的样本，求解HMM三要素。该问题可分为有监督学习和无监督学习。无监督学习需要使用EM算法求解，在之后的贝叶斯学习会涉及。下面是有监督学习问题的求解。

有监督学习，即观测序列和对应的状态序列已知。求解有监督学习问题常使用极大似然估计，该方法认为出现的观测往往对应着较高的概率，通过解析地表示出概率（似然函数），优化似然函数就可求得我们想要的某些概率。$D=\{(O^s, I^s\}_{s=1}^S,\lambda=(A,B,\pi)$，推导如下：

$$\begin{align*} L(D;\lambda)&=\text{ln}\prod_{s=1}^Sp(O^s,I^s)=\sum_{s=1}^S\text{ln}p(O^s,I^s)\\ &=\sum_{s=1}^S \text{ln}\left(\left(\prod_{t=1}^{T_s}p(O_t^s|I_t^s)\right)\times\left(\prod_{t=1}^{T_s-1}p(I_{t+1}^s|I_t^s)\right)\times p(I_1^s)\right)\\ &=\sum_{s=1}^S\sum_{t=1}^{T_s}\text{ln}p(O_t^s|I_t^s)+\sum_{s=1}^S\sum_{t=1}^{T_s-1}\text{ln}p(I_{t+1}^s|I_t^s)+\sum_{s=1}^S\text{ln}(p(I_1^s))\\ &=\sum_{s=1}^S\sum_{t=1}^{T_s}\sum_{j=1}^N\sum_{k=1}^M1\{I_t^s=j\wedge O_t^s=k\}\text{ln}b_{jk}+\sum_{s=1}^S\sum_{t=1}^{T_s-1}\sum_{i=1}^N\sum_{j=1}^N1\{I_t^s=i\wedge I_{t+1}^s=j\}\text{ln}a_{ij}\\&+\sum_{s=1}^S\sum_{i=1}^N1\{I_1^s=i\}\text{ln}(\pi_i) \end{align*}$$

上式为学习问题对应地似然函数，分为三部分，分别对应HMM中的三个参数。

对初始概率，可求得使似然函数最大的概率值（由上式第三项加上初始概率和为1的约束形成的拉格朗日函数）：

$$L(\pi,\alpha)=\sum_{s=1}^S\sum_{i=1}^N1\{I_t^s=i\}\text{ln}(\pi_i)+\alpha\left(1-\sum_{i=1}^N\pi_i\right)$$

对 $\pi_i$ 求偏导得：

$$\frac{\partial L}{\partial \pi_i}=\frac{1}{\pi_i}\sum_{s=1}^S1\{I_t^s=i\}-\alpha=0\Rightarrow\pi_i=\frac{\sum_{s=1}^S1\{I_t^s=i\}}{\alpha}$$

约束为：

$$\sum_{i=1}^N\pi_i=1 \Rightarrow\sum_{i=1}^N\left(\frac{\sum_{s=1}^S1\{I_t^s=i\}}{\alpha}\right)=\frac{\sum_{i=1}^N\sum_{s=1}^S1\{I_t^s=i\}}{\alpha}=1\Rightarrow\alpha=\sum_{i=1}^N\sum_{s=1}^S1\{I_t^s=i\}=S$$

所以，

$$\pi_i=\frac{\sum_{s=1}^S1\{I_1^s=i\}}{S}$$

直观来看，初始概率及统计初始状态中各个状态出现的频率。

同理，可推导转移概率矩阵得学习公式：

$$L(A,\beta)=\sum_{s=1}^S\sum_{t=1}^{T_s}\sum_{i=1}^N\sum_{i=1}^N1\{I_t^s=i\wedge I_{t+1}^s=j\}\text{ln}a_{ij}+\beta\left(1-\sum_{j=1}^Na_{ij}\right)$$

对 $a_{ij}$ 求偏导得：

$$\frac{\partial L}{\partial a_{ij}}=\sum_{s=1}^S\sum_{t=1}^{T_s}\frac{1\{I_t^s=1\wedge I_{t+1}^s=j\}}{a_{ij}}-\beta=0\Rightarrow a_{ij}=\frac{\sum_{s=1}^S\sum_{t=1}^{T_s-1}1\{I_t^s=1\wedge I_{t+1}^s=j\}}{\beta}$$

约束为：

$$\sum_{j=1}^Na_{ij}=1\Rightarrow\frac{\sum_{j=1}^N\sum_{s=1}^S\sum_{t=1}^{T_s-1}1\{I_t^s=i\wedge I_{t+1}^s=j\}}{\beta}=1$$

则，

$$\beta=\sum_{j=1}^N\sum_{s=1}^S\sum_{t=1}^{T_s-1}1\{I_t^s=i\wedge I_{t+1}^s=j\}=\sum_{s=1}^S\sum_{t=1}^{T_s-1}1\{I_t^s=i\}$$

代入，

$$a_{ij}=\frac{\sum_{s=1}^S\sum_{t=1}^{T_s-1}1\{I_t^s=i\wedge I_{t+1}^s=j\}}{\sum_{s=1}^S\sum_{t=1}^{T_s-1}1\{I_t^s=i\}}$$

直观解释就是，根据样本中不同状态之间的转移数目，计算某两个状态转移的频率，将其作为概率。

观测概率的学习公式为：

$$L(B,\gamma)=\sum_{s=1}^S\sum_{t=1}^{T_s}\sum_{i=1}^N\sum_{k=1}^M1\{I_t^s=j\wedge O_t^s=k\}\text{ln}b_{jk}+\gamma(1-\sum_{k=1}^Mb_{jk})$$

对 $b_{jk}$ 求偏导得：

$$\frac{\partial L}{\partial b_{jk}}=\sum_{s=1}^S\sum_{t=1}^{T_s}\frac{1\{I_t^s=j\wedge O_t^s=k\}}{b_{jk}}-\gamma=0\Rightarrow b_{jk}=\frac{\sum_{s=1}^S\sum_{t=1}^{T_s}1\{I_t^s=j\wedge O_t^s=k\}}{\gamma}$$

约束为：

$$\sum_{k=1}^Mb_{jk}=1\Rightarrow\frac{\sum_{k=1}^M\sum_{s=1}^S\sum_{t=1}^{T_s}1\{I_t^s=j\wedge O_t^s=k\}}{\gamma}=1$$

则，

$$\gamma=\sum_{k=1}^M\sum_{s=1}^S\sum_{t=1}^{T_s}1\{I_t^s=j\wedge O_t^s=k\}=\sum_{s=1}^S\sum_{t=1}^{T_s}1\{I_t^s=j\}$$

代入，

$$b_{jk}=\frac{\sum_{s=1}^S\sum_{t=1}^{T_s}1\{I_t^s=j\wedge O_t^s=k\}}{\sum_{s=1}^S\sum_{t=1}^{T_s}1\{I_t^s=j\}}$$

相似地，这可以被简单地归纳为统计频率的问题。这也比较符合通过频率逼近概率的思想。

贝叶斯学习

最大似然估计将待估计的参数看作是确定性的量，只是其值我们暂时还不知道；而贝叶斯估计则将待估计的参数看成是符合某种先验概率分布的随机变量。对样本进行观测的过程，就是把先验概率密度转化为后验概率密度，这样就利用样本的信息修正了对参数的初始估计值。这个过程也可以看成是贝叶斯学习的过程。

MLE

最大似然估计，用于估计模型参数。似然函数取最大时的参数值作为估计值。估计过程如下：

1. 假设有一批数据 $X$，且 $x$ 均服从某种参数为 $\theta$ 的分布：

   $$X=\{x_1,x_2,\dots,x_n\},\ x\sim P_\theta, p_\theta(x)=p(x|\theta)$$

2. 因为 $x$ 服从独立同分布，参数向量 $\theta$ 的MLE就是使 $P(X|\theta)$ 达到最大值的那个参数向量：

   $$p(X|\theta)=\prod_{x\in X}p(X=x|\theta)=\prod_{i=1}^n p(x_{i}|\theta)$$

3. 取对数：

   $$L(\theta|X)\equiv\text{log}p(X|\theta)=\sum_{t=1}^n\text{log}p(x_i|\theta)$$ $$\hat \theta=\text{arg max}_\theta L(\theta|X)=\text{arg max}_\theta \sum_{t=1}^n \text{log}p(x_i|\theta)$$

举例：伯努利分布

以抛硬币为例，正面朝上概率为 $p$，反面为 $1-p$。

参数估计：

$$\begin{align*} \theta^\ast&=\text{arg max}_\theta\sum_{t=1}^n \log p(x_i|\theta)\\ &=\text{arg max}_\theta\sum_{i=1}^n x_i\log p+(1-x_i)\log (1-p) \end{align*}$$

求导：

$$L(p)=\left(\sum_{i=1}^nx_i\right)\log p+\left(\sum_{i=1}^n(1-x_i)\right)\log (1-p)\\ \frac{dL(p)}{dp}=0 \Rightarrow p=\frac{\sum_{i=1}^nx_i}{n}=\frac{\sharp head}{n}$$

假设抛一枚硬币 $n$ 次，其中正面朝上的次数用 $\sharp head$ 表示，求正面朝上的概率为多少？很多人都会计算：$\sharp head/n$。

从上面的推导也可以看到正面朝上的概率 $p$ 通过MLE方法得到的值与我们日常中估计的方法一致。

举例：多项式分布

以抛骰子为例，

$$x\in\{1,2,\dots,K\},\ p(x=i)=q_i,\ \sum_{i=1}^Kq_i=1$$

则有：

$$p(x_i)=\prod_{j=1}^K(q_j)^{1\{x_i=j\}}\\ \textbf x_i=\{0,0,0,1,0,0\}\\ p(\textbf x_i)=\prod_{j=1}^K q_j^{x_{ij}}$$

取对数：

$$D=\{\textbf x_1, \textbf x_2,\dots, \textbf x_n\},\\ p(D)=\prod_{i=1}^n\prod_{j=1}^K=q_j^{x_{ij}}=\prod_{j=1}^Kq_j^{\sum_i x_{ij}}=\prod_{j=1}^Kq_j^{n_j},\\ \ln p(D)=\sum_{j=1}^Kn_j\ln q_j$$

拉格朗日函数为：

$$L(q,\lambda)=\sum_{j=1}^Kn_j\ln q_j+\lambda(1-\sum_{j=1}^Kq_j)$$

求导：

$$\frac{\partial L(q,\lambda)}{\partial q_j}=\frac{n_j}{q_j}-\lambda=0\Rightarrow q_j=\frac{n_j}{\lambda},\\ \sum_{j=1}^Kq_j=1\Rightarrow \sum_{j=1}^K\frac{n_j}{\lambda}=1\Rightarrow \lambda=\sum_{j=1}^Kn_j=n,\\ q_j=\frac{n_j}{n}$$

举例：高斯分布

已知：$D=\{x_1,x_2,\dots,x_n\}, x\sim N(\mu,\sigma^2)$，求 $\mu$。

假设有一批数据 $D$，且 $x$ 服从均值为 $\mu$，方差为 $\sigma^2$ 的高斯分布。数据 $D$ 是如下图 $x$ 轴坐标所示，那么哪一个模型比较好的刻画数据 $D$ 的分布呢？

很显然，中间的红色的曲线比较能够豪地刻画 $D$ 的分布，下面从MLE的角度分析一下。

取对数：

$$\begin{align*} \text{ln}L(D;\mu)&=\ln\prod_{i=1}^n\frac{1}{\sqrt{2\pi}\sigma}\text{exp}\left[-\frac{1}{2\sigma^2}(x_i-\mu)^2\right]\\ &=\sum_{i=1}^n\left[\ln\left(\frac{1}{\sqrt{2\pi}\sigma}\right)-\frac{1}{2\sigma^2}(x_i-\mu)^2\right] \end{align*}$$

求导：

$$\frac{d}{d\mu}\ln L(D;\mu)=\frac{1}{\sigma^2}\sum_{i=1}^n(x_i-\mu)=0 \Rightarrow \mu=\frac{1}{n}\sum_{i=1}^nx_i$$

得到的结果为所有样本的均值，对应于上图，即高斯分布的中轴线的位置，差不多红色曲线中轴线的位置与数据均值的位置最为相近。当然，此刻得到的参数只是对于真实值的一个估计，其对于真实值的接近程度是受训练样本数量限制的。如果训练样本数量越多，那么对于参数的估计值也就越接近于真实值。

MAP

最大后验估计。

$$\begin{align*} \theta_{MAP}&=\text{arg max}_\theta p(\theta|X)\\ &=\text{arg max}_\theta\frac{p(X|\theta)p(\theta)}{p(X)}\\ &=\text{arg max}_\theta p(X|\theta)p(\theta) \end{align*}\\ p_{MAP}(x|X)=p(x|\theta_{MAP})$$

其中 $p(\theta)$ 为参数 $\theta$ 的先验概率。即上述贝叶斯过程可以这样描述：

$$后验=\frac{样本似然\times先验}{证据}$$

所以，MLE可以说是当先验概率 $p(\theta)$ 为均匀分布时的MAP。

举例：单变量高斯分布的均值

简单起见，只考虑只有均值 $\mu$ 未知的情况，且参数 $\mu$ 也服从高斯分布（均值、方差已知）。

$$D=\{x_1,\dots,x_n\},\ x_i\in R,\ X\sim N(\mu,\sigma^2),\ \mu\sim N(\mu_0, 1)$$

采用MAP：

$$\begin{align*} \mu_{MAP}&=\text{arg max}_\mu \text{ln}[p(D|\mu)p(\mu)]\\ &=\text{arg max}_\mu \text{ln}\left[\prod_{i=1}^n\frac{1}{\sqrt{2\pi}\sigma}\text{exp}\left(-\frac{1}{2\sigma^2}(x_i-\mu)^2\right)\right]+\text{ln}\frac{1}{\sqrt{2\pi}}\text{exp}\left(-\frac{1}{2}(\mu-\mu_0)^2\right)\\ &=\text{arg max}_\mu \left\{-\frac{1}{2\sigma^2}\sum_{i=1}^n(x_i-\mu)^2-\frac{1}{2}(\mu-\mu_0)^2\right\}\\ &=L(\mu) \end{align*}$$

求导：

$$\begin{align*} &\frac{\partial}{\partial \mu}L(\mu)=0\Rightarrow \frac{1}{\sigma^2}\left(\sum_{i=1}^nx_i-n\mu\right)-(\mu-\mu_0)=0\\ \Rightarrow&\mu=\frac{\frac{1}{\sigma^2}\sum_{i=1}^nx_i+\mu_0}{\frac{n}{\sigma^2}+1}=\frac{n\mu_{MLE}+\sigma^2\mu_0}{n+\sigma^2}=\frac{n}{n+\sigma^2}\mu_{MLE}+\frac{\sigma^2}{n+\sigma^2}\mu_0 \end{align*}$$

可以看到求得的 $\mu$ 为MLE得到的 $\mu$ 以及 $\mu_0$ 的加权和，且权值之和为1。这种组合称为凸组合，当 $n$ 很大的时候，即训练样本充足，$\mu$ 的取值与 $\mu_{MLE}$ 接近，这说明此时采用极大似然估计比较可靠；反之，当 $n$ 很小时，参数原始的分布更可靠一些。

举例：回归模型

$$D=\{\langle x_i,y_i\rangle\}_{i=1}^n,\ x_i\in\mathbb R^d,\ y_i\in \mathbb R,\ y_i=w^Tx_i+\varepsilon, \varepsilon\sim N(0,1)\\ p(y_i)\sim N(w^Tx,1),\ w\sim N(0, \lambda^{-1}I)$$

其中，$\varepsilon$ 服从均值为0，方差为1的高斯分布，则 $p(y)$ 服从均值为 $w^Tx_i$，方差为1的高斯分布，且参数 $w$ 也服从一个高斯分布。

取对数：

$$\begin{align*} &\text{ln}p(D|w)=\text{ln}\prod_{i=1}^np(y_i)=\sum_{i=1}^n\left[\left(\text{ln}\frac{1}{\sqrt{2\pi}}\right)-\frac{1}{2}(y_i-w^Tx_i)^2\right]\\ &\text{ln}p(w)=-\frac{d}{2}\text{ln}(2\pi)+\frac{1}{2}\text{ln}\lambda-\frac{\lambda}{2}w^Tw\end{align*}$$

求导：

$$W_{MAP}=\text{arg max}_w[\text{ln}p(D|w)+\text{ln}p(w)]\propto w^TX^TXw-2w^TX^Ty+\lambda w^Tw\\ \Rightarrow X^TXw-X^Ty+\lambda w=0\Rightarrow w=(X^TX+\lambda I)^{-1}X^Ty$$

可以看到整个优化过程，从贝叶斯的角度来看是一种贝叶斯最大后验估计，正则化项变成了一种先验信息。

贝叶斯估计

在贝叶斯学习方法中，我们把参数向量 $\theta$ 本身看成时一个随机变量，已有的训练样本使我们能够对于 $\theta$ 的初始的密度估计转化为后验概率密度。

$$p(\theta|X)=\frac{p(X|\theta)p(\theta)}{p(X)}=\frac{p(X|\theta)p(\theta)}{\int p(X|\theta')p(\theta')d\theta'}$$

• 先验信息：看到样本前关于参数的信息；
• 样本似然：当分布的参数确定下来后，看到这些样本的可能性。

用概率图模型表示：

已知一组训练数据 $X$，这些样本都是从固定但未知的概率密度函数 $p(\theta)$ 中独立抽取的，要求根据这些样本估计 $p(x,X,\theta)$ ，这就是贝叶斯学习的核心问题。

问题也可以看成，用已知的训练集来预测未知的测试集，模型参数只是一个中间变量，即

$$\begin{align*} p(x|X)&=\frac{p(x,X)}{p(X)}=\frac{\int_\theta p(x,X,\theta)d\theta}{p(X)}\\ &=\frac{\int_\theta p(\theta)p(X|\theta)p(x|\theta)d\theta}{p(X)}\\ &=\int_\theta p(\theta|X)p(x|\theta)d\theta \end{align*}$$

其中，$p(\theta|X)=\frac{p(\theta)p(X|\theta)}{p(X)}$，而贝叶斯估计的方法是：

$$\theta_{Bayes}=E[\theta|X]=\int\theta p(\theta|X)d\theta.$$

由此也可以看出贝叶斯估计的优点：

• 贝叶斯推理需要一个先验假设；
• 最小化期望损失：多参数求平均；
• 方便进行模型选择；
• 避免过度拟合。

和缺点：

• 必须有一个先验假设；
• 准确计算非常困难。

三者的关系如下：

MLE	MAP	Bayesian
$\theta_{MLE}=\text{arg max}_{\theta}p(X	\theta)$	$\theta_{MAP}=\text{arg max}_\theta p(X	\theta)p(\theta)$	$\theta_{Bayes}=E[\theta	X]=\int \theta p(\theta	X)d\theta$

由MLE可以计算MAP，由MAP可以进行贝叶斯估计，MAP可以看成是贝叶斯估计中参数最大值时的一个特例，贝叶斯估计则考虑了参数取值的所有情况，换句话说，在进行贝叶斯估计的时候就进行了模型选择。

举例：估计分布的参数：离散变量（多值）

$$x_{ij}=\left\{\begin{align}&1, &&\text{if }x_i \text{ is in state }j\\&0,&&\forall j\neq i\end{align}\right.,\ \sum_{j=1}^K q_j=1, \ q=[q_1,q_2,\dots,q_K]^T$$

写出似然函数：

$$p(X|q)=\prod_{i=1}^n\prod_{j=1}^K(q_j)^{x_{ij}}=\sum_{j=1}^K(q_j)^{\sum_i x_{ij}}$$

而贝叶斯估计需要一个关于参数的先验假设。这个先验假设成什么分布好呢？

一些人，应该说一大堆人，到了这个地方，老是跟你说“刚好啊”，“凑巧啊”，“不妨啊”，“这个分布具有某些良好的性质啊”。就是没说最根本的原因，坑了另外一大堆人，在这里搞得半懂不懂的。

这里，涉及到一个重要的概念，共轭分布，定义为：$p(w)$ 与 $p(w|D)$ 同分布。就是说，它使得后验概率分布的函数形式与先验概率相同，因此使得贝叶斯分析得到了极大的简化。

简单来说，按照共轭分布去找，这样找下来有几个对应的：

• 伯努利分布的先验为Beta分布；
• 多项分布的先验为Dirichlet分布；
• 高斯分布的先验为高斯分布。

插一句：推荐马春鹏同学翻译的PRML。

Dirichlet分布

事先说明一下，这个分布实在是，老母猪戴胸罩——一套又一套。

Dirichlet分布如下：

$$\theta\sim Dirichlet(\theta|\alpha)=\frac{1}{B(\alpha)}\prod_{j=1}^n\theta_j^{\alpha_j-1}$$

其中：

$$\theta=(\theta_1,\dots,\theta_n), \ \alpha=(\alpha_1,\dots,\alpha_n)$$

$B(\alpha)$ 为Beta函数，其定义为：

$$\frac{1}{B(\alpha)}=\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)\dots\Gamma(\alpha_n)},\ \alpha_0=\sum_{j=1}^n\alpha_j$$

Gamma函数定义如下：

$$\Gamma(x)=\int_{0}^\infty u^{x-1}e^{-u}du$$

如果参数 $\theta$ 满足Dirichlet分布，考虑所有的参数 $\theta$ 的分布情况，将其代入公式（可以看作将所有的 $\theta$，根据其概率密度函数，给予其相应权重，得到的全局期望值）：

$$E[\theta_i]=\int \theta_i p(\theta)d\theta=\frac{\alpha_i}{\alpha_0}$$

公式推导过程如下：

$$\begin{align*} E[\theta_1]&=\int \theta_1p(\theta)d\theta\\ &=\int\theta_1\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)\dots\Gamma(\alpha_n)}\prod_{j=1}^n\theta_j^{\alpha_j-1}d\theta\\ &=\int\frac{\alpha_1}{\alpha_0}\frac{\Gamma(\alpha_0+1)}{\Gamma(\alpha_1+1)\dots\Gamma(\alpha_n)}\theta_1^{1+(\alpha_1-1)}\prod_{j=2}^n\theta_j^{\alpha_j-1}d\theta\\ &=\frac{\alpha_1}{\alpha_0} \end{align*}$$

PPT上这里讲得太粗糙了，关于最后一步，需要一些性质。

首先是Gamma函数的性质：

$$\Gamma(x+1)=x\Gamma(x)$$

所以有：

$$\frac{\Gamma(\alpha_0+1)}{\Gamma(\alpha_1+1)\dots\Gamma(\alpha_n)}=\frac{\alpha_0\Gamma(\alpha_0)}{\alpha_1\Gamma(\alpha_1)\dots\Gamma(\alpha_n)}$$

然后是概率密度函数的性质。概率密度函数的积分为1，所以Dirichlet分布有：

$$\int \frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)\dots\Gamma(\alpha_n)}\prod_{j=1}^n\theta_j^{\alpha_j-1}d\theta=1$$

把公式中的右半部分看成一个新的Dirichlet分布，即按以下对应：

$$\alpha_1'=\alpha_1+1\\ \alpha_0'=\alpha_0+1$$

所以有：

$$\int\frac{\Gamma(\alpha_0+1)}{\Gamma(\alpha_1+1)\dots\Gamma(\alpha_n)}\theta_1^{1+(\alpha_1-1)}\prod_{j=2}^n\theta_j^{\alpha_j-1}d\theta=1$$

这样就ok了。

---

现在将已知的满足Dirichlet分布的 $\theta$ 信息代入到贝叶斯公式中，得到后验概率公式：

$$\begin{align*} p(\theta|X)&\propto p(X|\theta)p(\theta|\alpha)\propto \prod_{j=1}^K\theta_j^{\sum_ix_{ij}}\times\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)\dots\Gamma(\alpha_K)}\prod_{j=1}^K\theta_j^{\alpha_j-1}\\ &=\frac{\Gamma(\alpha_0)}{\Gamma(\alpha_1)\dots\Gamma(\alpha_K)}\prod_{j=1}^K\theta_j^{\alpha_j+N_j-1}\propto \frac{\Gamma(\alpha_0+N)}{\Gamma(\alpha_1+N_1)\dots\Gamma(\alpha_K+N_K)}\prod_{j=1}^K\theta_j^{\alpha_j+N_j-1}\\ &=Dirichlet(\theta|\alpha+\textbf n) \end{align*}$$

其中，$\textbf n=[N_1,N_2,\dots,N_K]^T, N=\sum_i N_i, N_i$ 为 $x_i$ 出现的次数。这里能够看出结果是一个新的Dirichlet分布。这也是我们使用Dirichlet分布作为多项式分布的先验的原因，强大的聚合性。

之后根据Dirichlet分布的性质，有：

$$p(x|X)=\int_\theta p(x|\theta)p(\theta|X)d\theta=\int_\theta\theta_xp(\theta|X)d\theta=E[\theta_x|X]=\frac{n_x+a_x}{n+a_0}$$

通过上式可以很明显地看出，随着已知信息（先验概率）的不断更新，其对后验概率造成的直接影响。

举例：估计分布的参数：离散变量（二值）

$X=\{x_i\}_{i=1}^n, x_i\in\{0,1\}, X\sim Bern(\theta)$，其样本似然为：

$$p(X|\theta)=\prod_{i=1}^n\theta^{x_i}(1-\theta)^{1-x_i}=\theta^{n_1}(1-\theta)^{n_0},\ n_1=\sum_ix_i,\ n_0=n-n_1.$$

当参数 $\theta$ 服从Beta分布，Beta分布为：

$$Beta(\theta|a,b)=\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\theta^{a-1}(1-\theta)^{b-1}$$

Beta分布 $\theta \sim Beta(\theta|a,b)$ 有以下性质：

• $E[\theta]=\frac{a}{a+b}$；
• $\sigma^2(\theta)=\frac{ab}{(a+b)^2(a+b+1)}$；
• $\textbf{mode}=\frac{a-1}{a+b-2}$。这里的mode是指最大值，概念来源于聚类，一个高峰即一个类别即一个mode

同样的，先计算其后验概率：

$$\begin{align*} p(\theta|X)&\propto p(\theta)p(X|\theta)=\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\theta^{a-1}(1-\theta)^{b-1}\times\theta^{n_1}(1-\theta)^{n_0}\\ &=\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\theta^{n_1+a-1}(1-\theta)^{n_0+b-1}\\ &=Beta(\theta|n_1+a,n_0+b) \end{align*}$$

代入概率图模型得到：

$$\begin{align*} p(x=1|X)&=\int p(x=1|\theta)p(\theta|X)d\theta\\ &=\int \theta Beta(\theta|a+n_1,b+n_0)d\theta\\ &=E[\theta]=\frac{a+n_1}{a+b+n} \end{align*}$$

二项分布的MLE/MAP/贝叶斯估计三者比较

$$\begin{align*} &\theta_{MLE}=\text{arg max }p(X|\theta)\Rightarrow\theta_{MLE}=\frac{n_1}{n}\\ &\theta_{MAP}=\text{arg max }p(\theta|X)\Rightarrow\theta_{MAP}=\frac{a+n_1-1}{a+b+n-2}\\ &\theta_{Bayes}=\int p(x=1|\theta)p(\theta|X)d\theta\Rightarrow\theta_{Bayes}=\frac{a+n_1}{a+b+n} \end{align*}$$

通过比较我们可以看出在MLE中，对于概率的计算偏向于频率学派的观点，而MAP和Bayes中，充分吸收了Bayes学派的观点，考虑了先验信息。当已知信息有限时，MAP和Bayes的方法更加科学和有效。

其中Bayes的结果还可以写成：

$$\frac{a+n_1}{a+b+n}=\frac{a+b}{a+b+n}\times\frac{a}{a+b}+\frac{n}{a+b+n}\times\frac{n_1}{n}$$

通过上式可以更好地说明，当 $n$ （试验次数）充分大时，Bayes的结果与MLE相等。当 $n$ 有限时，则充分考虑先验信息。

估计分布的参数：单变量高斯，方差已知

$X=\{x_1,\dots,x_n\}, x_i\in \mathbb R, X\sim N(\mu, \sigma^2), \mu\sim N(\mu_0,\sigma_0^2)$，这里的似然为：

$$p(X|\mu)=\prod_{i=1}^n\frac{1}{\sqrt{2\pi}\sigma}\text{exp}\left[-\frac{(x_i-\mu)^2}{2\sigma^2}\right]=\left(\frac{1}{\sqrt{2\pi\sigma}}\right)^n\text{exp}\left[-\frac{1}{2\sigma^2}\sum_{i=1}^n(x_i-\mu)^2\right]$$

$\mu$ 的先验为：

$$p(\mu)\sim N(\mu_0,\sigma_0^2)=\frac{1}{\sqrt{2\pi}\sigma_0}\text{exp}\left[-\frac{1}{2\sigma_0^2}(\mu-\mu_0)^2\right]$$

根据后验概率公式：

$$\begin{align*} p(\mu|X)&\propto p(\mu)p(X|\mu)\\ &=\left\{\left(\frac{1}{\sqrt{2\pi}\sigma}\right)^n\text{exp}\left[-\frac{1}{2\sigma^2}\sum_{i=1}^n(x_i-\mu)^2\right]\right\}\times\left\{\frac{1}{\sqrt{2\pi}\sigma_0}\text{exp}\left[-\frac{(\mu-\mu_0)^2}{2\sigma_0^2}\right]\right\} \end{align*}$$

考虑概率密度函数的等比例性质：

$$f=\frac{g(x)}{\int g(x)dx}$$

其中等式左边为概率密度函数，证明如下：

$$\begin{align*} &f\propto g\Rightarrow g(x)=cf(x)\\ &\frac{g(x)}{\int g(x)dx}=\frac{cf(x)}{\int cf(x)dx}=f(x) \end{align*}$$

这里将概率密度函数的等比例性质用于从先验分布×似然分布=后验分布的过程中，多余的系数可以作为积分项省略，包括前面的Dirichlet分布和Beta分布。

即后验概率函数又可写成：

$$\begin{align*} p(\mu|X)&\propto p(\mu)p(X|\mu)\\ &=\left\{\left(\frac{1}{\sqrt{2\pi}\sigma}\right)^n\text{exp}\left[-\frac{1}{2\sigma^2}\sum_{i=1}^n(x_i-\mu)^2\right]\right\}\times\left\{\frac{1}{\sqrt{2\pi}\sigma_0}\text{exp}\left[-\frac{(\mu-\mu_0)^2}{2\sigma_0^2}\right]\right\}\\ &\propto \text{exp}\left[-\frac{1}{2\sigma^2}\sum_{i=1}^n(x_i-\mu)^2-\frac{1}{2\sigma_0^2}(\mu-\mu_0)^2\right]\\ &=\text{exp}\left[(\dots)\mu^2+(\dots)\mu+c\right]\\ &=\text{exp}\left[\left(-\frac{n}{2\sigma^2}-\frac{1}{2\sigma_0^2}\right)\mu^2+\left(\frac{\sum_i x_i}{\sigma^2}+\frac{\mu_0}{\sigma_0^2}\right)\mu+c\right]\\ &=\text{exp}\left[-\frac{1}{2}\left(\frac{n}{\sigma^2}+\frac{1}{\sigma_0^2}\right)\left(\mu-\frac{\frac{\sum_ix_i}{\sigma^2}+\frac{\mu_0}{\sigma_0^2}}{\frac{n}{\sigma^2}+\frac{1}{\sigma_0^2}}\right)^2+c'\right]\Leftarrow ax^2+bx+c=a\left(x+\frac{b}{2a}\right)^2+\frac{-V}{4ac}\\ &\propto \text{exp}\left[-\frac{1}{2\sigma_n^2}(\mu-\mu_n)^2\right],\ \frac{1}{\sigma_n^2}=\frac{1}{\sigma_0^2}+\frac{n}{\sigma^2},\ \mu_n=\frac{\sigma^2}{n\sigma_0^2+\sigma^2}\mu_0+\frac{n\sigma_0^2}{n\sigma_0^2+\sigma^2}\mu_{MLE} \end{align*}$$

从上式可知，先验信息（高斯分布的均值）会随着试验次数 $n$ 的增大，对后验信息的影响逐渐减小。

即后验的分布为：

$$p(\mu|X)=N(\mu_n,\sigma_n^2)$$

代入概率图为：

$$p(x|X)=\int p(x|\mu)p(\mu|X)d\mu=E_{\mu|X}[p(x|\mu)]$$

随着试验次数的增多，预测结果的方差会减小，即预测结果更加准确。

估计分布的参数：单变量高斯，方差未知

已知：

$$X=\{x_1,\dots,x_n\},\ x_i\in \mathbb R,\ X\sim N(\mu,\sigma^2),\ \textbf{precision:}\lambda=\frac{1}{\sigma^2}$$

且其中 $\lambda$ 满足 $\Gamma$ 函数分布：

$$p(\lambda)\sim Gamma(\lambda|a_0,b_0)=\frac{b_0^{a_0}\lambda^{a_0-1}\text{exp}(-b_0\lambda)}{\Gamma(a_0)}$$

似然函数为：

$$p(X|\lambda)=\prod_{i=1}^n\frac{\lambda^{1/2}}{\sqrt{2\pi}}\text{exp}\left[-\frac{\lambda(x_i-\mu)^2}{2}\right]=\lambda^{n/2}(2\pi)^{-n/2}\text{exp}\left[-\frac{\lambda}{2}\sum_{i=1}^n(x_i-\mu)^2\right]$$

计算后验概率：

$$\begin{align*} p(\lambda|X)&\propto p(\lambda)p(X|\lambda)\\ &\propto(\lambda^{a_0-1}\text{exp}(-b_0\lambda))\times\left(\lambda^{n/2}\text{exp}\left[-\frac{\lambda}{2}\sum_{i=1}^n(x_i-\mu)^2\right]\right)\\ &=\lambda^{a_0+n/2-1}\text{exp}\left\{-\left[b_0+\frac{1}{2}\sum_{i=1}^n(x_i-\mu)^2\right]\lambda\right\}\\ &\sim Gamma(a_n,b_n),\ a_n=a_0+\frac{n}{2},\ b_n=b_0+\frac{n}{2}s^2,\ s^2=\frac{1}{n}\sum_{i=1}^n(x_i-\mu)^2 \end{align*}$$

显然后验概率满足Gamma函数，其参数 $(a_n,b_n)$ 与试验次数相关。

同理，有：

$$\begin{align*} &p(\lambda|X)=Gamma(a_n,b_n)\\ &p(x|X)=\int p(x|\lambda)p(\lambda|X)d\lambda=E_{\lambda|X}[p(x|\lambda)] \end{align*}$$

估计模型的参数：回归

贝叶斯学习方法同样可以用在回归模型中，问题描述为：

$$D=\{\langle x_i,y_i\rangle\}_{i=1}^n,\ y=w^Tx+\varepsilon,\ \varepsilon\sim N(0,\beta^{-1})$$

则：

$$p(y|x,w)\sim N(w^Tx,\beta^{-1})$$

回顾最大似然估计，其目标函数为：

$$\begin{align*} L(X|w)&=\text{log}p(X|w)=\text{log}p(X,y|w)=\text{log}p(y|X,w)+\text{log}p(X) \end{align*}$$

其中：

$$\begin{align*} \text{log}p(y|X,w)&=\text{log}\prod_{i=1}^np(y_i|x_i,w)\\ &=\text{log}\prod_{i=1}^n\sqrt\frac{\beta}{ {2\pi} }\text{exp}\left[-\frac{\beta}{2}(y_i-w^Tx_i)^2\right]\\ &=-n\text{log}\sqrt{2\pi}+n\text{log}\sqrt \beta -\frac{\beta}{2}\sum_{i=1}^n(y_i-w^Tx_i)^2 \end{align*}$$

对参数的估计可转化为最小二乘法：

$$w^\ast=\text{arg max}_w\text{log}p(y|X,w)\Leftrightarrow\text{arg min}_w\frac{1}{2}\sum_{i=1}^n(y_i-w^Tx_i)^2$$

结果为：

$$w_{MLE}=(X^TX)^{-1}X^Ty$$

利用贝叶斯的逻辑，将充分考虑先验信息（结构风险最小化）：

$$p(w|X)\propto p(w)p(X|w)\propto \text{exp}\left[-\frac{1}{2}\left(\beta(y-Xw)^T(y-Xw)+\alpha w^Tw\right)\right]$$

接下来的公式推导目的是，后验分布为高斯分布：

$$\begin{align*} &\beta(y-Xw)^T(y-Xw)+\alpha w^Tw\\ =&w^T(\beta X^TX+\alpha I)w-2\beta w^TX^T y+\beta y^Ty \end{align*}$$

设 $\Sigma_n^{-1}=(\beta X^TX+\alpha I), \Sigma_n^{-1}\mu_n=\beta X^Ty$，则：

$$\begin{align*} &w^T(\beta X^TX+\alpha I)w-2\beta w^TX^T y+\beta y^Ty\\ =&(w-\mu_n)^T\Sigma_n^{-1}(w-\mu_n)+常数\sim N(\mu_n,\Sigma_n) \end{align*}$$

其中，$\mu_n=\beta\Sigma_nX^Ty, \Sigma_n=(\beta X^TX+\alpha I)^{-1}$。

以上推导需要 $\beta X^TX+\alpha I$ 可逆，证明：

因为 $X^TX$ 可逆，所以有：

$$\begin{align*} &(X^TX)u=\lambda u\Rightarrow(\beta X^TX+\alpha I)u=(\lambda \beta+\alpha)u,\\ &\lambda \ge 0, \beta \ge 0, \alpha \ge 0 \Rightarrow \lambda \beta+\alpha >0\Rightarrow \beta X^TX+\alpha I可逆. \end{align*}$$

考虑参数的所有取值情况，即代入概率图进行预测，有：

$$\begin{align*} y&=\int(w^Tx)p(w|X)dw\\ &=\left[\int wp(w|X)dw\right]^Tx\\ &=\left(E_{w|X}[w]\right)^Tx\\ &=\mu_n^Tx \end{align*}$$

同样先考虑参数的所有取值情况，由于参数符合高斯分布，且概率之和等于一，有：

$$\begin{align*} p(y|x,X)&=\int_w p(y,w|x,X)dw\\ &=\int_w p(y|x,w,X)p(w|x,X)dw\\ &=\int_w p(y|x,w)p(w|X)dw\\ &=\int_w p(y|w^Tx,\beta^{-1})N(w|\mu_n,\Sigma_n)dw\\ &=\int_w g(y)N(w|\mu,\Sigma)dw\Leftarrow我们的目标\\ &=g(y) \end{align*}$$

将两个高斯分布展开：

$$\begin{align*} p(y|x,X)&=\int_w N(y|w^Tx,\beta^{-1})N(w|\mu_n,\Sigma_n)dw\\ &\propto_y\int_w\text{exp}\left[-\frac{\beta}{2}(y-w^Tx)^2\right]\text{exp}\left[-\frac{1}{2}(w-\mu_n)^T\Sigma_n^{-1}(w-\mu_n)\right]dw\\ &\propto_y\int_w\text{exp}\left[-\frac{\beta}{2}(y^2-2w^Txy+w^Txx^Tw)-\frac{1}{2}(w^T\Sigma_n^{-1}w-2w^T\Sigma_n^{-1}\mu_n+\mu_n^T\Sigma_n^{-1}\mu_n)\right]dw \end{align*}$$

其中指数项：

$$\begin{align*} &\beta(y^2-2w^Txy+w^Txx^Tw)+w^T\Sigma_n^{-1}w-2w^T\Sigma_n^{-1}\mu_n\\ =&w^T(\beta xx^T+\Sigma_n^{-1})w-2w^T(\beta xy+\Sigma_n^{-1}\mu_n)+\beta y^2\\ =&(w-m)^TL(w-m)+\dots+\beta y^2\Leftarrow 目标\\ =&w^TLw-2w^TLm+m^TLm+\beta y^2-m^TLm \end{align*}$$

其中：

$$\left\{\begin{aligned}&L=\beta xx^T+\Sigma_n^{-1}\\&Lm=\beta xy+\Sigma_n^{-1}\mu_n\end{aligned}\right.\Rightarrow\left\{\begin{aligned}&L=\beta xx^T+\Sigma_n^{-1}\\&m=L^{-1}(\beta xy+\Sigma_n^{-1}\mu_n)\end{aligned}\right.$$

综上：

$$\begin{align*} p(y|x,X)&\propto_y\int_w\text{exp}\left[-\frac{1}{2}(w-m)^TL(w-m)\right]\text{exp}\left[\frac{1}{2}m^TLm-\frac{1}{2}\beta y^2\right]dw\\ &\propto_y\text{exp}\left[\frac{1}{2}m^TLm-\frac{1}{2}\beta y^2\right]\\ &=\text{exp}\left[-\frac{1}{2}\left(\beta y^2-(\beta xy+\Sigma_n^{-1}\mu_n)^TL^{-1}(\beta xy+\Sigma_n^{-1}\mu_n\right)\right]\\ &=\text{exp}\left[-\frac{1}{2}\left((\beta-\beta^2x^TL^{-1}x)y^2-2\beta yx^TL^{-1}\Sigma_n^{-1}\mu_n-(\Sigma_n^{-1}\mu_n)^TL^{-1}(\Sigma_n^{-1}\mu_n)\right)\right]\\ &\propto \text{exp}\left[-\frac{1}{2}\left(\lambda(y-u)^2\right)\right]\\ &=\text{exp}\left[-\frac{1}{2}(\lambda y^2-2\lambda uy+\lambda u^2)\right]\Leftarrow 我们的目标 \end{align*}$$

其中：

$$\left\{\begin{align}&L=\beta xx^T+\Sigma_n^{-1}\\&\lambda = \beta(1-\beta x^TL^{-1}x)\\&u=\frac{1}{\lambda} \beta x^T L^{-1} \Sigma_n^{-1}\mu_n\end{align}\right.$$

最后得到结果为：

$$p(y|x,X)=N\left(y\left|u,\frac{1}{\lambda}\right.\right)$$

根据S-M公式（舍曼和莫里森于1949年提出）有：

$$\begin{align*} &(A+uv^T)^{-1}=A^{-1}-\frac{A^{-1}uv^TA^{-1}}{1+v^TA^{-1}u}\\ \Rightarrow&L^{-1}=\left(\Sigma_n^{-1}+\beta xx^T\right)^{-1}=\Sigma_n-\frac{\beta \Sigma_nxx^T\Sigma_n}{1+\beta x^T\Sigma_n x} \end{align*}$$

代入 $L$，有：

$$\begin{align*} &x^TL^{-1}x=x^T\Sigma_n x-\frac{\beta x^T\Sigma_nxx^T\Sigma_nx}{1+\beta x^T\Sigma_nx}=\frac{x^T\Sigma_nx}{1+\beta x^T\Sigma_nx}\\ \Rightarrow&\lambda=\beta(1-\beta x^TL^{-1}x)=\frac{\beta}{1+\beta x^T\Sigma_nx}\Rightarrow\text{var}(y)=\frac{1}{\lambda}=\beta^{-1}+x^T\Sigma_nx\\ \Rightarrow&u=\frac{1}{\lambda}\beta x^TL^{-1}\Sigma_{n}^{-1}\mu_n=(\beta^{-1}+x^T\Sigma_nx)\beta x^T\left(\Sigma_n-\frac{\beta \Sigma_nxx^T\Sigma_n}{1+\beta x^T\Sigma_nx}\right)\Sigma_n^{-1}\mu_n=x^T\mu_n \end{align*}$$

估计模型的参数：分类

问题描述：

$$X=\{\langle x_i,y_i\rangle\}_{i=1}^n,\ \sigma_i=p(C_1|x_i)=sigmoid(w^Tx_i)$$

似然函数：

$$p(X|w)=\prod_{i=1}^n\sigma_i^{y_i}(1-\sigma_i)^{1-y_i}$$

先验信息：

$$p(w)=N(\textbf{m}_0,\textbf{S}_0)$$

参数的后验概率函数为：

$$\begin{align*} &p(w|X)\propto p(w)p(X|w)\\ \propto &\text{exp}\left[-\frac{1}{2}(w-\textbf m_0)^T\textbf S_0^{-1}(w-\textbf{m}_0)\right]\times \left[\prod_{i=1}^n\sigma_i^{y_i}(1-\sigma_i)^{1-y_i}\right] \end{align*}$$

解决方案：拉普拉斯近似，当要对某个分布进行近似的时候，可以考虑在它的一阶导数为0的点（mode）进行高斯函数的拟合。

$$\begin{align*} &\text{log}f(x)=\text{log}f(x_0)-\frac{1}{2}a(x-x_0)^2+\dots,\ a=-\frac{d}{dx^2}\text{log}f(x)|x=x_0\\ \Rightarrow &f(x)=f(x_0)\text{exp}\left[-\frac{1}{2}a(x-x_0)^2\right]\\ \Rightarrow &f(x)\propto \text{exp}\left[-\frac{1}{2}a(x-x_0)^2\right] \end{align*}$$

上式为泰勒展开式，其中一阶导项因为 $f’(x_0)=0$ ，为0。因为：

$$q(x)=\sqrt{\frac{a}{2\pi}}\text{exp}\left[-\frac{1}{2}a(x-x_0)^2\right]\sim N(x_0,\frac{1}{a})$$

多维情况为：

$$q(x)=\frac{|\textbf{A}|^{1/2}}{(2\pi)^{d/2}}\text{exp}\left[-\frac{1}{2}(\textbf{x}-\textbf x_0)^T\textbf A(\textbf x-\textbf x_0)\right]\sim N_d(\textbf x_0, \textbf A^{-1})$$

其中，$\textbf A=-\nabla \nabla\text{log}f(x)|_{x=x_0}$。

则：

$$p(w|X)\propto \text{exp}\left[-\frac{1}{2}(w-\textbf m_0)^T\textbf S_0^{-1}(w-\textbf m_0)\right]\times \left[\prod_{i=1}^n\sigma_i^{y_i}(1-\sigma_i)^{1-y_i}\right]$$

其中：

$$\begin{align*} &\textbf A=-\nabla\nabla\text{log}p(w|X)=\textbf S_0^{-1}+\sum_i \sigma_i(1-\sigma_i)x_ix_i^T=\textbf S_0^{-1}+\textbf X^T\pmb {\Lambda}\textbf X\\ \Rightarrow&q(w)\sim N(w_{MAP},\textbf S_n^{-1}), \ \textbf S_n=\textbf A. \end{align*}$$

估计模型的参数：朴素贝叶斯的贝叶斯估计

问题的形式化描述为：

$$D=\{\langle x^{(1)},y_1\rangle,\dots,\langle x^{(n)},y_n\rangle\},\ x^{(i)}=(x_1^{(i)},\dots,x_d^{(i)})^T,\ x_j^{(i)}\in A_j,\ y_i\in\{1,2,\dots,m\}\\ \langle x^{(1)},y_1\rangle,\dots,\langle x^{(n)},y_n\rangle\sim p(x,y|\theta)$$

贝叶斯规则：

$$p(x,y|\theta)=p(y|\theta)p(x|y,\theta)=p(y|\theta)\prod_{j=1}^dp(x_j|y,\theta)$$

其中：

$$\begin{align*} &p(x_j|y,\theta)=P(X_j=x_j|Y=y,\theta)=r_{jy}(x_j)\\ &p(y|\theta)=\pi(y),\ \pi=(\pi(1),\dots,\pi(m))\\ &\sum_{k\in A_j}r_{jy}(k)=1,\ \sum_y \pi(y)=1,\ j=1,2,\dots,d,\ y=1,2,\dots,m,\ x_j\in A_j \end{align*}$$

概率图描述为：

样本似然函数：

$$\begin{align*} p(X|\theta)&=\prod_{i=1}^np(x^{(i)},y_i|\theta)\\ &=\prod_{i=1}^np(y_i|\theta)p(x^{(i)}|y_i,\theta)\\ &=\prod_{i=1}^np(y_i|\theta)\prod_{j=1}^dp(x_j^{(i)}|y_i,\theta)\\ &=\prod_{i=1}^np(y_i|\theta)\prod_{j=1}^d\prod_{k=1}^mp(x_j^{(i)}|y_i,\theta)^{I\{y_i=k\}}\\ &=\prod_{i=1}^np(y_i|\theta)\prod_{j=1}^d\prod_{k=1}^mr_{jk}(x_j^{(i)})^{I\{y_i=k\}}\\ &=\prod_{i=1}^nCat(y_i|\pi)\prod_{j=1}^d\prod_{k=1}^mCat(x_j^{(i)}|r_{jk})^{I\{y_i=k\}} \end{align*}$$

同前面一样，对于似然函数是多项式分布的，我们取它的共轭分布——Dirichlet分布作为先验分布：

$$p(\theta)=Dir(\pi|\alpha)\prod_{jk}Dir(r_{jk}|\beta)$$

则后验函数为：

$$\begin{align*} p(\theta|X)&\propto p(\theta)p(X|\theta)\\ &=\left[Dir(\pi|\alpha)\prod_{jk}Dir(r_{jk}|\beta)\right]\times\left[\prod_{i=1}^nCat(y_i|\pi)\prod_{j=1}^d\prod_{k=1}^mCat(x_j^{(i)}|r_{jk})^{I\{y_i=k\})}\right]\\ &=\left[Dir(\pi|\alpha)\prod_{i=1}^nCat(y_i|\pi)\right]\times\left[\prod_{i=1}^n\prod_{j=1}^d\prod_{k=1}^mDir(r_{jk}|\beta)Cat(x_j^{(i)}|r_{jk})^{I\{y_i=k\}}\right]\\ &=Dir(\pi|\alpha+c)\prod_{jk}Dir(r_{jk}|\beta+d_{jk}) \end{align*}$$

其中，$c_k=\sharp\{i:y_i=k\}, d_{jk}(l)=\sharp\{i:x_j^{(i)}=l,y_i=k\}$。

预测分布：

$$\begin{align*} &p(x,y|\theta)=\pi(y)\prod_{j=1}^d\prod_{k=1}^mr_{jk}(x_j)^{I\{y=k\}}\\ &p(\theta|X)=Dir(\pi|\alpha+c)\prod_{jk}Dir(r_{jk}|\beta+d_{jk})\\ \end{align*}$$

所以：

$$\begin{align*} p(y|x,X)&\propto\int p(x,y|\theta)p(\theta|X)d\theta\\ &\propto \int\left[\pi(y)\prod_{jk}r_{jk}(x_j)^{I\{y=k\}}\right]\times\left[Dir(\pi|\alpha+c)\prod_{jk}Dir(r_{jk}|\beta+d_{jk})\right]d\theta\\ &=\int[\pi(y)Dir(\pi|\alpha+c)]d\pi\times\prod_{jk}\left[\int \left(r_{jk}(x_j)^{I\{y=k\}}Dir(r_{jk}|\beta+d_{jk})\right)dr_{jk}\right]\\ &=\left(\frac{\alpha_y+c_y}{\sum_{y=1}^m(\alpha_y+c_y)}\right)\prod_{j=1}^d\left(\frac{\beta_{x_j}+d_{jy}(x_j)}{\sum_{l=1}^N(\beta_l+d_{jy}(l))}\right)\\ &=\left(\frac{\alpha_y+c_y}{\alpha_0+n}\right)\prod_{j=1}^d\frac{\beta_{x_j}+d_{jy}(x_j)}{\beta_0+c_y} \end{align*}$$

其中，$\sum_{y=1}^m\alpha_y=\alpha_0, \sum_{y=1}^m=n, \sum_{l=1}^N\beta_1=\beta_0, \sum_{l=1}^Nd_{jy}(l)=c_y$。上式从积分推导到期望，其中的 $\pi(y)$和 $r_{jk}(x_j)^{I\{y=k\}}$ 都为概率，各自之和为1，所以积分后化去。

总结

以上的各种方法都是假设已知信息有效的情况下，实验数据偏少，充分利用先验信息，先给定一个初步预测，再通过实验的进行（已知信息的增加）不断对先验信息进行校正，使之不断地向精确偏移，很像某种程度上的拉普拉斯近似。在之前的推导中，可以发现随着实验次数的不断增大，贝叶斯学习的结果将无限趋近于MLE的结果。所以事实上当实验数据量足够大时，先验信息的取值误差对结果的影响很小。即当数据量有保障时，先验信息的取值不用过于谨慎。

EM

引入

当模型中的所有变量均是观测变量时，极大似然估计MLE可以很好地解决问题。但实际中往往会遇到未知地变量存在，称为“隐变量”。这时候，就需要EM（Expectation-Maximization）算法，直译为“期望最大化”算法，其实是一种迭代算法，本质上也是极大似然估计MLE。

极大似然估计为：

$$\theta = \text{arg max }L(\theta)=\text{arg max }\sum_{i=1}^m \text{log}p_\theta(x_i)$$

如果再除以 $m$，则变为期望：

$$\theta = \text{arg max }\left(\sum_{i=1}^m\text{log}p_\theta(x_i)\right)/m=\text{arg max E[log}p_\theta(x_i)]$$

混合高斯模型与EM

混合高斯模型，就是多个高斯模型混合在一起。假设有下图数据，如果只用一个高斯模型去拟合这堆数据，则拟合出的黑色曲线不能很好地解决问题。此时用两个高斯模型的话，对原始数据便能拟合得比较好。

设有随机变量 $\bf x$，则混合高斯模型可以表示为：

$$p(\textbf x)=\sum_{k=1}^K\pi_k\mathcal N(\textbf x|\pmb \mu_k,\pmb \Sigma_k)$$

其中，$\mathcal{N}(\textbf x|\pmb \mu_k,\Sigma_k)$ 称为混合模型中的第 $k$ 个分量，即第 $k$ 个高斯模型，$\pi_k$ 是每个高斯模型的系数，也称混合系数且满足：

$$\sum_{k=1}^K\pi_k=1$$

应用极大似然估计：

$$\text{ln}p(\textbf X|\pmb \pi,\pmb \mu, \pmb \Sigma)=\sum_{n=1}^N\text {ln}\left\{\sum_{k=1}^K\pi_k\mathcal{N}(\textbf {x}_n|\pmb \mu_k,\pmb \Sigma_k)\right\}$$

其中 $N$ 表示样本总数，需要求解的参数 $\{\pmb \pi_k,\pmb \mu_k,\pmb \Sigma_k\}$。对上式各个变量求偏导并令其等于0：

1. 对 $\mu_k$ 求导并令其等于0得到：

   $$0=-\sum_{n=1}^N\frac{\pi_k\mathcal{N}(\textbf{x}_n|\pmb \mu_k,\Sigma_k)}{\sum_j\pi_j\mathcal{N}(\textbf x_n|\pmb\mu_j,\pmb\Sigma_j)}\Sigma_k^{-1}(\textbf{x}_n-\pmb\mu_k)$$

   其中左项比较复杂，用 $\gamma(z_{nk})$ 表示。化简之后得到：

   $$\pmb\mu_k=\frac{1}{N_k}\sum_{n=1}^N\gamma(z_{nk})\textbf x_n\\ N_k=\sum_{n=1}^N\gamma(z_{nk})$$

2. 对 $\pmb\pi_k$ 和 $\pmb\Sigma_k$ 分别求偏导并令其等于0得到：

   $$\pmb\pi_k=\frac{N_k}{N}\\ \pmb\Sigma_k=\frac{1}{N_k}\sum_{n=1}^N\gamma(z_{nk})(\textbf x_n-\pmb\mu_k)(\textbf x_n-\pmb\mu_k)^T$$

   其中 $\pmb \pi_k$ 的求取，需要加入限制条件，$\sum_{k=1}^K\pmb\pi_k=1$，因此需要加入拉格朗日算子：

   $$\text{ln}p(\textbf x|\pmb\pi,\pmb\mu,\pmb\Sigma)+\lambda(\sum_{k=1}^K\pmb\pi_k-1)$$

   求上式关于 $\pmb\pi_k$ 的最大似然函数，得到：

   $$0=\sum_{n=1}^N\frac{\mathcal{N}(\textbf x_n|\pmb\mu_k,\pmb\Sigma_k)}{\sum_j\pmb\pi_j\mathcal{N}(\textbf x_n|\pmb\mu_j,\pmb\Sigma_j)}+\lambda$$

   上式两边同乘 $\pmb\pi_k$，可以得到 $\lambda=-N$，进而可以得到 $\pmb\pi_k$ 更简洁的表达式：

   $$\pmb\pi_k=\frac{N_k}{N}$$

这时候可以看到，要求的参数 $\{\pmb\pi_k,\pmb\mu_k,\pmb\Sigma_k\}$ 的表达式都包含了 $\gamma(z_{nk})$，而这一项的求解又依赖于 $\{\pmb\pi_k,\pmb\mu_k,\pmb\Sigma_k\}$ 的值，似乎陷入了循环。

所以可以先给 $\{\pmb\pi_k,\pmb\mu_k,\pmb\Sigma_k\}$ 一些初始值，之后计算 $\gamma(z_{nk})$ ，之后再计算 $\{\pmb\pi_k,\pmb\mu_k,\pmb\Sigma_k\}$ 。这就是EM算法求解GMM的流程了，把整个流程以算法的形式表示：

• E-step

  $$\gamma(z_{nk})=\frac{\pi_k\mathcal{N}(\textbf{x}_n|\pmb \mu_k,\Sigma_k)}{\sum_j\pi_j\mathcal{N}(\textbf x_n|\pmb\mu_j,\pmb\Sigma_j)}$$

• M-step

  $$\begin{align*} &\pmb\mu_k^{\text{new}}=\frac{1}{N_k}\sum_{n=1}^N\gamma(z_{nk})\textbf x_n\\ &\pmb\Sigma_k^{\text{new}}=\frac{1}{N_k}\sum_{n=1}^N\gamma(z_{nk})(\textbf x_n-\pmb\mu_k^{\text{new}})(\textbf x_n-\pmb\mu_k^{\text{new}})^T\\ &\pmb\pi_k^{\text{new}}=\frac{N_k}{N} \end{align*}$$

效果图如下：

GMM的贝叶斯理解

GMM模型表示为：

$$p(\textbf{x})=\sum_{k=1}^K\pmb\pi_k\mathcal N(\textbf x|\pmb\mu_k,\pmb\Sigma_k)$$

其中，$\pmb\pi_k$ 表示第 $k$ 个高斯模型的系数（权重），或者第 $k$ 类被选中的概率。我们引入一个 $k$ 维的随机变量 $z$。$z_k$ 表示它的第 $k$ 个维度，只能取1或0两个值。$z_k=1$ 表示第 $k$ 类被选中，$p(z_k=1)=\pmb\pi_k$。用概率图表示如下：

隐变量 $z$ 是独立同分布的，可以写出它的联合概率分布：

$$p(\textbf z)=p(z_1)p(z_2)\dots p(z_K)=\prod_{k=1}^K\pi_k^{z_k}$$

设样本被分为两类，则可以用条件概率来表示：

$$p(\textbf x|z_k=1)=\mathcal{N}(\textbf x|\pmb\mu_k,\pmb\Sigma_k)$$

综合以上两式，可以写成：

$$p(\textbf x|\textbf z)=\prod_{k=1}^K\mathcal{N}(\textbf x|\pmb\mu_k,\pmb\Sigma_k)^{z_k}$$

根据条件概率公式，求出 $p(\textbf x)$ 的形式：

$$p(\textbf x)=\sum_{\textbf z}p(\textbf z)p(\textbf x|\textbf z)=\sum_{k=1}^K\pmb\pi_k\mathcal N(\textbf x|\pmb\mu_k,\pmb\Sigma_k)$$

隐变量中的“隐”的意义是：我们知道数据可以分成两类，但是随机抽取一个数据点，不知道这个数据点属于第一类还是第二类，它的分类我们观察不到，因此引入一个隐含变量 $z$ 来描述这个现象。

注意到在贝叶斯的思想下， $p(\textbf z)$ 是先验概率，$p(\textbf x|\textbf z)$ 是似然概率，很自然我们会想到求出后验概率：

$$\begin{align*} \gamma(z_k)&=p(z_k=1|\textbf x)\\ &=\frac{p(z_k=1)p(\textbf{x}|z_k=1)}{p(\textbf x|\textbf z)}\\ &=\frac{p(z_k=1)p(\textbf x|z_k=1)}{\sum_{j=1}^Kp(z_j=1)p(\textbf x|z_j=1)}(全概率公式)\\ &=\frac{\pmb\pi_k\mathcal{N}(\textbf{x}|\pmb\mu_k,\pmb\Sigma_k)}{\sum_{j=1}^K\pmb\pi_j\mathcal{N}(\textbf x|\pmb\mu_j,\pmb\Sigma_j)} \end{align*}$$

上式的 $\gamma(z_k)$ 是对于一个样本而言。

完备数据

现在假设有 $N$ 个样本，$K$ 个类，则完备数据（有观测数据，即类别标签）的似然表示为：

$$p(\textbf X,\textbf Z|\pmb\mu, \pmb\Sigma, \pmb\pi)=\prod_{n=1}^N\prod_{k=1}^K\pi_k^{z_{nk}}\mathcal{N}(\textbf x_n|\pmb\mu_k,\pmb\Sigma_k)^{z_{nk}}$$

对数似然为：

$$\text{ln}p(\textbf X,\textbf Z|\pmb\mu, \pmb\Sigma, \pmb\pi)=\sum_{n=1}^N\sum_{k=1}^Kz_{nk}\{\text{ln}\pi_k+\text{ln}\mathcal{N}(\textbf x_n|\pmb\mu_k, \pmb\Sigma_k)\}$$

$z_{nk}$ 的期望值为：

$$\mathbb E[z_{nk}]=\frac{\pi_k\mathcal{N}(\textbf{x}_n|\pmb \mu_k,\Sigma_k)}{\sum_j\pi_j\mathcal{N}(\textbf x_n|\pmb\mu_j,\pmb\Sigma_j)}=\gamma(z_{nk})$$

则完备数据下对数似然的期望值为：

$$\mathbb{E}_{\textbf Z}[\text{ln }p(\textbf X,\textbf Z|\pmb\mu, \pmb\Sigma, \pmb\pi)]=\sum_{n=1}^N\sum_{k=1}^K\gamma(z_{nk})\{\text{ln} \pi_k+\text{ln}\mathcal{N}(\textbf x_n|\pmb\mu_k,\pmb\Sigma_k)\}$$

完备数据情况下，最大似然估计即可。

非完备数据

• 不引入隐变量
  非完备数据的对数似然表示为：

  $$\text{ln}p(\textbf X|\pmb\mu, \pmb\Sigma, \pmb\pi)=\sum_{n=1}^N\text{ln}\left\{\sum_{k=1}^K\pi_k\mathcal{N}(\textbf x_n|\pmb\mu_k,\pmb\Sigma_k)\right\}$$

  ​ 之后求偏导。

• 引入隐变量
  先给定隐变量 $\textbf Z$，之后按完备数据的最大似然求取。

以上两种方法都能推导出EM。

最大后验估计与EM

1. 初始化设置参数 $\pmb\theta^{\text{old}}$；

2. E-step计算 $p(\textbf Z|\textbf X,\pmb\theta^{\text{old}})$；

3. M-step计算

   $$\pmb\theta^{\text{new}}=\text{arg max}_{\pmb\theta}\mathcal{Q}(\pmb\theta,\pmb\theta^{\text{old}})$$

   其中，

   $$\mathcal{Q}(\pmb\theta,\pmb\theta^{\text{old}})=\sum_{\textbf Z}p(\textbf Z|\textbf X,\pmb\theta^{\text{old}})\text{ln}p(\textbf X,\textbf Z|\pmb\theta)+\text{ln}p(\pmb\theta)$$

k-means与EM

下图为生成的样本示意图，可以发现样本属于哪个分布不是hard的0-1，而是soft的[0, 1]概率。

可以发现上图其实与k-means有些相似，区别在于k-means是hard的0-1。所以说，k-means其实属于hard的EM算法。其中：

$$J=\sum_{n=1}^N\sum_{k=1}^Kr_{nk}||\textbf x_n-\pmb\mu_k||^2$$

E-step：

$$r_{nk}=\left\{\begin{aligned}&1&&\text{if }k=\text{arg min}_j||\textbf x_n-\pmb\mu_j||^2\\&0&&\text{otherwise.}\end{aligned}\right.$$

M-step：

$$\pmb\mu_k=\frac{\sum_nr_{nk}\textbf{x}_n}{\sum_nr_{nk}}$$

也可以将k-means从hard改为上述soft的混合高斯模型，即 $r_{nk}\in\{0,1\}$ 改为 $r_{nk}\in(0, 1)$，

$$p(\textbf x|\pmb\mu_k,\pmb\Sigma_k)=\frac{1}{(2\pi\epsilon)^{1/2}}\exp\left\{-\frac{1}{2\epsilon}||\textbf x-\pmb\mu_k||^2\right\}$$

隶属度为：

$$\gamma_{nk}=\frac{\pi_k\text{exp}\{-||\textbf x_n-\pmb\mu_k||^2/2\epsilon\}}{\sum_j\pi_j\text{exp}\{-||\textbf{x}_n-\pmb\mu_j||^2/2\epsilon\}}$$

期望为：

$$\mathbb{E}_\textbf Z[\text{ln}p(\textbf X,\textbf Z|\pmb\mu,\pmb\Sigma,\pmb\pi)]\rightarrow-\frac{1}{2}\sum_{n=1}^N\sum_{k=1}^Kr_{nk}||\textbf x_n-\pmb\mu_k||^2+\text{const}$$

算法推导（收敛性）

EM算法可以表示为下式：

其中 $\log p(X,Z|\theta)$ 可以看成一个随机变量的函数，而 $p(Z|X,\theta^{(t)})$ 是一个分布，这样积分后，就是期望，也就是E-step所做的；而M-step就是最大化该期望。

EM算法推导如下：

首先对于概率分布为 $p_\theta(x_i)$ 的 $n$ 次观测，我们将它添上隐变量 $z$，

$$p_\theta(x_i)\Rightarrow p_\theta(x_i,z)$$

之后取对数似然函数，

$$L(\theta)=\sum_{i=1}^n\ln p_\theta(x_i)$$

为了将隐变量 $z$ 强制暴露出来，将 $p_\theta(x_i)$ 化为 $p_\theta(x_i,z)$ 接着对 $z$ 求积分，可以得到：

$$L(\theta)=\sum_{i=1}^n\ln \sum_zp_\theta(x_i,z)$$

接下来要做的是求使 $L(\theta)$ 最大的 $\theta$ 值，

$$\theta^\ast=\text{arg max}_\theta\left\{\sum_{i=1}^n\text{ln}\sum_zp_\theta(x_i,z)\right\}$$

但由于此时 $L(\theta)$ 含有隐随机变量 $z$ ，直接找参数估计是困难的。我们的策略是建立 $L(\theta)$ 的下界，并求该下界的最大值，重复这个过程，直到收敛到局部最大值。

由上一步，已知不能直接最大化 $L(\theta)$，采用不断建立 $L(\theta)$ 的下界（E步），然后优化下界（M步）的方法。

具体如下：假定对于每一个 $x_i$ 样本，都有 $m$ 类的隐变量分布 $z_j$ ，并把隐变量 $z_j$ 的概率分布记为 $Q_i(z_j)$。已知，

$$L(\theta)=\sum_{i=1}^n\text{ln}\sum_{j=1}^mp_\theta(x_i,z_j)$$

为了将 $Q_i(z_j)$ 引入，就在求和公式里面乘以 $Q_i(z_j)$，再除以 $Q_i(z_j)$，

$$L(\theta)=\sum_{i=1}^n\text{ln}\sum_{j=1}^mQ_i(z_j)\frac{p_\theta(x_i,z_j)}{Q_i(z_j)}$$

为什么要这么做呢？根据Lazy Statistician规则：

Lazy Statistician规则：

设 $Y$ 是随机变量 $X$ 的函数，$Y=g(X)$（$g$ 是连续函数），$X$ 是离散型随机变量，它的分布为 $p(X=x_k)=p_k$，其中 $k=1,2,\dots$。

若 $\sum_{k=1}^\infty g(x_k)p_k$ 的值绝对收敛，则有，

$$E(Y)=E(g(X))=\sum_{k=1}^\infty g(x_k)p_k$$

所以，在引入 $Q_i(z_j)$ 以后，就得到了所要求的期望：

$$E\left[\frac{p_\theta(x_i,z_j)}{Q_i(z_j)}\right]=\sum_z Q_i(z_j)\frac{p_\theta(x_i,z_j)}{Q_i(z_j)}$$

接着，还需用到Jensen不等式：

Jensen不等式：

若 $f$ 是凹函数， $X$ 是随机变量，那么

$$E(f(X))\ge f(E(X))$$

特别地，如果 $f$ 是严格凹函数，那么当且仅当 $P(x=E(X))=1$ 时取等号，此时 $X$ 为常量，有：

$$E(f(X))=f(E(X))$$

另：Jensen不等式应用于凸函数时，不等号方向反向，

$$E(f(X))\le f(E(X))$$

根据Jensen不等式，分析如下，

$$\begin{align*} \text{ln}\frac{p_\theta(x_i,z_j)}{Q_i(z_j)}&\Rightarrow f(X)\\ \sum_{j=1}^mQ_i(z_j)\text{ln}\frac{p_\theta(x_i,z_j)}{Q_i(z_j)}&\Rightarrow E[f(X)]\\ \sum_{j=1}^mQ_i(z_j)\frac{p_\theta(x_i,z_j)}{Q_i(z_j)}=E\left[\frac{p_\theta(x_i,z_j)}{Q_i(z_j)}\right]&\Rightarrow E[X]\\ \text{ln}\sum_{j=1}^mQ_i(z_j)\frac{p_\theta(x_i,z_j)}{Q_i(z_j)}=\ln E\left[\frac{p_\theta(x_i,z_j)}{Q_i(z_j)}\right]&\Rightarrow f(E[X]) \end{align*}$$

由于对数函数 $\text{ln}$ 是凸函数，根据Jensen不等式，

$$\text{ln}\sum_zQ_i(z_j)\frac{p_\theta(x_i,z_j)}{Q_i(z_j)}\ge\sum_zQ_i(z_j)\ln \frac{p_\theta(x_i,z_j)}{Q_i(z_j)}$$

添上求和符号：

$$\sum_{i=1}^n\text{ln}\sum_{j=1}^mQ_i(z_j)\frac{p_\theta(x_i,z_j)}{Q_i(z_j)}\ge\sum_{i=1}^n\sum_{j=1}^mQ_i(z_j)\text{ln}\frac{p_\theta(x_i,z_j)}{Q_i(z_j)}$$

所以之前的最大似然为：

$$\theta^\ast=\text{arg max}_\theta\sum_{i=1}^n\sum_{j=1}^mQ_i(z_j)\text{ln}\frac{p(x_i,z_j)}{Q_i(z_j)}=\text{arg max}_\theta Q_i(z_j)\text{ln}p(x_i,z_j)$$

知道最大似然可以通过最大化其下界来求取，但是其中的 $Q_i(z_j)$ 具体是多少呢？

根据Jensen不等式，取等号的条件为：

$$\frac{p_\theta(x_i,z_j)}{Q_i(z_j)}=c$$

其中 $c$ 为常数，即

$$Q_i(z_j)=\frac{p(x_i,z_j)}{c}$$

对于每一个样本 $x_i$，$Q_i(z_j)$ 表示该样本隐含变量 $z$ 的某种分布，根据概率分布的性质可知，

$$\sum_{j=1}^mQ_i(z_j)=1$$

所以，

$$\begin{align*} \sum_{j=1}^m\frac{p(x_i,z_j)}{c}&=1\\ \frac{1}{c}\sum_{j=1}^mp(x_i,z_j)&=1\\ \sum_{j=1}^mp(x_i,z_j)&=c \end{align*}$$

则，

$$Q_i(z_j)=\frac{p(x_i,z_j)}{c}=\frac{p(x_i,z_j)}{\sum_{j=1}^mp(x_i,z_j)}=\frac{p(x_i,z_j)}{p(x_i)}=p(z_j|x_i)$$

又，

$$Q_i(z_j)=p(z_j|x_i)=\frac{p(x_i,z_j)}{\sum_{j=1}^mp(x_i,z_j)}=\frac{p(z_j)p(x_i|z_j)}{\sum_{j=1}^mp(z_j)p(x_i|z_j)}$$

$Q_i(z_j)$ 的值我们已经知道了，代回最大似然推导：

$$\begin{align*} \theta^\ast&=\text{arg max}_\theta\sum_{i=1}^n\sum_{j=1}^mQ_i(z_j)\ln p(x_i,z_j)\\ &=\text{arg max}_\theta\sum_{i=1}^nE_{Q_i(z)}[\ln p(x_i,z)] \end{align*}$$

之后才是这一节的正题，收敛性。因为，

$$p(x)=\frac{p(x,z)}{p(z|x)}$$

所以，

$$\ln p(x|\theta)=\ln p(x,z|\theta)-\ln p(z|x,\theta)$$

所以，

$$E_{p(z|x,\theta^{(t)})}[\ln p(x|\theta)]=E_{p(z|x,\theta^{(t)})}[\ln p(x,z|\theta)-\ln p(z|x,\theta)]$$

所以，

$$\begin{align*} \ln p(x|\theta)&=E_{p(z|x,\theta^{(t)})}[\ln p(x,z|\theta)-\ln p(z|x,\theta)]\\ &=\int_z[\ln p(x,z|\theta)]p(z|x,\theta^{(t)})dz-\int_z[\ln p(z|x,\theta)]p(z|x,\theta^{(t)})dz\\ &=Q(\theta,\theta^{(t)})-H(\theta,\theta^{(t)}) \end{align*}$$

因为，

$$\ln p(x|\theta^{(t+1)})\ge \ln p(x|\theta^{(t)})$$

所以，

$$Q(\theta^{(t+1)},\theta^{(t)})-H(\theta^{(t+1)},\theta^{(t)})\ge Q(\theta^{(t)},\theta^{(t)})-H(\theta^{(t)},\theta^{(t)})$$

所以我们要确定的是

$$\left\{\begin{align}&Q(\theta^{(t+1)},\theta^{(t)})\ge Q(\theta^{(t)},\theta^{(t)})\\&H(\theta^{(t+1)},\theta^{(t)})\le H(\theta^{(t)},\theta^{(t)})\end{align}\right.$$

由EM的迭代式，

$$\theta^{(t+1)}=\text{arg max}_\theta\int\ln p(X,Z|\theta)p(Z|X,\theta^{(t)})dZ=\text{arg max}_\theta Q(\theta,\theta^{(t)})$$

得，

$$Q(\theta^{(t+1)},\theta^{(t)})\ge Q(\theta,\theta^{(t)})$$

接下来证明 $H(\theta^{(t+1)},\theta^{(t)})\le H(\theta^{(t)},\theta^{(t)})$，

$$\text{If }H(\theta,\theta^{(t)})\le H(\theta^{(t)},\theta^{(t)})\ \forall\theta \text{, Then } H(\theta^{(t+1)},\theta^{(t)})\le H(\theta^{(t)},\theta^{(t)})$$

则，

$$\begin{align*} &H(\theta^{(t)},\theta^{(t)})-H(\theta,\theta^{(t)})\\ =&\int_z[\ln p(z|x,\theta^{(t)})]p(z|x,\theta^{(t)}) dz-\int_z[\ln p(z|x,\theta)]p(z|x,\theta^{(t)})dz\\ =&\int_z\left[\ln \frac{p(z|x,\theta^{(t)})}{p(z|x,\theta)}\right]p(z|x,\theta^{(t)})dz\\ =&\int_z-\left[\ln \left(\frac{p(z|x,\theta)}{p(z|x,\theta^{(t)})}\right)\right]p(z|x,\theta^{(t)})dz\\ \ge&-\ln \int_z\left(\frac{p(z|x,\theta)}{p(z|x,\theta^{(t)})}\right)p(z|x,\theta^{(t)})dz \Leftarrow \text{Jensen不等式}\\ =&-\ln 1\\ =&0 \end{align*}$$

所以，

$$\left\{\begin{align}&Q(\theta^{(t+1)},\theta^{(t)})\ge Q(\theta^{(t)},\theta^{(t)})\\&H(\theta^{(t+1)},\theta^{(t)})\le H(\theta^{(t)},\theta^{(t)})\end{align}\right.\Rightarrow\ln p(x|\theta^{(t+1)})\ge\ln p(x|\theta^{(t)})$$

硬币模型

有两个不均匀硬币A、B，每次实验抛10次，记录正面和反面次数。本例中记为 $H(head)$ 和 $T(tail)$，多次这样的实验后，预测硬币A和B各自的正面概率。

第一种情况是实验员记录了每次丢的是哪个硬币，没有隐变量，直接最大似然就可以算出A和B的正面概率 $\hat\theta_A,\hat\theta_B$。

A硬币丢了24+6=30次，其中24次正面 $\hat\theta_A=0.8$，$\hat\theta_B$ 同理可计算。

第二种情况就是忘记记录每次丢的是A硬币还是B硬币了。我们将观测到的结果记录成以下的表格，并使用EM算法的公式计算。

上表中给出的是不同样本中出现”正面朝上“和“反面朝上”的次数，EM算法在本例中的计算公式为：

• E-step

  $$Q_i(z_j)=\frac{p(z_j)p(x_i|z_j)}{\sum_{j=1}^mp(z_j)p(x_i|z_j)}$$

• M-step

  $$\theta^\ast=\text{arg max}_\theta\sum_{i=1}^n\sum_{j=1}^mQ_i(z_j)\text{ln}p(x_i,z_j)$$

计算过程为：

首先随机给定两个 $\hat\theta_A,\hat\theta_B$，假设为0.6和0.5，硬币个数 $m=2$，使用A硬币和B硬币的概率分别为0.5，进行E-step计算。

• E-step

$$\begin{align*} &Q_i(z_j)=\frac{p(z_j)p(x_i|z_j)}{\sum_{j=1}^mp(z_j)p(x_i|z_j)},\ \theta_A=0.6,\ \theta_B=0.5,\ m=2,\ z_1=A, z_2=B, p(z_1)=p(z_2)=0.5\\ &Q_1(z_1)=\frac{p(z_1)p(x_1|z_1)}{p(z_1)p(x_1|z_1)+p(z_2)p(x_1|z_2)}=\frac{0.5\times0.6^5\times(1-0.6)^5}{0.5\times0.6^5\times(1-0.6)^5+0.5\times0.5^5\times(1-0.5)^5}=0.45\\ &Q_2(z_1)=\frac{p(z_1)p(x_2|z_1)}{p(z_1)p(x_2|z_1)+p(z_2)p(x_2|z_2)}=\frac{0.5\times0.6^9\times(1-0.6)^1}{0.5\times0.6^9\times(1-0.6)^1+0.5\times0.5^9\times(1-0.5)^1}=0.8\\ &Q_3(z_1)=\frac{p(z_1)p(x_3|z_1)}{p(z_1)p(x_3|z_1)+p(z_2)p(x_3|z_2)}=\frac{0.5\times0.6^8\times(1-0.6)^2}{0.5\times0.6^8\times(1-0.6)^2+0.5\times0.5^8\times(1-0.5)^2}=0.73\\ &Q_4(z_1)=\frac{p(z_1)p(x_4|z_1)}{p(z_1)p(x_4|z_1)+p(z_2)p(x_4|z_2)}=\frac{0.5\times0.6^4\times(1-0.6)^6}{0.5\times0.6^4\times(1-0.6)^6+0.5\times0.5^4\times(1-0.5)^6}=0.35\\ &Q_5(z_1)=\frac{p(z_1)p(x_5|z_1)}{p(z_1)p(x_5|z_1)+p(z_2)p(x_5|z_2)}=\frac{0.5\times0.6^7\times(1-0.6)^3}{0.5\times0.6^7\times(1-0.6)^3+0.5\times0.5^7\times(1-0.5)^3}=0.65\\ \end{align*}$$

得：

• M-step

$$\begin{align*} \theta^\ast&=\text{arg max}_\theta\sum_{i=1}^n\sum_{j=1}^mQ_i(z_j)\text{ln}p(x_i,z_j)\\ &=\text{arg max}_\theta\sum_{i=1}^n\left[Q_i(z_1)\text{ln}p(x_i,z_1)+Q_i(z_2)\text{ln}p(x_i,z_2)\right]\\ &=0.45\text{ln}\left[p(z_1)\theta_A^5(1-\theta_A)^5\right]+0.55\text{ln}\left[p(z_2)\theta_B^5(1-\theta_B)^5\right]\\ &+0.8\text{ln}\left[p(z_1)\theta_A^9(1-\theta_A)^1\right]+0.2\text{ln}\left[p(z_2)\theta_B^9(1-\theta_B)^1\right]\\ &+0.73\text{ln}\left[p(z_1)\theta_A^8(1-\theta_A)^2\right]+0.27\text{ln}\left[p(z_2)\theta_B^8(1-\theta_B)^2\right]\\ &+0.35\text{ln}\left[p(z_1)\theta_A^4(1-\theta_A)^6\right]+0.65\text{ln}\left[p(z_2)\theta_B^4(1-\theta_B)^6\right]\\ &+0.65\text{ln}\left[p(z_1)\theta_A^7(1-\theta_A)^3\right]+0.35\text{ln}\left[p(z_2)\theta_B^7(1-\theta_B)^3\right]\\ \end{align*}$$

算完一轮得到下表，例如 $Q_1(z_1)=0.45$ 实际抛掷中有5次正面朝上，所以期望为 $0.45\times 5=2.25$，

coin A
1=2.2H, 2.2T
2=7.2H, 0.8T
3=5.9H, 1.5T
4=1.4H, 2.1T
5=4.5H, 1.9T
all=21.3H, 8.6T

所有的计算出来以后，求 $\hat\theta_A, \hat\theta_B$：0.596是怎么计算的？

$$\left\{\begin{aligned}\hat\theta_A^{(1)}=\frac{21.3}{21.3+8.6}=0.71\\\hat\theta_B^{(1)}=\frac{11.7}{11.7+8.4}=0.58\end{aligned}\right.\Rightarrow p(z_1)=0.596$$

迭代多次后正面向上概率收敛，即为预测结果：

$$\left\{\begin{aligned}\hat\theta_A^{(10)}=0.80\\\hat\theta_B^{(10)}=0.52\end{aligned}\right.$$

隐马尔可夫HMM与EM

在参数学习问题中，依题意，我们写出未知参数的一些性质：

$$\begin{align*} &\theta = (A,B,\Pi),\ A=(a_{ij})_{n\times n},\ B=(b_{jk})_{n\times m},\ \Pi=(\pi_i)_{1\times n}\\ &\sum_{j=1}^na_{ij}=1,\ \sum_{k=1}^mb_{jk}=1,\ \sum_{i=1}^n\pi_i=1 \end{align*}$$

还有几个重要计算公式：

$$\gamma_t(i)=p(i_t=i|O,\lambda)=\frac{P(i_t=i,O|\lambda)}{P(O|\lambda)}=\frac{\alpha_t(i)\beta_t(i)}{\sum_{j=1}^N\left[\alpha_t(j)\beta_t(j)\right]}$$ $$\begin{align*} \xi_t(i,j)&=P(i_t=i,i_{t+1}=j|O,\lambda)\\ &=\frac{P(i_t=i,i_{t+1}=j,O|\lambda)}{P(O|\lambda)}\\ &=\frac{\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)}{\sum_{i=1}^N\sum_{j=1}^N\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)} \end{align*}$$ $$\gamma_t(i)=P(i_t=i|O,\lambda)=\sum_jP(i_t=i,i_{t+1}=j|O,\lambda)=\sum_j\xi_t(i,j)$$

E-step：

已知：

$$P(X,Z)=P(Z_1)\prod_{t=1}^{T-1}P(Z_{t+1}|Z_t)\prod_{t=1}^TP(X_t|Z_t)$$

其中，$Z_t\in\{1,2,\dots,n\}$，$X_t\in\{1,2,\dots,m\}$。求状态序列的分布为：

$$Q(Z)=P(Z|X,\theta^{(t)})=\frac{P(X,Z)}{\sum_ZP(X,Z)}$$

M-step：

最大似然如下：

$$\theta^{(t+1)}=\text{arg max}_\theta L$$

其中，

$$\begin{align*} L&=\sum_ZQ(Z|\theta^{(t)})\text{ln}P(X,Z|\theta)\\ &=\sum_ZQ(Z|\theta^{(t)})\left(\text{ln}P(Z_1)+\sum_{t=1}^{T-1}\text{ln}P(Z_{t+1}|Z_t)+\sum_{t=1}^T\text{ln}P(X_t|Z_t)\right) \end{align*}$$

$L$ 分别对 $\Pi, A, B$ 分别求偏导，并令偏导等于零，则可求出 $\theta^{(t+1)}$。

• $\Pi$

$$\begin{align*} L(\pi,\lambda)&=\sum_{Z_{1,\dots,T}}Q(Z|\theta^{(t)})\text{ln}P(Z_1)+\lambda\left(1-\sum_{i=1}^n\pi_i\right)\\ &=\sum_{Z_1}Q(Z_1|\theta^{(t)})\text{ln}P(Z_1)+\lambda\left(1-\sum_{i=1}^n\pi_i\right)\\ &=\sum_{i=1}^nQ(Z_1=i|\theta^{(t)})\text{ln}P(Z_1=i)+\lambda\left(1-\sum_{i=1}^n\pi_i\right)\\ &=\sum_{i=1}^n\gamma_1(i)\text{ln}\pi_i+\lambda\left(1-\sum_{i=1}^n\pi_i\right) \end{align*}$$

求偏导：

$$\frac{\partial}{\partial \pi_i}L(\pi, \lambda)=\frac{\gamma_1(i)}{\pi_i}-\lambda=0\Rightarrow\pi_i=\frac{\gamma_1(i)}{\lambda}$$

因为，

$$\sum_{i=1}^n\pi_i=1\Rightarrow\lambda=\sum_{i=1}^n\gamma_1(i)\Rightarrow\pi_i=\frac{\gamma_1(i)}{\sum_{i=1}^n\gamma_1(i)},\ i=1,2,\dots,n$$

• $A$

$$\begin{align*} L(A,\lambda)&=\sum_{Z_{1,\dots,T}}Q(Z|\theta^{(t)})\sum_{t=1}^{T-1}\text{ln}P(Z_{t+1}|Z_t)+\lambda\left(1-\sum_{j=1}^na_{ij}\right)\\ &=\sum_{t=1}^{T-1}\sum_{i=1}^n\sum_{j=1}^nQ(Z_t=i,Z_{t+1}=j|\theta^{(t)})\text{ln}P(Z_{t+1}=j|Z_t=j)+\lambda\left(1-\sum_{j=1}^na_{ij}\right)\\ &=\sum_{t=1}^{T-1}\sum_{i=1}^n\sum_{j=1}^n\xi_t(i,j)\text{ln}a_{ij}+\lambda\left(1-\sum_{j=1}^na_{ij}\right) \end{align*}$$

求偏导：

$$\frac{\partial }{\partial a_{ij}}L(A,\lambda)=\frac{\sum_{t=1}^{T-1}\xi_t(i,j)}{a_{ij}}-\lambda=0\Rightarrow a_{ij}=\frac{\sum_{t=1}^{T-1}\xi_t(i,j)}{\lambda}$$

因为，

$$\sum_{j=1}^na_{ij}=1\Rightarrow\lambda=\sum_{t=1}^{T-1}\gamma_t(i)\Rightarrow a_{ij}=\frac{\sum_{t=1}^{T-1}\xi_t(i,j)}{\sum_{t=1}^{T-1}\gamma_t(i)},\ i,j=1,2,\dots,n$$

• $B$

$$\begin{align*} L(B,\lambda)&=\sum_{Z_{1,\dots,T}}Q(Z|\theta^{(y)})\sum_{t=1}^T\text{ln}P(X_t|Z_t)+\lambda\left(1-\sum_{k=1}^mb_{jk}\right)\\ &=\sum_{t=1}^T\sum_{j=1}^n\sum_{k=1}^mQ(Z_t=j,X_t=k|\theta^{(t)})\text{ln}P(X_t|Z_t=j)+\lambda\left(1-\sum_{k=1}^mb_{jk}\right)\\ &=\sum_{t=1}^T\sum_{j=1}^n\sum_{k=1}^mQ(Z_t=j,X_t=k|\theta^{(t)})\text{ln}b_{jk}+\lambda\left(1-\sum_{k=1}^mb_{jk}\right) \end{align*}$$

求偏导：

$$\frac{\partial}{\partial b_{jk}}L(B,\lambda)=\frac{\sum_{t=1}^TQ(Z_t=j,X_t=k|\theta^{(t)})}{b_{jk}}-\lambda=0\Rightarrow b_{jk}=\frac{\sum_{t=1}^TQ(Z_t=j,X_t=k|\theta^t)}{\lambda}$$

因为，

$$\sum_{k=1}^mb_{jk}=1\Rightarrow\lambda=\sum_{t=1}^TQ(Z_t=j|\theta^{(t)})=\sum_{t=1}^T\gamma_t(j)\Rightarrow b_{jk}=\frac{\sum_{t=1,X_t=k}^T\gamma_t(j)}{\sum_{t=1}^{T-1}\gamma_t(j)}$$

混合伯努利模型与EM

列出伯努利模型概率分布以及均值和协方差：

• 单个模型：

  $$\begin{align*} p(\textbf x|\pmb\mu)&=\prod_{i=1}^D\mu_i^{x_i}(1-\mu_i)^{(1-x_i)}\\ \mathbb{E}[\textbf x]&=\pmb\mu\\ \text{cov}[\textbf x]&=\text{diag}\{\mu_i(1-\mu_i)\} \end{align*}$$

• 混合模型：

  $$\begin{align*} p(\textbf x|\pmb\mu,\pmb\pi)&=\sum_{k=1}^K\pi_kp(\textbf x|\pmb\mu_k)\\ p(\textbf x|\pmb\mu_k)&=\prod_{i=1}^D\mu_{ki}^{x_i}(1-\mu_{ki})^{(1-x_i)}\\ \mathbb E[\textbf x]&=\sum_{k=1}^K\pi_k\pmb\mu_k\\ \text{cov}[\textbf x]&=\sum_{k=1}^K\pi_k\{\pmb\Sigma_k+\pmb\mu_k\pmb\mu_k^T\}-\mathbb E[\textbf x]E[\textbf x]^T \end{align*}$$

则对数似然函数为：

$$\text{ln}p(\textbf X|\pmb\mu,\pmb\pi)=\sum_{n=1}^N\text{ln}\left\{\sum_{k=1}^K\pi_kp(\textbf x_n|\pmb\mu_k)\right\}$$

引入隐变量：

$$p(\textbf x|\textbf z,\pmb\mu)=\prod_{k=1}^Kp(\textbf x|\pmb\mu_k)^{z_k}\\ p(\textbf z|\pmb\pi)=\prod_{k=1}^K\pi_k^{z_k}$$

完备数据下的对数似然为：

$$\text{ln}p(\textbf X, \textbf Z|\pmb\mu,\pmb\pi)=\sum_{n=1}^N\sum_{k=1}^Kz_{nk}\left\{\text{ln}\pi_k+\sum_{i=1}^D[x_{ni}\text{ln}\mu_{ki}+(1-x_{ni})\text{ln}(1-\mu_{ki})]\right\}$$

上式的期望：

$$\mathbb{E}[\text{ln}p(\textbf X, \textbf Z|\pmb\mu,\pmb\pi)]=\sum_{n=1}^N\sum_{k=1}^K\gamma(z_{nk})\left\{\text{ln}\pi_k+\sum_{i=1}^D[x_{ni}\text{ln}\mu_{ki}+(1-x_{ni})\text{ln}(1-\mu_{ki})]\right\}$$

E-step：

$$\gamma(z_{nk})=\mathbb E[z_{nk}]=\frac{\sum_{z_{nk}}z_{nk}[\pi_kp(\textbf x_n|\pmb\mu_n)]^{z_{nk}}}{\sum_{z_{nj}}[\pi_jp(\textbf x_n|\pmb\mu_j)]^{z_{nj}}}=\frac{\pi_kp(\textbf x_n|\pmb\mu_k)}{\sum_{j=1}^K\pi_jp(\textbf x_n|\pmb\mu_j)}$$

M-step：

$$\pmb\mu_k=\overline{\textbf{x}},\ \pi_k=\frac{N_k}{N},\ N_k=\sum_{n=1}^N\gamma(z_{nk}),\ \overline{\textbf{x}}_k=\frac{1}{N_k}\sum_{n=1}^N\gamma(z_{nk})\textbf{x}_n$$

广义EM算法（GEM）

广义期望最大算法（Generalized Expectation Maximization）是EM算法的推广。

EM算法中引入的是隐变量，如果将隐变量理解为隐分布，就得到了广义EM算法，它是在数据不完全或者存在缺失变量的情况下参数估计的迭代算法。

由EM算法中我们可以知道含有隐变量 $z$ 后事件分布规律可以写为：

$$p(\textbf X|\pmb\theta)=\frac{p(\textbf X,\textbf Z|\pmb\theta)}{p(\textbf Z|\textbf X,\pmb\theta)}$$

所以，

$$\begin{align*} &\text{ln}p(\textbf X|\pmb\theta)\\ =&\int q(\textbf Z)\text{ln}p(\textbf X|\pmb\theta)d\textbf Z\\ =&\int q(\textbf Z)\text{ln}\frac{p(\textbf X,\textbf Z|\pmb\theta)q(\textbf Z)}{p(\textbf Z|\textbf X,\pmb\theta)q(\textbf Z)}d\textbf Z\\ =&\int q(\textbf Z)\text{ln}\frac{p(\textbf X,\textbf Z|\pmb\theta)}{q(\textbf Z)}d\textbf Z-\int q(\textbf Z)\text{ln}\frac{p(\textbf Z|\textbf X,\pmb\theta)}{q(\textbf Z)}d\textbf Z\\ =&\mathcal L(q,\pmb\theta)+\text{KL}(q||p) \end{align*}$$

KL距离（Kullback-Leiblier Divergence）也叫相对熵（Relative Entropy）。它衡量的是相同事件空间里的两个概率分布的差异情况。其物理意义是：在相同事件空间里，概率分布 $P(x)$ 对应的每个事件，若用概率分布 $Q(x)$ 编码时，平均每个基本事件（符号）编码长度增加了多少比特。我们用 $D(P||Q)$ 表示KL距离，计算公式：

$$D(P||Q)=\sum_{x\in X}P(x)\text{ln}\frac{P(x)}{Q(x)}$$

可以看到对数似然被分解成了两部分，一个是下界 $\mathcal L(q,\pmb\theta)$，另一个是KL距离。如下图所示：

也就是说，只要令KL距离为0，那么 $\mathcal L(q,\pmb\theta)$ 就等于似然度的值了，此时再找到 $\pmb\theta$ 使得 $\mathcal L(q,\pmb\theta)$ 最大就好。

而最大化这个 $\mathcal L(q,\pmb\theta)$ 比较简单，首先，

$$q(\textbf Z)=p(\textbf Z|\textbf X,\pmb\theta^{\text{old}})$$

将其代入到 $\mathcal L(q,\theta)$，

$$\begin{align*} &\mathcal L(q,\pmb\theta)\\ =&\sum_{\textbf Z}p(\textbf Z|\textbf X,\pmb\theta^{\text{old}})\text{ln}p(\textbf X,\textbf Z|\pmb\theta)-\sum_{\textbf Z}p(\textbf Z|\textbf X,\pmb\theta^{\text{old}})\text{ln}p(\textbf Z|\textbf X,\pmb\theta^{\text{old}})\\ =&\mathcal Q(\pmb\theta,\pmb\theta^{\text{old}})+\text{const} \end{align*}$$

所以只要对 $\mathcal Q(\pmb\theta,\pmb\theta^{\text{old}})$ 最大化就可以了。

总体流程如下：

• E-step：

  

  初始化参数为 $\pmb\theta^{\text{old}}$，计算 $q(\textbf Z)$，

  $$q(\textbf Z)=p(\textbf Z|\textbf X,\pmb\theta^{\text{old}})$$

  把KL设为0，蓝色的线往上移使得 $\mathcal L(q,\theta)=$似然度。

• M-step：

  

  最大化 $\mathcal L(q,\theta)$ ，将参数更新成 $\pmb\theta^{\text{new}}$，

  $$\pmb\theta^{\text{new}}=\text{arg max}\mathcal Q(\pmb\theta,\pmb\theta^{\text{old}})\\ \mathcal Q(\pmb\theta,\pmb\theta^{\text{old}})=\sum_\textbf Zp(\textbf Z|\textbf X,\pmb\theta^{\text{old}})\text{ln}p(\textbf X,\textbf Z|\pmb\theta)$$

  此时对数似然函数的值也得到了上升，这样重复进行下去，可以收敛到一个局部极值。

总的效果图如下图所示：

VBEM

变分EM算法。当遇到比较复杂的、高维的概率图模型时，有可能是树状的、环状的，这时候求其后验分布就比较困难了。这时候有下面两种方法：

1. 采样，采样需要知道模型，之后可以通过积分求出；
2. 变分EM，基于独立性简化计算的假设，将复杂的高维空间分解成低维空间，将其变成一个最优化问题。

集成学习

集成学习是指通过训练多个分类器，然后将这些分类器组合起来，来获得比单个分类器更有的性能（比最好的那个分类器还要好）。如果每个分类器都是同种类型的（比如都是决策树或都是SVM等等），那么这些单个的分类器称为基学习器；如果继承中包含不同类型的分类器，这样的集成是异质的。需要注意的是，这些单个的分类器不一定要很好，只需要比随机猜测好就可以。在我们一般的经验中，如果把好的东西与坏的东西掺杂在一起，那么结果通常是比最坏的要好但比最好的要差一些，但是集成学习可以获得比最好的单一学习器更好的性能。

集成学习的发展历史：

决策树

CART框架

CART使用基尼指数来选择划分。基尼指数反映的是从数据集 $D$ 中随机抽取两个样本，其类别标签不一致的概率。因此基尼指数越少，数据集 $D$ 包含的信息量越少。CART在候选特征集中选取一个特征，使得划分后基尼指数最小。

• 分类举例：

  

  

  根据决策树，将下面的区域分成四块，分别是 $X_2\le0.7,X_1\le-1.4;X_2\le0.7,X_1>-1.4;X_2>0.7,X_1\le-0.6;X_2>0.7,X_1>0.6$。

• 回归举例：

  

  

  $$s^\ast=\text{arg min}_s\left\{\text{min}_{c_1}\sum_{x_i\le s}(y_i-c_1)^2+\text{min}_{c_2}\sum_{x_i>s}(y_i-c_2)^2\right\}$$

模型

分类决策树模型是一种描述对实例进行分类的树形结构，决策树由结点和有向边组成。结点有两种类型：内部结点和叶结点。内部结点表示一个特征或属性，叶节点表示一个类。用决策树分类，从根节点开始，对实例的某一特征进行测试，根据测试结果，将实例分配到其子结点；这时，每一个子结点对应着该特征的一个取值。如果递归地对实例进行测试并分配，直至达到叶结点。最后将实例分到叶结点的类中。

我们可以将决策树看成一个if-then规则的集合，转换成if-then规则的过程：由决策树的根结点到叶结点的每一条路径构建一条规则；路径上内部结点的特征对应着规则的条件，而叶结点的类对应着规则的结论。决策树的路径或其对应的if-then规则集合具有一个重要的性质：互斥并且完备，每一个实例都被一条路径或一条规则所覆盖，而且只被一条路径或一条规则所覆盖。

而分类决策树与回归决策树的区别是：分类构建的决策树的叶结点是一个类别，叶结点中哪个类别的训练样本数量最多，叶结点的类别就取这个类别（由于剪枝早停等优化，叶结点不一定纯，该方法即最大后验概率MAP）。而回归的叶结点是一个值，与目标取值集合是类别的分类问题不同的是，回归的目标取值集合是一个数值集合，对于一个叶结点有一个相对训练样本固定的输出值，这个值通常取叶结点中包含的所有训练样本的均值。特别的，回归在树划分的时候选取的准则是特有的方差度量。这里所谓覆盖是指实例的特征与路径上的特征一致或实例满足规则的条件。

优点

• 高精度、稳定
• 可解释性高
• 非线性
• 适用于解决分类或回归问题
• 既适用于分类的离散变量，也适用于连续的输入和输出变量

特征选择

特征选择，即每一步划分时选择哪个属性进行划分，选择的依据是什么正是解决这个问题的关键。定义每个阶段的划分纯不纯，从上个阶段到以某个属性划分后的下个阶段纯度增加得多还是少，这里列举三种常用准则（另外还有卡方检验等）：信息增益/信息增益率（分类）、基尼增益（分类）、方差度量（回归）。

信息增益/信息增益率（ID3/C4.5）

补充：熵、交叉熵、联合熵、条件熵、相对熵、互信息

• 设 $X$ 是一个取有限个值的离散随机变量，其概率分布为：

$$p(X=x_i)=p_i,\ i=1,2,\dots,n$$

则随机变量 $X$ 的熵定义为：

$$H(X)=-\sum_{i=1}^Np_i\text{log}p_i$$

又记，

$$H(p)=-\sum_{i=1}^Np_i\text{log}p_i=\sum_{i=1}^Np_i\text{log}\frac{1}{p_i}=E\left[\text{log}\frac{1}{p_i}\right]$$

• 对于 $p, q$ 两个分布，交叉熵的定义为：

$$H(p,q)=-\sum_{i=1}^Np_i\text{log}q_i$$

• 联合熵的定义为：

$$\begin{align*} &H(X,Y)\\ =&-\sum_{x\in X}\sum_{y\in Y}p(x,y)\text{log}p(x,y)\\ =&-E_{(X,Y)-p(x,y)}[\text{log}p(X,Y)] \end{align*}$$

条件熵 $H(Y|X)$ 表示在已知随机变量 $X$ 的条件下随机变量 $Y$ 的不确定性。随机变量 $X$ 给定的条件下随机变量 $Y$ 的条件熵（conditional entropy）$H(Y|X)$，定义为 $X$ 给定条件下 $Y$ 的条件概率分布的熵对 $X$ 的数学期望：

$$\begin{align*} &H(Y|X)\\ =&-\sum_{x\in X}p(X=x)H(Y|X=x)\\ =&-\sum_{x\in X}p(x)H(Y|X=x)\\ =&-\sum_{x\in X}p(x)\sum_{y\in Y}p(Y=y|X=x)\text{log}p(Y=y|X=x)\\ =&-\sum_{x\in X}p(x)\sum_{y\in Y}p(y|x)\text{log}p(y|x)\\ =&-\sum_{x\in X}\sum_{y\in Y}p(x,y)\text{log}p(y|x)\\ =&-E_{(X,Y)-p(x,y)}[\text{log}p(Y|X)]\\ =&-E_{p(x,y)}[\text{log}p(Y|X)] \end{align*}$$

当熵和条件熵中的概率由数据估计（特别是极大似然估计）得到时，所对应的熵与条件熵分别称为经验熵和经验条件熵。

• 联合熵与条件熵之间的关系

$$H(X,Y)=H(X)+H(Y|X)$$

证明：

$$\begin{align*} H(X,Y)&=-\sum_{x\in X}\sum_{y\in Y}p(x,y)\text{log}p(x,y)\\ &=-\sum_{x\in X}\sum_{y\in Y}p(x,y)\text{log}p(x)p(y|x)\\ &=-\sum_{x\in X}\sum_{y \in Y}p(x,y)\text{log}p(x)-\sum_{x\in X}\sum_{y\in Y}p(x,y)\text{log}p(y|x)\\ &=-\sum_{x\in X}p(x)\text{log}p(x)-\sum_{x\in X}\sum_{y\in Y}p(x,y)\text{log}p(y|x)\\ &=H(X)+H(Y|X) \end{align*}$$

• 相对熵（Kullack-Leibler距离）：度量两个分布之间的距离、

$$D(p||q)=\sum_{x\in X}p(x)\text{log}\frac{p(x)}{q(x)}=E_{p(x)}\left[\text{log}\frac{p(X)}{q(X)}\right]$$

• 互信息（Mutual Information）

$$\begin{align*} I(X;Y)&=\sum_{x\in X}\sum_{y\in Y}p(x,y)\text{log}\frac{p(x,y)}{p(x)p(y)}\\ &=E_{p(x,y)}\left[\text{log}\frac{p(X,Y)}{p(X)p(Y)}\right]\\ &=D(p(x,y)||p(x)p(y)) \end{align*}$$

• 互信息与熵 $$\begin{align*} I(X;Y)&=\sum_{x\in X}\sum_{y\in Y}p(x,y)\text{log}\frac{p(x,y)}{p(x)p(y)}\\ &=\sum_{x\in X}\sum_{y\in Y}p(x,y)\text{log}\frac{p(x|y)}{p(x)}\\ &=-\sum_{x\in X}\sum_{y\in Y}p(x,y)\text{log}p(x)+\sum_{x\in X}\sum_{y\in Y}p(x,y)\text{log}p(x|y)\\ &=H(X)-H(X|Y) \end{align*}$$

• $I(X;Y)=H(X)-H(X|Y)$
• $I(X;Y)=H(Y)-H(Y|X)$
• $I(X;Y)=H(X)+H(Y)-H(X,Y)$

---

对于一个样本数量为 $D$，类别数量为 $N$ 的集合 $D$，该集合的熵定义为：

$$H(D)=-\sum_{i=1}^Np_i\text{log}p_i$$

将该集合以某个属性 $A$ 划分一次，划分后的总熵为：

$$H(D|A)=\sum_{i\in\{L,R\}}\frac{|D^i|}{|D|}H(D^i)$$

假设将该集合以某个属性 $A$ 划分一次，得到两个样本数量分别为 $C_1$ 和 $C_2$ 的集合，则划分后的总熵为

$$H(D|A)=-\frac{C_1}{D}\sum_{i=1}^{N_1}p_i\text{log}p_i-\frac{C_2}{D}\sum_{i=1}^{N_2}p_i\text{log}p_i$$

即两个集合各自求本集合的熵，再计算加权平均和作为整个的熵，对于多分类，也是同理进行加权平均。那么，信息增益定义如下：

$$g(D;A)=H(D)-H(D|A)$$

• 决策树中的信息增益等价于训练数据集中类与特征的互信息。
• 熵随着分支的进行应该越来越小，因此以某个特征进行划分，不同的特征往往具有不同的信息增益，信息增益大的特征具有更强的分类能力。

假设某个属性 $X$，分裂前熵=10.0，分裂后=6.0，I=10.0-6.0=4.0；某个属性 $Y$，分裂前熵=4.0，分裂后=0.2，I=4.0-0.2=3.8。由上面这个例子可以看出，用信息增益的方式选择，存在偏向于选择取值比较多的特征的问题，因此便出现了信息增益率来作为一种补偿方案。

信息增益率的定义为：

$$IGR=\frac{H(D)-\sum_{i\in\{L,R\}}\frac{|D^i|}{|D|}H(D^i)}{-\sum_{i\in\{L,R\}}\frac{|D^i|}{|D|}H(D^i)}$$

Gini基尼指数（CART）

基尼指数的定义：

$$Gini(p)=\sum_{k=1}^Kp_k(1-p_k)=1-\sum_{k=1}^kp_k^2$$

类似于信息增益，将该集合以某个属性 $A$ 划分一次，划分后的Gini指数为：

$$Gini(D|A)=\frac{|D_1|}{|D|}Gini(D_1)+\frac{|D_2|}{|D|}Gini(D_2)$$

Gini指数 vs 熵

• Gini 指数的计算不需要对数运算，更加高效；
• Gini 指数更偏向于连续属性，熵更偏向于离散属性。

方差

利用方差定义数据之间的关系时显然有：方差越小，数据之间的关系越紧密，即这些数据之间相似度高。

假设现在按 $s$ 为间隔划分成两类：

$$\begin{align*} &R_1(j,s)=\{x|x_j\le s\},\\ &R_2(j,s)=\{x|x_j>s\}. \end{align*}$$

方差度量准则表达式为：

$$\text{min}_{j,s}\left[ \underbrace{\text{min}_{c_1}\sum_{x_i\in R_1(j,s)}(y_i-c_1)^2}_{\hat c_1=ave(y_i|x_i\in R_1(j,s))}+\underbrace{\text{min}_{c_2}\sum_{x_i\in R_2(j,s)}(y_i-c_2)^2}_{\hat c_2=ave(y_i|x_i\in R_2(j,s))}\right]$$

即每次划分要使划分的两类的均方差各自达到最小且两者的和也要最小。

树的生成

结点的取值情况？

最常见的决策树模型ID3是多叉树，即属性由几类取值就分为几叉，但二叉树更简单，CART模型就是采用的二叉树，实际上任何多叉树都可以转化为二叉树，即对于多于两个取值的属性，通过yes/no询问，一次只划分出某类和非某类两类，然后继续对非某类进行更细致的划分。

采用哪个属性进行测试？

• 不纯度准则

不纯度准则选择上，分类和回归各有特有准则，由于回归通常就用方差度量，这里对比一下分类问题使用的熵、基尼指数和基本的误分类不纯度（分类误差率）之间的关系。

其中，误分类的不纯度定义为：

$$i(N)=1-\text{max}_jp(C_j|N)$$

熵定义为：

$$i(N)=-\sum_jp(C_j|N)\text{log}_2p(C_j|N)$$

Gini不纯度定义为：

$$i(N)=\frac{1}{2}\left(1-\sum_jp^2(C_j|N)\right)$$

属性值类型

当属性值为标称值时，即固定的几个类别，则只需要遍历所有属性值，就知道该属性的所有类别。示例中年龄由三个类别，一目了然。

属性值为实数值时，对于该属性的类别，并不明显，此时一种解决方法是，先对该属性已经出现的值进行排序，然后在标签发生变化的地方进行切分。示例中体重是实数值，将训练样本中的体重值排序，注意该示例中目标分类（标签）为蓝色和红色两类，因此对排序序列从左向右看，在目标分类发生变化处切分，使同类的直接切到一个判断区间，这样在针对该属性进行划分时能尽可能地减少分支。

单属性vs多属性

即每次划分时判断是只考虑一个属性，还是同时考虑多个属性。对于如下图所示地训练数据分布，左边为单属性划分，右边为双属性划分，显然，双属性更加自然合理。但是多属性同时考虑效率并不高，现在换个角度看数据分布，有这样的分布，说明其 $x,y$ 两个维度直接存在纠缠关系，一次仅仅用一个属性（直线 $y=k$ 或者 $x=k$）根本分不开，此时如果旋转坐标系（斜着看）去相关性，又可以通过单属性进行划分了。也就说，需要用多属性同时进行划分地，可以对数据进行一定地映射变换去相关性，再用单属性进行划分，但是不管是去相关性还是直接多属性划分，没有绝对的谁更容易。

何时为叶结点

由于噪声地影响，把所有的数据都分类正确的树不一定是一棵好树，如下图所示，随着决策树越来越复杂，对训练集是越来越准确，但是泛化能力又变差了。

解决这样地问题，可以加入树分裂停止准则和剪枝进行优化。前者即预剪枝，后者即后剪枝。

预剪枝：树分裂时就实施措施。

1. 当树的深度达到预先设定的阈值停止当前分支的分裂；
2. 当分完之后，支持的数据变少，即事先划分验证集测试正确率。

后剪枝：树建好以后，删除掉一些无用的分支，即删除该分支有或无对测试结果没太大的影响的分支。

1. 降低错误剪枝 REP（Reduced Error Pruning）

   划分一个测试数据集来纠正它。该剪枝方法考虑将书上的每个节点作为修剪的候选对象，决定是否修剪这个结点有如下步骤组成：

   1. 删除以此结点为根的子树；
   2. 使其成为叶子结点；
   3. 赋予该结点关联的训练数据的最常见分类；
   4. 当修剪后的树对于验证集合的性能不会比原来的树差时，才真正删除该结点。

优点：

1. 最简单的后剪枝方法之一；
2. 计算复杂度是线性的；

3. 和原始决策树相比，修剪后的决策树对未来新事例的预测偏差较小。

   缺点：

4. 在数据量较少的情况下很少应用。REP方法趋于过拟合（overfitting），这是因为训练数据集中存在的特性在剪枝过程中都被忽略了，当剪枝数据集比训练数据集小得多时，这个问题特别值得注意。

1. 悲观错误剪枝 PEP（Pessimistic Error Pruning）

   该方法克服了REP需要独立剪枝集的缺点，基于训练数据的误差评估，因此不用单独找剪枝数据集。但训练数据也带来错分误差偏向于训练集，所以引入了统计学上连续修正的概念弥补REP中的缺陷，在评价子树的训练错误公式中添加了一个常数，假定每个叶子结点都自动对实例的某个部分进行错误的分类。

设原始决策树 $T$，叶结点 $z$，$z$ 结点训练实例个数为 $n_z$，其中错分个数为 $e_z$。

定义误差率且增加连续性校正：$pe_z=(e_z+0.5)/n_z$，

相应的误差率：$E(T)=\sum_{t\in T}\frac{e(t)+0.5}{N(t)}$，

剪枝后的误差率：$E(T’)=\sum_{t\in T,\text{except }K}\frac{e(t)+0.5}{N(t)}$

标准差错误：$SE(E(T’))=\sqrt{\frac{E(T’)\ast(N_t-E(T’))}{N_t}}$，其中，$N_t$ 是当前训练数据量。

剪枝条件：$SE(E(T’))\le E(T’)-E(T)$，符合此条件，剪掉 $t$ 。

优点：

1. PEP方法被认为是当前决策树后剪枝方法中精度较高的算法之一；
2. PEP 方法不需要分离的剪枝数据集, 这对于事例较少的问题非常有利；
3. 它的计算时间复杂性也只和未剪枝树的非叶节点数目成线性关系。

缺点：PEP是唯一使用自顶向下剪枝策略的事后剪枝方法, 这种策略会带来与事前剪枝方法出 现的同样问题, 那就是树的某个节点会在该节点的子孙根据同样准则不需要剪裁时也会被剪裁。

1. 基于错误剪枝 EBP（Error-based Pruning）

   1. 计算叶结点的错分样本率估计的置信区间上限 $U$；

   2. 计算叶结点的预测错分样本数

      • 叶结点的预测错分样本数=到达该叶结点的样本数$\times$该叶结点的预测错分样本率 $U$

   3. 判断是否剪枝及如何剪枝

      • 分别计算三种预测错分样本数：

        • 计算子树 $t$ 的所有叶结点预测错分样本数之和，记为 $E_1$
        • 计算子树 $t$ 被剪枝以叶结点代替时的预测错分样本数，记为 $E_2$
        • 计算子树 $t$ 的最大分支的预测错分样本数，记为 $E_3$

      • 比较 $E_1,E_2,E_3$，如下：

        • $E_1$ 最小时，不剪枝
        • $E_2$ 最小时，进行剪枝，以一个叶结点代替 $t$
        • $E_3$ 最小时，采用“嫁接”（grafting）策略，即用这个最大分支代替 $t$

2. 代价-复杂度剪枝 CCP（Cost-complexity Pruning）

   CART剪枝法，该算法为子树 $T_t$ 定义了代价（cost）和复杂度（complexity），以及一个可由用户设置的衡量代价与复杂度之间关系的参数 $α$，其中，代价指在剪枝过程中因子树 $T_t$ 被叶结点替代而增加的错分样本，复杂度表示剪枝后子树 $T_t$ 减少的叶结点数，$α$ 则表示剪枝后树的复杂度降低程度与代价间的关系，定义为：

   $$\alpha=\frac{R(t)-R(T_t)}{|N_1|-1}$$

   其中，$|N_1|$ 为子树 $T_t$ 中的叶结点数，$R(t)$ 为结点 $t$ 的错误代价，计算公式为 $R(t)=r(t)\times p(t)$，$r(t)$ 为结点 $t$ 的错分样本率，$p(t)$ 为落入结点 $t$ 的样本占所有样本的比例，$R(T_t)$ 子树 $T_t$ 错误代价，计算公式为 $R(T_t)=\sum R(i)$ ，$i$ 为子树 $T_t$ 的叶结点。

算法分为两个步骤：

1. 对于完全决策树 $T$ 的每个非叶结点计算 $α$ 值，循环剪掉具有最小 $α$ 值的子树，直到剩下根节点。在该步可得到一系列的剪枝树 $\{T_0,T_1,T_2,\dots,T_m\}$，其中 $T_0$ 为原有的完全决策树，$T_m$ 为根结点，$T_{i+1}$ 为对 $T_i$ 进行剪枝的结果；

2. 从子树序列中，根据真实的误差估计选择最佳决策树。

缺点：生成子树序列 $T(α)$ 所需要的时间和原决策树非叶节点的关系是二次的，这就意味着如果非叶节点的数目随着训练样本数线性增加，则CCP方法的运行时间和训练样本数的关系也是二次的。这就比其它剪枝方法所需要的时间长得多，因为其它剪枝方法的运行时间和非叶节点的关系是线性的。

1. 最小错误剪枝 MEP（Minimum Error Pruning）

   基本思路：采用自底向上的方式，对于树中每个非叶结点，

   1. 计算该结点的误差 $E(t)$；

      $$E_r(t)=\frac{n(t)-n_e(t)+k-1}{n(t)+k}$$

      其中，$n(t)$ 为结点 $t$ 中的样本总数，$n_e(t)$ 为 $t$ 中主类的样本数目，$k$ 为类数目。

   2. 计算该结点的各个分支的加权累计误差 $E’(t)$；

   3. 若 $E(t)>E’(t)$，保留；否则，剪枝。

优点：

1. MEP方法不需要独立的剪枝数据集，无论是初始版本，还是改进版本，在剪枝过程中，使用的信息都来自于训练样本集；
2. 它的计算时间复杂性也只和未剪枝树的非叶节点数目成线性关系。

缺点：

1. 类别平均分配的前提假设现实几率不大；
2. 对 $K$ 太过敏感。

比较项\剪枝方法	REP	PEP	CCP	MEP
独立剪枝集	需要	不需要	不需要（交叉验证）	不需要
剪枝方式	自底向上	自顶向下	自底向上	自底向上
误差估计	利用剪枝集	使用连续性校正	使用交叉验证或标准误差	基于概率估计
计算复杂性（$n$ 为非叶结点个数）	$O(n)$	$O(n)$	$O(n^2)$	$O(n)$

处理属性值丢失结点

• 丢弃这些样本；
• 赋给它结点 $t$ 所对应的训练实例中该属性的最常见值；
• 代理结点法：缺失的属性在测试阶段用其他的属性代替。

C4.5—缺失值

例如 $(X,y)$ 是样本集 $S$ 中的一个训练实例，$X=(F_1,F_2,\dots,F_n)$。但是其属性 $F_i$ 的值未知。

处理策略：为 $F_i$ 的每个可能值赋予一个概率。例如，给定一个布尔属性 $F_i$，如果结点 $t$ 包含6个已知 $F_i=1$ 和4个 $F_i=0$ 的实例，那么 $F_i=1$ 的概率是0.6，而 $F_i=0$ 的概率是0.4。于是，实例 $x$ 的60%被分配到 $F_i=1$ 的分支，40%被分配到另一个分支。这些片段样例（fractional examples）的目的是计算信息增益，

模型选择

1. CART：支持分类和回归，分类采用基尼指数，回归使用方差度量，二叉树，支持连续值处理，支持缺失值处理，利用独立的验证数据集进行剪枝。
2. ID3：仅分类，采用信息增益准则，多叉树，其他不支持。
3. C4.5：在ID3的基础上，仅分类，采用信息增益率准则，多叉树，支持连续值处理，利用基于错误剪枝方法进行剪枝。
4. C5.0：在C4.5基础上加入了Boosting技术，速度更快，内存使用更加有效，能生成更小的决策树，但算法细节并没有给公开太多，属于商业机密。

随机森林

随机森林最基本的思想为，在训练阶段树的生成过程中，制定随机部分，随机地生成一颗决策树，通过重复这样的过程，随机生成若干决策树组成森林，然后在预测的时候，综合所有决策树的判断结果得到最终预测结果。随机度越高，树的相关性越低，综合考虑结果的效果就越好。综合结果有加法模型和乘法模型，其中乘法模型有一个痛点就是如果其中某棵树预测值为0，那么整个结果就变为了0。

Randomized Decision Trees

随机进行属性选择。

特点：

• 不再寻找全局最优的分裂结点，随机分裂结点，训练更快，使用更简单的二值特征测试；
• 选择一个包含 $k$ 个属性的集合，从中依据某个准则函数来挑选最优的分裂结点，比如，最大化信息增益；
• 随机树的效果不如决策树，通过构建多棵树来弥补；
• 各棵树的投票结果求平均。

应用：

1. OCR；

2. 手写体识别；

3. 人脸识别；

4. 特定物体检测；

   

Random Forest

算法过程：

1. 有放回地随机选择 $N$ 个样本；
2. 随机从 $M$ 个属性中选取出 $m$ 个属性（$m<<M,m=\sqrt M$），依据某种策略（比如说信息增益）来选择一个属性作为分裂结点；
3. 按照步骤2，构建一棵决策树；
4. 按照步骤1~3建立大量的决策树，组成随机森林。

简单地说，行随机、列随机。

优点：

1. 比较适合做多分类问题；
2. 训练和预测速度快，能够有效地处理大的数据集；容易并行化；
3. 数据集中有大比例的数据缺失时，仍然可以保持精度不变；
4. 能够检测到特征之间的相互影响以及重要性程度；
5. 不易出现过度拟合。

分类森林：

回归森林：

随机度：

• 随机度高，树的相关性低；
• 随机度低，树的相关性高。

树的集成：

• 加法模型
• 乘法模型

应用：

1. Kinect上的姿态估计；

   

   

1. Hough Forest for object detection

   

   

   

2. 人脸关键点标定

   

   

Extremely Randomized trees

动机：随机森林中，每个维度上的最优切割点的方差较大。

算法过程：

1. 随机选择 $k$ 个属性值，计算出属性值的最小值和最大值；
2. 在 $[\text{min,max}]$ 区间上随机生成 $k$ 个切分点；
3. 选择 score 最好的切分点。

与随机森林的不同点：

1. 无bootstrap，使用全部的训练样本；
2. 在选择结点的分裂位置时，随机。

Ferns

算法：ET上各层上的结点都一样。

• 每个结点想象成一个特征；
• 其输出值为0或者1，表示该特征有无出现；
• 在图像块中，可以比较像素点值的大小

$$f_j=\left\{\begin{aligned}&1,&&\text{if }I(d_{j,1})<I(d_{j,2})\\&0,&&\text{otherwise}\end{aligned}\right.$$

应用：

1. 3D object detection
2. SLAM

Boosting

提升（Boosting）是一种常用得统计学习方法，在分类问题中，它通过改变训练样本的权重，学习多个分类器（一般是弱分类器），并将这些分类器线性组合，最终提高分类器的性能。而针对于这种提升方法而言，需要回答两个问题，一是在每一轮如何改变训练样本的权值或概率分布；二十如何将弱分类器组合成一个强分类器。

AdaBoost

Adaboost属于Boosting一种，针对第一个问题，Adaboost的做法是提高哪些被前一轮弱分类器错误分类样本的权值，而降低那些被正确分类的样本的权重。从而使得那些被错误分类的样本由于其权值被加大而受到后一轮弱分类器的更多的关注。针对第二个问题，Adaboost采取加权多数表决的方法，加大分类误差率小的弱分类器的权重，使其在最终的分类器表决中起较大作用，减小分类误差率大的弱分类器的权重，使其在表决中其较小的作用。

算法流程

算法流程如下：

1. 首先初始化各个样本的权重，为 $\frac{1}{m}$，$m$ 为样本数量，有些书本可能会这样初始化：每个正样本权重为 $\frac{1}{2u}$，每个负样本权重为 $\frac{1}{2v}$，$u+v$ 等于样本总数。这对于得到最终的强分类器并没有影响。

2. 开始进行迭代，在每一轮的迭代中：

   1. 寻找到一个弱分类器使得在当前样本权重分布下的错误率最小；
   2. 当确定了这样一个弱分类器后，计算错误率 $\varepsilon_j$ ；
   3. 根据误差率计算此弱分类器在最终的强分类器中所占的权重 $\alpha_t$；
   4. 根据公式更新训练样本的权值分布；
   5. 将此弱分类器加入到当前迭代次数下的强分类器中。

3. 直到迭代次数达到指定数值或者强分类器在训练集上的误差达到指定范围停止迭代。

例子

1. 第一次迭代我们给所有样本相同的权重，找到当前权重分布下，分类误差率最小的一个弱分类器；
2. 第二次迭代，经过上一轮弱分类器之后，有三个样本被分错，加大其权重，找到在当前权重分布下，误差率最小的一个弱分类器；
3. 以此类推；
4. 最终，将得到的一系列弱分类器加权组合起来就得到了一个强分类器（实黑线）。

应用于目标检测

软件资源

• Ensemble Learning Framework
• AdaBoost Matlab Toolbox
• Xgboost，陈天奇
• LightGBM

策略

AdaBoost算法是前向分步加法算法的特例，模型是由基本分类器组成的加法模型，损失函数是指数函数。

模型：

$$f(x)=\sum_{m=1}^M\alpha_mG_m(x)$$

策略：

$$L(y,f(x))=\text{exp}[-yf(x)]$$

$G_m(x)$ 在训练集上的分类误差率 $e_m$ 为：

$$e_m=\sum_{i=1}^Nw_{mi}I(G_m(x_i)\neq y_i)$$

$G_m(x)$的系数 $\alpha_m$ 为：

$$\alpha_m=\frac{1}{2}\text{log}\frac{1-e_m}{e_m}$$

在证明前，先介绍前向分布算法。给定一个加法模型：

$$f(x)=\sum_{m=1}^M\beta_mb(x;\gamma_m)$$

其中，$b$ 为基函数，$\gamma_m$ 为基函数的参数，$\beta_m$ 为基函数的系数。在给定训练数据基损失函数 $L(y,f(x))$ 的条件下，学习加法模型 $f(x)$ 成为经验风险最小化即损失函数最小化问题：

$$\text{min}_{\beta_m,\gamma_m}\sum_{i=1}^NL(y_i,\sum_{m=1}^M\beta_Mb(x_i,\gamma_m))$$

通常这是一个复杂的优化问题，一次性求解所有基函数的参数以及其系数比较困难。前向分布算法（forward stagewise algorithm）求解这一优化问题的思路是：因为学习的是加法模型，如果能够从后向前，每一步只学习一个基函数及其系数，逐步逼近其目标函数式，那么就可以简化优化的复杂度。具体地，每步只需优化如下损失函数：

$$\text{min}_{\beta,\gamma}\sum_{i=1}^NL(y_i,\beta b(x_i,\gamma))$$

了解了前向分布算法，下面就要证明上述结论了。从我们的模型入手：

$$f(x)=\sum_{m=1}^M\alpha_mG_m(x)$$

假设 $f_0(x)$ 表示第0此迭代后的模型，$f_m(x)$ 表示第 $m$ 次迭代后的模型，有：

$$\begin{align*} &f_0(x)=0\\ &f_1(x)=f_0(x)+\alpha_1G_1(x)=\sum_{m=1}^1\alpha_mG_m(x)\\ &f_2(x)=f_1(x)+\alpha_2G_2(x)=\alpha_1G_1(x)+\alpha_2G_2(x)=\sum_{i=1}^2\alpha_iG_i(x)\\ &\dots\\ &f_{m-1}(x)=f_{m-2}(x)+\alpha_{m-1}G_{m-1}(x)=\alpha_1G_1(x)+\dots+\alpha_{m-1}G_{m-1}(x)=\sum_{i=1}^{m-1}\alpha_iG_i(x)\\ &f_m(x)=f_{m-1}(x)+\alpha_mG_m(x)=\sum_{i=1}^m\alpha_iG_i(x) \end{align*}$$

那么经过 $m$ 次迭代所产生的最终模型 $f_m(x)$，表示为前 $m-1$ 次迭代产生的模型 $f_{m-1}(x)$ 加上第 $m$ 次迭代产生的加权基函数 $\alpha_mG_m(x)$。代入损失函数：

$$L(y,f(x))=\text{exp}(-yf(x))$$

即第 $m$ 轮迭代的基函数以及其系数必然满足下述极小化式：

$$(\alpha_m,G_m(x))=\text{arg min}_{\alpha,G}\sum_{i=1}^N\text{exp}[-y_i(f_{m-1}(x_i))+\alpha G(x_i)]$$

对上式变形：

$$\begin{align*} &\sum_{i=1}^N\text{exp}[-y_i(f_{m-1}(x_i)+\alpha G(x_i))]\\ =&\sum_{i=1}^N\text{exp}[-y_if_{m-1}(x_i)]\times\text{exp}[-y_i\alpha G(x_i)]\\ =&\sum_{i=1}^N\overline w_{mi}\times \text{exp}[-y_i\alpha G(x_i)] \end{align*}$$

其中，$\overline w_{mi}=\text{exp}[-y_if_{m-1}(x_i)]$，注意到这一项是与 $\alpha_m$ 以及 $G_m(x)$ 无关的。

即：

$$(\alpha_m,G_m(x))=\text{arg min}_{\alpha, G}\sum_{i=1}^N\overline w_{mi}\times \text{exp}[-y_i\alpha G(x_i)]$$

进一步展开，

$$\begin{align*} &\sum_{i=1}^N\overline w_{mi}\times\text{exp}[-y_i\alpha G(x_i)]\\ =&\sum_{i=1}^N\overline w_{mi}\times e^{-\alpha}\times I\{y_i=G(x_i)\}+\sum_{i=1}^N\overline w_{mi}\times e^\alpha \times I\{y_i\neq G(x_i)\}\\ =&e^{-\alpha}\sum_{i=1}^N\overline w_{mi}\times I\{y_i=G(x_i)\}+e^\alpha \sum_{i=1}^N\overline w_{mi}\times I\{y_i\neq G(x_i)\}\\ =&e^{-\alpha}\sum_{i=1}^N\overline w_{mi}\times I\{y_i=G(x_i)\}+e^\alpha \sum_{i=1}^N\overline w_{mi}\times I\{y_i\neq G(x_i)\}\\ +&e^{-\alpha}\sum_{i=1}^N\overline w_{mi}\times I\{y_i\neq G(x_i)\}-e^{-\alpha}\sum_{i=1}^N\overline w_{mi}\times I\{y_i\neq G(x_i)\}\\ =&e^{-\alpha}\sum_{i=1}^N\overline w_{mi}+(e^\alpha-e^{-\alpha})\sum_{i=1}^N\overline w_{mi}\times I\{y_i\neq G(x_i)\} \end{align*}$$

因为 $\overline w_{mi}$ 与 $\alpha_m$ 以及 $G_m(x)$ 无关，可以看到所求的 $G_m^\ast(x)$ 满足下列极小化式，

$$G_m^\ast(x)=\text{arg min}_G\sum_{i=1}^N\overline w_{mi}I(y_i\neq G(x_i))$$

假设我们已经求得了 $G_m^\ast(x)$，那接下来求它的系数 $\alpha_m^\ast$，对 $\alpha$ 求偏导，令其等于0：

$$\begin{align*} &\frac{\partial}{\partial \alpha}\left(e^{-\alpha}\sum_{i=1}^N\overline w_{mi}+(e^\alpha-e^{-\alpha})\sum_{i=1}^N\overline w_{mi}\times I\{y_i\neq G(x_i)\}\right)\\ =&-e^{-\alpha}\sum_{i=1}^N\overline w_{mi}+(e^\alpha+e^{-\alpha})\sum_{i=1}^N\overline w_{mi}\times I\{y_i\neq G(x_i)\}=0 \end{align*}$$

即得：

$$\frac{e^\alpha+e^{-\alpha}}{e^{-\alpha}}=\frac{\sum_{i=1}^N\overline w_{mi}}{\sum_{i=1}^N\overline w_{mi}\times I\{y_i\neq G(x_i)\}}\Rightarrow\alpha_m^\ast=\frac{1}{2}\text{log}\frac{1-e_m}{e_m}$$

其中，

$$e_m=\frac{\sum_{i=1}^N\overline w_{mi}\times I\{y_i\neq G(x_i)\}}{\sum_{i=1}^N\overline w_{mi}}$$

对 $e_m$ 进一步化简：

$$\begin{align*} e_m=&\frac{\sum_{i=1}^N\overline w_{mi}\times I\{y_i\neq G(x_i)\}}{\sum_{i=1}^N\overline w_{mi}}\\=&\sum_{i=1}^N\left(\frac{\overline w_{mi}}{\sum_{i=1}^N\overline w_{mi}}\right)I(y_i\neq G_m(x_i))\\=&\sum_{i=1}^Nw_{mi}I(y_i\neq G_m(x_i)) \end{align*}$$

最后推导出训练集样本权值的更新公式，由模型：

$$f_m(x)=f_{m-1}(x)+\alpha_mG_m(x)$$

以及权值：

$$\overline w_{mi}=\text{exp}(-y_if_{m-1}(x_i))$$

可以得到：

$$\overline w_{m+1,i}=\overline w_{m,i}\text{exp}(-y_i\alpha_mG_m(x))$$

决策规则为：

$$G(x)=\text{sgn}(f(x))=\text{sgn}\left(\sum_{m=1}^M\alpha_mG_m(x)\right)$$

损失函数与类概率的对数几率之间的关系

上述的损失函数可以写为：

$$f^\ast(x)=\text{arg min}_{f(x)}E_{Y|x\sim D}[e^{-Yf(x)}]$$

其中期望为：

$$E_{Y|x\sim D}[e^{-Yf(x)}]=p(Y=1|x)e^{-f(x)}+p(Y=-1|x)e^{f(x)}$$

对其求偏导，

$$\begin{align*} &\frac{\partial E_{Y|x\sim D}[e^{-Yf(x)}]}{\partial f(x)}=-p(Y=1|x)e^{-f(x)}+p(Y=-1|x)e^{f(x)}=0\\ \Rightarrow&f^\ast(x)=\frac{1}{2}\text{ln}\frac{p(Y=1|x)}{p(Y=-1|x)}\\ \end{align*}$$

所以，

$$\begin{align*} G(x)&=\text{sgn}(f(x))\\ &=\text{sgn}\left(\frac{1}{2}\text{ln}\frac{p(Y=1|x)}{p(Y=-1|x)}\right)\\ &=\left\{\begin{aligned}&1,&&p(Y=1|x)>p(Y=-1|x)\\&-1,&&p(Y=1|x)\le p(Y=-1|x)\end{aligned}\right.\\ &=\text{arg max}_{y\in\{-1,1\}}p(y|x) \end{align*}$$

可以看出其达到了贝叶斯最优错误率。

损失函数回顾

理想的0-1损失函数

$$L(y,f_\theta(x))=\left\{\begin{aligned}&0,&&y=f_\theta(x)\\&1,&&y\neq f_\theta(x)\end{aligned}\right.=I\{y\neq f_\theta(x)\}$$

这类损失函数并不常用，因为其有间断点，不可导。

基于最小二乘法的分类

$$\begin{align*} L(y,f_\theta(x))&=(y-f_\theta(x))^2\\ &=y^2\left(1-\frac{f_\theta(x)}{y}\right)^2\\ &=\left(1-\frac{f_\theta(x)}{y}\right)^2\\ &=(1-yf_\theta(x))^2\\ &=(1-m)^2 \end{align*}$$

这类损失函数的优点是有解析解，但它的缺点也很明显，就是对”坏点”非常敏感，并且当 $f(x)$ 的值很大时损失也很大，“妒忌心”很重。

SVM的Hinge损失函数

损失函数为：

$$\text{min}_{w,b}\sum_i^N[1-y_i(wx_i+b)]_++\lambda||w||^2$$

变形，令

$$[1-y_i(wx_i+b)]_+=\xi_i$$

于是，原式就变成了，

$$\text{min}_{w,b}\sum_i^N\xi_i+\lambda||w||^2$$

若进一步表示成，

$$\text{min}_{w,b}\frac{1}{C}\left(\frac{1}{2}||w||^2+C\sum_i^N\xi_i\right)$$

类似于下式：

$$\frac{1}{m}\sum_{i=1}^ml(wx_i+b,y_i)+||w||^2$$

前半部分中的就是Hinge损失函数，而后面相当于L2正则项，说明SVM损失函数是正则化了的Hinge损失函数。这类损失函数线性惩罚误分类的点，鲁棒性较强。

AdaBoost的指数损失函数

这类损失函数连续，可导，接近理想的0-1损失，对太正确的点惩罚低，收敛快，误分类点的损失指数增长，鲁棒性较差。

Logistic回归的指数损失函数

Logistic回归的指数损失，m>0时，类似指数损失；m<0时，损失是线性增长的。

AdaBoost的损失函数可以改成Logistic Regression的损失函数，从而更加鲁棒，这样的算法称作 GentleBoost。在每次迭代时，基于最小二乘去做一个加权回归，最后所有回归函数的和作为最终的分类器。

提升树（Boosting Trees）

提升树是以分类树或回归树为基本分类器的提升方法。提升树被认为是统计学习中性能最好的方法之一。

模型

提升方法实际采用加法模型（即基函数的线性组合）与前向分步算法。以决策树为基函数的提升方法成为提升树。对分类问题决策树是二叉分类树，对回归问题决策树是二叉回归树。只有一个条件的基本分类器可以看做是由一个根结点直接连接两个叶结点简单决策树，即所谓的决策树桩（decision stump）。提升树模型可以表示为决策树的加法模型：

$$f_M(x)=\sum_{m=1}^MT(x;\Theta_m)$$

其中，$T(x; \Theta_m)$ 表示决策树；$\Theta_m$ 为决策树的参数；$M$ 为树的个数。提升树中树之间没有权重，或者说它们的权重都是一样的，树之间是独立的，训练样本之间也没有权重的概念，这是提升树和随机森林、AdaBoost之间的区别。

算法

提升树算法采用前向分步算法。首先确定初始提升树 $f_0(x)=0$，第 $m$ 步的模型是：

$$f_m(x)=f_{m-1}(x)+T(x,\Theta_m)$$

其中，$f_{m-1}(x)$ 为当前模型，通过经验风险最小化确定下一决策树的参数 $\Theta_m$：

$$\hat \Theta=\arg \min_{\Theta_m}\sum_{i=1}^NL(y_i,f_{m-1}(x)+T(x_i;\Theta_m))$$

最终得到的模型为：

$$f_M(x)=\sum_{m=1}^MT(x;\Theta_m)$$

由于树的线性组合可以很好地拟合训练数据，即使数据中的输入与输出之间关系很复杂也是如此。所以提升树是一个表现很好的学习算法。

不同问题中的提升树学习算法的主要区别在于使用的损失函数不同。例如用平方损失函数的回归问题，用指数损失函数的分类问题，以及用一般损失函数的一般决策问题。

回归问题

回归问题采用平方误差损失函数，损失变为：

$$\begin{align*} &L(y_i,f_{m-1}(x_i)+T(x_i;\Theta_m))\\ =&[y_i-f_{m-1}(x_i)-T(x_i;\Theta_m))]^2\\ =&[r_i-T(x_i;\Theta_m)]^2 \end{align*}$$

其中，$r_i=y_i-f_{m-1}(x_i)$，是当前模型拟合数据的残差（residual）。所以，对回归问题的提升树算法来说，只需简单地拟合当前模型的残差。

梯度提升树（GBDT）

梯度提升树，Gradient Boosting，其学习流程与提升树类似，只是不再使用残差作为新的训练数据而是使用最速下降的近似方法，即使用损失函数的梯度作为新的新的训练数据的 $y$ 值，具体的来说就是使用损失函数对 $f(x)$ 求梯度然后带入 $f_{m-1}(x)$ 。$f(x)=f_{m-1}(x)+T(x;\Theta_m)$，损失为：

$$\begin{align*} &L(y,f(x))\\ \neq&L(y,f_{m-1}(x))+\left.\frac{\partial L(y,f(x))}{\partial f(x)}\right|_{f(x)=f_{m-1}(x)}\times(f(x)-f_{m-1}(x)) \end{align*}$$

负梯度为：

$$T(x;\Theta)=-\left.\frac{\partial L(y,f(x))}{\partial f(x)}\right|_{f(x)=f_{m-1}(x)}$$

GDBT与提升树之间的关系：

提升树模型每一次的提升都是靠上次的预测结果与训练数据的label值差值作为新的训练数据进行重新训练，GDBT则是将残差计算替换成了损失函数的梯度方向，将上一次的预测结果带入梯度中求出本轮的训练数据，也就是说这两种模型在生成新的训练数据时采用了不同的方法。在使用平方误差损失函数和指数损失函数时，提升树的残差求解比较简单，但是在使用一般的损失误差函数时，残差求解起来不是那么容易，所以就使用损失函数的负梯度在当前模型的值作为回归问题中残差的近似值。

上述给出的提升树是损失函数采用平方损失的版本，当采用其他损失函数的时候，可以将其细分为若干类算法。

正则化

1. 验证集

   为了防止过拟合，首先我们可以设置一个验证集，在验证集上调整叶子结点的个数 $J$ 和树的棵数 $M$。寻找最优 $M$ 的过程类似神经网络的早停策略。

1. Shrinkage

   对加入模型的树加入缩放因子 $v$ 进行缩放：

   $$f_m(x)=f_{m-1}(x)+v\times\sum_{j=1}^Jc_{mj}I\{x\in R_{mj}\}$$

   这里的缩放因子可以理解为学习率。那么类比学习率对模型的影响，我们不难得出这样的结论：相同的 $M$，较小的 $v$ 导致较大的训练风险；$v$ 较小，需要较大的 $M$ ；$v$ 和 $M$ 共同控制模型的泛化能力。在实际应用中通常采用使用较小的 $v$ ，通过早停选择 $M$。

1. Penalty

   一种“广泛撒网，重点选拔”的方法，即先训练出大量的树，然后选择好的树放到最终的模型中，

   $$\alpha^\ast(\lambda)=\arg \min_\alpha\left\{\sum_{i=1}^N\left(y_i-\sum_k\alpha_kT_k(x_i)\right)^2+\lambda J(\alpha)\right\}$$

   其中，正则化项 $J(\alpha)$ 可以选择岭回归或者Lasso。该方法即是在原本的策略中加一个罚项，从而训练并选出泛化能力更好的模型。

模型组合

得到一组分类器

交叉验证

同一个分类器在不同的数据集上训练。

Bagging

对训练集进行部分采样。主要目的是为了降低方差。

典型方法为随机森林。

组合分类器

Stacking

• 简单：只需利用已训练好的分类器，无需重新训练
• 分类器可以学习出分类器之间的相关性，比简单的朴素贝叶斯组合具有更好的效果
• 方便进行特征组合，比如：目标检测中的部位模型

Averaging

假设有 $h=1,\dots,H$ 个模型，每个模型的先验概率是 $p(h)$，则数据集的分布

$$p(X)=\sum_{h=1}^Hp(X|h)p(h)$$

解释：

• 只有一个模型能很好地解释数据；

• $p(h)$ 反映了各个模型之间的不确定性；

• 数据集增大，不确定性减少，$p(X|h)$ 收敛到某个模型。

Boosting

综合转载以下文章：

• 概率图模型(中)
• 概率图模型(下)

概率图

定义

概率图模型就是把概率和模型用图连接起来 ，图是一种数据结构，用来表示变量概率依赖关系的理论，结合概率论与图论的知识，利用图来表示与模型有关的变量的联合概率分布，通过图形化概率分布有助于求解复杂的模型，图还可以帮助人们做推理，表示和推理是图的两大主要功能。

形式有三种：有向图、无向图和混合图。

有向图

有向图 又称 贝叶斯网络 又称 信念网络。

例：多项式回归

$$D=\{\langle x_n,t_n \rangle\}_{n=1}^t\\ t_n=w^Tx_n+\varepsilon,\varepsilon\sim N(0,1)$$

回归模型是有监督的。在上图中，$ε$ 是一个均值为0，方差为1的高斯分布，加上一个常数，所以 $t_n$ 同样满足高斯分布，其中 $w$ 是未知的，为了求解 $w$ 可以使用最大似然估计进行求解，$t_n$ 依赖于 $x_n$ 和 $w$ ，所以高斯分布可以写成如下形式：

$$p(t_n|x_n,w)=N(t_n|w^Tx_n,1)=\frac{1}{\sqrt{2\pi}}\exp\left\{-\frac{(t_n-w^Tx_n)^2}{2}\right\}$$

这样写的好处，w可以理解成在极大似然估计中一个固定的点，当然w也有可能是一个随机变量。为了更好的解析一个模型，通常可以使用三个方法：极大似然估计、最大后验概率估计（正则化）、贝叶斯学习（把参数当成是随机变量）。

每个 $t_n$ 都会和 $w$ 相关，$w$ 满足某一分布，相当于很多个高斯分布叠加在一起，可以表现成如下形式：

对输入变量预测，假设输入值 $x$，得到根据观测数据为条件的对应的 $t$ 的概率分布，对于该问题的概率图模型如下：

它们有一个共同的参数 $α$，$t_n$ 和 $x_n$ 有关，$x_n$ 是一个标量，即一个固定的值，同时还和方差有关，使用空心圆来表示随机变量，确定的参数由小的实心圆表示，可得公式：

$$p(\textbf t,\textbf w|\textbf x,\alpha,\sigma^2)=p(\textbf w|\alpha)\prod_{n=1}^Np(t_n|\textbf w,x_n,\sigma^2)$$

该模型的联合概率密度函数如下：

$$p(\hat t,\textbf t,\textbf w|\hat x,\textbf x,\alpha,\sigma^2)=\left[\prod_{n=1}^Np(t_n|x_n,\textbf w,\sigma^2)\right]p(\textbf w|\alpha)p(\hat t|\hat x,\textbf w,\sigma^2)$$

例：Logistic Sigmoid函数

$x_1,x_2,\dots,x_M$ 都分别表示一个维度，每个维度之间并不相互独立。有：

$$p(y=1|x_1,\dots,x_M)=\sigma\left(w_0+\sum_{i=1}^Mw_ix_i\right)=\sigma(\textbf w^T\textbf x)$$

其中：

$$\sigma(\textbf w^T\textbf x)=\frac{1}{1+\exp(-\textbf w^T \textbf x)}$$

例：多元高斯模型

结点 $i$ 表示服从高斯分布的一元连续随机变量 $x_i$ ，该分布的均值是结点 $i$ 的父节点 $x_i$ 的状态的线性组合：

$$p(x_i|\text{pa}_i)=\cal N\left(x_i\left|\sum_{j \in \text{pa}_i}w_{ij}x_j+b_i,v_i\right.\right)$$

其中 $w_{ij}$ 和 $b_i$ 是控制均值的参数，$v_i$ 是 $x_i$ 的条件概率分布的方差。联合概率分布的对数为图模型中所有结点的条件分布的乘积的对数：

$$\begin{align*} \ln p(\textbf x)&=\sum_{i=1}^D\ln p(x_i|\text{pa}_i)\\ &=-\sum_{i=1}^D\frac{1}{2v_i}\left(x_i-\sum_{j\in\text{pa}_i}w_{ij}x_j-b_i\right)^2+\text{const} \end{align*}$$

随机变量 $x_i$ 的期望：

$$\mathbb E[x_i]=\sum_{j\in\text{pa}_i}w_{ij}\mathbb E[x_j]+b_i$$

它们的均值如下：

$$\pmb \mu=(b_1,b_2+w_{21}b_1,b_3+w_{32}b_2+w_{32}w_{21}b_1)^T$$

用一个图来表示三个高斯变量：

上述三个变量的协方差如下：

$$\pmb \Sigma=\left(\begin{aligned}&v_1&&w_{21}v_1&&w_{32}w_{21}v_1\\&w_{21}v_1 &&v_2+w_{21}^2v_1 &&w_{32}(v_2+w_{21}^2v_1)\\ &w_{32}w_{21}v_1 &&w_{32}(v_2+w_{21}^2v_1)&&v_3+w_{32}^2(v_2+w_{21}^2v_1) \end{aligned}\right)$$

例：Wet

一天早晨，Tracey离开她的房子，意识到她的草是湿的。是由 于连夜下雨，还是她昨晚忘了关掉洒水器?接着她注意到她邻居杰克的草也湿了。这就在一定程度上解释了洒水喷头已经关掉的可能性，因此她得出结论，可能昨晚下雨。 该模型的变量有 $R$、$S$、$J$、$T$，各自取值情况为{0,1}：

R=1：昨晚下雨了；	R=0：昨晚没有下雨；
S=1：水洒关掉了；	S=0：水洒没有关掉；
J=1：Jack家的草坪是湿的；	J=0：Jack家草坪没有湿；
T=1:Tracey家草坪是湿的；	T=0:Tracey家草坪没有湿；

保持该模型的联合概率分布一般性，得到分解式：

$$p(R,S,J,T)=\underbrace{p(T|S,R,J)}_{2^3=8}\underbrace{p(J|S,R)}_{2^2=4}\underbrace{p(R|S)}_2\underbrace{p(S)}_1$$

通过条件独立性，减少了对模型的变量的指定数目，使用概率图模型表示如下：

同时给出如下条件概率表：

根据条件概率表，以及相应的概率模型图，可以求解：

同理：

例：混合高斯模型

$$p(\textbf z)=\prod_{k=1}^K\pi_k^{z_k}$$ $$p(\textbf x|z_k=1)=\cal N(\textbf x|\pmb \mu_k,\pmb\Sigma_k)$$ $$p(\textbf x|\textbf z)=\prod_{k=1}^K\cal N(\textbf x|\pmb \mu_k,\pmb \Sigma_k)^{z_k}$$ $$p(\textbf x)=\sum_{\textbf z}p(\textbf z)p(\textbf x|\textbf z)=\sum_{k=1}^K\pi_k\cal N(\textbf x|\pmb\mu_k,\pmb\Sigma_k)$$

例：隐马尔可夫模型

把很多个时刻都关联起来，并且各个时刻之间还有相互关联，那么就变成了隐马尔科夫模型。

• 状态序列：$Z_t\in Q=\{q_1,q_2,\dots,q_N\}$

• 观测序列：$X_t\in V=\{v_1,v_2,\dots,v_M\}$

• $A$ —状态转移概率矩阵：$A=[a_{ij}]_{N\times N}, a_{ij}=P(Z_{t+1}=q_j|Z_t=q_i)$

• $B$—观测概率矩阵：$B=[b_j(k)]_{N\times M},b_j(k)=P(X_t=v_k|Z_t=q_j)$
• $\Pi$—初始状态概率矩阵：$\Pi=(\pi_i)_{1\times N},\pi_i=P(Z_1=q_i)$

这个模型的主要功能是进行推理，求其中的隐变量，在语音识别等领域用的比较多。

例：朴素贝叶斯

朴素贝叶斯假设条件独立性：

$$p(X_1,\dots,X_n|Y)=\prod_i p(X_i|Y)$$

推测接下来的概率：

$$p(Y=y_k|X_1,\dots,X_n)=\frac{p(Y=y_k)\prod_i p(X_i|Y=y_k)}{\sum_j p(Y=y_j)\prod_i p(X_i|Y=y_j)}$$

朴素贝叶斯主要的应用是垃圾邮件的分类等。

例：贝叶斯网络

不同的概率图可以将联合概率分布展开成不同的形式：

1. 

   $$p(x_1,x_2,x_3,x_4)=p(x_1|x_2,x_3,x_4)p(x_2|x_3,x_4)p(x_3|x_4)p(x_4)$$

2. 

   $$p(x_1,x_2,x_3,x_4)=p(x_3|x_1,x_2,x_4)p(x_4|x_1,x_2)p(x_1|x_2)p(x_2)$$

现在有概率图：

联合概率分布展开为：

$$p(x_1,x_2,x_3,x_4)=p(x_1|x_4)p(x_2|x_3,x_4)p(x_3)p(x_4)$$

证明在给定 $x_4$ 的条件下，考虑 $x_1$ 和 $x_2$ 是否独立：$p(x_1,x_2|x_4)=p(x_1|x_4)p(x_2|x_4)$。

$$\begin{align*} p(x_1,x_2|x_4)&=\frac{1}{p(x_4)}\sum_{x_3}p(x_1,x_2,x_3,x_4)\\ &=\frac{1}{p(x_4)}\sum_{x_3}p(x_1|x_4)p(x_2|x_3,x_4)p(x_3)p(x_4)\\ &=p(x_1|x_4)\sum_{x_3}p(x_2|x_3,x_4)p(x_3) \end{align*}$$

现在我们要将 $p(x_2|x_4)$ 转化成 $\sum_{x_3}p(x_2|x_3,x_4)p(x_3)$，得到：

$$\begin{align*} p(x_2|x_4)&=\frac{1}{p(x_4)}\sum_{x_1,x_3}p(x_1,x_2,x_3,x_4)\\ &=\frac{1}{p(x_4)}\sum_{x_1,x_3}p(x_1|x_4)p(x_2|x_3,x_4)p(x_3)p(x_4)\\ &=\sum_{x_3}p(x_2|x_3,x_4)p(x_3) \end{align*}$$

分析可知：$p(x_1,x_2|x_4)=p(x_1|x_4)p(x_2|x_4)$，所以给定 $x_4$ 的情况下 $x_1$ 和 $x_2$ 独立。

冲突变量

定义：给定一个路径，路径上的结点 $x$、$y$ 同时指向邻居结点 $z$，则 $z$ 结点为冲突变量。

上图中，(c)中给定 $z$ 时，$x$ 与 $y$ 不独立；(d)中给定 $z$ 时，$x$ 与 $y$ 不独立。

变化的规则：

1. 冲突且给定，保留；
2. 不冲突且不给定，保留；
3. 冲突且不给定，删除；
4. 不冲突且给定，删除。

几个例子：

无向图

定义

无向图模型（undirected graphical model），即马尔科夫图模型（Markov network），也称为马尔可夫随机场（Markov random field），含一组结点，每个结点都对应着一个变量或一组变量，链接是无向的，即不含有箭头。模型例子如下图所示，

形式上的定义为：它是由变量 $x_1,x_2 \dots,x_n$ 组成的概率分布 $p$，由 无向图 $G$ 定义，其中每个结点都对应着一个变量或一组变量，联合分布为

$$p(x_1,x_2 \dots,x_n) = \frac{1}{Z} \prod_{c \in C} \phi(\cal X_c)$$

其中 $c$ 代表 无向图 $G$ 中的 团块，$\cal X_c$ 代表团块中的变量，乘积包括了图中的所有团， 而团块，通俗来讲，它是无向图中的完全子图的集合，$\phi$ 代表势函数，势函数是一个关于参数的非负函数，参数的多少并不影响其非负的性质，比如 $\phi(x_1)$ 是关于 $x_1$ 的非负函数，$\phi(x_1,x_2 \dots,x_n)$ 同样也是一个非负函数，可以将势函数理解为输入变量的条件概率的分布的函数，对势函数求积分，理想状态为1，即 $\sum_{c \in C} \phi(x_c) =1$，而 $Z = \sum_{x_1, \dots, x_n} \prod_{c \in C} \phi_c(\cal X_c)$， 它是一个标准化常数，确保整个分布合为1。接下来举出一个简单的例子来更好的理解模型的定义。给定一个无向图 $G$ ，如下图。

它的团集合 $ C $ 是 $\{[a,c], [b,c]\}$，它的概率分布 $p(x_1,x_2 \dots,x_n) = \frac{1}{Z}\phi_{ac}(a,c) \phi_{bc}(b,c) $，$Z = \sum_{a,b,c} \phi_{ac}(a,c) \phi_{b,c}(b,c)$。

成对的马尔科夫网络：马尔可夫网络有一个特例，当每个团块只包含两个节点时，即每一个势函数中的变量只有两个变量，该情况下的马尔科夫网络叫做成对的马尔科夫网络（Pairwise Markov network）。

Gibbs分布：如果无向图模型能够表示成一系列在 $G$ 的最大团上的非负函数乘积的形式，这个无向图模型的概率分布 $p(X)$ 就称为Gibbs分布。

例： 玻尔兹曼机（Boltzmann Machine）

玻尔兹曼机（BM）是二值的马尔科夫随机场（Markov Random Filed），是一种基于能量的模型，其对应的联合概率分布为

$$p(\textbf x) = \frac{\exp(\frac{-E(\textbf x)}{T})}{Z}$$

玻尔兹曼机可以表示为带权重的无向图，

配分函数 $Z$ 通常难以计算，一般通过MCMC方法近似，生成一组服从 $p(\textbf x)$ 的样本，BM中的节点值是二值的（非0即1），且BM的节点是成对作用的，连接权重表示假设之间是相互支持还是互相抵消的关系，如果 $x_i$ 和 $x_j$ 的取值都为1，且 $w_{ij}>0$，则能量下降，概率变大，其能量函数为：

$$\begin{align*} E(\textbf x) &= E(\textbf X=\textbf x)\\ &= -(\sum_{i<j}w_{ij}x_ix_j + \sum_ib_ix_i)\\ \end{align*}$$

解释一下玻尔兹曼机基于能量模型的原因：首先，统计力学的结论表明，任何概率分布都可以转变成基于能量的模型，其次，对于一个给定的数据集，如果不知道其潜在的分布形式，那是非常难学习的，似然函数都写不出来。如果知道是高斯分布或者多项分布，那可以用最大化似然函数来学出需要学习的对应参数，但是如果分布的可能形式都不知道，这个方法就行不通。最后，能量函数能为无监督学习方法提供目标函数与目标解。

例：限制玻尔兹曼机（Restricted Boltzmann Machines）

限制玻尔兹曼机（RBM）包含两个层，分别为可见层（visible layer）和隐藏层（hidden layer）。神经元之间的连接具有如下特点：层内无连接，层间全连接，RBM对应的图是一个二部图。一般来说，可见层单元用来描述观察数据的一个方面或一个特征，而隐藏层单元的意义一般来说并不明确，可以看作特征提取层。RBM和BM的不同之处在于，BM允许层内神经元之间有连接，而RBM则要求层内神经元之间没有连接，它名字中的限制也就是指限制同一层内相互没有连接。

因此RBM的性质为：当给定可见层神经元的状态时，各隐藏层神经元的激活条件独立；反之当给定隐藏层神经元的状态是，可见层神经元的激活也条件独立。

其对应的联合概率分布为：

$$p(\textbf v,\textbf h) = \frac{\exp(-E(\textbf v,\textbf h))}{Z}$$

其能量函数为：

$$\begin{align*} E(\textbf v,\textbf h) &= -\sum_i a_iv_i - \sum_jb_jh_j - \sum_i\sum_jv_iw_{ij}h_j \\ &=-\textbf a^{\text{T}} \textbf v - \textbf b^{\text{T}}\textbf h - \textbf v^{\text{T}} W\textbf h \end{align*}$$

$Z$ 为配分函数，也称归一化因子

$$Z = \sum_{v,h}\exp(- E(v,h))$$

在RBM中，隐层神经元 $h_j$ 被激活的概率为

$$p(h_j|\textbf v) = \sigma(b_j + \sum_iW_{i,j}v_i)$$

由于是双向连接，显层的神经元同样能通过隐层神经元激活

$$p(v_i|\textbf h) = \sigma(c_i + \sum_jW_{i.j}h_j)$$

同一层神经元中具有独立性，得到

$$p(\textbf h|\textbf v) = \prod^{N_h}_{j=1} p(h_j|v)\\ p(\textbf v|\textbf h) = \prod^{N_h}_{i=1}p(v_i|h)$$

独立性推导过见受限玻尔兹曼机推导。

例：深度置信网络（Deep Belief Network）

深度信念网络，也称DBN，是神经网络的一种。既可以用于非监督学习，类似于一个自编码机；也可以用于监督学习，作为分类器来使用。作为神经网络，神经元自然是其必不可少的组成部分。DBN由若干层神经元构成，组成元件是限制玻尔兹曼机（RBM）。训练过程中，需要充分训练完上一层的RBM 之后才能训练当层的RBM，直到最后一层。

具体训练过程如下：分别单独无监督地训练每一层RBM网络，确保特征向量映射到不同特征空间中，尽可能保留特征信息，这个过程中，数据输入到可见层中，生成向量 $v$ 再通过权重 $W$ 传给 隐层 $h$ ，通过隐层激活单元与可是层输入之间的相关性差别作为权重更新的主要依据。

例：条件随机场（Conditional Random Field）

条件随机场是条件概率分布模型 $ p(Y|X)$， 表示的是给定一组输入随机变量X条件下，另一组输出随机变量 Y的马尔可夫随机场，CRF的特点是假设输出随机变量构成马尔可夫随机场，下图所示的是线性链条件随机场的模型

设 $X=(x_1,\dots x_n), Y = (y_1, \dots, y_n)$ 均为线性链表示的随机变量序列，在给定随机变量序列 $X$ 的条件下，随机变量序列 $Y$ 的条件概率分布 $p(Y|X)$ 构成条件随机场，同时满足：

$$p(y_i|X,y_1,\dots,y_{i-1}, y_{i+1},\dots, y_n) = p(y_i|X,y_{i-1},y_{i+1})$$

则称 $p(Y|X)$ 为线性链条件随机场。

特性

条件独立性

条件独立性

如果变量 $x, y$ 由未观测变量的路径连接，那么它们是相关的。但如果 $x$ 的邻居都被观测到，那么 $x$ 与其他变量无关。举个例子，下图中，在给定条件 $C$ 的情况下，$A, B$ 无关。

可得公式：

$$p(a,b|c) = p(a|c)p(b|c)$$

这与有向图模型中的 $A\rightarrow C\leftarrow B$ 相反。在有向图模型中，给定条件 $C$ 并且 $C$ 是一个冲突变量的情况下，$A$ 和 $B$ 不独立。

以点集为对象的马尔可夫性质（Global Markov Property）

给定三个互斥的变量集 $(A,B,S)$ ，如果 $S$ 将 $A$ 与 $B$ 分开，那么在给定条件 $S$ 的情况下，$A$ 与 $B$ 独立。

换个说法，所有由 $A$ 点集中的任何一个点和 $B$ 点集中的任何一个点之间的路径都需要经过 $S$ 点集中的一个点，则称 $S$ 隔离了 $A$ 和 $B$，即给定条件 $S$ 的情况下，$A$ 和 $B$ 是独立的。

以点为对象的马尔可夫性质（Local Markov Property）

在马尔可夫网络中，对于给定邻居结点的条件下，其条件概率的值等于给定除 $x$ 外所有结点条件下的条件概率的值， 即公式如下

$$p(x|\cal X \setminus x) = p(x |\text{ne} (x))$$

定义如下：$v$ 为无向图 $G$上的一个结点，$W$ 是与 $v$ 相连的所有结点，$O$ 为除 $W, v$ 以外的结点，则有：

$$p(v,O|W) = p(v|W)p(O|W)$$

给出一个小例子：

以节点4的所有邻居 $\{2,3,5,6\}$ 为条件，节点4与节点1，节点7相互独立。也就是说 $p(x_4| x_2,x_3,x_5,x_6)=p(x_4|x_1,x_2,x_3,x_5,x_6,x_7)$。

以对为对象的马尔可夫网络性质（Pairwise Markov Property）

设 $x,y$ 为无向图 $G$ 上任意两个不相连的结点，其他所有的结点为 $X \setminus \{x,y\}$，则有以下公式

$$p(x,y| X \setminus \{x,y\}) = P(x|X \setminus \{x,y\})p(y|X \setminus \{x,y\})$$

对于非负的势函数，局部、成对、以及全局的马尔可夫性质都是等价的。

马尔可夫随机场与贝叶斯网络对比

马尔可夫随机场与贝叶斯网络相比有几个优点：

• 它们可以应用于更广泛的问题，其中没有和变量依赖相关的天然方向性。

• 无向图可以简洁地表达某些依赖关系，贝叶斯网络不容易描述它们（尽管反过来也是如此）。

它们也有几个重要的缺点：

• 计算归一化常数 $Z$ 需要对可能为指数数量的赋值进行求和。 在一般情况下，我们会看到这将是NP难的；因此许多无向模型是棘手的并需要近似技术。
• 马尔可夫随机场可能难以解释。
• 从贝叶斯网络生成数据要容易得多，这在某些应用中很重要。

不难看出，贝叶斯网络是马尔可夫随机场的特例，带有非常特定类型的团聚因子（对应于条件概率分布并且暗示图中的有向非循环结构）和归一化常数 1。 特别是，如果我们取一个有向图 $G$，并向给定节点的所有父节点添加侧边（并消除它们的方向性），那么贝叶斯网络（看做变量及其祖先的因子）可以在生成的无向图上因式分解。 由此产生的过程被称为规范化（Moralization）。

Hammersley-clifford Theorem

Hammersley-Clifford定理指出：马尔科夫随机场和Gibbs分布是一致的。回顾一下马尔可夫随机场与Gibbs分布的概念：

• 马尔可夫随机场：对于一个无向图模型 $G$，对于其中的任意节点 $x_ i$，以除了它以外的所有点为条件的条件概率和以它的邻居节点为条件的条件概率相等，也就是满足以下公式，那么这个无向图就是马尔可夫随机场。

  $$p(x_i|X\setminus x_i) = p(x_i|\text{ne}(x_i))$$

• Gibbs分布：如果无向图模型能够表示成一系列在G的最大团（们）上的非负函数乘积的形式，也就是满足以下公式，这个无向图模型的概率分布 $ p(X)$ 就称为Gibbs分布。

  $$p(X) = \frac{1}{Z} \prod_{c \in C} \phi(\cal X_c)$$

以下例子，将马尔可夫随机场表示为多个势函数相乘的函数

从结点1开始：

$$p(x_1,\dots x_7) = p(x_1|x_2,x_3)p(x_2,x_3,x_4,x_5,x_6,x_7)$$

根据 $p(x| \cal X \setminus x) = p(x|\text{ne}(x))$可得

$$\begin{align*} p(x_1,x_2,x_3|x_4,x_5,x_6,x_7) &=p(x_1,x_2,x_3|x_4)\\ \Rightarrow \quad p(x_2,x_3|x_4,x_5,x_6,x_7) &= p(x_2,x_3|x_4) \end{align*}$$

可得：

$$\begin{align*} p(x_1,\dots x_7) &= p(x_1|x_2,x_3)p(x_2,x_3|x_4)p(x_4,x_5,x_6,x_7)\\ &= p(x_1|x_2,x_3)p(x_2,x_3|x_4)p(x_4|x_5,x_6)p(x_5,x_6| x_5)p(x_7)\\ \end{align*}$$

即得到的势函数（每个最大团块区域内变量）相乘的函数形式为：

$$p(x_1,\dots x_7) = \frac{1}{Z}\phi(x_1,x_2,x_3)\phi(x_2,x_3,x_4)\phi(x_4,x_5,x_6)\phi(x_5,x_6,x_7)$$

马尔可夫网络的分布是可表示的，它是在图中团块上定义的势函数的乘积，简单来说，马尔可夫网络与Gibbs分布的相互转换可以按以下形式理解，$G$ 表示马尔可夫网络，$F$ 为Gibbs分布。

• $G \Rightarrow F$：$F$ 表示一个连乘的式子，这些连乘的因子是由在定义在 $G$ 中的最大团块的势函数；
• $F \Rightarrow G$：给定一个由势函数组合的式子，通过对该分解式可得到马尔科夫网络。

证明Hammersley-Clifford定理的过程比较繁琐，就不再展开，详情可见Hammersley-Clifford定理证明。

条件独立判断

Ancestral Graph：Identify the ancestors $A$ of the nodes $X \cup Y \cup Z$. Retain the nodes $X \cup Y \cup Z$ but remove all other nodes which are not in $A$, together with any edges in or out of such nodes

Moralisation：Add link between any two remaining nodes which have a common child, but are not already connected by an arrow. Then remove remaining arrow heads.

Separation：Remove links neighbouring $Z$. In the undirected graph so constructed, look for a path which joins a node in $X$ to one in $Y$.

举个例子：

根据给定模型图，判断给定条件 $i,d$ 情况下，$a$ 和 $b$ 是否独立：

这个例子中，$X$ 为结点 $a$，$Y$ 为结点 $b$，$Z$ 为结点 $i, d$。(a)图像给定一个初始的BN图模型。首先根据Ancestral Graph，去掉无关结点，得到图(b)，然后根据Moralisation的规则去掉 $e$ 和 $f$ 分别到 $i$ 的有向连接，将 $e$ 和 $f$ 用无向线连接起来得到无向图(c)，最后通过图(c)，可以观察到在给定 $i,d$ 的条件下，$a$ 和 $b$ 是相互独立的。

Lattice Model（格子图模型）

设计一个模型，模型中的每个变量取值只有两种 $\{-1,1\}$。一共有9个变量 {$x_1,\dots, x_9$}，均匀分布的在模型中，相邻的变量是同样的状态，

$$p(x_1,x_2 \dots,x_n) = \frac{1}{Z} \prod_{i,j} \phi(x_i,x_j)$$

例：Ising Models（伊辛模型）

Ising模型根据定义势函数，其图模型同样是Lattice图模型，势函数定义如下：

$$\phi_{i,j}(x_i,x_j) = \exp(-\frac{1}{2T}(x_i-x_j)^2)$$

Ising模型是对应着物理学中的铁磁性物体的相变性质得到的，许多小铁磁性物质排成一列组成一维的Ising模型，随着温度增大，小磁性的物质之间的磁力就越小，此时没有出现全局磁化现象，比较该磁场相对不稳定；当T值减小则，小磁性的物质之间的磁力就越大，此时全部的小磁性物质稳定排列，出现全局磁化的现象。产生更大的磁力。考虑到二维的情况下，即Lattice模型下，将给定一个 $T$ 值，居里温度设为 $T(c)\sim2.269$，该温度下，Ising模型中的变量之间相互依赖，排列整齐，当温度增大或减少在一定程度上破坏了初始的稳定状态，假设变量 $M = |\sum^N_{i=1}x_i|/N$，该变量表明变量的平均整齐排列的值，通过改变 $T$ 值，$M$ 的曲线图像如下：

通过改变局部变量之间的关系从而对全部的变量产生影响。

例：图像去噪

上图中，左图为原始图片，中间是受到噪声污染的图像，我们的目标是从噪声污染的图片中恢复图像，如右图所示。

将图像的像素值标为 $\{+1,-1\}$，原始图像的第 $i$ 个像素记为 $x_i$，受到污染的噪声图像的第 $i$ 个像素记为 $y_i$，由于噪声等级比较小，因此 $x_i $ 和 $y_i$ 之间有着强烈的相关性。我们还知道图像中相邻像素 $x _i $ 和 $x_ j $ 的相关性很强。这种先验知识可以使用马尔可夫随机场模型进行描述，它的无向图如下图所示。

这个图中有两种类型的团块，每一种团块包含两个变量。形如 $\{x_i, y_i\}$ 的团块有一个关联的能量函数，表达了这些变量之间的相关性。对于这些团块，我们选择一个非常简单的能量函数 $−\beta x_iy_i$，其中 $\beta $ 是一个正的常数。这个能量函数的效果是：当 $x_i$ 和 $y_i$ 符号相同时，能量函数会给出一个较低的能量（即较高的概率），而当 $x_i$ 和 $y_i$ 符号相反时，能量函数会给出 一个较高的能量。

剩余的团块由变量 {$x_i,x_j$} 组成，其中 $i$ 和 $j$ 是相邻像素的下标。与之前一样，我们希望当两个像素符号相同时能量较低，当两个像素符号相反时能量较高，因此我们选择能量函数 $−\alpha x_i x_j$，其中 $\alpha$ 是一个正的常数。

由于势函数是最大团块上的一个任意的非负的函数，因此我们可以将势函数与团块的子集上的任意非负函数相乘，或者等价地，我们可以加上对应的能量。

模型的完整的能量函数的形式为：

$$E(\cal X,\cal Y) = - \alpha \sum_{i,j}x_ix_j - \beta \sum_ix_iy_i$$

联合概率分布形式为：

$$\begin{align*} p(\cal X,\cal Y) &= \frac{1}{Z}\exp(-E(\cal X,\cal Y)) \\ &= \frac{1}{Z}\exp(\alpha\sum_{i,j}x_ix_j + \beta\sum_ix_iy_i) \end{align*}$$

去噪的过程可以解释为给定被噪声污染的图像 $\cal Y$，求解隐变量 $\cal X $ ，这与MAP相似，于是有

$$\begin{align*} \cal X^* &= \underset{\cal X}{ \text{argmax}}\, p(\cal X|\cal Y) \\ &= \underset{\cal X}{\text{argmax}} \, \ln p(\cal X,\cal Y)\\ &= \underset{\cal X}{\text{argmax}} \, (\alpha\sum_{i,j}x_ix_j + \beta \sum_ix_iy_i) \end{align*}$$

以两个像素点为例，假设估计值 $\alpha$ =8, $\beta = 10$，观测值 $y_1=+1,y_2=-1$，则当 $x_1=+1,x_2=+1$ 时，$ \ln p(X,Y) = \alpha\sum_{i,j}x_ix_j + \beta\sum_ix_iy_i = 8$， 类似计算，当 $x_1=+1,x_2=-1$ 时，$ \ln p(X,Y) =12$， 当$x_1=-1,x_2=+1$ 时，$ \ln p(X,Y) =-28$，当 $x_1=-1,x_2=-1$ 时，$ \ln p(X,Y) =8$，可知，原始图像的像素值最可能为 $(+1,-1)$。

推理

基于概率图模型定义的联合概率分布，我们能够对目标变量的边缘分布或条件分布进行计算，这一过程被成为推理。

对于概率图模型，还需要确定具体分布的参数，通常使用极大似然估计或最大后验概率估计求解。但如果将参数视为待推理的变量，则参数估计和推理十分相似，可以“吸收”到推理问题中。

概率图模型的推理方法可分为两类，一类是精确推理，能够计算出目标变量的边缘分布或条件分布的精确值；另一类是近似推理，可以在较低的时间复杂度下得到近似解。

精确推理

因子图

引用wiki里对因子图的定义，将一个具有多个变量的全局函数因子分解，得到局部函数乘积，由此得到的双向图就是因子图（Factor Graph）。可以将因子图表示为三元组 $ G = (X,F,E)$， $X = x_1,\dots, x_n$ 为变量结点，$F = f_1,\dots, f_n$ 为因子结点， $E$ 为边的集合。设有一个全局函数 $p(x_1,\dots,x_n)$，分解完毕之后就可以将其表示为如下公式

$$p(x_1,\dots,x_n) = \prod^m_{j=1} f_j(S_j)$$

其中 $S_j \sube X $，其对应的因子图 $G =(X,F,E)$，因子节点 $F=\{f_{1},f_{2},\dots ,f_{n}\}$。因子图用圆圈表示概率分布中的每个变量，用正方形表示联合概率分布中的每个因子$f_s(S_j)$，各函数表述的是变量之间的关系，比如在有向图中，就是条件概率，而在马尔可夫随机场中就是势函数。在每个因子结点和因子所依赖的变量结点之间引入无向链接。简单理解一下因子图就是在有向图、无向图的基础上，引入了额外的结点来表示函数本身。下图就是一个因子图。

当联合概率分布是由无向图表示时，我们可以直接由无向图推出因子图。因子图的变量结点与无向图对应，然后构造额外的因子结点，对应于最大团，因子与最大团的势函数相等，如图所示。注意，不同的因子图可能对应相同的无向图。

类似的，对于有向图，先构造对应的变量结点，然后构造对应于条件分布的因子结点，最后添加上合适的连接，如下所示

以上三个都是 $p(x_1)p(x_2)p(x_3|x_1,x_2)$。

而同一个因子图可以转换为有向图也可以转换为无向图，如图所示

图(a)是因子图，而图(b)图(c)分别是无向图与有向图。

sum-product算法

在概率图中，求某个变量的边缘分布是常见的问题。这问题有很多求解方法，其中之一就是把贝叶斯网络或马尔科夫随机场转换成因子图，然后用sum-product算法求解。换言之，sum-product算法就是一个计算边缘概率的低复杂度的算法，求边缘分布需要大量的累加运算，而sum-product算法将大量的累加运算分配到乘积项内，从而降低复杂度。

sum-product算法可以使用信息的传递进行表示，我们将变量节点 $x$向因子结点$f$传递消息的过程记为 $u_{x \Rightarrow f}(x)$，因子结点$f$向变量结点 $x$ 传递消息记为 $u_{f \Rightarrow x}(x)$，变量的边缘分布于它接受的所有消息的乘积成正比，即

$$p(x) \propto \prod_{f \in ns(x)}u_{f \Rightarrow x}(x)$$

于是我们只要选择一个根结点，从所有叶子结点向根结点传递一次信息，根结点从叶子结点传递一次信息，就可以计算出所有变量的边缘分布。信息的传递分为四种情况：

1. 叶子结点为因子结点，传递的信息为：

   

2. 叶子结点为变量结点，传递的信息为：

   

3. 非叶子变量结点传递信息给因子结点时，传递的信息为：

   

4. 非叶子因子结点传递信息给变量结点时，传递的信息为 ：

   

   其中 $\chi_f$ 指的是 {$y_1,y_2,y_3, x$}， $ \sum_{\chi_f \setminus x}\phi_f(\chi_f)= \sum{y_1,y_2,y_3}f(y_1,y_2,y_3,x) $。

变量的边缘分布就可以通过这样的方式得到：

举一个例子，给定如下的因子图，现在欲求解 $x_2$ 的边缘分布。

本题中 $ x_2 $ 边缘分布的定义为联合概率分布对其他变量求和，即：

$$\tilde{p}(x_2) = \sum_{x_1}\sum_{x_3}\sum_{x_4} \tilde{p}(X)$$

此时，整个图未归一化的联合概率分布为：

$$\tilde{p}(X) = f_a(x_1,x_2)f_b(x_2,x_3)f_c(x_2,x_4)$$

使用sum-product算法对 $ x_ 2 $ 的边缘分布进行计算：

$$\begin{align*} \tilde{p}(x_2) &=u_{f_a \Rightarrow x_2}(x_2) * u_{f_b \Rightarrow x_2}(x_2) *u_{f_c \Rightarrow x_2}(x_2)\\ &=\sum_{x_1}f_a(x_1,x_2) \sum_{x_3}f_b(x_3,x_2) \sum_{x_4}f_a(x_4,x_2)\\ &= \sum_{x_1}\sum_{x_3}\sum_{x_4} \tilde{p}(X) \end{align*}$$

这与边缘分布的定义一致。求的如果是 $x_3$ 的边缘分布，那么

$$\begin{align*} \tilde{p}(x_3) &=u_{f_b \Rightarrow x_3}(x_3) \\ &= \sum_{x_2}f_b(x_2,x_3) u_{x_2 \Rightarrow f_b} \\ &= \sum_{x_2}f_b(x_2,x_3) u_{a \Rightarrow x_2} u_{c \Rightarrow x_2} \\ &= \sum_{x_2}f_b(x_2,x_3) \sum_{x_1}f_a(x_1,x_2)\sum_{x_4}f_c(x_2,x_4) \end{align*}$$

归一化因子的计算方法如下：

$$Z = \sum_x \prod_{f\in \text{ne}(x)}u_{f \Rightarrow x}(x)$$

关于更多sum-product算法的例子可见sum-product例子。

由于 $\tilde p(X)$ 是未归一化的，上面得到的 $\tilde p(x_2)$ 也是未归一化的，归一化因子的计算方法如下：

$$Z=\sum_x\prod_{f\in\text{ne}(x)}\mu_{f\rightarrow x}(x)$$

如果图比较大，概率连乘可能会变得非常小，使得难以精确计算。一种处理方法是使用对数，

$$\lambda = \log \mu$$

相应的有：

$$\log p(x)=\sum_{f \in \text{ne}(x)}\lambda_{f\rightarrow x}(x)$$

非叶子变量结点向因子结点传递消息修改为：

$$\lambda_{x\rightarrow f}(x)=\sum_{g\in \{\text{ne}(x)\setminus f\}}\lambda_{g\rightarrow x}(x)$$

非叶子因子结点向变量结点传递消息修改为：

$$\lambda_{f\rightarrow x}(x)=\log\left(\sum_{\cal X_f\setminus x}\phi_f(\cal X_f)\exp\left(\sum_{y\in \{\text{ne}(f)\setminus x\}}\lambda_{y\rightarrow f}(y)\right)\right)$$

这里，对数结果的幂运算可能会发生数值上涨，解决方法是在指数中减去一个常数。我们令：

$$\lambda_{y\rightarrow f}^\ast=\max_{y\in\{ne(f)\setminus x\}}\lambda_{y\rightarrow f}(y)$$

则非叶子因子结点向变量结点传递消息修改为：

$$\lambda_{f\rightarrow x}(x)=\lambda_{y\rightarrow f}^\ast+\log\left(\sum_{\cal X_f\setminus x}\phi_f(\cal X_f)\exp\left(\sum_{y\in \{\text{ne}(f)\setminus x\}}\lambda_{y\rightarrow f}(y)-\lambda_{y\rightarrow f}^\ast\right)\right)$$

即可避免这个问题。

近似推理

见另外一篇博客。