| Julia | 您所在的位置：网站首页 › 函数拟合法 › Julia |
Julia
|
分享Julia中使用LsqFit对带参数的隐函数方程进行拟合！ 之前分享了LsqFit对一般带参数的单值函数进行非线性拟合； 对于带参数的隐函数方程，也可以通过LsqFit进行拟合，需要对隐函数方程稍作处理； 1. LsqFit安装及基本用法# 安装 using Pkg Pkg.add("LsqFit") # 基本用法 # fit = curve_fit(circle, xdata, ydata, p0)2. 隐函数处理与绘制隐函数图像时类似处理，考虑二维隐函数 f(x,y)=0 ，该图像为平面 xOy 上的一条曲线； 令 z=f(x,y) ，该三维函数为空间中曲面，其与 z=0 (即平面 xOy)的交线为隐函数 f(x,y)=0 的图像； 对于包含参数的二维隐函数方程，假设参数为p(p为向量)，则 f(x,y,p)=0 ，已有的数据为一系列散点 (x_i,y_i) 对应该隐函数方程，当构造三维曲面 z=f(x,y,p) ，可以得到原隐函数方程对应一系列散点在三维空间中表达为 (x_i,y_i,0) ； 用矩阵的表达，令 X=[x \ y] （ x 与 y 为列向量）， Z=[0,0,...,0]' ( Z 维度与 x 维度一致)，则 Z=f(X,p) 通过这种方式处理，将隐函数拟合问题处理为更高维度的单值函数非线性拟合问题； 3. 隐函数拟合以下为两个例子，分别为圆的方程和双纽线方程； 过程中需要绘制这两个隐函数图像示意图，可以参考这篇文章： 3.1 圆的方程： (x-x_0)^2+(y-y_0)^2=r^2此处在处理隐函数方程时，得益于Julia的多重派发机制，构造了两个版本的circle(x,p)函数，当x为矩阵时，用于拟合，当x为向量时，用于绘图； 构造Fake data时假设圆心(2, 0)，半径2 # (x-x0)^2+(y-y0)^2=r^2 using LsqFit,GLMakie # 定义隐函数，输入为矩阵，用于参数拟合 function circle(x::Matrix,p) x0, y0, r = p return @. (x[:,1]-x0)^2 + (x[:,2]-y0)^2 - r^2 end # 多重派发，输入为向量的版本，用于绘图 function circle(x::Vector,p) x0, y0, r = p return @. (x[1]-x0)^2 + (x[2]-y0)^2 - r^2 end # fack data # 圆心(2,0),半径2 # x∈[0,4] x = 4*rand(120) y = vcat(@.(sqrt(4-(x[1:60]-2)^2)), @.(-sqrt(4-(x[61:120]-2)^2))) + 0.02*randn(length(x)) xdata = [x y] # 拟合输入参数X # fitting p0 = [0.5, 0.5, 0.3] # 初始值，任意猜测值 ydata = zeros(size(xdata,1)) # 拟合输入参数Y fit = curve_fit(circle, xdata, ydata, p0) p = fit.param # 绘图 fig,ax = scatter(x,y,label="Source") ax.aspect=DataAspect() xfit = range(0, stop=4, length=100) yfit = range(0, stop=2.5, length=100) zfit = [circle([x,y], fit.param) for x in xfit,y in yfit] # 绘制拟合的隐函数方程 contour!(xfit, yfit, zfit, levels=[0], color=:black, label="fitting") axislegend("Data & Fitting") fig输出参数p的拟合结果如下： 3-element Vector{Float64}: 1.999933801267914 0.0019273621817362085 1.9995800140760118与构造Fake data时设定的[2, 0, 2](即圆心(2, 0),半径2)基本一致，绘图结果如下：  3.2 双纽线方程： (x^2+y^2)^2=2a^2(x^2-y^2) 与上述圆的方程处理基本一致，仍构造两个版本的函数func(x,a) 构造Fake data时，假设a=2 # 双扭线 (x^2+y^2)^2=2a^2(x^2-y^2) using LsqFit,GLMakie function func(x::Matrix,a) return @. (x[:,1]^2+x[:,2]^2)^2-2*a[1]^2*(x[:,1]^2-x[:,2]^2) end function func(x::Vector,a) return @. (x[1]^2+x[2]^2)^2-2*a[1]^2*(x[1]^2-x[2]^2) end # fake data t=[-3/4:0.03:3/4;pi.+(-3/4:0.03:3/4)] # 假设a=2.0构造fake data x=@. a*(2*cos(2*t))^0.5*cos(t) y=@.(a*(2*cos(2*t))^0.5*sin(t)) + 0.022*randn(length(t)) # fitting a0 = [1.0] # 初始参数，参数数量为1时，括号不可省略，否则函数func输出类型会发生变化，后续绘图将出错 xdata = [x y] ydata = zeros(size(xdata,1)) fit = curve_fit(func, xdata, ydata, a0) a=fit.param # Souce data fig,ax=scatter(x,y,markersize=14,label="Source") # fitting data xi = range(-3, stop=3, length=400) yi = range(-2, stop=2, length=400) zi = [func([x,y], a) for x in xi,y in yi] contour!(xi,yi,zi, levels=[0], color=:black,label="fitting") ax.aspect=DataAspect() axislegend("Data & Fitting",titlesize=20,labelsize=24) fig输出参数a的拟合结果如下： 1-element Vector{Float64}: 1.9999136125892245与构造Fake data时设定的a=2基本一致，绘图结果如下： 
|
【本文地址】