## ----setup, include = FALSE---------------------------------------------------
library(XSLmodels)

knitr::opts_chunk$set(
  echo = TRUE,
  collapse = TRUE,
  comment = "#>",
  warning = FALSE
)

## ----data---------------------------------------------------------------------
labels <- vapply(xsl_datasets, `[[`, "", "label")
conds <- xsl_datasets[labels %in% c("201", "212", "215", "220", "filt3E_6L", "freq369-3x3hiCD")]
# DEoptim advises a population of at least 10 per parameter; this one is
# smaller to keep the vignette fast, so its warnings about that are hidden
de <- DEoptim::DEoptim.control(NP = 15, itermax = 8, trace = FALSE)
models <- c("baseline", "decay", "uncfam", "rescorla_wagner")

## ----sse----------------------------------------------------------------------
sse_fits <- lapply(models, \(m) get_group_model_fit(m, datasets = conds,
                                                    deoptim_control = de, seed = 1))
names(sse_fits) <- models
sapply(sse_fits, \(f) f$fit_result[[1]]$optim$bestval)

## ----binomial-----------------------------------------------------------------
binom <- xslControl(objective = "binomial")
ll_fits <- lapply(models, \(m) get_group_model_fit(m, datasets = conds, control = binom,
                                                   deoptim_control = de, seed = 1,
                                                   fit_lapse = TRUE))
names(ll_fits) <- models
ic <- \(fits) do.call(rbind, lapply(fits, \(f) {
  r <- f$fit_result[[1]]
  data.frame(k = r$k, loglik = r$loglik, AIC = r$aic, BIC = r$bic,
             lapse = unname(r$optim$bestmem["lapse"]))
}))
ic(ll_fits)

## ----foils--------------------------------------------------------------------
foil_cor <- function(mod, data) {
  pairs <- do.call(rbind, lapply(xsl_run(mod, data)$fits, \(f) {
    ch <- choice_probabilities(f$matrix, f$data)
    do.call(rbind, lapply(seq_len(nrow(ch)), \(i) {
      j <- match(ch$correct[i], ch$alternatives[[i]])
      q <- ch$p[[i]][-j]
      h <- ch$counts[[i]][-j]
      if (sum(h) == 0 || sum(q) == 0) return(NULL)
      data.frame(model = q / sum(q), people = h / sum(h))
    }))
  }))
  cor(pairs$model, pairs$people)
}
best <- \(m) ll_fits[[m]]$fit_result[[1]]$optim$bestmem
c(baseline = foil_cor(baseline(), conds),
  decay = foil_cor(decay(C = best("decay")[["C"]]), conds),
  uncfam = foil_cor(uncfam(X = best("uncfam")[["X"]], B = best("uncfam")[["B"]],
                           C = best("uncfam")[["C"]]), conds))

## ----multinomial--------------------------------------------------------------
multi <- xslControl(objective = "multinomial")
mn_fits <- lapply(c("decay", "uncfam"), \(m) get_group_model_fit(m, datasets = conds, control = multi,
                                                                  deoptim_control = de, seed = 1,
                                                                  fit_lapse = TRUE))
names(mn_fits) <- c("decay", "uncfam")
ic(mn_fits)

## ----cv-----------------------------------------------------------------------
cv <- lapply(c("decay", "uncfam"), \(m) get_crossvalidated_model_fit(
  m, n_folds = 3, datasets = conds, control = binom, deoptim_control = de,
  seed = 1, fit_lapse = TRUE)$cv_result)
names(cv) <- c("decay", "uncfam")
sapply(cv, \(r) c(held_out_nll = sum(r$test_nll), held_out_sse = mean(r$test_sse)))