機器學習第九篇：EM算法及其在高斯混合模型和K-means中的應用

1.EM算法：

EM算法産生的前提：當機率模型的變量都是觀測變量的時候可以直接用極大似然估計法求得我們的最優解，但是當模型含有隐變量的時候，就無法直接求得最優解，而此時就産生了我們的EM算法。

假設Y表示可觀測的變量，Z(跟Y是存在關聯的)表示不可觀測的資料而我們最終的目的是通過疊代來求解出

,每次疊代分兩步：E步，求期望；M步，求極大化。

算法：

• 選擇初始參數

  

  ，開始疊代；
• E步：第i次疊代參數

  

  的估計值，在第i+1次疊代的E步，計算：

• M步：求使

  

  極大化的

  

  ,确定i+1次疊代的估計值

  

• 重複前兩步直至收斂

這樣看的話可能比較抽象其實就是找出Q，對Q求解相應偏導即可，我們就以兩個例子來解釋這個算法。

2.EM算法在高斯混合模型學習中的應用：

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

其中，

是系數，

；

是高斯分布密度，

,

然後明确我們的隐變量，寫出對數似然函數：

根據這段話，直接寫出我們的Q函數：

然後利用Q分别對mu、sigma、alpha求偏導。

高斯混合模型參數估計的EM算法如下：

• 取參數初始值開始疊代
• E步：

   ,         j= 1,2,...,N ; k = 1,2,...,K

• M步：

• 重複前兩步直至收斂

算法就是這麼多，比起前面單說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中，這個原型就是每個簇的質心

;

從EM算法的觀點來看，K-means的參數就是每個簇的質心

，隐變量為每個樣本的所屬簇。如果事先已知每個樣本屬于哪個簇，則直接求平均即可得到

。但現實中不知道的情況下，則需要運用EM的思想：

假設要k個簇，先随機標明k個點作為質心

；

• 固定

  

  ，将樣本劃分到距離最近的

  

  所屬的簇中。若用

  

  表示第n個樣本所屬的簇，則

• 固定

  

  ，根據上一步的劃分結果重新計算每個簇的質心。由于我們的目标是最小化每個簇中樣本與質心的距離，可将目标函數表示為

要最小化J則對

求導得:

則有：

即簇中每個樣本的均值向量。

​上面兩步分别更新

和

就對應了EM算法中的E步和M步。和EM算法一樣，K-means每一步都最小化目标函數 J，因而可以保證收斂到局部最優值，但在非凸目标函數的情況下不能保證收斂到全局最優值。