拟合介绍:
所谓数据拟合是求一个简单的函数,例如是一个低次多项式,不要求通过已知的这些点,而是要求在整体上“尽量好”的逼近原函数。这时,在每个已知点上就会有误差,数据拟合就是从整体上使误差,尽量的小一些。
多项式拟合
n次多项式:
g(x)=c1xn+c2xn−1+⋯+cn+1
g
(
x
)
=
c
1
x
n
+
c
2
x
n
−
1
+
⋯
+
c
n
+
1
曲线与数据点的残差为:
ri=yi−g(xi),i=1,2,⋯,L
r
i
=
y
i
−
g
(
x
i
)
,
i
=
1
,
2
,
⋯
,
L
残差的平方和为:
R=∑i=1Lr2i
R
=
∑
i
=
1
L
r
i
2
为使其最小化,可令R关于
cj
c
j
的偏导数为零,即:
∂R∂cj=0,j=1,2,⋯,n+1
∂
R
∂
c
j
=
0
,
j
=
1
,
2
,
⋯
,
n
+
1
或:
∑j=1n+1(∑i=1Lx2n+2−j−ki)cj=∑i=1Lxn+1−kiyi,k=1,2,⋯,n+1
∑
j
=
1
n
+
1
(
∑
i
=
1
L
x
i
2
n
+
2
−
j
−
k
)
c
j
=
∑
i
=
1
L
x
i
n
+
1
−
k
y
i
,
k
=
1
,
2
,
⋯
,
n
+
1
或矩阵形式:
[∑i=1Lx2ni∑i=1Lx2n−1i.∑i=1Lxni ∑i=1Lx2n−1i..∑i=1Lxn−1i .... ∑i=1Lxni..∑i=1Lx0i]⎡⎣⎢⎢⎢c1c2.cn+1⎤⎦⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢∑i=1Lxniyi∑i=1Lxn−1iyi.∑i=1Lyi⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
[
∑
i
=
1
L
x
i
2
n
∑
i
=
1
L
x
i
2
n
−
1
.
∑
i
=
1
L
x
i
n
∑
i
=
1
L
x
i
2
n
−
1
.
.
∑
i
=
1
L
x
i
n
−
1
.
.
.
.
∑
i
=
1
L
x
i
n
.
.
∑
i
=
1
L
x
i
0
]
[
c
1
c
2
.
c
n
+
1
]
=
[
∑
i
=
1
L
x
i
n
y
i
∑
i
=
1
L
x
i
n
−
1
y
i
.
∑
i
=
1
L
y
i
]
多项式拟合MATLAB命令:polyfit
**格式:**p=polyfit(x,y,n)
例1:已知的数据点来自
f(x)=(x2−3x=5)e−5xsinx
f
(
x
)
=
(
x
2
−
3
x
=
5
)
e
−
5
x
sin
x
,用多项式拟合的方法在不同的阶次下进行拟合。
拟合该数据的3次多项式:
>> x0=0:.1:1; y0=(x0.^2-3*x0+5).*exp(-5*x0).*sin(x0);
>> p3=polyfit(x0,y0,3); vpa(poly2sym(p3),10) % 可以如下显示多项式
ans =
2.839962923*x^3-4.789842696*x^2+1.943211631*x+.5975248921e-1
绘制拟合曲线:
>> x=0:.01:1;
>>ya=(x.^2-3*x+5).*exp(-5*x).*sin(x);
>> y1=polyval(p3,x);
>>plot(x,y1,x,ya,x0,y0,'o')
就不同的次数进行拟合:
p4=polyfit(x0,y0,4); y2=polyval(p4,x);
p5=polyfit(x0,y0,5); y3=polyval(p5,x);
p8=polyfit(x0,y0,8); y4=polyval(p8,x);
plot(x,y0,x0,y0,'o',x,y2,'r',x,y3,'g',x,y4,'k')
拟合最高次数为8的多项式:
>> vpa(poly2sym(p8),5)
ans =
- 8.2586*x^8 + 43.566*x^7 - 101.98*x^6 + 140.22*x^5 - 125.29*x^4 + 74.45*x^3 - 27.672*x^2 + 4.9869*x + 4.2037e-7
Taylor幂级数展开:
不能运行
>> syms x; y=(x^2-3*x+5)*exp(-5*x)*sin(x);
>> vpa(taylor(y,9),5)
ans =
9.*x-29.*x^2+77.667*x^3-142.*x^4+192.17*x^5-204.96*x^6+179.13*x^7-131.67*x^8
多项式表示数据模型是不唯一的,即是两个多项式函数完全不同。在某一区域内其曲线可能特别近似。
例2:对
f(x)=1/(1+25x2),−1≤x≤1
f
(
x
)
=
1
/
(
1
+
25
x
2
)
,
−
1
≤
x
≤
1
进行多项式拟合,
多项式拟合的效果并不一定总是很精确的。
>> x0=-1+2*[0:10]/10; y0=1./(1+25*x0.^2);
>> x=-1:.01:1; ya=1./(1+25*x.^2);
>> p3=polyfit(x0,y0,3);
>> y1=polyval(p3,x);
>> p5=polyfit(x0,y0,5);
>> y2=polyval(p5,x);
>> p8=polyfit(x0,y0,8);
>> y3=polyval(p8,x);
>> p10=polyfit(x0,y0,10);
>> y4=polyval(p10,x);
>> plot(x,ya,x,y1,x,y2,'-.',x,y3,'--',x,y4,':')
用Taylor幂级数展开效果将更差
>> syms x; y=1/(1+25*x^2); p=taylor(y,x,10)
p =
1-25*x^2+625*x^4-15625*x^6+390625*x^8
多项式拟合效果
>> x1=-1:0.01:1;
>> ya=1./(1+25*x1.^2);
>> y1=subs(p,x,x1);
>> plot(x1,ya,'--‘,x1,y1)
函数线性组合的曲线拟合方法
已知某函数的线性组合为:
g(x)=c1f1(x)+c2f2(x)+c3f3(x)+....+cnfn(x)
g
(
x
)
=
c
1
f
1
(
x
)
+
c
2
f
2
(
x
)
+
c
3
f
3
(
x
)
+
.
.
.
.
+
c
n
f
n
(
x
)
其中
f1(x)f2(x)f3(x).....fn(x)
f
1
(
x
)
f
2
(
x
)
f
3
(
x
)
.
.
.
.
.
f
n
(
x
)
为已知函数
c1,c2c3.....cn
c
1
,
c
2
c
3
.
.
.
.
.
c
n
为待定系数
假设已经测出数据
(x1,y1),(x2,y2),.....(xM,yM)
(
x
1
,
y
1
)
,
(
x
2
,
y
2
)
,
.
.
.
.
.
(
x
M
,
y
M
)
则可建立如下的线性方程
Ac=y
A
c
=
y
例3:
直接拟合ci参数:
>> x=[0,0.2,0.4,0.7,0.9,0.92,0.99,1.2,1.4,1.48,1.5]';
>> y=[2.88;2.2576;1.9683;1.9258;2.0862;2.109; 2.1979;2.5409;2.9627;3.155;3.2052];
>> A=[ones(size(x)),exp(-3*x), cos(-2*x).*exp(-4*x) ,x.^2];
>> c=A\y; c1=c'
c1 =
1.2200 2.3397 -0.6797 0.8700
例4:
>> x=[1.1052,1.2214,1.3499,1.4918,1.6487,1.8221,2.0138,...
2.2255,2.4596,2.7183,3.6693];
>> y=[0.6795,0.6006,0.5309,0.4693,0.4148,0.3666,0.3241,...
0.2864,0.2532,0.2238,0.1546];
>> plot(x,y,x,y,'*')
分别对x,y进行对数变换:
>> x1=log(x); y1=log(y); plot(x1,y1)
>> A=[x1', ones(size(x1'))]; c=[A\y1']
c =
-1.2339 -0.2630
>> exp(c(2))
ans =
0.7687
例5:
>> x=[0:0.1:1]'; y=(x.^2-3*x+5).*exp(-5*x).*sin(x); n=8; A=[];
>> for i=1:n+1, A(:,i)=x.^(n+1-i); end
>> c=A\y; vpa(poly2sym(c),5)
ans =
-9.2586*x^8+43.566*x^7-101.98*x^6+140.22*x^5-129.29*x^4+74.450*x^3-27.672*x^2+4.9869*x+.42037e-6
最小二乘曲线拟合
格式: [a, jm]=lsqcurvefit(Fun,a0,x,y)
例6:
由下面语句生成一组数据,其中ai为待定系数,
>> x=0:.1:10;
>> y=0.12*exp(-0.213*x)+0.54*exp(-0.17*x).*sin(1.23*x);
>> f=inline('a(1)*exp(-a(2)*x)+a(3)*…
exp(-a(4)*x).*sin(a(5)*x)','a','x');
得出待定系数向量:
>> [xx,res]=lsqcurvefit(f,[1,1,1,1,1],x,y); xx',res
Optimization terminated successfully:
Relative function value changing by less than OPTIONS.TolFun
ans =
0.1200
0.2130
0.5400
0.1700
1.2300
res =
1.7927e-16
修改最优化选项:
不能运行
>> ff=optimset; ff.TolFun=1e-20; ff.TolX=1e-15; % 修改精度限制
>> [xx,res]=lsqcurvefit(f,[1,1,1,1,1],x,y,[],[],ff); xx‘,res % []变量界
Optimization terminated successfully:
Relative function value changing by less than OPTIONS.TolFun
ans =
0.1200
0.2130
0.5400
0.1700
1.2300
res =
9.5035e-021
例7:
>> x=0.1:0.1:1;
>> y=[2.3201,2.6470,2.9707,3.2885,3.6008,3.9090,4.2147,4.5191,
4.8232,9.1275];
function y=c8f3(a,x)
y=a(1)*x+a(2)*x.^2.*exp(-a(3)*x)+a(4);
求解
a=lsqcurvefit('c8f3',[1;2;2;3],x,y); a'
Maximum number of function evaluations exceeded;
increase options.MaxFunEvals
ans =
2.4575 2.4557 1.4437 2.0720
绘制曲线:
>> y1=c8f3(a,x); plot(x,y,x,y1,’o’)
跑出图片