高斯混合模型聚类_高斯混合模型及EM算法

EM(expectationmaximization algorithm)算法是一种迭代算法，1977年由Dempster等人总结提出， 用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望；M步，求极大，所以这一算法称为期望极大算法，简称EM算法。 一、EM算法的推导 用X=(x_1，x_2，…，x_n)表示观测数据，H=(h_1，h_2，…，h_n)表示隐变量数据，隐变量数据h_i都服从隐变量分布Z(z_1，z_2，…，z_n)，即为隐变量数据的全部取值集合，那么，X与Z一起成为完全数据。观测数据可以直接获得并使用，但是隐变量取值并不明确，所以要根据观测数据的对数似然函数极大化来求得最优参数。 1.E步

(1)计算观测数据的似然函数

(2)计算对数似然函数  

 

EM算法通过迭代逐步极大化对数似然函数 来逐步优化参数 。在第m次迭代得到参数的一个估计值 ，并且希望在第m+1次迭代得到的参数估计值 能使得 对 数似然函数 能再次增加， 就这样逐步使得 对 数似然函数 达到极大值。因为 对 数似然函数 的式子已经固定，现在考虑如何保证m+1次迭代时 对 数似然函数 可以再增加，EM算法用了一个比较巧妙的方法，那就是求 对 数似然函数 的下界 ，只要保证每次迭代时 对 数似然函数 的下界取极大，就可以保证 对 数似然函数 一直在增加。 (3)利用Jensen不等式求解第m+1次迭代时 对 数似然函数 的下界 先来简单介绍Jensen不等式：如果f(x)是凹函数， x_i(i=1，2，…)来自f(x)的定义域，λ_i(i=1，2，…) 大于等于0 且∑λ_i =1，则有如下不等式成立： 当且仅当x_i为常数时等号成立。如果f(x)是凸函数则不等号反向。 如果f(x)为凹函数时，将x_i看作是随机变量X得取值，将λ_i看作随机变量取值x_i得概率，则E[f(x)]大于等于 f[E(x)]，当且仅当随机变量X是常数时等号成立。当f(x)是凸函数时不等号反向。 将Jensen不等式应用于对数似然函数，构造λ_i。 因为Jensen不等式成立有如下约束条件：                                且要使等号成立则必须有：                                     所以：

故有：

2.M步

三、EM算法在高斯混合模型学习中的应用

EM算法的一个重要应用是高斯混合模型的参数估计。高斯混合模型(Gaussian Mixed Model，GMM)是一种常见的聚类算法，在图像分割、对象识别、视频分析等方面均有应用，对于任意给定的数据样本集合，根据其分布概率，可以计算每个样本数据向量的概率分布，从而根据概率分布对其进行分类。高斯混合模型是指具有如下形式的概率分布模型：

  

其中，α_k是系数，α_k大于等于0，∑α_k =1；Φ(x|θ_k)是高斯分布密度，θ_k=(μ_k，σ²_k)则可得第k个分模型： 由此看出，高斯混合模型混合的基本分布就是高斯分布。 假设观测数据x_1，x_2，…，x_N由高斯混合模型生成，其中，θ=(α_1，α_2，…，α_k；θ_1，θ_2，…，θ_k)，下面用EM算法估计高斯混合模型的参数θ。 1. 明确隐变量与完全数据的对数似然函数 设想观测数据x_j，j=1，2，…，N是这样产生的：首先依概率α_k选择第k个高斯分布分模型Φ(x|θ_k)，然后依第k个分模型的概率分布Φ(x|θ_k)生成观测数据x_j。这时的观测数据是已知的；但是反应观测数据x_j来自第k个分模型的数据是未知的，以隐变量γ_jk表示，其定义为： γ_jk是0-1随机变量。有了观测数据与未观测数据，那么完全数据是(x_j，γ_j1，γ_j2，…，γ_jk)，j=1，2，…，N。于是，可以得到完全数据的似然函数：

进一步得到完全数据的对数似然函数：

2.EM算法的E步—确定Q函数 是当前模型参数下第j个观测数据来自第k个分模型的概率，称为分模型k对观测数据x_j的响应度。 3. 确定EM算法的M步 四、高斯混合模型参数估计的EM算法 输入：观测数据x_j，j=1，2，…，N与高斯混合模型。 输出：高斯混合模型参数。 (1)取参数的初始值开始迭代。 (2)E步：依据当前模型参数，计算分模型k对观测数据x_j的响应度：

(3)M步：计算新一轮迭代的模型参数。 (4)重复第2、3步，直到收敛。 五、案例分析 已知观测数据-67，-48，6，8，14，16，23，24，28，29，41，49，56，60，75，估计两个分量的高斯混合模型的5个参数。 初值设置为σ₁=1，σ₂=1，μ₁=0.5，μ₂=0.5，α₁=0.5，α₂=0.5。

import numpy as npfrom scipy.stats import normy = np.array([-67, -48, 6, 8, 14, 16, 23, 24, 28, 29, 41, 49, 56, 60, 75])K = 2  # 两个高斯N = 15  # y有15个数据# 参数初始化mu = np.array([0.5, 0.5])sigma = np.array([1.0, 1.0]) * 10alpha = np.array([0.5, 0.5])for i in range(10):    gm = np.zeros((N, K))    # E 步    for j in range(N):        for k in range(K):            gm[j, k] = alpha[k] * norm(mu[k], sigma[k]).pdf(y[j]) #使用scipy实现高斯分布        gm[j, :] /= sum(gm[j, :])  # gm[j,:] = gm[j,:] /sum(gm[j,:])    # M 步    mu2 = y.dot(gm) / sum(gm)    alpha2 = sum(gm) / N    sigma2 = np.zeros((2,))    sigma2[0] = sum(gm[:, 0] * (y - mu[0]) ** 2) / sum(gm[:, 0])    sigma2[1] = sum(gm[:, 1] * (y - mu[1]) ** 2) / sum(gm[:, 1])    #判断是否收敛    if sum((mu - mu2) ** 2 + (sigma - sigma2) ** 2 + (alpha - alpha2) ** 2) < 0.01:        break    mu = mu2    sigma = sigma2    alpha = alpha2    print("第%d次迭代\nalpha_0=%f,mu_0=%f,sigma_0=%f\nalpha_1=%f,mu_1=%f,sigma_1=%f\n" %(i+1, alpha[0], mu[0], sigma[0], alpha[1], mu[1], sigma[1]))