1、普通高等院校計(jì)算機(jī)課程規(guī)劃教材,MATLAB數(shù)據(jù)分析方法,李柏年 吳禮斌 主編 張孔生 丁 華 參編,第5章 主成分與典型相關(guān)分析,主成分分析就是將原來(lái)指標(biāo)重新組合成一組新的互相無(wú)關(guān)的指標(biāo)來(lái)代替原來(lái)指標(biāo).這些綜合指標(biāo)就是原來(lái)指標(biāo)的線性組合，同時(shí)根據(jù)實(shí)際需要從中選取幾個(gè)較少的綜合指標(biāo)盡可能多地反映原來(lái)指標(biāo)的信息.,5.1主成分分析的基本原理,1.基本思想,主成分分析是一種數(shù)學(xué)降維的方法，找出幾個(gè)綜合變量來(lái)代替原來(lái)眾多的變量，使這些綜合變量能盡可能地代表原來(lái)變量的信息量，而且彼此之間互不相關(guān)。這種將把多個(gè)變量化為少數(shù)幾個(gè)互相無(wú)關(guān)的綜合變量的統(tǒng)計(jì)分析方法就叫做主成分分析或主分量分析。,2.主成分的數(shù)

2、學(xué)模型,設(shè)X1,X2,Xp，為實(shí)際問(wèn)題的p個(gè)n維隨機(jī)變量(p項(xiàng)指標(biāo))記X=(X1,X2,Xp)T，其協(xié)方差矩陣為,它是一個(gè)p階的非負(fù)定矩陣。設(shè)變量x1,x2,xp經(jīng)過(guò)線性變換后得到新的綜合變量Y1,Y2,Yp ，即,或,(5.1.1),其中系數(shù) 為常數(shù)向量。要求(5.1.1)滿足以下條件：,（1）系數(shù)向量是單位向量，即,（2）不同的主成分不相關(guān)，即,（3）各主成分的方差遞減，即,(5.1.2),(5.1.3),(5.1.4),于是，稱Y1為第一主成分，Y2為第二主成分，依此類推，Yp稱為第p個(gè)主成分。主成分又叫主分量。這里lij我們稱為主成分的系數(shù)。,3. 主成分的求法及性質(zhì),當(dāng)總體X=(X1

3、,X2,Xp)T的協(xié)方差矩陣=(ij)p已知時(shí)，我們可根據(jù)下面的定理求出主成分。,定理5.1 設(shè)p維隨機(jī)向量X的協(xié)方差矩陣的特征值滿足12 p 0，相應(yīng)的單位正交特征向量為e1,e2,ep，則X的第i個(gè)主成分為,(5.1.5),其中 ，且,(5.1.6),證明： 令 , 則P為正交矩陣，且,若 為X的第一主成分，其中 ，令,則 ，且,只有當(dāng)h1=(1,0,0)(標(biāo)準(zhǔn)單位向量)時(shí)等號(hào)成立，這時(shí),因此，X的第1個(gè)主成分為:,且方差 Var(Y1)=1, 達(dá)到最大 .,若 為X的第二主成分，其中 ，且,則 ，且,從而,只有當(dāng)h2=(0,1,0)=2時(shí)等號(hào)成立，這時(shí),因此的第2個(gè)主成分為：,且方差 V

4、ar(Y2)=2, 達(dá)到最大 .,類似可得其余主成分的表達(dá)式，且各主成份的方差等于相應(yīng)的特征值. 定理5.1表明：求X的主成分等價(jià)于求它的協(xié)方差矩陣的所有特征值及相應(yīng)的正交單位化特征向量.,推論：若記Y=(Y1,Y2,Yp)T為主成分向量，矩陣p=(e1,e2,ep),則 Y=pTX，且Y的協(xié)方差,主成分的總方差,證明: 由(5.1.5)式，顯然有Y=PTX，又由(5.1.6)式，有,又因?yàn)?此性質(zhì)表明主成分分析是將p個(gè)原始變量的總方差分解為p個(gè)不相關(guān)變量Y1,Y2,Yp的方差之和. 由于Var(Yk)=k，因此 描述了第k個(gè)主成分提取的信息占總信息的份額.,我們稱 為第個(gè)主成分的貢獻(xiàn)率，他表

5、示第個(gè)主成分提取的信息占總信息的百分比.,稱前m個(gè)主成分的貢獻(xiàn)率之和,為累計(jì)貢獻(xiàn)率，它表示前m個(gè)主成分綜合提供總信息的程度.通常mp，且累計(jì)貢獻(xiàn)率達(dá)到80%以上.,在實(shí)際應(yīng)用中，選擇了重要的主成分后，還要注意主成分實(shí)際含義解釋。主成分分析中一個(gè)很關(guān)鍵的問(wèn)題是如何給主成分賦予新的意義，給出合理的解釋。一般而言，這個(gè)解釋是根據(jù)主成分表達(dá)式的系數(shù)結(jié)合定性分析來(lái)進(jìn)行的。,在MATLAB中，運(yùn)用協(xié)方差矩陣進(jìn)行主成分分析的命令pcacov，調(diào)用格式： PC=pcacov(X) PC,latent,explained=pcacov(X),我們稱Yi與Xj的相關(guān)系數(shù)為因子載荷量，由于因子載荷量與主成分的系數(shù)向

