1200字范文,内容丰富有趣,写作的好帮手!
1200字范文 > 利用matlab实现非线性拟合(三维 高维 参数方程)

利用matlab实现非线性拟合(三维 高维 参数方程)

时间:2024-05-25 00:07:40

相关推荐

利用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 补充参考文献

.06 更新。补充了非线性拟合中,关于最小二乘定义的问题,放在了最后一章。

.12更新。应评论区建议,补充了3.3章绘图用的代码。

本文首发于 matlab爱好者 微信公众号,欢迎关注。

惯例声明:本人没有相关的工程应用经验,只是纯粹对相关算法感兴趣才写此博客。所以如果有错误,欢迎在评论区指正,不胜感激。本文主要关注于算法的实现,对于实际应用等问题本人没有任何经验,所以也不再涉及。

本人在学习matlab匿名函数时,作为练习有感而发,写下下面内容,非专业,难免有误。

0 前言

一般而言,通过已有的数据点去推导其它数据点,常见的方法有插值和拟合。插值适用性较广,尤其是线性插值或样条插值已被广泛的应用。但是通过已知的函数去拟合数据,是连接理论与实验重要的桥梁,这一点是插值无法替代的。

日常学习工作中,经常会遇到下面这种问题:想要用某个具体函数去拟合自己的数据,明明知道这个函数的具体形式,却不知道其中的参数怎么选取。本文就简单介绍一下matlab环境下,如何进行非线性拟合。

1 线性拟合

1.1 多项式拟合

多项式拟合就是利用下面形式的方程去拟合数据:

y=a+bx+cx2+dx3+...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阶进行拟合,结果如下。

可以看到拟合曲线与理论曲线基本一致,说明这种方法能够较好的拟合出原始数据的趋势。源代码如下:

x=0:0.5:10;y=0.03*x.^4-0.5*x.^3+2.0*x.^2-0*x-4+6*(rand(size(x))-0.5);p=polyfit(x,y,4);x2=0:0.05:10;y2=polyval(p,x2);figure();subplot(1,2,1)hold onplot(x,y,'linewidth',1.5,'MarkerSize',15,'Marker','.','color','k')plot(x,0.03*x.^4-0.5*x.^3+2.0*x.^2-0*x-4,'linewidth',1,'color','b')hold offlegend('原始数据点','理论曲线','Location','southoutside','Orientation','horizontal')legend('boxoff')box onsubplot(1,2,2)hold onplot(x2,y2,'-','linewidth',1.5,'color','r')plot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')hold offbox onlegend('拟合曲线','数据点','Location','southoutside','Orientation','horizontal')legend('boxoff')

对于多项式拟合,不是阶数越多越好。比如还是这个上面这个例子,阶数越多,曲线越能够穿过每一个点,但是曲线的形状与理论曲线偏离越大。所以选择最适合的才是最好的。

1.2 线性拟合

线性拟合就是能够把拟合函数写成下面这种形式的:

y=p1f1(x)+p2f2(x)+p3f3(x)+...y=p_1 f_1(x) +p_2 f_2(x)+p_3 f_{3}(x)+... y=p1​f1​(x)+p2​f2​(x)+p3​f3​(x)+...

其中f(x)是关于x的函数,其表达式是已知的。p是常数系数,这个是未知的。

对于这种形式的拟合,matlab内部有一个及其强悍的函数,可以自动输出p的解,并且满足最小二乘。这个函数就是\。没错,就是斜杠。这个符号通常用于求解方程AX=B的情况,我们用X=A\B可以求出未知数X。我们利用当A行和列不等时,输出X的最小二乘这个特性,就可以求出相应的最佳拟合。

还是举个例子

虽然看上去很非线性,但是x和y都是已知的,a、b、c是未知的,而且是线性的,所以可以被非常简单的拟合出来。假设a=2.5,b=0.5,c=-1,加入随机扰动。拟合的最终效果为:

最终得到的拟合参数为:a=2.47,b=0.47,c=-0.66。考虑到随机噪声的影响,与原始数据相差不大,源代码如下:

