一，介绍
学习混合高斯，先要了解几个概念：
1，协方差:
协方差是对两个随机变量联合分布线性相关程度的一种度量。两个随机变量越线性相关，协方差越大，完全线性无关，协方差为零。
根据数学期望的性质：
推导协方差为：
2，协方差矩阵：
对多维随机变量X=[X1,X2,X3,...,Xn]，我们往往需要计算各维度两两之间的协方差，这样各协方差组成了一个n×n的矩阵，称为协方差矩阵：
百度百科给出的实例：
3，服从高斯分布的概率密度函数：
可以看出高斯分布由均值向量μ和协方差矩阵∑这两个参数确定，故将概率密度函数记为：p(x|μ,∑)
可定义高斯混合分布为：
给定一组数据，假设该数据由多个高斯分布产生，现在我们要估计这些高斯分布的参数，以及每个样本属于每个高斯分布的概率，那么根据样本推测出元GMM的概率分布就可以了。具体求解借助EM算法可得。
二，公式推导
假设有K个样本，服从C个高斯分布，随机抽取一个样本X得到：
（1）
其中，
为协方差矩阵，
是聚类的中心（向量均值）
根据上式，可以写出关于样本值（第i个样本）的概率密度函数：
（2）
再根据贝叶斯定理：
该分布由c个混合成分组成，每个混合对应一个高斯分布，P（Cj）为对应的“混合系数”，且
=1。
我们定义m个观测样本的对数似然函数:
为了得到最大对数似然函数，需要对分别对和求偏导并令其为0：
（1）求:
将公式（2）代入，并且除去第K个分量其余可以看做常数，得到：
由于，只有第K个样本与有关，其他的都可以视为常数，上式变为：
再利用对数复合函数求导得到：
利用贝叶斯公式，整理获得：
（3）接着，求解，
：
求偏导：
（4）
令偏导为0，公式（4）代入公式（3），又知
不为0，可得：
(2)求
类似与求
利用矩阵求导的性质:
再继续求解得：
我们假设K的位置再行i和列j上，那么可知，∑除了在位置（i，j）上为1以外，其他的都为零，上式结果为：
代入上式，获得：
(3)求k个分量的比例的问题
在（1）和（2）中，我们已经获取了和，我们要求K个分量最佳比例值，相当于优化：
我们利用拉格朗日定理求解：
综上所述得到：
将其代入第K个分量，获得：
我们总结下算法流程：
第一步，初始化模型参数；
第二步，计算高斯混合分布
第三步，计算Ck的后验分布：
第四步，计算新的均值向量、协方差矩阵、混合系数
第五步，如果不满足条件，跳转到第二步继续计算。
第六步，满足条件，确定均值向量中心，获得分类模型
二，代码实现
数据：
0.697 0.46
0.774 0.376
0.634 0.264
0.608 0.318
0.556 0.215
0.403 0.237
0.481 0.149
0.437 0.211
0.666 0.091
0.243 0.267
0.245 0.057
0.343 0.099
0.639 0.161
0.657 0.198
0.36 0.37
0.593 0.042
0.719 0.103
0.359 0.188
0.339 0.241
0.282 0.257
0.748 0.232
0.714 0.346
0.483 0.312
0.478 0.437
0.525 0.369
0.751 0.489
0.532 0.472
0.473 0.376
0.725 0.445
0.446 0.459import numpy as np
import matplotlib.pyplot as plt
# 预处理数据
def loadDataSet(filename):
dataSet = []
fr = open(filename)
for line in fr.readlines():
curLine = line.strip().split('\t')
fltLine = list(map(float, curLine))
dataSet.append(fltLine)
return dataSet
# 高斯分布的概率密度函数
def prob(x, mu, sigma):
n = np.shape(x)[1]
expOn = float(-0.5 * (x - mu) * (sigma.I) * ((x - mu).T))
divBy = pow(2 * np.pi, n / 2) * pow(np.linalg.det(sigma), 0.5) # np.linalg.det 计算矩阵的行列式
return pow(np.e, expOn) / divBy
# EM算法
def EM(dataMat, maxIter=50,c=3):
m, n = np.shape(dataMat)
# 第一步，初始化参数
alpha = [1 / 3, 1 / 3, 1 / 3] # 初始化参数 alpha1=alpha2=alpha3=1/3
mu = [dataMat[5, :], dataMat[21, :], dataMat[26, :]] # 初始化聚类中心 mu1=x6,mu2=x22,mu3=x27
sigma = [np.mat([[0.1, 0], [0, 0.1]]) for x in range(3)] # 初始化协方差矩阵
gamma = np.mat(np.zeros((m, c)))
for i in range(maxIter):
for j in range(m):
sumAlphaMulP = 0
for k in range(c):
gamma[j, k] = alpha[k] * prob(dataMat[j, :], mu[k], sigma[k]) # 计算混合成分生成的后验概率
sumAlphaMulP += gamma[j, k] # 第二步，计算高斯混合分布
for k in range(c):
gamma[j, k] /= sumAlphaMulP # 第三步，计算后验分布
sumGamma = np.sum(gamma, axis=0)
for k in range(c):
mu[k] = np.mat(np.zeros((1, n)))
sigma[k] = np.mat(np.zeros((n, n)))
for j in range(m):
mu[k] += gamma[j, k] * dataMat[j, :]
mu[k] /= sumGamma[0, k] # 第四步，计算均值向量
for j in range(m):
sigma[k] += gamma[j, k] * (dataMat[j, :] - mu[k]).T *(dataMat[j, :] - mu[k])
sigma[k] /= sumGamma[0, k] # 第四步，计算协方差矩阵
alpha[k] = sumGamma[0, k] / m # 第四步，计算混合系数
return gamma
# init centroids with random samples
def initCentroids(dataMat, k):
numSamples, dim = dataMat.shape
centroids = np.zeros((k, dim))
for i in range(k):
index = int(np.random.uniform(0, numSamples))
centroids[i, :] = dataMat[index, :]
return centroids
def gaussianCluster(dataMat):
m, n = np.shape(dataMat)
centroids = initCentroids(dataMat, m)
clusterAssign = np.mat(np.zeros((m, 2)))
gamma = EM(dataMat)
for i in range(m):
clusterAssign[i, :] = np.argmax(gamma[i, :]), np.amax(gamma[i, :])
for j in range(m):
pointsInCluster = dataMat[np.nonzero(clusterAssign[:, 0].A == j)[0]]
centroids[j, :] = np.mean(pointsInCluster, axis=0) # 第五步，确定各均值中心，获得分类模型
return centroids, clusterAssign
def showCluster(dataMat, k, centroids, clusterAssment):
numSamples, dim = dataMat.shape
if dim != 2:
print("Sorry! I can not draw because the dimension of your data is not 2!")
return 1
mark = ['or', 'ob', 'og', 'ok', '^r', '+r', 'sr', 'dr', '<r', 'pr']
if k > len(mark):
print("Sorry! Your k is too large!")
return 1
for i in range(numSamples):
markIndex = int(clusterAssment[i, 0])
plt.plot(dataMat[i, 0], dataMat[i, 1], mark[markIndex])
mark = ['Dr', 'Db', 'Dg', 'Dk', '^b', '+b', 'sb', 'db', '<b', 'pb']
for i in range(k):
plt.plot(centroids[i, 0], centroids[i, 1], mark[i], markersize=12)
plt.show()
if __name__=="__main__":
dataMat = np.mat(loadDataSet('watermelon4.0.txt'))
centroids, clusterAssign = gaussianCluster(dataMat)
print(clusterAssign)
showCluster(dataMat, 3, centroids, clusterAssign)
