高斯混合模型GMM

高斯混合模型（Gaussian Mixed Model）是多个高斯分布函数的线性组合，理论上可以拟合出任意类型的分布，通常用于无监督聚类问题

设随机变量X，高斯混合模型如下所示
p(x)=\sum_{k=1}^K{\pi_k}N(x|\mu_k,\Sigma_k)
其中N(x|\mu_k,\Sigma_k)是混合模型中的第k个分量
\pi_k是每个分量的权重/概率

EM 估计GMM参数

由于\pi_k的引入，GMM无法通过极大似然计算模型参数，不过我们可以使用EM算法迭代求解

GMM中EM算法步骤

假定有C个数据 x_1,x_2,...,x_C，需要将其聚类到包含K个子模型的GMM中
0）首先初始化参数
\pi_k、\mu_k、\Sigma_k，k=1,2,...K

1）E-STEP
根据当前参数，计算每个数据在每个子模型的概率

\gamma_i^k=\frac{\pi_kN(x_i|\mu_k,\Sigma_k)}{\sum_{m=1}^K{\pi_mN(x_i|\mu_m,\Sigma_m)}},i=1,2,...,C;k=1,2,...,K \tag{1}

2）M-STEP
更新模型参数

\left\{\begin{aligned}\mu_k & = & \frac{\sum_i^C{\gamma_i^kx_i}}{\sum_i^C\gamma_i^k}\\ \Sigma_k & = & \frac{\sum_i^C\gamma_i^k(x_i-\mu_k)(x_i-\mu_k)^T}{\sum_i^C\gamma_i^k}\\ \pi_k & = & \frac{\sum_{i=1}^C\gamma_i^k}{C}\\ \end{aligned}\right. \tag{2}

重复E-STEP和M-STEP，直到参数收敛

EM算法推导

内容参考自小象学院的机器学习课程，讲的很好，推荐

EM算法一般推导

对于目标函数p(x:\theta)，取对数似然函数（连乘取对数）：

l(\theta)=\sum_{i=1}^M{log\ p(x;\theta)}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ =\sum_{i=1}^M{log{\sum_zp(x,z;\theta)}} \tag{3}

其中z是隐变量

下面我们看看GMM的对数似然函数：

l_{\pi,\mu,\Sigma}(x)=\sum_{i=1}^M{log\sum_{k=1}^K{\pi_kN(x_i|\mu_k,\Sigma_k)}} \tag{4}

上面的隐变量z就是GMM里的\pi_k的变量

由于隐随机变量的存在，无法直接估计参数，因此EM算法的基本策略是：计算对数似然函数的下界，求该下界的最大值，重复这个过程，直到收敛于局部最大值

图1. 策略示意（图片取自小象学院机器学习课程PPT）

图1紫色线假设是对数似然函数的曲线，现在因为隐变量z的存在，不方便直接求哪组\theta使得其取极值

现在给定一个初始值\theta_0（绿线与横轴交点）构造一个函数r(x|\theta)，该点取值A，并保证除了\theta_0外该函数的值均小于p(x|\theta)，即

r(x|\theta)\leq p(x|\theta) \tag{5}

接下来在r(x|\theta)上求极大值B（示意图红色线与蓝线交点），因为上式的约束条件，该点在p(x|\theta)的取值C（红线与紫线的交点）必然大于B，又因为B大于A，因此C必然大于A

重复上述过程就可以逐步逼近p(x|\theta)的局部最大值（从上述过程也可以知道，EM只能保证局部最优无法保证全局最优 ）

那么怎么构造这个辅助函数 r(x|\theta)呢
假设Q(z)是关于z的一个分布，式(1)可变为

l(\theta)=\sum_{i=1}^M{log\sum_{z^i}p(x^i,z^i;\theta)}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ =\sum_{i=1}^M{log\sum_z^iQ_i(z^i)\frac{p(x^i,z^i;\theta)}{Q_i(z^i)}}\\ \geq \sum_{i=1}^M{\sum_{z^i}Q_i(z^i)log\frac{p(x^i,z^i;\theta)}{Q_i(z^i)}}\ \ \tag{6}

式(4)中令\frac{p(x^i,z^i;\theta)}{Q_i(z^i)}=X，那么
log\sum_z^i{Q_i(z^i)X}就是对X关于Q_i(z^i)求期望再取对数 ：log{E_Q(X)}
\sum_{z^i}Q_i(z^i){logX}就是对X取对数后求期望 ：E_Q(logX)

又因为对数函数log是一个凹函数 ，根据Jensen不等式 ，对于凹函数f，总有

f(Ex)\geq Ef(x) \tag{7}

因此(6)式成立，其中的\sum_{i=1}^M{\sum_{z^i}Q_i(z^i)log\frac{p(x^i,z^i;\theta)}{Q_i(z^i)}}就是要找的辅助函数r(x|\theta)

再来看图1,需要保证找到的r(x|\theta)在某个位置与l(\theta)相切（6式取等号）（否则不能保证迭代的点往局部最大点接近）；
而log是个严格的凹函数，因此只有当\frac{p(x,z;\theta)}{Q_i(z)}=C取定值才可以，也就是说p(x,z;\theta)\propto Q_i(z)是正比关系，可以推得

Q_i(z^i)=\frac{p(x^i,z^i;\theta)}{\sum_z{p(x^i,z^i;\theta)}}\\ =\frac{p(x^i|z^i;\theta)}{p(x^i;\theta)}\\ =p(z^i|x;\theta) \ \ \ \tag{8}

也就是说，当Q(z)的分布取 给定样本时隐变量 z的条件概率 时等号可以成立

由此可以得到EM算法的步骤 ：
1）首先对所有数据，根据现有的参数求条件概率Q(z)，此为E-STEP；
2）然后将得到的Q(z)代入(4)式所述的辅助函数r(x|\theta)，求使得该辅助函数取最大值的\theta，更新参数，此为M-STEP；
重复1-2步，最终算法收敛于局部极值

EM算法在GMM中的应用推导

E-STEP：
对所有数据样本，应用现有的参数（\pi、\mu、\Sigma）,计算隐变量的条件概率

w_j^i=Q_i(z^i=j)=P(z^i=j|x^i;\phi,\mu,\Sigma) \tag{9}

对比式3和式4，w_j^i表示第i个样本属于第j个子模型的概率

M-STEP：
该步骤将条件概率w_j^i代入辅助函数，并求该函数的最大值

\sum_{i=1}^M{\sum_{z^i}Q(z^i)log{\frac{p(x^i,z^i;\phi,\mu,\Sigma)}{Q_i(z^i)}}}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ =\sum_{i=1}^M{\sum_{j=1}^K{Q_i(z^i=j)log{\frac{p(x^i|z^i=j;\mu,\Sigma)p(z^i=j;\phi)}{Q_i(z^i=j)}}}}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ =\sum_{i=1}^M{\sum_{j=1}^K{w_j^ilog\frac{\frac{1}{(2\pi)^{n/2|\Sigma_j|^{1/2}}}exp(-\frac{1}{2}(x^i-\mu_j)^T\Sigma_j^{-1}(x^i-\mu_j))\phi_j}{w_j^i}}} \tag{10}

上式(10)中，p(z^i=j;\phi)=\phi_j
而指定z^i=j后，p(x^i|z^i=j;\mu,\Sigma)就是一个关于\mu和\Sigma的高斯分布

M-STEP的目标就是找到使得式(10)取最大值的\phi、\mu、\Sigma
1) 对 \mu_l求偏导

