【Pythonのコードあり】『現代数学の捉え方「代数編」』にSIRモデルの紹介があったのでPythonでコードを書いた。ちなみに、SIRモデルには厳密解がある。

2020年10月28日 00:23

どうも、こんにちは（@t_kun_kamakiri）

本日は、「現代数学の捉え方「代数編」」を読んでいました。

サイエンス社が出している「数理科学」は毎月気になったものは買ってじっくり読んでいます。
が、今回は数学がメイントピックなので読みこなすのは無理でした( ﾉД`)ｼｸｼｸ…

その代わり、コラムでは感染者数の推移をモデル化するというのがあって読んでいまして・・・・

そしてこちらのツイート。

「現代数学の捉え方［代数編］」を読んでいました。

SIRモデルの厳密解について原論文の紹介がありました。

Exact analytical solutions of the Susceptible-Infected-Recovered (SIR) epidemic
model and of the SIR model with equal death and birth rateshttps://t.co/OJjT6yiDT9 https://t.co/3TcggGjL3D
— カマキリ🐲@Django勉強中 (@t_kun_kamakiri) October 27, 2020

感染者数の推移に関しては、SIRモデルというのが有名だそうですね。

本記事では、
SIRモデルについてはびっくりすくらいさらっと流して、Pythonのコードを紹介します。PythonのコードをGoogle Colabにペタペタ貼っていけば計算実行できます('ω')

SIRモデルの概要

簡単にまとめると以下です。

モデルの概要は理解したので、Pythonでコードを書いて感染者数の推移をグラフにします。

Google Colabを使っています。

PythonでSIRモデルを解く

以下にPythonのコードを書いていきます。

まずは初期設定をします。

import matplotlib.pyplot as plt
import numpy as np
import japanize_matplotlib

####################
##### 初期設定 #####
####################

#観測日数
t = 250

#分割数
nt = 10*t
dt = t/nt

t = np.linspace(0,t,nt)
dt

次に、パラメータの設定を行います。
こちらのパラメータは感染者数推移に効いてくるので、追々パラメータ化して条件を変えたスタディをしたいと思います（後日）

##########################
##### パラメータ設定 #####
##########################

N = 1E6
m = 10
p = 0.02
d = 14
beta = m*p/N
gamma = 1/d

##########################
######## 初期値 ##########
##########################

I0 = 100
R0 = 0
S0 = N - I0 - R0

次に、メインの連立微分方程式の部分を関数で用意します。

#微分方程式を関数で設定#####

dSdt = lambda S , I , R , t: -beta*S*I
dIdt = lambda S , I , R , t: beta*S*I - gamma*I
dRdt = lambda S , I , R , t: gamma*I

↑こちらのように関数に名前を付けなくても関数を設定することができます。
あまり慣れない方はこちらの記事をお読みください。

次に、微分方程式を積分する際のルンゲクッタ（4次）の部分のコードです。

#空のリストを用意
S = np.empty(nt)
I = np.empty(nt)
R = np.empty(nt)

S[0] = S0
I[0] = I0
R[0] = R0


for it in range(len(t)-1):
 kS1 = dt * dSdt( S[it], I[it], R[it], t[it])
 kI1 = dt * dIdt( S[it], I[it], R[it], t[it])
 kR1 = dt * dRdt( S[it], I[it], R[it], t[it])

 kS2 = dt * dSdt( S[it] + kS1/2, I[it] + kI1/2, R[it] + kR1/2, t[it])
 kI2 = dt * dIdt( S[it] + kS1/2, I[it] + kI1/2, R[it] + kR1/2, t[it])
 kR2 = dt * dRdt( S[it] + kS1/2, I[it] + kI1/2, R[it] + kR1/2, t[it])

 kS3 = dt * dSdt( S[it] + kS2/2, I[it] + kI2/2, R[it] + kR2/2, t[it])
 kI3 = dt * dIdt( S[it] + kS2/2, I[it] + kI2/2, R[it] + kR2/2, t[it])
 kR3 = dt * dRdt( S[it] + kS2/2, I[it] + kI2/2, R[it] + kR2/2, t[it])

 kS4 = dt * dSdt( S[it] + kS3, I[it] + kI3, R[it] + kR3, t[it])
 kI4 = dt * dIdt( S[it] + kS3, I[it] + kI3, R[it] + kR3, t[it])
 kR4 = dt * dRdt( S[it] + kS3, I[it] + kI3, R[it] + kR3, t[it])

 S[it + 1] = S[it] + (kS1 + 2*kS2 + 2*kS3 + kS4)/6 
 I[it + 1] = I[it] + (kI1 + 2*kI2 + 2*kI3 + kI4)/6
 R[it + 1] = R[it] + (kR1 + 2*kR2 + 2*kR3 + kR4)/6

これで、S,I,Rの時間推移が解けたことになるので、グラフ化します。

#グラフの作成
plt.figure(figsize=(10, 6)) 

plt.grid()
plt.xlabel('date(day)' , fontsize=16)
plt.ylabel('number', fontsize=16)
plt.xticks(fontsize=16)
plt.yticks(fontsize=16)

plt.plot(t, S, color = 'green', label='S(t)', linewidth = 2)
plt.plot(t, I, color = 'red', label='I(t)', linewidth = 2)
plt.plot(t, R, color = 'blue', label='R(t)', linewidth = 2)
plt.legend(fontsize=32)

SIRモデルを今更読んでいる。 pic.twitter.com/a217ijC8Dq
— カマキリ🐲@Django勉強中 (@t_kun_kamakiri) October 27, 2020

グラフ化ができました。

ちなみに、Pythonって何で学んだらいいの？って思っている方は、こちらのやさしいPythonの参考書を1週間くらいで読んだらだいたいわかります('ω')

あとは、Pythonは外部のライブラリを使って色々することが多いので、外部のライブラリを使う練習をした方が良いですね。

Excelをよく使うなら以下の参考書とかはお勧めです。

あとは、スクレイピングと機械学習と自然言語処理とかの超簡単な話は以下の参考書に載っています。
あまり深くはやっていないので、難しくはないですが本格的にやろうとすると物足りないです・・・・が、入門程度にはちょうど良いです。

SIRモデルは厳密解がある

「現代数学の捉え方「代数編」」で紹介されていましたが、SIRモデルは厳密解があるみたいですね。

2014年の論文がarXivに出ています。

タイトル
Exact analytical solutions of the Susceptible-Infected-Recovered (SIR) epidemic model and of the SIR model with equal death and birth rates.

SIRの解の求め方はそれほど難しくなく導出できます。

感染者数のデータ

※ちなみに何の感染者数のデータかは書いていません。
noteに警告が出るため。

感染者数のデータは以下にまとめられているので、SIRモデルと実際のデータを比較して分析をしたいですね。

本日は以上となります(^^♪

Twitter➡@t_kun_kamakiri
ブログ➡宇宙に入ったカマキリ（物理ブログ）
ココナラ➡物理の質問サポートサービス

この記事が気に入ったらサポートをしてみませんか？