A model-based approach to trial-by-trial P300 amplitude fluctuations

1. Introduction

The notion of a Bayesian brain is increasingly recognized as providing a distinctive framework for investigating cognitive brain functions (Kersten et al., 2004; Knill and Pouget, 2004; Friston, 2005; Doya et al., 2007). Predictive coding theories of cortical function provide a possible route to the Bayesian brain (Friston, 2002). According to the predictive coding approach, constraints from higher levels of a cortical hierarchy provide contextual guidance to lower levels of processing, providing a theory of how bottom-up evidence is combined with top-down priors to compute the most likely interpretation of sensory data. Specifically, predictive coding theory proposes that an internal representation of the world generates predictions that are compared with stimulus-driven activity to calculate the residual error between the predicted and the actual information. The residual error is then used to update the internal representation so as to minimize the residual error imposed by future stimuli (Friston, 2002, 2005; Spratling, 2010).

The general scheme of predictive coding as a ubiquitous mode of cortical processing offers an instrumental framework for analyzing functional correlates of the P300 event-related brain potential (Sutton et al., 1965; Kopp, 2008). It has long been recognized that fluctuations in P300 amplitude reflect the degree of surprise related to the processing of attended, but unforeseeable sensory events. In particular (Donchin, 1981) argued that P300 amplitude is not crucially determined by the inherent attributes of eliciting events. Instead of that, he ascertained that “surprising events elicit a large P300 component” (p. 498). Squires et al. (1976) had presented a model of P300 amplitude fluctuations, based on the concept of expectancy, which was thought to be determined by three factors: “(i) the memory for event frequency within the prior stimulus sequence, (ii) the specific structure of the prior sequence, and (iii) the global probability of the event” (p. 1144).

More recently, Mars et al. (2008) proposed a computational model of processes underlying the generation of the P300 in which trial-by-trial fluctuations in P300 amplitudes were explained in terms of a Bayesian observer keeping track of the global probabilities of sensory events. The subjective estimates of statistical regularities in the environment were thought to depend crucially on the integration of sensory data over long periods of time. However, the adequacy of this Bayesian observer model is limited, because it cannot account appropriately for the well-documented effects of the recent stimulus sequence on P300 amplitudes (e.g., Squires et al., 1976; Leuthold and Sommer, 1993).

Here we tested these two state-of-the-art models against a newly developed computational model of trial-by-trial P300 amplitude fluctuations by Bayesian model selection (Kass and Raftery, 1995; Raftery, 1995). The new model assumes three additive digital filtering processes, thereby integrating aspects of both state-of-the-art models. Specifically, subjective estimates of statistical regularities in sensory data are kept at short-term and long-term decay time parameters. Further it implements an alternation term (as Squires et al., 1976) as well as uniform initial prior probabilities (as Mars et al., 2008). Our findings show that this new approach provides a superior account of parietally distributed trial-by-trial P300 amplitudes compared to these two state-of-the-art models.

2. Materials and Methods

2.1. Participants, Experimental Design, and Data Acquisition

Sixteen healthy participants [fourteen women, mean age: 20 years; age range 18–23 years; mean handedness (Oldfield, 1971): 74; handedness range −76–100], all with normal or corrected-to-normal visual acuity participated in the experiment. All were recruited from introductory courses at the Department of Psychology at the Technische Universität Braunschweig in return for course credit. Experimental procedures were approved by the local ethics committee and in accordance with the Declaration of Helsinki.

Participants performed a simple two-choice response time [RT] task without feedback about response accuracy in which all stimuli had equal behavioral relevance. This feature of the experimental design constitutes an important difference between this and the classical oddball paradigm (Ritter and Vaughan, 1969) in which participants usually discriminate between task-relevant (target) and irrelevant (standard) stimuli.

The experiment was realized using the Presentation^® software (Neurobehavioral Systems, Albany, CA, USA). Visual stimuli were presented one at a time for 100 ms each, with a stimulus presentation rate of f_s = 2/3 Hz, i.e., one stimulus per 1.5 s. Stimuli were displayed at the center of a CRT monitor (FlexScan T766 19″; Eizo, Hakusan, Ishikawa, Japan) with a refresh rate of 100 Hz at a resolution of 1280 × 1024 pixels against a light gray background. Viewing distance amounted to 1.25 m. Two types of visual stimuli were presented: the stimulus event was either a red or a blue rectangle, each of which subtended approximately 2.75° × 2.25°.

Participants were required to respond to each stimulus with the previously associated button as quickly as possible but not at the expense of accuracy. They used the index finger of both hands (e.g., left button on response to the red rectangle, right button in response to the blue rectangle). Stimulus-response mapping (i.e., [red-left, blue-right] or [red-right, blue-left], respectively) was counterbalanced over participants.

Participants performed twelve blocks of N = 192 trials of the two-choice RT task. The probability of the occurrence of each stimulus event was manipulated between blocks such that the relative probabilities of events were either 0.5 for each event, across six consecutive blocks (1152 trials overall), or [0.3, 0.7], across the remaining six consecutive blocks (1152 trials overall). Stimulus-probability mapping was counterbalanced over participants (i.e., a stimulus color identified the rare (0.3) stimulus in fifty percent of the participants but the frequent stimulus (0.7) in the remaining participants).

The order of the probability manipulation was counterbalanced over participants (probability category [0.5, 0.5] prior to [0.3, 0.7] or vice versa) who were not informed about these probabilities. Participants were informed that the two different stimuli were randomly distributed across blocks. Between the blocks a break was scheduled, participants were free to initiate the subsequent block at their own pace.

A continuous electroencephalogram (EEG) was recorded using a QuickAmps-72 amplifier (Brain Products, Gilching, Germany) and the BrainVision Recorder^® Version 1.02 software (Brain Products, Gilching, Germany) from frontal (F7, F3, Fz, F4, F8), central (T7, C3, Cz, C4, T8), parietal (P7, P3, Pz, P4, P8), occipital (O1, O2), and mastoid (M1, M2) sites. Ag-AgCl EEG electrodes were used which were mounted on an EasyCap (EasyCap, Herrsching-Breitbrunn, Germany). Electrode impedance was kept below 10 kΩ. All EEG electrodes were referenced to average reference during the recording.

For each participant, the actual stimulus sequence of each probability category [0.5, 0.5], and [0.3, 0.7], respectively, was randomized only once in order to enhance the reliability of the sequential trial-by-trial P300 estimates (see below). Thus, each participant received solely one truly random arrangement of trials in each probability category. This arrangement was repeatedly presented across all six blocks of each probability category, unbeknownst to participants. In consequence, sequential P300 estimates could be averaged over the six sequence repetitions per probability category, thereby improving the notoriously low signal-to-noise ratio of single-trial EEG data. Task-related brain activity of a single trial is much more obscured by task-unrelated brain activity than is task-related activity averaged across trials (Blankertz et al., 2002).

