方差分析及MATLAB實現(xiàn)目標

方差分析及MATLAB實現(xiàn)目標第1章方差分析

在工農(nóng)業(yè)生產(chǎn)和科學(xué)研究中,經(jīng)常遇到這樣的問題:影響產(chǎn)品產(chǎn)量、質(zhì)量的因素很多,我們需要了解在這眾多的因素中,哪些因素對影響產(chǎn)品產(chǎn)量、質(zhì)量有顯著影響.為此,要先做試驗,然后對測試的結(jié)果進行分析.方差分析就是分析測試結(jié)果的一種方法.

在方差分析中,把在試驗中變化的因素稱為因子,用A、B、C、...表示;因子在試驗中所取的不同狀態(tài)稱為水平,因子A的r個不同水平用A1、A2、...、Ar表示.§1單因子方差分析§1.1基本概念與數(shù)學(xué)模型水平觀測值A(chǔ)1x11x12...x1n1A2x21x22…x2n2……………Arxr1xr2…xrnr

例:為尋求適應(yīng)本地區(qū)的高產(chǎn)油菜品種,今選了五種不同品種進行試驗,每一品種在四塊試驗田上得到在每一塊田上的畝產(chǎn)量如下:

我們要研究的問題是諸不同品種的平均畝產(chǎn)量是否有顯著差異.

試驗的目的就是要檢驗假設(shè)

H0:μ1=μ2=μ3=μ4=μ5是否成立.若是拒絕,那么我們就認為這五種品種的平均畝產(chǎn)量之間有顯著差異;反之,就認為各品種間產(chǎn)量的不同是由隨機因素引起的.方差分析就是檢驗假設(shè)的一種方法.

在本例中只考慮品種這一因子對畝產(chǎn)量的影響,五個不同品種就是該因子的五個不同水平.由于同一品種在不同田塊上的畝產(chǎn)量不同,我們可以認為一個品種的畝產(chǎn)量就是一個總體,在方差分析中總假定各總體獨立地服從同方差正態(tài)分布,即第i個品種的畝產(chǎn)量是一個隨機變量,它服從分布N(μi,σ2),i=1,2,3,4,5.

設(shè)在某試驗中,因子A有r個不同水平A1,A2,...,Ar,在Ai水平下的試驗結(jié)果Xi服從正態(tài)分布N(μi,σ2),i=1,2,...,r,且X1,X2,...,Xr間相互獨立.現(xiàn)在水平Ai下做了ni次試驗,獲得了ni個試驗結(jié)果Xij,j=1,2,...,ni這可以看成是取自Xi的一個容量為ni的樣本,i=1,2,...,r.

實際上,方差分析是檢驗同方差的若干正態(tài)總體均值是否相等的一種統(tǒng)計方法.

在實際問題中影響總體均值的因素可能不止一個.我們按試驗中因子的個數(shù),可以有單因子方差分析,雙因子分析,多因子分析等.例中是一個單因子方差分析問題.水平觀測值A(chǔ)1x11x12...x1n1A2x21x22…x2n2……………Arxr1xr2…xrnr

由于Xij~N(μi,σ2),故Xij與μi的差可以看成一個隨機誤差εij~N(0,σ2).這樣一來,可以假定Xij具有下述數(shù)據(jù)結(jié)構(gòu)式:

為了今后方便起見,把參數(shù)的形式改變一下,并記稱μ為一般平均,αi為因子A的第i個水平的效應(yīng).Xij=μi+εij,i=1,2,...,r;j=1,2,...,ni其中諸εij~N(0,σ2),且相互獨立.要檢驗的假設(shè)是

H0:μ1=μ2=…=μr

在這樣的改變下,單因子方差分析模型中的數(shù)據(jù)結(jié)構(gòu)式可以寫成:所要檢驗的假設(shè)可以寫成:

為了導(dǎo)出檢驗假設(shè)的統(tǒng)計量,下面我們分析一下什么是引起諸Xij

波動的原因.

平方和分解公式：引起諸Xij

波動的原因有兩個:一個是假設(shè)H0為真時,諸Xij的波動純粹是隨機性引起的;另一個可能是假設(shè)H0不真而引起的.因而我們就想用一個量來刻劃諸Xij之間的波動,并把引起波動的兩個原因用另兩個量表示出來,這就是方差分析中常用的平方和分解法.§1.2統(tǒng)計分析其中交叉乘積項下面我們來看各式的意義

