目录
- 聚类
- 基本原理
- 基本原则
- 步骤和过程
- 系统聚类
- K-Means
- DBSCAN
聚类
基本原理
多元统计分析中的聚类分析方法既可以对样本进行分类(记为 \(Q\) 型分类),也可以对反映事物特征的指标或变量(记为 \(R\) 型分类)进行分类。两种分类时对等的。在算法上没有任何区别。此处主要以 \(Q\)
人们往往根据事物之间的距离远近或相似程度来判定类别。个体与个体之间的距离越近,其相似性可能也越大,是同类的可能性也越大,聚在一起形成类别的可能性也就越大。因此有了巨累分析的基本原则。
基本原则
聚类过程所依据的距离主要与明氏距离、马氏距离等几大类。
设样本数据可以用如下矩阵形式表示
\[X= \left( \begin{array}{} x_{11} & x_{12} & \cdots & x_{1p} \\ x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots \\ x_{n1} & x_{n2} & \cdots & x_{np} \\ \end{array} \right),记为X=\{x_{ij}\}_{n\times p} \]
设 \(d_{ij}\) 表示第 \(i\) 个样本与第 \(j\) 个样本之间的距离。如果 \(d_{ij}\)满足以下4个条件,则称其为距离
- \(d_{ij} \geq 0\),对于一切 \(i,j\);
- \(d_{ij} = 0\),等价于 \(i=j\);
- \(d_{ij} = d_{ji}\),对于一切 \(i,j\);
- \(d_{ij} \leq d_{ik}+d_{kj}\),对于一切 \(i,j,k\);
第1个条件表明聚类分析中的距离是非负的;第2个条件表明个体自身与自身的距离为0;第3个条件表明距离的对等性,即A和B之间的距离与B和A之间的距离是一致的;第4个条件表明两点之间直线距离是最小的。
明氏距离是最常用的距离度量方法之一,其计算公式为
\[d_{ij}(q) = (\sum_{k=1}^{p}{|x_{ik}-x_{jk}|^q})^{1/q} \]
有如下几种典型情况
- 当 \(q=1\)时,\(d_{ij}(1) = \sum_{k=1}^{p}{|x_{ik}-x_{jk}|}\) 称为绝对距离
- 当 \(q=2\)时,\(d_{ij}(2) = (\sum_{k=1}^{p}{|x_{ik}-x_{jk}|^2})^{1/2}\)称为欧氏距离
- 当\(q=1\)时,\(d_{ij}(\infty) = \max_{1\leq k \leq p}{|x_{ik}-x_{jk}|}\)称为车比雪夫距离
但是明氏距离的大小与个体指标的观测单位有关,没有考虑指标之间的相关性。为克服此缺点,可以考虑马氏距离进行改造。马氏距离 是由协方差矩阵计算出来的相对距离,具体计算公式如下
\[d_{ij} = (X_i-X_j)^{'}\Sigma^{-1}(X_i-X_j) \]
其中,\(\Sigma\)
除了最短距离原则进行分类之外,还可以采用相关系数、相似系数、匹配系数等指标来衡量个体之间的相似性,以此为依据进行分类。
在分类过程中,为了便于分析,有如下3个重要原则:
- 同质性原则:同一类中个体之间有较大的相似性
- 互斥性原则:不同类中的个体差异很大
- 完备性原则:每个个体在同一次分类过程中,能且只能分在一个类别中
实际应用中,以最短距离原则进行系统聚类比较常用。
步骤和过程
实际数据分析中,一般会首先选取样本点之间的距离测算方式,然后按照最近或最像原则进行分类,当类别中有多个样本点的时候,就会设计到样本类别和类别之间距离确定,则个是不可避免的流程。
系统聚类
依据多个指标进行聚类分析时,可依照如下流程进行系统聚类分析
1.选择并计算距离
2.每个个体自成一类
3.按照一定原则把最近的两类合并为新类
4.计算当前类与新类的距离
5.判断是否仅有一类:若是,则依据谱系聚类图确定类的数目,若否,继续3/4/5的循环
当类别中多个样本点时,会涉及到类别与类别之间的距离定义方法。社会经济研究中有如下几种常用的类间距确定方法:
- 最短距离法(single linkage)
最短距离法又称为单连接聚类法。如果有两类 \(G_p,G_q\) 聚成新类 \(G_n\),在最短距离法中新类 \(G_n\)与其他的任意类 \(G_k\)
\[D_{kn} = Min(D_{kp}, D_{kq}) \]
其中, \(D_{kp},D_{kq}\) 是用来衡量原有类别 \(G_p,G_q\) 中各样本点与任意类\(G_k\)中各样本点的点间距离。
即如果新类与其他类别之间存在多个点与点之间的距离,则取这些距离当中最小者作为两类的距离,即在进行聚类的过程中应以最小点间剧里作为并类的依据。
- 最长距离法(complete method)
最长距离法也叫完全连接法。如果有两类 \(G_p,G_q\) 聚成新类 \(G_n\),在最长距离法中新类 \(G_n\)与其他的任意类 \(G_k\)
\[D_{kn} = Max(D_{kp}, D_{kq}) \]
即如果新类与其他类别之间存在多个点与点之间的距离,则取这些距离当中最大者作为两类的距离。
在及性能分类的过程中,首先以最小距离原则把最近的两个样本合并为一个新类,其余各样本自身自成一类。则刚合并的新类与其他类别之间按最大距离原则确定之后,再按照最小原则把类别之间距离最小的合并为一个新类,以此类推,直至把所有样本归为一类。
- 中间距离法(median method)
如果有两类 \(G_p,G_q\) 聚成新类 \(G_n\),在中间距离法中新类 \(G_n\)与其他的任意类 \(G_k\)
\[D_{kn} = \frac{D_{kp} + D_{kq}}{2}-\frac{D_{pq}}{4} \]
如果新类与其他类别之间存在多个点距离,按照上述公式计算的结果作为两类的距离,然后按照最小距离原则把类别之间距离最小的两类合并为一类,直至把所有样本归为一类。
- 重心法(centroid method)
在以上方法中,没有考虑每一类中所包含的样本数目。重心法可克服此缺点。
如果有两类 \(G_p,G_q\)
\[D_{pq} = ||\bar{X}_p-\bar{X}_q||^2 \]
其中,\(\bar{X}_p,\bar{X}_q\)
该距离即为两类重心(通常可用类内样本各指标的均值表示)之间的欧式距离平方。新类 \(G_n\)与其他的任意类 \(G_k\)
\[D_{kn}=\frac{N_pD_{kp}+N_qD_{kq}}{N_n}-\frac{N_pN_qD_{pq}}{N_n^2} \]
其中,\(N_i,i=p,q,n\)
即在并类的过程中,以最小重心距离作为依据并类,直至把所有样本归为一类。
- 类平均法(average linkage)
重心法虽有较好代表性,但并未充分利用各个样本的信息,可以用不同类的样本点两两之间的平均距离作为类间距离。如果有类\(G_p,G_q\),可以计算每类中每队样本点之间的平均距离,即
\[D_{pq} = \frac{1}{N}\sum_{i\in G_p}\sum_{j\in G_q}{d(x_i,x_j)} \]
若 \(d(x,y)=|x-y|^2\), 新类 \(G_n\)与其他的任意类 \(G_k\)
\[D_{kn}=\frac{N_pD_{kp}+N_qD_{kq}}{N_n} \]
即在并类的过程中,以类别样本点之间的平均距离作为依据并类,直至把所有样本归为一类。
- 离差平方和法(WARD)
离差平方和法又被称为WARD最小方差法。它的思想来源于方差分析,即如果聚类得恰当,类内样本点之间的理茬平方和应较小,而类间离差平方和应当大。该法要求样品间距离必须采用欧式距离。离差平方和法定义类间距离平方为
\[D_{pq}^2 = S_n^2-S_p^2-S_q^2 \]
其中,\(S_n^2\) 是类 \(G_p,G_q\)合并成 \(G_n\)
当观测距离 \(d(x,y)=||x-y||^2/2\) 时,则新类 \(G_n\)与其他的任意类 \(G_k\)
\[D_{kn}=\frac{(N_k+N_p)D_{kp}+(N_k+N_q)D_{kq}-N_kD_{pq}}{N_k+N_n} \]
使用离差平方和法进行分类时,先让每个样本自身各成一类,然后并类,每并一类离差平方和就要增大,选择使其增量最小的两类合并,直到所有的样品聚为一类。
使用哪种方法来对样本数据的类间距,可以参考经验分析及研究对象的特征和研究要求,进行经验判断。如可根据给定距离的阈值或是通过计算相应的统计量来判定,如Demirmen提出的一些在决定聚类方法取舍时应遵循的原则:
1.任何类必须在邻近的各类中是突出的,即各类重心(常取平均数衡量)之间应该有最大的距离
2.确定的类中,各类所包含的元素都不宜过分多
3.聚类数目应符合实际
4.当用许多方法进行分类时,应选出现次数最多的那种分类结果
实际应用中,通常是多用几种不同的类间距方法进行聚类,然后根据现有理论和研究要求挑出合适的分类类别
示例
import pandas as pd
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
# 通过专家和消费者打分形式,对13款鼠标做参数数据,进行聚类分析
# 指定为黑体中文字体,防止中文乱码
plt.rcParams["font.sans-serif"] = ["Heiti TC"]
# 解决保存图像是负号'-'显示为方块的问题
plt.rcParams['axes.unicode_minus'] = False
mc = pd.read_csv('./data/mouse_cluster.csv')
print(mc.head(5))
# brand Touch Chips Driver Compatibility Game
# 0 Brand1 7.5 17.5 7.0 8.0 8.0
# 1 Brand2 7.5 19.5 7.0 7.0 9.0
# 2 Brand3 8.5 18.0 8.5 8.0 9.5
# 3 Brand4 9.0 18.5 8.5 8.0 9.5
# 4 Brand5 7.0 14.0 6.5 7.0 7.5
# 系统聚类
import scipy.cluster.hierarchy as hc
# 采用类间距的类平均法
z = hc.linkage(mc.iloc[:, 1:], method='average')
# methord可以是:single,complete,average,weighted,centroid,median, ward等
# 链接矩阵z是一个展示并类过程的数组,每一行表示一次并类步骤,第三列表示并类距离
print(z)
# [[ 2. 3. 0.70710678 2. ]
# [ 0. 9. 0.8660254 2. ]
# [ 6. 7. 1. 2. ]
# [ 8. 11. 1. 2. ]
# [ 5. 16. 1.32287566 3. ]
# [14. 15. 1.43328552 4. ]
# [10. 17. 1.55409255 4. ]
# [12. 18. 1.64860349 5. ]
# [ 4. 19. 2.28324551 5. ]
# [13. 20. 2.48994247 7. ]
# [ 1. 22. 2.81192434 8. ]
# [21. 23. 3.34347471 13. ]]
# 方法一:绘制系统聚类谱系图
dd = hc.dendrogram(z, orientation='right', labels=list(mc.iloc[:, 0]))
# orientation可以选top,bottom,left,right分别表示对应方向的谱系图
plt.show()
# 在水平谱系图中任意地点画一条竖线(x=2.6),该直线与图中的横线有多少个交代呢就可以把样本分为多少类,
# 分类的依据便是这条竖线的直线所对应的横轴距离,与每个交点相连的样本同属一类。
# 根据连接线的长短,可以判定出各个样本并类的过程。
# 方法二:fcluster可对由原始数据计算的链接矩阵进行聚类
# 设置并类阈值来对计算出来的链接矩阵进行分组
res = hc.fcluster(z, 2.6, criterion='distance')
print(res) # [2 3 2 2 1 1 2 2 1 2 1 1 2]
# 数组元素分别表示各行数据对应的聚类分组,相同数字的行所代表的样本点是同一个类别
# 方法三:fclusterdata可对原始数据进行直接聚类
res = hc.fclusterdata(mc.iloc[:, 1:], 2.6, criterion='distance', metric='euclidean', method='average')
print(res) # [2 3 2 2 1 1 2 2 1 2 1 1 2]
# 方法四:sklearn
from sklearn.cluster import AgglomerativeClustering as AC
cm = AC(n_clusters=3, linkage='average', affinity='euclidean')
cl = cm.fit(mc.iloc[:, 1:])
print(cl.labels_) # [0 2 0 0 1 1 0 0 1 0 1 1 0]
# 由于有一个类别所包含的样本量仅为1,考虑使用其他类间距方法重新分裂
dd_ward = hc.dendrogram((hc.linkage(mc.iloc[:, 1:], method='ward')), orientation='right', labels=list(mc.iloc[:, 0]))
plt.show()
K-Means
在系统聚类中,研究者事先并不知道对样本据称多少类别。在有些情况下,研究者对于研究的对象事先知道分为几类,即已知类别的个数 \(k\),只是不知这些类别中的具体样本,此时可以考虑快速聚类方法进行聚类。
快速聚类一般用于大样本情况下的样品聚类。Anderberg提出了最近中心归类法为基本的快速聚类:
1.选择k各观测值组成初始类别并作为聚类种子;
2.找出聚类种子的中心;
3.把每一个观测值依据最小欧式距离(即观测值与聚类种子中心的欧式距离)原则归入各类,构成暂时的类别;
4.计算每个暂时类别中各个变量的均值,以此作为新的类别中心
5.再一次把每一个观测值根据最小欧式距离原则归入各类,购成新的暂时的类别
6.重复3~5的过程,中心的迭代标准达到要求时,聚类过程结束。
上述过程也被称为K-Means聚类。
示例
import pandas as pd
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
# 通过专家和消费者打分形式,对13款鼠标做参数数据,进行聚类分析
# 指定为黑体中文字体,防止中文乱码
plt.rcParams["font.sans-serif"] = ["Heiti TC"]
# 解决保存图像是负号'-'显示为方块的问题
plt.rcParams['axes.unicode_minus'] = False
mc = pd.read_csv('./data/mouse_cluster.csv')
print(mc.head(5))
# brand Touch Chips Driver Compatibility Game
# 0 Brand1 7.5 17.5 7.0 8.0 8.0
# 1 Brand2 7.5 19.5 7.0 7.0 9.0
# 2 Brand3 8.5 18.0 8.5 8.0 9.5
# 3 Brand4 9.0 18.5 8.5 8.0 9.5
# 4 Brand5 7.0 14.0 6.5 7.0 7.5
# k-means聚类
# 方法一:scipy的kmeans2
from scipy.cluster.vq import kmeans2
res = kmeans2(mc.iloc[:, 1:], 3)
print(res)
# (array([[ 7.2 , 15.6 , 6.3 , 7.5 , 7.3 ],
# [ 8.33333333, 18.66666667, 8. , 7.66666667, 9.33333333],
# [ 7.6 , 17.2 , 7.8 , 7.7 , 7.9 ]]),
# array([2, 1, 1, 1, 0, 0, 2, 2, 0, 2, 0, 0, 2], dtype=int32))
# 数组1表示最后一次迭代的聚类中心;数组2表示每行数据的类别标识,标识相同的表示同类别
# 方法二:scipy的kmeans
# 其要求在迭代中心值小于等于阈值时进行扭曲变换使得迭代过程得以停止。
# 注意:使用该函数进行kmeans聚类前,需要将原始数据进行洗白(即使得各列数据均具有标准方差,
# 亦即每列数据除以其在所有观测值中的标准差以给出单位方差)
from scipy.cluster.vq import kmeans, whiten
res = kmeans(whiten(mc.iloc[:, 1:]), 3)
# 给出最终的聚类中心即扭曲值
print(res)
# (array([[12.85063882, 13.22354999, 8.63685747, 14.90223347, 9.26392533],
# [14.74663471, 13.59604435, 9.32232235, 16.44384383, 10.83166654],
# [12.27883053, 11.97303319, 7.20717358, 15.8565637 , 8.55131569]]),
# 1.3410757207693236)
# 方法三:sklearn的KMeans
from sklearn.cluster import KMeans
mc_km = KMeans(n_clusters=3).fit(mc.iloc[:, 1:])
print(mc_km.labels_) # [2 1 1 1 0 0 2 2 0 2 0 0 2]
# 把brand变量与上述类别拼接起来查看每个品牌鼠标的聚类情况
cluster_by_Kmeans = pd.concat([mc['brand'], pd.DataFrame(mc_km.labels_, columns=['cluster'])], axis=1)
print(cluster_by_Kmeans.T)
# 0 1 2 3 ... 9 10 11 12
# brand Brand1 Brand2 Brand3 Brand4 ... Brand10 Brand11 Brand12 Brand13
# cluster 2 1 1 1 ... 2 0 0 2
DBSCAN
DBSCAN即具有噪声的基于密度的空间聚类方法(density-based spatial clustering of applications with noise)。该算法利用基于密度的聚类的概念,即要求聚类空间中的一定区域内所包含的对象(点或其他空间对象)的数目不小于某一给定阈值。DBSCAN算法的显著优点是聚类速度快且能够有效处理噪声点和发现任意形状的空间聚类。但是由于它直接对整个数据库进行操作且进行聚类时使用了一个全局性的表征密度的参数,因此也具有两个比较明显的弱点:
- 当数据量增大时,系统开销也很大;
- 当空间剧烈的密度不均匀、聚类间距离差想差很大时,聚类质量较差。
DBSCAN算法有2个输入参数:一个是半径(Eps),表示以给定 \(P\) 为中心的圆形邻域的范围;另一个是以点 \(P\) 为中心的邻域内最少点的数量(MinPts),这2个参数的计算都来自经验。如果满足:以点 \(P\) 为中心、半径为 \(Eps\) 的邻域内的点的个数不少于 \(MinPts\),则称点 \(P\)
其主要基本概念如下:
- \(\varepsilon\) 邻域:给定对象半径 \(\varepsilon\) 内的区域称为该对象的 \(\varepsilon\)
- 核心对象:如果给定对象 \(\varepsilon\) 邻域内的样本点数大于等于 \(MinPts\),则称该对象为核心对象。
- 直接密度可达:给定一个对象集合 \(D\),如果 \(p\) 在 \(q\) 的 \(\varepsilon\) 邻域内,且\(q\) 是一个核心对象,则对象 \(p\) 从对象 \(q\)
- 密度可达:对于样本集合 \(D\),如果存在一个对象链 \(p_1,p_2,\cdots,p_n;p_1=q,p_n=p\),对于 \(p_1\in D(1\leq i \leq n)\),\(p_{i+1}\) 是从 \(p_i\) 关于 \(\varepsilon\) 和 \(MinPts\) 直接密度可达,则对象 \(p\) 是从对象\(q\) 关于\(\varepsilon\) 和 \(MinPts\)
- 密度相连:如果存在对象 \(o \in D\),使对象 \(p,q\) 都是从 \(o\) 关于\(\varepsilon\) 和 \(MinPts\) 密度可达的,那么对象 \(p\) 到 \(q\) 是于\(\varepsilon\) 和 \(MinPts\)
可以发现,密度可达是直接密度可达的传递闭包,并且这种关系是非对称的。只有核心对象之间相互密度可达。密度相连是对称关系。DBSCAN的目的就是要找到密度相连对象的最大集合。
其算法的具体聚类过程如下:
扫描整个数据集,找到任意一个核心点,对该核心点进行扩充。扩充的方法是寻找从该核心点出发的所有密度相连的数据点。遍历该核心点的 \(\varepsilon\)
示例
import pandas as pd
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
# 通过专家和消费者打分形式,对13款鼠标做参数数据,进行聚类分析
# 指定为黑体中文字体,防止中文乱码
plt.rcParams["font.sans-serif"] = ["Heiti TC"]
# 解决保存图像是负号'-'显示为方块的问题
plt.rcParams['axes.unicode_minus'] = False
mc = pd.read_csv('./data/mouse_cluster.csv')
print(mc.head(5))
# brand Touch Chips Driver Compatibility Game
# 0 Brand1 7.5 17.5 7.0 8.0 8.0
# 1 Brand2 7.5 19.5 7.0 7.0 9.0
# 2 Brand3 8.5 18.0 8.5 8.0 9.5
# 3 Brand4 9.0 18.5 8.5 8.0 9.5
# 4 Brand5 7.0 14.0 6.5 7.0 7.5
# DBSCAN
from sklearn.cluster import DBSCAN
mc_DB = DBSCAN(eps=1.5, min_samples=1).fit(mc.iloc[:, 1:])
cluster_by_DBSCAN = pd.concat([mc['brand'], pd.DataFrame(mc_DB.labels_, columns=['cluster'])], axis=1)
print(cluster_by_DBSCAN.T)
# 0 1 2 3 ... 9 10 11 12
# brand Brand1 Brand2 Brand3 Brand4 ... Brand10 Brand11 Brand12 Brand13
# cluster 0 1 0 0 ... 0 3 3 0