[PYTHON] [Introduction au modèle SEIR] Essayez d'ajuster les données COVID-19 ♬

J'ai essayé d'ajuster les données d'infection de la nouvelle corona en utilisant le modèle SEIR. La méthode est la même que la précédente [Introduction pour minimiser] Analyse des données du modèle SEIR ♬. Les données d'infection ont été obtenues à partir du site de référence suivant. Les résultats suivants ont été obtenus dans [Introduction to matplotlib] Reading the end time from COVID-19 data ♬ pour les données de Wuhan. La valeur initiale est choisie et ajustée de manière à ce qu'une telle courbe puisse être obtenue, et le résultat est obtenu comme suit. Donc, cette fois, je voudrais faire des ajustements similaires pour chaque pays et examiner si de nouveaux résultats peuvent être obtenus. Les données sont obtenues à partir de ce qui suit, et sauf indication contraire, il s'agit de ** données en date du 24 mars 2020 **, mais ce qui précède est ** en date du 22 mars 2020 **. De plus, comme expliqué dans l'explication du code, veuillez noter que cette fois, l'ajustement est effectué pour les données d'infection et l'ajustement utilisant les données de guérison est exécuté par différents programmes. 【référence】 CSSEGISandData/COVID-19

Ce que j'ai fait

・ La fin du Japon ・ Explication du code ・ Adaptation des données nationales

・ La fin du Japon

J'ai essayé de l'écrire doucement comme mentionné ci-dessus, mais puisque le gouvernement métropolitain de Tokyo a dit que la ville serait fermée, je voudrais mentionner si la fin du Japon peut être vue en premier. De la conclusion, il semble que la fin sera sûrement vue dans un autre mois. En fait, dans les données ci-dessus du 22 mars, il y avait une situation où je me suis endormi pendant un moment et c'était le pic de l'infection, mais comme le nombre d'infections a augmenté linéairement ici, la saturation a été légèrement reportée. Donc, comme ajustement, j'ai obtenu le graphique suivant pour les données jusqu'au 24 mars.

L'appareillage détermine à la fois le nombre d'infections et la courbe de guérison. Les données pour le Japon sont un peu étranges du 53 au 59, mais elles évoluent généralement naturellement. Donc, j'ai une belle courbe comme suit. En d'autres termes, le résultat est qu'un pic apparaît dans environ 10 jours, une courbe de guérison culmine en environ un mois et le reste se termine doucement.

Les paramètres obtenus sont

beta_const=0.159
lp=1.97e-23
gamma=52.7
N=1474
R0 = 4.45

・ Explication du code

Le code complet est ci-dessous. Les deux raisons sont que la structure de données fournie a changé et que les données actuelles doivent être lues individuellement comme dans (2). En outre, dans (2), la fonction d'ajustement a été modifiée de la fonction normale la plus probable à l'erreur quadratique normale afin qu'elle s'adapte plus facilement à la courbe. ①collective_particles/fitting_COVID.py ②collective_particles/fitting_COVID2.py

Le code simple est expliqué ci-dessous. Utilisez le Lib. C'est la somme des articles des deux articles précédents.

#include package
import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize
import matplotlib.pyplot as plt
import pandas as pd

Ce qui suit est identique au dessin de données infectées.

#Lisez les données CSV avec les pandas.
data = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv')
data_r = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Recovered.csv')
data_d = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Deaths.csv')
    
confirmed = [0] * (len(data.columns) - 4)
confirmed_r = [0] * (len(data_r.columns) - 4)
confirmed_d = [0] * (len(data_d.columns) - 4)
days_from_22_Jan_20 = np.arange(0, len(data.columns) - 4, 1)

#City
city = "Hubei" #Données de Wuhan

#Données de processus
t_cases = 0

for i in range(0, len(data), 1):
    #if (data.iloc[i][1] == city): #for country/region
    if (data.iloc[i][0] == city):  #for province:/state  
        print(str(data.iloc[i][0]) + " of " + data.iloc[i][1])
        for day in range(4, len(data.columns), 1):
            confirmed[day - 4] += data.iloc[i][day] -  data_r.iloc[i][day] 
            confirmed_r[day - 4] += data_r.iloc[i][day]
            confirmed_d[day - 4] += data_d.iloc[i][day]

