小瓜讲机器学习——聚类算法（二）Mean Shift算法原理及Python代码实现

2. Mean Shift算法原理

2.1 核函数估计

统计推断是指获得样本数据（总体中的采样点）后通过分析，以推断总体特征，比如通过样本均值、方差推断总体的均值和方差是多少。
其中概率分布是统计推断中最重要的一个内容，通过样本来推断估计总体的概率分布有两种方法：1.假设一个有待定参数的概率分布形式，通过求解待定参数来确定总体分布称为参数估计，2.不使用假设，直接根据样本数据中的信息估计出总体的概率分布称为非参数估计。

核函数估计是一种非参数估计方法。

在$d$维特征空间上概率密度函数$f(X)$的核函数估计为
$$f(X)=\frac{{\sum }_{i=1}^{m}K(\frac{{\mathbf{X}}^{(\mathbf{i})}-\mathbf{X}}{\mathbf{H}})}{m\times h_1{h}_2...h_d}$$
其中$h_i$是第$i$维度上的带宽。如果是同带宽的情况下，核函数估计就变成
$$f(X)=\frac{{\sum }_{i=1}^{m}K(\frac{x^{(i)}-x}{h})}{m\times h^d}$$

• 注：更多非参数估计内容可以参见王星著的《非参数估计》

2.2 Mean Shift向量

对于$d$维特征的样本点集合 T = { X ( 1 ) , . . . , X ( m ) } T=\{X^{(1)}, ...,X^{(m)}\} T={X(1),...,X(m)}，其中$S_h(X^{(j)})$是一个超球体，定义如下
S h ( X ( j ) ) = { Y ∣ ∣ ∣ Y − X ∣ ∣ 2 2 ≤ h 2 } S_h(X^{(j)})=\{Y|\quad ||Y-X||_2^2\le h^2\} Sh​(X(j))={Y∣∣∣Y−X∣∣22​≤h2}
例如下图

$X^{(j)}$为蓝色点，所有黑边灰色点都是超球体内的其他样本点，红边白点为所有样本点的均值点，那么定义从$X^{(j)}$点到均值点的Mean Shift向量如下所示
$$M_h(X^{(j)})=\frac{1}{k}\sum _{X^{(i)\in S_h}}(X^{(i)}-X^{(j)})$$
其中$k$为所有超球体内的样本点。

其实在超球体内的样本点只是总体的抽样点，对于单峰值的概率密度函数来说，均值点处总是相对概率高的点。Mean Shift算法的本质是就是在每一步中将样本点$X^{(j)}$移动到概率高的均值点，最终形成聚类中心。具体的数学原理见下节。

• 这里也可见超球体的半径要取得合适。

2.3 Mean Shift算法原理（由概率密度函数导出）

对于某一$d$维聚类问题，其实存在对应$d$维特征空间上概率密度函数$f(X=(x_1,...,x_d))$，$m$个训练样本就是$d$维特征空间内的采样点。那么概率密度函数的核函数估计为
$$f(X)=\frac{{\sum }_{i=1}^{m}K(\frac{{\mathbf{X}}^{(\mathbf{i})}-\mathbf{X}}{\mathbf{H}})}{m\times h_1{h}_2...h_d}=\frac{{\sum }_{i=1}^{m}K(||\frac{x_j^{(i)}-x_j}{h_j}|{|}_2^2)}{m{h}_1{h}_2...h_d},j=1,...,l$$
其中$H=(h_1,...,h_d)$，当采用等带宽时上式就变成
$$f(X)=\frac{{\sum }_{i=1}^{m}K(||\frac{x_j^{(i)}-x_j}{h}|{|}_2^2)}{m{h}^d},j=1,...,l$$

