概述

EM算法是一种迭代算法，1977年由Dempster等人总结提出，用于含隐变量（Hidden variable）的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望（expectation）；M步，求极大（Maximization）。所以这一算法被称为期望极大算法（Expectation Maximization Algorithm），简称EM算法。本章首先叙述EM算法，然后讨论EM算法的收敛性；作为EM算法的应用，介绍高斯混合模型的学习；最后叙述EM算法的推广–GEM算法。

概率模型有时既含有观测变量（Observable Variable），又含有隐变量或潜在变量（Latent variable）。如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计法估计模型参数。但是，当模型含有隐变量时，就不能简单的使用这些估计方法。EM算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法。

EM算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。一般地，用Y表示观测随机变量的数据，Z表示隐随机变量的数据。Y和Z连在一起称为完全数据（complete-data），观测数据Y又称为不完全数据（incomplete-data）。假设给定观测数据Y，其概率分布是$P(Y|\theta )$，其中$\theta$是需要估计的模型参数，那么不完全数据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步，求极大化。

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 _{Z}logP(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})$$
（3）M步：求使$Q(\theta ,{\theta }^{(i)})$极大化的$\theta$，确定第i+1次迭代的参数的估计值$${\theta }^{(i+1)}=argma{x}_{\theta }Q(\theta ,{\theta }^{(i)})$$
（4）重复第2步和第4步直到收敛。

在构建具体的EM算法时，重要的是定义Q函数。每次迭代中，EM算法通过极大化Q函数来增大对数似然函数$L(\theta )$。

EM算法的收敛性

EM算法在每次迭代后，均提高观测数据的似然函数值，即
$$P(Y|{\theta }^{(i+1)})\ge P(Y|{\theta }^{(i)})$$
在一般条件下，EM算法是收敛的，但不能保证收敛到全局最优。

EM算法实现

import torch
import torch.nn as nn
import torch.nn.functional as F


class AlexNet(nn.Module):
    def __init__(self, num_classes):
        super(AlexNet, self).__init__()
        # AlexNet与LeNet一样也会同时使用卷积和池化层提取图像特征
        # 与LeNet不同的是激活函数换成了ReLU
        self.conv1 = nn.Conv2d(in_channels=3, out_channels=96, kernel_size=11, stride=4, padding=5)
        self.pool1 = nn.MaxPool2d(kernel_size=2, stride=2)

        self.conv2 = nn.Conv2d(in_channels=96, out_channels=256, kernel_size=5, stride=1, padding=2)
        self.pool2 = nn.MaxPool2d(kernel_size=2, stride=2)

        self.conv3 = nn.Conv2d(in_channels=256, out_channels=384, kernel_size=3, stride=1, padding=1)

        self.conv4 = nn.Conv2d(in_channels=384, out_channels=384, kernel_size=3, stride=1, padding=1)

        self.conv5 = nn.Conv2d(in_channels=384, out_channels=256, kernel_size=3, stride=1, padding=1)
        self.pool5 = nn.MaxPool2d(kernel_size=2, stride=2)

        self.fc1 = nn.Linear(in_features=12544, out_features=4096)
        self.drop_ratio1 = 0.5

        self.fc2 = nn.Linear(in_features=4096, out_features=4096)
        self.drop_ratio2 = 0.5

        self.fc3 = nn.Linear(in_features=4096, out_features=num_classes)

    def forward(self, x):
        in_size = x.size(0)
        x = F.relu(self.conv1(x))
        x = self.pool1(x)

        x = F.relu(self.conv2(x))
        x = self.pool2(x)

        x = F.relu(self.conv3(x))
        x = F.relu(self.conv4(x))

        x = x.view(in_size, -1)

        x = torch.relu(self.fc1(x))
        # 在全连接之后使用Dropout抑制过拟合
        x = nn.Dropout2d(x, self.drop_ratio1)

        x = torch.relu(self.fc2(x))
        # 在全连接之后使用Dropout抑制过拟合
        x = nn.Dropout2d(x, self.drop_ratio2)

        x = self.fc3(x)
        return x

K-Means与高斯混合模型

K-Means

• K-Means的目标是将样本集划分为K个簇，使得同一个簇内样本的距离尽可能的小，不同簇之间的样本距离尽可能大，即最小化每个簇中样本与质心的距离。K-means属于原型聚类（prototype-based clustering），原型聚类指聚类结构能通过一组原型刻画，而原型即为样本空间中具有代表性的点。在K-means中，这个原型就是每个簇的质心$\mu$。
• 从EM算法的观点来看，K-Means的参数就是每个簇的质心$\mu$，隐变量为每个样本的所属簇。如果事先已知每个样本属于哪个簇，则直接求平均即可得到$\mu$。但现实中不知道的情况下，则需要运用EM的思想。
• 假设要K个簇，先随机选定K个点作为质心 { μ 1 , μ 2 , μ 3 , ⋯   , μ k } \{\mu_1,\mu_2,\mu_3,\cdots,\mu_k\} {μ1​,μ2​,μ3​,⋯,μk​}

