SIR model with exponential waiting times in Julia

[mschauer]

How do you like something very simple, I think you don't have that yet:

https://gist.github.com/mschauer/be1196dedbf8afe7fb8a0b411672b4d6

[slwu89]

hi @mschauer, it looks like there are a few CTMC constructions of the SIR model in this repo, just listed as "jump processes":

1. https://github.com/epirecipes/sir-julia/blob/master/markdown/jump_process/jump_process.md
2. https://github.com/epirecipes/sir-julia/blob/master/markdown/jump_process_gillespie/jump_process_gillespie.md

[mschauer]

Think of it in the same fashion as https://discourse.julialang.org/t/improving-metropolis-hastings-code/67412/7, you could use Turing there too instead.

[sdwfrost]

I think this is already covered for now, unless @mschauer, you think there's something we're missing with the jump process implementations

[mschauer]

unless you think there's something we're missing with the jump process implementations

I would think so:

prob_jump = JumpProblem(prob,Direct(),infection_jump,recovery_jump);
sol_jump = solve(prob_jump,SSAStepper());

has a very different didactical value compared to

 function sir!(us, u, T, p)
    β, c, γ = p.β, p.c, p.γ
    while true
        t, S, I, R = u
       # @show t, u
        N = +(S, I, R)
        τ = randexp()/(S * β*c*I/N), randexp()/(I*γ) # renew both times, because state has changed
        t + min(τ[1], τ[2]) > T && return us
        if τ[1] < τ[2]
            t += τ[1]
            S -= 1
            I += 1
        else
            t += τ[2]
            I -= 1
            R += 1
        end
        u = (t, S, I, R)
        push!(us, u)
    end
end

[sdwfrost]

OK, I can add the plain Julia code, as I agree that there are didactic differences. This no-library approach was in the old repo, and was in an old blog post of mine too

[slwu89]

@sdwfrost if we are going to add plain Julia code, do you think its worth it to do plain Julia for different flavors of how to sample the CTMC? @mschauer has posted a nice Direct method, it might be worth it to show a First Reaction Method and Anderson's Next Reaction Method too.

[mschauer]

I think have something close to Algorithm 2 (the next reaction method) from the Anderson paper, except that I simply resample both times instead what is done in step (8) to reuse elements of the vector of random next event times.

[sdwfrost]

@slwu89 you probably have something along these lines hanging around too?

[slwu89]

Yeah, I'd love to write up an example with a direct method, first reaction method, and Anderson's NRM. I think I also coded Thanh's rejection sampler in R before.

[sdwfrost]

Hi @slwu89 - just seeing if you're still interested in this!