第八章插值 拟合与回归.docx
《第八章插值 拟合与回归.docx》由会员分享，可在线阅读，更多相关《第八章插值 拟合与回归.docx（24页珍藏版）》请在冰豆网上搜索。
第八章插值拟合与回归
第八章插值拟合与回归
在生产和实验中，由于函数
的表达式不便于计算或者没有表达式而只有在给定点的函数值（或其导数值），为此，我们希望建立一个简单的而且便于计算的近似函数
，来逼近函数
，这就用到插值和拟合方法。
8.1插值
插值就是在离散数据的基础上补插连续函数，使得这条连续曲线通过全部给定的离散数据点。
插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。
MATLAB中的插值函数主要有以下几个。
8.1.1interp1函数
MATLAB中用于一维数据插值的函数是interp1，其调用格式为：
yi=interp1（x,y,xi,'method,'extrap'）
该命令用于找出由参量x决定的一元函数y=y（x）在点xi处的值yi。
其中x,y为插值节点的横坐标和纵坐标，yi为在被插值点xi处的插值结果；x,y为向量，‘method’表示采用的插值方法，MATLAB提供的插值方法有几种：
‘nearest’：
最近邻点插值算法；
‘linear’：
线性插值（默认）；
‘spline’：
三次样条函数插值；
‘pchip’：
分段三次Hermite插值；
‘cubic’：
与’pchip’操作相同；
‘v5cubic’：
在MATLAB5.0中的三次插值。
缺省时表示线性插值。
‘extrap’表示对于超出x范围的xi中的分量将执行特殊的外插值法extrap。
注意：
所有的插值方法都要求x是单调的。
例1对函数
从0到10，步长为1进行采样得到插值节点，利用线性插值给出步长为0.25的插值函数图形。
在命令窗口输入：
>>x=0:
10;
y=x.*sin（x）;
xx=0:
.25:
10;
yy=interp1（x,y,xx）;
plot（x,y,'kd',xx,yy）
输出结果如图8.1所示。
例2在一天24小时内，从零点开始每间隔2小时测得的环境温度数据分别为
12，9，9，1，0，18，24，28，27，25，20，18，15，13
推测中午13时的温度。
我们采用三次样条插值，在命令窗口输入：
>>x=0:
2:
24;
a=13;
y=[129910182428272520181513];
y1=interp1（x,y,a,'spline'）
输出结果为：
y1=
27.8725
若要得到一天24小时的温度曲线，则输入：
xi=0:
1/3600:
24;
yi=interp1（x,y,xi,'spline'）;
plot（x,y,'o',xi,yi）
温度函数曲线图如图8.2所示。
图8.1线性插值函数曲线图图8.2温度曲线图
8.1.2interp2函数
MATLAB中用于二维数据插值的函数是interp2，其调用格式为：
zi=interp2（x,y,z,xi,yi,‘method’）
该命令用于找出由参量x,y决定的二元函数z=z（x,y）在点（xi,yi）处的值zi。
其中返回矩阵为zi，其元素为对应于参量xi与yi（可以是向量、或同型矩阵）的元素，若xi与yi中有在x与y范围之外的点，则相应地返回NaN（NotaNumber），method和interp1一样，常用的是‘cubic’（双三次插值），缺省为‘linear’（双线性插值算法）。
例3在命令窗口输入
[X,Y]=meshgrid（-3:
.25:
3）;
Z=peaks（X,Y）;%具有两个变量的采样函数，可产生一个凹凸有致的曲面，包含了三个局部极大点及三个局部极小点。
[XI,YI]=meshgrid（-3:
.125:
3）;
ZZ=interp2（X,Y,Z,XI,YI）;
surf（X,Y,Z）;holdon;
surf（XI,YI,ZZ+15）%为作比较，将插值曲面向上平移15单位。
holdoff
输出结果如图8.3
图8.3二维插值曲面
8.1.3interp3函数
MATLAB中用于三维数据插值的函数是interp3，其调用格式为：
vi=interp3（x,y,z,v,xi,yi,zi,‘method’）
该命令用于找出由参量x,y,z决定的三元函数v=v（x,y,z）在点（xi,yi,zi）处的值vi。
参量xi,yi,zi是同型阵列或向量。
若向量参量xi,yi,zi是不同长度，不同方向（行或列）的向量，这时输出参量vi与y1,y2,y3为同型矩阵。
其中y1,y2,y3为用命令meshgrid（XI,YI,ZI）生成的同型阵列。
若插值点（xi,yi,zi）中有位于点（x,y,z）之外的点，则相应地返回特殊变量值NaN。
method和interp1一样。
例4三维插值举例
命令窗口输入
[x,y,z,v]=flow（20）;
[xx,yy,zz]=meshgrid（.1:
.25:
10,-3:
.25:
3,-3:
.25:
3）;
vv=interp3（x,y,z,v,xx,yy,zz）;
slice（xx,yy,zz,vv,[69.5],[12],[-2.2]）;
shadinginterp;
colormapcool
输出图形如图8.4所示。
图8.4三维插值结果
8.1.4griddata函数
griddata也是一种常用的二维插值方法，其调用格式为：
zi=griddata（x,y,z,xi,yi,‘method’）
该命令用于找出由参量x,y决定的二元函数z=z（x,y）在点（xi,yi）处的值zi。
它和interp2的区别在于，interp2的插值数据必须是矩形域，即已知数据点（x,y）组成规则的矩阵，可使用meshgid生成。
而griddata函数的已知数据点（x,y）不要求规则排列，特别是对试验中随机没有规律采取的数据进行插值具有很好的效果。
method包括：
‘linear’（线性插值，为缺省算法）；‘natural’（自然邻点插值），‘nearest’（最近邻点插值）；‘cubic’（基于三角形的三次插值）和‘v4’（Matlab4中的griddata算法）。
例5在命令窗口输入：
rand（'seed',0）
x=rand（100,1）*4-2;y=rand（100,1）*4-2;
z=x.*exp（-x.^2-y.^2）;
ti=-2:
.125:
2;
[XI,YI]=meshgrid（ti,ti）;
ZI=griddata（x,y,z,XI,YI）;
mesh（XI,YI,ZI）,hold;
plot3（x,y,z,'o'）,holdoff;
输出结果如图8.5所示。
图8.5griddata插值结果
例6有一组散乱数据点矩阵如下：
A=[1.486,3.059,0.1;2.121,4.041,0.1;2.570,3.959,0.1;3.439,4.396,0.1;
4.505,3.012,0.1;3.402,1.604,0.1;2.570,2.065,0.1;2.150,1.970,0.1;
1.794,3.059,0.2;2.121,3.615,0.2;2.570,3.473,0.2;3.421,4.160,0.2;
4.271,3.036,0.2;3.411,1.876,0.2;2.561,2.562,0.2;2.179,2.420,0.2;
2.757,3.024,0.3;3.439,3.970,0.3;4.084,3.036,0.3;3.402,2.077,0.3;
2.879,3.036,0.4;3.421,3.793,0.4;3.953,3.036,0.4;3.402,2.219,0.4;
3.000,3.047,0.5;3.430,3.639,0.5;3.822,3.012,0.5;3.411,2.385,0.5;
3.103,3.012,0.6;3.430,3.462,0.6;3.710,3.036,0.6;3.402,2.562,0.6;
3.224,3.047,0.7;3.411,3.260,0.7;3.542,3.024,0.7;3.393,2.763,0.7];
x=A（:
1）;y=A（:
2）;z=A（:
3）;
[X,Y,Z]=griddata（x,y,z,linspace（min（x）,max（x）,20）',linspace（min（y）,max（y）,20）,'v4'）;
mesh（X,Y,Z）;holdon
plot3（x,y,z,'o'）;holdoff;
输出结果如图8.6所示。
图8.6散乱数据插值图
8.1.5spline函数
该函数是利用三次样条对数据进行插值，其调用格式为：
yy=spline（x,y,xx）
该命令用三次样条插值计算出由向量x与y确定的一元函数y=f（x）在点xx处的值。
若参量y是一矩阵，则以y的每一列和x配对，再分别计算由它们确定的函数在点xx处的值。
则yy是一个阶数为length（xx）*size（y,2）的矩阵。
例7对离散地分布在y=exp（x）sin（x）函数曲线上的数据点进行样条插值。
在命令窗口输入：
x=[024581212.817.219.920];
y=exp（x）.*sin（x）;
xx=0:
.25:
20;
yy=spline（x,y,xx）;
plot（x,y,'o',xx,yy）
输出结果如图8.7。
图8.7三次样条插值
8.2拟合
拟合就是用连续曲线近似地刻画或比拟平面上离散点组所表示的坐标之间的函数关系的一种数据处理方法。
如果已知某函数的若干离散函数值{f1,f2,…,fn}，通过调整该函数中若干待定系数f（λ1,λ2,…,λn），使得该函数与已知点集的差别（最小二乘意义）最小。
如果待定函数是线性，就叫线性拟合或者线性回归（主要在统计中），否则叫做非线性拟合或者非线性回归。
MATLAB中提供了线性最小二乘拟合和非线性最小二乘拟合函数。
8.2.1polyfit函数
ployfit函数是多项式拟合函数，其调用格式为：
p=polyfit（x,y,n）
其中x,y为长度相同的向量，n为拟合多项式的次数，返回值p是拟合多项式的系数向量，幂次由高到低，拟合多项式在x处的值可以通过y=polyval（p,x）来计算。
例8，2004年全国大学生数学建模竞赛C题，饮酒驾车问题，时间和血液中酒精浓度的函数关系，可以利用polyfit进行拟合。
可输入：
>>time=[0.25,0.50.75,11.522.533.544.55678910111213141516]';
vol=[2068758282776868585150413835282518151210774]';
plot（time,vol,'o'）
p=polyfit（time,vol,7）%7次多项式拟合
f=polyval（p,time）;%求多项式的值
holdon
plot（time,f）
holdoff
输出7多项次的系数向量为：
p=
0.0002-0.01130.2817-3.630725.5345-94.6710153.8757-1.2592
拟合曲线如图8.8所示。
图8.8血液酒精浓度拟合曲线
MATLAB提供了两个非线性最小二乘拟合函数：
lsqcurvefit和lsqnonlin。
两个命令都要先建立M-文件fun.m，在其中定义函数f（x），但两者定义f（x）的方式是不同的。
8.2.2lsqcurvefit函数
该函数用来进行非线性拟合，其调用格式为：
x=lsqcurvefit（'fun',x0,xdata,ydata,options）;
其中，fun为事先建立的拟合函数F（x,xdata），其中自变量x表示拟合函数中的待定参数，xdata为已知拟合节点的x坐标，x0为待定参数x的迭代初始值，xdata,ydata为已知数据点的x和y坐标，options是一些控制参数。
lsqcurvefit函数用来求含参数x（向量）的向量值函数
F（x,xdata）={f（x,xdata1）,f（x,xdata2）,…,f（x,xdatan）}
中的参数x（向量），使得
最小。
例9根据表1中的数据，利用lsqcurvefit函数拟合
。
表1已知数据点
X
100
200
300
400
500
600
700
800
900
1000
y×103
4.54
4.99
5.35
5.65
5.90
6.10
6.26
6.39
6.50
6.59
首先建立拟合函数的M文件myfit1.m，其内容如下：
functionf=myfit1（x,xdata）
f=x
（1）+x
（2）*exp（-0.02*x（3）*xdata）;
其中x
（1）,x
（2）,x（3）分别表示拟合曲线中的参数a,b,k。
然后在命令窗口输入：
>>xdata=100:
100:
1000;
ydata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];
x0=[0.2,0.05,0.05];
x=lsqcurvefit（'myfit1',x0,xdata,ydata）
f=myfit1（x,xdata）
输出结果为：
x=
0.0069-0.00290.0809
f=
Columns1through9
0.00440.00480.00510.00540.00560.00580.00600.00610.0063
Column10
0.0064%拟合函数在数据点的函数值
拟合曲线如图8.9。
>>plot（xdata,ydata,'o',xdata,f）
图8.9拟合曲线
8.2.3lsqnonlin函数
该函数用来进行非线性拟合，其调用格式为：
x=lsqnonlin（'fun',x0,options）;
其中fun为事先建立的拟合函数F（x），其中自变量x表示拟合函数中的待定参数，x0为待定参数x的迭代初始值，options是一些控制参数。
由于lsqnonlin中定义的拟合函数的自变量是x，所以已知参数xdata,ydata应写在该函数中。
lsqnonlin函数用来求含参量x（向量）的向量值函数
f（x）={f1（x）,f2（x）,…，fn（x）}
中的参量x，使得
最小。
其中fi（x）=f（x,xdatai,ydaytai）=F（x,xdatai）-ydatai。
例10根据表1中的数据，利用lsqnonlin函数拟合
。
首先建立拟合函数的M文件myfit2.m，其内容如下：
functionf=myfit2（x）
xdata=100:
100:
1000;
ydata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];
f=x
（1）+x
（2）*exp（-0.02*x（3）*xdata）-ydata;
然后在命令窗口输入：
>>x0=[0.2,0.05,0.05];
x=lsqnonlin（'myfit2',x0）
f=myfit2（x）
输出结果为：
x=
0.0069-0.00290.0809
f=
1.0e-03*
Columns1through9
-0.0971-0.1739-0.2164-0.2463-0.2666-0.2713-0.2651-0.2538-0.2436
Column10
-0.2313
该结果与lsqcurvefit函数拟合的结果相同。
8.3回归分析
MATLAB中提供了一些线性和非线性回归分析函数。
8.3.1regress函数
一般地，称由
确定的模型为一元线性回归模型，记为
固定的未知参数
、
称为回归系数，自变量x也称为回归变量。
一般地，称
为高斯—马尔柯夫线性模型（k元线性回归模型）。
，
，
，
MATLAB提供了多元线性回归函数regress，采用的是最小二乘估计，其调用格式有：
b=regress（y,x）
返回值为线性模型y=x*b的回归系数向量。
其中x为n×（k+1）矩阵，行对应于观测值，列对应于预测变量，y为n维向量，为因变量，一元线性回归可取k=1。
[b,bint,r,rint,stats]=regress（y,x,alpha）
其中bint是回归系数的区间估计，r是残差，rint是置信区间，stats是用于检验回归模型的统计量，有四个数值：
相关系数r^2，F值，与F对应的概率P和误差方差估计，alpha是显著性水平（缺省的时候为0.05）。
相关系数r^2越接近1，说明回归方程越显著；与F对应的概率P例11输入
>>y=[8.88188.94879.05419.15459.26939.42899.61609.81509.982510.155810.3193]';
k=[7.83817.91678.00488.10268.25568.58228.82879.07569.21759.41489.6198]';
l=[8.38718.38728.39358.39718.40258.40488.40798.41418.42618.43778.4444]';
t=[9.95519.905710.09729.95379.93709.94499.963610.129110.157310.294410.2093]';
x=[ones（size（k））klt];
[b,bint,r,rint,stats]=regress（y,x）
输出结果为
b=
-55.4988
0.5644
7.1254
0.0222
bint=
-83.6008-27.3968
0.46930.6594
3.574010.6769
-0.19140.2359
r=
-0.0262
-0.0032
0.0033
0.0261
0.0164
-0.0249
0.0006
0.0124
0.0137
-0.0101
-0.0081
rint=
-0.06200.0097
-0.04790.0415
-0.03060.0372
-0.01110.0632
-0.02180.0546
-0.06210.0124
-0.03850.0397
-0.01900.0438
-0.03040.0577
-0.04600.0259
-0.04200.0258
stats=
1.0e+003*
0.00102.10220.00000.0000
该结果说明y=-55.4988+0.5644k+7.1254l+0.0222t，stats中的数据说明r2=1，F=2102.2，p=0，由于p<0.05可知回归模型成立。
可以利用rcoplot函数画出残差及其置信区间，红色的表示超出期望值的数据，圆圈代表残差的值，竖线代表置信区间的范围。
输入
>>rcoplot（r,rint）
输出残差图见8.10。
图8.10残差图
8.3.2rstool函数
该函数是多元二项式回归函数，其调用格式为
rstool（x,y,'model',alpha）
其中x为n×m为矩阵，y为n维列向量，'model'为以下4种模型：
'linear'（线性，缺省）：
；
'interaction'（交叉）：
；
'quadratic'（完全二次）：
；
'purequadratic'（纯二次）：
。
alpha为显著性水平，默认值为0.05。
例13设某商品的需求量与消费者的平均收入、商品价格的统计数据如表2，建立多元二项式纯二次回归模型，并预测平均收入为1000、价格为6时的商品需求量。
表2需求量、平均收入和价格统计表
需求量
100
75
80
70
50
65
90
100
110
60
收入
1000
600
1200
500
300
400
1300
1100
1300
300
价格
5
7
6
6
8
7
5
4
3
9
可以直接使用多元二项式回归，在命令窗口输入：
>>x1=[10006001200500300400130011001300300];
x2=[5766875439];
y=[10075807050659010011060]';
x=[x1'x2'];
rstool（x,y,'purequadratic'）
其输出如图8.11。
图8.11多元二项式回归
在图8.11中x1上面的方框中输入1000，x2上面的方框中输入6，在图形框左侧的“PredictedY1”下方的数据变为88.47981，即预测出平均收入为1000，价格为6时的商品需求量为88.4791。
单击图形框左边Export，则出现图8.12对话框，可以将回归参数beta、剩余标准差rmse和残差residuals传送到MATALB的工作区中。
图8.12输出对话框
在命令窗口输入：
>>beta
rmse
residuals
输出结果为：
beta=
110.5313
0.1464
-26.5709
-0.0001
1.8475
rmse=
4.5362
residuals=
5.2724
-0.7162
-4.5158
-1.9390
-3.3315
3.4566
3.4843
-3.4452
-0.0976
1.8320
故回归模型为：
，剩余标准差为4.5362,说明此回归模型的显著性较好。
还可将该模型转化为多元线性回归模型，利用regress函数进求解。
可输入：
>>x1=[10006001200500300400130011001300300];
x2=[5766875439];
X=[ones（10,1）x1'x2'（x1.^2）'（x2.^2）'];
[b,bint,r,rint,stats]=regress（y,X）;
b,stats
输出结果为：
b=
110.5313
0.1464
-26.5709
升级会员 