真实的聚类中心（待求点）附近必然是样本点最密集的，即在待求解的聚类中心的邻域中是采样点最集中的区域，对应着在该邻域内也是概率密度函数$f(X=(x_1,...,x_l))$取到极大值的区域。那么有
$$\begin{array}{rl}\frac{\partial f(X)}{\partial x_j}& =\frac{{\sum }_{i=1}^m\partial K(||\frac{x_j^{(i)}-x_j}{h_j}|{|}_2^2)/\partial x_j}{m{h}_1{h}_2...h_d}\\ & =\frac{{\sum }_{i=1}^m\partial K(u)/\partial u\bullet \partial u/\partial x_j}{m{h}_1{h}_2...h_d}\Rightarrow u=||\frac{x_j^{(i)}-x_j}{h_j}|{|}_2^{2}j=1,...,l\\ & =\frac{1}{m{h}_1{h}_2...h_d}\sum _{i=1}^m\frac{\partial K(u)}{\partial u}\bullet 2\frac{x_j^{(i)}-x_j}{h_j}\bullet (-\frac{1}{h_j})\\ & =\frac{1}{m{h}_1{h}_2...h_d}\frac{2}{h_j^2}\sum _{i=1}^m[-\frac{\partial K(u)}{\partial u}]\bullet (x_j^{(i)}-x_j)\\ & =\frac{2}{h_j^2}\frac{{\sum }_{i=1}^m[-\frac{\partial K(u)}{\partial u}]}{m{h}_1{h}_2...h_d}\bullet \frac{{\sum }_{i=1}^m[-\frac{\partial K(u)}{\partial u}]\bullet (x_j^{(i)}-x_j)}{{\sum }_{i=1}^m[-\frac{\partial K(u)}{\partial u}]}\end{array}$$
其中由核函数的性质可知，下式也是核函数。
$$G(u)=-\frac{\partial K(u)}{\partial u}$$
那么有
$$\begin{array}{rl}\frac{\partial f(X)}{\partial x_j}& =\frac{2}{h_j^2}\frac{{\sum }_{i=1}^{m}G(||\frac{x_j^{(i)}-x_j}{h_j}|{|}_2^2)}{m{h}_1{h}_2...h_d}\bullet \frac{{\sum }_{i=1}^{m}G(||\frac{x_j^{(i)}-x_j}{h_j}|{|}_2^2)\bullet (x_j^{(i)}-x_j)}{{\sum }_{i=1}^{m}G(||\frac{x_j^{(i)}-x_j}{h_j}|{|}_2^2)}j=1,...,l\\ & =\frac{2}{h_j^2}\bullet f_g(X)\bullet M_h(x_j)\end{array}$$
概率密度函数$f_g(X)>0$，所以$\frac{\partial f(X)}{\partial x_j}$与$M_h(x_j)$同号，由$\frac{\partial f(X)}{\partial x_j}=0\Rightarrow M_h(x_j)=0$，这意味着，Mean Shift向量总是指向$f(X)$增加的方向，当向量为零时，概率密度函数取到极大值，而这就是聚类中心。

2.4 Mean Shift算法过程

Mean Shift向量更新公式如下
$$\begin{array}{rl}{M}_h(x_j)& =\frac{{\sum }_{i=1}^{m}G(||\frac{x_j^{(i)}-x_j}{h_j}|{|}_2^2)\bullet (x_j^{(i)}-x_j)}{{\sum }_{i=1}^{m}G(||\frac{x_j^{(i)}-x_j}{h_j}|{|}_2^2)}j=1,...,l\\ & =\frac{{\sum }_{i=1}^{m}G(||\frac{x_j^{(i)}-x_j}{h_j}|{|}_2^2)x_j^{(i)}}{{\sum }_{i=1}^{m}G(||\frac{x_j^{(i)}-x_j}{h_j}|{|}_2^2)}-x_j\\ & =m_h(x_j)-x_j=0\end{array}$$
写成矩阵形式就是$M_h(X^i)=m_h(X^i)-X^i=0$。
Mean Shift聚类算法的流程如下：
输入：
　　　训练样本为 T = { X ( 1 ) , X ( 2 ) , . . . , X ( m ) } T=\{X^{(1)},X^{(2)},...,X^{(m)}\} T={X(1),X(2),...,X(m)}，其中$X^{(i)}=(x_1^{(i)},x_2^{(i)},...,x_l^{(i)})$，选择核函数$K$，选择带宽参数$h_i$，迭代终止参数$\xi$。
过程：
　　　1.按顺序选择其中一个样本点$X^i$，计算与其他样本点的距离列阵$D$
　　　2.以距离列阵$D$，计算核函数$G(||\frac{x_j^{(i)}-x_j}{h_j}|{|}_2^2),j=1,...,l$
　　　3.计算其$m_h(X^i)$
　　　4.判断$||m_h(X^i)-X^i||<\xi$，若是，则$X^i$就是聚类中心，退出循环，若否，则做下一步
　　　5.将该样本点移动到$X_{new}^i=m_h(X^i)$
　　　6.重复上述动作，直到所有聚类中心确定。
输出：
　　　聚类中心。

2.5 Mean Shift算法Python代码实现

首先是建立距离列阵计算函数：

def distance(loop, data):
    m = int(data.shape[0])
    dis = np.zeros((m,1))

    delta = data-loop
    for loopi in range(m):
        dis[loopi] = np.dot(delta[loopi], delta[loopi])
    dis = np.sqrt(dis)
    return dis

Means Shift算法主要部分如下：