Utilisez le modèle SEIR.

#define differencial equation of seir model
def seir_eq(v,t,beta,lp,ip,N0):
    N=N0 #int(26749*1) #58403*4
    a = -beta*v[0]*v[2]/N
    b = beta*v[0]*v[2]/N-(1/lp)*v[1] #N
    c = (1/lp)*v[1]-(1/ip)*v[2]
    d = (1/ip)*v[2]
    return [a,b,c,d]

Ici, l'image obtenue change en fonction de la valeur de N0. C'est bien de l'adapter, mais cette fois, ce n'est pas fait. Lors de l'installation de Hubei, j'ai fait référence aux valeurs suivantes. Cependant, pour ** N et S0, le nombre réel de personnes infectées (valeurs déplacées) a été adopté en partant de l'idée que les résidents ne sont pas tous les participants du jeu mais la cible de l'infection **. 【référence】 ・ J'ai essayé de prédire le comportement du nouveau virus corona avec le modèle SEIR. Réglage de la valeur initiale

#solve seir model
N0 = 70939 #Hubei
N,S0,E0,I0=int(N0*1.),N0,int(318),int(389) #N is total population, S0 initial value of fresh peaple
ini_state=[S0,E0,I0,0] #initial value of diff. eq.
beta,lp,ip=1, 2, 7.4 
t_max=60 #len(days_from_22_Jan_20)
dt=0.01
t=np.arange(0,t_max,dt)
obs_i = confirmed

Une fonction qui renvoie une valeur d'évaluation pour un paramètre d'entrée. Le v [2] suivant; renvoie le nombre d'infections.

#function which estimate i from seir model func 
def estimate_i(ini_state,beta,lp,ip,N):
    v=odeint(seir_eq,ini_state,t,args=(beta,lp,ip,N))
    est=v[0:int(t_max/dt):int(1/dt)]
    return est[:,2]

La fonction objectif pour la minimisation ci-dessous (ce qui suit utilise la fonction la plus probable).

#define logscale likelihood function
def y(params):
    est_i=estimate_i(ini_state,params[0],params[1],params[2],params[3])
    return np.sum(est_i-obs_i*np.log(np.abs(est_i)))

Trouvez la valeur minimale avec minimiser @ scipy et calculez le nombre de reproduction de base.

#optimize logscale likelihood function
mnmz=minimize(y,[beta,lp,ip,N],method="nelder-mead")
print(mnmz)
#R0
beta_const,lp,gamma_const,N = mnmz.x[0],mnmz.x[1],mnmz.x[2],mnmz.x[3] #Taux d'infection, temps d'attente d'infection, taux d'élimination (taux de récupération)
print(beta_const,lp,gamma_const,N)
R0 = N*beta_const*(1/gamma_const)
print(R0)

Affichage graphique des résultats des calculs.

#plot reult with observed data
fig, (ax1,ax2) = plt.subplots(2,1,figsize=(1.6180 * 4, 4*2))
lns1=ax1.plot(obs_i,"o", color="red",label = "data")
lns2=ax1.plot(confirmed_r,"*", color="green",label = "data_r")
lns3=ax1.plot(estimate_i(ini_state,mnmz.x[0],mnmz.x[1],mnmz.x[2],mnmz.x[3]), label = "estimation")
lns0=ax1.plot(t,odeint(seir_eq,ini_state,t,args=(mnmz.x[0],mnmz.x[1],mnmz.x[2],mnmz.x[3])))
lns_ax1 = lns1+lns2 + lns3 + lns0
labs_ax1 = [l.get_label() for l in lns_ax1]
ax1.legend(lns_ax1, labs_ax1, loc=0)
ax1.set_ylim(0,N0)

Ce qui suit affiche un graphique à longue portée pour que la fin soit visible.

• Affiché jusqu'à 200 jours ici En conséquence, le résultat similaire à celui de Wuhan montré au début peut être obtenu. Cependant, l'image de Wuhan est une image jusqu'au 60e.

