matlab在数理统计中的应用1
matlab在数理统计中的应用1
一、直方图与经验分布函数图的绘制
hist(A,n) —— 对矩阵A按列作统计频数直方图, n为条形图的条数
hist(A,x0) —— 对矩阵A按列作以向量x0为划分区间中点的频数直方图 ni=hist(A,n) 或 ni=hist(A,x0) —— 对矩阵A按列得各划分区间内的统计频数 注意: 当A为向量时, 上述所有命令直接作用在向量上, 而不是列优先. bar(x,y) —— 绘制分别以x和y为横纵坐标的二维条形图 cdfplot(x) —— 绘制样本x的经验分布函数图
[Fn,x0]=ecdf(x) —— 得到样本x的经验分布函数值Fn, 当x中有m个不同的数 (记为向量x0) 时, 则Fn的个数为m+1个 例如:
>> x = [6 4 5 3 6 8 6 7 3 4]; >> [Fn,x0]=ecdf(x)
Fn = 0 0.2000 0.4000 0.5000 0.8000 0.9000 1.0000 x0 = 3 3 4 5 6 7 8 >>cdfplot(x)
【例1-1】
在齿轮加工中, 齿轮的径向综合误差是个随机变量, 今对200件同样的齿轮进行测量, 测得的数值 (mm) 如下, 求作的频率密度直方图, 并作出的经验分布函数图形.
16 25 19 20 25 33 24 23 20 24 25 17 15 21 22 26 15 23 22 24 20 14 16 11 14 28 18 13 27 31 25 24 16 19 23 26 17 14 30 21 18 16 18 19 20 22 19 22 18 26 26 13 21 13 11 19 23 18 24 28 13 11 25 15 17 18 22 16 13 12 13 11 09 15 18 21 15 12 17 13 14 12 16 10 08 23 18 11 16 28 13 21 22 12 08 15 21 18 16 16 19 28 19 12 14 19 28 28 28 13 21 28 19 11 15 18 24 18 16 28 19 15 13 22 14 16 24 20 28 18 18 28 14 13 28 29 24 28 14 18 18 18 08 21 16 24 32 16 28 19 15 18 18 10 12 16 26 18 19 33 08 11 18 27 23 11 22 22 13 28 14 22 18 26 18 16 32 27 25 24 17 17 28 33 16 20 28 32 19 23 18 28 15 24 28 29 16 17 19 18
·编写函数文件example1_6_1.m:
function [ni,frequency]=example1_6_1(Y,k) %Y为试验数据 %k为分组参数 n=length(Y); m=min(Y); M=max(Y);
%(1)下面作频率密度直方图 delta=(M-m)/k; x=m:delta:M; ni=hist(Y,x); frequency=ni/n;
density=frequency/delta;
bar(x,density) %作频率密度直方图 pause
%(2)下面作经验分布函数图
·编写命令文件example1_6_2.m:
F=[16 25 19 20 25 33 24 23 20 24 25 17 15 21 22 26 15 23 22 24.... 20 14 16 11 14 28 18 13 27 31 25 24 16 19 23 26 17 14 30 21.... 18 16 18 19 20 22 19 22 18 26 26 13 21 13 11 19 23 18 24 28.... 13 11 25 15 17 18 22 16 13 12 13 11 09 15 18 21 15 12 17 13.... 14 12 16 10 08 23 18 11 16 28 13 21 22 12 08 15 21 18 16 16.... 19 28 19 12 14 19 28 28 28 13 21 28 19 11 15 18 24 18 16 28.... 19 15 13 22 14 16 24 20 28 18 18 28 14 13 28 29 24 28 14 18.... 18 18 08 21 16 24 32 16 28 19 15 18 18 10 12 16 26 18 19 33.... 08 11 18 27 23 11 22 22 13 28 14 22 18 26 18 16 32 27 25 24.... 17 17 28 33 16 20 28 32 19 23 18 28 15 24 28 29 16 17 19 18]; [ni,frequency]=example1_6_1(F,9) ·运行命令文件example1_6_2.m: >> example1_6_2
ni = 5 16 21 34 44 25 23 22 4 6
frequency = 0.0250 0.0800 0.1050 0.1700 0.2200 0.1250 0.1150 0.1100 0.0200 0.0300
图1-1(a) 图1-1(b)
二、常见的概率分布
表1.1 概率分布分类表
三、MATLAB为常见分布提供的五类函数 1) 概率密度函数(pdf); 2) (累积)分布函数(cdf); 3) 逆(累积)分布函数(icdf); 4) 随机数发生器(random); 5) 均值和方差(stat). 1、概率密度函数
注意: Y=normpdf (X, MU, SIGMA)的SIGMA是指标准差 , 而非 . 【例1-2】绘制标准正态分布的概率密度图.
y=normpdf(x,0,1); plot(x,y)
title('N(0,1)的概率密度曲线图')
图1-2
2、累积分布函数
表1.3 累积分布函数(cdf)
常见的概率分布(matlab作图)
一、常见的概率分布
表1.1 概率分布分类表
二、MATLAB为常见分布提供的五类函数 1) 概率密度函数(pdf); 2) (累积)分布函数(cdf); 3) 逆(累积)分布函数(icdf); 4) 随机数发生器(random); 5) 均值和方差(stat). 1、概率密度函数
表1.2 概率密度函数(pdf)
注意【例1-2】绘制标准正态分布的概率密度图. x=-4:0.1:4;
y=normpdf(x,0,1);
title('N(0,1)的概率密度曲线图')
图1-2
2、累积分布函数
表1.3 累积分布函数(cdf)
【例1-3】求服从标准正态分布的随机变量落在区间[-2, 2]上的概率. >> P=normcdf ([-2, 2]) ans = 0.0228 0.9772 >>P(2)-P(1) ans = 0.9545
3、逆累积分布函数 (用于求分位点)
【例1-4】(书P22例1.13) 求下列分位数: (i) ; (ii) ; (iii) ; (iv) . >> u_alpha=norminv(0.9,0,1) u_alpha = 1.2816 >> t_alpha=tinv(0.25,4) t_alpha = -0.7407
>> F_alpha=finv(0.1,14,10) F_alpha = 0.4772
>> X2_alpha=chi2inv(0.025,50) X2_alpha = 32.3574 4、随机数发生函数
表1.5 随机数发生函数(random)
5、均值和方差
表1.6 常见分布的均值和方差函数(stat)
注意: 三、常用的统计量
表1.7 常用统计量
说明:
(1) y=var(X) ——计算X中数据的方差. .
y=var(X, 1) —— , 得到样本的二阶中心矩 (转动惯量).
(2) C=cov(X) ——返回一个协方差矩阵, 其中输入矩阵X的每列元素代表着一个随机变量的观测值. 如果X为n×m的矩阵, 则C为m×m的矩阵.
(3) var(X)=diag(cov(X)), std(X)=sqrt(diag(cov(X))).
参数估计(matlab)
参数估计包含两种常用方式: 点估计和区间估计.
Matlab统计工具箱给出了常用概率分布中参数的点估计 (采用最大似然估计法) 与区间估计, 另外还提供了部分分布的对数似然函数的计算功能.
由于点估计中的矩估计法的实质是求与未知参数相应的样本的各阶矩,统计工具箱提供了常用的求矩函数(见第一章), 读者可根据需要选择合适的矩函数进行点估计.
表2.1 统计工具箱中的参数估计函数 (fit / like)
说明: 调用格式只罗列了其中的一种. 需另外说明的是:
(1) unifit和normfit的格式与其它函数均不同, 此二者要求左边的输出变量必须将参数或分别列出.
(2) binofit (x,n,alpha)根据试验成功的次数x和总的试验次数n, 对中的p进行最大似然估计, 同时返回置信度为100(1-alpha)%的置信区间pci.
【例2-1】(书P692.3) 使用一测量仪器对同一值进行了12次独立测量, 其结果为 (单位: mm)
232.50, 232.48, 232.15, 232.52, 232.53, 232.30, 232.48, 232.05, 232.45, 232.60, 232.47, 232.30 试用矩法估计测量的真值和方差 (设仪器无系统误差). ·编写命令文件exercise2_3.m: %P66_2.3 mu与sigma^2的矩估计
x=[232.50, 232.48, 232.15, 232.52, 232.53, 232.30,... 232.48, 232.05, 232.45, 232.60, 232.47, 232.30]; mu_ju=mean(x) sigma2_ju=var(x,1) ·运行命令文件exercise2_3.m: >> exercise2_3 mu_ju = 232.4025 sigma2_ju = 0.0255
【例2-2】(书P692.22) 随机地从一批零件中抽取16个, 测得长度 (单位: cm) 为:
2.14, 2.10, 2.13, 2.15, 2.13, 2.12, 2.13, 2.10, 2.15, 2.12, 2.14, 2.10, 2.13, 2.11, 2.14, 2.11
设零件长度的分布为正态的, 试求总体均值的90%的置信区间: (1)若 (cm); (2) 若未知. (1)·编写函数文件zestimate.m:
%P69_2.22(1)sigma已知时mu的区间估计 function muci=zestimate(x,sigma,alpha) n=length(x); xhat=mean(x);
u_alpha=norminv(1-alpha/2,0,1); delta1=sigma/sqrt(n)*u_alpha; muci=[xhat-delta1,xhat+delta1]; ·调用函数文件zestimate.m:
>> x=[2.14, 2.10, 2.13, 2.15, 2.13, 2.12, 2.13, 2.10, 2.15, 2.12, 2.14, 2.10, 2.13, 2.11, 2.14, 2.11]; >>sigma=0.01; >>alpha=0.1;
>>muci=zestimate(x,sigma,alpha) muci = 2.1209 2.1291 (2)·编写命令文件exercise2_22_2.m:
%P69_2.22(1)sigma未知时mu的区间估计
x=[2.14, 2.10, 2.13, 2.15, 2.13, 2.12, 2.13, 2.10, 2.15, 2.12, 2.14, 2.10, 2.13, 2.11, 2.14, 2.11]; alpha=0.1;
[muhat,sigmahat,muci,sigmaci]= normfit(x,alpha); muci ·运行命令文件exercise2_22_2.m: >> exercise2_22_2
muci = 2.1175 2.1325
【例2-3】(书P66例2.31) 对一批产品, 欲通过抽样检查其合格率. 若产品不合格率在5%以下, 则该批产品可出厂. 检验时要求结果具有0.95的置信水平. 今抽取产品100件, 发现不合格品有4件, 问这批产品能否出厂? >> [phat,pci]=binofit(4,100,0.05) phat = 0.0400
pci = 0.0110 0.0993
由于置信区间的上限超出了规定指标(不合格率在5%以下), 因此不能出厂.
假设检验1
假设检验分为两种: 参数假设检验与分布假设检验. 一、正态总体的参数假设检验 表3.1的说明:
对一个正态总体 , 抽取样本 ;
对两个正态总体 , , 且X与Y独立,分别抽取样本与 .
表3.1 正态总体的参数假设检验
二、非参数假设检验
表3.2 非参数假设检验
三、统计工具箱中的假设检验
表3.3 统计工具箱中的假设检验 (test / rank)
1、jbtest, lillietest与kstest的比较:
(1) jbtest与lillietest均是检验样本是否来自正态分布, 而kstest可检验样本来自任意指定的分布; (2) jbtest是利用偏度峰度来检验, 适用于大样本; 而对于小样本, 则用lillietest来检验;
(3) lillietest与kstest的检验原理均是用x的经验分布函数与一个有相同均值与方差的正态分布的分布函数进行比较, 不同的是lisllietest中正态分布的参数是由x估计得来, 而kstest中正态分布的参数是事先指定的. 2、kstest2对应于斯米尔诺夫检验. 3、命令说明:
(1) [h,sig,ci,zval] = ztest(x,mu0,sigma,alpha,tail) 对已知方差的单个总体均值进行Z检验. 进行显著性水平为的Z假设检验, 以检验标准差为的正态分布样本的均值与的关系. 并可通过指定tail的值来控制备择假设的类型. tail 的取值及表示意义如下:
tail=0 备择假设为 (缺省值); tail=1备择假设为 ;
tail= -1备择假设为 . (原假设则为 ) ·输出变量含义:
h——如果h=0, 则接受 ; 如果h=1, 则拒绝而接受备择假设 ; sig——Z的观察值在下较大或统计意义上较大的概率值; ci——方差未知时均值的的置信区间. zval——Z统计量的观测值. ·单边检验对应单侧区间估计.
(2) [h,sig,ci,tval] = ttest(x,mu0,alpha,tail) 格式调用中无“tval”这个输出变量, 但可加上此项. tval——包含两个结果: tstat表示t统计量的值; df表示t分布的自由度.
(3) [h,p,jbstat,cv] = jbtest(x,alpha) 对“单个总体服从正态分布(未指定均值和方差)”假设进行显著水平为的
Jarque-Bera检验. 此检验基于x的偏度与峰度. 对于真实的正态分布, 样本偏度应接近于0, 样本峰度应接近于3. Jarque-Bera检验通过统计量来判定样本偏度和峰度是否与它们的期望值显著不同. ·输出变量含义:
h——如果h=0, 则接受“ : 认为x来自正态总体”; 如果h=1, 则接受备择假设“ : 认为x不是来自正态总体”; p——检验的概率p-值; jbstat——检验统计量的值;
cv——判断是否拒绝原假设的关键值.
(4) [h,p,ksstat,cv] = kstest(x,cdf,alpha,tail) 对“x的总体服从由两列矩阵cdf指定的分布G”假设进行显著水平为的Kolmogorov-Smirnov检验. 矩阵cdf的第一列包含可能的x值, 第二列包含相应的理论累积分布函数值G(x0). 在可能的情况下, 应定义cdf使每一列包含x中的值. 如果cdf=[ ], kstest( )将使用标准正态分布. (5) [h,p,ksstat] = kstest2(x,cdf,alpha,tail) 对“两个样本来自同一连续分布”假设进行显著水平为的
Kolmogorov-Smirnov检验. 对于大容量的样本来说, p-值将很精确, 一般来说, 当样本容量N1和N2满足时, p-值即可认为是精确的.
(6) normplot(x) 绘出x中数据的正态检验概率图. 如果x是一个矩阵, 则对每一列绘出一条线. 图中样本数据用符号„+‟来表示, 叠加在数据上的实线是数据的第一个与第三个四分位点之间的连线 (为样本顺序统计量的鲁棒线性拟
合). 这条线延伸到样本数据的两端, 以便估计数据的线性度. 如果数据是来自一个正态分布, 则„+‟线近似地在一直线上. 一般地, 中间的点离直线位置的偏差不能过大, 两头的点的偏差可以允许大一些. 当中间的点离直线位置偏差太大时, 就认为x来自其它分布.
(7) qqplot(x,y) 绘出两样本的分位数-分位数图. 图中样本数据用符号„+‟来表示, 叠加在数据上的实线是各分布的第一个与第三个四分位点之间的连线 (为两个样本顺序统计量的鲁棒线性拟合). 这条线延伸到样本数据的两端以便估计数据的线性度. 如果两个样本来源于同一个分布, 则„+‟线近似地在一直线上.
qqplot(x) 绘出样本x的分位数-正态分布的理论分位数图. 如x为正态分布, 则„+‟线近似地在一直线上.
【例3-1】(例3.4) 一台包装机装洗衣粉, 额定标准重量为500g, 根据以往经验, 包装机的实际装袋重量服从正态分布 , 其中 g, 为检验包装机工作是否正常, 随机抽取9袋, 称得洗衣粉净重数据如下 (单位: g):
497 506 518 524 488 517 510 515 516
若取显著性水平 , 问这包装机工作是否正常? >> x=[497,506,518,524,488,517,510,515,516]; >> [h,sig,ci,zval]=ztest(x,500,15,0.01,0) h = 0 %接受
sig = 0.0432 % 为真条件下P( )的值
ci = 497.2320 522.9903 % 未知时的置信水平为0.95的双侧置信区间 zval = 2.0222 %Z统计量的值. 所以认为包装机工作正常. 【例3-2】(例3.5) 某部门对当前市场的价格情况进行调查. 以鸡蛋为例, 所抽查的全省20个集市上, 售价分别为 (单位: 元/500克)
3.05, 3.31, 3.34, 3.82, 3.30, 3.16, 3.84, 3.10, 3.90, 3.18, 3.88, 3.22, 3.28, 3.34, 3.62, 3.28, 3.30, 3.22, 3.54, 3.30.
已知往年的平均售价一直稳定在3.25元/500克左右, 在显著性水平下, 能否认为全省当前的鸡蛋售价明显高于往年?
>> x=[3.05,3.31,3.34,3.82,3.30,3.16,3.84,3.10,3.90,3.18,... 3.88,3.22,3.28,3.34,3.62,3.28,3.30,3.22,3.54,3.30]; >> [h,sig,ci,tval]=ttest(x,3.25,0.025,1) h = 1
sig = 0.0114
ci = 3.2731 Inf
tval = tstat: 2.4763 df: 19
所以认为全省当前的鸡蛋售价明显高于往年.
假设检验2
【例3-3】(例3.6) 某工厂生产某种电器材料. 要检验原来使用的材料与一种新研制的材料的疲劳寿命有无显著性差异, 各取若干样品, 做疲劳寿命试验, 所得数据如下 (单位: 小时) : 原材料: 40, 110, 150, 65, 90, 210, 270,
新材料: 60, 150, 220, 310, 350, 250, 450, 110, 175.
一般认为, 材料的疲劳寿命服从对数正态分布, 并可以假定原材料疲劳寿命的对数与新材料疲劳寿命的对数有相同的方差, 即可设 , . 在显著性水平下, 能否认为两种材料的疲劳寿命没有显著性差异? >> x=[40,110,150,65,90,210,270];
>> y=[60,150,220,310,350,250,450,110,175]; >> [h,sig,ci,tval]=ttest2(log(x),log(y),0.05,0) h = 0
sig = 0.1001
ci = -1.2655 0.1244
tval = tstat: -1.7609 df: 14
所以认为两种材料的疲劳寿命没有显著差异.
【例3-4】(例3.18) 对一台设备进行寿命试验, 记录10次无故障工作时间, 并从小到大排列得 420, 500, 920, 1380, 1510, 1650, 1760, 2100,2300, 2350 问此设备的无故障工作时间X的分布是否服从的指数分布( )? >> x=[420, 500, 920, 1380, 1510, 1650, 1760, 2100,2300, 2350]; >> [h,p,ksstat,cv]=kstest(x,[x',expcdf(x',1500)],0.05,0) h = 0
p = 0.2700 ksstat = 0.3015 cv = 0.4093
所以认为此设备的无故障工作时间X服从的指数分布. 【例3-5】(例3.20) 在20天内, 从维尼纶正常生产时生产报表上看到的维尼纶纤度 (表示纤维粗细程度的一个量) 的情况, 有如下100个数据:
1.36, 1.49, 1.43, 1.41, 1.37, 1.40, 1.32, 1.42, 1.47, 1.39 1.41, 1.36, 1.40, 1.34, 1.42, 1.42, 1.45, 1.35, 1.42, 1.39 1.44, 1.42, 1.39, 1.42, 1.42, 1.30, 1.34, 1.42, 1.37, 1.36 1.37, 1.34, 1.37, 1.37, 1.44, 1.45, 1.32, 1.48, 1.40, 1.45 1.39, 1.46, 1.39, 1.53, 1.36, 1.48, 1.40, 1.39, 1.38, 1.40 1.36, 1.45, 1.50, 1.43, 1.38, 1.43, 1.41, 1.48, 1.39, 1.45 1.37, 1.37, 1.39, 1.45, 1.31, 1.41, 1.44, 1.44, 1.42, 1.47 1.35, 1.36, 1.39, 1.40, 1.38, 1.35, 1.42, 1.43, 1.42, 1.42 1.42, 1.40, 1.41, 1.37, 1.46, 1.36, 1.37, 1.27, 1.37, 1.38 1.42, 1.34, 1.43, 1.42, 1.41, 1.41, 1.44, 1.48, 1.55, 1.37 要求在显著性水平下检验假设
其中F(x)为纤度的分布函数, 为标准正态分布函数. >> x=[
1.36, 1.49, 1.43, 1.41, 1.37, 1.40, 1.32, 1.42, 1.47, 1.39... 1.41, 1.36, 1.40, 1.34, 1.42, 1.42, 1.45, 1.35, 1.42, 1.39... 1.44, 1.42, 1.39, 1.42, 1.42, 1.30, 1.34, 1.42, 1.37, 1.36... 1.37, 1.34, 1.37, 1.37, 1.44, 1.45, 1.32, 1.48, 1.40, 1.45... 1.39, 1.46, 1.39, 1.53, 1.36, 1.48, 1.40, 1.39, 1.38, 1.40... 1.36, 1.45, 1.50, 1.43, 1.38, 1.43, 1.41, 1.48, 1.39, 1.45... 1.37, 1.37, 1.39, 1.45, 1.31, 1.41, 1.44, 1.44, 1.42, 1.47... 1.35, 1.36, 1.39, 1.40, 1.38, 1.35, 1.42, 1.43, 1.42, 1.42... 1.42, 1.40, 1.41, 1.37, 1.46, 1.36, 1.37, 1.27, 1.37, 1.38... 1.42, 1.34, 1.43, 1.42, 1.41, 1.41, 1.44, 1.48, 1.55, 1.37]; %法一 (偏度峰度检验)
>> [h,p,jbstat,cv]= jbtest(x,0.01) h = 0
p = 0.4432 jbstat = 1.6276 cv = 9.2103
%法二 (Lilliefors检验)
>> [h,p,lstat,cv]= lillietest(x,0.01) h = 0
p = 0.0467
lstat = 0.0904 cv = 0.1103
所以不论是偏度峰度检验, 还是Lilliefors检验, 均认为维尼纶纤度服从正态分布.
【例3-6】(例3.23) 抽查用克矽平治疗的矽肺患者10名, 得他们治疗前后血红蛋白的差(g%)如下: 2.7, -1.2, -1.0, 0, 0.7, 2.0, 3.7, -0.6, 0.8, -0.3
试检验治疗前后血红蛋白的差是否服从正态分布( ). >> x=[2.7,-1.2,-1.0,0,0.7,.0,.7,-0.6,.8,-0.3];
>> [h,p,lstat,cv]= lillietest(x,0.05) %本题采用的是Lilliefors检验, 而非书上的W检验. h = 0 p = NaN
lstat = 0.1915 cv = 0.2580
所以认为治疗前后血红蛋白的差服从正态分布.
【例3-7】(例3.25) 以下是两个地区所种小麦的蛋白质含量检验数据: 地区1: 12.6, 13.4, 11.9, 12.8, 13.0
地区2: 13.1, 13.4, 12.8, 13.5, 13.3, 12.7, 12.4 问两地区小麦的蛋白质含量有无显著性差异( )? >> x=[12.6,13.4,11.9,12.8,13.0];
>> y=[13.1,13.4,12.8,13.5,13.3,12.7,12.4];
>> [p,h,stats] = ranksum(x,y,0.05) %秩和检验 p = 0.4066 h = 0
stats = ranksum: 27
>> h=kstest2(x,y,0.05) %斯米尔乐夫检验 h = 0
所以不论是秩和检验还是斯米尔乐夫检验, 均认为两地区小麦的蛋白质含量无显著差异.
【例
(注: 这批数据实际上是来自参数为100的指数分布) ·编写命令文件addexample3_1.m: %几种正态分布检验方法的比较
x=[0.05,4.29,5.22,6.10,6.93,6.96,11.11,11.78,14.18,16.12...
16.53,19.43,26.15,26.59,28.06,35.40,39.64,40.63,51.43,56.79... 57.62,64.99,69.63,78.73,98.07,98.11,98.84,114.16,123.62,124.20...
125.12,133.10,138.15,145.12,155.12,156.41,161.02,203.27,203.30,210.44... 14.51,228.69,234.95,251.246,260.79,272.85,276.26,300.59,301.50,306.54]; normplot(x) pause qqplot(x)
[h_jbtest,p,jbstat,cv]=jbtest(x,0.05) %偏度峰度检验 [h_lillitest,p,lstat,cv]= lillietest(x,0.05)%Lillifors检验 ·运行命令文件addexample3_1.m:
>> addexample3_1 h_jbtest = 0 p = 0.0742 jbstat = 5.2014 cv = 5.9915 h_lillitest = 1 p = 0.0288 lstat = 0.1416 cv = 0.1253
所以由正态概率纸、qq图、偏度峰度检验及Lilliefors检验均认为这批数据不是来自正态分布.
回归分析1
一、线性回归模型
线性模型是指因变量与一个或多个自变量之间的关系可由式(6.1)形式表示:
其中,
y——表示n×1的因变量观测值向量; X——表示n×(p+1)的由自变量决定的设计矩阵; ——表示(p+1)×1的参数向量; ——表示n×1的随机干扰向量, 相互独立且通常具有正态分布. 二、回归分析中研究的主要问题:
1) 确定因变量y与自变量间的定量关系表达式. 这种表达式称为回归方程. 2) 对求得的回归方程的可信度进行检验. 3) 判断自变量对y有无影响.
4) 利用所求得的回归方程进行预测和控制.
Matlab的统计工具箱提供了几种常用的回归分析方法. 三、分析函数
(1) [b,bint,r,rint,stats] = regress(y,X,alpha) 通过求解线性模型 (4.1), 返回y关于X的最小二乘估计. 其中 y——表示n×1的观测向量; X——表示n×(p+1)的设计矩阵;
(4.1)
b——表示的估计;
bint——表示的置信水平为的估计区间, 是(p+1)×2的矩阵; r——表示残差, 是n×1的矩阵;
rint——表示残差的置信水平为的估计区间, 是n×2的矩阵;
stats——包含回归中的 (复) 相关系数的平方R2统计量和F统计量以及概率p. (2) Malab提供了多项式拟合 (polyfit)、多项式评估 (polyval)等函数.
p=polyfit(x,y,n) 在最小二乘意义下, 将向量x,y进行n次多项式拟合, 并返回多项式的系数p. y_fit=polyval(p,x) 返回给定系数p的多项式在x处的预测值y_fit.
★p的第一个元素为的系数, 最后一个元素为常数项, 这与regress命令中的输出b的顺序刚好相反.
【例4-1】(书P
例4.1-4.4) 为研究温度对某个化学过程的生产量的影响, 收集到如下数据(规范化形式):
(1)利用最小二乘法求y对x的回归方程, 作拟合曲线图与观测数据的散点图; (2)在正态分布下利用F检验法检验回归方程效果是否显著 ( ); (3)求出回归系数的置信区间 ( );
(4)取 , 求的预测值与置信水平为的预测区间 ·编写命令文件example4_1.m: %书P138-159例4.1-4.4 alpha=0.05;
x=[-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5]'; n=length(x); X=[ones(n,1),x];
y=[1, 5, 4, 7, 10, 8, 9, 13, 14, 13, 18]'; plot(x,y,'b.') hold on
[b,bint,r,rint,stats] = regress(y,X,alpha); %利用最小二乘法进行线性回归 beta_hat=b %第(1)小题 y_hat=X*beta_hat; plot(x,y_hat,'r')
legend('散点图','回归直线') hold off
F_equation=stats(2:3) %第(2)小题 beta_ci=bint %第(3)小题
x0=3; %第(4)小题,y0的点估计 y0_hat=[1,x0]*beta_hat
Qe=r'*r; %第(4)小题,y0的区间估计 sigma_hat=sqrt(Qe/(n-2)); t_alpha=tinv(1-alpha/2,n-2); Lxx=(x-mean(x))'*(x-mean(x));
delta=sigma_hat*t_alpha*sqrt(1+1/n+(x0-mean(x))^2/Lxx); y0_ci=[y0_hat-delta,y0_hat+delta] ·运行命令文件example4_1.m: >> example4_1
回归分析2
beta_hat = 9.2727 1.4364 F_equation = 96.1798 0.0000 beta_ci = 8.2250 10.3204 1.1050 1.7677 y0_hat = 13.5818
y0_ci = 9.8188 17.3449 所以, 有:
(1)回归方程为 , 回归直线与观测数据散点图见图4-1;
(2)F检验法表明在显著性水平拒绝原假设 , 说明回归方程显著; (3) 的置信区间为[8.225, 10.3204], 的置信区间为[1.1050, 1.7677]; (4)当 , , 的置信水平为的预测区间为[9.8188, 17.3449]. 【例4-2】(例4.5-4.7) 在平炉炼钢中, 由于矿石与炉气的氧化作用, 铁水的总含碳量在不断降低, 一炉钢在冶炼初期总的去碳量y与所加的两种矿石的量及熔化时间有关. 经实测某号平炉的49组数据如下表4-6所列:
设y与之间有线性关系
(1)求y与的回归方程;
(2)检验回归方程和回归系数的显著性. 如有不显著的变量, 请剔除之并求剔除不显著的变量之后的回归方程 ( ); (3)在不剔除不显著变量的前提下, 求回归系数的置信区间 ( );
(4)在不剔除不显著变量的前提下, 若取 , 求的置信水平为95%的预测区间. ·编写命令文件example4_5.m:
%A为观测数据阵, 共49行4列, 此处省略. x=A(:,1:3); [n,p]=size(x); X=[ones(n,1),x]; Y=A(:,4); alpha_1=0.01;
p0=p; %用于记录删除某个自变量后剩余的自变量个数 for i=1:p
[b,bint,r,rint,stats]=regress(Y,X,alpha_1); %第(1)小题--多元线性回归 beta_hat=b %beta的LS估计
F_equation=stats(2:3) %第(2)小题:回归方程的显著性检验 Qe=r'*r; %回归系数的显著性检验 C=inv(X'*X); Cii=diag(C);
F_coeff=(beta_hat(2:p0+1).^2./Cii(2:p0+1))./(Qe/(n-p0-1)) F0=Finv(1-alpha_1,1,n-p0-1)
Hi=F_coeff>F0 %Hi=0,接受原假设H0i=0
k=find(Hi~=0) %剔除不显著变量再次建立回归方程并进行显著性检验 if length(k)
X=[ones(n,1),x]; else
break end end
x=A(:,1:3); %第(3)小题 [n,p]=size(x); X=[ones(n,1),x]; alpha_3=0.05;
[b,bint,r,rint,stats]=regress(Y,X,alpha_3);
beta_ci=bint(2:p+1);
x0=[5,10,50]; %第(4)小题 y0_hat=[1,x0]*b; Qe=r'*r;
sigma=sqrt(Qe/(n-p-1));
t_alpha=tinv(1-alpha_3/2,n-p-1);
y0_ci=[y0_hat-sigma*t_alpha,y0_hat+sigma*t_alpha] ·运行命令文件example4_5.m: >> example4_5
beta_hat = 0.6952 0.1606 0.1076 0.0359 %(1)小题 F_equation = 7.7011 0.0003 %(2)小题 F_coeff = 7.0935 8.2716 11.5683 F0 = 7.2339 Hi = 0 1 1 k = 2 3
beta_hat = 2.5150 0.0233 0.0364 F_equation = 7.0686 0.0021 F_coeff = 1.2050 10.4884 F0 = 7.2200 Hi = 0 1 k = 2
beta_hat = 2.6475 0.0393 F_equation = 12.8760 0.0008 F_coeff = 12.8760 F0 = 7.2068 Hi = 1 k = 1
beta_ci = 0.0392 0.2821 %(3)小题 0.0322 0.1829 0.0147 0.0572
y0_ci = 2.7359 6.0069 %(4)小题 所以, 有:
(1)回归方程为 ;
(2)F检验法表明在显著性水平拒绝原假设 , 说明回归方程显著; 第一次线性回归得回归方程为: , 经F检验知不显著;
剔除进行第二次线性回归, 得回归方程为: , 经F检验知不显著;
剔除进行第三次线性回归, 得回归方程为: , 经F检验知显著, 故最终的线性回归为 .
(3) 的置信区间为[0.0392, 0.2821], 的置信区间为[0.0322, 0.1829], 的置信区间为[0.0147, 0.0572]; (4)当 , 的置信水平为的预测区间为[2.7359, 6.0069].
回归分析3
【例4-3】(书P194例4.10) 某种半成品在生产过程中的废品率y(%)与它所含的某种化学成分x(0.01%)有关, 现将试验所得的16组数据记录如下:
求y对x·编写命令文件example4_10.m: %书P194例4.10
x=[34,36,37,38,39,39,39,40,40,41,42,43,43,45,47,48];
y=[1.30,1.00,0.73,0.90,0.81,0.70,0.60,0.50,0.44,0.56,0.30,0.42,0.35,0.40,0.41,0.60]; plot(x,y,'r*') hold on
p=polyfit(x,y,2)
y_fit=polyval(p,x); y的拟合值 Qe=(y-y_fit)*(y-y_fit)' %残差平方和 plot(x,y_fit,'k') hold off ·运行命令文件example4_10.m: >> example4_10
p = 0.0092 -0.8097 18.2642 Qe = 0.1304
回归分析4
出钢时所用的盛钢水的钢包, 由于钢水对耐火材料的浸蚀, 容积不断增大. 我们希望找出使用次数x与增大的容积y之间的关系. 试验数据列于下表4-7. 首先利用双曲线、倒指数曲线与对数曲线进行非线性回归并利用F检验法进行显著性检验, 然后作出散点图与拟合曲线图, 最后在显著性的模型中比较哪个模型为优. ( )
·编写命令文件example4_8.m: %书P190例4.8 x=2:16;
y=[6.42, 8.20, 9.58, 9.50, 9.70, 10.00, 9.93, 9.99, 10.49, 10.59, 10.60, 10.80, 10.60, 10.90, 10.76]; plot(x,y,'kO') hold on n=length(x);
y1=1./y; %(1)双曲线拟合 x1=1./x;
p1=polyfit(x1,y1,1); ab=[p1(2),p1(1)]
y1_fit=polyval(p1,x1);
U1=(y1_fit-mean(y))*(y1_fit-mean(y))'; %方程的显著性检验 Qe1=(y1-y1_fit)*(y1-y1_fit)'; F1=U1/Qe1/(n-2); F0=Finv(1-0.01,1,n-2); if F1>F0 h=1 else h=0 end
y_fit=1./y1_fit;
Qe=(y-y_fit)*(y-y_fit)' %残差平方和 plot(x,y_fit,'r')
y2=log(y); %(2)倒指数曲线拟合 x2=1./x;
p2=polyfit(x2,y2,1); ab=[exp(p2(2)),p2(1)] y2_fit=polyval(p2,x2);
U2=(y2_fit-mean(y))*(y2_fit-mean(y))'; %方程的显著性检验 Qe2=(y2-y2_fit)*(y2-y2_fit)'; F2=U2/Qe2/(n-2); F0=Finv(1-0.01,1,n-2); if F2>F0 h=1 else h=0 end
y_fit=exp(y2_fit);
Qe=(y-y_fit)*(y-y_fit)' %残差平方和 plot(x,y_fit,'b-.')
y3=y; %(3)对数曲线拟合 x3=log(x);
p3=polyfit(x3,y3,1); ab=[p3(2),p3(1)]
y3_fit=polyval(p3,x3);
U3=(y3_fit-mean(y))*(y3_fit-mean(y))'; %方程的显著性检验 Qe3=(y3-y3_fit)*(y3-y3_fit)'; F3=U3/Qe3/(n-2); F0=Finv(1-0.01,1,n-2); if F3>F0 h=1 else h=0 end
y_fit=y3_fit;
Qe=(y-y_fit)*(y-y_fit)' %残差平方和 plot(x,y_fit,'m:')
hold off
legend('散点图','1/y=a+b/x拟合','y=a*e^(b/x)拟合','y=a+b*lnx拟合') ·运行命令文件example4_8.m: >> example4_8
ab = 0.0823 0.1312 %第(1)模型的结果 h = 1
Qe = 1.4396
ab = 11.6789 -1.1107 %第(2)模型的结果 h = 1
Qe = 0.8913
ab = 6.2389 1.7761 %第(3)模型的结果 h = 0
Qe = 2.5592
所以, 有:
(1)双曲线回归方程为 , 此方程显著 ; 倒指数曲线回归方程为 , 此方程也显著 ; 对数曲线回归方程为 , 此方程不显著.
(2)在显著型模型中, 倒指数模型由于残差平方和最小, 所以比双曲线模型更优.
方差分析1
单因素方差与双因素方差分析均假定X中的数据具有以下特征: 独立、正态、方差齐性.
表5.1方差分析部分分析函数
(1) [p,table,stats] = anova1(X) 进行均衡的单因素方差分析 (如果数据缺失, 用NaN补齐). 检验数据X (m×n矩阵)中各列的均值是否相等, X的每一列 (一个水平) 表示独立样本, 每个样本包含m个相互独立的观察值. 默认返回两幅图表, 第一幅为标准anova表, 第二幅为X各列数据的盒形 (box) 图. 如果盒形图的中心线差别很大, 则对应的F值很大, 相应的概率p-值就小.
p——零假设 ( : 各样本具有相同的均值)的概率值. 如果p-值接近零则应对零假设置疑, 表明至少有一个样本的均值与其它样本的均值显著不同;
table——以单元数组的形式返回anova表;
stats——利用stats可接下来进行多重比较检验. 用户可以将stats结构作为输入利用multcompare函数进行这种检验. (2) [p,table,stats] = anova2(X,reps,'displayopt') 进行均衡的双因素方差分析 (对于不均衡设计, 使用函数anovan). 检验数据X中各列或各行的均值是否相等. 不同列中的数据表示因素A引起的变化情况, 不同行中的数据表示因素B引起的变化. 如果对于因素A和B的每一种水平组合都有超过一个的观察值 (这种情况又称重复试验双因素方差分析), 则输入reps表示每个单元 (cell) (对应一个水平组合) 中观察值的个数, 它必须为常数.
X数据的结构格式为一矩阵. 比如下面这个矩阵, 其中列因素A有两个水平, 行因素B有三个水平, 每个水平对有两个观察值 (reps=2). 下标分别表示行、列和每个水平对中的观察值.
当reps=1 (缺省值) 时, anova2返回的向量p中包含如下两个零假设的概率值:
: 因素A各样本 (即X中所有的列样本) 均来自相同的总体, 或指因素A的不同水平对试验结果具有相同的影响); : 因素B各样本 (即X中所有的行样本) 均来自相同的总体;
当reps>1时, anova2返回的向量p还包含下面这个零假设的概率值: : 因素A与B之间无交互作用.
(3) c = multcompare(stats,alpha,'displayopt') 均值或其它估计量的多重比较. 利用上面命令输出中的stats结构中的信息进行多重比较检验, 返回成对比较的结果矩阵c. c包含5列矩阵. 矩阵的每一行代表一次检验, 每一对均值比较对应一行. 行中的第一、二个数据表示比较组的编号, 第四个数据表示正在被比较的均值差的估计值, 第三、五个数据表示正在被比较的均值差的置信区间. 如果置信区间中包含0, 则说明在alpha的显著水平上差异不是显著的. 此命令结果也显示一个表示检验的交互式图形. 图中每组的均值用一个符号和符号周围的区间表示. 如果两个均值的区间不交叠, 说明它们显著不同; 如果两个区间交叠, 则说明它们不是显著不同的. 可以用鼠标选中任何一组, 图中其它任何与之显著不同的组将会高亮度显示.
【例5-1】(书P209例5.3) 对六种不同的农药在相同的条件下分别进行杀虫试验, 试验结果见下表. 问杀虫率是否因
·编写命令文件example5_3.m:
%书P209例5.3
X=[87,90,56,55,92,75; 85,88,62,48,99,72; 80,87,NaN,NaN,95,81;
NaN,94,NaN,NaN,91,NaN]; n=length(find(X>0)); r=6;
[p,table,stats] = anova1(X); F0=Finv(1-0.01,r-1,n-r) pause
c = multcompare(stats,0.01)
21 / 21