模式识别——EM算法

一、密度估计
  对于某些有一定意义的数据，其可能不服从于任何标准的概率分布，使用密度估计算法估计其概率密度。一般的，一个非标准的分布都可以使用多个高斯函数进行拟合，称为高斯混合模型【GMM】。
  考虑未知的随机变量$\mathbf{z}$，其可能是隐藏的，或者未知的。考虑${\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)}$的联合概率分布$$p({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)})=p({\mathbf{x}}^{(i)}|{\mathbf{z}}^{(i)})p({\mathbf{z}}^{(i)})$$且${\mathbf{z}}^{(i)}\sim B(k,\varphi )$代表k个高斯分布的概率，以及${\mathbf{x}}^{(i)}|{\mathbf{z}}^{(i)}=j\sim N({\mathbf{\mu }}_j,{\mathbf{\Sigma }}_j)$代表已知的第j个高斯分布的数据概率分布。
  如果${\mathbf{z}}^{(i)}$是已知的，则可以使用极大似然估计，形如$$lnL(\varphi ,\mathbf{\mu },\mathbf{\Sigma })=\sum _{i=1}^{m}logp({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\varphi ,\mathbf{\mu },\mathbf{\Sigma })$$可以得到ϕj=∑i=1m1{z(i)=j}/mμj=∑i=1m1{c(i)=j}x(i)/∑i=1m1{c(i)=j}\phi_j = \sum_{i=1}^m 1\{z^{(i)} = j\}/m \\ \bm\mu_j = \sum_{i = 1}^m 1\{c^{(i)} = j\}\bm{x}^{(i)}/ \sum_{i = 1}^m 1\{c^{(i)} = j\}ϕj​=i=1∑m​1{z(i)=j}/mμj​=i=1∑m​1{c(i)=j}x(i)/i=1∑m​1{c(i)=j}其称为高斯判别模型【GDA】。然而${\mathbf{z}}^{(i)}$是未知的，可以考虑尝试使用模型猜测${\mathbf{z}}^{(i)}$的值，使用极大似然拟合出更好的参数的值，再去猜测${\mathbf{z}}^{(i)}$的值，并进行迭代，该算法称为最大期望算法。

---

二、EM算法
  EM算法的步骤如下：
  （1）猜测未知的${\mathbf{z}}^{(i)}$的值，称为E步；
  （2）最大似然估计参数的值，称为M步。
  （3）迭代（1）-（2），直到收敛。

2.1 E步
  E步用于猜测未知的${\mathbf{z}}^{(i)}$的值，形如$$w_j^{(i)}=p({\mathbf{z}}^{(i)}=j|{\mathbf{x}}^{(i)};\varphi ,\mathbf{\mu },\mathbf{\Sigma })$$该步计算了${\mathbf{z}}^{(i)}=j$的概率，即$w_j^{(i)}$表示了${\mathbf{x}}^{(i)}$由第j个高斯分布生成占所有高斯分布生成的概率，根据贝叶斯公式$$p(B_i|A)=P(B_i)P(A|B_i)/\sum _{j=1}^{n}P(B_j)P(A|B_j)$$形如$${\hat{\mathbf{z}}}^{(j)}=p({\mathbf{x}}^{(i)}|{\mathbf{z}}^{(i)}=j)p({\mathbf{z}}^{(i)}=j)/\sum _{l=1}^{c}p({\mathbf{x}}^{(i)}|{\mathbf{z}}^{(i)}=l)p({\mathbf{z}}^{(i)}=l)$$
2.2 M步
  M步用于更新对参数$\varphi$的更新，形如$$\begin{gathered} {\varphi }_j=\sum _{i=1}^m{w}_j^{(i)}/m \\ {\mathbf{\mu }}_j=\sum _{i=1}^m{w}_j^{(i)}{\mathbf{x}}^{(i)}/\sum _{i=1}^m{\mathbf{x}}^{(i)} \\ {\mathbf{\Sigma }}_j=\sum _{i=1}^m{w}_j^{(i)}({\mathbf{x}}^{(i)}-{\mathbf{\mu }}_j)({\mathbf{x}}^{(i)}-{\mathbf{\mu }}_j{)}^T/\sum _{i=1}^m{\mathbf{x}}^{(i)} \end{gathered}$$  其与${\mathbf{z}}^{(i)}$已知情况的区别仅仅是第i个点不是由第j个高斯分布生成的指示函数，而是概率，这对应了${\mathbf{z}}^{(i)}$已知与${\mathbf{z}}^{(i)}$未知的情况。

---

三、EM算法的一般形式
  GMM只是EM算法的一个特例，接下来介绍EM算法的一般形式。

3.1 琴生不等式
  首先介绍琴生【Jensen】不等式。考虑凸函数$f$，取随机变量$\mathbf{x}$，那么有$$f(E[\mathbf{x}])\le E[f(\mathbf{x})]$$其直观理解为凸函数的割线在函数曲线上方。若$f$是严格凸的，即$d^{2}f/d{x}^2>0$，那么上述等号当且仅当$\mathbf{x}=E[\mathbf{x}]$时成立，即$\mathbf{x}$以1的概率取得某常数值。

3.2 最大似然下界理论
  考虑模型的概率分布$p(\mathbf{x},\mathbf{z};\mathbf{\theta })$，其中，仅有$\mathbf{x}$可以被观测，尝试最大似然$$\begin{array}{rl}lnL(\mathbf{\theta })& =\sum _{i=1}^{m}logp({\mathbf{x}}^{(i)};\mathbf{\theta })\\ & =\sum _{i=1}^{m}log\sum _{{\mathbf{z}}^{(i)}}p({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\mathbf{\theta })\end{array}$$EM算法就是极大似然估计的过程，但其并非直接求解最大似然导数，其迭代的对似然函数创造了一个下界并寻找下界的最大值，从而靠近似然函数的最大值。由于$\mathbf{z}$未被观测，因此似然函数的变量仅有$\mathbf{x}$，因此，有$$\begin{array}{rl}ma{x}_{\mathbf{\theta }}lnL(\mathbf{\theta })& =ma{x}_{\mathbf{\theta }}\sum _{i=1}^{m}logp({\mathbf{x}}^{(i)};\mathbf{\theta })\\ & =ma{x}_{\mathbf{\theta }}\sum _{i=1}^{m}log\sum _{{\mathbf{z}}^{(i)}}p({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\mathbf{\theta })\end{array}$$引入辅助函数$Q_i({\mathbf{z}}^{(i)})$，其是未知随机变量${\mathbf{z}}^{(i)}$的概率分布，有$$\begin{gathered} Q_i({\mathbf{z}}^{(i)})\ge 0 \\ \sum _{{\mathbf{z}}^{(i)}}{Q}_i({\mathbf{z}}^{(i)})=1 \end{gathered}$$再考虑数学期望公式$$EX=\sum _{k=1}^m{X}_k{p}_k$$那么有$$\begin{array}{rl}\sum _{i=1}^{m}log\sum _{{\mathbf{z}}^{(i)}}p({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\mathbf{\theta })& =\sum _{i=1}^{m}log\sum _{{\mathbf{z}}^{(i)}}{Q}_i({\mathbf{z}}^{(i)})p({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\mathbf{\theta })/Q_i({\mathbf{z}}^{(i)})\\ & =\sum _{i=1}^{m}logE[p({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\mathbf{\theta })/Q_i({\mathbf{z}}^{(i)})]\end{array}$$这里是对所有的${\mathbf{z}}^{(i)}$对其概率分布乘以其函数进行求和，即该函数的期望。
  对数函数$log(x)$是一个凹函数，那么根据琴生不等式，有$$log(EX)\ge E[log(X)]$$故$$\begin{array}{rl}\sum _{i=1}^{m}logE[p({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\mathbf{\theta })/Q_i({\mathbf{z}}^{(i)})]& \ge \sum _{i=1}^{m}E[logp({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\mathbf{\theta })/Q_i({\mathbf{z}}^{(i)})]\\ & =\sum _{i=1}^m\sum _{{\mathbf{z}}^{(i)}}{Q}_i({\mathbf{z}}^{(i)})log(p({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\mathbf{\theta })/Q_i({\mathbf{z}}^{(i)}))\end{array}$$即$$lnL(\mathbf{\theta })\ge \sum _{i=1}^m\sum _{{\mathbf{z}}^{(i)}}{Q}_i({\mathbf{z}}^{(i)})log(p({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\mathbf{\theta })/Q_i({\mathbf{z}}^{(i)}))$$为了取得下界的最优值，即等号成立，即需$$\sum _{i=1}^{m}logE[p({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\mathbf{\theta })/Q_i({\mathbf{z}}^{(i)})]=\sum _{i=1}^{m}E[logp({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\mathbf{\theta })/Q_i({\mathbf{z}}^{(i)})]$$又需要对任何数据均有等式成立，因此需要选择可能的分配$Q$来保证。根据等号成立的条件，需要$$p({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\mathbf{\theta })/Q_i({\mathbf{z}}^{(i)})=Const$$即对任何${\mathbf{z}}^{(i)}$都取得相同的值，即$$Q_i({\mathbf{z}}^{(i)})=\Theta (p({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\mathbf{\theta }))$$又$$\sum _{{\mathbf{z}}^{(i)}}{Q}_i({\mathbf{z}}^{(i)})=1$$因此$$\begin{array}{rl}{Q}_i({\mathbf{z}}^{(i)})& =p({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\mathbf{\theta })/\sum _{{\mathbf{z}}^{(i)}}p({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\mathbf{\theta })\\ & =p({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\mathbf{\theta })/p({\mathbf{x}}^{(i)};\mathbf{\theta })\\ & =p({\mathbf{z}}^{(i)}|{\mathbf{x}}^{(i)};\mathbf{\theta })\end{array}$$

3.3 EM算法理论
  在E步中猜测${\mathbf{z}}^{(i)}$的值，即是$$Q_i({\mathbf{z}}^{(i)})=p({\mathbf{z}}^{(i)}|{\mathbf{x}}^{(i)};\mathbf{\theta })$$这为当前的参数$\varphi$的似然函数确定了紧下界。
  在M步中更新$\varphi$的值，即是$$\varphi =argma{x}_{\varphi }\sum _{i=1}^m\sum _{{\mathbf{z}}^{(i)}}{Q}_i({\mathbf{z}}^{(i)})log(p({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\mathbf{\theta })/Q_i({\mathbf{z}}^{(i)}))$$这在当前的紧下界中找到似然函数下界中的最大值，从而更新参数$\varphi$。
  从凸优化的角度考虑，可以考虑$$J(\varphi ,Q)=\sum _{i=1}^m\sum _{{\mathbf{z}}^{(i)}}{Q}_i({\mathbf{z}}^{(i)})log(p({\mathbf{x}}^{(i)},{\mathbf{z}}^{(i)};\mathbf{\theta })/Q_i({\mathbf{z}}^{(i)}))$$其中，E步选择$Q$最优化$J$，M步选择$\varphi$最优化$J$，这是一种坐标上升法。