[PYTHON] [Einführung in das SEIR-Modell] Versuchen Sie, COVID-19-Daten anzupassen ♬

Ich habe versucht, die Infektionsdaten der neuen Korona mit dem SEIR-Modell anzupassen. Die Methode ist dieselbe wie bei der vorherigen [Einführung zur Minimierung] SEIR-Modelldatenanalyse (. Die Infektionsdaten wurden von der folgenden Referenzstelle erhalten. Die folgenden Ergebnisse wurden in [Einführung in matplotlib] Lesen der Endzeit aus COVID-19-Daten ( für die Daten von Wuhan erhalten. Der Anfangswert wird ausgewählt und angepasst, so dass eine solche Kurve erhalten werden kann, und das Ergebnis wird wie folgt erhalten. Daher möchte ich dieses Mal für jedes Land ähnliche Anpassungen vornehmen und prüfen, ob neue Erkenntnisse gewonnen werden können. Die Daten stammen aus dem Folgenden, und sofern nicht anders angegeben, handelt es sich um ** Daten zum 24. März 2020 **, die oben genannten Daten ** zum 22. März 2020 **. Beachten Sie außerdem, wie in der Code-Erläuterung erläutert, dass diesmal die Anpassung für die Infektionsdaten durchgeführt wird und die Anpassung, die auch die Heilungsdaten verwendet, von verschiedenen Programmen ausgeführt wird. 【Referenz】 CSSEGISandData/COVID-19

Was ich getan habe

・ Das Ende Japans ・ Codeerklärung ・ Anpassung nationaler Daten

・ Das Ende Japans

Ich habe versucht, es wie oben erwähnt sanft zu schreiben, aber da die Stadtregierung von Tokio angekündigt hat, dass die Stadt geschlossen wird, möchte ich erwähnen, ob das Ende Japans zuerst gesehen werden kann. Aus der Schlussfolgerung geht hervor, dass das Ende sicherlich in einem weiteren Monat zu sehen sein wird. Tatsächlich gab es in den obigen Daten vom 22. März eine Situation, in der ich eine Weile eingeschlafen bin und dies der Höhepunkt der Infektion war, aber da die Anzahl der Infektionen hier linear anstieg, wurde die Sättigung leicht verschoben. Als Anpassung habe ich die folgende Grafik für die Daten bis zum 24. März erhalten.

Die Anpassung bestimmt gleichzeitig die Anzahl der Infektionen und die Heilungskurve. Die Daten für Japan sind vom 53. bis zum 59. etwas seltsam, bewegen sich aber im Allgemeinen natürlich. Also habe ich eine schöne Kurve wie folgt. Mit anderen Worten, das Ergebnis ist, dass ein Peak in ungefähr 10 Tagen auftritt, eine Heilungskurve in ungefähr einem Monat einen Peak aufweist und der Rest sanft beendet wird.

Die erhaltenen Parameter sind

beta_const=0.159
lp=1.97e-23
gamma=52.7
N=1474
R0 = 4.45

・ Codeerklärung

Der gesamte Code ist unten. Die beiden Gründe sind, dass sich die bereitgestellte Datenstruktur geändert hat und die aktuellen Daten wie in (2) einzeln gelesen werden müssen. Außerdem wurde in (2) die Anpassungsfunktion von der normalsten wahrscheinlichsten Funktion in den normalen quadratischen Fehler geändert, damit sie leichter in die Kurve passt. ①collective_particles/fitting_COVID.py ②collective_particles/fitting_COVID2.py

Der einfache Code wird unten erklärt. Verwenden Sie die folgende Bibliothek. Dies ist die Summe der Artikel aus den beiden vorherigen Artikeln.

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

Das Folgende entspricht dem Zeichnen infizierter Daten.

#Lesen Sie CSV-Daten mit Pandas.
data = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv')
data_r = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Recovered.csv')
data_d = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Deaths.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)

#City
city = "Hubei" #Wuhan Daten

#Prozessdaten
t_cases = 0

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] -  data_r.iloc[i][day] 
            confirmed_r[day - 4] += data_r.iloc[i][day]
            confirmed_d[day - 4] += data_d.iloc[i][day]

Verwenden Sie das SEIR-Modell.

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

Hier ändert sich das erhaltene Bild in Abhängigkeit vom Wert von N0. Es ist gut, das zu passen, aber diesmal ist es nicht fertig. Bei der Anpassung von Hubei habe ich auf die folgenden Werte Bezug genommen. Für ** N und S0 wurde jedoch die tatsächliche Anzahl infizierter Personen (Werte verschoben) angenommen, basierend auf der Idee, dass die Bewohner nicht alle Teilnehmer des Spiels sind, sondern das Ziel der Infektion **. 【Referenz】 ・ Ich habe versucht, das Verhalten des neuen Koronavirus mit dem SEIR-Modell vorherzusagen. Anfangswerteinstellung

#solve seir model
N0 = 70939 #Hubei
N,S0,E0,I0=int(N0*1.),N0,int(318),int(389) #N is total population, S0 initial value of fresh peaple
ini_state=[S0,E0,I0,0] #initial value of diff. eq.
beta,lp,ip=1, 2, 7.4 
t_max=60 #len(days_from_22_Jan_20)
dt=0.01
t=np.arange(0,t_max,dt)
obs_i = confirmed

