I’m trying to put together a demonstration of a very simple particle filter for my repository https://github.com/epirecipes/sir-julia. This involves:
- Setting up an array of simulations.
- Simulating for a single timestep.
- Resampling the simulations that match the observed data
- Compute a partial likelihood as the proportion of simulations that are consistent with the data (as likelihoods are either 0 or 1).
- If there are more timepoints, go back to (2).
- Compute the
I need something that works for both in-place and out-of-place versions. My current code is below. A couple of questions:
- How can I reduce allocations/increase speed? Using
reinit!makes the code faster, but then I have to set
false, and then the in-place version doesn’t work.
- Why are the results different for the in-place and out-of-place versions?
julia> @time pfilter(prob_op, p, u0, X, 100000, 1234) 8.030784 seconds (46.90 M allocations: 3.463 GiB, 9.35% gc time) -113.48234265521722 julia> @time pfilter(prob_ip, p, u0, X, 100000, 1234) 8.026255 seconds (39.00 M allocations: 3.210 GiB, 8.44% gc time) -113.42214586062131
# Load packages using OrdinaryDiffEq using Random using Distributions using StatsBase # Set seed Random.seed!(1) # Define a model for a single time-step. This is an out-of-place version. # C is the number of cases per timestep, which is reset for each iteration. function sir_markov(u,p,t) (S,I,_) = u C=0 (β,γ)=p p δt=0.1 nsteps=10 for i in 1:nsteps ifrac = 1-exp(-β*I*δt) rfrac = 1-exp(-γ*δt) infection=rand(Binomial(S,ifrac)) recovery=rand(Binomial(I,rfrac)) S = S-infection I = I+infection-recovery C = C+infection end [S,I,C] end # In place version of the above function sir_markov!(du,u,p,t) (S,I,C) = u C=0 (β,γ)=p δt=0.1 nsteps=10 for i in 1:nsteps ifrac = 1-exp(-β*I*δt) rfrac = 1-exp(-γ*δt) infection=rand(Binomial(S,ifrac)) recovery=rand(Binomial(I,rfrac)) S = S-infection I = I+infection-recovery C = C+infection end @inbounds begin du = S du = I du = C end end ## Time, initial conditions, and parameter values tspan = (0,40) u0 = [990, 10, 0] # S, I, C p = [0.0005, 0.25] # β, γ ## Solve for a single run to generate simulated data prob_ip = DiscreteProblem(sir_markov!, u0, tspan, p, dt=1) prob_op = DiscreteProblem(sir_markov, u0, tspan, p, dt=1) sol_op = solve(prob_op, FunctionMap()) X = hcat(sol_op.u...)[3,2:end] # Simple particle filter function pfilter(prob, p, u0, X, trajectories=1000, seed=1234) prob = remake(prob, p=p, u0=u0) integrators = [init(prob, FunctionMap(),save_everystep=true) for i in 1:trajectories] Random.seed!(seed) liks = zeros(Float64,length(X)) weights = Weights(zeros(Float64,trajectories)) us = [copy(u0) for i in 1:trajectories] idx = collect(1:trajectories) @inbounds for t in 1:length(X) step!.(integrators) x = X[t] [us[i] = integrators[i].u for i in 1:trajectories] [weights[i]=Float64(us[i]==x) for i in 1:trajectories] liks[t] = mean(weights) if mean(weights)==0.0 return -Inf break end sample!(1:trajectories, weights, idx) [reinit!(integrators[i],us[idx[i]]) for i in 1:trajectories] end sum(log.(liks)) end pfilter(prob_op, p, u0, X, 100000, 1234) @time pfilter(prob_op, p, u0, X, 100000, 1234) @time pfilter(prob_ip, p, u0, X, 100000, 1234)