[PYTHON] [Analyse du modèle SIR] Transformez la formule pour déterminer γ et le nombre de reproduction effectif R ♬

Lorsque les dernières données sont naturellement extrapolées, la figure ci-dessous est obtenue. Je pense que les gens qui peuvent lire le graphique comprendront, mais je me suis rendu compte que la situation était plutôt mauvaise.

Donc, je voudrais parler de ce qui suit sans toucher à ce sujet. Les données traitées cette fois sont les données jusqu'au 3/4/2020 publiées sur le site de référence. Pour le nombre d'infections au Japon, la valeur lue à partir de la référence (2) a été utilisée. 【référence】 ①Novel Coronavirus (COVID-19) Cases, provided by JHU CSSE ② Nouveau virus Corona @ Nittele NEWS vu à partir des données et du graphique

Ce que j'ai fait

・ Transformation du coefficient du modèle SIR ・ Signification du modèle SIR ・ Explication des résultats ・ Explication du code

・ Transformation du coefficient du modèle SIR

Cette fois, je vais essayer de trouver le changement quotidien de $ \ gamma $ et le nombre de reproductions effectives en transformant le modèle SIR suivant comme suit.

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

Premièrement, la troisième équation peut être transformée comme suit, et puisque $ I $ et $ r $ ont des valeurs d'observation, $ \ gamma $ peut être obtenu.

{\begin{align}
\frac{1}{I}\frac{dr}{dt} &=  \gamma \\
\end{align} 
}

Ensuite, si vous remplacez ce $ \ gamma $ dans la deuxième équation

{\begin{align}
\frac{1}{\gamma I}\frac{dI}{dt} &= R-1\\
R &=  \frac{\beta}{\gamma} \frac{S}{N}\\
\end{align} 
}

Par conséquent, R est obtenu.

・ Signification du modèle SIR

En utilisant ces deux quantités, l'équation différentielle ci-dessus peut être réécrite comme suit.

{\begin{align}
\frac{dS}{dt} &= -\gamma R I \\
\frac{dI}{dt} &=  \gamma (R - 1) I \\
\frac{dr}{dt} &=  \gamma I \\
\end{align} 
}

Ici, $ \ gamma $ est un facteur lié à (taux de récupération, taux d'isolement et taux de mortalité) et signifie le taux auquel une personne infectée ne contribue pas à l'infection (supprimée). En revanche, $ R $ est appelé le nombre effectif de reproduction, qui est une quantité liée au nombre de personnes infectées par une personne infectée. Par conséquent, la première équation montre que plus $ R $ et $ \ gamma $ sont grands, plus l'infection progresse proportionnellement à la personne infectée. C'est la raison pour laquelle la propagation de l'infection augmente (de façon exponentielle) dans les calculs dits murins. D'un autre côté, la troisième formule signifie que plus le $ \ gamma $ est grand, plus il est guéri proportionnellement à la personne infectée. La deuxième équation est proportionnelle à $ (R-1) $, et il est montré que la valeur positive ou négative de cette valeur détermine si le nombre de personnes infectées augmente ou diminue. De plus, étant proportionnel au nombre de personnes infectées, on constate que le changement est exponentiel. Et ce changement est également proportionnel à $ \ gamma $.

・ Explication des résultats

La première chose qui attire l'attention cette fois est que le nombre d'infections (en rouge) a considérablement changé depuis le 64e jour (il y a 3 avril à 9 jours). Et le but principal de cette fois est de demander $ \ gamma $ et $ R $ pour cette tendance au Japon. Cependant, ce qui précède et surtout, les données quotidiennes fluctuent considérablement, le problème était donc de savoir comment y faire face. Si vous regardez les données, vous pouvez voir que la courbe du nombre de cures est terrible. J'ai pensé à cela comme un soi-disant lissage (3 points, 5 points, 7 points ...), mais c'est pire que cela, et dans son ensemble, cela change presque linéairement, donc je l'adapte avec une fonction exponentielle. J'ai décidé. De même, le nombre d'infections a également été ajusté, mais cela n'a été fait que pour la partie droite. Et cette fois, il n'a été utilisé que pour l'extrapolation présentée au début, pas pour ce calcul de coefficient. La courbe $ I / (R + D) $ représente le rapport entre les personnes infectées et les non-contributeurs. Ce nombre est le nombre auquel le taux de guérison augmente, et lorsque ce nombre diminue et tombe en dessous de 1, la fin est enfin visible. Au Japon, on s'attendait à ce que le nombre soit une fois de 3 ou moins, puis de passer à 1 tel quel, mais après le 64e jour ci-dessus, ce nombre a également commencé à augmenter de manière significative. De plus, il avait légèrement diminué jusque-là, mais il a commencé à augmenter linéairement, ce qui est vraiment surprenant. Le graphique du bas montre les $ \ gamma $ et $ R $ que je voulais trouver cette fois, le nombre de cures quotidiennes (moyenne sur 5 jours) utilisées dans le calcul et leurs courbes approximatives. $ \ Gamma $ divise le nombre de cicatrisation quotidien (différenciation de la courbe de cicatrisation) par le nombre d'infections. Cependant, étant donné que la courbe du nombre de guérison a une grande fluctuation, la courbe approximative a été utilisée ici. Même avec cela, il y avait une fluctuation considérable due à la fluctuation de la courbe du nombre d'infections. Cependant, des changements ont été observés chaque jour, et bien que ce nombre ait également augmenté régulièrement jusqu'au 64e jour, il a diminué de manière linéaire à partir de là. D'un autre côté, $ R $ est inversement proportionnel à ce $ \ gamma $, et le résultat fonctionne de cette façon. Cependant, la valeur la plus basse est d'environ 3 reproductions effectives et la position est légèrement décalée vers la gauche. Dans tous les cas, cette surface a complètement augmenté, et on peut voir que le nombre de reproductions effectives augmente à environ 10. Les dernières valeurs sont les suivantes. Le nombre de reproductions efficaces est une comparaison avec 1, qui est une valeur très grande (une valeur à laquelle l'infection se propage rapidement).

{\begin{align}
\gamma &= 1.6e^{-2} \\
R &= 10 \\
\end{align} 
}

La contribution réelle dépend de $ \ gamma (R-1) $, et le graphique ci-dessus montre que la contribution de la diminution de $ \ gamma $ et la contribution de l'augmentation de $ R-1 $ sont équilibrées. Par conséquent, lorsque ce montant est représenté graphiquement, il devient comme suit. En conséquence, la contribution de $ R $ était importante, et une fois qu'elle a diminué à environ 0,05, elle a augmenté après le 64e jour, et maintenant on peut voir qu'elle est passée à environ 0,14. On peut dire que la probabilité d'infection a triplé au cours des deux dernières semaines.

・ Explication du code Je mets tout le code ci-dessous collective_particles/fitting_gamma_R.py Cette fois, je n'expliquerai que les parties principales. Premièrement, la lecture et le traitement des données sont les suivants.

city = "Japan"

skd=5 #2 #1 #4 #3 #2 #slopes average factor
#Données de processus
t_cases = 0
t_recover = 0
t_deaths = 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]
            if day < 4+skd:
                day_confirmed_r[day-4] += data_r.iloc[i][day]
            else:
                day_confirmed_r[day-4] += (data_r.iloc[i][day] - data_r.iloc[i][day-skd])/(skd)
        t_recover += data_r.iloc[i][day]        
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]        
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] -confirmed_d[day-4]
            diff_confirmed[day - 4] += confirmed[day-4] /  (confirmed_r[day - 4]+confirmed_d[day-4])
            if day == 4:
                day_confirmed[day-4] += data.iloc[i][day]
            else:
                day_confirmed[day-4] += data.iloc[i][day] - data.iloc[i][day-1]

L'ajustement de la courbe de guérison est effectué ci-dessous.

t_max=len(confirmed)
dt=1
t=np.arange(0,t_max,dt)
t1=t

obs_i = confirmed_r[:]
#function which estimate i from seir model func 
def estimate_i(ini_state,r0,alpha):
    est = r0*np.exp(alpha*(t))
    return est

def y(params):
    est_i=estimate_i(ini_state,params[0],params[1])
    return np.sum((est_i-obs_i)*(est_i-obs_i))
r0=1
alpha = 1
ini_state=[5.70579672, 0.00755685]
#optimize logscale likelihood function
mnmz=minimize(y,ini_state,method="nelder-mead")
print(mnmz)
r0,alpha = mnmz.x[0],mnmz.x[1]
est=estimate_i(ini_state,r0,alpha)

Les quantités $ \ gamma $, $ R $, $ C = \ gamma * (R-1) $ sont calculées ci-dessous.

diff_est=[0] * (len(data.columns) - 4)
gamma_est=[0] * (len(data.columns) - 4)
R_est = [0] * (len(data_d.columns) - 4)
R_0 = [0] * (len(data_d.columns) - 4)
C = [0] * (len(data_d.columns) - 4)
for i in range(1,t_max):
    diff_est[i]=est[i]-est[i-1]
for i in range(0, len(confirmed), 1):        
    if confirmed[i] > 0:    
        gamma_est[i]=diff_est[i]/confirmed[i]
        R_est[i]= 1+day_confirmed[i]/diff_est[i] # diff_est=gamma*confirmed
        C[i]=gamma_est[i]*(R_est[i]-1)
    else:
        continue

Après cela, le résultat du calcul ci-dessus est dessiné.

Résumé

・ Γ et R ont été obtenus à partir de données réelles ・ On peut dire que la situation au Japon vient d'entrer dans une zone dangereuse.

・ Je ne cherche pas β, alors je vais l'essayer. ・ Je voudrais comparer chaque pays avec les différents paramètres obtenus, y compris ceux obtenus par l'équation différentielle.