素人が新型コロナ感染のシミュレーションをやってみる(日本編: python+numpy)

あと1か月くらい自粛モードやればなんとかなる？

厚生労働省のホームページに専門家会議の3/19のレポートが載っています。
https://www.mhlw.go.jp/content/10900000/000610566.pdf

これによると日本と状況としては、以下の様な感じの様です。

• 一人の感染者から他の人に感染する人数は2以下。(実効再生産数と呼ぶらしい)
  
• 累積感染者数1000人弱(検査で陽性判定された人数)

当ブログの前回のエントリーでは実効再生産数を2にした簡易モデル作ってシミュレーションしましたが、総人口を約1.2億として実効再生産数を2より少し下げた形で再シミュレーションしてみました。
rc30-popo.hatenablog.com

日本での感染シミュレーション(60日目から自粛モード)

実効再生産者数の平均値を1.85、最初の感染者が出た人をDay0として、60日目からは移動抑制策等を取り、実効再生産数を0.74(当初の再生産者数の半分以下)に調整しています。
day0を1月1日と仮置きとすると、3月20日近辺の感染検出者数が1000弱になって、まあまあ現実に近い雰囲気になりました。(Total Detectedの青い縦破線がday 80)

とりあえずday0から120日でシミュレーションしましたが、これだと120日目で完全に収束しきっていません。非専門家の単純シミュレーションですが、後1ヶ月くらいは自粛モードが必要な感じもします。5月の連休明けくらいまで今の抑制策が必要にも見える。(あくまでも素人の私見です。)
# これは日本だけの話で、世界レベルで見るともっとかかる気もします

ちなみに上記の専門家会議のレポートには大規模流行時の10万人あたりの新規感染者数と重篤患者数のシミュレーションを実施したグラフが載っています。

自分のシミュレーションを総人口10万にして、実効再生産者数は2のまま抑制無しとした場合のグラフが下記です。

10万人の大規模流行シミュレーション(私家版)

専門家会議のレポートに乗っているシミュレーション結果が下記です。
ピークのタイミングとかピーク時の1日あたりの感染者数の見積もりが1000人規模でずれていますが、オーダー的には割と近い感じになっています。

2020/3/19 専門家会議レポート抜粋

以下、最初のシミュレーションのソースです。

# -*- Coding: utf-8 -*-

# numpy
import numpy as np

# Matplotlib
import matplotlib as mpl
mpl.use('TkAgg')  # or whatever other backend that you want
import matplotlib.pyplot as plt

# シミュレーション日数
SIM_DAYS = 120
# 感染後、陰性になるまでの平均日数
AVG_RECOVERY_DAYS = 10
# 重症となる割合
SEVERE_PATIENT_RATE = 0.20
# 重症化するまでの平均日数
AVG_SEVERE_DAYS = 7
# 一人の陽性者が1日あたりに感染させる平均人数
AVG_INFECTION_RATE_PER_PERSON = 0.185
# 世界人口
#WORLD_POPULATION = 7700000000
WORLD_POPULATION = 120000000

# シミュレーション停止レシオ(対世界人口)
SIM_STOP_RATIO=0.9

# 移動制限等の対策パラメータ
# 対策有無
CONTROL_MEASURE=True
# 対策を開始する日
CONTROL_START_DAY=60
# 対策による感染率低減
CONTROL_FACTOR = 0.4

# グラフ上の現在日目安
CURRENT_DAY=80

def controlFactor(day):
    global CONTROL_MEASURE
    global CONTROL_START_DAY
    global CONTROL_FACTOR
    retvalue = 1.0
    if CONTROL_MEASURE:
        if day >= CONTROL_START_DAY:
            retvalue = CONTROL_FACTOR
    return retvalue

# 日別の新規感染者
new_infected = np.zeros(SIM_DAYS+1)
new_infected[0] = 1
# 日別の累積感染者
total_infected = np.zeros(SIM_DAYS+1)
total_infected[0] = 1
# 重症化して陽性判定される患者数
new_detected = np.zeros(SIM_DAYS+1)
# 重症化して陽性判定される患者数の累積
total_detected = np.zeros(SIM_DAYS+1)



sim_stop_day = SIM_DAYS
# シミュレーション
for day in range(1,SIM_DAYS + 1):
    # その日新規に感染する人数
    day_idx = (day - AVG_RECOVERY_DAYS) if day - AVG_RECOVERY_DAYS >= 0 else 0
    ctrl_factor = controlFactor(day)
    infected_cnt = np.sum(new_infected[day_idx:day] * AVG_INFECTION_RATE_PER_PERSON * ctrl_factor * (WORLD_POPULATION - total_infected[day - 1])/WORLD_POPULATION)
    new_infected[day] = infected_cnt
    total_infected[day] = total_infected[day - 1] + infected_cnt

    day_idx2 = day -  AVG_SEVERE_DAYS
    if day_idx2 >= 0:
        new_detected[day] = new_infected[day_idx2] * SEVERE_PATIENT_RATE
    total_detected[day] = total_detected[day - 1] + new_detected[day]

    print('day = %d, new_infected = %f, new_detected = %f,total_infected = %f,total_detected = %f' % (day,infected_cnt,new_detected[day],total_infected[day],total_detected[day] ))

    if total_infected[day] >= WORLD_POPULATION * SIM_STOP_RATIO: # 感染者が世界人口x停止レシオ超えたら停止
        print('Total infected count exceeds World population. x %f' % SIM_STOP_RATIO)
        sim_stop_day = day
        break



# グラフ描画
days = np.arange(sim_stop_day + 1)
fig = plt.figure()
ax1 = fig.add_subplot(221)
ax2 = fig.add_subplot(222)
ax3 = fig.add_subplot(223)
ax4 = fig.add_subplot(224)

ax1.plot(days[0:sim_stop_day + 1],total_infected[0:sim_stop_day + 1],label='Total Infected')
#ax1.hlines([WORLD_POPULATION],0,sim_stop_day,"red",linestyles='dashed')
ax2.bar(days[0:sim_stop_day + 1],new_infected[0:sim_stop_day + 1],label='Newly Infected')
ax3.plot(days[0:sim_stop_day + 1],total_detected[0:sim_stop_day + 1],label='Total Detected')
ax3.vlines(CURRENT_DAY,0,max(total_detected[0:sim_stop_day + 1]),"blue",linestyles='dashed')
ax4.bar(days[0:sim_stop_day + 1],new_detected[0:sim_stop_day + 1],label='Newly Detected')
ax4.vlines(CURRENT_DAY,0,max(new_detected[0:sim_stop_day + 1]),"blue",linestyles='dashed')

ax1.legend()
ax2.legend()
ax3.legend()
ax4.legend()

plt.show()