本文最后更新于：2021年10月28日晚上

K-Means&EM聚类

理论介绍

什么是聚类分析：

聚类分析是在数据中发现数据对象之间的关系，将数据进行分组，组内的相似性越大，组间的差别越大，则聚类效果越好。聚类是一种非监督学习。

K-Means算法：

优点：易于实现

缺点：容易陷于局部最小值

算法的伪代码可表示如下：

选择K个点作为初始质心
repeat
将每个点指派到最近的质心，形成K个簇
重新计算每个簇的质心
until 簇不发生变化或达到最大迭代次数

这里的重新计算每个簇的质心，如何计算的是根据目标函数得来的，因此在开始时我们要考虑距离度量和目标函数。

考虑欧几里得距离的数据，使用误差平方和（Sum of the Squared Error,SSE）作为聚类的目标函数，两次运行K均值产生的两个不同的簇集，我们更喜欢SSE最小的那个。

SSE=\Sigma^K_{i=1}\Sigma_{x\in C_i}dist(c_i, x)^2

$K$ 表示 $K$ 个聚类中心， $c_i$ 表示第几个中心， $dist$ 表示的是欧几里得距离。
这里有一个问题就是为什么，我们更新质心是让所有的点的平均值，这里就是SSE所决定的。

这里我们对第 $k$ 个质心 $c_k$ 求解，最小化公式(1)，也就是令导数=0.求解 $c_k$ ：

\frac{\partial}{\partial c_k}SSE = \frac{\partial}{\partial c_k}\Sigma^K_{i=1}\Sigma_{x\in C_i}(c_i - x)^2\\ =\Sigma^K_{i=1}\Sigma_{x\in C_i}\frac{\partial}{\partial c_k}(c_i - x)^2\\ =\Sigma_{x\in C_i}2(c_i - x)=0\\ \to m_kc_k=\Sigma_{x\in C_i}x_k=\frac{1}{m_k}\Sigma_{x\in C_i}x_k

所以，簇的最小化SSE的最佳质心是簇中各点中心的值。

EM聚类

最大期望算法（Expectation-maximization algorithm，又译期望最大化算法）在统计中被用于寻找，依赖于不可观察的隐性变量的概率模型中，参数的最大似然估计。

在统计计算中，最大期望（EM）算法是在概率（probabilistic）模型中寻找参数最大似然估计或者最大后验估计的算法，其中概率模型依赖于无法观测的隐藏变量（Latent Variable）。最大期望经常用在机器学习和计算机视觉的数据聚类（Data Clustering）领域。最大期望算法经过两个步骤交替进行计算，第一步是计算期望（E），利用对隐藏变量的现有估计值，计算其最大似然估计值；第二步是最大化（M），最大化在 E 步上求得的最大似然值来计算参数的值。M 步上找到的参数估计值被用于下一个 E 步计算中，这个过程不断交替进行。

下面介绍基于GMM(高斯混合模型)的EM聚类算法步骤：

初始化K个N维高斯分布参数 $\mu、\Sigma$ 和 $\pi$

repeat:

E步:

$\gamma(x_n,k) = \frac{\pi_k N(x_n|\mu_k,\Sigma_k)}{\Sigma^K_{j=1}\pi_j N(x_n|\mu_j,\Sigma_j)}, k = 1, 2\dots K$

$\pi_{k_{new}} = \Sigma^N_{j=1}\gamma(x_j,k) / N$

M步:

$\mu_{k_{new}}=\frac{1}{N_k}\Sigma_{j=1}^N \gamma(x_j,k)x_j$

$\Sigma_{k_{new}} = \frac{1}{N_k}\Sigma_{j=1}^N \gamma(x_j,k)(x_j-\mu_{k_{new}})(x_j-\mu_{k_{new}})^T$

更新参数

until 簇不发生变化或达到最大迭代次数

其中， $\pi_k$ 为第 $k$ 类中元素占全部元素的比例， $K$ 为要分类的种类数， $\gamma(x_n,k)$ 为 $x_n \in k类$ 的概率。

实验结果

实验中我们发现，基本上每次迭代次数达到60左右EM算法开始收敛，在下图的实验中，聚类成功率大概在54%以上

