[PYTHON] [SIR-Modellanalyse] Transformieren Sie die Formel, um γ und die effektive Reproduktionszahl R ♬ zu bestimmen

Wenn die neuesten Daten natürlich extrapoliert werden, wird die folgende Abbildung erhalten. Ich denke, Leute, die die Grafik lesen können, werden es verstehen, aber mir wurde klar, dass die Situation ziemlich schlecht war.

Daher möchte ich über Folgendes sprechen, ohne dieses Thema anzusprechen. Die Daten, die dieses Mal verarbeitet werden, sind die Daten bis zum 04.03.2020, die auf der Referenzseite veröffentlicht wurden. Für die Anzahl der Infektionen in Japan wurde der aus Referenz (2) abgelesene Wert verwendet. 【Referenz】 ①Novel Coronavirus (COVID-19) Cases, provided by JHU CSSE ② Neue Corona Virus @ Nittele-NACHRICHTEN aus Daten und Grafiken

Was ich getan habe

・ Koeffiziententransformation des SIR-Modells ・ Bedeutung des SIR-Modells ・ Erklärung des Ergebnisses ・ Codeerklärung

・ Koeffiziententransformation des SIR-Modells

Dieses Mal werde ich versuchen, die tägliche Änderung von $ \ gamma $ und die Anzahl der effektiven Reproduktionen zu ermitteln, indem ich das folgende SIR-Modell wie folgt transformiere.

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

Erstens kann die dritte Gleichung wie folgt transformiert werden, und da $ I $ und $ r $ Beobachtungswerte erhalten, kann $ \ gamma $ erhalten werden.

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

Als nächstes, wenn Sie dieses $ \ gamma $ in die zweite Gleichung einsetzen

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

Daher wird R erhalten.

・ Bedeutung des SIR-Modells

Unter Verwendung dieser beiden Größen kann die obige Differentialgleichung wie folgt umgeschrieben werden.

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

Hier ist $ \ gamma $ ein Faktor, der sich auf (Wiederherstellungsrate, Isolationsrate und Sterblichkeitsrate) bezieht und die Rate bedeutet, mit der eine infizierte Person nicht zur Infektion beiträgt (entfernt). Andererseits wird $ R $ als effektive Reproduktionsnummer bezeichnet. Dies ist eine Menge, die sich auf die Anzahl der von einer infizierten Person infizierten Personen bezieht. Daher zeigt die erste Gleichung, dass die Infektion im Verhältnis zur infizierten Person umso weiter voranschreitet, je größer $ R $ und $ \ gamma $ sind. Dies ist der Grund, warum die Ausbreitung der Infektion in sogenannten Mausberechnungen (exponentiell) zunimmt. Andererseits bedeutet die dritte Formel, dass je größer $ \ gamma $ ist, desto mehr geheilt wird im Verhältnis zur infizierten Person. Die zweite Gleichung ist proportional zu $ (R-1) $, und es wird gezeigt, dass der positive oder negative Wert dieses Wertes bestimmt, ob die Anzahl infizierter Personen zunimmt oder abnimmt. Darüber hinaus ist ersichtlich, dass die Änderung exponentiell ist, da sie proportional zur Anzahl der infizierten Personen ist. Und diese Änderung ist auch proportional zu $ \ gamma $.

・ Erklärung des Ergebnisses

