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-3rd-wave-five.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 = 14000000 # 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; E0 = 0; R0 = 0.0 # SIR differential equation# S = SEIR[0], E = SEIR[1], I = SEIR[2], R = SEIR[3]def SEIReq(SEIR, t, beta, ipusilon, gamma): dSdt = -beta*SEIR[0]*SEIR[2] dEdt = beta*SEIR[0]*SEIR[2] - ipusilon*SEIR[1] dIdt = ipusilon*SEIR[1] - gamma*SEIR[2] dRdt = gamma*SEIR[2] return [dSdt, dEdt, dIdt, dRdt] def I(t, beta, ipusilon, gamma): SEIRlist = odeint(SEIReq, (S0, E0, I0, R0), t, args = (beta, ipusilon, gamma)) return SEIRlist[:,2] p0 =[1, 1, 1]optparams, cov = optimize.curve_fit(I, t, normalized_cases,p0)print('R0=',optparams[0]/optparams[2])print(optparams[0])print(optparams[1])print(optparams[2])fitted = I(t, *optparams) plt.scatter(t, cases)plt.plot(t, fitted*Tokyo)plt.xlabel('the number of days from 2020/12/1')plt.ylabel('the number of confirmed cases in Tokyo')plt.show()f.close() # close the csv-file
0 コメント