pmwg: Particle Metropolis Within Gibbs.

Description

The pmwg package provides a general purpose implementation of the sampling techniques outlined in Gunawan et al. (2020) doi:10.1016/j.jmp.2020.102368. The user of this package is required to provide their own log likelihood function, but given this the functions provided can estimate model parameters, the full covariance matrix and subject random effects in a hierarchical Bayesian way.

Documentation

The documentation found at contains background information and motivation for the approach used in this package and several detailed examples of the package in action. It also includes a list of common problems and associated troubleshooting steps.

User input

The user is expected to provide a data source in a format that is compatible with R data.frame methods. This data must have at least one column named ‘subject' that has a unique identifier for each subject’s data.

Additionally the user should provide a function that when given a set of parameter estimates and the data for a single subject return the log of the likelihood of that data given the parameter estimates.

The final piece of required information is a list of the names of each parameter that should be estimated. There is also the capability to provide priors on the model parameters, start points for the model parameters and covariance matrix as well as options to fine tune the sampling process

Author(s)

Maintainer: Gavin Cooper gavin@gavincooper.net (Package creator and maintainer) [translator]

Authors:

• Reilly Innes Reilly.Innes@uon.edu.au

• Caroline Kuhne caroline.kuhne@newcastle.edu.au

• Jon-Paul Cavallaro jon-paul.cavallaro@uon.edu.au

• David Gunawan dgunawan@uow.edu.au (Author of original MATLAB code)

• Guy Hawkins guy.hawkins@newcastle.edu.au

• Scott Brown scott.brown@newcastle.edu.au (Original translation from MATLAB to R) [translator]

• Niek Stevenson niek.stevenson@gmail.com

References

Gunawan, D., Hawkins, G. E., Tran, M. N., Kohn, R., & Brown, S. D. (2020). New estimation approaches for the hierarchical Linear Ballistic Accumulator model. Journal of Mathematical Psychology, 96, 102368.

See Also

Useful links:

• https://github.com/university-of-newcastle-research/pmwg

• Report bugs at https://github.com/university-of-newcastle-research/pmwg/issues

An altered version of the utils:txtProgressBar that shows acceptance rate

Description

The progress bar displays several elements, the progress visually as a bar being filled and the percentage complete as per the standard utils::txtProgressBar and additionally the average across subjects of the rate of accepting newly generated particles.

Usage

accept_progress_bar(min = 0, max = 1)

Arguments

Value

A structure matching the structure of a txtProgresBar with additional info

Return the acceptance rate for new particles across all subjects

Description

Here the acceptance rate is defined as the rate of accepting newly generated particles for individuals random effects. That is the number of samples where a newly generated particle was accepted / the number of samples.

Usage

accept_rate(pmwgs, window_size = 200)

Arguments

Value

A vector with the acceptance rate for each subject for the last X samples

Return a CODA mcmc object with the required samples

Description

Given a sampler object and a specification of the samples required, return either an individual coda mcmc object, or a list of mcmc objects.

Usage

as_mcmc(sampler, selection = "theta_mu", filter = stages)

Arguments

Value

An mcmc object or list containing the selected samples.

Selecting sample types

The values that can be chosen for the selection argument can be one of the following list:

"theta_mu"

the model parameter estimate samples

"theta_sig"

the covariance matrix estimates, returns a list of mcmc objects, one for each model parameter.

"alpha"

the random effect estimates, returns a list of mcmc objects, one for each subject.

The default value for selection is "theta_mu"

Filtering samples

The filter argument can take one of two forms:

• An integer vector, usually a sequence of integers, that must fall within the range 1:end.

• A character vector, where each element corresponds to a stage of the sampling process, i.e. one or more of "init", "burn", "adapt" or "sample".

The default value for filter is all stages.

Examples

par_estimates <- as_mcmc(sampled_forstmann)
par_estimates_sample_stage <- as_mcmc(sampled_forstmann, filter = "sample")
rand_eff <- as_mcmc(
  sampled_forstmann,
  selection = "alpha",
  filter = "sample"
)

Augment existing sampler object to have subject specific epsilon storage

Description

Older sampler object will be missing subject specific scaling parameter (epsilon) storage, and running a stage with an updated pmwg will fail. To fix this you can run the augment_sampler_epsilon function to fill the appropriate array internals with NA values

Usage

augment_sampler_epsilon(sampler)

Arguments

Value

A pmwgs sampler with epsilon array set internally

Check whether the adaptation phase has successfully completed

Description

Check whether the adaptation phase has successfully completed

Usage

