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

薛亞宏，王 嘉，程喜林

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

0 引言

方差分析(Analysis of variance，ANOVA)是英國統計學家兼遺傳學家Gitbert C G提出的一種分析方法，在林業遺傳、科學試驗、醫學研究等眾多領域有極其廣泛的應用[1].方差是統計量分析中的一類，是假設檢驗與區間估計的推廣和延伸，隨著試驗樣本的成因方式不同，一般采用的分析方法也有所差異，在統計學界，通常用單因子(one-way)、復因子(double-way)或多因子(N-way)方法.

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

1 復合因子方差分析

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

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

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

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

算法設計：

基于以上監測數據，現構造出一個6×5型矩陣，命名為矩陣Q,對各組數據采用復因子方差分析，得出以下分析結果：

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

執行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);

在程序運行過程中，anova1( )會自動呈現兩個窗口，分別是盒式圖、分析表，同時顯示概率值p=0.005 51<α，其中α=0.03或0.04表示置信水平，顯然從結果來看應拒絕假設，即藥物對痊愈時間有顯著影響.

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

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

表2 甘肅省小隴山林區巖松等3類松樹生長觀測數據Tab.2 Observation data on the growth of three types of pine trees including Yansong in Xiaolongshan Forest Area, Gansu Province

數據來源：甘肅省小隴山林業實驗局林業科學研究所

算法設計：

根據表中觀測數據，構建矩陣H.然后調用anova2( )函數進行雙因子分析，對各組數據采用雙因子方差分析，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);

執行結果：

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

從結果來看，由于PA=0.013 93，所以應該拒絕H1假設.可以初步推斷，列數據對監測結果有顯著影響，即小隴山林區下轄黨川等4地3類松樹樹種對其胸徑有顯著影響.

以下計算均值，以反映不同樹種在同一林場生長差異：

》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')']

均值計算結果如下：

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

2 主成分分析

主成分分析是一種常見的多因素分析方法，在信息模擬、疾病預防、地理信息采集、工程造價測算、農作物產量分析等領域有著廣泛應用.通常采用SPSS、R等平臺進行分析，但受限于源數據類型的多樣性，輸出圖形的特征的兩極分化(異端非同步)現象較為普遍，經與實際監測比對出現較大偏差，結論不穩定，因此不具有代表性.在這種情況下，利用MATLAB在數據降維處理方面的精度、效度以及在圖象表現上的多維仿真優勢，通過調用Corr( )函數(Corr( )用來計算兩組列向量a和b的相關性，表達式為Corr(a,b))，建立協方差矩陣及特征向量、主成分貢獻率及累計貢獻率、建立變量指標方程組對數據進行降維處理，實現源數據分析從多維到低維的轉化.

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

假設某事件受N個因素(可記為k1,k2,…,kN)影響，監測數據有M組,可構建M×N監測數據矩陣K，再基于K建立協方差矩陣R,其矩陣元素為相關系數ri,j

其中：

根據R計算出特征向量ei及特征值，做反序處理；然后計算主成分貢獻率τi，累計貢獻率δi(基于主成分貢獻率的并行求和)；最后，通過轉換變量(降維)，構建平面坐標方程zi=azi+bzi+Ri,獲取主體變量影響因素間的關系構成，即主成分分析的基本表達式.

2.2 三元主成分降維分析(高程測量方向)

案例3：以甘肅省天水市李子園鉛鋅礦第四紀淺層地貌特征分析為例，某測量點三維坐標參數分別為x=ωcos2ω,y=ωsin2ω,z=0.887x+3.463y,現通過MATLAB生成一維數組，并輸出以2個測量數位為基本單位矢量模擬表達式.

算法設計如下：

》ω=[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)

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

執行以上命令，輸出結果為：

值得注意的是，結果中的d向量、e向量非測量高程測序排列，要通過fliplr()函數和real()函數執行反序和翻轉輸出，目地是使特征值按常規測序呈現，為RNSS測繪系統數據導入做必要的前期配置.

具體語句如下：

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

e=fliplr(e); %對矩陣e進行左右翻轉

D=[d';d';d']; %對矩陣d進行轉置

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

Z=S*M; %三維坐標矩陣與主成分矩陣M做積運算

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

主成分貢獻率γi和累計貢獻率δi可分別用以下公式求得：

本例主要側重對主成分降維分析，故轉換后的3*3矩陣提供二維數據(z列為0)如下：

故新坐標系可表示為：

該坐標方程實現了對三維高程測量數組的降維(縱向投影)，通過坐標轉化使RNSS源數據壓縮于二維平面上，一方面能準確表現該區域第四紀地貌分布特征，另一方面在同類型礦區主要作業區域地形圖繪制(表層、淺層)中提供了滿足繪制精度要求的一種新的計算途徑，其誤差范圍與多基點均勻采樣在同一水平[7].但其在數據生成原理、仿真形式以及中間變量轉換等多個方面集成了ArcGIS、C++的優勢，有效降低了測圖成本.

3 結語

基于MATLAB的方差分析與主成分分析數學原理清晰，算法邏輯性強，語法較為靈活.在實踐中，以數理統計基本實義為基礎，利用矩陣計算、坐標變換等方法，通過MATLAB實現對樣本數據的處理，能有效彌補SPSS、R等工具無法進行數據降維的不足.特別是在三維數字測圖、工程概預算、造價分析等領域能大幅降低數據交互，在一定精度范圍內能有效降低項目成本，具有較強的實用意義.