Load the ilmtools package...


In [10]:
include("ilmtools.jl")


Out[10]:
create_loglikelihood (generic function with 1 method)

Generate a population from a bi-variate uniform distribution...


In [11]:
pop_db1 = create_pop_db_univariates(100, Uniform(0,10), Uniform(0,10))


Out[11]:
ind_idxy
116.4182617205219921.445838370414274
228.3372500394220352.8406069386986066
338.9328501065182334.110438699717795
447.0746862385106726.088751472671898
559.2153445783370968.129404114502712
669.9647880272693145.58736355229446
773.0645097059542689.907423640128933
880.00453325405911142456.989576439410259
998.525039785034375.525692719498913
10101.5922951509521214.550163217556338
11115.0218555556858415.00398335676647
12123.9545550849504029.78530450070842
13133.65381039901042654.20223096376001
14144.6768891888476390.5412347803985229
15155.2616863901586511.6615794309246823
16168.2093736013342740.032678307031108744
17177.8554228627340083.217382582744357
18185.8147402157937144.32296386218005
19197.5082205124702935.093545343779899
20206.1327084673806656.678955952504357
21215.081192752663939.760072312951289
22226.0490622255551880.6319165607887012
23231.57523671810521074.68550047540613
24249.2976553879933572.2732830409997606
25256.7415765966408455.20796370657008
26265.40637157853788255.804786658776944
27278.3089802339045992.8389545440289154
28285.8171632224346651.3702372207750746
29290.426713518814927364.400693400628919
30304.30188513894064250.18116583630979788
&vellip&vellip&vellip&vellip

Generate the initial infection and propagate infection through population (SIR model) in continuous time with parameters $\alpha=2, \beta=6, \gamma=5$


In [12]:
evdb = infect_recover_loop(pop_db1, "continuous", "SIR", 2, 6, 5)


Out[12]:
edb(100x4 DataFrame
| Row | ind_id | s | i       | r       |
|-----|--------|---|---------|---------|
| 1   | 1      | 0 | 1.41448 | 36.1383 |
| 2   | 2      | 0 | 1.00501 | 1.74011 |
| 3   | 3      | 0 | 1.19872 | 6.83014 |
| 4   | 4      | 0 | 2.32629 | 8.44914 |
| 5   | 5      | 0 | 2.56692 | 2.77236 |
| 6   | 6      | 0 | 2.98095 | 9.69036 |
| 7   | 7      | 0 | 3.70584 | 5.34881 |
| 8   | 8      | 0 | 1.88059 | 2.75702 |
| 9   | 9      | 0 | 2.09324 | 3.23807 |
| 10  | 10     | 0 | 1.84361 | 9.93704 |
| 11  | 11     | 0 | 1.57463 | 5.59969 |
⋮
| 89  | 89     | 0 | 1.00075 | 2.98217 |
| 90  | 90     | 0 | 2.2662  | 5.65972 |
| 91  | 91     | 0 | 1.8723  | 17.4001 |
| 92  | 92     | 0 | 1.8418  | 2.20715 |
| 93  | 93     | 0 | 1.86421 | 7.6723  |
| 94  | 94     | 0 | 2.33476 | 5.85296 |
| 95  | 95     | 0 | 2.28791 | 11.7352 |
| 96  | 96     | 0 | 2.88563 | 12.0168 |
| 97  | 97     | 0 | 1.84779 | 14.3681 |
| 98  | 98     | 0 | 2.26694 | 8.66313 |
| 99  | 99     | 0 | 1.8453  | 3.63309 |
| 100 | 100    | 0 | 1.19141 | 3.60301 |,[1.0,1.00075,1.00349,1.00501,1.00501,1.00701,1.00709,1.00716,1.00916,1.01099  …  17.4001,18.2024,20.0489,20.0559,20.382,24.1985,24.2646,36.1383,Inf,Inf],"continuous")

Generate the efficient log likelihood function


In [13]:
ll = create_loglikelihood(pop_db1, evdb)


Out[13]:
cSIR_ll (generic function with 1 method)

Time the log likelihood function


In [14]:
ll((1, 7, 4)), ll((3, 7, 7)), ll((1, 7, 4)), @time ll((2, 6, 5))


elapsed time: 0.092694488 seconds (3106912 bytes allocated, 59.26% gc time)
Out[14]:
(-16.88799811085775,-98.0622516950068,-16.88799811085775,14.235248127465708)

Run MCMC


In [15]:
n = 20000
burnin = 10000
sim = Chains(n, 3, names = ["alpha", "beta", "gamma"])
theta = AMMVariate([2.0, 6.0, 5])
SigmaF = cholfact(eye(3))
@time for i in 1:n
    amm!(theta, SigmaF, ll, adapt = (i <= burnin))
    sim[i,:,1] = theta[1:3]
end
Mamba.describe(sim)


