Estimation with stochastic differential equations

estimation
#1

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!

#2

Unfortunately we currently don’t support fitting of models with stochastic ODEs.

#3

Thanks, and do you have any plans to implement this feature in the (nearer) future?

#4

While it’s something that we’d like to support eventually, it’s most likely not something that we’ll implement in the near future since there are many other things we’d like to do as well.

#5

We do support it. Just make sure to utilize an EnsembleProblem and do something similar to the DiffEqFlux docs

#6

I’ll probably get a tutorial on it, but not until September.

#7

Thanks, looking forward to your tutorial!