多元统计分析建模

目录一二三四五聚类分析主成分分析典型相关分析历史建模试题数据处理方法Box-Cox变换是统计建模中常用的一种数据变换，用于连续的响应变量不满足正态分布的情况。比如在使用线性回归的时候，由于残差不符合正态分布而不满足建模的条件，这时候要对响应变量进行变换，把数据变成正态的。1.boxcox变换(x>0)其中x为原始数据，y为变换后的数据.数据处理方法一在MATLAB中，上述变换的命令如下：[t,l]=boxcox(x)其中x是原始数据(列向量)，t是变换以后的数据，l是变换公式中参数

的数值例1.1949—1991淮河流域成灾面积说明如何利用上述的变换使得数据从不具备正态分布到符合正态分布.图11949—1991淮河流域成灾面积qq图QQ图的作用用于直观验证一组数据是否来自某个分布，或者验证某两组数据是否来自同一（族）分布。在教学和软件中常用的是检验数据是否来自于正态分布。把已知分布的分位数标在纵轴上，样本分位数标在横轴上，从图形上可以了解到分布的信息。图形是直线说明是正态分布。图形中有一段是直线，在两端存在弧度，说明峰度的情况。图形是曲线图，说明不对称。如果Q-Q图是直线，当该直线成45度角并穿过原点时，说明分布与给定的正态分布完全一样。如果是成45度角但不穿过原点，说明均值与给定的正态分布不同，如果是直线但不是45度角，说明均值与方差都与给定的分布不同。如果Q-Q图中间部分是直线，但是右边在直线下面，左边在直线上面，说明分布的峰度大于3，反之说明峰度小于3.从图上可以看出散点并不聚集在直线上，因此流域成灾面积（原始数据）不服从正态分布，这一点也可以通过jbtest检验来证实.但是通过变换以后的图形如图2所示，显然数据服从正态分布.图2流域成灾面积（变换后数据）qq图下面以人均粮食支出x1为例，进行正态检验。x1=[8.359.258.197.739.429.1610.069.099.418.706.938.679.986.778.147.677.907.188.826.2510.607.2713.4510.857.21,...7.687.787.948.2812.47];[h,p]=jbtest(x1)[y,t]=boxcox(x1');[h,p]=jbtest(y)数据类型hp是否拒绝正态分布原始数据10.0112是变换数据00.9766否表2.正态检验例2做出x1,x3的散点图，可否区别开？解：输入原始数据a=[8.35 23.53 7.51 8.62 17.42 10 1.04 11.21,…,12.47 76.39 5.52 11.24 14.52225.4625.5];

plot(a(:,3),'+'),holdon,plot(a(:,1),'or'),legend('x3','x1')图3x1,x3的散点图[y3,t3]=boxcox(a(:,3));[y1,t1]=boxcox(a(:,1));plot(y3,'+'),holdon,plot(y1,'or'),legend('y3','y1')图4变换后散点图练习：对1991年人均消费数据练习boxcox变换与正态分布检验2.其他变换公式标准化：设有数据x=(x1,x2,…,xn)MATLAB命令：y=zscore(x)规格化：设有数据x=(x1,x2,…,xn)[0,1]之间MATLAB:y=(x-min(x))./range(x)如果x为m行n列矩阵，列为指标，则命令为：y=(x-ones(m,1)*min(x))./[ones(m,1)*range(x)]此时y仍然是m行n列的一个矩阵.模糊交集与模糊并集设x=(x1,x2,…xn),y=(y1,y2,…,yn),xi,yi[0,1]爱因斯坦积：爱因斯坦和：注意：zi,wi仍然[0,1]，且有

zi<xi<wi，zi<yi<wi二聚类方法1.模糊C均值聚类定义目标函数为显然J(U,V)表示了各类中样本到聚类中心的加权距离平方和，权重是样本xk对第i类隶属度的m次方，聚类准则取为求的极小值：(min)｛J(U,V)｝。其中聚类中心为：其中在Matlab中(m=2)，我们只要直接调用如下程序即可：[center,U,fcn]=fcm(data,cluster_n)data：要聚类的数据集合，每一行为一个样本；cluster_n：聚类数（大于1）。Center：最终的聚类中心矩阵，其每一行为聚类中心的坐标值；U：最终的模糊分区矩阵；fcn：在迭代过程中的目标函数值注意：使用上述方法时，要根据中心坐标center的特点分清楚每一类中心代表的是实际中的哪一类，才能准确地将待聚类的各方案准确地分为各自所属的类别.

例3.2007年安徽省各地市工业企业效益指标如下表所示，请利用模糊C均值聚类方法分为三类。地区工业总产值工业增加值实收资本业务收入业务成本利润总额合肥市1099.82356.03191.241020.77823.4543.72淮北市239.27112.4279.80266.74202.254.07亳州市116.2344.7114.4165.8149.173.34宿州市109.0640.6414.18112.3497.32-1.10蚌埠市218.3882.2258.42202.28150.3512.39阜阳市192.0366.0033.02183.78133.8118.04淮南市300.09141.61120.09310.51239.8319.36滁州市248.4387.1747.61230.67185.7116.26六安市136.2457.8817.81128.5094.2610.15马鞍山687.38260.28180.80761.63653.1735.29巢湖市172.2852.7536.31166.52114.5812.78芜湖市674.33166.80100.15648.34566.8336.52宣城市160.9240.7922.72151.54126.449.89铜陵市513.95151.6963.52651.15571.0027.90池州市28.0111.268.2125.2119.143.68安庆市363.2676.7157.50372.97327.6115.90黄山市32.748.809.2231.4125.672.22表3.安徽工业企业数据b=[39.63,…];%输入数据作为一个矩阵[center,U,fcn]=fcm(b,3)；%模糊C均值聚类解：Matlab中计算的程序如下得到输出的结果为：

center=

746.0614237.5902141.0734769.5504653.425036.2240281.0558104.409275.5785287.9557230.567613.9213118.704341.051719.6115107.827883.16356.7159此时，center的每一行就是每一类最终的中心坐标，由效益型指标可知：第一行表示效益最好的一类，第三行表示效益最差的一类，第二行则介于两者之间。U=

0.80820.00410.00220.00030.01220.01030.00710.00800.00100.11540.94790.02480.00470.58270.31440.95790.85230.01610.07640.04800.97300.99500.40510.67530.03510.13960.9829

0.98610.00640.88960.00470.70690.01170.04670.01050.00910.14820.07520.09870.20460.09440.84570.08720.00480.84540.03530.89660.08850.89390.10760.9023U共有17列，每一列表示一个地区关于三类的隶属度，每一列最大值所在的行数即为该地区的类别。Y=sort(U);[Y,I]=sort(U);%排序L1=find(I(3,:)==1),L2=find(I(3,:)==2),L3=find(I(3,:)==3),%检索下面给出MATLAB检索类别的方法：此时，Li输出的是第i类的样本号码(即原始数据行号)表4.2001-2005年灾情数据年份死亡人口（人）紧急转移安置人口（万人）倒塌房屋（万间）直接经济损失（亿元）2005年24751570.3226.42042.12004年2250563.3155.01602.32003年2259707.3343.01884.22002年2384471.8189.51637.22001年2538211.192.21942.2上网查找2005年以后的数据，然后分为三类：灾害较重，灾害一般，灾害较轻.主成分分析三主成分分析的主要目的是希望用较少的变量去解释原来资料中的大部分变异，将我们手中许多相关性很高的变量转化成彼此相互独立或不相关的变量。通常是选出比原始变量个数少，能解释大部分资料中的变异的几个新变量，即所谓主成分，并用以解释资料的综合性指标。由此可见，主成分分析实际上是一种降维方法。主成分分析用于投资组合风险管理，企业效益的综合评价，图像特征识别，机械加工或传感器故障检测，灾害损失分析。如果将主成分分析技术与聚类分析、判别分析以及回归分析方法相结合则可以解决更多实际问题。1.主成分的基本思想设x1,x2,…,xp，为p个n维随机变量（p项指标）在统计学中，经常使用原始指标的线性组合所构成的综合指标来代替原有的指标，即要求Yi尽可能地反映原有P个变量的信息.这里的‘信息’用Yi方差来度量，即要求var(Yi)=l1T

l1达到最大，为此我们需要对系数向量加以限制即满足约束条件：求l1使var(Y1)取最大值，由此l1所确定的随机变量Y1称为随机变量的第一主成分.

如果第一主成分Y1还不足以反映原变量的信息，则进一步求Y2,为了使Y1和Y2所反映原变量的信息不相重叠，要求Y1和Y2不相关，即于是,在约束条件下，求l2使得var(Y2)达到最大，由此所确定的随机变量Y2称为的第二主成分.一般地，求第i个主成分Yi，则要求其系数及主成分满足以下条件：（1）系数向量是单位向量，即

（2）不同的主成分不相关，没有重叠信息,即

（3）各主成分的方差递减，重要性递减,即

Y1,Y2,…,Yp依次称为第一主成分，第二主成分，…，第p个主成分.①无量纲化2.主成分分析的计算步骤设有n个样本，p项指标的数据矩阵令则实现无量纲化在MATLAB中：Y=X./[ones(n,1)*std(X)]②构造矩阵Y的实对称矩阵，通常用协方差矩阵

或相关系数矩阵R.③计算

或R的特征值与相应的特征向量；④根据特征值计算累计贡献率（85%），确定主成分的个数，而特征向量就是主成分的系数向量.MATLAB:

=cov(Y),R=corrcoef(Y)[V,D]=eig(R)%D为对角矩阵,主对角为特征值，V的每一列为特征值对应的特征向量.⑤计算主成分的数值（即主成分得分）注意：利用第一主成分得分排序要满足两个条件：最大特征值对应的特征向量是正向量；贡献率>50%地区

x1

x2x3x4x5x6合肥市1932.271900.53653.83570.951810.70119.53淮北市367.05366.08186.16252.07395.4332.82亳州市86.8985.3840.8551.7183.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.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工业总产值，x2工业销售产值，x3流动资产年平均余额，x4固定资产净值年平均余额，x5业务收入,x6利润总额等六项指标进行主成分分析.(1)选取指标是否合适？（2）给出各市大中型工业企业排名。表5.安徽工业数据解：首先输入数据A=[data];%data即表中数据R=corrcoef(A);得到的相关系数矩阵为：由于r12=r21=1，表明指标x1,x2完全线性相关，故只需保留一个指标.A=A(:,2:6)./[ones(17,1)*std(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(-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淮南50.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.将各地区分成三类，比较主成分排名与分类是否具有一致性？3.因子载荷矩阵例4得到第一主成分公式为F1=0.4595x2+0.4552x3+0.4158x4+0.46x5+0.4441x6我们称主成分Yi与指标Xj的相关系数为Yi在Xj上的因子载荷量.因子载荷用于解释第j个变量对第i个主成分的重要程度.计算可得第一主成分与五个指标的载荷分别为：0.98670.97730.89270.98760.9535典型相关分析三四在实际问题中，经常遇到研究两组随机变量之间的相关性.比如工厂管理人员需要了解原料的主要质量指标与产品的主要质量指标之间的相关性，以便提高产品质量；医生要根据一组化验指标确定与一些疾病之间的关系；主教练排兵布阵要考虑自己的队员与对手之间的相生相克以便制定更好的对策，等等.受主成分分析的启发，对每组变量分别构造线性组合，将两组变量之间的相关性转化为两个变量之间的相关性进行研究.典型相关分析示意图如图4.7所示.图4.7典型相关分析示意图1.总体典型变量的定义设有两组随机变量(XT,YT)T=(X1,X2,…,Xp,Y1,Y2,…,Yq)T的协方差矩阵为

11=cov(X)，22=cov(Y),12=

T21=cov(X,Y)注意：的维数p+q阶方阵,11p阶方阵,22q阶方阵根据典型相关的思想，分别考虑X,Y的线性组合其中a1=(a11,a12,…,a1p)T,b1=(b11,b12,…,b1q)TVar(U1)=a1T

11a1,

Var(V1)=b1T

22b1,cov(U1,V1)=a1T

12b1由于则U1,V1的相关系数为在约束条件a1T

11a1=b1T

22b1=1下，求a1,b1，使得

u1,v1=a1T

12b1取得最大值.如果(U1,V1)还不足以反映X,Y之间的相关性，还可构造第二对线性组合：使得(U1,V1)与(U2,V2)不相关，即cov(u1,u2)=cov(u1,v2)=cov(u2,v1)=cov(v1,v2)=0在约束条件Var(u1)=Var(v1)=Var(u2)=Var(v2)=1下求a2,b2，使得

u2,v2=a2T

12b2取得最大值.一般地，若前k-1对典型变量还不足以反映X,Y之间的相关性，还可构造第k对线性组合：在约束条件Var(uk)=Var(vk)=1,及cov(uk,uj)=cov(uk,vj)=cov(vk,uj)=cov(vk,vj)=0,(1j<k)求ak,bk，使得

uk,vk=akT

12bk取得最大值.如此确定的(uk,vk)称为X,Y的第k对典型变量，相应的

uk,vk称为第k个典型相关系数.2.总体典型变量与典型相关系数的计算(1)计算矩阵(XT,YT)T的协方差矩阵(2)令求A,B的特征值

12,22…,

p2与对应的正交单位特征向量ek,fkk=1,…,p(3)X,Y的第k对典型相关变量为(4)X,Y的第k个典型相关系数为：

k(k=1,2,…p)3.典型相关分析的Matlab实现设X=(xij)np，Y=(yij)nq是取自总体的观测数据，典型相关分析的步骤如下①输入数据并计算协方差矩阵

a=[X,Y];%此前X,Y的数据应该已经输入[n,m]=size(a);R=cov(a);②计算典型相关系数R1=inv(R(1:p,1:p))*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④典型相关系数的显著性检验统计量其中检验程序如下：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及F统计量以及相应的概率.例4.选取1980-2008年安徽省人均粮食总产量（吨/人）、人均农业总产值(亿元/万人)、，人均粮食播种面积（千公顷/万人）、人均农业机械总动力（千瓦/人）、单位面积化肥施用（万吨/千公顷）、人均受灾面积（千公顷/万人）以及农业生产资料价格指数指标，分别记为：x1,x2,x3,y1,y2,y3,y4(1)按年度将安徽省粮食生产分成三类；(2)对安徽省粮食生产影响因素进行典型相关分析年份X1x2x3y1y2y3y419800.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.18720.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.01740.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.200019971.04780.26213.15500.83730.02850.224798.900019980.95330.25013.15150.93710.02960.142494.800019991.01720.25943.15021.01540.02980.333695.300020000.88360.24143.00