# 利用matlab实现非线性拟合（三维、高维、参数方程）

### 利用matlab实现非线性拟合[三维、高维、参数方程]

2021.06 更新。补充了非线性拟合中，关于最小二乘定义的问题，放在了最后一章。
2021.12更新。应评论区建议，补充了3.3章绘图用的代码。

# 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()函数进行多项式拟合。下面举一个小例子：

``````y=0.03 x^4 - 0.5 x^3 + 2 x^2 - 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 on
plot(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 off
legend('原始数据点','理论曲线','Location','southoutside','Orientation','horizontal')
legend('boxoff')
box on
subplot(1,2,2)
hold on
plot(x2,y2,'-','linewidth',1.5,'color','r')
plot(x,y,'LineStyle','none','MarkerSize',15,'Marker','.','color','k')
hold off
box on
legend('拟合曲线','数据点','Location','southoutside','Orientation','horizontal')
legend('boxoff')
``````

## 1.2 线性拟合

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)+...

``````%线性拟合
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 on
x2=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 off
box on
legend('拟合曲线','数据点','Location','southoutside','Orientation','horizontal')
legend('boxoff')
``````

``````clear;
clc;
close all
t=0:0.001:2*pi;%原函数
YS=sin(t);%基函数
N=21;
Yo=[];
for k=1:N
Yn=sawtooth(k*(t+pi/2),0.5);
Yo=[Yo,Yn'];
end
YS=YS';%拟合
a = Yo\YS;
%绘图
figure()
for k=1:N
clf
plot(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 ⁡ ( − b x ) y=a*\exp(-bx) y=aexp(bx)

lg ⁡ ( y ) = lg ⁡ ( a ) + b ∗ ( − x ) \lg(y)=\lg(a)+b*(-x) lg(y)=lg(a)+b(x)

``````clear
clc
close all

%方法1
x=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 on
x2=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函数来实现的拟合。

### 2.3.3 lsqnonlin()函数和lsqcurvefit()函数

lsqnonlin()也是matlab中自带的一个非线性拟合函数。它给出了两种计算非线性拟合的方法，一种是比较经典之前介绍过的LM方法，一种是信赖域方法。信赖域法（trust region reflective）是通过Hessian 矩阵，逐步试探邻域内的最小化，来求解问题的。相比较之前的那些雅克比相关的方法，信赖域法会占用更多内存和速度，所以适用于中小规模的矩阵。

lsqcurvefit()函数和lsqnonlin()内容上相似，只是引用格式上有所不同。

### 2.3.6 不同算法的对比效果

abcd
3.82.14.4-1.3

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算法')
``````

# 3 高维非线性方程组拟合

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

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

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

``````clear
clc
close 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] x
y=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);
end
rand_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 on
plot3(x2,y2,z2)
plot3(x1,y1,z1,'*');
hold off
view(3)

%% 演示2
%1 导入数据（这里用的是人工生成的数据）
%二维椭圆拟合
th=0:0.15:2*pi;
a=3.2;%椭圆轴1
b=4.8;%椭圆轴2
x0=-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 on
plot(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~=N
error('输入输出数量不一致')
end
for k=1:N
varargout{k}=f(k);
end
end

function 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:L
for n=1:N_L
XXm{n}=XX{n}(m);
end
Distance=FF_k(p(p_n:(p_n+F_p-1)),XXm);
Sum_N=Sum_N+Distance^2;
end
p_n=p_n+F_p;%更新函数索引
end
end
``````

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

``````x=sin(A*u)
y=cos(B*u)
z=sin(C*u)
``````

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

``````clear
clc
close 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 on
tt=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 off
view(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 on
plot3(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)
camlight
plot3([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 off
view(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>1
uu_new=ndgrid_h(uu{:});
uu=uu_new;
end
%循环每个点，求到直线的距离
%在假定参数p的情况下，计算假定函数
Point_NF=cell(N_F,1);
p_n=1;%参数的索引
for k=1:N_F
F_p=p_num{k};%该方程用到的参数数量
Point_NF{k}=FF{k}(p(p_n:(p_n+F_p-1)),uu);%计算假定函数各个点坐标
p_n=p_n+F_p;%更新函数索引
end
Sum_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;
end
Min_distance2=min(Point_distance2);%求出最小距离，即为点与假定函数之间的距离
Sum_N=Sum_N+Min_distance2;
end
end

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~=N
error('输入输出数量不一致')
end
for k=1:N
varargout{k}=f(k);
end
end

function X=rand_AB(AB,varargin)
%生成区间[A,B]之间的随机分布
%除了AB之外，其余输入与rand相同
[A,B]=Value(AB);
X=rand(varargin{1:end});
X=A+X*(B-A);
end

function S=ndgrid_h(varargin)
%来源于matlab自带的源代码。
%用法和ndgrid用法一样，但是将输出更改为cell
N=nargin;
S=cell(1,N);
if N==1
S{1}=varargin;
else
j = 1:N;
siz = cellfun(@numel,varargin);
if N == 2 % Optimized Case for 2 dimensions
x = full(varargin{1}(:));
y = full(varargin{2}(:)).';
S{1} = repmat(x,size(y));
S{2} = repmat(y,size(x));
else
for i=1:N
x = 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);
end
end
end
S2=cell(1,N);
for k=1:N
S2{k}=S{k}(:);
end
S=S2;
end
``````

## 3.3 补充

``````%构建原始数据
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({'原始数据','最小二乘','非线性拟合'})
``````

# 参考文献

原文作者：hyhhyh21
原文地址: https://blog.csdn.net/weixin_42943114/article/details/116803379
本文转自网络文章，转载此文章仅为分享知识，如有侵权，请联系博主进行删除。