t_max=200 #len(days_from_22_Jan_20)
t=np.arange(0,t_max,dt)

lns4=ax2.plot(obs_i,"o", color="red",label = "data")
lns5=ax2.plot(confirmed_r,"*", color="green",label = "recovered")
lns6=ax2.plot(t,odeint(seir_eq,ini_state,t,args=(mnmz.x[0],mnmz.x[1],mnmz.x[2],mnmz.x[3])))
ax2.legend(['data','data_r','Susceptible','Exposed','Infected','Recovered'], loc=1)
ax2.set_title('SEIR_b{:.2f}_ip{:.2f}_gamma{:.2f}_N{:.2f}_E0{:d}_I0{:d}_R0{:.2f}'.format(beta_const,lp,gamma_const,N,E0,I0,R0))
ax1.grid()
ax2.grid()
plt.savefig('./fig/SEIR_{}_b{:.2f}_ip{:.2f}_gamma{:.2f}_N{:.2f}_E0{:d}_I0{:d}_R0{:.2f}_.png'.format(city,beta_const,lp,gamma_const,N,E0,I0,R0)) 
plt.show()
plt.close()

En revanche, lors de l'ajustement simultané du nombre d'infections et du nombre de cures, la fonction d'évaluation et la fonction objective au milieu ont été modifiées comme suit.

#function which estimate i from seir model func 
def estimate_i(ini_state,beta,lp,ip,N):
    v=odeint(seir_eq,ini_state,t,args=(beta,lp,ip,N))
    est=v[0:int(t_max/dt):int(1/dt)] 
    return est[:,2],est[:,3] #v0-S,v1-E,v2-I,v3-R
    
#define logscale likelihood function
def y(params):
    est_i_2,est_i_3=estimate_i(ini_state,params[0],params[1],params[2],params[3])
    return np.sum((est_i_2-obs_i_2)*(est_i_2-obs_i_2)+(est_i_3-obs_i_3)*(est_i_3-obs_i_3))

Dans ce cas, il est possible d'ajuster uniquement le nombre d'infections ou uniquement le nombre de cures.

・ Adaptation des données nationales

· Iran L'Iran, qui a un taux de guérison relativement élevé, sature plus rapidement que le Japon et est susceptible de se terminer, mais le nombre d'infections continue d'augmenter et le taux de guérison est lent, donc cela peut prendre plus de temps.

·Corée D'une manière ou d'une autre, la distorsion est devenue plus perceptible par rapport aux données précédentes et la crédibilité des données a été brisée. Pourtant, le taux de guérison augmente globalement et la fin semble être relativement précoce.

·Italie Le taux de mortalité est extrêmement élevé et dangereux, mais les données sont très claires et sophistiquées. Par conséquent, il était facile à installer et je me suis bien mis. Le taux de guérison est toujours d'environ 10%, et il semble que cela prendra un certain temps (le pic de guérison est de 2 mois), mais si le pic est observé après cela, la fin est proche.

Ce qui suit est la période d'expansion initiale, et il semble que la précision de l'ajustement est médiocre et qu'elle n'est pas fiable, mais je la posterai. ·Espagne Il se lève de plus en plus, mais c'est difficile à égaler. ·Allemagne Cela n'est pas encore terminé et le nombre de traitements est trop petit en premier lieu, il semble donc que ce soit dans la période d'expansion initiale.

Résumé

・ J'ai essayé d'ajuster les données COVID-19 ・ Le Japon, l'Iran et la Corée du Sud devraient voir une fin après la Chine. ・ Les autres pays d'Europe et les États-Unis sont encore en période d'expansion et il est difficile de prévoir la fin.

・ J'ai parlé grossièrement cette fois, mais j'aimerais faire une analyse plus détaillée tout en regardant les mises à jour quotidiennes.

prime

Le montage simultané de Wuhan est le suivant Les données de Wuhan ne semblent pas être expliquées par ce modèle (> <) données brutes Ajuster uniquement les données d'infection Données de guérison uniquement adaptées Adaptation des données d'infection et de guérison