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

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

EM 估计GMM参数

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

GMM中EM算法步骤

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

1）E-STEP
根据当前参数，计算每个数据在每个子模型的概率
γ i k = π k N ( x i ∣ μ k , Σ k ) ∑ m = 1 K π m N ( x i ∣ μ m , Σ m ) , i = 1 , 2 , . . . , C ; k = 1 , 2 , . . . , K (1) \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} γik=∑m=1KπmN(xi∣μm,Σm)πkN(xi∣μk,Σk),i=1,2,...,C;k=1,2,...,K(1)
2）M-STEP
更新模型参数
{ μ k = ∑ i C γ i k x i ∑ i C γ i k Σ k = ∑ i C γ i k ( x i − μ k ) ( x i − μ k ) T ∑ i C γ i k π k = ∑ i = 1 C γ i k C (2) \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} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧μkΣkπk===∑iCγik∑iCγikxi∑iCγik∑iCγik(xi−μk)(xi−μk)TC∑i=1Cγik(2)

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

EM算法推导

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

EM算法一般推导

对于目标函数 p ( x : θ ) p(x:\theta) p(x:θ)，取对数似然函数（连乘取对数）：
l ( θ ) = ∑ i = 1 M l o g p ( x ; θ ) = ∑ i = 1 M l o g ∑ z p ( x , z ; θ ) (3) l(\theta)=\sum_{i=1}^M{log\ p(x;\theta)}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\=\sum_{i=1}^M{log{\sum_zp(x,z;\theta)}} \tag{3} l(θ)=i=1∑Mlog p(x;θ) =i=1∑Mlogz∑p(x,z;θ)(3)
其中 z z z是隐变量

下面我们看看GMM的对数似然函数：
l π , μ , Σ ( x ) = ∑ i = 1 M l o g ∑ k = 1 K π k N ( x i ∣ μ k , Σ k ) (4) l_{\pi,\mu,\Sigma}(x)=\sum_{i=1}^M{log\sum_{k=1}^K{\pi_kN(x_i|\mu_k,\Sigma_k)}} \tag{4} lπ,μ,Σ(x)=i=1∑Mlogk=1∑KπkN(xi∣μk,Σk)(4)
上面的隐变量 z z z就是GMM里的 π k \pi_k πk的变量

由于隐随机变量的存在，无法直接估计参数，因此EM算法的基本策略是：计算对数似然函数的下界，求该下界的最大值，重复这个过程，直到收敛于局部最大值
高斯混合模型GMM-LMLPHP
图1. 策略示意（图片取自小象学院机器学习课程PPT）

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

现在给定一个初始值 θ 0 \theta_0 θ0（绿线与横轴交点）构造一个函数 r ( x ∣ θ ) r(x|\theta) r(x∣θ)，该点取值A，并保证除了 θ 0 \theta_0 θ0外该函数的值均小于 p ( x ∣ θ ) p(x|\theta) p(x∣θ)，即
r ( x ∣ θ ) ≤ p ( x ∣ θ ) (5) r(x|\theta)\leq p(x|\theta) \tag{5} r(x∣θ)≤p(x∣θ)(5)
接下来在 r ( x ∣ θ ) r(x|\theta) r(x∣θ)上求极大值B（示意图红色线与蓝线交点），因为上式的约束条件，该点在 p ( x ∣ θ ) p(x|\theta) p(x∣θ)的取值C（红线与紫线的交点）必然大于B，又因为B大于A，因此C必然大于A

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

那么怎么构造这个辅助函数 r ( x ∣ θ ) r(x|\theta) r(x∣θ)呢
假设 Q ( z ) Q(z) Q(z)是关于 z z z的一个分布，式(1)可变为
l ( θ ) = ∑ i = 1 M l o g ∑ z i p ( x i , z i ; θ ) = ∑ i = 1 M l o g ∑ z i Q i ( z i ) p ( x i , z i ; θ ) Q i ( z i ) ≥ ∑ i = 1 M ∑ z i Q i ( z i ) l o g p ( x i , z i ; θ ) Q i ( z i ) (6) 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} l(θ)=i=1∑Mlogzi∑p(xi,zi;θ) =i=1∑Mlogz∑iQi(zi)Qi(zi)p(xi,zi;θ)≥i=1∑Mzi∑Qi(zi)logQi(zi)p(xi,zi;θ) (6)

式(4)中令 p ( x i , z i ; θ ) Q i ( z i ) = X \frac{p(x^i,z^i;\theta)}{Q_i(z^i)}=X Qi(zi)p(xi,zi;θ)=X，那么
l o g ∑ z i Q i ( z i ) X log\sum_z^i{Q_i(z^i)X} log∑ziQi(zi)X就是对 X X X关于 Q i ( z i ) Q_i(z^i) Qi(zi)求期望再取对数： l o g E Q ( X ) log{E_Q(X)} logEQ(X)
∑ z i Q i ( z i ) l o g X \sum_{z^i}Q_i(z^i){logX} ∑ziQi(zi)logX就是对 X X X取对数后求期望： E Q ( l o g X ) E_Q(logX) EQ(logX)

