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

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#!/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日火曜日

感染症の数理シミュレーション(9)

感染症の数理シミュレーション(8)からの続き

湖北省(武漢以外)を説明するSIIDR2モデルができたので,これを用いて典型的な2つのイメージを定量化して表現してみたい。以下の計算ではモデルを単純化して特徴を表わすために,感染確率の$\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}$と想定する。これは,厚生労働省の専門家会議における「新型コロナウイルス感染症の流行シナリオ」に準拠している。

A.短期集団免疫シナリオ
SARS-CoV-2の感染者が完全に免疫を獲得できるのかどうかがよく知らないのだけれど(まれに再感染するのか,単なる変異したウイルスへの感染か),英国のジョンソンドイツのメルケルが,国民の6〜7割が感染する可能性について言及しているようだ。果たしてそれが可能なのだろうか。SIIDR2モデルで$\beta=0.6,\ \nu=0.01,\ \lambda=10000,\ \tau=0,\ T=180$とする。$\lambda$を十分大きくとったのは$\bar{\beta} \approx \dfrac{8}{15}\beta$(一定)になるようにするためだ。


図1 短期集団免疫シナリオ(1万人当り)

グラフの基点は新規感染者累計が人口の1ppmに達した時点であり,感染者が増大している多くの国では条件を満足している。日本に当てはめると,ピークは3〜4ヶ月後(5〜6月)$u_1$:未感染数,$u_2$:軽症感染数(含む未発症),$u_3$:重症感染数(ピーク時400万人),$u_4$:死亡数(終期は60万人),$u_5$:免疫獲得回復数(終期は7000万人),$u_6$:重症感染数累計(終期は1400万人)である。したがってピーク時の重篤用必要設備(病床)数は20万床となり,完全に国内のキャパシティを越えてしまうものと思われる。当然オリンピックは不可能だ。

B.ピークシフトシナリオ
医療崩壊の本質は,PCR検査の数の問題でも,病床数の問題でもなく,重篤患者に対応するための設備と医者の数の問題であるといわれている。日本の厚生労働省がNHKを通じて広めている図も,ピークを低くしてその頂上の位置を先送りにする必要があるというものだった。SIIDR2モデルで$\beta=0.42,\ \nu=0.01,\ \lambda=10000,\ \tau=0,\ T=750$としたものが次の図である。図1とスケールのオーダーやファクターが違うことに注意しよう。

図2 ピークシフトシナリオ(1万人当り)

感染数のピークは10〜16ヶ月後,$u_2$:軽症感染数(含む未発症)(ピーク時50万人),$u_3$:重症感染数(ピーク時30万人),$u_4$:死亡数(終期は10万人),$u_6$:重症感染数累計(終期は300万人)である。したがってピーク時の重篤用必要設備(病床)数は1.5万床となり,かろうじて対応できるのかもしれないがよくわからない。大きな問題は2年にわたる長期の対応が必要なことである。ピークシフトを考える場合は,感染率$\beta$の抑制策との合わせ技が必須だ。

[1]ピークカット戦略(集団免疫戦略)地獄への道は善意で舗装されている(Sato Hiroshi)

2020年3月16日月曜日

日本の歴史

4月12日まで期間限定無料提供されている,小学館版学習漫画少年少女日本の歴史全24巻(22巻+付録2巻)をざっと読み終えた。第22巻の平成編はそれまでのものとトーンが違っており,全く面白くなかった。それ以外はとてもよかった。権力の移行過程に加えて,文化的なトピックスや民衆の動きなどが適当なバランスで記述されており,1巻の中を3〜4つのテーマに分けながら進めていくスタイルはとても理解しやすい。南京大虐殺や三光作戦や慰安婦などもちゃんと触れており,夏目漱石のところには米山保三郎も登場していた。第21巻あたりでは著作権の関係か背景がグレーに処理されていて新聞などの紙面?が再現されていないものがあるのがちょっと残念だった。

歴史といえば,大学時代に井上清(1915-2001)の岩波新書「日本の歴史(上・中・下)」を読んでいたことを思い出した。これもなかなか貴重な読書体験だった。史実の正確性や著者のマルクス主義的歴史観に対してはいろいろ批判はあるようだが,当時はとても新鮮に感じられた。

