スーパーシミュレーター

今朝，起床の直前に見た夢の話。夢言語から日本語へ翻訳すると概ね次のようなことだった。もちろん，下記のテキストには生成AIは一切関与していない。

登場人物：自分，良く知らない学生さん

場所：学校の教室らしきところ

学生さんに，最近のAIのトレンドを聞いていたら，新しいアプリを紹介してくれた。今年の8月に入ってから登場したアプリらしい。そのデモ画面では，机上の大きなスクリーンに展開した，子ども向けの図鑑のようなアイコンやシンボルが配列されたものだ。それは世界のスーパーシミュレーターになっている。我々の世界の経済生産活動や政治社会活動や文化芸術活動の様々なデータが地図上の操作によって簡単に取り出すことができる。さらに，その場所における社会のパラメタを少し変化させると，これまでくっきりしていたアイコンがぼやけて全体のフォーカスが崩れるのだが，しばらくすると世界は新しい条件の下に実現する状態に遷移して世界全体のシミュレーションが続いていく。試してみたのは中国の東北部なのか日本の東北地方なのかははっきりしないが，あっという間に結果を得ることができた。自然現象のシミュレーションも可能であり，日本列島を西から東へと鳥の様々な種類が変化しながら覆っていく様子がみえた。石川県の辺りからは，ピンクの朱鷺だとおもわれる鳥のアイコンが拡がっていった（夢から覚めた後で日経朝刊をみると，26面に「能登でトキ放鳥へ 26年度にも 復興支援で本州初」という記事がみつかったけれど・・・）。この新しいアプリは一体なんというものなのか，生成AI（LLM）を背景で使っているのだろうと学生さんに尋ねたけれど，どうやらそうではないらしい。さらに，URLを教えてよとお願いすると，教えてくれたので，チョークで黒板に，ボールペンでノートに書き留めようとするが，自分の電卓が見当たらないのでうまくいかない。あちらこちら探し回るが電卓がみつからないのだ。それはとても困る。別の部屋に行くと関数電卓が何台か置いてあったが，どれも自分のものではない。仕方がないので，それらをまとめて持って元の部屋に帰ると，学生さんが，書類を持っていた。彼はどうやら自分の研究室に所属する予定の学生だったようだ。改めてこのスーパーシミュレーターアプリのURLを教えてもらおうとするが。自分の学生だから，まあ安心だ。いつかは必ず聞き出すことができるだろう。それにしても，AIの進化のスピードははやい。あっという間に，かなり精度が高くて有用な世界のシミュレータが皆の手元に存在する時代になってしまった。これからいったいどうなるのだろう。

P. S. そうそう，たぶん自分の大好きなダイアスパーがモデルですよね。

図：DALL-E3のこのイメージは自分のPOVとは異なっている

［１］『都市と星』完成への道程とクラークの作家性

［２］都市と星（Manuke Station SF Review）

［３］JackIn Space: 一人称・三人称映像間の連続的な遷移を可 能にするテレプレゼンスシステム

ブランコの物理

Physical Reveiw E に，Initial phase and frequency modulations of pumping a playground swing という日本のグループの論文が掲載された。ブランコを漕いで揺らすコツを解明したというものだ，たぶん。アブストラクトをDeepLにかけてGPT-4で整えると次のようになる。

遊具のブランコは、ブランコ本体とそれを揺らす人間からなる動的な結合振動子系である。本研究では、上半身の自然な動きの初期段階がブランコの連続的なポンピングに与える影響を分析するモデルを提案し、3種類の長さのブランコで10人の参加者が実際にポンピングした際の運動データを用いて検証する。

研究の結果から、スイングが垂直位置（中点）にあり、振幅が小さい時に前方へ移動する際、最大リーンバックの初期段階が発生し、最大振幅のポンピングが予測された。振幅が大きくなるにつれて、最適な初期段階はサイクルの初期段階、すなわちスイングの軌道の後端に向かって徐々に移動する。モデルの予測に従って、参加者全員が、スイングの振幅が大きくなるにつれて上半身の動作の初期位相が早くなることが確認された。

この結果から、スイングする人は、上半身の動作の頻度と初期位相を適切に調整することで、ブランコを効果的にポンピングすることができると考えられる。

人間の体を3つの線質量要素で表現し，上半身と下肢の角度をそれぞれ変数として導入したラグランジアンを考えている。このとき，これらの身体の角度については，何通りかのあらかじめ定められたモデルの周期運動だとしている。

現実に近いのだけれど面倒なモデルなので今一つピンと来ない。そこで，人体を振り子につながって距離が固定された２質点で表し，その間で周期的な質量移動をすることでスウィングするというモデルを考えてみた。質量移動は，位相と角振動数の２つの変数で制御できる。速度に比例する減衰振動のある単振子をこの質量移動で励振するというものだ。

