第9章概率論與數(shù)理統(tǒng)計的MATLAB實現(xiàn)講稿

第9章概率論與數(shù)理統(tǒng)計的MATLAB實現(xiàn)-PAGE245-第9章概率論與數(shù)理統(tǒng)計的MATLAB實現(xiàn)MATLAB總包提供了一些進行數(shù)據(jù)統(tǒng)計分析的函數(shù)，但不完整。利用MATLAB統(tǒng)計工具箱，可以進行概率和數(shù)理統(tǒng)計分析，以及進行比較復雜的多元統(tǒng)計分析。9.1隨機變量及其分布利用統(tǒng)計工具箱提供的函數(shù)，可以比較方便地計算隨機變量的分布列（或密度函數(shù)）和分布函數(shù)。9.1.1常見離散型隨機變量的分布列如果隨機變量全部可能取到的不相同的值是有限個或可列無限多個，則稱為離散型隨機變量。MATLAB提供的計算常見離散型隨機變量分布列的函數(shù)及調(diào)用格式：函數(shù)調(diào)用格式(對應的分布)分布列y=binopdf(x,n,p)(二項分布)y=geopdf(x,p)(幾何分布)y=hygepdf(x,M,K,n)(超幾何分布)y=poisspdf(x,lambda)(泊松分布)y=unidpdf(x,n)(離散均勻分布)9.1.2常見連續(xù)型隨機變量的密度函數(shù)對于隨機變量的分布函數(shù)，如果存在非負函數(shù)，使對于任意實數(shù)有則稱為連續(xù)型隨機變量，其中函數(shù)稱為的密度函數(shù)。MATLAB提供的計算常見連續(xù)型隨機變量分布密度函數(shù)的函數(shù)及調(diào)用格式：函數(shù)調(diào)用格式(對應的分布)密度函數(shù)y=betapdf(x,a,b)(分布)y=chi2pdf(x,v)(卡方分布)y=exppdf(x,mu)(指數(shù)分布)y=fpdf(x,v1,v2)(F分布)y=gampdf(x,a,b)(伽馬分布)y=normpdf(x,mu,sigma)(正態(tài)分布)y=lognpdf(x,mu,sigma)(對數(shù)正態(tài)分布)y=raylpdf(x,b)(瑞利分布)y=tpdf(x,v)(學生氏t分布)y=unifpdf(x,a,b)(連續(xù)均勻分布)y=weibpdf(x,a,b)(威布爾分布)比如，用normpdf函數(shù)計算正態(tài)概率密度函數(shù)值。該函數(shù)的調(diào)用格式為：Y=normpdf(X,MU,SIGMA)計算數(shù)據(jù)X中各值處參數(shù)為MU和SIGMA的正態(tài)概率密度函數(shù)的值。參數(shù)SIGMA必須為正。正態(tài)概率密度函數(shù)的計算公式為：9.1.3用函數(shù)pdf計算隨機變量的分布列或概率除了用上述的函數(shù)計算服從相應分布的隨機變量的分布列或概率密度外，還可以用函數(shù)pdf計算隨機變量的分布列或概率密度。調(diào)用格式：Y=pdf('name',X,A1,A2,A3)返回服從參數(shù)為A1,A2,A3的'name'分布的隨機變量在X處的分布列或密度函數(shù)值。Y與X同型，分布函數(shù)名'name'常見的取值如下：'beta'或'Beta'：Beta分布'bino'或'Binomial'：二項分布'chi2'或'Chisquare'：卡方分布'exp'或'Exponential'：指數(shù)分布'f'或'F'：F分布'gam'或'Gamma'：GAMMA分布'geo'或'Geometric'：幾何分布'hyge'或'Hypergeometric'：超幾何分布'logn'或'Lognormal'：對數(shù)正態(tài)分布'nbin'或'NegativeBinomial'：負二項分布'ncf'或'NoncentralF'：非中心F分布'nct'或'Noncentralt'：非中心t分布'ncx2'或'NoncentralChi-square'：非中心卡方分布'norm'或'Normal'：正態(tài)分布'poiss'或'Poisson'：泊松分布'rayl'或'Rayleigh'：瑞利分布't'或'T'：T分布'unif'或'Uniform'：均勻分布'unid'或'DiscreteUniform'：離散均勻分布'weib'或'Weibull'：Weibull分布比如，計算自由度為8的卡方分布,在點2.18處的密度函數(shù)值的命令為：pdf('chi2',2.18,8)

