Matlab概率统计编程指南

第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的正态分布的随机数据,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
unifrnd ( A,B,m,n)
[A,B]上均匀分布(连续) 随机数
Unidrnd
unidrnd(N,m,n)
均匀分布(离散)随机数
Exprnd
exprnd(Lambda,m,n)
参数为Lambda的指数分布随机数
Normrnd
normrnd(MU,SIGMA,m,n)
参数为MU,SIGMA的正态分布随机数
chi2rnd
chi2rnd(N,m,n)
自由度为N的卡方分布随机数
Trnd
trnd(N,m,n)
自由度为N的t分布随机数
Frnd
frnd(N1, N2,m,n)
第一自由度为N1,第二自由度为N2的F分布随机数
gamrnd
gamrnd(A, B,m,n)
参数为A, B的分布随机数
betarnd
betarnd(A, B,m,n)
参数为A, B的分布随机数
lognrnd
lognrnd(MU, SIGMA,m,n)
参数为MU, SIGMA的对数正态分布随机数
nbinrnd
nbinrnd(R, P,m,n)
参数为R,P的负二项式分布随机数
ncfrnd
ncfrnd(N1, N2, delta,m,n)
参数为N1,N2,delta的非中心F分布随机数
nctrnd
nctrnd(N, delta,m,n)
参数为N,delta的非中心t分布随机数
ncx2rnd
ncx2rnd(N, delta,m,n)
参数为N,delta的非中心卡方分布随机数
raylrnd
raylrnd(B,m,n)
参数为B的瑞利分布随机数
weibrnd
weibrnd(A, B,m,n)
参数为A, B的韦伯分布随机数
binornd
binornd(N,P,m,n)
参数为N, p的二项分布随机数
geornd
geornd(P,m,n)
参数为 p的几何分布随机数
hygernd
hygernd(M,K,N,m,n)
参数为 M,K,N的超几何分布随机数
Pois-s-rnd
pois-s-rnd(Lambda,m,n)
参数为Lambda的泊松分布随机数
4.1.4 通用函数求各分布的随机数据
命令 求指定分布的随机数
函数 random
格式 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 常见分布函数表
name的取值
函数说明
\'beta\'
或
\'Beta\'
Beta分布
\'bino\'
或
\'Binomial\'
二项分布
\'chi2\'
或
\'Chisquare\'
卡方分布
\'exp\'
或
\'Exponential\'
指数分布
\'f\'
或
\'F\'
F分布
\'gam\'
或
\'Gamma\'
GAMMA分布
\'geo\'
或
\'Geometric\'
几何分布
\'hyge\'
或
\'Hypergeometric\'
超几何分布
\'logn\'
或
\'Lognormal\'
对数正态分布
\'nbin\'
或
\'Negative Binomial\'
负二项式分布
\'ncf\'
或
\'Noncentral F\'
非中心F分布
\'nct\'
或
\'Noncentral t\'
非中心t分布
\'ncx2\'
或
\'Noncentral Chi-square\'
非中心卡方分布
\'norm\'
或
\'Normal\'
正态分布
\'poiss\'
或
\'Poisson\'
泊松分布
\'rayl\'
或
\'Rayleigh\'
瑞利分布
\'t\'
或
\'T\'
T分布
\'unif\'
或
\'Uniform\'
均匀分布
\'unid\'
或
\'Discrete Uniform\'
离散均匀分布
\'weib\'
或
\'Weibull\'
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的密度函数值.
解:>> 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) %等同于, p — 每次试验事件A发生的概率;K—事件A发生K次;n—试验总次数
命令 泊松分布的概率值
函数 poisspdf
格式 poisspdf(k, Lambda) %等同于
命令 正态分布的概率值
函数 normpdf(K,mu,sigma) %计算参数为μ=mu,σ=sigma的正态分布密度函数在K处的值
专用函数计算概率密度函数列表如表4-3.
 
 