又因为对数函数 l o g log log是一个凹函数，根据Jensen不等式，对于凹函数 f f f，总有
f ( E x ) ≥ E f ( x ) (7) f(Ex)\geq Ef(x) \tag{7} f(Ex)≥Ef(x)(7)

因此(6)式成立，其中的 ∑ i = 1 M ∑ z i Q i ( z i ) l o g p ( x i , z i ; θ ) Q i ( z i ) \sum_{i=1}^M{\sum_{z^i}Q_i(z^i)log\frac{p(x^i,z^i;\theta)}{Q_i(z^i)}} ∑i=1M∑ziQi(zi)logQi(zi)p(xi,zi;θ)就是要找的辅助函数 r ( x ∣ θ ) r(x|\theta) r(x∣θ)

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

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

由此可以得到EM算法的步骤：
1）首先对所有数据，根据现有的参数求条件概率 Q ( z ) Q(z) Q(z)，此为E-STEP；
2）然后将得到的 Q ( z ) Q(z) Q(z)代入(4)式所述的辅助函数 r ( x ∣ θ ) r(x|\theta) r(x∣θ)，求使得该辅助函数取最大值的 θ \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 ; ϕ , μ , Σ ) (9) w_j^i=Q_i(z^i=j)=P(z^i=j|x^i;\phi,\mu,\Sigma) \tag{9} wji=Qi(zi=j)=P(zi=j∣xi;ϕ,μ,Σ)(9)
对比式3和式4， w j i w_j^i wji表示第 i i i个样本属于第 j j j个子模型的概率

M-STEP：
该步骤将条件概率 w j i w_j^i wji代入辅助函数，并求该函数的最大值
∑ i = 1 M ∑ z i Q ( z i ) l o g p ( x i , z i ; ϕ , μ , Σ ) Q i ( z i ) = ∑ i = 1 M ∑ j = 1 K Q i ( z i = j ) l o g p ( x i ∣ z i = j ; μ , Σ ) p ( z i = j ; ϕ ) Q i ( z i = j ) = ∑ i = 1 M ∑ j = 1 K w j i l o g 1 ( 2 π ) n / 2 ∣ Σ j ∣ 1 / 2 e x p ( − 1 2 ( x i − μ j ) T Σ j − 1 ( x i − μ j ) ) ϕ j w j i (10) \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} i=1∑Mzi∑Q(zi)logQi(zi)p(xi,zi;ϕ,μ,Σ) =i=1∑Mj=1∑KQi(zi=j)logQi(zi=j)p(xi∣zi=j;μ,Σ)p(zi=j;ϕ) =i=1∑Mj=1∑Kwjilogwji(2π)n/2∣Σj∣1/21exp(−21(xi−μj)TΣj−1(xi−μj))ϕj(10)
上式(10)中， p ( z i = j ; ϕ ) = ϕ j p(z^i=j;\phi)=\phi_j p(zi=j;ϕ)=ϕj
而指定 z i = j z^i=j zi=j后， p ( x i ∣ z i = j ; μ , Σ ) p(x^i|z^i=j;\mu,\Sigma) p(xi∣zi=j;μ,Σ)就是一个关于 μ \mu μ和 Σ \Sigma Σ的高斯分布

M-STEP的目标就是找到使得式(10)取最大值的 ϕ \phi ϕ、 μ \mu μ、 Σ \Sigma Σ
1) 对 μ l \mu_l μl求偏导
∂ μ l ∑ i = 1 M ∑ j = 1 K w j i l o g 1 ( 2 π ) n / 2 ∣ Σ j ∣ 1 / 2 e x p ( − 1 2 ( x i − μ j ) T Σ j − 1 ( x i − μ j ) ) ϕ j w j i = − ∂ μ l ∑ i = 1 M ∑ j = 1 K w j i 1 2 ( x i − μ j ) T Σ j − 1 ( x i − μ j ) = ∑ i = 1 M w l i ( Σ l − 1 x i − Σ l − 1 μ l ) = Σ l − 1 ∑ i = 1 M ( w l i x i − μ l ) (11) \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} ∂μli=1∑Mj=1∑Kwjilogwji(2π)n/2∣Σj∣1/21exp(−21(xi−μj)TΣj−1(xi−μj))ϕj=−∂μli=1∑Mj=1∑Kwji21(xi−μj)TΣj−1(xi−μj)=i=1∑Mwli(Σl−1xi−Σl−1μl) =Σl−1i=1∑M(wlixi−μl) (11)
上式(9)中，因为 Σ j − 1 \Sigma_j^{-1} Σj−1是对称阵，因此有 ∂ ( x T A x ) ∂ x = 2 A x \frac{\partial(x^TAx)}{\partial x}=2Ax ∂x∂(xTAx)=2Ax
令(11)式为0，可以推得
μ l = ∑ i = 1 M w l i x i ∑ i = 1 M w l i (12) \mu_l=\frac{\sum_{i=1}^Mw_l^ix^i}{\sum_{i=1}^Mw_l^i} \tag{12} μl=∑i=1Mwli∑i=1Mwlixi(12)