6、量成正比，與標(biāo)準(zhǔn)差成反比，因此因子載荷量的絕對(duì)值大小刻畫(huà)了該主成分的成因，可以解釋第j個(gè)變量對(duì)第i個(gè)主成分的重要程度.,其中輸入X是定理5.1中的協(xié)方差矩陣，輸出PC為定理5.1中的矩陣P,latent是協(xié)方差矩陣的從大到小排列的特征值向量,explained表示貢獻(xiàn)率向量即每個(gè)每個(gè)主成份的方差在觀測(cè)量總方差中所占的百分?jǐn)?shù)向量。,例5.1.1設(shè)隨機(jī)向量X=(X1,X2,X3)T的協(xié)方差矩陣為,求X的各主成份以及各主成份的貢獻(xiàn)率。,clear S=2,2,-2;2,5,-4;-2,-4,5;%S表示總體的協(xié)方差矩陣PC,vary,explained=pcacov(S)%總體主成分分析,解：因?yàn)橐?/p>

7、知總體的協(xié)方差矩陣，所以可調(diào)用主成分分析的命令pcacov，程序如下：,程序輸出結(jié)果： PC =-0.3333 0 0.9428 -0.6667 0.7071 -0.2357 0.6667 0.7071 0.2357 %主成份變換矩陣 vary =10.0000 1.0000 1.0000 %主成份方差向量 explained =83.3333 8.3333 8.3333 %各主成份貢獻(xiàn)率向量,由程序輸出結(jié)果以及公式(5.1.5)知，X的主成份為：,主成份的方差,第一主成份的貢獻(xiàn)率,前倆個(gè)主成份的累計(jì)貢獻(xiàn)率為 83.3333%+8.3333%=91.6666%， 因此，若用前兩個(gè)主成分代替原來(lái)

8、三個(gè)變量，其信息損失為8.3333%，是很小的.,定理5.2 設(shè)Y=(Y1,Y2,Yp)T為總體X=(X1,X2,Xp)T的主成分向量，則主成分Yi與變量Xj的相關(guān)系數(shù),(5.1.8),證明：由定理5.1及其推論，因?yàn)閅=PTX,所以X=PY從而,于是,所以Yi與Xj的相關(guān)系數(shù)為,Yi與Xj相關(guān)系數(shù)反映了主成分與原變量的關(guān)聯(lián)程度，它與Xj標(biāo)準(zhǔn)差成反比，與主成分的標(biāo)準(zhǔn)差成正比。,若記 ，則由代數(shù)學(xué)可以證明： (5.1.9) 其中diag()表示協(xié)方差矩陣的主對(duì)角線元組成的對(duì)角矩陣，P是主成份矩陣，是特征值對(duì)角矩陣。,例5.1.2 （續(xù)例5.1.1 ）求主成份與原變量的相關(guān)系數(shù).,解：由例5.1.

