《新编MATLAB自学一本通》课件第20章 多项式回归与数据插值

2025/4/30主要内容多项式回归

插值问题的数学描述一维插值二维插值高维插值2025/4/30第一节多项式回归2025/4/30一、多项式回归模型2025/4/30二、多项式回归的MATLAB实现[p,S,mu]=polyfit(x,y,n)自变量观测值向量因变量观测值向量多项式阶数1.polyfit函数的用法系数向量的估计值用于误差估计的结构体变量均值和标准差2025/4/30[y,delta]=polyval(p,x,S)用户指定的自变量取值向量因变量估计值向量2.polyval函数的用法多项式系数向量用于误差估计的结构体变量误差标准差的估计值向量2025/4/30r=poly2sym(p,v)多项式的符号表达式3.poly2sym函数的用法多项式系数向量自变量符号或取值2025/4/30三、多项式回归案例【例20.1-1】现有我国2007年1月至2011年11月的食品零售价格分类指数数据，如表20.1-1所示。数据来源：中华人民共和国国家统计局网站月度统计数据。试根据这59组统计数据研究全国食品零售价格分类指数y（上年同月=100）和时间x

之间的关系。序号统计月度上年同月=100上年同期=100全国城市农村全国城市农村12007年1月104.9104.4105.9104.9104.4105.922007年2月105.8105.2106.9105.3104.8106.432007年3月107.7107.4108.3106.1105.7107……………………572011年9月113.5113.4114112.6112.4113.1582011年10月111.9111.8112.2112.5112.3113592011年11月108.7108.8108.6112.2112112.62025/4/301.数据的散点图2025/4/302.四次多项式拟合假设y关于x的理论回归方程为调用polyfit函数求系数估计值>>p4=polyfit(x,y,4)p4=-0.00010.0096-0.39855.563594.2769>>r=poly2sym(p4);>>r=vpa(r,5)

r=-0.000074268*x^4+0.0096077*x^3-0.39845*x^2+5.5635*x+94.2772025/4/303.更高次多项式拟合调用polyfit函数作更高次（大于4次）多项式拟合，并把多次拟合的残差的模加以对比，评价拟合的好坏。>>[p5,S5]=polyfit(x,y,5);%5次多项式拟合>>S5.normr>>[p6,S6]=polyfit(x,y,6);%6次多项式拟合>>S6.normr>>[p7,S7]=polyfit(x,y,7);%7次多项式拟合>>S7.normr>>[p8,S8]=polyfit(x,y,8);%8次多项式拟合>>S8.normr>>[p9,S9]=polyfit(x,y,9);%9次多项式拟合>>S9.normr2025/4/304.拟合效果图在以上拟合结果的基础上，可以调用polyval函数计算给定自变量处的因变量的预测值，从而绘制拟合效果图，从拟合效果图上直观地看出拟合的准确性。代码略。。。2025/4/302025/4/30第二节插值问题的数学描述2025/4/30一、什么是插值所谓插值就是在已知离散数据的基础上补插连续函数，使得这条连续曲线（或曲面）通过全部已知的离散数据点，利用插值方法可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。

一维插值问题的数学描述为：

已知某一函数

y=g(x)（g(x)的解析表达式可能十分复杂,也可以是未知的）在区间[a,

b]上n+1个互异点xj处的函数值yj，j=0,1,…,n

，还知道g(x)在[a,

b]上有若干阶导数，如何求出g(x)在[a,

b]上任一点x的近似值。二、一维插值问题的数学描述１．数学描述

一维插值方法的基本思想是：根据g(x)在区间[a,b]上n+1个互异点

xj（称为节点）的函数值

yj，j=0,1,…,n，求一个足够光滑、简单便于计算的函数f(x)（称为插值函数）作为g(x)的近似表达式，使得f(xj

)=yj，j=0，1，…，n。（1）然后计算f(x)在区间[a,b](称为插值区间)上点x（称为插值点）的值作为原函数g(x)（称为被插函数）在此点的近似值。求插值函数f(x)的方法称为插值方法，（1）式称为插值条件。代数多项式比较简单，常用多项式作为插值函数。2．基本思想常用的一维插值方法有：分段线性插值、拉格朗日（Lagrange）多项式插值、牛顿（Newton）插值、Hermite插值、最近邻插值、三次样条插值和B样条插值等。3．常用一维插值方法

