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

Question

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-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 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/13')
plt.ylabel('the number of confirmed cases in Tokyo')
plt.show()
f.close() # close the csv-file

新型コロナウイルスの新規感染者数のデータをSIRモデルの数式を用いて
最小二乗法を実行して、一つのグラフに実際のデータとその曲線とを表したのですが
曲線がうまく表示されません。

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

/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/minpack.py:833: OptimizeWarning: Covariance of the parameters could not be estimated
  warnings.warn('Covariance of the parameters could not be estimated',

と表記されています。
### 該当のソースコード

optparams, cov = optimize.curve_fit(I, t, normalized_cases)
この辺りに原因があるのではないかとみています。

### 試したこと

ファイルもうまく開けているので何が原因なのかがわかっていません
他のファイルではうまくいくこともあります。

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

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

Comments

[melian]

COVID-first_wave.csv ですが、厚生労働省が発表しているオープンデータを加工したものでしょうか？具体的な内容を開示してもらえると、こちらでも動作確認ができるかと思います。

[SHIROPANDA]

厚生労働省が発表しているオープンデータのうち、2月中旬から5月中旬にかけてまでのデータをエクセルでまとめてcsvファイルとして保存しました。
まとめたと言いましても東京の情報だけをコピペしただけです。
プログラミング初心者のため色々と至らぬ点があると思いますのでたくさんの質問よろしくお願いいたします。

[meg_]

コードは「コードの挿入」で記入してください。

[jbpb0]

pythonのコードの一番最初の行のすぐ上に
```python
だけの行を追加してください
また、pythonのコードの一番最後の行のすぐ下に
```
だけの行を追加してください

または、
https://teratail.storage.googleapis.com/uploads/contributed_images/56957fe805d9d7befa7dba6a98676d2b.gif
を見て、そのようにしてみてください

現状、コードがとても読み辛いです
質問にコードを載せる際に上記をやってくれたら、他人がコードを読みやすくなり、コードの実行による現象確認もやりやすくなるので、回答されやすくなります

[SHIROPANDA]

質問のようにおこなってみました。見やすくなったと思われます。ありがとうございます

[SHIROPANDA]

「やってほしいことだけを記載した丸投げの質問」という指摘を受けました。
やってほしいことなど記載していないです。
日本語が読めない方がいたので再度記載します。
どのコードの部分が原因であるかを知れたら自分なりに対処するつもりです。まだプログラミング言語を学んで1ヶ月ですので優しくみてくれると嬉しいです。
よろしくお願いいたします。

[meg_]

> やってほしいことなど記載していないです。
> 日本語が読めない方がいたので再度記載します。
> どのコードの部分が原因であるかを知れたら自分なりに対処するつもりです。まだプログラミング言語を学んで1ヶ月ですので優しくみてくれると嬉しいです。

デバック作業は工数がかかることも多いです。そのため質問者ご自身でどこまでデバッグされたのか？を回答者としては知りたいところかと思われます。さらに、コードを検証するためには今回の質問の場合ですとデータを用意することから始めなければなりません。データの入手方法についてもこちらのコメント欄には書いていますが、質問本文では言及されていません。あとは「曲線がうまく表示されません。」について具体的にグラフの画像があると第三者に伝わりやすいと思います。（「うまく表示されません」だけでは”どうなるはず”が”どうだった”のか分からないので）

「どのコードの部分が原因であるかを知れたら自分なりに対処するつもりです。」は質問に追記されると良いかと思います。（質問は編集できます）プログラミング初心者については質問に付ける”初心者マーク”がありますので活用されると良いです。（コードをパッと見た感じでは初心者感はありませんでしたので）

[melian]

Teratail には「やってほしいことだけを記載した丸投げの質問」などの指摘をカジュアルに付けてしまうユーザが一定数存在する様です。気になさることもないでしょう。

[SHIROPANDA]

ものすごくわかりやすい質問ありがとうございます。
訳があって基礎的な部分を飛ばし、このようなレベルのプログラミングを扱っております。
megさんの指摘くださったものに関してはいますぐ解答できるような能力は自分にはありませんので、もう一度自分なりにまとめてから再度質問しようと思います。

このサイトを知ったのは今日で、このような指摘をする方が一定数いるのも初めて知りました。
それを知った上でこのサイトを有効に今後活用させていただきます。
親切にありがとうございます。

[SHIROPANDA]

melian さん　
ありがとうございます

[jbpb0]

> warnings.warn('Covariance of the parameters could not be estimated',

は、「scipy.optimize.curve_fit」の初期値が適切でない場合に出ることがあります

曲線とデータのグラフを見て、曲線がデータにそこそこ合うように曲線のパラメータを決めてから、そのパラメータを初期値として使ってみてください

参考
https://teratail.com/questions/53bfrobbpegf9r
https://qiita.com/kon2/items/6498e66af55949b41a99
の「いよいよ関数でフィッティング」
https://www.haya-programming.com/entry/2019/03/12/052716
の「対数の底」

[melian]

厚労省のオープンデータから 2020/2/10 〜 2020/5/20 までの東京の新規陽性者数をみると、jbpb0 さんの指摘の通り、初期値が 0 になっているのだと思います。

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']]
df_Tokyo.head(10)

Date Tokyo
25 2020-02-10 0
26 2020-02-11 0
27 2020-02-12 0
28 2020-02-13 0
29 2020-02-14 0
30 2020-02-15 0
31 2020-02-16 0
32 2020-02-17 0
33 2020-02-18 3
34 2020-02-19 3

[jbpb0]

melianさん

私がコメントに書いた「初期値」は、データの始まりの値ではなく、

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

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

ちなみに、melianさんのコメントを参照して、質問のコードを

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)

と変更して実行したら、

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

の結果は、

optparams：[1. 1.]
cov：
[[inf inf]
[inf inf]]

で、パラメータは初期値から変わってませんでした

[melian]

> jbpb0 さん
ごめんなさい、意味を取り違えていました。

[jbpb0]

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

の適切な初期値を探索しようと、melianさんのコメントを参考に、質問のコードを

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)

optparams, cov = optimize.curve_fit(I, t, normalized_cases)
↓ 変更
optparams = [1, 1]

と変更して、「optparams = [1, 1]」の数値「1, 1」をいろいろ変えてみましたが、

> fitted = I(t, *optparams)

の計算結果「fitted」は全部「0」で変わりませんでした

そこで、

df_Tokyo = df.loc[df['Date'].between('2020-02-10', '2020-05-20'), ['Date', 'Tokyo']]
↓ 変更
df_Tokyo = df.loc[df['Date'].between('2020-02-18', '2020-05-20'), ['Date', 'Tokyo']]

と日付を変えてみたら、「fitted」は「0」以外の値が計算されました

どうやら、melianさんが指摘してた、最初の日のデータが「0」だとダメなようです

> Date Tokyo
25 2020-02-10 0
26 2020-02-11 0
27 2020-02-12 0
28 2020-02-13 0
29 2020-02-14 0
30 2020-02-15 0
31 2020-02-16 0
32 2020-02-17 0
33 2020-02-18 3
34 2020-02-19 3

「SIRモデル」というものをよく知らないのですが、そういうものなのですかね

[jbpb0]

df_Tokyo = df.loc[df['Date'].between('2020-02-10', '2020-05-20'), ['Date', 'Tokyo']]
↓ 変更
df_Tokyo = df.loc[df['Date'].between('2020-02-18', '2020-05-20'), ['Date', 'Tokyo']]

と日付を変えて、「optparams = [1, 1]」の数値「1, 1」をいろいろ変えてみたところ、「optparams = [20.2, 20.1]」でそこそこ合ってるようにグラフで見えたので、下記のように変更したら、実行できました

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

【追記】
p0 = [10, 10]

でも実行できました

【追記2】
p0 = [1, 1]

ではダメでした
([1, 1]は初期値を指定しない場合と同じ)

[jbpb0]

以上をまとめると、下記の両方の条件が必要なようです
・「cases」(新規感染者数？)の最初の(日の)数値が「0」以外
・「scipy.optimize.curve_fit」の初期値を適切に設定

[SHIROPANDA]

解答の方ありがとうございます。初学者のため、jbpb0さんが変更してくださったcodeをすぐ理解できる能力がないので再度勉強しなしてからまた伺おうと思います。初期条件が0が無理ということでしたのでファイルで感染者が0じゃない箇所から始めたいと考えています。

[SHIROPANDA]

自分の描いたコードではそのままにし、初期条件をcsvファイルの方で感染者が０人でないところから始めたところうまくいきました。このまま続けてみます。本当にありがとうございます。

[jbpb0]

> 自分の描いたコードではそのままにし、初期条件をcsvファイルの方で感染者が０人でないところから始めたところうまくいきました。

https://teratail.com/questions/ut0nj7eedi8te6
を見ると、まだうまくいってないようなので、この質問の私の回答の通りに「scipy.optimize.curve_fit」の初期値を適切に設定してみてください
データを抽出する期間の日付と、パラメータ初期値の設定が適切なら、うまくいくはずです

【追記】
適切なパラメータ初期値は、データを抽出する期間の日付に依存するようです

Answer 1

下記の両方の条件が必要なようです
・「cases」(新規感染者数？)の最初の(日の)数値が「0」以外
・「scipy.optimize.curve_fit」の初期値を適切に設定

　
たとえば、
厚生労働省が発表しているオープンデータ
から2020-02-18〜2020-05-20のTokyoのデータを抽出した場合、

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

↓ 変更

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

としたら、質問のエラーは出ずに実行できました