Title: | Particle Metropolis Within Gibbs |
Version: | 0.2.7 |
Description: | Provides an R implementation of the Particle Metropolis within Gibbs sampler for model parameter, covariance matrix and random effect estimation. A more general implementation of the sampler based on the paper by Gunawan, D., Hawkins, G. E., Tran, M. N., Kohn, R., & Brown, S. D. (2020) <doi:10.1016/j.jmp.2020.102368>. An HTML tutorial document describing the package is available at https://university-of-newcastle-research.github.io/samplerDoc/ and includes several detailed examples, some background and troubleshooting steps. |
License: | GPL-3 |
URL: | https://github.com/university-of-newcastle-research/pmwg |
BugReports: | https://github.com/university-of-newcastle-research/pmwg/issues |
Depends: | R (≥ 3.6.0) |
Imports: | checkmate, coda, condMVNorm, MASS, mvtnorm, stats |
Suggests: | covr, rtdists, testthat, knitr, rmarkdown, V8 |
Encoding: | UTF-8 |
LazyData: | true |
RoxygenNote: | 7.2.3 |
VignetteBuilder: | knitr |
NeedsCompilation: | no |
Packaged: | 2024-01-31 04:20:35 UTC; gjc216 |
Author: | Gavin Cooper [aut, cre, trl] (Package creator and maintainer), Reilly Innes [aut], Caroline Kuhne [aut], Jon-Paul Cavallaro [aut], David Gunawan [aut] (Author of original MATLAB code), Guy Hawkins [aut], Scott Brown [aut, trl] (Original translation from MATLAB to R), Niek Stevenson [aut] |
Maintainer: | Gavin Cooper <gavin@gavincooper.net> |
Repository: | CRAN |
Date/Publication: | 2024-01-31 05:00:02 UTC |
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:
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
min |
The minimum of the value being updated for the progress bar |
max |
The maximum of the value being updated for the progress bar |
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
pmwgs |
The sampler object (containing random effects) with which we are working |
window_size |
The size of the window to calculate acceptance rate over |
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
sampler |
The pmwgs object containing samples to extract. |
selection |
The selection of sample types to return. |
filter |
A filter that defines which stage to draw samples from. |
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
sampler |
The sampler object to augment |
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
samples |
The subject mean samples with which we are working |
unq_vals |
The number of unique values for each subject |
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
mix_proportion |
A vector of floats between 0 and 1 and summing to 1 which give the proportion of particles to generate from the population level parameters, the individual random effects and the conditional parameters respectively |
efficient_mu |
The mu value for the efficient proposals |
efficient_sig2 |
The sigma value for the efficient proposals |
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
pmwgs |
A Particle Metropolis within Gibbs sampler which has been set up and initialised |
stage |
The sampling stage to run. Must be one of |
iter |
The number of iterations to run for the sampler. For
|
particles |
The default here is 1000 particles to be generated for each iteration, however during the sample phase this should be reduced. |
display_progress |
Display a progress bar during sampling. |
n_cores |
Set to more than 1 to use |
n_unique |
A number representing the number of unique samples to check
for on each iteration of the sampler (An initial test for the generation
of the proposal distribution). Only used during the |
epsilon |
A value between 0 and 1 that controls the extent to which the covariance matrix is scaled when generating particles from the previous random effect. The default will be chosen based on the number of random effects in the model. |
p_accept |
A value between 0 and 1 that will flexibly tune epsilon to achieve an acceptance ratio close to what you set p_accept to. The default is set at 0.8. |
mix |
A vector of floats that controls the mixture of different sources
for particles. The function |
pdist_update_n |
The number of iterations in the sample stage after which the proposal distribution will be recomputed. |
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
s |
current subject number |
samples |
A list containing previous samples |
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
x |
The current pmwgs object |
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
sampler |
The pmwgs object that we are adding the new samples to |
n_samples |
The number of new samples to increase by |
stage |
The name of the stage from which the new samples will be drawn |
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
sampler |
The pmwgs object containing the samples |
stage |
The stage, or list of stages from which you want the samples |
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
num_particles |
The total number of particles to generate using a combination of the three methods. |
mu |
A vector of means for the multivariate normal |
sig2 |
A covariate matrix for the multivariate normal |
particle |
A particle (re proposals for latent variables) |
mix_proportion |
A vector of floats between 0 and 1 and summing to 1 which give the proportion of particles to generate from the population level parameters, the individual random effects and the conditional parameters respectively |
epsilon |
Reduce the variance for the individual level samples by this factor |
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
sampler |
The pmwgs object from which to generate the new group parameters. |
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
pmwgs |
The pmwgs object for the current run. |
err_cond |
The original error condition that prompted this. |
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
pmwgs |
The sampler object that provides the parameters. |
start_mu |
An array of starting values for the group means |
start_sig |
An array of starting values for the group covariance matrix |
display_progress |
Display a progress bar during sampling |
particles |
The number of particles to generate in initialisation |
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
x |
An object to test |
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
store |
The list containing samples from which to grab the last. |
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
s |
A number - the index of the subject. For |
data |
A data.frame (or similar object) which contains the data against which the particles are assessed. The only strict requirement is that it contains a subject column named as such to allow for the splitting of the data by unique subject id. The provided log likelihood function is the only other contact with the data. |
num_particles |
The total number of particles to generate using a combination of the three methods. |
parameters |
A list containing:
|
efficient_mu |
The mu value for the efficient proposals |
efficient_sig2 |
The sigma value for the efficient proposals |
mix_proportion |
A vector of floats between 0 and 1 and summing to 1 which give the proportion of particles to generate from the population level parameters, the individual random effects and the conditional parameters respectively |
likelihood_func |
A likelihood function for calculating log likelihood
of samples. Usually provided internally in |
epsilon |
A scaling factor to reduce the variance on the distribution based on subject random effects when generating particles. |
subjects |
A list of unique subject ids in the order they appear in the data.frame |
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
pmwgs |
The pmwgs object for the current run. |
envir |
The environment of the function at this point in time. |
err_cond |
The original error condition that prompted this. |
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
mix_proportion |
A vector of floats between 0 and 1 and summing to 1 which give the proportion of particles to generate from the population level parameters, the individual random effects and the conditional parameters respectively |
num_particles |
The total number of particles to generate using a combination of the three methods. |
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
n |
number of observations |
mu |
mean vector |
covar |
covariance matrix |
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
subj |
The index of the subject where the error was detected. |
envir |
The enclosing environment of the function where the error occurred. |
err_cond |
The original error condition that prompted this. |
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
data |
A data frame containing empirical data to be modelled. Assumed
to contain at least one column called subject whose elements are unique
identifiers for each subject. Can be any of |
pars |
The list of parameter names to be used in the model |
ll_func |
A log likelihood function that given a list of parameter
values and a data frame (or other data store) containing subject data will
return the log likelihood of |
prior |
Specification of the prior distribution for the model
parameters. It should be a list with two elements, |
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
sampler |
The pmwgs object that we are relabelling samples from |
indices |
The sample iterations from burn-in to relabel |
from |
The stage you want to re-label from. Default is "burn" |
to |
The stage you want to relabel to. Default is "adapt" |
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
v |
Degrees of freedom (scalar). |
S |
Scale matrix |
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
pmwgs |
A Particle Metropolis within Gibbs sampler which has been set up and initialised |
stage |
The sampling stage to run. Must be one of |
iter |
The number of iterations to run for the sampler. For
|
particles |
The default here is 1000 particles to be generated for each iteration, however during the sample phase this should be reduced. |
display_progress |
Display a progress bar during sampling. |
n_cores |
Set to more than 1 to use |
n_unique |
A number representing the number of unique samples to check
for on each iteration of the sampler (An initial test for the generation
of the proposal distribution). Only used during the |
epsilon |
A value between 0 and 1 that controls the extent to which the covariance matrix is scaled when generating particles from the previous random effect. The default will be chosen based on the number of random effects in the model. |
p_accept |
A value between 0 and 1 that will flexibly tune epsilon to achieve an acceptance ratio close to what you set p_accept to. The default is set at 0.8. |
mix |
A vector of floats that controls the mixture of different sources
for particles. The function |
pdist_update_n |
The number of iterations in the sample stage after which the proposal distribution will be recomputed. |
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
v |
Degrees of freedom (scalar). |
S |
Inverse scale matrix |
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
par_names |
The names of each parameter as a character vector |
subject_ids |
The unique ID of each subjects as a character vector |
iters |
The number of iterations to be pre-allocated |
stage |
The stage for which the samples will be created. Should be one
of |
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
n_pars |
The number of parameters for the model |
epsilon |
The value of epsilon set in the call to 'run_stage' or NULL. |
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
stage |
The stage being run by the sampler. |
mix |
The value of mix set in the call to 'run_stage' or NULL. |
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
i |
The current iteration of the stage being run. |
stage |
The stage being run by the sampler. |
pmwgs |
The pmwgs object from which to attempt to create the proposal distribution |
pdist_update_n |
The number of iterations to run before recomputing the proposal distribution (NA to never update or for burnin/adaptation stages) |
Test that the sampler has successfully adapted
Description
Test that the sampler has successfully adapted
Usage
test_sampler_adapted(pmwgs, n_unique, i)
Arguments
pmwgs |
The full pmwgs object with all samples |
n_unique |
The number of unique samples to look for in random effects for each subject. |
i |
The number for the current iteration of the sampler |
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
sampler |
The pmwgs object that we are adding the new samples to |
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
var_matrix |
A variance matrix |
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
epsilon |
The scaling parameter for all subjects |
acc |
A boolean vector, TRUE if current sample != last sample |
p |
The target sample acceptance rate (0-1) |
i |
The current iteration for sampling |
d |
The number of parameters for the model |
alpha |
A hyperparameter for the epsilon tuning |
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
pb |
The progress bar object |
value |
The value to set the bar width to |
extra |
A value that represents the number of accepted particles |
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
var_vector |
A variance vector |
Value
The wound vector as a matrix