最小二乘法曲线拟合_原理及matlab实现
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
曲线拟合(curve-fitting ):工程实践中,用测量到的一些离散的数据},...2,1,0),,{(m i y x i i =求一个近似的函数)(x ϕ来拟合这组数据,要求所得的拟合曲线能最好的反映数据的基本趋势(即使)(x ϕ最好地逼近()x f ,而不必满足插值原则。因此没必要取)(i x ϕ=i y ,只要使i i i y x -=)(ϕδ尽可能地小)。
原理:
给定数据点},...2,1,0),,{(m i y x i i =。求近似曲线)(x ϕ。并且使得近似曲线与()x f 的偏差最小。近似曲线在该点处的偏差i i i y x -=)(ϕδ,i=1,2,...,m 。
常见的曲线拟合方法:
1.使偏差绝对值之和最小
2.使偏差绝对值最大的最小
3.使偏差平方和最小
最小二乘法:
按偏差平方和最小的原则选取拟合曲线,并且采取二项式方程为拟合曲线的方法,称为最小二乘法。
推导过程:
1. 设拟合多项式为:
k k x a x a a x +++=...)(10ϕ
2. 各点到这条曲线的距离之和,即偏差平方和如下:
3. 问题转化为求待定系数0a ...k a 对等式右边求i a 偏导数,因而我们得到了:
.......
4、把这些等式化简并表示成矩阵的形式,就可以得到下面的矩阵:
5. 将这个德蒙得矩阵化简后可得到:
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包括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函数数据输入的误差是独立正态的,并且方差为常数。则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.5614x^2+0.08287x+1.15560
0.51 1.52 2.53
图1 最小二乘法曲线拟合示例
对比检验拟合的有效性:
例:在[0,π]区间上对正弦函数进行拟合,然后在[0,2π]区间画出图形,比较拟合区间和非拟合区间的图形,考察拟合的有效性。
在MATLAB中输入如下代码:
clear
x=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 =
0.0000 0.0000 -0.0003 0.0002 0.0080 0.0002 -0.1668 0.0000 1.0000 0.0000
mu =
R: [10x10 double]
df: 22
normr: 1.6178e-07
MATLAB的最优化工具箱还提供了lsqcurvefit()函数命令进行最小二乘曲线拟合(Solve nonlinear curve-fitting (data-fitting) problems in least-squares sense)。
调用格式:
x = lsqcurvefit(fun,x0,xdata,ydata)
x = lsqcurvefit(fun,x0,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] = lsqcurvefit(...)
[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,xdata)-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
0.9000 0.8435
1.0000 0.6856
1.1000 0.6100
1.2000 0.5392
1.3000 0.3946
1.4000 0.3903