感染症の数理シミュレーション（８）

感染症の数理シミュレーション（７）からの続き

中間まとめ
いろいろと書き散らしてきたので，ここまでに得られた結果を整理してみる。

感染症の数理シミュレーション（３）において，感染対策時間因子を修正した図１の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)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

［１］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)。


感染症の数理シミュレーション（９）に続く

感染症の数理シミュレーション（７）

感染症の数理シミュレーション（６）からの続き

動機
岩田健太郎さんが，新型コロナウイルス感染症対策専門家会議(第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$である。

まとめ
湖北省（武漢以外）の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}


感染症の数理シミュレーション（８）に続く

小学館版学習漫画少年少女日本の歴史

新型コロナウイル感染症のため全国の学校が休校になり，さまざまな対応措置が講じられている。そのひとつが子ども向け学習・娯楽コンテンツの期間限定無料提供である。例えば，小学館版学習漫画少年少女日本の歴史全24巻（22巻+付録2巻）が2020年3月11日から4月12日まで無料公開されている。

さっそく読み始めました。いやーおもしろくてためになる。高校のときは理数科で世界史と日本史は片方のみが選択だった。それならば世界史でしょうということで，受験に不利かもしれない世界史を選択したため，日本史は小学校の知識でとどまっている。中学校社会の歴史は，途中で交代もあったかもしれないなんだかもうひとつの先生で（地理の米山先生は若くて授業に集中できたのに），試験の成績も全く振るわなかったのだ。

小学校で歴史を学ぶのは６年生である。５年生が地理。その泉野小学校６年の担任の前多光子先生が，４月の社会の授業の最初に１時間かけて何のために歴史を学ぶのかということをひたすら子どもたちに問い続けた。みんなで考え続ける時間がとても長かった。で，結論は「温故知新（ふるきをたずねてあたらしきをしる）」ということだった。学校生活で数少ない印象的だった授業の一つである。

さて，学習漫画のほうに話を戻す。
わかったことは，日本の政治は悪者が担ってきたということだ。どおりで最近の政権支持率が40%で高止まりしているわけだ。子どものときからみんながこの37年に渡るベストセラーの学習漫画で育っている。そんなわけで，今の50歳以下の世代は我々と違って政治の本質を見抜いてしまっており，なおかつこれが歴史というものだとあきらめて（達観して）いるのかもしれない。その土壌の上に新自由主義が猖獗するのもしょうがない。

各国の感染者比（２）

各国の感染者比（１）からの続き

新型コロナウイルス感染症の話。
感染者数が一定の数あるいは人口比になった日時とその値を共通の原点として，各国の感染者数や人口比を対数グラフに表して観察する方法について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%とその中間になっている。

ワン・モア・ヌーク：藤井太洋（２）

3月14日のエピローグも含めて，3月10日で丁度ワン・モア・ヌークを読了した。たいへん知的でスリリングで映画にしたらおもしろいなあと思いながらスイスイと読み進むことができた。テーマとしては9年目を迎える東日本大震災に伴う福島第一原子力発電所事故にピタリと照準をあわせてきた。原子爆弾に関するディテールとこれにイスラム国のテロがうまく調和した内容になっている。もちろん本当の専門家から見ればなんだこれはということになるのだろうけれど，我々素人からすれば十分にリアルに感じられる。

残念ながら著者の福島原発被害者の不安に寄り添おうとする視点には共感できなかった。これはほとんど例のアンチニセ科学グループの思想と同じだったから。具体的には65ページの「テレビ，新聞，雑誌やインターネットを通して死の恐怖を感じさせる言葉をひっきりなしに浴びせられた中で・・・」から67ページの「実はこの程度の被爆でがんの発症例が増えたという疫学的な証拠はない」あたりまで。著者の説明はほとんど間違ってはいないかもしれないが，同時に警鐘をならす人々をすべて飲み込んでしまう効果をもたらしている。そして今も処理水という名のトリチウム＋未処理放射性物質の海洋投棄や小児甲状腺がんの診断をめぐってその闘いは続いている。キクマコよりおしどりマコが，早野龍五より牧野淳一郎が信頼できる。

ワン・モア・ヌークはおもしろかったのでこれは評価したい。ただ，こころに残るかというとどうだろうか。この程度の科学的な精度を持ちながらより深い思索をもたらすものを読んでみたい。安部公房は残念ながらもういない。

各国の感染者比（１）

各国の新規確定感染者数の累計（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時点，縦軸はppmの常用対数）

（１）各国の1ppm以下のデータはバラつきが大きく，感染症の数理シミュレーションの基礎データとして用いるには不適当である。今後，さらに各国の比較を行う場合は感染者比1ppm以上のデータで議論することが望ましい。

（２）感染者数が多い国おける1ppm以上の感染者比の増加カーブの様子はこれまでのところよく似ている。韓国では中国と同様の飽和が始まっているようだ。韓国，イタリア，イランは現在まで湖北省と中国全土に挟まれたデータ領域で推移している。

（３）感染源近傍の湖北省（武漢以外）に比べ中国全土の感染者比の飽和値は約1/6となっている。中国における湖北省の封じ込めが成功しているのだと考えられる。各国や地域の社会構造（政治・経済・文化）や感染対策効果に依存して感染者比の飽和値が定まっていく。

（４）この中では日本のデータだけが他の国や地域と異なった振る舞いをしている。しかしまだ飽和に達している様子は見られない。なお，1ppmに達していない国のデータを混ぜて議論するべきではない。人口密度の低いオーストラリアなども日本と同様の傾向ではないかとの指摘があり，引き続き原因の分析が必要である。


（注１）日本のデータはWHOのデータの扱いと同様で，ダイヤモンド・プリンセス号の部分は除かれている。

（注２）湖北省（武漢以外）はWHOのデータからは直接得られないので，湖北省（武漢以外）のデータによる考察とモデル計算に基づいている。

（注３）WHOの日本のデータ2月5日分には奥村晴彦さんのCOVID-19が指摘するような誤記があるので修正した。なお，中国のデータも奥村さんのcsvファイルを利用させていただいた。

（注４）中国全土のデータの2月17日の段差は集計方法の変更によるものである。

（注５）＊1/30 -> 1/21 に訂正（グラフにはまだ反映されていない 3/15/2020）

［１］感染症の数理シミュレーション（６）

各国の感染者比（２）に続く

fudan-CCDCモデル

中国の上海にある復旦大学の研究グループが，新型コロナウイルス感染症（COVID-19）の流行を推測するために，統計的な時間遅延力学系モデル（fudan-CCDCモデル）を考案した。CCDC（中国疾病予防センター）は中国の防疫拠点であるがそこで得られたデータを再現している。彼らはこのモデルを用いて [1] の論文を2月21日に報告した。

日本でも，このモデルによる感染予測が始まった [2][3][4]。いずれも厚生労働省のデータをもとにしてパラメタを推定している。KEN_KEN2さんのpythonを使った計算によれば，3/10ごろに確定感染者が2000人となりそうだ。3月のはじめの数日の新規感染者数については厚生労働省の報告値の2倍以上の推計値を与えており，モデルが原因なのかデータが原因なのかはわからない。また，林祐輔さんのpythonを使った計算によれば，政府の観測にかからない潜在感染者数が 3月9日時点で約2万人あり，ピークアウトするのは4月末〜5月初旬という結果を得ている（注：この潜在感染者数は我々のSIIDR2モデルでは$u_2(t)$に対応する）。いずれもソースコードが公開された貴重な資料となっている。

