EM算法与高斯混合聚类

EM算法

用Y表示观测随机变量的数据，Z表示隐随机变量的数据，Y和Z连在一起成为完全数据，观测Y又称为不完全数据。假设给定观测数据Y，其概率分布是P(Y|θ)，其中θ是要估计的模型参数，完全数据的对数似然函数为logP(Y,Z|θ)，EM算法通过迭代求对数似然函数的极大似然估计，每次迭代包括两步：E步，求期望；M步，求极大化。

算法步骤如下：

1.选择参数的初值θ(0)，开始迭代

2.E步：记θ(i)为第i次迭代参数θ的估计值，在第i+1次迭代的E步，计算

Q(θ,θ(i))=EZ[logP(Y,Z)|Y,θ(i)]

3.M步：求使Q(θ,θ(i))极大化的θ，确定第i+1次的迭代参数估计值θ(i+1)

4.重复第2和第3步直到收敛。

要注意的是EM算法对初值是敏感的，算法不一定可能会收敛到局部最优值。

高斯混合聚类

上一篇中已经讲了高斯混合聚类的原理，这里主要说明高斯混合聚类用EM算法求解的过程。

首先是E步，写出Q函数的形式，然后计算函数中的E(γjk)

γ^=E(γjk|y,θ)=P(γjk=1|y,θ)=αk∗ϕ(yj|θk)∑Kk=1αkϕ(yj|θk)

随后是M步，M步要求Q函数的极大值，即

θ(i+1)=arg maxθQ(θ,θ(i))

求得μ,Σ,α的估计值，具体来说就是求偏导并令导数为0，可以得到

μ^k=∑Nj=1γ^jkyj∑Nj=1γ^jk

Σ^k=∑Nj=1γ^jk(xj−μk)(xj−μk)T∑Nj=1γ^jk

α^k=∑Nj=1γ^jkN

算法在达到最大迭代轮数或者更新量很小的时候停止。

高斯混合聚类代码

算法在西瓜数据集4.0_1上运行，令高斯混合成分的个数k=3，初始化时，将模型参数初始化为

α1=α2=α3=1/3

μ1=x6,μ2=x22,μ3=x27

Σ1=Σ2=Σ3=[0.1000.1]

默认最大迭代次数为50

from numpy import *

# 高斯混合聚类

# 预处理数据
def loadData(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 = shape(x)[1]
expOn = float(-0.5*(x-mu)*(sigma.I)*((x-mu).T))
divBy = pow(2*pi, n/2)*pow(linalg.det(sigma), 0.5)
return pow(e, expOn)/divBy

# EM算法
def EM(dataMat, maxIter=50):
m, n = shape(dataMat)
# 初始化各高斯混合成分参数
alpha = [1/3, 1/3, 1/3]
mu = [dataMat[5, :], dataMat[21, :], dataMat[26, :]]
sigma = [mat([[0.1, 0], [0, 0.1]]) for x in range(3)]
gamma = mat(zeros((m, 3)))
for i in range(maxIter):
for j in range(m):
sumAlphaMulP = 0
for k in range(3):
gamma[j, k] = alpha[k]*prob(dataMat[j, :], mu[k], sigma[k])
sumAlphaMulP += gamma[j, k]
for k in range(3):
gamma[j, k] /= sumAlphaMulP
sumGamma = sum(gamma, axis=0)
for k in range(3):
mu[k] = mat(zeros((1, n)))
sigma[k] = mat(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
#print(mu)
return gamma

def gaussianCluster(dataMat):
m, n = shape(dataMat)
# 每个样本的所属的簇，以及分到该簇对应的响应度
clusterAssign = mat(zeros((m, 2)))
gamma = EM(dataMat)
for i in range(m):
# amx返回矩阵最大值，argmax返回矩阵最大值所在下标
clusterAssign[i,:] = argmax(gamma[i, :]), amax(gamma[i, :])
return clusterAssign

dataMat = mat(loadData('watermelon4_1.txt'))
clusterAssign = gaussianCluster(dataMat)
#print(clusterAssign)