1 混合模型与高斯混合模型

混合模型在当前科学研究以及实际应用中都有广泛使用,作为分析复杂现象的一个灵活有力的模型,随着数据量的不断增大,混合模型相比单一模型有着明显的优势,所以它被广泛应用于医学、人脸识别、模式识别、图像处理等多个领域。本文用 EM 算法应用到混合模型的参数估计问题中。

混合模型是指随机变量 $Y$ 的概率密度函数为如下形式:
$$
P(y|\theta) = \sum_{k=1}^K \alpha_k\phi(y|\theta_k),\ \alpha_k\ge0\text{且}\sum_{k=1}^K \alpha_k=1
$$
此处的 $\theta = (\alpha_1,\cdots, \alpha_K,\theta_1,\cdots,\theta_K)​$ ,即混合模型有 $K​$ 个分支组成,每个分支的权值为 $\alpha_k​$ 。特别地,当每个分支服从的分布都是高斯分布时,则称混合分布为有 $K​$ 个分支的高斯混合模型(Gaussian miixture model,GMM)。

若为 GMM 模型,则第 $k​$ 个分支为:
$$
\phi(y|\theta_k) = \frac1{\sqrt{2\pi}\sigma_k}\exp{(-\frac{(y-\mu_k)^2}{2\sigma_k^2})},\ k=1,2,\cdots,K,\ \theta_k=(\mu_k,\sigma_k^2)
$$

一般混合模型可以由任意概率密度组合,在此仅介绍最常用的高斯混合模型。

2 高斯混合模型参数估计的 EM 算法

假设观测数据 $y_1,y_2, \cdots, y_N​$ 由高斯混合模型生产,
$$
P(y|\theta) = \sum_{k=1}^K \alpha_k\phi(y|\theta_k)
$$
其中, $\theta = (\alpha_1,\cdots, \alpha_K,\theta_1,\cdots,\theta_K)​$ 。我们用 EM 算法估计高斯混合模型的参数 $\theta​$。

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

我们的观测数据为 $y_1,y_2, \cdots, y_N$,这些数据的产生方法如下:

  1. 依据概率 $\alpha_k$ 选择第 $k$ 个高斯分布分模型 $\phi(y|\theta_k)$;
  2. 依据第 $k$ 个高斯分布分模型 $\phi(y|\theta_k)$ 生成观测数据 $y_j,j=1,2,\cdots,N$

因此,我们的观测数据 $y_1,y_2, \cdots, y_N​$ 是已知的。但是我们不知道观测数据 $y_j​$ 来自哪一个分模型,我们将 $y_j​$ 来自第 $k​$ 个分模型这一未能观测到的数据用隐变量 $\gamma_{jk}​$ 来表示,其定义如下:
$$
\gamma_{jk} = \begin{cases}
1,\text{第 $j$ 个观测来自第 $k$ 个分模型}\\
0,\text{其他}
\end{cases}
$$
其中,$j=1,2,\cdots,M;\ k=1,2,\cdots,K​$ 。

有了观测数据 $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|\theta_k)]^{\gamma_{jk}}\\
&= \prod_{k=1}^K \alpha_k^{n_k} \prod_{j=1}^N [\phi(y|\theta_k)]^{\gamma_{jk}}\\
&= \prod_{k=1}^K \alpha_k^{n_k} \prod_{j=1}^N \Big[\frac1{\sqrt{2\pi}\sigma_k}\exp{(-\frac{(y_j-\mu_k)^2}{2\sigma_k^2})}\Big]^{\gamma_{jk}}
\end{aligned}其中,$nk = \sum{j=1}^N\gamma_{jk}$ ,er
$$
其中,$n_k = \sum_{j=1}^N\gamma_{jk}​$ ,而 $ \sum_{k=1}^Kn_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 (\frac1{\sqrt{2\pi}})-\log\sigma_k - \frac1{2\sigma_k^2}(y_j-\mu_k)^2\Big] \Bigg\}
$$

2.2 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 (\frac1{\sqrt{2\pi}})-\log\sigma_k - \frac1{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 (\frac1{\sqrt{2\pi}})-\log\sigma_k - \frac1{2\sigma_k^2}(y_j-\mu_k)^2\Big] \Bigg\}
\end{aligned}
$$

这里最后一步用到了 $n_k = \sum_{j=1}^N\gamma_{jk}$

上式需要计算 $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)}{P(y_j|\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}^KP(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}
$$

由于 $\gamma_{jk}$ 是一个 0-1 随机变量,因此其期望就等于 $P(\gamma_{jk}=1|y,\theta)$,故第二个等式成立;
$\gamma_{jk}$ 只跟观测数据 $y_j$ 有关,但是观测数据 $y_j$ 跟 $\gamma_{j1},\gamma_{j2},\cdots,\gamma_{jK}$ 都有关,结合贝叶斯公式和边际分布可成立第三、四、五个等式。

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

将 $\hat \gamma_{jk} = E\gamma_{jk}$ 及 $n_k = \sum_{j=1}^N\gamma_{jk}$ 代入 $Q$ 函数中,即得:
$$
Q(\theta|\theta^{(i)}) = \sum_{k=1}^K\Bigg\{ n_k\log\alpha_k+\sum_{j=1}^N\hat\gamma_{jk}\Big[ \log (\frac1{\sqrt{2\pi}})-\log\sigma_k - \frac1{2\sigma_k^2}(y_j-\mu_k)^2 \Big] \Bigg\}
$$

2.3 确定 EM 算法的 M 步

迭代的 M 步是求函数 $Q(\theta,\theta^{(i)})$ 对 $\theta$ 的极大值,即求新一轮迭代的模型参数:
$$
\theta^{(i+1)} = arg\max_\theta Q(\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​$ 只需将 2.2 节最后得到的 $Q​$ 函数 分别对 $\hat\mu_k​$,$\hat\sigma_k^2​$ 求偏导数并令其为 0,即可得到;求 $\hat\alpha_k​$ 是在 $\sum_{k=1}^N\alpha_k=1​$ 条件下求偏导数并令其为 0 得到的。结果如下:
$$
\hat\mu_k = \frac{\sum_{j=1}^N\hat\gamma_{jk}y_j}{\sum_{j=1}^N\hat\gamma_{jk}},\ k=1,2,\cdots,K
$$

$$
\hat\sigma_k^2 = \frac{\sum_{j=1}^N\hat\gamma_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^N\hat\gamma_{jk}},\ k=1,2,\cdots,K
$$

$$
\hat\alpha_k = \frac{n_k}{N} = \frac{\sum_{j=1}^N\hat\gamma_{jk}}{N},\ k=1,2,\cdots,K
$$

重复以上计算,知道对数似然函数值不再有明显的变化为止。

2.4 高斯混合模型参数估计的 EM 算法步骤总结

输入:观测数据 $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 步:计算新一轮迭代的模型参数
$$
\hat\mu_k = \frac{\sum_{j=1}^N\hat\gamma_{jk}y_j}{\sum_{j=1}^N\hat\gamma_{jk}},\ k=1,2,\cdots,K
$$

$$
\hat\sigma_k^2 = \frac{\sum_{j=1}^N\hat\gamma_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^N\hat\gamma_{jk}},\ k=1,2,\cdots,K
$$

$$
\hat\alpha_k = \frac{n_k}{N} = \frac{\sum_{j=1}^N\hat\gamma_{jk}}{N},\ k=1,2,\cdots,K
$$

(4) 重复第 2 和第 3 步,直至收敛。

最后修改:2022 年 11 月 09 日
如果觉得我的文章对你有用,请随意赞赏