9.1.對于離散型隨機變量，設為任意實數(shù)，的分布函數(shù)為：對于連續(xù)型隨機變量，假設其概率密度函數(shù)為，則其分布函數(shù)為：對常見分布的隨機變量，MATLAB均提供了專門的函數(shù)來計算它們各自的分布函數(shù)，這些函數(shù)是具體如下：函數(shù)調(diào)用格式對應的分布F=betacdf(x,a,b)分布F=binocdf(x,n,p)二項分布F=chi2cdf(x,v)卡方分布F=expcdf(x,mu)指數(shù)分布F=fcdf(x,v1,v2)F分布F=gamcdf(x,a,b)伽馬分布F=geocdf(x,p)幾何分布F=hygecdf(x,M,K,n)超幾何分布F=normcdf(x,mu,sigma)正態(tài)分布F=logncdf(x,mu,sigma)對數(shù)正態(tài)分布F=poisscdf(x,lambda)泊松分布F=raylcdf(x,b)瑞利分布F=tcdf(x,v)學生氏t分布F=unidcdf(x,n)離散均勻分布F=unifcdf(x,a,b)連續(xù)均勻分布F=weibcdf(x,a,b)威布爾分布例如，用normcdf函數(shù)計算正態(tài)分布的分布函數(shù)。該函數(shù)的調(diào)用格式為：F=normcdf(X,MU,SIGMA)計算參數(shù)為MU和SIGMA的正態(tài)分布函數(shù)在數(shù)據(jù)X中每個值處的值。參數(shù)SIGMA必須為正。正態(tài)分布的分布函數(shù)為：結果為取自參數(shù)為和的正態(tài)分布總體的單個觀測量落在區(qū)間中的概率。另外，還可以用函數(shù)cdf計算隨機變量的分布函數(shù)。調(diào)用格式：F=cdf('name',X,A1,A2,A3)返回服從參數(shù)為A1,A2,A3的'name'分布的隨機變量在X處的分布函數(shù)值。分布函數(shù)名'name'常見的取值同函數(shù)pdf中的'name'。例9－1某儀器需安裝一個電子元件，需要電子元件的使用壽命不低于1000小時即可?，F(xiàn)有甲乙兩廠的電子元件可供選擇，甲廠生產(chǎn)的電子元件的壽命服從正態(tài)分布，乙廠生產(chǎn)的電子元件的壽命服從正態(tài)分布。問應選哪個工廠的產(chǎn)品呢？解：設，。則有：0.97720.9696因此，應選甲廠生產(chǎn)的產(chǎn)品。注：計算的命令為：1-normcdf(1000,1100,50)或1-cdf('norm',1000,1100,50)計算的命令為：1-normcdf(1000,1150,80)或1-cdf('norm',1000,1150,80)9.1.5分布函數(shù)的MATLAB中，常見分布的分布函數(shù)的逆函數(shù)及其調(diào)用格式：函數(shù)調(diào)用格式對應的分布x=betainv(P,a,b)分布x=binoinv(P,n,p)二項分布x=chi2inv(P,v)卡方分布x=expinv(P,mu)指數(shù)分布x=finv(P,v1,v2)F分布x=gaminv(P,a,b)伽馬分布x=geoinv(P,p)幾何分布x=hygeinv(P,M,K,n)超幾何分布x=norminv(P,mu,sigma)正態(tài)分布x=logninv(P,mu,sigma)對數(shù)正態(tài)分布x=poissinv(P,lambda)泊松分布x=raylinv(P,b)瑞利分布x=tcdfinv(P,v)學生氏t分布x=unidinv(P,n)離散均勻分布x=unifinv(P,a,b)連續(xù)均勻分布x=weibinv(P,a,b)威布爾分布在MATLAB中，還可以用函數(shù)icdf計算隨機變量的分布函數(shù)的逆函數(shù)。調(diào)用格式：X=icdf('name',P,A1,A2,A3)服從參數(shù)為A1,A2,A3的'name'分布的隨機變量的分布函數(shù)在X處值為P。分布函數(shù)名'name'常見的取值同函數(shù)pdf中的'name'。例9－2有同類設備300臺，各臺工作狀態(tài)相互獨立。已知每臺設備發(fā)生故障的概率為0.01，若一臺設備發(fā)生故障需要1人去處理，問至少需要多少工人，才能保證設備發(fā)生故障而不能及時維修的概率小于0.01？解：設表示同一時刻發(fā)生故障的設備臺數(shù)，則有。再設配備位維修人員，則有：即鍵入命令：x=binoinv(0.99,300,0.01)或x=icdf('bino',0.99,300,0.01)運行結果：x=8鍵入命令：F1=binocdf(8,300,0.01),F2=binocdf(7,300,0.01)運行結果：F1=0.9964，F(xiàn)2=0.9885。因此，至少需要8個工人，才能保證設備發(fā)生故障而不能及時維修的概率小于0.01。例9－3服從卡方分布的隨機變量的分布函數(shù)的逆函數(shù)的應用程序代碼：n=5;a=0.05;%n為自由度x_a=chi2inv(1-a,n);%x_a為臨界值x=linspace(0,20,1000);y_pdf=chi2pdf(x,n);%計算的概率密度函數(shù)值,供繪圖用.plot(x,y_pdf,'b')%繪密度函數(shù)圖形holdonxx=linspace(0,x_a,800);yy_pdf=chi2pdf(xx,n);%計算[0,x_a]上的密度函數(shù)值,供填色用fill([xx,x_a],[yy_pdf,0],'g')%填色,其中:點(x_a,0)使得填色區(qū)域封閉.text(x_a+0.01,0.005,num2str(x_a))%標注臨界值點text(10,0.10,['\fontsize{16}X~{\chi}^2(5)'])%圖中標注text(1.5,0.03,'\fontsize{16}1-alpha=0.95')%圖中標注運行結果見圖9―1。圖9―19.2多維隨機變量及其分布9.2.1用mvnpdf和mvncdf函數(shù)可以計算二維正態(tài)分布隨機變量在指定位置處的密度函數(shù)值和分布函數(shù)值。例9－4計算服從二維正態(tài)分布的隨機變量在指定范圍內(nèi)的概率密度值并繪圖。程序代碼：%二維正態(tài)分布的隨機變量在指定范圍內(nèi)的概率密度函數(shù)圖形mu=[00];sigma=[0.250.3;0.31];%協(xié)方差陣x=-3:0.1:3;y=-3:0.15:3;[x1,y1]=meshgrid(x,y);%將平面區(qū)域網(wǎng)格化取值f=mvnpdf([x1(:)y1(:)],mu,sigma);%計算二維正態(tài)分布概率密度函數(shù)值F=reshape(f,numel(y),numel(x));%矩陣重塑surf(x,y,F);%繪刻面圖caxis([min(F(:))-0.5*range(F(:)),max(F(:))]);%設置顏色的范圍%range(x)表示最大值與最小值的差，即極差。axis([-44-440max(F(:))+0.1]);%設置坐標軸范圍xlabel('x')ylabel('y')zlabel('ProbabilityDensity')運行結果見圖9－2。圖9－2二維正態(tài)分布的隨機變量的密度函數(shù)圖形

9.2.2若連續(xù)型隨機變量的密度函數(shù)為，則關于和的邊緣概率密度和分別為：例9－5設具有概率密度⑴確定常數(shù)；⑵求邊緣概率密度和。解：⑴由可得；計算程序代碼：clear;clc;symsxyCfxy=C*x^2*y;g=int(int(fxy,y,x*x,1),x,-1,1);C=double(solve(g-1))⑵程序代碼：clear;clc;symsxyfxy=5.25*x*x*y;fx=int(fxy,y,x*x,1)fy=int(fxy,x,-sqrt(y),sqrt(y))運行結果：fx=21/8*x^2*(1-x^4)fy=7/2*y^(5/2)因此，，

9.3隨機變量的數(shù)字特征在解決實際問題過程中，往往并不需要全面了解隨機變量的分布情況，而只需要知道它們的某些特征，這些特征通常稱為隨機變量的數(shù)字特征。常見的有數(shù)學期望、方差、相關系數(shù)和矩等。9.3.1⒈離散型隨機變量的數(shù)學期望設離散型隨機變量的分布律為：，如果絕對收斂，則稱的和為隨機變量的數(shù)學期望。例9－6設表示一張彩票的獎金額，的分布列如下：500000500005000500501000.0000010.0000090.000090.00090.0090.090.9試求。求解程序代碼：%離散型隨機變量的數(shù)學期望clear;clc;x=[50000050000500050050100]';p=[0.0000010.0000090.000090.00090.0090.090.9]';Ex=x'*p運行結果：Ex=3.2000例9－7設，，，求。求解程序：%離散型隨機變量的數(shù)學期望clear;clc;symspkEx=symsum(k*p*(1-p)^(k-1),k,1,inf)運行結果：Ex=1/p⒉連續(xù)型隨機變量的數(shù)學期望設連續(xù)型隨機變量的概率密度為，若積分絕對收斂，則稱該積分的值為隨機變量的數(shù)學期望。例9－8設的概率密度為：，求。求解程序代碼：%連續(xù)型隨機變量的數(shù)學期望clear;clc;symsxf1=x/1500^2;f2=(3000-x)/1500^2;Ex=double(int(x*f1,0,1500)+int(x*f2,1500,3000))運行結果：Ex=1500⒊隨機變量的函數(shù)的數(shù)學期望計算公式：例9－9設圓的直徑，求圓的面積的數(shù)學期望。求解程序代碼：%連續(xù)型隨機變量的函數(shù)的數(shù)學期望clear;clc;symsxabf=1/(b-a);g=pi*x^2/4;Ey=simplify(int(f*g,x,a,b))運行結果：Ey=1/12*(a^2+b*a+b^2)*pi所以，圓的面積的數(shù)學期望為。⒋二維隨機變量的函數(shù)的數(shù)學期望計算公式：例9－10設二維隨機變量的概率密度為，求。求解程序代碼：%二維連續(xù)型隨機變量的函數(shù)的數(shù)學期望clear;clc;symsxyf=x+y;Ex=double(int(int(x*y*f,y,0,1),0,1))運行結果：Ex=0.33339.3.2設是一個隨機變量，若存在，則稱為的方差，記為。即因此，隨機變量的方差實際上是隨機變量的函數(shù)的數(shù)學期望。這里不再舉例說明。