Participants were informed about the problem of non-cerebral artifacts, and they were encouraged to reduce the occurrence of movement artifacts (Picton et al., 2000). Ocular artifacts were monitored by means of bipolar pairs of electrodes positioned at the sub and supraorbital ridges (vertical electrooculogram, vEOG) and at the external ocular canthi (horizontal electrooculogram, hEOG). The EEG and EOG channels were subject to a bandpass of 0.01–30 Hz and digitized at 250 Hz sampling rate.

Off-line analysis of the EEG data was performed by means of the BrainVision Analyzer^® Version 2.0.1 software (Brain Products, Gilching, Germany). Careful manual artifact rejection was performed before averaging to discard trials during which eye movements, or any other non-cerebral artifact except blinks, had occurred. Deflections in the averaged EOG waveforms were small indicating that fixation was well maintained in those trials that survived the manual artifact rejection process. Semi-automatic blink detection and the application of an established method for blink artifact removal were employed for blink correction (Gratton et al., 1983). A digital high-pass filter was applied to the data (0.75 Hz cutoff frequency, 48 db/oct) in order to eliminate low-frequency variations in the EEG signal which were associated with the occasional occurrence of electro-dermal artifacts.

The EEG was then divided into epochs of 1000 ms duration, starting 100 ms before stimulus onset. Epochs were corrected using the interval [−100, 0 ms] before stimulus presentation as the baseline. As a start, event-related potential (ERP) waveforms were created (Luck, 2005). ERP waveforms were calculated as trial averages for each participant and for each event probability [i.e., 0.5, 0.3, 0.7], with the exception that those trials in which the participant selected the wrong behavioral response were excluded from averaging.

Thereafter, trial-by-trial P300s were estimated from the EEG data at electrode Pz, where this ERP component is traditionally reported to be maximal (Duncan-Johnson and Donchin, 1977). To estimate trial-by-trial P300 amplitudes, for each participant, the time point at which the averaged P300 waveforms at Pz were modulated maximally by relative stimulus frequency in the [0.3, 0.7] probability category was determined (M = 344 ms, SD = 48 ms; range 280–464 ms). Identifying the P300 in single trials is a notoriously difficult problem, due to the low signal-to-noise ratio of single-trial EEG data (Blankertz et al., 2002). In our study, for each event probability, trial-by-trial P300 estimates were extracted over a temporal window of ±60 ms around the individual time point of maximal modulation (Barceló et al., 2008), thereby completely ignoring latency variability across single trials (Luck, 2005). Albeit this drawback of the method, it was nevertheless chosen in order (1) to keep the testing environment as similar as possible to the procedures employed by Mars et al. (2008), and (2) to improve the reliability of trial-by-trial amplitude measures, in comparison to peak detection measures, akin to previous studies (Debener et al., 2005).

2.2. Conventional Data Analysis

Trial-by-trial P300 estimates, RTs, and error rates were averaged according to the three event probabilities [i.e., 0.3, 0.5, 0.7]. In the [0.5, 0.5] probability category, trial-by-trial P300 estimates were additionally averaged according to eight third-order stimulus sequences (denoted as aaaa, baaa, abaa, aaba, bbaa, abba, baba, bbba), four second-order stimulus sequences (aaa, baa, aba, bba), and two first-order sequences (aa, ba), with up to four consecutive trials (xxxx) = (trial n − 3, trial n − 2, trial n − 1, trial n = eliciting event). Please note that the symbol a simply denotes one of the two possible stimulus events while symbol b signifies the other one in this notation. For example, if a signifies the red rectangle, then b signifies the blue rectangle (and vice versa). In the [0.5, 0.5] probability category, sequential analysis could be collapsed across the two possible stimulus events since both stimuli were equally probable and task-relevant. The same kind of sequential analysis was performed in the [0.3, 0.7] probability category. However, in this experimental condition, a consistently denoted the rare stimulus in the [0.3, 0.7] probability category, whereas b signified the frequent stimulus in the [0.3, 0.7] probability category. We refrained from analyzing the eight third-order stimulus sequences in the [0.3, 0.7] probability category in order to obtain a sufficient number of trials entering the sequential P300 estimates. Thus, trial-by-trial P300 estimates were solely averaged according to eight second-order stimulus sequences (aaa, baa, aba, bba, bbb, abb, bab, aab) and four first-order sequences (aa, ba, bb, ab) in the [0.3, 0.7] probability category, separately for rare and frequent stimuli.

Individual medians of trial-by-trial P300 estimates, RTs, and error rates over the three event probabilities [i.e., 0.3, 0.5, 0.7] as well as the sequential P300 estimates, generated as described above, were submitted to repeated measures analysis of variance (ANOVAs), using the Greenhouse–Geisser correction. The results of the univariate tests are provided, using a format which gives the uncorrected degrees of freedom, and ε in order to compensate for violations of sphericity or equal covariance among all pairs of levels of the repeated measures (Picton et al., 2000). A measure of effect size, ${\eta }_p^2$ (partial eta squared), is also provided.

2.3. State of the Art Models

Let us call

an estimated subjective probability (henceforth simply called subjective probability) that event k ∈ {1, …, K} on trial n ∈ {1, …, N} will be observed, given a sequence ${\text{s}}_1^{n-1}=\left(s\left(1\right),s\left(2\right),\dots ,s\left(n-1\right)\right)$ of n − 1 former stimulus observations. While n is the discrete time index of the consecutive trials, the value N denotes the total number of trials in a block within an experimental probability category for one subject. Note that in (1) stimulus s(n) has not yet been observed, therefore, a subjective probability distribution P_k(n) for all possible stimuli k ∈ {1, …, K} on trial n is of interest. However, once the stimulus k = s(n) on trial n has been observed (which is only a single value k out of set {1, …, K}), the respective subjective probability P_k(n) can be used to calculate the degree of surprise (Shannon and Weaver, 1948; Strange et al., 2005)

Following Mars et al. (2008), we assume the trial-by-trial P300 estimate Y(n)[μV] to be proportional to the surprise I(n)[bit]:

Note that Squires et al. (1976) assumed direct proportionality between the so-called expectancy E_k(n) and the trial-by-trial P300 estimate:

In the following we briefly recapitulate these two well-known state-of-the-art approaches to compute the subjective probability P_k(n) or expectancy E_k(n), which play the role of a dynamically updated prior probability for learning statistical parameters of the stimulus sequence.

2.3.1. Approach by Mars et al. (MAR)

Mars et al. (2008) proposed a Bayesian observer model (henceforth called MAR) without forgetting according to

where

