本节书摘来自华章计算机《数学建模：基于R》一书中的第2章，第2.3节,作者：薛 毅 更多章节内容可以访问云栖社区“华章计算机”公众号查看。 2.3 判别分析判别分析是用以判别个体所属群体的一种统计方法,它产生于20世纪30年代.近年来,在许多现代自然科学的各个分支和技术部门中得到广泛的应用.例如,利用计算机对一个人是否有心脏病进行诊断时,可以取一批没有心脏病的人,测其p个指标的数据,然后再取一批已知患有心脏病的人,同样也测得p个相同指标的数据,利用这些数据建立一个判别函数,并求出相应的临界值,这时对于需要进行诊断的人,也同样测其p个指标的数据,将其代入判别函数,求得判别得分,再依判别临界值,就可以判断此人是属于有心脏病的那一群体,还是属于没有心脏病的那一群体.又如,在考古学中,对化石及文物年代的判断；在地质学中,判断是有矿还是无矿；在质量管理中,判断某种产品是合格品,还是不合格品;在植物学中,对于新发现的一种植物,判断其属于哪一科.总之,判别分析方法在很多学科中都有着广泛的应用.2.3.1 判别分析的基本原理所谓判别问题,就是将p维Euclid空间Rp划分成k个互不相交的区域R1,R2,…,Rk,即Ri∩Rj=,i≠j,i,j=1,2,…,k,∪kj=1Rj=Rp.当x∈Ri,i=1,2,…,k,就判定x属于总体Xi,i=1,2,…,k.特别地,当k=2时,就是两个总体的判别问题.关于判别方法,有距离判别、Fisher判别和Bayes判别等,这里以距离判别为例,说明判别分析的基本原理.距离判别的基本原理是:计算待判样本与总体之间的距离,样本距哪一类总体的距离最近,就判断待判样本属于哪一类.但在判别分析中,使用Euclid距离是不合适的,请见图2.6所示的例子.图2.6 不同方差的正态分布函数为简单起见,考虑p=1的情况.设X~N(0,1),Y~N(4,22),绘出相应的概率密度曲线,如图2.6所示.考虑图中的A点,A点距X的均值μ1=0较近,距Y的均值μ2=4较远.但从概率角度来分析问题,情况并非如此.经计算,A点的x值为1.66,也就是说,A点距μ1=0是1.66σ1,而A点距μ2=4却只有1.17σ2,因此,从概率分布的角度来讲,应该认为A点距μ2更近一点.所以,在定义距离时,要考虑随机变量方差的信息.设x,y是从均值为μ,协方差阵为Σ的总体X中抽取的两个样本,则总体X内两点x与y的Mahalanobis距离（马哈拉诺比斯距离,简称为马氏距离）定义为d(x,y)=(x－y)TΣ－1(x－y)（2.26）定义样本x与总体X的Mahalanobis距离为d(x,X)=(x－μ)TΣ－1(x－μ)（2.27）在这里,讨论两个总体的距离判别,分别讨论两总体协方差阵相同和协方差阵不同的情况.设总体X1和X2的均值向量分别为μ1和μ2,协方差阵分别为Σ1和Σ2,给定一个样本x,要判断x来自哪一个总体.如果两个总体X1和X2的协方差矩阵是相同的,即μ1≠μ2, Σ1=Σ2=Σ要判断x属于哪一个总体,需要计算x到总体X1和X2的Mahalanobis距离的平方d2(x,X1)和d2(x,X2),然后进行比较,若d2(x,X1)≤d2(x,X2),则判定x属于X1;否则判定x来自X2.由此得到如下判别准则R1={xd2(x,X1)≤d2(x,X2)}, R2={xd2(x,X1)>d2(x,X2)}（2.28）现在引进判别函数的表达式,考虑d2(x,X1)与d2(x,X2)之间的关系,有d=d2(x,X2)-d2(x,X1)=(x-μ2)TΣ-1(x-μ2)-(x-μ1)TΣ-1(x-μ1)=(xTΣ-1x-2xTΣ-1μ2+μT2Σ-1μ2)-(xTΣ-1x-2xTΣ-1μ1+μT1Σ-1μ1)=2xTΣ-1(μ1-μ2)+(μ1+μ2)TΣ-1（μ2-μ1）=2x-μ1+μ22TΣ-1(μ1-μ2)=2(x-μ)TΣ-1(μ1-μ2)(2.29)其中μ=μ1+μ22是两个总体均值的平均.令w(x)=(x－μ)TΣ－1(μ1－μ2)（2.30）称w(x)为两总体的距离判别函数.如果两个总体X1和X2的协方差阵是不同的,即μ1≠μ2, Σ1≠Σ2对于样本x,在协方差阵不同的情况下,判别函数为w(x)=(x－μ2)TΣ－12(x－μ2)－(x－μ1)TΣ－11(x－μ1)（2.31）无论是协方差矩阵相同的情况,还是不同的情况,判别准则(2.28)可改为R1={xw(x)≥0}, R2={xw(x)<0}（2.32）在实际计算中,总体的均值与协方差阵是未知的,因此,需要用样本的均值和样本的协方差阵来代替.设x(1)1,x(1)2,…,x(1)n1是来自总体X1的n1个样本,x(2)1,x(2)2,…,x(2)n2是来自总体X2的n2个样本,则样本均值为μi=x(i)=1ni∑nij=1x(i)j, i=1,2(2.33)如果认为两个总体的协方差矩阵相同,则样本协方差阵定义为Σ=1n1+n2－2∑2i=1∑nij=1(x(i)j－x(i))(x(i)j－x(i))T(2.34)对应的判别函数(2.30)改为w(x)=(x－x)TΣ－1(x(1)－x(2))（2.35）其中x=x(1)+x(2)2.如果认为两个总体的协方差矩阵不同,则样本协方差阵定义为Σi=1ni－1∑nij=1(x(i)j－x(i))(x(i)j－x(i))T, i=1,2（2.36）对应的判别函数(2.31)改为w(x)=(x－x(2))TΣ－12(x－x(2))－(x－x(1))TΣ－11(x－x(1))（2.37）判别准则(2.32)改为R1={xw(x)≥0}, R2={xw(x)<0}（2.38）2.3.2 判别分析的计算如果两个总体的协方差矩阵是相同时,则判别函数(2.30)或判别函数(2.35)是线性函数;当两个总体的协方差矩阵不同时,则判别函数(2.31})或判别函数(2.37)是二次函数.所以,距离判别属于线性判别或者二次判别.同样,可以证明:Fisher判别属于线性判别,Bayes判别本质属于线性判别或者二次判别.R并没有单独提供这三种判别,而是将判别分析的方法综合在一起,分别给出线性判别函数——lda()函数和二次判别函数——qda()函数.在使用这两个函数之前,需要加载MASS程序包,或使用命令library(MASS).两个函数的使用格式基本相同,一种是公式形式,其使用格式为lda(formula, data, ..., subset, na.action)qda(formula, data, ..., subset, na.action)参数formula为公式,形如groups ~ x1 + x2 + ....data为数据构成的数据框.subset为可选择向量,表示观察值的子集.na.action为函数,表示当数据中出现缺失数据(NA)的处理方法.另一种是矩阵或数据框形式,其使用格式为lda(x, grouping, prior = proportions, tol = 1.0e-4, method, CV = FALSE, nu, ...)qda(x, grouping, prior = proportions, method, CV = FALSE, nu, ...)参数x为矩阵或数据框,或者为包含解释变量的矩阵.grouping为指定样本属于哪一类的因子向量.prior为各类的先验概率,默认值是已有训练样本的计算结果.tol控制精度,当矩阵的方差小于tol^2时,表明矩阵奇异.method为字符串,表示估计方法,取"moment"表示均值和方差的标准估计,取"mle"表示极大似然估计,取"mve"表示使用cov.mve作估计,取"t"表示基于t分布的稳健估计.CV为逻辑变量,表示在返回值中是否包含留一交叉验证的结果,默认值为FALSE. nu为method = "t"的自由度....为附加参数.lda()函数的返回值是一个列表,其元素有:调用方法、先验概率、每一类样本的均值和线性判别系数等.qda()函数的返回值也是一个列表,其元素基本上与lda()函数相同.在lda()函数或qda()函数计算后,还需要用predict()函数给出待测样本的预测值或者训练样本的回代值.例2.15 表2.12是某气象站监测前14年气象的实际资料,有两项综合预报因子（气象含义略）,其中有春旱的是6个年份的资料,无春旱的是8个年份的资料.今年观测到两个指标的数据为(23.5,-1.6),试用lda()函数和qda()函数对数据作判别分析,并预报今年是否有春旱.
lda.sol <- lda(train, sp)计算待测样本的预测值和训练样本的回代值：> predict(lda.sol, tst)$class［1］NoLevels:Have No
table(sp,predict(lda.sol)$class)
sp Have No Have 5 1 No 0 8预测结果是无春旱,回代计算中,本有6个有春旱的年份,只判对了5个.使用table()函数是为了以表格形式列出.再看二次判别的计算结果：> qda.sol<-qda(train, sp); predict(qda.sol, tst)$class
［1］Have Levels:Have No >table(sp,predict(qda.sol)$class) sp Have No Have60No08这次得到的结果是有春旱.到底是有春旱还是无春旱呢?从回代结果来看,可能有春旱更合理一些,因为二次判别的回代正确率是100%.例2.16（Fisher Iris数据） Iris数据有4个属性,萼片的长度、萼片的宽度、花瓣长度和花瓣的宽度（单位：cm）.数据共150个样本,分为三类,前50个数据为第一类——Setosa,中间的50个数据为第二类——Versicolor,最后50个数据为第三类——Virginica.数据格式如表2.13所示.试用lda()函数（或qda()函数）对Iris数据作判别分析.表2.13 Fisher Iris数据序号萼片长度萼片宽度花瓣长度花瓣宽度种类15.13.51.40.2setosa24.93.01.40.2setosa（续）序号萼片长度萼片宽度花瓣长度花瓣宽度种类505.03.31.40.2setosa517.03.24.71.4versicolor526.43.24.51.5versicolor1005.72.84.11.3versicolor1016.33.36.02.5virginica1025.82.75.11.9virginica1505.93.05.11.8virginica解 R软件中提供了Iris数据(数据框iris),数据的前4列是数据的4个属性,第5列标明数据属于哪一类.用lda()函数作判别分析,读者可仿照本例完成qda()函数的判别分析.在150个样本中随机选取100个作为训练样本,余下的50个作为待测样本,先验概率各为1/3.程序(程序名:exam0216.R)为train <- sample(1:150, 100)z <- lda(Species ~ ., iris, prior = c(1,1,1)/3, subset = train)(class<-predict(z, iris［-train, ］)$class)程序中sample()函数为抽取样本,subset = train表示选择抽取的样本作为训练子集,iris［-train, ］表示在预测函数中使用其余的样本.计算结果略.看一下预测结果的准确性：> sum(class==iris$Species［-train］)［1］ 50全部正确.
相关资源：数学建模判别分析