緊急事態宣言の発出は本当は早すぎた？！ 新型コロナウイルス拡大防止の最適タイミングをSIRモデルを使って算出してみた

緊急事態宣言の発出は本当に遅すぎたのか？

厳かに新たな元号を迎えたその翌年、華々しくオリンピックイヤーを迎えるはずだった2020年が始まってすぐに、世界中に不穏なニュースが流れ始めた。
曰く、新種の病気が拡大しているのではないか、と……（画像はイメージです）。

事の起こりは中国の武漢と言われているこの新型コロナウイルスだが、すぐに感染者は世界各国で確認されるようになり、今や世界的大流行に発展してしまった。

とても強い感染力を持つ一方で死亡率は2％程度と言われており、4月14日時点での総死亡者数は世界で約11万人となっている。
年間の死亡者数で言えば、季節性インフルエンザの29万人〜64万人で大体は毎年推移するという数字にも今のところは劣るものであるが、なぜこの新型コロナウイルスがこんなにも騒がれているのか。
理由はいくつかあるようだが、主に以下の点が重要であると思われる。

• インフルエンザには特効薬があるが、コロナにはない
• インフルエンザ以上の死亡率がある
• インフルエンザと違い、潜伏感染する可能性がある（発症前の潜伏期間中も感染力を持つ）

インフルエンザと同程度の死亡率で特効薬も存在しないとくれば、その時点ですでにかなり危険なウイルスであると言える。
その上で、インフルエンザウイルスにはなく新型コロナウイルスにはある特徴として「潜伏感染」が挙げられる。

多くのウイルスは、発症してからでなければ他の個体への感染力が発生しない。
これは結核もインフルエンザも同じで、だから症状が出てから対策を打つのでも決して遅くはない（特効薬があるため）。
しかし、新型コロナウイルスは平均して1週間程度の潜伏期間があり、さらにその潜伏期間においても他の個体への感染力が発生するのだが、これが潜伏感染と呼ばれるものである。
気付かぬうちにウイルスをばらまいていたなんていう事が本当に起こりうるし、おそらくはすでにどこかで現実に発生している（感染経路不明の患者はおそらく全てこれではなかろうか）。

これが新型コロナウイルスの最も恐ろしい特性であると個人的に思う。

日本においては、初の感染者は2020年1月中旬に確認されている。
それから感染者数は今日に至るまで指数関数的に増加し、その勢いは止まる所を知らない。
そしていよいよ4月7日に緊急事態宣言が発出され、主な7都府県に対して外出自粛を求める発表がなされた。
首相官邸HP 新型コロナウイルス感染症対策本部（第２７回）
※ 以下、このHPに従い緊急事態宣言は「発令」ではなく「発出」されたものとする。

そもそもが他国おける宣言やロックダウン（都市封鎖）などとは異なり、この宣言は市民の活動に対する強制力を持たない。
しかし、国民は状況を真摯に捉えて真剣に新型コロナウイルスと向き合おうとしており、各企業は宣言の発出を受けて業務をリモートワークにするなり、自宅待機とするなり出来うる限りの対策を講じている。

同時に、巷でまことしやかにささやかれていることがある。
「緊急事態宣言の発出は遅すぎた」のではないかと。

前述したように、新型コロナウイルスには平均して1週間ほどの潜伏期間があるのだが、長い場合だと2週間ほどになるケースもあるという。
つまり、緊急事態宣言を発出して感染力を低く抑えようとしても、その効果はすぐには現れない。
宣言が発出されてから1週間前後はそれまでと同様か、もしくはそれ以上の新規感染者が出続けることになるだろう。
しかし緊急事態宣言の発出は、早すぎても効果が薄くなってしまうという議論があったりもする。
早すぎてもダメなら遅すぎてもダメ、つまり宣言の発出には絶好のタイミングを測る必要があったという事である。

従って今回は、日本におけるこの絶好のタイミングがいつだったのかに関して、既存のSIRモデルというデータモデルをもとに算出していきたいと思う。
その途中の計算において、「なぜ早すぎるとダメなのか」に対する回答も得られるはずである。

