実測値は Gompertz 曲線に当てはまるか？　その3

2020年1月に「ゴンペルツ曲線とは何か？」という連載を始めて1年たち、10章目となりました。
2021年もよろしくお願い申し上げます。

「2枚の三角定規の法則」の一般化

ゴンペルツ曲線とは何か？（8）では COVID-19 新規感染者数がゴンペルツ曲線に従う時、
流行開始からピークまでの発生数：ピークから終息までの発生数 ≒ 1 : √3
となり、2枚の三角定規に当てはまることを示しました。
45°、45°、90° の三角定規と 30°、60°、90° の三角定規を用いました。
2つの直角三角形は底辺の比率が 1 : √3 であれば以下のように高さを適当に変えることが可能です。

上り坂の角度を α 、下り坂の角度を β 、高さを h とすると以下が成り立ちます。

\begin{align}
\tan\alpha&=h\\
\tan\beta&=\frac{h}{\sqrt{3}}\\
\ &=\frac{\tan\alpha}{\sqrt{3}}\\
\ &=\frac{\sqrt{3}}{3}\tan\alpha\\
\therefore\beta&=\arctan\biggl(\frac{\sqrt{3}}{3}\tan\alpha\biggr)\\
\end{align}\\

よって下り坂の角度 β は上り坂の角度 α に規定されることがわかります。
※arctan: tan の逆関数。tan θ = y となる角度 θ （ラジアン）は、θ = arctan y である。

α：上り坂	β：下り坂
10°	5.8°
20°	11.9°
30°	18.4°
40°	25.8°
45°	30°
50°	34.5°
60°	45°
70°	57.8°
80°	73.0°
85°	81.4°

上り坂が急峻な時は下り坂も急峻になります。
以下に日本各地の第2波を示します。

\beta=\arctan\biggl(\frac{\sqrt{3}}{3}\tan\alpha\biggr)

が成り立つ（=ゴンペルツ曲線に従う）ことがわかります。

日本の第3波に『2枚の三角定規の法則』は当てはまるか

ゴンペルツ曲線とは何か？（8）で日本の第3波に2枚の三角定規を当てはめ、2021.1.20頃に終息すると予想しました。以下が2020.12.3時点での予測です。

結果は予測を大きく外れてしまいました。申し訳ございません。
2021.1.23の時点で以下のようになりました。

もう1つの新たな波が重なっているように見えます。
この傾向は北海道と大阪でも見られます。（参考のために同時期の札幌駅前と梅田駅前の人出¹も示します。）

いずれも

\beta=\arctan\biggl(\frac{\sqrt{3}}{3}\tan\alpha\biggr)

が成り立ち、ゴンペルツ曲線に従うことがわかります。
流行の拡大期と収束期で駅前の人出に変化がないことから、接触抑制の強弱はリアルタイムの実効再生産数よりむしろ、ゴンペルツ効果の強弱に関連しているのではないかと考えました。つまり、「マスク」「手洗い」「三密回避」などの行動変容をしているコホートで第○波の流行が始まった場合、行動変容の強度で規定されるゴンペルツ効果に基づいたゴンペルツ曲線を描く、ということです。新たな波の出現には、変異株などウイルス側の要因があるかも知れません。
上り坂が急峻な時は下り坂も急峻になるので終息は2021.2月中旬になる見込みです。（新たな波が来なければですが・・・。また外れたらすみません。）

2枚の三角定規からゴンペルツ曲線 G(t) を求める。

ゴンペルツ曲線（対数化していないもの）は

\ G(t)=ab^{e^{-kt}}

と二重指数関数になります。a が最大値、k が減衰係数となります。G(t) は累計感染者数を表し、新規感染者数は G(t) の1日毎の差分になります。ゴンペルツ曲線は最大値と減衰係数が決まれば b が 0 < b < 1 の範囲でどのような値になっても曲線の形状は同一です。（ただし x 軸方向の平行移動は生じます。）a を G_max とし、b を e に置き換えることも可能です。（この場合、e > 1 なので指数の係数にマイナス "–" をつけます。G(0) = G_max/e を通る位置に平行移動しています。）

\ G(t)=G_{max}e^{-e^{-kt}}

上式を時間 t で微分してみましょう。

