应用多元统计分析讲稿(朱建平)

第一节 引言

多元统计分析是运用数理统计方法来研究解决多指标问题的理论和方法。近30年来，随着计算机应用技术的发展和科研生产的迫切需要，多元统计分析技术被广泛地应用于地质、气象、水文、医学、工业、农业和经济等许多领域，已经成为解决实际问题的有效方法。然而，随着Internet的日益普及，各行各业都开始采用计算机及相应的信息技术进行管理和决策，这使得各企事业单位生成、收集、存储和处理数据的能力大大提高，数据量与日俱增，大量复杂信息层出不穷。在信息爆炸的今天，人们已经意识到数据最值钱的时代已经到来。

显然，大量信息在给人们带来方便的同时也带来一系列问题。比如：信息量过大，超过了人们掌握、消化的能力；一些信息真伪难辩，从而给信息的正确应用带来困难；信息组织形式的不一致性导致难以对信息进行有效统一处理等等，这种变化使传统的数据库技术和数据处理手段已经不能满足要求.Internet的迅猛发展也使得网络上的各种资源信息异常丰富，在其中进行信息的查找真如大海捞针。这样又给多元统计分析理论的发展和方法的应用提出了新的挑战。

多元统计分析起源于上世纪初，1928年Wishart发表论文《多元正态总体样本协差阵的精确分布》，可以说是多元分析的开端。20世纪30年代R.A. Fisher 、H.Hotelling、S.N.Roy、许宝騄等人作了一系列得奠基性工作，使多元分析在理论上得到了迅速得发展。20世纪40年代在心理、教育、生物等方面有不少得应用，但由于计算量大，使其发展受到影响，甚至停滞了相当长得时间。20世纪50年代中期，随着电子计算机得出现和发展，使多元分析方法在地质、气象、医学、社会学等方面得到广泛得应用。20世纪60年代通过应用和实践又完善和发展了理论，由于新的理论、新的方法不断涌现又促使它的应用范围更加扩大。20世纪70年代初期在我国才受到各个领域的极大关注，并在多元统计分析的理论研究和应用上也取得了很多显著成绩，有些研究工作已达到国际水平，并已形成一支科技队伍，活跃在各条战线上。

在20世纪末与本世纪初，人们获得的数据正以前所未有的速度急剧增加，产生了很多超大型数据库，遍及超级市场销售、银行存款、天文学、粒子物理、化学、医学以及政府统计等领域，多元统计与人工智能和数据库技术相结合，已在经济、商业、金融、天文等行业得到了成功的应用。

为了让人们更好的较为系统地掌握多元统计分析的理论与方法，本书重点介绍多元正态总体的参数估计和假设检验以及常用的统计方法。这些方法包括判别分析、聚类分析、主成分分析、因子分析、对应分析、典型相关分析、多维标度法以及多变量的可视化分析等。与此同时，我们将利用在我国广泛流行的SPSS统计软件来实现实证分析，做到在理论的学习中体会应用，在应用的分析中加深理论。

第二节 应用背景

二、多元统计分析方法的应用

这里我们要通过一些实际的问题，解释选择统计方法和研究目的之间的关系，这些问题以及本书中的大量案例能够使得读者对多元统计分析方法在各个领域中的广泛应用有一定的了解。多元分析方法从研究问题的角度可以分为不同的类，相应有具体解决问题的方法，参看表1.1。

多元统计分析方法在经济管理、农业、医学、教育学、体育科学、生态学、地质学、社会学、考古学、环境保护、军事科学、文学等方面都有广泛的应用，这里我们例举一些实际问题，进一步了解多元统计分析的应用领域，让读者从感性上加深对多元统计分析的认识。

问题 内容 方法 尽可能简单地表示所研究的现多元回归分析、聚类分析、象，但不损失很多有用的信息，主成分分析、因子分析、相数据或结构性化简 并希望这种表示能够很容易的应分析、多维标度法、可视解释。 化分析 分类和组合 基于所测量到的一些特征，给出判别分析、聚类分析、主成好的分组方法，对相似的对象或分分析、可视化分析 变量分组。 变量之间是否存在相关关系，相多元回归、典型相关、主成关关系又是怎样体现。 分分析、因子分析、相应分变量之间的相关关系 析、多维标度法、可视化分析 预测与决策 通过统计模型或最优准则，对未多元回归、判别分析、聚类来进行预见或判断。 分析、可视化分析 检验由多元总体参数表示的某多元总体参数估计、假设检假设的提出及检验 种统计假设，能够证实某种假设验 条件的合理性。 1、城镇居民消费水平通常用八项指标来描述，如人均粮食支出、人均副食支出、人均烟酒茶支出、人均衣着商品支出、人均日用品支出、人均燃料支出、人均非商品支出。这八项指标存在一定的线性关系。为了研究城镇居民的消费结构，需要将相关强的指标归并到一起，这实际就是对指标进行聚类分析。

2、在企业经济效益的评价中，涉及到的指标往往很多，如百元固定资产原值实现产值、百元固定资产原值实现利税、百元资金实现利税、百元工业总产值实现利税、百元销售收入实现利税、每吨标准煤实现工业产值、每千瓦时电力实现工业产值、全员劳动生产率、百元流动资金实现产值。如何将这些具有错综复杂关系的指标综合成几个较少的因子，既有利于对问题进行分析和解释，又能便于抓住主要矛盾做出科学的评价。可用主成分分析和因子分析法。

3、某一产品是用两种不同原料生产的，试问此两种原料生产的产品寿命有无显著差异？又比如，若考察某商业行业今年和去年的经营状况，这时需要看这两年经营指标的平均水平是否有显著差异以及经营指标之间的波动是否有显著差异。可用多元正态总体均值向量和协差阵的假设检验。

4、按现行统计报表制度，农村家庭纯收入是指农村常住居民家庭总收入中扣除从事生产和非生产经营用支出、税款和上交承包集体任务金额以后剩余的、可直接用于进行生产的、非生产性建设投资、生产性消费的那一部分收入。如果我们收集某年各个省、自治区、直辖市农民家庭人均纯收入的数据，可以用相应分析，揭示全国农民人均纯收入的特征以及各省、自治区、直辖市与各收入指标的关系。

5、某医院已有100个分别患有胃炎、肝炎、冠心病、糖尿病等的病人资料，记录了他们每个人若干项症状指标数据。如果对于一个新的病人，当也测得这若干项症状指标时，可以利用判别分析方法判定他患的是哪种病。

6、有100种酒，品尝家可以对每两种酒进行品尝对比，给出一种相近程度的得分（越相近得分越高，相差越远得分越低），希望用这些得分数据来了解这100种酒之间的结构关系。这样的问题就可以用多维标度法来解决。

7、在地质学中，常常要研究矿石中所含化学成分之间的关系。设在某矿体中采集了60个标本，对每个标本测得20个化学成分的含量。我们希望通过对这20个化学成分的分析，了解矿体的性质和矿体形成的主要原因。

8、对1000个类似的鱼类样本，如何根据测量的特征如体重、身长、鳍数、鳍长、头宽等，我们可以利用聚类分析方法将这类鱼分成几个不同品种。

9、考古学家对挖掘出来的人头盖骨的高、宽等特征来判断是男或女，根据挖掘出的动物牙齿的有关测试指标，判别它是属于哪一类动物牙齿、是哪一个时代的。

10、在高考招生工作中，我们知道每个考生的基本情况，通过分析我们不仅可以了解到学生喜欢学习的科目，还可以进一步从考生每门课程的成绩，分析出学生的逻辑思维能力、形象思维能力和记忆力等等对学习成绩的影响。

第二章 多元正态分布的参数估计

第一节 引言

多元统计分析涉及到的都是随机向量或多个随机向量放在一起组成的随机矩阵。例如在研究公司的运营情况时，要考虑公司的获利能力、资金周转能力、竞争能力以及偿债能力等财务指标；又如在研究国家财政收入时，税收收入、企业收入、债务收入、国家能源交通重点建设基金收入、基本建设贷款归还收入、国家预算调节基金收入、其他收入等都是需要同时考察的指标。显然，如果我们只研究一个指标或是将这些指标割裂开分别研究，是不能从整体上把握研究问题的实质的，解决这些问题就需要多元统计分析方法。为了更好的探讨这些问题，本章我们首先论述有关随机向量的基本概念和性质。

在实用中遇到的随机向量常常是服从正态分布或近似正态分布，或虽本身不是正态分布，但它的样本均值近似于正态分布。因此现实世界中许多实际问题的解决办法都是以总体服从正态分布或近似正态分布为前提的。在多元统计分析中， 多元正态分布占有很重要地位，本书所介绍的方法大都假定数据来之多元正态分布。为此，本章将要介绍多元正态分布的定义和有关性质。

然而在实际问题中，多元正态分布中均值向量和协差阵通常是未知的，一般的做法是由样本来估计。这是本章讨论的重要内容之一，在此我们介绍最常见的最大似然估计法对参数进行估计，并讨论其有关的性质。

第二节 基本概念

一、随机向量

我们所讨论的是多个变量的总体，所研究的数据是同时p个指标（变量），又进行了n次观测得到的，我们把这个p指标表示为X1 ，X2，L，Xp，常用向量X = (X1 ， X2 ， L ， XP)' 表示对同一个体观测的p个变量。这里我们应该强调，在多元统计分析中，仍然将所研究对象的全体称为总体，它是由许多（有限和无限）的个体构成的集合，如果构成总体的个体是具有p个需要观测指标的个体，我们称这样的总体为p维总体（或p元总体）。上面的表示便于人们用数学方法去研究p维总体的特性。这里“维”（或“元”）的概念，表示共有几个分量。若观测了n个个体，则可得到如表2.1的数据，称每一个个体的p个变量为一个样品，而全体n个样品组成一个样本。

表2.1 数据 变量 序号 1 2

X1 X2 Xp X11 X12 X22

X1p X2p

X21

n

在这里横看表2.1，记为 X()(X1,X2, Xj(X1j,X2j,Xn1 Xn2

Xnp

,Xp )， 1,2,n,

表示第个样品的观测值。竖看表2.1，第j列的元素

,Xnj) ， j1,2,p,

表示对第j个变量Xj的n次观测数值。

因此，表2.1所反映出的样本资料可用矩阵表示为

X11X21XXn1X12X22Xn2X1pX(1)XX2p(X,X，，X)(2) （2.1）

12pXnpX(n),Xp的整体称为p维随机向量，记为X(X1,X2,,Xp)。

简记为X。

定义2.1 将p个随机变量X1,X2,在对随机向量的研究仍然限于讨论离散型和连续型两类随机向量。 二、多元分布

先回顾一下一元统计中分布函数和密度函数的定义。

设X是一个随机变量，称F(x)p(Xx)为X的概率分布函数或简称为分布函数，记为X~F(x)。 若随机变量在有限或可列个值xk上取值，记P(Xxk)pk，(k1,2,散型随机变量，称P(Xxk)pk，(k1,2,)且pk1，则称X为离

k)为X的概率分布。

x设X~F(x)，若存在一个非负函数f(x)，使得一切实数x有：F(x)f(t)dt，则称f(x)为X的分布密度函数，简称为密度函数。一个函数f(x)能作为某个随机变量X的分布密度函数的重要条件是： （1）f(x)0，对一切实数x；

（2）

f(x)dx1。

,Xp)是p维随机向量，它的多元分布函数定义为

F(x)F(X1,X2,,Xp)P(X1x1,X2x2,,Xpxp) （2.2）

定义2.2 设X(X1,X2,记为X~F(x)，其中x(x1,x2,,xp)Rp，Rp表示p维欧氏空间。

多维随机向量的统计特性可用它的分布函数来完整地描述。

定义2.3 设X(X1,X2,,Xp)是p维随机向量，若存在有限个或可列个p维数向量x1,x2,,，记

P(Xxk)pk，(k1,2,)且满足p1p21，则称X为离散型随机向量，称P(Xxk)pk，(k1,2,)为X的概率分布。

设X~F(x)F(x1,x2,,xp)，若存在一个非负函数f(x1,x2,,xp)，使得对一切

x1xpx(x1,x2,,xp)Rp有F(x)F(x1,x2,,xp)f(t1,t2,,tp)dt1dtp（2.3）

则称X为连续型随机变量，称f(x1,x2,,xp)为分布密度函数，简称为密度函数或分布密度。

p一个p元函数f(x1,x2,,xp)能作为R中某个随机向量的密度函数的主要条件是：

p（1）f(x1,x2,,xp)0，(x1,x2,,xp)R；

12p（2）f(x,x,,x)dx1dxp1

离散型随机向量的统计性质可由它的概率分布完全确定，连续型随机向量的统计性质可由它的分布密度完全确定。

【例2.1】 试证函数

e(x1x2),x10,x20 f(x1,x2)

其它0,为随机向量X(X1,X2)密度函数。

证：只要验证满足密度函数两个条件即可

（1）显然，当x10,x20时有f(x1,x2)0

(x1x2)1 （2）edx1dx2ex2dx2ex2dx1dx2edx1dx2e000000定义2.4 设X(X1,X2,,Xp)是p维随机向量，称由它的q(p)个分量组成的子向量

(x1x2)(x1x2)X(i)(Xi1,Xi2,q,Xiq)的分布为X的边缘（或边际）分布，相对地把X的分布称为联合分布。通过变

换X中各分量的次序，总可假定X(1)正好是X的前q个分量，其余pq个分量为X(2)，则

x(1)X(1)X(2)，相应的取值也可分为两部分x(2)。

xXpq当X的分布函数是F(x1,x2,,xq)时，X(1)的分布函数即边缘分布函数为：

F(x1,x2,,xq)P(X1x1,,Xqxq)

,Xqxq,Xq1,,)

,Xp)