檢驗統(tǒng)計量的構(gòu)造:對于各組樣本有因此

一般,當F>F時,稱因子的影響高度顯著,記為“**”;當F>F≥F時,稱因子的影響顯著,記為“*”;當F＜F時,稱因子無顯著影響,即認為因子各水平間無差異.

檢驗過程：§1.3方差分析表

例:為尋求適應(yīng)本地區(qū)的高產(chǎn)油菜品種,今選了五種不同品種進行試驗,每一品種在四塊試驗田上得到在每一塊田上的畝產(chǎn)量如下:

我們要研究的問題是諸不同品種的平均畝產(chǎn)量是否有顯著差異.

解:先列表計算例:下面給出了隨機選取的,用于計算器的四種類型的電路的響應(yīng)時間(以毫秒計).

表:電路的響應(yīng)時間類型I類型II類型III類型IV15222018402133271617151826182219這里試驗的指標是電路的響應(yīng)時間.電路類型為因素.這一因素有四個水平,試驗的目的是要考察各類型電路對響應(yīng)時間的影響.設(shè)四種類型電路的響應(yīng)時間的總體均為正態(tài),且各總體方差相同,但參數(shù)均未知.又設(shè)各樣本相互獨立.解

分別以m1,m2,m3,m4記類型I,II,III,IV四種電路響應(yīng)時間總體的平均值.我們需檢驗(a=0.05)

H0:m1=m2=m3=m4,

H1:m1,m2,m3,m4不全相等.

現(xiàn)在n=18,s=4,n1=n2=n3=5,n4=3,試驗號12345和和平方類型I1915222018948836類型II204021332714119881類型III1617151826928464類型IV182219593481ST,SA,SE的自由度依次為17,3,14

表:方差分析表方差來源平方和自由度均方F值顯著性因素318.983106.333.76*誤差395.461428.25總和714.4417因F(3,14)=3.34<3.76<F,故認為各類型電路的響應(yīng)時間有顯著差異.1.4.單因素方差分析的Matlab實現(xiàn)單因素方差分析:anova1調(diào)用格式:(1)p=anova1(X)(2)p=anova1(X,group)(3)p=anova1(X,group,displayopt)(4)[p,table]=anova1(...)(5)[p,table,stats]=anova1(...)（2）p=anova1(X,group)輸入：X是一個向量，從第一個總體的樣本到第r個總體的樣本依次排列，group是與X有相同長度的向量，表示X中的元素是如何分組的.group中某元素等于i,表示X中這個位置的數(shù)據(jù)來自第i個總體.因此group中分量必須取正整數(shù)，從1直到r.（1）p=anova1(X)%比較X中各列數(shù)據(jù)的均值是否相等。此時輸出的p是零假設(shè)成立時，數(shù)據(jù)的概率，當稱差異是顯著的，當稱差異是高度顯著的.輸入X各列的元素相同，即各總體的樣本大小相等，稱為均衡數(shù)據(jù)的方差分析。不均衡時用下面的命令：Table——輸出anova表：stats——輸出boxplot圖：0.30352.87173.57744.98464.9438];>>p=anova1(X)例.某水產(chǎn)研究所為了比較四種不同配合飼料對魚的飼喂效果，選取了條件基本相同的魚20尾，隨機分成四組，投喂不同飼料，經(jīng)一個月試驗以后，各組魚的增重結(jié)果列于下表。表飼喂不同飼料的魚的增（單位：10g）飼料魚的增重（xij）A131.927.931.828.435.9A224.825.726.827.926.2A322.123.627.324.925.8A427.030.829.024.528.5四種不同飼料對魚的增重效果是否顯著？解：這是單因素均衡數(shù)據(jù)的方差分析，Matlab程序如下：27.0 30.8 29.0 24.5 28.5];%原始數(shù)據(jù)輸入B=A';%將矩陣轉(zhuǎn)置,Matlab中要求各列為不同水平p=anova1(B)

