【COVID19×Python】SIRモデルをアニメーションを作って理解する-Med Python

できあがり

感染病の数学予測モデル「SIRモデル」について、Pythonでシミュレーションをつくってみました。
さらに、アニメーションも作成しました。

YouTubeにて解説動画を公開中

こちらのコードを元に解説動画を公開しています！
Premiere Proを初めて使って、編集してみました。
10分ちょうどの動画ですが、編集するのはその何倍も時間がかかりました...

コロナ撲滅ちゃんねる - YouTube

ちなみに、こちらの知り合いの方のチャンネルにて他にもCOVID-19関連の動画が公開されています。

合わせてぜひご覧ください！

全体のコードはこちらにて

ライブラリはこちら

SIRモデルとは？

伝染病の流行プロセスを説明する基本モデルのこと。
1927年、W.O. KermackとA.G. McKendrickの論文にて提案された。
A contribution to the mathematical theory of epidemics | Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character

SIRとはそれぞれ、

感染可能者（Susceptible）

感染者（Infected）

免疫保持者または死亡者（Recovered）

を表す。

また、下のような微分方程式で立式される。

今回はこちらの式をもとに、Pythonでモデルを実装し、アニメーションの作成まで行う。

Pythonで実装してみよう

流れ

ライブラリの読み込み

パラメータを設定

微分方程式を立てて、解く

グラフをつくる

アニメーションをつくる

ライブラリの読み込み

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

パラメータを設定

必要なパラメータは、

日数（期間）
人数
感染率・回復率

の3つ。

# 日数の最大値
day_max = 150
# グラフにおいて、何日ごとに数値を求めるか
dt = 0.01

# 人数
S_0 = 99
I_0 = 1
R_0 = 0
N_total = S_0 + I_0 + R_0  # 人口
ini_state = [S_0, I_0, R_0]  #[S[0], I[0], R[0]]

# 感染率・回復率
infection_rate = 0.271438  # 感染率
recovery_rate = 0.031  # 回復率

# SIRモデルに適応できるように
beta = infection_rate / N_total
gamma = recovery_rate

人数は、初期の感染者数、健常者数、回復者数を定める。3つを定めると、全体の人口も決まる。

感染率、回復率については、

COVID-19情報共有 — COVID19-Information sharing – By Prof. Dr. Aki-Hiro Sato

を参考にさせていただき、今回は

感染率＝0.271438
回復率=0.031

と定めることにした。

微分方程式を立てて、解く

# SIRモデルの立式
def SIR_EQ(v, t, beta, gamma):
    '''
    [v[0], v[1], v[2]] = [S, I, R]

    dSdt = - beta * S * I
    dIdt = beta * S * I - gamma * I
    dRdt = gamma * I
    '''

    dSdt = - beta * v[0] * v[1]
    dIdt = beta * v[0] * v[1] - gamma * v[1]
    dRdt = gamma * v[1]

    return [dSdt, dIdt, dRdt]

# 微分方程式を解く
times = np.arange(0, day_max, dt)
args = (beta, gamma)  # scipy.odeintに渡す、微分方程式の追加の引数

result = odeint(SIR_EQ, ini_state, times, args)

微分方程式は、scipy.integrate.odeintを使って解を求める。

関数の変数において、

betaを感染率
gammaを回復率

とする。

resultが計算結果となっており、numpy.ndarray形式となっている。

result

array([[9.90000000e+01, 1.00000000e+00, 0.00000000e+00],
       [9.89973096e+01, 1.00238003e+00, 3.10369000e-04],
       [9.89946129e+01, 1.00476565e+00, 6.21476643e-04],
       ...,
       [1.83014271e-02, 1.81093375e+00, 9.81707648e+01],
       [1.83005276e-02, 1.81037334e+00, 9.81713261e+01],
       [1.82996285e-02, 1.80981311e+00, 9.81718873e+01]])

PandasのDataFrameに変換すれば、より見やすくなる。

import pandas as pd
df = pd.DataFrame(result, columns=["susceptible", "infectious", "removed"])

# 数値の表示形式を変更する
# floatの表示方法を、小数点以下3桁で設定する
pd.options.display.float_format = '{:.3f}'.format

df

susceptible    infectious  removed
0    99.000  1.000   0.000
1    98.997  1.002   0.000
2    98.995  1.005   0.001
3    98.992  1.007   0.001
4    98.989  1.010   0.001
...    ... ... ...
14995    0.018   1.812   98.170
14996    0.018   1.811   98.170
14997    0.018   1.811   98.171
14998    0.018   1.810   98.171
14999    0.018   1.810   98.172

グラフをつくる

susceptible, infectious, removedごとにplot

まずはS, I, Rごとの日数に伴う人数の推移をグラフにしてみる。

susceptibleについての計算結果は、

result[:, 0]

