基于MATLAB的復(fù)合因子方差分析與主成分分析算法

薛亞宏，王 嘉，程喜林

(1.甘肅工業(yè)職業(yè)技術(shù)學(xué)院電信學(xué)院，甘肅天水 741025； 2.蘭州大學(xué)信息科學(xué)與工程學(xué)院，甘肅蘭州 730000；3.西北大學(xué)城市與環(huán)境學(xué)院,陜西西安 741069)

0 引言

方差分析(Analysis of variance，ANOVA)是英國(guó)統(tǒng)計(jì)學(xué)家兼遺傳學(xué)家Gitbert C G提出的一種分析方法，在林業(yè)遺傳、科學(xué)試驗(yàn)、醫(yī)學(xué)研究等眾多領(lǐng)域有極其廣泛的應(yīng)用[1].方差是統(tǒng)計(jì)量分析中的一類，是假設(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)滿足達(dá)到精度即終止計(jì)算.

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

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

在MATLAB中，函數(shù)anova1()、anova2()可分別進(jìn)行單因子、復(fù)因子分析，并能得到較為精確的結(jié)果，以anova1()為例，其基本格式為[p,Tab,Stats]=anova1(Q)[3].其中，Q為需要分析的數(shù)據(jù)，該數(shù)據(jù)是一個(gè)k×w矩陣，其行對(duì)應(yīng)于分組號(hào)，運(yùn)算結(jié)果會(huì)返回概率p、方差表Tab、統(tǒng)計(jì)量Stats;該函數(shù)還將打開兩個(gè)主程序窗口viewer,power,分別以表式、盒式結(jié)構(gòu)呈現(xiàn).

1.1 單因子分析(療效分析方向)

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

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

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

t=[5.6,3.9,6.2322,5.2355,7.1112],p=[0.00551];

執(zhí)行anova1( )：

》Q=[3,5,8,5,2;8,7,6,8,4;5,1,3,5,6;6,8,4,5,6;3,4,4,4,6;3,5,7,4,4];

m=mean(Q);

[p,Tab,Stats]=anova1(Q);

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

1.2 雙因子分析(植物培育方向)

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

表2 甘肅省小隴山林區(qū)巖松等3類松樹生長(zhǎng)觀測(cè)數(shù)據(jù)Tab.2 Observation data on the growth of three types of pine trees including Yansong in Xiaolongshan Forest Area, Gansu Province

數(shù)據(jù)來(lái)源：甘肅省小隴山林業(yè)實(shí)驗(yàn)局林業(yè)科學(xué)研究所

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

根據(jù)表中觀測(cè)數(shù)據(jù)，構(gòu)建矩陣H.然后調(diào)用anova2( )函數(shù)進(jìn)行雙因子分析，對(duì)各組數(shù)據(jù)采用雙因子方差分析，anova2( )命令及矩陣列排列如下：

》H=[25,14,17,30,21,27,31,23,13,19,22,13,31,28,27,13,28,21,13,16,16,18,21,21;

12,23,24,21,16,16,16,23,19,14,21,23,22,23,31,17,29,16,14,19,26,25,14,21;

21,26,20,12,16,19,12,18,18,16,20,13,12,23,14,13,13,21,13,26,22,30,16,19];

