龔 靜 陳俊華

(1.湖南環(huán)境生物職業(yè)技術(shù)學(xué)院，湖南 衡陽 421005;2.四川省林業(yè)科學(xué)研究院，四川 成都 610081)

主成份分析(principal components analysis，PCA)是將多個(gè)觀測指標(biāo)因子變化為少數(shù)幾個(gè)新指標(biāo)的一種多元統(tǒng)計(jì)方法［1］，其分析結(jié)果可用于回歸分析、聚類分析、神經(jīng)網(wǎng)格分析等等。我國林業(yè)專家應(yīng)用主成份分析法在多個(gè)領(lǐng)域［2～6］作了研究。主成份分析的算法中，最麻煩的是矩陣特征根的計(jì)算，目前在一些統(tǒng)計(jì)軟件［7］中有集成模塊，當(dāng)然也可以通過編程來計(jì)算，但數(shù)據(jù)錄入、修改不方便。Microsoft Excel能很好地解決這些問題。

Microsoft Excel是美國微軟公司研制開發(fā)的用于個(gè)人財(cái)務(wù)分析與規(guī)劃、公司營運(yùn)管理與目標(biāo)設(shè)定和薪水管理等的一個(gè)出色的電子表格軟件，它不僅具有友好的界面、強(qiáng)大的數(shù)據(jù)計(jì)算，并且還具有統(tǒng)計(jì)分析功能。Excel進(jìn)行數(shù)據(jù)計(jì)算與統(tǒng)計(jì)分析的一個(gè)顯著特點(diǎn)是，它可以進(jìn)行公式編輯和插人函數(shù)，而且具有連環(huán)鏈鎖計(jì)算能力［8］。有關(guān)矩陣的特征根和特征向量的Excel算法，相關(guān)文獻(xiàn)中也有報(bào)道［9～11］，但對(duì)整個(gè)主成份分析的過程卻不見系統(tǒng)報(bào)道。下面，筆者詳細(xì)介紹如何在Excel 2003中建立數(shù)據(jù)進(jìn)行主成份分析，以與搞林業(yè)工作的同仁們分享。

1 原始數(shù)據(jù)的錄入及相關(guān)系數(shù)的計(jì)算

本文以參考文獻(xiàn)［7］第11章“因子分析法”的數(shù)據(jù)為例，建立矩陣來進(jìn)行主成份分析，以便與專業(yè)軟件(SPSS)進(jìn)行數(shù)據(jù)結(jié)果精度的比較。

在Excel2003中錄入數(shù)據(jù)進(jìn)行相關(guān)系數(shù)的計(jì)算，非常簡單。步驟如下:

(1)打開 Microsoft Excel 2003，選擇區(qū)域 A1:F9，輸入原始數(shù)據(jù)(圖1)。

圖1 原始數(shù)據(jù)的錄入

(2)相關(guān)系數(shù)矩陣的計(jì)算

單擊“工具(T)”下拉菜單中的“數(shù)據(jù)分析(D)…”，在彈出的“數(shù)據(jù)分析”窗口中，選擇“相關(guān)系數(shù)”，輸入要計(jì)算數(shù)據(jù)的區(qū)域等選項(xiàng)，單擊“確定”后，結(jié)果如圖2所示。

圖2 相關(guān)系數(shù)計(jì)算結(jié)果

2 Excel中矩陣的建立及命名

在Excel中，可以在一個(gè)單元格區(qū)域內(nèi)通過逐個(gè)輸入矩陣的各個(gè)元素來建立矩陣，還可以使用數(shù)組公式和數(shù)組常量更加方便地建立矩陣［12］。

(1)新建一個(gè)工作表，選擇單元格區(qū)域A1:F6;

(2)輸入公式 ={1，－0.557，－0.443，0.535，=0.614，0.107;……}，即(1)中計(jì)算出來的相關(guān)系數(shù)矩陣。(在Excel的數(shù)組公式中，將矩陣元素用大括號(hào){}括起來稱為數(shù)組常量，其中不同列的元素用逗號(hào)隔開，不同行的元素用分號(hào)隔開［12］);

(3)同時(shí)按下Ctrl+Alt+Enter鍵，完成數(shù)組的輸入，并在左上角框中將其命名為數(shù)組A，如圖3所示。

3 特征根的計(jì)算

(1)選擇單元格區(qū)域G1:L9，輸入6階單位矩陣并將其命名為I;

圖3 矩陣的建立及命名

(2)在單元格 C7 中輸入 MDETERM(A)［8］來計(jì)算矩陣A的行列式值，其結(jié)果為0.037248;

(3)預(yù)測矩陣A的特征值所在區(qū)間:

Gersgorin圓盤定理［13］認(rèn)為，設(shè)對(duì)稱方陣A=(aij)n×n，則A的每一個(gè)特征值必屬于下述某個(gè)圓盤之中:

理論上A的特征值包含在區(qū)間［0，6)內(nèi)。在區(qū)間［0，6)取若干個(gè)點(diǎn)作特征多項(xiàng)式f(λ)=det(A－λI)的“平滑散點(diǎn)圖”，步驟為:

(1)在工作表Sheet1的單元格A8中輸入0，按住Ctrl鍵，拖動(dòng)鼠標(biāo)直到單元格A68，在單元格B8中，輸入工式=A8/10，同樣拖動(dòng)鼠標(biāo)直到B68，即輸入0，0.1，0.2，…(每隔 0.1 取一個(gè)值)，在單元格C7中輸入公式=MDETERM(A－B8*I)，按Ctrl+Shift+Enter鍵，在C8中形成數(shù)組公式，則C8中的值0.037248為特征多項(xiàng)式f(λ)在λ=0的值;

(2)選定單元格C8，拖動(dòng)鼠標(biāo)直到C68，計(jì)算各特征根的值，可見有幾次符號(hào)改變，而每一次符號(hào)改變就說明此附近有一特征根，當(dāng)?shù)?.3以后符號(hào)不再改變，說明2.3附近為A的最大特征根;

(3)選定單元格B8:C31(即2.3)，作“XY散點(diǎn)圖”的“平滑線散點(diǎn)圖”類型曲線，見圖4。從圖4可以看出A的特征多項(xiàng)式與x軸有6個(gè)交點(diǎn)，還可以看出這6 個(gè)點(diǎn)分別落在區(qū)間［0，0.1］，(0.1，0.25)，(0.6，0.7)，(0.9，1)，(1.7，1.8)，(2.2，2.3)。

圖4 特征多項(xiàng)式f(λ)的平滑線散點(diǎn)圖

(4)求A的全部特征根

1)選擇工作表Sheet2里面的單元格A1:A6，即將這6個(gè)單元格作為求特征根的可變單元格;

2)在B1中輸入公式=MDETERM(A－A1*I)，按 Ctrl+Shift+Enter鍵，形成數(shù)組公式:{=MDETERM(A－A1*I)};

3)拖動(dòng)B1到B6;

4)單擊“工具”菜單中的“規(guī)劃求解”命令，打開“規(guī)劃求解參數(shù)”對(duì)話框(若沒有這項(xiàng)，可單擊“工具”菜單中的“加載宏(I)”，在彈出的窗口中將“規(guī)劃求解”勾選上);

5)在“設(shè)置目標(biāo)單元格(E)”中輸入B1，單擊“等于”欄中的“值為”輸入0，在“可變單元格(B)”中，輸入或選擇單元格區(qū)域A1:A6;

6)單擊“添加(A)”按鈕，在“添加約束”對(duì)話框中添加條件:A1＞=0，A1＜=0.1;

7)單擊“求解”按鈕完成操作。此時(shí)單元格A1中的值:0.067662就是實(shí)對(duì)稱矩陣A的第一個(gè)特征值;

8)重復(fù)步驟5)～7)(注意條件)，這樣就求出了全部特征根，如圖5所示，這與文獻(xiàn)［7］中用統(tǒng)計(jì)軟件SPSS求出來的值非常接近(SPSS求出的6個(gè)值分 別 為 2.285，1.793，0.954，0.700，0.201，0.068)［7］。

圖5 實(shí)對(duì)稱矩陣A的全部特征值

4 求累計(jì)方差貢獻(xiàn)率

方差貢獻(xiàn)率的計(jì)算公式較為簡單，例如，對(duì)于λ1=2.285，則方差貢獻(xiàn)率為=2.285/6=38.09%，以此類推。計(jì)算出來的方差貢獻(xiàn)率及累計(jì)方差貢獻(xiàn)率見圖5。

5 每個(gè)特征值對(duì)應(yīng)的特征向量及因子負(fù)荷矩陣的求解

(1)特征向量的解法

對(duì)應(yīng)于λ1的特征向量可由下面齊次方程組求出:

將矩陣A及λ的值代人方程組，令x6=1，用前面的5個(gè)方程聯(lián)立求解，就可以求得相對(duì)于λ的特征向量，這樣求特征向量問題就是求解多元一次線性方程組的問題。Excel解線性方程組的方法主要有逆矩陣法、迭代計(jì)算法以及規(guī)劃求解法等。本文以逆矩陣法為例來說明如何求特征向量。

1)選定Sheet3中的區(qū)域A1－F6，輸入公式=P－A(P為特征根形成的對(duì)角線矩陣)，按下Ctrl+Shift+Enter鍵，計(jì)算矩陣λ－A;

2)選擇區(qū)域A8:E12，輸入公式=Minverse(A1:E5)，按下Ctrl+Shift+Enter鍵，計(jì)算前面5個(gè)未知數(shù)形成的矩陣的逆矩陣;

3)在F8:F12中輸入=－F1:F5，作為常數(shù)項(xiàng)(因?yàn)榱顇6=1);

4)矩陣相乘

選取區(qū)域 H8:H12，輸入 =MMULT(A8:E12，F(xiàn)8:F12)，按下Ctrl+Shift+Enter鍵，并在H13中輸入1，求得方程的一組解;

5)對(duì)解向量進(jìn)行標(biāo)準(zhǔn)化