,xq,, P(X1x1, F(x1,x2,为：

(1)当X有分布密度f(x1,x2,,xp)时（亦称联合分布密度函数），则X也有分布密度，即边缘密度函数

f1(x1,x2,,xq)f(x1,,xp)dxq1,,dxp

【例2.2】对例2.1中的X(X1,X2)求边缘密度函数。

(x1x2)dx2ex1,x10e解： f(x1)f(x1,x2)dx2＝0

0,其它ex2,x20 同理f(x2)＝

0,其它定义2.5 若p个随机变量X1,X2,,Xp的联合分布等于各自的边缘分布的乘积，则称X1,X2,,Xp是

相互独立的。

【例2.3】 问例2.2中的X1与X2是否相互独立？

e(x1,x2),x10,x20 解：f(x1,x2)＝

其它0,ex1,x10ex2,x20 fx1(x1)＝ fx2(x2)＝

其它其它0,0,由于f(x1,x2)＝fx1(x1)fx2(x2)，故X1与X2相互独立。

这里我们应该注意，由X1,X2,定义

2.6

设

,Xp相互独立，可推知任何Xi与Xj(ij)独立，但反之不真。

X(X1,X2,,Xp)，若E(Xi)(i1,存p在且有限，则称

E(X)(E(X1),E(X2),为μ和i，即μ(1,2,,E(Xp))为X的均值（向量）或数学期望，有时也把E(X)和E(Xi)分别记,p)，容易推得均值（向量）具有以下性质：

（1）E(AX)AE(X) （2）E(AXB)AE(X)B

（3）E(AXBY)AE(X)BE(Y)

其中，X、Y为随机向量，A、B为大小适合运算的常数矩阵。

定义2.7 设X(X1,X2,,Xp)，Y(Y1,Y2,,Yp)，称D(X)E(XE(X))(XE(X))

Cov(X1,Xp)Cov(X1,X1)Cov(X1,X2)Cov(X,X)Cov(X,X)Cov(X,X)21222p（2.4） Cov(X,X)Cov(X,X)Cov(X,X)p1p2pp为X的方差或协差阵，有时把D(X)简记为Σ，Cov(Xi,Xj)简记为ij，从而有Σ(ij)pp；称随机向量X和Y的协差阵为

Cov(X1,Yp)Cov(X1,Y1)Cov(X1,Y2)Cov(X,Y)Cov(X,Y)Cov(X,Y)21222p（2.5） Cov(X,Y)E(XE(X))(YE(Y))Cov(X,Y)Cov(X,Y)Cov(X,Y)p1p2pp当X=Y时，即为D(X)。

若Cov(X,Y)0，则称X和Y不相关，由X和Y相互独立易推得Cov(X,Y)0，即X和Y不相关；但反过来，当X和Y不相关时，一般不能推知它们独立。

当A、B为常数矩阵时，由定义可以推出协方差阵有如下性质： （1）对于常数向量a，有D(Xa)D(X) （2）D(AX)AD(X)AAΣA （3）Cov(AX,BY)ACov(X,Y)B

（4）设X为n维随机向量，期望和协方差存在，记μE(X)，ΣD(X)，A为nn常数阵，则 E(XAX)tr(AΣ)μAμ

这里我们应该注意到，对于任何的随机向量X(X1,X2,,Xp)来说，其协差阵Σ都是对称阵，同时总

是非负定（半正定）的。大多数情况是正定的。

若X(X1,X2,,Xp)的协差阵存在，且每个分量的方差大于零，则称随机向量X的相关阵为

RCorr(X)(ij)pp，其中

ijCov(Xi,Xj)D(Xi)D(Xj)ijiijj i,j1,,p （2.6）

为Xi与Xj的相关系数。

在数据处理时，为了克服由于指标的量纲不同对统计分析结果带来的影响，往往在使用各种统计分析之前，常需要将每个指标“标准化”，即进行如下变换X那么由（2.7）构成的随机向量X(X1,X2,****jXjE(Xj)D(Xj)， j1,,p （2.7） ,pp)，有：

Cdiag(11,22,,X*p)。令，

X*C1(XE(X))

那么，标准化后的随机向量X均值和协差阵分别为

*E(X*)E[C1(XE(X))]C1E[(XE(X))]0

D(X*)DC[1X(EX(111))C]DX[(EX11C()1)]CD(X)CCΣCR即标准化数据的协差阵正好是原指标的相关阵。

第三节 多元正态分布

一、多元正态分布的定义

f(x)我们先来回顾一元正态分布的密度函数，即为上式可以改写为f(x)1e2(x)222,0

1121exp(x)()(x)（2.8） 1/221/2(2)()2由于（2.8）式中的x，均为一维的数字，可以用(x)代表(x)的转置。根据上面的表述形式，

我们可以将其推广，给出多元正态分布的定义。

定义2.8 若p维随机向量X(X1,X2,,Xp)的密度函数为：

f(x1,x2,其中x(x1,x2,,xp)1(2)p/21-1 （2.9） exp(x-μ)Σ(x-μ)122Σ,xp)，μ是p维随机向量，Σ是p阶正定阵，则称X服从p元正态分布，也称X为

p维正态随机向量，简记为X~NP(μ,Σ)，显然当p1时，即为一元正态分布密度函数。

可以证明μ为X的均值（向量），Σ为X的协差阵。

这里我们应该提及的是，当Σ0时，Σ1不存在，X也就不存在通常意义下的密度函数，然而可以形

式的给出一个表达式，是的有些问题可以利用这一形式对Σ0及Σ0的情况给出一个统一的处理。 当p2时，设X(X1,X2)服从二元正态分布，则

12111212，r1 Σ2221222122这里1，2分别是X1与X2的方差，是X1与X2的相关系数。即有

2|Σ|122(12)

Σ12212212(12)2112 12故X1与X2的密度函数为

(x11)211f(x1,x2)exp22212(12)1/22(1)12(x11)(x22)2(x22)2122对于0，那么X1与X2是相互独立的；若0，则X1与X2趋于正相关；若0，则X1与X2趋

于负相关。

定理2.1 设X~NP(μ,Σ)，则有E(X)μ，D(X)Σ。

关于这个定理的证明可以参考文献[1]，该定理将多元正态分布的参数μ和Σ赋予了明确的统计意义。 这里我们需要明确的是，多元正态分布的定义不止是一种，更广泛的可以采用特征函数来定义，也可以用一切线性组合均为正态的性质来定义。 二、多元正态分布的性质

在讨论多元统计分析的理论和方法时，经常用到多元正态变量的某些性质，利用这些性质可使得正态分布的处理变得容易一些。

1．若X(X1,X2,,Xp)~Np(μ,Σ)，Σ是对角阵，则X1,,Xp相互独立。 2．若X~Np(μ,Σ)，A为sp阶常数阵，d为s维常数向量，则 AXd~Ns(Aμd,AΣA) 即正态随机向量的线性函数还是正态的。

3．若X~Np(μ,Σ)，将X，μ，Σ作如下剖分

Xμ(1)Σ11 X(2) μ(2) ΣΣ21Xpqμpq则X(1)~Nq(μ(1),Σ11)，X(2)~Np-q(μ(2),Σ22)。

(1)qqqpqΣ12

Σ22pqq这里需要指出的是：第一，多元正态分布的任何边缘分布为正态分布，但反之不真。第二，由于

Σ12Cov(X(1),X(2))，故Σ120表示X(1)和X(2)不相关，因此可知，对于多元正态变量而言，X(1)和

X(2)的不相关与独立是等价的。

【例2.4】 若X(X1,X2,X3)~N3(μ,Σ)

其中，

11112 2212233132（1）

1310023 设 a(0,1,0) ，A，则 00133X1 aX(0,1,0)X2X2~N(aμ,aΣa) X3其中

111) aμ(0,1,0)22 aΣa(0,1,021331（2）

122232130231 03322X1100X1 AXX2X~N(Aμ,AΣA)

001X33其中

1111001100AΣA Aμ 2210013001331（3）

1222321310001113 2333013133X11(1)(1)X111213XμΣ11Σ1222 记X μ Σ21222 3Σ(2)(2)2122μX323331X33则 X其中 μ(1)(1)X1~N2(μX2(1),Σ11 )121 Σ1111 22122,Xp)服从p元正态分布，则它的每个分量必服从一元正态分

在此我们应该注意到，如果X(X1,X2,布，因此把某个分量的n个样品值作成直方图，如果断定不呈正态分布，则就可以断定随机向量

X(X1,X2,,Xp)也不可能服从p元正态分布。

第四节 多元正态分布的参数估计

一、多元样本的数字特征 设样本资料可用矩阵表示为

X11X21XXn1X12X22Xn2X1pX(1)XX2p(X,X，，X)(2)

12pXnpX(n) ,n。

在这里我们给出样本均值向量、样本离差阵、样本协差阵以及样本相关阵的定义。

定义2.9 设X(1),X(2),,X(n)为来自p元总体的样本，其中X(a)(Xa1,Xa2,,Xap)，a1,2,1nˆXX(a)(X1,X2,（1） 样本均值向量定义为μna1其中

,Xp)

1nX(a)na1X11X211X12X22nXX1p2pnXn1X11X21XXX1n22212nXnpX1pX2pXn1X1Xn2X2

XnpXp

（2）样本离差阵定义为

Spp这里，

(Xa1(a)X)(X(a)X)(sij)pp （2.11）

Xa1X1nnXX2(X(a)X)(X(a)X) a2(Xa1X1,Xa2X2,a1a1XXpap2(Xa1X1)(Xa1X1)(Xa2X2)n2(XX)(XX)(XX)a22a22a11a1(XX)(XX)(XX)(XX)pappa11a22aps11s12s1psss21222p(s) ijppsssp1p2pp11n（3）样本协差阵定义为VppS(X(a)X)(X(a)X)(vij)pp （2.12）

nna1,XapXp)

(Xa1X1)(XapXp)(Xa2X)(XapXp)22(XapXp)11n1n这里，S(X(a)X)(X(a)X)(XaiXi)(XajXj)vijnna1na1pppp

ˆr （2.13） （4）样本相关阵定义为Rppijpp 其中rijvijviivjjsijsiisjj

在此，我们应该提及的是，样本均值向量和离差阵也可用样本资料阵X直接表示如下： Xp11X1n 其中 1n(1,1,n,1)

X1111X12 由于 Xp1X1n nnX1p

那么，（2.11）式可以表示为：

nX21X22X2pXn11X11X21XX1Xn212212nXnp1X1pX2pXn1X1Xn2X2

XnpXpS(X(a)X)(X(a)X)XXnXXXXa111X1n1XX(I1n1nnn)X

nn（2.14）

1其中 In00 1二、均值向量与协差阵的最大似然估计

多元正态分布有两组参数，均值μ和协差阵Σ，在许多问题中它们是未知的，需要通过样本来估计。那么，通过样本来估计总体的参数叫做参数估计，参数估计的原则和方法是很多的，这里用最常见的且具有很多优良性质的最大似然法给出μ和Σ的估计量。 设X(1),X(2),,X(n)来自正态总体Np(μ,Σ)容量为n的样本，每个样品X(a)(Xa1,Xa2,,Xap)，

a1,2,,n，样本资料阵为（2.1）式表示，即

X1pX11X12XXX21222p XXXXn2npn11S （2.15） n实际上，最大似然法求估计量可以这样得到。针对X(1),X(2),,X(n)来自正态总体Np(μ,Σ)容量为n的ˆˆX ，Σ则可由最大似然法求出μ和Σ的估计量，即有μ样本，构造似然函数，即

L(μ,Σ)f(Xi,μ,Σ)i1n1(2)pn/21n1exp(Xμ)Σ(Xμ)i （2.16） in/22Σi1为了求出使（2.16）式取极值的μ和Σ的值，将（2.16）两边取对数，即

1n1nlnL(μ,Σ)pnln(2)lnΣ(Xiμ)Σ1(Xiμ) （2.17）

222i1因为对数函数是一个严格单调增函数，所以可以通过对lnL(μ,Σ)的极大值而得到μ和Σ的估计量。

(XAX)(XAX)2AX，XX，这里我们要注意到，根据矩阵代数理论，对于实对称矩阵A，有

XAlnAA1。 A那么，针对对数似然函数（2.17）分别对μ和Σ求偏导数，则有

lnL(μ,Σ)n1Σ(Xiμ)0μi1（2.18）由（2.18）式可以得到极大似nlnL(μ,Σ)nΣ11(Xμ)(Xμ)(Σ1)20iiΣ22i1然估计量分别为

1nˆXiXμni1 nΣˆ1(XX)(XX)1Siini1n由此可见，多元正态总体的均值向量μ的极大似然估计量就是样本均值向量，其协差阵Σ的极大似然

估计就是样本协差阵。

μ和Σ的估计量有如下基本性质： 1．E(X)μ，即X是μ的无偏估计； E(S)2．X，

1nn1111Σ，即S不是Σ的无偏估计，而E(S)Σ，即S是Σ的无偏估计； nnn1n11S分别是μ，Σ的有效估计； n1113．X，S（或。 S）分别是μ，Σ的一致估计（相合估计）

nn1样本均值向量和样本离差阵在多元统计推断中具有十分重要的作用，并有如下结论： 定理2.2 设X和S分别是正态总体Np(μ,Σ)的样本均值向量和离差阵，则 1．X~Np(μ,1Σ)； n2．离差阵S可以写为SZa1n1a 其中，Z1,Za,Zn1独立同分布于Np(0,Σ)；

3．X和S相互独立；

4．S为正定阵的充要条件是n三、Wishart分布

p。

1S来估计μ和Σ，前面已指出，均值向量X的分布仍为正态分布，n12而离差阵S的分布又是什么呢？为此给出维希特（Wishart）分布，并指出它是一元分布的推广，也是ˆ在实际应用中，常采用X和Σ构成其它重要分布的基础。

Wishart分布是Wishart在1928年推导出来的，而该分布的名称也即由此得来。

定义2.10 设X(a)(Xa1,Xa2,,Xap)~Np(μa,Σ)，a1,2,,n且相互独立，则由X(a)组成的随机矩阵：Wpp其中Z(a1,Xa1n(a)X(a) （2.19）的分布称为非中心Wishart分布，记为Wp(n,Σ,Z)。

，μa称为非中心参数；,an)μaμa当μa0时称为中心Wishart分布，

a1n,an)(a1,记为Wp(n,Σ)，当np，Σ0，Wp(n,Σ)有密度存在，其表达式为：

1(np1)112wexptrΣw2,当w为正定阵pni1f(w)np2p(p1)4n2 （2.20） 2Σ2i1其它0,222显然，当p1，Σ时，f(w)就是(n)的分布密度，此时（2.19）式为

nn1n22WX(a)X(a)X(a)，有2X(a)~2(n)。因此，Wishart分布是2分布在p维正态情况下

a1a1a1的推广。

下面给出Wishart分布的基本性质：

１．若X(a)~Np(μ,Σ)，a1,2,,n且相互独立，则样本离差阵

1nS(X(a)X)(X(a)X)~Wp(n1,Σ)，其中XX(a)。

na1a12．若Si~Wp(ni,Σ)，i1,,k，且相互独立，则

nSi1ki~Wp(ni,Σ)。

i1k3．若Xpp~Wp(n,Σ)，Cpp为非奇异阵，则CXC~Wp(n,CΣC)。

这里我们有必要说明一下什么是随机矩阵的分布。随机矩阵的分布有不同的定义，此处是利用已知向量分

布的定义给出矩阵分布的定义。

这里我们有必要说明一下什么是随机矩阵的分布。随机矩阵的分布有不同的定义，此处是利用已知向量分布的定义给出矩阵分布的定义。 设随机矩阵

X11X21 XXn1X12X22Xn2X1pX2p

Xnp将该矩阵的列向量（或行向量）一个接一个地连接起来，组成一个长的向量，即拉直向量：

(X11,X21,,Xn1,X12,X22,,Xn2,,X1p,X2p,,Xnp)的分布定义为该阵的分布。若X为对称阵时，由于XijX，

npn，故只取其下三角部分组成的拉直向量，即

1(X11,X2,1X,X,2n2,X,。n2, pX,)

第三章 多元正态分布均值向量和协差阵的检验

第一节 引言

在单一变量的统计分析中，已经给出了正态总体N（ ， 2） 的均值和方差2的各种检验。对于多变量的正态总体Np（ ， ∑ ） ，各种实际问题同样要求对和∑进行统计推断。

例如，我们要考察全国各省、自治区和直辖市的社会经济发展状况，与全国平均水平相比较有无显著性差异等，就涉及到多元正态总体均值向量的检验问题等。

本章类似单一变量统计分析中的各种均值和方差的检验，相应地给出多元统计分析中的各种均值向量和协差阵的检验。

其基本思想和步骤均可归纳为：

第一，提出待检验的假设H0和H1；

第二，给出检验的统计量及其服从的分布； 第三，给定检验水平，查统计量的分布表，确定相应的临界值，从而得到否定域； 第四，根据样本观测值计算出统计量的值，看是否落入否定域中，以便对待判假设做出决策（拒绝或接受）。

在检验的过程中，关键在于对不同的检验给出不同的统计量，而有关统计量的给出大多用似然比方法得到。由于多变量问题的复杂性，本章只侧重于解释选取统计量的合理性，而不给出推导过程，最后给出几个实例。

为了更好的说明检验过程中统计量的分布，本章还要介绍HotellingT2分布和Wilks分布的定义。

第二节 均值向量的检验

2

一、单一变量检验的回顾及HotellingT分布

为了对多元正态总体均值向量作检验，首先需要给出HotellingT2分布的定义。 在单一变量的检验问题中，设X1,X2,,Xn来自总体N(,2)的样本，我们要检验假设

H0:0;H1:0

当已知时，用统计量z2

(X0)n （3.1）

1n其中，XXi为样本均值。当假设成立时，统计量z服从正态分布z~N(0,1)，从而否定域为

ni1|z|z/2，z/2为N(0,1)的上/2分位点。

1n(XiX)2 （3.2） 当未知时，用 Sn1i1(X0)2

n （3.3） 作为的估计量，用统计量：tS来做检验。当假设成立时，统计量t服从自由度为n1的t分布，从而否定域为|t|t/2(n1)，t/2(n1)为自由度为n1的t分布上的/2分位点。

2

2这里我们应该注意到，（3.3）式可以表示为

n(X)2tn(X)(S2)1(X) （3.4） 2S2

对于多元变量而言，可以将t分布推广为下面将要介绍的HotellingT分布。

22-1定义3.1 设X~Np(μ,Σ)，S~Wp(n,Σ)且X与S相互独立，np，则称统计量TnXSX的

分布为非中心HotellingT2分布，记为T~T(p,n,μ)。当μ0时，称T2服从（中心）HotellingT2分布。记为T(p,n)。

由于这一统计量的分布首先由Harold Hotelling 提出来的，故称为HotellingT2分布，值得指出的是，我国著名统计学家许宝禄先生在1938年用不同方法也导出T2分布的密度函数，因表达式很复杂，故略去。 在单一变量统计分析中，若统计量t~t(n1)分布，则t~F(1,n1)分布，即把t分布的统计量转化为F统计量来处理，在多元统计分析中T2统计量也具有类似的性质。

21定理3.1 若X~Np(0,Σ)，S~Wp(n,Σ)且X与S相互独立，令TnXSX，则

2222

np12T~F(p,np1) （3.5） np,X(n)是来自p维正态总体Np(μΣ,在我们后面所介绍的检验问题中，经常会用到这一性质。 设X(1),X(2),n1n)的样本，且XX()，

n1S(X(a)X)(X(a)X)。

a1（一） 协差阵Σ已知时均值向量的检验 H0：μμ0（μ0为已知向量）H1：μμ0 假设H0成立，检验统计量为

212 T0n(Xμ0)Σ(Xμ0)~(p) （3.6）

2给定检验水平，查分布表使PT02，可确定出临界值，再用样本值计算出T0，若2T02，则否定H0，否则接受H0。

222这里要对统计量的选取做一些解释，为什么该统计量服从(p)分布。根据二次型分布定理知道，若

2X~Np(0,Σ)，则XΣ1X~2(p)。显然，

T02n(Xμ0)Σ1(Xμ0)n(Xμ0)Σ1n(Xμ0)YΣ1Y

n(Xμ0)~Np(0,)，因此，T02n(X0)Σ1(Xμ0)~2(p) （二）协差阵Σ未知时均值向量的检验

H0：μμ0（μ0为已知向量）H1：μμ0 假设H0成立，检验统计量为

(n1)p12T~F(p,np) （3.7）

(n1)p其中，Y其中，T2(n1)[n(Xμ0)S1n(Xμ0)] 给定检验水平，查F分布表，使P若

np22

可确定出临界值F，再用样本值计算出T，TF，

(n1)pnp2TF，则否定H0，否则接受H0。

(n1)p这里需要解释的是，当Σ未知时，自然想到要用样本协差阵估计量，而样本离差阵

S

n1S取代替Σ，因(n1)S1是Σ1的无偏n1(Xa1(a)X)X((a)X)W~pn(Σ1 ,)n(Xμ0)~Np(0,Σ)

由定义3.1知

T2(n1)[n(Xμ0)S1n(Xμ0)]~T2(p,np) 再根据Hotelling T2分布的性质，所以

(n1)p12T~F(p,np)

(n1)p在处理实际问题时，单一变量的检验和多变量检验可以联合使用，多元的检验具有概括和全面考察的特点，而一元的检验容易发现各变量之间的关系和差异，能给人们提供更多的统计分析信息。 三、两个正态总体均值向量的检验

（一）当协差阵相等时，两个正态总体均值向量的检验

设X(a)(Xa1,Xa2,,Xap)，a1,2,,n，为来自p维正态总体Np(μ1,Σ)的容量为n的样本；

Y(a)(Ya1,Ya2,，a1,2,,m，为来自p维正态总体Np(μ2,Σ)的容量为m的样本。两组样本,Yap)1n1m相互独立，np,mp，且XX(i)，YY(i)。

ni1mi11．针对有共同已知协差阵的情形

对假设 H0：μ1μ2 H1：μ1μ2 进行检验。 对此问题，假设H0成立时，所构造的检验统计量为

nm(XY)Σ1(XY)~2(p) （3.8） nm22给出检验水平，查(p)分布表使PT02，可确定出临界值2，再用样本值计算出T02，若

T022T02，则否定H0，否则接受H0。

这里，我们应该注意到，在单一变量统计中进行均值相等检验所给出的统计量为

zXY2nm(XY)2nmnm22(XY)显然z2(XY)(2)1(XY)~2(1) 22(nm)nmnm此式恰为上边统计量当p1时的情况，不难看出这里给出的检验统计量是单一变量检验情况的推广。

2．针对有共同的未知协差阵的情形

对假设H0：μ1μ2 H1：μ1μ2进行检验。 对此问题，假设H0成立时，所构造的检验统计量为

2~N(0,1)

F(nm2)p12T~F(p,nmp1) （3.9）

(nm2)p其中，

nm1nmT(nm2)(XY)S(XY)

nmnm2)SSxSy Sx(X(a)X)(X(a)X， X(X1,X2,a1n,Xp)

Sy(Y(a)Y)(Y(a)Y)， Y(Y1,Y2,a1n,Yp)

给定检验水平，查F分布表，使pFF，可确定出临界值F，再用样本值计算出F，若

FF，则否定H0，否则接受H0。

这里我们需要解释的是，当两个总体的协差阵未知时，自然想到用每个总体的样本协差阵

1Sx和n11Sy去代替，而 m1 Sx Sy(Xa1mn(a)X)(X(a)X)~Wp(n1,Σ) Y)(Y(a)Y)~Wp(m1,Σ)

(Ya1(a)从而SSxSy~Wp(nm2,Σ)。又由于

nm(XY)~Np(0,Σ) nm(nm2)p12T~F(p,nmp1) 所以

(nm2)p

下述假设检验统计量的选取和前边统计量的选取思路是一样的，以下只提出待检验的假设，然后给出统计量及其分布，为节省篇幅，不做重复解释。

（二）协差阵不等时，两个正态总体均值向量的检验

设从两个总体Np(μ1,Σ1)和Np(μ2,Σ2)中，分别抽取两个样本，即X(a)(Xa1,Xa2,,Xap)，

a1,2,,n；Y(a)(Ya1,Ya2,,Yap)，a1,2,,m，其容量分别为n和m，且两组样本相互独立，

np,mp，Σ10，Σ20。

对假设H0：μ1μ2 H1：μ1μ2进行检验。 1．针对nm的情形

令 Z(i)X(i)Y(i) i1,2,,n

1n ZZ(i)XY

ni1 S(Zi1n(i)Z)(Z(i)Z)(X(i)Y(i)XY)(X(i)Y(i)XY)

i1n假设H0成立时，构造检验统计量为

(np)nZS-1Z~F(p,np) （3.10） p2．针对nm的情形 F

在此，我们不妨假设nm，令

nn11nZ(i)X(i)Y(i)Y(i)Y(i) i1,2,,n mmj1nmj11nZZ(i)XY

ni1S(Z(i)Z)(Z(i)Z)i1nn1nn1n(X(i)X)(Y(i)Y(j))(X(i)X)(Y(i)Y(j))

mnmnj1i1j1假设H0成立时，构造检验统计量为

(np)nZS-1Z~F(p,np) Fpn四、多个正态总体均值向量的检验

解决多个正态总体均值向量的检验问题，实际上应用到多元方差分析的知识。多元方差分析是单因素方差分析直接的推广。为了容易理解多元方差分析方法，我们有必要先回顾单因素方差分析方法。 （一）单因素方差分析的基本思想及Wilks分布

22设k个正态总体分别为N(1,)，,N(k,)，从k个总体取ni个独立样本如下：

(1)(1) X1(1),X2,,Xn1 

(k)(k) X1(k),X2 ,,Xnk H0：12k H1：至少存在ij使ij 假设H0成立时，构造检验统计量为

FSSA(k1)~F(k1,nk) （3.11）

SSE(nk)这里SSAknin(Xii1kiX)称为组间平方和；SSE(X(ji)Xi)2 称为组内平方和；

2knii1j1SST(X(ji)X)2称为总平方和。其中

i1j11 XiniXj1n(i)j

1kni(i) XXj nn1nk

ni1j1给定检验水平，查F分布表，使pFF，可确定出临界值F，再用样本值计算出F值，若

FF，则否定H0，否则接受H0。

定义3.2 若X~Np(0,Σ)，则称协差阵的行列式Σ为X的广义方差。称

1S为样本广义方差。其中nS(X(a)X)(X(a)X)。

a1n定义3.3 若A1~Wp(n1,Σ),n1p，A2~Wp(n2,Σ),Σ0，且A1和A2相互独立，则称

A1 

A1A2为Wilks统计量，的分布称为Wilks分布，简记为~(p,n1,n2)，其中n1，n2为自由度。 这里我们需要说明的是，在实际应用中经常把统计量化为T统计量进而化为F统计量，利用F统计量来解决多元统计分析中有关检验问题。表3.1列举常见的一些情形。

表3.1 与F统计量的关系 p n1 n2 F统计量及分别 任意 任意 1 2

n1p11(p,n1,1)~F(p,n1p1) p(p,n1,1)任意 任意 2 n1p1(p,n1,2)~F(2p,2(n1p)) p(p,n1,2)n11(1,n1,n2)~F(n2,n1) n2(1,n1,n2)1 任意 任意 2 任意 任意 n111(2,n1,n2)~F(2n2,2(n11)) n2(2,n1,n2)2以上几个关系式说明对一些特殊的统计量可以化为F统计量，而当n22，p2时，可用统计量或F统计量来近似表示，后面给出。

（二）多元方差分析法

设有k个p维正态总体Np(μ1,Σ)，,Np(μk,Σ)，从每个总体抽取独立样本个数分别为n1,n2,,nk，

n1nkn，每个样品观测p个指标得观测数据如下：

(1),Xip)，i1,2,,n1 (2),Xip)，i1,2,,n2 (k),Xip)，i1,2,,nk

