itsdfish / differentialevolutionmcmc.jl Goto Github PK

ERROR: LoadError: UndefVarError: resample not defined

got unsupported keyword arguments "n_initial", "θsnooker"

using Plots, DifferentialEquations 
using DifferentialEvolutionMCMC, Distributions, Random, StatsBase, LabelledArrays

Random.seed!(42)

# ODE model
function SIR(du,u,p,t)

    # states
    S, I, R, C = u
    N = S + I + R
  
    # params
    β = p.β
    η = p.η
    φ = p.φ
    ω = p.ω
    μ = p.μ
    σ = p.σ
      
    # FOI
    βeff = β * (1.0+η*sin(2.0*π*(t-φ)/365.0))
    λ = βeff*I/N

    # change in states
    du[1] = μ*N - λ*S - μ*S + ω*R
    du[2] = λ*S - σ*I - μ*I
    du[3] = σ*I - μ*R - ω*R
    du[4] = σ*I # cumulative incidence
  
end
  
# Calculate incident cases from cumulative incidence
function get_incidence(sol)

    incidence = [sol.u[t].C - sol.u[t-1].C  for t = 2:length(sol.u)]
    return incidence

end

# Solver settings
tmin = 0
tmax = 10*365
tspan = (tmin, tmax)
abstol = 1.0e-12
saveat = 7.0
solver = AutoVern7(KenCarp4())

# Initiate ODE problem
p = @LArray [0.27,0.15,20,1.0/365,1.0/(80*365),1.0/4.98] (:β,:η,:φ,:ω,:μ,:σ)
u0 = @LArray [9999,1,0,1] (:S,:I,:R,:C)
problem = ODEProblem(SIR,u0,tspan,p)
sol = solve(problem, solver, abstol=abstol, saveat=saveat)
plot(sol)
incidence =get_incidence(sol)
plot(incidence)

# Fake some data from model
data = round.(incidence .* rand(Uniform(0.5,1.5), size(incidence)))
plot(incidence, legend = false); scatter!(data,legend = false)


# Fit model to fake data

# Priors
priors = (
    β = (Uniform(0.0,1.0)), 
    η = (Uniform(0.0,1.0)),
    ω = (Uniform(0.0,1.0/(5.0*365.0))),
    φ = (Uniform(0.0,364.0))
)

# Upper and lower bounds (is equivalent to argument of the Uniform)
bounds = (
    (0.0,1.0),
    (0.0,1.0),
    (0.0,1.0/(5.0*365.0)),
    (0.0,364.0)
)

# Log Likelihood function
function loglik(data, θ...)

    β, η, ω, φ = θ
     
    p_new = @LArray [β,η,ω,φ,1.0/(80.0*365.0),1.0/4.98] (:β,:η,:ω,:φ,:μ,:σ)
    # Update problem and solve ODEs
    problem_new = remake(problem, p=p_new)
    sol_new = solve(problem_new,solver, abstol=abstol, saveat=saveat) 
    # calculate incidence
    incsim = get_incidence(sol_new)

    return sum(logpdf.(Poisson.(incsim), data))
end


model = DEModel(; priors, model=loglik, data)

de = DE(; bounds, burnin=1000, priors)
n_iter = 2000
chains = sample(model, de, MCMCThreads(), n_iter, progress=true)
println(chains)

ERROR: LoadError: MethodError: no method matching iterate(::Uniform{Float64})
Closest candidates are:
  iterate(::ExponentialBackOff) at error.jl:252
  iterate(::ExponentialBackOff, ::Any) at error.jl:252
  iterate(::VSCodeServer.JuliaInterpreter.ExprSplitter) at c:\Users\mickey\.vscode\extensions\julialang.language-julia-1.1.38\scripts\packages\JuliaInterpreter\src\construct.jl:554

function loglik(θ, data)
    μ = θ[1]
    σ = θ[2]
    return sum(logpdf.(Normal(μ, σ), data))
end

using DifferentialEquations, StatsBase, StatsPlots, DifferentialEvolutionMCMC, Distributions, Random, StatsBase, LabelledArrays

Random.seed!(42)

# ODE model: simple SIR model with seasonally forced contact rate
function SIR!(du,u,p,t)

    # states
    S, I, R = u
    N = S + I + R
  
    # params
    β = p.β
    η = p.η
    φ = p.φ
    ω = p.ω
    μ = p.μ
    σ = p.σ
      
    # FOI
    βeff = β * (1.0+η*sin(2.0*π*(t-φ)/365.0))
    λ = βeff*I/N

    # change in states
    du[1] = μ*N - λ*S - μ*S + ω*R
    du[2] = λ*S - σ*I - μ*I
    du[3] = σ*I - μ*R - ω*R
    du[4] = σ*I # cumulative incidence
  
end
  
# Calculate incident cases from cumulative incidence
function get_incidence(sol)

    incidence = [sol.u[t].C - sol.u[t-1].C  for t = 2:length(sol.u)]
    return incidence

end

# Solver settings
tmin = 0.0
tmax = 10.0*365.0
tspan = (tmin, tmax)
abstol = 1.0e-8 
reltol = 1.0e-8 
maxiters = 1e7 
saveat = 7.0
solver = AutoVern7(Rodas5())

# Initiate ODE problem
p = (β=0.28,η=0.07,φ=20,ω=1.0/365,μ=1.0/(80*365),σ=1.0/5.0)
u0 = @LArray [9999.0,1.0,0.0,1.0] (:S,:I,:R,:C)

# Initiate ODE problem and run solution with dummy vars
problem = ODEProblem(SIR!,u0,tspan,p)
sol = solve(problem, 
            solver, abstol=abstol, 
            reltol=reltol, 
            maxiters=maxiters,
            isoutofdomain=(u,p,t)->any(x->x<0.0,u), 
            saveat=saveat)
#plot(sol)
incidence = get_incidence(sol)
plot(incidence)

# Fake some data from model
data = round.(incidence .* rand(Uniform(0.5,1.5), size(incidence)))
#plot(incidence, legend = false); scatter!(data,legend = false)
#data = @. round(incidence * rand(Uniform(0.5,1.5), size(incidence)))

# Fit model to fake data

# Priors
priors = (
    β = (Uniform(1e-3,1.0),), 
    η = (Uniform(0.0,1.0),),
    ω = (Uniform(1.0/(5.0*365.0),1),),
    φ = (Uniform(0.0,364.0),)
)

# Upper and lower bounds (in this case equivalent to arguments of the Uniform)
bounds = (
    (1e-3,1.0),
    (0.0,1.0),
    (1.0/(5.0*365.0),1.0),
    (0.0,364.0)
)

# Log Likelihood function
function loglik(data, problem, theta_fix, θ...)

    p_new = merge(theta_fix, θ) # does not work, because θ is not a named tuple
    
    # Update problem and solve ODEs
    problem_new = remake(problem, p=p_new)
    sol_new = solve(problem_new, 
                    solver, abstol=abstol, 
                    reltol=reltol, 
                    maxiters=maxiters,
                    isoutofdomain=(u,p,t)->any(x->x<0.0,u), 
                    saveat=saveat) 
    # calculate incidence
    incidence = get_incidence(sol_new)

    # due to numerical instability of the diff eqs, some states can become negative. In this case, return -Inf <-- this is obsolete actually, because of the isoutofdomain argument. 
    #if minimum(incidence) < 0.0
        #ll = -Inf
   # else
        ll = sum(logpdf.(Poisson.(incidence), data))
    #end
    return ll
end


# Run DE MCMC
burnin=0
n_iter = 5000
n_groups = 3 # default is 4
Np = length(priors)*3 # N of particles per group.Default is nparam*3
model = DEModel(problem, theta_fix; priors, model=loglik, data=data)
de = DE(; n_groups, Np, bounds, burnin=burnin, priors)
@time chains = sample(model, de, MCMCThreads(), n_iter, progress=true, discard_burnin=false)