Использование Lsq-Fit в Юлии - PullRequest
Новая
       25

Использование Lsq-Fit в Юлии

1 голос
/ 15 марта 2020

Я пытаюсь попрактиковаться в подборе функции Lsq-Fit в Джулии. Производная распределения Коши с параметрами \ gamma и x_0. Следуя этому руководству, я попробовал 

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[...], оставляя меня в замешательстве. Может кто-нибудь сказать мне, что я делаю не так?

1 Ответ

1 голос
/ 15 марта 2020

Есть несколько проблем с тем, что вы написали:

  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)

    Один из способов исправить это - использовать точечную запись для широковещания всех операторов, чтобы они работали поэлементно: 

    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

  2. Другая проблема заключается в том, что параметры, которые вы ищете, должны быть значениями с плавающей точкой. Но ваше первоначальное предположение p0 инициализируется целыми числами, что сбивает с толку curve_fit. Есть два способа исправить это. Либо поместите значения с плавающей точкой в ​​p0: 

    julia> p0 = [3349.0, 49.0]
2-element Array{Float64,1}:
 3349.0
   49.0

    или используйте типизированный инициализатор массива , чтобы явно указать тип элемента: 

    julia> p0 = Float64[3349, 49]
2-element Array{Float64,1}:
 3349.0
   49.0

  3. Это не совсем ошибка, но я бы посчитал более интуитивно понятным вычислять a/b вместо a*b^(-1). Кроме того, yarr можно вычислить с помощью простой трансляции, используя точечную запись вместо понимания.

Оборачивая все это вместе: 

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)

выход: 

julia> fit.param
2-element Array{Float64,1}:
 3349.999986535933  
   49.99999203625603
...