9、1的條件與結(jié)果，根據(jù) (5.1.9)式，可編寫(xiě)MATLAB程序如下：,clear S=2,2,-2;2,5,-4;-2,-4,5; %S表示總體的協(xié)方差矩陣PC,vary,explained=pcacov(S); %總體主成分分析 S1=diag(diag(S); %協(xié)方差矩陣的主對(duì)角線元組成的對(duì)角矩陣 SYX=inv(sqrt(S1)*PC*sqrt(diag(vary); %按(5.1.9)式計(jì)算 SYX,程序輸出結(jié)果： SYX = -0.7454 0 0.6667 -0.9428 0.3162 -0.1054 0.9428 0.3162 0.1054,所以，SYX的第一列元素依次為Y1與

10、X1,X2,X3的相關(guān)系數(shù)，即其余各列的元素類推。結(jié)果表明Y1與X2,X3高度相關(guān)，Y2與X1不相關(guān)等。,4. 標(biāo)準(zhǔn)化變量的主成分,實(shí)際問(wèn)題經(jīng)常遇到不同的指標(biāo)具有不同的量綱，有時(shí)會(huì)導(dǎo)致各指標(biāo)取值的分散程度較大，這樣在計(jì)算協(xié)方差矩陣時(shí)，可能出現(xiàn)總體的方差主要受方差較大的數(shù)據(jù)的控制，可能造成不合理的結(jié)果.為消除量綱的影響，我們對(duì)原始數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化.即令,(5.1.10),其中i=EXi，X*的協(xié)方差矩陣就是原始數(shù)據(jù)的相關(guān)系數(shù)矩陣 .計(jì)算標(biāo)準(zhǔn)化變量的主成分公式為：,其中X*是標(biāo)準(zhǔn)化以后的數(shù)據(jù)，(ei*)T是相關(guān)系數(shù)矩陣的特征值對(duì)應(yīng)的特征向量.,標(biāo)準(zhǔn)化變量的主成分具有以下性質(zhì) 性質(zhì)1 總體方差和等于向

11、量的維數(shù),其中i*是相關(guān)系數(shù)矩陣的特征值.,(5.1.11),性質(zhì)2 標(biāo)準(zhǔn)化變量的第個(gè)主成分的貢獻(xiàn)率為,標(biāo)準(zhǔn)化變量的前個(gè)主成分的累積貢獻(xiàn)率為,性質(zhì)3 主成分Yi*與標(biāo)準(zhǔn)化數(shù)據(jù)Xj*的相關(guān)系數(shù)為,值得注意的是同一個(gè)總體，分別從協(xié)方差矩陣和相關(guān)系數(shù)矩陣出發(fā)進(jìn)行主成份分析，所得的主成份的貢獻(xiàn)率可以不同，讀者可參考本章習(xí)題的第1題。,5.1.2樣本主成分分析,實(shí)際問(wèn)題中，總體的協(xié)方差矩陣一般是未知的，具有的資料只是來(lái)自于X的一個(gè)容量為n的樣本觀測(cè)數(shù)據(jù).設(shè),為取自總體的一個(gè)容量為n的簡(jiǎn)單隨機(jī)樣本，由第二章知樣本協(xié)方差矩陣及樣本相關(guān)矩陣分別為,其中,分別以S和R作為總體和的估計(jì)，然后按總體主成分分析的方法

12、作樣本主成分分析。關(guān)于樣本主成份, 有如下結(jié)論 :,設(shè)Spp為樣本協(xié)方差矩陣, 其特征值 , 相應(yīng)的單位正交化特征向量 , 第k個(gè)樣本主成份表為,當(dāng)依次代入觀測(cè)值xi時(shí), 便得到第k個(gè)樣本主成分的n個(gè)觀測(cè)值,稱其為第k個(gè)樣本主成份的得分.,同總體主成份分析一樣，為了消除量綱的影響，可利用第二章介紹的方法可以對(duì)樣本進(jìn)行標(biāo)準(zhǔn)化，即令,(5.1.12),由于標(biāo)準(zhǔn)化的樣本數(shù)據(jù)的協(xié)方差矩陣也即原始數(shù)據(jù)的樣本相關(guān)系數(shù)矩陣，由樣本相關(guān)系數(shù)矩陣出發(fā)作主成份分析即可。,在MATLAB中，運(yùn)用樣本數(shù)據(jù)矩陣進(jìn)行主成分分析的命令為princomp，調(diào)用格式： PC=princomp(X) PC,SCORE,laten

13、t,tsquare=princomp(X),其中輸入X是樣本數(shù)據(jù)矩陣，輸出PC為主成份矩陣,SCORE是樣本主成份的得分，latent是X的方差矩陣的特征值向量, tsquare是每個(gè)數(shù)據(jù)點(diǎn)的HotellingT2統(tǒng)計(jì)量。,綜上，由樣本觀測(cè)數(shù)據(jù)矩陣進(jìn)行主成分分析的步驟為： 第一步，對(duì)原始數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理； 第二步：計(jì)算樣本相關(guān)系數(shù)矩陣； 第三步：求相關(guān)系數(shù)矩陣的特征值和相應(yīng)的特征向量； 第四步：選擇重要的主成分，并寫(xiě)出主成分表達(dá)式。 第五步：計(jì)算主成分得分； 第六步：依據(jù)主成分得分的數(shù)據(jù)，進(jìn)一步從事統(tǒng)計(jì)分析.,例5.1.3 對(duì)10名男中學(xué)生的身高,胸圍和體重進(jìn)行測(cè)量，得數(shù)據(jù)如表5.1所示，對(duì)

14、其做主成分分析.,表5.1 10名男中學(xué)生的身離、胸圍及體重?cái)?shù)據(jù),解：表5.2中的數(shù)據(jù)可認(rèn)作為來(lái)自總體的樣本，由樣本主成份分析法，Matlab程序如下：,X=149.569.538.5;162.57755.5;162.778.550.8;162.287.565.5;156.574.549.0;156.174.545.5;172.076.551.0;173.281.559.5;159.574.543.5;157.77953.5; %輸入樣本數(shù)據(jù) P,SCORE,latent=princomp(X) %主成份分析,P =0.5592 0.8277 -0.0480 0.4213 -0.3335 -0

15、.8434 0.7140 -0.4514 0.5352 %正交單位化特征向量,輸出,latent=110.0041 25.3245 1.5680 %特征值,所以，樣本各主成分的貢獻(xiàn)率分別為,因?yàn)榍皟蓚€(gè)主成分的累計(jì)貢獻(xiàn)率已達(dá)88. 860%，實(shí)際應(yīng)用中可只取前兩個(gè)主成分，即,第一主成分是身高值、胸圍值和體重值的加權(quán)和，當(dāng)一個(gè)學(xué)生比較魁梧時(shí)，所對(duì)應(yīng)的值也較大，故第一主成分是反映學(xué)生身材是否魁梧的綜合指標(biāo)，稱之為“大小”因子.而在第二主成分的表達(dá)式中，身高前的系數(shù)為正，而胸圍和體重前的系數(shù)為負(fù)，瘦高型學(xué)生的值會(huì)較大，故是反映學(xué)生體型特征的綜合指標(biāo)，稱之為“形狀”因子。,例5.1.4 影響我國(guó)糧食安全

16、生產(chǎn)的主要因素有以下幾個(gè)方面: 有效灌溉面積，糧食播種面積，成災(zāi)面積，財(cái)政投入，農(nóng)業(yè)勞動(dòng)力 ，農(nóng)村用電量，農(nóng)業(yè)機(jī)械總動(dòng)力 及農(nóng)業(yè)化肥施用量，具體數(shù)據(jù)如表5.2所示。對(duì)糧食安全生產(chǎn)因素作主成分分析。,解：由于各個(gè)指標(biāo)的單位不同，且各指標(biāo)的方差相差很大，所以首先對(duì)樣本數(shù)據(jù)進(jìn)行無(wú)量綱的變換，然后對(duì)標(biāo)準(zhǔn)化的樣本數(shù)據(jù)進(jìn)行主成份分析。程序如下：,%將表5.3中各指標(biāo)值作為矩陣X輸入，省略號(hào)表示數(shù)據(jù)省略了 X=4740.31, 11346.60, 178.20, 221.76, 33336.00, 844.50, 28707.70, 647.58;5847.17, 10679.30, 222.80, 454

17、4.01, 20078.6, 5713.20, 82190.41,1309.19;,表5.2 影響我國(guó)糧食安全生產(chǎn)的主要因素表,X1=zscore(X); %按公式(5.1.12)對(duì)樣本數(shù)據(jù)標(biāo)準(zhǔn)化 pc,la,tent=princomp(X1); %主成份分析，pc是特征向量矩陣，la得分矩陣，tent特征值 tents=sum(tent); %特征值總和 gxl= tent/ tents; %各個(gè)主成份貢獻(xiàn)率,pc = %正交單位化特征向量 0.3933 -0.1518 -0.0544 0.3944 0.5494 0.3701 -0.1396 0.4533 -0.2821 0.3173 -0