SIRモデルとは

今回の検証に用いるのはSIRモデルというシンプルなモデルである。
参考：SIRモデル -wikipedia-

このモデルは時間の経過に応じて感染者数などがどう推移するかを導出するモデルであり、初めの関数の定義においてこそ簡単な常微分方程式の知識が必要にはなるが、数学的には比較的簡単なモデルであるため扱いやすい。
ただし、SIRモデルはそのままだと潜伏感染者を考慮に入れないので、潜伏感染者の存在を考慮に入れたSEIRモデルを今回は使用する。
参考：SEIRモデル -wikipedia-

SEIRモデルにおいては、全人口を以下のように

• Susceptible：感染前の人
• Exposed：感染した人（潜伏中の人）
• Infected：発症した人
• Recovered：治った人（もしくは亡くなった人）

の4つに分類し、時間の経過に応じてその数がどう推移するかを算出していく。
以下、Wikipediaからのそのままの抜粋。

\frac{dS}{dt}=m(N-S)-βSI　　　(1)

\frac{dE}{dt}=βSI-(m+α)E　　　(2)

\frac{dI}{dt}=αE-(m+γ)I　　　(3)

\frac{dR}{dt}=(m+γ)I　　　(4)

ここで使用されるS, E, I, Rの値は全て「全人口で割った数」なので、初期値においてはSはほぼ1、E, I, Rはほぼゼロと定義した上で計算を開始する。

この式のうちS, E, I, Rの各数を決定する上で重要なパラメータになるのは以下の3つである。

• β：感染率（感染力）
• α：発症率
• γ：回復率

「m：出生率、死亡率」に関しては今回の検証では省略（m=0とする）する（死亡者は回復者Rに含めて計算する）。
また上記に付け加えて「潜伏感染」という要素も式に含めることとする（発症済みの人と潜伏中の人のいずれからも感染するものとする）。

\frac{dS}{dt}=-βS(I+E)　　　(1)

\frac{dE}{dt}=βS(I+E)-αE　　　(2)

\frac{dI}{dt}=αE-γI　　　(3)

\frac{dR}{dt}=γI　　　(4)

「この数式から得られた値」と「実際の日本の感染者数」の推移を比較することで、β, α, γの値をまずは求める。

ちなみに、これ以外にも「完治後の再感染」を考慮に入れるSISモデルというのも存在するが、今回の新型コロナに関しては再感染している兆候は今の所認められないため、再感染は考慮しない。

日本の感染者数グラフの導出

グラフを導出するスクリプトをPythonにて実装する。
基本的な方法としては、β, α, γの値を0から1までの0.01刻みで総当たりでグラフを作成していき、全グラフの中と現実データを最小二乗法で比較した結果、最も現実データとの差異が少ないグラフが今回のベースとなる。
数学的かつPython的にもっとふさわしいやり方があるのかもしれないが、筆者の力不足ゆえ力技になってしまっている点ご容赦いただきたい。

感染者のデータはジョンズ・ホプキンス大学が独自に収集・提供しているデータを使用している（2020年4月14日時点のデータを使用）。
https://github.com/CSSEGISandData/COVID-19

※ 本データは「2020年1月22日」を起点としたデータセットとなっており、以下「〇〇日目」という場合は1月22日を起点とした日数を表すものとする。

また、日本の総人口は2019年度のものを使用する(126,180,643人)。
参考：日本の人口 -wikipedia-

1_analyzing_japanese_infection.py

import matplotlib.pyplot as plt
import numpy as np
import math