%线性拟合x=0:0.5:10;a=2.5;b=0.5;c=-1;%原函数y=a*sin(0.2*x.^2+x)+b*sqrt(x+1)+c+0.5*rand(size(x));%拟合部分yn1=sin(0.2*x.^2+x);yn2=sqrt(x+1);yn3=ones(size(x));X=[yn1;yn2;yn3];%拟合,得到系数Pn=X'\y';yn=Pn(1)*yn1+Pn(2)*yn2+Pn(3)*yn3;%绘图figure()subplot(1,2,1)plot(x,y,'linewidth',1.5,'MarkerSize',15,'Marker','.','color','k')legend('原始数据点','Location','southoutside','Orientation','horizontal')legend('boxoff')subplot(1,2,2)hold onx2=0:0.01:10;plot(x2,Pn(1)*sin(0.2*x2.^2+x2)+Pn(2)*sqrt(x2+1)+Pn(3),'-','linewidth',1.5,'color','r')plot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')hold offbox onlegend('拟合曲线','数据点','Location','southoutside','Orientation','horizontal')legend('boxoff')

事实上,其实常用的拟合方式中,有很多都是线性拟合,比如多项式拟合,傅里叶拟合等。对于多项式拟合,只需要把拟合的部分替换成x,x.^2,x.^3这种形式。对于傅里叶级数,只需要把拟合的部分替换为sin(x),sin(2*x),sin(3*x)这种形式就行。

下面展示一个利用线性拟合,进行不同频率的三角波级数拟合正弦函数的例子:

clear;clc;close allt=0:0.001:2*pi;%原函数YS=sin(t);%基函数N=21;Yo=[];for k=1:NYn=sawtooth(k*(t+pi/2),0.5);Yo=[Yo,Yn'];endYS=YS';%拟合a = Yo\YS;%绘图figure()for k=1:Nclfplot(t,Yo(:,1:k)*a(1:k,:),t,YS,'LineWidth',1)ylim([-1.3,1.3])xlim([0,6.3])pause(0.1)end

2 一维非线性拟合

2.1 简单的非线性拟合

前面介绍了线性的拟合方法。如果一个非线性方程,可以化为上面线性方程中公式中给出的样子,那么我们也可以套用线性的方法去求解。

比如下面的方程:

y=a∗exp⁡(−bx)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)

这样求出来之后,再反变换回去,就可以得到原方程的系数了。

代码如下:

clearclcclose all%方法1x=0:0.5:10;a=2.5;b=0.5;%原函数y=a*exp(-b*x);y=y+0.5*y.*rand(size(y));%加噪声%变形函数%Lg_y=Lg_a+b*(-x)Lg_y=log(y);%拟合部分yn1=ones(size(x));yn2=-x;X=[yn1;yn2];%拟合,得到系数Pn=X'\Lg_y';%反变换a_fit=exp(Pn(1));b_fit=Pn(2);%绘图figure()hold onx2=0:0.01:10;plot(x2,a_fit*exp(-b_fit*x2),'-','linewidth',1.5,'color','r')plot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')hold off

2.2 matlab中Curve Fitting App

matlab自带了一个Curve Fitting App,它在matlab集成的App里面。

界面左边输入数据,右侧选择方法。

界面里常用的拟合方式都有,而且直接展示拟合效果,非常方便。非常适合鼠标直接拖拖拽拽点点点的操作方式。

2.3 matlab中非线性拟合的实现

2.3.1 fit()函数

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,一个简单的高斯函数形式的非线性方程,其参数给定为:

在已知函数形式,求解这四个参数条件下,6种不同的函数非拟合效果如下:

可以看到,这几种方法都能够较好的拟合出想要的结果。

第二个例子是一个指数增长的正弦函数,在很多线性系统中都可以测量到这种信号。函数的形式为:

y=a*x+b*sin(c*x).*exp(d*x)+e 。其给定的参数为:

这个函数的拟合具有一定难度,拟合过程中会遇到非常多的局部解。

由于初始点比较随机,所以最终结果每次都会不一样。

其中前面的几种方法对于初始值的敏感度比较高,如果初始值选的比较接近原始解,也是可以得到较好的结果。其中nlinfit函数经常会报错,容错率较低。而PSO算法经常能够收敛到最优解(虽然不是每次都可以,偶尔也会陷入局部解)。

