第八章 聚类分析与判别分析 8.1 聚类分析 8.2 判别分析

8.1 聚类分析 8.1.1 聚类分析的一般概念 8.1.2 系统聚类法的基本思想和步骤 8.1 聚类分析 8.1.1 聚类分析的一般概念 8.1.2 系统聚类法的基本思想和步骤 8.1.3 用CLUSTER过程和TREE过程进行系统聚类 8.1.4 用VARCLUS过程进行变量聚类

8.1.1 聚类分析的一般概念 设有n个样品(多元观测值)，每个样品测得m项指标(变量)，得到观测数据xij（i=1,…,n；j=1,…,m），如表所示。 X1 X2 … Xm X(1) x11 x12 x1m X(2) x21 x22 x2m X(n) xn1 xn2 xnm

表8-1中数据又称为观测数据阵或简称为数据阵，其数学表示为： 其中列向量Xj = (x1j，x2j，…，xnj)'，表示第j项指标（j = 1，2，…，m），行向量X(i) = (xi1，xi2，…，xin)表示第i个样品。

1. 两种聚类分析 根据分类对象的不同，聚类分析分为两种： (1) 样品聚类：样品聚类是对样品（观测）进行的分类处理，又称为Q型分类，相当于对观测数据阵按行分类。 (2) 变量聚类：变量聚类是对变量（指标）进行的分类处理，又称为R型分类，相当于对观测数据阵按列分类。 两种聚类在形式上是对称的，处理方法也是相似的。

2. 聚类分析的方法 聚类方法大致可归纳如下： (1) 系统聚类法（谱系聚类） 先将l个元素（样品或变量）看成l类，然后将性质最接近（或相似程度最大）的2类合并为一个新类，得到l – 1类，再从中找出最接近的2类加以合并变成了l – 2类，如此下去，最后所有的元素全聚在一类之中。 (2) 分解法（最优分割法） 其程序与系统聚类相反。首先所有的元素均在一类，然后按照某种最优准则将它分成2类、3类，如此下去，一直分裂到所需的k类为止。

(3) 动态聚类法（逐步聚类法） 开始将l个元素粗糙地分成若干类，然后用某种最优准则进行调整，一次又一次地调整，直至不能调整为止。 (4) 有序样品的聚类 n个样品按某种因素（时间或年龄或地层深度等）排成次序，要求必须是次序相邻的样品才能聚在一类。 其他还有：有重叠聚类、模糊聚类、图论聚类等方法。

3. 聚类统计量 聚类分析实质上是寻找一种能客观反映元素之间亲疏关系的统计量，然后根据这种统计量把元素分成若干类。常用的聚类统计量有距离系数和相似系数两类。距离系数一般用于对样品分类，而相似系数一般用于对变量聚类。距离的定义很多，如马氏距离、明考斯基距离、兰氏距离、切比雪夫距离以及常见的欧氏距离： 等。相似系数有相关系数、夹角余弦、列联系数等。

当然，采用不同的分类方法会得到不同的分类结果，有时即使是同一种聚类方法，因距离的定义方法不同也会得到不同的分类结果。对任何观测数据都没有唯一“正确的”分类方法。实际应用中，常采用不同的分类方法对数据进行分类，可以提出多种分类意见，由实际工作者决定所需要的分类数和分类情况。

8.1.2 系统聚类法的基本思想和步骤 1. 系统聚类法的基本思想 下面以样品聚类为例介绍系统聚类法。 8.1.2 系统聚类法的基本思想和步骤 下面以样品聚类为例介绍系统聚类法。 1. 系统聚类法的基本思想 设有n个样品，每个样品测得m项指标（见表8-1）。系统聚类方法的基本思想是：首先定义样品间的距离（或相似系数）和类与类之间的距离。一开始将n个样品各自自成一类，这时类间的距离与样品间的距离是等价的；然后将距离最近的两类合并，并计算新类与其他类的类间距离，再按最小距离准则并类。这样每次减少一类，直到所有的样品都并成一类为止。这个并类过程可以用谱系聚类图形象地表达出来。

2. 系统聚类法的基本步骤 (1) 数据变换 为了便于比较或消除量纲的影响，在作聚类之前常常首先要对数据进行变换。变换的方法有中心化变换、标准化变换、极差标准化变换、极差正规化变换、对数变换等。 最常用的标准化变换为： 其中 ， ，j = 1，2，…，m。 变换后的数据，每个变量的样本均值为0，标准差为1，而且标准化变换后的数据{ }与变量的量纲无关。

(2) 计算n个样品两两间的距离 选择度量样品间距离的定义，计算n个样品两两间的距离，得样品间的距离矩阵D(0)。 (3) 聚类过程 首先n个样品各自构成一类，类的个数k = n：Gi = {X(i)}（i = 1，…，n），此时类间的距离就是样品间的距离（即D(1) = D(0)）。 令j = 2，…，n，执行如下并类过程： 1) 合并类间距离最小的两类为一新类（类间距离参见下文“系统聚类分析的方法”）。此时类的总个数k减少1类，即k = n – j + 1；

2) 计算新类与其他类的距离，得新的距离矩阵D(j)）。 若合并后类的总个数k仍大于1，重复1)和2)步，直到类的总个数为1时止。 (4) 画谱系聚类图 谱系图能明确清晰地描述各个样本点在不同层次上聚合分类的情况。 (5) 决定分类的个数及各类的成员

3. 系统聚类分析的方法 设有原始数据阵，如表8-1所示。G为在某一聚类水平上的类的个数，Ck是当前(水平G)的第k类，nk为Ck中的样品个数，DKL为第G水平的类CK和类CL之间的距离。 根据类间距离计算方法的不同，有11种不同的聚类方法： (1) 类平均法 类平均法(Average Linkage)用两类样品两两观测间距离的平均作为类间距离，即 类平均法是一种应用较广泛，聚类效果较好的方法。

