EM 算法 & 推广¶

EM 算法背景简介
EM 算法详解
高斯混合模型
Python 代码实现
EM 算法相关面经和其他资料
小结

EM 算法背景简介¶

EM 算法是在 1977 年由 Dempster 等人提出，是一种迭代算法。其用于还有隐变量（hidden variable）的概率模型参数的极大似然估计，或极大后验概率估计。EM 算法的每次迭代由两步组成：

E 步：求期望（expectation）
M 步：求极大（maximization）

所以这一算法称为期望极大算法。

概率模型有时既还有观测变量（observable variable），又含有隐变量或潜在变量（latent variable）。如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计，或者贝叶斯估计来估计参数。但是，当模型中含有隐变量时，这些简单的估计方法便不再适用。

EM 算法详解¶

E 和 M 是什么？

为了更形象的引入 EM 算法，不妨先介绍一个使用 EM 算法的例子。

假设有 3 枚硬币，分别记作 A，B，C。这些硬币正面出现的概率分别为 $\pi$，$p$ 和 $q$。进行如下掷硬币试验：

先抛硬币 A，根据其结果选出 B 或者 C，正面选 B，反面选 C。
然后抛选出的硬币，抛硬币的结果，出现正面记作 1，反面记作 0
独立重复 n 次试验（n=10），我们得到：

$$1,1,0,1,0,0,1,0,1,1$$

假设只能观测到抛硬币的结果，不能观测抛硬币的过程。问如何估计三个硬币正面出现的概率。

解：

三硬币模型可写作：

$$ \begin{aligned} P(y|\theta) &= \sum_{z}P(y,z|\theta) \\ &= \sum_{z}P(z|\theta)P(y|z,\theta) \\ &= \pi p^y(1-p)^{1-y} + (1-\pi)q^y(1-q)^{1-y} \end{aligned} $$

这里，随机变量 y 是观测变量，表示一次试验观测的结果是 1 或 0；随机变量 $z$ 是隐变量，表示未观测到的抛硬币 $A$ 的结果；$\theta = (\pi, p, q)$ 是模型参数。这一模型是以上数据的生成模型。注意，随机变量 $y$ 的数据可以观测，随机变量 $z$ 的数据不可观测。

将观测数据表示为 $Y = (Y_1,Y_2,\cdots,Y_n)^T$，未观测数据表示为 $Z = (Z_1,Z_2,\cdots,Z_n)^T$ 则观测数据的似然函数为

$$P(Y|\theta) = \sum_{Z}P(Z|\theta)P(Y|Z,\theta) $$

也就是

$$P(Y|\theta) = \prod_{j=1}^{n}[\pi p^{y_j}(1-p)^{1-y_j} + (1-\pi)q^{y_j}(1-q)^{1-y_j}]$$

那么我们要求参数 $\theta = (\pi,p,q)$ 的极大似然估计，也就是

$$\hat \theta = arg \max_{\theta} log P(Y|\theta)$$

这个问题没有解析解，只能用迭代的方法求解。所以 EM 算法首先初始化参数，记作 $\theta^{(0)} = (\pi^{(0)}, p^{(0)}, q^{(0)})$，然后通过迭代来计算参数的估计值，知道收敛。第 $i$ 次迭代参数的估计值为 $\theta^{(i)} = (\pi^{(i)}, p^{(i)}, q^{(i)})$。EM 算法的第 $i+1$ 次迭代如下：

E 步：计算在模型参数 $\pi^{(i)}, p^{(i)}, q^{(i)}$ 下观测数据 $y_j$ 来自抛硬币 $B$ 的概率
M 步：计算模型参数的新估计值

注意：

EM 算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。

一般，用 $Y$ 表示观测随机变量的数据，$Z$ 表示隐随机变量的数据。$Y$ 和 $Z$ 连在一起称为完全数据（complete-data），观测数据 $Y$ 又称为不完全数据。EM 算法通过迭代求 $L(\theta) = log P(Y|\theta)$ 的极大似然估计。每次迭代包含两步：E 步，求期望；M 步，求极大化。

EM 算法

输入：观测变量数据 $Y$，隐变量数据 $Z$，联合分布 $P(Y,Z|\theta)$，条件分布 $P(Z|Y,\theta)$；

输出：模型参数 $\theta$