check_adapted(samples, unq_vals = 20)

Arguments

Value

A boolean TRUE or FALSE depending on the result of the test

Check for efficient proposals if necessary

Description

Takes a mix proportion vector (3 x float) and the efficient proposal mu and sigma. If efficient proposals are to be used (mix_proportion[3] > 0) then test the efficient proposal values to see whether they are not null and appropriate.

Usage

check_efficient(mix_proportion, efficient_mu, efficient_sig2)

Arguments

Value

nothing, stops operation on incorrect combination of parameters.

Test the arguments to the run_stage function for correctness

Description

Takes the arguments to run_stage and checks them for completeness and correctness. Returns a list of cleaned/checked arguments to the caller.

Usage

check_run_stage_args(
  pmwgs,
  stage,
  iter,
  particles,
  display_progress,
  n_cores,
  n_unique,
  epsilon,
  subj_epsilon,
  p_accept,
  mix,
  pdist_update_n
)

Arguments

Obtain the efficent mu and sigma from the adaptation phase draws

Description

Obtain the efficent mu and sigma from the adaptation phase draws

Usage

conditional_parms(s, samples)

Arguments

Value

A list containing the conditional mean and variances for this subject

Create distribution parameters for efficient proposals

Description

From the existing samples, create a proposal distribution for drawing efficient samples from.

Usage

create_efficient(x)

Arguments

Value

A list containing the mu and sigma for the proposal distribution.

Extend the main data store with empty space for new samples

Description

Extend the main data store with empty space for new samples

Usage

extend_sampler(sampler, n_samples, stage)

Arguments

Value

The pmwgs object with the space for new samples added

Extract relevant samples from the list for conditional dist calc

Description

From the sampler, extract relevant samples for the creation of the proposal distribution.

Usage

extract_samples(sampler, stage = c("adapt", "sample"))

Arguments

Value

A list containing only appropriate samples (non init/burnin samples)

Forstmann et al.'s data

Description

A dataset containing the speed or accuracy manipulation for a Random Dot Motion experiment.

Usage

forstmann

Format

A data frame with 15818 rows and 5 variables:

subject

integer ID for each subject

rt

reaction time for each trial as a double

condition

Factor with 3 levels for Speed, Accuracy and Neutral

stim

Factor with 2 levels for Left and Right trials

resp

Factor with 2 levels for Left and Right responses

Details

Details on the dataset can be found in the following paper:

Striatum and pre-SMA facilitate decision-making under time pressure

Birte U. Forstmann, Gilles Dutilh, Scott Brown, Jane Neumann, D. Yves von Cramon, K. Richard Ridderinkhof, Eric-Jan Wagenmakers.

Proceedings of the National Academy of Sciences Nov 2008, 105 (45) 17538-17542; DOI: 10.1073/pnas.0805903105

Source

https://www.pnas.org/content/105/45/17538

Generate proposal particles

Description

Generates particles for the new_sample function

Usage

gen_particles(
  num_particles,
  mu,
  sig2,
  particle,
  ...,
  mix_proportion = c(0.5, 0.5, 0),
  prop_mu = NULL,
  prop_sig2 = NULL,
  epsilon = 1
)

Arguments

Details

Generate particles for a particular subject from a mix of population level (hierarchical) distribution, from the particles (containing individual level distribution) and/or from the conditional on accepted individual level particles, a more efficient proposal method.

This function is used in burnin, adaptation and sampling using various combinations of the arguments.

Value

The new proposals

Gibbs step of the Particle Metropolis within Gibbs sampler

Description

Samples new theta_mu and theta_sig using Gibbs sampling

Usage

gibbs_step(sampler)

Arguments

Value

A new sample for theta_mu, theta_sig and some new mixing weights in a list for use in the Particle Metropolis step.

Error handler for the gibbs_step call

Description

If an error was detected when generating new values in Gibbs step this function is called to generate the error message and save the state of the samples at that moment to help with debugging.

Usage

gibbs_step_err(pmwgs, err_cond)

Arguments

Initialise values for the random effects

Description

Initialise the random effects for each subject using MCMC.

Usage

init(
  pmwgs,
  start_mu = NULL,
  start_sig = NULL,
  display_progress = TRUE,
  particles = 100
)

Arguments

Details

Before sampling can start the Particle Metropolis within Gibbs sampler needs initial values for the random effects. The init function generates these values using a Monte Carlo algorithm. One alternative methods would be setting the initial values randomly.

Optionally takes starting values for the model parameters and the variance / covariance matrix. All arrays must match the appropriate shape.