(2) 重心法 重心法(Centroid Method)用两个类重心(均值)之间的（平方）欧氏距离定义类间距离，即 (3) 最长距离法 最长距离法(Complete Method)用两类观测间最远一对观测的距离定义类间距离，即 (4) 最短距离法 最短距离法(Single Linkage) 用两类观测间最近一对观测的距离定义类间距离，即

(5) Ward最小方差法（离差平方和法） Ward最小方差法(Ward's Mininum-Variance Method)也称Ward离差平方和法。类间距离定义为 Ward方法并类时总是使得并类导致的类内离差平方和增量最小。 其它的聚类方法还有最大似然法(EML)法、可变类平均法(Flexible-Beta Method)、McQuitty相似分析法(McQuitty's Similarity Analysis)、中间距离法(Median Method)、两阶段密度估计法(Two –Stage Density Linkage)等。 类平均法和Ward最小方差法使用最广泛。

4. 系统聚类类数的确定 (1) 由适当的阀值确定 选定某种聚类方法，按系统聚类的方法并类后，得到一张谱系聚类图，聚类图（又称谱系图）只反映样品间（或变量间）的亲疏关系，它本身并没有给出分类，需要给定一个临界相似尺度，用以分割谱系图而得到样品（或变量）的分类，如给定临界值（阀值）为d，那么，当样品间或已并类间距离小于d时，认为这些样品和类的关系密切，应该归属一类。

(2) 根据统计量确定分类个数 在SAS的CLUSTER过程中，提供一些统计量可以近似检验类个数如何选择更合适，用统计量决定类数的方法来自统计的方差分析思想，下面作一些介绍。 1) R2统计量 其中PG为分类数为G个类时的总类内离差平方和，T为所有样品或变量的总离差平方和。R2越大，说明分为G个类时每个类内的离差平方和都比较小，也就是分为G个类是合适的。但是，显然分类越多，每个类越小，R2越大，所以我们只能取G使得R2足够大，但G本身比较小，而且R2不再大幅度增加。

2) 半偏R2统计量 在把类CK和类CL合并为下一水平的类CM时，定义半偏相关 其中BKL=WM – (WK +WL)为合并类引起的类内离差平方和的增量，Wt为类Ct的类内离差平方和。半偏R2用于评价一次合并的效果，其值是上一步R2与该步R2的差值。其值越大，说明上一次合并的效果越好。

3) 伪F统计量 伪 伪F统计量评价分为G个类的效果。伪F统计量越大，表示分为G个类越合理。通常取伪F统计量较大而类数较小的聚类水平。 4) 伪t2统计量 用此统计量评价合并类CK和类CL的效果，该值大说明合并的两个类CK和CL是很分开的，不应合并这两个类，而应该取合并前的水平。

8.1.3 用CLUSTER过程和TREE过程进行系统聚类 PROC CLUSTER <选项列表>； VAR <聚类用变量>； COPY <复制变量>； ID <变量名>； RUN；

其中： 1) PROC CLUSTER语句为调用CLUSTERS过程的开始，其常用选项及功能见表8-2。 2) VAR语句指定用来聚类的数值型变量。如果缺省，则使用没有列在其他语句中的所有数值型变量。 3) COPY语句把指定的变量复制到OUTTREE = 的数据集中，以备后用。 4) ID语句中指定的变量用于区分聚类过程中的输出及OUTTREE数据集中的观测。

2. TREE过程 TREE过程可以把CLUSTER过程产生的OUTTREE = 数据集作为输入，画出聚类谱系图，并按照用户指定的聚类水平(类数)产生分类结果数据集。一般格式如下： PROC TREE <选项列表>； COPY <复制变量>； ID <变量>； RUN； 其中： 1) PROC TREE语句为调用TREE过程的开始，其常用选项及功能见表8-3。

2) COPY语句把输入数据集中的变量复制到输出数据集。 表8-3 PROC TREE语句的常用选项 2) COPY语句把输入数据集中的变量复制到输出数据集。 3) ID语句用于指定在输出树状图中的识别对象，ID变量可以是字符或数值变量。如果省略，TREE过程将使用变量_NAME_。 选项名称 功能说明 DATA = 数据集 指定从CLUSTER过程生成的OUTTREE数据集作为输入。 OUT = 数据集 指定包含最后分类结果(每一个观测属于哪一类，用一个CLUSTER变量区分)的输出数据集。 NCLUSTERS 指定最后把样本观测分为多少个类。 HORIZONTAL 横向画聚类谱系图。

3. 应用实例 【例8-1】表8-4是全国沿海10省市农民2004年支出情况的汇总资料，表中涉及生活消费支出情况的八个指标。 假定上述数据已经存放在数据集Mylib.nm10。 地区 食品 衣着 居住 家庭设备及服务 医疗保健 交通和通讯 文教、娱乐及服务 其他商品及服务 天津 1017.72 181.21 508.06 108.51 177.10 230.41 376.87 42.23 河北 780.09 127.06 340.88 80.42 115.97 176.60 182.56 31.33 辽宁 962.00 154.97 290.08 72.58 145.19 186.77 217.95 43.42 上海 2191.15 279.73 1446.15 344.41 424.55 720.37 805.55 116.94 江苏 1317.88 163.53 467.62 141.43 163.16 293.07 373.39 72.47 浙江 1838.57 258.58 798.88 242.09 326.12 496.86 597.96 100.05 福建 1408.54 159.60 430.14 154.43 136.40 306.06 313.09 107.32 山东 1000.13 139.18 365.97 110.12 155.85 221.93 298.23 97.85 广东 1581.68 116.83 494.89 131.17 153.17 350.27 314.34 98.44 广西 1047.58 64.42 311.48 64.87 83.64 140.13 178.83 37.65