选择参数的初值 $\theta^{(0)}$，开始迭代；
E 步：记 $\theta^{(i)}$ 为第 $i$ 次迭代参数 $\theta$ 的估计值，在第 $i+1$ 次迭代的 $E$ 步，计算

$$ \begin{aligned} Q(\theta,\theta^{(i)}) &= E_z[log P(Y,Z|\theta)|Y,\theta^{(i)}] \\ &= \sum_{Z}log P(Y, Z|\theta) P(Z|Y,\theta^{(i)}) \end{aligned}$$ 这里，$P(Z|Y,\theta^{(i)})$ 是在给定观测数据 $Y$ 和当前的参数估计 $\theta^{(i)}$ 下隐变量数据 $Z$ 的条件概率分布；
M 步：求使 $Q(\theta,\theta^{(i)})$ 极大化的 $\theta$，确定第 $i+1$ 次迭代的参数的估计值 $\theta^{(i+1)}$ $$\theta^{(i+1)} = arg \max_{\theta} Q(\theta,\theta^{(i)})$$
重复第 2 步和第 3 步，直到收敛。

那么，肯定会有人问这个 $Q(\theta,\theta^{(i)})$ 是什么，它就是 EM 算法的核心，称为 $Q$ 函数，下面来着重讲解它。

$Q$ 函数¶

$Q$ 函数定义：

完全数据的对数似然函数 $log P(Y,Z|\theta)$ 关于在给定观测数据 $Y$ 和当前参数 $\theta^{(i)}$ 下对未观测数据 $Z$ 的条件概率分布 $P(Z|Y,\theta^{(i)})$ 的期望称为 $Q$ 函数，即

$$Q(\theta,\theta^{(i)}) = E_Z[log P(Y,Z|\theta)|Y,\theta^{(i)}]$$

那么要注意以下几个问题：

参数的初值可以任意选择，但需要注意 EM 算法对初值是非常敏感的。
E 步求 $Q(\theta,\theta^{(i)})$。$Q$ 函数式中 $Z$ 是未观测数据，$Y$ 是观测数据。这里要注意上角标 $(i)$ 的含义。
- $\theta$ 表示要极大化的参数
- $\theta^{(i)}$ 表示参数的当前估计值每次迭代其实在求 $Q$ 函数及其极大！
M 步求 $Q(\theta,\theta^{(i)})$ 的极大化，得到 $\theta^{(i+1)}$，完成一次迭代 $\theta^{(i)} \rightarrow \theta^{(i+1)}$
那么什么时候停止迭代呢？一般是对较小的正数 $\epsilon_1$，$\epsilon_2$，若满足 $$||\theta^{(i+1)} - \theta^{(i)}|| < \epsilon_1$$ 或 $$||Q(\theta^{(i+1)} ,\theta^{(i)}) - Q(\theta^{(i)} ,\theta^{(i)})|| < \epsilon_2$$ 则停止迭代。

算法导出¶

EM 算法已经叙述完毕了，但是在理解一个算法的时候，要更深层次理解其原理。为什么 EM 算法能够近似实现对观测数据的极大似然估计呢？

对于一个含有隐变量的概率模型，目标是极大化观测数据（也就是不完全数据）$Y$ 关于参数 $\theta$ 的对数似然函数，也就是极大化

$$ \begin{aligned} L(\theta) &= log P(Y|\theta) = log \sum_{Z} P(Y,Z|\theta) \\ &= log \Big(\sum_{Z}P(Y|Z,\theta)P(Z|\theta) \Big) \end{aligned}$$

事实上，EM 算法是通过迭代逐步近似极大化 $L(\theta)$ 的。假设在第 $i$ 次迭代后 $\theta$ 的估计值是 $\theta^{(i)}$。我们希望新估计值 $\theta$ 能使 $L(\theta)$ 增加，也就是 $L(\theta) > L(\theta^{(i)})$，并逐步达到极大值。为此，考虑两者的差：

$$L(\theta) - L(\theta^{(i)}) = log \Bigg( \sum_{Z}P(Y|Z,\theta)P(Z|\theta) \Bigg) - log P(Y|\theta^{(i)})$$

我们利用 Jensen 不等式

$$log\sum_{j} \lambda_jy_j \geq \sum_{j} \lambda_j \log y_j, \ \lambda_j \geq 0, \ \sum_{j} \lambda_j = 1$$

得到其下界：

