人口比 (3?)

これの続きだと果たして言えるのやらどうやら。

これ、それにこれでの主張というのは、「質の良い見える化を皆で共有すべし、それをするのはメディアの責務だ、それによってコンセンサスを醸成せよ」「二週間前の行動の反映なのだから、二週間後予測に基づいて行動するのがデフォルトなんではないんかーい」の二点が主たるものである。さらに「人口比」という見出しで括って言いたかったのは「相応しい指標で思考すべきで、そのための「統一的なものの見方」が必要である」てこと。

後者の、「相応しい指標」は、ほんとは専門家会議や各地方自治体行政がきちんと考えて、報告してくれていたりはするわけである。むろん「個々が個々に」なのはマズいにせよ、各々ちゃんとやろうとはしているわけだ。でもメディアは「批判が仕事である」と勘違いしているので、「指標を示せ」と行政や政府に文句だけを言い続ける。いや…、「お前らにも出来るんだし、出てるだろが」と。そういうことをワタシは言いたいわけである。ワタシが「人口比」で主張したことは、東京都と政府の指標にちゃんと含まれてる。それをみてる。メディアがそこを軽視して「今日は新たに300人でしたきゃー」と、そこばかりを強調し、もっと議論しやすい移動平均で議論することを怠り、10万人あたりを伝えることも軽視し、なので「島根で10人感染者が出ることがどれほど危険なことなのか」が一向に国民に伝わらない。

まぁそんなことが、言いたかったことなわけなんだけれど、もう一つ、隠れたテーマとしては、「こうした可視化は、20年前と今とでは、簡単さが雲泥の差であるということを伝えること」でもあったりした。正直ワタシが学生の頃なんてのは、今では誰でも使える「スプライン補完」みたいな基礎的な平滑化ですら、インフラが今ほど整っていなかったので、まぁ言ってしまえば「大学のような恵まれた場所でのみ簡単」というシロモノだったわけである。けれども、今は、はっきりいって「パソコンをお持ちの全国民」が、少しの学習と若干の手間さえかければ、(基礎となる理論等々への理解の有無によらず)実施できるようになった。そう、「誰でも」だ。老婆だろうが小学生だろうが、誰でもだ。この手段は「ワタシの場合は Python + Matplotlib」ということにはなるんだけれど、この道具に関しては色々選べる。MS Excel で出来ることも増えたし、フリーのものが良ければ、LibreOffice もかなり強力なことが出来る。(ちなみに 8割オジサンでお馴染みの西浦先生は R 使いらしい。)

というわけで、「言いたかった主たること」はもう大方言い尽くしたのではあるけれど、もう少しだけ「見える化の実現、のサンプル」を続けようかと。今回は感染拡大のステージングを素直に可視化しようとするとどうなるんだろね、てやつ。ただ、一般人が簡単にアクセス出来るデータで出来ることは限られてるので、出来ることだけを。

まずは、改めて(少し拡充したりもしてるので)データ取得:

csvを自分の都合のいいデータに加工するスクリプト(NHKまとめのもの)→今回は本質的には使わないんだけど、依存はしてるので一応。
 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         "effrepro": np.zeros((1, 47)),  # 実効再生産数
