多元統(tǒng)計分析建模

目錄一二三四五聚類分析主成分分析典型相關分析歷史建模試題數(shù)據(jù)處理方法Box-Cox變換是統(tǒng)計建模中常用的一種數(shù)據(jù)變換，用于連續(xù)的響應變量不滿足正態(tài)分布的情況。比如在使用線性回歸的時候，由于殘差不符合正態(tài)分布而不滿足建模的條件，這時候要對響應變量進行變換，把數(shù)據(jù)變成正態(tài)的。1.boxcox變換(x>0)其中x為原始數(shù)據(jù)，y為變換后的數(shù)據(jù).數(shù)據(jù)處理方法一在MATLAB中，上述變換的命令如下：[t,l]=boxcox(x)其中x是原始數(shù)據(jù)(列向量)，t是變換以后的數(shù)據(jù)，l是變換公式中參數(shù)

的數(shù)值例1.1949—1991淮河流域成災面積說明如何利用上述的變換使得數(shù)據(jù)從不具備正態(tài)分布到符合正態(tài)分布.圖11949—1991淮河流域成災面積qq圖QQ圖的作用用于直觀驗證一組數(shù)據(jù)是否來自某個分布，或者驗證某兩組數(shù)據(jù)是否來自同一（族）分布。在教學和軟件中常用的是檢驗數(shù)據(jù)是否來自于正態(tài)分布。把已知分布的分位數(shù)標在縱軸上，樣本分位數(shù)標在橫軸上，從圖形上可以了解到分布的信息。圖形是直線說明是正態(tài)分布。圖形中有一段是直線，在兩端存在弧度，說明峰度的情況。圖形是曲線圖，說明不對稱。如果Q-Q圖是直線，當該直線成45度角并穿過原點時，說明分布與給定的正態(tài)分布完全一樣。如果是成45度角但不穿過原點，說明均值與給定的正態(tài)分布不同，如果是直線但不是45度角，說明均值與方差都與給定的分布不同。如果Q-Q圖中間部分是直線，但是右邊在直線下面，左邊在直線上面，說明分布的峰度大于3，反之說明峰度小于3.從圖上可以看出散點并不聚集在直線上，因此流域成災面積（原始數(shù)據(jù)）不服從正態(tài)分布，這一點也可以通過jbtest檢驗來證實.但是通過變換以后的圖形如圖2所示，顯然數(shù)據(jù)服從正態(tài)分布.圖2流域成災面積（變換后數(shù)據(jù)）qq圖下面以人均糧食支出x1為例，進行正態(tài)檢驗。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)數(shù)據(jù)類型hp是否拒絕正態(tài)分布原始數(shù)據(jù)10.0112是變換數(shù)據(jù)00.9766否表2.正態(tài)檢驗例2做出x1,x3的散點圖，可否區(qū)別開？解：輸入原始數(shù)據(jù)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年人均消費數(shù)據(jù)練習boxcox變換與正態(tài)分布檢驗2.其他變換公式標準化：設有數(shù)據(jù)x=(x1,x2,…,xn)MATLAB命令：y=zscore(x)規(guī)格化：設有數(shù)據(jù)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均值聚類定義目標函數(shù)為顯然J(U,V)表示了各類中樣本到聚類中心的加權距離平方和，權重是樣本xk對第i類隸屬度的m次方，聚類準則取為求的極小值：(min)｛J(U,V)｝。其中聚類中心為：其中在Matlab中(m=2)，我們只要直接調(diào)用如下程序即可：[center,U,fcn]=fcm(data,cluster_n)data：要聚類的數(shù)據(jù)集合，每一行為一個樣本；cluster_n：聚類數(shù)（大于1）。Center：最終的聚類中心矩陣，其每一行為聚類中心的坐標值；U：最終的模糊分區(qū)矩陣；fcn：在迭代過程中的目標函數(shù)值注意：使用上述方法時，要根據(jù)中心坐標center的特點分清楚每一類中心代表的是實際中的哪一類，才能準確地將待聚類的各方案準確地分為各自所屬的類別.

