【混合模型】一文搞懂混合模型、高斯混合模型（GMM）、EM算法

1. 聚类问题（Clustering Problem）

• 聚类是将一组对象根据它们的相似性进行分组，使得同一组内的对象相比于其他组的对象更加相似。

• 聚类是一种无监督学习任务，即在没有事先给定类别标签的情况下，根据相似性对对象进行分组。这与分类（有监督学习）问题不同，分类是基于已知的类别标签进行分类的。

• 聚类问题分为硬聚类（hard clustering）和软聚类（soft clustering），前者将每个对象明确地进行划分，如K-means算法，而后者每个对象可以以某种概率或权重属于多个簇。现实中，研究者可能更感兴趣的是一些具有连续性的聚类分配（软聚类），因为硬聚类的离散结果可能不符合实际需求。这就需要构建一种概率模型来处理具有多个簇的复杂数据。

2. 混合模型（Mixture Model）

• 混合模型不止使用一个单一参数模型，而是包含$K$个参数模型的概率密度函数。$p(y_i|\pi ,\theta )={\sum }_{j=1}^K{\pi }_{j}p(y_i|{\theta }_j)$。每一个样本$y_i$都从这样一个加权的概率分布中抽样而来，所以 $y$ 的边际似然函数为： ∏ i = 1 n p ( y i ∣ π , θ ) = ∏ i = 1 n { ∑ j = 1 K π j p ( y i ∣ θ j ) } \prod_{i=1}^np(y_i|\pmb{\pi},\pmb{\theta})=\prod_{i=1}^n\{\sum_{j=1}^K\pi_jp(y_i|\theta_j)\} ∏i=1n​p(yi​∣π,θ)=∏i=1n​{∑j=1K​πj​p(yi​∣θj​)}。

【举个例子】：下图为一元混合高斯模型示意图：

• 模型的概率密度函数（pdf）：$p(y_i|\pi ,\theta )$ 是$K$个独立的参数模型的加权组合。其中，$\{({\pi }_j,{\theta }_j),j=1,\cdots ,K\}$ 分别是$K$个模型的参数
• 混合模型引入隐变量 $z$ 来表示数据属于哪一个聚类。

3. 混合模型的生成过程

• Step 1：我们先生成分类隐变量 z ∈ { 1 , ⋯   , K } z\in\{1,\cdots,K\} z∈{1,⋯,K}，并且 $\mathrm{Pr}(z=j)={\pi }_j$

• Step 2：然后我们生成 $y|z\sim p(y|{\theta }_z)$。于是我们可以得到 $y,z$ 的联合分布：$p(y,z)=p(y|z)\mathrm{Pr}(z)$ 于是 $y$ 的边际似然函数：$p(y)={\sum }_{j=1}^{K}p(y,z=j)={\sum }_{j=1}^{K}p(y|z=j)\mathrm{Pr}(z=j)$

• $\mathrm{Pr}(z=j)={\pi }_j$ 表示在没有观察数据点时，数据点属于簇 $j$ 的先验概率

• 观察到数据点 $y$ 后，我们希望通过贝叶斯定理来计算数据点属于簇 $j$ 的后验概率：$\mathrm{Pr}(z=j|y)=\frac{p(y|z=j)\mathrm{Pr}(z=j)}{p(y)}=\frac{p(y|z=j)\mathrm{Pr}(z=j)}{{\sum }_{l=1}^{K}p(y|z=l)\mathrm{Pr}(z=l)}=\frac{p(y|z=j){\pi }_j}{{\sum }_{l=1}^{K}p(y|z=l){\pi }_l}$

4. 高斯混合模型（GMM）

4.1 已知参数 $\pi ,\mu ,{\sigma }^2$ 更新 $z$

• $p(y|\pi ,\mu ,{\sigma }^2)={\sum }_{j=1}^K{\pi }_j\mathcal{N}(y|{\mu }_j,{\sigma }_j^2)$：边际似然，其中 N ( y ∣ μ j , σ j 2 ) = 1 2 π σ 2 exp ⁡ { − ( y − μ ) 2 2 σ 2 } \mathcal{N}(y|\mu_j,\sigma_j^2)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\{-\frac{(y-\mu)^2}{2\sigma^2}\} N(y∣μj​,σj2​)=2πσ2 ​1​exp{−2σ2(y−μ)2​}

• $\mathrm{Pr}(z=j)={\pi }_j$：先验，从第 $j$ 个混合组件中抽样的概率

  $p(y|z=j)=\mathcal{N}(y|{\mu }_j,{\sigma }_j^2)$：条件似然

  $\Rightarrow p(z=j|y)=\frac{p(y|z=j)\mathrm{Pr}(z=j)}{p(y)}$

4.2 已知 $z$ 更新参数 $\pi ,\mu ,{\sigma }^2$

• 现在我们有观测数据 y = { y 1 , ⋯   , y n } \pmb{y}=\{y_1,\cdots,y_n\} y={y1​,⋯,yn​}，似然函数 L ( π , μ , σ 2 ∣ y ) = ∏ i = 1 n p ( y ∣ π , μ , σ 2 ) = ∏ i = 1 n { ∑ j = 1 K π j N ( y ∣ μ j , σ j 2 ) } L(\pmb{\pi},\pmb{\mu},\pmb{\sigma^2}|\pmb{y})=\prod_{i=1}^np(y|\pmb{\pi},\pmb{\mu},\pmb{\sigma^2})=\prod_{i=1}^n\{\sum_{j=1}^K\pi_j\mathcal{N}(y|\mu_j,\sigma_j^2)\} L(π,μ,σ2∣y)=∏i=1n​p(y∣π,μ,σ2)=∏i=1n​{∑j=1K​πj​N(y∣μj​,σj2​)} 。

  注意到，似然函数是一个求和累乘的形式，我们很难通过最大似然估计推导出${\hat{\pi }}_{MLE},{\hat{\mu }}_{MLE},{\hat{\sigma }}_{MLE}$

