# Validation and Comparison of Non-Stationary Cognitive Models: A Diffusion
Model Application

Lukas Schumacher
Institute of Psychology
Heidelberg University
&Martin Schnuerch
Institute of Psychology
University of Mannheim
&Andreas Voss
Institute of Psychology
Heidelberg University
&Stefan T. Radev
Department of Cognitive Science
Rensselaer Polytechnic Institute
For correspondence, please contact Lukas Schumacher<EMAIL_ADDRESS>
###### Abstract

Cognitive processes undergo various fluctuations and transient states across
different temporal scales. Superstatistics are emerging as a flexible
framework for incorporating such non-stationary dynamics into existing
cognitive model classes. In this work, we provide the first experimental
validation of superstatistics and formal comparison of four non-stationary
diffusion decision models in a specifically designed perceptual decision-
making task. Task difficulty and speed-accuracy trade-off were systematically
manipulated to induce expected changes in model parameters. To validate our
models, we assess whether the inferred parameter trajectories align with the
patterns and sequences of the experimental manipulations. To address
computational challenges, we present novel deep learning techniques for
amortized Bayesian estimation and comparison of models with time-varying
parameters. Our findings indicate that transition models incorporating both
gradual and abrupt parameter shifts provide the best fit to the empirical
data. Moreover, we find that the inferred parameter trajectories closely
mirror the sequence of experimental manipulations. Posterior re-simulations
further underscore the ability of the models to faithfully reproduce critical
data patterns. Accordingly, our results suggest that the inferred non-
stationary dynamics may reflect actual changes in the targeted psychological
constructs. We argue that our initial experimental validation paves the way
for the widespread application of superstatistics in cognitive modeling and
beyond.

## Introduction

The human brain operates in a perpetual state of activity, whether it is
focused on a particular task or wandering in the inner world of thoughts. This
activity reflects the non-stationary nature of neuronal dynamics, which are
characterized by a complex interplay between transient, evoked states, and
ongoing spontaneous fluctuations [1, 2]. The complex cognitive processes that
emerge from this neuronal activity also tend to exhibit non-stationary
dynamics [3, 4, 5, 6]. In other words, proverbial cognitive processes, such as
attention, memory, and decision-making, are not constant over time, but
instead undergo fluctuations, shifts, and alterations in their functions.

Lapses of attention are a canonical cause of such non-stationary dynamics.
Even when actively engaged in a task, our focus can drift or momentarily
falter [7]. Moreover, our capacity to sustain attention and concentrate may
vary, influenced by factors such as fatigue, motivation, and external
distractions [8, 9, 10]. These fluctuations can have a significant impact on
our cognitive functioning, but they are often overlooked or simplified in
traditional models of cognition. And while these often assume cognitive
processes to be stable and time-invariant, there has been a growing
recognition that traditional models do not fully capture the complexity and
variability of real-world cognition [11, 12, 13, 14, 15, 16, 17]. Common
approaches to address variability in the components of cognitive models can be
broadly classified into four categories: stationary variability, trial
binning, regression approach, and frontend-backend models.

Figure 1:  A conceptual illustration of the five main approaches to model
temporal variation in the parameters $\theta$ of a cognitive model
$\mathcal{G}$. A Stationary variability, also known as inter-trial
variability, assumes that parameter values fluctuate around a stable mean. B
Trial binning involves organizing the data into distinct bins and fitting a
cognitive model $\mathcal{G}$ to each bin individually. C Regression approach
employs time (and sometimes additional contextual variables) as predictors for
the parameters $\theta$. D Frontend-backend models employ a mechanistic model,
referred to as the frontend, to elucidate the dynamics of the parameter of the
cognitive model (i.e., the backend). E Superstatistics involve a superposition
of multiple stochastic processes operating on different temporal scales. They
comprise a low-level observation model $\mathcal{G}$ and a high-level
transition model $\mathcal{T}$ that specifies how the parameters $\theta_{t}$
evolve stochastically.

The first approach assumes random fluctuations around a stable mean, referred
to as stationary variability (see Figure 1A). A prominent example of this
approach is the “full” diffusion decision model (DDM), which allows for inter-
trial variability of its core parameters [18, 19]. However, stationary inter-
trial variability mainly improves in-sample model fit and cannot identify
systematic changes or sudden shifts in core model parameters. Moreover, the
resulting model family still treats behavioral data as independent and
identically distributed (IID) responses, making it unsuitable for
investigating systematic changes in cognitive constructs.

Another approach for detecting systematic changes in cognitive model
components is trial binning [13, 20, 21]. This method involves organizing data
into discrete bins and then applying a stationary model to each of these data
subsets separately (see Figure 1B). One can then examine variations in
parameter estimates across these bins. The challenge in employing this
approach is the selection of the number of time steps within each bin, which
introduces an unwelcome trade-off between temporal resolution and estimation
quality. For instance, if only a few time steps are chosen, the analysis can
yield relatively fine-grained, but very uncertain estimates due to the low
number of data points. A further shortcoming of trial binning is that
estimates within a specific bin are not informed by data from neighboring
bins. However, the appeal of dynamic modeling lies in the distinctive
capability to utilize both past and future data to constrain the estimated
parameter trajectories.

The third approach involves a generalized linear model (GLM) with time (and
possibly other contextual factors) as a predictor of model parameters [17,
22]. The GLM approach is more appealing than trial binning, as it can detect
linear or non-linear changes in model parameters without loss of resolution
(see Figure 1C). However, the underlying regression function makes strong
assumptions about the nature of the relationship between model parameters and
time. Thus, even though a modeler will typically fit and compare a few
plausible specifications (e.g., linear vs. exponential), it is often difficult
to determine all plausible specifications a priori, and the overall
flexibility of the GLM model as a process characterization remains severely
limited [16].

Differently, the frontend-backend approach aims to account for changes in
model parameters, while providing a mechanistic explanation for the dynamic
nature of the target system (see Figure 1D). Here, the backend model pertains
to the cognitive model which formalizes how the behavioral data is generated
(e.g., a DDM). The frontend constitutes a mechanistic model, elucidating how
the parameters of the backend model adapt over time, in different contexts and
in response to additional factors [23, 24, 25, 26]. This approach has several
advantages, as it not only accommodates the dynamic nature of the parameters,
but also provides a mechanistic description for their temporal variation
through a set of static parameters and deterministic functions. For instance,
there has been a recent trend to use reinforcement learning models as a
frontend model to inform changes in DDM parameters due to reward-based
learning [27, 28, 23]. Nevertheless, detailed frontend models are often
challenging to develop, estimate, and compare.

Recently, we proposed an alternative approach that infers non-stationary
parameter trajectories directly from the data, while imposing minimal
constraints on how parameters change over time [11]. Our approach leverages a
framework known as superstatistics [29, 30, 31], which involves a
superposition of multiple stochastic processes operating on distinct time
scales (see Figure 1E). At its core, this model comprises a low-level
observation model and a high-level transition model. The former describes how
data at a specific time point is generated, akin to the backend model. Like
the frontend model, the transition model characterizes how the parameters
change over time. However, from a superstatistics perspective, the transition
model is inherently a stochastic process, exemplified, for instance, by a
Gaussian random walk or a regime switching process.

The superstatistics approach effectively addresses the limitations of prior
methodologies. Unlike stationary models, superstatistical models can readily
generate non-stationary variations in the parameters of the low-level model,
facilitating gradual or sudden transitions between different states.
Furthermore, parameter estimates are contingent on past data points, thereby
treating the data no longer as IID. In contrast to the trial-binning approach,
models within the superstatistics framework leverage the entirety of available
data, mitigating concerns about insufficient data points for parameter
estimation. Different from GLM approaches, our superstatistics method imposes
minimal assumptions on potential parameter trajectories, making it
significantly less restrictive.

In contrast to frontend-backend models, superstatistics do not offer
mechanistic explanations for parameter dynamics but provide greater
flexibility in their estimation. Although mechanistic explanations are central
to psychological research, there are cases where suitable explanations are
lacking or are applicable only to specific parameters. Therefore, we consider
these two approaches as complementary. The superstatistical framework takes a
bottom-up, exploratory approach, functioning as a tool for generating
hypotheses. In subsequent stages, one could potentially formulate plausible
frontend models based on insights from parameter trajectories inferred with a
superstatistical model. Additionally, superstatistical models can serve as
benchmarks for testing and validating competing frontend-backend models by
comparing resulting parameter trajectories from both methods.

Having laid out the potential benefits of the superstatistics framework and
its applicability in the realm of cognitive process models [11], a pivotal
question arises: Do the inferred parameter trajectories genuinely reflect
shifts in the cognitive constructs they aim to represent, or are they merely a
modeling artefact? To address this inquiry, we embark on an experimental
validation study. In this study, we manipulate the experimental context in a
manner that allows us to confidently anticipate how individuals and,
consequently, their inferred cognitive constructs, will respond. In other
words, if the inferred parameter time series mirror the alterations in the
experimental context, we garner substantial evidence that these trajectories
indeed reflect changes in the psychological constructs.

Throughout, we employ the well-established $4$-parameter DDM [32] as a low-
level observation model. The DDM is a mathematical model that simultaneously
accounts for response time (RT) and choice data obtained from two-alternative
decision tasks. Fundamentally, it posits that, in forced-choice binary
decision task, individuals accumulate evidence for the decision alternatives
until a certain threshold is met, triggering a decision. Each of the DDM’s
four core parameters corresponds to a specific psychological construct: (i)
the drift rate $v$ signifies the average speed of information uptake; (ii) the
threshold $a$ serves as a proxy for decision caution; (iii) the relative
starting point $\beta$ represents a priori decision preferences; and (iv) the
additional constant $\tau$ accounts for the duration of all processes taking
place prior and following a decision, such as stimulus encoding or motor
action [33, but see].

A primary reason for our choice of the DDM as the observation model lies in
its rigorous prior validation [34, 35, 36]. These prior studies have
convincingly demonstrated that the DDM’s parameters are valid reflections of
the intended psychological constructs. Moreover, the manipulation of
experimental conditions leading to systematic alterations in specific DDM
parameters is well-documented and comprehensively understood [37]. For
example, varying the difficulty of an experimental task alters the drift rate
parameter, whereas providing verbal instructions to prioritize either speed or
accuracy during task-solving leads to observable shifts in the threshold
parameter and sometimes also in the non-decision time [38].

In this study, we focus on the aforementioned experimental manipulations
targeting the drift rate and the threshold parameters. We employed a color
discrimination task, which was also utilized in the validation study by [34].
During this task, individuals must decide whether there are more blue or more
orange pixels in a patch of pixels. The difficulty of the task can be easily
manipulated by adjusting the ratio of blue and orange pixels. The farther the
ratio is from 1:1, the easier the task becomes. Additionally, we manipulated
the emphasis on speed or accuracy by verbally instructing participants to
prioritize one over the other.

Systematic changes in cognitive model parameter can appear in different ways,
ranging from changing slowly and gradually to more rapid and large shifts. In
our experiment, we focus on two different types. Firstly, task difficulty
changes frequently to the next easier or harder level, imitating gradual
changes. Secondly, the speed-accuracy emphasis changes less regularly after
each trial block, resembling sudden shifts. The primary aim of our experiment
is to investigate whether the parameter trajectories inferred with a non-
stationary DDM (NSDDM) match these changing patterns of the experimental
conditions. Specifically, we expect the drift rate parameter to mirror the
gradual changes of the task difficulty. Additionally, the threshold parameter
should show sudden shifts when the priority switches between speed and
accuracy. It is crucial to understand that in this application, the NSDDM does
not have information about the experimental context and has to infer the
parameter trajectory solely on the behavioral data.

When dealing with different types of fluctuations, another crucial question
arises: What kind of transition model is most suitable for capturing the
expected dynamics? To address this question, we implement different NSDDMs
that solely differ in their transition model for the drift rate and the
threshold parameter. Specifically, we compare four different transition
models: (i) a Gaussian random walk; (ii) a mixture between a Gaussian random
walk and uniformly distributed regime changes; (iii) a Lévy flight; and (iv) a
regime switching function, where parameters either remain the same as in the
previous time step or shift uniformly. These four transition models differ in
their complexity (i.e., number of high-level parameters) and their ability to
account for different types of temporal shifts.

Performing Bayesian model comparison and parameter estimation with
superstatistical models can be computationally challenging [11]. Therefore, we
employ simulation-based inference [39, SBI,] as implemented in the BayesFlow
framework [40]. BayesFlow enables us to carry out a principled Bayesian
workflow utilizing simulation-based calibration [41, 42, SBC,] and other
validation methods [43, 44] that would otherwise be excessively time-
consuming. The contributions of the present study can be summarized as
follows:

  1. 1.

We perform an experimental validation of different non-stationary
instantiations of the diffusion decision model.

  2. 2.

We propose an amortized method for Bayesian model comparison of non-stationary
models via deep ensembles.

  3. 3.

We showcase the potential of amortized Bayesian inference for increasing the
aspirations of cognitive modeling.

## Materials and Methods

### Participants

A total of $14$ participants ($9$ female, $5$ male) were recruited for the
experiment. The participants had an average age of $23.14$ years
($\text{SD}=1.29$, $\text{Range}=[22,26]$). Every individual provided informed
consent to participate in the study, and the research protocol received
approval from the local ethics committee. The entire study was conducted in
accordance with the ethical principles outlined in the Helsinki Declaration.

### Task

The participants completed a total of $800$ trials in a color discrimination
task, including $32$ practice trials. In each trial, individuals were
presented with a rectangular patch containing blue and orange pixels and had
to determine whether there were more blue or orange pixels. Prior to the patch
presentation, a fixation cross was displayed for $300$ ms. All stimuli were
presented on a gray background.

Task difficulty was manipulated by varying the proportion of blue/orange
pixels in the patch. The following ratios were utilized: $50.5$:$49.5$;
$52.25$:$47.75$; $53.5$:$46.5$; and $55$:$45$. Half of the trials featured
orange as the dominant color, while the other half featured blue. The
difficulty level remained constant for either $8$ or $16$ trials before
transitioning to the next level of difficulty.

