EM算法(Expectation-Maximum)
原理篇
-
本质：概率模型，去估计一个密度函数，最大化对数似然函数去估计参数。无法直接用极大化对数似然函数得到模型分布的参数，用启发式跌代法，用于求解含有隐变量的最大似然估计、最大后延概率估计问题。
- 隐变量：对概率模型有一定的影响，但无法观测；
观测数据X+隐变量 = 完全数据
- 无法用极大似然的两种情况：数据缺失、含有未知隐变量
-
作用：经常存在数据缺失或者不可用的的问题，这时候直接处理数据比较困难，而数据添加办法有很多种，常用的有神经网络拟合、添补法、卡尔曼滤波法等等，但是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中的应用
![高斯分布分析]()
图中可知：
-
单个高斯拟合效果差，均值应该分布在数据密集处
-
混合高斯模型中的隐变量，同时，隐变量在概率模型中不能改变边缘分布，即：
[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}}
]
- 目标函数为：
[Θ_{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：统计学习方法