十大之——EM算法（Expectation-Maximum）

一. 为什么引入EM

1. 最大似然估计（MLE）

（1）已知：① 一堆观测数据X

② 数据服从的统计模型

估计：统计模型中的参数

（2）引入MLE： $\hat{\theta} _{MLE} =argMax\sum_{i=1}^n \lg P(x_{i};\theta)$

（3）非凸问题：无法求全局极值，只能求局部极值。

① 梯度下降法； ② 启发式方法； ③ EM算法

2. EM估计

（1）问题：① 观测数据X有缺失

② 难以得知 $P(X|\theta )$

方法：引入隐含变量Z，得 $∑P(X,Z|\theta )$ （可以简化模型，不改变原分布的边缘）

（2）优势：无需调参、编程简单、理论完整

二. EM算法（聚类）

1. 原理

（1）反复局部下限构造：巧妙的构造了一个 $P(X|\theta )$ 的下界，通过优化这个相对简单的下界来优化最终目标 $P(X|\theta )$

（2）启发式的迭代方法：先固定模型参数和观测值估计隐变量（E步）；再固定得到的隐变量，将其和观测数据一起通过极大化似然函数来估计模型参数（M步）。依次迭代直至模型参数变化小于指定阈值。

2. 推导

（1）凸函数：对于任意的实数，都有 $f{″}≥0$

对于由实数组成的向量，都有矩阵 $H″≥0$ （半正定）

Jensen不等式：E[f(X)] ≥ f(E[X])

（2）完整数据 $(X,Z)$ ：给定m个相互独立的观测数据 $x=(x^{(1)},..., x^{(m)})$

对应的未观测到的隐含变量值 $z=(z^{(1)},..., z^{(m)})$

（3）似然函数： $\hat{\theta}=argMax\sum_{i=1}^m \lg P(x^{(i)}|\theta)$

加入隐含变量： $\hat{\theta},\hat{z}=argMax l(\theta ,z)=argMax\sum_{i=1}^m \lg \sum_{z^{(i)}} P(x^{(i)},z^{(i)}|\theta)$

（4）设 $Q_i(z^{(i)})≥0$ 表示该样例隐含变量z的某种分布，且满足 $\sum_{z} Q_i(z)=1$ ，则：

$l(\theta,z)=\sum_{i=1}^m \lg \sum_{z^{(i)}} Q_i(z^{(i)})·\frac{P(x^{(i)},z^{(i)}|\theta)}{Q_i(z^{(i)})}≥\sum_{i=1}^m \sum_{z^{(i)}} Q_i(z^{(i)})·lg\frac{P(x^{(i)},z^{(i)}|\theta)}{Q_i(z^{(i)})}=J(z,Q)$

用到 Jensen 不等式： $lg[E(y)]≥E[lg(y)]$ （对数函数为凹函数）

（5） $J(z,Q)$ 中 $\sum_{z^{(i)}} Q_i(z^{(i)})·lg\frac{P(x^{(i)},z^{(i)}|\theta)}{Q_i(z^{(i)})}=E[\frac{P(x^{(i)},z^{(i)}|\theta)}{Q_i(z^{(i)})}]$

Lazy Statistician规则：

① X是离散随机变量，分布律为 $P(X=x_k)=p_k$ ，则： $E(Y)=E[g(X)]=\sum_{k=1}^∞g(x_k)p_k$

② X是连续型变量，他的概率密度为g(x)。则： $E(Y)=E[g(X)]=\int_{−∞}^{∞}g(x)f(x)dx$

（6）局部下限构造：由于 $l(θ,z)≥J(z,Q)$ ，只要不断最大化右侧下界（寻找 $Q_i(z)$ ），使得 $l(θ,z)$ 不断提高，就能在两边取等号时，使左侧达到最大值。若固定θ，那么 $l(θ,z)$ 的值就取决于 $Q_i(z^{(i)})$ 和 $p(x^{(i)},z^{(i)})$ 。

（7）等式条件：随机变量是常数

即： $g(X_k)= \frac{p(x ^{(i)} ,z ^{(i)} ∣θ)}{Q_i (z ^{(i)} )} =C$

两边同时累加： $\sum_{z}p(x ^{(i)} ,z ^{(i)} ∣θ)=C·\sum_{z} Q_i (z ^{(i)} )=C$

（8）E-step： $Q_i (z ^{(i)} )=\frac{p(x ^{(i)} ,z ^{(i)} ∣θ)}{ \sum_{z}p(x ^{(i)} ,z ^{(i)} ∣θ)}=\frac{p(x ^{(i)} ,z ^{(i)} ∣θ)}{ p(x_i|\theta )}=p(z^{(i)}|x^{(i)},\theta )$

联合概率： $P(X=x_i,Y=y_j)=p_{ij}$

