没看过上一篇的建议看一下前面的上篇。这一篇非线性拟合我就不废话,直接开始了。下面首先介绍几种matlab非线性拟合方法,之后将这几种方法进行对比研究。
如果你喜欢界面化的输入输出,那么可以尝试Curve Fitting App,它在matlab集成的App里面。
界面里常用的拟合方式都有,而且直接展示拟合效果,非常方便。非常适合鼠标直接拖拖拽拽点点点的操作方式。
除了界面拟合,下面介绍几种函数式拟合的方式。
1 fit()函数
matlab中,fit()函数是一个比较通用的用于函数拟合的函数。它的优点就是非常的全面,可以用于各种种类的拟合。上面的App里,很多拟合种类都是间接调用了fit函数来实现的拟合。
对于非线性拟合,可以使用fit()函数中的Nonlinear Least Squares方法。其大概原理为,首先确定一个初始的点,计算该点的雅可比矩阵J,来近似线性化的判断该点周围的趋势,并将这个点向更小的方向移动。
因此,这个方法的一个缺点在于,对于初始点的选取非常敏感,最终结果只能在初始点附近的局部最小值点上,而不能保证全局最小值。
2 nlinfit()函数
相比于前面的fit()函数,nlinfit()函数是matlab专门的非线性拟合函数。对于非稳健估计,采用的是Levenberg-Marquardt(LM)方法,也叫阻尼最小二乘法。对于稳健估计,采用的是Iteratively Reweighted Least Squares方法,也就是在Least Squares基础上,对每一个拟合点的权重进行调整的一种方法。这两者方法也都是基于雅克比矩阵的方法。
3 lsqnonlin()函数和lsqcurvefit()函数
lsqnonlin()也是matlab中自带的一个非线性拟合函数。它给出了两种计算非线性拟合的方法,一种是比较经典之前介绍过的LM方法,一种是信赖域方法。信赖域法(trust region reflective)是通过Hessian 矩阵,逐步试探邻域内的最小化,来求解问题的。相比较之前的那些雅克比相关的方法,信赖域法会占用更多内存和速度,所以适用于中小规模的矩阵。
lsqcurvefit()函数和lsqnonlin()内容上相似,只是引用格式上有所不同。
4 fsolve()函数
这也是一个求解非线性方程的函数,可以求解方程组或者矩阵形式,功能非常强大。默认的算法为trust-region-dogleg,俗称狗腿法,属于信赖域法。这里用到的功能比较基础,所以也不过多介绍了。
5 粒子群算法
说了那么多,发现逐渐从如何非线性拟合,陷入到了最优化的深坑里。而且前面的那么多方法,很多都解决不了陷入到局部最优解的问题。实际上,这种问题如果进入了最优化领域,很多智能算法也可以被考虑进来。所以我也把粒子群PSO算法加入到了里面,尝试将结果收敛到全局最优解。
前面介绍的这些方法究竟效果如何,下面用实际例子比试一下。
第一个例子是 y=a.*exp(-b*(x-c).^2)+d,一个简单的高斯函数形式的非线性方程,其参数给定为:
a | b | c | d |
3.8 | 2.1 | 4.4 | -1.3 |
在已知函数形式,求解这四个参数条件下,6种不同的函数非拟合效果如下:
可以看到,这几种方法都能够较好的拟合出想要的结果。
第二个例子是一个指数增长的正弦函数,在很多线性系统中都可以测量到这种信号。函数的形式为:
y=a*x+b*sin(c*x).*exp(d*x)+e 。其给定的参数为:
a | b | c | d | e |
-0.3 | 2.1 | 4.4 | 0.3 | 1.7 |
这个函数的拟合具有一定难度,拟合过程中会遇到非常多的局部解。
其中前面的几种方法对于初始值的敏感度比较高,如果初始值选的比较接近原始解,也是可以得到较好的结果。其中nlinfit函数经常会报错,容错率较低。而PSO算法经常能够收敛到最优解(虽然不是每次都可以,偶尔也会陷入局部解)。
上图展示了PSO算法逐渐收敛到全局最优解的过程。
6种算法对比的代码如下:
clear
clc
%函数大比拼
close all
%初始设置
x = 0:0.1:10;
a = -0.3;
b = 2.1;
c = 4.4;
d = 0.3;
e = 1.7;
y = a*x+b*sin(c*x).*exp(d*x)+e;
noise = 0.05*abs(y-1).*randn(size(x));
y = y+noise;%加噪声函数
figure();%plot(x,y)
y_lim = [-40,40];
%% 1 fit()函数 Least Squares
ft = fittype( 'a*x+b*sin(c*x).*exp(d*x)+e', 'independent', 'x', 'dependent', 'y' );
OP1 = fitoptions( 'Method', 'NonlinearLeastSquares' );
OP1.StartPoint = 5*rand(1,5);%初始值,越靠近真实值越好
OP1.Lower = [-2,0,2,0,0];%参数的最小边界
OP1.Upper = [1,3,5,3,3];%参数的最大边界
%拟合
fitobject = fit(x',y',ft,OP1);
a2=fitobject.a;
b2=fitobject.b;
c2=fitobject.c;
d2=fitobject.d;
e2=fitobject.e;
%展示结果
y1 = a2*x+b2*sin(c2*x).*exp(d2*x)+e2;
subplot(3,2,1)
hold on
plot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')
plot(x,y1,'-','linewidth',1.5,'color','r')
hold off
box on
ylim(y_lim)
title('fit函数')
%% 2 nlinfit()函数 Levenberg-Marquardt %容易报错
modelfun = @(p,x)( p(1)*x+p(2)*sin(p(3)*x).*exp(p(4)*x)+p(5) );
OP2 = statset('nlinfit');
%opts.RobustWgtFun = 'bisquare';
p0 = 5*rand(1,5);
p = nlinfit(x,y,modelfun,p0,OP2);
%拟合
y2 = p(1)*x+p(2)*sin(p(3)*x).*exp(p(4)*x)+p(5);
subplot(3,2,2)
hold on
plot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')
plot(x,y2,'-','linewidth',1.5,'color','r')
hold off
box on
ylim(y_lim)
title('nlinfit函数')
%% 3 lsqnonlin()函数 trust-region-reflective
modelfun = @(p)(p(1)*x+p(2)*sin(p(3)*x).*exp(p(4)*x)+p(5) -y) ;%这里和nlinfit()函数定义不一样
p0 = 5*rand(1,5);
OP3 = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective','MaxFunctionEvaluations',1e4,'MaxIterations',1e3);
[p,~] = lsqnonlin(modelfun,p0,[-2,0,2,0,0],[1,3,5,3,3],OP3);
y3 = p(1)*x+p(2)*sin(p(3)*x).*exp(p(4)*x)+p(5);
subplot(3,2,3)
hold on
plot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')
plot(x,y3,'-','linewidth',1.5,'color','r')
hold off
box on
ylim(y_lim)
title('lsqnonlin函数')
%% 4 lsqcurvefit()函数 trust-region-reflective
modelfun = @(p,x)(p(1)*x+p(2)*sin(p(3)*x).*exp(p(4)*x)+p(5)) ;%这里和其它函数的定义也不一样
p0 = 5*rand(1,5);
OP4 = optimoptions('lsqcurvefit','Algorithm','trust-region-reflective','MaxFunctionEvaluations',1e4,'MaxIterations',1e3);
[p,~] = lsqcurvefit(modelfun,p0,x,y,[-2,0,2,0,0],[1,3,5,3,3],OP4);
y4 = p(1)*x+p(2)*sin(p(3)*x).*exp(p(4)*x)+p(5);
subplot(3,2,4)
hold on
plot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')
plot(x,y4,'-','linewidth',1.5,'color','r')
hold off
box on
ylim(y_lim)
title('lsqcurvefit函数')
%% 5 fsolve()函数 %默认算法trust-region-dogleg
modelfun = @(p)(p(1)*x+p(2)*sin(p(3)*x).*exp(p(4)*x)+p(5) -y);
p0 = 5*rand(1,5);
OP5 = optimoptions('fsolve','Display','off');
p = fsolve(modelfun,p0,OP5);
y5 = p(1)*x+p(2)*sin(p(3)*x).*exp(p(4)*x)+p(5);
subplot(3,2,5)
hold on
plot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')
plot(x,y5,'-','linewidth',1.5,'color','r')
hold off
box on
ylim(y_lim)
title('fsolve函数')
%% 6 粒子群PSO算法
fun = @(p) ( norm(p(1)*x+p(2)*sin(p(3)*x).*exp(p(4)*x)+p(5) -y) );%这里需要计算误差的平方和
OP6 = optimoptions('particleswarm','InertiaRange',[0.4,1.2],'SwarmSize',100);
[p,~,~,~] = particleswarm(fun,5,[-5,-5,-5,-5],[5,5,5,5],OP6);%区间可以稍微放大一些,不怕
y6 = p(1)*x+p(2)*sin(p(3)*x).*exp(p(4)*x)+p(5);
subplot(3,2,6)
hold on
plot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')
plot(x,y6,'-','linewidth',1.5,'color','r')
hold off
box on
ylim(y_lim)
title('PSO算法')
参考资料:matlab官方文档
图片来源:由 hyhhyh21 在Pixabay上发布