# Performing maximum likelihood estimation for SIR Model

The code first produces a plot of the dataset, then simulates the SIR model using the parameter values we found in the manual calibration to the full dataset of the total number infected .

(*total size of population*) n = 763;

(*days*)
tmax = 14;

(*Flu dataset*)
numberinfected = {3, 8, 26, 76, 225, 298, 258, 233, 189, 128, 68, 29,
14, 4};

pointsPlot =
ListPlot[numberinfected, PlotStyle -> Purple,
PlotTheme -> "Detailed",
FrameLabel -> {"Time (days)", "Number Infected"},
PlotLegends -> {"Total cases"}]

*Set parameter values*)
{\[Beta], \[Gamma]} = {1.7, 0.45};

(*Define the SIR model equations*)
sireqns = {s'[t] == -\[Beta]*s[t]*i[t]/n,
i'[t] == \[Beta]*s[t]*i[t]/n - \[Gamma]*i[t],
r'[t] == \[Gamma]*i[t]};

initialConditions = {s[0] == n, i[0] == 1, r[0] == 0};

(*Solve the differential equations*)
sol = NDSolve[{sireqns, initialConditions}, {s, i, r}, {t, 0, tmax}];

(*Extract the infected values from the solution*)
infectedValues = i /. sol[[1]];
fitPlot =
Plot[infectedValues[t], {t, 1, 14}, PlotStyle -> Red,
PlotStyle -> Red, PlotRange -> All,
PlotLegends -> {"Simulated infected"}];
Show[pointsPlot, fitPlot, PlotRange -> All]

Now, we want to calculate the likelihood of the model with these specific parameter values, i . e . the probability of observing these numbers of reported cases given our simulated numbers of infected people .

We are building up to calibrating the SIR model to our flu outbreak data from previous exercises using likelihood as a measure of the divergence between the model projections and the data. This time, even though we are looking at the same outbreak, the dataset only shows the reported cases, and we know that 60 % of flu cases are reported.

infectedValues*0.6

(*Simulated reported cases*)

(*Likelihood calculation*)
likelihood =
Product[Exp[-\[Lambda]] \[Lambda]^k/k!, {\[Lambda],
simulatedReportedCases}];

(*Print the likelihood value*)
Print["Likelihood:", likelihood]

• Usually likelihoods are L~Exp[-(yi-y[x[i]])^2/2], so you're missing a minus sign there, but I haven't ran the code, perhaps it's also something else. Commented Aug 29, 2023 at 11:45
• @HansOlo Unfortunately it does not work Commented Aug 29, 2023 at 11:59
• @HansOlo Let’s take it a step back, can you explain that to me? Commented Aug 29, 2023 at 12:39
• The immediate cause of the \$RecursionLimit error is that k is undefined and {\[Lambda], simulatedReportedCases} does not provide a list of \[Lambda] values. But the underlying issue is that you haven't given a mathematical definition of the likelihood (or the log of the likelihood). I suggest asking the question about the likelihood of this model at CrossValidated (stats.stackexchange.com) and then back here for implementation.
– JimB
Commented Aug 29, 2023 at 14:28

• ParametricNDSolve and NonlinearModelFit work for this cases.
Clear["Global`*"];
n = 763;
tmax = 14;
numberinfected = {3, 8, 26, 76, 225, 298, 258, 233, 189, 128, 68, 29,
14, 4};
pointsPlot =
ListPlot[numberinfected, PlotStyle -> Purple,
PlotTheme -> "Detailed",
FrameLabel -> {"Time (days)", "Number Infected"},
PlotLegends -> {"Total cases"}];
sireqns = {s'[t] == -β*s[t]*
i[t]/n, i'[t] == β*s[t]*i[t]/n - γ*i[t],
r'[t] == γ*i[t]};

initialConditions = {s[0] == n, i[0] == 1, r[0] == 0};
sol = ParametricNDSolve[{sireqns, initialConditions}, {s, i, r}, {t,
0, tmax}, {β, γ}];
nlm = NonlinearModelFit[numberinfected,
i[β, γ]@t /. sol, {β, γ}, t]
Show[ListPlot[numberinfected],
Plot[nlm[t], {t, 0, 14}, PlotStyle -> Orange]]
nlm["BestFitParameters"]

{β -> 1.66747, γ -> 0.444057}

• my problem is how I could calculate the likelihood of reported cases in the above situation Commented Aug 29, 2023 at 13:08
• @AthanasiosParaskevopoulos Yes, I solve the original problem,the newly edited problem I think we need time to understand. Commented Aug 29, 2023 at 13:16
• To be clearer I want to calculate the likelihood and design the reported cases Commented Aug 29, 2023 at 13:18
(*Provided data*)
numberinfected = {3, 8, 26, 76, 225, 298, 258, 233, 189, 128, 68, 29,
14, 4};

(*Parameters*)
{\[Beta], \[Gamma]} = {1.7, 0.45};
n = 763; (*Total population*)
underreportingFactor = 0.6;
tmax = 14;

(*Define the SIR model equations with underreporting*)
sireqns = {s'[t] == -\[Beta]*s[t]*i[t]/n,
i'[t] == (\[Beta]*s[t]*i[t]/n - \[Gamma]*i[t]),
r'[t] == \[Gamma]*i[t]};

(*Initial conditions*)
initialConditions = {s[0] == n, i[0] == 1, r[0] == 0};

(*Solve the SIR model equations*)
sirSol =
NDSolve[{sireqns, initialConditions}, {s, i, r}, {t, 0, tmax - 1}];

(*Calculate reported cases from the model*)
reportedCasesModel[t_] := (i[t] /. sirSol)*underreportingFactor;

(*Plot reported cases from the model and the actual dataset*)
modelPlot =
Plot[reportedCasesModel[t], {t, 0, tmax - 1},
PlotStyle -> {Purple, Dashed, Thick},
PlotLegends -> {"Reported cases "}];

dataPlot =
ListPlot[numberinfected, PlotStyle -> Red,
PlotLegends -> {"Total cases"}];

Show[modelPlot, dataPlot,
FrameLabel -> {"Time (days)", "Reported Cases"},
PlotLabel -> "Reported Cases: Model vs Data", PlotRange -> All]

I am still struggling to calculate likelihood