1 简介
EM（Expectation-Maximum）算法也称期望最大化算法，它是为了解决在方程无法获得解析解的情况下，通过迭代给出数值解。
核心：EM算法是一种迭代算法，用于含有隐变量的概率模型参数的极大似然估计（因此在往下面看之前，我希望你对贝叶斯的基本理论有所了解）
2 极大似然估计
（1）问题背景
我们要去调查全校学生的身高分布，因此分别随机抽取男生、女生各100人，总共200人（假设全校的男生、女生一样多，如果不一样多，那就还需要根据男女比例进行抽取，但这并不是我们研究的重点）。我们想要去获取男生的身高分布
（2）问题假设
① 男生、女生的身高符合高斯正态分布（对于非高斯问题，则需要建模处理，它不是统计学科要解决的问题）
② 男生、女生的身高都是独立同分布（ i.i.d ）
（3）模型建立
现在我们已经有了一堆数据 X = {
x
1
,
x
2
,
.
.
.
,
x
n
{x_1,x_2,...,x_n}
x1,x2,...,xn} ，其中
x
i
x_i
xi表示第i个男生的身高，n为男生的个数。待估参数为
θ
=
\theta=
θ={
μ
,
σ
{\mu ,\sigma }
μ,σ} 。
我们从当前全校男生中抽取一个同学A，其身高为
y
A
y_A
yA的概率可以记为
P
(
x
A
∣
θ
)
P(x_A \mid \theta)
P(xA∣θ)。那么从全校男生中抽取100个同学，其身高为上述数据集合Y的概率可以用下列式子表示：
L
(
θ
)
=
L
(
x
1
,
⋯
,
x
n
;
θ
)
=
∏
i
=
1
n
p
(
x
i
;
θ
)
,
θ
∈
Θ
L(\theta)=L\left(x_1, \cdots, x_n ; \theta\right)=\prod_{i=1}^n p\left(x_i ; \theta\right), \theta \in \Theta
L(θ)=L(x1,⋯,xn;θ)=i=1∏np(xi;θ),θ∈Θ
上述式子反应了，在概率密度函数的参数是θ时，得到X这组样本的概率。上述式子中只有 θ 是未知的，因此称为参数θ的似然函数
L
(
θ
)
L(\theta)
L(θ) 。我们的目标就是寻找出最优的
θ
\theta
θ，使得
L
(
θ
)
L(\theta)
L(θ) 的数值最大，即：找到一个高斯分布，使得我随机抽100人，刚好就是上面的100个男生的概率最大。θ的最大似然估计量记为
θ
^
=
a
r
g
m
a
x
L
(
θ
)
\hat{\theta}=argmax L(\theta)
θ^=argmaxL(θ)
为便于分析，我们将
L
(
θ
)
L(\theta)
L(θ)取其对数，即：
ln
L
(
θ
)
=
log
∏
i
=
1
n
p
(
x
i
;
θ
)
=
∑
i
=
1
n
log
p
(
x
i
;
θ
)
\ln L(\theta)=\log \prod_{i=1}^n p\left(x_i ; \theta\right)=\sum_{i=1}^n \log p\left(x_i ; \theta\right)
lnL(θ)=logi=1∏np(xi;θ)=i=1∑nlogp(xi;θ)
由于p是高斯分布的概率密度函数（PDF），取对数后exp的指数会变为相加项。对其求导，令导数为0，得到似然方程；解似然方程即可得到
θ
M
L
E
\theta_{MLE}
θMLE，即
θ
\theta
θ最优估计。
3 EM算法推导
对于参数为θ且含有隐变量Z的概率模型，进行n次抽样。假设随机变量
x
x
x 的观察值为
X
X
X={
x
1
,
x
2
,
.
.
.
,
x
n
x_1,x_2,...,x_n
x1,x2,...,xn}，隐变量
Z
Z
Z的m个可能的取值为
Z
Z
Z={
z
1
,
z
2
,
.
.
.
,
z
m
z_1,z_2,...,z_m
z1,z2,...,zm}. （如果对
z
z
z存在疑惑，那么就先将其理解为不同模型的概率，其总和为1。后面在第四部分GMM，会有更加容易被理解的说明）
写出其似然函数：
EM算法的公式：
θ
(
t
+
1
)
=
arg
max
θ
∫
z
log
P
(
X
,
z
∣
θ
)
⋅
P
(
z
∣
X
,
θ
(
t
)
)
d
z
\theta^{(t+1)} = \arg\max_{\theta} \int_z \log P(X, z | \theta) \cdot P(z|X, \theta^{(t)}) dz
θ(t+1)=argθmax∫zlogP(X,z∣θ)⋅P(z∣X,θ(t))dz
- l o g P ( X , z ∣ θ ) log P(X, z | \theta) logP(X,z∣θ)是数据X的对数联合概率
- arg max θ ∫ z log P ( X , z ∣ θ ) ⋅ P ( z ∣ X , θ ( t ) ) d z \arg\max_{\theta} \int_z \log P(X, z | \theta) \cdot P(z|X, \theta^{(t)}) dz argmaxθ∫zlogP(X,z∣θ)⋅P(z∣X,θ(t))dz可以表示为 E z ∣ x , θ ( t ) [ log P ( x , z ∣ θ ) ] \mathbb{E}_{z|x,\theta^{(t)}}[\log P(x, z|\theta)] Ez∣x,θ(t)[logP(x,z∣θ)]
EM算法收敛性的证明
目标：证明下一次迭代的参数
θ
\theta
θ，要比当前时刻的
θ
\theta
θ好，即证明其收敛性，公式如下：
θ
(
t
+
1
)
←
θ
(
t
)
\theta^{(t+1)} \gets \theta^{(t)}
θ(t+1)←θ(t)
log
P
(
X
∣
θ
(
t
+
1
)
)
≥
log
P
(
X
∣
θ
(
t
)
)
（
1
）
\log P(X | \theta^{(t+1)})\ge\log P(X | \theta^{(t)})\text{ }（1）
logP(X∣θ(t+1))≥logP(X∣θ(t)) （1） 利用贝叶斯全概率公式可得：
P
(
X
∣
θ
)
=
P
(
X
,
θ
)
P
(
θ
)
=
P
(
X
,
θ
)
P
(
X
,
Z
,
θ
)
⋅
P
(
X
,
Z
,
θ
)
P
(
θ
)
=
P
(
X
∣
Z
,
θ
)
P
(
Z
∣
X
,
θ
)
\begin{align*} P(X|\theta) &= \frac{P(X, \theta)}{P(\theta)} \\ &= \frac{P(X, \theta)}{P(X, Z, \theta)} \cdot \frac{P(X, Z, \theta)}{P(\theta)} \\ &= \frac{P(X | Z, \theta)}{P(Z | X, \theta)} \end{align*}
P(X∣θ)=P(θ)P(X,θ)=P(X,Z,θ)P(X,θ)⋅P(θ)P(X,Z,θ)=P(Z∣X,θ)P(X∣Z,θ)因此：
log
P
(
X
∣
θ
)
=
log
P
(
X
,
z
∣
θ
)
P
(
z
∣
X
,
θ
)
=
log
P
(
X
,
z
∣
θ
)
−
P
(
z
∣
X
,
θ
)
（
2
）
\log P(X | \theta)=\log\frac{ P(X,z | \theta)}{P(z | X,\theta)}=\log{ P(X,z | \theta)}-{P(z | X,\theta)} \text{ }\text{ }（2）
logP(X∣θ)=logP(z∣X,θ)P(X,z∣θ)=logP(X,z∣θ)−P(z∣X,θ) （2） 对公式（2）两边同时求期望
E
z
∣
x
,
θ
(
t
)
[
log
P
(
X
∣
θ
)
]
=
E
z
∣
x
,
θ
(
t
)
[
log
P
(
X
,
z
∣
θ
)
−
log
P
(
z
∣
X
,
θ
)
]
\mathbb{E}_{z|x,\theta^{(t)}}[\log P(X|\theta)]=\mathbb{E}_{z|x,\theta^{(t)}}[\log P(X, z|\theta)-\log P(z |X,\theta)]
Ez∣x,θ(t)[logP(X∣θ)]=Ez∣x,θ(t)[logP(X,z∣θ)−logP(z∣X,θ)]
左边
=
∫
z
P
(
z
∣
X
,
θ
(
t
)
)
⋅
log
P
(
X
∣
θ
)
d
z
=
log
P
(
X
∣
θ
)
∫
z
P
(
z
∣
X
,
θ
(
t
)
)
d
z
=
log
P
(
X
∣
θ
)
\begin{align*}左边&= \int_z P(z|X, \theta^{(t)}) \cdot \log P(X|\theta) dz \\&= \log P(X|\theta) \int_z P(z|X, \theta^{(t)}) dz \\&= \log P(X|\theta) \\\end{align*}
左边=∫zP(z∣X,θ(t))⋅logP(X∣θ)dz=logP(X∣θ)∫zP(z∣X,θ(t))dz=logP(X∣θ)
注意：
∫
z
P
(
z
∣
X
,
θ
(
t
)
)
d
z
=
1
，可以将其理解为不同模型的占比，总的和为
1.
注意：\int_z P(z|X, \theta^{(t)}) dz=1，可以将其理解为不同模型的占比，总的和为1.
注意：∫zP(z∣X,θ(t))dz=1，可以将其理解为不同模型的占比，总的和为1.
右边
=
∫
z
P
(
z
∣
X
,
θ
(
t
)
)
⋅
log
P
(
X
,
z
∣
θ
)
d
z
−
∫
z
P
(
z
∣
X
,
θ
(
t
)
)
⋅
log
P
(
z
∣
X
,
θ
)
d
z
\begin{align*} 右边&= \int_z P(z|X, \theta^{(t)}) \cdot \log P(X, z|\theta) dz - \int_z P(z|X, \theta^{(t)}) \cdot \log P(z|X, \theta) dz & \end{align*}
右边=∫zP(z∣X,θ(t))⋅logP(X,z∣θ)dz−∫zP(z∣X,θ(t))⋅logP(z∣X,θ)dz 令
Q
(
θ
,
θ
(
t
)
)
=
∫
z
P
(
z
∣
X
,
θ
(
t
)
)
⋅
log
P
(
X
,
z
∣
θ
)
d
z
\text{Q}(\theta, \theta^{(t)})= \int_z P(z|X, \theta^{(t)}) \cdot \log P(X, z|\theta) dz
Q(θ,θ(t))=∫zP(z∣X,θ(t))⋅logP(X,z∣θ)dz ，
H
(
θ
,
θ
(
t
)
)
=
∫
z
P
(
z
∣
X
,
θ
(
t
)
)
⋅
log
P
(
z
∣
X
,
θ
)
d
z
\text{H}(\theta, \theta^{(t)})=\int_z P(z|X, \theta^{(t)}) \cdot \log P(z|X, \theta) dz
H(θ,θ(t))=∫zP(z∣X,θ(t))⋅logP(z∣X,θ)dz，则：
log
P
(
X
∣
θ
)
=
Q
(
θ
,
θ
(
t
)
)
−
H
(
θ
,
θ
(
t
)
)
(
3
)
\log P(X|\theta)=\text{Q}(\theta, \theta^{(t)})-\text{H}(\theta, \theta^{(t)}) \text{ }(3)
logP(X∣θ)=Q(θ,θ(t))−H(θ,θ(t)) (3)证明公式（3）成立，即证明：
Q
(
θ
(
t
)
,
θ
(
t
)
)
−
H
(
θ
(
t
)
,
θ
(
t
)
)
≤
Q
(
θ
(
t
+
1
)
,
θ
(
t
)
)
−
H
(
θ
(
t
+
1
)
,
θ
(
t
)
)
(
4
)
Q(\theta^{(t)}, \theta^{(t)}) - H(\theta^{(t)}, \theta^{(t)}) \leq Q(\theta^{(t+1)}, \theta^{(t)}) - H(\theta^{(t+1)}, \theta^{(t)})\text{ }(4)
Q(θ(t),θ(t))−H(θ(t),θ(t))≤Q(θ(t+1),θ(t))−H(θ(t+1),θ(t)) (4)
- 先证明 Q Q Q: Q ( θ ( t ) , θ ( t ) ) = ∫ z P ( z ∣ X , θ ( t ) ) ⋅ log P ( X , z ∣ θ ( t ) ) d z Q ( θ ( t + 1 ) , θ ( t ) ) = ∫ z P ( z ∣ X , θ ( t ) ) ⋅ log P ( X , z ∣ θ ( t + 1 ) ) d z Q(\theta^{(t)}, \theta^{(t)}) = \int_z P(z|X, \theta^{(t)}) \cdot \log P(X, z|\theta^{(t)}) dz \\ Q(\theta^{(t+1)}, \theta^{(t)}) = \int_z P(z|X, \theta^{(t)}) \cdot \log P(X, z|\theta^{(t+1)}) dz Q(θ(t),θ(t))=∫zP(z∣X,θ(t))⋅logP(X,z∣θ(t))dzQ(θ(t+1),θ(t))=∫zP(z∣X,θ(t))⋅logP(X,z∣θ(t+1))dz根据定义 θ ( t + 1 ) = arg max θ ∫ z P ( z ∣ X , θ ( t ) ) ⋅ log P ( X , z ∣ θ ) d z \theta^{(t+1)} = \arg\max_{\theta} \int_z P(z|X, \theta^{(t)}) \cdot \log P(X, z|\theta) dz θ(t+1)=argθmax∫zP(z∣X,θ(t))⋅logP(X,z∣θ)dz 因此 Q ( θ ( t + 1 ) , θ ( t ) ) ≥ Q ( θ ( t ) , θ ( t ) ) Q(\theta^{(t+1)}, \theta^{(t)}) \geq Q(\theta^{(t)}, \theta^{(t)}) Q(θ(t+1),θ(t))≥Q(θ(t),θ(t))
- 证明
H
H
H
H
(
θ
(
t
+
1
)
,
θ
(
t
)
)
−
H
(
θ
(
t
)
,
θ
(
t
)
)
=
∫
z
P
(
z
∣
X
,
θ
(
t
)
)
⋅
log
P
(
X
,
z
∣
θ
(
t
+
1
)
)
d
z
−
∫
z
P
(
z
∣
X
,
θ
(
t
)
)
⋅
log
P
(
X
,
z
∣
θ
(
t
)
)
d
z
=
∫
z
P
(
z
∣
X
,
θ
(
t
)
)
⋅
(
log
P
(
X
,
z
∣
θ
(
t
+
1
)
)
−
log
P
(
X
,
z
∣
θ
(
t
)
)
)
d
z
=
∫
z
P
(
z
∣
X
,
θ
(
t
)
)
⋅
log
P
(
X
,
z
∣
θ
(
t
+
1
)
)
P
(
X
,
z
∣
θ
(
t
)
)
d
z
\begin{align*} H(\theta^{(t+1)}, \theta^{(t)}) - H(\theta^{(t)}, \theta^{(t)}) &= \int_z P(z|X, \theta^{(t)}) \cdot \log P(X, z|\theta^{(t+1)}) dz - \int_z P(z|X, \theta^{(t)}) \cdot \log P(X, z|\theta^{(t)}) dz \\ &= \int_z P(z|X, \theta^{(t)}) \cdot (\log P(X, z|\theta^{(t+1)}) - \log P(X, z|\theta^{(t)})) dz \\ &= \int_z P(z|X, \theta^{(t)}) \cdot \log \frac{P(X, z|\theta^{(t+1)})}{P(X, z|\theta^{(t)})} dz \end{align*}
H(θ(t+1),θ(t))−H(θ(t),θ(t))=∫zP(z∣X,θ(t))⋅logP(X,z∣θ(t+1))dz−∫zP(z∣X,θ(t))⋅logP(X,z∣θ(t))dz=∫zP(z∣X,θ(t))⋅(logP(X,z∣θ(t+1))−logP(X,z∣θ(t)))dz=∫zP(z∣X,θ(t))⋅logP(X,z∣θ(t))P(X,z∣θ(t+1))dz
根据 [ J e n s e n ] ( h t t p s : / / z h u a n l a n . z h i h u . c o m / p / 39315786 ) [Jensen](https://zhuanlan.zhihu.com/p/39315786) [Jensen](https://zhuanlan.zhihu.com/p/39315786)不等式原理： 若 f ( x ) 是 c o n v e x f u n c t i o n （凸函数），则 E [ f ( x ) ] ≥ f ( E [ x ] ) 若f(x) 是convex function（凸函数），则 \mathbb{E}[f(x)] \geq f(\mathbb{E}[x]) 若f(x)是convexfunction（凸函数），则E[f(x)]≥f(E[x])因此： H ( θ ( t + 1 ) , θ ( t ) ) − H ( θ ( t ) , θ ( t ) ) = ∫ z P ( z ∣ X , θ ( t ) ) ⋅ log P ( X , z ∣ θ ( t + 1 ) ) P ( X , z ∣ θ ( t ) ) d z ≤ log ∫ z P ( z ∣ X , θ ( t ) ) ⋅ P ( X , z ∣ θ ( t + 1 ) ) P ( X , z ∣ θ ( t ) ) d z = log ∫ z P ( z ∣ X , θ ( t + 1 ) ) d z = log 1 = 0 \begin{align*} H(\theta^{(t+1)}, \theta^{(t)}) - H(\theta^{(t)}, \theta^{(t)}) &= \int_z P(z|X, \theta^{(t)}) \cdot \log \frac{P(X, z|\theta^{(t+1)})}{P(X, z|\theta^{(t)})} dz \\ &\leq \log \int_z P(z|X, \theta^{(t)}) \cdot \frac{P(X, z|\theta^{(t+1)})}{P(X, z|\theta^{(t)})} dz \\ &= \log \int_z P(z|X, \theta^{(t+1)}) dz \\ &= \log 1 \\ &= 0 \end{align*} H(θ(t+1),θ(t))−H(θ(t),θ(t))=∫zP(z∣X,θ(t))⋅logP(X,z∣θ(t))P(X,z∣θ(t+1))dz≤log∫zP(z∣X,θ(t))⋅P(X,z∣θ(t))P(X,z∣θ(t+1))dz=log∫zP(z∣X,θ(t+1))dz=log1=0综上所述，公式（4）得证。EM算法的收敛性得到证明。EM算法的流程为： - 设置
θ
(
0
)
\theta^{(0)}
θ(0)的初值。不同初值迭代出来的结果可能不同。可以观察下面的示意图，如果
θ
(
t
)
\theta^{(t)}
θ(t)在左边的峰值附近，EM最终就会迭代到左边的局部最优，无法发现右边更大的值，陷入局部最优。（注：该图像来源于:李航老师的《统计学习方法》）
- 更新 θ ( t ) \theta_{(t)} θ(t)。这一步要进行两种计算求期望E，求极大化M。
- 比较 θ t 与 θ ( t + 1 ) \theta^t与\theta^{(t+1)} θt与θ(t+1)的差异，若变化小于一定阈值则结束迭代；否则，返回第二步
4 高斯混合模型的应用
回到 1 极大似然估计的问题背景中，这时候我们不但想知道身高的分布，还想知道体重的分布情况，因此一个数据 x i x_i xi={ 身高、体重 身高、体重 身高、体重}，同时全校的男生、女生混在一起。
4.1 隐变量的引入
我们想一个问题：从全校学生中抽到一个身高为1.7m、体重为60kg的同学的概率是多少？如果明确告诉你这个同学是男生还是女生，那么将其带入高斯模型的中，利用概率密度公式求解，对应的概率就可以直接计算出来。但是，现在你并不知道这个同学的具体性别，这件事情就麻烦了。我们很容易想到用加权的思维去计算这个问题，在这个问题中男生、女生的概率刚好是50%。这里的概率，就是高斯混合模型中的隐变量。隐变量用数学语言做出以下定义： z i j z_{ij} zij表示第 x i x_i xi个数据属于第 j j j个高斯分布的概率 。
4.2 EM-GMM算法推导
- X ： o b s e r v e d d a t a ⟶ X：observed\text{ }data\longrightarrow X：observed data⟶ x 1 , x 2 , . . , x N x_1,x_2,..,x_N x1,x2,..,xN
- Z : l a t e n t d a t a ⟶ Z: latent\text{ }data\longrightarrow Z:latent data⟶ z i j z_{ij} zij
- x : o b s e r v e d v a r i a b l e x:observed\text{ }variable x:observed variable
-
z
:
l
a
t
e
n
t
v
a
r
i
a
b
l
e
z:latent\text{ }variable
z:latent variable
假设高斯混合模型混了 m m m个高斯分布，参数 θ \theta θ={ α 1 , α 2 , . . . , α m , μ 1 , μ 2 , . . . , μ m , Σ 1 , Σ 2 , . . . , Σ m {\alpha_1,\alpha_2,...,\alpha_m,\mu_1,\mu_2,...,\mu_m,\Sigma _1,\Sigma _2,...,\Sigma _m} α1,α2,...,αm,μ1,μ2,...,μm,Σ1,Σ2,...,Σm}，则整个概率密度为 ： P ( x ∣ θ ) = ∑ j = 1 m α j ϕ ( x ∣ μ j , Σ k ) ， w h e r e ∑ j = 1 m α j = 1 P(x|\theta)=\sum_{j=1}^{m} \alpha_j \phi(x|\mu_j,\Sigma_k)，where\sum_{j=1}^{m}\alpha _j=1 P(x∣θ)=j=1∑mαjϕ(x∣μj,Σk)，wherej=1∑mαj=1 对混合分布抽样n次得到 x 1 , . . . , x n {x_1,...,x_n} x1,...,xn，则在第k+1次迭代，待优化式为： max θ Q ( θ , θ k ) = max θ ∑ x ∈ X ∑ z ∈ Z P ( z ∣ y , θ k ) log P ( x , z ∣ θ ) = max θ ∑ x ∈ X ∑ z ∈ Z P ( z , y ∣ θ k ) P ( y ∣ θ k ) log P ( x , z ∣ θ ) = max θ ∑ i = 1 n ∑ j = 1 m α j k ϕ ( x i ∣ θ j k ) ∑ l = 1 m α l k ϕ ( x i ∣ θ l k ) log [ α j ϕ ( x i ∣ θ j ) ] = max θ ∑ i = 1 n ∑ j = 1 m α j k ϕ ( x i ∣ θ j k ) ∑ l = 1 m α l k ϕ ( x i ∣ θ l k ) log [ α j ( 2 π ) D / 2 ∣ Σ j ∣ 1 / 2 exp ( − 1 2 ( x i − μ ) T Σ j − 1 ( x i − μ ) ) ] = max θ ∑ j = 1 m ∑ i = 1 n α j k ϕ ( x i ∣ θ j k ) ∑ l = 1 m α l k ϕ ( x i ∣ θ l k ) [ log α j − 1 2 log Σ j − 1 2 ( x i − μ ) T Σ j − 1 ( x i − μ ) ] \begin{aligned} & \max _\theta Q\left(\theta, \theta^k\right) \\ = & \max _\theta \sum_{x \in {X}} \sum_{z \in{Z}} P\left(z \mid y, \theta^k\right) \log P(x, z \mid \theta) \\ = & \max _\theta \sum_{x \in {X}} \sum_{z \in {Z}} \frac{P\left(z, y \mid \theta^k\right)}{P\left(y \mid \theta^k\right)} \log P(x, z \mid \theta) \\ = & \max _\theta \sum_{i=1}^n \sum_{j=1}^m \frac{\alpha_j^k \phi\left(x_i \mid \theta_j^k\right)}{\sum_{l=1}^m \alpha_l^k \phi\left(x_i \mid \theta_l^k\right)} \log \left[\alpha_j \phi\left(x_i \mid \theta_j\right)\right] \\ = & \max _\theta \sum_{i=1}^n \sum_{j=1}^m \frac{\alpha_j^k \phi\left(x_i \mid \theta_j^k\right)}{\sum_{l=1}^m \alpha_l^k \phi\left(x_i \mid \theta_l^k\right)} \log \left [ \frac{\alpha_j}{(2 \pi)^{D / 2}|\Sigma_j|^{1 / 2}} \exp \left(-\frac{1}{2}(x_i-\mu)^T \Sigma_j^{-1}(x_i-\mu)\right) \right ] \\ = & \max _\theta \sum_{j=1}^m \sum_{i=1}^n \frac{\alpha_j^k \phi\left(x_i \mid \theta_j^k\right)}{\sum_{l=1}^m \alpha_l^k \phi\left(x_i \mid \theta_l^k\right)}\left[\log \alpha_j-\frac{1}{2}\log \Sigma_j-\frac{1}{2}(x_i-\mu)^T \Sigma_j^{-1}(x_i-\mu)\right] \end{aligned} =====θmaxQ(θ,θk)θmaxx∈X∑z∈Z∑P(z∣y,θk)logP(x,z∣θ)θmaxx∈X∑z∈Z∑P(y∣θk)P(z,y∣θk)logP(x,z∣θ)θmaxi=1∑nj=1∑m∑l=1mαlkϕ(xi∣θlk)αjkϕ(xi∣θjk)log[αjϕ(xi∣θj)]θmaxi=1∑nj=1∑m∑l=1mαlkϕ(xi∣θlk)αjkϕ(xi∣θjk)log[(2π)D/2∣Σj∣1/2αjexp(−21(xi−μ)TΣj−1(xi−μ))]θmaxj=1∑mi=1∑n∑l=1mαlkϕ(xi∣θlk)αjkϕ(xi∣θjk)[logαj−21logΣj−21(xi−μ)TΣj−1(xi−μ)]
4.2.1 求解 α \alpha α
记
p
j
=
α
j
k
ϕ
(
x
i
∣
θ
j
k
)
∑
l
=
1
m
α
l
k
ϕ
(
x
i
∣
θ
l
k
)
p_j= \frac{\alpha_j^k \phi\left(x_i \mid \theta_j^k\right)}{\sum_{l=1}^m \alpha_l^k \phi\left(x_i \mid \theta_l^k\right)}
pj=∑l=1mαlkϕ(xi∣θlk)αjkϕ(xi∣θjk)
构造拉格朗日方程
L
(
p
,
λ
)
=
∑
j
=
1
m
∑
i
=
1
n
p
j
∗
log
α
j
−
λ
(
∑
j
=
1
m
p
j
−
1
)
L(p,\lambda )= \sum_{j=1}^m \sum_{i=1}^np_j*\log \alpha_j-\lambda(\sum_{j=1}^{m}p_j-1)
L(p,λ)=j=1∑mi=1∑npj∗logαj−λ(j=1∑mpj−1)
对
p
j
p_j
pj求导：
∂
L
∂
p
j
=
∑
i
=
1
n
1
α
j
∗
p
j
−
λ
=
0
\frac{\partial L}{\partial p_j}= \sum_{i=1}^n\frac{1}{\alpha_j}*p_j-\lambda=0
∂pj∂L=i=1∑nαj1∗pj−λ=0
∂
L
∂
p
j
=
∑
i
=
1
n
p
j
−
α
j
∗
λ
=
0
\frac{\partial L}{\partial p_j}= \sum_{i=1}^np_j-\alpha_j*\lambda=0
∂pj∂L=i=1∑npj−αj∗λ=0对于所有的
j
j
j，即
j
j
j从
1
1
1到
m
m
m，对上式加和：
∑
j
=
1
m
∑
i
=
1
n
p
j
−
∑
i
=
j
m
α
j
∗
λ
=
0
**
\sum_{j=1}^m\sum_{i=1}^np_j-\sum_{i=j}^m\alpha_j*\lambda=0\text{ **}
j=1∑mi=1∑npj−i=j∑mαj∗λ=0 **
因为
∑
j
=
1
m
p
j
=
1
,
∑
j
=
1
m
α
j
=
1
\sum_{j=1}^mp_j=1,\sum_{j=1}^m\alpha_j=1
j=1∑mpj=1,j=1∑mαj=1
所以推出
λ
=
N
\lambda=N
λ=N,代入(**)式子，得
α
k
=
1
n
∑
i
=
1
n
α
j
k
ϕ
(
x
i
∣
θ
j
k
)
∑
l
=
1
m
α
l
k
ϕ
(
x
i
∣
θ
l
k
)
\alpha_k=\frac{1}{n}\sum_{i=1}^{n}\frac{\alpha_j^k \phi\left(x_i \mid \theta_j^k\right)}{\sum_{l=1}^m \alpha_l^k \phi\left(x_i \mid \theta_l^k\right)}
αk=n1i=1∑n∑l=1mαlkϕ(xi∣θlk)αjkϕ(xi∣θjk)
4.2.2 求解 μ 、 Σ \mu、\Sigma μ、Σ
对这两项的求导不进行展开，较为简单。只展示最终结果：
μ j ← ∑ i = 1 N p j x ( i ) ∑ i = 1 N p j Σ j ← ∑ i = 1 N p j ⋅ { ( x ( i ) − μ j ) ( x ( i ) − μ j ) T } ∑ i = 1 N p j \begin{aligned} & \mu_j \leftarrow \frac{\sum_{i=1}^Np_j x^{(i)}}{\sum_{i=1}^Np_j} \quad \Sigma_j \leftarrow \frac{\sum_{i=1}^N p_j \cdot\left\{\left(x^{(i)}-\mu_j\right)\left(x^{(i)}-\mu_j\right)^T\right\}}{\sum_{i=1}^Np_j} \end{aligned} μj←∑i=1Npj∑i=1Npjx(i)Σj←∑i=1Npj∑i=1Npj⋅{(x(i)−μj)(x(i)−μj)T} p j p_j pj在4.2.1第一个公式进行说明
python实现代码
import numpy as np
import matplotlib.pyplot as plt
from sklearn.mixture import GaussianMixture
# 初始化分布参数
MU1 = np.array([1, 2])
SIGMA1 = np.array([[1, 0], [0, 0.5]])
MU2 = np.array([-1, -1])
SIGMA2 = np.array([[1, 0], [0, 1]])
# 生成数据点
data1 = np.random.multivariate_normal(MU1, SIGMA1, 1000)
data2 = np.random.multivariate_normal(MU2, SIGMA2, 1000)
X = np.concatenate([data1, data2])
#对X进行随机打乱，此步复现时不可忽略
np.random.shuffle(X)
# Step 2：请在这里绘制二维散点图。
plt.scatter(data1[:, 0], data1[:, 1], label='Class 1')
plt.scatter(data2[:, 0], data2[:, 1], label='Class 2')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Scatter Plot of Data Points')
plt.legend()
plt.show()
# Step 3：GMMs的代码实现。
def my_fit_GMM(X, K, Max_iters):
N, D = X.shape # 获取数据X的维度
alpha = np.ones(K) / K # 初始化混合系数
mu = X[np.random.choice(N, K, False), :]
# 初始化COV
cov = [np.cov(X.T) for _ in range(K)]
E = np.zeros((N, K)) # 初始化后验概率矩阵
while Max_iters>0:
Max_iters-=1
oldmu = mu.copy() # E-Step
for i in range(N):
for j in range(K):
E[i, j] = alpha[j] * multivariate_gaussian_pdf(X[i], mu[j], cov[j])
E /= E.sum(axis=1, keepdims=True)
# M-Step
sum_E = E.sum(axis=0)
alpha = sum_E / N
for j in range(K):
mu[j] = E[:,j]@X / sum_E[j]
x_mu = X - mu[j]
cov[j] = (E[:, j, np.newaxis] * x_mu).T @ x_mu / sum_E[j]
if abs(np.linalg.det(mu-oldmu))<2.2204e-16:
break
return mu, cov,alpha
# 混合高斯的pdf
def multivariate_gaussian_pdf(x, mu, cov):
k = len(mu)
cov_det = np.linalg.det(cov)
cov_inv = np.linalg.inv(cov)
pdf = (1.0 / np.sqrt((2 * np.pi) ** k * cov_det))*np.exp(-0.5 * (x - mu) @ cov_inv @ (x - mu).T) # 高斯的公式
return pdf
# Step 4：请用GMMs拟合散点分布。
means, cov, weights = my_fit_GMM(X,2,100)
# Step 5：绘制GMMs所拟合分布的概率密度函数。画图函数
from scipy.stats import multivariate_normal
def draw_results(m,s,w):
# 定义两个二维高斯分布的参数
m1 = m[0]
s1 = s[0]
m2 = m[1]
s2 = s[1]
# 生成网格点
x, y = np.mgrid[-5:5:.01, -5:5:.01]
pos = np.dstack((x, y))
# 计算每个点的概率密度值
rv1 = multivariate_normal(m1, s1)
rv2 = multivariate_normal(m2, s2)
z1 = rv1.pdf(pos)
z2 = rv2.pdf(pos)
# 混合两个高斯分布的概率密度函数
z = w[0] * z1 + w[1] * z2 # 设置混合比例
# 绘制3D图
fig = plt.figure(figsize=(5, 10))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, z, cmap='viridis')
# 设置坐标轴标签和标题
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Probability Density')
ax.set_title('3D Plot of Mixture Gaussian Probability Density')
plt.show()
# Step 6：输出估计的均值和方差。 可以与Step 1的数据进行对比
print("Data1.Means:")
print(means[0])
print("Data2.Means:")
print(means[1])
print("Data1.Covariances:")
print(cov[0])
print("Data2.Covariances:")
print(cov[1])
print("Weights:")
print(weights)
draw_results(means, cov, weights)
参考链接
https://blog.csdn.net/zouxy09/article/details/8537620
https://www.cnblogs.com/qizhou/p/13100817.html
https://www.bilibili.com/video/BV13b411w7Xj/?p=4&vd_source=395b52d7c4e90a9c96a83181871f36bb