9.3.3MATLAB提供了常見分布的均值和方差的計算函數(shù)，其調(diào)用格式如下：函數(shù)調(diào)用格式對應的分布[M,V]=betastat(a,b)分布[M,V]=binostat(n,p)二項分布[M,V]=chi2stat(v)卡方分布[M,V]=expstat(mu)指數(shù)分布[M,V]=fstat(v1,v2)F分布[M,V]=gamstat(a,b)伽馬分布[M,V]=geostat(p)幾何分布[M,V]=hygestat(M,K,n)超幾何分布[M,V]=normstat(mu,sigma)正態(tài)分布[M,V]=lognstat(mu,sigma)對數(shù)正態(tài)分布[M,V]=poisstat(lambda)泊松分布[M,V]=raylstat(b)瑞利分布[M,V]=tstat(v)學生氏t分布[M,V]=unidstat(n)離散均勻分布[M,V]=unifstat(a,b)連續(xù)均勻分布[M,V]=weibstat(a,b)威布爾分布9.3.4協(xié)方差矩陣及相關系數(shù)隨機變量與的協(xié)方差：。隨機變量與的相關系數(shù)：。設與是容量均為的二個樣本，則有：這兩個樣本的樣本均值為和；這兩個樣本的協(xié)方差為；這兩個樣本的相關系數(shù)為。相關系數(shù)常常用來衡量兩個隨機變量之間的線性相關性，相關系數(shù)的絕對值越接近1，表示相關性越強，反之越弱。隨機向量的數(shù)學期望：隨機向量的協(xié)方差矩陣：設隨機變量的觀測值均有個，構成二維數(shù)組：記，(表示樣本均值向量)，其中。則其樣本的協(xié)方差矩陣為：。用cov函數(shù)計算樣本協(xié)方差矩陣。其簡單調(diào)用格式為：●C=cov(X)：X只能是一維數(shù)組(向量)或二維數(shù)組(矩陣)?！馛=cov(X,Y)=cov([X(:),Y(:)])：僅要求X與Y的元素個數(shù)相同，即僅要求X(:)與Y(:)長度相等。若X是一個向量，則cov(X)返回樣本X的方差。若X是矩陣(每一列為一個隨機變量的觀測值)，則cov(X)返回樣本協(xié)方差矩陣。cov函數(shù)的算法為：[n,p]=size(X);Y=X-ones(n,1)*mean(X);C=Y'*Y./(n-1)用corrcoef函數(shù)計算樣本數(shù)據(jù)的相關系數(shù)矩陣。其簡單調(diào)用格式為：●R=corrcoef(X)●R=corrcoef(x,y)例9－11程序代碼：%協(xié)方差陣的計算x=[10501038;11001089;11201118;12501256;12801300];a=cov(x)運行結果：a=1.0e+004*0.99501.12001.12001.26269.3.5MATLAB中可以利用moment函數(shù)計算樣本的中心矩。該函數(shù)的調(diào)用格式為：●m=moment(X,order)：若X是向量，則moment(X,order)返回X數(shù)據(jù)的指定階次中心矩；若X是矩陣，則moment(X,order)返回X數(shù)據(jù)的每一列的指定階次中心矩。注：一階中心矩為0，二階中心矩為用除數(shù)n（而非n-1）得到的方差，其中n為向量X的長度或是矩陣X的行數(shù)。9.4樣本描述采集到大量的樣本數(shù)據(jù)以后，常常需要用一些統(tǒng)計量來描述數(shù)據(jù)的集中程度和離散程度，并通過這些指標來對數(shù)據(jù)的總體特征進行歸納。9.4.1描述樣本數(shù)據(jù)集中趨勢的統(tǒng)計量有算術平均、中位數(shù)、眾數(shù)、幾何均值、調(diào)和均值和截尾均值等。⒈幾何均值樣本數(shù)據(jù)的幾何均值：。用geomean函數(shù)計算樣本數(shù)據(jù)的幾何均值。●m=geomean(X)：若X是向量，則geomean(X)返回數(shù)據(jù)X中元素的幾何均值；若X是矩陣，則geomean(X)返回一個行向量，包含每列數(shù)據(jù)的幾何均值。⒉調(diào)和均值樣本數(shù)據(jù)的調(diào)和平均值。用harmmean函數(shù)計算樣本數(shù)據(jù)的調(diào)和平均值?！駇=harmmean(X)：若X是向量，則harmmean(X)返回數(shù)據(jù)X的調(diào)和平均值；若X是矩陣，則harmmean(X)返回包含每列元素調(diào)和平均值的行向量。⒊算術平均值樣本數(shù)據(jù)的算術平均值。用mean函數(shù)計算樣本數(shù)據(jù)的算術平均值?！駇=mean(X)：若X是向量，則mean(X)返回X的算術平均值；若X是矩陣，則mean(X)返回包含X中每列元素算術平均值的行向量。⒋中值中值是樣本數(shù)據(jù)中心趨勢的穩(wěn)健估計，因為異常值的影響較小。用median函數(shù)計算樣本數(shù)據(jù)中值?！駇=median(X)：若X是向量，則median(X)返回X的中值；若X是矩陣，則median(X)返回包含X中每列元素中值的行向量。⒌截尾均值對樣本數(shù)據(jù)進行排序以后，去掉兩端的部分極值，然后對剩下的數(shù)據(jù)求算術平均值，得到截尾均值。截尾均值為樣本位置參數(shù)的穩(wěn)健性估計。若數(shù)據(jù)中異常值，截尾均值為數(shù)據(jù)中心的一個更有代表性的估計。用函數(shù)trimmean計算截尾均值●m=trimmean(X,percent)：若樣本數(shù)據(jù)X是向量，則剔除X中最大的和最小的各0.5*percent%以后，再計算算術平均值；若樣本數(shù)據(jù)X是矩陣，則X的每列均剔除其最大的和最小的各0.5*percent%以后，再分別計算算術平均值，構成一個行向量。例9－12程序代碼：%描述樣本的集中趨勢clear;clc;x=[0123456788991011121213141516171810002000];xbar=mean(x)m1=median(x)m2=trimmean(x,20)運行結果：xbar=133.3333m1=9.5000m2=9.95009.4.2描述樣本數(shù)據(jù)離散趨勢的統(tǒng)計量包括極差、平均差、平均絕對差、方差和標準差等。⒈均值絕對差用mad函數(shù)計算數(shù)據(jù)樣本的均值或中值絕對差（MAD）?！駓=mad(X)：若X為向量，則y用mean(abs(X-mean(X)))計算；如果X為矩陣，則y為包含X中每列數(shù)據(jù)均值絕對差的行向量?！駓=mad(X,0)：與y=mad(X)相同，使用均值?！駓=mad(X,1)：y=median(abs(X-median(X)))。⒉極差極差指的是樣本中最小值與最大值之間的差值?！駓=range(X)：若X是向量，則range(X)返回X中元素的極差；若X是矩陣，則range(X)返回包含X中列元素極差的行向量。⒊方差●y=var(X)：若X是向量，則var(X)返回樣本X的方差；若X是矩陣，var(X)返回包含X中每一列元素構成的樣本的方差的行向量?！駓=var(X,1)：除以n(n是樣本大小)，是樣本數(shù)據(jù)的二階矩。●y=var(X,w)：使用權重向量w計算方差。w中元素的個數(shù)必須等于矩陣X的行數(shù)。對于向量X，w和X必須在長度上匹配。w的每一個元素必須為正。注：令SS為向量X中元素與其均值之間的偏差的平方和，則var(X)＝SS/(n-1)為的最小方差無偏估計量，var(X,1)＝SS/n為的最大似然估計量。例9－13程序代碼：%說明方差的算法clear;clc;x=[-21;-15;13;27];w=[1;2;3;4];ss0=var(x)ss1=var(x,1)ssw=var(x,w)[m,n]=size(x);xbarw=zeros(1,n);forj=1:nxbarw(j)=x(:,j)'*w./sum(w);endssw1=zeros(1,n);forj=1:nssw1(j)=((x(:,j)-xbarw(j)).^2)'*w./sum(w);endssw1運行結果：ss0=3.33336.6667ss1=2.50005.0000ssw=2.01004.3600ssw1=2.01004.3600⒋標準差●y=std(X)：若X是向量，則std(X)返回X的標準差；若X是矩陣，則std(X)返回包含X中每一列標準差的行向量。9.4.3統(tǒng)計量的分布稱為抽樣分布。正態(tài)總體的幾個常用統(tǒng)計量的分布包括分布、分布和分布。下面給出用MATLAB進行這三大分布有關的計算表。表9－1用MATLAB進行分布有關的計算表函數(shù)名類別調(diào)用格式chi2pdf密度函數(shù)chi2pdf(x,n)chi2cdf分布函數(shù)chi2cdf(x,n)chi2inv分布函數(shù)的逆函數(shù)chi2inv(P,n)chi2stat數(shù)學期望與方差[M,V]=chi2stat(n)表9－2用MATLAB進行分布有關的計算表函數(shù)名類別調(diào)用格式tpdf密度函數(shù)tpdf(x,n)tcdf分布函數(shù)tcdf(x,n)tinv分布函數(shù)的逆函數(shù)tinv(P,n)tstat數(shù)學期望與方差[M,V]=tstat(n)表9－3用MATLAB進行分布有關的計算表函數(shù)名類別調(diào)用格式fpdf密度函數(shù)fpdf(x,n1,n2)fcdf分布函數(shù)fcdf(x,n1,n2)finv分布函數(shù)的逆函數(shù)finv(P,n1,n2)fstat數(shù)學期望與方差[M,V]=fstat(n1,n2)9.5參數(shù)估計參數(shù)估計的內(nèi)容包括點估計和區(qū)間估計。MATLAB統(tǒng)計工具箱提供了進行最大似然估計的函數(shù)，可以計算待估參數(shù)及其置信區(qū)間。利用專門的參數(shù)估計函數(shù)，可以估計服從不同分布的函數(shù)的參數(shù)。9.5.1點估計是用單個數(shù)值作為參數(shù)的估計，常用的方法有矩法和最大似然法。⒈矩法某些情況下，待估參數(shù)往往是總體原點矩或原點矩的函數(shù)，此時可以用取自該總體的樣本的原點矩或樣本原點矩的函數(shù)值作為待估參數(shù)的估計，這種方法稱為矩法。例如，樣本均值總是總體均值的矩估計量，樣本方差總是總體方差的矩估計量，樣本標準差總是總體標準差的矩估計量。用計算矩的函數(shù)moment(X,order)進行估計。例9－14對某型號的20輛汽車記錄其5L汽油的行駛里程（公里），觀測數(shù)據(jù)如下：29.827.628.327.930.128.729.928.027.928.728.427.229.528.528.030.029.129.829.626.9試求總體的均值和方差的矩法估計。求解程序代碼：%矩法估計clear;clc;x1=[29.827.628.327.930.128.729.928.027.928.7];x2=[28.427.229.528.528.030.029.129.829.626.9];x=[x1x2]';muhat=mean(x)sigma2hat=moment(x,2)%樣本二階中心矩,也可用var(x,1)。運行結果：muhat=28.6950sigma2hat=0.9185⒉最大似然法最大似然法是在待估參數(shù)的可能取值范圍內(nèi)進行挑選，使似然函數(shù)值（即樣本取固定觀察值鄰域內(nèi)的概率）最大的那個參數(shù)值即為最大似然估計量。由于最大似然估計法得到的估計量通常不僅僅滿足無偏性、有效性等基本條件，還能保證其為充分統(tǒng)計量，所以，在點估計和區(qū)間估計中，一般推薦使用最大似然法。用函數(shù)mle進行最大似然估計。該函數(shù)的常見調(diào)用格式為：phat=mle('dist',data)使用向量data中的樣本數(shù)據(jù)，返回dist指定的分布的最大似然估計(MLE)。'dist'常見取值：'beta'(Beta分布)，'bernoulli'(0-1分布)，'binomial'(二項分布)，'exponential'(指數(shù)分布)，'gamma'(GAMMA分布)，'geometric'(幾何分布)，'lognormal'(對數(shù)正態(tài)分布)，'normal'(正態(tài)分布)，'negativebinomial'(負二項分布)，'poisson'(泊松分布)，'rayleigh'(瑞利分布)，'discreteuniform'(離散均勻分布)，'uniform'(均勻分布)，'weibull'(Weibull分布)。'dist'取'normal'時可省略。例9－15對某型號的20輛汽車記錄其5L汽油的行駛里程（公里），觀測數(shù)據(jù)如下：29.827.628.327.930.128.729.928.027.928.728.427.229.528.528.030.029.129.829.626.9設行駛里程服從正態(tài)分布，試用最大似然估計法估計總體的均值和方差。求解程序代碼：%最大似然估計clear;clc;x1=[29.827.628.327.930.128.729.928.027.928.7];x2=[28.427.229.528.528.030.029.129.829.626.9];x=[x1x2];p=mle('norm',x);muhatmle=p(1)sigma2hatmle=p(2)^2運行結果：muhatmle=28.6950sigma2hatmle=0.91859.5.2求參數(shù)的區(qū)間估計，首先要求出該參數(shù)的點估計，然后構造一個含有該參數(shù)的隨機變量，并根據(jù)一定的置信水平求該估計值的范圍。函數(shù)mle除了可以用于求指定分布的參數(shù)的最大似然估計外，還可以用于求指定分布的參數(shù)的區(qū)間估計。其簡單的調(diào)用格式為：●[phat,pci]=mle('dist',data)：返回'dist'分布的參數(shù)的最大似然估計和置信度為95%的置信區(qū)間?！馵phat,pci]=mle('dist',data,alpha)：返回'dist'分布的參數(shù)的最大似然估計和置信度為100(1-alpha)%的置信區(qū)間。例9－16對某型號的20輛汽車記錄其5L汽油的行駛里程（公里），觀測數(shù)據(jù)如下：29.827.628.327.930.128.729.928.027.928.728.427.229.528.528.030.029.129.829.626.9設行駛里程服從正態(tài)分布，求平均行駛里程的95%置信區(qū)間。求解程序代碼：%正態(tài)總體參數(shù)的區(qū)間估計clear;clc;x1=[29.827.628.327.930.128.729.928.027.928.7];x2=[28.427.229.528.528.030.029.129.829.626.9];x=[x1x2]';[p,ci]=mle('norm',x,0.05)運行結果：p=28.69500.9584ci=28.23480.747829.15521.4361即平均行駛里程的95%置信區(qū)間為（28.2348，291552）。另外，0.9584，的95%置信區(qū)間(0.7478，1.4361)。9.5.3常見分布的參數(shù)估計除了使用mle函數(shù)求指定分布的參數(shù)估計量外，MATLAB統(tǒng)計工具箱還提供了求常見分布的參數(shù)的估計函數(shù)，如表9－4所示。表9－4常見分布的參數(shù)估計函數(shù)及其簡單調(diào)用格式分布調(diào)用格式貝塔分布phat=betafit(x)[phat,pci]=betafit(x,alpha)二項分布phat=binofit(x,n)[phat,pci]=binofit(x,n)[phat,pci]=binofit(x,n,alpha)指數(shù)分布muhat=expfit(x)[muhat,muci]=expfit(x)[muhat,muci]=expfit(x,alpha)伽馬分布phat=gamfit(x)[phat,pci]=gamfit(x)[phat,pci]=gamfit(x,alpha)正態(tài)分布[muhat,sigmahat,muci,sigmaci]=normfit(x)[muhat,sigmahat,muci,sigmaci]=normfit(x,alpha)泊松分布lambdahat=poissfit(x)[lambdahat,lambdaci]=poissfit(x)[lambdahat,lambdaci]=poissfit(x,alpha)均勻分布[ahat,bhat]=unifit(x)[ahat,bhat,aci,bci]=unifit(x)[ahat,bhat,aci,bci]=unifit(x,alpha)威布爾分布phat=weibfit(x)[phat,pci]=weibfit(x)[phat,pci]=weibfit(x,alpha)例如，用normfit函數(shù)對正態(tài)分布總體進行參數(shù)估計?！馵muhat,sigmahat,muci,sigmaci]=normfit(x)：對于給定的正態(tài)分布的數(shù)據(jù)x，返回參數(shù)的估計值muhat、的估計值sigmahat、的95%置信區(qū)間muci和的95%置信區(qū)間sigmaci。●[muhat,sigmahat,muci,sigmaci]=normfit(x,alpha)：進行參數(shù)估計并計算100(1-alpha)%置信區(qū)間。例9－17用normfit函數(shù)求解例9－16。解：由可得的置信區(qū)間：由可得的置信區(qū)間：求解程序代碼：%normfit函數(shù)應用舉例clear;clc;a=0.05;x1=[29.827.628.327.930.128.729.928.027.928.7];x2=[28.427.229.528.528.030.029.129.829.626.9];x=[x1x2]';[muhat,sigmahat,muci,sigmaci]=normfit(x,a)n=numel(x);muci1=[muhat-tinv(1-a/2,n-1)*sigmahat/sqrt(n),...muhat+tinv(1-a/2,n-1)*sigmahat/sqrt(n)]sigmaci1=[((n-1).*sigmahat.^2/chi2inv(1-a/2,n-1)).^0.5,...((n-1).*sigmahat.^2/chi2inv(a/2,n-1)).^0.5]運行結果：muhat=28.6950sigmahat=0.9833muci=28.234829.1552sigmaci=0.74781.4361muci1=28.234829.1552sigmaci1=0.74781.4361請讀者比較例9－16(利用mle函數(shù))和例9－17(利用normfit函數(shù))的結果。9.6假設檢驗在總體分布函數(shù)完全未知或只知分布形式，但不知其參數(shù)時，為了推斷總體的某些性質(zhì)，需要提出關于總體的假設。假設是否合理，需要檢驗。9.6.1在MATLAB中，對于方差已知的正態(tài)總體，關于均值的檢驗用ztest函數(shù)。其調(diào)用格式為：●h=ztest(x,m,sigma,alpha)：在顯著性水平alpha下進行z檢驗，以檢驗服從正態(tài)分布的樣本x是否來自均值為m的正態(tài)總體。sigma為標準差。若返回結果h=1，則可以在顯著性水平alpha下接受備擇假設(拒絕：)；若返回結果h=0，則在顯著性水平alpha下不能拒絕。在alpha為0.05時，可省略。●[h,sig,ci,zval]=ztest(x,m,sigma,alpha,tail)：tail用于指定是進行單側檢驗還是進行雙側檢驗。tail參數(shù)可以有下面幾個取值：tail=0或'both'(為默認設置)指定備擇假設為均值不等于m；tail=1或'right'指定備擇假設為均值大于m；tail=-1或'left'指定備擇假設為均值小于m。zval是標準正態(tài)分布統(tǒng)計量的值。sig為與統(tǒng)計量有關的p值。sig為能夠由統(tǒng)計量的值zval做出拒絕原假設的最小顯著性水平。即若tail=0，則sig=P{|U|>zval}；若tail=1，則sig=P{U>zval}；若tail=-1，則sig=P{U<zval}。ci為均值真值的1-alpha置信區(qū)間。例9－18下面列出的是某工廠隨機選取的20只零部件的裝配時間（分）：9.810.410.69.69.79.910.911.19.610.210.39.69.911.210.69.810.510.110.59.7設裝配時間的總體服從正態(tài)分布，標準差為0.4，是否可以認為裝配時間的均值在0.05的水平下不小于10分鐘。解：：：程序代碼：%正態(tài)總體的方差已知時的均值檢驗clear;clc;x1=[9.810.410.69.69.79.910.911.19.610.2];x2=[10.39.69.911.210.69.810.510.110.59.7];x=[x1x2]';m=10;sigma=0.4;a=0.05;[h,sig,muci]=ztest(x,m,sigma,a,1)運行結果：h=1sig=0.0127muci=10.0529Inf因此，在0.05的水平下，可以認為裝配時間的均值不小于10分鐘。9.6.2方差未知時在MATLAB中，對于均方差未知的正態(tài)總體，關于均值的檢驗用ttest函數(shù)。其調(diào)用格式為：●h=ttest(x,m)：在顯著性水平0.05下進行t檢驗，以檢驗服從正態(tài)分布(標準差未知)的樣本x是否來自均值為m的正態(tài)總體。(說明：當m=0時，可省略，即ttest(x)=ttest(x,0))●h=ttest(x,m,alpha)：在顯著性水平alpha下進行t檢驗，以檢驗服從正態(tài)分布(標準差未知)的樣本x是否來自均值為m的正態(tài)總體。若返回結果h=1，則可以在顯著性水平alpha下接受備擇假設(拒絕：)；若返回結果h=0，則在顯著性水平alpha下不能拒絕?！馵h,sig,ci,stats]=ttest(x,m,alpha,tail)：tail用于指定是進行單側檢驗還是進行雙側檢驗。tail參數(shù)可以有下面幾個取值：tail=0或'both'(為默認設置)指定備擇假設為均值不等于m；tail=1或'right'指定備擇假設為均值大于m；tail=-1或'left'指定備擇假設為均值小于m。sig為與樣本x有關的p值。sig為能夠由樣本x做出拒絕原假設的最小顯著性水平。ci為均值真值的1-alpha置信區(qū)間。結構數(shù)組stats中包含統(tǒng)計量的值、自由度和樣本標準差。例9－19某種電子元件的壽命(以小時計)服從正態(tài)分布，和均未知。現(xiàn)測得16只元件的壽命如下：159280101212224379179264222362168250149260485170問是否有理由認為元件的平均壽命大于225(小時)？解：：：程序代碼：%正態(tài)總體的方差求知時的均值檢驗clear,clc,x=[159280101212224379179264222362168250149260485170];m=225;a=0.05;[h,sig,muci,stats]=ttest(x,m,a,1)運行結果：h=0sig=0.2570muci=198.2321Infstats=tstat:0.6685df:15sd:98.7259由于sig=0.257，因此沒有充分的理由認為元件的平均壽命大于225小時。而對于：：程序代碼：%正態(tài)總體的方差求知時的均值檢驗clear;clc;x=[159280101212224379179264222362168250149260485170];m=225;a=0.05;[h,sig]=ttest(x,m,a,-1)運行結果：h=0sig=0.7430由于sig=0.743，因此更沒有充分的理由認為元件的平均壽命小于225小時。9.6.3兩個正態(tài)總體(方差未知但相等)在MATLAB中，對于兩個獨立正態(tài)總體(方差均未知但相等)，關于其均值的檢驗用ttest2函數(shù)。其簡單調(diào)用格式為：●h=ttest2(x,y)：x和y為取自兩個獨立正態(tài)分布(方差未知但相等)的兩個樣本，檢驗兩個正態(tài)總體的均值是否相等。若返回h=1時，則可以在0.05的水平下拒絕(均值相等)，即可認為兩個總體的均值不等；若返回h=0時，則不能在0.05水平下拒絕，即不能認為兩個總體均值不等。●[h,significance,ci]=ttest2(x,y,alpha,tail)：alpha為給定的顯著性水平，tail用于指定是進行單側檢驗還是雙側檢驗。若返回h=1時，則可以在alpha的水平下拒絕；若返回h=0時，則不能在alpha水平下拒絕。significance是與x,y有關的p值。即為能夠利用x,y做出拒絕的最小顯著性水平。ci為兩個總體均值差()的1-alpha置信區(qū)間。tail可以有下面3個取值：tail=0或'both'(為默認設置)指定備擇假設。tail=1或'right'指定備擇假設。tail=-1或'left'指定備擇假設。例9－20某廠鑄造車間為提高鑄件的耐磨性而試制了一種鎳合金鑄件以取代銅合金鑄件，為此，從兩種鑄件中各獨立地抽取一個容量分別為8和9的樣本，測得其硬度（一種耐磨性指標）為：鎳合金76.4376.2173.5869.6965.2970.8382.7572.34銅合金73.6664.2769.3471.3769.7768.1267.2768.0762.61根據(jù)專業(yè)經(jīng)驗，硬度服從正態(tài)分布，且方差保持不變，試在顯著性水平下判斷鎳合金的硬度是否有明顯提高。解：：：程序代碼：%正態(tài)總體的方差齊但求知時的均值差檢驗clear;clc;x=[76.4376.2173.5869.6965.2970.8382.7572.34]';y=[73.6664.2769.3471.3769.7768.1267.2768.0762.61]';a=0.05;[h,sig,ci]=ttest2(x,y,a,1)運行結果：h=1sig=0.0142ci=1.4148Inf因此，在顯著性水平下，可以認為鎳合金的硬度有明顯提高。在實際工作中，為了比較兩種方法或兩種產(chǎn)品的差異，常常進行對比試驗。這樣得到的數(shù)據(jù)具有成對的特點，即上述的x與y長度相同。在MATLAB中，分析這種數(shù)據(jù)，還可以用ttest函數(shù)進行檢驗。下面舉例說明。例9－21在不同蒸汽壓下保持了8小時的紅花苜蓿半穗中花密的含糖濃度數(shù)據(jù)如下：4.4mmHg62.565.267.669.969.470.167.867.068.562.49.9mmHg51.754.253.357.056.461.557.256.258.455.8根據(jù)經(jīng)驗，濃度服從正態(tài)分布，且方差保持不變，檢驗不同蒸汽壓對紅花苜蓿半穗中花密的含糖濃度是否有顯著影響。解：：：程序代碼：%基于成對數(shù)據(jù)的檢驗clear;clc;x=[62.565.267.669.969.470.167.867.068.562.4]';y=[51.754.253.357.056.461.557.256.258.455.8]';[h,sig]=ttest(x,y)運行結果：h=1sig=8.6831e-008由于sig=，因此有充分理由認為不同蒸汽壓對紅花苜蓿半穗中花密的含糖濃度有顯著影響。9.6.4分布擬合檢驗是統(tǒng)計分析中常常用到的方法，檢驗方法比較多，如圖示法、峰度-偏度法、以及一些非參數(shù)法等。下面介紹兩種比較簡單的方法，即q-q圖法和峰度-偏度法。⒈q-q圖q-q圖用變量數(shù)據(jù)分布的分位數(shù)與所指定分布的分位數(shù)之間的關系曲線來檢驗數(shù)據(jù)的分布。如果兩個樣本來自同一分布，則圖中數(shù)據(jù)點呈現(xiàn)直線關系，否則為曲線關系。用qqplot函數(shù)生成兩個樣本的q-q圖，其簡單調(diào)用格式如下：●qqplot(X)：若X是向量，則顯示X的樣本值與服從正態(tài)分布的理論數(shù)據(jù)之間的q-q圖(若X的分布為正態(tài)分布，則圖形接近直線。)；若X是矩陣，則顯示X的每列數(shù)據(jù)與服從正態(tài)分布的理論數(shù)據(jù)之間的q-q圖。●qqplot(X,Y)：若X和Y均是向量，則顯示兩個樣本的q-q圖(若樣本是來自于相同的分布，則圖形將是線性的。)；若X和Y不全是向量(X或Y是矩陣)時，q-q圖為一個配對列顯示分隔線。例9－22有關q-q圖的程序代碼：%分布擬合檢驗(q-q圖法)x=normrnd(0,1,100,1);%生成服從正態(tài)分布的隨機數(shù)y=normrnd(0.5,2,50,1);z=weibrnd(2,0.5,100,1);%生成服從威布爾分布的隨機數(shù)subplot(2,2,1)%說明生成子圖的位置qqplot(x)holdonsubplot(2,2,2)qqplot(x,y)holdonsubplot(2,2,3)qqplot(z)holdonsubplot(2,2,4)qqplot(x,z)holdoff運行結果見圖9－3。圖9－3q-q圖上面的q-q圖中第1個子圖用x的數(shù)據(jù)繪圖，因為服從正態(tài)分布，圖中數(shù)據(jù)點呈直線分布；第2個子圖用x數(shù)據(jù)和y數(shù)據(jù)（均服從正態(tài)分布），數(shù)據(jù)點的主體部分呈直線；第3個子圖用z數(shù)據(jù)繪圖，由于它服從威布爾分布，所以數(shù)據(jù)點不在一條直線上；第4個子圖是用x數(shù)據(jù)和z數(shù)據(jù)繪制的，因為它們不是同分布的，圖中數(shù)據(jù)點不呈直線分布。⒉峰度-偏度檢驗樣本階矩：樣本階中心矩：樣本偏度：樣本偏差反映了總體分布密度曲線的對稱性信息。如果數(shù)據(jù)完全對稱，則。樣本峰度：（茆詩松書中的定義為）峰度-偏度檢驗又稱為Jarque-Bera檢驗，評價給定數(shù)據(jù)服從未知均值和方差的正態(tài)分布的假設是否成立。該檢驗基于數(shù)據(jù)樣本的偏度和峰度。對于正態(tài)分布數(shù)據(jù)，樣本偏度接近于0，樣本峰度接近于3。Jarque-Bera檢驗確定樣本偏度和峰度是否與它們的期望值相差較遠。注意：Jarque-Bera檢驗不能用于小樣本的檢驗。在MATLAB中，用jbtest函數(shù)進行Jarque-Bera檢驗，測試數(shù)據(jù)對正態(tài)分布的似合程度，其調(diào)用格式為：●h=jbtest(X)：對輸入數(shù)據(jù)向量X進行Jarque-Bera檢驗，返回檢驗結果h。若h=1，則在顯著性水平0.05下拒絕X服從正態(tài)分布的假設；若h=0，我們可認為X服從正態(tài)分布。●h=jbtest(X,alpha)：在顯著性水平alpha下進行Jarque-Bera檢驗?！馵h,P,JBSTAT,CV]=jbtest(X,alpha)：還返回3個其他輸出。P為檢驗的值，JBSTAT為檢驗統(tǒng)計量，CV為確定是否拒絕零假設的臨界值。例9－23某氣象站收集了44個獨立的年降雨量數(shù)據(jù)，資料如下（已排序）：52055656161663566968669270470771171371471972773574074474575077677778678679179482182282683483785186287387988990090492292695296310561074試檢驗這些數(shù)據(jù)是否來自正態(tài)總體。求解程序代碼：%峰度－偏度檢驗clear;clcx1=[520556561616635669686692704707711];x2=[713714719727735740744745750776777];x3=[786786791794821822826834837851862];x4=[87387988990090492292695296310561074];x=[x3x2x1x4];[H,P,JBSTAT,CV]=jbtest(x)運行結果：H=0P=0.9356JBSTAT=0.1332CV=5.9915由于P=0.9356，因此有充分理由認為上述數(shù)據(jù)是來自正態(tài)總體。⒊秩和檢驗秩和檢驗是用來檢驗兩個總體的分布函數(shù)是否相等的，在MATLAB中用ranksum函數(shù)實現(xiàn)秩和檢驗。其調(diào)用格式為：●p=ranksum(x,y,alpha)：進行雙側秩和檢驗，零假設為向量x和y表示的兩個獨立樣本取自相同的分布。返回的p值表示由樣本x與y能拒絕零假設的最小顯著性水平?！馵p,h]=ranksum(x,y)：返回假設檢驗的結果到h中，若h為0，則在顯著性水平0.05下，可以認為x和y取自相同的分布；若h為1，則在顯著性水平0.05下，可以認為x和y取自不同的分布?！馵p,h]=ranksum(x,y,’alpha’,alpha)：在顯著性水平alpha下做檢驗。●[p,h]=ranksum(┅,’method’,method)：設置計算p值的算法。Method參數(shù)設置為’exact’時，使用精確算法；設置為’approximate’時，使用近似算法。如果忽略該參數(shù)，則對小樣本使用精確算法，對大樣本使用近似算法?！馵p,h,ststs]=ranksum(┅)：返回一個有一個或兩個字段的結構?！痳anksum’包含秩和統(tǒng)計量的值。如果樣本比較大，則p用近似算法計算，’zval’包含Z統(tǒng)計量的值。例9－24某醫(yī)院外科用兩種手術方法治療情況基本相同的肝癌患者11例?；颊卟捎秒S機方法分配到不同手術組。每例手術后生存月數(shù)在下面，試比較兩種手術方法的術后生存月數(shù)有否差別？甲法235612乙法37881011求解程序代碼：%秩和檢驗clear;clcx=[235612];y=[37881011];[p,h,ststs]=ranksum(x,y)運行結果：p=0.2727h=0ststs=ranksum:23.5000由于p=0.2727，因此可以認為兩種手術方法的術后生存月數(shù)沒有顯著差別。9.7方差分析事件的發(fā)生往往與多個因素有關，但各個因素對事件發(fā)生的影響可能是不一樣的，而且同一因素的不同水平對事件發(fā)生的影響也可能是不同的。通過方差分析，便可以研究不同因素以及因素的不同水平對事件發(fā)生的影響程度。根據(jù)自變量個數(shù)的不同，方差分析可以分為單因子方差分析和多因子方差分析。9.7.1用anova1函數(shù)進行單因子方差分析?！駊=anova1(X)：的第列是因子的第個水平的m個相互獨立觀測值，比較樣本數(shù)據(jù)的各列數(shù)據(jù)的均值。零假設為各水平的均值相等。若返回的p值接近0，則可拒絕，即可認為因子是“統(tǒng)計上顯著”。為了確定結果是否“統(tǒng)計上顯著”，需要確定p值。該值由自己確定。一般地，當p值小于0.05或0.01時，認為因子是顯著的。anova1函數(shù)還生成兩個圖表。第1個圖表為方差分析表，方差分析表中有6列：第1列顯示誤差的來源；第2列顯示每一個誤差來源的平方和（SS）；第3列顯示與每一個誤差來源相關的自由度（df）；第4列顯示均方和（MS），它是誤差來源平方和與自由度的比值，即SS/df；第5列顯示F統(tǒng)計量，它是均方和的比值；第6列顯示p值，p值是能夠拒絕零假設(為各水平的均值相等)的最小顯著性水平。第2個圖顯示X的每一列的箱形圖?！馻nova1(X,group)：若X是m行n列矩陣，則X的第列是因子的第個水平的m個相互獨立觀測值，group（字符型單元數(shù)組）的每個元素作為X中對應列中的數(shù)據(jù)的標簽，所以變量的長度必須等于X的列數(shù)。若X為向量，group用來標識向量X中的每個元素的水平，因此，group與X的長度必須相等。group中包含的標簽同樣用于箱形圖的標注。注：anova1(X,group)(X為向量)形式不需要每個水平的觀測值個數(shù)相同，所以它適用于不平衡(不等重復)單因子試驗?！駊=anova1(X,group,‘displayopt’)：當‘displayopt’參數(shù)設置為‘on’(默認設置)時，激活ANOVA表和箱形圖的顯示；‘displayopt’參數(shù)設置為‘off’時，不予顯示?！馵p,table]=anova1(…)：還返回ANOVA表（包含列標簽和行標簽）?！馵p,table,stats]=anova1(…)：還返回stats結構，可用于進行多重比較檢驗。當檢驗結果為因子顯著時，有時需要作進一步檢驗，確定哪對均值差異顯著，哪對均值差異不顯著。提供stats結構作為輸入，使用multcompare函數(shù)便可以進行此項檢驗。注：方差分析要求樣本數(shù)據(jù)滿足下面的假設條件：⑴所有樣本數(shù)據(jù)滿足正態(tài)分布條件；⑵所有樣本數(shù)據(jù)具有相等的方差；⑶所有觀測值相互獨立。在基本滿足前二個假設條件的情況下，一般認為ANOVA檢驗是穩(wěn)健的。例9－25某化工產(chǎn)品的產(chǎn)量是衡量經(jīng)濟效益的重要指標。為了考察反應溫度（因子）對該化工產(chǎn)品產(chǎn)量是否有顯著影響，我們選取因子的5個水平為：60℃，：65℃，：70℃，：75℃，：80℃。每個水平下各作5次重復試驗，試驗結果如下表所示：反映溫度觀察值（㎏）606570758090879291889791939592969293969584828683888481858682解：設該試驗的線性統(tǒng)計模型為：作方差分析的程序代碼：%等重復的單因子試驗的方差分析clear;clcy1=[9087929188]';y2=[9791939592]';y3=[9692939695]';y4=[8482868388]';y5=[8481858682]';y=[y1y2y3y4y5];p=anova1(y,{'60','65','70','75','80'})%y的列表示重復觀測值。運行結果見圖9－4和圖9－5。圖9－4方差分析表圖9－5箱形圖從方差分析表可見，反應溫度（因子）對該化工產(chǎn)品產(chǎn)量有顯著影響。