歴史といえば,小学校高学年のころか,父親が平泉澄(1895-1984)の「少年日本史」をどこからかもらってきて,これを読めと渡されたことがある。パラパラと見たところ直観的にこれはインチキだとすぐわかったのでほとんど読まなかった。父にこれ(皇国史観)はおかしいというと,叱られたような気がする。その頃は家にあった6巻ものの漫画ではない学研かどこかの日本の歴史を読んでいてこちらの方はまともでおもしろかった。


2020年3月15日日曜日

感染症の数理シミュレーション(8)

感染症の数理シミュレーション(7)からの続き

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

感染症の数理シミュレーション(3)において,感染対策時間因子を修正した図1のSIIDR2モデルで,湖北省(武漢以外)の新規感染者数累計と死亡数累計説明できるパラメタセットを得ることができた
図1 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}
湖北省(武漢以外)のデータを再現するパラメタセットは以下のとおりであり,これを用いたシミュレーション結果は図2のように与えられる。
$\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)$に対応している。


図2 湖北省(武漢以外)の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)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

[1]Coronavirus disease (COVID-2019) situation reports
[2]新型環状病毒肺炎疫情動態(tencentニュースサイト)
[3]新型コロナウイルス感染症対策専門家会議(第5回)(3月2日)
[4]特設サイト:新型コロナウイルス(NHK)
[5]新型コロナウイルス感染速報(Su Wei)
[6]Coronavirus Disease (COVID-19) – Statistics and Research(Our World in DATA)
[7]Databrew's COVID-19 data explorer(Databrew)
[8]COVID-19情報共有(佐藤彰洋)
[9]時間遅れを考慮した確率 SIR モデルの安定性解析(石川昌明)
[10]分布的時間遅れをもつ確率感染症モデルの安定性解析(石川昌明)

P. S. [6]のように,プロがシミュレーションをすると遅延SIRモデルみたいなことになるのか。感染率にステップ関数を導入しているので鋭い構造が出現しているようだ。遅延確率SIRモデルになるとちょっとついていけません(3/17/2020)。


感染症の数理シミュレーション(9)に続く

2020年3月14日土曜日

感染症の数理シミュレーション(7)

感染症の数理シミュレーション(6)からの続き

動機
岩田健太郎さんが,新型コロナウイルス感染症対策専門家会議(第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)$と重ねたものを図1に示す。

図1 ピークと変曲点の時点が$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$が得られるので,条件を満足している。

図2 $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$としてきた。ところで,各国の感染者比(1)各国の感染者比(2)から新規感染者数累計が人口の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}


感染症の数理シミュレーション(8)に続く


2020年3月13日金曜日

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

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

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

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

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


2020年3月12日木曜日

各国の感染者比(2)

各国の感染者比(1)からの続き

新型コロナウイルス感染症の話。
感染者数が一定の数あるいは人口比になった日時とその値を共通の原点として,各国の感染者数や人口比を対数グラフに表して観察する方法についてSNS上であれこれ話題になっている。最初の指摘は「他国と比べて日本だけが特異な振る舞いを示しているのは人為的な要因(検査の絞り込み等)の結果ではないか」というものであった。

これに対して,感染者数基準のプロットではほとんど重なって見えた他国・地域については,人口比にすれば分離して見えることを指摘した。それにもかかわらず,中国,湖北省(武漢以外),韓国,日本,イタリア,イランの比較では日本だけが異なった傾向を見せた(緩慢な指数関数型増加)。一方,症例の少ない多くの国を含めて比較すれば,日本だけが特異なのではないというクレームがついていた。ただし,無原則にあれもこれも混ぜてしまうのは問題の本質を見失う危険性があるので正しくない方法論だ。

そこで,アジア太平洋地域で一定の感染者が報告されている国々と日本,韓国について前回と同様の分析を行った。対象は,韓国,日本,台湾,香港,シンガポール,オーストラリアの6カ国とした。なお,今回の追加データは,ジョンズ・ホプキンズ大学のアーカイブから取得している。各国の人口,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の常用対数)

(1)韓国と台湾は,それぞれ異なった水準ではあるが収束に向かう様子が窺える。

(2)欧州の感染者急増国の多くは韓国と相似形の振る舞いをしている。

(3)取り上げたアジア太平洋諸国(韓国,台湾以外,日本を含む)は対数グラフでおおよそなだらかな1次関数で増加している(すなわち,底の小さな指数関数的増加)。したがって,まだ収束の様子はみられない。

