machine learning[GMM-EM]

---

介绍下EM算法和GMM模型，先简单介绍GMM的物理意义，而后给出最直接的迭代过程；而后再介绍EM。

1 高斯混合模型

高斯分布，是统计学中的模型，其输出值表示当前输入数据样本（一维标量，多维向量）的几率。

1.1 多元高斯分布

如高斯分布-笔记所述，多元高斯函数公式为：
\[p({\bf x})=\frac{1}{(2\pi)^{\frac{d}{2}}|\Sigma|^\frac{1}{2}}exp\{-\frac{1}{2}({\bf x-\mu})^T{\Sigma}^{-1}({\bf x-\mu})\} \]
其中\({\bf x}=[x_1,x_2,...,x_d]^T\)是\(d\)维的列向量；
\({\bf \mu}=[\mu_1,\mu_2,...,\mu_d]^T\)是\(d\)维均值的列向量；
\(\Sigma\)是\(d\times d\)维的协方差矩阵；
\({\Sigma}^{-1}\)是\(\Sigma\)的逆矩阵;
\(|\Sigma|\)是\(\Sigma\)的行列式；
\((\bf x-\mu)^T\)是\((\bf x-\mu)\)的转置，且
\[\mu=E(\bf x) \]
\[\Sigma=E\{(\bf x-\bf \mu)(\bf x - \mu)^T\}\]
其中\(\mu,\Sigma\)分别是向量\(\bf x\)和矩阵\((\bf x -\mu)(\bf x -\mu)^T\)的指望。

1.1 高斯混合模型定义

如图所示，假定咱们当前的样本都是2维特征的样本，即表如今特征空间中就是在2维坐标系上有多个数据点

图1.1.1 二维空间中的数据样本下，用一个高斯分布去拟合，或者2个高斯分布模型去拟合
从左图可看出，经过一个高斯模型去拟合，总以为很违和，由于所谓的高斯分布是中间密集，而后随着距离越远则越稀疏，即 密集和稀疏程度就表示着当前出现样本点的几率,如

图1.1.2 一维高斯图形,其中x轴表示输入值，y轴表示几率值

图 1.1.3二维高斯分布呈现，乃至球面高斯分布图形。
ps: 在最后一个图形中，z轴表示当前出现样本的几率，能够看出越靠近中心，样本出现的几率越大，而越远离，则样本出现的几率越小
回头看图1.1.1，能够得出，若是是这种数据，是不合适用单一高斯分布去拟合的。

接下来咱们看看两个例子：

• 假设paperback书籍的售价是服从均值为10.00，标准差为1.00的高斯分布；而hardback 书籍的售价服从均值为17，方差为1.5的高斯分布。那么随便挑一本书，其价格是服从单一高斯分布么？
  结果确定是否。直观来讲，单一高斯分布是中间密度大，而离中间越远则密度会快速降低就像单一山峰，可是这里随便挑一本的价格倒是服从双峰的。这里分布的中心的售价大约是13，而找到13块钱的一本书的几率却远低于10或者17的周边。
• 另外一个例子：略

图1.1.4 一维数据下的高斯混合模型的两个例子
这里咱们先简单拿出高斯混合模型的公式：

如公式所述，所谓多维高斯分布，就是多个高斯分布的线性组合，也就是经过将多个高斯分布乘以对应的权重来完成样本数据的拟合。在这其中
\(0 \leq \pi_k \geq 1\)，且 \(\sum{\pi_k} =1\)
即，这里表达的物理意义是：

• 在多维特征空间中，每个样本是基于多个高斯分布而生成的，好比基于当前\(\theta\)基础上，\(p(x_i| \theta)=0.5*G1+0.3*G2+0.2*G1\)
  意思就是该样本点50%几率来自高斯1分布，30%几率来自高斯2分布，20%几率来自高斯3分布。

上面介绍的都是2个高斯分布的状况，以下图，假设咱们引入三个高斯分布

图1.1.5 三个高斯分布在2维空间中的展现
如上图所述，咱们拿到的样本集合天然都是一些在特征空间中的点，所谓用 混合高斯模型去拟合，也就是为了求得每一个样本点来自对应的高斯分布的几率（\(\pi_k\)）以及每一个高斯分布的参数（均值和方差）。能够天然推广出，只要参与的高斯分布越多，那么就能表示越多的样本点，那么只要足够多的高斯分布进行线性组合，那么就能拟合整个特征维度的空间中表示的样本。

