人間には未来予知の能力はなく…ても (6)

汎用性皆無だが「今回のパンデミックにおける本日時点までは通用しそうな」話。

ちょっと前に「価値がないだろう」と切り捨てた話ではあるものの、個人的興味としてはあって、やってみた。

csvを自分の都合のいいデータに加工するスクリプト(2)(週平均、二週平均もついでに取ってる)
 1 # -*- coding: utf-8 -*-
 2 import io
 3 import sys
 4 import csv
 5 import datetime
 6 import numpy as np
 7 
 8 
 9 def _tonpa(fn):
10     # rownum=2020/1/16からの日数, colnum=都道府県コード-1
11     result = {
12         "newpat": np.zeros((1, 47)),
13         "newpatmovavg": np.zeros((1, 47)),  # 週平均
14         "newpatmovavg2": np.zeros((1, 47)),  # 2週平均
15         "newdeath": np.zeros((1, 47)),
16         "newdeathmovavg": np.zeros((1, 47)),
17         "newdeathmovavg2": np.zeros((1, 47)),
18         }
19     era = "2020/1/16"
20     era = datetime.datetime(*map(int, era.split("/")))
21     reader = csv.reader(io.open(fn, encoding="utf-8-sig"))
22     next(reader)
23     for row in reader:
24         dt, (p, c, d) = row[0], map(int, row[1::2])
25         dt = (datetime.datetime(*map(int, dt.split("/"))) - era).days
26         if dt > len(result["newpat"]) - 1:
27             for k in result.keys():
28                 result[k] = np.vstack((result[k], np.zeros((1, 47))))
29         result["newpat"][dt, p - 1] = c
30         result["newpatmovavg"][dt, p - 1] = result["newpat"][max(0, dt - 7):dt + 1, p - 1].mean()
31         result["newpatmovavg2"][dt, p - 1] = result["newpat"][max(0, dt - 14):dt + 1, p - 1].mean()
32         result["newdeath"][dt, p - 1] = d
33         result["newdeathmovavg"][dt, p - 1] = result["newdeath"][max(0, dt - 7):dt + 1, p - 1].mean()
34         result["newdeathmovavg2"][dt, p - 1] = result["newdeath"][max(0, dt - 14):dt + 1, p - 1].mean()
35     np.savez(io.open("nhk_news_covid19_prefectures_daily_data.npz", "wb"), **result)
36 
37 
38 if __name__ == '__main__':
39     fn = "nhk_news_covid19_prefectures_daily_data.csv"
40     _tonpa(fn)
2021-02-02追記:
スクリプトをマジメに読んだ人は気付いたかもしれない。「8日平均」と「15日平均」を計算しちゃってる。一連の「主張」の本題とはあまり関係ないとはいえ、さすがにもうちょっと慎重にやるべきでした、すまぬ。正しくは「dt – 6」「dt – 5」…「dt – 1」「dt」の7エレメントでの平均を取らなければならないので、「slice(dt – 6, dt + 1)」。8日平均が役に立たないわけではないけれど、やりたいこととやってることが違うので、間違いは間違い。グラフは滑らかさが若干減って、もうほんの少しだけガタガタになる。
上のスクリプトでセーブしたデータを読み込んで可視化を行うやーつ
 1 # -*- coding: utf-8 -*-
 2 import io
 3 import datetime
 4 import numpy as np
 5 import matplotlib as mpl
 6 import matplotlib.pyplot as plt
 7 import matplotlib.font_manager
 8 
 9 
10 def _vis(dat):
11     era = "2020/1/16"  # 木曜日
12     era = datetime.datetime(*map(int, era.split("/")))
13     npat = dat["newpat"][:,12]  # 12: 東京 (都道府県コード13)
14     npatma = dat["newpatmovavg"][:,12]
15     npatw = np.zeros(npat.shape)
16     for dow in range(7):  # eraが木曜なのでdow=0が木曜
17         wgt = npat[dow::7].sum() / npat.sum()
18         npatw[dow::7] = npatma[dow::7] * wgt * 7
19 
20         
21     fontprop = matplotlib.font_manager.FontProperties(
22         fname="c:/Windows/Fonts/meiryo.ttc")
23 
24     T = [era + datetime.timedelta(t) for t in range(len(npat))]
25     fig, ax1 = plt.subplots(tight_layout=True)
26     fig.set_size_inches(16.53 * 2, 11.69)
27     l1, = ax1.plot(T, npat, alpha=0.8)
28     l2, = ax1.plot(T, npatma, alpha=0.8)
29     l3, = ax1.plot(T, npatw, alpha=0.8)
30     ax1.legend((l1, l2, l3), [
31         "新規陽性者報告数", "新規陽性者報告数週平均", "新規陽性者報告数×曜日重み"],
32             shadow=True, prop=fontprop)
33     ax1.set_xlabel(
34         "東京の新規感染者数推移\n(「NHK新型コロナウイルス関連データ」を加工)",
35         fontproperties=fontprop)
36     fig.autofmt_xdate()
37     #
38     for evtxt, evdt in (
39             ("GoToトラベル全国一時停止", datetime.datetime(2020, 12, 15)),
40             ("1都3県に緊急事態宣言", datetime.datetime(2021, 1, 7)),
41             ("緊急事態宣言11都府県に拡大", datetime.datetime(2021, 1, 13)),
42     ):
43         ax1.annotate("{} ({})".format(evtxt, str(evdt)[:len("????-??-??")]),
44                      xy=(evdt, 0), xycoords='data',
45                      xytext=(
46                          0.5,
47                          (evdt - era).days / (T[-1] - era).days),
48                      textcoords='axes fraction',
49                      arrowprops=dict(facecolor='black', width=0.0001, alpha=0.1),
50                      horizontalalignment='right', verticalalignment='top',
51                      alpha=0.3,
52                      fontproperties=fontprop)
53         ax1.vlines([evdt], 0, 1, transform=ax1.get_xaxis_transform(), alpha=0.1)
54     #
55     #plt.show()
56     fig.savefig("nhk_news_covid19_prefectures_daily_data3.png")
57     plt.close(fig)
58 
59     
60 if __name__ == '__main__':
61     fn = "nhk_news_covid19_prefectures_daily_data.npz"
62     _vis(np.load(io.open(fn, "rb")))

nhk_news_covid19_prefectures_daily_data3

検査数が曜日ごとに違うが、検査数の比が最初の感染者判明の日から本日時点まで変わってなさげ、という感覚を「実際に確認してみた」てことね。検査数ではなく、陽性判明数の比は、上での計算で:

  • 木曜: 0.17
  • 金曜: 0.17
  • 土曜: 0.17
  • 日曜: 0.13
  • 月曜: 0.09
  • 火曜: 0.12
  • 水曜: 0.14

と、かなり有意な差となってあらわれてる。特に日・月曜ね。で、仮に陽性率に曜日ごとの差がないのだとすれば、この比は検査数の比とほとんど同義に違いないだろう、てわけだ。だとするならば、と。この比を重みとして使って、「週平均予測 x 重み」で、その日の感染者が概ね予測できそうだね、てことね。

ざっくりとは思った通りで「割と当たるのかもしらんなぁ」てことではあるね。そして。「で?」て思うことが大事ね。「特定の一日の感染者数予測」が役に立つの、立たないの、という意味では、これは全然役に立たないと言ってよいと思う。大事なのは「月曜が500人で本日木曜が1500人に違いないので今すぐやべー」なんて、こんなの予測でも計画でもなく、合理性のかけらもないんだから。そうじゃなくて、「一週間後には週平均が4000人であり、病床占有期間平均が2週間なので、入院待機人数がこのままだと~人ペースで増え続ける」みたいなことでしょ、「木曜1500人」という特定の日の感染者数そのものは、使い勝手が悪かろ?

そうなんだけど、こうやって実際に可視化してみるのって結構大事かもしれんなぁとは思った。「この比は検査数の比とほとんど同義に違いないだろう」が今でも成立し続けてるんだとしたら、これ、何をあらわしてるかと言えば、「検査体制の抜本的傾向は、今の今まで全然変わってない」てことだと思うのよね。つまり「週末は検査してくれない」という不安が初期にあったのだとして、これが改善したの、してないの、てのは、きっと改善してないんであろうと。


今回の本題とは全然関係ない話もちょっとしとく。

こういう精神的にツラい話ばっかり見聞きしてると心の健康には良くない、とは思うのよね。「病は気から」てのは実際は科学で証明されつつあって、事実なんだけど、心の健康を保つための術ってのも、平時とは違ってきちゃうよなぁとちょっと思っててね。

まず「深呼吸」の効果ってのは、はっきりいって劇的なのよ、本来は。けどさ、「飛沫だぁマイクロ飛沫だぁ」言うてるときに、普段通りの深呼吸って可能なのか、って思っちゃうよねぇ。「換気を十分に行った上で深呼吸しましょう」てことかねぇ?

あと、個人的には今年に入ってからのヒットはのんのんびよりだねぇ。かなり癒される、これは。まぁワタシの場合はアニメだけど、やっぱ何か癒しの手段を見つけておくことは必要だと思う。