[PYTHON] [Einführung in das SIR-Modell] Prognostizieren Sie die Endzeit jedes Landes mit der COVID-19-Datenanpassung ♬

Im vorherigen SEIR-Modell wurden unter den erhaltenen Parametern Werte wie $ lp \ fallenddotseq 1e ^ {-23} $ erhalten. In einem solchen Fall besteht die Möglichkeit, dass der Ziffernverlust oder die Differentialgleichung überhaupt nicht funktioniert. Seien Sie also vorsichtig. Also werde ich dieses Mal dieses Teil ändern und versuchen, mit dem einfachsten Modell zu passen.

Was ich getan habe

・ Eine kleine Theorie ・ Codeerklärung ・ Das Ende Japans ・ Ende jedes Landes ・ Eine kleine Diskussion

・ Eine kleine Theorie

Das SIR-Modell ist die folgende einfache Zeitentwicklungsgleichung.

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

Die erste Gleichung lautet, dass die Abnahme nicht infizierter Personen proportional zur Anzahl infizierter Personen und zur Anzahl nicht infizierter Personen und umgekehrt proportional zur Gesamtzahl der Probanden ist. Der Koeffizient β ist die Infektionsrate und die Anzahl nicht infizierter Personen ist das Verhältnis. Verringern. Als β / N kann es als Pro-Kopf-Infektionsrate angesehen werden. Die zweite Gleichung ist die Formel zur Erhöhung der Anzahl infizierter Personen. Der Anstieg ist der Unterschied zwischen dem Zustrom von nicht infizierten und den geheilten Personen. Die dritte Gleichung ist eine Formel, nach der ein bestimmter Prozentsatz infizierter Personen geheilt wird, und dieser Koeffizient γ ist eine Heilungsrate, die als Entfernungsrate bezeichnet wird. 【Referenz】 ・ [Verbreitet oder konvergiert diese Infektion: die physikalische Bedeutung und Bestimmung der Anzahl der Reproduktionen R](https://rad-it21.com/%E3%82%B5%E3%82%A4%E3%82%A8] % E3% 83% B3% E3% 82% B9 / kado-shinichiro_20200327 /)

Hier kann die zweite Gleichung wie folgt transformiert werden.

\frac{dI}{dt} =  \gamma (\frac{\beta}{\gamma} \frac{S}{N} - 1)I\\

Das heißt, wenn $ S / N $ als konstant angesehen wird, kann die folgende Lösung formal erhalten werden.

• Es kann davon ausgegangen werden, dass die zeitliche Entwicklung innerhalb einer ausreichend kurzen Zeit liegt.

I =  I_0\exp(\gamma (\frac{\beta}{\gamma} \frac{S}{N} - 1)t)\\

Zu Beginn der Infektion ist es $ S_0 \ fallenddotseq N $. In diesem Fall kann die folgende Lösung erhalten werden.

{\begin{align}
I &=  I_0\exp(\gamma (R_0 - 1)t)\\
Hier\\ 
R_0 &= \frac{\beta}{\gamma}
\end{align}
} 

Hier nimmt in der Formel von $ I $ die Anzahl der infizierten Personen exponentiell ab, wenn $ R_0 <1 $ ist, und steigt umgekehrt exponentiell an, wenn $ R_0> 1 $. Daher wird dieses ** $ R_0 $ als Basisreproduktionsnummer ** bezeichnet und dient als Leitfaden für die Bestimmung, ob eine Infektionsübertragung stattfinden wird oder nicht. Darüber hinaus wird in der obigen Referenz

R = R_0\frac{S}{N}

Ist definiert als ** die Anzahl der effektiven Reproduktionen **. Dieses $ R $ wird kleiner, wenn $ \ frac {S} {N} $ mit fortschreitender Infektion kleiner wird. Das heißt, wie Sie aus der obigen formalen Lösung sehen können, wenn $ R_0> 1 $, $ R> 1 $ bereits in der Anfangsphase der Expansion vorhanden ist und sich die Infektion ausbreitet, aber schließlich ist $ \ frac {S} {N} $ klein. Daher ist es ein Indikator dafür, dass das Ende beginnt, wenn $ R <1 $ realisiert wird. Wenn andererseits die linke Seite von $ \ frac {dI} {dt} = 0 $ = 0 ist, erreicht die Infektion ihren Höhepunkt und die Anzahl der nicht infizierten Personen zu diesem Zeitpunkt beträgt $ S_p $.

{\begin{align}
R_0 &= \frac{\beta}{\gamma}\\
& = \frac{N}{S_p}　\\
\end{align}
}

Ist zu diesem Zeitpunkt eingerichtet

{\begin{align}
R &= R_0\frac{S_p}{N}\\
&= 1
\end{align}
}

Es stellt sich heraus, dass dies zutrifft. Das heißt, ** die Anzahl der wirksamen Reproduktionen beträgt am Höhepunkt der Infektion R = 1, und die Übertragung der Infektion endet danach. ** ** ** Dieses Mal werde ich die Endvorhersage betrachten, indem ich auch diese ** effektive Reproduktionszahl ** aufzeichne.

・ Codeerklärung

Der gesamte Code ist unten. ・ Collective_particles /itting_SIR_COVID2.py Die verwendete Bibliothek lautet wie folgt wie beim letzten Mal.

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

Daten im Voraus von COVID-19-Site Wird heruntergeladen.

#Lesen Sie CSV-Daten mit Pandas.
data = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv')
data_r = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv')
data_d = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.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)

Die Aufzeichnungsreihenfolge der Daten wurde beim letzten Mal korrekt angeordnet, aber diesmal ist die Reihenfolge der Länder und Regionen der drei Daten nicht dieselbe. Holen Sie sich also wie folgt und berechnen Sie die Anzahl der vorhandenen Infektionen für die Infektionsnummer-Daten Um dies zu tun, versuche ich es zu bekommen, nachdem ich die Heilungszahldaten erhalten habe. Darüber hinaus, wenn Aussagen in jedem Land und jeder Region synonym verwendet werden. Außerdem werden die Gesamtzahl der Heilungen und die Gesamtzahl der Todesfälle, die für die Sterblichkeitsrate und die Heilungsrate erforderlich sind, dieses Mal nicht berechnet, sodass sie auskommentiert werden.

#City
city = "Japan"
#Prozessdaten
t_cases = 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]
            #t_recover += data_r.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]            
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]

Dieses Mal werden wir das SIR-Modell verwenden. Die Bewertungsfunktion gibt auch die nicht infizierte Nummer S, die infizierte Nummer I und die geheilte Nummer R zurück.

#define differencial equation of sir model
def sir_eq(v,t,beta,gamma):
    a = -beta*v[0]*v[1]
    b = beta*v[0]*v[1] - gamma*v[1]
    c = gamma*v[1]
    return [a,b,c]
def estimate_i(ini_state,beta,gamma):
    v=odeint(sir_eq,ini_state,t,args=(beta,gamma))
    est=v[0:int(t_max/dt):int(1/dt)] 
    return est[:,0],est[:,1],est[:,2] #v0-S,v1-I,v2-R

Die Zielfunktion ist eine Normalverteilung. Hier wird die lineare Summe verwendet, damit die Anzahl der Infektionen I und die Anzahl der Heilungen R gleichzeitig angepasst werden können.

#define logscale likelihood function
def y(params):
    est_i_0,est_i_1,est_i_2=estimate_i(ini_state,params[0],params[1])
    return np.sum(f1*(est_i_1-obs_i_1)*(est_i_1-obs_i_1)+f2*(est_i_2-obs_i_2)*(est_i_2-obs_i_2))

Geben Sie die Gesamtzahl der Infektionen in dem oben in N0 ausgewählten Land (oder der Region) ein. Außerdem bestimmen f1 und f2, auf welcher Kurve die Anpassung durchgeführt werden soll. N, S0, I0, R0 sind die Anfangswerte von jedem, aber dies sind auch Konstanten, die entsprechend den Anpassungsbedingungen angepasst werden sollten. Die Anfangswerte von Beta und Gamma sind in den meisten Ländern ähnlich, aber es scheint besser, die einmal erhaltenen Werte erneut einzugeben und rekursiv zu berechnen.

• Sie sollten auch diese Anfangswertanpassung durchführen, diesmal jedoch nicht

#solve seir model
N0 = 1517
f1,f2 = 1,1 #data & data_r fitting factor
N,S0,I0,R0=int(N0*1.5),int(N0*1.5),int(10),int(0)
ini_state=[S0,I0,R0]
beta,gamma = 1.6e-6, 1e-3 #Der Anfangswert konvergiert nicht, es sei denn, er liegt in gewissem Maße nahe

t_max=len(days_from_22_Jan_20) #Anpassung innerhalb eines bestimmten Datenbereichs
dt=0.01
t=np.arange(0,t_max,dt)

Ich werde versuchen, es unten anzuzeigen, einschließlich der Überprüfung des Anfangswertes.

plt.plot(t,odeint(sir_eq,ini_state,t,args=(beta,gamma)))
plt.legend(['Susceptible','Infected','Recovered'])

obs_i_1 = confirmed
obs_i_2 = confirmed_r

plt.plot(obs_i_1,"o", color="red",label = "data_1")
plt.plot(obs_i_2,"o", color="green",label = "data_2")
plt.legend()
plt.pause(1)
plt.close()

Machen Sie eine Anpassung. Auch diesmal habe ich minim @ scipy verwendet. Berechnen Sie die Grundreproduktionsnummer R0. Da sich der Variablenname mit der oben angegebenen anfänglichen Heilungsnummer überschneidet, wird er in Kleinbuchstaben geschrieben.

#optimize logscale likelihood function
mnmz=minimize(y,[beta,gamma],method="nelder-mead")
print(mnmz)
#R0
#N_total = S_0+I_0+R_0
#R0 = N_total*beta_const *(1/gamma_const)
beta_const,gamma = mnmz.x[0],mnmz.x[1] #Infektionsrate, Entfernungsrate (Heilungsrate)
print(beta_const,gamma)
r0 = N*beta_const*(1/gamma)
print(r0)

Stellen Sie es unten grafisch dar. Die Unterachse von ax3 und ax4 wird hinzugefügt, um die effektive Reproduktionsnummer zu zeichnen.

#plot reult with observed data
fig, (ax1,ax2) = plt.subplots(2,1,figsize=(1.6180 * 4, 4*2))
ax3 = ax1.twinx()
ax4 = ax2.twinx()

lns1=ax1.plot(confirmed,"o", color="red",label = "data_I")
lns2=ax1.plot(confirmed_r,"o", color="green",label = "data_R")
est_i_0,est_i_1,est_i_2=estimate_i(ini_state,mnmz.x[0],mnmz.x[1])
lns3=ax1.plot(est_i_1, label = "estimation_I")
lns0=ax1.plot(est_i_2, label = "estimation_R")
lns_1=ax3.plot((S0-est_i_1-est_i_2)*r0/N,".", color="black", label = "effective_R0")
ax3.set_ylim(0,r0+1)

lns_ax1 = lns1+lns2 + lns3 + lns0 +lns_1
labs_ax1 = [l.get_label() for l in lns_ax1]
ax1.legend(lns_ax1, labs_ax1, loc=0)
ax1.set_ylim(0,N0)
ax1.set_title('SIR_{} f1_{:,d} f2_{:,d};N={:.0f} S0={:.0f} I0={:.0f} R0={:.0f} R={:.2f}'.format(city,f1,f2,N,S0,I0,R0,(S0-est_i_1[-1]-est_i_2[-1])*r0/N))
ax1.set_ylabel("Susceptible, Infected, recovered ")
ax3.set_ylabel("effective_R0")

Die zukünftige Prognose ist unten grafisch dargestellt.

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

lns4=ax2.plot(confirmed,"o", color="red",label = "data")
lns5=ax2.plot(confirmed_r,"o", color="green",label = "data_r")
est_i_0,est_i_1,est_i_2=estimate_i(ini_state,mnmz.x[0],mnmz.x[1])
lns6=ax2.plot(est_i_0, label = "estimation_S")
lns7=ax2.plot(est_i_1, label = "estimation_I")
lns8=ax2.plot(est_i_2, label = "estimation_R")
lns_6=ax4.plot((S0-est_i_1-est_i_2)*r0/N,".", color="black", label = "effective_R0")

lns_ax2 = lns4+lns5 + lns6 + lns7+ lns8 +lns_6
labs_ax2 = [l.get_label() for l in lns_ax2]
ax2.legend(lns_ax2, labs_ax2, loc=0)
ax4.set_ylim(0,r0+1)
ax2.set_ylim(0,S0)
ax2.set_title('SIR_{};b={:.2e} g={:.2e} r0={:.2f}'.format(city,beta_const,gamma,r0))
ax2.set_ylabel("Susceptible, Infected, recovered ")
ax4.set_ylabel("effective_R0")

ax1.grid()
ax2.grid()
plt.savefig('./fig/SIR_{}f1_{:,d}f2_{:,d};b_{:.2e}g_{:.2e}r0_{:.2f}N_{:.0f}S0_{:.0f}I0_{:.0f}R0_{:.0f}.png'.format(city,f1,f2,beta_const,gamma,r0,N,S0,I0,R0)) 
plt.show()
plt.close()

・ Das Ende Japans

** Daten sind bis zum 27. März 2020 **. Die Daten sind wie folgt. Die Ergebnisse sind wie folgt. Dieses Ergebnis scheint vorerst passend zu sein, aber der heutige Bericht besagt, dass gestern am 28. März 200 Menschen infiziert wurden. Es liegt also außerhalb dieser Grafik und stimmt nicht überein. .. Ich habe dies als Bonus gepostet. Die erhaltenen Mengen sind wie folgt.

{\begin{align}
\beta &= 4.53e^{-5}\\
\gamma &= 1.58e^{-2}\\
N &= 2275\\
R0 &= 6.52\\
R(27.3.2020) &= 2.71\\
\end{align} 
}

・ Ende jedes Landes

Die Situation von Hubei

Die neuesten Daten (27.3.2020) lauten wie folgt. Die Daten von Wuhan können nicht wie folgt abgeglichen werden, selbst wenn Sie versuchen, die Daten sowohl der Anzahl der Infektionen als auch der Anzahl der Heilungen anzupassen. Wie Sie jedoch aus der Grafik sehen können, kann gesagt werden, dass die effektive Reproduktionszahl R = 0,01, 0,08 ist.

{\begin{align}
\beta &= 3.19e^{-6}\\
\gamma &= 6.53e^{-2}\\
N &= 106462\\
R0 &= 5.20\\
R(27.3.2020) &= 0.01\\
\end{align} 
}

Das Folgende sind Armaturen von beiden gleichzeitig, und die maximale Anzahl von Infektionen hat sich signifikant verschoben.

{\begin{align}
\beta &= 2.34e^{-6}\\
\gamma &= 5.33e^{-2}\\
N &= 106462\\
R0 &= 4.67\\
R(27.3.2020) &= 0.08\\
\end{align} 
}

Status von Korea, Süd

Die neuesten Daten (27.3.2020) lauten wie folgt.

Die koreanischen Daten sind nicht so schlecht wie in Wuhan, aber sie sind fast gleich. Selbst wenn Sie versuchen, sowohl die Daten zur Infektionsnummer als auch zur Heilungsnummer anzupassen, können sie nicht wie folgt abgeglichen werden. Wie Sie in der Grafik sehen können, endet jedoch die effektive Reproduktionszahl R = 0,13.

{\begin{align}
\beta &= 2.25e^{-5}\\
\gamma &= 2.94e^{-2}\\
N &= 11365\\
R0 &= 8.68\\
R(27.3.2020) &= 0.13\\
\end{align} 
}

Wenn Sie versuchen, beide anzupassen, ist die Anpassung der Daten zur Infektionsnummer lose und der Infektionspeak verschiebt sich ebenfalls nach rechts. Es kann möglich sein, es anzupassen, wenn es besser passt, aber dieses Mal werde ich es nicht mehr weiterverfolgen.

{\begin{align}
\beta &= 2.11e^{-5}\\
\gamma &= 2.14e^{-2}\\
N &= 11365\\
R0 &= 11.17\\
R(27.3.2020) &= 0.18\\
\end{align} 
}

・ Situation des Iran

Das nächste Land, das wahrscheinlich enden wird, ist der Iran. Die neuesten Daten (27.3.2020) lauten wie folgt. Als ich hierher kam und die Heilungsrate 30% überstieg, schien es jedoch zu stagnieren. Trotzdem beträgt die effektive Anzahl von Reproduktionen 1,81 und wird bald 1 sein, so dass gesagt werden kann, dass die maximale Anzahl von Infektionen nahe ist.

{\begin{align}
\beta &= 2.92e^{-6}\\
\gamma &= 4.94e^{-2}\\
N &= 62478\\
R0 &= 3.70\\
R(27.3.2020) &= 1.81\\
\end{align} 
}

・ Status von Italien

Das nächste betroffene Land ist Italien. Die neuesten Daten (27.3.2020) lauten wie folgt. Die Daten für dieses Land sind sehr solide und es fühlt sich gut an, sowohl die Infektions- als auch die Heilungskurve genau anzupassen. Das Ergebnis ist wie folgt: Die effektive Reproduktionszahl beträgt 4,07, ist jedoch stark gesunken und scheint in etwa 10 Tagen zu enden.

{\begin{align}
\beta &= 1.46e^{-6}\\
\gamma &= 1.91e^{-2}\\
N &= 143448\\
R0 &= 10.99\\
R(27.3.2020) &= 4.07\\
\end{align} 
}

Wenn ich mir diese Daten ansehe, habe ich den Eindruck, dass die medizinische Versorgung nicht zusammengebrochen ist und dass das Management wie die Auswirkungen der städtischen Blockade solide ist. (Kommentare haben keine andere wissenschaftliche Grundlage als Grafiken. Es ist ein individueller Eindruck.)

・ Situation in Spanien

Als nächstes werde ich die Situation in Spanien zeigen, in der die Todesfälle rapide zunehmen. Wenn ich die Schlussfolgerung schreibe, scheint es, dass sie sich noch in der frühen Expansionsphase befindet und die Anpassung unbestimmt ist.

• Situationen, die mit diesem Modell derzeit nicht entschieden werden können Die neuesten Daten (27.3.2020) lauten wie folgt. Das Folgende zeigt den Fall, in dem N kleiner ist, und den Fall, in dem N größer ist. N usw. sind unsicher, aber die Anzahl der effektiven Reproduktionen beträgt 7,81 bzw. 8,09, aber die Endzeit (Infektionspeak) beträgt etwa 80-90 Tage (2 Wochen bis 1 Monat später), und wenn der Peakwert ignoriert wird, Man kann sagen, dass sie sich an einem ähnlichen Ort befinden.

{\begin{align}
\beta &= 1.49e^{-6}\\
\gamma &= 1.37e^{-2}\\
N &= 127542\\
R0 &= 13.85\\
R(27.3.2020) &= 7.81\\
\end{align} 
}

{\begin{align}
\beta &= 5.39e^{-7}\\
\gamma &= 1.94e^{-2}\\
N &= 354285\\
R0 &= 9.86\\
R(27.3.2020) &= 8.09\\
\end{align} 
}

Eigentlich gilt das Gleiche für Deutschland, Frankreich und die Schweiz, aber ich habe versucht, es anzupassen, aber aufgrund des schnellen Anstiegs diesmal (Parameter sind undefiniert), also werde ich es erneut schreiben.

・ Eine kleine Diskussion

Ich bin mir nicht sicher, ob ich es wirklich akzeptieren kann, also schreibe ich ein wenig. Der größte Nachteil dieser Armatur ist, dass sie an N befestigt ist. Es scheint, dass die Situation der spontanen Übertragung einer Infektion bei einem geschlossenen Subjekt richtig simuliert werden kann. Viele Länder der Welt haben jedoch zumindest bis zu diesem Zeitpunkt keine vollständigen Sperren vorgenommen, und es gab (in Prozent) viele infizierte Menschen selbst. In Japan nimmt die Zahl solcher ambulanten Patienten zu, und es gibt auch Situationen, in denen sich von dort aus Cluster bilden. Was ich hier geschrieben habe, ist eine Anpassung des bisherigen Übergangs, und ich kann die zukünftige Situation nicht vorhersagen. Damit diese Vorhersagen korrekt sind, wird davon ausgegangen, dass sich die gleiche Situation fortsetzt, und es ist erforderlich, den Anteil der ambulanten Patienten extrem kontrollieren zu können. Obwohl wir ein Modell schreiben können, das diese Effekte berücksichtigt, scheint es unmöglich, mit diesem Modell vorherzusagen, da es schwierig ist, die Anzahl der ambulanten Patienten vorherzusagen, die sich von Moment zu Moment zufällig ändern. Ich denke also, dass eine Anpassung möglich sein wird, wenn ganz Europa und die Vereinigten Staaten sichtbar werden und die erhaltenen Parameter diskutiert werden können.

Zusammenfassung

・ Ich habe versucht, COVID-19-Daten an das SIR-Modell anzupassen. ・ Es wurde festgestellt, dass das Ende von Wuhan nahe war und dass der Höhepunkt der Infektion in Südkorea vorbei war. ・ Es ist zu erwarten, dass der Iran und Italien in etwa 10 Tagen relativ schnell ihre höchste Infektion erreichen. ・ In Japan ist der Infektionspeak relativ spät, da er langsam gesättigt ist, und es ist zu erwarten, dass der Infektionspeak je nach ambulanten Patienten in 20 Tagen eintreten wird. ・ Wenn sich die aktuelle Situation weiterentwickelt, wird davon ausgegangen, dass im nächsten Monat in vielen Ländern der Höhepunkt der Infektion eintreten wird.

・ Ich möchte täglich überwachen und Änderungen in der Situation sehen ・ Ich möchte die erhaltenen Parameter auswerten

Bonus

Wenn ich mir die Daten anschaue, werden sie bis zum 28. März 2020 wiedergegeben **, also habe ich versucht, sie anzupassen. Die Daten von gestern sind jedoch Ausreißer, und die Spitzenzeit der Infektion ändert sich nicht wesentlich. Die erhaltenen Parameter sind wie folgt und sie haben sich im Vergleich zu den obigen nicht wesentlich geändert. ** Die Anzahl der effektiven Reproduktionen R ist zwar gering, nimmt jedoch zu. Seien Sie also vorsichtig **.

{\begin{align}
\beta &= 3.91e^{-5}\\
\gamma &= 1.65e^{-2}\\
N &= 2617\\
R0 &= 6.22\\
R(28.3.2020) &= 2.79\\
\end{align} 
}

Ich würde gerne die zukünftigen Trends sehen.