Hi everyone,
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
set_u!instead ofreinit!makes the code faster, but then I have to setsave_everysteptofalse, 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[1] = S
du[2] = I
du[3] = 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][3]==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)