$$ \begin{aligned} L(\theta) - L(\theta^{(i)}) &= \log \Bigg( \sum_{Z} P(Z|Y,\theta^{(i)}) \frac{P(Y|Z,\theta)P(Z|\theta) }{P(Z|Y,\theta^{(i)})}\Bigg) - log P(Y|\theta^{(i)}) \\ &\geq \sum_{Z} P(Z|Y,\theta^{(i)}) \log \frac{P(Y|Z,\theta)P(Z|\theta) }{P(Z|Y,\theta^{(i)})} - log P(Y|\theta^{(i)}) \\ &= \sum_{Z} P(Z|Y,\theta^{(i)}) \log \frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})} \end{aligned} $$

令

$$B(\theta,\theta^{(i)}) \hat= L(\theta^{(i)}) + \sum_{Z} P(Z|Y,\theta^{(i)}) \log \frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})} \tag{1}$$

则

$$L(\theta) \geq B(\theta,\theta^{(i)})$$

所以函数 $B(\theta,\theta^{(i)})$ 是 $L(\theta)$ 的一个下界，且由 $(1)$ 式可知，

$$L(\theta^{(i)}) = B(\theta^{(i)}, \theta^{(i)})$$

因此，任何可以使 $B(\theta,\theta^{(i)})$ 增大的 $\theta$，也可以使 $L(\theta)$ 增大。为了使 $L(\theta)$ 有尽可能大的增长，我们应选择 $\theta^{(i+1)}$ 使 $B(\theta,\theta^{(i)})$ 达到极大，也就是

$$\theta^{(i+1)} = arg \max_{\theta} B(\theta,\theta^{(i)})$$

那么要求 $\theta^{(i+1)}$ 的表达式，就要提出对于 $\theta$ 而言是常数的项，从而得到：

$$ \begin{aligned} \theta^{(i+1)} &= arg \max_{\theta} \Bigg(L(\theta^{(i)}) + \sum_{Z} P(Z|Y,\theta^{(i)}) \log \frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})} \Bigg) \\ &= arg \max_{\theta} \Bigg( \sum_{Z} P(Z|Y,\theta^{(i)}) \log P(Y|Z,\theta)P(Z|\theta) \Bigg) \\ &= arg \max_{\theta} \Bigg( \sum_{Z} P(Z|Y,\theta^{(i)}) \log P(Y,Z|\theta)\Bigg) \\ &= arg \max_{\theta} Q(\theta,\theta^{(i)}) \end{aligned} $$

上式相当于 EM 算法的一次迭代，也就是求 $Q$ 函数及其极大化。EM 算法是通过不断求解下界的极大化逼近求解对数似然函数的极大化的算法。

下图给出了算法更直观的解释，上方的曲线为$L(\theta)$，下方曲线为 $B(\theta,\theta^{(i)})$。其为对数似然函数 $L(\theta)$ 的下界。两个函数在点 $\theta = \theta^{(i)}$ 处相等。EM 算法要找到下个点 $\theta^{(i+1)}$ 使得函数 $B(\theta,\theta^{(i)})$ 极大化，同时随着 $B(\theta,\theta^{(i)})$ 的增加，我们保证对数似然函数 $L(\theta)$ 也是在迭代过程中不断增加的。在迭代过程中，EM 算法不能保证找到全局最优，也就说明了为什么初始值的选取会很大程度上影响最后的结果！

EM 算法也可以用于生成模型的无监督学习。生成模型由联合概率分布 $P(X,Y)$ 表示，可认为无监督学习训练数据是联合概率分布产生的数据。

收敛性分析¶

我们已经了解了 EM 算法的基本原理和导出，知道它是一种近似计算含有隐变量概率模型的极大似然估计方法。那么我们怎么判断算法得到的估计序列是是否收敛呢？如果收敛，是否收敛到全局最大或者局部最大？

关于 EM 算法收敛性有两个定理，我们不妨了解一下：

