HIT机器学习Lab3

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()