边缘概率： $p_{i.}=\sum_{j=1}^∞ p_{ij}$

条件概率： $P(X=x_i|Y=y_j)=\frac{ P(X=x_i,Y=y_j)}{P(Y=y_j)}$

（9）M-step： $Q_i(z^{(i)})$ 已知，上式 $J(z,Q)$ 中只剩 $θ$ 一个未知量，因此可以直接使用极大似然法获得一个新的 $θ$ 值，再用这个新 $θ$ 重新进行E过程，直到收敛。

$\hat{\theta }=argMax \sum_{i=1}^m \sum_{z^{(i)}} Q_i(z^{(i)})·lg\frac{P(x^{(i)},z^{(i)}|\theta)}{Q_i(z^{(i)})}$

3. EM算法流程

（1）输入：观察数据 $x=(x^{(1)},..., x^{(m)})$ 、联合分布 $p(x,z|\theta )$ 、条件分布 $p(z|x,\theta )$

① 初始化参数的初值 $\theta ^{(0)}$ （任意，但敏感）

② E-step： $Q_i (z ^{(i)} )=p(z^{(i)}|x^{(i)},\theta )$

③ M-step： $\hat{\theta }=argMax \sum_{i=1}^m \sum_{z^{(i)}} Q_i(z^{(i)})·lg\frac{P(x^{(i)},z^{(i)}|\theta)}{Q_i(z^{(i)})}$

④ 重复E、M步骤直到 $\theta$ 收敛

（2）输出：模型参数 $θ$

三. EM在高斯混合模型中的应用（Gaussian Misture ,Model）

1. 高斯分布（Gaussian Distribution）

也被称为正态分布（normal distribution），是一种在自然界大量的存在、最为常见的分布。高斯分布的概率密度计算核心在于计算数据点到中心的距离，并且除以标准差将这个绝对距离转化为相对距离，然后通过距离平方的指数衰减计算概率密度。

（1）一维概率密度： $f(x|μ,σ^ 2)= \frac{1}{\sqrt{2π}σ} exp(− \frac{(x−μ)^ 2}{2σ^ 2})$

X∼N(μ,σ)，其中μ为均值，σ为标准差，σ²为方差

关于μ对称；总面积为1；在μ加减σ处为拐点（先内翻后外翻）

（2）多维联合概率密度： $f(x|\mu ,\Sigma )= \frac{1}{\sqrt{(2π)^d|Σ| }}exp[− \frac{1}{2} (x−\mu ) ^T Σ ^{−1} (x−\mu)]$

X∼N(μ,Σ)，其中μ为均值，Σ为协方差，|Σ|为行列式

协方差（Covariance）衡量两个变量的总体误差。方差是两个变量相同的特殊情况。协方差矩阵 $Σ=E[(x-\mu )(x-\mu)^T]$ ，由于假设了各个维度之间不相关，因此协方差矩阵只有在对角线的位置有值， $\sigma _1\sigma _2...\sigma _d=|\Sigma |^{\frac{1}{2} }$

2. 高斯混合模型（Gaussian Mixture Model）

（1）定义：通过求解K个高斯模型，并通过一定的权重将K个高斯模型融合成一个模型

（2）概率密度： $p(x)= \sum_{k=1}^K p(k)p(x|k)=\sum_{k=1}^K π_ kN(x|u _k ,Σ _k)$

$𝜋_𝑘$ ：权值因子，第k个模型的先验概率，且 $\sum_{k=1}^K 𝜋_𝑘=1$ （概率密度函数在其作用域内的积分之和必然为1）； $N(x|u_k,Σ_k)$ ：第k个高斯模型的概率密度函数

3. 求解GMM模型

（1）极大似然法：

① 似然函数：求样本集出现的联合概率，但无法得知样本到底来源于哪一类的高斯模型。

② 引入隐变量：一个K维二值随机变量 $\gamma =[\gamma _1,...\gamma _k]$

$γ_k\in$ {0,1}，且 $\sum_{k=1}^K γ_k=1$ ， $\gamma _k=1$ 表示样本x由分模型k抽样产生，共有K种可能，即每一次采样，选择第k个高斯模型的概率为： $p(γ _k =1)=π _k$

③ 根据 $γ$ 的取值，确定选择第k个高斯模型进行采样，满足： $p(x|γ _k =1)=N(x|μ _k ,Σ _k )$

④ 样本x的边缘分布： $p(x)=\sum_{\gamma } p(γ)·p(x|γ)= \sum_{k=1}^Kπ _k· N(x|μ_ k ,Σ _k )$

⑤ 样本集联合概率： $L(μ,Σ,π)=L(x_1...x_N;μ,Σ,π)=\prod_{n=1}^N p(x_ n ;μ,Σ,π)= \prod_{n=1}^N\sum_{k=1}^K π _kN(x _n|μ _k ,Σ_ k )$

