ディスクレーマー
私は感染学、免疫学の専門家でもなければ、セルオートマトンの専門家でもないため、用語遣いとか、前提の数字とかはめちゃめちゃ適当です。ご容赦ください。
あくまでも、数理モデルとpythonプログラミングに関する資料の一環としてご覧ください。
モチベーション
**接触８割減っていうけど、７割減だといけないの？**とか、**収束まであとどれくらいかかるの？という声が散見される今日この頃。
もちろん厳密な数値はPCR検査等で毎日確認されていますが、ざっくりした感染者数の広がり方を見積もりたければ、数理モデルの出番でしょう。
というわけで、かんたんなおもちゃモデル（toy model）**を作り、所与の人口密度、濃厚接触時に感染する確率、治癒率に対する、集団全体の感染者数の推移をシミュレーションしてみました。
※あくまでおもちゃモデルですので、これをもって政策がどうとか議論するのはナンセンスです。しつこいようですが、あくまで知的好奇心目的でご覧ください。
数理モデル
ここでは簡単のために、以下の仮定を置きます。
１. 海も山もない平らな2次元世界の住民とする（要はセルオートマトン）。
２. 全員徒歩移動で、新幹線や飛行機は使わないものとする（1時点分の移動距離は最大1で、斜め移動もあり）。
３. 濃厚接触（セルが隣接）した人は、20%の確率で次の日に感染するものとする。
言い換えると、濃厚接触しているクラスター内の感染者が4日で約2倍に増えるとする仮定です。
(1+0.20)^4 = 2.07... \approx 2
４. 感染している人は、2%の確率で次の日に治癒するものとする。
言い換えると、感染者の約半分が1ヶ月（31日）で治癒するとする仮定です。
1-(1-0.02)^{31} = 0.46... \approx 0.5
実際に作ってみた。
時刻 t=0（日）での感染者数は全体の5%（100人のうち5人が初めに感染している）とします。
左の図が開始から1年（365日）後の集団の様子を示しており、青が健康な人、赤が感染者で、全体的に真っ赤だとパンデミックという感じでしょうか。
また、右の図が感染率の時系列推移を示しています。横軸が時間、縦軸が感染率です。要はこれ↓です。
人口密度＝１０％のケース（大都会）
100個のセルに10人いるケースです。図を見ていただければわかりますが大都会のイメージですね。
1年経過時の集団は真っ赤っかのパンデミック状態になってしまいました。
自然治癒のロジックがあるため感染率100%には至らず、90%くらいで飽和していますが、これではとても収束が見込めません。
人口密度＝４％（接触６割減）のケース
ずいぶん自粛してると思うじゃないですか？とんでもない！
こんなもんじゃ全然足りないのです。
1年経過時の感染率は50%近くまで低下していますが、収束まではまだまだ遠い。
この世界では1年先延ばしにしたところでオリンピックが開催出来ません。
人口密度＝３％（接触７割減）のケース
8割じゃなくて7割じゃないとダメなんですか？
ええ、7割じゃダメなんでしょうね。
たしかに感染率はかなり抑えられていますが、ゆっくりではあるものの1年近く増加を続けています（終盤でややピークアウト）。
感染スピードは抑制できているものの、自然収束に至るまでに1年では足らないようです。
人口密度＝２％（接触８割減）のケース
新規感染者数と治癒者数がほとんど均衡しました。
新規感染率も治癒率も適当に決めた割に、ずいぶん綺麗な結果になりましたね。
まったく素人ですが、接触8割減という数字は、感染学において一種のrule of thumbなんでしょうかね？
人口密度＝１％（接触９割減）のケース
200日ちょっとで感染者数が0になりました。
現実世界も、はやくこうなってほしいものですね～
pythonコード
こちらリアルタイムで推移するアニメーションを描画する仕様となっております。
セルの値が0のとき誰もいない、-1は健康な人、+1は感染者であることを示しています。
import numpy as np
from numpy import random
import seaborn as sns
import matplotlib.pyplot as plt
sns.set()
# parameter
density = 0.03
prob_infect = 0.20
prob_heal = 0.02
infected_initial = 0.05
# initialize earth
earth = np.where((random.random(10000) > (1 - density)), -1, 0).reshape(100, 100)
for i in range(int(10000 * density * infected_initial)):
earth[random.randint(0, 100), random.randint(0, 100)] = 1
time = 0
infection_ratio = np.zeros(10 ** 6) * np.nan
infection_ratio[time] = 0
plt.close()
fig, ax = plt.subplots(figsize=(12, 6), nrows=1, ncols=2)
no_end = True
while no_end:
for (x, y) in [(x, y) for x in range(100) for y in range(100)]:
if earth[x, y] != 0:
# infection
(x_min, x_max) = (max(x - 1, 0), min(x + 1, 98) + 1)
(y_min, y_max) = (max(y - 1, 0), min(y + 1, 98) + 1)
if np.any(earth[x_min: x_max, y_min: y_max] == 1):
earth[x, y] = 1 if random.rand() < prob_infect else earth[x, y]
# healing
earth[x, y] = -1 if random.rand() < prob_heal else earth[x, y]
# random walk around
x_next = min(max(x + random.randint(-1, 2), 0), 99)
y_next = min(max(y + random.randint(-1, 2), 0), 99)
if ((x_next, y_next) != (x, y)) & (earth[x_next, y_next] == 0):
earth[x_next, y_next] = earth[x, y]
earth[x, y] = 0
time += 1
infection_ratio[time] = np.sum(np.where(earth == 1, 1, 0)) / (10000 * density) * 100
try:
subplot0.cla()
subplot1.cla()
except:
print('no need to refresh figure')
plt.suptitle('time: {}'.format(time))
subplot0 = sns.heatmap(earth, vmin=-1, vmax=1, center=0, cmap='coolwarm', ax=ax[0], cbar=False)
ax[0].axis('off')
ax[0].set_title('density: {}%'.format(int(density * 100)))
subplot1 = sns.lineplot(data=infection_ratio)
ax[1].set_title('time series of infection')
ax[1].set_xlabel('time')
ax[1].set_ylabel('infection ratio %')
ax[1].set_ylim([0, 100])
ax[1].set_xlim([0, time + 10])
plt.tight_layout()
plt.draw()
plt.pause(0.01)