(1) 第一个总体： Xi(1)(Xi(1),X1i2,(2)第二个总体： Xi(2)(Xi(2)1,Xi2,…… …… ……

k) 第k个总体： Xi(k)(Xi(1k),Xi(2,全部样品的总均值向量： X1p1knr(r)Xi(X1,X2,nr1i11nr(r)(r)Xi(X1,X(2r),ni1(r)j,Xp)

各总体样品的均值向量： X(r)1p,X(pr))，r1,2,,k

此处 X1nrXi1nr(r)ij

类似一元方差分析办法，将诸平方和变成了离差阵即： A E Tn(Xrr1knrr1i1k(r)X)(X(r)X) X(r))(Xi(r)X(r))

(X(r)i(Xr1i1knr(r)iX)(Xi(r)X)

这里，我们称A为组间离差阵；E为组内离差阵；T为总离差阵。很显然有 TAE。 我们的问题是检验假设

μk H1：至少存在ij使μiμj H0：μ1μ2用似然比原则构成的检验统计量为

EE~(p,nk,k1) （3.13） TAE给定检验水平，查Wilks分布表，确定临界值，然后作出统计判断。在这里我们特别要注意，Wilks分

2布表可用分布或F分布来近似。

2巴特莱特（Bartlett）提出了用分布来近似。设~(p,n,m)，令

1/t V(nm(pm1)2)lnln （3.14）

2则V近似服从(pm)分布。其中，tnm(pm1)2。 Rao后来又研究用F分布来近似。设~(p,n,m)，令

11/LtL2 R （3.15） 1/Lpm则R近似服从F(pm,tL2)，这里tL2不一定为整数，可用与它最近的整数来作为F的自由度，且min(p,m)2。其中，

p2m24pm2 tnm(pm1)2 L2 2pm54第三节 协差阵的检验

一、一个正态总体协差阵的检验 设X(a)(Xa1,Xa2,1/2,Xap)(a1,2,,n)来自p维正态总体Np(μ,Σ)的样本，Σ未知，且Σ0。

首先，我们考虑检验假设

H0：ΣIp H1：ΣIp

1n/2e所构造的检验统计量为exptrSS2n其中 Snp/2 （3.16）

(Xa1n(a) X)X(a()X)然后，我们考虑检验假设 H0：ΣΣ0Ip H1：ΣΣ0Ip 因为Σ00，所以存在D(D0)，使得DΣ0DIp。 令 Y(a)DX(a) a1,2,,n 则 Y(a)~Np(Dμ,DΣD)Np(μ,Σ) 因此，检验ΣΣ0等价于检验Σ*Ip

np/2**1n/2e此时构造检验统计量为exptrS*S*2n其中 S* （3.17）

(Ya1n(a)Y)(Y(a)Y)

222ln分布。因此当，由样本值计算出值，若npp(p1)/2给定检验水平，因为直接由分布计算临界值0很困难，所以通常采用的近似分布。 在H0成立时，2ln极限分布是即e2/2，则拒绝H0，否则接受H0。

k设有k个正态总体分别为Np(μ1,Σ1)，,Np(μk,Σk)，Σi0且未知，i1,,k。从k个总体分别取

ni个样本X(i)(a)(X,X,(i)a1(i)a2,X) i1,,k；a1,,ni这里nin为总样本容量。

(i)api1我们考虑检验假设 H0：Σ1Σ2Σk H1：Σi不全相等

np/2构造检验统计量为kn

kSi1(i)(a)kni/2iSn/2ni1kipni/2 （3.18）

1nii()(X)其中SSi Si(XX)X XX(a)

ni11i1，称为修正的统巴特莱特（Bartlett）建议 ，将ni改为ni1，从而n变为nk，变换以后的k记为k1近似分布2fp(p1)(k1) 计量，则2lnk。 其中(1D)f22p23p1k11 ninj，6(p1)(k1)i1ni1nk D2(2p3p1)(k1)n1n2nk6(p1)(nk),(i)i()(a)i()ni(i)第四章 判别分析

第一节 引言

在我们的日常生活和工作实践中，常常会遇到判别分析问题，即根据历史上划分类别的有关资料和某种最优准则，确定一种判别方法，判定一个新的样本归属哪一类。例如，某医院有部分患有肺炎、肝炎、冠心病、糖尿病等病人的资料，记录了每个患者若干项症状指标数据。现在想利用现有的这些资料找出一种方法，使得对于一个新的病人，当测得这些症状指标数据时，能够判定其患有哪种病。又如，在天气预报中，我们有一段较长时间关于某地区每天气象的记录资料（晴阴雨、气温、气压、湿度等），现在想建立一种用连续五天的气象资料来预报第六天是什么天气的方法。这些问题都可以应用判别分析方法予以解决。

把这类问题用数学语言来表达，可以叙述如下：设有n个样本，对每个样本测得p项指标（变量）的数据，已知每个样本属于k个类别（或总体）G1，G2， …，Gk中的某一类，且它们的分布函数分别为F1(x)，F2(x)， …，Fk(x)。我们希望利用这些数据，找出一种判别函数，使得这一函数具有某种最优性质，能把属于不同类别的样本点尽可能地区别开来，并对测得同样p项指标（变量）数据的一个新样本，能判定这个样本归属于哪一类。

判别分析内容很丰富，方法很多。判断分析按判别的总体数来区分，有两个总体判别分析和多总体判别分析；按区分不同总体所用的数学模型来分，有线性判别和非线性判别；按判别时所处理的变量方法不同，有逐步判别和序贯判别等。判别分析可以从不同角度提出问题，因此有不同的判别准则，如马氏距离最小准则、Fisher准则、平均损失最小准则、最小平方准则、最大似然准则、最大概率准则等等，按判别准则的不同又提出多种判别方法。本章仅介绍常用的几种判别分析方法：距离判别法、Fisher判别法、Bayes判别法和逐步判别法。

第二节 距离判别法

一、马氏距离的概念

设p维欧氏空间Rp中的两点X(X1,X2,2,Xp)和Y(Y1,Y2,2,Yp)，通常我们所说的两点之间的距

离，是指欧氏距离，即 d(X,Y)(X1Y1)2(XpYp)2 （4.1）

2在解决实际问题时，特别是针对多元数据的分析问题，欧氏距离就显示出了它的薄弱环节。

第一、设有两个正态总体，X~N(1,)和Y~N(2,4)，现有一个样品位于如图4.1所示的A点，距总体X的中心2远，距总体Y的中心3远，那么，A点处的样品到底离哪一个总体近呢？若按欧氏

距离来量度，A点离总体X要比离总体Y “近一些”。但是，从概率的角度看，A点位于1右侧的2x处，而位于2左侧1.5y处，应该认为A点离总体Y“近一些”。显然，后一种量度更合理些。

第二、设有量度重量和长度的两个变量X与Y，以单位分别为kg和cm得到样本A(0,5)，B(10,0)，C(1,0)，

D(0,10)。今按照欧氏距离计算，有

AB10252125； CD12102101如果我们将长度单位变为mm，那么，有

AB1025022600； CD1210021000 1量纲的变化，将影响欧氏距离计算的结果。

为此，我们引入一种由印度著名统计学家马哈拉诺比斯（Mahalanobis, 1936）提出的“马氏距离”的概念。 设X和Y是来自均值向量为μ，协方差为Σ(0)的总体G中的p维样本，则总体G内两点X与Y之间的马氏距离定义为

D(X,Y)(XY)Σ(XY) (4.2) 定义点X到总体G的马氏距离为

D(X,G)(Xμ)Σ(Xμ) (4.3)

这里应该注意到，当ΣI（单位矩阵）时，即为欧氏距离的情形。 二、距离判别的思想及方法

21211、两个总体的距离判别问题

问题：设有协方差矩阵∑相等的两个总体G1和G2，其均值分别是1和 2，对于一个新的样品X，要判断它来自哪个总体。

一般的想法是计算新样品X到两个总体的马氏距离D2（X，G1）和D2（X，G2），并按照如下的判别规则进行判断

XG1,XG2,如果如果D2(X,G1)D2(X,G2) （4.4） 22D(X,G1)D(X,G2)这个判别规则的等价描述为：求新样品X到G1的距离与到G2的距离之差，如果其值为正，X属于G2；否则X属于G1。 我们考虑

D2(X,G1)D2(X,G2)

(Xμ1)Σ1(Xμ1)(Xμ2)Σ1(Xμ2)1Σ1μ1(XΣ1X2XΣ1μ2μXΣ1X2XΣ1μ1μ12Σμ2)1Σ1μ1μ2XΣ1(μ2μ1)μ12Σμ22XΣ1(μ2μ1)(μ1μ2)Σ1(μ1μ2)μ1μ21 2XΣ(μ1μ2)22(Xμ)α2α(Xμ)11其中μ(μ1μ2)是两个总体均值的平均值，αΣ(μ1μ2)，记

2W(X)α(Xμ) （4.5）

则判别规则（4.4）式可表示为

XG1,如果W(X)0 （4.6） XG,如果W(X)02这里称W(X)为两总体距离判别的判别函数，由于它是X的线性函数，故又称为线性判别函数，α称为判

别系数。

在实际应用中，总体的均值和协方差矩阵一般是未知的，可由样本均值和样本协方差矩阵分别进行估计。设X1,(1)(2),X(1)n1来自总体G1的样本，X1,(1),X(2)n2是来自总体G2的样本，μ1和μ2的一个无偏估计分别为

1n1(1)1n2(2)(2)XXi 和 XXi

n1i1n2i11ˆ(S1S2) Σ的一个联合无偏估计为 Σn1n22这里 S(Xi1n()iX())(Xi()X()),1,2

ˆ(X)αˆ(XX) 此时，两总体距离判别的判别函数为 W其中X1(1)ˆ1(X(1)X(2))。这样，判别规则为 ˆΣ(XX(2))，α2XG1,XG2,如果如果ˆ(X)0W （4.7）

ˆW(X)022 这里我们应该注意到：

2（1） 当p1，G1和G2的分布分别为N(1,)和N(2,)时，1,2,均为已知，且12，

120，判别函数为W(x)(x) 2xxG1,如果判别规则为

xxG2,如果则判别系数为（2） 当μ1μ2，Σ1Σ2时，我们采用（4.4）式作为判别规则的形式。选择判别函数为

11W*(X)D2(X,G1)D2(X,G2)(Xμ1)Σ1(Xμ1)(Xμ2)Σ2(Xμ2)

XG1,它是X的二次函数，相应的判别规则为XG2,2、多个总体的距离判别问题

如果如果W*(X)0

W*(X)0问题：设有k个总体G1,G2,,Gk，其均值和协方差矩阵分别是μ1,μ2,,μk和Σ1,Σ2,,Σk，而且

Σ1Σ2ΣkΣ。对于一个新的样品X，要判断它来自哪个总体。

该问题与两个总体的距离判别问题的解决思想一样。计算新样品X到每一个总体的距离，即 D2(X,G)(Xμ)Σ1(Xμ)

Σ1XμΣ1μXΣ1X2μXC)XΣ1X2(I111这里IΣμ，CμΣμ，1,2,,k。

2由（4.8）式，可以取线性判别函数为

4.8）

XC， 1,2,,k W(X)I相应的判别规则为

XC) （4.9） XGi 如果 Wi(X)max(I1k针对实际问题，当μ1,μ2,,μk和Σ均未知时，可以通过相应的样本值来替代。设X1,体G中的样本（1,2,,k），则μ（1,2,,k）和Σ可估计为

())是来自总,X(n1kS ， 其中 nn1n2nk XXnk1i1同样，我们注意到，如果总体G1,G2,,Gk的协方差矩阵分别是Σ1,Σ2,,Σk，而且它们不全相等，则计算X到各总体的马氏距离，即

D2(X,G)(Xμ)Σ1(Xμ) 1,2,,k

()1nn()iˆ， 1,2,,k和 Σ则判别规则为

22 XGi 如果 D(X,Gi)minD(X,G) （4.10）

1k当μ1,μ2,,μk和Σ1,Σ2,,Σk均未知时，μ（1,2,,k）的估计同前，Σ（1,2,,k）的

ˆ估计为Σ1S， 1,2,,k n1三、判别分析的实质

我们知道，判别分析就是希望利用已经测得的变量数据，找出一种判别函数，使得这一函数具有某种最优性质，能把属于不同类别的样本点尽可能地区别开来。为了更清楚的认识判别分析的实质，以便能灵活的应用判别分析方法解决实际问题，我们有必要了解“划分”这样概念。

设R1，R2，…，Rk是p维空间R p的k个子集，如果它们互不相交，且它们的和集为R p，则称R1，R2， …，Rk为R p的一个划分。

在两个总体的距离判别问题中，利用W(X)α(Xμ)可以得到空间Rp的一个划分

R1{X:W(X)0} （4.11） R2{X:W(X)0}新的样品X落入R1推断XG1，落入R2推断XG2

这样我们将会发现，判别分析问题实质上就是在某种意义上，以最优的性质对p维空间R p构造一个“划分”，这个“划分”就构成了一个判别规则。这一思想将在后面的各节中体现的更加清楚。

第三节 贝叶斯（Bayes）判别法

从上节看距离判别法虽然简单，便于使用。但是该方法也有它明显的不足之处。 第一，判别方法与总体各自出现的概率的大小无关； 第二，判别方法与错判之后所造成的损失无关。Bayes判别法就是为了解决这些问题而提出的一种判别方法。

一、Bayes判别的基本思想

问题：设有k个总体G1,G2,,Gk，其各自的分布密度函数f1(x),f2(x),,fk(x)互不相同的，假设k个总体各自出现的概率分别为q1,q2,,qk（先验概率），qi0，

kqi1i1。假设已知若将本来属于Gi总体的样

品错判到总体Gj时造成的损失为C(j|i)，i,j1,2,,k。在这样的情形下，对于新的样品X判断其来自哪个总体。

下面我们对这一问题进行分析。首先应该清楚C(i|i)0、C(j|i)0，对于任意的i,j1,2,,k成立。设k个总体G1,G2,,Gk相应的p维样本空间为R1,R2,,Rk ，即为一个划分，故我们可以简记一个判别规则为R(R1,R2,,Rk)。从描述平均损失的角度出发，如果原来属于总体Gi且分布密度为fi(x)的样品，正好取值落入了Rj，我们就将会错判为属于Gj。 故在规则R下，将属于Gi的样品错判为Gj的概率为

P(j|i,R)fi(x)dx i,j1,2,,kRjij

如果实属Gi的样品，错判到其它总体G1,,Gi1,Gi1,Gk所造成的损失为

C(1|i),,C(i1|i),C(i1|i),C(k|i)，则这种判别规则R对总体Gi而言，样品错判后所造成的平

均损失为

r(i|R)[C(j|i)P(j|i,R)] i1,2,,k

j1k其中C(i|i)0

由于k个总体G1,G2,,Gk出现的先验概率分别为q1,q2,,qk，则用规则R来进行判别所造成的总平均损失为

g(R)qir(i,R)qiC(j|i)P(j|i,R) （4.12）

i1i1j1kkk所谓Bayes判别法则，就是要选择R1,R2,,Rk，使得（4.12）式表示的总平均损失g(R)达到极小。 二、Bayes判别的基本方法

设每一个总体Gi的分布密度为fi(x)，i1,2,,k，来自总体Gi的样品X被错判为来自总体Gj（i,j1,2,,k）时所造成的损失记为C(j|i)，并且C(i|i)0。那么，对于判别规则

R(R1,R2,,Rk)产生的误判概率记为P(j|i,R)，有P(j|i,R)fi(x)dx

Rj如果已知样品X来自总体Gi 的先验概率为qi，则在规则R下，由（4.12）式知，误判的总平均损失为

g(R)qiC(j|i)P(j|i,R)qiC(j|i)fi(x)dx(qiC(j|i)fi(x))dx （4.13）

i1j1i1j1Rjj1Rji1kkkkkk令

（4.13）式为 g(R)qC(j|i)f(x)h(x)，那么，

iiji1kkj1Rjhj(x)dx

*如果空间R有另一种划分R(R,R,,R)，则它的总平均损失为g(R)那么，在两种划分下的总平均损失之差为g(R)g(R)*p**1*2*kj1kR*jhj(x)dx

i1j1kkRiR*j[hi(x)hj(x)]dx （4.14）

由Ri的定义，在Ri上hi(x)hj(x)对一切j成立，故（4.14）式小于或等于零，这说明R1,R2,,Rk确能使总平均损失达到极小，它是Bayes判别的解。这样，我们以Bayes判别的思想得到的划分

R(R1,R2,,Rk)为Ri{x|hi(x)minhj(x)} i1,2,,k （4.15）

1jk具体说来，当抽取了一个未知总体的样本值X，要判断它属于哪个总体，只要前计算出k个按先验分布加权的误判平均损失

hj(x)qiC(j|i)fi(x)) j1,2,,k （4.16）

i1k然后比较这k个误判平均损失h1(x),h2(x),,hk(x)的大小，选取其中最小的，则判定样品X来自该总体。 这里我们看一个特殊情形，当k2时，由（4.16）式得h1(x)q2C(1|2)f2(x))

h2(x)q1C(2|1f)x( ))1从而R1{x|q2C(1|2)f2(x)q1C(2|1)f1(x)} R2{x|q2C(1|2)f2(x)q1C(2|1)f1(x)} 若令 V(x)f1(x)qC(1|2)， d2 f2(x)q1C(2|1)则判别规则可表示为