array([9.90000000e+01, 9.89973096e+01, 9.89946129e+01, ...,
       1.83014271e-02, 1.83005276e-02, 1.82996285e-02])

のように取り出せるので、グラフは、

plt.plot(times, result[:, 0], color="blue")
plt.legend("Susceptible")  # 凡例
plt.title("Susceptible")  # タイトル
plt.xlabel("Day")  # x軸ラベル
plt.ylabel("Nubmer")  # y軸ラベル
plt.ylim([0, 100])  # y軸の表示範囲

このようにplotすることができる。

他の2つについても、同様に作成すると、

plt.plot(times, result[:, 1], color="orange")
plt.legend("Infectious")
plt.title("Infectious")
plt.xlabel("Day")
plt.ylabel("Nubmer")
plt.ylim([0, 100])

plt.plot(times, result[:, 2], color="green")
plt.legend("Removed")
plt.title("Removed")
plt.xlabel("Day")
plt.ylabel("Nubmer")
plt.ylim([0, 100])

全体をまとめてplot

3つのグラフをまとめると、

plt.plot(times, result)
plt.legend(['Susceptible','Infectious', 'Removed'])
plt.title("SIR model")
plt.xlabel("Day")
plt.ylabel("Nubmer")

パラメータを変えて、比較してみる

vulの値を変える

vulという新たな変数を定め、日本国内における、感染防止対策のレベルを感染率に反映させる。

(感染率) * vul

を感染率として更新する。

全く対策をしない自然の場合が、vul=1.0の状態であり、vulの値が小さいほど対策が成されている状態を表す。

同じく、

数値シミュレーションによる検討 – COVID-19情報共有 — COVID19-Information sharing

を参考にさせていただいた。

これまでのパラメータなどをまとめたコードは、

# パラメータ
day_max = 150  # 日数の最大値
dt = 0.01  # グラフにおいて、何日ごとに数値を求めるか

S_0 = 99
I_0 = 1
R_0 = 0
N_total = S_0 + I_0 + R_0  # 人口
ini_state = [S_0, I_0, R_0]  #[S[0], I[0], R[0]]

infection_rate = 0.271438  # 感染率
recovery_rate = 0.031  # 回復率

vul = 1.0
beta = infection_rate / N_total * vul
gamma = recovery_rate

times = np.arange(0, day_max, dt)
args = (beta, gamma)  # scipy.odeintに渡す、微分方程式の追加の引数

result = odeint(SIR_EQ, ini_state, times, args)
#plot
plt.plot(times, result)
plt.legend(['Susceptible','Infectious', 'Removed'])
plt.title("SIR model")
plt.xlabel("Day")
plt.ylabel("Nubmer")

vulの値を1.0, 0.5, 0.2, 0.1に変えてplotをつくると、

day_max = 150  # 日数の最大値
dt = 0.01  # グラフにおいて、何日ごとに数値を求めるか

S_0 = 99
I_0 = 1
R_0 = 0
N_total = S_0 + I_0 + R_0  # 人口
ini_state = [S_0, I_0, R_0] #[S[0], I[0], R[0]]

infection_rate = 0.271438  # 感染率
recovery_rate = 0.031  # 回復率

# vulの値を変えて
for vul in [1.0, 0.5, 0.2, 0.1]:
    # SIRモデルに適応できるように
    beta = infection_rate / N_total * vul
    gamma = recovery_rate

    #numerical integration
    times = np.arange(0, day_max, dt)
    args = (beta, gamma)  # scipy.odeintに渡す、微分方程式の追加の引数

    #Numerical Solution using scipy.integrate
    #Solver SIR model
    result = odeint(SIR_EQ, ini_state, times, args)
    #plot
    plt.plot(times,result)
    plt.legend(['Susceptible','Infectious', 'Removed'])
    plt.title("vul = " + str(vul))
    plt.xlabel("Day")
    plt.ylabel("Nubmer")
    plt.show()

日数を変えてみる

day_max = 300  # 日数の最大値
dt = 0.01  # グラフにおいて、何日ごとに数値を求めるか

S_0 = 99
I_0 = 1
R_0 = 0
N_total = S_0 + I_0 + R_0  # 人口
ini_state = [S_0, I_0, R_0] #[S[0], I[0], R[0]]

infection_rate = 0.271438  # 感染率
recovery_rate = 0.031  # 回復率

# vulの値を変えて
for vul in [1.0, 0.5, 0.2, 0.1]:
    # SIRモデルに適応できるように
    beta = infection_rate / N_total * vul
    gamma = recovery_rate

    #numerical integration
    times = np.arange(0, day_max, dt)
    args = (beta, gamma)  # scipy.odeintに渡す、微分方程式の追加の引数

    #Numerical Solution using scipy.integrate
    #Solver SIR model
    result = odeint(SIR_EQ, ini_state, times, args)
    #plot
    plt.plot(times,result)
    plt.legend(['Susceptible','Infectious', 'Removed'])
    plt.title("vul = " + str(vul))
    plt.xlabel("Day")
    plt.ylabel("Nubmer")
    plt.show()

