[PYTHON] (Premier article) Une histoire sur le calcul numérique de la grippe et du nouveau coronavirus de la pneumonie avec Tensorflow
introduction
Abstrait
- Grippe
- Virus Corona
J'ai essayé d'analyser en utilisant un modèle mathématique.
Contexte
Cette fois, j'ai essayé de calculer numériquement le modèle des maladies infectieuses, mais comprenez bien que l'auteur n'est pas un professionnel des modèles mathématiques tels que le modèle des maladies infectieuses. Depuis que je me suis spécialisé en physique à l'université et à l'école supérieure, je n'ai abordé que dans une certaine mesure l'informatique. Le contexte de cet article était que j'étais curieux de savoir si je pouvais analyser le coronavirus, qui est un sujet brûlant en ce moment, de manière quantitative. Il y a eu le SRAS et le MERS dans les cas passés de maladies infectieuses, mais j'ai pensé qu'il serait utile d'étudier le type de modèle mathématique disponible pour prédire le temps de convergence des maladies infectieuses.
Modèle de maladie infectieuse
Modèle SIR
Il existe un modèle SIR comme équation modèle qui décrit le processus épidémique des maladies infectieuses. Il existe trois variables: sensible (suspecté d'être infecté), infecté (infecté) et rétabli (porteur d'immunité). Ce modèle a été publié dans un article de 1927 («A Contribution to the Mathematical Theory of Epidemics»). En particulier, il est souvent utilisé dans des domaines traitant de phénomènes familiers tels que la physique sociale.
S+I \rightarrow 2I \\
I \rightarrow R \\
Le processus élémentaire est constitué des deux processus ci-dessus. Le haut est le processus d'infection et le bas est le processus de récupération. Pensez-y comme une formule de réaction chimique. Lors de son incorporation dans une équation différentielle, la même logique de calcul que la théorie de la vitesse de réaction en chimie doit être utilisée.
Formule
Le modèle SIR est formulé comme suit.
\frac{dS(t)}{dt} = -\beta S(t)I(t) \\
\frac{dI(t)}{dt} = \beta S(t)I(t)-\gamma I(t)\\
\frac{dR(t)}{dt} = \gamma I(t) \\
Pour $ t $ à un moment donné, $ S (t) $ est le nombre d'infections suspectées, $ I (t) $ est le nombre de personnes infectées et $ R (t) $ est le nombre de porteurs immunitaires.
- Formule supérieure ➤ Le nombre de personnes soupçonnées d'être infectées diminue proportionnellement au nombre de contacts avec la personne infectée.
- Formule intermédiaire ➤ Le nombre de personnes infectées est proportionnel au nombre de contacts, mais il diminuera à mesure que le nombre de personnes qui guérissent augmente.
- Formule inférieure ➤ Si le nombre de personnes en convalescence augmente, le nombre de personnes infectées diminuera.
Ici, $ \ beta $ représente le taux d'infection et $ \ gamma $ représente le taux de récupération et le taux d'isolement. Aussi
-\beta I(t) \\
Équivalents à l'infectiosité.
Surtout pourquoi il peut être décrit comme l'équation différentielle ci-dessus http://www2.meijo-u.ac.jp/~tnagata/education/react/2019/react_03_slides.pdf Je pense que vous devriez vous y référer. C'est la même logique que lors de la formulation de la formule de vitesse de réaction pour la chimie du lycée.
Ici, je résous les équations différentielles simultanées ci-dessus, et quand j'ai vu cela, je me suis souvenu des équations de Lorenz que j'avais apprises à l'université. L'équation de Lorenz est une équation différentielle ordinaire (communément appelée ODE) dont on parle dans le domaine du chaos en physique. Je pense que je suis assez maniaque, donc je n'entrerai pas dans les détails, mais je vais juste écrire la formule.
\frac{dx}{dt} = -px+py \\
\frac{dy}{dt} = -xz+rx-y\\
\frac{dz}{dt} = xy-bz \\
J'avais l'habitude de faire cette simulation numérique pour une mission universitaire, alors j'ai pensé que je devrais faire la même simulation qu'à l'époque. À ce moment-là, je me souviens de l'avoir résolu avec Rungekutta, mais cette fois, j'essaierai d'utiliser TensorFlow pour étudier.
Qu'est-ce que TensorFlow
TensorFlow est une bibliothèque de logiciels open source pour le Machine Learning développée par Google. En un mot, c'est une bibliothèque Python pour l'analyse numérique à grande vitesse. Je ne ferai pas de Machine Learning cette fois, mais je l'utiliserai pour le calcul numérique. Je ne pense pas que vous puissiez entendre le tenseur, mais je comprends que c'est une extension de la matrice. Puisqu'il est hors de portée de cet article, je le posterai quelque part. La caractéristique de TensorFlow est que le code source peut être écrit facilement. D'autre part, la plupart du traitement des calculs est comme une boîte noire, de sorte que le traitement des calculs n'est pas clair.
Simulation numérique réelle (1)
Code source
Les paramètres, etc. seront expliqués dans la partie résultat.
sir_model.py
import numpy as np
import tensorflow.compat.v1 as tf
tf.disable_v2_behavior()
#Cela peut ne pas fonctionner selon la version de TensorFlow(Probablement à cause de la compatibilité)
#import tensorflow as tf
import matplotlib.pyplot as plt
import time
from scipy.integrate import solve_ivp
%matplotlib inline
def main():
x=tf.placeholder(dtype=tf.float64,shape=(3),name='x')
class SIR(object):
def __init__(self, beta=0.00218, gamma=0.452, **kwargs):
self.beta = float(beta)
self.gamma = float(gamma)
def __call__(self, t, x):
dx_dt = -self.beta * x[0]*x[1]
dy_dt = self.beta*x[0] *x[1] - self.gamma*x[1]
dz_dt = self.gamma*x[1]
dX_dt = tf.stack([dx_dt, dy_dt, dz_dt])
return dX_dt
f_sir = SIR()
fx =f_sir(None, x) # ODE
x0 = np.array([762, 1, 0], dtype=np.float64) # initial value
sess = tf.Session()
sess.run(tf.initializers.global_variables())
print('###value of f(x)###')
print(sess.run(fx, feed_dict={x:x0}))
dt=0.1
tstart=0.0
tend=30.0
ts=np.arange(tstart, tend+dt, dt) # 0.Sortie 30 jours en une seule étape
def f_sir_tf(t,xt):
return sess.run(fx, feed_dict={x: xt})
start_time =time.time() # measurement of time
sol_sir = solve_ivp(fun=f_sir_tf,
t_span=[tstart, tend], y0=x0, t_eval=ts)
integration_time_tf = time.time() - start_time
t_lo = sol_sir['t'] #Obtenez l'heure de chaque étape
x_lo = sol_sir['y'] #X de chaque étape(t)Obtenez la valeur de
#Calculer le temps de traitement
print("processing time (tf) %.5f" % integration_time_tf)
fig=plt.figure(figsize=(14,8))
ax=fig.add_subplot(111)
ax.plot(t_lo,x_lo[0],'*',color='blue',alpha=0.2,label="Susceptible",markersize=10)
ax.plot(t_lo,x_lo[1],'^',color='red',alpha=0.2,label="Infected",markersize=10)
ax.plot(t_lo,x_lo[2],'.',color='green',alpha=0.2,label="Recovered",markersize=10)
ax.grid()
ax.legend()
ax.set_xlabel('day')
ax.set_ylabel('number')
if __name__ == "__main__":
main()
résultat de la simulation
Pour cette simulation, je me suis référé à la p14 de la diapositive suivante.
http://www.actuaries.jp/lib/meeting/reikai20-7-siryo.pdf
Il s'agit d'un cas de grippe dans un internat enquêté en 1978. Sur un total de 763 personnes, $ S = 762 $ personnes et $ I = 1 $ personnes sont les conditions initiales. Aussi, pour le taux d'infection $ \ beta $ et le taux de récupération $ \ gamma $
À propos, le graphique de la relation entre $ S $ et $ I $ à ce moment est le suivant. Il converge vers l'origine avec l'image allant du bas à droite vers la gauche.
Interprétation
Le nombre de $ R $ verts récupérés n'est pas très important ici. Le plus important est le nombre de personnes infectées par le $ I $ rouge. Dans ce cas, à partir de l'état initial $ I = 1 $, $ S $ en contact avec la personne infectée devient $ I $, augmente rapidement, culmine une semaine plus tard et converge vers 0 après environ 14 à 15 jours. Vous pouvez voir comment ça se passe. Vous pouvez voir qu'il est en fait bien ajusté comme le graphique en cercle blanc sur la diapositive de référence p14.
Simulation numérique réelle (2)
Maintenant que nous avons calculé les valeurs numériques de la grippe, la suivante est "Coronavirus". Je n'ai pas expliqué la grippe, mais j'expliquerai à partir du "nombre de reproduction de base $ R_0 $" qui est un paramètre important lorsqu'il s'agit de maladies infectieuses.
Numéro de reproduction de base
La définition est «le nombre attendu d'infections secondaires produites par une personne infectée pendant toute la période d'infection». Par exemple, si $ R_0 = 2
Je vais citer l'URL suivante pour ce $ R_0 $. https://wired.jp/2020/01/26/scientists-predict-wuhans-outbreak-will-get-much-worse/
Un facteur d'incertitude majeur est l'infectiosité du 2019-nCoV. En fait, à quel point est-il contagieux? Dans le modèle de Reed, le nombre de personnes qu'une personne infectée peut infecter (reproduction de base du virus) est estimé à 3,6-4,0. En comparaison, le SRAS (syndrome respiratoire aigu sévère) est de 2 à 5 et la rougeole la plus contagieuse chez l'homme est de 12 à 18.
Calculons ce virus corona avec le nombre de reproduction de base $ R_0 = 3,4 $. Bien que le calcul soit omis, le taux d'infection $ \ beta $ peut être calculé par la formule suivante.
Considérons maintenant les conditions initiales. La ville de Wuhan, où le virus corona est répandu, est une grande ville de Chine avec 11,08 millions d'habitants.
Ensuite, vous obtiendrez $ \ beta = 2.19 \ times10 ^ {-8} $. Ici, j'ai mis $ \ gamma = 0.071428 $ (1/14 équivaut à 2 semaines). $ S (0) $ est 11 079 999.
résultat
Le code source est fondamentalement le même que pour la grippe. Seuls les paramètres ont été modifiés.
Interprétation
Les conditions initiales supposent le moment où la première personne infectée est trouvée. En réalité, je pense que c'est le 1er décembre 2019. Ce n'est que le début du mois de février, donc si vous pensez qu'environ deux mois se sont écoulés, c'est 60 jours. C'est vraiment effrayant quand j'y pense. .. .. Cette fois, j'ai utilisé un assez grand nombre de R = 3,4, donc je suppose le pire des cas. En fait, selon l'OMS (Organisation mondiale de la santé), il est de $ R_0 $ = 2,2 ~ 2,8, donc il devrait converger plus tôt. Le point clé est de savoir quand trouver la première personne infectée. En supposant qu'il s'agit du pic de personnes infectées (environ 100 jours dans le graphique), il est prévu qu'il faudra environ 40 jours pour que le nombre de personnes infectées converge vers 0 à partir d'ici. Comme beaucoup d'entre vous l'ont peut-être remarqué, ce modèle ne considère pas les morts. Par conséquent, le système S + I + R est conservé. Si vous voulez penser à un modèle plus réaliste, vous pouvez utiliser le modèle SEIR.
Résumé
Ce post
J'ai posté à Qiita pour la première fois. Cette fois, nous avons calculé numériquement le modèle SIR, qui est le plus élémentaire des «modèles d'infection». Les paramètres réellement utilisés ont été tirés de cas antérieurs, mais des mathématiques plus avancées sont nécessaires pour effectuer une estimation détaillée des paramètres. Si j'ai le temps, j'écrirai un tel article à l'avenir. En tant qu'individu, je m'intéresse à l'apprentissage automatique et à la programmation de concours, alors je pense créer un référentiel qui comprendra des mémorandums.
Perspective
- Calcul numérique du modèle SEIR (modèle incluant la période de latence dans SIR)
- SIR dans le modèle de réseau
- Jouez avec TensorFlow
- Lavez-vous sérieusement les mains
Document de référence
- https://www.ism.ac.jp/editsec/toukei/pdf/54-2-461.pdf
- http://www.orsj.or.jp/archive2/or56-12/or56-12_728.pdf
Recommended Posts