Оценка SIR Поиск причины ошибки при оценке параметров

Question

Я пытаюсь подогнать модель распространения эпидемий SIR к текущим новым данным о случаях заболевания в странах. Для этого я использовал работу здесь: https://github.com/epimath/param-estimation-SIR. Основная идея заключалась в том, чтобы наилучшим образом подогнать кривую инфицированных SIR к новым данным случаев для этой конкретной c страны и рассчитать общее прогнозируемое количество случаев и дни, которые относятся к% 98 и% 95 от общего числа случаев. Проблема в том, что я выбираю Бразилию, Мексику или США. Это показывает, что это никогда не закончится. Мне любопытна причина. Мы будем благодарны за любую помощь о том, что можно сделать для решения этих несходящихся случаев.

Измените переменную selected_country с «Испания» на одну из этих трех стран (Бразилия, Мексика или США), чтобы воспроизвести результат, который заставляет меня спросить здесь.

PS Я знаю ограничения работы. Например, новый номер дела привязан к количеству тестов et c. Пожалуйста, игнорируйте эти ограничения. Я хотел бы увидеть, что необходимо для получения результата из следующего кода.

Вот некоторые результаты:

Испания (пример ожидаемого результата)

Турция (пример ожидаемого результата)

Франция (пример ожидаемого результата)

США (пример неожиданного результата)

Бразилия (пример неожиданного результата)

Я подозреваю, что параметр гаммы (скорость восстановления) слишком мал, что приводит к одинаковому количеству случаев для каждого день. Но я не смог go дальше и выяснил, чем это вызвано. (Я понял это, проверив переменную параметров, распечатав и изучив ее значения.)

Здесь вы можете найти мой код ниже.

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

df = pd.read_csv('https://github.com/owid/covid-19-data/raw/master/public/data/owid-covid-data.csv')  

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

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,
                                       'adaptive': False})
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()