In [1]:

    

include("ilmtools.jl")


    

    












    












    
Out[1]:





create_loglikelihood (generic function with 1 method)

In [2]:

    

pop_db1 = create_pop_db_univariates(100, Uniform(0,10), Uniform(0,10))


    

    
Out[2]:




ind_id
x
y
1
1
4.1959307141148505
8.345040165635586
2
2
0.757855949502162
2.8627572180673333
3
3
0.6867987378469431
4.994887655952017
4
4
7.296195861300429
5.160499597359864
5
5
1.6084775590068712
7.937516460440543
6
6
7.18159614845665
0.8084079030469238
7
7
8.450015049946607
1.9899835909128916
8
8
5.367445910080182
3.502961934477753
9
9
1.9279391723095651
2.9981330856735533
10
10
0.8633187128834874
3.89319942952856
11
11
0.07271448682305781
6.708507256143637
12
12
8.567441477733267
9.509432856860327
13
13
0.9020116978575965
3.399215417995618
14
14
5.25415517879662
2.4701519813766093
15
15
2.0709212942856103
6.147927749582873
16
16
9.451541410561251
2.4957134143648507
17
17
7.409850410675558
4.660549787000139
18
18
9.881765118108408
7.136547301380782
19
19
3.921458954700787
2.2803619566965283
20
20
1.690996321431042
8.476806412018693
21
21
1.3468285911243982
6.956183225178854
22
22
9.697256005757964
5.878702209204876
23
23
8.059190536416365
8.462256451002052
24
24
1.1020420520140206
0.674768351108157
25
25
8.88032790275243
5.206400278134051
26
26
3.1026552770039895
5.543952647485863
27
27
6.990753821503426
9.517107980418231
28
28
0.07787412870816413
4.214877739481945
29
29
2.173887588510137
7.213152348451124
30
30
4.355545894327227
3.969627194635257
&vellip
&vellip
&vellip
&vellip

In [3]:

    

evdb = infect_recover_loop(pop_db1, "discrete", "SIR", 2, 6, 5)


    

    
Out[3]:





edb(100x4 DataFrame
| Row | ind_id | s | i    | r    |
|-----|--------|---|------|------|
| 1   | 1      | 0 | 10.0 | 11.0 |
| 2   | 2      | 0 | 8.0  | 9.0  |
| 3   | 3      | 0 | 8.0  | 9.0  |
| 4   | 4      | 0 | 4.0  | 10.0 |
| 5   | 5      | 0 | 9.0  | 15.0 |
| 6   | 6      | 0 | 4.0  | 7.0  |
| 7   | 7      | 0 | 3.0  | 6.0  |
| 8   | 8      | 0 | 4.0  | 6.0  |
| 9   | 9      | 0 | 6.0  | 7.0  |
| 10  | 10     | 0 | 7.0  | 11.0 |
| 11  | 11     | 0 | 9.0  | 14.0 |
⋮
| 89  | 89     | 0 | 6.0  | 19.0 |
| 90  | 90     | 0 | 3.0  | 10.0 |
| 91  | 91     | 0 | 4.0  | 8.0  |
| 92  | 92     | 0 | 2.0  | 7.0  |
| 93  | 93     | 0 | 9.0  | 12.0 |
| 94  | 94     | 0 | 8.0  | 11.0 |
| 95  | 95     | 0 | 4.0  | 9.0  |
| 96  | 96     | 0 | 8.0  | 20.0 |
| 97  | 97     | 0 | 8.0  | 12.0 |
| 98  | 98     | 0 | 7.0  | 8.0  |
| 99  | 99     | 0 | 11.0 | 31.0 |
| 100 | 100    | 0 | 10.0 | 15.0 |,[1.0,2.0,2.0,2.0,2.0,2.0,3.0,3.0,3.0,3.0  …  18.0,18.0,19.0,19.0,19.0,20.0,20.0,26.0,31.0,41.0],"discrete")

In [4]:

    

ll = create_loglikelihood(pop_db1, evdb)


    

    
Out[4]:





dSIR_ll (generic function with 1 method)

In [5]:

    

ll((2, 6, 5))
@time ll((2, 6, 5))


    

    




elapsed time: 0.000945815 seconds (306736 bytes allocated)






    
Out[5]:





-334.73176378026545

In [6]:

    

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: 50.989012571 seconds (11789386316 bytes allocated, 24.36% 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.48586  0.324372  0.00229365   0.00666812  6879.46  
 "beta"   5.41494  0.489484  0.00346117   0.01091     6344.96  
 "gamma"  5.29738  0.485908  0.00343589   0.0120808   5688.17  

Quantiles:
4x6 Array{Any,2}:
 ""        "2.5%"    "25.0%"   "50.0%"   "75.0%"   "97.5%"
 "alpha"  0.949197  1.24934   1.45223   1.68763   2.19406 
 "beta"   4.53964   5.05828   5.39032   5.74214   6.42832 
 "gamma"  4.46216   4.96115   5.25849   5.59741   6.33642 

In [7]:

    

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


    

    
Out[7]:

In [8]:

    

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


    

    
Out[8]:

In [9]:

    

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


    

    
Out[9]: