高斯混合模型

阅读量：2173 次

发布时间：2019-05-01

本文共 1438 字，大约阅读时间需要 4 分钟。

定义

先给出高斯混合模型的定义，高斯混合模型是指具有如下形式的概率分布模型：

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

其中， $\alpha_k$ 是系数，且 $\alpha_k \geq 0$ ， $\sum_{k=1}^{K}\alpha_k=1$ ，而 $\phi(y|\theta_k)$ 是高斯分布密度， $\theta_k=(\mu_k, \sigma_k^2)$ ，对于随机变量y是一维数据时，

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

称为第k个分模型。理论上只要分模型足够多，并且各分模型的系数设置合理，就能够产生任意分布的样本。

分析

高斯混合模型属于生成模型，可以设想观测数据 $y_j$ ， $j=1,2,...,N$ ，是这样生成的：首先以概率 $\alpha_k$ 选择第k个分模型，然后由第k个分模型的概率分布生成观测数据 $y_j$ 。观测数据是能直接观测到的，已知的；而反映观测数据 $y_j$ 来自第k个分模型的数据是未知的，称为隐随机变量。一般地，用y表示观测随机变量的数据；z表示隐随机变量的数据。y称为不完全数据，而y和z连在一起称为完全数据。

为了求出模型的参数，我们先采用最大似然估计，似然函数可以写为:

$l(y|\theta)=\prod _{j=1}^{N}p(y_j|\theta)=\prod _{j=1}^N \sum_{k=1}^K\alpha_k\phi(y_j|\theta_k)=\prod _{j=1}^N\sum_{k=1}^Kp(z_j|\theta_k)p(y_j|z_j;\theta_k)$ （3)

对数化后得到：

$logl(y|\theta)=\sum_{j=1}^N log\sum_{k=1}^K p(z_j|\theta_k)p(y_j|z_j;\theta_k)$ (4)

然而(4)式无法直接通过求导并令导数为0的方法来求得参数的极大似然估计，但是如果知道了每个样本具体来自哪个分模型，即 $z_j$ ，那么（4）式便可以简化为：

$logl(y|\theta)=\sum_{j=1}^N log(p(z_j|\theta_k)p(y_j|z_j;\theta_k))=\sum_{j=1}^N(log\ p(z_j|\theta_k)+log\ p(y_j|z_j;\theta_k))$ (5)

此时便可以通过令偏导为零的方法求其最大似然估计。但实际上 $z_j$ 是无法获知的，既然无法获知那么能不能使用某种方法对其进行一个估计呢？因此这里需要采用EM算法，EM算法中使用关于z的期望来猜测z。

推导

为了方便推导，这里隐变量用 $\gamma _j$ 来表示， $\gamma _j$ 是一个k维的随机变量，当第j个观测来自第k个分模型时 $\gamma _{jk}=1$ ，否则 $\gamma_{jk}=0$ ， $\gamma_{jk}$ 是0-1随机变量。很明显 $\sum_{k=1}^K \gamma_{jk}=1$ 。

有了观测数据 $y_j$ 和未观测数据 $\gamma_{jk}$ ，那么完全数据就是: $(y_j,\gamma_{j1},\gamma_{j2},...,\gamma_{jk}),j=1,2,...,N$ 。于是，可以写出完全数据的似然函数：

$p(y,\gamma|\theta)=\prod _{j=1}^{N}p(y_j,\gamma_{j1},\gamma_{j2},...,\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 [\phi(y_j|\theta_k)]^{\gamma_{jk}}$ (6)

其中 $n_k=\sum _{j=1}^N \gamma_{jk}$ 表示N个数据中有多少个来自于第k个分模型，另外，很明显 $\sum_{k=1}^K n_k = N$ 。与(6)相比，(3)中没有出现隐变量，可以理解为不完全数据的似然函数，而(6)中限定了 $y_j$ 只会来源于某个分模型。

那么完全数据的对数似然函数为：

$log\ p(y,\gamma,|\theta)=\sum_{k=1}^K n_klog\alpha_k+\sum_{j=1}^N\gamma_{jk}\phi(y_j|\theta_j)$ (7)

(7)中存在隐变量 $\gamma$ ，使得我们无法直接使用极大似然估计来求解最优参数。既然 $\gamma$ 无法具体获知，那么我们也许可以对其进行一个估计，由此引入EM算法！

EM算法求解参数

接下来使用EM算法来求解参数，EM算法分为E步，求期望；和M步，求期望最大化时的参数值。

E-step：确定Q函数

既然无法直接对(7)进行似然估计，那么就转为对Q函数进行似然估计。Q函数表示完全数据的对数似然函数 $log\ p(y,z|\theta)$ 关于在给定观测数据y和当前参数 $\theta^{(i)}$ 下对未观测数据z的条件概率分布 $p(z|y,\theta^{(i)})$ 的期望。

$Q(\theta,\theta^{(i)})=E[log\ p(y,\gamma|\theta)|y,\theta^{(i)}]=E[\sum_{k=1}^K n_klog\alpha_k+\sum_{j=1}^N\gamma_{jk}\phi(y_j|\theta_j)] =\sum_{k=1}^K\left \{ \sum_{j=1}^N E(\gamma_{jk}|y,\theta)log\alpha_k+\sum_{j=1}^NE(\gamma_{jk}|y,\theta)\phi(y_j|\theta_j) \right \}$ (8)

这里将值不确定的随机变量 $\gamma_{jk}$ 换成了值确定的 $E(\gamma_{jk=1}|y_j,\theta)$ 。Q函数表达的其实就是先对样本 $y_j$ 做一个最有可能的划分，即以 $p(\gamma_{jk}|y_j,\theta^{(i)})$ 的概率确定样本 $\gamma_j$ 来自第k个分模型，再描述产生这个样本的可能性。有了对 $\gamma_{jk}$ 的估计后，似然估计中就不含隐变量，从而可以最大似然估计便可派上用场了。

将 $E(\gamma_{jk}|y,\theta)$ 记为 $\hat{\gamma_{jk}}$ ，这里计算 $\hat{\gamma_{jk}}$ 使用了贝叶斯公式：

$\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}^Kp(\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)}$ (9)

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

M-step

通过最大化Q函数便可以求新一轮的模型参数：

$\theta^{(i+1)}=\underset{\theta}{argmax}Q(\theta,\theta^{(i)})$ (10)

用 $\hat\mu_k$ ， $\hat \sigma^2$ 和 $\hat \alpha_k$ ,k=1,2,...,K，表示 $\theta^{(i+1)}$ 的各参数，用Q函数分别对 $\hat\mu_k$ ， $\hat \sigma^2$ 和 $\hat \alpha_k$ 求偏导，并令其为0得到第（i+1)轮迭代的参数值：