Add some examples for documentation

.github/workflows/document-EpiAware.yaml

@@ -1,7 +1,7 @@
 name: Documenter
 on:
   push:
-    branches: main
+    branches: [main, master]
     tags: [v*]
   pull_request:
 

EpiAware/docs/Project.toml

@@ -1,4 +1,16 @@
 [deps]
 Changelog = "5217a498-cd5d-4ec6-b8c2-9b85a09b6e3e"
+DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-EpiAware = "b2eeebe4-5992-4301-9193-7ebc9f62c855"
+DynamicPPL = "366bfd00-2699-11ea-058f-f148b4cae6d8"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
+Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"

EpiAware/test/predictive_checking/discretized_pmfs.jl → EpiAware/docs/examples/discretized_pmfs.jl

@@ -81,7 +81,6 @@ plt1 = let
         xlabel = "Days",
         ylabel = "Probability")
 end
-savefig(plt1, joinpath(@__DIR__(), "assets/", "discrete_pmf_daily.png"))
 
 # For hourly censoring the difference is not noticable.
 
@@ -100,4 +99,3 @@ plt2 = let
         xlabel = "Days",
         ylabel = "Probability")
 end
-savefig(plt2, joinpath(@__DIR__(), "assets/", "discrete_pmf_hourly.png"))

EpiAware/test/predictive_checking/toy_model_log_infs_RW.jl → EpiAware/docs/examples/toy_model_log_infs_RW.jl

@@ -68,7 +68,6 @@ using Random
 using DynamicPPL
 using Statistics
 using DataFramesMeta
-using CSV # For outputting the MCMC chain
 
 Random.seed!(0)
 
@@ -142,7 +141,8 @@ truth_data = random_epidemic.y_t
 model = make_epi_aware(truth_data, time_horizon, ; epi_model = epi_model,
     latent_model_model = rwp, observation_model = obs_model,
     pos_shift = 1e-6)
