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.
pmwgs(data, pars, ll_func, prior = NULL)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 data.frame,
data.table or tibble, or any other data frame like object
that can have subsets created in an identical way.
The list of parameter names to be used in the model
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 data given x.
Specification of the prior distribution for the model
parameters. It should be a list with two elements, theta_mu_mean and
theta_mu_var which fully specify the prior distribution. If left as
the default (NULL) then the theta_mu_mean will be all zeroes and
theta_mu_var will be 1 on the diagonal and zero elsewhere.
A pmwgs object that is ready to be initialised and sampled.
# 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)
#> MESSAGE: Sampling Initial values for random effects
#>
|
| | 0%
|
|==== | 5%
|
|======= | 11%
|
|=========== | 16%
|
|=============== | 21%
|
|================== | 26%
|
|====================== | 32%
|
|========================== | 37%
|
|============================= | 42%
|
|================================= | 47%
|
|===================================== | 53%
|
|========================================= | 58%
|
|============================================ | 63%
|
|================================================ | 68%
|
|==================================================== | 74%
|
|======================================================= | 79%
|
|=========================================================== | 84%
|
|=============================================================== | 89%
|
|================================================================== | 95%
|
|======================================================================| 100%
sampler = run_stage(sampler, stage="burn", iter=10, particles=10)
#> MESSAGE: Epsilon has been set to 0.5 based on number of parameters
#> MESSAGE: mix has been set to c(0.5, 0.5, 0) based on the stage being run
#> Phase 1: Burn in
#>
| | 0% | New( 0%)
|====== | 10% | New( 0%)
|============ | 20% | New( 0%)
|================== | 30% | New( 97%)
|======================== | 40% | New( 95%)
|============================== | 50% | New( 96%)
|==================================== | 60% | New( 96%)
|========================================== | 70% | New( 96%)
|================================================ | 80% | New( 95%)
|====================================================== | 90% | New( 95%)
|============================================================| 100% | New( 95%)