For example, with 5 parameters and 10 subjects, the model parameter start means must be a vector of length 5 and the covariance matrix must be an array of 5 x 5.

If the start_mu and start_sig arguments are left at the default (NULL) then start_mu will be sampled from a normal distribution with mean as the prior mean for eac variable and sd as the square of the variance from the prior covariance matrix. start_sig by default is sampled from an inverse wishart (IW) distribution. For a model with the number of parameters N the degrees of freedom of the IW distribution is set to N*3 and the scale matrix is the identity matrix of size NxN.

Value

The sampler object but with initial values set for theta_mu, theta_sig, alpha and other values for the first sample.

Examples

lba_ll <- function(x, data) {
  x <- exp(x)
  if (any(data$rt < x["t0"])) {
    return(-1e10)
  }
  sum(
    log(
      rtdists::dLBA(
        rt = data$rt,
        response = data$correct,
        A = x["A"],
        b = x["A"] + x[c("b1", "b2", "b3")][data$condition],
        t0 = x["t0"],
        mean_v = x[c("v1", "v2")],
        sd_v = c(1, 1),
        silent = TRUE
      )
    )
  )
}
sampler <- pmwgs(
  forstmann,
  c("b1", "b2", "b3", "A", "v1", "v2", "t0"),
  lba_ll
)
sampler <- init(sampler, particles=10)

Test whether object is a pmwgs

Description

Checks whether object is a Particle Metropolis with Gibbs sampler

Usage

is.pmwgs(x)

Arguments

Value

logical, whether object inherits from pmwgs

Examples

if (is.pmwgs(sampled_forstmann)) {
  print("sampled_forstmann object is a pmwgs")
}

Create a list with the last samples in the pmwgs object

Description

Create a list with the last samples in the pmwgs object

Usage

last_sample(store)

Arguments

Value

A list containing the last sample of group mean and variance and subject means.

Generate particles and select one to be the new sample

Description

Generate a new sample for a particular subject given their data and the new model parameter estimates. This should not be called directly, rather it is used internally to run_stage.

Usage

new_sample(
  s,
  data,
  num_particles,
  parameters,
  efficient_mu = NULL,
  efficient_sig2 = NULL,
  mix_proportion = c(0.5, 0.5, 0),
  likelihood_func = NULL,
  epsilon = 1,
  subjects = NULL
)

Arguments

Details

The function that controls the generation of new samples for the Particle Metropolis within Gibbs sampler. It generates samples from either the initial proposal or from the last iteration of the sampler. This function should not usually need to be called, as the run_stage function uses this internally.

The way it selects a new sample is by generating proposal particles from up to three different distributions (according to a mixing proportion).

The first distribution is based on the current model parameter sample values. The second distribution is based on the last random effects for the subject. The third distribution is only used in the final sampling phase and is based on the conditional distribution built from accepted particles in the adapt phase of the sampler.

Value

A single sample from the new proposals

Error handler forany error in new_sample function call(s)

Description

If an error was detected when generating new samples. Save the state of the samples and particles at that moment to help with debugging.

Usage

new_sample_err(pmwgs, envir, err_cond)

Arguments

Check and normalise the number of each particle type from the mix_proportion

Description

Takes a mix proportion vector (3 x float) and a number of particles to generate and returns a vector containing the number of each particle type to generate.

Usage

numbers_from_proportion(mix_proportion, num_particles = 1000)

Arguments

Value

The wound vector as a matrix

Generate a cloud of particles from a multivariate normal distribution

Description

Takes the mean and variance for a multivariate normal distribution, as well as the number of particles to generate and return random draws from the multivariate normal if the numbers of particles is > 0, otherwise return NULL. At least one of mean or sigma must be provided.

Usage

particle_draws(n, mu, covar)

Arguments

Value

If n > 0 returns n draws from the multivariate normal with mean and sigma, otherwise returns NULL

Error handler for the particle selection call

Description

If an error was detected when selecting the winning particle, save the state of the samples and particles at that moment to help with debugging.

Usage

particle_select_err(subj, envir, err_cond)

Arguments

Create a PMwG sampler and return the created object

Description

This function takes a few necessary elements for creating a PMwG sampler. Each pmwgs object is required to have a dataset, a list of parameter names, a log likelihood function and optionally a prior for the model parameters.

Usage

pmwgs(data, pars, ll_func, prior = NULL)

Arguments

Value

A pmwgs object that is ready to be initialised and sampled.

Examples

