はじめに
新型コロナウイルス感染症(COVID-19)は、世界中で猛威を振るっており、中国・欧州・米国・アジアなど様々な地域で都市封鎖、いわゆるロックダウンが実行されています。日本においても、4月4日現在、国内での感染者は3360名で、東京都は891名となっており、週末の外出自粛や特定飲食業の回避を知事が要請する事態となっています。また、東京・大阪など都心部だけでなく、全国的な広がりを見せつつあり、まさにオーバーシュートが生じる緊急事態に限りなく近い状況と言えるでしょう。そのような緊急事態を回避するため、政府の専門家会議などでは、３密空間（密閉・密集・密接）を避けるよう、盛んに呼び掛けをしています。
では、このような３密空間をどれだけ避ければいいのでしょうか？本記事では、感染拡大を防ぐために必要な３密空間の避け具合を数値として定量化してみたいと思います。
３密空間と基本再生産数
基本再生産数は、
- 人口集団の典型的な１人の感染者が、その全感染期間において再生産する２次感染者数の期待値
と定義されます。つまり、感染者が他の人に新型コロナウイルスを感染させてしまう強度と言えます。また、その強度が強いと、多くの人に一度に感染させてしまい、クラスターが発生します。厚生労働省が公開している全国クラスターマップによると、クラスターを構成する要素として以下が挙げられています。
- 場所：病院、福祉施設、屋形船、ライブハウス、ジム、パブ、展示会
- 行為：合唱、懇親会、スポーツスクール
これらの共通点が、３密空間なのでしょう。また、クラスター対策班の推計で、基本再生産数の環境別における分析結果が示されています。
この図では、感染源が換気の悪い環境にいた場合とその他の環境にいた場合で、基本再生産数の頻度分布が異なることが示唆されています。つまり、３密空間であるか否かによって、基本再生産数の確率分布が大きく変わるわけです。そこで、環境の違いを表すパラメータを次で定義したいと思います。
- 環境を示すパラメータer: 感染源が非３密空間（良い環境）にいる確率(0-1)
そして、上記のパラメータerによって、基本再生産数を近似的に計算するコードを前回の記事で示しました。条件分岐と一様分布を組み合わせて上の図を近似的に表現しています。
"""
er: good environment ratio (0.-1.)
"""
def COVID19R0(er):
if np.random.rand() < er:
# good environment
if np.random.rand() < 0.8:
R0 = 0
else:
R0 = np.random.randint(1,4)*1.0
else:
# bad environment
R0 = np.random.randint(0,12)*1.0
return R0
では、環境パラメータerがどの程度の値であれば、基本再生産数を１未満（感染を収束させる方向）にできるのでしょうか？以下、統計の道具を使って検討したいと思います。
中心極限定理
中心極限定理とは、
- 観測値$X_1, X_2, \cdots, X_n$が互いに独立な確率変数で全て同じ分布に従う。
- 全てのiについて、$E[X_i]=\mu, V[X_i]=\sigma$。
であるとき、観測値の標本平均$\bar{X}_n$について、
\frac{\bar{X}_n - \mu}{\left( \frac{\sigma}{\sqrt{n}} \right)} \sim N(0,1)\ ({\rm as}\ n \rightarrow \infty) \\
\bar{X}_n = \frac{1}{n} \sum_{i=1}^n X_i
が成り立つことでした。ここで重要なのは、母集団の確率分布については、平均と分散の真値が存在すること以外には仮定がありません。従って、先ほどの基本再生産数の確率分布についてもこの定理を適用することができます。
母平均の信頼区間推定
では、環境パラメータerごとに、基本再生産数の母平均がどうなるかを推定するにはどうすればいいでしょうか。統計の教科書によれば、母集団の分散が未知の場合の、母平均の推定値（信頼区間が$1-\alpha$）は、下記で与えられます。
\left[ \bar{X}_n - \frac{U}{\sqrt{n}} t_{\alpha / 2} (n-1), \bar{X}_n + \frac{U}{\sqrt{n}} t_{\alpha / 2} (n-1) \right]\\
U^2 = \frac{1}{n-1} \sum_{i=1}^n (X_i - \bar{X}_n)^2
ここで、$U$は不偏標本分散、$t_{\alpha }(n-1)$は自由度$n-1$のt分布の$\alpha /2$点（確率分布の裾面積が$\alpha /2$となる確率変数の値）です。教科書を読むと、この信頼区間は、中心極限定理に基づき、母平均が正規分布、母分散が$χ^２$分布に従うことにより推定して、組み合わせた結果導かれるようです。
Pythonで計算してみる
では、Pythonを使って、環境パラメータer依存の基本再生産数R0の母平均と信頼区間を計算してみましょう。
まず、ライブラリをインポートします。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
さきほどの基本再生産数の生成関数です。
"""
er: good environment ratio (0.-1.)
"""
def COVID19R0(er):
if np.random.rand() < er:
# good environment
if np.random.rand() < 0.8:
R0 = 0
else:
R0 = np.random.randint(1,4)*1.0
else:
# bad environment
R0 = np.random.randint(0,12)*1.0
return R0
与えられた確率分布$f$に従って、N個のサンプリングを行い、標本平均と不偏分散を計算する関数です。
def EstFromNsample(f, N):
# 確率分布fにしたがって、N個サンプリング
sample = [f() for i in range(N)]
# 標本平均
Eave = np.average(sample)
# 不偏標本分散
Esig2 = 1./(np.float(N)-1.) * np.sum([pow(s - Eave, 2) for s in sample])
return Eave, Esig2
環境パラメータerを0から1までメッシュで刻んで、それぞれでサンプリングを行い、母平均と母分散の推定値を計算する関数です。
def makeERsamples(keys):
ediv = keys['ediv']
N = keys['N']
#
l = []
#
for e in np.linspace(0, 1, ediv):
f = lambda : COVID19R0(e) # 特定の環境ﾊﾟﾗﾒｰﾀでのｻﾝﾌﾟﾘﾝｸﾞ関数
Eave, Esig2 = EstFromNsample(f, N)
l.append([e, Eave, Esig2])
return np.array(l)
信頼区間を計算する関数です。ここで、alphaは上記の$\alpha$とは意味が逆になっていることに注意です。scipy.statsによれば、"Endpoints of the range that contains alpha percent of the distribution"なので、t_dist.intervalはt分布のalphaパーセントを含む区間の端点を返す関数になっています。
def calcBand01(ave, sigma, N, alpha):
t_dist = stats.t(loc=0,scale=1,df=N-1)
b, u = t_dist.interval(alpha=alpha)
b = ave + b * sigma / np.sqrt(N)
u = ave + u * sigma / np.sqrt(N)
return b, u
母平均推定値とその信頼区間を表示する関数です。
def showBandEstimation(sample, keys):
#
alpha = keys['alpha']
N = keys['N']
#
x = sample[:,0]
y = sample[:,1]
s2 = sample[:,2]
s = np.sqrt(s2)
#
yb, yu = calcBand01(y, s, N, alpha)
#
fig, ax = plt.subplots(1,1,figsize=(5,4), dpi=200)
ax.plot([0,1],[1,1], 'r--')
ax.plot(x, y, label='ave')
ax.plot(x, yu, '.', label='upper')
ax.plot(x, yb, '.', label='lower')
ax.fill_between(x, yu, yb, facecolor='r',alpha=0.2)
ax.set_xlabel('er')
ax.set_ylabel('R0')
ax.set_xticks(np.linspace(0,1,11))
ax.grid(b=True, which='both', axis='x')
ax.legend()
plt.show()
最後に、サンプリングを行って母平均の推定値と信頼区間を表示する部分です。
keys = {'N':100, 'ediv':100, 'alpha':0.95 }
sample = makeERsamples(keys)
fig = showBandEstimation(sample, keys)
計算結果
それでは、いよいよ計算結果を見てみましょう。
サンプリング数(N=100)の場合
サンプリング数(N=1000)の場合
サンプリング数(N=10000)の場合
$N=100$ではかなりばらつきがありますが、$N=10000$では殆ど収束しているように見えますね。ですので、10000サンプルほどあれば相当信頼性の高い推定が出来そうです。
大事なポイントとして、この推定は、１．中心極限定理により標本平均が正規分布へ収束、２．ｔ分布による信頼区間推定、という２段構えで成立しているということです。なので、どの程度の$N$があればこの推定が成り立つのに十分かは、見極めが大事ですね。
考察
以上から、環境パラメータer(３密空間の避け具合）に対する基本再生産数R0の母平均の推定に関して、次の傾向がシミュレーションから導けます。
- 基本再生産数の母平均を１未満にするには、$er > 0.89$ (95%)が必要です。つまり、国民の89%以上が３密空間にいない状態や、１人の人間が１日のうち、89%以上３密空間を避ける状態でないと、感染は収束に向かいません。
- 最も悪い状態（すべての人が３密空間を避けない）のとき、基本再生産数は5.5程度です。
- ちなみに、下のように日本全体の県別実効再生産数ランキングを出してみると、12を超えることもあります。これは、基本再生産数の分散が大きいためと思われます。
さらに言えば・・・
- 国民全体の89%以上が３密空間を避ける必要があるということは、ほぼロックダウンに等しいことかもしれません。しかし、在宅勤務やオンライン教育など、ITを駆使すれば地域別では出来ないこともないでしょう。
- 避けるべきなのはあくまで３密空間なので、そういう空間を徹底的に対策し感染源を隔離したうえで、敢えて経済活動を止めないという日本モデルも成立させることもできるかもしれません。
- 新型コロナとの闘いは、ワクチンや特効薬が普及するまでの長期戦であることを考えると、ロックダウンを継続することは息を止め続けることに等しくいずれ破綻します。経済を殺さず、人を生かすための人工呼吸器のような政策が必要です。
参考リンク
下記のページを参考にさせて頂きました。