巧妙地进行非线性拟合——非线性拟合转化为线性拟合

1.  前言(不在意来龙去脉的可忽略不看)
对于多项式函数,可以用最小二乘法求得精确的拟合结果,使得拟合函数具有全局最优的拟合误差;对于某些非线性函数,如指数函数y=e^(ax+b),也可以对函数转化后,求得精确的拟合结果,如上述指数函数可转化为x=(ln y)/a -b/a,同样可以求得具有全局最优拟合误差的拟合函数。上述函数都可以用MATLAB的regress函数或者polyfit函数求得最优的拟合结果,或者可以用广义逆矩阵的最小二乘解计算而得。
但是,对于大多数非线性函数,难以求得拟合误差全局最优的拟合结果,形如y=c1*e^(px)+c2*e^(qx)。这一类函数,一般做法是首先预估拟合函数的参数,以此作为初始值,采用迭代优化的方法求出局部最优的拟合参数,MATLAB中用lsqcurvefit或者nonlinfit函数求出拟合函数,拟合结果的准确性由选取的初始点决定,初始点的选取极为困难。


2. 非线性拟合转化为线性拟合
对于大多数指数函数、三角函数、多项式函数通过四则运算或者复合得到的函数,通常可以用线性微分方程来表示,形如y''+ay'+by=0的常系数线性微分方程则可以表示以下形式的函数:①y=C1*e^(p*x)+C2*e^(q*x);②y=C1*(1+C2*x)*e^(w*x);③y=e^(p*x)*(C1*sin (w*x)+C2*cos (w*x))。
对于上述非线性函数,则可以通过拟合线性微分方程的系数,再计算出对应的非线性函数的参数,并利用样本数据作为初始条件,拟合出函数中的其他常数项。
微分方程需要转化为差分方程:①y'(n)=(y(n+1)-y(n-1))/(x(n+1)-x(n-1));②y''(n)=(y(n+1)-2*y(n)+y(n-1))/((x(n+1)-x(n-1))/2)^2。
通过拟合方程y''+ay'+by=0则可得到a和b,再通过求解二次方程r^2+a*r+b=0,得到两个解r1和r2。
当r1和r2不相等时,则可得到函数y=C1*e^(r1*x)+C2*e^(r2*x),再利用样本数据线性拟合则可得到C1和C2。
当r1=r2时,可得到函数y=C1*e^(r1x)+C2*x*e^(r1x),再利用样本数据线性拟合则可得到C1和C2。
当r1和r2为共轭复数时,可得函数y=e^(a*x)*(C1*sin (b*x)+C2*cos (b*x)),其中a和b分别为r1的实部和虚部,再利用样本数据线性拟合则可得到C1和C2。


事例MATLAB代码如下:
% 产生样本数据
x = linspace(0,98,50)';
y = 0.9*exp(-0.2*x)-0.3*exp(-0.06*x)+10.57;
plot(x,y)
% 计算差分
dyx = (y(3:end) - y(1:end-2))./(x(3:end) - x(1:end-2));
ddyx = (y(3:end) - 2*y(2:end-1) + y(1:end-2))./(x(3:end) - x(1:end-2)).^2*4;
% 拟合微分方程系数
A = [dyx,y(2:end-1),-1*ones(length(dyx),1)];
b = -ddyx;
a = (A.'*A)\A.'*b;
a(3) = a(3)/a(2);
% 求解二次方程
syms r
r0 = double(solve(r^2+a(1)*r+a(2)))
% 线性拟合得到常数项C1,C2
A = [exp(r0(1)*x),exp(r0(2)*x)];
C = (A.'*A)\A.'*(y-a(3))
% 作图
hold on
plot(x,C(1)*exp(r0(1)*x)+C(2)*exp(r0(2)*x)+a(3),'r.')
最后,还可以用lsqcurvefit函数进一步优化拟合结果,采用上述方法得到的参数作为初始点,可以修正差分方程拟合导致的误差。
f = @(b,x)b(1)*exp(b(2)*x)+b(3)*exp(b(4)*x)+b(5);
b0 = [C(1);r0(1);C(2);r0(2);a(3)];
b1 = lsqcurvefit(f,b0,x,y)
hold on
plot(x,f(b1,x),'k.')

阅读更多

没有更多推荐了,返回首页