SIRモデルを扱った新型コロナウイルスの新規感染者シミュレーション

Python

import numpy as np import matplotlib.pyplot as plt from scipy.integrate import odeint from scipy import optimize import csv # loading csv-filef = open("./COVID-first_wave.csv", "r", encoding="UTF-8", errors="", newline="" )fcsv = csv.reader(f, delimiter=",", doublequote=True, lineterminator="\r\n", quotechar='"', skipinitialspace=True) next(f) # skip to the header of the csv-file cases = [] for row in fcsv: cases.append(int(row[1])) Tokyo = 13999568 # the population of Tokyo in 2020normalized_cases = np.array(cases, dtype = float)/Tokyo days = len(cases)t = np.arange(days)# initial valuesI0 = normalized_cases[0]; S0 = 1.0 - I0; R0 = 0.0 # SIR differential equation# S = SIR[0], I = SIR[1], R = SIR[2]def SIReq(SIR, t, beta, gamma): dSdt = -beta*SIR[0]*SIR[1] dIdt = beta*SIR[0]*SIR[1] - gamma*SIR[1] dRdt = gamma*SIR[1] return [dSdt, dIdt, dRdt] def I(t, beta, gamma): SEIRlist = odeint(SIReq, (S0, I0, R0), t, args = (beta, gamma)) return SEIRlist[:,1] optparams, cov = optimize.curve_fit(I, t, normalized_cases)print('R0=',optparams[0]/optparams[1])fitted = I(t, *optparams) plt.scatter(t, cases)plt.plot(t, fitted*Tokyo)plt.xlabel('the number of days from 2020/2/13')plt.ylabel('the number of confirmed cases in Tokyo')plt.show()f.close() # close the csv-file