Implement the MechBayes SEIRD ODE and SDE model as a test case

[ChrisRackauckas]

The ODE is needed for a test case too:

a. We are interested in determining how detection rate can affect the accuracy and uncertainty in our forecasts. In particular, suppose we can improve the baseline detection rate by 20%, and the detection rate stays constant throughout the duration of the forecast. Assuming no additional interventions (ignoring Question 3), does that increase the amount of cumulative forecasted cases and deaths after six weeks? How does an increase in the detection rate affect the uncertainty in our estimates? Can you characterize the relationship between detection rate and our forecasts and their uncertainties, and comment on whether improving detection rates would provide decision-makers with better information (i.e., more accurate forecasts and/or narrower prediction intervals)? b. Compute the accuracy of the forecast assuming no mask mandate (ignoring Question 3) in the same way as you did in Question 1 and determine if improving the detection rate improves forecast accuracy.

[ChrisRackauckas]

@variables t
D = Differential(t)
sts = @variables S(t) E(t) I(t) R(t) D1(t) D2(t) β(t) logβ(t)
ps = @parameters σ ρ γ λ N
@brownian α
eqs = Equation[
               β ~ logβ #exp(logβ)
               D(logβ) ~ sqrt(0.2) * α
               D(S) ~ -β * S * I / N
               E ~ N - (S + I + R + D1 + D2)
               D(I) ~ σ * E - γ * I
               D(R) ~ (1 - ρ) * γ * I
               D(D1) ~ ρ * γ * I - λ * D1
               D(D2) ~ λ * D1]
@named mecha_bayes = System(eqs, t)
mecha_bayes = structural_simplify(mecha_bayes)
dE = 0.4
dI = 2.0
Rh = 3.0
NN = 1
using Random
Random.seed!(1)

E0 = rand(Uniform(0, 0.02NN))
I0 = rand(Uniform(0, 0.02NN))
R0 = rand(Uniform(0, 0.02NN))
D10 = rand(Uniform(0, 0.02NN))
D20 = rand(Uniform(0, 0.02NN))
S0 = 1 - (E0 + I0 + R0 + D10 + D20)
u0 = [
      S =>  S0
      E =>  E0
      I =>  I0
      R =>  R0
      D1 => D10
      D2 => D20
      β => 0.0
      logβ => log(rand(Gamma(1, dI / Rh)))
     ]

p = [N => NN,
     σ => rand(Gamma(100, 100dE)),
     ρ => rand(Beta(1, 99)),
     γ => rand(Gamma(100, 100dI)),
     λ => rand(Gamma(10, 10*25))]
prob = SDEProblem(mecha_bayes, u0, (0.0, 10.0), p)
@test_nowarn sol = solve(prob, SRIW1(), reltol=1e-3, abstol=1e-3, seed = 1);
eqs = Equation[
               β ~ 2.6
               D(S) ~ -β * S * I / N
               E ~ N - (S + I + R + D1 + D2)
               D(I) ~ σ * E - γ * I
               D(R) ~ (1 - ρ) * γ * I
               D(D1) ~ ρ * γ * I - λ * D1
               D(D2) ~ λ * D1]
@named mecha_bayes = System(eqs, t)
mecha_bayes = structural_simplify(mecha_bayes)
prob = ODEProblem(mecha_bayes, u0, (0.0, 1e5), p)
sol = solve(prob, Tsit5())