Forcing one of the parameters to a specific value

[AlexisRenchon]

sorry to open 2 issues in a row, but they are separate so I thought it would be best this way.

I have a non-linear model with 6 parameters, and 2 independent variables, that you see here https://github.com/AlexisRenchon/DAMMmodel.jl/blob/master/src/DAMM.jl
that I want to fit to data using LsqFit.jl

My problem is that I often know the value of the parameter 5 (p[5]) before doing the fitting.

I have read before that doing

# porosity is my known value (Float64) of p[5]
lb = [0.0, 0.0, 0.0, 0.0, porosity, 0.0] # params can't be negative
ub = [Inf, Inf, Inf, Inf, porosity, Inf]

is not a good idea.

so what should I do instead?

[Gkreindler]

AlexisRenchon not sure if this is still relevant, but you can hack this as so: (not tested)

model_with_one_param_fixed = ((x,p) -> model(x, vcat(p[1:4], porosity, p[5])))

and then only estimate a 5-element vector omitting porosity

[AlexisRenchon]

Thanks @Gkreindler I will test this and let you know if it works, I think you meant p[6] and not p[5]

model_with_one_param_fixed = ((x,p) -> model(x, vcat(p[1:4], porosity, p[6])))

[Gkreindler]

@AlexisRenchon I mean p[5]. The idea is that you only want to estimate 5 parameters (original positions 1,2,3,4,6) because porosity is fixed. So you "skip" the 5th parameter in p, and the 6th parameter "moves" to the 5th place. Makes sense?

[AlexisRenchon]

you solution was correct @Gkreindler , also, @rafael-guerra-www wrote a nice MWE on Julia discourse https://discourse.julialang.org/t/mwe-fitting-a-model-with-a-known-parameter/73308/14

using LsqFit
@. model(t, p) = p[1] + p[2]*sin(2π*p[3]*t)*exp(-t*p[4])

tdata = LinRange(0,4,1000)
ydata = model(tdata, [1.0 2.0 30.0 0.5]) + 0.01*randn(length(tdata))
lb = [0, 0.0, 0.0, 0.0]
ub = [Inf, Inf, Inf,Inf]
p0 = [1.2, 1.2, 27.0, 0.1]
fit_bounds = curve_fit(model, tdata, ydata, p0, lower=lb, upper=ub)
fit_bounds.param    # 1.0003, 2.0001, 30.000, 0.4996

# Fix 3rd parameter (which disappears from constraints):
lb = [0, 0.0, 0.0]
ub = [Inf, Inf, Inf]
p0 = [1.2, 1.2, 0.1]
modelfixparam3 = (t,p) -> model(t, (p[1], p[2], 30.0, p[3]))
fit_bounds2 = curve_fit(modelfixparam3, tdata, ydata, p0, lower=lb, upper=ub)
fit_bounds2.param  # 1.0000, 2.0013, 0.5004