From c43a52b196452f8a6b08b8397f33aae829d7648f Mon Sep 17 00:00:00 2001 From: "Documenter.jl" Date: Fri, 8 Mar 2024 13:05:15 +0000 Subject: [PATCH] build based on d10b1e5 --- previews/PR123/.documenter-siteinfo.json | 2 +- previews/PR123/checklist/index.html | 2 +- previews/PR123/contributing/index.html | 2 +- previews/PR123/examples/getting_started.jl | 46 ++++++- .../PR123/examples/getting_started/index.html | 122 +++++++++++------- previews/PR123/index.html | 2 +- previews/PR123/lib/internals/index.html | 23 +++- previews/PR123/lib/public/index.html | 2 +- previews/PR123/man/guide/index.html | 2 +- previews/PR123/objects.inv | Bin 785 -> 867 bytes previews/PR123/release-notes/index.html | 2 +- previews/PR123/search_index.js | 2 +- 12 files changed, 142 insertions(+), 65 deletions(-) diff --git a/previews/PR123/.documenter-siteinfo.json b/previews/PR123/.documenter-siteinfo.json index 520121a96..2ea0182c4 100644 --- a/previews/PR123/.documenter-siteinfo.json +++ b/previews/PR123/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.2","generation_timestamp":"2024-03-08T12:39:30","documenter_version":"1.3.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.2","generation_timestamp":"2024-03-08T13:05:11","documenter_version":"1.3.0"}} \ No newline at end of file diff --git a/previews/PR123/checklist/index.html b/previews/PR123/checklist/index.html index d0a3cf201..fb8dbb4cf 100644 --- a/previews/PR123/checklist/index.html +++ b/previews/PR123/checklist/index.html @@ -20,4 +20,4 @@ Either of those should automatically publish a new version to the Julia registry. - Once registered, the `TagBot.yml` workflow should create a tag, and rebuild the documentation for this tag. - - These steps can take quite a bit of time (1 hour or more), so don't be surprised if the new documentation takes a while to appear. + - These steps can take quite a bit of time (1 hour or more), so don't be surprised if the new documentation takes a while to appear. diff --git a/previews/PR123/contributing/index.html b/previews/PR123/contributing/index.html index 9ce136c40..376bda799 100644 --- a/previews/PR123/contributing/index.html +++ b/previews/PR123/contributing/index.html @@ -1,2 +1,2 @@ -Contributing · EpiAware.jl
GitHub

Contributing

This page details the some of the guidelines that should be followed when contributing to this package. It is adapted from Documenter.jl.

Branches

release-* branches are used for tagged minor versions of this package. This follows the same approach used in the main Julia repository, albeit on a much more modest scale.

Please open pull requests against the master branch rather than any of the release-* branches whenever possible.

Backports

Bug fixes are backported to the release-* branches using git cherry-pick -x by a EpiAware member and will become available in point releases of that particular minor version of the package.

Feel free to nominate commits that should be backported by opening an issue. Requests for new point releases to be tagged in METADATA.jl can also be made in the same way.

release-* branches

  • Each new minor version x.y.0 gets a branch called release-x.y (a protected branch).
  • New versions are usually tagged only from the release-x.y branches.
  • For patch releases, changes get backported to the release-x.y branch via a single PR with the standard name "Backports for x.y.z" and label "Type: Backport". The PR message links to all the PRs that are providing commits to the backport. The PR gets merged as a merge commit (i.e. not squashed).
  • The old release-* branches may be removed once they have outlived their usefulness.
  • Patch version milestones are used to keep track of which PRs get backported etc.

Style Guide

Follow the style of the surrounding text when making changes. When adding new features please try to stick to the following points whenever applicable. This project follows the SciML style guide.

Tests

Unit tests

