Hi,

I am currently exploring Pumas and I tried to estimate parameters based on an SDE model. While simulations work fine, I am not able to run parameter estimation successfully. I have to say that my understanding of SDEs is limited, so maybe it’s some obvious mistake.

This is my model definition:

```
p = ParamSet((θ = VectorDomain(3, lower = zeros(3)), Σ = RealDomain(lower = 0), Ω = RealDomain(lower = 0)))
function rfx(p)
ParamSet((η = Normal(0, p.Ω), ))
end
function pre_f(params, randoms, subject)
function (t)
Σ = params.Σ
θ = params.θ
η = randoms.η
(CL = θ[1]*exp(η[1]),
Vc = θ[2],
KA = θ[3],
Σ = Σ)
end
end
function f(du, u, p, t)
depot, cent = u
du[1] = -p.KA*depot
du[2] = p.KA*depot - (p.CL/p.Vc)*cent
end
function g(du, u, p, t)
du[1] = 0
du[2] = 0.5u[2]
end
prob_sde = SDEProblem(f, g, nothing, nothing)
init_f = (col, t) -> [0.0, 0.0]
function derived_f(col, sol, obstimes, subject, param, randeff)
Vc = getfield.(col.(obstimes),:Vc)
Σ = getfield.(col.(obstimes),:Σ)
central = sol(obstimes, idxs = 2)
conc = @. central / Vc
dv = @. Normal(conc, Σ)
(dv = dv, conc = conc, )
end
function observed_f(col, sol, obstimes, samples, subject)
(dv = samples.dv)
end
stoch_model_sde = Pumas.PumasModel(p, rfx, pre_f, init_f, prob_sde, derived_f, observed_f)
```

Then, I simulate data as follows.

```
N = 20
ka = 1*exp.(randn((1, N)))
CL = 10*exp.(randn((1, N)))
V = 40*exp.(randn((1, N)))
k = CL./V
t = 0:0.5:24
Cp = (1000 ./V).*(ka./(ka-k)).*(exp.(-k.*t) - exp.(-ka.*t))
data = DataFrame()
for i = 1:N
tmp = DataFrame(id = i, time = t, amt = 0, evid = 0, cmt = 2, dv = Cp[:, i])
tmp.:dv = convert(Array{Union{Missing, Float64}}, tmp.:dv)
tmp[tmp.:time .== 0, :amt] = 1000
tmp[tmp.:time .== 0, :evid] = 1
tmp[tmp.:time .== 0, :cmt] = 1
tmp[tmp.:time .== 0, :dv] = missing
append!(data, tmp)
end
timeu = u"hr"
concu = u"mg/L"
amtu = u"mg"
data.route = "ev"
ds = read_pumas(data; covariates = Symbol[], observations = Symbol[:dv], id = :id, time = :time, evid = :evid, amt = :amt)
```

Finally, I try to run parameter estimation:

```
param = (θ = [250.0, 3021.0, 0.3], Σ = 0.1, Ω = 1)
f2 = fit(stoch_model_sde, ds, param, Pumas.FOCEI(), alg = SOSRI())
```

Which results in the following warning and, after a few minutes, error:

```
Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
**gradient element 3 of θ is exactly zero. This indicates that θ isn't identified.**
```

Is there any obvious error in this code or any documentation available on parameter estimation with SDEs? E.g. a code example would probably help a lot.

Thanks a lot!