EM算法与高斯混合模型

EM算法：
解决有隐变量的优化问题。如：
设想第k个观测数据\(y_k\)是这样产生的，首先依据概率选择第k个高斯分布模型，然后依第k个分模型的概率分布生成数据\(y_K\)。

1.高斯混合模型
一般地，用Y表示观测随机变量的数据，Z表示隐随机变量的数据。Y和Z连在一起称为完全数据，观测数据Y又称为不完全数据，假设给定观测数据Y，其概率分布是\(P(Y|\theta)\),对数似然函数\(L(\theta)=logP(Y|\theta);\)假设Y和Z的联合概率分布是P(Y,Z|\theta),那么完全数据的对数似然函数就是\(logP(Y,Z|\theta)\)
EM算法通过迭代求\(L(\theta)=logP(Y|\theta)\)的极大似然估计，每次迭代包含两步：E步，求期望；M步，求极大化

例子(三硬币模型)
假设有3枚硬币，分别记作A,B,C，这些硬币正面出现的概率分别是\(\pi,p,q\)，进行如下抛硬币试验:先抛硬币A，若为正面，再抛B，否则抛C；然后掷选出的硬币，掷硬币的结果，出现正面记为1，出现反面为0，独立地重复n此试验，观测结果为:1,1,0,1,0,0,1,0,1,1
假设只能观测到掷硬币的结果，不能观测掷硬币的过程。问如何估计三硬币正面出现的概率\(\pi,p,q\)

解:三硬币模型可以写作
\(p(y|\theta)=\sum_zP(y,z|\theta)=\sum_zP(z|\theta)P(y|z,\theta)=\pi p^y(1-p)^{1-y}+(1-\pi)q^y(1-q)^{1-y}\)
这里随机变量y是观测变量，表示一次试验观测的结果是1或0.随机变量z是隐变量，表示未观测到的抛硬币A的结果；\(\theta=(\pi,p,q)\)是模型参数。随机变量y的数据是可以观测的，随机变量z的数据不可观测
考虑求模型参数\(\theta=(\pi,p,q)\)的极大似然估计，即\(\hat{\theta}=argmax_{\theta}logP(Y|\theta)\)
该问题没有解析解，只有通过迭代的方法求解，EM算法是其中一种迭代求解方法
EM算法首先取参数的初值，记为\(\theta^{(0)}=(\pi^{(0)},p^{(0)},q^{(0)})\),然后通过下面的步骤迭代计算参数的估计值，直至收敛。第i次迭代参数的估计值为\(\theta^{(i)}=(\pi^{(i)},p^{(i)},q^{(i)})\),EM算法的第i+1次迭代如下
E步：计算在模型参数\(\pi^{(i)},p^{(i)},q^{(i)}\)下观测数据\(y_j\)来自掷硬币B的概率(来自各子模型概率)
\(\mu^{(i+1)}=\frac{\pi^{(i)}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j}}{\pi^{(i)}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j}+(1-\pi)^{(i)}(q^{(i)})^{y_j}(1-q^{(i)})^{1-y_j}}\)
M步：计算模型参数的新估计值
\(\pi^{(i+1)}=\frac{1}{n}\sum_{j=1}^n\mu_j^{(i+1)}\\ p^{(i+1)}=\frac{\sum \mu_j^{(i+1)}y_j}{\sum\mu_j^{(i+1)}} \\ q^{(i+1)}=\frac{\sum (1-\mu_j^{(i+1)})y_j}{\sum(1-\mu_j^{(i+1)})}\)

EM算法
输入：观测变量数据Y,隐变量数据Z，联合分布\(P(Y,Z|\theta)\),条件分布\(P(Z|Y,\theta)\)
输出：模型参数\(\theta\)
1.选择参数的初值\(\theta^{(0)}\),开始迭代
2.E步:记\(\theta^{(i)}\)为第i次迭代参数\(\theta\)的估计值，在第i+1次迭代的E步，计算

\[Q(\theta,\theta^{i})=E_Z[logP(Y,Z|\theta)|Y,\theta^{(i)}]\\ =\sum_ZlogP(Y,Z|\theta)P(Z|Y，\theta^{(j)}) \]

3.M步：求使\(Q(\theta,\theta^{i})\)极大化的\(\theta\),确定第i+1次迭代的参数的估计值\(\theta^{(i+1)}\)=\(argmax_{\theta}Q(\theta,\theta^{(i)})\)
4.重复2,3直至收敛

Q函数
完全数据的对数似然函数\(logP(Y,Z|\theta)\)关于在给定观测数据Y和当前参数\(\theta^{(i)}\)下对未观测数据Z的条件概率分布\(P(Z|Y,\theta^{(i)})\)的期望称为Q函数，即\(Q(\theta,\theta^{(i)})=E_Z[logP(Y,Z|\theta)|Y,\theta^{(i)}]=\sum logP(Y,Z|\theta)P(Z|Y,\theta^{(i)})\)

EM算法的导出
当我们面对一个含有隐变量的概率模型，目标是极大化观测数据(不完全数据)Y关于参数\(\theta\)的对数似然函数，即极大化\(L(\theta)=logP(Y|\theta)=log\sum_ZP(Y,Z|\theta)=log(\sum_ZP(Y|Z,\theta)P(Z|\theta))\)极大化这一式子是很困难的.事实上，EM算法是通过迭代逐步近似极大化\(L(\theta)\)的，假设在第i次迭代后\(\theta\)的估计值是\(\theta^{(i)}\)。我们希望新估计值\(\theta\)能使\(L(\theta)\)增加，即\(L(\theta)>L(\theta^{(i)})\)，并逐步达到极大值，为此，我们考虑两者的差：
\(L(\theta)-L(\theta^{(i)})=logP(Y|\theta)-logP(Y|\theta^{(i)})=log(\sum_ZP(Y|Z,\theta)P(Z|\theta))-logP(Y|\theta^{(i)})\)