In addition to manipulating task difficulty, participants received two types
of instructions which changed every $48$ trials. In the “accuracy” condition,
individuals were instructed to prioritize accuracy in their responses.
Conversely, in the “speed” condition, participants were directed to emphasize
speed while maintaining a reasonable level of accuracy. Feedback was provided
after each trial to make participants aware of their performance: a green
cross for correct responses, a red minus for incorrect responses, and a red
clock for responses slower than 700 ms in the speed condition.

### Superstatistics Framework

To represent non-stationary changes in DDM parameters, we adopt a
superstatistics framework [29, 31]. Within this framework, each generative
model comprises (at least) a low-level observation model $\mathcal{G}$
characterized by time-dependent local parameters $\theta_{t}\in\mathbb{R}^{K}$
that vary according to a high-level transition model $\mathcal{T}$ with static
high-level parameters $\eta\in\mathbb{R}^{D}$. These models simulate
parameters and observable data $x_{t}\in\mathcal{X}$ according to the
following general recurrent system

$\displaystyle\theta_{t}$
$\displaystyle=\mathcal{T}(\theta_{0:t-1},\eta,\xi_{t})\quad\,\text{with}\quad\xi_{t}\sim
p(\xi\,|\,\eta),\,\,\theta_{0}\sim p(\theta)$ (1) $\displaystyle x_{t}$
$\displaystyle=\mathcal{G}(x_{1:t-1},\theta_{t},z_{t})\quad\text{with}\quad
z_{t}\sim p(z\,|\,\theta_{t}),$

where $\mathcal{T}$ represents an arbitrary high-level transition function
parameterized by $\eta$, and $\mathcal{G}$ is a (non-linear) transformation
that encapsulates the functional assumptions of the low-level model. The
random variates $\xi_{t}$ and $z_{t}$ govern the stochastic nature of the two
model components through noise outsourcing. The initial parameter
configuration $\theta_{0}$ adheres to a prior distribution $\theta_{0}\sim
p(\theta)$ encoding the available information about feasible starting
parameter values.

The above formulation is very abstract and general, highlighting the
flexibility of the superstatistics framework. Moreover, it does not assume
that the corresponding transition or likelihood densities, given by

$\displaystyle\mathbb{T}(\theta_{t}\,|\,\eta,\theta_{0:t-1})$
$\displaystyle=\int p(\theta_{t},\xi\,|\,\eta,\theta_{0:t-1})\,d\xi$
$\displaystyle(\text{implied transition density})$ (2) $\displaystyle
p(x_{t}\,|\,\theta_{t},x_{1:t-1})$ $\displaystyle=\int
p(x_{t},z\,|\,\theta_{t},x_{1:t-1})\,dz$ $\displaystyle(\text{implied
likelihood density}),$ (3)

are tractable or available in closed-form, situating our approach in the
context of simulation-based inference [39, SBI,]. Here, we build on SBI with
neural networks [45, 46, 47] as a principled approach to perform fully
Bayesian inference by using only samples from the generative system defined by
Equation 1. Importantly, our estimation methods overcome key limitations of
previous approaches related to the curse of dimensionality [31].

### Low-Level Model

In this work, we use the same standard DDM implementation as a low-level
observation model $\mathcal{G}$ for all NSDDMs. The low-level dynamics of the
evidence accumulation process are described by the following stochastic
ordinary differential equation:

$\displaystyle\mathrm{d}x_{n}=v\mathrm{d}t_{s}+z\sqrt{\mathrm{d}t_{s}}\quad\text{with}\quad
z\sim\mathcal{N}(0,1)$ (4)

Accordingly, the evidence $x_{n}$ on a given trial $n$ follows a random walk
with drift $v$ and Gaussian noise $z$, where $t_{s}$ represents time on a
continuous time scale. The core assumption of the DDM is that evidence is
accumulated with a fixed rate $v$ until one of two thresholds, $a$ or $0$, is
reached, and the corresponding decision $D_{n}$ is made:

$\displaystyle D_{n}=\begin{cases}1,&\text{if }x_{n}\geq a\\\ 0,&\text{if
}x_{n}\leq 0\end{cases}.$ (5)

Furthermore, the DDM incorporates an additive constant $\tau$, which
represents the time allocated to all non-decisional processes (i.e., stimulus
encoding and motor action). Consequently, the DDM encompasses three distinct
free parameters, namely $\theta=(v,a,\tau)$. We fixed the starting point of
the evidence accumulation process at $a/2$, since, in our case, the two
boundaries of the accumulation process correspond to correct and incorrect
responses, respectively. Thus, it is unwarranted to estimate any potential a
priori bias towards either of these boundaries [48].

Figure 2: An example illustration of the four high-level (transition) models
considered in our study, governing the temporal variation of a hypothetical
low-level model parameter.

### High-Level Models

We formulate and compare four different high-level transition models, denoted
as $\mathbb{T}_{1}$, …, $\mathbb{T}_{4}$, which govern the trial-by-trial
changes in local DDM parameters $\theta_{1:T}$. These transition models vary
in terms of their flexibility in allowing changes to the low-level parameters
and their underlying complexity, including the number of high-level parameters
involved (see Figure 2 for an exemplar trajectory generated by each transition
model). To ensure that the low-level parameters remain within plausible
ranges, we impose both lower and upper bounds on their trajectories.
Specifically, we set the upper bounds for the parameters $v$, $a$, and $\tau$
to $8$, $6$, and $4$, respectively. Additionally, since negative parameter
values are not meaningful for our DDM specification, we set the lower bounds
for all parameters to $0$. For all transition models, we assume independence
between the trajectories of the local DDM parameters.

#### Random Walk

The first transition model ($\mathbb{T}_{1}$) convolves the low-level model’s
parameters with a Gaussian distribution, resulting in a gradual change that
follows a random walk:

$\displaystyle\mathbb{T}_{1}(\theta_{k,t}\,|\,\theta_{k,t-1},\sigma_{k})=\mathcal{N}(\theta_{k,t}\,|\,\theta_{k,t-1},\sigma_{k})$
(6)

According to this transition model, the current value of each parameter
$\theta_{k,t}$ is only influenced by its previous value $\theta_{k,t-1}$,
generating more or less auto-correlated and gradual changes.

#### Mixture Random Walk

The second transition model ($\mathbb{T}_{2}$) corresponds to a mixture
distribution between a random walk (cf. Equation 6) and uniformly distributed
shifts:

$\displaystyle\mathbb{T}_{2}(\theta_{k,t}\,|\,\theta_{k,t-1},\rho_{k},\sigma_{k},a_{k},b_{k})$
$\displaystyle=\rho_{k}\,\mathcal{N}(\theta_{k,t}\,|\,\theta_{k,t-1},\sigma_{k})+(1-\rho_{k})\,\mathcal{U}(a_{k},b_{k})$
(7)

where $\rho$ indicates the probability of the type of change (gradual change
or shift) as a mixing coefficient for the two states. The upper and lower
bounds of the uniform distribution, denoted as $a$ and $b$, are set to cover
plausible parameter ranges and are not treated as free parameters themselves.

#### Lévy-Flight

The Lévy flight transition model ($\mathbb{T}_{3}$) is similar to the Gaussian
random walk. However, instead of assuming normally distributed noise, it
assumes an alpha-stable transition for each component of $\theta$:

$\displaystyle\mathbb{T}_{3}(\theta_{k,t}\,|\,\theta_{k,t-1},\sigma_{k},\alpha_{k})=\text{Alpha-
Stable}(\theta_{k,t}\,|\,\theta_{k,t-1},\sigma_{k},\beta=0,\alpha_{k})$ (8)

where $0<\alpha\leq 2$ governs the heaviness of the noise distribution’s
tails. If $\alpha_{k}=2$ then the distribution is equivalent to a Gaussian
distribution. Notably, as the value of $\alpha$ decreases, the distribution’s
tails get heavier, allowing for larger shifts in the parameter values. When
simulating from the Lévy flight transition model, we use a scale of
$\sigma_{k}/\sqrt{2}$, such that the corresponding Gaussian distribution for
$\alpha_{k}=2$ has a standard deviation of $\sigma_{k}$.

#### Regime Switching

Finally, the regime switching transition model ($\mathbb{T}_{4}$) is a simpler
version of the mixture random walk. The parameter’s trajectory adheres to one
of two possibilities: it either maintains its previous value or undergoes a
uniform shift:

$\displaystyle\mathbb{T}_{4}(\theta_{k,t}\,|\,\theta_{k,t-1},\rho_{k},a_{k},b_{k})=\rho_{k}\,\delta(\theta_{k,t}-\theta_{k,t-1})+(1-\rho_{k})\,\mathcal{U}(a_{k},b_{k}),$
(9)

where $\delta(\cdot)$ is the Dirac delta distribution indicating that the
parameter either does not change at all with probability $\rho$ or undergoes a
sudden change with probability $1-\rho$.

Strictly speaking, some of the above transition models can effectively be
transformed into others by employing specific high-level parameter
configurations. For instance, the mixture random walk with $\sigma=0$ reduces
to the regime switching transition function. Conversely, when $\rho=0$ it
reduces to a simple Gaussian random walk. Also, the Lévy flight transition
model with $\alpha=2$ turns into a random walk transition function. The
mixture random walk and the Lévy flight transition function both have two
high-level parameters and can thus be regarded as more complex and more
flexible than the other two transition models, which only have a single high-
level parameter. Notably, the random walk transition model is the only one
that cannot generate relatively large sudden shifts in parameter values.

### Model Comparison Setup

One of the major aims of this study is to compare four NSDDMs sharing the same
low-level diffusion model but differing in their assumptions about the type of
stochastic variation of the drift rate ($v$) and threshold ($a$) parameters.
All four NSDDMs employ the same Gaussian random walk model $\mathbb{T}_{1}$
for the non-decision time parameter ($\tau$). We base this decision on
previous research [11] and the rationale of our experimental manipulations,
which should not imply sudden large shifts in the $\tau$ parameter. For
$\mathcal{M}_{1}$, the drift rate and threshold parameter also follow a
Gaussian random walk, resulting in three high-level parameters,
$\eta=(\sigma_{v},\sigma_{a},\sigma_{\tau})$. In $\mathcal{M}_{2}$, both $v$
and $a$ follow a mixture between a Gaussian random walk and uniform shifts
($\mathbb{T}_{2}$), which results in a total of five high-level parameters,
$\eta=(\sigma_{v},\sigma_{a},\sigma_{\tau},\rho_{v},\rho_{a})$. In contrast,
$\mathcal{M}_{3}$ introduces a trajectory for the drift rate and threshold
parameters characterized by a Lévy flight ($\mathbb{T}_{3}$), which has five
free high-level parameters,
$\eta=(\sigma_{v},\sigma_{a},\sigma_{\tau},\alpha_{v}$, $\alpha_{a})$. Lastly,
for $\mathcal{M}_{4}$, the two parameters $v$ and $a$ either remain the same
as in the previous time point or shift uniformly ($\mathbb{T}_{4}$). This
model has a total of three high-level parameters,
$\eta=(\sigma_{\tau},\rho_{v},\rho_{a})$. A listing of the weakly informative
prior distributions assigned to the model parameters can be found in the
Appendix.

Figure 3: A conceptual illustration of our amortized Bayesian inference
training setup. A Parameter estimation A recurrent summary network processes
the synthetic time series $x_{1:T}$ and learns maximally informative temporal
summary statistics $\tilde{x}_{1:T}$. An inference network (i.e., normalizing
flow) learns to approximate the joint posterior distribution of time-varying
low-level parameters $\theta_{1:T}$ and static high-level parameters $\eta$
given the learned summaries. B Model comparison A transformer summary network
consumes time series simulated from competing models and learns maximally
informative summary vectors $\tilde{x}$. An inference network (i.e., a
probabilistic classifier) learns to approximate posterior model probabilities
(PMPs) given the summary vectors. Once trained, the networks can be
efficiently validated using principled Bayesian methods and applied to the
observed data.

### Amortized Bayesian Inference

Amortized Bayesian inference (ABI) is a flexible framework for estimating,
comparing, and validating complex models through simulation-based training of
specialized neural networks [40]. ABI consists of (i) a training phase where
the networks learn a surrogate distribution, and (ii) an inference phase where
the networks infer the target quantities (e.g., model parameters or model
posterior probabilities) in real-time for any new data set supported by the
model(s). The neural networks are trained purely on simulations from the
generative model and do not require an explicit likelihood or numerical
integration. Thus, ABI re-casts expensive Bayesian inference into a neural
network prediction task, such that sampling from the target posterior and
model refits happen almost instantaneously.

#### Amortized Parameter Estimation

Our deep learning approach for jointly estimating time-varying and static
parameters follows [11], who extend ideas from ABI with static parameters [47,
49] to non-stationary Bayesian models. Accordingly, our goal is not only to
infer the trajectories of all three model parameters
$\\{\theta_{t}\\}^{T}_{t=1}$, but also to estimate the posterior distribution
for the static high-level parameters $\eta$ of the transition model. Thus, we
are interested in recovering the full joint posterior
$p(\theta_{1:T},\eta\,|\,x_{1:T})$ from the observed time series
$\\{x_{t}\\}_{t=1}^{T}$:

$p(\theta_{1:T},\eta\,|\,x_{1:T})\propto
p(\eta,\theta_{0})\,p(x_{1}\,|\,\theta_{1})\,\prod_{t=2}^{T}p(x_{t}\,|\,\theta_{t},x_{1:t-1})\,\prod_{t=1}^{T}\mathbb{T}(\theta_{t}\,|\,\eta,\theta_{0:t-1})$
(10)

where $p(\eta,\theta_{0})$ is the joint prior over high-level parameters and
initial low-level parameter values. The joint prior typically factorizes as
$p(\eta,\theta_{0})=p(\eta)p(\theta_{0})$, assuming that $\eta$ and
$\theta_{0}$ are independent in the absence of any information. Even though
our SBI method is applicable to any model of the general form in Eq. 10, our
low-level (Low-Level Model) and high-level (High-Level Models) specifications
lead to a simplified formulation