19         }
20     era = "2020/1/16"
21     era = datetime.date(*map(int, era.split("/")))
22     reader = csv.reader(io.open(fn, encoding="utf-8-sig"))
23     next(reader)
24     for row in reader:
25         dt, (p, c, d) = row[0], map(int, row[1::2])
26         dt = (datetime.date(*map(int, dt.split("/"))) - era).days
27         if dt > len(result["newpat"]) - 1:
28             for k in result.keys():
29                 result[k] = np.vstack((result[k], np.zeros((1, 47))))
30         result["newpat"][dt, p - 1] = c
31         result["newpatmovavg"][dt, p - 1] = result[
32             "newpat"][max(0, dt-6):dt+1, p - 1].mean()
33         result["newpatmovavg2"][dt, p - 1] = result[
34             "newpat"][max(0, dt-13):dt+1, p - 1].mean()
35         result["newdeath"][dt, p - 1] = d
36         result["newdeathmovavg"][dt, p - 1] = result[
37             "newdeath"][max(0, dt-6):dt+1, p - 1].mean()
38         result["newdeathmovavg2"][dt, p - 1] = result[
39             "newdeath"][max(0, dt-13):dt+1, p - 1].mean()
40     # 実効再生産数は、現京都大学西浦教授提案による簡易計算(を一般向け説明にまとめてくれ
41     # ている東洋経済ONLINE)より:
42     #   (直近7日間の新規陽性者数/その前7日間の新規陽性者数)^(平均世代時間/報告間隔)
43     #   平均世代時間は=5日、報告間隔=7日の仮定。
44     # のつもりだが、東洋経済ONLINEで GitHub にて公開(閉鎖予定)されてる計算結果の csv とは
45     # 値が違う。違うことそのものは致し方ないこと。計算タイミングによって結果は変わるし、
46     # NHK まとめと東洋経済ONLINEのデータはシンクロしてないから。思いっきり間違ってるって
47     # ことはなさそう。近い値にはなってる。
48     for dt in range(len(result["newpat"])):
49         # 実効再生産数(都道府県別)
50         sum1 = result["newpat"][max(0, dt-6):dt+1, :].sum(axis=0)
51         sum2 = result["newpat"][max(0, dt-6-7):max(0, dt-6)+1, :].sum(axis=0)
52         result["effrepro"][dt, :] = np.power(sum1 / sum2, 5. / 7.)
53     np.savez(io.open("nhk_news_covid19_prefectures_daily_data.npz", "wb"), **result)
54 
55 
56 if __name__ == '__main__':
57     fn = "nhk_news_covid19_prefectures_daily_data.csv"
58     _tonpa(fn)
csvを自分の都合のいいデータに加工するスクリプト(東洋経済オンラインのもの)
  1 # -*- coding: utf-8-unix -*-
  2 import io
  3 import csv
  4 import datetime
  5 import numpy as np
  6 
  7 
  8 _PREFCODE = {
  9     '北海道': 1,
 10     '青森県': 2,
 11     '岩手県': 3,
 12     '宮城県': 4,
 13     '秋田県': 5,
 14     '山形県': 6,
 15     '福島県': 7,
 16     '茨城県': 8,
 17     '栃木県': 9,
 18     '群馬県': 10,
 19     '埼玉県': 11,
 20     '千葉県': 12,
 21     '東京都': 13,
 22     '神奈川県': 14,
 23     '新潟県': 15,
 24     '富山県': 16,
 25     '石川県': 17,
 26     '福井県': 18,
 27     '山梨県': 19,
 28     '長野県': 20,
 29     '岐阜県': 21,
 30     '静岡県': 22,
 31     '愛知県': 23,
 32     '三重県': 24,
 33     '滋賀県': 25,
 34     '京都府': 26,
 35     '大阪府': 27,
 36     '兵庫県': 28,
 37     '奈良県': 29,
 38     '和歌山県': 30,
 39     '鳥取県': 31,
 40     '島根県': 32,
 41     '岡山県': 33,
 42     '広島県': 34,
 43     '山口県': 35,
 44     '徳島県': 36,
 45     '香川県': 37,
 46     '愛媛県': 38,
 47     '高知県': 39,
 48     '福岡県': 40,
 49     '佐賀県': 41,
 50     '長崎県': 42,
 51     '熊本県': 43,
 52     '大分県': 44,
 53     '宮崎県': 45,
 54     '鹿児島県': 46,
 55     '沖縄県': 47,
 56 }
 57 _COLSNM = (
 58     "testedPositive",
 59     "peopleTested",
 60     "hospitalized",
 61     "serious",
 62     "discharged",
 63     "deaths",
 64     "effectiveReproductionNumber",
 65 )
 66 def _tonpa(fn):
 67     # 元データ prefectures.csv の著作権:
 68     # 『東洋経済オンライン「新型コロナウイルス 国内感染の状況」制作:荻原和樹』
 69     #
 70     # 1. 開始日が NHK まとめと違う
 71     # 2. 更新のない都道府県は省略してある
 72     # 3. 日付は一応抜けなく連続している、今のところ
 73     #
 74     # NHKまとめのほうで密行列として扱った同じノリにしたい、ここでは面倒だけど。
 75     # つまり、列方向に都道府県、行方向に日付。
 76     result = {
 77         k: np.zeros((1, 47))
 78         for k in _COLSNM}
 79     reader = csv.reader(io.open(fn, encoding="utf-8-sig"))
 80     next(reader)
 81     def _f(d):
 82         if not d or d == "-":
 83             # 不明とか未報告の意味だと思う。
 84             return float("nan")
 85         return float(d)
 86 
 87     era = datetime.date(2020, 1, 16)  # NHKまとめのほうの起点日
 88     curdate = era
 89     for line in reader:
 90         date = datetime.date(*map(int, line[:3]))
 91         dadv = (date - curdate).days
 92         data = list(map(_f, line[5:]))
 93         p = _PREFCODE[line[3]]
 94         for i, k in enumerate(_COLSNM):
 95             for _ in range(dadv):
 96                 result[k] = np.vstack((result[k], np.zeros((1, 47))))
 97             result[k][-1, p - 1] = data[i]
 98         curdate = date
 99     # 東洋経済オンラインのデータ更新がおよそ一日遅く、ズレてると扱いにくいので、
100     # 末尾を埋めておく。
101     nhkdat = np.load(io.open("nhk_news_covid19_prefectures_daily_data.npz", "rb"))
102     days = len(nhkdat["newpat"])
103     for d in range(days - len(result[_COLSNM[0]])):
104         for i, k in enumerate(_COLSNM):
105             result[k] = np.vstack((result[k], np.full((1, 47), np.nan)))
106 
107     # 元データについて、どうやら厚労省のオープンデータそのものを転記しただけの
108     # 項目は累積でのみ管理しているらしく、この csv 全体では「累積だったりそうで
109     # なかったり」と大変使いにくい。ゆえ、累積のものは読み替えておく。
110     for k in ("testedPositive", "peopleTested", "discharged", "deaths",):
111         result[k][1:] = (
112             result[k] - np.roll(result[k], 1, axis=0))[1:]
113     #
114     # peopleTested の最新2日分はゴミ:
115     # ---
116     # 2021,1,28,東京都,Tokyo,  97571,  1303823.0,  14957.0,  150,  81767.0,  847,  0.74
117     # 2021,1,29,東京都,Tokyo,  98439,  1309999.0,  13807.0,  147,  83768.0,  864,  0.76
118     # 2021,1,30,東京都,Tokyo,  99208,  1309999.0,  13383.0,  141,  84942.0,  883,  0.77
119     # 2021,1,31,東京都,Tokyo,  99841,  1309999.0,  13257.0,  140,  85698.0,  886,  0.78
120     # ---
121     # 集計が3日遅れ、てこと、だろう。
122     result["peopleTested"][np.where(result["peopleTested"] == 0.0)] = float("nan")
123     # 移動平均
124     for k in ("testedPositive", "peopleTested", "hospitalized", "serious", "deaths"):
125         k1, k2 = k, k + "_movavg"
126         result[k2] = np.zeros(result[k1].shape)
127         for dt in range(len(result[k1])):
128             #result[k2][dt, :] = result[k1][
129             #    max(0, dt - 6):dt + 1, :].mean(axis=0)
130             result[k2][dt, :] = np.nanmean(result[k1][
131                 max(0, dt - 6):dt + 1, :], axis=0)
132     #
133     np.savez(io.open("toyokeizai_online_covid19.npz", "wb"), **result)
134 
135 
136 if __name__ == '__main__':
137     _tonpa("prefectures.csv")

で、人口比の計算に必要なデータ:

pp.txt
 1 1	東京都	13000-1	13,515,271	13,960,236	+3.29	2021年1月1日
 2 2	神奈川県	14000-7	9,126,214	9,216,009	+0.98	2020年9月1日
 3 3	大阪府	27000-8	8,839,469	8,815,082	-0.28	2020年11月1日
 4 4	愛知県	23000-6	7,483,128	7,538,701	+0.74	2020年11月1日
 5 5	埼玉県	11000-1	7,266,534	7,342,682	+1.05	2020年12月1日
 6 6	千葉県	12000-6	6,222,666	6,281,869	+0.95	2021年1月1日
 7 7	兵庫県	28000-3	5,534,800	5,434,645	-1.81	2021年1月1日
 8 8	北海道	01000-6	5,381,733	5,231,685	-2.79	2020年11月30日
 9 9	福岡県	40000-9	5,101,556	5,108,038	+0.13	2020年9月1日
10 10	静岡県	22000-1	3,700,305	3,613,788	-2.34	2021年1月1日
11 11	茨城県	08000-4	2,916,976	2,852,499	-2.21	2021年1月1日
12 12	広島県	34000-6	2,843,990	2,793,470	-1.78	2020年11月1日
13 13	京都府	26000-2	2,610,353	2,566,341	-1.69	2021年1月1日
14 14	宮城県	04000-2	2,333,899	2,290,915	-1.84	2021年1月1日
15 15	新潟県	15000-2	2,304,264	2,195,068	-4.74	2021年1月1日
16 16	長野県	20000-0	2,098,804	2,031,795	-3.19	2021年1月1日
17 17	岐阜県	21000-5	2,031,903	1,975,397	-2.78	2020年9月1日
18 18	栃木県	09000-0	1,974,255	1,930,000	-2.24	2021年1月1日
19 19	群馬県	10000-5	1,973,115	1,924,116	-2.48	2021年1月1日
20 20	岡山県	33000-1	1,921,525	1,880,772	-2.12	2021年1月1日
21 21	福島県	07000-9	1,914,039	1,820,949	-4.86	2021年1月1日
22 22	三重県	24000-1	1,815,865	1,768,632	-2.60	2020年9月1日
23 23	熊本県	43000-5	1,786,170	1,734,231	-2.91	2021年1月1日
24 24	鹿児島県	46000-1	1,648,177	1,587,785	-3.66	2021年1月1日
25 25	沖縄県	47000-7	1,433,566	1,460,427	+1.87	2021年1月1日
26 26	滋賀県	25000-7	1,412,916	1,412,095	-0.06	2021年1月1日
27 27	山口県	35000-1	1,404,729	1,339,003	-4.68	2021年1月1日
28 28	愛媛県	38000-8	1,385,262	1,323,851	-4.43	2021年1月1日
29 29	長崎県	42000-0	1,377,187	1,308,277	-5.00	2021年1月1日
30 30	奈良県	29000-9	1,364,316	1,321,250	-3.16	2021年1月1日
31 31	青森県	02000-1	1,308,265	1,228,730	-6.08	2020年12月1日
32 32	岩手県	03000-7	1,279,594	1,209,457	-5.48	2021年1月1日
33 33	大分県	44000-1	1,166,338	1,124,309	-3.60	2020年11月1日
34 34	石川県	17000-3	1,154,008	1,129,362	-2.14	2020年12月1日
35 35	山形県	06000-3	1,123,891	1,062,239	-5.49	2021年1月1日
36 36	宮崎県	45000-6	1,104,069	1,062,180	-3.79	2021年1月1日
37 37	富山県	16000-8	1,066,328	1,033,981	-3.03	2020年11月1日
38 38	秋田県	05000-8	1,023,119	948,964	-7.25	2021年1月1日
39 39	香川県	37000-2	976,263	949,357	-2.76	2020年9月1日
40 40	和歌山県	30000-4	963,579	912,364	-5.32	2021年1月1日
41 41	山梨県	19000-4	834,930	805,339	-3.54	2021年1月1日
42 42	佐賀県	41000-4	832,832	808,074	-2.97	2021年1月1日
43 43	福井県	18000-9	786,740	762,272	-3.11	2020年11月1日
44 44	徳島県	36000-7	755,733	721,721	-4.50	2020年9月1日
45 45	高知県	39000-3	728,276	688,583	-5.45	2021年1月1日
46 46	島根県	32000-5	694,352	665,702	-4.13	2021年1月1日
47 47	鳥取県	31000-0	573,441	550,651	-3.97	2021年1月1日

の上で、「10万人当たりの新規感染者数に基づくステージング」:

可視化するやーつ
  1 # -*- coding: utf-8-unix -*-
  2 import io
  3 import csv
  4 import datetime
  5 import numpy as np
  6 import matplotlib.pyplot as plt
  7 import matplotlib.font_manager
  8 import matplotlib.ticker as ticker
  9 from matplotlib.colors import BoundaryNorm
 10 
 11 
 12 _PREFCODE = {
 13     '北海道': 1,
 14     '青森県': 2,
 15     '岩手県': 3,
 16     '宮城県': 4,
 17     '秋田県': 5,
 18     '山形県': 6,
 19     '福島県': 7,
 20     '茨城県': 8,
 21     '栃木県': 9,
 22     '群馬県': 10,
 23     '埼玉県': 11,
 24     '千葉県': 12,
 25     '東京都': 13,
 26     '神奈川県': 14,
 27     '新潟県': 15,
 28     '富山県': 16,
 29     '石川県': 17,
 30     '福井県': 18,
 31     '山梨県': 19,
 32     '長野県': 20,
 33     '岐阜県': 21,
 34     '静岡県': 22,
 35     '愛知県': 23,
 36     '三重県': 24,
 37     '滋賀県': 25,
 38     '京都府': 26,
 39     '大阪府': 27,
 40     '兵庫県': 28,
 41     '奈良県': 29,
 42     '和歌山県': 30,
 43     '鳥取県': 31,
 44     '島根県': 32,
 45     '岡山県': 33,
 46     '広島県': 34,
 47     '山口県': 35,
 48     '徳島県': 36,
 49     '香川県': 37,
 50     '愛媛県': 38,
 51     '高知県': 39,
 52     '福岡県': 40,
 53     '佐賀県': 41,
 54     '長崎県': 42,
 55     '熊本県': 43,
 56     '大分県': 44,
 57     '宮崎県': 45,
 58     '鹿児島県': 46,
 59     '沖縄県': 47,
 60 }
 61 def _vis(dattko):
 62     reader = csv.reader(io.open("pp.txt", encoding="cp932"), delimiter="\t")
 63     fac = []
 64     for row in reader:
 65         fac.append((_PREFCODE[row[1]], int(row[4].replace(",", ""))))
 66     fac.sort()
 67     fac = np.array([(100000. / n) for pc, n in fac])  #10万人換算
 68     #fac = np.array([(fac[12][1] / n) for pc, n in fac])  #東京換算
 69     #
 70     era = datetime.date(2020, 1, 16)
 71     #
 72     fontprop = matplotlib.font_manager.FontProperties(
 73         fname="c:/Windows/Fonts/meiryo.ttc")
 74 
 75     era = "2020/1/16"  # 木曜日
 76     era = datetime.date(*map(int, era.split("/")))
 77     #
 78     k, kt = "testedPositive_movavg", "新規陽性者数(10万人あたり)"
 79     #k, kt = "hospitalized", "入院治療等を要する者(10万人あたり)"
 80     #k, kt = "deaths", "死者数(10万人あたり)"
 81     #k, kt = "serious_movavg", "重症者数(10万人あたり)"
 82     #k, kt = "serious", "重症者数(10万人あたり)"
 83     #
 84     skip = 7*4*7
 85     Z = dattko[k] * fac
 86     Z = Z[skip:]
 87     Z[np.where(np.isnan(Z))] = 0
 88     #
 89     T = np.array([era + datetime.timedelta(t + skip) for t in range(len(Z))])
 90     Y = np.arange(48)
 91     X, Y = np.meshgrid(T, Y)
 92     #
 93     fig, ax1 = plt.subplots(tight_layout=True)
 94     fig.set_size_inches(16.53 * 1.2, 11.69)
 95     ax1.invert_yaxis()
 96     ax1.tick_params(labelright=True)
 97     ax1.yaxis.set_major_formatter(
 98         ticker.FuncFormatter(
 99             lambda x, pos=None: {
100                 pc: pn for (pn, pc) in _PREFCODE.items()}.get(x + 1, "")))
101     for lab in ax1.get_yticklabels():
102         lab.set_fontproperties(fontprop)
103     CS = ax1.pcolor(
104         X, Y, Z.T * 7,
105         norm=BoundaryNorm([0, 15, 25, 50], 400), cmap="coolwarm")
106     ax1.set_yticks(range(47))
107 
108     cbar = fig.colorbar(CS)
109     ax1.grid(True)
110     ax1.set_xlabel(
111         "{}\n「東洋経済オンライン「新型コロナウイルス 国内感染の状況」」を加工".format(kt),
112         fontproperties=fontprop)
113 
114     fig.savefig("{}_per100000_{}.png".format(k, str(T[-1])))
115     plt.close(fig)
116 
117 
118 if __name__ == '__main__':
119     _vis(np.load(io.open("toyokeizai_online_covid19.npz", "rb")))

Matplotlib の実例としては「BoundaryNorm例」のつもり。最新データを使うとこんな絵になる:
testedPositive_movavg_per100000_2021-02-19
東京がステージ3基準、大阪がステージ2基準になってるのがわかるだろう。(なお大阪問題に関しては、ワタシはこれまで主張してきた通り、重傷者数と死者数の減りが非常に鈍くいまだにワーストな点から、大阪だけ先走っていいという意見には素直に賛成出来るわけではなく「判断が難しいところ」と思ってる。)

もう一つは、「ステージング」には直結出来ないんだけれど、そのステージングの指標の一つの「前週比」を可視化したらどんななんだろうかと、やってみた:

可視化するやーつ
  1 # -*- coding: utf-8-unix -*-
  2 import io
  3 import datetime
  4 import numpy as np
  5 import matplotlib.pyplot as plt
  6 import matplotlib.font_manager
  7 import matplotlib.ticker as ticker
  8 from matplotlib.colors import BoundaryNorm
  9 
 10 
 11 _PREFCODE = {
 12     '北海道': 1,
 13     '青森県': 2,
 14     '岩手県': 3,
 15     '宮城県': 4,
 16     '秋田県': 5,
 17     '山形県': 6,
 18     '福島県': 7,
 19     '茨城県': 8,
 20     '栃木県': 9,
 21     '群馬県': 10,
 22     '埼玉県': 11,
 23     '千葉県': 12,
 24     '東京都': 13,
 25     '神奈川県': 14,
 26     '新潟県': 15,
 27     '富山県': 16,
 28     '石川県': 17,
 29     '福井県': 18,
 30     '山梨県': 19,
 31     '長野県': 20,
 32     '岐阜県': 21,
 33     '静岡県': 22,
 34     '愛知県': 23,
 35     '三重県': 24,
 36     '滋賀県': 25,
 37     '京都府': 26,
 38     '大阪府': 27,
 39     '兵庫県': 28,
 40     '奈良県': 29,
 41     '和歌山県': 30,
 42     '鳥取県': 31,
 43     '島根県': 32,
 44     '岡山県': 33,
 45     '広島県': 34,
 46     '山口県': 35,
 47     '徳島県': 36,
 48     '香川県': 37,
 49     '愛媛県': 38,
 50     '高知県': 39,
 51     '福岡県': 40,
 52     '佐賀県': 41,
 53     '長崎県': 42,
 54     '熊本県': 43,
 55     '大分県': 44,
 56     '宮崎県': 45,
 57     '鹿児島県': 46,
 58     '沖縄県': 47,
 59 }
 60 def _vis(dattko):
 61     #
 62     era = datetime.date(2020, 1, 16)
 63     #
 64     fontprop = matplotlib.font_manager.FontProperties(
 65         fname="c:/Windows/Fonts/meiryo.ttc")
 66 
 67     era = "2020/1/16"  # 木曜日
 68     era = datetime.date(*map(int, era.split("/")))
 69     #
 70     k, subj = "testedPositive_movavg", "新規陽性判明者数週移動平均"
 71     #
 72     skip = 7*4*7
 73     Z = dattko[k]
 74     # 「前週比」の計算で、分母が1未満だと発散してウザいので、シンプルに1としておく。
 75     Z[np.where(Z < 1)] = 1
 76     Z /= np.roll(Z, 7, axis=0)
 77     Z *= 100
 78     Z = Z[skip:]
 79     Z[np.where(np.isnan(Z))] = 0
 80     #
 81     T = np.array([era + datetime.timedelta(t + skip) for t in range(len(Z))])
 82     Y = np.arange(48)
 83     X, Y = np.meshgrid(T, Y)
 84     #
 85     fig, ax1 = plt.subplots(tight_layout=True)
 86     fig.set_size_inches(16.53 * 1.2, 11.69)
 87 
 88     ax1.invert_yaxis()
 89     ax1.tick_params(labelright=True)
 90     ax1.yaxis.set_major_formatter(
 91         ticker.FuncFormatter(
 92             lambda x, pos=None: {
 93                 pc: pn for (pn, pc) in _PREFCODE.items()}.get(x + 1, "")))
 94     for lab in ax1.get_yticklabels():
 95         lab.set_fontproperties(fontprop)
 96     #
 97     norm = BoundaryNorm([0, 70, 90, 100, 200], 220)
 98     CS1 = ax1.pcolor(X, Y, Z.T, norm=norm, cmap="coolwarm")
 99     ax1.set_yticks(range(47))