二维插值问题的数学描述为：

已知某二元函数

z=G(x,y)（解析表达式可能十分复杂,也可以是未知的）在平面区域D上N个互异点(xi,yi)处的函数值zi，i=0,1,…,N

，求一个足够光滑、简单便于计算的插值函数f(x,y)。由插值函数可以计算原函数在平面区域上任意点处的近似值。常用的二维插值方法有：分片线性插值、双线性插值、最近邻插值、三次样条插值和B样条插值等。三、二维插值问题的数学描述

是一个细的、可弯曲的木制或塑料条，在飞机或轮船等的设计制造过程中为描绘出光滑的外形曲线(放样)所用的工具

从物理上讲，样条满足插值点的约束，同时使势能达到最小

三次样条本质上是一段一段的三次多项式拼合而成的曲线，在拼接处,不仅函数是连续的,且一阶和二阶导数也是连续的，1946年,Schoenberg将样条引入数学,即所谓的样条函数1.什么是样条？四、三次样条插值的数学描述2.三次样条插值函数3.三次样条插值原理在n个小区间构造S(x)，共有n个三次多项式，需确定4n个参数在所有节点上n+1个方程在除端点外的节点上3(n-1)个方程这样就得到4n-2个方程，为保证待定参数的唯一性，还差两个方程。为此，常用的方法是对边界节点除函数值外附加要求，这就是所谓的边界条件。根据实际问题的不同，三次样条插值常用到下列三类边界条件。

周期边界条件：

当y=g(x)是以

b-a=x0-xn

为周期的周期函数时，要求S(x)也是周期函数，故端点处要满足此条件称为周期条件。

m边界条件：

即给定端点处的一阶导数值。

M边界条件：

即给定端点处的二阶导数值。2025/4/30第三节一维插值一、Lagrange（拉格朗日）插值1．拉格朗日线性插值插值多项式：直线的两点式表达式分别称为节点x0和x1的一次插值基函数插值函数为基函数的线性组合，组合系数就是对应节点上的函数值在上，过点2.拉格朗日二次插值二次插值多项式：对于，可假定：基函数满足：基函数如何确定？从而过点

构造一组插值基函数（n次多项式）缺点：当n

比较大时，插值多项式Ln(x)的收敛性与稳定性变

差，逼近效果不理想。

li(x)是n次多项式，且满足

由插值基函数li(x)构造n次拉格朗日插值多项式如下：3.拉格朗日n次插值