$p(\theta_{1:T},\eta\,|\,x_{1:T})\propto
p(\eta,\theta_{0})\,\prod_{t=1}^{T}p(x_{t}\,|\,\theta_{t})\,\prod_{t=1}^{T}\mathbb{T}(\theta_{t}\,|\,\eta,\theta_{t-1}).$
(11)

The simplified formulation follows from the fact that our transition models
share the Markov property and the DDM likelihood depends on time only through
the current parameter $\theta_{t}$ in the latent trajectory $\theta_{1:T}$.

Following the typical ABI offline training setting (see Figure 3A for a
conceptual illustration), we generate a data set of simulated data sets,
$\smash{\mathcal{D}=\\{\eta^{(b)},\theta^{(b)}_{1:T},x^{(b)}_{1:T}\\}_{b=1}^{B}}$,
and use the simulated data to train a specialized neural network,
$F_{\psi}(\theta_{1:T},\eta;\,x_{1:T})$, which approximates the full joint
posterior [50, i.e., a normalizing flow, see]. In particular, we minimize the
following loss in expectation over the full non-stationary generative model
(i.e., the right hand-side of Eq. 10)

$\mathcal{L}(\psi)=\mathbb{E}_{(\eta,\theta_{1:T},x_{1:T})\sim\mathcal{D}}\left[-\log
q_{\psi}(\theta_{1:T},\eta\,|\,x_{1:T})\right],$ (12)

where we approximate the expectation over
$p(\theta_{0})\,p(\eta,\theta_{1:T},x_{1:T})$ via our training set
$\mathcal{D}$ and regularize against overfitting with standard techniques,
such as dropout and weight decay. It is also possible to run the simulator(s)
indefinitely and perform online training using on-the-fly simulation [47]. In
fact, this approach should be preferred for fast simulators, as it makes
overfitting hardly possible. Thus, online learning is the approach we pursue
for estimating the parameters of our NSDDMs.

In the context of dynamic Bayesian models, we have many choices on how to
factorize the joint posterior [51]. The two most common choices are to
approximate the filtering distribution or the smoothing distribution [31]. The
filtering distribution corresponds to an online analysis, where the low-level
parameters $\theta_{t}$ at time step $t$ are only informed by past data
points. Differently, the smoothing distribution conditions the posterior of
$\theta_{t}$ on all past and future data points, and provides potentially
sharper estimates. Thus, in this study, we exclusively target the approximate
smoothing distribution due to its superior parameter recoverability in an
offline analysis.111Note, that [11] focused exclusively on the filtering
distribution in their benchmarking experiments. In practice, we employ
unidirectional or bidirectional long-short term memory (LSTM) networks [52]
with many-to-many input-output relationships as a backbone for approximating
the filtering or the smoothing distribution, respectively. We then train four
separate neural approximators, such that each network becomes an “expert” in
inferring the smoothing distribution of the corresponding NSDDM. The Appendix
contains more details on the neural network settings and training
hyperparameters.

#### Amortized Model Comparison

To conduct a comparative analysis of the four NSDDMs, we focus on Bayes
factors (BFs) and posterior model probabilities (PMPs), which can be
classified as prior predictive methods embodying Occam’s razor [53, 54]. The
efficacy of these measures has been demonstrated in a wide range of
psychological modeling studies [55]. Nevertheless, an ongoing debate surrounds
the preference between the two [56, 57]. Since BFs and posterior odds (i.e.,
ratios between PMPs) are equivalent when all models are assumed to be equally
likely a priori, we estimate and analyse both quantities in our study.

Following the common Bayesian terminology [54], we can refer to the four
competing models through an index set
$\mathcal{M}=\\{\mathcal{M}_{1},\mathcal{M}_{2},\mathcal{M}_{3},\mathcal{M}_{4}\\}$.
The aim of prior predictive Bayesian model comparison is to find the simplest
most plausible model within $\mathcal{M}$. To this end, we can compute PMPs
for each of the competing models

$p(\mathcal{M}_{j}\,|\,x_{1:T})=\frac{p(x_{1:T}\,|\,\mathcal{M}_{j})\,p(\mathcal{M}_{j})}{\mathbb{E}_{p(\mathcal{M})}\left[p(x_{1:T}\,|\,\mathcal{M})\right]},$
(13)

where $p(\mathcal{M})$ refers to the prior distribution over the discrete
model space. The marginal likelihood $p(x_{1:T}\,|\,\mathcal{M}_{j})$ plays a
crucial role in Equation 13, and can be expressed by integrating out all
parameters of the joint model,

$p(x_{1:T}\,|\,\mathcal{M}_{j})=\int
p(\eta,\theta_{0})\,\prod_{t=1}^{T}p(x_{t}\,|\,\theta_{t},\mathcal{M}_{j})\,\prod_{t=1}^{T}\mathbb{T}_{j}(\theta_{t}\,|\,\eta,\theta_{t-1})\,d\eta\,d\theta_{0},\dots,d\theta_{T}.$
(14)

Importantly, since the marginal likelihood averages the likelihood over the
joint prior, it automatically incorporates a probabilistic Occam’s razor,
favoring models with constrained prior predictive flexibility. When comparing
a pair of competing models, $\mathcal{M}_{j}$ and $\mathcal{M}_{i}$, we can
compute the ratio between their respective marginal likelihood,

$\text{BF}_{ji}=\frac{p(x_{1:T}\,|\,\mathcal{M}_{j})}{p(x_{1:T}\,|\,\mathcal{M}_{i})}.$
(15)

This ratio is referred to as the Bayes factor (BF). Consequently, a
$\text{BF}_{ji}>1$ signifies a relative preference for model $j$ over model
$i$ based on the given data $x_{1:T}$ [53].

Unfortunately, the marginal likelihood is notoriously hard to approximate [58]
and even doubly intractable for mechanistic models with unknown or
unnormalized likelihoods. To circumvent this intractability, we follow the
neural method of [59, 60] which enables amortized Bayesian model comparison
for arbitrary computational models (see Figure 3B for a graphical
illustration). This method involves the simultaneous training of two neural
networks with different roles: a summary network and an inference network. The
summary network learns maximally informative summary statistics from the raw
data (e.g., behavioral time series). The inference network approximates the
PMPs for the candidate models, $q_{\phi}(\mathcal{M}\,|\,x_{1:T})$ given the
outputs of the summary network. Here, we subsume all trainable network
parameters under $\phi$ and refer to the composition of the two networks as an
evidential network.

The training data for the evidential network consists of all simulations from
the candidate models together with the corresponding model index,
$\smash{\mathcal{D}(\mathcal{M})=\\{x_{1:T}^{(b)},\mathcal{M}_{j}^{(b)}\\}_{b=1}^{B^{\prime}}}$,
where $B^{\prime}$ denotes the total number of simulations from all models.
Together, the two networks minimize the standard cross-entropy loss,

$\mathcal{L}(\phi)=\mathbb{E}_{(\mathcal{M}_{j},x_{1:T})\sim\mathcal{D}(\mathcal{M})}\Big{[}-\sum_{j=1}^{J}\mathbb{I}_{\mathcal{M}_{j}}\log
q_{\phi}(\mathcal{M}_{j}\,|\,x_{1:T})\Big{]},$ (16)

and we approximate the expectation over $p(\eta,\theta_{1:T},x_{1:T})$ by our
training set $\mathcal{D}(\mathcal{M})$, and $\mathbb{I}_{\mathcal{M}_{j}}$
denotes an indicator function (i.e., one-hot encoding) for the true model
index. In principle, we could use online learning for amortized model
comparison as well, but we found offline training to yield sufficiently
accurate results.

More recently, [61] demonstrated the importance of gauging the sensitivity of
amortized neural approximators, especially in the context of model comparison.
The authors suggest to train an ensemble of multiple evidential networks,
instead of relying on a single network. Accordingly, we can measure the (lack
of) agreement between ensemble members and obtain a hint at the robustness of
the approximate PMPs. Here, we trained an ensemble of ten evidential networks
and computed the mean and standard deviation of the estimated PMPs across all
ten networks. For more details regarding the neural network architecture and
training settings, we refer the reader to the Appendix.

#### Code Availability

Complete code for reproducing the results reported in this manuscript is
available in the project’s GitHub repository https://github.com/bayesflow-
org/Non-Stationary-DDM-Validation.

## Results

### Model Comparison

Figure 4: In silico model comparison and sensitivity results. A Calibration
curves of all four NSDDMs aggregated across the neural approximator ensemble.
Additionally, the expected calibration error ($\widehat{{\mathrm{{ECE}}}}$) is
annotated within each subfigure. The gray histograms depict the relative
frequencies of the predicted model probabilities. B Confusion matrix between
true data generating model and predicted model. The proportion values were
averaged across the ten neural approximator within the ensemble.

As a first step, we assess the closed-world (i.e., in silico) performance of
our model comparison method in terms of computational faithfulness and
accuracy of model recovery. To assess the former, we perform simulation-based
calibration [41, 42, SBC;] based on $10\,000$ synthetic data sets per model.
Figure 4A shows the calibration curves for each NSDDM averaged across the ten
evidential networks in our deep ensemble. We observe excellent calibration
with very minimal expected calibration errors ($\widehat{{\mathrm{{ECE}}}}$)
across all models. Thus, we conclude that the approximate posterior
probabilities are well-calibrated in the closed-world setting.

Next, we assess the accuracy of our model comparison networks in terms of
their ability to correctly identify the ground-truth data-generating model. To
this end, we apply the deep ensemble to the $40\,000$ synthetic data sets we
have already simulated for assessing calibration. In Figure 4B, we present the
resulting confusion matrix, which illustrates the agreement between true and
predicted models averaged across the ten approximators. Among the four models,
the random walk DDM is the only one that rarely gets confused with the other
models. A possible explanation is that it is the only transition model not
capable of generating sudden shifts in parameter values. The remaining models
are susceptible to more frequent misclassifications. For example, the mixture
random walk DDM is correctly identified only $54\%$ of the time, and it is
often confused with the regime switching model, occurring $43\%$ of the time.
Notably, the Lévy flight DDM is prone to mimicry with the random walk DDM (on
average $30\%$ of the time).

It is essential to emphasize that these results do not imply a deficiency in
your model comparison method, but rather underscore the fact that certain
pairs of models, such as the mixture random walk and the regime switching DDM,
can generate remarkably similar data patterns. For instance, a significant
portion of the prior distribution’s mass for the $\alpha$ parameter of the
Lévy flight transition model centers around $2$. If $\alpha\approx 2$, then
the Lévy alpha-stable distribution closely resembles a Gaussian distribution,
with equality in the case of $\alpha=2$. Consequently, simulating the Lévy
flight DDM would often yield data patterns that could have just as plausibly
originated from the simpler random walk DDM.

Similarly, a substantial portion of the prior mass for the $\sigma$ priors of
the mixture random walk transition model clusters around $0$, which
subsequently transforms it into a regime switching transition model, resulting
in large overlap in synthetic data sets. Interestingly, the mixture random
walk and the Lévy flight DDM are seldom confused, even though both models can
produce subtle local changes and large sudden shifts. This implies that these
two transition models generate qualitatively similar but quantitatively easy
to distinguish parameter trajectories. In summary, the observation of
occasional model confusion is not a limitation of our method; rather, it
underscores our method’s effectiveness in discerning when two models generate
highly similar data, making them less straightforward to differentiate from
each other. Moreover, the amortization property of our method enables us to
easily conduct such simulation studies prior to analyzing real data –
estimating $40\,000$ posterior model probabilities would have been infeasible
for any other method.

Figure 5: Empirical model comparison results. A Aggregate posterior model
probabilities (PMP) across the ensemble and the $14$ individual participants.
Points depict the mean, stars the median, and the error bars indicate the
$75\%$ credibility interval (CI). B Heatmap of average $\text{log}_{10}$ Bayes
factors (BF). Both metrics agree on favoring the Lévy flight DDM over the
other models.

After successfully validating our model comparison method, we apply the deep
ensemble to the empirical data of the $14$ participants. Each approximator in
the ensemble was used to infer posterior model probabilities (PMP) for each
model, considering each individual’s data separately. Subsequently, we
calculated the mean (points), median (stars), and $75\%$ credibility interval
(CI) for the PMPs for all approximators the $14$ individuals (Figure 5A). The
analysis reveals that the Lévy flight DDM is the most plausible model with an
average PMP of approximately $60\%$. It was the most plausible model for $9$
out of the $14$ participants. In contrast, the mixture random walk model
collects an average PMP of less than $30\%$. Nevertheless, it was estimated to
be the most plausible model for $5$ participants. The random walk DDM and
regime switching DDM were consistently less plausible than the other models
and did not emerge as superior for any of the participants.

In addition to PMPs, we computed $\text{log}_{10}$ Bayes factors (BF). Figure
5B depicts a heatmap of BFs for all one-to-one comparisons between our four
NSDDMs, averaged across the participants and the evidential networks of the
ensemble. Following [53], an absolute value of $\text{log}_{10}(\text{BF})>2$
indicates decisive evidence, absolute values between $1$ to $2$ signify
strong, and between $0.5$ to $1$ substantial evidence. An absolute value of
$\text{log}_{10}(\text{BF})<0.5$ is labeled as not worth more than a bare
mention. The BF patterns in Figure 5B align with the PMP findings, implying
strong evidence for the Lévy flight DDM over the random walk DDM and
substantial evidence over the other NSDDMs. Also, both the mixture random walk
and the regime switching DDM have substantial evidence over the random walk
model. Interestingly, there is little evidence favoring the mixture random
walk DDM over the regime-switching model, suggesting comparable performance.

These findings offer two substantive insights. First, the ability of
transition models to generate sudden shifts in parameters seems essential, as
seen in the random walk DDM’s lower plausibility. Moreover, the regime
switching DDM, allowing for occasional shifts, but neglecting small gradual
changes, performed less effectively than the more complex models. This result
underscores the importance of accommodating both gradual as well as sharp
changes in model parameters for achieving optimal fit. Consequently, the more
complex NSDDMs, particularly the mixture random walk DDM and Lévy flight DDM,
emerged as more plausible than their simpler counterparts, despite the
implicit penalty for prior complexity imposed by Bayesian model comparison.

