多元统计分析建模

1、统计与应用数学学院统计与应用数学学院STATISTICS & APPLIED MATHEMATICS多元统计分析建模方法多元统计分析建模方法 李柏年李柏年 刘德志刘德志一一二二三三四四五五聚类分析聚类分析主成分分析主成分分析典型相关分析典型相关分析历史建模试题历史建模试题数据处理方法数据处理方法 Box-Cox变换是统计建模中常用的一种变换是统计建模中常用的一种数据变换，用于连续的响应变量不满足正态数据变换，用于连续的响应变量不满足正态分布的情况。比如在使用线性回归的时候，分布的情况。比如在使用线性回归的时候，由于残差不符合正态分布而不满足建模的条由于残差不符合正态分布而不满足建模的条

2、件，这时候要对响应变量进行变换，把数据件，这时候要对响应变量进行变换，把数据变成正态的。变成正态的。 1. boxcox变换变换0 xlog01xy )(/ )( (x0)其中其中x为原始数据，为原始数据，y为变换后的数据为变换后的数据.数据处理方法数据处理方法一一在在MATLAB中，上述变换的命令如下：中，上述变换的命令如下： t,l=boxcox(x)其中其中 x是原始数据是原始数据(列向量列向量)，t是变换以后的数据，是变换以后的数据，l是变换公式中参数是变换公式中参数 的数值的数值例例1. 19491991淮河流域成灾面积说明如何利用上述淮河流域成灾面积说明如何利用上述的变换使得数据从

3、不具备正态分布到符合正态分布的变换使得数据从不具备正态分布到符合正态分布.图图1 19491991淮河流域成灾面积淮河流域成灾面积qq图图 QQ图的作用用于直观验证一组数据是否来自图的作用用于直观验证一组数据是否来自某个分布，或者验证某两组数据是否来自同一某个分布，或者验证某两组数据是否来自同一（族）分布。在教学和软件中常用的是检验数据（族）分布。在教学和软件中常用的是检验数据是否来自于正态分布。是否来自于正态分布。 把已知分布的分位数标在纵轴上，样本分位数把已知分布的分位数标在纵轴上，样本分位数标在横轴上，从图形上可以了解到分布的信息。标在横轴上，从图形上可以了解到分布的信息。 图形是直线说

4、明是正态分布。图形中有一段图形是直线说明是正态分布。图形中有一段是直线，在两端存在弧度，说明峰度的情况。图是直线，在两端存在弧度，说明峰度的情况。图形是曲线图，说明不对称。形是曲线图，说明不对称。 如果如果Q-Q图是直线，当该直线成图是直线，当该直线成45度角并穿过度角并穿过原点时，说明分布与给定的正态分布完全一样。原点时，说明分布与给定的正态分布完全一样。如果是成如果是成45度角但不穿过原点，说明均值与给定度角但不穿过原点，说明均值与给定的正态分布不同，如果是直线但不是的正态分布不同，如果是直线但不是45度角，说度角，说明均值与方差都与给定的分布不同。如果明均值与方差都与给定的分布不同。如果

5、Q-Q图图中间部分是直线，但是右边在直线下面，左边在中间部分是直线，但是右边在直线下面，左边在直线上面，说明分布的峰度大于直线上面，说明分布的峰度大于3，反之说明峰，反之说明峰度小于度小于3.从图上可以看出散点并不聚集在直线上，因此流从图上可以看出散点并不聚集在直线上，因此流域成灾面积（原始数据）不服从正态分布，这一域成灾面积（原始数据）不服从正态分布，这一点也可以通过点也可以通过jbtest检验来证实检验来证实.但是通过变换以后但是通过变换以后的图形如图的图形如图2所示，显然数据服从正态分布所示，显然数据服从正态分布.图图2 流域成灾面积（变换后数据）流域成灾面积（变换后数据）qq图图X1X

6、2X3X4X5X6X7X88.3523.537.518.6217.42101.0411.219.2523.756.619.1917.7710.481.7210.518.1930.54.729.7816.287.62.5210.327.7329.25.429.4319.298.492.52109.4227.98.28.1416.179.421.559.769.1627.989.019.3215.999.11.8211.3510.0628.6410.5210.0516.188.391.9610.819.0928.127.49.6217.2611.122.4912.659.4128.25.7710.

7、816.3611.561.5312.178.728.127.2110.5319.4513.31.6611.966.9329.854.549.4916.6210.651.8813.618.6736.057.317.7516.6711.682.3812.889.9837.697.018.9416.1511.080.8311.676.7738.696.018.8214.7911.441.7413.238.1437.759.618.4913.159.761.2811.287.6735.718.048.3115.137.761.4113.257.939.778.4912.9419.2711.052.04

8、13.297.1840.917.328.9417.612.751.1414.88.8233.77.5910.9818.8214.731.7810.16.2535.024.726.2810.037.151.9310.3910.652.417.79.9812.5311.72.3114.697.2752.653.849.1613.0315.261.9814.5713.4555.855.57.459.559.522.2116.310.8544.687.3214.5117.1312.081.2611.577.2145.797.6610.3616.5612.862.2511.697.6850.3711.3

9、513.319.2514.592.7514.877.7848.44820.5122.1215.731.1516.617.9439.6520.9720.8222.5212.411.757.98.2864.34822.2220.0615.120.7222.8912.4776.395.5211.2414.52225.4625.5表表1. 1991各地区居民月人均消费数据各地区居民月人均消费数据下面以人均粮食支出下面以人均粮食支出x1为例，进行正态检验。为例，进行正态检验。x1=8.35 9.25 8.19 7.73 9.42 9.16 10.06 9.09 9.41 8.70 6.93 8.67 9

10、.98 6.77 8.14 7.67 7.90 7.18 8.82 6.25 10.60 7.27 13.45 10.85 7.21,.7.68 7.78 7.94 8.28 12.47;h,p=jbtest(x1)y,t=boxcox(x1); h,p=jbtest(y)数据类型数据类型hp是否拒绝正态分布是否拒绝正态分布原始数据原始数据10.0112是是变换数据变换数据00.9766否否表表2. 正态检验正态检验例例2 做出做出x1,x3的散点图，可否区别开？的散点图，可否区别开？解：输入原始数据解：输入原始数据a=8.3523.537.51 8.6217.42101.04 11.21,

11、12.4776.395.52 11.24 14.52 22 5.46 25.5;0510152025300510152025 x3x1 plot(a(:,3),+),hold on,plot(a(:,1),or),legend(x3,x1)图图3 x1,x3的散点图的散点图0510152025300.811.21.41.6 y3y1y3,t3=boxcox(a(:,3);y1,t1=boxcox(a(:,1); plot(y3,+),hold on, plot(y1,or),legend(y3,y1)图图4 变换后散点图变换后散点图练习：对练习：对1991年人均消费数据练习年人均消费数据练习b

12、oxcox变变换与正态分布检验换与正态分布检验2. 其他变换公式其他变换公式标准化：设有数据标准化：设有数据x=(x1,x2,xn)niiniiiixxssxxy1211)(,MATLAB命令：命令：y=zscore(x)规格化：设有数据规格化：设有数据x=(x1,x2,xn)0,1之间之间minmaxminiiiiixxxxyMATLAB: y=(x-min(x)./range(x)如果如果x为为m行行n列矩阵，列为指标，则命令为：列矩阵，列为指标，则命令为：y=(x-ones(m,1)*min(x)./ones(m,1)*range(x)此时此时y仍然是仍然是m行行n列的一个矩阵列的一个矩

13、阵.模糊交集与模糊并集模糊交集与模糊并集设设x=(x1,x2,xn),y=(y1,y2,yn),xi,yi 0,1爱因斯坦积：爱因斯坦积：iiiiiiiyxyxyxz2爱因斯坦和：爱因斯坦和：iiiiiyxyxw1注意：注意：zi,wi仍然仍然 0,1，且有，且有 zixiwi，ziyi50%1x2x3x4x5x6x地地 区区 x1 x2 x3 x4 x5 x6合合 肥肥 市市1932.271900.53653.83570.951810.70119.53淮淮 北北 市市367.05366.08186.16252.07395.4332.82亳亳 州州 市市86.8985.3840.8551.71

14、83.268.95宿宿 州州 市市154.27147.0730.6857.96146.30-1.27蚌蚌 埠埠 市市197.21193.28104.5690.15182.607.85阜阜 阳阳 市市244.17231.5556.37121.96224.0426.49淮淮 南南 市市497.74483.69206.80501.37496.5927.76滁滁 州州 市市308.91296.99118.6576.90277.4219.32六六 安安 市市191.77189.0570.1962.31191.9823.08马鞍山市马鞍山市905.32894.61351.52502.991048.0253

15、.88巢巢 湖湖 市市254.99242.38106.6675.48234.7619.65芜芜 湖湖 市市867.07852.34418.82217.76806.9437.01宣宣 城城 市市219.36207.0782.5854.74192.7411.02铜铜 陵陵 市市570.33563.33224.23190.77697.9120.61池池 州州 市市59.1157.3216.9740.3356.566.03安安 庆庆 市市430.58426.25103.08147.05442.040.79黄黄 山山 市市65.0364.3628.388.5860.482.88例例4. 根据根据x1工业

16、总产值，工业总产值，x2工业销售产值，工业销售产值，x3流动资产年平流动资产年平均余额，均余额，x4固定资产净值年平均余额，固定资产净值年平均余额，x5业务收入业务收入,x6利润总利润总额等六项指标进行主成分分析额等六项指标进行主成分分析.(1)选取指标是否合适？（选取指标是否合适？（2）给出各市大中型工业企业排名。）给出各市大中型工业企业排名。表表5. 安徽工业数据安徽工业数据解：首先输入数据解：首先输入数据A=data; % data即表中数据即表中数据R=corrcoef(A);得到的相关系数矩阵为：得到的相关系数矩阵为：1.0000 1.0000 0.9754 0.8231 0.991

17、4 0.93751.0000 1.0000 0.9758 0.8236 0.9920 0.93690.9754 0.9758 1.0000 0.8245 0.9712 0.91270.8231 0.8236 0.8245 R 1.0000 0.8502 0.80200.9914 0.9920 0.9712 0.8502 1.0000 0.92120.9375 0.9369 0.9127 0.8020 0.9212 1.0000由于由于r12=r21=1，表明指标，表明指标x1,x2完全线性相关，故完全线性相关，故只需保留一个指标只需保留一个指标.A=A(:,2:6)./ones(17,1)*s

18、td(A(:,2:6); % 消除量纲消除量纲d,v=eig(corrcoef(A); % 计算特征值与特征向量计算特征值与特征向量w=sum(d)/sum(sum(d); % 计算贡献率计算贡献率F=A-ones(17,1)*mean(A)*v(:,5); % 计算主成分得分计算主成分得分F1,I1=sort(F,descend); % I1给出各名次的序号给出各名次的序号F2,I2=sort(I1); % I2给出各市排名给出各市排名特征值特征值特征向量特征向量贡献率贡献率4.6100(0.4595, 0.4552, 0.4158, 0.4600, 0.4441)0.92200.2475(

19、-0.2517, -0.2103, 0.9054,-0.1315,-0.2354)0.04950.1050(0.1926, 0.3702, -0.0390, 0.3029, -0.8559)0.02100.0322(-0.3510, 0.7779, 0.0275, -0.5153, 0.0738)0.00640.0053(0.7518, -0.0803, 0.0719, -0.6434, -0.0965)0.0011表表6. 特征值、特征向量及贡献率特征值、特征向量及贡献率 地地 区区得分得分排排 名名地地 区区排排 名名得分得分地地 区区得分得分排排 名名合合 肥肥18.671淮淮 南南 5

20、0.642宣宣 城城 -2.64711淮淮 北北 1.6364滁滁 州州 10-1.017铜铜 陵陵 -0.7638亳亳 州州 -3.05412六六 安安 7-0.278池池 州州 -3.62814宿宿 州州 -5.06317马鞍山马鞍山25.774安安 庆庆 -4.65816蚌蚌 埠埠 -3.27113巢巢 湖湖 9-0.952黄黄 山山 -4.24715阜阜 阳阳 0.9326芜芜 湖湖 32.459表表7. 各市第一主成分得分排名各市第一主成分得分排名 练习：练习：1. 根据软件输出结果，写出第一、第二主成分根据软件输出结果，写出第一、第二主成分的公式；的公式；2. 将各地区分成三类，比

21、较主成分排名与分将各地区分成三类，比较主成分排名与分类是否具有一致性？类是否具有一致性？3. 因子载荷矩阵因子载荷矩阵例例4 得到第一主成分公式为得到第一主成分公式为F1=0.4595x2+0.4552x3+0.4158x4+0.46x5+0.4441x6我们称主成分我们称主成分Yi与指标与指标Xj的相关系数的相关系数( ,)( ,)/( )()iji ijijiijjjijijjCov Y XeY XeVar YVar X为为Yi在在Xj上的因子载荷量上的因子载荷量.因子载荷用于解释第因子载荷用于解释第j个变量对第个变量对第i个主成分的重要程度个主成分的重要程度.计算可得第一主成分与五个指标

22、的载荷分别为：计算可得第一主成分与五个指标的载荷分别为： 0.9867 0.9773 0.8927 0.9876 0.9535典型相关分析典型相关分析三三四四 在实际问题中，经常遇到研究两组随机变量在实际问题中，经常遇到研究两组随机变量之间的相关性之间的相关性.比如工厂管理人员需要了解原料的比如工厂管理人员需要了解原料的主要质量指标主要质量指标 与产品的主要质量指标之间的相关与产品的主要质量指标之间的相关性，以便提高产品质量；医生要根据一组化验指性，以便提高产品质量；医生要根据一组化验指标确定与一些疾病之间的关系；主教练排兵布阵标确定与一些疾病之间的关系；主教练排兵布阵要考虑自己的队员与对手之

23、间的相生相克以便制要考虑自己的队员与对手之间的相生相克以便制定更好的对策，等等定更好的对策，等等.受主成分分析的启发，对每组变量分别构造线性受主成分分析的启发，对每组变量分别构造线性组合，将两组变量之间的相关性转化为两个变量组合，将两组变量之间的相关性转化为两个变量之间的相关性进行研究之间的相关性进行研究.典型相关分析示意图如图典型相关分析示意图如图4. 7所示所示.图图4.7 典型相关分析示意图典型相关分析示意图1.总体典型变量的定义总体典型变量的定义设有两组随机变量设有两组随机变量 12(,) ,TpXXXX12( ,)TqYY YY11122122 (XT,YT)T=(X1,X2,Xp,

24、Y1,Y2,Yq)T的协方差矩阵为的协方差矩阵为 11=cov(X)， 22=cov(Y), 12= T21=cov(X,Y)注意：注意： 的维数的维数p+q阶方阵阶方阵, 11 p阶方阵阶方阵, 22 q阶方阵阶方阵根据典型相关的思想，分别考虑根据典型相关的思想，分别考虑X,Y的线性组合的线性组合1111112211111 11221TppTqqUa Xa Xa Xa XVb Yb Yb Yb Y其中其中 a1=(a11,a12,a1p)T, b1=(b11,b12,b1q)TVar(U1)=a1T 11a1, Var(V1)=b1T 22b1,cov(U1,V1)= a1T 12b1由于由

25、于则则U1,V1的相关系数为的相关系数为11112 1,111 1122 1TU VTTabaabb在约束条件在约束条件 a1T 11a1=b1T 22b1=1下，求下，求a1,b1，使得，使得 u1,v1=a1T 12b1取得最大值取得最大值. 如果如果(U1,V1)还不足以反映还不足以反映X,Y之间的相关性，之间的相关性，还可构造第二对线性组合：还可构造第二对线性组合：2221122222221 12222TppTqqUa Xa Xa XaXVb Yb Yb Yb Y使得使得(U1,V1)与与(U2,V2)不相关，即不相关，即cov(u1,u2)=cov(u1,v2)=cov(u2,v1)

26、=cov(v1,v2)=0在约束条件在约束条件Var(u1)=Var(v1)=Var(u2)=Var(v2)=1下下求求a2,b2，使得，使得 u2,v2=a2T 12b2取得最大值取得最大值. 一般地，若前一般地，若前k -1对典型变量还不足以反映对典型变量还不足以反映X,Y之之间的相关性，还可构造第间的相关性，还可构造第k对线性组合：对线性组合：11221 122TkkkkkppTkkkkkqqUa Xa XaXa XVb Yb Yb Yb Y在约束条件在约束条件 Var(uk)=Var(vk)=1,及及cov(uk,uj)=cov(uk,vj)=cov(vk,uj)=cov(vk,vj)

27、=0,(1 jk)求求ak,bk，使得，使得 uk,vk=akT 12bk取得最大值取得最大值.如此确定的如此确定的(uk,vk)称为称为X,Y的第的第k对典型变量，相对典型变量，相应的应的 uk,vk称为第称为第k个典型相关系数个典型相关系数.2. 总体典型变量与典型相关系数的计算总体典型变量与典型相关系数的计算(1) 计算矩阵计算矩阵(XT,YT)T的协方差矩阵的协方差矩阵11122122 AB-1-111122221-1-122211112= = (2) 令令 求求A,B的特征值的特征值 12, 22, p2与对应的正交单位特征向量与对应的正交单位特征向量ek,fkk=1,p(3) X,

28、Y的第的第k对典型相关变量为对典型相关变量为0.511,TkkUeX0.522TkkVfY(4) X,Y的第的第k个典型相关系数为：个典型相关系数为： k(k=1,2,p)3. 典型相关分析的典型相关分析的Matlab实现实现设设X=(xij)n p， Y=(yij)n q是取自总体的观测数据，典是取自总体的观测数据，典型相关分析的步骤如下型相关分析的步骤如下 输入数据并计算协方差矩阵输入数据并计算协方差矩阵 a=X, Y; % 此前此前X,Y的数据应该已经输入的数据应该已经输入n,m=size(a); R=cov (a); 计算典型相关系数计算典型相关系数R1=inv(R(1:p,1:p)*

29、R(1:p,p+1:p+q)*inv(R(p+1:p+q, p+1:p+q)*R(p+1:p+q,1:p);d=sort(eig(R1),descend);p=sqrt(d); %典型相关系数典型相关系数 计算典型相关向量计算典型相关向量X=X./ones(n,1)*std(X); % n为具体的样本容量数为具体的样本容量数Y=Y./ones(n,1)*std(Y);A,B = canoncorr(X,Y);U=(X-ones(n,1)*mean(X)*AV=(Y-ones(n,1)*mean(Y)*B 典型相关系数的显著性检验典型相关系数的显著性检验221(3) ln(1) ( )2pjii

30、jQnpqf 统计量统计量其中其中(1)(1)fpjqj检验程序如下：检验程序如下：D=1-d;f1=fliplr(D); % 矩阵左右翻转矩阵左右翻转f2=cumprod(f1); % 向量累积乘积向量累积乘积d1k=(p-k+1).*(q-k+1);Qk=-n-0.5*(p+q+3).*(log(fliplr(f2);1-chi2cdf(Qk,d1k)注意：如果使用下面的命令，则、可一步实现注意：如果使用下面的命令，则、可一步实现A,B,r,U,V,stats = canoncorr(X,Y)其中其中r是三个典型相关系数，是三个典型相关系数，stats 包括包括Wilks 、chisq及及

31、F统计量以及相应的概率统计量以及相应的概率.例例4. 选取选取1980-2008年安徽省人均粮食总产量（吨年安徽省人均粮食总产量（吨/人）人）、人均农业总产值、人均农业总产值(亿元亿元/万人万人)、，人均粮食播种面、，人均粮食播种面积（千公顷积（千公顷/万人）、人均农业机械总动力（千瓦万人）、人均农业机械总动力（千瓦/人）、单位面积化肥施用（万吨人）、单位面积化肥施用（万吨/千公顷）、人均受千公顷）、人均受灾面积（千公顷灾面积（千公顷/万人）以及农业生产资料价格指数万人）以及农业生产资料价格指数指标，分别记为：指标，分别记为：x1,x2,x3,y1,y2,y3,y4 (1) 按年度将安徽省粮食

32、生产分成三类；按年度将安徽省粮食生产分成三类；(2) 对安徽省粮食生产影响因素进行典型相关分析对安徽省粮食生产影响因素进行典型相关分析年份年份X1 x2x3y1y2y3y419800.87040.04114.63320.39790.00710.2628102.100019811.05380.05684.56640.39290.00910.6130101.700019821.08180.05864.48080.40470.01140.2349101.300019831.08980.06044.26000.41470.01150.1339102.800019841.15760.06644.1872

33、0.41910.01270.4036107.000019851.09830.07364.14700.42230.01390.2131101.700019861.16490.08174.00890.45030.01410.3617102.100019871.16700.09024.02260.49780.01440.2895112.800019881.06610.09923.76960.52970.01550.6894118.600019891.08800.10603.69690.54920.01670.3333121.700019901.09500.11343.61230.56800.0174

34、0.4821103.900019910.74100.08743.47200.58470.01760.4576102.300019920.96280.10783.35270.59700.01910.3754102.500019931.03740.14273.30300.62030.02150.5717112.900019940.92860.19953.24990.66210.02300.1626122.800019951.02330.24613.22280.70830.02430.3637128.000019961.03120.26113.19300.77020.02970.2261107.20

35、0019971.04780.26213.15500.83730.02850.224798.900019980.95330.25013.15150.93710.02960.142494.800019991.01720.25943.15021.01540.02980.333695.300020000.88360.24143.00881.06360.03010.347098.200020010.88620.24382.91911.12170.03180.262897.900020020.97310.25072.95801.18680.03220.613099.900020030.77390.21592.94381.23870.03340.2349100.200020040.94240.28932.96561.30010.03220.13