为了进行系统聚类并产生帮助确定类数的统计量，使用如下代码： 试利用汇总资料对10个地区进行分类 (1) 使用CLUSTER过程 为了进行系统聚类并产生帮助确定类数的统计量，使用如下代码： proc cluster data = mylib.nm10 standard method = ward outtree = otree pseudo; copy group; run; 结果中首先给出数据相关系数矩阵特征值方面的信息。

结果的最后部分为聚类分析的完整过程(Cluster History)，如图所示，内容按列划分从左到右依次为： “NCL”为类别数量，表示新类别形成后类别的总数； “--Clusters Joined---”为合并的类别，指明这一步合并了哪两个类，有两列。其中OBxxx表示某一个原始样品，而CLxxx表示在某一个聚类水平上产生的类。 “FREQ”表示这次合并得到的类有多少个样品。 “SPRSQ”是半偏R2，“RSQ”是R2，“PSF”为伪F统计量，“PST2”为伪t2统计量，“Tie”指示距离最小的候选类对是否有多对，本例全无。

(2) 确定分类个数 1) R2统计量（列标题为RSQ）用于评价每次合并成NCL个类时的聚类效果。R2越大说明NCL个类越分开，故聚类的效果好。R2的值总是在0和1之间，而且R2的值总是随着分类个数NCL的减少而变小。通过查看R2值的变化，可以确定n个样品分为几类最为合适。本例中，分为3个类之前(NCL>3)的并类过程中R2的减少是逐渐的，改变不大；当分为3类时的R2 = 0.838，而下一次合并后分为2类时R2下降较多(R2 = 0.721)。这时通过分析R2统计量可得出分为3个类是较合适的。

查看R2变化的大小也可以由合并类时的半偏R2(列标题为SPRSQ)得到。半偏R2的值是上一步R2与该步R2的差值，故某步的半偏R2值越大，说明上一步合并的效果好。本例中半偏R2最大和次大分别为NCL = 1和2，说明根据半偏R2准则分为两个类或三个类是较合适的。 2) 伪F统计量(列标题为PSF)用于评价分为NCL个类的聚类效果。伪F值越大表示这些观测样品可显著地分为NCL个类。本例中伪F最大和次大依次为NCL = 2和5（局部），说明根据伪F准则分为两个类、五个类是较合适的。

4) 伪t2统计量用以评价此步合并类的效果。由该统计量的定义知伪t2值大表明上一次合并的两个类是很分开的，也就是上一次聚类的效果是好的。本例中伪t2最大和次大分别为NCL = 1和2，说明根据伪t2准则分为两个类或三个类是较合适的。 由此看出：R2准则支持分为两类和三类；伪F统计量支持分为两类和五类；伪t2统计量支持分为两类和三类。综合分析认为，用Ward法10个地区分为两类或三类较合适。

从图中可见，分为两类的结果为：G1 = {上海,浙江}，G2 = {天津,江苏,福建,广东,山东,河北,辽宁,广西}。 使用如下代码画出谱系图（如图所示）： proc tree data = otree horizontal; ID Group; run; 从图中可见，分为两类的结果为：G1 = {上海,浙江}，G2 = {天津,江苏,福建,广东,山东,河北,辽宁,广西}。 分为三类的结果为：G1 = {上海,浙江}，G2 = {天津,江苏,福建,广东,山东}，G3 = {河北,辽宁,广西}。

相仿地，可以使用类平均法、中间距离法、可变类平均法等方法。不同的聚类方法得到的聚类结果或多或少会有些差别，在实际应用中，应综合各种计算结果，提出合适的分类个数。

8.1.4 用VARCLUS过程进行变量聚类 1. VARCLUS过程 常用的VARCLUS语句格式为： 其中： PROC VARCLUS <选项列表>; VAR <变量列表>； RUN； 其中： 1) PROC VARCLUS语句为调用VARCLUS过程的开始，其常用选项及功能见表8-5。 2) VAR语句指定要分类的变量，如省略VAR语句，则使用没有在其他语句中列出的所有数值型变量。

2. VARCLUS过程变量聚类的步骤 如果没有为VARCLUS过程提供初始分类的情况，VARCLUS过程开始把所有变量看成一个类，然后它重复以下步骤： 1) 首先挑选一个将被分裂的类。通常这个被选中的类的类分量所解释的方差百分比最小（选项PRECENT=）或者同第二主成分有关的特征值为最大（选项MAXETGH=）。 2) 把选中的类分裂成两个类。首先计算前两个主成分，再进行斜交旋转，并把每个变量分配到旋转分量对应的类里，分配的原则是使变量与这个主成分的相关系数为最大。

3) 变量重新归类。通过多次反复循环，变量被重新分配到这些类里，使得由这些类分量所解释的方差为最大。 当每一类满足用户规定的准则时，VARCLUS过程停止。所谓准则，或是每个类分量所解释的方差的百分比，或是每一类的第二个特征值达到预设定的标准为止。如果没有规定准则，则当每个类只有一个特征值大于1时，VARCLUS过程停止。

3. 应用实例 【例8-2】为研究人脑老化的严重程度，有人测定了不同年龄的60名正常男性10项有关指标的数据，见表8-6。各变量的含义如下：AGE为年龄、TJ为图片记忆、SG为数字广度记忆、TS为图形顺序记忆、XX为心算位数、XS为心算时间、CK为规定时间内穿孔数、BJ为步距、JJ为步行时双下肢夹角、BS步速。 假定上述数据已经保存在逻辑库Mylib.rnlh中，试对这些指标作变量聚类分析。