设 $P(Y|\theta)$ 为观测数据的似然函数，$\theta^{(i)}$ 为 EM 算法得到的参数估计序列，$P(Y|\theta^{(i)})$ 为对应的似然函数序列，则 $P(Y|\theta^{(i)})$ 是单调递增的，也就是 $$P(Y|\theta^{(i+1)}) \geq P(Y|\theta^{(i)})$$
设 $L(\theta) = \log P(Y|\theta)$ 为观测数据的对数似然函数，$\theta^{(i)}$ 为 EM 算法得到的参数估计序列，$L(\theta^{(i)})$ 为对应的对数似然函数序列。
- 如果 $P(Y|\theta)$ 有上界，则 $L(\theta^{(i)}) = \log P(Y|\theta^{(i)})$ 收敛到某一值 $L^*$
- 在函数 $Q(\theta,\theta')$ 与 $L(\theta)$ 满足一定条件下，由 EM 算法得到的参数估计序列 $\theta^{(i)}$ 的收敛值 $\theta^*$ 是 $L(\theta)$ 的稳定点。

高斯混合模型¶

EM 算法的一个非常重要的应用是高斯混合模型的参数估计。高斯混合模型应用广泛，那么先来介绍一下高斯混合模型吧。

高斯混合模型是指具有如下形式的概率分布模型：

$$P(y|\theta) = \sum_{k=1}^{K}\alpha_k\phi(y|\theta_k)$$

其中，$\alpha_k$ 是系数，$\alpha_k \geq 0$，$\sum_{k=1}^{K}\alpha_k = 1$；$\phi(y|\theta_k)$ 是高斯分布密度，$\theta_k = (\mu_k, \sigma_k^2)$，

$$\phi(y|\theta_k) = \frac{1}{\sqrt{2\pi}\sigma_k}exp \Bigg(-\frac{(y-\mu_k)^2}{2\sigma_k^2} \Bigg)$$

称为第 $k$ 个模型。

下面再介绍高斯混合模型参数估计的 EM 算法！

假设观测数据 $y_1,y_2,\cdots,y_N$ 由高斯混合模型生成，

$$P(y|\theta) = \sum_{k=1}^{K}\alpha_k\phi(y|\theta_k)$$

其中，$\theta = (\alpha_1, \alpha_2, \cdots, \alpha_K; \theta_1,\theta_2,\cdots,\theta_K)$。我们用 EM 算法估计高斯混合模型的参数 $\theta$。

明确隐变量，写出完全数据的对数似然函数

可以设想观测数据 $y_j$，$j=1,2,\cdots,N$，是这样产生的：首先依概率 $\alpha_k$ 选择第 $k$ 个高斯分布分模型 $\phi(y|\theta_k)$，然后依第 $k$ 个分模型的概率分布 $\phi(y|\theta_k)$ 生成观测数据 $y_j$。这时观测数据 $y_j$，$j=1,2,\cdots,N$ 是已知的；反映观测数据 $y_j$ 来自第 $k$ 个分模型的数据是未知的，$k=1,2,\cdots,K$，以隐变量 $\gamma_{jk}$ 表示，其定义如下：

$$ \gamma_{jk}=\begin{cases} 1, \ 第\ j\ 个观测来自第 \ k\ 个分模型\\ 0, \ 否则\end{cases} $$ $$j=1,2,\cdots,N; \ k = 1,2,\cdots,K$$ $\gamma_{jk}$ 是 0-1 随机变量。有了观测数据 $y_j$ 及未观测数据 $\gamma_{jk}$，那么完全数据是 $$(y_j,\gamma_{j1}, \gamma_{j2},\cdots,\gamma_{jK}), \ j = 1,2,\cdots,N$$ 于是，可以写出完全数据的似然函数：

$$ \begin{aligned} P(y,\gamma|\theta) &= \prod_{j=1}^{N}P(y_j,\gamma_{j1}, \gamma_{j2},\cdots,\gamma_{jK}|\theta) \\ &= \prod_{k=1}^{K}\prod_{j=1}^{N}[\alpha_k\phi(y_j|\theta_k)]^{\gamma_{jk}} \\ &= \prod_{k=1}^{K} \alpha_k^{n_k} \prod_{j=1}^{N}[\phi(y_j|\theta_k)]^{\gamma_{jk}} \\ &= \prod_{k=1}^{K} \alpha_k^{n_k} \prod_{j=1}^{N} \Bigg[\frac{1}{\sqrt{2\pi}\sigma_k}exp \Bigg(-\frac{(y_j-\mu_k)^2}{2\sigma_k^2} \Bigg) \Bigg]^{\gamma_{jk}} \end{aligned} $$

式中，$n_k = \sum_{j=1}^{N}\gamma_{jk}$，$\sum_{k=1}^{K}n_k = N$。那么，完全数据的对数似然函数为

$$\log P(y,\gamma|\theta) = \sum_{k=1}^{K} \Bigg \{ n_k \log \alpha_k + \sum_{j=1}^{N} \gamma_{jk} \Big[ log \Big( \frac{1}{\sqrt{2 \pi}} \Big) - \log \sigma_k - \frac{1}{2 \sigma_k^2} (y_j - \mu_k)^2 \Big] \Bigg \}$$
EM 算法的 E 步：确定 $Q$ 函数

$$ \begin{aligned} Q(\theta,\theta^{(i)}) &= E[\log P(y,\gamma|\theta)|y,\theta^{(i)}] \\ &= E\Bigg \{\sum_{k=1}^{K} \Bigg \{ n_k \log \alpha_k + \sum_{j=1}^{N} \gamma_{jk} \Big[ log \Big( \frac{1}{\sqrt{2 \pi}} \Big) - \log \sigma_k - \frac{1}{2 \sigma_k^2} (y_j - \mu_k)^2 \Big] \Bigg \} \Bigg \} \\ &= E\Bigg \{\sum_{k=1}^{K} \Bigg \{ \sum_{j=1}^{N}\gamma_{jk} \log \alpha_k + \sum_{j=1}^{N} \gamma_{jk} \Big[ log \Big( \frac{1}{\sqrt{2 \pi}} \Big) - \log \sigma_k - \frac{1}{2 \sigma_k^2} (y_j - \mu_k)^2 \Big] \Bigg \} \Bigg \} \\ &= \sum_{k=1}^{K} \Bigg \{ \sum_{j=1}^{N} (E(\gamma_{jk})) \log \alpha_k + \sum_{j=1}^{N} (E(\gamma_{jk})) \Big[ log \Big( \frac{1}{\sqrt{2 \pi}} \Big) - \log \sigma_k - \frac{1}{2 \sigma_k^2} (y_j - \mu_k)^2 \Big] \Bigg \} \Bigg \} \end{aligned} $$ 这里需要计算 $E(\gamma_{jk} | y,\theta)$，记为 $\hat \gamma_{jk}$。