1.2 求每一个样本属于不一样高斯分布的几率

如1.1中混合高斯的公式

咱们将其简写成

这里表示的就是每一个样本的几率来自全部基于选定第 \(k\)个高斯分布基础上该高斯分布的几率之和。那么对应的后验几率 \(p(k|\mathbf{x})\)咱们表示成 \(\gamma_k( \mathbf{x} )\):

对应的含义：

• 基于当前参数集合固定的基础上，给定一个样本，该样原本自第k个高斯分布的几率，从而完成了先混合高斯聚类，而后进行简单分类的过程。

1.3 用EM进行估计GMM的参数的直接过程

GMM模型如上述解释的，如今，咱们关注点就到了如何去用GMM拟合样本集，也就是训练过程。这其中，K是须要预先人为指定的，就和K-means同样。那么在指定了K以后，整个GMM模型的参数一共为\(\pi,\mu,\Sigma\).而一般用来训练GMM的算法叫作EM算法，也就是(expectation maximization,EM) (Dempster et al. 1977; Meng and van Dyk 1997; McLachlan and Krishnan 1997)。下面会详细介绍EM算法，而这里为了简洁，就直接给出基于EM算法的GMM求解过程：
众所周知，EM包含了两个阶段：

• E阶段：求指望函数；
• M阶段：对指望函数最大化。

在这其中，咱们须要计算最大似然估计，也就是,假设咱们有\(N\)个样本，每一个样本表示为\(\mathbf{x}_i\)，且假设有\(K\)个高斯分布，则对应的最大似然估计就是
\[p(\mathbf{X}|\pi, \mu, \Sigma)=\prod_n^{N} p(x_n|\pi, \mu, \Sigma)\]
左右取对数得：

而对于GMM模型来讲，对应的阶段为：

• 初始化阶段：先随机初始化\(\pi, \mu, \Sigma\) 。
• E阶段：使用当前的参数，基于每一个样本计算各自的高斯响应值(就是以前介绍的后验几率\(p(k|\mathbf{x})\))
  
  ps:一共\(N\)个样本，每一个样本都有\(K\)次计算，则当前E阶段一共须要计算\(N \times K\)次
• M阶段：更新参数，基于全部样本，计算各自对应的高斯参数。咱们能够想象这样一个网格，其中行表示\(N\)个样本数量，列表示\(K\)个高斯模型，则E阶段是计算网格全部的部分，而M阶段是按照列，每次计算一列中\(N\)个值的均值，要计算\(K\)个模型次数。
  
  这里
  
• 评估阶段：计算log似然
  
  并观察log似然和参数\(\pi, \mu, \Sigma\) 是否收敛，若是收敛则中止迭代，不然，接着计算E阶段和M阶段。

从上面能够看出，基于EM的GMM计算须要的计算量有些大，故而若是样本量较多，该方法一次迭代所须要的计算次数就不少。

2 EM算法

上面基于GMM给出了最终的简洁结果，这里咱们详细介绍下EM的原理及概念。值得注意的是，EM算法对初始值十分敏感，且容易陷入局部最优。

2.1 EM通用计算流程

一般在ML中，给的数据无外乎样本和标签，若是咱们采用的模型是直接基于这些计算的，那么就很容易了，好比k-means。而若是采用的模型是包含了隐含变量的，那么对于整个算法来讲，就还须要考虑到缺失的数据信息。也就是所谓的潜在变量。例如：

• 李航书里面介绍的例子：有三个硬币A，B，C。咱们先抛掷硬币A，若是朝上则抛掷硬币B，不然抛掷硬币C。而后将硬币B或者C的结果进行收集，看成观测结果，那么假设一共收集了十个[1,1,0,1,0,0,1,0,1,1]。但是咱们并不知道每次结果是来自B仍是来自C，其中包含着硬币A的未观测数据。因此咱们将观测的结果视为观测变量，而硬币A视为潜在变量。
• GMM模型：给定一个样本集合，其中每一个样本点是观测到的数据，而该数据来自哪一个高斯分布是潜在变量（虽然\(\pi_k\)表示几率，但是咱们若是取最大几率的为生成该样本的高斯分布，则能够将其视为one-hot形式）