運行后得到一表一圖，表是方差分析表（重要）；圖是各列數(shù)據(jù)的盒子圖，離盒子圖中心線較遠的對應(yīng)于較大的F值，較小的概率p.Source方差來源SS平方和df自由度MS均方差F統(tǒng)計量P值Columns(因素A組間)SSAr-1SS/(r-1)7.140.0029Error誤差（組內(nèi)）SSEn-rSS/(n-r)Total總和SSTn-1表中所列出的各項意義如下：因為，故不同飼料對魚的增重效果極為顯著.如果沒有給出概率。四種不同飼料對魚的增重效果極為顯著，那么哪一種最好呢？請看下圖

此時，第一個圖對應(yīng)第一種飼料且離盒子圖中心線較遠，效果最突出。如果從原始數(shù)據(jù)中去掉第一種飼料的試驗數(shù)據(jù)，得到的結(jié)果為各種飼料之間對魚的增重效果不顯著.p=anova1(B(:,2:4))

例.為比較同一類型的三種不同食譜的營養(yǎng)效果，將19支幼鼠隨機分為三組，各采用三種食譜喂養(yǎng).12周后測得體重，三種食譜營養(yǎng)效果是否有顯著差異？食譜體重增加量甲164190203205206214228257乙185197201231丙187212215220248265281解：這是單因素非均衡數(shù)據(jù)的方差分析A=[164190203205206214228257185197201231187212215220248265281];group=[ones(1,8),2*ones(1,4),3*ones(1,7)];p=anova1(A,group)方差分析表:

均值盒子圖由于概率比較大，故認為三種食料沒有顯著差異.(3)多重比較的MATLAB實現(xiàn)

為了便于解決實際問題，我們給出多重比較的MATLAB命令。c=multcompare(s)其中輸入s，由[p,c,s]=anova1(B);得到輸出C共有6列，每一行給出均值差的置信區(qū)間例.四個實驗室試制同一型號紙張，為了比較光滑度每個實驗室測量了8張紙，進行方差分析

實驗室紙張光滑度A138.741.543.844.545.54647.758A239.239.339.741.441.842.943.345.8A33435394043434445A43434.834.835.437.237.841.242.8解：A=[38.7,41.5,43.8,44.5,45.5,46,47.7,5834,35,39,40,43,43,44,4534,34.8,34.8,35.4,37.2,37.8,41.2,42.8];%輸入數(shù)據(jù)B=A’;%MATLAB只對各列進行分析[p,c,s]=anova1(b);%方差分析c=multcompare(s)%多重比較從方差分析表可知：四個實驗室生產(chǎn)有差異，那么如何比較？軟件輸出c如下所示：1,2列表示比較的實驗室號碼，3，5列分別為置信區(qū)間左右端點，第4列是均值差的統(tǒng)計量觀測值.

1.00002.0000-1.47534.03759.55031.00003.0000-0.17535.337510.85031.00004.00002.94978.462513.97532.00003.0000-4.21281.30006.81282.00004.0000-1.08784.42509.93783.00004.0000-2.38783.12508.6378若置信區(qū)間包含原點則無顯著差異，可見只有1,4實驗室有顯著差異.另外，軟件輸出一幅圖形，告知1，4有顯著差異.§2雙因子方差分析§2.1數(shù)學(xué)模型§2.2無交互影響的雙因子方差分析SE表示試驗的隨機波動引起的誤差,稱為誤差平方和;SA除了反映了試驗的隨機波動引起的誤差外,還反映了因子A的效應(yīng)間的差異,稱為因子A的偏差平方和;SB除了反映了試驗的隨機波動引起的誤差外,還反映了因子B的效應(yīng)間的差異,稱為因子B的偏差平方和.

具體計算時可用計算表和方差分析表:

一般,當F>F時,稱因子的影響高度顯著,記為“**”;當F>F≥F時,稱因子的影響顯著,記為“*”;當F＜F時,稱因子無顯著影響,即認為因子各水平間無差異.

例:為了考察蒸餾水的pH值和硫酸銅溶液濃度對化驗血清中白蛋白與球蛋白的影響,對蒸餾水的pH值(A)取了4個不同水平,對硫酸銅溶液濃度(B)取了3個不同水平,在不同水平組合(Ai,Bj)下各測一次白蛋白與球蛋白之比,其結(jié)果列于計算表的左上角.試檢驗兩因子對化驗結(jié)果有無顯著差異.解查F-分布表得:F(3,6)=4.76,F(2,6)=5.14,F(3,6)=9.78,F(2,6)=10.9,由此可知FA>F(3,6);FB>F(2,6).所以因子A及因子B的不同水平對化驗結(jié)果有高度顯著影響.§2.3有交互作用的雙因子方差分析其中n=rst

