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

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

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

9.2.2若連續(xù)型隨機(jī)變量的密度函數(shù)為，則關(guān)于和的邊緣概率密度和分別為：例9－5設(shè)具有概率密度⑴確定常數(shù)；⑵求邊緣概率密度和。解：⑴由可得；計(jì)算程序代碼：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))運(yùn)行結(jié)果：fx=21/8*x^2*(1-x^4)fy=7/2*y^(5/2)因此，，

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

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

例9－26研究6種農(nóng)藥對殺蟲效果的影響，試驗(yàn)所得數(shù)據(jù)如下表：農(nóng)藥編號殺蟲量12345687.485.080.290.588.587.394.356.262.455.048.292.099.295.391.575.272.381.3解：設(shè)該試驗(yàn)的線性統(tǒng)計(jì)模型為：其中，。作方差分析的程序代碼：%不等重復(fù)的單因子試驗(yàn)的方差分析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)運(yùn)行結(jié)果的方差分析表如下：可見，這6種農(nóng)藥對殺蟲效果有顯著的影響。運(yùn)行結(jié)果的stats結(jié)構(gòu)如下：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)運(yùn)行結(jié)果：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兩列標(biāo)識水平，第4列表示均值差的估計(jì)（的估計(jì)），第3列和第5列表示的置信區(qū)間，若置信區(qū)間中包含0，則說明在顯著性水平（這里為0.05）下，相應(yīng)的與沒有顯著差異；若置信區(qū)間中不包含0，則說明在顯著性水平（這里為0.05）下，相應(yīng)的與有顯著差異。在本例中，第5種農(nóng)藥與第1、3、4、6四種農(nóng)藥有顯著差異，而第5種農(nóng)藥與第2種農(nóng)藥沒有顯著差異。

9.7.2用anova2函數(shù)進(jìn)行平衡雙因子方差分析?！駊=anova2(X,reps)：進(jìn)行平衡(等重復(fù)reps次)雙因子方差分析，的不同列中的數(shù)據(jù)代表列因子A(個(gè)水平)的變化。不同行中的數(shù)據(jù)代表行因子B(個(gè)水平)的變化。其中。當(dāng)reps=1(默認(rèn)值可省略)時(shí)，anova2函數(shù)返回兩個(gè)p值到p向量中：◣零假設(shè)H0A的p值。零假設(shè)H0A為列因子A的所有樣本（即X中的所有列樣本）取自相同的總體。◣零假設(shè)H0B的p值。零假設(shè)H0B為行因子B的所有樣本（如X中的所有行樣本）取自相同的總體。當(dāng)reps＞1時(shí)，anova2在向量中還返回第3個(gè)值：◣零假設(shè)H0AB的p值。零假設(shè)H0AB為因子A和因子B之間沒有交互效應(yīng)。如果任意一個(gè)p值接近于0，則認(rèn)為相應(yīng)的零假設(shè)不成立。對于零假設(shè)H0A，一個(gè)足夠小的p值表示至少有一個(gè)列樣本均值明顯地不同于其他列樣本均值，即列因子A存在主效應(yīng)；對于零假設(shè)H0B，一個(gè)足夠小的p值表示至少有一個(gè)行樣本均值明顯地不同于其他行樣本均值，即行因子B存在主效應(yīng)；對于零假設(shè)H0AB，一個(gè)足夠小的p值表示因子A與因子B之間存在交互效應(yīng)。為了決定結(jié)果是否是“統(tǒng)計(jì)上顯著的”，需要確定p值。一般地，當(dāng)p值小于0.05或0.01時(shí)，認(rèn)為結(jié)果是顯著的。anova2函數(shù)還顯示一個(gè)含方差分析表的圖形。該方差分析表中包含6列：第1列顯示誤差的來源；第2列顯示來源于每一個(gè)誤差來源的平方和（SS）；第3列為與每一個(gè)誤差來源相關(guān)的自由度（df）；第4列為均方和（MS），它是誤差平方和與自由度的比值，即SS/df；第5列為F統(tǒng)計(jì)量，它是均方和的比值；第6列為p值?！駊=anova2(X,reps,‘displayopt’)：進(jìn)行平衡(等重復(fù)reps次)雙因子方差分析。當(dāng)‘displayopt’為‘on’（默認(rèn)設(shè)置）時(shí)，激活A(yù)NOVA表的顯示；‘displayopt’為‘off’時(shí)，不予顯示。●[p,table]=anova2(…)：還返回ANOVA表（包含列標(biāo)簽和行標(biāo)簽）?！馵p,table,stats]=anova2(…)：返回stats結(jié)構(gòu)，用于進(jìn)行列因子均值的多重比較檢驗(yàn)。例9－27為了考察高溫合金中碳的含量（因子）和銻與鋁的含量之和（因子）對合金強(qiáng)度的影響。因子取3個(gè)水平0.03、0.04、0.05（上述數(shù)字表示碳的含量占合金總量的百分比），因子取4個(gè)水平3.3、3.4、3.5、3.6（上述數(shù)字的意義同上）。在每個(gè)水平組合下各作一次試驗(yàn)，試驗(yàn)結(jié)果如下：（銻與鋁的含量之和）3.33.43.53.6（碳的含量）0.030.040.0563.163.965.666.865.166.467.869.067.271.071.973.5設(shè)上表中的數(shù)據(jù)符合，請作方差分析。解：作方差分析程序代碼：%平衡雙因子方差分析(不重復(fù))clear;clcy=[63.163.965.666.8;65.166.467.869.0;67.271.071.973.5];p=anova2(y)運(yùn)行結(jié)果：可見，列因子B和行因子A對合金強(qiáng)度的影響均是顯著的。例9－28為了考察某種電池的最大輸出電壓受板極材料與使用電池的環(huán)境溫度的影響，材料類型(因子)取3個(gè)水平(即3種不同的材料)，溫度(因子)也取3個(gè)水平，每個(gè)水平組合下重復(fù)4次試驗(yàn)，數(shù)據(jù)如下：溫度152535材料類型113015517418034408075207082582150188159126136122106115257058453138110168160174120150139961048260解作方差分析程序代碼(材料類型作為列因子)：%兩因子等重復(fù)試驗(yàn)（方差分析）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)運(yùn)行結(jié)果(方差分析表)：可見，列因子A、行因子B以及其交互作用均顯著。運(yùn)行結(jié)果(stats的結(jié)構(gòu))：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)運(yùn)行結(jié)果：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