(1) SAS程序 输入如下程序： OPTIONS PS=500; /*要求输出的结果中每页包括500行内容，可避免不必要的SAS标题反复出现。*/ PROC VARCLUS data = Mylib.rnlh; VAR age tj sg xx xs ts ck bj jj bs; RUN; 说明：过程语句中没有任何选择项，默认的聚类方法为主成分聚类法。过程步最终会聚成多少类，将由默认的临界值来决定，即当每个类只有一个特征根大于1时，VARCLUS过程停止。

(2) 结果分析 1) 第1步：如图8-4所示。这是用分解法思想进行斜交主成分聚类的第1步，把全部10个变量聚成一类，能解释的方差为4.146566，占总方差10的41.47％，第2特征值为1.4695，并预告这一类将被分裂。

2) 第2步将1类分裂成2类，分别含3个和7个变量，图8-5给出聚类概要。 图中表明各类中的Variation Explained（解释方差，即第一特征值）、Proportion Explained（解释方差占本类总方差的百分比）、Second Eigenvalue（类中的第2特征值）等。此时能解释的方差为5.403654，占总方差10的54.04％。

图8-6给出相关系数的平方。 其中第3列R-Squared With Own Cluster是指每个变量与所属类分量之间相关系数的平方R2，如：变量BJ在第1类中，它与第1类分量（相当于主成分分析中的第1主成分）之间的R2是0.8166，同理可理解该列中的其他相关系数的含义；

第4列R-squared with Next Closest是指每个变量与相邻类的类分量之间的相关系数的平方R2，如：BJ与第2类分量之间的相关系数平方R2为0.2121，该值越小，说明分类越合理。 第5列R-squared with 1-R**2 Ratio是由同一横行的数据求得： 如：(1 – 0.8166)/(1 – 0.2121) = 0.2327，此值越小，表明分类越合理。从此列可看出，很多比值较大，说明这10个变量分成2类是不太合适的。

图8-7给出从标准化变量预测类分量的标准回归系数。 若设C1、C2分别为第1和第2类分量，则有： C1 = 0.356608 BJ + 0.377719 JJ – 0.353237 BS C2 = – 0.281106 AGE + 0.278996 TJ + 0.239620 SG – 0.120883 XX + 0.1567777XS + 165498TS + 0.260054BS

图8-8给出类结构。 类结构相当于因子分析中的因子模型，即每个标准化变量可以表示成全部类分量的线性组合。如：AGE = – 0 .574313C1 – 0.806650C2。 类内相关（图8-9）就是类分量之间的相关系数阵。 这里预告第2类将被分裂。

3) 第3步。 首先给出将2类分裂成3类时的聚类概要。3个类分别含3个、4个和3个变量，此时能解释的方差为6.392112，占总方差10的63.92％，如图8-10所示。 接着给出各变量与类间的相关系数平方（图8-11）

标准回归系数（上图）以及类分量之间的相关系数阵（下图）。 此时已达到默认的停止分裂的临界值（即每个类中只有一个特征值大于1），停止分裂。

最后，给出整个聚类过程的汇总信息（图8-14）。 ① ② ③ ④ ⑤ ⑥ 第①列表示分成1类、2类或3类时分别能解释的总方差量； 第②列表示分成1类、2类或3类时分别能解释的方差占全部10个变量的总方差的百分比； 第③列表示分成1类、2类或3类时由1个类成分能解释的方差占全部10个变量的总方差的最小百分比 第④列为各类中最大的第2特征值；

第⑤列为各类中1个变量与其所在类的类分量的最小相关系数的平方R2； 第⑥列为各类中(1-R2)own/(1-R2)next的最大比值。参见前面输出结果及其解释。 说明：如果设定分类的个数，或是设定其他停止分裂的准则，则每次分裂按第2特征值最大选择分裂的类。

(3) 结论 结合专业知识发现：第1类中含BJ(步距)、JJ(步行时双下肢夹角)、BS(步速)这三个与走步有关的变量；第2类中含TJ(图片记忆)、XX(心算位数)、AGE(年龄)、CK(穿孔)这4个与视力、记忆和协调能力有关的变量；第3类中含SG(数字广度记忆)、XS(心算时间)、TS(图形顺序记忆)这三个与记忆、计算有关的变量。 分类的结果将有助于研究者对影响人脑老化本质的认识，为进一步从事这方面的研究提供了一些线索。 说明：对于一批给定的数据，究竟应聚成几类合适，没有统一的规则。可先将数据聚成类，然后结合专业知识和各类能解释总方差的百分比来权衡。

8.2 判别分析 8.2.1 判别分析的基本概念 8.2.2 SAS中作判别分析的过程 8.2.3 分析实例

8.2.1 判别分析的基本概念 1. 距离判别法 距离判别的基本思想是：样品和哪个总体的距离最近，就判断它属于哪个总体。 (1) 两总体情况 8.2.1 判别分析的基本概念 1. 距离判别法 距离判别的基本思想是：样品和哪个总体的距离最近，就判断它属于哪个总体。 (1) 两总体情况 设有两个总体G1和G2，若定义样品x到G1和G2的距离分别为d2(x,G1)和d2(x,G2)，则距离判别的判别规则是： 即未知所属总体的样品x离哪个总体较近，就判x属于哪个总体。

当总体G1和G2为正态总体，距离选用马氏距离，公式为： 当1 = 2 = 时，可以证明： d2(x,G1) – d2(x,G2) = – 2[x – (μ1+μ2)/2] ' –1 (μ1–μ2)

称W(x)为判别函数。由于它是x的线性函数，因此又称为线性判别函数。线性判别的应用最为广泛。 如果1与2不等，则判别函数W(x)为： 令 于是判别规则可表示为 称W(x)为判别函数。由于它是x的线性函数，因此又称为线性判别函数。线性判别的应用最为广泛。 如果1与2不等，则判别函数W(x)为： W(x) = d2(x，G1) – d2(x，G2) = 即W(x)是x的二次函数。 ，