Das erste, was diesmal auffällt, ist, dass sich die Anzahl der Infektionen (rot) seit dem 64. Tag (3. bis 9. April) erheblich geändert hat. Und der Hauptzweck dieser Zeit ist es, nach $ \ gamma $ und $ R $ für diesen Trend in Japan zu fragen. Vor allem aber schwanken die täglichen Daten stark, so dass das Problem darin bestand, wie man damit umgeht. Wenn Sie sich die Daten ansehen, können Sie sehen, dass die Kurve der Anzahl der Heilungen schrecklich ist. Ich habe das als sogenannte Glättung betrachtet (3 Punkte, 5 Punkte, 7 Punkte ...), aber es ist schlimmer als das, und insgesamt ändert es sich fast linear, also passe ich es mit einer Exponentialfunktion an. Ich habe mich entschieden. In ähnlicher Weise wurde auch die Anzahl der Infektionen angepasst, dies wurde jedoch nur für den geraden Teil durchgeführt. Und dieses Mal wurde es nur für die eingangs dargestellte Extrapolation verwendet, nicht für diese Koeffizientenberechnung. Die Kurve $ I / (R + D) $ zeigt das Verhältnis von infizierten und nicht beitragenden Personen. Diese Zahl ist eine Zahl, die die Heilungsrate erhöht, und wenn diese Zahl abnimmt und unter 1 fällt, ist das Ende endlich sichtbar. In Japan wurde erwartet, dass die Zahl einmal 3 oder weniger betragen und dann so wie sie ist auf 1 übergehen würde, aber nach dem 64. Tag oben begann diese Zahl ebenfalls signifikant zu steigen. Außerdem war es bis dahin schwach zurückgegangen, aber es begann linear zu steigen, was wirklich überraschend ist. Das untere Diagramm zeigt $ \ gamma $ und $ R $, die ich dieses Mal finden wollte, die Anzahl der täglichen Heilungen (5-Tage-Durchschnitt), die für die Berechnung verwendet wurden, und ihre ungefähren Kurven. $ \ Gamma $ dividiert die tägliche Heilungszahl (Differenzierung der Heilungskurve) durch die Anzahl der Infektionen. Da die Heilungszahlkurve jedoch stark schwankt, wurde hier die ungefähre Kurve verwendet. Trotzdem gab es aufgrund der Schwankung der Infektionszahlkurve erhebliche Schwankungen. Es wurden jedoch jeden Tag Veränderungen beobachtet, und obwohl diese Zahl bis zum 64. Tag ebenfalls stetig anstieg, fiel sie von dort linear ab. Andererseits ist $ R $ umgekehrt proportional zu diesem $ \ gamma $, und das Ergebnis funktioniert auf diese Weise. Der niedrigste Wert liegt jedoch bei etwa 3 effektiven Reproduktionen, und die Position ist leicht nach links verschoben. In jedem Fall hat sich dieser Bereich vollständig vergrößert, und es ist ersichtlich, dass die Anzahl der effektiven Reproduktionen auf etwa 10 ansteigt. Die neuesten Werte sind wie folgt. Die Anzahl der effektiven Reproduktionen ist ein Vergleich mit 1, was ein sehr großer Wert ist (ein Wert, bei dem sich die Infektion schnell ausbreitet).

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

Der tatsächliche Beitrag hängt von $ \ gamma (R-1) $ ab, und die obige Grafik zeigt, dass der Beitrag der Abnahme von $ \ gamma $ und der Beitrag der Zunahme von $ R-1 $ ausgeglichen sind. Wenn dieser Betrag grafisch dargestellt wird, wird er daher wie folgt. Infolgedessen war der Beitrag von $ R $ groß, und sobald er auf ungefähr 0,05 abnahm, stieg er nach dem 64. Tag an, und jetzt ist zu sehen, dass er auf ungefähr 0,14 gestiegen ist. Man kann sagen, dass sich die Infektionswahrscheinlichkeit in den letzten zwei Wochen verdreifacht hat.

・ Codeerklärung Ich habe den gesamten Code unten angegeben collective_particles/fitting_gamma_R.py Dieses Mal werde ich nur die Hauptteile erklären. Erstens sind das Lesen und Verarbeiten der Daten wie folgt.

city = "Japan"

skd=5 #2 #1 #4 #3 #2 #slopes average factor
#Prozessdaten
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]

Die Anpassung der Heilungskurve wird unten durchgeführt.

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)

Die Mengen $ \ gamma $, $ R $, $ C = \ gamma * (R-1) $ werden unten berechnet.

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

Danach wird das obige Berechnungsergebnis gezeichnet.

Zusammenfassung

・ Γ und R wurden aus realen Daten erhalten ・ Man kann sagen, dass die Situation in Japan gerade in eine gefährliche Zone eingetreten ist.

・ Ich suche nicht nach β, also werde ich es versuchen. ・ Ich werde jedes Land mit den verschiedenen erhaltenen Parametern vergleichen, einschließlich der durch die Differentialgleichung erhaltenen.