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

§1. 不思議の国の γ

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

ClearAll["Global`*"]
(*時系列リストdata={v1,v2,...vn}の変域tをbrに従って分割し,dataの区分近似直線関数を作る．関数のリストを作成する関数がfuncs,
傾きのリストを作成する関数がslopes. funcsを組み合わせdataの折れ線近似を作る関数が multifit. 
slopesを組み合わせdataの傾きの棒グラフを作る関数が multislope. 共にtの関数.*)
zip[A_List,B_List]:=Partition[Riffle[A,B],2]
slopes[data_,br_]:=Module[{a,b},Table[a/.FindFit[zip[Range[br[[k]],br[[k+1]]],Take[data,{br[[k]],br[[k+1]]}]],a x+b,{a,b},x],{k,1,Length[br]-1}]]
pieces[data_,br_]:=Table[Fit[zip[Range[br[[k]],br[[k+1]]],Take[data,{br[[k]],br[[k+1]]}]],{1,x},x],{k,1,Length[br]-1}]
multifit[data_,br_]:=Which@@Flatten[Table[{br[[k]]≤x≤br[[k+1]],pieces[data,br][[k]]},{k,1,Length[br]-1}]]/.{x->t}
multislope[data_,br_]:=Which@@Flatten[Table[{br[[k]]≤x≤br[[k+1]],slopes[data,br][[k]]},{k,1,Length[br]-1}]]/.{x->t}

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)

1. Johns Hopkins 大学 (GitHub) から DownLoad

Johns Hopkinsからのデータでは recovered は 死者を含まないが， SIR理論では 死者も含む．以下
「r」 または「rec」と書くと SIR理論の回復者の意味とする．日本語で「回復者」と書くと，死者を含まないとする.
ConfirmedCaseは感染者の積算データなので，infected= active case（感染者.以下「inf」) は,
それから回復者と死者を引いたものになる. さらに(回復者＋死者)が50人になった日から改めて数える．
以下,day1→recが50人以上となった日のindex. dayl→ inf<10 または最新データ獲得日のindex.
inf0,rec0,deaths0,cases0,days0→全てのデータ.  inf,rec,deaths,cases,days→day1からのデータ,  
period→day1からdaylまでの日数.

casesGlobal=Import["https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv"];
deathsGlobal=Import["https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv"];
recGlobal=Import["https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv"];
countries=Drop[casesGlobal[[All,2]],1];

nation="Japan";(*上のリストから. United States は US. UK は United Kingdom.*)

