机器学习--EM算法

引入

当我们假设了一个概率分布之后往往需要做的是进行参数的估计，比如我们假设了一个概率分布P(y)他表示随机变量y的概率分布，我们此时还有一些样本$y_1,y_2,y_3...y_n$我们需要对这个概率分布的参数$\theta$进行概率估计。
我们最常使用的概率估计是极大似然估计也叫做MLE，它的做法就是找到一个$\theta$让所有的样本的出现概率尽可能大，即：
$$\hat{\theta }=\underset{\theta }{\mathrm{argmax}}(\prod _{i=1}^{n}p(y_i))$$
为了方便表示究竟是那个参数下的概率分布，我们把当前的概率分布写作$p(y|\theta )$ 表示在参数$\theta$下的$y$的概率分布，那么上式可以写成
$$\hat{\theta }=\underset{\theta }{\mathrm{argmax}}(\prod _{i=1}^{n}p(y_i|\theta ))$$
在进一步，我们记$P(Y|\theta )=\prod _{i=1}^{n}p(y_i|\theta )$其中Y就表示了所有非样本。并且我们在进行求解上述式子时我们常常会加入一个log函数来简化运算。于是上式就变成了下面的式子
$$\hat{\theta }=\underset{\theta }{\mathrm{argmax}}log[P(Y|\theta )]$$
对于复杂的概率分布来说想要直接解出这个式子是十分困难的，可能甚至没有解析解，那么此时我们就需要寻找另一种方法来进行求解

含有因变量的参数估计

什么是含有因变量的参数估计呢，举一下李航老师书上的例子。
我们依次抛两次硬币，分别是A,B,C，我们先抛出A硬币后依据A硬币的结果来决定下一次是抛B还是C，最后我们求出一个第二次抛硬币的结果的概率分布P(y)这里y表示第二次抛硬币的结果，而我们知道第二次抛硬币的结果依赖于第一次抛硬币的结果，但是我们并不知道第第一次抛的是什么，第一次抛硬币结果的概率分布也是我们计算P(y)的一个关键。

可以直接把第一次的抛硬币结果也记做一个随机变量，叫做z，我们就得到一个新的概率分布P(y, z)表示在第一次结果为z第二次结果为y时的概率分布。根据边缘分布的计算式子，我们可以得到
$$p(y)=\sum _{z}p(y,z)$$
其中这个z就是因变量，简单的来说可以把z当成一个种子，对应的z会生成对应的y。

根据上式，我们知道极大似然估计可以写作
$$\hat{\theta }=\underset{\theta }{\mathrm{argmax}}log[P(Y|\theta )]$$
那么对于上述含有隐变量的概率分布，我们就可以强制的把因变量加上去。
$$\hat{\theta }=\underset{\theta }{\mathrm{argmax}}log\sum _{Z}P(Y,Z|\theta )$$
可以看到这里Z是大写的，也就是意味着隐变量可以不止一个。

然后根据概率的乘法公式$p(x,y)=p(x)p(y|x)$我们可以把上式进行改写
$$\hat{\theta }=\underset{\theta }{\mathrm{argmax}}log\sum _{Z}P(Y|Z,\theta )P(Z|\theta )$$
我们设
$$L(\theta )=log\sum _{Z}P(Y|Z,\theta )P(Z|\theta )$$
我们希望极大化这个函数，而EM算法给出了一个方法进行极大化，那就是通过不断地迭代求解$\theta$

首先记${\theta }^{(1)}$为我们初始的$\theta$，即我们初始化的参数，注意，一旦有了$\theta$什么都能计算出来了，下面有${\theta }^{(1)}$的式子其实都是常数，因为可以计算出来。

有了${\theta }^{(1)}$我们就可以计算一下我们需要极大化$L(\theta )$和我们现在的$L({\theta }^{(1)})$之间的差距。
$$L(\theta )-L({\theta }^{(1)})=log\sum _{Z}P(Y|Z,\theta )P(Z|\theta )-log\sum _{Z}P(Y|Z,{\theta }^{(1)})P(Z|{\theta }^{(1)})$$

然后我们使用一个trick，我们让第一个式子中同时乘以除以$P(Z|Y,{\theta }^{(1)})$，方便我们使用jensen不等式，同时我们把后一项写成$P(Y|{\theta }^{(1)})$。
$$L(\theta )-L({\theta }^{(1)})=log\sum _{Z}P(Z|Y,{\theta }^{(1)})\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(1)})}-logP(Y|{\theta }^{(1)})$$
然后我们合并这两项
$$L(\theta )-L({\theta }^{(1)})=log\sum _{Z}P(Z|Y,{\theta }^{(1)})\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(1)})P(Y|{\theta }^{(1)})}$$