$$ \begin{aligned} \hat \gamma_{jk} &= E(\gamma_{jk} | y,\theta) = P(\gamma_{jk} = 1|y,\theta) \\ &= \frac{P(\gamma_{jk} = 1,y_j|\theta)}{\sum_{k=1}^{K}P( \gamma_{jk} = 1,y_j | \theta)} \\ &= \frac{P(y_j|\gamma_{jk} = 1,\theta)P(\gamma_{jk}=1 | \theta)}{\sum_{k=1}^{k}P(y_j|\gamma_{jk} = 1,\theta)P(\gamma_{jk}=1 | \theta)} \\ &= \frac{\alpha_k \phi(y_j|\theta_k)}{\sum_{k=1}^{K}\alpha_k \phi(y_j|\theta_k)} , \ j = 1,2,\cdots, N; \ k = 1,2, \cdots, K \end{aligned} $$

$\hat \gamma_{jk}$ 是在当前模型参数下第 $j$ 个观测数据来自第 $k$ 个分模型的概率，称为分模型 $k$ 对观测数据 $y_j$ 的响应度。所以得到：

$$Q(\theta,\theta^{(i)}) = \sum_{k=1}^{K} \Bigg \{ n_k \log \alpha_k + \sum_{j=1}^{N} \hat \gamma_{jk} \Big[ log \Big( \frac{1}{\sqrt{2 \pi}} \Big) - \log \sigma_k - \frac{1}{2 \sigma_k^2} (y_j - \mu_k)^2 \Big] \Bigg \} \tag{1}$$

确定 EM 算法的 M 步 迭代的 M 步是求函数 $Q(\theta,\theta^{(i)})$ 对 $\theta$ 的极大值，即求新一轮迭代的模型参数： $$\theta^{(i+1)} = \arg \max_{\theta}(\theta,\theta^{(i)})$$