18、.7669 0.4403 0.0325 -0.1907 -0.0141 -0.0121 -0.0479 -0.8701 -0.4379 -0.1998 -0.0337 -0.0838 0.0257 -0.0115 0.3856 0.1952 -0.2572 -0.3364 0.2998 0.0210 0.7181 -0.1668 -0.3605 -0.2487 0.3850 0.4799 0.2640 -0.3102 0.5109 -0.0512 0.4083 0.0358 0.0265 -0.0014 -0.2123 -0.7233 0.0088 0.5128 0.4070 -0.0643

19、0.0422 0.1814 0.2707 -0.3433 -0.3636 -0.6872 0.3905 -0.1175 -0.0151 0.4841 -0.6466 0.2861 0.2659 -0.1692,由于第一、二主成分的累計(jì)貢獻(xiàn)率為88.04%?？扇∏皟蓚€(gè)主成分，即,gxl = 0.7388 0.1416 0.0765 0.0355 0.0043 0.0027 0.0005 0.0001 %樣本各主成分的貢獻(xiàn)率,其中,e1T=( 0.3933,-0.2821,-0.0479,0.3856,-0.3605,0.4083,0.4070,0.3905),e2T=(-0.1518,0.317

20、3,-0.8701,0.1952,-0.2487,0.0358,-0.0643, -0.1175).,第一主成分反映了我國(guó)在農(nóng)業(yè)生產(chǎn)基礎(chǔ)設(shè)施投入情況，第二主成分反映了我國(guó)糧食生產(chǎn)抗災(zāi)能力與勞動(dòng)力情況。,5.2 主成分分析的應(yīng)用,主成分分析的應(yīng)用范圍非常廣泛，諸如投資組合風(fēng)險(xiǎn)管理，企業(yè)效益的綜合評(píng)價(jià)，圖像特征識(shí)別，機(jī)械加工或傳感器故障檢測(cè)，災(zāi)害損失分析等.將主成分分析與聚類分析、判別分析以及回歸分析方法相結(jié)合，還可以解決更多實(shí)際問(wèn)題.,5.2.1 主成分分析用于綜合評(píng)價(jià),主成分分析用于綜合評(píng)價(jià)的一般步驟：, 若各指標(biāo)的屬性不同(成本型、利潤(rùn)型、適度型)，則將原始數(shù)據(jù)矩陣A統(tǒng)一趨勢(shì)化，得到屬性一致

21、的指標(biāo)矩陣B（具體過(guò)程參見(jiàn)第二章的數(shù)據(jù)變換一節(jié)）., 計(jì)算的協(xié)方差矩陣,或相關(guān)系數(shù)矩陣R(當(dāng)?shù)牧烤V不同，或矩陣主對(duì)角元素差距過(guò)大時(shí)，用相關(guān)系數(shù)矩陣)., 計(jì)算或R的特征值與相應(yīng)的特征向量., 根據(jù)特征值計(jì)算累計(jì)貢獻(xiàn)率，確定主成分的個(gè)數(shù)，而特征向量就是主成分的系數(shù)向量。, 計(jì)算主成分的數(shù)值（即主成分得分）。若利用協(xié)方差矩陣計(jì)算特征值與特征向量，則主成分得分為,若利用相關(guān)系數(shù)矩陣R計(jì)算特征值與特征向量，則主成分得分為：,其中，V是特征向量矩陣，B*是將矩陣標(biāo)準(zhǔn)化以后的矩陣（即zscore(B)）, 計(jì)算綜合評(píng)價(jià)值，進(jìn)行排序.若為效益型矩陣，則評(píng)價(jià)值越大排名越靠前；若為成本型矩陣，則評(píng)價(jià)值越小排名越

22、靠前.,通常計(jì)算綜合評(píng)價(jià)值的公式為:Z=FW 其中F是主成分得分矩陣，W是將特征值歸一化后得到的權(quán)向量.,例5.2.1 根據(jù)2008年安徽統(tǒng)計(jì)年鑒資料，選擇x1：工業(yè)總產(chǎn)值(現(xiàn)價(jià))，X2：工業(yè)銷售產(chǎn)值(當(dāng)年價(jià))，X3：流動(dòng)資產(chǎn)年平均余額，X4：固定資產(chǎn)凈值年平均余額， X5 ：業(yè)務(wù)收入,X6：利潤(rùn)總額等六項(xiàng)指標(biāo)進(jìn)行主成分分析. (1)選取指標(biāo)是否合適？（2）給出各市大中型工業(yè)企業(yè)排名。,解：首先輸入數(shù)據(jù) A=data; % data即表5.3中數(shù)據(jù) R=corrcoef(A); 得到的相關(guān)系數(shù)矩陣為：,，,表 5. 3安徽省各市大中型工業(yè)企業(yè)主要經(jīng)濟(jì)指標(biāo)(單位：億元),由于r12=r21=1，

