混合高斯模型Gaussian Mixture Model(GMM)的EM(Expectation Maxium)求解代码

最新推荐文章于 2025-05-13 07:30:00 发布

转载最新推荐文章于 2025-05-13 07:30:00 发布 · 2.2k 阅读

文章标签：

#matlab #c

机器学习专栏收录该内容

28 篇文章

订阅专栏

本文介绍了高斯混合模型（GMM）的基本概念及其在描述复杂数据分布中的应用，并详细阐述了如何使用期望最大化（EM）算法来求解GMM中的参数估计问题。

原帖请见：http://www.matrixq.net/2011/09/10218.html

与GMM有关的matlab和c代码可以参见：

(matlab+c)https://engineering.purdue.edu/~bouman/software/cluster/

(matlab toolbox)http://www.kostaskyriakoulis.com/gmmgui.html

(matlab)http://blog.pluskid.org/?p=39

(matlab)http://www.cnblogs.com/cfantaisie/archive/2011/08/20/2147075.html

高斯分布有很多重要的性质，但是用它来描述现实中的数据的话，它还是有很多局限的。如果将这些简单的分布线性组合，就可以更好的描述实际数据的性质了，这样的模型便被称为混合模型。最常用和最流行的混合模型就是高斯混合模型GMM(Gaussian Mixture Model )。如果用足够多的高斯分布，调整其期望和协方差矩阵，以及线性组合的系数，就可以精确的表述任何连续分布。

1). 有K个高斯分布的高斯混合模型，可以表示为：

(1) $\begin{equation*} p(x)=\sum_{k=1}^{K}\pi_k \mathcal{N}(x | u_k, \Sigma_k) \end{equation*}$

其中每个高斯分布 $\mathcal{N}(x | u_k, \Sigma_k)$ 被称为一个component，参数 $\pi_k$ 被称为混合系数，为了计算和分析的方便 $\pi_k$ 通常被正则化为：

(2) $\begin{equation*} \sum_{k=1}^K \pi_k = 1 \end{equation*}$

如果把 $\pi_k=p(k)$ 看做是能取到第k个component的先验概率， $\matchcal{N}(x | u_k, \Sigma_k)=p(x|k)$ 看做是第k个component中发生x的条件概率，那么根据贝叶斯定理，发生x之后，推断其属于第k个component的后验概率为：

(3) $\begin{equation*} \gamma_k(x) \equiv p(k|x) =\frac{p(k)p(x|k)}{\sum_l p(l) p(x|l)} \end{equation*}$

2) 高斯混合模型的最大似然
想要确定 $p(X|\pi,u,\sigma)$ 中的参数 $\pi_k$ $u_k$ $\sigma_k$ ，很容易想到用最大似然法。

(4) $\begin{equation*} \begin{split} \ln p(X|\pi,u,\sigma) &= \ln \prod_{i=1}^{N}p(x_i)\\ &=\sum_{i=1}^{N} \ln \p(x_i) \\ &= \sum_{i=1}^{N} \ln \sum_{k=1}^{K}\pi_k \mathcal{N}(x_n|u_k,\Sigma_k) \end{split} \end{equation*}$

但是这个log likelihood方程并不是很好解，因为需要在对数里面求和。于是一个优雅而强大的解决此问题的方法诞生了：期望最大化法EM（Expectation-Maximization）

3) 用EM法求高斯混合模型的参数
3.1）首先求解 $u_k$ 。把 $\ln p(X|\pi,u,\Sigma)$ 看做是 $u_k$ 的方程，令式（4）的导数等于0，可得：

(5) $\begin{equation*} 0=-\sum_{n=1}^{N}\underbrace{ \frac{\pi_k \mathcal{N}(x_n | u_k,\sigma_k) }{\sum_j \pi_j \mathcal{N}(x_n | u_j, \sigma_j)}}_{\gamma (z_{nk})} \Sigma_{k} (x_n-u_k) \end{equation*}$

两边同乘以 $\Sigma_k^{-1}$ ,可得：

(6) $\begin{equation*} u_k = \frac{1}{N_k}\sum \gamma (z_{nk})x_n \end{equation*}$

其中 $N_k=\sum \gamma (z_{nk})$

3.2) 同理将 $\ln p(X|\pi,u,\Sigma)$ 看做是 $\Sigma_k$ 的方程，令式（4）的导数等于0

(7) $\begin{equation*} \Sigma_k=\frac{1}{N_k}\sum \gamma(z_nk)(x_n-u_k)(x_n-u_k)^T \end{equation*}$

3.3) 求 $\pi_k$ . 将 $\ln p(X|\pi,u,\Sigma)$ 看做是 $\pi_k$ 的方程，由于 $\sum_{k=1}^{K}\pi_k=1$ ，在这里引入拉格朗日乘子，令
$\ln p(X|\pi,u,\Sigma)+\lambda(\sum_{k=1}^{K}-1)$ 的导数等于0，可得

(8) $\begin{equation*} 0=\sum_{n=1}^{N} \frac{\pi_k \mathcal{N}(x_n | u_k,\sigma_k) } {\sum_j \pi_j \mathcal{N}(x_n | u_j, \sigma_j)} + \lambda \end{equation*}$