# 日本の感染者データ（2020/01/22~）
jp_confirm = np.array([2,2,2,2,4,4,7,7,11,15,20,20,20,22,22,22,25,25,26,26,26,28,28,29,43,59,66,74,84,94,105,122,147,159,170,189,214,228,241,256,274,293,331,360,420,461,502,511,581,639,639,701,773,839,839,878,889,924,963,1007,1101,1128,1193,1307,1387,1468,1693,1866,1866,1953,2178,2495,2617,3139,3139,3654,3906,4257,4667,5530,6005,6748,7370])
jp_death = np.array([0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,2,4,4,5,6,6,6,6,6,6,6,6,10,10,15,16,19,22,22,27,29,29,29,33,35,41,42,43,45,47,49,52,54,54,56,57,62,63,77,77,85,92,93,94,99,99,108,123])
jp_recoverd = np.array([0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,4,9,9,9,9,12,12,12,13,18,18,22,22,22,22,22,22,22,22,32,32,32,43,43,43,46,76,76,76,101,118,118,118,118,118,144,144,144,150,191,232,235,235,285,310,359,372,404,424,424,424,472,472,514,514,514,575,592,622,632,685,762,762,784])
jp_confirm = (jp_confirm - jp_death - jp_recoverd) / 126180643
t = len(jp_confirm)

(ba, bb, bc) = (0, 1, 0.01)
(aa, ab, ac) = (0, 1, 0.01)
(ga, gb, gc) = (0, 1, 0.01)
beta = np.arange(ba, bb, bc)
alfa = np.arange(aa, ab, ac)
gamma = np.arange(ga, gb, gc)

min_n = float('inf')
min_c = '0,0,0'
for b in range(len(beta)):
  for a in range(len(alfa)):
    for g in range(len(gamma)):
      i = np.array([jp_confirm[0]] + [0] * (t - 1), dtype = 'float128')
      e = np.array([jp_confirm[0]] + [0] * (t - 1), dtype = 'float128')
      r = np.array([0] + [0] * (t - 1), dtype = 'float128')
      s = np.array([1 - r[0] - i[0]] + [0] * (t - 1), dtype = 'float128')

      for j in range(t - 1):
        s[j + 1] = s[j] - beta[b] * s[j] * (i[j] + e[j])
        e[j + 1] = e[j] + beta[b] * s[j] * (i[j] + e[j]) - alfa[a] * e[j]
        i[j + 1] = i[j] + alfa[a] * e[j] - gamma[g] * i[j]
        r[j + 1] = r[j] + gamma[g] * i[j]

      estimate = math.sqrt(np.sum((i - jp_confirm) ** 2))
      if min_n > estimate:
        min_n = estimate
        min_c = str(b) + ',' + str(a) + ',' + str(g)

(minb, mina, ming) = map(int, min_c.split(','))

print('β = ' + str(ba + minb * bc))
print('α = ' + str(aa + mina * ac))
print('γ = ' + str(ga + ming * gc))

これを実行すると、以下のようになる（正確に回帰分析を実施すればもっと細かい値になる）。

β = 0.18
α = 0.24
γ = 0.16

このパラメータをもとにSEIRモデルに基づくグラフを作成する。

2-1_japanese_sir_chart.py

import matplotlib.pyplot as plt
import numpy as np
import math

beta = 0.18
alfa = 0.24
gamma = 0.16

population = 126180643

init = 1 / population
t = 400

i = np.array([init] + [0] * (t - 1), dtype = 'float128')
e = np.array([init] + [0] * (t - 1), dtype = 'float128')
r = np.array([0] + [0] * (t - 1), dtype = 'float128')
s = np.array([1 - r[0] - i[0]] + [0] * (t - 1), dtype = 'float128')

for j in range(t - 1):
  s[j + 1] = s[j] - beta * s[j] * (i[j] + e[j])
  e[j + 1] = e[j] + beta * s[j] * (i[j] + e[j]) - alfa * e[j]
  i[j + 1] = i[j] + alfa * e[j] - gamma * i[j]
  r[j + 1] = r[j] + gamma * i[j]

print(r[-1] * 100)

plt.plot(s) 
plt.plot(i + e) 
plt.plot(r) 
plt.savefig("2-1_japanese_sir_chart.png")

青：未感染者
橙：感染者 ＋ 発症者
緑：回復者（発症後に完治もしくは死亡した人、隔離された人）
※ 緑は「回復者」としているが、最終的には全体の感染者の総数となるためその意味で使用することもある