counts the number of occurrences of event k until trial n − 1. The time sequence d_k(ν) holds d_k(ν) = 1 if s(ν) = k, ν = 1,2, …, otherwise d_k(ν) = 0. Note that ${\sum }_{k=1}^K{\tilde{c}}_{\text{L},k}\left(n\right)=n-1.$ As can be easily seen in (5), the subjective probability for event k on trial n = 1 equals a uniform initial prior P_k(1) = 1/K. After many trials (n ≫ K, and ${\tilde{c}}_{\text{L},k}\left(n\right)\gg 1$), the subjective probability approximates $P_k\left(n\right)\approx \frac{{\tilde{c}}_{\text{L},k}\left(n\right)}{n-1},$ i.e., the relative frequency of event k until trial n − 1. Note that the index “L” of the count function ${\tilde{c}}_{\text{L,}k}\left(n\right)$ expresses the long-term memory character of Mars’ model.

2.3.2. Approach by Squires et al. (SQU)

Unlike Mars et al. (2008), Squires et al. (1976) did not formulate a strict computational model to compute the subjective probability P_k(n). Moreover, having investigated solely a K = 2 case, they use the notion of expectancy¹E_k(n) for stimulus k on trial n. While Squires et al. (1976) have described their model, hence called SQU, partly in math, partly in words, in the following we present a complete analytical formulation of their approach, which is straightforward to implement in software. Their empirical formulation of expectancy that event k ∈ {1, 2} will be observed on trial n ∈ {1, …, N} is given as

with three expectancy contributions, namely the assumed-to-be-known global probability P_k, a count function for the short-term memory “S”, and a count function for the alternation expectancy “Δ” (and an additive constant).

2.3.2.1. Short-term memory. The short-term count function is defined as

which is different to (6), since only a limited memory span of N_depth = 5 is covered, and an exponential forgetting factor ${\gamma }_{\text{S}}=e^{-\frac{1}{{\beta }_{\text{S}}}}$ with 0 ≤ γ_S ≤ 1 and time constant 0 ≤ β_S < ∞ is introduced, with γ_S = 0.6 for all probability categories (i.e., β_S = 1.96). Note that the count function (8) depends only on stimulus observations in the recent past.

2.3.2.2. Global probability. The term P_k = P(s(n) = k) in (7) denotes the true global probability of the stimulus being k. It is nothing else but the relative frequency of the stimulus in the current experimental probability category which must be made known to this model.

2.3.2.3. Alternation expectancy. In contrast, the term ${\stackrel{⌣}{c}}_{\Delta ,k}\left(n\right)\in \left\{-3,-2,0,2,3\right\}$ denotes the expectancy w.r.t. alternating stimuli, and how this expectancy is met by the present stimulus s(n). The sign of ${\stackrel{}{c}}_{\Delta ,k}\left(n\right)$ is negative if the stimulus s(n) violates the alternation expectation [i.e., s(n) and s(n − 1) are identical] and positive if the alternation expectation is met [i.e., s(n) and s(n − 1) differ from each other]. The amplitude of ${\stackrel{}{c}}_{\Delta ,k}\left(n\right)$ depends on the number of previous alternations in a row. The formulas for calculating ${\stackrel{}{c}}_{\Delta ,k}\left(n\right)$ are detailed in the Appendix.

2.3.3. Explanatory notes

In summary, (7) provides a model for expectancy that is linearly composed of three contributions: Firstly, the relative frequency P_k of event k which equals the correct global probability throughout probability categories. Note that, in contrast to Mars et al. (2008), the relative frequency P_k is not learned sequentially by experience but assumed to be known by participants. Secondly, a purely predictive limited length (N_depth = 5) exponentially decaying short-term memory [cf. count function ${\stackrel{}{c}}_{\text{S},k}\left(n\right)$ in (8)]. Thirdly, an expectancy contribution in the range $-3\le {\stackrel{⌣}{c}}_{\Delta ,k}\left(n\right)\le +3$ depending on the extent to which a first-order alternation (aa or ba) expectation has been build up and then met/violated within the latest observed N_depth = 5 trials.

2.4. Proposed Digital Filter Model

In this section we present our newly proposed model, inspired by both Mars et al. (2008) and Squires et al. (1976). Our aim is to unify the learned relative frequency estimation of Mars et al. (2008; long-term) with the exponentially decaying short-term memory and alternation expectation capabilities of Squires et al. (1976), and to express the result in terms of a simple new digital filter (DIF) model. Besides an additive probability-normalizing constant 1/C (see Appendix for details), it consists of three additive contributions to subjective probability: a long-term contribution (“L”), a short-term one (“S”), and one term capturing alternations (“Δ”) as depicted in Figure 1:

There are three different count functions used, each represented by a digital filter transfer function H(f ) applied to the common input signal g_k(n), which is given as:

which implicitly contains an initial prior of 1/K at the start of a block of trials², and a “1” wherever a past stimulus s(ν) equals the current stimulus s(n) = k, otherwise a “0.” Note that in contrast to the sequence d_k(ν) as used in Mars et al. (2008) (6) and Squires et al. (1976) (8) we define a model-exciting infinite length signal g_k(ν), ν ∈ {−∞, …, n − 2, n − 1}. The digital filter model yields an output signal P_k(n) as given in (9). The weighting parameters α_L, α_S, α_Δ hold α_L + α_S + α_Δ = 1 and 0 ≤ α_i ≤ 1, i ∈ {L,S,Δ}.

2.4.1. Short-term memory

The block diagram of the infinite impulse response (IIR) digital filter is shown in Figure 2. The respective short-term memory count function can be expressed as

with some normalizing constant C_S and an exponential forgetting factor ${\gamma }_{\text{S}}={\text{e}}^{-\frac{1}{{\beta }_{\text{S}}}}$ with 0 ≤ β_S < ∞, as with count function ${\stackrel{⌣}{c}}_{\text{s},k}\left(n\right)$ in (8). The transfer function of the short-term digital filtering process as described by (11) is depicted in Figure 1 as H_S(f ), and is plotted in Figure 4 as dashed curve, revealing a smooth (i.e., weak) low-pass characteristic. Note that the short-term memory count function (11) can be expressed mathematically equivalent in a recursive form according to (see the Appendix)

initialized with the uniform initial prior c_S,k(0) = 1/K, which in (11) was contained in the values g_k(ν) = 1/K for ν ≤ 0 of (10). The recursive character of (12) becomes apparent by substituting the right hand side of (12) in itself for calculating c_S,k(n − 1). Figure 2 further illustrates the updating process inherent in (12). The input signal defined in (10) and the weights (1 − γ_S) and γ_S guarantee 0 ≤ c_S,k(n) ≤ 1. The equivalence of (11) and (12) and the derivation of C_S = γ_S/(1 − γ_S) are shown in the Appendix.

2.4.2. Long-term memory

The long-term memory count function can be expressed as

with the time-dependent (i.e., dynamic) exponential forgetting factor ${\gamma }_{\text{L},n}\left(\nu \right)={\prod }_{\upsilon =\nu +1}^n{\gamma }_{\text{L},\upsilon }\frac{1-{\gamma }_{\text{L},\upsilon -1}}{1-{\gamma }_{\text{L},\upsilon }}$$\left(\text{using}{\gamma }_{\text{L},\upsilon }=e^{-\frac{1}{{\beta }_{\text{L},\upsilon }}}\right),$ the dynamic normalizing value C_L,n, and the same model-exciting signal g_k(ν) as before (10). The formulas for calculating γ_L,n(ν) and C_L,n are derived in the Appendix. The transfer function of the long-term digital filtering process as described by (13) is depicted in Figure 1 as H_L(f,n). Analog to (12) a recursive function with the same behavior as (13) can be defined as

with the same initial value C_L,k(0) = 1/K. The dynamics of the forgetting factor of the long-term memory are detailed in the Appendix.

Figure 3 illustrates the updating process inherent in (14). The long-term transfer function H_L(f,n) is plotted in Figure 4 as dash-dotted curve (for n = 1) and as solid curve (for n = N = 192), respectively. Inspection of Figure 4 reveals an initially moderate low-pass characteristic which becomes much sharper when the number of trials increases.

2.4.3. Alternation expectation

Finally, our model comprises a count function capturing alternations:

with some normalizing constant C_Δ and the same model-exciting signal g_k(ν) as before (10). In contrast to the short- and long-term IIR filters which both have a low-pass characteristic, this finite impulse response (FIR) filter reveals a high-pass characteristic. Its transfer function is plotted in Figure 4 as ascending dashed curve. The coefficients γ_Δ,n−ν are specified in more detail in the Appendix.

2.4.4. Explanatory notes

There are two important differences to Squires et al. (1976): We propose to use two terms with different time parameters β_S and β_L,n, one accounting for the short-term memory [as in (8)], the other one accounting for a dynamically adapted long-term memory. Secondly, we allow for N_depth → ∞. Moreover, the role of negative trial indices in our model is to define the initial subjective probability distribution, which is P_k(1) = 1/K, k ∈ {1, …, K} in (9), reflecting that participants were not informed about the actual relative frequencies of events over a block of trials.

In summary, the DIF model expresses both a long-term memory contribution, and a short-term memory contribution by exponential decay processes, with uniform initial subjective prior probabilities, and a contribution of alternation expectation. Though there are similarities between (7) and (9), the model of Squires et al. (1976) uses information about the experimental design (P_k) which was actually unknown to the participants. In contrast, the DIF model uses only the information contained in the stimulus sequence as observed by participants, and it always starts with a uniform initial subjective probability distribution P_k(1) = 1/K, k ∈ {1, …, K}, regardless of the actual relative frequencies of events over a block of trials. Finally it should be noted that our new model yields a conceptually well-defined subjective probability, as opposed to an expectancy as in Squires et al. (1976).

2.4.5. Evaluation methods

Following Mars et al. (2008) we compared the DIF to the MAR and SQU models using the log-Bayes factor based on the model evidences. The evidences were approximated using the variational free energy which consists of an accuracy and complexity term, thus enabling the comparison and selection of competing models (Penny et al., 2004; Friston et al., 2007; Penny, 2012). We employed the same three-level hierarchical general linear model (GLM) as Mars et al. (2008). For model fitting and calculation of the model evidences we used parametric empirical Bayes (PEB) from the spm_PEB.m function of the Statistical Parametric Mapping (SPM8) software (Friston et al., 2002, 2007).

The different models (DIF, MAR, and SQU) generate the model-specific surprise³I_ℓ(n) or expectancy E_k,ℓ(n) values as regressors, with the subscript ℓ ∈ {1, …, L = 16} denoting the individual participants and n ∈ {1, …, N = 192} being the discrete time index of the consecutive trials within one block. The first level of the GLM models the measured trial-by-trial P300 estimates Y_ℓ(n) [μV] as a linear function of the surprise I_ℓ(n) [bit] with the intercept ${\theta }_{\ell }^{\left(1\right)}\left[\mu \text{V}\right],$ the slope ${\vartheta }_{\ell }^{\left(1\right)}\left[\mu \text{V}∕\text{bit}\right],$ and an error ${ϵ}_{n,\ell }^{\left(1\right)}\left[\mu \text{V}\right]:$

Note that the fitted model-based P300 estimates then follow

The second level models the participant-specific parameters ${\theta }_{\ell }^{\left(1\right)}\text{and}{\vartheta }_{\ell }^{\left(1\right)}$ as deviations from the corresponding group parameters θ⁽²⁾ [μV] and ϑ⁽²⁾[μV/bit]:

The third level functions as a shrinkage prior on the group parameters. In matrix notation the GLM structure is

Y = [Y_{ℓ=1, …,} Y_ℓ=L]^T ∈ ℝ^{2NL × 1} is a vector concatenating the trial-by-trial P300 estimates for all participants, with []^T being the transpose, making Y a column vector. The participant-specific vector Y_ℓ = [Y_ℓ(n = 1), …, Y_ℓ(n = 2N)] ∈ ℝ^{1 × 2N} contains the trial-by-trial P300 estimates for one participant, averaged over the six sequence repetitions, for both probability categories. The first level design matrix X⁽¹⁾ ∈ ℝ^{2NL × 2L} is block-diagonal with L partitions ${\mathbf{\text{X}}}_{\ell }^{\left(1\right)}=\left[{\mathbf{\text{1}}}_{2N}{\mathbf{\text{I}}}_{\ell }\right]\in {ℝ}^{2N\times 2},$ each of which contains an all-one column vector 1_2N of length 2N, and surprise values I_ℓ = [I_ℓ(n = 1), …, I_ℓ(n = 2N)]^T ∈ ℝ^{2N × 1} as explanatory variables. The second level design matrix X⁽²⁾ = 1 _L⊗I_{2 × 2} = [I_{2 × 2,ℓ=1, …}, I_{2 × 2,ℓ=L}]^T ∈ ℝ^{2L × 2} is the Kronecker product of an all-one column vector 1 _L of length L and an identity matrix I_{2 × 2}. The third level design matrix X⁽³⁾ shall have all-zero elements. The unknown level-one parameters ${\theta }_{\ell }^{\left(1\right)}\text{and}{\vartheta }_{\ell }^{\left(1\right)}$ are assembled in the parameter vector ${\Theta }^{\left(1\right)}={\left[{\theta }_{\ell =1}^{\left(1\right)},{\vartheta }_{\ell =1}^{\left(1\right)},\dots ,{\theta }_{\ell =L}^{\left(1\right)},{\vartheta }_{\ell =L}^{\left(1\right)}\right]}^T\in {ℝ}^{2L\times 1}$. Likewise, the second level parameters θ⁽²⁾ and ϑ⁽²⁾ are assembled in the vector Θ⁽²⁾ ∈ ℝ^{2 × 1}. All errors are assumed to be normally distributed ${\mathbf{\text{E}}}^{\left(j\right)}\sim \mathcal{N}\left(0,{\Sigma }_{ϵ}^{\left(j\right)}\right).$ The covariance is parameterized following ${\Sigma }_{ϵ}^{\left(j\right)}={\lambda }^{\left(j\right)}{\mathbf{\text{I}}}^{\left(j\right)}$, with I^(j) as an identity matrix with the same dimension as the number of rows of the design matrix of the corresponding level X^(j). The hyperparameters λ^(j) are the free parameters of the hierarchical linear model and are estimated using an EM algorithm for maximum likelihood estimation.