例3.2007年安徽省各地市工業(yè)企業(yè)效益指標如下表所示，請利用模糊C均值聚類方法分為三類。地區(qū)工業(yè)總產(chǎn)值工業(yè)增加值實收資本業(yè)務收入業(yè)務成本利潤總額合肥市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.安徽工業(yè)企業(yè)數(shù)據(jù)b=[39.63,…];%輸入數(shù)據(jù)作為一個矩陣[center,U,fcn]=fcm(b,3)；%模糊C均值聚類解：Matlab中計算的程序如下得到輸出的結(jié)果為：

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列，每一列表示一個地區(qū)關于三類的隸屬度，每一列最大值所在的行數(shù)即為該地區(qū)的類別。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類的樣本號碼(即原始數(shù)據(jù)行號)表4.2001-2005年災情數(shù)據(jù)年份死亡人口（人）緊急轉(zhuǎn)移安置人口（萬人）倒塌房屋（萬間）直接經(jīng)濟損失（億元）2005年24751570.3226.42042.12004年2250563.3155.01602.32003年2259707.3343.01884.22002年2384471.8189.51637.22001年2538211.192.21942.2上網(wǎng)查找2005年以后的數(shù)據(jù)，然后分為三類：災害較重，災害一般，災害較輕.主成分分析三主成分分析的主要目的是希望用較少的變量去解釋原來資料中的大部分變異，將我們手中許多相關性很高的變量轉(zhuǎn)化成彼此相互獨立或不相關的變量。通常是選出比原始變量個數(shù)少，能解釋大部分資料中的變異的幾個新變量，即所謂主成分，并用以解釋資料的綜合性指標。由此可見，主成分分析實際上是一種降維方法。主成分分析用于投資組合風險管理，企業(yè)效益的綜合評價，圖像特征識別，機械加工或傳感器故障檢測，災害損失分析。如果將主成分分析技術與聚類分析、判別分析以及回歸分析方法相結(jié)合則可以解決更多實際問題。1.主成分的基本思想設x1,x2,…,xp，為p個n維隨機變量（p項指標）在統(tǒng)計學中，經(jīng)常使用原始指標的線性組合所構(gòu)成的綜合指標來代替原有的指標，即要求Yi盡可能地反映原有P個變量的信息.這里的‘信息’用Yi方差來度量，即要求var(Yi)=l1T

l1達到最大，為此我們需要對系數(shù)向量加以限制即滿足約束條件：求l1使var(Y1)取最大值，由此l1所確定的隨機變量Y1稱為隨機變量的第一主成分.

如果第一主成分Y1還不足以反映原變量的信息，則進一步求Y2,為了使Y1和Y2所反映原變量的信息不相重疊，要求Y1和Y2不相關，即于是,在約束條件下，求l2使得var(Y2)達到最大，由此所確定的隨機變量Y2稱為的第二主成分.一般地，求第i個主成分Yi，則要求其系數(shù)及主成分滿足以下條件：（1）系數(shù)向量是單位向量，即

（2）不同的主成分不相關，沒有重疊信息,即

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

Y1,Y2,…,Yp依次稱為第一主成分，第二主成分，…，第p個主成分.①無量綱化2.主成分分析的計算步驟設有n個樣本，p項指標的數(shù)據(jù)矩陣令則實現(xiàn)無量綱化在MATLAB中：Y=X./[ones(n,1)*std(X)]②構(gòu)造矩陣Y的實對稱矩陣，通常用協(xié)方差矩陣

或相關系數(shù)矩陣R.③計算

或R的特征值與相應的特征向量；④根據(jù)特征值計算累計貢獻率（85%），確定主成分的個數(shù)，而特征向量就是主成分的系數(shù)向量.MATLAB:

=cov(Y),R=corrcoef(Y)[V,D]=eig(R)%D為對角矩陣,主對角為特征值，V的每一列為特征值對應的特征向量.⑤計算主成分的數(shù)值（即主成分得分）注意：利用第一主成分得分排序要滿足兩個條件：最大特征值對應的特征向量是正向量；貢獻率>50%地區(qū)

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.根據(jù)x1工業(yè)總產(chǎn)值，x2工業(yè)銷售產(chǎn)值，x3流動資產(chǎn)年平均余額，x4固定資產(chǎn)凈值年平均余額，x5業(yè)務收入,x6利潤總額等六項指標進行主成分分析.(1)選取指標是否合適？（2）給出各市大中型工業(yè)企業(yè)排名。表5.安徽工業(yè)數(shù)據(jù)解：首先輸入數(shù)據(jù)A=[data];%data即表中數(shù)據(jù)R=corrcoef(A);得到的相關系數(shù)矩陣為：由于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.特征值、特征向量及貢獻率