\partial_{\mu_l}\sum_{i=1}^M{\sum_{j=1}^K{w_j^ilog\frac{\frac{1}{(2\pi)^{n/2|\Sigma_j|^{1/2}}}exp(-\frac{1}{2}(x^i-\mu_j)^T\Sigma_j^{-1}(x^i-\mu_j))\phi_j}{w_j^i}}}\\ =-\partial_{\mu_l}\sum_{i=1}^M{\sum_{j=1}^K{w_j^i\frac{1}{2}(x^i-\mu_j)^T\Sigma_j^{-1}(x^i-\mu_j)}}\\ =\sum_{i=1}^M{w_l^i(\Sigma_l^{-1}x^i-\Sigma_l^{-1}\mu_l)}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ =\Sigma_l^{-1}{\sum_{i=1}^M(w_l^ix^i-\mu_l)}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \tag{11}

上式(9)中，因为\Sigma_j^{-1}是对称阵，因此有\frac{\partial(x^TAx)}{\partial x}=2Ax
令(11)式为0，可以推得

\mu_l=\frac{\sum_{i=1}^Mw_l^ix^i}{\sum_{i=1}^Mw_l^i} \tag{12}

2）对 \Sigma_l求偏导
同上面的方法，可以求得

\Sigma_l=\frac{\sum_{i=1}^M{w_j^i(x^i-\mu_j)(x^i-\mu_j)^T}}{\sum_{i=1}^Mw_j^i} \tag{13}

