基于MATLAB的復(fù)合因子方差分析與主成分分析算法設(shè)計(jì)

薛亞宏，王 嘉，程喜林

（1.甘肅工業(yè)職業(yè)技術(shù)學(xué)院 經(jīng)管學(xué)院，甘肅 天水 7410252；2.蘭州大學(xué) 信息科學(xué)與工程學(xué)院，蘭州 7300003；3.西北大學(xué) 城市與環(huán)境學(xué)院，西安 741069）

方差分析（Analysis of variance，ANOVA）是英國(guó)統(tǒng)計(jì)學(xué)家兼遺傳學(xué)家Gitbert提出的一種分析方法，在林業(yè)遺傳、科學(xué)試驗(yàn)、醫(yī)學(xué)研究等眾多領(lǐng)域有極其廣泛的應(yīng)用[1]。方差是統(tǒng)計(jì)量分析中的一類(lèi)，是假設(shè)檢驗(yàn)與區(qū)間估計(jì)的推廣和延伸，由于試驗(yàn)樣本的成因方式不同，一般采用的分法方法也有所差異，在統(tǒng)計(jì)學(xué)界，通常用單因子（one-way）、復(fù)因子（double-way）或多因子（N-way）方法。

主成分分析是根據(jù)已有數(shù)據(jù)推斷假設(shè)數(shù)據(jù)受某主要因素影響程度的一種分析方法，其本質(zhì)仍為假設(shè)檢驗(yàn)[2]。通常需要提供理論樣本結(jié)果和實(shí)測(cè)樣本結(jié)果，兩者通過(guò)矩陣排列得到回歸模型，最終得到影響總體數(shù)據(jù)分布的主要因素及數(shù)值，一般要采用二維或三維曲線進(jìn)行二次以上模擬，在誤差允許范圍內(nèi)滿(mǎn)足達(dá)到精度即終止計(jì)算。

1 復(fù)合因子方差分析

復(fù)因子分析即有兩個(gè)影響因子。通常應(yīng)用于醫(yī)學(xué)療效的驗(yàn)證，其一般的分析步驟為：首先將病人隨機(jī)（一般、等距、整群）分為k組，每組有x人，將每位病人的療效監(jiān)測(cè)指標(biāo)記為tm,n，其中下標(biāo)m、n分別表示第m組，（i=1,2,…,k），n表示某組內(nèi)病人的診療編號(hào)，（n=1,2,…,w），則第m組的第n個(gè)病人的監(jiān)測(cè)指標(biāo)即為tm,n。按照MATLAB特有的表達(dá)式記號(hào)，記第k組所有病人的所有監(jiān)測(cè)指標(biāo)為tm，或各組的第w位病人的監(jiān)測(cè)標(biāo)為tn，則這樣計(jì)算出向量，這涉及算術(shù)平均數(shù)的范籌，如此則構(gòu)造出標(biāo)準(zhǔn)方差分析表，并根據(jù)所給出的監(jiān)測(cè)數(shù)據(jù)找出效果分析數(shù)據(jù)。

以上所采用藥物作為分組的依據(jù)，稱(chēng)為復(fù)因子（Complex factor），它們的差異均值稱(chēng)為復(fù)水平。其中m值與p值較為關(guān)鍵，直接影響概率值p＜a的置信度及拒絕假設(shè)目標(biāo)H0，否則假設(shè)不成立，療效分析驗(yàn)證為假。

在MATLAB中，anoval0、anova2可分別進(jìn)行單因子、復(fù)因子分析，并有效地給出較為精確的結(jié)果，其基本格式為[p,Tab,Stats]=anova1(Q)[3]，其中，Q為需要分析的數(shù)據(jù)，該數(shù)據(jù)是一個(gè)k×w矩陣，其行對(duì)應(yīng)于分組號(hào)，運(yùn)算結(jié)果會(huì)返回檢驗(yàn)值s，以及檢測(cè)數(shù)據(jù)表Tab；該函數(shù)還將打開(kāi)兩個(gè)主程序窗口viewer和power，分別以表式、盒式結(jié)構(gòu)呈現(xiàn)。

1.1 單因子分析（療效分析方向）

