From 0dbf48d8faca768798e7d8604700ecc6d1fb42df Mon Sep 17 00:00:00 2001 From: annabellasd Date: Mon, 25 Sep 2023 20:20:43 -0800 Subject: [PATCH] updated lectures --- lectures/continuous_time/covid_sde.md | 121 +++++++++++++------------ lectures/continuous_time/seir_model.md | 81 +++++++++-------- 2 files changed, 103 insertions(+), 99 deletions(-) diff --git a/lectures/continuous_time/covid_sde.md b/lectures/continuous_time/covid_sde.md index c0576465..b13b6b85 100644 --- a/lectures/continuous_time/covid_sde.md +++ b/lectures/continuous_time/covid_sde.md @@ -230,40 +230,41 @@ First, construct our $F$ from {eq}`dfcvsde` and $G$ from {eq}`dG` ```{code-cell} julia function F(x, p, t) - s, i, r, d, R₀, δ = x - (;γ, R̄₀, η, σ, ξ, θ, δ_bar) = p - - return [-γ*R₀*s*i; # ds/dt - γ*R₀*s*i - γ*i; # di/dt - (1-δ)*γ*i; # dr/dt - δ*γ*i; # dd/dt - η*(R̄₀(t, p) - R₀);# dR₀/dt - θ*(δ_bar - δ); # dδ/dt - ] + s, i, r, d, R_0, delta = x + (; gamma, R_bar_0, eta, sigma, xi, theta, delta_bar) = p + + return [-gamma * R_0 * s * i; # ds/dt + gamma * R_0 * s * i - gamma * i; # di/dt + (1 - delta) * gamma * i; # dr/dt + delta * gamma * i; # dd/dt + eta * (R_bar_0(t, p) - R_0); # dR_0/dt + theta * (delta_bar - delta); # ddelta/dt + ] end function G(x, p, t) - s, i, r, d, R₀, δ = x - (;γ, R̄₀, η, σ, ξ, θ, δ_bar) = p + s, i, r, d, R_0, delta = x + (; gamma, R_bar_0, eta, sigma, xi, theta, delta_bar) = p - return [0; 0; 0; 0; σ*sqrt(R₀); ξ*sqrt(δ * (1-δ))] + return [0; 0; 0; 0; sigma * sqrt(R_0); xi * sqrt(delta * (1 - delta))] end ``` Next create a settings generator, and then define a [SDEProblem](https://docs.sciml.ai/stable/tutorials/sde_example/#Example-2:-Systems-of-SDEs-with-Diagonal-Noise-1) with Diagonal Noise. ```{code-cell} julia -p_gen = @with_kw (T = 550.0, γ = 1.0 / 18, η = 1.0 / 20, - R₀_n = 1.6, R̄₀ = (t, p) -> p.R₀_n, δ_bar = 0.01, - σ = 0.03, ξ = 0.004, θ = 0.2, N = 3.3E8) -p = p_gen() # use all defaults +p_gen(; T = 550.0, gamma = 1.0 / 18, eta = 1.0 / 20, + R_0_n = 1.6, R_bar_0 = (t, p) -> p.R_0_n, delta_bar = 0.01, + sigma = 0.03, xi = 0.004, theta = 0.2, N = 3.3E8) = (; T, gamma, eta, R_0_n, R_bar_0, delta_bar, sigma, xi, theta, N) + +p = p_gen() # use all defaults i_0 = 25000 / p.N r_0 = 0.0 d_0 = 0.0 s_0 = 1.0 - i_0 - r_0 - d_0 -R̄₀_0 = 0.5 # starting in lockdown -δ_0 = p.δ_bar -x_0 = [s_0, i_0, r_0, d_0, R̄₀_0, δ_0] +R_bar_0_0 = 0.5 # starting in lockdown +delta_0 = p.delta_bar +x_0 = [s_0, i_0, r_0, d_0, R_bar_0_0, delta_0] prob = SDEProblem(F, G, x_0, (0, p.T), p) ``` @@ -370,29 +371,29 @@ We will shut down the shocks to the mortality rate (i.e. $\xi = 0$) to focus on Consider $\eta = 1/50$ and $\eta = 1/20$, where we start at the same initial condition of $R_0(0) = 0.5$. ```{code-cell} julia -function generate_η_experiment(η; p_gen = p_gen, trajectories = 100, - saveat = 1.0, x_0 = x_0, T = 120.0) - p = p_gen(η = η, ξ = 0.0) +function generate_eta_experiment(eta; p_gen, trajectories = 100, saveat = 1.0, x_0, T = 120.0) + p = p_gen(eta = eta, xi = 0.0) ensembleprob = EnsembleProblem(SDEProblem(F, G, x_0, (0, T), p)) - sol = solve(ensembleprob, SOSRI(), EnsembleThreads(), - trajectories = trajectories, saveat = saveat) + sol = solve(ensembleprob, SOSRI(), EnsembleThreads(), trajectories = trajectories, saveat = saveat) return EnsembleSummary(sol) end # Evaluate two different lockdown scenarios -η_1 = 1/50 -η_2 = 1/20 -summ_1 = generate_η_experiment(η_1) -summ_2 = generate_η_experiment(η_2) -plot(summ_1, idxs = (4,5), +eta_1 = 1 / 50 +eta_2 = 1 / 20 +summ_1 = generate_eta_experiment(eta_1) +summ_2 = generate_eta_experiment(eta_2) + +plot(summ_1, idxs = (4, 5), title = ["Proportion Dead" L"R_0"], ylabel = [L"d(t)" L"R_0(t)"], xlabel = L"t", legend = [:topleft :bottomright], - labels = L"Middle 95% Quantile, $\eta = %$η_1$", + labels = L"Middle 95% Quantile, $\eta = %$eta_1$", layout = (2, 1), size = (900, 900), fillalpha = 0.5) -plot!(summ_2, idxs = (4,5), + +plot!(summ_2, idxs = (4, 5), legend = [:topleft :bottomright], - labels = L"Middle 95% Quantile, $\eta = %$η_2$", size = (900, 900), fillalpha = 0.5) + labels = L"Middle 95% Quantile, $\eta = %$eta_2$", size = (900, 900), fillalpha = 0.5) ``` While the the mean of the $d(t)$ increases, unsurprisingly, we see that the 95% quantile for later time periods is also much larger - even after the $R_0$ has converged. @@ -411,24 +412,24 @@ Since empirical estimates of $R_0(t)$ discussed in {cite}`NBERw27128` and other We start the model with 100,000 active infections. ```{code-cell} julia -R₀_L = 0.5 # lockdown -η_experiment = 1.0/10 -σ_experiment = 0.04 -R̄₀_lift_early(t, p) = t < 30.0 ? R₀_L : 2.0 -R̄₀_lift_late(t, p) = t < 120.0 ? R₀_L : 2.0 +R_0_L = 0.5 # lockdown +eta_experiment = 1.0 / 10 +sigma_experiment = 0.04 -p_early = p_gen(R̄₀ = R̄₀_lift_early, η = η_experiment, σ = σ_experiment) -p_late = p_gen(R̄₀ = R̄₀_lift_late, η = η_experiment, σ = σ_experiment) +R_bar_0_lift_early(t, p) = t < 30.0 ? R_0_L : 2.0 +R_bar_0_lift_late(t, p) = t < 120.0 ? R_0_L : 2.0 +p_early = p_gen(R_bar_0 = R_bar_0_lift_early, eta = eta_experiment, sigma = sigma_experiment) +p_late = p_gen(R_bar_0 = R_bar_0_lift_late, eta = eta_experiment, sigma = sigma_experiment) # initial conditions i_0 = 100000 / p_early.N r_0 = 0.0 d_0 = 0.0 s_0 = 1.0 - i_0 - r_0 - d_0 -δ_0 = p_early.δ_bar +delta_0 = p_early.delta_bar -x_0 = [s_0, i_0, r_0, d_0, R₀_L, δ_0] # start in lockdown +x_0 = [s_0, i_0, r_0, d_0, R_0_L, delta_0] # start in lockdown prob_early = SDEProblem(F, G, x_0, (0, p_early.T), p_early) prob_late = SDEProblem(F, G, x_0, (0, p_late.T), p_late) ``` @@ -438,13 +439,13 @@ Simulating for a single realization of the shocks, we see the results are qualit ```{code-cell} julia sol_early = solve(prob_early, SOSRI()) sol_late = solve(prob_late, SOSRI()) -plot(sol_early, vars = [5, 1,2,4], +plot(sol_early, vars = [5, 1, 2, 4], title = [L"R_0" "Susceptible" "Infected" "Dead"], layout = (2, 2), size = (900, 600), ylabel = [L"R_0(t)" L"s(t)" L"i(t)" L"d(t)"], xlabel = L"t", legend = [:bottomright :topright :topright :topleft], label = ["Early" "Early" "Early" "Early"]) -plot!(sol_late, vars = [5, 1,2,4], +plot!(sol_late, vars = [5, 1, 2, 4], legend = [:bottomright :topright :topright :topleft], label = ["Late" "Late" "Late" "Late"], xlabel = L"t") ``` @@ -535,24 +536,26 @@ We will redo the "Ending Lockdown" simulation from above, where the only differe ```{code-cell} julia function F_reinfect(x, p, t) - s, i, r, d, R₀, δ = x - (;γ, R̄₀, η, σ, ξ, θ, δ_bar, ν) = p - - return [-γ*R₀*s*i + ν*r; # ds/dt - γ*R₀*s*i - γ*i; # di/dt - (1-δ)*γ*i - ν*r # dr/dt - δ*γ*i; # dd/dt - η*(R̄₀(t, p) - R₀);# dR₀/dt - θ*(δ_bar - δ); # dδ/dt - ] + s, i, r, d, R_0, delta = x + (;gamma, R_bar_0, eta, sigma, xi, theta, delta_bar, nu) = p + + return [ + -gamma * R_0 * s * i + nu * r; # ds/dt + gamma * R_0 * s * i - gamma * i; # di/dt + (1 - delta) * gamma * i - nu * r; # dr/dt + delta * gamma * i; # dd/dt + eta * (R_bar_0(t, p) - R_0); # dR_0/dt + theta * (delta_bar - delta); # ddelta/dt + ] end -p_re_gen = @with_kw ( T = 550.0, γ = 1.0 / 18, η = 1.0 / 20, - R₀_n = 1.6, R̄₀ = (t, p) -> p.R₀_n, - δ_bar = 0.01, σ = 0.03, ξ = 0.004, θ = 0.2, N = 3.3E8, ν = 1/360) +p_re_gen(; T = 550.0, gamma = 1.0 / 18, eta = 1.0 / 20, + R_0_n = 1.6, R_bar_0 = (t, p) -> p.R_0_n, + delta_bar = 0.01, sigma = 0.03, xi = 0.004, theta = 0.2, + N = 3.3E8, nu = 1/360) = (; T, gamma, eta, R_0_n, R_bar_0, delta_bar, sigma, xi, theta, N, nu) -p_re_early = p_re_gen(R̄₀ = R̄₀_lift_early, η = η_experiment, σ = σ_experiment) -p_re_late = p_re_gen(R̄₀ = R̄₀_lift_late, η = η_experiment, σ = σ_experiment) +p_re_early = p_re_gen(R_bar_0 = R_bar_0 _lift_early, eta = eta_experiment, sigma = sigma_experiment) +p_re_late = p_re_gen(R_bar_0 = R_bar_0 _lift_late, eta = eta_experiment, sigma = sigma_experiment) trajectories = 400 saveat = 1.0 diff --git a/lectures/continuous_time/seir_model.md b/lectures/continuous_time/seir_model.md index c36fdb06..14c98d0c 100644 --- a/lectures/continuous_time/seir_model.md +++ b/lectures/continuous_time/seir_model.md @@ -171,13 +171,13 @@ We begin by implementing a simple version of this model with a constant $R_0$ an First, define the system of equations ```{code-cell} julia -function F_simple(x, p, t; γ = 1/18, R₀ = 3.0, σ = 1/5.2) +function F_simple(x, p, t; gamma = 1/18, R_0 = 3.0, sigma = 1/5.2) s, e, i, r = x - return [-γ*R₀*s*i; # ds/dt = -γR₀si - γ*R₀*s*i - σ*e;# de/dt = γR₀si -σe - σ*e - γ*i; # di/dt = σe -γi - γ*i; # dr/dt = γi + return [-gamma * R_0 * s * i; # ds/dt + gamma * R_0 * s * i - sigma * e; # de/dt + sigma * e - gamma * i; # di/dt + gamma * i; # dr/dt ] end ``` @@ -195,7 +195,7 @@ tspan = (0.0, 350.0) # ≈ 350 days prob = ODEProblem(F_simple, x_0, tspan) ``` -With this, choose an ODE algorithm and solve the initial value problem. A good default algorithm for non-stiff ODEs of this sort might be `Tsit5()`, which is the Tsitouras 5/4 Runge-Kutta method). +With this, choose an ODE algorithm and solve the initial value problem. A good default algorithm for non-stiff ODEs of this sort might be `Tsit5()`, which is the Tsitouras 5/4 Runge-Kutta method. ```{code-cell} julia sol = solve(prob, Tsit5()) @@ -242,7 +242,7 @@ While we could integrate the deaths given the solution to the model ex-post, it This is a common trick when solving systems of ODEs. While equivalent in principle to using the appropriate quadrature scheme, this becomes especially convenient when adaptive time-stepping algorithms are used to solve the ODEs (i.e. there is not a regular time grid). Note that when doing so, $d(0) = \int_0^0 \delta \gamma i(\tau) d \tau = 0$ is the initial condition. -The system {eq}`seir_system` and the supplemental equations can be written in vector form $x := [s, e, i, r, R₀, c, d]$ with parameter tuple $p := (\sigma, \gamma, \eta, \delta, \bar{R}_0(\cdot))$ +The system {eq}`seir_system` and the supplemental equations can be written in vector form $x := [s, e, i, r, R_0, c, d]$ with parameter tuple $p := (\sigma, \gamma, \eta, \delta, \bar{R}_0(\cdot))$ Note that in those parameters, the targeted reproduction number, $\bar{R}_0(t)$, is an exogenous function. @@ -289,35 +289,36 @@ First, construct our $F$ from {eq}`dfcv` ```{code-cell} julia function F(x, p, t) - s, e, i, r, R₀, c, d = x - (;σ, γ, R̄₀, η, δ) = p - - return [-γ*R₀*s*i; # ds/dt - γ*R₀*s*i - σ*e; # de/dt - σ*e - γ*i; # di/dt - γ*i; # dr/dt - η*(R̄₀(t, p) - R₀);# dR₀/dt - σ*e; # dc/dt - δ*γ*i; # dd/dt + s, e, i, r, R_0, c, d = x + (;sigma, gamma, R_0_bar, eta, delta) = p + + return [-gamma * R_0 * s * i; # ds/dt + gamma * R_0 * s * i - sigma * e; # de/dt + sigma * e - γ*i; # di/dt + gamma * i; # dr/dt + eta * (R_bar_0(t, p) - R_0); # dR_0/dt + sigma * e; # dc/dt + delta * gamma * i; # dd/dt ] end; + ``` This function takes the vector `x` of states in the system and extracts the fixed parameters passed into the `p` object. -The only confusing part of the notation is the `R̄₀(t, p)` which evaluates the `p.R̄₀` at this time (and also allows it to depend on the `p` parameter). +The only confusing part of the notation is the `R_0(t, p)` which evaluates the `p.R_bar_0` at this time (and also allows it to depend on the `p` parameter). ### Parameters - +#####HERE ENDED The baseline parameters are put into a named tuple generator (see previous lectures using [Parameters.jl](https://github.com/mauro3/Parameters.jl)) with default values discussed above. ```{code-cell} julia -p_gen = @with_kw ( T = 550.0, γ = 1.0 / 18, σ = 1 / 5.2, η = 1.0 / 20, - R₀_n = 1.6, δ = 0.01, N = 3.3E8, - R̄₀ = (t, p) -> p.R₀_n); +p_gen(T = 550.0, gamma = 1.0 / 18, sigma = 1 / 5.2, eta = 1.0 / 20, + R_0_n = 1.6, delta = 0.01, N = 3.3E8, + R_0_bar = (t, p) -> p.R_0_n) = (T, gamma, sigma, eta, R_0_n, delta, N, R_0_bar) ``` -Note that the default $\bar{R}_0(t)$ function always equals $R_{0n}$ -- a parameterizable natural level of $R_0$ used only by the `R̄₀` function +Note that the default $\bar{R}_0(t)$ function always equals $R_{0n}$ -- a parameterizable natural level of $R_0$ used only by the `R_hat_0` function Setting initial conditions, we choose a fixed $s, i, e, r$, as well as $R_0(0) = R_{0n}$ and $m(0) = 0.01$ @@ -328,7 +329,7 @@ i_0 = 1E-7 e_0 = 4.0 * i_0 s_0 = 1.0 - i_0 - e_0 -x_0 = [s_0, e_0, i_0, 0.0, p.R₀_n, 0.0, 0.0] +x_0 = [s_0, e_0, i_0, 0.0, p.R_0_n, 0.0, 0.0] tspan = (0.0, p.T) prob = ODEProblem(F, x_0, tspan, p) ``` @@ -372,9 +373,9 @@ Let's start with the case where $\bar{R}_0(t) = R_{0n}$ is constant. We calculate the time path of infected people under different assumptions of $R_{0n}$: ```{code-cell} julia -R₀_n_vals = range(1.6, 3.0, length = 6) -sols = [solve(ODEProblem(F, x_0, tspan, p_gen(R₀_n = R₀_n)), - Tsit5(), saveat=0.5) for R₀_n in R₀_n_vals]; +R_0_n_vals = range(1.6, 3.0, length = 6) +sols = [solve(ODEProblem(F, x_0, tspan, p_gen(R_0_n = R_0_n)), + Tsit5(), saveat=0.5) for R_0_n in R_0_n_vals]; ``` Here we chose `saveat=0.5` to get solutions that were evenly spaced every `0.5`. @@ -384,7 +385,7 @@ Changing the saved points is just a question of storage/interpolation, and does Let's plot current cases as a fraction of the population. ```{code-cell} julia -labels = permutedims([L"R_0 = %$r" for r in R₀_n_vals]) +labels = permutedims([L"R_0 = %$r" for r in R_0_n_vals]) infecteds = [sol[3,:] for sol in sols] plot(infecteds, label=labels, legend=:topleft, lw = 2, xlabel = L"t", ylabel = L"i(t)", title = "Current Cases") @@ -413,21 +414,21 @@ In the simple case, where $\bar{R}_0(t) = R_{0n}$ is independent of the state, t We will examine the case where $R_0(0) = 3$ and then it falls to $R_{0n} = 1.6$ due to the progressive adoption of stricter mitigation measures. -The parameter `η` controls the rate, or the speed at which restrictions are +The parameter `eta` controls the rate, or the speed at which restrictions are imposed. We consider several different rates: ```{code-cell} julia -η_vals = [1/5, 1/10, 1/20, 1/50, 1/100] -labels = permutedims([L"\eta = %$η" for η in η_vals]); +eta_vals = [1/5, 1/10, 1/20, 1/50, 1/100] +labels = permutedims([L"\eta = %$eta" for eta in eta_vals]); ``` Let's calculate the time path of infected people, current cases, and mortality ```{code-cell} julia x_0 = [s_0, e_0, i_0, 0.0, 3.0, 0.0, 0.0] -sols = [solve(ODEProblem(F, x_0, tspan, p_gen(η=η)), Tsit5(), saveat=0.5) for η in η_vals]; +sols = [solve(ODEProblem(F, x_0, tspan, p_gen(eta=eta)), Tsit5(), saveat=0.5) for eta in eta_vals]; ``` Next, plot the $R_0$ over time: @@ -468,11 +469,11 @@ The parameters considered here start the model with 25,000 active infections and 75,000 agents already exposed to the virus and thus soon to be contagious. ```{code-cell} julia -R₀_L = 0.5 # lockdown -R̄₀_lift_early(t, p) = t < 30.0 ? R₀_L : 2.0 -R̄₀_lift_late(t, p) = t < 120.0 ? R₀_L : 2.0 -p_early = p_gen(R̄₀= R̄₀_lift_early, η = 10.0) -p_late = p_gen(R̄₀= R̄₀_lift_late, η = 10.0) +R_0_L = 0.5 # lockdown +R_bar_0_lift_early(t, p) = t < 30.0 ? R_0_L : 2.0 +R_bar_0_lift_late(t, p) = t < 120.0 ? R_0_L : 2.0 +p_early = p_gen(R_bar_0= R_bar_0_lift_early, eta = 10.0) +p_late = p_gen(R_bar_0= R_bar_0_lift_late, eta = 10.0) # initial conditions @@ -480,9 +481,9 @@ i_0 = 25000 / p_early.N e_0 = 75000 / p_early.N s_0 = 1.0 - i_0 - e_0 -x_0 = [s_0, e_0, i_0, 0.0, R₀_L, 0.0, 0.0] # start in lockdown +x_0 = [s_0, e_0, i_0, 0.0, R_0_L, 0.0, 0.0] # start in lockdown -# create two problems, with rapid movement of R₀(t) towards R̄₀(t) +# create two problems, with rapid movement of R_0(t) towards R_hat_0(t) prob_early = ODEProblem(F, x_0, tspan, p_early) prob_late = ODEProblem(F, x_0, tspan, p_late) ``` @@ -501,7 +502,7 @@ plot!(sol_late, vars = [7], label = "Lift Late", xlabel = L"t") Next we examine the daily deaths, $\frac{d D(t)}{dt} = N \delta \gamma i(t)$. ```{code-cell} julia -flow_deaths(sol, p) = p.N * p.δ * p.γ * sol[3,:] +flow_deaths(sol, p) = p.N * p.delta * p.gamma * sol[3,:] plot(sol_early.t, flow_deaths(sol_early, p_early), title = "Flow Deaths", label = "Lift Early") plot!(sol_late.t, flow_deaths(sol_late, p_late), label = "Lift Late", xlabel = L"t")