23、表明指標(biāo)x1,x2完全線性相關(guān)，故只需保留一個(gè)指標(biāo)x2.,A1=A(:,2:6)./ones(17,1)*std(A(:,2:6); % 消除量綱 v,d=eig(corrcoef(A1); % 計(jì)算特征值與特征向量 w=sum(d)/sum(sum(d); % 計(jì)算貢獻(xiàn)率 F=A1-ones(17,1)*mean(A1)*d(:,5);% 計(jì)算主成分得分 F1,I1=sort(F,descend); % I1各名次的序號(hào) F2,I2=sort(I1); % I2給出各市排名,表5.4 特征值、特征向量及貢獻(xiàn)率,結(jié)果見(jiàn)表5.4、5.5,表5.5 各市第一主成分得分排名,5.2.3 主成分分析用

24、于分類,利用主成分分析可以計(jì)算出主成分的得分，如果主成分的貢獻(xiàn)率較大，則其提取了原始數(shù)據(jù)的主要信息，因此可以利用主成分分析進(jìn)行分類。,例5.2.2做出例4.1.1中蠓蟲(chóng)原始數(shù)據(jù)圖與主成分得分?jǐn)?shù)據(jù)圖。,解：作圖程序如下,apf=1.14,1.78;1.18,1.96;1.20,1.86;1.26,2.;1.28,2;1.30,1.96; af=1.24,1.72;1.36,1.74;1.38,1.64;1.38,1.82;1.38,1.90;1.40,1.70;1.48,1.82;1.54,1.82;1.56,2.08; x=1.24,1.8;1.28,1.84; 1.4,2.04; % 輸入原

25、始數(shù)據(jù) subplot(2,1,1) plot(apf(:,1),apf(:,2),*,af(:,1),af(:,2),or,x(:,1),x(:,2),p) % 原始數(shù)據(jù)圖形 c1,s1,l1,t1=princomp(apf;af;x); % 計(jì)算主成分得分s1 subplot(2,1,2), plot(1:6,s1(1:6,2),*) hold on plot(7:15,s1(7:15,2),or) hold on plot(16:18,s1(16:18,2),p) legend(apf,af,x) % 主成分得分圖形 hold on plot(0:18,0*ones(1,19),-) %

26、畫(huà)直線,圖5.1蠓蟲(chóng)原始數(shù)據(jù)圖(上方)與主成分得分?jǐn)?shù)據(jù)圖(下方),從圖5.1上方圖形可以看出原始數(shù)據(jù)的圖形中不同蠓蟲(chóng)觸長(zhǎng)與翅長(zhǎng)兩個(gè)指標(biāo)無(wú)法做到楚漢分明，但是下方主成分得分的圖形則顯示Af類蠓蟲(chóng)的第一個(gè)主成分小于Apf的得分,且基本小于零；未知的三個(gè)蠓蟲(chóng)第一主成分的得分大于Af而小于Apf.依此為判據(jù)可以得到與第4章馬氏距離判別同樣的結(jié)果.,由于例4.1.1蠓蟲(chóng)的原始數(shù)據(jù)是二維向量可以做出圖形，如果原始數(shù)據(jù)維數(shù)大于3，借助主成分分析，我們可以選取23個(gè)主成分，做出主成分得分的散點(diǎn)圖.,例5.2.3 瑞士銀行紙幣數(shù)據(jù)(見(jiàn)附錄A：Swiss Bank Notes)為2006的矩陣，其中100行是真紙

