From 141dbcce4187d92753364ee6584c56d0558482d9 Mon Sep 17 00:00:00 2001 From: Rich FitzJohn Date: Thu, 26 Sep 2024 16:29:57 +0100 Subject: [PATCH 1/6] Fix doc rendering issues --- R/sample-manual.R | 4 ++++ _pkgdown.yml | 12 ++++++++++++ man/monty_sample_manual_collect.Rd | 4 +--- man/monty_sample_manual_prepare_continue.Rd | 5 +---- 4 files changed, 18 insertions(+), 7 deletions(-) diff --git a/R/sample-manual.R b/R/sample-manual.R index 70dd0b83..e0a97895 100644 --- a/R/sample-manual.R +++ b/R/sample-manual.R @@ -159,6 +159,8 @@ monty_sample_manual_run <- function(chain_id, path, progress = NULL) { ##' [monty_sample_manual_prepare] and [monty_sample_manual_run]. If ##' any chain has not completed, we will error. ##' +##' @title Collect manually run samples +##' ##' @inheritParams monty_sample_manual_run ##' ##' @param samples Samples from the parent run. You need to provide @@ -273,6 +275,8 @@ sample_manual_info_chain <- function(complete) { ##' [monty_sample_manual_prepare] is to [monty_sample]. The original ##' set of samples do not need to have been run manually. ##' +##' @title Prepare to continue sampling with manual scheduling +##' ##' @param save_samples Control over saving samples into the inputs. ##' The choices here are `hash` (the default) where we save a hash ##' and validate that in [monty_sample_manual_collect], `value` diff --git a/_pkgdown.yml b/_pkgdown.yml index d42df802..f80156ba 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -77,3 +77,15 @@ reference: - monty_rng_distributed_pointer - monty_rng_distributed_state - monty_rng_pointer + +articles: + - title: Basics + navbar: Basics + contents: + - dsl + - samples + - title: Details + navbar: Details + contents: + - dsl-errors + - samplers diff --git a/man/monty_sample_manual_collect.Rd b/man/monty_sample_manual_collect.Rd index e626f28c..72a868f4 100644 --- a/man/monty_sample_manual_collect.Rd +++ b/man/monty_sample_manual_collect.Rd @@ -2,9 +2,7 @@ % Please edit documentation in R/sample-manual.R \name{monty_sample_manual_collect} \alias{monty_sample_manual_collect} -\title{Collect samples from chains that have been run with -\link{monty_sample_manual_prepare} and \link{monty_sample_manual_run}. If -any chain has not completed, we will error.} +\title{Collect manually run samples} \usage{ monty_sample_manual_collect(path, samples = NULL, restartable = FALSE) } diff --git a/man/monty_sample_manual_prepare_continue.Rd b/man/monty_sample_manual_prepare_continue.Rd index b54ad262..56fc3d01 100644 --- a/man/monty_sample_manual_prepare_continue.Rd +++ b/man/monty_sample_manual_prepare_continue.Rd @@ -2,10 +2,7 @@ % Please edit documentation in R/sample-manual.R \name{monty_sample_manual_prepare_continue} \alias{monty_sample_manual_prepare_continue} -\title{Prepare to continue sampling from a model, with manual chain -orchestration. This function is to \link{monty_sample_continue} what -\link{monty_sample_manual_prepare} is to \link{monty_sample}. The original -set of samples do not need to have been run manually.} +\title{Prepare to continue sampling with manual scheduling} \usage{ monty_sample_manual_prepare_continue( samples, From 0d4efccdf7d1af3ca143c53f09759d6fb238f2ca Mon Sep 17 00:00:00 2001 From: Rich FitzJohn Date: Thu, 26 Sep 2024 16:31:46 +0100 Subject: [PATCH 2/6] Bump version --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 5ed6714c..9146d764 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: monty Title: Monte Carlo Models -Version: 0.2.6 +Version: 0.2.7 Authors@R: c(person("Rich", "FitzJohn", role = c("aut", "cre"), email = "rich.fitzjohn@gmail.com"), person("Wes", "Hinsley", role = "aut"), From 4d80cfc260240185a3077d8fd68c11d2b4671185 Mon Sep 17 00:00:00 2001 From: Rich FitzJohn Date: Mon, 30 Sep 2024 18:58:35 +0100 Subject: [PATCH 3/6] Doc improvements --- vignettes/monty.Rmd | 2 +- vignettes/samplers.Rmd | 72 +++++++++++++++++------------------------- 2 files changed, 30 insertions(+), 44 deletions(-) diff --git a/vignettes/monty.Rmd b/vignettes/monty.Rmd index 65dde092..29ad479b 100644 --- a/vignettes/monty.Rmd +++ b/vignettes/monty.Rmd @@ -68,7 +68,7 @@ likelihood <- monty_model_function(fn, fixed = list(data = data)) likelihood ``` -The prior we'll make much nicer to work with in the future, but here we construct the density by hand as a sum of normally distributed priors on `a` and `b`, and a weak uniform prior on `sigma`. +We construct a prior for this model using the monty DSL (`vignette("dsl")`), using normally distributed priors on `a` and `b`, and a weak uniform prior on `sigma`. ```{r} prior <- monty_dsl({ diff --git a/vignettes/samplers.Rmd b/vignettes/samplers.Rmd index 6638cc21..88c0e369 100644 --- a/vignettes/samplers.Rmd +++ b/vignettes/samplers.Rmd @@ -25,45 +25,13 @@ library(monty) ## The bendy banana -This example shows HMC outperforming a random walk on a two dimensional banana-shaped function. Our model takes two parameters `alpha` and `beta`, and is based on two successive simple draws, with the one conditional on the other: +This example shows HMC outperforming a random walk on a two dimensional banana-shaped function. Our model takes two parameters `alpha` and `beta`, and is based on two successive simple draws, with the one conditional on the other, so $\beta \sim Normal(1,0)$ and $\alpha \sim Normal(\beta^2, \sigma)$, with $\sigma$ the standard deviation of the conditional draw. -\[ -\beta \sim Normal(1,0) \\ -\alpha \sim Normal(\beta^2, \sigma) -\] - -with $\sigma$ the standard deviation of the conditional draw. - -We'll use R's `dnorm` for the likelihood calculations and then differentiate the density by hand (using the formula of the gaussian density) to derive the gradient by hand. We skip the details here as this will be automated in an upcoming version of the package. +We include this example withn the package; here we create a model with $\sigma = 0.5$ ```{r} -banana_model <- function(sigma = 0.5) { - monty_model(list( - parameters = c("alpha", "beta"), - direct_sample = function(rng) { - beta <- rng$random_normal(1) - alpha <- rng$normal(1, beta^2, sigma) - c(alpha, beta) - }, - density = function(x) { - alpha <- x[[1]] - beta <- x[[2]] - dnorm(beta, log = TRUE) + dnorm((alpha - beta^2) / sigma, log = TRUE) - }, - gradient = function(x) { - alpha <- x[[1]] - beta <- x[[2]] - c((beta^2 - alpha) / sigma^2, - -beta + 2 * beta * (alpha - beta^2) / sigma^2) - }, - domain = rbind(c(-Inf, Inf), c(-Inf, Inf)))) -} -``` - -Let's create a model with $\sigma = 0.5$ - -```{r create_banana_model} -m <- banana_model(0.5) +m <- monty_example("banana", sigma = 0.5) +m ``` We can plot a greyscale visualisation of its density by computing the density over a grid. Normally this is not possible of course: @@ -71,7 +39,9 @@ We can plot a greyscale visualisation of its density by computing the density ov ```{r} a <- seq(-2, 6, length.out = 1000) b <- seq(-2, 2, length.out = 1000) -z <- outer(a, b, function(a, b) dnorm(b) * dnorm((a - b^2) / 0.5)) +z <- outer(a, b, function(alpha, beta) { + exp(monty_model_density(m, rbind(alpha, beta))) +}) image(a, b, z, xlab = "alpha", ylab = "beta") ``` @@ -104,24 +74,40 @@ points(s[1, ], s[2, ], pch = 19, col = "#00000055") It is not generally possible to directly sample from a density (otherwise MCMC and similar methods would not exist!). In these cases we need to use a sampler based on the density and if available possibly the gradient of the density. -We can start with a basic random-walk sampler +We can start with a basic random-walk sampler: ```{r RW_sampling} sampler_rw <- monty_sampler_random_walk(vcv = diag(2) * 0.01) -res_rw <- monty_sample(m, sampler_rw, 2000) -plot(res_rw$pars, pch = 19, col = "#ff222277") +res_rw <- monty_sample(m, sampler_rw, 1000) +plot(t(drop(res_rw$pars)), type = "l", col = "#0000ff66", + xlim = range(a), ylim = range(b)) lines(z95[, 1], z95[, 2]) ``` -As we can see this is not great. We can probably improve the samples here by finding a better variance covariance matrix (VCV), but a single VCV will not hold well over the whole surface because it is not very similar to a multivariate normal (that is, the appropriate VCV will change depending on the position in parameter space) +As we can see this is not great, exhibiting strong random walk behaviour as it slowly explore the surface (this is over 1,000 steps). Another way to view this is the parameters varying over steps: + +``` +matplot(t(drop(res_rw$pars)), lty = 1, type = "l", col = c(2, 4), + xlab = "Step", ylab = "Value") +``` + +We can probably improve the samples here by finding a better variance covariance matrix (VCV), but a single VCV will not hold well over the whole surface because it is not very similar to a multivariate normal (that is, the appropriate VCV will change depending on the position in parameter space) Let's try the Hamiltonian Monte Carlo (HMC) sampler, which uses the gradient to move efficiently in parameter space: ```{r HMC_sampling} sampler_hmc <- monty_sampler_hmc(epsilon = 0.1, n_integration_steps = 10) -res_hmc <- monty_sample(m, sampler_hmc, 2000) -plot(res_hmc$pars, pch = 19, col = "#ff222277") +res_hmc <- monty_sample(m, sampler_hmc, 1000) +plot(t(drop(res_hmc$pars)), type = "l", col = "#0000ff33", + xlim = range(a), ylim = range(b)) lines(z95[, 1], z95[, 2]) ``` +or viewed over steps: + +```{r} +matplot(t(drop(res_hmc$pars)), lty = 1, type = "l", col = c(2, 4), + xlab = "Step", ylab = "Value") +``` + Clearly better! From 92deddd867c0e26b6f84b9f7057656ee9a2f5371 Mon Sep 17 00:00:00 2001 From: Rich FitzJohn Date: Mon, 30 Sep 2024 18:59:11 +0100 Subject: [PATCH 4/6] Fix spelling --- vignettes/samplers.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/samplers.Rmd b/vignettes/samplers.Rmd index 88c0e369..c136eed5 100644 --- a/vignettes/samplers.Rmd +++ b/vignettes/samplers.Rmd @@ -27,7 +27,7 @@ library(monty) This example shows HMC outperforming a random walk on a two dimensional banana-shaped function. Our model takes two parameters `alpha` and `beta`, and is based on two successive simple draws, with the one conditional on the other, so $\beta \sim Normal(1,0)$ and $\alpha \sim Normal(\beta^2, \sigma)$, with $\sigma$ the standard deviation of the conditional draw. -We include this example withn the package; here we create a model with $\sigma = 0.5$ +We include this example within the package; here we create a model with $\sigma = 0.5$ ```{r} m <- monty_example("banana", sigma = 0.5) From 3f3dcb428f5d57e4f6580d049d926050898907bd Mon Sep 17 00:00:00 2001 From: Rich FitzJohn Date: Tue, 1 Oct 2024 22:27:14 +0100 Subject: [PATCH 5/6] Update vignettes/samplers.Rmd Co-authored-by: Wes Hinsley --- vignettes/samplers.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/samplers.Rmd b/vignettes/samplers.Rmd index c136eed5..632433f0 100644 --- a/vignettes/samplers.Rmd +++ b/vignettes/samplers.Rmd @@ -84,7 +84,7 @@ plot(t(drop(res_rw$pars)), type = "l", col = "#0000ff66", lines(z95[, 1], z95[, 2]) ``` -As we can see this is not great, exhibiting strong random walk behaviour as it slowly explore the surface (this is over 1,000 steps). Another way to view this is the parameters varying over steps: +As we can see this is not great, exhibiting strong random walk behaviour as it slowly explores the surface (this is over 1,000 steps). Another way to view this is the parameters varying over steps: ``` matplot(t(drop(res_rw$pars)), lty = 1, type = "l", col = c(2, 4), From 355200a75c26e625735d5a8a304ec0a85061f524 Mon Sep 17 00:00:00 2001 From: Rich FitzJohn Date: Wed, 2 Oct 2024 07:34:19 +0100 Subject: [PATCH 6/6] Additional fixes --- vignettes/samplers.Rmd | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/vignettes/samplers.Rmd b/vignettes/samplers.Rmd index 632433f0..eecf6e8d 100644 --- a/vignettes/samplers.Rmd +++ b/vignettes/samplers.Rmd @@ -34,7 +34,7 @@ m <- monty_example("banana", sigma = 0.5) m ``` -We can plot a greyscale visualisation of its density by computing the density over a grid. Normally this is not possible of course: +We can plot a visualisation of its density by computing the density over a grid. Normally this is not possible of course: ```{r} a <- seq(-2, 6, length.out = 1000) @@ -54,7 +54,7 @@ image(a, b, z, xlab = "alpha", ylab = "beta") points(s[1, ], s[2, ], pch = 19, col = "#00000055") ``` -It is also possible to compute the 95% confidence interval of the distribution using the relationship between the standard bivariate normal distribution and the banana shaped distribution as defined above. We can check than roughly 10 samples (out of 200) are out of this 95% CI contour. +It is also possible to compute the 95% confidence interval of the distribution using the relationship between the standard bivariate normal distribution and the banana shaped distribution as defined above. We can check that roughly 10 samples (out of 200) are out of this 95% CI contour. ```{r} theta <- seq(0, 2 * pi, length.out = 10000)