Dale Black (Sep 08 2023 at 18:38): Dale Black (Sep 08 2023 at 18:38):I am porting my labs code from Matlab to Julia and we have a Gamma Variate fitting function that I cannot seem to make work in Julia.
The matlab code is simple enough
gamma = @(gm,time)( ((aif_vec_new(end) - gm(3)) ./ (((time_vec_new(end)./gm(1)).^gm(2)).*(exp(gm(2).*(1-(time_vec_new(end)./gm(1))))))).*((time./gm(1)).^gm(2)).*(exp(gm(2).*(1-(time./gm(1))))) + gm(3) )
gamma_coef = lsqcurvefit(gamma,gm0,time_vec_new,aif_vec_new,gmlb,gmub);
And the Julia equivalent looks something like
xdata = [
0
1.1
2.2
3.3
4.4
5.5
6.6
7.7
8.8
9.9
11.0
12.1
13.2
14.3
]
ydata = [
122.8
115.2
92.4
103.6
161.4
156.2
196.9
188.0
286.1
244.4
196.1
363.6
420.8
395.7
]
p0 = [14.3, 0, 104]
lb = [14.3, -100, 0.0]
ub = [14.3, 100, 200]
function model(x, p; x_end = 14.3, y_end = 395)
r1 = (y_end - p[3]) / (((x_end ./ p[1]).^p[2]) * (exp(p[2] * (1 - (x_end ./ p[1])))))
r2 = ((x ./ p[1]).^p[2]) .* (exp.(p[2] .* (1 .- (x ./ p[1]))))
return r1 .* r2 .+ p[3]
end
gamma_coef = curve_fit(model, xdata, ydata, p0; lower = lb, upper = ub)
Where model = gamma, xdata = time_vec_new, ydata = aif_vec_new, p0 = gm0, lb = gmlb, and up = gmub
Unfortunately the Matlab code snippet fits the gamma function properly but the Julia curve_fit does not.
ArgumentError: matrix contains Infs or NaNs
chkfinite@lapack.jl:86[inlined]
var"#getrf!#1"(::Bool, ::typeof(LinearAlgebra.LAPACK.getrf!), ::Matrix{Float64})@lapack.jl:559
getrf!@lapack.jl:557[inlined]
#lu!#170@lu.jl:81[inlined]
lu!@lu.jl:80[inlined]
#lu#176@lu.jl:299[inlined]
lu@lu.jl:298[inlined]
\(::Matrix{Float64}, ::Vector{Float64})@generic.jl:1115
var"#levenberg_marquardt#10"(::Float64, ::Float64, ::Int64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Bool, ::Vector{Float64}, ::Vector{Float64}, ::Nothing, ::typeof(LsqFit.levenberg_marquardt), ::NLSolversBase.OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, ::Vector{Float64})@levenberg_marquardt.jl:126
levenberg_marquardt@levenberg_marquardt.jl:34[inlined]
#lmfit#15@curve_fit.jl:68[inlined]
lmfit@curve_fit.jl:67[inlined]
var"#lmfit#14"(::Symbol, ::Base.Pairs{Symbol, Vector{Float64}, Tuple{Symbol, Symbol}, NamedTuple{(:lower, :upper), Tuple{Vector{Float64}, Vector{Float64}}}}, ::typeof(LsqFit.lmfit), ::LsqFit.var"#18#20"{typeof(Main.var"workspace#118".model), Vector{Float64}, Vector{Float64}}, ::Vector{Float64}, ::Vector{Float64})@curve_fit.jl:64
lmfit@curve_fit.jl:46[inlined]
var"#curve_fit#16"(::Bool, ::Base.Pairs{Symbol, Vector{Float64}, Tuple{Symbol, Symbol}, NamedTuple{(:lower, :upper), Tuple{Vector{Float64}, Vector{Float64}}}}, ::typeof(LsqFit.curve_fit), ::typeof(Main.var"workspace#118".model), ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64})@curve_fit.jl:118
top-level scope@Local: 1[inlined]
The only thing that is potentially different between these two is the underlying mechanics of curve_fit vs lsqcurvefit from Matlab, as far as I can tell. Does anyone have familiarity with the LsqFit that might be able to point out the issue?
Dale Black (Sep 08 2023 at 18:57):I also found Gamma from Distributions.jl. I am not sure if this could be useful instead of the LsqFit approach
Dale Black (Sep 08 2023 at 19:35):Actually it looks like this is a simple enough approach that should work
data = DataFrame(X = xdata, Y = ydata)
gamma_model = glm(@formula(Y ~ X), data, Gamma(), LogLink())
gamma_fit = predict(gamma_model)
Edit, does not quite work
Dale Black (Sep 08 2023 at 22:52):I guess, what is the simplest way to fit a curve like this
Dale Black (Sep 08 2023 at 22:53):Where Bolus Tracking are the data points to be fit
Dale Black (Sep 08 2023 at 23:36):Kind of fixed it. The original problem was division by zero
begin
function gamma_model2(x, p; gm1 = 14.307, x_end = 14.307, y_end = 395)
r1 = (y_end - p[2]) / (((x_end / (gm1 + eps()))^p[1]) * exp(p[1] * (1 - x_end / (gm1 + eps()))))
r2 = ifelse.(x .== 0, 0, (x / (gm1 + eps())).^p[1]) .* exp.(p[1] .* (1 .- x / (gm1 + eps())))
return r1 .* r2 .+ p[2]
end
# Data points
x_data = [0.0, 0.35, 0.7, 1.05, 1.4, 1.75, 2.1, 2.45, 2.8, 3.15, 3.5, 3.85, 4.2, 14.307]
y_data = [122.8, 115.2, 92.4, 103.6, 161.4, 156.2, 196.9, 188.0, 286.1, 244.4, 196.1, 363.6, 420.8, 395.7]
# Initial guess and bounds for the remaining parameters
initial_guess = [0.0, 104.0]
lb = [-100.0, 0.0]
ub = [100.0, 200.0]
# Curve fitting
fit = curve_fit(gamma_model2, x_data, y_data, initial_guess, lower=lb, upper=ub)
gamma_fit2 = gamma_model2(x_data, fit.param)
end
And this is a little closer to the goal
image.png
Romeo Valentin (Feb 09 2024 at 16:46):Not quite your question, but I tried out different lsq solvers last year and found (the pretty new) NonlinearSolve.jl the best and most flexible. It's part of the SciML ecosystem so good bus factor.
Last updated: Nov 06 2024 at 04:40 UTC