Figure 6: Aggregated results from all models fitted to the empirical data. The
top row illustrates posterior re-simulations as a measure of the model’s
generative performance and absolute goodness-of-fit to the data. The bottom
row depicts parameter estimates of the drift rate and the threshold parameter
from the non-stationary diffusion decision models (NSDDM). A Empirical and re-
simulated RTs for each difficulty level and both conditions. B Empirical and
re-simulated proportions of correct choices (accuracy) for each difficulty
level and both conditions separately. C Posterior estimates of the drift rate
parameter for each difficulty level and both conditions separately. D
Posterior estimates of the threshold parameter for each difficulty level and
both conditions separately. Points indicate medians and the error bars
represent the median absolute deviations (MAD) across individuals and re-
simulations.

### Posterior Re-simulation

Subsequently, we fit all four variants of the NSDDM to each of the $14$ data
sets, evaluating the absolute goodness-of-fit of each model. To achieve this,
we conducted $500$ re-simulations with randomly sampled posterior parameter
trajectories for each individual data set. In Figure 6A, we present the median
and median absolute deviation (MAD) of response times (RT) across all
individuals and re-simulations. We provide these aggregates for each NSDDM,
categorized by task difficulty level and the two experimental conditions.
Notably, an initial observation reveals that the experimental manipulations
were effective on average: empirical median RTs increased with task
difficulty, and individuals tended to respond faster in the speed condition
compared to the accuracy condition. Remarkably, all four variants of the NSDDM
demonstrated an outstanding fit to these empirical data patterns. Solely, RTs
in the accuracy condition with the highest task difficult level consistently
are underestimated by all NSDDM variants.

The empirical and re-simulated proportion of correct choices (accuracy) are
aggregated and presented in the same way as the RTs (see Figure 6B). Again,
the empirical data mirror the anticipated patterns resulting from our
experimental manipulations. As expected, accuracy diminishes with increasing
task difficulty. Individuals are generally less accurate in the speed
condition compared to the accuracy condition. Although NSDDMs successfully
reproduce the general patterns in the choice data, we observe notably worse
re-simulation compared to that of the RTs data. In both accuracy and speed
conditions, re-simulated accuracies exhibit a less pronounced decline as a
function of difficulty than observed in the empirical data. Further, the
difference in accuracy between the two experimental conditions is less
pronounced in the re-simulated data compared to the behavioral data. Notably,
the random walk DDM underperforms relative to the other three NSDDMs in this
analysis.

It is important to highlight that, unlike conventional approaches, the models
did not receive any information regarding the specific experimental context an
individual faced at any given moment. From these analyses, we conclude that
all NSDDM implementations successfully capture the general patterns in the
empirical RT data. Individual participant analyses, detailed in the Appendix,
affirm the same conclusions.

In addition to the analysis of the absolute model fit on the aggregate level,
we evaluated the fit across the RT time series. For every participant, we
generated $250$ posterior re-simulations for the first $700$ trials with the
corresponding best fitting NSDDM. The remaining $68$ data points were left out
for predictive analysis. Employing a one-step-ahead prediction approach, we
iteratively forecasted the subsequent data point, followed by a re-fitting of
the model in each step.

Figure 7 illustrates the empirical and re-simulated RT time series for two
exemplary participants. Results for the remaining $12$ participants can be
seen in the Appendix. The colored lines depict the median and the shaded bands
represent to $90\%$ credibility intervals (CI) across the $250$ re-
simulations. Both the empirical data (solid black lines) and the re-
simulated/predicted RTs were smoothed using a simple moving average (SMA) with
a period of $5$. Yellow shaded regions highlight trials where speed was
emphasised over accuracy, whereas blank white areas denote instances where the
opposite emphasis was applied. Overall, RTs were slower and more variable in
the accuracy condition. Notably, the NSDDM not only closely replicated the
empirical time series but also effectively predicted future data points. This
suggests that the model does not overfit the data.

Figure 7: Model fit to response time (RT) time series. The empirical RT time
series of two exemplar individuals are shown in black. From trial $1$ to
$700$, the posterior re-simulations (aka retrodictive checks) using the best
fitting non-stationary diffusion decision model (NSDDM) for the specific
individual are shown in blue and red, respectively. In this instance, the left
column showcases results from a Lévy flight DDM, while the right column
displays parameter trajectories from a mixture random walk DDM. For the
remaining trials, one-step-ahead posterior predictions from the NSDDMs are
depicted in cyan and orange, respectively. Solid lines correspond to the
median and shaded bands to $90\%$ credibility intervals (CI). The empirical,
re-simulated, and predicted RT time series were smoothed via a simple moving
average (SMA) with a period of $5$. The yellow shaded regions indicate trials
where speed was emphasised over accuracy, while blank white areas denote
instances where the opposite emphasis was applied.

### Parameter Estimates

At the heart of the current validation study are the inferred parameters,
prompting a crucial question: Do these parameter dynamics align with the
sequence of experimental manipulations? We address this question by examining
both the time-averaged and time-varying estimates.

#### Aggregate Analysis

We initially examine the parameter estimates averaged across individuals for
each difficulty level and condition separately. This provides a comprehensive
overview of average effects on model parameters in different experimental
contexts, at first, without delving into the temporal aspect. The bottom panel
of Figure 6 illustrates the posterior medians and MADs collapsed onto the
different experimental contexts for the drift rate (Figure 6C) and threshold
parameter (Figure 6D).

Analyzing the aggregated drift rate estimates reveal an anticipated pattern.
On average, the drift rate decreases as task difficulty increases, observed in
both the accuracy and speed conditions. Additionally, slightly higher overall
values are estimated in the speed condition compared to the accuracy
condition. While all four NSDDMs yield fairly similar parameter values, the
distinctions in average parameter values between difficulty levels are less
pronounced when estimated with the random walk DDM.

With the second experimental manipulation - namely, the instruction to
emphasize speed or accuracy - we aimed to manipulate the participants’
decision caution, which is assumed to be captured by the threshold parameter.
Examining the aggregated estimates of the threshold parameter in Figure 6D, we
observe generally increased values in the accuracy condition compared to the
speed condition. Interestingly, in the accuracy condition, the threshold
parameter also slightly increases with growing task difficulty — a pattern not
observed in the speed condition. A comparison between the estimates of the
four NSDDMs reveals that the mixture random walk DDM and the Lévy flight DDM
yield higher threshold estimates in the accuracy condition compared to the
other two NSDDMs. Conversely, all four NSDDMs seem to converge in their
threshold parameter estimates in the speed condition.

Figure 8: Estimated parameter trajectories of two exemplar individuals
corresponding to the respective best-fitting non-stationary diffusion decision
model (NSDDM). In this instance, the left column showcases results from a Lévy
flight DDM, while the right column displays parameter trajectories from a
mixture random walk DDM. Each low-level parameter (drift rate, threshold, and
non-decision time) is displayed on a separate row. The solid lines are color-
coded (blue for the Lévy flight DDM and red for the mixture random walk DDM)
to represent the posterior medians, while the shaded regions mark the median
absolute deviation (MAD). The yellow shaded regions indicate trials where
speed was emphasised over accuracy, while blank white areas denote instances
where the opposite emphasis was applied. The sequences of task difficulty
levels are depicted with black lines and overlaid with the drift rate in the
top panels.

#### Parameter Trajectories

For a more fine-grained analysis, particularly considering temporal aspects,
we present the complete inferred parameter trajectories of the three low-level
parameters of a NSDDM for two exemplary individuals in Figure 8. The Appendix
contains the inferred parameter trajectories of the remaining $12$
participants. Each participant’s trajectory is depicted with the posterior
median (solid lines) and the median absolute deviation (MAD, shaded bands)
across all $768$ experimental trials, estimated with the model with the
highest posterior model probability for that specific individual. The
trajectory of participant $11$ corresponds to a Lévy flight DDM, whereas the
trajectory of participant $6$ comes from a mixture random walk DDM. Shaded
blocks along the timeline denote the experimental condition at a given trial,
with yellow indicating an emphasis on speed.

The top panel illustrates the estimated trajectories of the drift rate
parameter alongside the sequences of task difficulty levels (depicted by black
line). Here, $0$ corresponds to the most difficult level, while $6$ represents
the easiest. It is important to note that the absolute values of the
difficulty conditions hold no intrinsic meaning. As observed, the drift rates
of both participants align with the overarching trend of the difficulty
condition sequence. They decrease when the difficulty is high and increase as
the task becomes easier.

Regarding the trajectory of the threshold parameter (middle panel), we
anticipated that a shift from an accuracy instruction to a speed instruction
would lead to a decrease in the threshold parameter, and vice versa. This
hypothesized pattern is clearly evident when examining the estimated threshold
parameter trajectories of the two participants in the middle panel of Figure
8. For instance, the threshold parameter estimated for participant $11$
oscillates around an approximate value of $1$ in the speed condition.
Moreover, it consistently rises whenever a switch in the accuracy condition
takes place. Intriguingly, the parameter’s value during accuracy emphasis is
not as uniform compared to the speed condition. In some blocks, it fluctuates
around $2$, while in others it hovers around $1.5$ or even lower. Similarly,
participant $6$ displays pronounced shifts in the threshold parameter when a
change in the condition occurs, with these shifts being more pronounced in the
first half of the experiment and diminishing in the second half.

Finally, the bottom panel of Figure 8 illustrates the trajectory of the non-
decision time parameter. Although our experimental manipulations did not
systematically target the dynamics of this parameter, it is sometimes assumed
that the manipulation of speed and accuracy instructions may also influence it
[34, 36]. While both individuals exhibit some fluctuations in $\tau$, no
systematic differences between the two conditions are apparent.

Upon reviewing the parameter trajectories of the remaining participants in the
Appendix, similar patterns emerge. In summary, both the inferred means and
trajectories of the drift rate and threshold parameters align with the
sequence of experimental manipulations, as predicted by our design. Moreover,
our NSDDMs were able to estimate these trajectories directly from the
behavioral data, getting no explicit information whatsoever about the
experimental context. Thus, our validation study suggests that NSDDMs are able
to detect genuine changes in cognitive constructs.

## Discussion

Psychology and cognitive science are witnessing a growing interest in
incorporating dynamic aspects into mechanistic models that seek to formalize
and explain cognitive processes. In a previous study, we explored a method to
estimate plausible trajectories of cognitive process model parameters directly
from behavioral data [11]. Nevertheless, an empirical validation of this
modeling approach was lacking. Thus, the current study sought to bridge this
gap by experimentally examining the validity of the inferred diffusion
decision model (DDM) parameter dynamics.

### Experimental validation

The present study posed to the following core question: Can non-stationary
DDMs (NSDDM) effectively detect experimentally induced changes in cognitive
constructs from behavioral data alone? If this holds true, our findings can
provide the first substantial evidence for the validity of the superstatistics
framework as applied to cognitive models. Notably, our results demonstrated
that the NSDDMs indeed reliably identified the sequence of two experimental
manipulations, despite the absence of any contextual information. Moreover,
posterior re-simulation revealed an outstanding fit to the response data, both
on an aggregate level as well as on the level of the raw time series. This
performance stands as compelling evidence supporting the validity of NSDDMs.

The trajectory of the drift rate parameter for all individuals closely
mirrored the sequence of the task difficulty manipulation. Specifically, the
drift rate parameter decreased when task difficulty increased, and conversely,
increased as task difficulty decreased. This not only confirms the anticipated
impact of the manipulation, but also highlights the NSDDMs’ ability to discern
these variations directly from the behavioral data, agnostic to additional
contextual information.

Interestingly, drift rates increased throughout the experiment, although this
was not the case for trials with the highest task difficulty. This observation
suggests a practice effect among participants, where task performance
generally improved with experience, except under the most challenging
condition. Practice effects are a widely recognized phenomenon in various
decision-making and memory paradigms [62, 63, 64, 65, 66]. In fact, practice
effects have been studied with various dynamic cognitive modeling approaches
[21, 16, 20, 22]. A notable contribution to this field comes from [16], who
conducted a comprehensive re-analysis of three datasets derived from widely
cited articles. Their study compared three dynamic models: (i) a smooth
polynomial trend, (ii) a non-smooth autoregressive process, and (iii) a regime
switching model instantiated by a hidden Markov model (HMM) with two different
states.

In their study, [16] employed a low-level model similar to the DDM, namely the
linear ballistic accumulator model [67, LBA;]. However, their transition
models, specifically the polynomial trend and the autoregressive process,
differed in that they allowed LBA parameters to change only from block-to-
block, neglecting trial-to-trial parameter fluctuations (except for the HMM).
Their findings indicated that the HMM outperformed the other two dynamic model
instantiations. This superiority can possibly be attributed to the model’s
capacity to flexibly change parameters from trial-to-trial, in contrast to
changes occurring only from block-to-block. Even though the trial-by-trial
specification of the HMM captures the microstructure of the decision-making
process, it is still less flexible than the models we examined in the current
study. HMMs assume a pre-defined number of possible states, whereas this is
not the case with the implementation of our regime switching model. The
advantage of not fixing the number of distinct states beforehand is
particularly evident when the exact latent quantity is unknown prior to
investigation. Moreover, results from our model comparison clearly favored
transition models that account for both, gradual changes as well as sudden
shifts. This suggests that regime-switching models may fall short in certain
fields of application. Nevertheless, both models have their merits, and the
choice between them should be guided by the specific research question at hand
and formal model comparison.

As our study focused on experimentally validating parameter trajectories
estimated with NSDDMs, we deliberately refrained from further analysing
practice effects. However, we suggest that our flexible framework could be a
promising alternative for investigating practice effects. Unlike pure regime
switching models, it has the capacity to reveal a mixture of practice-related
changes, ranging from abrupt shifts to gradual changes. When exploring
substantive research questions, such as practice effects, with
superstatistical models, it is imperative to depart from the approach taken in
the current study. That is, one should always incorporate contextual
information from the experimental setting when estimating parameter
trajectories. Here the question arises, how to incorporate this information?
In a previous study, we simply assumed separate low-level parameters for each
experimental condition [11]. This approach is particularly appropriate when
conditions randomly change from trial to trial. However, future research could
explore alternative ways of including experimental context information with
the goal to further inform the parameters.