elapsed time: 167.104514631 seconds (88319132372 bytes allocated, 55.09% gc time)
Iterations = 1:20000
Thinning interval = 1
Chains = 1
Samples per chain = 20000

Empirical Posterior Estimates:
4x6 Array{Any,2}:
 ""        "Mean"   "SD"      "Naive SE"   "MCSE"         "ESS"
 "alpha"  1.92849  0.216126  0.00152824   0.0051368   5950.16  
 "beta"   6.17221  0.129917  0.000918654  0.00285728  6430.28  
 "gamma"  5.65376  0.578209  0.00408855   0.0144666   5652.4   

Quantiles:
4x6 Array{Any,2}:
 ""        "2.5%"   "25.0%"   "50.0%"   "75.0%"   "97.5%"
 "alpha"  1.52376  1.77957   1.92054   2.07426   2.36573 
 "beta"   5.91418  6.08666   6.17488   6.2602    6.41556 
 "gamma"  4.63311  5.24117   5.62026   6.02057   6.89933 

Produce trace plots of $\alpha$, $\beta$, and $\gamma$


In [16]:
plot(x=burnin:n, y=sim.value[burnin:n, 1, 1], Geom.line)


Out[16]:
x -5.0×10³ 0 5.0×10³ 1.0×10⁴ 1.5×10⁴ 2.0×10⁴ 2.5×10⁴ 3.0×10⁴ 3.5×10⁴ 0 5.00×10² 1.00×10³ 1.50×10³ 2.00×10³ 2.50×10³ 3.00×10³ 3.50×10³ 4.00×10³ 4.50×10³ 5.00×10³ 5.50×10³ 6.00×10³ 6.50×10³ 7.00×10³ 7.50×10³ 8.00×10³ 8.50×10³ 9.00×10³ 9.50×10³ 1.00×10⁴ 1.05×10⁴ 1.10×10⁴ 1.15×10⁴ 1.20×10⁴ 1.25×10⁴ 1.30×10⁴ 1.35×10⁴ 1.40×10⁴ 1.45×10⁴ 1.50×10⁴ 1.55×10⁴ 1.60×10⁴ 1.65×10⁴ 1.70×10⁴ 1.75×10⁴ 1.80×10⁴ 1.85×10⁴ 1.90×10⁴ 1.95×10⁴ 2.00×10⁴ 2.05×10⁴ 2.10×10⁴ 2.15×10⁴ 2.20×10⁴ 2.25×10⁴ 2.30×10⁴ 2.35×10⁴ 2.40×10⁴ 2.45×10⁴ 2.50×10⁴ 2.55×10⁴ 2.60×10⁴ 2.65×10⁴ 2.70×10⁴ 2.75×10⁴ 2.80×10⁴ 2.85×10⁴ 2.90×10⁴ 2.95×10⁴ 3.00×10⁴ 0 1×10⁴ 2×10⁴ 3×10⁴ 0 1.0×10³ 2.0×10³ 3.0×10³ 4.0×10³ 5.0×10³ 6.0×10³ 7.0×10³ 8.0×10³ 9.0×10³ 1.0×10⁴ 1.1×10⁴ 1.2×10⁴ 1.3×10⁴ 1.4×10⁴ 1.5×10⁴ 1.6×10⁴ 1.7×10⁴ 1.8×10⁴ 1.9×10⁴ 2.0×10⁴ 2.1×10⁴ 2.2×10⁴ 2.3×10⁴ 2.4×10⁴ 2.5×10⁴ 2.6×10⁴ 2.7×10⁴ 2.8×10⁴ 2.9×10⁴ 3.0×10⁴ -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5.0 -2 0 2 4 6 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 y

In [17]:
plot(x=burnin:n, y=sim.value[burnin:n, 2, 1], Geom.line)


