注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

点点滴滴

感悟人生,享受精彩。

 
 
 

日志

 
 

第4章 概率统计  

2009-03-30 15:48:20|  分类: 默认分类 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
第4章  概率统计
本章介绍MATLAB在概率统计中的若干命令和使用格式,这些命令存放于MatlabR12ToolboxStats中。
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的超几何分布随机数 
Poissrnd poissrnd(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。
表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
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处的值。
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, [' ontsize{16}X~{chi}^2(4)'])
   %图中标注
text(1.5,0.05, ' ontsize{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
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
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的对角元素为相应参数的渐近方差。
命令  ?
  评论这张
 
阅读(264)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017