Do[

 s0 = Pi/6*0; w0 = Pi; m1 = 1; m2 = 1;

 λ = Pi/12; ν = 0.05/Pi; ε = 0.95 + i*0.01; δ = Pi/8*0;

 m[t_] := (m1 + m2)/20 * Cos[ε w0 t + δ];

 sol = NDSolve[

       {(m1 + m2) s''[t] == -w0^2 * 

        ( (m1 + m[t]) * (Sin[s[t] + λ + ν * s'[t]) 

        + (m2 - m[t]) * (Sin[s[t] - λ + ν * s'[t])), 

        s[0] == s0, s'[0] == 0}, s, {t, 0, 60}];

  s /. sol[[1]];

  f[t_] := s[t] /. sol[[1]];

  g[t_] := s'[t] /. sol[[1]];

  h[i] = Plot[f[t], {t, 0, 30}, AspectRatio -> 1/2, 

   PlotRange -> {-0.6, 0.6}, PlotStyle -> Hue[0.08*i]],

 {i, 0, 7}]

もとの振動子の初期位相は s0 であり，固有角振動数は w0 (=√g/ℓ)，２質点の距離は角度 λ で表現している。速度に比例した減衰係数は ν，質量移動の角振動数と固有角振動数の比が ε，質量移動の位相のズレが δ である。以下の図は，ε = 0.95 + i*0.01 としたもの（上）と，δ = Pi/8*i; としたもの（下）である。固有角振動数比が 0.97のときが最も励振が大きく，2%ふえるだけで減衰に移行してしまう。

図：振幅の質量移動振動への依存性，振動数（上）と位相（下）

弾道ミサイルの軌道（６）

弾道ミサイルの軌道（４）からの続き

H3ロケットの発射がフェイルセーフで中止された次の日，なんだか間の悪いことに北朝鮮のICBM級ミサイルが平壌近傍から発射され，北海道渡島大島西200kmに落下した。

飛行時間は66分，飛行距離は900km，最高高度は5700km，落下速度はマッハ20程度と推測される。ペイロードを調整すれば，射程が14,000kmでアメリカ本土が射程に入ると報道された。

それでは，以前作ったのプログラムで確認してみよう。一生懸命Juliaのフォルダを探していたが，Mathematicaのコードだった。昨年11月のコードのパラメータで，飛行時間と飛行距離を微調整した結果が次の通り。

g = 0.0098; R = 6350; τ= 86; p = 0.75; 

a = 0.0446; s = 86.5 Degree; T = 3960; 

(* g = 0.0098; R = 6350; τ = 87; p = 0.75; a

 = 0.0446; s = 45 Degree; T = 5100; *)

fr[t_, τ_] := a*Sin[s]*HeavisideTheta[τ - t]

ft[t_, τ_] := a*Cos[s]*r[t]*HeavisideTheta[τ - t]

fm[t_, τ_] := -p/(τ - p*t)*HeavisideTheta[τ - t]

sol = NDSolve[{r''[t] == -fm[t, τ]*r'[t] + h[t]^2/r[t]^3 - 

 g R^2/r[t]^2 + fr[t, τ], r[0] == R, r'[0] == 0, 

 h'[t] == -fm[t, τ]*h[t] + ft[t, τ], h[0] == 0}, {r, h},

 {t, 0, T}]

f[t_] := r[t] /. sol[[1, 1]]

d[t_] := h[t] /. sol[[1, 2]]

Plot[{6350, f[t]}, {t, 0, T}]

Plot[{f[t + 1] - f[t], d[t]*R/f[t]^2, d[t]/f[t]}, 

{t, 0, T}, PlotRange -> {-4, 8}]

tyx[T_] := {T, f[T] - R, NIntegrate[R d[t]/f[t]^2 , {t, 0, T}]}

v[T_] := Sqrt[(f[T] - f[T - 1])^2 + (R d[T]/f[T]^2)^2]

{tyx[T], v[T]/0.340}

g0 = ParametricPlot[{NIntegrate[R d[t]/f[t]^2 , {t, 0, tt}], 

  f[tt] - R}, {tt, 0, T}]

燃焼時間τが86秒，燃料重量比が0.75，燃焼加速度が 0.0446km/s^2，投射角が86.5度である。これで飛行時間を与えると，飛行距離と最高高度と落下時速度が概ね再現できる。万有引力は距離の2乗に反比例し，コリオリ力や空気抵抗は無視，地球の形を考慮して2次元極座標で質量が変化する1段ロケットの微分方程式を解いている。

燃焼時間τを87秒にして，投射角を45度にすれば，飛行時間が5100秒で飛行距離は14,600km，落下時速度はマッハ24になる。このときペイロードはほとんど変化させていない。燃料重量比は同じで燃焼時間を1秒のばしただけだ。

図：投射角86度のロフテッド軌道と投射角45度の弾道軌道

SIDRモデル（２）

 SIDRモデル（１）からの続き

牧野さんの例では，感染期間/免疫消失期間が1/100 と小さくなる場合が示されていた。そこで，これに対応するケースの計算をしてみた。

R = 5; \[Gamma] = 13; \[Beta] = R/\[Gamma]; \[Alpha] = 130*10;
sol = NDSolve[{
x'[t] == \[Beta] x[t] (1 - x[t] - y[t] - z[t]) - x[t]/\[Gamma],
y'[t] == 0.9987*x[t]/\[Gamma] - y[t]/\[Alpha],
z'[t] == 0.0013*x[t]/\[Gamma],
x[0] == 0.01, y[0] == 0, z[0] == 0}, {x, y, z}, {t, 0, 3000}];
fx[t_] := x[t] /. sol[[1, 1]]
fy[t_] := y[t] /. sol[[1, 2]]
fz[t_] := z[t] /. sol[[1, 3]]
Plot[fx[t], {t, 0, 3000}, PlotRange -> {0, 0.5}]
Plot[fz[t], {t, 0, 3000}, PlotRange -> {0, 0.01}]

基本再生産数と免疫消失期間を前回のそれぞれ2.5倍，10倍にしている。

この結果，定常状態における全人口に対する感染者割合は0.77%，一日当たり死亡数は，97人/日となる。年間死亡数は3.5万人なので，季節性インフルエンザの3.5倍程度に収まることになる。

図１：SIDRモデルにおける感染者割合の推移

図２：SIDRモデルにおける累計死亡数の推移

（注）西浦さんの公式だと，$\frac{1-1/R}{1+\alpha / \gamma}\cdot \frac{n}{\gamma} =100$人/日となるので，上記の結果とよく対応している。

SIDRモデル（１）

 エンデミックからの続き　（参考：感染症の数理シミュレーション（２））

西浦博さんが紹介している既知のSIRSモデルは，年齢構造や平均余命まで考慮した精緻なモデルだが，一般ピープルがその振る舞いを試してみるには大層なので，単純なSIDRモデルを考えてみた。本質的に，牧野淳一郎さんが最近計算したモデルと等価であり，念のため死亡者数をちょっと取り出しただけだ。

図１：SIDRの概念図

$\dfrac{du_1}{dt} = -\beta \dfrac{ u_1 \cdot u_2}{n} + \dfrac{u_3}{\alpha}$

$\dfrac{du_2}{dt} = \beta \dfrac{ u_1 \cdot u_2}{n} -\dfrac{u_2}{\gamma_1} -\dfrac{u_2}{\gamma_2}$

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

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

ここで，$u_1+u_2+u_3+u_4=n$は定数で，全国民数である。また，$\alpha$が免疫保持期間（day），$\gamma$は感染期間（day）である。CFRを$c$とすると，$\frac{1}{\gamma}=\frac{1}{\gamma_1} + \frac{1}{\gamma_2}$が成り立ち，$\frac{1}{\gamma_1}=\frac{1-c}{\gamma}，\frac{1}{\gamma_2}=\frac{c}{\gamma}$である。

微分方程式の両辺を$n$で割って，$x=\frac{u_2}{n}，y=\frac{u_3}{n}，z=\frac{u_4}{n}$とする。また，$s=\frac{u_1}{n}$を，$s=1-(x+y+z)$で消去すると次の3変数の微分方程式系が得られる。

$\dfrac{d x}{dt} = \beta (1-x-y-z)\ x -\dfrac{x}{\gamma}$

$\dfrac{d y}{dt} = \dfrac{(1-c)\ x}{\gamma} -\dfrac{y}{\alpha}$

$\dfrac{d z}{dt} = \dfrac{ c\ x}{\gamma}$

この式で$c=0$とすれば，牧野さんの式となる。次のMathematicaコードで検算したところ，彼の結果は再現できた。ここでは，次の値で計算してみる。$R=\beta \gamma =2，\gamma=13，\alpha=130，c=0.0013，x[0]=0.01$

R = 2; \[Gamma] = 13; \[Beta] = R/\[Gamma]; \[Alpha] = 130;
sol = NDSolve[{
x'[t] == \[Beta] x[t] (1 - x[t] - y[t] - z[t]) - x[t]/\[Gamma],
y'[t] == 0.9987 * x[t]/\[Gamma] - y[t]/\[Alpha],
z'[t] == 0.0013 * x[t]/\[Gamma],
x[0] == 0.01, y[0] == 0, z[0] == 0}, {x, y, z}, {t, 0, 1000}]
fx[t_] := x[t] /. sol[[1, 1]]
fy[t_] := y[t] /. sol[[1, 2]]
fz[t_] := z[t] /. sol[[1, 3]]
Plot[fx[t], {t, 0, 1000}, PlotRange -> {0, 0.2}]

その結果は図２のようなものであり，振動の後に平衡状態となり，感染者割合は4.5%に収束する。

図２：SIDRモデルにおける感染者割合の推移

また，400日目を過ぎると，死亡数は図３のように線形で増加するが，このときの1日あたりの死亡者割合は，4.5 x 10^-6 となり，540人/日に相当するので，年間死亡数は19万人に達する。これは季節性インフルエンザの超過死亡数（1万人）の20倍のオーダーである。

図３：SIDRモデルにおける累計死亡数の推移

手計算による第7波ピーク予測（１）

新型コロナウィルスオミクロンBA.5株（BA.4はどこに行った？）による感染者数がこのところ急増している。その一方で，近鉄大和西大寺や祇園祭は人であふれているらしい。

京大の西浦博さんは過去の経験で懲りたのか，ポルトガルの分析グラフを解説するだけで，本邦の数値予測をマスコミには見せていない。KEKの野尻美保子さんも仕事先のヨーロッパで感染して忙しいらしく，以前のような感染者の重症度重み付きグラフを出していない。

名古屋工業大学のグループによるAI予測の結果が，7月8日にNHKで取り上げられていた。それによれば，東京では7月25日に18,000人のピークを迎えるという。実際には，7月15日に19,000人であり，たった1週間先もまともに推定できていなかった。感染症数理モデルの専門家をよんでこないとだめじゃないの。タイトルにAIがあればよいという発想がアウトだ。

そんなわけで，素人が，費用0円，手計算で第7波のピークを予測することにする。

図：新規感染者数の増減率の常用対数値（2021.7.7 - 2022.7.13）

NHKによる都道府県別新規感染者数の日々データから，前週値との新規感染者数の比率を増減率とする。これを1週間移動平均した値の対数値をプロットしたものが上図になる。図のピークは新規感染者数グラフの変曲点を表わし，正の値が続いている間は感染者数が増加している。第6波の増加期間は約60日であり（コロラド博士の9週間=63日ピーク説と符合），この山の増加率の常用対数値はピークを挟んでほぼ1次関数で増減している。

第7波が始まったのは約1ヶ月前の6月中旬であり，7月中旬にはすでに増加率はピークを越えてしまった。増加率の絶対値は前回より小さい（4割程度）が，開始からピークまでが30日というのは前回と同じだ。増加率の減少過程も前回同様だと仮定すれば，8月13日のお盆前まで新規感染者数は増加が続くことになる。

そこで東京について，3つの場合について計算してみる。7月中旬のピーク以降，(A) 3週間で増加率の常用対数値が一様に減少して0になる。(B) 4週間で増加率の常用対数値が一様に減少して0になる。(C) 5週間で増加率の常用対数値が一様に減少して0になる。例えば，7/17 の1.4万というのは前週の新規感染者数の平均値を表わしている。

(A) では，7/17 1.4万，7/24 2.1 万，7/31 2.5万，8/7 2.6万 → 8月第1週でピークの2.6万人/日

(B) では，7/24 2.2万，7/31 2.9 万，8/ 7  3.4万，8/14 3.5万 → 8月第2週でピークの3.5万人/日

(C) では，7/24 2.2万，7/31 3.1 万，8/ 7  3.9万，8/14-21 4.4万 → 8月第3週でピークの4.4万人/日

ということになる。全国に換算する場合は7倍して，それぞれ，(A) 18万人/日，(B) 24万人/日(C) 31万人/日が第7波のピークの値ということになる。

まとめれば，使用する元データは，NHKの新規感染者数から得られた対前週増減率の7日間移動平均だけ。仮定する理論は，その対数値が幅60日の三角波になっている第6波の経験値だけ。その結果は，8月の初旬〜中旬に全国20〜30万人/日でピークとなるということ。

（付録）計算方法：　7/17（東京1.4万人/日）における一日当たり新規感染者数の対前週増加率は，10^0.25=1.8倍である。そこで，(A)では，0.25/3≒0.08として，10^0.17=1.48，10^0.09=1.23，10^0.01=1.02と毎週減少する，(B)では，0.25/4≒0.06として，10^0.19=1.55，10^0.13=1.35，10^0.07=1.17，10-0.01=1.02と毎週減少する，(C)では，0.25/5≒0.05として，10^0.20=1.58，10^0.15=1.41，10^0.10=1.26，10^0.05=1.12，10^0.00=1.00と毎週減少するとして，上記の予測値を得た。

越中大門

 越中大門は富山県射水市にある北陸本線の駅だ。北陸新幹線の開通後は，並行在来線の廃止にともなって第3セクターのあいのかぜ富山鉄道が運営している。

ところで，Unreal Engine 5 は，エピックゲームズというアメリカのコンピュータゲーム・ソフトウェアの開発会社がつくった，ゲームエンジンである。3次元コンピュータグラフィックス生成用のプログラムを3Dエンジンとよんでいたが，この延長上の概念に相当する。

Unreal Engine 5 は2020年に発表された。この機能を利用して，Lorentz Dragoというイタリアの3D環境アーティストが個人でモデリング，テクスチャリング，ライティング，アニメーションのすべてを手がけたものが，越中大門の駅の3DCGアニメーションだ。YouTube で2分49秒のコンテンツだけれど，見た目にはまったく現実の駅を撮影したものといわれてもわからない。途中で昼のシーンが急に夜に変わったところではじめてコンピュータによる合成画像だと気付く。

この技術がメタバースに組み込まれると，視覚上はほとんど現実と区別のつかない仮想世界体験が得られるかもしれない。

越中大門は，子どものころ母の実家の滑川まで里帰りするとき，蒸気機関車が牽引する鈍行列車で通過したなじみ深い駅だ。当時の北陸本線の駅を確認すると次のようになっていて，越中大門はちょうど金沢と滑川の中間あたりにある。

（現ＩＲいしかわ鉄道）金沢−東金沢−森本−津幡−倶利迦羅−（現あいの風とやま鉄道）石動−福岡−西高岡−高岡−越中大門−小杉−呉羽−富山−東富山−水橋−東滑川−滑川

ものごころついてから小学校5-6年ごろまで通ったのだが，西高岡，東富山，東滑川はあまり記憶に残っていない。東滑川はまだ開業していなかったかもしれない。北陸線の電化が進んで蒸気機関車が廃止になるころには，里帰りについていくこともなくなってしまった。

何度かトンネルをくぐるのだが，そのたびに蒸気機関車の煙が入らないよう窓を閉めていた。高岡駅のホームなどで母がアイスクリームを買いに降りるのだが，時間内に帰ってこられるかどうかドキドキしながら待っている小学生であった。

写真：越中大門駅のホーム（あいの風とやま鉄道から引用）

太陽自転

 千葉大学の堀田英之准教授と名古屋大学の草野完也教授が，スーパーコンピュータ富岳を用いて太陽内部の熱対流・磁場を精密に計算することによって，太陽では赤道が北極・南極（極地方）よりも速く自転するという基本自転構造（差動回転） を，世界で初めて人工的な仮説を用いずに再現することに成功した。

これまでのスーパーコンピュータ京では計算可能な解像度が１億点だったのが，富岳では54億点の超高解像度計算が可能となった。これにより，太陽内部の磁場のエネルギーが乱流のエネルギーの最大２倍以上という従来の常識（乱流エネルギーに比べて磁場エネルギーは十分小さい）を破る結果が得られ，同時に，赤道が極より速く回転する差動回転を初めて再現することができた。

京と富岳ではスーパーコンピュータ高速化のためのプログラミング技術がかなり異なるのでそれを克服したのが重要だった。その上で，計算メッシュ数を増やすことによって本質的に異なった物理が得られ，「熱対流の難問」と呼ばれる長年の謎を解決したという興味深い結果だ。

図：解像度と自転速度の関係（千葉大学プレスリリースより引用）

MCMCへの道（３）

 MCMCへの道（２）からの続き

次は，MCMC法#3マルコフ連鎖である。昼食のメニューが３種類に限られていて，なおかつ当日のメニュー選択の確率が，前日のメニューだけによって定まるという例があげられていた。はい，この例は非常によくわかる。

問題は，Juliaの配列と行列の扱いだ。MahthematicaはすべてListと考えればよいのでわかりやすかったが，Juliaでは型として，ArrayもMatrixもある。見様見真似でとりあえずコードを書いたのだけれど，まだ十分理解できていない。このため，3x3行列の配列であるbをベクトルvにかけて得られる配列の転置ができなかったので，手動で転置している。

#using LinearAlgebra
using Plots

x=[zeros(3) for i in 1:11]
y=[zeros(11) for i in 1:3]
#
a=[0.2 0.1 0.3 ; 0.2 0.6 0.5 ; 0.6 0.3 0.2]
b=[zeros(3,3) for i in 1:11]
b[1]=[1 0 0 ; 0 1 0 ; 0 0 1]
v=[0.3,0.2,0.5]
#
for i in 2:11
  b[i]=a*b[i-1]
end

for j in 1:11
  x[j]=b[j]*v
end
#
for i in 1:11
  for j in 1:3
    y[j][i]=x[i][j]
  end
end
#
plot(y,legend=:bottom, xlim=(0,10),ylim=(0,0.6))

MCMCへの道（２）

 MCMCへの道（１）からの続き

さて，MCMC法#2棄却サンプリングである。手順は次の通り。

(1) 目標分布 p(x) を事後分布ともいう。これは既知の関数。

(2) 乱数発生計算が容易な分布 q(x) で， k q(x) ≧ p(x) を満たすものを選ぶ。 kは正定数。例えば，一様分布とか正規分布。

(3) q(x) にしたがう乱数 x' を発生する。

(4) [0, k q(x)] の一様分布にしたがう乱数 y' を発生する。

(5) y' ≦ p(x') ならx' を採用する。そうでなければ棄却する。

この(3) から(5) を必要なだけ繰り返す。

using Plots

function beta(x,a,b)
  return x^(a-1)*(1-x)^(b-1)
end

x= [k for k in 0.0:0.001:1.0]
y= beta.(x,10.2,5.8)
z= [0.00013 for k in 0.0:0.001:1.0]
scatter(x,y,xlim=(0,1.0),ylim=(0,0.00015),
  markershape = :circle,
  markersize = 1,
  markeralpha = 0.75,
  markercolor = :blue,
  markerstrokewidth = 0.1,
  markerstrokealpha = 0.1,
  markerstrokecolor = :white)

g1=scatter!(x,z,xlim=(0,1.0),ylim=(0,0.00015),
  markershape = :circle,
  markersize = 1,
  markeralpha = 0.75,
  markercolor = :orange,
  markerstrokewidth = 0.1,
  markerstrokealpha = 0.1,
  markerstrokecolor = :white)

n=3000
p=rand(n)
r=0.00013*rand(n)
q=ones(n)
for i in 1:n
  f=beta(p[i],10.2,5.8)
  q[i]=ifelse(f <= r[i],0,r[i])
end

g2=scatter(p,q,xlim=(0,1.0),ylim=(0,0.00015),
  markershape = :circle,
  markersize = 1,
  markeralpha = 0.75,
  markercolor = :green,
  markerstrokewidth = 0.1,
  markerstrokealpha = 0.1,
  markerstrokecolor = :white)

plot(g1,g2,legend=:left,aspect_ratio=6000)

図：ベータ分布 p(x) と一様分布 q(x) と定数 k=0.00013 による棄却法の例

MCMCへの道（１）

 そういえば，機械学習をまともに勉強していなかった。いや，そもそも統計学も確率論もあれもこれもである。先日の水素原子波動関数の可視化についても進んでいないのだけれど，どうやら，メトロポリス・ヘイスティング法を使うのが望ましいらしい（by tsujimotter）。そういえば，モンテカルロ法もまともに勉強してこなかった。

そこで，反省してMCMC（マルコフ連鎖モンテカルロ法）の全体像をマスターすべく，いろいろ探したが，おこちゃま向けの適当なテキストがない。こうなるとYouTubeだのみだ。機械学習基礎理論独習というサイトがあったので，このPythonコードをJuliaで再現しながら学ぶことにした。

その１回目はMCMC法#1モンテカルロ法なので，これを再現してみた。

using BenchmarkTools
using Random
using Plots

function findpi(n)
  rng = MersenneTwister(0)
  count = 0
  for i in 1:n
    count += ifelse(rand(rng)^2+rand(rng)^2<=1.0,1,0)
  end
  return 4*count/n
end

#@btime findpi(10^6)
#@benchmark findpi(10^8)

function circle(n)
  x=rand(n)
  y=rand(n)
  c=zeros(n)
  for i in 1:n
    c[i]=ifelse(x[i]^2+y[i]^2<=1.0,0.0,1.0)
  end
  scatter(x,y,
  marker_z = c,
  markershape = :circle,
  markersize = 1,
  markeralpha = 0.2,
  markercolor = :grey,
  markerstrokewidth = 0.1,
  markerstrokealpha = 0.1,
  markerstrokecolor = :white,
  markerstrokestyle = :dot,
  aspect_ratio = 1.0,
  xlim=(0.0,1.0),
  ylim=(0.0,1.0),
  legend=:topright
  )
end

circle(10000)

jupyterlabを使っているが，多数のデータの描画で表示が重くなってかたまる現象がみられる。黒木さんは400msで10^8点のデータが計算できるとしていたが，こちらではメルセンヌツイスターアルゴリズムを用いて300msくらいかかってしまう。

図：単純なモンテカルロ法による円周率の計算

緊急事態宣言解除（３）

緊急事態宣言解除（２）からの続き

昨日，残っていた北海道，東京都，神奈川県，埼玉県，千葉県の緊急事態宣言が解除された。記念に東京都のモデル計算を行ってみた。なんだかグダグダな対応のうちに，よく理由はわからないが，欧米とは異なり，他の東アジアやオセアニア地域なみの水準で第１波はおよそ終息したように思われる。

モデル計算で使用したパラメタは，$\alpha_1=5/0.8,\ \alpha_2=5/0.2,\ \beta=0.505,\ \gamma_1=15/0.94,\ \gamma_2=15/0.06,\ \lambda=35,\ \tau=14,\ \nu=0.1$である。時刻の基準は，新規感染数累計と東京の人口比が10ppmになる時点の3月23日としており，日本の場合とちがって，新規感染数の初期値は正しい値となっている。

図１　東京の新規感染数累計（u6），死亡数累計（u4），回復数累計（u5）

図２　モデル計算に用いた東京（オレンジ）と日本（水色）の実効再生産数

NHKなどで毎日報道されているデータを東京の人口で割って１万人あたりの数に換算したものが上記グラフのサークルであり，下記の値を用いた。なお日時の原点は3/22としている（無駄に有効数字が多くて気持ち悪いのだが面倒なのでそのままにした結果）。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
xt=[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,50,51,52,53,54,55,
56,57,58,59,60,61,62]
yt=[0.1100,0.1221,0.1514,0.1850,0.2136,0.2571,0.3057,
0.3164,0.3721,0.4193,0.4886,0.5521,0.6364,0.7386,
0.7964,0.8529,0.9557,1.0850,1.2179,1.3586,1.4771,
1.5414,1.6564,1.7471,1.8536,1.9957,2.1250,2.2014,
2.2743,2.3621,2.4564,2.5521,2.6664,2.7400,2.7914,
2.8193,2.8993,2.9329,2.9657,3.0836,3.1979,3.2629,
3.3250,3.3657,3.3914,3.4079,3.4357,3.4614,3.4771,
3.5421,3.5621,3.5693,3.5907,3.5971,3.6071,3.6107,
3.6179,3.6214,3.6250,3.6664,3.6686,3.6700,3.6800]
zt=[0.286,0.357,0.357,0.357,0.357,0.357,0.357,
0.857,0.857,0.857,1.000,1.357,1.643,2.143,
2.214,2.500,2.571,2.857,2.857,3.000,3.000,
3.357,3.786,4.000,4.500,4.857,5.071,5.500,
5.786,5.786,6.214,6.643,7.143,7.143,7.571,
7.714,8.357,8.571,9.000,10.071,10.357,10.714,
10.714,11.071,11.429,12.214,12.857,12.857,13.500,
14.000,14.500,15.143,15.643,16.429,16.929,17.214,
17.429,17.643,18.286,18.786,18.786,18.786,18.800]/100
plot(xt,yt,st=:scatter,label="Confirmed-tokyo")
plot!(xt,zt,st=:scatter,label="Deaths-tokyo")
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

緊急事態宣言解除（２）

東京高等検察庁の黒川弘務検事長が軽い訓告（懲戒処分ではない）を受けつつ辞表を提出した日，大阪府・兵庫県・京都府の緊急事態宣言が解除された（北海道・東京都・神奈川県・埼玉県・千葉県は未解除）。

無理無理日本のデータにあわせた単純なモデル計算を実行してみると次のようになった。全回復数累計（無症状感染からの回復を含む）は１万人当り６人強にしかならず，例の0.6%とはどうしても１桁違うのであった。なお，使用したパラメタは，$\alpha_1=5/0.8,\ \alpha_2=5/0.2,\ \beta=0.67,\ \gamma_1=15/0.94,\ \gamma_2=15/0.06,\ \lambda=63,\ \tau=14,\ \nu=0.001$である。$\nu$は想定値の0.01より１桁小さく取ってはじめてこの程度の一致をみているので，真の感染数や死亡数は現在報告された値より10倍程度までの範囲で変化することがありうるかもしれない。そのときには，0.6%が再現できるのことになるのだが・・・

図１　日本の新規感染数累計（u6）と死亡数累計（u4）

図２　日本の総回復数累計（u5）

図３　このモデルパラメタでの有効再生産数（$R_{\rm eff} = \beta(t) \alpha$）

このモデルの$\beta(t)$では，４月中旬に有効再生産数 $R_{\rm eff}$ が１を切るようなパラメタを選択したことになっている。なお，報告時点ベースの観測値を再現するようにモデルパラメタを設定しているので，実際には４月の初旬ということになるだろう。

WHOに報告されているデータを日本の人口で割って，１万人あたりの数に換算したものが上記グラフのサークルであり，下記の値を用いている。なお日時の原点は3/1としている。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

#japan-data(start=3/1)

xj=[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,50,51,52,53,54,55,

56,57,58,59,60,61,62]

yj=[0.0190,0.0202,0.0213,0.0225,0.0252,0.0277,0.0324,

0.0361,0.0387,0.0408,0.0451,0.0492,0.0536,0.0568,

0.0619,0.0646,0.0658,0.0658,0.0693,0.0754,0.0790,

0.0830,0.0864,0.0895,0.0947,0.103,0.110,0.119,

0.134,0.148,0.155,0.173,0.189,0.208,0.232,

0.260,0.290,0.310,0.338,0.378,0.424,0.477,

0.536,0.576,0.607,0.643,0.681,0.728,0.777,

0.822,0.853,0.882,0.912,0.946,0.983,1.018,

1.046,1.062,1.078,1.099,1.118,1.133,1.154]

zj=[0.040,0.048,0.048,0.048,0.048,0.048,0.048,

0.048,0.056,0.071,0.095,0.119,0.151,0.167,

0.175,0.190,0.222,0.222,0.230,0.262,0.278,

0.286,0.325,0.333,0.341,0.357,0.365,0.389,

0.413,0.429,0.444,0.452,0.452,0.516,0.548,

0.556,0.579,0.635,0.643,0.675,0.698,0.746,

0.778,0.810,0.865,0.944,1.079,1.175,1.222,

1.278,1.357,1.476,2.198,2.278,2.516,2.651,

2.762,2.786,2.984,3.087,3.294,3.429,3.603]/100

plot(xj,yj,st=:scatter,label="Confirmed-japan")

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

小田垣さんのSIQR

朝日新聞に，九州大学名誉教授の小田垣さんの計算として「PCR検査を倍にすれば、接触「5割減」でも収束可能？」という記事がでた。twitterでも注目を集めていた。

小田垣さんのホームページにその論考があった[1]。SIQRモデルということで，感染者を2段階に分けていた。我々のSIIDR2モデルと本質的に同じではないか。感受性保持者S(t)から直接我々の重症患者I2(t)，すなわち小田垣さんの隔離感染者Q(t)に遷移する項があって，ここは違うのだが，議論が始まる前の段階でこの項を落としているので結局同じです。違うのはこちらには死亡数D(t)への遷移を含んでいることくらいである。

そこで両方のモデルで使用しているパラメタを比較してみた。左が小田垣さんのSIQR，右が我々のSIIDR2である。
\begin{equation}
\begin{aligned}
N\beta \quad (0.07) &= \beta \quad (0.4-0.6)  \\
p \quad (0.096) &= \dfrac{1}{\alpha_2} \quad (0.04) \\
\gamma \quad (0.04) &= \dfrac{1}{\alpha_1} \quad (0.16) \\
\gamma \quad (0.04) &= \dfrac{1}{\gamma_1}\quad  (0.064) \\
noparameter (0) &= \dfrac{1}{\gamma_2}\quad   (0.00267) \\
\end{aligned}
\end{equation}
感染率の$\beta$が一桁違うことがわかる。その他はファクターが異なるくらいである。I(t)からR(t)への遷移時間とQ(t)からR(t)への遷移時間が同じとしていることもやや疑問に感じる。また，小田垣氏は論考の中で，基本再生産数（実効再生産数）は$N\beta/p$でなく$N\beta/(p+\gamma)$とすべきだと主張しているが，これはどうなのだろうか。

とりあえず，感染初期のS(t)=Nのときにあてはめて，接触5割減の議論に持ち込んでいる部分は興味深いが，実際にこのパラメータを使って日本の感染データの全体像を説明できるのだろうか。我々のコードを少し修正して彼らのモデルとパラメタを使った計算を再現してみた（初期値の細かな調整は我々のモデルを前提としているがとりあえずはそこはそのままにしておく）。

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

sky = @ode_def SIQR_model begin
  du0 = 1 # u0:time
  du1 = -β*u1*u2/n # u1:Noimmunity(Susceptible)
  du2 = β*u1*u2/n -u2/α2 -u2/α1 # u2:(Infected)
  du3 = u2/α2 -u3/γ1 # u3:(Quarantined)
  du4 = u3/γ2 # u4:Dead(not in use)
  du5 = u2/α1 +u3/γ1 # u5:Recovered
  du6 = u2/α2 # u6:Accumulated Quarantined
  du7 = u3/γ1 # u7:Accululated Recovered
end n α1 α2 β γ1 γ2 λ τ

function epidm(β,ν,λ,τ,T)
n=10000.0 #total number of population
α1=1/0.04   #5.0/0.8 #latent to recovery (days/%)
α2=1/0.096  #5.0/0.2 #latent to onset (days/%)
β=0.07  #0.45  #infection rate (/day・person)
γ1=1/0.04   #15.0/0.96 #onset to recovery (days/%)
γ2=15.0/0.04   #onset to death (days/%) (not in use)
u0 = [0.0,n-11ν,4ν,2ν,0.0,5ν,ν,0.0] #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

#japan-data(start=3/1)
xj=[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,50,51,52,53,54,55,
56,57,58,59,60,61,62]
yj=[0.0190,0.0202,0.0213,0.0225,0.0252,0.0277,0.0324,
0.0361,0.0387,0.0408,0.0451,0.0492,0.0536,0.0568,
0.0619,0.0646,0.0658,0.0658,0.0693,0.0754,0.0790,
0.0830,0.0864,0.0895,0.0947,0.103,0.110,0.119,
0.134,0.148,0.155,0.173,0.189,0.208,0.232,
0.260,0.290,0.310,0.338,0.378,0.424,0.477,
0.536,0.576,0.607,0.643,0.681,0.728,0.777,
0.822,0.853,0.882,0.912,0.946,0.983,1.018,
1.046,1.062,1.078,1.099,1.118,1.133,1.154]
zj=[0.040,0.048,0.048,0.048,0.048,0.048,0.048,
0.048,0.056,0.071,0.095,0.119,0.151,0.167,
0.175,0.190,0.222,0.222,0.230,0.262,0.278,
0.286,0.325,0.333,0.341,0.357,0.365,0.389,
0.413,0.429,0.444,0.452,0.452,0.516,0.548,
0.556,0.579,0.635,0.643,0.675,0.698,0.746,
0.778,0.810,0.865,0.944,1.079,1.175,1.222,
1.278,1.357,1.476,2.198,2.278,2.516,2.651,
2.762,2.786,2.984,3.087,3.294,3.429,3.603]/100
plot(xj,yj,st=:scatter,label="Confirmed-japan")
#plot!(xj,zj,st=:scatter,label="Deaths-japan")

β=0.07
ν=0.01
T=60

@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))
- - - - - - - - - - - - - - - - - - - - - - - - - - - - -

結果は以下の通りである。確かに初期の段階ではデータを説明しているが，その後の振る舞いは説明できない。このモデル（のパラメタ）はあまりよろしくないのかもしれない。

図　日本の新規感染数累計の推移（3/1-5/1）u6をデータと比較する

［１］新型コロナウイルスの蔓延に関する一考察（小田垣孝，2020.5.5）
［２］隔離と市中の感染者を分ける SIR モデル（佐野雅己，2020.4.29）
［３］3.11以後の科学リテラシー No. 89（牧野淳一郎，2020 科学５月号）
［４］感染症の数理シミュレーション（８）（2020.3.15）

米国の集団免疫率（３）

米国の集団免疫率（２）からの続き

タイトルはもう変更したほうがいいかもしれない。というのも，抗体検査の結果，ニューヨーク州では12.3%が抗体を持っている（感染済）という結果がでているからだ。こちらの結果とはほぼ１桁違うので，我々のモデルの前提や仮定のどれかががまったく間違っているのではないか。しかし，モデルを検討する余力がないので（遠隔授業の準備で手いっぱい），そこは放置したまま，米国の新しいデータに基づいたパラメタ推定を行う。というのも，これまでの6万人から6.5万人という発表に代わり，再び死亡者が10万人を越えるという予想が出ているからだ。前回と同様に，HEMLのCOVID-19 projections のページを見れば，確かに米国全体の死亡数は8月には13.4万人になりそうだとある。

前回同様のSIIDR2モデルで計算する。使用するWHOのデータ（人口１万人当り）はこれ。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
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.280,1.576,1.929,2.074,2.587,3.136,3.722,
4.268,4.953,5.684,6.483,7.335,8.310,9.327,
10.13,11.03,11.99,12.93,14.00,14.96,15.92,
16.81,17.55,18.33,19.20,20.19,21.10,21.96,
22.80,23.58,24.31,25.19,26.12,27.29,28.28,
29.16,29.85,30.47,31.42,32.39,33.20,34.16]
za=[0.0006,0.0008,0.0009,0.0011,0.0012,0.0012,0.0012,
0.0018,0.0018,0.0030,0.0046,0.0061,0.0061,0.0122,
0.0143,0.0204,0.0268,0.0301,0.0377,0.0506,0.0641,
0.0728,0.0865,0.117,0.146,0.178,0.213,0.254,
0.290,0.329,0.387,0.445,0.504,0.562,0.620,
0.667,0.712,0.785,0.856,0.922,0.984,1.038,
1.089,1.141,1.216,1.284,1.337,1.402,1.456,
1.492,1.532,1.591,1.679,1.742,1.792,1.842]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - -

パラメタはこれ，前回と少し変えている。

#β=0.61,ν=0.12,λ=49,τ=16,α2=5.0/0.20,γ2†=15.0/0.10 

結果はこれ，確かに15〜16万人くらいになりそうだ。

図　米国の感染カーブ（u3=重症感染数，u4=死亡数，u6=新規感染数累計，1万人当）

P. S. 日本でも，神戸の病院の外来患者の約3%に新型コロナウイルス感染症の抗体が検出されたとある。なかなか大きな数字だ。これが1%だとしても，大都市部の数十万人が感染済みということ。

自分で計算できる実効再生産数

新型コロナウイルス感染症専門家会議が５月１日に出した提言・状況分析における実効再生産数のグラフについてあれこれいわれている。こういう場合は自分で確かめておよその感じをつかみたい。いやいやこんな緊急事態に素人が世情を惑わせるよけいな計算をしてはいけないという声も，かつてよけいなことをしまくったその分野の素人（別の分野の専門家）から聞こえてくるのである。ここでComputational Thinkingを標榜するのであれば，計算の自由（Freedom of Computation）を宣言しておきたい。もちろんそのためには最低限データとアルゴリズムを明らかにしておく必要がある（専門家会議はそうしていない）。

①使用するデータ：日本の新規感染数の日次統計（WHOに報告されたもの）WHOのSituation Reportsのデータ（確定感染者数 Confirmed）は次のようになった。ただし，2月22日から5月1日まで（WHOが公表している日付を基準としたもの）の70点のデータである。
{12,27,12,13,7,22,24,20,9,15,14,16,33,32,59,48,33,26,54,52,55,41,64,34,15,15,44,77,46,50,43,39,65,98,96,112,194,173,87,225,206,233,303,351,383,252,351,511,579,658,743,507,390,455,482,585,628,566,390,367,378,423,469,441,353,203,191,276,236,193}

②使用するモデル：単純なSIRモデルを仮定し，感染が人口の数%を越えて拡がっていないものとし，未感染者数（感受性保持者 Susceptible $S(t)$）が全人口（$N$ 定数）とほぼ等しいとする。なお，感染者（Infected $I(t)$）と未感染者の接触による１人１日当りの感染率 $\beta$は対策効果を含めて時間の関数$\beta(t)$とした。さらに，感染期間（日）を$\alpha$として，実効再生産数を$R_t=\alpha \beta(t)$で定義する。このとき，感染者数$I(t)$は次の微分方程式を満足し，$R_t$は$I(t)$とその時間微分から求まる。
\begin{equation}
\begin{aligned}
\dfrac{d I(t) }{dt} = \beta(t) S(t) I(t) / N - I(t) / \alpha \approx \dfrac{ R_t - 1}{ \alpha} I(t) \\
\therefore R_t = 1 + \alpha  \dfrac{d I(t) }{dt} / I(t)
\end{aligned}
\end{equation}
③計算方法：誤差はあるけれど時間の単位を１日とする差分式に直して，エクセルで計算する。もとの新規感染数データのままではゆらぎが大きいので，５日移動平均を求めて$I(t)$とする。これから中心差分（$I(t+1)-I(t-1))/2$）で１日当りの新規感染数の変化分を求める。さらに，このゆらぎを緩和するためにこの５日移動平均をもとめて$d I(t)$とする。これから②の式を用いて$R_t$を求めた。

④結果：図のとおりである。

図１　全国の新規感染数の推移（2/22-5/1）

図２　全国の実効再生産数の推移（2/24-4/29）

図２には社会的な事象を書き加えているが，WHOへの報告公表時点とは時間的なずれ（１週間〜10日）があることに注意する。→（注）12日ぐらいか（P. S. 参照）

⑤ 結論：１ヶ月前に1.5〜2.0近くまであった$R_t$が現在0.5〜1.0の範囲にまで落ちてきたことがみてとれ，これは専門家会議の結果とおおよそ近いものである。しかし，細かく見ればそのピークの位置や高さは異なっている。とくに，2月末から3月中旬までの結果は待ったく違う。そもそも出発点となる新規感染数のデータもWHOに報告されたものは，１週間の大きな周期構造を持っており，牧野さんのいうように（ちょっと意味が違うかもしれないけれど），専門家会議が確定感染者数（Confirmed）からどのように処理して新規感染数の推定値を導いているのかがはっきりしないのでもやもやが残る。

P. S. 高橋健太郎さん（@kentarotakahash）によれば，「潜伏期間+発症〜報告までの日数の平均を12日と考えてみたら、かなり納得できました。15日目あたりから一度、下降してRe=1を切りますが、2/24+15-12ですので、2月27日の休校要請の効果に見えます。その下降のピークは3月4日頃。一週間で緩んで、再び上昇が始まった。3月10日頃にRe=2を越えるピーク。」だそうだ。またさらに，「その後、Re=1.5以上の状態が続きますが、三連休の人手の反省があり、3月25日の都知事緊急記者会見を経て、3月27日頃から下降に転じる。3月31日にはRe=1を切り、専門家会議の4月1日には1を下回ったという分析と合致します。」

そして，「その後一度、上昇して、1を越えますが、4月5日頃から本格的な下降が始まる。これは緊急事態宣言が出るというムードの先取り。が、4月14日あたりから再び上昇の兆しというところでしょうか。」なるほどなのだった。

［１］新型コロナウイルス感染症対策専門家会議の見解等（牧野淳一郎，2020.04.18）

［２］新型コロナウイルス感染症対策専門家会議の見解等その２（牧野淳一郎，2020.05.12）
［３］感染症数理モデル（北海道大学医学統計学教室のSqquential SEIRモデル）
［４］Rt-COVID-19 Japan （都道府県別新型コロナウイルスの実効再生産数）
［５］山中伸弥による新型コロナウイルス情報発信
［６］A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics（Cori et al.）
［７］新型コロナ対策専門家会議が判断の拠り所にしている『実効再生産数・倍加時間』の算出方法に関する考察（@makirin1230  2020.05.06）

米国の集団免疫率（２）

米国の集団免疫率（１）からの続き

トランプは4/17の会見で米国におけるCOVID-19の死亡数は6〜6.6万人にとどまるとした。これは勝手な想像ではなく，IHMEの最新の予測である。前回の10〜25万人から減少させたことを自分の政策の成果であるかのようにアピールしつつ，ロックダウンを解消して経済回復を誘導しようという意図に基づくものだろう。

前回のSIIDR2モデルの適用はちょうど2週間前だったので，あらためてこのモデルと上記の主張を組み合わることで米国の集団免疫率を推定してみる。前回のように1ppm到達の基準日3/10から4/17までの39日分の新規感染数累計と死亡数累計の人口比データを示す。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
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,25,26,27,
28,29,30,31,32,33,34,35,36,37,38]
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.280,1.576,1.929,2.074,2.587,3.136,3.722,
4.268,4.953,5.684,6.483,7.335,8.310,9.327,
10.13,11.03,11.99,12.93,14.00,14.96,15.92,
16.81,17.55,18.33,19.20]
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,17.8,21.3,25.4,
29.0,32.9,38.7,44.5,50.4,56.2,62.0,
66.7,71.2,78.5,85.6]/100
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

