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')运行结果：0246810125101520253035HourDegreesCelsius2.二维插值（自变量是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）XY12001600200024002800320036004000120011301250128012301040900500700160013201450142014001300700900850200013901500150014009001100106095024001500120011001350145012001150101028001500120011001550160015501380107032001500155016001550160016001600155036001480150015501510143013001200980试作出该山区的地貌图和等高线图。代码：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）里的哪些地方船要避免进入．xyz1297.54140141.58103.5236881478185.522.56195137.5810585.58157.5-6.59107.5-81977388156.58162-66.59162844117.5-33.59作海底曲面图，作出水深小于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)pausez1(z1=5)=nan;%将水深大于5的置为nan，这样绘图就不会显示出来meshc(x1,y1,z1)运行结果：50100150200-5005010015024681012二、数据拟合拟合最常用的方法为最小二乘法。原理是：因变量的实际值与拟合值之差，称为残差（拟合误差），将所有残差平方后相加，即得误差平方和，“最好”拟合效果就是使误差平方和最小，于是运用极值原理，将“求最好拟合问题”转换为“求误差平方和最小”问题。Matlab实现：1.多项式拟合f(x)=a1xm+…+amx+am+1语法：50100150200-500501001503.544.55a=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年我国人口数据资料如下：年份x01949195419591964196919741979198419891994人口数y05.46.06.77.08.19.19.810.311.311.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,'*')holdonplot(x1,y1,'--r')holdonplot(x1,y2,'-k')legend('原数据散点图','线性拟合曲线','指数拟合曲线')运行结果：y1_2005=13.5080y2_2005=15.0596注：（1）实际2005年人口是13.3亿，可见线性拟合效果更好；（2）人口指数增长模型为：bxyae，两边取ln得lny=lna+bx令Y=lny，A=lna，B=b，则转化为直线方程：Y=A+Bx故可以作线性拟合得到A，B值，从而得到ABxyee.一般可以转化为直线方程的曲线有：19401950196019701980199020005678910111213原数据散点图线性拟合曲线指数拟合曲线（Logistic生长曲线，就是S型曲线）2.多元线性回归011nnyxxβ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.5234.57.59.110.512体重5.66.67.27.810.110.813.516.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.15750.8950——回归方程：y=4.1575+0.8950xbint=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.10.65250.3575-0.3850-0.7701-1.5021-0.05511.6024rint=-2.12832.3283——0.1000∈[-2.1283,2.3283]-1.52822.8332-2.00742.7223-2.84452.0744-3.14011.5999-3.29350.2893-