As is conventional for Julia packages, unit tests are located at test/*.jl with the entrypoint test/runtests.jl.

End to end testing

Tests that build example package docs from source and inspect the results (end to end tests) are located in /test/examples. The main entry points are test/examples/make.jl for building and test/examples/test.jl for doing some basic checks on the generated outputs.

+Contributing · EpiAware.jl
GitHub

Contributing

This page details the some of the guidelines that should be followed when contributing to this package. It is adapted from Documenter.jl.

Branches

release-* branches are used for tagged minor versions of this package. This follows the same approach used in the main Julia repository, albeit on a much more modest scale.

Please open pull requests against the master branch rather than any of the release-* branches whenever possible.

Backports

Bug fixes are backported to the release-* branches using git cherry-pick -x by a EpiAware member and will become available in point releases of that particular minor version of the package.

Feel free to nominate commits that should be backported by opening an issue. Requests for new point releases to be tagged in METADATA.jl can also be made in the same way.

release-* branches

  • Each new minor version x.y.0 gets a branch called release-x.y (a protected branch).
  • New versions are usually tagged only from the release-x.y branches.
  • For patch releases, changes get backported to the release-x.y branch via a single PR with the standard name "Backports for x.y.z" and label "Type: Backport". The PR message links to all the PRs that are providing commits to the backport. The PR gets merged as a merge commit (i.e. not squashed).
  • The old release-* branches may be removed once they have outlived their usefulness.
  • Patch version milestones are used to keep track of which PRs get backported etc.

Style Guide

Follow the style of the surrounding text when making changes. When adding new features please try to stick to the following points whenever applicable. This project follows the SciML style guide.

Tests

Unit tests

As is conventional for Julia packages, unit tests are located at test/*.jl with the entrypoint test/runtests.jl.

End to end testing

Tests that build example package docs from source and inspect the results (end to end tests) are located in /test/examples. The main entry points are test/examples/make.jl for building and test/examples/test.jl for doing some basic checks on the generated outputs.

diff --git a/previews/PR123/examples/getting_started.jl b/previews/PR123/examples/getting_started.jl index 1ce170fc0..555145e1e 100644 --- a/previews/PR123/examples/getting_started.jl +++ b/previews/PR123/examples/getting_started.jl @@ -27,6 +27,7 @@ begin using Statistics using DataFramesMeta using LinearAlgebra + using Transducers end # ╔═╡ 3ebc8384-f73d-4597-83a7-07a3744fed61 @@ -243,14 +244,14 @@ However, we now treat the generated data as `truth_data` and make inference with We do the inference by MCMC/NUTS using the `Turing` NUTS sampler with default warm-up steps. " -# ╔═╡ c8ce0d46-a160-4c40-a055-69b3d10d1770 -truth_data = generated_obs - # ╔═╡ 4a4c6e91-8d8f-4bbf-bb7e-a36dc281e312 md" The observation model supports partially complete data. To test this we set some of the generated observations to be `missing`. " +# ╔═╡ 525aa98c-d0e5-4ffa-b808-d90fc986204c +truth_data = generated_obs + # ╔═╡ 259a7042-e74f-43c7-aeb4-97a3beeb7776 let truth_data = Union{Int, Missing}[truth_data...] @@ -270,12 +271,43 @@ inference_mdl = fix( (rw_init = 0.0,) ) +# ╔═╡ 9222b436-9445-4039-abbf-25c8cddb7f63 +md" +### Initialising inference + +It is possible for the default warm-up process for NUTS to get stuck in low probability or otherwise degenerate regions of parameter space. + +To make NUTS more robust we provide `manypathfinder`, which is built on pathfinder variational inference from [Pathfinder.jl](https://mlcolab.github.io/Pathfinder.jl/stable/). `manypathfinder` runs `nruns` pathfinder processes on the inference problem and returns the pathfinder run with maximum estimated ELBO. + +`manypathfinder` differs from `Pathfinder.multipathfinder`; `multipathfinder` is aimed at sampling from a potentially non-Gaussian target distribution which is first approximated as a uniformly weighted collection of normal approximations from pathfinder runs. `manypathfinder` is aimed at moving rapidly to a 'good' part of parameter space, and is robust to runs that fail. +" + +# ╔═╡ 197a4fbb-b71a-475a-bb78-28ff613e3094 +best_pf = manypathfinder(inference_mdl, 10; nruns = 20, executor = Transducers.ThreadedEx()); + +# ╔═╡ 073a1d40-456a-450e-969f-11b23eb7fd1f +md" +We can use draws from the best pathfinder run to initialise NUTS. +" + +# ╔═╡ 0379b058-4c35-440a-bc01-aafa0178bdbf +best_pf.draws_transformed + +# ╔═╡ a7798f71-9bb5-4506-9476-0cc11553b9e2 +init_params = collect.(eachrow(best_pf.draws_transformed.value[1:4, :, 1])) + +# ╔═╡ 4deb3a51-781d-48c4-91f6-6adf2b1affcf +md" +**NB: We are running this inference run for speed rather than accuracy as a demonstration. Use a higher target acceptance and more samples in a typical workflow.** +" + # ╔═╡ 3eb5ec5e-aae7-478e-84fb-80f2e9f85b4c chn = sample(inference_mdl, NUTS(; adtype = AutoReverseDiff(true)), MCMCThreads(), 250, 4; + init_params, drop_warmup = true) # ╔═╡ 30498cc7-16a5-441a-b8cd-c19b220c60c1 @@ -407,11 +439,17 @@ end # ╠═a04f3c1b-7e11-4800-9c2a-9fc0021de6e7 # ╟─f68b4e41-ac5c-42cd-a8c2-8761d66f7543 # ╟─b5bc8f05-b538-4abf-aa84-450bf2dff3d9 -# ╠═c8ce0d46-a160-4c40-a055-69b3d10d1770 # ╟─4a4c6e91-8d8f-4bbf-bb7e-a36dc281e312 +# ╠═525aa98c-d0e5-4ffa-b808-d90fc986204c # ╠═259a7042-e74f-43c7-aeb4-97a3beeb7776 # ╟─32638954-2c99-4d4e-8e03-52154030c657 # ╠═b4033728-b321-4100-8194-1fd9fe2d268d +# ╟─9222b436-9445-4039-abbf-25c8cddb7f63 +# ╠═197a4fbb-b71a-475a-bb78-28ff613e3094 +# ╠═073a1d40-456a-450e-969f-11b23eb7fd1f +# ╠═0379b058-4c35-440a-bc01-aafa0178bdbf +# ╠═a7798f71-9bb5-4506-9476-0cc11553b9e2 +# ╟─4deb3a51-781d-48c4-91f6-6adf2b1affcf # ╠═3eb5ec5e-aae7-478e-84fb-80f2e9f85b4c # ╟─30498cc7-16a5-441a-b8cd-c19b220c60c1 # ╠═e9df22b8-8e4d-4ab7-91ea-c01f2239b3e5 diff --git a/previews/PR123/examples/getting_started/index.html b/previews/PR123/examples/getting_started/index.html index 43af2f465..129034ae1 100644 --- a/previews/PR123/examples/getting_started/index.html +++ b/previews/PR123/examples/getting_started/index.html @@ -25,7 +25,7 @@ @@ -49,6 +49,7 @@ using Statistics using DataFramesMeta using LinearAlgebra + using Transducers end @@ -130,67 +131,67 @@ 
Model{typeof(make_epi_aware), (:y_t, :time_steps), (:epi_model, :latent_model, :observation_model, :pos_shift), (:y_t,), Tuple{Missing, Int64}, Tuple{DirectInfections{Normal{Float64}}, RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, Float64}, DynamicPPL.FixedContext{@NamedTuple{rw_init::Float64, init_incidence::Float64}, DefaultContext}}(EpiAware.make_epi_aware, (y_t = missing, time_steps = 100), (epi_model = DirectInfections{Normal{Float64}}(EpiData{Float64, typeof(exp)}([0.056722565637507084, 0.0944813112842689, 0.11668723746973605, 0.12773814268679387, 0.130948497124998, 0.1287976856448402, 0.12312384970427602, 0.11527489380975135, 0.10622581663782862], 9, exp), Distributions.Normal{Float64}(μ=4.605170185988092, σ=1.0)), latent_model = RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(Distributions.Normal{Float64}(μ=0.0, σ=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.02); lower=0.0, upper=Inf)), observation_model = DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(sparse(Integer[1, 2, 3, 4, 2, 3, 4, 5, 3, 4  …  97, 98, 99, 100, 98, 99, 100, 99, 100, 100], Integer[1, 1, 1, 1, 2, 2, 2, 2, 3, 3  …  97, 97, 97, 97, 98, 98, 98, 99, 99, 100], [0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25  …  0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25], 100, 100), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.12533141373155002); lower=0.0, upper=Inf)), pos_shift = 1.0e-6), FixedContext((rw_init = 0.0, init_incidence = 4.605170185988092), DynamicPPL.DefaultContext()))
random_epidemic = rand(cond_generative_model)
-
(ϵ_t = [0.7541364299786495, 0.6228295155052502, -0.08437805815079774, -1.6027945048950885, 0.5288638281567213, -1.6896912280141976, -1.6475552716129747, -1.4734033626550678, 0.2292953790785949, -1.4453965284513017  …  -0.5319988979469747, -0.1702093780381867, 0.7423616167719554, 0.5084938565160623, -0.7552829192371855, -0.12043077985555731, -1.7538972893049496, 0.5572895108833601, -0.7266998274413684, 1.0439092851436422], σ²_RW = 0.030509604813344984, neg_bin_cluster_factor = 0.003822359489741684, var"y_t[1]" = 27, var"y_t[2]" = 63, var"y_t[3]" = 75, var"y_t[4]" = 121, var"y_t[5]" = 126, var"y_t[6]" = 99, var"y_t[7]" = 84, var"y_t[8]" = 59, var"y_t[9]" = 62, var"y_t[10]" = 42, var"y_t[11]" = 60, var"y_t[12]" = 36, var"y_t[13]" = 52, var"y_t[14]" = 47, var"y_t[15]" = 47, var"y_t[16]" = 47, var"y_t[17]" = 35, var"y_t[18]" = 58, var"y_t[19]" = 55, var"y_t[20]" = 62, var"y_t[21]" = 66, var"y_t[22]" = 88, var"y_t[23]" = 73, var"y_t[24]" = 78, var"y_t[25]" = 73, var"y_t[26]" = 68, var"y_t[27]" = 81, var"y_t[28]" = 64, var"y_t[29]" = 78, var"y_t[30]" = 68, var"y_t[31]" = 71, var"y_t[32]" = 80, var"y_t[33]" = 95, var"y_t[34]" = 74, var"y_t[35]" = 79, var"y_t[36]" = 86, var"y_t[37]" = 101, var"y_t[38]" = 94, var"y_t[39]" = 89, var"y_t[40]" = 93, var"y_t[41]" = 77, var"y_t[42]" = 90, var"y_t[43]" = 98, var"y_t[44]" = 108, var"y_t[45]" = 129, var"y_t[46]" = 130, var"y_t[47]" = 115, var"y_t[48]" = 126, var"y_t[49]" = 99, var"y_t[50]" = 121, var"y_t[51]" = 96, var"y_t[52]" = 111, var"y_t[53]" = 82, var"y_t[54]" = 80, var"y_t[55]" = 67, var"y_t[56]" = 80, var"y_t[57]" = 91, var"y_t[58]" = 87, var"y_t[59]" = 103, var"y_t[60]" = 118, var"y_t[61]" = 129, var"y_t[62]" = 138, var"y_t[63]" = 133, var"y_t[64]" = 150, var"y_t[65]" = 163, var"y_t[66]" = 190, var"y_t[67]" = 219, var"y_t[68]" = 216, var"y_t[69]" = 225, var"y_t[70]" = 335, var"y_t[71]" = 362, var"y_t[72]" = 447, var"y_t[73]" = 469, var"y_t[74]" = 460, var"y_t[75]" = 465, var"y_t[76]" = 458, var"y_t[77]" = 457, var"y_t[78]" = 470, var"y_t[79]" = 473, var"y_t[80]" = 516, var"y_t[81]" = 631, var"y_t[82]" = 606, var"y_t[83]" = 593, var"y_t[84]" = 563, var"y_t[85]" = 505, var"y_t[86]" = 451, var"y_t[87]" = 453, var"y_t[88]" = 427, var"y_t[89]" = 352, var"y_t[90]" = 310, var"y_t[91]" = 243, var"y_t[92]" = 277, var"y_t[93]" = 247, var"y_t[94]" = 238, var"y_t[95]" = 248, var"y_t[96]" = 251, var"y_t[97]" = 229, var"y_t[98]" = 234, var"y_t[99]" = 187, var"y_t[100]" = 193)
+
(ϵ_t = [-0.20112065205042204, 0.3243818279907073, -0.9802665412256135, 1.5899093334448624, 1.1765040928990662, -1.097368025541141, 0.21404530796180515, 0.036040772035425335, 0.19674066706454033, -1.1183113744954276  …  2.925721551631621, -0.2811334030236948, -1.0839466539012816, -0.01945211239319567, -0.525340621257463, -1.5377346782206847, -1.1521273618730166, 0.9841152135889327, 0.03839535544638948, -0.9408477622299963], σ²_RW = 0.007678264825513905, neg_bin_cluster_factor = 0.07657768923618172, var"y_t[1]" = 28, var"y_t[2]" = 59, var"y_t[3]" = 73, var"y_t[4]" = 99, var"y_t[5]" = 94, var"y_t[6]" = 118, var"y_t[7]" = 93, var"y_t[8]" = 90, var"y_t[9]" = 95, var"y_t[10]" = 96, var"y_t[11]" = 94, var"y_t[12]" = 110, var"y_t[13]" = 87, var"y_t[14]" = 86, var"y_t[15]" = 88, var"y_t[16]" = 80, var"y_t[17]" = 88, var"y_t[18]" = 74, var"y_t[19]" = 87, var"y_t[20]" = 77, var"y_t[21]" = 100, var"y_t[22]" = 103, var"y_t[23]" = 87, var"y_t[24]" = 167, var"y_t[25]" = 117, var"y_t[26]" = 126, var"y_t[27]" = 152, var"y_t[28]" = 153, var"y_t[29]" = 179, var"y_t[30]" = 184, var"y_t[31]" = 212, var"y_t[32]" = 160, var"y_t[33]" = 155, var"y_t[34]" = 151, var"y_t[35]" = 170, var"y_t[36]" = 172, var"y_t[37]" = 192, var"y_t[38]" = 205, var"y_t[39]" = 185, var"y_t[40]" = 168, var"y_t[41]" = 169, var"y_t[42]" = 149, var"y_t[43]" = 161, var"y_t[44]" = 168, var"y_t[45]" = 184, var"y_t[46]" = 187, var"y_t[47]" = 215, var"y_t[48]" = 255, var"y_t[49]" = 254, var"y_t[50]" = 209, var"y_t[51]" = 221, var"y_t[52]" = 190, var"y_t[53]" = 256, var"y_t[54]" = 257, var"y_t[55]" = 240, var"y_t[56]" = 245, var"y_t[57]" = 222, var"y_t[58]" = 274, var"y_t[59]" = 221, var"y_t[60]" = 229, var"y_t[61]" = 264, var"y_t[62]" = 324, var"y_t[63]" = 323, var"y_t[64]" = 280, var"y_t[65]" = 281, var"y_t[66]" = 235, var"y_t[67]" = 222, var"y_t[68]" = 262, var"y_t[69]" = 248, var"y_t[70]" = 271, var"y_t[71]" = 232, var"y_t[72]" = 243, var"y_t[73]" = 262, var"y_t[74]" = 253, var"y_t[75]" = 297, var"y_t[76]" = 250, var"y_t[77]" = 304, var"y_t[78]" = 283, var"y_t[79]" = 289, var"y_t[80]" = 303, var"y_t[81]" = 345, var"y_t[82]" = 301, var"y_t[83]" = 230, var"y_t[84]" = 226, var"y_t[85]" = 264, var"y_t[86]" = 243, var"y_t[87]" = 265, var"y_t[88]" = 259, var"y_t[89]" = 301, var"y_t[90]" = 306, var"y_t[91]" = 253, var"y_t[92]" = 257, var"y_t[93]" = 337, var"y_t[94]" = 292, var"y_t[95]" = 246, var"y_t[96]" = 275, var"y_t[97]" = 256, var"y_t[98]" = 259, var"y_t[99]" = 212, var"y_t[100]" = 226)
true_infections = generated_quantities(cond_generative_model, random_epidemic).I_t
100-element Vector{Float64}:
- 114.07945592233989
- 127.19035360323254
- 125.32952857511957
-  94.72579974976007
- 103.89315740245587
-  77.3410367647116
-  58.00018680718289
+  98.2531043136177
+ 101.08593876245777
+  92.76547032111293
+ 106.63282760687878
+ 118.21244895560318
+ 107.37482705017699
+ 109.40774047821078
    ⋮
- 248.0613826703851
- 242.89775255149056
- 178.80348692965416
- 197.0838192629199
- 173.59004678270438
- 208.31181943121163
+ 271.51631386943313 + 237.2886030084489 + 214.50244950617125 + 233.82078619031455 + 234.6087824059075 + 216.04288275076027 
generated_obs = generated_quantities(cond_generative_model, random_epidemic).generated_y_t
100-element Vector{Int64}:
-  27
-  63
-  75
- 121
- 126
+  28
+  59
+  73
   99
-  84
+  94
+ 118
+  93
    ⋮
- 248
- 251
- 229
- 234
- 187
- 193
+ 246 + 275 + 256 + 259 + 212 + 226 - +

Inference

Fixing $Z_0 = 0$ for the random walk was based on inference principles; in this model $Z_0$ and $\log I_0$ are non-identifiable.

However, we now treat the generated data as truth_data and make inference without fixing any other parameters.

We do the inference by MCMC/NUTS using the Turing NUTS sampler with default warm-up steps.

+ +

The observation model supports partially complete data. To test this we set some of the generated observations to be missing.

+ 
truth_data = generated_obs
100-element Vector{Int64}:
-  27
-  63
-  75
- 121
- 126
+  28
+  59
+  73
   99
-  84
+  94
+ 118
+  93
    ⋮
- 248
- 251
- 229
- 234
- 187
- 193
- - -

The observation model supports partially complete data. To test this we set some of the generated observations to be missing.

+ 246 + 275 + 256 + 259 + 212 + 226 
let
     truth_data = Union{Int, Missing}[truth_data...]
@@ -220,15 +221,38 @@ 
Inference
-
Model{typeof(make_epi_aware), (:y_t, :time_steps), (:epi_model, :latent_model, :observation_model, :pos_shift), (), Tuple{Vector{Int64}, Int64}, Tuple{DirectInfections{Normal{Float64}}, RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, Float64}, DynamicPPL.FixedContext{@NamedTuple{rw_init::Float64}, DefaultContext}}(EpiAware.make_epi_aware, (y_t = [27, 63, 75, 121, 126, 99, 84, 59, 62, 42  …  243, 277, 247, 238, 248, 251, 229, 234, 187, 193], time_steps = 100), (epi_model = DirectInfections{Normal{Float64}}(EpiData{Float64, typeof(exp)}([0.056722565637507084, 0.0944813112842689, 0.11668723746973605, 0.12773814268679387, 0.130948497124998, 0.1287976856448402, 0.12312384970427602, 0.11527489380975135, 0.10622581663782862], 9, exp), Distributions.Normal{Float64}(μ=4.605170185988092, σ=1.0)), latent_model = RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(Distributions.Normal{Float64}(μ=0.0, σ=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.02); lower=0.0, upper=Inf)), observation_model = DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(sparse(Integer[1, 2, 3, 4, 2, 3, 4, 5, 3, 4  …  97, 98, 99, 100, 98, 99, 100, 99, 100, 100], Integer[1, 1, 1, 1, 2, 2, 2, 2, 3, 3  …  97, 97, 97, 97, 98, 98, 98, 99, 99, 100], [0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25  …  0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25], 100, 100), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.12533141373155002); lower=0.0, upper=Inf)), pos_shift = 1.0e-6), FixedContext((rw_init = 0.0,), DynamicPPL.DefaultContext()))
+
Model{typeof(make_epi_aware), (:y_t, :time_steps), (:epi_model, :latent_model, :observation_model, :pos_shift), (), Tuple{Vector{Int64}, Int64}, Tuple{DirectInfections{Normal{Float64}}, RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, Float64}, DynamicPPL.FixedContext{@NamedTuple{rw_init::Float64}, DefaultContext}}(EpiAware.make_epi_aware, (y_t = [28, 59, 73, 99, 94, 118, 93, 90, 95, 96  …  253, 257, 337, 292, 246, 275, 256, 259, 212, 226], time_steps = 100), (epi_model = DirectInfections{Normal{Float64}}(EpiData{Float64, typeof(exp)}([0.056722565637507084, 0.0944813112842689, 0.11668723746973605, 0.12773814268679387, 0.130948497124998, 0.1287976856448402, 0.12312384970427602, 0.11527489380975135, 0.10622581663782862], 9, exp), Distributions.Normal{Float64}(μ=4.605170185988092, σ=1.0)), latent_model = RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(Distributions.Normal{Float64}(μ=0.0, σ=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.02); lower=0.0, upper=Inf)), observation_model = DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(sparse(Integer[1, 2, 3, 4, 2, 3, 4, 5, 3, 4  …  97, 98, 99, 100, 98, 99, 100, 99, 100, 100], Integer[1, 1, 1, 1, 2, 2, 2, 2, 3, 3  …  97, 97, 97, 97, 98, 98, 98, 99, 99, 100], [0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25  …  0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25], 100, 100), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.12533141373155002); lower=0.0, upper=Inf)), pos_shift = 1.0e-6), FixedContext((rw_init = 0.0,), DynamicPPL.DefaultContext()))
+ + +

Initialising inference

It is possible for the default warm-up process for NUTS to get stuck in low probability or otherwise degenerate regions of parameter space.

To make NUTS more robust we provide manypathfinder, which is built on pathfinder variational inference from Pathfinder.jl. manypathfinder runs nruns pathfinder processes on the inference problem and returns the pathfinder run with maximum estimated ELBO.

manypathfinder differs from Pathfinder.multipathfinder; multipathfinder is aimed at sampling from a potentially non-Gaussian target distribution which is first approximated as a uniformly weighted collection of normal approximations from pathfinder runs. manypathfinder is aimed at moving rapidly to a 'good' part of parameter space, and is robust to runs that fail.

+ +
best_pf = manypathfinder(inference_mdl, 10; nruns = 20, executor = Transducers.ThreadedEx());
+ + + +

We can use draws from the best pathfinder run to initialise NUTS.

+ +
best_pf.draws_transformed
+
iterationchainϵ_t[1]ϵ_t[2]ϵ_t[3]ϵ_t[4]ϵ_t[5]ϵ_t[6]...
1111.745171.3581-0.0387770.228527-0.4804070.259776
2211.739621.966440.2950680.6561230.2273570.221215
3311.404741.62267-0.1437840.913016-0.212709-0.0593607
4411.989611.71673-0.1428980.558126-0.3282210.278743
5511.839921.609580.3763271.052640.1482740.358493
6611.568591.3840.2471130.647590.08802540.372194
7711.802511.857780.1909470.6891620.1628090.01183
8812.126171.599590.1070831.17275-0.2168740.267504
9911.517171.88996-0.2935890.3183850.146870.335505
101011.622831.57171-0.2107330.3872-0.5130320.182996
+ +
init_params = collect.(eachrow(best_pf.draws_transformed.value[1:4, :, 1]))
+
4-element Vector{Vector{Float64}}:
+ [1.7451655099039896, 1.3581032532363655, -0.03877701856907928, 0.2285266700718785, -0.4804068463802501, 0.2597760274093325, -0.72329152773823, -0.004892028183825725, 0.30371317937314624, -0.32360713324738793  …  -0.5689987466999031, -1.069569400690149, -0.006438779107851761, 0.12884728013030794, -0.4672454241686336, 0.46663119811384574, -0.5012766699018397, 0.0008443002919903706, 4.503333034127224, 0.10534955177982713]
+ [1.7396239623829572, 1.9664442888461848, 0.29506847368478406, 0.6561230931785763, 0.2273567486516539, 0.22121457862999117, -0.4457969493025312, 0.4933162930390505, 0.5387079340846253, 0.27978737948037563  …  -0.26300000639899557, -1.0805768157412714, -0.2542903280341239, 0.3450712794593742, -0.41739949283166994, 0.5189459878777011, -0.2407079378968638, 0.0013653332228509899, 4.40762882420073, 0.12133329079454364]
+ [1.404744492096187, 1.6226685394934721, -0.14378388615214202, 0.9130160620170369, -0.21270896000629097, -0.059360707272928415, -0.38729412142681957, 0.3423067106403553, 0.4423491622359259, 0.03856933348504296  …  -0.6845978102747797, -1.0692267124647084, -0.1627366413969122, 0.4973267181861608, -0.2835918876396665, 0.7286810493507772, -0.14793772810228983, 0.0012724611829520575, 4.364726695263788, 0.11203225705238223]
+ [1.9896074833972643, 1.7167261947293577, -0.14289800447970843, 0.5581263080306341, -0.3282212207720035, 0.27874272220508667, -0.9184680486663799, 0.5964059365665253, 0.5400965831532827, 0.5955039812253591  …  -0.6508965991738329, -1.2498544646684877, -0.35170979199800345, 0.3039492244589093, -0.5363086467209464, 0.37978414791959775, 0.2275674079496028, 0.0015740315635789245, 4.358921604618429, 0.09758372594224174]
+ + +

NB: We are running this inference run for speed rather than accuracy as a demonstration. Use a higher target acceptance and more samples in a typical workflow.

chn = sample(inference_mdl,
     NUTS(; adtype = AutoReverseDiff(true)),
     MCMCThreads(),
     250,
     4;
+    init_params,
     drop_warmup = true)
-
iterationchainϵ_t[1]ϵ_t[2]ϵ_t[3]ϵ_t[4]ϵ_t[5]ϵ_t[6]...
11261-1.236280.23389-0.485811-0.494185-0.2627-2.06568
212711.25850.289792-0.0274798-0.715777-0.759128-1.7599
31281-2.127861.34623-1.200970.124472-2.265070.198596
41291-2.086831.00838-1.834050.228425-2.24750.219736
513010.884139-0.0395667-0.8535551.94183-1.763630.482344
61311-0.5916391.949591.0209-1.98888-1.21572-2.21748
71321-0.2620560.8361030.315522-0.723914-1.51473-0.888812
81331-0.6898560.812266-1.08283-0.331231-1.17112-0.520497
913411.033671.28017-0.821226-1.472-0.90994-0.244272
101351-1.893661.18913-0.578158-1.30672-1.0126-2.69586
...
+
iterationchainϵ_t[1]ϵ_t[2]ϵ_t[3]ϵ_t[4]ϵ_t[5]ϵ_t[6]...
112610.871897-0.130904-0.919508-0.0402790.4041450.142873
21271-1.10953-1.335251.81422-1.51668-0.190370.143744
312810.6702260.648365-0.791261-0.151059-0.2185021.85856
412910.417116-0.272337-0.671253-0.6251520.595989-1.02231
513012.452411.239250.745333-1.320331.2241-1.87632
613111.973160.2079611.01515-2.481431.1235-1.28877
71321-0.590797-0.778022-0.1580630.748415-0.172206-0.517546
81331-1.402431.63989-0.339024-0.8965040.2595890.507557
91341-0.915431.201580.554195-0.549210.4315130.349227
1013510.510275-1.12904-0.3463790.104358-0.576478-0.813278
...

Predictive plotting

We can spaghetti plot generated case data from the version of the model which hasn't conditioned on case data using posterior parameters inferred from the version conditioned on observed data. This is known as posterior predictive checking, and is a useful diagnostic tool for Bayesian inference (see here).

Because we are using synthetic data we can also plot the model predictions for the unobserved infections and check that (at least in this example) we were able to capture some unobserved/latent variables in the process accurate.

@@ -265,7 +289,7 @@

Inference - +

As well as checking the posterior predictions for latent infections, we can also check how well inference recovered unknown parameters, such as the random walk variance or the cluster factor of the negative binomial observations.

@@ -286,7 +310,7 @@

Inference - +

Reproductive number back-calculation

As mentioned at the top, we don't directly use the concept of reproductive numbers in this note. However, we can back-calculate the implied $\mathcal{R}(t)$ values, conditional on the specified generation interval being correct.

Here we spaghetti plot posterior sampled time-varying reproductive numbers against the actual.

@@ -313,6 +337,6 @@

- + - + diff --git a/previews/PR123/index.html b/previews/PR123/index.html index aff5904bc..23ebeceb1 100644 --- a/previews/PR123/index.html +++ b/previews/PR123/index.html @@ -1,2 +1,2 @@ -EpiAware.jl: Real-time epidemic monitoring · EpiAware.jl
GitHub

EpiAware.jl

Infectious disease situational awareness modelling toolkit for Julia.

A package for building and fitting situational awareness models for infectious diseases. The package is designed to be flexible and extensible, and to provide a consistent interface for fitting and simulating models.

Package Features

  • Flexible: The package is designed to be flexible and extensible, and to provide a consistent interface for fitting and simulating models.
  • Modular: The package is designed to be modular, with a clear separation between the model and the data.
  • Extensible: The package is designed to be extensible, with a clear separation between the model and the data.
  • Consistent: The package is designed to provide a consistent interface for fitting and simulating models.
  • Efficient: The package is designed to be efficient, with a clear separation between the model and the data.

See the Index for the complete list of documented functions and types.

Manual Outline

    Library Outline

    Index

      +EpiAware.jl: Real-time epidemic monitoring · EpiAware.jl
      GitHub

      EpiAware.jl

      Infectious disease situational awareness modelling toolkit for Julia.

      A package for building and fitting situational awareness models for infectious diseases. The package is designed to be flexible and extensible, and to provide a consistent interface for fitting and simulating models.

      Package Features

      • Flexible: The package is designed to be flexible and extensible, and to provide a consistent interface for fitting and simulating models.
      • Modular: The package is designed to be modular, with a clear separation between the model and the data.
      • Extensible: The package is designed to be extensible, with a clear separation between the model and the data.
      • Consistent: The package is designed to provide a consistent interface for fitting and simulating models.
      • Efficient: The package is designed to be efficient, with a clear separation between the model and the data.

      See the Index for the complete list of documented functions and types.

      Manual Outline

        Library Outline

        Index

          diff --git a/previews/PR123/lib/internals/index.html b/previews/PR123/lib/internals/index.html index ae56876a0..b7a1c4e97 100644 --- a/previews/PR123/lib/internals/index.html +++ b/previews/PR123/lib/internals/index.html @@ -1,8 +1,23 @@ -Internals · EpiAware.jl
          GitHub

          internal Documentation

          Documentation for EpiAware.jl's internal interface.

          Contents

            Index

            Internal API

            EpiAware.NegativeBinomialMeanClust — Method

            NegativeBinomialMeanClust NegativeBinomialMeanClust(μ, α)

            Compute the mean-cluster factor negative binomial distribution.

            Arguments

            • μ: The mean of the distribution.
            • α: The clustering factor parameter.

            Returns

            A NegativeBinomial distribution object.

            Signatures

            NegativeBinomialMeanClust(
+Internals · EpiAware.jl
            GitHub
            internal Documentation
            Documentation for EpiAware.jl's internal interface.
            Contents
              Index
              Internal API
              EpiAware.NegativeBinomialMeanClust — Method
              NegativeBinomialMeanClust     NegativeBinomialMeanClust(μ, α)
              Compute the mean-cluster factor negative binomial distribution.
              Arguments
              • μ: The mean of the distribution.
              • α: The clustering factor parameter.
              Returns
              A NegativeBinomial distribution object.
              Signatures
              NegativeBinomialMeanClust(
     μ,
     α
 ) -> Distributions.NegativeBinomial
-
              Methods
              NegativeBinomialMeanClust(μ, α)
              defined at /home/runner/work/Rt-without-renewal/Rt-without-renewal/src/utilities.jl:210.
              source
              EpiAware.generate_observation_kernel — Method
              generateobservationkernel     generateobservationkernel(delayint, timehorizon)
              Generate an observation kernel matrix based on the given delay interval and time horizon.
              Arguments
              • delay_int::Vector{Float64}: The delay PMF vector.
              • time_horizon::Int: The number of time steps of the observation period.
              Returns
              • K::SparseMatrixCSC{Float64, Int}: The observation kernel matrix.
              Signatures
              generate_observation_kernel(delay_int, time_horizon) -> Any
-
              Methods
              generate_observation_kernel(delay_int, time_horizon)
              defined at /home/runner/work/Rt-without-renewal/Rt-without-renewal/src/utilities.jl:229.
              source
              
+

              Methods

              NegativeBinomialMeanClust(μ, α)

              defined at /home/runner/work/Rt-without-renewal/Rt-without-renewal/src/utilities.jl:210.

              source
              EpiAware._continue_manypathfinder! — Method

              continuemanypathfinder! Continue running the pathfinder algorithm until a pathfinder succeeds or the maximum number of tries is reached.

              Arguments

              • pfs: An array of pathfinder objects.
              • mdl::DynamicPPL.Model: The model to perform inference on.
              • max_tries: The maximum number of tries to run the pathfinder algorithm. Default is Inf.
              • nruns: The number of times to run the pathfinder function.
              • kwargs...: Additional keyword arguments passed to pathfinder.

              Returns

              • pfs: The updated array of pathfinder objects.

              Signatures

              _continue_manypathfinder!(
+    pfs,
+    mdl::DynamicPPL.Model;
+    max_tries,
+    nruns,
+    kwargs...
+)
+

              Methods

              _continue_manypathfinder!(
+    pfs,
+    mdl;
+    max_tries,
+    nruns,
+    kwargs...
+)

              defined at /home/runner/work/Rt-without-renewal/Rt-without-renewal/src/inference-methods.jl:42.

              source
              EpiAware._run_manypathfinder — Method

              runmanypathfinder

              Run pathfinder multiple times and store the results in an array. Fails safely.

              Arguments

              • mdl::DynamicPPL.Model: The Turing model to be used for inference.
              • nruns: The number of times to run the pathfinder function.
              • kwargs...: Additional keyword arguments passed to pathfinder.

              Returns

              An array of PathfinderResult objects or Symbol values indicating success or failure.

              Signatures

              _run_manypathfinder(mdl::DynamicPPL.Model; nruns, kwargs...)
+

              Methods

              _run_manypathfinder(mdl; nruns, kwargs...)

              defined at /home/runner/work/Rt-without-renewal/Rt-without-renewal/src/inference-methods.jl:13.

              source
              EpiAware.generate_observation_kernel — Method

              generateobservationkernel generateobservationkernel(delayint, timehorizon)

              Generate an observation kernel matrix based on the given delay interval and time horizon.

              Arguments

              • delay_int::Vector{Float64}: The delay PMF vector.
              • time_horizon::Int: The number of time steps of the observation period.

              Returns

              • K::SparseMatrixCSC{Float64, Int}: The observation kernel matrix.

              Signatures

              generate_observation_kernel(delay_int, time_horizon) -> Any
+

              Methods

              generate_observation_kernel(delay_int, time_horizon)

              defined at /home/runner/work/Rt-without-renewal/Rt-without-renewal/src/utilities.jl:229.

              source
              diff --git a/previews/PR123/lib/public/index.html b/previews/PR123/lib/public/index.html index 989536627..0a467faca 100644 --- a/previews/PR123/lib/public/index.html +++ b/previews/PR123/lib/public/index.html @@ -1,2 +1,2 @@ -Public API · EpiAware.jl
              GitHub
              +Public API · EpiAware.jl
              GitHub
              diff --git a/previews/PR123/man/guide/index.html b/previews/PR123/man/guide/index.html index d534b6187..188bf8013 100644 --- a/previews/PR123/man/guide/index.html +++ b/previews/PR123/man/guide/index.html @@ -1,2 +1,2 @@ -Guide · EpiAware.jl
              GitHub
              +Guide · EpiAware.jl
              GitHub
              diff --git a/previews/PR123/objects.inv b/previews/PR123/objects.inv index ab7c86d765b104813d035873b1e28cc1e88c6532..df6dec51b9c9f426c107d7d05ddb00d69ea64899 100644 GIT binary patch delta 747 zcmV_ep@p)451 zW!4wfPIz?Ld=gV}Fr_Ux8P3Fht@4*NRP08Fc>>CEE3_)pB@aLjm;C&NHb2x!BG0N` zDsuB(0hrIoQh&XI0TqIa9PL=0RBOd9oKV<6%Z5EcW-~Xz=~YSoj0R-(G_t_a5=U-T z#H?me3>awaFq(>ZYBhyosPLq+)v4mwCPa$bhHl*#8?NFrirzVdF1DPlr5vp7KtNH} z$ARI-7k5p!;wr~B{>-fu3I1#fBasH9E>eE8cLUc@ltO0q%i?T!y(kVB%xvy30l@ zHri(+i+}!|dqm!%aYTP3U8gW|C*IF!-h;sMnBAyX>-pBK?HN-L`>Xpq_LSqs|IXCz z8=VE2H9WL05aL_+kK$XHfZgbb4sNi+G?S}25>2Ly$@IHRHo3a#9&q2}%EvczmvXy# zEu3U(AvZh4@xfiMA1+op0(qW_Q9e}EXY5PCmSp<@UUxpoK>GthxS=yO94D%4<1QY9 z0>dr%I>^Ahm0t3znXd1gN(<^sJuiav7$1>!m<%N!kPf)~)E;UIW?YTrIK`$aL3W z*nc1W87m)6u#;`_@#*Q2lC@3o6$(nPQkQ`I0&ug| z>q7xtFU7+CM1PMq!4P@0V{K8ZmAi04VF#@n_J9QyVL{NVn*5OmP||5+g@YO=Zfq*a zNvs|m=*Hn_DdT}0VD)J5q8Y1G=C?!06zv#>zE7W(&CW1=a0o+eoo%HY^mb1`W$XRk zz#U&a9J*U&OK#)$l8m;vr=>|ln`ZQit&&Q%tUD(6r+=ufte~Je0b-(()!O;f?f*d= zgL2;AJ4c~4u)T!FX?lk3cdS8WC}*gs^t)D7o;CpJ0DJ{%D|LT~4zY_W+Mo3gGWN={ z1#j-JLys=TTYOIeT}VvM;Hoxsk8AJJG9mrK8!ht~y^MJeGEpWe&4D3)>s$}npkaeu zHZk43bAR8EH|-{3cO%0|kz%G_&uHFP&Gz`s(OGTt{jm0D%wo2@ozLV$$s7MOwwsSI zjtlZCp(jlEU-uXIUzj63yA<|ru){Q?%Lx;WE~cZ4&mr09@_Kl{GqmUx;nKwg`j^_n z<%-GC1cL15JQXMT0%iOX%K3(ws@glz*5z86TM|@uYRZECOUwzjgHrwhb%Obe)$vYT diff --git a/previews/PR123/release-notes/index.html b/previews/PR123/release-notes/index.html index fdcc68321..0cb4d6d1c 100644 --- a/previews/PR123/release-notes/index.html +++ b/previews/PR123/release-notes/index.html @@ -1,2 +1,2 @@ -Release notes · EpiAware.jl
              GitHub
              +Release notes · EpiAware.jl
              GitHub
              diff --git a/previews/PR123/search_index.js b/previews/PR123/search_index.js index 564dc502b..fccecafee 100644 --- a/previews/PR123/search_index.js +++ b/previews/PR123/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"checklist/#Checklists","page":"Release checklist","title":"Checklists","text":"","category":"section"},{"location":"checklist/","page":"Release checklist","title":"Release checklist","text":"The purpose of this page is to collate a series of checklists for commonly performed changes to the source code of EpiAware. It has been adapted from Documenter.jl.","category":"page"},{"location":"checklist/","page":"Release checklist","title":"Release checklist","text":"In each case, copy the checklist into the description of the pull request.","category":"page"},{"location":"checklist/#Making-a-release","page":"Release checklist","title":"Making a release","text":"","category":"section"},{"location":"checklist/","page":"Release checklist","title":"Release checklist","text":"In preparation for a release, use the following checklist. These steps should be performed on a branch with an open pull request, either for a topic branch, or for a new branch release-1.y.z (\"Release version 1.y.z\") if multiple changes have accumulated on the master branch since the last release.","category":"page"},{"location":"checklist/","page":"Release checklist","title":"Release checklist","text":"## Pre-release\n\n - [ ] Change the version number in `Project.toml`\n * If the release is breaking, increment MAJOR\n * If the release adds a new user-visible feature, increment MINOR\n * Otherwise (bug-fixes, documentation improvements), increment PATCH\n - [ ] Update `CHANGELOG.md`, following the existing style (in particular, make sure that the change log for this version has the correct version number and date).\n - [ ] Run `make changelog`, to make sure that all the issue references in `CHANGELOG.md` are up to date.\n - [ ] Check that the commit messages in this PR do not contain `[ci skip]`\n - [ ] Run https://github.com/JuliaDocs/Documenter.jl/actions/workflows/regression-tests.yml\n using a `workflow_dispatch` trigger to check for any changes that broke extensions.\n\n## The release\n\n - [ ] After merging the pull request, tag the release. There are two options for this:\n\n 1. [Comment `[at]JuliaRegistrator register` on the GitHub commit.](https://github.com/JuliaRegistries/Registrator.jl#via-the-github-app)\n 2. Use [JuliaHub's package registration feature](https://help.juliahub.com/juliahub/stable/contribute/#registrator) to trigger the registration.\n\n Either of those should automatically publish a new version to the Julia registry.\n - Once registered, the `TagBot.yml` workflow should create a tag, and rebuild the documentation for this tag.\n - These steps can take quite a bit of time (1 hour or more), so don't be surprised if the new documentation takes a while to appear.","category":"page"},{"location":"release-notes/","page":"Release notes","title":"Release notes","text":"EditURL = \"https://github.com/JuliaDocs/Documenter.jl/blob/master/CHANGELOG.md\"","category":"page"},{"location":"release-notes/#Release-notes","page":"Release notes","title":"Release notes","text":"","category":"section"},{"location":"release-notes/","page":"Release notes","title":"Release notes","text":"The format is based on Keep a Changelog, and this project adheres to Semantic Versioning.","category":"page"},{"location":"release-notes/#Unreleased","page":"Release notes","title":"Unreleased","text":"","category":"section"},{"location":"release-notes/#Added","page":"Release notes","title":"Added","text":"","category":"section"},{"location":"release-notes/#Changed","page":"Release notes","title":"Changed","text":"","category":"section"},{"location":"release-notes/#Fixed","page":"Release notes","title":"Fixed","text":"","category":"section"},{"location":"lib/internals/#internal-Documentation","page":"Internals","title":"internal Documentation","text":"","category":"section"},{"location":"lib/internals/","page":"Internals","title":"Internals","text":"Documentation for EpiAware.jl's internal interface.","category":"page"},{"location":"lib/internals/#Contents","page":"Internals","title":"Contents","text":"","category":"section"},{"location":"lib/internals/","page":"Internals","title":"Internals","text":"Pages = [\"internal.md\"]\nDepth = 2:2","category":"page"},{"location":"lib/internals/#Index","page":"Internals","title":"Index","text":"","category":"section"},{"location":"lib/internals/","page":"Internals","title":"Internals","text":"Pages = [\"internals.md\"]","category":"page"},{"location":"lib/internals/#Internal-API","page":"Internals","title":"Internal API","text":"","category":"section"},{"location":"lib/internals/","page":"Internals","title":"Internals","text":"Modules = [EpiAware]\nPublic = false","category":"page"},{"location":"lib/internals/#EpiAware.NegativeBinomialMeanClust-Tuple{Any, Any}","page":"Internals","title":"EpiAware.NegativeBinomialMeanClust","text":"NegativeBinomialMeanClust NegativeBinomialMeanClust(μ, α)\n\nCompute the mean-cluster factor negative binomial distribution.\n\nArguments\n\nμ: The mean of the distribution.\nα: The clustering factor parameter.\n\nReturns\n\nA NegativeBinomial distribution object.\n\n\n\nSignatures\n\nNegativeBinomialMeanClust(\n μ,\n α\n) -> Distributions.NegativeBinomial\n\n\n\n\nMethods\n\nNegativeBinomialMeanClust(μ, α)\n\ndefined at /home/runner/work/Rt-without-renewal/Rt-without-renewal/src/utilities.jl:210.\n\n\n\n\n\n","category":"method"},{"location":"lib/internals/#EpiAware.generate_observation_kernel-Tuple{Any, Any}","page":"Internals","title":"EpiAware.generate_observation_kernel","text":"generateobservationkernel generateobservationkernel(delayint, timehorizon)\n\nGenerate an observation kernel matrix based on the given delay interval and time horizon.\n\nArguments\n\ndelay_int::Vector{Float64}: The delay PMF vector.\ntime_horizon::Int: The number of time steps of the observation period.\n\nReturns\n\nK::SparseMatrixCSC{Float64, Int}: The observation kernel matrix.\n\n\n\nSignatures\n\ngenerate_observation_kernel(delay_int, time_horizon) -> Any\n\n\n\n\nMethods\n\ngenerate_observation_kernel(delay_int, time_horizon)\n\ndefined at /home/runner/work/Rt-without-renewal/Rt-without-renewal/src/utilities.jl:229.\n\n\n\n\n\n","category":"method"},{"location":"lib/internals/#EpiAware.neg_MGF-Tuple{Any, AbstractVector}","page":"Internals","title":"EpiAware.neg_MGF","text":"negMGF negMGF(r, w::AbstractVector)\n\nCompute the negative moment generating function (MGF) for a given rate r and weights w.\n\nArguments\n\nr: The rate parameter.\nw: An abstract vector of weights.\n\nReturns\n\nThe value of the negative MGF.\n\n\n\nSignatures\n\nneg_MGF(r, w::AbstractVector) -> Any\n\n\n\n\nMethods\n\nneg_MGF(r, w)\n\ndefined at /home/runner/work/Rt-without-renewal/Rt-without-renewal/src/utilities.jl:132.\n\n\n\n\n\n","category":"method"},{"location":"contributing/#Contributing","page":"Contributing","title":"Contributing","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"This page details the some of the guidelines that should be followed when contributing to this package. It is adapted from Documenter.jl.","category":"page"},{"location":"contributing/#Branches","page":"Contributing","title":"Branches","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"release-* branches are used for tagged minor versions of this package. This follows the same approach used in the main Julia repository, albeit on a much more modest scale.","category":"page"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"Please open pull requests against the master branch rather than any of the release-* branches whenever possible.","category":"page"},{"location":"contributing/#Backports","page":"Contributing","title":"Backports","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"Bug fixes are backported to the release-* branches using git cherry-pick -x by a EpiAware member and will become available in point releases of that particular minor version of the package.","category":"page"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"Feel free to nominate commits that should be backported by opening an issue. Requests for new point releases to be tagged in METADATA.jl can also be made in the same way.","category":"page"},{"location":"contributing/#release-*-branches","page":"Contributing","title":"release-* branches","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"Each new minor version x.y.0 gets a branch called release-x.y (a protected branch).\nNew versions are usually tagged only from the release-x.y branches.\nFor patch releases, changes get backported to the release-x.y branch via a single PR with the standard name \"Backports for x.y.z\" and label \"Type: Backport\". The PR message links to all the PRs that are providing commits to the backport. The PR gets merged as a merge commit (i.e. not squashed).\nThe old release-* branches may be removed once they have outlived their usefulness.\nPatch version milestones are used to keep track of which PRs get backported etc.","category":"page"},{"location":"contributing/#Style-Guide","page":"Contributing","title":"Style Guide","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"Follow the style of the surrounding text when making changes. When adding new features please try to stick to the following points whenever applicable. This project follows the SciML style guide.","category":"page"},{"location":"contributing/#Tests","page":"Contributing","title":"Tests","text":"","category":"section"},{"location":"contributing/#Unit-tests","page":"Contributing","title":"Unit tests","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"As is conventional for Julia packages, unit tests are located at test/*.jl with the entrypoint test/runtests.jl.","category":"page"},{"location":"contributing/#End-to-end-testing","page":"Contributing","title":"End to end testing","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"Tests that build example package docs from source and inspect the results (end to end tests) are located in /test/examples. The main entry points are test/examples/make.jl for building and test/examples/test.jl for doing some basic checks on the generated outputs.","category":"page"},{"location":"examples/getting_started/","page":"Getting started","title":"Getting started","text":"\n\n\n\n\n\n\n\n

              Getting stated with EpiAware

              This tutorial introduces the basic functionality of EpiAware. EpiAware is a package for making inferences on epidemiological case/determined infection data using a model-based approach.

              EpiAware models

              The models we consider are discrete-time $t = 1,\\dots, T$ with a latent random process, $Z_t$ generating stochasticity in the number of new infections $I_t$ at each time step. Observations are treated as downstream random variables determined by the actual infections and a model of infection to observation delay.

              Mathematical definition

              $$\\begin{align}\nZ_\\cdot &\\sim \\mathcal{P}(\\mathbb{R}^T) | \\theta_Z, \\\\\nI_0 &\\sim f_0(\\theta_I), \\\\\nI_t &\\sim g_I(\\{I_s, Z_s\\}_{s < t}, \\theta_{I}), \\\\\ny_t &\\sim f_O(\\{I_s\\}_{s \\leq t}, \\theta_{O}).\n\\end{align}$$

              Where, $\\mathcal{P}(\\mathbb{R}^T) | \\theta_Z$ is a parametric process on $\\mathbb{R}^T$. $f_0$ and $f_O$ are parametric distributions on, respectively, the initial number of infections and the observed case data conditional on underlying infections. $g_I$ is distribution of new infections conditional on infections and latent process in the past. Note that we assume that new infections are conditional on the strict past, whereas new observations can depend on infections on the same time step.

              Code structure outline

              An EpiAware model in code is created from three modular components:

              • A LatentModel: This defines the model for $Z_\\cdot$. This chooses $\\mathcal{P}(\\mathbb{R}^T) | \\theta_Z$.

              • An EpiModel: This defines a generative process for infections conditional on the latent process. This chooses $f_0(\\theta_I)$, and $g_I(\\{I_s, Z_s\\}_{s < t}, \\theta_{I})$.

              • An ObservationModel: This defines the observation model. This chooses $f_O({I_s}_{s \\leq t}, \\theta_{O})$

              Reproductive number

              EpiAware models do not need to specify a time-varying reproductive number $\\mathcal{R}_t$ to generate $I_\\cdot$, however, this is often a quantity of interest. When not directly used we will typically back-calculate $\\mathcal{R}_t$ from the generated infections:

              $$\\mathcal{R}_t = {I_t \\over \\sum_{s \\geq 1} g_s I_{t-s} }.$$

              Where $g_s$ is a discrete generation interval. For this reason, even when not using a reproductive number approach directly, we ask for a generation interval.

              \n\n
              begin\n    using EpiAware\n    using Turing\n    using Distributions\n    using StatsPlots\n    using Random\n    using DynamicPPL\n    using Statistics\n    using DataFramesMeta\n    using LinearAlgebra\nend
              \n\n\n\n

              Random walk LatentModel

              As an example, we choose the latent process as a random walk with parameters $\\theta_Z$:

              • $Z_0$: Initial position.

              • $\\sigma^2_{Z}$: The step-size variance.

              Conditional on the parameters the random walk is then generated by white noise:

              $$\\begin{align}\nZ_t &= Z_0 + \\sigma_{RW} \\sum_{t= 1}^T \\epsilon_t, \\\\\n\\epsilon_t &\\sim \\mathcal{N}(0,1).\n\\end{align}$$

              In EpiAware we provide a constructor for random walk latent models with priors for $\\theta_Z$. We choose priors,

              $$\\begin{align}\nZ_0 &\\sim \\mathcal{N}(0,1),\\\\\n\\sigma^2_Z &\\sim \\text{HalfNormal}(0.01).\n\\end{align}$$

              \n\n
              rwp = EpiAware.RandomWalk(Normal(),\n    truncated(Normal(0.0, 0.02), 0.0, Inf))
              \n
              RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(Distributions.Normal{Float64}(μ=0.0, σ=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.02); lower=0.0, upper=Inf))
              \n\n\n

              Direct infection EpiModel

              This is a simple model where the unobserved log-infections are directly generated by the latent process $Z$.

              $$\\log I_t = \\log I_0 + Z_t.$$

              As discussed above, we still ask for a defined generation interval, which can be used to calculate $\\mathcal{R}_t$.

              \n\n
              truth_GI = Gamma(2, 5)
              \n
              Distributions.Gamma{Float64}(α=2.0, θ=5.0)
              \n\n\n

              The EpiData constructor performs double interval censoring to convert our continuous estimate of the generation interval into a discretized version. We also implement right truncation using the keyword D_gen.

              \n\n
              model_data = EpiData(truth_GI, D_gen = 10.0)
              \n
              EpiData{Float64, typeof(exp)}([0.056722565637507084, 0.0944813112842689, 0.11668723746973605, 0.12773814268679387, 0.130948497124998, 0.1287976856448402, 0.12312384970427602, 0.11527489380975135, 0.10622581663782862], 9, exp)
              \n\n\n

              We can supply a prior for the initial log_infections.

              \n\n
              log_I0_prior = Normal(log(100.0), 1.0)
              \n
              Distributions.Normal{Float64}(μ=4.605170185988092, σ=1.0)
              \n\n\n

              And construct the EpiModel.

              \n\n
              epi_model = DirectInfections(model_data, log_I0_prior)
              \n
              DirectInfections{Normal{Float64}}(EpiData{Float64, typeof(exp)}([0.056722565637507084, 0.0944813112842689, 0.11668723746973605, 0.12773814268679387, 0.130948497124998, 0.1287976856448402, 0.12312384970427602, 0.11527489380975135, 0.10622581663782862], 9, exp), Distributions.Normal{Float64}(μ=4.605170185988092, σ=1.0))
              \n\n\n

              Delayed Observations ObservationModel

              The observation model is a negative binomial distribution parameterised with mean $\\mu$ and 'successes' parameter $r$. The standard deviation relative to the mean $\\sigma_{\\text{rel}} = \\sigma / \\mu$ for negative binomial observations is,

              $$\\sigma_{\\text{rel}} =(1/\\sqrt{\\mu}) + (1 / \\sqrt{r}).$$

              It is standard to use a half-t distribution for standard deviation priors (e.g. as argued in this paper); we specialise this to a Half-Normal prior and use an a priori assumption that a typical observation fluctuation around the mean (when the mean is $\\sim\\mathcal{O}(10^2)$) would be 10%. This implies a standard deviation prior,

              $$1 / \\sqrt{r} \\sim \\text{HalfNormal}\\Big(0.1 ~\\sqrt{{\\pi \\over 2}}\\Big).$$

              The $\\sqrt{{\\pi \\over 2}}$ factor ensures the correct prior mean (see here).

              The expected observed cases are delayed infections. Delays are implemented as the action of a sparse kernel on the infections $I(t)$.

              $$y_t \\sim \\text{NegBinomial}\\Big(\\mu = \\sum_{s\\geq 0} K[t, t-s] I(s), r\\Big). \\\\$$

              \n\n\n

              We also set up the inference to occur over 100 days.

              \n\n
              time_horizon = 100
              \n
              100
              \n\n\n

              We choose a simple observation model where infections are observed 0, 1, 2, 3 days later with equal probability.

              \n\n
              obs_model = EpiAware.DelayObservations(\n    fill(0.25, 4),\n    time_horizon,\n    truncated(Normal(0, 0.1 * sqrt(pi) / sqrt(2)), 0.0, Inf)\n)
              \n
              DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(sparse(Integer[1, 2, 3, 4, 2, 3, 4, 5, 3, 4  …  97, 98, 99, 100, 98, 99, 100, 99, 100, 100], Integer[1, 1, 1, 1, 2, 2, 2, 2, 3, 3  …  97, 97, 97, 97, 98, 98, 98, 99, 99, 100], [0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25  …  0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25], 100, 100), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.12533141373155002); lower=0.0, upper=Inf))
              \n\n\n

              Generate cases from the EpiAware model

              Having chosen an EpiModel, LatentModel and ObservationModel, we can implement the model as a Turing model using make_epi_aware.

              By giving missing to the first argument, we indicate that case data will be generated from the model rather than treated as fixed.

              \n\n
              full_epi_aware_mdl = make_epi_aware(missing, time_horizon;\n    epi_model = epi_model,\n    latent_model = rwp,\n    observation_model = obs_model)
              \n
              Model{typeof(make_epi_aware), (:y_t, :time_steps), (:epi_model, :latent_model, :observation_model, :pos_shift), (:y_t,), Tuple{Missing, Int64}, Tuple{DirectInfections{Normal{Float64}}, RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, Float64}, DefaultContext}(EpiAware.make_epi_aware, (y_t = missing, time_steps = 100), (epi_model = DirectInfections{Normal{Float64}}(EpiData{Float64, typeof(exp)}([0.056722565637507084, 0.0944813112842689, 0.11668723746973605, 0.12773814268679387, 0.130948497124998, 0.1287976856448402, 0.12312384970427602, 0.11527489380975135, 0.10622581663782862], 9, exp), Distributions.Normal{Float64}(μ=4.605170185988092, σ=1.0)), latent_model = RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(Distributions.Normal{Float64}(μ=0.0, σ=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.02); lower=0.0, upper=Inf)), observation_model = DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(sparse(Integer[1, 2, 3, 4, 2, 3, 4, 5, 3, 4  …  97, 98, 99, 100, 98, 99, 100, 99, 100, 100], Integer[1, 1, 1, 1, 2, 2, 2, 2, 3, 3  …  97, 97, 97, 97, 98, 98, 98, 99, 99, 100], [0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25  …  0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25], 100, 100), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.12533141373155002); lower=0.0, upper=Inf)), pos_shift = 1.0e-6), DefaultContext())
              \n\n\n

              We choose some fixed parameters:

              • Initial incidence is 100.

              • In the direct infection model, the initial incidence and in the initial value of the random walk form a non-identifiable pair. Therefore, we fix $Z_0 = 0$.

              \n\n
              fixed_parameters = (rw_init = 0.0, init_incidence = log(100.0))
              \n
              (rw_init = 0.0, init_incidence = 4.605170185988092)
              \n\n\n

              We fix these parameters using fix, and generate a random epidemic.

              \n\n
              cond_generative_model = fix(full_epi_aware_mdl, fixed_parameters)
              \n
              Model{typeof(make_epi_aware), (:y_t, :time_steps), (:epi_model, :latent_model, :observation_model, :pos_shift), (:y_t,), Tuple{Missing, Int64}, Tuple{DirectInfections{Normal{Float64}}, RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, Float64}, DynamicPPL.FixedContext{@NamedTuple{rw_init::Float64, init_incidence::Float64}, DefaultContext}}(EpiAware.make_epi_aware, (y_t = missing, time_steps = 100), (epi_model = DirectInfections{Normal{Float64}}(EpiData{Float64, typeof(exp)}([0.056722565637507084, 0.0944813112842689, 0.11668723746973605, 0.12773814268679387, 0.130948497124998, 0.1287976856448402, 0.12312384970427602, 0.11527489380975135, 0.10622581663782862], 9, exp), Distributions.Normal{Float64}(μ=4.605170185988092, σ=1.0)), latent_model = RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(Distributions.Normal{Float64}(μ=0.0, σ=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.02); lower=0.0, upper=Inf)), observation_model = DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(sparse(Integer[1, 2, 3, 4, 2, 3, 4, 5, 3, 4  …  97, 98, 99, 100, 98, 99, 100, 99, 100, 100], Integer[1, 1, 1, 1, 2, 2, 2, 2, 3, 3  …  97, 97, 97, 97, 98, 98, 98, 99, 99, 100], [0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25  …  0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25], 100, 100), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.12533141373155002); lower=0.0, upper=Inf)), pos_shift = 1.0e-6), FixedContext((rw_init = 0.0, init_incidence = 4.605170185988092), DynamicPPL.DefaultContext()))
              \n\n
              random_epidemic = rand(cond_generative_model)
              \n
              (ϵ_t = [0.7541364299786495, 0.6228295155052502, -0.08437805815079774, -1.6027945048950885, 0.5288638281567213, -1.6896912280141976, -1.6475552716129747, -1.4734033626550678, 0.2292953790785949, -1.4453965284513017  …  -0.5319988979469747, -0.1702093780381867, 0.7423616167719554, 0.5084938565160623, -0.7552829192371855, -0.12043077985555731, -1.7538972893049496, 0.5572895108833601, -0.7266998274413684, 1.0439092851436422], σ²_RW = 0.030509604813344984, neg_bin_cluster_factor = 0.003822359489741684, var\"y_t[1]\" = 27, var\"y_t[2]\" = 63, var\"y_t[3]\" = 75, var\"y_t[4]\" = 121, var\"y_t[5]\" = 126, var\"y_t[6]\" = 99, var\"y_t[7]\" = 84, var\"y_t[8]\" = 59, var\"y_t[9]\" = 62, var\"y_t[10]\" = 42, var\"y_t[11]\" = 60, var\"y_t[12]\" = 36, var\"y_t[13]\" = 52, var\"y_t[14]\" = 47, var\"y_t[15]\" = 47, var\"y_t[16]\" = 47, var\"y_t[17]\" = 35, var\"y_t[18]\" = 58, var\"y_t[19]\" = 55, var\"y_t[20]\" = 62, var\"y_t[21]\" = 66, var\"y_t[22]\" = 88, var\"y_t[23]\" = 73, var\"y_t[24]\" = 78, var\"y_t[25]\" = 73, var\"y_t[26]\" = 68, var\"y_t[27]\" = 81, var\"y_t[28]\" = 64, var\"y_t[29]\" = 78, var\"y_t[30]\" = 68, var\"y_t[31]\" = 71, var\"y_t[32]\" = 80, var\"y_t[33]\" = 95, var\"y_t[34]\" = 74, var\"y_t[35]\" = 79, var\"y_t[36]\" = 86, var\"y_t[37]\" = 101, var\"y_t[38]\" = 94, var\"y_t[39]\" = 89, var\"y_t[40]\" = 93, var\"y_t[41]\" = 77, var\"y_t[42]\" = 90, var\"y_t[43]\" = 98, var\"y_t[44]\" = 108, var\"y_t[45]\" = 129, var\"y_t[46]\" = 130, var\"y_t[47]\" = 115, var\"y_t[48]\" = 126, var\"y_t[49]\" = 99, var\"y_t[50]\" = 121, var\"y_t[51]\" = 96, var\"y_t[52]\" = 111, var\"y_t[53]\" = 82, var\"y_t[54]\" = 80, var\"y_t[55]\" = 67, var\"y_t[56]\" = 80, var\"y_t[57]\" = 91, var\"y_t[58]\" = 87, var\"y_t[59]\" = 103, var\"y_t[60]\" = 118, var\"y_t[61]\" = 129, var\"y_t[62]\" = 138, var\"y_t[63]\" = 133, var\"y_t[64]\" = 150, var\"y_t[65]\" = 163, var\"y_t[66]\" = 190, var\"y_t[67]\" = 219, var\"y_t[68]\" = 216, var\"y_t[69]\" = 225, var\"y_t[70]\" = 335, var\"y_t[71]\" = 362, var\"y_t[72]\" = 447, var\"y_t[73]\" = 469, var\"y_t[74]\" = 460, var\"y_t[75]\" = 465, var\"y_t[76]\" = 458, var\"y_t[77]\" = 457, var\"y_t[78]\" = 470, var\"y_t[79]\" = 473, var\"y_t[80]\" = 516, var\"y_t[81]\" = 631, var\"y_t[82]\" = 606, var\"y_t[83]\" = 593, var\"y_t[84]\" = 563, var\"y_t[85]\" = 505, var\"y_t[86]\" = 451, var\"y_t[87]\" = 453, var\"y_t[88]\" = 427, var\"y_t[89]\" = 352, var\"y_t[90]\" = 310, var\"y_t[91]\" = 243, var\"y_t[92]\" = 277, var\"y_t[93]\" = 247, var\"y_t[94]\" = 238, var\"y_t[95]\" = 248, var\"y_t[96]\" = 251, var\"y_t[97]\" = 229, var\"y_t[98]\" = 234, var\"y_t[99]\" = 187, var\"y_t[100]\" = 193)
              \n\n
              true_infections = generated_quantities(cond_generative_model, random_epidemic).I_t
              \n
              100-element Vector{Float64}:\n 114.07945592233989\n 127.19035360323254\n 125.32952857511957\n  94.72579974976007\n 103.89315740245587\n  77.3410367647116\n  58.00018680718289\n   ⋮\n 248.0613826703851\n 242.89775255149056\n 178.80348692965416\n 197.0838192629199\n 173.59004678270438\n 208.31181943121163
              \n\n
              generated_obs = generated_quantities(cond_generative_model, random_epidemic).generated_y_t
              \n
              100-element Vector{Int64}:\n  27\n  63\n  75\n 121\n 126\n  99\n  84\n   ⋮\n 248\n 251\n 229\n 234\n 187\n 193
              \n\n\n\n\n","category":"page"},{"location":"examples/getting_started/#Inference","page":"Getting started","title":"Inference","text":"","category":"section"},{"location":"examples/getting_started/","page":"Getting started","title":"Getting started","text":"
              \n

              Fixing $Z_0 = 0$ for the random walk was based on inference principles; in this model $Z_0$ and $\\log I_0$ are non-identifiable.

              However, we now treat the generated data as truth_data and make inference without fixing any other parameters.

              We do the inference by MCMC/NUTS using the Turing NUTS sampler with default warm-up steps.

              \n\n
              truth_data = generated_obs
              \n
              100-element Vector{Int64}:\n  27\n  63\n  75\n 121\n 126\n  99\n  84\n   ⋮\n 248\n 251\n 229\n 234\n 187\n 193
              \n\n\n

              The observation model supports partially complete data. To test this we set some of the generated observations to be missing.

              \n\n
              let\n    truth_data = Union{Int, Missing}[truth_data...]\n    truth_data[vcat([3, 5], 10:20)] .= missing\nend
              \n
              13-element view(::Vector{Union{Missing, Int64}}, [3, 5, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]) with eltype Union{Missing, Int64}:\n missing\n missing\n missing\n missing\n missing\n missing\n missing\n missing\n missing\n missing\n missing\n missing\n missing
              \n\n\n

              We now make the model but fixing the initial condition of the random walk to be 0.

              \n\n
              \ninference_mdl = fix(\n    make_epi_aware(truth_data, time_horizon; epi_model = epi_model,\n        latent_model = rwp, observation_model = obs_model),\n    (rw_init = 0.0,)\n)
              \n
              Model{typeof(make_epi_aware), (:y_t, :time_steps), (:epi_model, :latent_model, :observation_model, :pos_shift), (), Tuple{Vector{Int64}, Int64}, Tuple{DirectInfections{Normal{Float64}}, RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, Float64}, DynamicPPL.FixedContext{@NamedTuple{rw_init::Float64}, DefaultContext}}(EpiAware.make_epi_aware, (y_t = [27, 63, 75, 121, 126, 99, 84, 59, 62, 42  …  243, 277, 247, 238, 248, 251, 229, 234, 187, 193], time_steps = 100), (epi_model = DirectInfections{Normal{Float64}}(EpiData{Float64, typeof(exp)}([0.056722565637507084, 0.0944813112842689, 0.11668723746973605, 0.12773814268679387, 0.130948497124998, 0.1287976856448402, 0.12312384970427602, 0.11527489380975135, 0.10622581663782862], 9, exp), Distributions.Normal{Float64}(μ=4.605170185988092, σ=1.0)), latent_model = RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(Distributions.Normal{Float64}(μ=0.0, σ=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.02); lower=0.0, upper=Inf)), observation_model = DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(sparse(Integer[1, 2, 3, 4, 2, 3, 4, 5, 3, 4  …  97, 98, 99, 100, 98, 99, 100, 99, 100, 100], Integer[1, 1, 1, 1, 2, 2, 2, 2, 3, 3  …  97, 97, 97, 97, 98, 98, 98, 99, 99, 100], [0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25  …  0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25], 100, 100), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.12533141373155002); lower=0.0, upper=Inf)), pos_shift = 1.0e-6), FixedContext((rw_init = 0.0,), DynamicPPL.DefaultContext()))
              \n\n
              chn = sample(inference_mdl,\n    NUTS(; adtype = AutoReverseDiff(true)),\n    MCMCThreads(),\n    250,\n    4;\n    drop_warmup = true)
              \n
              iterationchainϵ_t[1]ϵ_t[2]ϵ_t[3]ϵ_t[4]ϵ_t[5]ϵ_t[6]...
              11261-1.236280.23389-0.485811-0.494185-0.2627-2.06568
              212711.25850.289792-0.0274798-0.715777-0.759128-1.7599
              31281-2.127861.34623-1.200970.124472-2.265070.198596
              41291-2.086831.00838-1.834050.228425-2.24750.219736
              513010.884139-0.0395667-0.8535551.94183-1.763630.482344
              61311-0.5916391.949591.0209-1.98888-1.21572-2.21748
              71321-0.2620560.8361030.315522-0.723914-1.51473-0.888812
              81331-0.6898560.812266-1.08283-0.331231-1.17112-0.520497
              913411.033671.28017-0.821226-1.472-0.90994-0.244272
              101351-1.893661.18913-0.578158-1.30672-1.0126-2.69586
              ...
              \n\n\n

              Predictive plotting

              We can spaghetti plot generated case data from the version of the model which hasn't conditioned on case data using posterior parameters inferred from the version conditioned on observed data. This is known as posterior predictive checking, and is a useful diagnostic tool for Bayesian inference (see here).

              Because we are using synthetic data we can also plot the model predictions for the unobserved infections and check that (at least in this example) we were able to capture some unobserved/latent variables in the process accurate.

              \n\n
              let\n    post_check_mdl = fix(full_epi_aware_mdl, (rw_init = 0.0,))\n    post_check_y_t = mapreduce(hcat, generated_quantities(post_check_mdl, chn)) do gen\n        gen.generated_y_t\n    end\n\n    predicted_I_t = mapreduce(hcat, generated_quantities(inference_mdl, chn)) do gen\n        gen.I_t\n    end\n\n    p1 = plot(post_check_y_t, c = :grey, alpha = 0.05, lab = \"\")\n    scatter!(p1, truth_data,\n        lab = \"Observed cases\",\n        xlabel = \"Time\",\n        ylabel = \"Cases\",\n        title = \"Post. predictive checking: cases\",\n        ylims = (-0.5, maximum(truth_data) * 1.5),\n        c = :green)\n\n    p2 = plot(predicted_I_t, c = :grey, alpha = 0.05, lab = \"\")\n    scatter!(p2, true_infections,\n        lab = \"Actual infections\",\n        xlabel = \"Time\",\n        ylabel = \"Unobserved Infections\",\n        title = \"Post. predictions: infections\",\n        ylims = (-0.5, maximum(true_infections) * 1.5),\n        c = :red)\n\n    plot(p1, p2,\n        layout = (1, 2),\n        size = (700, 400))\nend
              \n\n\n\n

              As well as checking the posterior predictions for latent infections, we can also check how well inference recovered unknown parameters, such as the random walk variance or the cluster factor of the negative binomial observations.

              \n\n
              let\n    parameters_to_plot = (:σ²_RW, :neg_bin_cluster_factor)\n\n    plts = map(parameters_to_plot) do name\n        var_samples = chn[name] |> vec\n        histogram(var_samples,\n            bins = 50,\n            norm = :pdf,\n            lw = 0,\n            fillalpha = 0.5,\n            lab = \"MCMC\")\n        vline!([getfield(random_epidemic, name)], lab = \"True value\")\n        title!(string(name))\n    end\n    plot(plts..., layout = (2, 1))\nend
              \n\n\n","category":"page"},{"location":"examples/getting_started/#Reproductive-number-back-calculation","page":"Getting started","title":"Reproductive number back-calculation","text":"","category":"section"},{"location":"examples/getting_started/","page":"Getting started","title":"Getting started","text":"
              \n

              As mentioned at the top, we don't directly use the concept of reproductive numbers in this note. However, we can back-calculate the implied $\\mathcal{R}(t)$ values, conditional on the specified generation interval being correct.

              Here we spaghetti plot posterior sampled time-varying reproductive numbers against the actual.

              \n\n
              let\n    n = epi_model.data.len_gen_int\n    Rt_denom = [dot(reverse(epi_model.data.gen_int), true_infections[(t - n):(t - 1)])\n                for t in (n + 1):length(true_infections)]\n    true_Rt = true_infections[(n + 1):end] ./ Rt_denom\n\n    predicted_Rt = mapreduce(hcat, generated_quantities(inference_mdl, chn)) do gen\n        _It = gen.I_t\n        _Rt_denom = [dot(reverse(epi_model.data.gen_int), _It[(t - n):(t - 1)])\n                     for t in (n + 1):length(_It)]\n        Rt = _It[(n + 1):end] ./ _Rt_denom\n    end\n\n    plt = plot((n + 1):time_horizon, predicted_Rt, c = :grey, alpha = 0.05, lab = \"\")\n    plot!(plt, (n + 1):time_horizon, true_Rt,\n        lab = \"true Rt\",\n        xlabel = \"Time\",\n        ylabel = \"Rt\",\n        title = \"Post. predictions: reproductive number\",\n        c = :red,\n        lw = 2)\nend
              \n\n\n","category":"page"},{"location":"examples/getting_started/","page":"Getting started","title":"Getting started","text":"EditURL = \"https://github.com/CDCgov/Rt-without-renewal/blob/main/docs/src/examples/getting_started.jl\"","category":"page"},{"location":"lib/public/#Public-Documentation","page":"Public API","title":"Public Documentation","text":"","category":"section"},{"location":"lib/public/","page":"Public API","title":"Public API","text":"Documentation for EpiAware.jl's public interface.","category":"page"},{"location":"lib/public/","page":"Public API","title":"Public API","text":"See the Internals section of the manual for internal package docs covering all submodules.","category":"page"},{"location":"lib/public/#Contents","page":"Public API","title":"Contents","text":"","category":"section"},{"location":"lib/public/","page":"Public API","title":"Public API","text":"Pages = [\"public.md\"]\nDepth = 2:2","category":"page"},{"location":"lib/public/#Index","page":"Public API","title":"Index","text":"","category":"section"},{"location":"lib/public/","page":"Public API","title":"Public API","text":"Pages = [\"public.md\"]","category":"page"},{"location":"lib/public/#Public-API","page":"Public API","title":"Public API","text":"","category":"section"},{"location":"#EpiAware.jl","page":"EpiAware.jl: Real-time epidemic monitoring","title":"EpiAware.jl","text":"","category":"section"},{"location":"","page":"EpiAware.jl: Real-time epidemic monitoring","title":"EpiAware.jl: Real-time epidemic monitoring","text":"Infectious disease situational awareness modelling toolkit for Julia.","category":"page"},{"location":"","page":"EpiAware.jl: Real-time epidemic monitoring","title":"EpiAware.jl: Real-time epidemic monitoring","text":"A package for building and fitting situational awareness models for infectious diseases. The package is designed to be flexible and extensible, and to provide a consistent interface for fitting and simulating models.","category":"page"},{"location":"#Package-Features","page":"EpiAware.jl: Real-time epidemic monitoring","title":"Package Features","text":"","category":"section"},{"location":"","page":"EpiAware.jl: Real-time epidemic monitoring","title":"EpiAware.jl: Real-time epidemic monitoring","text":"Flexible: The package is designed to be flexible and extensible, and to provide a consistent interface for fitting and simulating models.\nModular: The package is designed to be modular, with a clear separation between the model and the data.\nExtensible: The package is designed to be extensible, with a clear separation between the model and the data.\nConsistent: The package is designed to provide a consistent interface for fitting and simulating models.\nEfficient: The package is designed to be efficient, with a clear separation between the model and the data.","category":"page"},{"location":"","page":"EpiAware.jl: Real-time epidemic monitoring","title":"EpiAware.jl: Real-time epidemic monitoring","text":"See the Index for the complete list of documented functions and types.","category":"page"},{"location":"#Manual-Outline","page":"EpiAware.jl: Real-time epidemic monitoring","title":"Manual Outline","text":"","category":"section"},{"location":"","page":"EpiAware.jl: Real-time epidemic monitoring","title":"EpiAware.jl: Real-time epidemic monitoring","text":"Pages = [\n \"man/contributing.md\",\n]\nDepth = 1","category":"page"},{"location":"#Library-Outline","page":"EpiAware.jl: Real-time epidemic monitoring","title":"Library Outline","text":"","category":"section"},{"location":"","page":"EpiAware.jl: Real-time epidemic monitoring","title":"EpiAware.jl: Real-time epidemic monitoring","text":"Pages = [\"lib/public.md\", \"lib/internals.md\"]","category":"page"},{"location":"#main-index","page":"EpiAware.jl: Real-time epidemic monitoring","title":"Index","text":"","category":"section"},{"location":"","page":"EpiAware.jl: Real-time epidemic monitoring","title":"EpiAware.jl: Real-time epidemic monitoring","text":"Pages = [\"lib/public.md\"]","category":"page"}] +[{"location":"checklist/#Checklists","page":"Release checklist","title":"Checklists","text":"","category":"section"},{"location":"checklist/","page":"Release checklist","title":"Release checklist","text":"The purpose of this page is to collate a series of checklists for commonly performed changes to the source code of EpiAware. It has been adapted from Documenter.jl.","category":"page"},{"location":"checklist/","page":"Release checklist","title":"Release checklist","text":"In each case, copy the checklist into the description of the pull request.","category":"page"},{"location":"checklist/#Making-a-release","page":"Release checklist","title":"Making a release","text":"","category":"section"},{"location":"checklist/","page":"Release checklist","title":"Release checklist","text":"In preparation for a release, use the following checklist. These steps should be performed on a branch with an open pull request, either for a topic branch, or for a new branch release-1.y.z (\"Release version 1.y.z\") if multiple changes have accumulated on the master branch since the last release.","category":"page"},{"location":"checklist/","page":"Release checklist","title":"Release checklist","text":"## Pre-release\n\n - [ ] Change the version number in `Project.toml`\n * If the release is breaking, increment MAJOR\n * If the release adds a new user-visible feature, increment MINOR\n * Otherwise (bug-fixes, documentation improvements), increment PATCH\n - [ ] Update `CHANGELOG.md`, following the existing style (in particular, make sure that the change log for this version has the correct version number and date).\n - [ ] Run `make changelog`, to make sure that all the issue references in `CHANGELOG.md` are up to date.\n - [ ] Check that the commit messages in this PR do not contain `[ci skip]`\n - [ ] Run https://github.com/JuliaDocs/Documenter.jl/actions/workflows/regression-tests.yml\n using a `workflow_dispatch` trigger to check for any changes that broke extensions.\n\n## The release\n\n - [ ] After merging the pull request, tag the release. There are two options for this:\n\n 1. [Comment `[at]JuliaRegistrator register` on the GitHub commit.](https://github.com/JuliaRegistries/Registrator.jl#via-the-github-app)\n 2. Use [JuliaHub's package registration feature](https://help.juliahub.com/juliahub/stable/contribute/#registrator) to trigger the registration.\n\n Either of those should automatically publish a new version to the Julia registry.\n - Once registered, the `TagBot.yml` workflow should create a tag, and rebuild the documentation for this tag.\n - These steps can take quite a bit of time (1 hour or more), so don't be surprised if the new documentation takes a while to appear.","category":"page"},{"location":"release-notes/","page":"Release notes","title":"Release notes","text":"EditURL = \"https://github.com/JuliaDocs/Documenter.jl/blob/master/CHANGELOG.md\"","category":"page"},{"location":"release-notes/#Release-notes","page":"Release notes","title":"Release notes","text":"","category":"section"},{"location":"release-notes/","page":"Release notes","title":"Release notes","text":"The format is based on Keep a Changelog, and this project adheres to Semantic Versioning.","category":"page"},{"location":"release-notes/#Unreleased","page":"Release notes","title":"Unreleased","text":"","category":"section"},{"location":"release-notes/#Added","page":"Release notes","title":"Added","text":"","category":"section"},{"location":"release-notes/#Changed","page":"Release notes","title":"Changed","text":"","category":"section"},{"location":"release-notes/#Fixed","page":"Release notes","title":"Fixed","text":"","category":"section"},{"location":"lib/internals/#internal-Documentation","page":"Internals","title":"internal Documentation","text":"","category":"section"},{"location":"lib/internals/","page":"Internals","title":"Internals","text":"Documentation for EpiAware.jl's internal interface.","category":"page"},{"location":"lib/internals/#Contents","page":"Internals","title":"Contents","text":"","category":"section"},{"location":"lib/internals/","page":"Internals","title":"Internals","text":"Pages = [\"internal.md\"]\nDepth = 2:2","category":"page"},{"location":"lib/internals/#Index","page":"Internals","title":"Index","text":"","category":"section"},{"location":"lib/internals/","page":"Internals","title":"Internals","text":"Pages = [\"internals.md\"]","category":"page"},{"location":"lib/internals/#Internal-API","page":"Internals","title":"Internal API","text":"","category":"section"},{"location":"lib/internals/","page":"Internals","title":"Internals","text":"Modules = [EpiAware]\nPublic = false","category":"page"},{"location":"lib/internals/#EpiAware.NegativeBinomialMeanClust-Tuple{Any, Any}","page":"Internals","title":"EpiAware.NegativeBinomialMeanClust","text":"NegativeBinomialMeanClust NegativeBinomialMeanClust(μ, α)\n\nCompute the mean-cluster factor negative binomial distribution.\n\nArguments\n\nμ: The mean of the distribution.\nα: The clustering factor parameter.\n\nReturns\n\nA NegativeBinomial distribution object.\n\n\n\nSignatures\n\nNegativeBinomialMeanClust(\n μ,\n α\n) -> Distributions.NegativeBinomial\n\n\n\n\nMethods\n\nNegativeBinomialMeanClust(μ, α)\n\ndefined at /home/runner/work/Rt-without-renewal/Rt-without-renewal/src/utilities.jl:210.\n\n\n\n\n\n","category":"method"},{"location":"lib/internals/#EpiAware._continue_manypathfinder!-Tuple{Any, DynamicPPL.Model}","page":"Internals","title":"EpiAware._continue_manypathfinder!","text":"continuemanypathfinder! Continue running the pathfinder algorithm until a pathfinder succeeds or the maximum number of tries is reached.\n\nArguments\n\npfs: An array of pathfinder objects.\nmdl::DynamicPPL.Model: The model to perform inference on.\nmax_tries: The maximum number of tries to run the pathfinder algorithm. Default is Inf.\nnruns: The number of times to run the pathfinder function.\nkwargs...: Additional keyword arguments passed to pathfinder.\n\nReturns\n\npfs: The updated array of pathfinder objects.\n\n\n\nSignatures\n\n_continue_manypathfinder!(\n pfs,\n mdl::DynamicPPL.Model;\n max_tries,\n nruns,\n kwargs...\n)\n\n\n\n\nMethods\n\n_continue_manypathfinder!(\n pfs,\n mdl;\n max_tries,\n nruns,\n kwargs...\n)\n\ndefined at /home/runner/work/Rt-without-renewal/Rt-without-renewal/src/inference-methods.jl:42.\n\n\n\n\n\n","category":"method"},{"location":"lib/internals/#EpiAware._get_best_elbo_pathfinder-Tuple{Any}","page":"Internals","title":"EpiAware._get_best_elbo_pathfinder","text":"getbestelbopathfinder Selects the pathfinder with the highest ELBO estimate from a list of pathfinders.\n\nArguments\n\npfs: A list of pathfinders results or Symbol values indicating failure.\n\nReturns\n\nThe pathfinder with the highest ELBO estimate.\n\n\n\nSignatures\n\n_get_best_elbo_pathfinder(pfs) -> Any\n\n\n\n\nMethods\n\n_get_best_elbo_pathfinder(pfs)\n\ndefined at /home/runner/work/Rt-without-renewal/Rt-without-renewal/src/inference-methods.jl:72.\n\n\n\n\n\n","category":"method"},{"location":"lib/internals/#EpiAware._run_manypathfinder-Tuple{DynamicPPL.Model}","page":"Internals","title":"EpiAware._run_manypathfinder","text":"runmanypathfinder\n\nRun pathfinder multiple times and store the results in an array. Fails safely.\n\nArguments\n\nmdl::DynamicPPL.Model: The Turing model to be used for inference.\nnruns: The number of times to run the pathfinder function.\nkwargs...: Additional keyword arguments passed to pathfinder.\n\nReturns\n\nAn array of PathfinderResult objects or Symbol values indicating success or failure.\n\n\n\nSignatures\n\n_run_manypathfinder(mdl::DynamicPPL.Model; nruns, kwargs...)\n\n\n\n\nMethods\n\n_run_manypathfinder(mdl; nruns, kwargs...)\n\ndefined at /home/runner/work/Rt-without-renewal/Rt-without-renewal/src/inference-methods.jl:13.\n\n\n\n\n\n","category":"method"},{"location":"lib/internals/#EpiAware.generate_observation_kernel-Tuple{Any, Any}","page":"Internals","title":"EpiAware.generate_observation_kernel","text":"generateobservationkernel generateobservationkernel(delayint, timehorizon)\n\nGenerate an observation kernel matrix based on the given delay interval and time horizon.\n\nArguments\n\ndelay_int::Vector{Float64}: The delay PMF vector.\ntime_horizon::Int: The number of time steps of the observation period.\n\nReturns\n\nK::SparseMatrixCSC{Float64, Int}: The observation kernel matrix.\n\n\n\nSignatures\n\ngenerate_observation_kernel(delay_int, time_horizon) -> Any\n\n\n\n\nMethods\n\ngenerate_observation_kernel(delay_int, time_horizon)\n\ndefined at /home/runner/work/Rt-without-renewal/Rt-without-renewal/src/utilities.jl:229.\n\n\n\n\n\n","category":"method"},{"location":"lib/internals/#EpiAware.neg_MGF-Tuple{Any, AbstractVector}","page":"Internals","title":"EpiAware.neg_MGF","text":"negMGF negMGF(r, w::AbstractVector)\n\nCompute the negative moment generating function (MGF) for a given rate r and weights w.\n\nArguments\n\nr: The rate parameter.\nw: An abstract vector of weights.\n\nReturns\n\nThe value of the negative MGF.\n\n\n\nSignatures\n\nneg_MGF(r, w::AbstractVector) -> Any\n\n\n\n\nMethods\n\nneg_MGF(r, w)\n\ndefined at /home/runner/work/Rt-without-renewal/Rt-without-renewal/src/utilities.jl:132.\n\n\n\n\n\n","category":"method"},{"location":"contributing/#Contributing","page":"Contributing","title":"Contributing","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"This page details the some of the guidelines that should be followed when contributing to this package. It is adapted from Documenter.jl.","category":"page"},{"location":"contributing/#Branches","page":"Contributing","title":"Branches","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"release-* branches are used for tagged minor versions of this package. This follows the same approach used in the main Julia repository, albeit on a much more modest scale.","category":"page"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"Please open pull requests against the master branch rather than any of the release-* branches whenever possible.","category":"page"},{"location":"contributing/#Backports","page":"Contributing","title":"Backports","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"Bug fixes are backported to the release-* branches using git cherry-pick -x by a EpiAware member and will become available in point releases of that particular minor version of the package.","category":"page"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"Feel free to nominate commits that should be backported by opening an issue. Requests for new point releases to be tagged in METADATA.jl can also be made in the same way.","category":"page"},{"location":"contributing/#release-*-branches","page":"Contributing","title":"release-* branches","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"Each new minor version x.y.0 gets a branch called release-x.y (a protected branch).\nNew versions are usually tagged only from the release-x.y branches.\nFor patch releases, changes get backported to the release-x.y branch via a single PR with the standard name \"Backports for x.y.z\" and label \"Type: Backport\". The PR message links to all the PRs that are providing commits to the backport. The PR gets merged as a merge commit (i.e. not squashed).\nThe old release-* branches may be removed once they have outlived their usefulness.\nPatch version milestones are used to keep track of which PRs get backported etc.","category":"page"},{"location":"contributing/#Style-Guide","page":"Contributing","title":"Style Guide","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"Follow the style of the surrounding text when making changes. When adding new features please try to stick to the following points whenever applicable. This project follows the SciML style guide.","category":"page"},{"location":"contributing/#Tests","page":"Contributing","title":"Tests","text":"","category":"section"},{"location":"contributing/#Unit-tests","page":"Contributing","title":"Unit tests","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"As is conventional for Julia packages, unit tests are located at test/*.jl with the entrypoint test/runtests.jl.","category":"page"},{"location":"contributing/#End-to-end-testing","page":"Contributing","title":"End to end testing","text":"","category":"section"},{"location":"contributing/","page":"Contributing","title":"Contributing","text":"Tests that build example package docs from source and inspect the results (end to end tests) are located in /test/examples. The main entry points are test/examples/make.jl for building and test/examples/test.jl for doing some basic checks on the generated outputs.","category":"page"},{"location":"examples/getting_started/","page":"Getting started","title":"Getting started","text":"\n\n\n\n\n\n\n\n

              Getting stated with EpiAware

              This tutorial introduces the basic functionality of EpiAware. EpiAware is a package for making inferences on epidemiological case/determined infection data using a model-based approach.

              EpiAware models

              The models we consider are discrete-time $t = 1,\\dots, T$ with a latent random process, $Z_t$ generating stochasticity in the number of new infections $I_t$ at each time step. Observations are treated as downstream random variables determined by the actual infections and a model of infection to observation delay.

              Mathematical definition

              $$\\begin{align}\nZ_\\cdot &\\sim \\mathcal{P}(\\mathbb{R}^T) | \\theta_Z, \\\\\nI_0 &\\sim f_0(\\theta_I), \\\\\nI_t &\\sim g_I(\\{I_s, Z_s\\}_{s < t}, \\theta_{I}), \\\\\ny_t &\\sim f_O(\\{I_s\\}_{s \\leq t}, \\theta_{O}).\n\\end{align}$$

              Where, $\\mathcal{P}(\\mathbb{R}^T) | \\theta_Z$ is a parametric process on $\\mathbb{R}^T$. $f_0$ and $f_O$ are parametric distributions on, respectively, the initial number of infections and the observed case data conditional on underlying infections. $g_I$ is distribution of new infections conditional on infections and latent process in the past. Note that we assume that new infections are conditional on the strict past, whereas new observations can depend on infections on the same time step.

              Code structure outline

              An EpiAware model in code is created from three modular components:

              • A LatentModel: This defines the model for $Z_\\cdot$. This chooses $\\mathcal{P}(\\mathbb{R}^T) | \\theta_Z$.

              • An EpiModel: This defines a generative process for infections conditional on the latent process. This chooses $f_0(\\theta_I)$, and $g_I(\\{I_s, Z_s\\}_{s < t}, \\theta_{I})$.

              • An ObservationModel: This defines the observation model. This chooses $f_O({I_s}_{s \\leq t}, \\theta_{O})$

              Reproductive number

              EpiAware models do not need to specify a time-varying reproductive number $\\mathcal{R}_t$ to generate $I_\\cdot$, however, this is often a quantity of interest. When not directly used we will typically back-calculate $\\mathcal{R}_t$ from the generated infections:

              $$\\mathcal{R}_t = {I_t \\over \\sum_{s \\geq 1} g_s I_{t-s} }.$$

              Where $g_s$ is a discrete generation interval. For this reason, even when not using a reproductive number approach directly, we ask for a generation interval.

              \n\n
              begin\n    using EpiAware\n    using Turing\n    using Distributions\n    using StatsPlots\n    using Random\n    using DynamicPPL\n    using Statistics\n    using DataFramesMeta\n    using LinearAlgebra\n    using Transducers\nend
              \n\n\n\n

              Random walk LatentModel

              As an example, we choose the latent process as a random walk with parameters $\\theta_Z$:

              • $Z_0$: Initial position.

              • $\\sigma^2_{Z}$: The step-size variance.

              Conditional on the parameters the random walk is then generated by white noise:

              $$\\begin{align}\nZ_t &= Z_0 + \\sigma_{RW} \\sum_{t= 1}^T \\epsilon_t, \\\\\n\\epsilon_t &\\sim \\mathcal{N}(0,1).\n\\end{align}$$

              In EpiAware we provide a constructor for random walk latent models with priors for $\\theta_Z$. We choose priors,

              $$\\begin{align}\nZ_0 &\\sim \\mathcal{N}(0,1),\\\\\n\\sigma^2_Z &\\sim \\text{HalfNormal}(0.01).\n\\end{align}$$

              \n\n
              rwp = EpiAware.RandomWalk(Normal(),\n    truncated(Normal(0.0, 0.02), 0.0, Inf))
              \n
              RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(Distributions.Normal{Float64}(μ=0.0, σ=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.02); lower=0.0, upper=Inf))
              \n\n\n

              Direct infection EpiModel

              This is a simple model where the unobserved log-infections are directly generated by the latent process $Z$.

              $$\\log I_t = \\log I_0 + Z_t.$$

              As discussed above, we still ask for a defined generation interval, which can be used to calculate $\\mathcal{R}_t$.

              \n\n
              truth_GI = Gamma(2, 5)
              \n
              Distributions.Gamma{Float64}(α=2.0, θ=5.0)
              \n\n\n

              The EpiData constructor performs double interval censoring to convert our continuous estimate of the generation interval into a discretized version. We also implement right truncation using the keyword D_gen.

              \n\n
              model_data = EpiData(truth_GI, D_gen = 10.0)
              \n
              EpiData{Float64, typeof(exp)}([0.056722565637507084, 0.0944813112842689, 0.11668723746973605, 0.12773814268679387, 0.130948497124998, 0.1287976856448402, 0.12312384970427602, 0.11527489380975135, 0.10622581663782862], 9, exp)
              \n\n\n

              We can supply a prior for the initial log_infections.

              \n\n
              log_I0_prior = Normal(log(100.0), 1.0)
              \n
              Distributions.Normal{Float64}(μ=4.605170185988092, σ=1.0)
              \n\n\n

              And construct the EpiModel.

              \n\n
              epi_model = DirectInfections(model_data, log_I0_prior)
              \n
              DirectInfections{Normal{Float64}}(EpiData{Float64, typeof(exp)}([0.056722565637507084, 0.0944813112842689, 0.11668723746973605, 0.12773814268679387, 0.130948497124998, 0.1287976856448402, 0.12312384970427602, 0.11527489380975135, 0.10622581663782862], 9, exp), Distributions.Normal{Float64}(μ=4.605170185988092, σ=1.0))
              \n\n\n

              Delayed Observations ObservationModel

              The observation model is a negative binomial distribution parameterised with mean $\\mu$ and 'successes' parameter $r$. The standard deviation relative to the mean $\\sigma_{\\text{rel}} = \\sigma / \\mu$ for negative binomial observations is,

              $$\\sigma_{\\text{rel}} =(1/\\sqrt{\\mu}) + (1 / \\sqrt{r}).$$

              It is standard to use a half-t distribution for standard deviation priors (e.g. as argued in this paper); we specialise this to a Half-Normal prior and use an a priori assumption that a typical observation fluctuation around the mean (when the mean is $\\sim\\mathcal{O}(10^2)$) would be 10%. This implies a standard deviation prior,

              $$1 / \\sqrt{r} \\sim \\text{HalfNormal}\\Big(0.1 ~\\sqrt{{\\pi \\over 2}}\\Big).$$

              The $\\sqrt{{\\pi \\over 2}}$ factor ensures the correct prior mean (see here).

              The expected observed cases are delayed infections. Delays are implemented as the action of a sparse kernel on the infections $I(t)$.

              $$y_t \\sim \\text{NegBinomial}\\Big(\\mu = \\sum_{s\\geq 0} K[t, t-s] I(s), r\\Big). \\\\$$

              \n\n\n

              We also set up the inference to occur over 100 days.

              \n\n
              time_horizon = 100
              \n
              100
              \n\n\n

              We choose a simple observation model where infections are observed 0, 1, 2, 3 days later with equal probability.

              \n\n
              obs_model = EpiAware.DelayObservations(\n    fill(0.25, 4),\n    time_horizon,\n    truncated(Normal(0, 0.1 * sqrt(pi) / sqrt(2)), 0.0, Inf)\n)
              \n
              DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(sparse(Integer[1, 2, 3, 4, 2, 3, 4, 5, 3, 4  …  97, 98, 99, 100, 98, 99, 100, 99, 100, 100], Integer[1, 1, 1, 1, 2, 2, 2, 2, 3, 3  …  97, 97, 97, 97, 98, 98, 98, 99, 99, 100], [0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25  …  0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25], 100, 100), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.12533141373155002); lower=0.0, upper=Inf))
              \n\n\n

              Generate cases from the EpiAware model

              Having chosen an EpiModel, LatentModel and ObservationModel, we can implement the model as a Turing model using make_epi_aware.

              By giving missing to the first argument, we indicate that case data will be generated from the model rather than treated as fixed.

              \n\n
              full_epi_aware_mdl = make_epi_aware(missing, time_horizon;\n    epi_model = epi_model,\n    latent_model = rwp,\n    observation_model = obs_model)
              \n
              Model{typeof(make_epi_aware), (:y_t, :time_steps), (:epi_model, :latent_model, :observation_model, :pos_shift), (:y_t,), Tuple{Missing, Int64}, Tuple{DirectInfections{Normal{Float64}}, RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, Float64}, DefaultContext}(EpiAware.make_epi_aware, (y_t = missing, time_steps = 100), (epi_model = DirectInfections{Normal{Float64}}(EpiData{Float64, typeof(exp)}([0.056722565637507084, 0.0944813112842689, 0.11668723746973605, 0.12773814268679387, 0.130948497124998, 0.1287976856448402, 0.12312384970427602, 0.11527489380975135, 0.10622581663782862], 9, exp), Distributions.Normal{Float64}(μ=4.605170185988092, σ=1.0)), latent_model = RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(Distributions.Normal{Float64}(μ=0.0, σ=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.02); lower=0.0, upper=Inf)), observation_model = DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(sparse(Integer[1, 2, 3, 4, 2, 3, 4, 5, 3, 4  …  97, 98, 99, 100, 98, 99, 100, 99, 100, 100], Integer[1, 1, 1, 1, 2, 2, 2, 2, 3, 3  …  97, 97, 97, 97, 98, 98, 98, 99, 99, 100], [0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25  …  0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25], 100, 100), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.12533141373155002); lower=0.0, upper=Inf)), pos_shift = 1.0e-6), DefaultContext())
              \n\n\n

              We choose some fixed parameters:

              • Initial incidence is 100.

              • In the direct infection model, the initial incidence and in the initial value of the random walk form a non-identifiable pair. Therefore, we fix $Z_0 = 0$.

              \n\n
              fixed_parameters = (rw_init = 0.0, init_incidence = log(100.0))
              \n
              (rw_init = 0.0, init_incidence = 4.605170185988092)
              \n\n\n

              We fix these parameters using fix, and generate a random epidemic.

              \n\n
              cond_generative_model = fix(full_epi_aware_mdl, fixed_parameters)
              \n
              Model{typeof(make_epi_aware), (:y_t, :time_steps), (:epi_model, :latent_model, :observation_model, :pos_shift), (:y_t,), Tuple{Missing, Int64}, Tuple{DirectInfections{Normal{Float64}}, RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, Float64}, DynamicPPL.FixedContext{@NamedTuple{rw_init::Float64, init_incidence::Float64}, DefaultContext}}(EpiAware.make_epi_aware, (y_t = missing, time_steps = 100), (epi_model = DirectInfections{Normal{Float64}}(EpiData{Float64, typeof(exp)}([0.056722565637507084, 0.0944813112842689, 0.11668723746973605, 0.12773814268679387, 0.130948497124998, 0.1287976856448402, 0.12312384970427602, 0.11527489380975135, 0.10622581663782862], 9, exp), Distributions.Normal{Float64}(μ=4.605170185988092, σ=1.0)), latent_model = RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(Distributions.Normal{Float64}(μ=0.0, σ=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.02); lower=0.0, upper=Inf)), observation_model = DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(sparse(Integer[1, 2, 3, 4, 2, 3, 4, 5, 3, 4  …  97, 98, 99, 100, 98, 99, 100, 99, 100, 100], Integer[1, 1, 1, 1, 2, 2, 2, 2, 3, 3  …  97, 97, 97, 97, 98, 98, 98, 99, 99, 100], [0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25  …  0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25], 100, 100), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.12533141373155002); lower=0.0, upper=Inf)), pos_shift = 1.0e-6), FixedContext((rw_init = 0.0, init_incidence = 4.605170185988092), DynamicPPL.DefaultContext()))
              \n\n
              random_epidemic = rand(cond_generative_model)
              \n
              (ϵ_t = [-0.20112065205042204, 0.3243818279907073, -0.9802665412256135, 1.5899093334448624, 1.1765040928990662, -1.097368025541141, 0.21404530796180515, 0.036040772035425335, 0.19674066706454033, -1.1183113744954276  …  2.925721551631621, -0.2811334030236948, -1.0839466539012816, -0.01945211239319567, -0.525340621257463, -1.5377346782206847, -1.1521273618730166, 0.9841152135889327, 0.03839535544638948, -0.9408477622299963], σ²_RW = 0.007678264825513905, neg_bin_cluster_factor = 0.07657768923618172, var\"y_t[1]\" = 28, var\"y_t[2]\" = 59, var\"y_t[3]\" = 73, var\"y_t[4]\" = 99, var\"y_t[5]\" = 94, var\"y_t[6]\" = 118, var\"y_t[7]\" = 93, var\"y_t[8]\" = 90, var\"y_t[9]\" = 95, var\"y_t[10]\" = 96, var\"y_t[11]\" = 94, var\"y_t[12]\" = 110, var\"y_t[13]\" = 87, var\"y_t[14]\" = 86, var\"y_t[15]\" = 88, var\"y_t[16]\" = 80, var\"y_t[17]\" = 88, var\"y_t[18]\" = 74, var\"y_t[19]\" = 87, var\"y_t[20]\" = 77, var\"y_t[21]\" = 100, var\"y_t[22]\" = 103, var\"y_t[23]\" = 87, var\"y_t[24]\" = 167, var\"y_t[25]\" = 117, var\"y_t[26]\" = 126, var\"y_t[27]\" = 152, var\"y_t[28]\" = 153, var\"y_t[29]\" = 179, var\"y_t[30]\" = 184, var\"y_t[31]\" = 212, var\"y_t[32]\" = 160, var\"y_t[33]\" = 155, var\"y_t[34]\" = 151, var\"y_t[35]\" = 170, var\"y_t[36]\" = 172, var\"y_t[37]\" = 192, var\"y_t[38]\" = 205, var\"y_t[39]\" = 185, var\"y_t[40]\" = 168, var\"y_t[41]\" = 169, var\"y_t[42]\" = 149, var\"y_t[43]\" = 161, var\"y_t[44]\" = 168, var\"y_t[45]\" = 184, var\"y_t[46]\" = 187, var\"y_t[47]\" = 215, var\"y_t[48]\" = 255, var\"y_t[49]\" = 254, var\"y_t[50]\" = 209, var\"y_t[51]\" = 221, var\"y_t[52]\" = 190, var\"y_t[53]\" = 256, var\"y_t[54]\" = 257, var\"y_t[55]\" = 240, var\"y_t[56]\" = 245, var\"y_t[57]\" = 222, var\"y_t[58]\" = 274, var\"y_t[59]\" = 221, var\"y_t[60]\" = 229, var\"y_t[61]\" = 264, var\"y_t[62]\" = 324, var\"y_t[63]\" = 323, var\"y_t[64]\" = 280, var\"y_t[65]\" = 281, var\"y_t[66]\" = 235, var\"y_t[67]\" = 222, var\"y_t[68]\" = 262, var\"y_t[69]\" = 248, var\"y_t[70]\" = 271, var\"y_t[71]\" = 232, var\"y_t[72]\" = 243, var\"y_t[73]\" = 262, var\"y_t[74]\" = 253, var\"y_t[75]\" = 297, var\"y_t[76]\" = 250, var\"y_t[77]\" = 304, var\"y_t[78]\" = 283, var\"y_t[79]\" = 289, var\"y_t[80]\" = 303, var\"y_t[81]\" = 345, var\"y_t[82]\" = 301, var\"y_t[83]\" = 230, var\"y_t[84]\" = 226, var\"y_t[85]\" = 264, var\"y_t[86]\" = 243, var\"y_t[87]\" = 265, var\"y_t[88]\" = 259, var\"y_t[89]\" = 301, var\"y_t[90]\" = 306, var\"y_t[91]\" = 253, var\"y_t[92]\" = 257, var\"y_t[93]\" = 337, var\"y_t[94]\" = 292, var\"y_t[95]\" = 246, var\"y_t[96]\" = 275, var\"y_t[97]\" = 256, var\"y_t[98]\" = 259, var\"y_t[99]\" = 212, var\"y_t[100]\" = 226)
              \n\n
              true_infections = generated_quantities(cond_generative_model, random_epidemic).I_t
              \n
              100-element Vector{Float64}:\n  98.2531043136177\n 101.08593876245777\n  92.76547032111293\n 106.63282760687878\n 118.21244895560318\n 107.37482705017699\n 109.40774047821078\n   ⋮\n 271.51631386943313\n 237.2886030084489\n 214.50244950617125\n 233.82078619031455\n 234.6087824059075\n 216.04288275076027
              \n\n
              generated_obs = generated_quantities(cond_generative_model, random_epidemic).generated_y_t
              \n
              100-element Vector{Int64}:\n  28\n  59\n  73\n  99\n  94\n 118\n  93\n   ⋮\n 246\n 275\n 256\n 259\n 212\n 226
              \n\n\n\n\n","category":"page"},{"location":"examples/getting_started/#Inference","page":"Getting started","title":"Inference","text":"","category":"section"},{"location":"examples/getting_started/","page":"Getting started","title":"Getting started","text":"
              \n

              Fixing $Z_0 = 0$ for the random walk was based on inference principles; in this model $Z_0$ and $\\log I_0$ are non-identifiable.

              However, we now treat the generated data as truth_data and make inference without fixing any other parameters.

              We do the inference by MCMC/NUTS using the Turing NUTS sampler with default warm-up steps.

              \n\n\n

              The observation model supports partially complete data. To test this we set some of the generated observations to be missing.

              \n\n
              truth_data = generated_obs
              \n
              100-element Vector{Int64}:\n  28\n  59\n  73\n  99\n  94\n 118\n  93\n   ⋮\n 246\n 275\n 256\n 259\n 212\n 226
              \n\n
              let\n    truth_data = Union{Int, Missing}[truth_data...]\n    truth_data[vcat([3, 5], 10:20)] .= missing\nend
              \n
              13-element view(::Vector{Union{Missing, Int64}}, [3, 5, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]) with eltype Union{Missing, Int64}:\n missing\n missing\n missing\n missing\n missing\n missing\n missing\n missing\n missing\n missing\n missing\n missing\n missing
              \n\n\n

              We now make the model but fixing the initial condition of the random walk to be 0.

              \n\n
              \ninference_mdl = fix(\n    make_epi_aware(truth_data, time_horizon; epi_model = epi_model,\n        latent_model = rwp, observation_model = obs_model),\n    (rw_init = 0.0,)\n)
              \n
              Model{typeof(make_epi_aware), (:y_t, :time_steps), (:epi_model, :latent_model, :observation_model, :pos_shift), (), Tuple{Vector{Int64}, Int64}, Tuple{DirectInfections{Normal{Float64}}, RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}, Float64}, DynamicPPL.FixedContext{@NamedTuple{rw_init::Float64}, DefaultContext}}(EpiAware.make_epi_aware, (y_t = [28, 59, 73, 99, 94, 118, 93, 90, 95, 96  …  253, 257, 337, 292, 246, 275, 256, 259, 212, 226], time_steps = 100), (epi_model = DirectInfections{Normal{Float64}}(EpiData{Float64, typeof(exp)}([0.056722565637507084, 0.0944813112842689, 0.11668723746973605, 0.12773814268679387, 0.130948497124998, 0.1287976856448402, 0.12312384970427602, 0.11527489380975135, 0.10622581663782862], 9, exp), Distributions.Normal{Float64}(μ=4.605170185988092, σ=1.0)), latent_model = RandomWalk{Normal{Float64}, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(Distributions.Normal{Float64}(μ=0.0, σ=1.0), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.02); lower=0.0, upper=Inf)), observation_model = DelayObservations{Float64, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}(sparse(Integer[1, 2, 3, 4, 2, 3, 4, 5, 3, 4  …  97, 98, 99, 100, 98, 99, 100, 99, 100, 100], Integer[1, 1, 1, 1, 2, 2, 2, 2, 3, 3  …  97, 97, 97, 97, 98, 98, 98, 99, 99, 100], [0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25  …  0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25], 100, 100), Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.12533141373155002); lower=0.0, upper=Inf)), pos_shift = 1.0e-6), FixedContext((rw_init = 0.0,), DynamicPPL.DefaultContext()))
              \n\n\n

              Initialising inference

              It is possible for the default warm-up process for NUTS to get stuck in low probability or otherwise degenerate regions of parameter space.

              To make NUTS more robust we provide manypathfinder, which is built on pathfinder variational inference from Pathfinder.jl. manypathfinder runs nruns pathfinder processes on the inference problem and returns the pathfinder run with maximum estimated ELBO.

              manypathfinder differs from Pathfinder.multipathfinder; multipathfinder is aimed at sampling from a potentially non-Gaussian target distribution which is first approximated as a uniformly weighted collection of normal approximations from pathfinder runs. manypathfinder is aimed at moving rapidly to a 'good' part of parameter space, and is robust to runs that fail.

              \n\n
              best_pf = manypathfinder(inference_mdl, 10; nruns = 20, executor = Transducers.ThreadedEx());
              \n\n\n\n

              We can use draws from the best pathfinder run to initialise NUTS.

              \n\n
              best_pf.draws_transformed
              \n
              iterationchainϵ_t[1]ϵ_t[2]ϵ_t[3]ϵ_t[4]ϵ_t[5]ϵ_t[6]...
              1111.745171.3581-0.0387770.228527-0.4804070.259776
              2211.739621.966440.2950680.6561230.2273570.221215
              3311.404741.62267-0.1437840.913016-0.212709-0.0593607
              4411.989611.71673-0.1428980.558126-0.3282210.278743
              5511.839921.609580.3763271.052640.1482740.358493
              6611.568591.3840.2471130.647590.08802540.372194
              7711.802511.857780.1909470.6891620.1628090.01183
              8812.126171.599590.1070831.17275-0.2168740.267504
              9911.517171.88996-0.2935890.3183850.146870.335505
              101011.622831.57171-0.2107330.3872-0.5130320.182996
              \n\n
              init_params = collect.(eachrow(best_pf.draws_transformed.value[1:4, :, 1]))
              \n
              4-element Vector{Vector{Float64}}:\n [1.7451655099039896, 1.3581032532363655, -0.03877701856907928, 0.2285266700718785, -0.4804068463802501, 0.2597760274093325, -0.72329152773823, -0.004892028183825725, 0.30371317937314624, -0.32360713324738793  …  -0.5689987466999031, -1.069569400690149, -0.006438779107851761, 0.12884728013030794, -0.4672454241686336, 0.46663119811384574, -0.5012766699018397, 0.0008443002919903706, 4.503333034127224, 0.10534955177982713]\n [1.7396239623829572, 1.9664442888461848, 0.29506847368478406, 0.6561230931785763, 0.2273567486516539, 0.22121457862999117, -0.4457969493025312, 0.4933162930390505, 0.5387079340846253, 0.27978737948037563  …  -0.26300000639899557, -1.0805768157412714, -0.2542903280341239, 0.3450712794593742, -0.41739949283166994, 0.5189459878777011, -0.2407079378968638, 0.0013653332228509899, 4.40762882420073, 0.12133329079454364]\n [1.404744492096187, 1.6226685394934721, -0.14378388615214202, 0.9130160620170369, -0.21270896000629097, -0.059360707272928415, -0.38729412142681957, 0.3423067106403553, 0.4423491622359259, 0.03856933348504296  …  -0.6845978102747797, -1.0692267124647084, -0.1627366413969122, 0.4973267181861608, -0.2835918876396665, 0.7286810493507772, -0.14793772810228983, 0.0012724611829520575, 4.364726695263788, 0.11203225705238223]\n [1.9896074833972643, 1.7167261947293577, -0.14289800447970843, 0.5581263080306341, -0.3282212207720035, 0.27874272220508667, -0.9184680486663799, 0.5964059365665253, 0.5400965831532827, 0.5955039812253591  …  -0.6508965991738329, -1.2498544646684877, -0.35170979199800345, 0.3039492244589093, -0.5363086467209464, 0.37978414791959775, 0.2275674079496028, 0.0015740315635789245, 4.358921604618429, 0.09758372594224174]
              \n\n\n

              NB: We are running this inference run for speed rather than accuracy as a demonstration. Use a higher target acceptance and more samples in a typical workflow.

              \n\n
              chn = sample(inference_mdl,\n    NUTS(; adtype = AutoReverseDiff(true)),\n    MCMCThreads(),\n    250,\n    4;\n    init_params,\n    drop_warmup = true)
              \n
              iterationchainϵ_t[1]ϵ_t[2]ϵ_t[3]ϵ_t[4]ϵ_t[5]ϵ_t[6]...
              112610.871897-0.130904-0.919508-0.0402790.4041450.142873
              21271-1.10953-1.335251.81422-1.51668-0.190370.143744
              312810.6702260.648365-0.791261-0.151059-0.2185021.85856
              412910.417116-0.272337-0.671253-0.6251520.595989-1.02231
              513012.452411.239250.745333-1.320331.2241-1.87632
              613111.973160.2079611.01515-2.481431.1235-1.28877
              71321-0.590797-0.778022-0.1580630.748415-0.172206-0.517546
              81331-1.402431.63989-0.339024-0.8965040.2595890.507557
              91341-0.915431.201580.554195-0.549210.4315130.349227
              1013510.510275-1.12904-0.3463790.104358-0.576478-0.813278
              ...
              \n\n\n

              Predictive plotting

              We can spaghetti plot generated case data from the version of the model which hasn't conditioned on case data using posterior parameters inferred from the version conditioned on observed data. This is known as posterior predictive checking, and is a useful diagnostic tool for Bayesian inference (see here).

              Because we are using synthetic data we can also plot the model predictions for the unobserved infections and check that (at least in this example) we were able to capture some unobserved/latent variables in the process accurate.

              \n\n
              let\n    post_check_mdl = fix(full_epi_aware_mdl, (rw_init = 0.0,))\n    post_check_y_t = mapreduce(hcat, generated_quantities(post_check_mdl, chn)) do gen\n        gen.generated_y_t\n    end\n\n    predicted_I_t = mapreduce(hcat, generated_quantities(inference_mdl, chn)) do gen\n        gen.I_t\n    end\n\n    p1 = plot(post_check_y_t, c = :grey, alpha = 0.05, lab = \"\")\n    scatter!(p1, truth_data,\n        lab = \"Observed cases\",\n        xlabel = \"Time\",\n        ylabel = \"Cases\",\n        title = \"Post. predictive checking: cases\",\n        ylims = (-0.5, maximum(truth_data) * 1.5),\n        c = :green)\n\n    p2 = plot(predicted_I_t, c = :grey, alpha = 0.05, lab = \"\")\n    scatter!(p2, true_infections,\n        lab = \"Actual infections\",\n        xlabel = \"Time\",\n        ylabel = \"Unobserved Infections\",\n        title = \"Post. predictions: infections\",\n        ylims = (-0.5, maximum(true_infections) * 1.5),\n        c = :red)\n\n    plot(p1, p2,\n        layout = (1, 2),\n        size = (700, 400))\nend
              \n\n\n\n

              As well as checking the posterior predictions for latent infections, we can also check how well inference recovered unknown parameters, such as the random walk variance or the cluster factor of the negative binomial observations.

              \n\n
              let\n    parameters_to_plot = (:σ²_RW, :neg_bin_cluster_factor)\n\n    plts = map(parameters_to_plot) do name\n        var_samples = chn[name] |> vec\n        histogram(var_samples,\n            bins = 50,\n            norm = :pdf,\n            lw = 0,\n            fillalpha = 0.5,\n            lab = \"MCMC\")\n        vline!([getfield(random_epidemic, name)], lab = \"True value\")\n        title!(string(name))\n    end\n    plot(plts..., layout = (2, 1))\nend
              \n\n\n","category":"page"},{"location":"examples/getting_started/#Reproductive-number-back-calculation","page":"Getting started","title":"Reproductive number back-calculation","text":"","category":"section"},{"location":"examples/getting_started/","page":"Getting started","title":"Getting started","text":"
              \n

              As mentioned at the top, we don't directly use the concept of reproductive numbers in this note. However, we can back-calculate the implied $\\mathcal{R}(t)$ values, conditional on the specified generation interval being correct.

              Here we spaghetti plot posterior sampled time-varying reproductive numbers against the actual.

              \n\n
              let\n    n = epi_model.data.len_gen_int\n    Rt_denom = [dot(reverse(epi_model.data.gen_int), true_infections[(t - n):(t - 1)])\n                for t in (n + 1):length(true_infections)]\n    true_Rt = true_infections[(n + 1):end] ./ Rt_denom\n\n    predicted_Rt = mapreduce(hcat, generated_quantities(inference_mdl, chn)) do gen\n        _It = gen.I_t\n        _Rt_denom = [dot(reverse(epi_model.data.gen_int), _It[(t - n):(t - 1)])\n                     for t in (n + 1):length(_It)]\n        Rt = _It[(n + 1):end] ./ _Rt_denom\n    end\n\n    plt = plot((n + 1):time_horizon, predicted_Rt, c = :grey, alpha = 0.05, lab = \"\")\n    plot!(plt, (n + 1):time_horizon, true_Rt,\n        lab = \"true Rt\",\n        xlabel = \"Time\",\n        ylabel = \"Rt\",\n        title = \"Post. predictions: reproductive number\",\n        c = :red,\n        lw = 2)\nend
              \n\n\n","category":"page"},{"location":"examples/getting_started/","page":"Getting started","title":"Getting started","text":"EditURL = \"https://github.com/CDCgov/Rt-without-renewal/blob/main/docs/src/examples/getting_started.jl\"","category":"page"},{"location":"lib/public/#Public-Documentation","page":"Public API","title":"Public Documentation","text":"","category":"section"},{"location":"lib/public/","page":"Public API","title":"Public API","text":"Documentation for EpiAware.jl's public interface.","category":"page"},{"location":"lib/public/","page":"Public API","title":"Public API","text":"See the Internals section of the manual for internal package docs covering all submodules.","category":"page"},{"location":"lib/public/#Contents","page":"Public API","title":"Contents","text":"","category":"section"},{"location":"lib/public/","page":"Public API","title":"Public API","text":"Pages = [\"public.md\"]\nDepth = 2:2","category":"page"},{"location":"lib/public/#Index","page":"Public API","title":"Index","text":"","category":"section"},{"location":"lib/public/","page":"Public API","title":"Public API","text":"Pages = [\"public.md\"]","category":"page"},{"location":"lib/public/#Public-API","page":"Public API","title":"Public API","text":"","category":"section"},{"location":"#EpiAware.jl","page":"EpiAware.jl: Real-time epidemic monitoring","title":"EpiAware.jl","text":"","category":"section"},{"location":"","page":"EpiAware.jl: Real-time epidemic monitoring","title":"EpiAware.jl: Real-time epidemic monitoring","text":"Infectious disease situational awareness modelling toolkit for Julia.","category":"page"},{"location":"","page":"EpiAware.jl: Real-time epidemic monitoring","title":"EpiAware.jl: Real-time epidemic monitoring","text":"A package for building and fitting situational awareness models for infectious diseases. The package is designed to be flexible and extensible, and to provide a consistent interface for fitting and simulating models.","category":"page"},{"location":"#Package-Features","page":"EpiAware.jl: Real-time epidemic monitoring","title":"Package Features","text":"","category":"section"},{"location":"","page":"EpiAware.jl: Real-time epidemic monitoring","title":"EpiAware.jl: Real-time epidemic monitoring","text":"Flexible: The package is designed to be flexible and extensible, and to provide a consistent interface for fitting and simulating models.\nModular: The package is designed to be modular, with a clear separation between the model and the data.\nExtensible: The package is designed to be extensible, with a clear separation between the model and the data.\nConsistent: The package is designed to provide a consistent interface for fitting and simulating models.\nEfficient: The package is designed to be efficient, with a clear separation between the model and the data.","category":"page"},{"location":"","page":"EpiAware.jl: Real-time epidemic monitoring","title":"EpiAware.jl: Real-time epidemic monitoring","text":"See the Index for the complete list of documented functions and types.","category":"page"},{"location":"#Manual-Outline","page":"EpiAware.jl: Real-time epidemic monitoring","title":"Manual Outline","text":"","category":"section"},{"location":"","page":"EpiAware.jl: Real-time epidemic monitoring","title":"EpiAware.jl: Real-time epidemic monitoring","text":"Pages = [\n \"man/contributing.md\",\n]\nDepth = 1","category":"page"},{"location":"#Library-Outline","page":"EpiAware.jl: Real-time epidemic monitoring","title":"Library Outline","text":"","category":"section"},{"location":"","page":"EpiAware.jl: Real-time epidemic monitoring","title":"EpiAware.jl: Real-time epidemic monitoring","text":"Pages = [\"lib/public.md\", \"lib/internals.md\"]","category":"page"},{"location":"#main-index","page":"EpiAware.jl: Real-time epidemic monitoring","title":"Index","text":"","category":"section"},{"location":"","page":"EpiAware.jl: Real-time epidemic monitoring","title":"EpiAware.jl: Real-time epidemic monitoring","text":"Pages = [\"lib/public.md\"]","category":"page"}] }