技術標籤：聚類python機器學習

1 實驗目的
實現一個k-means演算法和混合高斯模型，並且用EM演算法估計模型中的引數。

2實驗要求
測試：
用高斯分佈產生k個高斯分佈的資料（不同均值和方差）（其中引數自己設定）。
（1）用k-means聚類，測試效果；
（2）用混合高斯模型和你實現的EM演算法估計引數，看看每次迭代後似然值變化情況，考察EM演算法是否可以獲得正確的結果（與你設定的結果比較）。

應用：可以UCI上找一個簡單問題資料，用你實現的GMM進行聚類。

3實驗環境
Windows 10,
Python 3.8.0,
Visual Studio Code

4.1實驗原理
本實驗分為兩部分：K-means和GMM，這兩部分都是屬於EM演算法，而EM演算法主要分為兩步：E步：求期望,M步：求最大似然。E步是調整分佈，M步是根據E步得到的分佈求得當前分佈下取到最大似然時的引數，然後更新引數，再次進入E步根據得到的引數調整新的分佈，如此往復迴圈，直到引數收斂。

4.2 K-means
原理及具體實現：
k-means聚類就是根據某種度量方式(常用歐氏距離，如歐氏距離越小，相關性越大)，將相關性較大的一些樣本點聚集在一起，一共聚成k個堆，每一個堆我們稱為一“類”。k-means的過程為：先在樣本點中選取k個點作為暫時的聚類中心，然後依次計算每一個樣本點與這k個點的距離，將每一個與距離這個點最近的中心點聚在一起，這樣形成k個類“堆”，求每一個類的期望，將求得的期望作為這個類的新的中心點。一直不停地將所有樣本點分為k類，直至中心點不再改變停止。
給定訓練樣本X={x1,x2,…,xm}，和劃分聚類的數量 k kk，給出一個簇劃分C = C 1 , C 2 , . . . , C k C=C_1,C_2,…,C_kC=C1,C2,…,Ck，E=i=1∑kx∈Ci∑∣∣x−μi∣∣22
其中，μ i = 1 ∣ C i ∣ ∑ x ∈ C i x i \mu_i=\frac 1 {|C_i|}\sum_{x\in C_i}x_iμi=∣Ci∣1∑x∈Cixi，它是簇C i C_iCi的均值向量。E EE刻畫了簇內樣本圍繞簇的均值向量的緊密程度，E EE越小簇內樣本的相似度越高。
具體迭代過程如下：使得該劃分的平方誤差E EE最小化，即使下式取最小值：根據輸入的超引數K KK首先初始化一些向量（可以從現有的向量中挑選），作為各簇的均值向量。根據初始化的均值向量給出訓練樣本的一個劃分，計算各個訓練樣本到各個均指向量的距離，找出距離最近的均值向量，並將該樣本分至該均值向量所代表的簇。根據新的簇劃分，重新計算每個簇的均值向量，如果新的均值向量與舊的均值向量差小於ε，則認為演算法收斂；否則，更新均值向量，回到第2步重新迭代求解。

4.3 GMM
原理及具體實現：初始化響應度矩陣γ, 協方差矩陣，均值和α，E步：定義響應度矩陣γ，其中γjk表示第j個樣本屬於第k類的概率，計算式如下

M步：將響應度矩陣視為定值，更新均值，協方差矩陣和α

5.1 結論
1.與K-means相比，GMM-EM似乎更容易陷入到區域性最優解當中，可能因為我此次實驗中的EM演算法採用的是隨機中心點才導致了該現象，我推測GMM-EM演算法對初值的敏感程度更高，正因如此我認為將K-Means得到的結果作為GMM-EM方法的初值會使得GMM-EM演算法更精準。
2.正如結論一所說將K-Means得到的結果作為GMM-EM方法的初值會使得GMM-EM演算法更精準。所以GMM-MM適合作為優化演算法，即適合對K-Means的分類結果進行進一步優化。

import numpy as np
import random
import matplotlib.pyplot as plt

weidu = 2

def rand_params(dim, k):
    """
    隨機初始引數
    """
    mu = np.random.rand(k, dim).reshape(k, dim)
    sigma = np.array([np.eye(dim)] * k).reshape(k, dim, dim)
    alpha = (np.ones(k) * (1.0 / k)).reshape(k, 1)
    return mu, sigma, alpha


def get_probability(x, mu, sigma, threshold=1e-8):
    """
    計算概率
    """
    n = mu.shape[1]
    if np.linalg.det(sigma) == 0:
        for i in range(sigma.shape[0]):
            sigma[i, i] += threshold
    p = np.exp(-0.5 * np.dot(np.dot(x - mu, np.linalg.pinv(sigma)), (x - mu).T))
    p = p / (np.power(2 * np.pi, n / 2.0) * np.sqrt(np.linalg.det(sigma)))
    return p


def gmm(x, k):  #k=3中心點的數量
    """
    GMM
    """
    x_size = np.shape(x)[0]   #點數（矩陣的行數）
    mu, sigma, alpha = rand_params(weidu, k)
    gamma = np.zeros((x_size, k))
    lld = np.array([])
    last_l_theta = 1e9
    t = 0
    for times in range(1000):
        prob = np.zeros((x_size, k))
        for i in range(x_size):
            for j in range(k):
                prob[i, j] = get_probability(x[i, :].reshape(1, -1), mu[j, :].reshape(1, -1), sigma[j])
        # E步
        for i in range(k):
            gamma[:, i] = alpha[i, 0] * prob[:, i]
        # 計算似然值
        l_theta = np.sum(np.log(np.sum(gamma, axis=1)))
        if np.abs(last_l_theta - l_theta) < 1e-10:
            t += 1
        else:
            t = 0
        if t == 10:
           # print('已收斂，迭代次數{}'.format(times + 1))
            break
        last_l_theta = l_theta
       # print(l_theta)
        lld = np.append(lld, l_theta)
        for i in range(x_size):
            gamma[i, :] /= np.sum(gamma[i, :])
        # M步
        alpha = (np.sum(gamma, axis=0) / x_size).reshape(k, 1)
        for i in range(k):
            nk = np.sum(gamma[:, i])
            mu[i, :] = np.dot(gamma[:, i].reshape(1, x_size), x) / nk
            tmp = np.zeros((weidu, weidu))
            for j in range(x_size):
                v = (x[j, :] - mu[i, :]).reshape(-1, 1)
                tmp += (gamma[j, i] * np.dot(v, v.T))
            sigma[i, :, :] = tmp / nk
    return gamma