(4)香港とシンガポールは1次関数のスロープが原点から2週目以降でよりゆるやかに変化している。なんらかの抑制策が効果を発揮したのかもしれない。シンガポールは0日から18日まで平均15.5%,18日から42日まで平均3.5%で増加,香港は0日から15日まで平均11.0%,15日から43日まで平均4%で増加している。なお,日本は,0日から16日まで平均9%とその中間になっている。

2020年3月11日水曜日

ワン・モア・ヌーク:藤井太洋(2)

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

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

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

2020年3月10日火曜日

各国の感染者比(1)

各国の新規確定感染者数の累計(Confirmed)が100名または200名となった時点を1としてその前後の感染者数累計を対数プロットしたものがtwitter上でいくつか紹介され,日本だけが異なった振る舞いを示していることが指摘されている。日本の感染者数の統計が不自然ではないかということはこれまでも話題になってきた。

ところで国によって人口は大きく異なるのにもかかわらず基準を100名などの固定した人数にするのはどうかと思う。そこで,新規確定感染者数の累計(Confirmed)を各国/地域の人口で割った感染者比が1ppm(百万分の一)となるところを基準にしてそれらの時間経過を比較してみよう。データはWHOのCoronavirus disease (COVID-2019) situation reportsを用いる。感染者数が多い次の6つの国や地域(日本を含む),(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の常用対数)

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

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

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

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


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

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

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

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

(注5)1/30 -> 1/21 に訂正(グラフにはまだ反映されていない 3/15/2020)

[1]感染症の数理シミュレーション(6)

各国の感染者比(2)に続く

2020年3月9日月曜日

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のシミュレーション部分はミスリーディングである。

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

2020年3月8日日曜日

感染症の数理シミュレーション(6)

感染症の数理ミュレーション(5)からの続き

湖北省(武漢以外)のデータでは,新規/累計感染者数 f(t), g(t) と新規/累計死亡数 h(t), i(t) を簡単な近似関数で表現した。これを用いると,感染症の数理シミュレーション(3)の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人(1万人当り)[4]。
⑥同上の新規感染者数f(t)は12日のピークに0.25/日,40日目の累計で3.70人(1万人当り)[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)→感染症の数理シミュレーション(7)

各関数の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)→感染症の数理シミュレーション(7)
\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)$を上式に代入したものと,感染症の数理シミュレーション(5)で考えたモデル感染対策時間因子で$\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$と矛盾しないような感染対策時間因子のパラメタセットが存在することがわかる。

[1]新型コロナウイルス感染症対策の見解(厚生労働省)
[2]新型コロナウイルス(COVID-19)感染症(川名・三笠・泉川)
[3]新型コロナウイルス肺炎COVID-19(下畑享良)
[4]tenetニュースサイト
[5]Coronavirus disease (COVID-2019) situation reports(WHO)
[6]COVID-19(奥村晴彦)

感染症の数理シミュレーション(7)へ続く

2020年3月7日土曜日

感染症の数理シミュレーション(5)

感染症の数理シミュレーション(4)からの続き

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

感染症の数理シミュレーション(3)では,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$で割って,感染対策時間因子だけを取り出したものをグラフにすると次のようになる。横軸はシミュレーション開始からの日数を表す。
図1 感染対策時間因子β(t) λ=5 左からτ=5,10,15

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

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






2020年3月6日金曜日

ワン・モア・ヌーク:藤井太洋(1)

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

2020年3月5日木曜日

湖北省(武漢以外)のデータ

新型コロナウイルス感染症の事例数が最も多い地位は中国の湖北省であるが,武漢のデータは特異的な値を示している(高い死亡率)。そこで,モデルの妥当性を検討するためには湖北省(武漢以外)のデータを用いるのが適当だ。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 など,観測データをおおむね再現している。

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


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

湖北省(武漢以外)での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日目にはおおむね収束に近づいている。これをグラフにすると次のようになる。

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

2020年3月4日水曜日

あぶり餅

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

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



写真:あぶり餅の一和(2020年3月3日撮影)




2020年3月3日火曜日

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

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

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

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

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


2020年3月2日月曜日

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

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

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

まず,生活時間を調べよう。5年おきに実施している最新の平成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人でいるものとする。そこで感染機会Cに寄与するのは複数で過ごす時間である。
C(家庭) = 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
職場は非常に多様な環境であり,えいやっと1人1坪のオフィスや工場を考える。
C(職場) = 0.47*1人/3.3㎡*8時間 = 1.14

