丸一年思い続けてたことである。
新型コロナウィルスが問題になり始めた昨年一月には問題意識を持っていたオレ、かっけー。
なわけはなく、一回目の緊急事態宣言の半月ほど前からようやく危機意識を持ち始めた、要するにごく普通の凡人である。先見の明なんかありやしない。3月中旬だったと思うなぁ、ちゃんと「これはやべー」と思ったのって。
ただ、その頃からずっと感じ、今でもそうだと思っていることがある。というかいくつもあるんだけどね、そのうちの一つについて。それは、ある違和感についてである。
最初にでっかいクエスチョンマークが頭に飛び交ったのは、一回目の緊急事態宣言を出す出さないでモタモタしてた際に、「今は緊急事態宣言を出す状態にない」であるとか「オーバーシュートには至っていない」といった説明が飛び交った際である。なんで「んんん????」と思ったか。これは、ワタシは日々の新規陽性者報告の推移を聞き、「あぁ今日も順調に予定通り増えとるわ、やべーー」と毎日思っていたことと真っ向から対立する主張だったからだ。そしてこの感覚は今でも絶賛継続中である。
確か3月中旬にはもうメディアは「二週間前の我々の行動の反映が今日の陽性判明者数である」と散々説明していたと記憶している。であるから、当然「科学」は新規陽性判明者について「最低でも2週間先の予測」を目指す必要があるということになるし、少なくとも報道を真面目にウォッチしている一般人も、そのことを理解しているわけだ。
が、どうにも現在の分科会、かつての専門家会議の説明、それに政府も、「最低でも2週間先の予測」に基づいた発言に全く思えないと感じ続けた一年間だったのだ、てことをね、ちゃんと言いたくなったわけさ。「遅い、後手後手だ」という日本語に至る前の根拠にあたる部分としてね、「先を見越した対策を打とうとしてない」と見え続けた、てこと。もっとヒドい言い方がいい? 「二週間前の我々の行動の反映」にのみ基づいてるように見え続けたんだよ、ぶっちゃけ。平たくいえば「予測してない」。
じゃぁこの「予測」って、難しいの? 簡単なの?
当たり前だが、まず「人間は予知能力なんて持ってない」ので、「絶対に外れない予測」というのはほとんど存在しない。科学・数学を駆使すれば、100%未満の範囲内で「良い予測」と「良くない予測」があり、一定以上になると、精度の良しあしによらず「そこそこは難解」でもあり、なかなかに大学レベルになりがちであり、そしてそして、一般人が行う素人予測のほとんどはさほどあてにならない、ということが「一般的には」言える。
話はここから。「一般的には」といった。「一般人が行う素人予測」って、実際なにを指すかというとこれは要するに「線形の外挿」であり、非常に感覚的に理解しやすく、また誰もが日常的に利用していて、「当てるための予測」に使うには難があったとしても、「数日後にびっくりしないための覚悟」目的には実は結構適ってるんじゃないのか、と、ワタシは思ったのだ。
そりゃね、教科書的には最低でも最小二乗法だのほか関数近似だのの「大学生レベル」の数学を試みる方が普通は良い予測になるし、感染症においてはモデルに基づくシュミレーションの手法が比較的確立しているようであるから、まかせられるのならば専門家の予測を待つのが最も正確な理解が出来る。けれども「専門家の予測がないと一般人はまったく予測が不可能なので、昨日感染者が1000人だったので来週は1570人になるかもしれない、などと想像出来るはずはないし、しちゃいけない」となるかといえば…、これはならない。うん、ならなかった。
やってみたわけよ。毎日感覚的に「リニアに外挿してなんとなくの予測」してたのを、ちゃんと機械的に。個人的にはもっと高等な数学を知ってるけどそうじゃなくて「史上最もバカに属する予測」の実力ってどんなよ、てのが知りたくてね。
元にするデータは、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 "newdeath": np.zeros((1, 47)),
15 "newdeathmovavg": np.zeros((1, 47)),
16 }
17 era = "2020/1/16"
18 era = datetime.datetime(*map(int, era.split("/")))
19 reader = csv.reader(io.open(fn, encoding="utf-8-sig"))
20 next(reader)
21 for row in reader:
22 dt, (p, c, d) = row[0], map(int, row[1::2])
23 dt = (datetime.datetime(*map(int, dt.split("/"))) - era).days
24 if dt > len(result["newpat"]) - 1:
25 for k in result.keys():
26 result[k] = np.vstack((result[k], np.zeros((1, 47))))
27 result["newpat"][dt, p - 1] = c
28 result["newpatmovavg"][dt, p - 1] =
29 result["newpat"][max(0, dt - 7):dt + 1, p - 1].mean()
30 result["newdeath"][dt, p - 1] = d
31 result["newdeathmovavg"][dt, p - 1] =
32 result["newdeath"][max(0, dt - 7):dt + 1, p - 1].mean()
33 np.savez(io.open("nhk_news_covid19_prefectures_daily_data.npz", "wb"), **result)
34
35
36 if __name__ == '__main__':
37 fn = "nhk_news_covid19_prefectures_daily_data.csv"
38 _tonpa(fn)
スクリプトをマジメに読んだ人は気付いたかもしれない。「8日平均」を計算しちゃってる。一連の「主張」の本題とはあまり関係ないとはいえ、さすがにもうちょっと慎重にやるべきでした、すまぬ。正しくは「dt – 6」「dt – 5」…「dt – 1」「dt」の7エレメントでの平均を取らなければならないので、「slice(dt – 6, dt + 1)」。8日平均が役に立たないわけではないけれど、やりたいこととやってることが違うので、間違いは間違い。グラフは滑らかさが若干減って、もうほんの少しだけガタガタになる。
1 # -*- coding: utf-8 -*-
2 # 13:00記: すまん、npatmadif について、二重に間違えた結果偶然正しい、ってダメなやつだった。これ、
3 # 「傾き」のつもりなので「7日前との差」を取ったなら7で割って使う必要がある「つもりで作った」はず
4 # なのだが、これを使う部分で「7日後」を計算するのに7をかけていないミスをしてる。つまり「結果的に
5 # 正しい」。スクリプトは修正しないが、なんらか同じようなことをしようと思う場合はコードの見かけに
6 # 騙されないように注意。
7 import io
8 import datetime
9 import numpy as np
10 import matplotlib as mpl
11 import matplotlib.pyplot as plt
12 import matplotlib.font_manager
13
14
15 def _vis(dat):
16 era = "2020/1/16"
17 era = datetime.datetime(*map(int, era.split("/")))
18 npat = dat["newpat"][:,12] # 12: 東京 (都道府県コード13)
19 npatma = dat["newpatmovavg"][:,12]
20 #npat = dat["newpat"].sum(axis=1)
21 #npatma = dat["newpatmovavg"].sum(axis=1)
22 npatmadif = npatma - np.roll(npatma, 7)
23
24 fontprop = matplotlib.font_manager.FontProperties(fname="c:/Windows/Fonts/meiryo.ttc")
25
26 T = [era + datetime.timedelta(t) for t in range(len(npat))]
27 #for t0 in range(0, len(T), 28):
28 for t0 in range(28, len(T)):
29 tlast = min(len(T) - 1, t0 + 28)
30 cspan = slice(t0, tlast + 1)
31 pspan = slice(tlast, tlast + 7)
32 T2 = [era + datetime.timedelta(t) for t in range(pspan.start, pspan.stop)]
33 dtbeg, dtend = era + datetime.timedelta(t0), era + datetime.timedelta(tlast)
34 fig, ax1 = plt.subplots(tight_layout=True)
35 fig.set_size_inches(16.53, 11.69)
36 l1, = ax1.plot(T[cspan], npat[cspan], color="blue")
37 l2, = ax1.plot(T[cspan], npatma[cspan], color="orange")
38 # 数ある「予測」のなかでも最も単純(バカ)な部類に属する予測をあえて。
39 # たぶん中学生程度の学力で十分に理解でき、かつ普通の素人がかなり真っ先に思いつくやつ。
40 # こむずかしい用語なら「線形外挿」ということにはなるが、よーは「びろーんと延ばす」。
41 # ここでは週平均の傾きを延長している。(傾きは、週平均とさらに一週前週平均の差を採用。)
42 last = npat[cspan][-1]
43 cdif = npatmadif[cspan.stop - 7:cspan.stop]
44 for pr in (cdif.min(), cdif.mean(), cdif.max()):
45 prd0, prd1 = last, last + pr
46 ax1.plot(T2, np.linspace(prd0, prd1 + 1, len(T2)), linestyle="dotted", color="black")
47 tit = "{}の新規陽性判明数(東京)の素人予測: {}人n(予測実施日: {})".format(
48 str(dtend + datetime.timedelta(7)).partition(" ")[0],
49 int(last + cdif.mean()), str(dtend).partition(" ")[0])
50 ax1.set_title(tit, fontproperties=fontprop)
51 l3, = ax1.plot(T2[:len(npat[pspan])], npat[pspan], color="blue", alpha=0.5)
52 ax1.legend(
53 (l1, l2), ("新規陽性判明数", "新規陽性判明数(週平均)"),
54 shadow=True, prop=fontprop)
55 ax1.set_xlabel(
56 "「NHK新型コロナウイルス関連データ」を加工", fontproperties=fontprop)
57 plt.axvspan(T2[0], T2[-1], facecolor="lightgray", alpha=0.5)
58 fig.autofmt_xdate()
59 fig.savefig("result/{:04d}.png".format(t0))
60 #plt.show()
61 plt.close(fig)
62
63
64 if __name__ == '__main__':
65 fn = "nhk_news_covid19_prefectures_daily_data.npz"
66 _vis(np.load(io.open(fn, "rb")))
作った絵の一枚:
結合してビデオにしたんだけど、なんかアップロードが調子悪い。うまくいったら追記で貼り付ける。
「リニアに外挿」といっている実際の処理は、「週平均の傾きを素直に延長」としてる。点線がそれ。リニアに、なので、「指数関数的に」となればこの予測から上に振れるので「指数関数的に増えたーっ」とビックリすることが出来るし、そうでなくてこのアホ予測の範囲内に収まれば、驚かないことが出来る。さてさて。「分科会とやら」は、この「アホ予測」が予測できないような「予測困難性」に基づいてびっくりしてたようにみえるかい? いや、ワタシにはそうな全然みえない。「わかってたじゃん」と。最低でも一週間前にあらかじめビックリしておくことは、このアホ予測ですら結構可能であり、少なくとも「ギリギリ持ちこたえている」と頭の良いお方々がたがたが言い続けていた時分は「いやぁ、もうダメじゃね?」と素人が理解することが出来たはずだ、てのがワタシの偽らざる感覚ね。だってさぁ、「一週間後予測」でさえ「一週間前の我々の行動の反映」なんだぞ? 「実況値だけみて議論してねー?」と思い続けたのは、そういう意味。
なお、こうした機械予測は「自分でやらんといかん」てことはないしそれを推奨しているのでもなく、実際 google の予測があるし、感染症学の専門家が出す予測も比較的簡単に手に入るので、普通はそうしたものを能動的に見に行くのがベストである。けれどもワタシが言いたいのは、今の政治家どもが「今すぐ緊急事態宣言を出す状況にない」だの「ギリギリ持ちこたえている」だのが「予測に基づいているようには到底思えない」と思わされてる以上は、願いたいのは「専門家による正しい予測がなくても感覚的な予測は割と簡単に出来るのだから、毎日テレビが報じたらいいんじゃね?」てことなのよ。現況だけでなくシンプルでも一週間後の予測が周知されるだけで、世論の誘導が今よりずっとやりやすくなると思うし、頭のオカシイ政治家の発言の良しあしも「中学生が指摘できる」ようになると思うわけよね。「西浦教授がこんなどどーん予測を出しましたー、びっくりびっくり」と、3か月に一回だけびっくりしてたってアカンと思うんだわ。こういうのさ、テレビ会社には「科学班」があるはずよね、そうしたところが自前でやる手もあるんじゃないかな。高等なやつは確かに結構知識やスキルが必要だし、「あまり科学的でない」という但し書きが常に必要なことは注意は必要だけれどもね、専門家任せにしない以上は。
感染症学の専門家の予測がなかなか出ないという場合に、たぶんだけれど、経済学の専門家に頼るといいんじゃないのかな、とちょっと思ってる。株価予想につかう数学の応用も、それほど大外ししないんではないのかしら、って気がしてる。
14:20追記:
ビデオのアップロードが出来た。
あと二つ目のスクリプトに訂正コメントを追加済だが、そこに追加したのとは別の2つの問題がある。一つは確信犯なのに説明しなかったことで、「負の予測値」に対する措置をしてないこと。説明すりゃいっしょ、とわかってたのに書き忘れたの。もう一つは終了処理で、意図としてはダウンロードしたデータの日付の7日後までの処理のつもりだったのに、28日後までの処理になっちゃってる。ビデオはこの余分な末尾の絵を取り除いてあるが、スクリプトの問題はそのままなので注意。
あと念のためもう一度注意しておくけれど、「自力でやるべきである」てことを言いたいのではないよ。マスメディア、特にテレビがなんらか毎日予測を市民に提供して欲しい、てのが言いたいこと。