100     cbar1 = fig.colorbar(CS1)
101     ax1.set_xlabel(
102         "{}({})\n「東洋経済オンライン「新型コロナウイルス 国内感染の状況」」を加工".format(
103             subj, "前週比"),
104         fontproperties=fontprop)
105     ax1.grid(True)
106     fig.savefig("{}_{}_{}.png".format(k, "cpw", str(T[-1])))
107     #plt.show()
108     plt.close(fig)
109 
110 
111 if __name__ == '__main__':
112     _vis(np.load(io.open("toyokeizai_online_covid19.npz", "rb")))

これも最新データだとこんな絵:
testedPositive_movavg_cpw_2021-02-19
まぁさすがにこれは「こういう質の良い見える化をしたらいいのに」とは言えんね、わかりやすいとは思えん。が、まぁ全体の傾向を一瞥するのには悪くないかもしれんね、今まさに「下がっていってる最中」だということが、最近の「真っ青っぷり」で表現出来てる。

倍加時間」もしくは実効再生産数を計算して推定するのと「前週比」を計算するのは、「過去履歴の傾き」を見ることなので、本質的にはほとんど同じ意味。正直ワタシ的には「やってみた」以上の意味はないのだが、なんというか、この「前週比」を指標にしてるのって、「理解しやすさ」以上の意味があるんだろうか、って思う。しかも「1倍 or not」って、ラフ過ぎるだろ、と。

ほんとはさぁ…、「老婆でも小学生でも」の主張のために、道具を Python + Matplotlib に縛らずに、たとえば Google Drive のワークシートを使うとか、LibreOffice でやるならどんなだとか、あるいは「お仕事としてこういうことをしなければならない人」向けには、WEB で(つまりおそらく javascript で)実現する方法、とかも示せたら良かったのかもしれない、とはずっと思ってるんだよね。まぁ、気が向いたらそっち方面もやってみるかもしれないし、やらないかもしれない。「猿でも出来る」わけではないからね、学習と手間は必要だ、そのための情報を提供できるスキルがワタシにあるのならば、それを伝えるのにはやぶさかではない。


2021-02-21追記:
こういうのはつきものだし、どういう場合であっても「単一の指標だけで判断することは危険」なのだとはいえ、前週比のほう、「一人が一人になりました、前週比1未満を満たしておらず危険ですきゃー」となっちゃってる見せ方はさすがにヒド過ぎるよな、と。「Z[np.where(Z < 1)] = 1」としてただけの部分を以下のように:

変更箇所のみ
1     # 1.「100%」となるケースにおいて、純然たるゴミといえるパターン:
2     #     前週も今週も(限りなく)0
3     #   このケースは「無意味」ということで nan としておく。としないと、まったく問題が
4     #   ない「一人が一人になった!」として「前週比1倍未満ではない」ものとしてピックアップ
5     #   することになり、馬鹿げている。
6     # 2.「前週比」の計算で、分母が1未満だと発散してウザいので、シンプルに1としておく。
7     Z[np.where(np.logical_and(Z < 1, np.roll(Z, 7, axis=0) < 1))] = np.nan
8     Z[np.where(Z < 1)] = 1

とするといいかなと:
testedPositive_movavg_cpw_2021-02-19改
狼少年的になっちゃうのが一番よくないんだよね、こういうのって。