用 $\hat \mu_k$，$\hat \sigma_k^2$ 以及 $\hat \alpha_k$，$k=1,2,\cdots,K$，表示 $\theta^{(i+1)}$ 的各参数。求 $\hat \mu_k$，$\hat \sigma_k^2$ 只需将式（1）分别对 $\hat \mu_k$，$\hat \sigma_k^2$ 求偏导数并令其为 0，即可得到；求 $\hat \alpha_k$ 是在 $\sum_{k=1}^{K} \alpha_k = 1$ 条件下（拉格朗日乘子）求偏导数并令其为 0 得到的。结果如下：

重复以上计算，直到对数似然函数值不再有明显变化为止。

算法：

输入：观测数据 $y_1,y_2,\cdots,y_N$，高斯混合模型；

输出：高斯混合模型参数

（1）取参数的初始值开始迭代

（2）E 步：依据当前模型参数，计算分模型 $k$ 对观测数据 $y_j$ 的响应度 $$\hat \gamma_{jk} = \frac{\alpha_k \phi(y_j|\theta_k)}{\sum_{k=1}^{K}\alpha_k \phi(y_j|\theta_k)} , \ j = 1,2,\cdots, N; \ k = 1,2, \cdots, K$$

（3）M 步：计算新一轮迭代的模型参数

（4）重复第 2 步和第 3 步，直到收敛

Python 代码实现¶

Numpy¶

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from scipy.stats import multivariate_normal
plt.style.use('seaborn')

# 生成数据
def generate_X(true_Mu, true_Var):
    # 第一簇的数据
    num1, mu1, var1 = 400, true_Mu[0], true_Var[0]
    X1 = np.random.multivariate_normal(mu1, np.diag(var1), num1)
    # 第二簇的数据
    num2, mu2, var2 = 600, true_Mu[1], true_Var[1]
    X2 = np.random.multivariate_normal(mu2, np.diag(var2), num2)
    # 第三簇的数据
    num3, mu3, var3 = 1000, true_Mu[2], true_Var[2]
    X3 = np.random.multivariate_normal(mu3, np.diag(var3), num3)
    # 合并在一起
    X = np.vstack((X1, X2, X3))
    # 显示数据
    plt.figure(figsize=(10, 8))
    plt.axis([-10, 15, -5, 15])
    plt.scatter(X1[:, 0], X1[:, 1], s=5)
    plt.scatter(X2[:, 0], X2[:, 1], s=5)
    plt.scatter(X3[:, 0], X3[:, 1], s=5)
    plt.show()
    return X

#### E step ####

# 更新W也就是隐变量gamma_jk --> 第j个观测量属于第k类的概率
# gamma_jk = alpha*phi / sum alpha*phi
# alpha*phi = pdfs
def update_W(X, Mu, Var, Pi):
    n_points, n_clusters = len(X), len(Pi)
    pdfs = np.zeros(((n_points, n_clusters)))
    for k in range(n_clusters):
        pdfs[:, k] = Pi[k] * multivariate_normal.pdf(X, Mu[k], np.diag(Var[k]))
    W = pdfs / pdfs.sum(axis=1).reshape(-1, 1)
    return W


# 更新pi也就是alpha,也就是先验概率，每一类的概率
def update_Pi(W):
    Pi = W.sum(axis=0) / 2000 # W.sum()
    return Pi

#### M step ####
    
# 计算log似然函数 Q函数
def logLH(X, Pi, Mu, Var):
    n_points, n_clusters = len(X), len(Pi)
    pdfs = np.zeros(((n_points, n_clusters)))
    for k in range(n_clusters):
        pdfs[:, k] = Pi[k] * multivariate_normal.pdf(X, Mu[k], np.diag(Var[k]))
    return np.mean(np.log(pdfs.sum(axis=1)))
# loglh = E(log P(y|theta))


# 画出聚类图像
def plot_clusters(X, Mu, Var, Mu_true=None, Var_true=None):
    colors = ['b', 'g', 'r']
    n_clusters = len(Mu)
    plt.figure(figsize=(10, 8))
    plt.axis([-10, 15, -5, 15])
    plt.scatter(X[:, 0], X[:, 1], s=5)
    ax = plt.gca()
    for i in range(n_clusters):
        plot_args = {'fc': 'None', 'lw': 2, 'edgecolor': colors[i], 'ls': ':'}
        ellipse = Ellipse(Mu[i], 3 * Var[i][0], 3 * Var[i][1], **plot_args)
        ax.add_patch(ellipse)
    if (Mu_true is not None) & (Var_true is not None):
        for i in range(n_clusters):
            plot_args = {'fc': 'None', 'lw': 2, 'edgecolor': colors[i], 'alpha': 0.5}
            ellipse = Ellipse(Mu_true[i], 3 * Var_true[i][0], 3 * Var_true[i][1], **plot_args)
            ax.add_patch(ellipse)         
    plt.show()


