アジアの状況からの続き
欧米の状況を見ると新規感染数累計は人口比で10ppmを越えているところがある。スペインの30ppmは湖北省（武漢以外）の10倍近い水準であり，下記の国々の新規感染数累計をはすべて湖北省（武漢以外）以上の値となっている。ただし，対数グラフ上は上に凸となっていて増加率は減少に向いつつある。
図１ 人口当りの新規感染数累計（単位100ppm，基準日は1ppm達成時）
図２ 人口当りの新規感染数累計の対数（単位100ppm，基準日は1ppm達成時）
2020年4月12日日曜日
アジアの状況
新型コロナウイルス感染症の感染者数の増加が5/6には収まるように考えている人が多いのかもしれない。うまくいけば７月には一端終息に向ったようにみえる可能性もある。しかし，集団免疫が獲得できずワクチンもない現状では，緊急事態宣言レベルの制限を継続するか，断続的に緩めたり強めたりすることの繰り替えしかの二択ではないだろうか。中国以外で終息に近い状態を実現しているのは台湾だけだ。それに近いのは韓国。香港もシンガポールも完全にはおさまっていない。どこまで耐えられることか。
図１ 人口当りの新規感染数累計（単位100ppm，基準日は1ppm達成時）
図２ 人口当りの新規感染数累計の対数（100 ppm，基準日は1ppm達成時）
注：上記は武漢を除いた湖北省の値であり，370ppmに収束している。武漢を含めた湖北省の収束値は1150ppmであり，上記の3.1倍に相当する。湖北省の全体イメージは湖北省（武漢以外）を全体に3倍程度スケールしたものと考えられることに注意する。
2020年4月11日土曜日
2020年4月10日金曜日
BCG
Wikipediaによれば，BCG とは次のようなものだった。
自分が小学生低学年の間は陰性が続いていて，ほぼ毎年のようにBCGを注射していたような気がする。1回くらいのBCGでは抗体ができなかったということか。で，BCGも初めのうちは今のようなスタンプ型ではなかった。それがいつの間にか（日本では1974年から）定期化され全員接種に切り替わっていた。
そのBCGが新型コロナウイルス感染症に効果があるかどうかということで議論になっている。乳児向けのBCGが不足することがないように祈るだけである。
BCGは，実験室で長期間培養を繰り返すうちにヒトに対する毒性が失われて抗原性だけが残った結核菌であり，BCGワクチンはBCGを人為的にヒトに接種して感染させることで，結核に罹患することなく結核菌に対する免疫を獲得させることを目的としたものである。小学生の時には，毎年ツベルクリン反応を調べる注射が恒例行事だった。ツベルクリン（独: Tuberkulin）とは，結核菌感染の診断に用いられる抗源である。前腕の内側に注射して数日後にこれが赤く腫れている（陽性）か変化なし（陰性）かを透明の物差しで測って確認する。たぶん1cm以上のサイズなら陽性なのでこれでOKだ。それ未満なら陰性や疑陽性，すなわち結核菌に対する抗体がないので，BCGを注射しなければならない。したがって小学生にとってはこれはなかなかドキドキする検査なのである。
自分が小学生低学年の間は陰性が続いていて，ほぼ毎年のようにBCGを注射していたような気がする。1回くらいのBCGでは抗体ができなかったということか。で，BCGも初めのうちは今のようなスタンプ型ではなかった。それがいつの間にか（日本では1974年から）定期化され全員接種に切り替わっていた。
そのBCGが新型コロナウイルス感染症に効果があるかどうかということで議論になっている。乳児向けのBCGが不足することがないように祈るだけである。
2020年4月9日木曜日
緊急事態宣言と接触制限モデル（３）
緊急事態宣言と接触制限モデル（２）からの続き
4月7日のtwitterで牧野淳一郎さんがいうには，
上に示されている数字は行動抑制の因子であり，モデルの感染率にかける係数だと思われる。ここでは，どれが正しいかについては考察せずに，ドイツなど欧州の基本再生産数$R_0=2.5$（7日増倍率が7倍）に基づいた西浦さんの計算（我々の評価では$R_0=2.3$だったが）を，東京の現状に合わせて再計算した結果を報告する。
前回述べたように，現在の東京では，7日増倍率が2.5倍（$R_0=1.65$）程度なので，西浦さんが用いている値よりかなり小さいのだ。単純なSIRモデルを用いてこれを再現するパラメタセットを探し，その場合に接触制限の効果がどうなるかを調べよう。
やり方は前回と同じなので，パラメタを示す。β=0.28，初期値はν=0.005とすれば上記の増倍率が再現できた。つまり，さきのMathematicaコードにおいて，
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
f[n_, β_, α_, ν_, p_]:=
NDSolve[{S'[t]==-(p+(10-p)*Tanh[2*(20-t)])/10*β/n*J[t]*S[t],
J'[t]==(p+(10-p)*Tanh[2*(20-t)])/10*β/n*J[t]*S[t]-J[t]/α,
R'[t]==J[t]/α, S[0]==n, J[0]==ν, R[0]==0}, {S,J,R}, {t,0,100}];
n=1400; β=0.28; ν=0.005; p=10;
sol = f[1400, β, 7, 0.005, p];
s[t_] := S[t] /. sol[[1, 1]];
i[t_] := J[t] /. sol[[1, 2]];
cft1[t_]:= (p+(10-p)*Tanh[2*(20-t)])/10*β/n*i[t]*s[t];
・・・
Plot[{cf1[t], cf2[t], cf3[t], cf4[t], cf5[t], cf6[t]},
{t, 10, 30}, PlotRange -> {0, 0.1},
GridLines -> {{19, 20, 21}, {0.005, 0.01, 0.015, 0.02, 0.025}}]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
などとすると計算結果は次のようになった。
4月7日のtwitterで牧野淳一郎さんがいうには，
一応どれも現実をモデルしているはずのところで、そこまで行動抑制しないといけない割合が、西浦氏:0.2、大橋氏「0.45](R_0>1でいくとして)、佐藤氏 0.02 で3人で20倍違う時点で少なくとも2人は間違っているのは明らかであろう。西浦氏は，緊急事態宣言と接触制限モデル（１）で取り上げた西浦博さん（北海道大学），大橋氏は，新型コロナウイルス感染症の 流行予測の大橋順さん（東京大学），佐藤氏は，COVID-19情報共有の佐藤彰洋さん（横浜国立大学）である。
上に示されている数字は行動抑制の因子であり，モデルの感染率にかける係数だと思われる。ここでは，どれが正しいかについては考察せずに，ドイツなど欧州の基本再生産数$R_0=2.5$（7日増倍率が7倍）に基づいた西浦さんの計算（我々の評価では$R_0=2.3$だったが）を，東京の現状に合わせて再計算した結果を報告する。
前回述べたように，現在の東京では，7日増倍率が2.5倍（$R_0=1.65$）程度なので，西浦さんが用いている値よりかなり小さいのだ。単純なSIRモデルを用いてこれを再現するパラメタセットを探し，その場合に接触制限の効果がどうなるかを調べよう。
やり方は前回と同じなので，パラメタを示す。β=0.28，初期値はν=0.005とすれば上記の増倍率が再現できた。つまり，さきのMathematicaコードにおいて，
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
f[n_, β_, α_, ν_, p_]:=
NDSolve[{S'[t]==-(p+(10-p)*Tanh[2*(20-t)])/10*β/n*J[t]*S[t],
J'[t]==(p+(10-p)*Tanh[2*(20-t)])/10*β/n*J[t]*S[t]-J[t]/α,
R'[t]==J[t]/α, S[0]==n, J[0]==ν, R[0]==0}, {S,J,R}, {t,0,100}];
n=1400; β=0.28; ν=0.005; p=10;
sol = f[1400, β, 7, 0.005, p];
s[t_] := S[t] /. sol[[1, 1]];
i[t_] := J[t] /. sol[[1, 2]];
cft1[t_]:= (p+(10-p)*Tanh[2*(20-t)])/10*β/n*i[t]*s[t];
・・・
Plot[{cf1[t], cf2[t], cf3[t], cf4[t], cf5[t], cf6[t]},
{t, 10, 30}, PlotRange -> {0, 0.1},
GridLines -> {{19, 20, 21}, {0.005, 0.01, 0.015, 0.02, 0.025}}]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
などとすると計算結果は次のようになった。
図 SIRモデル＋シグモイド関数制限措置による新規感染者数の変化（東京のデータ）
（上から順に，青：制限なし，80%，60%，50%，40%，20%に制限）
増加と減少の臨界点は50%あたりにある。すなわち，数字だけでいえば大橋氏の値と近い結果が得られたのかな。大橋氏の資料に対する牧野さんのもろもろの批判にはうなずけるものもあるけれど。そう，このブログも拡散してはいけない情報の仲間である。
まとめ
新型コロナウイルス感染症専門家会議（というか西浦博さん）から出てきた接触8割削減の前提条件は，欧州並の$R_0=2.5$に対応する場合ということだった。それを現在の東京の$R_0=1.65$に対応させると，上記のように接触6割削減で新規感染数は減少させることができる。いずれにせよモデル計算の結果なので，他地域の今後の動向や安全率も勘案して，目標値として接触8割削減を掲げることは意味があるだろう。
P. S. 佐藤彰洋さんのシミュレーションに対する疑義がでてきた。ちょっと安心した。Delay Differential Equationによるモデリングは，それでもまだ自分には理解できていない。
P. P. S. 4/10に牧野淳一郎さんの解説が出てきたので問題点がよく理解できるようになった。ただ，6割削減でも新規感染数が減少に転ずるという牧野さんの解説は，7割削減が必要であるというこちらの結果とは食い違っていた。
［１］人との接触７～８割減、効果は 専門家「感染抑制できる」、全員やれば６割減でも 新型コロナ（朝日新聞）→タイトルと本文があまり整合していない記事
［２］佐藤彰洋氏(横浜市大)の新型コロナ感染予測シミュレーションに関する疑義について（本堂毅・佐野雅己・松下貢）
［３］いろいろなモデル計算（牧野淳一郎，2020.4.10）
［４］交通整理（牧野淳一郎，2020.4.12）
［５］公開質問状に対するコメント（佐藤彰洋，2020.4.13）
［６］交通整理の続き（牧野淳一郎，2020.4.12）
［７］交通整理の続きその２（牧野淳一郎，2020.4.16）
［８］新型コロナウイルス感染症対策専門家会議の見解等（牧野淳一郎，2020.4.18）
2020年4月8日水曜日
緊急事態宣言と接触制限モデル（２）
緊急事態宣言と接触制限モデル（１）からの続き
緊急事態宣言がでた翌朝（4/8）の日本経済新聞にも再び西浦さんの図が「感染拡大阻止 接触８割減が必要」という記事とともに掲載されていた。安倍晋三も同内容を専門家の見解として強調していた。
緊急事態宣言がでた翌朝（4/8）の日本経済新聞にも再び西浦さんの図が「感染拡大阻止 接触８割減が必要」という記事とともに掲載されていた。安倍晋三も同内容を専門家の見解として強調していた。
図 日本経済新聞（2020.4.8 朝刊３面）から引用
日経の記事を批判的に読み解いてみる。
①「接触８割減」は政府の専門家会議メンバーで，感染者数の予測を数理モデルで解析している北海道大学の西浦博教授がはじいた数字だ 。
→ 西村さんは新型コロナウイルス感染症対策専門家会議の当初の構成員ではなかった。クラスター対策班のメンバーか。
②１人の感染者が平均で何人に感染させたかを示す「実効再生産数」は3月21〜30日の改定データで推定1.7。
→新規感染者数累計の1週間の増倍率を$k$とすると，実効再生産数$R_0$との間に，$R_0=1+5/7*\log k $の関係がある（倍加時間と基本再生産数）。東京（3/20-3/31）のデータは，平均で$k=2.3$であることから，$R_0=1.6$となった（注：$R_0=1.7$に対応するのは$k=2.7$）。
③その時点では数万人の感染者が出ていたドイツ（2.5）を下回っていたが，
→ドイツ（3/20-3/31）のデータは，平均で$k=3.9$であることから，$R_0=2.0$となった（注：$R_0=2.5$に対応するのは$k=8$，3/10以前には$k=8$を超えていたが・・・）。
④4月に入っても感染増が続き，実効再生産数は3を超えてドイツを上回った可能性があるという。
→実効再生産数がを3を超えるとは。$R_0=3$に対応するのは$k=16$である。日本でそのようなタイミングがあったのだろうか。福岡で9倍，福井で20倍くらいのタイミングは一時的に見られるが，これらはいずれも感染数が少なくて揺らぎが効いてしまう段階の話である。
⑤このままだと1日あたりの新規感染が米ニューヨークのように数千人に達する。
→仮定が④であれば正しいし，指数関数的な増大であれば常に正しい言明ではあるが，事態を強調するためにやや話を盛っている印象を受けた。
ただし，新規感染数累計について単純な指数関数型増加傾向が続くと仮定すると（接触制限の効果がない場合），直近の増倍率は東京で1.15/日，大阪で1.10/日となっているので，
東京：4/14：3千人，4/21：8千人，4/28：2万人，5/5：5万人
大阪：4/14：1千人，4/21：2千人，4/28：4千人，5/5：8千人
である。1日あたりの数ではなく累計であることに注意。
図 新規感染数の１週間増倍率$k$と実効再生産数$R_0$
2020年4月7日火曜日
緊急事態宣言と接触制限モデル（１）
北海道大学の西浦博さんが4月3日にマスコミで示していた図がある。対策がない状態の新規感染者数が指数関数的な増大をしているときに，人の接触を20%減らした場合と80%減らした場合の新規感染者数の変化を示した図だ。これにより彼は欧米に近い外出制限の必要性を訴えていた。残念ながら，これは専門家会議や政府の共通のコンセンサスとはならず，やがて4月7日の7都府県に対する5月6日までの緊急事態宣言につながっていく。
ここではその図がどのようにして得られたものかを，素人が持っている簡単な道具と知識で理解することを目指す。報道されたニュースや伝わってくる情報は鵜呑みにはせず，自分で考えることが必要だから。集団としては人口1400万人の東京都を想定する。日本経済新聞におけるこの内容に関する報道の文脈でも東京を対象としていることがうかがえる。東京の4月3日における新規感染者数は14人（新規感染数累計は753人）であり，人口の1ppm条件がちょうど達成された頃なので，これ以降をSIRモデルなどの単純な微分方程式系で扱うことは可能だろう。
今，探してもみあたらないのだが，西浦さんはこの報道の後に解説のための動画を公開していた。それによればしばらく前のドイツなみの再生産数を仮定し，現時点から15日目以降に制限措置を行った場合をシミュレートしているとした。なお，彼のモデルでは，計算開始から15日目が現在であり，30日目が制限措置の開始時点となる。
今から3週間前のドイツでは新規感染数累計の7日増倍率が7倍になっていた。仮に6.3倍（1日に1.3倍）を仮定すると，基本再生産数は$R_0 = 1 + 5/7 * \log 6.3 \approx 2.3$となる*。一方現在の東京では新規感染数累計の7日増倍率が2.5倍（1日に1.14倍）程度なので，基本再生産数は$R_0 = 1 + 5/7 * \log 2.5 \approx 1.65$（倍加時間と基本再生産率を参照）である。だからこのドイツ並の仮定はちょっとどうかと思うが，いちおうこれを採用する。
（*西浦さんはドイツなどの欧州の平均基本再生産数が$R_0=2.5$であると話していた）
つまり，現在100人の新規感染者数はこのまま放置すると15日後のt=30（1.3^15≒50）に5000人以上に達するという状況を西浦さんは想定していると思われる。今の東京の水準をそのままあてはめた場合は，15日後のt=30（1.14^15≒7）に東京の新規感染者数は700人程度であることに注意しよう。
SIRモデルに接触制限措置の効果を含めた微分方程式系をMathematicaで書くと次のようになる。nは集団人口（1400，以下人数の単位は1万人当りとなる），βは感染率（0.44），αは感染期間（5），νは感染数の初期値（0.0002），pは制限措置の程度（p=10→制限なし，p=9→80%に制限，p=6→20%に制限）である。なお，接触制限措置は，βにかけるシグモイド関数（幅が1日のオーダー）によってモデル化して階段関数は用いない。t=15が4月3日現在（新規感染者数=100）に対応する。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
f[n_,β_,α_,ν_,p_]:=
NDSolve[{S'[t]==-(p+(10-p)*Tanh[2*(30-t)])/10*β/n*J[t]*S[t],
J'[t]==(p+(10-p)*Tanh[2*(30-t)])/10*β/n*J[t]*S[t]-J[t]/α, R'[t]==J[t]/α, S[0]==n,J[0]==ν,R[0]==0}, {S,J,R}, {t,0,100}]
β=0.44 sol = f[1400, β, 7, 0.0002, 10];
s[t_] := S[t] /. sol[[1, 1]];
i[t_] := J[t] /. sol[[1, 2]];
r[t_] := R[t] /. sol[[1, 3]];
cf1[t_] = (10+0*Tanh[2*(30-t)])/10*β*i[t]*s[t]/1400;
sol = f[1400, β, 7, 0.0002, 9];
s[t_] := S[t] /. sol[[1, 1]];
i[t_] := J[t] /. sol[[1, 2]];
r[t_] := R[t] /. sol[[1, 3]];
cf2[t_] = (9+1*Tanh[2*(30-t)])/10*β*i[t]*s[t]/1400;
sol = f[1400, β, 7, 0.0002, 6];
s[t_] := S[t] /. sol[[1, 1]];
i[t_] := J[t] /. sol[[1, 2]];
r[t_] := R[t] /. sol[[1, 3]];
cf3[t_] = (6+4*Tanh[2*(30-t)])/10*β*i[t]*s[t]/1400;
Plot[{cf1[t], cf2[t], cf3[t]}, {t, 15, 45},
PlotRange -> {0, 2}, GridLines ->
{{29, 30, 31}, {0.1, 0.4, 0.5, 0.6, 0.7, 0.8}}]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
計算の結果，次のグラフが得られた。西浦さんのものと概ね一致する。
ここではその図がどのようにして得られたものかを，素人が持っている簡単な道具と知識で理解することを目指す。報道されたニュースや伝わってくる情報は鵜呑みにはせず，自分で考えることが必要だから。集団としては人口1400万人の東京都を想定する。日本経済新聞におけるこの内容に関する報道の文脈でも東京を対象としていることがうかがえる。東京の4月3日における新規感染者数は14人（新規感染数累計は753人）であり，人口の1ppm条件がちょうど達成された頃なので，これ以降をSIRモデルなどの単純な微分方程式系で扱うことは可能だろう。
今，探してもみあたらないのだが，西浦さんはこの報道の後に解説のための動画を公開していた。それによればしばらく前のドイツなみの再生産数を仮定し，現時点から15日目以降に制限措置を行った場合をシミュレートしているとした。なお，彼のモデルでは，計算開始から15日目が現在であり，30日目が制限措置の開始時点となる。
今から3週間前のドイツでは新規感染数累計の7日増倍率が7倍になっていた。仮に6.3倍（1日に1.3倍）を仮定すると，基本再生産数は$R_0 = 1 + 5/7 * \log 6.3 \approx 2.3$となる*。一方現在の東京では新規感染数累計の7日増倍率が2.5倍（1日に1.14倍）程度なので，基本再生産数は$R_0 = 1 + 5/7 * \log 2.5 \approx 1.65$（倍加時間と基本再生産率を参照）である。だからこのドイツ並の仮定はちょっとどうかと思うが，いちおうこれを採用する。
（*西浦さんはドイツなどの欧州の平均基本再生産数が$R_0=2.5$であると話していた）
つまり，現在100人の新規感染者数はこのまま放置すると15日後のt=30（1.3^15≒50）に5000人以上に達するという状況を西浦さんは想定していると思われる。今の東京の水準をそのままあてはめた場合は，15日後のt=30（1.14^15≒7）に東京の新規感染者数は700人程度であることに注意しよう。
SIRモデルに接触制限措置の効果を含めた微分方程式系をMathematicaで書くと次のようになる。nは集団人口（1400，以下人数の単位は1万人当りとなる），βは感染率（0.44），αは感染期間（5），νは感染数の初期値（0.0002），pは制限措置の程度（p=10→制限なし，p=9→80%に制限，p=6→20%に制限）である。なお，接触制限措置は，βにかけるシグモイド関数（幅が1日のオーダー）によってモデル化して階段関数は用いない。t=15が4月3日現在（新規感染者数=100）に対応する。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
f[n_,β_,α_,ν_,p_]:=
NDSolve[{S'[t]==-(p+(10-p)*Tanh[2*(30-t)])/10*β/n*J[t]*S[t],
J'[t]==(p+(10-p)*Tanh[2*(30-t)])/10*β/n*J[t]*S[t]-J[t]/α, R'[t]==J[t]/α, S[0]==n,J[0]==ν,R[0]==0}, {S,J,R}, {t,0,100}]
β=0.44 sol = f[1400, β, 7, 0.0002, 10];
s[t_] := S[t] /. sol[[1, 1]];
i[t_] := J[t] /. sol[[1, 2]];
r[t_] := R[t] /. sol[[1, 3]];
cf1[t_] = (10+0*Tanh[2*(30-t)])/10*β*i[t]*s[t]/1400;
sol = f[1400, β, 7, 0.0002, 9];
s[t_] := S[t] /. sol[[1, 1]];
i[t_] := J[t] /. sol[[1, 2]];
r[t_] := R[t] /. sol[[1, 3]];
cf2[t_] = (9+1*Tanh[2*(30-t)])/10*β*i[t]*s[t]/1400;
sol = f[1400, β, 7, 0.0002, 6];
s[t_] := S[t] /. sol[[1, 1]];
i[t_] := J[t] /. sol[[1, 2]];
r[t_] := R[t] /. sol[[1, 3]];
cf3[t_] = (6+4*Tanh[2*(30-t)])/10*β*i[t]*s[t]/1400;
Plot[{cf1[t], cf2[t], cf3[t]}, {t, 15, 45},
PlotRange -> {0, 2}, GridLines ->
{{29, 30, 31}, {0.1, 0.4, 0.5, 0.6, 0.7, 0.8}}]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
計算の結果，次のグラフが得られた。西浦さんのものと概ね一致する。
図１ SIRモデル＋シグモイド関数制限措置による新規感染者数の変化
（縦軸の単位は万人，横軸の単位は日＝3/20を基準とした経過日数,
青：制限なし，オレンジ：80%に制限，緑：20%に制限）
図２ 上記に，60%，40%，30%に制限などの場合を加えたもの
青：制限なし，オレンジ：80%に制限，緑：20%に制限）
図２ 上記に，60%，40%，30%に制限などの場合を加えたもの
なるほど，中途半端な制限ではだめだということか。30%に制限するあたりが増加かどうかの臨界点かもしれない。この度の緊急事態宣言はどの程度の効果があるのだろうか。
2020年4月6日月曜日
奈良コンベンションセンター
４月１日，近鉄奈良線新大宮駅から徒歩10分の大宮通りと三条通りに挟まれた区域に奈良コンベンションセンターがオープンした。近くのホテルマリオットはまだオープンしていなかったが，おしゃれな蔦屋書店は中川政七商店とコラボレーションしながら開店していた。新型コロナウイルス感染症で不要不急の外出を自粛させられている中，用事があったのと，まあ人はいないだろうという予測のもとに探索に出かけた。実際お客さんはほとんどいなかった。地下の駐車場もガラガラだった。
蔦屋書店がおしゃれであることは認めるけれど，本屋の価値の本質はおしゃれにあるわけではないので，あまり好きではない。見掛け倒しのために洋書を壁面に糊付けするセンスには耐えられないのだった（それはここじゃないけど）。もう一つの大問題。理工書がほとんどないのである。しかも数学と物理は妙にかけ離れた場所に存在していた。なんだかなあ。
もちろんビレッジバンガード（こちらはきらいじゃない）とか丸善などの本のアレンジだって好き嫌いはあるのだろうけれど，ジュンク堂にがっつり並んだ理工書のボリュームを見ると安心するのであった。
そんなわけで，文句をいいながらも，無人支払機を試してみるべく久しぶりに文庫本を一冊買って帰った。親切に使い方を教えていただいてありがとうございました。
蔦屋書店がおしゃれであることは認めるけれど，本屋の価値の本質はおしゃれにあるわけではないので，あまり好きではない。見掛け倒しのために洋書を壁面に糊付けするセンスには耐えられないのだった（それはここじゃないけど）。もう一つの大問題。理工書がほとんどないのである。しかも数学と物理は妙にかけ離れた場所に存在していた。なんだかなあ。
もちろんビレッジバンガード（こちらはきらいじゃない）とか丸善などの本のアレンジだって好き嫌いはあるのだろうけれど，ジュンク堂にがっつり並んだ理工書のボリュームを見ると安心するのであった。
そんなわけで，文句をいいながらも，無人支払機を試してみるべく久しぶりに文庫本を一冊買って帰った。親切に使い方を教えていただいてありがとうございました。
2020年4月5日日曜日
ニューヨーク州の集団免疫率
米国の集団免疫率からの続き
前回に続いて，IHME（Institute for Health Metrics and Ecaluation）におけるニューヨーク州（人口1900万人）の新型コロナウイルス感染症についての予測について考える。4月4日のWHOのデータによれば，米国の新規感染数累計は24万人，死亡数累計は5800人である。ニューヨーク州ではそれぞれが10万人と2900人であることから，米国全体の40%と50%をしめている。
IHMEのニューヨーク州の予測を，我々のSIIDR2モデルで追試してみる。ニューヨーク州のデータは上記の観察から，前回求めた米国のデータyaを2.5で，zaを2.0で割ることで求める（アバウトすぎるがオーダーがわかればよいという立場なので）。時間データxaは前回と同じ基準日（3月10日）である。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
xa=[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]
ya=[0.014,0.021,0.030,0.038,0.051,0.051,0.051,0.106,0.107,0.215,
0.317,0.462,0.462,0.958,1.28,1.58,1.93,2.07,2.59,3.14,
3.72,4.27,4.95,5.68,6.48]/2.5
za=[0.06,0.08,0.09,0.11,0.12,0.12,0.12,0.18,0.18,0.30,
0.46,0.61,0.61,1.22,1.43,2.04,2.68,3.01,3.77,5.06,
6.41,7.28,8.65,11.7,14.6]/100/2.0
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
SIIDR2モデル計算の下図の結果は次のパラメタで与えられる。
$\beta = 0.89, \nu = 0.01, \lambda = 21, \tau = 14, \alpha_1 = 5.0/0.80 ,\alpha_2 = 5.0/0.20, \gamma_1 = 15/0.92, \gamma_2 = 15/0.08$
$\gamma_2$は米国の場合と比べて2/3程度にとっていることに注意する。
前回に続いて，IHME（Institute for Health Metrics and Ecaluation）におけるニューヨーク州（人口1900万人）の新型コロナウイルス感染症についての予測について考える。4月4日のWHOのデータによれば，米国の新規感染数累計は24万人，死亡数累計は5800人である。ニューヨーク州ではそれぞれが10万人と2900人であることから，米国全体の40%と50%をしめている。
IHMEのニューヨーク州の予測を，我々のSIIDR2モデルで追試してみる。ニューヨーク州のデータは上記の観察から，前回求めた米国のデータyaを2.5で，zaを2.0で割ることで求める（アバウトすぎるがオーダーがわかればよいという立場なので）。時間データxaは前回と同じ基準日（3月10日）である。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
xa=[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]
ya=[0.014,0.021,0.030,0.038,0.051,0.051,0.051,0.106,0.107,0.215,
0.317,0.462,0.462,0.958,1.28,1.58,1.93,2.07,2.59,3.14,
3.72,4.27,4.95,5.68,6.48]/2.5
za=[0.06,0.08,0.09,0.11,0.12,0.12,0.12,0.18,0.18,0.30,
0.46,0.61,0.61,1.22,1.43,2.04,2.68,3.01,3.77,5.06,
6.41,7.28,8.65,11.7,14.6]/100/2.0
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
SIIDR2モデル計算の下図の結果は次のパラメタで与えられる。
$\beta = 0.89, \nu = 0.01, \lambda = 21, \tau = 14, \alpha_1 = 5.0/0.80 ,\alpha_2 = 5.0/0.20, \gamma_1 = 15/0.92, \gamma_2 = 15/0.08$
$\gamma_2$は米国の場合と比べて2/3程度にとっていることに注意する。
図１ ニューヨーク州の感染カーブ（u3=重症感染数，u4=死亡数，u6=新規感染数累計）
図２ ニューヨーク州の感染カーブ（同上，u5=回復（免疫獲得）数）
IHME予測の定性的な振る舞いを再現することができたりできなかったりしている。
① 4月10日の中旬（t=30）に重症感染数はピークアウトする（図１u3のピーク位置）
は10日ほど後にずれている。
② ピーク時の必要病床数（図１ u3のピーク値）が7.5万床→6.5万床程度。
③ 6月上旬（t=90）段階の死亡数は1.6万人に達する（図１のu4の収束値）はほぼ再現。
④ 6月上旬（t=90）にはおおむね終息している（ただしu3はテイルが長い）。
⑤ 終息後のニューヨーク州の集団免疫率は4%にになる（図２のu6の値400/1万=4%）
主観的意見：たぶん報道の様子からすると，米国の様々な感染症対策はこのIHMTのデータに基づいて立案されているのではないかと推察される。しかし，日本にはこれに対応するシミュレーションが存在しない，若しくは存在していてもオープンにされていない，若しくは存在しているがインプットデータの信頼性が著しく低いので正しく利用されていない（印象操作には用いられている）。
［１］基本再生産数と集団免疫率（2020.3.27）
［２］倍加時間と基本再生産数（2020.3.29）
［３］新規感染数累計の増倍率（2020.3.31）
［４］湖北省の集団免疫率（2020.3.28）
［５］韓国の集団免疫率（2020.4.3）
［６］感染症の数理シミュレーション（８）（2020.3.15）
2020年4月4日土曜日
米国の集団免疫率（１）
韓国の集団免疫率からの続き
ワシントン大学にある研究所のIHME（Institute for Health Metrics and Ecaluation）では，COVID-19 Projections として，米国の新型コロナウイルス感染症についての予測を行っている。入院病床，ICU病床，人工呼吸器などの必要数がシミュレーションされている。それによればこれらは4月の中旬にピークアウトし，ピーク時の必要病床数波26万床（必要ICU病床数は4万床）である。また，7月の初旬には終息するとされている。また7月段階での死亡数は9万人以上に達する。
これは，3月27日付の論文 "Forecasting COVID-19 impact on hospital bed-days, ICU-days, ventilator- days and deaths by US state in the next 4 months" によるものであり，Githubには彼らのGeneric curve fitting package with nonlinear mixed effects model の説明とそのコードもある。
これを我々のSIIDR2モデルで追試してみる。まず米国のデータをWHOのSituation Reportsからとって人口で正規化すると1万人当りの新規感染数累計ykと死亡数累計zkが得られる。なお，新規感染数累計が人口の1ppmを超えた基準日は3/10であり，3/10から4/3までの25日分のデータを与える。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
xa=[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]
ya=[0.014,0.021,0.030,0.038,0.051,0.051,0.051,0.106,0.107,0.215,
0.317,0.462,0.462,0.958,1.28,1.58,1.93,2.07,2.59,3.14,
3.72,4.27,4.95,5.68,6.48]
za=[0.06,0.08,0.09,0.11,0.12,0.12,0.12,0.18,0.18,0.30,
0.46,0.61,0.61,1.22,1.43,2.04,2.68,3.01,3.77,5.06,
6.41,7.28,8.65,11.7,14.6]/100
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
SIIDR2モデル計算の下図の結果は次のパラメタで与えられる。
$\beta = 0.72, \nu = 0.03, \lambda = 21, \tau = 14, \alpha_1 = 5.0/0.80 ,\alpha_2 = 5.0/0.20, \gamma_1 = 15/0.87, \gamma_2 = 15/0.13$
これまでの韓国や湖北省の計算では，$\gamma_1 = 15/0.96, \gamma_2 = 15/0.04$に固定していたが，それでは死亡数が再現できなかったので約3倍の値に設定していることに注意する。これによって1万人当たり（米国の人口は3.27億人なので下記を3270倍すると実人数になる）の値が得られる。
ワシントン大学にある研究所のIHME（Institute for Health Metrics and Ecaluation）では，COVID-19 Projections として，米国の新型コロナウイルス感染症についての予測を行っている。入院病床，ICU病床，人工呼吸器などの必要数がシミュレーションされている。それによればこれらは4月の中旬にピークアウトし，ピーク時の必要病床数波26万床（必要ICU病床数は4万床）である。また，7月の初旬には終息するとされている。また7月段階での死亡数は9万人以上に達する。
これは，3月27日付の論文 "Forecasting COVID-19 impact on hospital bed-days, ICU-days, ventilator- days and deaths by US state in the next 4 months" によるものであり，Githubには彼らのGeneric curve fitting package with nonlinear mixed effects model の説明とそのコードもある。
これを我々のSIIDR2モデルで追試してみる。まず米国のデータをWHOのSituation Reportsからとって人口で正規化すると1万人当りの新規感染数累計ykと死亡数累計zkが得られる。なお，新規感染数累計が人口の1ppmを超えた基準日は3/10であり，3/10から4/3までの25日分のデータを与える。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
xa=[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]
ya=[0.014,0.021,0.030,0.038,0.051,0.051,0.051,0.106,0.107,0.215,
0.317,0.462,0.462,0.958,1.28,1.58,1.93,2.07,2.59,3.14,
3.72,4.27,4.95,5.68,6.48]
za=[0.06,0.08,0.09,0.11,0.12,0.12,0.12,0.18,0.18,0.30,
0.46,0.61,0.61,1.22,1.43,2.04,2.68,3.01,3.77,5.06,
6.41,7.28,8.65,11.7,14.6]/100
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
SIIDR2モデル計算の下図の結果は次のパラメタで与えられる。
$\beta = 0.72, \nu = 0.03, \lambda = 21, \tau = 14, \alpha_1 = 5.0/0.80 ,\alpha_2 = 5.0/0.20, \gamma_1 = 15/0.87, \gamma_2 = 15/0.13$
これまでの韓国や湖北省の計算では，$\gamma_1 = 15/0.96, \gamma_2 = 15/0.04$に固定していたが，それでは死亡数が再現できなかったので約3倍の値に設定していることに注意する。これによって1万人当たり（米国の人口は3.27億人なので下記を3270倍すると実人数になる）の値が得られる。
図１ 米国の感染カーブ（u3=重症感染数，u4=死亡数，u6=新規感染数累計）
図２ 米国の感染カーブ（同上，u5=回復（免疫獲得）数）
我々のモデルでもIHME予測の定性的な振る舞いを再現することができる。
① 4月の中旬（t=35）に重症感染数はピークアウトする（図１u3のピーク位置）。
② ピーク時の必要病床数は26万床に達する（図１ u3のピーク値）。
③ 7月上旬（t=110）段階の死亡数は9万人に達する（図１のu4の収束値）。
④ 7月上旬（t=110）には終息している。
⑤ 終息後の米国の集団免疫率は1%にすぎない（図２のu6の値100/1万=1%）
［１］基本再生産数と集団免疫率（2020.3.27）
［２］倍加時間と基本再生産数（2020.3.29）
［３］新規感染数累計の増倍率（2020.3.31）
［４］湖北省の集団免疫率（2020.3.28）
［５］韓国の集団免疫率（2020.4.3）
［６］感染症の数理シミュレーション（８）（2020.3.15）
2020年4月3日金曜日
韓国の集団免疫率
湖北省の集団免疫率からの続き
中国以外の多くの国では感染拡大が終息していないが，韓国では中国に続いて感染拡大が収まりつつあるように見える。そこで，韓国では現在どの程度の集団免疫が獲得されたかを，前回と同様にSIIDR2モデルで推定してみよう。
中国以外の多くの国では感染拡大が終息していないが，韓国では中国に続いて感染拡大が収まりつつあるように見える。そこで，韓国では現在どの程度の集団免疫が獲得されたかを，前回と同様にSIIDR2モデルで推定してみよう。
まず，データはいつものようにWHOのSituation Report（South Korea）を用いる。新規感染数累計（Confirmed）と死亡数累計（Deaths）の日次統計を韓国の人口5180万人でわったものから，1万人当りの人口比データを取り出した。新規感染数累計が人口の1ppm（百万分の一）を越えた時点を基準日とすると，韓国の基準日は2020年2月19日となる。このデータを示しておく。xkが基準日からの日数（2/19〜4/2までの44日間），ykがConfirmedの人口比（/万人），zkがDeathsの人口比（/万人（zk全体を100で割っていることに注意する）。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
xk=[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]
yk=[0.010,0.020,0.039,0.067,0.116,0.147,0.189,0.243,0.341,0.451,
0.608,0.721,0.813,0.929,1.029,1.113,1.213,1.306,1.377,1.425,
1.450,1.497,1.519,1.540,1.561,1.576,1.590,1.606,1.606,1.624,
1.670,1.699,1.718,1.730,1.745,1.764,1.784,1.802,1.830,1.850,
1.865,1.889,1.909,1.926]
zk=[0.00,0.02,0.02,0.04,0.10,0.14,0.19,0.23,0.25,0.25,
0.33,0.35,0.42,0.54,0.62,0.68,0.81,0.85,0.97,0.98,
1.04,1.16,1.27,1.27,1.39,1.45,1.45,1.56,1.56,1.62,
1.81,1.97,2.01,2.14,2.32,2.43,2.53,2.68,2.78,2.93,
3.05,3.13,3.19,3.26]/100
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
xk=[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]
yk=[0.010,0.020,0.039,0.067,0.116,0.147,0.189,0.243,0.341,0.451,
0.608,0.721,0.813,0.929,1.029,1.113,1.213,1.306,1.377,1.425,
1.450,1.497,1.519,1.540,1.561,1.576,1.590,1.606,1.606,1.624,
1.670,1.699,1.718,1.730,1.745,1.764,1.784,1.802,1.830,1.850,
1.865,1.889,1.909,1.926]
zk=[0.00,0.02,0.02,0.04,0.10,0.14,0.19,0.23,0.25,0.25,
0.33,0.35,0.42,0.54,0.62,0.68,0.81,0.85,0.97,0.98,
1.04,1.16,1.27,1.27,1.39,1.45,1.45,1.56,1.56,1.62,
1.81,1.97,2.01,2.14,2.32,2.43,2.53,2.68,2.78,2.93,
3.05,3.13,3.19,3.26]/100
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
SIIDR2モデルのパラメタを次のようにとると韓国の上記データがおおむね再現できる。
$\beta=0.95, \nu=0.03, \lambda=8, \tau=5, \alpha_1=5.0/0.80, \alpha_2=5.0/0.20, \gamma_1=15.0/0.96, \gamma_2=15.0/0.04$
なお，下記の図では，u3が重症感染者，u4が死亡数累計（Deaths），u5が回復(免疫獲得)者，u6が新規感染数累計（Confirmed）を与えている。
図１ 韓国の感染カーブ
もちろん他のパラメタセットの範囲でも同じ結果を与えるが，その場合でも最終的な目標値の集団免疫率の推定値は同じオーダーに収まっている。集団免疫率は集団全体に対する回復（免疫獲得）者の割合なので，現時点での図からわかるように，韓国における集団免疫率は10人/1万人〜0.1%のオーダーに過ぎない。
なお，最近（25〜40日）のu6は落ち着いているとはいうものの完全な収束ではなく，直線状に増加しているようにみえるので，我々のモデルではこの定性的な振る舞いを説明できていない。これは毎日，人口の0.2%程度にあたる一定数の感染者が引き続き発生していることを意味している。
参考のため，このモデルの結果を与えている実効再生産数（パラメタセットで定まる）$R_{\rm eff} = \alpha*\beta(t)$を以下に示しておく。
図２ 韓国の実効再生産数カーブ
2020年4月2日木曜日
「オーバーシュート」の見分け方
新型コロナウイルス感染症対策本部は2020年1月30日に閣議決定により内閣官房に設置された。内閣総理大臣を本部長とし，官房長官や他の国務大臣が構成員である。その後3月26日に改めて，「新型インフルエンザ等対策特別措置法(平成24年法律第31号) 第15条第1項の規定に基づいて新型コロナウイルス感染症対策本部を設置した」とされた。
2月16日には，この対策本部のもとに，国立感染症研究所の脇田隆字所長を座長として医学・医療関係者からなる新型コロナウイルス感染症対策専門家会議が設置され第１回の会合が開かれている。なお，テレビでよくみかける感染症数理が専門の北海道大学の西浦博教授は座長が出席を求める関係者だった。
3月19日開催の第8回会議における「新型コロナウイルス感染症の対策の状況分析・提言」において，第1回から第7回までまったく登場していなかった「オーバーシュート」というキーワードが突如として次の文脈で デビューした（強調は筆者）。
○オーバーシュートは指数関数的増大でなく過渡現象の振れ過ぎ部分を指すものだ。
○爆発的患者急増はアウトブレイクというのではなかったのか（図１）。
○手元の辞書にはそんな意味が載っていない。
○何故一般向けにわかりにくい英語の専門用語を使うのか。
○そもそもこれは学術専門用語ではない（日本環境感染学会用語集に存在しない）。
○欧米の感染症情報にもほとんど出てこないしあってもまったく意味が違う。
○指数関数的増大のどこからがオーバーシュートなのかわからない。
それにもかかわらず，マスコミや政治家によってオーバーシュートは我が国に定着することになった。さすがにちょっと気が引けたのか，専門家会議は4月1日の専門家会議の後の会見で，オーバーシュートの定義を次のように説明することになった。
前提として，SIR的なモデルを仮定し感染期間（感染者数の崩壊定数）を$\gamma$=5日とする。累積患者数は 新規感染数累計，$J(t)=\int_{0}^{t} I(t') dt'$のことであるとして，感染者数が非常に少ない状態から指数関数的に増大する場面を設定する。すでに基本的な関係式は下記参考資料で導出しているので，ここでは結果のみ整理する。時刻$t$日の$J(t)$に対する$\tau$日後の$J(t+\tau)$の比率を$r$とおく。なお，$R_0$は基本再生産数（あるいはその時点での実効再生産数）とする（図２）。
\begin{equation}
\begin{aligned}
r &= \dfrac{J(t+\tau)}{J(t)} = \dfrac{e^{\frac{R_0-1}{\gamma}(t+\tau)}-1}{e^{\frac{R_0-1}{\gamma} t}-1} = \dfrac{e^{\frac{R_0-1}{\gamma}\tau}-e^{-\frac{R_0-1}{\gamma} t}}{1 - e^{- \frac{R_0-1}{\gamma} t}} \approx e^{\frac{R_0-1}{\gamma}\tau}\\
R_0 &= 1 + \dfrac{\gamma}{\tau} \log r
\end{aligned}
\end{equation}
① 3日で2倍すなわち$\tau=3, r=2$とすると，$R_0$=2.15となる。
② $R_0$=2.15の場合，7日で何倍かというと，$r=e^{7 * 1.15 /5}$=5.00倍になる。
③ 東京では，7日で3倍になっているので，$R_0$=1.78となる。
つまり，日次報告されている新規感染数累計（Confirmed）$J(t)$が 1週間で5倍程度になる状況が続けば，いわゆる「オーバーシュート」の状態であるということになる。東京ではこれが3倍程度でありそこまでは至っていない。欧州や以前の韓国や湖北省ではオーバーシュートの状態が生じていたが，現在は収まっており$r$も減少している。北米ではオーバーシュート近傍にあるが，$r$は減少傾向にある。日本や東京ではオーバーシュート状態ではないが，$r$は増加傾向にあって予断を許さない（下記の各国の感染者比で$r$を確認）。
$r$は比率なので，報告されている感染数の把握度が不十分であっても（見逃している感染数がかなり存在しているかもしれないとしても），把握度の時間依存性が大きくない限りこの結論には影響しない。
［１］新規感染数累計の増倍率（2020.3.31）
［２］倍加時間と基本再生産数（2020.3.29）
［３］基本再生産数と集団免疫率（2020.3.27）
［４］感染症の数理シミュレーション（８）（2020.3.15）
［５］各国の感染者比（２）（2020.3.12）
［６］各国の感染者比（１）（2020.3.10）
［６］感染症の数理シミュレーション（２）（2020.2.25）
［７］感染症の数理シミュレーション（１）（2020.2.24）
2月16日には，この対策本部のもとに，国立感染症研究所の脇田隆字所長を座長として医学・医療関係者からなる新型コロナウイルス感染症対策専門家会議が設置され第１回の会合が開かれている。なお，テレビでよくみかける感染症数理が専門の北海道大学の西浦博教授は座長が出席を求める関係者だった。
3月19日開催の第8回会議における「新型コロナウイルス感染症の対策の状況分析・提言」において，第1回から第7回までまったく登場していなかった「オーバーシュート」というキーワードが突如として次の文脈で デビューした（強調は筆者）。
特に，気付かないうちに感染が市中 に拡がり，あるときに突然爆発的に患者が急増(オーバーシュート(爆発的患者急増))す ると，医療提供体制に過剰な負荷がかかり，それまで行われていた適切な医療が提供でき なくなることが懸念されます。SNSでは理系の人々を中心に違和感ありまくりでブーイングが起こっている。
○オーバーシュートは指数関数的増大でなく過渡現象の振れ過ぎ部分を指すものだ。
○爆発的患者急増はアウトブレイクというのではなかったのか（図１）。
○手元の辞書にはそんな意味が載っていない。
○何故一般向けにわかりにくい英語の専門用語を使うのか。
○そもそもこれは学術専門用語ではない（日本環境感染学会用語集に存在しない）。
○欧米の感染症情報にもほとんど出てこないしあってもまったく意味が違う。
○指数関数的増大のどこからがオーバーシュートなのかわからない。
図１ outbreak(赤) と overshoot(青) の google serarch trend (30日)
上は世界全体，下は日本
副座長の尾身茂氏（地域医療機能推進機構・理事長）は「オーバーシュートは欧米で見られるように爆発的な患者数の増加を示すが，2日ないし3日のうちに累積患者数が倍増し，しかもそのスピードが継続的にみられる状態を指すと私たちは定義した」と述べた。実は，「医療提供体制に過剰な負荷がかかり，それまで行われていた適切な医療が提供できなくなる」ことが彼らの強調したいことであって，そのためにあえて新しい用語の選択によって注意喚起したのではないかとも推察できるが，これをグダグダいっても仕方がないので， 2から3日のうちに累積患者数が倍増しという定義を受け入れた場合のオーバーシュートの見分け方について考える。
前提として，SIR的なモデルを仮定し感染期間（感染者数の崩壊定数）を$\gamma$=5日とする。累積患者数は 新規感染数累計，$J(t)=\int_{0}^{t} I(t') dt'$のことであるとして，感染者数が非常に少ない状態から指数関数的に増大する場面を設定する。すでに基本的な関係式は下記参考資料で導出しているので，ここでは結果のみ整理する。時刻$t$日の$J(t)$に対する$\tau$日後の$J(t+\tau)$の比率を$r$とおく。なお，$R_0$は基本再生産数（あるいはその時点での実効再生産数）とする（図２）。
\begin{equation}
\begin{aligned}
r &= \dfrac{J(t+\tau)}{J(t)} = \dfrac{e^{\frac{R_0-1}{\gamma}(t+\tau)}-1}{e^{\frac{R_0-1}{\gamma} t}-1} = \dfrac{e^{\frac{R_0-1}{\gamma}\tau}-e^{-\frac{R_0-1}{\gamma} t}}{1 - e^{- \frac{R_0-1}{\gamma} t}} \approx e^{\frac{R_0-1}{\gamma}\tau}\\
R_0 &= 1 + \dfrac{\gamma}{\tau} \log r
\end{aligned}
\end{equation}
① 3日で2倍すなわち$\tau=3, r=2$とすると，$R_0$=2.15となる。
② $R_0$=2.15の場合，7日で何倍かというと，$r=e^{7 * 1.15 /5}$=5.00倍になる。
③ 東京では，7日で3倍になっているので，$R_0$=1.78となる。
図２ 1週間での$J(t)$の増倍率$r$と基本再生産数$R_0$
つまり，日次報告されている新規感染数累計（Confirmed）$J(t)$が 1週間で5倍程度になる状況が続けば，いわゆる「オーバーシュート」の状態であるということになる。東京ではこれが3倍程度でありそこまでは至っていない。欧州や以前の韓国や湖北省ではオーバーシュートの状態が生じていたが，現在は収まっており$r$も減少している。北米ではオーバーシュート近傍にあるが，$r$は減少傾向にある。日本や東京ではオーバーシュート状態ではないが，$r$は増加傾向にあって予断を許さない（下記の各国の感染者比で$r$を確認）。
$r$は比率なので，報告されている感染数の把握度が不十分であっても（見逃している感染数がかなり存在しているかもしれないとしても），把握度の時間依存性が大きくない限りこの結論には影響しない。
［１］新規感染数累計の増倍率（2020.3.31）
［２］倍加時間と基本再生産数（2020.3.29）
［３］基本再生産数と集団免疫率（2020.3.27）
［４］感染症の数理シミュレーション（８）（2020.3.15）
［５］各国の感染者比（２）（2020.3.12）
［６］各国の感染者比（１）（2020.3.10）
［６］感染症の数理シミュレーション（２）（2020.2.25）
［７］感染症の数理シミュレーション（１）（2020.2.24）
2020年4月1日水曜日
インターネット授業への道
昨日，非常勤先の大阪教育大学の方針が変更された。当初，予定されていた学年暦どおりに4月8日から授業を開始するとのことであったのだが，やはりこの状況なので，通常授業は困難であると判断され，4月8日から4月19日までは休講となった。一安心。
4月20日からインターネット授業（インターネットを活用した授業），5月11日からは対面授業を含めて実施ということである。昨日はzoomによる緊急のFDが開催されて，moodleを活用したインターネット授業の典型例，実施のための体制はどうすればよいか，ということが議論された。
昨年度の大阪教育大学のmoodleの活用率は1/6程度らしいのでなかなか険しい道のりなのである。
学習院大学の田崎さんのGIGAにやさしい遠隔授業を参考にしようと思うが，あと3週間で準備できるだろうか。
4月20日からインターネット授業（インターネットを活用した授業），5月11日からは対面授業を含めて実施ということである。昨日はzoomによる緊急のFDが開催されて，moodleを活用したインターネット授業の典型例，実施のための体制はどうすればよいか，ということが議論された。
昨年度の大阪教育大学のmoodleの活用率は1/6程度らしいのでなかなか険しい道のりなのである。
学習院大学の田崎さんのGIGAにやさしい遠隔授業を参考にしようと思うが，あと3週間で準備できるだろうか。
2020年3月31日火曜日
新規感染数累計の増倍率
ある集団の単位時間（1日）当りの新規感染数を$f(t)$，新規感染数累計（Confirmed）を$g(t)=\int_{t_0}^t f(t') dt'$とする。$t_0$は基準日である。このとき，$r(t)=g(t+\tau)/g(t)$は期間$\tau$の間に新規感染数累計が何倍に増えたかを示す増倍率を表す。感染が終息すれば
$f(t)$は0となり，$g(t)$は一定に近づくので$r(t) \rightarrow 1$となる。
新規感染数が指数関数的に増大する場合は，$f(t)= a e^{t/\lambda}$となるので，$g(t) = a \lambda (e^{t / \lambda} - 1)$ となる。この場合，$r(t) = \dfrac{e^{(t+\tau)/ \lambda}-1}{e^{ t/\lambda}-1}$であることから，$\lim_{t \to \infty} r(t) = e^{\tau / \lambda}$ に漸近する。
そこで，WHOが毎日報告しているCOVID-19のSituation Reportsの各国の新規感染数累計（Confirmed）の日次統計から $\tau$=7日間としたときの$r(t)$を求めてその特徴を比較する。日次統計には揺らぎがあるので，$\bar{r}(t) = \frac{1}{5} \int_{t-2}^{t+2} r(t') dt'$という5日移動平均をExcelによって求めてグラフを描いた。なお基準日はWHOの統計が始まった2020年1月21日をt=1(日)とする。図１のアジア・太平洋地域の国々はt=41(3月1日)〜t=68(3月28日)を，図２の欧州・北米の国々はt=47(3月7日)〜t=68(3月28日)の範囲で計算した。
① アジア・太平洋地域の直近で，オーストラリアの4倍〜韓国の1倍までに対して，欧州・北米地域でアメリカ合衆国の6倍〜イラン・イタリアの2倍の範囲にあり，それぞれ漸減中である。ただし，東京・日本は増加傾向にあり，カナダはっきりした減少傾向が見えない。
② 初期にインフレーションを起した韓国・イラン・イタリアはその後単純な減少カーブを描いている。それ以外の国で複数の山が観測される。例えば，当初収束しているといわれていた台湾，香港，シンガポールに第二の緩やかなピークが生じ，マレーシアやオーストラリア，スペインやアメリカにも大きな第二のピークが現れている。
$f(t)$は0となり，$g(t)$は一定に近づくので$r(t) \rightarrow 1$となる。
新規感染数が指数関数的に増大する場合は，$f(t)= a e^{t/\lambda}$となるので，$g(t) = a \lambda (e^{t / \lambda} - 1)$ となる。この場合，$r(t) = \dfrac{e^{(t+\tau)/ \lambda}-1}{e^{ t/\lambda}-1}$であることから，$\lim_{t \to \infty} r(t) = e^{\tau / \lambda}$ に漸近する。
そこで，WHOが毎日報告しているCOVID-19のSituation Reportsの各国の新規感染数累計（Confirmed）の日次統計から $\tau$=7日間としたときの$r(t)$を求めてその特徴を比較する。日次統計には揺らぎがあるので，$\bar{r}(t) = \frac{1}{5} \int_{t-2}^{t+2} r(t') dt'$という5日移動平均をExcelによって求めてグラフを描いた。なお基準日はWHOの統計が始まった2020年1月21日をt=1(日)とする。図１のアジア・太平洋地域の国々はt=41(3月1日)〜t=68(3月28日)を，図２の欧州・北米の国々はt=47(3月7日)〜t=68(3月28日)の範囲で計算した。
図１ アジア・太平洋地域の$\bar{r}(t)$
図２ 欧州・北米地域の$\bar{r}(t)$
① アジア・太平洋地域の直近で，オーストラリアの4倍〜韓国の1倍までに対して，欧州・北米地域でアメリカ合衆国の6倍〜イラン・イタリアの2倍の範囲にあり，それぞれ漸減中である。ただし，東京・日本は増加傾向にあり，カナダはっきりした減少傾向が見えない。
② 初期にインフレーションを起した韓国・イラン・イタリアはその後単純な減少カーブを描いている。それ以外の国で複数の山が観測される。例えば，当初収束しているといわれていた台湾，香港，シンガポールに第二の緩やかなピークが生じ，マレーシアやオーストラリア，スペインやアメリカにも大きな第二のピークが現れている。
2020年3月30日月曜日
深紅の碑文：上田早夕里
上田早夕里は，しっかりした世界観に基づいた骨太のSF作品を産み出し続けている。「深紅の碑文」はオーシャンクロニクルシリーズの「華竜の宮」に続く長編作品だ。ホットプルームによる海底隆起で陸地が水没し，人類が陸上民と遺伝子操作された海上民・魚舟に分離する世界で，さらなる環境激変（大異変）にどう対応するかが，複雑な社会関係と政治的駆け引きを書き込みながら語られる。
小松左京と比較されることが多いが，一定の科学的論理に裏付けられて大きな物語が構築される様子は似ているのだろう。ただ，小松左京には，その短編にも見られるような関西のお笑いのノリ（初期の漫才台本原稿の効果？）というのか，長編にも明るさが感じられるのだが，上田早夕里はどこかきまじめでその語りにも甘さより複雑な苦さや暗さが先にくる。日本酒の分類表みたいなものであり，どれも味わい深いのだけれど。最も上田早夕里の他のシリーズは読んだことがないのであくまでも個人の感想です。
ハヤカワJA文庫下巻の渡邊利道さんの解説によれば，長編パートはこれで終了であとは中短編による連作らしいが，海と地球だけでなく深宇宙の物語も早く読んでみたい。
小松左京と比較されることが多いが，一定の科学的論理に裏付けられて大きな物語が構築される様子は似ているのだろう。ただ，小松左京には，その短編にも見られるような関西のお笑いのノリ（初期の漫才台本原稿の効果？）というのか，長編にも明るさが感じられるのだが，上田早夕里はどこかきまじめでその語りにも甘さより複雑な苦さや暗さが先にくる。日本酒の分類表みたいなものであり，どれも味わい深いのだけれど。最も上田早夕里の他のシリーズは読んだことがないのであくまでも個人の感想です。
ハヤカワJA文庫下巻の渡邊利道さんの解説によれば，長編パートはこれで終了であとは中短編による連作らしいが，海と地球だけでなく深宇宙の物語も早く読んでみたい。
2020年3月29日日曜日
倍加時間と基本再生産数
マスコミを巧妙に活用した情宣作戦というのか，"やってる感パフォーマンス"の演出が大好きなのは維新と安倍だが，土曜日の夕方にまたまた"記者会見もどき"の宣伝工作を行っていた。内容はないようだったが，いつものように滑舌が悪くて"指示"が"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}
さて，あいかわらず醜悪なポエムに包まれた左右首振りプロンプター官僚作文の中で，一点だけ定量的な発言として興味を引いたのが，感染者数が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月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 は，ウエブのスレッドということでした。以前のみらい翻訳と比べてどうだろうか。一長一短かもしれない。
登録:
投稿 (Atom)