GMM算法的实现

GMM算法的实现

程序结构

main.py

import os
from numpy import testing
import numpy as np
import DataReader, GMM
import matplotlib.pyplot as plt
from matplotlib import animation


def main():

    # # 读入数据并分割成训练集和测试集
    # dataReader = DataReader.DataReader()
    # X = dataReader.readDataX("cluster.dat")
    # trainData, testData = dataReader.sampleSplit(X)

    # 直接从文件中读入分割好的训练集和测试集
    dataReader = DataReader.DataReader()
    trainData, testData = dataReader.readDataX("train.dat"), dataReader.readDataX("test.dat")

    # 构建GMM模型，指定类簇数和数据维度，进行训练
    model = GMM.GMM(n_components=4, n=2)
    model.train(trainData)

    # 绘制train过程中的对数似然变化图
    model.plotLogLikelihood()
    plt.savefig("./img/LogLikelihood.png")
    plt.close('all')

    # 绘制训练集的结果图
    labels = model.classify(trainData)
    model.plotAndSave(trainData, labels, "./img/trainDataCluster.png")

    # 绘制测试集的结果图
    labels = model.classify(testData)
    model.plotAndSave(testData, labels, "./img/testDataCluster.png")

    # 输出所有样本的聚类标签
    X = dataReader.readDataX("cluster.dat")
    labels = model.classify(X)
    np.savetxt("labels.csv", np.concatenate((X, labels), axis=1), delimiter=',')

    # 将迭代过程中保存下来的中间过程的图片制作成gif, 可以看到GMM模型的变化过程
    imgs = []
    fig = plt.figure()

    for i in range(model.iterations):
        path = "./img/iterations/iteration" + str(i) + ".png"
        img = plt.imread(path)
        img = plt.imshow(img, animated=True)
        imgs.append([img])

    # fig = plt.figure()
    # plt.close(2)
    """plt.figure()函数的作用是新建一张图像，这个图像是一个类似于相框一样的东西，上面可以放图片。这里animation中一张张的图片就是要依次出现在一个figure上。每调用一次函数，就会新建一个figure. 
    最开始写代码的时候，我在for循环添加img之后又调用了一次plt.figure(), 程序结束show的时候就会出现两张figure. 动画展示在figure1上，而figure2上什么都没有，因为对figure2进行任何操作。最后保存的时候，保存的也是最新的figure2. 这里可以认为所有对figure的操作都是在操作最新创建的figure.
    也可以调用plt.close()函数关闭figure. 调用plt.close(2)就是关闭了figure2, 这时候就能够正常显示和保存gif了(gif显示在figure1中)。"""

    ani = animation.ArtistAnimation(fig,
                                    imgs,
                                    interval=500,
                                    blit=True,
                                    repeat_delay=1000)
    ani.save("./img/iterationAnime.gif", writer="pillow")
    # plt.show()


if __name__ == "__main__":
    main()

GMM.py

from io import StringIO
from math import log
from re import S
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from numpy.core.fromnumeric import shape
from scipy.stats import multivariate_normal
from sklearn.cluster import KMeans


