お詫びと訂正
本記事で、数理モデルによる新型コロナウイルス感染者数の推移の分析において、感染者数の計算に不備があることが、ユーザーの方からのご指摘で分かりました。
このため設計通りの分析結果が得られていない状態で情報を提供しておりました。
私たちが直面している、非常に関心の高い内容にも関わらず、十分な検証を実施せず情報提供をしていましたことを深くお詫び申し上げます。
申し訳ございません。
感染者数の計算処理を修正し、分析結果、グラフ及びプログラム（github）を訂正致しました。
また、タイトルとサムネイルだけをご覧になって、誤解される方もいらっしゃるかもしれませんので、数理モデルが推定した収束時期は削除しました。
さらに、感染者データのCSVファイルが更新されていましたので、3月11日までの感染者数データをダウンロードして使用しています。
お気づきの点等ございましたら、ご指摘いただければ幸いでございます。
なお、本記事で使用している数理モデルは、実際の新型コロナウイルスの感染症流行を適切に記述しているわけではなく、得た観察を近似する論理的な説明に過ぎません。
従いました、数理モデルの結果は、実際とは異なることを改めてご理解頂きたく存じます。
最新のグラフをここからご覧いただけます。
四谷ラボのメンバーが、有志で作ってくれました。
新型コロナウイルス COVID-19の今後の感染予測 - 四谷ラボ
使用したCSVファイル
3/11までの感染者データを使用して、グラフを更新しています。
本編
こんにちは、四谷ラボのやましんです。 お気づきでしょうか？四谷ラボのロゴが完成しました！！ とても気に入ってます。みなさんにも気に入ってもらえると嬉しいです。
さて、新型コロナウイルス感染拡大防止のため、学校は臨時休業になってしまいましたね。いろんな意見があると思うけど、日本の感染状況を可視化したら、やむを得ない気がします。
もしかすると、在宅勤務しながら、このブログをご覧になっている人も多いのではないでしょうか？ 私は今日もマスクをして出勤しまーす。
そうそう、前回このタイトルで記事を書いたところ・・・。
新型ウイルスの感染者数を可視化して気づいたこと - 四谷ラボ公式ブログ
「何も気づいてねぇじゃねーか」という暖かいご声援を頂きました。ありがとうございます。
ごもっともです。・・・とほほ。その反省を踏まえて・・・
今回挑戦するのは・・・
日本の発症者数のピークと収束時期を数理モデルを使って推定します。
むろん、数理モデルは現実とは違うので、推定された結果を「信じるか信じないかはあなた次第・・」ではなく「信じないでください」。
もちろん、厚労省の発表を信じるかは「あなた次第です」
いきなり結論！！
今後の日本の新型コロナウイルス発症者数の推定
数理モデル曰く・・・
「日本の発症者数のピークは3月16日」
「ピーク時の発症者数は478人に至る（※）」
「完全に収束するのは9月以降」
（※）1日に478人発症するという意味ではなく、その日の発症者数が478人存在する。という意味です。
実際の観測データは、棒グラフで表現し、赤棒は発症者数。青棒は回復者数。黒棒は死亡者数。 一方、推定データは、線グラフで表現し、赤線は発症者数。青線は回復者数。黒線は死亡者数。
良くメディアで見る、感染者数のグラフは、回復者数、死亡者数を減算していないですが、このグラフでは感染者数から回復者数と死亡者数を引いたものを発症者数としています。
言い換えると、このグラフは、その日、何人の発症者、回復者、死亡者が存在しているかを意味しています。
海外からの感染者の流入が０であることを前提にしたモデルなので、実際にはピークはこの日より前になることは考えにくいと思います。
過去の観測データだけのグラフからもわかるように、日本の感染の広がりは始まったばかりで、今後の感染防止対策の効果を期待するばかりです。
今後の中国の新型コロナウイルス発症者の推定
多くの感染者を出している、中国ではどうでしょうか？
数理モデル曰く・・・
「完全に収束するのは8月以降」
グラフの見方は、日本と同じ。
過去の観測データからもわかるように、ピークは過ぎていると思われる。
もし、理系のパスポートを高校に置いてきた方は、高校に取りに戻って、先を読み進めるか、この記事に興味を持ってくれそうな方へ「シェア」してちょーだい。
ちなみに、僕は理系のパスポートを中二まで取りに戻りました。文系のパスポートは小学校に放置してます。
現在の感染状況を再確認
それでは順を追って、新型コロナウイルスの感染状況と数理モデルの適用について確認しましょう。 まずは、3月11日時点の日本と中国の感染状況のグラフを確認しておきましょう。
今回のデータはkaggleというサイトからダウンロードしました。 ちなみにKaggleは、世界中の機械学習・データサイエンスが集まるコミニティーです。
Novel Corona Virus 2019 Dataset | Kaggle
いろんなデータを参考にしてわかったけど、メディアによってデータがずいぶん異なります。
もちろん、国家レベルで感染者数を隠蔽している場合も考えられるけど、一旦、このデータを信用します。
日本の新型ウイルス発症者数の推移（観測データ）
まずは、日本の発症者数の状況を見てみよう。
赤棒は発症者数。青棒は回復者数。黒棒は死亡者数。
上述の「結論」で紹介したグラフと違って、過去の実際のデータのみで描画しています。
ん？なんか目がチカチカするぞ。グラフの色がどぎつ過ぎた。
う～ん、下の中国のグラフと比べると、発症者数は圧倒的に少ないものの、まだまだ増加傾向にあるように見えるね。
TOKYO2020は無事開催されるのだろうか？政府が慌てて対策するのも納得。
中国の新型ウイルス発症者数の推移（観測データ）
次は、感染者数の最も多い、中国のグラフはこちら。
同様に、赤棒は発症者数。青棒は回復者数。黒棒は死亡者数。
ふむふむ、中国では発症者数のピークは過ぎていて、減少傾向であるようにみえる。
それでは、今後の日本の発症者数の推移を数理モデルを使って、推定しよう！
数理モデル
感染症流行の数理モデルに、SEIRモデルというのがあることがわかった。今回はこのSEIRモデルを使って、日本の新型コロナウイルスのピークと収束時期を推定する。
他の感染症流行の数理モデルとして、潜伏期間を考慮しないSIRモデルや、コンピュータ上に仮想都市を設計して、その中に暮らす住人を動作させることでシミュレーションを行うMAS（Multi Agent Simulation）などもある。Unityのml-agentsを使ってMASを実装できそうですねぇ。
SEIRモデル
SEIRは下記の４つの頭文字をとっている。
Exposedは感染はしているけど、人に感染させる能力を持っていない状態。Infectiousは発症しており、人に感染させる能力を持っている状態。
今回のモデルでは、SEIRの状態に加えて「死亡」を含めた下記の５つの状態を常微分方程式で表現する。
パラメータは、それぞれ単位時間当たりの、Beta：伝達係数、Kappa：遷移数（E→I）、Gamma：回復率、Tau：死亡率。
遷移数、回復率、死亡率は、過去のデータから単純に計算する。これらを既知のパラメータ値として、常微分方程式に代入し、伝達係数を常微分方程式の解として導き出す。
導かれた伝達係数を使って、Infectious（発症者）の推移を推定し、実際の発症者数の推移との差が最小になる値を求めた。
厄介なのは、Susceptible（感染症に対して免疫を持たない者）の初期値がわからいこと。これも同様にSusceptibleの値を変化させて、実際の観測者数との差が最小になる値を探した。 （よりスマートな他の手法も試したが、これが一番うまくいった）
冒頭に書いたけど、これはあくまでも数理モデル。実際の世界とは違います。 この数理モデルでは、回復した人は二度と感染しないし、海外からの感染者の流入は考慮に入れてません。
プログラミング
プログラム
「三度の飯よりプログラム」という方は、こちらからクローンして、動かしちゃってください。 github.com
以降、プログラムの主なコードを紹介します。
CSVファイル読み込み
Kaggleから取得できるZIPファイルに格納されているCSVデータは５つある。今回はcovid_19_data.csvというファイルを使います。
他のCSVファイルには、感染者発生都市の緯度経度などの情報を含むものもあるので、感染ルートを地理的な観点から分析することもできそうだ。
実際のCSVファイルの書式は、こんな感じです。今回は使わないけど、都市名もあるよ。
| SNo | ObservationDate | Province/State | Country/Region | Last Update | Confirmed | Deaths | Recovered |
|---|---|---|---|---|---|---|---|
| 1 | 01/22/2020 | Anhui | Mainland China | 1/22/2020 17:00 | 1 | 0 | 0 |
| 2 | 01/22/2020 | Beijing | Mainland China | 1/22/2020 17:00 | 14 | 0 | 0 |
| 3 | 01/22/2020 | Chongqing | Mainland China | 1/22/2020 17:00 | 6 | 0 | 0 |
| 4 | 01/22/2020 | Fujian | Mainland China | 1/22/2020 17:00 | 1 | 0 | 0 |
| 5 | 01/22/2020 | Gansu | Mainland China | 1/22/2020 17:00 | 0 | 0 | 0 |
| 6 | 01/22/2020 | Guangdong | Mainland China | 1/22/2020 17:00 | 26 | 0 | 0 |
| 7 | 01/22/2020 | Guangxi | Mainland China | 1/22/2020 17:00 | 2 | 0 | 0 |
| 8 | 01/22/2020 | Guizhou | Mainland China | 1/22/2020 17:00 | 1 | 0 | 0 |
| 9 | 01/22/2020 | Hainan | Mainland China | 1/22/2020 17:00 | 4 | 0 | 0 |
CSVデータは定期的に更新されてるようなので、ご自身で最新のデータをダウンロードして、より精度の高い推定をしてみてはいかがでしょうか？
- CSVファイルを読み込む関数
def read_csv(filename): data = {} with open(filename) as f: reader = csv.reader(f) header = next(reader) for row in reader: for cnt, name in enumerate(header): if name not in data: data[name] = [] data[name].append(row[cnt]) return convert_count_by_country(data)
最初にヘッダーを読み込んで、各フィールド名（ObservationDate等）が、キーとなるディクショナリを構築してます。
- 扱いやすいデータに変換する関数
def convert_count_by_country(csv_data): data = {} for date, country, confirmed ,deaths, recovered in zip(csv_data['ObservationDate'], csv_data['Country/Region'],csv_data['Confirmed'],csv_data['Deaths'],csv_data['Recovered']): date = re.sub("[0-9]{2}:[0-9]{2}:[0-9]{2}","",date) if country not in data: data.setdefault(country,{date:[0, 0, 0]}) if date not in data[country]: data[country].setdefault(date, [0, 0, 0]) data[country][date][0] += float(confirmed) data[country][date][1] += float(deaths) data[country][date][2] += float(recovered) return data
「国」ｘ「日付」でアクセス可能な、ディクショナリを再構築します。後々データをこねくり回せるように、こういう構造のデータにしました。
SEIRモデルの実装
SEIRモデルと可視化機能をクラスEstimationInfectedPeopleとして実装した。以降、主要なメソッドを紹介する。
- 常微分方程式を定義する関数
def SEIR(self, vec, time, Beta): # vec[0]: S:Susceptible(未感染) # vec[1]: E:Exposed(潜伏感染) # vec[2]: I:Infected(発症) # vec[3]: R:Recovered(回復) # vec[4]: D:Died(死亡) # # Beta:伝達係数 # Kappa:遷移率（E→I） # Gamma:回復率 # Tau:死亡率 S = vec[0] E = vec[1] I = vec[2] R = vec[3] D = vec[4] N = S+E+I+R+D Kappa = 1/self.latent_period Gamma = self.recovery_rate Tau = self.mortality_rate return [-Beta*S*I/N, Beta*S*I/N-Kappa*E, Kappa*E - Gamma*I - Tau*I, Gamma*I, Tau*I]
遷移率、回復率、死亡率は、このクラスのコンストラクタ内で、既知のデータから計算して、常微分方程式に代入しました。
- 伝達係数Betaを与えて、S・E・I・R・Dのそれぞれの時系列データを取得する関数
def estimate(self, Beta): vec=odeint(self.SEIR,self.initParams,self.time,args=(Beta,)) est=vec[0:int(self.max/self.dt):int(1/self.dt)] return est
常微分方程式の解を求めるため、scipyのodeint()を使用した。
- 過去の観測データと推定データの差の合計を返却する関数
def func(self, params): est_i=self.estimate(params[0]) return np.sum(est_i[:,2] - self.infected * np.log(est_i[:,2]))
est_i[:,2]は、推定した発症者数の時系列データ。self.infectedは、実際の発症者数の観測データ。 計算しやすいようLogをとって、差の合計を計算している。
- susceptible（感染症に対して免疫を持たない者）の最適値を推定する関数
def getEstimatedParams(self): no_new_record_cnt = 0 max_fun = 0 bounds = [(0, None)] initParams = [0.001] for susceptible in range(self.confirmed[len(self.confirmed)-1], self.population ): self.initParams = [susceptible, 0, np.min(self.confirmed), 0, 0] estimatedParams = minimize(self.func,initParams,method="L-BFGS-B", bounds=bounds) if estimatedParams.success == True: if max_fun < -estimatedParams.fun: no_new_record_cnt = 0 max_fun = -estimatedParams.fun best_population = population print('Susceptible:',susceptible, ' Score:',max_fun) self.bestEstimatedParams = estimatedParams self.bestInitParams = self.initParams else: no_new_record_cnt += 1 if no_new_record_cnt > 1000: break return self.bestEstimatedParams
得られる推定発症者数が、実際の発症者数と近くなるよう、minimize()の戻り値が最小になるsusceptibleの値を探した。 常微分方程式の解として求める遷移率Betaは０より大きいため、パラメータの上限・下限を指定することができるL-BFGS-B法を指定して、minimize()を使った。
no_new_record_cnt は、処理時間の短縮のために設けた。最小値が連続で1000回更新されない場合は、その先、もう見つからないと決めつけている。このため局所解に陥る可能性を持っているが、最終的にはグラフを目視確認するので、これで良しとする。
可視化の実装
グラフ表示用に、次の関数を用意した。estimate()とほぼ同じだけど、過去の観測データの期間を６倍に延長してS・E・I・R・Dの未来の予測値を得れるようにした。
６倍は適当で、発症者数が収束する日が、グラフに収まるように調整して得た数値。
def estimate4plot(self, Beta): multiple = 6 v=odeint(self.SEIR,self.bestInitParams,np.arange(0, self.max * multiple, self.dt),args=(Beta,)) est=v[0:int(self.max * multiple/self.dt):int(1/self.dt)] return est
過去のデータを棒グラフ表示する関数。matplotlibを使用してます。使い方はググってください。
def plot_bar(self, ax):
未来の発症者数、回復者数、死亡者数の推定値を線グラフで描画する関数。
def plot_estimation(self, ax, estimatedParams):
クラスの呼び出し
population = 120000000 fig = plt.figure(figsize=(20,10),dpi=200) ax = fig.add_subplot(1, 1, 1) fig.suptitle('Infections of a new coronavirus in Japan', fontsize=30) Japan = EstimationInfectedPeople('Japan', population, data['Japan']) Japan.plot_bar(ax) Japan.save_plot('obcerved') estParams = Japan.getEstimatedParams() print(estParams) Japan.print_estimation(estParams) Japan.plot(ax, estParams) Japan.save_plot('estimation') ax.clear()
EstimationInfectedPeople()で、観測データ、日本人口を引数にインスタンスを生成する。 plot_bar()で過去データを棒グラフで表示。getEstimatedParams()で常微分方程式の最適解を求め、plot()で過去の観測データと未来の推定データを同じスケールで描画する。
すごく参考になったサイト
https://www.ism.ac.jp/editsec/toukei/pdf/54-2-461.pdf
感染症数理モデル事始め PythonによるSEIRモデルの概要とパラメータ推定入門 - Qiita
最後に
感染された方をご家族に持つ方の心情を想像すると、心が痛みます。
いち早く、新型ウイルスの感染拡大を抑えることができること。被害を最小限に留めることを心から願います。
私にできることがあればという思いで執筆しています。