通勤通学時間は,NHKによる2015年の国民生活時間調査と国勢調査が参考になる。徒歩(7%),自転車等(15%),自家用車(47%),電車バス等(31%)である。また,電車の1人当り立席専有面積は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 だけ減少するので,ほぼ全校休校措置と同程度の効果をもたらすことになる。かもしれない。

[1]新型コロナウイルス感染症対策本部(首相官邸)
[2]新型ウイルス研究者試算 対策組み合わせで入院患者6割以上減(NHK)


2020年3月1日日曜日

感染症の数理シミュレーション(4)

感染症の数理シミュレーション(3)からの続き

新型コロナウイルス感染症に対するWHOなどの共同調査報告書[1]が公表された。同時に,新型コロナウイルス感染症の世界的なリスクは4段階の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)とは,特定の疾患に感染した母集団に対する死亡数の割合である。ここで,疾患に感染していることが判明している必要があるので,今回のようにほぼ無症状や他の疾患と区別できない軽症の感染者がいる場合には取り扱いに注意がいる。

感染症の数理シミュレーション(3)では,東大の黒田産業医のデータをもとにパラメタを定めたが,今回の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$の3要素が必要である。いずれにせよ,合理的な判断ができなくて肝腎の検査がシステマティックに行われず,感染状況に関するデータが完全に欠落している日本では推定が難しい。

[1]Report of the WHO-China Joint Mission on Coronavirus Disease 2019 (COVID-19)

2020年2月29日土曜日

昭和のプログラミング教育(4)

昭和のプログラミング教育(3)からの続き

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

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

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

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


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


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

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

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


2020年2月28日金曜日

感染症の数理シミュレーション(3)

感染症の数理シミュレーション(2)からの続き

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

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

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

東京大学の黒田玲子産業医が個人的な見解としてまとめたメモによれば,新型コロナウイルス感染症の潜伏期は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%とは矛盾しない。こうして4つのパラメタは合理的な範囲で収められる。

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

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

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

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

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

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

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

図3 小さな山/丘(λ=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]のように,中国からはよりまともな遅延ダイナミカル・シミュレーションの結果が報告されている。日本では誰か計算していないのだろうか・・・

[1]緊急寄稿新型コロナウィルス(COVID-19)感染症(川名明彦・三笠桂一・泉川公一)
[2]2019-nCoVについてのメモとリンク(中澤港)
[3]人口と感染症の数理(稲葉寿)
[4]数理モデルによる感染症流行制御の検証(大森亮介)
[5]新興感染症の国際的伝播を予測する数理モデル(西浦博・木下諒)
[6]CoVID-19 in Japan: What could happen in the future?(medRxivから)
[7]The reproductive number R0 of COVID-19 based on estimate of a statistical time delay dynamical system(medRxivから)
[8]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))
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

感染症の数理シミュレーション(4)に続く

2020年2月27日木曜日

昭和のプログラミング教育(3)

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

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

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


写真:CASIO AL-1000 の雄姿(Vintage Calculators Web Museumから引用)

2020年2月26日水曜日

特集−量子コンピュータ(5)

特集−量子コンピュータ(4)からの続き

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

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

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

2020年2月25日火曜日

感染症の数理シミュレーション(2)

感染症の数理シミュレーション(1)からの続き

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

ところで,日本の医療資源の限界とはどの程度のもので,感染によってこれがどの程度枯渇することが想定されているのか,あるいは,ピークシフトの時間スケールはどうなのかという半定量的な情報はほとんど出されていない。重症化した感染者が病院で病床について治療されることを考えてみる。厚生労働省の医療施設動態調査(平成30年2月)によれば,一般診療所を除く日本全国の病院の療養病床を除く病床数は約123万床なので国民100人に1床の割合だ。なお月末病床利用率というのも厚生労働省から報告されており,全国平均で約80%である。したがって,コロナウィルスの対応可能病床数は100人に0.2床ということになる。実際にはすべての病床が対応可能ではないから,100人に0.1床(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人の初期状態から約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$ はすべて集団1万人あたりの人数になる。そこで,重症感染者数 $u_3$と1万人当たりの受け入れ可能病床数(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))
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

結果
上記モデルを用いたピークシフトのふるまいは,次のグラフのような結果となる。
図:1万人当たりの重症感染者数の感染率β依存性,左からβ=0.30, 0.20, 0.15, 0.12, 0.10

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

