前言
MATLAB语言一直是国际科学界应用和影响最广泛的三大计算机数学语言之一。从某种意义上讲,在纯数学以外的领域中,MATLAB语言有着其他两种计算机数学语言Mathematica和Maple无法比拟的优势和适用性。在很多领域,MATLAB语言是科学研究者首选的计算机数学语言。目前国内外关于MATLAB语言和应用书籍数以千计,但从其覆盖面和应用水平来说,往往难以满足日益增长的MATLAB语言使用者的要求。已出版的著作从涵盖面及深度与广度上缺乏高层次、全面系统介绍高等应用数学问题各个分支的计算机求解的书籍。随着计算机的发展与普及,概率与数理统计已成为处理信息、进行决策的重要理论和方法。在科学研究中,用概率和数理统计方法从数据中获取信息和判别初步规律,往往成为重大科学发现的先导。概率与数理统计是数学方法与实际相结合,应用最为广泛、最为重要的方式之一。因此,现代科研人员和工程技术人员都应该掌握概率与数理统计的基础知识。同时概率与数理统计在自然科学、工程技术、管理科学及人文社会科学中得到越来越广泛和深入的应用,其研究的内容也随着科学技术和经济与社会的不断发展而逐步扩大。为了更好地满足高等学校培养高等技术应用型人才的需要,提高学生的基本素质和教学质量,解决高等学校概率与数理统计理论课与实践课相结合的问题,根据高等院校对数学教学的基本要求,应用数学与专业相融,基础数学为专业服务和以应用为目的,以必需、够用为度的基本原则,在多年从事高等教育教学实践的基础上,以MATLAB R2015b为平台编写了本书。本书编写时力求做到以下几点。1 数学软件命令的介绍符合学生的知识水平、浅显易理解。本书以新的MATLAB R2015b为平台,将数学理论与软件命令介绍有机地结合,使学生学会数学软件的使用方法,培养学生运用软件求解实际问题模型的能力。2 注重思想方法介绍。在阐述某一概率统计方法时,有的从具体实例开始引出相关内容的背景,有的从概念上开始,以实例总结。3 注重应用性。概率与数理统计是一门应用性很强的学科,其应用几乎遍及各个领域,成为解决实际问题的重要工具,因此,本书充实了许多应用性内容,以适应读者解决实际问题的需要。4 重视MATLAB应用于概率与数理统计方法时的简单性、实用性和可操作性,就可得到计算与分析的结果。本书主要介绍了概率与数理统计在各领域中的应用。全书共11章,各章的主要内容介绍如下。第1章对MATLAB R2015b进行概述,主要包括MATLAB的功能及发展史、MATLAB R2015b开发环境、MATLAB的语言基础等内容。第2章对概率与数理统计进行概述,主要包括概率论基础、随机变量与随机分布等内容。第3章介绍统计估计,主要包括点估计、区间估计、核密度估计和统计作图等内容。第4章介绍假设检验,主要包括正态总体参数的假设检验、其他检验等内容。第5章介绍方差分析,主要包括对方差分析进行概述、单因素方差分析和双因子方差分析等内容。第6章介绍回归分析,主要包括一元线性回归分析、多元线性回归分析、非线性回归分析和逐步回归分析等内容。第7章介绍正交实验,主要包括正交表、无交互作用的正交实验和交互作用正交实验等内容。第8章介绍主成分分析,主要包括主成分分析的概述、实现步骤及MATLAB实现等内容。第9章介绍因子分析,主要包括因子分析的概述、R型因子、Q型因子分析和目标因子分析等内容。第10章介绍判别分析,主要包括判别分析的概述、距离判别分析、Fisher判别法和Bayes判别法等内容。第11章介绍聚类分析,主要包括聚类分析的概述、距离与相似系数、K均值聚类法和模糊C均值聚类等内容。本书主要由邓奋发编写,此外参加编写的还有栾颖、周品、曾虹雁、邓俊辉、邓秀乾、邓耀隆、高泳崇、李嘉乐、李旭波、梁朗星、梁志成、刘超、刘泳、卢佳华、张棣华、张金林、钟东山、詹锦超、叶利辉、杨平和许兴杰。
本书在编写过程中,参考了大量的资料文献,在此对其作者表示感谢。本书可作为工科硕士研究生应用概率与数理统计学课程的基础教材、本科生相关专业的专业基础教材或实验教材,也可作为科研人员、工程技术人员的工具书或理论参考书。由于我们水平有限,书中难免存在不足之外,敬请读者批评指正。
作者2016年7月
第3章统计估计
统计估计是统计推断的主要内容,包括两个方面的任务。(1) 变量的分布形态未知,根据样本数据对变量的分布形态做出推测(估计)。(2) 变量的分布形态已知,即已知其概率分布函数(或概率分布律,或概率密度函数)的数学表达式,但是某些参数(或数字特征)未知,根据样本数据对未知的参数(或未知参数的函数)做出估计。3.1点估计点估计(point estimation)是用样本统计量来估计总体参数,因为样本统计量为数轴上某一点值,估计的结果也以一个点的数值表示,所以称为点估计。点估计和区间估计属于总体参数估计问题。何为总体参数统计,当在研究中从样本获得一组数据后,如何通过这组信息,对总体特征进行估计,也就是如何从局部结果推论总体的情况,称为总体参数估计。点估计也称定估计,它是以抽样得到的样本指标作为总体指标的估计量,并以样本指标的实际值直接作为总体未知参数的估计值的一种推断方法。点估计的方法有矩估计法、顺序统计量法、最大似然法、最小二乘法等。3.1.1矩估计在统计学中,矩是指以期望为基础而定义的数字特征,一般分为原点矩和中心矩。设X为随机变量,对任意正整数k,称EXk为随机变量k阶原点矩,记为
mk=EXk
当k=1时,m1=EX=。可见一阶原点矩为随机变量X的数学期望。我们把ck=E[X-EXk]称为以EX为中心的k阶中心矩。显然,当k=2时,c2=E[X-EX2]=2。可见,二阶中心矩为随机变量X的方差。【例31】已知某种灯泡的寿命X~N,2,其中,、2都是未知的,今随机取得4只灯泡,测得寿命(单位: 小时)为1052、1453、1367、1650,试估计和。解: 因为为全体灯泡的平均寿命,x-为样本的平均寿命,很自然地会想到用x-去估计; 同理用S去估计。由于
x-=141502 1453 1367 1650=1493
S2=1502-14932 1453-14932 1367-14932 1650-149324-1=14068.7
S=118.61
因此,和的估计值分别为1493小时及118.61小时。矩估计简便、直观、比较常用,但是矩估计法也有其局限性。首先,它要求总体的k阶原点矩存在,如果不存在则无法估计; 其次,矩估计法不能充分地利用估计时已掌握的有关总体分布形式的信息。通常设为总体X的待估计参数,一般用样本X1,X2,,Xn构成一个统计量^=^X1,X2,,Xn来估计,则称^为的估计量。对于样本的一组数值x1,x2,,xn,估计量^的值^x1,x2,,xn称为的估计值。于是点估计即是寻求一个作为待估计参数的估计量^x1,x2,,xn的问题。但是必须注意,对于样本的不同数值,估计值是不相同的。如在上例中,分别用样本平均数和样本修正方差来估计总数数学期望和总体均方差,即有
^=^X1,X2,,Xn=1nni=1Xi=
^=^X1,X2,,Xn=ni=1Xi-2n-1=S
其中对应于给定的估计值=x-=1493h,^=S=118.61h。3.1.2极大似然估计极大似然估计方法(Maximum Likelihood Estimate,MLE)也称为最大概率似然估计或最大似然估计,是求估计的另一种方法,1821年首先由德国数学家C.F.Gauss(高斯)提出,但是这个方法通常被归功于英国的统计学家R.A.Fisher(罗纳德费希尔),他在1922年的论文On the mathematical foundations of theoretical statistics, reprinted in Contributions to Mathematical Statisticsby R.A.Fisher,1950, J. Wiley & Sons, New York中再次提出了这个思想,并且首先探讨了这种方法的一些性质,极大似然估计这一名称也是费希尔给的。这是一种目前仍然得到广泛应用的方法。下面分别讨论X为离散型和连续型随机变量时总体中某些参数的最大似然估计。设总体X是离散型随机变量,其分布律P{X=x}=p{x;},其中是未知参数,如果取得样本观测值为x1,x2,,xn,则表示随机事件X1=x1,X2=x2,,Xn=xn发生了。考虑n个事件X1=x1,X2=x2,,Xn=xn的交点的概率,注意到X1,X2,,Xn的独立性,即有
L=P{X1=x1,X2=x2,,xn=xn}
=P{X1=x1}P{X2=x2}P{Xn=xn}
=p{x1;}p{x2;}p{xn;}
=ni=1p{x;}(31)
函数L称为似然函数,对于已给定的x1,x2,,xn,它是未知参数的函数。按极大似然估计法的直观想法是: 若抽样的结果得到样本观测值x1,x2,,xn,则应当这样选取参数L的值,使这组样本观测值出现的可能性最大,也就是使似然函数L达到最大值,从而求得参数的估计值^,利用极大似然估计法求得的参数估计值称为极大似然估计值。极大似然估计值的问题,就是求似然函数L的最大值问题,这个问题可以通过解下面的方程
dLd=0(32)
来解决。因为lnL是L的增函数,所以lnL与L在的同一值处取得最大值。因此,也可将方程(32)换成下面的方程
dlnLd=0(33)
解方程(32)或式(33)得到的^就是参数的最大似然估计值,而从后一方程求解往往比较方便,式(33)称为对数似然方程。【例32】设有甲、乙两个布袋,甲袋中有99个白球和1个黑球,乙袋中有1个白球和99个黑球。由于某种原因已不能识别哪一个是甲袋,哪一个是乙袋。你能否用统计的方法识别出来?解: 下面对这一问题进行数学描述与分析。不妨设变量X表示袋中的白球数,则X~199p1-p,p是未知的分布参数,其取值依赖于变量X代表的是甲袋中的白球数还是乙袋中的白球数。显然,变量X代表的是甲袋中的白球数与p=99100是等价的,变量X的代表是乙袋中的白球数与p=1100是等价的。可以通过抽样(任取一袋,从该袋中任取一球,观察其颜色)的方法来确定p=99100还是p=1100。设事件A表示取出的一袋为甲袋,事件B表示从袋子中取出的是白球,则
PA=0.5,PB|A=99100,PB|=1100
假定取出的是白球。在已知取出的是白球的条件下,判断该球来自甲袋还是乙袋的问题,可由贝叶斯公式,通过比较概率PB|A和P|B的大小来做出判断。由于在一次试验中大概率事件容易发生,因此,若PA|BP|B,则该球来自甲袋; 如果PA|BP|B,p=1100等价于PA|BP|B。通过计算可知,PA|BP|B,因此p=99100,即现在取出的这一袋是甲袋。概括这里的思想方法,就可以得到极大似然估计法的数学原理大概率原理: 大概率事件在一次试验中容易发生。或者说,在一次试验中已经发生的事件具有较大的概率,而变量的分布参数有助于关于该变量的大概率事件的发生。【例33】设X~N,2,求和2的极大似然估计。解: 正态样本N,2的密度函数是12e-x-222,则似然函数为
L,2=ni=112e-xi-222=122n2eni=1xi-222
将其取对数,并令关于、2的一阶导数为零,则得
lnL,2=12ni=1xi-=0
lnL,22=-n22 1222ni=1xi-2=0
解此关于,2的方程组,得驻点
=x-=1nni=1xi,2=1nni=1xi-2
又可求得对数似然函数的二阶导函数矩阵是非正定矩阵,因此驻点处即为似然函数的极大值点处,并将的样本表达式代入2的驻点表达式,得与2的极大似然估计为
^=x-=1nni=1xi,^2=1nni=1xi-x-2
表31所示函数的返回值为数据向量x的参数最大似然估计值,以及置信度为(1-a)100%的置信区间。a的默认值为0.05,即置信度为95%。
表31参数估计函数
函数名调 用 格 式函 数 说 明
binofitphat = binofitx,n二项分布的概率最大似然估计[phat,pci] = binofitx,n置信度为95%的参数估计和置信区间[phat,pci] = binofitx,n,alpha返回水平a的参数估计和置信区间
续表
函数名调 用 格 式函 数 说 明
poissfitlambdshat = poissfitdata泊松分布的参数的最大似然估计[lambdahat,lambdaci] = poissfitdata置信度为95%的参数估计和置信区间[lambdahat,lambdaci] = poissfitdata,alpha返回水平a的参数估计和置信区间normfit[muhat,sigmahat] = normfitdata正态分布的最大似然估计, 置信度为95%[muhat,sigmahat,muci,sigmaci] = normfitdata,alpha返回水平a的期望、方差值和置信区间betafitphat = betafitdata返回分布参数a和b的最大似然估计[phat,pci] = betafitdata,alpha返回最大似然估计值和水平a的置信区间unifit[ahat,bhat] = unifitdata均匀分布参数的最大似然估计[ahat,bhat,ACI,BCI] = unifitdata置信度为95%的参数估计和置信区间[ahat,bhat,ACI,BCI] = unifitdata,alpha返回水平a的参数估计和置信区间expfitmuhat = expfitdata指数分布参数的最大似然估计[muhat,muci] = expfitdata置信度为95%的参数估计和置信区间[muhat,muci] = expfitdata,alpha返回水平a的参数估计和置信区间gamfitphat = gamfitdatar分布参数的最大似然估计[phat,pci] = gamfitdata置信度为95%的参数估计和置信区间[phat,pci] = gamfitdata,alpha返回最大似然估计值和水平a的置信区间wblfitparmhat = wblfitdata韦伯分布参数的最大似然估计[parmhat,parmci] = wblfitdata置信度为95%的参数估计和置信区间[parmhat,parmci] = wblfitdata,alpha返回水平a的参数估计及其区间估计mle[phat=mle''dist'', data分布函数名为dist的最大似然估计[phat, pci]=mle''dist'', data置信度为95%的参数估计和置信区间[phat, pci]=mle''dist'', data, alpha返回水平a的最大似然估计值和置信区间[phat, pci]=mle''dist'', data, alpha, p1仅用于二项分布,p1为试验总次数【例34】随机产生100个服从正态分布N2,0.52的样本数据X,并用这些数据估计总体N,2中的参数、,求出参数的最大似然估计值和置信水平为99%的置信区间。分析: 随机产生的100个数据可视为总体中抽出容量为100的样本,样本的观测值就是这具体的100个数据,可用命令normfitX, alpha求出参数、的估计。其MATLAB代码编程如下。
clear all;
X=normrnd2,0.5,100,1;%产生100个样本数据
[muhat,sigmahat,muci,sigmaci]=normfitX,0.01运行程序,输出如下。
muhat =
2.0240
sigmahat =
0.4343
muci =
1.9099
2.1380
sigmaci =
0.3665
0.5298说明: 参数、的估计最大似然值分别为2.0240、0.4343,参数,的置信水平为99%的置信区间分别为[1.9099, 2.1380]、[0.3665, 0.5298]。这一估计结果和总体N,2中的参数真实数值=2,=0.5是非常接近的。可以概括出求极大似然估计值的一般步骤如下。(1) 明确变量的分布律和密度函数;(2) 写出似然函数L;(3) 求似然函数L的最大值点,得^MLE;(4) 应用问题中,将样本数据代入^MLE,求出具体的估计值。值得注意的是,求解对数似然方程组是在假定其可导并且导数变号的基础上,如果不满足这一条件,需针对似然函数L1,2,,k的单调性,利用极大似然估计的基本原理直接进行L1,2,,k的最大值问题的讨论。极大似然估计量有一个简单而有用的性质: 设的函数g=g是上的实值函数,且有唯一反函数。如果^是的极大似然估计量,则g^也是g的极大似然估计量。这个性质称为极大似然估计的不变性。根据这一性质可以使一些复杂结构的参数的极大似然估计问题简单化。极大似然估计法是在变量分布类型已知的情况下使用的一种参数估计法。一般地,用极大似然法所得的估计的性质比用矩估计法所得的要好,故通常多用极大似然法。在MATLAB中,提供了mle函数进行极大似然估计,函数的调用格式如下。
phat = mledata
[phat,pci] = mledata
[] = mledata,''distribution'',dist
[] = mledata,,name1,val1,name2,val2,
[] = mledata,''pdf'',pdf,''cdf'',cdf,''start'',start,
[] = mledata,''logpdf'',logpdf,''logsf'',logsf,''start'',start,
[] = mledata,''nloglf'',nloglf,''start'',start,
[phat, pci]=mledata, ''distribution'', dist, ''alpha'', a, ''ntrials'', n其中,输出参数phat是指定分布的参数的极大似然估计值(多参数时为行向量),pci是参数的区间估计的置信上限和下限(与参数对应的二维列向量,可以缺省)。输入参数data是样本数据向量(不可缺省)。引用参数''distribution''及其取值dist设置变量的分布类型(应用中dist要用具体的分布名称字符串替换并用单引号引起),二者要成对出现(可以同时缺省,缺省时分布类型默认为正态分布)。引用参数''alpha''及其取值a设置区间估计的显著性水平,二者成对出现(可以同时缺省,缺省时默认为0.05,即置信水平为0.95)。引用参数''ntrials''及其取值n仅在分布类型为二项分布时引用(对于其他分布可以缺省),设置二项分布中试验的次数。dist的取值包括: Beta,Bernoulli,Binomial,Discrete Uniform,Exponential,Extreme Value,Gamma,Geometric,Lognormal,Negative Binomial,Normal,Poisson,Rayleigh,Uniform,Weibull。【例35】引用常数的测定值服从均值为、标准差为的正态分布。某人在实验中使用金球测定引力常数,6次测定观察值为: 6.683,6.681,6.676,6.678,6.679,6.672。试用极大似然估计法对未知参数和做出估计。其MATLAB代码编程如下:
clear all;
x=[6.683,6.681,6.676,6.678,6.679,6.672];
phat=mlex,''distribution'',''norm'',''alpha'',0.05运行程序,输出如下:
phat =
6.67820.0035
即金球测定的估计值为6.6782,的估计值为0.0035。其实,此例计算中mle函数的调用可以简化为p=mlex。3.1.3顺序统计量1. 统计量法的定义
设 1, 2,, n是总体 的样本,将其按大小排列为 *1 *2 *n,则称 *1, *2,, *n为顺序统计量。明显地, *1与 *n分别为样本的最小值与最大值。称 -= *k 1,n=2k 1
*k*k 12,n=2k为样本中位数。样本中位数的取值规则为: 将样本值x1,x2,,xn从小至大排成x*1x*2x*n,当n=2k 1时, -取居中的数据x*k 1为其观测值; 当n=2k时, -取居中的两个数据的平均值 *k*k 12为其观测值,中位数 -带来了总体 取值的平均数的信息,因此用 -估计总体 的数学期望是合适的。用样本中位数 -估计总体 的数学期望的方法称为数学期望的顺序统计量估计法。顺序统计量估计法的优点是计算简便,且不易受个别异常数据的影响。如果一组样本值某一数据异常(如过小或过大),则这个异常数据可能是总体 的随机性造成的,也可能是受外来干扰造成的(如工作人员粗心,记录错误),当原因属于后者,用样本平均值x-估计Ex显然受到影响,但用样本中位数 -估计Ex时,由于一个(甚至几个)异常的数据不易改变中位数取值,所以估计值不易受到影响。即称R= *n- *1为样本极差。由于样本极差带来总体样本取值离散程度的信息,因此可以用R作为对总体 的标准差的估计(R与量纲相同)。用样本极差对总体 的标准差做估计的方法称为极差估计法。极差估计法的优点是计算简便,但不如用S可靠,n越大两者可靠的程度差别越大,这时一般不用极差估计。2. 顺序统计量法主要适用范围顺序统计量法主要适用于正态总体,当总体不是正态分布,但是连续型且分布密度对称时,也常用样本中位数来估计总体的期望。3.1.4最小二乘法在科学实验数据处理中,往往要根据一组给定的实验数据xi,yii=0,1,,m,求出自变量x与因变量y的函数关系y=sx,a0,,ann<m,这是ai为待定参数,由于观测数据总有误差,且待定参数ai的数量比给定数据点的数量少(即n<m),因此它不同于插值问题。这类问题不要求y=sx=sx,a0,,an通过点xi,yii=0,1,,m,而只要求在给定点xi上的误差i=sxi-yii=0,1,,m的平方和mi=02i最小。当Sxspan{0,1,,n}时,即
sx=a00x a11x annx(34)
这里0x,1x,,nxC[a,b]是线性无关的函数族,假定在[a,b]上给出一组数据{xi,yi,i=0,1,,m},axib以及对应的一组权{i}m0,这里i0为权系数,要求sx=span{0,1,,n}使Ia0,a1,,an最小,其中
Ia0,a1,,an=mi=0i[sxi-yi]2(35)
这就是最小二乘逼近,得到的拟合曲线为y=sx,这种方法称为曲线拟合的最小二乘法。式(35)中Ia0,a1,,an实际上是关于a0,a1,,an的多元函数,求I的最小值就是求多元函数I的极值,由极值必要条件,可得
Iak=2mi=0i[a00xi a11xi annxi-yi]kxi,k=0,1,,n(36)
根据内积定义引入相应带权内积记号
j,k=mi=0ijxikxi
y,k=mi=0iyikxi(37)
则(36)可改写为
0,ka0 1,ka1 n,kan=y,k,k=0,1,,n
这是关于参数a0,a1,,an的线性方程组,用矩阵表示为
0,00,10,n
1,01,11,n
n,0n,1n,na0a1an=y,0y,1y,n(38)
式(38)称为法方程。当{jx;j=0,1,,n}线性无关,且在点集X={x0,x1,,xm}mn上至多只有n个不同零点,则称0,1,,n在X上满足Haar条件,此时式(38)的解存在唯一。记式(38)的解为
ak=a*k,k=0,1,,n
从而得到最小二乘拟合曲线
y=s*x=a*00x a*11x a*nnx(39)
可以证明对a0,a1,,anTRn 1,有
Ia*0,a*1,,a*nIa0,a1,,an
故式(39)得到的s*x即为所求的最小二乘解。它的平方误差为
‖‖22=mi=0i[s*xi-yi]2(310)
均方误差为
‖‖2=mi=0i[s*xi-yi]2
在最小二乘逼近中,若取kx=xkk=0,1,,n,则sxspan{1,x,,xn},表示为
sx=a0 a1x anxn(311)
此时关于系数a0,a1,,an的方程(38)是病态方程,通常当n3时都不直接取kx=xk作为基。3.1.5点估计的优良性准则样本统计量,如样本均值,样本标准差S,样本成数如何用于对相应总体参数、和p的点估计值。直观上,这些样本统计量对相应总体参数的点估计值是很有吸引力的。然而,在用一个样本统计量作为点估计量之前,统计学应检验说明这些样本统计量是否具有某些与好的点估计量相联系的性质。本节讨论点估计量的性质: 无偏性、有效性和一致性。1. 无偏性设^=^X1,X2,,Xn的数学期望等于,即
E^=
则称^是参数的无偏估计量; 如果样本观测值为x1,x2,,xn,则称^x1,x2,,xn为参数的无偏估计值。在科学技术中E^-称为以^作为的估计的系统误差,无偏估计的实际意义就是无系统误差。无偏性是对估计量的一个最重要、最常见的要求,它的实际意义在于,当这个估计量经常使用时,在多次重复的平均意义下,给出了接近于真值的估计,在此应当指出,同一个参数的无偏估计量不是唯一的。例如,有EXi=,这表明任一样本的每一分量Xii=1,2,,n都是总体均值的无偏估计量,在参数的许多无偏估计中,当然是以对的平均偏差较小者为好,即较好的估计量应当有尽可能小的方差。因此,便有了第二个评选标准。下面列举出关于无偏性的几个重要结论。1 无论变量X服从何种分布,样本的k阶原点矩Ak=1nni=1Xkii=1,2,,n是变量X的k阶原点矩EXk的无偏估计。自然,是EX的无偏估计。2 无论变量X服从何种分布,样本(修正)方差S2=1n-1ni=1Xi-2是变量X的方差2的无偏估计。3 样本方差(二阶中心矩)B2不是变量方差2的无偏估计,但是limnEB2=2,所以B2是2的渐近无偏估计。4 样本标准差S=1n-1ni=1Xi-2不是变量X的标准差的无偏估计。但是,在变量的正态性假设下,可将样本标准差修正为^S=CnS,^S是的无偏估计,其中Cn=n-12n-12n2称为正态标准差的无偏系数。由于limnCn=1,所以S是的渐近无偏估计。无偏性准则是对估计量的一个基本要求。无偏性估计的统计意义是指估计量不产生系统性的偏差。例如,用样本均值作为变量均值的估计时,由于是随机变量,因此在一次估计中的实现值与其真值之间存在偏差-。这种偏差是随机的,虽无法说明一次估计所产生的偏差,但是对同一统计问题大量重复使用估计时,实际产生的偏差-随机地在0的周围波动,不会产生系统的偏大(小)于的情况。渐近无偏差是指估计量存在系统性的偏差,但是这种系统性偏差随着样本容量的增加而趋向于消失。【例36】设总体X~X2n,X1,X2,,X20为来自总体的简单随机样本,想要估计总体均值(注意n未知),比较以下三个点估计量的好坏: ^1=101X1-100X2,^2=12X10 X11,^3=。实例中给出了利用MSE评价点估计量的随机模拟方法。由于X2n的总体均值为n,因此可以先取定一个固定值,例如n=0=5,然后在这个参数已知且固定的总体中抽取容量为20的样本,分别用样本值依照三种方法分别计算估计值,看看哪种方法误差大,哪种方法误差小。一次估计的比较一般不能说明问题,正如低手射击也可能命中10环,高手射击也可能命中9环,如果连续射击10000次,比较总环数,多者一定是高手。同理,如果抽取容量为20的样本N=10000次,分别计算: MSE^i1Nnk=1[^ik-0]2,值小者为好。其MATLAB代码编程如下:
clear all;
N=10000;
m=5;n=20;
mse1=0;mse2=0;mse3=0;
for k=1:N
x=chi2rndm,1,n;
m1=101*x1-100*x2;
m2=medianx;
m3=meanx;
mse1=mse1 m1-m^2;
mse2=mse2 m2-m^2;
mse3=mse3 m3-m^2;
end
mse1=mse1N
mse2=mse2N
mse3=mse2N运行程序,输出如下:
mse1 =
1.9716e 05
mse2 =
0.9909
mse3 =
9.9087e-052. 有效性一个未知参数的估计量^仅有无偏性是不够的。因为一方面,无偏性仅反映估计量在参数真值周围波动,而没有反映出集中的程度; 另一方面,一个参数的无偏估计量可能不止一个,对于数学期望,样本均值是它的无偏估计量,样本的第一个观测值X1也是它的无偏估计量(因EX1=),那么哪个更好呢?仅有无偏性一个标准是不能确定的。一个自然的想法是进一步比较它们的方差,方差越小,表示^越集中在的附近,从这个意义上讲方差越小的无偏估计量越好。设^1,^2都是的无偏估计量,若D^10,有
limnP{|^-|}=0
而且这对的一切可能取的值都成立,则称^是参数的一个一致性估计。一致性准则是对一个估计量最基本的要求。它说明,随着样本容量的增大,一个好的估计量^应该越来越靠近参数的真值,使绝对偏差|^-|较大的概率越来越小。如果一个估计量没有一致性,那么,不论样本取多大,我们也不可能把未知参数估计到预定的精度。这种估计量显然是不可取的。下面给出一致性估计的几个重要结论。1 一致性估计具有不变性。即当^1,^2,,^k分别是1,2,,k的一致性估计时,如果g1,2,,k为连续函数,则g^1,^2,,^k是g1,2,,k的一致性估计。2 样本的k阶原点矩Ak=1nni=1Xki是变量X的k阶原点矩EXk的一致性估计,因此样本均值是变量均值的一致性估计。3 样本的二阶中心矩B2=1nni=1Xi-2是变量X的方差2的一致性估计。4 样本方差S2=1n-1ni=1Xi-2是变量的方差2的一致性估计,样本标准差S=1n-1ni=1Xi-2是变量的标准差的一致性估计。5 事件发生的频率是其概率的一致性估计。6 极大似然估计量往往具有一致性。3.2区间估计上一节讨论了参数点估计,它是用样本算得的一个值去估计未知参数,但是,点估计值仅仅是未知参数的一个近似值,它没有反映出这个近似值的误差范围,使用起来把握不大,区间估计(interval estimation)正好弥补了点估计的这个缺陷。3.2.1区间估计简介区间估计就是以一定的概率保证估计包含总体参数的一个值域,即根据样本指标和抽样平均误差推断总体指标的可能范围。它包括两部分内容: 一是这一可能范围的大小; 二是总体指标落在这个可能范围内的概率。区间估计既说清估计结果的准确程度,又同时表明这个估计结果的可靠程度,所以区间估计是比较科学的。用样本指标来估计总体指标,要达到100%的准确而没有任何误差,几乎是不可能的,所以在估计总体指标时就必须同时考虑估计误差的大小。从人们的主观愿望上看,总是希望花较少的钱取得较好的效果,也就是说希望调查费用和调查误差越小越好。但是,在其他条件不变的情况下,缩小抽样误差就意味着增加调查费用,它们是一对矛盾。因此,在进行抽样调查时,应该根据研究目的和任务以及研究对象的标志变异程度,科学确定允许的误差范围。区间估计必须同时具备3个要素。即具备估计值、抽样极限误差和概率保证程度3个基本要素。抽样误差范围决定抽样估计的准确性,概率保证程度决定抽样估计的可靠性,二者密切联系,但同时又是一对矛盾,所以,对估计的精确度和可靠性的要求应慎重考虑。3.2.2区间估计的含义区间估计就是根据样本来确定统计量X1,X2,,Xn和-X1,X2,,Xn,使
PX1,X2,,Xn clear all;
x=[1.1 2.2 3.3 4.4 5.5];
n=lengthx;
m=meanx;
c=2.3sqrtn;
d=c*norminv0.975;
a1=m-d;
b1=m d;
[a1,b1]运行程序,输出如下:
ans =
1.28405.31602) 2为未知时,均值的置信区间这时,自然会想到以样本标准差S代替总体均方差,即知选取统计量
T=-Sn~tn-1
对给定的数,由
P|T| clear all;
x=[1.1 2.2 3.3 4.4 5.5];
n=lengthx;
m=meanx;
S=stdx;
dd=S*tinv0.975,4sqrtn;
a2=m-dd;
b2=m dd;
[a2,b2]运行程序,输出如下:
ans =
1.14045.45963) 单正态方差的区间估计设总体X~N,2,X1,X2,,Xn是X的样本,求2的1-置信区间。由2分布知选取统计量
n-1S22~2n-1
对给定的,取2分布分位点22n和21-2n,使
21-2n-1 clear all;
x=[1.1 2.2 3.3 4.4 5.5];
n=lengthx;
c1=chi2inv0.025,4;
c2=chi2inv0.975,4;
T=n-1*varx;
a3=Tc2;
b3=Tc1;
[a3,b3]运行程序,输出如下:
ans =
1.085924.97844) 两正态总体均值差的置信区间当方差已知时,设X1,X2,,Xm~N1,21,Y1,Y2,,Ym~N2,22,两样本独立,此时1-2的置信区间为
--u221m 22n,- u221m 22n
在此已经知道u2可用norminv0.975求得。当方差未知但相等时,此时1-2的置信区间为
--t2C,- t2C
其中,C=1m 1nm-1S21 n-1S22m n-2,而t2依照自由度m n-2计算。5) 两正态总体方差比的置信区间查自由度为m-1,n-1的F分布临界值表使得,
Pc1 clear all;
x=[-0.12 -0.80 -0.05 -0.04 -0.01 0.05 0.07 0.21];
y=[-1.50 -0.80 -0.40 -0.10 0.20 0.61 0.82 1.24];
v1=varx;
v2=vary;
c1=finv0.025,7,7;
c2=finv0.975,7,7;
a4=v1v2c1;
b4=v1v2c1;
[a4,b4]运行程序,输出如下:
ans =
0.57200.5720方差比小于1的概率至少达到95%,说明车床A的精度明显高。2. 单侧置信区间上面的讨论中,对于未知参数,给出两个统计量和-,得到的置信区间为,-的形式。但在有些实际应用中,常常只关心参数的上限或下限。例如,对于设备、元件的寿命来说,我们只关心平均寿命至少是多少(的下限); 与之相反,在考虑化学药品中杂质含量时,我们关心的却是平均杂质含量最多是多少(的上限)。这就引出了单侧置信区间的概念。对于给定值0=1-
则称随机区间X1,X2,,Xn, 是的置信水平为1-的下侧置信区间,称X1,X2,,Xn是置信水平为1-的置信下限。又如果统计量-X1,X2,,Xn满足对任意有
P-tan-1Sn=1-
可得的置信度为1-的单侧置信下限为
-tn-1Sn
由所得数据计算,有
x-=1160,s=99.75,n=5,=0.05
查表得t0.054=2.14,从而平均寿命的置信度为95%的置信下限为
x--tn-1sn=1064.56
也就是说,该批灯泡的平均寿命至少在1064.56h以上,可靠程度为95%。其MATLAB代码编程如下:
clear all;
x=[1050,1100,1120,1250,1280];
N=lengthx;
muEST=meanx
muLOWER=muEST-tinv0.95,N-1*sqrtvarxN上述指令的运行结果是:
muEST =
1160
muLOWER =
1.0649e 003计算结果表明,这批灯泡的平均寿命约为1160h,以95%的概率保证这批灯泡的平均寿命不低于1065h。3.2.5区间估计函数在MATLAB的函数工具箱中,也提供了相关函数用于实现区间估计,下面分别对这些函数给予介绍。1. nlinfit函数在MALTAB中,提供了nlinfit函数用于求解高斯牛顿法的非线性最小二乘数据拟合。函数的调用格式如下。beta = nlinfitX,y,fun,beta0: 返回在fun中描述的非线性函数的系数。fun为用户提供的形如y^=f,x的函数,该函数返回已给初始参数估计值和自变量x的y的预测值y^。[beta,r,J,COVB,mse] = nlinfitX,y,fun,beta0: 同时返回的beta为拟合系数,r为残差,J为jacobi矩阵,COVB为评估的协方差矩阵,mse为误差的方差。输入参数beta0为初始预测值。[] = nlinfitX,y,fun,beta0,options: 指定控制参数后返回值。参数options包括MaxIter、TolFun、TolX、Display、DerivStep等。当X为矩阵时,则X的每一列为自变量的取值,y是一个相应的列向量。如果fun中使用了@,则表示函数的句柄。【例314】使用nlinfit函数求高斯牛顿法的非线性最小二乘数据拟合。其MATLAB代码编程如下:
clear all;
S = load''reaction'';
X = S.reactants;
y = S.rate;
beta0 = S.beta;
%利用nlinfit函数求非线性最小二乘数据拟合
beta = nlinfitX,y,@hougen,beta0运行程序,输出如下:
beta =
1.2526
0.0628
0.0400
0.1124
1.19142. nlparci函数在MATLAB中,提供了nlparci函数用于求解非线性模型的参数估计的置信区间。函数的调用格式为ci = nlparcibeta,resid,''covar'',sigma:返回置信度为95%的置信区间,beta为非线性最小二乘法估计的参数值,resid为残差,sigma为协方差矩阵系数。ci = nlparcibeta,resid,''jacobian'',J:返回置信度为95%的置信区间,beta为非线性最小二乘法估计的参数值,resid为残差,J为Jacobian矩阵。ci = nlparci,''alpha'',alpha:返回1001-alpha%的置信区间。
【例315】利用nlparci函数求求非线性模型yj=1 2exp-3xj j的参数估计的置信区间。
clear all;
%用函数句柄表示模型
mdl = @a,xa1 a2*exp-a3*x;
rng9845,''twister'' %可重复性
a = [1;3;2];
x = exprnd2,100,1; %指数分布
epsn = normrnd0,0.1,100,1; %正态分布
y = mdla,x epsn;
%数据拟合模型为随机的
a0 = [2;2;2];
[ahat,r,J,cov,mse] = nlinfitx,y,mdl,a0;
ahat
%检查在95%的置信区间是否[1 3 2]是使用雅可比参数nlparci
ci1 = nlparciahat,r,''Jacobian'',J
%使用的协方差参数
ci2 = nlparciahat,r,''covar'',cov
运行程序,输出如下:
ahat =
1.0153
3.0229
2.1070
ci1 =
0.98691.0438
2.94013.1058
1.99632.2177
ci2 =
0.98691.0438
2.94013.1058
1.99632.2177
3. nlintool函数在MATLAB中,提供了nlintool函数用于求解非线性拟合并显示交互图形。函数的调用格式如下。nlintoolX,y,fun,beta0: 返回数据X,y的非线性曲线的预测图形,它用两条红色曲线预测全局置信区间。beta0为参数的初始预测值,默认值为0.05,即置信度为95%。nlintoolX,y,fun,beta0,alpha: 将置信度设置为1-alpha100%。nlintoolX,y,fun,beta0,alpha,''xname'',''yname'': 给X和y的变量分别赋予变量名xname和yname。【例316】使用nlintool函数求非线性拟合并显示交互图形。其MATLAB代码编程如下:
clear all;
load reaction
nlintoolreactants,rate,@hougen,beta,0.01,xn,yn运行程序,效果如图31所示。
图31非线性拟合交互图形
4. nlpredci函数在MATLAB中,提供了nlpredci函数用于求解非线性模型置信区间预测。函数的调用格式如下。[ypred,delta] = nlpredcimodelfun,x,beta,resid,''jacobian'',J: 返回预测值ypred,fun与前面相同,beta为给出的适当参数,resid为残差,J为Jacobi矩阵,x为非线性函数中的独立变量的矩阵值。返回的delta为非线性最小二乘法估计的置信区间长度的一半。当resid长度超过beta的长度,并且J的列满秩时,置信区间的计算才是有效的,[ypred-delta,ypred delta]为置信度,是95%的不同步置信区间。[ypred,delta] = nlpredcimodelfun,x,beta,resid,''covar'',sigma: 参数sigma为协方差矩阵系数。[] = nlpredci,param1,val1,param2,val2,: 设置多个参数的名称及其对应的值,参数包括: ''alpha''、''mse''、''predopt''和''simopt''。【例317】使用nlpredci函数求非线性最小二乘预测置信区间。其MATLAB代码编程如下:
clear all;
S = load''reaction'';
X = S.reactants;
y = S.rate;
beta0 = S.beta;
[beta,R,J] = nlinfitX,y,@hougen,beta0;
[ypred,delta] = nlpredci@hougen,meanX,beta,R,''Jacobian'',J运行程序,输出如下:
ypred =
5.4622
delta =
0.19213.3参数估计实例下面通过一个实例来演示怎样利用参数估计实现工程实际领域中的应用。【例318】分别使用金球和铂球测定引力常数(单位: 10-11m3kg-1s-2)。(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的置信区间,并求的置信水平为0.9的置信区间。其MATLAB代码编程如下:
clear all;
data1=[6.683 6.681 6.676 6.678 6.679 6.672];
alpha=0.1;
[muhat1,sigmahat1,muci1,sigmaci1]=normfitdata1,alpha
[phat1,pci1]=mledata1,''distribution'',''normal'',''alpha'',alpha
data2=[6.661 6.661 6.667 6.667 6.664];
[muhat2,sigmahat2,muci2,sigmaci2]=normfitdata2,alpha
[phat2,pci2]=mledata1,''distribution'',''normal'',''alpha'',alpha运行程序,输出如下:
muhat1 =
6.6782
sigmahat1 =
0.0039
muci1 =
6.6750
6.6813
sigmaci1 =
0.0026
0.0081
phat1 =
6.67820.0035
pci1 =
6.67500.0026
6.68130.0081
muhat2 =
6.6640
sigmahat2 =
0.0030
muci2 =
6.6611
6.6669
sigmaci2 =
0.0019
0.0071
phat2 =
6.67820.0035
pci2 =
6.67500.0026
6.68130.00813.4核密度估计核密度估计(Kernel Density Estimation)是在概率论中用来估计未知的密度函数,属于非参数检验方法之一,由Rosenblatt 1955和Emanuel Parzen1962提出,又名Parzen窗(Parzen Window)。Ruppert和Cline基于数据集密度函数聚类算法提出修订的核密度估计方法。3.4.1核密度估计的概述对于一组关于X和Y观测数据{xi,yi}ni=1,假设它们存在关系yi=mxi i,通常我们的目的在于估计mx的形式。在样本数量有限的情况下,我们无法准确估计mx的形式。这时,可以采用非参数方法,在非参数方法中,并不假定也不固定mx的形式,仅假设mx满足一定的光滑性,函数在每一点的值都由数据决定。显然,由于随机扰动的影响数据有很大的波动,极不光滑,因此要去除干扰使图形光滑。最简单最直接的方法就是取多点平均,也就是每一点mx的值都由离x最近的多个数据点所对应的y值的平均值得到。显然,如果用来平均的点越多,所得的曲线越光滑。当然,如果用n个数据点来平均,则mx为常数,这时它最光滑,但失去了大量的信息,拟合的残差也很大。所以说,这就存在了一个平衡的问题,也就是说,要决定每个数据点在估计mx的值时要起到的作用问题。直观上,和x点越近的数据对决定mx的值所起作用越大,这就需要加权平均。因此,如何选择权函数来光滑及光滑到何种程度即是我们这里所关心的核心问题。3.4.2核密度估计的形式对于数据x1,x2,,xn,核密度估计的形式如下。
fhx=1nhni=1kx-xih
这是一个加权平均,而核函数(kernal function)K为一个权函数,核函数的形状和值域控制着用来估计fx在点x的值时所用数据点的个数和利用的程度,直观来看,核密度估计的好坏依赖于核函数和带宽h的选取。通常考虑的核函数关于原点对称且其积分为1,下面四个函数为最为常用的权函数:Uniform:
12I|t|1Epanechikov:
34I1-t2I|t| clear all;
rng default
x = [randn30,1; 5 randn30,1];
[f,xi] = ksdensityx;
figure
plotxi,f;运行程序,效果如图32所示。
图32核密度估计曲线图
【例320】核密度估计的案例分析。其MATLAB代码编程如下:
%对MATLAB自带的数据绘制其直方图
clear all;
cars = load''carsmall'',''MPG'',''Origin'';
MPG = cars.MPG;
histMPG%效果如图33所示
xlabel''样本'';ylabel''直方图'';
setgetgca,''Children'',''FaceColor'',[.8 .8 1]
%绘制不同的窗宽核密度估计曲线,效果如图34所示
[f,x,u] = ksdensityMPG;
plotx,f
title''MPG的核密度估计''
hold on
[f,x] = ksdensityMPG,''width'',u3;
plotx,f,'':r'';
[f,x] = ksdensityMPG,''width'',u*3;
plotx,f,''--g'';
legend''默认窗宽'',''默认13窗宽'',''默认3倍窗宽'';
xlabel''x'';ylabel''f'';
hold off
%设置不同核密度估计参数,绘制相应曲线,效果如图35所示
hname = {''normal'' ''epanechnikov'' ''box'' ''triangle''};
colors = {''r'' ''b'' ''g'' ''m''};
for j=1:4
[f,x] = ksdensityMPG,''kernel'',hname{j};
plotx,f,colors{j};
hold on;
end
legendhname{:};
xlabel''x'';ylabel''f'';
hold off
%比较密度估计,显示了从不同的原产地国家的汽车省油分布曲线,效果如图36所示
Origin = cellstrcars.Origin;
I = strcmp''USA'',Origin;
J = strcmp''Japan'',Origin;
K = ~I|J;
MPG_USA = MPGI;
MPG_Japan = MPGJ;
MPG_Europe = MPGK;
[fI,xI] = ksdensityMPG_USA;
plotxI,fI,'':b''
hold on
[fJ,xJ] = ksdensityMPG_Japan;
plotxJ,fJ,''r''
[fK,xK] = ksdensityMPG_Europe;
plotxK,fK,''--g''
legend''USA'',''Japan'',''Europe''
xlabel''样本'';ylabel''核密度估计''
hold off
图33数据频率直方图
图34不同窗宽下的核密度估计曲线
图35不同参数设置核密度估计曲线
图36各个国家汽车省油核密度估计曲线图
3.5统计作图用图形表达样本数据的统计特征具有生动、直观的特点,MATLAB提供了多种常用的统计图形函数来完成统计图绘制。3.5.1直方图直方图又称质量分布图,它是表示资料变化情况的一种主要工具。用直方图可以解析出资料的规则性,比较直观地看出产品质量特性的分布状态,对于资料分布状况一目了然,便于判断其总体质量分布情况。在制作直方图时,牵涉统计学的概念,首先要对资料进行分组,因此如何合理分组是其中的关键问题。按组距相等的原则进行的两个关键数位是分组数和组距,是一种几何形图表,它是根据从生产过程中收集来的质量数据分布情况,画成以组距为底边、以频数为高度的一系列连接起来的直方型矩形图。在MATLAB中,提供了hist函数用于实现绘制直方图。函数的调用格式如下。histx: 表示把矩阵x中的数据等距地划分为10个区间进行统计,并将每一区间内的数据个数作为返回值矢量的元素,最后画出10个柱形,如果x是矩阵,则将矩阵x的每一列进行统计。histx,nbins: 表示nbins是一个常量且指定了统计的区间个数。histx,xbins: 表示x是要统计的数据,nbins为一个矢量,矢量xbins的长度指定了统计的区间数,并以该矢量的各元素为中心进行统计。histax,: 表示在ax指定的坐标系中画出直方图。counts = hist: 表示只返回数据的频数。[counts,centers] = hist: 表示返回矢量counts和centers,分别表示频数和各个区间的位置。【例321】利用函数hist绘制randn概率分布图。其MATLAB代码编程如下:
clear all;
x = randn1000,1;
subplot3,1,1
xbins1 = -4:4;
histx,xbins1
subplot3,1,2
xbins2 = -2:2;
histx,xbins2
subplot3,1,3
xbins3 = [-4 -2.5 0 0.5 1 3];
histx,xbins3运行程序,效果如图37所示。
图37随机数据的直方图
3.5.2频数表在观察值个数较多时,为了解一组同质观察值的分布规律和便于指标的计算,可编制频数分布表,简称频数表。频数表是统计描述中经常使用的基本工具之一。1. 频数分布的特征由频数表可看出频数分布的两个重要特征: 集中趋势(Central Tendency)和离散程度(Dispersion)。身高有高有矮,但多数人身高集中在中间部分组段,以中等身高居多,此为集中趋势; 由中等身高到较矮或较高的频数分布逐渐减少,反映了离散程度。对于数值变量资料,可从集中趋势和离散程度两个侧面去分析其规律性。2. 频数分布的类型频数分布有对称分布和偏态分布之分。对称分布是指多数频数集中在中央位置,两端的频数分布大致对称。偏态分布是指频数分布不对称,集中位置偏向一侧,若集中位置偏向数值小的一侧,称为正偏态分布; 集中位置偏向数值大的一侧,称为负偏态分布,如冠心病、大多数恶性肿瘤等慢性病患者的年龄分布为负偏态分布。临床上正偏态分布资料较多见。不同的分布类型应选用不同的统计分析方法。3. 频数表的用途频数表可以揭示资料分布类型和分布特征,以便选取适当的统计方法; 便于进一步计算指标和统计处理; 便于发现某些特大或特小的可疑值。4. MATLAB实现在MATLAB中,提供了tabulate函数用于绘制频数表。函数的调用格式如下。tb1 = tabulatex: 表示对矢量x中的数据绘制频数表,返回值tb1的第1列是矢量x中的唯一值,第2列是每一个值出现的次数,第3列是每一个值出现的百分比例,如果x是一个数值型数组,则tb1是一个数值型矩阵,如果x的每一个元素都是非负整数,则tb1包含0到不包含在x中的从1到maxx的整数; 如果x是一个分类变量、字符数组或字符串单元数组,则tb1是一个单元数组。tabulatex: 表示不返回频数表。例如,在命令窗口中输入:
clear all;
tb1=tabulate[1 2 4 4 3 4]运行程序,输出如下:
tb1 =
1.00001.000016.6667
2.00001.000016.6667
3.00001.000016.6667
4.00003.000050.00003.5.3箱形图箱形图可以比较清晰地表示数据的分布特征,MATLAB提供了boxplot函数来绘制箱形图,它由5个部分组成:1 箱形上、下横线为样本的25%和75%分位数,箱形顶部和底部的差值为内四分位极值。2 箱形中间的横线为样本的中值,如果该横线没在箱形中央,则说明存在偏度。3 箱形向上或向下延伸的直线称为触须,如果没有异常值,样本的最大值为上触须的顶部,样本最小值为下触须的底部。默认情况下,距离箱形顶部或底部大于1.5倍同四分极值的值为异常值。4 图中顶部的加号表示该处数据为一异常值,该值的异常可能是输入错误、测量失误或系统误差引起的。5 箱形两侧的V形槽口对应于样本中值的置信区间。默认情况下,箱形图没有V形槽口。boxplot函数的调用格式如下。boxplotX: 对X中的每列数据绘制一个箱形图。boxplotX,notch: 当notch=1,得到一个有凹口的盒子图; notch=0,得到一个矩形箱形图。boxplotX,notch,''sym'': ''sym''为标记符号,缺省符号为 。boxplotX,nocth,''sym'',vert,whis: 参数vert控制箱形图水平放置还是垂直放置。当vert=0时,箱形图水平放置; 当vert=1时(缺省),箱形图垂直放置; whis定义虚线的长度,为内四分位间距(IQR)的函数(缺省情况为15*IQR)。如果whis=0,则box图用''sym''规定的标记显示箱子外所有的数据。【例322】根据参数的设置不同,绘制对应的样本的盒子图。其MATLAB代码编程如下:
clear all;
%产生正态分布的样本
%样本长度
N=1024;
x1=normrnd5,1,N,1;
x2=normrnd6,1,N,1;
x=[x1 x2];
%参数
figure1;
sym1=''*'';
notch1=1;%凹口
boxplotx,notch1,sym1;
figure2;
notch2=0;%矩形
boxplotx,notch2;
figure3;
vert=0;%水平
boxplotx,notch1,'' '',vert;运行程序,效果如图38~图310所示。
图38垂直、带凹口的盒子图
图39垂直、矩形的盒子图
图310水平、带凹口的盒子图
3.5.4经验累加分布图称函数
Fnx=0,xx1
ii=1fk,xix clear all;
rng default;%设置重复性
y = evrnd0,3,100,1;
[h,stats]=cdfploty
hold on
x = -20:0.1:10;
f = evcdfx,0,3;
plotx,f,''m:''
legend''经验分布曲线'',''理论上分布曲线'',''Location'',''NW''运行程序,输出如下,效果如图311所示。
h =
Line 具有属性:
Color: [0 0.4470 0.7410]
LineStyle: ''-''
LineWidth: 0.5000
Marker: ''none''
MarkerSize: 6
MarkerFaceColor: ''none''
XData: [1x202 double]
YData: [1x202 double]
ZData: [1x0 double]
显示所有属性
stats =
min: -10.5349
max: 4.4659
mean: -2.0350
median: -1.5294
std: 3.7186
图311经验累积分布函数图
3.5.5误差条图误差条图通常用于统计或科学数据,显示潜在的误差或相对于系列中每个数据标志的不确定程度。误差条图可以用标准差(平均偏差)或标准误差来表示,它们的区别如下。1 概念不同。标准差是离均差平方和平均后的方根,标准误差是标准误差定义为各测量值误差的平方和的平均值的平方根;2 用途不同。标准差与均数结合估计参考值范围,计算变异系数,计算标准误差等。标准误差用于估计参数的可信区间,进行假设检验等;3 它们与样本含量的关系不同。当样本含量n足够大时,标准差趋向稳定; 而标准误等随n的增大而减小,甚至趋于0。误差条形图类型的序列具有三个Y值。虽然可以手动将这些值分配给每个点,但在大多数情况下,是从其他序列中的数据来计算这些值。Y值的顺序十分重要,因为值数组中的每个位置都表示误差条形图上的一个值。在MATLAB统计工具箱中提供了errorbar函数用于绘制误差条图。函数调用格式如下。errorbarY,E: 表示绘制Y,以及对Y的每个元素绘制误差条,误差条的上半部分和下半部分都是长为Ei的对称条。errorbarX,Y,E: X、Y与E必须有相同的大小,如果X与Y是矢量,则误差条以(Xi,Yi)为中心,上下各长为Ei的线段条,如果X与Y是矩阵,则误差条是以(Xi,j,Yi,j)为中心,上下各长为Ei,j的线段条。errorbarX,Y,L,U: 表示由L和U指定误差条的上下界,在此X、Y、L、U必须有相同的长度,如果X是矢量,则在误差条是以(Xi,Yi)为中心,下长为Li上长为Ui的线段条,如果X是矩阵则误差条是以(Xi,j,Yi,j)为中心,下长为Li,j上长为Ui,j的线段条。errorbar,LineSpec: 表示由字符串LineSpec指定的线颜色、线类型来绘制误差条图。errorbarax,: 指定当前指定的坐标轴ax上绘制误差条图。h = errorbar: 表示返回误差条对象的一个句柄h。【例324】对所给定数据绘制其误差条图。其MATLAB代码编程如下:
clear all;
X = 0:pi10:pi;
Y = sinX;
E = stdY*onessizeX;
errorbarX,Y,E
xlabel''数据'';ylabel''误差条图'';运行程序,效果如图312所示。
图312误差条图
3.5.6交互等值线图交互等值线图即是指既可用代码实现绘图也可以通过手工来实现绘图。在MATLAB统计工具箱中提供了fsurfht函数用于绘制交互等值线图。函数的调用格式如下。fsurfhtfun,xlims,ylims: 表示生成由变量fun指定的交互式等值线图,x轴的限制由xlims=[xmin,xmax]来指定,y轴的限制由ylims=[ymin,ymax]来指定。fsurfhtfun,xlims,ylims,p1,p2,p3,p4,p5: 表示允许函数fun提供5个选项参数,函数fun的前面两个变量分别为x轴变量和y轴变量。图中有垂直参照线和水平参照线,两者的交点对应于当前点的x值和y值,可以通过拖拉这些带点的白色参考线来查看计算的z值(在图形上方)。另外,也可以通过在x轴和y轴的文本框中输入z值来得到指定的z值。【例325】绘制由gas.mat文件提供数据的高斯似然函数图形。其MATLAB代码编程如下:
function z = gauslikemu,sigma,p1
n = lengthp1;
z = onessizemu;
for i = 1:n
z = z .* normpdfp1i,mu,sigma;
end调用fsurfht函数绘制高斯似然函数图形,代码为
clear al;
load gas
fsurfht''gauslike'',[112 118],[3 5],price1%求似然函数图形
mumax = meanprice1%求price1的均值
sigmamax = stdprice1*sqrt1920 %求price1的标准差运行程序,输出如下,效果如图313所示。
mumax =
115.1500
sigmamax =
3.7719
图313交互等值线图
3.5.7散点图散点图(Scatter Diagram)是指在回归分析中,数据点在直角坐标系平面上的分布图。其是表示因变量随自变量而变化的大致趋势,据此可以选择合适的函数对数据点进行拟合。在MATLAB中,提供了gscatter函数用于实现散点图的绘制。函数的调用格式如下。gscatterx,y,group: 表示创建x和y的散点图,用group进行分组,其中x和y是矢量,且它们具有相同的大小,group可以是矢量、字符串数组或字符串单元数组,具有相同group值的点分在一组,在图中用相同的标记和颜色来表示,另外,group可以是包含一些分组变量(如[G1,G2,G3])的单元数组。gscatterx,y,group,clr,sym,siz: 表示指定每组的颜色、标记类型和大小,默认是,clr=''bgrcmyk'',sym是可以被函数plot识别的字符串数组,其默认值为.,siz是数组大小组成的矢量,其默认值由''DefaultLineMarkerSize''属性指定。gscatterx,y,group,clr,sym,siz,doleg: 表示由doleg指定是否在图中显示图例,当deleg=''on''时,表示显示图例,当doleg=''off''时表示不显示图例,默认值为''on''。gscatterx,y,group,clr,sym,siz,doleg,xnam,ynam: 表示由xnam和ynam指定x轴和y轴的名称,如果x和y的输入为简单的变量名,而且xnam和ynam被忽略,则函数gscatter用变量名标示坐标轴。h=gscatter: 表示返回图中直线的句柄数组。【例326】比较三种不同类型汽车的重量和里程数。其MATLAB代码编程如下:
clear all;
%装载数据
load carsmall
%比较不同类型汽车的重量和里程数
gscatterWeight,MPG,Model_Year,'''',''xos'';
xlabel''重量'';
ylabel''里程数'';运行程序,效果如图314所示。
图314三种不同类型汽车的重量和里程数散度图
由图314可以看出,1982年生产的汽车的里程数和重量明显区别于其他两种汽车。3.5.8最小二乘拟合线最小二乘法(又称最小平方法)是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据,并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用于曲线拟合。其他一些优化问题也可通过最小化熵或最大化熵用最小二乘法来表达。在MATLAB中,提供了lsline函数用于添加最小二乘拟合线。函数的调用格式如下。lsline: 表示在当前轴中每一直线对象上添加最小二乘直线。lslineax: 在指定的坐标轴ax中添加最小二乘拟合线。h=lsline: 返回直线对象的句柄h。【例327】利用lsline函数绘制最小二乘拟合线。其MATLAB代码编程如下:
clear all;
x = 1:10;
rng default;%设置重复性
figure;
y1 = x randn1,10;
scatterx,y1,25,''b'',''*''
hold on
y2 = 2*x randn1,10;
plotx,y2,''mo''
y3 = 3*x randn1,10;
plotx,y3,''rx:''运行程序,在一个图形中得到3条拟合曲线,效果如图315所示。
图315三条拟合线
在拟合曲线上添加最小二乘线,在命令窗口中输入:
lsline运行程序,效果如图316所示。
图316添加最小二乘直线
3.5.9正态概率图正态概率图用于检查一组数据是否服从正态分布,是实数与正态分布数据之间函数关系的散点图。如果这组实数服从正态分布,正态概率图将是一条直线。通常,概率图也可以用于确定一组数据是否服从任一已知分布,如二项分布或泊松分布。在MATLAB统计工具箱中提供了normplot函数用于绘制图形化正态性检验的正态概率图。函数调用格式如下。h = normplotX: 显示数据X的正态概率图,如果X为矩阵,则为X的每一列生成一条直线,该图中的样本数据用图形标记 显示,并在图中添加X中每列数据14和34处的连线,该线可以看做样本次序统计量的稳健性直线拟合,它可帮助评价数据的线性特征,如果数据源于正态分布,则图形呈现直线形,否则为曲线。【例328】利用normplot函数对给定的正态数据绘制概率图。其MATLAB代码编程如下:
clear al;
%生成正态分布数据
M=100;N=1;
x=normrnd0,1,M,N;
%生成均匀分布
y=randM,N;
z=[x,y];
%绘制正态概率图
h=normplotz;
xlabel''数据'';ylabel''概率'';
title''正态概率图'';
legend''正态分布数据'',''均匀分布数据'';
grid on;运行程序,效果如图317所示。
图317正态概率分布图
在正态概率图中有三个图形元素: 号表示每一个样本点数值的经验概率; 实线连接了数据的第25个和第75个百分点,表示一个线性拟合; 点画线将实线延伸到样本的两端。在正态概率图中,如果所有的样本点都在直线附近,则假设样本服从正态分布是合理的; 否则,如果样本不是正态分布的,则 号构成了一条曲线。通过观察图317中的两种不同分布样本的概率图可以验证这一点。3.5.10QQ图由两个样本的分位数绘制成的效果图称为QQ图,QQ图亦称为分位数图。在MATLAB中提供了qqplot函数用于绘制QQ图,其调用格式如下。qqplotX: 显示一个分位数分位数图。如果绘制分位数图的样本X源于正态分布,则绘制的QQ图近似于直线。qqplotX,Y: 显示两个样本的分位数分位数。如果两个样本来源于同一分布,那么,图中的曲线为直线。如果X与Y为乱阵,则为它们的每列数据绘制单独的曲线。图中样本数据以 符号表示,并将位于第一分位数和第三分位数间的数据拟合绘制成一条线(这是两个样本顺序统计量的鲁棒性拟合)。此线外推到样本数据的两端,以帮助用户评估数据的线性程度。qqplotX,Y,pvec: 函数可在pvec矢量中规定分位数。h=qqplotX,Y,pvec: 返回线段的句柄值h。【例329】绘制样本的QQ图。其MATLAB代码编程如下:
clear all;
%生成正态分布数据
M=100;N=1;
x=normrnd0,1,M,N;
%生成均匀分布
y=randM,N;
z=[x,y];
%绘制QQ图
subplot221;
h1=qqplotz;
xlabel''标准正态样本'';ylabel''输入样本'';title''QQ图'';
legend''正态分布数据'',''均匀分布数据'';
grid on;
%生成两个正态分布样本
x=normrnd0,1,100,1;
y=normrnd0.5,2,50,1;
subplot222
h2=qqplotx,y;
xlabel''输入样本x'';ylabel''输入样本y '';title''QQ图'';
grid on;
%生成两个不同分布的样本
x=normrnd5,1,100,1;
y=weibrnd2,0.5,100,1;
subplot223
h3=qqplotx,y;
xlabel''输入样本x'';ylabel''输入样本y'';title''QQ图'';
grid on;
subplot224
%生成一个正态分布的样本
x=normrnd10,1,100,1;
subplot224
qqplotx;
xlabel''输入样本x'';ylabel''输入样本x'';title''QQ图'';
grid on;运行程序,效果如图318所示。
图318QQ图
3.5.11帕累托图帕累托图又叫排列图、主次图,是按照发生频率大小顺序绘制的直方图,表示有多少结果是由已确认类型或范畴的原因所造成。它是将出现的质量问题和质量改进项目按照重要程度依次排列而采用的一种图表。可以用来分析质量问题,确定产生质量问题的主要因素。按等级排序的目的是指导如何采取纠正措施,项目班子应首先采取措施纠正造成最多数量缺陷的问题。从概念上说,帕累托图与帕累托法则一脉相承,该法则认为相对来说数量较少的原因往往造成绝大多数的问题或缺陷。帕累托法则往往称为二八原理,即百分之八十的问题是百分之二十的原因所造成的。帕累托图在项目管理中主要用来找出产生大多数问题的关键原因,用来解决大多数问题。在帕累托图中,不同类别的数据根据其频率降序排列的,并在同一张图中画出累积百分比图。帕累托图可以体现帕累托原则: 数据的绝大部分存在于很少类别中,极少剩下的数据分散在大部分类别中。这两组经常被称为至关重要的极少数和微不足道的大多数。帕累托图能区分微不足道的大多数和至关重要的极少数,从而方便人们关注于重要的类别。帕累托图是进行优化和改进的有效工具,尤其应用在质量检测方面。在MATLAB的统计工具箱中提供了pareto函数用于绘制帕累托图。函数调用格式如下。paretoY: 将矢量Y中的每个元素按元素数值递减顺序绘成直方条,并以其Y中的索引号进行标记。各直方条上方的折线显示累积频率。paretoY,names: 以字符串names中的名称对Y中相应的元素所绘的直方条进行标记。paretoY,X: 根据给定的X值对直方条进行标记。H = pareto: 返回帕累托图的句柄值H。【例330】根据给定的一组生产工的生产情况,绘制帕累托图。其MATLAB代码编程如下:
clear all;
%给定生产力
codelines = [200 120 555 608 1024 101 57 687];
%生产工名
coders ={''Fred'',''Ginger'',''Norman'',''Max'',''Julia'',''Wally'',''Heidi'',''Pat''};
paretocodelines, coders%绘制帕累托图
title''生产力制帕累托图''
xlabel''数据'';ylabel''效果图'';运行程序,效果如图319所示。
图319帕累托图
3.5.12频率直方图将样本观测值x1,x2,,xn从小到大排序并去除多余的重复值,得到x1 clear all;
y = exprnd10,50,1; %随机故障次数
d = exprnd20,50,1; %随机丢失次数
t = miny,d;%最低次数
censored = yd;%观察是否受失败
%计算经验分布并绘制频率直方图
[f,x] = ecdft,''censoring'',censored;
ecdfhistf,x
setgetgca,''Children'',''FaceColor'',[.8 .8 1]
hold on
% Superimpose a plot of the known population pdf
xx = 0:.1:maxt;
yy = exp-xx1010;
plotxx,yy,''r-'',''LineWidth'',2
hold off运行程序,效果如图320所示。
图320随机数据的频率直方图
|