# Specify the log likelihood function
lba_loglike <- function(x, data) {
  x <- exp(x)
  if (any(data$rt < x["t0"])) {
    return(-1e10)
  }
  # This is faster than "paste".
  bs <- x["A"] + x[c("b1", "b2", "b3")][data$condition]

  out <- rtdists::dLBA(
    rt = data$rt, # nolint
    response = data$stim,
    A = x["A"],
    b = bs,
    t0 = x["t0"],
    mean_v = x[c("v1", "v2")],
    sd_v = c(1, 1),
    distribution = "norm",
    silent = TRUE
  )
  bad <- (out < 1e-10) | (!is.finite(out))
  out[bad] <- 1e-10
  out <- sum(log(out))
  out
}

# Specify parameter names and priors
pars <- c("b1", "b2", "b3", "A", "v1", "v2", "t0")
priors <- list(
  theta_mu_mean = rep(0, length(pars)),
  theta_mu_var = diag(rep(1, length(pars)))
)

# Create the Particle Metropolis within Gibbs sampler object
sampler <- pmwgs(
  data = forstmann,
  pars = pars,
  ll_func = lba_loglike,
  prior = priors
)

sampler = init(sampler, particles=10)
sampler = run_stage(sampler, stage="burn", iter=10, particles=10)

Relabel requested burn-in samples as adaptation

Description

Given a sampler object and a vector of sample indices, relabel the given samples to be adaptation samples. This will allow them to be used in the calculation of the conditional distribution for efficient sampling.

Usage

relabel_samples(sampler, indices, from = "burn", to = "adapt")

Arguments

Value

The pmwgs object with re-labelled samples

Further information

This should not usually be needed, however if you have a model that is slow to fit, and upon visual inspection and/or trace analysis you determine that during burn-in the samples had already approached the posterior distribution then you can use this function to re-label samples from that point onwards to be classed as adaptation samples.

This will allow them to be used in tests that check for the number of unique samples, and in the building of the conditional distribution (which is used for efficient sampling)

If all old samples do not match 'from' then an error will be raised.

Examples

new_pmwgs <- relabel_samples(sampled_forstmann, 17:21)

The Inverse Wishart Distribution

Description

Random generation from the Inverse Wishart distribution.

Usage

riwish(v, S)

Arguments

Details

The mean of an inverse Wishart random variable with v degrees of freedom and scale matrix S is (v-p-1)^{-1}S.

Value

riwish generates one random draw from the distribution.

Examples

## Not run: 
draw <- riwish(3, matrix(c(1,.3,.3,1),2,2))

## End(Not run)

Run a stage of the PMwG sampler

Description

Run one of burnin, adaptation or sampling phase from the PMwG sampler. Each stage involves slightly different processes, so for the full PMwG sampling we need to run this three times.

Usage

run_stage(
  pmwgs,
  stage,
  iter = 1000,
  particles = 100,
  display_progress = TRUE,
  n_cores = 1,
  n_unique = ifelse(stage == "adapt", 100, NA),
  epsilon = NULL,
  p_accept = 0.8,
  mix = NULL,
  pdist_update_n = ifelse(stage == "sample", 50, NA)
)

Arguments

Details

The burnin stage by default selects 50 parameter sample (selected through a Gibbs step) and 50 the previous random effect of each subject. It assesses each particle with the log-likelihood function and samples from all particles weighted by their log-likelihood.

The adaptation stage selects and assesses particle in the same was as burnin, however on each iteration it also checks whether each subject has enough unique random effect samples to attempt to create a conditional distribution for efficient sampling in the next stage. If the attempt at creating a conditional distribution fails, then the number of unique samples is increased and sampling continues. If the attempt succeeds then the stage is finished early.

The final stage (sampling) by default samples predominantly from the conditional distribution created at the end of adaptation. This is more efficient and allows the number of particles to be reduced whilst still getting a high enough acceptance rate of new samples.

Once complete each stage will return a sampler object with the new samples stored within it.

