(Premier article) Une histoire sur le calcul numérique de la grippe et du nouveau coronavirus de la pneumonie avec Tensorflow

introduction

Abstrait

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.

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

image.png

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 $ $\beta=0.00218,\gamma=0.452$ Nous simulons comme. Dans un peu plus de détails sur les paramètres, le taux d'infection $ \ beta $ est le taux d'infection par jour et le taux de récupération $ \ gamma $ équivaut au temps qu'il faut pour ne plus être infecté.

À 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.

image.png

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 , une personne infectera deux personnes en moyenne. Je n'entrerai pas dans les détails, mais il est généralement connu que l'épidémie se terminera si $ R_0 <1 $$. Il peut être calculé en calculant $ côté droit = 0 $ dans la deuxième formule du modèle SIR.

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.

\beta=\frac{\gamma R_0}{S(0)}

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. $(S,I,R)=(11079999,1,0)$ Je vais calculer comme.

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.

image.png

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

Document de référence

Recommended Posts

(Premier article) Une histoire sur le calcul numérique de la grippe et du nouveau coronavirus de la pneumonie avec Tensorflow
Construire un environnement de calcul numérique avec pyenv et miniconda3
Une histoire sur la prédiction des préfectures à partir des noms de villes avec Jubatus
J'ai remplacé le calcul numérique de Python par Rust et comparé la vitesse
Une histoire sur l'automatisation du mahjong en ligne (Jakutama) avec OpenCV et l'apprentissage automatique
L'histoire de la création d'une caméra sonore avec Touch Designer et ReSpeaker
Créer un nouveau projet de calcul numérique Python
Une histoire sur le calcul de la vitesse d'une petite balle tombant tout en recevant la résistance de l'air avec Python et Sympy
Une histoire sur l'apprentissage automatique avec Kyasuket
Une histoire sur Python pop and append
Une histoire sur l'obtention du champ Atom (télégramme XML) de l'Agence météorologique avec une tarte à la râpe et de le tweeter
Une histoire accro aux variables globales et à la portée de Go
Une histoire sur l'implémentation d'un écran de connexion avec django
Une histoire sur la modification de Python et l'ajout de fonctions
Une histoire sur une erreur lors du chargement d'un modèle TensorFlow créé avec Google Colab localement