xG1,xG2,V(x)f1(x) f2(x)当V(x)d当V(x)d （4.17）

如果在此，f1(x)与f2(x)分别为N(μ1,Σ)和N(μ2,Σ)，那么

11exp(xμ1)Σ1(xμ1)(xμ2)Σ1(xμ2)22

exp[x(μ1μ2)/2]Σ1(μ1μ2)expW(x)其中W(x)由（4.5）所定义。于是，判定样品X来自该总体时，判别规则（4.17）成

XG1,XG2,如果如果W(X)lnd （4.18）

W(X)lnd对比判别规则（4.6），唯一的差别仅在于阈值点，（4.6）用0作为阈值点，而这里用lnd。当q1q2，

C(1|2)C(2|1)时，d1，lnd0，则（4.6）与（4.18）完全一致。

第四节 费歇（Fisher）判别法

Fisher判别法是1936年提出来的，该方法的主要思想是通过将多维数据投影到某个方向上，投影的原则是将总体与总体之间尽可能的放开，然后再选择合适的判别规则，将新的样品进行分类判别。

一、Fisher判别的基本思想

从k个总体中抽取具有p个指标的样品观测数据，借助方差分析的思想构造一个线性判别函数

U(X)u1X1u2X2upXpuX （4.19）

其中系数u(u1,u2,,up)确定的原则是使得总体之间区别最大，而使每个总体内部的离差最小。有了线性判别函数后，对于一个新的样品，将它的p个指标值代入线性判别函数（4.19）式中求出U(X)值，然后根据判别一定的规则，就可以判别新的样品属于哪个总体。 二、Fisher判别函数的构造 1、针对两个总体的情形 假设有两个总体G1,G2，其均值分别为μ1和μ2，协方差矩阵为Σ1和Σ2。当XGi时，我们可以求出uX的均值和方差，即

