【机器学习】EM算法在高斯混合模型学习中的应用

前言

EM算法，此博客介绍了$EM$算法相关理论知识，看本篇博客前先熟悉$EM$算法。

本篇博客打算先从单个高斯分布说起，然后推广到多个高斯混合起来，最后给出高斯混合模型参数求解过程。

单个高斯分布

假如我们有一些数据，这些数据来自同一个高斯分布（独立同分布），那个我们如何根据这些数据估计出这个高斯分布的参数呢？我们知道只要知道高斯分布的参数$\theta=\{\mu,\sigma^2\}$就能确定此高斯分布。

从上图中，我们要想知道数据是来自哪个高斯分布，我们就要知道高斯分布的参数，直观上可以认定数据来自参数为 $\theta_1$的高斯分布，然而这毕竟是我们直观上的，那么我们应该如何根据数据估计高斯分布的参数呢？

假设数据为$X = \{x_1,x_2,...,x_N\}，x_i$ 独立同分布 $p(X|\theta)$，其中$\theta=\{\mu,\sigma^2\}$;

由贝叶斯公式知：

$$p(\theta|X) \propto p(X|\theta)p(\theta)$$
$p(\theta|X)$为后验概率， $p(X|\theta)$为似然度， $p(\theta)$为先验概率。

关于贝叶斯公式的，可参考我之前的博客，里面有提到：链接

不加上先验概率，就是极大似然估计；
加上先验概率，就是极大后验概率估计
这里我们只介绍极大似然估计，极大后验概率估计类似

（1）写出对数似然函数

$$\mathcal{L}(\theta|X) = log\ [p(X|\theta)] \\ = \sum_{i=1}^N log\ p(x_i|\theta)\\ = \sum_{i=1}^N log\ [\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(x_i - \mu)^2}{2\sigma^2})]$$

（2）求极大似然估计

分别对$\mu,\sigma$求偏导，并令为0：

$$\frac{\partial\mathcal{L}(\theta|X)}{\partial\mu} = \frac{\partial(\sum_{i=1}^N log\ [\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(x_i - \mu)^2}{2\sigma^2})])}{\partial\mu} \\ =\frac{\partial(\sum_{i=1}^Nlog[exp(-\frac{(x_i - \mu)^2}{2\sigma^2})])}{\partial \mu} \\ =\frac{\partial(\sum_{i=1}^N-\frac{(x_i - \mu)^2}{2\sigma^2})}{\partial \mu} = - \sum_{i=1}^N\frac{(x_i - \mu)}{\sigma^2} = 0$$
由此得到 $\mu$的似然估计为：
$$\mu_{MLE} = \frac{1}{N} \sum_{i=1}^Nx_i$$
$$\frac{\partial\mathcal{L}(\theta|X)}{\partial\sigma} = \frac{\partial(\sum_{i=1}^N log\ [\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(x_i - \mu_{MLE})^2}{2\sigma^2})])}{\partial\sigma}\\ = \frac{\partial\sum_{i=1}^Nlog\frac{1}{\sqrt{2\pi}\sigma}}{\partial\sigma}+\frac{\partial(\sum_{i=1}^Nlog[exp(-\frac{(x_i - \mu_{MLE})^2}{2\sigma^2})])}{\partial\sigma} \\ =\sum_{i=1}^N-(\frac{1}{\sqrt{2\pi}\sigma})\sqrt{2\pi} + \frac{\partial(\sum_{i=1}^N-\frac{(x_i - \mu_{MLE})^2}{2\sigma^2})}{\partial \sigma} \\ =-N+\frac{\sum_{i=1}^N(x_i - \mu_{MLE})^2}{\sigma^2} = 0$$
其中用到 $\frac{d \ log(\frac{1}{x})}{dx} = -\frac{1}{x}$
得到 $\sigma^2$的似然估计为：
$$\sigma^2 = \frac{\sum_{i=1}^N(x_i - \mu_{MLE})^2}{N}$$

综上我们可以表述为：

θ=argmaxθ[∑i=1Nlog N(xi|μ,σ2)] θ = a r g max θ [ ∑ i = 1 N l o g   N ( x i | μ , σ 2 ) ]

$$\theta = arg \max_{\theta} [\sum_{i=1}^{N}log \ \mathcal{N}(x_i|\mu,\sigma^2)]$$
$$\mu \rightarrow \frac{\partial \mathcal{L}(\mu, \sigma^2|X)}{\partial \mu} \\ \sigma^2 \rightarrow \frac{\partial \mathcal{L}(\mu, \sigma^2|X)}{\partial \sigma}$$

以上问题是只有一个高斯分布的，如果不止一个高斯分布呢？

高斯混合模型

同样，我们看下图：

