EM算法与高斯混合聚类

最新推荐文章于 2024-12-28 10:18:57 发布

原创最新推荐文章于 2024-12-28 10:18:57 发布 · 1.7k 阅读

4 ·

CC 4.0 BY-SA版权

文章标签：

#机器学习 #聚类 #高斯混合聚类 #EM算法 #Python

机器学习专栏收录该内容

18 篇文章

订阅专栏

本文详细介绍了EM算法及其在高斯混合聚类中的应用。首先阐述了EM算法的基本思想，强调了其对初值的敏感性和可能存在的局部最优问题。接着，讲解了在高斯混合聚类中如何运用EM算法，包括E步和M步的具体计算过程。最后，文章通过在西瓜数据集上的实际操作，展示了如何设置和运行高斯混合聚类模型，其中设置的混合成分个数为3。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

EM算法

用Y表示观测随机变量的数据，Z表示隐随机变量的数据，Y和Z连在一起成为完全数据，观测Y又称为不完全数据。假设给定观测数据Y，其概率分布是 $P(Y|\theta)$ ，其中 $\theta$ 是要估计的模型参数，完全数据的对数似然函数为 $logP(Y,Z|\theta)$ ，EM算法通过迭代求对数似然函数的极大似然估计，每次迭代包括两步：E步，求期望；M步，求极大化。
算法步骤如下：
1.选择参数的初值 $\theta^{(0)}$ ，开始迭代
2.E步：记 $\theta{(i)}$ 为第i次迭代参数 $\theta$ 的估计值，在第i+1次迭代的E步，计算
$Q(\theta,\theta^{(i)})=E_Z[logP(Y,Z)|Y,\theta^{(i)}]$
3.M步：求使 $Q(\theta,\theta^{(i)})$ 极大化的 $\theta$ ，确定第i+1次的迭代参数估计值 $\theta^{(i+1)}$
4.重复第2和第3步直到收敛。
要注意的是EM算法对初值是敏感的，算法不一定可能会收敛到局部最优值。

高斯混合聚类

上一篇中已经讲了高斯混合聚类的原理，这里主要说明高斯混合聚类用EM算法求解的过程。
首先是E步，写出Q函数的形式，然后计算函数中的 $E(\gamma_{jk})$
$\hat{\gamma} = E(\gamma_{jk}|y,\theta)\\ =P(\gamma_{jk}=1|y,\theta)\\={\alpha_k*\phi(y_j|\theta_k)\over \sum_{k=1}^K\alpha_k\phi(y_j|\theta_k)}$
随后是M步，M步要求Q函数的极大值，即
$\theta^{(i+1)}=arg\ max_{\theta}Q(\theta,\theta{(i)})$
求得 $\mu,\Sigma,\alpha$ 的估计值，具体来说就是求偏导并令导数为0，可以得到
$\hat \mu_k={\sum_{j=1}^N \hat \gamma_{jk}y_j \over \sum_{j=1}^N \hat \gamma_{jk}}$
$\hat\Sigma_k={\sum_{j=1}^N \hat \gamma_{jk}(x_j-\mu_k)(x_j-\mu_k)^T\over \sum_{j=1}^N \hat \gamma_{jk}}$
$\hat\alpha_k={\sum_{j=1}^N \hat \gamma_{jk}\over N}$
算法在达到最大迭代轮数或者更新量很小的时候停止。

高斯混合聚类代码

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

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

$\alpha_1=\alpha_2=\alpha_3=1/3$

μ 1 = x 6, μ 2 = x 22, μ 3 = x 27

$\mu_1=x_6,\mu_2=x_{22},\mu_3=x_{27}$

Σ 1 = Σ 2 = Σ 3 = [0.1 0 0 0.1]

$\Sigma_1=\Sigma_2=\Sigma_3=\begin{bmatrix} 0.1&0\\ 0&0.1\end{bmatrix}$
默认最大迭代次数为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)