ところで，復旦大学グループの論文[1]の結論部分には次のようなコメント「上海と同じ防疫措置を2/20までに実行すれば潜在感染者数は15万人に押さえられるが，2/29にのびれば45万人になるだろう」がある。彼らは日本の感染者数の初期データとしてダイヤモンド・プリンセス号を含めたものを採用し，武漢とほぼ同じ値で増加していることを冒頭で示している。また，これを用いて複数のシナリオに基づくシミュレーションを行っている。残念ながら，これにはダイヤモンド・プリンセス号という特異事例の感染者数の値，すなわち2月10日に135人（日本国内のみ26人）から2月19日に542人（日本国内のみ73人）を含んでおり，第5節のPossible future scenarios for COVID-19 in Japanのシミュレーション部分はミスリーディングである。

［１］CoVID-19 in Japan: What could happen in the future（Fudan University）
［２］中国の最新論文の方式で日本のコロナウイルスの感染者を予想してみた（KEN_KEN2）
［３］ Fudan_CCDC_Japan_Optuna.ipynb（林祐輔）
［４］COVID-19 in Japan（林祐輔）
［４］Coronavirus disease (COVID-2019) situation reports（WHO）

感染症の数理シミュレーション（６）

感染症の数理ミュレーション（５）からの続き

湖北省（武漢以外）のデータでは，新規/累計感染者数 f(t), g(t) と新規/累計死亡数 h(t), i(t) を簡単な近似関数で表現した。これを用いると，感染症の数理シミュレーション（３）のSIIDR2モデルにおけるパラメータを絞り込むことができるのではないだろうか。

そこで，前回示したSIIDR2の方程式系と今回の実測データ近似関数の関係を整理してみよう。
\begin{equation}
\begin{aligned}
S: \dfrac{du_1}{dt} &= -\dfrac{\beta}{n} u_1 u_2 = -\dfrac{\beta}{n}\alpha_2 f u_1 \\
I_1: \dfrac{du_2}{dt} &= \dfrac{\beta}{n} u_1 u_2 -\dfrac{u_2}{\alpha_1} -\dfrac{u_2}{\alpha_2} = \dfrac{\beta}{n}\alpha_2 f u_1 - (1 + \dfrac{\alpha_2}{\alpha_1}) f \\
I_2: \dfrac{du_3}{dt} &= \dfrac{u_2}{\alpha_2} - \dfrac{u_3}{\gamma_1} -\dfrac{u_3}{\gamma_2} = f - (1 + \dfrac{\gamma_2}{\gamma_1}) h \\
D: \dfrac{du_4}{dt} &= \dfrac{u_3}{\gamma_2} = h\\
R: \dfrac{du_5}{dt} &= \dfrac{u_2}{\alpha_1} + \dfrac{u_3}{\gamma_1} = \dfrac{\alpha_2}{\alpha_1} f + \dfrac{\gamma_2}{\gamma_1} h\\
I_a: \dfrac{du_6}{dt} &= \dfrac{u_2}{\alpha_2} = f
\end{aligned}
\end{equation}
これから次の関係式が得られる。
\begin{equation}
\begin{aligned}
u_1 &= n - u_2 - u_3 - u_4 - u_5\\
u_2 &= \alpha_2 f\\
u_3 &= \gamma_2 h\\
u_4 &= \int h dt = i\\
u_5 &= \dfrac{\alpha_2}{\alpha_1} g + \dfrac{\gamma_2}{\gamma_1}i \\
u_6 &= \int f dt = g
\end{aligned}
\end{equation}
さらに，次の前提条件が成り立つとする。
①感染者の2割が発症する[1]。
②致命率（発症者の死亡率）は0.7%である[2]。
③発症までの平均潜伏期間は$\alpha=5$とする$(1/\alpha=1/\alpha_1+1/\alpha_2)$ [3]。
④基本再生産数$R_o$=2.2(1.4〜3.9)より，$\beta=R_o/\lambda=0.46\pm0.12$  である [2]。
⑤湖北省（武漢以外）の新規死亡数h(t)は20日のピークに0.006/日，40日目の累計で0.12人（１万人当り）[4]。
⑥同上の新規感染者数f(t)は12日のピークに0.25/日，40日目の累計で3.70人（１万人当り）[4]。

$I_2$の$u_3$に対する微分方程式$du_3/dt + u_3/\gamma = f$を解いて求めた$u_3$の関数形が$\gamma_2 h$に近づくように定めると，$\gamma=5,\   \gamma_1=5/0.965,  \gamma_2=5/0.035$が得られる。ただし，$(1/\gamma=1/\gamma_1+1/\gamma_2)$である（訂正有り 3/13/2020）→感染症の数理シミュレーション（７）へ

各関数のT日目の値がわかれば，非感染者の減少分$\Delta u_1$は次式で与えられる。
\begin{equation}
\Delta u_1=n-u_1(T)=u_2(t)+u_3(T)+u_4(T)+u_5(T)
\end{equation}
さらに湖北省（武漢以外）のT=40日のデータを代入して，$u_2(T)=\alpha_2 f(T)≒0$, $u_3(T)=\gamma_2 h(T)≒0$より，$\Delta u_1= \dfrac{\alpha_2}{\alpha_1} g(T) + \dfrac{\gamma_2}{\gamma_1}i(T)+i(T) =14.8+3.4=18.2 $である。ここで，$\alpha_1=5/0.8, \alpha_2=5/0.2$という①の仮定を用いた。

$u_1(t)$において，$\beta$が定数であると仮定して微分方程式を解けば，$u_1(t)=n \exp (-\frac{\beta}{n} \alpha_2 g(t))$となることから，$\beta$と$\alpha_2$を与えれば，$u_1(T)=9982$が得られる。先ほどの式と組み合わせると$\beta=0.2$となる。

まとめると，湖北省（武漢以外）のデータを再現するSIIDR2モデルのパラメタのセットは次のようになる。（訂正有り 3/13/2020）→感染症の数理シミュレーション（７）へ
\begin{equation}
\alpha_1=\dfrac{5}{0.8},\  \alpha_2=\dfrac{5}{0.2},\  \beta=0.2,\  \gamma_1=\dfrac{5}{0.965},\   \gamma_2=\dfrac{5}{0.035}
\end{equation}

SIIDR2モデルでは，感染率を時間の関数$\beta(t)$として感染対策が行われることを表している。上記の$\beta$はそれを平均したものと考えることができる。ところで，上記の$I_2$の$u_2$に対する微分方程式は，$du_2/dt + u_2/\alpha = \beta u_2$と近似できる（$\therefore u_1/n ≒ 1$）ことから，次のように与えられる。
\begin{equation}
\beta(t) = \dfrac{du_2/dt}{u_2} + \dfrac{1}{\alpha}
\end{equation}
$u_2(t)=\alpha_2 f(t)$を上式に代入したものと，感染症の数理シミュレーション（５）で考えたモデル感染対策時間因子で$\beta=0.91, \lambda=7, \tau=7$と置いたものを下記の図に示す。