Concerning the second experimental manipulation, that is, the emphasis on
speed or accuracy, their effect on the threshold parameter is more diverse
across individuals. While a majority of participants demonstrated shifts in
the threshold parameter in response to instructional changes, the consistency
and magnitude of these changes varied significantly among individuals. Some
participants exhibited only a few adjustments in the threshold parameters,
seemingly overlooking the change in instruction on certain occasions. In
contrast, others consistently heightened their threshold parameter during
accuracy-focused tasks, followed by a subsequent decrease when transitioning
to speed-oriented conditions. Meanwhile, some participants displayed rather
unsystematic changes in decision caution, suggesting that the participants of
reacted differently to the speed-accuracy manipulation.

[15] introduced a dynamic LBA incorporating a hidden Markov transition model
with two states, akin to the model proposed by [16]. Their focus centered on
scrutinizing the speed-accuracy trade-off, exploring the hypothesis that
individuals dynamically switch between different operating states under
varying instruction conditions. By fitting their model to previously collected
data, they provided evidence that individuals tend to oscillate between two
stable states: a deliberative, stimulus-driven mode emphasizing accuracy and
sacrificing speed, and a guessing mode characterized by random and relatively
faster choices.

However, our approach for estimating parameter trajectories reveals a more
intricate scenario, challenging the assumed binary operational shift. Contrary
to expectations, individuals manifest more than two discernible states. At
times, they exhibit an extreme adaptation to a change in condition, while at
other times, they display little or no reaction to the altered condition. This
complexity underscores the necessity for more flexible transition models, as
employed in our study. Failing to utilize such adaptive models could
potentially obscure the complex unfolding of individuals’ cognition and
behavior over time.

### Model comparison

When implementing non-stationary models, a modeler encounters a myriad of
options, ranging from various transition models to decisions about which
parameter follows which transition model. In this study, we limited our
choices to a small subset of the possibility space. Based on our experimental
manipulations we anticipated that the DDM parameters, particularly the
threshold parameter, would not only undergo gradual changes, but also manifest
more abrupt shifts in response to changing conditions. Consequently, we tested
different implementations accommodating such shifts (mixture random walk, Lévy
flight, regime switching) against a transition model that does not, namely,
the simple Gaussian random walk.

The inferred posterior model probabilities (PMPs) and Bayes factors (BFs)
consistently favored the Lévy flight and occasionally the mixture random walk
transition models. However, in terms of the absolute goodness-of-fit, as
assessed through posterior re-simulations, the performance of all four NSDDMs
showed remarkable similarity. This leads to two notable conclusions. First,
even the models with lower PMPs demonstrated a good fit to the data, likely
owing to the inherent flexibility of the superstatistical framework. Second,
our Bayesian model comparison method could reliably detect the most favorable
model even when the absolute differences were marginal.

### Limitations

Psychological research is usually interested in some group or overall estimate
of parameters. Thus, it would have been informative to compute and inspect
“average” parameter trajectories. Unfortunately, our experiment was designed
in a way that the difficulty and the speed-accuracy instruction manipulation
was randomized across participants. This made it impossible to average the
individual trajectories directly. Instead, we collapsed the estimates by the
different experimental conditions and provided an aggregate view across
individuals. Although this is certainly a limitation of this study, we argue
that the current analysis is sufficient to address our specific research
question.

Moreover, despite using many default settings from the BayesFlow software
[40], the configuration and training of neural approximators for both
parameter inference and model comparison for non-stationary models can still
be a challenge. A basic understanding of deep learning principles and
simulation-based inference is an essential prerequisite. These requirements
may pose obstacles to the adoption of our method, highlighting the necessity
for improved software and tutorials addressing these intricacies.

### Outlook

Going forward, we see the relevance of our superstatistics framework as
twofold. First, superstatistics could become a powerful tool in the
methodological toolkit of the researcher interested in temporal changes in
cognitive constructs. It is a general framework and provides large
flexibility. Thus far, we only used the DDM as a low-level observation model.
However, there are many other cognitive process models that could benefit from
such an framework. For instance, reinforcement learning model parameters, such
as the learning rate or the softmax temperature parameter,likely change over
time [14]. Second, even when the temporal evolution of cognitive parameters is
not a central research question, the adoption of non-stationary models may
bring advantages over their stationary counterparts [11]. Our analysis of
estimated trajectories vividly illustrates discernible changes in parameters.
Assuming stationarity would have led to misleading substantive conclusions.

With great flexibility comes a great plethora of choices. In this study, we
compared different transition models guided by the contrast between gradual
and sudden changes. However, there are more degrees of freedom when
implementing superstatistical models, or Bayesian models in general [44]. [61]
advocates for the crucial role of sensitivity analysis, illustrating a potent
methodology to facilitate informed decisions regarding factors such as the
type and shape of prior distributions, neural network architectures, and other
pivotal elements. We believe that using such an approach in the context of
superstatistics could provide better guidelines for their implementation.

Up to this point, we focused on the estimates of the low-level parameter
trajectories. Yet, it is crucial to note that we also obtain posterior
distributions for the static high-level parameters. These estimates can also
yield valuable insights into individuals’ behavior and cognition. Depending on
the chosen transition model, these estimates can offer indications of the
frequency with which individuals transition between distinct operational
states or the variability inherent in their cognitive constructs. Thus,
analyzing these high-level parameters could constitute a compelling avenue for
future research.

### Conclusion

In conclusion, the experimental validation of non-stationary diffusion
decision models presented in this study represents a significant step forward
in the field of cognitive modeling. Our results provide compelling evidence
that the estimated parameter trajectories genuinely reflect tangible changes
in the targeted psychological constructs. We hope that our validation opens
the door to widespread applications of non-stationary models in future
modeling endeavors, offering a more nuanced understanding of cognitive
processes across varying time scales.

### Acknowledgments

We thank Steffen Ernst for his efforts in programming the experiment and
collecting the data. L.S., M.S., and A.V. were supported by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation; grant number GRK 2277
“Statistical Modeling in Psychology”).

### Author contributions

L.S. and M.S conceived the initial idea and the experimental design. L.S. and
S.T.R. created and applied the models, and wrote the initial draft of the
manuscript. A.V. supervised the project. All authors reviewed and refined the
initial draft of the manuscript and agreed to its current version.

### Data and Code Availability

All models, data, and scripts for reproducing the results of this work are
publicly available in the project’s GitHub repository
https://github.com/bayesflow-org/Non-Stationary-DDM-Validation. The neural
superstatistics method is implemented in the BayesFlow Python library for
amortized Bayesian workflows [40].