例9－26研究6種農(nóng)藥對殺蟲效果的影響，試驗所得數(shù)據(jù)如下表：農(nóng)藥編號殺蟲量12345687.485.080.290.588.587.394.356.262.455.048.292.099.295.391.575.272.381.3解：設該試驗的線性統(tǒng)計模型為：其中，。作方差分析的程序代碼：%不等重復的單因子試驗的方差分析clear;clcy1=[87.485.080.2];y2=[90.588.587.394.3];y3=[56.262.4];y4=[55.048.2];y5=[92.099.295.391.5];y6=[75.272.381.3];y=[y1y2y3y4y5y6];A1=ones(numel(y1),1);A2=2*ones(numel(y2),1);A3=3*ones(numel(y3),1);A4=4*ones(numel(y4),1);A5=5*ones(numel(y5),1);A6=6*ones(numel(y6),1);A=[A1;A2;A3;A4;A5;A6];[p,table,stats]=anova1(y,A)運行結果的方差分析表如下：可見，這6種農(nóng)藥對殺蟲效果有顯著的影響。運行結果的stats結構如下：stats=gnames:{6x1cell}n:[342243]source:'anova1'means:[84.200090.150059.300051.600094.500076.2667]df:12s:3.8470從means:[84.200090.150059.300051.600094.500076.2667]可以看出，第5種農(nóng)藥的殺蟲效果最好。作多重比較的命令：B=multcompare(stats,0.05)運行結果：B=1.00002.0000-15.8193-5.95003.91931.00003.000013.104024.900036.69601.00004.000020.804032.600044.39601.00005.0000-20.1693-10.3000-0.43071.00006.0000-2.61747.933318.48402.00003.000019.659330.850042.04072.00004.000027.359338.550049.74072.00005.0000-13.4872-4.35004.78722.00006.00004.014113.883323.75263.00004.0000-5.22197.700020.62193.00005.0000-46.3907-35.2000-24.00933.00006.0000-28.7627-16.9667-5.17064.00005.0000-54.0907-42.9000-31.70934.00006.0000-36.4627-24.6667-12.87065.00006.00008.364118.233328.1026在multcompare函數(shù)的輸出B中，第1、2兩列標識水平，第4列表示均值差的估計（的估計），第3列和第5列表示的置信區(qū)間，若置信區(qū)間中包含0，則說明在顯著性水平（這里為0.05）下，相應的與沒有顯著差異；若置信區(qū)間中不包含0，則說明在顯著性水平（這里為0.05）下，相應的與有顯著差異。在本例中，第5種農(nóng)藥與第1、3、4、6四種農(nóng)藥有顯著差異，而第5種農(nóng)藥與第2種農(nóng)藥沒有顯著差異。