Jensen不等式的式子是
$$log\sum _j{\lambda }_j{y}_j\ge \sum _j{\lambda }_{j}log(y_j)$$
其中$\lambda$是离散概率分布
我们就可以使用jensen不等式得到
$$\begin{gathered} L(\theta )-L({\theta }^{(1)})=log\sum _{Z}P(Y|Z,{\theta }^{(1)})\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Y|Z,{\theta }^{(1)})P(Y|{\theta }^{(1)})} \\ \ge \sum _{Z}P(Y|Z,{\theta }^{(1)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Y|Z,{\theta }^{(1)})P(Y|{\theta }^{(1)})} \end{gathered}$$
我们最后得到
$$L(\theta )-L({\theta }^{(1)})\ge \sum _{Z}P(Z|Y,{\theta }^{(1)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(1)})P(Y|{\theta }^{(1)})}$$
我们进行移项
$$L(\theta )\ge \sum _{Z}P(Y|Z,{\theta }^{(1)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Y|Z,{\theta }^{(1)})P(Y|{\theta }^{(1)})}+L({\theta }^{(1)})$$
于是我们神奇的得到了一个$L(\theta )$的下界
也就是如果我们想让$L(\theta )$增大，我们就可以让它的下界增大，因为$L(\theta )$大于等于它的下界。
我们记$B(\theta ,{\theta }^{(1)})=\sum _{Z}P(Z|Y,{\theta }^{(1)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(1)})P(Y|{\theta }^{(1)})}+L({\theta }^{(1)})$
那我们现在的任务就是最大化这个下界，假设我们最大化的结果是${\theta }^{(2)}$
那么就有
$${\theta }^{(2)}=\underset{\theta }{\mathrm{argmax}}B(\theta ,{\theta }^{(1)})$$
我们就可以由此得到${\theta }^{(2)}$，那么它就是最终的答案么，显然不是，我们可以把它当成${\theta }^{(1)}$对待然后求出${\theta }^{(3)}$可见我们可以得到一个一般的式子
$${\theta }^{(i+1)}=\underset{\theta }{\mathrm{argmax}}B(\theta ,{\theta }^{(i)})$$

我们再对最优化的式子进行一下化简，使我们更方便计算
$${\theta }^{(i+1)}=\underset{\theta }{\mathrm{argmax}}[P(Z|Y,{\theta }^{(i)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})}+L({\theta }^{(i)})]$$
上式我们最优化的$\theta$对于$L({\theta }^{(i)})$只是一个已知的常数，所以我们不需要对它进行考虑，于是得到
$${\theta }^{(i+1)}=\underset{\theta }{\mathrm{argmax}}[P(Z|Y,{\theta }^{(i)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})}]$$

同样，我们把对数除法拆开，也可以得到一个常数
$${\theta }^{(i+1)}=\underset{\theta }{\mathrm{argmax}}[P(Z|Y,{\theta }^{(i)})log[P(Y|Z,\theta )P(Z|\theta )]-P(Z|Y,{\theta }^{(i)})log[P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})]]$$
我们也可以把它省掉
$${\theta }^{(i+1)}=\underset{\theta }{\mathrm{argmax}}P(Z|Y,{\theta }^{(i)})log[P(Y|Z,\theta )P(Z|\theta )]$$
最后运用乘法公式进行合并我们得到了一个更简洁的式子
$${\theta }^{(i+1)}=\underset{\theta }{\mathrm{argmax}}P(Z|Y,{\theta }^{(i)})logP(Y,Z|\theta )$$
我们一般把后面这个式子叫做$Q(\theta ,{\theta }^{(i)})$ 也就是
$$Q(\theta ,{\theta }^{(i)})=P(Z|Y,{\theta }^{(i)})logP(Y,Z|\theta )$$

至此EM算法的推导也就结束了，关于EM算法的收敛性证明就不证了，其大概想法就是如果这个迭代让$P(Y|\theta )$不断增大，那么这个算法就收敛，因为概率不会大于1，单调递增必收敛。
还有广义EM算法也暂时先不写，因为还没看，等以后有时间看了再回来补上吧。

应用

隐马尔科夫参数估计

隐马尔科中，我们做出了假设，即每个状态由其对应的一个隐状态生成。
这里就可以看出，隐马尔科夫就是一个含有隐状态的模型。
那么按照EM算法，就有两步需要进行走。

E步：
求$Q(\lambda ,\hat{\lambda })$，其中$\lambda$是参数
$$\begin{gathered} Q(\lambda ,\hat{\lambda })=\sum _{I}logP(O,I|\lambda )P(I|O,\hat{\lambda }) \\ =\sum _{I}logP(O,I|\lambda )\frac{P(I,O|\hat{\lambda })}{P(O|\hat{\lambda })} \end{gathered}$$
上式中$I,Q$分别代表，隐状态序列和观测态序列。
设$\lambda =(\pi ,A,B)$，就有
$$P(O,I|\lambda )={\pi }_{I_1}{B}_{I_1,O_1}\prod _{t=2}^T[A_{I_{t-1},I_t}{B}_{I_t,O_t}]$$
带入第一个式子
$$\begin{gathered} Q(\lambda ,\hat{\lambda })=\sum _{I}logP(O,I|\lambda )\frac{P(I,O|\hat{\lambda })}{P(O|\hat{\lambda })} \\ =\sum _{I}log{\pi }_{I_1}{B}_{I_1,O_1}\prod _{t=2}^T[A_{I_{t-1},I_t}{B}_{I_t,O_t}]\frac{P(I,O|\hat{\lambda })}{P(O|\hat{\lambda })} \end{gathered}$$

