pmwgs: Create a PMwG sampler and return the created object

View source: R/particle_sampler.R

pmwgsR Documentation

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 data.frame, data.table or tibble, or any other data frame like object that can have subsets created in an identical way.

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 data given x.

prior

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.

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)

pmwg documentation built on May 29, 2024, 12:09 p.m.