【ML】高斯混合模型GMM（生成模型）

  高斯混合模型（Gaussian Mixture Model）是一种业界广泛使用的聚类算法，该方法使用了高斯分布作为参数模型，并使用了EM算法进行训练。

高斯分布

  高斯分布是一种常用的连续变量分布的模型。若单个随机变量$x$服从均值为$\mu$，方差为${\sigma }^2$的高斯分布，记为$x\mathcal{N}(\mu ,{\sigma }^2)$，则其概率密度函数为：
$$p\left(x|\mu ,{\sigma }^2\right)=\frac{1}{{\left(2\pi {\sigma }^2\right)}^{1/2}}\exp \left\{-\frac{1}{2{\sigma }^2}(x-\mu {)}^2\right\}$$

  对于一个$d$维向量$x$，若其各元素服从均值为向量$\mu$，协方差矩阵为$\Sigma$的多元高斯分布，记为$x\mathcal{N}(\mu ,\Sigma )$，则其概率密度函数为：
$$p(\mathbf{x}|\mu ,\Sigma )=\frac{1}{(2\pi {)}^{D/2}}\frac{1}{|\Sigma {|}^{1/2}}\exp \left\{-\frac{1}{2}(\mathbf{x}-\mu {)}^{\mathbf{T}}{\Sigma }^{-1}(\mathbf{x}-\mu )\right\}$$

高斯分布函数求导

知识准备：
$$\frac{\partial }{\partial x}(x^{T}a)=\frac{\partial }{\partial x}(a^{T}x)=a$$

$$\frac{\partial }{\partial x}(AB)=\frac{\partial A}{\partial x}B+A\frac{\partial B}{\partial x}$$

$$\frac{\partial }{\partial x}(A^{-1})=-A^{-1}\frac{\partial A}{\partial x}{A}^{-1}$$

$$\frac{\partial }{\partial x}\ln |A|=(A^{-1}{)}^T$$

高斯分布的对数似然函数为：
$$\ln p(\mathbf{x}|\mu ,\Sigma )=-\frac{D}{2}\ln (2\pi )-\frac{1}{2}\ln |\Sigma |-\frac{1}{2}(\mathbf{x}-\mu {)}^{\mathbf{T}}{\Sigma }^{-1}(\mathbf{x}-\mu )$$
关于均值及协方差矩阵求导：
$$\frac{\partial }{\partial \mu }\ln p(\mathbf{x}|\mu ,\Sigma )={\Sigma }^{-1}(\mathbf{x}-\mu )$$

$$\frac{\partial }{\partial \Sigma }\ln p(\mathbf{x}|\mu ,\Sigma )=-\frac{1}{2}{\Sigma }^{-1}+\frac{1}{2}{\Sigma }^{-1}(\mathbf{x}-\mu )(\mathbf{x}-\mu {)}^T{\Sigma }^{-1}$$
则，
$$\frac{\partial }{\partial \mu }p(\mathbf{x}|\mu ,\Sigma )=p(\mathbf{x}|\mu ,\Sigma ){\Sigma }^{-1}(\mathbf{x}-\mu )$$

$$\frac{\partial }{\partial \Sigma }p(\mathbf{x}|\mu ,\Sigma )=p(\mathbf{x}|\mu ,\Sigma )\left[-\frac{1}{2}{\Sigma }^{-1}+\frac{1}{2}{\Sigma }^{-1}(\mathbf{x}-\mu )(\mathbf{x}-\mu {)}^T{\Sigma }^{-1}\right]$$

$$\frac{\partial }{\partial \Sigma }\ln |\Sigma |=({\Sigma }^{-1}{)}^T={\Sigma }^{-1}$$

