第1讲插值、曲线拟合与回归分析

1、第9讲 插值、曲线拟合和回归分析研究生E题，罗伦滋曲线拟合。多项式曲线拟合x=0:0,1:5;y=x.A2+rand(1,length(x);p=polyfit(x,y,2);%2次多项式曲线拟合f=poly2sym(p,t);%以多项式系数p生成多项式函数f=inline(f); %f(x)相当匿名函数plot(x,f(x)syms x;sym2poly(f(x);%返回多项式函数的系数向量p在实际中，常常要处理有实验或测量所得到的一批离散数据。插值与拟合方法就是要通过这些数据去决定莫一类已知函数的参数或寻找莫个近似函数， 使得近似函数与已知数据有较高的拟合精度。如果要求这个近似函数(曲线或

2、曲面)经过所有已知数据点，则这类问题称为插值问题。这种寻找函数的方法称为插值i 方法。如果不要求这个近似函数(曲线或曲面)经过所有已知数据点，而是要求它能较好的反映数据的整体变化趋势，称解决这类问题的方法为数据拟合。共同点：都是根据一组已知数据来构造反映数据变化规律 的近似函数的方法。不同点：由于对近似函数要求不同，二者在数学方法上有很大差异。回归分析研究的是随机变量与普通变量的统计关系，如体 重与身高的关系等。一插值方法定义已知某未知函数y = f (x)的一组观测(或试验)数据 (Xi,yi)(i=1,2n),要寻找一个函数 欠x),使外x) = V = f (x),则称 此类问题为插值问

3、题。并称Wx)为f(x)的插值函数，并称Xi,X2，Xn 为样本点；称(x) = X为插值条件，可得：(x)*f(x)常见的插值方法1、线性插值，linear；2、拉格朗日插值方法；lagri3、牛顿插值方法；newton4、最近点插值方法，nearst5、cubic 插值方法，cubic 或 pchip6、三次分段样条插值，spline7、B样条函数插值方法(k次等距B样条函数插值函数方 法)8、人工神经网络插值方法9、分形插值方法详细算法参考数学建模相关书籍，本文只讲给定样本点， 如何利用数学软件实现各类插值的计算。一维插值方法的 matlab软件实现一般格式形如：y = interp1(

4、%,y0,x,lagrT),其中A4为样本数据序列，X是待插值点。lagr1是插值方 法选项，上述7种插值方法都可以选择。 方法缺省，表示线性 插值。例1.1利用y=sin(x),在0,2冗上取10个点作为样本点，在取 12个点作为插值点，利用最近点插值和线性插值分别计算， 结果如下图显示。三次样条插值三次样条插值原理见高等应用数学的Matlab求解P286二次分段多项式：s（x） =Gi（x-为）3 +q2（x x ）2 +63（乂一为）+q4以y =sin（x）为例x0=linspace(0,2*pi,10);y0=sin(x0);sp=csapi(x0,y0);%sp由(x0,y0)确定

5、具体的样条函数，主要确定插值系数x=linspace(0,2*pi,20);% 待插值点y=fnval(sp,x);%三次样条插值B样条插值以y =sin(x)为例x0=linspace(0,2*pi,10);y0=sin(x0);sp=spapi(k,x0,y0);%sp由(x0,y0)确定具体的样条函数，k表示 样条阶次，一般 4, 5 x=linspace(0,2*pi,20);% 待插值点 y=fnval(sp,x);%三次样条插值二维插值问题观测数据格式如下：1小乙”*2 丫2 Z2 : : : sin(x.*y)0.84150.90930.9093-0.75680.1411-0.2

6、794三维网格点数据：复杂些，以三维矩阵形式给由O假设自变量x取值1, 2, y取值1, 2, 3, z取值为1, 2, 3, 4,则自变 量在空间上形成2M3=24个点，这24个点以三维矩阵形式给 出,x网格点数据形式：x(:,:,1), x(:,:,2), x(:,:,3), x(:,:,4),并 且x(:,:,i)均是相同的平面网格点数据，y,z也依此类似，其中 x(:,:,1)2 TOC o 1-5 h z 2 y(：,：,i)ii23 z(:,:,1)111111 v=sin(xx.*yy.*zz)v(：,：,1)=0.84150.90930.9093-0.75680.1411-0.

7、2794v(:,:,2), v(:,:,3), v(:,:,4)与 v(:,:,1)类似，没有列由。高维插值计算 1: vv=interp3(x,y,z,v,a,b,c, spline)其中，a,b,c是插值点，可以是列向量或行向量，vv也是对应的列（行）向量。插值方法还有:nearest - nearest neighbor interpolation linear - linear interpolationspline - spline interpolationcubic - cubic interpolation as long as the data is uniformlyspa

