All infected individuals go through an incubation period, which lasts on average for 4 days. During this time, individuals are not infectious, nor do they have any excess mortality risk.
All infected individuals eventually develop symptoms, and the mean duration of symptoms before either recovery or death is 5 days. Symptomatic individuals are infectious, as well as having a 3% case fatality rate. Those who survive the infection are thought to have long-term immunity.
In the source country, the peak prevalence (i.e. maximum number of symptomatic people during the epidemic) was observed to be 8% of the population.
The total population of this country is N=11.000.000
\begin{align*} \text{Case Fatality Rate (CFR)} &= 3\% = 0.03 \\ \text{Mean duration of symptoms (}D\text{)} &= 5 \text{ days} \end{align*}
The rate of recovery ($\gamma$) can be calculated using the formula:
$$ \gamma = \frac{1 - \text{CFR}}{D} $$
$$\gamma = \frac{0.97}{5} = 0.194 \text{ per day}$$
Then ($\mu$) is $$CFR = \frac{\mu}{\mu + \gamma}$$ So $$\mu=0.006$$ The estimated value of $(I_{peak})$ can be calculated as: $$I_{\text{peak}} = \text{Peak Prevalence} \times N = 0.08 \times 11,000,000 = 880,000$$
The differential equations for this model is \begin{align*} \frac{dS}{dt}&= -\beta \cdot S(t)\cdot \frac{I(t)}{N} - \mu \cdot S(t)\\ \frac{dI}{dt} &= \beta \cdot S(t) \cdot \frac{I(t)}{N} - \gamma \cdot I(t) - \mu \cdot I(t)\\ \frac{dR}{dt} &= \gamma \cdot I(t) - \mu \cdot R(t) \end{align*} What estimate do get for $\beta$ which corresponds to a peak prevalence (of symptomatic infection) of 0.08(i.e. 8%)
(*total size of population*)n = 11000000;
times = 100;
reportedData = n*0.08
sireqns = {s'[t] == -\[Beta]*s[t]*i[t]/n - \[Mu]*s[t],
i'[t] == \[Beta]*s[t]*i[t]/n - \[Gamma]*i[t] - \[Mu]*i[t],
r'[t] == \[Gamma]*i[t] - \[Mu]*r[t]};
initialConditions = {s[0] == n, i[0] == 1, r[0] == 0};
sol = ParametricNDSolve[{sireqns, initialConditions}, {s, i, r}, {t,
0, times}, {\[Beta], \[Gamma], \[Mu]}]
(*Extract the infected values from the solution*)
infectedValues = i /. sol;
I have tried the following code to Define a function to calculate the Poisson log-likelihood but it does not work.
loglik[\[Beta]_, \[Gamma]_, \[Mu]_, data_] :=
Module[{output, modelI, \[Lambda], loglik},
output = {i} /. sol[\[Beta], \[Gamma], \[Mu]];
modelI =
Interpolation[Transpose[{Range[0, times], output[[1]]}],
InterpolationOrder -> 1];
\[Lambda] = 0.6*modelI[data[[All, 1]]];
loglik =
Total[Log[PoissonDistribution[\[Lambda]] /@ data[[All, 2]]]];
loglik]
result =
NMaximize[{loglik[\[Beta], \[Gamma], \[Mu], reportedData],
0 <= \[Beta] <= 1, 0 <= \[Gamma] <= 1,
0 <= \[Mu] <= 1}, {\[Beta], \[Gamma], \[Mu]}]
estimatedParameters = {\[Beta], \[Gamma], \[Mu]} /. result[[2]]
I want to define a function that simulates the model for a given combination of parameters and calculates the Poisson log-likelihood for the epidemic curve of reported cases as R code shows us:
# DISTANCE FUNCTION
loglik_function <- function(parameters, dat) { # takes as inputs the parameter values and dataset
beta <- parameters[1] # extract and save the first value in the "parameters" input argument as beta
gamma <- parameters[2] # extract and save the second value in the "parameters" input argument as gamma
mu <- parameters[2] # extract and save the second value in the "parameters" input argument as mu
# Simulate the model with initial conditions and timesteps defined above, and parameter values from function call
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = c(beta = beta, # ode() takes the values for beta and gamma extracted from
gamma = gamma,
mu=mu))) # the "parameters" input argument of the loglik_function()
So far I have figured out how to code it. My question is how do we generate this code in Mathematica, do I need Poisson Distribution? Could some one explain to me how could i code that?
LL <- sum(dpois(x = dat$number_reported, lambda = 0.6 * output$I[output$time %in% dat$time], log = TRUE))
return(LL)
}
# OPTIMISATION
optim(par = c(1.7, 0.1), # starting values for beta and gamma - you should get the same result no matter
# which values you choose here
fn = loglik_function, # the distance function to optimise
dat = reported_data, # the dataset to fit to ("dat" argument is passed to the function specified in fn)
control = list(fnscale=-1)) # tells optim() to look for the maximum number instead of the minimum (the default)
Also, I know that the result of $\beta=0.38$
If someone needs it this the whole code I used in R
# Load required packages
library(deSolve)
# Define initial conditions
initial_state_values <- c(S = 11000000,
I = 1,
R = 0)
# Define time points
times <- seq(from = 0, to = 100, by = 0.1)
# Define SIR model function
sir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + R
lambda <- beta * I / N
dS <- -lambda * S - mu * S
dI <- lambda * S - gamma * I - mu * I
dR <- gamma * I - mu * R
return(list(c(dS, dI, dR)))
})
}
# Known values
gamma <- 0.194
mu <- 0.006
Ipeak <- 880000
# Define a range of possible beta values
possible_beta_values <- seq(0, 1, by = 0.01)
# Find beta that results in Ipeak
for (beta in possible_beta_values) {
parameters <- c(beta = beta, gamma = gamma, mu = mu)
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters))
peak_time <- times[which.max(output$I)]
if (output$I[peak_time / 0.1 + 1] >= Ipeak) {
cat("Found beta:", beta, "\n")
break
}
}
datais not defined. Is it the same as in your previous question? $\endgroup$