机器学习小白笔记系列——EM算法

本系列是针对于DataWhale学习小组的笔记，从一个对统计学和机器学习理论基础薄弱的初学者角度出发，在小组学习资料的基础上，由浅入深地对知识进行总结和整理，今后有了新的理解可能还会不断完善。由于水平实在有限，不免产生谬误，欢迎读者多多批评指正。如需要转载请与博主联系，谢谢

EM算法核心概念

---

什么是EM算法？

EM算法是最大期望算法的英文缩写。它是一种利用迭代来完成极大似然估计或极大后验概率估计的算法，主要用于对包含隐变量或缺失数据的概率模型进行参数估计。实际上，EM算法更接近于一种算法思路而非模型，它的应用往往需要与某些模型相结合，比如对高斯混合分布的估计。

EM算法的基本运用场景就是存在不可测量的未知变量的概率模型，如同统计学习方法书中举的例子，有三个硬币A,B和C，其各抛一次为正面朝上的概率分别为π,p和q（这里确实有点奇怪哈，竟然不是0.5）或写作θ=(π,p,q)，先抛硬币A，由A的正反来决定接下来抛B还是C，第二次抛出的硬币则决定输出（正为1，反为0）。这样经过几次实验后可以获得输出值组成的序列（如1,1,0,1,0,0,1,1等），问假设只能观测到掷硬币的结果，而不能观测掷硬币的过程，如何估计三个硬币分别出现正面的概率π,p和q？从这个问题中我们可以发现，我们的样本仅知道观测变量Y=0或1（由第二轮B或C抛出的结果），但不知道其中包含的隐变量Z（由第一轮A抛出的结果）。在这里我们将实验中的观测数据（样本中的已知变量取值）和未观测数据（样本中的隐变量取值）加起来被称为完全数据。(P(Y|θ))代表当我们的模型参数θ——也就是三个硬币抛出后正面朝上的概率确定时，实验中出现训练集中所有观测样本的概率（也就是当前模型参数下所有单个样本出现概率(P(y_i|θ))的乘积），我们要做的就是通过不断尝试不同的θ让这个概率最大化。这里也很好理解，这个过程就相当于我们通过不断调整自己的判断倾向（即参数θ）使得自己心中的预测方案与已有的事实Y尽可能贴近。根据这种思路，我们可以写出：

[egin{aligned} P(Y| heta)&=sum_{Z}P(Z| heta)P(Y|Z, heta)\ &=prod_{i=0} ( pi p^{y_i}(1-p)^{1-y_{i}}+(1-pi)q^{y_{i}}(1-q)^{1-y_{i}}) end{aligned} ]

我们的目标是使这个概率最大化：

[hat{ heta}=argmax_{ heta}logP(Y| heta) ]

这个问题通常没有解析解，因此我们选择利用迭代的方式找到θ的最佳取值。

EM算法基本步骤

EM算法的核心就是通过迭代求出$L(θ)=logP(Y|θ)的极大似然估计。
算法输入：可观测变量Y，隐变量Z，联合分布P(Y,Z|θ)，条件分布P(Z|Y,θ)。
算法输出：模型参数θ。
算法步骤：

1. 随机设置一组参数的初值( heta ^{left ( 0 ight )})，但注意EM算法对于初值敏感，可能导致获得局部最优解，解决办法就是多取几组初值看最终结果，比较选择出最优解。
2. E步骤（Expectation）：令( heta ^{left ( i ight )})为第i次迭代后参数的估计值，则在第i+1次迭代中，利用当前参数估计值来计算Q函数（Q函数根据书上解释是完整数据的对数似然函数(logP(Y,Z| heta))关于在给定观测数据Y和当前参数( heta ^{left ( i ight )})下对未观测数据Z的条件概率分布(P(Z|Y, heta ^{left ( i ight )}))的期望，由于直接对含有隐变量/未观测数据的对数似然函数求极大值很困难，我们通过极大化Q函数来代替，可以证明后者是对前者的逼近，在此暂不讨论证明过程）：

[egin{aligned} Q( heta, heta^i)&=E_{Z}[logP(Y,Z| heta)|Y, heta^i]\ &=sum_{Z}logP(Y,Z| heta)P(Z|Y, heta^i) end{aligned} ]

3. M步骤（Maximization）：求使Q函数极大化的参数θ作为下一次迭代的参数估计( heta ^{left ( i+1 ight )})，可以证明每次经过此过程都相当于使似然函数增大：