-@time chn = sample(model,
+
+chn = sample(model,
     NUTS(; adtype = AutoReverseDiff(true)),
     MCMCThreads(),
     250,
@@ -189,5 +189,4 @@ We can use `spread_draws` to convert the MCMC chain into a tidybayes format.
 =#
 
 df_chn = spread_draws(chn)
-save_path = joinpath(@__DIR__, "assets/toy_model_log_infs_RW_draws.csv")
-CSV.write(save_path, df_chn)
+first(df_chn, 10)

EpiAware/docs/make.jl

@@ -1,12 +1,16 @@
-using Documenter
+using Documenter, Changelog
 using EpiAware
 
 include("changelog.jl")
 include("pages.jl")
 
 makedocs(; sitename = "EpiAware.jl",
     authors = "Samuel Brand, Zachary Susswein, Sam Abbott, and contributors",
-    clean = true, doctest = true, linkcheck = true,
+    clean = true, doctest = true, #linkcheck = true,
     warnonly = [:docs_block, :missing_docs],
     modules = [EpiAware],
-    pages = pages)
+    pages = pages,
+    format = Documenter.HTML(
+        prettyurls = get(ENV, "CI", nothing) == "true"
+    )
+)

EpiAware/docs/pages.jl

@@ -1,3 +1,8 @@
 pages = [
-    "EpiAware.jl: Real-time epidemic monitoring" => "index.md"
+    "EpiAware.jl: Real-time epidemic monitoring" => "index.md",
+    "Examples" => [
+        "Double censoring" => "man/discretized_pmfs.md",
+        "Calculating exponential growth rate from reproduction number" => "man/R_to_r.md",
+        "Running inference for a toy model" => "man/toy_model_log_infs_RW.md"
+    ]
 ]

EpiAware/docs/src/man/R_to_r.md

@@ -0,0 +1,86 @@
+```@meta
+EditURL = "../../examples/R_to_r.jl"
+```
+
+# Fast approximation for `r` from `R₀`
+I use the negative moment generating function (MGF).
+
+Let
+```math
+G(r) = \sum_{i=1}^{\infty} w_i e^{-r i}.
+```
+and
+```math
+f(r, \mathcal{R}_0) = \mathcal{R}_0 G(r) - 1.
+```
+then the connection between `R₀` and `r` is given by
+```math
+f(r, \mathcal{R}_0) = 0.
+```
+Given an estimate of $\mathcal{R}_0$ we implicit solve for $r$ using a root
+finder algorithm. In this note, I test a fast approximation for $r$ which
+should have good autodifferentiation properties. The idea is to start from the
+small $r$ approximation to the solution of $f(r, \mathcal{R}_0) = 0$ and then
+apply one step of Newton's method. The small $r$ approximation is given by
+```math
+r_0 = { \mathcal{R}_0 - 1 \over  \mathcal{R}_0 \langle W \rangle }.
+```
+where $\langle W \rangle$ is the mean of the generation interval.
+
+To rapidly improve the estimate for `r` we use Newton steps given by
+```math
+r_{n+1} = r_n - {\mathcal{R}_0 G(r) - 1\over \mathcal{R}_0 G'(r)}.
+```
+
+### Test mode
+
+Run this script in Test environment mode. If not, run the following command:
+
+```julia
+using TestEnv
+TestEnv.activate()
+```
+
+````@example R_to_r
+using EpiAware
+using Distributions
+using StatsPlots
+````
+
+Create a discrete probability mass function (PMF) for a negative binomial distribution
+with left truncation at 1.
+
+````@example R_to_r
+w = create_discrete_pmf(NegativeBinomial(2, 0.5), D = 20.0) |>
+    p -> p[2:end] ./ sum(p[2:end])
+
+#
+
+jitter = 1e-17
+idxs = 0:10
+doubling_times = [1.0, 3.5, 7.0, 14.0]
+
+errors = mapreduce(hcat, doubling_times) do T_2
+    true_r = log(2) / T_2 # 7 day doubling time
+    R0 = EpiAware.r_to_R(true_r, w)
+
+    return map(idxs) do ns
+        @time r = EpiAware.R_to_r(R0, w, newton_steps = ns)
+        abs(r - true_r) + jitter
+    end
+end
+
+plot(idxs,
+    errors,
+    yscale = :log10,
+    xlabel = "Newton steps",
+    ylabel = "abs. Error",
+    title = "Fast approximation for r",
+    lab = ["T_2 = 1.0" "T_2 = 3.5" "T_2 = 7.0" "T_2 = 14.0"],
+    yticks = [0.0, 1e-15, 1e-10, 1e-5, 1e0] |> x -> (x .+ jitter, string.(x)),
+    xticks = 0:2:10)
+````
+
+---
+
+*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*

EpiAware/docs/src/man/discretized_pmfs.md

@@ -0,0 +1,118 @@
+```@meta
+EditURL = "../../examples/discretized_pmfs.jl"
+```
+
+# Discretized PMFs
+
+## Analytical PMF for the Exponential distribution
+
+For unit testing it is useful to have an analytically solvable example of double interval censoring. Easiest distribution we
+could solve analytically but was also not completely trivial was $X \sim \text{Exp}(1)$ day delay with daily interval censoring.
+W.l.o.g. Primary censored obs time is $t = 0$, and we go from the double-censored interval with uniform on interval primary
+approximation as per [here](https://www.medrxiv.org/content/10.1101/2024.01.12.24301247v1).
+
+For above example the secondary censored obs time can be $s = 0, 1, 2,...$ days. The probability mass function is:
+
+```math
+P_S(s) = \int_0^1 \int_s^{s+1} f_X(y-x)dy dx.
+```
+
+This splits into two cases: $s = 0$ and $s \geq 1$.
+
+**Case 1:** $s=0$
+
+```math
+P_S(0) = \int_0^1 \int_x^1 \exp(-(y - x)) dy dx = \exp(-1).
+```
+
+_NB: the density is zero for negative values._
+
+**Case 2:** $s \geq 1$
+
+```math
+P_S(s) = \int_0^1 \int_s^{s+1} \exp(-(y - x)) dy dx = (1 - \exp(-1)) (\exp(1) - 1) \exp(-s).
+```
+
+we can directly check that the above is a discrete prob distribution. First, non-negativity is obvious. Second,
+normalisation to 1 can be directly calculated,
+
+```math
+\begin{align}
+ \sum_{s \geq 1} P_S(s)&= (1 - \exp(-1)) (\exp(1) - 1) \sum_{s \geq 1}  \exp(-s)  \\
+&=  (1 - \exp(-1)) (\exp(1) - 1)  {\exp(-1) \over 1 - \exp(-1)} \\
+& = 1 -  \exp(-1).
+\end{align}
+```
+Therefore,
+
+```math
+ \sum_{s \geq 0} P_S(s) = 1.
+```
+
+## Predictive checking for the `create_discrete_pmf` function
+
+This predictive checking shows the difference between the two methods of the `create_discrete_pmf` function
+for creating a discrete PMF from a continuous distribution with a given discretization interval `Δd` and upper bound `D`.
+
+The default method is double censoring based on censoring intervals of width `Δd`. The basic method is based on the
+same but with the assumption that the primary event happens at the edge of the censoring interval. The left edge implies that
+the discrete PMF starts at `0`, the right edge implies that the discrete PMF starts at `Δd`.
+
+````@example discretized_pmfs
+using EpiAware
+using StatsPlots
+using Distributions
+````
+
+Example distribution is a Gamma distribution with shape 2 and scale 3/2 (mean = 3 days, std = √4.5 days) with an upper bound of 21 days.
+
+````@example discretized_pmfs
+cont_dist = Gamma(2, 3.0 / 2)
+D = 21.0
+````
+
+For daily censoring there is a fairly big difference between the two methods, as well as left/right interval endpointing.
+
+````@example discretized_pmfs
+plt1 = let
+    Δd = 1
+    ts = (0.0:Δd:(D - Δd)) |> collect
+    pmf1 = create_discrete_pmf(cont_dist, Val(:single_censored); Δd = Δd, D = D)
+    pmf2 = create_discrete_pmf(cont_dist; Δd = Δd, D = D)
+
+    bar(ts,
+        [pmf1;; pmf2],
+        fillalpha = 0.5,
+        lw = 0,
+        title = "Discrete PMF with Δd = 1 day",
+        label = ["Single censoring (midpoint primary)" "Double Censoring"],
+        xlabel = "Days",
+        ylabel = "Probability")
+end
+# savefig(plt1, joinpath(@__DIR__(), "assets/", "discrete_pmf_daily.png"))
+````
+
+For hourly censoring the difference is not noticable.
+
+````@example discretized_pmfs
+plt2 = let
+    Δd = 1 / 24
+    ts = (0.0:Δd:(D - Δd)) |> collect
+    pmf1 = create_discrete_pmf(cont_dist, Val(:single_censored); Δd = Δd, D = D)
+    pmf2 = create_discrete_pmf(cont_dist; Δd = Δd, D = D)
+
+    bar(ts,
+        [pmf1;; pmf2],
+        fillalpha = 0.5,
+        lw = 0,
+        title = "Discrete PMF with Δd = 1 hour",
+        label = ["Single censoring (midpoint primary)" "Double Censoring"],
+        xlabel = "Days",
+        ylabel = "Probability")
+end
+# savefig(plt2, joinpath(@__DIR__(), "assets/", "discrete_pmf_hourly.png"))
+````
+
+---
+
+*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*