Matlab插值与拟合教程

MATLAB插值与拟合§1曲线拟合实例：温度曲线问题气象部门观测到一天某些时刻的温度变化数据为：t012345678910T1315171416192624262729试描绘出温度变化曲线。曲线拟合就是计算出两组数据之间的一种函数关系，由此可描绘其变化曲线及估计非采集数据对应的变量信息。曲线拟合有多种方式，下面是一元函数采用最小二乘法对给定数据进行多项式曲线拟合，最后给出拟合的多项式系数。1.1.线性拟合函数：regress()调用格式：b=regress(y,X)[b,bint,r,rint,stats]=regress(y,X)[b,bint,r,rint,stats]=regress(y,X,alpha)说明：b=regress(y,X)返回X处y的最小二乘拟合值。该函数求解线性模型：y=Xβ+εβ是p1的参数向量；ε是服从标准正态分布的随机干扰的n1的向量；y为n1的向量；X为np矩阵。bint返回β的95%的置信区间。r中为形状残差，rint中返回每一个残差的95%置信区间。Stats向量包含R2统计量、回归的F值和p值。例1：设y的值为给定的x的线性函数加服从标准正态分布的随机干扰值得到。即y=10+x+ε；求线性拟合方程系数。程序：x=[ones(10,1)(1:10)’]y=x*[10;1]+normrnd(0,0.1,10,1)[b,bint]=regress(y,x,0.05)结果：x=111213141516171819110y=10.956711.833413.012514.028814.885416.119117.118917.996219.032720.0175b=9.92131.0143bint=9.788910.05370.99301.0357即回归方程为：y=9.9213+1.0143x2.2.多项式曲线拟合函数：polyfit()调用格式：p=polyfit(x,y,n)[p,s]=polyfit(x,y,n)说明：x,y为数据点，n为多项式阶数，返回p为幂次从高到低的多项式系数向量p。矩阵s用于生成预测值的误差估计。(见下一函数polyval)例2：由离散数据x0.1.2.3.4.5.6.7.8.91y.3.511.41.61.9.6.4.81.52拟合出多项式。程序：x=0:.1:1;y=[.3.511.41.61.9.6.4.81.52]n=3;p=polyfit(x,y,n)xi=linspace(0,1,100);z=polyval(p,xi);%多项式求值plot(x,y,’o’,xi,z,’k:’,x,y,’b’)legend(‘原始数据’,’3阶曲线’)结果：p=16.7832-25.745910.9802-0.0035多项式为：16.7832x3-25.7459x2+10.9802x-0.0035曲线拟合图形：也可由函数给出数据。例3：x=1:20,y=x+3*sin(x)程序：x=1:20;y=x+3*sin(x);p=polyfit(x,y,6)xi=1inspace(1,20,100);z=poyval(p,xi);%多项式求值函数plot(x,y,’o’,xi,z,’k:’,x,y,’b’)legend(‘原始数据’,’6阶曲线’)结果：p=0.0000-0.00210.0505-0.59713.6472-9.729511.3304再用10阶多项式拟合程序：x=1:20;y=x+3*sin(x);p=polyfit(x,y,10)xi=linspace(1,20,100);z=polyval(p,xi);plot(x,y,'o',xi,z,'k:',x,y,'b')legend('原始数据','10阶多项式')结果：p=Columns1through70.0000-0.00000.0004-0.01140.1814-1.806511.2360Columns8through11-42.086188.5907-92.815540.2671可用不同阶的多项式来拟合数据，但也不是阶数越高拟合的越好。3.3.多项式曲线求值函数：polyval()调用格式：y=polyval(p,x)[y,DELTA]=polyval(p,x,s)说明：y=polyval(p,x)为返回对应自变量x在给定系数P的多项式的值。[y,DELTA]=polyval(p,x,s)使用polyfit函数的选项输出s得出误差估计YDELTA。它假设polyfit函数数据输入的误差是独立正态的，并且方差为常数。则YDELTA将至少包含50%的预测值。4.4.多项式曲线拟合的评价和置信区间函数：polyconf()调用格式：[Y,DELTA]=polyconf(p,x,s)[Y,DELTA]=polyconf(p,x,s,alpha)说明：[Y,DELTA]=polyconf(p,x,s)使用polyfit函数的选项输出s给出Y的95%置信区间YDELTA。它假设polyfit函数数据输入的误差是独立正态的，并且方差为常数。1-alpha为置信度。例4：给出上面例1的预测值及置信度为90%的置信区间。程序：x=0:.1:1;y=[.3.511.41.61.9.6.4.81.52]n=3;[p,s]=polyfit(x,y,n)alpha=0.05;[Y,DELTA]=polyconf(p,x,s,alpha)结果：p=16.7832-25.745910.9802-0.0035s=R:[4x4double]df:7normr:1.1406Y=Columns1through7-0.00350.85381.29701.42661.34341.14800.9413Columns8through110.82380.89631.25942.0140DELTA=Columns1through71.36391.15631.15631.15891.13521.12021.1352Columns8through111.15891.15631.15631.36395.5.稳健回归函数：robust()稳健回归是指此回归方法相对于其他回归方法而言，受异常值的影响较小。调用格式：b=robustfit(x,y)[b,stats]=robustfit(x,y)[b,stats]=robustfit(x,y,’wfun’,tune,’const’)说明：b返回系数估计向量；stats返回各种参数估计；’wfun’指定一个加权函数；tune为调协常数；’const’的值为’on’(默认值)时添加一个常数项；为’off’时忽略常数项。例5：演示一个异常数据点如何影响最小二乘拟合值与稳健拟合。首先利用函数y=10-2x加上一些随机干扰的项生成数据集，然后改变一个y的值形成异常值。调用不同的拟合函数，通过图形观查影响程度。程序：x=(1:10)’;y=10-2*x+randn(10,1);y(10)=0;bls=regress(y,[ones(10,1)x])%线性拟合brob=robustfit(x,y)%稳健拟合scatter(x,y)holdonplot(x,bls(1)+bls(2)*x,’:’)plot(x,brob(1)+brob(2)*x,’r‘)结果：bls=8.4452-1.4784brob=10.2934-2.0006分析：稳健拟合(实线)对数据的拟合程度好些，忽略了异常值。最小二乘拟合(点线)则受到异常值的影响，向异常值偏移。6.6.向自定义函数拟合对于给定的数据，根据经验拟合为带有待定常数的自定义函数。所用函数：nlinfit()调用格式：[beta,r,J]=nlinfit(X,y,’fun’,betao)说明：beta返回函数’fun’中的待定常数；r表示残差；J表示雅可比矩阵。X,y为数据；‘fun’自定义函数；beta0待定常数初值。例6：在化工生产中获得的氯气的级分y随生产时间x下降，假定在x≥8时，y与x之间有如下形式的非线性模型：)8()49.0(xbeaay现收集了44组数据，利用该数据通过拟合确定非线性模型中的待定常数。xyxyxy80.49160.43280.4180.49180.46280.40100.48180.45300.40100.47200.42300.40100.48200.42300.38100.47200.43320.41120.46200.41320.40120.46220.41340.40120.45220.40360.41120.43240.42360.36140.45240.40380.40140.43240.40380.40140.43260.41400.36160.44260.40420.39160.43260.41首先定义非线性函数的m文件：fff6.mfunctionyy=model(beta0,x)a=beta0(1);b=beta0(2);yy=a+(0.49-a)*exp(-b*(x-8));程序：x=[8.008.0010.0010.0010.0010.0012.0012.0012.0014.0014.0014.00...16.0016.0016.0018.0018.0020.0020.0020.0020.0022.0022.0024.00...24.0024.0026.0026.0026.0028.0028.0030.0030.0030.0032.0032.00...34.0036.0036.0038.0038.0040.0042.00]';y=[0.490.490.480.470.480.470.460.460.450.430.450.430.430.440.43...0.430.460.420.420.430.410.410.400.420.400.400.410.400.410.41...0.400.400.400.380.410.400.400.410.380.400.400.390.39]';beta0=[0.300.02];betafit=nlinfit(x,y,'sta67_1m',beta0)结果：betafit=0.38960.1011即：a=0.3896，b=0.1011拟合函数为：)8(1011.0)3896.049.0(3896.0xey§2插值问题在应用领域中，由有限个已知数据点，构造一个解析表达式，由此计算数据点之间的函数值，称之为插值。实例：海底探测问题某公司用声纳对海底进行测试，在5×5海里的坐标点上测得海底深度的值，希望通过这些有限的数据了解更多处的海底情况。并绘出较细致的海底曲面图。一、一、一元插值一元插值是对一元数据点(xi,yi)进行插值。1．1．线性插值：由已知数据点连成一条折线，认为相临两个数据点之间的函数值就在这两点之间的连线上。一般来说，数据点数越多，线性插值就越精确。调用格式：yi=interp1(x,y,xi,’linear’)%线性插值zi=interp1(x,y,xi,’spline’)%三次样条插值wi=interp1(x,y,xi,’cubic’)%三次多项式插值说明：yi、zi、wi为对应xi的不同类型的插值。x、y为已知数据点。例1：已知数据：x0.1.2.3.4.5.6.7.8.91y.3.511.41.61.9.6.4.81.52求当xi=0.25时的yi的值。程序：x=0:.1:1;y=[.3.511.41.61.6.4.81.52];yi0=interp1(x,y,0.025,'linear')xi=0:.02:1;yi=interp1(x,y,xi,'linear');zi=interp1(x,y,xi,'spline');wi=interp1(x,y,xi,'cubic');plot(x,y,'o',xi,yi,'r+',xi,zi,'g*',xi,wi,'k.-')legend('原始点','线性点','三次样条','三次多项式')结