拉格朗日n次插值的误差估计拉格朗日插值产生的截断误差为Rn(x)，则4.自编拉格朗日插值函数lagrangefunctiony=lagrange(x0,y0,x)ifnumel(x0)~=numel(y0)error('原始坐标点应等长')endx0=x0(:);ifany(diff(sort(x0))==0)error('插值节点不能有重复值')endn=numel(x0);m=numel(x);y0=y0(:);x=x(:)';y=zeros(n,m);fori=1:ny(i,:)=y0(i)*prod(repmat(x,n-1,1)-...repmat(x0(1:n~=i),1,m))/prod(x0(i)-x0(1:n~=i));endy=sum(y);2025/4/30【例20.3-1】Runge（龙格）现象：用函数在区间[-1,1]上产生11个等距节点，然后调用自编lagrange函数作拉格朗日插值。>>x0=linspace(-1,1,11);>>y0=1./(1+25*x0.^2);>>x=linspace(-1,1,100);>>f=1./(1+25*x.^2);>>y=lagrange(x0,y0,x);>>plot(x0,y0,'ko');>>holdon;>>plot(x,f,'k','linewidth',2);>>plot(x,y,'k:','linewidth',2);>>xlabel('X');>>ylabel('$$f(x)=\frac{1}{1+25x^2}$$','Interpreter','latex');>>legend('插值节点','原函数图像','Lagrange插值')2025/4/30二、interp1函数yi=interp1(x,Y,xi,method)自变量观测值向量因变量观测值向量用户指定的插值点横坐标1.interp1函数的用法计算得到的近似函数值指定所用的插值方法2025/4/30【例20.3-1续】针对例20.3-1中的数据，调用interp1函数作一维插值。>>ylin=interp1(x0,y0,x);>>yspl=interp1(x0,y0,x,'spline');>>plot(x0,y0,'ko');>>holdon;>>plot(x,f,'k','linewidth',2);>>plot(x,ylin,':','linewidth',2);>>plot(x,yspl,'r-.','linewidth',2);>>xlabel('X');>>ylabel('$$f(x)=\frac{1}{1+25x^2}$$','Interpreter','latex');>>legend('插值节点','原函数图像','分段线性插值','三次样条插值')2025/4/30【例20.3-2】在加工机翼的过程中，已有机翼断面轮廓线上的20组坐标点数据，如表20.3-2所列，其中(x,y1)和(x,y2)分别对应轮廓线的上下线。假设需要得到x坐标每改变0.1时的y坐标，试通过插值方法计算加工所需的全部数据，并绘制机翼断面轮廓线，求加工断面的面积。03579111213141501.82.22.73.03.12.92.52.01.601.21.72.02.12.01.81.21.01.62025/4/30x0=[0,3,5,7,9,11,12,13,14,15];y01=[0,1.8,2.2,2.7,3.0,3.1,2.9,2.5,2.0,1.6];y02=[0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6];x=0:0.1:15;ysp1=interp1(x0,y01,x,'spline');pp=interp1(x0,y02,'spline','pp');ysp2=ppval(pp,x);xx=[x,fliplr(x)];ysp=[ysp1,fliplr(ysp2)];plot([x0,x0],[y01,y02],'o');holdon;plot(xx,ysp,'r','linewidth',2);xlabel('X');ylabel('Y');legend('插值节点','三次样条插值','location','northwest');%截面面积：S1=trapz(x,ysp1)-trapz(x,ysp2)S2=trapz(xx,ysp)2025/4/30三、spline函数yy=spline(x,Y,xx)自变量观测值向量因变量观测值向量用户指定的插值点横坐标1.spline函数的用法计算得到的近似函数值2025/4/30四、csape和csapi函数pp=csape(x,y,conds)自变量观测值向量因变量观测值向量边界条件参数1.csape函数的用法分段多项式形式的插值结果conds参数取值说

明'complete'或'clamped'给定端点处的一阶导数值，默认为拉格朗日边界条件'not-a-knot'非纽结边界条件，csapi函数使用的就是这种边界条件'periodic'周期边界条件'second'给定端点处的二阶导数值。默认为[0,0]，同'variational'情形'variational'设定端点处的二阶导数值为02025/4/30values=csapi(x,y,xx)2.csapi函数的用法自变量观测值向量因变量观测值向量用户指定的插值点横坐标计算得到的近似函数值2025/4/302025/4/30fun=@(x)sin(pi*x/2).*(x>=-1&x<1)+x.*exp(1-x.^2).*(x>=1|x<-1);%%----------------区间[0,1]上的三次样条插值------------------x01=linspace(0,1,6);y01=fun(x01);x1=linspace(0,1,20);pp1=csape(x01,[1,y01,0],'complete');y1=fnval(pp1,x1);%%----------------区间[1,3]上的三次样条插值------------------x02=linspace(1,3,8);y02=fun(x02);x2=linspace(1,3,30);pp2=csape(x02,[0,y02,0.01],[1,2]);y2=fnval(pp2,x2);%%-----------------------绘图---------------------plot([x01,x02],[y01,y02],'ko');holdon;plot([x1,x2],fun([x1,x2]),'k','linewidth',2);plot([x1,x2],[y1,y2],'--','linewidth',2);xlabel('X');ylabel('Y=f(x)');legend('插值节点','原函数图像','三次样条插值');2025/4/30fun=@(x)sin(pi*x/2).*(x>=-1&x<1)+x.*exp(1-x.^2).*(x>=1|x<-1);%%----------------区间[0,3]上的三次样条插值------------------x0=[linspace(0,1,6),linspace(1.1,3,8)];y0=fun(x0);x=linspace(0,3,61);y=csapi(x0,y0,x);%%-----------------------绘图---------------------plot(x0,y0,'ko');holdon;plot(x,fun(x),'k','linewidth',2);plot(x,y,'--','linewidth',2);xlabel('X');ylabel('Y=f(x)');legend('插值节点','原函数图像','三次样条插值');2025/4/30五、spapi函数（B样条插值）sp=spapi(k,x,y)自变量观测值向量因变量观测值向量控制节点坐标或B样条阶数1.spapi函数的用法B样条函数2025/4/302025/4/30fun=@(x)sin(pi*x/2).*(x>=-1&x<1)+x.*exp(1-x.^2).*(x>=1|x<-1);%%----------------区间[0,3]上的三次B样条插值------------------x0=[linspace(0,1,6),linspace(1.1,3,8)];y0=fun(x0);x=linspace(0,3,61);sp3=spapi(4,x0,y0);sp4=spapi(8,x0,y0);%%-----------------------绘图---------------------plot(x0,y0,'ko');holdon;plot(x,fun(x),'k','linewidth',2);fnplt(sp3,2,'--');fnplt(sp4,2,'r:');xlabel('X');ylabel('Y=f(x)');legend('插值节点','原函数图像','三次B样条插值','七次B样条插值');2025/4/30六、其他一维插值函数

