Incorporating advanced language models into the P300 speller using particle filtering

SWLDA is a classification algorithm that selects a set of signal features to include in a discriminant function [23]. Signals are assigned labels based on two classes: those corresponding to flashes containing the attended character and those without the attended character. The algorithm uses ordinary least-squares regression to predict class labels for the training set. It then adds the most significant features in the forward stepwise analysis and removes the least significant features in the backward analysis step. These steps are repeated until either the target number of features was met or it reached a state where no features were added or removed [10].

The score for each flash in the test set can then be computed as the dot product of the feature weight vector, w, with the features from that trial's signal, $y_{t}^{i}.$ Traditionally, the score for each possible next character, ${{x}_{t}},$ is computed as the sum of the individual scores for flashes that contain that character:

$$\begin{eqnarray}&&g\left( {{x}_{t}} \right)=\mathop \sum \limits_{i:{{x}_{t}}\in A_{t}^{i}} y_{t}^{{{i}^{T}}}w,\end{eqnarray} \tag{ 1 }$$

where $A_{t}^{i}$ is the set of characters illuminated for the $i$th flash for character $t$ in the sequence. After a predetermined set of stimuli, the character with the highest score is selected.

It has been shown that scores can be approximated as independent samples from a Gaussian distribution given the target character [12]:

$$\begin{eqnarray}&&p\left( {{y}_{t}}|{{x}_{t}} \right)\propto \mathop \prod \limits_{i} f\left( y_{t}^{i}|x_{t}^{(L)} \right)\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}f\left( y_{t}^{i}|{{x}_{t}} \right)=\left\{ \begin{array}{ccccccccccccccc} \frac{1}{\sqrt{2\pi \sigma _{a}^{2}}}{{{\rm e}}^{-\frac{1}{2\sigma _{a}^{2}}{{\left( y_{t}^{i}-{{\mu }_{a}} \right)}^{2}}}} & {\rm if}\;{{x}_{t}}\in A_{t}^{i}, \\ \frac{1}{\sqrt{2\pi \sigma _{n}^{2}}}{{{\rm e}}^{-\frac{1}{2\sigma _{n}^{2}}{{\left( y_{t}^{i}-{{\mu }_{n}} \right)}^{2}}}} & {\rm if}\;{{x}_{t}}\notin A_{t}^{i}, \\ \end{array} \right.\end{eqnarray} \tag{ 3 }$$

where ${{\mu }_{a}},$ $\sigma _{a}^{2},$ ${{\mu }_{n}}$ and $\sigma _{n}^{2}$ are the means and variances of the distributions for the attended and non-attended flashes, respectively. The conditional probability of a target given the EEG signal and typing history can then be found:

$$\begin{eqnarray}&&p\left( {{x}_{t}}|{{y}_{t}},{{x}_{t-1}},\ldots ,{{x}_{1}} \right)\propto p\left( {{y}_{t}}|{{x}_{t}} \right)p\left( {{x}_{t}}|{{x}_{t-1}},\ldots ,{{x}_{1}} \right),\end{eqnarray} \tag{ 4 }$$

where $p\left( {{x}_{i}}|{{x}_{i-1}},\ldots ,{{x}_{1}} \right)$ is the prior probability of character ${{x}_{i}}$ determined from a language model. In the simplest case, a uniform prior probability is used, simplifying the posterior to

$$\begin{eqnarray}&&p\left( {{x}_{t}}|{{y}_{t}},{{x}_{t-1}},\ldots ,{{x}_{1}} \right)=p\left( {{x}_{t}}|{{y}_{t}} \right)\propto p\left( {{y}_{t}}|{{x}_{t}} \right).\end{eqnarray} \tag{ 5 }$$

Dynamic classification was implemented by setting a threshold probability, ${{p}_{{\rm thresh}}},$ to determine when a decision should be made. The program flashed characters until either ${{{\rm max} }_{{{x}_{t}}}}\;p\left( {{x}_{t}}|{{y}_{t}} \right)\geqslant {{p}_{{\rm thresh}}}$ or the number of sets of flashes reached the maximum (15). The classifier then selected the character that satisfied ${\rm argma}{{{\rm x}}_{{{x}_{t}}}}\;p\left( {{x}_{t}}|{{y}_{t}} \right).$ In offline analysis, the speeds, accuracies, and bit rates were found for values of ${{p}_{{\rm thresh}}}$ between 0 and 1 in increments of 0.01 and the threshold probability that maximized the bit rate was chosen for each subject.

N-gram models are simple representations of language that produce probability distributions based on relative character frequency in a training corpus. States in these models represent characters in the alphabet with transitions based on the probability of their co-occurrence in text. For example, a state with the character 'Q' would have a high probability transition to a state with the character 'U' because the string 'QU' is relatively common in English text while other strings beginning with 'Q' are much less common. In this study, the second order Markov assumption is used, so prior probabilities can determined from relative trigram counts in the Brown English language corpus [24]. The prior probability of character, ${{x}_{t}},$ given the previous characters, ${{x}_{1:t-1}},$ is

$$\begin{eqnarray}&&p\left( {{x}_{t}}|{{x}_{1:t-1}} \right)=p\left( {{x}_{t}}|{{x}_{t-2}},{{x}_{t-1}} \right)=\frac{c\left( {{x}_{t-2}},{{x}_{t-1}},{{x}_{t}} \right)}{c\left( {{x}_{t-2}},{{x}_{t-1}} \right)},\end{eqnarray} \tag{ 6 }$$

where c('a', 'b', 'c') denotes the number of times the string 'abc' occurs in the corpus.

In the Viterbi algorithm, a typed word is simply a sequence of states of the Markov process, $x=({{x}_{0}},\ldots ,{{x}_{n}}).$ The goal of the algorithm is to search all possible state sequences to determine the most probable output given the known language model and the observed EEG signal. At any given time $t,$ the Viterbi algorithm finds the highest probability path to each possible value for the target character, ${{x}_{t}}.$ It then determines the highest probability of those paths and returns it to the user. The lower probability paths are retained for calculation of the probability distribution for the next time step.

The EEG response at time $t$ is dependent only on the current target character and governed by the conditional probability, $p\left( {{y}_{t}}|{{x}_{t}} \right).$ At each time step, the joint probability of being in state ${{x}_{t}}$ and previously being in state ${{x}_{t-1}}$ can then be computed by

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} p\left( {{x}_{t}},{{x}_{t-1}}|{{y}_{t}},\ldots ,{{y}_{1}} \right)\propto p\left( {{y}_{t}}|{{x}_{t}} \right)\mathop \sum \limits_{{{x}_{t-2}}} p\left( {{x}_{t}}|{{x}_{t-2}},{{x}_{t-1}} \right) \\ \quad \times p\left( {{x}_{t-1}},{{x}_{t-2}}|{{y}_{t-1}},\ldots ,{{y}_{1}} \right), \\ \end{array}\end{eqnarray} \tag{ 7 }$$

