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"}]
```
[![enter image description here][1]][1]

```
*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]
```
[![enter image description here][2]][2]


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.
```
(*Adjust infected values for reporting rate*)adjustedInfectedValues = 
 infectedValues*0.6

(*Simulated reported cases*)
simulatedReportedCases = Round[adjustedInfectedValues]

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

(*Print the likelihood value*)
Print["Likelihood:", likelihood]
```
[![enter image description here][3]][3]


  [1]: https://i.sstatic.net/RaMEU.png
  [2]: https://i.sstatic.net/XoyCP.png
  [3]: https://i.sstatic.net/WvO9U.png