第四讲数据分析方法.docx
《第四讲数据分析方法.docx》由会员分享,可在线阅读,更多相关《第四讲数据分析方法.docx(30页珍藏版)》请在冰豆网上搜索。
第四讲数据分析方法
第四讲数据分析方法
第一节、数据拟合
问题:
给定一批数据点(输入变量与输出变量的数据),需确定满足特定要求的曲线或
曲面。
如果输入变量和输出变量都只有一个,则属于一元函数的拟合和插值;而若输入变量
有多个,则为多元函数的拟合和插值(有点回归分析的意思)
解决方案:
(1)若要求所求曲线(面)通过所给所有数据点,就是插值问题;
(2)若不要求曲线(面)通过所有数据点,而是要求它反映对象整体的变化趋势,这就
是数据拟合,又称曲线拟合或曲面拟合。
注意:
插值和拟合都是要根据一组数据构造一个函数作为近似,由于近似的要求不同,
二者的数学方法上是完全不同的。
而面对一个实际问题,究竟应该用插值还是拟合,有时
容易确定,有时则并不明显。
例1:
下面数据是某次实验所得,希望得到X和f之间的关系?
x
f
1
2
4
7
9
12
13
15
17
1.5
3.9
6.6
11.7
15.6
18.8
19.6
20.6
21.1
曲线拟合问题最常用的解法——最小二乘法的基本思路
第一步:
确定拟合的函数类型y=f(x;a1,a2,",am),其中a1,a2,",am为待定系数。
(函数类型的确定可以根据内在的规律确定,如果无现成的规则,则可以通过散点图,联系
曲线的形状进行分析)
第二步:
确定a1,a2,",am的最小二乘准则:
要求n个已知点(xi,yi)与曲线y=f(x)
n
∑
的距离di的平方和(yi−fx2最小。
())
i
i=1
用MATLAB作拟合
1.多项式拟合。
作多项式y=a0x
m
+a1xm−1+"+am拟合,可利用
a=polyfit(x,y,m)—其中x,y为给出的数据,m为多项式的次数。
多项式在x处的值y可用以下命令计算:
y=polyval(a,x)
2.用MATLAB作非线性最小二乘拟合
Matlab的提供了两个求非线性最小二乘拟合的函数:
lsqcurvefit和lsqnonlin。
两个命令
都要先建立M-文件fun.m,在其中定义函数f(x)。
(1)x=lsqcurvefit(‘fun’,x0,xdata,ydata);
(2)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);
(3)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options,’grad’);
(4)[x,options]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);
(5)[x,options,funval]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);
1
(6)[x,options,funval,Jacob]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);
(7)x=lsqnonlin(‘fun’,x0);
(8)x=lsqnonlin(‘fun’,x0,options);
(9)x=lsqnonlin(‘fun’,x0,options,‘grad’);
(10)[x,options]=lsqnonlin(‘fun’,x0,…);
(11)[x,options,funval]=lsqnonlin(‘fun’,x0,…);
在使用lsqcurvefit与lsqnonlin命令时,共同的问题是要先知道函数的类型,而拟合其实是决
定函数中的待定系数。
第二节插值
插值的基本问题:
给出n个数对,(Pi,f(Pi)),i=1,2,",n,求点P处对应的函数值f(P)。
一、一维插值
已知n+1个节点(xi,yi),i=0,1,",n,求任意点*处的函数值*。
常用的插值方法
x
y
有拉格朗日多项式插值、牛顿插值、分段线性插值、Hermite插值和三次样条插值。
9分段线性插值:
将各数据点用折线连接起来
9多项式插值:
求一个多项式通过所有数据点,可以假设出多项式的系数,最后通过
求解方程得到每个系数(拉格朗日插值,用n次多项式描述n+1个点)
9样条插值:
分段多项式的光滑连接(三次样条插值)
9牛顿插值:
利用节点之间的各阶差商和差分构造多项式
9
Hermite插值:
对插值函数,不仅要求它在节点处与函数同值,而且要求它与函数有
相同的一阶、二阶甚至更高阶的导数值
(1)MATLAB命令:
y=interp1(x0,y0,x,'method')
method指定插值的方法,默认为线性插值。
其值可为:
'nearest'最近项插值
'linear'
'spline'
'cubic'
线性插值
立方样条插值
立方插值。
所有的插值方法要求x0是单调的。
当x0为等距时可以用快速插值法,使用快速插值法的格式为'*nearest'、'*linear'、'*spline'、
'*cubic'。
(2)三次样条插值在Matlab中的实现
在Matlab中数据点称之为断点。
如果三次样条插值没有边界条件,最常用的方法,就
是采用非扭结(not-a-knot)条件。
这个条件强迫第1个和第2个三次多项式的三阶导数相
等。
对最后一个和倒数第2个三次多项式也做同样地处理。
Matlab中三次样条插值也有现成的函数:
y=interp1(x0,y0,x,'spline');
y=spline(x0,y0,x);
pp=csape(x0,y0,conds),
pp=csape(x0,y0,conds,valconds),y=ppval(pp,x)。
其中x0,y0是已知数据点,x是插值点,y是插值点的函数值。
对于三次样条插值,我们提倡使用函数csape,csape的返回值是pp形式,要求插值点
2
的函数值,必须调用函数ppval。
pp=csape(x0,y0):
使用默认的边界条件,即Lagrange边界条件。
pp=csape(x0,y0,conds,valconds)中的conds指定插值的边界条件,其值可为:
'complete'
边界为一阶导数,一阶导数的值在valconds参数中给出,若忽略valconds
参数,则按缺省情况处理。
'not-a-knot'非扭结条件
'periodic'
'second'
周期条件
边界为二阶导数,二阶导数的值在valconds参数中给出,若忽略valconds
参数,二阶导数的缺省值为[0,0]。
'variational'设置边界的二阶导数值为[0,0]。
对于一些特殊的边界条件,可以通过conds的一个1×2矩阵来表示,conds元素的取值
为0,1,2。
conds(i)=j的含义是给定端点i的j阶导数,即conds的第一个元素表示左边界的条件,
第二个元素表示右边界的条件,conds=[2,1]表示左边界是二阶导数,右边界是一阶导数,对
应的值由valconds给出。
例2机床加工
待加工零件的外形根据工艺要求由一组数据(x,y)给出(在平面情况下),用程控铣床
加工时每一刀只能沿x方向和y方向走非常小的一步,这就需要从已知数据得到加工所要求
的步长很小的(x,y)坐标。
表中给出的x,y数据位于机翼断面的下轮廓线上,假设需要得到x坐标每改变0.1时的
y坐标。
试完成加工所需数据,画出曲线,并求出x=0处的曲线斜率和13≤x≤15范围
内y的最小值。
x
y
0
3
5
7
9
11
2.0
12
1.8
13
1.2
14
1.0
15
1.6
0
1.2
1.7
2.0
2.1
要求用分段线性和三次样条计算。
解编写以下程序:
x0=[035791112131415];
y0=[01.21.72.02.12.01.81.21.01.6];
x=0:
0.1:
15;
y2=interp1(x0,y0,x);
y3=interp1(x0,y0,x,'spline');
pp1=csape(x0,y0);
y4=ppval(pp1,x);
pp2=csape(x0,y0,'second');
y5=ppval(pp2,x);
[x',y1',y2',y3',y4',y5']
subplot(2,2,1)
plot(x0,y0,'+',x,y2)
title('Piecewiselinear')
subplot(2,2,2)
plot(x0,y0,'+',x,y3)
3
title('Spline1')
subplot(2,2,3)
plot(x0,y0,'+',x,y4)
title('Spline2')
dx=diff(x);
dy=diff(y3);
dy_dx=dy./dx;
dy_dx0=dy_dx
(1)
ytemp=y3(131:
151);
ymin=min(ytemp);
index=find(y3==ymin);
xmin=x(index);
[xmin,ymin]
二、二维插值:
前面讲述的都是一维插值,即节点为一维变量,插值函数是一元函数(曲线)。
若节点
是二维的,插值函数就是二元函数,即曲面。
如在某区域测量了若干点(节点)的高程(节
点值),为了画出较精确的等高线图,就要先插入更多的点(插值点),计算这些点的高程(插
值)。
(1)
(2)
网格节点:
已知m×n个节点(xi,yj,zij),i=1,2,",m;j=1,2,",n,构造
一个二元函数z=f(x,y)通过全部已知节点,即zij=f(xi,yj),再利用
z=f(x,y)插值。
散乱节点:
已知n个节点(xi,yi,zi),i=1,2,",n,构造一个二元函数
z=f(x,y)通过全部已知节点,即zi=f(xi,yi),再利用z=f(x,y)插值。
二维插值方法:
(1)最邻近插值:
二维或高维情形的最邻近插值,与被插值点最邻近的节点的函数值
即为所求。
(2)分片线性插值:
将四个插值点(矩形的四个顶点)处的函数值依次简记为:
f(xi,yj)=f1,f(xi+1,yj)=f2,f(xi+1,yj+1)=f3,f(xi,yj+1)=f4
分两片的函数表达式如下:
第一片(下三角形区域):
f(x,y)=f1+(f2−f1)(x−xi)+(f3−f2)(y−yj)
第二片(上三角形区域):
f(x,y)=f1+(f4−f1)(y−yj)+(f3−f4)(x−xi)
(3)双线性插值:
双线性插值是一片一片的空间二次曲面构成。
双线性插值函数的形
式如下:
f(x,y)=(ax+b)(cy+d)。
其中有四个待定系数,利用该函数在矩形的四个顶点
(插值节点)的函数值,得到四个代数方程,正好确定四个系数。
用MATLAB作网格节点数据的插值:
4
z=interp2(x0,y0,z0,x,y,’method’)—method可以选择‘nearest’:
最邻近插值;‘linear’:
双线性插值;‘cubic’:
双三次插值;缺省时,双线性插值。
如果是三次样条插值,可以使用命令
pp=csape({x0,y0},z0,conds,valconds),z=fnval(pp,{x,y})
其中x0,y0分别为m维和n维向量,z0为m×n维矩阵,z为矩阵,它的行数为x的维数,
列数为y的维数,表示得到的插值,具体使用方法同一维插值。
例3在一丘陵地带测量高程,x和y方向每隔100米测一个点,得高程如下表,试插值
一曲面,确定合适的模型,并由此找出最高点和该点的高程。
x
y
100
200
300
400
100
200
697
712
674
626
300
624
630
598
552
400500
478450
478420
412400
334310
636
698
680
662
解编写程序如下:
x=100:
100:
500;
y=100:
100:
400;
z=[636697624478450;
698712630478420;
680674598412400;
662626552334310];
pp=csape({x,y},z')
xi=100:
10:
500;yi=100:
10:
400
cz1=fnval(pp,{xi,yi})
cz2=interp2(x,y,z,xi,yi','spline')
[i,j]=find(cz1==max(max(cz1)))
x=xi(i),y=yi(j),zmax=cz1(i,j)
用MATLAB作散点数据的插值计算
cz=griddata(x,y,z,cx,cy,‘method’)—method可以选择‘nearest’:
最邻近插值;‘linear’:
双线性插值;‘cubic’:
双三次插值;'v4':
Matlab提供的插值方法;缺省时,双线性插值
例4在某海域测得一些点(x,y)处的水深z由下表给出,在矩形区域(75,200)×(-50,150)
内画出海底曲面的图形。
x129
y7.5
z4
140
141.5
8
103.5
23
6
88
147
8
185.5
22.5
6
195
137.5
8
105
157.5
–6.5
9
107.5
-81
9
77
81
56.5
8
162
–66.5
9
162
117.5
-33.5
9
85.5
3
84
8
8
4
解编写程序如下:
x=[129,140,103.5,88,185.5,195,105,157.5,107.5,77,81,162,162,117.5];
y=[7.5,141.5,23,147,22.5,137.5,85.5,-6.5,-81,3,56.5,-66.5,84,-33.5];
z=[4,8,6,8,6,8,8,9,9,8,8,9,4,9];
xi=75:
1:
200;
5
yi=-50:
1:
150;
zi=griddata(x,y,z,xi,yi','cubic')
subplot(1,2,1)
plot(x,y,'*')
subplot(1,2,2)
mesh(xi,yi,zi)
第三节、回归分析方法
“回归”问题最早来源于生物界,英国生物学家兼统计学家高尔顿(Galton,1822-1911)
发现同一种族中儿子的平均高度介于其父亲的高度与种族平均高度之间。
儿子的身高有返归
于种族平均身高的趋势,即回归于种族的平均身高。
回归分析是指对具有相关关系的现象,根据其关系形态,选择一个合适的数学模型,用
来近似地表示变量间的平均变化关系的一种统计方法。
回归分析的分类:
按照回归模型中变量个数分(一元回归,多元回归);按照回归曲线
的形态分(线性回归,非线性回归);按照是否要求总体分布类型已知分(参数回归,非参
数回归)
一元线性回归特点
¾两个变量中,一个是自变量,一个是因变量
¾回归方程不是抽象的数学模型,而随机方程,可以进行实证
¾因果关系不明显时,应同时作两个回归方程
¾回归系数具有较强的经济含义
¾作为回归模型的因变量是随机变量,而自变量是确定性变量,即可控变量
1.一元线性回归模型
一般地,称由y=β0+β1x+ε确定的模型为一元线性回归模型,记为
⎧y=β0+β1x+ε
⎨
Eε=0,Dε=σ
2
⎩
其中固定的未知参数β0,β1称为回归系数,自变量x也称为回归变量,Y=β0+β1x称
为y对x的回归直线方程。
一元线性回归分析的主要任务
1.用试验值(样本值)对β0,β1和σ作点估计;
2.对回归系数β0,β1作假设检验;
3.在x=x0处对y作预测,对y作区间估计。
回归系数的最小二乘估计
有n组独立观测值,(xi,yi),i=1,2,",n。
⎧yi=β0+βx1+εi,i=1,2,...,n
设⎩⎨Eεi=0,Dεi=σ2且ε1ε2,...,εn相互独立
6
n
n
∑∑
ε
=(yi−β0−β1xi)2。
最小二乘法就是选择
2
i
记Q=Q(β0,β1)=
β0和β1的估
i=1
i=1
计β0,β1使得Q(βˆ0,βˆ1)=minQ(β0,β1)。
计算得到
ˆ
ˆ
β,β
0
1
⎧ˆ
β0=y−βˆ1x
⎪
⎨
xy−xy
ˆ
β=
⎪
⎩
1
x−x
22
1
n
1
n
1
n
1
n
∑
∑
∑
∑
其中x=
xi,y=n
yi,x2=
xi2,xy=
xiyi。
(经验)回归方程为:
n
n
n
i=1
i=1
i=1
i=1
yˆ=βˆ0+βˆ1x=y+βˆ1(x−x)。
记Q=Q(βˆ0,βˆ)
=(yi−0−1xi)2=(yi−y2,称e为残差平方和或剩
∑
n
β
ˆ
β
ˆ
∑
n
Q
ˆ)
i
e
1
i=1
i=1
σ
余平方和。
2的无偏估计为
σˆe
2
=Qe(n−2)。
回归方程的显著性检验
对回归方程Y=β0+β1x的显著性检验,归结为对假设H0:
β1=0;H1:
β1≠0进行
检验。
假设H0:
β1=0被拒绝,则回归显著,认为y与x存在线性关系,所求的线性回归方
程有意义;否则回归不显著,y与x的关系不能用一元线性回归模型来描述,所得的回归方
程也无意义。
U
Qe/(n−2)
n
∑
F检验法:
当H0成立时,F=
~F(1,n−2),其中U=(yˆi−y)2(回
i=1
归平方和)。
若F>F1−α(1,n−2),拒绝H0,否则就接受H0。
回归系数的置信区间
β0和β1置信水平为1−α的置信区间分别为
⎡
⎤
⎥
⎦
1
+x,βˆ0+tα(n−2)σˆe
2
1
+x2
ˆ
⎢β0−tα(n−2)σˆe
和
n
Lxx
n
Lxx
⎢
⎣
1−
1−
⎥
2
2
⎡
⎢
⎣
⎤
⎥
β−t(n−2)σˆe/Lxx,βˆ1+t(n−2)σˆe/L;
ˆ
1
1−α2
1−α2
xx
⎦
7
⎡
⎢
⎢
⎤
⎥
⎥
Qe
Qe
σ
2
的置信水平为1−α的置信区间为
χ2(n−2),χ(n−2)
。
2
α
⎢
⎥
α
⎣
1−2
2
⎦
预测
用y的回归值yˆ0=βˆ0+βˆ1x
0
0作为0的预测值。
0的置信水平为
y
y
1−α的预测区间为
(x0−x)
Lxx
2
[yˆ0−δ(x0),yˆ0+δ(x0)],其中δ(x0)=σˆet(n−2)1+1+
。
特别,当n很
1−α
n
2
大且x0在x附近取值时,
y的置信水平为1−α的预测区间近似为
⎡
⎢
⎣
⎤
⎥
⎦
yˆ−σˆeu,yˆ+σˆeu
。
1−α2
1−α2
2.可线性化的一元非线性回归(需要配曲线)
先对两个变量x和y作n次试验观察得(xi,yi),i=1,2,...,n画出散点图,根据散点图确
定须配曲线的类型.然后由n对试验数据确定每一类曲线的未知参数a和b。
采用的方法是通
过变量代换把非线性回归化成线性回归,即采用非线性回归线性化的方法。
通常选择的六类曲线如下
(1)双曲线1=a+
b
y
x
(2)幂函数曲线y=axb,其中x>0,a>0
(3)指数曲线y=aebx其中参数a>0
(4)倒指数曲线y=aeb/x其中a>0
(5)对数曲线y=a+blnx,x>0
1
(6)S型曲线y=a+be−x
3.多元线性回归
⎧Y=Xβ+ε
一般称⎨
为高斯—马尔柯夫线性模型(k元线性回归模
⎩E(ε)=0,COV(ε,ε)=σ2
In
型),并简记为(Y,Xβ,σ2In)。
8
⎡y1⎤
⎡1x11x12...x1k⎤
⎡β0⎤
⎡ε⎤
1
⎢
⎥
⎥
⎢
⎥
⎥
⎢
⎢
⎥
⎥
⎢⎥
ε
...
1
x21x22...x2k
β
其中Y=⎢
,
X=⎢
,
β=
,
ε=⎢⎥
,
1
2
⎢...⎥
⎢...............⎥
⎢...⎥
⎢...⎥
⎢
⎣
⎥
⎦
⎢
⎣
⎥
⎦
⎢
⎣
⎥
⎦
⎢⎥
y
1
xn1xn2...xnk
β
ε
n
⎣⎦
n
k
y=β0+β1x1+...+βkxk称为回归平面方程。
线性模型(Y,Xβ,σ
2In)考虑的主要问题是:
(1)用试验值(样本值)对未知参数β和σ
作点估计和假设检验,从而建立y与x1,x2,...,xk
2
之间的数量关系;
(2)在x1=x01,x2=x02,...,xk=x0k处对y的值作预测与控制,即对y作区间估计。
对βi和σ
2
作估计,用最小二乘法求β0,...,βk的估计量:
作离差平方和
n
Q=
∑(
yi−β0−β1xi1−...−β
kxik)2
i=1
选择β0,...,βk使Q达到最小。
解得估计值βˆ=(XTX)−1(XTY),得到的βˆ代入回归平面
i
方程得y=βˆ0+βˆ1x1+...+βˆkxk,称为经验回归平面方程.βˆi称为经验回归系数
4.多元非线性回归