The conditional means of the first level parameters ${\mu }_{\Theta |Y}^{\left(1\right)}$ of the posterior densities $\mathcal{N}\left({\mu }_{\Theta |Y}^{\left(j\right)},{\Sigma }_{\Theta |Y}^{\left(j\right)}\right)$ were used as maximum a posteriori point estimates of the parameters for the model fitting for the Figures 6–8 and the calculation of the mean squared error (MSE) and fraction of variance explained (FVE) between ${Ŷ}_{\ell }\left(n\right)$ and Y_ℓ(n) in Table 4 (Friston et al., 2002). The log-evidences or marginal log-likelihoods of the models F = ln(p(Y | ℳ)), with p(Y | ℳ) being the likelihood of the data Y given the model ℳ, were used for model comparison via the log-Bayes factor ln(BF) which is the natural logarithm of the quotient of the model likelihoods or the difference in log-evidence (Kass and Raftery, 1995; Penny et al., 2004; Friston et al., 2007):

If the log-evidence is calculated this way, positive values reflect evidence in favor of the DIF model and negative values in favor of the XXX model (being SQU or MAR), respectively. Values larger than five are considered “very strong” evidence (Kass and Raftery, 1995; Penny et al., 2004). To summarize, in our evaluation method we closely followed (Mars et al., 2008).

2.4.6. DIF model parameter identification

The values for the free model parameters, namely α_L in (9), τ₁ and τ₂ (both (A12)), β_S [for (11)], α_Δ in (9), and γ_Δ,2 in (15), have to be trained on the measured data, i.e., on the trial-by-trial P300 estimates. The model-based P300 estimates are calculated with the same model parameters for all participants and then used for maximization of the DIF model evidence, which is our optimization criterion, as described in Section 2.4.5.

The calculation of the model evidence for the whole range of possible combinations of the parameters with a reasonable resolution is computationally too expensive. For this reason only subsets of parameter combinations were optimized simultaneously with a resolution of 100 values per parameter, which results in 30,000 possible parameter combinations for one iteration. In the first iteration, while optimizing one set of parameters, the not yet optimized parameters were fixed to the center of their respective intervals. In the following iteration, parameters not currently optimized were fixed to the optimal values from the last iteration. In these two iterations a total of 60,000 parameter combinations have been evaluated, and the set with the highest evidence was considered optimal. Note that only a locally optimal parameter combination can be found using this procedure, as many iterations may be necessary for convergence toward the global optimum, if it can be found at all. Table 1 gives an overview over the searched parameter space.

3. Results

3.1. Behavioral Results

RTs showed clear dependence on stimulus probability (0.3, M = 373.4 ms, SE = 6.9 ms; 0.5, M = 353.5 ms, SE = 7.7 ms; 0.7, M = 319.8 ms, SE = 6.4 ms). The slowdown of responding to less probable stimuli was confirmed by an ANOVA on RTs as a function of probability, F_(2,30) = 80.98, p < 0.001, ${\eta }_p^2=0.84,$ ∊ = 0.95. Polynomial contrasts revealed a linear trend, F_(1,15) = 166.79, p < 0.001, ${\eta }_p^2=0.92,$ in the absence of a quadratic trend, F_(1,15) = 3.37, p > 0.05.

Error rates similarly showed clear dependence on stimulus probability (0.3, M = 9.6%, SE = 1.6%; 0.5, M = 4.7%, SE = 0.9%; 0.7, M = 2.1%, SE = 0.4%). The enhanced error proneness in response to less probable stimuli was confirmed by an ANOVA on arcsin-transformed error rates as a function of probability, F_(2,30) = 25.19, p < 0.001, ${\eta }_p^2$ = 0.63, ∊ = 0.65. Polynomial contrasts revealed a linear trend, F_(1,15) = 28.94, p < 0.001, ${\eta }_p^2=0.66$ as well as a quadratic trend, F_(1,15) = 4.86, p < 0.05, ${\eta }_p^2=0.25$.

3.2. Conventional ERP Results

Figure 5 depicts grand-average ERP waveforms (upper panels) and topographic maps (lower panels). Left panels illustrate ERP waveforms at Pz that were obtained in the [0.3, 0.7] probability category. Right panels show third-order sequence effects on ERP waveforms at Pz that were obtained in the [0.5, 0.5] probability category. Note that sequences of four successive stimuli are illustrated, in temporal order (trial n − 3, trial n − 2, trial n − 1, trial n = eliciting event); a signifies a particular stimulus, b the other one. For example, aaaa gives a description of stimulus a being repeated across four consecutive trials (shown as green dashed curve), whereas bbba represents the presentation of stimulus a after having stimulus b repeated across the three immediately preceding trials (shown as black dashed curve).

As can be seen from Figure 5A, trial-by-trial P300 estimates showed clear dependence on stimulus probability (0.3, M = 4.84 μV, SE = 0.72 μV; 0.5, M = 3.51 μV, SE = 0.58 μV; 0.7, M = 2.00 μV, SE = 0.52 μV). P300 augmentation over stimulus improbability was confirmed by an ANOVA on P300 amplitudes as a function of probability, F_(2,30) = 39.88, p < 0.001, ${\eta }_p^2=0.73,$ ∊ = 0.82. Polynomial contrasts revealed a linear trend, F_(1,15) = 55.49, p < 0.001, ${\eta }_p^2=0.79,$ in the absence of a quadratic trend, F_(1,15) = 0.17, p > 0.05.

Sequential P300 estimates in the [0.5, 0.5] probability category (Figure 5B) yielded main effects of first-, F_(1,15) = 6.72, p < 0.05, ${\eta }_p^2=0.31,$ second-, F_(1,15) = 21.04, p < 0.001, ${\eta }_p^2=0.58,$ and third-order sequences, F_(1,15) = 6.89, p < 0.05, ${\eta }_p^2=0.32,$ as well as a significant three-way first- by second- by third-order sequence interaction, F_(1,15) = 6.70, p < 0.05, ${\eta }_p^2=0.31.$ First-order alternations (ba; M = 3.75 μV, SE = 0.57 μV) were associated with enhanced P300 amplitudes compared to first-order repetitions (aa; M = 3.26 μV, SE = 0.59 μV). Likewise, second-order alternations (bxa; M = 3.98 μV, SE = 0.64 μV) were associated with enhanced P300 amplitudes compared to second-order repetitions (axa; M = 3.03 μV, SE = 0.51 μV). Finally, third-order alternations (bxxa; M = 3.64 μV, SE = 0.58 μV) were associated with enhanced P300 amplitudes compared to third-order repetitions (axxa; M = 3.37 μV, SE = 0.57 μV). Separate ANOVAs on sequential P300 estimates were performed in each second-order sequence condition to further parse the three-way interaction. These ANOVAs revealed that the two-way first- by third-order sequence interaction was not significant when the second-order sequence consisted of stimulus repetitions (i.e., when xaxa sequences were included), F_(1,15) = 0.94, p > 0.05, whereas the two-way first- by third-order sequence interaction was significant when the second-order sequence consisted of stimulus alternations (i.e., when xbxa sequences were included), F_(1,15) = 6.34, p < 0.05, ${\eta }_p^2=0.30.$ Further comments on these data are deferred to the Discussion.

