三维点云课程---GMM
由于KMeans缺乏不确定性的指标,GMM算法可以解决这个问题,现在就GMM算法的原理进行解释
1.前言
单变量高斯分布
N
(
x
∣
μ
,
σ
2
)
=
1
(
2
π
σ
2
)
1
/
2
exp
{
−
1
2
σ
2
(
x
−
μ
)
2
}
\mathcal{N}\left(x \mid \mu, \sigma^{2}\right)=\frac{1}{\left(2 \pi \sigma^{2}\right)^{1 / 2}} \exp \left\{-\frac{1}{2 \sigma^{2}}(x-\mu)^{2}\right\}
N(x∣μ,σ2)=(2πσ2)1/21exp{−2σ21(x−μ)2}
多变量的高斯分布,D是维数
N
(
x
∣
μ
,
Σ
)
=
1
(
2
π
)
D
/
2
1
∣
Σ
∣
1
/
2
exp
{
−
1
2
(
x
−
μ
)
T
Σ
−
1
(
x
−
μ
)
}
\mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma})=\frac{1}{(2 \pi)^{D / 2}} \frac{1}{|\boldsymbol{\Sigma}|^{1 / 2}} \exp \left\{-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\mathrm{T}} \boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})\right\}
N(x∣μ,Σ)=(2π)D/21∣Σ∣1/21exp{−21(x−μ)TΣ−1(x−μ)}
2.GMM原理
2.1E步推导
对于一个混合高斯模型,可以写成多个单个高斯模型通过不同的权重进行叠加,其中
π
k
\pi_k
πk表示权重
p
(
x
)
=
∑
k
=
1
K
π
k
N
(
x
∣
μ
k
,
Σ
k
)
p(x)=\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(x \mid \mu_{k}, \Sigma_{k}\right)
p(x)=k=1∑KπkN(x∣μk,Σk)
现在引入一个K维二进制变量z,表达形式如下
z
k
∈
{
0
,
1
}
,
Σ
k
z
k
=
1
z_k\in \{0,1\},\Sigma_{k} z_{k}=1
zk∈{0,1},Σkzk=1
这
p
(
z
k
=
1
)
p(z_k=1)
p(zk=1)是高斯分布
N
(
x
∣
μ
k
,
Σ
k
)
\mathcal{N}\left(x \mid \mu_{k}, \Sigma_{k}\right)
N(x∣μk,Σk)的先验概率,表达如下
p
(
z
k
=
1
)
=
π
k
p(z_k=1)=\pi_k
p(zk=1)=πk
先验的直观理解就是,不告诉我一个数据点在什么位置,只知道高斯模型的概率。以上式子有
0
≤
π
k
≤
1
,
∑
k
=
1
K
π
k
=
1
0 \leq \pi_{k} \leq 1, \sum_{k=1}^{K} \pi_{k}=1
0≤πk≤1,∑k=1Kπk=1的约束,联合起来,
p
(
z
)
p(z)
p(z)的表达如下
p
(
z
)
=
∏
k
=
1
K
π
k
z
k
p(z)=\prod_{k=1}^{K} \pi_{k}^{z_{k}}
p(z)=k=1∏Kπkzk
其中
z
=
[
z
1
,
.
.
.
,
z
k
,
.
.
.
,
z
K
]
z=[z_1,...,z_k,...,z_K]
z=[z1,...,zk,...,zK]。那么
p
(
x
∣
z
k
=
1
)
p(x|z_k=1)
p(x∣zk=1)就坍塌成一个高斯分布了。
p
(
x
∣
z
k
=
1
)
=
N
(
x
∣
μ
k
,
Σ
k
)
p\left(x \mid z_{k}=1\right)=\mathcal{N}\left(x \mid \mu_{k}, \Sigma_{k}\right)
p(x∣zk=1)=N(x∣μk,Σk)
所以
p
(
x
∣
z
)
p(x|z)
p(x∣z)表示如下
p
(
x
∣
z
)
=
∏
k
=
1
K
N
(
x
∣
μ
k
,
Σ
k
)
z
k
p(x \mid z)=\prod_{k=1}^{K} \mathcal{N}\left(x \mid \mu_{k}, \Sigma_{k}\right)^{z_{k}}
p(x∣z)=k=1∏KN(x∣μk,Σk)zk
根据概率中的联合密度和边缘概率的关系
p
(
x
)
=
∑
z
p
(
x
,
z
)
p(x)=\sum_{z} p(x, z)
p(x)=∑zp(x,z),可以得到
p
(
x
)
=
∑
z
p
(
x
,
z
)
=
∑
z
p
(
z
)
p
(
x
∣
z
)
=
∑
z
∏
k
=
1
K
π
k
z
k
∏
k
=
1
K
N
(
x
∣
μ
k
,
Σ
k
)
z
k
=
∑
k
=
1
K
π
k
N
(
x
∣
μ
k
,
Σ
k
)
p(x)=\sum_{z} p(x, z)=\sum_{z}p(z)p(x|z)=\sum_{z}\prod_{k=1}^{K} \pi_{k}^{z_{k}}\prod_{k=1}^{K} \mathcal{N}\left(x \mid \mu_{k}, \Sigma_{k}\right)^{z_{k}}=\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(x \mid \mu_{k}, \Sigma_{k}\right)
p(x)=z∑p(x,z)=z∑p(z)p(x∣z)=z∑k=1∏Kπkzkk=1∏KN(x∣μk,Σk)zk=k=1∑KπkN(x∣μk,Σk)
根据贝叶斯概率分布公式可以得到
p
(
z
∣
x
)
=
p
(
z
,
x
)
p
(
x
)
=
p
(
z
)
p
(
x
∣
z
)
∑
j
=
1
K
p
(
x
∣
z
j
)
p
(
z
j
)
p(\mathbf{z} \mid \mathbf{x})=\frac{p(\mathbf{z},\mathbf{x})}{p(\mathbf{x})}=\frac{p(\mathbf{z}) p(\mathbf{x} \mid \mathbf{z})}{\sum_{j=1}^{K} p\left(\mathbf{x} \mid z_{j}\right) p\left(z_{j}\right)}
p(z∣x)=p(x)p(z,x)=∑j=1Kp(x∣zj)p(zj)p(z)p(x∣z)
令
γ
(
z
k
)
=
p
(
z
k
=
1
)
\gamma(z_k)=p(z_k=1)
γ(zk)=p(zk=1),那么
γ
(
z
k
)
≡
p
(
z
k
=
1
∣
x
)
=
p
(
z
k
=
1
)
p
(
x
∣
z
k
=
1
)
∑
j
=
1
K
p
(
z
j
=
1
)
p
(
x
∣
z
j
=
1
)
=
π
k
N
(
x
n
∣
μ
k
,
Σ
k
)
∑
j
=
1
K
π
j
N
(
x
n
∣
μ
j
,
Σ
j
)
\begin{aligned} \gamma\left(z_{k}\right) \equiv p\left(z_{k}=1 \mid \mathbf{x}\right) &=\frac{p\left(z_{k}=1\right) p\left(\mathbf{x} \mid z_{k}=1\right)}{\sum_{j=1}^{K} p\left(z_{j}=1\right) p\left(\mathbf{x} \mid z_{j}=1\right)} \\ &=\frac{\pi_{k} \mathcal{N}\left(\mathbf{x}_{n} \mid \mu_{k}, \Sigma_{k}\right)}{\sum_{j=1}^{K} \pi_{j} \mathcal{N}\left(\mathbf{x}_{n} \mid \mu_{j}, \Sigma_{j}\right)} \end{aligned}
γ(zk)≡p(zk=1∣x)=∑j=1Kp(zj=1)p(x∣zj=1)p(zk=1)p(x∣zk=1)=∑j=1KπjN(xn∣μj,Σj)πkN(xn∣μk,Σk)
物理意义
γ n k = p ( z n k = 1 ∣ x n ) \gamma_{nk}=p(z_{nk}=1 | x_n) γnk=p(znk=1∣xn):一个点属于哪个类的可能性
2.2.M步推导
上式当我们知道
{
π
k
,
μ
k
,
Σ
k
}
\{\pi_k,\mu_k,\Sigma_k\}
{πk,μk,Σk}和
x
x
x,我们可以得到
γ
(
z
k
)
\gamma(z_k)
γ(zk)。但是实际是我只知道数据点
x
1
,
.
.
.
,
x
N
{x_1,...,x_N}
x1,...,xN和聚类的个数,让我们来估计
{
π
k
,
μ
k
,
Σ
k
}
,
γ
(
z
k
)
\{\pi_k,\mu_k,\Sigma_k\},\gamma(z_k)
{πk,μk,Σk},γ(zk),此时我们需要最大似然估计
p
(
x
)
=
∑
k
=
1
K
π
k
N
(
x
∣
μ
k
,
Σ
k
)
ln
p
(
X
∣
π
,
μ
,
Σ
)
=
∑
n
=
1
N
ln
{
∑
k
=
1
K
π
k
N
(
x
n
∣
μ
k
,
Σ
k
)
}
\begin{aligned} p(\mathbf{x}) &=\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\mathbf{x} \mid \mu_{k}, \Sigma_{k}\right) \\ \ln p(\mathbf{X} \mid \pi, \mu, \Sigma) &=\sum_{n=1}^{N} \ln \left\{\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\mathbf{x}_{\mathbf{n}} \mid \mu_{k}, \Sigma_{k}\right)\right\} \end{aligned}
p(x)lnp(X∣π,μ,Σ)=k=1∑KπkN(x∣μk,Σk)=n=1∑Nln{k=1∑KπkN(xn∣μk,Σk)}
其中
X
=
[
x
1
,
.
.
.
,
x
N
]
X=[x_1,...,x_N]
X=[x1,...,xN],因此最大似然表达如下
π
,
μ
,
Σ
=
arg
max
π
,
μ
,
Σ
ln
p
(
X
∣
π
,
μ
,
Σ
)
=
arg
max
π
,
μ
,
Σ
∑
n
=
1
N
ln
{
∑
k
=
1
K
π
k
N
(
x
n
∣
μ
k
,
Σ
k
)
}
\pi, \mu, \boldsymbol{\Sigma}=\underset{\pi, \mu, \boldsymbol{\Sigma}}{\arg \max } \ln p(\mathbf{X} \mid \pi, \mu, \boldsymbol{\Sigma})=\underset{\pi, \mu, \boldsymbol{\Sigma}}{\arg \max } \sum_{n=1}^{N} \ln \left\{\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\mathbf{x}_{n} \mid \mu_{k}, \boldsymbol{\Sigma}_{k}\right)\right\}
π,μ,Σ=π,μ,Σargmaxlnp(X∣π,μ,Σ)=π,μ,Σargmaxn=1∑Nln{k=1∑KπkN(xn∣μk,Σk)}
对于求解多个变量的最优问题,大多数是固定其中多个变量,只求解一个,进行迭代求解即可。
2.2.1固定 π k , Σ k \pi_k,\Sigma_k πk,Σk,求解 μ k \mu_k μk
对上式进行求导可知,省略求导过程
0
=
−
∑
n
=
1
N
π
k
N
(
x
n
∣
μ
k
,
Σ
k
)
∑
j
π
j
N
(
x
n
∣
μ
j
,
Σ
j
)
Σ
K
(
x
n
−
μ
k
)
⇔
0
=
−
∑
n
=
1
N
γ
(
z
n
k
)
Σ
K
(
x
n
−
μ
k
)
0 = - \sum\limits_{n = 1}^N {\frac{{{\pi _k}{\mathcal N}\left( {{x_n}\mid {\mu _k},{\Sigma _k}} \right)}}{{\sum\nolimits_j {{\pi _j}{\mathcal N}\left( {{x_n}\mid {\mu _j},{\Sigma _j}} \right)} }}} {\Sigma _K}({x_n} - {\mu _k})\\ \Leftrightarrow 0 = - \sum\limits_{n = 1}^N {\gamma ({z_{nk}})} {\Sigma _K}({x_n} - {\mu _k})
0=−n=1∑N∑jπjN(xn∣μj,Σj)πkN(xn∣μk,Σk)ΣK(xn−μk)⇔0=−n=1∑Nγ(znk)ΣK(xn−μk)
化简得
μ
k
=
∑
n
=
1
N
γ
(
z
n
k
)
x
n
N
k
,
N
k
=
∑
n
=
1
N
γ
(
z
n
k
)
{\mu _k} = \frac{{\sum\limits_{n = 1}^N {\gamma ({z_{nk}})} {x_n}}}{{{N_k}}},{N_k} = \sum\limits_{n = 1}^N {\gamma ({z_{nk}})}
μk=Nkn=1∑Nγ(znk)xn,Nk=n=1∑Nγ(znk)
物理意义
N
k
N_k
Nk:分配给聚类K的有效点数,即有多个点属于k这个类
μ k \mu_k μk:所有点在k这个类上的加权平均,这个不同于KMeans简单的进行求平均,在这里引入了权重的概念
2.2.2固定 π k , μ k \pi_k,\mu_k πk,μk,求解 Σ k \Sigma_k Σk
求导过程省略
Σ
k
=
1
N
k
∑
n
=
1
N
γ
(
z
n
k
)
(
x
n
−
μ
k
)
(
x
n
−
μ
k
)
T
,
N
k
=
∑
n
=
1
N
γ
(
z
n
k
)
{\Sigma _k} = \frac{1}{{{N_k}}}\sum\limits_{n = 1}^N {\gamma ({z_{nk}})({x_n} - {\mu _k})} {({x_n} - {\mu _k})^T},{N_k} = \sum\limits_{n = 1}^N {\gamma ({z_{nk}})}
Σk=Nk1n=1∑Nγ(znk)(xn−μk)(xn−μk)T,Nk=n=1∑Nγ(znk)
物理意义
Σ k \Sigma_k Σk:是以 μ k \mu_k μk为中心的所有点方差的加权平均值,权重为后验概率 γ ( z n k ) \gamma(z_{nk}) γ(znk)
2.2.3固定 μ k , Σ k \mu_k,\Sigma_k μk,Σk,求解 π k \pi_k πk
解这个问题时,存在一个限制条件$\sum\nolimits_{k = 1}^K {{\pi _k} = 1} $。关于有限制条件的求导,使用拉格朗日求导法
ln
p
(
X
∣
π
,
μ
,
Σ
)
+
λ
(
∑
k
=
1
K
π
k
−
1
)
\ln p(\mathbf{X} \mid \pi, \mu, \Sigma) +\lambda(\sum \limits_{k=1}^K \pi_k-1)
lnp(X∣π,μ,Σ)+λ(k=1∑Kπk−1)
对上式进行关于
μ
k
\mu_k
μk的求导
∑
n
=
1
N
N
(
x
n
∣
μ
k
,
Σ
k
)
∑
j
π
j
N
(
x
n
∣
μ
j
,
Σ
j
)
+
λ
=
0
\sum\limits_{n = 1}^N {\frac{{{\mathcal N}\left( {{x_n}\mid {\mu _k},{\Sigma _k}} \right)}}{{\sum\nolimits_j {{\pi _j}{\mathcal N}\left( {{x_n}\mid {\mu _j},{\Sigma _j}} \right)} }}}+\lambda=0
n=1∑N∑jπjN(xn∣μj,Σj)N(xn∣μk,Σk)+λ=0
首先求解
λ
\lambda
λ,等式两边同时乘以
∑
k
=
1
K
π
k
\sum \limits_{k=1}^K \pi_k
k=1∑Kπk,化简得
∑
k
=
1
K
π
k
∑
n
=
1
N
N
(
x
n
∣
μ
k
,
Σ
k
)
∑
j
π
j
N
(
x
n
∣
μ
j
,
Σ
j
)
+
∑
k
=
1
K
π
k
λ
=
0
\sum \limits_{k=1}^K \pi_k {\sum\limits_{n = 1}^N {\frac{{{\mathcal N}\left( {{x_n}\mid {\mu _k},{\Sigma _k}} \right)}}{{\sum\nolimits_j {{\pi _j}{\mathcal N}\left( {{x_n}\mid {\mu _j},{\Sigma _j}} \right)} }}} }+\sum \limits_{k=1}^K \pi_k\lambda=0
k=1∑Kπkn=1∑N∑jπjN(xn∣μj,Σj)N(xn∣μk,Σk)+k=1∑Kπkλ=0
∑
k
=
1
K
∑
n
=
1
N
π
k
N
(
x
n
∣
μ
k
,
Σ
k
)
∑
j
π
j
N
(
x
n
∣
μ
j
,
Σ
j
)
+
∑
k
=
1
K
π
k
λ
=
0
\sum \limits_{k=1}^K {\frac{{\sum\limits_{n = 1}^N {\pi_k}{\mathcal N}\left( {{x_n}\mid {\mu _k},{\Sigma _k}} \right)}}{{\sum\nolimits_j {{\pi _j}{\mathcal N}\left( {{x_n}\mid {\mu _j},{\Sigma _j}} \right)} }}}+\sum \limits_{k=1}^K \pi_k\lambda =0
k=1∑K∑jπjN(xn∣μj,Σj)n=1∑NπkN(xn∣μk,Σk)+k=1∑Kπkλ=0
(
∑
n
=
1
N
1
)
+
λ
=
0
λ
=
−
N
(\sum \limits_{n=1}^N 1)+\lambda=0\\ \lambda=-N
(n=1∑N1)+λ=0λ=−N
再次将
λ
\lambda
λ代入,进而求解
π
k
\pi_k
πk
∑
n
=
1
N
N
(
x
n
∣
μ
k
,
Σ
k
)
∑
j
π
j
N
(
x
n
∣
μ
j
,
Σ
j
)
+
λ
=
0
1
π
k
∑
n
=
1
N
π
k
N
(
x
n
∣
μ
k
,
Σ
k
)
∑
j
π
j
N
(
x
n
∣
μ
j
,
Σ
j
)
−
N
=
0
1
π
k
∑
n
=
1
N
γ
n
k
−
N
=
0
π
k
=
N
k
N
\sum\limits_{n = 1}^N {\frac{{{\mathcal N}\left( {{x_n}\mid {\mu _k},{\Sigma _k}} \right)}}{{\sum\nolimits_j {{\pi _j}{\mathcal N}\left( {{x_n}\mid {\mu _j},{\Sigma _j}} \right)} }}}+\lambda=0 \\ \frac{1}{\pi_k} \sum\limits_{n = 1}^N {\frac{{{\pi_k}{\mathcal N}\left( {{x_n}\mid {\mu _k},{\Sigma _k}} \right)}}{{\sum\nolimits_j {{\pi _j}{\mathcal N}\left( {{x_n}\mid {\mu _j},{\Sigma _j}} \right)} }}}-N=0\\ \frac {1}{\pi_k} \sum \limits_{n=1}^N\gamma_{nk}-N=0\\ \pi_k=\frac {N_k}{N}
n=1∑N∑jπjN(xn∣μj,Σj)N(xn∣μk,Σk)+λ=0πk1n=1∑N∑jπjN(xn∣μj,Σj)πkN(xn∣μk,Σk)−N=0πk1n=1∑Nγnk−N=0πk=NNk
物理意义
π k \pi_k πk:表示的是实际上有多少个点属于我这个类$\div $所有点的数量
3.GMM总结
-
初始化每一个高斯模型的权重 π k \pi_k πk,均值 μ k \mu_k μk和方差 Σ k \Sigma_k Σk
# π(alpha),初始化权重值,即初始时,平分权重 alpha = np.ones(self.n_clusters) / self.n_clusters # μ(mu),随机初始化期望值 mu = np.array([[.403, .237], [.714, .346], [.532, .472]]) # ∑(sigma),这是一个k*data.shape[1]*data.shape[1]三维矩阵,目的是一次性的将方差储存起来,其中初值为对角矩阵, # 这里以3*2*2,∑为[[[0.1,0],[0,0.1]],[[0.1,0],[0,0.1]],[[0.1,0],[0,0.1]]],第一块(sigma[0])为k=0对应的方差,同理第二块,第三块也是这样 sigma = np.full((self.n_clusters, n_features, n_features), np.diag(np.full(n_features, 0.1)))
其中 self.n_clusters表示为聚类的个数,n_features表示数据的列数,这里是2,数据是以N*2表示的
-
E步:通过权重,均值,方差估计后验概率 γ n k = p ( z n k = 1 ∣ x n ) \gamma_{nk}=p(z_{nk}=1 | x_n) γnk=p(znk=1∣xn),表示一个点属于哪个类的可能性
gamma = self.computeGamma(data, mu, sigma, alpha, self.multiGaussian) #维数N*k
计算 γ ( z n k ) \gamma(z_{nk}) γ(znk),调用computeGamma函数和multiGaussian函数,函数如下
def computeGamma(self, X, mu, sigma, alpha, multiGaussian): n_samples = X.shape[0] n_clusters = len(alpha) gamma = np.zeros((n_samples, n_clusters)) p = np.zeros(n_clusters) g = np.zeros(n_clusters) for i in range(n_samples): for j in range(n_clusters): p[j] = multiGaussian(X[i], mu[j], sigma[j]) #对应公式中的 N(xn|μj,∑j) g[j] = alpha[j] * p[j] #对应公式中的πj*N(xn|μj,∑j) for k in range(n_clusters): if np.sum(g)!=0: gamma[i, k] = g[k] / np.sum(g) #对应公式的γnk的求解 return gamma
其中 p [ j ] p[j] p[j]对应 N ( x ∣ μ k , Σ k ) \mathcal{N}\left(x \mid \mu_{k}, \Sigma_{k}\right) N(x∣μk,Σk), g [ j ] g[j] g[j]对应 π k N ( x ∣ μ k , Σ k ) \pi_{k} \mathcal{N}\left(x \mid \mu_{k}, \Sigma_{k}\right) πkN(x∣μk,Σk)
def multiGaussian(self, x, mu, sigma): return 1 / ((2 * np.pi) * pow(np.linalg.det(sigma), 0.5)) * np.exp(-0.5 * (x - mu).dot(np.linalg.pinv(sigma)).dot((x - mu).T))
对应于前言的多元高斯分布函数
-
M步:计算最新的权重,均值,和方差
{ μ k = ∑ n = 1 N γ ( z n k ) x n N k Σ k = 1 N k ∑ n = 1 N γ ( z n k ) ( x n − μ k ) ( x n − μ k ) T π k = N k N N k = ∑ n = 1 N γ ( z n k ) N = ∑ n = 1 N 1 \left\{ \begin{array}{l} {\mu _k} = \frac{{\sum\limits_{n = 1}^N {\gamma ({z_{nk}})} {x_n}}}{{{N_k}}}\\ {\Sigma _k} = \frac{1}{{{N_k}}}\sum\limits_{n = 1}^N {\gamma ({z_{nk}})({x_n} - {\mu _k})} {({x_n} - {\mu _k})^T}\\ {\pi _k} = \frac{{{N_k}}}{N}\\ {N_k} = \sum\limits_{n = 1}^N {\gamma ({z_{nk}})} \\ N = \sum\limits_{n = 1}^N 1 \end{array} \right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧μk=Nkn=1∑Nγ(znk)xnΣk=Nk1n=1∑Nγ(znk)(xn−μk)(xn−μk)Tπk=NNkNk=n=1∑Nγ(znk)N=n=1∑N1
计算权重alpha = np.sum(gamma, axis=0) / n_samples
对应公式 π k = N k N {\pi _k} = \frac{{{N_k}}}{N} πk=NNk
计算均值
#更新均值 μ(mu) gamma_xn=data * gamma[:, i].reshape((n_samples, 1)) #维数N*k mu_molecule=np.sum(gamma_xn, axis=0) #维数(k,) mu_denominator=np.sum(gamma, axis=0)[i] mu[i]=mu_molecule/mu_denominator
其中gamma_xn表示 γ ( z n k ) x n \gamma ({z_{nk}}){x_n} γ(znk)xn,mu_molecule表示 ∑ n = 1 N γ ( z n k ) x n \sum \limits_{n=1}^N \gamma(z_{nk})x_n n=1∑Nγ(znk)xn,mu_denominator表示 N k N_k Nk
计算协方差矩阵
sigma[i] = 0 for j in range(n_samples): #更新方差∑(sigma) xn_muk=data[j].reshape((1, n_features))-mu[i] #维数 1*k sigma[i]+=(xn_muk).T.dot(xn_muk)*gamma[j,i] #维数 k*k sigma[i] = sigma[i] / (np.sum(gamma, axis=0)[i])
其中xn_muk表示为 x n − μ k x_n-\mu_k xn−μk。注意这里和公式不一样,是 ( x n − μ k ) T ( x n − μ k ) (x_n-\mu_k)^T(x_n-\mu_k) (xn−μk)T(xn−μk)这样才可以表示为协方差矩阵,维数为2*2。
-
评估log函数的,如果收敛,则停止,否则返回到E步
效果如下,不同于KMeans的是,数据点有个渐变的过程,并不是非0即1
4.GMM完整代码
import numpy as np
from numpy import *
import pylab
import random, math
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from scipy.stats import multivariate_normal
plt.style.use('seaborn')
class GMM(object):
def __init__(self, n_clusters, max_iter=50):
self.n_clusters = n_clusters #聚类个数
self.ITER = max_iter #最大迭代次数
self.mu = 0 #初始化均值
self.sigma = 0 #初始化协方差矩阵
self.alpha = 0 #初始化权重
#计算多元高斯分布函数,注意这里是(x-mu)*sigma^(-1)*(x-mu)^T公式没有错,因为这里算出来的是一个数,并且不能直接运用在三维中,缺了个系数
def multiGaussian(self, x, mu, sigma):
return 1 / ((2 * np.pi) * pow(np.linalg.det(sigma), 0.5)) * np.exp(-0.5 * (x - mu).dot(np.linalg.pinv(sigma)).dot((x - mu).T))
#计算后验概率γnk,维数为data.shape[0]*k,这里以k=2为例,第一列放的是所有属于k=0类的概率,第二列放的是所有属于k=1类的概率
def computeGamma(self, X, mu, sigma, alpha, multiGaussian):
n_samples = X.shape[0]
n_clusters = len(alpha)
gamma = np.zeros((n_samples, n_clusters))
p = np.zeros(n_clusters)
g = np.zeros(n_clusters)
for i in range(n_samples):
for j in range(n_clusters):
p[j] = multiGaussian(X[i], mu[j], sigma[j]) #对应公式中的 N(xn|μj,∑j)
g[j] = alpha[j] * p[j] #对应公式中的πj*N(xn|μj,∑j)
for k in range(n_clusters):
if np.sum(g)!=0:
gamma[i, k] = g[k] / np.sum(g) #对应公式的γnk的求解
return gamma
#GMM的核心代码
def fit(self, data):
n_samples = data.shape[0]
n_features = data.shape[1]
'''
mu=data[np.random.choice(range(n_samples),self.n_clusters)]
'''
#初始化初值
# π(alpha),初始化权重值,即初始时,平分权重
alpha = np.ones(self.n_clusters) / self.n_clusters
# μ(mu),随机初始化期望值
mu = np.array([[.403, .237], [.714, .346], [.532, .472]])
# ∑(sigma),这是一个k*data.shape[1]*data.shape[1]三维矩阵,目的是一次性的将方差储存起来,其中初值为对角矩阵,
# 这里以3*2*2,∑为[[[0.1,0],[0,0.1]],[[0.1,0],[0,0.1]],[[0.1,0],[0,0.1]]],第一块(sigma[0])为k=0对应的方差,同理第二块,第三块也是这样
sigma = np.full((self.n_clusters, n_features, n_features), np.diag(np.full(n_features, 0.1)))
for i in range(self.ITER):
#E步,更新后验概率γnk,维数N*k
gamma = self.computeGamma(data, mu, sigma, alpha, self.multiGaussian)
#M步,更新三大参数:权重,均值,方差
#计算权重π(alpha),维数1*k
alpha = np.sum(gamma, axis=0) / n_samples
for i in range(self.n_clusters):
#计算均值 μ(mu)
gamma_xn=data * gamma[:, i].reshape((n_samples, 1)) #维数N*k
mu_molecule=np.sum(gamma_xn, axis=0) #维数(k,)
mu_denominator=np.sum(gamma, axis=0)[i]
mu[i]=mu_molecule/mu_denominator
sigma[i] = 0
for j in range(n_samples):
#计算协方差矩阵∑(sigma)
xn_muk=data[j].reshape((1, n_features)) #维数 1*k
#维数 k*k
sigma[i]+=(xn_muk-mu[i]).T.dot(xn_muk-mu[i])*gamma[j,i]
sigma[i] = sigma[i] / (np.sum(gamma, axis=0)[i])
self.mu = mu #更新权重π(alpha)
self.sigma = sigma #更新均值μ(mu)
self.alpha = alpha #更新协方差矩阵∑(sigma)
#GMM算法中的最终预测点的分类函数
def predict(self, data):
pred = self.computeGamma(data, self.mu, self.sigma, self.alpha, self.multiGaussian)
cluster_results = np.argmax(pred, axis=1)
return cluster_results
##初始化画布为4行2列
fig, ax = plt.subplots(4, 2)
## 可视化函数
## 输入:x 原始数的路径
# i 画布的行
# j 画布的列
# n_clusters 聚类的个数,认为指定
def show_result(x,i,j,n_clusters):
ax[i][j].scatter(x[:, 0], x[:, 1], s=20, c="b", marker='o')
ax[i][j].set_xlabel('x')
ax[i][j].set_ylabel('y')
gmm = GMM(n_clusters)
gmm.fit(x)
cat = gmm.predict(x)
list_max = max(cat)
if list_max == 1:
for idx, value in enumerate(cat):
if value == 0:
ax[i][j+1].scatter(x[idx, 0], x[idx, 1], s=20, c="r", marker='*')
elif value == 1:
ax[i][j+1].scatter(x[idx, 0], x[idx, 1], s=20, c="g", marker='+')
elif list_max == 2:
for idx, value in enumerate(cat):
if value == 0:
ax[i][j+1].scatter(x[idx, 0], x[idx, 1], s=20, c="r", marker='*')
elif value == 1:
ax[i][j+1].scatter(x[idx, 0], x[idx, 1], s=20, c="g", marker='+')
elif value == 2:
ax[i][j+1].scatter(x[idx, 0], x[idx, 1], s=20, c="y", marker='^')
ax[i][j+1].set_xlabel('x')
ax[i][j+1].set_ylabel('y')
## 读取数据的函数
def read_txt(path):
filename = path # txt文件和当前脚本在同一目录下,所以不用写具体路径
pos = []
Efield = []
with open(filename, 'r') as file_to_read:
while True:
lines = file_to_read.readline() # 整行读取数据
if not lines:
break
pass
p_tmp, E_tmp = [float(i) for i in lines.split(",")] # 将整行数据分割处理,如果分割符是空格,括号里就不用传入参数,如果是逗号, 则传入‘,'字符。
pos.append(p_tmp) # 添加新读取的数据
Efield.append(E_tmp)
pass
pos = np.array(pos) # 将数据从list类型转换为array类型。
Efield = np.array(Efield)
x=np.vstack(pos)
y=np.vstack(Efield)
X=np.concatenate((x,y),axis=1)
return X
if __name__ == '__main__':
x_circle = read_txt("circle.txt")
show_result(x_circle, 0, 0, 2)
x_moons = read_txt("moons.txt")
show_result(x_moons, 1, 0, 2)
x_blobs = read_txt("varied.txt")
show_result(x_blobs, 2, 0, 3)
x_aniso = read_txt("aniso.txt")
show_result(x_aniso, 3, 0, 3)
plt.show()
仿真结果如下
数据集下载参见之前的KMeans文章
5.GMM的一些缺点
最大似然公式存在奇异性的问题:假设GMM的协方差矩阵为
Σ
k
=
σ
k
2
I
\Sigma_k=\sigma^2_{k}I
Σk=σk2I,
I
I
I为单位向量,有一个数据点
x
n
=
μ
j
x_n=\mu_j
xn=μj,那么高斯分布为
N
(
x
∣
x
n
,
σ
j
2
I
)
=
1
(
2
π
)
1
/
2
1
σ
j
\mathcal{N}\left(x \mid x_n, \sigma_{j}^2I\right)=\frac{1}{(2\pi)^{1/2}} \frac{1}{\sigma_j}
N(x∣xn,σj2I)=(2π)1/21σj1
那么如果
σ
j
\sigma_j
σj非常小,高斯分布函数就非常大,系统可能会崩溃,例如如果一个点就是一个高斯分布,那么
σ
=
0
\sigma=0
σ=0,N->0,
∏
p
(
x
i
)
\prod p(x_i)
∏p(xi)->0。因此,在GMM最大似然公式中,我们将高斯点放在单个点上,这会导致很大的问题。
解决方法:
1.采用MAP的方法,初始考虑 μ k , Σ k \mu_k,\Sigma_k μk,Σk,假如高斯分布的方差不能太小,就不会掉到奇点问题,但比较复杂
2.如果出现某个高斯点的方差很小,随机初始化另外的一个值。
其次GMM存在和KMeans同样的通病,需要人工给定聚类数。