The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Pooling Accumulators

This vignette shows how to define a response from a pool of accumulators. Here response A is produced by a 2-of-3 pool: the response is observed once any two of A1, A2, and A3 have finished. Response B is generated by a single competing accumulator.

library(AccumulatR)

Define the model

We use three lognormal accumulators feeding a pooled response A, plus one lognormal accumulator feeding response B.

model <- race_spec() |>
  add_accumulator("A1", "lognormal") |>
  add_accumulator("A2", "lognormal") |>
  add_accumulator("A3", "lognormal") |>
  add_accumulator("B", "lognormal") |>
  add_pool("A_pool", c("A1", "A2", "A3"), k = 2L) |>
  add_outcome("A", "A_pool") |>
  add_outcome("B", "B") |>
  set_parameters(separate = list(m = TRUE, s = TRUE)) |>
  finalize_model()

true_params <- c(
  A1.m = log(0.28),
  A1.s = 0.16,
  A2.m = log(0.28),
  A2.s = 0.16,
  A3.m = log(0.28),
  A3.s = 0.16,
  B.m = log(0.28),
  B.s = 0.18
)

Simulate data

We simulate a response and response time for each trial.

set.seed(123456)

n_trials <- 1500
params_df <- build_param_matrix(model, true_params, n_trials = n_trials)

sim <- simulate(model, params_df)

data_df <- data.frame(
  trials = sim$trials,
  R = factor(sim$R),
  rt = sim$rt,
  stringsAsFactors = FALSE
)

table(data_df$R)
## 
##   A   B 
## 751 749

When response A is observed, its response time is the finishing time of the second finisher of the pool.

Estimate parameters with optim()

We estimate a shared location and spread for the three pool members, plus the location and spread of accumulator B. The spread parameters are optimized on the log scale.

prepared <- prepare_data(model, data_df)
ctx <- make_context(model)
neg_loglik <- function(theta) {
  est <- true_params
  est[c("A1.m", "A2.m", "A3.m")] <- theta[["A_pool.m"]]
  est[c("A1.s", "A2.s", "A3.s")] <- exp(theta[["log_A_pool.s"]])
  est["B.m"] <- theta[["B.m"]]
  est["B.s"] <- exp(theta[["log_B.s"]])
  params_df <- build_param_matrix(
    model,
    est,
    n_trials = n_trials
  )
  ll <- log_likelihood(ctx, prepared, params_df)
  -as.numeric(ll)
}

start <- c(
  A_pool.m = log(0.24),
  log_A_pool.s = log(0.12),
  B.m = log(0.34),
  log_B.s = log(0.12)
)

fit <- optim(
  start,
  neg_loglik,
  method = "Nelder-Mead",
  control = list(maxit = 4000, reltol = 1e-9)
)

fit_params <- c(
  A_pool.m = fit$par[["A_pool.m"]],
  A_pool.s = exp(fit$par[["log_A_pool.s"]]),
  B.m = fit$par[["B.m"]],
  B.s = exp(fit$par[["log_B.s"]])
)
target <- c(
  A_pool.m = true_params[["A1.m"]],
  A_pool.s = true_params[["A1.s"]],
  B.m = true_params[["B.m"]],
  B.s = true_params[["B.s"]]
)

data.frame(true = target, recovered = fit_params, miss = abs(target - fit_params))

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.