本篇文章是之前期望极大算法(EM算法)文章的后续,有需要可以先看看那篇文章关于EM算法的推导。
高斯混合模型
高斯混合模型是研究算法的人避不开的一个东西,其在非深度学习的远古时代经常被用到,比如图像处理任务的前背景提取,点云处理任务的点云聚类等等等等。
具体的,高斯混合模型是指具有如下形式的概率分布模型:
P
(
y
∣
θ
)
=
∑
k
=
1
K
α
k
ϕ
(
y
∣
θ
k
)
P(y \mid \theta)=\sum_{k=1}^{K} \alpha_{k} \phi\left(y \mid \theta_{k}\right)
P(y∣θ)=k=1∑Kαkϕ(y∣θk)
其中,
α
k
\alpha_{k}
αk是系数,
α
k
⩾
0
\alpha_{k} \geqslant 0
αk⩾0,
∑
k
=
1
K
α
k
=
1
\quad \sum_{k=1}^{K} \alpha_{k}=1
∑k=1Kαk=1;
ϕ
(
y
∣
θ
k
)
\phi\left(y \mid \theta_{k}\right)
ϕ(y∣θk)是高斯分布密度,
θ
k
=
(
μ
k
,
σ
k
2
)
\theta_{k}=\left(\mu_{k}, \sigma_{k}^{2}\right)
θk=(μk,σk2),
ϕ
(
y
∣
θ
k
)
=
1
2
π
σ
k
exp
(
−
(
y
−
μ
k
)
2
2
σ
k
2
)
\phi\left(y \mid \theta_{k}\right)=\frac{1}{\sqrt{2 \pi} \sigma_{k}} \exp \left(-\frac{\left(y-\mu_{k}\right)^{2}}{2 \sigma_{k}^{2}}\right)
ϕ(y∣θk)=2πσk1exp(−2σk2(y−μk)2)
称为第
k
k
k个分模型。
Q Q Q函数的一般表达
在算法处理的过程中,将问题建模成高斯混合模型后,往往需要去解模型中的参数,这个时候就需要用到EM算法。
在期望极大算法(EM算法)文章的分析中我们已经知道,要进行EM算法,得先得到
Q
Q
Q函数:
Q
(
θ
,
θ
(
i
)
)
=
∑
Z
log
P
(
Y
,
Z
∣
θ
)
P
(
Z
∣
Y
,
θ
(
i
)
)
Q\left(\theta, \theta^{(i)}\right)=\sum_{Z} \log P(Y, Z \mid \theta)P\left(Z \mid Y, \theta^{(i)}\right)
Q(θ,θ(i))=Z∑logP(Y,Z∣θ)P(Z∣Y,θ(i))
在随后的抛硬币问题中我们也介绍了一种求解这个
Q
Q
Q函数的方法.但如果看过李老师的书的人会发现,本人没有用书中给出的看起来更具总结性的
Q
Q
Q函数形式来求解问题,原因在于当时本人也不明白那个式子是怎么来的。。。
在李航老师的书上,
Q
Q
Q函数是这样定义的:完全数据的对数似然函数
log
P
(
Y
,
Z
∣
θ
)
\log P(Y, Z \mid \theta)
logP(Y,Z∣θ)关于在给定观测数据
Y
Y
Y和当前参数
θ
(
i
)
\theta^{(i)}
θ(i)下对未观测数据
Z
Z
Z的条件概率分布
P
(
Z
∣
Y
,
θ
(
i
)
)
P\left(Z \mid Y, \theta^{(i)}\right)
P(Z∣Y,θ(i))的期望称为
Q
Q
Q函数,即:
Q
(
θ
,
θ
(
i
)
)
=
E
Z
[
log
P
(
Y
,
Z
∣
θ
)
∣
Y
,
θ
(
i
)
]
Q\left(\theta, \theta^{(i)}\right)=E_{Z}\left[\log P(Y, Z \mid \theta) \mid Y, \theta^{(i)}\right]
Q(θ,θ(i))=EZ[logP(Y,Z∣θ)∣Y,θ(i)]
第一次看到这个定义和下面的公式感觉整个人都不好了!实在不知道他们之间为什么是个相等的关系。在经过一两个小时的发呆、无助、掉头发后突然看懂了!
假设
log
P
(
Y
,
Z
∣
θ
)
\log P(Y, Z \mid \theta)
logP(Y,Z∣θ)是只与变量
Z
Z
Z相关的函数,则可以把其写成
f
(
Z
)
f(Z)
f(Z),当
Z
Z
Z取得一个定值的时候,其就是一个固定的数值。如果这个时候对它取期望,就有:
E
Z
(
f
(
Z
)
)
=
∑
Z
[
f
(
Z
)
P
(
Z
)
]
E_Z(f(Z))=\sum_{Z}\left[f(Z)P\left(Z\right)\right]
EZ(f(Z))=Z∑[f(Z)P(Z)]
而如果Z的取值本身受到别的参数
x
x
x,
y
y
y影响,而这些参数都已经给出,则原式可以写成:
E
(
f
(
Z
)
∣
x
,
y
)
=
∑
Z
[
f
(
Z
)
P
(
Z
∣
x
,
y
)
]
E(f(Z)\mid x,y)=\sum_{Z}\left[f(Z)P\left(Z\mid x,y\right )\right]
E(f(Z)∣x,y)=Z∑[f(Z)P(Z∣x,y)]
代入我们已知的量,等式得证。
由证明的过程也可以知道,这里的给定观测数据 Y Y Y和当前参数两个已知量影响的只有 Z Z Z这一变量。由于他们与 log P ( Y , Z ∣ θ ) \log P(Y, Z \mid \theta) logP(Y,Z∣θ)中的两个参数看起来似乎有关系,因此才大大增加了理解的难度。
有了 Q Q Q函数的一般表达,从式子的形式我们知道关键的一步是把高斯混合模型的 log P ( Y , Z ∣ θ ) \log P(Y, Z \mid \theta) logP(Y,Z∣θ)给列出来,也就是把其完全数据的似然函数的对数列出来。
高斯混合模型参数估计的EM算法
假设数据
y
1
,
y
2
,
⋯
,
y
N
y_{1},y_{2}, \cdots, y_{N}
y1,y2,⋯,yN由高斯混合模型生成,
P
(
y
∣
θ
)
=
∑
k
=
1
K
α
k
ϕ
(
y
∣
θ
k
)
P(y \mid \theta)=\sum_{k=1}^{K} \alpha_{k} \phi\left(y \mid \theta_{k}\right)
P(y∣θ)=k=1∑Kαkϕ(y∣θk)
其中,
θ
=
(
α
1
,
α
2
,
⋯
,
α
K
;
θ
1
,
θ
2
,
⋯
,
θ
K
)
\theta=\left(\alpha_{1}, \alpha_{2}, \cdots, \alpha_{K} ; \theta_{1}, \theta_{2}, \cdots, \theta_{K}\right)
θ=(α1,α2,⋯,αK;θ1,θ2,⋯,θK),求高斯混合模型的参数,就是用EM算法估计高斯混合模型的参数
θ
\theta
θ。
回顾观测数据 y j , j = 1 , 2 , ⋯ , N , y_{j},j=1,2, \cdots, N, yj,j=1,2,⋯,N, 的产生过程:
- 依概率 α k \alpha_{k} αk选择第 k k k个高斯分布分模型 ϕ ( y ∣ θ k ) \phi\left(y \mid \theta_{k}\right) ϕ(y∣θk);
- 依第 k k k个分模型的概率分布 ϕ ( y ∣ θ k ) \phi\left(y \mid \theta_{k}\right) ϕ(y∣θk)生成观测数据 y j y_{j} yj。
在高斯混合模型建模中,结果通常是已知的—观测数据 y j , j = 1 , 2 , ⋯ , N , y_{j},j=1,2, \cdots, N, yj,j=1,2,⋯,N, 是已知的;
而结果产生自哪个分模型未知—反映观测数据
y
j
y_{j}
yj来自第
k
k
k个分模型的数据未知,
k
=
1
,
2
,
⋯
,
K
,
k=1,2, \cdots, K ,
k=1,2,⋯,K,以隐变量
γ
j
k
\gamma_{j k}
γjk表示,可定义如下:
γ
j
k
=
{
1
,
第
j
个观测来自第
k
个分模型
0
,
否则
\gamma_{j k}=\left\{\begin{array}{ll} 1, & \text { 第 } j \text { 个观测来自第 } k \text { 个分模型 } \\ 0, & \text { 否则 } \end{array}\right.
γjk={1,0, 第 j 个观测来自第 k 个分模型 否则
j = 1 , 2 , ⋯ , N ; k = 1 , 2 , ⋯ , K j=1,2, \cdots, N ; \quad k=1,2, \cdots, K j=1,2,⋯,N;k=1,2,⋯,K
有了观测数据
y
j
y_j
yj及未观测数据
γ
j
k
\gamma_{j k}
γjk,那么完全数据是
(
y
j
,
γ
j
1
,
γ
j
2
,
⋯
,
γ
j
K
)
,
j
=
1
,
2
,
⋯
,
N
\left(y_{j}, \gamma_{j 1}, \gamma_{j 2}, \cdots, \gamma_{j K}\right), \quad j=1,2, \cdots, N
(yj,γj1,γj2,⋯,γjK),j=1,2,⋯,N
于是,可以写出完全数据的似然函数:
P
(
y
,
γ
∣
θ
)
=
∏
j
=
1
N
P
(
y
j
,
γ
j
1
,
γ
j
2
,
⋯
,
γ
j
K
∣
θ
)
=
∏
k
=
1
K
∏
j
=
1
N
[
α
k
ϕ
(
y
j
∣
θ
k
)
]
γ
j
k
=
∏
k
=
1
K
α
k
n
∏
j
=
1
N
[
ϕ
(
y
j
∣
θ
k
)
]
γ
j
k
=
∏
k
=
1
K
α
k
n
k
∏
j
=
1
N
[
1
2
π
σ
k
exp
(
−
(
y
j
−
μ
k
)
2
2
σ
k
2
)
]
γ
j
k
\begin{array}{l} P(y, \gamma \mid \theta)=\prod_{j=1}^{N} P\left(y_{j}, \gamma_{j 1}, \gamma_{j 2}, \cdots, \gamma_{j K} \mid \theta\right)\\ =\prod_{k=1}^{K} \prod_{j=1}^{N}\left[\alpha_{k} \phi\left(y_{j} \mid \theta_{k}\right)\right]^{\gamma_{jk}} \\ =\prod_{k=1}^{K} \alpha_{k}^{n} \prod_{j=1}^{N}\left[\phi\left(y_{j} \mid \theta_{k}\right)\right]^{\gamma_{jk}} \\ =\prod_{k=1}^{K} \alpha_{k}^{n_{k}} \prod_{j=1}^{N}\left[\frac{1}{\sqrt{2 \pi} \sigma_{k}} \exp \left(-\frac{\left(y_{j}-\mu_{k}\right)^{2}}{2 \sigma_{k}^{2}}\right)\right]^{\gamma_{jk}} \end{array}
P(y,γ∣θ)=∏j=1NP(yj,γj1,γj2,⋯,γjK∣θ)=∏k=1K∏j=1N[αkϕ(yj∣θk)]γjk=∏k=1Kαkn∏j=1N[ϕ(yj∣θk)]γjk=∏k=1Kαknk∏j=1N[2πσk1exp(−2σk2(yj−μk)2)]γjk
式中,
n
k
=
∑
j
=
1
N
γ
j
k
n_{k}=\sum_{j=1}^{N} \gamma_{j k}
nk=∑j=1Nγjk,
∑
k
=
1
K
n
k
=
N
\sum_{k=1}^{K} n_{k}=N
∑k=1Knk=N。
那么,完全数据的对数似然函数为:
log
P
(
y
,
γ
∣
θ
)
=
∑
k
=
1
K
n
k
log
α
k
+
∑
j
=
1
N
γ
j
k
[
log
(
1
2
π
σ
k
)
−
(
y
j
−
μ
k
)
2
2
σ
k
2
]
=
∑
k
=
1
K
n
k
log
α
k
+
∑
j
=
1
N
γ
j
k
[
log
(
1
2
π
)
−
log
σ
k
−
(
y
j
−
μ
k
)
2
2
σ
k
2
]
\log P(y, \gamma \mid \theta)=\sum_{k=1}^{K} n_{k} \log \alpha_{k}+\sum_{j=1}^{N} \gamma_{j k}\left[\log \left(\frac{1}{\sqrt{2 \pi}\sigma_{k}}\right)-\frac{\left(y_{j}-\mu_{k}\right)^{2}}{2 \sigma_{k}^{2}}\right]\\ =\sum_{k=1}^{K} n_{k} \log \alpha_{k}+\sum_{j=1}^{N} \gamma_{j k}\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log{\sigma_{k}-}\frac{\left(y_{j}-\mu_{k}\right)^{2}}{2 \sigma_{k}^{2}}\right]\\
logP(y,γ∣θ)=k=1∑Knklogαk+j=1∑Nγjk[log(2πσk1)−2σk2(yj−μk)2]=k=1∑Knklogαk+j=1∑Nγjk[log(2π1)−logσk−2σk2(yj−μk)2]
整个
Q
Q
Q函数为:
Q
(
θ
,
θ
(
i
)
)
=
E
[
log
P
(
y
,
γ
∣
θ
)
∣
y
,
θ
(
i
)
]
=
E
{
∑
k
=
1
K
n
k
log
α
k
+
∑
j
=
1
N
γ
j
k
[
log
(
1
2
π
)
−
log
σ
k
−
(
y
j
−
μ
k
)
2
2
σ
k
2
]
}
=
∑
k
=
1
K
{
∑
j
=
1
N
(
E
γ
j
k
)
log
α
k
+
∑
j
=
1
N
(
E
γ
j
k
)
[
log
(
1
2
π
)
−
log
σ
k
−
1
2
σ
k
2
(
y
j
−
μ
k
)
2
]
}
\begin{aligned} Q\left(\theta, \theta^{(i)}\right) &=E\left[\log P(y, \gamma \mid \theta) \mid y, \theta^{(i)}\right] \\ &=E\left\{\sum_{k=1}^{K} n_{k} \log \alpha_{k}+\sum_{j=1}^{N} \gamma_{j k}\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log \sigma_{k}-\frac{\left(y_{j}-\mu_{k}\right)^{2}}{2 \sigma_{k}^{2}}\right]\right\} \\ &=\sum_{k=1}^{K}\left\{\sum_{j=1}^{N}\left(E \gamma_{j k}\right) \log \alpha_{k}+\sum_{j=1}^{N}\left(E \gamma_{j k}\right)\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log \sigma_{k}-\frac{1}{2 \sigma_{k}^{2}}\left(y_{j}-\mu_{k}\right)^{2}\right]\right\} \end{aligned}
Q(θ,θ(i))=E[logP(y,γ∣θ)∣y,θ(i)]=E{k=1∑Knklogαk+j=1∑Nγjk[log(2π1)−logσk−2σk2(yj−μk)2]}=k=1∑K{j=1∑N(Eγjk)logαk+j=1∑N(Eγjk)[log(2π1)−logσk−2σk21(yj−μk)2]}
这里出现了
E
γ
j
k
E \gamma_{j k}
Eγjk,依据我们前面的推导它是已知观测数据和当前参数的情况下,隐函数的似然。可以写成:
E
(
γ
j
k
∣
y
,
θ
)
E\left(\gamma_{j k} \mid y, \theta\right)
E(γjk∣y,θ),记为
γ
^
j
k
\hat{\gamma}_{j k}
γ^jk。有
γ
^
j
k
=
E
(
γ
j
k
∣
y
,
θ
)
=
P
(
γ
j
k
=
1
∣
y
,
θ
)
=
P
(
γ
j
k
=
1
,
y
j
∣
θ
)
∑
k
=
1
K
P
(
γ
j
k
=
1
,
y
j
∣
θ
)
=
P
(
y
j
∣
γ
j
k
=
1
,
θ
)
P
(
γ
j
k
=
1
∣
θ
)
∑
k
=
1
K
P
(
y
j
∣
γ
j
k
=
1
,
θ
)
P
(
γ
j
k
=
1
∣
θ
)
=
α
k
ϕ
(
y
j
∣
θ
k
)
∑
k
=
1
K
α
k
ϕ
(
y
j
∣
θ
k
)
,
j
=
1
,
2
,
⋯
,
N
;
k
=
1
,
2
,
⋯
,
K
\begin{aligned} \hat{\gamma}_{j k} &=E\left(\gamma_{j k} \mid y, \theta\right)=P\left(\gamma_{j k}=1 \mid y, \theta\right) \\ &=\frac{P\left(\gamma_{j k}=1, y_{j} \mid \theta\right)}{\sum_{k=1}^{K} P\left(\gamma_{j k}=1, y_{j} \mid \theta\right)} \\ &=\frac{P\left(y_{j} \mid \gamma_{j k}=1, \theta\right) P\left(\gamma_{j k}=1 \mid \theta\right)}{\sum_{k=1}^{K} P\left(y_{j} \mid \gamma_{j k}=1, \theta\right) P\left(\gamma_{j k}=1 \mid \theta\right)} \\ &=\frac{\alpha_{k} \phi\left(y_{j} \mid \theta_{k}\right)}{\sum_{k=1}^{K} \alpha_{k} \phi\left(y_{j} \mid \theta_{k}\right)}, \quad j=1,2, \cdots, N ; \quad k=1,2, \cdots, K \end{aligned}
γ^jk=E(γjk∣y,θ)=P(γjk=1∣y,θ)=∑k=1KP(γjk=1,yj∣θ)P(γjk=1,yj∣θ)=∑k=1KP(yj∣γjk=1,θ)P(γjk=1∣θ)P(yj∣γjk=1,θ)P(γjk=1∣θ)=∑k=1Kαkϕ(yj∣θk)αkϕ(yj∣θk),j=1,2,⋯,N;k=1,2,⋯,K
推到这可以知道
γ
^
j
k
\hat{\gamma}_{j k}
γ^jk等于当前模型参数下第
j
j
j个观测数据来自第
k
k
k个分模型的概率,称为分模型
k
k
k对观测数据
y
j
y_j
yj的响应度。代入
Q
Q
Q函数可以得到:
Q
(
θ
,
θ
(
i
)
)
=
∑
k
=
1
K
{
n
k
log
α
k
+
∑
j
=
1
N
γ
^
j
k
[
log
(
1
2
π
)
−
log
σ
k
−
1
2
σ
k
2
(
y
j
−
μ
k
)
2
]
}
Q\left(\theta, \theta^{(i)}\right) =\sum_{k=1}^{K}\left\{n_{k} \log \alpha_{k}+\sum_{j=1}^{N}\hat{\gamma}_{j k}\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log \sigma_{k}-\frac{1}{2 \sigma_{k}^{2}}\left(y_{j}-\mu_{k}\right)^{2}\right]\right\}
Q(θ,θ(i))=k=1∑K{nklogαk+j=1∑Nγ^jk[log(2π1)−logσk−2σk21(yj−μk)2]}
到此就得到了只含有模型参数的
Q
Q
Q函数,真的不容易啊!要是没有李航老师的书做参考,估计推到吐血都推不出来!
有了
Q
Q
Q函数相当于
E
E
E步就有了,
M
M
M步很简单,就是
Q
Q
Q函数取相应模型参数的偏导然后求其极值点也就是等于0的点即可。求得结果如下:
μ
^
k
=
∑
j
=
1
N
γ
^
j
k
y
j
∑
j
=
1
N
γ
^
j
k
,
k
=
1
,
2
,
⋯
,
K
\hat{\mu}_{k}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k} y_{j}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2, \cdots, K
μ^k=∑j=1Nγ^jk∑j=1Nγ^jkyj,k=1,2,⋯,K
σ ^ k 2 = ∑ j = 1 N γ ^ j k ( y j − μ k ) 2 ∑ j = 1 N γ ^ j k , k = 1 , 2 , ⋯ , K \hat{\sigma}_{k}^{2}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}\left(y_{j}-\mu_{k}\right)^{2}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2, \cdots, K σ^k2=∑j=1Nγ^jk∑j=1Nγ^jk(yj−μk)2,k=1,2,⋯,K
α ^ k = n k N = ∑ j = 1 N γ ^ j k N , k = 1 , 2 , ⋯ , K \hat{\alpha}_{k}=\frac{n_{k}}{N}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}}{N}, \quad k=1,2, \cdots, K α^k=Nnk=N∑j=1Nγ^jk,k=1,2,⋯,K
可以看到其只包含有 γ ^ j k \hat{\gamma}_{j k} γ^jk这一与当前参数相关的变量,因此在实际计算过程中,我们只需要设定初值,然后在 E E E步计算出 γ ^ j k \hat{\gamma}_{j k} γ^jk,在M步将计算得到的 γ ^ j k \hat{\gamma}_{j k} γ^jk代入求得新的参数,一直重复直到收敛即可。
真的是推导要老命,编程3分钟啊!
高斯混合模型使用EM算法解决实际问题
已知观测数据 -67,-48,6,8,14,16,23,24,28,29,41,49,56,60,75 试估计两个分量的高斯混合模型的5个参数。
由上面的推导已经知道了所有需要的信息,因此只要设定初值然后直接代公式即可,代码如下:
#include <iostream>
#include <cmath>
#define N 15
#define pi 3.1415926535898
class theta
{
public:
double alpha;
double mu;
double sigma;
void print()
{
std::cout << "--------------" << std::endl;
std::cout << "alpha:" << alpha << std::endl;
std::cout << "mu:" << mu << std::endl;
std::cout << "sigma:" << sigma << std::endl;
std::cout << "sigma平方" << sigma* sigma << std::endl;
}
};
theta mtheta[2];//两个高斯分模型参数
double gamma[2][N];//E步结果gamma
double y[N] = { -67,-48,6,8,14,16,23,24,28,29,41,49,56,60,75 };//观测结果
double phi(theta& mtheta, double yj)
{
return 1 / (sqrt(2 * pi) * mtheta.sigma) * exp(-pow((yj - mtheta.mu), 2) / (2 * pow(mtheta.sigma, 2)));
}
void EStep()
{
for (size_t j = 0; j < N; j++)
{
gamma[0][j] = mtheta[0].alpha * phi(mtheta[0], y[j]);
gamma[1][j] = mtheta[1].alpha * phi(mtheta[1], y[j]);
double sum = gamma[0][j] + gamma[1][j];
gamma[0][j] = gamma[0][j] / sum;
gamma[1][j] = gamma[1][j] / sum;
//std::cout << "gamma0:" << gamma[0][j] << std::endl;
//std::cout << "gamma1:" << gamma[1][j] << std::endl;
}
}
void MStep()
{
for (size_t k = 0; k < 2; k++)
{
double mu = 0;
double sigma = 0;
double sumgamma = 0;
for (size_t j = 0; j < N; j++)
{
mu += gamma[k][j] * y[j];
sigma += gamma[k][j] * pow((y[j] - mtheta[k].mu), 2);
sumgamma += gamma[k][j];
}
mtheta[k].mu = mu / sumgamma;
mtheta[k].sigma = sqrt(sigma / sumgamma);
mtheta[k].alpha = sumgamma / N;
}
}
int main()
{
//初始化高斯分模型参数变量
mtheta[0].alpha = 0.5;
mtheta[0].mu = 30;
mtheta[0].sigma = sqrt(500);
mtheta[1].alpha = 0.5;
mtheta[1].mu = -30;
mtheta[1].sigma = sqrt(500);
for (size_t k = 0; k < 2; k++)
{
mtheta[k].print();
}
//迭代10次
for (size_t i = 0; i < 10; i++)
{
EStep();
MStep();
for (size_t k = 0; k < 2; k++)
{
mtheta[k].print();
}
}
}
如果求出了模型参数后想知道观察结果在这一参数下的分类,再求一次 γ ^ j k \hat{\gamma}_{j k} γ^jk即可!
参考文章
https://zhuanlan.zhihu.com/p/32508410