class GMM():
    def __init__(self, n_components=1, n=1, eps=1e-2):
        """
        模型本身具有的参数有：
        n_components 聚类的类别数目
        n 特征维度
        mu 模型中心点
        covariance 协方差矩阵
        class_prior 各个模型的先验概率
        eps 停止迭代的阈值。如果前后两次EM计算出的LogLikelihood差值小于eps, 就停止迭代
        LogLikelihood 存储迭代过程中所有的对数似然，绘图用
        iterations 保存最后的迭代次数，绘图用
        """
        self.n_components = n_components
        self.n = n
        self.mu = np.zeros((self.n_components, self.n))
        self.class_prior = np.zeros(self.n_components)
        self.covariance = np.zeros((self.n_components, self.n, self.n))
        self.eps = eps
        self.LogLikelihood = []
        self.iterations = 0

    def train(self, X):
        """
        使用数据X训练模型，没有返回值，更新模型自身的参数
        """
        m = X.shape[0]
        self._paramInitialization(X, m)
        try:
            self._EM(X, m)
        except:
            print(
                "error in EM algorithm. Didn't get a successful model. Need to try again."
            )
            raise
        # self._EM(X, m)

    def classify(self, X):
        """
        使用训练好的模型对X进行聚类，返回聚类得到的标签        
        """
        m = X.shape[0]
        pdfs = np.zeros((m, self.n_components))
        for k in range(self.n_components):
            pdfs[:, k] = self._gaussianPDF(X, self.mu[k], self.covariance[k])
        labels = np.argmax(pdfs, axis=1)
        return labels.reshape(-1, 1)

    def plot(self, X, labels=None):
        """
        将数据集X中的数据点和模型中所有的高斯分布的等高线图绘制在figure上。如果给出了数据点的labels, 那么不同label的数据点有不同的颜色，否则都是灰色。
        """
        if labels is not None:
            plt.scatter(X[:, 0], X[:, 1], s=8, c=labels, alpha=0.5)
        else:
            plt.scatter(X[:, 0], X[:, 1], s=8, c="grey", alpha=0.5)

        for k in range(self.n_components):
            self._plotGaussianContour(self.mu[k], self.covariance[k], X)

    def plotAndSave(self, X, labels=None, savePath='./img/picture.png'):
        """
        调用plot函数，绘制图像，并且按照savePath的路径名保存figure.
        """
        self.plot(X, labels)
        plt.savefig(savePath)
        plt.close('all')

    def plotLogLikelihood(self):
        """
        画出迭代过程中对数似然的变化图象。
        """
        length = len(self.LogLikelihood)
        X = np.linspace(0, length, length)
        Y = self.LogLikelihood
        plt.xlabel('iteration')
        plt.ylabel('LogLikelihood')
        plt.plot(X, Y, color='red')

    def _paramInitialization(self, X, m):
        """
        模型参数初始化。对于mu和covariance的初始化做了一些特殊处理，使得模型中心点尽可能靠近样本中心，协方差矩阵尽可能合理，减小迭代过程中出现奇异矩阵的可能性。
        """

        # """
        # mu是scale过的，保证不超出样本点的范围
        # """
        # self.mu = np.random.random((self.n_components, self.n))
        # minCol = np.min(X, axis=0)
        # maxCol = np.max(X, axis=0)
        # self.mu = minCol + self.mu * (maxCol - minCol)
        """
        使用k-means进行预聚类，获取聚类中心点。
        """
        estimator = KMeans(n_clusters=self.n_components)  #构造聚类器
        estimator.fit(X)  #聚类
        centroids = estimator.cluster_centers_  #获取聚类中心
        self.mu = centroids
        """
        dist是一个(m, n_components)的矩阵，dist[i, k]是第i个样本点到第k个中心点的距离。这里使用了一些神奇的操作使代码简洁。两点间距离的计算公式为 (x-mu)^2 , 这里将公式展开为 x^2 + mu^2 - 2*x*mu。 
        """
        dist = np.tile(
            np.sum(X * X, axis=1).reshape(-1, 1),
            (1, self.n_components)) + np.tile(
                np.sum(self.mu * self.mu, axis=1).reshape(1, -1),
                (m, 1)) - 2 * np.dot(X, self.mu.T)
        labels = np.argmin(dist, axis=1)
        for k in range(self.n_components):
            cluster = X[labels == k, :]
            self.class_prior[k] = cluster.shape[0] / m
            self.covariance[k, :, :] = np.cov(cluster.T)

    def _EM(self, X, m):
        """
        执行EM算法，学习参数
        """
        pdfs = np.zeros((m, self.n_components))  # 产生的概率分布
        gamma = np.zeros((m, self.n_components))  # 归一化之后的概率分布，即后验概率

        # 需要设定算法终止条件
        # 1. likelihood变化非常微小
        # 2. 迭代轮数达到上限
        self.iterations = totalIteration = 100
        preLogLikelihood = LogLikelihood = 0
        for nowIter in range(totalIteration):

            # e-step, 归一化之后得到后验概率
            for k in range(self.n_components):
                try:
                    pdfs[:, k] = self.class_prior[k] * self._gaussianPDF(
                        X, self.mu[k], self.covariance[k])
                except:
                    print(
                        "covariance matrix {} is singular matrix.\n something wrong in producing probabilities. iteration {}"
                        .format(k + 1, nowIter))
                    raise
            gamma = pdfs / np.sum(pdfs, axis=1).reshape(-1, 1)

            # m-step， 依据公式更新参数
            class_prior = np.sum(gamma, axis=0) / np.sum(gamma)
            mu = np.zeros((self.n_components, self.n))
            covariance = np.zeros((self.n_components, self.n, self.n))
            for k in range(self.n_components):
                mu[k] = np.average(X, axis=0, weights=gamma[:, k])
                cov = np.zeros((self.n, self.n))
                for i in range(m):
                    tmp = (X[i] - mu[k]).reshape(-1, 1)
                    cov += gamma[i, k] * np.dot(tmp, tmp.T)
                covariance[k, :, :] = cov / np.sum(gamma[:, k])

            # 更新参数
            self.class_prior = class_prior
            self.mu = mu
            self.covariance = covariance

            savePath = "./img/iterations/iteration" + str(nowIter) + ".png"
            self.plotAndSave(X, savePath=savePath)

            # 计算likelihood, 判断是否能够停止
            try:
                LogLikelihood = self._computeLogLikelihood(X)
            except:
                print(
                    "something wrong in computing LogLikelihood. iteration {}".
                    format(nowIter + 1))
                raise

            self.LogLikelihood.append(LogLikelihood)
            if abs(preLogLikelihood - LogLikelihood) < self.eps:
                self.iterations = nowIter
                break
            preLogLikelihood = LogLikelihood

    def _computeLogLikelihood(self, X):
        """
        计算对数似然
        """
        m = X.shape[0]
        pdfs = np.zeros((m, self.n_components))
        for k in range(self.n_components):
            try:
                pdfs[:, k] = self.class_prior[k] * self._gaussianPDF(
                    X, self.mu[k], self.covariance[k])
            except:
                print("covariance matrix {} is singular matrix".format(k + 1))
                raise
        return np.sum(np.log(np.sum(pdfs, axis=1)))

    def _plotGaussianContour(self, mu, covariance, X):
        """
        在图像上绘制高斯分布的等高线图，同时标注中心点。
        """
        xmin, xmax = np.min(X[:, 0]), np.max(X[:, 0])
        ymin, ymax = np.min(X[:, 1]), np.max(X[:, 0])
        M = 200
        X, Y = np.meshgrid(np.linspace(xmin - 1, xmax + 1, M),
                           np.linspace(ymin - 1, ymax + 1, M))
        d = np.dstack((X, Y))
        Z = multivariate_normal.pdf(d, mu, covariance)
        plt.contour(X, Y, Z, levels=7)

        plt.scatter(mu[0], mu[1], s=100, marker='x', c='red')

    def _gaussianPDF(self, X, mu, covariance):
        """
        按照给定的mu和covariance生成数据集X的概率分布
        """
        # 加一个微小的单位矩阵应对矩阵不可逆的情况
        covdet = np.linalg.det(covariance + np.eye(self.n) * 0.01)  #协方差矩阵的行列式
        covinv = np.linalg.inv(covariance + np.eye(self.n) * 0.01)  #协方差矩阵的逆
        xdiff = X - mu  # (m,n)的矩阵
        #概率密度
        prob = 1.0 / np.power(2 * np.pi, self.n / 2) / np.sqrt(
            np.abs(covdet)) * np.exp(
                -1.0 / 2 * np.diag(np.dot(np.dot(xdiff, covinv), xdiff.T)))
        # xdiff * covinv * xdiff.T  是 (1000, 1000) 的矩阵，对角线上是 (xi - mu)^T * mu^-1 * (xi - mu) , 所以用 np.diag取得对角线上的值
        return prob

