第六章MATLAB解曲线拟合问题曲线拟合就是计算出两组数据之间的一种函数关系,由此可描绘其变化曲线及估计非采集数据对应的变量信息。
曲线拟合有多种方式,下面是一元函数采用最小二乘法对给定数据进行多项式曲线拟合,最后给出拟合的多项式系数。
§1 线性拟合一、数学模型y=Xβ+εβ是p⨯1的参数向量;ε是服从标准正态分布的随机干扰的n⨯1的向量;y为n⨯1的向量;X为n⨯p矩阵。
二、求解线性拟合函数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的最小二乘拟合值。
bint返回β的95%的置信区间。
r 中为形状残差,rint中返回每一个残差的95%置信区间。
Stats向量包含R2统计量、回归的F 值和p值。
三、举例例1:设y的值为给定的x的线性函数加服从标准正态分布的随机干扰值得到。
即y=10x+ε;求线性拟合方程系数。
程序:编写M文件mainG.m如下:x=[ones(10,1) (1:10) '] %ones(m,n)产生一个m行n列的元素全为1的矩阵%(1:10)产生一个10行1列的元素值从1~10的矩阵% [A B] 将矩阵A和B拼接成新矩阵或者写成[A;B]y=x*[10;1]+normrnd(0,0.1,10,1)% R = normrnd(MU,SIGMA,m,n) 依据参数MU、SIGMA生成一个随机数,m和n是R的行数和列数.[b,bint]=regress(y,x,0.05)结果:x = 1 11 21 31 41 51 61 71 81 91 10x=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.0021 0.0505 -0.5971 3.6472 -9.7295 11.3304图 6阶曲线图10阶曲线例4:再用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 = Columns 1 through 70.0000 -0.0000 0.0004 -0.0114 0.1814 -1.8065 11.2360 Columns 8 through 11-42.0861 88.5907 -92.8155 40.2671说明:可用不同阶的多项式来拟合数据,但也不是阶数越高拟合的越好。
§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 多项式曲线拟合的评价和置信区间函数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 为置信度。
例5:给出上面例1的预测值及置信度为90%的置信区间。
程序:x=0:.1:1;y=[.3 .5 1 1.4 1.6 1.9 .6 .4 .8 1.5 2]n=3;[p,s]=polyfit(x,y,n)alpha=0.05;[Y,DELTA]=polyconf(p,x,s,alpha)结果:p = 16.7832 -25.7459 10.9802 -0.0035s = R: [4x4 double]df: 7normr: 1.1406Y = Columns 1 through 7-0.0035 0.8538 1.2970 1.4266 1.3434 1.1480 0.9413Columns 8 through 110.8238 0.8963 1.2594 2.0140DELTA = Columns 1 through 71.3639 1.1563 1.1563 1.1589 1.1352 1.1202 1.1352Columns 8 through 111.1589 1.1563 1.1563 1.3639§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 ’时忽略常数项。
例6:演示一个异常数据点如何影响最小二乘拟合值与稳健拟合。
首先利用函数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)hold onplot(x,bls(1)+bls(2)*x,’:’)plot(x,brob(1)+brob(2)*x,’r‘)结果:bls = 8.4452-1.4784brob = 10.2934-2.0006分析:稳健拟合(实线)对数据的拟合程度好些,忽略了异常值。
最小二乘拟合(点线)则受到异常值的影响,向异常值偏移。
§6 向自定义函数拟合对于给定的数据,根据经验拟合为带有待定常数的自定义函数。
二、向自定义函数拟合所用函数:nlinfit( )调用格式: [beta,r,J]=nlinfit(X,y,’fun’,betao)说明:beta 返回函数’fun’中的待定常数;r 表示残差;J 表示雅可比矩阵。
X,y 为数据;‘fun’自定义函数;beta0待定常数初值。
例7:在化工生产中获得的氯气的级分y 随生产时间x 下降,假定在x ≥8时,y 与x 之间有如下形式的非线性模型:)8()49.0(---+=x b e a a y现收集了44组数据,利用该数据通过拟合确定非线性模型中的待定常数。
x y x y x y8 0.49 16 0.43 28 0.418 0.49 18 0.46 28 0.4010 0.48 18 0.45 30 0.4010 0.47 20 0.42 30 0.4010 0.48 20 0.42 30 0.3810 0.47 20 0.43 32 0.4112 0.46 20 0.41 32 0.4012 0.46 22 0.41 34 0.4012 0.45 22 0.40 36 0.4112 0.43 24 0.42 36 0.3614 0.45 24 0.40 38 0.4014 0.43 24 0.40 38 0.4014 0.43 26 0.41 40 0.3616 0.44 26 0.40 42 0.3916 0.43 26 0.41首先定义非线性函数的m 文件:fff6.mfunction yy=model(beta0,x)a=beta0(1);b=beta0(2);yy=a+(0.49-a)*exp(-b*(x-8));程序:x=[8.00 8.00 10.00 10.00 10.00 10.00 12.00 12.00 12.00 14.00 14.00 14.00...16.00 16.00 16.00 18.00 18.00 20.00 20.00 20.00 20.00 22.00 22.00 24.00...24.00 24.00 26.00 26.00 26.00 28.00 28.00 30.00 30.00 30.00 32.00 32.00...34.00 36.00 36.00 38.00 38.00 40.00 42.00]';y=[0.49 0.49 0.48 0.47 0.48 0.47 0.46 0.46 0.45 0.43 0.45 0.43 0.43 0.44 0.43...0.43 0.46 0.42 0.42 0.43 0.41 0.41 0.40 0.42 0.40 0.40 0.41 0.40 0.41 0.41...0.40 0.40 0.40 0.38 0.41 0.40 0.40 0.41 0.38 0.40 0.40 0.39 0.39]';beta0=[0.30 0.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---+=x e y。