# 更新Mu 也就是对于数据关于隐变量的加权平均
def update_Mu(X, W):
    n_clusters = W.shape[1]
    Mu = np.zeros((n_clusters, 2))
    for i in range(n_clusters):
        Mu[i] = np.average(X, axis=0, weights=W[:, i])
    return Mu


# 更新Var
def update_Var(X, Mu, W):
    n_clusters = W.shape[1]
    Var = np.zeros((n_clusters, 2))
    for i in range(n_clusters):
        Var[i] = np.average((X - Mu[i]) ** 2, axis=0, weights=W[:, i])
    return Var


if __name__ == '__main__':
    # 生成数据
    true_Mu = [[0.5, 0.5], [5.5, 2.5], [1, 7]]
    true_Var = [[1, 3], [2, 2], [6, 2]]
    X = generate_X(true_Mu, true_Var)
    # 初始化
    n_clusters = 3
    n_points = len(X)
    Mu = [[0, -1], [6, 0], [0, 9]]
    Var = [[1, 1], [1, 1], [1, 1]]
    Pi = [1 / n_clusters] * 3
    W = np.ones((n_points, n_clusters)) / n_clusters 
    Pi = W.sum(axis=0) / W.sum()
    # 迭代
    loglh = []
    for i in range(5):
        plot_clusters(X, Mu, Var, true_Mu, true_Var)
        loglh.append(logLH(X, Pi, Mu, Var))
        W = update_W(X, Mu, Var, Pi)
        Pi = update_Pi(W)
        Mu = update_Mu(X, W)
        print('log-likehood:%.3f'%loglh[-1])
        Var = update_Var(X, Mu, W)

log-likehood:-8.166

log-likehood:-4.719

log-likehood:-4.715

log-likehood:-4.715

log-likehood:-4.714

Sklearn¶

import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
import numpy as np
from matplotlib.patches import Ellipse

#产生实验数据
from sklearn.datasets.samples_generator import make_blobs
X, y_true = make_blobs(n_samples=700, centers=4,
                       cluster_std=0.5, random_state=2019)
X = X[:, ::-1] #方便画图

from sklearn.mixture import GaussianMixture as GMM
gmm = GMM(n_components=4).fit(X) #指定聚类中心个数为4
labels = gmm.predict(X)
plt.scatter(X[:, 0], X[:, 1], c=labels, s=50, cmap='viridis')

probs = gmm.predict_proba(X)
print(probs[:10].round(2))

#给定的位置和协方差画一个椭圆
def draw_ellipse(position, covariance, ax=None, **kwargs):
    ax = ax or plt.gca()
    #将协方差转换为主轴
    if covariance.shape == (2, 2):
        U, s, Vt = np.linalg.svd(covariance)
        angle = np.degrees(np.arctan2(U[1, 0], U[0, 0]))
        width, height = 2 * np.sqrt(s)
    else:
        angle = 0
        width, height = 2 * np.sqrt(covariance)

    #画出椭圆
    for nsig in range(1, 4):
        ax.add_patch(Ellipse(position, nsig * width, nsig * height,
                             angle, **kwargs))
#画图
def plot_gmm(gmm, X, label=True, ax=None):
    ax = ax or plt.gca()
    labels = gmm.fit(X).predict(X)
    if label:
        ax.scatter(X[:, 0], X[:, 1], c=labels, s=4, cmap='viridis', zorder=2)
    else:
        ax.scatter(X[:, 0], X[:, 1], s=4, zorder=2)
    ax.axis('equal')
    w_factor = 0.2 / gmm.weights_.max()
    for pos, covar, w in zip(gmm.means_, gmm.covariances_  , gmm.weights_):
        draw_ellipse(pos, covar, alpha=w * w_factor)
        
        
# 椭圆拟合数据       
rng = np.random.RandomState(13)
X_stretched = np.dot(X, rng.randn(2, 2))
gmm = GMM(n_components=4, covariance_type='full', random_state=42)
plot_gmm(gmm, X_stretched)
plt.show()