$$\begin{array}{rl}\frac{\partial }{\partial {\Sigma }_{ij}}(\mathbf{x}-\mu {)}^{\mathbf{T}}{\Sigma }^{-1}(\mathbf{x}-\mu )& =(\mathbf{x}-\mu {)}^{\mathbf{T}}\frac{\partial {\Sigma }^{-1}}{\partial {\Sigma }_{ij}}(\mathbf{x}-\mu )\\ & =-(\mathbf{x}-\mu {)}^{\mathbf{T}}{\Sigma }^{-1}\frac{\partial \Sigma }{\partial {\Sigma }_{ij}}{\Sigma }^{-1}(\mathbf{x}-\mu )\\ & =-B^T\frac{\partial \Sigma }{\partial {\Sigma }_{ij}}B=-B_i{B}_j=-(B{B}^T{)}_{ij}\end{array}$$
其中$B={\Sigma }^{-1}(\mathbf{x}-\mu )$为$n\times 1$列向量，$B_i$为B的第$i$个元素；$\frac{\partial \Sigma }{\partial {\Sigma }_{ij}}$除了在第$i$行第$j$列为1外，其他位置处的值为0。
$$\frac{\partial }{\partial \Sigma }(\mathbf{x}-\mu {)}^{\mathbf{T}}{\Sigma }^{-1}(\mathbf{x}-\mu )=-B{B}^T=-{\Sigma }^{-1}(\mathbf{x}-\mu )(\mathbf{x}-\mu {)}^T{\Sigma }^{-1}$$

高斯混合模型

  高斯混合模型（Gaussian Mixed Model）指的是多个高斯分布函数的线性组合，多个高斯分布的线性叠加能拟合非常复杂的密度函数：通过足够多的高斯分布叠加，并调节它们的均值，协方差矩阵，以及线性组合的系数，可以精确地逼近任意连续密度。

  我们考虑$K$个高斯分布的线性叠加，这个高斯混合分布(Gaussian mixture distiburion)的概率密度函数为：
$$p(\mathbf{x})=\sum _{k=1}^K{\pi }_{k}p\left(\mathbf{x}|{\mu }_{\mathbf{k}},{\Sigma }_{\mathbf{k}}\right)$$
其中$p\left(\mathbf{x}|{\mu }_{\mathbf{k}},{\Sigma }_{\mathbf{k}}\right)$表示参数为${\mu }_{\mathbf{k}}$，${\Sigma }_{\mathbf{k}}$的高斯分布的概率密度。每个高斯密度函数称为混合模型的一个分模型(component)，有自己的均值${\mu }_{\mathbf{k}}$和协方差矩阵${\Sigma }_{\mathbf{k}}$。式中的参数${\pi }_k$是模型的混合系数(mixing coefficients)，满足：
$$\sum _{k=1}^K{\pi }_k=1,0\le {\pi }_k\le 0$$

  高斯混合模型可以等价的写成：
$$p(x)=\sum _{k=1}^{K}p(k)p(x|k)$$
其中，$p(k)={\pi }_k$代表选择第$k$个分模型的先验概率；$p(x|k)=p\left(\mathbf{x}|{\mu }_{\mathbf{k}},{\Sigma }_{\mathbf{k}}\right)$是$x$对$k$的条件概率密度；$p(k|x)$是指在观测到$x$后，其来自第$k$个分模型的后验概率，称为第k个分模型的响应度。

  设样本$x$由一个高斯混合模型产生，即其是从$K$个高斯分布中的某一个采样得到的，但具体是哪一个高斯分布不知道。引入一个$K$维的二值型随机变量$z=[z_1,...,z_K]$，来表示样本$x$由哪一分模型产生，$z$满足条件： z k ∈ { 0 , 1 } z_{k}\in\{0,1\} zk​∈{0,1}，且${\sum }_{k=1}^K{z}_k=1$，$z_k=1$表示样本$x$由分模型$k$抽样得到，和为1表示样本$x$只可能从一个高斯分布产生。$z$一共有$K$种可能的状态，其边缘分布由混合系数给出：$p(z_k=1)={\pi }_k$，于是$z$的分布可写作：
$$p(z)=\prod _{k=1}^K{\pi }_k^{z_k},\sum _{k=1}^K{z}_k=1$$
  来自高斯混合模型的样本$x$可以看作由以下方式产生：

1. 先以离散分布$p(z)$抽样得到变量$z$；
2. 设根据$z$的取值选择了第$k$个分模型，则以高斯分布$p(x|{\mu }_k,{\Sigma }_k)$抽样得到$x$。

  这个过程可由如下的概率图模型表示：

  给定$z$后$x$的条件概率密度为：