2021-02-22追記:
「人口比 (3?)」への追記というよりは(2)(1)への追記としてのほうが相応しいのかもしれんけれど:

  1 # -*- coding: utf-8-unix -*-
  2 import io
  3 import datetime
  4 import numpy as np
  5 import matplotlib.pyplot as plt
  6 import matplotlib.font_manager
  7 import matplotlib.ticker as ticker
  8 from matplotlib.colors import LogNorm, PowerNorm, LightSource
  9 
 10 
 11 _PREFCODE = {
 12     '北海道': 1,
 13     '青森県': 2,
 14     '岩手県': 3,
 15     '宮城県': 4,
 16     '秋田県': 5,
 17     '山形県': 6,
 18     '福島県': 7,
 19     '茨城県': 8,
 20     '栃木県': 9,
 21     '群馬県': 10,
 22     '埼玉県': 11,
 23     '千葉県': 12,
 24     '東京都': 13,
 25     '神奈川県': 14,
 26     '新潟県': 15,
 27     '富山県': 16,
 28     '石川県': 17,
 29     '福井県': 18,
 30     '山梨県': 19,
 31     '長野県': 20,
 32     '岐阜県': 21,
 33     '静岡県': 22,
 34     '愛知県': 23,
 35     '三重県': 24,
 36     '滋賀県': 25,
 37     '京都府': 26,
 38     '大阪府': 27,
 39     '兵庫県': 28,
 40     '奈良県': 29,
 41     '和歌山県': 30,
 42     '鳥取県': 31,
 43     '島根県': 32,
 44     '岡山県': 33,
 45     '広島県': 34,
 46     '山口県': 35,
 47     '徳島県': 36,
 48     '香川県': 37,
 49     '愛媛県': 38,
 50     '高知県': 39,
 51     '福岡県': 40,
 52     '佐賀県': 41,
 53     '長崎県': 42,
 54     '熊本県': 43,
 55     '大分県': 44,
 56     '宮崎県': 45,
 57     '鹿児島県': 46,
 58     '沖縄県': 47,
 59 }
 60 
 61 
 62 def _get_scales():
 63     import csv
 64     reader = csv.reader(io.open("pp.txt", encoding="cp932"), delimiter="\t")
 65     fac = []
 66     for row in reader:
 67         fac.append((_PREFCODE[row[1]], int(row[4].replace(",", ""))))
 68     fac.sort()
 69     return (
 70         (lambda z: z, "", ""),
 71         (lambda z: z * np.array([(100000. / n) for pc, n in fac]), "(10万人換算)", "_per100000"),
 72         (lambda z: z * np.array([(fac[12][1] / n) for pc, n in fac]), "(東京換算)", "_asTokyo"),
 73         (lambda z: z / np.nanmax(z, axis=0), "(都道府県ごと最大値との比)", "_emphPeaks"),
 74         )
 75 
 76 
 77 def _vis(dattko):
 78     scales = _get_scales()
 79     #
 80     era = datetime.date(2020, 1, 16)
 81     #
 82     fontprop = matplotlib.font_manager.FontProperties(
 83         fname="c:/Windows/Fonts/meiryo.ttc")
 84 
 85     era = "2020/1/16"  # 木曜日
 86     era = datetime.date(*map(int, era.split("/")))
 87     #
 88     skip = 7*4*8 + 7*2
 89     jmeanings = {
 90         "testedPositive": "新規陽性者数",
 91         "testedPositive_movavg": "新規陽性者数7日移動平均",
 92         "hospitalized": "入院治療等を要する者",
 93         "hospitalized_movavg": "入院治療等を要する者7日移動平均",
 94         "serious": "重症者数",
 95         "serious_movavg": "重症者数7日移動平均",
 96         "deaths": "死者数",
 97         "deaths_movavg": "死者数7日移動平均",
 98         }
 99     #
100     for i, (k0, fac, norm) in enumerate((
101             ("serious", scales[0], PowerNorm(0.4)),
102             ("serious", scales[1], PowerNorm(0.4)),
103             ("serious", scales[2], PowerNorm(0.4)),
104             ("serious", scales[3], PowerNorm(3)),
105             ("deaths", scales[0], PowerNorm(0.4)),
106             ("deaths", scales[1], PowerNorm(0.4)),
107             ("deaths", scales[2], PowerNorm(0.4)),
108             ("deaths", scales[3], PowerNorm(3)),
109             ("hospitalized", scales[0], PowerNorm(0.4)),
110             ("hospitalized", scales[1], PowerNorm(0.4)),
111             ("hospitalized", scales[2], PowerNorm(0.4)),
112             ("hospitalized", scales[3], PowerNorm(3)),
113     )):
114         for k in (k0, k0 + "_movavg"):
115             Z = fac[0](dattko[k])
116             Z = Z[skip:]
117             Z[np.where(np.isnan(Z))] = 0
118             Z[np.where(Z < 0)] = 0
119             T = np.array([era + datetime.timedelta(t + skip) for t in range(len(Z))])
120             Y = np.arange(48)
121             #
122             X, Y = np.meshgrid(T, Y)
123             #
124             fig, ax1 = plt.subplots(tight_layout=True)
125             fig.set_size_inches(16.53 * 1.2, 11.69)
126     
127             ax1.invert_yaxis()
128             ax1.tick_params(labelright=True)
129             ax1.yaxis.set_major_formatter(
130                 ticker.FuncFormatter(
131                     lambda x, pos=None: {pc: pn for (pn, pc) in _PREFCODE.items()}.get(x + 1, "")))
132             for lab in ax1.get_yticklabels():
133                 lab.set_fontproperties(fontprop)
134             #
135             CS1 = ax1.pcolor(X, Y, Z.T, norm=norm, cmap="coolwarm")
136             ax1.set_yticks(range(47))
137             cbar1 = fig.colorbar(CS1)
138             ax1.set_xlabel(
139                 "{}{}\n「東洋経済オンライン「新型コロナウイルス 国内感染の状況」」を加工".format(
140                     jmeanings[k], fac[1]),
141                 fontproperties=fontprop)
142             ax1.grid(True)
143             fig.savefig("{}{}_{}.png".format(k, fac[2], str(T[-1])))
144             #plt.show()
145             plt.close(fig)
146 
147 
148 if __name__ == '__main__':
149     _vis(np.load(io.open("toyokeizai_online_covid19.npz", "rb")))

色んな指標ごとに都度都度やってたらダルくなってきて、の、少しだけ整理したヤツね。最新データを使えばこんな絵になる:
hospitalized_2021-02-21hospitalized_asTokyo_2021-02-21hospitalized_per100000_2021-02-21hospitalized_emphPeaks_2021-02-21
hospitalized_movavg_2021-02-21hospitalized_movavg_asTokyo_2021-02-21hospitalized_movavg_per100000_2021-02-21hospitalized_movavg_emphPeaks_2021-02-21

serious_2021-02-21serious_asTokyo_2021-02-21serious_per100000_2021-02-21serious_emphPeaks_2021-02-21
serious_movavg_2021-02-21serious_movavg_asTokyo_2021-02-21serious_movavg_per100000_2021-02-21serious_movavg_emphPeaks_2021-02-21

deaths_2021-02-21deaths_asTokyo_2021-02-21deaths_per100000_2021-02-21deaths_emphPeaks_2021-02-21
deaths_movavg_2021-02-21deaths_movavg_asTokyo_2021-02-21deaths_movavg_per100000_2021-02-21deaths_movavg_emphPeaks_2021-02-21
ここ何日か大阪問題について言ってきたが、改めて。注目して欲しいのは「重症者数(東京換算)」と「死亡者数(東京換算)」。後者にあまり差がないことにより、「大阪の重傷者数が多いのは東京と基準が異なるからである」というところの説得力が失われ、すなわち「本当に重症者が多い」ということがわかる、という点。ワタシがずっと主張してるのが、つまりそういうこと。「皆で同じ可視化を共有し、同じ問題意識を持つべきである」てこと。こういうちゃんとした分析的視点をこそ、メディアが率先して伝えるべきなのではないのか、てことね。