def meanshift(data, h, delta):
    m = data.shape[0]
    max_dis = 1
    min_dis = delta
    dealed = np.zeros((m,1))
    ones = np.ones((m,1))
    data_new = data.copy()
    # do when ||m_h(X^i)-X^i||<delta
    while max_dis>delta:
        max_dis = 0
        for loopi in range(int(data.shape[0])):
        # select undealed node
            if(dealed[loopi]==0):
            # dis is distance vector
            # guass is kernel value of all node in X^i 
                dis = distance(data_new[loopi], data)
                guass = 1/(math.sqrt(2*math.pi)*h)*np.exp(-dis**2/h**2)
                sum1 = guass.T.dot(data)
                sum2 = guass.T.dot(ones)
                mh = sum1/sum2
                vector = mh-data_new[loopi]
                # get max value
                if(max_dis<np.linalg.norm(vector)):
                    max_dis = np.linalg.norm(vector)
                # if ||vector||<min_dis, do nothing
                if(np.linalg.norm(vector)<min_dis):
                    dealed[loopi] = 1
                else:
                    data_new[loopi] = mh

    return data_new

返回的data_new是将所有的样本点都移动到聚类中心，然后需要将原来的data标记上类别标签

def group(data_new, delta):
    m = data_new.shape[0]
    label = np.zeros((m,1))
    number = 1
    
    for loopi in range(m):
        if(label[loopi]==0):
            deltadata = data_new - data_new[loopi]
            dis = (deltadata * deltadata).sum(axis=1)
            label[dis<delta] = number
            number += 1

    return label

结果如下

附录：Python源代码

python源代码和数据集如下

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import math as math

def meanshift(data, h, delta):
    m = data.shape[0]
    max_dis = 1
    min_dis = delta
    dealed = np.zeros((m,1))
    ones = np.ones((m,1))
    data_new = data.copy()
    while max_dis>delta:
        max_dis = 0
        for loopi in range(int(data.shape[0])):
            if(dealed[loopi]==0):
                dis = distance(data_new[loopi], data)
                guass = 1/(math.sqrt(2*math.pi)*h)*np.exp(-dis**2/h**2)
                sum1 = guass.T.dot(data)
                sum2 = guass.T.dot(ones)
                mh = sum1/sum2

                vector = mh-data_new[loopi]
                if(max_dis<np.linalg.norm(vector)):
                    max_dis = np.linalg.norm(vector)
                if(np.linalg.norm(vector)<min_dis):
                    dealed[loopi] = 1
                else:
                    data_new[loopi] = mh

    return data_new

def distance(loop, data):
    m = int(data.shape[0])
    dis = np.zeros((m,1))

    delta = data-loop
    for loopi in range(m):
        dis[loopi] = np.dot(delta[loopi], delta[loopi])
    #dis = delta.dot(delta.T)
    dis = np.sqrt(dis)
    return dis

def group(data_new, delta):
    m = data_new.shape[0]
    label = np.zeros((m,1))
    number = 1
    for loopi in range(m):
        if(label[loopi]==0):
            deltadata = data_new - data_new[loopi]
            dis = (deltadata * deltadata).sum(axis=1)
            label[dis<delta] = number
            number += 1

    return label


if __name__=='__main__':
    data = []
    with open(r'x\data.txt') as f:
        for loop in f.readlines():
            line = loop.strip().split('\t')
            temp = [float(line[0]), float(line[1])]
            data.append(temp)
    data1 = np.array(data)
    data_new = meanshift(data1, 2.8, 1e-6)
    label = group(data_new, 1e-6)

    data2 = pd.DataFrame(data, columns=['x', 'y'])
    label2 = pd.DataFrame(label, columns=['label'])

    data_set = pd.concat([data2, label2], axis=1)

    colors = ['black','green', 'red', 'blue','yellow', 'hotpink']

    for loopi in range(int(label.max())+1):
        plt.scatter(data_set[data_set.label==loopi].x, data_set[data_set.label==loopi].y, color=colors[loopi])
    plt.show()

数据集