代码如下:

clearclc%函数大比拼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 Squaresft = 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 onplot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')plot(x,y1,'-','linewidth',1.5,'color','r')hold offbox onylim(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 onplot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')plot(x,y2,'-','linewidth',1.5,'color','r')hold offbox onylim(y_lim)title('nlinfit函数')%% 3 lsqnonlin()函数 trust-region-reflectivemodelfun = @(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 onplot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')plot(x,y3,'-','linewidth',1.5,'color','r')hold offbox onylim(y_lim)title('lsqnonlin函数')%% 4 lsqcurvefit()函数 trust-region-reflectivemodelfun = @(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 onplot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')plot(x,y4,'-','linewidth',1.5,'color','r')hold offbox onylim(y_lim)title('lsqcurvefit函数')%% 5 fsolve()函数 %默认算法trust-region-doglegmodelfun = @(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 onplot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')plot(x,y5,'-','linewidth',1.5,'color','r')hold offbox onylim(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 onplot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')plot(x,y6,'-','linewidth',1.5,'color','r')hold offbox onylim(y_lim)title('PSO算法')

下图展示了PSO求解过程中逐渐收敛到全局最优解的过程。

3 高维非线性方程组拟合

3.1 一般形式高维方程或方程组的拟合

之前的文章中的数据具有一 一对应的特点,所以严格来讲并不是普遍的二维拟合。对于一些复杂的二维函数,比如椭圆,可能原本的拟合就会失效。

对于这种一般性质的方程或方程组,将原本输入方程整理为f(x,y,z,…)=0的形式。衡量拟合程度的优化函数,就直接取函数f(xi,yi,zi,…)的值即可。

下面演示最终的两个例子:

第一个是三维直线,采用两平面式描述。

Ax+By+Cz-1=0Dx+Ey+Fz-1=0

总共2个方程,维度为3维,第一个方程有3个参数,第二个方程也有3个参数。离散点已知的条件下,三维直线的两平面表达式不唯一。

最终拟合效果如下:

第二个是二维椭圆,方程为:

x^2+Axy+By^2+Cx+Dy+E=0

总共1个方程,维度为2维。方程共有5个参数。

最终拟合效果如下:

clearclcclose all%% 演示1%1 导入数据(这里用的是人工生成的数据)%三维直线拟合,函数表示%1.0*x+1.9*y+3.0*z=1;%1.2*x-0.4*y-1.7*z=1;x=0:0.04:1;%求解出y和z% [ 1.0, 3.0 [ y [1.0% -0.4,-1.7] * z] = 1 - 1.2] xy=zeros(size(x));z=y;for k=1:length(x)R=[1.9,3.0;-0.4,-1.7]\[1-1.0*x(k);1-1.2*x(k)];[y(k),z(k)]=Value(R);endrand_n=randperm(length(x));%生成随机的初始输入数据x1=x(rand_n)+0.05*randn(size(x));y1=y(rand_n)+0.05*randn(size(x));z1=z(rand_n)+0.05*randn(size(x));%2 整理格式,将数据和要拟合的函数进行格式整理%准备数据XX={x1,y1,z1};%准备函数F1=@(p1,XX1) p1(1)*XX1{1}+p1(2)*XX1{2}+p1(3)*XX1{3}-1;F2=@(p2,XX2) p2(1)*XX2{1}+p2(2)*XX2{2}+p2(3)*XX2{3}-1;FF={F1,F2};%3 生成最终优化函数,带入到优化方程中求解fun=@(p) Dis(p,{3,3},XX,FF);%p参数,{3,3}第一个方程3个参数,第二个方程3个参数。XX离散点。FF函数表达式。OP=optimoptions('particleswarm','FunctionTolerance',0);[p,fval,~,~]=particleswarm(fun,6,[-5,-5,-5,-5,-5,-5],[5,5,5,5,5,5],OP);%4 对比拟合效果figure()x2=0:0.01:1;y2=zeros(size(x2));z2=y2;for k=1:length(x2)R=[p(2),p(3);p(5),p(6)]\[1-p(1)*x2(k);1-p(4)*x2(k)];[y2(k),z2(k)]=Value(R);end%系数可能不同。因为直线的两平面表示不唯一hold onplot3(x2,y2,z2)plot3(x1,y1,z1,'*');hold offview(3)%% 演示2%1 导入数据(这里用的是人工生成的数据)%二维椭圆拟合th=0:0.15:2*pi;a=3.2;%椭圆轴1b=4.8;%椭圆轴2x0=-1.9;y0=-4.1;beta=1.1;%椭圆旋转角度%绘制椭圆x=a*cos(th);y=b*sin(th);%旋转beta角度x_r=x*cos(beta)-y*sin(beta);y_r=x*sin(beta)+y*cos(beta);rand_n=randperm(length(x));%生成随机的初始输入数据x1=x_r(rand_n)+0.1*randn(size(x));y1=y_r(rand_n)+0.1*randn(size(x));%2 整理格式,将数据和要拟合的函数进行格式整理%准备数据XX={x1,y1};%准备函数F1=@(p1,XX1) XX1{1}.*XX1{1}+p1(1)*XX1{1}.*XX1{2}+p1(2)*XX1{2}.*XX1{2}+p1(3)*XX1{1}+p1(4)*XX1{2}+p1(5);FF={F1};%3 生成最终优化函数,带入到优化方程中求解fun=@(p) Dis(p,{5},XX,FF);OP=optimoptions('particleswarm','FunctionTolerance',0);[p,fval,~,~]=particleswarm(fun,5,[-20,-20,-20,-20,-20],[20,20,20,20,20],OP);%4 对比拟合效果figure()hold onplot(x1,y1,'*')Fun_Plot=@(xp,yp) xp.*xp+p(1)*xp.*yp+p(2)*yp.*yp+p(3)*xp+p(4)*yp+p(5);fimplicit(Fun_Plot,[-6 6 -6 6],'LineWidth',1)hold off%% 用到的函数function varargout=Value(f)%多元素赋值,例子:%[a,b,c]=Value([1,2,3]);%多变量赋值%[xy(1),xy(2),t]=Value([2,5,3]);%复合赋值%[b,a]=Value([a,b]);%元素交换%来源:hyhhhyh21,matlab爱好者N=numel(f);varargout=cell(1,N);if nargout~=Nerror('输入输出数量不一致')endfor k=1:Nvarargout{k}=f(k);endendfunction Sum_N=Dis(p,p_num,XX,FF)%用函数值评价拟合程度%输入:要拟合的参数p%输入:p_num cell格式,每个方程的参数数量%输入:XX 数据,以cell形式储存%输入:FF 拟合函数,以cell形式储存N_F=numel(FF);%要联立的方程数量L=length(XX{1});%离散点的数量N_L=numel(XX);%离散点维度Sum_N=0;XXm=cell(1,N_L);%直接计算函数的值p_n=1;%参数的索引for k=1:N_F%对每一个方程进行计算FF_k=FF{k};%方程F_p=p_num{k};%该方程用到的参数数量for m=1:Lfor n=1:N_LXXm{n}=XX{n}(m);endDistance=FF_k(p(p_n:(p_n+F_p-1)),XXm);Sum_N=Sum_N+Distance^2;endp_n=p_n+F_p;%更新函数索引endend

3.2 一般形式参数方程的拟合

高维参数方程的拟合比较困难。因为原本方程中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的值不唯一,但是它们的比例肯定唯一。

第二个例子是一个三维旋转曲面。参数方程为:

x= A*u.*sin(v+B)y=-C*u.*cos(v+D)z=v

方程为三维参数方程,有2个参数u、v。第一个方程有2个未知量AB,第二个方程有2个未知量CD,第三方程有0个未知量。

最终拟合效果如下:

这两个例子的演示代码如下。这里参数方程的Dis_P()由于频繁的计算点与点之间的距离,所以运算速度比第一章单纯函数的Dis()较慢。