admin 2007-11-29 20:43
表4-3 专用函数计算概率密度函数表
函数名
调用形式
注 释
Unifpdf
unifpdf (x, a, b)
[a,b]上均匀分布(连续)概率密度在X=x处的函数值
unidpdf
Unidpdf(x,n)
均匀分布(离散)概率密度函数值
Exppdf
exppdf(x, Lambda)
参数为Lambda的指数分布概率密度函数值
normpdf
normpdf(x, mu, sigma)
参数为mu,sigma的正态分布概率密度函数值
chi2pdf
chi2pdf(x, n)
自由度为n的卡方分布概率密度函数值
Tpdf
tpdf(x, n)
自由度为n的t分布概率密度函数值
Fpdf
fpdf(x, n1, n2)
第一自由度为n1,第二自由度为n2的F分布概率密度函数值
gampdf
gampdf(x, a, b)
参数为a, b的分布概率密度函数值
betapdf
betapdf(x, a, b)
参数为a, b的分布概率密度函数值
lognpdf
lognpdf(x, mu, sigma)
参数为mu, sigma的对数正态分布概率密度函数值
nbinpdf
nbinpdf(x, R, P)
参数为R,P的负二项式分布概率密度函数值
Ncfpdf
ncfpdf(x, n1, n2, delta)
参数为n1,n2,delta的非中心F分布概率密度函数值
Nctpdf
nctpdf(x, n, delta)
参数为n,delta的非中心t分布概率密度函数值
ncx2pdf
ncx2pdf(x, n, delta)
参数为n,delta的非中心卡方分布概率密度函数值
raylpdf
raylpdf(x, b)
参数为b的瑞利分布概率密度函数值
weibpdf
weibpdf(x, a, b)
参数为a, b的韦伯分布概率密度函数值
binopdf
binopdf(x,n,p)
参数为n, p的二项分布的概率密度函数值
geopdf
geopdf(x,p)
参数为 p的几何分布的概率密度函数值
hygepdf
hygepdf(x,M,K,N)
参数为 M,K,N的超几何分布的概率密度函数值
poisspdf
poisspdf(x,Lambda)
参数为Lambda的泊松分布的概率密度函数值
例4-6 绘制卡方分布密度函数在自由度分别为1,5,15的图形
>> x=0:0.1:30;
>> y1=chi2pdf(x,1); plot(x,y1,\':\')
>> hold on
>> 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,\'+\')
2.卡方分布
例4-8
>> x = 0:0.2:15;
>>y = chi2pdf(x,4);
>>plot(x,y)

图4-2
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)

图4-3
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,\'-\')

图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\'))

图4-5
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)

图4-6
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)