(2) 多总体情况 设有m个总体：G1，G2，…，Gm（m > 2），它们的均值、协差阵分别为μi，i，(i = 1，2，…，m)。对任意给定的样品x，要判断它来自哪个总体。 多总体按距离最近的准则对x进行判别归类时，首先计算样品x到m个总体的马氏距离di2(x)（i = 1,2,…,m），然后进行比较，把x判归距离最小的那个总体。 即若dh2(x) = min{di2(x) | i = 1,2,…,m}，则x  Gh。 计算马氏距离di2(x)（i = 1,2,…,m）时，类似地可考虑1 = 2 = … = m或i不全相等的两种情况，并用样本统计量作为μi和i的估计.

2. 费歇尔判别（Fisher） 下面以两个总体为例说明费歇尔判别的思想。 设有两个总体G1、G2，其均值分别为1和2，协方差阵分别1和2，并假定1 = 2 = ，考虑线性组合：y = L'x。通过寻求合适的L向量，使得来自两个总体的数据间的距离较大，而来自同一个总体数据间的差异较小。为此，可以证明，当选L=c–1(1–2)，其中c  0时，所得的投影即满足要求。从而称c = 1时的线性函数： y = L'x = (1 – 2)'–1x 为费歇尔线性判别函数。

其判别规则为： 其中，m为两个总体均值在投影方向上的中点，即 当1、2和未知时，可由总体G1和G2中分别抽出n1和n2个样品，计算相应的样本均值和协方差阵作为1、2和的估计。

3. 贝叶斯判别（Bayes） (1) 贝叶斯判别法 设有m个总体G1,…,Gm，假定它们各自的分布密度分别为f1(x)，f2(x)，…，fm(x)，各自的先验概率（可以根据经验事先给出或估出）分别为q1,q2,…,qm，显然 贝叶斯判别的方法是：当抽取了一个未知总体的样品x，要判断它属于哪个总体，可用著名的贝叶斯公式计算x属于第j个总体的后验概率：

当 时，判断x属于第h个总体。或者计算按先验概率加权的误判平均损失： 然后再比较这m个误判平均损失的h1(x)，h2(x)，…，hm(x)的大小，选取其中最小的，就可以判定样品x来自该总体。 上式中C(j | i)为假定本来属于Gi的样品被判为属于Gj时造成的损失。当然C(i | i) = 0，C(j | i) 0（i，j = 1，…，m）。

在实际问题中，错判的损失可以给出定性的分析，但很难用数值来表示，但应用贝叶斯判别的方法时，要求定量给出C(j | i)，C(j | i)的赋值常用以下两种方法： ● 根据经验人为赋值； ● 假定各种错判的损失都相等。

(2) 错判概率 当样品xGi，用判别法D判别时，把x判归Gj (i≠j)，出现错判。用P(j|I;D)（或简记为P(j|i)）表示实属Gi的样品错判为Gj的概率，错判概率的估计方法有下面几种： ● 利用训练样本作为检验集，即用判别方法对已知样品进行回判，统计错判的个数以及错判的比率，作为错判率的估计。此法得出的估计一般较低。 ● 当训练样本的大小足够大时，可留出一些已知类别的样品不参加建立判别准则，而是作为检验集，并把错判的比率作为错判率的估计。此法当检验集较小时估计的方差大。

● 设一法（或称为交叉确认法），每次留出一个已知类别的样品，而用其他n – 1个样品建立判别准则，然后对留出的这一个已知类别的样品进行判别归类。对训练样本中n个样品逐个处理后把错判的比率作为错判率的估计。 (3) 广义平方距离判别法 在正态总体的假定下，按贝叶斯判别的思想，在错判造成的损失认为相等的情况下得到的判别函数其实就是马氏距离判别在考虑先验概率及协方差阵是否相等情况下的推广，故在SAS的DISCRIM过程中称为广义平方距离判别法。

4. 逐步判别 (1) 各变量判别能力的检验 当检验k个类的均值向量是否全都相等(即检验H0：1 = 2 = … = k)时，否定了这个假设H0(即表明各总体的均值向量有显著性差异)，也并不能保证其各分量的均值有显著差异，若第i个分量间没有显著差异时，说明相应的变量Xi对判别分类不起作用，应该剔除。关于各变量判别能力的检验问题是筛选判别变量的理论基础，也是逐步判别的理论基础。

变量判别能力的度量通常采用删去该变量后考察判别能力的变化，即考察该变量对区分k个类是否能提供更多的附加信息，然后由附加信息构造F统计量进行检验。 利用F统计量对假设H0(i)(第i个变量在k个总体中的均值相等)作统计检验。若否定H0(i)，表示变量Xi对区分k个总体的判别能力是显著的(在显著水平α下)。否则，变量Xi对区分k个总体的判别能力不能提供附加信息，这个变量应剔除。

(2) 逐步判别法的基本思想 前面讨论了用全部m个变量：X1,X2,…,Xm来建立判别函数对样品进行判别归类的几种方法。在这m个变量中，有的变量对区分k个总体的判别能力可能很强，有的可能很弱。如果不加区别地把m个变量全部用来建立判别函数，则必增加大量的计算，还可能因为变量间的相关性引起计算上的困难(病态或退化等)及计算精度的降低。另一方面由于一些对区分k个总体的判别能力很小的变量的引入，产生干扰，致使建立的判别函数不稳定，反而影响判别效果，因此自然提出一个变量的选择问题。即如何从m个变量中挑选出对区分k个总体有显著判别能力的变量，来建立判别函数，用以判别归类。