利用Jensen不等式得其下界：
\(L(\theta)-L(\theta^{(i)})=log(\sum_ZP(Y|Z,\theta^{(i)})\frac{P(Y|Z,\theta)P(Z|\theta))}{P(Y|Z,\theta^{(i)})})-logP(Y|\theta^{(i)})\\ \geq \sum_ZP(Y|Z,\theta^{(i)})log\frac{P(Y|Z,\theta)P(Z|\theta))}{P(Y|Z,\theta^{(i)})}-logP(Y|\theta^{(i)}\\ =\sum_ZP(Y|Z,\theta^{(i)})log\frac{P(Y|Z,\theta)P(Z|\theta))}{P(Y|Z,\theta^{(i)})P(Y|\theta^{(i)})}\)
令\(B(\theta,\theta^{(i)})=L(\theta^{(i)})+\sum_ZP(Y|Z,\theta^{(i)})log\frac{P(Y|Z,\theta)P(Z|\theta))}{P(Y|Z,\theta^{(i)})P(Y|\theta^{(i)})}\)
则\(L(\theta)\geq B(\theta,\theta^{(i)})\),因此任何可以使\(B(\theta,\theta^{(i)})\)增大的\(\theta\)，也可以使\(L(\theta)\)增大，为了使其进可能增大，选择\(\theta^{(i+1)}\)使\(B(\theta,\theta^{(i)})\)达到极大，即

\[\theta^{(i+1)}=argmax_{\theta}B(\theta,\theta^{(i)}) \]

而\(\theta^{(i)}=argmax_{\theta}(\sum P(Z|Y,\theta^{(i)})log(P(Y|Z,\theta)P(Z|\theta) )\\ =argmax_{\theta}(\sum P(Z|Y,\theta^{(i)})log(P(Y,Z|\theta) \\ =argmax_{\theta}Q(\theta,\theta^{(i)})\)

EM算法那在高斯混合模型学习中的应用
高斯混合模型是指具有如下形式的概率分布模型：
设想第k个观测数据\(y_k\)是这样产生的，首先依据概率选择第k个高斯分布模型，然后依第k个分模型的概率分布生成数据\(y_K\)。
\(P(y|\theta)=\sum a_k\phi(y|\theta_k)\),其中\(\alpha_k\geq 0 , \sum \alpha =1 ;\phi(y|\theta_k)\)是高斯分布密度,\(\theta_k=(\mu_k,\sigma_k^2)\)
\(\phi(y|\theta_k)=\frac{1}{\sqrt{2\pi} \sigma_k}exp(-\frac{(y-mu_k)^2}{2\sigma_k^2})\)称为第k个分模型
一般混合模型可以由任意概率分布密度代替上式中的高斯分布密度

我们最终给出的只是某一个观测是由某一个高斯分量唯一生成的概率值，而不是确定下来的属于某一类。每个高斯分量其实都可以产生这个观测数据只是输出不同而已，即产生观测数据的概率不同。最后，根据每个高斯分量产生观测数据的可能性不同，结合其权值汇总出整个GMM产生这个观测数据的概率值。

高斯混合模型参数估计的EM算法
输入：观测数据\(y_1,y_2,\cdots,y_N\),高斯混合模型
输出：高斯混合模型参数
(1)取参数的初始值开始迭代
(2)E步：依据当前模型参数，计算分模型k对观测数据\(y_j\)的响应度
\(\hat{\gamma}_{jk}=\frac{\alpha_k\phi(y_j|\theta_k)}{\sum\alpha_k\phi(y_j|\theta_k)}\)
(3)M步：计算新一轮迭代的模型参数
\(\hat{\mu}_k=\frac{\sum_{j=1}^N \hat{\gamma}_{jk}y_j }{\sum_{j=1}^N\hat{\gamma}_{jk}}\)
\(\hat{\sigma}_k^2=\frac{\sum_{j=1}^N \hat{\gamma}_{jk}(y_j -\mu_k)^2 }{\sum_{j=1}^N\hat{\gamma}_{jk}},k=1,2,\cdots,K\)
\(\hat{\alpha}_k=\frac{\sum_{j=1}^N \hat{\gamma}_{jk} }{N}\)
(4)重复(2)和(3)，直至收敛

GMM和K-means直观对比

最后我们比较GMM和K-means两个算法的步骤。

GMM：
先计算所有数据对每个分模型的响应度
根据响应度计算每个分模型的参数
迭代

K-means：
先计算所有数据对于K个点的距离，取距离最近的点作为自己所属于的类
根据上一步的类别划分更新点的位置（点的位置就可以看做是模型参数）
迭代

可以看出GMM和K-means还是有很大的相同点的。GMM中数据对高斯分量的响应度就相当于K-means中的距离计算，GMM中的根据响应度计算高斯分量参数就相当于K-means中计算分类点的位置。然后它们都通过不断迭代达到最优。不同的是：GMM模型给出的是每一个观测点由哪个高斯分量生成的概率，而K-means直接给出一个观测点属于哪一类。