Sequential P300 estimates in the [0.3, 0.7] probability category yielded significant main effects of stimulus probability, F_(1,15) = 44.08, p < 0.001, ${\eta }_p^2=0.75$ (0.3, M = 4.59 μV, SE = 0.70 μV > 0.7, M = 2.20 μV, SE = 0.53 μV), of first-order sequence, F_(1,15) = 5.80, p < 0.05, ${\eta }_p^2=0.28$ (ba, M = 3.68 μV, SE = 0.58 μV > aa, M = 3.10 μV, SE = 0.63 μV), and of second-order sequence, F_(1,15) = 11.20, p < 0.01, ${\eta }_p^2=0.43$ (bxa, M = 3.81 μV, SE = 0.66 μV > axa, M = 2.98 μV, SE = 0.54 μV), but without interactions between these factors.

3.3. Model-Based Trial-By-Trial Analysis

The maximum log-Bayes factors in favor of the DIF model over the MAR and SQU models are shown in Table 2. In both cases, this is considered a very strong evidence in favor of the DIF model (Kass and Raftery, 1995; Penny et al., 2004). Table 3 shows the free parameters of the DIF model which were used for calculating the log-Bayes factors in Table 2. Figure 6 illustrates how the log-Bayes factors vary in dependence on the model parameters, with the tip of the “V” marking the parameter combination with the highest evidence. It is important to note the relatively flat tops of the contours implying good generalization capability of the DIF model. Due to computational complexity, only two parameters were optimized simultaneously, as described in Section 2.4.6. The relatively high value of α_L = 0.83 shows that the subjective probability mainly follows the long-term memory. With the identified values for τ_l and τ₂ we get β_L,1 = 40.3 and β_L,192 = 11787, which is further illustrated in Figure A2 in the Appendix. With a short-term memory time constant of β_S = 1.82 and a weight of α_S = 0.12 the influence of recent events to the subjective probability is captured. While the weight of the filter modeling alternation expectancy α_Δ = 0.05 appears to be small, Figure 6C clearly shows the importance of this contribution.

Figure 7 shows timeline plots for one exemplary participant. Figures 7A,B show plots of the expectancy E_k=1(n) and the subjective probability P_k=1(n) of seeing stimulus k = 1 on trial n for the [0.5, 0.5] and [0.3, 0.7] probability category for all competing models. Figures 7C,D show plots of the expectancy E_k=s(n)(n) and the subjective probability P_k=s(n)(n) of seeing the actually occurring stimulus k = s(n) on trial n for both probability categories and all models. The transition from Figures 7A,B to Figures 7C,D illustrates that the subjective probability is traced as a distribution for all possible events k ∈ {1, …, K} simultaneously over all trials and that at the moment of seeing a new stimulus k = s(n) only the corresponding subjective probability of that event k is relevant for the surprise I_ℓ(n) and consequently for the model-based P300 estimate ${\hat{Y}}_{\ell }\left(n\right).$ Figures 7E,F show plots of both the measured and the model-based trial-by-trial P300 estimates Y_ℓ(n) and ${\hat{Y}}_{\ell }\left(n\right)$ for both probability categories, where the latter is calculated according to (17). It is visible that the DIF model estimates are smoother over trials than those of the SQU model but not as undynamic as the estimates of the MAR model, which loses its initial dynamic and becomes almost binary over increasing trial number n. This effect is especially prominent in Figures 7E,F. Similar consecutive trials n elicit a descent in the measured trial-by-trial P300 estimate. For small n all models show this behavior, but for increasing n the MAR model yields nearly constant estimates.

Furthermore the MAR model does not account appropriately for the well-documented sequence effects (Squires et al., 1976). Figure 8 shows the tree diagrams of the measured (Yℓ(n)) and model-based $\left({\hat{Y}}_{\ell }\left(n\right)\right)$ P300 estimates as a function of the preceding stimuli sequences for the different probability categories. The DIF and SQU model are both capable of estimating the envelope and general tree structure quite well, but the SQU model fans out too much for higher order effects for the frequent stimulus in the [0.3, 0.7] and in general in the [0.5, 0.5] probability category, while for the MAR model higher-order effects are nearly non-existent.

As an additional measure of goodness-of-fit of the models Table 4 shows the mean squared error (MSE) and the fraction of variance explained (FVE) of the fitted model predictions ${\hat{Y}}_{\ell }\left(n\right).$ Although differences appear to be somewhat smaller, still the superiority of the DIF model is evident, supporting the log-Bayes factors presented in Table 2.

4. Discussion

We tested three computational models of trial-by-trial P300 amplitudes using Bayesian model selection (Kass and Raftery, 1995; Raftery, 1995). Trial-by-trial P300 amplitude estimates at Pz were obtained in a two-choice RT task. Behavioral data indicated that on average participants reacted slower and committed more errors when they responded to rarely occurring stimuli, consistent with many earlier reports (Miller, 1998). P300 amplitudes showed the expected relationships with stimulus probability. Further, they were influenced by the immediately preceding stimulus sequence, a finding which is also consistent with earlier reports. Thus, our data replicate two of the most ubiquitous P300 findings, notably probabilistic and sequential effects on P300 amplitudes.

The DIF model (9) possesses important advantages over previous models of P300 amplitude fluctuations. It relies completely on mathematical notations and definitions, unlike the notions of expectancy (Squires et al., 1976), global vs. local probability (Squires et al., 1976), temporal probability (Gonsalvez and Polich, 2002), or context updating (Donchin and Coles, 1988). It is a formal model, akin to the MAR model (Mars et al., 2008), but it offers a more parsimonious explanation of trial-by-trial P300 fluctuations. The competitive advantage of the DIF model over the MAR model stems, in large part, from the non-negligible contribution of the short-term traces to subjective estimates of event probabilities, as evidenced by the scarce sensitivity of the MAR model to the sequential effects on P300 amplitudes (Figures 7 and 8).

The SQU model of P300 amplitude fluctuations (Squires et al., 1976) can be considered as a precursor of our DIF model insofar as it comprised memory for event frequencies within the prior stimulus sequence (equivalent to the short-term trace) and event probabilities (loosely related to the long-term trace). Yet, the long-term contribution to the DIF model does not incorporate global event probabilities, P_k, themselves, but rather subjective estimates, c_L,k(n), of these probabilities, these being based on counting observed events in continuously larger samples. Thus, subjective probability estimates are constantly revised while evidence is accumulating, and the DIF model is a pure model of subjective statistical parameters, rather than a mixture of subjective and objective parameters such as the SQU model.