逐步判别的基本思想和逐步回归是类似的。逐个引入变量，每次把一个判别能力最强的变量引入判别式，每引入一个新变量，对判别式中的老变量逐个进行检验，如其判别能力因新变量的引入而变得不显著，应把它从判别式中剔除。这种通过逐步筛选变量使得建立的判别函数中仅保留判别能力显著的变量的方法，就是逐步判别法。

(3) 逐步判别法的基本步骤 1) 逐步筛选变量：根据各变量对区分k个总体的判别能力的大小，按基本思想所介绍的过程来筛选变量。SAS中的STEPDISC过程专用于筛选变量子集。该过程利用向前选入、向后剔除或逐步筛选的方法来选择区分k个总体的最佳变量子集。 2) 判别归类：对已选出变量子集，使用以上几节介绍的判别方法对样品进行判别归类。

8.2.2 SAS中作判别分析的过程 1. DISCRIM过程简介

如果假设每个组内样本的分布为多元正态分布，基于多元正态分布理论的参数法将导出一个线性或二次的判别函数；否则，将采用不基于任何分布的假定的非参数方法，来估计类别的密度从而实现分类的功能。此类非参数检验的方法包括核密度法(Kernel Method)和K最邻近法(K-Nearest-Neighbor-Method)，其中核密度法又可因核的确定方法不同而分为多种不同的判别方法。

其中PROC DISCRIM语句和CLASS语句是DISCRIM过程的必选语句，其它均为可选项。 CLASS <变量>；BY <变量列表>； ID <变量>； PRIORS <先验概率值>； TESTCLASS <变量>； TESTID <变量>； VAR <变量列表>； RUN; 其中PROC DISCRIM语句和CLASS语句是DISCRIM过程的必选语句，其它均为可选项。

(2) PROC DISCRIM语句 PROC DISCRIM语句为调用DISCRIM过程的开始，语句中可设置的选项及其功能见表8-7。 (3) CLASS语句 该语句在PROC DISCRIM中不可省略，分类变量的不同值决定判别分析的组。类水平由CLASS变量的格式化值决定，指定的变量可以是数值的或字符的。

(4) PRIORS语句 指定各组中成员出现的先验概率。若觉得各类别可能出现的概率相等，则使用PRIORS EQUAL语句。若规定先验概率为各组样本出现的比例，则使用PRIORS PROPORTIONAL(或PRIORS PROP语句)。若知道各个组的先验概率，则可以在此指令中列出各类别的概率。 注：概率的总和必须是1。 (5) VAR语句 列出所有判别分析可能用到的变量名称。若此语句省略，则数据集中未在其他语句中使用的全部数值型变量被作为VAR变量来使用。

2. STEPDISC过程简介 针对具有一个分类变量和若干数值型变量(指标变量)的数据集，STEPDISC过程执行逐步判别分析(stepwise discriminant analysis):从指定的指标变量(VAR变量)中筛选出一组变量，以用于随后的判别分析。逐步判别分析要求指标变量在各组内服从多元正态分布，并且具有相同的协方差矩阵。STEPDISC过程选择变量的方法包括前进法(forward selection)、后退法(backward elimination)以及逐步法(stepwise selection)等。

在使用DISCRIM过程进行判别分析之前，一般情况下均需先调用STEPDISC过程对指标变量进行筛选，以避免无关变量对判别分析结果的影响或指标变量过多导致判别结果的不稳定。

其中PROC STEPDISC语句和CLASS语句为STEPDISC过程运行所必需的语句，其余语句均为可选项。 BY <变量列表>； VAR <变量列表>； RUN； 其中PROC STEPDISC语句和CLASS语句为STEPDISC过程运行所必需的语句，其余语句均为可选项。

(2) PROC STEPDISC语句 PROC STEPDISC语句为调用STEPDISC过程的开始，语句中可设置的选项及其功能见表8-8。 (3) CLASS语句 CLASS语句为STEPDISC过程所必需的语句，其语法与前述的DISCRIM过程的同名语句完全相同。 (4) VAR语句 指定筛选变量的范围，需要考察的变量必须全部在VAR语句中指定，且须为数值型变量。如果忽略VAR语句，则数据集中未在其他语句中使用的全部数值型变量将被作为VAR变量来使用。

8.2.3 分析实例 【例8-3】中小企业的破产模型。为了研究中小企业的破产模型，选定4个经济指标：x1总负债率（现金收益/总负债）、x2收益性指标（纯收入/总财产）、x3短期支付能力（流动资产/流动负债）、x4生产效率性指标（流动资产/纯销售额）。 对17个破产企业（1类）和21个正常运行企业（2类）进行了调查，调查资料见表8-9。 假定上述数据已经保存在逻辑库Mylib.pcmx中，试对这些指标作判别分析。

1. 使用DISCRIM过程作距离判别 (1) SAS程序 OPTIONS PS=500; proc discrim data = Mylib.PCMX distance listerr; class Type; var x1-x4; run; 说明：其中distance选项表示可根据类间距离来判断各类样品有无显著差异，从而说明判别归类问题有无意义。

(2) 结果分析 DISCRIM过程的分析结果首先给出观测（样品）个数、变量个数、分类个数以及自由度（如图所示）。 “Class Level Information”为各类别信息，依次为：分类变量水平值（Type）、哑变量名（Variable Name）、各类别观测数、各类别权重、各类别占总样本数的比例、各类别的先验概率（Prior Probability）。 其中先验概率的默认值为0.5。

“Squared Distance to Type”（图8-16左）给出组间距离平方：组间的距离平方及其F检验值、F检验的p值。 结果表明：两组间的距离（马氏距离）平方为：3.17425，检验H0：μ1 = μ2的F统计量为6.83404，相应的p值为0.0004 < 0.05，这说明两组企业的模型（破产）特征有显著差异，讨论判别归类问题是有意义的。

