# SIR Estimation Finding Cause of Error While Estimating Parameters With Python

Hello I had a question on stackoverflow and one user redirected me to here as this question is more related with statistics.

I'm trying to fit SIR Epidemics Spread Model to the current new case data of the countries. In order to do that I used the work here: https://github.com/epimath/param-estimation-SIR . Main Idea was to fit best possible SIR's Infected curve to the new case data for that specific country, and calculate total predicted case number and the days that belong to %98 and %95 of total cases. The problem is, when I select Brazil, Mexico or United States. It shows that it will never end. I am curious about the reason. Any help about what can be done to deal with this non converging cases would be appreciated.

Please change the selected_country variable from "Spain" to one of those three countries(Brazil, Mexico or United States) to reproduce the result that leads me to ask here.

P.S. I know the limitations of the work. For example, new case number is bound to the number of tests etc. Please ignore those limitations. I'd like to see what is needed to produce a result out of the following code.

Here are some outputs:

Spain (Expected Output Example)

Turkey (Expected Output Example)

France (Expected Output Example)

USA (Unexpected Output Example)

Brazil (Unexpected Output Example)

I suspect something that cause gamma(the rate of recovering) parameter too small which leads the same amount of cases for each day. But I couldn't go further and found out what causing that. (I understood that by checking paramests variable by printing and examining it's values.)

Here you can find my code below.

import scipy.optimize as optimize
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import poisson
from scipy.stats import norm
import json
from scipy.integrate import odeint as ode

import pandas as pd

from datetime import datetime
time_start = datetime.timestamp(datetime.now())
output = {"result": "error"}
error = False

def model(ini, time_step, params):
Y = np.zeros(3)  # column vector for the state variables
X = ini
mu = 0
beta = params[0]
gamma = params[1]

Y[0] = mu - beta * X[0] * X[1] - mu * X[0]  # S
Y[1] = beta * X[0] * X[1] - gamma * X[1] - mu * X[1]  # I
Y[2] = gamma * X[1] - mu * X[2]  # R

return Y

def x0fcn(params, data):
S0 = 1.0 - (data[0] / params[2])
I0 = data[0] / params[2]
R0 = 0.0
X0 = [S0, I0, R0]

return X0

def yfcn(res, params):
return res[:, 1] * params[2]

# cost function for the SIR model for python 2.7
# Marisa Eisenberg (marisae@umich.edu)
# Yu-Han Kao (kaoyh@umich.edu) -7-9-17

def NLL(params, data, times):  # negative log likelihood
params = np.abs(params)
data = np.array(data)
res = ode(model, x0fcn(params, data), times, args=(params,))
y = yfcn(res, params)
nll = sum(y) - sum(data * np.log(y))
# note this is a slightly shortened version--there's an additive constant term missing but it
# makes calculation faster and won't alter the threshold. Alternatively, can do:
# nll = -sum(np.log(poisson.pmf(np.round(data),np.round(y)))) # the round is b/c Poisson is for (integer) count data
# this can also barf if data and y are too far apart because the dpois will be ~0, which makes the log angry

# ML using normally distributed measurement error (least squares)
# nll = -sum(np.log(norm.pdf(data,y,0.1*np.mean(data)))) # example WLS assuming sigma = 0.1*mean(data)
# nll = sum((y - data)**2)  # alternatively can do OLS but note this will mess with the thresholds
#                             for the profile! This version of OLS is off by a scaling factor from
#                             actual LL units.
return nll

selected_location = 'Spain'
selected_df = df[df.location == selected_location].reset_index()
selected_df.date = pd.to_datetime(selected_df.date)

selected_df.date = pd.to_datetime(selected_df.date)
selected_df = selected_df[['date', 'new_cases']]
print(selected_df)

df = selected_df

optimizer = optimize.minimize(NLL, params, args=(data, times), method='Nelder-Mead',
options={'disp': False, 'return_all': False, 'xatol': 3.1201, 'fatol': 0.0001,
paramests = np.abs(optimizer.x)
iniests = x0fcn(paramests, data)

print('Paramests:')
print(paramests)

times_long = range(0, int(len(times) * 10))
start_day = df['date'][0]
dates_long = []
for i in range(0, int(len(times) * 10)):
dates_long.append(start_day + (np.timedelta64(1, 'D') * i))
# print(df)
# print(dates_long)
# sys.exit()
#### Re-simulate and plot the model with the final parameter estimates ####
xest = ode(model, iniests, times_long, args=(paramests,))
# print(xest)
est_measure = yfcn(xest, paramests)

# plt.plot(times, data, 'k-o', linewidth=1, label='Data')

json_dict = {}
time_end = datetime.timestamp(datetime.now())
json_dict['duration'] = time_end - time_start
json_df = pd.DataFrame()
json_df['dates'] = dates_long
json_df['new_cases'] = df['new_cases']
json_df['prediction'] = est_measure
json_df = json_df.fillna("")
json_df['cumulative'] = json_df['prediction'].cumsum()

json_df = json_df[json_df['prediction'] >= 1]

if error == True:
json_dict['result'] = 'error'
json_dict['message'] = error_message
json_dict['timestamp'] = datetime.timestamp(datetime.now())
json_dict['chart_data'] = json_df.drop(columns=['prediction'], axis=1)

else:
json_dict['result'] = 'success'
json_dict['day_for_95_percent_predicted_cases'] = \
json_df[json_df['cumulative'] > (json_df['cumulative'].iloc[-1] * 0.95)]['dates'].reset_index(drop=True)[0]
json_dict['day_for_98_percent_predicted_cases'] = \
json_df[json_df['cumulative'] > (json_df['cumulative'].iloc[-1] * 0.98)]['dates'].reset_index(drop=True)[0]

# json_dict['timestamp'] = str(f"{datetime.now():%Y-%m-%d %H:%M:%S}")
json_dict['timestamp'] = datetime.timestamp(datetime.now())

json_dict['chart_data'] = json_df.to_dict()

json_string = json.dumps(json_dict, default=str)
print(json_string)
output = json_string  # json string

plt.plot(json_df['dates'], json_df['prediction'], 'r-', linewidth=3, label='Predicted New Cases')
plt.bar(df['date'], data)
plt.axvline(x=json_dict['day_for_95_percent_predicted_cases'], label='(95%) '+str(json_dict['day_for_95_percent_predicted_cases'].date()),color='red')
plt.axvline(x=json_dict['day_for_98_percent_predicted_cases'], label='(98%) '+str(json_dict['day_for_98_percent_predicted_cases'].date()),color='green')
plt.xlabel('Time')
plt.ylabel('Individuals')
plt.legend()
plt.show()