1. 固定${\mu }_k$，将样本划分到距离最近的${\mu }_k$所属的簇中，若用$r_{nk}$表示第n个样本所属的簇，则
   $$r_{nk}=\{\begin{array}{ll}1,& ifk=argmi{n}_j||X_n-{\mu }_j|{|}^2\\ 2,& otherwise\end{array}$$
2. 固定$r_{nk}$，根据上一步的划分结果重新计算每个簇的质心。由于我们的目标是最小化每个簇中样本与质心的距离，可将目标函数表示为$$J=\sum _{n=1}^N{r}_{nk}||X_n-{\mu }_k|{|}^2$$,要最小化J则对${\mu }_k$求导得$$2\sum _{n=1}^N{r}_{nk}(X_n-{\mu }_k)=0$$，则$${\mu }_k=\frac{{\sum }_n{r}_{nk}{X}_n}{{\sum }_n{r}_{nk}}$$，即簇中每个样本的均值向量。

上面两步分别更新$r_{nk}$和${\mu }_k$就对应了EM算法中的E步和M步，和EM算法一样，K-Means每一步都最小化目标函数J，因而可以保证收敛到局部最优值，但是在非凸目标函数的情况下，不能保证收敛到全局最优值。

另外，K-means对每个样本进行硬分配（Hard Assignment），即只归为一个簇，然而某些样本可能处于簇与簇的边界处，将这些样本强行分到其中一个簇可能并不能反映确信度，高斯混合模型可以进行软分配（Soft Assignment），即对每个样本划分的簇进行概率估计。

最后总结一下K-Means算法的优缺点：
优点：

1. 可解释性比较强
2. 调参的参数仅为簇数K
3. 相对于高斯混合模型而言，收敛速度快，因而常用于高斯混合模型的初始值选择，K-Means的时间复杂度为$O(N\cdot K\cdot I)$，簇数K和迭代次数I通常远小于N，所以可优化为$O(N)$，效率较高。

缺点：

• 对离群点敏感。
• K值难以事先选取，交叉验证不大合适，因为簇越多，目标函数$$\sum _{n=1}^N{r}_{nk}||X_n-{\mu }_k|{|}^2$$就越小。常采用的方法有：

一、拐点法，如下图K=3就是一个拐点。

二、加入正则化系数$\lambda$，使得$$\sum _{n=1}^N(r_{nk}||X_n-{\mu }_k|{|}^2)+\lambda K$$最小。

• 无法保证收敛到全局最优值，常使用不同的初始值进行多次实验。也可以通过K-means++算法进行优化，核心思想是选取与已有质心距离较远的点作为初始值。
• 只能发现球状的簇。
• 由于采用欧氏距离，无法直接计算类别型变量。

高斯混合模型

高斯混合模型同样用于聚类，与K-means不同的是高斯混合模型采用概率模型来表达聚类原型。
首先，高斯混合模型（Gaussian Distribution）：对随机变量x，其概率密度函数可表示为：$$N(x|\mu ,{\sigma }^2)=\frac{1}{\sqrt{(}2\pi {\sigma }^2)}exp(-\frac{(x-\mu {)}^2}{2{\sigma }^2})$$
若x为n维随机变量，则多元高斯分布（Multivariate Gaussian Distribution）为：$$N(x|\mu ,\Sigma )=\frac{1}{(2\pi {)}^{n/2}|\Sigma {|}^{1/2}}exp(-\frac{1}{2}(x-\mu {)}^T{\Sigma }^{-1}(x-\mu ))$$
其中$\mu$为n维均值向量，$\Sigma$为$n\ast n$的协方差矩阵，$|\Sigma |$为$\Sigma$的行列式。
很多时候我们发现单个高斯分布无法很好地描述数据的性质，如下图数据分为两个簇，如果使用两个高斯分布明显能更好地描述数据结构。

因此沿着这个思路就诞生了高斯混合模型（Mixture of Gaussian），本质上是K个高斯分布的线性组合，这样灵活性大增，能描述更多样的分布：$$p(x)=\sum _{k=1}^K{\pi }_{k}N(x|{\mu }_k,{\Sigma }_k)$$
其中，$0\le {\pi }_k\le 1$为混合系数，满足$$\sum _{k=1}^K{\pi }_k=1$$

下面以EM算法的角度看待高斯混合模型。EM算法是一种对含有隐变量的概率模型的极大似然估计，通常隐变量Z未知，而实际知晓的只有观测变量X。而对于高斯混合模型来说，通过采样得到的样本，却不知道每个样本来自哪个样本，这就是高斯混合模型的隐变量。

图（a）是依照完全数据的联合分布$p(x,z)$作图，每个点依照其所属于的第K个类别标记颜色，这样本的聚类属于硬分配。图（b）则不考虑隐变量Z，仅依照不完全数据x直接作图；图（c）则考虑了每个样本每个样本来自于各个类别的后验概率，这样一些在簇中间得的点的颜色是三种颜色的混合，表明这些点来自于每个簇的概率比较接近，即软分配。

高斯混合模型的优缺点：
优点：

• 相对于K-Means更具一般性，能形成各种不同大小和形状的簇。K-Means可视为高斯混合聚类中每个样本仅指派给一个混合成分的特例，且各混合成分协方差相等，均为对角矩阵${\sigma }^{2}I$。
• 仅使用少量的参数，就能较好地描述数据的特性。

缺点：

• 高斯混合模型的计算量较大且收敛较慢，因此先对样本集进行K-Means聚类，依据得到的各个簇来确定高斯混合模型的初始值，其中质心即为均值向量，协方差矩阵为每个簇中样本的协方差矩阵，混合系数为每个簇中样本占总体样本的比例。
• 分模型数量难以事先选择，但可以通过划分验证集来比较。
• 对异常点敏感。
• 数据量少时效果不好。

参考文献
1.https://blog.csdn.net/weixin_37763870/article/details/103012009
2.https://www.cnblogs.com/massquantity/p/9416109.html