E(uX)E(uX|Gi)uE(X|Gi)uμii， i1,2

D(uX)D(uX|Gi)uD(X|Gi)uuΣiui2， i1,2

在求线性判别函数时，尽量使得总体之间差异大，也就是要求uμ1uμ2尽可能的大，即12变大；同时要求每一个总体内的离差平方和最小，即12，则我们可以建立一个目标函数

22(12) （4.20） 2212这样，将问题转化为，寻找u使得目标函数(u)达到最大。从而可以构造出所要求的线性判别函数。

(u)

2、针对多个总体的情形

假设有k个总体G1,G2,,Gk，其均值和协方差矩阵分别为μi

和Σi(0)（i1,2,,k）。同样，我们考虑线性判别函数uX，在XGi的条件下，有

E(uX)E(uX|Gi)uE(X|Gi)uμi i1,2,,k D(uX)D(uX|Gi)uD(X|Gi)uuΣiu i1,2,,k

令

b(uμiuμ) euΣiuu(Σi)uuEu

2i1i1i1k1k其中μμi，EΣi。这里b相当于一元方差分析中的组间差e相当于组内差，应用方差分析

ki1i1b的思想，选择u使得目标函数(u)（4.21）达到极大。

e这里我们应该说明的是，如果我们得到线性判别函数uX，对于一个新的样品X可以这样构造一个判别规则，如果uXuμjminuXuμi（4.22）则判定X来自总体Gj。

kkk1ik三、线性判别函数的求法

针对多个总体的情形，我们讨论使目标函数（4.21）式达到极大的求法。设X为p维空间的样品，那么

1k1μμiM1其中

ki1k1112 M1k21221kpkkμp1μ11p1μ2 11 1注意到 MMμ1μ2μ1μ2kμkμiμi

i1μk从而 b(uμi1ki1kiuμ)2

ku(μiμ)(μiμ)uu[μiμikμμ]ui111u(MMM11M)uuM(IJ)MuuBukk111。 这里，BM(IJ)M，Ipp 为pp的单位阵，Jk11uBu即有(u) （4.23）求使得（4.23）式达到极大的u。

uEu为了确保解的唯一性，不妨设uEu1，这样问题转化为，在uEu1的条件下，求u使得uBu式达到

极大。

考虑目标函数 (u)uBu(uEu1) （4.24） 对（4.24）式求导，有

2(BE)u0(4.25)u uEu10(4.26)对（4.25）式两边同乘u，有 uBuuEu

11从而，uBu的极大值为。再用E左乘（4.25）式，有(EBI)u0 （4.27）

11由（4.27）式说明为EB特征值，u为EB的特征向量。在此最大特征值所对应的特征向量u(u1,u2,,up)为我们所求结果。

这里值得注意的是，本书有几处利用极值原理求极值时，只给出了不要条件的数学推导，而有关充分条件的论证省略了，因为在实际问题中，往往根据问题本身的性质就能肯定有最大值（或最小值），如果所求的驻点只有一个，这时就不需要根据极值存在的充分条件判定它是极大还是极小而就能肯定这唯一的驻点就是所求的最大值（或最小值）。为了避免用较多的数学知识或数学上的推导，这里不追求数学上的完整性。

在解决实际问题时，当总体参数未知，需要通过样本来估计，我们仅对k2的情形加以说明。设样本分别为X1,X2,(1)(1)(2)(2)和X(1)X,Xn112,X(2)n2，则