大阪問題、といってるけれど、正直2月5日頃に吉村知事が言い出した頃と今とではだいぶ状況が変わってる。そこで「300を切ったら」と言い出した時にはさすがに絶句した。というか、300という数字そのものより、2月7日の週にでも解除できる、という姿勢が見受けられたことに絶句したんだけどね。けど、今は一応重症者数も緩やかながら減少に転じ始めてるし、言ってることも「2月末」と随分慎重になった。「2月末」の妥当性についても、正直重症者数の数からいえばもろ手で賛成とは言わないけど、絶対論外だというほどでもない。

本日時点だと、大阪より東京の方が気になる傾向にあるやもしれん。死者数の推移がね、今ちょっと…:
nhk_news_covid19_prefectures_daily_data5_2_9_東京都_0221
報道でも伝えられる通り、新規感染者数の下げ止まりも気にはなるけれど、死者数が若干上昇傾向なのは、注視しておく必要がありそうだと思う。大阪はこんなね:
nhk_news_covid19_prefectures_daily_data5_2_9_大阪府_0221


2021-03-05追記:
今回の追記も、「良い可視化を共有せよ」という主張を本題とするなら蛇足、その手法の方についての追記。スクリプトの改善は「前週比」の方だけで示す:

  1 # -*- coding: utf-8-unix -*-
  2 import io
  3 import datetime
  4 import numpy as np
  5 import matplotlib.pyplot as plt
  6 import matplotlib.font_manager
  7 import matplotlib.ticker as ticker
  8 from matplotlib.colors import BoundaryNorm, Normalize
  9 
 10 
 11 _PREFCODE = {
 12     '北海道': 1,
 13     '青森県': 2,
 14     '岩手県': 3,
 15     '宮城県': 4,
 16     '秋田県': 5,
 17     '山形県': 6,
 18     '福島県': 7,
 19     '茨城県': 8,
 20     '栃木県': 9,
 21     '群馬県': 10,
 22     '埼玉県': 11,
 23     '千葉県': 12,
 24     '東京都': 13,
 25     '神奈川県': 14,
 26     '新潟県': 15,
 27     '富山県': 16,
 28     '石川県': 17,
 29     '福井県': 18,
 30     '山梨県': 19,
 31     '長野県': 20,
 32     '岐阜県': 21,
 33     '静岡県': 22,
 34     '愛知県': 23,
 35     '三重県': 24,
 36     '滋賀県': 25,
 37     '京都府': 26,
 38     '大阪府': 27,
 39     '兵庫県': 28,
 40     '奈良県': 29,
 41     '和歌山県': 30,
 42     '鳥取県': 31,
 43     '島根県': 32,
 44     '岡山県': 33,
 45     '広島県': 34,
 46     '山口県': 35,
 47     '徳島県': 36,
 48     '香川県': 37,
 49     '愛媛県': 38,
 50     '高知県': 39,
 51     '福岡県': 40,
 52     '佐賀県': 41,
 53     '長崎県': 42,
 54     '熊本県': 43,
 55     '大分県': 44,
 56     '宮崎県': 45,
 57     '鹿児島県': 46,
 58     '沖縄県': 47,
 59 }
 60 
 61 
 62 class MidpointNormalize(Normalize):
 63     def __init__(self, vmin=None, vmax=None, vcenter=None, clip=False):
 64         self.vcenter = vcenter
 65         Normalize.__init__(self, vmin, vmax, clip)
 66 
 67     def __call__(self, value, clip=None):
 68         # I'm ignoring masked values and all kinds of edge cases to make a
 69         # simple example...
 70         x, y = [self.vmin, self.vcenter, self.vmax], [0, 0.5, 1]
 71         return np.ma.masked_array(np.interp(value, x, y))
 72 
 73 
 74 def _vis(dattko):
 75     #
 76     era = datetime.date(2020, 1, 16)
 77     #
 78     fontprop = matplotlib.font_manager.FontProperties(
 79         fname="c:/Windows/Fonts/meiryo.ttc")
 80 
 81     era = "2020/1/16"  # 木曜日
 82     era = datetime.date(*map(int, era.split("/")))
 83     #
 84     k, subj = "testedPositive_movavg", "新規陽性判明者数週移動平均"
 85     #
 86     skip = 7*4*7
 87     Z = dattko[k]
 88     # 1.「100%」となるケースにおいて、純然たるゴミといえるパターン:
 89     #     前週も今週も(限りなく)0
 90     #   このケースは「無意味」ということで nan としておく。としないと、まったく問題が
 91     #   ない「一人が一人になった!」として「前週比1倍未満ではない」ものとしてピックアップ
 92     #   することになり、馬鹿げている。
 93     # 2.「前週比」の計算で、分母が1未満だと発散してウザいので、シンプルに1としておく。
 94     Z[np.where(np.logical_and(Z < 1, np.roll(Z, 7, axis=0) < 1))] = np.nan
 95     Z[np.where(Z < 1)] = 1
 96     Z /= np.roll(Z, 7, axis=0)
 97     Z *= 100
 98     Z = Z[skip:]
 99     Z[np.where(np.isnan(Z))] = 0