3）对 \phi求偏导

\partial_{\phi}\sum_{i=1}^M{\sum_{j=1}^K{w_j^ilog\frac{\frac{1}{(2\pi)^{n/2|\Sigma_j|^{1/2}}}exp(-\frac{1}{2}(x^i-\mu_j)^T\Sigma_j^{-1}(x^i-\mu_j))\phi_j}{w_j^i}}}\\ =\partial_{\phi}\sum_{i=1}^M{\sum_{j=1}^Kw_j^ilog\phi_j}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \tag{14}

注意到\phi有约束条件\sum_{j=1}^K\phi_j=1，使用拉格朗日乘子法，建立目标函数：

L(\phi)=\sum_{i=1}^M{\sum_{j=1}^Kw_j^ilog\phi_j}+\beta(\sum_{j=1}^K\phi_j-1) \tag{15}

式(15)对\phi_l求偏导：

\frac{\partial}{\partial \phi_j}L(\phi)=\sum_{i=1}^M{\frac{w_j^i}{\phi_j}}+\beta \tag{16}

式子(16)对每一个j都要为0（取极值），不妨将所有j累加，有：

\sum_{j=1}^K{\sum_{i=1}^Mw_j^i}+\beta\sum_{j=1}^K\phi_j=0 \tag{17}

注意到\sum_{j=1}^K\phi_j=1，而\sum_{j=1}^K{\sum_{i=1}^Mw_j^i}=\sum_{i=1}^M{\sum_{j=1}^Kw_j^i}，于是式(17)等价为

\sum_{i=1}^M{\sum_{j=1}^Kw_j^i}=-\beta \tag{18}

显然，式(18)的\sum_{j=1}^Kw_j^i=1对每个样本i都成立（某个样本属于所有类的概率之和为1），因此可以推得

\beta=-M \tag{19}

最后将式(19)代入式(17)，可以求得\phi

\phi_l=\frac{1}{M}\sum_{i=1}^Mw_j^i \tag{20}

至此，GMM中EM算法的迭代公式全部推导完成，式子(9)、(12)、(13)和(20)分别对应式(1)和(2)。

应用于GMM中的EM算法代码

python中的sklearn包提供了EM算法接口

    from sklearn import mixture
    gmm = mixture.GaussianMixture(n_components=nclass, covariance_type='full').fit(sample_)
    labels = gmm.predict(sample_)

当然，E-STEP和M-STEP的内容其实很简单很容易编程，自己写也没多少代码量
需要注意的是有可能由于初始化不合适的原因，E-STEP的\sum_{i=1}^Mw_j^i有可能为0，而该项在M-STEP中被作为分母计算，需要考虑保护。

–以上–