n2n1n1X(1)n2X(2)(1)2)(1)(X(1)X(2)) X(2)X(X(X) XXXn1n2n1n2n1n2ˆn(X(1)X)(X(1)X)n(X(2)X)(X(2)X)n1n2(X(1)X(2))(X(1)X(2)) 那么 B12n1n2当μ1,μ2,,μk和Σ1,Σ2,,Σk均未知时，μ（1,2,,k）的估计同前，Σ（1,2,,k）的

ˆ1S， 1,2,,k估计为Σn1第五章 聚类分析

第一节 引言

“物以类聚，人以群分”。对事物进行分类，是人们认识事物的出发点，也是人们认识世界的一种重要方法。因此，分类学已成为人们认识世界的一门基础科学。

在生物、经济、社会、人口等领域的研究中，存在着大量量化分类研究。例如：在生物学中，为了研究生物的演变，生物学家需要根据各种生物不同的特征对生物进行分类。在经济研究中，为了研究不同地区城镇居民生活中的收入和消费情况，往往需要划分不同的类型去研究。在地质学中，为了研究矿物勘探，需要根据各种矿石的化学和物理性质和所含化学成分把它们归于不同的矿石类。在人口学研究中，需要构造人口生育分类模式、人口死亡分类状况，以此来研究人口的生育和死亡规律。

但历史上这些分类方法多半是人们主要依靠经验作定性分类，致使许多分类带有主观性和任意性，不能很好地揭示客观事物内在的本质差别与联系；特别是对于多因素、多指标的分类问题，定性分类的准确性不好把握。为了克服定性分类存在的不足，人们把数学方法引入分类中，形成了数值分类学。后来随着多元统计分析的发展，从数值分类学中逐渐分离出了聚类分析方法。随着计算机技术的不断发展，利用数学方法研究分类不仅非常必要而且完全可能，因此近年来，聚类分析的理论和应用得到了迅速的发展。

聚类分析就是分析如何对样品（或变量）进行量化分类的问题。通常聚类分析分为Q型聚类和R型聚类。Q型聚类是对样品进行分类处理，R型聚类是对变量进行分类处理。

第二节 相似性的量度

一、样品相似性的度量

 在聚类之前，要首先分析样品间的相似性。Q型聚类分析，常用距离来测度样品之间的相似程度。

每个样品有p个指标（变量）从不同方面描述其性质，形成一个p维的向量。如果把n个样品看成p维空间中的n个点，则两个样品间相似程度就可用p维空间中的两点距离公式来度量。两点

距离公式可以从不同角度进行定义，令dij 表示样品Xi与Xj的距离，存在以下的距离公式：

1．明考夫斯基距离dij(q)(k1pXikXjk)1/q明考夫斯基距离简称明氏距离，按的取值不同又可分成：

pq（1）绝对距离（q1） dij(1)XiXj k (5.2) k（2）欧氏距离（q2） dij(2)(XiXkk1kp121/2jk) (5.3)

（3）切比雪夫距离（q） dij()max XikXjk1kp欧氏距离是常用的距离，大家都比较熟悉，但是前面已经提到，在解决多元数据的分析问题时，欧氏距离就显示出了它的不足之处。一是它没有考虑到总体的变异对“距离”远近的影响，显然一个变异程度大的总体可能与更多样品近些，既使它们的欧氏距离不一定最近；另外，欧氏距离受变量的量纲影响，这对多元数据的处理是不利的。为了克服这方面的不足，可用“马氏距离”的概念。 2．马氏距离

设Xi与Xj是来自均值向量为 ，协方差为∑ =（＞0）的总体G中的p维样品，则两个样品间的马氏距

21离为dij(M)(XiXj)Σ(XiXj)

马氏距离又称为广义欧氏距离。显然，马氏距离与上述各种距离的主要不同就是它考虑了观测变量之间的相关性。如果各变量之间相互独立，即观测变量的协方差矩阵是对角矩阵，则马氏距离就退化为用各个观测指标的标准差的倒数作为权数的加权欧氏距离。马氏距离还考虑了观测变量之间的变异性，不再受各指标量纲的影响。将原始数据作线性变换后，马氏距离不变。

3．兰氏距离 dij(L)1pXikXjk pk1XikXjk它仅适用于一切Xij>0的情况，这个距离也可以克服各个指标之间量纲的影响。这是一个自身标准化的量，由于它对大的奇异值不敏感，它特别适合于高度偏倚的数据。虽然这个距离有助于克服明氏距离的第一个

缺点，但它也没有考虑指标之间的相关性。 4．距离选择的原则

一般说来，同一批数据采用不同的距离公式，会得到不同的分类结果。产生不同结果的原因，主要是由于不同的距离公式的侧重点和实际意义都有不同。因此我们在进行聚类分析时，应注意距离公式的选择。通常选择距离公式应注意遵循以下的基本原则：

（1）要考虑所选择的距离公式在实际应用中有明确的意义。如欧氏距离就有非常明确的空间距离概念。马氏距离有消除量纲影响的作用。

（2）要综合考虑对样本观测数据的预处理和将要采用的聚类分析方法。如在进行聚类分析之前已经对变量作了标准化处理，则通常就可采用欧氏距离。

（3）要考虑研究对象的特点和计算量的大小。样品间距离公式的选择是一个比较复杂且带有一定主观性的问题，我们应根据研究对象的特点不同做出具体分折。实际中，聚类分析前不妨试探性地多选择几个距离公式分别进行聚类，然后对聚类分析的结果进行对比分析，以确定最合适的距离测度方法。 二、变量相似性的度量

多元数据中的变量表现为向量形式，在几何上可用多维空间中的一个有向线段表示。在对多元数据进行分析时，相对于数据的大小，我们更多地对变量的变化趋势或方向感兴趣。因此，变量间的相似性，我们可以从它们的方向趋同性或“相关性”进行考察，从而得到“夹角余弦法”和“相关系数”两种度量方法。 1、夹角余弦

两变量Xi与Xj看作p维空间的两个向量，这两个向量间的夹角余弦可用下式进行计算

cosijXk1pk1pikXjkp 显然，∣cos  ij∣  1。

2(Xik)(X2jk)k12．相关系数

相关系数经常用来度量变量间的相似性。变量Xi与Xj的相关系数定义为

rij(Xk1pikXi)(XjkXj)p(Xk1p显然也有，∣rij∣  1。

ikXi)2(XjkXj)2k1

无论是夹角余弦还是相关系数，它们的绝对值都小于1，作为变量近似性的度量工具，我们把它们统记为cij。当∣cij∣= 1时，说明变量Xi与Xj完全相似；当∣cij∣近似于1时，说明变量Xi与Xj非常密切；当∣cij∣ = 0时，说明变量Xi与Xj完全不一样；当∣cij∣近似于0时，说明变量Xi与Xj差别很大。 据此，我们把比较相似的变量聚为一类，把不太相似的变量归到不同的类内。 在实际聚类过程中，为了计算方便，我们把变量间相似性的度量公式作一个变换为

dij = 1  ∣cij∣ (5.9)

或者 dij2=1  cij2 (5.10)

用表示变量间的距离远近，小则与先聚成一类，这比较符合人们的一般思维习惯。

第三节 系统聚类分析法

一、系统聚类的基本思想

系统聚类的基本思想是：距离相近的样品（或变量）先聚成类，距离相远的后聚成类，过程一直进行下去，每个样品（或变量）总能聚到合适的类中。系统聚类过程是：假设总共有n个样品（或变量），第一步将每个样品（或变量）独自聚成一类，共有n类；第二步根据所确定的样品（或变量）“距离”公式，把距离较近的两个样品（或变量）聚合为一类，其它的样品（或变量）仍各自聚为一类，共聚成n 1类；第三步将“距离”最近的两个类进一步聚成一类，共聚成n 2类；……，以上步骤一直进行下去，最后将所有的样品（或变量）全聚成一类。为了直观地反映以上的系统聚类过程，可以把整个分类系统画成一张谱系图。所以有时系统聚类也称为谱系分析。除系统聚类法外，还有有序聚类法、动态聚类法、图论聚类法、模糊聚类法等，限于篇幅，我们只介绍系统聚类方法。 二、类间距离与系统聚类法

在进行系统聚类之前，我们首先要定义类与类之间的距离，由类间距离定义的不同产生了不同的系统聚类法。常用的类间距离定义有8种之多，与之相应的系统聚类法也有8种，分别为最短距离法、最长距离法、中间距离法、重心法、类平均法、可变类平均法、可变法和离差平方和法。它们的归类步骤基本上是一致的，主要差异是类间距离的计算方法不同。以下用dij表示样品Xi与Xj之间距离，用Dij表示类Gi与Gj之间的距离。 1. 最短距离法

定义类与之间的距离为两类最近样品的距离，即为DijXiGi,XjGjmindij

设类与合并成一个新类记为，则任一类与的距离为

DkrXiGk,XjGrmindijmin{XiGk,XjGpmindij,minxiGk,xjGqdij}min{Dkp,Dkq}

最短距离法进行聚类分析的步骤如下：

（1）定义样品之间距离，计算样品的两两距离，得一距离 阵记为D（0） ，开始每个样品自成一类，

显然这时Dij = dij。

（2）找出距离最小元素，设为Dpq，则将Gp和Gq合并成一个新类，记为Gr，即Gr = ｛Gp，Gq｝。

（3）按（5.12）计算新类与其它类的距离。 （4）重复（2）、（3）两步，直到所有元素。并成一类为止。如果某一步距离最小的元素不止一个，则对应这些最小元素的类可以同时合并。

【例5.1】设有六个样品，每个只测量一个指标，分别是1，2，5，7，9，10，试用最短距离法将它们分类。

（1）样品采用绝对值距离，计算样品间的距离阵D（0） ，见表5.1

G1 G2 G3 G4 G1 0 1 4 6 G2 0 3 5 G3 0 2 G4 0 G5 G6 G5 G6 8 9 7 8 4 5 2 3 0 1 0 （2）D（0）中最小的元素是D12＝D56＝1，于是将G1和G2合并成G7，G5和G6合并成G8，并利用（5.12）式计算新类与其它类的距离D（1） ，见表5.2

G7 G3 G4 G8 G7 0 3 5 7 G3 0 2 4 G4 0 2 G8 0 （3）在D（1）中最小值是D34＝D48＝2，由于G4与G3合并，又与G8合并，因此G3、G4、G8合并成一个新类G9，其与其它类的距离D（2） ，见表5.3

G7 G9 G7 0 3 G9 0 （4）最后将G7和G9合并成G10，这时所有的六个样品聚为一类，其过程终止。

上述聚类的可视化过程见图5.1所示，横坐标的刻度表示并类的距离。这里我们应该注意，聚类的个数要以实际情况所定，其详细内容将在后面讨论。

2. 最长距离法

定义类Gi与Gj之间的距离为两类最远样品的距离，即为DpqXiGp,XjGqmaxdij (5.13)

最长距离法与最短距离法的并类步骤完全一样，也是将各样品先自成一类，然后将距离最小的两类合并。将类Gp与Gq合并为Gr，则任一类Gk与Gr的类间距离公式为

DkrXiGk,XjGrmaxmaxdij,maxdij}max{Dkp,Dkq} dijmax{XG,XGxG,xGikjpjikjq再找距离最小两类并类，直至所有的样品全归为一类为止。可以看出最长距离法与最短距离法只有两点不

同：一是类与类之间的距离定义不同；另一是计算新类与其它类的距离所用的公式不同。 3. 中间距离法

最短、最长距离定义表示都是极端情况，我们定义类间距离可以既不采用两类之间最近的距离也不采用两类之间最远的距离，而是采用介于两者之间的距离，称为中间距离法。

中间距离将类Gp与Gq类合并为类Gr，则任意的类Gk和Gr的距离公式为

2Dkr12122 (1／4    0) (5.15) DkpDkqDpq22设Dkq＞Dkp，如果采用最短距离法，则Dkr = Dkp，如果采用最长距离法，则Dkr = Dkq。如图5.2所示，(5.15)式就是取它们（最长距离与最短距离）的中间一点作为计算Dkr的根据。 特别当 =  1／4，它表示取中间点算距离，公式为Dkr121212DkpDkpDpq 224图5.2 中间距离法

4. 重心法

重心法定义类间距离为两类重心（各类样品的均值）的距离。重心指标对类有很好的代表性，但利用各样本的信息不充分。

设Gp与Gq分别有样品np，nq个，其重心分别为Xp和Xq，则Gp与Gq之间的距离定义为Xp和Xq之间的

2距离，这里我们用欧氏距离来表示，即Dpq(XpXq)(XpXq) (5.17)

设将Gp和Gq合并为Gr，则Gr内样品个数为nrnpnq，它的重心是Xr2重心是Xk，那么依据（5.17）式它与新类Gr的距离为Dkr1类Gk的(npXpnqXq)，

nrnpnr2Dkpnqnr2Dkqnpnqn2r2 （5.18） Dpq

这里我们应该注意，实际上（5.18）式表示的类Gk与新类Gr的距离为：

2Dkr(XkXr)(XkXr)[Xk11(npXpnqXq)][Xk(npXpnqXq)] nrnrXk2XknpnrXp2XknqnrXqXk

122(n2pXpXp2npnqXpXqnqXqXq)nrXk利用Xk1XknqXkXk)代入上式，有 (npXknrD2krnpnrXk2XkXpX(XkpXp)nqnrXk2XkXqXqXq)(Xknpnqnr(XpXp2XpXqXqXq)

npnrD2kpnqnrD2kqnpnqn2r2 （5.19） Dpq【例5.2】针对例5.1的数据，试用重心法将它们聚类。

（1）样品采用欧氏距离，计算样品间的平方距离阵D2（0），见表5.4所示。

G1 G2 G3 G4 G5 G6 G1 0 1 16 36 64 81 G2 0 9 25 49 64 G3 0 4 16 25 G4 0 4 9 G5 0 1 G6 0 （2）D2（0）中最小的元素是D212＝D256＝1，于是将G1和G2合并成G7，G5和G6合并成G8，并利用（5.18）式计算新类与其它类的距离得到距离阵D2（1） ，见表5.5：

G1 G2 G3 G4 2D37G1 0 12.25 30.25 64 G2 0 4 20.25 G3 0 6.25 G4 0 12121121111D31D32D12169112.25其它结果类似可以求得 22222222（3）在D2（1）中最小值是D234＝4，那么G3与G4合并一个新类G9，其与与其它类的距离D2（2） ，见表5.6：

G7 G9 G8 G7 G10 G7 0 20.25 64 G7 0 39.0625 G9 0 12.5 G8 0 G10 0 （4）在中最小值是＝12.5，那么与合并一个新类，其与与其它类的距离，见表5.7：

（5）最后将G7和G10合并成G11，这时所有的六个样品聚为一类，其过程终止。 上述重心法聚类的可视化过程见图5.3所示，横坐标的刻度表示并类的距离。 5．类平均法

类平均法定义类间距离平方为这两类元素两两之间距离平方的平均数，即为

2Dpq1npnqXiGpXjGjd2ij (5.20)

设聚类的某一步将Gp和Gq合并为Gr，则任一类类Gk与Gr的距离为：

2 Dkr1nknrXiGkXjGr2dijn2nq2122(dijdij)pDkpDkq （5.21） nknrXiGkXjGpnrnrXiGkXjGq类平均法的聚类过程与上述方法完全类似，这里就不在详述了。

6. 可变类平均法

由于类平均法中没有反映出Gp和Gq之间的距离Dpq的影响，因此将类平均法进一步推广，如果将Gp

2和Gq合并为新类Gr，类Gk与新并类Gr的距离公式为：Dkr(1)(npnr2Dkpnqnr22 Dkq)Dpq其中是可变的且 <1，称这种系统聚类法为可变类平均法。

7．可变法

针对于中间法而言，如果将中间法的前两项的系数也依赖于，那么，如果将Gp和Gq合并为新类Gr，

1222 （5.23） (DkpDkq)Dpq2nn1其中是可变的，且1。显然在可变类平均法中取pq，即为可变法。可变类平均法与可变法nrnr2的分类效果与的选择关系很大，在实际应用中常取负值。

2类Gk与新并类Gr的距离公式为：Dkr8. 离差平方和法

该方法是Ward提出来的，所以又称为Ward法。该方法的基本思想来自于方差分析，如果分类正确，同类样品的离差平方和应当较小，类与类的离差平方和较大。具体做法是先将n个样品各自成一类，然后每次缩小一类，每缩小一类，离差平方和就要增大，选择使方差增加最小的两类合并，直到所有的样品归为一类为止。

设将n个样品分成k类G1，G2，…，Gk，用Xit表示Gt中的第I个样品，nt表示Gt中样品的个数，

是Gt的重心，则Gt的样品离差平方和为St(Xt1ntitXt)(XitXt)

如果Gp和Gq合并为新类Gr类内离差平方和分别为

Sp(XipXp)(XipXp) Sq(XiqXq)(XiqXq) Sr(XirXr)(XirXr)

i1i1i1npnqnr它们反映了各自类内样品的分散程度，如果Gp和Gq这两类相距较近，则合并后所增加的离散平方和

2SrSpSq应较小；否则，应较大。于是定义Gp和Gq之间的平方距离为：DpqSrSpSq 2其中GrGpGq，可以证明类间距离的递推公式为DkrnknpnrnkD2kpnknqnrnk2Dkqnk2 Dpqnrnk这种系统聚类法称为离差平方和法或Ward方法。下面论证离差平方和法的距离递推（5.26）式。

Sr(XirXr)(XirXr)(XirXpXpXr)(XirXpXpXr)

i1i1nrnr(XirXp)(XirXp)(XirXp)(XpXr)(XpXr)(XirXp)(XpXr)(XpXr)i1i1i1i1nrnrnrnr

(XipXp)(XipXp)(XiqXp)(XiqXp)2(XpXr)(XirXp)nr(XpXr)(XpXr)i1i1i1npnqnr

Sp(XiqXqXqXp)(XiqXqXqXp)nr(XpXr)(XpXr) Sp(XiqXq)(XiqXq)nq(XpXq)(XpXq)nr(Xpi1nqi1nqnpXpnqXqnr)(XpnpXpnqXqnr)

SpSqnq(XpXq)(XpXq)n2pSpSqnq(XpXq)(XpXq)从而，由（5.25）式知

nrnqnp(XpXq)(XpXq)

nr(XpXq)(XpXq)

D2pqnqnpnr(XpXq)(XpXq) （5.27）

那么，由（5.27）式和（5.19）式，可以得到离差平方和法的平方距离的递推公式为：

2Dkrnrnk(XrXk)(XrXk)

nrnknrnknrnknqnpnqnp (XX)(XX)(XX)(XX)(XX)(XX)kpkpkqkqpqpq2nrnrnrnknpnrnknpnknknp(XkXp)(XkXp)nknqnrnknqnknknq(XkXq)(XkXq)

npnqnk(XpXq)(XpXq)nrnknrnknp2nknq2nk2 DkpDkqDpqnrnknrnknrnk上述八种系统聚类法的步骤完全一样，只是距离的递推公式不同。兰斯（Lance）和威廉姆斯（Williams）于1967年给出了一个统一的公式。

222222 (5.28) DkrpDkpqDkqDpqDkpDkq其中ap、aq、  、 是参数，不同的系统聚类法，它们取不同的数，详见表5.8。

这里应该注意，不同的聚类方法结果不一定完全相同，一般只是大致相似。如果有很大的差异，则应

该仔细考查，找到问题所在；另外，可将聚类结果与实际问题对照，看哪一个结果更符合经验。

表5.8 系统聚类法参数表

方 法 最短距离法 最长距离法 中间距离法 重心法 类平均法 可变类平均法 可变法 离差平方和法 p 1/2 1/2 1/2 q 1/2 1/2 1/2  0 0 -1/4  -1/2 1/2 0 0 0 0 0 0 npnrnpnrnqnrnqnrpq 0 (1)npnr (1)nqnr (1) (1) (1)/2 (1)/2 (npnk)(nrnk) (nqnk)(nrnk) nk(nknr) 第四节 K均值聚类分析

系统聚类法需要计算出不同样品或变量的距离，还要在聚类的每一步都要计算“类间距离”，相应的计算量自然比较大；特别是当样本的容量很大时，需要占据非常大的计算机内存空间，这给应用带来一定的困难。而K—均值法是一种快速聚类法，采用该方法得到的结果比较简单易懂，对计算机的性能要求不高，因此应用也比较广泛。

K均值法是麦奎因（MacQueen，1967）提出的，这种算法的基本思想是将每一个样品分配给最近中心（均值）的类中，具体的算法至少包括以下三个步骤：

1．将所有的样品分成K个初始类；2．通过欧氏距离将某个样品划入离中心最近的类中，并对获得样品与失去样品的类，重新计算中心坐标；3．重复步骤2，直到所有的样品都不能再分配时为止。

K均值法和系统聚类法一样，都是以距离的远近亲疏为标准进行聚类的，但是两者的不同之处也是明显的：系统聚类对不同的类数产生一系列的聚类结果，而K—均值法只能产生指定类数的聚类结果。具体类数的确定，离不开实践经验的积累；有时也可以借助系统聚类法以一部分样品为对象进行聚类，其结果作为K—均值法确定类数的参考。

下面通过一个具体问题说明K均值法的计算过程。

【例5.3】假定我们对A、B、C、D四个样品分别测量两个变量和得到结果见表5.9。

样品 X1 变量 X2 A 5 3 B -1 1 C 1 -2 D -3 -2 试将以上的样品聚成两类。 第一步：按要求取K=2，为了实施均值法聚类，我们将这些样品随意分成两类，比如（A、B）和（C、D），然后计算这两个聚类的中心坐标，见表5.10所示。

聚类 （A、B） （C、D） 中心坐标 X1 2 -1 X2 2 -2 表5.10中的中心坐标是通过原始数据计算得来的，比如（A、B）类的,X15(1)2等等。 2第二步：计算某个样品到各类中心的欧氏平方距离，然后将该样品分配给最近的一类。对于样品有变动的类，重新计算它们的中心坐标，为下一步聚类做准备。先计算A到两个类的平方距离：

d2(A,(AB))(52)2(32)210， d2(A,(CD))(51)2(32)261

由于A到（A、B）的距离小于到（C、D）的距离，因此A不用重新分配。计算B到两类的平方距离：

d2(B,(AB))(12)2(12)210，d2(B,(CD))(11)2(12)29

由于B到（A、B）的距离大于到（C、D）的距离，因此B要分配给（C、D）类，得到新的聚类是（A）和（B、C、D）。更新中心坐标如表5.11所示。

聚类 （A） （B、C、D） 中心坐标 X1 5 -1 样品到中心的距离平方 B C 40 41 4 5 X2 3 -1 第三步：再次检查每个样品，以决定是否需要重新分类。计算各样品到各中心的距离平方，得结果见表5.12。

聚类 （A） （B、C、D） A 0 52 D 89 5 到现在为止，每个样品都已经分配给距离中心最近的类，因此聚类过程到此结束。最终得到K=2的聚类结果是A独自成一类，B、C、D聚成一类。

第五节 有序样品的聚类分析法

以上的系统聚类和K—均值聚类中，样品的地位是彼此独立的，没有考虑样品的次序。但在实际应用中，有时样品的次序是不能变动的，这就产生了有序样品的聚类分析问题。例如对动植物按生长的年龄段进行分类，年龄的顺序是不能改变的，否则就没有实际意义了；又例如在地质勘探中，需要通过岩心了解地层结构，此时按深度顺序取样，样品的次序也不能打乱。

如果用X（1） ， X（2） ， …，X（n）表示n个有序的样品，则每一类必须是这样的形式，即X。在同一类中的（i） ， X（i+1)，…，X（j） ，其中1  r  n，且j  n，简记为Gi = ｛i，i+1，…，j｝样品是次序相邻的。这类问题称为有序样品的聚类分析。