Eine Funktion, die einen Auswertungswert für einen Eingabeparameter zurückgibt. Das folgende v [2] gibt die Anzahl der Infektionen zurück.

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

Die Zielfunktion zur Minimierung unten (im Folgenden wird die wahrscheinlichste Funktion verwendet).

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

Finden Sie den Minimalwert mit minim @ scipy und berechnen Sie die Basisreproduktionsnummer.

#optimize logscale likelihood function
mnmz=minimize(y,[beta,lp,ip,N],method="nelder-mead")
print(mnmz)
#R0
beta_const,lp,gamma_const,N = mnmz.x[0],mnmz.x[1],mnmz.x[2],mnmz.x[3] #Infektionsrate, Wartezeit der Infektion, Entfernungsrate (Wiederherstellungsrate)
print(beta_const,lp,gamma_const,N)
R0 = N*beta_const*(1/gamma_const)
print(R0)

Grafische Darstellung der Berechnungsergebnisse.

#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(confirmed_r,"*", color="green",label = "data_r")
lns3=ax1.plot(estimate_i(ini_state,mnmz.x[0],mnmz.x[1],mnmz.x[2],mnmz.x[3]), label = "estimation")
lns0=ax1.plot(t,odeint(seir_eq,ini_state,t,args=(mnmz.x[0],mnmz.x[1],mnmz.x[2],mnmz.x[3])))
lns_ax1 = lns1+lns2 + lns3 + lns0
labs_ax1 = [l.get_label() for l in lns_ax1]
ax1.legend(lns_ax1, labs_ax1, loc=0)
ax1.set_ylim(0,N0)

Im Folgenden wird ein Diagramm mit großer Spannweite angezeigt, damit das Ende sichtbar ist.

• Bis zu 200 Tage hier angezeigt Infolgedessen kann das zu Beginn gezeigte ähnliche Ergebnis wie das von Wuhan erhalten werden. Das Bild von Wuhan ist jedoch ein Bild bis zum 60 ..

t_max=200 #len(days_from_22_Jan_20)
t=np.arange(0,t_max,dt)

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

Andererseits wurden bei der gleichzeitigen Anpassung der Anzahl von Infektionen und der Anzahl von Heilungen die Bewertungsfunktion und die Zielfunktion in der Mitte wie folgt geändert.

#function which estimate i from seir model func 
def estimate_i(ini_state,beta,lp,ip,N):
    v=odeint(seir_eq,ini_state,t,args=(beta,lp,ip,N))
    est=v[0:int(t_max/dt):int(1/dt)] 
    return est[:,2],est[:,3] #v0-S,v1-E,v2-I,v3-R
    
#define logscale likelihood function
def y(params):
    est_i_2,est_i_3=estimate_i(ini_state,params[0],params[1],params[2],params[3])
    return np.sum((est_i_2-obs_i_2)*(est_i_2-obs_i_2)+(est_i_3-obs_i_3)*(est_i_3-obs_i_3))

In diesem Fall ist es möglich, nur die Anzahl der Infektionen oder nur die Anzahl der Heilungen anzupassen.

・ Anpassung nationaler Daten

· Iran Der Iran, der eine relativ hohe Heilungsrate aufweist, sättigt schneller als Japan und wird wahrscheinlich enden, aber die Anzahl der Infektionen steigt immer noch und die Heilungsrate ist träge, so dass es möglicherweise länger dauert.

·Korea Irgendwie wurde die Verzerrung im Vergleich zu den vorherigen Daten deutlicher und die Glaubwürdigkeit der Daten wurde gebrochen. Dennoch steigt die Heilungsrate insgesamt an und das Ende scheint relativ früh zu sein.

·Italien Die Sterblichkeitsrate ist extrem hoch und gefährlich, aber die Daten sind sehr sauber und anspruchsvoll. Daher war es einfach zu montieren und ich kam sauber voran. Die Heilungsrate beträgt immer noch etwa 10%, und es scheint, dass es einige Zeit dauern wird (der Heilungspeak beträgt 2 Monate), aber wenn der Peak danach gesehen wird, ist das Ende nahe.

Das Folgende ist die anfängliche Expansionsperiode, und es scheint, dass die Anpassungsgenauigkeit schlecht und nicht zuverlässig ist, aber ich werde sie veröffentlichen. ·Spanien Es steht immer mehr auf, aber es ist schwer zu erreichen. ·Deutschland Dies ist noch nicht abgeschlossen, und die Anzahl der Heilungen ist in erster Linie zu gering, so dass es den Anschein hat, dass es sich in der anfänglichen Expansionsphase befindet.

Zusammenfassung

・ Ich habe versucht, COVID-19-Daten anzupassen ・ Japan, Iran und Südkorea werden nach China voraussichtlich ein Ende haben. ・ Andere Europa und die Vereinigten Staaten befinden sich noch in der Expansionsphase, und es ist schwierig, das Ende vorherzusagen.

・ Ich habe diesmal ungefähr gesprochen, möchte aber eine detailliertere Analyse durchführen, während ich mir die täglichen Updates ansehe.

Bonus

Die gleichzeitige Anpassung von Wuhan ist wie folgt Wuhans Daten scheinen durch dieses Modell nicht erklärt zu werden (> <) Originale Daten Nur Infektionsdaten anpassen Heilungsdaten passen nur Anpassen von Infektions- und Heilungsdaten