From 7f168cf5ed8a2983177347cac3fe557b108aafb6 Mon Sep 17 00:00:00 2001 From: Samuel Brand Date: Fri, 8 Mar 2024 21:31:22 +0000 Subject: [PATCH 1/6] Add the _make_halfnormal_prior internal function with unit tests and included in getting started --- EpiAware/docs/src/examples/getting_started.jl | 18 +++++++++------ EpiAware/src/latentmodels/randomwalk.jl | 3 +-- EpiAware/src/utils/dist.jl | 16 +++++++++++-- EpiAware/test/test_latent-models.jl | 2 +- EpiAware/test/test_observation-models.jl | 2 +- EpiAware/test/test_utilities.jl | 23 +++++++++++++++++++ 6 files changed, 51 insertions(+), 13 deletions(-) diff --git a/EpiAware/docs/src/examples/getting_started.jl b/EpiAware/docs/src/examples/getting_started.jl index 555145e1e..61b4caa1a 100644 --- a/EpiAware/docs/src/examples/getting_started.jl +++ b/EpiAware/docs/src/examples/getting_started.jl @@ -89,14 +89,14 @@ In `EpiAware` we provide a constructor for random walk latent models with priors ```math \begin{align} Z_0 &\sim \mathcal{N}(0,1),\\ -\sigma^2_Z &\sim \text{HalfNormal}(0.01). +\sigma_{RW} &\sim \text{HalfNormal}(0.1 * \sqrt{\pi} / \sqrt{2})). \end{align} ``` " # ╔═╡ 56ae496b-0094-460b-89cb-526627991717 rwp = EpiAware.RandomWalk(Normal(), - truncated(Normal(0.0, 0.02), 0.0, Inf)) + EpiAware._make_halfnormal_prior(0.1)) # ╔═╡ 767beffd-1ef5-4e6c-9ac6-edb52e60fb44 md" @@ -179,7 +179,7 @@ We choose a simple observation model where infections are observed 0, 1, 2, 3 da obs_model = EpiAware.DelayObservations( fill(0.25, 4), time_horizon, - truncated(Normal(0, 0.1 * sqrt(pi) / sqrt(2)), 0.0, Inf) + EpiAware._make_halfnormal_prior(0.1) ) # ╔═╡ e49713e8-4840-4083-8e3f-fc52d791be7b @@ -298,12 +298,15 @@ 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.** +**NB: We are running this inference run for speed rather than accuracy as a demonstration. You might want to use a higher target acceptance and more samples in a typical workflow.** " +# ╔═╡ 946b1c43-e750-40c9-9f14-79da9735e437 +target_acc_prob = 0.8 + # ╔═╡ 3eb5ec5e-aae7-478e-84fb-80f2e9f85b4c chn = sample(inference_mdl, - NUTS(; adtype = AutoReverseDiff(true)), + NUTS(target_acc_prob; adtype = AutoReverseDiff(true)), MCMCThreads(), 250, 4; @@ -360,7 +363,7 @@ As well as checking the posterior predictions for latent infections, we can also # ╔═╡ 10d8fe24-83a6-47ac-97b7-a374481473d3 let - parameters_to_plot = (:σ²_RW, :neg_bin_cluster_factor) + parameters_to_plot = (:σ_RW, :neg_bin_cluster_factor) plts = map(parameters_to_plot) do name var_samples = chn[name] |> vec @@ -449,7 +452,8 @@ end # ╠═073a1d40-456a-450e-969f-11b23eb7fd1f # ╠═0379b058-4c35-440a-bc01-aafa0178bdbf # ╠═a7798f71-9bb5-4506-9476-0cc11553b9e2 -# ╟─4deb3a51-781d-48c4-91f6-6adf2b1affcf +# ╠═4deb3a51-781d-48c4-91f6-6adf2b1affcf +# ╠═946b1c43-e750-40c9-9f14-79da9735e437 # ╠═3eb5ec5e-aae7-478e-84fb-80f2e9f85b4c # ╟─30498cc7-16a5-441a-b8cd-c19b220c60c1 # ╠═e9df22b8-8e4d-4ab7-91ea-c01f2239b3e5 diff --git a/EpiAware/src/latentmodels/randomwalk.jl b/EpiAware/src/latentmodels/randomwalk.jl index d3fe8bd36..20a2ae7f3 100644 --- a/EpiAware/src/latentmodels/randomwalk.jl +++ b/EpiAware/src/latentmodels/randomwalk.jl @@ -10,9 +10,8 @@ end @model function generate_latent(latent_model::RandomWalk, n) ϵ_t ~ MvNormal(ones(n)) - σ²_RW ~ latent_model.var_prior + σ_RW ~ latent_model.var_prior rw_init ~ latent_model.init_prior - σ_RW = sqrt(σ²_RW) rw = Vector{eltype(ϵ_t)}(undef, n) rw[1] = rw_init + σ_RW * ϵ_t[1] diff --git a/EpiAware/src/utils/dist.jl b/EpiAware/src/utils/dist.jl index 3bae62382..e759c95a9 100644 --- a/EpiAware/src/utils/dist.jl +++ b/EpiAware/src/utils/dist.jl @@ -150,8 +150,6 @@ function r_to_R(r, w::AbstractVector) end """ - NegativeBinomialMeanClust(μ, α) - Compute the mean-cluster factor negative binomial distribution. # Arguments @@ -167,3 +165,17 @@ function NegativeBinomialMeanClust(μ, α) r = μ^2 / ex_σ² return NegativeBinomial(r, p) end + +""" +Create a half-normal prior distribution with the specified mean. + +# Arguments +- `prior_mean::AbstractFloat`: The mean of the prior distribution. + +# Returns +- `Truncated{Normal}`: The half-normal prior distribution. + +""" +function _make_halfnormal_prior(prior_mean::AbstractFloat) + return truncated(Normal(0.0, prior_mean * sqrt(pi) / sqrt(2)), 0.0, Inf)::Distribution +end diff --git a/EpiAware/test/test_latent-models.jl b/EpiAware/test/test_latent-models.jl index f1357722c..9bab27770 100644 --- a/EpiAware/test/test_latent-models.jl +++ b/EpiAware/test/test_latent-models.jl @@ -8,7 +8,7 @@ rw_process = EpiAware.RandomWalk(Normal(0.0, 1.0), truncated(Normal(0.0, 0.05), 0.0, Inf)) model = EpiAware.generate_latent(rw_process, n) - fixed_model = fix(model, (σ²_RW = 1.0, init_rw_value = 0.0)) #Fixing the standard deviation of the random walk process + fixed_model = fix(model, (σ_RW = 1.0, init_rw_value = 0.0)) #Fixing the standard deviation of the random walk process n_samples = 1000 samples_day_5 = sample(fixed_model, Prior(), n_samples) |> chn -> mapreduce(vcat, generated_quantities(fixed_model, chn)) do gen diff --git a/EpiAware/test/test_observation-models.jl b/EpiAware/test/test_observation-models.jl index df6b4f7c9..5437176c4 100644 --- a/EpiAware/test/test_observation-models.jl +++ b/EpiAware/test/test_observation-models.jl @@ -1,6 +1,6 @@ @testitem "Testing delay obs against theoretical properties" begin using DynamicPPL, Turing, Distributions - using HypothesisTests#: ExactOneSampleKSTest, pvalue + using HypothesisTests # Set up test data with fixed infection I_t = [10.0, 20.0, 30.0] diff --git a/EpiAware/test/test_utilities.jl b/EpiAware/test/test_utilities.jl index d7cf01c7b..b39be39dc 100644 --- a/EpiAware/test/test_utilities.jl +++ b/EpiAware/test/test_utilities.jl @@ -196,3 +196,26 @@ end @test df == expected_df end end + +@testitem "Testing _make_halfnormal_prior" begin + using Distributions, HypothesisTests + @testset "Check distribution type" begin + prior_mean = 10.0 + prior_dist = EpiAware._make_halfnormal_prior(prior_mean) + @test typeof(prior_dist) <: Distribution + end + + @testset "Check distribution properties" begin + prior_mean = 2.0 + prior_dist = EpiAware._make_halfnormal_prior(prior_mean) + #Check Distributions.jl mean function + @test mean(prior_dist) ≈ prior_mean + samples = rand(prior_dist, 10_000) + #Check mean from direct sampling of folded distribution and ANOVA and Variance F test comparisons + direct_samples = randn(10_000) * prior_mean * sqrt(pi) / sqrt(2) .|> abs + mean_pval = OneWayANOVATest(samples, direct_samples) |> pvalue + @test mean_pval > 1e-6 #Very unlikely to fail if the model is correctly implemented + var_pval = VarianceFTest(samples, direct_samples) |> pvalue + @test var_pval > 1e-6 #Very unlikely to fail if the model is correctly implemented + end +end From 4229320317fd68b4f2df9342af989c8426b20ab1 Mon Sep 17 00:00:00 2001 From: Samuel Brand Date: Fri, 8 Mar 2024 21:31:43 +0000 Subject: [PATCH 2/6] tweak manypathfinder to default do more samples for elbo estimate --- EpiAware/src/inferencemethods/manypathfinder.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/EpiAware/src/inferencemethods/manypathfinder.jl b/EpiAware/src/inferencemethods/manypathfinder.jl index 7fe546f63..1bbe42516 100644 --- a/EpiAware/src/inferencemethods/manypathfinder.jl +++ b/EpiAware/src/inferencemethods/manypathfinder.jl @@ -94,8 +94,8 @@ largest ELBO estimate. - `best_pfs::PathfinderResult`: Best pathfinder result by estimated ELBO. """ function manypathfinder(mdl::DynamicPPL.Model, ndraws; nruns = 4, - maxiters = 50, max_tries = 100, kwargs...) + maxiters = 50, max_tries = 100, ndraws_elbo = 10, kwargs...) _run_manypathfinder(mdl; nruns, ndraws, maxiters, kwargs...) |> - pfs -> _continue_manypathfinder!(pfs, mdl; max_tries, nruns, kwargs...) |> + pfs -> _continue_manypathfinder!(pfs, mdl; max_tries, nruns, ndraws_elbo, kwargs...) |> pfs -> _get_best_elbo_pathfinder(pfs) end From fb774e534289279ae254b805b3e5f2c208ffddd9 Mon Sep 17 00:00:00 2001 From: Samuel Brand Date: Mon, 11 Mar 2024 15:19:14 +0000 Subject: [PATCH 3/6] Revert "tweak manypathfinder to default do more samples for elbo estimate" This reverts commit 4229320317fd68b4f2df9342af989c8426b20ab1. --- EpiAware/src/inferencemethods/manypathfinder.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/EpiAware/src/inferencemethods/manypathfinder.jl b/EpiAware/src/inferencemethods/manypathfinder.jl index 1bbe42516..7fe546f63 100644 --- a/EpiAware/src/inferencemethods/manypathfinder.jl +++ b/EpiAware/src/inferencemethods/manypathfinder.jl @@ -94,8 +94,8 @@ largest ELBO estimate. - `best_pfs::PathfinderResult`: Best pathfinder result by estimated ELBO. """ function manypathfinder(mdl::DynamicPPL.Model, ndraws; nruns = 4, - maxiters = 50, max_tries = 100, ndraws_elbo = 10, kwargs...) + maxiters = 50, max_tries = 100, kwargs...) _run_manypathfinder(mdl; nruns, ndraws, maxiters, kwargs...) |> - pfs -> _continue_manypathfinder!(pfs, mdl; max_tries, nruns, ndraws_elbo, kwargs...) |> + pfs -> _continue_manypathfinder!(pfs, mdl; max_tries, nruns, kwargs...) |> pfs -> _get_best_elbo_pathfinder(pfs) end From 1f46019e16e3156951f91d6339a894a32a14a27e Mon Sep 17 00:00:00 2001 From: Samuel Brand Date: Mon, 11 Mar 2024 16:18:51 +0000 Subject: [PATCH 4/6] Rename `var_prior` to `std_prior` --- EpiAware/src/latentmodels/randomwalk.jl | 4 ++-- EpiAware/test/test_latent-models.jl | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/EpiAware/src/latentmodels/randomwalk.jl b/EpiAware/src/latentmodels/randomwalk.jl index 20a2ae7f3..1dda224b2 100644 --- a/EpiAware/src/latentmodels/randomwalk.jl +++ b/EpiAware/src/latentmodels/randomwalk.jl @@ -1,6 +1,6 @@ struct RandomWalk{D <: Sampleable, S <: Sampleable} <: AbstractLatentModel init_prior::D - var_prior::S + std_prior::S end function default_rw_priors() @@ -10,7 +10,7 @@ end @model function generate_latent(latent_model::RandomWalk, n) ϵ_t ~ MvNormal(ones(n)) - σ_RW ~ latent_model.var_prior + σ_RW ~ latent_model.std_prior rw_init ~ latent_model.init_prior rw = Vector{eltype(ϵ_t)}(undef, n) diff --git a/EpiAware/test/test_latent-models.jl b/EpiAware/test/test_latent-models.jl index 9bab27770..0c16430b1 100644 --- a/EpiAware/test/test_latent-models.jl +++ b/EpiAware/test/test_latent-models.jl @@ -33,8 +33,8 @@ end end @testset "Testing RandomWalk constructor" begin init_prior = Normal(0.0, 1.0) - var_prior = truncated(Normal(0.0, 0.05), 0.0, Inf) - rw_process = RandomWalk(init_prior, var_prior) + std_prior = truncated(Normal(0.0, 0.05), 0.0, Inf) + rw_process = RandomWalk(init_prior, std_prior) @test rw_process.init_prior == init_prior - @test rw_process.var_prior == var_prior + @test rw_process.std_prior == std_prior end From b844043f260e49b8378302810f12f410a16e77ee Mon Sep 17 00:00:00 2001 From: Samuel Brand Date: Mon, 11 Mar 2024 16:19:08 +0000 Subject: [PATCH 5/6] removed mistaken type tag --- EpiAware/src/utils/dist.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/EpiAware/src/utils/dist.jl b/EpiAware/src/utils/dist.jl index e759c95a9..8a2baac39 100644 --- a/EpiAware/src/utils/dist.jl +++ b/EpiAware/src/utils/dist.jl @@ -177,5 +177,5 @@ Create a half-normal prior distribution with the specified mean. """ function _make_halfnormal_prior(prior_mean::AbstractFloat) - return truncated(Normal(0.0, prior_mean * sqrt(pi) / sqrt(2)), 0.0, Inf)::Distribution + return truncated(Normal(0.0, prior_mean * sqrt(pi) / sqrt(2)), 0.0, Inf) end From e36b550a0f7a82642e6c02666f2eac6178895e42 Mon Sep 17 00:00:00 2001 From: Samuel Brand Date: Mon, 11 Mar 2024 16:19:14 +0000 Subject: [PATCH 6/6] Update getting_started.jl --- EpiAware/docs/src/examples/getting_started.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/EpiAware/docs/src/examples/getting_started.jl b/EpiAware/docs/src/examples/getting_started.jl index 61b4caa1a..5140c14f7 100644 --- a/EpiAware/docs/src/examples/getting_started.jl +++ b/EpiAware/docs/src/examples/getting_started.jl @@ -449,10 +449,10 @@ end # ╠═b4033728-b321-4100-8194-1fd9fe2d268d # ╟─9222b436-9445-4039-abbf-25c8cddb7f63 # ╠═197a4fbb-b71a-475a-bb78-28ff613e3094 -# ╠═073a1d40-456a-450e-969f-11b23eb7fd1f +# ╟─073a1d40-456a-450e-969f-11b23eb7fd1f # ╠═0379b058-4c35-440a-bc01-aafa0178bdbf # ╠═a7798f71-9bb5-4506-9476-0cc11553b9e2 -# ╠═4deb3a51-781d-48c4-91f6-6adf2b1affcf +# ╟─4deb3a51-781d-48c4-91f6-6adf2b1affcf # ╠═946b1c43-e750-40c9-9f14-79da9735e437 # ╠═3eb5ec5e-aae7-478e-84fb-80f2e9f85b4c # ╟─30498cc7-16a5-441a-b8cd-c19b220c60c1