[Introduction à minimiser] Analyse des données avec le modèle SEIR ♬

En préparation de l'analyse des données du COVID-19, nous expliquerons comment générer le graphique suivant. Cette figure est une figure dans laquelle le modèle SEIR de maladie infectieuse est ajusté aux données données en utilisant minimiser @ scipy. SEIR_b6.88_ip1.22_gamma2.01_N762_E01_I00_R03.41_.png La méthode est presque comme ci-dessous. Les deux ont été utiles à bien des égards. 【référence】 ① Début du modèle mathématique des maladies infectieuses: aperçu du modèle SEIR par Python et introduction à l'estimation des paramètresModèle de prédiction mathématique des maladies infectieuses (modèle SIR): étude de cas (1)J'ai essayé de prédire le comportement du nouveau virus corona avec le modèle SEIR.

Ce que j'ai fait

・ Explication du code ・ Essayez le montage avec le modèle SIR et le modèle SEIR ・ Autres modèles

・ Explication du code

Le code complet est ci-dessous. ・ Collective_particles / minimiser_params.py Je laisserai l'explication à la référence ci-dessus, et ici je vais expliquer le code utilisé cette fois. Les données ont été utilisées à partir de la référence (2). Aussi, pour le code de minimisation, j'ai fait référence à Reference ①. La Lib à utiliser est la suivante. Cette fois, il est implémenté dans l'environnement keras-gpu de l'environnement conda de Windows 10, mais il est précédemment installé y compris scipy.

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

L'équation différentielle du modèle SEIR est définie ci-dessous. Ici, ce qu'il faut faire avec N sera décidé comme il convient, mais dans cette proposition, N = 763 personnes de l'article publié dans Reference (2), qui a été corrigé. L'équation différentielle du modèle SEIR est la suivante Ici, dans la référence (1), N est poussé dans la variable bêta et n'est pas exprimé explicitement, mais la formule suivante est plus facile à comprendre en termes de Wan. Ce qui suit est cité de la référence ③.

{\begin{align}
\frac{dS}{dt} &= -\beta \frac{SI}{N} \\
\frac{dE}{dt} &=  \beta \frac{SI}{N}  -\epsilon E \\
\frac{dI}{dt} &=  \epsilon E -\gamma I \\
\frac{dR}{dt} &=  \gamma I \\
\end{align} 
}

Si vous regardez ce qui suit tout en regardant ce qui précède, vous pouvez le comprendre car c'est juste.

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

Tout d'abord, une solution est obtenue sur la base de la valeur initiale sous forme d'équation différentielle ordinaire.

#solve seir model
N,S0,E0,I0=762,0,1,0
ini_state=[N,S0,E0,I0]
beta,lp,ip=6.87636378, 1.21965986, 2.01373496 #2.493913  , 0.95107715, 1.55007883
t_max=14
dt=0.01
t=np.arange(0,t_max,dt)
plt.plot(t,odeint(seir_eq,ini_state,t,args=(beta,lp,ip))) #0.0001,1,3
plt.legend(['Susceptible','Exposed','Infected','Recovered'])
plt.pause(1)
plt.close()

Les résultats sont les suivants: Si vous estimez et modifiez les paramètres ici, vous pouvez ajuster approximativement, mais c'est difficile pour les paramètres compliqués. SEIR_calc_.png Comme indiqué dans la référence (2), remplacez les données et affichez-les dans le graphique comme indiqué ci-dessous.

#show observed i
#obs_i=np.loadtxt('fitting.csv')
data_influ=[3,8,28,75,221,291,255,235,190,125,70,28,12,5]
data_day = [1,2,3,4,5,6,7,8,9,10,11,12,13,14]
obs_i = data_influ

plt.plot(obs_i,"o", color="red",label = "data")
plt.legend()
plt.pause(1)
plt.close()

Le résultat est le suivant, et comme il s'agit d'une sortie pour confirmation en cours de route, l'en-tête, etc. SEIR_data_.png Maintenant, le suivant est la partie estimation. Cette fonction d'évaluation calcule avec une équation différentielle en utilisant la valeur de l'argument et renvoie le résultat de l'évaluation.

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