在單元格I8中，輸入公式=H8/SQRT(SUMSQ(MYMHMYM8:MYMHMYM13))(SUMSQ()函數(shù)為求平方和)，拖動(dòng)到I13得到特征向量的值;

將另外的λ代入方程組，求得其他特征根所對(duì)應(yīng)的特征向量。

(2)因子負(fù)荷矩陣的求法

例如，對(duì)于 λ1=2.285，在 J8中輸入 =SQRT(2.285)*I8，拖動(dòng)至J13可求得λ1所對(duì)應(yīng)的主成份因子負(fù)荷向量。同樣求得其他特征根所對(duì)應(yīng)的因子負(fù)荷向量，全部向量構(gòu)成因子負(fù)荷矩陣。

6 結(jié)論

主成份分析法在林業(yè)中應(yīng)用非常廣泛［2～6］。目前有些專業(yè)統(tǒng)計(jì)軟件，如SPSS、SAS等可以進(jìn)行計(jì)算，但這些軟件屬于“傻瓜”型軟件，其操作過程仿佛是在一個(gè)“暗箱”內(nèi)進(jìn)行的，對(duì)于使用者而言，可能對(duì)計(jì)算結(jié)果無法進(jìn)行合理解釋。本文利用Excel進(jìn)行主成份的計(jì)算，使讀者理解其計(jì)算過程，如果再結(jié)合SPSS進(jìn)行計(jì)算結(jié)果的分析，就會(huì)感到理解透徹些。此外，作者寫這篇文章的目的還為了提供一些Excel的基本操作，尤其是有關(guān)矩陣(數(shù)組)的建立及操作等。Excel計(jì)算精度控制方便，操作簡單，且在各類計(jì)算機(jī)上隨處可見，使用者眾多［14－16］，這也是寫這篇文章的目的所在。

［1］ 陳華豪，洪偉等.林業(yè)應(yīng)用數(shù)理統(tǒng)計(jì)［M］.遼寧:大連海運(yùn)學(xué)院出版社，1988，8.第一版:243 ～251.

［2］ 陳松河，鄭清芳.黃甜竹筍用林豐產(chǎn)培育技術(shù)模式的研究［J］.竹子研究匯刊，2001，20(1):61 ～67.

［3］ 高永，孟海生，高國雄，等.吉蘭泰鹽湖引種樟子松的試驗(yàn)研究［J］.內(nèi)蒙古林學(xué)院學(xué)報(bào)(自然科學(xué)版)，1999，21(2):12 ～18.

［4］ 朱春云，趙越，劉霞，等.錦雞兒等旱生樹種抗旱生理的研究［J］.干早區(qū)研究，1996，13(1):59 ～63.

［5］ 謝會(huì)成，姜志林，葉鏡中.麻櫟光合作用的特性及其對(duì)CO2倍增的響應(yīng)［J］.南京林業(yè)大學(xué)學(xué)報(bào)(自然科學(xué)版)，2002，26(4):67～70.

［6］ 陳存及，梁一池，邱爾發(fā)，等.毛竹種源多性狀綜合選擇的研究［J］.林業(yè)科學(xué)，2001，37(1)，Sp.:18 ～23.

［7］ 劉大海，李寧，晁陽編著.SPSS15.0統(tǒng)計(jì)分析送入門到精通［M］.清華大學(xué)出版社，北京，2008:298～317.

［8］ 桂紅義，樊榮，張建森，等.EXCE1 97中文版入門與提高［M］.北京:清華大學(xué)出版社，1997.

［9］ 楊明波，盧建立.Excel中矩陣的命名與特征值特征向量的計(jì)算［J］.電腦知識(shí)與技術(shù)，2007(3):1295～1296.

［10］ 曹永卿，吳宏斌.相關(guān)矩陣的特征值及特征向量的算法［J］.湖南城建高等?？茖W(xué)校學(xué)報(bào)，2001，10(4):64～65.

［11］ 邵榮.用Excel VBA進(jìn)行矩陣特征值的數(shù)值計(jì)算［J］.江南大學(xué)學(xué)報(bào)(自然科學(xué)版)，2009，8(4):418 ～422.

［12］ 盧秋根.中文版Office應(yīng)用基礎(chǔ)教程［M］.上海:上海科學(xué)普及出版社，2005年9月第1版，151～152.

［13］ 李慶楊，王能超，易大義.數(shù)值分析［M］.4版.北京:清華大學(xué)出版社，2001:295～296.

［14］ 陳彥光編著.基于Excel的地理數(shù)據(jù)分析［M］，科學(xué)出版社，北京，2008.

［15］ 陳俊華，慕長龍，朱志芳.Excel在物元模型及層次分析法(AHP)中的應(yīng)用［J］.四川林業(yè)科技，2009，30(5):58 ～62.

［16］ 陳俊華，文吉富，王國良，向成華，等.Excel在計(jì)算群落生物多樣性指數(shù)中的應(yīng)用［J］.四川林業(yè)科技，2009，30(3):88～90，60.