表:感染率βと1万人当たりの重症感染者のピーク及び1年後の1万人当たり死亡者数
 β  ピーク日 ピーク病床  死亡者
 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を超えてしまい,医療パニックが発生する。これが,感染開始1ヶ月から半年の間に起こる(当然オリンピックの時期を含む)。非常に激しければ短期(3ヶ月)で収束するが混乱は限度を超える。感染率を押さえこんだとしても影響が長期化することは避けられない。

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

感染症の数理シミュレーション(3)に続く


2020年2月24日月曜日

感染症の数理シミュレーション(1)

新型環状病毒肺炎からの続き

あの神戸大学の岩田健太郎さんが新型コロナウィルス感染症の数理シミュレーションの論文を書いていた[1]。調べてみると基本的なSEIRモデルを使った簡単な微分方程式系だ。感染症に対して免疫を持たない者(Susceptible),感染症が潜伏期間中の者(Exposed),発症者(Infectious),感染症から回復し免疫を獲得した者(Recovered)から構成され,その頭文字をとってSEIRモデルと呼ばれる。このうち集団全体の出生率と死亡率を0に置いたモデルで計算していた。

このモデルの詳細については,Compartmental models in epidemiologyに詳しいが,ここでは上の日本語版WikippediaのSEIRの記事で十分だ。早速,Juliaで追計算してみた。下の例では,n=1000人の集団に1人の感染者がいるところからはじまる。感染率βは1日に1人で,潜伏期間αは12日,発症期間γは8.5日である。このモデルでは潜伏期間には感染しない。コロナウイルスは潜伏期間にも感染するようなので,この部分は簡単にモデルを修正できる。以下のコードではδとして導入した項だ(ピークがかなり前倒しになる)。

αとγは制御できないので,感染率βを如何にコントロールするかが対策のキモとなる。潜伏期間に感染するδ=1.0のモデルでは,β=0.2で集団の9割以上が感染してしまうが,β=0.1で8割,β=0.07では1年かけて5割強,β=0.05では1年後にも1割以下となり,非常にセンシティブだ。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
using DifferentialEquations
using ParameterizedFunctions
using Plots; gr()

sky = @ode_def SEIR_model begin
  du1 = -β*u1*(δ*u2+u3)/n # u1:Susceptible
  du2 = β*u1*(δ*u2+u3)/n -u2/α # u2:Exposed
  du3 = u2/α -u3/γ # u3:Infectious
  du4 = u3/γ # u4:Recovered
end n α β γ δ

function dif()

# hayashi-model
#n=1000.0 #total number of population
#α=5.5 #latent period (days)
#β=1.0 #disease transmission rate
#γ=7.0 #infectious period (days)
#δ=0.0 #latent transmission rate

# iwata-model
n=1000.0 #total number of population
α=12.0 #latent period (days)
β=1.0 #disease transmission rate
γ=8.5 #infectious period (days)
δ=0.0 #latent transmission rate

u0 = [999.0,1.0,0,0] #initial values
p = (n,α,β,γ,δ) #parameters
tspan = (0.0,100.0) #time span in days
prob = ODEProblem(sky,u0,tspan,p)
sol = solve(prob)
plot(sol,vars=(0,1))
plot!(sol,vars=(0,2))
plot!(sol,vars=(0,3))
plot!(sol,vars=(0,4))
end

@time dif()
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
岩田さんが例にあげているアウトブレイクするパラメタを代入してみると,よく似たグラフが得られたが,発症者のピーク位置や感染者と発症者の山の重なりが微妙に違う。彼の論文ではパラメタの領域しか提示されていないため,厳密な比較ができない。そこで,林祐輔さんのシミュレーション[2]と比較した。こちらのほうとは,ピタリと一致したのでたぶん上のコードは正しいのだろう。

図 JuliaによるSEIRモデルの確認


[1]A Simulation on Potential Secondary Spread of Novel Coronavirus in an Exported Country Using a Stochastic Epidemic SEIR Model(K. Iwata, C. Miyakoshi)
[2]COVID-19/stochasticSEIRmodel.ipynbY. Hayashi
[3]面向新冠疫情的数据可视化分析与模拟预测(陈宝权 北京大学前沿计算研究中心)
[4]corona_virus_model fitting_by optina(H. Maruyama)

 感染症の数理シミュレーション(2)に続く

2020年2月23日日曜日