Given the hereby documented superiority of the DIF model over its competitors, we will shortly consider some of its cornerstones. To begin with, it is important to view the DIF model in the context of the processing of event frequencies (Sedlmeier and Betsch, 2002). In particular, the reliable encoding of the frequency with which events occur (Underwood, 1969; Hintzman, 1976) led to the claim that event frequency is automatically encoded in memory, placing only minimal demands on attentional resources (Hasher and Zacks, 1984; Zacks and Hasher, 2002). A variety of representation modes for memory for event frequency can be envisaged. According to multiple-trace views (Hintzman, 1976), a record of individual events is stored such that each attended occurrence of an event results in an independent memory trace. In contrast, according to strength views, each attended event occurrence produces an increment in the strength of a single memory trace or a frequency counter (Alba et al., 1980), supporting the event frequency counter assumption which is inherent in the DIF model. Both, the short-term, c_S,k(n), and the long-term, c_L,k(n), memory traces are frequency counters (11) and (13).

The retention functions describing the short-term and long-term traces in the DIF model are of exponential-decay nature (Lu et al., 2005), differing mainly with regard to their decay half-lives. The dual decay rate assumption is compatible with the fact that short-term and long-term memory functions depend on dissociable neuronal processes (Jonides et al., 2008). Further, recent functional brain imaging data suggest different distributions of cortical responses for short-term and long-term decay functions (Harrison et al., 2011).

The optimal short-term time constant β_S approximated the value of two. Note that the DIF model does not allow for variations in β_S, implying that the short-term contribution to the subjective probability for the appearance of event k on trial n, P_k=s(n)(n), occurs stable over the progression of n. In contrast, β_L,n varies as a function of n (A12), with γ_L,n gradually approximating the value of one. β_L,n (and hence γ_L,n) are relatively low during early trials when compared to late trials within blocks of trials (Figure A2 in Appendix). On the one hand, the long-term quality of the long-term contribution to P_k=s(n)(n) is gradually increasing as a function of n, as revealed by the dynamics of the long-term low-pass filter (Figure 4). In other words, the decay half-life of the long-term trace gradually increases when the observer experiences more and more trials. On the other hand, the recursive formulation of the long-term contribution in the DIF model (14) reveals that the balance between the most recently experienced stimuli (which occurred on trial n − 1, weighted by (1 − γ_L,n−1)) and the counted frequency (weighted by γ_L,n−1) is biased toward recent stimuli during early trials (when γ_L,n−1 is relatively low), but biased toward the counted frequency during late trials (when γ_L,n−1 is relatively high). Thus, the DIF model postulates that the decay half-life of the long-term trace evolves dynamically with the amount of experience. The observer is modeled to rely more and more on environmental experience rather than on prior assumptions, possibly reflecting the fact that the exploitation of statistical redundancy becomes gradually more reliable with progressing time (Barlow, 2001).

Visual working memory (vWM, Baddeley, 2003) can maintain representations of only around four objects at any given moment (Cowan, 2005). The surprisingly limited vWM capacity offers a rationale for the assumption of a capacity-limited alternation term in the DIF model (15), c_Δ,k(n). Its finite impulse response (FIR) characteristic resembles the alternation term in the SQU model (Squires et al., 1976). This FIR high-pass filter searches for alternation patterns over short sequences of trials (such as those in abab and in baba sequences). The discovery of such patterns leads one to expect the completion of the pattern in the upcoming trial, an expectation which will be confirmed in the ababa sequence, but will be disconfirmed in the babaa sequence. It is important that alternation expectation appears conditional upon the detection of alternation patterns in vWM, as revealed by the larger P300 amplitudes in response to ba sequences compared to aa sequences (Figure 8). Thus, the detection of alternation patterns in vWM entrains alternation expectation, as evidenced by our data (Figure 5B, see also Jentzsch and Sommer, 2001; Ford et al., 2010). Further, the effects of pattern completion (such as baba sequences) vs. pattern violation (such as abaa sequences) might underlie the first- by second- by third-order sequence interaction which we identified in the [0.5, 0.5] probability category.

While the effects of alternation expectation are non-negligibly measurable at Pz, visual inspection of our data at more anterior electrode sites suggested that these alternation expectancy effects might possess a more anterior, i.e., P3a-like scalp topography than the proper event frequency effects which showed the typical, P3b-like scalp topography (Polich, 2007; Duncan et al., 2009). While these initial observations ask for multi-channel data analyses, the present modeling work should be mainly considered as a model for P3b generation, since the task procedures (all stimuli required a button press) and the ERP waveforms [i.e., their scalp distribution (cf. Figure 5) and peak latencies] favor an interpretation in terms of predominant P3b-potentials.

The DIF model offers a digital filtering account of multiple memory systems in the brain (Figures 1–4). Specifically, the DIF model characterizes frequency memory as two digital first-order infinite impulse response (IIR) low-pass filters, one filter with an experience-invariant short-term exponential-decay function (Figure 2), and another filter with an experience-dependent long-term exponential-decay function, such that the low-pass characteristic becomes progressively apparent as the amount of experience increases (Figure 3). Moreover, vWM is conceptualized as an additional fourth-order finite impulse response (FIR) high-pass filter (Figure 9). The input signal g_k(n) in (10) to all three filters is a binary representation of the stimulus sequence, with all samples prior to the first trial filled with the uniform initial prior.

Our theory of variation in trial-by-trial P300 amplitudes bears implications on the nature of cortical processing. It is in agreement with the predictive coding approach (Friston, 2002, 2005; Spratling, 2010). Viewed from this perspective, predictive surprise – and hence trial-by-trial P300 amplitude – is proportional to the residual error between top-down priors and bottom-up sensory evidence. Predictive coding theory is a successful guiding model for functional neuroimaging and electrophysiological studies of sensory cortical processing (Summerfield et al., 2006; Garrido et al., 2009; Summerfield and Egner, 2009; Egner et al., 2010; Rauss et al., 2011; Winkler and Czigler, 2011). Further, the DIF model is a Bayesian model of cortical processing (Knill and Pouget, 2004; Friston, 2005). It conceives performance on our two-choice RT task as sequential Bayesian learning (MacKay, 2003), with initial prior knowledge being conceptualized as a uniform prior probability distribution (10), consistent with Laplace’s Principle of Indifference.

It is important to note that we do not claim that the observed P300 modulations were exclusively related to predictive surprise over sensory input, since in the present task design the probabilities of sensory events were mirrored on probabilities of motor responses, as each stimulus was mapped onto a distinct motor response. Thus, particular stimuli also called for particular motor programs, and it could be that the observed P300 modulations are related to predictive surprise over motor responses. We deliberately leave it open whether the observed P300 modulations were due to surprise conveyed by the visual stimuli, or whether they were related to surprise associated with the selection of a motor response, given a visual stimulus on each trial (Barceló et al., 2008; O’Connell et al., 2012).