27、幣數(shù)據(jù)，100行是假幣數(shù)據(jù).六項(xiàng)指標(biāo)為：紙幣長(zhǎng)度，左、右側(cè)紙幣高度，上、下圖廓內(nèi)骨架距離以及對(duì)角線長(zhǎng)度，分別用表示.(1)利用協(xié)方差矩陣進(jìn)行主成分分析，此時(shí)可否利用第一主成分得分進(jìn)行排名？(2)利用R矩陣進(jìn)行主成分分析，此時(shí)可否利用第一主成分得分進(jìn)行排名？(3)選擇兩個(gè)主成分的得分做出平面圖形，能否從圖形上分別真假紙幣？,解：(1)首先輸入原始數(shù)據(jù)，然后利用協(xié)方差矩陣進(jìn)行主成分分析，程序如下：,a=data; %輸入原始數(shù)據(jù) M,N=size(a); v,d=eig(cov(a); %樣本協(xié)方差矩陣的特征值 輸出結(jié)果顯示最大特征值對(duì)應(yīng)的不是正向量，所以不能用第一主成分得分進(jìn)行排名. (2)利用

28、R矩陣進(jìn)行主成分分析，程序如下： a=data; %輸入原始數(shù)據(jù) M,N=size(a); %計(jì)算原始數(shù)據(jù)維數(shù) for i=1:N for j=1:N R(i,j)=2*dot(a(:,i),a(:,j)./sum(a(:,i).2)+ sum(a(:,j).2); %計(jì)算R矩陣 end end v,d=eig(R); %R矩陣的特征值與特征向量 q=sum(d)/sum(sum(d); %計(jì)算貢獻(xiàn)率,輸出結(jié)果顯示最大特征值對(duì)應(yīng)的是正向量，且其貢獻(xiàn)率達(dá)到60，所以可以用第一主成分得分進(jìn)行排名.,(3)分別選擇第一，二兩個(gè)主成分，第二，三主成分以及第一，三兩個(gè)主成分得分做出平面圖形，程序如下:,

29、a=data; %輸入原始數(shù)據(jù) M,N=size(a); v,d=eig(cov(a); %樣本協(xié)方差矩陣的特征值 q=sum(d)/sum(sum(d); %計(jì)算貢獻(xiàn)率 F=a*v; %計(jì)算主成分得分 subplot(2,2,1) plot(F(1:100,6),F(1:100,5),o,F(101:200,6),F(101:200,5),+) %第1、2兩個(gè)主成分 title(w-pc1-pc2) subplot(2,2,2) plot(F(1:100,5),F(1:100,4),or,F(101:200,5),F(101:200,4),+) %第2、3兩個(gè)主成分 title(w-pc2-

30、pc3) subplot(2,2,3) plot(F(1:100,6),F(1:100,4),or,F(101:200,6),F(101:200,4),+) %第1、3兩個(gè)主成分 title(w-pc1-pc3) subplot(2,2,4) plot(1:6,fliplr(sum(d),-or) %從大到小特征值 legend(eigenvalues ) title(從大到小特征值),從圖5.2可以看出，第一，第二兩個(gè)主成分以及第一，第三兩個(gè)主成分的得分做出平面圖形還可以區(qū)分真假紙幣，而第二，第三主成分的得分做出的圖形（右上）則呈現(xiàn)出混亂現(xiàn)象.為了使得從圖形上更好地區(qū)別真假紙幣，我們可以計(jì)算

31、加權(quán)主成分得分，將上面程序中計(jì)算主成分得分的命令行“F=a*v”改寫(xiě)成如下,圖 5.2瑞士銀行紙幣主成分圖形,w=q/sum(q); %計(jì)算權(quán)向量 F=a.*(ones(M,1)*w)*v; %計(jì)算加權(quán)主成分得分,繪圖命令如“title (pc1-pc2) 改寫(xiě)成title(w-pc1-pc2)，加了字符“w-”,此時(shí)輸出圖形如圖5.3所示，其效果明顯優(yōu)于圖5.2.,圖 5.3 瑞士銀行紙幣加權(quán)主成分圖形,5.2.3 主成分分析用于信號(hào)分離,在電波、聲波訊號(hào)傳輸過(guò)程中由于各種干擾較多,給實(shí)際數(shù)據(jù)帶來(lái)了一定的噪聲污染，主成分分析利用真實(shí)信號(hào)在混合信號(hào)中占主導(dǎo)地位的特性，可以進(jìn)行信號(hào)分離與提取。,

32、例5.2.4 設(shè) (1)構(gòu)造信號(hào) 并做出圖像； (2)利用MZ算法(Marsaglias Ziggurat algorithm)生成隨機(jī)噪聲，將三個(gè)信號(hào)混合為：,(3)使用主成分分析與獨(dú)立成分分析將混合信號(hào)分離.,解：(1) 程序如下：,clear i = 1:0.01:10 * pi; % 輸入t的取值 dummy index = sort(sin(i); s1(index,1) = i/31; s1 = s1 - mean(s1); % 生成s1 s2 = abs(cos(1.89*i); s2 = s2 - mean(s2); % 生成s2 s3 = sin(3.43*i); % 生成s

33、3 subplot(311), plot(s1), ylabel(s_1), title(Raw signals) subplot(312), plot(s2), ylabel(s_2) subplot(313), plot(s3), ylabel(s_3),圖 5.4 三個(gè)原始信號(hào)圖形,(2) 生成三個(gè)混合信號(hào)程序如下：,randn(state,1); %生成隨機(jī)噪聲 y1=0.1*s1 + 0.8*s2 + 0.01*randn(length(i),1); %生成混合信號(hào)y1 y2= 0.4*s1 + 0.3*s2 + 0.01*randn(length(i),1); %生成混合信號(hào)y2

34、y3=0.1*s1 + s3 + 0.02*randn(length(i),1); %生成混合信號(hào)y3 y=y1,y2,y3; subplot(311), plot(y(:,1), ylabel(y_1), title(Mixed signals) subplot(312), plot(y(:,2), ylabel(y_2) subplot(313), plot(y(:,3), ylabel(y_3),圖 5.5 混合信號(hào)圖形,(3) 分離信號(hào)程序,coeff, score,latent = princomp(y) sPCA=score; sPCA = sPCA./repmat(std(sPC

35、A),length(sPCA),1); subplot(321), plot(sPCA(:,1) ylabel(s_PCA1), title(Separated signals - PCA) subplot(323), plot(sPCA(:,2), ylabel(s_PCA2) subplot(325), plot(sPCA(:,3), ylabel(s_PCA3),rand(state,1); div = 0;B = orth(rand(3, 3)-0.5);BOld = zeros(size(B); while (1 - div) eps B = B * real(inv(B * B)(

36、1/2);div = min(abs(diag(B * BOld); BOld = B; B = (sPCA*( sPCA * B). 3)/length(sPCA) -3*B; sICA = sPCA*B; end,% 獨(dú)立成分分析,subplot(322), plot(sICA(:,1), ylabel(s_ICA1), title(Separated signals - ICA) subplot(324), plot(sICA(:,2), ylabel(s_ICA2) subplot(326), plot(sICA(:,3), ylabel(s_ICA3),圖 5.6 分離信號(hào)圖形,5.

37、3 典型相關(guān)分析,在對(duì)經(jīng)濟(jì)和管理問(wèn)題的研究中，不僅經(jīng)常需要考察兩個(gè)變量之間的相關(guān)程度，而且還經(jīng)常需要考察多個(gè)變量與多個(gè)變量之間即兩組變量之間的相關(guān)性。比如工廠管理人員需要了解原料的主要質(zhì)量指標(biāo)x1,x2,xp 與產(chǎn)品的主要質(zhì)量指標(biāo)Y1,Y2,Yq 之間的相關(guān)性，以便提高產(chǎn)品質(zhì)量；醫(yī)生要根據(jù)病人的一組體檢化驗(yàn)指標(biāo)與一些疾病之間的相關(guān)性，以便確定治療方法等等.典型相關(guān)分析就是測(cè)度兩組變量之間相關(guān)程度的一種多元統(tǒng)計(jì)方法，它是兩個(gè)隨機(jī)變量之間的相關(guān)性在兩組變量之下的推廣,5.3.1典型相關(guān)分析的基本原理,對(duì)于兩組隨機(jī)變量(X1,X2,Xp)和(Y1,Y2,Yq) ，象主成分分析那樣，考慮(X1,X2,

38、Xp)一個(gè)線性組合U及的(Y1,Y2,Yq)一個(gè)線性組合V，希望找到的U和V之間有最大可能的相關(guān)系數(shù)，以充分反映兩組變量間的關(guān)系。這樣就把研究?jī)山M隨機(jī)變量間相關(guān)關(guān)系的問(wèn)題轉(zhuǎn)化為研究?jī)蓚€(gè)隨機(jī)變量間的相關(guān)關(guān)系。如果一對(duì)變量(U,V)還不能完全刻劃兩組變量間的相關(guān)關(guān)系時(shí)，可以繼續(xù)找第二對(duì)變量，希望這對(duì)變量在與第一對(duì)變量(U,V)不相關(guān)的情況下也具有盡可能大的相關(guān)系數(shù)。直到進(jìn)行到找不到相關(guān)變量對(duì)時(shí)為止。這便引導(dǎo)出典型相關(guān)變量的概念。,1.總體典型相關(guān)變量,設(shè)有兩組隨機(jī)變量,(XT,YT)T=(X1,X2,Xp,Y1,Y2,Yq)T的協(xié)方差矩陣為,其中，11=cov(X)， 22=cov(Y), 12=

39、 T21=cov(X,Y),根據(jù)典型相關(guān)思想是要尋找 的線性組合,(pq),將兩組合并成一組向量,（5.3.1）,使U1,V1的相關(guān)系數(shù)(U1,V1)達(dá)到最大，這里,由（5.3.1）式，,所以U1,V1的相關(guān)系數(shù)為,（5.3.2）,又由于相關(guān)系數(shù)與量綱無(wú)關(guān)，因此可設(shè)約束條件,（5.3.3）,滿足約束條件（5.3.3）的相關(guān)系數(shù)的最大值稱為第一典型相關(guān)系數(shù)，U1,V1稱為第一對(duì)典型相關(guān)變量.,典型相關(guān)分析在約束條件,a1T11a1=b1T22b1=1下，求,a1,b1，使得u1,v1=a1T12b1取得最大值.,如果(U1,V1)還不足以反映X,Y之間的相關(guān)性，還可構(gòu)造第二對(duì)線性組合：,使得(U

40、1,V1)與(U2,V2)不相關(guān)，即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=a2T12b2取得最大值.,一般地，若前k -1對(duì)典型變量還不足以反映X,Y之間的相關(guān)性，還可構(gòu)造第k對(duì)線性組合：,在約束條件 Var(uk)=Var(vk)=1,及cov(uk,uj)=cov(uk,vj)=cov(vk,uj)=cov(vk,vj)=0,(1jk) 求ak,bk，使得uk,vk=akT12bk取得最大值.,如此確定的(uk,vk)稱為X,Y的第

41、k對(duì)典型變量，相應(yīng)的uk,vk稱為第k個(gè)典型相關(guān)系數(shù).,2. 總體典型變量與典型相關(guān)系數(shù)的計(jì)算,(1) 計(jì)算矩陣(XT,YT)T的協(xié)方差矩陣或相關(guān)系數(shù)矩陣,(2) 令 或,求A,B的特征值12, 22, p2與對(duì)應(yīng)的正交單位特征向量ek,fk,k=1,p,(4) X,Y的第k個(gè)典型相關(guān)系數(shù)為：k,k=1,2,p,X=data; %輸入?yún)f(xié)方差矩陣X p=c1; q=c2; %c1 ,c2分別表示X，Y向量的維數(shù) R11=X(1:p,1:p); R12=X(1:p,p+1:p+q); %讀取11,12 R21=X(p+1:p+q,1:p); R22=X(p+1:p+q,p+1:p+q);%讀取21

42、,22 v1,d1=eig(R11); %計(jì)算R11的特征值與單位正交向量 v2,d2=eig(R22); %計(jì)算R22的特征值與單位正交向量 p1=inv(v1*sqrt(d1)*v1); p2=inv(v2*sqrt(d2)*v2); % p1，p2表示的平方根矩陣的逆 A=p1*R12*inv(R22)*R21*p1; %計(jì)算矩陣A B=p2*R21*inv(R11)*R12*p2; %計(jì)算矩陣B va,da=eig(A ), %計(jì)算A的特征值與特征向量 vb,db=eig(B), %計(jì)算B的特征值與特征向量 A1=p1*va, %計(jì)算典型相關(guān)變量U的系數(shù) B1=p2*vb, %計(jì)算典型

43、相關(guān)變量V的系數(shù) r=sqrt(sum(da), %計(jì)算典型相關(guān)系數(shù),以上過(guò)程的MATLAB實(shí)現(xiàn)程序如下：,例5.3.1 設(shè)樣本的相關(guān)系數(shù)矩陣為,計(jì)算典型相關(guān)系數(shù)與典型相關(guān)變量.,解：由于給出了相關(guān)系數(shù)矩陣，計(jì)算程序如下,R11=1,0.505;0.505,1; R12=0.569,0.602;0.422,0.467; R21=0.569,0.422;0.602,0.467; R22=1,0.926;0.926,1; v1,d1=eig(R11); %計(jì)算R11的特征值與單位正交向量 v2,d2=eig(R22); %計(jì)算R22的特征值與單位正交向量,p1=inv(v1*sqrt(d1)*v1

44、); %p1就是平方根矩陣的逆 p2=inv(v2*sqrt(d2)*v2); %p2就是平方根矩陣的逆 T1=p1*R12*inv(R22)*R21*p1; %計(jì)算矩陣A T2=p2*R21*inv(R11)*R12*p2; %計(jì)算矩陣B va,da=eig(T1), vb,db=eig(T2), A1=p1*va, %計(jì)算典型相關(guān)變量U的系數(shù) B1=p2*vb, %計(jì)算典型相關(guān)變量V的系數(shù) r=sqrt(sum(da), %計(jì)算典型相關(guān)系數(shù),典型相關(guān)系數(shù)為:,典型變量為：,5.3.2樣本的典型變量與典型相關(guān)系數(shù),在實(shí)際問(wèn)題中(XT,YT)T的協(xié)方差矩陣(或相關(guān)系數(shù)矩陣R)一般是未知的，我們

45、所具有的資料通常是關(guān)于X和Y的n組觀測(cè)數(shù)據(jù):,同主成分分析一樣，利用觀測(cè)數(shù)據(jù)的樣本協(xié)方差矩陣或者相關(guān)系數(shù)矩陣,或,作為或的估計(jì)，其中,以S代替或R代替所求得的典型變量和典型相關(guān)系數(shù)分別稱為樣本典型變量和樣本典型相關(guān)系數(shù)。此時(shí)樣本典型變量和典型相關(guān)系數(shù)計(jì)算方法同總體典型變量和典型相關(guān)系數(shù)的計(jì)算方法一樣。,在MATLAB中，樣本典型相關(guān)分析的命令canoncorr，其調(diào)用格式為 A,B,r,U,V,stats = canoncorr(X,Y),其中輸入X表示第一組向量的觀測(cè)矩陣，Y表示第二組向量的觀測(cè)矩陣，輸出A,B是典型相關(guān)變量的系數(shù)矩陣；r表示典型相關(guān)系數(shù)；U,V表示典型相關(guān)變量的得分；輸出s

46、tats包括wilks 、chisq及F統(tǒng)計(jì)量以及相應(yīng)的概率。,例5.3.2 某康復(fù)俱樂(lè)部對(duì)20名中年人測(cè)量了三項(xiàng)生理指標(biāo)：體重(weight)、腰圍(waist)、脈搏(pulse )和三項(xiàng)訓(xùn)練指標(biāo)：引體向上(chins )、起坐次數(shù)(situps )、跳躍次數(shù)(jumps )。其數(shù)據(jù)列于表5.6。試分析這兩組變量間的相關(guān)性。,解：三項(xiàng)生理指標(biāo)作為第一組向量X，三項(xiàng)訓(xùn)練指標(biāo)作為第二組向量Y，表5.6中的數(shù)據(jù)作為樣本數(shù)據(jù)，調(diào)用典型相關(guān)分析命令.程序如下:,DATA=; %將表5.6中的數(shù)據(jù)輸入DATA X=DATA(:,1:3); %第一組向量觀測(cè)值 Y=DATA(:,4:6); %第二組向量

47、觀測(cè)值 A,B,r,U,V,stats = canoncorr(X,Y),表5.6 某康復(fù)俱樂(lè)部測(cè)量的生理指標(biāo)和訓(xùn)練指標(biāo),輸出結(jié)果為：,A = -0.0314 -0.0763 0.0077 0.4932 0.3687 -0.1580 -0.0082 -0.0321 -0.1457 B = -0.0661 -0.0710 0.2453 -0.0168 0.0020 -0.0198 0.0140 0.0207 0.0082 r = 0.7956 0.2006 0.0726,5.3.3典型相關(guān)系數(shù)的顯著性檢驗(yàn),設(shè)總體X,Y的各對(duì)典型相關(guān)系數(shù)為1 2 p 0首先提出檢驗(yàn)原假設(shè)與備擇假設(shè),若不能拒絕原假

48、設(shè)，則 1= 2 = = k=0，此時(shí)不能做典型相關(guān)分析；若拒絕H0(1),繼續(xù)如下檢驗(yàn),若不能拒絕H0(2),表明只有第一對(duì)典型變量顯著相關(guān)外，其余變量均不顯著，實(shí)際應(yīng)用只需考慮第一對(duì)典型變量；若拒絕H0(2),則需檢驗(yàn)3是否為零,以此類推，若假設(shè)k-1 =0被拒絕，則檢驗(yàn),1.檢驗(yàn)方法,若不能拒絕H0(k),則只需考慮前k-1對(duì)典型相關(guān)變量,否則繼續(xù)檢驗(yàn)，直至檢驗(yàn)p是否為零.,在總體服從p+q維正態(tài)分布條件下，可用如下的似然比統(tǒng)計(jì)量進(jìn)行檢驗(yàn),對(duì)于給定的，計(jì)算概率,若pk ，即認(rèn)為第k對(duì)典型變量顯著相關(guān).上述檢驗(yàn)依次對(duì)k=1,2,p進(jìn)行，若對(duì)某個(gè)k檢驗(yàn)概率首次大于，則檢驗(yàn)停止，即認(rèn)為只有前k

49、-1對(duì)典型變量顯著相關(guān).,2.典型相關(guān)分析檢驗(yàn)的Matlab實(shí)現(xiàn),設(shè) 是取自總體的觀測(cè)數(shù)據(jù)，利用 MATLAB軟件進(jìn)行典型相關(guān)分析的步驟如下：, 輸入數(shù)據(jù)并計(jì)算協(xié)方差矩陣或相關(guān)系數(shù)矩陣,a=X, Y; % 此前X,Y的數(shù)據(jù)應(yīng)該已經(jīng)輸入 n,m=size(a); R=cov (a);, 計(jì)算典型相關(guān)系數(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); xgxs=sqrt(d);, 計(jì)算典型相關(guān)向量,X=X./ones(n,1)*std(X); Y=Y./ones(n,1)*std(Y); A,B = canoncorr(X,Y); U=(X-ones(n,1)*mean(X)*A V=(Y-ones(n,1)*mean(Y)*B, 典型相關(guān)系數(shù)的顯著性檢驗(yàn),其中,檢驗(yàn)程序如下： D=1-d; f1=fliplr(D);