条件	结果
EM算法陷入局部最小值
EM算法陷入局部最小值（2）
EM算法成功聚类

对于K-Means，我们惊人地发现，它在这个同样的数据集下，迭代10次即可收敛，且聚类成功率为100%

条件	结果
K-Means成功聚类（2）
K-Means成功聚类（2）

在多次试验后我们发现，K-Means最终都收敛与同一分布情况。

随后我们还进行了多次实验，普遍发现EM算法收敛慢，容易陷入局部最小值，但是也往往能够应对各种不同的情况，而K-Means收敛快，往往只会收敛于几个确定的分布。我们认为这是因为K-Means只能处理很多点聚集在一起的情况，一但聚类的形状与混合情况复杂，K-Means就会完全失去处理能力。而基于GMM的EM由于 $\Sigma$ 参数的存在，可以处理聚类形状和分布更加复杂的情况。

代码实现

本代码在去耦和的基础上，提升了其可视化性能，并针对EM算法实现了由点的属于三个类的概率到RGB三种颜色的映射，使得我们能够通过点的颜色窥察该点身上的概率分布情况

首先是K-Means代码实现:

K-Means对数据初始化最优的方法应该是K-Means++方法，但是由于这里数据较少，直接随机抽取也能取得较好的结果。

class KMeansCluster:
    def __init__(self, data: np.ndarray):
        self.dataSet = data
        self.K = 0
        self.centerPos = None

    def addData(self, data: np.ndarray):
        self.dataSet = np.append(self.dataSet, data, axis=0)

    def reInit(self, K=2, init='random'):
        indices = np.random.choice(self.dataSet.shape[0], K, replace=False)
        self.centerPos = self.dataSet[indices]
        self.K = K

    def train(self, n=1000, plt=None):
        for i in range(n):
            # 对每个点，找到距离它最近的中心归类
            classes = [[] for k in range(self.K)]
            for data in self.dataSet:
                distanceWithCenter = np.linalg.norm(data - self.centerPos, ord=2, axis=1, keepdims=False)
                nearCenterI = np.argmin(distanceWithCenter)
                classes[nearCenterI].append(data)
            classes = [np.array(c) for c in classes]
            classes = np.array(classes)

            # 对每个类，找出其新中心
            self.centerPos = [c.mean(axis=0) for c in classes]

            print(i)

            if plt is not None:
                plt.subplot(1, 2, 2, frameon=False)
                plt.cla()
                plt.scatter(classes[0].T[0], classes[0].T[1], marker='o', c='r', alpha=0.7, label='red')
                plt.scatter(classes[1].T[0], classes[1].T[1], marker='o', c='b', alpha=0.7, label='blue')
                plt.scatter(classes[2].T[0], classes[2].T[1], marker='o', c='g', alpha=0.7, label='green')
                for center in self.centerPos:
                    plt.scatter(center[0], center[1], marker='o', c='black', alpha=0.7, label='center')
                plt.legend()
                plt.pause(0.01)

EM代码实现：

EM算法中包含大量求和等运算，我们尽可能地将其转化为基于NumPy的矩阵运算以提高其运行速度

import numpy as np
from scipy.stats import multivariate_normal