9.7.2用anova2函數(shù)進行平衡雙因子方差分析。●p=anova2(X,reps)：進行平衡(等重復reps次)雙因子方差分析，的不同列中的數(shù)據(jù)代表列因子A(個水平)的變化。不同行中的數(shù)據(jù)代表行因子B(個水平)的變化。其中。當reps=1(默認值可省略)時，anova2函數(shù)返回兩個p值到p向量中：◣零假設H0A的p值。零假設H0A為列因子A的所有樣本（即X中的所有列樣本）取自相同的總體。◣零假設H0B的p值。零假設H0B為行因子B的所有樣本（如X中的所有行樣本）取自相同的總體。當reps＞1時，anova2在向量中還返回第3個值：◣零假設H0AB的p值。零假設H0AB為因子A和因子B之間沒有交互效應。如果任意一個p值接近于0，則認為相應的零假設不成立。對于零假設H0A，一個足夠小的p值表示至少有一個列樣本均值明顯地不同于其他列樣本均值，即列因子A存在主效應；對于零假設H0B，一個足夠小的p值表示至少有一個行樣本均值明顯地不同于其他行樣本均值，即行因子B存在主效應；對于零假設H0AB，一個足夠小的p值表示因子A與因子B之間存在交互效應。為了決定結果是否是“統(tǒng)計上顯著的”，需要確定p值。一般地，當p值小于0.05或0.01時，認為結果是顯著的。anova2函數(shù)還顯示一個含方差分析表的圖形。該方差分析表中包含6列：第1列顯示誤差的來源；第2列顯示來源于每一個誤差來源的平方和（SS）；第3列為與每一個誤差來源相關的自由度（df）；第4列為均方和（MS），它是誤差平方和與自由度的比值，即SS/df；第5列為F統(tǒng)計量，它是均方和的比值；第6列為p值?！駊=anova2(X,reps,‘displayopt’)：進行平衡(等重復reps次)雙因子方差分析。當‘displayopt’為‘on’（默認設置）時，激活ANOVA表的顯示；‘displayopt’為‘off’時，不予顯示?！馵p,table]=anova2(…)：還返回ANOVA表（包含列標簽和行標簽）?！馵p,table,stats]=anova2(…)：返回stats結構，用于進行列因子均值的多重比較檢驗。例9－27為了考察高溫合金中碳的含量（因子）和銻與鋁的含量之和（因子）對合金強度的影響。因子取3個水平0.03、0.04、0.05（上述數(shù)字表示碳的含量占合金總量的百分比），因子取4個水平3.3、3.4、3.5、3.6（上述數(shù)字的意義同上）。在每個水平組合下各作一次試驗，試驗結果如下：（銻與鋁的含量之和）3.33.43.53.6（碳的含量）0.030.040.0563.163.965.666.865.166.467.869.067.271.071.973.5設上表中的數(shù)據(jù)符合，請作方差分析。解：作方差分析程序代碼：%平衡雙因子方差分析(不重復)clear;clcy=[63.163.965.666.8;65.166.467.869.0;67.271.071.973.5];p=anova2(y)運行結果：可見，列因子B和行因子A對合金強度的影響均是顯著的。例9－28為了考察某種電池的最大輸出電壓受板極材料與使用電池的環(huán)境溫度的影響，材料類型(因子)取3個水平(即3種不同的材料)，溫度(因子)也取3個水平，每個水平組合下重復4次試驗，數(shù)據(jù)如下：溫度152535材料類型113015517418034408075207082582150188159126136122106115257058453138110168160174120150139961048260解作方差分析程序代碼(材料類型作為列因子)：%兩因子等重復試驗（方差分析）clear;clcy11=[130155174180];y12=[34408075];y13=[20708258];y1=[y11y12y13]';y21=[150188159126];y22=[136122106115];y23=[25705845];y2=[y21y22y23]';y31=[138110168160];y32=[174120150139];y33=[961048260];y3=[y31y32y33]';y=[y1y2y3];[p,table,stats]=anova2(y,4)運行結果(方差分析表)：可見，列因子A、行因子B以及其交互作用均顯著。運行結果(stats的結構)：stats=source:'anova2'sigmasq:502.9907colmeans:[91.5000108.3333125.0833]coln:12rowmeans:[153.1667107.583364.1667]rown:12inter:1pval:8.0678e-004df:27◣再作列因子的多重比較，命令如下：c=multcompare(stats,0.05)運行結果：c=1.00002.0000-39.5348-16.83335.86811.00003.0000-56.2848-33.5833-10.88192.00003.0000-39.4515-16.75005.9515對于不平衡雙因子方差分析和多因子方差分析，在MATLAB中用anovan函數(shù)。anov