第4章 概率统计

更新时间:2023-12-13 20:02:01 阅读量: 教育文库 文档下载

说明:文章内容仅供预览,部分内容可能不全。下载后的文档,内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的,是否完整无缺。

MATLAB6.0数学手册 第4章 概率统计

本章介绍MATLAB在概率统计中的若干命令和使用格式,这些命令存放于MatlabR12\\Toolbox\\Stats中。

4.1 随机数的产生

4.1.1 二项分布的随机数据的产生

命令 参数为N,P的二项随机数据 函数 binornd

格式 R = binornd(N,P) %N、P为二项分布的两个参数,返回服从参数为N、P的二

项分布的随机数,N、P大小相同。

R = binornd(N,P,m) %m指定随机数的个数,与R同维数。 R = binornd(N,P,m,n) %m,n分别表示R的行数和列数

例4-1

>> R=binornd(10,0.5) R = 3

>> R=binornd(10,0.5,1,6) R =

8 1 3 7 6 4 >> R=binornd(10,0.5,[1,10]) R =

6 8 4 6 7 5 3 5 6 2 >> R=binornd(10,0.5,[2,3]) R =

7 5 8 6 5 6 >>n = 10:10:60;

>>r1 = binornd(n,1./n) r1 =

2 1 0 1 1 2 >>r2 = binornd(n,1./n,[1 6]) r2 =

0 1 2 1 3 1

4.1.2 正态分布的随机数据的产生

命令 参数为μ、σ的正态分布的随机数据 函数 normrnd

格式 R = normrnd(MU,SIGMA) %返回均值为MU,标准差为SIGMA的正态分布的

134 第4章 概率统计

随机数据,R可以是向量或矩阵。

R = normrnd(MU,SIGMA,m) %m指定随机数的个数,与R同维数。 R = normrnd(MU,SIGMA,m,n) %m,n分别表示R的行数和列数

例4-2

>>n1 = normrnd(1:6,1./(1:6)) n1 =

2.1650 2.3134 3.0250 4.0879 4.8607 6.2827 >>n2 = normrnd(0,1,[1 5]) n2 =

0.0591 1.7971 0.2641 0.8717 -1.4462 >>n3 = normrnd([1 2 3;4 5 6],0.1,2,3) %mu为均值矩阵 n3 =

0.9299 1.9361 2.9640 4.1246 5.0577 5.9864

>> R=normrnd(10,0.5,[2,3]) %mu为10,sigma为0.5的2行3列个正态随机数 R =

9.7837 10.0627 9.4268 9.1672 10.1438 10.5955

4.1.3 常见分布的随机数产生

常见分布的随机数的使用格式与上面相同

表4-1 随机数产生函数表

函数名 Unifrnd Unidrnd Exprnd Normrnd chi2rnd Trnd Frnd gamrnd betarnd lognrnd nbinrnd ncfrnd nctrnd ncx2rnd raylrnd weibrnd binornd geornd hygernd Poissrnd 调用形式 unifrnd ( A,B,m,n) unidrnd(N,m,n) exprnd(Lambda,m,n) normrnd(MU,SIGMA,m,n) chi2rnd(N,m,n) trnd(N,m,n) frnd(N1, N2,m,n) gamrnd(A, B,m,n) betarnd(A, B,m,n) lognrnd(MU, SIGMA,m,n) nbinrnd(R, P,m,n) ncfrnd(N1, N2, delta,m,n) nctrnd(N, delta,m,n) ncx2rnd(N, delta,m,n) raylrnd(B,m,n) weibrnd(A, B,m,n) binornd(N,P,m,n) geornd(P,m,n) hygernd(M,K,N,m,n) poissrnd(Lambda,m,n) 注 释 [A,B]上均匀分布(连续) 随机数 均匀分布(离散)随机数 参数为Lambda的指数分布随机数 参数为MU,SIGMA的正态分布随机数 自由度为N的卡方分布随机数 自由度为N的t分布随机数 第一自由度为N1,第二自由度为N2的F分布随机数 参数为A, B的?分布随机数 参数为A, B的?分布随机数 参数为MU, SIGMA的对数正态分布随机数 参数为R,P的负二项式分布随机数 参数为N1,N2,delta的非中心F分布随机数 参数为N,delta的非中心t分布随机数 参数为N,delta的非中心卡方分布随机数 参数为B的瑞利分布随机数 参数为A, B的韦伯分布随机数 参数为N, p的二项分布随机数 参数为 p的几何分布随机数 参数为 M,K,N的超几何分布随机数 参数为Lambda的泊松分布随机数 4.1.4 通用函数求各分布的随机数据

命令 求指定分布的随机数 函数 random

135 MATLAB6.0数学手册 格式 y = random('name',A1,A2,A3,m,n) %name的取值见表4-2;A1,A2,A3为分

布的参数;m,n指定随机数的行和列

例4-3 产生12(3行4列)个均值为2,标准差为0.3的正态分布随机数

>> y=random('norm',2,0.3,3,4) y =

2.3567 2.0524 1.8235 2.0342 1.9887 1.9440 2.6550 2.3200 2.0982 2.2177 1.9591 2.0178

4.2 随机变量的概率密度计算

4.2.1 通用函数计算概率密度函数值

命令 通用函数计算概率密度函数值 函数 pdf

格式 Y=pdf(name,K,A)

Y=pdf(name,K,A,B) Y=pdf(name,K,A,B,C)

说明 返回在X=K处、参数为A、B、C的概率密度值,对于不同的分布,参数个数是不同;name为分布函数名,其取值如表4-2。

表4-2 常见分布函数表

'beta' 'bino' 'chi2' 'exp' 'f' 'gam' 'geo' 'hyge' 'logn' 'nbin' 'ncf' 'nct' 'ncx2' 'norm' 'poiss' 'rayl' 't' 'unif' 'unid' 'weib' 或 或 或 或 或 或 或 或 或 或 或 或 或 或 或 或 或 或 或 或 name的取值 'Beta' 'Binomial' 'Chisquare' 'Exponential' 'F' 'Gamma' 'Geometric' 'Hypergeometric' 'Lognormal' 'Negative Binomial' 'Noncentral F' 'Noncentral t' 'Noncentral Chi-square' 'Normal' 'Poisson' 'Rayleigh' 'T' 'Uniform' 'Discrete Uniform' 'Weibull' 函数说明 Beta分布 二项分布 卡方分布 指数分布 F分布 GAMMA分布 几何分布 超几何分布 对数正态分布 负二项式分布 非中心F分布 非中心t分布 非中心卡方分布 正态分布 泊松分布 瑞利分布 T分布 均匀分布 离散均匀分布 Weibull分布 例如二项分布:设一次试验,事件A发生的概率为p,那么,在n次独立重复试验中,事件A恰好发生K次的概率P_K为:P_K=P{X=K}=pdf('bino',K,n,p)

例4-4 计算正态分布N(0,1)的随机变量X在点0.6578的密度函数值。

136 第4章 概率统计

解:>> pdf('norm',0.6578,0,1)

ans =

0.3213

例4-5 自由度为8的卡方分布,在点2.18处的密度函数值。 解:>> pdf('chi2',2.18,8)

ans =

0.0363

4.2.2 专用函数计算概率密度函数值

命令 二项分布的概率值 函数 binopdf

格式 binopdf (k, n, p) %等同于pdf(?bino?K,n,p), p — 每次试验事件A发生的概率;K—事件A发生K次;n—试验总次数 命令 泊松分布的概率值 函数 poisspdf

格式 poisspdf(k, Lambda) %等同于pdf(?poiss?,K,Lamda) 命令 正态分布的概率值

函数 normpdf(K,mu,sigma) %计算参数为μ=mu,σ=sigma的正态分布密度函数在K处的值

专用函数计算概率密度函数列表如表4-3。

表4-3 专用函数计算概率密度函数表

函数名 Unifpdf unidpdf Exppdf normpdf chi2pdf Tpdf Fpdf gampdf betapdf lognpdf nbinpdf Ncfpdf Nctpdf ncx2pdf raylpdf weibpdf binopdf geopdf hygepdf poisspdf 调用形式 unifpdf (x, a, b) Unidpdf(x,n) exppdf(x, Lambda) normpdf(x, mu, sigma) chi2pdf(x, n) tpdf(x, n) fpdf(x, n1, n2) gampdf(x, a, b) betapdf(x, a, b) lognpdf(x, mu, sigma) nbinpdf(x, R, P) ncfpdf(x, n1, n2, delta) nctpdf(x, n, delta) ncx2pdf(x, n, delta) raylpdf(x, b) weibpdf(x, a, b) binopdf(x,n,p) geopdf(x,p) hygepdf(x,M,K,N) poisspdf(x,Lambda) 注 释 [a,b]上均匀分布(连续)概率密度在X=x处的函数值 均匀分布(离散)概率密度函数值 参数为Lambda的指数分布概率密度函数值 参数为mu,sigma的正态分布概率密度函数值 自由度为n的卡方分布概率密度函数值 自由度为n的t分布概率密度函数值 第一自由度为n1,第二自由度为n2的F分布概率密度函数值 参数为a, b的?分布概率密度函数值 参数为a, b的?分布概率密度函数值 参数为mu, sigma的对数正态分布概率密度函数值 参数为R,P的负二项式分布概率密度函数值 参数为n1,n2,delta的非中心F分布概率密度函数值 参数为n,delta的非中心t分布概率密度函数值 参数为n,delta的非中心卡方分布概率密度函数值 参数为b的瑞利分布概率密度函数值 参数为a, b的韦伯分布概率密度函数值 参数为n, p的二项分布的概率密度函数值 参数为 p的几何分布的概率密度函数值 参数为 M,K,N的超几何分布的概率密度函数值 参数为Lambda的泊松分布的概率密度函数值 例4-6 绘制卡方分布密度函数在自由度分别为1、5、15的图形

>> x=0:0.1:30;

>> y1=chi2pdf(x,1); plot(x,y1,':') >> hold on

137 MATLAB6.0数学手册 >> y2=chi2pdf(x,5);plot(x,y2,'+') >> y3=chi2pdf(x,15);plot(x,y3,'o')

