[精品课程]chap7判别分析.docx

第七章 判别分析判别问题在日常生活中很常见：在河里逮住一条怪鱼，要判别它是哪种鱼，这是判别问题；一个非婚生子女，判别他的父亲是否某人，这是判别问题；从今年春天气象资料预报今年夏天是偏热、偏冷还是正常，也是判别问题。判别问题的特点是，事先已知若干个类，今有一个样本，需要判别它属于哪一类。解决判别问题有多种方法，例如观测鱼的特点，dna鉴定等。多元统计分析的判别分析方法则是从可观测的多个指标出发。例如已知某地第二年汛期雨量偏大、偏小与正常这三种情况取决于第一年10个气象指标。我们想建立从10个气象指标预报第二年雨量是否正常的公式，于是我们把第二年雨量分为三类，建立判别法则：用前一年10个气象指标判定第二年属于哪一类。这10个气象指标也可称为10个预报因子，第二年雨量类型作为预报对象。从预报因子预报对象的类型来看，判别分析的预报因子是数值变量，预报对象是属性变量。从多元统计分析的角度来看，每个预报因子作为一个随机变量，所有预报因子作为随机向量，每一类作为一个总体，由均值的不同确定不同类。7.1数学模型与判别方法例7.1 已知长江中游6月份降水分为三级：偏少、偏多、正常，分别用1，2，3来表示这三个母体（用变量species=1,2,3表示）。现利用气象因子x1,x2,x3,x4为预报因子（随机向量）作判别，1951-1975年观测数据(已分好类)见表7.1。1976-1979 x1,x2,x3,x4年数据已知，不知分类（待分类）见表7.2。试利用1951-1975年资料所得判别方法判定1976-1979年长江中游6月份降水级别。表7.1 1951-1975长江中游6月份降水及气象因子x1,x2,x3,x4的值yearx1x2x3x4species19510.5882.044.040.6119520.4083.018.043.0219530.5585.036.030.7219540.4085.036.040.7219550.4888.049.043.0219560.4182.035.078.6319570.6580.029.033.2119580.4582.032.033.1319590.3981.027.046.5319600.3485.028.041.7319610.4284.038.020.4319620.5286.038.00.2119630.4688.025.056.7219640.4883.046.013.6119650.5384.041.032.3119660.6581.031.028.9119670.6683.038.046.6119680.5380.042.093.1319690.5685.018.016.3319700.4583.037.023.9319710.3480.042.026.3319720.4179.038.040.8319730.5383.023.061.3319740.4884.019.023.2219750.3085.027.017.52表7.2 1976-1979长江中游x1,x2,x3,x4的值yearx1x2x3x419760.4281.021.052.219770.5281.038.045.819780.3682.034.034.919790.4384.034.060.5多元统计判别把预报因子看成随机向量，每一类看成一个总体。判别分析的数学模型是：设有k个总体，（这k个总体通常称为k个类）。它们分别服从多元正态分布，其中，是已知或未知向量，.是已知或未知矩阵，每次观测是一个样品。今有一观测向量x（样品），需要判定x来自k个总体中的哪一个（属于哪一类）。 实际问题中， 总是未知的，但我们知道每个总体的若干个样品：第一类样品，第k类样品，。从而可以分别计算每类样本均值（也称为j类均值），和样本方差阵（也称为i类样本方差） （7.1）用它们代替各总体的均值和方差阵。判别分析的方法很多，诸如距离判别法，最大概率判别法等，由于，都是随机变量的观测值因而存在误判问题。我们这儿对误判问题不作讨论，仅介绍这些方法。1 欧氏距离判别:计算x与第j个总体均值的欧氏距离，也称为x与第j个总体的欧氏距离；判定x属于距离最小的一类。一般情况总是未知的，用样本均值 代替，即计算x与第j个总体样本均值的欧氏距离，判定x属于距离最小的一类，即x来自距离最小的总体。欧氏距离判别计算简单，不用sas软件也能很快编程(如c语言编程)，非常方便。但由于欧氏距离与样品各分量单位有关，例如长度分量由cm为单位，改为mm为单位，长度分量数值扩大10倍，其影响扩大10倍，而重量分量由g改为kg，重量分量的影响缩小1000倍；而且有些相互联系的分量都作为一个加数计算，出现了重复计算，不合理，以下马氏距离判别法能避免这些缺点。 2.马氏距离判别：计算x与第i个总体均值的马氏距离，也称为x与第i类的距离；判定x属于距离最小的一类。一般情况总是未知的，用样本均值代替，即计算x与第i个总体样本均值的马氏距离。由于平方根函数是增函数，通常考虑而不是其平方根；判定x属于马式距离最小的类，即x来自距离最小的总体。马氏距离判别法的优点是：马式距离与单位无关，且相关的随机变量不被统计。简证如下：将单位改变或改变其相关性，这相当于x做变换y=cx,（无妨设c是可逆阵），则每类样本均值变为，样本方差阵变为=。y与第j个总体均值的马式距离= 用sas简单过程计算马氏距离很方便，但考虑到常用的随机向量都服从正态分布，sas中discrim过程采用最大概率判别，当各个总体协方差阵相同时，最大概率判别就是马氏距离判别。3.最大概率判别 若x是来自第i个总体，则其概率密度为用样本均值，和样本方差阵代替各总体的均值和方差阵。得到来自第i个总体的近似概率密度 （7.2）若这些值中第j个最大，则判定x属于第j个总体等价的做法是比较的对数的大小。sas discrim过程中比较简化了的，即将舍去共同的，并乘以-2，得到，定义它为x与第i类广义（平方）距离，这样的广义距离称为判别函数。若x与第j类广义距离最小，则判定x属于第j类。当假定各总体协差阵相同时，用联合协差阵 （7.3）代替。其中。这时只需比较的大小，即马氏距离判别。而且，为了简化，可以略去共同的，判别函数件化为线性函数4贝叶斯判别方法 当总体有一定先验概率时，先验概率对判别有较大影响，例如某地历史上多年统计的结果，旱涝年先验概率很小，常年先验概率很大，当气象指标（样品）既接近常年又接近涝年时，应当判为常年，这就是贝叶斯方法。sas的计算原理是：设第i个总体的先验概率为，则由贝叶斯定理，样品x属于第i个母体的后验概率，近似为其中由(7.3)定义。若其中最大，则判定x属于第i个总体。值得注意的是，我们这儿介绍的贝叶斯判别与2、3中贝叶斯判别不同，不涉及损失函数。sas的discrim过程能实施以上介绍的贝叶斯方法。 5投影降维问题 从几何角度考虑，判别就是把分为k个部分：；每一部分对应一个类，例如对应第i类,i=1,，k；当观测值x落在时，就判定x属于第i类。一般判别问题中随机向量维数p很大。很不直观，能否把这p维空间投影到较低维的1、2维空间，使投影后较直观？一般总希望原有的观测值（对应1，2维空间中的点）在新坐标系中尽可能分开。 2、3中介绍的fisher判别分析是一种方法，sas中则使用另一种方法一典型判别，它借助于主分量分析和典型相关分析形成投影方向。sas的candisc过程执行典型判别分析。6模型选择 实际问题中随机向量维数p有时很大。特别在作多项式判别时，可观测变量的多项式项数很多，p也就很大但每一变量都对判别有用，选出合适的随机变量。以保证每一选出的随机变量都对判别作出显著贡献。这就是模型选择。sas包含的stepdisc过程就能有效地实施逐步判别法等方法做模型选择。7.2 用discrim过程实施最大概率判别和贝叶斯判别 sas用discrim过程实施最大概率判别和贝叶斯判别。discrim过程主要有三个语句：proc discrim语句，class语句和var语句。1 proc discrim语句。用以调用discrim过程。其中data选项用以指定训练数据集，即已知类别的若干个样品所成数据集（各个总体的已知样品所成数据集），由于这些样品所属总体已知，可以用来建立判别法则。pool=选项用以说明按各个总体协差阵相等与否的模型计算：当poolyes，指定总体协差阵相同，按(72)式计算协差阵估计值；当poolno时，指定总体协差阵不同，按(71)式计算协差阵估计值。若此选项缺省，按poolyes处理。list选项用以指式sas打印每个样品分类结果。选项out=指示sas把分类结果存入指定的数据库。2class语句用以说明分类变量，该分类变量的值用以决定母体的类别。例如class a；指示sas把a作为分类变量。a取不同值，表示不同的类。3var语句用以规定作判别分析的随机变量,即样品的分量。例如var x1-x3 u1 v2；以x1,x2,x3,u1,v2作为预报因子，它们形成5维随机向量。4.当我们需要作贝叶斯判别时，可以加上priors语句。priors语句有如下形式：(1)priors prop；此语句说明先验概率等于训练数据集中各类样品个数所占比例，即先验概率为。(2)priors 1= 03 2=05 302；或priors a=0.1 b=0.3 c=0.5 d=0.1; 这里用空格作为区分符，前一语句说明分类变量取1，2，3的类的先验概率分别是03，05，02；后一语句说明分类变量取a,b,c,d类的先验概率分别为0.1,0.3,0.5,0.1。当使用priors语句时， discrim过程计算样品属于某一类的后验概率，同时给出判别函数。当poolyes时，判别函数是线性的，当poolno时，判别函数是二次函数。对例7.1使用discrim过程的过程： 用19511975的数据，建立分类规则(判别函数)，以19761979的数据，作为待判数据，将这4年分类为此，我们可以采用缺省技巧，将这4年分类信息（species）取为“.”。与1951-1975年资料一同输入计算机，再用list选项列出29年分类结果，包括原来已分好类的25年数据的回代分类；为此采用sas程序data rainfall;input year x1-x4 species;cards;1951 0.58 82.0 44.0 40.6 11952 0.40 83.0 18.0 43.0 21953 0.55 85.0 36.0 30.7 21954 0.40 85.0 36.0 40.7 21955 0.48 88.0 49.0 43.0 21956 0.41 82.0 35.0 78.6 31957 0.65 80.0 29.0 33.2 11958 0.45 82.0 32.0 33.1 31959 0.39 81.0 27.0 46.5 31960 0.34 85.0 28.0 41.7 31961 0.42 84.0 38.0 20.4 31962 0.52 86.0 38.0 0.2 11963 0.46 88.0 25.0 56.7 21964 0.48 83.0 46.0 13.6 11965 0.53 84.0 41.0 32.3 11966 0.65 81.0 31.0 28.9 11967 0.66 83.0 38.0 46.6 11968 0.53 80.0 42.0 93.1 31969 0.56 85.0 18.0 16.3 31970 0.45 83.0 37.0 23.9 31971 0.34 80.0 42.0 26.3 31972 0.41 79.0 38.0 40.8 31973 0.53 83.0 23.0 61.3 31974 0.48 84.0 19.0 23.2 21975 0.30 85.0 27.0 17.5 21976 0.42 81.0 21.0 52.2 .1977 0.52 81.0 38.0 45.8 .1978 0.36 82.0 34.0 34.9 .1979 0.43 84.0 34.0 60.5 .;proc discrim out=wu list;class species;var x1-x4;id year;run;其中语句proc discrim out=wu list;功能是调用discrim过程，将分类结果存入数据库wu，并打印分类结果。语句class species;指示变量species是分类变量。语句id year;指示将变量year作为id变量。执行此程序得到的输出是 class level information output prior species sas name frequency weight proportion probability 1 _1 7 7.0000 0.280000 0.333333 2 _2 7 7.0000 0.280000 0.333333 3 _3 11 11.0000 0.440000 0.333333以上是1951-1975（25个已分好类）的简单信息，用它们的数据为训练数据，计算协差阵，判别函数。每类先验概率相等都是1/3。 discriminant analysis pooled covariance matrix information covariance natural log of the determinant matrix rank of the covariance matrix4 6.2109918以上说明pool选项缺省，按poolyes处理用(8.2)式计算联合协差阵。其秩为4，行列式的自然对数为6.2109918。 pairwise generalized squared distances between groups 2 _ _ -1 _ _ d (i|j) = (x - x ) cov (x - x ) i j i j generalized squared distance to species from species 1 2 3 1 0 10.55484 7.47484 2 10.55484 0 2.84068 3 7.47484 2.84068 0以上按计算各类间的距离的平方。 discriminant analysis linear discriminant function _ -1 _ -1 _ constant = -.5 x cov x coefficient vector = cov x j j j species 1 2 3 constant -960.94226 -1018 -948.80619 x1 76.64880 39.25276 40.56251 x2 22.37432 23.33004 22.48088 x3 0.30314 0.08712 0.13713 x4 0.54023 0.64151 0.63430以上给出线性判别函数的具体形式判为第1类的函数为判为第2类的函数为判为第3类的函数为 classification results for calibration data: work.rainfall resubstitution results using linear discriminant function generalized squared distance function: 2 _ -1 _ d (x) = (x-x ) cov (x-x ) j j j posterior probability of membership in each species: 2 2 pr(j|x) = exp(-.5 d (x) / sum exp(-.5 d (x) j k k以上说明判别观测值x为第j类的算法：先计算他与第j类的马氏距离，再计算他判为j类的概率。classification results for calibration data: work.rainfallresubstitution results using linear discriminant functionposterior probability of membership in species:year from classified species into species 1 2 31951 1 1 0.9692 0.0027 0.02811952 2 2 0.0002 0.5139 0.48581953 2 1 * 0.6735 0.2085 0.11801954 2 2 0.0030 0.6959 0.30111955 2 2 0.0564 0.8812 0.06241956 3 3 0.0004 0.1978 0.80191957 1 1 0.9885 0.0004 0.01111958 3 3 0.0681 0.1527 0.77921959 3 3 0.0012 0.1135 0.88531960 3 2 * 0.0001 0.7898 0.21021961 3 3 0.1088 0.3828 0.50841962 1 1 0.9120 0.0667 0.02141963 2 2 0.0000 0.9815 0.01851964 1 1 0.9258 0.0119 0.06241965 1 1 0.7907 0.0794 0.12991966 1 1 0.9936 0.0004 0.00601967 1 1 0.9891 0.0024 0.00851968 3 3 0.0339 0.0283 0.93781969 3 2 * 0.2472 0.5982 0.15461970 3 3 0.2431 0.1895 0.56751971 3 3 0.0195 0.0229 0.95761972 3 3 0.0355 0.0121 0.95241973 3 3 0.0123 0.4368 0.55101974 2 2 0.0211 0.6338 0.34511975 2 2 0.0001 0.7775 0.22241976 . 3 * 0.0007 0.1477 0.85161977 . 3 * 0.4842 0.0301 0.48561978 . 3 * 0.0033 0.1676 0.82911979 . 2 * 0.0016 0.5469 0.4514* misclassified observation以上给出历年原来所属类（第2列），判归的类（第3列，加*号的表示误判）以及判归各类的概率。1976年判为第3类，1977年判为第3类，1978年判为第3类，1979年判为第3类，sas认为它们原来属于待判的类，判归第2，3类都是错判，故加*号。 discriminant analysisclassification summary for calibration data: work.rainfallresubstitution summary using linear discriminant function generalized squared distance function: 2 _ -1 _ d (x) = (x-x ) cov (x-x ) j j j posterior probability of membership in each species: 2 2 pr(j|x) = exp(-.5 d (x) / sum exp(-.5 d (x) j k k以上总结判别方法：马氏距离与判为第j类的计算公式。 classification summary for calibration data: work.rainfall resubstitution summary using linear discriminant function number of observations and percent classified into species:from species 1 2 3 total. 0 1 3 4 0.00 25.00 75.00 100.00 1 7 0 0 7 100.00 0.00 0.00 100.00 2 1 6 0 7 14.29 85.71 0.00 100.00 3 0 2 9 11 0.00 18.18 81.82 100.00 total 8 9 12 29percent 27.59 31.03 41.38 100.00priors 0.3333 0.3333 0.3333以上总结判别情况：将原来属于待判，1，2，3类的年份判为1，2，3类的个数及百分比。因为没有给出先验概率，先验概率按相等概率计算。这等价于不考虑先验概率的最大概率原则。若上例改用bayes判别，由25年资料决定先验概率，则可改用程序proc discrim outstat=info method=normal list;class species;var x1-x4;priors prop;id year;run;得到的输出中与非bayes判别不同的是 class level information priorspecies frequency weight proportion probability1 7 7.0000 0.280000 0.280000 2 7 7.0000 0.280000 0.2800003 11 11.0000 0.440000 0.440000给出先验概率0.28，0.28，0.44 pairwise generalized squared distances between groups 2 _ _ -1 _ _ d (i|j) = (x - x ) cov (x - x ) - 2 ln prior i j i j j generalized squared distance to species from species 1 2 3 1 2.54593 13.10077 9.11680 2 13.10077 2.54593 4.48264 3 10.02077 5.38661 1.64196以上按计算各类间的距离的平方。discriminant analysis linear discriminant function_ -1 _ -1 _constant = -.5 x cov x + ln prior coefficient vector = cov x j j j j species 1 2 3constant -962.21523 -1019 -949.62717 x1 76.64880 39.25276 40.56251 x2 22.37432 23.33004 22.48088 x3 0.30314 0.08712 0.13713 x4 0.54023 0.64151 0.63430以上给出线性判别函数的具体形式判为1，2，3类的判别函数分别为对于给定x的观测值，计算上述判别函数，x属于函数值最小的类。 discriminant analysis classification results for calibration data: work.rainfall resubstitution results using linear discriminant function generalized squared distance function: 2 _ -1 _ d (x) = (x-x ) cov (x-x ) - 2 ln prior j j j j posterior probability of membership in each species: 2 2 pr(j|x) = exp(-.5 d (x) / sum exp(-.5 d (x) j k kyear from classified species into species 1 2 31951 1 1 0.9539 0.0026 0.04341952 2 3 * 0.0002 0.4023 0.59751953 2 1 * 0.6310 0.1953 0.17371954 2 2 0.0026 0.5938 0.40371955 2 2 0.0544 0.8508 0.09481956 3 3 0.0003 0.1356 0.86411957 1 1 0.9823 0.0004 0.01741958 3 3 0.0471 0.1056 0.84721959 3 3 0.0008 0.0754 0.92381960 3 2 * 0.0001 0.7051 0.29491961 3 3 0.0843 0.2966 0.61911962 1 1 0.9010 0.0659 0.03321963 2 2 0.0000 0.9712 0.02871964 1 1 0.8939 0.0114 0.09461965 1 1 0.7361 0.0739 0.19001966 1 1 0.9903 0.0004 0.00931967 1 1 0.9844 0.0024 0.01321968 3 3 0.0221 0.0184 0.95951969 3 2 * 0.2271 0.5496 0.22331970 3 3 0.1835 0.1431 0.67341971 3 3 0.0126 0.0148 0.97261972 3 3 0.0230 0.0078 0.96921973 3 3 0.0093 0.3322 0.65851974 2 2 0.0176 0.5294 0.45291975 2 2 0.0001 0.6899 0.31001976 . 3 * 0.0005 0.0993 0.90021977 . 3 * 0.3790 0.0236 0.59741978 . 3 * 0.0022 0.1137 0.88411979 . 3 * 0.0013 0.4348 0.5639 * misclassified observation classification summary for calibration data: work.rai