统计的回归分析 --matlab实现
前面我们讲过曲线拟合问题。曲线拟合问题的特点是,根据得到的若干有关变量的 一组数据,寻找因变量与(一个或几个)自变量之间的一个函数,使这个函数对那组数 据拟合得好。通常,函数的形式可以由经验、先验知识或对数据的直观观察决定,要 作的工作是由数据用小二乘法计算函数中的待定系数。从计算的角度看,问题似乎已 经完全解决了,还有进一步研究的必要吗?
目录
1.1 样本空间 、样本均值、样本协方差矩阵、样本相关系数矩阵
(1)数据的中心化处理 (2)数据的无量纲化处理 (3)标准化处理
2.1 小二乘估计方法 、正规方程组 2.2 回归方程的估计
2.4 拟合效果分析 2.4.1 残差的样本方差 2.4.2 判定系数(拟合优度)
3.4 回归模型的假设检验 3.5 回归系数的假设检验和区间估计
4.2.1 一元多项式回归 polyfit 4.2.2 多元二项式回归 rstool
从数理统计的观点看,这里涉及的都是随机变量,我们根据一个样本计算出的那些 系数,只是它们的一个(点)估计,应该对它们作区间估计或假设检验,如果置信区间 太大,甚至包含了零点,那么系数的估计值是没有多大意义的。另外也可以用方差分析方法对模型的误差进行分析,对拟合的优劣给出评价。简单地说,回归分析就是对拟合 问题作的统计分析。
具体地说,回归分析在一组数据的基础上研究这样几个问题:
(i)建立因变量 y 与自变量 之间的回归模型(经验公式);
(ii)对回归模型的可信度进行检验;
(iii)判断每个自变量 对 y 的影响是否显著;
(iv)诊断回归模型是否适合这组数据;
(v)利用回归模型对 y 进行预测或控制。
1 数据表的基础知识
1.1 样本空间 、样本均值、样本协方差矩阵、样本相关系数矩阵
1.2 数据的标准化处理
(1)数据的中心化处理
该变换可以使样本的均值变为 0,而这样的变换既不改变样本点间的相互位置,也 不改变变量间的相关性。但变换后,却常常有许多技术上的便利。
(2)数据的无量纲化处理
在实际问题中,不同变量的测量单位往往是不一样的。为了消除变量的量纲效应, 使每个变量都具有同等的表现力,数据分析中常用的消量纲的方法,是对不同的变量进行所谓的压缩处理,即使每个变量的方差均变成 1,即
(3)标准化处理
所谓对数据的标准化处理,是指对数据同时进行中心化-压缩处理,即
2 一元线性回归模型
2.1 小二乘估计方法 、正规方程组
2.2 回归方程的估计
2.3.1
、
的性质
2.3.2 其它性质
用小二乘法拟合的回归方程还有一些值得注意的性质:
2.4 拟合效果分析
当根据一组观测数据得到小二乘拟合方程后,必须考察一下,是否真的能由所得的模型( )来较好地拟合观测值
?用
能否较好地反映 (或者说解释)
值的取值变化?回归方程的质量如何?误差多大?对这些,都必须 予以正确的评估和分析。
2.4.1 残差的样本方差
一个好的拟合方程,其残差总和应越小越好。残差越小,拟合值与观测值越接近, 各观测点在拟合直线周围聚集的紧密程度越高,也就是说,拟合方程 解释 y 的能力越强。 另外,当
越小时,还说明残差值
的变异程度越小。由于残差的样本均值为零, 所以,其离散范围越小,拟合的模型就越为精确。
2.4.2 判定系数(拟合优度)
对应于不同的 值,观测值
的取值是不同的。建立一元线性回归模型的目的, 就是试图以x 的线性函数(
)来解释 y 的变异。那么,回归模型
究竟能以多大的精度来解释 y 的变异呢?又有多大部分是无法用这个回归方程来解释 的呢?
从上式可以看出,y 的变异是由两方面的原因引起的;一是由于x 的取值不同,而 给 y 带来的系统性变异;另一个是由除 x以外的其它因素的影响。
注意到对于一个确定的样本(一组实现的观测值),SST 是一个定值。所以,可解释变异SSR 越大,则必然有残差SSE 越小。这个分解式可同时从两个方面说明拟合方 程的优良程度:
(1)SSR 越大,用回归方程来解释 变异的部分越大,回归方程对原数据解释得 越好;
(2)SSE 越小,观测值 绕回归直线越紧密,回归方程对原数据的拟合效果越好。
因此,可以定义一个测量标准来说明回归方程对原始数据的拟合程度,这就是所谓的判定系数,有些文献上也称之为拟合优度。 判定系数是指可解释的变异占总变异的百分比,用 表示,有
从判定系数的定义看, 有以下简单性质:
(1) 0 ≤ ≤ 1 ;
(2)当 =1 时,有 SSR = SST,也就是说,此时原数据的总变异完全可以由拟 合值的变异来解释,并且残差为零( SSE =0 ),即拟合点与原数据完全吻合;
(3)当 =0 时,回归方程完全不能解释原数据的总变异,y 的变异完全由与 x无关的因素引起,这时 SSE = SST 。
测定系数时一个很有趣的指标:一方面它可以从数据变异的角度指出可解释的变异占总变异的百分比,从而说明回归直线拟合的优良程度;另一方面,它还可以从相关性 的角度,说明原因变量 y 与拟合变量 的相关程度,从这个角度看,拟合变量
与原 变量 y 的相关度越大,拟合直线的优良度就越高。
2.5 显著性检验
2.5.1 回归模型的线性关系检验
在拟合回归方程之前,我们曾假设数据总体是符合线性正态误差模型的,也就是说, y 与 x之间的关系是线性关系,即
然而,这种假设是否真实,还需进行检验。 对于一个实际观测的样本,虽然可以用判定系数 说明 y 与
的相关程度,但是, 样本测度指标具有一定的随机因素,还不足以肯定 y 与 x的线性关系。
假设 y 与 x之间存在线性关系,则总体模型为
如果 ,则称这个模型为全模型。 用小二乘法拟合全模型,并求出误差平方和为
也即,全模型的误差总是小于(或等于)选模型的误差的。其原因是在全模型中有 较多的参数,可以更好地拟合数据。
假若在某个实际问题中,全模型的误差并不比选模型的误差小很多的话,这说明 假设成立,即
近似于零。因此,差额 (SST − SSE )很少时,表明
成立。若这个差额很大,说明增加了 x的线性项后,拟合方程的误差大幅度减少,则应否定
, 认为总体参数
显著不为零。
就是一个恰当的回归模型,事实上,当 假设被拒绝后,只能说明 y 与x之间存在显 著的线性关系,但很有可能在模型中还包括更多的回归变量,而不仅仅是一个回归变量 x 。
一般地,回归方程的假设检验包括两个方面:
一个是对模型的检验,即检验自变量与因变量之间的关系能否用一个线性模型来表示,这是由F 检验来完成的;
另一个检验是关于回归参数的检验,即当模型检验通过后,还要具体检验每一个自变量对因变量 的影响程度是否显著。这就是下面要讨论的t检验。在一元线性分析中,由于自变量的个数只有一个,这两种检验是统一的,它们的效果完全是等价的。 但是,在多元线性回归分析中,这两个建议的意义是不同的。从逻辑上说,一般常在F 检验通过后,再进一步进行t检验。
2.5.2 回归系数的显著性检验
回归参数的检验是考察每一个自变量对因变量的影响是否显著。换句话说,就是要检验每一个总体参数是否显著不为零。 首先看对 的检验。
代表
变化一个单位对
的影响程度。对
的检验 就是要看这种影响程度与零是否有显著差异。
3 多元线性回归模型
3.1 参数估计
3.3 统计分析
不加证明地给出以下结果:
3.4 回归模型的假设检验
3.5 回归系数的假设检验和区间估计
3.6 利用回归模型进行预测
4 Matlab 中的回归分析
4.1 多元线性回归
Matlab 统计工具箱用命令 regress 实现多元线性回归,用的方法是小二乘法,用 法是:
b=regress(Y,X)
[b,bint,r,rint,stats]=regress(Y,X,alpha)
rcoplot(r,rint)
残差及其置信区间可以用 rcoplot(r,rint)画图。
例 1 合金的强度 y 与其中的碳含量 x有比较密切的关系,今从生产中收集了一批 数据如下表 1。
试先拟合一个函数 y (x),再用回归分析对它进行检验。
解 先画出散点图:
x=0.1:0.01:0.18;
y=[42,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0];
plot(x,y,'+')
可知 y 与 x大致上为线性关系。 设回归模型为
用 regress 和 rcoplot 编程如下:
观察命令 rcoplot(r,rint)所画的残差分布,除第 8 个数据外其余残差的置信区间均包 含零点,第 8 个点应视为异常点,将其剔除后重新计算,可得
应该用修改后的这个结果。
编写如下程序:
x1=[120 140 190 130 155 175 125 145 180 150]';
x2=[100 110 90 150 210 150 250 270 300 250]';
y=[102 100 120 77 46 93 26 69 65 85]';
x=[ones(10,1),x1,x2];
[b,bint,r,rint,stats]=regress(y,x);
b,bint,stats
4.2 多项式回归
如果从数据的散点图上发现 y 与 x呈较明显的二次(或高次)函数关系,或者用线 性模型(20)的效果不太好,就可以选用多项式回归。
4.2.1 一元多项式回归 polyfit
一元多项式回归可用命令 polyfit 实现。
例 3 将 17 至 29 岁的运动员每两岁一组分为 7 组,每组两人测量其旋转定向能力, 以考察年龄对这种运动能力的影响。现得到一组数据如表 3。
试建立二者之间的关系。 解 数据的散点图明显地呈现两端低中间高的形状,所以应拟合一条二次曲线。 选用二次模型
编写如下程序:
x0=17:2:29;x0=[x0,x0];
y0=[20.48 25.13 26.15 30.0 26.1 20.3 19.35...
24.35 28.11 26.3 31.4 26.92 25.7 21.3];
[p,s]=polyfit(x0,y0,2);
p
上面的s是一个数据结构,用于计算函数值,如
[y,delta]=polyconf(p,x0,s);
y
得到 y 的拟合值,及预测值 y 的置信区间半径delta。
4.2.2 多元二项式回归 rstool
统计工具箱提供了一个作多元二项式回归的命令rstool,它也产生一个交互式画面, 并输出有关信息,用法是
rstool(x,y,model,alpha)
编程如下:
x1=[120 140 190 130 155 175 125 145 180 150]';
x2=[100 110 90 150 210 150 250 270 300 250]';
y=[102 100 120 77 46 93 26 69 65 85]';
x=[x1 x2];
rstool(x,y,'purequadratic')
图的左下方有两个下拉式菜单,一个菜单Export用以向Matlab工作区传送数据,包 括beta(回归系数),rmse(剩余标准差),residuals(残差)。
模型(41)的回归系数和剩余 标准差为
5 非线性回归和逐步回归
本节介绍怎样用Matlab统计工具箱实现非线性回归和逐步回归。
5.1 非线性回归
非线性回归是指因变量 y 对回归系数 (而不是自变量)是非线性的。 Matlab统计工具箱中的nlinfit,nlparci,nlpredci,nlintool,不仅给出拟合的回归系数, 而且可以给出它的置信区间,及预测值和置信区间等。下面通过例题说明这些命令的用 法。
例4 在研究化学动力学反应过程中,建立了一个反应速度和反应物含量的数学模 型,形式为
解 首先,以回归系数和自变量为输入变量,将要拟合的模型写成函数文件 huaxue.m:
function yhat=huaxue(beta,x);
yhat=(beta(4)*x(:,2)-x(:,3)/beta(5))./(1+beta(1)*x(:,1)+...
beta(2)*x(:,2)+beta(3)*x(:,3));
然后,用nlinfit计算回归系数,用nlparci计算回归系数的置信区间,用nlpredci 计算预测值及其置信区间,编程如下:
clc,clear
x0=[ 1 8.55 470 300 10
2 3.79 285 80 10
3 4.82 470 300 120
4 0.02 470 80 120
5 2.75 470 80 10
6 14.39 100 190 10
7 2.54 100 80 65
8 4.35 470 190 65
9 13.00 100 300 54
10 8.50 100 300 120
11 0.05 100 80 120
12 11.32 285 300 10
13 3.13 285 190 120];
x=x0(:,3:5);
y=x0(:,2);
beta=[0.1,0.05,0.02,1,2]'; %回归系数的初值,任意取的
[betahat,r,j]=nlinfit(x,y,@huaxue,beta); %r,j是下面命令用的信息
betaci=nlparci(betahat,r,'jacobian',j);
betaa=[betahat,betaci] %回归系数及其置信区间
[yhat,delta]=nlpredci(@huaxue,x,betahat,r,'jacobian',j) %y的预测值及其置信区间的半径,置信区间为yhat±delta。
用nlintool得到一个交互式画面,左下方的Export可向工作区传送数据,如剩余标准差等。使用命令
nlintool(x,y,'huaxue',beta)
可看到画面,并传出剩余标准差rmse= 0.1933。
4.2 逐步回归
实际问题中影响因变量的因素可能很多,我们希望从中挑选出影响显著的自变量来 建立回归模型,这就涉及到变量选择的问题,逐步回归是一种从众多变量中有效地选择 重要变量的方法。以下只讨论线性回归模型(1)式的情况。 变量选择的标准,简单地说就是所有对因变量影响显著的变量都应选入模型,而影 响不显著的变量都不应选入模型,从便于应用的角度应使模型中变量个数尽可能少。
逐步回归是实现变量选择的一种方法,基本思路为,先确定一初始子集,然后每次 从子集外影响显著的变量中引入一个对 y 影响大的,再对原来子集中的变量进行检 验,从变得不显著的变量中剔除一个影响小的,直到不能引入和剔除为止。使用逐步 回归有两点值得注意:
一是要适当地选定引入变量的显著性水平 和剔除变量的显著 性水平
,显然,
越大,引入的变量越多;
越大,剔除的变量越少。
二是由 于各个变量之间的相关性,一个新的变量引入后,会使原来认为显著的某个变量变得不 显著,从而被剔除,所以在初选择变量时应尽量选择相互独立性强的那些。
在Matlab统计工具箱中用作逐步回归的是命令stepwise,它提供了一个交互式画面,通过这个工具你可以自由地选择变量,进行统计分析,其通常用法是:
stepwise(x,y,inmodel,alpha)
其中x是自变量数据,y是因变量数据,分别为 n× m 和n ×1 矩阵,inmodel是矩阵x的 列数的指标,给出初始模型中包括的子集(缺省时设定为空),alpha为显著性水平。
Stepwise Regression 窗口,显示回归系数及其置信区间,和其它一些统计量的信 息。绿色表明在模型中的变量,红色表明从模型中移去的变量。在这个窗口中有Export 按钮,点击Export产生一个菜单,表明了要传送给Matlab工作区的参数,它们给出了统 计计算的一些结果。
下面通过一个例子说明stepwise的用法。
例5 水泥凝固时放出的热量 y 与水泥中4种化学成分 有关,今测得一 组数据如表5,试用逐步回归来确定一个线性模型
编写程序如下:
clc,clear
x0=[1 7 26 6 60 78.5
2 1 29 15 52 74.3
3 11 56 8 20 104.3
4 11 31 8 47 87.6
5 7 52 6 33 95.9
6 11 55 9 22 109.2
7 3 71 17 6 102.7
8 1 31 22 44 72.5
9 2 54 18 22 93.1
10 21 47 4 26 115.9
11 1 40 23 34 83.8
12 11 66 9 12 113.3
13 10 68 8 12 109.4];
x=x0(:,2:5);
y=x0(:,6);
stepwise(x,y)
得到图3所示的图形界面。
习 题
1. 某人记录了21天每天使用空调器的时间和使用烘干器的次数,并监视电表以计 算出每天的耗电量,数据见表6,试研究耗电量(KWH)与空调器使用的小时数(AC)和烘 干器使用次数(DRYER)之间的关系,建立并检验回归模型,诊断是否有异常点。
2. 在一丘陵地带测量高程,x和 y 方向每隔100米测一个点,得高程如下表,试拟 合一曲面,确定合适的模型,并由此找出高点和该点的高程。
3. 一矿脉有13个相邻样本点,人为地设定一原点,现测得各样本点对原点的距离 x,与该样本点处某种金属含量 y 的一组数据如下,画出散点图观测二者的关系,试建 立合适的回归模型,如二次曲线、双曲线、对数曲线等。