>> axis([0,30,0,0.2]) %指定显示的图形区域

则图形为图4-1。

4.2.3 常见分布的密度函数作图

1.二项分布 例4-7

>>x = 0:10;

>>y = binopdf(x,10,0.5); >>plot(x,y,'+')

图4-1

2.卡方分布 例4-8

>> x = 0:0.2:15; >>y = chi2pdf(x,4); >>plot(x,y)

0.250.20.150.20.150.10.10.05002468100.05 图4-2

0051015

3.非中心卡方分布 例4-9

>>x = (0:0.1:10)';

>>p1 = ncx2pdf(x,4,2); >>p = chi2pdf(x,4); >>plot(x,p,'--',x,p1,'-')

4.指数分布 例4-10

>>x = 0:0.1:10; >>y = exppdf(x,2); >>plot(x,y)

0.20.50.150.40.10.30.20.050.100246810 图4-3

00246810

138

第4章 概率统计

5.F分布 例4-11

>>x = 0:0.01:10; >>y = fpdf(x,5,3); >>plot(x,y)

6.非中心F分布 例4-12

>>x = (0.01:0.1:10.01)'; >>p1 = ncfpdf(x,5,20,10); >>p = fpdf(x,5,20); >>plot(x,p,'--',x,p1,'-')

0.80.80.60.60.40.40.20.2002468100

024681012

图4-4

7.Γ分布 例4-13

>>x = gaminv((0.005:0.01:0.995),100,10); >>y = gampdf(x,100,10); >>y1 = normpdf(x,1000,100); >>plot(x,y,'-',x,y1,'-.')

8.对数正态分布 例4-14

>>x = (10:1000:125010)';

>>y = lognpdf(x,log(20000),1.0); >>plot(x,y)

>>set(gca,'xtick',[0 30000 60000 90000 120000])

>>set(gca,'xticklabel',str2mat('0','$30,000','$60,000',? '$90,000','$120,000'))

x 103.5-53x 105-32.54231.52110.507008009001000110012001300

00 $30,000 $60,000 $90,000 $120,000

图4-5

139 MATLAB6.0数学手册 9.负二项分布 例4-15

>>x = (0:10);

>>y = nbinpdf(x,3,0.5); >>plot(x,y,'+')

10.正态分布 例4-16

>> x=-3:0.2:3;

>> y=normpdf(x,0,1); >> plot(x,y)

0.20.40.150.30.10.20.050.100246810 图4-6

0-3-2-10123

11.泊松分布 例4-17

>>x = 0:15;

>>y = poisspdf(x,5); >>plot(x,y,'+')

12.瑞利分布 例4-18

>>x = [0:0.01:2]; >>p = raylpdf(x,0.5); >>plot(x,p)

0.21.50.1510.10.50.050051015 图4-7

000.511.52

13.T分布 例4-19

>>x = -5:0.1:5; >>y = tpdf(x,5);

>>z = normpdf(x,0,1); >>plot(x,y,'-',x,z,'-.')

14.威布尔分布

140 第4章 概率统计

例4-20

>> t=0:0.1:3;

>> y=weibpdf(t,2,2); >> plot(y)

0.4

1.50.310.20.50.10-505 图4-8

005101520253035

4.3 随机变量的累积概率值(分布函数值)

4.3.1 通用函数计算累积概率值

命令 通用函数cdf用来计算随机变量X?K的概率之和(累积概率值) 函数 cdf

格式 cdf(?name?,K,A)

cdf(?name?,K,A,B)

cdf(?name?,K,A,B,C)

说明 返回以name为分布、随机变量X≤K的概率之和的累积概率值,name的取值见表4-1 常见分布函数表

例4-21 求标准正态分布随机变量X落在区间(-∞,0.4)内的概率(该值就是概率统计教材中的附表:标准正态数值表)。

解:

>> cdf('norm',0.4,0,1) ans =

0.6554

例4-22 求自由度为16的卡方分布随机变量落在[0,6.91]内的概率

>> cdf('chi2',6.91,16) ans =

0.0250

4.3.2 专用函数计算累积概率值(随机变量X?K的概率之和)

命令 二项分布的累积概率值 函数 binocdf

格式 binocdf (k, n, p) %n为试验总次数,p为每次试验事件A发生的概率,k为n

次试验中事件A发生的次数,该命令返回n次试验中事件A恰好发生k次的概率。

141 MATLAB6.0数学手册 命令 正态分布的累积概率值 函数 normcdf

格式 normcdf(x,mu,sigma) %返回F(x)=???p(t)dt的值,mu、sigma为正态分布的

两个参数

例4-23 设X~N(3, 22)

(1)求P{2?X?5},P{?4?X?10},P{X?2},P{X?3} (2)确定c,使得P{X?c}?P{X?c} 解(1) p1=P{2?X?5} p2=P{?4?X?10}

p3=P{X?2}?1?P{X?2} p4=P{X?3}?1?P{X?3}

则有:

>>p1=normcdf(5,3,2)-normcdf(2,3,2) p1 =

0.5328

>>p2=normcdf(10,3,2)-normcdf(-4,3,2) p2 =

0.9995

>>p3=1-normcdf(2,3,2)-normcdf(-2,3,2) p3 =

0.6853

>>p4=1-normcdf(3,3,2) p4 =

0.5000

x专用函数计算累积概率值函数列表如表4-4。

表4-4 专用函数的累积概率值函数表

函数名 unifcdf unidcdf expcdf normcdf chi2cdf tcdf fcdf gamcdf betacdf logncdf nbincdf ncfcdf nctcdf ncx2cdf raylcdf weibcdf binocdf geocdf hygecdf poisscdf 142 调用形式 unifcdf (x, a, b) unidcdf(x,n) expcdf(x, Lambda) normcdf(x, mu, sigma) chi2cdf(x, n) tcdf(x, n) fcdf(x, n1, n2) gamcdf(x, a, b) betacdf(x, a, b) logncdf(x, mu, sigma) nbincdf(x, R, P) ncfcdf(x, n1, n2, delta) nctcdf(x, n, delta) ncx2cdf(x, n, delta) raylcdf(x, b) weibcdf(x, a, b) binocdf(x,n,p) geocdf(x,p) hygecdf(x,M,K,N) poisscdf(x,Lambda) 注 释 [a,b]上均匀分布(连续)累积分布函数值 F(x)=P{X≤x} 均匀分布(离散)累积分布函数值 F(x)=P{X≤x} 参数为Lambda的指数分布累积分布函数值 F(x)=P{X≤x} 参数为mu,sigma的正态分布累积分布函数值 F(x)=P{X≤x} 自由度为n的卡方分布累积分布函数值 F(x)=P{X≤x} 自由度为n的t分布累积分布函数值 F(x)=P{X≤x} 第一自由度为n1,第二自由度为n2的F分布累积分布函数值 参数为a, b的?分布累积分布函数值 F(x)=P{X≤x} 参数为a, b的?分布累积分布函数值 F(x)=P{X≤x} 参数为mu, sigma的对数正态分布累积分布函数值 参数为R,P的负二项式分布概累积分布函数值 F(x)=P{X≤x} 参数为n1,n2,delta的非中心F分布累积分布函数值 参数为n,delta的非中心t分布累积分布函数值 F(x)=P{X≤x} 参数为n,delta的非中心卡方分布累积分布函数值 参数为b的瑞利分布累积分布函数值 F(x)=P{X≤x} 参数为a, b的韦伯分布累积分布函数值 F(x)=P{X≤x} 参数为n, p的二项分布的累积分布函数值 F(x)=P{X≤x} 参数为 p的几何分布的累积分布函数值 F(x)=P{X≤x} 参数为 M,K,N的超几何分布的累积分布函数值 参数为Lambda的泊松分布的累积分布函数值 F(x)=P{X≤x} 第4章 概率统计

说明 累积概率函数就是分布函数F(x)=P{X≤x}在x处的值。

4.4 随机变量的逆累积分布函数

MATLAB中的逆累积分布函数是已知F(x)?P{X?x},求x。 逆累积分布函数值的计算有两种方法

4.4.1 通用函数计算逆累积分布函数值

命令 icdf 计算逆累积分布函数

