[PYTHON] [Introduction au modèle SIR] Prédire l'heure de fin de chaque pays avec l'ajustement des données COVID-19 ♬
Dans le modèle SEIR précédent, parmi les paramètres obtenus, des valeurs telles que $ lp \ falldotseq 1e ^ {-23} $ ont été obtenues. Dans un tel cas, il est possible que la perte de chiffres ou l'équation différentielle ne fonctionne pas en premier lieu, alors soyez prudent. Donc, cette fois, je vais changer cette pièce et essayer de l'adapter au modèle le plus simple.
Ce que j'ai fait
・ Un peu de théorie ・ Explication du code ・ La fin du Japon ・ Fin de chaque pays ・ Une petite discussion
・ Un peu de théorie
Le modèle SIR est l'équation simple d'évolution temporelle suivante.
{\begin{align}
\frac{dS}{dt} &= -\beta \frac{SI}{N} \\
\frac{dI}{dt} &= \beta \frac{SI}{N} -\gamma I \\
\frac{dR}{dt} &= \gamma I \\
\end{align}
}
La première équation est que la diminution du nombre de personnes non infectées est proportionnelle au nombre de personnes infectées et au nombre de personnes non infectées, et est inversement proportionnelle au nombre total de personnes cibles. Diminution. En tant que β / N, il peut être considéré comme le taux d'infection par habitant. La deuxième équation est la formule pour augmenter le nombre de personnes infectées. L'augmentation est la différence entre l'afflux d'individus non infectés et d'individus guéris. La troisième équation est une formule selon laquelle un certain pourcentage de personnes infectées sont guéris, et ce coefficient γ est un taux de guérison appelé taux d'élimination. 【référence】 ・ [Cette infection se propage-t-elle ou converge-t-elle: la signification physique et la détermination du nombre de reproductions R](https://rad-it21.com/%E3%82%B5%E3%82%A4%E3%82%A8] % E3% 83% B3% E3% 82% B9 / kado-shinichiro_20200327 /)
Ici, la deuxième équation peut être transformée comme suit.
\frac{dI}{dt} = \gamma (\frac{\beta}{\gamma} \frac{S}{N} - 1)I\\
Autrement dit, si $ S / N $ est considéré comme constant, la solution suivante peut être formellement obtenue.
- On peut considérer que l'évolution du temps se fait dans un temps suffisamment court.
I = I_0\exp(\gamma (\frac{\beta}{\gamma} \frac{S}{N} - 1)t)\\
Au début de l'infection, c'est $ S_0 \ Fallingdotseq N $, auquel cas la solution suivante peut être obtenue.
{\begin{align}
I &= I_0\exp(\gamma (R_0 - 1)t)\\
ici\\
R_0 &= \frac{\beta}{\gamma}
\end{align}
}
Ici, dans la formule de $ I $, le nombre de personnes infectées diminue exponentiellement lorsque $ R_0 <1 $, et inversement augmente exponentiellement lorsque $ R_0> 1 $. Par conséquent, ce ** $ R_0 $ est appelé le nombre de reproduction de base ** et sert de guide pour déterminer si une transmission d'infection se produira ou non. De plus, dans la référence ci-dessus,
R = R_0\frac{S}{N}
Est défini comme ** le nombre de reproductions effectives **. Ce $ R $ devient plus petit à mesure que $ \ frac {S} {N} $ devient plus petit à mesure que l'infection progresse. Autrement dit, comme vous pouvez le voir à partir de la solution formelle ci-dessus, lorsque $ R_0> 1 $, $ R> 1 $ est toujours présent même au stade initial de l'expansion et que l'infection se propage, mais finalement $ \ frac {S} {N} $ est petit. Par conséquent, c'est un indicateur que la fin commence lorsque $ R <1 $ est réalisé. D'un autre côté, si le côté gauche de $ \ frac {dI} {dt} = 0 $ = 0, l'infection culmine et le nombre de personnes non infectées à ce moment est $ S_p $,
{\begin{align}
R_0 &= \frac{\beta}{\gamma}\\
& = \frac{N}{S_p} \\
\end{align}
}
Est établi, à ce moment
{\begin{align}
R &= R_0\frac{S_p}{N}\\
&= 1
\end{align}
}
Il s'avère que c'est vrai. Autrement dit, ** le nombre de reproductions efficaces est R = 1 au pic de l'infection, et la transmission de l'infection se terminera par la suite. ** ** Cette fois, je vais regarder la prédiction finale en traçant ce ** nombre de reproduction effectif ** également.
・ Explication du code
Le code complet est ci-dessous. ・ Collective_particles / fit_SIR_COVID2.py Le Lib utilisé est le suivant identique à la dernière fois.
import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize
import matplotlib.pyplot as plt
import pandas as pd
Données à l'avance à partir du site COVID-19 Est en téléchargement.
#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_covid19_confirmed_global.csv')
data_r = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv')
data_d = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.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)
L'ordre d'enregistrement des données a été correctement organisé la dernière fois, mais cette fois, l'ordre des pays et des régions des trois données n'est pas le même, alors obtenez chacun comme suit et calculez le nombre d'infections existantes pour les données sur le nombre d'infections Pour ce faire, j'essaie de l'acquérir après avoir acquis les données du numéro de guérison. En outre, si les déclarations sont utilisées de manière interchangeable dans chaque pays et région. De plus, le nombre total de guérisons et le nombre total de décès requis pour le taux de mortalité et le taux de guérison ne sont pas calculés cette fois, ils sont donc commentés.
#City
city = "Japan"
#Données de processus
t_cases = 0
for i in range(0, len(data_r), 1):
if (data_r.iloc[i][1] == city): #for country/region
#if (data_r.iloc[i][0] == city): #for province:/state
print(str(data_r.iloc[i][0]) + " of " + data_r.iloc[i][1])
for day in range(4, len(data.columns), 1):
confirmed_r[day - 4] += data_r.iloc[i][day]
#t_recover += data_r.iloc[i][day]
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] - confirmed_r[day - 4]
for i in range(0, len(data_d), 1):
if (data_d.iloc[i][1] == city): #for country/region
#if (data_d.iloc[i][0] == city): #for province:/state
print(str(data_d.iloc[i][0]) + " of " + data_d.iloc[i][1])
for day in range(4, len(data.columns), 1):
confirmed_d[day - 4] += data_d.iloc[i][day] #fro drawings
#t_deaths += data_d.iloc[i][day]
Cette fois, nous utiliserons le modèle SIR. La fonction d'évaluation renvoie également le nombre non infecté S, le nombre infecté I et le nombre guéri R.
#define differencial equation of sir model
def sir_eq(v,t,beta,gamma):
a = -beta*v[0]*v[1]
b = beta*v[0]*v[1] - gamma*v[1]
c = gamma*v[1]
return [a,b,c]
def estimate_i(ini_state,beta,gamma):
v=odeint(sir_eq,ini_state,t,args=(beta,gamma))
est=v[0:int(t_max/dt):int(1/dt)]
return est[:,0],est[:,1],est[:,2] #v0-S,v1-I,v2-R
La fonction objectif est une distribution normale. Ici, la somme linéaire est utilisée pour que le nombre d'infections I et le nombre de guérisons R puissent être ajustés en même temps.
#define logscale likelihood function
def y(params):
est_i_0,est_i_1,est_i_2=estimate_i(ini_state,params[0],params[1])
return np.sum(f1*(est_i_1-obs_i_1)*(est_i_1-obs_i_1)+f2*(est_i_2-obs_i_2)*(est_i_2-obs_i_2))
Entrez le nombre total d'infections dans le pays (ou la région) sélectionné ci-dessus dans N0. De plus, f1 et f2 déterminent sur quelle courbe l'ajustement doit être effectué. N, S0, I0, R0 sont les valeurs initiales de chacun, mais ce sont également des constantes qui doivent être ajustées en fonction de la condition d'ajustement. Les valeurs initiales de bêta et de gamma sont similaires dans la plupart des pays, mais il semble préférable de ressaisir les valeurs une fois obtenues et de calculer récursivement.
- Vous devriez également faire cet ajustement de valeur initiale, mais cette fois nous ne le ferons pas
#solve seir model
N0 = 1517
f1,f2 = 1,1 #data & data_r fitting factor
N,S0,I0,R0=int(N0*1.5),int(N0*1.5),int(10),int(0)
ini_state=[S0,I0,R0]
beta,gamma = 1.6e-6, 1e-3 #La valeur initiale ne converge que si elle est proche dans une certaine mesure
t_max=len(days_from_22_Jan_20) #S'adapter à une certaine plage de données
dt=0.01
t=np.arange(0,t_max,dt)
Je vais essayer de l'afficher ci-dessous, notamment en vérifiant la valeur initiale.
plt.plot(t,odeint(sir_eq,ini_state,t,args=(beta,gamma)))
plt.legend(['Susceptible','Infected','Recovered'])
obs_i_1 = confirmed
obs_i_2 = confirmed_r
plt.plot(obs_i_1,"o", color="red",label = "data_1")
plt.plot(obs_i_2,"o", color="green",label = "data_2")
plt.legend()
plt.pause(1)
plt.close()
Faites un essayage. J'ai également utilisé minimiser @ scipy cette fois. Calculez le nombre de reproduction de base R0. Puisque le nom de la variable chevauche le numéro de guérison initial ci-dessus, il est en minuscules.
#optimize logscale likelihood function
mnmz=minimize(y,[beta,gamma],method="nelder-mead")
print(mnmz)
#R0
#N_total = S_0+I_0+R_0
#R0 = N_total*beta_const *(1/gamma_const)
beta_const,gamma = mnmz.x[0],mnmz.x[1] #Taux d'infection, taux d'élimination (taux de guérison)
print(beta_const,gamma)
r0 = N*beta_const*(1/gamma)
print(r0)
Représentez-le ci-dessous. Le sous-axe de ax3 et ax4 est ajouté pour tracer le nombre de reproduction effectif.
#plot reult with observed data
fig, (ax1,ax2) = plt.subplots(2,1,figsize=(1.6180 * 4, 4*2))
ax3 = ax1.twinx()
ax4 = ax2.twinx()
lns1=ax1.plot(confirmed,"o", color="red",label = "data_I")
lns2=ax1.plot(confirmed_r,"o", color="green",label = "data_R")
est_i_0,est_i_1,est_i_2=estimate_i(ini_state,mnmz.x[0],mnmz.x[1])
lns3=ax1.plot(est_i_1, label = "estimation_I")
lns0=ax1.plot(est_i_2, label = "estimation_R")
lns_1=ax3.plot((S0-est_i_1-est_i_2)*r0/N,".", color="black", label = "effective_R0")
ax3.set_ylim(0,r0+1)
lns_ax1 = lns1+lns2 + lns3 + lns0 +lns_1
labs_ax1 = [l.get_label() for l in lns_ax1]
ax1.legend(lns_ax1, labs_ax1, loc=0)
ax1.set_ylim(0,N0)
ax1.set_title('SIR_{} f1_{:,d} f2_{:,d};N={:.0f} S0={:.0f} I0={:.0f} R0={:.0f} R={:.2f}'.format(city,f1,f2,N,S0,I0,R0,(S0-est_i_1[-1]-est_i_2[-1])*r0/N))
ax1.set_ylabel("Susceptible, Infected, recovered ")
ax3.set_ylabel("effective_R0")
Les prévisions futures sont présentées ci-dessous.
t_max=200 #len(days_from_22_Jan_20)
t=np.arange(0,t_max,dt)
lns4=ax2.plot(confirmed,"o", color="red",label = "data")
lns5=ax2.plot(confirmed_r,"o", color="green",label = "data_r")
est_i_0,est_i_1,est_i_2=estimate_i(ini_state,mnmz.x[0],mnmz.x[1])
lns6=ax2.plot(est_i_0, label = "estimation_S")
lns7=ax2.plot(est_i_1, label = "estimation_I")
lns8=ax2.plot(est_i_2, label = "estimation_R")
lns_6=ax4.plot((S0-est_i_1-est_i_2)*r0/N,".", color="black", label = "effective_R0")
lns_ax2 = lns4+lns5 + lns6 + lns7+ lns8 +lns_6
labs_ax2 = [l.get_label() for l in lns_ax2]
ax2.legend(lns_ax2, labs_ax2, loc=0)
ax4.set_ylim(0,r0+1)
ax2.set_ylim(0,S0)
ax2.set_title('SIR_{};b={:.2e} g={:.2e} r0={:.2f}'.format(city,beta_const,gamma,r0))
ax2.set_ylabel("Susceptible, Infected, recovered ")
ax4.set_ylabel("effective_R0")
ax1.grid()
ax2.grid()
plt.savefig('./fig/SIR_{}f1_{:,d}f2_{:,d};b_{:.2e}g_{:.2e}r0_{:.2f}N_{:.0f}S0_{:.0f}I0_{:.0f}R0_{:.0f}.png'.format(city,f1,f2,beta_const,gamma,r0,N,S0,I0,R0))
plt.show()
plt.close()
・ La fin du Japon
** Les données sont jusqu'au 27 mars 2020 **.
Les données sont les suivantes.
Les résultats sont les suivants.
Ce résultat semble convenir pour le moment, mais le rapport d'aujourd'hui indique que 200 personnes ont été infectées hier le 28 mars, il est donc hors de ce graphique et ne correspond pas. ..
J'ai posté ceci en bonus.
Les quantités obtenues sont les suivantes.
{\begin{align}
\beta &= 4.53e^{-5}\\
\gamma &= 1.58e^{-2}\\
N &= 2275\\
R0 &= 6.52\\
R(27.3.2020) &= 2.71\\
\end{align}
}
・ Fin de chaque pays
La situation du Hubei
Les dernières données (27.3.2020) sont les suivantes.
Les données de Wuhan ne peuvent pas être comparées comme suit, même si vous essayez d'ajuster les données à la fois du nombre d'infections et du nombre de traitements.
Cependant, comme vous pouvez le voir sur le graphique, on peut dire qu'il s'est terminé avec le nombre de reproduction effectif R = 0,01, 0,08.
{\begin{align}
\beta &= 3.19e^{-6}\\
\gamma &= 6.53e^{-2}\\
N &= 106462\\
R0 &= 5.20\\
R(27.3.2020) &= 0.01\\
\end{align}
}
Ce qui suit sont des raccords des deux en même temps, et le nombre maximal d'infections a considérablement changé.
{\begin{align}
\beta &= 2.34e^{-6}\\
\gamma &= 5.33e^{-2}\\
N &= 106462\\
R0 &= 4.67\\
R(27.3.2020) &= 0.08\\
\end{align}
}
Statut de la Corée du Sud
Les dernières données (27.3.2020) sont les suivantes.
Les données coréennes ne sont pas aussi mauvaises que celles de Wuhan, mais c'est presque la même chose, et même si vous essayez d'ajuster à la fois le nombre d'infection et les données de numéro de traitement, elles ne peuvent pas être comparées comme suit.
Cependant, comme vous pouvez le voir sur le graphique, le nombre de reproduction effectif R = 0,13 se termine.
{\begin{align}
\beta &= 2.25e^{-5}\\
\gamma &= 2.94e^{-2}\\
N &= 11365\\
R0 &= 8.68\\
R(27.3.2020) &= 0.13\\
\end{align}
}
Si vous essayez d'ajuster les deux, l'ajustement des données du nombre d'infection sera lâche et le pic d'infection se déplacera également vers la droite. Il sera peut-être possible de l'adapter s'il s'adapte mieux, mais cette fois je ne le poursuivrai plus.
{\begin{align}
\beta &= 2.11e^{-5}\\
\gamma &= 2.14e^{-2}\\
N &= 11365\\
R0 &= 11.17\\
R(27.3.2020) &= 0.18\\
\end{align}
}
・ Situation de l'Iran
Le prochain pays qui est susceptible de se terminer est l'Iran.
Les dernières données (27.3.2020) sont les suivantes.
Cependant, lorsque je suis arrivé ici et que le taux de guérison dépassait 30%, il semblait stagnant.
Même ainsi, le nombre effectif de reproductions est de 1,81, et il sera bientôt de 1, on peut donc dire que le nombre maximal d'infections est proche.
{\begin{align}
\beta &= 2.92e^{-6}\\
\gamma &= 4.94e^{-2}\\
N &= 62478\\
R0 &= 3.70\\
R(27.3.2020) &= 1.81\\
\end{align}
}
・ Statut de l'Italie
Le prochain pays préoccupant est l'Italie.
Les dernières données (27.3.2020) sont les suivantes.
Les données pour ce pays sont très solides et il est bon d'ajuster exactement les courbes d'infection et de guérison. Le résultat est le suivant, le nombre de reproduction effectif est de 4,07, mais il a fortement chuté et il semble qu'il se terminera dans environ 10 jours.
{\begin{align}
\beta &= 1.46e^{-6}\\
\gamma &= 1.91e^{-2}\\
N &= 143448\\
R0 &= 10.99\\
R(27.3.2020) &= 4.07\\
\end{align}
}
En regardant ces données, j'ai l'impression que les soins médicaux ne se sont pas effondrés et qu'une gestion telle que l'effet du blocus urbain est solide. (Les commentaires n'ont aucune base scientifique autre que des graphiques. Il s'agit d'une impression individuelle.)
・ Situation de l'Espagne
Ensuite, je montrerai la situation en Espagne où les décès augmentent rapidement. Conclusion Quand je l'écris, cela ressemble encore à une période d'expansion précoce, et le montage est indéfini.
- Situations qui ne peuvent pas être décidées maintenant avec ce modèle
Les dernières données (27.3.2020) sont les suivantes.
Ce qui suit montre le cas où N est plus petit et le cas où N est plus grand.
N etc. sont incertains, mais le nombre de reproductions efficaces est respectivement de 7,81 et 8,09, mais l'heure de fin (pic d'infection) est d'environ 80 à 90 jours (2 semaines à 1 mois plus tard), et si la valeur de pic est ignorée, On peut dire qu'ils sont dans un endroit similaire.
{\begin{align}
\beta &= 1.49e^{-6}\\
\gamma &= 1.37e^{-2}\\
N &= 127542\\
R0 &= 13.85\\
R(27.3.2020) &= 7.81\\
\end{align}
}
{\begin{align}
\beta &= 5.39e^{-7}\\
\gamma &= 1.94e^{-2}\\
N &= 354285\\
R0 &= 9.86\\
R(27.3.2020) &= 8.09\\
\end{align}
}
En fait, la même chose est vraie pour l'Allemagne, la France et la Suisse, mais j'ai essayé de l'adapter, mais en raison de l'augmentation rapide, cette fois (les paramètres ne sont pas définis), je vais donc l'écrire à nouveau.
・ Une petite discussion
Je ne sais pas si je peux vraiment l'accepter, alors j'écrirai un peu. Le plus gros inconvénient de ce raccord est qu'il est fixé sur N. Il semble qu'il soit possible de simuler correctement la situation où l'infection se transmet spontanément chez un sujet fermé. Cependant, de nombreux pays dans le monde n'ont pas pris de verrouillage complet, du moins jusqu'à présent, et il y a eu (en pourcentage) de nombreuses personnes infectées elles-mêmes. Au Japon, le nombre de ces patients ambulatoires augmente et il existe également des situations où des grappes se forment à partir de là. Donc, ce que j'ai écrit ici est un ajustement de la transition jusqu'à présent, et je ne peux pas prédire la situation future. Pour que ces prévisions soient correctes, il est supposé que la même situation perdurera, et il est nécessaire de pouvoir contrôler extrêmement la proportion de patients ambulatoires. En effet, s'il est possible d'écrire un modèle prenant en compte ces effets, il semble impossible de prédire avec ce modèle car il est difficile de prédire le nombre de patients ambulatoires qui change au hasard d'un moment à l'autre. Donc, je pense que l'ajustement sera possible si toute l'Europe et les Etats-Unis commencent à être vus, et il sera possible de discuter des paramètres obtenus.
Résumé
・ J'ai essayé d'ajuster les données COVID-19 avec le modèle SIR. ・ Il a été constaté que la fin de Wuhan était proche et que le pic d'infection était terminé en Corée du Sud. ・ On peut s'attendre à ce que l'Iran et l'Italie atteignent leur pic d'infection dans environ 10 jours relativement rapidement. ・ Au Japon, le pic d'infection est relativement tardif car il est bâclé, et on peut s'attendre à ce que le pic d'infection survienne dans 20 jours, selon les patients ambulatoires. ・ Si la situation actuelle évolue, on pense que le pic d'infection commencera à apparaître dans de nombreux pays en un mois.
・ Je veux surveiller quotidiennement et voir les changements de situation ・ Je souhaite évaluer les paramètres obtenus
prime
Maintenant, quand je regarde les données, elles ont été reflétées ** jusqu'au 28 mars 2020, alors j'ai essayé de les ajuster.
Cependant, les données d'hier sont aberrantes et la période de pointe de l'infection ne change pas beaucoup.
Les paramètres obtenus sont les suivants, et ils n'ont pas beaucoup changé par rapport à ce qui précède.
Cependant, ** le nombre de reproductions effectives R est petit, mais il augmente, alors soyez prudent **.
{\begin{align}
\beta &= 3.91e^{-5}\\
\gamma &= 1.65e^{-2}\\
N &= 2617\\
R0 &= 6.22\\
R(28.3.2020) &= 2.79\\
\end{align}
}
Je voudrais voir les tendances futures.
Recommended Posts