cases0=casesFromJH=Plus@@(Drop[#,4]&)/@Select[casesGlobal,Part[#,2]==nation&];
recFromJH=Plus@@(Drop[#,4]&)/@Select[recGlobal,Part[#,2]==nation&];
deaths0=deathsFromJH=Plus@@(Drop[#,4]&)/@Select[deathsGlobal,Part[#,2]==nation&];
rec0=recFromJH+deathsFromJH;
inf0=cases0-rec0;

day1=FirstPosition[ rec0-50,_?Positive][[1]];(*rec>50からカウント*)
dayInfLess10=Flatten[Position[ inf0-10,_?Negative]];
dayl=Min[Select[dayInfLess10,#>day1&]-1,Length[cases0]];(*inf<10になると中止*)
datesGlobal=Drop[casesGlobal[[1]],4];
period=dayl-day1+1;
inf=Take[ inf0,{day1,dayl}];
rec=Take[ rec0,{day1,dayl}];
cases=Take[ cases0,{day1,dayl}];
deaths=Take[ deaths0,{day1,dayl}];
dates=Take[datesGlobal,{day1,dayl}];

2. Log(rec)の傾き (Subscript[m, r])を求める.

「移動平均+差分を使う方法」と「全体のFit関数の微分(oneFit)」と「下のbreaksに従って期間を分割して
Fitを使う方法(multifit)」を比較する.  また breaksは，Default では「始点はt=1, 幅5，最後のみ少し長め」に設定.
しかし グラフを目で確認しながら設定可能.(但し最初は1,最後はperiodにする)

haba=5;(*←幅を設定*)
breaks=Table[1+haba k,{k,0,(period-1)/haba-1}]//Append[period];(*←グラフを目で見て直線部分を見つけると，更に適切*)
mhaba=Round[Mean[Take[RotateLeft[breaks]-breaks,{1,Length[breaks]-1}]]];
Print[Style["breaks(解析期間の分割)="<>ToString[breaks],Larger,Blue]]

logRmultifit[t_]=multifit[Log[rec],breaks];(*Log(rec)の折れ線近似*)
mr[t_]=multislope[Log[rec],breaks];(*(Subscript[m, r])の棒グラフ*)
mrByOneFit=D[Fit[Log[rec],Table[t^k,{k,0,10}],t],t];(*全体にFitを使用し微分-onefit*)
mrByMovingAverage=Differences[MovingAverage[Log[rec],5]];
GraphicsRow[{DateListPlot[TimeSeries[rec,{dates}],PlotLabel->nation<>": recovered (recが50人以上になった日から)"],
DateListLogPlot[TimeSeries[cases,{dates}],PlotLabel->nation<>": Log(rec) (recが50人以上になった日から)"]},ImageSize->700]//Quiet
GraphicsGrid[{{Show[{ListLinePlot[Log[rec],PlotStyle->{Thin,Blue}],Plot[logRmultifit[t],{t,1,period},PlotStyle->Red]},PlotLabel->"recの対数グラフ"<>",青:data,赤:multifit(平均区間幅"<>ToString[mhaba]<>")"],ListLinePlot[{Log[rec],MovingAverage[Log[rec],5],MovingAverage[Log[rec],10]},PlotLabel->"Log[rec]の移動平均(0,5,10日)"]},
{Plot[{Tooltip[mr[t],"multifit"],Tooltip[Evaluate[mrByOneFit],"onefit"]},{t,1,period},AxesOrigin->{0,0},ExclusionsStyle->Lighter[Purple],
PlotLabel->"Log[rec]の傾き (Subscript[m, r]):"<>"　multifit(平均区間幅"<>ToString[mhaba]<>")とoneFit"],
ListLinePlot[mrByMovingAverage,PlotLabel->"Log[rec]の傾き (Subscript[m, r]): 5日の移動平均と差分"]}},ImageSize->700]

3つのグラフの傾向は同じだが，移動平均は,やや変動が大きく,かつ変域が左にずれる. onefit と multifit を比べると 
onefit の方は両端が大きくハズレる．次章で扱う log(cases)の傾き (Subscript[m, c]) と違い　はっきり言ってどちらを採用してもよいのだが
 (Subscript[m, c]) に合わせて「階段関数 multifit 」を (Subscript[m, r]) に採用したい．即ち (Subscript[m, r])=mr[t] とする．

3. (Subscript[m, r])からγを求める.

SIR理論では
                               ds/dt=-βs*i,    di/dt= βs*i-γi,   dr/dt =γi　　（β,γは定数)
一番右の式から  ((d(log r))/dt)=(1/r)・(dr/dt)=((γ*i)/r).　よって y=log r のグラフの傾きを (Subscript[m, r]) とすると, γ=(Subscript[m, r])・(r/i)

DateListPlot[TimeSeries[#,{dates}]&/@{inf/cases,rec/cases,rec/inf},
PlotLegends->{"(inf/cases)","(rec/cases)","(rec/inf)"},PlotRange->{0,1},PlotLabel->"i,r,c の比"]

γ=(Subscript[m, r])*r/i だから, これからγのグラフが描ける.
r/i は短い期間では一定としてよいはずなのだが，最近の急上昇を見るとそうも言えない．
したがって r/i は平均で考えずに時間tにおける値で考える．

γTable=Table[mr[t],{t,1,period}]*rec/inf;(*γ=(Subscript[m, r])*r/i. 区分関数*)
DateListPlot[TimeSeries[γTable,{dates}],PlotLabel->"γのグラフ"<>"　(平均区間幅"<>ToString[mhaba]<>")",AxesLabel->"γ"]

本来は γは 「回復日数の逆数＝感染期間の逆数」 なので，定数と考えられる （統計のとり方の違いで国には依るだろう) 
しかし上のグラフを見ると，定数と見ることは難しい．(5/11では定数と見ることができたが 5/19にはもはや無理.)
とりあえず直近を除いて平均をとってみると，

γ=TrimmedMean[γTable,{0.2,0.2}]//N;
Print[Style["単純平均によると {γ, 回復（感染)日数}={"<>ToString[γ]<>","<>ToString[1/γ]<>"日}",Larger,Blue]]

上の回復(感染)日数は長過ぎる. コロナに2ヶ月も罹っているとは考えられない．
一方，直近の値では γ=0.13(回復日数=8日)
すごい違いがある．この様に rの変化だけでは γ の値を決めようがないが，
次章で微分方程式と組み合わせて γと (Subscript[R, t]) を一度に求めてみたい．

[注]この章を最初に書いた5/11時点では γ の急上昇がなかったので γ=0.016だった．そのとき
「実際には r が,統計データより大きい可能性が高い．実際,回復者のデータが何日も更新されない国もあり
worldodata.org では回復者のデータは置いてないので 十分可能性がある」と書いて
その理由をもう一つ上げた．以下は γ=0.016(5/11時点)に対して, 書いた文である．

4.  r は統計データより多いはず (致死率との関係)

( γ=0.016(5/11時点)に対して, 書いた文である．回復者は死者を含めないとする．但しrec は 死者を含めるとする. )

DateListPlot[TimeSeries[{(rec-deaths)/deaths},{dates}],PlotLabel->"回復者/死者 "]
TrimmedMean[N[(rec-deaths)/deaths],{0.2,0.2}];
Print[Style["回復者/死者="<>ToString[%],Blue,Larger]]

上を見ると 回復者は死者の6.5倍しかないが，日本の致死率 は 2% ぐらいである．
ところが回復者が死者の6.5倍なら，致死率はもっと高いはずである．(1/7.5=13%)
これは 回復するのに掛かる平均日数を p日, 死亡するのに掛かる平均日数をq 
としたとき（p=回復日数=1/γ ), 十分に p>q なら説明できる.

sol=FindFit[Log[cases],m t+n,{m,n},t];
g[t_]=(m t+n)/.sol;
a=Exp[m]/.sol;
GraphicsRow[{Import["http://mixedmoss.com/miscellaneous/MathematicaSIRModel/deathRatio.jpg"],
DateListPlot[TimeSeries[#,{dates}]&/@{Log[cases],Table[g[t],{t,1,period}]},PlotLabel->"y=log(cases)と直線近似y=g[t]"]},
ImageSize->600]
Print[Style["log(cases)="<>ToString[g[t]]<>", cases=c*(a^t)(a="<>ToString[a]<>")",Blue,Larger]]

Confirmed cases は基本的に指数関数なので cases=c*(a^t)=f(t) (a,c定数)  とおける．(5/11時点では)
p>q のとき，回復者/死者が6.5, 死亡率が2%となることから,
                              6.5= 回復者/死者=  ((f(t-p)*0.98)/(f(t-q)*0.02))=49*(a^(q-p))                 
                                    故に    p-q=(Subscript[log, a])(49/6.5) 
                            
 casesのデータから a を概算で求めると  a=1.07 で 
                                     p=1/γ=1/0.0166=60日,
　　　　　　　　　　　　　q=p-(Subscript[log, 1.07])(49/6.5)=30 日
となる．これは おかしい．故にγは 0.0166より大きくないといけない.   

 
 ［結論］基本的な考え方は今でも正しいと思っているが，5/19時点ではもはや「cases=c*(a^t)」が成り立たない．
(右上図で大きく直線からずれている)したがって上の式そのものは改良しないと成り立たない．

何れにせよ γの値は0.016より大きいことは間違いない．正確な値は微分方程式を解き，総合的判断も入れて 次章で求めたい．

5. 次章のため「nation, dates , breaks ,cases,rec,deaths」 を Export

Export["nation.txt",nation];
Export["numbers.csv",{dates, breaks, cases,rec,deaths}];