La fonction de vraisemblance logarithmique est définie comme suit.

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

Utilisez la fonction de réduction de Scipy ci-dessous pour trouver la valeur maximale.

#optimize logscale likelihood function
mnmz=minimize(y,[beta,lp,ip],method="nelder-mead")
print(mnmz)

La sortie est la suivante

>python minimize_params.py
 final_simplex: (array([[6.87640764, 1.21966435, 2.01373196],
       [6.87636378, 1.21965986, 2.01373496],
       [6.87638203, 1.2196629 , 2.01372646],
       [6.87631916, 1.21964456, 2.0137297 ]]), array([-6429.40676091, -6429.40676091, -6429.40676091, -6429.4067609 ]))
           fun: -6429.406760912483
       message: 'Optimization terminated successfully.'
          nfev: 91
           nit: 49
        status: 0
       success: True
             x: array([6.87640764, 1.21966435, 2.01373196])

Retirez beta_const, lp, gamma_const du résultat du calcul comme suit, et calculez R0 ci-dessous. Il s'agit du nombre de reproductions de base, qui est un guide pour la propagation de l'infection.Si R0 <1, il ne se propage pas, mais si R0> 1, il se propage et est un indice indiquant le degré.

#R0
beta_const,lp,gamma_const = mnmz.x[0],mnmz.x[1],mnmz.x[2] #Taux d'infection, temps d'attente d'infection, taux d'élimination (taux de récupération)
print(beta_const,lp,gamma_const)
R0 = beta_const*(1/gamma_const)
print(R0)

Cette fois, c'est comme suit, et depuis R0 = 3,41, l'infection s'est propagée.

6.876407637532918 1.2196643491443309 2.0137319643699927
3.4147581501415165

Le résultat ainsi obtenu est affiché sur le graphique.

#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(estimate_i(ini_state,mnmz.x[0],mnmz.x[1],mnmz.x[2]), label = "estimation")
lns_ax1 = lns1+lns2
labs_ax1 = [l.get_label() for l in lns_ax1]
ax1.legend(lns_ax1, labs_ax1, loc=0)

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

Le résultat est comme ci-dessus.

・ Autres modèles

Lorsque je faisais ce calcul, j'ai trouvé un modèle SIR légèrement plus simple et un modèle légèrement plus compliqué, donc je vais les résumer. Le modèle SIR est le suivant car c'est un modèle qui ignore Exposed de S-E-I-R.

{\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} 
}

Vous pouvez comprendre la signification de $ R_0 $ en réécrivant comme suit.

{\begin{align}
\frac{dS}{dt} &= -R_0 \frac{S}{N} \frac{dR}{dt} \\
R_0 &= \frac{\beta}{\gamma}\\ 
\frac{dI}{dt} &= 0 \La personne infectée culmine à\\
R_0 &=  \frac{\beta}{\gamma}=  \frac{N}{S_0} \\
\end{align} 
}

En revanche, dans une situation compliquée telle qu'une transition entre la mort et la mort sans notification, ce sera comme suit. On suppose que les morts ne contribuent pas à l'infection.

Résumé

-L'adaptation des données peut maintenant être effectuée automatiquement par minimiser @ scipy. ・ J'ai pu ajuster les données réelles avec le modèle SEIR

・ La prochaine fois, j'aimerais analyser les données de chaque pays du COVID-19.

Recommended Posts

