第四讲 插值、拟合与回归分析在生产实践和科学研究中，常常有这样的问题：由实验或测量得到变量间的一批离散样本点，要求得到变量之间的函数关系或得到样本点之外的数据。

解决此类问题的方法一般有插值、拟合和回归分析等。

设有一组实验数据0011(,),(,),(,)n n x y x y x y ，当原始数据精度较高，要求确定一个简单函数()y x ϕ=（一般为多项式或分段多项式）通过各数据点，即(),0,,i i y x i n ϕ== ，称为插值问题。

另一类是拟合问题，当我们已经有了函数关系的类型，而其中参数未知或原始数据有误差时，我们确定的初等函数()y x ϕ=并不要求经过数据点，而是要求在某种距离度量下总体误差达到最小，即(),0,,i i i y x i n ϕε=+= ，且20ni i ε=∑达到最小值。

对同一组实验数据，可以作出各种类型的拟合曲线，但拟合效果有好有坏，需要进行有效性的统计检验，这类问题称为回归分析。

一、插值(interpolation)常用的插值方法有分段线性插值、分段立方插值、样条插值等。

1、一元插值yi=interp1(x,y,xi,method)对给定数据点(x,y)，按method 指定的方法求出插值函数在点（或数组）xi 处的函数值yi 。

其中method 是字符串表达式，可以是以下形式：'nearest' ——最邻近点插值'linear' ——分段线性插值（也是缺省形式）'spline' ——分段三次样条插值'cubic' 分段立方插值例：在一天24小时内，从零点开始每间隔2小时测得环境温度数据分别为（℃）：12，9，9，10，18，24，28，27，25，20，18，15，13用不同的插值方法估计中午1点（即13点）的温度，并绘出温度变化曲线。

>> x=0:2:24;>> y=[12 9 9 10 18 24 28 27 25 20 18 15 13];>>y_linear=interp1(x,y,13),y_nearest=interp1(x,y,13,'nearest')>>y_cubic=interp1(x,y,13,'cubic'),y_spline=interp1(x,y,13,'spline')>> y1=interp1(x,y,xx); y2=interp1(x,y,xx,'nearest');>> y3=interp1(x,y,xx,'cubic');y4=interp1(x,y,xx,'spline');>> subplot(2,2,1),plot(x,y,'or',xx,y1)>> subplot(2,2,2),plot(x,y,'or',xx,y2)>> subplot(2,2,3),plot(x,y,'or',xx,y3)>> subplot(2,2,4),plot(x,y,'or',xx,y4)2、二元插值zi=interp2(X,Y,Z,xi,yi,method)已知数据点（X，Y，Z），求插值函数在（xi,yi）处的函数值zi,插值方法method同interp1。

这里要求X,Y,Z是同维矩阵，且X，Y是网格矩阵，或者X是与Z列数相同的行向量，Y是与Z行数相同的列向量。

例：测得平板表面5 3网格点处的温度分别为试作出平板表面的温度分布图>> x=1:5;y=1:3;z=[82 81 80 82 84;79 63 61 65 81;84 84 82 85 86];>> xx=1:0.1:5;yy=1:0.1:3;yy=yy';>> zz=interp2(x,y,z,xx,yy,'cubic');>> mesh(xx,yy,zz)3、不规则点的插值若数据是不规则的，即数据不能构成矩阵形式，从而不能用interp2函数进行插值。

zi=griddata(x,y,z,xi,yi,method)这里,x,y,z为同维向量，表示已知数据点的坐标，xi,yi是行向量和列向量，返回值zi为在meshgrid(xi,yi)网格矩阵处的函数值。

method 可选择’linear’,’nearest’,’cubic’。

例：假如上例中的数据残缺不全>> x=[3 4 5 1 3 4 1 2 5];y=[1 1 1 2 2 2 3 3 3];>>z=[80 82 84 79 61 65 84 84 86];>> xx=1:0.1:5;yy=1:0.1:3;>> zz=griddata(x,y,z,xx,yy','cubic');>> mesh(xx,yy',zz)二、拟合(Fit)1、多项式拟合p=polyfit(x,y,n) 用n次多项式拟合向量数据(x,y)。