## References

  * [1] J.. Galadí, S. Silva Pereira, Y. Sanz Perl, M.. Kringelbach, I. Gayte, H. Laufs, E. Tagliazucchi, J.. Langa and G. Deco  “Capturing the Non-Stationarity of Whole-Brain Dynamics Underlying Human Brain States”  In _NeuroImage_ 244, 2021, pp. 118551  DOI: 10.1016/j.neuroimage.2021.118551
  * [2] Alexandre Melanson, Jorge F. Mejias, James J. Jun, Leonard Maler and André Longtin  “Nonstationary Stochastic Dynamics Underlie Spontaneous Transitions between Active and Inactive Behavioral States”  In _eNeuro_ 4.2 Society for Neuroscience, 2017  DOI: 10.1523/ENEURO.0355-16.2017
  * [3] Guy C. Van Orden, John G. Holden and Michael T. Turvey  “Self-Organization of Cognitive Performance”  In _Journal of Experimental Psychology: General_ 132.3 US: American Psychological Association, 2003, pp. 331–350  DOI: 10.1037/0096-3445.132.3.331
  * [4] Eric-Jan Wagenmakers, Simon Farrell and Roger Ratcliff  “Estimation and Interpretation of $1/f^{\alpha}$ Noise in Human Cognition”  In _Psychonomic Bulletin & Review_ 11.4, 2004, pp. 579–615  DOI: 10.3758/BF03196615
  * [5] Rob R. Sebastian Castro-Alvarez and Jorge N. Tendeiro  “A Time-Varying Dynamic Partial Credit Model to Analyze Polytomous and Multivariate Time Series Data”  In _Multivariate Behavioral Research_ 0.0 Routledge, 2023, pp. 1–20  DOI: 10.1080/00273171.2023.2214787
  * [6] Peter F. Craigmile, Mario Peruggia and Trisha Van Zandt  “Hierarchical Bayes Models for Response Time Data”  In _Psychometrika_ 75.4, 2010, pp. 613–632  DOI: 10.1007/s11336-010-9172-6
  * [7] Daniel H Weissman, KC Roberts, KM Visscher and MG Woldorff  “The neural bases of momentary lapses in attention”  In _Nature neuroscience_ 9.7 Nature Publishing Group US New York, 2006, pp. 971–978
  * [8] Michael Esterman and David Rothlein  “Models of sustained attention”  In _Current opinion in psychology_ 29 Elsevier, 2019, pp. 174–180
  * [9] Roger Ratcliff and Hans P.. Van Dongen  “Diffusion Model for One-Choice Reaction-Time Tasks and the Cognitive Effects of Sleep Deprivation”  In _Proceedings of the National Academy of Sciences_ 108.27 Proceedings of the National Academy of Sciences, 2011, pp. 11285–11290  DOI: 10.1073/pnas.1100483108
  * [10] Matthew M. Walsh, Glenn Gunzelmann and Hans P.. Van Dongen  “Computational Cognitive Modeling of the Temporal Dynamics of Fatigue from Sleep Loss”  In _Psychonomic Bulletin & Review_ 24.6, 2017, pp. 1785–1807  DOI: 10.3758/s13423-017-1243-6
  * [11] Lukas Schumacher, Paul-Christian Bürkner, Andreas Voss, Ullrich Köthe and Stefan T. Radev  “Neural Superstatistics for Bayesian Estimation of Dynamic Cognitive Models”  In _Scientific Reports_ 13.1 Nature Publishing Group, 2023, pp. 13778  DOI: 10.1038/s41598-023-40278-3
  * [12] Randall D Beer  “On the proper treatment of dynamics in cognitive science”  In _Topics in cognitive science_ Wiley Online Library, 2023
  * [13] Nathan J. Evans and Scott D. Brown  “People Adopt Optimal Policies in Simple Decision-Making, after Practice and Guidance”  In _Psychonomic Bulletin & Review_ 24.2, 2017, pp. 597–606  DOI: 10.3758/s13423-016-1135-1
  * [14] Jing-Jing Li, Chengchun Shi, Lexin Li and Anne Collins  “Dynamic noise estimation: A generalized method for modeling noise in sequential decision-making behavior”  In _bioRxiv_ Cold Spring Harbor Laboratory, 2023  DOI: 10.1101/2023.06.19.545524
  * [15] Šimon Kucharský, N.-Han Tran, Karel Veldkamp, Maartje Raijmakers and Ingmar Visser  “Hidden Markov Models of Evidence Accumulation in Speeded Decision Tasks”  In _Computational Brain & Behavior_ 4.4, 2021, pp. 416–441  DOI: 10.1007/s42113-021-00115-0
  * [16] David Gunawan, Guy E Hawkins, Robert Kohn, Minh-Ngoc Tran and Scott D Brown  “Time-evolving psychological processes over repeated decisions.”  In _Psychological review_ 129.3 American Psychological Association, 2022, pp. 438
  * [17] Aaron Cochrane, Chris R. Sims, Vikranth R. Bejjanki, C. Green and Daphne Bavelier  “Multiple Timescales of Learning Indicated by Changes in Evidence-Accumulation Processes during Perceptual Decision-Making”  In _npj Science of Learning_ 8.1 Nature Publishing Group, 2023, pp. 1–10  DOI: 10.1038/s41539-023-00168-9
  * [18] Roger Ratcliff and Jeffrey N. Rouder  “Modeling Response Times for Two-Choice Decisions”  In _Psychological Science_ 9.5 SAGE Publications Inc, 1998, pp. 347–356  DOI: 10.1111/1467-9280.00067
  * [19] Roger Ratcliff and Francis Tuerlinckx  “Estimating Parameters of the Diffusion Model: Approaches to Dealing with Contaminant Reaction Times and Parameter Variability”  In _Psychonomic Bulletin & Review_ 9.3, 2002, pp. 438–481  DOI: 10.3758/BF03196302
  * [20] Nathan J. Evans and Guy E. Hawkins  “When Humans Behave like Monkeys: Feedback Delays and Extensive Practice Increase the Efficiency of Speeded Decisions”  In _Cognition_ 184, 2019, pp. 11–18  DOI: 10.1016/j.cognition.2018.11.014
  * [21] Michael J. Kahana, Eash V. Aggarwal and Tung D. Phan  “The Variability Puzzle in Human Memory”  In _Journal of Experimental Psychology: Learning, Memory, and Cognition_ 44.12 US: American Psychological Association, 2018, pp. 1857–1863  DOI: 10.1037/xlm0000553
  * [22] Nathan J. Evans, Scott D. Brown, Douglas J.. Mewhort and Andrew Heathcote  “Refining the Law of Practice”  In _Psychological Review_ 125.4 US: American Psychological Association, 2018, pp. 592–605  DOI: 10.1037/rev0000105
  * [23] Laura Fontanesi, Sebastian Gluth, Mikhail S. Spektor and Jörg Rieskamp  “A Reinforcement Learning Diffusion Decision Model for Value-Based Decisions”  In _Psychonomic Bulletin & Review_ 26.4, 2019, pp. 1099–1121  DOI: 10.3758/s13423-018-1554-2
  * [24] Adam F. Osth, Anna Jansson, Simon Dennis and Andrew Heathcote  “Modeling the Dynamics of Recognition Memory Testing with an Integrated Model of Retrieval and Decision Making”  In _Cognitive Psychology_ 104, 2018, pp. 106–142  DOI: 10.1016/j.cogpsych.2018.04.002
  * [25] Lukas Schumacher and Andreas Voss  “Duration Discrimination: A Diffusion Decision Modeling Approach”  In _Attention, Perception, & Psychophysics_ 85.2, 2023, pp. 560–577  DOI: 10.3758/s13414-022-02604-1
  * [26] Scott D. Brown, A… Marley, Christopher Donkin and Andrew Heathcote  “An Integrated Model of Choices and Response Times in Absolute Identification”  In _Psychological Review_ 115.2 US: American Psychological Association, 2008, pp. 396–425  DOI: 10.1037/0033-295X.115.2.396
  * [27] Samuel D. McDougle and Anne G.. Collins  “Modeling the Influence of Working Memory, Reinforcement, and Action Uncertainty on Reaction Time and Choice during Instrumental Learning”  In _Psychonomic Bulletin & Review_ 28.1, 2021, pp. 20–39  DOI: 10.3758/s13423-020-01774-z
  * [28] Steven Miletić, Russell J Boag, Anne C Trutti, Niek Stevenson, Birte U Forstmann and Andrew Heathcote  “A New Model of Decision Processing in Instrumental Learning Tasks”  In _eLife_ 10 eLife Sciences Publications, Ltd, 2021, pp. e63055  DOI: 10.7554/eLife.63055
  * [29] C. Beck and E… Cohen  “Superstatistics”  In _Physica A: Statistical Mechanics and its Applications_ 322, 2003, pp. 267–275  DOI: https://doi.org/10.1016/S0378-4371(03)00019-0
  * [30] Christian Beck  “Superstatistics: theory and applications”  In _Continuum mechanics and thermodynamics_ 16 Springer, 2004, pp. 293–304
  * [31] Christoph Mark, Claus Metzner, Lena Lautscham, Pamela L. Strissel, Reiner Strick and Ben Fabry  “Bayesian model selection for complex dynamic systems”  In _Nature Communications_ 9.1, 2018, pp. 1803  DOI: 10.1038/s41467-018-04241-5
  * [32] Roger Ratcliff  “A Theory of Memory Retrieval”  In _Psychological Review_ 85.2 US: American Psychological Association, 1978, pp. 59–108  DOI: 10.1037/0033-295X.85.2.59
  * [33] Stijn Verdonck, Tim Loossens and Marios G. Philiastides  “The Leaky Integrating Threshold and Its Impact on Evidence Accumulation Models of Choice Response Time (RT)”  In _Psychological Review_ 128.2 US: American Psychological Association, 2021, pp. 203–221  DOI: 10.1037/rev0000258
  * [34] Andreas Voss, Klaus Rothermund and Jochen Voss  “Interpreting the Parameters of the Diffusion Model: An Empirical Validation”  In _Memory & Cognition_ 32.7, 2004, pp. 1206–1220  DOI: 10.3758/BF03196893
  * [35] Veronika Lerche and Andreas Voss  “Experimental Validation of the Diffusion Model Based on a Slow Response Time Paradigm”  In _Psychological Research_ 83.6, 2019, pp. 1194–1209  DOI: 10.1007/s00426-017-0945-8
  * [36] Nina R. Arnold, Arndt Bröder and Ute J. Bayen  “Empirical Validation of the Diffusion Model for Recognition Memory and a Comparison of Parameter-Estimation Methods”  In _Psychological Research_ 79.5, 2015, pp. 882–898  DOI: 10.1007/s00426-014-0608-y
  * [37] Roger Ratcliff and Gail McKoon  “The Diffusion Decision Model: Theory and Data for Two-Choice Decision Tasks”  In _Neural computation_ 20.4, 2008, pp. 873–922  DOI: 10.1162/neco.2008.12-06-420
  * [38] Veronika Lerche and Andreas Voss  “Speed–Accuracy Manipulations and Diffusion Modeling: Lack of Discriminant Validity of the Manipulation or of the Parameter Estimates?”  In _Behavior Research Methods_ 50.6, 2018, pp. 2568–2585  DOI: 10.3758/s13428-018-1034-7
  * [39] Kyle Cranmer, Johann Brehmer and Gilles Louppe  “The Frontier of Simulation-Based Inference”  In _Proceedings of the National Academy of Sciences_ 117.48 Proceedings of the National Academy of Sciences, 2020, pp. 30055–30062  DOI: 10.1073/pnas.1912789117
  * [40] Stefan T. Radev, Marvin Schmitt, Lukas Schumacher, Lasse Elsemüller, Valentin Pratz, Yannik Schälte, Ullrich Köthe and Paul-Christian Bürkner  “BayesFlow: Amortized Bayesian Workflows With Neural Networks”  In _Journal of Open Source Software_ 8.89 The Open Journal, 2023, pp. 5702  DOI: 10.21105/joss.05702
  * [41] Sean Talts, Michael Betancourt, Daniel Simpson, Aki Vehtari and Andrew Gelman  “Validating Bayesian Inference Algorithms with Simulation-Based Calibration”, 2020  DOI: 10.48550/arXiv.1804.06788
  * [42] Teemu Säilynoja, Paul-Christian Bürkner and Aki Vehtari  “Graphical Test for Discrete Uniformity and Its Applications in Goodness-of-Fit Evaluation and Multiple Sample Comparison”  In _Statistics and Computing_ 32.2, 2022, pp. 32  DOI: 10.1007/s11222-022-10090-6
  * [43] Daniel J Schad, Michael Betancourt and Shravan Vasishth  “Toward a principled Bayesian workflow in cognitive science.”  In _Psychological methods_ 26.1 American Psychological Association, 2021, pp. 103
  * [44] Andrew Gelman, Aki Vehtari, Daniel Simpson, Charles C Margossian, Bob Carpenter, Yuling Yao, Lauren Kennedy, Jonah Gabry, Paul-Christian Bürkner and Martin Modrák  “Bayesian workflow”  In _arXiv preprint arXiv:2011.01808_ , 2020
  * [45] Lynton Ardizzone, Jakob Kruse, Sebastian Wirkert, Daniel Rahner, Eric W Pellegrini, Ralf S Klessen, Lena Maier-Hein, Carsten Rother and Ullrich Köthe  “Analyzing inverse problems with invertible neural networks”  In _arXiv preprint arXiv:1808.04730_ , 2018  DOI: 10.48550/arXiv.1808.04730
  * [46] David Greenberg, Marcel Nonnenmacher and Jakob Macke  “Automatic posterior transformation for likelihood-free inference”  In _International Conference on Machine Learning_ , 2019, pp. 2404–2414  PMLR
  * [47] Stefan T Radev, Ulf K Mertens, Andreas Voss, Lynton Ardizzone and Ullrich Köthe  “BayesFlow: Learning complex stochastic models with invertible neural networks”  In _IEEE transactions on neural networks and learning systems_ 33.4 IEEE, 2020, pp. 1452–1466
  * [48] Andreas Voss, Markus Nagler and Veronika Lerche  “Diffusion models in experimental psychology”  In _Experimental psychology_ Hogrefe Publishing, 2013
  * [49] Pedro J Gonçalves, Jan-Matthis Lueckmann, Michael Deistler, Marcel Nonnenmacher, Kaan Öcal, Giacomo Bassetto, Chaitanya Chintaluri, William F Podlaski, Sara A Haddad and Tim P Vogels  “Training deep neural density estimators to identify mechanistic models of neural dynamics”  In _Elife_ 9 eLife Sciences Publications, Ltd, 2020, pp. e56261
  * [50] George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed and Balaji Lakshminarayanan  “Normalizing flows for probabilistic modeling and inference”  In _The Journal of Machine Learning Research_ 22.1 JMLRORG, 2021, pp. 2617–2680
  * [51] Simo Särkkä  “Bayesian Filtering and Smoothing”, Institute of Mathematical Statistics Textbooks  Cambridge: Cambridge University Press, 2013  DOI: 10.1017/CBO9781139344203
  * [52] Felix A Gers, Jürgen Schmidhuber and Fred Cummins  “Learning to forget: Continual prediction with LSTM”  In _Neural computation_ 12.10 MIT Press, 2000, pp. 2451–2471
  * [53] Robert E. Kass and Adrian E. Raftery  “Bayes Factors”  In _Journal of the American Statistical Association_ 90.430 Taylor & Francis, 1995, pp. 773–795  DOI: 10.1080/01621459.1995.10476572
  * [54] David JC MacKay  “Information theory, inference and learning algorithms”  Cambridge university press, 2003
  * [55] Daniel W. Heck, Udo Boehm, Florian Böing-Messing, Paul-Christian Bürkner, Koen Derks, Zoltan Dienes, Qianrao Fu, Xin Gu, Diana Karimova, Henk A.. Kiers, Irene Klugkist, Rebecca M. Kuiper, Michael D. Lee, Roger Leenders, Hidde J. Leplaa, Maximilian Linde, Alexander Ly, Marlyne Meijerink-Bosman, Mirjam Moerbeek, Joris Mulder, Bence Palfi, Felix D. Schönbrodt, Jorge N. Tendeiro, Don van den Bergh, Caspar J. Van Lissa, Don van Ravenzwaaij, Wolf Vanpaemel, Eric-Jan Wagenmakers, Donald R. Williams, Mariëlle Zondervan-Zwijnenburg and Herbert Hoijtink  “A Review of Applications of the Bayes Factor in Psychological Research”  In _Psychological Methods_ 28.3, 2023, pp. 558–579  DOI: 10.1037/met0000454
  * [56] Don van Ravenzwaaij and Eric-Jan Wagenmakers  “Advantages Masquerading as “Issues” in Bayesian Hypothesis Testing: A Commentary on Tendeiro and Kiers (2019)”  In _Psychological Methods_ 27.3 US: American Psychological Association, 2022, pp. 451–465  DOI: 10.1037/met0000415
  * [57] Jorge N. Tendeiro and Henk A.. Kiers  “A Review of Issues about Null Hypothesis Bayesian Testing”  In _Psychological Methods_ 24.6, 2019, pp. 774–795  DOI: 10.1037/met0000221
  * [58] Quentin F Gronau, Alexandra Sarafoglou, Dora Matzke, Alexander Ly, Udo Boehm, Maarten Marsman, David S Leslie, Jonathan J Forster, Eric-Jan Wagenmakers and Helen Steingroever  “A tutorial on bridge sampling”  In _Journal of mathematical psychology_ 81 Elsevier, 2017, pp. 80–97
  * [59] Stefan T. Radev, Marco D’Alessandro, Ulf K. Mertens, Andreas Voss, Ullrich Köthe and Paul-Christian Bürkner  “Amortized Bayesian Model Comparison With Evidential Deep Learning”  In _IEEE Transactions on Neural Networks and Learning Systems_ 34.8, 2020, pp. 4903–4917  DOI: 10.1109/TNNLS.2021.3124052
  * [60] Lasse Elsemüller, Martin Schnuerch, Paul-Christian Bürkner and Stefan T. Radev  “A Deep Learning Method for Comparing Bayesian Hierarchical Models”, 2023  arXiv:2301.11873 [stat.ML]
  * [61] Lasse Elsemüller, Hans Olischläger, Marvin Schmitt, Paul-Christian Bürkner, Ullrich Köthe and Stefan T. Radev  “Sensitivity-Aware Amortized Bayesian Inference”  arXiv, 2023  DOI: 10.48550/arXiv.2310.11122
  * [62] M. Healey and Michael J. Kahana  “Is Memory Search Governed by Universal Principles or Idiosyncratic Strategies?”  In _Journal of Experimental Psychology: General_ 143.2 US: American Psychological Association, 2014, pp. 575–596  DOI: 10.1037/a0033715
  * [63] M. Healey and Michael J. Kahana  “A Four-Component Model of Age-Related Memory Change”  In _Psychological Review_ 123.1 US: American Psychological Association, 2016, pp. 23–69  DOI: 10.1037/rev0000015
  * [64] Birte U. Forstmann, Gilles Dutilh, Scott Brown, Jane Neumann, D. von Cramon, K. Ridderinkhof and Eric-Jan Wagenmakers  “Striatum and Pre-SMA Facilitate Decision-Making under Time Pressure”  In _Proceedings of the National Academy of Sciences_ 105.45 Proceedings of the National Academy of Sciences, 2008, pp. 17538–17542  DOI: 10.1073/pnas.0805903105
  * [65] Eric-Jan Wagenmakers, Roger Ratcliff, Pablo Gomez and Gail McKoon  “A Diffusion Model Account of Criterion Shifts in the Lexical Decision Task”  In _Journal of Memory and Language_ 58.1, 2008, pp. 140–159  DOI: 10.1016/j.jml.2007.04.006
  * [66] Sarah K.. Wynton and Jeromy Anglim  “Abrupt Strategy Change Underlies Gradual Performance Change: Bayesian Hierarchical Models of Component and Aggregate Strategy Use”  In _Journal of Experimental Psychology: Learning, Memory, and Cognition_ 43.10 US: American Psychological Association, 2017, pp. 1630–1642  DOI: 10.1037/xlm0000404
  * [67] Scott D. Brown and Andrew Heathcote  “The Simplest Complete Model of Choice Response Time: Linear Ballistic Accumulation”  In _Cognitive Psychology_ 57.3, 2008, pp. 153–178  DOI: 10.1016/j.cogpsych.2007.12.002
  * [68] Qingsong Wen, Tian Zhou, Chaoli Zhang, Weiqi Chen, Ziqing Ma, Junchi Yan and Liang Sun  “Transformers in Time Series: A Survey”  arXiv, 2023  DOI: 10.48550/arXiv.2202.07125

Appendix

## S1 Appendix. Prior distributions

In the following we list the prior distributions we used for all four NSDDM’s.

### DDM Starting Values

For the starting values of the parameter trajectories we used half-normal
distributions with a mean $\mu$ and a standard deviation $\sigma$ denoted as
$\mathcal{HN}(\mu,\sigma)$:

$\displaystyle v_{0}$ $\displaystyle\sim\mathcal{HN}(2.0,2.0)$ $\displaystyle
a_{0}$ $\displaystyle\sim\mathcal{HN}(2.0,1.5)$ $\displaystyle\tau_{0}$
$\displaystyle\sim\mathcal{HN}(0.3,1.0)$