地區(qū)得分排名地區(qū)排名得分地區(qū)得分排名合肥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.根據(jù)軟件輸出結(jié)果，寫出第一、第二主成分的公式；2.將各地區(qū)分成三類，比較主成分排名與分類是否具有一致性？3.因子載荷矩陣例4得到第一主成分公式為F1=0.4595x2+0.4552x3+0.4158x4+0.46x5+0.4441x6我們稱主成分Yi與指標Xj的相關系數(shù)為Yi在Xj上的因子載荷量.因子載荷用于解釋第j個變量對第i個主成分的重要程度.計算可得第一主成分與五個指標的載荷分別為：0.98670.97730.89270.98760.9535典型相關分析三四在實際問題中，經(jīng)常遇到研究兩組隨機變量之間的相關性.比如工廠管理人員需要了解原料的主要質(zhì)量指標與產(chǎn)品的主要質(zhì)量指標之間的相關性，以便提高產(chǎn)品質(zhì)量；醫(yī)生要根據(jù)一組化驗指標確定與一些疾病之間的關系；主教練排兵布陣要考慮自己的隊員與對手之間的相生相克以便制定更好的對策，等等.受主成分分析的啟發(fā)，對每組變量分別構(gòu)造線性組合，將兩組變量之間的相關性轉(zhuǎn)化為兩個變量之間的相關性進行研究.典型相關分析示意圖如圖4.7所示.圖4.7典型相關分析示意圖1.總體典型變量的定義設有兩組隨機變量(XT,YT)T=(X1,X2,…,Xp,Y1,Y2,…,Yq)T的協(xié)方差矩陣為

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

T21=cov(X,Y)注意：的維數(shù)p+q階方陣,11p階方陣,22q階方陣根據(jù)典型相關的思想，分別考慮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的相關系數(shù)為在約束條件a1T

11a1=b1T

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

u1,v1=a1T

12b1取得最大值.如果(U1,V1)還不足以反映X,Y之間的相關性，還可構(gòu)造第二對線性組合：使得(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之間的相關性，還可構(gòu)造第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個典型相關系數(shù).2.總體典型變量與典型相關系數(shù)的計算(1)計算矩陣(XT,YT)T的協(xié)方差矩陣(2)令求A,B的特征值

12,22…,

p2與對應的正交單位特征向量ek,fkk=1,…,p(3)X,Y的第k對典型相關變量為(4)X,Y的第k個典型相關系數(shù)為：

k(k=1,2,…p)3.典型相關分析的Matlab實現(xiàn)設X=(xij)np，Y=(yij)nq是取自總體的觀測數(shù)據(jù)，典型相關分析的步驟如下①輸入數(shù)據(jù)并計算協(xié)方差矩陣

a=[X,Y];%此前X,Y的數(shù)據(jù)應該已經(jīng)輸入[n,m]=size(a);R=cov(a);②計算典型相關系數(shù)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);%典型相關系數(shù)③計算典型相關向量X=X./[ones(n,1)*std(X)];%n為具體的樣本容量數(shù)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④典型相關系數(shù)的顯著性檢驗統(tǒng)計量其中檢驗程序如下：D=1-d;f1=fliplr(D');%矩陣左右翻轉(zhuǎn)f2=cumprod(f1);%向量累積乘積d1k=(p-k+1).*(q-k+1);Qk=-[n-0.5*(p+q+3)].*(log(fliplr(f2)));1-chi2cdf(Qk,d1k)注意：如果使用下面的命令，則③、④可一步實現(xiàn)[A,B,r,U,V,stats]=canoncorr(X,Y)其中r是三個典型相關系數(shù)，stats包括Wilks、chisq及F統(tǒng)計量以及相應的概率.例4.選取1980-2008年安徽省人均糧食總產(chǎn)量（噸/人）、人均農(nóng)業(yè)總產(chǎn)值(億元/萬人)、，人均糧食播種面積（千公頃/萬人）、人均農(nóng)業(yè)機械總動力（千瓦/人）、單位面積化肥施用（萬噸/千公頃）、人均受災面積（千公頃/萬人）以及農(nóng)業(yè)生產(chǎn)資料價格指數(shù)指標，分別記為：x1,x2,x3,y1,y2,y3,y4(1)按年度將安徽省糧食生產(chǎn)分成三類；(2)對安徽省糧食生產(chǎn)影響因素進行典型相關分析年份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