This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License

using CovidSEIR,  Plots, DataFrames, JLD2, StatsPlots, MCMCChains
Plots.pyplot()
jmddir = normpath(joinpath(dirname(Base.find_package("CovidSEIR")),"..","docs","jmd"))
covdf = covidjhudata();

China

using Dates
dayt0 = Dates.Date("2020-01-22") - Dates.Day(65)
china = CountryData(covdf, "China", 65)
CovidSEIR.CountryData{Float64,Int64}(1.39273e9, [65, 66, 67, 68, 69, 70, 71
, 72, 73, 74  …  141, 142, 143, 144, 145, 146, 147, 148, 149, 150], [17.0, 
18.0, 26.0, 42.0, 56.0, 82.0, 131.0, 133.0, 171.0, 213.0  …  3335.0, 3337.0
, 3339.0, 3340.0, 3343.0, 3343.0, 3345.0, 3345.0, 3346.0, 3346.0], [28.0, 3
0.0, 36.0, 39.0, 49.0, 58.0, 101.0, 120.0, 135.0, 214.0  …  77410.0, 77567.
0, 77679.0, 77791.0, 77877.0, 77956.0, 78039.0, 78200.0, 78311.0, 78401.0],
 [503.0, 595.0, 858.0, 1325.0, 1970.0, 2737.0, 5277.0, 5834.0, 7835.0, 9375
.0  …  1973.0, 1905.0, 1865.0, 1810.0, 1794.0, 1835.0, 1829.0, 1761.0, 1699
.0, 1656.0])
using Turing
mdl = CovidSEIR.TimeVarying.countrymodel(china)
cc = Turing.psample(mdl, NUTS(0.65), 5000, 4)
import JLD2
JLD2.@save "$jmddir/china_$(Dates.today()).jld2" cc dayt0
JLD2.@load "$jmddir/china_dhmc_2020-04-13.jld2" cc dayt0;
cc = MCMCChains.Chains(collect(cc.value.data), replace.(cc.name_map.parameters, r"([^\[])([1-9])" => s"\1[\2]"))
Object of type Chains, with data of type 5000×15×4 Array{Float64,3}

Iterations        = 1:5000
Thinning interval = 1
Chains            = 1, 2, 3, 4
Samples per chain = 5000
parameters        = τ, sigD, sigC, sigRc, a, pE0, p[1], p[2], β[1], β[2], β
[3], γ[1], γ[2], ρ[1], ρ[2]

2-element Array{MCMCChains.ChainDataFrame,1}

Summary Statistics
  parameters        mean        std  naive_se      mcse       ess    r_hat
  ──────────  ──────────  ─────────  ────────  ────────  ────────  ───────
           τ      0.0000     0.0000    0.0000    0.0000   80.3213   1.5862
        sigD     90.1865   108.1679    0.7649    7.6652   80.3213  16.4095
        sigC  12749.7844  4659.6012   32.9484  322.0747   80.3213   4.5433
       sigRc   8238.5858  3466.7500   24.5136  242.3239   80.3213   4.9360
           a      0.7721     0.1378    0.0010    0.0069  120.2959   1.1108
         pE0      0.0000     0.0000    0.0000    0.0000   80.3213   1.6474
        p[1]      0.0004     0.0005    0.0000    0.0000   80.3213   7.1795
        p[2]      0.4376     0.4398    0.0031    0.0309   80.3213   6.8857
        β[1]      0.3384     0.3067    0.0022    0.0213   80.3213   2.0684
        β[2]      0.7464     0.6525    0.0046    0.0298  124.5414   1.1639
        β[3]      1.4304     1.4707    0.0104    0.1042   80.3213   9.1300
        γ[1]      1.2668     1.2386    0.0088    0.0878   80.3213  13.3651
        γ[2]      0.1599     0.1464    0.0010    0.0081   86.9155   1.3707
        ρ[1]      0.8703     0.1652    0.0012    0.0117   80.3213   2.2865
        ρ[2]     87.7994     1.4004    0.0099    0.0850   80.3213   1.5699

