マスコミを巧妙に活用した情宣作戦というのか，"やってる感パフォーマンス"の演出が大好きなのは維新と安倍だが，土曜日の夕方にまたまた"記者会見もどき"の宣伝工作を行っていた。内容はないようだったが，いつものように滑舌が悪くて"指示"が"Siri"になってしまうため，どこでもDIGA視聴中のiPadがよけいな応答をして画面が止まってしまうのだった...orz
さて，あいかわらず醜悪なポエムに包まれた左右首振りプロンプター官僚作文の中で，一点だけ定量的な発言として興味を引いたのが，感染者数が2週間で30倍以上になるというくだりだった。さっそく前回のメモ基本再生産数と集団免疫率の式を使って考えてみよう。
とりあえず，簡単のためここで登場した感染者数はSIRモデルの$I(t)$のことだとしよう。初期条件が，$S(0) = N, I(0) \ll 1, R(0)=0$の場合，$t=0$近傍の方程式は，次のように近似される。
\begin{equation}
\begin{aligned}
\dfrac{dI}{dt} &=\beta \dfrac{S}{N} I -\dfrac{I}{\gamma} \\
& \approx \bigl( \beta -\dfrac{1}{\gamma} \bigr) I \\
&= \dfrac{R_0 - 1}{\gamma} I
\end{aligned}
\end{equation}
ただし，$R_0 = \beta \gamma $を用いた。これから次の解が得られる。
\begin{equation}
I(t) = I(0) \ e^{\frac{R_0 - 1}{\gamma} t}
\end{equation}
$k$倍加時間を$\tau_k$として，
\begin{equation}
\begin{aligned}
k &= \dfrac{I(\tau_k)}{I(0)} = e^{\frac{R_0 - 1}{\gamma} \tau_k} \\
\log k &= \frac{R_0 - 1}{\gamma} \tau_k \\
\therefore R_0 &= 1 + \dfrac{\gamma \log k}{\tau_k}
\end{aligned}
\end{equation}
そこで，2週間で30倍≒2倍加時間が3日であることから，$k=2, \gamma=5, \tau=3$をあてはめると，この場合の$R_0=1+ \dfrac{5*0.6931}{3} = 2.15$となる。
付録A：
感染者数というのが$I(t)$でなく，通常示されている新規感染数累計（Confirmed）（→$J(t)$と定義）のことだと解釈するとどうなるだろうか。
\begin{equation}
J(t) = \int_0^t I(t') dt' \approx \int_0^t I(0) e^{\frac{R_0 - 1}{\gamma} t'} dt' = \dfrac{I(0) \gamma}{R_0-1} \Bigl\{ e^{\frac{R_0-1}{\gamma} t } - 1 \Bigr\}
\end{equation}
２週間で30倍ということは，$30 = \frac{J(t_0+14)}{J(t_0)}$であるが，これを数値計算すると$t_0$=3.8日となる。まあ，そういう解がありうるという以上の情報はないので，これは$R_0$に対する制約条件にはならないということか。
追伸（2020.4.1）
付録B：
新規感染数累計の増倍率では，感染数の指数関数的な増大$f(t)= a e^{t/\lambda}$が観測される場合の一定期間における累計の増倍率$r(t)$を次のように定義した。
\begin{equation}
r(t) = \dfrac{e^{(t+\tau)/ \lambda}-1}{e^{ t/\lambda}-1} = \dfrac{e^{\tau/ \lambda}-e^{-t/\lambda}}{1-e^{ -t/\lambda}} \approx e^{\tau/ \lambda}
\end{equation}
これを上記の$J(t)$と組み合わせると，$\lambda=\dfrac{\gamma}{R_0-1}$とすればよい。
したがって，
\begin{equation}
\lim_{t \to \infty}r(t) = e^{\frac{\tau (R_0-1)}{ \gamma}}
\end{equation}
2020年3月29日日曜日
2020年3月28日土曜日
湖北省の集団免疫率
基本再生産数と集団免疫率からの続き
昨日の基本再生産数と集団免疫率のところで，湖北省（武漢以外）における集団免疫率をSIIDR2モデルで推定したところ，0.2%という結果が得られた。ところで湖北省の感染中心である武漢はどうだろうか。そこで，武漢を含む湖北省について同様の計算を行った。
パラメタを決定するための新規感染数累計と死亡数累計のデータとしては，WHOのSituation Reportを使う。ただし，データが存在しているのは，2月1日のNo.12から3月15日のNo.55までであり，死亡数データはNo.25からNo.55までに限定されている。また2月13日のNo.24から2月14日No.25以降の新規感染数累計のデータには，確定感染者（Confirmed）の定義変更に伴うギャップが存在する。そこで，次のようにデータを加工した。
①定義変更ギャップについては変更前後のデータ倍率を古いデータにも当てはめて新しい定義にそろえる。
②欠落している初期部分については指数関数的な振る舞いを仮定して，未欠落部分につががるように補完する。こちらの方は，データ数もそれほど多くないので，後のモデルパラメタ推定にはそれほど大きな影響を与えない。
2月1日を基準日0として50日分の以下のjuliaコード用のデータが得られた。xhは時間（日）
yhは新規感染数累計，zhは死亡数累計である。なお湖北省の人口を5900万人として，1万人当りの数に換算している。
xh=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,
22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,
43,44,45,46,47,48,49]
yh=[0.370,0.478,0.612,0.777,0.979,1.223,1.517,1.944,2.420,2.957,
3.684,4.387,4.982,5.677,6.225,6.872,7.428,7.885,8.318,8.811,
9.221,9.534,9.861,10.168,10.455,10.514,10.621,10.755,10.862,
10.896,10.981,11.049,11.118,11.172,11.244,11.340,11.373,11.393,
11.412,11.435,11.456,11.469,11.476,11.482,11.485,11.487,11.488,
11.489,11.490,11.491]
zh=[0.028,0.031,0.035,0.039,0.043,0.049,0.054,0.060,0.067,0.075,
0.084,0.094,0.105,0.117,0.130,0.145,0.162,0.181,0.202,0.223,
0.247,0.271,0.287,0.303,0.326,0.344,0.363,0.381,0.398,0.423,
0.434,0.443,0.448,0.455,0.462,0.468,0.475,0.480,0.487,0.492,
0.497,0.502,0.506,0.510,0.513,0.516,0.518,0.519,0.521,0.523]
SIIDR2モデルのパラメタを次のようにとると，上記データがほぼ再現できる。
$\alpha_1=5/0.80, \alpha_2=5/0.20, \gamma_1=15/0.96, \gamma_2=15/0.04, \beta=0.88, \nu=0.30, \lambda=7, \tau=5$。もちろんパラメタには不定性があるが，データの再現性を制約条件とした場合に，感染数や回復数（免疫獲得数）などの結果が大きく変わることはない。以下の図では，u2:軽症感染数（無症状含む），u3:重症感染数（隔離済），u4:死亡数累計=Deaths，u5:新規感染数累計=Confirmed，u6:回復数（免疫獲得数）である。
昨日の基本再生産数と集団免疫率のところで，湖北省（武漢以外）における集団免疫率をSIIDR2モデルで推定したところ，0.2%という結果が得られた。ところで湖北省の感染中心である武漢はどうだろうか。そこで，武漢を含む湖北省について同様の計算を行った。
パラメタを決定するための新規感染数累計と死亡数累計のデータとしては，WHOのSituation Reportを使う。ただし，データが存在しているのは，2月1日のNo.12から3月15日のNo.55までであり，死亡数データはNo.25からNo.55までに限定されている。また2月13日のNo.24から2月14日No.25以降の新規感染数累計のデータには，確定感染者（Confirmed）の定義変更に伴うギャップが存在する。そこで，次のようにデータを加工した。
①定義変更ギャップについては変更前後のデータ倍率を古いデータにも当てはめて新しい定義にそろえる。
②欠落している初期部分については指数関数的な振る舞いを仮定して，未欠落部分につががるように補完する。こちらの方は，データ数もそれほど多くないので，後のモデルパラメタ推定にはそれほど大きな影響を与えない。
2月1日を基準日0として50日分の以下のjuliaコード用のデータが得られた。xhは時間（日）
yhは新規感染数累計，zhは死亡数累計である。なお湖北省の人口を5900万人として，1万人当りの数に換算している。
xh=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,
22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,
43,44,45,46,47,48,49]
yh=[0.370,0.478,0.612,0.777,0.979,1.223,1.517,1.944,2.420,2.957,
3.684,4.387,4.982,5.677,6.225,6.872,7.428,7.885,8.318,8.811,
9.221,9.534,9.861,10.168,10.455,10.514,10.621,10.755,10.862,
10.896,10.981,11.049,11.118,11.172,11.244,11.340,11.373,11.393,
11.412,11.435,11.456,11.469,11.476,11.482,11.485,11.487,11.488,
11.489,11.490,11.491]
zh=[0.028,0.031,0.035,0.039,0.043,0.049,0.054,0.060,0.067,0.075,
0.084,0.094,0.105,0.117,0.130,0.145,0.162,0.181,0.202,0.223,
0.247,0.271,0.287,0.303,0.326,0.344,0.363,0.381,0.398,0.423,
0.434,0.443,0.448,0.455,0.462,0.468,0.475,0.480,0.487,0.492,
0.497,0.502,0.506,0.510,0.513,0.516,0.518,0.519,0.521,0.523]
SIIDR2モデルのパラメタを次のようにとると，上記データがほぼ再現できる。
$\alpha_1=5/0.80, \alpha_2=5/0.20, \gamma_1=15/0.96, \gamma_2=15/0.04, \beta=0.88, \nu=0.30, \lambda=7, \tau=5$。もちろんパラメタには不定性があるが，データの再現性を制約条件とした場合に，感染数や回復数（免疫獲得数）などの結果が大きく変わることはない。以下の図では，u2:軽症感染数（無症状含む），u3:重症感染数（隔離済），u4:死亡数累計=Deaths，u5:新規感染数累計=Confirmed，u6:回復数（免疫獲得数）である。
図１ 湖北省（2/1〜3/15）の１万人当り感染状況
つまり，人口5900万人の湖北省では1万人当り約60人の水準（0.6%）で免疫獲得数が飽和している。これを前回の結果（武漢以外の湖北省で1万人当り約20人）と組み合わせると，武漢市の免疫獲得者は5900*60-4800*20=258,000人となり，武漢市の集団免疫率は2.6%になる。これでもまだ少ないのでいつでも再燃する可能性がある。
このパラメタにおける実効再生産数 $R_{\rm eff}= \alpha \beta(t) = \alpha \dfrac{\beta}{15}\{8+7 \tanh(-\frac{t-\tau}{\lambda}) \} $は次のようになる。$t$=20日までの平均で$R_{\rm eff}=1.415$，$t$=30日までの平均は$R_{\rm eff}=1.044$となる。
図２ 湖北省データを再現する実効再生産数 R_eff（t=0〜49）
2020年3月27日金曜日
基本再生産数と集団免疫率
感染症数理の基本的なモデルであるSIRモデルについて自分は十分理解できていなかった。奈良女子大の高須さん[1]や神戸大の牧野さん[2]や京大の門さん[3]の資料が大変わかりやすくて勉強になった。
きっかけはツイッター@nagayaさんの次のコメント。
SIRモデルは，感受性保持者（Susceptible，$S(t)$）・感染者（Infected，$I(t)$）・免疫保持者（Recovered，$R(t)$ ）からなる力学系モデルであり次の微分方程式を満たす（$\beta, \gamma$の定義に注意 ）。ここで$t$は時間であり，最も単純なバージョンでは集団の人数を$N$として$S(t)+I(t)+R(t)=N$は保存量である。なお，$\beta$は感染率（１人単位時間当りの感染数），$\gamma$は感染期間である。
\begin{equation}
\begin{aligned}
\dfrac{dS}{dt} &= -\beta \dfrac{S}{N} I\\
\dfrac{dI}{dt} &= \beta \dfrac{S}{N} I - \dfrac{I}{\gamma} \\
\dfrac{dR}{dt} &= - \dfrac{I}{\gamma}\\
\end{aligned}
\end{equation}
これからわかること。集団中の感受性保持者（免疫非保持者）の割合を$p(t)=\dfrac{S}{N} $とすると，$p(t) > \dfrac{1}{\beta \gamma}$で感染拡大，$p(t) < \dfrac{1}{\beta \gamma}$で感染縮小となる。そこで，基本再生産数を$R_0 = \beta \gamma$と定義すると，$p(0)=1$の初期状態において，$R_0>1$で感染拡大，$R_0<1$で感染縮小となる。なお，実効再生産数が$R_{\rm eff}=R_0 p(t)$ と定義される。
また集団免疫率を$\pi(t) = 1-p(t)$と定義すると，感染拡大と縮小の分岐点$\tau$では，$1-\pi(\tau) = \dfrac{1}{R_0}$である。これが@nagayaさんの示した関係式だった。
次に十分時間がたって$I(\infty)=0$となるときの集団免疫率$\pi(\infty)$と$R_0$の関係を求めてみよう。もとの微分方程式系の第１式と第３式から次の微分方程式が得られる。
\begin{equation}
\begin{aligned}
\dfrac{dS}{dR} &= -R_0 \dfrac{S}{N}\\
\dfrac{dS}{S} &= -R_0 \dfrac{dR}{N}\\
\log S + C &= -R_0 \dfrac{R}{N} = -R_0 (1-\dfrac{I}{N} - p(t) ) = -R_0 (\pi(t) - i(t) )\\
\end{aligned}
\end{equation}
ただし，$i(t)=\dfrac{I}{N}$であり，$i(\infty)=0$。そこで，$S(0)=N$として$\pi(\infty)=\pi$とすると，
\begin{equation}
\begin{aligned}
\dfrac{S(t)}{S(0)} &= e^{-R_0\{ \pi(t) - i (t) \}} \\
p(t) = 1-\pi(t) &= e^{-R_0\{ \pi(t) - i (t) \}} \\
1- \pi &= e^{ -R_0 \pi } \\
R_0 &= - \dfrac{1}{\pi}\log (1- \pi )
\end{aligned}
\end{equation}
分岐点における$\pi(t)$と$R_0$の関係式とは異なった関係を与えていることに注意しよう。
付録：juliaにおける上記プロットのコード
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
using Plots
gr()
x = 0.05:0.05:0.95
f(x) = .- 1 ./ x .* log.(1 .- x)
y = f(x)
plot(x, y, label="")
savefig("sir.png")
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
追伸：
以上は簡単なSIRモデルの場合だが，湖北省（武漢以外）をおおむね説明したSIIDR2モデルにたいして，十分な時間がたったときの湖北省（武漢以外）の集団免疫率はどのように推定できるだろうか。パラメタセット，$\alpha_1=5/0.80, \alpha_2=5/0.20, \gamma_1=15/0.96, \gamma_2=15/0.04, \beta=0.915, \nu=0.025, \lambda=7, \tau=7$について計算すると，Recovered $u_5$は1万人当り19程度に収束し，集団免疫率はたかだか0.2%にどどまる。
［１］伝染病のモデル−大域情報学（高須夫悟）
［２］コロナウィルス(SARS-CoV-2)の(単純な) 数学モデル（牧野淳一郎）
［３］この感染は拡大か収束か：再生産数 R の物理的意味と決定（門信一郎）
きっかけはツイッター@nagayaさんの次のコメント。
R0＝2なら国民の50%、3なら66%、逆に1.5なら33%の免疫獲得で理論上は収束に向かいます（この理論値を超えて免疫獲得した分が「オーバーシュート」）。日本は自粛とクラスター追跡、検査を絞って1.2程度にして医療崩壊を避けようとしたんでしょう。無理です。長期化しますし、指数関数的に潜在します。基本再生産数$R_0$と集団免疫率に関係あったっけ？ということで高須さんの資料で学び直し。
SIRモデルは，感受性保持者（Susceptible，$S(t)$）・感染者（Infected，$I(t)$）・免疫保持者（Recovered，$R(t)$ ）からなる力学系モデルであり次の微分方程式を満たす（$\beta, \gamma$の定義に注意 ）。ここで$t$は時間であり，最も単純なバージョンでは集団の人数を$N$として$S(t)+I(t)+R(t)=N$は保存量である。なお，$\beta$は感染率（１人単位時間当りの感染数），$\gamma$は感染期間である。
\begin{equation}
\begin{aligned}
\dfrac{dS}{dt} &= -\beta \dfrac{S}{N} I\\
\dfrac{dI}{dt} &= \beta \dfrac{S}{N} I - \dfrac{I}{\gamma} \\
\dfrac{dR}{dt} &= - \dfrac{I}{\gamma}\\
\end{aligned}
\end{equation}
これからわかること。集団中の感受性保持者（免疫非保持者）の割合を$p(t)=\dfrac{S}{N} $とすると，$p(t) > \dfrac{1}{\beta \gamma}$で感染拡大，$p(t) < \dfrac{1}{\beta \gamma}$で感染縮小となる。そこで，基本再生産数を$R_0 = \beta \gamma$と定義すると，$p(0)=1$の初期状態において，$R_0>1$で感染拡大，$R_0<1$で感染縮小となる。なお，実効再生産数が$R_{\rm eff}=R_0 p(t)$ と定義される。
また集団免疫率を$\pi(t) = 1-p(t)$と定義すると，感染拡大と縮小の分岐点$\tau$では，$1-\pi(\tau) = \dfrac{1}{R_0}$である。これが@nagayaさんの示した関係式だった。
次に十分時間がたって$I(\infty)=0$となるときの集団免疫率$\pi(\infty)$と$R_0$の関係を求めてみよう。もとの微分方程式系の第１式と第３式から次の微分方程式が得られる。
\begin{equation}
\begin{aligned}
\dfrac{dS}{dR} &= -R_0 \dfrac{S}{N}\\
\dfrac{dS}{S} &= -R_0 \dfrac{dR}{N}\\
\log S + C &= -R_0 \dfrac{R}{N} = -R_0 (1-\dfrac{I}{N} - p(t) ) = -R_0 (\pi(t) - i(t) )\\
\end{aligned}
\end{equation}
ただし，$i(t)=\dfrac{I}{N}$であり，$i(\infty)=0$。そこで，$S(0)=N$として$\pi(\infty)=\pi$とすると，
\begin{equation}
\begin{aligned}
\dfrac{S(t)}{S(0)} &= e^{-R_0\{ \pi(t) - i (t) \}} \\
p(t) = 1-\pi(t) &= e^{-R_0\{ \pi(t) - i (t) \}} \\
1- \pi &= e^{ -R_0 \pi } \\
R_0 &= - \dfrac{1}{\pi}\log (1- \pi )
\end{aligned}
\end{equation}
図１ t=∞における集団免疫率π（横軸）と基本再生産数R0（縦軸）
分岐点における$\pi(t)$と$R_0$の関係式とは異なった関係を与えていることに注意しよう。
付録：juliaにおける上記プロットのコード
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
using Plots
gr()
x = 0.05:0.05:0.95
f(x) = .- 1 ./ x .* log.(1 .- x)
y = f(x)
plot(x, y, label="")
savefig("sir.png")
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
追伸：
以上は簡単なSIRモデルの場合だが，湖北省（武漢以外）をおおむね説明したSIIDR2モデルにたいして，十分な時間がたったときの湖北省（武漢以外）の集団免疫率はどのように推定できるだろうか。パラメタセット，$\alpha_1=5/0.80, \alpha_2=5/0.20, \gamma_1=15/0.96, \gamma_2=15/0.04, \beta=0.915, \nu=0.025, \lambda=7, \tau=7$について計算すると，Recovered $u_5$は1万人当り19程度に収束し，集団免疫率はたかだか0.2%にどどまる。
図２ 湖北省（武漢以外）のSIIDR2モデルよる集団免疫率の評価
［２］コロナウィルス(SARS-CoV-2)の(単純な) 数学モデル（牧野淳一郎）
［３］この感染は拡大か収束か：再生産数 R の物理的意味と決定（門信一郎）
2020年3月26日木曜日
パーセク
天文学における距離の単位に，光年とならんでパーセクがあるが実感としてわかっていなかった。もちろん定義も。
どこかで，１天文単位（太陽地球間の距離）を見込む角が１秒になる距離がパーセクであるという定義をみかけた。散歩中に暗算してみたが，その過程で太陽半径（地球半径の100倍，64万km）は地球−月軌道半径（38万km）より大きいことに驚いてしまった。まあ，いわれてみればそうである。で，太陽を見込む視角と比較しようとしたが，暗算力が落ちていてよくわからなかった。
まず，１秒をラジアンに直す。１秒 = 1/3600度 = 1/3600 * π/180 ラジアン = 1/(6^3*10^3) ラジアン ≒ 1/(2*10^5) ラジアンかな。１パーセク× tan 1秒 = １天文単位なので，１パーセク =
１天文単位 ÷ {1/(2*10^5)}，１天文単位は1.5×10^8 kmなので，１パーセクは3×10^13 kmだ。１光年は 3×10^5 km /s * π×10^7 s = 10^13 km なので，１パーセクは３光年くらいかな。しらんけど。
どこかで，１天文単位（太陽地球間の距離）を見込む角が１秒になる距離がパーセクであるという定義をみかけた。散歩中に暗算してみたが，その過程で太陽半径（地球半径の100倍，64万km）は地球−月軌道半径（38万km）より大きいことに驚いてしまった。まあ，いわれてみればそうである。で，太陽を見込む視角と比較しようとしたが，暗算力が落ちていてよくわからなかった。
まず，１秒をラジアンに直す。１秒 = 1/3600度 = 1/3600 * π/180 ラジアン = 1/(6^3*10^3) ラジアン ≒ 1/(2*10^5) ラジアンかな。１パーセク× tan 1秒 = １天文単位なので，１パーセク =
１天文単位 ÷ {1/(2*10^5)}，１天文単位は1.5×10^8 kmなので，１パーセクは3×10^13 kmだ。１光年は 3×10^5 km /s * π×10^7 s = 10^13 km なので，１パーセクは３光年くらいかな。しらんけど。
2020年3月25日水曜日
DeepL翻訳
googleからスピンオフしたチームによる「DeepL翻訳」のクオリティが高いという評判だったので試してみた。
One hundred years have passed since Ryunosuke Akutagawa published "The Spider's Thread". The theme of the class costume procession in the autumn school festival in the second year of high school was the spider thread. He cut the bamboo of the Buddha's paradise tower from the neighborhood and drew a picture of the flames and smoke of hell. When he was a hell-raiser flocking to the spider thread (the mountaineer's red foil) after Gubata.
カンダタがGubataになっているものの，なんとなくうまくできているようだ。タイトルの On a thread of the web は，ウエブのスレッドということでした。以前のみらい翻訳と比べてどうだろうか。一長一短かもしれない。
One hundred years have passed since Ryunosuke Akutagawa published "The Spider's Thread". The theme of the class costume procession in the autumn school festival in the second year of high school was the spider thread. He cut the bamboo of the Buddha's paradise tower from the neighborhood and drew a picture of the flames and smoke of hell. When he was a hell-raiser flocking to the spider thread (the mountaineer's red foil) after Gubata.
カンダタがGubataになっているものの，なんとなくうまくできているようだ。タイトルの On a thread of the web は，ウエブのスレッドということでした。以前のみらい翻訳と比べてどうだろうか。一長一短かもしれない。
2020年3月24日火曜日
翔んで兵庫
先日，大阪府の吉村知事は，春分の日を含む三連休に大阪府と兵庫県の間での不要不急の往来自粛を要請した。法的根拠を指摘されると，厚生労働省から非公開の文書が出ていたのでそれを重視して府知事の判断で要請を発したと弁明した。
ところがこの「大阪府・兵庫県における緊急対策の提案(案)（3月16日，北海道大学西浦博教授）」が公開されるとSNS上ではさっそく突っ込みが入り，危険区域は大阪府・兵庫県であり，「大阪府・兵庫県内外の不要不急な往来の自粛を呼びかける」とは大阪府と兵庫県の間だけの遮断を意味するのではないだろう，読解力はないのかという声があがった。
吉村知事の判断は，半分はリーズナブルで，半分はいつものような維新の粗っぽい決断力顕示欲（パフォーマンス）が見られるような気がする（本気ならば関西広域連合の定例会にも出席し，兵庫県と打ち合わせした上で決断するべきだろう）。それはさておき，ここではシミュレーションによってその判断の合理性を検討する（ただし，吉村も大阪府もまともに分析・吟味はしていないだろう，していれば説明できるはず）。
問題を，「感染者が急増している地域とそれほどでもない地域の間の交流の遮断が，感染者の増加率やピーク日にどのような影響があるのか」というふうに設定しよう。簡単のためにSIRモデルを用いる。大阪府の人口は880万人，兵庫県への移動人口は11万人（1.25%→pとおく），兵庫県の人口は550万人，大阪府への移動人口は33万人（6.00%→qとおく）であり，感染確率は大阪府的A地域で$\beta_a$=0.20，兵庫県的B地域で$\beta_b$=0.25とし，B地域の方がより感染が拡大していると仮定する（本当のところは，大阪府のPCR検査拒否率がトップであることを考えるとよくわからない）。なお，3/22現在の感染者数は，大阪府125人（14ppm），兵庫県107人（19ppm）であるが，初期条件はともに10ppmから計算をはじめる。（注：この移動人口は平日の場合だろうと思われるが，これを適用する。）
２つの地域A,Bに対して，SIRモデルをあてはめ，感染者の移動人口分も感染に寄与すると仮定すると次式が成り立つ。計算はExcelで1000人を1単位，時間は1日を1単位として差分（以下では微分方程式で表現）により実行する。したがって，人口の初期値はA地域で$n_a$=8800，B地域で$n_b$=5500とする。また平均回復日数を$\gamma$=10とする。$p,q$は各地域の移動人口の割合である。
\begin{equation}
\begin{aligned}
\dfrac{dS_a}{dt} &= - \beta_a \dfrac{S_a}{n_a} \{ I_a - (1-p) q I_b \} \\
\dfrac{dI_a}{dt} &= \beta_a \dfrac{S_a}{n_a} \{ I_a - (1-p) q I_b \} - \dfrac{I_a}{\gamma} \\
\dfrac{dR_a}{dt} &= \dfrac{I_a}{\gamma} \\
\dfrac{dS_b}{dt} &= - \beta_b \dfrac{S_b}{n_b} \{ I_b - (1-q) p I_a \} \\
\dfrac{dI_b}{dt} &= \beta_b \dfrac{S_b}{n_b} \{ I_b - (1-q) p I_a \} - \dfrac{I_b}{\gamma} \\
\dfrac{dR_b}{dt} &= \dfrac{I_b}{\gamma}
\end{aligned}
\end{equation}
ところがこの「大阪府・兵庫県における緊急対策の提案(案)（3月16日，北海道大学西浦博教授）」が公開されるとSNS上ではさっそく突っ込みが入り，危険区域は大阪府・兵庫県であり，「大阪府・兵庫県内外の不要不急な往来の自粛を呼びかける」とは大阪府と兵庫県の間だけの遮断を意味するのではないだろう，読解力はないのかという声があがった。
吉村知事の判断は，半分はリーズナブルで，半分はいつものような維新の粗っぽい決断力顕示欲（パフォーマンス）が見られるような気がする（本気ならば関西広域連合の定例会にも出席し，兵庫県と打ち合わせした上で決断するべきだろう）。それはさておき，ここではシミュレーションによってその判断の合理性を検討する（ただし，吉村も大阪府もまともに分析・吟味はしていないだろう，していれば説明できるはず）。
問題を，「感染者が急増している地域とそれほどでもない地域の間の交流の遮断が，感染者の増加率やピーク日にどのような影響があるのか」というふうに設定しよう。簡単のためにSIRモデルを用いる。大阪府の人口は880万人，兵庫県への移動人口は11万人（1.25%→pとおく），兵庫県の人口は550万人，大阪府への移動人口は33万人（6.00%→qとおく）であり，感染確率は大阪府的A地域で$\beta_a$=0.20，兵庫県的B地域で$\beta_b$=0.25とし，B地域の方がより感染が拡大していると仮定する（本当のところは，大阪府のPCR検査拒否率がトップであることを考えるとよくわからない）。なお，3/22現在の感染者数は，大阪府125人（14ppm），兵庫県107人（19ppm）であるが，初期条件はともに10ppmから計算をはじめる。（注：この移動人口は平日の場合だろうと思われるが，これを適用する。）
２つの地域A,Bに対して，SIRモデルをあてはめ，感染者の移動人口分も感染に寄与すると仮定すると次式が成り立つ。計算はExcelで1000人を1単位，時間は1日を1単位として差分（以下では微分方程式で表現）により実行する。したがって，人口の初期値はA地域で$n_a$=8800，B地域で$n_b$=5500とする。また平均回復日数を$\gamma$=10とする。$p,q$は各地域の移動人口の割合である。
\begin{equation}
\begin{aligned}
\dfrac{dS_a}{dt} &= - \beta_a \dfrac{S_a}{n_a} \{ I_a - (1-p) q I_b \} \\
\dfrac{dI_a}{dt} &= \beta_a \dfrac{S_a}{n_a} \{ I_a - (1-p) q I_b \} - \dfrac{I_a}{\gamma} \\
\dfrac{dR_a}{dt} &= \dfrac{I_a}{\gamma} \\
\dfrac{dS_b}{dt} &= - \beta_b \dfrac{S_b}{n_b} \{ I_b - (1-q) p I_a \} \\
\dfrac{dI_b}{dt} &= \beta_b \dfrac{S_b}{n_b} \{ I_b - (1-q) p I_a \} - \dfrac{I_b}{\gamma} \\
\dfrac{dR_b}{dt} &= \dfrac{I_b}{\gamma}
\end{aligned}
\end{equation}
図 SIRモデルシミュレーションの対数プロット
その結果，次のことがわかった。なお，感染ピークは感染者数 $I(t)$ のピーク日，感染倍率は $I(t)$ を等比級数に近似した場合の公比/日である。ピーク日が前倒しされ，感染倍率が増えるほうが危険度が増すということになる。
①両地域（$\beta_a$=0.20, $\beta_b$=0.25）を分離した場合はそれぞれ次のようになる。
（感染ピークＡ 94日，感染倍率A 1.099，感染ピークB 68日，感染倍率B 1.148）
②①の両地域がp=1.25% q=6.00%で結合した場合は次のようになる。
（感染ピークＡ 81日，感染倍率A 1.116，感染ピークB 67日，感染倍率B 1.151）
③②の条件で地域Bの感染率が低い場合（$\beta_b=0.15$）は次のようになる。
（感染ピークＡ 93日，感染倍率A 1.102，感染ピークB 130日，感染倍率B 1.066）
したがって，自分の地域よりも感染率の高い地域は遮断したほうがよい（感染者数が倍増する日数が7.5日から6.5日になる程度）。感染率の低い地域を遮断してもあまり影響はない（それでも若干危険度は増大する）。ただし現実の大阪・兵庫がこのモデルのケースに当てはまっているかどうかは断定できない。モデルが単純すぎるかもしれないし，基礎データの信頼度が低いと思われることによる。
2020年3月18日水曜日
久々のPerl
久しぶりにPerlのプログラムを書いた。20年ぶりとか13年ぶりとか8年ぶりとかの感じ。目的はWHOの "Coronavirus disease (COVID-2019) situation reports" から各国別のConfirmedとDeathsの継時データを取り出すことである。一応形が整ったのであるが，WHOはpdfファイル中の表のデータ構造をコロナウィルスのようにどんどん変化させているので何だかすっきり行かない。結局，一部の古いものについては年度ごとにpdfから取り出したテキストデータ側を手動で調整する必要がある（これらすべてに対応させるほうが面倒なのだ）。
なお，pdfからのテキストファイルの取出にはpdftotextを使った。
（例）$ pdftotext -f 3 -l 6 20200312-sitrep-52-covid-19.pdf 52.txt
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
手元のMacBook Proの Catalina(10.15.3)に附属のperlは2013年バージョン5.18.4（This is perl 5, version 18, subversion 4 (v5.18.4) built for darwin-thread-multi-2level）だった。最新は，5.30.2のようで，Perl6はどうなっているのかと思えば，そうそう楽土になっていた。
なお，pdfからのテキストファイルの取出にはpdftotextを使った。
（例）$ pdftotext -f 3 -l 6 20200312-sitrep-52-covid-19.pdf 52.txt
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#!/usr/bin/perl
# /usr/local/bin/perl
# 03/17/2020 K. Koshigiri
# extract data from WHO covid-19 reports
# https://www.who.int/emergencies/diseases/
# novel-coronavirus-2019/situation-reports/
# usage:: ./corona.pl Japan ??.txt
$country = shift(@ARGV);
print("$country:\n");
foreach $file (@ARGV) {
open(IN, \$file) || die "\$file: ";
$file =~ /^([0-9]+).txt/;
\$no=\$1;
print($no,',');
while($line = <IN>) {
chomp($line);
if(\$line =~ /\$country\$/) {
$flg='y';
} elsif(\$flg eq 'y' && \$line =~ /([0-9]+)/) {
print("$1,");
} else {
$flg='n';
}
}
print("\n");
}
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -手元のMacBook Proの Catalina(10.15.3)に附属のperlは2013年バージョン5.18.4（This is perl 5, version 18, subversion 4 (v5.18.4) built for darwin-thread-multi-2level）だった。最新は，5.30.2のようで，Perl6はどうなっているのかと思えば，そうそう楽土になっていた。
2020年3月17日火曜日
感染症の数理シミュレーション（９）
感染症の数理シミュレーション（８）からの続き
湖北省（武漢以外）を説明するSIIDR2モデルができたので，これを用いて典型的な２つのイメージを定量化して表現してみたい。以下の計算ではモデルを単純化して特徴を表わすために，感染確率の$\beta$は一定だと仮定する。なお，次のパラメタは固定しておく（$\alpha_1=\dfrac{5}{0.80}, \alpha_2=\dfrac{5}{0.20}, \gamma_1=\dfrac{15}{0.96}, \gamma_2=\dfrac{15}{0.04}$）。また，重篤な患者に対して必要な設備（病床）数を，発症して隔離されている感染者数$u_3$の$\dfrac{1}{20}$と想定する。これは，厚生労働省の専門家会議における「新型コロナウイルス感染症の流行シナリオ」に準拠している。
Ａ．短期集団免疫シナリオ
SARS-CoV-2の感染者が完全に免疫を獲得できるのかどうかがよく知らないのだけれど（まれに再感染するのか，単なる変異したウイルスへの感染か），英国のジョンソンやドイツのメルケルが，国民の６〜７割が感染する可能性について言及しているようだ。果たしてそれが可能なのだろうか。SIIDR2モデルで$\beta=0.6,\ \nu=0.01,\ \lambda=10000,\ \tau=0,\ T=180$とする。$\lambda$を十分大きくとったのは$\bar{\beta} \approx \dfrac{8}{15}\beta$（一定）になるようにするためだ。
グラフの基点は新規感染者累計が人口の1ppmに達した時点であり，感染者が増大している多くの国では条件を満足している。日本に当てはめると，ピークは3〜4ヶ月後（5〜6月）$u_1$:未感染数，$u_2$:軽症感染数(含む未発症)，$u_3$:重症感染数（ピーク時400万人），$u_4$:死亡数（終期は60万人），$u_5$:免疫獲得回復数（終期は7000万人），$u_6$:重症感染数累計（終期は1400万人）である。したがってピーク時の重篤用必要設備（病床）数は20万床となり，完全に国内のキャパシティを越えてしまうものと思われる。当然オリンピックは不可能だ。
Ｂ．ピークシフトシナリオ
医療崩壊の本質は，PCR検査の数の問題でも，病床数の問題でもなく，重篤患者に対応するための設備と医者の数の問題であるといわれている。日本の厚生労働省がNHKを通じて広めている図も，ピークを低くしてその頂上の位置を先送りにする必要があるというものだった。SIIDR2モデルで$\beta=0.42,\ \nu=0.01,\ \lambda=10000,\ \tau=0,\ T=750$としたものが次の図である。図１とスケールのオーダーやファクターが違うことに注意しよう。
感染数のピークは10〜16ヶ月後，$u_2$:軽症感染数(含む未発症)（ピーク時50万人），$u_3$:重症感染数（ピーク時30万人），$u_4$:死亡数（終期は10万人），$u_6$:重症感染数累計（終期は300万人）である。したがってピーク時の重篤用必要設備（病床）数は1.5万床となり，かろうじて対応できるのかもしれないがよくわからない。大きな問題は２年にわたる長期の対応が必要なことである。ピークシフトを考える場合は，感染率$\beta$の抑制策との合わせ技が必須だ。
［１］ピークカット戦略（集団免疫戦略）地獄への道は善意で舗装されている（Sato Hiroshi）
湖北省（武漢以外）を説明するSIIDR2モデルができたので，これを用いて典型的な２つのイメージを定量化して表現してみたい。以下の計算ではモデルを単純化して特徴を表わすために，感染確率の$\beta$は一定だと仮定する。なお，次のパラメタは固定しておく（$\alpha_1=\dfrac{5}{0.80}, \alpha_2=\dfrac{5}{0.20}, \gamma_1=\dfrac{15}{0.96}, \gamma_2=\dfrac{15}{0.04}$）。また，重篤な患者に対して必要な設備（病床）数を，発症して隔離されている感染者数$u_3$の$\dfrac{1}{20}$と想定する。これは，厚生労働省の専門家会議における「新型コロナウイルス感染症の流行シナリオ」に準拠している。
Ａ．短期集団免疫シナリオ
SARS-CoV-2の感染者が完全に免疫を獲得できるのかどうかがよく知らないのだけれど（まれに再感染するのか，単なる変異したウイルスへの感染か），英国のジョンソンやドイツのメルケルが，国民の６〜７割が感染する可能性について言及しているようだ。果たしてそれが可能なのだろうか。SIIDR2モデルで$\beta=0.6,\ \nu=0.01,\ \lambda=10000,\ \tau=0,\ T=180$とする。$\lambda$を十分大きくとったのは$\bar{\beta} \approx \dfrac{8}{15}\beta$（一定）になるようにするためだ。
図１ 短期集団免疫シナリオ（１万人当り）
グラフの基点は新規感染者累計が人口の1ppmに達した時点であり，感染者が増大している多くの国では条件を満足している。日本に当てはめると，ピークは3〜4ヶ月後（5〜6月）$u_1$:未感染数，$u_2$:軽症感染数(含む未発症)，$u_3$:重症感染数（ピーク時400万人），$u_4$:死亡数（終期は60万人），$u_5$:免疫獲得回復数（終期は7000万人），$u_6$:重症感染数累計（終期は1400万人）である。したがってピーク時の重篤用必要設備（病床）数は20万床となり，完全に国内のキャパシティを越えてしまうものと思われる。当然オリンピックは不可能だ。
Ｂ．ピークシフトシナリオ
医療崩壊の本質は，PCR検査の数の問題でも，病床数の問題でもなく，重篤患者に対応するための設備と医者の数の問題であるといわれている。日本の厚生労働省がNHKを通じて広めている図も，ピークを低くしてその頂上の位置を先送りにする必要があるというものだった。SIIDR2モデルで$\beta=0.42,\ \nu=0.01,\ \lambda=10000,\ \tau=0,\ T=750$としたものが次の図である。図１とスケールのオーダーやファクターが違うことに注意しよう。
図２ ピークシフトシナリオ（１万人当り）
［１］ピークカット戦略（集団免疫戦略）地獄への道は善意で舗装されている（Sato Hiroshi）
2020年3月16日月曜日
日本の歴史
4月12日まで期間限定無料提供されている，小学館版学習漫画少年少女日本の歴史全24巻（22巻+付録2巻）をざっと読み終えた。第22巻の平成編はそれまでのものとトーンが違っており，全く面白くなかった。それ以外はとてもよかった。権力の移行過程に加えて，文化的なトピックスや民衆の動きなどが適当なバランスで記述されており，１巻の中を３〜４つのテーマに分けながら進めていくスタイルはとても理解しやすい。南京大虐殺や三光作戦や慰安婦などもちゃんと触れており，夏目漱石のところには米山保三郎も登場していた。第21巻あたりでは著作権の関係か背景がグレーに処理されていて新聞などの紙面？が再現されていないものがあるのがちょっと残念だった。
歴史といえば，大学時代に井上清（1915-2001）の岩波新書「日本の歴史（上・中・下）」を読んでいたことを思い出した。これもなかなか貴重な読書体験だった。史実の正確性や著者のマルクス主義的歴史観に対してはいろいろ批判はあるようだが，当時はとても新鮮に感じられた。
歴史といえば，小学校高学年のころか，父親が平泉澄（1895-1984）の「少年日本史」をどこからかもらってきて，これを読めと渡されたことがある。パラパラと見たところ直観的にこれはインチキだとすぐわかったのでほとんど読まなかった。父にこれ（皇国史観）はおかしいというと，叱られたような気がする。その頃は家にあった６巻ものの漫画ではない学研かどこかの日本の歴史を読んでいてこちらの方はまともでおもしろかった。
歴史といえば，大学時代に井上清（1915-2001）の岩波新書「日本の歴史（上・中・下）」を読んでいたことを思い出した。これもなかなか貴重な読書体験だった。史実の正確性や著者のマルクス主義的歴史観に対してはいろいろ批判はあるようだが，当時はとても新鮮に感じられた。
歴史といえば，小学校高学年のころか，父親が平泉澄（1895-1984）の「少年日本史」をどこからかもらってきて，これを読めと渡されたことがある。パラパラと見たところ直観的にこれはインチキだとすぐわかったのでほとんど読まなかった。父にこれ（皇国史観）はおかしいというと，叱られたような気がする。その頃は家にあった６巻ものの漫画ではない学研かどこかの日本の歴史を読んでいてこちらの方はまともでおもしろかった。
2020年3月15日日曜日
感染症の数理シミュレーション（８）
感染症の数理シミュレーション（７）からの続き
中間まとめ
いろいろと書き散らしてきたので，ここまでに得られた結果を整理してみる。
感染症の数理シミュレーション（３）において，感染対策時間因子を修正した図１のSIIDR2モデルで，湖北省（武漢以外）の新規感染者数累計と死亡数累計を説明できるパラメタセットを得ることができた。
\begin{equation}
\begin{aligned}
S: \dfrac{du_1}{dt} &= -\dfrac{\beta(t)}{n} u_1 u_2\\
I_1: \dfrac{du_2}{dt} &= \dfrac{\beta(t)}{n} u_1 u_2 -\dfrac{u_2}{\alpha_1} -\dfrac{u_2}{\alpha_2}\\
I_2: \dfrac{du_3}{dt} &= \dfrac{u_2}{\alpha_2} - \dfrac{u_3}{\gamma_1} -\dfrac{u_3}{\gamma_2}\\
D: \dfrac{du_4}{dt} &= \dfrac{u_3}{\gamma_2}\\
R: \dfrac{du_5}{dt} &= \dfrac{u_2}{\alpha_1} + \dfrac{u_3}{\gamma_1}\\
I_a: \dfrac{du_6}{dt} &= \dfrac{u_2}{\alpha_2}\\
\beta(t)&= \dfrac{\beta}{15}\{8+7 \tanh(-\frac{t-\tau}{\lambda})
\end{aligned}
\end{equation}
湖北省（武漢以外）のデータを再現するパラメタセットは以下のとおりであり，これを用いたシミュレーション結果は図２のように与えられる。
$\alpha_1=5/0.80, \alpha_2=5/0.20=25$, $\beta=0.915$, $\gamma_1=15/0.96, \gamma_2=15/0.04=375$, $\lambda=\tau=7, \nu=0.025$，$u_2$が軽症感染者数（発症無），$u_3$が重症感染者数（発症有），$\dfrac{u_2}{\alpha_2}$が新規重症感染者数$f(t)$，$\dfrac{u_3}{\gamma_2}$が新規死亡数$h(t)$，$u_4$が死亡数累計$i(t)=\int h(t) dt$，$u_6$が新規重症感染者累計$g(t)=\int f(t) dt$，である。WHOで報告されているデータ[1]では，感染者累計（Confirmed）が$g(t)$，死亡数累計（Deaths）が$i(t)$に対応している。
半分パラメタフィッティングなので，$\nu,\tau,\lambda$を調整して勝手に象の鼻を描いている感は否めないのが微妙なところである。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
using DifferentialEquations
using ParameterizedFunctions
using Plots; gr()
sky = @ode_def SIIDR2_model begin
du0 = 1 # u0:time
du1 = -β/15*(8+7*tanh(-(u0-τ)/λ))*u1*u2/n # u1:Noimmunity(Susceptible)
du2 = β/15*(8+7*tanh(-(u0-τ)/λ))*u1*u2/n -u2/α1 -u2/α2 # u2:Mild(Infected-a)
du3 = u2/α2 -u3/γ1 -u3/γ2 # u3:Serious(Infected-b)
du4 = u3/γ2 # u4:Dead
du5 = u2/α1 +u3/γ1 # u5:Recovered
du6 = u2/α2 # u6:Accumulated Infected-b
end n α1 α2 β γ1 γ2 λ τ
function epidm(β,ν,λ,τ,T)
n=10000.0 #total number of population
α1=5.0/0.80 #latent to recovery (days/%)
α2=5.0/0.20 #latent to onset (days/%)
#β=0.45 #infection rate (/day・person)
γ1=15.0/0.96 #onset to recovery (days/%)
γ2=15.0/0.04 #onset to death (days/%)
u0 = [0.0,n-11ν,4ν,2ν,0.0,5ν,ν] #initial values
p = (n,α1,α2,β,γ1,γ2,λ,τ) #parameters
tspan = (0.0,T) #time span in days
prob = ODEProblem(sky,u0,tspan,p)
sol = solve(prob)
return sol
end
β=0.915 #infection rate
ν=0.025 #inital value of accumulated infected-b
λ=7 #pandemic supression range (days)
τ=7 #pandemic supression start (day)
T=40 #period of simulation
xc=[0,3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48]
yc=[0.003,0.03,0.23,0.67,1.27,2.00,2.52,2.88,
3.44,3.58,3.71,3.69,3.71,3.71,3.71,3.71,3.71]
zc=[0.0,0.0,0.0,0.01,0.02,0.03,0.04,0.05,
0.07,0.09,0.10,0.11,0.11,0.12,0.12,0.13,0.13]
# kohoku-bukan model
# β=0.915,ν=0.025,λ=7,τ=7,α2=5.0/0.20,γ2=15.0/0.04
@time so=epidm(β,ν,λ,τ,T)
#plot(so,vars=(0,2))
plot(so,vars=(0,3))
plot!(so,vars=(0,4))
plot!(so,vars=(0,5))
plot!(so,vars=(0,6))
plot!(so,vars=(0,7))
plot!(xc,yc,st=:scatter)
plot!(xc,zc,st=:scatter)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
［１］Coronavirus disease (COVID-2019) situation reports
［２］新型環状病毒肺炎疫情動態（tencentニュースサイト）
［３］新型コロナウイルス感染症対策専門家会議(第5回)（3月2日）
［４］特設サイト：新型コロナウイルス（NHK）
［５］新型コロナウイルス感染速報（Su Wei）
［６］Coronavirus Disease (COVID-19) – Statistics and Research（Our World in DATA）
［７］Databrew's COVID-19 data explorer（Databrew）
［８］COVID-19情報共有（佐藤彰洋）
［９］時間遅れを考慮した確率 SIR モデルの安定性解析（石川昌明）
［１０］分布的時間遅れをもつ確率感染症モデルの安定性解析（石川昌明）
P. S. [6]のように，プロがシミュレーションをすると遅延SIRモデルみたいなことになるのか。感染率にステップ関数を導入しているので鋭い構造が出現しているようだ。遅延確率SIRモデルになるとちょっとついていけません(3/17/2020)。
感染症の数理シミュレーション（９）に続く
中間まとめ
いろいろと書き散らしてきたので，ここまでに得られた結果を整理してみる。
感染症の数理シミュレーション（３）において，感染対策時間因子を修正した図１のSIIDR2モデルで，湖北省（武漢以外）の新規感染者数累計と死亡数累計を説明できるパラメタセットを得ることができた。
図１ SIIDR2モデルの概念図，$n=u_1+u_2+u_3+u_4+u_5$, $p+q=1$
\begin{equation}
\begin{aligned}
S: \dfrac{du_1}{dt} &= -\dfrac{\beta(t)}{n} u_1 u_2\\
I_1: \dfrac{du_2}{dt} &= \dfrac{\beta(t)}{n} u_1 u_2 -\dfrac{u_2}{\alpha_1} -\dfrac{u_2}{\alpha_2}\\
I_2: \dfrac{du_3}{dt} &= \dfrac{u_2}{\alpha_2} - \dfrac{u_3}{\gamma_1} -\dfrac{u_3}{\gamma_2}\\
D: \dfrac{du_4}{dt} &= \dfrac{u_3}{\gamma_2}\\
R: \dfrac{du_5}{dt} &= \dfrac{u_2}{\alpha_1} + \dfrac{u_3}{\gamma_1}\\
I_a: \dfrac{du_6}{dt} &= \dfrac{u_2}{\alpha_2}\\
\beta(t)&= \dfrac{\beta}{15}\{8+7 \tanh(-\frac{t-\tau}{\lambda})
\end{aligned}
\end{equation}
湖北省（武漢以外）のデータを再現するパラメタセットは以下のとおりであり，これを用いたシミュレーション結果は図２のように与えられる。
$\alpha_1=5/0.80, \alpha_2=5/0.20=25$, $\beta=0.915$, $\gamma_1=15/0.96, \gamma_2=15/0.04=375$, $\lambda=\tau=7, \nu=0.025$，$u_2$が軽症感染者数（発症無），$u_3$が重症感染者数（発症有），$\dfrac{u_2}{\alpha_2}$が新規重症感染者数$f(t)$，$\dfrac{u_3}{\gamma_2}$が新規死亡数$h(t)$，$u_4$が死亡数累計$i(t)=\int h(t) dt$，$u_6$が新規重症感染者累計$g(t)=\int f(t) dt$，である。WHOで報告されているデータ[1]では，感染者累計（Confirmed）が$g(t)$，死亡数累計（Deaths）が$i(t)$に対応している。
図２ 湖北省（武漢以外）のSIIDR2モデルによるシミュレーション
（○はConfirmedとDeathsの観測値，WHOとtencentニュースサイトのデータ）
半分パラメタフィッティングなので，$\nu,\tau,\lambda$を調整して勝手に象の鼻を描いている感は否めないのが微妙なところである。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
using DifferentialEquations
using ParameterizedFunctions
using Plots; gr()
sky = @ode_def SIIDR2_model begin
du0 = 1 # u0:time
du1 = -β/15*(8+7*tanh(-(u0-τ)/λ))*u1*u2/n # u1:Noimmunity(Susceptible)
du2 = β/15*(8+7*tanh(-(u0-τ)/λ))*u1*u2/n -u2/α1 -u2/α2 # u2:Mild(Infected-a)
du3 = u2/α2 -u3/γ1 -u3/γ2 # u3:Serious(Infected-b)
du4 = u3/γ2 # u4:Dead
du5 = u2/α1 +u3/γ1 # u5:Recovered
du6 = u2/α2 # u6:Accumulated Infected-b
end n α1 α2 β γ1 γ2 λ τ
function epidm(β,ν,λ,τ,T)
n=10000.0 #total number of population
α1=5.0/0.80 #latent to recovery (days/%)
α2=5.0/0.20 #latent to onset (days/%)
#β=0.45 #infection rate (/day・person)
γ1=15.0/0.96 #onset to recovery (days/%)
γ2=15.0/0.04 #onset to death (days/%)
u0 = [0.0,n-11ν,4ν,2ν,0.0,5ν,ν] #initial values
p = (n,α1,α2,β,γ1,γ2,λ,τ) #parameters
tspan = (0.0,T) #time span in days
prob = ODEProblem(sky,u0,tspan,p)
sol = solve(prob)
return sol
end
β=0.915 #infection rate
ν=0.025 #inital value of accumulated infected-b
λ=7 #pandemic supression range (days)
τ=7 #pandemic supression start (day)
T=40 #period of simulation
xc=[0,3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48]
yc=[0.003,0.03,0.23,0.67,1.27,2.00,2.52,2.88,
3.44,3.58,3.71,3.69,3.71,3.71,3.71,3.71,3.71]
zc=[0.0,0.0,0.0,0.01,0.02,0.03,0.04,0.05,
0.07,0.09,0.10,0.11,0.11,0.12,0.12,0.13,0.13]
# kohoku-bukan model
# β=0.915,ν=0.025,λ=7,τ=7,α2=5.0/0.20,γ2=15.0/0.04
@time so=epidm(β,ν,λ,τ,T)
#plot(so,vars=(0,2))
plot(so,vars=(0,3))
plot!(so,vars=(0,4))
plot!(so,vars=(0,5))
plot!(so,vars=(0,6))
plot!(so,vars=(0,7))
plot!(xc,yc,st=:scatter)
plot!(xc,zc,st=:scatter)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
［２］新型環状病毒肺炎疫情動態（tencentニュースサイト）
［３］新型コロナウイルス感染症対策専門家会議(第5回)（3月2日）
［４］特設サイト：新型コロナウイルス（NHK）
［５］新型コロナウイルス感染速報（Su Wei）
［６］Coronavirus Disease (COVID-19) – Statistics and Research（Our World in DATA）
［７］Databrew's COVID-19 data explorer（Databrew）
［８］COVID-19情報共有（佐藤彰洋）
［９］時間遅れを考慮した確率 SIR モデルの安定性解析（石川昌明）
［１０］分布的時間遅れをもつ確率感染症モデルの安定性解析（石川昌明）
P. S. [6]のように，プロがシミュレーションをすると遅延SIRモデルみたいなことになるのか。感染率にステップ関数を導入しているので鋭い構造が出現しているようだ。遅延確率SIRモデルになるとちょっとついていけません(3/17/2020)。
感染症の数理シミュレーション（９）に続く
2020年3月14日土曜日
感染症の数理シミュレーション（７）
感染症の数理シミュレーション（６）からの続き
動機
岩田健太郎さんが，新型コロナウイルス感染症対策専門家会議(第5回)の参考資料「新型コロナウイルス感染症の流行シナリオ(2 月 29 日時点)」を紹介していた。そこでは，基本再生産数$R_o=1.7$として，発病間隔を4日としていた。$\beta=1.7/4=0.42$であり，これらはこれまでの想定と矛盾しない。ところが入院期間や重症期間を15日と設定しており，これまでの$\gamma=5$という想定とは異なっている。そこで，前回の考察を見直したところ，導出根拠が不正確な部分があったので再検討する。
方針
中国の湖北省（武漢以外，人口4800万人）の感染データは感染者数が最も大きくほぼ収束しているデータである。新規確定感染者数 , 累計 $f(t) , g(t)=\int f(t) dt$と新規死亡数 , 累計 $h(t) , i(t)=\int h(t) dt$ が観測値として得られており，武漢のような特異性（突出した致死率）をもたない。そこで，これのデータを簡単な関数で近似して，SIIDR2モデルにあてはめることにより，モデルのパラメタである$\beta, \gamma_1, \gamma_2$を観測値から推定しようというものである。なお，前回の議論から$\alpha_1 = \dfrac{5}{0.8} = 6.25,\ \alpha_2 = \dfrac{5}{0.2} = 25,\ \dfrac{\alpha_2}{\alpha_1}=4$は前提とする。
SIIDR2モデル式（$u_1 \approx n$, $\beta$=一定とした場合）
\begin{equation}
\begin{aligned}
u_1 &= n - \int \beta u_2 dt = n - 25 \beta g\\
u_2 &= 25 \beta g -5 g\\
u_3 &= \gamma_2 h\\
u_4 &= i\\
u_5 &= 4 g + \dfrac{\gamma_2}{\gamma_1} i \\
u_6 &= g
\end{aligned}
\end{equation}
これを$u_1+u_2+u_3+u_4+u_5=n$に代入して次式が得られる。
\begin{equation}
g(t) = \gamma_2 h(t) + (1+ \dfrac{\gamma_2}{\gamma_1}) i(t)
\end{equation}
この左辺と右辺がほぼ等しくなる$\gamma_1, \gamma_2$を探す。
観測値の近似式
湖北省（武漢以外）の観測値（1万人当りに換算）を再現するものとして前回と同様に次の近似関数を採用する。ただし，$h(t)$は前回から修正した。
\begin{equation}
\begin{aligned}
f(t) = t^4 \exp(-t/3) / 1600\\
g(t) = \int f(t) dt ; g(0)=0\\
h(t) = 1.4 t \exp(-(t - 20)^2/12^2) /4800\\
i(t) = \int h(t) dt ; i(0)=0
\end{aligned}
\end{equation}
$h(t), i(t)$を24倍($\gamma_2/\gamma_1$)したグラフを$f(t), g(t)$と重ねたものを図１に示す。
最適値の検討
$r(t) = \{ \gamma_2 h(t) + (1+ \dfrac{\gamma_2}{\gamma_1}) i(t) \}/g(t)$を$t=10$から$t=30$の範囲でプロットして$r(t)$の値がほぼフラットであり平均して1前後になる$\gamma_1, \gamma_2$を探す。上記専門家会議の議論を踏まえて，$\frac{1}{\gamma}=\frac{1}{\gamma_1}+\frac{1}{\gamma_2}=\frac{1}{15}$とする。なお，これまでの$\gamma=5$を用いると$r(t)$が0.6あたりまで減ってしまう。結局，
\begin{equation}
\gamma_1 = \dfrac{15}{0.96}, \quad \gamma_2 = \dfrac{15}{0.04}
\end{equation}
このとき，$r(t)$の平均値として$\dfrac{1}{20} \int_{10}^{30} r(t) dt = 1.05$が得られるので，条件を満足している。
まとめ
湖北省（武漢以外）の1/20/2020から2/29/2020のデータを再現するSIIDR2モデルのパラメタとしては，次のセットが推奨される。
\begin{equation}
\alpha_1 = \dfrac{5}{0.8} ,\quad \alpha_2 = \dfrac{5}{0.2} ,\quad \bar{\beta}=0.4 \sim 0.5 ,\quad \gamma_1 = \dfrac{15}{0.96} ,\quad \gamma_2 = \dfrac{15}{0.04}
\end{equation}
付録 初期値$\nu$の取り扱い
これまでは，$u_2$の初期値を$\nu$としてきた。ところで，各国の感染者比（１）と各国の感染者比（２）から新規感染者数累計が人口の1ppmを越えた$t_0$を始点とするモデル化が必要があることがわかった。そこで，$t_0$における初期条件を再度検討する。はじめに$h(t_0)=i(t_0)=0$と近似した上で，上記関係式より，$u_2(t_0)=25 f(t_0), u_5(t_0)=4 g(t_0) , u_6(t_0)=g(t_0)$となる。ここで，$g(t_0)=\nu$と再定義する。$\nu$=1ppmが始点になる。
ところで，$f(t)=k t^m$と仮定すると
\begin{equation}
\begin{aligned}
g(t) &=\int f(t) dt = \dfrac{k t^{m+1}}{m+1}=\dfrac{t f(t)}{m+1}\\
\therefore f(t_0) &= \dfrac{m+1}{t_0} g(t_0) = \dfrac{m+1}{t_0} \nu
\end{aligned}
\end{equation}
したがって，$u_2(t_0)=\dfrac{25}{t_0}(m+1) \nu \approx 4\nu$となる。ただし，日本の1ppmに至るまでの観測値から，$t_0 \approx 25, m \approx 3$とおいた。
追加（3/17/2020）
次に，$u_3(t_0)$について考える。これは次の微分方程式を満たす。
\begin{equation}
\dfrac{du_3}{dt}=f -\dfrac{u_3}{\gamma}
\end{equation}
ここで，$u_3(t)=k t$と仮定して代入すると，$k = f -\dfrac{k t}{\gamma}$である。そこで，
$k=\dfrac{f(t_0)}{1+t_0/\gamma}$となる。上の結果から，$u_2(t_0)=25 f(t_0)$，したがって，次式が成り立つ。
\begin{equation}
u_3(t_0)= \dfrac{f(t_0)}{1+t_0/\gamma} t_0 = u_2(t_0) \dfrac{25}{t_0 (1+t_0/\gamma)}
\approx u_2(t_0) \dfrac{1}{(1+t_0/\gamma)} \sim \dfrac{1}{2} u_2(t_0) = 2\nu
\end{equation}
さらに，$h(t_0)=i(t_0)=0$という条件をゆるめてみる。$u_3=\dfrac{u_2}{2}=375 h$より，$\int u_2 dt = 750 \int h dt = 750 i$，一方，$g = \int f dt = \dfrac{1}{25} \int u_2 dt$。したがって，$i = \dfrac{1}{30} g$となる。これから，$u_5(t_0)=4 g(t_0) + \dfrac{24}{30} g(t_0) \approx 5 \nu$となる。なお，$u_4(t_0)=i(t_0)= \dfrac{g(t_0)}{30} = \dfrac{\nu}{30} $なので，これはゼロのままとする。
この結果$u_1(t_0)=n-u_2(t_0)-u_3(t_0)-u_4(t_0)-u_5(t_0)=n-4\nu-2\nu-5\nu=n-11\nu$となり，新規感染者数累計が1ppmになる時点には，約11ppmの軽症感染者・重症感染者・死亡者・免疫獲得回復者が存在することになる。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
まとめると，$\nu$=1ppm時点の初期値を次のように設定する。
\begin{equation}
u_1(t_0)=n-11\nu, \quad u_2(t_0)=4\nu, \quad u_3(t_0)=2\nu, \quad u_4(t_0)=0, \quad u_5(t_0)=5\nu, \quad u_6(t_0)=\nu
\end{equation}
感染症の数理シミュレーション（８）に続く
動機
岩田健太郎さんが，新型コロナウイルス感染症対策専門家会議(第5回)の参考資料「新型コロナウイルス感染症の流行シナリオ(2 月 29 日時点)」を紹介していた。そこでは，基本再生産数$R_o=1.7$として，発病間隔を4日としていた。$\beta=1.7/4=0.42$であり，これらはこれまでの想定と矛盾しない。ところが入院期間や重症期間を15日と設定しており，これまでの$\gamma=5$という想定とは異なっている。そこで，前回の考察を見直したところ，導出根拠が不正確な部分があったので再検討する。
方針
中国の湖北省（武漢以外，人口4800万人）の感染データは感染者数が最も大きくほぼ収束しているデータである。新規確定感染者数 , 累計 $f(t) , g(t)=\int f(t) dt$と新規死亡数 , 累計 $h(t) , i(t)=\int h(t) dt$ が観測値として得られており，武漢のような特異性（突出した致死率）をもたない。そこで，これのデータを簡単な関数で近似して，SIIDR2モデルにあてはめることにより，モデルのパラメタである$\beta, \gamma_1, \gamma_2$を観測値から推定しようというものである。なお，前回の議論から$\alpha_1 = \dfrac{5}{0.8} = 6.25,\ \alpha_2 = \dfrac{5}{0.2} = 25,\ \dfrac{\alpha_2}{\alpha_1}=4$は前提とする。
SIIDR2モデル式（$u_1 \approx n$, $\beta$=一定とした場合）
\begin{equation}
\begin{aligned}
u_1 &= n - \int \beta u_2 dt = n - 25 \beta g\\
u_2 &= 25 \beta g -5 g\\
u_3 &= \gamma_2 h\\
u_4 &= i\\
u_5 &= 4 g + \dfrac{\gamma_2}{\gamma_1} i \\
u_6 &= g
\end{aligned}
\end{equation}
これを$u_1+u_2+u_3+u_4+u_5=n$に代入して次式が得られる。
\begin{equation}
g(t) = \gamma_2 h(t) + (1+ \dfrac{\gamma_2}{\gamma_1}) i(t)
\end{equation}
この左辺と右辺がほぼ等しくなる$\gamma_1, \gamma_2$を探す。
観測値の近似式
湖北省（武漢以外）の観測値（1万人当りに換算）を再現するものとして前回と同様に次の近似関数を採用する。ただし，$h(t)$は前回から修正した。
\begin{equation}
\begin{aligned}
f(t) = t^4 \exp(-t/3) / 1600\\
g(t) = \int f(t) dt ; g(0)=0\\
h(t) = 1.4 t \exp(-(t - 20)^2/12^2) /4800\\
i(t) = \int h(t) dt ; i(0)=0
\end{aligned}
\end{equation}
$h(t), i(t)$を24倍($\gamma_2/\gamma_1$)したグラフを$f(t), g(t)$と重ねたものを図１に示す。
図１ ピークと変曲点の時点が$t=12$の$f(t),\ g(t)$及び$t=24$の$h(t),\ i(t)$
最適値の検討
$r(t) = \{ \gamma_2 h(t) + (1+ \dfrac{\gamma_2}{\gamma_1}) i(t) \}/g(t)$を$t=10$から$t=30$の範囲でプロットして$r(t)$の値がほぼフラットであり平均して1前後になる$\gamma_1, \gamma_2$を探す。上記専門家会議の議論を踏まえて，$\frac{1}{\gamma}=\frac{1}{\gamma_1}+\frac{1}{\gamma_2}=\frac{1}{15}$とする。なお，これまでの$\gamma=5$を用いると$r(t)$が0.6あたりまで減ってしまう。結局，
\begin{equation}
\gamma_1 = \dfrac{15}{0.96}, \quad \gamma_2 = \dfrac{15}{0.04}
\end{equation}
このとき，$r(t)$の平均値として$\dfrac{1}{20} \int_{10}^{30} r(t) dt = 1.05$が得られるので，条件を満足している。
図２ $t=10$から$t=30$における$r(t)$の振る舞い
また前回の，$\beta(t)=\dfrac{df}{dt} / f + \dfrac{1}{\alpha}$についての議論はそのままであり変更されていないことに注意する。SIIDR2モデルを採用すれば，武漢（湖北以外）では$\beta(t)$を急速に小さくしたことが考えられる。なお，基本再生産率$R_o$と$\alpha$の値から$\beta$を一定とした場合の平均値は$\bar{\beta}=0.4 \sim 0.5$である。
また前回の，$\beta(t)=\dfrac{df}{dt} / f + \dfrac{1}{\alpha}$についての議論はそのままであり変更されていないことに注意する。SIIDR2モデルを採用すれば，武漢（湖北以外）では$\beta(t)$を急速に小さくしたことが考えられる。なお，基本再生産率$R_o$と$\alpha$の値から$\beta$を一定とした場合の平均値は$\bar{\beta}=0.4 \sim 0.5$である。
まとめ
湖北省（武漢以外）の1/20/2020から2/29/2020のデータを再現するSIIDR2モデルのパラメタとしては，次のセットが推奨される。
\begin{equation}
\alpha_1 = \dfrac{5}{0.8} ,\quad \alpha_2 = \dfrac{5}{0.2} ,\quad \bar{\beta}=0.4 \sim 0.5 ,\quad \gamma_1 = \dfrac{15}{0.96} ,\quad \gamma_2 = \dfrac{15}{0.04}
\end{equation}
付録 初期値$\nu$の取り扱い
これまでは，$u_2$の初期値を$\nu$としてきた。ところで，各国の感染者比（１）と各国の感染者比（２）から新規感染者数累計が人口の1ppmを越えた$t_0$を始点とするモデル化が必要があることがわかった。そこで，$t_0$における初期条件を再度検討する。はじめに$h(t_0)=i(t_0)=0$と近似した上で，上記関係式より，$u_2(t_0)=25 f(t_0), u_5(t_0)=4 g(t_0) , u_6(t_0)=g(t_0)$となる。ここで，$g(t_0)=\nu$と再定義する。$\nu$=1ppmが始点になる。
ところで，$f(t)=k t^m$と仮定すると
\begin{equation}
\begin{aligned}
g(t) &=\int f(t) dt = \dfrac{k t^{m+1}}{m+1}=\dfrac{t f(t)}{m+1}\\
\therefore f(t_0) &= \dfrac{m+1}{t_0} g(t_0) = \dfrac{m+1}{t_0} \nu
\end{aligned}
\end{equation}
したがって，$u_2(t_0)=\dfrac{25}{t_0}(m+1) \nu \approx 4\nu$となる。ただし，日本の1ppmに至るまでの観測値から，$t_0 \approx 25, m \approx 3$とおいた。
追加（3/17/2020）
次に，$u_3(t_0)$について考える。これは次の微分方程式を満たす。
\begin{equation}
\dfrac{du_3}{dt}=f -\dfrac{u_3}{\gamma}
\end{equation}
ここで，$u_3(t)=k t$と仮定して代入すると，$k = f -\dfrac{k t}{\gamma}$である。そこで，
$k=\dfrac{f(t_0)}{1+t_0/\gamma}$となる。上の結果から，$u_2(t_0)=25 f(t_0)$，したがって，次式が成り立つ。
\begin{equation}
u_3(t_0)= \dfrac{f(t_0)}{1+t_0/\gamma} t_0 = u_2(t_0) \dfrac{25}{t_0 (1+t_0/\gamma)}
\approx u_2(t_0) \dfrac{1}{(1+t_0/\gamma)} \sim \dfrac{1}{2} u_2(t_0) = 2\nu
\end{equation}
さらに，$h(t_0)=i(t_0)=0$という条件をゆるめてみる。$u_3=\dfrac{u_2}{2}=375 h$より，$\int u_2 dt = 750 \int h dt = 750 i$，一方，$g = \int f dt = \dfrac{1}{25} \int u_2 dt$。したがって，$i = \dfrac{1}{30} g$となる。これから，$u_5(t_0)=4 g(t_0) + \dfrac{24}{30} g(t_0) \approx 5 \nu$となる。なお，$u_4(t_0)=i(t_0)= \dfrac{g(t_0)}{30} = \dfrac{\nu}{30} $なので，これはゼロのままとする。
この結果$u_1(t_0)=n-u_2(t_0)-u_3(t_0)-u_4(t_0)-u_5(t_0)=n-4\nu-2\nu-5\nu=n-11\nu$となり，新規感染者数累計が1ppmになる時点には，約11ppmの軽症感染者・重症感染者・死亡者・免疫獲得回復者が存在することになる。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
まとめると，$\nu$=1ppm時点の初期値を次のように設定する。
\begin{equation}
u_1(t_0)=n-11\nu, \quad u_2(t_0)=4\nu, \quad u_3(t_0)=2\nu, \quad u_4(t_0)=0, \quad u_5(t_0)=5\nu, \quad u_6(t_0)=\nu
\end{equation}
感染症の数理シミュレーション（８）に続く
2020年3月13日金曜日
小学館版学習漫画少年少女日本の歴史
新型コロナウイル感染症のため全国の学校が休校になり，さまざまな対応措置が講じられている。そのひとつが子ども向け学習・娯楽コンテンツの期間限定無料提供である。例えば，小学館版学習漫画少年少女日本の歴史全24巻（22巻+付録2巻）が2020年3月11日から4月12日まで無料公開されている。
さっそく読み始めました。いやーおもしろくてためになる。高校のときは理数科で世界史と日本史は片方のみが選択だった。それならば世界史でしょうということで，受験に不利かもしれない世界史を選択したため，日本史は小学校の知識でとどまっている。中学校社会の歴史は，途中で交代もあったかもしれないなんだかもうひとつの先生で（地理の米山先生は若くて授業に集中できたのに），試験の成績も全く振るわなかったのだ。
小学校で歴史を学ぶのは６年生である。５年生が地理。その泉野小学校６年の担任の前多光子先生が，４月の社会の授業の最初に１時間かけて何のために歴史を学ぶのかということをひたすら子どもたちに問い続けた。みんなで考え続ける時間がとても長かった。で，結論は「温故知新（ふるきをたずねてあたらしきをしる）」ということだった。学校生活で数少ない印象的だった授業の一つである。
さて，学習漫画のほうに話を戻す。
わかったことは，日本の政治は悪者が担ってきたということだ。どおりで最近の政権支持率が40%で高止まりしているわけだ。子どものときからみんながこの37年に渡るベストセラーの学習漫画で育っている。そんなわけで，今の50歳以下の世代は我々と違って政治の本質を見抜いてしまっており，なおかつこれが歴史というものだとあきらめて（達観して）いるのかもしれない。その土壌の上に新自由主義が猖獗するのもしょうがない。
さっそく読み始めました。いやーおもしろくてためになる。高校のときは理数科で世界史と日本史は片方のみが選択だった。それならば世界史でしょうということで，受験に不利かもしれない世界史を選択したため，日本史は小学校の知識でとどまっている。中学校社会の歴史は，途中で交代もあったかもしれないなんだかもうひとつの先生で（地理の米山先生は若くて授業に集中できたのに），試験の成績も全く振るわなかったのだ。
小学校で歴史を学ぶのは６年生である。５年生が地理。その泉野小学校６年の担任の前多光子先生が，４月の社会の授業の最初に１時間かけて何のために歴史を学ぶのかということをひたすら子どもたちに問い続けた。みんなで考え続ける時間がとても長かった。で，結論は「温故知新（ふるきをたずねてあたらしきをしる）」ということだった。学校生活で数少ない印象的だった授業の一つである。
さて，学習漫画のほうに話を戻す。
わかったことは，日本の政治は悪者が担ってきたということだ。どおりで最近の政権支持率が40%で高止まりしているわけだ。子どものときからみんながこの37年に渡るベストセラーの学習漫画で育っている。そんなわけで，今の50歳以下の世代は我々と違って政治の本質を見抜いてしまっており，なおかつこれが歴史というものだとあきらめて（達観して）いるのかもしれない。その土壌の上に新自由主義が猖獗するのもしょうがない。
2020年3月12日木曜日
各国の感染者比（２）
各国の感染者比（１）からの続き
新型コロナウイルス感染症の話。
感染者数が一定の数あるいは人口比になった日時とその値を共通の原点として，各国の感染者数や人口比を対数グラフに表して観察する方法についてSNS上であれこれ話題になっている。最初の指摘は「他国と比べて日本だけが特異な振る舞いを示しているのは人為的な要因（検査の絞り込み等）の結果ではないか」というものであった。
これに対して，感染者数基準のプロットではほとんど重なって見えた他国・地域については，人口比にすれば分離して見えることを指摘した。それにもかかわらず，中国，湖北省（武漢以外），韓国，日本，イタリア，イランの比較では日本だけが異なった傾向を見せた（緩慢な指数関数型増加）。一方，症例の少ない多くの国を含めて比較すれば，日本だけが特異なのではないというクレームがついていた。ただし，無原則にあれもこれも混ぜてしまうのは問題の本質を見失う危険性があるので正しくない方法論だ。
そこで，アジア太平洋地域で一定の感染者が報告されている国々と日本，韓国について前回と同様の分析を行った。対象は，韓国，日本，台湾，香港，シンガポール，オーストラリアの６カ国とした。なお，今回の追加データは，ジョンズ・ホプキンズ大学のアーカイブから取得している。各国の人口，1ppm達成日，3/9時点での新規確定感染数累計（Confirmed）/人口（Population）は，オーストラリア（2460万人，2月29日，3.7ppm），日本（12600万人，2月22日，4.1ppm），韓国（5180万人，2月19日，140ppm），シンガポール（560万人，1月27日，15ppm），台湾（2360万人，2月20日，1.9ppm），香港（740万人，1月26日，8.1ppm）である。
新型コロナウイルス感染症の話。
感染者数が一定の数あるいは人口比になった日時とその値を共通の原点として，各国の感染者数や人口比を対数グラフに表して観察する方法についてSNS上であれこれ話題になっている。最初の指摘は「他国と比べて日本だけが特異な振る舞いを示しているのは人為的な要因（検査の絞り込み等）の結果ではないか」というものであった。
これに対して，感染者数基準のプロットではほとんど重なって見えた他国・地域については，人口比にすれば分離して見えることを指摘した。それにもかかわらず，中国，湖北省（武漢以外），韓国，日本，イタリア，イランの比較では日本だけが異なった傾向を見せた（緩慢な指数関数型増加）。一方，症例の少ない多くの国を含めて比較すれば，日本だけが特異なのではないというクレームがついていた。ただし，無原則にあれもこれも混ぜてしまうのは問題の本質を見失う危険性があるので正しくない方法論だ。
そこで，アジア太平洋地域で一定の感染者が報告されている国々と日本，韓国について前回と同様の分析を行った。対象は，韓国，日本，台湾，香港，シンガポール，オーストラリアの６カ国とした。なお，今回の追加データは，ジョンズ・ホプキンズ大学のアーカイブから取得している。各国の人口，1ppm達成日，3/9時点での新規確定感染数累計（Confirmed）/人口（Population）は，オーストラリア（2460万人，2月29日，3.7ppm），日本（12600万人，2月22日，4.1ppm），韓国（5180万人，2月19日，140ppm），シンガポール（560万人，1月27日，15ppm），台湾（2360万人，2月20日，1.9ppm），香港（740万人，1月26日，8.1ppm）である。
図 各国の感染者比（基準日は1ppm時点，縦軸はppmの常用対数）
（１）韓国と台湾は，それぞれ異なった水準ではあるが収束に向かう様子が窺える。
（２）欧州の感染者急増国の多くは韓国と相似形の振る舞いをしている。
（３）取り上げたアジア太平洋諸国（韓国，台湾以外，日本を含む）は対数グラフでおおよそなだらかな１次関数で増加している（すなわち，底の小さな指数関数的増加）。したがって，まだ収束の様子はみられない。
（４）香港とシンガポールは１次関数のスロープが原点から２週目以降でよりゆるやかに変化している。なんらかの抑制策が効果を発揮したのかもしれない。シンガポールは0日から18日まで平均15.5%，18日から42日まで平均3.5%で増加，香港は0日から15日まで平均11.0%，15日から43日まで平均4%で増加している。なお，日本は，0日から16日まで平均9%とその中間になっている。
2020年3月11日水曜日
ワン・モア・ヌーク：藤井太洋（２）
3月14日のエピローグも含めて，3月10日で丁度ワン・モア・ヌークを読了した。たいへん知的でスリリングで映画にしたらおもしろいなあと思いながらスイスイと読み進むことができた。テーマとしては9年目を迎える東日本大震災に伴う福島第一原子力発電所事故にピタリと照準をあわせてきた。原子爆弾に関するディテールとこれにイスラム国のテロがうまく調和した内容になっている。もちろん本当の専門家から見ればなんだこれはということになるのだろうけれど，我々素人からすれば十分にリアルに感じられる。
残念ながら著者の福島原発被害者の不安に寄り添おうとする視点には共感できなかった。これはほとんど例のアンチニセ科学グループの思想と同じだったから。具体的には65ページの「テレビ，新聞，雑誌やインターネットを通して死の恐怖を感じさせる言葉をひっきりなしに浴びせられた中で・・・」から67ページの「実はこの程度の被爆でがんの発症例が増えたという疫学的な証拠はない」あたりまで。著者の説明はほとんど間違ってはいないかもしれないが，同時に警鐘をならす人々をすべて飲み込んでしまう効果をもたらしている。そして今も処理水という名のトリチウム＋未処理放射性物質の海洋投棄や小児甲状腺がんの診断をめぐってその闘いは続いている。キクマコよりおしどりマコが，早野龍五より牧野淳一郎が信頼できる。
ワン・モア・ヌークはおもしろかったのでこれは評価したい。ただ，こころに残るかというとどうだろうか。この程度の科学的な精度を持ちながらより深い思索をもたらすものを読んでみたい。安部公房は残念ながらもういない。
残念ながら著者の福島原発被害者の不安に寄り添おうとする視点には共感できなかった。これはほとんど例のアンチニセ科学グループの思想と同じだったから。具体的には65ページの「テレビ，新聞，雑誌やインターネットを通して死の恐怖を感じさせる言葉をひっきりなしに浴びせられた中で・・・」から67ページの「実はこの程度の被爆でがんの発症例が増えたという疫学的な証拠はない」あたりまで。著者の説明はほとんど間違ってはいないかもしれないが，同時に警鐘をならす人々をすべて飲み込んでしまう効果をもたらしている。そして今も処理水という名のトリチウム＋未処理放射性物質の海洋投棄や小児甲状腺がんの診断をめぐってその闘いは続いている。キクマコよりおしどりマコが，早野龍五より牧野淳一郎が信頼できる。
ワン・モア・ヌークはおもしろかったのでこれは評価したい。ただ，こころに残るかというとどうだろうか。この程度の科学的な精度を持ちながらより深い思索をもたらすものを読んでみたい。安部公房は残念ながらもういない。
2020年3月10日火曜日
各国の感染者比（１）
各国の新規確定感染者数の累計（Confirmed）が100名または200名となった時点を１としてその前後の感染者数累計を対数プロットしたものがtwitter上でいくつか紹介され，日本だけが異なった振る舞いを示していることが指摘されている。日本の感染者数の統計が不自然ではないかということはこれまでも話題になってきた。
ところで国によって人口は大きく異なるのにもかかわらず基準を100名などの固定した人数にするのはどうかと思う。そこで，新規確定感染者数の累計（Confirmed）を各国／地域の人口で割った感染者比が1ppm（百万分の一）となるところを基準にしてそれらの時間経過を比較してみよう。データはWHOのCoronavirus disease (COVID-2019) situation reportsを用いる。感染者数が多い次の６つの国や地域（日本を含む），(1) 中国（13億9500万人）(2) 湖北省（武漢以外）（4800万人），(3) 韓国（5200万人），(4) 日本（1億2700万人），(5) イタリア（6000万人），(6) イラン（8100万人）を選んだ。
それぞれの国の現在（グラフの終端 3/9/2020）の感染者比と感染者が1ppmになる日時はおおむね次のようになっている。(1) 中国（58ppm，1月25日），(2) 湖北省（武漢以外）（370ppm，1月21日*），(3) 韓国（140ppm，2月19日），(4) 日本（3.8ppm，2月23日），(5) イタリア（120ppm，2月27日），(6) イラン（81ppm，2月25日）その結果は次の図で示される。
（１）各国の1ppm以下のデータはバラつきが大きく，感染症の数理シミュレーションの基礎データとして用いるには不適当である。今後，さらに各国の比較を行う場合は感染者比1ppm以上のデータで議論することが望ましい。
（２）感染者数が多い国おける1ppm以上の感染者比の増加カーブの様子はこれまでのところよく似ている。韓国では中国と同様の飽和が始まっているようだ。韓国，イタリア，イランは現在まで湖北省と中国全土に挟まれたデータ領域で推移している。
（３）感染源近傍の湖北省（武漢以外）に比べ中国全土の感染者比の飽和値は約1/6となっている。中国における湖北省の封じ込めが成功しているのだと考えられる。各国や地域の社会構造（政治・経済・文化）や感染対策効果に依存して感染者比の飽和値が定まっていく。
（４）この中では日本のデータだけが他の国や地域と異なった振る舞いをしている。しかしまだ飽和に達している様子は見られない。なお，1ppmに達していない国のデータを混ぜて議論するべきではない。人口密度の低いオーストラリアなども日本と同様の傾向ではないかとの指摘があり，引き続き原因の分析が必要である。
（注１）日本のデータはWHOのデータの扱いと同様で，ダイヤモンド・プリンセス号の部分は除かれている。
（注２）湖北省（武漢以外）はWHOのデータからは直接得られないので，湖北省（武漢以外）のデータによる考察とモデル計算に基づいている。
（注３）WHOの日本のデータ2月5日分には奥村晴彦さんのCOVID-19が指摘するような誤記があるので修正した。なお，中国のデータも奥村さんのcsvファイルを利用させていただいた。
（注４）中国全土のデータの2月17日の段差は集計方法の変更によるものである。
（注５）＊1/30 -> 1/21 に訂正（グラフにはまだ反映されていない 3/15/2020）
［１］感染症の数理シミュレーション（６）
各国の感染者比（２）に続く
ところで国によって人口は大きく異なるのにもかかわらず基準を100名などの固定した人数にするのはどうかと思う。そこで，新規確定感染者数の累計（Confirmed）を各国／地域の人口で割った感染者比が1ppm（百万分の一）となるところを基準にしてそれらの時間経過を比較してみよう。データはWHOのCoronavirus disease (COVID-2019) situation reportsを用いる。感染者数が多い次の６つの国や地域（日本を含む），(1) 中国（13億9500万人）(2) 湖北省（武漢以外）（4800万人），(3) 韓国（5200万人），(4) 日本（1億2700万人），(5) イタリア（6000万人），(6) イラン（8100万人）を選んだ。
それぞれの国の現在（グラフの終端 3/9/2020）の感染者比と感染者が1ppmになる日時はおおむね次のようになっている。(1) 中国（58ppm，1月25日），(2) 湖北省（武漢以外）（370ppm，1月21日*），(3) 韓国（140ppm，2月19日），(4) 日本（3.8ppm，2月23日），(5) イタリア（120ppm，2月27日），(6) イラン（81ppm，2月25日）その結果は次の図で示される。
図 各国の感染者比（基準日は1ppm時点，縦軸はppmの常用対数）
（１）各国の1ppm以下のデータはバラつきが大きく，感染症の数理シミュレーションの基礎データとして用いるには不適当である。今後，さらに各国の比較を行う場合は感染者比1ppm以上のデータで議論することが望ましい。
（３）感染源近傍の湖北省（武漢以外）に比べ中国全土の感染者比の飽和値は約1/6となっている。中国における湖北省の封じ込めが成功しているのだと考えられる。各国や地域の社会構造（政治・経済・文化）や感染対策効果に依存して感染者比の飽和値が定まっていく。
（４）この中では日本のデータだけが他の国や地域と異なった振る舞いをしている。しかしまだ飽和に達している様子は見られない。なお，1ppmに達していない国のデータを混ぜて議論するべきではない。人口密度の低いオーストラリアなども日本と同様の傾向ではないかとの指摘があり，引き続き原因の分析が必要である。
（注１）日本のデータはWHOのデータの扱いと同様で，ダイヤモンド・プリンセス号の部分は除かれている。
（注２）湖北省（武漢以外）はWHOのデータからは直接得られないので，湖北省（武漢以外）のデータによる考察とモデル計算に基づいている。
（注３）WHOの日本のデータ2月5日分には奥村晴彦さんのCOVID-19が指摘するような誤記があるので修正した。なお，中国のデータも奥村さんのcsvファイルを利用させていただいた。
（注４）中国全土のデータの2月17日の段差は集計方法の変更によるものである。
（注５）＊1/30 -> 1/21 に訂正（グラフにはまだ反映されていない 3/15/2020）
［１］感染症の数理シミュレーション（６）
各国の感染者比（２）に続く
登録:
投稿 (Atom)