Out[17]:
x -5.0×10³ 0 5.0×10³ 1.0×10⁴ 1.5×10⁴ 2.0×10⁴ 2.5×10⁴ 3.0×10⁴ 3.5×10⁴ 0 5.00×10² 1.00×10³ 1.50×10³ 2.00×10³ 2.50×10³ 3.00×10³ 3.50×10³ 4.00×10³ 4.50×10³ 5.00×10³ 5.50×10³ 6.00×10³ 6.50×10³ 7.00×10³ 7.50×10³ 8.00×10³ 8.50×10³ 9.00×10³ 9.50×10³ 1.00×10⁴ 1.05×10⁴ 1.10×10⁴ 1.15×10⁴ 1.20×10⁴ 1.25×10⁴ 1.30×10⁴ 1.35×10⁴ 1.40×10⁴ 1.45×10⁴ 1.50×10⁴ 1.55×10⁴ 1.60×10⁴ 1.65×10⁴ 1.70×10⁴ 1.75×10⁴ 1.80×10⁴ 1.85×10⁴ 1.90×10⁴ 1.95×10⁴ 2.00×10⁴ 2.05×10⁴ 2.10×10⁴ 2.15×10⁴ 2.20×10⁴ 2.25×10⁴ 2.30×10⁴ 2.35×10⁴ 2.40×10⁴ 2.45×10⁴ 2.50×10⁴ 2.55×10⁴ 2.60×10⁴ 2.65×10⁴ 2.70×10⁴ 2.75×10⁴ 2.80×10⁴ 2.85×10⁴ 2.90×10⁴ 2.95×10⁴ 3.00×10⁴ 0 1×10⁴ 2×10⁴ 3×10⁴ 0 1.0×10³ 2.0×10³ 3.0×10³ 4.0×10³ 5.0×10³ 6.0×10³ 7.0×10³ 8.0×10³ 9.0×10³ 1.0×10⁴ 1.1×10⁴ 1.2×10⁴ 1.3×10⁴ 1.4×10⁴ 1.5×10⁴ 1.6×10⁴ 1.7×10⁴ 1.8×10⁴ 1.9×10⁴ 2.0×10⁴ 2.1×10⁴ 2.2×10⁴ 2.3×10⁴ 2.4×10⁴ 2.5×10⁴ 2.6×10⁴ 2.7×10⁴ 2.8×10⁴ 2.9×10⁴ 3.0×10⁴ 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8 4.60 4.65 4.70 4.75 4.80 4.85 4.90 4.95 5.00 5.05 5.10 5.15 5.20 5.25 5.30 5.35 5.40 5.45 5.50 5.55 5.60 5.65 5.70 5.75 5.80 5.85 5.90 5.95 6.00 6.05 6.10 6.15 6.20 6.25 6.30 6.35 6.40 6.45 6.50 6.55 6.60 6.65 6.70 6.75 6.80 6.85 6.90 6.95 7.00 7.05 7.10 7.15 7.20 7.25 7.30 7.35 7.40 7.45 7.50 7.55 7.60 4 5 6 7 8 4.6 4.7 4.8 4.9 5.0 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6.0 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 7.0 7.1 7.2 7.3 7.4 7.5 7.6 y

In [18]:
plot(x=burnin:n, y=sim.value[burnin:n, 3, 1], Geom.line)


Out[18]:
x -5.0×10³ 0 5.0×10³ 1.0×10⁴ 1.5×10⁴ 2.0×10⁴ 2.5×10⁴ 3.0×10⁴ 3.5×10⁴ 0 5.00×10² 1.00×10³ 1.50×10³ 2.00×10³ 2.50×10³ 3.00×10³ 3.50×10³ 4.00×10³ 4.50×10³ 5.00×10³ 5.50×10³ 6.00×10³ 6.50×10³ 7.00×10³ 7.50×10³ 8.00×10³ 8.50×10³ 9.00×10³ 9.50×10³ 1.00×10⁴ 1.05×10⁴ 1.10×10⁴ 1.15×10⁴ 1.20×10⁴ 1.25×10⁴ 1.30×10⁴ 1.35×10⁴ 1.40×10⁴ 1.45×10⁴ 1.50×10⁴ 1.55×10⁴ 1.60×10⁴ 1.65×10⁴ 1.70×10⁴ 1.75×10⁴ 1.80×10⁴ 1.85×10⁴ 1.90×10⁴ 1.95×10⁴ 2.00×10⁴ 2.05×10⁴ 2.10×10⁴ 2.15×10⁴ 2.20×10⁴ 2.25×10⁴ 2.30×10⁴ 2.35×10⁴ 2.40×10⁴ 2.45×10⁴ 2.50×10⁴ 2.55×10⁴ 2.60×10⁴ 2.65×10⁴ 2.70×10⁴ 2.75×10⁴ 2.80×10⁴ 2.85×10⁴ 2.90×10⁴ 2.95×10⁴ 3.00×10⁴ 0 1×10⁴ 2×10⁴ 3×10⁴ 0 1.0×10³ 2.0×10³ 3.0×10³ 4.0×10³ 5.0×10³ 6.0×10³ 7.0×10³ 8.0×10³ 9.0×10³ 1.0×10⁴ 1.1×10⁴ 1.2×10⁴ 1.3×10⁴ 1.4×10⁴ 1.5×10⁴ 1.6×10⁴ 1.7×10⁴ 1.8×10⁴ 1.9×10⁴ 2.0×10⁴ 2.1×10⁴ 2.2×10⁴ 2.3×10⁴ 2.4×10⁴ 2.5×10⁴ 2.6×10⁴ 2.7×10⁴ 2.8×10⁴ 2.9×10⁴ 3.0×10⁴ -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8 8.0 8.2 8.4 8.6 8.8 9.0 9.2 9.4 9.6 9.8 10.0 10.2 10.4 10.6 10.8 11.0 11.2 11.4 11.6 11.8 12.0 12.2 12.4 12.6 12.8 13.0 13.2 13.4 13.6 13.8 14.0 14.2 14.4 14.6 14.8 15.0 -5 0 5 10 15 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 15.0 y