案例1：以非嗎啡類(lèi)中樞型鎮(zhèn)痛藥物鹽酸曲馬多（Tramadol）為例[4]，現(xiàn)將40個(gè)病人（醫(yī)學(xué)低于30為小樣本）樣本分為6組，每組5人，患者（patient）使用同一藥物（假定無(wú)其他輔助藥物），記錄從用藥到痊愈時(shí)間（h），觀測(cè)所用藥物的療效是否存在顯著差異，觀測(cè)數(shù)據(jù)如表1所示。

表1 痊愈時(shí)間觀測(cè)數(shù)據(jù)表

算法設(shè)計(jì)：

基于以上監(jiān)測(cè)數(shù)據(jù)，現(xiàn)構(gòu)造出一個(gè)5×6型矩陣，命名為矩陣Q，對(duì)各組數(shù)據(jù)采用復(fù)因子方差分析，得出以下分析結(jié)果：

在程序運(yùn)行過(guò)程中，anoval（）會(huì)自動(dòng)呈現(xiàn)兩個(gè)窗口，分別是盒式圖、分析表，同時(shí)顯示概率值p，其中a=0.03或0.04，表示置信水平，顯然從結(jié)果來(lái)看應(yīng)拒絕假設(shè)，即藥物對(duì)痊愈時(shí)間有顯著影響。

1.2 雙因子分析（植物培育方向）

案例2：以巖松、油松、赤松3種松樹(shù)樹(shù)種在甘肅省天水市小隴山林區(qū)黨川、利橋、草川、草灘4地（林場(chǎng)）的生長(zhǎng)情況為例，每地每類(lèi)樹(shù)種選擇6株，測(cè)量其胸徑，并進(jìn)行雙因子方差分析，觀測(cè)數(shù)據(jù)如表2所示：

表2 甘肅省小隴山林區(qū)巖松等3類(lèi)松樹(shù)生長(zhǎng)觀測(cè)數(shù)據(jù)

Q算法設(shè)計(jì)：

基于以上監(jiān)測(cè)數(shù)據(jù)，現(xiàn)調(diào)用anoval2（）函數(shù)，命名為矩陣H，對(duì)各組數(shù)據(jù)采用雙因子方差分析，anoval2（）命令及矩陣列排列如下：

從結(jié)果來(lái)看，由于PA=0.01393，所以應(yīng)該拒絕H1假設(shè)。可以初步推斷，列數(shù)據(jù)對(duì)監(jiān)測(cè)結(jié)果有顯著影響，即小隴山林區(qū)下轄黨川等4地3類(lèi)松樹(shù)樹(shù)種對(duì)其胸徑有顯著影響。

以下計(jì)算均值，以反映不同樹(shù)種在同一林場(chǎng)生長(zhǎng)差異：

根據(jù)結(jié)果分析，赤松胸徑明顯大于巖松和油松。PH與PHQ的差距較大，從而判斷假設(shè)為真，故接受假設(shè)[5]。即黨川等4地各自對(duì)3類(lèi)松樹(shù)樹(shù)種的胸徑有輕微影響，不同區(qū)域（林場(chǎng)）對(duì)不同松樹(shù)樹(shù)種胸徑成長(zhǎng)無(wú)顯著影響。

2 主成分分析

主成分分析是一種常見(jiàn)的多因素分析方法，在信息模擬、疾病預(yù)防、地理信息采集、工程造價(jià)測(cè)算、農(nóng)作物產(chǎn)量分析等領(lǐng)域有著廣泛應(yīng)用[6]。通常采用SPSS、R等平臺(tái)進(jìn)行分析，但由于原數(shù)據(jù)類(lèi)型的多樣性，輸出圖形特征的兩極分化（異端非同步）現(xiàn)象較為普遍，經(jīng)與實(shí)際監(jiān)測(cè)比對(duì)出現(xiàn)較大偏差，結(jié)論不穩(wěn)定，因此不具有代表性。在這種情況下，利用MATLAB在數(shù)據(jù)降維處理方面的精度、效度以及在圖象表現(xiàn)上的多維仿真優(yōu)勢(shì)，通過(guò)調(diào)用Corr（）函數(shù)，建立協(xié)方差矩陣及特征向量、主成分貢獻(xiàn)率及累計(jì)貢獻(xiàn)率、建立變量指標(biāo)方程組對(duì)數(shù)據(jù)進(jìn)行降維處理，實(shí)現(xiàn)源數(shù)據(jù)分析從多維到低維的轉(zhuǎn)化。