例：拟合下列数据>> x=[0.1 0.2 0.15 0 -0.2 0.3];y=[0.95 0.84 0.86 1.06 1.50 0.72];>> p=polyfit(x,y,2);>> xx=-0.2:0.01:0.3;yy=polyval(p,xx);>> plot(x,y,'or',xx,yy)2、曲线拟合当经验函数不是多项式，而是其它类型的函数时，可以用lsqcurvefit 函数对拟合函数中的未知参数进行估计。

c=lsqcurvefit(fun,c0,xdata,ydata)fun 是经验拟合函数，含有未知参数，即具有形式fun(c,x)，c0是未知参数的预估计值，(xdata,ydata)是已知实验数据。

例：已知数据表用适当的曲线进行数据拟合。

先画散点图，根据散点图确定拟合曲线为对数函数ln b t y a += >> t=1:16;>> y=[4 6.4 8 8.4 9.28 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 ... 10.55 10.58 10.6]; >> plot(t,y,'or')>> f=inline('c(1)+c(2)*log(t)','c','t') %建立拟合函数 >> c=lsqcurvefit(f,[1,1],t,y) %求未知参数 >> tt=1:0.1:16;yy=f(c,tt); >> hold on >> plot(tt,yy) 3、拟合工具箱Matlab 中的拟合工具箱是一个更方便、更直观进行曲线拟合的图形界面，用cftool 指令打开拟合工具箱。

拟合效果主要看2个参数：SSE （误差平方和）和R-Square ，SSE 越接近0，R-Square 越接近1，拟合效果越好。

三、多元线性回归问题：设有因变量y 和p 个自变量12,,p x x x ，它们具有某种线性关系 1122p p y x x x βββε=+++ 其中12,,p βββ 为待定系数，ε为随机误差。

现有容量为n 观测数据,,1,2,,1,2,i ij y x i n j p == ，怎样确定待定系数12,,p βββ ，并进行有效性检验？将样本代入关系式，得Y X βε=+，其中1111211221222212,,p p p n n np p y x x x y x x x Y X y x x x ββββ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦Matlab 求解：[b,bint,r,rint,stats]=regress(Y ,X,alpha)其中alpha 为显著性水平，省略时为默认值0.05；输出变量中，b 为β的参数估计值，bint 为b 的置信度为1-alpha 的置信区间，r 为残差向量Y X β-，也即ε，rint 为ε的置信区间，stats 是包含3个元素的检验统计量，分别是R-square ：相关系数R 的平方，F-统计量和p 值。

回归效果：R-Square 越接近1，p 值越接近0（一般要求p<0.05）。

例：某种水泥在凝固时放出的热量（单位:卡/克）Y 与水泥中下列4种化学成分所占的百分比有关：x1:233CaO Al O ⋅; x2: 23CaO SiO ⋅; x3:23234.CaO Al O Fe O ⋅; x4:22CaO SiO ⋅ 现测得13组数据，见表，要求建立热量与水泥化学成分之间的经验回归关系式。

输入数据，可以先建立2个全零矩阵x=zeros(13,4);y=zeros(13,1);然后将表中的数据直接复制、粘贴到相应位置。

>> [b,bint,r,rint,stats]=regress(y,x) 最后得到的回归方程为：12342.193 1.15330.75850.4863y x x x x =+++如果回归方程是形式：01122p p y x x x ββββε=++++ ，相当于增加一个变量0x ，001122p p y x x x x ββββε=++++ ，而01x ≡。

如上例>> x=[ones(13,1),x];>> [b,bint,r,rint,stats]=regress(y,x) 得回归关系式：123462.4054 1.55110.51020.10190.1441y x x x x =+++-上机练习1、在1-12的11小时内，每隔1小时测量一次温度，测得的温度依次为：5，8，9，15，25，29，31，30，22，25，27，24。

试估计每隔1/10小时的温度值。

2、已知飞机下轮廓线上数据如下，求x 每改变0.1时的y 值。

X 035791112131415Y0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.63、在某海域测得一些点(x,y)处的水深z 由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）*（-50，150）里的哪些地方船要避免进入。

4、山区地貌：在某山区测得一些地点的高程如下表：(平面区域1200<=x<=4000,1200<=y<=3600)，试作出该山区的地貌图和等高线图，并对几种插值方法进行比较。