このグラフのうち、X軸が2020年1月22日からの経過日数で、縦軸が人口に占める各人数の割合である。
何も対策を講じなければ、最終的に日本の全人口のうち約76.5%が新型コロナへ感染する計算となる。
2020年4月14日現在(83日目)だと、グラフ上はまだ何も始まってすらいないようにすら見えるが、150日目（2020年6月20日）くらいから感染拡大が始まり、7月下旬にピークを迎えることになる。
自然災害の発生も考慮に入れると、今年の夏はさながら地獄の様相を呈することになるだろう。
当然この予測に関しては、データのグラフへの回帰方法が雑であったり前提が色々と限定的だったりするため、決して正確なものではないと信じたいが、何かしらの目安にはなるはずなので一旦はこの値を信用した上で次の議論へと進む。

また、各算出結果の値の比較を行う場合は、常にこのグラフとの比較であるものとする。

緊急事態宣言の絶好タイミングを算出する

上記で得られたβ, α, γの値をもとに、緊急事態宣言発出を想定した数値シミュレーションを行い、最適なタイミングとなる日付を算出していく。
目標とするのは、以下の図の通り「最も流行のピークを遅くでき」、かつ「ピーク時の最大数を最小化できる」タイミングの算出である。

これにつけ加うるに、「最終的な総感染者数の最小化」も視野に入れて議論を進めていく。

目標

• 流行のピークを最も遅らせる
• ピーク時の最大数を最小化する
• 最終的な総感染者数を最小化する

緊急事態宣言の期間は一旦は30日間で固定とする。
また、感染率の軽減効果に関しては4月14日の感染者数が約半数にまで軽減されたことから、半分にまで軽減できるものと仮定して計算を行う。

計算の方法は単純で、指定された期間の間は緊急事態宣言が発出されているという事で、ある程度の感染率抑制効果が出ているものとみなし、その期間だけ感染率を半分に下げる。

現実の通り、2020年4月7日に宣言発出した場合

まずは、緊急事態宣言が現実の通り2020年4月7日（76日目）に発出された前提のグラフを作成する。

2-2_one_delayed_sir_chart.py

import matplotlib.pyplot as plt
import numpy as np
import math

beta = 0.18
alfa = 0.24
gamma = 0.16

population = 126180643

init = 1 / population
t = 400

i = np.array([init] + [0] * (t - 1), dtype = 'float128')
e = np.array([init] + [0] * (t - 1), dtype = 'float128')
r = np.array([0] + [0] * (t - 1), dtype = 'float128')
s = np.array([1 - r[0] - i[0]] + [0] * (t - 1), dtype = 'float128')

k = 76

for j in range(t - 1):
  beta_t = beta
  if k <= j and j < k + 30:
    beta_t = beta / 2
  s[j + 1] = s[j] - beta_t * s[j] * (i[j] + e[j])
  e[j + 1] = e[j] + beta_t * s[j] * (i[j] + e[j]) - alfa * e[j]
  i[j + 1] = i[j] + alfa * e[j] - gamma * i[j]
  r[j + 1] = r[j] + gamma * i[j]

plt.plot(s) 
plt.plot(i + e) 
plt.plot(r) 
plt.savefig("2-2_one_delayed_sir_chart.png")

約30日ほどピークの時期を遅らせることができるようであるが、ピーク時の最大感染者数の抑制効果はあまり見られない。

最大感染者数（E + I）の、未対策の場合との比較

緊急事態宣言を発出はしたが、どうやら8月中の急激な感染拡大は避けられそうにない（ピークの頂点が8月25日）。
なお今は期間を30日で計算しているが、これの期間を伸ばしてもピークが後ろ倒しになるだけでピーク時の最大感染者数に変化はない。
この場合、経済への打撃の方が深刻になるとも思われるため期間延長は現実的ではないだろう。
参考：期間を90日にした場合の、最大感染者数（I）の未対策の場合との比較

流行のピークを最も遅らせる

流行のピークを最も遅らせることができたであろう発出日を求める。
各発出日ごとにSEIRモデルのグラフを計算し、「最も感染者数が多い日」が一番早く来る発出日を算出することで、この日付を求める。

3-1_delay_max.py

import matplotlib.pyplot as plt
import numpy as np
import math

beta = 0.18
alfa = 0.24
gamma = 0.16

population = 126180643

init = 1 / population
t = 400
dist = 30
emg = [0 for _ in range(t - dist)]

max_delay = 0
delay_k = 0
for k in range(t - dist):
  i = np.array([init] + [0] * (t - 1), dtype = 'float128')
  e = np.array([init] + [0] * (t - 1), dtype = 'float128')
  r = np.array([0] + [0] * (t - 1), dtype = 'float128')
  s = np.array([1 - r[0] - i[0]] + [0] * (t - 1), dtype = 'float128')

  tmp_max = 0
  max_t = 0
  for j in range(t - 1):
    beta_t = beta
    if k <= j and j < k + dist:
      beta_t = beta / 2
    s[j + 1] = s[j] - beta_t * s[j] * (i[j] + e[j])
    e[j + 1] = e[j] + beta_t * s[j] * (i[j] + e[j]) - alfa * e[j]
    i[j + 1] = i[j] + alfa * e[j] - gamma * i[j]
    r[j + 1] = r[j] + gamma * i[j]

    if tmp_max < i[j + 1]:
      tmp_max = i[j + 1]
      max_t = j # 感染者数がピークの時点

  if max_delay < max_t:
    max_delay = max_t
    delay_k = k
  emg[k] = max_t
print(delay_k)
plt.plot(emg)
plt.savefig("3-1_delay_max.png")

結果

158

これは、発出日ごとの「最も感染者数が多い日」の分布である。
縦軸が上であるほど、ピークの日付を遅らせることができていることを表す。

これをもとにSEIRモデルのグラフを作成すると

山が二つになる。
つまり仮に158日目（2020年6月28日）に宣言を発出していたとしたら、ピークを遅らせるというよりはピークを2回に分散させることができていたことになる。

最大感染者数（E + I）の、未対策の場合との比較

ピーク時の最大数を最小化する

ピーク時の最大感染者数を最小化することができたであろう発出日を求める。
各発出日ごとにSEIRモデルのグラフを計算し、「最も感染者数が多い日」の感染者数が一番少なくなる発出日を算出することで、この日付を求める。

3-2_minimize_moment_infection.py

import matplotlib.pyplot as plt
import numpy as np
import math

beta = 0.18
alfa = 0.24
gamma = 0.16

population = 126180643

init = 1 / population
t = 400
dist = 30
emg = [0 for _ in range(t - dist)]

min_max = float('inf')
min_c = 0
for k in range(t - dist):
  i = np.array([init] + [0] * (t - 1), dtype = 'float128')
  e = np.array([init] + [0] * (t - 1), dtype = 'float128')
  r = np.array([0] + [0] * (t - 1), dtype = 'float128')
  s = np.array([1 - r[0] - i[0]] + [0] * (t - 1), dtype = 'float128')

  tmp_max = 0
  for j in range(t - 1):
    beta_t = beta
    if k <= j and j < k + dist:
      beta_t = beta / 2
    s[j + 1] = s[j] - beta_t * s[j] * (i[j] + e[j])
    e[j + 1] = e[j] + beta_t * s[j] * (i[j] + e[j]) - alfa * e[j]
    i[j + 1] = i[j] + alfa * e[j] - gamma * i[j]
    r[j + 1] = r[j] + gamma * i[j]

    tmp_max = max(tmp_max, i[j + 1])

  if min_max > tmp_max:
    min_max = tmp_max
    min_c = k
  emg[k] = tmp_max
print(min_max * 100)
print(min_c)
plt.plot(emg)
plt.savefig("3-2_minimize_moment_infection.png")

結果

4.455430750546662306
159

これは、発出日ごとの「ピーク時の最大数」の推移を表す。
日付自体は、ピークを最も遅らせるケースとほぼ同じとなったため、グラフの作成と比較は省略する。

総感染者数を最小化する

最終的な総感染者数を最小化することができたであろう発出日を求める。
各発出日ごとにSEIRモデルのグラフを計算し、最終的な総感染者数が最も少なくなる発出日を算出することで、この日付を求める。

3-3_minimize_total_infection.py

import matplotlib.pyplot as plt
import numpy as np
import math

beta = 0.18
alfa = 0.24
gamma = 0.16

population = 126180643

init = 1 / population
t = 400
dist = 30
emg = [0 for _ in range(t - dist)]

min_n = float('inf')
min_c = 0
for k in range(t - dist):
  i = np.array([init] + [0] * (t - 1), dtype = 'float128')
  e = np.array([init] + [0] * (t - 1), dtype = 'float128')
  r = np.array([0] + [0] * (t - 1), dtype = 'float128')
  s = np.array([1 - r[0] - i[0]] + [0] * (t - 1), dtype = 'float128')

  for j in range(t - 1):
    beta_t = beta
    if k <= j and j < k + dist: # ← この期間内は感染率が半分になるものとする
      beta_t = beta / 2
    s[j + 1] = s[j] - beta_t * s[j] * (i[j] + e[j])
    e[j + 1] = e[j] + beta_t * s[j] * (i[j] + e[j]) - alfa * e[j]
    i[j + 1] = i[j] + alfa * e[j] - gamma * i[j]
    r[j + 1] = r[j] + gamma * i[j]

  if min_n > r[-1]:
    min_n = r[-1]
    min_c = k
  emg[k] = r[-1]
print(min_n * 100)
print(min_c)
plt.plot(emg)
plt.savefig("img2.png")

グラフにすると以下のようになる。

このグラフは、発出日ごとの総感染者数の推移を表す。

この大きく落ちくぼんでいる個所の値は

59.260485603852678493
173

173日目、2020年7月13日となる。
これを、SEIRモデルのグラフに当てはめると以下のようになる。

総感染者数（R）の、未対策の場合との比較

減らすことはできているが、それでも総感染者数は数割程度しか減らなかった。

計算結果総括

ピークを遅らせるという目的とピーク時の最大数を減らすという目的に関しては、いずれも「山を二つに分散させる」ことで効果は最大化される。
一方で、ピーク直前に宣言を発出して「残りの未感染者への感染を防ぐ」ことで、最終的な総感染者数を減らす効果を最大化することができる。

いずれにせよ、データから導き出したグラフをもとに算出した絶好のタイミングは、現実よりもはるかに後ろであったことがわかる。
どの効果を狙うにしても、少なくともある程度の感染者がすでに存在している段階でなければ軽減の効果は薄い。
すなわち、完全に計算上の妙によって今回の結果ははじき出されたものではあるが、逆に言えば計算上は緊急事態宣言の発出は決して遅くないどころか逆に早すぎたということになる。

ならばどうするか、緊急事態宣言をもう一回発出してしまえば良いのである。

麻疹の予防接種も2度打つではないか、それと同じことである。

再発出の絶好タイミングを算出する

4月7日から5月6日までの緊急事態宣言はこのまま固定で置くとして、そのあと再宣言を発出するタイミングを再度、2通り求める。

• ピーク時の最大数を抑制する（≒ ピークを遅らせる）
• 総感染者数を抑制する

ピークの最大数を抑制する

4-1_redeclaration_momemt_minimize.py

import matplotlib.pyplot as plt
import numpy as np
import math

beta = 0.18
alfa = 0.24
gamma = 0.16

population = 126180643

init = 1 / population
t = 400
dist = 30
emg = [0 for _ in range(t - dist)]

min_max = float('inf')
min_c = 0
for k in range(t - dist):
  i = np.array([init] + [0] * (t - 1), dtype = 'float128')
  e = np.array([init] + [0] * (t - 1), dtype = 'float128')
  r = np.array([0] + [0] * (t - 1), dtype = 'float128')
  s = np.array([1 - r[0] - i[0]] + [0] * (t - 1), dtype = 'float128')

  tmp_max = 0
  for j in range(t - 1):
    beta_t = beta
    if (76 <= j and j < 106) or (k <= j and j < k + dist):
      beta_t = beta / 2
    s[j + 1] = s[j] - beta_t * s[j] * (i[j] + e[j])
    e[j + 1] = e[j] + beta_t * s[j] * (i[j] + e[j]) - alfa * e[j]
    i[j + 1] = i[j] + alfa * e[j] - gamma * i[j]
    r[j + 1] = r[j] + gamma * i[j]

    tmp_max = max(tmp_max, i[j + 1])

  if min_max > tmp_max:
    min_max = tmp_max
    min_c = k
  emg[k] = tmp_max