8、ced, otherwise the same as spline高维插值计算 2: vv=griddata3(x,y,z,v,a,b,c,linear)其中插值方法还可选择nearest；其它方法不行。1.8普通数据的高维插值问题V=griddatan(x0,v0,x);其中x0, v0构成样本点，x0是自变量，每列表示一个自 变量，v0是因变量，是列向量。类似I xl1xlnv V1x2ix2n 丫2x0 = ： : :, v0 =: xm1xmn fkvm /同时，插值点x也必须与x0有相同结构。其中插值方法还可 选择nearest,其它方法不行。1.9地貌插值问题已知某处山区地形选点测

9、量坐标数据为：x=0 0,5 1 1.5 2 2.5 3 3.5 4 4.5 5y=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 海拔高度数据为： TOC o 1-5 h z z=89 90 878592 91 96 93 90878292 96 989995 91 89 86 84828496 98 959290 88 85 84 83818580 81 828995 96 93 92 89868682 85 879899 96 97 88 85828382 85 899495 93 92 91 86848888 92 939495 89 87 86 83819

10、292 96 979896 93 95 84 82818485 85 818280 80 81 85 90939584 86 819899 98 97 96 95848780 81 858283 84 87 90 95868880 82 818485 86 83 82 81808287 88 899899 97 96 98 949287其地貌图为:对数据插值加密形成地貌图。程序：x=0:.5:5;y=0:.5:6;z=89 90 87 85 92 91 96 93 90 87 8292 96 98 99 95 91 89 86 84 82 8496 98 95 92 90 88 85 84 8

11、3 81 8580 81 82 89 95 96 93 92 89 86 8682 85 87 98 99 96 97 88 85 82 8382 85 89 94 95 93 92 91 86 84 8888 92 93 94 95 89 87 86 83 81 9292 96 97 98 96 93 95 84 82 81 8485 85 81 82 80 80 81 85 90 93 95 TOC o 1-5 h z 84 8681 98 99 98 979695848780 8185 82 83 84 879095868880 8281 84 85 86 838281808287 88

12、89 98 99 97 9698949287;mesh(x,y,z) 腔原始数据图xi=linspace(0,5,50);%加密横坐标数据到 50个yi=linspace(0,6,80);呦口密纵坐标数据到 60个xii,yii=meshgrid(xi,yi);%& 成网格数据zii=interp2(x,y,z,xii,yii,cubic); %插值mesh(xii,yii,zii)吸口密后的地貌图hold on %保持图形xx,yy=meshgrid(x,y);%成网格数据plot3(xx,yy,z+0.1, ob )% 原始数据用O绘出File Edit Tools Windw Help口

13、令。昌|，A?/彦-JDI x| 0101.20二元非等距插值调用格式：zi=griddata(x,y,z,xi,yi, 指定插值方法) 插值方法有： linear %线性插值(默认)bilinear%双线性插值cubic%三次插值bicubic%双三次插值nearest%最近邻域插值例：用随机数据生成地貌图再进行插值 程序： x=rand(100,1)*4-2; y=rand(100,1)*4-2;z=x.*exp(-x.A2-y.A2);ti=-2:.25:2;xi,yi=meshgrid(ti,ti);% 加密数据zi=griddata(x,y,z,xi,yi);% 线性插值mesh(x

14、i,yi,zi) hold on plot3(x,y,z,o)11该例中使用的数据是随机形成的，故函数griddata 可以处理无规则的数据。二数据拟合方法数据拟合 方法不要求函数（曲线或曲面）经过所有已知数 据点，而是要求它能较好的反映数据的整体变化趋势，称解决这类问题的方法为数据拟合。常用方法：最小二乘拟合方法，分为线性最小二乘拟合方 法，非线性最小二乘拟合方法。 思想：数据点的误差平方和最 小。多项式数据拟合数学原理：无穷级数，泰勒定理12多项式函数形式：y = PnXn + PnXnJL+ PiX1 + Po , p= Pn, Pn,，R, P0Syntax:P = Polyfit(x

15、,y,n)数据拟合P,S = Polyfit(x,y,n)生成多项式函数：f=Poly2sym(P)多项式函数插值点计算：polyval(p,xi) , p是多项式系数。一般非线性曲线拟合拟合 Cobb-Douglas生产函数：Q(K,L) =aK： L： 0 :二一：1式中a,% P要由经济统计数据确定。现有美国马萨诸塞州19001926年上述三个经济指数的统计数据，如下表，试用数据拟合的方法，求生(*)式中的参数民,U),a。TQKLTQKL19001.051.041.0519131.952.821.6819011.181.061.0819142.013.241.6519021.291.1

16、61.1819152.003.241.6219031.301.221.2219162.093.611.8619041.301.271.1719171.964.101.931319051.421.371.3019182.204.361.9619061.501.441.3919192.124.771.9519071.521.531.4719202.164.751.9019081.461.571.3719212.084.541.5819091.602.051.4319222.244.541.6719101.692.511.5819232.564.581.8219111.812.631.5919242

