あと1か月くらい自粛モードやればなんとかなる?
厚生労働省のホームページに専門家会議の3/19のレポートが載っています。
https://www.mhlw.go.jp/content/10900000/000610566.pdf
これによると日本と状況としては、以下の様な感じの様です。
- 一人の感染者から他の人に感染する人数は2以下。(実効再生産数と呼ぶらしい)
- 累積感染者数1000人弱(検査で陽性判定された人数)
当ブログの前回のエントリーでは実効再生産数を2にした簡易モデル作ってシミュレーションしましたが、総人口を約1.2億として実効再生産数を2より少し下げた形で再シミュレーションしてみました。
rc30-popo.hatenablog.com
実効再生産者数の平均値を1.85、最初の感染者が出た人をDay0として、60日目からは移動抑制策等を取り、実効再生産数を0.74(当初の再生産者数の半分以下)に調整しています。
day0を1月1日と仮置きとすると、3月20日近辺の感染検出者数が1000弱になって、まあまあ現実に近い雰囲気になりました。(Total Detectedの青い縦破線がday 80)
とりあえずday0から120日でシミュレーションしましたが、これだと120日目で完全に収束しきっていません。非専門家の単純シミュレーションですが、後1ヶ月くらいは自粛モードが必要な感じもします。5月の連休明けくらいまで今の抑制策が必要にも見える。(あくまでも素人の私見です。)
# これは日本だけの話で、世界レベルで見るともっとかかる気もします
ちなみに上記の専門家会議のレポートには大規模流行時の10万人あたりの新規感染者数と重篤患者数のシミュレーションを実施したグラフが載っています。
自分のシミュレーションを総人口10万にして、実効再生産者数は2のまま抑制無しとした場合のグラフが下記です。
専門家会議のレポートに乗っているシミュレーション結果が下記です。
ピークのタイミングとかピーク時の1日あたりの感染者数の見積もりが1000人規模でずれていますが、オーダー的には割と近い感じになっています。
以下、最初のシミュレーションのソースです。
# -*- 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()