本文主要参考周志华《机器学习》的9.4.3章节,对高斯混合聚类的原理做简单介绍,并使用numpy实现GMM。
要想很好得理解掌握高斯混合聚类算法,以我的学习经验来看,需要掌握两方面背景知识。
- 多维正态分布
- EM算法
关于上述两方面知识,我只做简单的介绍。
多维正态分布
首先,什么是多维正态分布?就是多变量的正态分布。我们所熟知的正态分布往往是一维的,但在现实中,我们所获得的数据往往是多维的。这就需要用到多维正态分布。因为二维的正态分布是两个独立同分布的联合概率分布,所以从一维到2维的推广:
p
(
x
,
y
)
=
p
(
x
)
p
(
y
)
=
1
2
π
exp
(
−
x
2
+
y
2
2
)
p(x, y)=p(x) p(y)=\frac{1}{2 \pi} \exp \left(-\frac{x^{2}+y^{2}}{2}\right)
p(x,y)=p(x)p(y)=2π1exp(−2x2+y2)
用一维向量表示变量为:
V
=
[
x
,
y
]
T
V = [x,y]^T
V=[x,y]T
则可写成:
p
(
x
,
y
)
=
1
2
π
exp
(
−
V
V
T
2
)
p(x,y) = \frac{1}{2 \pi} \exp \left(-\frac{VV^T}{2}\right)
p(x,y)=2π1exp(−2VVT)
进一步对V做线性变换,令V=A(V-μ),则:
p
(
x
)
=
∣
d
e
t
(
A
)
∣
2
π
exp
[
−
1
2
(
x
−
μ
)
T
A
T
A
(
x
−
μ
)
]
p(\mathbf{x})=\frac{|det{(A)}|}{2 \pi} \exp \left[-\frac{1}{2}(\mathbf{x}-\mu)^{T} A^{T} A(\mathbf{x}-\mu)\right]
p(x)=2π∣det(A)∣exp[−21(x−μ)TATA(x−μ)]
观察上面的公式发现,
Σ
=
(
A
T
A
)
−
1
是
多
维
高
斯
分
布
的
协
方
差
Σ={(A^TA)}^{-1}是多维高斯分布的协方差
Σ=(ATA)−1是多维高斯分布的协方差,那么上式可以进一步写成:
p
(
x
)
=
1
2
π
∣
Σ
∣
1
/
2
exp
[
−
1
2
(
x
−
μ
)
T
Σ
−
1
(
x
−
μ
)
]
p(\mathbf{x})=\frac{1}{2 \pi|Σ|^{1/2}} \exp \left[-\frac{1}{2}(\mathbf{x}-\mu)^{T}Σ^{-1}(\mathbf{x}-\mu)\right]
p(x)=2π∣Σ∣1/21exp[−21(x−μ)TΣ−1(x−μ)]
以上便是多维高斯分布的简单介绍。
EM算法
EM算法主要是针对模型中具有隐变量时的参数估计。例如估计
λ
=
a
r
g
m
a
x
(
L
(
x
∣
y
,
λ
)
)
λ = argmax(L(x|y,λ))
λ=argmax(L(x∣y,λ))
其中y即为隐变量,这种似然函数很难估计,所以采用两步的方法。
主要分为E步和M步。
E步:求关于隐变量的似然函数的期望。一般我们成为Q函数。
M步:通过拉格朗日乘子法最大化Q函数。
高斯混合模型
该模型认为,需要聚类的样本是符合k个多维高斯分布的函数的组合。用公式表示如下:
p
M
(
x
)
=
∑
i
=
1
k
α
i
⋅
p
(
x
∣
μ
i
,
Σ
i
)
p_{\mathcal{M}}(\boldsymbol{x})=\sum_{i=1}^{k} \alpha_{i} \cdot p\left(\boldsymbol{x} | \boldsymbol{\mu}_{i}, \mathbf{\Sigma}_{i}\right)
pM(x)=i=1∑kαi⋅p(x∣μi,Σi)
其中α为各个分布的权重,和为1。后面的部分是一个高斯分布。因为该模型就存在隐变量,所以求解使用EM算法。算法的具体推导详见《统计学习方法》李航。
求解该模型的算法如下:
numpy实现的代码如下:
import numpy as np
from sklearn.datasets import load_iris
import clustering as cl
def multi_norm(x, mu, sigma):
"""
return the probability of Multidimensional gaussian distribution, and there is a better implementation in scipy.stats
返回多维高斯分布的结果,该方法在scipi库中有更好的实现。
:param x: x
:param mu: mean vector
:param sigma: covariance matrix
:return: the probability of Multidimensional gaussian distribution
"""
det = np.linalg.det(sigma)
inv = np.matrix(np.linalg.inv(sigma))
x_mu = np.matrix(x - mu).T
const = 1 / (((2 * np.pi) ** (len(x) / 2)) * (det ** (1 / 2)))
exp = -0.5 * x_mu.T * inv * x_mu
return float(const * np.exp(exp))
def distance(a, b, p):
a = np.array(a)
b = np.array(b)
return np.sum((a - b) ** p) ** (1 / p)
def pairing(data, truth, label):
datatemp = data.copy()
centerTruth = []
center = []
new_label = np.zeros(label.shape) - 1
for i in range(0, np.max(truth) + 1):
centerTruth.append(list(np.mean(datatemp[np.argwhere(truth == i)], axis=0)))
center.append(list(np.mean(datatemp[np.argwhere(label == i)], axis=0)))
for i in range(0, np.max(truth) + 1):
temp = []
for j in range(0, np.max(truth) + 1):
temp.append(distance(centerTruth[i], center[j], 2))
number = temp.index(min(temp))
print(number)
new_label[label == number] = i
return new_label
def calcu_acc(truth, label):
temp = np.zeros(label.shape)
temp[np.argwhere(label == truth)] = 1
return np.sum(temp) / label.shape[0]
class GMM(object):
def __init__(self, n_clusters, max_iter, init_params="random"):
"""
初始化变量:聚类数目以及最大迭代数
:param n_clusters: 最大迭代数
"""
self.max_iter = max_iter
self.n_clusters = n_clusters
self.init_params = init_params
self.num = None
self.dim = None
self.X = None
self.Q = None
self.weight = None
self.covar = None
self.mu = None
self.labels = None
def _initialize_params(self, X):
"""
初试化模型参数,
:param X: 分类的数据集
:return:
"""
self.X = X # 分类的数据集
self.num = X.shape[0] # 样本数目
self.dim = X.shape[1] # 特征维度
self.Q = np.zeros((self.num, self.n_clusters)) # 初始化各高斯分布对观测数据的响应度矩阵
if self.init_params == "random":
self.weight = [1 / self.n_clusters] * self.n_clusters # 初始化各高斯分布的权重为聚类数目分之一
self.mu = np.random.uniform(0, 1, (self.n_clusters, self.dim)) * np.max(X, axis=0) # 随机产生均值向量
self.covar = np.array([np.identity(self.dim) for _ in range(self.n_clusters)]) # 随机产生协方差矩阵
if self.init_params == "kmeans":
kmeanmodel = cl.Mykmean(self.n_clusters, 20)
kmeanmodel.fit(self.X)
self.mu = kmeanmodel.center
self.weight = []
# print(kmeanmodel.labels)
for i in range(self.n_clusters):
temp = np.zeros(kmeanmodel.labels.shape)
temp[kmeanmodel.labels == i] = 1
self.weight.append(np.sum(temp) / self.num)
self.covar = np.zeros((self.n_clusters, self.dim, self.dim))
for i in range(self.n_clusters):
temp = self.X.copy()
# print(np.argwhere(kmeanmodel.labels == i).T.tolist())
temp = temp[np.argwhere(kmeanmodel.labels == i).T.tolist()[0], :]
self.covar[i, :, :] = np.cov(temp, rowvar=False)
# print(self.mu)
# print(self.covar)
# print(self.weight)
def e_step(self):
"""
e步,更新分模型对数据的响应度矩阵Q, 计算公式为
:return:
"""
for i in range(self.num):
q_i = []
for k in range(0, self.n_clusters):
postProb = multi_norm(self.X[i, :], self.mu[k, :], self.covar[k, :, :])
q_i.append(self.weight[k] * postProb)
self.Q[i, :] = np.array(q_i) / np.sum(q_i)
def m_step(self):
# update weight 更新权值矩阵
self.weight = np.mean(self.Q, axis=0)
# update mu 更新均值向量
temp = []
for k in range(self.n_clusters):
up = np.zeros(self.dim)
for j in range(self.num):
up += self.Q[j, k] * np.array(self.X[j, :])
down = np.sum(self.Q[:, k])
temp.append(up / down)
self.mu = np.array(temp)
# update covar
for k in range(self.n_clusters):
up = np.zeros((self.dim, self.dim))
for j in range(self.num):
x_mu = np.matrix(self.X[j, :] - self.mu[k, :])
# print(x_mu.T*x_mu)
up += self.Q[j, k] * (x_mu.T * x_mu)
# print(up)
down = np.sum(self.Q[:, k])
var = np.array(up / down)
self.covar[k, :, :] = var
def fit(self, X):
self.X = X
self._initialize_params(X)
while self.max_iter > 0:
# 初始化变量
# e-step
self.e_step()
# m-step
self.m_step()
self.max_iter -= 1
self.labels = np.argmax(self.Q, axis=1)
if __name__ == '__main__':
iris = load_iris()
# # 测试post_prob
# from scipy.stats import multivariate_normal
# mean = [0, 0,2]
# cov = [[1, 0,0], [0, 1,0],[0,0,1]]
# x = [18, 2, 1]
# var = multivariate_normal(mean, cov)
# print(var.pdf([18, 2,1]))
# print(multi_norm(np.array(x), np.array(mean), np.array(cov)))
# 鸢尾花的数据对比
from sklearn import mixture
gmm = mixture.GaussianMixture(n_components=3, init_params="kmeans").fit(iris.data)
labels = gmm.predict(iris.data)
label = pairing(iris.data, iris.target, labels)
acc = calcu_acc(iris.target, label)
print("GMM acc is :", acc)
mygmmmodel = GMM(3, 100, init_params="kmeans")
mygmmmodel.fit(iris.data, )
label = pairing(iris.data, iris.target, mygmmmodel.labels)
acc = calcu_acc(iris.target, label)
print("MYGMM acc is :", acc)
# 生成数据进行对比
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# Generate some data
from sklearn.datasets.samples_generator import make_blobs
X, y_true = make_blobs(n_samples=400, centers=4,
cluster_std=0.60, random_state=0)
X = X[:, ::-1] # flip axes for better plotting
gmm = mixture.GaussianMixture(n_components=4, init_params="kmeans").fit(X)
labels = gmm.predict(X)
fig, axs = plt.subplots(1,2)
ax1 = axs[0]
ax1.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis')
ax1.set_title("Scipy GMM")
model = GMM(4, 50, init_params="kmeans")
model.fit(X, )
my_label = model.labels
ax2 = axs[1]
ax2.scatter(X[:, 0], X[:, 1], c=my_label, s=40, cmap='viridis')
ax2.set_title("My GMM")
plt.show()
参数初始化问题
原型聚类一般都有一个局部极值的问题,因此初始化参数的选择就十分重要。在本例中,我选择了以kmeans聚类得到的聚类中心,以及计算每一类的协方差矩阵作为GMM的初试参数。在scipy库中,默认是采用了kmean的方法,当时没考虑到这个,用鸢尾花数据测试时浪费了不少时间。在都用了kmean初始化之后,可以看到结果基本一致。您也可以试试用random初试化得结果,可以init_params改为“random”。
附加:上传的代码中的kmeans聚类初试化参数是自已实现的,读者直接拷贝代码直接运行会报错。如果您想自己实现,可以写一个kmeans聚类或者参考我的另一篇博文kmeans聚类,或者您可以将_initialize_params方法中cl.Mykmeans()函数改为scipy库中的函数。为了保证numpy实现,我就自己不改啦,hah。
鸢尾花数据数据运行结果:
scipy gmm的准确率为:0.9667
自己实现的gmm的准确率为:0.9667
生成数据如下:
tips:ML新手,如果哪里写得不好,大家多多包涵。如果哪里有疑问,欢迎留言讨论。此外,我将继续有间断用numpy实现周志华《机器学习》的算法,欢迎大家关注。