2）对 Σ l \Sigma_l Σl求偏导
同上面的方法，可以求得
Σ l = ∑ i = 1 M w j i ( x i − μ j ) ( x i − μ j ) T ∑ i = 1 M w j i (13) \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} Σl=∑i=1Mwji∑i=1Mwji(xi−μj)(xi−μj)T(13)

3）对 ϕ \phi ϕ求偏导
∂ ϕ ∑ i = 1 M ∑ j = 1 K w j i l o g 1 ( 2 π ) n / 2 ∣ Σ j ∣ 1 / 2 e x p ( − 1 2 ( x i − μ j ) T Σ j − 1 ( x i − μ j ) ) ϕ j w j i = ∂ ϕ ∑ i = 1 M ∑ j = 1 K w j i l o g ϕ j (14) \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} ∂ϕi=1∑Mj=1∑Kwjilogwji(2π)n/2∣Σj∣1/21exp(−21(xi−μj)TΣj−1(xi−μj))ϕj=∂ϕi=1∑Mj=1∑Kwjilogϕj (14)

注意到 ϕ \phi ϕ有约束条件 ∑ j = 1 K ϕ j = 1 \sum_{j=1}^K\phi_j=1 ∑j=1Kϕj=1，使用拉格朗日乘子法，建立目标函数：
L ( ϕ ) = ∑ i = 1 M ∑ j = 1 K w j i l o g ϕ j + β ( ∑ j = 1 K ϕ j − 1 ) (15) L(\phi)=\sum_{i=1}^M{\sum_{j=1}^Kw_j^ilog\phi_j}+\beta(\sum_{j=1}^K\phi_j-1) \tag{15} L(ϕ)=i=1∑Mj=1∑Kwjilogϕj+β(j=1∑Kϕj−1)(15)
式(15)对 ϕ l \phi_l ϕl求偏导：
∂ ∂ ϕ j L ( ϕ ) = ∑ i = 1 M w j i ϕ j + β (16) \frac{\partial}{\partial \phi_j}L(\phi)=\sum_{i=1}^M{\frac{w_j^i}{\phi_j}}+\beta \tag{16} ∂ϕj∂L(ϕ)=i=1∑Mϕjwji+β(16)
式子(16)对每一个 j j j都要为0（取极值），不妨将所有 j j j累加，有：
∑ j = 1 K ∑ i = 1 M w j i + β ∑ j = 1 K ϕ j = 0 (17) \sum_{j=1}^K{\sum_{i=1}^Mw_j^i}+\beta\sum_{j=1}^K\phi_j=0 \tag{17} j=1∑Ki=1∑Mwji+βj=1∑Kϕj=0(17)
注意到 ∑ j = 1 K ϕ j = 1 \sum_{j=1}^K\phi_j=1 ∑j=1Kϕj=1，而 ∑ j = 1 K ∑ i = 1 M w j i = ∑ i = 1 M ∑ j = 1 K w j i \sum_{j=1}^K{\sum_{i=1}^Mw_j^i}=\sum_{i=1}^M{\sum_{j=1}^Kw_j^i} ∑j=1K∑i=1Mwji=∑i=1M∑j=1Kwji，于是式(17)等价为
∑ i = 1 M ∑ j = 1 K w j i = − β (18) \sum_{i=1}^M{\sum_{j=1}^Kw_j^i}=-\beta \tag{18} i=1∑Mj=1∑Kwji=−β(18)
显然，式(18)的 ∑ j = 1 K w j i = 1 \sum_{j=1}^Kw_j^i=1 ∑j=1Kwji=1对每个样本 i i i都成立（某个样本属于所有类的概率之和为1），因此可以推得
β = − M (19) \beta=-M \tag{19} β=−M(19)
最后将式(19)代入式(17)，可以求得 ϕ \phi ϕ
ϕ l = 1 M ∑ i = 1 M w j i (20) \phi_l=\frac{1}{M}\sum_{i=1}^Mw_j^i \tag{20} ϕl=M1i=1∑Mwji(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的 ∑ i = 1 M w j i \sum_{i=1}^Mw_j^i ∑i=1Mwji有可能为0，而该项在M-STEP中被作为分母计算，需要考虑保护。

–以上–