clearclcclose all%% 演示1%三维李萨如图形拟合%1 导入数据(这里用的是人工生成的数据)t=0:0.01:4*pi;x=sin(4*t);y=cos(5*t);z=sin(6*t);rand_n=randperm(length(t));x1=x(rand_n)+0.02*randn(size(t));y1=y(rand_n)+0.02*randn(size(t));z1=z(rand_n)+0.02*randn(size(t));%2 整理格式,将数据和要拟合的函数进行格式整理%输入参数方程XX={x1,y1,z1};%离散点输入F1=@(p1,uu1) sin(p1(1)*uu1{1});F2=@(p2,uu1) cos(p2(1)*uu1{1});F3=@(p3,uu1) sin(p3(1)*uu1{1});FF={F1,F2,F3};%方程输入u1=0:0.05:13;%设置参数的最大最小范围以及精度,能达到绘图精度即可uu={u1};%参数输入%3 生成最终优化函数,带入到优化方程中求解fun=@(p) Dis_P(p,{1,1,1},uu,XX,FF);%使得DisP函数输出的Sum_N最小%p参数,{1,1,1}代表有3个方程,每个方程的代求参数p个数为1。uu为参数输入。XX为离散点输入。FF为方程输入。OP=optimoptions('particleswarm','FunctionTolerance',0,'InertiaRange',[0.3,1.2],'MaxStallIterations',200);[pp,fval,~,~]=particleswarm(fun,3,[1,1,1],[10,10,10],OP);%4 对比拟合效果figure()hold ontt=0:0.01:4*pi;%pp=pp/pp(1)*4;%这里不一定必须是4,5,6,只需要比例为4:5:6就行[a2,b2,c2]=Value(pp);x2=sin(a2*tt);y2=cos(b2*tt);z2=sin(c2*tt);plot3(x2,y2,z2);plot3(x1,y1,z1,'*');hold offview(3)%% 演示2%三维螺旋面拟合%1 导入数据(这里用的是人工生成的数据)F1=@(p1,uu1) p1(1).*uu1{1}.*sin(uu1{2}+p1(2));F2=@(p2,uu1) -p2(1).*uu1{1}.*cos(uu1{2}+p2(2));F3=@(p3,uu1) uu1{2};u_rand=rand_AB([-4,4],100,1);v_rand=rand_AB([-5,5],100,1);x=F1([2,0.1],{u_rand,v_rand});y=F2([1,0.1],{u_rand,v_rand});z=F3([],{u_rand,v_rand});rand_n=randperm(length(x));x1=x(rand_n)+0.01*randn(size(x));y1=y(rand_n)+0.01*randn(size(x));z1=z(rand_n)+0.01*randn(size(x));%2 整理格式,将数据和要拟合的函数进行格式整理%输入参数方程XX={x1,y1,z1};%离散点输入FF={F1,F2,F3};%方程输入u1=-4:0.1:4;%设置参数的最大最小范围以及精度,能达到绘图精度即可v1=-5:0.1:5;%设置参数的最大最小范围以及精度,能达到绘图精度即可uu={u1,v1};%参数输入%3 生成最终优化函数,带入到优化方程中求解fun=@(p) Dis_P(p,{2,2,0},uu,XX,FF);%使得DisP函数输出的Sum_N最小OP=optimoptions('particleswarm','FunctionTolerance',0,'InertiaRange',[0.3,1.2],'MaxStallIterations',200);[pp,fval,~,~]=particleswarm(fun,4,[0.1,0,0.1,0],[10,2*pi,10,2*pi],OP);%4 对比拟合效果figure()hold onplot3(x1,y1,z1,'*');funx = @(u,v) pp(1)*u.*sin(v+pp(2));funy = @(u,v) -pp(3)*u.*cos(v+pp(4));funz = @(u,v) v;fsurf(funx,funy,funz,[-4 4 -5 5],'LineStyle','none','FaceAlpha',0.5)xlim([-8,8]);ylim([-8,8]);zlim([-5,5]);colormap(hsv)camlightplot3([0,8],[0,0],[0,0],'linewidth',2,'color','k')plot3([0,0],[0,8],[0,0],'linewidth',2,'color','k')plot3([0,0],[0,0],[0,5],'linewidth',2,'color','k')hold offview(3)function Sum_N=Dis_P(p,p_num,uu,XX,FF)%用距离曲线的距离评价拟合程度(参数方程)%输入:p 要拟合的参数%输入:p_num 数组,每个方程的参数数量%输入:uu 参数方程中的参数,以cell形式储存%输入:XX 数据,以cell形式储存%输入:FF 拟合函数,以cell形式储存N_F=numel(FF);%要联立的方程数量L=length(XX{1});%离散点的数量N_L=numel(XX);%拟合参数p的数量N_u=numel(uu);%参数方程中参数的数量if N_u>1uu_new=ndgrid_h(uu{:});uu=uu_new;end%循环每个点,求到直线的距离%在假定参数p的情况下,计算假定函数Point_NF=cell(N_F,1);p_n=1;%参数的索引for k=1:N_FF_p=p_num{k};%该方程用到的参数数量Point_NF{k}=FF{k}(p(p_n:(p_n+F_p-1)),uu);%计算假定函数各个点坐标p_n=p_n+F_p;%更新函数索引endSum_N=0;for k=1:L%分别求每个假定函数的点,与真实离散点之间距离的平方和Point_distance2=0;for m=1:N_F%读取真实点坐标Point_distance2=Point_distance2+(Point_NF{m}-XX{m}(k)).^2;endMin_distance2=min(Point_distance2);%求出最小距离,即为点与假定函数之间的距离Sum_N=Sum_N+Min_distance2;endendfunction varargout=Value(f)%多元素赋值,例子:%[a,b,c]=Value([1,2,3]);%多变量赋值%[xy(1),xy(2),t]=Value([2,5,3]);%复合赋值%[b,a]=Value([a,b]);%元素交换%来源:hyhhhyh21,matlab爱好者N=numel(f);varargout=cell(1,N);if nargout~=Nerror('输入输出数量不一致')endfor k=1:Nvarargout{k}=f(k);endendfunction X=rand_AB(AB,varargin)%生成区间[A,B]之间的随机分布%除了AB之外,其余输入与rand相同[A,B]=Value(AB);X=rand(varargin{1:end});X=A+X*(B-A);endfunction S=ndgrid_h(varargin)%来源于matlab自带的源代码。%用法和ndgrid用法一样,但是将输出更改为cellN=nargin;S=cell(1,N);if N==1S{1}=varargin;else j = 1:N;siz = cellfun(@numel,varargin);if N == 2 % Optimized Case for 2 dimensionsx = full(varargin{1}(:));y = full(varargin{2}(:)).';S{1} = repmat(x,size(y));S{2} = repmat(y,size(x));elsefor i=1:Nx = full(varargin{j(i)});s = ones(1,N);s(i) = numel(x);x = reshape(x,s);s = siz;s(i) = 1;S{i} = repmat(x,s);endendendS2=cell(1,N);for k=1:NS2{k}=S{k}(:);endS=S2;end

