机器学习第九篇：EM算法及其在高斯混合模型和K-means中的应用

最新推荐文章于 2024-07-14 08:21:19 发布

一个小路人

最新推荐文章于 2024-07-14 08:21:19 发布

阅读量630

点赞数 3

分类专栏：机器学习文章标签： kmeans EM 高斯混合模型

本文链接：https://blog.youkuaiyun.com/h_q_c_/article/details/89889069

版权

机器学习专栏收录该内容

12 篇文章

订阅专栏

1.EM算法：

EM算法产生的前提：当概率模型的变量都是观测变量的时候可以直接用极大似然估计法求得我们的最优解，但是当模型含有隐变量的时候，就无法直接求得最优解，而此时就产生了我们的EM算法。

假设Y表示可观测的变量，Z(跟Y是存在关联的)表示不可观测的数据而我们最终的目的是通过迭代来求解出 $L(\theta )=logP(Y|\theta)$ ,每次迭代分两步：E步，求期望；M步，求极大化。

算法：

选择初始参数 $\theta$ ，开始迭代；
E步：第i次迭代参数 $\theta$ 的估计值，在第i+1次迭代的E步，计算：

$Q(\theta,\theta^i) = E_Z[logP(Y,Z|\theta)|Y,\theta^i]$

$=\sum_ZlogP(Y,Z|\theta)P(Z|Y,\theta^i)$

M步：求使 $Q(\theta,\theta^i)$ 极大化的 $\theta$ ,确定i+1次迭代的估计值 $\theta^{(i+1)}$

$\theta^{i+1} = argmaxQ(\theta,\theta^i)$

重复前两步直至收敛

这样看的话可能比较抽象其实就是找出Q，对Q求解相应偏导即可，我们就以两个例子来解释这个算法。

2.EM算法在高斯混合模型学习中的应用：

高斯混合模型是指具有如下形式的概率分布模型：

$P(y|\theta) = \sum^K_{k=1}\alpha _k\o (y|\theta_k)$

其中， $\alpha_k$ 是系数， $\alpha_k\geq 0,\sum\alpha_k = 1$ ； $\o (y|\theta_k)$ 是高斯分布密度， $\theta_k = (\mu _k,\sigma _k^2)$ ,

$\o (y|\theta_k) = \frac{1}{\sqrt{2\pi}\sigma _k}exp(-\frac{(y-\mu_k)^2}{2\sigma_k^2})$

然后明确我们的隐变量，写出对数似然函数：

根据这段话，直接写出我们的Q函数：

$Q(\theta,\theta^i) = E_Z[logP(Y,Z|\theta)|Y,\theta^i]$

$= \sum_{k=1}^K\{n_klog\alpha_k+\sum_{k=1}^N \hat{\gamma}_{jk}[log(\frac{1}{\sqrt{2\pi}})-log\sigma_k -\frac{1}{2\sigma_k^2}(y_i-\mu_k)^2]\}$

然后利用Q分别对mu、sigma、alpha求偏导。

高斯混合模型参数估计的EM算法如下：

取参数初始值开始迭代
E步：

$\hat{\gamma _{jk}} = \frac{\alpha_k\o(y_j|\theta_k)}{\sum_{k=1}^K\alpha_k\o(y_j|\theta_k)}$ , j= 1,2,...,N ; k = 1,2,...,K

M步：

$\hat{\mu}_k = \frac{\sum_{j=1}^N\hat{\gamma}_{jk}y_j}{\sum_{j=1}^N\hat{\gamma}_{jk}}$

$\hat{\sigma}_k^2 = \frac{\sum_{j=1}^N\hat{\gamma}_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^N\hat{\gamma}_{jk}}$

$\hat{\alpha_k} = \frac{\sum_{j=1}^N\hat{\gamma}_{jk}}{N}$

重复前两步直至收敛

算法就是这么多，比起前面单说EM算法要容易理解的多，下面用程序具体实现一下这样的一个过程：

from numpy import *
import matplotlib.pyplot as plt
import math



def get_set_data(file_name):
    fd = open(file_name)
    set_data = []

    for line in fd.readlines():
        line_data = line.strip().split()
        set_data.append([float(line_data[0]), float(line_data[1])])

    return set_data

def gaussion_distribution(x, miu, sigma):
    n = shape(x)[1]
    sigma_mat = mat(sigma)
    expOn = float(-0.5*(x - miu)*(sigma_mat.I)*((x - miu).T))
    divBy = pow(2*pi, 1/2)*pow(linalg.det(sigma_mat), 0.5)
    return pow(e, expOn)/divBy

def gaussion_em(data_in, alpha, miu, sigma, iterator_num):
    K = shape(alpha)[1]
    N,n = shape(data_in)
    gama = mat(zeros((K, N)))
    temp_value = mat(zeros((K, N)))

    for s in range(iterator_num):   #迭代iterator_num次
        for j in range(N):
            for k in range(K):
               temp_value[k, j] =  alpha[:, k] * gaussion_distribution(data_in[j, :], miu[k], sigma[k])
        for j in range(N):
            sum = 0
            for k in range(K):
                sum += temp_value[k, j]
            for k in range(K):
                gama[k, j] = temp_value[k, j] / sum

        for k in range(K):
            sum0 = mat(zeros((1, n)))
            sum2 = mat(zeros((n, n)))
            sum1 = 0
            for j in range(N):
                sum0[0, :] += gama[k, j] * data_in[j, :]
                sum1 += gama[k, j]
            miu[k] = sum0[0, :] / sum1
            for j in range(N):
                sum2 += gama[k, j] * (data_in[j, :] - miu[k]).T * (data_in[j, :] - miu[k])

            sigma[k] = sum2 / sum1
            alpha[:, k] = sum1 / N
        
    return alpha,miu,sigma,gama

def get_classify(gama):
    label = 0

    K,N = shape(gama)
    labels = mat(zeros((N, 1)))
    for i in range(N):
        index_max = argmax(gama[:, i])  #计算每一个观测变量最有可能属于第几个高斯分量
        labels[i, 0] = index_max
    return labels

def show_experiment_plot(data_mat, labels):
    N = shape(data_mat)[0]

    for i in range(N):
        if(labels[i, 0] == 0):
            plt.plot(data_mat[i, 0], data_mat[i, 1], "ob")
        elif(labels[i, 0] == 1):
            plt.plot(data_mat[i, 0], data_mat[i, 1], "or")
        else:
            plt.plot(data_mat[i, 0], data_mat[i, 1], "og")

    plt.xlabel("X")
    plt.ylabel("Y")
    plt.show()


if __name__ == "__main__":
    set_data = get_set_data("1.txt")
    alpha_mat = mat([1/3, 1/3, 1/3])
    data_mat = mat(set_data)
    miu_mat = [data_mat[5, :], data_mat[21, :], data_mat[26, :]]
    sigma_mat = [mat([[0.1, 0], [0, 0.1]]) for x in range(3)]
    #print(data_mat)
    alpha,miu,sigma,gama = gaussion_em(data_mat, alpha_mat, miu_mat, sigma_mat, 50)
    print("alpha = ", alpha)
    print("\n")
    print("miu = ", miu)
    print("\n")
    print("sigma = ", sigma)
    print("\n")
    print("gama = ", gama)
    print("\n")
    labels = get_classify(gama)
    show_experiment_plot(data_mat, labels)



#实验数据
0.697 0.460
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.360 0.370
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.459

程序还是很简单的，采用了西瓜书kmeans的数据来测试的。

3.K-means算法：

在介绍EM的kmeans之前，先简单介绍一下kmeans算法：

kmeans是很简单的无监督聚类算法，其伪代码如下：

创建k个点作为起始质心（经常是随机选择）

当任意一个点的簇分配结果发生改变时

对数据集中的每个数据点

对每个质心

计算质心与数据点直接的距离

将数据点分配到距其最近的簇

对每一个簇，计算簇中所有点的均值并将均值作为质心

from numpy import *

def loadDataSet(fileName):      #general function to parse tab -delimited floats
    dataMat = []                #assume last column is target value
    fr = open(fileName)
    for line in fr.readlines():
        curLine = line.strip().split('\t')
        fltLine = [float(i) for i in curLine] #map all elements to float()
        dataMat.append(fltLine)
    return dataMat

def distEclud(vecA, vecB):
    return sqrt(sum(power(vecA - vecB, 2))) #la.norm(vecA-vecB)

def randCent(dataSet, k):
    n = shape(dataSet)[1]
    centroids = mat(zeros((k,n)))#create centroid mat
    for j in range(n):#create random cluster centers, within bounds of each dimension
        minJ = min(dataSet[:,j]) 
        rangeJ = float(max(dataSet[:,j])-minJ)
        centroids[:,j] = mat(minJ + rangeJ * random.rand(k,1))
    return centroids
    
def kMeans(dataSet, k, distMeas=distEclud, createCent=randCent):
    m = shape(dataSet)[0]
    clusterAssment = mat(zeros((m,2)))#create mat to assign data points 
                                      #to a centroid, also holds SE of each point
    centroids = createCent(dataSet, k)
    clusterChanged = True
    while clusterChanged:
        clusterChanged = False
        for i in range(m):#for each data point assign it to the closest centroid
            minDist = inf; minIndex = -1
            for j in range(k):
                distJI = distMeas(centroids[j,:],dataSet[i,:])
                if distJI < minDist:
                    minDist = distJI; minIndex = j
            if clusterAssment[i,0] != minIndex: clusterChanged = True
            clusterAssment[i,:] = minIndex,minDist**2
        print (centroids)
        for cent in range(k):#recalculate centroids
            ptsInClust = dataSet[nonzero(clusterAssment[:,0].A==cent)[0]]#get all the point in this cluster
            centroids[cent,:] = mean(ptsInClust, axis=0) #assign centroid to mean 
    return centroids, clusterAssment

kmeans与EM算法：

K-means的目标是将样本集划分为K个簇，使得同一个簇内样本的距离尽可能小，不同簇之间的样本距离尽可能大，即最小化每个簇中样本与质心的距离。K-means属于原型聚类(prototype-based clustering)，原型聚类指聚类结构能通过一组原型刻画，而原型即为样本空间中具有代表性的点。在K-means中，这个原型就是每个簇的质心 $\mu$ ;