聚类——MeanShift算法以及Python实现

均值漂移算法(MeanShift)是一种旨在发现团(blobs)的聚类算法

核心思想

寻找核密度极值点并作为簇的质心，然后根据最近邻原则将样本点赋予质心

算法简介

核密度估计

根据样本分布估计在样本空间的每一点的密度。估计某点的密度时，核密度估计方法会考虑该点邻近区域的样本点的影响，邻近区域大小由带宽h决定，该参数对最终密度估计的影响非常大。通常采用高斯核：$$N(x)=\frac{1}{\sqrt{2\pi }h}{e}^{\frac{x^2}{2h^2}}$$

均值漂移

算法初始化一个质心(向量表示)，每一步迭代都会朝着当前质心领域内密度极值方向漂移，方向就是密度上升最大的方向，即梯度方向（参考资料）。求导即得漂移向量的终点就是下一个质心：$$x_k^{i+1}=\frac{{\sum }_{x_j\in N}K(x_j-x_k^i)x_j}{{\sum }_{x_j\in N}K(x_j-x_k^i)}$$N是当前质心的领域内样本集合。这个公式可以从另一方面理解：当前质心的邻域内样本点以核密度为权重的均值就是更新后的质心。

算法流程

• Input: 高斯核带宽，bin_seeding(是否对数据粗粒化)， min_fre(可以作为起始质心的样本点领域内的最少样本数)，阈值epsilon(漂移最小长度)
• Output: 样本簇标签
• Step1: 获取可以作为起始质心的样本点
• Step2: 对每个起始质心进行漂移，漂移终止条件就是漂移距离小于epsilon。若漂移结束最终的质心与已存在的质心距离小于带宽则合并
• Step3: 分类。将样本点归属到距离最近的质心中

代码

"""
meanshift聚类算法
核心思想：
寻找核密度极值点并作为簇的质心，然后根据最近邻原则将样本点赋予质心
"""
from collections import defaultdict
import numpy as np
import math


class MeanShift:
    def __init__(self, band_width=2.0, min_fre=3, epsilon=None, bin_seeding=False, bin_size=None):
        self.epsilon = epsilon if epsilon else 1e-3 * band_width
        self.bin_size = bin_size if bin_size else self.band_width
        self.band_width = band_width
        self.min_fre = min_fre  # 可以作为起始质心的球体内最少的样本数目
        self.bin_seeding = bin_seeding
        self.radius2 = self.band_width ** 2  # 高维球体半径的平方

        self.N = None
        self.labels = None
        self.centers = []

    def init_param(self, data):
        # 初始化参数
        self.N = data.shape[0]
        self.labels = -1 * np.ones(self.N)
        return

    def get_seeds(self, data):
        # 获取可以作为起始质心的点（seed）
        if not self.bin_seeding:
            return data
        seed_list = []
        seeds_fre = defaultdict(int)
        for sample in data:
            seed = tuple(np.round(sample / self.bin_size))  # 将数据粗粒化，以防止非常近的样本点都作为起始质心
            seeds_fre[seed] += 1
        for seed, fre in seeds_fre.items():
            if fre >= self.min_fre:
                seed_list.append(np.array(seed))
        if not seed_list:
            raise ValueError('the bin size and min_fre are not proper')
        if len(seed_list) == data.shape[0]:
            return data
        return np.array(seed_list) * self.bin_size

    def euclidean_dis2(self, center, sample):
        # 计算均值点到每个样本点的欧式距离（平方）
        delta = center - sample
        return delta @ delta

    def gaussian_kel(self, dis2):
        # 计算高斯核
        return 1.0 / self.band_width * (2 * math.pi) ** (-1.0 / 2) * math.exp(-dis2 / (2 * self.band_width ** 2))

    def shift_center(self, current_center, data):
        # 计算下一个漂移的坐标
        denominator = 0  # 分母
        numerator = np.zeros_like(current_center)  # 分子, 一维数组形式
        for sample in data:
            dis2 = self.euclidean_dis2(current_center, sample)
            if dis2 <= self.radius2:
                d = self.gaussian_kel(dis2)
                denominator += d
                numerator += d * sample
        if denominator > 0:
            return numerator / denominator
        else:
            return None

    def classify(self, data):
        # 根据最近邻将数据分类到最近的簇中
        center_arr = np.array(self.centers)
        for i in range(self.N):
            delta = center_arr - data[i]
            dis2 = np.sum(delta * delta, axis=1)
            self.labels[i] = np.argmin(dis2)
        return

    def fit(self, data):
        # 训练主函数
        self.init_param(data)
        seed_list = self.get_seeds(data)
        for seed in seed_list:
            bad_seed = False
            current_center = seed
            # 进行一次独立的均值漂移
            while True:
                next_center = self.shift_center(current_center, data)
                if next_center is None:
                    bad_seed = True
                    break
                delta_dis = np.linalg.norm(next_center - current_center, 2)
                if delta_dis < self.epsilon:
                    break
                current_center = next_center
            if not bad_seed:
                # 若该次漂移结束后，最终的质心与已存在的质心距离小于带宽，则合并
                for i in range(len(self.centers)):
                    if np.linalg.norm(current_center - self.centers[i], 2) < self.band_width:
                        break
                else:
                    self.centers.append(current_center)
        self.classify(data)
        return


if __name__ == '__main__':
    from sklearn.datasets import make_blobs

    data, label = make_blobs(n_samples=500, centers=5, cluster_std=1.2, random_state=7)
    MS = MeanShift(band_width=3, min_fre=3, bin_size=4, bin_seeding=True)
    MS.fit(data)
    labels = MS.labels
    print(MS.centers, np.unique(labels))
    import matplotlib.pyplot as plt
    from itertools import cycle


    def visualize(data, labels):
        color = 'bgrymk'
        unique_label = np.unique(labels)
        for col, label in zip(cycle(color), unique_label):
            partial_data = data[np.where(labels == label)]
            plt.scatter(partial_data[:, 0], partial_data[:, 1], color=col)
        plt.show()
        return


    visualize(data, labels)

我的GitHub
注：如有不当之处，请指正。