10.91079039	8.389412017
9.875001645	9.9092509
7.8481223	10.4317483
8.534122932	9.559085609
10.38316846	9.618790857
8.110615952	9.774717608
10.02119468	9.538779622
9.37705852	9.708539909
7.670170335	9.603152306
10.94308287	11.76207349
9.247308233	10.90210555
9.54739729	11.36170176
7.833343667	10.363034
10.87045922	9.213348128
8.228513384	10.46791102
12.48299028	9.421228147
6.557229658	11.05935349
7.264259221	9.984256737
4.801721592	7.557912927
6.861248648	7.837006973
13.62724419	10.94830031
13.6552565	9.924983717
9.606090699	10.29198795
12.43565716	8.813439258
10.0720656	9.160571589
8.306703028	10.4411646
8.772436599	10.84579091
9.841416158	9.848307202
15.11169184	12.48989787
10.2774241	9.85657011
10.1348076	8.892774944
8.426586093	11.30023345
9.191199877	9.989869949
5.933268578	10.21740004
9.666055456	10.68814946
5.762091216	10.12453436
5.224273746	9.98492559
10.26868537	10.31605475
10.92376708	10.93351512
8.935799678	9.181397458
2.978214427	3.835470435
4.91744201	2.674339991
3.024557256	4.807509213
3.019226157	4.041811881
4.131521545	2.520604653
0.411345842	3.655696597
5.266443567	5.594882041
4.62354099	1.375919061
5.67864342	2.757973123
3.905462712	2.141625079
8.085352646	2.58833713
6.852035583	3.610319053
4.230846663	3.563377115
6.042905325	2.358886853
4.20077289	2.382387946
4.284037893	7.051142553
3.820640884	4.607385052
5.417685111	3.436339164
8.21146303	3.570609885
6.543095544	-0.150071185
9.217248861	2.40193675
6.673038102	3.307612539
4.043040861	4.849836388
3.704103266	2.252629794
4.908162271	3.870390681
5.656217904	2.243552275
5.091797066	3.509500134
6.334045598	3.517609974
6.820587567	3.871837206
7.209440437	2.853110887
2.099723775	2.256027992
4.720205587	2.620700716
6.221986574	4.665191116
5.076992534	2.359039927
3.263027769	0.652069899
3.639219475	2.050486686
7.250113206	2.633190935
4.28693774	0.741841034
4.489176633	1.847389784
6.223476314	2.226009922
2.732684384	4.026711236
6.704126155	1.241378687
6.406730922	6.430816427
3.082162445	3.603531758
3.719431124	5.345215168
6.190401933	6.922594241
8.101883247	4.283883063
2.666738151	1.251248672
5.156253707	2.957825121
6.832208664	3.004741194
-1.523668483	6.870939176
-6.278045454	5.054520751
-4.130089867	3.308967776
-2.298773883	2.524337553
-0.186372986	5.059834391
-5.184077845	5.32761477
-5.260618656	6.373336994
-4.067910691	4.56450199
-4.856398444	3.94371169
-5.169024046	7.199650795
-2.818717016	6.775475264
-3.013197129	5.307372667
-1.840258223	2.473016216
-3.806016495	3.099383642
-1.353873198	4.60008787
-5.422829607	5.540632064
-3.571899549	6.390529804
-4.037978273	4.70568099
-1.110354346	4.809405537
-3.8378779	6.029098753
-6.55038578	5.511809253
-5.816344971	7.813937668
-4.626894927	8.979880178
-3.230779355	3.295580582
-4.333569224	5.593364339
-3.282896829	6.590185797
-7.646892109	7.527347421
-6.461822847	5.62944836
-6.368216425	7.083861849
-4.284758729	3.842576327
-2.29626659	7.288576999
1.101278199	6.548796127
-5.927942727	8.655087775
-3.954602311	5.733640188
-3.160876539	4.267409415

文章导引列表：

机器学习

1. 小瓜讲机器学习——分类算法（一）logistic regression(逻辑回归)算法原理详解
2. 小瓜讲机器学习——分类算法（二）支持向量机（SVM）算法原理详解
3. 小瓜讲机器学习——分类算法（三）朴素贝叶斯法（naive Bayes）
4. 小瓜讲机器学习——分类算法（四）K近邻法算法原理及Python代码实现
5. 小瓜讲机器学习——分类算法（五）决策树算法原理及Python代码实现
6. 小瓜讲机器学习——聚类算法（一）K-Means算法原理Python代码实现
7. 小瓜讲机器学习——聚类算法（二）Mean Shift算法原理及Python代码实现
8. 小瓜讲机器学习——聚类算法（三）DBSCAN算法原理及Python代码实现

---

数据分析

1. 小呆学数据分析——使用pandas中的merge函数进行数据集合并
2. 小呆学数据分析——使用pandas中的concat函数进行数据集堆叠
3. 小呆学数据分析——pandas中的层次化索引
4. 小呆学数据分析——使用pandas的pivot进行数据重塑
5. 小呆学数据分析——用duplicated/drop_duplicates方法进行重复项处理
6. 小呆学数据分析——缺失值处理（一）
7. 小呆学数据分析——异常值判定与处理（一）
8. 小瓜讲数据分析——数据清洗

---

数据可视化

1. 小瓜讲数据分析——数据可视化工程（matplotlib库使用基础篇）
2. 小瓜讲matplotlib高级篇——坐标轴设置（坐标轴居中、坐标轴箭头、刻度设置、标识设置）