def get_y(gamma):
    return np.argmax(gamma, axis=1)   #每行的最大值


def main():
        data = np.loadtxt('bezdekIris.csv', delimiter=',')
        dataset = data[:,0:weidu]
        print("*************************************")
        print('real data:')
        print(np.shape(data[np.where(data[:,-1]==1)])[0])
        print(np.shape(data[np.where(data[:,-1]==2)])[0])
        print(np.shape(data[np.where(data[:,-1]==3)])[0])
        print("*************************************")
        x = data[:, 0:weidu]
        y = data[:,-1]
        a=gmm(x, 3)
        b=gmm(x, 3)
        gmm_y = get_y(a)
        gmm_y = get_y(b)
        cluPoint=dataset[np.where(gmm_y==0)] 
        plt.scatter(cluPoint[:,0],cluPoint[:,1], c='red')
        cluPoint=dataset[np.where(gmm_y==1)] 
        plt.scatter(cluPoint[:,0],cluPoint[:,1], c='green')
        cluPoint=dataset[np.where(gmm_y==2)] 
        plt.scatter(cluPoint[:,0],cluPoint[:,1], c='black')
        print("*************************************")
        print('gmm result:')
        print(np.sum(gmm_y == 0))
        print(np.sum(gmm_y == 1))
        print(np.sum(gmm_y == 2))
        print("*************************************")
        plt.show()

if __name__ == '__main__':
    main()


    

import argparse
from math import inf

import matplotlib.pyplot as plt
import numpy as np

w = 2 #維數

# 生成K=3箇中心點
def centerGen(dataSet, k):
    n = np.shape(dataSet)[1]   #dataSet的列數
    cpoints = np.zeros((k, n)) #3*2
    # create centroid mat
    for j in range(0,w):
        # create random cluster centers, within bounds of each dimension
        min = np.min(dataSet[:, j])      #第J列的最小值
        cha = float(np.max(dataSet[:, j]) - min)    #第j列最大值與最小值之差
        cpoints[:, j] = min + cha*np.random.rand(k)  #生成中心的橫縱座標,np.random.rand(k)生成K=3個隨機資料
    return cpoints

# 求解距離
def distMeasure(A, B):
    return np.linalg.norm(A - B)   #求解歐氏距離

# K-means main method.
def kmeans(dataSet, k,centers):
    m = np.shape(dataSet)[0]  # 點數
    cluResult = np.zeros((m, 2))  # 30*2全0矩陣
    cluChanged = True
    while cluChanged:    #用於迭代，找到對應點
        cluChanged = False
        for i in range(m):
            minDist = inf
            minIndex = -1
            for j in range(0,k):
                # Step 1: 為每個點找到最近的中心
                y = j+1
                newDist = distMeasure(centers[j, :], dataSet[i, :])
                if newDist < minDist:
                    minDist = newDist   #找到距離點的最近的中心點距離
                    minIndex = y        #最近的點
            if cluResult[i, 0] != minIndex:   #如果點的標籤改變了
                cluChanged = True
                cluResult[i, :] = minIndex, minDist   #附上最新的標籤和距離中心點的距離
        # Step 2: 重新評估中心點
        for cent in range(0,k):
            ptsInClust = dataSet[np.where(cluResult[:, 0] == (cent+1))]
            centers[cent, :] = np.mean(ptsInClust, axis=0)   #np.mean計算均值
    return centers, cluResult   #返回K=3箇中心點，和每個點到中心點距離


def main():
    data = np.loadtxt('test.csv', delimiter=',')
    dataSet = data[:, 0:w]
    print("*************************************")
    print(np.shape(data[np.where(data[:,-1]==1)])[0])
    print(np.shape(data[np.where(data[:,-1]==2)])[0])
    print(np.shape(data[np.where(data[:,-1]==3)])[0])
    print("*************************************")
    k=3
    centers = centerGen(dataSet, k)  #初始的K=3箇中心點
    centers, cluResult = kmeans(dataSet, k,centers)
    print('Centers:', centers)  #中心點的座標
    #print('Clustering result:', cluResult)   #改過標籤的點集
    # Drawing
    colors=['red', 'pink', 'green']
    print("*************************************")
    print(np.shape(cluResult[np.where(cluResult[:,0]==1)])[0])
    print(np.shape(cluResult[np.where(cluResult[:,0]==2)])[0])
    print(np.shape(cluResult[np.where(cluResult[:,0]==3)])[0])
    print("*************************************")
    for i in range(0,3):
        center=centers[i]
        cluPoint=dataSet[np.where(cluResult[:,0]==(i+1))]  #np.where(cluResult[:,0]==i)返回標籤為I的資料的序號
        plt.scatter(center[0],center[1],c=colors[i],marker='*')
        plt.scatter(cluPoint[:,0],cluPoint[:,1], c=colors[i])
    plt.show()


if __name__ == '__main__':
    main()