新型コロナウイルス感染者数予測をSIRモデルで表現（Python）

$\frac{dy(t)}{dt} = \beta x(t)y(t)-\gamma y(t)$

新型の病気のため、人々は免疫を持っていない
モデルとする人間は、平均的に不特定多数の人と接触している
人口の流入、流出はない（ロックダウンされた国や都市を想定したモデル）
接触によって、一定の確率 $p$ で感染する
感染者は快復後に免疫を獲得し、再び罹患者と接触しても再発しない、未感染者に罹患させることもない。

ここでは、S,I,Rを用いて
S：(Susceptable)未感染者
I : (Infected)感染者
R : (Recovered)快復者
を表すこととする。

総人口を $N$ 、未感染者 $(S)$ が $x$ 人、感染者 $(I)$ が $y$ 人、快復した人 $(R)$ が $z$ 人と置く。また、経過時間を変数 $t$ （日）で表すことにする。集団の中で、一人当たり平均 $m$ 人と接触するとする。そのうち感染している人の割合は $y(t)/N$ であるから、未感染者 $x(t)$ 人が一日あたりに接触する感染者の人数は $m\frac{y(t)}{N} x(t)$ 人と表せる。さらに、未感染者が感染者と接触した際に感染する確率を $p$ とすれば、 $\Delta t$ 日の間に生じる新たな感染者は、 $m\frac{y(t)}{N} x(t)p\Delta t$ となる。

$\beta = \frac{mp}{N}$ と置くと、 $t$ 日目から $t+\Delta t$ 日目の間に、新たな感染者によって非感染者の数が

$x(t+\Delta t)-x(t) = -\beta x(t)y(t)\Delta t$
だけ変化することになる。この式の $\Delta t$ を左辺の分母に移行し、
$\frac{x(t+\Delta t)-x(t)}{\Delta t} = -\beta x(t)y(t)$
と書き直して、 $\Delta t \rightarrow 0$ の極限をとると
$\frac{dx(t)}{dt} = -\beta x(t)y(t)$
を得る。次に、感染者は一日当たり一定の率 $\gamma$ で快復してゆくとすると、回復者の増加分は
$\frac{dz(t)}{dt} = \gamma y(t)$
人口は一定であるため、(非感染者)+(感染者)+(快復者)＝ $N$
$x(t)+y(t)+z(t) = N$
両辺を $t$ で微分すると
$\begin{align} \frac{dy}{dt} &= -\frac{dx(t)}{dt}-\frac{dz(t)}{dt}\\ &= \beta x(t)y(t)-\gamma y(t)\cdots (1) \end{align}$

これをPythonで数値シミュレーションする。4次のルンゲクッタ法を用いて微分方程式からグラフ描写を行う。
$\begin{align} &k_1 \leftarrow f(x, y(x)) \\ &k_2 \leftarrow f(x+h/2, y(x)+k_1 h/2)\\ &k_3 \leftarrow f(x+h/2, y(x)+k_2 h/2)\\ &k_4 \leftarrow f(x+h, y(x)+k_3 h)\\ &y(x+h) \leftarrow y(x)+(k_1 +2k_2+2k_3+k_4)h/6 \end{align}$

import math
import matplotlib.pyplot as plt

N = 1000000
x = N-10
y = 10
z = 0

taxis = []
xaxis = []
yaxis = []
zaxis = []

beta = 0.2/N
gamma = 0.07
dt = 0.001

t = 0
cnt = 0
while t < 180:
    if cnt%100 == 0:
        taxis.append(t)
        xaxis.append(x)
        yaxis.append(y)
        zaxis.append(z)
# step 1
    kx1 = -beta*x*y
    ky1 = beta*x*y - gamma*y
# step 2
    x2 = x+kx1*dt/2
    y2 = y+ky1*dt/2
    kx2 = -beta*x2*y2
    ky2 = beta*x2*y2 - gamma*y2
# step3
    x3 = x+kx2*dt/2
    y3 = y+ky2*dt/2
    kx3 = -beta*x3*y3
    ky3 = beta*x3*y3 - gamma*y3
# step4
    x4 = x+kx3*dt
    y4 = y+ky3*dt
    kx4 = -beta*x4*y4
    ky4 = beta*x4*y4 - gamma*y4

# update
    x = x+(kx1+2*kx2+2*kx3+kx4)*dt/6
    y = y+(ky1+2*ky2+2*ky3+ky4)*dt/6
    z = N-x-y
    t = t+dt
    cnt+=1

plt.title("SIR MODEL")
plt.plot(taxis,xaxis, color=(0.0,1,0.0), linewidth=1.0, label='S')
plt.plot(taxis,yaxis, color=(1.0,0,0.0), linewidth=1.0, label='I')
plt.plot(taxis,zaxis, color=(0.0,0,1.0), linewidth=1.0, label='R')
plt.xlim(0,180)
plt.legend()
plt.xlabel('DAY')
plt.ylabel('X, Y, Z')
plt.grid(True)

f:id:nya__nya:20210418005530p:plain

参考文献

情報基礎「Pythonプログラミング」（ステップ６：SIRモデル）
情報基礎「Pythonプログラミング」（ステップ６：ルンゲクッタ法）

メモ帳

楽しいアウトプットの場所

新型コロナウイルス感染者数予測をSIRモデルで表現（Python）

参考文献