特集−量子コンピュータ(4)

特集−量子コンピュータ(3)からの続き

13.大黒岳彦(哲学)
   量子力学・情報科学・社会システム論(22p ☆☆☆☆)
   量子情報科学の思想的地平
著者は次のようにいう。量子情報科学(量子計算−量子暗号−量子通信の三位一体)は次世代の情報社会(N. ルーマンの社会システム論=コミュニケーションの連鎖的接続)の要をなす技術であり,パラダイム変革的ではなくパラダイム完成的テクノロジーであると。その後,情報科学の誕生から情報科学と現象学の対比をへて,物的世界像を逆転させた情報的世界観を提示する。さらに,物理学におけるモノと情報の相克が量子力学に現れることから,量子情報科学の誕生からその思想的地平までに及ぶ。哲学者の議論は往々にして深い迷い道に誘い込まれるのだが,大黒さんはとても明解である。ルーマンの社会システム理論への過剰な肩入れが気になるのだが(ここは自分の勉強不足か?),それを除けばとても実り多い論考だった。

14.河島茂生(情報倫理)
   未来技術の倫理(13p ☆)
「倫理(ethics)」とはもともと「慣習(ethos)」を表す言葉である。ここまでは良かったが,オートポイエティック・システム(自分自身で自分をつくるシステム)やアロポイエティック・システム(他のものによってつくられ他のものをつくるシステム)あたりからついていけなくなり,転調してから,社会−技術システムと功利主義的計算が何なのかさっぱり分からず,すっ飛ばしながら終了した。チーン。

15.斉藤賢爾(インターネットと社会)
   デジタル署名,ブロックチェーンと量子アルゴリズム(10p ☆☆☆)
   新たなるサイファーパンクの夜明け
フィンテックが専門の斉藤賢爾さんは以前から独特の発想がおもしろくてときどきウォッチしてきた。「サイファーパンク」とは「強力な暗号技術やプライバシー強化技術を社会的・政治的変化への道として広く利用することを提唱する活動家」のことであり,斉藤さん自身のことかもしれない。暗号技術の本質は記録の真正性を保つことにあるとして,ブロックチェーンの話につなげていく。ただ,量子計算との関係はあまり深まらなかったのが少し残念だった。

16.羽根次郎(中国近現代史)
   チャイバースペース(Chyberspace)の出現について(10p ☆☆)
   中国の「サイバー主権」論の背景にあるもの
話は2000年の九州・沖縄サミットの沖縄IT憲章から始まる。結論は,サイバースペースの「国家からの独立性」が存在しないのは,中国のみに限定される話では最初からない,である。もうひとつのキーポイントは,次の部分だ。確率上稀にしか起せないイノベーションを,「ベンチャー支援」という名の「失敗時のコスト最小化作業」を通じてまで支援する所以が・・・(割に合わないイノベーションにかくも期待するのは)・・・イノベーションを起すための不可欠な条件として新自由主義的諸改革への白紙委任を正当化するためだ,というところ。

2020年2月22日土曜日

昭和のプログラミング教育(2)

昭和のプログラミング教育(1)からの続き

楠先生から,北國新聞2020年2月18日の紙面コピーが送られてきた。タイトルは「昭和40〜50年代のプログラミング,最先端教育記憶つなぐ,元泉丘高教諭・楠さん 教え子から感想文収集」であった。楠先生は,93歳で,1958〜77年に泉丘高で勤務しているので,我々が高校生だった50年前には,43歳ぐらいだった。

クスリのアオキの青木桂生さんが,日本経済新聞の2016年2月22日の記事「パチ先生の教え」で,楠先生(あだながパチ)の思い出と,実家の薬屋を継ぐように諭されたことを書いている。青木さんは1942年生まれくらいなので,我々よりほぼ10歳上,楠先生が30歳で泉丘高校に着任したころの生徒だ。クスリのアオキの会社沿革には,1985年に株式会社クスリのアオキを設立したときの本社所在地は,金沢市泉野出町4丁目になっている。泉丘高校のすぐ南になる。昔はほとんど田んぼだったが,円光寺線の道路沿いには,何かあったのかもしれない。

