利用matlab实现非线性拟合(三维、高维、参数方程) |
您所在的位置:网站首页 › origin曲线拟合求公式 › 利用matlab实现非线性拟合(三维、高维、参数方程) |
利用matlab实现非线性拟合[三维、高维、参数方程]
0 前言1 线性拟合1.1 多项式拟合1.2 线性拟合
2 一维非线性拟合2.1 简单的非线性拟合2.2 matlab中Curve Fitting App2.3 matlab中非线性拟合的实现2.3.1 fit()函数2.3.2 nlinfit()函数2.3.3 lsqnonlin()函数和lsqcurvefit()函数2.3.4 fsolve()函数2.3.5 粒子群算法2.3.6 不同算法的对比效果
3 高维非线性方程组拟合3.1 一般形式高维方程或方程组的拟合3.2 一般形式参数方程的拟合3.3 补充
参考文献
2021.06 更新。补充了非线性拟合中,关于最小二乘定义的问题,放在了最后一章。 2021.12更新。应评论区建议,补充了3.3章绘图用的代码。 本文首发于 matlab爱好者 微信公众号,欢迎关注。 惯例声明:本人没有相关的工程应用经验,只是纯粹对相关算法感兴趣才写此博客。所以如果有错误,欢迎在评论区指正,不胜感激。本文主要关注于算法的实现,对于实际应用等问题本人没有任何经验,所以也不再涉及。 本人在学习matlab匿名函数时,作为练习有感而发,写下下面内容,非专业,难免有误。 一般而言,通过已有的数据点去推导其它数据点,常见的方法有插值和拟合。插值适用性较广,尤其是线性插值或样条插值已被广泛的应用。但是通过已知的函数去拟合数据,是连接理论与实验重要的桥梁,这一点是插值无法替代的。 日常学习工作中,经常会遇到下面这种问题:想要用某个具体函数去拟合自己的数据,明明知道这个函数的具体形式,却不知道其中的参数怎么选取。本文就简单介绍一下matlab环境下,如何进行非线性拟合。 1 线性拟合 1.1 多项式拟合多项式拟合就是利用下面形式的方程去拟合数据: y = a + b x + c x 2 + d x 3 + . . . y=a +bx+cx^2+dx^3+... y=a+bx+cx2+dx3+... matlab中可以用polyfit()函数进行多项式拟合。下面举一个小例子: 对于已有的数据点,我们采用4阶多项式拟合。其中已知函数的的表达式为 y=0.03 x^4 - 0.5 x^3 + 2 x^2 - 4在此基础上添加了一些噪声点。拟合曲线依然采用4阶进行拟合,结果如下。
对于多项式拟合,不是阶数越多越好。比如还是这个上面这个例子,阶数越多,曲线越能够穿过每一个点,但是曲线的形状与理论曲线偏离越大。所以选择最适合的才是最好的。 线性拟合就是能够把拟合函数写成下面这种形式的: y = p 1 f 1 ( x ) + p 2 f 2 ( x ) + p 3 f 3 ( x ) + . . . y=p_1 f_1(x) +p_2 f_2(x)+p_3 f_{3}(x)+... y=p1f1(x)+p2f2(x)+p3f3(x)+... 其中f(x)是关于x的函数,其表达式是已知的。p是常数系数,这个是未知的。 对于这种形式的拟合,matlab内部有一个及其强悍的函数,可以自动输出p的解,并且满足最小二乘。这个函数就是\。没错,就是斜杠。这个符号通常用于求解方程AX=B的情况,我们用X=A\B可以求出未知数X。我们利用当A行和列不等时,输出X的最小二乘这个特性,就可以求出相应的最佳拟合。 还是举个例子
事实上,其实常用的拟合方式中,有很多都是线性拟合,比如多项式拟合,傅里叶拟合等。对于多项式拟合,只需要把拟合的部分替换成x,x.^2,x.^3这种形式。对于傅里叶级数,只需要把拟合的部分替换为sin(x),sin(2*x),sin(3*x)这种形式就行。 下面展示一个利用线性拟合,进行不同频率的三角波级数拟合正弦函数的例子: 前面介绍了线性的拟合方法。如果一个非线性方程,可以化为上面线性方程中公式中给出的样子,那么我们也可以套用线性的方法去求解。 比如下面的方程:
y
=
a
∗
exp
(
−
b
x
)
y=a*\exp(-bx)
y=a∗exp(−bx) 经过取对数变换,可以变为下面等效的线性形式:
lg
(
y
)
=
lg
(
a
)
+
b
∗
(
−
x
)
\lg(y)=\lg(a)+b*(-x)
lg(y)=lg(a)+b∗(−x) 这样求出来之后,再反变换回去,就可以得到原方程的系数了。 matlab自带了一个Curve Fitting App,它在matlab集成的App里面。 matlab中,fit()函数是一个比较通用的用于函数拟合的函数。它的优点就是非常的全面,可以用于各种种类的拟合。上面的App里,很多拟合种类都是间接调用了fit函数来实现的拟合。 对于非线性拟合,可以使用fit()函数中的Nonlinear Least Squares方法。其大概原理为,首先确定一个初始的点,计算该点的雅可比矩阵J,来近似线性化的判断该点周围的趋势,并将这个点向更小的方向移动。 因此,这个方法的一个缺点在于,对于初始点的选取非常敏感,最终结果只能在初始点附近的局部最小值点上,而不能保证全局最小值。 2.3.2 nlinfit()函数相比于前面的fit()函数,nlinfit()函数是matlab专门的非线性拟合函数。对于非稳健估计,采用的是Levenberg-Marquardt(LM)方法,也叫阻尼最小二乘法。对于稳健估计,采用的是Iteratively Reweighted Least Squares方法,也就是在Least Squares基础上,对每一个拟合点的权重进行调整的一种方法。这两者方法也都是基于雅克比矩阵的方法。 2.3.3 lsqnonlin()函数和lsqcurvefit()函数lsqnonlin()也是matlab中自带的一个非线性拟合函数。它给出了两种计算非线性拟合的方法,一种是比较经典之前介绍过的LM方法,一种是信赖域方法。信赖域法(trust region reflective)是通过Hessian 矩阵,逐步试探邻域内的最小化,来求解问题的。相比较之前的那些雅克比相关的方法,信赖域法会占用更多内存和速度,所以适用于中小规模的矩阵。 lsqcurvefit()函数和lsqnonlin()内容上相似,只是引用格式上有所不同。 2.3.4 fsolve()函数这也是一个求解非线性方程的函数,可以求解方程组或者矩阵形式,功能非常强大。默认的算法为trust-region-dogleg,俗称狗腿法,属于信赖域法。这里用到的功能比较基础,所以也不过多介绍了。 2.3.5 粒子群算法说了那么多,发现逐渐从如何非线性拟合,陷入到了最优化的深坑里。而且前面的那么多方法,很多都解决不了陷入到局部最优解的问题。实际上,这种问题如果进入了最优化领域,很多智能算法也可以被考虑进来。所以我也把粒子群PSO算法加入到了里面,尝试将结果收敛到全局最优解。 2.3.6 不同算法的对比效果第一个例子是 y=a.*exp(-b\*(x-c).^2)+d,一个简单的高斯函数形式的非线性方程,其参数给定为: abcd3.82.14.4-1.3在已知函数形式,求解这四个参数条件下,6种不同的函数非拟合效果如下: 第二个例子是一个指数增长的正弦函数,在很多线性系统中都可以测量到这种信号。函数的形式为: y=a*x+b*sin(c*x).*exp(d*x)+e 。其给定的参数为: abcde-0.32.14.40.31.7这个函数的拟合具有一定难度,拟合过程中会遇到非常多的局部解。
代码如下: 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算法')下图展示了PSO求解过程中逐渐收敛到全局最优解的过程。 之前的文章中的数据具有一 一对应的特点,所以严格来讲并不是普遍的二维拟合。对于一些复杂的二维函数,比如椭圆,可能原本的拟合就会失效。 对于这种一般性质的方程或方程组,将原本输入方程整理为f(x,y,z,…)=0的形式。衡量拟合程度的优化函数,就直接取函数f(xi,yi,zi,…)的值即可。 下面演示最终的两个例子: 第一个是三维直线,采用两平面式描述。 Ax+By+Cz-1=0 Dx+Ey+Fz-1=0总共2个方程,维度为3维,第一个方程有3个参数,第二个方程也有3个参数。离散点已知的条件下,三维直线的两平面表达式不唯一。 最终拟合效果如下: 第二个是二维椭圆,方程为: x^2+Axy+By^2+Cx+Dy+E=0总共1个方程,维度为2维。方程共有5个参数。 最终拟合效果如下: 高维参数方程的拟合比较困难。因为原本方程中x、y、z的坐标点都是已知的。但是参数方程中,x、y、z的坐标点已知,但是与参数u、v往往未知。所以相当于原本的方程中引入了额外的未知信息。 但是基本思路和普通方程是一样的。离散点距离假想方程的距离,需要用该点距方程上每个点的距离中的最小值,来近似判断。 依然是展示两个例子。 第一个是计算三维李萨如图形。参数方程为: x=sin(A*u) y=cos(B*u) z=sin(C*u)方程为三维参数方程,只有一个参数u。第一个方程有1个未知量A,第二个方程有1个未知量B,第三方程有1个未知量C。 最终拟合效果如下。由于李萨如图形中,只要频率的比例固定,图案就会固定。所以最终ABC的值不唯一,但是它们的比例肯定唯一。
方程为三维参数方程,有2个参数u、v。第一个方程有2个未知量AB,第二个方程有2个未知量CD,第三方程有0个未知量。 最终拟合效果如下:
主要函数就是Dis和Dis_P,准备工作就是把方程和离散点整理成可以输入的形式。优化用到的函数就是PSO(particleswarm),需要更改未知参数数量和范围就可以。 3.3 补充这一部分的拟合优度或者置信区间就没办法拿出来说了,毕竟有试凑的嫌疑。适合单纯的验证。 另外这里和最小二乘求出来的结果会有所不同。这主要是定义的问题。其中最小二乘指的是距离y方向上最小的二乘距离。而这一章节中定义的最小二乘,是指的点到拟合曲线距离的最小二乘(有点类似于主成分分析)。虽然一般情况下两者并不会有太大差别,但是如果有误差的话,还请注意这一点, 如下图,这里是一个椭圆分布的离散点。 单独在后面补充了这么一个章节,希望在应用中注意。 上面的图,绘图代码如下: %构建原始数据 t1=rand(3000,1);t2=rand(3000,1); x=2*sqrt(t2).*cos(2*pi*t1); y=sqrt(t2).*sin(2*pi*t1); XY=[x,y]; XY2=XY*[1,1;-1,1];%变形 x2=XY2(:,1);y2=XY2(:,2); x=x2;y=y2; %% 演示1 %线性最小二乘拟合,以|y-yi|^2作为衡量指标 yn1=x; yn2=ones(size(x)); X=[yn1,yn2]; %拟合,得到系数 Pn=X\y; yn=Pn(1)*yn1+Pn(2)*yn2; %% 演示2 %二维直线非线性拟合,以|ax+by+c|^2作为衡量指标 %准备数据 XX={x,y}; %准备函数 F1=@(p1,XX1) p1(1)*XX1{1}-p1(2)*XX1{2}+p1(3);%ax-by+c=0 FF={F1}; %3 生成最终优化函数,带入到优化方程中求解 fun=@(p) Dis(p,{3},XX,FF);%p参数,{2}第一个方程2个参数。XX离散点。FF函数表达式。 % Sum_N=Dis([1,0],{2},XX,FF) OP=optimoptions('particleswarm','FunctionTolerance',0); [p,fval,~,~]=particleswarm(fun,3,[1,1,-1],[5,5,1],OP); %% 对比 figure() hold on plot(x,y,'.') plot(x,yn,'linewidth',4) plot(x,(p(1)*x+p(3))/p(2),'linewidth',4) hold off legend({'原始数据','最小二乘','非线性拟合'}) 参考文献信赖域法 https://zhuanlan.zhihu.com/p/205555114 |
今日新闻 |
推荐新闻 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |