itsdfish / differentialevolutionmcmc.jl Goto Github PK
View Code? Open in Web Editor NEWA Julia package for Differential Evolution MCMC
License: MIT License
A Julia package for Differential Evolution MCMC
License: MIT License
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.4.5
@JuliaRegistrator register()
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.4.0
@JuliaRegistrator register()
Hi @itsdfish
I have just updated to v0.5.3 and tried to run your example (https://github.com/itsdfish/DifferentialEvolutionMCMC.jl/blob/master/Examples/Multivariate_Guassian_Example.jl) for the snooker update, and it throws an error:
ERROR: LoadError: UndefVarError: resample not defined
When I omit the sample=resample argument, it throws another error:
got unsupported keyword arguments "n_initial", "θsnooker"
Thanks
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.3.2
@JuliaRegistrator register()
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.5.4
Hey
I can replicate your example without any errors. But when I try to fit an ODE model, it throws an error I cannot explain. The model specs should be fine, I can't make any sense of this. Thanks in advance! Bw, Fabienne
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
Hi
I would like to use DifferentialEvolutionMCMC for fitting a complex ODE model (with DifferentialEquations). I like that the priors and the log-likelihood function are defined separately (opposed to for example Turing). However, if there are multiple parameters to be estimated, it seems a bit annoying to list each of them as a separate argument in the loglik function. Is there a way to have a theta vector as one argument, which is matched to the priors in the DE model? The following example (based on your example) does not work:
function loglik(θ, data)
μ = θ[1]
σ = θ[2]
return sum(logpdf.(Normal(μ, σ), data))
end
returns
ERROR: LoadError: MethodError: no method matching loglik(::Float64, ::Float64, ::Array{Float64,1})
I assume the DE function expects n+1 arguments from the loglik function for a prior function with n parameters. Is there a way to adapt this?
Thanks, Fabienne
@JuliaRegistrator register()
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.3.2
@JuliaRegistrator register()
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.4.1
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.4.2
@JuliaRegistrator register()
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.3.0
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.3.1
Hi
I have a follow-up question for my ODE model. I want to estimate only parts of the parameters and fix the rest. I have adjusted my loglik (model) function accordingly. Unfortunately, the DE sampler does not return a Named Tuple for the parameter draws, which is why my attempt below does not work. I was hoping for this solution to work, so that I don't have to write out all the estimated parameters names. That's partly laziness, but also convenience if I want to change which params are estimated. Is there any way I can retrieve the parameter names? Thanks.
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)
Hi @itsdfish
I have been toying with the DEzs MCMC now for a while. Initially I had poor mixing for all model parameters, but I tuned the sampler parameters a bit and achieved acceptable (albeit not perfect) chain mixing for all of them except one. The results below are achieved with burnin=250000, n_iter=300000, n_groups=3, Np=3, θsnooker=0.1, burnin, α=0.2, β=0.2, ϵ=0.001, σ=0.1, κ=0.5
. They are also not fully converged, so may need to run longer. I will continue to improve the mixing also for ρ₀, but do you have any tips in the meanwhile? Could it be that the mixing of ρ₀ is so poor because its scale is much smaller than the other params? TIA
Hei
I am running your example with 10,000 iterations and no burn-in. I fit it repeatedly and noticed that occasionally, some chains converge only after several thousands of iterations. I have added one example. The burn-in of 1000 and iterations of 2000 in your example would thus be way to few. Is this something to be worried about? Thanks, Fabienne
Hi
Is the calculation of rhat and ESS based on the chains with or without the burnin? Does it matter for the calculation whether I specify discard_burnin=false
or discard_burnin=true
? Thanks
@JuliaRegistrator register()
Hey @itsdfish
I have a question regarding the arguments of the DE object. What does the n_groups argument stand for? I understand that the standard recommendation is to have 2-3 times as many chains as parameters, which is set with Np. What does the n_groups argument do, is it related to block updating? The default is 4, which leads to 12 times as many chains as parameters. If, for the sake of speed, I wanted to reduce the number of chains, should I reduce Np or n_groups? TIA
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.5.0
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.2.1
Hi @itsdfish
I was wondering what the official citation for your library is (for scientific publications)? Does it have DOI?
@JuliaRegistrator register()
Hi,
I was trying to fit a model with DEzs (θsnooker=0.1, generate_proposal=random_gamma
) and 3 independent groups of chains (n_groups = 3
and Np=3
). I set α=0.0
so that there would be no exchange between the groups, and I expected the three groups to be independent samples. However, the trace looks like they are not truly independent. How is this possible? This was fit with nthreads=3
. Thanks for your insight.
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.4.4
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.3.3
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.4.3
@JuliaRegistrator register()
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.5.2
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.2.0
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.3.4
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.5.3
@JuliaRegistrator register()
This issue is used to trigger TagBot; feel free to unsubscribe.
If you haven't already, you should update your TagBot.yml
to include issue comment triggers.
Please see this post on Discourse for instructions and more details.
If you'd like for me to do this for you, comment TagBot fix
on this issue.
I'll open a PR within a few hours, please be patient!
@JuliaRegistrator register() DifferentialEvolutionMCMC v0.5.1
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.