我正在尝试在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[...]
,让我感到困惑。可以告诉我我在做什么错吗?
您写的内容有一些问题:
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
但是如果这太繁琐,您可以让@.
宏为您完成工作:
@.
另一个问题是,您要查找的参数是浮点值。但是您最初的猜测# 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
这并不是真正的错误,但是我发现计算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)