前提
Pythonで数理モデルを用いて新型コロナウイルスの感染者のデータをフィッテイングさせています。
その際に最小二乗法を用いています。
実現したいこと
SEIRモデルにおいて最小二乗法を用いて新規感染者Iのフィッティングをしたいと考えています。
発生している問題・エラーメッセージ
今回は何かエラーが出ているということではありません。気になる点がいくつかありましたのでそこを伺いたいです。 1. 今回のこのコードで新規感染者Iのフィッティングができているのか。 後述しますが特にこの部分が見よう見真似でやったため少し自信がありません。 自分でコードを書いたのですがここがどうなっているのかを詳細に教えていただけると助かります。 自分なりにodeint関数も少し勉強したのですがイマイチわかっていないと思われます。 def I(t, beta, ipusilon, gamma): SEIRlist = odeint(SEIReq, (S0, E0, I0, R0), t, args = (beta, ipusilon, gamma)) return SEIRlist[:,2] 2. (解決済み)今回フィッティングさせて際に t = 1 あたりで曲線が緩やかになっていないのが気になります。 今までSIRモデルでフィッティングさせた際に緩やかな曲線を描いていましたのでなぜこのようになるのか わかる方がいらっしゃったら教えていただきたいです。 尚、初期条件を色々いじっても曲線は折れました。
該当のソースコード
Python
1import numpy as np 2import matplotlib.pyplot as plt 3from scipy.integrate import odeint 4from scipy import optimize 5import csv 6 7# loading csv-file 8f = open("./COVID_second_wave.csv", "r", encoding="UTF-8", errors="", newline="" ) 9fcsv = csv.reader(f, delimiter=",", doublequote=True, lineterminator="\r\n", quotechar='"', skipinitialspace=True) 10 11next(f) # skip to the header of the csv-file 12 13cases = [] 14for row in fcsv: 15 cases.append(int(row[1])) 16 17Tokyo = 13999568 # the population of Tokyo in 2020 18normalized_cases = np.array(cases, dtype = float)/Tokyo 19days = len(cases) 20t = np.arange(days) 21# initial values 22I0 = normalized_cases[0]; S0 = 1.0 - I0; E0 = 0; R0 = 0.0 23 24# SIR differential equation 25# S = SEIR[0], E = SEIR[1], I = SEIR[2], R = SEIR[3] 26def SEIReq(SEIR, t, beta, ipusilon, gamma): 27 dSdt = -beta*SEIR[0]*SEIR[2] 28 dEdt = beta*SEIR[0]*SEIR[2] - ipusilon*SEIR[1] 29 dIdt = ipusilon*SEIR[1] - gamma*SEIR[2] 30 dRdt = gamma*SEIR[2] 31 32 return [dSdt, dEdt, dIdt, dRdt] 33 34def I(t, beta, ipusilon, gamma): 35 SEIRlist = odeint(SEIReq, (S0, E0, I0, R0), t, args = (beta, ipusilon, gamma)) 36 return SEIRlist[:,2] 37 38optparams, cov = optimize.curve_fit(I, t, normalized_cases, ) 39print('R0=',optparams[0]/optparams[2]) 40fitted = I(t, *optparams) 41 42plt.scatter(t, cases) 43plt.plot(t, fitted*Tokyo) 44plt.xlabel('the number of days from 2020/7/15') 45plt.ylabel('the number of confirmed cases in Tokyo') 46plt.show() 47f.close() # close the csv-file
試したこと
p0 = [10, 10, 10]等、初期条件を変えてみました。
上記にも述べましたが曲線は折れました。
補足情報(FW/ツールのバージョンなど)
あなたの回答
tips
プレビュー