EM算法就是一个迭代算法，经过E计算潜在变量的指望值；并在M阶段基于潜在变量指望值去计算其余模型变量的指望值。
咱们假设全部的观测数据为\(\mathbf{X}=\{\mathbf{x}_1,\mathbf{x}_2,\mathbf{x}_3,...,\mathbf{x}_N\}\)，其中，该模型包含潜在变量\(\mathbf{Z}=\{\mathbf{z}_1,\mathbf{z}_2,\mathbf{z}_3,...,\mathbf{z}_M\}\),并假设全部的模型参数为\(\theta\),那么该模型的log似然函数为：

等效于下面式子

其中 \(\mathbf{x}_i\)表示在数据集中第 \(i\)个观测样本， \(\mathbf{z}_i\)表示该观测样本下涉及的潜在变量，上式子中 \(p(\mathbf{x}_i|\theta)=\sum_{\mathbf{z}_i} p(\mathbf{x}_i,\mathbf{z}_i|\theta)\),即边缘分布。
从上面式子能够看出log没法更进一步的塞到求和操做里面，而是先求和而后在log，而后再求和，这就致使计算十分困难。这里咱们将 \(\{\mathbf{X},\mathbf{Z}\}\)称为彻底数据，而将 \(\mathbf{X}\)称为不彻底数据 。然而基于不彻底数据基础上，咱们须要进行模型的参数迭代，因此咱们先要基于 \(\theta^{old}\)计算 \(p(\mathbf{Z}|\mathbf{X},\theta^{old})\)，而后计算 \(\theta^{new}\).

这里一个重要函数\(Q(\theta,\theta^{old})\)（在M阶段须要用到），其是基于给定观测数据\(\mathbf{X}\)和当前参数\(\theta^{old}\)基础上，对未观测数据\(\mathbf{Z}\)的条件几率分布\(p(\mathbf{Z}|\mathbf{X},\theta^{old})\)的指望（李航的观点）.也能够理解成基于给定观测数据和当前参数基础上对彻底数据的指望（mlapp观点)。其中彻底数据的极大似然log函数为\(\log p(\mathbf{Z},\mathbf{X}|\theta)\)

\[\begin{eqnarray} Q(\theta,\theta^{old}) &=& E_z [\log p(\mathbf{Z},\mathbf{X}|\theta)|\mathbf{X},\theta^{old} ] \\ &=& \sum_\mathbf{Z} \log p(\mathbf{X}, \mathbf{Z} | \theta)p(\mathbf{Z} | \mathbf{X}, \theta^{old}) \\ &=& \sum_{\mathbf{Z}}p(\mathbf{Z}|\mathbf{X},\theta^{old})\log p(\mathbf{X},\mathbf{Z}|\theta) \end{eqnarray}\]
ps：上式是示意式子，具体的过程能够看基于GMM的EM过程更具体的了解这部分的变化。
通用EM流程为：

• 初始化：先初始化参数\(\theta^{old}\)
• E阶段：计算\(p(\mathbf{Z}|\mathbf{X},\theta^{old})\)
• M阶段：计算\(\theta^{new}=\underset{\theta}{argmax}Q(\theta,\theta^{old})\)
  这里\(Q(\theta,\theta^{old})=\sum_{\mathbf{Z}}p(\mathbf{Z}|\mathbf{X},\theta^{old})\log p(\mathbf{X},\mathbf{Z}|\theta)\)，所谓最大化也就是最大化\(\sum_{\mathbf{Z}}log p(\mathbf{X},\mathbf{Z}|\theta)\)这部分，即基于当前参数的彻底数据的指望
• 评估阶段：检查当前log似然函数和参数是否收敛，更换参数：\(\theta^{old}=\theta^{new}\)
  若是收敛则中止迭代，不然接着循环E阶段和M阶段

2.2 基于GMM的EM计算过程

这里咱们从新回头看GMM是如何用EM来求解的。
首先，假设观测数据\({\mathbf{x}_1,\mathbf{x}_2,...\mathbf{x}_N}\)是由高速混合模型生成的：
\[p(\mathbf{X}|\theta)=\sum_{k=1}^{K}\pi_k\phi(\mathbf{X}|\theta_k)\]
其中，\(\theta={\pi_1,\pi_2,...\pi_K;\theta_1,\theta_2,...\theta_K}\)。

1 先明确潜在变量，写出彻底数据的对数似然函数