図：湖北省（武漢以外）のデータからSIIDR2で導いたβ(t)(下)とモデル化したβ(t)(上)

また，このモデル化された感染対策時間因子$\beta(t)$をt=0からt=40まで積分したものを40日で割って平均すると（$T=40, \beta=0.91, \lambda=7, \tau=7$），
\begin{equation}
\bar{\beta} = \dfrac{1}{T} \int_0^T \dfrac{\beta}{15} (8 + 7 \tanh (-(t-\tau)/\lambda ) dt = 0.24
\end{equation}
となって，先ほど湖北省（武漢以外）のデータからSIIDR2モデルを経由して導いた$\beta=0.2$と矛盾しないような感染対策時間因子のパラメタセットが存在することがわかる。

［１］新型コロナウイルス感染症対策の見解（厚生労働省）
［２］新型コロナウイルス(COVID-19)感染症（川名・三笠・泉川）
［３］新型コロナウイルス肺炎COVID-19（下畑享良）
［４］tenetニュースサイト
［５］Coronavirus disease (COVID-2019) situation reports（WHO）
［６］COVID-19（奥村晴彦）

感染症の数理シミュレーション（７）へ続く

感染症の数理シミュレーション（５）

感染症の数理シミュレーション（４）からの続き

ここ１，２週間が瀬戸際といわれ続けて，２週間を迎えようとしている。そこで，対策の効果や時期を感染率$\beta(t)$という１つの時間依存のパラメタに代表させ，これが感染者の増加にどのような影響を及ぼすかを見えるようにしたい。

感染症の数理シミュレーション（３）では，SIIDR2モデルとして，感染対策によって感染率を下げる効果を，感染シミュレーション開始時を基点とする指数関数的な因子$\beta(t) = \beta \exp(-t/\lambda) $によって導入した。これは，感染対策時間$\lambda$の数倍で感染率をほとんどゼロに近づけてしまうが，それは難しいのではないだろうか。

そこで，感染率を1桁程度減少させるとともに，対策導入日$\tau$を新たなパラメータとして加えた新しい対策時間因子を含む感染率を定義する。
\begin{equation}
\beta(t) =\dfrac{\beta}{15} \{ 8 + 7 \tanh \dfrac{-(t - τ)}{\lambda} \}
\end{equation}
感染対策の前後にはそれぞれ$\beta$及び$\beta/15$という一定の感染率を与える（15にはとくに根拠はない，パラメータの定量的な精度は求めておらず，オーダやファクターの評価ができればよいという立場である）。もし，対策開始日が$\tau=0$であれば，シミュレーション開始日の$t=0$の$\beta(0)=8\beta/15$となるので，これまでの$\beta$の定義の半分の値になっていることに注意する。

$\beta(t)$を$\beta$で割って，感染対策時間因子だけを取り出したものをグラフにすると次のようになる。横軸はシミュレーション開始からの日数を表す。

図１　感染対策時間因子β(t) λ=5 左からτ=5,10,15

図２　感染対策時間因子β(t)  τ=20 急峻な方からλ=5,10,15

どうやら政府は情報公開に力を入れるのではなく，マスコミやインターネットにおける情報抑制に全力を投球し始めた。ますます，正確な情報が見えなくなってしまう。

感染症の数理シミュレーション（６）に続く

ワン・モア・ヌーク：藤井太洋（１）

2020年3月6日を迎えたので，藤井太洋のワン・モア・ヌークを小説内の時間経過と同期しながら読み始めることにする。藤井太洋はオービタルクラウドとGeneMapperを読んだはずなのだが，どんな話だったのかまったく印象にない。本棚を探してみると，前者はみつからない。退職時の大掃除のときに捨てたのかもしれない。それでも藤井太洋には注目しているのは，そのITへの造形の深さと中国SFブームへの対応姿勢からだろうか。

湖北省（武漢以外）のデータ

新型コロナウイルス感染症の事例数が最も多い地位は中国の湖北省であるが，武漢のデータは特異的な値を示している（高い死亡率）。そこで，モデルの妥当性を検討するためには湖北省（武漢以外）のデータを用いるのが適当だ。tencentのニュースサイトには湖北省（武漢以外）のデータが報告されている。

湖北省（武漢以外）の新規感染者数データの実測値 K(t) にはバラつきがあるが，次のような傾向を持っている。
・K(t)は40日の範囲で0から増加してその後減衰して0に近づく。
・そのピークはt=14日のK(t)=1200人である。
・この40日間の感染者数の累計数は$\int K(t) dt$ = 17,800人である。

この性質をみたす簡単な近似関数を求めておけば，感染モデルとの比較が容易になる。K(t)の解析的な近似関数をf(t)（一日に報告される新規感染者の数）とし，その積分をg(t)（新規感染者の累計数）とする。

目の子で探すと，tを報告開始からの日数として次のようになる。
$f(t) = 3 t^4 \exp(-t/3)$
f'(12)=0, f(12)=1140, g(40)=17400 など，観測データをおおむね再現している。

感染症の数理シミュレーション（３）で用いたSIIDR2モデルでは，$u_6$がg(t)に対応する。なお，$u_6(t)$は重症感染者数$u_3(t)$の変化分から減少項を除いたものを積分したものである。

図１　湖北省（武漢以外）の新規/累計感染者数の近似関数，上がg(t)，下がf(t)（横軸は日，縦軸は１万人当たり）

湖北省（武漢以外）での1万人当りの累計感染者数（t=40）は3.7人，新規感染者数（t=peak）は0.25人，死亡数（t=40）は0.12人である。

P. S.
なお，死亡数についても20日をピークとする滑らかな関数となっているので，新規/
累計死亡数の近似関数を求めることができる。それぞれh(t), i(t)とすると，
$h(t) = 30 \sin^2[ \pi / 40 t]$で与えられる。ピーク時の一日あたりの死亡数は30人であり，40日目にはおおむね収束に近づいている。これをグラフにすると次のようになる。

図２　湖北省（武漢以外）の新規/累計死亡数の近似関数，上がi(t)，下がh(t)（横軸は日，縦軸は１万人当たり）

あぶり餅

京都の紫野今宮神社には名物のあぶり餅屋がある。境内前に２軒ありそのうちの一文字屋和輔，一和のほうに入った。新型コロナウイルス感染症の影響で，京都は観光客が少なく，あぶり餅屋も手持ちぶさた。庭のみえるお座敷に通されたが，我々夫婦の他には誰もいない。名物のあぶり餅はおいしかった。その縁起によれば，一条天皇（980-1011）の頃に悪疫がはやり，今宮神社で悪疫退散の祈願をした。そのときに初代が，この阿ぶり餅を神前に供え，そのお下がりを喰う人が疫病をのがれたとのこと。一和の創業は長保二年（1000）とある，すごすぎる。というわけで，疫病退散の今宮神社と一和を訪れたのはまことにタイムリーだったわけである。

疫病対策に大仏を建立したり，祇園祭を始めたり，あぶり餅を普及させたりと，奈良や京都では昔から工夫をこらし，その成果が観光資源として現在まで続いている。疫病の社会に対するインパクトはこれほどまでに大きい。このたびの新型コロナ疫でも後世に残るような文化資産につながる対策を考案する必要があるのではないかいな。

　　　

写真：あぶり餅の一和（2020年3月3日撮影）

休校中の高校生に贈る「情報」の課題Ａ

［課題Ａ］
新型コロナウイルス感染症の拡大を防ぐためという理由で，3月2日から全国ほぼ一斉に小中高等学校が休校となりました。これが本当に効果があるのかどうかを検討し，より良い対策がないかをグループで相談しながら考えてください。なお，国会中継がある日は必ずそれを見るようにしてください。

①日本における人々全体の平均的な１日における感染機会Ｃ(活動)は，(1)その活動をする人数の割合（分母は日本の人口），(2)その活動時間（分母は24時間），(3)その活動中の平均的な人の密集度（単位は人/㎡），の３つの積をすべての人について加えたものに比例すると仮定します。このとき，全国の学校を休校にした効果，すなわちＣ(平常時)-C(休校時)はいくらになるでしょうか，考えて下さい。
（ヒント）ある人のブログに「全国の小中高等学校の休業の効果」という計算例があるので，どうしてもわからない場合にはこれも参考にして下さい。

②全国の学校を休校にする以外で，社会を大きく混乱させずに感染機会を減らすにはどうすればよいか，それぞれについてC(活動)を計算して，最も効果的な対策を考えなさい。
（ヒント）通勤ラッシュの解消，カラオケ，居酒屋，パチンコ，スポーツジムにいかない，等。ただし，必ず数値的な根拠を引用データによって示してください。

③上記のシミュレーションの前提条件は，どの程度正しいか，どこが間違っている可能性があるのかについて友達と議論してください。
（ヒント）友達がいない場合は家族の方々でも結構です。

全国の小中高等学校の休業の効果

新型コロナウイルス感染症の感染が拡大する中，2月27日(木)の夕方，安倍首相は全国すべての小学校・中学校・高校・特別支援学校等に，3月2日から春休みまでの期間，臨時休業とするよう要請する考えを表明した。そして今日からその休業が始まった。

ここではその是非ではなくて，感染機会を減少させる上で休業の効果がどの程度ありうるかについて考えてみたい。人間の１日の生活時間において感染の機会は，その活動における他人との接触による（間接的なモノを通しての感染もあるとは思う）。このとき，人との距離が近い密集した空間のほうがその危険性は大きいとNHKでも宣伝していた。そこで，１日の活動時間とその活動時の単位面積当りの他者の人数の積和を全国民に渡って平均すれば，感染機会の大きさを定量的に評価できるのではないかと考えた。

まず，生活時間を調べよう。５年おきに実施している最新の平成28年社会生活基本調査のデータをみよう・・・うーん。面倒なので勝手に考えることにする。睡眠7時間，家庭・食事5時間，通勤・通学等1時間，仕事・学校・塾8時間，その他活動3時間。

また，次に平成27年国勢調査と学校基本調査を参考にしよう・・・うーん。まあだいたい，1億2700万人の国民が，幼稚園・幼保連携型認定こども園・保育園に400万人(3%)，学校（小〜大）に1700万人(13%)，職場に5900万人(47%)，老人・乳幼児等が4500万人(35%)，病院・介護施設等に200万人(2%)となる。国勢調査によれば世帯数は5300万世帯で，世帯人数が1人(34%)，2人(28%)，3人(18%)，4人(13%)，5人以上(7%)となっている。

まず，家庭での時間を考える。15㎡の部屋に家族が集まって食事や団欒の時間を過ごす。その時間が3時間として，睡眠時間を含め9時間は1人でいるものとする。そこで感染機会Ｃに寄与するのは複数で過ごす時間である。
Ｃ(家庭) = 0人/㎡*9時間 + (0*0.34 + 1*0.28 + 2*0.18 + 3*0.13 + 4*0.07)人/15㎡*3時間 = 0.21

次に，学校と職場である。教室の面積を9*7=60㎡として，平均収容人数を30人とすると，
C(幼保+学校) = (0.03 + 0.13)*30人/60㎡*8時間 = 0.64
職場は非常に多様な環境であり，えいやっと１人１坪のオフィスや工場を考える。
C(職場) = 0.47*1人/3.3㎡*8時間 = 1.14

通勤通学時間は，NHKによる2015年の国民生活時間調査と国勢調査が参考になる。徒歩(7%)，自転車等(15%)，自家用車(47%)，電車バス等(31%)である。また，電車の１人当り立席専有面積は0.3㎡らしいので，3.3人/㎡を用いることにする。
C(通勤+通学) =(0.03+0.13+0.47){ (0.07+0.15+0.47)*0人/㎡*1時間 + 0.31*3.3人/㎡*1時間 }= 0.64

その他はよくわからないので，ざっくりと職場と学校の平均0.4人/㎡をとる。
C(その他) = 0.98* 0.4人/㎡ *3時間 = 1.18

病院・介護施設もよくわからないが2%なので家庭に含めてしまった。したがって，
C(全体) = C(家庭) + C(幼保+学校) + C(職場) + C(通勤通学) + C(その他) = 0.21+0.64+1.14+0.64+1.18 = 3.81

学校が休業になるとC(学校)とC(通学)の分は減るが，半分は学童保育などで残るだろう。
なお，小中高等学校だけなので，全人口に対する寄与率は0.13ではなく，高等教育機関などを除いた0.10である。そこで学校一斉休業による減少項は，
C(休校) = {0.10*30/60*8 + 0.10*0.31*3.3 } * 0.5 = 0.25

つまり，この度の全国一斉学校休業による感染機会Cの減少効果は，0.25/3.81≒7%ということになる。まあこんなものだろうか。これからみると，さらに努力することができる部分は家庭や仕事以外のその他の部分である。

P. S. 時差出勤やテレワークなどによって，通勤時の混雑緩和をはかり，乗車率が50%になれば，C(通勤)= {0.47*0.31*3.3 } = 0.24 だけ減少するので，ほぼ全校休校措置と同程度の効果をもたらすことになる。かもしれない。

［１］新型コロナウイルス感染症対策本部（首相官邸）
［２］新型ウイルス研究者試算 対策組み合わせで入院患者６割以上減（NHK）

感染症の数理シミュレーション（４）

感染症の数理シミュレーション（３）からの続き

新型コロナウイルス感染症に対するWHOなどの共同調査報告書［１］が公表された。同時に，新型コロナウイルス感染症の世界的なリスクは４段階のHighからVery Highに引き上げられた。

報告書によれば，感染すると平均で5〜6日後に症状が出る。感染者のおよそ80％は症状が比較的軽く，肺炎の症状がみられない場合もあり，呼吸困難などを伴う重症患者は全体の13.8％，呼吸器の不全や敗血症，多臓器不全など命に関わる重篤な症状の患者は6.1％だとされた。

逆に子どもの感染例は少なく，症状も比較的軽い。19歳以下の感染者は全体の2.4％にとどまり，重症化する人はごくわずかだ。子どもの感染について報告書では主に家庭内で大人から感染しており，子どもから大人への感染例は調査をした中では確認されていないとしている。

ということは，安倍が週末にある世論調査の支持率回復を念頭に慌てて実施した小中高等学校の全国一斉休校要請は十分な根拠がないということだろう。危機的な状況を避けるための適切な判断だという説も多いが，適切な分析やシミュレーションには裏付けられていない（危機感の喚起には一定の意味があるかもしれないが，他がじゃじゃ漏れな状態で一部だけに注力するとよけいな混乱と感染拡大を招きかねない）。

2月20日までのデータによれば，5万5924人の感染者のうち死亡したのは2114人で，全体の致死率は3.8％である。ただし，感染拡大が最も深刻な湖北省武漢だけが致死率が5.8％であり，1月の最初の10日間に発病した患者の致死率は17.3%に達しているので，ここは特異的状況として除いて考えるべきだ。その他の地域の最近の致死率は0.7％である。疫学における致死率（致命率＝case fatality rate）とは，特定の疾患に感染した母集団に対する死亡数の割合である。ここで，疾患に感染していることが判明している必要があるので，今回のようにほぼ無症状や他の疾患と区別できない軽症の感染者がいる場合には取り扱いに注意がいる。

感染症の数理シミュレーション（３）では，東大の黒田産業医のデータをもとにパラメタを定めたが，今回のWHO報告とほぼ整合している。ウイルスの伝搬可能な感染者全体のうち症状があって報告に計上される感染者が1/3である（残りの2/3はほぼ無症状や他の疾患と区別できない軽症の感染者）と仮定すれば，上記の13.8%の1/3の4.6%が重症な感染者であり，モデルの$\alpha_2 / \lambda$=5%とほぼ一致する。また，致死率0.7%もウイルスの伝搬可能な感染者全体を分母にすると1/3になるので0.23%であり，$\gamma_2 / \lambda $=5%としたときの，$\alpha_2 /\lambda  \cdot \gamma_2 / \lambda $ = 0.25%とほぼ一致する。

前回のSIIDR2のシミュレーションでは，武漢を除く湖北省のデータに近くて日本に妥当するモデルとして，$\beta=0.20, \lambda=100$を例示した。これは，感染者の初期値が1万人に1人という仮定を含んでいた。感染者が急増した武漢を取り巻く湖北省全体では，それぞれの地域に飛び火した数千人の初期感染者から計算を開始することは，一定の合理性があるが，日本に換算すると，全国で1万人の初期感染者を仮定することになって不自然である。

つまり，SIIDR2を用いた日本の現状のモデルシミュレーションには，1万人当りの感染者の初期値$u_2(0)$，感染率の初期値$\beta$，対策による感染率の減少係数$\lambda$の３要素が必要である。いずれにせよ，合理的な判断ができなくて肝腎の検査がシステマティックに行われず，感染状況に関するデータが完全に欠落している日本では推定が難しい。

［１］Report of the WHO-China Joint Mission on Coronavirus Disease 2019 (COVID-19)

昭和のプログラミング教育（４）

昭和のプログラミング教育（３）からの続き

金沢泉丘高校理数科の２年になって，カシオのAL-1000の次に取り組んだのがアナログ計算機である。

野田中学校３年のときの技術の先生が黒丸の眼鏡をかけた辻先生といって，やんちゃな生徒まで含めて全員が恐れるたいへん怖い先生だった。電気や電子回路が専門のようで，技術準備室には分解・組立中のテレビなどが山積していた。電気の授業では，電気工事士の実技試験のような設定で，木の板に碍子を固定して電線を貼るというものがあって，クラスから２名選ばれた。ドライバーなどの工具を与えられ持ち時間10分くらいで取り組んだ。もう一人は誰だったかはっきりしないが，なんと自分も選ばれてしまったのだ。練習も何もなしで「始め」といわれてもウロウロするだけで碍子一つ取り付けれずに終了した。チーン。もう一人は途中までできていたが，自分はあまりにも不器用だった。

３年の技術の授業の後半にはいよいよ真空管式ラジオの回路図が登場した。五球スーパーか何かのキットを買ってもらって真空管式ラジオを実際に組み立てたのは高校生になってからだったが，それまでも，はやくから電子工作のおもちゃは家にあり，ハンダごてを使ってトランシーバーのキットに挑戦したりと，小学生の時の鉱石ラジオから始まった電子回路への興味はとても大きかった。授業にも真剣に取り組み（ラジオの実習はなかったと思う），試験前には教科書の回路図を完全に記憶していた。試験では微分回路と積分回路が出て，そこだけが間違っていた。採点の方が誤りだったので，恐る恐る先生に申し出たところ，無事に訂正してもらえてよくできたと褒められた。

話が長くなったが，その積分回路と微分回路を使うと，高校数学で出てくるような微分方程式（当時は高等学校の数学でも簡単な微分方程式を扱っていた）を解くことができる。この原理を利用したのがアナログ計算機であり高専などでも使われていた。20年ほど前に，インターネットで探してみると神戸高専のホームページに丁度自分たちが高校時代に使っていたものの写真が見つかった。それからしばらくするとそのページも無くなってしまい，今では完全に過去の遺物となってしまった。

図１　理数科で使ったアナログ計算機（引用：神戸高専の古いウェブサイトより）


北國新聞の記事にある楠先生の収集資料の中には，昭和48年9月に全国理数科高等学校長会によって発行された理数科紀要の第１集があり，その一部のデータを送っていただいた。当時，非常にがんばってこのアナログ計算機の実験に取り組んでレポートを書いたので，楠先生に取り上げられて，あちこちの発表に利用されたようだ。同級生の北原さんと連名になっているが，ほとんど私の作品だ。昔はよく勉強していたことがわかる。

図２　理数科紀要の中で報告された当時の私のレポート

アナログ計算機は理科の準備室かどこかに設置されていて，放課後一人で残って黙々と楽しい作業を続けた。微分方程式を解いた計算結果は振動や過渡現象のグラフになって出力される。あるとき，配線を間違えてスイッチをいれると警告ブザーが鳴り響き，学校に１台しかない高価な機械を壊してしまったかと肝を冷やした。アナログ計算機は無事だった。なお，図２のレポートをワープロで書いたのは楠先生だ。自分が書いたレポートは手書きであり，ちゃんと返却してもらったのだが，50年のうちに紛失してしまった。

感染症の数理シミュレーション（３）

感染症の数理シミュレーション（２）からの続き

一昨日のニュース。2020年2月26日，基本方針が全く具体的でないと厳しく指摘されてあせる政府の要請により，PerfumeやEXILEのコンサートが中止になった。また，3月11日まで国立劇場の公演は延期や中止，国立博物館・美術館が閉館となり，社会の空気が相転移した。一方，大阪でただ1名の感染報告者のガイドさんが，回復・退院後再び症状が現れてしまった。さらに，昨日は安倍官邸が全国の小中学校の休校を要請するに至った。SNSでは根拠が不確かな情報がひろがっているということなので，このblogもあまり信用しないほうがいい（ソースコードは提供しており反証は可能だ）。

さて，パラメタが４つあれば象の絵が描けるといわれているので，５つもパラメタを含んだモデルで遊ぶのは慎んだほうがよいのだ。しかし，精緻な理論化に取り組む前に，現象の本質につながるトイモデルであれこれ試すのは，物理学あるいはコンピュータ・シンキング（プログラミング的思考じゃないダス）のキモである。前回のモデルにおけるパラメタセットの問題点を指摘し，それを除く修正モデルを提案して今後の感染状況の推移の全体像の可能性を探りたい。

前回のSIIDRモデルには６つのパラメタが現れる。このうち集団の人口$n$は固定されており，これは問題ない。感染率（１人１日あたりの感染人数）$\beta$や軽症感染者からの平均離脱日$\alpha_1, \alpha_2$，重症感染者からの平均離脱日$\gamma_1, \gamma_2$の５つの範囲をどのように合理的に推定するかがポイントだ。

東京大学の黒田玲子産業医が個人的な見解としてまとめたメモによれば，新型コロナウイルス感染症の潜伏期は1〜12.5日（多くは5〜6日）で，感染力は1人→2〜3人，予後はほぼ治癒だが，ICU入院が5%弱で致死率が0.7%とある。したがって，軽症感染者の平均離脱日を7日とするのはまあまあだ。また，重症感染者への移行率は5%とした。黒田メモの分母は感染者だが，我々のモデルの感染者には発症していないものも含む（5%より小さくなる）。その一方でICUまで至らない重症感染者も考えられる（5%より大きくなる）ので，とりあえず5%を採用する。そこで，$\alpha_1=7/0.95, \alpha_2=7/0.05$とする。また，重症感染者についても同じ7日を仮定し，重症感染者からの死亡率を5%とした。つまり，$\gamma_1=7/0.95, \gamma_2=7/0.05$である。その結果，このモデルの感染者致死率は0.25%となるが，ここでの軽症感染者には発症していないものを含めているので，上記の致死率0.7%とは矛盾しない。こうして４つのパラメタは合理的な範囲で収められる。

最も重要な不確定因子は感染率$\beta$である。これは，ウイルスとホストの生物学的な特性でほぼ決まるが，同時に，ホストである人間集団の社会構造（人口密度，経済活動，生活習慣，文化的要因）にも強く影響されるので，想定する集団によって大きく変わりうる。一方で，感染力（基本再生産数 Ro）は2.2と報告されているため，軽症感染者の平均離脱日数7日で割れば，$\beta=0.3$あたりがこの値の想定範囲となる。

とりあえず，前回のモデルをダイヤモンドプリンセスに当てはめて計算してみた。初期値として5%が感染したところから初めた（2/28/2020日経朝刊の対談より）。$\beta=0.3$のときに3700人あたり4人の死亡数（30日目）が再現できる。このときの累積重症感染者は3700人に換算して590人であり，報告されている700人とほぼ対応している。

次に，このパラメタセットで中国の湖北省の状況を説明できるかどうかを確認したい。モデルの軽症感染者には感染力を持つ非発症者を含んでいるため，報告されている感染者数とは直接比較できない。一方，軽症感染者数の時間構造と重症感染者数の時間構造は１週間ほどのタイムラグがあるがほぼ同形である。また，死亡者数についてもあいまいさが少ない。そこで，我々は，(1)モデルの重症感染者数（1万人当り）の感染者数飽和時期と(2) 現時点の死亡数（1万人当り）の２要素を湖北省のデータと比較する。2019年12月8日に最初の報告があってから，ほぼ３ヶ月たった現時点で湖北省の感染者数は飽和しそうである。また，湖北省の死亡数は約2700人だ（湖北省総人口は5900万人で人口密度は日本と同程度。また，死亡数の8割が集中する武漢は人口が1100万人で人口密度は愛知県なみ）。湖北省の感染者数飽和時の1万人当りの死亡数は0.5人のオーダーと推定される（特異的な武漢を除く湖北省では0.1人）。

図１　β=0.3(左の山/丘)とβ=0.2(右の山/丘)の重症感染者数(u4)と死亡数(u4)

図１のように，上記のパラメタで3ヶ月後(90日目)に感染者数飽和を与える$\beta$は0.30前後だが，このときの死亡数は20人となり1桁以上違う。一方，$\beta=0.20$として，3ヶ月後の死亡数を0.5人あたりに近づけると感染者数飽和時期（右の山の累積値の飽和）が7ヶ月後になり，このモデルでは湖北省データをうまく説明できない。その原因は，中国の強力な対策によって感染率$\beta$を制御できたことにあるのではないかと考え，$\beta$が時間$t$の関数として減少する図２のSIIDR2モデルを導入する。パラメタ数をなるべく減らしたいので，ここでは，$\beta(t)=\beta \exp(-t / \lambda)$という１パラメタの指数関数を導入する。$\lambda$は感染対策時間因子であり１〜数ヶ月のオーダーだ。短期間の変化は１次関数で近似される。

図２　感染率がβ(t)として時間に依存するSIIDR2モデル

図３のように，$\beta=0.30, \lambda=65$で，湖北省データと矛盾しない結果が得られる。これは感染初期に3ヶ月で感染率を1/4まで減らしたことに対応し，実際中国は武漢封鎖などの強力な措置によってこれを達成したのかもしれない。このときの飽和感染者数は70人（1万人当り）であり，武漢に換算すると8万人となる（非発症者を含む）。なお，この感染対策時間因子が$\lambda=100$にとどまれば，重症感染者数の飽和数は数倍に跳ね上がる。

図３　小さな山/丘(λ=65)が湖北省データに対応，大きな山/丘はλ=100の場合

この湖北省で得られたパラメタをこのまま日本に当てはめると，30日目の現在の重症感染者数は0.65人/1万人，0.05人/1万人である。つまり現時点における全国の重症感染者数はu3=6500人，死亡数はu4=500人に達することになる。しかし，たいへん残念なことに日本には正確なデータがない。通常の肺炎の死亡者は，年間10万人（日平均300人）なので，埋もれて観測されていない可能性も排除できないがそれはないと思いたい。例えば，$\beta=0.20, \lambda=100$とすると（武漢を除く湖北省に近いケース），30日目の現時点においてそれぞれu3=1000人とu4=100人となり（観測されていないものがあるという考え），よりマイルドな結果を与える。この場合，3〜4ヶ月後に飽和する重症感染者数は全国で7〜8万人，最終死亡数は500人のオーダーとなる（パラメタ次第で，結果を象の鼻を動かすように変えられるのがいたい・・・orz）。

P. S. 下記の参考文献[6][7][8]のように，中国からはよりまともな遅延ダイナミカル・シミュレーションの結果が報告されている。日本では誰か計算していないのだろうか・・・

［１］緊急寄稿新型コロナウィルス(COVID-19)感染症（川名明彦・三笠桂一・泉川公一）
［２］2019-nCoVについてのメモとリンク（中澤港）
［３］人口と感染症の数理（稲葉寿）
［４］数理モデルによる感染症流行制御の検証（大森亮介）
［５］新興感染症の国際的伝播を予測する数理モデル（西浦博・木下諒）
［６］CoVID-19 in Japan: What could happen in the future?（medRxivから）
［７］The reproductive number R0 of COVID-19 based on estimate of a statistical time delay dynamical system（medRxivから）
［８］A Time Delay Dynamical Model for Outbreak of 2019-nCoV and the Parameter Identification（arxivから）

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
using DifferentialEquations
using ParameterizedFunctions
using Plots; gr()

sky = @ode_def SIIDR2_model begin
  du0 = 1 # u0:time
  du1 = -β*exp(-u0/λ)*u1*u2/n # u1:Noimmunity(Susceptible)
  du2 = β*exp(-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 = u3 # u6:Accumulated Infection
end n α1 α2 β γ1 γ2 λ

n=10000.0 #total number of population
α1=7.0/0.95 #latent to recovery (days/%)
α2=7.0/0.05 #latent to onset (days/%)
β=0.3 #infection rate
γ1=7.0/0.95 #onset to recovery (days/%)
γ2=7.0/0.05 #onset to death (days/%)
λ=65 #pandemic supression time (days)

function epidm(n,α1,α2,β,γ1,γ2,λ)
u0 = [0.0,9999.0,1.0,0.0,0.0,0.0,0.0] #initial values
p = (n,α1,α2,β,γ1,γ2,λ) #parameters
tspan = (0.0,30.0) #time span in days
prob = ODEProblem(sky,u0,tspan,p)
sol = solve(prob)
return sol
end

β=0.20
λ=100
@time so=epidm(n,α1,α2,β,γ1,γ2,λ)
plot(so,vars=(0,4))
plot!(so,vars=(0,5))
#plot!(so,vars=(0,7))
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

感染症の数理シミュレーション（４）に続く

昭和のプログラミング教育（３）

1969年，金沢泉丘高校の我々理数科二期生のプログラミング教育は，CASIO AL-1000（1967年発売）からはじまった。Nigel ToutによるVintage Calculators Web MuseumのDesk Electric Calculatorsのトップにきているのが，このCASIO AL-1000 である。楠先生によれば，泉丘高校に４台が導入されたとのことだ。当時の価格で30万円だときかされていた。今の物価に換算するとその数倍に相当するのではないか。

AL-1000の入力はテンキーの数字と四則演算記号などで，出力は14桁のニキシー管による。プログラムステップは30とあるが，メモリに４つ使うため実質は26ステップだったのではないか。楠先生はこの４メモリあるのがすごい（値段が高い原因でもある）と強調していた。アセンブリ言語のような命令を使うのだが，フラッグへの条件分岐を多用する。与えられた課題の一つが与えられた自然数の階乗を計算するというものだった。ようやくできて喜んでいたら，同級生の土田純一君（物理部）がクラスで最小のステップ数で実現したと楠先生に褒められていた。残念，はやくも負けてしまったのだった。

その後，高校３年になると新しいオリベッティの卓上計算機が１台導入された。たぶん100ステップぐらいの計算ができ，メモリも10くらいになっていたのではないかと思うが，入試モードの我々は姿を見ただけで使わせてはもらえなかった。その後，大学に入り，大学院に進学する頃には自分の手元にAL-1000の機能を遥かに上回るカシオのプログラミング電卓がやってくることになる。時代はあっという間に進んでしまった。

写真：CASIO AL-1000 の雄姿（Vintage Calculators Web Museumから引用）

特集−量子コンピュータ（５）

特集−量子コンピュータ（４）からの続き

１７．河西棟馬（技術史）
　　　情報概念の形成（15p ☆☆☆）
　　　1920年代における物理学と工学の接近
科学的な情報概念の誕生の契機として，統計学，通信工学，物理学（統計力学）をあげている。大阪教育大学に教養学科を設置したとき，情報科学専攻が設けられ，そのカリキュラム試案が提示されたことがある。木立先生が，なんでこんなところに統計力学が入っているんだとディスっていて削除されたのだが，たぶん，どこかの大学をモデルにしていたのだろう。まあ，入っていると物理の誰かが担当しなければならないので，正解だったかもしれないが，後に物理出身の藤田修先生が着任されたので，入っていても良かったかな。河西さんはハートリーやナイキストに通信工学における定量的情報概念の形成の出発点を見出し，シラードエンジンからシャノンまでの歴史を辿っている。

１８．鈴木真奈（論理学／科学史）
　　　ホームオートメーション再考（12p ☆☆）
　　　1980年代の日本が描いた21世紀の情報化社会
たぶん，本書のテーマから最も距離のある論考だ。確かに，2010年代のIoTと1980年代のホームオートメーションは同じなのだ。奈良県の東生駒でCATVが始まったのはよくおぼえているし，それがゆえに，自宅のネットワーク回線はいまだにKCNだったりする。ISDNも使っていたが，一瞬でなくなってしまった。CDの誕生から消滅も一瞬といえるのだろうか。1990年代に大学のネットワークを含んだ将来構想を考えていたときは，この1980年代の古いマルチメディア社会をベースにしながら，それを越えるものとしてイメージしていたはずだ。

１９．加藤夢三（日本近代文学）
　　　偶然性・平行世界・この私（11p ☆☆☆）
　　　量子力学と文学をめぐる諸問題
東浩紀のクオンタム・ファミリーズはおもしろくなかった。東浩紀が嫌いだからではない（津田大介の件では最低でしたね）。文学的な量子飛躍は，文学の標準機能として備わっているものだから，クオンタムと称して物語を展開すること自身がなんだかなあになってしまう。さて，偶然文学論の中河輿一は古い人なのかと思ったが，まだ著作権は切れていない天の夕顔の人だった。まったく議論は深まっていかないのだけれど，量子情報が本質的な役割を果たすサイエンスフィクションを読んでみたい。なお，最初にＳＦの多世界解釈を見たのは，旺文社の学年別学習雑誌の付録のＳＦ読み物だ。ヒューゴー・ガーンズバックのラルフ１２４Ｃ４１＋もエドモンド・ハミルトンのフェッセンデンの宇宙もエドワード・スミスの宇宙のスカイラークもこの短い小冊子で読んだ。そこに，著者もタイトルも忘れたけれど，自動車事故で夫が死ぬ世界と妻が死ぬ世界が交錯する物語を読んだのをはっきり憶えている。

感染症の数理シミュレーション（２）

感染症の数理シミュレーション（１）からの続き

動機
政府は，NHKテレビを通じて「ここ１，２週間が瀬戸際」で「感染のピークシフトを」と訴えている。ところが，具体策としては，多少具合が悪くても自宅で寝ておけというもので，それ以上の（中国のような）大胆な政策出動には及び腰だ。北海道や和歌山県などでは感染者が報告されてきたが，大阪や京都で感染者がほとんどいないのはどうみても不自然だ。SNSのうわさでは，保健所や病院でたらい回しにされながら検査を拒否されているという。ごく一部の例で感染拡大をアピールしながら，厚生労働省の検査基準と各地方自治体の保健所の連係プレーで，医療資源保護や東京五輪のために検出の押さえ込みを図っていると見えてしまう。

ところで，日本の医療資源の限界とはどの程度のもので，感染によってこれがどの程度枯渇することが想定されているのか，あるいは，ピークシフトの時間スケールはどうなのかという半定量的な情報はほとんど出されていない。重症化した感染者が病院で病床について治療されることを考えてみる。厚生労働省の医療施設動態調査（平成30年2月）によれば，一般診療所を除く日本全国の病院の療養病床を除く病床数は約123万床なので国民100人に１床の割合だ。なお月末病床利用率というのも厚生労働省から報告されており，全国平均で約80%である。したがって，コロナウィルスの対応可能病床数は100人に0.2床ということになる。実際にはすべての病床が対応可能ではないから，100人に0.1床（１万人に10人）を上回る重症者が発生すれば危険というオーダーか。

そこで，今回の新型コロナウイルス感染症の特徴を満たすSEIRをもとにした数理モデルを作成して，重症者や死亡者のシミュレーションを行ってみる。なお，PFNの丸山宏さんが，林祐輔さんが公開したSEIRモデルを参考に，ダイアモンド・プリンセスにおけるCOVID-19発症日別報告数を観測データとして，最適化ツールOptunaを用いてパラメターフィッティングを試みている。丸山さんは，SIRモデルを仮定し（感染しない潜伏期間がないもの），N=3711人の集団に対し，初期の感染者（発症していない人を含む）I0，感染率β（Nで割っていない），感染者がある日に発症する率α，発症期間γをパラメータ推定した。発症率が4%で，初期感染者が209名，一人当たり感染数（Nβ）は41だ（βの定義に注意）。

モデル
SEIRモデルを少し修正する。免疫を持たないものが感染する部分は同じである。この軽症感染者は一定期間をへて，回復者又は重症感染者にいたる。重症感染者は，一定期間をへて回復者または死亡者にいたる。重症感染者は病院または自宅で隔離されるためそれ以上の感染源とはならない。死亡者も同様だ。軽症感染者は無症状かもしくは軽い風邪程度の症状なので普通に活動し，免疫非保持者を感染させる可能性がある。この重症感染者の数を病床数と比較することで，医療資源の破綻の程度が判断できる。

図　SEIRモデルを修正したSIIDRモデルの概念図

$\dfrac{du_1}{dt} = -\beta \dfrac{ u_1*u_2}{n}$

$\dfrac{du_2}{dt} = \beta \dfrac{ u_1*u_2}{n} -\dfrac{u_2}{\alpha_1} -\dfrac{u_2}{\alpha_2}$

$\dfrac{du_3}{dt} = \dfrac{u_2}{\alpha_2} -\dfrac{u_3}{\gamma_1} -\dfrac{u_3}{\gamma_2}$

$\dfrac{du_4}{dt} = \dfrac{u_3}{\gamma_2}$

$\dfrac{du_5}{dt} = \dfrac{u_2}{\alpha_1} +\dfrac{u_3}{\gamma_1}$

計算
$n$=10,000人の集団を考える。免疫非保持者$u_1$が9,999人，軽症感染者$u_2$が1人の初期状態から約１年間の推移をシミュレートする。軽症感染者の感染期間を14日として，その96%が回復し，4%が重症化すると仮定する。すなわち，$\alpha_1=14/0.96$ ，$\alpha_2=14/0.04$。さらに，重症感染者の感染期間を7日として，その90%が回復し，10%が死亡するとする。感染者の死亡率は0.4%になる。すなわち，$\gamma_1=7/0.90$ ，$\gamma_2=7/0.10$ 。感染率$\beta$をパラメタとして，0.1から0.3まで変化させた。ここで得られる結果$u_1, \cdots u_5$ はすべて集団１万人あたりの人数になる。そこで，重症感染者数 $u_3$と１万人当たりの受け入れ可能病床数（10床）を比較して，医療資源の枯渇の様子を調べる。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
using DifferentialEquations
using ParameterizedFunctions
using Plots; gr()

sky = @ode_def SIIDR_model begin
  du1 = -β*u1*u2/n # u1:Noimmunity(Susceptible)
  du2 = β*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
end n α1 α2 β γ1 γ2

n=10000.0 #total number of population
α1=14.0/0.97 #latent to recovery (days/%)
α2=14.0/0.03 #latent to onset (days/%)
β=1.0 #infection rate
γ1=7.0/0.90 #onset to recovery (days/%)
γ2=7.0/0.10 #onset to death (days/%)

function epidm(n,α1,α2,β,γ1,γ2)
u0 = [9999.0,1.0,0.0,0.0,0.0] #initial values
p = (n,α1,α2,β,γ1,γ2) #parameters
tspan = (0.0,360.0) #time span in days
prob = ODEProblem(sky,u0,tspan,p)
sol = solve(prob)
return sol
end

β=0.3
@time so=epidm(n,α1,α2,β,γ1,γ2)
println(β," ",so[4,45]," ",so[3,19])
plot(so,vars=(0,3))
β=0.2
@time so=epidm(n,α1,α2,β,γ1,γ2)
println(β," ",so[4,40]," ",so[3,19])
plot!(so,vars=(0,3))
β=0.15
@time so=epidm(n,α1,α2,β,γ1,γ2)
println(β," ",so[4,34]," ",so[3,20])
plot!(so,vars=(0,3))
β=0.12
@time so=epidm(n,α1,α2,β,γ1,γ2)
println(β," ",so[4,30]," ",so[3,22])
plot!(so,vars=(0,3))
β=0.1
@time so=epidm(n,α1,α2,β,γ1,γ2)
println(β," ",so[4,28]," ",so[3,24])
plot!(so,vars=(0,3))
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

結果
上記モデルを用いたピークシフトのふるまいは，次のグラフのような結果となる。

図：１万人当たりの重症感染者数の感染率β依存性，左からβ=0.30, 0.20, 0.15, 0.12, 0.10

重症感染者数$u_3$のグラフから得られた結果と死亡者数$u_4$の値をまとめると次のようになる。

表：感染率βと１万人当たりの重症感染者のピーク及び１年後の１万人当たり死亡者数
　β　　ピーク日　ピーク病床 　死亡者
　0.30 　 　　 50 　　　54　　 　30
　0.20  　　　 80 　　　38　　 　28
　0.15 　　　125　　　 23　 　　25
　0.12 　　　185　　　 13　 　　20
　0.10 　　　280　　　 7　 　　13

重症感染者のふるまいの感染率$\beta$への依存性は非常に大きい。$\beta$が0.12を越えると，必要なピーク病床数は簡単に受け入れ可能病床数の10を超えてしまい，医療パニックが発生する。これが，感染開始１ヶ月から半年の間に起こる（当然オリンピックの時期を含む）。非常に激しければ短期（３ヶ月）で収束するが混乱は限度を超える。感染率を押さえこんだとしても影響が長期化することは避けられない。

付記
肺炎による死亡は，悪性新生物や心疾患についで日本人の死因の第３位になっており，年々死亡率（10万人当たり死者数）が上昇している。１万人当りに換算すると，年に10人が肺炎で死亡していることになる。上記のβ=0.10の場合の13人はこれと同じオーダーであることから，もし中間的な過程が隠蔽されたとしても，１年後には明らかに新型コロナウィルス肺炎による肺炎死亡者の急増が観測されるはずである。

感染症の数理シミュレーション（３）に続く