SIIDR2の計算において次のパラメタを用いると上記の米国のデータが再現できる。
$\beta = 0.60, \nu =0.12, \lambda=28, \tau=16$, $\alpha_1=5/0.80, \alpha_2 = 5/0.20$, $\gamma_1 = 15/0.95, \gamma_2 = 15/0.06$。前回と異なり，$\gamma_2$の値は，中国や韓国などを説明した値の方にややに戻している。

　図１　米国の感染カーブ（u3=重症感染数，u4=死亡数，u6=新規感染数累計）

図２　米国の感染カーブ（同上，u5=回復（免疫獲得）数）

① 米国では重症感染数のピークを迎えている。

② 最終的な死亡数は6〜7万人程度になる。

③ 第１回目の終息が想定される２ヶ月後の集団免疫率は1〜2%のオーダーである。

緊急事態宣言と接触制限モデル（３）

緊急事態宣言と接触制限モデル（２）からの続き

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）

緊急事態宣言と接触制限モデル（２）

緊急事態宣言と接触制限モデル（１）からの続き

緊急事態宣言がでた翌朝（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$

緊急事態宣言と接触制限モデル（１）

北海道大学の西浦博さんが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}}]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

計算の結果，次のグラフが得られた。西浦さんのものと概ね一致する。

図１　SIRモデル＋シグモイド関数制限措置による新規感染者数の変化

　（縦軸の単位は万人，横軸の単位は日＝3/20を基準とした経過日数,
青：制限なし，オレンジ：80%に制限，緑：20%に制限）

図２　上記に，60%，40%，30%に制限などの場合を加えたもの

なるほど，中途半端な制限ではだめだということか。30%に制限するあたりが増加かどうかの臨界点かもしれない。この度の緊急事態宣言はどの程度の効果があるのだろうか。