# 特殊分布形式
from sklearn.datasets import make_moons
Xmoon, y = make_moons(100, noise=.04, random_state=0)
plt.scatter(Xmoon[:, 0], Xmoon[:, 1]);
plt.show()

# 2个高斯分布拟合
gmm2 = GMM(n_components=2, covariance_type='full', random_state=0)
plot_gmm(gmm2, Xmoon)
plt.show()

# 10个高斯分布拟合
gmm10 = GMM(n_components=10, covariance_type='full', random_state=0)
plot_gmm(gmm10, Xmoon, label=False)
plt.show()

# 200个高斯分布拟合
Xnew = gmm10.sample(200)[0]
plt.scatter(Xnew[:, 0], Xnew[:, 1])
plt.show()

# 最优类别个数确定(k)
n_components = np.arange(1, 21)
models = [GMM(n, covariance_type='full', random_state=0).fit(Xmoon)
          for n in n_components]

plt.plot(n_components, [m.bic(Xmoon) for m in models], label='BIC')
plt.plot(n_components, [m.aic(Xmoon) for m in models], label='AIC')
plt.legend(loc='best')
plt.xlabel('n_components')

C:\Users\gzjgz\anaconda3\lib\site-packages\sklearn\utils\deprecation.py:144: FutureWarning: The sklearn.datasets.samples_generator module is  deprecated in version 0.22 and will be removed in version 0.24. The corresponding classes / functions should instead be imported from sklearn.datasets. Anything that cannot be imported from sklearn.datasets is now part of the private API.
  warnings.warn(message, FutureWarning)

[[0.   0.   1.   0.  ]
 [0.96 0.   0.   0.04]
 [1.   0.   0.   0.  ]
 [1.   0.   0.   0.  ]
 [0.   0.   1.   0.  ]
 [0.   1.   0.   0.  ]
 [1.   0.   0.   0.  ]
 [1.   0.   0.   0.  ]
 [0.   0.   0.   1.  ]
 [0.   1.   0.   0.  ]]

Text(0.5, 0, 'n_components')

EM 算法相关面经和其他资料¶

面经¶

EM 算法的流程是什么？ E 是为了求出类别对隐变量的响应程度，也就是观测数据属于某一类的概率，M 是为了最大化对数似然函数，并且更新模型参数，通过 EM 的不断迭代，直到参数和似然函数收敛。

GMM 为什么要用 EM 算法求解： 高斯混合模型也就是用多个高斯分布密度函数的加权平均来分析数据分布。为了估计其参数，一般选用最大似然估计，但是得到的函数没有解析解，所以使用 EM 算法，EM算法解决这个的思路是使用启发式的迭代方法，既然我们无法直接求出模型分布参数，那么我们可以先猜想隐含数据（EM算法的E步），接着基于观察数据和猜测的隐含数据一起来极大化对数似然，求解我们的模型参数（EM算法的M步)。由于我们之前的隐藏数据是猜测的，所以此时得到的模型参数一般还不是我们想要的结果。不过没关系，我们基于当前得到的模型参数，继续猜测隐含数据（EM算法的E步），然后继续极大化对数似然，求解我们的模型参数（EM算法的M步)。以此类推，不断的迭代下去，直到模型分布参数基本无变化，算法收敛，找到合适的模型参数。

K-Means 和 GMM 的区别
- kmeans是GMM的一种特例
- kmeans使用了集群中心的平均自，针对扭曲混合型的数据难以发挥作用，而GMM假设数据点服从高斯分布，假设限制较少
- kmeans是硬分配（要么属于，要么不属于），GMM是软分配（计算数据点属于某一类的概率）
- GMM比K-Means在群协方面更灵活。由于标准偏差参数，集群可以采取任何椭圆形状，而不是限于圆形
- GMM计算较为复杂，效果更好，一般用看means做粗略估计，再将得到的结果（质心）带入GMM做进一步的细致估计

其他资料¶

小结¶

EM 算法是还有隐变量的概率模型极大似然估计的迭代算法。含有隐变量的概率模型的数据表示为 $P(Y,Z|\theta)$。包括 E 步，求期望，M 步，求极大。
一般条件下，EM 算法是收敛的，但不能保证全局最优。
EM 算法很广泛，GMM（高斯混合模型）和 HMM（隐藏马尔科夫模型）的参数估计都是 EM 的重要应用，HMM 接下来会介绍。