where $p\left( {{y}_{t}}|{{x}_{t}} \right)$ is computed as in equation (2) and $p\left( {{x}_{t}}|{{x}_{t-2}},{{x}_{t-1}} \right)$ is determined by the language model (equation (6)). The total probability of the target character ${{x}_{t}}$ can then be computed by summing over all possible previous states ${{x}_{t-1}}$

$$\begin{eqnarray}&&p\left( {{x}_{t}}|{{y}_{t}},\ldots ,{{y}_{1}} \right)\propto \mathop \sum \limits_{{{x}_{t-1}}} p\left( {{x}_{t}},{{x}_{t-1}}|{{y}_{t}},\ldots ,{{y}_{1}} \right).\end{eqnarray} \tag{ 8 }$$

A back pointer to a previous state $\left\langle {{x}_{t-1}},{{x}_{t-2}} \right\rangle$ is saved representing the highest probability path to the state $\left\langle {{x}_{t}},{{x}_{t-1}} \right\rangle ,$ determined using the Viterbi algorithm

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} {{V}_{t}}\left( {{x}_{t}},{{x}_{t-1}} \right)=\mathop{{\rm max} }\limits_{{{x}_{t-2}}}p\left( {{y}_{t}}|{{x}_{t}} \right)p\left( {{x}_{t}}|{{x}_{t-2}},{{x}_{t-1}} \right) \\ \quad \quad \quad \quad \quad \;\;\times {{V}_{t-1}}\left( {{x}_{t-1}},{{x}_{t-2}} \right). \\ \end{array}\end{eqnarray} \tag{ 9 }$$

As in the previous method, characters are selected when a threshold probability is reached (i.e., ${{{\rm max} }_{{{x}_{t}}}}\;p\left( {{x}_{t}}|{{y}_{t}},\ldots ,{{y}_{0}} \right)\geqslant {{p}_{{\rm thresh}}})$ and the character ${{x}_{t}}$ is selected that satisfies ${\rm argma}{{{\rm x}}_{{{x}_{t}},{{x}_{t-1}}}}\;p\left( {{x}_{t}},{{x}_{t-1}}|{{y}_{t}},\ldots ,{{y}_{0}} \right).$ The back pointers are then followed from the selected character to find the optimal string up to the current time step. This string is then presented as the output of the system up to the current time. If this string differs from the previous output, the previous characters are assumed to be errors and are replaced by those in the current string.

The probabilistic automaton models the English language by creating states for every substring that starts a word in the corpus. Thus, the word 'the' would result in three states: 't,' 'th,' and 'the'. The start state is the root node, ${{x}_{0}},$ which corresponds to a blank string. Each state then links to every state that has a string that is a superstring that is one character longer. Thus, the state 't' will link to the states 'th' and 'to' (figure 1).

Zoom In Zoom Out Reset image size

States that represent completed words contain links back to ${{x}_{0}}$ to begin a new word. The state 'the,' for instance, links to the terminal state 'the_' which includes the valid English word 'the' and a space indicating the intent to start a new word. The state 'the' also contains links to nonterminal states such as 'them' because it begins other English words in addition to being a complete word.

