質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
86.02%
Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

Q&A

解決済

Python SIRモデル 感染者予測

SHIROPANDA
SHIROPANDA

総合スコア2

Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

3回答

0グッド

0クリップ

333閲覧

投稿2022/10/31 08:13

編集2022/10/31 08:19

前提

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

実現したいこと

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

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

厚生労働省がサイトで出している新規感染者のデータを抽出しPythonでファイルの読み込みをしてプロットまではできていますが、フィットさせようとしてる曲線が全く立ち上がらず一直線になってしまっています。画像がうまくアップロードできないので後付けでそちらの画像をアップロードします。
イメージ説明
青いプロットが今回の実際のデータとなっており、下の真っ直ぐに伸びている直線が最小二乗法によって導出した曲線となっています。

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

該当のソースコード

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_1st_wave1.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; R0 = 0.0 23 24# SIR differential equation 25# S = SIR[0], I = SIR[1], R = SIR[2] 26def SIReq(SIR, t, beta, gamma): 27 dSdt = -beta*SIR[0]*SIR[1] 28 dIdt = beta*SIR[0]*SIR[1] - gamma*SIR[1] 29 dRdt = gamma*SIR[1] 30 31 return [dSdt, dIdt, dRdt] 32 33def I(t, beta, gamma): 34 SEIRlist = odeint(SIReq, (S0, I0, R0), t, args = (beta, gamma)) 35 return SEIRlist[:,1] 36 37optparams, cov = optimize.curve_fit(I, t, normalized_cases) 38print('R0=',optparams[0]/optparams[1]) 39fitted = I(t, *optparams) 40 41plt.scatter(t, cases) 42plt.plot(t, fitted*Tokyo) 43plt.xlabel('the number of days from 2020/2/18') 44plt.ylabel('the number of confirmed cases in Tokyo') 45plt.show() 46f.close() # close the csv-file

試したこと

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

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

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

以下のような質問にはグッドを送りましょう

  • 質問内容が明確
  • 自分も答えを知りたい
  • 質問者以外のユーザにも役立つ

グッドが多くついた質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

気になる質問をクリップする

クリップした質問は、後からいつでもマイページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

下記のような質問は推奨されていません。

  • 間違っている
  • 質問になっていない投稿
  • スパムや攻撃的な表現を用いた投稿

適切な質問に修正を依頼しましょう。

SHIROPANDA

2022/10/31 08:20

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

2022/10/31 09:44

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

2022/10/31 12:55

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

2022/11/01 00:12

> 用いている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ファイルをアップロードしなくても、コードで使ってるデータ全体を回答者と共有できると思います

回答3

2

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

 

python

1f = open("./COVID_1st_wave1.csv", "r", encoding="UTF-8", errors="", newline="" ) 2fcsv = csv.reader(f, delimiter=",", doublequote=True, lineterminator="\r\n", quotechar='"', skipinitialspace=True) 3next(f) # skip to the header of the csv-file

↓ 変更

python

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

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

イメージ説明

 
同じデータで、

python

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

↓ 変更

python

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

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

イメージ説明

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

イメージ説明

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

投稿2022/11/01 01:11

jbpb0

総合スコア7521

SHIROPANDA❤️を押しています
melian😍を押しています

良いと思った回答にはグッドを送りましょう。
グッドが多くついた回答ほどページの上位に表示されるので、他の人が素晴らしい回答を見つけやすくなります。

下記のような回答は推奨されていません。

  • 間違っている回答
  • 質問の回答になっていない投稿
  • スパムや攻撃的な表現を用いた投稿

このような回答には修正を依頼しましょう。

回答へのコメント

SHIROPANDA

2022/11/04 04:37

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

1

ベストアンサー

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

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

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

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

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

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

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

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

投稿2022/10/31 09:39

編集2022/10/31 09:54
PondVillege

総合スコア1066

can110👍を押しています

良いと思った回答にはグッドを送りましょう。
グッドが多くついた回答ほどページの上位に表示されるので、他の人が素晴らしい回答を見つけやすくなります。

下記のような回答は推奨されていません。

  • 間違っている回答
  • 質問の回答になっていない投稿
  • スパムや攻撃的な表現を用いた投稿

このような回答には修正を依頼しましょう。

回答へのコメント

PondVillege

2022/10/31 10:01

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

2022/10/31 13:01

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

2022/11/01 00:52 編集

最小二乗法では理想の関数 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

2022/10/31 23:52 編集

質問者さん > 何の初期値なのか 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

2022/11/04 03:10

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

2022/11/04 03:17

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

2022/11/04 09:06

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

2022/11/04 09:23 編集

> optparamsは変数(β, γ)のリスト optparams, cov = optimize.curve_fit(I, t, normalized_cases) のすぐ下に print(type(optparams)) を追加して実行したら分かりますが、「リスト」ではなくて「numpy配列」(<class 'numpy.ndarray'>)です

1

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

投稿2022/10/31 08:30

dark-eater-kei

総合スコア1248

SHIROPANDA😍を押しています

良いと思った回答にはグッドを送りましょう。
グッドが多くついた回答ほどページの上位に表示されるので、他の人が素晴らしい回答を見つけやすくなります。

下記のような回答は推奨されていません。

  • 間違っている回答
  • 質問の回答になっていない投稿
  • スパムや攻撃的な表現を用いた投稿

このような回答には修正を依頼しましょう。

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
86.02%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問

同じタグがついた質問を見る

Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。