### Random Walk Transition Model

Half-normal distributions were used for the standard deviations of the
Gaussian random walk transition model:

$\displaystyle\sigma_{v}$ $\displaystyle\sim\mathcal{HN}(0.0,0.1)$
$\displaystyle\sigma_{a}$ $\displaystyle\sim\mathcal{HN}(0.0,0.1)$
$\displaystyle\sigma_{\tau}$ $\displaystyle\sim\mathcal{HN}(0.0,0.01)$

We decided to use a relatively narrower prior on $\sigma_{\tau}$ because the
non-decision time parameter is not expected to fluctuate as heavily as the
other two parameters.

### Mixture Random Walk Transition Model

The mixture random walk transition model used the same prior for the Gaussian
random walk as described above. Additionally, Uniform distributions denoted as
$\mathcal{U}$ were used for the mixture proportion parameter $\rho$:

$\displaystyle\rho_{v}$ $\displaystyle\sim\mathcal{U}(0.0,0.2)$
$\displaystyle\rho_{a}$ $\displaystyle\sim\mathcal{U}(0.0,0.1)$

### Levy Flight Transition Model

The Levy flight transition model uses an alpha stable distribution instead of
a Gaussian distribution for the transition. We used the same priors for the
standard deviations as in the random walk and the mixture random walk. The
alpha stable distribution has an additional parameter $\alpha$, which
determines the fatness of the tails. This parameter is bound between $1$ and
$2$. Therefore, we used a Beta distribution denoted as B and added $1$ to the
sampled values:

$\displaystyle\tilde{\alpha}_{v}$ $\displaystyle\sim\text{B}(1.5,1.5)$
$\displaystyle\tilde{\alpha}_{a}$ $\displaystyle\sim\text{B}(2.5,1.5)$
$\displaystyle\alpha_{v}$ $\displaystyle=\tilde{\alpha}_{v}+1$
$\displaystyle\alpha_{a}$ $\displaystyle=\tilde{\alpha}_{a}+1$

### Regime Switching Transition Model

The same prior distributions as for the mixture random walk were used for the
mixture probabilities of the regime switching transition model.

## S2 Appendix. Neural network architectures and training setups

In the following, we outline our implementation of the neural approximators
and the training setup used for model comparison and parameter estimation.

### Model comparison

For model comparison we trained an ensemble of ten neural approximators. Each
approximator consists of a summary network and an inference network. The
summary network is a many-to-one transformer architecture for time series
encoding [68]. The time series transformer has $128$ template and $64$ summary
dimensions. For inference, we use a network that approximates posterior model
probabilities (PMPs) as employed in [60].

We performed offline training for each of the ten neural approximators
separately. The training data consisted of $25\,000$ simulations per model.
Training was performed with $25$ epochs and a batch size of $16$ starting with
an initial learning rate of $0.0005$. The learning rate was adjusted with a
cosine decay from its initial value to $0$.

### Parameter estimation

For parameter estimation we trained one neural approximator for each of the
four NSDDM implementations. Each approximator consists of a hierarchical
summary network as employed in [60] and two inference networks. Three
bidirectional long-short term memory (LSTM) networks were used for the
hierarchical summary network. The number of hidden units were $512$, $256$,
and $128$ respectively.

For inference, we use a composition of two invertible neural networks [47],
one for the low-level and one for the high-level parameters. The network for
the low-level parameters has $8$ coupling layers with an interleaved affine
and spline internal coupling design. The network for the high-level parameters
only differs from the former in its number of coupling layers which is $6$.

Since our simulators can be run fast, the training of the four neural
approximators was performed online, with $75$ epochs, $1\,000$ iterations per
epoch, and a batch size of $16$. Thus, each approximator was trained on
$N=1\,200\,000$ simulated data sets. The initial learning rate was set to
$0.0005$ and was reduced with a cosine decay function to $0$.

## S3 Appendix. Individual analyses

The following section shows the individual specific posterior re-simulations
and parameter estimates for each difficulty level and both conditions
separately. The visualizations are constructed in the vain of Figure 6 in the
main text.

Figure 9: Aggregate results from all models fitted to the data from
participant $1$. The top row illustrates posterior re-simulations as a measure
of the model’s generative performance and absolute goodness-of-fit to the
data. The bottom row depicts parameter estimates of the drift rate and the
threshold parameter from the non-stationary diffusion decision models (NSDDM).
A Empirical and re-simulated response times for each difficulty level and both
conditions. B Empirical and re-simulated proportions of correct choices
(accuracy) for each difficulty level and both conditions separately. C
Posterior estimates of the drift rate parameter for each difficulty level and
both conditions separately. D Posterior estimates of the threshold parameter
for each difficulty level and both conditions separately. Points indicate
medians and the error bars represent the median absolute deviations (MAD)
across individual data and re-simulations. Figure 10: Aggregate results from
all models fitted to the data from participant $2$. The top row illustrates
posterior re-simulations as a measure of the model’s generative performance
and absolute goodness-of-fit to the data. The bottom row depicts parameter
estimates of the drift rate and the threshold parameter from the non-
stationary diffusion decision models (NSDDM). A Empirical and re-simulated
response times for each difficulty level and both conditions. B Empirical and
re-simulated proportions of correct choices (accuracy) for each difficulty
level and both conditions separately. C Posterior estimates of the drift rate
parameter for each difficulty level and both conditions separately. D
Posterior estimates of the threshold parameter for each difficulty level and
both conditions separately. Points indicate medians and the error bars
represent the median absolute deviations (MAD) across individual data and re-
simulations. Figure 11: Aggregate results from all models fitted to the data
from participant $3$. The top row illustrates posterior re-simulations as a
measure of the model’s generative performance and absolute goodness-of-fit to
the data. The bottom row depicts parameter estimates of the drift rate and the
threshold parameter from the non-stationary diffusion decision models (NSDDM).
A Empirical and re-simulated response times for each difficulty level and both
conditions. B Empirical and re-simulated proportions of correct choices
(accuracy) for each difficulty level and both conditions separately. C
Posterior estimates of the drift rate parameter for each difficulty level and
both conditions separately. D Posterior estimates of the threshold parameter
for each difficulty level and both conditions separately. Points indicate
medians and the error bars represent the median absolute deviations (MAD)
across individual data and re-simulations. Figure 12: Aggregate results from
all models fitted to the data from participant $4$. The top row illustrates
posterior re-simulations as a measure of the model’s generative performance
and absolute goodness-of-fit to the data. The bottom row depicts parameter
estimates of the drift rate and the threshold parameter from the non-
stationary diffusion decision models (NSDDM). A Empirical and re-simulated
response times for each difficulty level and both conditions. B Empirical and
re-simulated proportions of correct choices (accuracy) for each difficulty
level and both conditions separately. C Posterior estimates of the drift rate
parameter for each difficulty level and both conditions separately. D
Posterior estimates of the threshold parameter for each difficulty level and
both conditions separately. Points indicate medians and the error bars
represent the median absolute deviations (MAD) across individual data and re-
simulations. Figure 13: Aggregate results from all models fitted to the data
from participant $5$. The top row illustrates posterior re-simulations as a
measure of the model’s generative performance and absolute goodness-of-fit to
the data. The bottom row depicts parameter estimates of the drift rate and the
threshold parameter from the non-stationary diffusion decision models (NSDDM).
A Empirical and re-simulated response times for each difficulty level and both
conditions. B Empirical and re-simulated proportions of correct choices
(accuracy) for each difficulty level and both conditions separately. C
Posterior estimates of the drift rate parameter for each difficulty level and
both conditions separately. D Posterior estimates of the threshold parameter
for each difficulty level and both conditions separately. Points indicate
medians and the error bars represent the median absolute deviations (MAD)
across individual data and re-simulations. Figure 14: Aggregate results from
all models fitted to the data from participant $6$. The top row illustrates
posterior re-simulations as a measure of the model’s generative performance
and absolute goodness-of-fit to the data. The bottom row depicts parameter
estimates of the drift rate and the threshold parameter from the non-
stationary diffusion decision models (NSDDM). A Empirical and re-simulated
response times for each difficulty level and both conditions. B Empirical and
re-simulated proportions of correct choices (accuracy) for each difficulty
level and both conditions separately. C Posterior estimates of the drift rate
parameter for each difficulty level and both conditions separately. D
Posterior estimates of the threshold parameter for each difficulty level and
both conditions separately. Points indicate medians and the error bars
represent the median absolute deviations (MAD) across individual data and re-
simulations. Figure 15: Aggregate results from all models fitted to the data
from participant $7$. The top row illustrates posterior re-simulations as a
measure of the model’s generative performance and absolute goodness-of-fit to
the data. The bottom row depicts parameter estimates of the drift rate and the
threshold parameter from the non-stationary diffusion decision models (NSDDM).
A Empirical and re-simulated response times for each difficulty level and both
conditions. B Empirical and re-simulated proportions of correct choices
(accuracy) for each difficulty level and both conditions separately. C
Posterior estimates of the drift rate parameter for each difficulty level and
both conditions separately. D Posterior estimates of the threshold parameter
for each difficulty level and both conditions separately. Points indicate
medians and the error bars represent the median absolute deviations (MAD)
across individual data and re-simulations. Figure 16: Aggregate results from
all models fitted to the data from participant $8$. The top row illustrates
posterior re-simulations as a measure of the model’s generative performance
and absolute goodness-of-fit to the data. The bottom row depicts parameter
estimates of the drift rate and the threshold parameter from the non-
stationary diffusion decision models (NSDDM). A Empirical and re-simulated
response times for each difficulty level and both conditions. B Empirical and
re-simulated proportions of correct choices (accuracy) for each difficulty
level and both conditions separately. C Posterior estimates of the drift rate
parameter for each difficulty level and both conditions separately. D
Posterior estimates of the threshold parameter for each difficulty level and
both conditions separately. Points indicate medians and the error bars
represent the median absolute deviations (MAD) across individual data and re-
simulations. Figure 17: Aggregate results from all models fitted to the data
from participant $9$. The top row illustrates posterior re-simulations as a
measure of the model’s generative performance and absolute goodness-of-fit to
the data. The bottom row depicts parameter estimates of the drift rate and the
threshold parameter from the non-stationary diffusion decision models (NSDDM).
A Empirical and re-simulated response times for each difficulty level and both
conditions. B Empirical and re-simulated proportions of correct choices
(accuracy) for each difficulty level and both conditions separately. C
Posterior estimates of the drift rate parameter for each difficulty level and
both conditions separately. D Posterior estimates of the threshold parameter
for each difficulty level and both conditions separately. Points indicate
medians and the error bars represent the median absolute deviations (MAD)
across individual data and re-simulations. Figure 18: Aggregate results from
all models fitted to the data from participant $10$. The top row illustrates
posterior re-simulations as a measure of the model’s generative performance
and absolute goodness-of-fit to the data. The bottom row depicts parameter
estimates of the drift rate and the threshold parameter from the non-
stationary diffusion decision models (NSDDM). A Empirical and re-simulated
response times for each difficulty level and both conditions. B Empirical and
re-simulated proportions of correct choices (accuracy) for each difficulty
level and both conditions separately. C Posterior estimates of the drift rate
parameter for each difficulty level and both conditions separately. D
Posterior estimates of the threshold parameter for each difficulty level and
both conditions separately. Points indicate medians and the error bars
represent the median absolute deviations (MAD) across individual data and re-
simulations. Figure 19: Aggregate results from all models fitted to the data
from participant $11$. The top row illustrates posterior re-simulations as a
measure of the model’s generative performance and absolute goodness-of-fit to
the data. The bottom row depicts parameter estimates of the drift rate and the
threshold parameter from the non-stationary diffusion decision models (NSDDM).
A Empirical and re-simulated response times for each difficulty level and both
conditions. B Empirical and re-simulated proportions of correct choices
(accuracy) for each difficulty level and both conditions separately. C
Posterior estimates of the drift rate parameter for each difficulty level and
both conditions separately. D Posterior estimates of the threshold parameter
for each difficulty level and both conditions separately. Points indicate
medians and the error bars represent the median absolute deviations (MAD)
across individual data and re-simulations. Figure 20: Aggregate results from
all models fitted to the data from participant $12$. The top row illustrates
posterior re-simulations as a measure of the model’s generative performance
and absolute goodness-of-fit to the data. The bottom row depicts parameter
estimates of the drift rate and the threshold parameter from the non-
stationary diffusion decision models (NSDDM). A Empirical and re-simulated
response times for each difficulty level and both conditions. B Empirical and
re-simulated proportions of correct choices (accuracy) for each difficulty
level and both conditions separately. C Posterior estimates of the drift rate
parameter for each difficulty level and both conditions separately. D
Posterior estimates of the threshold parameter for each difficulty level and
both conditions separately. Points indicate medians and the error bars
represent the median absolute deviations (MAD) across individual data and re-
simulations. Figure 21: Aggregate results from all models fitted to the data
from participant $13$. The top row illustrates posterior re-simulations as a
measure of the model’s generative performance and absolute goodness-of-fit to
the data. The bottom row depicts parameter estimates of the drift rate and the
threshold parameter from the non-stationary diffusion decision models (NSDDM).
A Empirical and re-simulated response times for each difficulty level and both
conditions. B Empirical and re-simulated proportions of correct choices
(accuracy) for each difficulty level and both conditions separately. C
Posterior estimates of the drift rate parameter for each difficulty level and
both conditions separately. D Posterior estimates of the threshold parameter
for each difficulty level and both conditions separately. Points indicate
medians and the error bars represent the median absolute deviations (MAD)
across individual data and re-simulations. Figure 22: Aggregate results from
all models fitted to the data from participant $14$. The top row illustrates
posterior re-simulations as a measure of the model’s generative performance
and absolute goodness-of-fit to the data. The bottom row depicts parameter
estimates of the drift rate and the threshold parameter from the non-
stationary diffusion decision models (NSDDM). A Empirical and re-simulated
response times for each difficulty level and both conditions. B Empirical and
re-simulated proportions of correct choices (accuracy) for each difficulty
level and both conditions separately. C Posterior estimates of the drift rate
parameter for each difficulty level and both conditions separately. D
Posterior estimates of the threshold parameter for each difficulty level and
both conditions separately. Points indicate medians and the error bars
represent the median absolute deviations (MAD) across individual data and re-
simulations.