\begin{align}
\ G(t)&=G_{max}e^{-e^{-kt}}\\
\ U=-e^{-kt}&,　V=-kt　とすると\\
\ G(t)&=G_{max}e^U\\
\ U&=-e^V\\
\frac{dG}{dt}&=\frac{dG}{dU}\frac{dU}{dV}\frac{dV}{dt}\\
\frac{dG}{dU}&=G_{max}e^U\\
\frac{dU}{dV}&=-e^V\\
\frac{dV}{dt}&=-k\\
\frac{dG}{dt}&=G_{max}e^U(-e^V)(-k)\\
\ &=kG_{max}e^Ue^V\\
\ &=kG_{max}e^{-e^{-kt}}e^{-kt}\\
\ G'(0)&=\frac{kG_{max}}{e}\\
\therefore k&=\frac{eG'(0)}{G_{max}}
\end{align}\\

G(0) = G_max/e を通る位置に平行移動しており、ゴンペルツ曲線は最大値の 1/e の時にピークとなるので、G'(0) はピーク時の1日あたり新規発生数と等しくなります。G_max は流行開始からピークまでの累計発生数を e 倍したものに等しいので、下図のように最初の三角形の面積の e 倍となります。

上記を踏まえて減衰係数 k を求めると、

\begin{align}
\ k&=\frac{eG'(0)}{G_{max}}\\
\ &=\frac{e\timesピーク時発生数}{開始〜ピークまでの日数\timesピーク時発生数\times\frac{1}{2}\times e}\\
\ &=\frac{2}{開始〜ピークまでの日数}
\end{align}\\

このように k と G_max は「開始〜ピークまでの日数」と「ピーク時発生数」から算出できます。求めた値を

\ G(t)=G_{max}e^{-e^{-kt}}

に代入すれば良いです。
日本の第2波で開始〜ピークまでの日数（33日）、ピーク時発生数（1600）より
k = 0.0606、G_max = 71763 となり、

\ G(t)=71763e^{-e^{-0.0606t}}

が得られました。新規感染者数は G(t) の1日毎の差分であり、 G(t) の1日毎の差分は dG/dt に等しくなります。

\begin{align}
\frac{dG}{dt}&=kG_{max}e^{-e^{-kt}}e^{-kt}\\
\ &=kG_{max}e^{-kt-e^{-kt}}\\
\ &=eG'(0)e^{-kt-e^{-kt}}\\
\ &=e\times ピーク時発生数 \times e^{-kt-e^{-kt}}\\
\ &=4349e^{-0.0606t-e^{-0.0606t}}\\
\end{align}\\

と表せます。

・・・良く fit しています。
e^-0.0606 = 0.9412 なので日本の第2波では、累計感染者数を対数表示した時の傾きが、2020.6.30〜2020.9.16 までの期間は毎日 0.9412 倍に減衰していたことになります。（一定減衰仮定）

上り坂の途中でピークを予測できるか？

上記にて「開始〜ピークまでの日数」と「ピーク時発生数」という2つのパラメータから流行の全体像を描けることを示しました。ではピークアウトする前の上り坂の途中でピークを予測することは可能でしょうか。
以下は米国の第2波です。

上り坂の途中（黄色の矢印）の時点で小さな赤い三角定規のコースになるか大きな青い三角定規のコースになるかは予測できません。
以下はどちらもゴンペルツ曲線ですが

流行初期の経過は重なっていますがピーク値に大きな違いがあります。
東京海洋大学准教授　勝川俊雄先生は流行初期のデータにゴンペルツ曲線を当てはめて未来予測をすることの問題点を以下のように述べています。

もし仮に、感染者数がゴンペルツ曲線に従うとしても、ノイズが多い感染初期のデータだけを用いて、曲線の全体図を正確に把握するのは、ほぼ不可能です。
『K値による予測を使うべきでない理由』²より

したがってピークアウトする前の上り坂の途中でピークを予測することは出来ないと思われます。
しかし、

\beta=\arctan\biggl(\frac{\sqrt{3}}{3}\tan\alpha\biggr)

から、上り坂の途中でも下り坂の角度を予測することが可能です。

大阪の第2波を検証する

比較的きれいにゴンペルツ曲線が当てはまる大阪の第2波で検証してみました。第2波の始まりからピークまでの日数（30日）とピーク時発生数（200人）を当てはめました。

\begin{align}
\ k&=\frac{2}{開始〜ピークまでの日数}\\
\ &=\frac{2}{30}\\
\ &=0.0667\\
\ G_{max}&=開始〜ピークまでの日数\times ピーク時発生数\times \frac{1}{2}\times e\\
&=30\times 200\times \frac{1}{2}\times e\\
&=8155\\
\ G(t)&=G_{max}e^{-e^{-kt}}\\
\ &=8155e^{-e^{-0.0667t}}\\
\ G'(t)&=e\times ピーク時発生数 \times e^{-kt-e^{-kt}}\\
\ &=544e^{-0.0667t-e^{-0.0667t}}\\
\end{align}\\

わずか2つのパラメータで流行の始まりから終息まできれいに fit しています。二重指数関数の減衰係数は k = 0.0667 と1つの値に決定することができます。梅田駅前の人出に大きな変化がなかったことから、この期間の行動変容の強度は定常であったと考えられます。ゴンペルツ曲線では定常の行動変容下における拡大から収束までのトレンド変化を以下のように解釈します。

「マスク」「手洗い」「三密回避」等の行動様式を持つ集団の中で流行がスタートする。
ウイルス側の要因と人間側の要因で規定されるゴンペルツ効果に基づき、累計感染者数 G(t) を対数表示した時の傾き d ln G(t)/dt が日々 e^-kt 倍に減衰する。（一定減衰仮定）

指数関数の当てはめ

大阪の第2波を二重指数関数（ゴンペルツ曲線）ではなく指数関数（SIRモデル）で当てはめると以下のようになります。

全経過を6つのタイムゾーンに区切って6つの指数係数を使います。時間の経過と共に指数係数が漸減していることがわかります。実効再生産数は指数係数の一次関数であり、平均罹病期間を14日とすると 実効再生産数 = 14 × 指数係数 + 1 で表せます。

タイムゾーン	指数係数	実効再生産数
① 7/8~7/17	0.293	5.102
② 7/18~7/27	0.118	2.652
③ 7/28~8/6	0.028	1.392
④ 8/7~8/16	-0.018	0.748
⑤ 8/17~8/26	-0.042	0.412
⑥ 8/27~	-0.054	0.244

大阪の第2波では、山の始まりの時点から実効再生産数は減り続けています。

指数関数（SIRモデル）で当てはめることの利点は、

感染者数の推移がゴンペルツ曲線に当てはまらない時でも使える。
タイムゾーン毎の実効再生産数の変化がわかる。

ということです。

ゴンペルツ曲線とは何か？（10）

実測値は Gompertz 曲線に当てはまるか？ その3