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

**URL:** <https://discourse.julialang.org/t/cant-fit-dose-response-curve-could-you-tell-me-why-and-how/79996>\
**Category:** General Usage\
**Tags:** question, optim\
**Created:** [April 25, 2022, 3:33pm UTC](https://discourse.julialang.org/t/cant-fit-dose-response-curve-could-you-tell-me-why-and-how/79996 "2022-04-25T15:33:26Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![MEEMO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/meemo/32/35689_2.png) [@MEEMO](https://discourse.julialang.org/u/MEEMO)\
**Post date:** [April 25, 2022, 3:33pm UTC](https://discourse.julialang.org/t/cant-fit-dose-response-curve-could-you-tell-me-why-and-how/79996/1 "2022-04-25T15:33:26Z")

</div>

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;

```julia
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;

```julia
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
   ⋮

```

---

<div class="post-metadata">

**Author:** ![evan-wehi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evan-wehi/32/32692_2.png) [@evan-wehi](https://discourse.julialang.org/u/evan-wehi)\
**Post date:** [June 23, 2023, 12:38am UTC](https://discourse.julialang.org/t/cant-fit-dose-response-curve-could-you-tell-me-why-and-how/79996/2 "2023-06-23T00:38:31Z")

</div>

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.

```julia
[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.

---

<div class="post-metadata">

**Author:** ![sijo](https://avatars.discourse-cdn.com/v4/letter/s/da6949/32.png) [@sijo](https://discourse.julialang.org/u/sijo)\
**Post date:** [June 23, 2023, 8:43am UTC](https://discourse.julialang.org/t/cant-fit-dose-response-curve-could-you-tell-me-why-and-how/79996/3 "2023-06-23T08:43:04Z")

</div>

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:

```julia
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:

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

```

![image](https://global.discourse-cdn.com/julialang/original/3X/b/5/b56ad30092f47dd918f6b6d6ff8d21f917c1e8a0.png)