2.1 基于MATLAB的主成分分析算法原理

首先建立協(xié)方差矩陣

再根據(jù)R計(jì)算出特征向量ei及特征值，做反序處理；

最后，通過(guò)轉(zhuǎn)換變量（降維），構(gòu)建平面坐標(biāo)方程zi=azi+bzi+Ri，獲取主體變量影響因素間的關(guān)系構(gòu)成，即主成分分析的基本表達(dá)式。

2.2 三元主成分分析（高程測(cè)量方向）

案例3：以甘肅省天水市李子園鉛鋅礦第四紀(jì)淺層地貌特征分析為例，某測(cè)量點(diǎn)三維坐標(biāo)參數(shù)分別為 x=ωcos2ω, y=ωsin2ω, z=0.887x+3.463y，現(xiàn)通過(guò)MATLAB生成一維數(shù)組，并輸出以2個(gè)測(cè)量數(shù)位為基本單位矢量模擬表達(dá)式。

算法設(shè)計(jì)如下：

顯然，基于對(duì)降維矢量輸出原理的分析，進(jìn)一步利用空間坐標(biāo)變換，對(duì)三維原數(shù)據(jù)做放樣投影，得到二維數(shù)組[7]。

執(zhí)行以上命令，輸出結(jié)果為：

值得注意的是，結(jié)果中的σ向量、e向量非測(cè)量高程測(cè)序排列，要通過(guò)fliplr（）函數(shù)和real（）函數(shù)執(zhí)行反序和翻轉(zhuǎn)輸出，目地是使特征值按常規(guī)測(cè)序呈現(xiàn)，為RNSS測(cè)繪系統(tǒng)數(shù)據(jù)導(dǎo)入做必要的前期配置。

具體語(yǔ)句如下：

轉(zhuǎn)換后的3×3矩陣提供二維數(shù)據(jù)（z列為0）如下：

故新坐標(biāo)系可表示為：

該坐標(biāo)方程實(shí)現(xiàn)了對(duì)三維高程測(cè)量數(shù)組的降維（縱向投影），通過(guò)坐標(biāo)轉(zhuǎn)化使RNSS源數(shù)據(jù)壓縮于二維平面上，一方面能準(zhǔn)確表現(xiàn)該區(qū)域第四紀(jì)地貌分布特征，另一方面在同類(lèi)型礦區(qū)主要作業(yè)區(qū)域地形圖繪制（表層、淺層）中提供了滿(mǎn)足繪制精度要求的一種新的計(jì)算途徑，其誤差范圍與多基點(diǎn)均勻采樣在同一水平[8]，但其在數(shù)據(jù)生成原理、仿真形式以及中間變量轉(zhuǎn)換等多個(gè)方面集成了ArcGIS、C++的優(yōu)勢(shì)，有效降低了測(cè)圖成本。

3 結(jié)語(yǔ)

基于MATLAB的復(fù)合因子方差分析與主成分分析計(jì)算原理清晰，算法邏輯性強(qiáng)，語(yǔ)法調(diào)用靈活。在實(shí)踐中，以應(yīng)用統(tǒng)計(jì)學(xué)基本計(jì)算理論為基礎(chǔ)，結(jié)合矩陣運(yùn)算、坐標(biāo)變換等數(shù)學(xué)手段，最終通過(guò)MATLAB實(shí)現(xiàn)對(duì)樣本數(shù)據(jù)的終端處理，有效彌補(bǔ)了SPSS、R等分析工具在圖像擬合優(yōu)度與坐標(biāo)維度無(wú)法兼顧的不足，特別是三維數(shù)字測(cè)圖、工程概預(yù)算、造價(jià)分析等領(lǐng)域內(nèi)能大大地降低數(shù)據(jù)交叉，有效降低項(xiàng)目成本，有較強(qiáng)的實(shí)用意義。