第九章 聚类分析与判别分析在实际工作中，我们经常遇到分类问题.若事先已经建立类别，则使用判别分析，若事先没有建立类别，则使用聚类分析.聚类分析主要是研究在事先没有分类的情况下，如何将样本归类的方法.聚类分析的内容包含十分广泛，有系统聚类法、动态聚类法、分裂法、最优分割法、模糊聚类法、图论聚类法、聚类预报等多种方法.§9.1 聚类分析基本知识介绍在MA TLAB 软件包中，主要使用的是系统聚类法.系统聚类法是聚类分析中应用最为广泛的一种方法.它的基本原理是：首先将一定数量的样品(或指标)各自看成一类，然后根据样品(或指标)的亲疏程度，将亲疏程度最高的两类合并，然后重复进行，直到所有的样品都合成一类.衡量亲疏程度的指标有两类：距离、相似系数.一、常用距离1）欧氏距离假设有两个n 维样本和),,,(112111n x x x x =),,,(222212n x x x x =，则它们的欧氏距离为∑=-=nj j jx xx x d 122121)(),(2）标准化欧氏距离假设有两个n 维样本),,,(112111n x x x x =和),,,(222212n x x x x =，则它们的标准化欧氏距离为12(,)sd x x ==其中：D 表示n 个样本的方差矩阵，),,,(22221n diag D σσσ =，2j σ表示第j 列的方差,即每个指标的方差。
若每个指标的均值相等，方差相同，则有12(,)sd x x ==3）马氏距离假设共有n 个指标，第i 个指标共测得m 个数据(要求n m >)：⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=im i i i x x x x M 21于是，我们得到n m ⨯阶的数据矩阵),,,(21n x x x X =，每一行是一个样本数据.n m ⨯阶数据矩阵X 的n n ⨯阶协方差矩阵记做)(X Cov .两个n 维样本),,,(112111n x x x x =和),,,(222212n x x x x =的马氏距离如下：T x x X Cov x x x x mahal )())()((),(2112121--=-马氏距离考虑了各个指标量纲的标准化，是对其它几种距离的改进.马氏距离不仅排除了量纲的影响，而且合理考虑了指标的相关性.缺点：夸大了微小变量的作用，因为协方差阵的存在，有一定的不稳定性 4）布洛克距离（曼哈顿距离） 求方形建筑的最短距离。
两个n 维样本),,,(112111n x x x x =和),,,(222212n x x x x =的布洛克距离如下：∑=-=nj j j x x x x b 12121||),(5）切比雪夫距离按照国际象棋中国王走法，求任意格A 到格B 的最少步数|max ),(21|21j j x x x x c -=6）闵可夫斯基距离（明氏距离）两个n 维样本),,,(112111n x x x x =和),,,(222212n x x x x =的闵可夫斯基距离如下：pn j pj j x x x x m 112121||),(⎪⎪⎭⎫ ⎝⎛-=∑= 注：1=p 时是布洛克距离；2=p 时是欧氏距离，∞→p 时是切比雪夫距离 缺点：没有考虑量纲问题，同时没有考虑各个分量的分布可能不同7）余弦距离1212(,)11cos ,T d x x x x ⎛⎫==-<> ⎝ 这是受相似性几何原理启发而产生的一种标准，在识别图像和文字时，常用夹角余弦为标准.8）相似距离（用1减去皮尔森相关系数）关于样本的协方差以及相关系数1212(,)111T x x x x d x x ρ===-二、MATLAB 中常用的计算距离的函数假设我们有n m ⨯阶数据矩阵),,,(21n x x x x =，每一行是一个样本数据. 在MATLAB 中计算样本点之间距离的内部函数为：y=pdist(x，‘euclid’) 计算样本点之间的欧氏距离y=pdist(x,'seuclid') 计算样本点之间的标准化欧氏距离y=pdist(x,'mahal') 计算样本点之间的马氏距离y=pdist(x,'cityblock') 计算样本点之间的布洛克距离y=pdist(x,'minkowski') 计算样本点之间的闵可夫斯基距离y=pdist(x,'minkowski',p) 计算样本点之间的参数为p的闵可夫斯基距离y=pdist(x,'cosine') 计算样本点之间的余弦距离y=pdist(x,'correlation') 计算样本点之间的相似距离y为距离向量，共有m*（m-1）/2个距离元素。
另外，内部函数yy=squareform(y)表示将样本点之间的距离用矩阵的形式输出.三、常用的聚类方法常用的聚类方法主要有以下几种：最短距离法、最长距离法、中间距离法、质心法、内平方距离法等等.四、创建系统聚类树假设已经得到样本点之间的距离y，可以用linkage函数创建系统聚类树，格式为z=linkage(y).其中：z为一个包含聚类树信息的(m-1) 3的矩阵.例如：z=2.000 5.000 0.23.0004.000 1.28则z的第一行表示第2、第5样本点连接为一个类，它们距离为0.2；则z的第二行表示第3、第4样本点连接为一个类，它们距离为1.28.在MA TLAB中创建系统聚类树的函数为z=linkage(y，‘single’) 表示用最短距离法创建系统聚类树z=linkage(y,'complete') 表示用最长距离法创建系统聚类树z=linkage(y,'average') 表示用平均距离法创建系统聚类树z=linkage(y,'centroid') 表示用质心距离法创建系统聚类树z=linkage(y,'ward') 表示用内平方距离法创建系统聚类树§9.2 聚类分析示例例1 在MA TLAB中写一个名为opt_linkage_1的M文件：x=[3 1.7;1 1;2 3; 2 2.5; 1.2 1; 1.1 1.5; 3 1];y=pdist(x,'mahal');yy=squareform(y)；z=linkage(y,'centroid')h=dendrogram(z)存盘后按F5键执行，得到结果如下：yy =0 2.3879 2.1983 1.6946 2.1684 2.2284 0.88952.3879 0 2.6097 2.0616 0.2378 0.6255 2.37782.1983 2.6097 0 0.6353 2.5522 2.0153 2.98901.69462.0616 0.6353 0 1.9750 1.5106 2.41722.1684 0.2378 2.5522 1.9750 0 0.6666 2.14002.2284 0.6255 2.0153 1.5106 0.6666 0 2.45170.8895 2.3778 2.9890 2.4172 2.1400 2.4517 0z =2.0000 5.0000 0.23786.0000 8.0000 0.63533.00004.0000 0.63531.0000 7.0000 0.88959.0000 10.0000 2.106311.0000 12.0000 2.0117按重心距离法得到的系统聚类树为其中：h=dendrogram(z)表示输出聚类树形图的冰状图.一、根据系统聚类树创建聚类假设已经求出系统聚类树z，我们根据z来创建聚类，使用cluster函数.例2 在MA TLAB中写一个名为opt_cluster_1的M文件：x=[3 1.7;1 1;2 3; 2 2.5; 1.2 1; 1.1 1.5; 3 1];y=pdist(x,'mahal');yy=squareform(y)z=linkage(y,'centroid')h=dendrogram(z)t=cluster(z,3) %%%% cluster(Z,'cutoff',c) cluster(Z,'maxclust',n)其中：“t=cluster(z,3)”表示分成3个聚类，需要分成几个由人工选择.存盘后按F5键执行，得到结果如下：t =3122113即第1、第7样本点为第3类，第2、第5、第6样本点为第1类，第3、第4样本点为第2类.二、根据原始数据创建分类在MA TLAB软件包中，内部函数clusterdata对原始数据创建分类，格式有两种：1）clusterdata(x,a)，%%%%%clusterdata(x,’cutoff’,a)其中0<a<2，表示在系统聚类树中inconsistency coefficient小于a的样本点归结为一类；2）clusterdata(x,b)，%%%%%clusterdata(x,’maxclust’,b)其中b>=2(b>1)是整数，表示将原始数据x 分为b类.%%%%T = clusterdata(X,cutoff) is the same asY = pdist(X,'euclid');Z = linkage(Y,'single');T = cluster(Z,'cutoff',cutoff);%%%%T = clusterdata(X,maxclust) is the same asY = pdist(X,'euclid');Z = linkage(Y,'single');T = cluster(Z,'maxclust',maxclust);例3 在MA TLAB中写一个名为opt_clusterdata_1的M文件：x=[3 1.7;1 1;2 3; 2 2.5; 1.2 1; 1.1 1.5; 3 1];t= clusterdata(x,0.5)z= clusterdata(x,3)存盘后按F5键执行，得到结果如下：t =4322314z =2311332其中：t的结果表示距离小于0.5的样本点归结为一类，这样，共有四类，第1类：样本点6；第2类：样本点3、4；第3类：样本点2、5；第4类：样本点1、7.而z的结果表示首先约定将原始数据x分为3类，然后计算，结果如下：第1类：样本点3、4；第2类：样本点1、7；第3类：样本点2、5、6.利用内部函数clusterdata对原始数据创建分类，其缺点是不能更改距离的计算法.比较好的方法是分步聚类法.三、分步聚类法假设有样本数据矩阵x，第一步对于不同的距离，利用pdist函数计算样本点之间的距离：y1=pdist(x)y2=pdist(x,'seuclid')y3=pdist(x,'mahal')y4=pdist(x,'cityblock')第二步计算系统聚类树以及相关信息：z1=linkage(y1)z2=linkage(y2)z3=linkage(y3)z4=linkage(y4)第三步利用cophenet函数计算聚类树信息与原始数据的距离之间的相关性，这个值越大越好：t1=cophenet(z1,y1)t2=cophenet(z2,y2)t3=cophenet(z3,y3)t4=cophenet(z4,y4)注：z在前，y在后，顺序不能颠倒.第四步选择具有最大的cophenet值的距离进行分类.%%%%%利用函数cluster(x,a)对数据x进行分类，其中0<a<2，表示在系统聚类树中inconsistency coefficient 小于a的样本点归结为一类.例4 在MA TLAB中写一个名为opt_cluster_2的M文件：x=[3 1.7;1 1;2 3; 2 2.5; 1.2 1; 1.1 1.5; 3 1];y1=pdist(x);y2=pdist(x,'seuclid');y3=pdist(x,'mahal');y4=pdist(x,'cityblock');z1=linkage(y1);z2=linkage(y2);z3=linkage(y3);z4=linkage(y4);t1=cophenet(z1,y1)t2=cophenet(z2,y2)t3=cophenet(z3,y3)t4=cophenet(z4,y4)存盘后按F5键执行，得到结果如下：t1 =0.9291t2 =0.9238t3 =0.9191t4 =0.9242结果中t1=0.9291最大，可见此例利用欧式距离最合适.于是，在MA TLAB中另写一个名为opt_cluster_3的M文件：x=[3 1.7;1 1;2 3; 2 2.5; 1.2 1; 1.1 1.5; 3 1];y1=pdist(x);z1=linkage(y1)w=inconsistent(z1)；存盘后按F5键执行，得到结果如下：z1 =2.0000 5.0000 0.20003.00004.0000 0.50006.0000 8.0000 0.50991.0000 7.0000 0.70009.0000 11.0000 1.280610.0000 12.0000 1.3454矩阵z1的第1行表示样本点2、5为一类，在系统聚类树上的距离为0.2，其它类推.在MA TLAB中另写一个名为opt_cluster_4的M文件：x=[3 1.7;1 1;2 3; 2 2.5; 1.2 1; 1.1 1.5; 3 1];y1=pdist(x);z1=linkage(y1)b1=cluster(z1,0.5)存盘后按F5键执行，得到结果如下：b1 =4322314结果表示将原始数据x分为4类，第1类：样本点6；第2类：样本点3、4；第3类：样本点2、5；第4类：样本点1、7.例5下表是1999年中国省、自治区的城市规模结构特征的一些数据，可通过聚类分析将这些省、自治区进行分类，具体过程如下：省、自治区首位城市规模（万人）城市首位度四城市指数基尼系数城市规模中位值（万人）京津冀699.70 1.437 10.936 40.780 410.880山西179.46 1.898 2 1.000 60.587 011.780内蒙古111.13 1.418 00.677 20.515 817.775辽宁389.60 1.918 20.854 10.576 226.320吉林211.34 1.788 0 1.079 80.456 919.705黑龙江259.00 2.305 90.341 70.507 623.480浙江 139.29 1.871 2 0.885 8 0.453 6 12.670 安徽 102.78 1.233 3 0.532 6 0.379 8 27.375 福建 108.50 1.729 1 0.932 5 0.468 7 11.120 江西 129.20 3.245 4 1.193 5 0.451 9 17.080 山东 173.35 1.001 8 0.429 6 0.450 3 21.215 河南 151.54 1.492 7 0.677 5 0.473 8 13.940 湖北 434.46 7.132 8 2.441 3 0.528 2 19.190 湖南 139.29 2.350 1 0.836 0 0.489 0 14.250 广东 336.54 3.540 7 1.386 3 0.402 0 22.195 广西 96.12 1.228 8 0.638 2 0.500 0 14.340 海南 45.43 2.191 5 0.864 8 0.413 6 8.730 川渝 365.01 1.680 1 1.148 6 0.572 0 18.615 云南 146.00 6.633 3 2.378 5 0.535 9 12.250 贵州 136.22 2.827 9 1.291 8 0.598 4 10.470 西藏 11.79 4.151 4 1.179 8 0.611 8 7.315 陕西 244.04 5.119 4 1.968 2 0.628 7 17.800 甘肃 145.49 4.751 5 1.936 6 0.580 6 11.650 青海 61.36 8.269 5 0.859 8 0.809 8 7.420 宁夏 47.60 1.507 8 0.958 7 0.484 3 9.730 新疆 128.67 3.853 5 1.621 6 0.490 1 14.470§9.3 判别分析基本知识介绍一、协方差矩阵设(,)ξη是一个二维随机变量，又|()()|E E E ξξηη--<+∞，则称()()E E E ξξηη--为ξ与η的协方差，记作(,)Cov ξη.这里我们用到的是样本的协方差，样本的协方差系数.样本协方差 )()(111y y x x n i ni i ---∑=()())nT iix x y y x y --=∑例如，在MATLAB 软件包中写一个名为opt_cov_1的M 文件：x=[1 2 3]; y=[3 2 1]; cov(x,y)存盘后按F5键执行，得到结果如下：ans =1 -1 -1 1设12(,,,)n ξξξ是n 维随机变量，定义()()ij i i j j E E E σξξξξ=--，则称111212122212n n n n nn σσσσσσσσσ⎛⎫⎪ ⎪⎪⎪⎝⎭为12(,,,)n ξξξ的协方差矩阵.如果要进行相关系数分析，则在MA TLAB 软件包中写一个名为opt_corrcoef_1的M 文件： x=[1 2 3]; y=[3 2 1]; corrcoef(x,y)存盘后按F5键执行，得到结果如下：ans =1 -1 -1 1注：两个随机变量的相关系数(,)r x y 满足1(,)1r x y -≤≤.二、基本数学原理判别分析是利用原有的分类信息，得到判别函数(判别函数是这种分类的函数关系式，一般是与分类相关的若干个指标的线性关系式)，然后利用该函数去判断未知样品属于哪一类.因此这是一个学习和预测的过程.常用的判别分析法有距离判别法、费歇尔判别法、贝叶斯判别法等. 1）距离判别法距离判别法有欧式距离法和马氏距离法等.其中，欧式距离法比较粗糙，MA TLAB 软件包中采用的是马氏距离法.假设共有n 个指标，第i 个指标共测得m 个数据(要求m n >)：⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=im i i i x x x x 21于是，我们得到n m ⨯阶的数据矩阵),,,(21n x x x X =，每一行是一个样本数据.n m ⨯阶数据矩阵X 的n n ⨯阶协方差矩阵记做)(X Cov .a ）求n 维向量12(,,,)n r r r r =到n m ⨯矩阵X 的马氏距离.定义1122(,,,)n n r X r x r x r x -=---，则n 维向量12(,,,)n r r r r =到n m ⨯矩阵X 的马氏距离的公式如下：1(,)()()()T mahal r X r X Cov X r X -=--其中：i x 表示第i 个指标i x 的算术平均值.mahal X X.做b）矩阵X到自身的马氏距离，相当于矩阵X的每一行到X的马氏距离，记做(,)这种运算可以求出矩阵X的异常值.注：马氏距离法没有考虑样本指标值的分布.例1 有8 3数据矩阵如下：X=[0.052 0.084 0.021;0.037 0.0071 0.022;0.041 0.055 0.11;0.11 0.021 0.0073;0.030 0.112 0.072;0.16 0.056 0.021;0.074 0.083 0.105;0.19 0.02 1]在MA TALB软件包中输入以下程序：X=[0.052 0.084 0.021;0.037 0.0071 0.022;0.041 0.055 0.11;0.11 0.021 0.0073;0.030 0.112 0.072;0.16 0.056 0.021;0.074 0.083 0.105;0.19 0.02 1];d=mahal(X,X)执行后得到结果：d =0.69983.90500.83101.88192.76724.24240.59916.0735猜测矩阵X中第二行、第五行、第六行、第八行数据异常.阈值2）费歇尔判别法该法是以费歇尔准则为标准来评选判别函数的.所谓费歇尔判别法，指的是较优的判别函数应该能根据待判对象的n个指标最大程度地将它所属的类与其它类区分开来.一般采用线性函数作为判别函数.基本方法是：首先假定一个线性的判别函数，然后根据已知信息对判别函数进行训练，得到函数关系式的系数，从而最终确定判别函数.该法有时会使误判次数增加，但由于使用的是线性函数，所以使用起来也比较方便.3）贝叶斯判别法贝叶斯判别法是一种概率方法，它的好处是可以充分利用先验信息，可以考虑专家的意见.应用此方法，需要事先假定样本指标值的分布(例如，多元正态分布等).在MA TLAB软件包中，将已经分类的m个数据(长度为n)作为行向量，得到一个矩阵training，每行都属于一个分类类别，分类类别构成一个整数列向量g(共有m行)，待分类的k个数据(长度为n)作为行向量，得到一个矩阵sample，然后利用classify函数进行线性判别分析(默认).它的格式为classify(sample,training,group)其中：sample与training必须具有相同的列数；group与training必须具有相同的行数；group是一个整数向量.MA TLAB内部函数classify的功能是将sample的每一行进行判别，分到training指定的类中.进一步，较复杂的格式为[class,err]=classify(sample,training,group,type)其中：class返回分类表；err返回误差比例信息；sample是样本数据矩阵；training是已有的分类数据矩阵；group是分类列向量；type有3种选择：1）type=linear(默认设置)：表示进行线性判别分析；2）type=quadratic：表示进行二次判别分析；3）type=mahalanobis：表示用马氏距离进行判别分析.§9.4 判别分析示例例1 某地大气样品污染分类表如下：在此地某大型化工厂的厂区及临近地区挑选4个有代表性的大气样本取样点，获取数据如下：求它们的污染分类.在MA TLAB软件包中写一个名为opt_class_1的M文件：training=[0.056 0.084 0.031 0.038 0.0081 0.0220.04 0.055 0.1 0.11 0.022 0.00730.05 0.074 0.041 0.048 0.0071 0.020.045 0.05 0.11 0.1 0.025 0.00630.038 0.13 0.079 0.17 0.058 0.0430.03 0.11 0.07 0.16 0.05 0.0460.034 0.095 0.058 0.16 0.2 0.0290.03 0.09 0.068 0.18 0.22 0.0390.084 0.066 0.029 0.32 0.012 0.0410.085 0.076 0.019 0.3 0.01 0.040.064 0.072 0.02 0.25 0.028 0.0380.054 0.065 0.022 0.28 0.021 0.040.048 0.089 0.062 0.26 0.038 0.0360.045 0.092 0.072 0.2 0.035 0.0320.069 0.087 0.027 0.05 0.089 0.021];group=[1;1;1;1;2;2;1;1;2;2;2;2;2;2;1];sample=[0.052 0.084 0.021 0.037 0.0071 0.0220.041 0.055 0.11 0.11 0.021 0.00730.03 0.112 0.072 0.16 0.0056 0.0210.074 0.083 0.105 0.19 0.02 1];class=classify(sample,training,group)存盘后按F5键执行，得到结果如下：class =1122即样品1、样品2为1类污染，样品3、样品4为2类污染.。
