Matlab在统计中的应用

数理统计建模

Matlab在统计中的应用山西财经大学应用数学学院高崇山一、概率分布及有关函数字符分布unif均匀分布exp指数分布norm正态分布chi2chi方分布tt分布ff分布bino二项分布poisspoisson分布字符功能调用格式pdf概率密度namepdf(x,参数）cdf分布函数inv逆概率分布stat均值与方差ran随机数生成调用格式为：分布命令符功能命令符(x,参数)[m,v]=分布stat(x,参数)rand产生[0,1]上的随机数，randn产生标准正态分布随机数。y=normpdf(x,mu,sigma)或y=pdf(‘norm’,x,mu,sigma)%正态分布N(mu,sigma2)在x处的概率密度；y=normcdf(x,mu,sigma)或y=cdf(‘norm’,x,mu,sigma)%正态分布N(mu,sigma2)在x处的分布函数；y=norminv(alpha,mu,sigma)%正态分布N(mu,sigma2)在对应于alfa的分位数。即[m,v]=normstat(mu,sigma)%正态分布N(mu,sigma2)的期望和方差；y=exprnd(lamda,[m,n])或random(‘exp’,lamda,[m,n])%产生一个m×n的服从参数为lamda的指数分布的随机矩阵二、描述性统计描述性统计就是搜集、整理、加工和分析统计数据，使之系统化，以显示出数据的趋势、特征和数量关系。函数名描述函数名描述max求向量或矩阵列的最大值sort升序排列min求向量或矩阵列的最小值sum求向量或矩阵列的和mean求向量或矩阵列的平均值cumsum累计求和median求向量或矩阵列的中间值cov求协方差std求标准差corrcoef求相关系数var求方差kurtosis计算样本峰度rang样本极差skewness计算样本偏度2.1样本均值mean和中值median它们都是样本数据在数据分布线上中心位置的度量.A=[1244;3466;5688;5688];mean(A)%计算矩阵每列的均值，相当于mean(A,1)mean(A,2)%计算矩阵每行的均值median(A)%计算矩阵每列的中值〔中位数〕，相当于median(A,1)median(A,2)%计算矩阵每行的中值〔中位数〕2.2方差var、标准差std、极差range和协方差cov它们都是描述样本中的数据偏离其中心值的程度X=rand(4,5);std(X)%计算矩阵X每列的标准差var(X)%计算矩阵X每列的方差range(X)%计算矩阵X每列的极差cov(X)%计算协方差var(X)=diag(cov(X))’std(X)=sqrt(diag(cov(X)))’X假设为向量，cov(X)=var(X)；假设X为矩阵，X的每一列表示一个变量而行元素为观察值。对于二维随机向量(X,Y)，x为X的观察值，y为Y的观察值(x,y为同维向量)，那么有：cov(x,y)=cov([x,y])2.3百分位数及其图形描述百分位数(percentile)是把数据按从小到大的顺序排列后，位于p%位置的值称为第p百分位数。第25百分位数由叫做四分之一分位数(下四分位数)，75百分位数由叫做四分之三分位数(上四分位数)，第50百分位数就是median中数。最小值是第0百分位数，最大值是第100百分位数。百分位数是用于反映样本数据形态信息的数据统计量，它也可以刻划数据的位置和散布特征。Y=prctile(X,p)返回样本X中大于p%(0<p<100)的值。如果X是向量，那么返回X中p百分位数，假设X为矩阵，那么返回一个关于每列元素的p百分位数行向量。注意：p也可以是一个向量，此时，返回一组百分位数。eg.x=100*rand(1,10),y=prctile(x,0:10:100),subplot(1,2,1),boxplot(x),subplot(1,2,2),bar(x)x=[61.543279.193792.181373.820717.626640.570693.547091.690441.027089.3650]y=[17.626629.098640.798851.285167.682076.507284.279390.527791.935992.864193.5470]2.4相关系数

相关系数反映两个随即变量之间线性相依程度的变量。R=corrcoef(X)R=corrcoef(x,y)含义同协方差。同协方差之间有如下关系2.5样本峰度和偏度偏度描述的是分布的对称性，其定义为：当f>0时，表示数据在均值右边的比左边的多；f<0正好相反；f接近于0，那么表示分布是对称的。峰度描述的是分布曲线的陡缓程度，定义为：它是以正态分布为标准，比较两侧极端数据分布的情况的指标。g较大，那么表示样本中有许多远离均值的数据。上述公式中，s是样本标准差。f=skewness(X)g=kurtosis(X)三、参数估计函数名功能binofit二项分布参数估计expfit指数分布的参数估计normfit正态分布的参数估计poissfitpoisson分布的参数估计unifit均匀分布的参数估计mle最大似然估计……1〕fit函数的调用方法类似，以正态分布说明之。格式为：[mu,sigma,muci,sigmaci]=normfit(x,alpha)说明：x是样本(矩阵或向量)；alpha是显著性水平(默认值为0.05)；mu是总体均值的点估计值；sigma是总体方差的点估计值；muci是总体均值的区间估计；sigmaci是总体方差的区间估计。2〕mle的调用格式[phat,pci]=mle(‘dist’,data,alpha,pl)说明：dist是所给的分布名(如:norm,exp…)；data是样本数据；alpha为可选项，表示显著性水平；pl仅用于二项分布，表示试验的次数；phat为返回的点估计值；pci为返回相应置信区间。eg.x=normrnd(2,4,100,1);[mu,sigma,muci,sigmaci]=normfit(x)[mu,sigma,muci,sigmaci]=normfit(x,0.1)[phat,pci]=mle('norm',x)[phat,pci]=mle('norm',x,0.1)四、假设经验4.1单个样本的t检验功能：进行样本均值的t检验格式：h=ttest(x,m)；%在0.05的显著性水平下进行t检验，以确定在标准差未知的情况下取自正态分布的均值是否为m，假设输出h=0，那么接受零假设，h=1那么否认零假设；h=ttest(x,m,alpha)；%alpha为给定的显著性水平。[h,sig,ci]=ttest(x,m,alpha,tail)%假设原假设为μ=μ0，那么取tail=0(可省略)；假设原假设为μ>μ0，那么取tail=1；假设原假设为μ<μ0，那么取tail=-1。sig表示在零假设下，样本均值出现的概率，当sig>alpha时不能否认零假设，一般sig越大零假设越可信。ci为均值真值的1-alpha置信区间。4.2单个样本的z检验功能：在给定方差的条件下进行z检验格式：h=ztest(x,m,sigma)%sigma正态总体的标准差，alpha=0.05h=ztest(x,m,sigma,alpha)[h,sig,ci]=ztest(x,m,sigma,alpha,tail)4.3两个样本的t检验功能：两个服从正态总体样本均值差异的t检验(σ12=σ22均未知)格式：[h,significance,ci]=ttest2(x,y)%默认alpha=0.05[h,significance,ci]=ttest2(x,y,alpha)[h,significance,ci]=ttest2(x,y,alpha,tail)五、统计绘图5.1box图boxplot(X,notch,’sym’,vert,whis)X：样本数据；notch=1有切口,notch=0无切口,(默认notch=0)；‘sym’野值标记符号,默认‘+’；vert=0,box图是水平放置,vert=1是垂直放置(默认)。whis定义虚线的长度，一般用缺省值.eg:x1=normrnd(4,1,200,1);x2=normrnd(8,1,200,1);x3=normrnd(6,2,200,1);x=[x1,x2,x3];boxplot(x,1)图见下页说明：1.盒子的上下两条线分别为样本的75%和25%分位线，中间为样本中位数；2.虚线表示样本的其余局部，位于盒子的上下两侧；3.‘+’表示野值〔奇异值〕，位于虚线的上方和下方；4.‘切口’表示样本中位数的置信区间。默认状态下无切口。5.2正态概率图正态概率图用于判断样本数据是否服从正态分布。格式：normplot(X)X：数据.假设X为矩阵，那么为X的每列显示一条线。图形以符号‘+’显示样本数据。如果数据服从正态分布，那么图形呈现直线，否那么会表现不同程度的曲线。x1=normrnd(4,1,200,1);x2=normrnd(8,1,200,1);x3=normrnd(6,2,200,1);x=[x1,x2,x3];normplot(x,1)5.3分位数—分位数图分位数—分位数图用于比较两个样本的分布.格式:qqplot(X,Y,pvec)其中，X,Y分别是两个样本的数据。如果两个样本来自同一分布，那么绘制的曲线为直线。假设X,Y为矩阵，那么为他们每一列显示一条直线。图形以符号‘+’显示样本数据。参数pvec是可选项，用于规定分位数。x1=normrnd(4,1,200,1);x2=normrnd(0,1,200,1);qqplot(x1,x2)六、分布检验6.1Jarque-bera检验该检验评价X服从未知均值和方差的正态分布的假设是否成立。该检验基于X的样本偏度和峰度。对于正态分布数据，偏度接近于0，峰度接近于3。Jarque-bera检验就是确定样本偏度、峰度是否与它们的期望值相差较远。功能：测试数据对正态分布的拟合程度。格式:h=jbtest(X)%当h=1,拒绝X服从正态分布;否那么h=0。(默认alpha=0.05)h=jbtest(X,alpha)[h,P,jbstat,cv]=jbtest(X,alpha)%P为检验的p值，jbstat为检验的统计量，cv为确定是否拒绝零假设的临界值。当jbstat<cv时，同样接受零假设。注意：该检验不能用于小样本的检验，只能用于大样本。对于小样本，用lillietest检验较适宜。6.2Lilliefors检验该检验评价X服从未知均值和方差的正态分布的假设是否成立，对应的备择假设为X不服从正态分布。本检验比较X的经验分布与具有相同均值和方差的正态分布。格式：H=lillietest(X)%假设H=0，那么接受X服从正态分布；否那么，H=1。(默认alpha=0.05)H=lillietest(X,alpha)[H,P,LSTAT,CV]=lillietest(X,alpha)%P为检验的p值，通过在一系列由Lilliefors创立的表中进行插值得到；LSTAT为检验统计量的值；CV为确定是否拒绝零假设的临界值。如果LSTAT的值位于Lilliefors表之外，那么P返回NaN，但H显示是否拒绝假设。当LSTAT>CV时，同样拒绝零假设。七、回归分析7.1多元线性回归分析数学模型一元回归模型为：y=β0+β1x+ε其中ε服从N(0,1)多元回归模型为：

y=β0+β1x1+β2x2+…+βmxm+ε

其中ε服从N(0,σ2)

回归问题就是求出xi的系数βi，并求出误差σ2的估计，回归系数β的区间估计和假设检验，模型的有效性检验及对给定的x做出y的预测。预测分点预测和区间预测：其中点预测将x代入模型中即可，区间预测需要编一个小程序。命令为：b=regress(y,x)[b,bint,r,rint,s]=regress(y,x,alpha)说明：

输入y(因变量，列向量)；x(第一列全为1，第二列为x1的观察值，第三列为x2的观察值，…)；alpha为显著性水平α(默认值为0.05)；输出b为(β0,β2,…β2,βm)的估计值；bint是(β0,β2,…β2,βm)的置信区间；r是残差(观察值与预测值之差，为列向量，主要用于探测模型假设的合理性)，rint是残差的置信区间；s包含3个统计量：第一个是决定系数R2(其值越大，说明自变量对因变量的所起的作用也越大，但无明确界限说明模型是否有效)，第二个是F值，第三个是F(1,N-2)分布大于F值的概率p，p<alpha时，回归模型有效。y=[144 215 138 145 162 142 170 124 158 154162 150 140 110 128 130 135 114 116 124136 142 120 120 160 158 144 130 125 175];x=[39 47 45 47 65 46 67 42 67 5664 56 59 34 42 48 45 18 20 1936 50 39 21 44 53 63 29 25 69];n=length(y);X=[ones(n,1)x'];[b,bint,r,rint,s]=regress(y',X);b,bint,s,s2=sum(r.^2)/(n-2)rcoplot(r,rint)pausey=[y(1)y(3:30)];x=[x(1)x(3:30)];n=length(y);X=[ones(n,1)x'];[b,bint,r,rint,s]=regress(y',X);b,bint,s,s2=sum(r.^2)/(n-2)rcoplot(r,rint)pausey0=b(1)+b(2)*50;%预测y(x=50)xb=mean(x);sxx=sum((x-xb).^2);a=sqrt((50-xb)^2/sxx+1/n+1);t=tinv(0.975,n-2);d=t*a*sqrt(s2);y1=y0-d;y2=y0+d;%预测y(x=50)区间〔t分布〕[y0y1y2]d1=norminv(0.975)*sqrt(s2);y3=y0-d1;y4=y0+d1;[y0y3y4]%预测y(x=50)区间〔N分布〕7.2多项式回归1)多项式曲线拟合多项式回归的模型为：p(x)=p1xn+p2xn-2+…+pnx+pn+1格式为：[p,s]=polyfit(x,y,n)其中：n为拟合次数；x,y分别是自变量和因变量；s是一个矩阵，用于polyval函数，可进行预测的误差估计；p为系数向量(p1,p2,…,pn,pn+1)的估计值。2〕多项式预测和置信区间的评估[y,delta]=polyconf(p,x,s,alpha)其中：p,s是拟合输出的结果，x是要预测的点，alpha是置信度，输出的是1-alpha的置信区间y±delta.说明：命令polytool(x,y,n,alpha)作用类似于polyfit.他是一个交互式画面。eg.y=[1035624108410521015106670496099010508391030985855];x=[6.00002.50007.50008.500010.00007.00003.000011.50005.50006.50004.00009.000011.000012.5000];plot(x,y,'+'),pausex2=x.^2;X=[ones(14,1)x'x2'];[b,bi,r,ri,s]=regress(y',X);b,bi,s,pausexx=2:.1:13;yy=b(1)+b(2)*xx+b(3)*xx.^2;plot(x,y,'+',xx,yy),grid