为了方便求参数，我们把每个因变量按照$\pi ,A,B$的方式进行合并同类项
$$\begin{gathered} Q(\lambda ,\hat{\lambda })=\sum _I\frac{P(I,O|\hat{\lambda })}{P(O|\hat{\lambda })}log{\pi }_{I_1}+\sum _I\frac{P(I,O|\hat{\lambda })}{P(O|\hat{\lambda })}\sum _{t=2}^{T}log[A_{I_{t-1},I_t}] \\ +\sum _I\frac{P(I,O|\hat{\lambda })}{P(O|\hat{\lambda })}\sum _{t=1}^{T}log[B_{I_t,O_t}] \end{gathered}$$
E步就结束了
接下来是M步：
M我们需要极大化$Q(\lambda ,\hat{\lambda })$，由于分母上的$P(O|\hat{\lambda })$对最优化没有影响，所以我们舍去。
于是我们就需要极大化：
$$\begin{gathered} Q(\lambda ,\hat{\lambda })=\sum _{I}P(I,O|\hat{\lambda })log{\pi }_{I_1}+\sum _{I}P(I,O|\hat{\lambda })\sum _{t=2}^{T}log[A_{I_{t-1},I_t}] \\ +\sum _{I}P(I,O|\hat{\lambda })\sum _{t=1}^{T}log[B_{I_t,O_t}] \end{gathered}$$
我们可以使用坐标轴下降法，依次选定每一项然后固定其他项，进行极大化。

1. 选定$\pi$项，然后其他项与其都是加和关系所以求导之后会变成0，所以不用关。
   于是我们就有
   $$\underset{\pi }{\mathrm{argmax}}\sum _{I}P(I,O|\hat{\lambda })log{\pi }_{I_1}$$
   我们可以对其进行改写
   $$\underset{\pi }{\mathrm{argmax}}\sum _{k=1}^{N}P(I_1=k,O|\hat{\lambda })log{\pi }_k$$
   这里解释一下为什么这么写，我们看上面的式子${\pi }_{I_i}$对应$\pi$中的一个参数，但是这么表示不方便我们求导，因为我们不知道$I_i$到底是多少，所以我们把以$I$为循环的累加，改为以$k$为循环的累加，其中k表示$\pi$中的第几个参数。
   那这么改为什么正确呢？
   首先我们看${\pi }_{I_1}$，由于我们把循环变量换成了k，而一个k对应着多个$I_1$。
   所以就有
   $$log{\pi }_k\sum _{I_1=k}P(I,O|\hat{\lambda })=log{\pi }_{k}P(I_1=k,O|\hat{\lambda })$$
   解释完毕后，我们继续看。
   我们注意到有一个限制条件，那么就是${\sum }_{i=1}^N{\pi }_i=1$
   所以我们用拉格朗日乘子法来进行限制。
   于是我们的优化就变成了。
   $$\sum _{k=1}^{N}log{\pi }_{k}P(O,I_1=k|\lambda )+\gamma (\sum _{i=1}^N{\pi }_i-1)$$
   我们进行求导，然后令其等于零，就得到了
   $${\pi }_k=\frac{P(O,I_1=k|\hat{\lambda })}{P(O|\hat{\lambda })}$$
2. 选定A项
   $$\sum _{I}P(I,O|\hat{\lambda })\sum _{t=2}^{T}log[A_{I_{t-1},I_t}]$$
   同样道理我们可以进行改写，不过这里我么可以注意到，对于$A_{ij}$有两个变量i，j所以改写后会多出来两个循环。
   $$\sum _{I}P(I,O|\hat{\lambda })\sum _{t=2}^{T}log{A}_{I_{t-1},I_t}=\sum _{k=1}^N\sum _{l=1}^N\sum _{t=2}^{T}log{A}_{kl}P(O,I_{t-1}=k,I_t=l|\hat{\lambda })$$
   同样的，我们也需要对概率进行约束。
   ∑ k = 1 N { [ ∑ l = 1 N ∑ t = 2 T l o g A k l P ( O , I t − 1 = k , I t = l ∣ λ ^ ) ] + γ ( ∑ i = 1 N A k j − 1 ) } \sum_{k = 1}^N\{[\sum_{l = 1}^N\sum_{t=2}^TlogA_{kl}P(O,I_{t - 1}=k,I_{t}=l|\hat \lambda)] +\gamma(\sum_{i = 1}^NA_{kj} - 1)\} k=1∑N​{[l=1∑N​t=2∑T​logAkl​P(O,It−1​=k,It​=l∣λ^)]+γ(i=1∑N​Akj​−1)}
   最后求导然后也可以得出结果。
   对于B也是同样道理的，这里就不再多写了。