机器学习-EM算法

定义

按照[wiki百科的定义，EM算法被用于寻找，依赖于不可观察的隐形变量的模型中参数的最大似然估计。在统计计算中，最大期望（EM）算法是在概率模型中寻找参数最大似然估计或者最大后验估计的算法，其中概率模型依赖于无法观测的隐性变量。最大期望算法经常用在机器学习和计算机视觉的数据聚类（Data Clustering）领域。最大期望算法经过两个步骤交替进行计算，第一步是计算期望（E），利用对隐藏变量的现有估计值，计算其最大似然估计值；第二步是最大化（M），最大化在E步上求得的最大似然值来计算参数的值。M步上找到的参数估计值被用于下一个E步计算中，这个过程不断交替进行。

极大似然估计

首先是最大估计的原理，wiki定义如下，进行通俗解释如下

问题引入

给定了一个样本集D，现在需要进行估计这个样本集的数据对于一个特征的分布(比如估计一个学校的学生的身高分布)，一般情况下D是很大的，因此这时候可以随机抽取一定数量的样本(N)进行样本集数据分布的估计。通常先假设样本集$D$服从一个分布，比如服从高斯分布$N(\mu,\delta^{2})$（对于极大似然估计，一定是要先假设数据服从某一个分布，如果不知道数据的总体分布则无法使用极大似然估计），那么这个分布中的参数$\mu$和$\delta$就是需要进行估计并计算的参数。

为了得到总体的样本分布，我们独立的按照概率密度$p(x|\theta)$抽取了N个样本，并组成了样本集$X=x_{1},x_{2},…,x_{N}$，其中$x_{i}$表示抽取到的第$i$个人的特征值(例如身高值)，我们通过样本集$X$来估计总体的位置参数$\theta$，这里的概率密度$p(x| \theta)$服从提前假设的高斯分布$N(\mu,\delta^{2})$，其中未知参数为$\theta=[\mu,\delta]^{T}$，极大似然的任务也就是去估计这个参数$\theta$

参数估计

得到样本$X$的概率是多大？

由于抽取到的样本$X$，是独立抽取到的，也就是说抽取每个的概率为$p(x_{i}|\theta)$，由于样本之间相互独立，因此得到样本集$X$的联合概率为
$$
L(\theta)=L(x_{1},x_{2},…,x_{N};\theta)=\prod_{i=1}^{N}p(x_{i}|\theta),::::: \theta \in \Theta
$$
上面这个式子表明了在概率密度参数为$\theta$时，得到$X$这组样本集的概率，这是一个关于$\theta$的函数。

这个式子反应了在不同概率密度参数取值情况下，得到当前样本集的可能性，因此称作参数$\theta$相对于样本集$X$的似然函数，记作$L(\theta)$，对$L$取对数，可以将原来的连乘变作连加，称作对数似然函数，如下
$$
H_{\theta}=\ln L(\theta)=\ln \prod_{i=1}^{N}p(x_{i}|\theta)=\sum_{i=1}^{N} \ln p(x_{i}|\theta)
$$
- 为什么要使用对数似然函数
  
  之所以使用对数似然函数，是因为把连乘变为连加，求导更加方便；概率累计会产生非常小的数值，有可能会小于计算机可以识别的精度，取对数之后，非常小的数值可以变得相对较大（例如概率为$1e-50$，对其取对数之后可以得到为$-50$，共容易计算）
为什么随机抽样得到的是这个样本$X$

整个样本集$D$那么大，为什么随机抽取$N$个样本得到的是这个样本集$X$的这些样本，而不是其他的样本，是否可以认为抽取$N$个样本，着$N$个样本在整体的样本集$D$中出现的概率是最大的，也就是对应的似然函数$L(\theta)$取得极大值
$$
\hat {\theta}= argmax L(\theta)
$$
那么这个参数$\hat \theta$就是$\theta$的极大似然估计量，也就是要求的值
如何求解极大似然函数

对于似然函数$L(\theta)$对所有的参数求导，并让偏导数为$0$，如果有$n$个参数吗，那么就得到了$n$个方程组，这些方程组的解就是似然函数的极值点，从而得到$\hat \theta$

极大似然估计总结

实际上可以看做一个反推的过程，在我们不知道样本数据分布的前提下，首先根据我们的先验来假设一个样本分布，之后寻找使得这个结果出现的可能性最大的条件，并将满足这个极值的条件的相关参数作为整个样本的分布。已知某个样本符合某种概率分布，但是具体的参数我们是不知道的，通过若干次抽样及实验，利用得到的结果推出参数的大概值，这个参数就可以标书整体的样本分布。当一个参数可以使得这个样本出现的概率最大时，我们显然不会再去选择其他小概率的样本，所以就用这个参数当做整个样本的概率密度分布的参数。

求解的一般步骤

写出似然函数，首先当然要假设样本服从某个概率分布
求解对数似然函数
求导并使得导数为0，得到似然方程组
求解似然方程组，最终得到我们要估计的参数

极大似然估计应用

回归问题的极小化平方差
- 最小二乘法思想
  
  对于回归问题的优化目标一般是使用的MSE的方法，即估计值和真实值之间的差的平方，并通过梯度下降的方法不断迭代最终得到参数值，这也是常用的最小二乘估计思想
$$
J(\theta)=\sum {i=1}^{n}(h{\theta}(x_{i})-y_{i})^{2}
$$

其中$h_{\theta}$是预测函数，$x_{i}$为输入的样本，$y_{i}$为输入的样本对应的label，
- 极大似然估计的思想
  
  假设样本分布满足正太分布$X \thicksim N(\mu,\delta^{2})$，其概率密度函数为
  $$
  f(x)=\frac{1}{\delta \sqrt{2\pi}} \exp({-\frac{(x-\mu)^2}{2\delta ^{2}}})
  $$
  可以得到抽取每个样本的满足的概率密度分布为$y_{i} \in N(\theta x_{i},\delta^2)$，得到每每个样本的概率为
  $$
  p(y_{i}|x_{i};\theta)=\frac{1}{\delta \sqrt{2\pi}} \exp({-\frac{(y_{i}- \theta x_{i})^2}{2\delta ^{2}}})
  $$
  可以得到似然函数为
  $$
  \begin{align}
  L(\theta)&=\prod_{i=1}^{n}p(y_{i}|x_{i};\theta) \
  &=\prod_{i=1}^{n}\frac{1}{\delta \sqrt{2\pi}} \exp({-\frac{(y_{i}- \theta x_{i})^2}{2\delta ^{2}}})
  \end{align}
  $$
  对数似然函数为：
  $$
  \begin{align}
  H(\theta)=\log L(\theta) &= \log \prod_{i=1}^{n}\frac{1}{\delta \sqrt{2\pi}} \exp({-\frac{(y_{i}- \theta x_{i})^2}{2\delta ^{2}}}) \
  &=-\frac{1}{2 \delta^{2}} \sum_{i=1}^{n}(y_{i}-\theta x_{i})^{2} -n\ln (\delta \sqrt{2 \pi}) \ &\Rightarrow \sum_{i=1}^{n}(y_{i}-\theta x_{i})^{2}
  \end{align}
  $$
  可以得到最大化对数似然函数等价于最小化最小二乘法优化目标
  $$
  arg \max_{\theta}H(\theta)\Leftrightarrow arg \min _{\theta}J(\theta)
  $$
  虽然最后的表达式一致，但这只是巧合，两种优化方法的原理截然不同
分类问题极小化交叉熵损失函数

实际上似然函数的极大化就是交叉熵的本质

对于逻辑斯蒂模型，设$p(y=1|x)=p,:: p(y=0|x)=1-p$，由此可以得到$p(y|x,\theta)=p^{y}(1-p)^{1-y}$

对应的似然函数为
$$
L(\theta)=\prod_{i=1}^{n}p(y|x,\theta)=\prod_{i=1}^{n}p^{y}(1-p)^{1-y}
$$
对数似然函数为
$$
H(\theta)=\log L(\theta)=\sum_{i=1}^{n} y_{i} \log p +(1-y_{i}) \log (1-p)
$$
极大化对数似然函数就是对应的交叉熵优化目标

EM算法

对于极大似然估计，我们可以假设整个样本服从一个分布，但是现实情况下一个样本不一定会服从一个分布，比如一个学校的学生的身高，对于男生和女生是满足两个正太分布的(EM算法在此处与极大似然估计相同，都要先假设数据分布，不同的是，极大似然估计针对的是一个分布，而EM是针对的多种分布)。在这种情况下，估计样本的分布，我们从整个样本空间中抽取一个样本，有两个问题需要估计，首先这个人是男是女，其次是男生和女生满足的正太分布的参数，这两个是相互依赖的问题，如果我们知道了男女，那就可以很显然的使用极大似然估计的方法来估计男女的正太分布；反过来，我们知道了男女的分布，那么也知道了抽取到的样本是男女的哪个的可能性更大。

对于例子中的男女身高，这个人属于男生还是女生称之为隐含参数，男女的分布参成为模型参数。我们的目标是求出模型的参数，现在存在隐含参数对模型参数存在着影响。由此EM算法为：首先猜测隐含参数(EM中的E)，接着基于观察数据和猜测的隐含参数一起来极大化对数似然函数(EM中的M)。由于之前的隐含参数是猜测得到的，所以这个时候的模型参数是不可靠的，接下来基于得到的模型参数继续猜测隐含参数(EM中的E)，继续极大化对数似然函数，不断迭代下去直至算法收敛，找到合适的模型参数。

对于机器学习算法中的Kmeans算法就是使用的EM的思想，首先先假设K个样本值为最终的聚类簇的中心点(EM中的E)，之后计算样本中的其他点距离这些点的距离(EM中的M)，并将属于这个簇的所有点的平均值作为新的簇的中心，之后再进行距离计算，直至簇中心坐标不再变化位置。相关具体实现细节可以参考另一边文章机器学习-聚类

EM数学知识准备

凸函数

对于定义在实数域的函数，如果对于任意的实数，其二阶导数都存在$f^{‘’}(x) \geq 0$，这个函数称作凸函数。

Jensen不等式

因为该不等式在EM算法的推导过程中总会被用到，先了解一下该不等式

如果$f$是凸函数，$X$是随机变量那么
$$
E[f(x)]\ge f(E(x))
$$
如果$f$是严格的凸函数，那么$E[f(x)]=f(E(x))$，当且仅当$p(x=E[x])=1$时，也就是说随机变量$X$是一个常量

在图中可以看到对应着凸函数$f$，随机变量$X$有0.5的概率取到a，有0.5的概率取到b，则a,b的中值为对应的期望值，图中所示，满足Jensen不等式$E[f(x)] \ge f(E(x))$

同理，如果f是凹函数，则Jensen不等式反号$E[f(x)] \le f(E(x))$

期望

wiki定义

对于离散变量x的概率分布$p_{i}=p[X=x_{i}]$，数学期望$E(X)$为$E(X)=\sum_{i}x_{i}p_{i}$

其中$p_{i}$是权值满足条件$0 \le p_{i} \le 1$，$\sum_{i}p_{i}=1$

对于连续的随机变量，满足的概率密度函数为$f(x)$，则数学期望$E(X)$为
$$
E(X)=\int_{-\infty}^{\infty}xf(x)dx
$$
如果$Y=g(X)$，如果$X$是离散随机变量，则$Y$的期望值为
$$
E(Y)=\sum_{i}g(x_{i})p_{i}
$$
如果$X$是连续速记变量，则$Y$的期望值为
$$
E(Y)=\int _{-\infty}^{\infty}g(x)f(x)dx
$$

EM公式推导

抽取到的样本为$x=(x_{1},x_{2},…x_{m})$，隐含参数为$z=(z_{1},z_{2},…,z_{m})$，对应的模型参数为$\theta$，如果是极大似然函数的话，由于没有隐含参数的干扰，因此可以得到抽取个每个样本的概率为$p(x_{i}|\theta)$，对应的似然函数及优化目标为
$$
\theta: arg \max_{\theta} (H_{\theta})=\ln L(\theta)=\ln \prod_{i=1}^{m}p(x_{i}|\theta)=\sum_{i=1}^{m} \ln p(x_{i}|\theta)
$$
但是现在存在着隐含变量，因此对应的抽取每个样本的概率为$p(x_{i},z_{i}|\theta)$，其对应的似然函数及优化目标可以写成
$$
\theta, : z: arg \max_{\theta,z} (H_{\theta})=\ln L(\theta)=\ln \prod_{i=1}^{m}p(x_{i},z_{i}|\theta)=\sum_{i=1}^{m} \ln \sum_{z} p(x_{i},z_{i}|\theta)
$$
引入中间变量$Q(z_{i})$表示$z_{i}$的分布，可以对上面的式子进行整理
$$
\theta, : z: arg \max_{\theta,z} (H_{\theta})=\sum_{i=1}^{m} \ln \sum_{z} p(x_{i},z_{i}|\theta)=\sum_{i=1}^{m} \ln \sum_{z} Q(z_{i}) \frac{p(x_{i},z_{i}|\theta)}{Q(z_{i})}=\ln E(y)
$$
并且满足条件$1 \le Q(z_{i}) \le 1$及$\sum_z Q(z_{i})=1$

结合上面的数值期望的计算方法，可以将上面的式子看成$\ln E(y)$其中$y=Q(z_{i}) \frac{p(x_{i},z_{i}|\theta)}{Q(z_{i})}$

再结合Jensen不等式的定义，对于凸函数存在着不等式$E[f(x)]\ge f(E(x))$，此处的$f()=\ln()$，由此可以再次对上面的式子进行整理
$$
\begin{align}
\theta, : z: arg \max_{\theta,z} (H_{\theta})=\sum_{i=1}^{N} \ln \sum_{z_{i}} p(x_{i},z_{i}|\theta) &=\sum_{i=1}^{N} \ln \sum_{z_{i}} Q(z_{i}) \frac{p(x_{i},z_{i}|\theta)}{Q(z_{i})} =\ln E(y) \
& \ge \sum_{i=1}^{m}\sum_{z_{i}}Q(z_{i}) \ln (\frac{p(x_{i},z_{i}|\theta)}{Q(z_{i})})=E(\ln(y))
\end{align}
$$
由此，利用Jensen不等式，对优化目标函数进行了优化，得到了EM算法中的E步，求期望，找到了优化目标函数的下界。

为了找到使等号成立的条件，同样根据Jensen定义中得到等号成立的条件为随机变量是常数
$$
\begin{align}
\ln(y) = &\ln (\frac{p(x_{i},z_{i}|\theta)}{Q(z_{i})})=C \
\Longrightarrow & p(x_{i},z_{i}|\theta)=CQ(z_{i}) \
\Longrightarrow & \sum_{z}p(x_{i},z_{i}|\theta) = C\sum_{z}Q(z_{i})=C

\end{align}
$$
由此，可以得到
$$
Q(z_{i})=\frac{p(x_{i},z_{i}|\theta)}{C}=\frac{p(x_{i},z_{i}|\theta)}{\sum_{z}p(x_{i},z_{i}|\theta)} =\frac{p(x_{i},z_{i}|\theta)}{p(x_{i}|\theta)}=p(z_{i}|x_{i},\theta)
$$
其中$p(x_{i}|\theta)$是条件概率

至此，得到了隐含参数和观察变量及模型参数的关系，由此固定参数$\theta$，就可以根据观察变量来确定$Q(z_{i})$，由此确定Jensen不等式等号成立的条件，也就是优化目标函数的下界，这就是EM中的E步，然后根据固定的参数$Q(z_{i})$，来极大化优化目标$\ln (y)$调整模型参数$\theta$，当参数$Q(z_{i})$固定之后原始的优化目标变为
$$
\begin{align}
\theta, : z: arg \max_{\theta,z} (H_{\theta})=\sum_{i=1}^{N} \ln \sum_{z_{i}} p(x_{i},z_{i}|\theta) &=\sum_{i=1}^{N} \ln \sum_{z_{i}} Q(z_{i}) \frac{p(x_{i},z_{i}|\theta)}{Q(z_{i})} =\ln E(y) \
& \ge \sum_{i=1}^{m}\sum_{z_{i}}Q(z_{i}) \ln (\frac{p(x_{i},z_{i}|\theta)}{Q(z_{i})})=E(\ln(y)) \
& \Longrightarrow arg \max {\theta} \sum{i=1}^{m}\sum_{z_{i}}Q(z_{i}) \ln (p(x_{i},z_{i}|\theta))
\end{align}
$$

EM流程

根据给定的数据$x=(x_{1},x_{2},…,x_{m})$，求取联合概率分布$p(x,z|\theta)$和条件分布$p(z|x,\theta)$，极大迭代次数$J$

初始化模型参数$\theta$为$\theta^{0}$
for i from 1 to J:
- E步，计算联合概率分布的条件概率期望
  $$
  Q(z_{i}):=P(z_{i}|x_{i},\theta)
  $$
- M步，极大化对数似然函数
  $$
  \theta :arg \max {\theta} \sum{i=1}^{m}\sum_{z_{i}}Q(z_{i}) \ln (p(x_{i},z_{i}|\theta))
  $$
- 重复E,M直至$\theta$收敛
得到模型参数$\theta$

应用

如果我们从算法思想的角度来思考EM算法，我们可以发现我们的算法里已知的是观察数据，未知的是隐含数据和模型参数，在E步，我们所做的事情是固定模型参数的值，优化隐含数据的分布，而在M步，我们所做的事情是固定隐含数据分布，优化模型参数的值。EM的应用包括：

支持向量机的SMO算法
混合高斯模型
K-means
隐马尔可夫模型

举个栗子

本例子来自于参考文献1。扔硬币实验，存在两种硬币A,B，并进行5次实验，每次实验扔硬币10次(结果中H为正面，T为背面)，求硬币出现正面的概率。实验结果如下所示。

对于情况a，事先知道所选择的硬币进行实验，并得到5次实验的结果
对于情况b，实现不知道所选则的硬币，同样得到了5次实验的结果。

求解

对于情况a，显然可以利用极大似然估计的方法求得
$$
\begin{align}
arg \max {\theta{A},\theta_{B}} H(\theta) &=\log (\prod {i=1}^{5}p(x{i},\theta)) \
&=\log [\theta_{B}^{5}(1-\theta_{B})^5\theta_{A}^{9}(1-\theta_A)\theta_{A}^{8}(1-\theta_A)^2\theta_{B}^{4}(1-\theta_{B})^6\theta_{A}^{7}(1-\theta_A)^3] \
&= \log[\theta_{B}^{9}(1-\theta_{B})^{11}*\theta_{A}^{24}(1-\theta_{B})^6]
\end{align}
$$
对数似然函数求偏导
$$
\begin{align}
\frac{\partial H(\theta)}{\partial \theta_{B}}&=\frac{9}{\theta_{B}}+\frac{-11}{1-\theta_{B}}=0 \
\frac{\partial H(\theta)}{\partial \theta_{A}}&=\frac{24}{\theta_{A}}+\frac{-6}{1-\theta_{A}}=0
\end{align}
$$
由此可以求得
$$
\begin{align}
\theta_{A}&=\frac{24}{24+6}=0.8\
\theta_{B}&=\frac{9}{11+9}=0.45
\end{align}
$$
对于情况b，这是一个EM的应用，由于事先并不知道选择的硬币，因此按照EM算法，先假设初始的模型参数$\hat \theta_{A}=0.6,\hat \theta_{B}=0.5$，对于每次实验，先计算概率分布的条件概率$Q(z_{i})$
- 针对每次实验得到的结果分别计算对应的条件概率
  $$
  Q_{z_{i}}=p(z_{i} |x_{i},\theta)=\frac{p(z_{i},x_{i}|\theta)}{\sum p(z_{i},x_{i}|\theta)}
  $$
  根据示例中得到的实验结果，对每次实验分别都分别计算对应的条件概率
  $$
  Q_{z_{i}=A}=\frac {(\theta_{A})^{5}(1-\theta_{A})^{5}}{(\theta_{A})^{5}(1-\theta_{A})^{5}+(\theta_{B})^{5}(1-\theta_{B})^{5}}\
  Q_{z_{i}=B}=\frac {(\theta_{b})^{5}(1-\theta_{B})^{5}}{(\theta_{A})^{5}(1-\theta_{A})^{5}+(\theta_{B})^{5}(1-\theta_{B})^{5}}
  $$
  例如
  $$
  \begin{align}
  Q_{z_{1=A}}&=\frac{0.6^5*(1-0.6)^5}{0.6^5*(1-0.6)^5+0.5^5*(1-0.5)^5}=0.45\
  Q_{z_{1=B}}&=1-Q_{z_{1=A}}=0.55
  \end{align}
  $$
  由此可以得到每次实验的条件概率分布如下
  $$
  \begin{align}
  Q_{z_{1}=A}& =0.45,:::: & Q_{z_{1}=B}& =0.55 \
  Q_{z_{2}=A}& =0.80,:::: & Q_{z_{2}=B}& =0.20 \
  Q_{z_{3}=A}& =0.73,:::: & Q_{z_{3}=B}& =0.27\
  Q_{z_{4}=A}& =0.35,:::: & Q_{z_{4}=B}& =0.65\
  Q_{z_{4}=A}& =0.65,:::: & Q_{z_{4}=B}& =0.35\
  \end{align}
  $$
  至此完成了EM 算法的E步，接下啦进入M步，
  $$
  \begin{align}
  & \theta :arg \max {\theta} \sum{i=1}^{m}\sum_{z_{i}}Q(z_{i}) \ln (p(x_{i},z_{i}|\theta)) \
  \Longrightarrow & \sum_{i=1}^5 Q_{z_{i}=A}\log[(\hat \theta_{A})^{k_{i}}(1-\hat \theta_{A})^{10-{k_{i}}}]+ Q_{z_{i}=B}\log[(\hat \theta_{B})^{k_{i}}(1-\hat \theta_{B})^{10-{k_{i}}}]
  \end{align}
  $$
  其中$k_{i}$为第$i$次实验中正面朝上的次数。上面的优化目标函数求偏导并取值为0，由此可以得到更新之后的$\hat \theta_{A}$及$\hat \theta_{B}$，不断迭代直至模型参数收敛。