MATLAB样条插值工具箱中还提供了其他一些可以作一维插值的MATLAB函数，例如csaps、spaps、spap2和cscvn等函数，其中csaps、spaps和spap2函数还可以作二维和高维插值。函数名调用格式功能及参数说明csapspp=csaps(x,y,p)values=csaps(x,y,p,xx)三次光滑样条插值spapssp=spaps(x,y,tol)[sp,values]=spaps(x,y,tol)光滑B样条插值spap2sp=spap2(knots,k,x,y)最小二乘B样条近似cscvncurve=cscvn(points)具有自然边界条件或周期边界条件的三次样条插值2025/4/30【例20.3-5】产生区间[0,2p]上的带有噪声的正弦函数值，然后分别调用csaps、spaps和spap2函数作样条插值，绘制插值效果图。>>x0=linspace(0,2*pi,15);>>y0=sin(x0)+0.3*(rand(size(x0))-0.5);>>pp=csaps(x0,y0,0.9);>>sp1=spaps(x0,y0,1e-3);>>sp2=spap2(3,4,x0,y0);>>plot(x0,y0,'ko');>>holdon>>fnplt(pp,2,'r:')>>fnplt(sp1,2,'k-.')>>fnplt(sp2,2,'--')>>xlabel('X');ylabel('Y=sin(x)');>>legend('插值节点','三次光滑样条插值','光滑B样条插值','最小二乘B样条近似');2025/4/30【例20.3-6】MATLAB创始者CleveMoler博士写的《MATLAB数值计算》一书中有一个有趣的例子，就是利用样条插值方法描绘手的轮廓线。思路就是首先运行[x,y]=ginput命令（ginput为取点函数），然后把手放在屏幕前，用鼠标沿着手的轮廓点击取点，不用点太密，轮廓变化大的地方（例如指尖和两指相连处）可多取几个点。点击完毕后按Enter键退出取点模式，MATLAB会记录下这些点的坐标。对这些点调用cscvn函数进行顺序节点插值，最后就可画出手的轮廓。2025/4/30>>figure('position',get(0,'screensize'));>>axes('position',[0011]);>>[x,y]=ginput;>>curve=cscvn([x';y']);>>plot(x,y,'ko');>>holdon;>>fnplt(curve);>>xlabel('X');>>ylabel('笔者的手');2025/4/30第四节二维插值2025/4/30一、网格节点插值1.网格节点插值图解2025/4/30前面介绍的csape、csapi、spapi、csaps、spaps和spap2函数均可以作二维和高维网格节点插值，除此之外，MATLAB多项式函数工具箱中还提供了interp2函数，用来作二维网格节点插值。2.网格节点插值的MATLAB函数2025/4/30

csape、csapi、spapi、csaps、spaps和spap2函数用法函数名调用格式功能及参数说明csapepp=csape({x1,x2},Y)pp=csape({x1,x2},Y,{conds1,conds2})使用指定边界条件作三次样条插值csapipp=csapi({x1,x2},Y)values=csapi({x1,x2},Y,{xx1,xx2})使用非纽结边界条件作三次样条插值spapisp=spapi({knork1,knork2},{x1,x2},Y)指定阶次的B样条插值csapspp=csaps({x1,x2},Y,{p1,p2})values=csaps({x1,x2},Y,{p1,p2},{xx1,xx2})三次光滑样条插值spapssp=spaps({x1,x2},Y,{tol1,tol2})[sp,values]=spaps({x1,x2},Y,{tol1,tol2})光滑B样条插值spap2sp=spap2({knork1,knork2},k,{x1,x2},Y)最小二乘B样条近似2025/4/30（1）ZI=interp2(X,Y,Z,XI,YI)（2）ZI=interp2(Z,XI,YI)（3）ZI=interp2(Z,ntimes)（4）ZI=interp2(X,Y,Z,XI,YI,method)节点Y坐标节点X坐标插值点处近似函数值