Similar to the trigram model, the transition probabilities are determined by the relative frequencies of words starting with the states' substrings in the Brown English language corpus [24]

$$\begin{eqnarray}&&p\left( {{x}_{0:t}}|{{x}_{0:t-1}} \right)=\frac{c\left( {{x}_{1}},\ldots ,{{x}_{t-1}},{{x}_{t}} \right)}{c\left( {{x}_{1}},\ldots ,{{x}_{t-1}} \right)},\end{eqnarray} \tag{ 10 }$$

where c('a', 'b') denotes the number of occurrences of a word that starts with the string 'ab' in the corpus. When $t=1,$ $c(\;)$ represents the total number of words in the corpus and $p\left( {{x}_{0:1}}|{{x}_{0}} \right)$ is the fraction of words in the copus that begin with character 'x₁'. Similarly, the probability that a word ends and the model transitions back to the root, ${{x}_{0}},$ is the ratio of the number of occurrences of complete words consisting of a string to the total number of occurrences of words beginning with that string

$$\begin{eqnarray}&&p\left( {{x}_{0}}|{{x}_{0:t-1}} \right)=\frac{c\left( {{x}_{1}},\ldots ,{{x}_{t-1}}, \right)}{c\left( {{x}_{1}},\ldots ,{{x}_{t-1}} \right)},\end{eqnarray} \tag{ 11 }$$

where c('a', 'b', ' ') is the number of occurrences of the word 'ab' in the corpus.

The SIR method estimates the probability distribution over possible outputs by sampling a batch of possible realizations of the model called particles. Each of these particles moves through the language model independently, based on the model transition probabilities. Low probability realizations are removed and replaced by more likely realizations by resampling the particles based on weights derived from the observed EEG response. The algorithm estimates the probability distribution of the possible output strings by computing a histogram of the particles after they have moved through the model.

Each particle $j$ consists of a link to a state in the language model, ${{x}^{(j)}};$ a string consisting of the particle's state history; and a weight, ${{w}^{(j)}}.$ When the system begins, a set of $P$ particles is generated and each is associated with the root node, ${{x}_{0}},$ with an empty history and a weight equal to $1/P.$ At the start of a new character $t,$ a sample ${{x}_{t}}^{(j)}$ is drawn for each particle, $j,$ from the proposal distribution defined by the language model's transition probabilities from the particle's history, ${{x}_{0:t-1}}^{\left( j \right)}$

$$\begin{eqnarray}&&{{x}_{0:t}}^{(j)}\sim \;p\left( \left. {{x}_{0:t}} \right|{{x}_{0:t-1}}^{(j)} \right),\end{eqnarray} \tag{ 12 }$$

where $p\left( {{x}_{0:t}}|{{x}_{0:t-1}}^{(j)} \right)$ is provided by the language model (equation (10)). When a particle transitions between states, its pointer changes from the previous state in the model, ${{x}_{0:t-1}}^{(j)},$ to the new state ${{x}_{0:t}}.$ The history for each particle, ${{x}_{0:t}}^{(j)},$ is stored to represent the output character sequence associated with that particle. After each stimulus response, the probability weight is computed for each of the particles

$$\begin{eqnarray}&&{{w}_{t}}^{(j)}\propto \;p\left( {{y}_{t}}|{{x}_{t}}^{(j)} \right)\propto \mathop \prod \limits_{i} f\left( y_{t}^{i}|{{x}_{t}}^{(j)} \right),\end{eqnarray} \tag{ 13 }$$

where $f\left( y_{t}^{i}|{{x}_{t}}^{(j)} \right)$ is computed as in equation (2). The weights are then normalized and the probability of the current character is found by summing the weights of all particles that end in that character

$$\begin{eqnarray}&&p\left( {{x}_{0:t}}|{{y}_{1:t}} \right)=\mathop \sum \limits_{k} w_{t}^{\left( k \right)}\delta _{{{x}_{t}}}^{{{x}_{t}}^{\left( k \right)}},\end{eqnarray} \tag{ 14 }$$

where δ is the Kronecker delta. If the maximum probability is above a threshold, the particle with the maximum weight is selected and its history is used as the output text. As in the HMM algorithm, if characters in this output differ from the previous output text, the previous characters are assumed to be errors and are replaced by those in the current text. A new batch of particles, ${\boldsymbol{x}} _{t}^{*},$ are then sampled from the current particles, ${{{\boldsymbol{x}} }_{t}},$ based on the weight distribution, ${{{\boldsymbol{w}} }_{t}}.$ Each of the new particles are then assigned an equal weight $w_{t}^{*\left( j \right)}=1/P.$ The subject then moves on to the next character and the process then repeats with the new batch of particles.

The main concern with using this method is the number of particles to use. Using more particles increases the processing necessary for estimating the distributions. However, a low number of particles could result in undersampling the distributions and missing important possible sequences. Sensitivity analysis is performed on the number of particles by determining the offline performance using 10, 100, 1000, 10 000 and 100 000 particles.