[Introduction à minimiser] Analyse des données avec le modèle SEIR ♬
[Introduction au modèle SEIR] Essayez d'ajuster les données COVID-19 ♬
Note de lecture: Introduction à l'analyse de données avec Python
20200329_Introduction à l'analyse de données avec Python 2nd Edition Personal Summary
Introduction à la modélisation statistique pour l'analyse des données Sélection du modèle GLM
Introduction à l'analyse de données avec Python P32-P43 [ch02 3.US Baby Names 1880-2010]
Introduction à l'analyse de données par Python P17-P26 [ch02 1.usa.gov données de bit.ly]
Analyse de données avec python 2
J'ai essayé l'analyse de données IRMf avec python (Introduction au décodage des informations cérébrales)
Analyse de données avec Python
[Livre technique] Introduction à l'analyse de données avec Python -1 Chapitre Introduction-
Introduction à RDB avec sqlalchemy Ⅰ
Introduction à l'apprentissage automatique avec scikit-learn - De l'acquisition de données à l'optimisation des paramètres
[Introduction à Python] Comment obtenir des données avec la fonction listdir
Introduction à RDB avec sqlalchemy II
Comment augmenter les données avec PyTorch
Analyse de données à partir de python (visualisation de données 1)
Introduction à l'analyse d'image opencv python
Analyse de données à partir de python (visualisation de données 2)
[Analyse du cours de l'action] Apprenez les pandas avec la moyenne Nikkei (004: Changer les données lues en moyenne Nikkei)
Liens vers des personnes qui commencent tout juste l'analyse de données avec python
Introduction à la modélisation statistique pour le modèle linéaire généralisé d'analyse de données (GLM)
De l'introduction de JUMAN ++ à l'analyse morphologique du japonais avec Python
"Introduction à l'analyse de données par modélisation statistique bayésienne à partir de R et Stan" implémenté en Python
[Introduction à WordCloud] Jouez avec le scraping ♬
Introduction à la modélisation statistique pour l'analyse des données Test de rapport de ressemblance GLM et asymétrie de test
[Introduction à Python3, jour 17] Chapitre 8 Destinations de données (8.1-8.2.5)
[Introduction à Python3, jour 17] Chapitre 8 Destinations de données (8.3-8.3.6.1)
Introduction au remplissage d'image Python Remplissage d'image à l'aide d'ImageDataGenerator
J'ai essayé l'analyse factorielle avec des données Titanic!
Convertir des données Excel en JSON avec python
[Introduction à Python] Utilisons foreach avec Python
[Introduction à Python3 Jour 19] Chapitre 8 Destinations de données (8.4-8.5)
Envoyer des données à l'API DRF avec Vue.js
Convertissez des données FX 1 minute en données 5 minutes avec Python
Essayez de convertir en données ordonnées avec les pandas
[Introduction à Pytorch] J'ai joué avec sinGAN ♬
Analyse de données à partir de python (pré-traitement des données-apprentissage automatique)
Comment lire les données de problème avec Paiza
J'ai essayé de prédire le comportement du nouveau virus corona avec le modèle SEIR.
[Introduction au Data Scientist] Bases de Python ♬
Introduction à la modélisation statistique pour l'analyse des données Élargissement de la gamme d'applications de GLM
[Didacticiel d'analyse Python dans la base de données avec SQL Server 2017] Étape 2: importer des données dans SQL Server à l'aide de PowerShell
Introduction à l'analyse des séries temporelles ~ Modèle d'ajustement saisonnier ~ Implémenté en R et Python
Une introduction à l'analyse de données à l'aide de Python - Pour augmenter le nombre de vues vidéo -
[Introduction à Python] Comment obtenir l'index des données avec l'instruction for
Comment créer des exemples de données CSV avec hypothèse
Markov Chain Artificial Brainless avec Python + Janome (1) Introduction à Janome
Chaîne de Markov artificielle sans cervelle avec Python + Janome (2) Introduction à la chaîne de Markov
Essayez d'agréger les données de musique doujin avec des pandas
Convertissez les données avec la forme (nombre de données, 1) en (nombre de données,) avec numpy.
J'ai essayé de sauvegarder les données avec discorde
[Introduction à StyleGAN2] Apprentissage indépendant avec 10 visages d'anime ♬
Introduction à Tornado (1): Framework Web Python démarré avec Tornado
J'ai essayé d'analyser les principaux composants avec les données du Titanic!
Enregistrez les données pour flasher avec la carte Nucleo STM32
J'ai essayé d'obtenir des données CloudWatch avec Python