anova2(H',6);

執(zhí)行結(jié)果：

Totalmm=0.139 3,ColumnSS=358.1,RowsMS=16.55,Interactiondf=0.475 5

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

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

》D=[];

for m=1:4,for n=1:3,

R(m,n)=mean(H(m,[1:6]+(n-1)*6));

end,end

R=[R;mean(H)];

R=[R mean(R')']

均值計(jì)算結(jié)果如下：

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

2 主成分分析

主成分分析是一種常見的多因素分析方法，在信息模擬、疾病預(yù)防、地理信息采集、工程造價(jià)測(cè)算、農(nóng)作物產(chǎn)量分析等領(lǐng)域有著廣泛應(yīng)用.通常采用SPSS、R等平臺(tái)進(jìn)行分析，但受限于源數(shù)據(jù)類型的多樣性，輸出圖形的特征的兩極分化(異端非同步)現(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ù)(Corr( )用來(lái)計(jì)算兩組列向量a和b的相關(guān)性，表達(dá)式為Corr(a,b))，建立協(xié)方差矩陣及特征向量、主成分貢獻(xiàn)率及累計(jì)貢獻(xiàn)率、建立變量指標(biāo)方程組對(duì)數(shù)據(jù)進(jìn)行降維處理，實(shí)現(xiàn)源數(shù)據(jù)分析從多維到低維的轉(zhuǎn)化.

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

假設(shè)某事件受N個(gè)因素(可記為k1,k2,…,kN)影響，監(jiān)測(cè)數(shù)據(jù)有M組,可構(gòu)建M×N監(jiān)測(cè)數(shù)據(jù)矩陣K，再基于K建立協(xié)方差矩陣R,其矩陣元素為相關(guān)系數(shù)ri,j

其中：

根據(jù)R計(jì)算出特征向量ei及特征值，做反序處理；然后計(jì)算主成分貢獻(xiàn)率τi，累計(jì)貢獻(xiàn)率δi(基于主成分貢獻(xiàn)率的并行求和)；最后，通過(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ì)如下：

》ω=[0:0.2:3*pi]';

x=ω.*cos(2*ω);

y=ω.*sin(2*ω);

z=0.887*x+3.463*y;

S=(x y h);

R=corr(S);

[e,d]=eig(R),d=diag(d);

plot3(x,y,z)

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

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

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

具體語(yǔ)句如下：

》d=(end:-1:1); %對(duì)矩陣d進(jìn)行反序處理

e=fliplr(e); %對(duì)矩陣e進(jìn)行左右翻轉(zhuǎn)

D=[d';d';d']; %對(duì)矩陣d進(jìn)行轉(zhuǎn)置

M=real(sqrt(d)).*e; %對(duì)矩陣d求平方根取實(shí)部

Z=S*M; %三維坐標(biāo)矩陣與主成分矩陣M做積運(yùn)算

plot(Z(:,1),Z(:,2)) %第1、2列散點(diǎn)輸出

主成分貢獻(xiàn)率γi和累計(jì)貢獻(xiàn)率δi可分別用以下公式求得：

本例主要側(cè)重對(duì)主成分降維分析，故轉(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ì)地貌分布特征，另一方面在同類型礦區(qū)主要作業(yè)區(qū)域地形圖繪制(表層、淺層)中提供了滿足繪制精度要求的一種新的計(jì)算途徑，其誤差范圍與多基點(diǎn)均勻采樣在同一水平[7].但其在數(shù)據(jù)生成原理、仿真形式以及中間變量轉(zhuǎn)換等多個(gè)方面集成了ArcGIS、C++的優(yōu)勢(shì)，有效降低了測(cè)圖成本.

3 結(jié)語(yǔ)

基于MATLAB的方差分析與主成分分析數(shù)學(xué)原理清晰，算法邏輯性強(qiáng)，語(yǔ)法較為靈活.在實(shí)踐中，以數(shù)理統(tǒng)計(jì)基本實(shí)義為基礎(chǔ)，利用矩陣計(jì)算、坐標(biāo)變換等方法，通過(guò)MATLAB實(shí)現(xiàn)對(duì)樣本數(shù)據(jù)的處理，能有效彌補(bǔ)SPSS、R等工具無(wú)法進(jìn)行數(shù)據(jù)降維的不足.特別是在三維數(shù)字測(cè)圖、工程概預(yù)算、造價(jià)分析等領(lǐng)域能大幅降低數(shù)據(jù)交互，在一定精度范圍內(nèi)能有效降低項(xiàng)目成本，具有較強(qiáng)的實(shí)用意義.