## S4 Appendix. Response time time series

In the following, we present the model fit to the whole response time time
series for the remaining $12$ participants.

Figure 23: Model fit to response time (RT) time series. The empirical RT time
series of participant $1$ is shown in black. From trial $1$ to $700$, the
posterior re-simulation (aka retrodictive check) using the best fitting non-
stationary diffusion decision model (NSDDM) for this specific individual are
shown in blue. In this instance, the results stem from a Lévy flight DDM. For
the remaining trials, one-step-ahead predictions are depicted in cyan. Solid
lines correspond to the median and shaded bands to $90\%$ credibility
intervals (CI). The empirical, re-simulated, and predicted RT time series were
smoothed via a simple moving average (SMA) with a period of $5$. The yellow
shaded regions indicate trials where speed was emphasised over accuracy, while
blank white areas denote instances where the opposite emphasis was applied.
Figure 24: Model fit to response time (RT) time series. The empirical RT time
series of participant $2$ is shown in black. From trial $1$ to $700$, the
posterior re-simulation (aka retrodictive check) using the best fitting non-
stationary diffusion decision model (NSDDM) for this specific individual are
shown in blue. In this instance, the results stem from a Lévy flight DDM. For
the remaining trials, one-step-ahead predictions are depicted in cyan. Solid
lines correspond to the median and shaded bands to $90\%$ credibility
intervals (CI). The empirical, re-simulated, and predicted RT time series were
smoothed via a simple moving average (SMA) with a period of $5$. The yellow
shaded regions indicate trials where speed was emphasised over accuracy, while
blank white areas denote instances where the opposite emphasis was applied.
Figure 25: Model fit to response time (RT) time series. The empirical RT time
series of participant $3$ is shown in black. From trial $1$ to $700$, the
posterior re-simulation (aka retrodictive check) using the best fitting non-
stationary diffusion decision model (NSDDM) for this specific individual are
shown in red. In this instance, the results stem from a mixture random walk
DDM. For the remaining trials, one-step-ahead predictions are depicted in
orange. Solid lines correspond to the median and shaded bands to $90\%$
credibility intervals (CI). The empirical, re-simulated, and predicted RT time
series were smoothed via a simple moving average (SMA) with a period of $5$.
The yellow shaded regions indicate trials where speed was emphasised over
accuracy, while blank white areas denote instances where the opposite emphasis
was applied. Figure 26: Model fit to response time (RT) time series. The
empirical RT time series of participant $4$ is shown in black. From trial $1$
to $700$, the posterior re-simulation (aka retrodictive check) using the best
fitting non-stationary diffusion decision model (NSDDM) for this specific
individual are shown in red. In this instance, the results stem from a mixture
random walk DDM. For the remaining trials, one-step-ahead predictions are
depicted in orange. Solid lines correspond to the median and shaded bands to
$90\%$ credibility intervals (CI). The empirical, re-simulated, and predicted
RT time series were smoothed via a simple moving average (SMA) with a period
of $5$. The yellow shaded regions indicate trials where speed was emphasised
over accuracy, while blank white areas denote instances where the opposite
emphasis was applied. Figure 27: Model fit to response time (RT) time series.
The empirical RT time series of participant $5$ is shown in black. From trial
$1$ to $700$, the posterior re-simulation (aka retrodictive check) using the
best fitting non-stationary diffusion decision model (NSDDM) for this specific
individual are shown in red. In this instance, the results stem from a mixture
random walk DDM. For the remaining trials, one-step-ahead predictions are
depicted in orange. Solid lines correspond to the median and shaded bands to
$90\%$ credibility intervals (CI). The empirical, re-simulated, and predicted
RT time series were smoothed via a simple moving average (SMA) with a period
of $5$. The yellow shaded regions indicate trials where speed was emphasised
over accuracy, while blank white areas denote instances where the opposite
emphasis was applied. Figure 28: Model fit to response time (RT) time series.
The empirical RT time series of participant $7$ is shown in black. From trial
$1$ to $700$, the posterior re-simulation (aka retrodictive check) using the
best fitting non-stationary diffusion decision model (NSDDM) for this specific
individual are shown in red. In this instance, the results stem from a mixture
random walk DDM. For the remaining trials, one-step-ahead predictions are
depicted in orange. Solid lines correspond to the median and shaded bands to
$90\%$ credibility intervals (CI). The empirical, re-simulated, and predicted
RT time series were smoothed via a simple moving average (SMA) with a period
of $5$. The yellow shaded regions indicate trials where speed was emphasised
over accuracy, while blank white areas denote instances where the opposite
emphasis was applied. Figure 29: Model fit to response time (RT) time series.
The empirical RT time series of participant $8$ is shown in black. From trial
$1$ to $700$, the posterior re-simulation (aka retrodictive check) using the
best fitting non-stationary diffusion decision model (NSDDM) for this specific
individual are shown in blue. In this instance, the results stem from a Lévy
flight DDM. For the remaining trials, one-step-ahead predictions are depicted
in cyan. Solid lines correspond to the median and shaded bands to $90\%$
credibility intervals (CI). The empirical, re-simulated, and predicted RT time
series were smoothed via a simple moving average (SMA) with a period of $5$.
The yellow shaded regions indicate trials where speed was emphasised over
accuracy, while blank white areas denote instances where the opposite emphasis
was applied. Figure 30: Model fit to response time (RT) time series. The
empirical RT time series of participant $9$ is shown in black. From trial $1$
to $700$, the posterior re-simulation (aka retrodictive check) using the best
fitting non-stationary diffusion decision model (NSDDM) for this specific
individual are shown in blue. In this instance, the results stem from a Lévy
flight DDM. For the remaining trials, one-step-ahead predictions are depicted
in cyan. Solid lines correspond to the median and shaded bands to $90\%$
credibility intervals (CI). The empirical, re-simulated, and predicted RT time
series were smoothed via a simple moving average (SMA) with a period of $5$.
The yellow shaded regions indicate trials where speed was emphasised over
accuracy, while blank white areas denote instances where the opposite emphasis
was applied. Figure 31: Model fit to response time (RT) time series. The
empirical RT time series of participant $10$ is shown in black. From trial $1$
to $700$, the posterior re-simulation (aka retrodictive check) using the best
fitting non-stationary diffusion decision model (NSDDM) for this specific
individual are shown in blue. In this instance, the results stem from a Lévy
flight DDM. For the remaining trials, one-step-ahead predictions are depicted
in cyan. Solid lines correspond to the median and shaded bands to $90\%$
credibility intervals (CI). The empirical, re-simulated, and predicted RT time
series were smoothed via a simple moving average (SMA) with a period of $5$.
The yellow shaded regions indicate trials where speed was emphasised over
accuracy, while blank white areas denote instances where the opposite emphasis
was applied. Figure 32: Model fit to response time (RT) time series. The
empirical RT time series of participant $7$ is shown in black. From trial $1$
to $700$, the posterior re-simulation (aka retrodictive check) using the best
fitting non-stationary diffusion decision model (NSDDM) for this specific
individual are shown in red. In this instance, the results stem from a mixture
random walk DDM. For the remaining trials, one-step-ahead predictions are
depicted in orange. Solid lines correspond to the median and shaded bands to
$90\%$ credibility intervals (CI). The empirical, re-simulated, and predicted
RT time series were smoothed via a simple moving average (SMA) with a period
of $5$. The yellow shaded regions indicate trials where speed was emphasised
over accuracy, while blank white areas denote instances where the opposite
emphasis was applied. Figure 33: Model fit to response time (RT) time series.
The empirical RT time series of participant $13$ is shown in black. From trial
$1$ to $700$, the posterior re-simulation (aka retrodictive check) using the
best fitting non-stationary diffusion decision model (NSDDM) for this specific
individual are shown in blue. In this instance, the results stem from a Lévy
flight DDM. For the remaining trials, one-step-ahead predictions are depicted
in cyan. Solid lines correspond to the median and shaded bands to $90\%$
credibility intervals (CI). The empirical, re-simulated, and predicted RT time
series were smoothed via a simple moving average (SMA) with a period of $5$.
The yellow shaded regions indicate trials where speed was emphasised over
accuracy, while blank white areas denote instances where the opposite emphasis
was applied. Figure 34: Model fit to response time (RT) time series. The
empirical RT time series of participant $14$ is shown in black. From trial $1$
to $700$, the posterior re-simulation (aka retrodictive check) using the best
fitting non-stationary diffusion decision model (NSDDM) for this specific
individual are shown in blue. In this instance, the results stem from a Lévy
flight DDM. For the remaining trials, one-step-ahead predictions are depicted
in cyan. Solid lines correspond to the median and shaded bands to $90\%$
credibility intervals (CI). The empirical, re-simulated, and predicted RT time
series were smoothed via a simple moving average (SMA) with a period of $5$.
The yellow shaded regions indicate trials where speed was emphasised over
accuracy, while blank white areas denote instances where the opposite emphasis
was applied.

## S5 Appendix. Parameter trajectories

In the following, we present the inferred parameter trajectories for the
remaining participants. For each visualisation the model with the highest
posterior model probability for that specific individual was used.

Figure 35: Posterior parameter trajectory inferred with the best fitting NSDDM
of participant $1$ (a Lévy flight DDM in this case) for all three DDM
parameters (drift rate, threshold, and non-decision time) separately. The
yellow shaded areas indicate trials where speed was emphasised over accuracy
and blank white area indicated where the opposite was asked for. In the top
panel, the task difficulty levels sequence is depicted in black lines. Figure
36: Posterior parameter trajectory inferred with the best fitting NSDDM of
participant $2$ (a Lévy flight DDM in this case) for all three DDM parameters
(drift rate, threshold, and non-decision time) separately. The yellow shaded
areas indicate trials where speed was emphasised over accuracy and blank white
area indicated where the opposite was asked for. In the top panel, the task
difficulty levels sequence is depicted in black lines. Figure 37: Posterior
parameter trajectory inferred with the best fitting NSDDM of participant $3$
(a mixture random walk DDM in this case) for all three DDM parameters (drift
rate, threshold, and non-decision time) separately. The yellow shaded areas
indicate trials where speed was emphasised over accuracy and blank white area
indicated where the opposite was asked for. In the top panel, the task
difficulty levels sequence is depicted in black lines. Figure 38: Posterior
parameter trajectory inferred with the best fitting NSDDM of participant $4$
(a mixture random walk DDM in this case) for all three DDM parameters (drift
rate, threshold, and non-decision time) separately. The yellow shaded areas
indicate trials where speed was emphasised over accuracy and blank white area
indicated where the opposite was asked for. In the top panel, the task
difficulty levels sequence is depicted in black lines. Figure 39: Posterior
parameter trajectory inferred with the best fitting NSDDM of participant $5$
(a mixture random walk DDM in this case) for all three DDM parameters (drift
rate, threshold, and non-decision time) separately. The yellow shaded areas
indicate trials where speed was emphasised over accuracy and blank white area
indicated where the opposite was asked for. In the top panel, the task
difficulty levels sequence is depicted in black lines. Figure 40: Posterior
parameter trajectory inferred with the best fitting NSDDM of participant $7$
(a Lévy flight DDM in this case) for all three DDM parameters (drift rate,
threshold, and non-decision time) separately. The yellow shaded areas indicate
trials where speed was emphasised over accuracy and blank white area indicated
where the opposite was asked for. In the top panel, the task difficulty levels
sequence is depicted in black lines. Figure 41: Posterior parameter trajectory
inferred with the best fitting NSDDM of participant $8$ (a Lévy flight DDM in
this case) for all three DDM parameters (drift rate, threshold, and non-
decision time) separately. The yellow shaded areas indicate trials where speed
was emphasised over accuracy and blank white area indicated where the opposite
was asked for. In the top panel, the task difficulty levels sequence is
depicted in black lines. Figure 42: Posterior parameter trajectory inferred
with the best fitting NSDDM of participant $9$ (a Lévy flight DDM in this
case) for all three DDM parameters (drift rate, threshold, and non-decision
time) separately. The yellow shaded areas indicate trials where speed was
emphasised over accuracy and blank white area indicated where the opposite was
asked for. In the top panel, the task difficulty levels sequence is depicted
in black lines. Figure 43: Posterior parameter trajectory inferred with the
best fitting NSDDM of participant $10$ (a mixture random walk DDM in this
case) for all three DDM parameters (drift rate, threshold, and non-decision
time) separately. The yellow shaded areas indicate trials where speed was
emphasised over accuracy and blank white area indicated where the opposite was
asked for. In the top panel, the task difficulty levels sequence is depicted
in black lines. Figure 44: Posterior parameter trajectory inferred with the
best fitting NSDDM of participant $12$ (a Lévy flight DDM in this case) for
all three DDM parameters (drift rate, threshold, and non-decision time)
separately. The yellow shaded areas indicate trials where speed was emphasised
over accuracy and blank white area indicated where the opposite was asked for.
In the top panel, the task difficulty levels sequence is depicted in black
lines. Figure 45: Posterior parameter trajectory inferred with the best
fitting NSDDM of participant $13$ (a Lévy flight DDM in this case) for all
three DDM parameters (drift rate, threshold, and non-decision time)
separately. The yellow shaded areas indicate trials where speed was emphasised
over accuracy and blank white area indicated where the opposite was asked for.
In the top panel, the task difficulty levels sequence is depicted in black
lines. Figure 46: Posterior parameter trajectory inferred with the best
fitting NSDDM of participant $14$ (a Lévy flight DDM in this case) for all
three DDM parameters (drift rate, threshold, and non-decision time)
separately. The yellow shaded areas indicate trials where speed was emphasised
over accuracy and blank white area indicated where the opposite was asked for.
In the top panel, the task difficulty levels sequence is depicted in black
lines.