Matlab学习系列13 数据插值与拟合.docx
《Matlab学习系列13 数据插值与拟合.docx》由会员分享，可在线阅读，更多相关《Matlab学习系列13 数据插值与拟合.docx（21页珍藏版）》请在冰豆网上搜索。
Matlab学习系列13数据插值与拟合
13.数据插值与拟合
实际中，通常需要处理实验或测量得到的离散数据（点）。
插值与拟合方法就是要通过离散数据去确定一个近似函数（曲线或曲面），使其与已知数据有较高的拟合精度。
1.如果要求近似函数经过所已知的所有数据点，此时称为插值问题（不需要函数表达式）。
2.如果不要求近似函数经过所有数据点，而是要求它能较好地反映数据变化规律，称为数据拟合（必须有函数表达式）。
插值与拟合都是根据实际中一组已知数据来构造一个能够反映数据变化规律的近似函数。
区别是：
【插值】不一定得到近似函数的表达形式，仅通过插值方法找到未知点对应的值。
【拟合】要求得到一个具体的近似函数的表达式。
因此，当数据量不够，但已知已有数据可信，需要补充数据，此时用【插值】。
当数据基本够用，需要寻找因果变量之间的数量关系（推断出表达式），进而对未知的情形作预测，此时用【拟合】。
一、数据插值
根据选用不同类型的插值函数，逼近的效果就不同，一般有：
（1）拉格朗日插值（lagrange插值）
（2）分段线性插值
（3）Hermite
（4）三次样条插值
Matlab插值函数实现：
（1）interp1（）一维插值
（2）intep2（）二维插值
（3）interp3（）三维插值
（4）intern（）n维插值
1.一维插值（自变量是1维数据）
语法：
yi=interp1（x0,y0,xi,‘method’）
其中，x0,y0为原离散数据（x0为自变量，y0为因变量）；xi为需要插值的节点，method为插值方法。
注：
（1）要求x0是单调的，xi不超过x0的范围；
（2）插值方法有‘nearest’——最邻近插值；‘linear’——线性插值；‘spline’——三次样条插值；‘cubic’——三次插值；
默认为分段线性插值。
例1从1点12点的11小时内，每隔1小时测量一次温度，测得的温度的数值依次为：
5，8，9，15，25，29，31，30，22，25，27，24．试估计每隔1/10小时的温度值。
代码：
hours=1:
12;
temps=[589152529313022252724];
h=1:
0.1:
12;
t=interp1（hours,temps,h,'spline'）;
plot（hours,temps,'+',h,t,hours,temps,'r:
'）
%作图，原散点数据用+标记
xlabel（'Hour'）,ylabel（'DegreesCelsius'）
运行结果：
2.二维插值（自变量是2维数据）
语法：
zi=interp2（x0,y0,z0,xi,yi,‘method’）
其中，x0,y0,z0为原离散数据（x0,y0为自变量，z0为因变量）；xi,yi为需要插值的节点，method为插值方法。
注：
（1）要求x0,y0是单调的，xi,yi不超过x0,y0的范围，注意：
要求xi取为行向量，yi取列向量；
（2）插值方法有‘nearest’——最邻近插值；‘linear’——双线性插值；‘cubic’——双三次插值；默认为双线性插值。
例2山区地貌：
在某山区测得一些地点的高度如下表：
（平面区域为：
1200<=x<=4000,1200<=y<=3600）
X
Y
1200
1600
2000
2400
2800
3200
3600
4000
1200
1130
1250
1280
1230
1040
900
500
700
1600
1320
1450
1420
1400
1300
700
900
850
2000
1390
1500
1500
1400
900
1100
1060
950
2400
1500
1200
1100
1350
1450
1200
1150
1010
2800
1500
1200
1100
1550
1600
1550
1380
1070
3200
1500
1550
1600
1550
1600
1600
1600
1550
3600
1480
1500
1550
1510
1430
1300
1200
980
试作出该山区的地貌图和等高线图。
代码：
x0=1200:
400:
4000;
y0=1200:
400:
3600;
height=[11301250128012301040900500700;
13201450142014001300700900850;
139015001500140090011001060950;
15001200110013501450120011501010;
15001200110015501600155013801070;
15001550160015501600160016001550;
1480150015501510143013001200980];
mesh（x0,y0,height）
title（'原数据散点图'）;
pause%等待，按任意键切换到新图形
xi=1200:
5:
4000;
yi=1200:
5:
3600;
zi=interp2（x0,y0,height,xi',yi,'cubic'）;
mesh（xi,yi,zi）
title（'cubic插值图'）;
运行结果：
另外，二维插值函数还有griddata函数，格式类似
cz=griddata（x0,y0,z0,cx,cy,‘method’）
插值方法有‘nearest’——最邻近插值；‘linear’——双线性插值；‘cubic’——双三次插值；‘v4’——Matlab提供的插值方法；默认为双线性插值。
与interp2区别是，interp2的插值数据必须是矩形域，即已知数据点（x,y）组成规则的矩阵（或称之为栅格），可使用meshgid生成。
而griddata函数的已知数据点（X，Y）不要求规则排列，特别是对试验中随机没有规律采取的数据进行插值具有很好的效果。
例3在某海域测得一些点（x,y）处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）×（-50，150）里的哪些地方船要避免进入．
x
y
z
129
7.5
4
140
141.5
8
103.5
23
6
88
147
8
185.5
22.5
6
195
137.5
8
105
85.5
8
157.5
-6.5
9
107.5
-81
9
77
3
8
81
56.5
8
162
-66.5
9
162
84
4
117.5
-33.5
9
作海底曲面图，作出水深小于5的海域范围,即z=5的等高线。
代码：
x=[129.0140.0103.588.0185.5195.0105.5157.5107.577.081.0162.0162.0117.5];
y=[7.5141.523.0147.022.5137.585.5-6.5-813.056.5-66.584.0-33.5];
z=[48686889988949];
x1=75:
1:
200;
y1=-50:
1:
150;
[x1,y1]=meshgrid（x1,y1）;
z1=griddata（x,y,z,x1,y1,'v4'）;
meshc（x1,y1,z1）
pause
z1（z1>=5）=nan;%将水深大于5的置为nan，这样绘图就不会显示出来
meshc（x1,y1,z1）
运行结果：
二、数据拟合
拟合最常用的方法为最小二乘法。
原理是：
因变量的实际值与拟合值之差，称为残差（拟合误差），将所有残差平方后相加，即得误差平方和，“最好”拟合效果就是使误差平方和最小，于是运用极值原理，将“求最好拟合问题”转换为“求误差平方和最小”问题。
Matlab实现：
1.多项式拟合
f（x）=a1xm+…+amx+am+1
语法：
a=polyfit（x0,y0,m）;
y=polyval（a,x）;
注：
（1）x，y的长度相同，m为拟合多项式次数，m=1时，一次多项式对应一条直线，也称为线性拟合（一元线性回归）；
（2）a返回的是拟合多项式的系数：
[a1,a2,…,am];
（3）polyval函数用来计算拟合多项式在x点处的取值；
（4）拟合多项式的次数一般为1-3次，不宜过高，否则有龙格现象（两端出现很大误差）。
例41949年—1994年我国人口数据资料如下：
年份x0
1949
1954
1959
1964
1969
1974
1979
1984
1989
1994
人口数y0
5.4
6.0
6.7
7.0
8.1
9.1
9.8
10.3
11.3
11.8
分析我国人口增长的规律，预报2005年我国人口数。
代码：
x0=[1949195419591964196919741979198419891994];
y0=[5.46.06.77.08.19.19.810.311.311.8];
a=polyfit（x0,y0,1）;%线性拟合
x1=[1949:
0.5:
1994];
y1=polyval（a,x1）;
%或者直接用多项式计算：
y1=a
（1）*x1+a
（2）;
y1_2005=polyval（a,2005）
b=polyfit（x0,log（y0）,1）;%指数函数拟合
y2=exp（b
（2））*exp（b
（1）*x1）;
y2_2005=exp（b
（2））*exp（b
（1）*2005）
plot（x0,y0,'*'）
holdon
plot（x1,y1,'--r'）
holdon
plot（x1,y2,'-k'）
legend（'原数据散点图','线性拟合曲线','指数拟合曲线'）
运行结果：
y1_2005=13.5080
y2_2005=15.0596
注：
（1）实际2005年人口是13.3亿，可见线性拟合效果更好；
（2）人口指数增长模型为：
，两边取ln得
lny=lna+bx
令Y=lny，A=lna，B=b，则转化为直线方程：
Y=A+Bx
故可以作线性拟合得到A，B值，从而得到
.一般可以转化为直线方程的曲线有：
（Logistic生长曲线，就是S型曲线）
2.多元线性回归
βi为回归系数，ε为残差。
n=1时为线性回归。
Matlab实现：
[b,bint,r,rint,stats]=regress（y,x,α）;
注：
（1）整体是一个假设检验（是否符合多元线性关系）：
原假设H0：
βi=0（即不存在线性关系）
备择假设H1：
βi≠0（回归式成立）；
（2）若有常数项β0，则x的第一列必须全为1；
（3）α为显著性水平的概率，默认是0.05;
（4）b返回β0,β1,…βn值的点估计，
bint返回β0,β1,…βn值的区间估计，
r返回残差的点估计,
rint返回残差的区间估计；
（5）stats返回4个数：
相关系数平方r2（越大越好，若小于0.5结果无价值）
F统计量（用于计算概率）
拒绝原假设的概率（<α时拒绝）（越小越好）
剩余方差=残差平方和/[数据个数-（n+1）]（越小越好）
其中，n+1为回归系数β0,β1,…βn的个数。
例5已知某动物的进食量与体重数据为
进食量
1.5
2
3
4.5
7.5
9.1
10.5
12
体重
5.6
6.6
7.2
7.8
10.1
10.8
13.5
16.5
拟合直线y=β0+β1x.
代码：
food=[1.5234.57.59.110.512];
n=length（food）;
x=[ones（n,1）,food'];
y=[5.66.67.27.810.110.813.516.5]';
[b,bint,r,rint,stat]=regress（y,x,0.05）
运行结果：
b=4.1575
0.8950——>回归方程：
y=4.1575+0.8950x
bint=2.46925.8458——>4.1575∈[2.4692,5.8458]
0.66440.6644——>0.8950∈[0.6644,0.6644]
r=0.1000——>残差1:
5.6-（4.1575+0.8950*1.5）=0.1
0.6525
0.3575
-0.3850
-0.7701
-1.5021
-0.0551
1.6024
rint=-2.12832.3283——>0.1000∈[-2.1283,2.3283]
-1.52822.8332
-2.00742.7223
-2.84452.0744
-3.14011.5999
-3.29350.2893
-2.35192.2417
0.48462.7201
stat=0.937690.18710.00011.0220
——>r2=0.9376>0.5（说明结果很好）
0.0001<α=0.05（拒绝原假设，接受备择假设）
【残差r的平方和/[数据个数-（n+1）]】：
（0.10002+0.65252+…+1.60242）/[8-（1+1）]=1.0220
3.一般曲线拟合（非线性回归）
曲线形式y=f（a,x）已知，但其中a=[a1,a2,…]为未知的待定系数。
例如，已知曲线形式为y=f（a,x）=a1x2+a2sinx+a3x3，以及原始数据x1,…,xn，y1,…,yn，利用最小二乘法原理，求出让误差平方和
[F（a,x1）-y1]2+…+[F（a,xn）-yn]2
取最小值的a1，a2，a3.从而得到拟合曲线方程。
Matlab实现：
lsqcurvefit和lsqnonlin函数。
（1）lsqcurvefit函数
语法：
先定义形式函数f（a,x）的M-文件fun.m;
再调用
[x,resnorm,residual]=lsqcurvefit（‘fun’,x0,xdata,ydata,lb,ub,options）;
注：
（1）返回值x为拟合出的待定系数，resnorm返回残差的平方和（越小越好），residual返回残差；
（2）‘fun’为预定义的形式函数，xdata,ydata为已知数据，x0为迭代初始值（初始值选取的不同将影响拟合结果），lb，ub是可选参数，为x指定范围“下界和上界”；
（3）options是可选参数，指定采用的优化参数，语法：
OPTIONS=OPTIMSET（'PARAM1',VALUE1,'PARAM2',VALUE2,...）
示例：
options=optimset（'TolFun',1e-3）;
注意：
形式函数是以待定系数a和已知数据x作为参数，lsqcurvefit函数的参数中包含已知数据x，y.
例6用下面一组数据拟合
中的待定系数a，b，k.
tj
100
200
300
400
500
600
700
800
900
1000
cj（*10-3）
4.54
4.99
5.35
5.65
5.90
6.10
6.26
6.39
6.50
6.59
该问题即解最优化问题：
满足下式的的a，b，k:
代码：
（lsqcurvefit函数）
%%%%%%%%%%函数单独存为m文件%%%%%%%%%%%%%%%%%%%%%%%%
functionf=curvefun1（x,tdata）%tdata是向量
f=x
（1）+x
（2）*exp（-0.02*x（3）*tdata）;
%其中x
（1）=a;x
（2）=b;x（3）=k;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tdata=[100:
100:
1000]';
cdata=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,resnorm,residual]=lsqcurvefit（'curvefun1',x0,tdata,cdata）
C_fit=curvefun1（x,tdata）
plot（tdata,cdata,'o',tdata,f,'r'）
运行结果：
x=0.0069-0.00290.0809
——>拟合曲线：
resnorm=5.3938e-007（残差平方和）
residual=1.0e-003*（残差）
-0.0971-0.1739-0.2164-0.2463-0.2666-0.2713-0.2651-0.2538-0.2436-0.2313
c_fit=0.00440.00480.00510.00540.00560.00580.00600.00610.00630.0064
（2）lsqnonlin函数
与lsqcurvefit函数用法类似，不同在于形式函数只以待定系数a作为参数，但函数体里用到已知数据x,y.而在调用拟合函数时不再需要已知数据x,y作为参数：
[x,resnorm,residual]=lsqnonlin（‘fun’,x0,lb,ub,options）;
对于例6，调用lsqnonlin函数实现的代码如下：
（运行结果一样）
代码：
（lsqnonlin函数）
%%%%%%%%%%%%%%%%函数单独存为m文件%%%%%%%%%%%%%%%%%%
functionf=curvefun2（x）
tdata=[100:
100:
1000]';cdata=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）*tdata）-cdata;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x0=[0.2,0.05,0.05]';
[x,resnorm,residua]=lsqnonlin（'curvefun2',x0）
（3）nlinfit函数
非线性拟合的通用函数，适用面比lsqcurvefit函数广，例如可以做加权最小二乘拟合。
语法：
回归：
[beta,r,J]=nlinfit（x,y,‘model’,beta0）;
预测：
[Y,DELTA]=nlpredci（‘model’,x,beta,r,J）;
注：
（1）beta返回拟合的待定系数；r返回残差；J返回拟合函数的雅可比矩阵；
（2）x、y为已知数据，分别为n×m矩阵和n维列向量，对一元非线性回归，x为n维列向量；
（3）‘model’为预先定义的形式函数；beta0为指定的待定系数初始值；
（4）nlpredic函数用来求nlinfit所得的拟合函数在x处的预测值Y及预测值的显著性为1-alpha的置信区间Y±DELTA.
例7混凝土的抗压强度随养护时间的延长而增加，现将一批混凝土作成12个试块，下表记录了养护日期x（日）及抗压强度y（kg/cm2）的数据：
养护时间x
2
3
4
5
7
9
12
14
17
21
28
56
抗压强度y（+r）
35
42
47
53
59
65
68
73
76
82
86
99
建立非线性回归模型：
对其中的待定系数a，k1，k2，m进行拟合（回归），并对结果做检验。
注：
表中的+r代表加上一个[-0.5,0.5]之间的随机数。
代码：
（使用了内联函数，也可以如前例将函数定义为m文件）
x=[234579121417212856];
r=rand（1,12）-0.5;
y1=[354247535965687376828699];
y=y1+r;
myfunc=inline（'beta
（1）+beta
（2）*exp（beta（4）*x）+beta（3）*exp（-beta（4）*x）','beta','x'）;
[beta,r,J]=nlinfit（x,y,myfunc,[0.50.50.50.5]）
[y_fit,delta]=nlpredci（myfunc,x',beta,r,J）
plot（x,y,'o',x,y_fit,'r'）
运行结果：
beta=87.57900.0283-62.72140.1070
r=-1.76570.13760.03782.28821.17550.9553-2.6909-0.6823-1.11950.62541.0755-0.0368
J=1.00001.23860.8074101.3489
1.00001.37840.7255136.6213
1.00001.53410.6518163.7126
1.00001.70730.5857183.9240
1.00002.11470.4729208.0387
1.00002.61920.3818216.1866
1.00003.61050.2770209.6925
1.00004.47190.2236198.1339
1.00006.16430.1622175.9447
1.00009.45660.1057144.9124
1.000019.99780.0500103.6939
1.0000399.91160.0025643.7199
y_fit=36.975142.117146.738350.891557.979663.707370.309973.680777.579381.215085.009898.7606
delta=2.66
升级会员 