“Generalized Squared Distance to Type”(图左)给出的广义组间距离平方与组间的马氏距离平方相同。 “Linear Discriminant Function for Type”(图右)给出线性判别函数： W1(x) = –3.04357 – 3.11878 x1 – 0.44945 x2 + 1.58110 x3 + 9.05966 x4 W2(x) = –6.05818 + 2.22679 x1 – 2.75270 x2 +2.91289 x3 + 8.71392 x4

“Posterior Probability of Membership in Type” 是由选项Listerr要求列出的回判结果：组中成员的后验概率，此部分为判别分类中被误判的样品及其有关信息。由于代码中忽略了ID语句，此处以观测号对样品进行标识。 各列从左到右依次为观测号、原来所属的类别、被归入的新类别以及各样品依据判别函数所得的归属各类别的后验概率，DISCRIM过程依据此后验概率对每个样品进行判别归类

图8-19为已知类别的样品被判别函数归入各类别的频数和百分比。 图8-20为各类别的错误分类率(Error Count Estimates for type)，可见第二类(Type = 2)的错分率较高，为0.2381 = 5/21，各类别的总错分率为0.1779 = 0.1176*0.5 + 0.2381*0.5。

2. 使用DISCRIM过程作贝叶斯判别 (1) SAS程序 proc discrim data = Mylib.PCMX crosslist; class Type; var x1-x4; run; 说明：用贝叶斯判别法（广义平方距离判别法）建立判别准则，并对各样品进行判别归类。其中crosslist要求输出交叉验证后的分类结果。

(2) 结果分析 除了给出一般的结果之外，还有如下输出。 1) 交叉验证判别分类的样品及其有关的信息。各列从左至右依次为样品号、原来所属的类别（From Type）、被归入的新类别（Classified into Type）以及对各样品依据判别函数所得的归属各类别的后验概率（Posterior Probability of Membership），如图所示。

2) 来自已知类别的样品被判别函数归入各类别的频数和百分比，以及各类别的错误分类率(Error Count Estimates for type)，如图8-22所示。 与训练样本回代的判别分类结果相比，交叉验证的错分率明显增大，但交叉验证的分类结果更加真实和客观，在没有检验数据集时交叉验证的结果是评价判别函数之判别效能的可靠指标。

3. 使用DISCRIM过程作非参数判别 (1) SAS程序 修改上述代码如下： proc discrim data = Mylib.PCMX Method = npar k = 3 crosslisterr; class Type; var x1-x4; run; 说明：Method = npar指定采用非参数判别法，k = 3指定采用非参数判别法中的K最邻近法，且近邻个数为3（可经过多次尝试，得到最佳的判别效果）。

(2) 结果分析 因篇幅所限，此处仅给出交叉验证的分类结果。 从结果可见非参数方法的判别效果比参数法大大提高，错分率估计值为0.1891。

4. 综合实例分析 【例8-4】为了明确诊断小儿肺炎的三种类型，某研究单位测得32名结核性、13名化脓性和18名细菌性肺炎患儿的7项生理、生化指标，见表8-10。试建立判别函数并进行判别归类。 将表8-10中的数据存为数据集Mylib.fylx，然后编写SAS程序如下(分为变量筛选和判别分析两个部分)。

(1) 变量筛选SAS程序 Proc stepdisc slentry = 0.1 slstay = 0.1; class type; var x1 - x7; RUN; PROC STEPDISC语句中，选项“slentry = 0.1”和“slstay = 0.1”设定变量进入模型和变量留在模型中的显著性水平均为0.1，CLASS语句指定变量“type”为类别标志的变量，VAR语句指定从变量x1～x7中选择构建判别模型的指标变量。此处未设置“method =”选项，STEPDISC过程将以默认的变量选择方法(即逐步法)进行变量筛选。

将以上程序提交运行，结果如下。逐步判别分析在默认情况下的输出结果篇幅较大，其中包含了变量筛选的每一个步骤的分析结果。为全面了解逐步判别分析过程，此处将输出结果全部列出，并对每一部分进行简要说明。 首先给出的是有关数据集的结构情况和逐步判别分析的参数设置情况，如图8-24所示。

数据结构：样品数63，指标变量数7，分类变量水平数3。参数设置：变量筛选方法为逐步法，变量进入模型和变量留在模型中的显著性水平均为0.1。 标题为“Class Level Information”的结果列表列出与各类别有关的信息，各列从左到右依次为分类变量水平值(type)、哑变量名(Variable Name，STEPDISC过程运行时以此哑变量进行运算)、各类别样品数(Frequency)、各类别权重(Weight)以及各类别占总样品数的比例(Proportion)。

图8-25给出变量筛选第一步的结果。 首先对模型外的变量进行协方差分析(此处因为模型中尚无已选入的变量，协方差分析即为方差分析)，标题为“Statistics for Entry，DF：2，60 (F检验的自由度为：2，60)”的部分即为分析的结果。从结果中的F值(F Value)一栏，可判断x1与分类变量的关系最为密切，相应的p值也最小并达到变量进入模型的标准(p < 0.10)。从复相关系数(R Square)也可以做出同样的判断。因此，变量x1被选入判别模型。

“Multivariate statistics”为对现有模型判别效能大小进行估计的统计量(本质上为对判别模型的多元方差分析)，包括Wilks'值和Pillai's值等。

图8-26为变量筛选第二步到第七步的结果。基本内容和第一步没有太多的区别。 第二步中首先对模型内的变量分别进行协方差分析(此处模型内仅有一个变量，所作分析为方差分析)，标题为“Statistics for Removal，DF = 2，60”的部分即为分析结果。结果可见x1的R2为0.483，F值为28.03，p < 0.0001，符合留在模型中的标准。

…

然后再对模型外剩余的6个变量进行协方差分析，标题为"Statistics for Entry，DF = 2，59”的部分即为分析的结果。从F值一栏，可判断x5与分类变量的关系最为密切，相应的p值也最小并达到变量进入模型的标准(p < 0.10)。变量x5被选入判别模型。

