EM算法

EM算法(Expectation-Maximum)

原理篇

• 本质：概率模型，去估计一个密度函数，最大化对数似然函数去估计参数。无法直接用极大化对数似然函数得到模型分布的参数，用启发式跌代法，用于求解含有隐变量的最大似然估计、最大后延概率估计问题。

  1. 隐变量：对概率模型有一定的影响，但无法观测；

  ​ 观测数据X+隐变量 = 完全数据

  2. 无法用极大似然的两种情况：数据缺失、含有未知隐变量

• 作用：经常存在数据缺失或者不可用的的问题，这时候直接处理数据比较困难，而数据添加办法有很多种，常用的有神经网络拟合、添补法、卡尔曼滤波法等等，但是EM算法之所以能迅速普及主要源于它算法简单，稳定上升的步骤能非常可靠地找到“最优的收敛值”。

• ⽬标：是使包含隐变量的数据集的后验概率或似然函数最⼤化， 进⽽得到最优的参数估计

• 思想：我们可以发现我们的算法里已知的是观察数据，未知的是隐含数据和模型参数，在E步，我们所做的事情是固定模型参数的值，优化隐含数据的分布，而在M步，我们所做的事情是固定隐含数据分布，优化模型参数的值。

主要步骤：

已知：概率分布、随机抽取的样本；

未知：分类、模型参数

• E-step：猜想隐含数据，更新隐含数据和模型参数

• M-step：基于观察数据和猜测的隐含数据求极大对数似然，求解模型K

• E-step：基于前面得到的模型K，继续猜测隐含数据，继续极大化对数似然

• M-step：求模型参数

• 直至模型分布无明显变化，算法收敛

  EM算法步骤如下：

  输入：观测变量数据X，隐变量Z，联合分布p(X,Z|Θ)

  输出：模型参数Θ

  (1) 选择初始参数Θ_0

  (2)E步：记Θi为第i次迭代参数Θ的估计值，在第i+1次迭代的E步，计算Q(Θ,Θg);

  (3)M步：确定第i+1次迭代的参数的估计值Θ(i+1),即为：

  [Θ^{(i+1)}=argmax_ΘQ(Θ,Θ^{(i)}) ]

  (4)重复(2)和(3)步，直至收敛

EM算法的引入篇

有三枚硬币A、B、C

三枚硬币的模型：

[p(y| heta) = sum_zP(y,z| heta) = sum_zp(z| heta)p(y|z, heta)\ =pi p ^y(1-p)^{1-y} + (1-pi)q^y(1-q)^{1-y} ]

N枚硬币的模型：

[P(Y|Theta) = sum_zP(Z|Theta)P(Y|Z,Theta)\ = prod_{j = 1}^n[pi p^{y_j}](1-p)^{1-y_j} + (1-pi)q^{y_j}(1-q)^{q-y_j} ]

求模型的参数 Θ= (π,p,q)的极大似然估计，即目标函数为：

[Theta = argmaxlogp(Y|Theta) = argmaxlogsum_zP(Y|Z,Theta)p(Z,Theta) ]

面对含隐变量的概率模型，目标是极大化观测数据（不完全数据）Y关于参数Θ的对数似然函数，即极大化(取对数方便计算)：

[L(Theta) = logP(Y|Theta) = logsum_zP(Y,Z|Theta) \ = log(sum_zP(Y|Z,Theta)P(Z|Theta)) ]

EM算法通过迭代逐步近似极大化L(θ)的 + Jensen不等式

[L(Theta) - L(Theta^{i}) = log(sum_Z P(Y|Z,Theta)P(Z|Theta)) - logP(Y|Theta^{(i)})\ =log[sum_ZP(Z|Y,Theta^{(i)})frac{P(Y|Z)P(Z|Theta)}{P(Z|Y,Theta^{i})}] - logP(Y|Theta^{(i)})\ 由Jensen不等式且 sum_zP(Z|Y,Theta^{(i)}) = 1\ geq sum_ZP(Z|Y,Theta^{(i)})logfrac{P(Y|Z,Theta)P(Z|Theta)}{P(Z|Y,Theta^{(i)})} - frac{logP(Y|Theta^{(i)})*sum_ZP(Z|Y,Theta^{(i)})}{sum_ZP(Z|Y,Theta^{(i)})}\ = sum_Z P(Z|Y,Theta^{(i)})*logfrac{P(Y|Z,Theta)P(Z|Theta)}{P(Z|Y,Theta^{(i)})P(Y|Theta^{(i)})} ]