一、有序样品可能的分类数目

n个有序样品分成k类，则一切可能的分法有Cn1种。实际上，n个有序样品共有（n 1）个间隔，分成k类相当于在这（n 1）个间隔中插入k 1根“棍子”。由于不考虑棍子的插入顺序，是一个组合问题，共有Cn1种插法。

k1k1

这就是n个有序样品分成k类的一切可能分法。因此，对于有限的n和k，有序样品的所有可能分类结果是有限的，可以在某种损失函数意义下，求得最优解。所以有序样品聚类分析又称为最优分割，该算法是费希尔（Fisher）最先提出来的，故也称之为费希尔最优求解法。 二、费希尔最优求解法 2．定义分类的损失函数。费希尔最优求解法定义的分类损失函数的思想类似于系统聚类分析中的Ward法，即要求分类后产生的离差平方和的增量最小。用b(n,k)表示将n个有序样品分为k类的某一种分法：

G1{i1,i11,,i21}，G2{i2,i21,,i31}，

其中i11i1k，Gk{ik,ik1,,n}

ikn。定义上述分类法的损失函数为

L[b(n,k)]D(it,it11) (5.31)

t1上式中的ik1n1。

对于固定的n和k，L[(b(n,k)]越小，表示各类的离差平方和越小，分类就是越有效的。因此，要求寻找一种分法b(n,k)，使分类的损失函数L[(b(n,k)]最小，这种最优分类法记为p(n,k)。

3．求最优分类法的递推公式。具体计算最优分类的过程是通过递推公式获得的。先考虑k2的情形对所有的j考虑使得,L[(b(n,2)]D(1,j)D(j,n)最小的j。得到最优分类p(n,2)：G1{1,2,,j*1}，G2{j*,,n}。

*图5.5 k2时的情形

进一步考虑对于k，求p(n,k)。

这里需要注意，若要寻找将n个样品分为k类的最优分割，则对于任意的j（k  j  n）,先将前面j 1个样品最优分割为k 1类，得到p（j  1，k  1），否则从j到n这最后一类就不可能构成k类的最优分割，

参见图5.6。再考虑使L[b（n，k）]最小的j＊，得到p（n，k）。

因此我们得到费希尔最优求解法的递推公式为

L[p(n,2)]min{D(1,j1)D(j,n)}2jn L[p(n,k)]min{L[p(j1,k1)]D(j,n)}kjn4．费希尔最优求解法的实际计算。从递推公式(5.32)可知，要得到分点jk，使得

L[(p(n,k)]L[p(jk1,k1)]D(jk,n)

从而获得第k类：Gk{jk,,n}，必须先计算jk1使得

L[(p(jk1,k1)]L[p(jk11,k2)]D(jk1,jk1)

从而获得第k1类：Gk1{jk1,,jk1}。

依此类推，…，要得到分点j3，使得L[(p(j41,3]L[p(j31,2)]D(j3,j41)

从而获得第3类：G3{j3,,j41}，必须先计算j2L[(p(j31,2)]min{D(1,j1)D(j,j31)}

2jj31G1,G2,,Gk。从而获得第2类：这时自然获得G1{1,,j21}。最后获得最优分割： G2{j2,,j31}。

因此，实际计算过程中是从计算j2开始的，一直到最后计算出jk为止。

总之，为了求最优解，主要是计算{D(i,j),1ijn}和{L[p(l,k)]，3ln,2kl,kn1}。 三、一个典型例子

【例5.4】为了了解儿童的生长发育规律，今随机抽样统计了男孩从出生到11岁每年平均增长的重量数据表5.13，试问男孩发育可分为几个阶段？ 2 3 年龄（岁） 1 增重（公斤） 9.3 1.8 1.9 在分析这是一个有序样品的聚类问题时，我们通过图形可以看到男孩增重随年龄顺序变化的规律，从图5.6中发现男孩发育确实可以分为几个阶段。

下面通过有序样品的聚类分析确定男孩发育分成几个阶段较合适。步骤如下：

（1）计算直径{D(i,j)}，结果如表5.14。例如计算D(1,2)，此类包含两个样品{9.3，1.8}，故有：

122XG(9.31.8)=5.55， D(1,2)(9.35.55)(1.85.55)=28.125

2其它依此计算，其结果见表5.14。 i j 2 3 1 28.125 37.007 2 0.005 3 4 5 6 7 8 9 10 4 5 6 7 8 9 10 11 42.208 45.992 49.128 51.100 51.529 51.980 52.029 52.182 0.020 0.088 0.232 0.280 0.417 0.467 0.802 0.909 0.020 0.080 0.200 0.232 0.393 0.454 0.800 0.909 0.020 0.080 0.088 0.308 0.393 0.774 0.895 0.020 0.020 0.290 0.388 0.773 0.889 0.005 0.287 0.370 0.708 0.793 0.180 0.207 0.420 0.452 0.005 0.087 0.088 0.080 0.080 0.020 （2）计算最小分类损失函数{L[p(l,k)]}，结果如表5.15。 k l 3 4 5 6 7 8 9 10 11 2 0.005(2) 0.020(2) 0.088(2) 0.232(2) 0.280(2) 0.417(2) 0.469(2) 0.802(2) 0.909(2) 3 0.005(4) 0.020(5) 0.040(5) 0.040(5) 0.280(8) 0.285(8) 0.367(8) 0.368(8) 4 0.005(5) 0.020(6) 0.025(6) 0.040(8) 0.045(8) 0.127(8) 0.128(8) 5 0.005(6) 0.010(6) 0.025(8) 0.030(8) 0.045(10) 0.065(10) 6 0.005(6) 0.010(8) 0.015(8) 0.030(10) 0.045(11) 7 0.005(8) 0.010(3) 0.015(10) 0.030(11) 8 0.005(8) 0.010(10) 0.015(11) 9 0.005(8) 0.010(11) 10 0.005(11) 首先计算{L[p(l,2)],3l11}（即表中的k2列），例如计算： L[p(3,2)]min{D(1,j1)D(j,3)}min{D(1,1)D(1,3), D(1,2)D(3,3)}

2j3min{00.005, 28.1250}0.005

极小值是在j2处达到，故记L[p(3,2)]0.005(2)，其它类似计算。 再计算{L[p(l,3)],4l11}（即表中的k3列），例如计算：

L[p(4,3)]min{L[p(2,2)]D(3,4),L[p(3,2)]D(4,4)} min{00.02, 0.0050}0.005(4) 表5.15中其它数值同样计算，括弧内的数字表示最优分割处的序号。

（3）分类个数的确定。如果能从生理角度事先确定k当然最好；有时不能事先确定k时，可以从L[p（l，k）]随k的变化趋势图中找到拐点处，作为确定k的根据。当曲线拐点很平缓时，可选择的k很多，这时需要用其它的办法来确定，比如均方比和特征根法，限于篇幅此略，有兴趣的读者可以查看其它资料。

本例从表5.15中的最后一行可以看出k =3，4处有拐点，即分成3类或4类都是较合适的，从图5.8中可以更明显看出这一点。

（4）求最优分类。例如我们把儿童生长分成4个阶段，即可查表5.15中k4例的最后一行（即l11行）得L[p(11,4)]0.128(8)，说明最优损失函数值为0.128，最后的最优分割在第8个元素处，因此G4{8~11}或G4{2.0,1.9,2.3,2.1}。进一步从表中查L[p(7,3)]0.040(5)，因此G3{5~7}或G3{1.5,1.3,1.4}，再从表中查得L[p(4,2)]0.020(2)最后G2{2~4}或G2{1.8,1.9,1.7}，剩下的G1{9.3}。

第六章 主成分分析

第一节 引言

多元统计分析处理的是多变量（多指标）问题。由于变量较多，增加了分析问题的复杂性。但在实际问题中，变量之间可能存在一定的相关性，因此，多变量中可能存在信息的重叠。人们自然希望通过克服相关性、重叠性，用较少的变量来代替原来较多的变量，而这种代替可以反映原来多个变量的大部分信息，这实际上是一种“降维”的思想。

主成分分析也称主分量分析，是由Hotelling于1933年首先提出的。由于多个变量之间往往存在着一定程度的相关性。人们自然希望通过线性组合的方式，从这些指标中尽可能快地提取信息。当第一个线性组合不能提取更多的信息时，再考虑用第二个线性组合继续这个快速提取的过程，……，直到所提取的信息与原指标相差不多时为止。这就是主成分分析的思想。一般说来，在主成分分析适用的场合，用较少的主成分就可以得到较多的信息量。以各个主成分为分量，就得到一个更低维的随机向量；因此，通过主成分既可以降低数据“维数”又保留了原数据的大部分信息。

我们知道，当一个变量只取一个数据时，这个变量（数据）提供的信息量是非常有限的，当这个变量取一系列不同数据时，我们可以从中读出最大值、最小值、平均数等信息。变量的变异性越大，说明它对各种场景的“遍历性”越强，提供的信息就更加充分，信息量就越大。主成分分析中的信息，就是指标的变异性，用标准差或方差表示它。

主成分分析的数学模型是，设p个变量构成的p维随机向量为X = （X1，…，Xp）′。对X作正交变换，令Y = T′X，其中T为正交阵，要求Y的各分量是不相关的，并且Y的第一个分量的方差是最大的，第二个分量的方差次之，……，等等。为了保持信息不丢失，Y的各分量方差和与X的各分量方差和相等。

第二节 主成分的几何意义及数学推导

一、主成分的几何意义

主成分分析数学模型中的正交变换，在几何上就是作一个坐标旋转。因此，主成分分析在二维空间中有明显的几何意义。假设共有n个样品，每个样品都测量了两个指标（X1，X2），它们大致分布在一个椭圆内如图6.1所示。事实上，散点的分布总有可能沿着某一个方向略显扩张，这个方向就把它看作椭圆的长轴方向。显然，在坐标系x1Ox2中，单独看这n个点的分量X1和X2，它们沿着x1方向和x2方向都具有较大的离散性，其离散的程度可以分别用的X1方差和X2的方差测定。如果仅考虑X1或X2中的任何一个分量，那么包含在另一分量中的信息将会损失，因此，直接舍弃某个分量不是“降维”的有效办法。 如果我们将该坐标系按逆时针方向旋转某个角度变成新坐标系y1Oy2，这里y1是椭圆的长轴方向，y2是

YX1cosX2sin椭圆的短轴方向。旋转公式为1 (6.1）

YXsinXcos212我们看到新变量Y1和Y2是原变量X1和X2的线性组合，它的矩阵表示形式为：

Y1cosYsin2sinX1TX (6.2） Xcos21其中，T为旋转变换矩阵，它是正交矩阵，即有TT或TTI。

易见，n个点在新坐标系下的坐标Y1和Y2几乎不相关。称它们为原始变量X1和X2的综合变量，n个点y1在轴上的方差达到最大，即在此方向上包含了有关n个样品的最大量信息。因此，欲将二维空间的点投影到某个一维方向上，则选择y1轴方向能使信息的损失最小。我们称Y1为第一主成分，称Y2为

第二主成分。第一主成分的效果与椭圆的形状有很大的关系，椭圆越是扁平，n个点在y1轴上的方差就相对越大，在y2轴上的方差就相对越小，用第一主成分代替所有样品所造成的信息损失也就越小。 考虑两种极端的情形：

一种椭圆的长轴与短轴的长度相等，即椭圆变成圆，第一主成分只含有二维空间点的约一半信息，若仅用这一个综合变量，则将损失约50％的信息，这显然是不可取的。造成它的原因是，原始变量X1和X2的相关程度几乎为零，也就是说，它们所包含的信息几乎不重迭，因此无法用一个一维的综合变量来代替。

另一种是椭圆扁平到了极限，变成y1轴上的一条线，第一主成分包含有二维空间点的全部信息，仅用这一个综合变量代替原始数据不会有任何的信息损失，此时的主成分分析效果是非常理想的，其原因是，第二主成分不包含任何信息，舍弃它当然没有信息损失。 二、主成分的数学推导

设X(X1,,Xp)为一个p维随机向量，并假定存在二阶矩，其均值向量与协差阵分别记为： μE(X), ΣD(X) (6.3） 考虑如下的线性变换

YtXtXtXTX111122pp111Yt2X11tX222tpXp2TX2 2 (6.4）

tppXpTpXYptp1X1tpX22用矩阵表示为 YTX，其中Y(Y1,Y2,Yp)，T(T1,T2,,Tp)。

,Xp的信息，

我们希望寻找一组新的变量Y1,（mp），这组新的变量要求充分地反映原变量X1,,Ym,Ym有

而且相互独立。这里我们应该注意到，对于Y1,D(Yi)D(TiX)TiD(X)TiTiΣTi i1,2,m,

Cov(Yi,Yk)Cov(TiX,TkX)TiCov(X,X)TkTiΣTk i,k1,2,m,

这样，我们所要解决的问题就转化为，在新的变量Y1,,Ym相互独立的条件下，求Ti使得D(Yi)TiΣTi，

i1,2,,m,达到最大。

我们下面将借助投影寻踪（Projection Pursuit）的思想来解决这一问题。首先应该注意到，使得D(Yi)达到最大的线性组合，显然用常数乘以Ti后，D(Yi)也随之增大，为了消除这种不确定性，不妨假设Ti满足TiTi1或者T1。那么，问题可以更加明确。

第一主成分为，满足T1T11，使得D(Y1)T1ΣT1达到最大的Y1T1X。 (2Y,1)CovT(第二主成分为，满足T2T21，且CovY2XT,1)X0，使得D(Y2)T2ΣT2达到最大的Y2T2X。

TkX,TiX)0（ik）一般情形，第k主成分为，满足TkTk1，且Cov(Yk,Yi)Cov(，使得D(Yk)TkΣTk达

到最大的YkTkX。

求第一主成分，构造目标函数为：1(T1,)T1ΣT1(T1T11) (6.5)

对目标函数1(T1,)求导数12ΣT12T10(6.6)，即(ΣI)T10 (6.7)

T1由6.7式两边左乘T1得到T1ΣT1 (6.8)

由于X的协差阵Σ为非负定的，其特征方程(6.7)的根均大于零，不妨设12道Y1的方差为。那么，Y1的最大方差值为1，其相应的单位化特征向量为T1。

在求第二主成分之前，我们首先明确，由(6.6)知Cov(Y2,Y1)T2ΣT1 T2T1。那么，如果Y2与Y1相互独立，

p0。由(6.8)知

即有T2T10或T1T20。这时，我们可以构造求第二主成分的目标函数，即

2(T2,,)T2ΣT2(T2T21)2(T1T2) (6.9)

对目标函数2(T2,,)求导数有：22ΣT22T22T10 (6.10)

T2用T1左乘(6.10)式有 T1ΣT2T1T2T1T10

由于T1ΣT20，T1T20，那么，T1T10，即有0。从而(ΣI)T20 (6.11) 而且 T2ΣT2 (6.12)

这样说明，如果X的协差阵Σ的特征根为12大特征根2，其相应的单位化的特征向量为T2。

针对一般情形，第k主成分应该是在TkTk1且TkTi0或TiTk0（ik）的条件下，使得D(Yk)TkΣTk达到最大的YkTkX。这样我们构造目标函数为

p0。由(6.12)知道Y2的最大方差值为第二

k(Tk,,i)TkΣTk(TkTk1)2i(TiTk) (6.13)

i1k1k1k对目标函数k(Tk,,i)求导数有：2ΣTk2Tk2iTi0 (6.14)

i1Tk用Ti左乘(6.14)式有 TiΣTkTiTkTi(iTi)0

i1k1即有iTiTi0，那么，i0（i1,2,。从而(ΣI)Tk0 (6.15)，而且TkΣTk (6.16) k1）

对于X的协差阵Σ的特征根12p0。由(6.15)和(6.16)知道Yk的最大方差值为第k大特

,Xp)的协差阵为Σ，其特征根为

征根k，其相应的单位化的特征向量为Tk。综上所述，设X(X1,12Y2T2X，

p0，相应的单位化的特征向量为T1,T2,,Tp。那么，由此所确定的主成分为Y1T1X，

，YmTmX，其方差分别为Σ的特征根。

第三节 主成分的性质

一、主成分的一般性质 设Y(Y1,Y2,1,Yp)是X的主成分，由Σ的所有特征根构成的对角阵为Λ00 (6.17） p主成分可表示为YTX (6.18）

性质1 主成分的协方差矩阵是对角阵。

证明：实际上，由（6.3）式知E(Y)E(TX)Tμ ，D(Y)TD(X)TTΣTΛ (6.19） 性质2 主成分的总方差等于原始变量的总方差。

证明：由矩阵“迹”的性质知tr(Λ)tr(TΣT)tr(ΣTT)tr(Σ)

所以

ii1i1ppii (6.20） 或 D(Yi)D(Xi) (6.21）

i1i1ppktki (6.22）并称之为因子负荷量（或因子性质3 主成分Yk与原始变量Xi的相关系数为(Yk,Xi)ii载荷量）。

证明：事实上(Yk,Xi)其中的ei(0,Cov(Yk,Xi)D(Yk)D(Xi)Cov(TkX,eiX)kii

,0,1,0,,0)，它是除第i个元素为1外其他元素均为0的单位向量。而

Cov(TkX,eiX)TkΣeiei(ΣTk)ei(kTk)keTikktki

所以(Yk,Xi)性质4

ktki。 iii1p2（k1,2,,p）。 (Yk,Xi)iik，

证明：只须将(6.22）代入左边式子整理化简即可。

二、主成分的方差贡献率

由性质2可以看出，主成分分析把p个原始变量X1,X2,p,Xp的总方差tr(Σ)分解成了p个相互独立

的变量Y1,Y2,,Yp的方差之和k。主成分分析的目的是减少变量的个数，所以一般不会使用所有p个

k1主成分的，忽略一些带有较小方差的主成分将不会给总方差带来太大的影响。这里我们称kk(6.23）为第

k1pk

k个主成分Yk的贡献率。第一主成分的贡献率最大，这表明Y1T1X综合原始变量

,Yp的综合能力依次递减。若只取m(p)个主成分，则称

,Ym的累计贡献率，累计贡献率表明Y1,,Ym综合X1,X2,pX1,X2,mkk1m,Xp的能力最强，而Y2,Y3,k1k (6.24）为主成分Y1,,Xp的

能力。通常取m，使得累计贡献率达到一个较高的百分数（如85％以上）。

第四节 主成分方法应用中应注意的问题

一、实际应用中主成分分析的出发点

我们前面讨论的主成分计算是从协方差矩阵Σ出发的，其结果受变量单位的影响。不同的变量往往有不同的单位，对同一变量单位的改变会产生不同的主成分，主成分倾向于多归纳方差大的变量的信息，对于方差小的变量就可能体现得不够，也存在“大数吃小数”的问题。为使主成分分析能够均等地对待每一个原始变量，消除由于单位的不同可能带来的影响，我们常常将各原始变量作标准化处理，即令

XE(Xi)Xi*i i1,,p (6.25）

D(Xi)显然，X*(X1*,,X*实际应用中，X的相关系数矩阵R可以p)的协方差矩阵就是X的相关系数矩阵R。

通过（2.13）式，利用样本数据来估计。

这里我们需要进一步强调的是，从相关阵求得的主成分与协差阵求得的主成分一般情况是不相同的。实际表明，这种差异有时很大。我们认为，如果各指标之间的数量级相差悬殊，特别是各指标有不同的物理量纲的话，较为合理的做法是使用R代替∑。对于研究经济问题所涉及的变量单位大都不统一，采用R代替∑后，可以看作是用标准化的数据做分析，这样使得主成分有现实经济意义，不仅便于剖析实际问题，又可以避免突出数值大的变量。

因此，在实际应用中，主成分分析的具体步骤可以归纳为： 1. 将原始数据标准化；

2. 建立变量的相关系数阵；

****3. 求R的特征根为1*p0，相应的特征向量为T1,T2,,Tp；

4. 由累积方差贡献率确定主成分的个数（m），并写出主成分为Yi(Ti*)X， i1,2,,m

二、如何利用主成分分析进行综合评价

人们在对某个单位或某个系统进行综合评价时都会遇到如何选择评价指标体系和如何对这些指标进行综合的困难。一般情况下，选择评价指标体系后通过对各指标加权的办法来进行综合。但是，如何对指标加权是一项具有挑战性的工作。指标加权的依据是指标的重要性，指标在评价中的重要性判断难免带有一定的主观性，这影响了综合评价的客观性和准确性。由于主成分分析能从选定的指标体系中归纳出大部分信息，根据主成分提供的信息进行综合评价，不失为一个可行的选择。这个方法是根据指标间的相对重要性进行客观加权，可以避免综合评价者的主观影响，在实际应用中越来越受到人们的重视。

对主成分进行加权综合。我们利用主成分进行综合评价时，主要是将原有的信息进行综合，因此，要充分的利用原始变量提供的信息。将主成分的权数根据它们的方差贡献率来确定，因为方差贡献率反映了各个主成分的信息含量多少。

设Y1,Y2,,Yp是所求出的p个主成分，它们的特征根分别是1,2,,p，将特征根“归一化”即有

wii i1,2,p,

i1mi记为W(w1,w2,wp)，由YTX，构造综合评价函数为

Zw1Y1w2Y2wpYpWYWTX(TW)X (6.26）

*令TWwk，并代入(6.26）式，有Z(w*)X (6.27) 1这里我们应该注意，从本质上说综合评价函数是对原始指标的线性综合，从计算主成分到对之加权，经过两次线性运算后得到综合评价函数。