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

§2. 拡大再生産数 (Subscript[R, t]) と γ を求める

テーマはあくまでも「Mathematicaによるデータ分析の試み」です．Mathematica に対しても
 Covid-19 に関しても 専門家ではないので,  間違いがあるかもしれません.    ^_^;   　　
2020年5月13日(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. 前章で Export したデータ(nation,breaks,case,rec,deaths)の Import

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

inf=cases- rec;
period=Length[cases];
mhaba=Round[Mean[Take[RotateLeft[breaks]-breaks,{1,Length[breaks]-1}]]];
Print[Style[nation<>": period(データ期間)="
<>ToString[period]<>"日間,"<>" breaks="<>ToString[breaks], Larger,Blue]]

2．ConfirmedCases の対数グラフの傾き (Subscript[m, c]) の値を求める．

前章の繰り返しになるが， 「移動平均+差分を使う方法」と「全体のFit関数の微分(oneFit)」と
「breaksに従って期間を分割してFitを使う方法(multifit)」を比較する.

casesTL=TimeSeries[cases,{dates}]//Quiet;
GraphicsRow[{DateListPlot[casesTL,PlotLabel->nation<>": cases \n(recが50人以上になった日から)",PlotRange->All],
DateListLogPlot[casesTL,PlotLabel->nation<>": cases (対数グラフ)\n(recが50人以上になった日を第1日とする)",PlotRange->All]},ImageSize->700]

logCmultifit[t_]=multifit[Log[cases],breaks];(*Log(Case)の折れ線近似*)
mc[t_]=multislope[Log[cases],breaks];(*(Subscript[m, c])の関数．multiFit*)
mr[t_]=multislope[Log[rec],breaks];(*(Subscript[m, r])の関数．multiFit*)
mcByOneFit=D[Fit[Log[cases],Table[t^k,{k,0,10}],t],t];(*全体にFitを使用し微分-onefit*)
mcByMovingAverage=Differences[MovingAverage[Log[cases],5]];(*(Subscript[m, c])の5日の移動平均*)

GraphicsGrid[{{Show[{ListLinePlot[Log[cases],PlotStyle->{Thin,Blue}],Plot[logCmultifit[t],{t,1,period},PlotStyle->Red]},PlotLabel->"Log[cases]."<>"青:data,赤:multifit(平均区間幅"<>ToString[mhaba]<>")"],ListLinePlot[{Log[cases],MovingAverage[Log[cases],5],MovingAverage[Log[cases],10]},PlotLabel->"Log[cases]の移動平均(0,5,10日)"]},
{Plot[{Tooltip[mc[t],"multifit"],Tooltip[Evaluate[mcByOneFit],"onefit"]},{t,1,period},PlotLabel->"Log[cases]の傾き (Subscript[m, c]):"<>"　multifit(平均区間幅"<>ToString[mhaba]<>")と oneFit"],
ListLinePlot[mcByMovingAverage,PlotLabel->"Log[cases]の傾き (Subscript[m, c]) : 5日の移動平均と差分"]}},ImageSize->700]

3つのグラフの傾向は同じだが，移動平均の方は,やや変動が大きく,かつ変域が左にずれる.　また onefit は両端が
大きく外にずれていく．  (Subscript[R, t])=((c*(Subscript[m, c]))/(i*γ)) を使うと (Subscript[R, t]) も同じく両端が外にずれることになる. onefitは表示には使えないが，
微分方程式に使うことは出来る.　両方に使うため multifit を使う．即ち  (Subscript[m, c])=mc[t],  (Subscript[m, r])=mr[t] (前章) とする．

3. (Subscript[m, c]) と γ から 再生産数 (Subscript[R, t]) を求める. (Subscript[R, t])<1 のときは収束,  (Subscript[R, t])>1 のときは拡大

(＃3), (#4)より γ=(Subscript[m, r])*(r/i) ,    (Subscript[R, t]) =((Sβ)/γ)=((c*(Subscript[m, c]))/(i*γ)) =((c*(Subscript[m, c]))/(r*(Subscript[m, r])))  (=(dc/dr))

γが固定された時の (Subscript[R, t]) の式

(Subscript[m, c])=mc[t] と決めたが，(Subscript[R, t]) の計算方法は  (Subscript[R, t])= ((c*(Subscript[m, c]))/(i*γ)) , (Subscript[R, t]) = ((c*(Subscript[m, c]))/(r*(Subscript[m, r])))  の2つある．数学的には同じだが，実際には
データの変動により， 2つの式の値は大きく違う．  ここでは基本方針として，

　　　　「できるだけ変動が少なく，かつその (Subscript[R, t]) で計算した cases が　データのc と一致する」

ように選ぶ．まず変動の少ないのは c/i と γ を平均値で考え，(Subscript[R, t])= ((c*(Subscript[m, c]))/(i*γ))　を使った式だろう．

cOveri=TrimmedMean[N[cases/inf],{0.2,0.2}];
γ=γ0=TrimmedMean[rec/inf*Table[mr[t],{t,1,period}],{0.2,0.2}]//N;
Rt[t_]=mc[t]*cOveri/γ;
RtTL=TimeSeries[Table[Rt[t],{t,1,period}],{dates}]//Quiet;
constTL=TimeSeries[1^Range[period],{dates}]//Quiet;

Show[{DateListPlot[{RtTL,constTL},PlotStyle->{Line,{Gray,Dashed}},
PlotLabel->nation<>": (Subscript[R, t]) by c*(Subscript[m, c])/(i*γ) ,γ="<>ToString[γ]]},ImageSize->300]

嬉しいことに 最近では (Subscript[R, t])<1と なっているようだ. 実際には γ　の値はもっと大きいと思われる(前章)．
そして, (Subscript[R, t])= ((c*(Subscript[m, c]))/(i*γ)) より，(Subscript[R, t])の値はもっと下がるはず．
(しかしc を変えないでγ のみを大きくするとiは減るので反比例ではない.でもこの値より小さくなることは間違いない)．

4.微分方程式の解と比較することにより， (Subscript[R, t]) と γ を求める

4-1.  (Subscript[R, t])=Rt[t] として良いか?

先のRt(=mc[t]*c/(i*γ)) は　「丸め」や「近似」を繰り返して得た式である．正しい保証はない．
前の sub chapter で述べた基本方針に従って「この Rt を微分方程式(#2)に使ったとき
得られる結果式の c(t),i(t),r(t) がデータに合っているかどうか」を調べる．以下 γ として
データから得られた平均値を利用する．

cOveri=TrimmedMean[N[cases/inf],{0.2,0.2}];
γ=γ0=Round[TrimmedMean[rec/inf*Table[mr[t],{t,1,period}],{0.2,0.2}],0.001];
Rt[t_]=mc[t]*cOveri/γ0;(*以上3つは直前のsub chapter と同じ*)
u=0.98;
R0[t_]:=u Rt[t];
cByData=Tooltip[cases,"cases by Data"];
dByData=Tooltip[deaths,"deathes by Data"];
iByData=Tooltip[inf,"inf by Data"];
rByData=Tooltip[rec,"rec by Data"];

Module[{sol},
sol=NDSolve[(*←微分方程式*)
{i'[t]==R0[t]i[t]*γ0-γ0*i[t], r'[t]==γ0*i[t],
i[1]==inf[[1]],r[1]==rec[[1]]},{i,r},{t,1,period}];
Show[
Plot[{Evaluate[i[t]/.sol]+Evaluate[r[t]/.sol],
Evaluate[i[t]/.sol],
Evaluate[r[t]/.sol]},{t,1,period}, PlotStyle-> {{Purple,Dashed},{Red,Dashed},{Gray,Dashed}},PlotLegends->" Line-微分方程式の解　\n",PlotRange->All],
ListLinePlot[{cByData,iByData,rByData},
PlotStyle-> {Purple,Red,Gray},PlotLegends->{"confirmed cases(感染者積算)","infected=Active Case(感染者)","recovered (回復者,死者含む)"}]]]

上のプログラムで見るように, (Subscript[R, t])=0.98*Rt[t]  とすると，cについてはほぼ完全に一致する．
i,rについてはかなりずれているが，これは 計算されたγ の値が　本当はもっと大きいからだと
推定している．(前章と次のsub chapter参照)

従って 日本については (Subscript[R, t]) =Rt[t] として良い．そして実は世界の殆どの国に対し
「Rt[t]を使って解いた方程式の解」と「データ」は, c についてはほぼ一致する.
(この章の最後のManipulateで確認できます.)

4-2. γ と (Subscript[R, t]) の正しい値を求める

(この節では 回復者は死亡者を除く．ただし recやr は 死亡者も含むとする)
上の点線と実線のグラフを比較すると r のグラフはやや不自然である．また
γ= (Subscript[γ, 0]) =0.021(5/19現在)というのは小さすぎる (回復期間=50日は長すぎる).　かつ
回復者/死亡者=6.5なのに，死亡率がたった2%というのもおかしい．(13%のはず?)
それは本当の回復者はもっと多いからではないか? (前章)

前章より 「cases=C*(Subscript[a, t])」 とざっくり近似したとすると a=1.06

a=Module[{m,n,t,sol},sol=FindFit[Log[cases],m t+n,{m,n},t];Exp[m]/.sol];
Print[Style["概算で cases=c (a^t), 但し a="<>ToString[a],Blue,Larger]]

平均の回復日数をp(=1/γ), 死亡者が死亡するまでの平均期間をq日,死亡率を2%とすると 
                          (a^(p-q))= (0.98/0.02)* 死者/回復者=49 * 死者/回復者
故に                        　　　　   p-q= (Subscript[log, a]) (49 * 死者/回復者)
                           　　　∴      q=1/γ- (Subscript[log, a]) (49 * 死者/回復者) 　　···(#5)
回復者/死亡者=6.5 を代入して計算すると,  p=48日,  q=16日 (5/19現在) となり　おかしい．
しかし回復者の率が大きくなるとp， qは小さくなる．

以上も考慮し　ここでは 「 Confirmed Cases(感染者積算)とdeaths(死亡者積算)の2つは正しいと仮定し，
γ のみを動かして，p, qの値と rec,inf さらに Rt の推移を見てみたい．」

(#1),(#3)　より                      (dc/dt)= c*(Subscript[m, c])=　β s*i  =(Subscript[R, t])*γ*i 
Confirmed cases は正しいと仮定するので　
　　　　　　　　　　　　(dc/dt)= c*(Subscript[m, c]) で方程式を立て，その解の c[t],i[t]　を使って     (Subscript[R, t])= c*mc[t]/ i/γ 
とする．(先程の Rt=mc[t]*cOveri/γ　と比較すると,共に, c/i を平均値に取るか，時間tの関数と考えるかの違い)

ClearAll[γ]

DynamicModule[{i,c,r,q,sol,rt1Table,rt2Table,cAssoc, cTable, γTable, dOverr},
γTable=Table[γ,{γ,γ0,0.1,0.001}];
cTable=Quiet[Table[
sol=NDSolve[(*微分方程式*)
{i'[t]==0.985(i[t]+r[t])mc[t]-γ*i[t], r'[t]==γ*i[t],
i[1]==inf[[1]],r[1]==rec[[1]]},{i,r},{t,1,period}];
{Evaluate[i[t]/.sol]+Evaluate[r[t]/.sol],
Evaluate[i[t]/.sol],Evaluate[r[t]/.sol]},{γ,γ0,0.1,0.001}]];
cAssoc=AssociationThread[γTable->cTable];
Manipulate[
c=cAssoc[γ][[1]];i=cAssoc[γ][[2]];r=cAssoc[γ][[3]];
dOverr=Transpose[Table[deaths[[t]]/r,{t,1,period}]];
q=Round[1/γ-First[Log[a,49*TrimmedMean[Table[deaths[[t]]/(r-deaths[[t]]),{t,1,period}],{0.2,0.2}]]]];
rt1Table=Table[Rt[t],{t,1,period}];(*γ=(Subscript[γ0に対するR, t])*)
rt2Table=Transpose[Table[c*mc[t]/i/γ,{t,1,period}]]//First;(*Manipulating γ (Subscript[に対するR, t])*)
GraphicsGrid[{
{
Show[Plot[Evaluate[cAssoc[γ]],{t,1,period},
PlotStyle->{{Purple,Dashed},{Red,Dashed},{Gray,Dashed},{Black,Dashed}},PlotRange->All],
ListLinePlot[{cByData,iByData,rByData, dByData},PlotStyle->{Purple,Red,Gray,Black}],PlotLabel->nation<>": 点線-統計データ, 実線-微分方程式の解,\n 紫:cases, 赤:inf, 灰色:rec, 黒:deaths"],

ListLinePlot[ dOverr,PlotRange->{0,0.2},
PlotLabel->" 青:deaths/rec のグラフ\n現在のγでは,死亡者の平均入院日数 q= "<>ToString[q]
<>"日\n 回復者の平均日数p(1/γ)="<>ToString[Round[1/γ]]<>"日"]
},(*ここまでRow1*)
{
ListLinePlot[ {rt1Table,rt2Table,Table[1,{period}]},
PlotRange->{{0,period},{0,8}},PlotStyle->{Red,Lighter[Purple],{Gray,Dashed}},
PlotLabel->"Rt: 赤:γ="<>ToString[γ0]<>",Rt[t]利用, 青:Manipulating γ で計算 "]
}(*ここまでRow2*)
},ImageSize->700],
{γ,γ0,0.1,0.001,Appearance->"Labeled"},
SaveDefinitions->True,FrameLabel->"Decide γ and (Subscript[R, t])"]]

［結論］ まず (Subscript[R, t])= c*mc[t]/ i/γ...(*)  と置いて微分方程式を解くと， 
c/iを平均値で近似した Rt[t] の式を使った時より c[t] とcはさらに一致している．(左上図)
次に左下図を見ると γ=(Subscript[γ0における2つのR, t]) のグラフはほぼ一致している．
即ち (Subscript[R, t])= c*mc[t]/ i/γ  において c/i は平均値を使っても時間tの式で計算しても大差ない!
（r1Tableとr2Tableの定義を参照．この結果は少し驚きであるが...)

次にγの妥当な値を決定する．これは「総合的判断」が必要になる．また現在日時にも依る. 
5/19現在, 死亡者の平均入院日数を見ると γ＝0.027ぐらいが妥当だと思う(9日間)．
しかしこれでは 5/19現在　「データのrec(統計rと書く)」>「解曲線のrec(計算rと書く)」 となりおかしい．
(5/13現在では 入院日数からγ=0.033とすると，「統計rec<計算rec」 だったので 丁度良かった.)
そして 5/19現在で「統計rec=計算rec」 とすると　γ=0.05前後となる．
しかしこの時　死亡者の平均入院日数は(-1)日となってしまう．

この理由は最近のcのカーブがなだらかに減少(aが減少)しているからだと思われる．
故にcのカーブが「ざっくり見ても指数関数」とは言えなくなる．
先の式は「ざっくり見て指数関数」というのに依存しているので，カーブが平らになると
式を改良しないと，この方法では適切なγを求めることはできないだろう．
(そんなに難しくないと思うので，時間があればやってみたい.)

たった一週間でこんなに変わって驚きだが,　5/19現在では  γ=0.05 としたい.
そうすると現時点における 計算r と 統計r が一致する．そして過去に「計算r > 統計r」だったのは
適切に統計されていなかったためと考えることが出来る．

以上から　5/19現在に置いて 妥当なγは　γ=0.05　としたい．そして そのγ　と方程式の解に基づいて計算した
　　　　　　　　　　　　　　　   (Subscript[R, t])= c*mc[t]/ i/γ 　
を採用したい． (c/i は時間平均でも，瞬間値でも(Subscript[R, t])の値は大差ない)

(Subscript[R, t])は 残念であるが，γ がほぼ2倍になっても　反比例して1/2倍にはならない．
(Subscript[R, t])>>1のときはかなり変化が大きいが, (Subscript[R, t])≤1のときは，あまり変化しない.

＊なお私が γを決めた時の Manipulate の実行画面を下に置いておきます.

Hyperlink["5月13日の実行画面","http://mixedmoss.com/miscellaneous/MathematicaSIRModel/decideGamma.png",ImageSize->700]
Hyperlink["5月19日の実行画面","http://mixedmoss.com/miscellaneous/MathematicaSIRModel/decideGamma2.png",ImageSize->700]

5. (Subscript[R, t]) の表示関数の作成

いろいろ分析したので長くなりましたが，以上から 統計データを完全に信ずるなら 
(γ=γ0=0.021(5/19現在)で良いなら) (Subscript[R, t ])を求めるのは簡単です．
r の統計データを信じないのなら(γ=γ0は良くないと思うなら) ，
微分方程式を解かないと (Subscript[R, t]) は求まりません．でもMathematica では比較的簡単に求められます． 
なお，使うには「zip,slopes,multislope 」(notebook最上部)が必要です．

以下のプログラムは 
     (1)casesが正しいと仮定して，mc[t]を求める． 
　(2)  cases を動かさずに  γを与えて(データから計算されたγの値は使わないで)微分方程式を解いて r,i,c を求める. 
(3)  (Subscript[R, t])= c*mc[t]/ i/γ  で (Subscript[R, t]) を求める.( c/i　は平均値ではなく,時間tの関数).

の手順で作成しています．最後の引数は ラベル(PlotLabel)です.必要ない時は入力不要です.
casesやrecやdates は整形前でも大丈夫です．また見やすいように最後に Fit で平滑化しています．

あまり複雑な関数を作ったことがないので使いづらいと思います．プログラムを見て適当に書き直して下さい.
(こんな作り方をしてたら，何種類も作らないといけないので，なにか高度な良い方法があるはず.)

入力が 「γ , cases ,rec, dates (,Label) 」 のとき， (Subscript[R, t]) の表示関数の作成．

showRt[gamma_,cases0_List,rec0_List,dates0_List,label_String:""]:=Module[
{day1,dayl, dates, cases,rec,inf0, dayInfLess10, cc,ii, haba,period,breaks,mc,sol,rtTable,rtTL,constTL,
rtTableFit},
day1=First[FirstPosition[ rec0-10,_?Positive]];(*rec>10になってから解析*)
inf0=cases0-rec0;
dayInfLess10=Flatten[Position[ inf0-10,_?Negative]];
dayl=Min[Select[dayInfLess10,#>day1&]-1,Length[cases0]];(*inf<10になると中止*)
cases=Take[ cases0,{day1,dayl}];
rec=Take[ rec0,{day1,dayl}];
dates=Take[ dates0,{day1,dayl}];
period=Length[cases];haba=5;(*期間分割幅*)
breaks=Table[1+haba k,{k,0,(period-1)/haba-1}]//Append[period];
mc[t_]=multislope[Log[cases],breaks];
sol=NDSolve[(*微分方程式*)
{i'[t]==0.99(i[t]+r[t])mc[t]-gamma*i[t], r'[t]==gamma*i[t],
i[1]==Round[cases[[1]]/cOveri],r[1]==Round[cases[[1]](1-1/cOveri)]},{i,r},{t,1,period}];
cc=Evaluate[i[t]/.sol]+Evaluate[r[t]/.sol];
ii=Evaluate[i[t]/.sol];
rtTable=First[Transpose[Table[mc[t]cc/ii/gamma,{t,1,period}]]];
(*rtTL=TimeSeries[rtTable,{dates}]//Quiet;*)(*平滑化していない*)
rtTableFit=Fit[rtTable,Table[t^k,{k,1,Floor[period/4]}],t]//Evaluate;
rtTL=TimeSeries[Table[rtTableFit,{t,1,period}],{dates}]//Quiet;(*平滑化*)
constTL=TimeSeries[Table[1,{period}],{dates}]//Quiet;
DateListPlot[{rtTL,constTL},PlotStyle->{Line,Dashed},Epilog->{Lighter[Purple],PointSize[Small],Tooltip[Point[{Last[dates],Last[rtTable]}],ToString[Round[Last[rtTable],0.1]]],White,PointSize[Medium],Tooltip[Point[{Last[dates],0}],ToString[Last[dates]]]},PlotLabel->label]]

6. 東京,北海道,大阪と日本全体の (Subscript[R, t]) を図示

表示関数(showLocalRt)の作成

東京のデータは荻原さんがGitHubに公開されておられる素データをDownload します．
(https://github.com/kaz-ogiwara/covid19/blob/master/data/prefectures.csv)
こちらは3月11日からのデータになります. 
先程の showRt を使って γと地方名を入力すると, (Subscript[R, t])が描けるような関数　showLocalRt を作ります．

japanData=Import["https://raw.githubusercontent.com/kaz-ogiwara/covid19/master/data/prefectures.csv"];
japanDates=Select[japanData,Part[#,5]=="Tokyo"&][[All,1;;3]];

showLocalRt[gamma_,region_String]:=Module[
{localData,localCases,localDischarged,localDeaths,localRec,localInf},
localData=Select[japanData,Part[#,5]==region&];
localCases=localData[[All,6]]/.{""->0};
localDischarged=localData[[All,8]]/.{""->0};
localDeaths=localData[[All,9]]/.{""->0};
localRec=localDischarged+localDeaths;
localInf=localCases-localRec;
showRt[gamma,localCases,localRec,japanDates,"(Subscript[R, t]):"<>region<>" : (Subscript[γ, input])="<>ToString[gamma]]]

東京,北海道,大阪と日本全体の(Subscript[R, t]) の比較

γ = 0.05 で計算します．また最近の R の値を Tooltip　しました．

γ1=0.05

GraphicsGrid[{{showLocalRt[γ1,"Tokyo"],showLocalRt[γ1,"Osaka"]},
{showLocalRt[γ1,"Hokkaido"],showRt[γ1,cases,rec,dates,"(Subscript[R, t]):"<>nation<>" : (Subscript[γ, input])="<>ToString[γ1]]}}]

3月11日からなので,北海道の2月のピークは見えません．傾向は意外なほど違いますが， 
どちらも最近は余裕で R < 1 となっているようです．

7. 世界の (Subscript[R, t]) を表示 (with Manipulate)

(Subscript[R, t]) の表示関数:showGlobalRt を作成　(入力は γ 　と　国名 )

showGlobalRt[γ,国名]で，Rtのグラフのみ, 
showGlobalRt[1,γ,国名]で，Epidemic curve も描きます. 
showGlobalRt[2,γ,国名]で，(Subscript[γ, input]) により計算された微分方程式の解も重ねて図示します. 
showGlobalRt[3...]とshowGlobalRt[4...]は，mode1,mode2の対数グラフです. 
さらに mode が 5～8の時は，Rtの代わりに, データから計算したγのグラフを描きます.

r>50となった日(day1)から計算を始め, 「i<10または最新データ獲得の日」(dayl)まで (Subscript[R, t]) を計算します．
従って一度完全に収束して(i<10になって)，その後第2波が来た時は，計算されません. 

見やすくするために, (Subscript[R, t ])の グラフは Fit を使って平滑化しています.
showRt を利用して作ると悲しいことに遅くなったので，最初から作ってます．

入力されたγの値(Subscript[γ, input]) と, データから計算した (OverBar[(Subscript[γ, t])]) の値を表示します．
故に (Subscript[γ, input])=   (OverBar[(Subscript[γ, t])])　 とすると，データ曲線と解曲線は近くなります．このとき γ の値が妥当なものならOKですが，
そうでない時はどこかが間違っていることになります．

casesGlobal0=Import["https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv"];
recGlobal0=Import["https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv"];
deathsGlobal0=Import["https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv"];

datesGlobal=Drop[casesGlobal0[[1]],4];
casesWorld=Join[{ "","World","","" },Plus@@(Drop[#,4]&)/@ Drop[casesGlobal0,1]];
recWorld=Join[{ "","World","","" },Plus@@(Drop[#,4]&)/@ Drop[recGlobal0,1]];
deathsWorld=Join[{ "","World","","" },Plus@@(Drop[#,4]&)/@ Drop[deathsGlobal0,1]];

Clear[p,q]
rule={{"Hong Kong","China",q__}->{" ","Hong Kong",q},{p_,"China",q__}->{p,"China(Others)",q},{"Hubei","China(Others)",q__}->{"Hubei","China(Hubei)",q}};
casesGlobalChinaModified=casesGlobal0//.rule;
recGlobalChinaModified=recGlobal0//.rule;
deathsGlobalChinaModified=deathsGlobal0//.rule;

casesGlobal=Append[casesGlobalChinaModified,casesWorld];
recGlobal=Append[recGlobalChinaModified,recWorld];
deathsGlobal=Append[deathsGlobalChinaModified,deathsWorld];

countriesAll=Union[Drop[casesGlobalChinaModified[[All,2]],1]];
cOver1000Q[x_]:=Last[Flatten[Select[casesGlobal,Part[#,2]==x&]]]≥1000;
rOver200Q[x_]:=Last[Flatten[Select[recGlobal,Part[#,2]==x&]]]≥200;
countries=Append[Select[countriesAll,(cOver1000Q[#]&&rOver200Q[#]&&#≠"Tajikistan")&],"World"];
(*Tajikistanのデータはおかしいので除く*)


showGlobalRt[mode_Integer:0,gamma_,nation_String]:=Module[
{day1,dayl,dates0,dates,dayInfLess50,cases0,deaths0,cases,rec,deaths,rec0,inf0,inf,haba,period,breaks,mr,mc,sol,γculculated,
rtTable,dayInfLess10,rtFit,γFit,γTL,γCurve,csolTL,isolTL,rsolTL,constTL,rtTL,rtCurve,caseTL,recTL,deathsTL,infTL,myPlot,plotmode,cir,cirsol,cc,ii,rr},
dates0=datesGlobal;
cases0=Plus@@(Drop[#,4]&)/@Select[casesGlobal,Part[#,2]==nation&];
rec0=Plus@@(Drop[#,4]&)/@Select[recGlobal,Part[#,2]==nation&]+Plus@@(Drop[#,4]&)/@Select[deathsGlobal,Part[#,2]==nation&];
deaths0=Plus@@(Drop[#,4]&)/@Select[deathsGlobal,Part[#,2]==nation&];
inf0=cases0-rec0;
day1=First[FirstPosition[ rec0-30,_?Positive]];(*rec>30になってから解析*)
dayInfLess10=Flatten[Position[ inf0-10,_?Negative]];
dayl=Min[Select[dayInfLess10,#>day1&]-1,Length[cases0]];(*inf<10になると中止*)
dates=Take[ dates0,{day1,dayl}];
cases=Take[ cases0,{day1,dayl}];(*day1>=dayl となると問題!*)
deaths=Take[ deaths0,{day1,dayl}];
rec=Take[ rec0,{day1,dayl}];
inf=cases-rec;
period=Length[cases];haba=5;(*期間分割幅*)
breaks=Table[1+haba k,{k,0,(period-1)/haba-1}]//Append[period];

mc[t_]=multislope[Log[cases],breaks];
mr[t_]=multislope[Log[rec],breaks];
γculculated=Round[TrimmedMean[rec/inf*Table[mr[t],{t,1,period}],{0.1,0.1}],0.0001];
sol=NDSolve[
{i'[t]==(i[t]+r[t])mc[t]-gamma*i[t], r'[t]==gamma*i[t],
i[1]==inf[[1]],r[1]==rec[[1]]},{i,r},{t,1,period}];
ii=Evaluate[i[t]/.sol];
rr=Evaluate[r[t]/.sol];
cc=ii+rr;

(*rtTL:=Tooltip[TimeSeries[zip[dates,rtTable]],"(Subscript[R, t])"]//Quiet;平滑化していない曲線*)
rtTable:=First[Transpose[Table[mc[t]*cc/ii/gamma,{t,1,period}]]];(*Rt*)
rtFit:=Fit[rtTable,Table[t^k,{k,1,Floor[period/4]}],t]//Evaluate;(*Rt*)
rtTL:=Tooltip[TimeSeries[Table[rtFit,{t,1,period}],{dates}],"(Subscript[R, t])"]//Quiet;(*Rt*)
rtCurve:=DateListPlot[{rtTL,constTL},PlotStyle->{Line,Dashed},
PlotLabel->nation<>": (Subscript[R, t])"<>", (Subscript[γ, input])="<>ToString[gamma]<>",  (OverBar[(Subscript[γ, t])])="<>ToString[γculculated],
Epilog->{Lighter[Purple],PointSize[Small],Tooltip[Point[{Last[dates],Last[rtTable]}],ToString[Round[Last[rtTable],0.1]]],
White,PointSize[Medium],Tooltip[Point[{Last[dates],0}],ToString[Last[dates]]]},ImageSize->Medium];

γFit:=Fit[rec/inf*Table[mr[t],{t,1,period}],Table[t^k,{k,1,Floor[period/4]}],t]//Evaluate;(*γfit*)
γTL:=Tooltip[TimeSeries[Table[γFit,{t,1period}],{dates}],"(Subscript[γ, t]) from data"]//Quiet;(*γfit曲線*)
γCurve:=DateListPlot[γTL,PlotRange->All,PlotLabel->nation<>": (Subscript[γ, t]) from data, (OverBar[(Subscript[γ, t])])="<>ToString[γculculated]];(*γfit曲線*)

csolTL:=Tooltip[TimeSeries[Table[cc,{t,1,period}],{dates}],"感染者積算(Confirmed Cases)"]//Quiet;(*解曲線*)
isolTL:=Tooltip[TimeSeries[Table[ii,{t,1,period}],{dates}],"感染者(Active cases)"]//Quiet;(*解曲線*)
rsolTL:=Tooltip[TimeSeries[Table[rr,{t,1,period}],{dates}],"回復者(死者含む)"]//Quiet;(*解曲線*)
cirsol:=myPlot[{csolTL,isolTL,rsolTL},
PlotStyle->{{Purple,Dashed},{Red,Dashed},{Gray,Dashed}}];

caseTL=Tooltip[TimeSeries[cases,{dates}],"感染者積算(Confirmed Cases)"]//Quiet;(*データ曲線*)
infTL=Tooltip[TimeSeries[inf,{dates}],"感染者(Active cases)"]//Quiet;(*データ曲線*)
recTL=Tooltip[TimeSeries[rec,{dates}],"回復者(死者含む)"]//Quiet;(*データ曲線*)
constTL=TimeSeries[Table[1,{period}],{dates}]//Quiet;
cir:=myPlot[{caseTL,infTL,recTL},
,PlotStyle->{Purple,Red,Gray},
Epilog->{Purple,PointSize[Small],Tooltip[Point[{Last[dates],Last[cases]}],ToString[Last[cases]]],
Red,Tooltip[Point[{Last[dates],Last[inf]}],ToString[Last[inf]]],
Gray,Tooltip[Point[{Last[dates],Last[rec]}],ToString[Last[rec]]],
White,PointSize[Medium],Tooltip[Point[{Last[dates],0}],ToString[Last[dates]]]
},
PlotRange->All,PlotLabel->"["<>plotmode<>"] 直線:データ, 点線:ODEの解.  ☑ToolTip"
];

Which[
mode==0,rtCurve,
mode==1,myPlot=DateListPlot;plotmode="Linear"; GraphicsRow[{ cir,rtCurve}],
mode==2,myPlot=DateListPlot;plotmode="Linear";GraphicsRow[{Show[{cir,cirsol}],rtCurve}],
mode==3,myPlot=DateListLogPlot;plotmode="Log"; GraphicsRow[{ cir,rtCurve}],
mode==4,myPlot=DateListLogPlot;plotmode="Log";GraphicsRow[{Show[{cir,cirsol}],rtCurve}],
mode==5,myPlot=DateListPlot;plotmode="Linear"; GraphicsRow[{ cir,γCurve}],
mode==6,myPlot=DateListPlot;plotmode="Linear"; GraphicsRow[{Show[{cir,cirsol}],γCurve}],
mode==7,myPlot=DateListLogPlot;plotmode="Log"; GraphicsRow[{ cir,γCurve}],
mode==8,myPlot=DateListLogPlot;plotmode="Log";GraphicsRow[{Show[{cir,cirsol}],γCurve}]
]
]

US, Italy, Brazil, Sweden の (Subscript[R, t ])を見る

γ1=0.05;
GraphicsGrid[{{
showGlobalRt[γ1,"US"],
showGlobalRt[γ1,"Italy"]},
{showGlobalRt[γ1,"Brazil"],
showGlobalRt[γ1,"Sweden"]}},ImageSize->600]

USはいま一歩.  Italyは 余裕で1 より小さいです.   逆にBrazilは「2.2前後」，
Swedenは 徐々に下がってきているもののまだ1より大きいです．

ただこれらは γ=0.05 という前提です．γの値は世界中でそんなに変わらないはずですが，
多少は医療設備や統計のとり方により変わると思います．
デー(Subscript[タから単純に計算したγの値がγ, cul]) ですが，国によって本当にまちまちです．
やはり私は 統計データのrecはあまり信頼できないと思います．
同じであるべき γ が国によってこんなに違うからです.

Italyを見ると (Subscript[R, t ]) の値が最近30日間でも変わっているのがお分かりと思います．(SpainやGermanyなども同じ) 
SIR理論に依ると β,γ は時間に依らないので 初期の段階では(s/人口=1 に近い間は) (Subscript[R, t ])は定数のはずです．
もちろん それは何も対策を講じなかった時の話で実際には各国で対策をするのでそれによって (Subscript[R, t]) は変わります．
ブラジルは何も対策してないようなので最近50日間は (Subscript[R, t]) は一定です．ブラジルはSIR理論に忠実です.
しかしイタリアは3月初旬からロックダウンをしているので，2~3週間は徐々に (Subscript[R, t ])が直線的に下がるのは理論通りです．
しかしその後も緩やかに下がり続けているのは何故でしょう？　
スウェーデンは何も対策してないはずなのに何故下がっているのでしょう？
素人考えですが，これは SIR理論だけでは説明できないと思います．

(Subscript[R, t]) (& epidemic curve) の　Manipulate

Confirmed cases≥1000かつ回復者≥200 となっている国と, 全世界 の (Subscript[R, t]) を表示できます．
(但しカジキスタンはデータが少しおかしいので除きました.)
香港は中国の一部になっていましたが，独立に扱えるように変えました．
Hubei省と他の省も分けました．なお台湾はサンプルが小さいので(Subscript[R, t])　は計算できませんでした．
さらに 全世界(World)も付け加えました．
dayl(daylast) の s,i,r,Rtの値は Tooltipで見ることができます．

showRt は,中で微分方程式を解いているので重いです．
重すぎて γ をスライダーで動かす事はできないので PopupMenu で動かします．
さらにGitHub(Johns Hopkins 大学)から データを Download しているので
スタート時は余計に時間がかかります．また外部からDownloadするので
 警告メッセージ（動的機能を有効化) が出るかも知れません．
[cdf playerは遅いので 特にスタート時は時間がかかります．10~15秒ぐらいお待ち下さい!]

Manipulate[
showGlobalRt[1+Boole[!noODE]+2Boole[!linear]+4Boole[!hideγ],γ,country]
,{{linear,True,Style["描画モード",Smaller]},{True->"Linear",False->"Log"},ControlType->RadioButton},
{{noODE,True,Style["微分方程式の解",Smaller]},{True->"Hide",False->"Show"},ControlType->RadioButton},
{{hideγ,True,Style["右側の曲線",Smaller]},{True->"(Subscript[R, t]) ",False->"(Subscript[γ, t])"},ControlType->RadioButton},
{{γ,0.05},Range[100]*0.001,ControlType->PopupMenu},
{{country,"Japan"},countries,ControlType->PopupMenu},
FrameLabel->"Epidemic and (Subscript[R, t])/γ curves with Mathematica "
, SaveDefinitions->True
]