The progress bar (which is displayed by default) shows the number of iterations out of those requested which have been completed. It also contains additional information at the end about the number of newly generated particles that have been accepted. This is show as New(XXX average across subjects of newly sampled random effects accept rate. See accept_rate for more detail on getting individual accept rate values per subject.

Value

A pmwgs object with the newly generated samples in place.

Examples

library(rtdists)
sampled_forstmann$data <- forstmann
run_stage(sampled_forstmann, "sample", iter = 1, particles = 10)

The Wishart Distribution

Description

Random generation from the Wishart distribution.

Usage

rwish(v, S)

Arguments

Details

The mean of a Wishart random variable with v degrees of freedom and inverse scale matrix S is vS.

Value

rwish generates one random draw from the distribution.

Examples

## Not run: 
draw <- rwish(3, matrix(c(1,.3,.3,1),2,2))

## End(Not run)

Create a new list for storage samples in the pmwgs object

Description

Create a new list for storage samples in the pmwgs object

Usage

sample_store(par_names, subject_ids, iters = 1, stage = "init")

Arguments

Value

A list containing the conditional mean and variances for this subject

A sampled object of a model of the Forstmann dataset

Description

A pmwgs object with a limited number of samples of the Forstmann dataset.

Usage

sampled_forstmann

Format

A pmwgs object minus the data. A pmwgs opbject is a list with a specific structure and elements, as outlined below.

par_names

A character vector containing the model parameter names

n_pars

The number of parameters in the model

n_subjects

The number of unique subject ID's in the data

subjects

A vector containing the unique subject ID's

prior

A list that holds the prior for theta_mu (the model parameters). Contains the mean (theta_mu_mean), covariance matrix (theta_mu_var) and inverse covariance matrix (theta_mu_invar)

ll_func

The log likielihood function used by pmwg for model estimation

samples

A list with defined structure containing the samples, see the Samples Element section for more detail

Details

The pmwgs object is missing one aspect, the pmwgs$data element. In order to fully replicate the full object (ie to run more sampling stages) you will need to add the data back in, via sampled_forstmann$data <- forstmann

Samples Element

The samples element of a PMwG object contains the different types of samples estimated by PMwG. These include the three main types of samples theta_mu, theta_sig and alpha as well as a number of other items which are detailed here.

theta_mu

samples used for estimating the model parameters (group level), an array of size (n_pars x n_samples)

theta_sig

samples used for estimating the parameter covariance matrix, an array of size (n_pars x n_pars x n_samples)

alpha

samples used for estimating the subject random effects, an array of size (n_pars x n_subjects x n_samples)

stage

A vector containing what PMwG stage each sample was drawn in

subj_ll

The winning particles log-likelihood for each subject and sample

a_half

Mixing weights used during the Gibbs step when creating a new sample for the covariance matrix

last_theta_sig_inv

The inverse of the last samples covariance matrix

idx

The index of the last sample drawn

Source

https://www.pnas.org/content/105/45/17538

Set default values for epsilon

Description

Takes the number of parameters and the epsilon arg value and sets a default if necessary.

Usage

set_epsilon(n_pars, epsilon)

Arguments

Set default values for mix

Description

Takes the current stage and the mix arg value and sets a default if necessary.

Usage

set_mix(stage, mix)

Arguments

Setup the proposal distribution arguments (if in sample stage)

Description

Takes the current stage and the mix arg value and sets a default if necessary.

Usage

set_proposal(i, stage, pmwgs, pdist_update_n)

Arguments

Test that the sampler has successfully adapted

Description

Test that the sampler has successfully adapted

Usage

test_sampler_adapted(pmwgs, n_unique, i)

Arguments

Value

A list containing a string representing successful/unsuccessful adaptation and an optional message. The string representing the success or failure can be one of c("success", "continue", "increase")

Trim the unneeded NA values from the end of the sampler

Description

Trim the unneeded NA values from the end of the sampler

Usage

trim_na(sampler)

Arguments

Value

The pmwgs object without NA values added during extend_sampler

Unwinds variance matrix to a vector

Description

Takes a variance matrix and unwind to a vector via Cholesky decomposition then take the log of the diagonal.

Usage

unwind(var_matrix, ...)

Arguments

Value

The unwound matrix as a vector

Update the subject specific scaling parameters (epsilon)

Description

Update the subject specific scaling parameter according to procedures outlined in P. H. Garthwaite, Y. Fan & S. A. Sisson (2016) Adaptive optimal scaling of Metropolis–Hastings algorithms using the Robbins–Monro process, Communications in Statistics - Theory and Methods, 45:17, 5098-5111, DOI: 10.1080/03610926.2014.936562

Usage

update_epsilon(epsilon, acc, p, i, d, alpha)

Arguments

Value

A vector with the new subject specific epsilon values

A function that updates the accept_progress_bar with progress and accept rate

Description

A function that updates the accept_progress_bar with progress and accept rate

Usage

update_progress_bar(pb, value, extra = 0)

Arguments

Value

The old value that was present before updating

Winds a variance vector back to a vector

Description

The reverse of the unwind function, takes a variance vector and windows back into matrix

Usage

wind(var_vector, ...)

Arguments

Value

The wound vector as a matrix