• 假设我们知道每一个观测变量 $y_i$ 的类别 $z_i$，$(y_i,z_i)$ 的联合分布有如下形式：$p(y_i,z_i|\pi ,\mu ,{\sigma }^2)=\mathrm{Pr}(z_i)p(y_i|z_i)={\pi }_{z_i}\mathcal{N}(y_i|{\mu }_{z_i},{\sigma }_{z_i}^2)$

  联合对数似然：
  log ⁡ p ( y , z ) = ∑ i = 1 n log ⁡ p ( y i , z i ) = ∑ i = 1 n { log ⁡ π z i + log ⁡ N ( y i ∣ μ z i , σ z i 2 ) } = ∑ i = 1 n { log ⁡ π z i − ( y i − μ z i ) 2 2 σ z i 2 − 1 2 log ⁡ ( 2 π σ z i 2 ) } \begin{aligned} \log p(\pmb{y},\pmb{z})&=\sum_{i=1}^n\log p(y_i,z_i)\\ &=\sum_{i=1}^n\{\log \pi_{z_i}+\log \mathcal{N}(y_i|\mu_{z_i},\sigma_{z_i}^2)\}\\ &=\sum_{i=1}^n\{\log \pi_{z_i}-\frac{(y_i-\mu_{z_i})^2}{2\sigma_{z_i}^2}-\frac{1}{2}\log (2\pi\sigma_{z_i}^2)\} \end{aligned} logp(y,z)​=i=1∑n​logp(yi​,zi​)=i=1∑n​{logπzi​​+logN(yi​∣μzi​​,σzi​2​)}=i=1∑n​{logπzi​​−2σzi​2​(yi​−μzi​​)2​−21​log(2πσzi​2​)}​

• 记 $\mathrm{T}=[{\tau }_{ij}{]}_{i=1,j=1}^{n,K}\in {\mathbb{R}}^{n\times K}$，${\tau }_{ij}$ 表示第 $i$ 个观测数据属于第 $j$ 个高斯组件的概率， 定义指示函数：
  $$I(z_i=j)=\{\begin{array}{l}1,\arg \underset{l}{\max }{z}_{il}=j\\ 0,\arg \underset{l}{\max }{z}_{il}\ne j\end{array}$$

• 则：

$$\begin{array}{rl}\log p(y,z)& =\sum _{i=1}^n\sum _{j=1}^{K}I(z_i=j)\cdot (\log {\pi }_j+\log \mathcal{N}(y_i|{\mu }_j,{\sigma }_j^2))\end{array}$$

• 通过最大联合对数似然，我们可以得到：

$$\begin{array}{rl}{\pi }_j^{\ast }& =\frac{{\sum }_{i=1}^{n}I(z_i=j)}{n}\\ {\mu }_j^{\ast }& =\frac{{\sum }_{i=1}^n{y}_i\cdot I(z_i=j)}{{\sum }_{i=1}^{n}I(z_i=j)}\\ {\sigma }_j^{2\ast }& =\frac{{\sum }_{i=1}^n(y_i-{\mu }_j^{\ast }{)}^2\cdot I(z_i=j)}{{\sum }_{i=1}^{n}I(z_i=j)}\end{array}$$

• ${\pi }_j^{\ast }$：$n$ 个观测数据中，属于第 $j$ 个高斯组件的样本个数

  ${\mu }_j^{\ast }$：第 $j$ 个高斯组件的样本均值

  ${\sigma }_j^{2\ast }$：第 $j$ 个高斯组件的样本方差

4.3 EM算法

由4.1和4.2我们知道，如果我们知道混合高斯模型的参数，我们能够使用贝叶斯定理轻松地进行软聚类；如果我们提前知道每一个样本的类别标签 $z$，我们能够轻松地估计混合高斯模型的参数 $\pi ,\mu ,{\sigma }^2$。但问题是，我们啥也不知道！然而，我们也许可以构建一个迭代算法，来反复进行4.1和4.2中的过程知道算法收敛。这就是大名鼎鼎的Expectation-Maximization算法，简称EM算法。

• EM算法:
  • 初始化参数 $(\pi ,\mu ,{\sigma }^2)$
  • E-step：计算软聚类 ${\tau }_{ij}=\mathrm{Pr}(z_i=j|\pi ,\mu ,{\sigma }^2)=\frac{{\pi }_j\mathcal{N}(y_i|{\mu }_j,{\sigma }_j^2)}{{\sum }_{l=1}^K{\pi }_l\mathcal{N}(y_i|{\mu }_l,{\sigma }_l^2)}$
  • M-step：更新GMM参数：
    • ${\pi }_j^{(new)}=\frac{{\sum }_{i=1}^n{\tau }_{ij}}{n}$
    • ${\mu }_j^{(new)}=\frac{{\sum }_{i=1}^n{y}_i\cdot {\tau }_{ij}}{{\sum }_{i=1}^n{\tau }_{ij}}$
    • $({\sigma }_j^2{)}^{(new)}=\frac{{\sum }_{i=1}^n{\tau }_{ij}(y_i-{\mu }_j^{(new)}{)}^2}{{\sum }_{i=1}^n{\tau }_{ij}}$
  • 重复E-step和M-step直到收敛。

EM算法是否收敛到MLE的证明（略）。