interp2函数用法节点Z坐标插值点X坐标插值点Y坐标指定所用的插值方法2025/4/30【例20.4-1】在一丘陵地带测量高程，x和y方向每隔100米测一个点，得高程数据如表20.4-3所列，试用插值方法拟合出地形曲面。3.网格节点插值举例高程1002003004005001004504786246976362004204786307126983004004125986746804003103345526266622025/4/30x=100:100:500;y=100:100:400;[X,Y]=meshgrid(x,y);Z=[450478624697636420478630712698400412598674680310334552626662];xd=100:20:500;yd=100:20:400;[Xd1,Yd1]=meshgrid(xd,yd);[Xd2,Yd2]=ndgrid(xd,yd);figure;%新建图形窗口%--------------调用interp2函数作三次样条插值-------------------Zd1=interp2(X,Y,Z,Xd1,Yd1,'spline');subplot(2,3,1);surf(Xd1,Yd1,Zd1);xlabel('X');ylabel('Y');zlabel('Z');title('interp2')2025/4/30二、散乱节点插值1.散乱节点插值图解2025/4/30

MATLAB多项式函数工具箱中提供了TriScatteredInterp和griddata函数，其中前者用来作二维或三维散乱节点插值，后者只能用来作二维散乱节点插值，前者比后者的运行效率高。2.散乱节点插值的MATLAB函数函数名调用格式功能及参数说明griddataZI=griddata(x,y,z,XI,YI)[XI,YI,ZI]=griddata(x,y,z,XI,YI)[...]=griddata(...,method)返回插值点XI,YI处的近似函数值矩阵ZI。输入参数x,y是节点坐标向量，z是相应的原函数值向量。XI，YI是用户指定的插值点坐标，可以是同型矩阵，也可以是向量（XI为行向量，YI为列向量）。method参数用来指定所用的插值方法，其可用取值为：'linear'线性插值（默认）'cubic'立方插值'nearest'最近邻插值'v4'MATLAB4中用到的插值方法TriScatteredInterpF=TriScatteredInterp(X,V)返回TriScatteredInterp类变量F。输入参数X是m行n列的矩阵，其中m为节点个数，n为2或3，表示维数，V是节点处原函数值，是长度为m的列向量。此时由VI=F(XI)计算插值点XI处的近似函数值VI。F=TriScatteredInterp(X,Y,V)X,Y,V是等长的列向量，X,Y为节点坐标，V为相应的原函数值。此时由VI=F(XI,YI)计算插值点XI,YI处的近似函数值VI。F=TriScatteredInterp(X,Y,Z,V)X,Y,Z,V是等长的列向量，X,Y,Z为节点坐标，V为相应的原函数值。此时由VI=F(XI,YI,ZI)计算插值点XI,YI,ZI处的近似函数值VI。F=TriScatteredInterp(...,method)用method参数指定所用的插值方法，其可用取值为：'natural'自然近邻插值'linear'线性插值（默认）'nearest'最近邻插值2025/4/30【例20.4-2】2011高教社杯全国大学生数学建模竞赛A题中给出了某城市城区土壤地质环境调查数据，包括采样点的位置、海拔高度及其所属功能区等信息数据，以及8种主要重金属元素在采样点处的浓度、8种主要重金属元素的背景值数据。试根据调查数据中给出的采样点坐标和8种主要重金属元素的浓度数据绘制Cd元素的空间分布图。3.散乱节点插值举例2025/4/30xyz=xlsread('cumcm2011A.xls',1,'B4:D322');Cd=xlsread('cumcm2011A.xls',2,'C4:C322');x=xyz(:,1);y=xyz(:,2);z=xyz(:,3);xd=linspace(min(x),max(x),60);yd=linspace(min(y),max(y),60);[Xd,Yd]=meshgrid(xd,yd);%调用griddata