上图是一个高斯混合模型，此处只画了两个高斯分布，可以是多个高斯分布。

如果我们知道每一个数据属于哪一个高斯分布，就会很容易求解，但是我们不可能都知道的，这时随便一个数据点，我们应该如何判断它是哪个高斯分布产生的呢？

这里我们引入$\alpha_k$表示属于第$k$个高斯分布的权重，并满足$\sum_{k=1}^{K} \alpha_k = 1$

这样我们可以得到：

$$p(X|\theta) = \sum_{k=1}^{K} \alpha_k \mathcal{N}(\mu_k, \sigma^2_k) \\ \sum_{i=k}^{K} \alpha_k= 1$$

其中$\mathcal{N}(\mu_k, \sigma^2_k)$是高斯分布，$\alpha_k$是系数

给出定义：

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

$$p(x|\theta) = \sum_{k=1}^{K}\alpha_k \phi(x|\theta_k)$$
其中， $\alpha_k$是系数， $\alpha_k\geqslant 0 \sum_{k=1}^{K} \alpha_k= 1$； $\phi(x|\theta_k)$是高斯分布密度， $\theta_k = (\mu_k, \sigma^2_k)$，
$$\phi(x|\theta_k) = \frac{1}{\sqrt{2\pi}\sigma_k}exp(-\frac{(x - \mu_k)^2}{2\sigma_k^2})$$
称为第 $k$个分模型。

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

假设观测数据$x_1,x_2,...,x_N$由高斯混合模型生成，

$$p(x|\theta) = \sum_{k=1}^K\alpha_k \phi(x|\theta_k)$$
其中， $\theta = (\alpha_1,\alpha_2,..,\alpha_k;\theta_1,\theta_2,..,\theta_k)$，我们用 $EM$算法估计高斯混合模型的参数 $\theta$.

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

我们设想数据是这样产生的：

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

$$\gamma_{jk}=\begin{cases} 1,\quad 第j个观测来自第k个分模型 \\ 0, \quad 否则 \\ \end{cases} \\ j = 1,2,...,N; k = 1,2,...,K$$

有了观测数据$x_j$及未观测数据$\gamma_{jk}$，那么完全数据是：$(x_j,\gamma_{j1},\gamma_{j2},...,\gamma_{jK})， j = 1,2,...,N$

于是，可以写出完全数据的似然函数：

$$p(x,\gamma|\theta) = \prod_{j=1}^Np(x_j,\gamma_{j1},\gamma_{j2},...,\gamma_{jK}|\theta)\\ = \prod_{k=1}^K\prod_{j=1}^N[\alpha_k \phi(x_j|\theta_k)]^{\gamma_{jk}} \\ = \prod_{k=1}^K\alpha_k^{n_k}\prod_{j=1}^N[\phi(x_j|\theta_k)]^{\gamma_{jk}}\\ = \prod_{k=1}^K\alpha_k^{n_k}\prod_{j=1}^N[\frac{1}{\sqrt{2\pi}\sigma_k}exp(-\frac{(x_j - \mu_k)^2}{2\sigma_k^2})]^{\gamma_{jk}}$$
式中： $n_k = \sum_{j=1}^N \gamma_{jk},\sum_{k=1}^Kn_k = N$

那么，对数似然函数为：

$$log \ p(x,\gamma|\theta) = \sum_{k=1}^K \Big\{ 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^2_k}(x_j-\mu_k)^2\Big]\Big\}$$

2. $EM$算法$E$步：确定Q函数

$$Q(\theta, \theta^{(i)}) = E[log \ p(x,\gamma | \theta)|x,\theta^{(i)}] \\ = E\Big\{ \sum_{k=1}^K \Big\{ 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^2_k}(x_j-\mu_k)^2\Big]\Big\}\Big\}\\ = \sum_{k=1}^K\Big\{ \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^2_k}(x_j-\mu_k)^2\Big]\Big\}$$
这里需要计算$E(\gamma_{jk} | x, \theta)$，记为$\hat{\gamma_{jk}}$

γjk^=E(γjk|x,θ)=p(γjk=1|x,θ)=p(γjk=1,xj|θ)∑Kk=1p(γjk=1,xj|θ)=p(xj|γjk=1,θ)p(γjk=1|θ)∑Kk=1p(xj|γjk=1,θ)p(γjk=1|θ)=αkϕ(xj|θk)∑Kk=1αkϕ(xj|θk),j=1,2,...,N;k=1,2,...,K γ j k ^ = E ( γ j k | x , θ ) = p ( γ j k = 1 | x , θ ) = p ( γ j k = 1 , x j | θ ) ∑ k = 1 K p ( γ j k = 1 , x j | θ ) = p ( x j | γ j k = 1 , θ ) p ( γ j k = 1 | θ ) ∑ k = 1 K p ( x j | γ j k = 1 , θ ) p ( γ j k = 1 | θ ) = α k ϕ ( x j | θ k ) ∑ k = 1 K α k ϕ ( x j | θ k ) , j = 1 , 2 , . . . , N ; k = 1 , 2 , . . . , K