⑥ 对数似然函数： $lnL(μ,Σ,π)=\sum_{n=1}^Nln[ \sum_{k=1}^K π _kN(x _n|μ _k,Σ_ k )]$

（2）使用EM算法求解GMM参数

① 单个完全数据的似然函数：设 $\gamma _{t,k}=1$ （第t个观测变量来自第k个高斯分量）

$p(x_j,γ|μ,Σ,π)=\prod_{k=1}^K [π_ kN(x_j ;μ _k,Σ _k )]^{γ _{j,k}}=π_tN(x_j ;μ _t,Σ _t)$

② 所有观测数据的完全数据似然函数：

$p(x,γ|μ,Σ,π)=\prod_{j=1}^N \prod_{k=1}^K [π_ kN(x_j ;μ _k,Σ _k )]^{γ _{j,k}}p(y,γ∣μ,Σ,π)=\prod_{k=1}^K π_k^{\sum_{j=1}^N γ_{j,k}}\prod_{j=1}^N N(x_j;μ _k,Σ _k)^{γ _{j,k}}$

③ 完全数据的对数似然函数：

$l(x;Θ)=lnp(x,\gamma|\mu ,\Sigma,\pi )=\sum_{j=1}^N ln\sum_{k=1}^KQ_j(\gamma _{j,k}) \frac{p(x_j|\gamma_{j,k} =k;\mu ,\Sigma,\pi)p(\gamma _{j,k}=k)}{Q_j(\gamma _{j,k})}$

$lnp(x,\gamma|\mu ,\Sigma,\pi )=\sum_{j=1}^N\sum_{k=1}^K[γ_{j,k}ln\pi _k+γ_{j,k}lnN(x_j|μ_k,Σ_k)]$

$lnN(x_j|μ_k,Σ_k)=−\frac{D}{2} ln(2π)−\frac{1}{2} ln|Σ_k|−\frac{1}{2}(x_j−μ_k)^TΣ_k^{−1}(x_j−μ_k)$

$lnp(x,γ|μ,Σ,π)= \sum_{k=1}^K ( \sum_{j=1}^N γ _{j,k} )lnπ_k+ \sum_{j=1}^N γ _{j,k} [−\frac{D}{2} ln(2π)− \frac{1}{2} ln|Σ|− \frac{1}{2} (x_j−μ) ^T (Σ)^{−1}(x_j −μ)]$

④ 设Q函数表示隐含变量 $\gamma$ 的高斯分布： $Q_j(\gamma _{j,k})=P(\gamma _{j,k}=k|x_j;μ,Σ,\pi )$

$l(x;Θ)≥ \sum_{j=1}^N \sum_{k=1}^K Q _j(\gamma _{j,k}) ln\frac{p(x_j|\gamma_{j,k} =k;\mu ,\Sigma,\pi)p(\gamma _{j,k}=k)}{Q _j(\gamma _{j,k})} =J(\gamma ,Q)$

⑤ 最大化 J 使得上式取等号：

$J(Θ,Θ^i)=\sum_{j=1}^N \sum_{k=1}^KQ_j(\gamma _{j,k})ln\frac{\pi_kN(x_j|Θ_k)}{Q_j(\gamma _{j,k})} =\sum_{k=1}^K(\sum_{j=1}^N\hat{\gamma } )ln\pi_k+ \sum_{j=1}^N\hat{\gamma } (−ln(2π)− \frac{1}{2}ln|Σ_k|− \frac{1}{2}(x_j−μ_k)^T(Σ_k) ^{−1}(x_j−μ_k))$

⑥ E-step：假设第j个观测来自第k个component的概率，对 $\gamma$ 进行估计

$\hat{\gamma } =E(γ_{j,k}|x;Θ)=P(γ_{j,k}=1|x,Θ)=\frac{P(γ_{j,k}=1,x_j|Θ)}{p(x_j)}=\frac{P(γ_{j,k}=1,x_j|Θ)}{\sum_{k=1}^K P(γ_{j,k}=1,x_j|Θ)}=\frac{\pi_kN(x_j|Θ_k)}{\sum_{k=1}^K\pi_kN(x_j|Θ_k)}$

⑥ M-step：假设我们已经知道隐变量 $γ^{j,k}$ 的取值，求偏导后参数估计为

$μ_k=\frac{ \sum_{j=1}^N γ_{j,k}x_j}{\sum_{j=1}^Nγ_{j,k}}$ ， $Σ_k=\frac{\sum_{j=1}^Nγ_{j,k}(x_j−μ_k)(x_j−μ_k)^T }{ \sum_{j=1}^N γ_{j,k}}$ ， $\pi _k=\frac{\sum_{j=1}^N γ^{j,k}}{N}$