标题为“Multivariate statistics”的部分为对新模型判别效能大小进行估计的统计量，可见与第一步所建立的模型相比，Wilks'值和Pillai's迹均有所变化，判别效能增加。 变量筛选第三步到第六步，STEPDISC过程分别对模型中和模型外的变量进行相应的分析，判断其是否符合剔除或进入的标准，以此进行变量的剔除和选入。此过程中变量x6，x7，x2，x3被分别选入模型，没有变量被剔除。

变量筛选第七步，从协方差分析结果中可见x1，x2，x3，x5，x6，x7均符合留在模型中的标准，因此没有变量被从模型中剔除。对模型外剩余的1个变量x4进行协方差分析，判断x4与分类变量的关系不够密切，相应的p值达不到进入模型的标准，变量x4不能选入判别模型。 此时模型中没有变量可剔除，模型外也没有变量可进入，变量选择过程停止。 结果的最后部分(标题为“Stepwise Selection Summary”的结果列表)是对变量选择过程和结果的汇总，如图8-27所示。

各列从左到右依次为步骤序号，模型中所含变量数，每一步所选入模型的变量，每一步从模型中剔除的变量，所考查变量的平方偏相关系数，对所考查变量进行协方差分析的F值及相应p值，每一步完成后模型的Wilks'值及相应的p值，平方典则相关的平均及相应的p值。

(2) 判别分析SAS程序 Proc discrim data = Mylib.fylx pool = test crosslisterr; Class type; Var x1 - x3 x5 - x7; Run; 在PROC DISCRIM语句中，“data = Mylib.fylx”用来指定判别分析的训练数据集，“pool = test”要求在决定是否用合并协方差矩阵进行广义平方距离的计算以前，先对各类别内协方差矩阵进行齐性检验，根据检验结果决定所用的协方差矩阵。“crosslisterr”要求在输出结果中给出交叉验证结果中被错误分类的样品，此选项隐含了对“crossvalidate”的设置，要求DISCRIM过程对判别函数进行交叉验证。

程序中未设置“method=”选项，DISCRIM过程采用默认的设置即参数法进行判别分析的过程。程序中忽略了priors语句，因此各类别先验概率值将被设置为相等。 提交上述程序，分析结果如下。

图8-28给出的是有关数据集的结构情况。 数据一般情况：样品数63，指标变量数6，分类变量水平数3，总自由度(DF Total)62，类内自由度(DF Within Class)60，类间自由度(DF Between Class)2。

“Class Level Information”：各列从左到右依次为Type、Variable Name、Frequency、Weight以及Proportion、Prior Probability等，由于先验概率设置为各类别相等，因此各类的先验概率均为0.333333。 “Within Covariance Matrix Information”：各列从左到右依次为Type、Covariance Matrix Rank、Natural Log of the Determinant of the Covariance Matrix。此部分的最后一行为合并类内协方差矩阵的有关信息。

图为各类内协方差矩阵齐性检验的结果(Test of Homogeneity of Within Covariance Matrices)，其中的大部分内容是有关符号的定义及有关指标计算方法的介绍。 此处2检验的结果表明各类内协方差不齐(p < 0.0001)，因此各类内的协方差矩阵将被用于广义平方距离的计算，从而将得到一个二次判别函数(quadratic discriminant function)而非线性判别函数(linear discriminant function)。

图8-30为各类间的广义平方距离，可见类内平方距离小于类间平方距离。

图8-31所示为训练样本回代的判别分类汇总结果，上半部分为来自各类别的样品被判别函数归入各类别的频数和百分比，下半部分为各类别的错误分类率，可见第三类的错分率最高，为0.1667，各类别的总错分率为0.1124。

图8-32为交叉验证判别分类中被错误分类的样品及其有关的信息。各列从左到右依次为样品号、原来所属的类别、被归入的新类别，后三列为对各条样品依据判别函数所得的归属各类别的后验概率，discrim过程依据此后验概率对每条样品进行判别归类。

最后部分（图8-33）给出的是交叉验证的判别分类汇总结果，上半部分为来自各类别的样品被判别函数归入各类别的频数和百分比，下半部分为各类别的错误分类率，可见第三类的错分率最高，为0.3889，各类别的总错分率为0.3516。与训练样本回代的判别分类结果相比，交叉验证的错分率明显增大，但交叉验证的分类结果更加真实和客观，在没有检验数据集时交叉验证的结果是评价判别函数之判别效能的可靠指标。

由于判别函数为二次方程，输出结果中没有关于判别函数的内容，只有判别函数为线性方程时，输出结果中才会给出有关的模型和相应的参数。

将程序(2)中PROC DISCRIM语句修改如下： (3) 非参数判别法 从本例结果来看，所建立判别函数的判别效能不太理想，错分率估计值达到0.3516，这可能与数据不符合参数法判别分析的条件(各变量在各类别内服从多元正态分布)有关，应采用的判别分析方法为非参数法。 将程序(2)中PROC DISCRIM语句修改如下： Proc discrim data = Mylib.fylx method = npar k = 3 crosslisterr; “method=npar”指定采用非参数方法进行判别分析，“k = 3”指定采用非参数方法中的K最近邻法，且近邻个数为3(实际操作中此数值需经过多次尝试，方可获得最佳的判别效能)。

从结果可见非参数方法的判别效能比参数法大大提高，错分率估计值为0.2064。 修改后的代码如下： Proc discrim data = Mylib.fylx method = npar k = 3 crosslisterr; Class type; Var x1 - x3 x5 - x7; Run; 程序(2)提交得结果如图。 从结果可见非参数方法的判别效能比参数法大大提高，错分率估计值为0.2064。