[ heta^{i+1}=arg max limits_{ heta}Q( heta, heta^{i}) ]

4. 判断是否停止迭代，判断依据是当前参数θ或者Q函数与其上轮循环数值之差的模是否小于一个阈值，即：

[|| heta^{i+1}- heta^{i} || < varepsilon_1 ]

或者：

[||Q( heta^{i+1}, heta^{i})-Q( heta^{i}, heta^{i})|| <varepsilon_2 ]

EM算法在高斯混合模型上的应用

EM算法是处理高斯混合模型的有效方法，高斯混合分布模型顾名思义就是由多个参数不同的高斯分布按一定的权重混合而成，其概率分布如下：

[P(y| heta)=sum_{k=1}^{K}a_kphi(y| heta_{k}) \ phi(y| heta_{k})=frac{1}{sqrt{2pi}delta_{k}}exp(-frac{(y-mu_{k})^2}{2 delta_{k}^{2}}) \ a_k>0,sum a_k =1 ]

其中系数(a_k)大于等于0且总和归一化为1，(phi(y| heta_{k}))为第k个高斯分布密度，( heta _{k}=left ( u_{k},sigma _{2}^{k} ight ))。
为了利用EM算法求解高斯混合模型参数，我们可以将模型这样理解，第一步利用概率(a_k)选择第k个单一高斯分布模型(phi(y| heta_{k}))，第二步根据选中的第k个分模型的概率分布(phi(y| heta_{k}))来生成观察数据(y_j)，其中j=1,2……,N;此时任意观测数据(y_j)都对应一个未观测的第一步输出k=1,2……,K，以隐变量(y_{jk})表示，完全数据如下：

[left ( y_{j},gamma _{j1},gamma _{j2},cdots,gamma _{jK} ight ), j=1,2,cdots,N ]

为了计算已知y及θ时(gamma _{jk})的期望E，我们定义(hat{gamma_{jk}})为当前模型参数下第j个观测数据来自第k个分模型的概率，称为分模型k对观测数据(y_j)的响应度。

至此我们可以通过推导得出EM方法在高斯混合模型中关键步骤的计算形式（推导过程可参考《统计学习方法》中相关章节），即我们前面写过的第二步——E步骤（求期望）及第三步——M步骤（求极大），其余为参数赋初值及迭代终止阈值的设置的步骤仍根据实际要求完成：
E步骤：依据当前模型参数，计算响应度

[hat{gamma_{jk}} = E(gamma_{jk}|y, heta)=frac{a_kphi(y_i| heta_{k})}{sum_{k=1}^{K}a_kphi(y_i| heta_{k}) }\ j=1,2,..,N；k=1,2,...,K\ n_k=sum_{j=i}^{N}Egamma_{jk} ]

M步骤：构造及化简Q函数，对Q函数求导，令偏导等于0得到三个参数的迭代公式（这里写出只是为了从数学上便于理解函数迭代的意义，但在实际计算机编程中不需要写，直接对下面的三个参数进行更新即可）：

[egin{aligned} Q( heta. heta^i) &= E[log(P(y,gamma| heta))]\ &=sum_{K}^{k=1}igg(sum_{j=1}^{N}(Egamma_{jk}) log(a_k)+sum_{j=1}^{N}(Egamma_{jk})igg[log(frac{1}{sqrt{2pi}})-log(delta_{k})-frac{(y_i-mu_{k})^2}{2 delta_{k}^{2}}igg]igg) end{aligned} ]

θ的三个参数更新如下：

[hat{mu_k}=frac{sum_{j=1}^{N}hat{gamma_{jk}}y_i}{sum_{j=1}^{N}hat{gamma_{jk}}} \ \ hat{delta_k}=frac{sum_{j=1}^{N}hat{gamma_{jk}}(y_i-mu_k)^2}{sum_{j=1}^{N}hat{gamma_{jk}}} \ \ \ hat{a_k}=frac{sum_{j=1}^{N}hat{gamma_{jk}} }{N} ]

然后重复这两个过程，直至满足迭代终止判定即可。

参考资料：

1. 《统计学习方法》李航著
2. https://github.com/datawhalechina/team-learning/tree/master/机器学习算法基础 DataWhale小组学习资料
3. https://baike.baidu.com/item/最大期望算法/10180861?fromtitle=em算法&fromid=1866163&fr=aladdin 百度百科
4. https://www.jianshu.com/p/c57ef1508fa7 你真的了解EM算法吗