print(min_max * 100)
print(min_c)
plt.plot(emg)
plt.savefig("4-1_redeclaration_momemt_minimize.png")

実行すると以下のようになり、191日目（2020年07月31日）に再宣言をすることでピーク時の最大感染者数はかなり抑えられるようになる。

4.3781021262811077857
191

これをSEIRモデルに当てはめると

最大感染者数（E + I）の、未対策の場合との比較

やはり最大数は半分以下に抑えることができる。

総感染者数を抑制する

4-2_redeclaration_total_minimize.py

import matplotlib.pyplot as plt
import numpy as np
import math

beta = 0.18
alfa = 0.24
gamma = 0.16

population = 126180643

init = 1 / population
t = 400
dist = 30
emg = [0 for _ in range(113, t - dist)]

min_n = float('inf')
min_c = 0
for k in range(t - dist):
  i = np.array([init] + [0] * (t - 1), dtype = 'float128')
  e = np.array([init] + [0] * (t - 1), dtype = 'float128')
  r = np.array([0] + [0] * (t - 1), dtype = 'float128')
  s = np.array([1 - r[0] - i[0]] + [0] * (t - 1), dtype = 'float128')

  for j in range(t - 1):
    beta_t = beta
    if (83 <= j and j < 113) or (k <= j and j < k + dist):
      beta_t = beta / 2
    s[j + 1] = s[j] - beta_t * s[j] * (i[j] + e[j])
    e[j + 1] = e[j] + beta_t * s[j] * (i[j] + e[j]) - alfa * e[j]
    i[j + 1] = i[j] + alfa * e[j] - gamma * i[j]
    r[j + 1] = r[j] + gamma * i[j]

  if min_n > r[-1]:
    min_n = r[-1]
    min_c = k
  emg[k - 113] = r[-1]
print(min_n * 100)
print(min_c)
plt.plot(emg)
plt.savefig("img2.png")

実行すると以下のようになり、205日目（2020年08月14日）に再宣言をすることで総感染者数はかなり抑えられる。

59.185455818817687235
205

SEIRモデルに当てはめると

総感染者数（R）の、未対策の場合との比較

相当抑制できる。

結論

結論としては、感染率を減らすための対策は今回だけではなく二段構えで実施するのがよく、今このタイミングで緊急事態宣言が発出されているのは「ワクチンの開発が完了するまでの時間稼ぎ」という意味で考えれば決して悪くはないし、タイミングとしては遅くなくむしろ全然早い。
そして緊急事態宣言の再発出をする場合は、タイミングとしては8月上旬から中旬にかけてが望ましいという結果となった。

当然、それまでの間にワクチンなどが開発され感染拡大が収束すれば、宣言の再発出は不要となる。

所感

一番初めに導出したSEIRのグラフが予想以上にそれっぽくて、その時点ですでにかなりビビった（その前にSIRで試行錯誤しまくって全然うまくいっていなかったので心は半分以上折れていた）。

この記事における計算結果はあくまでも机上の計算であるので（あえて空論とは言わない）、現実通りであるかどうかは別の話である（特に全体的な人数に関して）。
また、今回の計算においては「緊急事態宣言を発出することにより、感染率が半分になる」という前提に基づいているということもあり、正直この前提もかなり疑わしいと言える。
ベースとなる関数の回帰も力技でガバガバである。
そのため、今回の結論は一つの目安として留めておいて頂けると幸いである。

なおこの分析は、これまでデータサイエンスとはあまり縁のなかった筆者が、自宅待機となって暇を持て余しすぎたがため行うに至ったものであること、ご承知置き頂きたく。

今回の分析に使用したソースは以下のリポジトリに入っているので、自由に参照されたし。
https://github.com/mashinosatoshi/corona_japanese_infection_analyzation