即：

[Theta^{(i+1)} = argmax [L(Theta^{(i)})+ sum_ZP(Z|Y,Theta^{(i)})*logfrac{P(Y|Z,Theta)P(Z|Theta)}{P(Z|Y,Theta^{(i)})P(Y|Theta^{(i)})}] ]

去掉与Θ无关的变量的式子：

[Theta^{(i+1)} = argmax [sum_ZP(Z|Y,Theta^{(i)})*logP(Y|Z,Theta)P(Z|Theta)]\ =argmaxQ(Theta,Theta^{(i)})= E[logP(Z|Y,Theta^{(i)})]....(1) ]

观测数据 + 隐变量 = 完全数据

Jensen不等式：

[log(sum alpha_i varphi(x_i)) >= sumalpha_i logvarphi(x_i)\ sum alpha_i = 1 且 alpha_i >= 0 ]

EM算法的收敛性

证明：P(Y|θ)为观测数据的似然函数，且是递增的，即：

[P(Y| heta^{(i+1)}) geq P(Y| heta^{(i)}) ]

证明如下：

[P(Y| heta) = frac{P(Y,Z| heta)}{P(Z|Y, heta)}......有贝叶斯公式展开得\ logP(Y| heta) = logP(Y.Z| heta) - logP(Z|Y, heta) ....取对数\ ]

[Q( heta, heta^{(i)}) = sum_z logP(Y,Z| heta)P(Z|Y, heta^{(i)}) ...由（1）得\ ]

构造下式（因为取对数方便相减和相除，同时构造了贝叶斯公式）：

[H( heta, heta^{(i)}) = sum_zlogP(Z|Y, heta)P(Z|Y, heta^{(i)})\ logP(Y| heta) = Q( heta, heta^{(i)}) - H( heta, heta^{(i)}).....(2) ]

在式（2）中分别取θ为 θi 和θ i+1，并相减：

[logP(Y| heta^{(i+1)}) - logP(Y| heta^{(i)}) = \ [Q( heta^{(i+1)}, heta^{(i)})- Q(( heta^{(i)}, heta^{(i)}))] - [H( heta^{(i+1)}, heta^{(i)}) - H(( heta^{(i)}, heta^{(i)}))] ]

其中对H：

