在Julia中使用Lsq-Fit

问题描述 投票:1回答:1

我正在尝试在Julia中练习Lsq-Fit-fit的拟合。参数\ gamma和x_0的柯西分布的导数。按照我尝试过的this手册

f(x, x_0, γ) = -2*(x - x_0)*(π * γ^3 * (1 + ((x - x_0)/γ)^2)^2)^(-1)
x_0 = 3350
γ = 50
xarr = range(3000, length = 5000, stop = 4000)
yarr = [f(x, x_0, γ) for x in xarr]

using LsqFit
# p ≡ [x_0, γ]
model(x, p) = -2*(x - p[1])*(π * (p[2])^3 * (1 + ((x - p[1])/p[2])^2)^2)^(-1)
p0 = [3349, 49]
curve_fit(model, xarr, yarr, p0)
param = fit.param

...它不起作用,给出MethodError: no method matching -(::StepRangeLen[...],让我感到困惑。可以告诉我我在做什么错吗?

julia curve-fitting
1个回答
1
投票

您写的内容有一些问题:

  1. model函数的第一个参数(x)是自变量的完整向量,而不仅仅是一个值。这是您提到的错误来源:

    julia> model(x, p) = -2*(x - p[1])*(π * (p[2])^3 * (1 + ((x - p[1])/p[2])^2)^2)^(-1);
    julia> p0 = [3349, 49];
    julia> model(xarr, p0);
    ERROR: MethodError: no method matching -(::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}, ::Float64)
    

    解决此问题的一种方法是对所有运算符使用点符号broadcast,以便它们逐个元素地工作:

    julia> model(x, p) = -2*(x .- p[1]) ./ (π * (p[2])^3 * (1 .+ ((x .- p[1])/p[2]).^2).^2);
    julia> model(xarr, p0); # => No error
    

    但是如果这太繁琐,您可以让@.宏为您完成工作:

    @.
  2. 另一个问题是,您要查找的参数是浮点值。但是您最初的猜测# just put @. in front of the expression to transform every # occurrence of a-b into a.-b (and likewise for all operators) # which means to compute the operation elementwise julia> model(x, p) = @. -2*(x - p[1])*(π * (p[2])^3 * (1 + ((x - p[1])/p[2])^2)^2)^(-1); julia> model(xarr, p0); # => No error 用整数初始化,这使p0感到困惑。有两种解决方法。可以将浮点值放在curve_fit

    p0

    或使用julia> p0 = [3349.0, 49.0] 2-element Array{Float64,1}: 3349.0 49.0 明确指定元素类型:

    typed array initializer
  3. 这并不是真正的错误,但是我发现计算julia> p0 = Float64[3349, 49] 2-element Array{Float64,1}: 3349.0 49.0 而不是a/b更直观。同样,可以使用a*b^(-1)而不是理解来通过简单的广播来计算yarr

将所有这些包装在一起:

dot notation

产量:

f(x, x_0, γ) = -2*(x - x_0)*(π * γ^3 * (1 + ((x - x_0)/γ)^2)^2)^(-1)
(x_0, γ) = (3350, 50)

xarr = range(3000, length = 5000, stop = 4000);
# use dot-notation to "broadcast" f and map it
# elementwise to elements of xarr
yarr = f.(xarr, x_0, γ);

using LsqFit
model(x, p) = @. -2*(x - p[1]) / (π * (p[2])^3 * (1 + ((x - p[1])/p[2])^2)^2)
p0 = Float64[3300, 10]

fit = curve_fit(model, xarr, yarr, p0)
© www.soinside.com 2019 - 2024. All rights reserved.