4.6.6 样本数据的盒图
函数 boxplot
格式 boxplot(X) %产生矩阵X的每一列的盒图和"须"图,"须"是从盒的尾部延伸出来,并表示盒外数据长度的线,如果"须"的外面没有数据,则在"须"的底部有一个点.
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
格式 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.6.11 在指定的界线之间画正态密度曲线
函数 normspec
格式 p = normspec(specs,mu,sigma) %specs指定界线,mu,sigma为正态分布的参数p 为样本落在上,下界之间的概率.
例4-60
>>normspec([10 Inf],11.5,1.25)
图4-19 
 

admin 2007-11-29 20:50
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
说明 估计值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分别为置信区间,其置信度为;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
设测定值总体为,μ和σ为未知.对(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 =
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指定分布的最大似然估计值
[phat, pci]=mle %置信度为95%
[phat, pci]=mle %置信度由alpha确定
[phat, pci]=mle %仅用于二项分布,pl为试验次数.
说明 dist为分布函数名,如:beta(分布),bino(二项分布)等,X为数据样本,alpha为显著水平α,为置信度.
例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 参数估计函数表
函数名
调 用 形 式
函 数 说 明
binofit
PHAT= binofit(X, N)
[PHAT, PCI] = binofit(X,N)
[PHAT, PCI]= binofit (X, N, ALPHA)
二项分布的概率的最大似然估计
置信度为95%的参数估计和置信区间
返回水平α的参数估计和置信区间
poissfit
Lambdahat=poissfit(X)
[Lambdahat, Lambdaci] = poissfit(X)
[Lambdahat, Lambdaci]= poissfit (X, ALPHA)
泊松分布的参数的最大似然估计
置信度为95%的参数估计和置信区间
返回水平α的λ参数和置信区间
normfit
[muhat,sigmahat,muci,sigmaci] = normfit(X)
[muhat,sigmahat,muci,sigmaci] = normfit(X, ALPHA)
正态分布的最大似然估计,置信度为95%
返回水平α的期望,方差值和置信区间
betafit
PHAT =betafit (X)
[PHAT, PCI]= betafit (X, ALPHA)
返回β分布参数a和 b的最大似然估计
返回最大似然估计值和水平α的置信区间
unifit
[ahat,bhat] = unifit(X)
[ahat,bhat,ACI,BCI] = unifit(X)
[ahat,bhat,ACI,BCI]=unifit(X, ALPHA)
均匀分布参数的最大似然估计
置信度为95%的参数估计和置信区间
返回水平α的参数估计和置信区间
expfit
muhat =expfit(X)
[muhat,muci] = expfit(X)
[muhat,muci] = expfit(X,alpha)
指数分布参数的最大似然估计
置信度为95%的参数估计和置信区间
返回水平α的参数估计和置信区间
gamfit
phat =gamfit(X)
[phat,pci] = gamfit(X)
[phat,pci] = gamfit(X,alpha)
γ分布参数的最大似然估计
置信度为95%的参数估计和置信区间
返回最大似然估计值和水平α的置信区间
weibfit
phat = weibfit(X)
[phat,pci] = weibfit(X)
[phat,pci] = weibfit(X,alpha)
韦伯分布参数的最大似然估计
置信度为95%的参数估计和置信区间
返回水平α的参数估计及其区间估计
Mle
phat = mle('dist',data)
[phat,pci] = mle('dist',data)
[phat,pci] = mle('dist',data,alpha)
[phat,pci] = mle('dist',data,alpha,p1)
分布函数名为dist的最大似然估计
置信度为95%的参数估计和置信区间
返回水平α的最大似然估计值和置信区间
仅用于二项分布,pl为试验总次数
说明 各函数返回已给数据向量X的参数最大似然估计值和置信度为(1-α)×100%的置信区间.α的默认值为0.05,即置信度为95%.
4.7.2 非线性模型置信区间预测
命令 高斯—牛顿法的非线性最小二乘数据拟合
函数 nlinfit
格式 beta = nlinfit(X,y,FUN,beta0) %返回在FUN中描述的非线性函数的系数.FUN为用户提供形如的函数,该函数返回已给初始参数估计值β和自变量X的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
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
格式 ypred = nlpredci(FUN,inputs,beta,r,J) % ypred 为预测值,FUN与前面相同,beta为给出的适当参数,r为残差,J为Jacobian矩阵,inputs为非线性函数中的独立变量的矩阵值.
[ypred,delta] = nlpredci(FUN,inputs,beta,r,J) �lta为非线性最小二乘法估计的置信区间长度的一半,当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的条件下使得最小的向量x,其中A和b必须为实矩阵或向量.
x = nnls(A,b,tol) % tol为指定的误差
[x,w] = nnls(A,b) %当x中元素时,,当时.
[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
格式 x = lsqnonneg(C,d) %返回在x≥0的条件下使得最小的向量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中输入的参数是极大似然估计值,那么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的对角元素为相应参数的渐近方差.
说明 威布尔分布的负对数似然函数定义为
例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 已知,单个正态总体的均值μ的假设检验(U检验法)
函数 ztest
格式 h = ztest(x,m,sigma) % x为正态总体的样本,m为均值μ0,sigma为标准差,显著性水平为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下,可以拒绝原假设.
原假设:,
若tail=0,表示备择假设:(默认,双边检验);
tail=1,表示备择假设:(单边检验);
tail=-1,表示备择假设:(单边检验).
例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
问机器是否正常
解:总体μ和σ已知,该问题是当为已知时,在水平下,根据样本值判断μ=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,说明在水平下,可拒绝原假设,即认为包装机工作不正常.
4.8.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下,可以拒绝原假设.
原假设:,
若 tail=0,表示备择假设:(默认,双边检验);
tail=1,表示备择假设:(单边检验);
tail=-1,表示备择假设:(单边检验).
例4-75 某种电子元件的寿命X(以小时计)服从正态分布,,σ2均未知.现测得16只元件的寿命如下
159 280 101 212 224 379 179 264 222 362 168 250
149 260 485 170
问是否有理由认为元件的平均寿命大于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表示在水平下应该接受原假设,即认为元件的平均寿命不大于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下,可以拒绝原假设.
原假设:, (为X为期望值,为Y的期望值)
若 tail=0,表示备择假设:(默认,双边检验);
tail=1,表示备择假设:(单边检验);
tail=-1,表示备择假设:(单边检验).
例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
设这两个样本相互独立,且分别来自正态总体和,,,均未知.问建议的新操作方法能否提高产率 (取α=0.05)
解:两个总体方差不变时,在水平下检验假设::,:
>> 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表示在水平下,应该拒绝原假设,即认为建议的新操作方法提高了产率,因此,比原方法好.
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两个公司的商品次品率总体的均值.则该问题为在水平α=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];
>> [p,h,stats]=ranksum(A,B,0.05)
结果为:
p =
0.8041
h =
0
stats =
zval: -0.2481
ranksum: 116
结果表明:一方面,两样本总体均值相等的概率为0.8041,不接近于0;另一方面,H=0也说明可以接受原假设,即认为两个公司的商品的质量无明显差异.
4.8.5 两个总体中位数相等的假设检验——符号秩检验
函数 signrank
格式 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,则可以拒绝是正态分布的原假设;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
说明 h=1表示拒绝正态分布的假设;p = 0.0175表示服从正态分布的概率很小;统计量的值l = 0.1062大于接受假设的临界值cv =0.0886,因而拒绝假设(测试水平为5%).
>>hist(Y)
从图中看出,数据Y不服从正态分布.
4.8.9 单个样本分布的 Kolmogorov-Smirnov 测试
函数 kstest
格式 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,表示应接受原假设;若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')
结果为:
P =
1.3431e-005
还有两个图,即图4-22和图4-23.

图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

图4-24 图4-25
说明 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行,如:
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)
结果为:
P =
0.0035 0.0260 0.0001
显示方差分析图为图4-26.
图4-26
Matlab6.0数学手册