[[H( heta^{(i+1)}, heta^{(i)}) - H(( heta^{(i)}, heta^{(i)}))] = sum_z(logP(Z|Y, heta^{(i+1})) =\ sum_z(logfrac{P(Z|Y, heta^{(i+1})}{P(Z|Y, heta^{(i)})}P(Z|Y, heta^{(i)})leq\ log(sum_Zfrac{P(Z|Y, heta^{(i+1})}{P(Z|Y, heta^{(i)})}P(Z|Y, heta^{(i)}))...Jensen不等式得\ log(sum_ZP(Z|Y, heta^{(i+1)}))=log1 = 0 ]

对Q，由于Q的i+1项已经达到极大，所以有：

[[Q( heta^{(i+1)}, heta^{(i)})- Q(( heta^{(i)}, heta^{(i)}))]geq0 ]

最后证明得到：

[logP(Y| heta^{(i+1)})geq logP(Y| heta^{(i)}) ]

应用篇--GMM

EM在GMM中的应用

图中可知：

1. 单个高斯拟合效果差，均值应该分布在数据密集处

2. 混合高斯模型中的隐变量，同时，隐变量在概率模型中不能改变边缘分布，即：

   [p(x_i)=int_{z_i}p(x_i|z_i)*p(z_i)dz_i = alpha_z = sum_{z_i=1}^kalpha_{z_i}N(x_i|mu_z,Sigma_z) ]

   每个数据都有一个隐变量，告诉你在哪个高斯模型中(由两个高斯扩展到n个高斯)

   [P(z_i = z_1|x_i,Theta^{g}) = frac{a}{a+b} ]

[P(x) = sum_{l=1}^{k}alpha_l*N(X|mu_l,Sigma_l) quadquadquad sum^k_{l=1}alpha_l = 1 ]

[Θ={α_1 ,…,α_k,μ_1,…,μ_k,Σ_1,…,Σ_{k-1}} ]

4. 目标函数为：
   [Θ_{MLE}=argmax_{Θ}*L(Θ∣X) =argmax_Θ(sum_{l=1}^nlog*sum_{l=1}^kα_lN(X∣μ_l,Σ_l)) ]
   该式子包含和（或积分）的对数，不能像单个高斯模型那样直接求导，再令导数为0来求解。这时我们需要利用 EM 算法通过迭代逐步近似极大化L(Θ∣X)来求解。

EM算法在GMM中的应用：

高斯混合模型的概率分布模型如下：

[P(Y| heta) = sum^{K}_{K=1}alpha_kphi(Y| heta)\ 其中：sum_{K=1}^{K}alpha_k = 1 , heta_k = (mu_k,sigma_k^2)\ heta = (alpha_1,alpha_2,...,alpha_k, heta_1, heta_2,..., heta_k)\ phi(Y| heta) = frac{1}{sqrt{2pi}sigma_k}exp(-frac{(y-mu_k)^2}{2sigma_k^2}) ]

用 γ 作为隐变量，去确定是哪个模型，γ =1/0.

完整数据的似然函数为：

[P(y,gamma| heta) = prod_{j=1}^NP(y_j,gamma_{j1},gamma_{j2},...,gamma_{jK}| heta)\ = prod_{K=1}^Kprod_{j=1}^N[alpha_kphi(y_j| heta_k)]^{gamma^{jK}}\ ]

完全数据的对数似然为：

[logP(y,gamma| heta) = sum_{K=1}^K[sum_{j=1}^Nlogalpha_k + sum_{j=1}^Ngamma_{jk}[log(frac{1}{sqrt{2pi}})-logsigma_k - frac{1}{2sigma^2_K}(y_j-mu_k)^2]] ]

EM算法中的E步：确定Q函数，对每一个隐变量求期望

根据当前模型参数，计算分模型k对观测数据y_j的响应度

[Q( heta, heta^{(i)}) = E[logP(y,gamma| heta)|y, heta^{(i)}]\ =E[sum_{K=1}^K[sum_{j=1}^Nlogalpha_k + sum_{j=1}^Ngamma_{jk}[log(frac{1}{sqrt{2pi}})-logsigma_k - frac{1}{2sigma^2_K}(y_j-mu_k)^2]]]\ =sum_{K=1}^K[sum_{j=1}^N(E_{gamma_{jk}})logalpha_K+sum_{j=1}^N(Egamma_{jk})[log(frac{1}{sqrt{2pi}})-logsigma_K-]frac{1}{2sigma_k^2}(y_j-mu_k)2] ]

计算E

[gamma_{jk} = E(gamma_{jk}|y, heta)\ =P(gamma_{jk}=1|y, heta)*1 + P(gamma_{jk}=0, heta)*0\ =frac{P(gamma_{jk}=1,y_j| heta)}{P(y_j| heta)} = frac{P(gamma_{jk}=1,y_j| heta)}{sum_{K=1}^KP(gamma=1,y_j| heta)}\ =...=frac{alpha_kphi(y_j| heta_k)}{sum_{K=1}^Kalpha_kphi(y_j| heta_k)} ]

最终求得Q

[Q( heta, heta^{(i)}) = sum_{K=1}^K[sum_{j=1}^Nlogalpha_k+sum_{j=1}^Ngamma_{jk}[log(frac{1}{sqrt{2pi}})-logsigma_k-frac{1}{2sigma_k^2}(y_j-mu_k)^2]] ]

M步,进行迭代模型的参数

[ heta^{(i+1)} = argmax_ heta Q( heta, heta^{(i)}) ]

[mu_k = frac{sum_{j=1}^Ngamma_{jk}y_j}{sum_{j=1}^Ngamma_{jk}}\ \ sigma^2_k = frac{sum_{j=1}^Ngamma_{jk}(y_j-mu_k)^2}{sum_{j=1}^Ngamma_{jk}}\ \ alpha_k= frac{sum_{j=1}^Ngamma_{jk}}{N} ]

参考1：知乎关于EM

参考2：GMM应用

参考3：EM九个境界

参考4：统计学习方法