class EM_Cluster:
    def __init__(self, data: np.ndarray, K):
        self.dataSet = data
        self.K = K
        self.rank = self.dataSet.shape[1]
        self.PiK = np.ones(K) / 3.0
        self.centerNormal = [(np.random.uniform(-3, 3, size=self.rank), np.identity(self.rank) * 3) for i in range(K)]

    def addData(self, data: np.ndarray):
        self.dataSet = np.append(self.dataSet, data, axis=0)

    def train(self, n=1000, plt=None):
        N = self.dataSet.shape[0]
        for i in range(n):
            # E步，根据当前分布计算每个点属于各个分布的概率
            gammas = []
            for data in self.dataSet:
                sumGamma = 0
                gamma = []
                # 对每个类别，计算这个点的pdf
                m = 0
                for mean, sigma in self.centerNormal:
                    g = self.PiK[m] * multivariate_normal.pdf(data, mean=mean, cov=sigma)
                    sumGamma += g
                    gamma.append(g)
                    m += 1
                # 归一化概率
                gamma = np.array(gamma)
                gamma /= sumGamma
                gammas.append(gamma)
            gammas = np.array(gammas)
            self.PiK = gammas.sum(axis=0) / gammas.sum()

            # M步根据每个点的分布情况更新参数
            mu_new = gammas.T.dot(self.dataSet) / np.array([self.PiK, self.PiK]).T / N
            sigmas_new = []
            for j in range(self.K):
                m = 0
                sigma_new = np.array([[0.0, 0.0], [0.0, 0.0]])
                for data in self.dataSet:
                    sigma_new += (data - mu_new[j])[None].T.dot((data - mu_new[j])[None]) * gammas[m][j]
                    m += 1
                sigmas_new.append(sigma_new)
            sigmas_new /= self.PiK[:, None, None] * N

            self.centerNormal = [(mu_new[m], sigmas_new[m]) for m in range(self.K)]

            # print(mu_new)
            # print(sigmas_new)

            print('\r trun:%d/%d' % (i + 1, n), end='')

            if plt is not None:
                m = 0
                plt.subplot(1, 2, 2, frameon=False)
                plt.cla()
                for data in self.dataSet:
                    plt.scatter(data[0], data[1], marker='o', color=gammas[m], alpha=0.7)
                    m += 1
                for k in range(self.K):
                    plt.scatter(self.centerNormal[k][0][0], self.centerNormal[k][0][1], marker='o', c='black', alpha=0.7)
                plt.pause(0.01)

测试程序部分：

from GaussianDistribution import *
from KMeansCluster import *
from EM_Cluster import *
import time

if __name__ == '__main__':
    plt.figure("cluster")
    plt.ion()
    np.random.seed(1553215005)
    redData, _ = Gaussian_Distribution(N=2, M=30, m=[0, 0], sigma=[[1, 0], [0, 1]])
    redLabel = np.ones(30) * -1
    blueData, _ = Gaussian_Distribution(N=2, M=20, m=[1, 2], sigma=[[2, 0], [0, 1]])
    blueLabel = np.ones(20) * -1
    greenData, _ = Gaussian_Distribution(N=2, M=20, m=[2, 0], sigma=[[1, 0.3], [0.3, 1]])
    greenLabel = np.ones(20) * -1

    plt.subplot(1, 2, 1, frameon=False)
    plt.scatter(redData.T[0], redData.T[1], marker='o', c='r', alpha=0.7, label='red')
    plt.scatter(blueData.T[0], blueData.T[1], marker='o', c='b', alpha=0.7, label='blue')
    plt.scatter(greenData.T[0], greenData.T[1], marker='o', c='g', alpha=0.7, label='green')
    plt.legend()
    plt.pause(2)

    np.random.seed(int(time.time()))
    kCluster = KMeansCluster(redData)
    kCluster.addData(blueData)
    kCluster.addData(greenData)
    kCluster.reInit(K=3)
    kCluster.train(plt=plt)

    # emCluster = EM_Cluster(redData, 3)
    # emCluster.addData(blueData)
    # emCluster.addData(greenData)
    # emCluster.train(plt=plt)

高斯分布数据生成部分

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
from mpl_toolkits.mplot3d import Axes3D


def Gaussian_Distribution(N=2, M=1000, m: list = None, sigma: list = None):
    """
    :param sigma: 样本方差
    :param M: 样本数
    :param m: 样本均值
    :param N: 维度

    :return data:  shape(M, N), M 个 N 维服从高斯分布的样本
    :return Gaussian:  高斯分布概率密度函数
    """
    mean = np.array(m)  # or np.zeros(N)  # 均值矩阵，每个维度的均值都为 m
    cov = np.array(sigma)  # or np.eye(N)  # 协方差矩阵，每个维度的方差都为 sigma

    # 产生 N 维高斯分布数据
    data = np.random.multivariate_normal(mean, cov, M,)
    # N 维数据高斯分布概率密度函数
    Gaussian = multivariate_normal(mean=mean, cov=cov)

    return data, Gaussian

机器学习

理论与代码

本博客所有文章除特别声明外，均采用 CC BY-SA 4.0 协议，转载请注明出处！

扩增你的KiCAD库上一篇

Hexo博客搭建后的细节优化下一篇