# SIR Parameter Estimation in Mathematica

• 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
}
}

• "... but it does not work" - how does it fail exactly? Where is the failure? Commented Sep 2, 2023 at 1:54
• Parameter data is not defined. Is it the same as in your previous question? Commented Sep 2, 2023 at 11:44
• You need something like Total[Map[Log[PDF[PoissonDistribution[\Lambda], #]&, data[[All,2]] ].]. This is equivalent to using the density dpois. The PDF is missing in what you did.
– Asim
Commented Sep 2, 2023 at 13:10
• @MarcoB I get the following message when I run the code, you will see it in my updated post Commented Sep 2, 2023 at 14:14

It looks like we try to estimate parameter $$\beta$$ with known $$\gamma, \mu$$ using one point only reportedData = n*0.08. In this case we can compute $$\beta$$ as follows

n = 11000000;
times = 100;
reportedData = n*0.08;
model[\[Beta]_?NumberQ] := (model[\[Beta]] =
First[i /.
NDSolve[{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], s[0] == n, i[0] == 1,
r[0] == 0} /. {\[Gamma] -> 0.194, \[Mu] -> 0.006},
i, {t, times}]]); sol =
NMinimize[{(Max[model[b][#] & /@ Range[0, times]] -
reportedData)^2, {b > 0}}, {b}, MaxIterations -> 1000]

(*Out[]= {1.95698*10^-17, {b -> 0.510274}}*)


Visualization

Plot[model[b][t] /. sol[[2]], {t, 0, 100}]


We also can try to estimate $$\beta, \gamma, \mu$$ using one point only

model1[\[Beta]_?NumberQ, \[Gamma]_?NumberQ, \[Mu]_?
NumberQ] := (model1[\[Beta], \[Gamma], \[Mu]] =
First[i /.
NDSolve[{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], s[0] == n, i[0] == 1,
r[0] == 0}, i, {t, times}]]); sol1 =
NMinimize[{(Max[
model1[\[Beta], \[Gamma], \[Mu]][#] & /@ Range[0, times]] -
reportedData)^2, {0 <= \[Beta] <= 1, 0 <= \[Gamma] <= 1,
0 <= \[Mu] <= 1}}, {\[Beta], \[Gamma], \[Mu]},
MaxIterations -> 1000]

(*Out[]= {4.87078*10^-11, {\[Beta] -> 0.88748, \[Gamma] ->
0.035164, \[Mu] -> 0.0405108}}*)


Visualization

Plot[
model1[\[Beta], \[Gamma], \[Mu]][t] /. sol1[[2]], {t, 0, 100}]


What is the advance to use LogLikelihood[] or loglik?

• The professor had said that the poisson log-likelihood gives better accuracy Commented Sep 2, 2023 at 13:20
• Because the estimate value of $\beta$=0.38 Commented Sep 2, 2023 at 14:29
• Try model with $\beta=0.38$ at $\gamma= 0.194, \mu -> 0.006$, You will get Max[model[0.38][#] & /@ Range[0, times]] of about 5178.92, or 0.047%. No, I don't believe your professor :) Commented Sep 2, 2023 at 14:43