利用matlab实现非线性拟合（补）

    之前在群里看有人问过三维拟合的问题。回去思考了一下，感觉和之前的非线性拟合还是有很多共同之处的。所以，这次将之前PSO方法的非线性拟合代码改动了一下，将其更改为适用性更广的高维拟合。

    没看过前面两篇文章的强烈建议回看一下。之前的一些应用背景和方法就不再重复说了。

利用matlab实现非线性拟合（上）

利用matlab实现非线性拟合（下）

1 高维方程或方程组拟合

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

    对于这种一般性质的方程或方程组，将原本输入方程整理为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个参数。

    最终拟合效果如下：

两个例子的代码如下：

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

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()较慢。

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

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

    临时补了之前没想到的一小点内容，可能这一部分的拟合优度或者置信区间就没办法拿出来说了，毕竟有试凑的嫌疑。

参考资料：

matlab官方文档

图片来源：hyhhyh21 由  在Pixabay上发布