17、.344.581.6019121.932.741.6619252.454.581.61X/19262.584.541.64Matlab软件计算：f=(c,x)c(1)*x(:,1),Ac(2).*x(:,2),Ac(3);y=1.05 1.18 1.29 1.30 1.30 1.42 1.50 1.52 1.46 1.60 1.691.81 1.93 1.95 2.01 2.00 2.09 1.96 2.20 2.12 2.16 2.08 2.242.56 2.34 2.45 2.58;x=1.04 1.06 1.16 1.22 1.27 1.37 1.44 1.53 1.57 2.05 2.

18、512.63 2.74 2.82 3.24 3.24 3.61 4.10 4.36 4.77 4.75 4.54 4.544.58 4.58 4.58 4.54;1.05 1.08 1.18 1.22 1.17 1.30 1.39 1.47 1.31 1.43 1.58 1.591.66 1.68 1.65 1.62 1.86 1.93 1.96 1.95 1.90 1.58 1.67 1.82141.60 1.61 1.64;c0=0,1,0.1,0.2;c=nlinfit(x,y,f,c0)运行结果： c=1.2246 0.4612 -0.1277 eviews软件计算：(1)输入x1 ,

19、 x2, y序列；(2)在 equation 窗口里输入：y=c(1)*x1Ac(2)*x2Ac(3)(3)运行结果如下：Dependent VariablesMethod: Least SquaresDate： 05/23/09 Tima： 10:27Sample: 1 27Incljded observations： 27Convergence achieved after 7 iterationsY=C(1AC(2)|t|95常 conf. interval/cl1.224636.057060621.460.0001.1068691.342403/C2.461A56W,05521778.

20、50.000.3471926.57512心-.1276689.1725373-0.740.467-4酩7684.2284306说明：简单问题，上述三个都可以，但是对于复杂问题，stata功能更全、更强大。三回归分析回归模型的计算与曲线拟合没有区别，只是解释不同罢163.1多元线性回归 2、 y =汾 P1X1 P2X2PnXn； N(0, r )Syntaxb = regress(y,X)b,bint,r,rint,stats = regress(y,X)b,bint,r,rint,stats = regress(y,X,alpha)yJ y=yn)(1)原始数据 TOC o 1-5 h z

21、11 X11X1m：X = 1 1Xn1Xnm，其中解释变量第1列的系数是P0O(2) p=regress(y,x)3.2非线性回归非线性模型：y = G , (0.49 - g )* ef(x,) 观察数据如下：17XyXyXyXy80.49140.43240.42320.480.49160.44240.4340.4100.48160.43240.4360.41100.47160.43260.41360.38100.48180.46260.4380.4100.47180.45260.41380.4120.46200.42280.41400.39120.46200.42280.4420.391

22、20.45200.43300.4120.43220.41300.4140.45220.41300.38140.43220.4320.41使用nlinfit命令f=(c,x)c(1)+(0.49-c(1)*exp(-c(2)*(x(:,1)-8);c0=0.1,0.01;%赋参数初值c=nlinfit(x(:,1),x(:,2),f,c0)运行结果：c = 0.3901400079615900.101632683988975使用lsqcurvefit命令f=(c,x)c(1)+(0.49-c(1)*exp(-c(2)*(x(:,1)-8);c=lsqcurvefit(f, 0.1,0.01,x(

23、:,1),x(:,2)运行结果：c= 0.3901400079615900.101632683988975注意：lsqcurve与nlinfit结合使用，lsqcurve给由结果不确定，与初值有很大关系，往往给由满足默认精度要求的结果； nlinfit对初值要求高，如果初值与真值偏离很大，会提示不能运行。实际中，可以把 lsqcurve 给由结果作为nlinfit 的 初值。还可以用eviews软件和stata软件计算。18例2007赛题：中国人口预测生育率(r岁育龄妇女千人中生育的比例)曲线可用如下模型 刻画：(r -15)：e一三,二()为了拟合由模型中参数，则对函数进行简化：hr =G(

24、r -15)气(,5) r =1549Matlab 程序：f=(c,x)c(1)*(x-15)Ac(2).*exp(-c(3)*(x-15);x=15:49;y=0.02,0.88,1.16,2.83,15.23,33.05,97.53,161.11,200.7,176.2,179.37,148.09,112.65,97.64,79.58,59.66,59.65,48.57,35.18,27.23,16.83,13.67,7.21,6.51,5.68,3.89,5.97,2.01,0.29,1.5,1.37,0.41,0.45,1.23,0.64;c=nlinfit(x,y,f,0.0211 7.5051 0.8248)y2=c(