Gaussian Mixed Model
应用
-
聚类
K-means无法处理两个聚类中心点相同的类。比如 A ∼ N ( μ , σ 1 2 ) , B ∼ N ( μ , σ 2 2 ) A\sim N(\mu,\;\sigma_1^2),\;B\sim N(\mu,\sigma_2^2) A∼N(μ,σ12),B∼N(μ,σ22) 是无法用k-means进行聚类的。
-
密度估计
-
新数据的生成
原理
我们认为数据空间是由某些高斯分布生成的，但对于某一具体的样本单元，我们只能观测到它的数值，而无法观测到它具体是由哪个高斯分布生成的，所以我们对样本单元背后的高斯分布进行建模的时候，认为样本单元以一定的权重（也可以理解为概率）服从某个具体的高斯分布。
符号说明
| 符号 | 解释 | |
|---|---|---|
| n n n | 高斯分布的个数（总体聚类的类个数） | |
| N N N | 样本数 | |
| x i x_i xi | i = 1 , 2 , … , N i=1,2,\dots,N i=1,2,…,N | 第 i i i个样本单元 |
| z i z_i zi | i = 1 , 2 , … , N i=1,2,\dots,N i=1,2,…,N | 第 i i i个样本单元的隐变量 |
| c j c_j cj | j = 1 , 2 , … , n j=1,2,\dots,n j=1,2,…,n | z i z_i zi可能的取值 |
| π i j \pi_{ij} πij | i = 1 , 2 , … , N j = 1 , 2 , … , n i=1,2,\dots,N\\j=1,2,\dots,n i=1,2,…,Nj=1,2,…,n | p ( z i = c j ) p(z_i=c_j) p(zi=cj) |
| π i j ( t ) \pi_{ij}^{(t)} πij(t) | i = 1 , 2 , … , N , j = 1 , 2 , … , n , t ∈ N + i=1,2,\dots,N,\\j=1,2,\dots,n,\\t\in\mathbb{N}^+ i=1,2,…,N,j=1,2,…,n,t∈N+ | 隐变量 z i z_i zi在第 t t t步时取值为 c j c_j cj的先验概率 |
| γ t ( z i = c j ) \gamma_t(z_i=c_j) γt(zi=cj) | i = 1 , 2 , … , N , j = 1 , 2 , … , n , t ∈ N + i=1,2,\dots,N,\\j=1,2,\dots,n,\\t\in\mathbb{N}^+ i=1,2,…,N,j=1,2,…,n,t∈N+ | 隐变量 z i z_i zi在第 t t t步时取值为 c j c_j cj的后验概率 |
| π i = ( π i 1 , π i 2 , … , π i n ) \pi_i=(\pi_{i1},\pi_{i2},\dots,\pi_{in}) πi=(πi1,πi2,…,πin) | i = 1 , 2 , … , N i=1,2,\dots,N i=1,2,…,N | z i z_i zi的分布 |
| 30 | ||
| μ j \mu_j μj | j = 1 , 2 , … , n j=1,2,\dots,n j=1,2,…,n | 第 j j j个高斯分布的均值 |
| μ j ( t ) \mu_j^{(t)} μj(t) | j = 1 , 2 , … , n t ∈ N + j=1,2,\dots,n\\t\in\mathbb{N}^+ j=1,2,…,nt∈N+ | 第 t t t步时，第 j j j个高斯分布的均值 |
| μ = ( μ 1 , μ 2 , … , μ n ) \mu=(\mu_1,\mu_2,\dots,\mu_n) μ=(μ1,μ2,…,μn) | 所有高斯分布的均值 | |
| Σ j \Sigma_j Σj | j = 1 , 2 , … , n j=1,2,\dots,n j=1,2,…,n | 第 j j j个高斯分布的协方差矩阵 |
| Σ j ( t ) \Sigma_j^{(t)} Σj(t) | j = 1 , 2 , … , n t ∈ N + j=1,2,\dots,n\\t\in\mathbb{N}^+ j=1,2,…,nt∈N+ | 第 t t t步时，第 j j j个高斯分布的协方差矩阵 |
| Σ = ( Σ 1 , Σ 2 , … , Σ n ) \Sigma=(\Sigma_1,\Sigma_2,\dots,\Sigma_n) Σ=(Σ1,Σ2,…,Σn) | 所有高斯分布的协方差矩阵 | |
| θ i = ( μ i , Σ i ) \mathbf{\theta_{i}}=(\mu_i,\;\Sigma_i) θi=(μi,Σi) | i = 1 , 2 , … , n i=1,2,\dots,n i=1,2,…,n | 第 i i i个高斯分布参数的值 |
| θ i ( t ) = ( μ i ( t ) , Σ i ( t ) ) \mathbf{\theta_{i}^{(t)}}=(\mu_i^{(t)},\;\Sigma_i^{(t)}) θi(t)=(μi(t),Σi(t)) | i = 1 , 2 , … , n , t ∈ N + i=1,2,\dots,n,\\t\in\mathbb{N}^+ i=1,2,…,n,t∈N+ | 第 t t t步时，第 i i i个高斯分布参数的值 |
| θ = ( θ 1 , θ 2 , … , θ n ) \mathbf{\theta}=(\mathbf{\theta_1},\mathbf{\theta_2},\dots,\mathbf{\theta_n}) θ=(θ1,θ2,…,θn) | 所有高斯分布参数的值 | |
| θ ( t ) = ( θ 1 ( t ) , θ 2 ( t ) , … , θ n ( t ) ) \mathbf{\theta^{(t)}}=(\mathbf{\theta_1^{(t)}},\mathbf{\theta_2^{(t)}},\dots,\mathbf{\theta_n^{(t)}}) θ(t)=(θ1(t),θ2(t),…,θn(t)) | 第 t t t步时，所有高斯分布参数的值 | |
| 40 | ||
| l ( θ ) l(\mathbf{\theta}) l(θ) | 似然函数的下界函数 | |
| l t ( θ ) l_t(\mathbf{\theta}) lt(θ) | 似然函数的下界函数在 θ ( t ) \mathbf{\theta^{(t)}} θ(t)下的值 |
隐变量（latent variable）
个体到底怎样服从于这
n
n
n个正态分布呢？我们引入隐变量
z
z
z来表示这一点。每一个个体的背后都有一个
z
z
z，
z
z
z的每一个取值都对应着一个高斯分布，
p
(
z
i
=
c
j
)
,
i
=
1
,
2
,
…
,
n
p(z_i=c_j),\;i=1,2,\dots,n
p(zi=cj),i=1,2,…,n即表示第
i
i
i个样本服从第
j
j
j个正态分布的概率，也即第
j
j
j个高斯分布对于第
i
i
i个样本的权重。
| z i z_i zi | c 1 c_1 c1 | c 2 c_2 c2 | ⋯ \cdots ⋯ | c n c_n cn |
|---|---|---|---|---|
| p ( z i = c j ) p(z_i=c_j) p(zi=cj) | π i 1 \pi_{i1} πi1 | π i 2 \pi_{i2} πi2 | ⋯ \cdots ⋯ | π i n \pi_{in} πin |
其中：
∑
j
=
1
n
π
i
j
=
1
,
∀
i
=
1
,
2
,
…
,
N
\sum_{j=1}^n\pi_{ij}=1,\;\forall\;i=1,2,\dots,N
j=1∑nπij=1,∀i=1,2,…,N
算法具体过程
我们使用EM算法来求得参数
θ
\mathbf{\theta}
θ。
似然函数的构建
给定一个数据集，假设样本间相互独立，则可以给出如下似然函数：
log
L
(
θ
)
=
∑
i
=
1
N
log
p
(
x
i
∣
θ
)
=
∑
i
=
1
N
log
(
∑
j
=
1
n
p
(
x
i
,
z
i
=
c
j
∣
θ
)
)
70
\begin{align*} \log L(\mathbf{\theta}) &=\sum_{i=1}^N\log p(x_i|\mathbf{\theta}) \\ &=\sum_{i=1}^N\log \left(\sum_{j=1}^np(x_i,z_i=c_j|\mathbf{\theta})\right) 70 \end{align*}
logL(θ)=i=1∑Nlogp(xi∣θ)=i=1∑Nlog(j=1∑np(xi,zi=cj∣θ))70
step1
定义高斯分布个数
n
n
n，对每个高斯分布设置初始参数值
θ
i
(
0
)
=
(
μ
i
(
0
)
,
Σ
i
(
0
)
)
,
i
=
1
,
2
,
…
,
n
\mathbf{\theta_i^{(0)}}=(\mu_i^{(0)},\;\Sigma_i^{(0)}),\;i=1,2,\dots,n
θi(0)=(μi(0),Σi(0)),i=1,2,…,n，并对所有
π
i
j
\pi_{ij}
πij设置初始参数值
π
i
j
(
0
)
\pi_{ij}^{(0)}
πij(0)。一般通过k-means算法计算初始值，即先使用k-means聚出
n
n
n类，对每一类计算它们的均值和协方差矩阵作为初始值。
step2 E-step
因为对数函数是个凹函数，由Jensen不等式的期望形式（可见(Jensen’s inequality - Wikipedia)），将
p
(
x
i
,
z
i
∣
θ
)
π
i
j
\frac{p(x_i,z_i|\mathbf{\theta})}{\pi_{ij}}
πijp(xi,zi∣θ)看作为
g
(
z
i
)
g(z_i)
g(zi)，
p
(
x
i
,
z
i
=
c
j
∣
θ
)
π
i
j
\frac{p(x_i,z_i=c_j|\mathbf{\theta})}{\pi_{ij}}
πijp(xi,zi=cj∣θ)即为
g
(
z
i
=
c
j
)
g(z_i=c_j)
g(zi=cj)，
π
i
j
\pi_{ij}
πij看作为
p
(
g
(
z
i
=
c
j
)
)
p(g(z_i=c_j))
p(g(zi=cj))，就有：
log
L
(
θ
)
=
∑
i
=
1
N
log
∑
j
=
1
n
p
(
x
i
,
z
i
=
c
j
∣
θ
)
=
∑
i
=
1
N
log
∑
j
=
1
n
π
i
j
p
(
x
i
,
z
i
=
c
j
∣
θ
)
π
i
j
=
∑
i
=
1
N
log
E
[
g
(
z
i
)
]
⩾
∑
i
=
1
N
E
[
log
g
(
z
i
)
]
=
∑
i
=
1
N
∑
j
=
1
n
π
i
j
log
p
(
x
i
,
z
i
=
c
j
∣
θ
)
π
i
j
\begin{align*} \log L(\mathbf{\theta}) &=\sum_{i=1}^N\log \sum_{j=1}^np(x_i,z_i=c_j|\mathbf{\theta}) \\ &=\sum_{i=1}^N\log \sum_{j=1}^n\pi_{ij}\frac{p(x_i,z_i=c_j|\mathbf{\theta})}{\pi_{ij}} \\ &=\sum_{i=1}^N\log E[g(z_i)] \\ &\geqslant\sum_{i=1}^N E[\log g(z_i)] \\ &=\sum_{i=1}^N\sum_{j=1}^n\pi_{ij}\log\frac{p(x_i,z_i=c_j|\mathbf{\theta})}{\pi_{ij}} \end{align*}
logL(θ)=i=1∑Nlogj=1∑np(xi,zi=cj∣θ)=i=1∑Nlogj=1∑nπijπijp(xi,zi=cj∣θ)=i=1∑NlogE[g(zi)]⩾i=1∑NE[logg(zi)]=i=1∑Nj=1∑nπijlogπijp(xi,zi=cj∣θ)
因此我们可以通过提高
∑
i
=
1
N
∑
j
=
1
n
π
i
j
log
p
(
x
i
,
z
i
=
c
j
∣
θ
)
π
i
j
\sum_{i=1}^N\sum_{j=1}^n\pi_{ij}\log\frac{p(x_i,z_i=c_j|\mathbf{\theta})}{\pi_{ij}}
∑i=1N∑j=1nπijlogπijp(xi,zi=cj∣θ) 的下界来提高似然函数值。由Jensen不等式等号成立的条件，当
g
(
z
i
)
=
E
[
g
(
z
i
)
]
g(z_i)=E[g(z_i)]
g(zi)=E[g(zi)]时等号成立。令
E
[
g
(
z
i
)
]
=
C
E[g(z_i)]=C
E[g(zi)]=C，于是此时：
p
(
x
i
,
z
i
=
c
j
∣
θ
)
π
i
j
=
C
,
∀
j
∈
{
1
,
2
,
…
,
n
}
∑
j
=
1
n
p
(
x
i
,
z
i
=
c
j
∣
θ
)
=
∑
j
=
1
n
C
π
i
j
∑
j
=
1
n
p
(
x
i
,
z
i
=
c
j
∣
θ
)
=
C
100
\begin{gather*} \frac{p(x_i,z_i=c_j|\mathbf{\theta})}{\pi_{ij}}=C,\;\forall\;j\in\{1,2,\dots,n\} \\ \sum_{j=1}^np(x_i,z_i=c_j|\mathbf{\theta})=\sum_{j=1}^nC\pi_{ij} \\ \sum_{j=1}^np(x_i,z_i=c_j|\mathbf{\theta})=C \end{gather*} 100
πijp(xi,zi=cj∣θ)=C,∀j∈{1,2,…,n}j=1∑np(xi,zi=cj∣θ)=j=1∑nCπijj=1∑np(xi,zi=cj∣θ)=C100
也即：
π
i
j
=
p
(
x
i
,
z
i
=
c
j
∣
θ
)
∑
k
=
1
n
p
(
x
i
,
z
i
=
c
k
∣
θ
)
=
p
(
x
i
,
z
i
=
c
j
∣
θ
)
p
(
x
i
∣
θ
)
=
p
(
z
i
=
c
j
∣
x
i
,
θ
)
\begin{align*} \pi_{ij} &=\frac{p(x_i,z_i=c_j|\mathbf{\theta})}{\sum_{k=1}^np(x_i,z_i=c_k|\mathbf{\theta})} \\ &=\frac{p(x_i,z_i=c_j|\mathbf{\theta})}{p(x_i|\mathbf{\theta})} \\ &=p(z_i=c_j|x_i,\mathbf{\theta}) \end{align*}
πij=∑k=1np(xi,zi=ck∣θ)p(xi,zi=cj∣θ)=p(xi∣θ)p(xi,zi=cj∣θ)=p(zi=cj∣xi,θ)
这说明在固定高斯分布参数的时候，使似然函数下界达到最大值的
π
i
j
\pi_{ij}
πij的计算公式即为
z
i
=
c
j
z_i=c_j
zi=cj的后验概率。
根据当前参数
θ
(
t
)
=
(
θ
1
(
t
)
,
θ
2
(
t
)
,
…
,
θ
n
(
t
)
)
\mathbf{\theta^{(t)}}=(\mathbf{\theta_1^{(t)}},\mathbf{\theta_2^{(t)}},\dots,\mathbf{\theta_n^{(t)}})
θ(t)=(θ1(t),θ2(t),…,θn(t))和
π
i
j
(
t
)
,
i
=
1
,
2
,
…
,
N
,
j
=
1
,
2
,
…
,
n
\pi_{ij}^{(t)},\;i=1,2,\dots,N,\;j=1,2,\dots,n
πij(t),i=1,2,…,N,j=1,2,…,n，计算每一个隐变量的后验概率分布：
γ
t
(
z
i
=
c
j
)
=
π
i
j
(
t
)
N
(
x
i
∣
μ
j
(
t
)
,
Σ
j
(
t
)
)
∑
k
=
1
n
π
i
k
(
t
)
N
(
x
i
∣
μ
k
(
t
)
,
Σ
k
(
t
)
)
\gamma_t(z_i=c_j)=\frac{\pi_{ij}^{(t)}N(x_i|\mu_j^{(t)},\Sigma_j^{(t)})}{\sum\limits_{k=1}^n\pi_{ik}^{(t)}N(x_i|\mu_k^{(t)},\Sigma_k^{(t)})}
γt(zi=cj)=k=1∑nπik(t)N(xi∣μk(t),Σk(t))πij(t)N(xi∣μj(t),Σj(t))
令
π
i
j
(
t
+
1
)
=
γ
t
(
z
i
=
c
j
)
\pi_{ij}^{(t+1)}=\gamma_t(z_i=c_j)
πij(t+1)=γt(zi=cj)。
step3 M-step
我们在前面使用隐变量的后验概率分布提高了似然函数能达到的理论下界，接下来的工作就是优化参数 θ \mathbf{\theta} θ来提高现在似然函数实际下界的值，那么我们就提高了似然函数的值。也即：
θ
(
t
+
1
)
=
arg max
θ
(
t
+
1
)
l
(
θ
)
=
arg max
θ
(
t
+
1
)
∑
i
=
1
N
∑
j
=
1
n
π
i
j
log
p
(
x
i
,
z
i
=
c
j
∣
θ
)
π
i
j
=
arg max
θ
(
t
+
1
)
∑
i
=
1
N
∑
j
=
1
n
π
i
j
log
N
(
x
i
∣
μ
j
,
Σ
j
)
π
i
j
\begin{align*} \mathbf{\theta^{(t+1)}} &=\argmax_\mathbf{\theta^{(t+1)}}l(\mathbf{\theta}) \\ &=\argmax_\mathbf{\theta^{(t+1)}}\sum_{i=1}^N\sum_{j=1}^n\pi_{ij}\log\frac{p(x_i,z_i=c_j|\mathbf{\theta})}{\pi_{ij}} \\ &=\argmax_\mathbf{\theta^{(t+1)}}\sum_{i=1}^N\sum_{j=1}^n\pi_{ij}\log\frac{N(x_i|\mu_j,\Sigma_j)}{\pi_{ij}} \end{align*}
θ(t+1)=θ(t+1)argmaxl(θ)=θ(t+1)argmaxi=1∑Nj=1∑nπijlogπijp(xi,zi=cj∣θ)=θ(t+1)argmaxi=1∑Nj=1∑nπijlogπijN(xi∣μj,Σj)
对于这个优化问题，我们就可以使用一般的极大似然估计（MLE）去做了。
对
l
(
θ
)
l(\mathbf{\theta})
l(θ)中的各参数求偏导，令偏导为
0
0
0，即可求得
θ
\mathbf{\theta}
θ的更新公式为：
μ
j
(
t
+
1
)
=
∑
i
=
1
N
π
i
j
(
t
+
1
)
x
i
∑
i
=
1
N
π
i
j
(
t
+
1
)
Σ
j
(
t
+
1
)
=
∑
i
=
1
N
π
i
j
(
t
+
1
)
(
x
i
−
μ
j
(
t
)
)
(
x
i
−
μ
j
(
t
)
)
T
∑
i
=
1
N
π
i
j
(
t
+
1
)
\begin{align*} \mu_j^{(t+1)}&=\frac{\sum_{i=1}^N\pi_{ij}^{(t+1)}x_i}{\sum_{i=1}^N\pi_{ij}^{(t+1)}} \\ \Sigma_j^{(t+1)}&=\frac{\sum_{i=1}^N\pi_{ij}^{(t+1)}(x_i-\mu_j^{(t)})(x_i-\mu_j^{(t)})^T}{\sum_{i=1}^N\pi_{ij}^{(t+1)}} \end{align*}
μj(t+1)Σj(t+1)=∑i=1Nπij(t+1)∑i=1Nπij(t+1)xi=∑i=1Nπij(t+1)∑i=1Nπij(t+1)(xi−μj(t))(xi−μj(t))T
step4 重复E-step和M-step直到收敛
由于似然函数值有上界
1
1
1，而EM算法会不断提高似然函数值，所以算法最终会收敛。但是要注意，迭代一定会收敛，但不一定会收敛到真实的参数值，因为可能会陷入局部最优。所以 EM 算法的结果很受初始值的影响。
可设置
ε
\varepsilon
ε，当
∣
l
t
+
1
(
θ
)
−
l
t
(
θ
)
∣
<
ε
|l_{t+1}(\mathbf{\theta})-l_{t}(\mathbf{\theta})|<\varepsilon
∣lt+1(θ)−lt(θ)∣<ε时，认为已经达到收敛，停止训练。
超参数 n n n的选择
使用AIC与BIC去选择合适的
n
n
n。
代码实现
这里使用sklearn中的函数进行建模。
from sklean.mixture import GaussianMixture
GaussianMixture(n_components=1, *, covariance_type='full', tol=0.001, reg_covar=1e-06, max_iter=100, n_init=1, init_params='kmeans', weights_init=None, means_init=None, precisions_init=None, random_state=None, warm_start=False, verbose=0, verbose_interval=10)[source]
具体参数解释与使用请参考sklearn官网文档。GaussianMixture — scikit-learn 1.6.0 documentation
结果
我们最终得到了三组参数，分别为
π
=
(
π
1
,
π
2
,
…
,
π
N
)
,
μ
=
(
μ
1
,
μ
2
,
…
,
μ
n
)
,
Σ
=
(
Σ
1
,
Σ
2
,
…
,
Σ
n
)
\pi=(\pi_1,\pi_2,\dots,\pi_N),\;\mu=(\mu_1,\mu_2,\dots,\mu_n),\;\Sigma=(\Sigma_1,\Sigma_2,\dots,\Sigma_n)
π=(π1,π2,…,πN),μ=(μ1,μ2,…,μn),Σ=(Σ1,Σ2,…,Σn)
-
目的是聚类
对于每一个样本单元产生了一个 π i = ( π i 1 , π i 2 , … , π i n ) \pi_i=(\pi_{i1},\pi_{i2},\dots,\pi_{in}) πi=(πi1,πi2,…,πin)，即样本单元以怎样的概率服从 n n n个具体的高斯分布，取 j = arg max j π i j j=\argmax_j\pi_{ij} j=argmaxjπij，可认为该样本单元属于第 j j j个高斯分布,也即第 j j j类。
GMM优缺点
-
优点：
-
GMM聚出的类可以呈现出椭圆形，优于k-means的圆形
-
GMM得到的最终结果是样本属于每一类的概率，而不像k-means那样必须把样本分给某一类。
-
-
缺点：
-
对大规模数据和多维高斯分布，计算量大，迭代速度慢
-
如果初始值设置不当，收敛过程的计算代价会非常大。
-
EM算法求得的是局部最优解而不一定是全局最优解。
-
References
【机器学习笔记11】高斯混合模型（GMM）【上篇】原理与推导_高斯混合模型推导-CSDN博客
【机器学习】EM——期望最大（非常详细）
高斯混合模型（GMM）