アニメーションをつくる

matplotlibのanimationを用いて、冒頭のようなアニメーションを作ります。

アニメーション作成のイメージ

動かしたいパラメータを決める

それぞれのパラメータでの、静止画をplotする

フレームごとに動かす

今回、動かすパラメータはvulに設定し、対策の度合いを見ていきたいと思います。

from matplotlib import animation
from IPython.display import HTML

JupyterNotebook上で出力したい場合は、HTML形式で見ることができます。

まずはシンプルにアニメーションを作ってみます。

fig = plt.figure()
ims = []

day_max = 150  # 日数の最大値
dt = 0.01  # グラフにおいて、何日ごとに数値を求めるか

S_0 = 99
I_0 = 1
R_0 = 0
N_total = S_0 + I_0 + R_0  # 人口
ini_state = [S_0, I_0, R_0] #[S[0], I[0], R[0]]

infection_rate = 0.271438  # 感染率
recovery_rate = 0.031  # 回復率

for vul in [1.0, 0.5, 0.2, 0.1]:
    # SIRモデルに適応できるように
    beta = infection_rate / N_total * vul
    gamma = recovery_rate

    #numerical integration
    times =np.arange(0,day_max, dt)
    args = (beta, gamma)  # scipy.odeintに渡す、微分方程式の追加の引数
    # 計算
    result = odeint(SIR_EQ, ini_state, times, args)

    img = plt.plot(times,result)
    plt.title("Infection rate transition")
    plt.legend(['Susceptible','Infectious', 'Removed'])

    ims.append(img) # グラフを配列に追加

# 100枚のプロットを 100ms ごとに表示するアニメーション
ani = animation.ArtistAnimation(fig, ims, interval=300)
# plt.show()
HTML(ani.to_jshtml())

できました！

が、フレームごとに色が変わってしまっているので、統一したいと思います。

colorを統一する場合、

fig = plt.figure()
ims = []
flag_legend = True  # 凡例描画のフラグ

day_max = 150  # 日数の最大値
dt = 0.01  # グラフにおいて、何日ごとに数値を求めるか

S_0 = 99
I_0 = 1
R_0 = 0
N_total = S_0 + I_0 + R_0  # 人口
ini_state = [S_0, I_0, R_0] #[S[0], I[0], R[0]]

infection_rate = 0.271438  # 感染率
recovery_rate = 0.031  # 回復率

for vul in [1.0, 0.5, 0.2, 0.1]:
    # SIRモデルに適応できるように
    beta = infection_rate / N_total * vul
    gamma = recovery_rate

    #numerical integration
    times =np.arange(0,day_max, dt)
    args = (beta, gamma)  # scipy.odeintに渡す、微分方程式の追加の引数
    # 計算
    result = odeint(SIR_EQ, ini_state, times, args)

    # 3つに分けてplot
    plt.title("Infection rate transition")
    im1 = plt.plot(times, result[:, 0], label="Susceptible", color="blue")
    im2 = plt.plot(times, result[:, 1], label="Infectious", color="orange")
    im3 = plt.plot(times, result[:, 2], label="Removed", color="green")

    if flag_legend:  # 一回のみ凡例を描画
        plt.legend(loc="best")
        flag_legend = False

    ims.append(im1 + im2 + im3)  # グラフを配列に追加

# 100枚のプロットを 100ms ごとに表示するアニメーション
ani = animation.ArtistAnimation(fig, ims, interval=300)
# plt.show()
HTML(ani.to_jshtml())

MP4またはgifとして保存する

作成したアニメーションはMP4やgif形式で保存ができます。

ani.save('sir-10_100_5-color.mp4', writer="ffmpeg")
ani.save('sir_test-10_100_5-color.gif', writer='imagemagick', fps=20)

フレームレートfpsの値を変えることができます。

ひとこと

作成したmp4ファイルは動画編集ソフトを使えば、いろいろ説明を加えたり、見栄えをよくできます。
これを期に動画を作ってみてはいかがでしょうか。

記事の感想をリアクションでお願いします！

【COVID19×Python】SIRモデルをアニメーションを作って理解する

できあがり

YouTubeにて解説動画を公開中

全体のコードはこちらにて

SIRモデルとは？

Pythonで実装してみよう

流れ

ライブラリの読み込み

パラメータを設定

微分方程式を立てて、解く

グラフをつくる

susceptible, infectious, removedごとにplot

全体をまとめてplot

パラメータを変えて、比較してみる

vulの値を変える

日数を変えてみる

アニメーションをつくる

アニメーション作成のイメージ

MP4またはgifとして保存する

ひとこと

Search

Archive

カテゴリ

記事一覧

Favorites

Twitter

LINK

Favorites

About