The DIF model specifies how predictive surprise determines trial-by-trial P300 amplitudes, seemingly representing barely more than a re-iteration of Donchin’s (1981) surprise hypothesis of P300 amplitude fluctuations. However, one should not confuse predictive surprise, as defined by the DIF model, with Bayesian surprise (Ostwald et al., 2012). Bayesian surprise numeralizes the divergence between P₁(n), …, P_K(n) and P₁(n + 1), …, P_K(n + 1), i.e., the divergence between probability distributions across successive trials which can be computed using the Kullback–Leibler metric (Baldi and Itti, 2010). Bayesian surprise thus quantifies the revision of the internal model of the world, given stimulus s(n), whereas predictive surprise I(n) in (2) refers to the unpredictability of s(n), given the internal model immediately before observing s(n). To conclude, we propose a formal computational model of predictive surprise, along with a strategy for testing the model’s ability to predict trial-by-trial P300 amplitudes.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Thanks are due Sandra Mann for her help in data collection and analysis. We also thank Herbert Hoijtink, University of Utrecht, NL, for his valuable suggestions and advice concerning Bayesian model selection.

Footnotes

References

Appendix

Alternation Expectation of the SQU Model

This Appendix gives computational details on the expectation w.r.t. alternating stimuli ${\stackrel{⌣}{c}}_{\Delta ,k}\left(n\right)$ as required in (7) and described in words in Squires et al. (1976). The expectation w.r.t. alternating stimuli, and how this expectancy is met by the present stimulus s(n) is given as:

with N_depth = 5. The use of the minimum function ensures that an expectation for alternation requires at least N_alt = 2 consecutive previous stimulus alternations. Following Squires et al. (1976), its sign

is negative (u_k(n) = −1) if stimulus s(n) violates the alternation expectation [i.e., s(n) and s(n − 1) are identical]. On the other hand, if the alternation expectation is met [i.e., s(n) and s(n − 1) differ from each other], the sign is positive: u_k(n) = +1. The number of previous alternations in a row constitutes the amplitude of the expectation, which is given by

with⁴ $\mathcal{N}=\left\{1,2,\dots ,\left(N_{\text{depth}}-2\right)\right\}.$

Alternation Expectation of the DIF Model

This Appendix gives details on the coefficients γ_Δ,n−ν and the probability normalizing constants C and C_Δ. Figure 9 shows the block diagram of this fourth order finite impulse response filter. In order to reduce the complexity of the model by keeping the number of independent parameters small, only the coefficient γ_Δ,2 is chosen freely within the range 0 ≤ γ_Δ,2 ≤ γ_Δ,max. The effect of the multiplicative constant C_Δ and additive constant C (as shown in Figure 1) is merely to ensure 0 ≤ c_Δ,k(n) + 1/C ≤ 1. Thus, γ_Δ,max can be set to an arbitrary value if some constraints are met: The filter coefficients are set following γ_Δ,1 = − γ_Δ,2, γ_Δ,3 = − γ_Δ,4, and γ_Δ,4 = γ_Δ,max − γ_Δ,2. The normalizing constants have to be set according to: C_Δ = γ_Δ,max + γ_Δ,2 + γ_Δ,4 = 2γ_Δ,max and $C=\frac{{\gamma }_{\Delta ,max}+{\gamma }_{\Delta ,2}+{\gamma }_{\Delta ,4}}{{\gamma }_{\Delta ,max}}=2.$ We chose γ_Δ,max = 1.

On Count Functions and Digital Filters (DIF Model)

In this Appendix we present some mathematical details of the DIF model, its constants, and its relation to the count functions. We will show that the count functions (11) for c_S,k(n) and (13) for c_L,k(n) are equivalent to their recursive formulations (12) and (14), respectively. To this end we will transform the block diagram of Figure 3 to obtain the count functions from the new resulting block diagrams.

Figure A1A shows the block diagram of Figure 3 according to (14) with some notation omitted. A delay of one trial is represented by , the input signal is defined as g_k(n), the output signal as c_L,k(n), (1 − γ_L,n−1) and γ_L,n−1 are time-dependent values. The multiplier (1 − γ_L,n−1) can be moved to the right yielding the block diagram shown in Figure A1B. Figure A1C shows the block diagram after moving (1 − γ_L,n−1) even further to the right. Finally, when moving the multiplier γ_L,n−1/(1 − γ_L,n−1) to the right of the lower branch delay unit, the time dependency has to be accounted for as is shown in Figure A1D. For the long-term filter we obtain an effective filter coefficient

while for the short-term filter this simplifies to

The time-varying discrete-time impulse response of the long-term filter depicted in Figure A1D is given as

with n(>ν) being the currently observed trial, and ν being the trial when the stimulus initiated the impulse response. We make use of the step function

Substituting (A4) in (A6) yields

with the dynamic normalizing value $C_{\text{L},n}=\frac{{\gamma }_{\text{L},n}}{1-{\gamma }_{\text{L},n}}.$ The input-output relation of this time-varying linear system is given by (c.f. Claasen and Mecklenbräuker, 1982; Prater and Loeffler, 1992):

Due to the step function (A7) in the impulse response (A6) only stimuli at trials ν ≤ n − 1 contribute to the output at trial n and we obtain

which is exactly the long-term count function (13) with ${\gamma }_{\text{L},n}\left(\nu \right)={\prod }_{\upsilon =\nu +1}^n{\gamma }_{\text{L},\upsilon }^{\prime }={\prod }_{\upsilon =\nu +1}^n{\gamma }_{\text{L},\upsilon }\frac{1-{\gamma }_{\text{L},\upsilon -1}}{1-{\gamma }_{\text{L},\upsilon }},$ and g_k(ν), as before. Note that for the short-term filter we obtain ${\prod }_{\upsilon =\nu +1}^n{\gamma }_{\text{S}}={\gamma }_{\text{S}}^{n-\nu }$ and C_S = γ_S/(1 − γ_S), which simplifies (A10) to the discrete-time convolution

which is the short-term count function (11). Thus, equivalence has been shown between the count functions (11) and (13), and their simple recursive formulations (12) and (14), respectively. The behavior of the dynamic long-term time value β_L,n follows:

with normalized time constants τ₁, τ₂, controlling the speed of transition from reliance on prior assumptions to experience. The time value 0 ≤ β_L,n < ∞ holds β_L,n > β_S. The effect of this dynamic formulation of β_L,n is further illustrated in Figure A2 which shows the values of β_L,n and the corresponding forgetting factor ${\gamma }_{\text{L},n}=e^{-\frac{1}{{\beta }_{\text{L},n}}}$ for trials n ∈ {1, … N} with N = 192.

Footnote