图4-7
13.T分布
例4-19
>>x = -5:0.1:5;
>>y = tpdf(x,5);
>>z = normpdf(x,0,1);
>>plot(x,y,\'-\',x,z,\'-.\')
14.威布尔分布
例4-20
>> t=0:0.1:3;
>> y=weibpdf(t,2,2);
>> plot(y)

图4-8 
 

admin 2007-11-29 20:43
4.3 随机变量的累积概率值(分布函数值)
4.3.1 通用函数计算累积概率值
命令 通用函数cdf用来计算随机变量的概率之和(累积概率值)
函数 cdf
格式
说明 返回以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 专用函数计算累积概率值(随机变量的概率之和)
命令 二项分布的累积概率值
函数 binocdf
格式 binocdf (k, n, p) %n为试验总次数,p为每次试验事件A发生的概率,k为n次试验中事件A发生的次数,该命令返回n次试验中事件A恰好发生k次的概率.
命令 正态分布的累积概率值
函数 normcdf
格式 normcdf() %返回F(x)=的值,mu,sigma为正态分布的两个参数
例4-23 设X~N(3, 22)
(1)求
(2)确定c,使得
解(1) p1=
p2=
p3=
p4=
则有:
>>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
专用函数计算累积概率值函数列表如表4-4.
表4-4 专用函数的累积概率值函数表
函数名
调用形式
注 释
unifcdf
unifcdf (x, a, b)
[a,b]上均匀分布(连续)累积分布函数值 F(x)=P{X≤x}
unidcdf
unidcdf(x,n)
均匀分布(离散)累积分布函数值 F(x)=P{X≤x}
expcdf
expcdf(x, Lambda)
参数为Lambda的指数分布累积分布函数值 F(x)=P{X≤x}
normcdf
normcdf(x, mu, sigma)
参数为mu,sigma的正态分布累积分布函数值 F(x)=P{X≤x}
chi2cdf
chi2cdf(x, n)
自由度为n的卡方分布累积分布函数值 F(x)=P{X≤x}
tcdf
tcdf(x, n)
自由度为n的t分布累积分布函数值 F(x)=P{X≤x}
fcdf
fcdf(x, n1, n2)
第一自由度为n1,第二自由度为n2的F分布累积分布函数值
gamcdf
gamcdf(x, a, b)
参数为a, b的分布累积分布函数值 F(x)=P{X≤x}
betacdf
betacdf(x, a, b)
参数为a, b的分布累积分布函数值 F(x)=P{X≤x}
logncdf
logncdf(x, mu, sigma)
参数为mu, sigma的对数正态分布累积分布函数值
nbincdf
nbincdf(x, R, P)
参数为R,P的负二项式分布概累积分布函数值 F(x)=P{X≤x}
ncfcdf
ncfcdf(x, n1, n2, delta)
参数为n1,n2,delta的非中心F分布累积分布函数值
nctcdf
nctcdf(x, n, delta)
参数为n,delta的非中心t分布累积分布函数值 F(x)=P{X≤x}
ncx2cdf
ncx2cdf(x, n, delta)
参数为n,delta的非中心卡方分布累积分布函数值
raylcdf
raylcdf(x, b)
参数为b的瑞利分布累积分布函数值 F(x)=P{X≤x}
weibcdf
weibcdf(x, a, b)
参数为a, b的韦伯分布累积分布函数值 F(x)=P{X≤x}
binocdf
binocdf(x,n,p)
参数为n, p的二项分布的累积分布函数值 F(x)=P{X≤x}
geocdf
geocdf(x,p)
参数为 p的几何分布的累积分布函数值 F(x)=P{X≤x}
hygecdf
hygecdf(x,M,K,N)
参数为 M,K,N的超几何分布的累积分布函数值
poisscdf
poisscdf(x,Lambda)
参数为Lambda的泊松分布的累积分布函数值 F(x)=P{X≤x}
说明 累积概率函数就是分布函数F(x)=P{X≤x}在x处的值. 
 

admin 2007-11-29 20:44
4.4 随机变量的逆累积分布函数
MATLAB中的逆累积分布函数是已知,求x.
逆累积分布函数值的计算有两种方法
4.4.1 通用函数计算逆累积分布函数值
命令 icdf 计算逆累积分布函数
格式
说明 返回分布为name,参数为,累积概率值为P的临界值,这里name与前面表4.1相同.
如果,则
例4-24 在标准正态分布表中,若已知=0.975,求x
解:>> x=icdf(\'norm\',0.975,0,1)
x =
1.9600
例4-25 在分布表中,若自由度为10,=0.975,求临界值Lambda.
解:因为表中给出的值满足,而逆累积分布函数icdf求满足的临界值.所以,这里的取为0.025,即
>> Lambda=icdf(\'chi2\',0.025,10)
Lambda =
3.2470
例4-26 在假设检验中,求临界值问题:
已知:,查自由度为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 设,确定c使得.
解:由得,=0.5,所以
>>X=norminv(0.5, 3, 2)
X=
3
关于常用临界值函数可查下表4-5.
表4-5 常用临界值函数表
函数名
调用形式
注 释
unifinv
x=unifinv (p, a, b)
均匀分布(连续)逆累积分布函数(P=P{X≤x},求x)
unidinv
x=unidinv (p,n)
均匀分布(离散)逆累积分布函数,x为临界值
expinv
x=expinv (p, Lambda)
指数分布逆累积分布函数
norminv
x=Norminv(x,mu,sigma)
正态分布逆累积分布函数
chi2inv
x=chi2inv (x, n)
卡方分布逆累积分布函数
tinv
x=tinv (x, n)
t分布累积分布函数
finv
x=finv (x, n1, n2)
F分布逆累积分布函数
gaminv
x=gaminv (x, a, b)
分布逆累积分布函数
betainv
x=betainv (x, a, b)
分布逆累积分布函数
logninv
x=logninv (x, mu, sigma)
对数正态分布逆累积分布函数
nbininv
x=nbininv (x, R, P)
负二项式分布逆累积分布函数
ncfinv
x=ncfinv (x, n1, n2, delta)
非中心F分布逆累积分布函数
nctinv
x=nctinv (x, n, delta)
非中心t分布逆累积分布函数
ncx2inv
x=ncx2inv (x, n, delta)
非中心卡方分布逆累积分布函数
raylinv
x=raylinv (x, b)
瑞利分布逆累积分布函数
weibinv
x=weibinv (x, a, b)
韦伯分布逆累积分布函数
binoinv
x=binoinv (x,n,p)
二项分布的逆累积分布函数
geoinv
x=geoinv (x,p)
几何分布的逆累积分布函数
hygeinv
x=hygeinv (x,M,K,N)
超几何分布的逆累积分布函数
poissinv
x=poissinv (x,Lambda)
泊松分布的逆累积分布函数
例4-28 公共汽车门的高度是按成年男子与车门顶碰头的机会不超过1%设计的.设男子身高X(单位:cm)服从正态分布N(175,36),求车门的最低高度.
解:设h为车门高度,X为身高
求满足条件的h,即,所以
>>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); %计算的概率密度函数值,供绘图用
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.5 随机变量的数字特征
4.5.1 平均值,中值
命令 利用mean求算术平均值
格式 mean(X) %X为向量,返回X中各元素的平均值
mean(A) %A为矩阵,返回A中各列元素的平均值构成的向量
mean(A,dim) %在给出的维数内的平均值
说明 X为向量时,算术平均值的数学含义是,即样本均值.
例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
命令 忽略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 =
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中各列元素的几何平均数构成的向量.
说明 几何平均数的数学含义是,其中:样本数据非负,主要用于对数正态分布.
例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中各列元素的调和平均值构成的向量.
说明 调和平均值的数学含义是,其中:样本数据非0,主要用于严重偏斜分布.
例4-35
>> B=[1 3 4 5]
B =
1 3 4 5
>> M=harmmean(B)
M =
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
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) %计算样本均值
则结果如下:
ans =
15.0600
命令 由分布律计算均值
利用sum函数计算
例4-40 设随机变量X的分布律为:
X
-2
-1
0
1
2
P
0.3
0.1
0.2
0.1
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
格式 D=var(X) %var(X)=,若X为向量,则返回向量的样本方差.
D=var(A) %A为矩阵,则D为A的列向量的样本方差构成的行向量.
D=var(X, 1) %返回向量(矩阵)X的简单方差(即置前因子为的方差)
D=var(X, w) %返回向量(矩阵)X的以w为权重的方差
命令 求标准差
函数 std
格式 std(X) %返回向量(矩阵)X的样本标准差(置前因子为)即:
std(X,1) %返回向量(矩阵)X的标准差(置前因子为)
std(X, 0) %与std (X)相同
std(X, flag, dim) %返回向量(矩阵)中维数为dim的标准差值,其中flag=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(默认)表示偏斜不纠正.
说明 偏斜度样本数据关于均值不对称的一个测度,如果偏斜度为负,说明均值左边的数据比均值右边的数据更散;如果偏斜度为正,说明均值右边的数据比均值左边的数据更散,因而正态分布的偏斜度为 0;偏斜度是这样定义的:
其中:μ为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 
 

admin 2007-11-29 20:48
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为二项分布的两个参数,可为标量也可为向量或矩阵.
例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
[M,V]=unifstat ( a, b)
均匀分布(连续)的期望和方差,M为期望,V为方差
unidstat
[M,V]=unidstat (n)
均匀分布(离散)的期望和方差
expstat
[M,V]=expstat (p, Lambda)
指数分布的期望和方差
normstat
[M,V]=normstat(mu,sigma)
正态分布的期望和方差
chi2stat
[M,V]=chi2stat (x, n)
卡方分布的期望和方差
tstat
[M,V]=tstat ( n)
t分布的期望和方差
fstat
[M,V]=fstat ( n1, n2)
F分布的期望和方差
gamstat
[M,V]=gamstat ( a, b)
分布的期望和方差
betastat
[M,V]=betastat ( a, b)
分布的期望和方差
lognstat
[M,V]=lognstat ( mu, sigma)
对数正态分布的期望和方差
nbinstat
[M,V]=nbinstat ( R, P)
负二项式分布的期望和方差
ncfstat
[M,V]=ncfstat ( n1, n2, delta)
非中心F分布的期望和方差
nctstat
[M,V]=nctstat ( n, delta)
非中心t分布的期望和方差
ncx2stat
[M,V]=ncx2stat ( n, delta)
非中心卡方分布的期望和方差
raylstat
[M,V]=raylstat ( b)
瑞利分布的期望和方差
Weibstat
[M,V]=weibstat ( a, b)
韦伯分布的期望和方差
Binostat
[M,V]=binostat (n,p)
二项分布的期望和方差
Geostat
[M,V]=geostat (p)
几何分布的期望和方差
hygestat
[M,V]=hygestat (M,K,N)
超几何分布的期望和方差
Poisstat
[M,V]=poisstat (Lambda)
泊松分布的期望和方差
4.5.6 协方差与相关系数
命令 协方差
函数 cov
格式 cov(X) %求向量X的协方差
cov(A) %求矩阵A的协方差矩阵,该协方差矩阵的对角线元素是A的各列的方差,即:var(A)=diag(cov(A)).
cov(X,Y) %X,Y为等长列向量,等同于cov([X Y]).
例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 Y]).
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 正整数的频率表
命令 正整数的频率表
函数 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
>> 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 %样本标准差
4.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
4.6.4 绘制正态分布概率图形
函数 normplot
格式 normplot(X) %若X为向量,则显示正态分布概率图形,若X为矩阵,则显示每一列的正态分布概率图形.
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)
图4-13
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) %data为所给样本数据,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) %data为向量,返回直方图
和正态曲线.
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) %beta为拟合系数,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) %delta为非线性最小二乘法估计的置信区间长度的一半,当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函数产生两个图: