第二十六章多元分析

第二十六章多元分多元分析（multivariateysis）是多变量的统计分析方法，是数理统计中应用广 将认识对象进行分类是人类认识世界的一种重要方法，比关世界的时间进程 1）d(x,y)0，x,yd(xy0xyd(x,y)d(y,x)，x,yd(xyd(xzd(zyx,yz分析中，对于定量变量，最常用的是Minkowski距离 qdq(x,y)xkk

，qq1,2qpd1(x,y)xkk

2d2(x,y)xkk

Chebyshevd(x,y)maxxk

Minkowski距离中，最常用的是欧氏距离，它的主要优点是当坐标轴进行正交值得注意的是在采用Minkowski距离时，一定要采用相同量纲的变量。如果变量算距离。在采用Minkowski距离时，还应尽可能地避免变量的多重相关性（multicollinearity由于Minkowski距离的这些缺点，一种改进的距离就是马氏距离，定义如下马氏（Mahalanobis）距(xy)T1(x(xy)T1(xxy为来自p维总Z的样本观测Z的协方差矩阵，实际中往往是不如果有两个样本类G1和G2，我们可以用下面的一系列方法度量它们间的距离最短距离法（nearestneighbororsinglelinkageD(G1,G2)min{d(xi,yj y最长距离法（farthestneighbororcompletelinkageD(G1,G2)max{d(xi,yj y重心法（centroidD(G1,G2)d(x,y) x,y分别为G1G2的重心类平均法（groupaverageD(G,G) d(x,x) n 12xiG1x它等于G1G2中两两样本点距离的平均，式中n1,n2分别为G1G2中的样本点个数D(xx)T(xx)，D(xx)T(xx)

xD12

x)T(xkx)x

xi，x2 xj，xn

1 1

1 2x 2D(G1,G2)D12D1 事实上，若G1G2内部点与点距离很小，则它们能很好地各自聚为一类，并且这两类又能够充分分离（D12很大），这时必然有DD12D1D2很大。因此，按定义可以认为，两类G1G2Ward1936年提出，Orloci1976年发展起来的，故又称为Ward方法。（b）1记{w1w2…w7}，聚类结果如下：当距离值为f5时，分为G1w1w2w3w4w5w6w7}；距离值为f4分为两类：G1w1w2w3}G2w4w5w6w7}；距离值为f3分为三类：G1{w1w2w3}G2{w4w5w6}G3{w7}；距离值为f2分为四类：G1{w1,w2,w3}，G2{w4,w5}，G3{w6}，G4{w7距离值为f1分为六类G1{w4,w5}，G2{w1}，G3{w2}，G4{w3}，G5{w6}，G6距离小于f1分为七类，每一个点自成一类怎样才能生成这样的聚类图呢？步骤如下：设{w1,w2,…,首先构造 个类，每一个类中只包含一个样本点，每一类的平台高度均为零法（又称最近邻法），最先由Florek等人1951年和Sneath1957年引入。下面举例说明例1设有5个销售w1w2w3w4w5，他们的销售业绩由二维变量(v1v2描述，见表1。1v1（销售量）v2（回收款项）1011324325记销售员wi(i1,2,3,4,5)的销售业绩为(vi1vi2。如果使用绝对值距离来测量点2d(wi,wj)vikvk

，D(Gp,Gq)min{d(wi,wj由距 w1w2w3w4 5 5

4第一步，所有的元素自成一类H1{w1w2w3w4w5。每一个类的平台高度为f(wi)0(i1,2,3,4,5)D(GpGq)d(wpwq。第二步，取新类的平台高度为1w1w2合成一个新类h6，此时的H2{h6,w3,w4,第三步，取新类的平台高度为2w3w4合成一个新类h7，此时的分类情况H3{h6,h7,第四步，取新类的平台高度为3，把h6h7合成一个新类h8，此时的分类情况H4{h8,第五步，取新类的平台高度为4，把h8w5合成一个新类h9，此时的分类情况H52b有了聚类图，就可以按要求进行分类。可以看出，在这五个中 的工作绩最w3w4的工作成绩较好，而w1w2的工作成绩fori=1:mfornd=nonzeros(d);%去掉d中的零元素，非零元素按列排列nd=union(nd,nd)去掉重复的非零元素fori=1:m-[row,col]=find(d==nd_min);tm=union(row,col);%row和col归为一类tm=reshape(tm,1,length(tm));%把数组tm变成行向量iflength(nd)==0y=pdist(a,'cityblock')求a的两两行向量间的绝对值距离yc=squareform(y)%变换成距离方阵z=linkage(y)%产生等级聚类树[h,t]=dendrogram(z)%画聚类图T=cluster(z,'maxclust',3)%把对象划分成3类fori=1:3tm=find(T==i);%求第i类的对象Y=pdist(X)计算mn矩阵X（看作mn维行向量）中两两对象间的欧氏距离。对于有m个对象组成的数据集，共有(m1)m/2个两两对象组合。输出Y是包含距离信息的长度为(m1m2的向量。可用squareform函数将此向表2’metric’含欧氏距离（缺省马氏距离（Mahalanobis距离闵氏距离（Minkowski距离输出的(m1)m/2维距离行向量。表3’method’含最短距离（缺省离差平方和方法（Ward方法输出Z为包含聚类树信息的(m1)3矩阵。聚类树上的叶节点为原始数据集中的对象，由1到m。它们是单元素的类，级别更高的类都由它们生成。对应于Zj每个新生成的类，其索引为mjm为初始叶节点的数量。表4cutoff 进行了量化。如果接的不一致系数大于阈值，则clusterxxijxssjxjsj是矩X(xij)mn每一列的均值和标准差和Y中的距离信息（由pdist()函数产生）进行比较。Z为(m1)3含在第三列。Y是(m1)m/2维的行向量例如，给定距离为Y的一组对象{1,2,…,m}，函数linkage()生成聚类树。((yijy)(ziji(yijy)2(zijiicyij为Y中对象ij间的距离；zij为Z(:,3)中对象ij间的距离；yz分别为Yx的取值(xx,…x)TRnj1,2,…m。则可以用两变x 1 2 n(xijxj)(xikxkrjk

12

))

(xijxj (xikxk) rjk

1

2 2xxxx rjk1jkrjkrkjjk

R(G1,G2)max{djk x 1 或d21r2，这时，R(G,G)与两类中相似性最小的两变量间 R(G1,G2)min{djk xxk 1 或d21r2，这时，R(G,G)与两类中相似性最大的两个变量 例2服装标准制定中的变量聚类法。量资料，获得各因间的相关系数表（见表5）。表5成年女子各部位相关系数 111111111111 x1上体长，x2长，x3，x4颈围，x5总肩围，x6总胸宽，x7后背宽，x8前腰节高，x9后腰节高，x10总体长，x11身高，x12长，x13腰围，x14臀围。用最大系数法对这14个变量进行系统聚类，分类结果如图3。图3成年女子14d=tril(d);%提出d矩阵的下三角部分b=b';%化成行向量 ind1=find(y==1);ind1=ind1'ind2=find(y==2);ind2=ind2'dendrogram(z);%画聚类图体长，长，前腰节高，后腰节高，总体长，身高，长；另一类是反映胖瘦§2指标的原始数据取自《中计年鉴，1995》和《中国教育统计年鉴，1995》除以各地区相应的人口数得到十项指标值见表6其中：x1为每百万人口高等院校数；x2为每十万人口高等院校毕业生数；x3为每十万人口高等院校招生数；x4为每十万人口高等院校在校生数x5为每十万人口高等院校教职工数x6为每十万人口高等院校专职教师数；x7为高级占专职教师的比例；x8为平均每所高等院校的在校生数；x9为国家财政预算内普通高教经费占国内生产总值的；x10为生均教育经费。图4高等教育的十项评价指标表6内种想法，运用软件计算十个指标之间的相关系数，相关系数矩阵如表6所示。表6

5load %计算相关系数矩 %取出矩阵d的下三角元素 fori=1:6 tm=reshape(tm,1,length(tm));%变成fprintf('第%d类的有%s\n',i,int2str(tm));%显示分类

x1：每百万人口高等院校数x7：高级占专职教师的比例； ：国家财政预算内普通高教经费占国内生产总值的x10：生均教育经费

6loadgj.txt%把原始数据保存在纯文本文件gj.txtgj=zscore(gj);%数据标准化z=linkage(y,'average');%按类平均法聚类dendrogram(z);%画聚类图forT=cluster(z,'maxclust',k);%把样本点划分成k类fori=1:ktm=find(T==i);%求第i类的对象if第一类：；第二类：；第三类：其他地区第一类：；第二类：；第三类：，；第四类：其他地区。国的政治、经济与文化中心的地位是吻合的。和作为另外两个较早的直辖市，源相对匮乏。作为一个非常特殊的民族地区，其高等教育状况具有和其它地区不同施对、、青海和地区进行扶持，促进当地高等教育事业的发展。§3主成分分析主成分分析（principalcomponentysis）是1901年Pearson对非随量引如果x1x2·xp表示p门课程c1c2·cp表示各门课程的权重，那么sc1x1c2x2·cp 绩，记为s1s2,·,snn为学生人数。如果这些值很分散，表明区分得好，即是说，需要寻找这样的，能使s1s2·,sn尽可能的分散，下面来看它的统计定义。设X1,X2,·,Xp表示以x1,x2,·,xp为样本观测值的随量，如果能找c1c2·cp，使Var(c1X1c2X2·cpXp p个c2c2·c2 Zi表示第ii1,2,·pZ1c11X1c12X2·c1pXZ2

X1c22X2·c2pXZcX X·c p1

p

其中对每一个i，均有c2c2·c21，且(cc,·c使得Var(Z i 11 1 到最大；(c21c22·c2p)不仅垂直于(c11c12·c1p)，而且使Var(Z2的值达到最大；(c31c32·c3p)同时垂直于(c11c12·c1p)和(c21c22·c2p)，并使Var(Z3的值达到最大；以此类推可得全部p个主成分，这项工作用手做是很繁琐的，但借助于计主成分估计（principalcomponentestimate）是Massy在1965年，它是回时表现出的不稳定性而。设有p个回归（自）变量x1,x2·,xp，它在第i次试验中的取xi1xi2·xip（i1,2,·n）

·x1p·x2pX xx

xn

xnpY01X~N(0,2I)， 其中Y(Y,Y,·,Y)Tn1为未知参数，1为所有元素均为1的n维列向 均已标准化，如X未标准化，需要作变量的标准化变换(xjxj)/sj，其中xj,sj分别为X的第j列的均值和标准差），此时ˆY1n0ni10

zcxc

…cx，c2 p11 22 pp z(i)c1xi1c2xi2…cpxip（i1,2,…,n） znz(i)ncjxijncjxij wcc,…c)T

i1 M*2

z)2

nz21(Xw)T(

(i 2M新变量可以去掉。反之，2M2zi的取值与ci的选取有关。因此，我们总是希望所选择的ci(i1,2,…,p2M*2XTX的特征值12p,,…,M*Xw)TXwnw时达到，且最大值为n zxxx,…x)TzxT zixTi（i1,2,…,k

0，i1,2,…,k，wTw12在条件1）之M*达到最大2k1zk1xTk1zixTi（i1,2,…p）现在回到线性模型（19）x1x2…xp变换为主z1z2…zp之后再z11 …z2ZX(1,2,…,p) zz

zn

znp记Q1,2,…,pppQZXQ，引入新参数QT，或者Q，则Y01ZQT01Z 0 0ZTZQTXTXQQT(XTX)Q

p剔除变量zi。如果r1p0zr1zr2.zpr个分量1,2.,r，设它的最小二乘估计为ˆ1,ˆ2,.,ˆr，而后面的pr个分量则以0作为它们的估计，然后由关系式Q即可确定的估计，我们称的主成分估计，先将Q,1Q(Q1,Q2)， 2 其中Q1pr 1 Q (Q1,Q2)0 1 理论上表明：主成分估计在设计阵时优于LS估计，但（31）在特征值为1的附定义1若存在1rp，使r1r1Adiag(11,.,r1 ,., r 这里p,1)为平稳参数，我们称ˆQAQTQˆ的单参数主成分1例3Hald水泥问题，含如下四种化学成 x23CaO·SiO2的含量x34CaO·Al2O3·Fe2O3的含量（％），x42CaO·SiO2的含量组，见表7，对数据实施标准化，则XTX/12就是样本相关系数阵（见表8）。表7Haldy1762138485766973681924198表8Hald1--1----1--11T12T23T3yˆ62.40541.5511x10.5102x20.102x3yˆ85.74331.3119x10.2694x20.1428x3loadsn.txt%把原始的x1,x2,x3,x4,y的数据保存在纯文本文件sn.txtr=corrcoef(x0)%计算相关系数矩阵xb=zscore(x0);yb=zscore(y0);%对y0进行标准化处理 contr=cumsum(t)/sum(t)%计算累积贡献率，第i个分量表示前i个主成分的贡献率 fori=1:n-hg=s(:,1:num)\yb;%主成分变量的回归系数hg=c(:,1:num)*hg;%标准化变量的回归方程系数fori=1:n-1何筛选这些特征值？一个实用的方法是删去r1,r2,…,p后，这些删去的特征值之和占整个特征值之和i的15％以下，换句话说，余下的特征值所占的 论很多，有很多比较成方法，这里不一一介绍。变量xi的贡献值为r r(zjxizjxi的相关系数。例4设x(x,x,x XTX 0 则可算得15.8284，20.1716，如果我们仅取第一个主成分，由于其累积贡r2(z,x) 可见，第一个主成分对第三个变量的贡献值为0，这是因为x3和x1,x2都不相关。由于x3的信息，这时只选择一个主成分就不够了，需要再例5研究纽约市场上五种的周回升率。这里，周回升率＝（本星期五场收盘价－上星期五市场收盘价）/上星期五市场收盘价。从年月到年月，因此，不同周末回升率是彼此相关的。x1x2…x5分别为五只xT(0.0054,0.0048,0.0057,0.0063,RR12.857，20.809，30.540，40.452，50.3431T(0.464,0.457,0.470,0.421,12T(0.240,0.509,0.260,0.526,2x0.457x0.470x0.421x0.421 20.240x0.509x0.260x0.526x 512100%5系数为正的三只都是化学工业上市企业）和石油（在z2中系数为负的两只股§4假设进行主成分分析的指标变量有mx1x2.xmn个评价对象，第ijaa~， ~ j，（i1,2,.,n；j1,2,.,msjs21 2jnaijsjn1(aijj)j1,2,.mjsjj xxii，（i,,.,mii~nki~nrij ，（i,j1,2,.,mnR的特征值12m0，及对应的特征向量u1u2.,um，其中uju1ju2j,.umj)T，由特征向量组成my1

~u 11

21x2.xy2

~u 12

22x2.um2

u1m

2.umm~1my1是第1主成分，y2是第2主成分，…，ymm主成分~1mp（pm）①计算特征值jj1,2,.m的信息贡献率和累积贡献率。kbjk

（j1,2,·,mkppmpkmk为主成分y1y2·yp的累积贡献率，当p接近于1（p0.850.900.95）时，则选择p个指标变量y1y2·,yp作为p个主成分，代替原来m个指标变量，从而可对p个主成分进行综合分析。ppZbjyj其中bjj表7主成分分析结果123456表8标准化变量的前41234567890------------------y0.3497x0.359x·0.2452x 12y0.1972~x0.0343x·0.28612

y

0.1084x·0.8637 y0.1022~x0.2266x·0.2457 Z0.7502y10.1577y20.0536y30.0206表9和综合评价结123456789--0--------内--0--------loadgj.txt%把原始数据保存在纯文本文件gj.txt中gj=zscore(gj);%数据标准化r=corrcoef(gj)[x,y,z]=pcacov(r)%y为r的特征值，z为各个主成分的贡献率x=x.*f;%修改特征向量的正负号，每个特征向量乘以所有分量和的符号函数值num=4;%num为选取的主成分的个数df=gj*x(:,1:num);%计算各个主成分的得分tf=df*z(1:num)/100;%计算综合得分[stf,ind]=sort(tf,'descend');%把得分按照从高到低的次序排列stf=stf',ind=ind'方面。陕西和东北三省高等教育发展水平也比较高。、广西、、等地区高等教育发展水平比较，这些地区的高等教育发展需要政策和的扶持。值得一提的是、、等经济不发达地区的高等教育发展水平居于中上游水平，可能是§5因子分析因子分析（factorysis）是由英国心理学家Spearman在1904年提出来的，他部分。初学因子分析的最大在于理解它的模型，我们先看如下几个例子。因子来确定，即可设想第iXiF1F2…F5的Xiiai1F1ai2F2…ai5F5Ui（i1,2,…,N）i是总平均，Ui是第i个学生的能力和知识不能被这五个因子包含的部分，称为特殊因子，常假定Ui~N(0,2，不难发现，这个模型与回归模型i在形式上是很相似的，但这里F1,F2…,F5的值却是未知的，有关参数的意义也有很大 例7诊断时，医生检测了的五个生理指标：收缩压、舒张压、心跳间隔、呼例8Holjinger和Swineford在芝加哥郊区对145名七、八年级学生进行了24个心理推理因子和因子。到新变量zi，则可以建立因子分析模型如下：ziai1F1ai2F2…aimFmUi（i1,2,…,p）， Fj(j1,2,…m)出现在每个变量的表达式中，称为公共因子，它们的含义要根2,…,p，j1,2,…m）称为因子载荷，A(aij称为载荷矩阵。zAFz(z,z,…,z)T，F(F,F,…,F)T UU1,U2,…,Up)TAaijpm。Cov(U)diag(2,2,…,2 Cov(F,U)

各公共因子都是均值为0，方差为1的独立正态随量，其协方差矩阵为单位ImF~N(0,Im)。当因子F的各个分量相关时Cov(F)不再是对角阵，这样im个公共因子对第i个变量方差的贡献称为第i共同度，记为h2ih2a2a2… i i而特殊因子的方差称为特殊方差或者特殊值（即（37）式中的2i1,2,…pi从而第iVarzh22，i1,2,…, 设12p为样本相关系数R的特征值，1,2,…,p为相应的标准A( 22 mm) 21 因子分析主成分解见表10，对m2RAATCov(U)为 0.0173

0

0

0.0118 表10因子分析主成分解1-2-3-45这表明AAT产生的数比R中对应元素（相关系数）要大。第一个因子F1代表了一般经济条件，称为市场因子，所有在这个因子上的载0.4620.3220.4260.523 a2=a(:,[1,2])%提出两个因子的载荷矩阵tcha2=diag(r-a2*a2')%计算两个因子的特殊方差ccha2=r-a2*a2'-diag(tcha2)%求两个因子时的残差矩阵 上面主成分解是不唯一的，因为对A作任何正交变换都不会改变原来的AAT，即设Q为m阶正交矩阵，BAQBBTAAT，载荷矩阵的这种不唯一性表面看是B中有尽可能多的元A(aij)，i1,2,·,p，j

sinQ

sin cosBAQ(bij)，i1,2,·,p，jzB(QTF)

的两个部分，因此，要求(b2b2,·b2)和(b2b2,·b2)T这两列数据分别求得的 p1pb2

1pb2V V h2 h2i1i i1i 的影响，正交旋转的目的是为了使总方差VV1V2d0应满

tan4 D02A0B0/C(A2B2)/

A0ui B0 C0(uivi D0 a

2auii1i2 h h

i1i i i 当m2时，还可以通过图解法，凭 公共因子数m2可以每次考虑不同的两个因子的旋转，从m个因子中每次选两个旋转，共有只会变大k次循环后的V(k)与上一次循环的V(k1)比较变化不大时停止旋转。例10 1/ 2/R1/ 02/ 1解1）R由特征方程det(RI)0可得三个特征值，依大小次序记为11.74542130.2546，由于前面两个特征值的累积方差贡献率已达91.51，因而只1T(0.7071,0.3162,12T(0,0.8944,2 A 0.4472 Q 0.3625，A 0.9320

0.1139 r=[1-1/32/3;-1/310;2/30[vec,val,con]=pcacov(r)%val为r的特征值，con为各个主成分的贡献率 num=2;%选择两个主因子[b,t]=rotatefactors(a(:,1:num),'method'varimax

R

1

1

1表11因子分析表QT1QT212-34-5--AATCov(UR比较接近，所以从直观上，我们可以认为两个因子的模很明显，变量2，4，5在QTF上有大载荷，而在QTF上的载荷较小或 反，变量1，3在QTF2上有大载荷，而在QTF1上的载荷却是可以忽略。因此，我们有理由称QTF为营养因QTF为滋味因子。旋转的效果一目了然。 loadr.txtr.txt中 num=2;%num为因子的个数 ccha=r-a1*a1'-diag(tcha)%求残差矩阵 是对不可观测的随机向量Fi取值的估计。通常可以用最小二乘法和回归法来估计§6因子分析案例因子分析(factorysis)是一种数据简化的技术。它通过研究众多变量之间的Xiiai1F1…aimFmi（mp 或X1 1 ·a1mF1 1X · F 22222m: : :::: :aX aX

·

p或

p p pmm pXAFX1 1 ·a1m 1X · X 2,2

，A

2m

，

2: : : : Xp

p

ap

·apm

pF1,F2·,Fp为公共因子，是不可观测的变量，它们的系数称为载荷因子。i是特殊因子，是不能被前m个公共因子包含的部分。并且满足E(F)0，E()0，Cov(F)ImD(Cov()diag(2,2,·,2cov(F)0 XAF，得CovXACov(FATCov(Cov(X)AATdiag(2,2,·,2 2,2,·,2的值越小，则公共因子共享的成分 p p设为一

AT

T

X~~aijX的共同度是因子载荷矩阵的第ih2

a2iVar(X)a2Var(F)·a2Var(F)Var(

即m1 a2

jiipSjFjj1,2,·m对所有的XiFjpSja2 RAATR*AATRiR*为约相关系数矩R*对角线上的元素是h2，而不是1i ·rr

·r1pR*RD 2p ; ·hˆ2 p p直接R*的前p个特征值和对应的正交特征向量。得到如下的2*pA[* * u*2*p1 R***·*，对应的正交特征向量为u*u*,·u* ii取hˆ2R2R2xxx 的p1xj的线性组合联系起来的。取hˆ2max

取hi

ppi取hˆ21/rii，其中riiR1的对角元素i 1/ 1/5 1/ 2/ 2/ 解特征值为11.546420.853630.60.4597 0.8881 u0.628，u0.3251，u

A 2 3u3

x10.5717F1

x20.7809F10.3003F2x30.7809F10.3003F2r=[11/5-1/5;1/51-2/5;-1/5-2/5[vec,val,con]=pcacov(r)%val为r的特征值，con为各个主成分的贡献率num=input('请选择公共因子的个数：');%交互式选取主因子的个数 aa=a(:,1:num);%提出两个主因子的载荷矩阵 s2=sum(aa.^2,2)%计算共同度11/1/51/12/ 2/ 解假定用hˆ2maxrh2h21h22h22i1/R*1/

1/2/

1/52/

1 1 2/ 2/5特征值为10.912320.087730uu

0.9294 A

0.2752 0.0773旋转。目的是使因子载荷阵的结构简化，使载荷矩阵每列或行的元素平方值向0和1两级分化。有三种主要的正交旋转法，方差最、四次方最和等量最。方差最从简化因子载荷矩阵的每一列出发，使和每个因子有关的载荷的平方的方差最大。当数几个变量在某个因子上有较高的载荷时，对因子的解释最简单。方等量最把四次方最和方差最结合起来，求它们的平均最大6.21.X1 .a1mF1 1X . F 2 2m22: :: :X Xa

.

p p pmm pFjj1X1.jpXp，j1,2,.,（1）巴特莱特因子得分(最小二乘法 .a1m . 2m : ap

.apm 11 12 1m x 11 12 1m xaf f. fi

21 22 2m xippap1f1ap2f2.apmfm [(xi)(aˆ ˆ. ˆ)]/ i

fˆ.fˆ xAF

(xAF)TD1(xAF 1 1D

2 pATD1FATD1A(x

Fˆ(ATD1A)1ATD1(x（2）回归方X1 ·a1mF1 1X · F 2 2m22: :: :X Xa

·

p p pmm piaijXFE(XiFj)E[Xi(bj1X1·bjpXpibj1b·b

·

j2j1 jp i ip:b bjp ·1pbj1 a1j · b a 2pj22ja :: :a

bp

·bppbjp

pj ·1p bj1 a1j · b a 2p，j2，2j : : : bp

bpp

bjpb

apja ·bm1 · m2 R : b1 b2 ·bmp(ˆ XR1ij Fˆ为第ijFXnm 表12上市公司数据12347雅px1x2.xpn个评价对象，第i个jxxx， x xj，（i,,.,n；j1,2,.,p）sjs21 2xjnxijsjn1(xijxj)j1,2,.,pxjsjj xxixi，（i1,2,.,pii相关R(rijpnxkinrij ，（i,j1,2,.,pnR的特征值12p0，及对应的特征向量u1u2.up，其中uju1ju2j,.unj)T，初等载荷矩A[ 2 pup选择m（mp）~xbF. 11 1mxpbp1F1.bpm表1312表主因子主因子-Fˆb

·b

，j1,2,·, j1

jp记第ijFjˆb

b

·b

（i1,2,·,n，j1,2,·,m j1

j2i

jp ·bm1 · m2 R ; b1且

b2 ·bmp Xij 0.531x0.1615x0.1831~0.5015x 0.5151x0.581x F44.49F1计算出家上市公司赢利能力的综合得分见表。 45678F-F---F雅F-------F----F-------我们通过相关分析，得出赢利能力F与资产负债率x之间的相关系数为-0.6987，这F0.829loaddata.txt x=zscore(x);%数据标准化r=cov(x)%求标准化数据的协方差阵，即求相关系数矩阵[vec,val,con]=pcacov(r)%进行主成分分析的相关计算num=input('请选择主因子的个数：');%交互式选择主因子的个数vec=vec.*f1;%特征向量正负号转换a=vec.*f2%求初等载荷矩阵am=a(:,1:num);%提出num个主因子的载荷矩阵[b,t]=rotatefactors(am,'method','varimax')%旋转变换,b为旋转后的载荷阵 weight=rate/sum(rate)%计算得分的权重 disy=[score(ind,:)';STscore';ind']%显示排序结果 d,stats%显示回归系数，和相关统计量的值6.5率的影响因素分对率进行分析。均国民收入。表是年中国个省、、直辖市的数据。16率有关数表17特征根与各因子的贡献1表18因子分析表QT1QT2因子因子1--2----3-4-5-loadsy.txt sy=zscore(sy);%数据标准化r=cov(sy);%求标准化数据的协方差阵，即求相关系数矩阵[vec,val,con]=pcacov(r)%进行主成分分析的相关计算num=input('请选择主因子的个数：');%交互式选择主因子的个数 am=a(:,1:num);%提出num个主因子的载荷矩阵[b,t]=rotatefactors(am,'method','varimax')%旋转变换,b为旋转后的载荷阵 6.6m，单个主成分与综合主成分的分析评价、单因子19Fia1ix1a2ix2·apiaT i1,2,·,ixjbj1F1bj2F2·bjmFmj1,2,·,R为相关系数矩阵，ip(1a1,…,mam为初等因子载荷矩阵（iai同左）C为正交旋转矩ATAI（A为正交矩阵BTBI（B为非正交阵12,…m互不相同aij唯i0,i协方cov(FiFjiijij1,i0,i1,ii（特征值）为主成分Fi的方pvb2Fx k主成Fjx确定Fi是不可观测主成分函数(FF,…F)TAT 因子得分函数(FF,…F) pFx的系数平方和a2 kmb22h221h2称为共同 度2称为特殊jmF(ip)FimpmF(vi/p)Fimp 判别分析（discriminantysis）是根据所研究的的观测指标来推断该所来确定该地区属于哪一种经济类型地区等等。该方法于1921年Pearson的种族相似系数法，1936年Fisher提出线性判别函数，并形成把一个样本归类到两个总体之一判别问题用统计的语言来表达，就是q个总体X1,X2…,Xq，它们的分布函数分别为F1xF2x),Fq(x，每Fi(xp维函数。对于给定的样本X，要判判别和Fisher判别。Mahalanobis通常我们定义的距离是Euclid距离（简称欧氏距离67p1X~N(0,1)Y~N(4,22)7上题，情况并非如此。Ax1.66，也就是说A10是1.661，A24却只有1.172A2更近一点。2xy，协方差为AA内x与y的Mahalanobis距离（简称马氏距离）定义为d(x,y)d(x,A)设总AB的均值向量分别为12，协方差阵分别为1和2，今给一个样本x，要判断x来自哪一个总体。12，12d(xB，然后进行比较d(xAd(xBxA；否则判xB。x

d(x,A)d(x,d(x,A)d(x,现在引进判别函数的表达式 d2(x,A)与d2(x,B)之间的关系，d2(x,B)d2(x,A)(x2)T1(x2)(x1)T1(x 2(x)T1( 122令w(x)(x)T1( x

w(x)w(x)样本的均值与协方差来代替x(1x(1·x(1)是来自An个样本点， x(2x(2·x(2)是来自总体B的n个样本点，则样本的均值与协方差 ˆx(i)i

1nij

x(i)j

j1,

2 (i (i) n (x )n (S1S2

j Si

j

(x(i)x(i))(x(i)x(i))T，i1, ˆ(x)(xx)Tˆ1(x(1)x(2))xx(1)x2x

ˆ(ˆ(12，1w(x)(x)T1(x)(x)T1(x 样本的均值与协方差来代替。因此，对于待测样本x，判别函数定义为ˆ2

(i (i) j n1(xj )(xj )n1Sii1, j Fisher判Fisher判别的基本思想是投影，即将表面上不易分类的数据通过投影到某个方向pX1X2，且都有二阶矩存在。Fisher均值向量分12（均为p维且有公共的协方差矩阵（0。那么线性组合yaTx的均值为 E(y|xX)aT E(y|xX)aT yy

2Var(y)aT( [aT( (aT y aT aTy定理 x为p维 量设yaTx当选取ac1，c0为常数时（54特别当c1yaTx()T 称为Fisher线性判别函数。K1 )1(aTaT)1()T1( 定理 利用上面的记号，取aT(12)T1，则yK0，y

Ky yxX 当x使得()T1x1 12)xxX2

2 W(x)()T1xK(x1())T1(

x 当x使得W(x xX 当x使得W(x 当总体的参数未知时，我们用样本对1,2及进行估计，注意到这里的Fisher判BayesBayes判别和Bayes估计的思想方法是一样的，即假定对研究的对象已经有一设有两个总体X1X2，根据某一个判别规则，将实际上为X1的判为X2或者将实际上为X2的判为X1的概率就是误判概率，一个好的判别规则应该使误判概X1的误判到X2的损失比X2的误判到X1严重得多，则人们在作前一种判断分别具有密度函数f1(x)与f2(x)，其x为p维向量。记x的所有可能观测值的全体，称它为样本空间，R1为根据我们的规则要判为X1的那些x的全体，而X1，但被X2的概率为P(2|1)P(xR2|X1).f1X2，但被判为X1的概率P(1|2)P(xR1|X2).f2类似地，来自X1被判X1的概率，来自X2被判X2的概率分别P(1|1)P(xR1|X1).P(2|2)P(xR2|X2).f2P(正确地判为X1P(来自X1被判为X1P(xR1|X1PX1P(1|1p1P(误判到X1P(来自X2被判为X1P(xR1|X2PX2P(1|2)p2P(正确地判为X2P(2|2)p2P(误判到X2)P(2|1)p1起的损失，并规定L(1|1)L(2|2)0。将上述的误判概率与误判损失结合起来，定义平均误判损失（expectedcostofECM(R1,R2)L(2|1)P(2|1)p1L(1|2)P(1|2)p2 两总体的Bayes定理 极小平均损失（56）的区域R1和R2Rx:f1(x)L(1|2)p2 f(x) L(2|1)p 1Rx:f1(x)L(1|2)p2 f(x) L(2|1)p 1（f1(x)L(1|2)p2xRRf2 L(2 便就将它归入R1xX 当x使得f1x)L(1|2)

L(2

xX 当x使得f1x)L(1|2) f2 L(2 1）x0x01x02,·x0p)Tf1(x0)/f2(x02）损失L(1|2L(2|1p2p11）p2p1xX 当x使得f1(x)L(1|

L(2|

xX 当x使得f1(x)L(1| L(2|2）当L(1|2)L(2|1)1xX 当x使得f(x)1

xX 当x使得f1(x)2 2

3）p1p2L(1|2)L(2|11xX 当x使得f1x)

xX 当x使得f1(x)2 2

（58样如误判损失或者其比值都是难以确定，此时就利用规则（59），如果上述两者都难以f1xf2xxX1xX2。BayesXi~Np(ii（i1,21）12（0Xif(x)(2)p/2||1/2exp{1(x)T1(x 定理4 设总体Xi~Np(i,)（i1,2，其中0，则使平均误判损失极小R1{x:W(x)R{x:W(x) W(x)[x1()]T1( lnL(1|2)L(2|1)

（如果总体的1,2和未知，用式（52）和（53）算出总体样ˆ1ˆ2和ˆ，来代12和，得到的判别函数1W(x)[x

1称为Anderson线性判别函数，判别的规则x 当x使得W(x xX 当x使得W(x) 这里应该，总体参数用其估计来代替，所得到的规则，仅仅只是最优（在平2）12（10，20f1x)/f2x)ln(f1x)/f2xi1对于12时，由定理3可得判别区域：R1{x:W(x)R{x:W(x) W(x)1xT(11)x(T1T L(1|2)p2 T TKlnL(2|1)p2 1

(111222 2显然，判别函数W(x)x的二次函数，它比12ii(i1,2)未知，仍可采用其估计来代替1420是某气象站预报有无春旱的实际资料，x1x2（气象含义从略68个年份的资料，它们的先验概率分别用6/14和8/14来估并设误判损失相建立Anderson线性判别函20序号12345678 - - - - - -W(x1,x2 - --------W(x1,x2-------- 0.3109

， 将上述计算结果代入Anderson线性判别函数W(x)W(x1,x2)2.0893x13.3165x2栏目中，错判的只有一个，即春旱中的第4号，与历史资料的拟合率达93％。26.6 25.5-2.0-- b=[22.1 22.0 22.7 -0.7-1.4 -0.8-1.6 -1.5-1.0-1.2-1.3]';sig1=cov(a);sig2=cov(b);%计算两个总体样本的协方差矩阵symsx1x2x=[x1x2];wx=(x-0.5*(mu1+mu2))*inv(sig)*(mu1-mu2)';%构造判别函数wx=vpa(wx,6)%显示判别函数ahat=subs(wx,{x1,x2},{a(:,1),a(:,2)})'%计算总体1bhat=subs(wx,{x1,x2},{b(:,1),b(:,2)})'%计算总体2样本的判别函数值sol1=(ahat>beta),sol2=(bhat<beta) 12的 26.6-2.0- - b=[22.1 22.0 22.7 -0.7-1.4 -0.8-1.6 -1.5-1.0-1.2-1.3]';cov1=cov(a);cov2=cov(b);%计算两个总体样本的协方差矩阵0.5*(mu1*inv(cov1)*mu1'-mu2*inv(cov2)*mu2Ksymsx1x2x=[x1 ahat=subs(wx,{x1,x2},{a(:,1),a(:,2)})'1bhat=subs(wx,{x1,x2},{b(:,1),b(:,2)})'2sol1=(ahat>=k),sol2=(bhat<k)分类正确率为100％。 26.6 25.5-2.0- - b=[22.1 22.0 22.7 -0.7-1.4 -0.8-1.6 -1.5-1.0-1.2-1.3]';train=[a;b];%train为已知样本prior=[p1;p2];%已知样本的先验概率%线性分类品的式样，包装和耐久性进行了评估后，得分资料见表21。21生产厂家的数据1234567896234247685376463525111111122222（6，4，5（8，1，3（2，4，5788987436242 853764635sample=[645;813;24group=[ones(7,1);2*ones(5,1)];%已知样本的分类[x2,y2]=classify(sample,train,group,'linear')%线性分类[x3,y3]=classify(sample,train,group,'quadratic')%二次分类 典型相关分析（Canonicalcorrelation(x1,x2,·,xp)，(y1,y2,·,yqu1a11x1a21x2·v1b11v

b21

·bq1vbyby…by 12 22 q2u2与u1、v2与v1不相关，但u2和v2相关。如此继续下去，直至进行到r步，两组变量的相关性被提取完为止，可以得到r组变量，这里rmin(pq)。Holing将简单相关系数推广到多个随量与多个随量之间的相关关系的讨论 复相关系数描述两组随量Xx1x2…xp与Yy1y2…yp之间的相 uaTXax,vbTYb i

juv的相关系数。由于uv与投影向量abruvab有关，rr(a,b)我们取在aT a1和bTb1的条件下使r达到最大的a,b作 ruvaTXXbT

X Var(X Cov(X,Y)

XYCovYCov(Y,X Var(Y) Cov(aTX,bTY

YYT D(bTY

aXY aTXXa1和bTYYb1aTXYb的极大值。S(a,b)aT b

(aT

a1)

YYb 的极大值，其中,是Lagrange乘数。

b

a

YXaYYb将上二式分别左乘aT与bTaT

baT a

bTYXabTYYb注意 T，所

aTXY XYbXXaYXaYYb以1左乘（78）第二式得b1 YYb 1 YY

(

YY

(1

)b YX记M1

1

11 M

XYYY

YYYXXXM1a2a，M2b 说明2

M2的特征根非负，均在0和1之间，非零特征根的个数等于min(p,q)，不妨设q。Ma2a的特征根排序为22·2pq个特征根为0 称,·Ma2aaa,·a Mb2b解出的特征向量为bb,·bq uaTX，vbTY，i1,2,·, 还可以证明，当ij时 ,u)Cov(aTX,aTX)aT a XX

Cov(v,v)Cov(bTY,bTY)bTb YY表示一切典型变量都是不相关的，并且其方差为Cov(ui,uj)E(uiuj)Cov(vi,vj)E(vivj)

ij

ii

X与Y的同一对典型变量uivii，不同对的典型变量uiv（ij）之间不相关，也就是说协方差为0iCov(ui,vj)E(uivj)

ii

当总体的均值向量和协差阵未知时，无法求总体的典型相关系数和典型变量，X(1…,X(n和Y(1…,Y(n)为来自总体容量n的样本，这时协方差阵的

n1

(X(i

X)(

(i

X

n1

(Y(i

Y

Y

ˆT n(

X

Y

n

(i

(i

1 1 ˆXn

X(i),Y

Y(i)，用代替并按（81）和（82）求出iab 系数矩阵R取代协方差阵，计算过程是一样的。是描述一个随量y与多个随量（一组随量）X(x,x,…,x)T之间 量ucTX

c i再研究y与u的相关系数。由于u与投影向量c有关，所以ryu与c有关，ryuryu(c)。我们取在cT c1的条件下使r达到最大的c作为投影向量得到的相关系数为偏相

ryucTXX

ryu R 22X ·a1r ·aA ·a 2r r

ap

·apr ·b1r ·bB ·b 2r r p

p

bqrcov(xi,uj)cov(xi,akjxk)akjcov(xi,xkk kxi与uj

(xi,uj)akjcov(xi,xk)k

D(xiq(xi,vj)bkjcov(xi,yk)k

D(xip(yi,uj)akjcov(yi,xk)k

D(yiq(yi,vj)bkjcov(yi,yk)k

D(yiX组原始变量被uipm2(u,

)/p kX组原始变量被vipm2(v,

)/ kY组原始变量被uiqn2(u,

)/ kY组原始变量被viqn2(v,

)/ k ·x1 ·y1q

x2

y2p ·

· nqxyyx11 ·x1p y11 ·y1qyqxyy

x21

x2

yqZ

xn1 ·xnp yn1 ·y2qyq

11

ZTZ

n

YY整体检验（H0XY0H1XY0H0:12·r0（rmin(p,q) Iˆ1

p1

YYYX (1i若1小，则支持H1在原假设为1下，检验的Qn1 (pq1)ln pq2分布。在给定的显著水平下，如果Q2pq H023r0H1:23,…r至少有一非H0kk1r0H1kk1,…r至少有一非零k

(12)，Q[nki

(pq1)]ln2近似服从自由度为pk1)(qk1)2分布。在给定的显著水平下，如果Q2pk1)(qk1，则原假设，认为至少第k对典型变量之间的相关22Xx1用户反馈 x2任务重要性 x3任务多样性，x4任务特殊x5Yy5y6y723y表 X组的典型变u-----------25u---------y--------26-----------u---------y-------表2712345可以看出，所有五个表示职业特性的变量与u1有大致相同的相关系数，u1视为形容职业特性的指标。第一对典型变量的第二个成员v1y1y2y5y6有较大的相关系而u1和v1之间的相关系数0.5537。u1和v1解释的本组原始变量的比率vum0.5818，vu

X组的原始变量被u1到u5解释了10％，Y组的原始变量被v1到5解释了803％。的：loadr.txt s12=r(1:n1,n1+1:end);%提出X与Y的相关系数s21=s12';%提出Y与X的相关系数s2=r(n1+1:end,n1+1:end);%提出Y与Y的相关系数[vec1,val1]=eig(m1);%求M1的特征向量和特征值forvec1(:,i)=vec1(:,i)/sqrt(vec1(:,i)'*s1*vec1(:,i));%特征向量归一化，满足a's1a=1 flag=1;%把计算结果写到Excel中的行计数变量 forvec2(:,i)=vec2(:,i)/sqrt(vec2(:,i)'*s2*vec2(:,i));%特征向量归一化，满足b's2b=1 %y,v的相关系数 %x,v的相关系数 flag=flag+2;str=char(['A',int2str(flag)]); WTO，作为区域中心的城市在区域经济发展中的作用69个具体指标；沈正平、马晓冬、戴先杰和翟仁祥（2002）构建了测度城市竞争力的指标体系，并用因子分析、聚类分析等方法对新亚欧大陆桥统计分析中，我们用简单相关系数反映两个变量之间的线性相关关系。1936年Holing关本思想是仿照主成分分析法中把多变量与多变量之间的相关化为两个变量之间相关的做x11 Z

y1qx1x1·x2y·2q · · nq

R12 第一步，计算相关系数阵R，并将R剖分为R ，其中R11，R22 22 RT为第一组与第二组变量的相

RRRRRR

的特2iaMR1RR1R的特征根2，特征向量b 222111 uaTX,vbTY;uaTX,vbTY;·;uaTX,vbTY(tmin(p, 记Uuu,·u)T，Vvv,·v)T 第三步，典型相关系数 的显著性检验典型结构指原始变量与典型变量之间的相关系数阵RX,U)，RX,V)R(Y,UR(Y,V，据此可以计算任一个典型变量uk或vkX(或Y)变差的百分RdX;uk（Rd(Yvk)。同时可求得前t个典型变量u1,·ut(或v1,·vtXY)总变差的累计百分比RdX;u1,·ut)或Rd(Yv1,·vt)。组典型变量uk解释第二组变量Y总变差的百分比Rd(Y;uk)（或第二组中典型变量解释城市竞争要取决于产业经济效益、对外开放程度、基础设施、市民素质、政体系：市场占有率、GDP增长率、劳动生产率和居民人均收入。城市基础设施指标体本设施指数（由城市能源、交通、道路、住房等具体指标综合而成每百人拥有 （由医院个数、万人医院床位数综合构成据如表28、表29。28y2y3y3--------200329数数-------------------------------------------------------2003 301234其结果如表31所示。3113.7608e-23.3963e-3843注：表中的e-007表示107典型相关模型，如下表32所示：321u0.1535x*0.3423x*0.4913x*0.3372x*0.1149x* v0.1395y*0.7185y*0.427y*0.0285 2u0.2134x*0.2637x*0.3953x*0.869x*0.2429x* v0.1322y*0.7361y*0.772y*0.0059 32第一组典型相关方程可知，基础设施方面的主要因素是x2,x3,x4数分别为0.34230.49130.3372设（2（3（4变量v1与y2呈高度相关，说明在城市竞争力中，市场占有率（y2）占有主要地位。根x4（技术设施指数）是基础设施方面的主要因素，而居民人0.869型变量占有信息量较大，所以总体上基础设施方面的主要因素按重要程度依次是x3x2x4，反映城市竞争力的主要指标是y2y3。33结构分析(相关系数----------由表33x1x2x3x4与“基础设施组”的第一典型变量u1均呈高度相关，说明对外x4与基础设施组的第二典型变量和竞争力组的第二典型变在反映城市竞争力中占有主导地位。y3与v1呈较高相关，与v2呈高相关，但v2凝聚的y2则与“影响组”的第一典型变量也呈高度相关。这种一致基础设施方面的u2与x5相关系数为正值（0.4688，而典型系数却为负值（-0.2429。x5为抑制变量(Suppressor)。由表33的相关系数还可以看出，“影响组”的第一典型变量u1对y2（市场占有率）有相当高的预测能力，相关系数为0.8137，而对y4（长期经济增长率）预测能力较差，相关系数仅为0.1625。表34被典型变量解释的X表35被典型变量解释的Y3435可以看出，两对典型变量u1、u2和v1、v2均较好地预测了对应典型变量u1、u2解释的比例和为64.04%；来自“基础设施组”的方差被“竞争力组”典型变量v1、v2解释的方差比例和为56.81%。城市竞争力变量组被其自身及其对立典市场占有率是企业竞争力大小的最直接表现，它反映一个城市的产品在全部城市占主每人数人市、、城市居民人均收入和长期经济增长率综合反映了城市在域内和域外创造价值的状PP增长率反映了城市价值扩展的速度和潜力。因此，居民人均收入可以综合反映出竞争已成为竞争的主要对象和，占有人才便控制了城市竞争的制高点，也就决定了设施等决定着城市利用资源和对人才的。因此，城市基础设施中的对内设施建设数和每百人数与居民人均收入和长期经济增长率反方向增长，设施和方面的投资在一定程度上影响了城市利用资源、创造价值的水平。因为设施和投资必然要占用GDP的增长。昌，2002计算的程序如下loadx.txt%原始的x组的数据保存在纯文本文件x.txt中loady.txt%原始的y组的数据保存在纯文本文件y.txt中x=zscore(x);y=zscore(y);%标准化数据n=size(x,1);%%下面做典型相关分析，a,b返回的是典型变量的系数，r返回的是典型相关系%下面修正a,b每一列的正负号,使得a,b每一列的系数和为 %典型相关系数的平方，M1或M2矩阵的非零特征值§9对应分析 服样品量大时作Q型因子分析所带来计算上的，且把R型和Q型因子分析统一起来，Z将两者有ZTZZZT有相同的非零特征值，记为12m0SR对应的标准化特征向量为vi，则SQ的特征值i对应的标准化特征向1ui

由SR的特征值和特征向量即可写出R型因子分析的因子载荷矩阵（记为AR）和型因子分析的因子载荷矩阵（AQ Av

．

v v vp ． mu

mA

m

u u ;

; 1 un1

un2

·

m由于SRSQ具有相同的非零特征值，而这些特征值又正是各个公共因子的方差，设有np ·x1pX21

·x2p· ;· n npP1X(p ijn T其中Txij，pij xij（i1,2,·,n，j1,2,·,p。不难看出