仍然用平方和分解的思想來給出檢驗用的統(tǒng)計量，先引入下述記號：由此可知總的偏差平方和可作如下的分解:其中各偏差平方和表達式如下:

各偏差平方和的意義:SE表示試驗的隨機波動引起的誤差,稱為誤差平方和;SA除了反映了試驗的隨機波動引起的誤差外,還反映了因子A的效應(yīng)間的差異,稱為因子A的偏差平方和;SB除了反映了試驗的隨機波動引起的誤差外,還反映了因子B的效應(yīng)間的差異,稱為因子B的偏差平方和;SA×B除了反映了試驗的隨機波動引起的誤差外,還反映了交互效應(yīng)的差異所引起的波動,稱為交互作用的偏差平方和.同無交互作用的情況類似可得:

檢驗統(tǒng)計量及顯著性檢驗:

這就是用來檢驗假設(shè)H01,H02,H03,的統(tǒng)計量.按照顯著性假設(shè)檢驗程序,對給定的顯著性水平α,

當FA>F1-α(r-1,rs(t-1))時拒絕H01;

當FB>F1-α(s-1,rs(t-1))時拒絕H02;

當

FA×B>F1-α((r-1)(s-1),rs(t-1))時拒絕H03.

具體的計算過程,各偏差平方和的計算也可用下面簡化的表達式,且可列成一張計算表和方差分析表.

一般,當F>F時,稱因子的影響高度顯著,記為“**”;當F0.99>F≥F時,稱因子的影響顯著,記為“*”;當F＜F時,稱因子無顯著影響,即認為因子各水平間無差異.

例:在某化工生產(chǎn)中為了提高收率,選了三種不同濃度,四種不同溫度做試驗.在同一濃度與同一溫度組合下各做二次試驗,其收率數(shù)據(jù)如下而計算表所列(數(shù)據(jù)均已減去75).試檢驗不同濃度,不同溫度以及它們間的交互作用對收率有無顯著影響.解:查表知F(2,12)=3.89,F(2,12)=6.93;F(3,12)=3.49,F(3,12)=5.95;F(6,12)=3.00,F(6,12)=4.81.由此知F<FA<F,而FB<F,FA×B<F.故濃度不同將對收率產(chǎn)生顯著影響;而溫度和交互作用的影響都不顯著.2.4.雙因素的方差分析的MATLAB實現(xiàn)

在Matlab中雙因素的方差分析命令如下：雙因素方差分析:anova2調(diào)用格式:(1)p=anova2(X)(2)p=anova2(X,reps)(3)p=anova2(X,reps,displayopt)(4)[p,table]=anova1(...)(5)[p,table,stats]=anova1(...)

在Matlab中雙因素有交互作用的方差分析命令如下：[p,t,s]=anova2(X,resp)其中輸入X是一個矩陣；resp表示試驗的重復(fù)次數(shù)輸出的p值有三個，分別為各行、各列以及交互作用的概率.若p<0.05,有顯著差異若p<0.01,有高度顯著差異t是方差分析表，s用于各因素均值估計與比較.82

雙因素方差分析:anova2

調(diào)用格式:(1)p=anova2(X)(2)p=anova2(X,reps)(3)p=anova2(X,reps,displayopt)(4)[p,table]=anova1(...)(5)[p,table,stats]=anova1(...)7.05.04.5];>>p=anova2(X,3)p=2023年6月23日MATLAB和R軟件83例一火箭使用了4種燃料，3種推進器作射程試驗，每種燃料與每種推進器的組合各發(fā)射火箭2次，得到結(jié)果如下：B1B2B3A158.2,52.656.2,41.265.3,60.8A249.1,42.854.1,50.551.6,48.4A360.1,58.370.9,73.239.2,40.7A475.8,71.558.2,51.048.7,41.4試在水平下，檢驗不同燃料（因素A）、不同推進器（因素B）下的射程是否有顯著差異？交互作用是否顯著？84解編寫程序如下：clc,clear75.8,71.558.2,51.048.7,