理数科2期生の我々8組のクラス担任は,1年が松川先生,2年が楠先生,3年が物理の谷口先生だった。松川先生は幼なじみの松川隆行君の父上で,東大の理学部数学科だったと思う。寺町2丁目の電車通りに面して,横山酒店のとなりが松川金物店で,泉野小学校の1,2年のときは同級生だった松川君のうちに時々遊びに行った。担任の川村先生に指名されて,2人で毎日クラス絵日記をかかされていた。その後,雪の市立工業高校のグラウンドをころげまわって帰っていた。その転がり方を○○流だといいあっていたのだが,彼が泉丘流だといったのが,はじめて「泉丘」という言葉を知ったときだ。

松川君にはお姉さんがいて,弟共々よくできるひとだった。お姉さんもお茶の水だったのではないか。3人でダイヤモンドゲームをして,自分がこてんぱんに負けた記憶がある。松川君は東大に進みその後,建設省だかどこかの高級官僚になったと思う。しばらく前は,日本木材住宅産業協会の専務理事の足跡がある。

かつての横山酒店の向いは,今は北國銀行の支店だが,当時は協栄のパン屋があった。その間の路地を入って100mほどで寺町のおばあちゃんの家に至る。途中に大杉和人君というのがいて,彼が泉野小学校の我々の学年で最も優秀だった。卒業式の答辞も読んだ。彼も東大に進学している。泉野小学校からは,他にも赤井勇一,西田荘平,山本賢正,山田晴美など計6人が東大に,片岡稔,山田良平,平崎鉄外喜など3人が京大に進んでいる。他にも忘れている人がいるかもしれない。昭和のプログラミング教育の話にたどりつかない。

P. S. 当時の泉野小学校は,1学年150人ほど,1年から5年までは3クラス,6年ではじめて4クラスになった。野田中学校は1学年450人で10クラス。

2020年2月21日金曜日

特集−量子コンピュータ(3)

特集−量子コンピュータ(2)からの続き

9.田中成典(計算生物学)
  量子と生命(11p ☆☆☆☆)
量子生物学というキーワードは1969年に大木幸介さんのブルーバックスを買って以来だ。最近注目を浴びているのは2007年の光合成系の量子コヒーレンスを契機としている。以前取り上げた鳥の磁気コンパスの話題もそれで,生物の感覚センサー機能における量子現象が注目されているとのこと。著者は量子多体系の第一原理計算についての素養があるので,生命現象の本質と量子計算の関りについて,生命量子計算や量子知性というキーワードで照射しようとしている。予想以上に興味深い内容だった。

10.郡司ペギオ幸夫(生命基礎論)
   「わたし」に向って一般化される量子コンピューティング(14p ☆☆)
元神戸大学の原俊雄先生にいわせると,郡司ペギオ幸夫さんは天才だが(なので?),普通の人の理解の枠を越えている。目次や表題をみると非常におもしろそうなのだが,読み出すとほとんどついていけなくなるのが常である。今回のキーワードは認知的非局所性だったが,そのスタートのロジック「XであってXでない」につまずいてしまった。

11.全卓樹(理論物理学)
   量子力学と現代の思潮(10p ☆☆☆☆)
全さんは東大の有馬研の出身だ。昔,森田研でポスドクをされていた。兵庫県の中山間部であった深山良徳君の葬儀の帰りに車で送っていただいたことなども思い出す。全さんの論説を読むためにこの雑誌を買ったのだ。期待を裏切らず,自分の知りたいと思っていた量子力学の「観測問題」をすっきりと整理してくれていた。東北大学の堀田昌寛さんは,コペンハーゲン解釈は量子力学の認識論的解釈であると言い放って,多世界解釈を簡単に斬って捨て,量子力学に観測問題などないと断定している。なんか違うのよね。哲学を含む人文科学一般について深い知識を持っている全さんは,スタンフォード哲学百科のミルヴォルトの整理に準拠して問題を解きほぐしてくれるところから始め,時代思潮と結びつけて着地する。この雑誌の特集の趣旨を十分に理解した論説だ。

12.内井惣七(科学哲学)
   無時間,無空間からの出発(8p ☆☆)
内井さんの本も何冊か挑戦しているのだけれど,まったくピンと来ないことが多い。なぜだろう。量子重力のロヴェリを最初に持ってくるところで引っかかるのかしら。いや,ライプニッツのモナドとの関連は確かにいいところをついていると思うのだ。でも確率のコード解釈がどうのこうのという結論へ向って,まったく筋が追えなくなってしまう。科学哲学の人の書いたもの,読める場合と読めない場合がある。認知の抗体反応を起してしまうのかもしれない。