Python SIRモデル 感染者予測

Question

前提

Pythonで新型コロナウイルス感染者の予測をおこなっています。
プログラム初学者のためSIRモデルを用いて過去のデータで数式に当てはめています。

実現したいこと

最小二乗法を用いて実際の感染者のデータにフィットするように曲線を描きたいと考えています。

発生している問題・エラーメッセージ

厚生労働省がサイトで出している新規感染者のデータを抽出しPythonでファイルの読み込みをしてプロットまではできていますが、フィットさせようとしてる曲線が全く立ち上がらず一直線になってしまっています。画像がうまくアップロードできないので後付けでそちらの画像をアップロードします。

青いプロットが今回の実際のデータとなっており、下の真っ直ぐに伸びている直線が最小二乗法によって導出した曲線となっています。

エラーメッセージは出ておりません

該当のソースコード

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy import optimize
import csv

# loading csv-file
f = open("./COVID_1st_wave1.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 2020
normalized_cases = np.array(cases, dtype = float)/Tokyo
days = len(cases)
t = np.arange(days)
# initial values
I0 = 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/18')
plt.ylabel('the number of confirmed cases in Tokyo')
plt.show()
f.close() # close the csv-file

試したこと

今回は2月18日から5月23日までの範囲で実際のデータにフィットさせようとしましたがうまくいかなかったのですが、1月24日から5月23日までの範囲ではうまくいきました。
また以前も少し似たような質問をしたことがあり、データの範囲の設定の際に必ず新規感染者が1人はいるところから始めないといけないことがわかっています。
よろしくおねがいいたします。

補足情報（FW/ツールのバージョンなど）

ここにより詳細な情報を記載してください。

Comments

[SHIROPANDA]

初学者なので伺いたいのですがこのサイトでは用いているcsvファイルのアップロードも可能だったりするのでしょうか。

[ps_aux_grep]

> 初学者なので伺いたいのですがこのサイトでは用いているcsvファイルのアップロードも可能だったりするのでしょうか。
残念ながら，それは不可能です．
かわりに，表を描くことはできます．また，コピペで解答者側が再現しやすいよう，コードブロックにCSVを書くのが主流です．いずれにしろ全データを示す必要はないです(今回は必要かもだけど)

[SHIROPANDA]

情報ありがとうございます。機会がまたあると思うのでその時載せてみようと思います。

[jbpb0]

> 用いているcsvファイルのアップロード

https://teratail.com/questions/ifoq8glk78r6rm
の「質問へのコメント」に、melianさんがcsvファイルを使わずにコードを実行する方法を書いてて、その後に私が、それを参考にして下記のコード変更をコメントに書いてます

f = 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
↓ 変更
import pandas as pd
url = 'https://covid19.mhlw.go.jp/public/opendata/newly_confirmed_cases_daily.csv'
df = pd.read_csv(url, parse_dates=['Date'])
df_Tokyo = df.loc[df['Date'].between('2020-02-10', '2020-05-20'), ['Date', 'Tokyo']]
fcsv = np.array(df_Tokyo)

この質問で使ってるデータの場合は、コードを上記のように変更して、抽出する期間の日付を適切に設定すれば、csvファイルをアップロードしなくても、コードで使ってるデータ全体を回答者と共有できると思います

Answer 1

1月24日から5月23日までの範囲ではうまくいきました。

　

f = open("./COVID_1st_wave1.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

↓ 変更

import pandas as pd
url = 'https://covid19.mhlw.go.jp/public/opendata/newly_confirmed_cases_daily.csv'
df = pd.read_csv(url, parse_dates=['Date'])
df_Tokyo = df.loc[df['Date'].between('2020-01-24', '2020-05-23'), ['Date', 'Tokyo']]
fcsv = np.array(df_Tokyo)

として、
厚生労働省が発表しているオープンデータ
から2020-01-24〜2020-05-23のTokyoのデータを抽出して、質問のコードのまま「optimize.curve_fit」の「初期値」を設定しなかった場合の結果

　
同じデータで、

optparams, cov = optimize.curve_fit(I, t, normalized_cases)

↓ 変更

p0 = [20, 20]
optparams, cov = optimize.curve_fit(I, t, normalized_cases, p0)

として、「optimize.curve_fit」の「初期値」を「p0 = [20, 20]」に設定した場合の結果

　
同じデータで、「optimize.curve_fit」の「初期値」を「p0 = [15, 15]」に設定した場合の結果

　
初期値にかなり依存してますよね

Comments

[SHIROPANDA]

初期値の重要性に気付けました。ありがとうございます。

Answer 2

以前他の方の回答で申し上げたことがあるのですが，そのような結果でも上記データを目的関数でfittingした結果で間違いありません．

最小二乗法で得られる解において，与えた初期状態p0

Initial guess for the parameters (length N). If None, then the initial values will all be 1

が悪かったため，局所的最適解を得ています．

現状，「初期値が全て1で初期化された局所的最適解は望んだ関数を描かない」という情報が得られたことになるので，これ以外の初期値で，試行すると良いでしょう．

もし自身で当たっていると思しきbeta, gammaを知っているなら，その値をp0に入れてやってください．運が良ければ最適解にシフトするものと思われます．

値が何もわからないのであれば，引数の数だけランダムに初期値を生成して与えると良いでしょう．

初期値が全て1でフィットする関数なんてたかが知れており，正弦波ですらフィットしないソルバを扱う中で，解に近い値を入れざるを得ないのは当然の事態と思います．

Comments

[ps_aux_grep]

https://teratail.com/questions/ifoq8glk78r6rm
前回の質問とその回答をみましたが，同様の内容を答えることになってしまいました．
解決方法はこれ以上でもこれ以下でもありません．適切な初期値を与えてください．

[SHIROPANDA]

すいません。数学があまり得意ではないのですがこの初期値と呼ばれるものは何の初期値なのかがわかっていません。もしよろしければ詳しく解説していただいてもよろしいでしょうか。不躾な質問であったら申し訳ございません。

[ps_aux_grep]

最小二乗法では理想の関数
y = f(x, β,γ)
に対して誤差e(x)が乗ったデータxi, yi
yi = f(xi, β,γ) + e(xi)
において，移項して得られる誤差を
e(xi,β,γ) = yi - f(xi, β,γ)
としてこの誤差関数e(xi,β,γ)の二乗和
E(β,γ) = Σ_i e(xi,β,γ)^2 =Σ_i {yi - f(xi, β,γ)}^2
が尤も小さくなるようなβ,γを探索します．

演算が簡単な関数であれば，二乗和した関数E(β,γ)を微分して0になるβ,γを手計算で導出するのが普通です．
ググって出てくるのも大抵は y = a * x + b における a, b の最適化で，教科書的な内容なので必ず目を通しましょう．
https://analytics-notty.tech/derivation-of-least-squares-equation/

curve_fitは，与えた初期値p0=[β0,γ0]（デフォルトで[1, 1]）を起点に誤差関数が尤も小さくなる場所の探索を始めるようなアルゴリズムになっています．誤差関数の値を小さくしたい．という現状の要件では，関数の傾きが負の方向に[β,γ]を更新しまくると良いでしょう．これを繰り返せばいつかは理想の関数にフィットするはずです．

更新アルゴリズムのイメージとしては，霧がかって見晴らしの悪い山を降りるときに，初期地点[β0,γ0]から「足場の傾きだけをみて下山する」ような処理です．足を一歩だけ下山する方向に踏み出して[β1,γ1]，さらにその足場の傾きから下山できる方向を考え足を踏み出して[β2,γ2]と進みます．これを繰り返すのがcurve_fitです．

式に当てはめると
E(β0,γ0) > E(β1,γ1) > E(β2,γ2) > ... > E(βn,γn)
となるように値を更新していきます．この誤差関数がこれ以上小さくならないような[βn,γn]を得たら更新終了，理想の関数にフィットしたと考えて最適解optparamsに書き出されます．

では，実は下山してるのではなく，山腹に存在してしまっている窪みに向かっているとしたらどうでしょう．カルデラみたいなイメージです．初期位置が悪かったがために，下山方向を誤ってしまったことになりますよね．窪みにハマってしまったら最後，前後左右どこに行っても登ることになります．窪みの中から脱出して下山するには，誤差を最小にする．というはずだったアルゴリズムに対して一度は誤差が大きい方向に行く(登山する)ことになり矛盾します．なので窪みから抜け出すことなくアルゴリズムは停止し，最適解としてoptparamsが吐き出されます．ここで得られた最適解β,γを，局所的最適解と呼びます．

これを回避すべく，初期位置は「ちゃんと下山できる場所からスタートしないといけない」．というのがセオリーになります．ちゃんと下山できたとき，すなわち関数の定義域全域で誤差の2乗和関数E(β,γ)が最小になったときのβ,γを大域的最適解と呼びます．

基本的に求めなければならないのは大域的最適解ですので，局所的最適解が得られた(と思った)場合は導出アルゴリズムや目的関数，初期値等のハイパーパラメータを見直す必要があります．

導出アルゴリズム，すなわち山の下り方も1パターンではありません． https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html#scipy.optimize.least_squares 効率の良い下山をする他のアルゴリズムの使用だって可能です．

[jbpb0]

質問者さん

> 何の初期値なのか

https://teratail.com/questions/ifoq8glk78r6rm
の「質問へのコメント」に私が、「初期値」は

> optparams, cov = optimize.curve_fit(I, t, normalized_cases)

で決定される(はずの)曲線のパラメータ「optparams」の初期値のことで、その質問のコードでは初期値を指定してないので全て「1」になってて、それが適切ではないのかもしれない、と書いてます

「optimize.curve_fit」は、「初期値」からスタートして、そこからちょっとずつ値を変えながら、「optparams」の最適値を探索します

【追記】
「初期値」をいろいろ変えたらどうなるかも、
https://teratail.com/questions/ifoq8glk78r6rm
の「質問へのコメント」に私が書いてます

[SHIROPANDA]

数弱で申し訳ないです。psさんのおっしゃってることが何となくは把握できました。
初期値の設定をうまく行うことによってうまく関数を導いてみたいと思います。
ありがとうございます。

[SHIROPANDA]

すいません
もう一つ質問なんですがoptparamsは変数(β, γ)のリストとして機能しているのでしょうか

[ps_aux_grep]

optparamsの0番目に最適解β，1番目にγが入っています．
fitted = I(t, *optparams)
としている箇所で，入っている値はこれらoptparams = [β, γ]というリストなので
fitted = I(t, *[β, γ])
になっており，インタプリンタが展開して
fitted = I(t, β, γ)
にしてくれます．

[jbpb0]

> optparamsは変数(β, γ)のリスト

optparams, cov = optimize.curve_fit(I, t, normalized_cases)

のすぐ下に

print(type(optparams))

を追加して実行したら分かりますが、「リスト」ではなくて「numpy配列」(<class 'numpy.ndarray'>)です

Answer 3

ちゃんとfittingできていると思います