EM算法、GMM应用及代码实现，如何一网打尽？

EM算法是一种迭代算法，用于含有隐变量的概率模型参数的极大似然估计。 使用EM算法的原因 　　首先举李航老师《统计学习方法》中的例子来说明为什么要用EM算法估计含有隐变量的概率模型参数。 　　假设有三枚硬币，分别记作A， B， C。这些硬币正面出现的概率分别是$\pi,p,q$。进行如下掷硬币试验：先掷硬币A，根据其结果选出硬币B或C，正面选硬币B，反面边硬币C；然后掷选出的硬币，掷硬币的结果出现正面记作1，反面记作0；独立地重复$n$次试验，观测结果为$\{y_1,y_2,...,y_n\}$。问三硬币出现正面的概率。 　　三硬币模型，也就是第二枚硬币为正面或反面的概率（$y=1$表示正面，$y=0$表示反面），或者说观测变量的概率，可以写作 $ \begin{aligned} &P(y|\pi,p,q) \\ =&\sum\limits_z P(y,z|\pi,p,q)\\ =&\sum\limits_z P(y|z,\pi,p,q)P(z|\pi,p,q)\\ =&\pi p^y(1-p)^{1-y}+(1-\pi)q^y(1-q)^{1-y} \end{aligned} $ 　　其中$z$表示硬币A的结果，也就是前面说的隐变量。为了求得参数$\pi,p,q$，我们通常会使用极大似然估计，即最大化似然函数 $ \begin{gather}\begin{aligned}&\max\limits_{\pi,p,q}\prod\limits_{i=1}^n P(y_i|\pi,p,q) \\ =&\max\limits_{\pi,p,q}\prod\limits_{i=1}^n[\pi p^{y_i}(1-p)^{1-y_i}+(1-\pi)q^{y_i}(1-q)^{1-y_i}]\\ =&\max\limits_{\pi,p,q}\sum\limits_{i=1}^n\log[\pi p^{y_i}(1-p)^{1-y_i}+(1-\pi)q^{y_i}(1-q)^{1-y_i}]\\ =&\max\limits_{\pi,p,q}L(\pi,p,q)\end{aligned}\label{}\end{gather} $ 　　分别对$\pi,p,q$求偏导并等于0，求解方程组来估计这三个参数。但是，由于它是带有隐变量的，在计算最终的概率之前有一个分支选择的过程，导致这个$\log$的内部是加和的形式，不但计算导数十分困难，待求解的方程组还不是线性方程组。当复杂度一高，解这种方程组几乎成为不可能的事。以下推导EM算法，它以迭代的方式来求解这些参数，它包含了一种“贪心”的思想。 算法导出与理解 　　对于参数为$\theta$且含有隐变量$Z$的概率模型，进行$n$次抽样。假设随机变量$Y$的观察值为$\mathcal{Y} = \{y_1,y_2,...,y_n\}$，隐变量$Z$的$m$个可能的取值为$\mathcal{Z}=\{z_1,z_2,...,z_m\}$。 　　写出似然函数： $ \begin{aligned} L(\theta) &= \sum\limits_{Y\in\mathcal{Y}}\log P(Y|\theta)\\ &=\sum\limits_{Y\in\mathcal{Y}}\log \sum\limits_{Z\in \mathcal{Z}} P(Y,Z|\theta)\\ \end{aligned} $ 　　EM算法首先初始化参数$\theta = \theta^0$，然后每一步迭代都会使似然函数增大，即$L(\theta^{k+1})\ge L(\theta^k)$。如何做到不断变大呢？考虑第$k+1$步迭代似然函数（这一步很重要！）： $ \begin{gather}\begin{aligned} L(\theta)=&\sum\limits_{Y\in \mathcal{Y}} \log\sum\limits_{Z\in \mathcal{Z}} P(Y,Z|\theta)\\ =&\sum\limits_{Y\in \mathcal{Y}} \log\sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta^k)\frac{P(Y,Z|\theta)}{P(Z|Y,\theta^k)}\\ \end{aligned}\label{} \end{gather} $ 　　至于上式的第二个等式为什么取出$P(Z|Y,\theta^k)$而不是别的，正向的原因我想不出来，马后炮原因在后面记录。