Mathematica(ver11)+SIR理論 による 
Covid-19 のデータの分析

§3 . 微分方程式で 未来を見る

テーマはあくまでも「Mathematicaによるデータ分析の試み」です．Mathematica 
に対しても Covid-19 に関しても専門家ではないので,  間違いがあるかもしれません.  ^_^;   　　
2020年5月13日(5月19日改定).    生越 茂樹.  oh.shigeki@gmail.com

0. SIR理論

s=Susceptible(感染可能人数), i=Infected(感染者), r=Recovered(回復者，死者含む), c=Confirmed Case=i+r とすると　
(この r と (Subscript[R, t]) は全く異なる.また s,i,r,c　をプログラム中では sus,inf,rec,cases　とも書く) ,SIR理論では,
                                          (ds/dt)=-βs*i,    (di/dt)=βs*i-γi,    (dr/dt)=γi
Nを全人口とすると i+r+s=c+s=N だから， 
                                   (ds/dt)=-βs*i,   (dc/dt)=βs*i,    (di/dt)=βs*i-γi,    (dr/dt)=γi           ···(#1)
故に (Subscript[R, t])=((sβ)/γ)と定義すると ( これから (Subscript[R, 0])=((Nβ)/γ), (Subscript[R, t])=(Subscript[R, 0])*(s/N) となる),   
                           (ds/(d(γt)))=-(Subscript[R, t])*i,    (dc/(d(γt)))=(Subscript[R, t])*i,    (di/(d(γt)))=(Subscript[R, t])*i-i=((Subscript[R, t])-1)i,     (dr/(d(γt)))=i     ···(#2)
さらに 時間の単位を1/γ日に変えて γt=t’とすると(例えばγ=0.1のときは 10日ごとの変化を考えると)， 
                                       (ds/(dt'))=-(Subscript[R, t])*i,     (dc/(dt'))=(Subscript[R, t])*i,    (di/(dt'))=((Subscript[R, t])-1)i,     (dr/(dt'))=i    
(Subscript[左の2式からR, t]) は 感染期間(回復期間) 毎に, 1人の感染者が2次感染させる人数(即ち,1人の感染者が2次感染させる人数)となる事がわかる．
(#2より (Subscript[R, t])=dc/dr となることからも同じ事が分かる)．
また右から2番めの式から  (Subscript[R, t])>1のとき I は増加，(Subscript[R, t])<1のとき I は減少と分かる．

さらに (Subscript[m, r])=((d(log r))/dt),(Subscript[m, c])=((dLog[c])/dt) と定義すると,

            (Subscript[m, r])=((d(log r))/dt)=(1/r)*(dr/dt)=((γ*i)/r)     ∴ γ= ((r*(Subscript[m, r]))/i),          (Subscript[m, c])=((dLog[c])/dt)=(1/c)*(dc/dt)       ∴  (dc/dt)= c*(Subscript[m, c]).   ···(#3) 
故に                                               
                                                  (Subscript[R, t])=((sβ)/γ) (=(1/(i*γ))(dc/dt))=((c*(Subscript[m, c]))/(i*γ)) =((c*(Subscript[m, c]))/(r*(Subscript[m, r]))) (=(dc/dr))   ···(#4)

オーバーシュートが起こると (Subscript[R, t])= (Subscript[R, 0])s において, sの変化は無視できなくなる．即ち sが減少するに連れ (Subscript[R, t ]) の値も小さくなり 流行は収束する．
しかし (Subscript[R, t])=(Subscript[R, 0])s はsを含み tの式では与えられない．従って, (Subscript[R, t ])も含めた微分方程式を, (Subscript[R, t])’[t]= (Subscript[R, 0])*s’[t], (Subscript[R, t])[0]= (Subscript[R, 0]) という条件を加えて解く． 
また 上の(#1)の式は両辺をN(全人口)で割り s/N, i/N, r/N を, 改めてs,i,r としても成り立つ．以下 s,i,r はこの様な比を表す．

1. (Subscript[R, 0])の値を変えて 微分方程式を解く

前章で述べたように γ=0.05(回復期間20日)とします． 
γ の値はすぐ下の行で簡単に変えられます.    しかしそもそも(#3)から 「時間単位を1/γ に変えると方程式はRの値のみで決まる」ので, 
1つのγでだけチェックすれば十分と言えます.　　またデフォルトでは対象国は日本で，s,i,r などの初期値は5月8日現在のデータ．
 データを変える場合は，下のプログラムの「caseLatest,deathLatest,recLatest, population」 を変更して下さい．
 (Mathematica notebookの場合. cdfでは固定．しかし 後述してますが，R≥1.1 の時は 初期値には殆ど無関係です．)
 ちなみに (Subscript[R, 0])=2 からのスタートです．

DynamicModule[{sol,Rt,s,i,r,contentR,contents,contenti,contentr,content1,myrange},
Manipulate[Quiet[
sol=NDSolve[{Rt'[t]==Rt0*s'[t],s'[t]==-γ*Rt[t]*i[t],i'[t]==γ*(Rt[t]-1)*i[t], r'[t]==γ*i[t],
Rt[0]==Rt0,
s[0]==(population-caseLatest)/population,
i[0]==(caseLatest-recLatest)/population,
r[0]==recLatest/population},{Rt,s,i,r},{t,1,n}];
contentR=If[showR&&!ionly,Evaluate[Rt[t]/.sol]];
contents=If[shows&&!ionly,Evaluate[s[t]/.sol]];
contenti=If[showi&&!ionly,Evaluate[i[t]/.sol],If[ionly,population*Evaluate[i[t]/.sol]]];
contentr=If[showr&&!ionly,Evaluate[r[t]/.sol]];
content1=If[!ionly,1];
myrange=If[ionly,All,Automatic];
Plot[{content1,Tooltip[contentR,"再生産数 (Subscript[R, t])"],Tooltip[contentr,"回復者&死者"],Tooltip[contents,"感染可能な人"],
Tooltip[contenti,"感染者"]},{t,0,n},PlotStyle->{Dashed,Line,Line,Line,Line,Line},
PlotLegends->{"y=1","再生産数 (Subscript[R, t])","r(回復者&死者)","s(感染可能な人)","i(感染者)"},ImageSize->600,PlotRange->myrange]],

{{ionly,False,"i のみを人数表示"},{True,False},ControlType->RadioButton},Delimiter,
{{showR,True,"(Subscript[R, t]) を表示"},{True,False}},{{shows,True,"s (人口比) を表示"},{True,False}},
{{showi,True,"i (人口比) を表示"},{True,False}},{{showr,True,"r (人口比) を表示"},{True,False}},
{{Rt0,2,"(Subscript[R, 0])"},0.5,5,0.01,Appearance->"Labeled"},{{n,Floor[20/γ],"期間"},Floor[5/γ],Floor[200/γ],1,Appearance->"Labeled"},
Initialization:>(
γ=0.05;
caseLatest=15575;(*この後3つの変数は日本の5月8日現在のデータ*)
deathLatest=590;
recLatest=5146;
population=127484450
),FrameLabel->"(Subscript[R, 0])の変化 と s,i,rの増減 "]]

2.  (Subscript[R, 0])>1のとき, 「感染者の最大値」 を求める

(Subscript[R, t])=((Sβ)/γ) と定義すると (これから (Subscript[R, 0])=((Nβ)/γ), (Subscript[R, t])= (Subscript[R, 0])*(S/N) = (Subscript[R, 0])*s  となる),   
                              (ds/dt)=-(Subscript[R, t])*i*γ···(1),     (di/dt)=((Subscript[R, t])*i-i)*γ=((Subscript[R, t])-1)i*γ···(2),      (dr/dt)=i*γ···(3)

この式はs, i, r が全て人口に対する比率のときも成り立つ. (1),(2)より 
                                                            di/ds=-1+ (1/(Subscript[R, 0]))×(1/s)
これを積分して,                                  i+s -(1/(Subscript[R, 0]))log s=一定
「t=0の時の i,s,rの値を(Subscript[i, 0]),(Subscript[s, 0]),(Subscript[r, 0]) とすると,     i+s-(1/(Subscript[R, 0]))log s=(Subscript[i, 0])+(Subscript[s, 0])-(1/(Subscript[R, 0]))log (Subscript[s, 0]) ···(4)
(2)より (Subscript[R, t]) >1 のとき, i が最大⇔(Subscript[R, t])=1 (グラフからも読み取れる)であり，「(Subscript[R, t]) =(Subscript[R, 0]) s 」だから,　
                                                                 (Subscript[s, (iが最大の時)]) =  (1/(Subscript[R, 0]))   
(4)とあわせて   i の最大値は,
                                                   (Subscript[i, Max])=1- (1/(Subscript[R, 0]))(1+log(Subscript[R, 0])) +δ     ( δ=-(Subscript[r, 0])-(1/(Subscript[R, 0]))log (Subscript[s, 0]) ) 
 初期の段階では (Subscript[s, 0])=1, (Subscript[r, 0])=0 と考えて良いから, その近似値は 
                                                    M((Subscript[R, 0]))= 1- (1/(Subscript[R, 0]))(1+log(Subscript[R, 0]))

GraphicsRow[{Plot[{Tooltip[(1-1/r(1+Log[r])),"M((Subscript[R, 0]))"],Tooltip[Log[r]/r^2,"M'((Subscript[R, 0]))"]},{r,1,5},
PlotLabel->"R>１の時, M((Subscript[R, 0]))=1-(1+(Subscript[logR, 0]))/(Subscript[R, 0])(青) と, その微分(緑)",AxesLabel->{"(Subscript[R, 0])","M"},PlotLegends->{},  AxesOrigin->{0,0}],Plot[Evaluate[Table[Tooltip[Exp[0.05(r-1)t],"(Subscript[R, 0])="<>ToString[r]],{r,1,0.3,-0.1}]],{t,0,5/γ},PlotRange->{0,1},
AxesLabel->{"t 日後",Style["i(t)/i(0)",Larger]},PlotLabel->"(Subscript[R, 0])≤1 の時の i の変化 (R=0.3~1)"]},ImageSize->700]

例えば　現在（5月13日), 私の計算では 世界の (Subscript[R, 0])は　1.3～1.4 でした．　
                          M(1)=0,     M(1.1)=0.004 ，M(1.2)=0.015,  M(1.3)=0.029,  M(1.4)=0.045

Clear[r]

Manipulate[Round[1-(1+Log[r])/r,0.001],{{r,2," (Subscript[R, 0]) "},1,5,0.01,Appearance->"Labeled"},FrameLabel->"M((Subscript[R, 0]))の値"]

また 世界人口は77億人, Confirmed Cases(感染者積算)=460万,　死者は30万, 回復者(死者含まず) は 170万，なので, 
　　　　　　　　　　　　　　　 (Subscript[i, 0])=0.00034,  (Subscript[r, 0])=0.00026, (Subscript[s, 0])=0.9994 
故に  1<(Subscript[R, 0]) の時は， δ<0.00075 です． 従って  1.1≤(Subscript[R, 0]) の時は, δ は無視できます.

M((Subscript[R, 0]))の値から，感染者は 
            (Subscript[R, 0])=1.4 のとき 3.5億人,   (Subscript[R, 0])=1.3 のとき 2.2億人,   (Subscript[R, 0])=1.2 のとき 1.1億人,   (Subscript[R, 0])=1.1 のとき 0.3億人
仮に致死率2%とすると, 死者は 
              (Subscript[R, 0])=1.4 のとき 700万人,   (Subscript[R, 0])=1.3 のとき 400万人,   (Subscript[R, 0])=1.2 のとき 200万人,   (Subscript[R, 0])=1.1 のとき 65万人

M((Subscript[R, 0]))は(Sqrt[e])≐1.648で極大となるので，(Subscript[R, 0])=1.65前後では, (Subscript[R, 0])の削減の効果が特に大きいです.

ですから，世界全体にとっては今が頑張りどきと言えるでしょう．

一方 (Subscript[R, 0])<1 の時は i　は減少関数なので,t=0のときiは最小です．そしてこの時 s≐1のままだから (Subscript[R, t])≐(Subscript[R, 0 ])で(2)より,
                                               (di/dt)=((Subscript[R, 0])-1)i*γ
故に                                     i[t]=i[0]*exp[γ((Subscript[R, 0])-1)t] (右上図)

そして 6割の削減で R=1,7割の削減で R=0.75, 8割の削減で R=0.5　となり，確かに8割の削減では急速に収束します．
しかし 以上はあくまでもSIR理論による理論値です．

3.  果たして SIR理論は正しいのか？

まず私は専門家ではないので，以下の事は単なる個人的意見です．批判的に読んで下さい．
先ず感染初期では　sの値は人口比でほぼ1なので　(Subscript[R, t]) =(Subscript[R, 0])=((Nβ)/γ) で一定のはずです．
しかし実際には各国も対策をするので一定でなくても構いません．しかし (Subscript[R, t])　が突然変わると
過渡期の数週間は1(Subscript[次関数的にR, t])が下がり，その後(2~3週間)で (Subscript[R, t]) の値は一定になるはずです．
しかしドイツ，イタリアなどがロックダウンを始めてから，1月以上立ちますがまだ　(Subscript[R, t]) は一定になっていません．
なだらかに減少してます．これはSIR理論では説明できないと思います．


SIR理論というのはそもそも長期的な理論で，短期的な事を解決する理論ではないと感じています．
感染初期でs=1とみなして良い時は (#1)は
　　　　　　　　　　　　　　 (di/dt)=(β-γ)*i,    (dr/dt)=γi (β,γ定数)
となるので「 iもrも指数関数という事」以外には何も言っていません．(β,γのこの世での意味を除けば)
「短期では指数関数的に増える」というのは,生物の増加なら当たり前で高校でも教えているはずです．
したがって短期では SIR理論は何者でもありません．長期的な対象に対して初めて
威力を発揮する理論だと思います．

ある国で (Subscript[R, 0])を抑えるのに失敗して「放ったらかし」
(或いは何か対策をとっているとしてもそれが  (Subscript[R, 0])<1 にするのに成功していない)
となった時，SIR理論が正しいか否かは，確実にわかると思います．
もちろん そんな事期待する訳には行かないので，大掛かりな動物実験でもするか，
マウスの大増殖などの観察しかないと思います．
(これを書いた後で Wikiを読むと, 
「1905–06年のボンベイにおけるペスト流行のデータをうまく再現することが知られている」
とありました．長期の場合は検証されているようです.)

「私の私見」SIR理論は　ボンベイでのペストの流行で正しいことが証明されているようだが
長期的な理論なので，このCovid19の流行を分析する”だけ”では，現時点においては
正しいか否かの結論は出せない．また初期を説明する理論ではなく，
ロックダウン1ヵ月後でも (Subscript[R, t])が減少することの説明もできない．
しかしとても簡潔でわかりやすくきれいなモデルだと思う．



最後は本当に個人的なことです．3月末にオーストラリアから帰国してから，
（普段見ない)　朝のワイドショーなどを見ていると，本当に恐怖にかられて
しまっていました．思えば 3/11 の時も同じでしたが「恐怖を煽って視聴率を稼ぐ」手口です．
一言でいうと「恐怖洗脳」されていました．

その時 SIR理論を知り，Mathematica を使って Covid19の分析を進めることで
再び自分で自分の感情をコントロール出来るようになりました．

その点では本当に Mathematica　には感謝しています．
この後も Mathematicaや他の数学ソフトを使って
いろいろな現象を分析，可視化してみたいと思っています．
 （ちなみに私の一番使っている数学のソフトは GeoGebraです．
Mathematica は2番めです．また現在孤軍奮闘中の数学ソフトは Coq です.)


p.s. Wolfram　Japan の方へ
もしこれを読んでいてくださったら，こんな素晴らしいソフトを日本で販売してくださる事に感謝します．
ただ日本語で長い Text を入力すると　2秒!! 前後の遅延がでます．これを是非なんとかして頂きたいです．
Editorで書いて貼り付けないといけないなんて，信じられません．
また 「pdf ソフト」（縮小したいのでソフトを使うしかない)を使うと「超頻繁に」(50%以上の確率で)
クラッシュします．更にいうと，印刷も何もしないのに突然クラッシュすることも有ります．
Mathematica は イタリア車を目指しているのでしょうか？　日本車の良さも身に着けて欲しいです．

この2点の改良を”是非”お願いします．