-
-
+
@@ -106,118 +108,182 @@
- +-
Using Likelihood Blocks
-Bob Verity and Pete Winskill
+Bob Verity and +Pete Winskill
-2021-12-16
+2024-06-26
- Source:vignettes/blocks.Rmd+ Source:vignettes/blocks.Rmdblocks.RmdOnce we are familiar with the basics, there are some more advanced drjacoby techniques that can be useful when running computationally intensive likelihoods, or when fitting multi-level models. These rely on the fact that drjacoby works by breaking the likelihood down into a series of independent blocks, which are combined to produce the overall likelihood.
-This vignette demonstrates the use of likelihood blocks to fit a simple random-effects model.
----Problem motivation
-For this example we will work with the
-chickwtsdataset (one of the datasets that comes with R), which lists the weights of various chicks broken down by the feed they were given.--
head(chickwts, 12) -#> weight feed -#> 1 179 horsebean -#> 2 160 horsebean -#> 3 136 horsebean -#> 4 227 horsebean -#> 5 217 horsebean -#> 6 168 horsebean -#> 7 108 horsebean -#> 8 124 horsebean -#> 9 143 horsebean -#> 10 140 horsebean -#> 11 309 linseed -#> 12 229 linseedWe will pass this into drjacoby as a list over groups:
+
+#> Registered S3 method overwritten by 'GGally': +#> method from +#> +.gg ggplot2Once we are familiar with the +basics, there are some more advanced drjacoby techniques +that can be useful when running computationally intensive likelihoods, +or when fitting multi-level models. These rely on the fact that +drjacoby works by breaking the likelihood down into a series of +independent blocks, which are combined to produce the overall +likelihood.
+This vignette demonstrates the use of likelihood blocks to fit a +simple random-effects model.
+++Problem motivation +
+For this example we will work with the
chickwtsdataset +(one of the datasets that comes with R), which lists the weights of +various chicks broken down by the feed they were given.--
data_list <- split(chickwts$weight, f = chickwts$feed) -data_list -#> $casein -#> [1] 368 390 379 260 404 318 352 359 216 222 283 332 -#> -#> $horsebean -#> [1] 179 160 136 227 217 168 108 124 143 140 -#> -#> $linseed -#> [1] 309 229 181 141 260 203 148 169 213 257 244 271 -#> -#> $meatmeal -#> [1] 325 257 303 315 380 153 263 242 206 344 258 -#> -#> $soybean -#> [1] 243 230 248 327 329 250 193 271 316 267 199 171 158 248 -#> -#> $sunflower -#> [1] 423 340 392 339 341 226 320 295 334 322 297 318There are 6 feed groups in total, and we will assume that values within a group are normally distributed each with a distinct unknown mean, leading to 6 parameters \(\mu_1\) to \(\mu_6\). For simplicity we will assume the same standard deviation \(\sigma\) in all groups. To make this a random effects framework, we will assume that \(\mu_1\) to \(\mu_6\) are themselves random draws from a normal distribution with mean \(\phi\) and standard deviation \(\tau\). Finally, we will put diffuse priors on \(\sigma\), \(\phi\) and \(\tau\) so that the data can “speak for itself”. For the more mathematically minded, the complete model can be written as follows:
+head(chickwts, 12) +#> weight feed +#> 1 179 horsebean +#> 2 160 horsebean +#> 3 136 horsebean +#> 4 227 horsebean +#> 5 217 horsebean +#> 6 168 horsebean +#> 7 108 horsebean +#> 8 124 horsebean +#> 9 143 horsebean +#> 10 140 horsebean +#> 11 309 linseed +#> 12 229 linseedWe will pass this into drjacoby as a list over groups:
+++
data_list <- split(chickwts$weight, f = chickwts$feed) +data_list +#> $casein +#> [1] 368 390 379 260 404 318 352 359 216 222 283 332 +#> +#> $horsebean +#> [1] 179 160 136 227 217 168 108 124 143 140 +#> +#> $linseed +#> [1] 309 229 181 141 260 203 148 169 213 257 244 271 +#> +#> $meatmeal +#> [1] 325 257 303 315 380 153 263 242 206 344 258 +#> +#> $soybean +#> [1] 243 230 248 327 329 250 193 271 316 267 199 171 158 248 +#> +#> $sunflower +#> [1] 423 340 392 339 341 226 320 295 334 322 297 318There are 6 feed groups in total, and we will assume that values +within a group are normally distributed each with a distinct unknown +mean, leading to 6 parameters \(\mu_1\) +to \(\mu_6\). For simplicity we will +assume the same standard deviation \(\sigma\) in all groups. To make this a +random effects framework, we will assume that \(\mu_1\) to \(\mu_6\) are themselves random draws from a +normal distribution with mean \(\phi\) +and standard deviation \(\tau\). +Finally, we will put diffuse priors on \(\sigma\), \(\phi\) and \(\tau\) so that the data can “speak for +itself”. For the more mathematically minded, the complete model can be +written as follows:
\[ \begin{aligned} - x_{i,j} &\sim Normal(\mu_i, \sigma^2) \hspace{5mm} \textrm{for } i \in 1:6, \\ - \mu_i &\sim Normal(\phi, \tau^2) \hspace{5mm} \textrm{for } i \in 1:6, \\ + x_{i,j} &\sim Normal(\mu_i, \sigma^2) \hspace{5mm} \textrm{for } i +\in 1:6, \\ + \mu_i &\sim Normal(\phi, \tau^2) \hspace{5mm} \textrm{for } i \in +1:6, \\ \sigma &\sim Gamma(0.01, 100), \\ \phi &\sim Normal(0, 1000), \\ \tau &\sim Gamma(0.01, 100) \end{aligned} -\] One way to fit this model in drjacoby is to write a likelihood function that simply loops through all the data and calculates the probability as a function of the parameters. However, this is very wasteful - if we are updating the parameter \(\mu_1\) then we only need to know the observed chick weights in feed group 1, along with the current estimates of \(\sigma\), \(\phi\) and \(\tau\). The same is true for all other \(\mu_i\) values, which only depend on a small subset of the data. In this example the inefficiency probably doesn’t matter too much, because the dataset is small and our likelihood is fast enough that we can brute force the problem, but it’s easy to imagine situations where this wouldn’t be the case. A better way to approach this problem is through likelihood blocks.
+\] One way to fit this model in drjacoby is to write a +likelihood function that simply loops through all the data and +calculates the probability as a function of the parameters. However, +this is very wasteful - if we are updating the parameter \(\mu_1\) then we only need to know the +observed chick weights in feed group 1, along with the current estimates +of \(\sigma\), \(\phi\) and \(\tau\). The same is true for all other +\(\mu_i\) values, which only depend on +a small subset of the data. In this example the inefficiency probably +doesn’t matter too much, because the dataset is small and our likelihood +is fast enough that we can brute force the problem, but it’s easy to +imagine situations where this wouldn’t be the case. A better way to +approach this problem is through likelihood blocks.-@@ -327,11 +406,13 @@-Defining blocks
-You can imagine likelihood blocks as a vector of values, each of which gives the log-likelihood of one small part of the model. The overall log-likelihood is simply the sum over this vector. We are free to define as many blocks as we like, and the exact number we need will depend on the model design. A good way of figuring out how many blocks you need is to sketch out a “conditional dependence” table like the one shown below.
++-Defining blocks +
+You can imagine likelihood blocks as a vector of values, each of +which gives the log-likelihood of one small part of the model. The +overall log-likelihood is simply the sum over this vector. We are free +to define as many blocks as we like, and the exact number we need will +depend on the model design. A good way of figuring out how many blocks +you need is to sketch out a “conditional dependence” table like the one +shown below.
-
In the first column we list all the data and free parameters of our model, and in the second column we list every parameter on which the first column depends (“depends” in this context means we need to know this parameter value in order to write down the probability). Don’t put data or fixed values in the second column, only free parameters, i.e. those we are trying to infer. We need one block for every unique combination of parameters in the second column. For example, \(\mu_1\) and \(\mu_2\) depend on the same combination of parameters (\(\phi\) and \(\tau\)) therefore they belong in the same block. We find that we need seven blocks in this example.
-We define blocks in drjacoby using the optional “block” column of the parameters dataframe. For each parameter you should list all the blocks in which it can be found - for example \(\mu_1\) can be found in blocks 1 and 7 in the table above:
-+-
# define parameters dataframe -df_params <- define_params(name = "mu1", min = 0, max = Inf, init = 1, block = c(1, 7), - name = "mu2", min = 0, max = Inf, init = 1, block = c(2, 7), - name = "mu3", min = 0, max = Inf, init = 1, block = c(3, 7), - name = "mu4", min = 0, max = Inf, init = 1, block = c(4, 7), - name = "mu5", min = 0, max = Inf, init = 1, block = c(5, 7), - name = "mu6", min = 0, max = Inf, init = 1, block = c(6, 7), - name = "sigma", min = 0, max = Inf, init = 1, block = 1:6, - name = "phi", min = -Inf, max = Inf, init = 0, block = 7, - name = "tau", min = 0, max = Inf, init = 1, block = 7) - -df_params -#> name min max init block -#> 1 mu1 0 Inf 1 1, 7 -#> 2 mu2 0 Inf 1 2, 7 -#> 3 mu3 0 Inf 1 3, 7 -#> 4 mu4 0 Inf 1 4, 7 -#> 5 mu5 0 Inf 1 5, 7 -#> 6 mu6 0 Inf 1 6, 7 -#> 7 sigma 0 Inf 1 1, 2, 3, 4, 5, 6 -#> 8 phi -Inf Inf 0 7 -#> 9 tau 0 Inf 1 7In the first column we list all the data and free parameters of our +model, and in the second column we list every parameter on which the +first column depends (“depends” in this context means we need to know +this parameter value in order to write down the probability). Don’t put +data or fixed values in the second column, only free parameters, +i.e. those we are trying to infer. We need one block for every +unique combination of parameters in the second column. For +example, \(\mu_1\) and \(\mu_2\) depend on the same combination of +parameters (\(\phi\) and \(\tau\)) therefore they belong in the same +block. We find that we need seven blocks in this example.
+We define blocks in drjacoby using the optional “block” +column of the parameters dataframe. For each parameter you should list +all the blocks in which it can be found - for example \(\mu_1\) can be found in blocks 1 and 7 in +the table above:
++
# define parameters dataframe +df_params <- define_params(name = "mu1", min = 0, max = Inf, init = 1, block = c(1, 7), + name = "mu2", min = 0, max = Inf, init = 1, block = c(2, 7), + name = "mu3", min = 0, max = Inf, init = 1, block = c(3, 7), + name = "mu4", min = 0, max = Inf, init = 1, block = c(4, 7), + name = "mu5", min = 0, max = Inf, init = 1, block = c(5, 7), + name = "mu6", min = 0, max = Inf, init = 1, block = c(6, 7), + name = "sigma", min = 0, max = Inf, init = 1, block = 1:6, + name = "phi", min = -Inf, max = Inf, init = 0, block = 7, + name = "tau", min = 0, max = Inf, init = 1, block = 7) + +df_params +#> name min max init block +#> 1 mu1 0 Inf 1 1, 7 +#> 2 mu2 0 Inf 1 2, 7 +#> 3 mu3 0 Inf 1 3, 7 +#> 4 mu4 0 Inf 1 4, 7 +#> 5 mu5 0 Inf 1 5, 7 +#> 6 mu6 0 Inf 1 6, 7 +#> 7 sigma 0 Inf 1 1, 2, 3, 4, 5, 6 +#> 8 phi -Inf Inf 0 7 +#> 9 tau 0 Inf 1 7--The Likelihood
-Next we need to access these blocks inside the likelihood function. The current block is stored within the
-miscobject, inside an element called “block”. We can access this value in a C++ function usingmisc["block"], or in an R function usingmisc$block. Once we know the current block we can vary the likelihood accordingly. Our objective is to calculate and return just the log-likelihood of this block.Take a look at the following likelihood function, which takes this approach:
+++The Likelihood +
+Next we need to access these blocks inside the likelihood function. +The current block is stored within the
+miscobject, inside +an element called “block”. We can access this value in a C++ function +usingmisc["block"], or in an R function using +misc$block. Once we know the current block we can vary the +likelihood accordingly. Our objective is to calculate and return just +the log-likelihood of this block.Take a look at the following likelihood function, which takes this +approach:
-
-
- Unpack each of the parameters from the
paramsobject
+ - Unpack each of the parameters from the
params+object - Find out which block we are in using
misc["block"] - Split based on block:
-
-
- For blocks 1:6, extract data corresponding to this block and calculate likelihood -
- For block 7, calculate likelihood over all
muparameters
+ - For blocks 1:6, extract data corresponding to this block and +calculate likelihood +
- For block 7, calculate likelihood over all
mu+parameters
stop("cpp function %i not found", function_name); } -
If we were to run this likelihood for blocks 1 through 7 and sum the values, we would obtain the full model log-likelihood.
-For the prior function we don’t need to worry about blocks, we just need to apply the appropriate distributions:
+If we were to run this likelihood for blocks 1 through 7 and sum the +values, we would obtain the full model log-likelihood.
+For the prior function we don’t need to worry about blocks, we just +need to apply the appropriate distributions:
We can then run the MCMC as normal:
---
# run MCMC -mcmc <- run_mcmc(data = data_list, - df_params = df_params, - loglike = "loglike", - logprior = "logprior", - burnin = 1e3, - samples = 1e4, - chains = 5, - silent = TRUE)Assuming we are happy with MCMC performance (convergence, ESS etc.), we can explore the posterior credible intervals of all parameters:
+-
plot_credible(mcmc)# run MCMC +mcmc <- run_mcmc(data = data_list, + df_params = df_params, + loglike = "loglike", + logprior = "logprior", + burnin = 1e3, + samples = 1e4, + chains = 5, + silent = TRUE)Assuming we are happy with MCMC performance (convergence, ESS etc.), +we can explore the posterior credible intervals of all parameters:
++
plot_credible(mcmc)
-
Here we obtain a different mean estimate for each group, with the global mean (
-phi) being somewhere near the centre. If you compare themuestimates to the raw mean of the data in each group you will find that they are “shrunk” towards the global mean. This is the impact of the random effects model - by sharing information between groups we end up pulling everything towards the centre.Likelihood blocks take some getting used to, but for complex models they can make life much easier and be much more efficient. In reality drjacoby always uses blocks internally, it just places all parameters into block 1 if not specified by the user. There are no checks performed on the blocks that you define, so make sure your blocks are independent and cover the complete likelihood of the model.
+Here we obtain a different mean estimate for each group, with the +global mean (
+phi) being somewhere near the centre. If you +compare themuestimates to the raw mean of the data in +each group you will find that they are “shrunk” towards the global mean. +This is the impact of the random effects model - by sharing information +between groups we end up pulling everything towards the centre.Likelihood blocks take some getting used to, but for complex models +they can make life much easier and be much more efficient. In reality +drjacoby always uses blocks internally, it just places all +parameters into block 1 if not specified by the user. There are no +checks performed on the blocks that you define, so make sure your blocks +are independent and cover the complete likelihood of the model.
@@ -340,5 +421,7 @@
+ + diff --git a/docs/articles/blocks_files/figure-html/unnamed-chunk-9-1.png b/docs/articles/blocks_files/figure-html/unnamed-chunk-9-1.png index e7c629f..8fe0721 100644 Binary files a/docs/articles/blocks_files/figure-html/unnamed-chunk-9-1.png and b/docs/articles/blocks_files/figure-html/unnamed-chunk-9-1.png differ diff --git a/docs/articles/checks_double_well.html b/docs/articles/checks_double_well.html index ba58d59..bafaf6f 100644 --- a/docs/articles/checks_double_well.html +++ b/docs/articles/checks_double_well.html @@ -6,7 +6,7 @@
Double well • drjacoby - + @@ -19,6 +19,8 @@ + +@@ -31,14 +33,14 @@@@ -229,11 +237,13 @@- - + Checks @@ -91,7 +93,7 @@
-
-
+
@@ -106,114 +108,120 @@
- +-+-
Double well
-Bob Verity and Pete Winskill
+Bob Verity and +Pete Winskill
-2021-12-16
+2024-06-26
- Source:vignettes/checks_double_well.Rmd+ Source:vignettes/checks_double_well.Rmdchecks_double_well.RmdPurpose: to compare drjacoby results for a challenging problem involving a multimodal posterior, both with and without temperature rungs.
--+-Model
-We assume a single parameter
+mudrawn from a double well potential distribution, defined by the formula:
+## Registered S3 method overwritten by 'GGally': +## method from +## +.gg ggplot2Purpose: to compare drjacoby results for a +challenging problem involving a multimodal posterior, both with and +without temperature rungs.
++-Model +
+We assume a single parameter
mudrawn from a double well +potential distribution, defined by the formula:\[ \begin{aligned} \mu &\propto exp\left(-\gamma(\mu^2 - 1)^2\right) \end{aligned} -\] where \(\gamma\) is a parameter that defines the strength of the well (higher \(\gamma\) leads to a deeper valley and hence more challenging problem). NB, there is no data in this example, as the likelihood is defined exactly by these parameters.
+\] where \(\gamma\) is a +parameter that defines the strength of the well (higher \(\gamma\) leads to a deeper valley and hence +more challenging problem). NB, there is no data in this example, as the +likelihood is defined exactly by these parameters.Likelihood and prior:
Parameters dataframe:
---
L <- 2 -gamma <- 30 -df_params <- define_params(name = "mu", min = -L, max = L, - name = "gamma", min = gamma, max = gamma)-+-Single temperature rung (no Metropolis coupling)
+-
mcmc <- run_mcmc(data = list(x = -1), - df_params = df_params, - loglike = "loglike", - logprior = "logprior", - burnin = 1e3, - samples = 1e5, - chains = 1, - rungs = 1, - silent = TRUE)L <- 2 +gamma <- 30 +df_params <- define_params(name = "mu", min = -L, max = L, + name = "gamma", min = gamma, max = gamma)+Single temperature rung (no Metropolis coupling) +
--
# trace plot -plot_par(mcmc, show = "mu")
+
mcmc <- run_mcmc(data = list(x = -1), + df_params = df_params, + loglike = "loglike", + logprior = "logprior", + burnin = 1e3, + samples = 1e5, + chains = 1, + rungs = 1, + silent = TRUE)+-
# extract posterior draws -output_sub <- subset(mcmc$output, phase == "sampling") -mu_draws <- output_sub$mu - -# get analytical solution -x <- seq(-L, L, l = 1001) -fx <- exp(-gamma*(x^2 - 1)^2) -fx <- fx / sum(fx) * 1/(x[2]-x[1]) - -# overlay plots -hist(mu_draws, breaks = seq(-L, L, l = 201), probability = TRUE, main = "", col = "black") -lines(x, fx, col = 2, lwd = 2)# trace plot +plot_trace(mcmc, show = "mu")
++
# extract posterior draws +output_sub <- subset(mcmc$output, phase == "sampling") +mu_draws <- output_sub$mu + +# get analytical solution +x <- seq(-L, L, l = 1001) +fx <- exp(-gamma*(x^2 - 1)^2) +fx <- fx / sum(fx) * 1/(x[2]-x[1]) + +# overlay plots +hist(mu_draws, breaks = seq(-L, L, l = 201), probability = TRUE, main = "", col = "black") +lines(x, fx, col = 2, lwd = 2)
-+-Multiple temperature rungs
---
mcmc <- run_mcmc(data = list(x = -1), - df_params = df_params, - loglike = "loglike", - logprior = "logprior", - burnin = 1e3, - samples = 1e5, - chains = 1, - rungs = 11, - alpha = 2, - pb_markdown = TRUE)
-## MCMC chain 1 -## burn-in -## - | - |======================================================================| 100% -## acceptance rate: 21.7% -## sampling phase -## - | - |======================================================================| 100% -## acceptance rate: 22% -## chain completed in 4.397968 seconds
-## total MCMC run-time: 4.4 seconds--
# trace plot -plot_par(mcmc, show = "mu")
+
++Multiple temperature rungs +
+++
mcmc <- run_mcmc(data = list(x = -1), + df_params = df_params, + loglike = "loglike", + logprior = "logprior", + burnin = 1e3, + samples = 1e5, + chains = 1, + rungs = 11, + alpha = 2, + pb_markdown = TRUE)
+## MCMC chain 1 +## burn-in +## | |======================================================================| 100% +## acceptance rate: 21.7% +## sampling phase +## | |======================================================================| 100% +## acceptance rate: 22% +## chain completed in 2.187709 seconds## total MCMC run-time: 2.19 seconds--
# coupling acceptance plot -plot_mc_acceptance(mcmc)
+
# trace plot +plot_trace(mcmc, show = "mu")+-
# extract posterior draws -output_sub <- subset(mcmc$output, phase == "sampling") -mu_draws <- output_sub$mu - -# overlay plots -hist(mu_draws, breaks = seq(-L, L, l = 201), probability = TRUE, main = "", col = "black") -lines(x, fx, col = 2, lwd = 2)# coupling acceptance plot +plot_mc_acceptance(mcmc)
++
# extract posterior draws +output_sub <- subset(mcmc$output, phase == "sampling") +mu_draws <- output_sub$mu + +# overlay plots +hist(mu_draws, breaks = seq(-L, L, l = 201), probability = TRUE, main = "", col = "black") +lines(x, fx, col = 2, lwd = 2)
@@ -242,5 +252,7 @@
+ +
- Unpack each of the parameters from the