DataReader.py

import numpy as np  

class DataReader():
    """
    负责数据的读入和预处理    
    """
    def __init__(self):
        pass

    def readDataX(self, filename):
        """
        从文件中读入只有X的数据，返回数据集X，格式为m*n的二位矩阵，m个样本点，n维特征。
        """
        X = []
        with open(filename, "r", encoding="utf-8") as inputFile:
            for line in inputFile.readlines():
                line = line.strip('\n').split(' ')
                line = [eval(x) for x in line]
                X.append(line)
        X = np.array(X)
        return X

    def sampleSplit(self, data, ratio=0.2):
        """
        将数据随机划分为训练集和测试集，ratio是测试集数据比例
        """
        np.random.shuffle(data)
        n = np.int(data.shape[0]*ratio) 
        testData, trainData = data[:n], data[n:]
        return trainData, testData

经验教训

Ⅰ

main程序中最后绘制动画的时候，一开始怎么都不对，使用plt.show()总是出现两张figure, figure1上面是正确的动画，但是figure2上面是一片空白，执行 ani.save("./img/iterationAnime.gif", *writer*="pillow") 保存的时候也是保存了一片空白，看来是保存了figure2.

经过多次实验，最终找到了原因：

plt.figure()函数的作用是新建一张图像，这个图像是一个类似于相框一样的东西，上面可以放图片。这里animation中一张张的图片就会依次出现在一个figure上。每调用一次plt.figure()，就会新建一个figure。每个figure都有数字编号，从1开始递增。

最开始写代码的时候，我在for循环添加img之后又调用了一次plt.figure(), 程序结束plt.show()的时候就会出现两张figure. 动画展示在figure1上，而figure2上什么都没有，因为并没有对figure2进行任何操作。最后保存的时候，保存的却是最新的figure2. 想必是所有对figure的操作都是在操作最新创建的figure，而保存animation也不例外，当然是默认保存最新的figure.

	imgs = []
    fig = plt.figure()

    for i in range(model.iterations):
        path = "./img/iterations/iteration" + str(i) + ".png"
        img = plt.imread(path)
        img = plt.imshow(img, animated=True)
        imgs.append([img])

    # fig = plt.figure()  就是这一行搞的鬼
    # plt.close(2)  加上这个就没事了
    
    ani = animation.ArtistAnimation(fig,
                                    imgs,
                                    interval=500,
                                    blit=True,
                                    repeat_delay=1000)
    ani.save("./img/iterationAnime.gif", writer="pillow")

调用plt.close()函数可以关闭figure. 调用plt.close(2)就是关闭了figure2, 这时候就能够正常显示和保存gif了(gif绘制在figure1中)。

Ⅱ

编写程序的初期总是出现 singular matrix 的问题。这是在调用 scipy.stats.multivariate_normal.pdf 的时候出现的，因为需要生成数据点的概率分布，而其中会有求协方差矩阵的逆矩阵的操作。这个函数不能处理矩阵不可逆（奇异矩阵）的情况，所以就会报错。而如果GMM一开始的聚类中心点是随机生成的话，就很容易迭代着迭代着就有一个协方差矩阵不可逆了。

后面是参考了一个能够处理逆矩阵的高斯分布概率生成函数。

        # 加一个微小的单位矩阵应对矩阵不可逆的情况
        covdet = np.linalg.det(covariance + np.eye(self.n) * 0.01)  #协方差矩阵的行列式
        covinv = np.linalg.inv(covariance + np.eye(self.n) * 0.01)  #协方差矩阵的逆

这样做确实是有用的，起码程序不会总是抛出异常。虽然最后可能会出现这样的情况：

其中有一个分布中间就缩没了，协方差矩阵明显成了奇异矩阵了，但是中间运行过程不会报错， 只会提示

后来试着把单位矩阵的系数改小，最大程度上减小影响，发现改到1e-4的时候出现了singular matrix 的情况。

后来改进了一下求中心点的初始化，用了一个KMeans先求一下中心点。这样即使是使用原来的概率分布函数也基本不会出现有singular matrix的情况了。

		"""
        使用k-means进行预聚类，获取聚类中心点。
        """
        estimator = KMeans(n_clusters=self.n_components)  #构造聚类器
        estimator.fit(X)  #聚类
        centroids = estimator.cluster_centers_  #获取聚类中心
        self.mu = centroids

Ⅲ

对程序中一些比较难理解的代码做一下解释

_gaussianPDF

prob = 1.0 / np.power(2 * np.pi, self.n / 2) / np.sqrt(
            np.abs(covdet)) * np.exp(
                -1.0 / 2 * np.diag(np.dot(np.dot(xdiff, covinv), xdiff.T)))

首先，多元高斯分布的概率密度函数形式为

这里传进来的参数是一整个数据集$X$ , 不想写循环，就直接用$X-\mu$去和${\Sigma }^{-1}$相乘了。

np.dot(xdiff, covinv), xdiff.T)这个结果是一个$m\times m$的矩阵($m$是数据点的个数)，用$A$表示，那么$A_{ij}=({\mathbf{x}}_i-\mu {)}^T{\Sigma }^{-1}({\mathbf{x}}_j-\mu )$，那么对角线上的$m$个数据就是我们想要的数据$({\mathbf{x}}_i-\mu {)}^T{\Sigma }^{-1}({\mathbf{x}}_i-\mu )$。之后的操作就按照公式来就行了。

_paramInitialization