Quantiles
  parameters       2.5%       25.0%       50.0%       75.0%       97.5%
  ──────────  ─────────  ──────────  ──────────  ──────────  ──────────
           τ     0.0000      0.0000      0.0000      0.0000      0.0000
        sigD    23.6088     26.6456     28.8155     91.4797    292.9900
        sigC  4388.5359  10075.1023  14738.3582  15853.7267  17853.3371
       sigRc  4613.6964   5290.9410   6950.5976  10492.7804  14948.0661
           a     0.4875      0.6742      0.7763      0.8873      0.9793
         pE0     0.0000      0.0000      0.0000      0.0000      0.0000
        p[1]     0.0000      0.0000      0.0001      0.0007      0.0014
        p[2]     0.0021      0.0033      0.1933      0.8908      0.9903
        β[1]     0.0269      0.1482      0.2338      0.3377      1.1417
        β[2]     0.0000      0.2806      0.5279      1.1032      2.3880
        β[3]     0.0530      0.1431      0.8618      2.9038      3.8533
        γ[1]     0.0819      0.0852      0.6694      2.5336      3.0000
        γ[2]     0.0336      0.0461      0.0849      0.2376      0.5292
        ρ[1]     0.4234      0.7741      0.9711      0.9943      0.9995
        ρ[2]    85.6977     86.6899     87.6515     88.7356     90.8120

Estimates

plot(cc)

describe(cc)
2-element Array{MCMCChains.ChainDataFrame,1}

Summary Statistics
  parameters        mean        std  naive_se      mcse       ess    r_hat
  ──────────  ──────────  ─────────  ────────  ────────  ────────  ───────
           τ      0.0000     0.0000    0.0000    0.0000   80.3213   1.5862
        sigD     90.1865   108.1679    0.7649    7.6652   80.3213  16.4095
        sigC  12749.7844  4659.6012   32.9484  322.0747   80.3213   4.5433
       sigRc   8238.5858  3466.7500   24.5136  242.3239   80.3213   4.9360
           a      0.7721     0.1378    0.0010    0.0069  120.2959   1.1108
         pE0      0.0000     0.0000    0.0000    0.0000   80.3213   1.6474
        p[1]      0.0004     0.0005    0.0000    0.0000   80.3213   7.1795
        p[2]      0.4376     0.4398    0.0031    0.0309   80.3213   6.8857
        β[1]      0.3384     0.3067    0.0022    0.0213   80.3213   2.0684
        β[2]      0.7464     0.6525    0.0046    0.0298  124.5414   1.1639
        β[3]      1.4304     1.4707    0.0104    0.1042   80.3213   9.1300
        γ[1]      1.2668     1.2386    0.0088    0.0878   80.3213  13.3651
        γ[2]      0.1599     0.1464    0.0010    0.0081   86.9155   1.3707
        ρ[1]      0.8703     0.1652    0.0012    0.0117   80.3213   2.2865
        ρ[2]     87.7994     1.4004    0.0099    0.0850   80.3213   1.5699

Quantiles
  parameters       2.5%       25.0%       50.0%       75.0%       97.5%
  ──────────  ─────────  ──────────  ──────────  ──────────  ──────────
           τ     0.0000      0.0000      0.0000      0.0000      0.0000
        sigD    23.6088     26.6456     28.8155     91.4797    292.9900
        sigC  4388.5359  10075.1023  14738.3582  15853.7267  17853.3371
       sigRc  4613.6964   5290.9410   6950.5976  10492.7804  14948.0661
           a     0.4875      0.6742      0.7763      0.8873      0.9793
         pE0     0.0000      0.0000      0.0000      0.0000      0.0000
        p[1]     0.0000      0.0000      0.0001      0.0007      0.0014
        p[2]     0.0021      0.0033      0.1933      0.8908      0.9903
        β[1]     0.0269      0.1482      0.2338      0.3377      1.1417
        β[2]     0.0000      0.2806      0.5279      1.1032      2.3880
        β[3]     0.0530      0.1431      0.8618      2.9038      3.8533
        γ[1]     0.0819      0.0852      0.6694      2.5336      3.0000
        γ[2]     0.0336      0.0461      0.0849      0.2376      0.5292
        ρ[1]     0.4234      0.7741      0.9711      0.9943      0.9995
        ρ[2]    85.6977     86.6899     87.6515     88.7356     90.8120

Fit

sdf = simtrajectories(cc, china, 1:150)
f = plotvars(sdf, china, dayt0=dayt0)
plot(f.fit, xlim=nothing, ylim=(0, maximum(china.active)*2))

Implications

for fig in f.trajectories
  display(plot(fig, xlim=nothing))
end