$$\hat{\gamma_{jk}} = E(\gamma_{jk} | x, \theta) = p(\gamma_{jk} = 1 |x,\theta) \\ = \frac{p(\gamma_{jk} = 1, x_j |\theta)}{\sum_{k=1}^Kp(\gamma_{jk}=1,x_j | \theta)} \\ = \frac{p(x_j|\gamma_{jk} = 1,\theta)p(\gamma_{jk}=1|\theta)}{\sum_{k=1}^Kp(x_j|\gamma_{jk} = 1,\theta)p(\gamma_{jk}=1|\theta)}\\ = \frac{\alpha_k \phi(x_j|\theta_k)}{\sum_{k=1}^K \alpha_k \phi(x_j|\theta_k)}, j = 1,2,...,N; k = 1,2,...,K$$

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

将$\hat{\gamma_{jk}} = E\gamma_{jk} 及 n_k = \sum_{j=1}^NE\gamma_{jk}$代入之前式子即得：

$$Q(\theta, \theta^{(i)}) =\sum_{k=1}^K \Big \{n_klog \ \alpha_k + \sum_{k=1}^K \hat{\gamma_{jk}} \Big[ log\Big(\frac{1}{\sqrt{2\pi}}\Big)-log \sigma_k - \frac{1}{2\sigma^2_k}(x_j-\mu_k)^2\Big]\Big\}$$

3.$EM$算法$M$步：

迭代的$M$步是求函数$Q(\theta, \theta^{(i)})$对$\theta$的极大值，即求新一轮迭代的模型参数：

θ(i+1)=argmaxθQ(θ,θ(i)) θ ( i + 1 ) = a r g max θ Q ( θ , θ ( i ) )

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

求 $\hat{\mu_k}, \hat{\sigma^2_k}$只需要将Q函数对$\mu_k,\sigma^2_k$求偏导并令为0，即可得到；求$\hat{\alpha_k}$是在$\sum_{i=k}^{K} \alpha_k= 1$条件下求偏导数并令为0得到的，用到拉格朗日乘子法。
结果如下：

$$\hat{\mu_k} = \frac{\sum_{j=1}^N \hat{\gamma_{jk}}x_j}{\sum_{j=1}^N \hat{\gamma_{jk}}}, k = 1,2,...,K \\ \hat{\sigma_k^2 }= \frac{\sum_{j=1}^N \hat{\gamma_{jk}}(x_j - \mu_k)^2}{\sum_{j=1}^N \hat{\gamma_{jk}}}, k= 1,2,...,K\\ \hat{\alpha_k} =\frac{n_k}{N} = \frac{\sum_{j=1}^N \hat{\gamma_{jk}}}{N}, k = 1,2,...,K$$
重复以上计算，直到对数似然函数值不再有明显变化为止。

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

输入：观测数据$x_1,x_2,...,x_N$，高斯混合模型

输出：高斯混合模型参数

（1）取参数的初始值开始迭代
（2）E步：依据当前模型参数，计算分模型$k$对观测数据$x_i$的响应度

$$\hat{\gamma_{jk}} =\frac{\alpha_k \phi(x_j|\theta_k)}{\sum_{k=1}^K \alpha_k \phi(x_j|\theta_k)}, j = 1,2,...,N; k = 1, 2,...,K$$
（3）M步：计算新一轮迭代的模型参数
$$\hat{\mu_k} = \frac{\sum_{j=1}^N \hat{\gamma_{jk}}x_j}{\sum_{j=1}^N \hat{\gamma_{jk}}}, k = 1,2,...,K \\ \hat{\sigma_k^2 }= \frac{\sum_{j=1}^N \hat{\gamma_{jk}}(x_j - \mu_k)^2}{\sum_{j=1}^N \hat{\gamma_{jk}}}, k= 1,2,...,K\\ \hat{\alpha_k} =\frac{n_k}{N} = \frac{\sum_{j=1}^N \hat{\gamma_{jk}}}{N}, k = 1,2,...,K$$
（4）重复第（2）步和第（3）步，直到收敛。

至此，高斯混合模型的$EM$算法我们就介绍完了，关于算法 M M <script type="math/tex" id="MathJax-Element-14350">M</script>步中的计算，这里没有给出，求解方法会使用到带约束条件的拉格朗日乘子法
参考资料：
李航《统计学习方法》，徐亦达机器学习