《最小二乘法曲线拟合原理及matlab实现》由会员分享，可在线阅读，更多相关《最小二乘法曲线拟合原理及matlab实现（10页珍藏版）》请在金锄头文库上搜索。
1、曲线拟合curve-fitting：工程实践中，用测量到的一些离散的数据求一个近似的函数来拟合这组数据，要求所得的拟合曲线能最好的反映数据的根本趋势即使最好地逼近，而不必满足插值原如此。因此没必要取=，只要使尽可能地小。原理：给定数据点。求近似曲线。并且使得近似曲线与的偏差最小。近似曲线在该点处的偏差，i=1,2,.,m。 常见的曲线拟合方法: 1.使偏差绝对值之和最小 2.使偏差绝对值最大的最小 3.使偏差平方和最小 最小二乘法： 按偏差平方和最小的原如此选取拟合曲线，并且采取二项式方程为拟合曲线的方法,称为最小二乘法。推导过程： 1. 设拟合多项式为： 2. 各点到这条曲线的距离之和，即偏
2、差平方和如下： 3. 问题转化为求待定系数.对等式右边求偏导数，因而我们得到了： . 4、 把等式化简并表示成矩阵的形式，就可以得到下面的矩阵： 5. 将这个X德蒙得矩阵化简后可得到: 6. 也就是说X*A=Y，那么A = (X*X)-1*X*Y，便得到了系数矩阵A，同时，我们也就得到了拟合曲线。MATLAB实现：MATLAB提供了polyfit函数命令进展最小二乘曲线拟合。调用格式：p=polyfit(x,y,n) p,s= polyfit(x,y,n) p,s,mu=polyfit(x,y,n)x,y为数据点，n为多项式阶数，返回p为幂次从高到低的多项式系数向量p。x必须是单调的。矩阵s包
3、括R对x进展QR分解的三角元素、df(自由度、normr(残差用于生成值的误差估计。 p,s,mu=polyfit(x,y,n)在拟合过程中，首先对x进展数据标准化处理，以在拟合中消除量纲等影响，mu包含标准化处理过程中使用的x的均值和标准差。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函数数据输入的误差是独立正态
4、的，并且方差为常数。如此Y DELTA将至少包含50%的值。如下给定数据的拟合曲线：x=0.5,1.0,1.5,2.0,2.5,3.0，y=1.75,2.45,3.81,4.80,7.00,8.60。解：MATLAB程序如下：x=0.5,1.0,1.5,2.0,2.5,3.0;y=1.75,2.45,3.81,4.80,7.00,8.60;p=polyfit(x,y,2)x1=0.5:0.05:3.0;y1=polyval(p,x1);plot(x,y,*r,x1,y1,-b)运行结果如图1计算结果为：p =0.5614 0.8287 1.1560即所得多项式为y=0.5614x2+0.082
5、87x+1.15560 图1 最小二乘法曲线拟合示例比照检验拟合的有效性：例：在0,区间上对正弦函数进展拟合，然后在0,2区间画出图形，比拟拟合区间和非拟合区间的图形，考察拟合的有效性。在MATLAB中输入如下代码：clearx=0:0.1:pi;y=sin(x);p,mu=polyfit(x,y,9)x1=0:0.1:2*pi;y1=sin(x1);%实际曲线y2=polyval(p,x1);%根据由区间0到pi上进展拟合得到的多项式计算0到2pi上的函数值， %需要注意的是polyval返回的函数值在pi到2pi上并没有进展拟合plot(x1,y2,k*,x1,y1,k-)运行结果：p =
6、 0.0000 0.0000 -0.0003 0.0002 0.0080 0.0002 -0.1668 0.0000 1.0000 0.0000mu = R: 10x10 double df: 22 normr: 1.6178e-07MATLAB的最优化工具箱还提供了lsqcurvefit函数命令进展最小二乘曲线拟合(Solve nonlinear curve-fitting (data-fitting) problems in least-squares sense)。调用格式：x = lsqcurvefit(fun,x0,xdata,ydata)x = lsqcurvefit(fun,x0
7、,xdata,ydata,lb,ub)x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)x = lsqcurvefit(problem)x,resnorm = lsqcurvefit(.)x,resnorm,residual = lsqcurvefit(.)x,resnorm,residual,exitflag = lsqcurvefit(.)x,resnorm,residual,exitflag,output = lsqcurvefit(.)x,resnorm,residual,exitflag,output,lambda = lsqcurve
8、fit(.)x,resnorm,residual,exitflag,output,lambda,jacobian = x0为初始解向量;xdata，ydata为满足关系ydata=F(x, xdata)的数据;lb、ub为解向量的下界和上界 ，假如没有指定界，如此lb= ，ub= ;options为指定的优化参数;fun为拟合函数，其定义方式为:x = lsqcurvefit(myfun,x0,xdata,ydata)，其中myfun已定义为 function F = myfun(x,xdata)F = % 计算x处拟合函数值fun的用法与前面一样;resnorm=sum (fun(x,xda
9、ta)-ydata).2)，即在x处残差的平方和;residual=fun(x,xdata)-ydata，即在x处的残差;exitflag为终止迭代的条件;output为输出的优化信息;lambda为解x处的Lagrange乘子;jacobian为解x处拟合函数fun的jacobian矩阵。例：lsqcurvefit()优化程序Data = . 0.0000 5.8955 0.1000 3.5639 0.2000 2.5173 0.3000 1.9790 0.4000 1.8990 0.5000 1.3938 0.6000 1.1359 0.7000 1.0096 0.8000 1.0343
10、0.9000 0.8435 1.0000 0.6856 1.1000 0.6100 1.2000 0.5392 1.3000 0.3946 1.4000 0.3903 1.5000 0.5474 1.6000 0.3459 1.7000 0.1370 1.8000 0.2211 1.9000 0.1704 2.0000 0.2636;t = Data(:,1);y = Data(:,2);% axis(0 2 -0.5 6) plot(t,y,ro)title(Data points)%We would like to fit the function y = c(1)*exp(-lam(1)
11、*t) + c(2)*exp(-lam(2)*t) to the data%The lsqcurvefit function solves this type of problem easily.%To begin, define the parameters in terms of one variable x:%x(1) = c(1)%x(2) = lam(1)%x(3) = c(2)%x(4) = lam(2)%Then define the curve as a function of the parameters x and the data t:F = (x,xdata)x(1)*
12、exp(-x(2)*xdata) + x(3)*exp(-x(4)*xdata);x0 = 1 1 1 0;x,resnorm,exitflag,output = lsqcurvefit(F,x0,t,y)hold onplot(t,F(x,t)hold offFsumsquares = (x)sum(F(x,t) - y).2);opts = optimset(LargeScale,off);xunc,ressquared,eflag,outputu = .fminunc(Fsumsquares,x0,opts)fprintf(There were %d iterations using f
13、minunc, . and %d using lsqcurvefit.n, . outputu.iterations,output.iterations)fprintf(There were %d function evaluations using fminunc, . and %d using lsqcurvefit., . outputu.funcCount,output.funcCount)type fitvectorx02 = 1 0;F2 = (x,t) fitvector(x,t,y);x2,resnorm2,exitflag2,output2 = lsqcurvefit(F2,x02,t,y)fprintf(There were %d function evaluations using the 2-d . formulation, and %d using the 4-d formulation., . output2.funcCount,output.funcCount)x0bad = 5 1 1 0;xbad,resnormbad,exitflagbad,outpu
