Convert R code to Mathematica code of Poisson log-likelihood for the epidemic curve

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:

# INPUT
initial_state_values <- c(S = 762,
I = 1,
R = 0)

times <- seq(from = 0, to = 14, by = 0.1)

# SIR MODEL FUNCTION
sir_model <- function(time, state, parameters) {

with(as.list(c(state, parameters)), {

N <- S+I+R

lambda <- beta * I/N

# The differential equations
dS <- -lambda * S
dI <- lambda * S - gamma * I
dR <- gamma * I

# Output
return(list(c(dS, dI, dR)))
})
}



The translated Mathematica code

data = {{10, 19}, {17, 4}, {24, 5}, {32, 47}, {36, 38}, {38, 54}, {41,
84}, {42, 55}, {47, 116}, {50, 102}, {52, 126}, {58, 153}, {59,
164}, {62, 172}, {63, 177}, {67, 142}, {68, 158}, {69, 151}, {75,
155}, {80, 125}, {81, 116}, {84, 101}, {87, 106}, {88, 116}, {94,
85}, {98, 59}, {102, 67}, {105, 57}, {114, 59}, {117, 24}, {118,
39}, {120, 31}, {128, 36}, {134, 8}, {140, 7}, {142, 3}, {147,
20}, {150, 7}, {155, 2}, {159, 25}, {163, 2}, {186, 10}, {191,
17}, {200, 7}};

dataTable =
TableForm[data,
TableHeadings -> {None, {"Time", "Number Infected"}}]

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

times = data[[All, 1]];

(*Flu dataset*)
numberinfected = data[[All, 2]]

ClearAll[\[Beta], \[Gamma]]
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};
sol = ParametricNDSolve[{sireqns, initialConditions}, {s, i, r}, {t,
0, times}, {\[Beta], \[Gamma]}];

(*Extract the infected values from the solution*)
infectedValues = i /. sol;


Now I am trying to translate the following code from R to Mathematica, so I could calculate the parameters. But I face difficulties how to do it.

My main difficulty how can I create this code that calculates the Poisson log-likelihood

 LL <- sum(dpois(x = dat$$number_reported, lambda = 0.6 * output$$I[output$$time %in% dat$$time], log = TRUE))



So then I will try to find through this method the parameters with greater accuracy than with the (sum-of-squares) method.

# 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

# 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)))   # the "parameters" input argument of the loglik_function()

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)


• Since not all of use know R, we should post the ideas and the process of R instead of only post the code of R. Sep 1, 2023 at 8:47
• "But I face difficulties..." Specifying the difficulties would likely get you more and specific help.
– JimB
Sep 1, 2023 at 14:21

It looks like we need to compute $$\beta, \gamma$$ using optimization with NMinimize in the form

data = {{10, 19}, {17, 4}, {24, 5}, {32, 47}, {36, 38}, {38, 54}, {41,
84}, {42, 55}, {47, 116}, {50, 102}, {52, 126}, {58, 153}, {59,
164}, {62, 172}, {63, 177}, {67, 142}, {68, 158}, {69, 151}, {75,
155}, {80, 125}, {81, 116}, {84, 101}, {87, 106}, {88, 116}, {94,
85}, {98, 59}, {102, 67}, {105, 57}, {114, 59}, {117, 24}, {118,
39}, {120, 31}, {128, 36}, {134, 8}, {140, 7}, {142, 3}, {147,
20}, {150, 7}, {155, 2}, {159, 25}, {163, 2}, {186, 10}, {191,
17}, {200, 7}};

dataTable =
TableForm[data,
TableHeadings -> {None, {"Time", "Number Infected"}}];

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

times = data[[All, 1]];

(*Flu dataset*)
numberinfected = data[[All, 2]];

ClearAll[\[Beta], \[Gamma]]
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};
infVal =
ParametricNDSolveValue[{sireqns, initialConditions},
i[#] & /@ times, {t, 0, Last[times]}, {\[Beta], \[Gamma]}];

func[b_?NumericQ, g_?NumericQ] := Norm[infVal[b, g] - numberinfected];

sol1 = NMinimize[{func[b, g], {b > 0, g > 0}}, {b, g}]

(*Out[]= {84.3981, {b -> 0.180905, g -> 0.0748557}}*)


Visualization

Show[ListLinePlot[
Transpose[{times, Evaluate[infVal[b, g] /. sol1[[2]]]}]],
ListPlot[data, PlotStyle -> Red]]


• First of all, thank you for your answer. Is this way like the Poisson log-likelihood? Sep 1, 2023 at 20:10
• @AthanasiosParaskevopoulos If you need exactly the Poisson log-likelihood, then we have function LogLikelihood[] and PoissonDistribution[\[Mu]] to optimize parameters. Sep 2, 2023 at 6:02