利用matlab实现非线性拟合（下）
没看过上一篇的建议看一下前面的上篇。这一篇非线性拟合我就不废话，直接开始了。下面首先介绍几种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上发布
腾讯云开发者