$$p(x|z_k=1)=p(x_{|}{\mu }_k,{\Sigma }_k)$$
或者写成：
$$p(x|z)=\prod _{k=1}^{K}p(x|{\mu }_k,{\Sigma }_k{)}^{z_k}$$
$x$的边缘分布为联合概率分布$p(x,z)$对$z$的所有可能状态求和：
$$\begin{array}{rl}p(x)& =\sum _{z}p(x,z)=\sum _{z}p(z)p(x|z)\\ & =\sum _z\left(\prod _{k=1}^K{\pi }_k^{z_k}\prod _{k=1}^{K}p{\left(x|{\mu }_k,{\Sigma }_k\right)}^{z_k}\right)\\ & =\sum _z\left(\prod _{k=1}^K{\left[{\pi }_{k}p\left(x|{\mu }_k,{\Sigma }_k\right)\right]}^{z_k}\right)\\ & =\sum _{k=1}^K{\pi }_{k}p\left(x|{\mu }_k,{\Sigma }_k\right)\end{array}$$

  观测到$x$后，其来自第$k$个分模型的后验概率(posterior responsibility)为：
$$\begin{array}{rl}p\left(z_k=1|x\right)& =\frac{p\left(z_k=1\right)p\left(x|z_k=1\right)}{p(x)}\\ & =\frac{{\pi }_{k}p\left(x|{\mu }_k,{\Sigma }_k\right)}{{\sum }_{j=1}^K{\pi }_{j}p\left(x|{\mu }_j,{\Sigma }_j\right)}\end{array}$$

  将上式定义为：$\gamma (z_k)=p(z_k=1|x)$，称为第$k$个分模型对$x$的响应度(responsibility)。

高斯混合模型的参数估计

  假设我们有数据集 X = { x 1 , . . . , x N } \mathbf{X}=\{x_1,...,x_N\} X={x1​,...,xN​}, 数据集$\mathbf{X}$由一个高斯混合模型产生。要估计这个模型的参数： π = { π 1 , . . . , π K } \pi=\{\pi_1,...,\pi_K\} π={π1​,...,πK​}， μ = { μ 1 , . . . , μ K } \mu=\{\mu_1,...,\mu_K\} μ={μ1​,...,μK​}， Σ = { Σ 1 , . . . , Σ K } \Sigma=\{\Sigma_1,...,\Sigma_K\} Σ={Σ1​,...,ΣK​}。记$x_n$对应的隐变量为$z_n=[z_{n1},...,z_{nK}]$，则：
$$p(z_n)=\prod _{k=1}^K{\pi }_k^{z_{nk}}$$

$$p(x_n|z_n)=\prod _{k=1}^{K}p(x_n|{\mu }_k,{\Sigma }_k{)}^{z_{nk}}$$

$$p(x_n)=\sum _{k=1}^K{\pi }_{k}p\left(x_n|{\mu }_k,{\Sigma }_k\right)$$

$$\gamma (z_{nk})=p\left(z_{nk}=1|x_n\right)=\frac{{\pi }_{k}p\left(x_n|{\mu }_k,{\Sigma }_k\right)}{{\sum }_{j=1}^K{\pi }_{j}p\left(x_n|{\mu }_j,{\Sigma }_j\right)}$$

最大似然估计

  样本集$X$的似然函数为：
