---
title: "Pooling Accumulators"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Pooling Accumulators}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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.

```{r}
library(AccumulatR)
```

## Define the model

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

```{r}
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.

```{r}
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)
```

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.

```{r eval=FALSE}
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))
```