matlab之数据拟合

拟合介绍：

所谓数据拟合是求一个简单的函数,例如是一个低次多项式,不要求通过已知的这些点,而是要求在整体上“尽量好”的逼近原函数。这时,在每个已知点上就会有误差,数据拟合就是从整体上使误差,尽量的小一些。

多项式拟合

n次多项式：g(x)=c1xn+c2xn−1+⋯+cn+1g(x)=c1xn+c2xn−1+⋯+cn+1g(x) = {c_1}{x^n} + {c_2}{x^{n - 1}} + \cdots + {c_{n + 1}}
曲线与数据点的残差为：ri=yi−g(xi),i=1,2,⋯,Lri=yi−g(xi),i=1,2,⋯,L{r_i} = {y_i} - g({x_i}),\quad i = 1,2, \cdots ,L
残差的平方和为：R=∑i=1Lr2iR=∑i=1Lri2R = \sum\limits_{i = 1}^L {r_i^2}
为使其最小化，可令R关于cjcj{c_j}的偏导数为零，即：∂R∂cj=0,j=1,2,⋯,n+1∂R∂cj=0,j=1,2,⋯,n+1\frac{{\partial R}}{{\partial {c_j}}} = 0,\quad j = 1,2, \cdots ,n + 1
或:
∑j=1n+1(∑i=1Lx2n+2−j−ki)cj=∑i=1Lxn+1−kiyi,k=1,2,⋯,n+1∑j=1n+1(∑i=1Lxi2n+2−j−k)cj=∑i=1Lxin+1−kyi,k=1,2,⋯,n+1\sum\limits_{j = 1}^{n + 1} {(\sum\limits_{i = 1}^L {x_i^{2n + 2 - j - k}} } ){c_j} = \sum\limits_{i = 1}^L {x_i^{n + 1 - k}} {y_i},\;k = 1,2, \cdots ,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=1Lxi2n∑i=1Lxi2n−1.∑i=1Lxin∑i=1Lxi2n−1..∑i=1Lxin−1....∑i=1Lxin..∑i=1Lxi0][c1c2.cn+1]=[∑i=1Lxinyi∑i=1Lxin−1yi.∑i=1Lyi]\left[ {\begin{array}{*{20}{c}} {\sum\limits_{i = 1}^L {x_i^{2n}} }^L {x_i^{2n - 1}} }&.^L {x_i^n} }\ {\sum\limits_{i = 1}^L {x_i^{2n - 1}} }&.&.^L {x_i^{n - 1}} }\ .&.&.&.\ {\sum\limits_{i = 1}^L {x_i^n} }&.&.^L {x_i^0} } \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{c_1}}\\ {{c_2}}\\ .\\ {{c_{n + 1}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\sum\limits_{i = 1}^L {x_i^n} {y_i}}\\ {\sum\limits_{i = 1}^L {x_i^{n - 1}{y_i}} }\\ .\\ {\sum\limits_{i = 1}^L {y_i^{}} } \end{array}} \right]

多项式拟合MATLAB命令：polyfit
**格式：**p=polyfit(x,y,n)

例1：已知的数据点来自f(x)=(x2−3x=5)e−5xsinxf(x)=(x2−3x=5)e−5xsin⁡xf(x) = ({x^2} - 3x = 5){e^{ - 5x}}\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≤1f(x)=1/(1+25x2),−1≤x≤1f(x) = 1/(1 + 25x2), - 1 \le x \le 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)=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)f1(x)f2(x)f3(x).....fn(x){f_1}(x){f_2}(x){f_3}(x).....{f_n}(x)为已知函数
c1,c2c3.....cnc1,c2c3.....cn{c_1},{c_2}{c_3}.....{c_n}为待定系数
假设已经测出数据(x1,y1),(x2,y2),.....(xM,yM)(x1,y1),(x2,y2),.....(xM,yM)({x_1},{y_1}),({x_2},{y_2}),.....({x_M},{y_M})
则可建立如下的线性方程
Ac=yAc=yAc = 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.TolFunans =0.12000.21300.54000.17001.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.12000.21300.54000.17001.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’)

跑出图片