dist = np.tile(
            np.sum(X * X, axis=1).reshape(-1, 1),
            (1, self.n_components)) + np.tile(
                np.sum(self.mu * self.mu, axis=1).reshape(1, -1),
                (m, 1)) - 2 * np.dot(X, self.mu.T)

这个算完之后是一个$m\ast n\_components$的矩阵，记为$D$吧。那么$D_{ij}$表示的就是$x_i$到${\mu }_j$的距离。之后我们要从每一行中找最小值，$x_i$离哪个中心点近，就是哪一类。统计个数，确定一下先验概率。说一下里面的各个操作。

np.sum(X * X, axis=1)这个算完之后得到一个$m$维的向量，第$i$个数据是$x_i^2$ 。 np.tile是做一下扩展，np.tile(np.sum(X * X, axis=1).reshape(-1, 1), (1, self.n_components))首先把那个$m$维的向量重构成$m\ast 1$的矩阵，然后按照(1, self.n_components)的模式进行复制，在第0维上复制1次，在第1维上复制self.n_components次。这样就得到了一个$m\ast n\_components$的矩阵。

np.tile(np.sum(self.mu * self.mu, axis=1).reshape(1, -1),(m, 1))的操作模式基本一样，只不过计算的是${\mu }_k^2$，而且扩展的方向不一样。

2 * np.dot(X, self.mu.T)计算的东西很明显。

最后合起来，得到矩阵$D$。$D_{ik}=x_i^2+{\mu }_k^2-2x_i{\mu }_k=(x_i-{\mu }_k{)}^2$

Ⅳ

关于公式部分的编程是很重要的部分，要经常拿出来复习一下，尤其是怎么把代码写的简洁高效。

def _EM(self, X, m):
for k in range(self.n_components):
                try:
                    pdfs[:, k] = self.class_prior[k] * self._gaussianPDF(
                        X, self.mu[k], self.covariance[k])
                except:
                    print(
                        "covariance matrix {} is singular matrix.\n something wrong in producing probabilities. iteration {}"
                        .format(k + 1, nowIter))
                    raise
            gamma = pdfs / np.sum(pdfs, axis=1).reshape(-1, 1)

            # m-step， 依据公式更新参数
            class_prior = np.sum(gamma, axis=0) / np.sum(gamma)
            mu = np.zeros((self.n_components, self.n))
            covariance = np.zeros((self.n_components, self.n, self.n))
            for k in range(self.n_components):
                mu[k] = np.average(X, axis=0, weights=gamma[:, k])
                cov = np.zeros((self.n, self.n))
                for i in range(m):
                    tmp = (X[i] - mu[k]).reshape(-1, 1)
                    cov += gamma[i, k] * np.dot(tmp, tmp.T)
                covariance[k, :, :] = cov / np.sum(gamma[:, k])
                
                
                
                
                
def _computeLogLikelihood(self, X):
        """
        计算对数似然
        """
        m = X.shape[0]
        pdfs = np.zeros((m, self.n_components))
        for k in range(self.n_components):
            try:
                pdfs[:, k] = self.class_prior[k] * self._gaussianPDF(
                    X, self.mu[k], self.covariance[k])
            except:
                print("covariance matrix {} is singular matrix".format(k + 1))
                raise
        return np.sum(np.log(np.sum(pdfs, axis=1)))
    
    
    
    
    
    
    
    
    
def _gaussianPDF(self, X, mu, covariance):
        """
        按照给定的mu和covariance生成数据集X的概率分布
        """
        # 加一个微小的单位矩阵应对矩阵不可逆的情况
        covdet = np.linalg.det(covariance + np.eye(self.n) * 1e-3)  #协方差矩阵的行列式
        covinv = np.linalg.inv(covariance + np.eye(self.n) * 1e-3)  #协方差矩阵的逆
        # covdet = np.linalg.det(covariance)  
        # covinv = np.linalg.inv(covariance)  
        xdiff = X - mu  # (m,n)的矩阵
        #概率密度
        prob = 1.0 / np.power(2 * np.pi, self.n / 2) / np.sqrt(
            np.abs(covdet)) * np.exp(
                -1.0 / 2 * np.diag(np.dot(np.dot(xdiff, covinv), xdiff.T)))
        # xdiff * covinv * xdiff.T  是 (1000, 1000) 的矩阵，对角线上是 (xi - mu)^T * mu^-1 * (xi - mu) , 所以用 np.diag取得对角线上的值
        return prob

附录:理论基础

$${\pi }_k=\frac{N_k}{N}$$