100     #
101     T = np.array([era + datetime.timedelta(t + skip) for t in range(len(Z))])
102     Y = np.arange(48)
103     X, Y = np.meshgrid(T, Y)
104     #
105     for i, norm in enumerate((
106             BoundaryNorm([0, 70, 90, 100, 200], 220),
107             MidpointNormalize(vcenter=90, vmax=200),
108     )):
109         fig, ax1 = plt.subplots(tight_layout=True)
110         fig.set_size_inches(16.53 * 1.2, 11.69)
111     
112         ax1.invert_yaxis()
113         ax1.tick_params(labelright=True)
114         ax1.yaxis.set_major_formatter(
115             ticker.FuncFormatter(
116                 lambda x, pos=None: {
117                     pc: pn for (pn, pc) in _PREFCODE.items()}.get(x + 1, "")))
118         for lab in ax1.get_yticklabels():
119             lab.set_fontproperties(fontprop)
120         #
121         CS1 = ax1.pcolor(X, Y, Z.T, norm=norm, cmap="coolwarm")
122         ax1.set_yticks(range(47))
123         cbar1 = fig.colorbar(CS1)
124         # 最大値詐称していることをカラーバーで主張(「> 200」という具合に)しておく。としないと誤解を招く。
125         cbticklbls = [tlobj.get_text() for tlobj in cbar1.ax.get_yticklabels()]
126         cbticklbls[-1] = "> " + cbticklbls[-1]
127         cbar1.ax.set_yticklabels(cbticklbls)
128         #
129         ax1.set_xlabel(
130             "{}({})\n「東洋経済オンライン「新型コロナウイルス 国内感染の状況」」を加工".format(
131                 subj, "前週比"),
132             fontproperties=fontprop)
133         ax1.grid(True)
134         fig.savefig("{}_{}_{}_{}.png".format(k, "cpw", str(T[-1]), i + 1))
135         #plt.show()
136         plt.close(fig)
137 
138 
139 if __name__ == '__main__':
140     _vis(np.load(io.open("toyokeizai_online_covid19.npz", "rb")))

ポイントは二つで、一つはカラーバーのティックラベルの変更、もう一つは「Custom normalization」の実例。前者は些細だけどプレゼンテーションとしては結構重要で、後者は実は前からやりたくてやり方を知らなくて、初めてやってみたやつ。絵はこうなる:
testedPositive_movavg_cpw_2021-03-04_1
testedPositive_movavg_cpw_2021-03-04_2
「後者」の目的には、最新の matplotlib.colors には TwoSlopeNorm があるのでそれを使えばいいのだが、ワタシがインストールしてるバージョンにはまだ含まれてないので、まさに Custom normalization に書かれてるそのままを使った。何をしたかったのかはたぶん絵を見ればわかってもらえると思う…んだけど日本語で説明しづらいんだよな…。滑らかな階調で見たくて、なおかつ、その色の傾きが可変、てことなのよ。例えば標高の可視化を例にすると、「500m以上はほとんど茶色でいいが、500m以下の傾きが詳細にわかればいいなぁ」みたいなときに使うやつね。今ワタシがやってみせたヤツだと「前週比が90%以下まではあんまし重要でないのでほとんど真っ青でよろし」という制御をしてる、てこと。


2021-03-27追記:
「良い見える化を共有すべし、いまや小学生でも実現できるほどに容易なのだからなおさら」てのが、一連の話の主題の半分だった、てことはいい?

そうなんだけれど、実際問題 matplotlib はさすがに多少多めの鍛錬が必要で、さすがにそこいらの一般人が初見で怯んでしまう程度には難解なのは問題で、そういう意味で、ワタシの主張を弱めてしまう理由にもなっている。

「matplotlib だから」である、これは。matplotlib の名前の由来はおそらく「matrix plotting library」であろう。これからわかるように、最初に念頭に置かれたのは「行列の可視化」であり、すなわち「NumPy/SciPyとセットで」使われることが前提になっている。となると、これはもうユーザの層が自ずと偏ってきて、つまり「科学技術系」(特に自然科学系)に属する学者や技術者が好むものとなった。「難しいのが当たり前」の世界にいるもんだから、matplotlib に対しても「もっと簡単にすべきだ」なんて流れにはならん気がする。実際ワタシ自身も結構そっち寄り。

で、「ワタシには matplotlib があるのでおk (あるいは R とか gnuplot でもいい)」で思考停止状態だったとも言えるけれど、今回このようなシリーズを書いた手前、ちったぁ「もっと誰にでも」ネタがあればいいなぁ、とは思っていたの。そんなタイミングで、たまたま glances の(optional な)依存パッケージであるところの「pygal」に気付いた。おそらく WEB ベース前提のプロジェクトと思われる。

一つには結構古参でなおかつ、「WEBで」の方向性が SVG である点は、「いまどきの HTML5 時代」にはそぐわないのかもしれんけれど、いやまぁ、なにせ簡単だ、これは:

 1 >>> import pygal
 2 >>> from math import cos
 3 >>> xy_chart = pygal.XY()
 4 >>> xy_chart.title = 'XY Cosinus'
 5 >>> xy_chart.add('x = cos(y)', [(cos(x / 10.), x / 10.) for x in range(-50, 50, 5)])
 6 <pygal.graph.xy.XY object at 0x000001F221B730F0>
 7 >>> xy_chart.add('y = cos(x)', [(x / 10., cos(x / 10.)) for x in range(-50, 50, 5)])
 8 <pygal.graph.xy.XY object at 0x000001F221B730F0>
 9 >>> xy_chart.add('x = 1',  [(1, -5), (1, 5)])
10 <pygal.graph.xy.XY object at 0x000001F221B730F0>
11 >>> xy_chart.add('x = -1', [(-1, -5), (-1, 5)])
12 <pygal.graph.xy.XY object at 0x000001F221B730F0>
13 >>> xy_chart.add('y = 1',  [(-5, 1), (5, 1)])
14 <pygal.graph.xy.XY object at 0x000001F221B730F0>
15 >>> xy_chart.add('y = -1', [(-5, -1), (5, -1)])
16 <pygal.graph.xy.XY object at 0x000001F221B730F0>
17 >>> xy_chart.render_to_png("pygal_example_xychart.png")
18 Traceback (most recent call last):
19   File "<stdin>", line 1, in <module>
20   File "C:\Program Files\Python35\lib\site-packages\pygal\graph\public.py", line 118, in render_to_png
21     import cairosvg
22 ImportError: No module named 'cairosvg'
23 >>> #xy_chart.render_to_file("pygal_example_xychart.svg")
24 >>> xy_chart.render_data_uri()
25    ...

cairo って Windows 版あったかなぁ? あれば png 化も出来ることになる。まぁ WEB 相手だけなら SVG だけでもありがたかろうね。結果はこんなね: