Can't fit dose response curve, could you tell me Why and How

General Usage
,
1

Hello there !

I just have started programming and Julia these days, and I am trying fitting curve (dose response) by Optim.jl
But I can’t obtain better parameters and predictions (the all predictions have same value eventhough I use parameters that are calculated by optimize()) .
So, I would like you to tell me which is/are wrong and how to fix it.
I just have started programming and Julia

Here is my codes;

xdata = [0.1, 0.3, 1.0, 3.0, 10.0]
ydata = [138.48, 150.23, 192.14, 199.78, 233.76]
xbase = collect(range(minimum(xdata), maximum(xdata), 100))

@. DoseRes(x, p) = p[1] + ((p[2] - p[1]) / (1 + 10 ^ ((p[3] - x) * p[4])))

function cost(x, y, p)
    err = 0.0
    for i in 1:length(x)
        err = sum(abs2(y[i] - DoseRes(x[i], p...))) / (2 * length(x))
    end
    return err
end
result = optimize(p -> cost(xdata, ydata, p), [0.0, 0.0, 0.0, 0.0]) 
pbest = Optim.minimizer(result)

DoseRes(xbase, pbest)

result is like this;

100-element Vector{Float64}:
 233.75995328637563
 233.75995328637563
 233.75995328637563
 233.75995328637563
 233.75995328637563
 233.75995328637563
 233.75995328637563
 233.75995328637563
 233.75995328637563
 233.75995328637563
 233.75995328637563
 233.75995328637563
 233.75995328637563
   ⋮
2

Wow, over 12 months and no reply.

You need a better starting guess. Most of the optimisers take your starting guess and roll down hill from there. I think in your example, the error surface is very flat and there is no where for the optimiser to go.

[ydata[1], 1.0, 1.0, 1.0]

produced a different result (I didn’t check if it was “good”). For this sort of relatively simple problem you should develop a heuristic to determine an initial guess.

This is not a Julia problem, the optimisation libraries in other languages will behave the same.

3

Apart from the bad initial guess, there’s also an error in the definition of the cost function: every pass of the for loop takes the “sum” of only one value and overwrites the previous err with the result. So in the end cost returns the error for the last point only. Then it’s not surprising that the optimization returns a flat function.

Here’s a version that seems to work:

using CairoMakie
using Optim

xdata = [0.1, 0.3, 1.0, 3.0, 10.0]
ydata = [138.48, 150.23, 192.14, 199.78, 233.76]
xbase = range(minimum(xdata), maximum(xdata), 100)

DoseRes(x; p) = p[1] + ((p[2] - p[1]) / (1 + 10 ^ ((p[3] - x) * p[4])))

cost(x, y, p) = sum(abs2.(y - DoseRes.(x; p))) / (2 * length(x))

p0 = [100, 240, 0.0, 1.0]
result = optimize(p -> cost(xdata, ydata, p), p0) 
pbest = Optim.minimizer(result)
@show cost(xdata, ydata, p0)
@show cost(xdata, ydata, pbest)
@show pbest

fig = scatter(xdata, ydata, figure=(; resolution=(600, 400)))
lines!(xbase, DoseRes.(xbase, p=p0), label="Initial guess")
lines!(xbase, DoseRes.(xbase, p=pbest), label="Optimized")
axislegend(position=:rb, framevisible=false)
fig

Output:

cost(xdata, ydata, p0) = 629.4901852908467
cost(xdata, ydata, pbest) = 42.599002665390294
pbest = [-14075.614358714607, 227.81615863795133, -8.098909537951354, 0.2707927218725344]

image

1 Like