$$p(X|\pi ,\mu ,\Sigma )=\prod _{n=1}^N[\sum _{k=1}^K{\pi }_{k}p(x_n|{\mu }_k,{\Sigma }_k)]$$
似然函数中的连乘求导比较麻烦，取自然对数将其转换成对数的和，得到对数似然函数:
$$lnp(X|\pi ,\mu \Sigma )=\sum _{n=1}^{N}ln[\sum _{k=1}^K{\pi }_{k}p(x_n|{\mu }_k,{\Sigma }_k)]$$
其中，
$$p(x_n|{\mu }_k,{\Sigma }_k)=\frac{1}{(2\pi {)}^{d/2}}\frac{1}{|{\Sigma }_k{|}^{1/2}}\exp \left\{-\frac{1}{2}(x_n-{\mu }_k{)}^{\mathbf{T}}{\Sigma }^{-1}(x_n-{\mu }_k)\right\}$$
用最大似然估计来计算高斯混合模型的参数，即求解如下的优化问题：
$$\{\begin{array}{l}{\max }_{\pi ,\mu ,\Sigma }lnp(X|\pi ,\mu ,\Sigma )\\ \\ s.t.{\sum }_{k=1}^K{\pi }_k=1\end{array}$$

采用拉格朗日乘子法来求极值，拉格朗日函数为：
$$L(\pi ,\mu ,\Sigma )=\ln p(X|\pi ,\mu ,\Sigma )+\lambda \left(\sum _{k=1}^K{\pi }_k-1\right)$$

  $L(\pi ,\mu ,\Sigma )$对${\mu }_k$求梯度：
$$\begin{array}{rl}\frac{\partial L}{\partial {\mu }_k}& =\frac{\partial }{\partial {\mu }_k}\ln p(X|\pi ,\mu ,\Sigma )\\ & =\frac{\partial }{\partial {\mu }_k}\left(\sum _{n=1}^N\ln \left[\sum _{k=1}^K{\pi }_{k}p\left(x_n|{\mu }_k,{\Sigma }_k\right)\right]\right)\\ & =\sum _{n=1}^N\frac{\partial }{\partial {\mu }_k}\ln \left[\sum _{k=1}^K{\pi }_{k}p\left(x_n|{\mu }_k,{\Sigma }_k\right)\right]\\ & =\sum _{n=1}^N\left[\frac{{\pi }_{k}p\left(x_n|{\mu }_k,{\Sigma }_k\right)}{{\sum }_{j=1}^K{\pi }_{j}p\left(x_n|{\mu }_j,{\Sigma }_j\right)}\left[{\Sigma }_k^{-1}\left(x_n-{\mu }_k\right)\right]\right]\\ & =\sum _{n=1}^N[\gamma (z_{nk}){\Sigma }_k^{-1}(x_n-{\mu }_k)]\end{array}$$
令$\frac{\partial L}{\partial {\mu }_k}=0$，即${\sum }_{n=1}^N[\gamma (z_{nk}){\Sigma }_k^{-1}(x_n-{\mu }_k)]=0$，左乘${\Sigma }_k$，得：
$$\sum _{n=1}^N\left[\gamma \left(z_{nk}\right)x_n\right]=\sum _{n=1}^N\left[\gamma \left(z_{nk}\right){\mu }_k\right]={\mu }_k\sum _{n=1}^N\gamma \left(z_{nk}\right)$$
定义$N_k={\sum }_{n=1}^N\gamma (z_{nk})$，可理解为被分配到第$k$个分模型(聚类)的“有效“的样本数。则：

$${\mu }_k=\frac{1}{N_k}\sum _{n=1}^N[\gamma (z_{nk})x_n]$$

  对${\Sigma }_k$中各元素求偏导：
$$\begin{array}{rl}\frac{\partial L}{\partial {\Sigma }_k}& =\frac{\partial }{\partial {\Sigma }_k}\ln p(X|\pi ,\mu ,\Sigma )\\ & =\frac{\partial }{\partial {\Sigma }_k}\left(\sum _{n=1}^N\ln \left[\sum _{k=1}^K{\pi }_{k}p\left(x_n|{\mu }_k,{\Sigma }_k\right)\right]\right)\\ & =\sum _{n=1}^N\frac{\partial }{\partial {\Sigma }_k}\ln \left[\sum _{k=1}^K{\pi }_{k}p\left(x_n|{\mu }_k,{\Sigma }_k\right)\right]\\ & =\sum _{n=1}^N\left[\frac{{\pi }_k}{{\sum }_{j=1}^K{\pi }_{j}p\left(x_n|{\mu }_j,{\Sigma }_j\right)}\frac{\partial }{\partial {\Sigma }_k}p\left(x_n|{\mu }_k,{\Sigma }_k\right)\right]\end{array}$$
令$\frac{\partial L}{\partial {\Sigma }_k}=0$，

$${\Sigma }_{\mathrm{k}}=\frac{1}{N_k}\sum _{n=1}^N\gamma \left(z_{nk}\right)\left(x_{\mathrm{n}}-{\mu }_{\mathrm{k}}\right){\left(x_{\mathrm{n}}-{\mu }_{\mathrm{k}}\right)}^T$$

  对${\pi }_k$求偏导并令其为0，
$$\begin{array}{rl}\frac{\partial L}{\partial {\pi }_k}& =\frac{\partial }{\partial {\pi }_k}\ln p(X|\pi ,\mu ,\Sigma )\\ & =\frac{\partial }{\partial {\pi }_k}\left[\sum _{n=1}^N\ln \left[\sum _{k=1}^K{\pi }_{k}p\left(x_{\mathrm{n}}|{\mu }_{\mathrm{k}},{\Sigma }_{\mathrm{k}}\right)\right]+\lambda \left(\sum _{k=1}^K{\pi }_k-1\right)\right]\\ & =\sum _{n=1}^N\frac{p\left(x_{\mathrm{n}}|{\mu }_{\mathrm{k}},{\Sigma }_{\mathrm{k}}\right)}{{\sum }_{j=1}^K{\pi }_{j}p\left(x_{\mathrm{n}}|{\mu }_{\mathrm{j}},{\Sigma }_{\mathrm{j}}\right)}+\lambda =0\end{array}$$
左右两边乘以${\pi }_k$，
$$\begin{array}{rl}\sum _{n=1}^N\frac{{\pi }_{k}p\left(x_n|{\mu }_k,{\Sigma }_k\right)}{{\sum }_{j=1}^K{\pi }_{j}p\left(x_n|{\mu }_j,{\Sigma }_j\right)}+\lambda {\pi }_k& =\sum _{n=1}^N\gamma \left(z_{nk}\right)+\lambda {\pi }_k\\ & =N_k+\lambda {\pi }_k=0\end{array}$$
上式对$k$求和得，
$$\sum _{k=1}^K\left(N_k+\lambda {\pi }_k\right)=\sum _{k=1}^K{N}_k+\lambda \sum _{k=1}^K{\pi }_k=\sum _{k=1}^K{N}_k+\lambda =0$$
又
$$\begin{array}{rl}\sum _{k=1}^K{N}_k& =\sum _{k=1}^K\sum _{n=1}^N\gamma \left(z_{nk}\right)=\sum _{n=1}^N\sum _{k=1}^K\gamma \left(z_{nk}\right)\\ & =\sum _{n=1}^N\sum _{k=1}^K\frac{{\pi }_{k}p\left(x|{\mu }_{\mathrm{k}},{\Sigma }_{\mathrm{k}}\right)}{{\sum }_{j=1}^K{\pi }_{j}p\left(x|{\mu }_{\mathrm{j}},{\Sigma }_{\mathrm{j}}\right)}\\ & =\sum _{n=1}^{N}1=N\end{array}$$

所以$\lambda =-N$，进而

$${\pi }_k=-\frac{N_k}{\lambda }=\frac{N_k}{N}$$

  综上，对数似然函数$lnp(X|\pi ,\mu \Sigma )$的极值点满足的条件为：
$$\{\begin{array}{l}{\mu }_{\mathrm{k}}=\frac{1}{N_k}{\sum }_{n=1}^N\left[\gamma \left(z_{nk}\right){\mathbf{x}}_{\mathrm{n}}\right]\\ \\ {\mathbf{\Sigma }}_{\mathrm{k}}=\frac{1}{N_k}{\sum }_{n=1}^N\gamma \left(z_{nk}\right)\left({\mathbf{x}}_{\mathrm{n}}-{\mu }_{\mathrm{k}}\right){\left({\mathbf{x}}_{\mathrm{n}}-{\mu }_{\mathrm{k}}\right)}^T\\ \\ {\pi }_k=\frac{N_k}{N}\end{array}$$

需要注意的是，上式并未给出高斯混合模型的解析解/闭式解，不过可使用迭代算法来计算模型参数。

求解高斯混合模型的EM算法

1. 初始化：给参数均值向量${\mu }_k$，协方差矩阵${\Sigma }_k$和混合系数${\pi }_k$赋初值；并计算对数似然函数的初值；(可以用K均值聚类(K-means clustering)算法来得到初始参数。)
2. E步(Expectation step)：用当前参数值${\mu }_k$，${\Sigma }_k$，${\pi }_k$计算后验概率(响应度)$\gamma (z_{nk})$：
   $$\gamma \left(z_{nk}\right)\leftarrow \frac{{\pi }_{k}p\left(x_n|{\mu }_k,{\Sigma }_k\right)}{{\sum }_{j=1}^K{\pi }_{j}p\left(x_n|{\mu }_j,{\Sigma }_j\right)}$$
3. M步(Maximization step)：用当前响应度$\gamma (z_{nk})$，重新估计参数：
   $$\begin{array}{rl}{\mu }_k^{\text{new }}& \leftarrow \frac{1}{N_k}\sum _{n=1}^N\left[\gamma \left(z_{nk}\right)x_n\right]\\ {\Sigma }_k^{\text{new }}& \leftarrow \frac{1}{N_k}\sum _{n=1}^N\gamma \left(z_{nk}\right)\left(x_n-{\mu }_k^{\text{new }}\right){\left(x_n-{\mu }_k^{\text{new }}\right)}^T\\ {\pi }_k^{\text{new }}& \leftarrow \frac{N_k}{N}\end{array}$$
   其中，$N_k={\sum }_{n=1}^N\gamma (z_{nk})$。
   在M步中，先计算${\mu }_k^{\text{new }}$的值，然后用${\mu }_k^{\text{new }}$来计算新的协方差矩阵${\Sigma }_k^{\text{new }}$。
4. 用新的参数${\mu }_k^{\text{new }}$、${\Sigma }_k^{\text{new }}$、${\pi }_k^{\text{new }}$计算对数似然函数：
   $$lnp(X|\pi ,\mu \Sigma )=\sum _{n=1}^{N}ln[\sum _{k=1}^K{\pi }_{k}p(x_n|{\mu }_k,{\Sigma }_k)]$$
   然后检查是否收敛。(当然，也可直接检查参数是否收敛。)
   如果收敛，则得到了模型参数；如果还没收敛，则返回第2步继续迭代。

EM算法：换个角度看高斯混合模型的极大似然估计

  在前面的极大似然估计中，我们的目标函数是不完全数据集 X = { x 1 , . . . , x N } \mathbf{X}=\{x_1,...,x_N\} X={x1​,...,xN​}的对数似然函数$lnp(X|\pi ,\mu \Sigma )$。现在我们换个角度来看：假如我们也有隐变量的观测值 Z = { z 1 , . . . , z N } \mathbf{Z}=\{z_1,...,z_N\} Z={z1​,...,zN​}，即知道$x_n$来自高斯混合中的哪一个分模型，换句话说，我们有完全数据集 { X , Z } \{\mathbf{X,Z}\} {X,Z}，那么我们可以计算完全数据的对数似然函数$lnp(X,Z|\pi ,\mu \Sigma )$。
  首先，一组完全数据 { x n , z n } \{x_n,z_n\} {xn​,zn​}的似然函数为：
$$p\left(x_n,z_n|\pi ,\mu ,\Sigma \right)=\prod _{k=1}^K{\pi }_k^{z_{nk}}p{\left(x_n|{\mu }_k,{\Sigma }_k\right)}^{z_{nk}}$$

所以，整个完全数据集 { X , Z } \{\mathbf{X,Z}\} {X,Z}的似然函数（假设每组数据独立）为：
$$p\left(\mathbf{X},\mathbf{Z}|\pi ,\mu ,\Sigma \right)=\prod _{n=1}^N\prod _{k=1}^K{\pi }_k^{z_{nk}}p{\left(x_n|{\mu }_k,{\Sigma }_k\right)}^{z_{nk}}$$

取对数得，
$$lnp\left(\mathbf{X},\mathbf{Z}|\pi ,\mu ,\Sigma \right)=\sum _{n=1}^N\sum _{k=1}^K{z}_{nk}[ln{\pi }_k+lnp\left(x_n|{\mu }_k,{\Sigma }_k\right)]$$

可以看出，和不完全数据的对数似然函数$p(X|\pi ,\mu \Sigma )={\prod }_{n=1}^N[{\sum }_{k=1}^K{\pi }_{k}p(x_n|{\mu }_k,{\Sigma }_k)]$相比，${\sum }_{k=1}^K$和$ln$换了顺序。

  然而实际情况是，我们并不知道隐变量的观测值，所以我们无法直接计算$lnp\left(\mathbf{X},\mathbf{Z}|\pi ,\mu ,\Sigma \right)$，因此，我们考虑对$lnp\left(\mathbf{X},\mathbf{Z}|\pi ,\mu ,\Sigma \right)$在隐变量 $z_{nk}$的后验分布下求期望，利用求期望的过程，将$z_{nk}$消掉。

  $lnp\left(\mathbf{X},\mathbf{Z}|\pi ,\mu ,\Sigma \right)$在隐变量 $z_{nk}$的后验分布下的期望：
E Z ln ⁡ p ( X , Z ∣ π , μ , Σ ) = E ∑ n = 1 N ∑ k = 1 K { z n k [ ln ⁡ π k + ln ⁡ p ( x n ∣ μ k , Σ k ) ] } = ∑ n = 1 N ∑ k = 1 K { E ( z n k ) [ ln ⁡ π k + ln ⁡ p ( x n ∣ μ k , Σ k ) ] } \begin{aligned} E_{Z} \ln p(X, Z | \pi, \mu, \Sigma) &=E \sum_{n=1}^{N} \sum_{k=1}^{K}\left\{z_{n k}\left[\ln \pi_{k}+\ln p\left(x_{n} | \mu_{k}, \Sigma_{k}\right)\right]\right\} \\ &=\sum_{n=1}^{N} \sum_{k=1}^{K}\left\{E\left(z_{n k}\right)\left[\ln \pi_{k}+\ln p\left(x_{n} | \mu_{k}, \Sigma_{k}\right)\right]\right\} \end{aligned} EZ​lnp(X,Z∣π,μ,Σ)​=En=1∑N​k=1∑K​{znk​[lnπk​+lnp(xn​∣μk​,Σk​)]}=n=1∑N​k=1∑K​{E(znk​)[lnπk​+lnp(xn​∣μk​,Σk​)]}​
其中，$E\left(z_{nk}\right)$表示$z_{nk}$在后验分布下的期望$E\left(z_{nk}|X,\theta \right)$。因 z n k ∈ { 0 , 1 } z_{n k}\in\{0,1\} znk​∈{0,1}，所以
$$\begin{array}{rl}E\left(z_{nk}|X,\theta \right)& =0\times p\left(z_{nk}=0|X,\theta \right)+1\times p\left(z_{nk}=1|X,\theta \right)\\ & =p\left(z_{nk}=1|X,\theta \right)\\ & =\gamma \left(z_{nk}\right)\end{array}$$
进而，
E Z ln ⁡ p ( X , Z ∣ π , μ , Σ ) = ∑ n = 1 N ∑ k = 1 K { γ ( z n k ) [ ln ⁡ π k + ln ⁡ p ( x n ∣ μ k , Σ k ) ] } E_{Z} \ln p(X, Z | \pi, \mu, \Sigma)=\sum_{n=1}^{N} \sum_{k=1}^{K}\left\{\gamma\left(z_{n k}\right)\left[\ln \pi_{k}+\ln p\left(x_{n} | \mu_{k}, \Sigma_{k}\right)\right]\right\} EZ​lnp(X,Z∣π,μ,Σ)=n=1∑N​k=1∑K​{γ(znk​)[lnπk​+lnp(xn​∣μk​,Σk​)]}
现在将$\gamma \left(z_{nk}\right)$当作常数，求参数$\pi$，$\mu$，$\Sigma$，使$E_Z\ln p(X,Z|\pi ,\mu ,\Sigma )$最大化，即
$$\{\begin{array}{l}{\max }_{\pi ,\mu ,\Sigma }{E}_Z\ln p(X,Z|\pi ,\mu ,\Sigma )\\ s.t.{\sum }_{k=1}^K{\pi }_k=1\end{array}$$
采用拉格朗日乘子法，拉格朗日函数为：
$$L_E=E_Z\ln p(X,Z|\pi ,\mu ,\Sigma )+\lambda (\sum _{k=1}^K{\pi }_k-1)$$
  先对${\mu }_k$求梯度，
$$\begin{array}{rl}\frac{\partial L_E}{\partial {\mu }_k}=0& \Rightarrow \sum _{n=1}^N\gamma \left(z_{nk}\right)\left[{\Sigma }_k^{-1}\left(x_n-{\mu }_k\right)\right]=0\\ & \Rightarrow \sum _{n=1}^N\gamma \left(z_{nk}\right)\left(x_n-{\mu }_k\right)=0\\ & \Rightarrow \sum _{n=1}^N\gamma \left(z_{nk}\right)x_n=\sum _{n=1}^N\gamma \left(z_{nk}\right){\mu }_k\end{array}$$
同之前极大化$lnp(X|\pi ,\mu ,\Sigma )$一样，可推得
$${\mu }_k=\frac{1}{N_k}\sum _{n=1}^N[\gamma (z_{nk})x_n]$$
  对${\Sigma }_k$求梯度，
$$\begin{array}{rl}\frac{\partial L_E}{\partial {\Sigma }_{\mathrm{k}}}& =\sum _{n=1}^N\left[\frac{\partial {\sum }_{k=1}^K\gamma \left(z_{nk}\right)\ln p\left(x_n|{\mu }_k,{\Sigma }_k\right)}{\partial {\Sigma }_k}\right]\\ & =\sum _{n=1}^N\left[\gamma \left(z_{nk}\right)\frac{\partial \ln p\left(x_n|{\mu }_k,{\Sigma }_k\right)}{\partial {\Sigma }_k}\right]\\ & =\sum _{n=1}^N\left[\frac{\gamma \left(z_{nk}\right)}{p\left(x_n|{\mu }_k,{\Sigma }_k\right)}\frac{\partial p\left(x_n|{\mu }_k,{\Sigma }_k\right)}{\partial {\Sigma }_k}\right]\end{array}$$
求$lnp(X|\pi ,\mu ,\Sigma )$的极值点时，有
$$\frac{\partial L}{\partial {\Sigma }_{\mathrm{k}}}=\sum _{n=1}^N\left[\frac{{\pi }_k}{{\sum }_{j=1}^K{\pi }_{j}p\left(x_n|{\mu }_j,{\Sigma }_j\right)}\frac{\partial p\left(x_n|{\mu }_k,{\Sigma }_k\right)}{\partial {\Sigma }_k}\right]$$
结合响应度$\gamma (z_{nk})$的公式，可发现两式是一致的，这意味着两种方法计算出的${\Sigma }_k$ 也相同。
  再来看${\pi }_k$，
$$\begin{array}{rl}\frac{\partial L_E}{\partial {\pi }_k}& =\sum _{n=1}^N\left[\frac{\partial {\sum }_{k=1}^K\gamma \left(z_{nk}\right)\ln {\pi }_k}{\partial {\pi }_k}\right]+\lambda \\ & =\sum _{n=1}^N\left[\gamma \left(z_{nk}\right)\frac{\partial \ln {\pi }_k}{\partial {\pi }_k}\right]+\lambda \\ & =\sum _{n=1}^N\left[\gamma \left(z_{nk}\right)\frac{1}{{\pi }_k}\right]+\lambda =\frac{1}{{\pi }_k}\sum _{n=1}^N\gamma \left(z_{nk}\right)+\lambda \\ & =\frac{1}{{\pi }_k}{N}_k+\lambda \end{array}$$
可以发现${\pi }_k$也和$lnp(X|\pi ,\mu ,\Sigma )$极值点的参数一致。

  综上，通过极大化$E_Z\ln p(X,Z|\pi ,\mu ,\Sigma )$得到的模型参数，和$lnp(X|\pi ,\mu ,\Sigma )$的极大似然估计结果是一致的。在EM算法的每一次迭代中，我们都在对$E_Z\ln p(X,Z|\pi ,\mu ,\Sigma )$求极大值。