主要函数就是Dis和Dis_P,准备工作就是把方程和离散点整理成可以输入的形式。优化用到的函数就是PSO(particleswarm),需要更改未知参数数量和范围就可以。

3.3 补充

这一部分的拟合优度或者置信区间就没办法拿出来说了,毕竟有试凑的嫌疑。适合单纯的验证。

另外这里和最小二乘求出来的结果会有所不同。这主要是定义的问题。其中最小二乘指的是距离y方向上最小的二乘距离。而这一章节中定义的最小二乘,是指的点到拟合曲线距离的最小二乘(有点类似于主成分分析)。虽然一般情况下两者并不会有太大差别,但是如果有误差的话,还请注意这一点,

如下图,这里是一个椭圆分布的离散点。

如果将方程的形式定义为y=ax+b,要求距离y最小,则最终结果为红色的线。如果将方程定义为ax+by+c=0,要求距离直线的距离最小,则最终拟合结果为黄色的线。这里无所谓对错,只是要求的最小定义不同。

单独在后面补充了这么一个章节,希望在应用中注意。

上面的图,绘图代码如下:

%构建原始数据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=0FF={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 onplot(x,y,'.')plot(x,yn,'linewidth',4)plot(x,(p(1)*x+p(3))/p(2),'linewidth',4)hold offlegend({'原始数据','最小二乘','非线性拟合'})

参考文献

信赖域法 /p/205555114

本内容不代表本网观点和政治立场,如有侵犯你的权益请联系我们处理。
网友评论
网友评论仅供其表达个人看法,并不表明网站立场。