这里假设观测数据\(\mathbf{x}_i\)是这样产生的：首先按照几率\(\pi_k\)选择第k个高斯分布，而后根据该高斯分布生成该观测数据。这时候，观测数据\(\mathbf{x}_i\)是已知的，而该观测数据来自哪一个分布是未知的，咱们假设为\(\gamma_{jk}\),表示第\(j\)个样本关于第\(k\)个高斯分布的选择：
\[\gamma_{jk}=\left\{\begin{array}{cc} 1, & 第j个样原本自第k个分模型\\ 0, & other\ values \end{array}\right.\]
其中\(\gamma_{jk}\)是个\({0,1}\)二值变量。
此时 基于每一个样本的彻底数据是\((\mathbf{x}_j,\gamma_{j1},\gamma_{j2},...,\gamma_{jK})\),其中\(j=1,2,...,N\)
则对应的彻底数据的似然函数为：
\[\begin{eqnarray}p(\mathbf{X},\mathbf{\gamma}|\theta) &=& \prod_{j=1}^{N}p(\mathbf{x}_j,\gamma_{j1},\gamma_{j2},...,\gamma_{jK}|\theta) \\ &=& \prod_{k=1}^{K}\prod_{j=1}^{N}\left [\alpha_k \phi(\mathbf{x}_j|\theta_k)\right ]^{\gamma_{jk}}\\ &=& \prod_{k=1}^K\alpha_k^{n_k}\prod_{j=1}^N\left [ \frac{1}{(2\pi)^{\frac{d}{2}}|\Sigma_k|^\frac{1}{2}}exp\{-\frac{1}{2}({\mathbf{x}_j-\mu_k})^T{\Sigma_k}^{-1}({\mathbf{x}_j-\mu_k})\}\right ]^{\gamma_{jk}} \end{eqnarray}\]
其中，\(n_k=\sum_{j=1}^N\gamma_{jk}\)，\(\sum_{k=1}^Kn_k=N\)
则，彻底数据的对数似然函数为：
\[\log p(\mathbf{X},\mathbf{\gamma}|\theta)=\sum_{k=1}^K n_k \log \alpha_k+\sum_{j=1}^N\gamma_{jk}\left [ \log\left(\frac{1}{(2\pi)^{\frac{d}{2}}}\right)-\log\left(\frac{1}{(|\Sigma_k|)^{\frac{1}{2}}}\right)-\frac{1}{2}({\mathbf{x}_j-\mu_k})^T{\Sigma_k}^{-1}({\mathbf{x}_j-\mu_k})\right ]\]

2 E阶段，肯定Q函数

\[\begin{eqnarray}Q(\theta,\theta^{old}) &=& E\left [\log p(\mathbf{X},\mathbf{\gamma}|\theta)|\mathbf{X},\theta^{old}\right ]\\ &=& E\left (\sum_{k=1}^K n_k \log \alpha_k+\sum_{j=1}^N\gamma_{jk}\left [\log\left(\frac{1}{(2\pi)^{\frac{d}{2}}}\right)-\log\left(\frac{1}{(|\Sigma_k|)^{\frac{1}{2}}}\right)-\frac{1}{2}({\mathbf{x}_j-\mu_k})^T{\Sigma_k}^{-1}({\mathbf{x}_j-\mu_k})\right ] \right ) \\ &=& \sum_{k=1}^K\left( \sum_{j=1}^N(E\gamma_{jk})\log\alpha_k+\sum_{j=1}^N(E\gamma_{jk})\left [ \log\left(\frac{1}{(2\pi)^{\frac{d}{2}}}\right)-\log\left(\frac{1}{(|\Sigma_k|)^{\frac{1}{2}}}\right)-\frac{1}{2}({\mathbf{x}_j-\mu_k})^T{\Sigma_k}^{-1}({\mathbf{x}_j-\mu_k})\right ]\right)\\ &=& \sum_{k=1}^K\sum_{j=1}^N(E\gamma_{jk})\left(\log\alpha_k+\left [ \log\left(\frac{1}{(2\pi)^{\frac{d}{2}}}\right)-\log\left(\frac{1}{(|\Sigma_k|)^{\frac{1}{2}}}\right)-\frac{1}{2}({\mathbf{x}_j-\mu_k})^T{\Sigma_k}^{-1}({\mathbf{x}_j-\mu_k})\right ]\right) \\ \end{eqnarray}\]
ps:最后一个式子，就等于2.1部分的表现形式\(\sum_{\mathbf{Z}}p(\mathbf{Z}|\mathbf{X},\theta^{old})\log p(\mathbf{X},\mathbf{Z}|\theta)\)
能够看到，对整个彻底数据的指望，最后就须要求解关于潜在变量的指望，这里只须要计算\(E(\gamma_{jk}|\mathbf{X},\theta)\),这里记为\(\hat \gamma_{jk}\)
则
\[\begin{eqnarray}\hat \gamma_{jk} &=& E(\gamma_{jk}|\mathbf{X},\theta^{old})=p(\gamma_{jk}=1|\mathbf{X},\theta^{old})\\ &=& \frac{p(\gamma_{jk}=1,\mathbf{x}_j|\theta^{old})}{\sum_{k=1}^Kp(\gamma_{jk}=1,\mathbf{x}_j|\theta^{old})}\\ &=& \frac{p(\gamma_{jk}=1,\mathbf{x}_j|\theta)p(\gamma_{jk}=1|\theta)}{\sum_{k=1}^Kp(\gamma_{jk}=1,\mathbf{x}_j|\theta)p(\gamma_{jk}=1|\theta)}\\ &=& \frac{\alpha_k\phi(\mathbf{x}_j|\theta_k)}{\sum_{k=1}^K\alpha_k\phi(\mathbf{x}_j|\theta_k)} \end{eqnarray}\]
其中\(j=1,2,...,N;k=1,2,...K\)。
将\(\hat \gamma_{jk}=E(\gamma_{jk}|\mathbf{X},\theta)\)和\(n_k=\sum_{j=1}^NE\hat \gamma_{jk}=\sum_{j=1}^N\gamma_{jk}\)带入得：
\[Q(\theta,\theta^{old})=\sum_{k=1}^K n_k \log \alpha_k+\sum_{j=1}^N\hat \gamma_{jk}\left [ \log\left(\frac{1}{(2\pi)^{\frac{d}{2}}}\right)-\log\left(\frac{1}{(|\Sigma_k|)^{\frac{1}{2}}}\right)-\frac{1}{2}({\mathbf{x}_j-\mu_k})^T{\Sigma_k}^{-1}({\mathbf{x}_j-\mu_k})\right ] \\ \]
ps:值得注意的是，上面式子中的中括号部分，能够替换成任意其余分模型，不必定是高斯模型

3 肯定M阶段

M阶段就是为了求\(Q(\theta,\theta^{old})\)的关于\(\theta\)的极大值，即：
\[\theta^{new}=\underset{\theta}{argmax}Q(\theta,\theta^{old})\]
这里用\(\hat \mu_k\),\(\hat \sigma_k^2\),\(\hat \alpha_k\),\(k=1,2,3,...,K\)表示\(\theta^{new}\)中的各参数，求\(\hat \mu_k\),\(\hat \sigma_k^2\)只须要将\(Q(\theta,\theta^{old})\)求偏导数，而后等于0，求出其值便可，而在求\(\hat \alpha_k\)时，由于有约束条件\(\sum_{k=1}^K\alpha_k=1\)，因此须要用拉格朗日乘子法求偏导。
\[\begin{eqnarray}\frac{dQ}{d\mu_k} &=& \sum_{j=1}^K\hat \gamma_{jk}\left[\Sigma_k(\mathbf{x}_j-\mu_k)\right]\\ \end{eqnarray}\]
令其等于0，并左右乘以\(\Sigma_k^{-1}\)(咱们假设该协方差矩阵是非奇异的)。则\(\sum_{j=1}^N\hat \gamma_{jk}\mathbf{x}_j=\sum_{j=1}^N\hat \gamma_{jk}\mu_k=E\mu_k\sum_{j=1}^N\hat \gamma_{jk}\),因此
\[\hat \mu_k=\frac{\sum_{j=1}^N\hat \gamma_{jk}\mathbf{x}_j}{\sum_{j=1}^N\hat \gamma_{jk}}\]
同理：
\[\frac{dQ}{d\Sigma_k}=\]

2.3 EM的证实过程 李航

参考文献：

• 高斯分布-笔记
• Gaussian mixture models and the EM algorithm
• gmm
• Pattern.Recognition.and.Machine.Learning,.Christopher.M..Bishop,.Springer,.2006
• Machine learning a probabilistic perspective
• 李航统计学习方法
• 高斯混合模型（GMM）及其EM算法的理解
• matrix_cookbook
• Matrix_calculus