f?name?,P,a1,a2,a3) 格式 icd(说明 返回分布为name,参数为a1,a2,a3,累积概率值为P的临界值,这里name与前

面表4.1相同。

如果P?cdf(?name?,x,a1,a2,a3),则x?icdf(?name?,P,a1,a2,a3) 例4-24 在标准正态分布表中,若已知?(x)=0.975,求x 解:>> x=icdf('norm',0.975,0,1)

x =

1.9600

例4-25 在?2分布表中,若自由度为10,?=0.975,求临界值Lambda。

解:因为表中给出的值满足P{?2??}??,而逆累积分布函数icdf求满足P{?2??}??的临界值?。所以,这里的?取为0.025,即

>> Lambda=icdf('chi2',0.025,10) Lambda = 3.2470

例4-26 在假设检验中,求临界值问题:

已知:??0.05,查自由度为10的双边界检验t分布临界值

>>lambda=icdf('t',0.025,10) lambda =

-2.2281

4.4.2 专用函数-inv计算逆累积分布函数

命令 正态分布逆累积分布函数 函数 norminv

格式 X=norminv(p,mu,sigma) %p为累积概率值,mu为均值,sigma为标准差,X

为临界值,满足:p=P{X≤x}。 例4-27 设X~N(3,22),确定c使得P{X?c}?P{X?c}。 解:由P{X?c}?P{X?c}得,P{X?c}?P{X?c}=0.5,所以

>>X=norminv(0.5, 3, 2) X= 3

关于常用临界值函数可查下表4-5。

143

MATLAB6.0数学手册 表4-5 常用临界值函数表

函数名 unifinv unidinv expinv norminv chi2inv tinv finv gaminv betainv logninv nbininv ncfinv nctinv ncx2inv raylinv weibinv binoinv geoinv hygeinv poissinv 调用形式 x=unifinv (p, a, b) x=unidinv (p,n) x=expinv (p, Lambda) x=Norminv(x,mu,sigma) x=chi2inv (x, n) x=tinv (x, n) x=finv (x, n1, n2) x=gaminv (x, a, b) x=betainv (x, a, b) x=logninv (x, mu, sigma) x=nbininv (x, R, P) x=ncfinv (x, n1, n2, delta) x=nctinv (x, n, delta) x=ncx2inv (x, n, delta) x=raylinv (x, b) x=weibinv (x, a, b) x=binoinv (x,n,p) x=geoinv (x,p) x=hygeinv (x,M,K,N) x=poissinv (x,Lambda) 注 释 均匀分布(连续)逆累积分布函数(P=P{X≤x},求x) 均匀分布(离散)逆累积分布函数,x为临界值 指数分布逆累积分布函数 正态分布逆累积分布函数 卡方分布逆累积分布函数 t分布累积分布函数 F分布逆累积分布函数 ?分布逆累积分布函数 ?分布逆累积分布函数 对数正态分布逆累积分布函数 负二项式分布逆累积分布函数 非中心F分布逆累积分布函数 非中心t分布逆累积分布函数 非中心卡方分布逆累积分布函数 瑞利分布逆累积分布函数 韦伯分布逆累积分布函数 二项分布的逆累积分布函数 几何分布的逆累积分布函数 超几何分布的逆累积分布函数 泊松分布的逆累积分布函数 例4-28 公共汽车门的高度是按成年男子与车门顶碰头的机会不超过1%设计的。设男子身高X(单位:cm)服从正态分布N(175,36),求车门的最低高度。

解:设h为车门高度,X为身高

求满足条件P{X?h}?0.01的h,即P{X?h}?0.99,所以

>>h=norminv(0.99, 175, 6) h =

188.9581

例4-29 卡方分布的逆累积分布函数的应用 在MATLAB的编辑器下建立M文件如下:

n=5; a=0.9; %n为自由度,a为置信水平或累积概率 x_a=chi2inv(a,n); %x_a 为临界值

x=0:0.1:15; yd_c=chi2pdf(x,n); %计算?2(5)的概率密度函数值,供绘图用 plot(x,yd_c,'b'), hold on %绘密度函数图形

xxf=0:0.1:x_a; yyf=chi2pdf(xxf,n); %计算[0,x_a]上的密度函数值,供填色用 fill([xxf,x_a], [yyf,0], 'g') %填色,其中:点(x_a, 0)使得填色区域封闭 text(x_a*1.01,0.01, num2str(x_a)) %标注临界值点 text(10,0.10, ['\\fontsize{16}X~{\\chi}^2(4)'])

%图中标注

text(1.5,0.05, '\\fontsize{22}alpha=0.9' ) %图中标注 结果显示如图4-9。

图4-9

144 第4章 概率统计

4.5 随机变量的数字特征

4.5.1 平均值、中值

命令 利用mean求算术平均值

格式 mean(X) %X为向量,返回X中各元素的平均值

mean(A) %A为矩阵,返回A中各列元素的平均值构成的向量 mean(A,dim) %在给出的维数内的平均值 说明 X为向量时,算术平均值的数学含义是x?1?xi,即样本均值。

ni?1例4-30

>> A=[1 3 4 5;2 3 4 6;1 3 1 5] A =

1 3 4 5 2 3 4 6 1 3 1 5 >> mean(A) ans =

1.3333 3.0000 3.0000 5.3333 >> mean(A,1) ans =

1.3333 3.0000 3.0000 5.3333

n命令 忽略NaN计算算术平均值

格式 nanmean(X) %X为向量,返回X中除NaN外元素的算术平均值。

nanmean(A) %A为矩阵,返回A中各列除NaN外元素的算术平均值向量。 例4-31

>> A=[1 2 3;nan 5 2;3 7 nan] A =

1 2 3 NaN 5 2 3 7 NaN >> nanmean(A) ans =

2.0000 4.6667 2.5000

命令 利用median计算中值(中位数)

格式 median(X) %X为向量,返回X中各元素的中位数。

median(A) %A为矩阵,返回A中各列元素的中位数构成的向量。 median(A,dim) %求给出的维数内的中位数 例4-32

>> A=[1 3 4 5;2 3 4 6;1 3 1 5] A =

1 3 4 5 2 3 4 6 1 3 1 5 >> median(A) ans =

145 MATLAB6.0数学手册 1 3 4 5

命令 忽略NaN计算中位数

格式 nanmedian(X) %X为向量,返回X中除NaN外元素的中位数。

nanmedian(A) %A为矩阵,返回A中各列除NaN外元素的中位数向量。 例4-33

>> A=[1 2 3;nan 5 2;3 7 nan] A =

1 2 3 NaN 5 2 3 7 NaN >> nanmedian(A) ans =

2.0000 5.0000 2.5000

命令 利用geomean计算几何平均数

格式 M=geomean(X) %X为向量,返回X中各元素的几何平均数。

M=geomean(A) %A为矩阵,返回A中各列元素的几何平均数构成的向量。 说明 几何平均数的数学含义是M?(?xi)n,其中:样本数据非负,主要用于对数正

i?1n1态分布。

例4-34

>> B=[1 3 4 5] B =

1 3 4 5 >> M=geomean(B) M =

2.7832

>> A=[1 3 4 5;2 3 4 6;1 3 1 5] A =

1 3 4 5 2 3 4 6 1 3 1 5 >> M=geomean(A) M =

1.2599 3.0000 2.5198 5.3133

命令 利用harmmean求调和平均值

格式 M=harmmean(X) %X为向量,返回X中各元素的调和平均值。

M=harmmean(A) %A为矩阵,返回A中各列元素的调和平均值构成的向量。 说明 调和平均值的数学含义是M?1?xii?1nn,其中:样本数据非0,主要用于严重偏斜

分布。

例4-35

>> B=[1 3 4 5] B =

1 3 4 5 >> M=harmmean(B) M =

146 第4章 概率统计

2.2430

>> A=[1 3 4 5;2 3 4 6;1 3 1 5] A =

1 3 4 5 2 3 4 6 1 3 1 5 >> M=harmmean(A) M =

1.2000 3.0000 2.0000 5.2941

4.5.2 数据比较

命令 排序

格式 Y=sort(X) %X为向量,返回X按由小到大排序后的向量。

Y=sort(A) %A为矩阵,返回A的各列按由小到大排序后的矩阵。

[Y,I]=sort(A) % Y为排序的结果,I中元素表示Y中对应元素在A中位置。

sort(A,dim) %在给定的维数dim内排序 说明 若X为复数,则通过|X|排序。 例4-36

>> A=[1 2 3;4 5 2;3 7 0] A =

1 2 3 4 5 2 3 7 0 >> sort(A) ans =

1 2 0 3 5 2 4 7 3 >> [Y,I]=sort(A) Y =

1 2 0 3 5 2 4 7 3 I =

1 1 3 3 2 2 2 3 1

命令 按行方式排序 函数 sortrows

格式 Y=sortrows(A) %A为矩阵,返回矩阵Y,Y按A的第1列由小到大,以

行方式排序后生成的矩阵。

Y=sortrows(A, col) %按指定列col由小到大进行排序 [Y,I]=sortrows(A, col) % Y为排序的结果,I表示Y中第col列元素在A中位置。 说明 若X为复数,则通过|X|的大小排序。 例4-37

>> A=[1 2 3;4 5 2;3 7 0] A =

1 2 3 4 5 2

147 MATLAB6.0数学手册 3 7 0 >> sortrows(A) ans =

1 2 3 3 7 0 4 5 2 >> sortrows(A,1) ans =

1 2 3 3 7 0 4 5 2 >> sortrows(A,3) ans =

3 7 0 4 5 2 1 2 3 >> sortrows(A,[3 2]) ans =

3 7 0 4 5 2 1 2 3 >> [Y,I]=sortrows(A,3) Y =

3 7 0 4 5 2 1 2 3 I = 3 2 1

命令 求最大值与最小值之差 函数 range

格式 Y=range(X) %X为向量,返回X中的最大值与最小值之差。

Y=range(A) %A为矩阵,返回A中各列元素的最大值与最小值之差。 例4-38

>> A=[1 2 3;4 5 2;3 7 0] A =

1 2 3 4 5 2 3 7 0 >> Y=range(A) Y =

3 5 3

4.5.3 期望

命令 计算样本均值 函数 mean

格式 用法与前面一样

例4-39 随机抽取6个滚珠测得直径如下:(直径:mm)

14.70 15.21 14.90 14.91 15.32 15.32

试求样本平均值

解:>>X=[14.70 15.21 14.90 14.91 15.32 15.32];

>>mean(X) %计算样本均值

148

第4章 概率统计

则结果如下:

ans =

15.0600

命令 由分布律计算均值 利用sum函数计算

例4-40 设随机变量X的分布律为:

X P -2 0.3 -1 0.1 0 0.2 1 0.1 2 0.3 求E (X) E(X2-1)

解:在Matlab编辑器中建立M文件如下:

X=[-2 -1 0 1 2];

p=[0.3 0.1 0.2 0.1 0.3]; EX=sum(X.*p) Y=X.^2-1

EY=sum(Y.*p)

运行后结果如下:

EX = 0 Y =

3 0 -1 0 3 EY =

1.6000

4.5.4 方差

命令 求样本方差 函数 var

n1格式 D=var(X) %var(X)=s?则返回向量的样本方差。 (xi?X)2,若X为向量,?n?1i?12 D=var(A) %A为矩阵,则D为A的列向量的样本方差构成的行向量。

D=var(X, 1) %返回向量(矩阵)X的简单方差(即置前因子为1的方差)

nD=var(X, w) %返回向量(矩阵)X的以w为权重的方差 命令 求标准差 函数 std

格式 std(X) %返回向量(矩阵)X的样本标准差(置前因子为1)即:

n?1std?1nx?X in?1?i?1std(X,1) %返回向量(矩阵)X的标准差(置前因子为

1) nstd(X, 0) %与std (X)相同

std(X, flag, dim) %返回向量(矩阵)中维数为dim的标准差值,其中flag=0

时,置前因子为1;否则置前因子为1。

nn?1

149 MATLAB6.0数学手册 例4-41 求下列样本的样本方差和样本标准差,方差和标准差

14.70 15.21 14.90 15.32 15.32

解:

>>X=[14.7 15.21 14.9 14.91 15.32 15.32]; >>DX=var(X,1) %方差 DX =

0.0559

>>sigma=std(X,1) %标准差 sigma =

0.2364

>>DX1=var(X) %样本方差 DX1 =

0.0671

>>sigma1=std(X) %样本标准差 sigma1 = 0.2590

命令 忽略NaN的标准差 函数 nanstd

格式 y = nanstd(X) %若X为含有元素NaN的向量,则返回除NaN外的元素的标准

差,若X为含元素NaN的矩阵,则返回各列除NaN外的标准差构成的向量。

例4-42

>> M=magic(3) %产生3阶魔方阵 M =

8 1 6 3 5 7 4 9 2

>> M([1 6 8])=[NaN NaN NaN] %替换3阶魔方阵中第1、6、8个元素为NaN M =

NaN 1 6 3 5 NaN 4 NaN 2

>> y=nanstd(M) %求忽略NaN的各列向量的标准差 y =

0.7071 2.8284 2.8284

>> X=[1 5]; %忽略NaN的第2列元素

>> y2=std(X) %验证第2列忽略NaN元素的标准差 y2 =

2.8284

命令 样本的偏斜度 函数 skewness

格式 y = skewness(X) %X为向量,返回X的元素的偏斜度;X为矩阵,返回X各

列元素的偏斜度构成的行向量。

y = skewness(X,flag) %flag=0表示偏斜纠正,flag=1(默认)表示偏斜不纠正。 说明 偏斜度样本数据关于均值不对称的一个测度,如果偏斜度为负,说明均值左边的数据比均值右边的数据更散;如果偏斜度为正,说明均值右边的数据比均值左边的数据更散,

E(x??)3因而正态分布的偏斜度为 0;偏斜度是这样定义的:y?

?3 150 第4章 概率统计

其中:μ为x的均值,σ为x的标准差,E(.)为期望值算子 例4-43

>> X=randn([5,4]) X =

0.2944 0.8580 -0.3999 0.6686 -1.3362 1.2540 0.6900 1.1908 0.7143 -1.5937 0.8156 -1.2025 1.6236 -1.4410 0.7119 -0.0198 -0.6918 0.5711 1.2902 -0.1567 >> y=skewness(X) y =

-0.0040 -0.3136 -0.8865 -0.2652 >> y=skewness(X,0) y =

-0.0059 -0.4674 -1.3216 -0.3954

4.5.5 常见分布的期望和方差

命令 均匀分布(连续)的期望和方差 函数 unifstat

格式 [M,V] = unifstat(A,B) %A、B为标量时,就是区间上均匀分布的期望和方差,

A、B也可为向量或矩阵,则M、V也是向量或矩阵。

例4-44

>>a = 1:6; b = 2.*a; >>[M,V] = unifstat(a,b) M =

1.5000 3.0000 4.5000 6.0000 7.5000 9.0000 V =

0.0833 0.3333 0.7500 1.3333 2.0833 3.0000

命令 正态分布的期望和方差 函数 normstat

格式 [M,V] = normstat(MU,SIGMA) %MU、SIGMA可为标量也可为向量或矩阵,

则M=MU,V=SIGMA2。

例4-45

>> n=1:4;

>> [M,V]=normstat(n'*n,n'*n) M =

1 2 3 4 2 4 6 8 3 6 9 12 4 8 12 16 V =

1 4 9 16 4 16 36 64 9 36 81 144 16 64 144 256

命令 二项分布的均值和方差 函数 binostat

格式 [M,V] = binostat(N,P) %N,P为二项分布的两个参数,可为标量也可为向量

或矩阵。

151 MATLAB6.0数学手册 例4-46

>>n = logspace(1,5,5) n =

10 100 1000 10000 100000 >>[M,V] = binostat(n,1./n) M =

1 1 1 1 1 V =

0.9000 0.9900 0.9990 0.9999 1.0000 >>[m,v] = binostat(n,1/2) m =

5 50 500 5000 50000 v =

1.0e+04 *

0.0003 0.0025 0.0250 0.2500 2.5000

常见分布的期望和方差见下表4-6。

表4-6 常见分布的均值和方差

函数名 unifstat unidstat expstat normstat chi2stat tstat fstat gamstat betastat lognstat nbinstat ncfstat nctstat ncx2stat raylstat Weibstat Binostat Geostat hygestat Poisstat 调用形式 [M,V]=unifstat ( a, b) [M,V]=unidstat (n) [M,V]=expstat (p, Lambda) [M,V]=normstat(mu,sigma) [M,V]=chi2stat (x, n) [M,V]=tstat ( n) [M,V]=fstat ( n1, n2) [M,V]=gamstat ( a, b) [M,V]=betastat ( a, b) [M,V]=lognstat ( mu, sigma) [M,V]=nbinstat ( R, P) [M,V]=ncfstat ( n1, n2, delta) [M,V]=nctstat ( n, delta) [M,V]=ncx2stat ( n, delta) [M,V]=raylstat ( b) [M,V]=weibstat ( a, b) [M,V]=binostat (n,p) [M,V]=geostat (p) [M,V]=hygestat (M,K,N) [M,V]=poisstat (Lambda) 注 释 均匀分布(连续)的期望和方差,M为期望,V为方差 均匀分布(离散)的期望和方差 指数分布的期望和方差 正态分布的期望和方差 卡方分布的期望和方差 t分布的期望和方差 F分布的期望和方差 ?分布的期望和方差 ?分布的期望和方差 对数正态分布的期望和方差 负二项式分布的期望和方差 非中心F分布的期望和方差 非中心t分布的期望和方差 非中心卡方分布的期望和方差 瑞利分布的期望和方差 韦伯分布的期望和方差 二项分布的期望和方差 几何分布的期望和方差 超几何分布的期望和方差 泊松分布的期望和方差 4.5.6 协方差与相关系数

命令 协方差 函数 cov

格式 cov(X) %求向量X的协方差

cov(A) %求矩阵A的协方差矩阵,该协方差矩阵的对角线元素是A的各列的

方差,即:var(A)=diag(cov(A))。

cov(X,Y) %X,Y为等长列向量,等同于cov([X Y])。

152 第4章 概率统计

例4-47

>> X=[0 -1 1]';Y=[1 2 2]';

>> C1=cov(X) %X的协方差 C1 = 1

>> C2=cov(X,Y) %列向量X、Y的协方差矩阵,对角线元素为各列向量的方差 C2 =

1.0000 0 0 0.3333 >> A=[1 2 3;4 0 -1;1 7 3] A =

1 2 3 4 0 -1 1 7 3

>> C1=cov(A) %求矩阵A的协方差矩阵 C1 =

3.0000 -4.5000 -4.0000 -4.5000 13.0000 6.0000 -4.0000 6.0000 5.3333

>> C2=var(A(:,1)) %求A的第1列向量的方差 C2 = 3

>> C3=var(A(:,2)) %求A的第2列向量的方差 C3 = 13

>> C4=var(A(:,3)) C4 =

5.3333

命令 相关系数 函数 corrcoef

格式 corrcoef(X,Y) %返回列向量X,Y的相关系数,等同于corrcoef([X corrcoef (A) %返回矩阵A的列向量的相关系数矩阵 例4-48

>> A=[1 2 3;4 0 -1;1 3 9] A =

1 2 3 4 0 -1 1 3 9

>> C1=corrcoef(A) %求矩阵A的相关系数矩阵 C1 =

1.0000 -0.9449 -0.8030 -0.9449 1.0000 0.9538 -0.8030 0.9538 1.0000

>> C1=corrcoef(A(:,2),A(:,3)) %求A的第2列与第3列列向量的相关系数矩阵 C1 =

1.0000 0.9538 0.9538 1.0000

4.6 统计作图

4.6.1 正整数的频率表

命令 正整数的频率表

Y])。

153

MATLAB6.0数学手册 函数 tabulate

格式 table = tabulate(X) %X为正整数构成的向量,返回3列:第1列中包含X的值

第2列为这些值的个数,第3列为这些值的频率。

例4-49

>> A=[1 2 2 5 6 3 8] A =

1 2 2 5 6 3 8 >> tabulate(A)

Value Count Percent 1 1 14.29% 2 2 28.57% 3 1 14.29% 4 0 0.00% 5 1 14.29% 6 1 14.29% 7 0 0.00% 8 1 14.29%

4.6.2 经验累积分布函数图形

函数 cdfplot

格式 cdfplot(X) %作样本X(向量)的累积分布函数图形

h = cdfplot(X) %h表示曲线的环柄

[h,stats] = cdfplot(X) %stats表示样本的一些特征 例4-50

Empirical CDF1>> X=normrnd (0,1,50,1); >> [h,stats]=cdfplot(X) h =

3.0013 stats =

min: -1.8740 %样本最小值 max: 1.6924 %最大值 mean: 0.0565 %平均值 median: 0.1032 %中间值

std: 0.7559 %样本标准差

0.90.80.70.6F(x)0.50.40.30.20.10-2-1.5-1-0.50x0.511.524.6.3 最小二乘拟合直线

函数 lsline

格式 lsline %最小二乘拟合直线 h = lsline %h为直线的句柄 例4-51

>> X = [2 3.4 5.6 8 11 12.3 13.8 16 18.8 19.9]'; >> plot(X,'+') >> lsline

201816141210864图4-10

4.6.4 绘制正态分布概率图形

20246810函数 normplot

格式 normplot(X) %若X为向量,则显示正态分布概率图形,若X为矩阵,则显

示每一列的正态分布概率图形。

154 图4-11

第4章 概率统计

h = normplot(X) %返回绘图直线的句柄

说明 样本数据在图中用“+”显示;如果数据来自正态分布,则图形显示为直线,而其它分布可能在图中产生弯曲。

例4-53

>> X=normrnd(0,1,50,1); >> normplot(X)

图4-12

4.6.5 绘制威布尔(Weibull)概率图形

函数 weibplot

格式 weibplot(X) %若X为向量,则显示威布尔(Weibull)概率图形,若X为矩阵,

则显示每一列的威布尔概率图形。

h = weibplot(X) %返回绘图直线的柄

说明 绘制威布尔(Weibull)概率图形的目的是用图解法估计来自威布尔分布的数据X,如果X是威布尔分布数据,其图形是直线的,否则图形中可能产生弯曲。

例4-54

>> r = weibrnd(1.2,1.5,50,1); >> weibplot(r)

Weibull Probability Plot0.99 0.96 0.90 0.75 0.50 Probability0.25 0.10 0.05 0.02 0.01 -101010Data

图4-13

4.6.6 样本数据的盒图

函数 boxplot

格式 boxplot(X) %产生矩阵X的每一列的盒图和“须”图,“须”是从盒的尾部延

155 MATLAB6.0数学手册 伸出来,并表示盒外数据长度的线,如果“须”的外面没有数据,则在“须”的底部有一个点。

boxplot(X,notch) %当notch=1时,产生一凹盒图,notch=0时产生一矩箱图。 boxplot(X,notch,'sym') %sym表示图形符号,默认值为“+”。

boxplot(X,notch,'sym',vert) %当vert=0时,生成水平盒图,vert=1时,生成竖

直盒图(默认值vert=1)。

boxplot(X,notch,'sym',vert,whis) %whis定义“须”图的长度,默认值为1.5,若

whis=0则boxplot函数通过绘制sym符号图来显示盒外的所有数据值。

例4-55

>>x1 = normrnd(5,1,100,1); >>x2 = normrnd(6,1,100,1); >>x = [x1 x2];

>> boxplot(x,1,'g+',1,0)

图4-14

4.6.7 给当前图形加一条参考线

函数 refline

格式 refline(slope,intercept) % slope表示直线斜率,intercept表示截距

refline(slope) slope=[a b],图中加一条直线:y=b+ax。 例4-56

>>y = [3.2 2.6 3.1 3.4 2.4 2.9 3.0 3.3 3.2 2.1 2.6]'; >>plot(y,'+') >>refline(0,3)

图4-15

4.6.8 在当前图形中加入一条多项式曲线

函数 refcurve

156 第4章 概率统计

格式 h = refcurve(p) %在图中加入一条多项式曲线,h为曲线的环柄,p为多项式系

数向量,p=[p1,p2, p3,…,pn],其中p1为最高幂项系数。

例4-57 火箭的高度与时间图形,加入一条理论高度曲线,火箭初速为100m/秒。

>>h = [85 162 230 289 339 381 413 437 452 458 456 440 400 356]; >>plot(h,'+')

>>refcurve([-4.9 100 0])

图4-16

4.6.9 样本的概率图形

函数 capaplot

格式 p = capaplot(data,specs) úta为所给样本数据,specs指定范围,p表示在指定

范围内的概率。

说明 该函数返回来自于估计分布的随机变量落在指定范围内的概率 例4-58

>> data=normrnd (0,1,30,1); >> p=capaplot(data,[-2,2]) p =

0.9199

图4-17

4.6.10 附加有正态密度曲线的直方图

函数 histfit

格式 histfit(data) úta为向量,返回直方图

和正态曲线。

histfit(data,nbins) % nbins指定bar的个数,

缺省时为data中数据个数的平方根。

例4-59

>>r = normrnd (10,1,100,1); >>histfit(r)

图4-18

157 MATLAB6.0数学手册 4.6.11 在指定的界线之间画正态密度曲线

函数 normspec

格式 p = normspec(specs,mu,sigma) %specs指定界线,mu,sigma为正态分布的参数p

为样本落在上、下界之间的概率。

例4-60

>>normspec([10 Inf],11.5,1.25)

图4-19

4.7 参数估计

4.7.1 常见分布的参数估计

命令 β分布的参数a和b的最大似然估计值和置信区间 函数 betafit

格式 PHAT=betafit(X)

[PHAT,PCI]=betafit(X,ALPHA)

说明 PHAT为样本X的β分布的参数a和b的估计量

PCI为样本X的β分布参数a和b的置信区间,是一个2×2矩阵,其第1例为参数a的置信下界和上界,第2例为b的置信下界和上界,ALPHA为显著水平,(1-α)×100%为置信度。

例4-61 随机产生100个β分布数据,相应的分布参数真值为4和3。则4和3的最大似然估计值和置信度为99%的置信区间为:

解:

>>X = betarnd (4,3,100,1); %产生100个β分布的随机数

>>[PHAT,PCI] = betafit(X,0.01) %求置信度为99%的置信区间和参数a、b的估计值

结果显示

PHAT =

3.9010 2.6193 PCI =

2.5244 1.7488 5.2776 3.4898

158

第4章 概率统计

说明 估计值3.9010的置信区间是[2.5244 5.2776],估计值2.6193的置信区间是[1.7488 3.4898]。

命令 正态分布的参数估计 函数 normfit

格式 [muhat,sigmahat,muci,sigmaci] = normfit(X)

[muhat,sigmahat,muci,sigmaci] = normfit(X,alpha)

说明 muhat,sigmahat分别为正态分布的参数μ和σ的估计值,muci,sigmaci分别为置信区间,其置信度为(1??)?100%;alpha给出显著水平α,缺省时默认为0.05,即置信度为95%。

例4-62 有两组(每组100个元素)正态随机数据,其均值为10,均方差为2,求95%的置信区间和参数估计值。

解:>>r = normrnd (10,2,100,2); %产生两列正态随机数据

>>[mu,sigma,muci,sigmaci] = normfit(r)

则结果为

mu =

10.1455 10.0527 %各列的均值的估计值 sigma =

1.9072 2.1256 %各列的均方差的估计值 muci =

9.7652 9.6288 10.5258 10.4766 sigmaci =

1.6745 1.8663 2.2155 2.4693

说明 muci,sigmaci中各列分别为原随机数据各列估计值的置信区间,置信度为95%。 例4-63 分别使用金球和铂球测定引力常数

(1)用金球测定观察值为:6.683 6.681 6.676 6.678 6.679 6.672 (2)用铂球测定观察值为:6.661 6.661 6.667 6.667 6.664

设测定值总体为N(?,?2),μ和σ为未知。对(1)、(2)两种情况分别求μ和σ的置信度为0.9的置信区间。

解:建立M文件:LX0833.m

X=[6.683 6.681 6.676 6.678 6.679 6.672]; Y=[6.661 6.661 6.667 6.667 6.664];

[mu,sigma,muci,sigmaci]=normfit(X,0.1) %金球测定的估计 [MU,SIGMA,MUCI,SIGMACI]=normfit(Y,0.1) %铂球测定的估计

运行后结果显示如下:

mu =

6.6782 sigma =

0.0039 muci =

6.6750 6.6813 sigmaci =

0.0026 0.0081 MU =

159 MATLAB6.0数学手册 6.6640 SIGMA =

0.0030 MUCI =

6.6611 6.6669 SIGMACI =

0.0019 0.0071

由上可知,金球测定的μ估计值为6.6782,置信区间为[6.6750,6.6813]; σ的估计值为0.0039,置信区间为[0.0026,0.0081]。

泊球测定的μ估计值为6.6640,置信区间为[6.6611,6.6669]; σ的估计值为0.0030,置信区间为[0.0019,0.0071]。 命令 利用mle函数进行参数估计 函数 mle

格式 phat=mle(?dist?,X) %返回用dist指定分布的最大似然估计值

[phat, pci]=mle(?dist?,X) %置信度为95% [phat, pci]=mle(?dist?,X,alpha) %置信度由alpha确定

[phat, pci]=mle(?dist?,X,alpha,pl) %仅用于二项分布,pl为试验次数。

说明 dist为分布函数名,如:beta(?分布)、bino(二项分布)等,X为数据样本,alpha为显著水平α,(1??)?100%为置信度。

例4-64

>> X=binornd(20,0.75) %产生二项分布的随机数 X = 16

>> [p,pci]=mle('bino',X,0.05,20) %求概率的估计值和置信区间,置信度为95% p =

0.8000 pci =

0.5634 0.9427

常用分布的参数估计函数

表4-7 参数估计函数表

调 用 形 式 PHAT= binofit(X, N) binofit [PHAT, PCI] = binofit(X,N) [PHAT, PCI]= binofit (X, N, ALPHA) Lambdahat=poissfit(X) poissfit [Lambdahat, Lambdaci] = poissfit(X) [Lambdahat, Lambdaci]= poissfit (X, ALPHA) [muhat,sigmahat,muci,sigmaci] = normfit(X) normfit [muhat,sigmahat,muci,sigmaci] = normfit(X, ALPHA) PHAT =betafit (X) betafit [PHAT, PCI]= betafit (X, ALPHA) [ahat,bhat] = unifit(X) unifit [ahat,bhat,ACI,BCI] = unifit(X) [ahat,bhat,ACI,BCI]=unifit(X, ALPHA) muhat =expfit(X) expfit [muhat,muci] = expfit(X) [muhat,muci] = expfit(X,alpha) gamfit phat =gamfit(X) 160 函数名 函 数 说 明 二项分布的概率的最大似然估计 置信度为95%的参数估计和置信区间 返回水平α的参数估计和置信区间 泊松分布的参数的最大似然估计 置信度为95%的参数估计和置信区间 返回水平α的λ参数和置信区间 正态分布的最大似然估计,置信度为95% 返回水平α的期望、方差值和置信区间 返回β分布参数a和 b的最大似然估计 返回最大似然估计值和水平α的置信区间 均匀分布参数的最大似然估计 置信度为95%的参数估计和置信区间 返回水平α的参数估计和置信区间 指数分布参数的最大似然估计 置信度为95%的参数估计和置信区间 返回水平α的参数估计和置信区间 γ分布参数的最大似然估计 第4章 概率统计

[phat,pci] = gamfit(X) [phat,pci] = gamfit(X,alpha) phat = weibfit(X) weibfit [phat,pci] = weibfit(X) [phat,pci] = weibfit(X,alpha) phat = mle('dist',data) [phat,pci] = mle('dist',data) Mle [phat,pci] = mle('dist',data,alpha) [phat,pci] = mle('dist',data,alpha,p1) 置信度为95%的参数估计和置信区间 返回最大似然估计值和水平α的置信区间 韦伯分布参数的最大似然估计 置信度为95%的参数估计和置信区间 返回水平α的参数估计及其区间估计 分布函数名为dist的最大似然估计 置信度为95%的参数估计和置信区间 返回水平α的最大似然估计值和置信区间 仅用于二项分布,pl为试验总次数 说明 各函数返回已给数据向量X的参数最大似然估计值和置信度为(1-α)×100%的置信区间。α的默认值为0.05,即置信度为95%。

4.7.2 非线性模型置信区间预测

命令 高斯—牛顿法的非线性最小二乘数据拟合 函数 nlinfit

格式 beta = nlinfit(X,y,FUN,beta0) %返回在FUN中描述的非线性函数的系数。FUN

??f(?,X)的函数,为用户提供形如y该函数返回已给初始参数估计

?。 值β和自变量X的y的预测值y[beta,r,J] = nlinfit(X,y,FUN,beta0) ?ta为拟合系数,r为残差,J为Jacobi矩

阵,beta0为初始预测值。

说明 若X为矩阵,则X的每一列为自变量的取值,y是一个相应的列向量。如果FUN中使用了@,则表示函数的柄。

例4-65 调用MATLAB提供的数据文件reaction.mat

>>load reaction

>>betafit = nlinfit(reactants,rate,@hougen,beta) betafit = 1.2526 0.0628 0.0400 0.1124 1.1914

命令 非线性模型的参数估计的置信区间 函数 nlparci

格式 ci = nlparci(beta,r,J) %返回置信度为95%的置信区间,beta为非线性最小二乘

法估计的参数值,r为残差,J为Jacobian矩阵。nlparci可以用nlinfit函数的输出作为其输入。

例4-66 调用MATLAB中的数据reaction。

>>load reaction

>>[beta,resids,J] = nlinfit(reactants,rate,'hougen',beta) beta =

1.2526 0.0628 0.0400 0.1124 1.1914 resids = 0.1321 -0.1642 -0.0909

161 MATLAB6.0数学手册 0.0310 0.1142 0.0498 -0.0262 0.3115 -0.0292 0.1096 0.0716 -0.1501 -0.3026 J =

6.8739 -90.6536 -57.8640 -1.9288 0.1614 3.4454 -48.5357 -13.6240 -1.7030 0.3034 5.3563 -41.2099 -26.3042 -10.5217 1.5095 1.6950 0.1091 0.0186 0.0279 1.7913 2.2967 -35.5658 -6.0537 -0.7567 0.2023 11.8670 -89.5655 -170.1745 -8.9566 0.4400 4.4973 -14.4262 -11.5409 -9.3770 2.5744 4.1831 -41.7896 -16.8937 -5.7794 1.0082 11.8286 -51.3721 -154.1164 -27.7410 1.5001 9.1514 -25.5948 -76.7844 -30.7138 2.5790 3.3373 0.0900 0.0720 0.1080 3.5269 9.3663 -102.0611 -107.4327 -3.5811 0.2200 4.7512 -24.4631 -16.3087 -10.3002 2.1141 >>ci = nlparci(beta,resids,J) ci =

-0.7467 3.2519 -0.0377 0.1632 -0.0312 0.1113 -0.0609 0.2857 -0.7381 3.1208

命令 非线性拟合和显示交互图形 函数 nlintool

格式 nlintool(x,y,FUN,beta0) %返回数据(x,y)的非线性曲线的预测图形,它用2条红

色曲线预测全局置信区间。beta0为参数的初始预测值,置信度为95%。

nlintool(x,y,FUN,beta0,alpha) %置信度为(1-alpha)×100% 例4-67 调用MATLAB数据

>> load reaction

>> nlintool(reactants,rate,'hougen',beta)

图4-20

命令 非线性模型置信区间预测 函数 nlpredci

162 第4章 概率统计

格式 ypred = nlpredci(FUN,inputs,beta,r,J) % ypred 为预测值,FUN与前面相同,beta

为给出的适当参数,r为残差,J为Jacobian矩阵,inputs为非线性函数中的独立变量的矩阵值。

[ypred,delta] = nlpredci(FUN,inputs,beta,r,J) Tlta为非线性最小二乘法估计

的置信区间长度的一半,当r长度超过beta的长度并且J的列满秩时,置信区间的计算是有效的。[ypred-delta,ypred+delta]为置信度为95%的不同步置信区间。

ypred = nlpredci(FUN,inputs,beta,r,J,alpha,'simopt','predopt') %控制置信区间的

类型,置信度为100(1-alpha)%。'simopt' = 'on' 或'off' (默认值)分别表示同步或不同步置信区间。'predopt'='curve' (默认值) 表示输入函数值的置信区间, 'predopt'='observation' 表示新响应值的置信区间。nlpredci可以用nlinfit函数的输出作为其输入。

例4-68 续前例,在[100 300 80]处的预测函数值ypred和置信区间一半宽度delta

>> load reaction

>> [beta,resids,J] = nlinfit(reactants,rate,@hougen,beta);

>> [ypred,delta] = nlpredci(@hougen,[100 300 80],beta,resids,J)

结果为:

ypred = 10.9113 delta = 0.3195

命令 非负最小二乘

函数 nnls(该函数已被函数lsnonneg代替,在6.0版中使用nnls将产生警告信息) 格式 x = nnls(A,b) %最小二乘法判断方程A×x=b的解,返回在x≥0的条件下使得

||A?x?b||最小的向量x,其中A和b必须为实矩阵或向量。 x = nnls(A,b,tol) % tol为指定的误差

[x,w] = nnls(A,b) %当x中元素xi?0时,wi?0,当xi?0时wi?0。 [x,w] = nnls(A,b,tol) 例4- 69

>> A =[0.0372 0.2869;0.6861 0.7071;0.6233 0.6245;0.6344 0.6170]; >> b=[0.8587 0.1781 0.0747 0.8405]'; >> x=nnls(A,b)

Warning: NNLS is obsolete and has been replaced by LSQNONNEG. NNLS now calls LSQNONNEG which uses the following syntax: [X,RESNORM,RESIDUAL,EXITFLAG,OUTPUT,LAMBDA] =lsqnonneg(A,b,X0, Options) ;

Use OPTIMSET to define optimization options, or type 'edit nnls' to view the code used here. NNLS will be

removed in the future; please use NNLS with the new syntax. x =

0 0.6929

命令 有非负限制的最小二乘 函数 lsqnonneg

163

MATLAB6.0数学手册 格式 x = lsqnonneg(C,d) %返回在x≥0的条件下使得||C?x?d||最小的向量x,其中C和d必须为实矩阵或向量。

x = lsqnonneg(C,d,x0) % x0为初始点,x0≥0

x = lsqnonneg(C,d,x0,options) %options为指定的优化参数,参见options函数。 [x,resnorm] = lsqnonneg(?) %resnorm表示norm(C*x-d).^2的残差 [x,resnorm,residual] = lsqnonneg(?) %residual表示C*x-d的残差 例4- 70

>> A =[0.0372 0.2869;0.6861 0.7071;0.6233 0.6245;0.6344 0.6170]; >> b=[0.8587 0.1781 0.0747 0.8405]'; >> [x,resnorm,residual] = lsqnonneg(A,b) x =

0 0.6929 resnorm = 0.8315 residual = 0.6599 -0.3119 -0.3580 0.4130

4.7.3 对数似然函数

命令 负?分布的对数似然函数 函数 Betalike

格式 logL=betalike(params,data) %返回负?分布的对数似然函数,params为向量[a,

b],是?分布的参数,data为样本数据。

[logL,info]=betalike(params,data) %返回Fisher逆信息矩阵info。如果params 中

输入的参数是极大似然估计值,那么info的对角元素为相应参数的渐近方差。 说明 betalike是?分布最大似然估计的实用函数。似然函数假设数据样本中,所有的元素相互独立。因为betalike返回负?对数似然函数,用fmins函数最小化betalike与最大似然估计的功能是相同的。

例4-71 本例所取的数据是随机产生的?分布数据。

>>r = betarnd(3,3,100,1);

>>[logL,info] = betalike([2.1234,3.4567],r) logL =

-12.4340 info =

0.1185 0.1364 0.1364 0.2061

命令 负?分布的对数似然估计 函数 Gamlike

格式 logL=gamlike(params,data) %返回由给定样本数据data确定的?分布的参数为params(即[a,b])的负对数似然函数值

[logL,info]=gamlike(params,data) %返回Fisher逆信息矩阵info。如果params

164 第4章 概率统计

中输入的参数是极大似然估计值,那么info的对角元素为相应参数的渐近方差。

说明 gamlike是?分布的最大似然估计函数。因为gamlike返回?对数似然函数值,故用fmins函数将gamlike最小化后,其结果与最大似然估计是相同的。 例4-72

>>r=gamrnd(2,3,100,1);

>>[logL,info]=gamlike([2.4212, 2.5320],r) logL =

275.4602 info =

0.0453 -0.0538 -0.0538 0.0867

命令 负正态分布的对数似然函数 函数 normlike

格式 logL=normlike(params,data) %返回由给定样本数据data确定的、负正态分布的、

参数为params(即[mu,sigma])的对数似然函数值。

[logL,info]=normlike(params,data) %返回Fisher逆信息矩阵info。如果params

中输入的参数是极大似然估计值,那么info的对角元素为相应参数的渐近方差。

命令 威布尔分布的对数似然函数 函数 Weiblike

格式 logL = weiblike(params,data) %返回由给定样本数据data确定的、威布尔分布

的、参数为params(即[a,b])的对数似然函数值。

[logL,info]=weiblike(params,data) %返回Fisher逆信息矩阵info。如果params中

输入的参数是极大似然估计值,那么info的对角元素为相应参数的渐近方差。

说明 威布尔分布的负对数似然函数定义为

?logL??log?f(a,b|xi)???logf(a,b|xi)

i?1i?1nn例4-73

>>r=weibrnd(0.4,0.98,100,1);

>>[logL,info]=weiblike([0.1342,0.9876],r) logL =

237.6682 info =

0.0004 -0.0002 -0.0002 0.0078

4.8 假设检验

4.8.1 ?2已知,单个正态总体的均值μ的假设检验(U检验法)

函数 ztest

格式 h = ztest(x,m,sigma) % x为正态总体的样本,m为均值μ0,sigma为标准差,

165 MATLAB6.0数学手册 显著性水平为0.05(默认值)

h = ztest(x,m,sigma,alpha) %显著性水平为alpha

[h,sig,ci,zval] = ztest(x,m,sigma,alpha,tail) %sig为观察值的概率,当sig为小概

率时则对原假设提出质疑,ci为真正均值μ的1-alpha置信区间,zval为统计量的值。

说明 若h=0,表示在显著性水平alpha下,不能拒绝原假设; 若h=1,表示在显著性水平alpha下,可以拒绝原假设。 原假设:H0:???0?m,

若tail=0,表示备择假设:H1:???0?m(默认,双边检验); tail=1,表示备择假设:H1:???0?m(单边检验);

tail=-1,表示备择假设:H1:???0?m(单边检验)。 例4-74 某车间用一台包装机包装葡萄糖,包得的袋装糖重是一个随机变量,它服从正态分布。当机器正常时,其均值为0.5公斤,标准差为0.015。某日开工后检验包装机是否正常,随机地抽取所包装的糖9袋,称得净重为(公斤)

0.497, 0.506, 0.518, 0.524, 0.498, 0.511, 0.52, 0.515, 0.512

问机器是否正常?

解:总体μ和σ已知,该问题是当?2为已知时,在水平??0.05下,根据样本值判断μ=0.5还是??0.5。为此提出假设:

原假设: H0:???0?0.5 备择假设:H1:??0.5

>> X=[0.497,0.506,0.518,0.524,0.498,0.511,0.52,0.515,0.512]; >> [h,sig,ci,zval]=ztest(X,0.5,0.015,0.05,0)

结果显示为

h = 1 sig =

0.0248 %样本观察值的概率 ci =

0.5014 0.5210 %置信区间,均值0.5在此区间之外 zval =

2.2444 %统计量的值

结果表明:h=1,说明在水平??0.05下,可拒绝原假设,即认为包装机工作不正常。

4.8.2 ?2未知,单个正态总体的均值μ的假设检验( t检验法)

函数 ttest

格式 h = ttest(x,m) % x为正态总体的样本,m为均值μ0,显著性水平为0.05

h = ttest(x,m,alpha) %alpha为给定显著性水平

[h,sig,ci] = ttest(x,m,alpha,tail) %sig为观察值的概率,当sig为小概率时则对原

假设提出质疑,ci为真正均值μ的1-alpha置信区间。

说明 若h=0,表示在显著性水平alpha下,不能拒绝原假设; 若h=1,表示在显著性水平alpha下,可以拒绝原假设。 原假设:H0:???0?m,

166 第4章 概率统计

若 tail=0,表示备择假设:H1:???0?m(默认,双边检验); tail=1,表示备择假设:H1:???0?m(单边检验); tail=-1,表示备择假设:H1:???0?m(单边检验)。

例4-75 某种电子元件的寿命X(以小时计)服从正态分布,?、σ2均未知。现测得16只元件的寿命如下

159 280 101 212 224 379 179 264 222 362 168 250 149 260 485 170

问是否有理由认为元件的平均寿命大于225(小时)?

解:未知?2,在水平??0.05下检验假设:H0:???0?225,H1:??225

>> X=[159 280 101 212 224 379 179 264 222 362 168 250 149 260 485 170]; >> [h,sig,ci]=ttest(X,225,0.05,1)

结果显示为:

h = 0 sig =

0.2570 ci =

198.2321 Inf %均值225在该置信区间内

结果表明:H=0表示在水平??0.05下应该接受原假设H0,即认为元件的平均寿命不大于225小时。

4.8.3 两个正态总体均值差的检验(t检验)

两个正态总体方差未知但等方差时,比较两正态总体样本均值的假设检验 函数 ttest2

格式 [h,sig,ci]=ttest2(X,Y) %X,Y为两个正态总体的样本,显著性水平为0.05 [h,sig,ci]=ttest2(X,Y,alpha) %alpha为显著性水平

[h,sig,ci]=ttest2(X,Y,alpha,tail) %sig为当原假设为真时得到观察值的概率,当

sig为小概率时则对原假设提出质疑,ci为真正均值μ的1-alpha置信区间。

说明 若h=0,表示在显著性水平alpha下,不能拒绝原假设; 若h=1,表示在显著性水平alpha下,可以拒绝原假设。

原假设:H0:?1??2, (?1为X为期望值,?2为Y的期望值)

若 tail=0,表示备择假设:H1:?1??2(默认,双边检验); tail=1,表示备择假设:H1:?1??2(单边检验); tail=-1,表示备择假设:H1:?1??2(单边检验)。

例4-76 在平炉上进行一项试验以确定改变操作方法的建议是否会增加钢的产率,试验是在同一只平炉上进行的。每炼一炉钢时除操作方法外,其他条件都尽可能做到相同。先用标准方法炼一炉,然后用建议的新方法炼一炉,以后交替进行,各炼10炉,其产率分别为

(1)标准方法:78.1 72.4 76.2 74.3 77.4 78.4 76.0 75.5 76.7 77.3 (2)新方法: 79.1 81.0 77.3 79.1 80.0 79.1 79.1 77.3 80.2 82.1

设这两个样本相互独立,且分别来自正态总体N(?1,?2)和N(?2,?2),?1、?2、?2均未知。问建议的新操作方法能否提高产率?(取α=0.05)

167 MATLAB6.0数学手册 解:两个总体方差不变时,在水平??0.05下检验假设:H0:?1??2,H1:?1??2

>> X=[78.1 72.4 76.2 74.3 77.4 78.4 76.0 75.5 76.7 77.3]; >>Y=[79.1 81.0 77.3 79.1 80.0 79.1 79.1 77.3 80.2 82.1]; >> [h,sig,ci]=ttest2(X,Y,0.05,-1)

结果显示为:

h = 1 sig =

2.1759e-004 %说明两个总体均值相等的概率很小 ci =

-Inf -1.9083

结果表明:H=1表示在水平??0.05下,应该拒绝原假设,即认为建议的新操作方法提高了产率,因此,比原方法好。

4.8.4 两个总体一致性的检验——秩和检验

函数 ranksum

格式 p = ranksum(x,y,alpha) %x、y为两个总体的样本,可以不等长,alpha为显著

性水平

[p,h] = ranksum(x,y,alpha) % h为检验结果,h=0表示X与Y的总体差别不显

著h=1表示X与Y的总体差别显著

[p,h,stats] = ranksum(x,y,alpha) %stats中包括:ranksum为秩和统计量的值以及

zval为过去计算p的正态统计量的值

说明 P为两个总体样本X和Y为一致的显著性概率,若P接近于0,则不一致较明显。 例4-77 某商店为了确定向公司A或公司B购买某种商品,将A和B公司以往的各次进货的次品率进行比较,数据如下所示,设两样本独立。问两公司的商品的质量有无显著差异。设两公司的商品的次品的密度最多只差一个平移,取α=0.05。

A:7.0 3.5 9.6 8.1 6.2 5.1 10.4 4.0 2.0 10.5

B:5.7 3.2 4.1 11.0 9.7 6.9 3.6 4.8 5.6 8.4 10.1 5.5 12.3 解:设?A、?B分别为A、B两个公司的商品次品率总体的均值。则该问题为在水平α=0.05下检验假设:H0:?A??B,H1:?A??B

>> A=[7.0 3.5 9.6 8.1 6.2 5.1 10.4 4.0 2.0 10.5];

>> B=[5.7 3.2 4.1 11.0 9.7 6.9 3.6 4.8 5.6 8.4 10.1 5.5 12.3]; >> [p,h,stats]=ranksum(A,B,0.05)

结果为:

p =

0.8041 h = 0 stats =

zval: -0.2481 ranksum: 116

结果表明:一方面,两样本总体均值相等的概率为0.8041,不接近于0;另一方面,H=0也说明可以接受原假设H0,即认为两个公司的商品的质量无明显差异。

4.8.5 两个总体中位数相等的假设检验——符号秩检验

函数 signrank

168

第4章 概率统计

格式 p = signrank(X,Y,alpha) % X、Y为两个总体的样本,长度必须相同,alpha为

显著性水平,P两个样本X和Y的中位数相等的概率,p接近于0则可对原假设质疑。

[p,h] = signrank(X,Y,alpha) % h为检验结果:h=0表示X与Y的中位数之差不

显著,h=1表示X与Y的中位数之差显著。

[p,h,stats] = signrank(x,y,alpha) % stats中包括:signrank为符号秩统计量的值以

及zval为过去计算p的正态统计量的值。

例4-78 两个正态随机样本的中位数相等的假设检验

>> x=normrnd(0,1,20,1); >> y=normrnd(0,2,20,1);

>> [p,h,stats]=signrank(x,y,0.05) p =

0.3703 h = 0 stats =

zval: -0.8960 signedrank: 81

结果表明:h=0表示X与Y的中位数之差不显著

4.8.6 两个总体中位数相等的假设检验——符号检验

函数 signtest

格式 p=signtest(X, Y, alpha) % X、Y为两个总体的样本,长度必须相同,alpha为显

著性水平,P两个样本X和Y的中位数相等的概率,p接近于0则可对原假设质疑。

[p, h]=signtest(X, Y, alpha) % h为检验结果:h=0表示X与Y的中位数之差不显

著,h=1表示X与Y的中位数之差显著。

[p,h,stats] = signtest(X,Y,alpha) % stats中sign为符号统计量的值 例4-79 两个正态随机样本的中位数相等的假设检验

>> X=normrnd(0,1,20,1); >> Y=normrnd(0,2,20,1);

>> [p,h,stats]=signtest(X,Y,0.05) p =

0.2632 h = 0 stats = sign: 7

结果表明:h=0表示X与Y的中位数之差不显著

4.8.7 正态分布的拟合优度测试

函数 jbtest

格式 H = jbtest(X) %对输入向量X进行Jarque-Bera测试,显著性水平为0.05。

H = jbtest(X,alpha) %在水平alpha而非5%下施行 Jarque-Bera 测试,alpha在

0和1之间。

[H,P,JBSTAT,CV] = jbtest(X,alpha) %P为接受假设的概率值,P越接近于0,则

169 MATLAB6.0数学手册 可以拒绝是正态分布的原假设;JBSTAT为测试统计量的值,CV为是否拒绝原假设的临界值。

说明 H为测试结果,若H=0,则可以认为X是服从正态分布的;若X=1,则可以否定X服从正态分布。X为大样本,对于小样本用lillietest函数。

例4-80 调用MATLAB中关于汽车重量的数据,测试该数据是否服从正态分布

>> load carsmall

>> [h,p,j,cv]=jbtest(Weight) h = 1 p =

0.0267 j =

7.2448 cv =

5.9915

说明 p=2.67%表示应该拒绝服从正态分布的假设;h=1也可否定服从正态分布;统计量的值j = 7.2448大于接受假设的临界值cv =5.9915,因而拒绝假设(测试水平为5%)。

4.8.8 正态分布的拟合优度测试

函数 lillietest

格式 H = lillietest(X) %对输入向量X进行Lilliefors测试,显著性水平为0.05。

H = lillietest(X,alpha) %在水平alpha而非5%下施行Lilliefors测试,alpha在

0.01和0.2之间。

[H,P,LSTAT,CV] = lillietest(X,alpha) %P为接受假设的概率值,P越接近于0,

则可以拒绝是正态分布的原假设;LSTAT为测试统计量的值,CV为是否拒绝原假设的临界值。

说明 H为测试结果,若H=0,则可以认为X是服从正态分布的;若X=1,则可以否定X服从正态分布。

例4-81

>> Y=chi2rnd(10,100,1); >> [h,p,l,cv]=lillietest(Y) h = 1 p =

0.0175 l =

0.1062 cv =

0.0886

图4-21

说明 h=1表示拒绝正态分布的假设;p = 0.0175表示服从正态分布的概率很小;统计量的值l = 0.1062大于接受假设的临界值cv =0.0886,因而拒绝假设(测试水平为5%)。

>>hist(Y)

从图中看出,数据Y不服从正态分布。

4.8.9 单个样本分布的 Kolmogorov-Smirnov 测试

函数 kstest

170 第4章 概率统计

格式 H = kstest(X) %测试向量X是否服从标准正态分布,测试水平为5%。

H = kstest(X,cdf) %指定累积分布函数为cdf的测试(cdf=[ ]时表示标准正态分

布),测试水平为5%

H = kstest(X,cdf,alpha) % alpha为指定测试水平

[H,P,KSSTAT,CV] = kstest(X,cdf,alpha) %P为原假设成立的概率,KSSTAT为测

试统计量的值,CV为是否接受假设的临界值。

说明 原假设为X服从标准正态分布。若H=0则不能拒绝原假设,H=1则可以拒绝原假设。

例4-82 产生100个威布尔随机数,测试该随机数服从的分布

>> x=weibrnd(1,2,100,1);

>> [H,p,ksstat,cv]=kstest(x,[x weibcdf(x,1,2)],0.05) %测试是否服从威布尔分布 H = 0 p =

0.3022 ksstat = 0.0959 cv =

0.1340

说明 H=0表示接受原假设,统计量ksstat小于临界值表示接受原假设。 >> [H,p,ksstat,cv]=kstest(x,[x expcdf(x,1)],0.05) %测试是否服从指数分布

H = 1 p =

0.0073 ksstat = 0.1653 cv =

0.1340

说明 H=1表明拒绝服从指数分布的假设。

>> [H,p,ksstat,cv]=kstest(x,[ ],0.05) %测试是否服从标准正态分布 H = 1 p =

3.1285e-026 ksstat = 0.5380 cv =

0.1340

说明 H=1表明不服从标准正态分布。

4.8.10 两个样本具有相同的连续分布的假设检验

函数 kstest2

格式 H = kstest2(X1,X2) %测试向量X1与X2是具有相同的连续分布,测试水

平为5%。

H = kstest2(X1,X2,alpha) % alpha为测试水平

[H,P,KSSTAT] = kstest(X,cdf,alpha) %与指定累积分布cdf相同的连续分布,P

为假设成立的概率,KSSTAT为测试统计量的值。

说明 原假设为具有相同连续分布。测试结果为H,若H=0,表示应接受原假设;若

171 MATLAB6.0数学手册 H=1,表示可以拒绝原假设。这是Kolmogorov-Smirnov测试方法。

例4-83

>> x=-1:1:5;

>> y=randn(20,1); >> [h,p,k]=kstest2(x,y) h = 1 p =

0.0444 k =

0.5643

说明 h=1表示可以认为向量x与y的分布不相同,相同的概率只有4.4%。

4.9 方差分析

4.9.1 单因素方差分析

单因素方差分析是比较两组或多组数据的均值,它返回原假设——均值相等的概率 函数 anova1

格式 p = anova1(X) %X的各列为彼此独立的样本观察值,其元素个数相同,p为各

列均值相等的概率值,若p值接近于0,则原假设受到怀疑,说明至少有一列均值与其余列均值有明显不同。

p = anova1(X,group) %X和group为向量且group要与X对应

p = anova1(X,group,'displayopt') % displayopt=on/off表示显示与隐藏方差分析

表图和盒图

[p,table] = anova1(?) % table为方差分析表

[p,table,stats] = anova1(?) % stats为分析结果的构造 说明 anova1函数产生两个图:标准的方差分析表图和盒图。

方差分析表中有6列:第1列(source)显示:X中数据可变性的来源;第2列(SS)显示:用于每一列的平方和;第3列(df)显示:与每一种可变性来源有关的自由度;第4列(MS)显示:是SS/df的比值;第5列(F)显示:F统计量数值,它是MS的比率;第6列显示:从F累积分布中得到的概率,当F增加时,p值减少。

例4-84 设有3台机器,用来生产规格相同的铝合金薄板。取样测量薄板的厚度,精确至?厘米。得结果如下:

机器1:0.236 0.238 0.248 0.245 0.243 机器2:0.257 0.253 0.255 0.254 0.261 机器3:0.258 0.264 0.259 0.267 0.262

检验各台机器所生产的薄板的厚度有无显著的差异? 解:

>> X=[0.236 0.238 0.248 0.245 0.243; 0.257 0.253 0.255 0.254 0.261;?

0.258 0.264 0.259 0.267 0.262]; >> P=anova1(X')

172 第4章 概率统计

结果为:

P =

1.3431e-005

还有两个图,即图4-22和图4-23。

0.2650.260.255Values0.250.2450.240.235

12Column Number3

图4-22 图4-23

例4-85 建筑横梁强度的研究:3000磅力量作用在一英寸的横梁上来测量横梁的挠度,钢筋横梁的测试强度是:82 86 79 83 84 85 86 87;其余两种更贵的合金横梁强度测试为合金1:74 82 78 75 76 77;合金2:79 79 77 78 82 79]。

检验这些合金强度有无明显差异? 解:

>> strength = [82 86 79 83 84 85 86 87 74 82 78 75 76 77 79 79 77 78 82 79]; >>alloy = {'st','st','st','st','st','st','st','st', 'al1','al1','al1','al1','al1','al1',? 'al2','al2','al2','al2','al2','al2'};

>> [p,table,stats] = anova1(strength,alloy,'on')

结果为

p =

1.5264e-004 table =

'Source' 'SS' 'df' 'MS' 'F' 'Prob>F' 'Groups' [184.8000] [ 2] [92.4000] [15.4000] [1.5264e-004] 'Error' [102.0000] [17] [ 6.0000] [ ] [ ] 'Total' [286.8000] [19] [ ] [ ] [ ] stats =

gnames: {3x1 cell} n: [8 6 6] source: 'anova1' means: [84 77 79] df: 17 s: 2.4495

868482Values80787674

stal1al2

图4-24 图4-25

173

MATLAB6.0数学手册 说明 p值显示,3种合金是明显不同的,盒图显示钢横梁的挠度大于另两种合金横梁的挠度。

4.9.2 双因素方差分析

函数 anova2

格式 p = anova2(X,reps)

p = anova2(X,reps,'displayopt') [p,table] = anova2(?)

[p,table,stats] = anova2(…)

说明 执行平衡的双因素试验的方差分析来比较X中两个或多个列(行)的均值,不同列的数据表示因素A的差异,不同行的数据表示另一因素B的差异。如果行列对有多于一个的观察点,则变量reps指出每一单元观察点的数目,每一单元包含reps行,如:

A?1A?2?x111?x?121?x211??x221?x311??x321x112??B?1x122??x212???B?2 x222?x312???B?3x322?reps=2

其余参数与单因素方差分析参数相似。

例4-86 一火箭使用了4种燃料,3种推进器作射程试验,每种燃料与每种推进器的组合各发射火箭2次,得到结果如下:

推进器(B) B1 B2 B3 A1 58.2000 56.2000 65.3000 52.6000 41.2000 60.8000 A2 49.1000 54.1000 51.6000 燃料A 42.8000 50.5000 48.4000 A3 60.1000 70.9000 39.2000 58.3000 73.2000 40.7000 A4 75.8000 58.2000 48.7000 71.5000 51.0000 41.4000

考察推进器和燃料这两个因素对射程是否有显著的影响? 解:建立M文件

X=[58.2000 56.2000 65.3000 52.6000 41.2000 60.8000 49.1000 54.1000 51.6000 42.8000 50.5000 48.4000 60.1000 70.9000 39.2000 58.3000 73.2000 40.7000 75.8000 58.2000 48.7000 71.5000 51.0000 41.4000]; P=anova2(X,2)

结果为:

174 第4章 概率统计

P =

0.0035 0.0260 0.0001

显示方差分析图为图4-26。

图4-26

175

本文来源:https://www.bwwdw.com/article/wre5.html

Top