Wahlpflichtfach “Kognitions- und Klinisch-Psychologische Forschungsmethoden”
Wintersemester 2025/2026

Teil 1: Kognitive Modellierung und Diffusionsmodelle

Willkommen bei Stats by Randolph. Hier geht es zur Hauptseite mit weiteren Informationen und Inhalten.

Autor:innen dieser Seite: An den Inhalten dieser Seite haben Markus Janczyk und Valentin Koob mitgearbeitet. Der Inhalt dieses Textes wird in der Lehre in den Studiengängen Psychologie von der AG Psychologische Forschungsmethoden und Kognitive Psychologie an der Universität Bremen verwendet, steht aber allen Interessierten zur Verfügung. Rückmeldungen/Fehler/Vorschläge können gesendet werden an randolph@uni-bremen.

Versionshistory:

• v2.02: minimal korrigierte Version (24.10.2025)
• v2.01: minimal korrigierte Version (29.9.2025)
• v2.0: zweite, überarbeitete Version (7.9.2025)
• v1.0: erste Version (14.10.2024)

1 Einführung

Dieser Teil bietet eine kurze Einführung in die mathematische Modellierung kognitiver Prozesse mit Random Walk und Diffusionsmodellen. Detailierte und weiterreichende Informationen und Themen dazu, können unserem Online-Buch Modeling Experimental Data with (Time-Dependent) Diffusion Models entnommen werden.

1.1 Warum mathematische Modelle in der Psychologie?

Wir wollen in diesem Teil uns weniger mit statistischen Modellen, sondern mit Prozessmodellen beschäftigen, die Verhalten von Versuchspersonen vorhersagen. Wir fokussieren dabei auf Modelle zur Vorhersage von Reaktionszeiten (RTs) und Fehlerraten in einfachen RT-Experimenten, die grundsätzliche Idee lässt sich aber auf viele Anwendungsbereiche übertragen (z.B. Entscheidungspsychologie, Lernen von Kategorien, visuelle Suche, motorische Programmierung, Arbeitsgedächtnis, …). Wir wollen uns nun zunächst der Frage widmen, welchen Vorteil die Mathematisierung psychologischer Theorien denn überhaupt bietet und warum sich psychologische Theoriebildung solcher Methoden bedienen sollte und hierfür drei Gründe anführen.

1.1.1 Interpretation von Daten

Der vielleicht wichtigste Grund für formale, d.h. mathematische, Modelle in der Wissenschaft (und daher auch in der Psychologie) ist die Einsicht, dass Daten i.d.R. „nicht für sich selber sprechen“ und sich so einer adäquaten Interpretation entziehen. Dieser Gedanke ist unter statistikaffinen Personen eigentlich relativ anerkannt und klar, wie sich in diesem Zitat zeigt: “Although we often hear that data speak for themselves, their voices can be soft and sly” (Mosteller et al., 1983, S. 234).

Als ein illustratives Beispiel dafür werden oft die retrograden Bewegungen von Planeten (siehe hier für eine Abbildung) angeführt (siehe auch Lewandowsky & Farrell, 2013): Beobachtet gegen einen feststehenden Hintergrund der Sterne scheinen Planeten manchmal plötzlich ihre Richtung zu ändern, bevor sie dann wieder die ursprüngliche Bahn aufnehmen. Die reine Beobachtung, also die Daten, gibt es schon sehr lange. Dennoch haben diese Daten lange Zeit für Rätsel gesorgt, da sie ohne ein dahinterstehendes Modell nicht wirklich Sinn ergeben.

Tatsächlich gab es bereits mit dem geozentrischen Weltbild von Claudius Ptolemäus (verfasst in etwa 150 nach Christus) ein recht erfolgreiches Modell, welches die Daten gar nicht so schlecht erklären konnte: Während die Planeten um die Erde kreisen, rotieren sie auch noch um einen Kreis entlang ihres Orbits. Wenn dann die Erde noch leicht aus dem Mittelpunkt versetzt angenommen wird, können die Daten mit einer Diskrepanz von etwa 1° (für den Mars) beschrieben werden.

Warum hat sich dann aber das heliozentrische Weltbild von Kopernikus, also dass sich die Planeten inkl. der Erde um die Sonne drehen, so schnell und nachhaltig durchgesetzt? Ein Teil der Antwort ist, dass die Vorhersagequalität des Modells besser war, der Fehler also ein wenig kleiner als im ptolemäischen Weltbild. Der wichtigere Punkt ist die Einfachheit und Eleganz des neuen Modells, welches einfache kreisförmige Planetenbahnen um die Sonne zur Beschreibung benötigt (siehe auch hier für mehr Informationen dazu).

Quantitative Genauigkeit ist also ein wichtiger Faktor, Einfachheit ist aber ein anderer wichtiger Faktor, der über Erfolg und Misserfolg eines Modells entscheidet und immer berücksichtigt werden sollte. Das wichtigste ist aber hier: Manche Daten ergeben auf den ersten Blick, und für sich betrachtet, keinen Sinn; nur adäquate Modelle ermöglichen dann eine adäquate Interpretation. (Kepler hat übrigens die kreisförmigen Bahnen durch elliptische Bahnen ersetzt und damit einen Vorhersagefehler von nahezu Null erreicht.)

1.1.2 Kommunikation von Annahmen

Der zweite wichtige Punkt betrifft die Kommunikation von Annahmen. Die typische “verbale Theorie” bietet einen nicht unerheblichen Spielraum für die “Interpretationsfreiheit”. Ein einfaches Beispiel dafür ist es, wenn jemand sagt “Die Aktivation nimmt mit der Zeit ab.” Klar ist dann: Zu Beginn eines betrachteten Zeitraums sollte mehr Aktivation im (kognitiven) System vorhanden sein als zu einem späteren Zeitpunkt. Unklar ist allerdings der genaue Verlauf der Aktivationsreduktion. Abbildung 1.1 stellt drei mögliche Reduktionsverläufe dar (linear, exponentiell, logistisch), die allesamt mit der qualitativen Aussage konform sind. Allerdings unterscheiden sich die Verläufe doch sehr. Die entsprechenden Formeln lauten (\(t\) ist die Zeit und \(Y\) die Aktivierung):

• linear: \(Y = b \cdot t + a\)
• exponentiell: \(Y=b \cdot e^{k\cdot t}\)
• logistisch: \(Y=\frac{b}{1 + e^{k\cdot(t-x_0)}}\)

Abbildung 1.1: Verschiedene Varianten, wie Aktivation über die Zeit abnehmen kann.

Klar: Je präziser ein Modell formuliert wird, desto eher finden sich Daten, die gegen das Modell sprechen und es somit falsifizieren. Dies ist allerdings durchaus ein Qualitätsmerkmal von Modellen.

Wir werden gleich auch sehen, dass selbst sehr ähnlich ausschauende Funktionen durchaus verschiedene psychologische Implikationen haben können. An dieser Stelle halten wir erst einmal fest, dass eine Theorie umso “besser” ist, je präziser ihre Vorhersage ist. Sagt eine Theorie also nicht nur vorher “Die Aktivation nimmt mit der Zeit ab.”, sondern “Die Aktivation nimmt mit der Zeit einer Exponentialfunktion folgend ab.”, ist dies eine präzisere Vorhersage, die auch leichter falsifizierbar ist.

1.1.3 Fehlerhafte Intuition

Der dritte Punkt ist, dass unsere Intuition oft auch einfach falsch ist. Gerade bei komplexen Systemen, bei denen Prozesse voneinander abhängig sind, sind intuitive Vorhersagen manchmal wirklich schwer und oft sogar falsch. An dieser Stelle hilft es dann, die Systeme mathematisch zu formulieren und dann Vorhersagen durch Simulationen zu generieren. Tatsächlich ist der Zugang durch Simulationen eng verknüpft mit mathematischen Modellen und auch wir werden im Verlaufe dieses Themas Modellvorhersagen durch Simulationen gewinnen.

Das sog. Geburtstagsparadox ist ein bekanntes Beispiel aus der Wahrscheinlichkeitstheorie für problematische Intuitionen. Probieren Sie es zunächst selber und beantworten Sie folgende Frage: Wieviele Personen müssen sich in einem Raum befinden, damit die Wahrscheinlichkeit größer als 50% dafür ist, dass zwei Personen am gleichen Tag Geburtstag haben (das Jahr ignorieren wir dabei)?

Für dieses Problem gibt es zwar eine halbwegs einfache analytische Lösung; man kann aber auch sehr leicht in R ein Programm schreiben, welches die Situation für eine gegebene Anzahl von Personen, \(n\), simuliert und dann die jeweilige Wahrscheinlichkeit durch die relative Häufigkeit des Eintreffens des Ereignisses schätzt. Beide Varianten zur Lösung finden Sie in dieser Ergänzung.

1.2 Ein einfaches Beispiel: Die Lernfunktion

Ein erstes Beispiel, wo wir versuchen, den Zusammenhang zweier Variablen zu beschreiben, wenn dieser nicht linear ist (wie bei einer einfachen, linearen Regression), betrifft die typische Lernkurve. Wird die zur Bearbeitung einer Aufgabe benötigte Zeit (\(RT\)) als Funktion der Anzahl der Durchgänge die bearbeitet wurden (\(T\)) betrachtet, sieht man, dass diese Zeit immer kürzer wird. Häufig wird die Lernkurve durch eine Power-Funktion \[RT = a_P+b_P\cdot (T+1)^{-\beta}\] beschrieben (vgl. auch Rosenbaum & Janczyk, 2019). Heathcote et al. (2000) haben allerdings (auf Basis der Analyse individueller Daten anstelle von aggregierten Gruppendaten) argumentiert, dass der typische Verlauf auch/besser durch eine Exponentialfunktion \[RT = a_E+b_E\cdot e^{-\alpha\cdot T}\] beschrieben werden kann. Abbildung 1.2 zeigt mit einer Power-Funktion simulierte “empirische” Daten (die schwarzen Punkte) sowie “die beste” Power- sowie Exponentialfunktion zur Beschreibung dieser Daten (wie diese Funktionen gefunden werden, wird im Wahlpflichtfach in einer praktischen Übung erarbeitet).

Abbildung 1.2: Simulierte Daten einer Powerfunktion sowie gefittete Power- und Exponentialfunktionen.

Beide Funktionen passen tatsächlich relativ gut zu den empirischen Daten. Allerdings haben sie Detail unterschiedliche Implikationen für Lernprozesse und die Abnahme bzw. Konstanz der Lernrate über die Zeit. Dies können wir mit den einfachsten Formen einer Power- und einer Exponentialfunktion demonstrieren, indem wir das Verhältnis der vorhergesagten \(RT\) von Durchgang \(T\) zu \(T+1\) betrachten:

Trial <- c(1:10)
Power <- (Trial)^(-1) # einfache Power-Funktion
Exp <- exp(-Trial) # einfache Exponentialfunktion

# nun schauen wir die Lernrate an, also das Verhältnis der Funktionswerte
# eines Zeitpunkts T und des folgenden Zeitpunktes T+1
rateExp <- Exp[1:9] / Exp[2:10]
ratePower <- Power[1:9] / Power[2:10]
rateExp # bleibt konstant

## [1] 2.718282 2.718282 2.718282 2.718282 2.718282 2.718282 2.718282 2.718282
## [9] 2.718282

ratePower # wird kleiner

## [1] 2.000000 1.500000 1.333333 1.250000 1.200000 1.166667 1.142857 1.125000
## [9] 1.111111

Die Rate des Lernzuwachses bleibt bei einer Exponentialfunktion also konstant, d.h. man würde relativ gesehen mit jeder Übung mehr einen gleichen Anteil lernen und entsprechend Forschritte machen. Demgegenüber wird die Lernrate bei einer Powerfunktion mit der Zeit kleiner, man würde also mit der Zeit immer weniger lernen bis kein Fortschritt mehr erkennbar ist. Die Kenntnis der “richtigen” Funktion zur Beschreibung der Lernkurve, hat also Implikationen für Theorien des Lernens.

1.3 Arten und Korrektheit von Modellen

Wenn wir uns also darauf einigen, dass Modelle wichtig und hilfreich sind: Was ist eigentlich ein Modell und auf welche Arten von Modellen können wir uns beziehen? Ohne zu sehr ins Detail zu gehen, fassen wir ein Modell als eine vereinfachte, abstrakte Struktur auf, die darauf abzielt, etwas Komplizierteres zu beschreiben, vorherzusagen und/oder zu erklären (in unserem Fall ist dieses „Etwas“ die menschliche Kognition). Mit anderen Worten: Modelle “seek to retain essential features of the system while discarding unnecessary details“ (Lewandowsky & Farrell, 2013, S. 11).

Es lassen sich mehrere Unterscheidungen zwischen verschiedenen Klassen von Modellen treffen (siehe auch Luce, 1995, für eine ausführliche Darstellung) und wir besprechen hier nur kurz eine vereinfachte Version (wie sie von Lewandowsky & Farrell, 2013, vorgeschlagen wurde):

• Deskriptive Modelle beschreiben und fassen Daten zusammen, jedoch ohne notwendigerweise ein psychologisches Element. Beispiele hierfür sind die Power- oder Exponentialfunktion, wie im vorherigen Abschnitt besprochen: Beide Funktionen haben an sich keine psychologische Bedeutung. Dennoch können sie, wie wir gesehen haben, Implikationen für psychologische Theorien haben. Auch statistische Modell (wie z.B. die Regression, testtheoretische Messmodelle, Strukturgleichungsmodelle, …) sind in erster Linie deskriptive Modelle.

• Prozesscharakterisierende Modelle zielen darauf ab, die Zwischenschritte, die einem Verhalten zugrunde liegen, genauer zu erfassen. Das Ergebnis kann z.B. die Erkenntnis sein, dass bestimmte Prozesse unterschieden werden können/sollten und wie ein Prozess auf einem anderen aufbaut. Wie diese Prozesse genau funktionieren, ist dabei allerdings noch nicht Bestandteil des Modells.

• Prozesserklärende Modelle versuchen eine detaillierte Erklärung dafür zu liefern, wie ein Prozess tatsächlich funktioniert. Die Modelle, mit denen wir uns hier beschäftigen, Random Walks und Diffusionsmodelle, gehören zu dieser Klasse von Modellen.

Zuletzt noch eine Anmerkung: Vereinfachung in Modellen ist immer notwendig, einfach deshalb, weil wir das Modell am Ende verstehen wollen. Der Einsatz übermäßig komplexer Modelle untergräbt dieses Verständnis nur wieder. Vereinfachung bedeutet jedoch auch, dass (selbst gute) Modelle niemals wahr, sondern immer falsch sind (für einige Überlegungen dazu, was ein gutes Modell ausmacht, siehe Shiffrin et al., 2008; vgl. auch MacCallum, 2003). Das schließt jedoch nicht aus, dass es „bessere“ Modelle geben kann oder unplausible Modelle, die sogar verworfen werden können. Dennoch tragen gut formulierte, mathematische Modelle zur Weiterentwicklung des Faches bei, umso mehr, wenn sie unerwartete Vorhersagen machen. Mindestens aber verbinden sie Theorien mit Vorhersagen durch formale Herleitungen und/oder Simulationen und ermöglichen so eine strengere Überprüfung von Hypothesen auf Basis empirischer Daten (siehe auch Ulrich, 2009; Oberauer & Lewandowsky, 2019).

Prägnant werden Modelle manchmal auch mit dem Satz charakterisiert: “All models are wrong, but some are useful“ (siehe dazu auch Box, 1976).

1.4 Übersicht

Wir werden uns hier i.W. mit einem prozesserklärenden Modell befassen, welches die Reaktionswahlen und die dazugehörigen Reaktionszeiten (“response times”; RTs) erklären soll: das sog. Diffusionsmodell. Am Ende werden wir auch eine allgemeinere Form dieses Modells betrachten, welches verschiedene empirische Phänomene bei Konfliktaufgaben (z.B. Simon- oder Eriksen Flanker-Aufgaben) erklären kann.

Auf dem Weg dahin benötigen wir allerdings verschiedene Schritte. Zunächst zeichnen sich alle Modelle dadurch aus, dass sie eine Reihe von Parametern besitzen, deren Werte das genaue Verhalten des Modells festlegen. Wir suchen dann zunächst diejenigen Parameter, die das Modell so “ausrichten”, dass seine Vorhersagen möglichst nah an beobachtete, empirische Daten herankommen. Dieses Prinzip haben wir in Statistik I bei der einfachen, linearen Regression bereits eingeführt. Im Allgemeinen sind die “besten Werte” der Parameter allerdings oft nicht analytisch bestimmbar. In Kapitel 2 werden wir daher das grundlegende Prinzip numerischer Optimierung behandeln und dort auch Begriffe wie Nelder-Mead oder Differential Evolution Algorithmen und Maximum Likelihood Schätzung einführen. Schätzen wir Parameter auf Basis empirischer Daten, dann machen wir genau das, was gemeint ist, wenn jemand sagt “ein Modell sei an Daten gefittet worden”.

In Kapitel 3 führen wir dann mit Random Walks ein erstes stochastisches Modell zur Erzeugung von RTs und Reaktionswahlen ein und erweitern dieses dann in Kapitel 4 zum Diffusionsmodell, einem der am weitest verbreiteten Modelle in der psychologischen Forschung. Kapitel 5 schließlich widmet sich dann einem speziellen Diffusionsmodell, welches für die Analyse von Daten aus Konfliktaufgaben entwickelt wurde. An diesem werden wir dann auch eine Anwendung im klinisch-psychologischen Kontext betrachten.

2 Numerische Optimierung und Parameterschätzung

Ein wichtiger Teil zum Verständnis mathematischer Modellierung ist das Verständnis darüber, was es bedeutet “ein Modell an Daten zu fitten”. Anders formuliert ist dies die Frage: Wie müssen Parameter eines Modells gewählt werden, damit das Modell Vorhersagen macht, die möglichst gut mit der empirischen Realität übereinstimmen?

Grundsätzlich gibt es dazu zwei verschiedene Methoden:

1. Man kann das zugrunde liegende Problem analytisch lösen und so eine exakte Lösung für die besten Parameter angeben. Einen Fall, auf den dies zutrifft, haben wir mit der einfachen, linearen Regression bereits in Statistik 1 kennengelernt.
2. Man kann das zugrunde liegende Problem nicht analytisch lösen. In diesem Fall können die Parameter oft–zumindest annähernd–durch Computeralgorithmen bestimmt werden. Die ist z.B. bei der logistischen Regression der Fall gewesen, trifft aber auch auf andere statistische Verfahren zu, wie z.B. Mixed Effect Modelle.

Wir werden in diesem Kapitel das Beispiel der einfachen, linearen Regression aufgreifen und an ihr das grundsätzliche Vorgehen im zweiten Fall beschreiben. Wenngleich wir hier dann später auf (prozesserklärende) Diffusionsmodelle fokussieren werden sei angemerkt, dass die hier dargestellten Methoden und Verfahren in anderen Bereichen mathematischer bzw. statistischer Modellierung ebenso zum Einsatz kommen. Insbesondere die Maximum-Likelihood Schätzung (Abschnitt 2.6) findet auch zur Parameterbestimmung in Mess- und Strukturmodellen und damit auch der Faktorenanalyse und Linearen Strukturgleichungsmodellen ihre Anwendung. Insofern führen wir hier etwas sehr grundsätzliches, und nicht spezifisch für Diffusionsmodelle geltendes, ein.

2.1 Die einfache, lineare Regression als (statistisches) Modell

Die einfache, lineare Regression ist ein statistisches Modell welches beschreibt, dass zwischen zwei (hier intervallskalierten) Variablen ein linearer Zusammenhang besteht. Das heißt, es liegen von \(n\) Personen zwei Variablen \(X\) und \(Y\) vor, die üblicherweise als Scatterplot (Punktewolke) visualisiert werden. Dann wird diejenige Gerade gesucht, die die Daten bestmöglich beschreibt. Da Regressionen auch zur Vorhersage von Werten auf einer Variablen auf Basis (nicht erhobener Werte) auf der anderen Variablen verwendet werden, wird \(X\) auch als Prädiktor und \(Y\) als Kriterium bezeichnet, wenn \(X\) zur Vorhersage von \(Y\) genutzt werden soll. Die allgemeine Gleichung einer Regressionsgeraden lautet dann \[ \hat{Y}=bX+a\,. \] Hierbei sind dann:

• \(X\) der Prädiktor und \(\hat{Y}\) die von der Regressionsgeraden, d.h. vom statistischen Modell, vorhergesagten Werte des Kriteriums für die vorhandenen Prädiktorwerte.

• Ferner sind \(b\) und \(a\) die (zwei) unbekannten Parameter des Modells, wobei \(b\) die Steigung und \(a\) den Achsenabschnitt der Geraden meint. Sie sollen so gewählt werden, dass die Gerade möglichst gut zu den empirischen Daten passt.

“Möglichst gut zu empirischen Daten passen” die vom Modell vorhergesagten Daten dann, wenn die vorhergesagten und die empirischen Daten möglichst wenig voneinander abweichen. Man sagt auch, dass ihre Diskrepanz möglichst klein werden soll. Die Abweichung kann formalisiert werden als die Differenz \(E=Y-\hat{Y}\), die auch Residuum genannt wird. Die Idee ist nun, \(b\) und \(a\) so zu wählen, dass die Summe aller quadratischen Residuen insgesamt möglichst klein wird. Diese Summe bezeichnen wir mit \(Q\) und verstehen sie als Funktion der beiden Parameter \(b\) und \(a\). Nun suchen wir diejenigen Werte für \(b\) und \(a\), bei denen die Funktion \(Q\) \[\begin{align*} Q(b,a) &= \sum_{i=1}^n(y_i-\hat{y}_i)^2 \\ &= \sum_{i=1}^n(y_i-(bx_i+a))^2 \\ &= \sum_{i=1}^n (y_i-bx_i-a)^2 \end{align*}\] ihr Minimum annimmt. Dieses Vorgehen haben wir als die Methode der Kleinsten Quadrate kennengelernt (siehe hier für eine interaktive ShinyApp). Die analytische Möglichkeit zur Bestimmung der Werte für \(a\) und \(b\), die \(Q\) möglichst klein werden lassen, besteht in der Lösung eines Gleichungssystems partieller Differentialgleichungen. Es ergeben sich dann geschlossene Lösungen zur Berechnung von \(a\) und \(b\) mit beliebigen Daten \(X\) und \(Y\), nämlich \[b=\frac{\text{Kov}(X,Y)}{S_X^2}\text{ und }a=M_Y-bM_X\,.\] Das praktische Vorgehen bei einer Regression ist mit R sehr einfach. Zur kurzen Wiederholung schaffen wir uns dazu einen Beispieldatensatz, in welchem von 10 Versuchspersonen Werte auf den Variablen \(X\) (Prädiktor) und \(Y\) (Kriterium) vorhanden sind. Mit lm() berechnen wir dann eine einfache lineare Regression und geben uns die Parameter \(a\) und \(b\) mit coef() aus:

# Anlegen des Datensatzes
daten <- data.frame(
  vp = c(1:10),
  X = c(4, 6, 5, 9, 6, 5, 4, 10, 4, 7),
  Y = c(10, 9, 9, 11, 10, 10, 12, 11, 9, 13)
)
head(daten, 2)

reg.modell <- lm(Y ~ X, # Y als Kriterium, X als Prädiktor
  data = daten # DataFrame daten wird verwendet
)
coef(reg.modell) # Koeffizienten (Parameter) anzeigen

## (Intercept)           X 
##       9.050       0.225

Gibt es keine geschlossenen Gleichungen zur Bestimmung des Minimums einer Funktion müssen wir anders verfahren, um das Minimum von \(Q\) und damit die besten Parameter zu finden. Wie dies funktioniert, erläutern wir nun am gleichen Beispiel der einfachen, linearen Regression.

2.2 Numerische Optimierung

Natürlich könnten wir einfach für viele (alle?) mögliche Kombinationen von \(a\) und \(b\) den Wert von \(Q\) ausrechnen und dann das Minimum heraussuchen. Wenn dies systematisch durchgeführt wird, wird dies manchmal als Grid Search bezeichnet. Effizienter sind allerdings i.d.R. Verfahren aus der Numerischen Mathematik, die bestimmte Algorithmen benutzen, um ein Minimum zu suchen. Das praktische Vorgehen wird nun in den nächsten Abschnitten eingeführt.

Im Prinzip benötigen wir zwei Zutaten. Zunächst müssen wir die Abweichung (oder Diskrepanz) der vorhergesagten von den empirischen Daten durch eine einzige Zahl ausdrücken (“Kostenfunktion”; manchmal auch als Diskrepanzfunktion, objective function, … bezeichnet). Dann benötigen wir einen Algorithmus, der möglichst effizient diejenigen Parameter findet, für die die Kostenfunktion ihr Minimum annimmt. Dies sind dann automatisch diejenigen Parameter, bei denen die vorhergesagten und die empirischen Daten am wenigsten voneinander abweichen.

2.2.1 Definition einer Kostenfunktion

Die von uns bereits definierte Funktion \(Q\) ist tatsächlich bereits eine Kostenfunktion und wir verwenden sie an dieser Stelle daher auch zunächst weiter. Letztlich ist die Methode der Kleinsten Quadrate, die wir zur Definition von \(Q\) genutzt haben, aber nur eine Möglichkeit zur Formulierung einer Kostenfunktion.

Tatsächlich benötigen wir auch in R eine Funktion, die den Wert von \(Q\) auf Basis vorhergesagter (auf aktuell geschätzten Parameterwerten basierend) und empirischer Daten berechnet und an eine aufrufende Funktion zurückliefert. Bei den Argumenten und Parametern, die der Funktion übergeben werden, ist es wichtig, dass das erste Argument die Parameter enthält (das wird von der Funktion optim(), die wir gleich verwenden werden, vorausgesetzt). In der Kostenfunktion wird zunächst eine Funktion aufgerufen, die für aktuelle Parameterwerte \(a\) und \(b\) (gespeichert im Vektor parameter) die Modellvorhersagen generiert (also die \(\hat{Y}\)-Werte. Dann wird \(Q\) berechnet als quadrierte Differenz der vorhersagten und beobachteten Werte. Für unser Beispiel könnte dies dann so aussehen:

# Kostenfunktion: berechnet Q
Q_costFunction <- function(parameter, # Vektor mit den Parametern
                           obsData) { # empirische ("observed") Daten
  # Berechnung der vorhergesagten Werte für die aktuellen Parameter die in
  # parameter gespeichert sind
  predData <- predictRegression(
    parameter,
    obsData
  )

  # Berechnung von Q: quadrierte Residuen aufsummieren
  Q <- sum((obsData$Y - predData)^2)

  # Rückgabe von Q
  return(Q) # Rückgabe
}

2.2.2 Prinzip der numerischen Optimierung mit optim()

Das generelle Vorgehen bei einer Optimierung ist in Abbildung 2.1 dargestellt.

1. Den Ausgangspunkt bilden empirische Daten. Es sollen diejenigen Parameter gesucht werden, mit denen Modellvorhersagen generiert werden können, die möglichst nah an diese empirischen Daten herankommen, d.h. bei denen der Wert der Kostenfunktion möglichst klein wird. Diese Daten werden i.d.R. in das Skript in R geladen (auch wenn wir das Beispiel oben einfach direkt erzeugt haben).
2. Die empirischen Daten werden gemeinsam mit Startwerten für die gesuchten Parameter an die Funktion optim() übergeben. Zusätzlich wird der Funktion optim() der Name der Kostenfunktion übergeben, damit optim() weiß, welcher Wert minimiert werden soll. Die Funktion optim() stellt mehrere Algorithmen zur Verfügung, die das Minimum der Kostenfunktion finden können. Als Default wird der sog. Simplex-Algorithmus nach Nelder und Mead (1965) verwendet (siehe Abschnitt 2.5.1 für mehr Informationen).
3. Die Kostenfunktion ruft eine Funktion auf, die mit den aktuellen Parameterwerten die Modellvorhersagen berechnet. Im Beispiel haben wir bereits eine Funktion predictRegression() dafür eingebaut. Danach berechnet die Kostenfunktion den Wert von \(Q\) und gibt diesen an optim() zurück.
4. Ist der Wert von \(Q\) kleiner als in der letzten Iteration, passt optim() nach einem bestimmten Algorithmus die Parameterwerte an und der Ablauf beginnt von vorn.¹
5. Wenn der Wert von \(Q\) nicht mehr kleiner wird (bzw. genauer: sich nur noch in einem festgelegten Toleranzbereich bewegt), dann wird die Optimierung beendet und optim() springt in das aufrufende Programm zurück und liefert u.a. diejenigen Parameterwerte zurück, für die \(Q\) minimal war.

Im Detail und abhängig von den verwendeten Algorithmen, gibt es noch ein paar Punkte auf die es zu achten gilt (wir kommen darauf später noch einmal zurück). Wir bleiben aber nun bei der Funktion optim() und setzen das beschriebene Vorgehen in entsprechenden R Code um, um die praktische Durchführung kennen zu lernen.

Abbildung 2.1: Illustration des grundlegenden Ablaufs bei numerischer Parameterbestimmung.

2.2.3 Praktische Umsetzung in R

Die Kostenfunktion Q_costFunction() haben wir bereits weiter oben beschrieben. Als nächstes benötigen wir die in der Kostenfunktion aufgerufene Funktion predictRegression(), die mit einem Satz von Parametern \(a\) und \(b\) (die im Vektor parameter stehen) die für die empirischen Prädiktorwerte \(X\) vorhergesagten Werte \(\hat{Y}\) berechnet. Wir haben diese Funktion hier zum einfacheren Verständnis relativ ausführlich gehalten und extrahieren zunächst aus dem Parametervektor parameter Werte für \(a\) und \(b\) (wir einigen uns darauf, dass der erste Wert des Vektors \(b\) meint und der zweite Wert meint \(a\)). Danach wird dann mit diesen Parametern und den empirischen Werten des Prädiktors \(X\) eine Regression berechnet und die vorhergesagten Werte werden in Y_Dach gespeichert und zurückgegeben:

# Funktion zur Berechnung der vorhergesagten Werte mit aktuellen Parametern
predictRegression <- function(parameter,
                              obsData) {
  # wir extrahieren die beiden Parameter aus dem Vektor parameter
  b <- parameter[1] # extrahieren von b
  a <- parameter[2] # extrahieren von a

  # dann verwenden wir eine normale Regression und berechnen Y-Dach als
  # aktuelle Modellvorhersage
  Y_Dach <- b * obsData$X + a

  # der resultierende Vektor wird dann an die aufrufende Funktion (bei uns
  # die Kostenfunktion Q_costFunction) zurückgegeben
  return(Y_Dach)
}

Damit haben wir beide Bestandteile zusammen, die Modellvorhersagen generieren (predictRegression()) und das Modell an den empirischen Daten evaluieren (Q_costFunction()). Schließlich fehlt noch das zentrale Stück Code, in welchem die Funktion optim() aufgerufen wird. Wir benötigen dazu noch den Vektor mit den Startwerten für die Parameter (damit wird in der ersten Iteration eine Modellvorhersage generiert) und rufen dann optim() direkt auf und speichern das Ergebnis in result_regression:

# wir definieren zunächst einen Vektor mit Startwerten für die Parameter
parameter_start <- c(0, 0) # c(b, a)

# dann rufen wir optim() auf
result_regression <- optim(
  par = parameter_start, # Startwerte der Parameter
  fn = Q_costFunction, # Kostenfunktion: wird minimiert
  obsData = daten # empirische Daten
)

Die Ergebnisse sind nun in result_regression gespeichert und mit dem $-Operator können wir auf verschiedene Elemente zugreifen. Die wichtigsten sind dabei nun $par und $value:

# unter $value finden wir den kleinsten Wert der Kostenfunktion (also von Q)
# der von optim() gefunden werden konnte:
round(result_regression$value, 3)

## [1] 14.375

# unter $par finden wir dann schließlich den Vektor derjenigen Parameter,
# die zu diesem kleinsten Q geführt haben
round(result_regression$par, 3)

## [1] 0.225 9.049

Wie wir sehen, stimmen diese Werte ziemlich gut mit den analytisch bestimmten Werten von oben überein. Dass sie nicht ganz genau übereinstimmen liegt eben daran, dass hier keine analytische Bestimmung durchgeführt wurde.

Zur Illustration betrachten wir noch zwei weitere Abbildungen. Abbildung 2.2 zeigt die empirischen Daten als Scatterplot sowie die analytisch bestimmte Regressionsgerade in rot. Zusätzlich sind in (dunkler werdenden) Grautönen diejenigen 81 Regressionsgeraden eingezeichnet, die für die verschiedenen Zwischenwerte der Parameter auf dem Weg zur finalen Lösung entstanden sind.

Abbildung 2.2: Illustration verschiedener Regressionsgeraden (grau) auf dem Weg zur besten numerischen Lösung (rot).

Die Parameter auf dem Weg zur finalen Lösung und die jeweils resultierenden Werte für \(Q\) sind in Abbildung 2.3 als Funktion der Iteration dargestellt. Im linken Teil wird ersichtlich, dass beide Parameter zunächst eher grob variiert werden, um sich dann auf einen kleinen Wertebereich festzulegen. Im rechten Teil der Abbildung ist zu sehen, wie die Kostenfunktion zunächst stark reduziert wird und am Ende nicht mehr kleiner wird. Man kann aber auch sehen, dass ab und an Parameterwerte ausprobiert werden, die die Kostenfunktion wieder leicht größer werden lassen.

Abbildung 2.3: Entwicklung der Parameter und der Kostenfunktion als Funktion der Iterationen während der Optimierung.

2.3 Zusammenfassung des generellen Vorgehens

Wir haben am Beispiel der einfachen, linearen Regression nun kennengelernt, wie Modellparameter mit numerischen Methoden so bestimmt werden können, dass die Modellvorhersagen möglichst wenig von empirischen, beobachteten Daten abweichen. Diese Parameterbestimmung ist das, was mit “ein Modell an Daten fitten” gemeint ist.

Das Vorgehen dabei umfasste i.W. die folgenden Schritte:

1. Das (statistische) Modell macht für bestimmte Parameterwerte Vorhersagen. Im Fall der einfachen, linearen Regression wurden die \(\hat{Y}\)-Werte auf Basis (aktueller) Werte von \(b\) und \(a\) berechnet.
2. Die vorhergesagten Werte werden mit den empirischen Daten verglichen und die Abweichung wird in einer Zahl, hier \(Q\), ausgedrückt. Dies wird allgemeiner als Kostenfunktion bezeichnet.
3. Ein Optimierungsalgorithmus sucht nach bestimmten Prinzipien diejenigen Werte für die Parameter, für die die Kostenfunktion ihr Minimum annimmt.

In den folgenden Kapiteln werden wir uns dann weniger mit statistischen Modellen als eher mit prozesserklärenden Modellen befassen, die das Verhalten von Versuchspersonen bei bestimmten Aufgaben in Form von RTs und Reaktionswahlen vorhersagen. Das Prinzip wird dann aber sehr ähnlich sein:

1. Das Modell macht für bestimmte Parameterwerte Vorhersagen. Hier werden dann Verteilungen von RTs vorhergesagt.
2. Die vorhergesagten Verteilungen werden mit empirischen Verteilungen verglichen und die Abweichung wird in einer (Kosten-)Funktion ausgedrückt.
3. Ein Optimierungsalgorithmus sucht nach denjenigen Werten für die Modellparameter, für die die (Kosten-)Funktion ihr Minimum annimmt.

2.4 Varianten von Kostenfunktionen

Wir haben im Beispiel das Kriterium der Kleinsten Quadrate als Kostenfunktion \(Q\) genutzt, um im theoretischen Rahmen der analytischen Herleitung zu bleiben. Es gibt aber eine Vielzahl anderer möglicher Varianten, die Abweichungen der vorhergesagten von den beobachteten Daten auszudrücken.

Eine einfache und weitverbreitete Variante für kontinuierliche, intervallskalierte Daten, die wir auch in den folgenden Kapiteln zur weiteren Illustration verwenden werden, ist der sog. Root Mean Square Error oder kurz: RMSE. Haben wir \(J\)-viele Datenpunkte und meinen mit \(d_j\) einen empirischen Datenpunkt \(j\) und mit \(p_j\) den vom Modell vorhergesagten Datenpunkt, dann ist \[ RMSE=\sqrt{\frac{1}{J}\cdot \sum_{j=1}^J(d_j-p_j)^2}\,. \] Klar sollte sein, dass der Teil \(\sum_{j=1}^J(d_j-p_j)^2\) auch als die Summe der quadrierten Residuen verstanden werden kann, wenn wir im Kontext der Regression denken.

Eine weitere öfters verwendete Kostenfunktion beruht auf einer \(\chi^2\)-Statistik. Hierbei werden für bestimmte Wertebereiche (z.B Quantil-Bereiche von RTs) die beobachteten und die vom Modell vorhergesagten Häufigkeiten miteinander verrechnet (ganz so wie wir es in Statistik 1 für die \(\chi^2\)-Statistik kennengelernt haben). Je ähnlicher sich beide Häufigkeiten sind, desto kleiner wird der resultierende \(\chi^2\)-Wert.

Basierend auf einer angenommenen Verteilung kann auch die Likelihood berechnet und als Ausgangspunkt der Kostenfunktion genutzt werden. In vielen Fällen ist dies die “beste” Möglichkeit, allerdings ist die Berechnung der Likelihood komplexer (und manchmal nicht möglich). Abschnitt 2.6 widmet sich der Likelihood-Methode daher etwas mehr im Detail.

2.5 Optimierungs-Algorithmen

2.5.1 Simplex-Algorithmus nach Nelder-Mead

Wir haben bisher die Funktion optim() zur Bestimmung der Parameter verwendet. Unterschlagen haben wir dabei, dass die Funktion mit der Defaulteinstellung method = "Nelder-Mead" arbeitet. Damit ist ein sog. Simplex-Algorithmus nach Nelder und Mead (1965) gemeint, der einen “allgemein-verwendbaren Algorithmus” darstellt, der auch in vielen Arbeiten aus der Psychologie verwendet wird. Im Folgenden wollen wir die Funktionsweise des Simplex-Algorithmus anhand einer einfachen linearen Regression darstellen.

Wie funktioniert der Simplex-Algorithmus nun? Zunächst ist ein Simplex eine geometrische Figur, die aus \(p+1\) Punkten in einem \(p\)-dimensionalen Parameterraum besteht. Wenn, wie im Falle der einfachen linearen Regression, \(p=2\) ist, ist ein Simplex ein Dreieck. Abbildung 2.4 illustriert die Kostenfunktion \(Q\) (d.h., die Sum of Squares) in Abhängigkeit von \(b\) und \(a\) und in Rot einzeichnet ist der initiale Simplex. Dieser ergibt sich wie folgt: Eine Ecke des Dreiecks ergibt sich aus den Startwerten, die der Funktion optim() übergeben werden. In den ersten \(p=2\) Schritten verändert der Algorithmus dann jeweils einen der Parameter um einen bestimmten Wert und generiert damit die beiden anderen Ecken des Dreiecks.

Abbildung 2.4: Illustration eines Simplex-Durchlaufs auf dem Weg ins globale Minimum im Falle einer Regression.

Nun kann der Wert der Kostenfunktion \(Q\) für jede Ecke des Dreiecks bestimmt werden. Der Algorithmus funktioniert dann so, dass der Simplex im Parameterraum so verschoben und umgeformt wird, damit in jedem Schritt der schlechteste Wert der Kostenfunktion zu verbessern versucht wird. Dabei sind die wichtigsten Operationen die folgenden:

1. Reflektion/Spiegelung: Der schlechteste Wert der Kostenfunktion wird durch den Zentroiden des Simplex reflektiert.

2. Der neue Wert der Kostenfunktion wird dann evaluiert und mit den anderen Werten verglichen:

   • Expansion/Streckung: Ist der neue Wert besonders gut, da er ein neues Minimum liefert, wird der Simplex in diese Richtung weiter expandiert, d.h. gestreckt.
   • Akzeptanz: Ist der neue Punkt weder der niedrigste noch der höchste Wert der Kostenfunktion, dann wird er so belassen, also akzeptiert.
   • Kontraktion: Ist der neue Wert besonders schlecht, da er ein neues Maximum liefert, wird der Simplex zum Zentroiden hin kontraktiert.
   • Zusammenziehen: Wenn keiner der vorangegangenen Schritte zu einer Verbesserung führt, wird der Simplex zu dessen bestem Punkt zusammengezogen.

3. Wenn alle Ecken des Simplex ähnliche Werte auf der Kostenfunktion liefern, dann stoppt der Algorithmus. Ansonsten wird mit Schritt 1 weitergemacht.

In Abbildung 2.4 ist in blau ein exemplarischer Pfad auf dem Weg zum (globalen) Minimum (grüner Punkt) eingezeichnet. So einfach und gut der Simplex-Algorithmus im Prinzip funktioniert, er hat auch eine Reihe von Problemen, die gerade bei komplexen Modellen seine Funktion erschweren können. Zwei dieser Probleme erwähnen wir hier daher direkt:

1. Der Wertebereich der Parameter ist in der Standardversion erstmal nicht eingeschränkt. Wir kommen aber auch in Situationen, in denen z.B. Wahrscheinlichkeiten geschätzt werden müssen und für diese gilt nun einmal \(p\in [0,1]\). In solchen Fällen ist es wenig hilfreich, wenn optim() mit Werten außerhalb dieses Intervalls versucht zu arbeiten. In anderen Situationen weiß man vielleicht vorab schon, in welchem Wertebereich ein Parameter plausibel wäre. In solchen Fällen können wir z.B. die Funktion nmkb() aus dem Paket dfoptim nutzen. Sie funktioniert quasi wie optim(), nutzt auch den Simplex-Algorithmus, kann aber mit lower und upper auch untere und obere Grenzen für die Parameter entgegen nehmen:

# Angabe der Startwerte sowie der Wertebereiche, innerhalb derer
# die Parameterwerte bleiben sollen während der Optimierung
parameter_start <- c(0.6, 30) # Startwerte
parameter_lower <- c(0.5, 1) # untere Grenze
parameter_upper <- c(1, 60) # obere Grenze

# Aufruf ähnlich wie optim()
result <- nmkb(
  par = start.pars, # Startwerte
  lower = lower.pars, # untere Grenze
  upper = upper.pars, # obere Grenze
  fn = Q, # Kostenfunktion
  emp.data = ... # empirische Daten die gefittet werden
)

2. Der Simplex-Algorithmus hat eine gewisse Neigung dazu, in lokalen Minima “hängenzubleiben” und damit das globale Minimum nicht zu finden. Im obigen Beispiel mit zwei Parametern einer einfachen, linearen Regression tritt dies wegen der gut strukturierten Oberfläche von \(Q\) zwar nicht auf, Modelle die wir später behandeln werden, sind aber oft nicht so einfach strukturiert. Hier gibt dann z.B. die folgenden zwei Auswege. (1) Wird mit dem Simplex-Algorithmus gearbeitet, wird die Optimierung oft mehrere Male mit zufällig gezogenen Startwerten wiederholt. Kommt der Algorithmus immer wieder zu ähnlichen Ergebnissen, liegt es zumindest nahe, dass dann tatsächlich das globale Minimum gefunden wurde. (2) Andere Algorithmen haben generell weniger Probleme mit lokalen Minima, funktionieren aber auch grundsätzlich anders und sind mitunter auch etwas rechenintensiver. Ein häufig verwendetes Beispiel ist der Differential Evolution Algorithmus.

2.5.2 Differential Evolution Algorithmus

Differential Evolution (DE) ist ein populationsbasierter Optimierungsalgorithmus, der von Storn und Price (1997) eingeführt wurde. Der Algorithmus beginnt mit der zufälligen Erzeugung von \(NP\) Parameterkombinationen, die typischerweise gleichmäßig über den Parameterraum verteilt sind, dessen untere und obere Grenzen für jeden Parameter dem Algorithmus vorgegeben werden (was auch das gerade genannte Problem 1 eines einfachen Simplex-Algorithmus mitberücksichtigt). Diese Parameterkombinationen werden typischerweise als jeweils ein Vektor \(\theta\) zusammengefasst und bilden die Populationselemente der Anfangsgeneration.

In jeder Iteration erzeugt DE eine neue Generation, indem bestehende Mitglieder rekombiniert und mutiert werden. Aufgrund dieser biologischen Metapher wird DE manchmal als genetischer Algorithmus bezeichnet und seine Iterationen werden in Begriffen wie „Generationen“, „Mutation“, „Kinder“ usw. beschrieben.

Eine gängige Strategie zur Aktualisierung der Population ist wie folgt: Für jeden aktuellen Parametervektor \(\theta_{i,g}\) in Generation \(g\) (wobei \(i \in \{1, \ldots, NP\}\)), wählt der Algorithmus drei verschiedene andere Mitglieder der Population aus: \(x_{r_0,g}\), \(x_{r_1,g}\) und \(x_{r_2,g}\). Dann wird ein skalierter Differenzvektor zwischen zwei von ihnen zum dritten addiert, um einen Vorschlagsvektor (eine „Mutation“) \[\begin{equation} \theta_i^* = x_{r_0,g} + F \cdot (x_{r_1,g} - x_{r_2,g}) \tag{2.1} \end{equation}\] zu bilden. Dieser Vorschlag wird anschließend vor dem Hintergrund der Kostenfunktion bewertet. Falls der Vorschlag eine bessere Anpassung als das aktuelle Mitglied \(\theta_{i,g}\) liefert (also einen kleineren Wert), ersetzt er dieses in der nächsten Generation. Dieser Prozess wird für alle Mitglieder der aktuellen Generation \(g\) wiederholt, wodurch eine neue Generation \(g+1\) entsteht. Der Faktor \(F\) steuert das Ausmaß der Mutation und wird typischerweise auf einen Wert um \(0.8\) gesetzt (wobei er zwischen 0 und 2 liegen muss).

Abbildung 2.5 veranschaulicht die Grundidee anhand derselben Oberfläche von \(Q\), die weiter oben für den Simplex-Algorithmus bereits benuzt wurde. Anders als die vorherige 3D-Darstellung zeigt diese aber eine 2D-Konturkarte, was die Visualisierung erleichtert. Die schwarzen Kreise repräsentieren die aktuellen Mitglieder der Population, nummeriert von 1-20 (zur besseren Beschreibung). Der blaue Punkt markiert ein ausgewähltes Zielmitglied \(\theta_{1,g}\), welches eine Mutation gemäß Gleichung (2.1) durchläuft.

Abbildung 2.5: Illustration einer Mutation im Rahmen eines Differential Evolution Algorithmus.

Zwei Mutationsvorschläge für das erste Mitglied werden dargestellt, einer davon erfolgreich, der andere nicht:

• Im ersten Fall wählt der Algorithmus zufällig \(x_{19,g}\), \(x_{3,g}\) und \(x_{9,g}\) (grün hervorgehoben). Er berechnet die Differenz zwischen den Mitgliedern 3 und 9 (grün gestrichelter Pfeil), skaliert diese mit \(F\) und addiert sie zu Mitglied 19 (grüner Pfeil). Das Ergebnis ist der Mutationsvorschlag \(\theta_1^*\), dargestellt als orangefarbener Punkt mit der Bezeichnung „1*“. Da dieser Vorschlag einen niedrigeren Wert auf der Kostenfunktion liefert, wird er akzeptiert und \(\theta_{1,g}\) wird in der nächsten Generation \(g+1\) durch \(\theta_1^*\) ersetzt. In der Sprache genetischer Algorithmen könnte man sagen, dass das „Kind“ fitter ist als der „Elternteil“.

• Im zweiten Fall wählt der Algorithmus \(x_{14,g}\), \(x_{8,g}\) und \(x_{5,g}\) (rot hervorgehoben). Die Mutation verläuft nach demselben Prinzip, aber der resultierende Vorschlag hat einen höheren Wert auf der Kostenfunktion. Daher wird der Vorschlag abgelehnt und \(\theta_{1,g}\) bleibt für die nächste Generation \(g+1\) unverändert.

Die hier dargestellte Version eines DE Algorithmus bildet quasi die Basis; letzlich handelt es sich um eine Familie von Algorithmen mit verschiedenen Verfeinerungen, auf die wir hier nicht eingehen. Eine Implementation für R liegt im Paket DEoptim mit der gleichnamigen Funktion DEoptim() vor. Der generelle Aufruf der Funktion ist ähnlich wie in den vorherigen Beispielen und wir werden in der Praxis auch mit dieser Funktion später arbeiten.

2.5.3 Performanzvergleich: Simplex vs. Differential Evolution

Um die Unterschiede in der Robustheit gegenüber lokalen Minima zwischen DE und Simplex zu veranschaulichen, betrachten wir die Rastrigin-Funktion, eine mathematische Funktion, die speziell dafür entworfen wurde, Optimierungsalgorithmen herauszufordern. Abbildung 2.6 veranschaulicht, warum diese Funktion schwer zu optimieren ist: Sie besitzt viele lokale Minima, in denen Optimierungsalgorithmen leicht stecken bleiben und dadurch suboptimale Lösungen finden können. Nichtsdestotrotz hat diese Funktion bei \(x=y=0\) ein wohldefiniertes globales Minimum.

Abbildung 2.6: Beispiel einer Rastrigin-Funktion zur Performanzprüfung von Optimierungsalgrithmen.

Nun setzen wir einen DE Algorithmus darauf an, das globale Minimu zu finden. Abbildung 2.7 zeigt mehrere Generationen des Durchlaufs, um dessen Prinzip zu veranschaulichen. Hierfür verwenden wir eine 2D-Projektion der Rastrigin-Funktion, um zu sehen, wie sich die Population (in Blau dargestellt) über die Generationen hinweg entwickelt. Innerhalb jeder Generation wird die Parameterkombination mit dem niedrigsten Wert der Kostenfunktion in Orange hervorgehoben.

Wir sehen, dass sich die gesamte Population bereits in den ersten 30 Generationen schnell in Richtung des globalen Minimums bewegt (die oberen vier Teilabbildungen). Danach sammeln sich die meisten Mitglieder der Population entweder direkt im globalen Minimum oder in benachbarten lokalen Minima (untere linke Teilabbildung). Schließlich haben sich bis zur 50. Generation alle Mitglieder stark zusammengezogen und sich um das globale Minimum herum stabilisiert (untere rechte Teilabbildung).

Abbildung 2.7: Beispielhafte Generation eines Differential Evolution Algorithmus bei der Suche nach dem globalen Minimum der Rastrigin-Funktion.

Zum Vergleich nun, wie sich ein Nelder-Mead Algorithmus verhält. Wie Abbildung 2.8 entnommen werden kann, hat dieser offensichtlich Schwierigkeiten, das Minimum zu finden. Die Abbildung zeigt die Suchlinie des Nelder-Mead-Verlaufs, wobei der Startpunkt in Blau und der Konvergenzpunkt in Orange dargestellt ist. Tatsächlich ist es so, dass selbst bei einem Startpunkt in der Nähe des globalen Minimums, die Funktion nmkb() dieses nicht gefunden hat, und sogar schlechter abschneidet als es der Startpunkt war. Es ist deutlich zu erkennen, dass Nelder-Mead die Oberfläche nicht besonders gründlich erkundet hat. Stattdessen entfernte sich der Algorithmus vom globalen Minimum und landete in einem weit entfernten lokalen Minimum.

Abbildung 2.8: Beispielhafter Verlauf eines Nelder-Mead Simplex Algorithmus bei der Suche nach dem globalen Minimum der Rastrigin-Funktion.

Mehr Details und Informationen (auch zu Bayesianischer Parameterschätzung) können hier gefunden werden.

2.6 Maximum-Likelihood Schätzung

Die Maximum-Likelihood Methode spielt eine zentrale Rolle in der statistischen Modellierung und hat eine lange Tradition sowohl im frequentistischen als auch im bayesianischen Rahmen. Grob gesagt versucht man bei der Maximum-Likelihood-Schätzung (MLE; von engl. Maximum-Likelihood Estimation), jene Parameterwerte zu finden, unter denen die beobachteten Daten am plausibelsten (“wahrscheinlichsten”) sind. MLE ist tief in der statistischen Theorie begründet und liefert eine Reihe wünschenswerter Eigenschaften für die geschätzten Parameter, oft unabhängig vom konkreten Modell. Sofern bestimmte Regularitätsbedingungen erfüllt sind (siehe z.B. Henze, 2024), sind die Schätzer…

• konsistent: Die Schätzwerte werden mit zunehmender Beobachtungszahl immer genauer.
• asymptotisch normalverteilt: Ihre Verteilung nähert sich mit wachsender Stichprobengröße einer Normalverteilung an.
• effizient: Sie erreichen die geringstmögliche Varianz unter den unverzerrten Schätzern.

Darüber hinaus bildet die Likelihood-Theorie eine Grundlage zur Beurteilung der Modellgüte und zum Vergleich alternativer Modelle anhand etablierter Kriterien wie dem Likelihood-Ratio-Test, dem AIC und dem BIC. Insgesamt stellt die Maximum-Likelihood Methode ein leistungsfähiges und fundiertes Verfahren zur Parameterschätzung dar, welches in vielfältigsten Kontexten zur Anwendung kommt. Eine generelle Einführung in die MLE (am Beispiel der Schätzung der Parameter einer Normalverteilung) findet sich hier und in einer passenden ShinyApp.

Um mit MLE zu arbeiten, benötigen wir (bei stetigen Variablen) eine Dichtefunktion, die probability density function (PDF). In unserem Kontext werden wir uns mit Modellen für RTs befassen und daher wird die PDF von RTs benötigt. Tatsächlich sagen die Modelle genau diese PDF vorher, wir werden uns aber im nächsten Abschnitt generell erst einmal damit befassen, wie RTs verteilt sind und was für statistische Modelle für ihre Verteilung auch herangezogen werden.

2.6.1 Die Verteilung von Reaktionszeiten

2.6.1.1 Daten eines Simon-Experiments

Wir verwenden hier Daten von \(n=32\) Versuchspersonen, die eine Simon-Aufgabe bearbeitet haben (vgl. auch den Teil des Expras im Wintersemester; die Daten können hier heruntergeladen werden). Grob gesagt, sollte auf die Buchstaben X bzw. H mit einem linken bzw. rechten Tastendruck reagiert werden, wobei der Buchstabe jeweils auf der linken oder rechten Seite des Bildschirms erscheinen kann. Erscheint er auf der Seite der korrekten Reaktion, redet man von einem kongruenten Trial, ansonsten von einem inkongruenten Trial (synomym werden dafür die Begriffe kompatibel und inkompatibel verwendet). Typischerweise sind die RTs in kongruenten Trials kürzer als in inkongruenten Trials. Hier würde man dann von einem Simon-Effekt sprechen, der anzeigt, dass die Versuchspersonen die Lokation nicht ignorieren können, auch wenn sie für die Bearbeitung der Aufgabe keine Rolle spielt.

simon_daten <- read.table("./Daten/daten_Simon.txt",
  header = TRUE,
  sep = ";"
)
head(simon_daten, 2)

length(unique(simon_daten$vpid))

## [1] 32

Die einzelnen Variablen bedeuten dabei:

• vpid: Versuchspersonen-Nummer zur Identifikation zusammengehörender Daten
• SR: Variante der Ausbalancierung “Stimulus zu Reaktion”
• blocknr: Nummer des aktuellen Blocks, 0 ist dabei der Übungsblock
• trialnr: Nummer des Trials im aktuellen Block
• stimulus: Stimulus des aktuellen Trials
• req_response: erforderte (d.h. korrekte) Reaktion
• location: Lokation des Stimulus im aktuellen Trial
• cong: Kongruenz (d.h. Übereinstimmung von Stimuluslokation und erforderter Reaktion)
• response: gegebene Reaktion
• rt: Reaktionszeit (in Millsekunden)
• error: 0 = korrekte Reaktion, 1= Fehler
• errortype: 0 = alles richtig im Trial, sonst: zu schnell reagiert, zu langsam reagiert, …. (also generelle Fehler)

Wir werten die RTs nun einmal mit den Standardmethoden aus, um das Vorhandensein eines Simon-Effektes zu evaluieren:

# Filtern der Daten: Auschließen der Übungsblöcke, fehlerhafter Trials,
# und zu kurzer/langer RTs ("Ausreißer")
f_simon_daten <- subset(
  simon_daten,
  blocknr > 0 & error == 0 & errortype == 0 &
    rt >= 200 & rt <= 1200
)

# faktorisieren der Daten
f_simon_daten$vpid <- as.factor(f_simon_daten$vpid)
f_simon_daten$cong <- as.factor(f_simon_daten$cong)

# die gemittelten RTs sind kürzer in kongruenten als in inkongruenten
# Durchgängen = Simon-Effekt
ezStats(
  data = f_simon_daten,
  dv = rt,
  wid = vpid,
  within = cong
)

## Warning: Collapsing data to cell means. *IF* the requested effects are a subset
## of the full design, you must use the "within_full" argument, else results may
## be inaccurate.

# Varianzanalyse mit einem 2-stufigem within-subject Faktor "cong":
# der RT Unterschied ist signifikant
ez_result <- ezANOVA(
  data = f_simon_daten,
  dv = rt,
  wid = vpid,
  within = cong,
  detailed = TRUE
)

## Warning: Collapsing data to cell means. *IF* the requested effects are a subset
## of the full design, you must use the "within_full" argument, else results may
## be inaccurate.

anova_out(ez_result)

## $`--- ANOVA RESULTS     ------------------------------------`
##        Effect       MSE df1 df2       F     p petasq getasq
## 1 (Intercept) 4627.6835   1  31 3797.48 0.000   0.99   0.99
## 2        cong  362.4898   1  31    9.99 0.004   0.24   0.02
## 
## $`--- SPHERICITY TESTS  ------------------------------------`
## [1] "N/A"
## 
## $`--- FORMATTED RESULTS ------------------------------------`
##        Effect                                    Text
## 1 (Intercept) F(1, 31) = 3797.48, p < .001, np2 = .99
## 2        cong F(1, 31) =    9.99, p = .004, np2 = .24
## 
## $`NOTE:`
## [1] "Reporting unadjusted p-values."

2.6.1.2 Verteilungsmodelle

In Experimenten wie dem Simon-Experiment von gerade liegen i.d.R. etliche RTs einer Versuchsperson für jede Bedingung vor. Eine wichtige Frage ist nun: Welcher Verteilung folgen die RTs, wenn man sie als Zufallsvariable auffassen würde?

Naheliegend, weil so oft angenommen, wäre eine Normalverteilung. Allerdings ist auch klar, dass (a) RTs\(\leq 0\) nicht möglich sind und, dass (b) man zwar nicht beliebig schnell reagieren kann, aber prinzipiell beliebig langsam. Wir visualisieren in Abbildung 2.9 daher einmal als Histogramm die RTs aller Versuchspersonen (links) sowie exemplarisch die RTs der Versuchsperson 5 (rechts):

Abbildung 2.9: Histogramme von RTs von Versuchspersonen in einem Simon-Experiment.

Wir können sehen, dass RTs nicht symmetrisch, sondern linkssteil bzw. rechtsschief verteilt sind: Ein etwas größerer Anteil der RTs liegt im “kurzen” Bereich, während der Anteil hin zu den längeren RTs immer kleiner wird. Zur formalen Beschreibung von RTs wurden verschiedene Verteilungen vorgeschlagen, u.a. die Wald-Verteilung (auch inverse Gauß-Verteilung genannt) oder die exGauß-Verteilung:

• Die Wald-Verteilung ist eine Dichtefunktion, die auf \([0,\infty]\) definiert ist und zwei Parameter hat: \(\mu>0\) ist der Erwartungswert und \(\lambda>0\) ein Skalierungsparameter.
• Die exGauß-Verteilung ist eine Faltung (“Addition”) einer normalverteilten und einer exponentiell-verteilten Zufallsvariablen mit drei Parametern: \(\mu\) und \(\sigma\) sind der Erwartungswert und die Standardabweichung der Normalverteilung, \(\lambda\) der Parameter der Exponentialverteilung. Insgesamt beschreibt \(\mu\) damit die Lokation, \(\sigma\) die Variabilität und \(\lambda\) die Schiefe der Verteilung.

Mitunter werden die Parameter dieser Verteilungen auch psychologisch interpretiert, was aber nicht unumstritten ist (siehe Matzke & Wagenmakers, 2009). Abbildung 2.10 illustriert beide Verteilungen für ausgewählte Parameterwerte:

Abbildung 2.10: Illustration von Wald-Verteilungen (links) und exGauß-Verteilungen (rechts).

Betrachten wir die Wald-Verteilung, wird klar, dass wir sie entlang der x-Achse noch um eine bestimmte Konstante “nach rechts verschieben” müssen, um in den Bereich realistischer RTs zu gelangen. Tut man dies, spricht man auch von einer “shifted Wald-Verteilung”. Auf die Wald-Verteilung werden wir später noch einmal kurz zurückkommen. Im Folgenden gehen wir zunächst davon aus, RTs folgen einer exGauß-Verteilung.

2.6.2 Likelihood-Funktion

Angenommen, die rote Kurve in Abbildung 2.11 visualiert eine (von einem Modell vorhergesagte) Verteilung von RTs. Wir sehen, dass die RTs zu einem großen Teil zwischen 300 und 1000 ms liegen. Wichtig ist dabei: Da wir es mit einer PDF zu tun haben, stellen die Werte auf der y-Achse keine direkten Wahrscheinlichkeiten dar, sondern Dichten. Dies liegt daran, dass Zeit eine kontinuierliche Variable ist, was bedeutet, dass es eine unendliche Anzahl möglicher RTs gibt, die theoretisch beobachtet werden könnten. Deshalb ergibt es keinen Sinn, einer einzelnen Zeit eine bedeutungsvolle Wahrscheinlichkeit zuzuweisen (die Wahrscheinlichkeit wäre 0). Stattdessen müssen Wahrscheinlichkeiten als Flächen unter der Kurve berechnet werden, also durch Integration der Dichte über ein Intervall.

Dennoch besagen die Werte auf der PDF, welche Werte “plausibler” sind als andere, und dies bringt uns dem Konzept der Likelihood bereits näher. Wichtig ist zu sehen, dass jede PDF bedingt ist, da sie von den Parametern des Modells (und tatsächlich auch vom Modell selbst) abhängt. Im Falle der exGauß-Verteilung also von \(\mu\), \(\sigma\) und \(\lambda\). Verändern sich die Parameter, verändert sich Form und Lage der PDF, und damit verändert sich auch der Dichtewert für eine fixe/gegebene RT. Mit \(\theta\) bezeichnet man im Allgemeinen einen Vektor der Parameter, hier also \(\theta=(\mu,\sigma,\lambda)\). Mit \[ p(y|\theta) \] meinen wir dann die PDF des Modells, wobei \(y\) ein Platzhalter für beobachtete RTs ist. Die Schreibweise betont allerdings, dass der Wert auf der PDF bedingt ist auf einen bestimmten Parametersatz \(\theta\).

Abbildung 2.11 illustriert für zwei verschiedene Parametervektoren, \(\theta_{rot}=(400, 50, 200)\) und \(\theta_{blau}=(500, 60, 150)\), die resultierenden PDFs. Zusätzlich sind als roter und blauer Punkt die PDF-Werte für eine exemplarische RT von 450 ms eingezeichnet. Ersichtlich wird, dass diese RT unter dem Modell mit den Parametern \(\theta_{rot}\) plausibler ist, als unter dem Modell mit den Parametern \(\theta_{blau}\).

Abbildung 2.11: Illustration zweier PDFs und der resultierenden Likelihood-Werte für eine fixierte RT.

Bei der Parameterschätzung mittels MLE möchten wir nun diejenigen Parameterwerte, also ein bestimmtes \(\theta\), finden, mit denen das Modell die Daten am besten beschreibt. Im Kontext von PDFs stellt sich damit die Frage: Welche Parameterwerte \(\theta\) liefern den höchsten Dichtewert, gegeben die Daten? Anders ausgedrückt: Wie müssen wir die PDF des Modells in Abhängigkeit der Parameter in \(\theta\) gestalten, um die maximalen Wahrscheinlichkeitsdichtewerte für die vorliegenden Daten zu erhalten?

Die Funktion, die uns dabei hilft, wird Likelihood-Funktion \[ \mathcal{L}(\theta|y) = p(y|\theta) \] genannt. Sie beschreibt, wie plausibel es ist, dass ein bestimmter Datensatz \(y\) von einem Modell mit dem entsprechenden Parametersatz \(\theta\) erzeugt wurde. Wichtig ist dabei: Die Likelihood-Funktion führt nichts Neues ein, sie verweist letztlich zurück auf die PDF. Sie zeigt lediglich auf, wie sich die PDF-Werte für einen festen Datensatz \(y\) in Abhängigkeit von \(\theta\) verändern.

Nun liegen i.d.R. mehrere RTs vor, statt nur eine wie im Beispiel von gerade. Allerdings könnnen wir, unter Annahme der Unabhängigkeit, die gemeinsame Likelihood als Produkt der einzelnen Likelihoods berechnen. Liegen als RTs \[ y=(450, 600, 800) \] vor, dann würden die einzelnen Likelihoods sich wie in Abbildung 2.12 bestimmen.

Abbildung 2.12: Illustration der Likelihood-Werte bei drei fixierten RTs.

Die Gesamt-Likelihood ergibt sich dann als \[\begin{align*} \mathcal{L}(\theta|y)&=p(y_1|\theta)\cdot p(y_2|\theta)\cdot p(y_3|\theta)\\ &= 0.0031\cdot 0.0019\cdot 0.0007\\ &= 4.123\cdot 10^{-09}\,. \end{align*}\]

Variieren wir nun z.B. systematisch den Parameter \(\mu\) in \(\theta\) und berechnen dann jeweils für die drei RTs die gemeinsame Likelihood, erhalten wir eine Funktion der Likelihood in Abhängigkeit von \(\mu\), wie Abbildung 2.13 dargestellt. Bei \(\mu=434\) erreicht diese Funktion ihr Maximum und der Wert wäre daher der Maximum-Likelihood Schätzer für \(\mu\); also derjenige Wert für \(\mu\) unter dem die Gesamtdaten am plausibelsten sind.

Abbildung 2.13: Likelihood in Abhängigkeit von \(\mu\).

Natürlich liegen in der Regel mehrere Datenpunkte in \(y\) vor. Bei \(n\) Datenpunkten (also z.B. \(n\)-vielen RTs) ergibt sich ganz allgemein dann die Likelihood-Funktion als \[\begin{equation} \mathcal{L}(\theta|y)=\prod_{i=1}^np(y_i|\theta)\,. \tag{2.2} \end{equation}\] In der Realität wird diese Schätzung simultan natürlich alle Parameter in \(\theta\) berücksichtigen. Klar ist aber auch hier: Es wird ein Extremwert einer Funktion gesucht und dies können am Ende die weiter oben erwähnten Algorithmen wie Simplex oder Differential Evolution für uns tun. Allerdings führen wir im nächsten Abschnitt noch einen technischen Kniff ein, der die Suche nach dem Extremwert erleichtert.

2.6.3 (negative) log-Likelihood

Die in Gleichung (2.2) dargestellte Likelihood-Funktion ist prinzipiell zur Schäzung der Parameter in \(\theta\) geeignet. Allerdings beinhaltet die Likelihood-Funktion das Produkt vieler PDF-Werte, und genau dies kann für einen Computer problematisch werden: Likelihood-Werte sind fast immer entweder sehr groß oder sehr klein, so extrem, dass Computer Schwierigkeiten haben können, sie genau darzustellen. In Abbildung 2.13 z.B. sind die resultierenden Werte extrem klein (siehe die Skalierung der y-Achse).

Daher betrachtet man i.d.R. den Logarithmus der Likelihood-Funktion, also die sog. log-Likelihood \[\begin{align*} \ell(\theta|y) &= \log[\mathcal{L}(\theta|y)] \\ &= \log \left[ \prod_{i=1}^np(y_i|\theta) \right]\\ &= \sum_{i=1}^n \log[p(y_i|\theta)] \,. \end{align*}\]

Dies kleine Umformung bringt drei Vorteile mit sich:

• Der Logarithmus verwandelt Produkte in Summen, was rechnerisch viel einfacher zu handhaben ist. Dies geht auf eine Grundregel der Rechnung mit Logarithmen zurück: \(\log(x\cdot y) = \log(x) + \log(y)\), solange \(x > 0\) und \(y > 0\).
• Der Logarithmus kleiner Zahlen ergibt Werte mit größerem absoluten Betrag, z.B. \(\log(0.0001) = –9.21\). Umgekehrt führt der Logarithmus großer Zahlen zu Werten mit kleinerem absoluten Betrag, z.B. \(\log(10000) = 9.21\).
• Der Logarithmus ist eine streng monoton steigende Funktion, sodass das Maximum der log-Likelihood an genau derselben Stelle liegt wie das Maximum der ursprünglichen Likelihood-Funktion. Das bedeutet, dass man durch Maximierung der log-Likelihood dieselben Parameterschätzungen erhält. Dies illustriert Abbildung 2.14, die nun die log-Likelihood für das obige Beispiel als Funktion von \(\mu\) zeigt.

Abbildung 2.14: log-Likelihood in Abhängigkeit von \(\mu\).

Der letzte Schritt ist noch: Die Optimierungsalgorithmen wie Simplex oder DE sind üblicherweise darauf ausgelegt, ein Minimum zu finden. Wir sind hier aber am Maximum der (log)-Likelihood-Funktion interessiert. Der Trick hier ist ganz einfach: Die log-Likelihood wird mit -1 multipliziert und aus Maxima werden dadurch Minima. Daher ist der Ansatz zur Bestimmung der Werte in \(\theta\), das Minimum der negativen log-Likelihood \[ -\ell(\theta|y) \] zu bestimmen. In der Regel übernehmen dies wieder die oben besprochenen Algorithmen wie Simplex oder DE. Die MLE stellt daher im Prinzip eine weitere Art von Kostenfunktion zur Verfügung, aus der allerdings, vermöge ihrer Verankerung in der statistischen Theorie, eine Reihe von Vorteilen erwachsen.

2.6.4 Likelihood-Ratio Tests, AIC und BIC

Ein Likelihood-Ratio Test (LRT) basiert auf der (skalierten) Differenz der Likelihoods (bzw. der Deviance-Statistik) zweier genesteter Modelle (siehe auch bereits hier im Kontext der logistischen Regression). Unter der Deviance-Statistik versteht man \[ \text{Deviance}=-2\cdot \ell\,, \] also die zweifache negative log-Likelihood.

Konkret wird hierbei ein Obermodell mit allen Parametern mit einem Untermodell, bei dem die interessierenden Parameter weggelassen werden, auf Basis der (maximierten) Likelihoods verglichen. Der Name “Ratio” rührt dabei daher, dass ein Bruch der beiden Likelihoods betrachtet wird, der aber durch eine Logarithmierung oft in eine Differenz überführt wird. Die resultierende Teststatistik ist dabei approximativ \(\chi^2\)-verteilt mit so vielen Freiheitsgraden wie Parameter im Vergleich vom Ober- zum Untermodell weggelassen wurden: \[\begin{align*} -2\log\left( \frac{\mathcal{L}_{\text{Untermodell}}}{\mathcal{L}_{\text{Obermodell}}} \right) &= [-2\log(\mathcal{L}_{\text{Untermodell}})]-[-2\log(\mathcal{L}_{\text{Obermodell}})]\\ &=2\cdot[\log(\mathcal{L}_{\text{Obermodell}})-\log(\mathcal{L}_{\text{Untermodell}})]\underset{app.}{\sim}\chi^2_{df = \Delta p}\,. \end{align*}\] Ausgedrückt durch die Deviance ergibt sich \[ \text{Deviance}_\text{Untermodell} - \text{Deviance}_\text{Obermodell} \underset{app.}{\sim}\chi^2_{df = \Delta p}\,. \]

LRTs können in R (v.a. für Regressionsmodelle) durch eine Übergabe von Ober- und Untermodell der Funktion anova() erhalten werden.

Um (nicht-genestete) Modelle mit verschiedenen Parametern zu vergleichen, oder auch, um zu berücksichtigen, dass Modelle mit mehr Parametern oft welchen mit weniger Parametern überlegen sind, wurde vorgeschlagen, die log-Likelihood noch für die Anzahl der Parameter zu korrigieren. Ein sehr einflussreicher Vorschlag dazu stammt von Akaike (1974) und wird heute Akaike Information Criterion (AIC) genannt: \[AIC=-2\ell+2p\] wobei \(p\) die Anzahl der Parameter darstellt. Für diese Formel gibt es noch verschiedene Korrekturen für die Stichprobengröße, wobei die vielleicht wichtigste Alternative das Bayesian Information Criterion (BIC) ist (Schwarz, 1978): \[BIC = -2\ell+p\log(n)\] wobei (streng genommen) \(n\) der sogenannte “limitierende Stichprobenumfang” ist, was aber bei kontinuierlich verteilten Daten der Anzahl der Datenpunkte \(n\) entspricht.

Da die log-Likelihood mit einem negativen Vorzeichen in die Berechnung beider Maße eingeht, sprechen nun wieder kleinere Werte für einen besseren Fit.

2.7 Zusammenfassung

In diesem Kapitel haben wir uns mit den technischen Voraussetzung für die Schätzung von Parametern in mathematischen Modellen befasst. Der Einfachheit halber haben wir auf das bekannte Problem der einfachen, linearen Regression zurückgegriffen. Wichtig ist festzuhalten, dass Kostenfunktion (\(Q\), \(RMSE\), neg. log-Likelihood, …) und Optimierungsalgorithmus im Prinzip beliebig miteinander kombinierbar sind. Wir werden aus diaktischen Gründen im nächsten Kapitel 3 auch mit dem \(RMSE\) arbeiten, für die Diffusionsmodelle in Kapitel 4 und 5 verwenden wir dann ausschließlich MLE zur Bestimmung der Parameter. Hierbei greifen wir dann auch auf unser R-Paket dRiftDM zurück, worin derlei Dinge direkt implementiert sind.

3 Random Walks als Modelle der Reaktionsauswahl

In diesem Kapitel führen wir mit Random Walks ein erstes stochastisches Modell ein, welches Reaktionswahlen und RTs erklären kann. Wir beginnen dazu mit einem verbalen Modell der Reaktionsauswahl und überführen intutive Annahmen dann in ein mathematisches Modell. Im Anschluss behandeln wir den Random Walk erst klassisch per Simulation und kommen dann zu einer Darstellung, mittels derer relevante Größen der Reaktionswahlen und RTs sehr effizient und generalisiert bestimmt werden können. Der einfachste Random Walk ist dann auch die Ausgangsbasis für das Diffusionsmodell, welches im folgenden Kapitel dann eingeführt wird.

3.1 Reaktionsauswahl in psychologischen Stadienmodellen

Wie können Random Walks uns helfen, die Entstehung einer Reaktion und der damit zusammenhängenden RT zu verstehen? Aus der Allgemeinen bzw. Kognitiven Psychologie stammen Vorstellungen der seriellen Abläufe mehrerer Stufen zur Bearbeitung typischer Aufgabe (siehe Abb. 3.1).

Die Stufe der Reaktionsauswahl soll also eine Repräsentation des wahrgenommenen Stimulus “übersetzen” in eine Repräsentation einer motorischen Efferenz, in Folge derer dann eine Taste gedrückt wird und das Ende des RT-Intervalls gemessen wird. Wie kann diese Übersetzung verstanden und idealerweise formalisiert werden? Eine Idee hierbei ist, dass der Reaktionsauswahlprozess mit jedem Zeitschritt Evidenz für die eine oder die andere Reaktion ansammelt (was sich z.B. im Aktivationsniveau von Neuronen(verbänden) niederschlägt). Hat das kognitive System genug Evidenz für die eine (oder eine andere) Reaktion gesammelt, gilt die Reaktion als gewählt und die motorische Efferenz wird programmiert. Hierbei müssen nun zwei weitere Aspekte beachtet werden:

1. Stimuli können sich dahingehend unterscheiden, wie “gut” die aus ihnen extrahierte Evidenz für eine Reaktion ist: Ein sehr gut sichtbar präsentierter Stimulus liefert mehr Evidenz als ein schlecht sichtbarer Stimulus.
2. Neurophysiologisch ist nicht davon auszugehen, dass die Evidenzsammlung deterministisch abläuft, d.h., dass mit jedem Zeitschritt einfach der gleiche Evidenzwert addiert wird. Vielmehr ist von einem zusätzlichen Rauschen auszugehen, welches zufällige Schwankungen der Evidenzsammlung bewirkt. Andernfalls wäre die Zeit bis zum Auswählen einer Reaktion auch immer identisch. Dass RTs aber klar einer Variabilität unterliegen, haben wir aber weiter oben bereits festgestellt.

Nun sollte nicht ein vollkommen beliebiges Ensemble aus Systematik und Zufall als Modell dieser Beschreibung verwendet werden, sondern zu bevorzugen wäre eines, welches systematisch untersucht werden kann und damit auch in seinem Verhalten verstanden werden kann. Random Walks sind ein Beispiel dafür und eine ihrer Anwendungen liegt nun genau darin, sie als Modell für die Übersetzung eines Stimulus in eine Reaktion anzusehen.

Abbildung 3.1: Illustration eines typischen Stufenmodells aus der kognitiven Psychologie.

3.2 Stochastische Prozesse und Markov-Ketten

Stochastische Prozesse bilden die mathematische Grundlage von Evidenzakkumulationsmodellen, zu denen Random Walks, aber eben auch Diffusionsmodelle, gehören. In unserem Kontext ist ein stochastischer Prozess eine Sammlung von Zufallsvariablen \[ \{\mathbf{X}_t\}_{t\in T}\,, \] die Werte aus einer Menge \(S\) annehmen, dem sog. state space, der entweder diskret oder kontinuierlich sein kann. Die Indexmenge \(T\) wird in unserem Zusammenhang als Zeit bzw. time space bezeichnet und kann ebenfalls entweder diskret oder kontinuierlich sein.

Grundsätzlich entwickelt sich ein stochastischer Prozess also über die Zeit und nimmt verschiedene Werte aus \(S\) an. Dies geschieht mit Unsicherheit, das heißt auf stochastische Weise abhängig vom Zufall. Zu jedem Zeitpunkt (bzw. Index) \(t\) hängt der Zustand nur vom Zustand zum Zeitpunkt \(t-1\) ab, nicht jedoch davon, wie der Zustand bei \(t-1\) erreicht wurde. Wir betrachten hier nun den einfachsten Fall, bei dem sowohl \(T\) als auch \(S\) diskret sind, also \(t\in \mathbb{N}_0\) und \(X_t\in \mathbb{Z}\). Ein solcher Prozess ist auch als Markov-Kette bekannt.

Betrachten wir zunächst ein sehr einfaches Beispiel zur Wettervorhersage und führen dabei einige Konzepte ein, denen wir später wieder begegnen werden. Der Einfachheit halber ist das Wetter entweder sonnig (\(s\)) oder regnerisch (\(r\)). Die Zufallsvariablen \(\mathbf{X}_0\), \(\mathbf{X}_1\), \(\mathbf{X}_2\), … nehmen jeweils den Wert \(s\) oder \(r\) an und repräsentieren das Wetter im Zeitverlauf. Wenn wir außerdem annehmen, dass das Wetter von morgen nur vom heutigen Wetter abhängt, fügen wir die Markov-Eigenschaft hinzu. In diesem Fall können wir die Wetterentwicklung durch eine Transitionsmatrix

\[ P= \begin{bmatrix} 0.7 & 0.3 \\ 0.2 & 0.8 \end{bmatrix}\,. \] darstellen. In dieser Matrix stellen die Zeilen den aktuellen Zustand am Tag \(t\) dar, also den Zustand, aus dem der Prozess kommt, und die Spalten den Zustand am Tag \(t+1\), also den Zustand, in den der Prozess übergeht. Die erste Zeile und Spalte stehen für einen sonnigen Tag, die zweite Zeile und Spalte für einen regnerischen Tag.

Das bedeutet: Wenn es heute sonnig ist, wird es mit Wahrscheinlichkeit \(p = 0.7\) auch morgen sonnig sein (der Prozess bleibt im gleichen Zustand) und mit Wahrscheinlichkeit \(p = 0.3\) regnerisch (der Prozess wechselt den Zustand). Ebenso gilt: Wenn es heute regnet, wird es mit Wahrscheinlichkeit \(p = 0.8\) auch morgen regnen und mit Wahrscheinlichkeit \(p = 0.2\) sonnig.

Formaler gesagt, beschreibt jeder Eintrag \(p_{ij}\in P\) eine Übergangswahrscheinlichkeit, nämlich \[ p_{ij} = P(X_{t+1} = j \mid X_t = i)\,. \] Die Matrix \(P\) wird als stochastische Matrix bezeichnet, wenn alle Einträge nicht-negativ sind und die Einträge jeder Zeile sich jeweils zu 1 summieren. Der Random Walk ist ein Beipiel eines Markov-Kette und wir werden ihn gegen Ende dieses Kapitels auch in Form einer Transitionsmatrix \(P\) schreiben.

3.3 Grundlagen von Random Walk Modellen

Mit diesem Vorwissen entwickeln wir nun ein Prozessmodell, welches (a) plausible Annahmen über die Entstehung einer Reaktion und damit der dazugehörigen RT macht und (b) in der Lage ist, die spezifische Verteilungsform der RTs vorherzusagen. Dabei beginnen wir mit sog. Random Walk Modellen, die uns dann im nächsten Kapitel zum Diffusionsmodell bringen werden.

3.3.1 Einfacher, symmetrischer Random Walk

Random Walks sind einfache Beispiele eines stochastischen Prozesses. Eine Zufallsvariable \({\bf X}_t\) beschreibt dabei die Position des “Walks” im state space \(S\) zum Zeitpunkt \(t=0, 1,2,\ldots\). Für die späteren Schritte \(t>0\) ergibt sich dann \[{\bf X}_t = {\bf X}_{t-1}+{\bf Y}_t\,,\] wobei die \({\bf Y}_t\) unabhängig voneinander und Bernoulli-verteilt seien. Diese Darstellung nennt man auch eine stochastische Differenzgleichung; “einfach” meint in diesem Kontext, dass die Ausgangsbedingung \({\bf X}_0=0\) ist (d.h. initial ist die Position 0). “Symmetrisch” meint, dass in jedem Schritt die Wahrscheinlichkeit eines Sprungs nach oben gleich der Wahrscheinlichkeit eines Sprungs nach unten ist, also \[\begin{align} P({\bf Y}_t = 1) = P({\bf Y}_t = -1) = 0.5\,. \end{align}\] In Worten bedeutet dies: Der Random Walk startet an der Position 0. Zum Zeitpunkt \(t=1\) springt die Position mit gleicher Wahrscheinlichkeit nach oben oder nach unten. Mögliche Positionen sind dann also 1 und -1. Von dort aus geht es wiederum zum Zeitpunkt \(t=2\) mit gleicher Wahrscheinlichkeit einen Schritt nach oben oder nach unten. Die nun möglichen Positionen sind also 2, 0 und -2. Diesen Gedanken können wir nun weiter fortführen. Ein solcher Random Walk kann aber auch sehr leicht in R simuliert werden, wobei wir uns hier auf die ersten \(t=20\) Zeitpunkte beschränken, aber bereits einige typische Aspekte des Programmierens hier einführen. Der folgende Code produziert den linken Teil der Abbildung 3.2, die einen einzelnen Random Walk visualisiert. Im rechten Teil von Abbildung 3.2 sind 5 Realisierungen gleichzeitig dargestellt:

# wir setzen einen seed, damit das Ergebnis der folgenden Zufallsziehung
# replizierbar ist... wird der seed verändert, erhalten wir auch ein anderes
# Ergebnis
set.seed(2)

# Definition einiger Parameter der folgenden Simulation
t_max <- 20 # maximale Anzahl an Zeitschritten
p <- 0.5 # Wahrscheinlichkeit eines Schritts nach oben
X_0 <- 0 # Ausgangsbedingung zu t = 0

# Anlegen eines Vektors mit den diskreten Zeitschritten 0, ..., t_max
t <- c(0:t_max)

# Hier werden nun 20-mal die Werte 1 oder -1 mit gleicher Wahrscheinlichkeit
# und mit Zurücklegen gezogen. Anschließend wird die kumulative Summe dieses
# Vektors gebildet, was entsprechend die Positionen dieser Realisierung des
# Random Walks sind
X_t <- cumsum(sample(
  x = c(-1, 1),
  size = t_max,
  prob = c(p, 1 - p),
  replace = TRUE
))
# Der Vollständigkeit halber müssen nun an die Position t = 0 noch den
# Wert der Ausgangsbedingung einfügen
X_t <- c(X_0, X_t)

# Dann können wir uns den Random Walk plotten
# (linker Plot der nachfolgenden Abbildung)
plot(X_t ~ t,
  type = "b",
  ylim = c(-5, 5),
  xlab = expression(plain("Zeitschritt ") ~ italic(t)),
  ylab = expression(bold(X)[t]),
  main = "Einfacher, symmetrischer Random Walk"
)

Abbildung 3.2: Beispiele für einfache, symmetrische Random Walks.

Wir können hier bereits zwei Aspekte eines Random Walks festhalten:

• Der Random Walk ist ein stochastischer Prozess in diskreter Zeit, da die Zeitschritte nur Werte aus \(\mathbb{Z}_{\geq 0}\) sind.
• Der Random Walk kann auch nur diskrete Positionen im state space \(S\) annehmen, in unserem Fall prinzipiell alle ganzen Zahlen, also \(\mathbb{Z}\). Die Position wird in der mathematischen Theorie auch als State bezeichnet.

Jeder einzelne Random Walk entspricht hier dem Pfad einer Kugel eines Galton-Bretts, welches auch zur Illustration der Entwicklung einer Normalverteilung genutzt werden kann (da es allerdings immer endlich viele Kugeln in einer Demonstration sind, wird dann eigentlich eine Binomial-Verteilung demonstriert).

Die Zufallsvariable \({\bf X}_t\) ist also die Summe aus \(t\)-vielen unabhängigen Zufallsvariablen \({\bf Y}_t\), und mit \({\bf X_0} = 0\) also \[{\bf X}_t = {\bf X}_0 + \sum_{i = 1}^t{\bf Y}_i = \sum_{i = 1}^t{\bf Y}_i\,.\] Angenommen nun, wir würden unendlich viele Random Walks erzeugen und dann zu jedem Zeitpunkt \(t\) die mittlere Position sowie die Varianz der möglichen Positionen berechnen wollen. Anders gesagt, interessieren wir uns für den Erwartungswert und die Varianz von \({\bf X}_t\). Solche Fragen sind typisch für die Analyse stochastischer Prozesse und im vorliegenden Fall ist die Analyse noch sehr einfach (siehe auch hier):

• Der Erwartungswert von \({\bf Y}_t\) ist wegen \(P({\bf Y}_t = 1) = P({\bf Y}_t = -1) = 0.5\) leicht zu berechnen als \(\mathbb{E}({\bf Y}_t) = 0.5\cdot 1 + 0.5\cdot -1 = 0\). Da sich \({\bf X}_t\) als Summe dieser Zufallsvariablen ergibt, folgt daraus auch \[\mathbb{E}({\bf X}_t) = 0\,.\]
• Zunächst berechnet sich die Varianz von \({\bf Y}_t\) als \[\begin{align*} V({\bf Y}_t) &= \mathbb{E}[({\bf Y}_t-\mathbb{E}({\bf Y}_t))^2] \\ &= (-1 - 0)^2\cdot 0.5+(1-0)^2\cdot 0.5 \\ &= 1\,. \end{align*}\] Die Varianz von \({\bf X}_t\) ist dann die Summe der Varianzen von \({\bf Y}_t\) und wegen derer Unabhängigkeit gilt \[\begin{align*} V({\bf X}_t) &= \sum_{i = 1}^t V({\bf Y}_t) \\ &= \sum_{i = 1}^t 1\\ &= t\,. \end{align*}\] Diese Ergebnisse werden im Kontext Diffusionsmodelle noch einmal wichtig werden.

3.3.2 Einfacher Random Walk mit Drift

Als nächstes erweitern wir den Random Walk, indem wir Wahrscheinlichkeiten von \(p\neq 0.5\) für die \({\bf Y}_t\) zulassen, also \[\begin{align} P({\bf Y}_t = 1) = p\quad\text{und}\quad P({\bf Y}_t = -1) = 1-p\,. \end{align}\] In Worten bedeutet dies also, dass eine Versetzung in der Position in jedem Zeitschritt mit einer Wahrscheinlichkeit von \(p\) nach oben erfolgt (1) und mit der Gegenwahrscheinlichkeit \(1-p\) nach unten erfolgt (-1). In Bezug auf die oben angesprochene Formalisierung von Reaktionsauswahl könnte man sich bereits vorstellen, dass bei einem Stimulus \(S_1\) tendenziell ein Sprung nach oben (also ein positiver State) induziert wird, während der andere Stimulus \(S_2\) eher einen Sprung nach unten, also einen negativen State induziert. Wie verhalten sich dann Erwartungswert und Varianz von \({\bf X}_t\)?

• Der Erwartungswert einer Versetzung der Position von \(t-1\) zu \(t\) ergibt sich nun als \[\begin{align*} \mathbb{E}({\bf Y}_t) &= [1 p] + [(-1)\cdot (1-p)] \\ &= p -1 + p\\ &= 2 p -1\,. \end{align*}\] Da sich \({\bf X}_t\) als Summe \(t\)-vieler unabhängiger Zufallsvariablen \({\bf Y}_t\) ergibt, folgt \[\begin{align*} \mathbb{E}({\bf X}_t) &= \mathbb{E}\left( \sum_{i=1}^t {\bf Y}_t \right) \\ &= \sum_{i=1}^t \mathbb{E}({\bf Y}_t) \\ &= t (2 p-1)\,. \end{align*}\]
• Die Berechnung der Varianz von \({\bf Y}_t\) ist etwas umfangreicher, aber nicht kompliziert, und es ergibt sich \[\begin{align*} V({\bf Y}_t) &= \mathbb{E}[({\bf Y}_t-\mathbb{E}({\bf Y}_t))^2] \\ &= (1-[2p-1])^2\cdot p + (-1-[2p-1])^2\cdot (1-p) \\ &= (1-2p+1)^2\cdot p + (-1-2p+1)^2\cdot (1-p) \\ &= (2-2p)^2\cdot p + (2p)^2\cdot(1-p) \\ &= (4-8p+4p^2)\cdot p + 4p^2\cdot (1-p) \\ &= 4p-8p^2+4p^3+4p^2-4p^3 \\ &= 4p - 4p^2 \\ &= 4p(1-p)\,. \end{align*}\] Die Varianz von \({\bf X}_t\) ist dann wieder die Summe der Varianzen von \({\bf Y}_t\) und wegen derer Unabhängigkeit gilt \[\begin{align*} V({\bf X}_t) &= \sum_{i = 1}^t V({\bf Y}_t) \\ &= \sum_{i = 1}^t 4p(1-p)\\ &= 4pt(1-p)\,. \end{align*}\]

Wir visualisieren die Situation für \(p=0.7\) nun noch einmal in der gleichen Art wie bereits oben für den symmetrischen Fall (siehe Abb. 3.3). Dabei sehen wir nun, dass die einzelnen Random Walks eine gewisse Tendenz haben mit zunehmenden Zeitschritten positive States anzunehmen (also würde ggf. Stimulus \(S_1\) vorliegen). Wir können aber auch nicht ausschließen, dass zufallsbedingt vereinzelte Random Walks auch negative States annehmen.

Abbildung 3.3: Beispiele für Random Walks mit Drift.

3.3.3 Zusammenfassung einfacher Random Walk Modelle

Wir haben nun Random Walks als einfache stochastische Prozesse kennengelernt. Random Walks sind diskret sowohl in der Zeit als auch im State. Während wir hier noch relativ einfache Analysen des Prozesses (z.B. Erwartungswert und Varianz zu jedem Zeitschritt \(t\) durchführen konnten), werden solche Analysen bei komplexeren stochastischen Prozessen schnell komplizierter. Wir werden in den folgenden Abschnitten Random Walks zunächst kurz simulationsbasiert explorieren, bevor wir uns danach anderen Lösungen zuwenden; mehr formale Details zu Random Walks (und vielen anderen stochastischen Prozessen) finden sich in Cox und Miller (1965), Schwarz (2022), Diederich und Busemeyer (2003) sowie Diederich and Mallahi-Karai (2018).

3.4 First-Passage Time Verteilung von Random Walks

Die grundlegende Idee bei binären Entscheidungen (bspw. die Entscheidung über “linke” Taste vs. “rechte” Taste) ist nun: Die Sammlung von Evidenz für eine Reaktion beginnt bei 0 als initiale Bedingung. “Reagiert” das Modell auf einen Stimulus \(S_1\), der die Reaktion \(R_1\) verlangt, dann ist die Wahrscheinlichkeit für eine Stateveränderung nach oben \(p>0.5\). Verlangt der Stimulus \(S_2\) hingegen nach der anderen Reaktion \(R_2\), dann ist \(p<0.5\). Damit erfolgt die Sammlung von Evidenz nicht deterministisch, denn in jedem Zeitschritt kann die Evidenz auch in die “falsche” Richtung gehen. Allerdings weist die Sammlung von Evidenz je nach Stimulus einen “Trend” auf.

Bisher haben wir keinerlei Begrenzung in der zeitlichen Entwicklung der stochastischen Prozesse betrachtet. Da aber Versuchspersonen ja durchaus reagieren, muss dieser Prozess, und damit die Entscheidung für \(R_1\) oder \(R_2\), zu irgendeinem Zeitpunkt beendet werden. Die Idee ist, dass es einen bestimmten “Evidenz”-Wert gibt, der als “ausreichend Evidenz für eine Reaktion” zählt. Dieser Wert wird in der Literatur oft Boundary oder Threshold genannt und wir nennen ihn als Modellparameter \(b\) (obere Reaktion, \(R_1\)) bzw. \(-b\) (untere Reaktion, \(R_2\)). Der Zeitschritt \(t\), zu dem zum ersten Mal gilt \[{\bf X}_t \geq b\quad\text{oder}\quad{\bf X}_t \leq -b\,,\] wird als First-Passage Time bezeichnet und umfasst dann denjenigen Anteil einer RT, der auf die Entscheidung, d.h. auf die Reaktionsauswahl entfällt (man sagt dann auch oft, dies sei die Decision Time). Durch die stochastische Bestimmung der Stateveränderung kann es im Extremfall auch sein, dass bei \(p>0.5\) der stochastische Prozess dennoch zuerst \(-b\) unterschreitet und umgekehrt. In diesem Fall würde entsprechend die falsche Reaktion ausgewählt werden, was durchaus sinnvoll vor dem Hintergrund empirischer Daten ist. In der mathematischen Literatur wird diese Situation mit zwei Boundaries auch two absorbing barriers genannt, da der stochastische Prozesse von der Boundary “absorbiert” wird und nicht mehr weiterläuft.

Die Bestimmung der First-Passage Time bedeutet mathematisch, dass analyisiert werden muss, wann \({\bf X}_t\) zum ersten Mal den Zustand \(b\) oder \(-b\) erreicht. Die analytische Bestimmung ist aber leider nicht mehr so einfach wie im Fall der Erwartungswerte und Varianzen weiter oben. Wir werden daher zunächst einen anderen Ansatz wählen und die (erwartete) First-Passage Time als Mittelwert der First-Passage Times einer großen Zahl simulierter Random Walks bestimmen. Dies ermöglicht zugleich, den Anteil beider Reaktionen \(R_1\) und \(R_2\) bzw. anders ausgedrückt: den Anteil an Fehlern zu bestimmen und eine Idee der Verteilung der einzelnen First-Passage Times zu bekommen.² Mit der Simulationsmethode greifen wir auch nochmal die Parameterschätzung am Beispiel des Random Walks auf. Im Anschluss betrachten wir noch kurz eine generelle Methode zur Bestimmung von First-Passage Times und Reaktionswahlen basierend auf Transitionsmatrizen.

3.4.1 Implementation eines Random Walks

Wir benötigen nun zunächst eine Funktion, die eine große Zahl Random Walks simulieren kann und ein DataFrame zurückliefert, in welchem für jeden Random Walk die First-Passage Time zu finden ist und zusätzlich, ob der stochastische Prozess zuerst \(b\) oder \(-b\) über- bzw. unterschritten hat. Zur Verdeutlichung des Prinzips und der Erhöhung der Lesbarkeit arbeiten wir hier mit einer for-Schleife. Dies geht zwar zu Lasten der Geschwindigkeit, fällt aber hier nicht wirklich auf:

# Funktion zur Simulation vieler Random Walks
simulate_RW <- function(t_max = 200, # Wieviele Zeitschritte?
                        n_trials = 10, # Wieviele Trials / Random Walks?
                        p_up = 0.5, # Wahrscheinlichkeit nach oben zu gehen
                        boundary = 10) { # Boundary

  # zwei Vektoren vorbereiten, in die die First-Passage Time und die
  # Reaktion geschrieben werden
  firstPassageTimes <- numeric(n_trials)
  response <- numeric(n_trials)

  for (i in c(1:n_trials)) { # n_trial viele Random Walks
    # Generierung eines Random Walks
    X <- cumsum(c(0, sample(
      x = c(1, -1),
      size = t_max,
      replace = TRUE,
      prob = c(p_up, 1 - p_up)
    )))
    # Ermittlung des Zeitschritts wenn First-Passage erfolgt und welche Boundary
    # dabei über-/unterschritten wurde. Dabei heißt 1 = b überschritten und
    # -1 = -b unterschritten
    firstPassage <- which(abs(X) > boundary)[1] # Zeitschritt
    response[i] <- sign(X[firstPassage]) # welche Boundary?
    firstPassageTimes[i] <- firstPassage - 1
  }

  # dann werden Reaktionen und First-Passage Times zu einem DataFrame zusammen-
  # gefasst, wobei noch alle NAs entfernt werden, d.h. die Instanzen, wo der
  # Random Walk bis zum Zeitpunkt t_max weder b noch -b über-/unterschritten hat
  RW <- data.frame(
    response,
    firstPassageTimes
  )
  RW <- na.omit(RW)

  # Schließlich wird der DataFrame an eine aufrufende Funktion zurückgegeben
  return(RW)
}

Diese Funktion kann nun wie jede andere Funktion in R aufgerufen und verwendet werden. Simulieren wir bspw. 6 Random Walks mit \(p=0.7\) und \(b=40\) (und maximal 200 Zeitschritten), dann würde der Aufruf wie folgt aussehen:

myRW <- simulate_RW(
  t_max = 200, # Wieviele Zeitschritte?
  n_trials = 6, # Wieviele Trials / Random Walks?
  p_up = 0.7, # Wahrscheinlichkeit nach oben zu gehen
  boundary = 40
)
myRW # DataFrame anzeigen

3.4.2 Analyse der (simulierten) First-Passage Times

Zur Analyse und Exploration des Modells benötigen wir nun noch Zusammenfassungen der Ergebnisse einer Simulation. Dazu schreiben wir uns eine kleine Funktion analyze_sim(), die die mittleren First-Passage Times für beide Reaktionen sowie deren relative Häufigkeiten berechnet und ausgibt. In der Funktion müssen wir ein wenig aufpassen, dass wir sicherstellen, ob auch wirklich Reaktionen oben und unten überhaupt vorliegen. Daher werden in der folgenden Funktion zwei Fälle unterschieden:

analyze_sim <- function(data) {
  # Gibt es Reaktionen durch die obere Boundary (1)?
  if (any(data$response == 1)) {
    # Relative Häufigkeit und Mittelwert berechnen
    mean_resp_1 <- mean(data$firstPassageTime[data$response == 1])
    perc_resp_1 <- mean(data$response == 1)
    # Formatieren und anzeigen
    upper <- paste0(
      "Rel. Häufigkeit obere Boundary: ", round(perc_resp_1, 3),
      "; M = ", round(mean_resp_1)
    )
    print(upper)
  }

  # Gibt es Reaktionen durch die untere Boundary (-1)?
  if (any(data$response == -1)) {
    # Relative Häufigkeit und Mittelwert berechnen
    mean_resp_2 <- mean(data$firstPassageTime[data$response == -1])
    perc_resp_2 <- mean(data$response == -1)
    # Formatieren und anzeigen
    lower <- paste0(
      "Rel. Häufigkeit untere Boundary: ", round(perc_resp_2, 3),
      "; M = ", round(mean_resp_2)
    )
    print(lower)
  }
}

Nun simulieren wir eine etwas realistischere Anzahl Random Walks, indem wir n_trials = 10000 setzen. Beim Aufruf sehen wir auch, dass dann auch mehr Rechenzeit benötigt wird. Da wir hier einen symmetrischen Random Walk mit p_up = 0.5 (also ohne Drift) simulieren, ist davon auszugehen, dass die einzelnen Random Walks um den State 0 mäandrieren, ohne schnell die obere oder untere Boundary zu über- bzw. unterschreiten. Wir setzen daher auch die maximale Anzahl der Zeitschritte mit t_max = 2000 auf einen eher großen Wert, um den Random Walks zumindest in vielen Fällen eine Chance geben, eine Boundary durch Zufall zu über- bzw. unterschreiten:

set.seed(1)
myRW <- simulate_RW(
  t_max = 2000, # Wieviele Zeitschritte?
  n_trials = 10000, # Wieviele Trials / Random Walks?
  p_up = 0.5, # Wahrscheinlichkeit nach oben zu gehen
  boundary = 10
)
analyze_sim(myRW)

## [1] "Rel. Häufigkeit obere Boundary: 0.505; M = 121"
## [1] "Rel. Häufigkeit untere Boundary: 0.494; M = 121"

Erwartungsgemäß treten beide Reaktionen beinahe gleich häufig auf und auch die mittleren First-Passage Times sind relativ ähnlich. Nun schauen wir uns noch die Verteilung der First-Passage Times an, indem wir Histogramme mit relativen Häufigkeiten erstellen (siehe Abb. 3.4). Hier können wir sehen, dass die oben beschriebene linkssteile, rechtsschiefe Form tasächlich erzeugt wird.

par(mfrow = c(1, 2))
# linkes Histogramm
hist(myRW$firstPassageTimes[myRW$response == 1],
  breaks = 40,
  freq = FALSE,
  xlab = "First-Passage Time",
  ylab = "rel. Häufigkeit",
  main = "Reaktion oben",
  ylim = c(0, 0.01)
)
# rechtes Histogramm
hist(myRW$firstPassageTimes[myRW$response == -1],
  breaks = 40,
  freq = FALSE,
  xlab = "First-Passage Time",
  ylab = "rel. Häufigkeit",
  main = "Reaktion unten",
  ylim = c(0, 0.01)
)

Abbildung 3.4: Histogramme der First-Passage Times oberer und unterer Reaktionen bei einem einfachen, symmetrischen Random Walk.

Aufgabe zum Ausprobieren: Hier ist es nun instruktiv, das Verhalten des Modells zu explorieren und festzustellen, was sich für verschiedene Parameterkonstellationen ergibt. Achten Sie dabei besonders auf die mittleren First-Passage Times, wenn \(p\neq 0.5\) gewählt wird. Welches unintuitive Verhalten des Modells fällt dabei auf? Variieren Sie darüber hinaus auch den Wert für die Boundary und versuchen Sie dabei die Systematik der Veränderungen zu erkennen!

3.5 Parameterschätzung für ein einfaches Random Walk Modell

Zur Illustration wollen wir nun die Inhalte aus Kapitel 2 zur Parameterschätzung auf ein Random Walk Modell anwenden. Die beiden Parameter, die es dabei zu bestimmen gilt, sind die Boundary \(b\) sowie \(p\). Mangels “echter” empirischer Daten (und da wir bisher auch nur über den Teil “Reaktionsauswahl” geredet haben und Wahrnehmung sowie Motorik bisher ignoriert haben), erzeugen wir uns quasi die beobachteten (empirischen) Daten synthetisch selber mit dem Modell und versuchen im Anschluss die erzeugenden Parameter “zurückzuschätzen”:

# Erzeugung der "beobachteten" empirischen Daten
set.seed(1)
obsData <- simulate_RW(
  t_max = 200, # Wieviele Zeitschritte?
  n_trials = 1000, # Wieviele Trials / Random Walks?
  p_up = 0.6, # Wahrscheinlichkeit nach oben zu gehen
  boundary = 12 # Boundary, die erreicht werden muss
)
analyze_sim(obsData)

## [1] "Rel. Häufigkeit obere Boundary: 0.996; M = 64"
## [1] "Rel. Häufigkeit untere Boundary: 0.004; M = 69"

Da hier das verwendete p_up = 0.6 für eine obere Reaktion spricht, werten wir obere Reaktionen entsprechend als “korrekte Reaktionen” und untere Reaktionen als “Fehler”. Als nächstes müssen wir die Daten irgendwie adäquat beschreiben, damit die beobachteten und die dann vom Modell vorhergesagten (d.h. die jeweils simulierten) Daten in einer Kostenfunktion miteinander verglichen werden können. Man könnte zwar prinzipiell die gerade errechneten Mittelwerte und relativen Häufigkeiten der richtigen und fehlerhaften Reaktionen verwenden, eine “reichhaltigere” Beschreibung, die auch die Form der Verteilung berücksichtigt, erweist sich aber als sinnvoller. Dies liegt u.a. daran, dass es viele Parameterkonstellationen geben kann, die die beobachteten Mittelwerte erzeugen würden.

Betrachten wir dazu zunächst die korrekten First-Passage Times. Zunächst halten wir fest, dass wir durch Simulationen zwar eine Art Verteilung der First-Passage Times erhalten, diese aber keine kontinuierliche Funktion (wie eine exGauß- oder eine Wald-Verteilung) liefern. Wir müssen also Datenpunkte auswählen, die wir zur Beschreibung nutzen wollen. Eine häufig genutzte Möglichkeit ist dabei, auf Quantile der Verteilungsfunktion zurückzugreifen, die wir mit der Funktion quantile() sehr einfach bestimmen können. Wir verwenden hier die 0.1, 0.3, 0.5, 0.7 und 0.9 Quantile. Abbildung 3.5 illustriert, dass diese Quantile die empirische Verteilungsfunktion (cumulative distribution function, CDF) gut zusammenfassen können:

# nur korrekte First-Passage Times verwenden
correct_FPT <- subset(
  obsData,
  response == 1
)
# Plotten der empirischen Verteilungsfunktion
plot(ecdf(correct_FPT$firstPassageTimes),
  cex = 0.1,
  main = "Verteilungsfunktion der korrekten First-Passage Times",
  xlab = "First-Passage Time",
  ylab = "F"
)

# Berechnen der (beobachteten) Quantile
obsQuantiles <- quantile(correct_FPT$firstPassageTimes,
  probs = c(0.1, 0.3, 0.5, 0.7, 0.9)
)
# Einzeichnen als schwarz-gefüllte Punkte
points(c(0.1, 0.3, 0.5, 0.7, 0.9) ~ obsQuantiles,
  pch = 19,
  col = "black"
)

Abbildung 3.5: CDF und Quantile der beobachteten (simulierten) Daten bei einem Random Walk.

Zusätzlich sollten wir auch die Fehlerraten (also den Anteil unterer Reaktionen) mit berücksichtigen. Eine einfache Möglichkeit hierfür wäre es, die \(-1\) für fehlerhafte Reaktionen durch eine \(0\) zu ersetzen und dann den Mittelwert der Variablen response zu berechnen, der dann dem Anteil korrekter Reaktionen entspricht:

obsData$response[obsData$response == -1] <- 0
mean(obsData$response) # rel. Häufigkeit korrekter Reaktionen

## [1] 0.9959718

Nun haben wir alle Bestandteile zusammen, die wir in Kapitel 2 mehrfach angesprochen hatten:

1. Das (Prozess-)Modell macht für bestimmte Parameterwerte Vorhersagen. Damit sind dann die vom Random Walk Modell vorhergesagten Verteilungen von First-Passage Times (ausgedrückt über die Quantile) und Fehlerraten (ausgedrückt über den Anteil korrekter Reaktionen) gemeint.
2. Die vorhergesagten Verteilungen und Fehlerraten werden mit den entsprechenden empirischen Werten (Quantilen und Anteil korrekter Reaktionen) verglichen und die Abweichung wird in einer Kostenfunktion ausgedrückt. Dazu werden wir hier nun den RMSE verwenden und insgesamt 6 Datenpunkte verrechnen: 5 Quantile und 1 Anteil korrekter Reaktionen.
3. Ein Optimierungsalgorithmus soll dann diejenigen Parameterwerte suchen, für die RMSE das Minimum annimmt.

Die Funktion, die die Modellvorhersagen generiert, haben wir mit simulate_RW() bereits geschrieben. Es fehlt also noch die Kostenfunktion, die diese Funktion dann aufruft und RMSE berechnet:

# Kostenfunktion basierend auf RMSE
RMSE_costFunction <- function(parameter, # Vektor mit den Parametern
                              obsData) { # empirische ("observed") Daten
  # Auslesen und Benennen der Parameter (zur besseren Lesbarkeit)
  my_p <- parameter[1]
  my_boundary <- parameter[2]

  # dann Modellvorhersagen für diese Parameter erzeugen
  predData <- simulate_RW(
    t_max = 200,
    n_trials = 1000,
    p_up = my_p,
    boundary = my_boundary
  )

  # Quantile der vorhergesagten und beobachteten korrekten First-Passage Times
  # berechnen und speichern
  obsQuantiles <- quantile(obsData$firstPassageTimes[obsData$response == 1],
    probs = c(0.1, 0.3, 0.5, 0.7, 0.9)
  )
  predQuantiles <- quantile(predData$firstPassageTimes[predData$response == 1],
    probs = c(0.1, 0.3, 0.5, 0.7, 0.9)
  )

  # Anteil korrekter Reaktionen in vorhergsagten und beobachteten Daten bestimmen
  obsData$response[obsData$response == -1] <- 0
  predData$response[predData$response == -1] <- 0
  obsPropCorrect <- mean(obsData$response)
  predPropCorrect <- mean(predData$response)

  # Nun können wir RMSE aus den 6 berücksichtigten Werten berechnen. Der besseren
  # Lesbarkeit halber machen wir das in mehreren Schritten
  squaredDeviationFirstPassageTimes <- sum((obsQuantiles - predQuantiles)^2)
  squaredDeviationAccuary <- (obsPropCorrect - predPropCorrect)^2
  RMSE <- sqrt(squaredDeviationFirstPassageTimes + squaredDeviationAccuary)

  # Wert der Kostenfunktion zurückgeben
  return(RMSE)
}

Nun versuchen wir die Parameter zu schätzen. Da wir hier Wahrscheinlichkeiten verwenden, ist es natürlich sinnvoll, mit eingeschränkten Wertebereichen und damit mit der Funktion nmkb() zu arbeiten, die einen Nelder-Mead Simplex Algorithmus implementiert, aber sog. box constraints erlaubt, also untere und obere Grenzen für die Parameter. Wir einigen uns darauf, dass der erste Wert im Parametervektor die Wahrscheinlichkeit für einen Schritt nach oben, also \(p\), darstellt, und der zweite Wert die Boundary \(b\) ist.

parameter_start <- c(0.70, 7) # Startwerte
parameter_lower <- c(0.5, 5) # untere Grenze
parameter_upper <- c(1, 20) # obere Grenze

result <- nmkb(
  par = parameter_start, # Startwerte
  lower = parameter_lower, # untere Grenze
  upper = parameter_upper, # obere Grenze
  fn = RMSE_costFunction, # Kostenfunktion
  obsData = obsData # empirische Daten die gefittet werden
)

Nun schauen wir uns die besten Parameterwerte und den kleinsten erzielten Wert für RMSE an. Wir können sehen, dass die Werte durchaus nahe an die zur Erzeugung der beobachteten Daten verwendeten Parameterwerte (p_up = 0.6 und boundary = 12) herankommen, aber diese (natürlich) auch nicht perfekt widerspiegeln:

result$par

## [1]  0.5983166 12.9529761

result$value

## [1] 0.006072813

Im aktuellen Beispiel haben wir den Vorteil, dass wir die “wahren” Parameter, die den beobachteten Daten zugrunde liegen, kennen. Daher wissen wir natürlich, wie gut oder schlecht die Paramerschätzung funktioniert hat, was im Normalfall aber nicht der Fall ist. Was könnten wir nun prinzipiell zur Verbesserung der Parameterschätzung tun?

• Im Beispielcode haben wir für jede Vorhersage nur n_trials = 1000 Random Walks simuliert. In der Regel werden die Parameter präziser geschätzt, wenn die Anzahl der simulierten Random Walks (bzw. allgemeiner: der stochastischen Prozesse) größer wird. Dies können Sie leicht ausprobieren, indem Sie den entsprechenden Parameter im Funktionsaufruf ändern. Probieren Sie auch andere Varianten der Startwerte aus, wenn die initiale Schätzung eher nicht gut ist.
• Die Parameterschätzung wird auch in der Regel besser, wenn die Daten genauer beschrieben werden. Dies könnte z.B. die Verwendung von mehr Quantilen sein. Auch die fehlerhaften Reaktionen können anders berücksichtigt werden: Zum einen können auch hier Quantile der First-Passage Times benutzt werden. Zum anderen ist es üblicher, sog. Conditional Accuracy Functions zu bestimmen, die den Anteil korrekter Reaktionen in separaten Quantilen der First-Passage Times angeben (siehe Abschnitt 5.1.2).

Das, was wir gerade durchgeführt haben, ist tatsächlich ein sehr wichtiger Bestandteil der Modellentwicklung und Modellevaluation, nämlich eine Ministudie zur Parameter Recovery des Modells: Für ein gutes Modell muss gezeigt werden, dass die Parameter des Modells auch gut geschätzt werden können. Dazu werden dann z.B. sehr viele synthetische Datensätze mit verschiedenen (meist zufällig gezogenen) Parametern simuliert und anschließend wird–so wie wir es gerade getan haben–das Modell “an diese Daten gefittet”. Hat das Modell gute Eigenschaften, entsprechen sich die erzeugenden und die geschätzten Parameter weitgehend, was sich über alle Datensätze hinweg dann als sehr hohe Korrelation zeigt. Dass dies nicht immer so ist und Modelle sich in dieser Eigenschaft stark unterscheiden können, zeigt ein Beispiel, welches bei Janczyk et al. (2025a) beschrieben ist.

Ein letzter Schritt wäre nun noch eine Viualisierung dessen, inwiefern die empirischen Daten und die “besten” Modellvorhersagen auch tatsächlich zu einander passen. Dazu simulieren wir nun mit den geschätzten Parametern viele Random Walks und visualisieren das Ergebnis gemeinsam mit den “empirischen” Daten (siehe Abb. 3.6). Wir sehen, dass die vorhergesagten Werte nicht perfekt mit den beobachteten Werten übereinstimmen. Allerdings sind derartige Abweichungen durchaus im Rahmen des Üblichen.

bestRW <- simulate_RW(
  t_max = 2000,
  n_trials = 50000, # Simulation mit vielen Durchgängen
  p_up = result$par[1], # bester geschätzter Wert
  boundary = result$par[2]
) # bester geschätzter Wert

### zunächst die beobachteten Daten
# nur korrekte First-Passage Times verwenden
correct_FPT <- subset(
  obsData,
  response == 1
)
# Berechnen der (beobachteten) Quantile
obsQuantiles <- quantile(correct_FPT$firstPassageTimes,
  probs = c(0.1, 0.2, 0.3, 0.5, 0.5, 0.6, 0.7, 0.8, 0.9)
)
# Plotten als schwarz-gefüllte Punkte
plot(c(0.1, 0.2, 0.3, 0.5, 0.5, 0.6, 0.7, 0.8, 0.9) ~ obsQuantiles,
  pch = 19,
  col = "black",
  ylab = "Quantile",
  xlab = "First-Passage Time",
  xlim = c(0, 120),
  ylim = c(0, 1)
)

### nun die besten vorhergesagten Daten
# nur korrekte First-Passage Times verwenden
correct_bestRW <- subset(
  bestRW,
  response == 1
)
# Berechnen der (beobachteten) Quantile
bestQuantiles <- quantile(correct_bestRW$firstPassageTimes,
  probs = c(0.1, 0.2, 0.3, 0.5, 0.5, 0.6, 0.7, 0.8, 0.9)
)
# Plotten als schwarz-gefüllte Punkte
points(c(0.1, 0.2, 0.3, 0.5, 0.5, 0.6, 0.7, 0.8, 0.9) ~ bestQuantiles,
  pch = 19,
  type = "l",
  col = "red"
)

# Legende
text(10,
  c(0.9, 0.8),
  c("beobachtet", "vorhergesagt"),
  pos = 4
)
points(7, 0.9, pch = 19)
segments(5, 0.8,
  9, 0.8,
  col = "red"
)

Abbildung 3.6: Quantile der CDF der beobachteten (simulierten) Daten sowie CDF der Vorhersage mit den besten Parameterwerten bei einem Random Walk.

3.6 Random Walks via Transitionsmatrix

Bevor wir zu Diffusionsmodellen übergehen, betrachten wir Random Walks nun noch einmal anders. Wir hatten mit der Feststellung begonnen, dass Random Walks ein Beispiel für Markov-Ketten sind. Das heißt, wir können diese prinzipiell auch als Transitionsmatrix \(P\) schreiben. Tatsächlich eröffnen sich damit weitreichende Möglichkeiten der Bestimmung von First-Passage Times und Reaktionswahlen (siehe Diederich & Busemeyer, 2003; Diederich & Mallahi-Karai, 2018, für mehr Details) und diese Methode gibt auch einen kleinen Einblick darin, wie das R-Paket dRiftDM intern arbeitet, welches wir ab dem nächsten Kapitel dann nutzen werden.

3.6.1 Transitionsmatrix und Markov-Ketten Evolution

Der Einfachheit halber nehmen wir an, dass die Boundaries bei \(b = 3\) und \(-b = -3\) liegen, und wir nehmen an, dass \(P({\mathbf Y}_i = 1) = p = 0.6\) und \(P({\mathbf Y}_i = -1) = 1 - p = 0.4\) sind. Wenn die Zeilen und Spalten die möglichen Zustände darstellen, die der Random Walk \({\mathbf X}_t \in \{-3, -2, -1, 0, 1, 2, 3\}\) annehmen kann, dann ergibt sich als Transitionsmatrix \[ P= \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0.4 & 0 & 0.6 & 0 & 0 & 0 & 0 \\ 0 & 0.4 & 0 & 0.6 & 0 & 0 & 0 \\ 0 & 0 & 0.4 & 0 & 0.6 & 0 & 0 \\ 0 & 0 & 0 & 0.4 & 0 & 0.6 & 0 \\ 0 & 0 & 0 & 0 & 0.4 & 0 & 0.6 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix}\,. \] Wie im Beispiel zu Beginn dieses Kapitels stellen die Zeilen den aktuellen Zustand zur Zeit \(t\) dar, also den Zustand, aus dem der Prozess „kommt“, und die Spalten stellen den Zustand zur Zeit \(t+1\) dar, also den Zustand, in den der Prozess übergeht. Wenn der Prozess zum Beispiel bereits absorbiert ist, das heißt, er befindet sich in \(-3\) oder \(+3\), dann bleibt er mit Wahrscheinlichkeit \(p=1\) in diesem Zustand (weil die Boundaries absorbierend sind). Ist der aktuelle Zustand \(-2\) (zweite Zeile), so beträgt die Wahrscheinlichkeit, nach \(-3\) überzugehen, \(p=0.4\) und die Wahrscheinlichkeit, nach \(-1\) überzugehen, ist \(p=0.6\). Die Zustände „zwischen“ den absorbierenden Zuständen (von \(-2\) bis \(2\) in unserem Beispiel) nennt man transiente Zustände. Um die Matrix \(P\) zu erstellen, können wir die folgende Funktion verwenden:

# Erstellung der Transitionsmatrix P mit absorbing states b and -b,
# und der Wahrscheinlichkeit p für einen Schritt nach oben 
create_transition_matrix <- function(b, p) {
  n <- 2 * b + 1 # Wie viele states gibt es?
  P <- matrix(0, # Passende Matrix voller Nullen 
    nrow = n,
    ncol = n
  )
  for (i in 2:(n - 1)) { # Transiente states
    P[i, i - 1] <- 1 - p
    P[i, i + 1] <- p
  }
  P[1, 1] <- 1 # Absorbing states
  P[n, n] <- 1
  return(P)
}

# Transitionsmatrix P anlegen
b <- 3 # Absorbing state
p <- 0.6 # p für einen Schritt nach oben
(P <- create_transition_matrix(b, p)) # anlegen und zeigen

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]  1.0  0.0  0.0  0.0  0.0  0.0  0.0
## [2,]  0.4  0.0  0.6  0.0  0.0  0.0  0.0
## [3,]  0.0  0.4  0.0  0.6  0.0  0.0  0.0
## [4,]  0.0  0.0  0.4  0.0  0.6  0.0  0.0
## [5,]  0.0  0.0  0.0  0.4  0.0  0.6  0.0
## [6,]  0.0  0.0  0.0  0.0  0.4  0.0  0.6
## [7,]  0.0  0.0  0.0  0.0  0.0  0.0  1.0

Die Kenntnis der Matrix \(P\) eröffnet nun den Zugang zu leistungsfähigeren Methoden zur Analyse von Random Walks, die dennoch auf (relativ einfacher) Matrizenalgebra basiert. Angenommen, ein Zeilenvektor \(X_0\) beschreibt den Zustand eines Random Walks zum Zeitpunkt \(t = 0\) in Form von Wahrscheinlichkeiten. Wenn der Random Walk in der Mitte zwischen zwei Boundaries \(b = 3\) und \(-b = -3\) starten soll, dann wäre dies also \[ X_0 = \begin{bmatrix}0 & 0 & 0 & 1 & 0 & 0 & 0 \end{bmatrix}\,. \] Mit anderen Worten: Der Random Walk befindet sich zum Zeitpunkt \(t=0\) mit Wahrscheinlichkeit 1 im Zustand 0. Bleiben wir bei unserem Beispiel und nehmen an, dass \(p = 0.6\) gilt. Dann erwarten wir, dass sich der Random Walk zum Zeitpunkt \(t = 1\) mit Wahrscheinlichkeit \(0.6\) im Zustand \(+1\) und mit Wahrscheinlichkeit \(0.4\) im Zustand \(-1\) befindet, also \[ X_1 = \begin{bmatrix}0 & 0 & 0.4 & 0 & 0.6 & 0 & 0 \end{bmatrix}\,. \] Das ist genau das Ergebnis, wenn \(X_0\) matrix-multipliziert wird mit \(P\): \[\begin{align*} X_1 &= X_0\cdot P \\ &= \begin{bmatrix}0 & 0 & 0 & 1 & 0 & 0 & 0 \end{bmatrix} \cdot \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0.4 & 0 & 0.6 & 0 & 0 & 0 & 0 \\ 0 & 0.4 & 0 & 0.6 & 0 & 0 & 0 \\ 0 & 0 & 0.4 & 0 & 0.6 & 0 & 0 \\ 0 & 0 & 0 & 0.4 & 0 & 0.6 & 0 \\ 0 & 0 & 0 & 0 & 0.4 & 0 & 0.6 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix} \\ &= \begin{bmatrix}0 & 0 & 0.4 & 0 & 0.6 & 0 & 0 \end{bmatrix}\,. \end{align*}\] Dies können wir in R leicht prüfen:

# X_0 anlegen
X0 <- rep(0, b * 2 + 1) # Alle Elemente als Nullen zunächst...
X0[b + 1] <- 1 # ...und dann eine 1 in die Mitte
X0

## [1] 0 0 0 1 0 0 0

# Dann X0 mit P matrix-multiplizieren
X1 <- X0 %*% P
X1

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]    0    0  0.4    0  0.6    0    0

Diese Idee generalisiert auch weiter. Um z.B. die Wahrscheinlichkeitsverteilung der States für \(t=2\) zu bestimmen, multiplizieren wir \(X_1\) noch einmal mit \(P\), \[ X_2 = X_1\cdot P = X_0\cdot P \cdot P = X_0\cdot P^2\,, \] bzw. ganz allgemein gilt für den Zeitpunkt \(t\) \[ X_t = X_0\cdot P^t\,. \] Damit können wir, ausgehend von einem initialen Zustand \(X_0\), die Wahrscheinlichkeitsverteilung zu allen Zeitpunkten \(t>0\) recht einfach berechnen. Das Paket expm stellt den Operator %^% zur Verfügung, mit dem Matrizen zur Potenz erhoben werden können (hier für \(t=10\)):

library(expm)
X10 <- X0 %*% (P %^% 10)
X10

##           [,1]       [,2] [,3]      [,4] [,5]       [,6]     [,7]
## [1,] 0.1671455 0.04299817    0 0.1289945    0 0.09674588 0.564116

3.6.2 Reaktionswahlen und First-Passage Times (ohne Simulation)

Im nächsten Schritt nutzen wir nun \(P\) um die Verteilung der First-Passage Times zu bestimmen. Weiter oben haben wir dies bereits simulationsbasiert getan. Um allerdings ein “genaues” Ergebnis zu erhalten, müssten wir unendlich viele Simulationen durchführen; idealer wäre daher natürlich eine Methode, die ohne Simulation auskommt.

Unter Verwendung der Transitionsmatrix \(P\) gehen wir dazu konzeptuell wie folgt vor:

• Wir beginnen bei \(t = 0\) mit einer Anfangsverteilung der Wahrscheinlichkeiten, zusammengefasst im Vektor \(X_0\). Diese Wahrscheinlichkeitsverteilung ist häufig auf einen einzigen Zustand in der Mitte der Boundaries konzentriert, jedoch sind auch andere Verteilungen möglich.
• Den nächsten Zustand zum Zeitpunkt \(t + 1\) bestimmen wir, indem wir den aktuellen Zustand mit der Matrix \(P\) multiplizieren, d.h., durch das zuvor eingeführte Matrixprodukt \[X_{t+1} = X_t \cdot P\,.\]
• Anschließend ermitteln wir, wie stark die Wahrscheinlichkeit, in \(-b\) oder \(b\) absorbiert zu werden, im Vergleich zum vorherigen Zustand zugenommen hat. Dies ist deshalb wichtig, weil die Wahrscheinlichkeit eines First-Passage-Ereignisses (also absorbiert zu werden) zu einem bestimmten Zeitpunkt die Wahrscheinlichkeiten früherer First-Passages (also bereits früher absorbiert worden zu sein) nicht mit einbeziehen darf. Anders ausgedrückt: Die Wahrscheinlichkeit eines First-Passage-Ereignisses zu einem bestimmten Zeitpunkt ist gerade der inkrementelle Anstieg der Gesamtwahrscheinlichkeit an der Boundary.

Formaler meinen wir: Sei \(X_0\) die Anfangsverteilung der Wahrscheinlichkeiten zum Zeitpunkt \(t = 0\), und \(P\) die zugehörige Übergangsmatrix einer Markov-Kette mit absorbierenden Boundaries bei \(-b\) und \(b\). Dann gilt für jeden Zeitpunkt \(t > 0\) \[ X_{t+1} = X_{t} \cdot P\,. \] Die Verteilungen der First-Passage Times \(f_t^+\) und \(f_t^-\) durch die Absorption in \(b\) und \(-b\) sind definiert als \[\begin{align*} f_t^+ &= P(\text{erste Absorption in }b\text{ zum Zeitpunkt }t)\\ \text{und }f_t^- &= P(\text{erste Absorption in }-b\text{ zum Zeitpunkt }t) \,. \end{align*}\]

Wir schreiben nun \[ f_t^+ = X_{t;2b+1}\text{ und }f_t^- = X_{t;1}\,, \] wobei der zweite Index auf das letzte bzw. das erste Element des Vektors \(X_t\) verweist. Dann können wir die gesuchten Wahrscheinlichkeiten leicht berechnen als \[\begin{equation*} f_t^+ = X_{t;2b+1} - X_{t-1;2b+1}\text{ und }f_t^- = X_{t;1} - X_{t-1;1} \,. \end{equation*}\]

Die Berechnung der Verteilungen mit R ist sehr einfach. Als erstes Beispiel benutzen wir die Transitionsmatrix \(P\) und den initialen Vektor \(X_0\) von oben. Die Verteilungen für die Zeitschritte \(t\in \{1,\ldots,t_{max}\}\) berechnen wir dann als Markov-Ketten Evolution und extrahieren die Verteilungen der First-Passage Times als Differenzen der Wahrscheinlichkeitsmassen an beiden absorbierenden Boundaries zwischen zwei sukzessiven Zeitschritten (mittels der R-Funktion diff()):

t_max <- 10 # Wie viele Zeitschritte?
n_states <- 2 * b + 1 # Wie viele States insgesamt?
state <- X0 # Initialer State
allStates <- state # ...wird auch die erste Zeile einer späteren Matrix 

# "Vorwärts-Berechnung" der States zu t = 1, ..., t_max
for (i in 1:t_max) {
  new_state <- state %*% P # Markov Ketten-Evolution
  allStates <- rbind(allStates, new_state) # Zu Matrix hinzufügen
  state <- new_state # Vorbereitung für nächste Iteration
}

# Extrahieren der First-Passage Wahrscheinlichkeiten in ein Dataframe
FPT <- data.frame(
  t = 1:t_max,
  FPT_lower = diff(allStates[, 1]),
  FPT_upper = diff(allStates[, n_states])
)

# Runden und Anzeigen
round(FPT, 3)

Abbildung 3.7 visualisiert die First-Passage Verteilungen für die untere und die obere Boundary noch einmal für \(b = 10\), \(p=0.55\) und \(t_{max}=300\). Hier wird die angestrebte Form der Verteilung schon sehr deutlich.

Abbildung 3.7: Visualisierung der First-Passage Time Verteilung eines Random Walks mit Drift.

Wie können wir aus diesen Informationen nun die relativen Häufigkeiten der oberen (korrekten) und unteren (fehlerhaften) Reaktionen sowie die (erwartete) First-Passage Time berechnen? Die Wahlhäufigkeiten sind zunächst sehr einfach, da sie sich aus der Summe der First-Passage Wahrscheinlichkeiten für alle Zeischritte \(t\) ergeben, getrennt für die untere und obere Boundary:

\[ P_{+b} = \sum_t P_{+b}(t) \quad \text{und} \quad P_{-b} = \sum_t P_{-b}(t) \,. \] Mit R geht das schnell so (wobei wir mit dem Beispiel der letzten Abbildung arbeiten):

prob_lower <- sum(FPT$FPT_lower) # Relative Häufigkeit unten
prob_upper <- sum(FPT$FPT_upper) # Relative Häufigkeit oben
cat("Rel. Häufigkeit unten: ", prob_lower, "; Rel. Häufigkeit oben: ", prob_upper)

## Rel. Häufigkeit unten:  0.117618 ; Rel. Häufigkeit oben:  0.8749346

Instruktiv ist auch zu sehen, dass sich beide Wahrscheinlichkeiten fast zu 1 aufaddieren:

prob_lower + prob_upper

## [1] 0.9925526

Dies ist auch der Grund, warum wir in dem Beispiel t_max <- 300 gesetzt haben. Würden wir hier “zuwenig” Zeitschritte verwenden, verpassen wir quasi einige der späteren First-Passages und insgesamt verfälscht dies dann die weiteren Informationen. Dies trifft dann insbesondere auf die erwarteten First-Passage Times zu (die ansonsten auf unendlich vielen simulierten Random Walks beruhen würden), die wir als nächstes berechnen.

Im Wesentlichen handelt es sich hier um Erwartungswerte und daher die Summe der Produkte von Zeitschritt \(t\) (zu dem eine First-Passage stattfand) und der Wahrscheinlichkeit einer First-Passage, also: \[\begin{equation*} \begin{aligned} \mathbb{E}(f^+_t) & =\sum_t t\cdot P_{+b}(t) \\ \text{und }\mathbb{E}(f^-_t) & =\sum_t t\cdot P_{-b}(t)\,. \end{aligned} \end{equation*}\] Die Berechnung mit R ergibt schnell:

fpt_lower <- sum(FPT$t * FPT$FPT_lower)
fpt_upper <- sum(FPT$t * FPT$FPT_upper)
cat("Untere FPT: ", fpt_lower, "; obere FPT: ", fpt_upper)

## Untere FPT:  8.725253 ; obere FPT:  64.90524

Mmh… während diese Berechnung eigentlich klar ist, scheint dennoch irgendetwas falsch zu sein. Weiter oben hatten wir bereits festgestellt, dass die mittleren First-Passage Times für die obere und untere Boundary gleich sind. Was hier noch fehlt ist, dass wir die Werte bedingen müssen auf diejenigen Random Walks, die tatsächlich in der oberen bzw. unteren Boundary absorbiert werden. Das heißt, wir berechnen \[\begin{equation} \begin{aligned} \mathbb{E}(f^+_t|+b) &= \frac{\sum_t t\cdot P_{+b}(t)}{\sum_t P_{+b}(t)} \\ \text{und }\mathbb{E}(f^-_t|-b) &= \frac{\sum_t t\cdot P_{-b}(t)}{\sum_t P_{-b}(t)}\,, \end{aligned} \end{equation}\] wobei der Nenner nichts anderes ist als die relative Häufigkeiten in der oberen bzw. unteren Boundary absorbiert zu werden. Damit ergibt sich:

expected_fpt_lower <- sum(FPT$t * FPT$FPT_lower) / prob_lower
expected_fpt_upper <- sum(FPT$t * FPT$FPT_upper) / prob_upper
cat("Untere FPT: ", expected_fpt_lower, "; obere FPT: ", expected_fpt_upper)

## Untere FPT:  74.18297 ; obere FPT:  74.18297

Wir halten fest:

1. Diese Methoden erlauben uns eine recht einfache und flexible Berechnung der First-Passage Time Verteilungen sowie daraus der relativen Häufigkeiten und der erwarteten First-Passage Times. Diese Werte sind nahezu “optimal”. Im folgenden Exkurs lernen wir noch eine Methode kennen, wie diese Werte in bestimmten Situationen noch einfacher und genauer berechnet werden können.
2. Die maximale Zeit \(t_{max}\), für die die Wahrscheinlichkeiten mittels Markov-Ketten Evolution berechnet werden, ist sehr wichtig. In der Praxis müssen wir sicherstellen, dass dieser Wert ausreichend groß gewählt wird, um (1) keine Reaktionen zu “verpassen” und (2) die erwarteten First-Passage Times nicht zu unterschätzen. Ab dem nächsten Kapitel werden wir mit dem R-Paket dRiftDM arbeiten und sehen, dass uns dort der Parameter t_max direkt wieder begegnet. Tatsächlich verwendet das Paket i.W. die hier dargestellten Methoden (allerdings in teilweise komplexerer, aber dafür präziserer, Form).

3.6.3 Exkurs: Bestimmung der Reaktionswahlen und First-Passage Times ohne Markov Ketten-Evolution

Es ist möglich, die erwartete First-Passage Time und die relative Häufigkeit der Reaktionswahlen basierend auf der Transitionsmatrix \(P\) und etwas Matrizenalgebra zu berechnen, ohne explizit die Verteilungen der First-Passage Times zu evaluieren. Dazu sortieren wir \(P\) etwas um, indem wir die letzte Zeile zur zweiten Zeile und die letzte Spalte zur zweiten Spalte machen. Nehmen wir das einführende Beispiel von oben, sieht das Resultat dann wie folgt aus:

\[ P= \left[ \begin{array}{cc|ccccc} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ \hline 0.4 & 0 & 0 & 0.6 & 0 & 0 & 0 \\ 0 & 0 & 0.4 & 0 & 0.6 & 0 & 0 \\ 0 & 0 & 0 & 0.4 & 0 & 0.6 & 0 \\ 0 & 0 & 0 & 0 & 0.4 & 0 & 0.6 \\ 0 & 0.6 & 0 & 0 & 0 & 0.4 & 0 \\ \end{array} \right]\,. \] Diese umstrukturierte Matrix wird manchmal auch die kanonische Form von \(P\) genannt und besteht aus vier Teilmatrizen, angezeigt durch die vertikale und die horizontale Linie:

\[ P = \left[ \begin{array}{c|c} P_I & 0 \\\hline R & Q \end{array} \right]\,. \] \(P_I\) ist dabei eine \(2\times 2\)-Einheitsmatrix (da wir zwei absorbierende States haben). Die Null-Matrix \(0\) enthält nur Wahrscheinlichkeiten von 0, weil der Random Walk sich nicht mehr weiterbewegt, sobald er in einer Boundary absorbiert wurde. Die \(5\times 2\)-Matrix \(R\) umfasst die transienten States, die im nächsten Schritt in einer absorbierenden Boundary enden würden und die \(5\times 5\)-Matrix \(Q\) enthält die Wahrscheinlichkeiten der verbleibenden transienten States.

Zusätzlich benötigen wir noch einen initialen State, also einen Vektor \(X_0\), dessen Anzahl an Elementen sich aus der Anzahl der transienten States ergibt. Wenn der Random Walk z.B. im State 0 beginnen sollen, dann verwenden wir \[ X_0= \begin{bmatrix} 0 & 0 & 1 & 0 & 0 \end{bmatrix}\,. \] Um die Matrizen \(R\) und \(Q\), die wir gleich benötigen, aus \(P\) zu extrahieren, nutzen wir eine kleine Funktion. Zusätzlich bereiten wir den Vektor \(X_0\) und die Einheitsmatrix mit der gleichen Dimension wie \(Q\) vor:

# Extrahieren der Matrizen aus P
extract_submatrices <- function(P) {
  # P umsortieren
  P <- P[c(1, nrow(P), 2:(nrow(P) - 1)), ]
  P <- P[, c(1, ncol(P), 2:(ncol(P) - 1))]

  # R and Q extrahieren
  dimP <- nrow(P)
  R <- P[3:dimP, 1:2]
  Q <- P[3:dimP, 3:dimP]

  # RÜckgabe
  return(list(R = R, Q = Q))
}

# Hier ein überschaubares Beispiel 
b <- 3
p <- 0.6
# Transitionsmatrix P erstellen
P <- create_transition_matrix(b, p)
# R und Q extrahieren
submats <- extract_submatrices(P)
(R <- submats$R)

##      [,1] [,2]
## [1,]  0.4  0.0
## [2,]  0.0  0.0
## [3,]  0.0  0.0
## [4,]  0.0  0.0
## [5,]  0.0  0.6

(Q <- submats$Q)

##      [,1] [,2] [,3] [,4] [,5]
## [1,]  0.0  0.6  0.0  0.0  0.0
## [2,]  0.4  0.0  0.6  0.0  0.0
## [3,]  0.0  0.4  0.0  0.6  0.0
## [4,]  0.0  0.0  0.4  0.0  0.6
## [5,]  0.0  0.0  0.0  0.4  0.0

# X0 und I anlegen
(X0 <- c(rep(0, times = b - 1), 1, rep(0, times = b - 1)))

## [1] 0 0 1 0 0

(I <- diag(nrow(Q)))

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    0    0    0    0
## [2,]    0    1    0    0    0
## [3,]    0    0    1    0    0
## [4,]    0    0    0    1    0
## [5,]    0    0    0    0    1

Diese Matrizen können wir nun sehr leicht nutzen, um die relativen Häufigkeiten der Reaktionswahlen ud die erwarteten First-Passage Times zu berechnen (siehe auch Diederich & Busemeyer, 2003; Diederich & Mallahi-Karai, 2018). Wir beginnen mit den Reaktionswahlen. Die \(2\times 1\)-Matrix \(C\) wird berechnet als \[\begin{equation*} C=\begin{bmatrix}P(\text{"unten"}) & P(\text{"oben"})\end{bmatrix}= X_0\cdot (I-Q)^{-1}\cdot R\,, \end{equation*}\] wobei das erste Element \(C_{1,1}\) die relative Häufigkeit einer “unteren” Reaktion enthält und das zweite Element \(C_{1,2}\) enthält die relative Häufigkeit einer “oberen” Reaktion. Dabei ist \((I-Q)^{-1}\) die Inverse der Matrix \((I-Q)\) mit \((I-Q)\cdot (I-Q)^{-1}=I\). Die Berechnung mit R ist sehr einfach, da wir die Inverse einer Matrix mit der R-Funktion solve() berechnen können:

# Berechnung der Reaktionswahlen
(C <- X0 %*% solve(I - Q) %*% R)

##           [,1]      [,2]
## [1,] 0.2285714 0.7714286

Nun folgen die (bedingten) erwarteten First-Passage Times. Wenn \(\mathbf{T}\) als Zufallsvariable die First-Passage Times beschreibt, dann wird die \(2\times 1\)-Matrix \(\mathbb{E}_\mathbf{T}\) berechnet als \[\begin{equation*} \begin{aligned} \mathbb{E}_\mathbf{T} &= \begin{bmatrix}\mathbb{E}(t|\text{"lower"}) & \mathbb{E}(t|\text{"upper"})\end{bmatrix} \\ &= [X0\cdot(I-Q)^{-2}\cdot R] \oslash C \,. \end{aligned} \end{equation*}\] Hierbei ist \((I-Q)^{-2}\) das Matrixprodukt \((I-Q)^{-1}\cdot (I-Q)^{-1}\) und \(\oslash\) meint eine Hadamard Division, d.h., eine element-weise Division der Matrizen. Mit R kann die Berechnung wir folgt durchgeführt werden:

# Berechnung der erwarteten First-Passage Times
(E_T <- (X0 %*% (solve(I - Q) %*% solve(I - Q)) %*% R) / C)

##          [,1]     [,2]
## [1,] 8.142857 8.142857

Nun wenden wir diese Methode noch einmal auf das Beispiel von oben (siehe Abb. 3.7) an:

b <- 10 # Absorbing state
p <- 0.55 # p für einen Schritt nach oben (+1)
t_max <- 300

# Transitionsmatrix P erstellen
P <- create_transition_matrix(b, p)
# R und Q extrahieren
submats <- extract_submatrices(P)
R <- submats$R
Q <- submats$Q
# X0 und I anlegen
X0 <- c(rep(0, times = b - 1), 1, rep(0, times = b - 1))
I <- diag(nrow(Q))

# Berechnung der Reaktionswahlen
(C <- X0 %*% solve(I - Q) %*% R)

##           [,1]      [,2]
## [1,] 0.1185005 0.8814995

# Berechnung der erwarteten First-Passage Times
(E_T <- (X0 %*% (solve(I - Q) %*% solve(I - Q)) %*% R) / C)

##          [,1]     [,2]
## [1,] 76.29989 76.29989

Zum Vergleich hier noch einmal die Werte von oben:

## Rel. Häufigkeit unten:  0.117618 ; Rel. Häufigkeit oben:  0.8749346

## Untere FPT:  74.18297 ; obere FPT:  74.18297

Hier können wir sehen, dass die oben berechneten Werte (sowohl für die Reaktionswahlen als auch für die erwarteten First-Passage Times) die “wahren” Werte leicht unterschätzen. Dies rührt daher, dass auch bei \(t_{max} = 300\) einige der First-Passages zu späteren Zeitpunkten nicht erfasst werden. Je höher \(t_{max}\) gewählt wird, desto mehr nähern sich diese Werte allerdings an.

3.7 Zusammenfassung Random Walk Modelle

In diesem Kapitel haben wir mit Random Walks stochastische Prozesse kennengelernt, die eine erste Annäherung an die Modellierung von Reaktionsauswahl bieten. Wir haben dabei mit dem einfachen, symmetrischen Random Walk begonnen, der auch gleich im Kontext der Diffusionsmodelle eine wichtige Rolle spielen wird. Anhand von Random Walks mit Drift und mit absorbierenden Boundaries konnten wir dann Verteilungen von First-Passage Times bestimmen und diese als Grundlage zur Berechnung der Reaktionswahlen und der (erwarteten) First-Passage Times nutzen.

Ein Nachteil der bisher behandelten Random Walks ist, dass sie sowohl im state space als auch in der Zeit jeweils diskret sind, also nur bestimmte Werte annehmen können. Die Differenz jener Werte kann auch zusammengefasst werden mit \(\Delta t=1\) und \(\Delta x = 1\). Diffusionsmodelle sind quasi die kontinuierliche Version eines Random Walks, d.h., es können alle potenziellen Zwischenschritte angenommen werden.

4 Diffusionsmodelle

Wir haben in Kapitel 3 ausführlich das Prinzip des Random Walks beschrieben, einschließlich einer Anwendung des Grundprinzips der Parameterschätzung für solche Modelle. Wenngleich Random Walks schon recht gut reale First-Passage Times vorhersagen können, kommen sie auch an ihre Grenzen. Wegen ihres Grundprinzips des Sammelns von Evidenz sind Random Walks auch ein Beispiel für die größere Klasse sog. Evidenzakkumulationsmodelle (engl. evidence accumulation models).

Wir werden in diesem Kapitel das, in der Psychologie, vielleicht am weitesten verbreitete Modell dieser Klasse kennen lernen: das Diffusionsmodell. Nach einer grundlegenden Einführung und einer kurzen Exploration des Modellverhaltens, kommen wir dann auch auf Anwendungen zu sprechen. Im folgenden Kapitel 5 greifen wir dann speziellere Formen von Diffusionsmodellen wieder auf.

Im Random Walk Modell aus Kapitel 3 waren die Zeitschritte ganzzahlig und ebenso galt dies für die möglichen Stateversetzungen (um +1 oder -1). In anderen Worten war der Unterschied im state space immer \(\Delta x = 1\) und die betrachtete Zeitschrittweite war \(\Delta t =1\). Das Modell war also sowohl in der Zeit als auch im State diskret. Der Übergang zum Diffusionsmodell besteht nun daraus, den Fall zu untersuchen wenn sowohl \(\Delta x\rightarrow 0\) als auch \(\Delta t \rightarrow 0\) gilt. Das heißt, wir betrachten den Fall, dass sowohl Zeit als auch state space kontinuierlich sind und die entsprechenden stochastischen Prozesse werden dann Diffusionsprozesse genannt.

Der Übergang mag auf den ersten Blick verkomplizierend wirken, hat aber zwei entscheidende Vorteile:

• Zum einen betrachten wir damit einen generellen Rahmen von Modellierungsmöglichkeiten, der uns mehr Flexibilität in Bezug auf Anwendungen und Erweiterungen erlaubt.
• Zum anderen ist in vielen Fällen das Modellverhalten mit Standardmethoden der Mathematik, wie z.B. (partiellen) Differentialgleichungen gut untersuchbar und auch verstanden. Wir werden hier wenig in die Tiefe gehen, sondern eher auf die Konzepte und die praktische Verwendung des R-Pakets dRiftDM (v. 0.2.2.9000) fokussieren. Allerdings nutzen die dahinterstehenden Methoden des öfteren kompliziertere mathematische Verfahren.

Wir beginnen nun damit, den unsystematischen, stochastischen Anteil eines Diffusionsprozesses herzuleiten bevor wir dann später wieder einen Drift hinzufügen (weswegen auch vom Drift-Diffusionsmodell gesprochen wird).

4.1 Der stochastische Anteil: Brown’sche Bewegung

Der unsystematische oder stochastische Anteil eines Diffusionsprozesses wird i.d.R. durch eine Brown’sche Bewegung modelliert, die zufällige und ruckartige Bewegungen beschreibt. Ursprünglich vom Botaniker Robert Brown anhand der Bewegung kleiner Teilchen in Gasen und Flüssigkeiten entdeckt, wurde das mathematische Modell vom Mathematiker Norbert Wiener beschrieben und in der Folge auch als Wiener-Prozess bezeichnet. Etwas vereinfacht gesagt ergeben sich die Zuwächse von jedem \(t\) zu \(t+\Delta t\) durch stochastisch unabhängige normalverteilte Zufallsvariablen.

4.1.1 Definition der Brown’schen Bewegung bzw. eines Wiener-Prozesses

Eine Brown’sche Bewegung \(B(t)\) ist für \(t\geq 0\) üblicherweise über die folgenden vier Eigenschaften definiert:

1. Die initiale Bedingung ist \(B(0)=0\).
2. Für alle \(0\leq s < t\) sind die Inkremente \(B(t)-B(s)\) unabhängig von ihrer Vergangenheit, d.h. von \(B(u), \forall u \leq s\).
3. Die Verteilung von \(B(t)-B(s)\) hängt nicht vom absoluten Zeitpunkt ab, sondern nur von \(t-s\) und ist \[B(t)-B(s)\sim\mathcal{N}(0,\sigma^2(t-s))\,,\] mit \(\sigma^2>0\). Hierbei ist zu beachten, dass \(\sigma^2\) nur eine Skalierungskonstante ist, die die Varianz der Brown’schen Bewegung bestimmt.
4. Die Funktion \(t\rightarrow B(t)\) ist kontinuierlich.

Aus dieser Definition folgt \[\begin{equation*} B(t)\sim \mathcal{N}(0,\sigma^2\cdot t)\,. \end{equation*}\] Wenn \(\sigma^2 = 1\) gewählt wird, wird oft von einer Standard Brown’schen Bewegung gesprochen.

Bemerkenswert ist ürigens, dass die Funktion \(t\rightarrow B(t)\) nicht differenzierbar ist, obwohl sie kontinuierlich ist: Das liegt daran, dass sie immer unregelmäßig und zufällig bleibt, auch wenn man beliebig weit in sie “herein-zoomt” (sie ist in dieser Hinsicht daher ähnlich zur Weierstrass Funktion).

Auch wenn dies alles abstrakt klingt, hilft es, sich die Brown’sche Bewegung als die kumulative Summe normalverteilter Zufallsvariablen vorzustellen. Tatsächlich reicht uns am Ende eine Approximation, die als Grenzwert in die Brown’sche Bewegung übergeht: Computer können generell nicht wirklich mit kontinuierlichen Variablen umgehen und wir müssen stattdessen wieder die Zeit diskretisieren, d.h., die Software arbeitet ohnehin mit diskreten Varianten. Wir müssen aber sicherstellen, dass die verwendete Approximation Eigenschaften der Brown’schen erfüllt, wenn die Diskretisierung feiner und feiner wird. Dazu werden wir nun auf ein mathematisches Modell hinarbeiten, welches die wichtigen Eigenschaften erfüllt. Dabei sehen wir, dass wir \(\Delta t\) und \(\Delta x\) nicht unabhängig voneinander betrachten dürfen. Dies zu sehen, ist die wichtige Erkenntnis des nächsten längeren Abschnitts.

4.1.2 Herleitung aus dem einfachen, symmetrischen Random Walk

Beim Random Walk war der state space \(S = \mathbb{Z}\) und die Zeit war definiert auf \(T=\mathbb{N}_0\). Da der Random Walk auch immer nur um \(1\) nach oben oder unten versetzt werden konnte, war \(\Delta x=1\). Ganz ähnlich waren die Zeitschritte immer separiert um \(\Delta t=1\). Nun zerlegen wir diese grobe Bewegung in viele kleine Schritte, und dies sowohl im State als auch in der Zeit, und betrachten schließlich den Fall, dass \(\Delta t \rightarrow 0\) und \(\Delta x\rightarrow 0\). Wir werden sehen, dass wir dabei beide Größen nicht unabhängig voneinander betrachten können, sondern die eine Größe von der anderen abhängt. Allerdings wird die Bewegung mit kleiner werdenden Zeitschritten immer glatter, aber immer noch zufällig und ohne klare Richtung. Das Ergebnis ist dann eine Standard Brown’sche Bewegung. Weil Markov-Prozesse die kontinuierlich in State und Zeit verlaufen auch Diffusionsprozesse genannt werden, ist dies dann die einfachste Form eines solchen Diffusionsprozesses.

Wir beginnen nun damit, das interessante Zeitintervall \([0,t_{max}]\) in \(n\)-viele kleinere Schritte der Größe \(\Delta t = \frac{t_{max}}{n}\) zu zerlegen. Dann modifizieren wir den einfachen, symmetrischen Random Walk so, dass er seine “Sprünge” zu den Zeitpunkten \(t \in \{0\cdot \Delta t, 1\cdot \Delta t, 2\cdot \Delta t, \ldots, n\cdot \Delta t \}\) (mit \(n\in\mathbb{N}\)) macht. In anderen Worten, wir betrachten die Zeitschritte nun als Vielfache von \(\Delta t\), d.h., als \(k\cdot \Delta t\), \(0\leq k \leq n\). In R könnten diese Zeitschritte sehr leicht angelegt werden:

t_max <- 10 # Maximale Zeit (z.B. in ms)
# 1) Mit dt = 1 (wie vorher)
dt <- 1.0
(t <- seq(0, t_max, by = dt))

##  [1]  0  1  2  3  4  5  6  7  8  9 10

# 2) Kleinere Zeitschritte, z.B. dt = 0.5
dt <- 0.5
(t <- seq(0, t_max, by = dt))

##  [1]  0.0  0.5  1.0  1.5  2.0  2.5  3.0  3.5  4.0  4.5  5.0  5.5  6.0  6.5  7.0
## [16]  7.5  8.0  8.5  9.0  9.5 10.0

Wir lassen außerdem zu, dass die Inkremente nach oben oder unten allgemein der Größe \(+\Delta x\) bzw. \(-\Delta x\) sind (statt nur \(1\) oder \(-1\)). Der state space \(S\) wird dadurch

\[S=\{\ldots,-2\cdot\Delta x, -1\cdot\Delta x, 0, + 1\cdot \Delta x,+2\cdot \Delta x,\ldots\}\,.\] Die stochastische Differenzgleichung dieses Random Walks ist dann \[\begin{equation*} \mathbf{X}(t+\Delta t) = \mathbf{X}(t) + \Delta x\cdot \mathbf{Y}_t\text{ mit }\mathbf{X}(0)=0 \,, \end{equation*}\] und sie beschreibt wieder die Evolution des Prozesses in der Zeit. Da wir auf einen kontinuierlichen Prozess hinarbeiten, schreiben wir nun \(\mathbf{X}(t)\), \(\mathbf{X}(t+\Delta t)\), …, auch wenn in manchen Büchern hierfür auch \(\mathbf{X}_t\), \(\mathbf{X}_{t+\Delta t}\), … geschrieben würde. Wie vorher auch sind die \(\mathbf{Y}_i\) unabhängige Zufallsvariablen mit \(P(\mathbf{Y}_i = +1) = P(\mathbf{Y}_i=-1) = 0.5\).

Für die nächsten Schritte ist es instruktiv, sich dies für einige Schritte von \(k\) explizit aufzuschreiben: \[\begin{equation} \begin{aligned} \mathbf{X}(0) &= 0 \\ \mathbf{X}(1\cdot\Delta t) &= \mathbf{X}(0) + \Delta x\cdot \mathbf{Y}_1 \\ \mathbf{X}(2\cdot\Delta t) &= \mathbf{X}(1\cdot \Delta t) + \Delta x\cdot \mathbf{Y}_2 \\ &= \mathbf{X}(0) + \Delta x\cdot \mathbf{Y}_1 + \Delta x\cdot \mathbf{Y}_2 \\ &= \Delta x\cdot (\mathbf{Y}_1+\mathbf{Y}_2)\\ \mathbf{X}(3\cdot\Delta t) &= \mathbf{X}(2\cdot \Delta t) + \Delta x\cdot \mathbf{Y}_3 \\ &= \mathbf{X}(0) + \Delta x\cdot \mathbf{Y}_1 + \Delta x\cdot \mathbf{Y}_2 + \Delta x\cdot \mathbf{Y}_3\\ &= \Delta x\cdot (\mathbf{Y}_1+\mathbf{Y}_2+\mathbf{Y}_3)\\ \vdots \quad&\phantom{=} \quad\quad\vdots \\ \mathbf{X}(k\cdot\Delta t) &= \mathbf{X}((k-1)\cdot \Delta t) + \Delta x\cdot \mathbf{Y}_k \\ &= \mathbf{X}(0) + \Delta x\cdot \mathbf{Y}_1 + \Delta x\cdot \mathbf{Y}_2 + \Delta x\cdot \mathbf{Y}_3+\ldots+\Delta x\cdot \mathbf{Y}_k\\ &= \Delta x\cdot (\mathbf{Y}_1+\mathbf{Y}_2+\mathbf{Y}_3+\ldots+\mathbf{Y}_k)\\ &= \Delta x\cdot \sum_{i=1}^k \mathbf{Y}_i\,. \end{aligned} \end{equation}\] Da die Zeit beim Schritt \(k\) als \(t=k\cdot \Delta t \Leftrightarrow k = \frac{t}{\Delta t}\) geschrieben werden kann, können wir die letzte Gleichung auch schreiben als \[\begin{equation*} \begin{aligned} \mathbf{X}(k\cdot \Delta t) = \mathbf{X}(t) &= \Delta x \cdot (\mathbf{Y}_1+\mathbf{Y}_2+\mathbf{Y}_3+\ldots+\mathbf{Y}_{\frac{t}{\Delta t}}) \\ &= \Delta x \cdot \sum_{i = 1}^{t /\Delta t}\mathbf{Y}_i\,. \end{aligned} \end{equation*}\] Was sind nun der Erwartungswert und die Varianz von \(\mathbf{X}(t)\)? Aus Kapitel 3 erinnern wir uns, dass \(\mathbb{E}(\mathbf{Y}_i)=0\) war. Damit können wir sehr leicht den Erwartungswert bestimmen als \[\begin{equation*} \begin{aligned} \mathbb{E}[\mathbf{X}(t)] &= \mathbb{E}\left[\Delta x \cdot \sum_{i = 1}^{t /\Delta t}\mathbf{Y}_i\right] \\ &= \Delta x \cdot \sum_{i = 1}^{t /\Delta t}\mathbb{E}(\mathbf{Y}_i) \\ &= 0\,. \end{aligned} \end{equation*}\] Die Varianz berechnet sich nun als \[\begin{equation*} \begin{aligned} V[\mathbf{X}(t)] &= V\left[ \Delta x \cdot \sum_{i = 1}^{t /\Delta t}\mathbf{Y}_i \right]\\ &= (\Delta x)^2\cdot V\left[\sum_{i = 1}^{t /\Delta t}\mathbf{Y}_i \right] \\ &= (\Delta x)^2\cdot \sum_{i = 1}^{t /\Delta t}V(\mathbf{Y}_i) \\ &= (\Delta x)^2\cdot \frac{t}{\Delta t}\,, \end{aligned} \end{equation*}\]

was wir zunächst im Kopf behalten. Nützlich ist nun, dass \(\mathbf{X}(t)\) die Summe vieler unabhängiger und identisch verteilter Zufallsvariablen \(\mathbf{Y}_i\) ist. Gemäß des Zentralen Grenzwertsatzes konvergiert diese Summe daher gegen eine Normalverteilung wenn \(\Delta t \rightarrow0\). Damit bleibt festzuhalten, dass \[\begin{equation*} \mathbf{X}(t) \overset{app.}{\sim} \mathcal{N}\left(0, (\Delta x)^2\cdot \frac{t}{\Delta t}\right)\,. \end{equation*}\] Denken wir nun zurück an die Definition der Brown’schen Bewegung, sehen wir, warum dies nützlich ist: \(\mathbf{X}(t)\) folgt genau wie \(B(t)\) einer Normalverteilung zu jedem Zeitpunkt \(t\)!

Abbildung 4.1 illustriert, dass für eine wachsende Anzahl simulierter \(\mathbf{Y}_i\), deren Summe gegen eine Normalverteilung konvergiert.

Abbildung 4.1: Illustration, dass \(\mathbf{X}(t)\) gegen eine Normalverteilung konvergiert, je mehr \(\mathbf{Y}_i\) in die Summe \(\mathbf{X}(t)\) eingehen. Die schwarzen Kurven viualisieren die theoretisch erwarteten Normalverteilungen. Für jede Abbildung wurden 1000 \(\mathbf{X}(t)\) realisiert.

Der offene Punkt aus der obigen Herleitung ist nun noch, dass die Varianz von \(\mathbf{X}(t)\) von \(\Delta x\) und \(\Delta t\) abhängt. An den x-Achsen von Abbildung 4.1 sehen wir auch, dass die Varianz von \(\mathbf{X}(t)\) größer wird, je mehr \(\mathbf{Y}_i\) in die Summe eingehen, d.h., je kleiner \(\Delta t\) gewählt wird (\(\Delta x\) war in der Simulation fixiert).

Damit kommen wir zum letzten Teil: Wir wollen zwar \(\Delta t \rightarrow 0\) und \(\Delta x\rightarrow 0\), aber in einer derartigen Weise, dass die Varianz von \(\mathbf{X}(t)\) konstant zur Varianz aus der oben definierten Brown’schen Bewegung ist. In anderen Worten, wollen wir, dass \[ V[\mathbf{X}(t)] = (\Delta x)^2\cdot\frac{t}{\Delta t}=\sigma^2\cdot t = V[B(t)] \] ist. Lösen wir dies nach \(\Delta x\) auf, erhalten wir \[ \begin{aligned} & & \; (\Delta x)^2\cdot\frac{t}{\Delta t} &= \sigma^2\cdot t \\ & \Leftrightarrow & \; \frac{(\Delta x)^2}{\Delta t} &= \sigma^2 \\ & \Leftrightarrow & \; (\Delta x)^2 &= \sigma^2\cdot \Delta t \\ & \Rightarrow & \; \Delta x &= \sigma\sqrt{\Delta t} \end{aligned} \] um die Varianz \(V[\mathbf{X}(t)]\) konstant zur Varianz einer Brown’schen Bewegung zu halten: \[ \begin{aligned} V[\mathbf{X}(t)] &= (\Delta x)^2\cdot\frac{t}{\Delta t} \\ &= (\sigma\sqrt{\Delta t})^2\cdot\frac{t}{\Delta t} \\ &= \sigma^2\cdot t \\ &= V[B(t)]\,. \end{aligned} \] Abbildung 4.2 visualisiert die Trajektorien approximierter Brown’scher Bewegungen für verschiedene Werte von \(\Delta t\). Je kleiner dieser Wert gewählt wird, desto glatter wird die Bewegung, wenngleich sie immer zufällig und nicht-deterministisch verläuft.

Abbildung 4.2: Illustration verschiedener Brown’scher Bewegungen mit verschiedenen Werten für \(\Delta t\) (\(\sigma=1\) wurde in allen Fällen verwendet).

Abbildung 4.3 illustriert die Sinnhaftigkeit der oben abgeleiteten Skalierung noch einmal anders. Dort sind 1000 simulierte Trajektorien für \(\Delta t = 1\) (oben) und \(\Delta t = 0.1\) (unten; jeweils mit \(\sigma=1\)) abgebildet. Zu den Zeitpunkten \(t=15\) und \(t=30\) wurde dann die Varianz der Werte auf \(\mathbf{X}(t)\) berechnet. Diese Werte sind nahe an den theoretisch erwarteten Werten 15 und 30 (und unterliegen auch den einzeichneten Normalverteilungen).

Abbildung 4.3: Illustration der Verteilung und Varianz von \(\mathbf{X}(t)\) für zwei verschiedene Werte von \(\Delta t\).

Insgesamt bleibt also festzuhalten, dass die Skalierung \(\Delta x=\sigma\sqrt{\Delta t}\) den Zweck erfüllt, \(V[\mathbf{X}(t)]\) korrekt zur theoretischen Varianz einer Brown’schen Bewegung zu halten, und dies für verschiedene Werte von \(\Delta t\). Was passiert, wenn \(\Delta x\) falsch gewählt wird?

• Wird \(\Delta x\) zu groß gewählt, dann geht \(V[\mathbf{X}(t)]\rightarrow \infty\) wenn \(\Delta t\rightarrow 0\), d.h., die Varianz “explodiert” .
• Wenn \(\Delta x\) zu klein gewählt wird, dann geht \(V[\mathbf{X}(t)]\rightarrow 0\) wenn \(\Delta t\rightarrow 0\), d.h., die Varianz verschwindet und der Prozess konvergiert gegen einen konstanten Wert.

Nun können wir also eine Brown’sche Bewegung adäquat approximieren. Natürlich wird dies nie eine “echte” Brown’sche Bewegung werden, da dies unendlich kleine Zeitschritte \(\Delta t\) benötigen würde, was ein Computer nicht kann.

Bisher hat unser Diffusionsprozess zudem noch keinerlei systematische Auf- oder Abwärtsbewegung, da ja \(\mathbb{E}[\mathbf{X}(t)]=0\) ist. Genau wie wir in Kapitel 3 eine Tendenz nach oben oder unten in einen Random Walk eingebaut haben, wenn \(p\neq 0.5\) war, implementieren wir im nächsten Abschnitt eine Driftrichtung in die Brown’sche Bewegung.

4.2 Brown’sche Bewegung mit Drift

Die Brown’sche Bewegung mit Drift wird aus einer Standard Brown’schen Bewegung abgeleitet, indem sie mit einer (nicht-stochastischen) Funktion \(m(t)\) superimponiert wird (d.h., beide werden addiert). Während prinzipiell \(m(t)\) frei gewählt werden kann, ist die typische Form eine lineare Funktion \[m(t)=\mu\cdot t\,,\] die dann zur Brown’schen Bewegung addiert wird (manchmal mit einem initialen Startwert \(x_0\), den wir vorerst bei \(x_0=0\) belassen): \[\begin{equation} \mathbf{X}(t)= x_0 + \mu\cdot t + \sigma B(t) = \mu\cdot t + \sigma B(t)\,. \end{equation}\] In dieser Gleichung wird \(\mu\) typischerweise als Drift Rate und \(\sigma\) als Diffusionskonstante bezeichnet. Dies ist quasi auch die einfachste Form eines Diffusionsmodells, wie es von Ratcliff (1978) verwendet wurde und das, was in der Psychologie, oft als Diffusionsprozess/-modell verstanden wird.

Mit \(x_0=0\) ist der Erwartungswert von \(\mathbf{X}(t)\) nun eine lineare Funktion von \(t\), also \[ \mathbb{E}[\mathbf{X}(t)] = \mathbb{E}[x_0 + \mu\cdot t + \sigma B(t)] = \mu\cdot t\,. \] Dies liegt natürlich daran, dass der Erwartungswert der Brown’schen Bewegung Null ist. Die Varianz verändert sich durch den addierten Drift nicht und sie ist daher immer noch \[ V[\mathbf{X}(t)] = \sigma^2\cdot t\,. \] Anders herum betrachtet, ist die Drift Rate nichts anders als die erste Ableitung des erwarteten Verlaufs der akkumulierten Evidenz nach der Zeit \(t\) (siehe auch Cox & Miller, 1965; Schwarz, 2022), also, \[ \frac{d\mathbb{E}[\mathbf{X}(t)]}{dt}=\frac{d(\mu\cdot t)}{dt}=\mu\,. \] In diesem Beispiel ist die Drift Rate \(\mu\) also unabhängig von der Zeit \(t\) und damit konstant über die Zeit. Daher wird dieses einfache Modell auch zeit-unabhängig genannt. Die Drift Rate indiziert also die Steigung des erwarteten Zeitverlaufs der Akkumulation und sie ist daher die zeitpunktbezogene Änderungsrate (die aber hier für alle Zeitpunkte eben identisch ist). Beispiele für Brown’sche Bewegungen, erwartete Zeitverläufe und die resultierenden Diffusionsprozesse sind in Abbildung 4.4 visualisiert.

Abbildung 4.4: Illustration von Brown’schen Bewegungen (rot), erwarteten Zeitverläufen der Evidenzakkumulation (blau) und den resultierenden Diffusionsprozessen (schwarz).

Insgesamt kann ein Diffusionsprozess also geschrieben werden als \[\mathbf{X}(t)=\mu\cdot t+B(t)\cdot\sigma\,,\] wobei \(\mu\cdot t\) auch als der systematische Anteil bezeichnet wird, und \(B(t)\) der unsystematische/stochastische Anteil ist. Um das Verhalten des Diffusionsprozesses für einen Zeitschritt \(\Delta t\) zu approximieren, können wir dies explizit als stochastische Differenzgleichung herleiten und gelangen damit zur diskreten Approximation des eigentlich kontinuierlichen Prozesses: \[ \begin{aligned} & & \mathbf{X}(t+\Delta t)-\mathbf{X}(t) &= \mu\cdot \Delta t + \sigma\cdot(\underset{\Delta B(t)\sim\mathcal{N}(0, \Delta t)}{\underbrace{B(t+\Delta t) - B(t)}}) \\\\ & \Leftrightarrow & \mathbf{X}(t+\Delta t) &= \mathbf{X}(t) + \mu\cdot \Delta t + \sigma\cdot \Delta B(t) \\\\ & \Rightarrow & \mathbf{X}(t+\Delta t) &= \mathbf{X}(t) + \mu\cdot \Delta t + \sigma\cdot\sqrt{\Delta t}\cdot \mathbf{Z}_t \end{aligned} \] mit \(\mathbf{Z}_t\sim\mathcal{N}(0,1)\) und \(\mathbf{X}(0)=0\).

Schließlich, aber hier ohne Herleitung (siehe hier für mehr Informationen), können wir auch eine andere Notation betrachten, die den kontinuierlichen Charakter des Diffusionsprozesses hervorhebt. Wenn wir \(\Delta t\rightarrow 0\) gehen lassen, gelangen wir zur stochastischen Differentialgleichung \[\begin{equation*} d\mathbf{X}(t) = \mu dt +\sigma dB(t)\,. \end{equation*}\] Diese Gleichung bedeutet grob, dass der Prozess \(\mathbf{X}(t)\) sich in jedem (unendlich) kleinen Zeitschritt \(\Delta t\) um einen systematischen, festen Anteil \(\mu dt\) und um einen stochastischen Anteil \(\sigma dB(t)\) ändert, der aus einer Brown’schen Bewegung stammt. Dieser letzte Term ist allerdings nicht im Sinne klassischer Differentialrechnung zu verstehen, sondern im Sinne von It\(\overline{\text{o}}\)-Integration, da die Funktion nicht glatt verläuft. It\(\overline{\text{o}}\)-Integration ist quasi der stochastische Gegenpart zur Riemann-Integration. Insgesamt ist diese Gleichung eine Kurzform dessen, wie eine Brown’sche Bewegung mit Drift geschrieben werden und daher wiederum das einfachste Diffusionsmodell.

4.3 Boundary und Residualzeit

Schließlich führen wir noch, wie beim Random Walk Modell auch schon, eine Boundary ein: Überschreitet der Diffusionsprozess zuerst die obere Boundary, also \(\mathbf{X}(t)>b\), dann gilt die obere Reaktion als ausgewählt; unterschreitet er die untere Boundary, also \((\mathbf{X}(t)<-b)\), dann gilt die untere Reaktion als ausgewählt. Die Zeit \(t\) zu der dies passiert wird auch wieder First-Passage Time genannt. Aufgrund der stochastischen Komponente erreicht der Diffusionsprozess bei gleichbleibendem Wert \(\mu\) eine Boundary immer zu einem anderen Zeitpunkt, wodurch Varianz in der First-Passage Time erzeugt wird. Zusätzlich kann es auch passieren, dass der Diffusionsprozess durch Zufall bei \(\mu>0\) zuerst den Zustand \(\mathbf{X}(t)<-b\) und bei \(\mu<0\) zuerst den Zustand \(\mathbf{X}(t)>b\) erreicht. Dies würde dann wiederum einer falschen Reaktion entsprechen.

Häufig wird, vor allem im hier betrachteten Kontext, davon ausgegangen, dass der Betrag der Drift Rate \(\mu\) für die Stimuli identisch ist. Dann wird i.d.R. ein positiver Wert für \(\mu\) angenommen und eine “obere” Reaktion ist dann eine korrekte Reaktion, während eine untere Reaktion einen Fehler darstellt.

Anmerkung: Wir haben in Kapitel 3 die Wald-Verteilung als eine mögliche Verteilung zur Beschreibung von RTs kennengelernt. Würde man nur von einer oberen Boundary ausgehen, die jeder Diffusionsprozess bei \(\mu>0\) irgendwann erreicht (“single absorbing barrier”), würde dies tatsächlich eine Wald-Verteilung erzeugen. Sie wird bei Feller (1968) hergeleitet und eine Methode zur Approximation wird bei Blurton et al. (2012) beschrieben. Bei zwei Boundaries, so wie wir sie benutzen, ist die Verteilung der First-Passage Times allerdings keine Wald-Verteilung mehr.

Damit haben wir im Prinzip die Komponenten eines Diffusionsmodells beisammen, mit welchen die Reaktionswahlen und ihre RTs modelliert werden. Was wir bisher noch nicht bedacht haben ist, dass eine “echte” RT auch noch die Dauern der Wahrnehmung und der Motorik umfasst. Um diesen Anteil abzudecken, wird zur First-Passage Time noch eine Residualzeit \(t0\) addiert (auch non-decision time oder \(t_{\text{er}}\) (für encoding und responding) genannt).

Die Residualzeit \(t0\) kann natürlich einfach eine Konstante sein, um die die Verteilung der First-Passage Times verschoben wird. Häufiger wird allerdings angenommen, dass die \(t0\) einer Verteilung folgt, z.B., dass sie gleichverteilt sei, \[t0\sim U(t0, S_{t0})\,,\] wobei sich das Minimum und Maximum der Gleichverteilung dann als \(t0 \pm S_{t0}\) ergibt. In anderen Fällen wird auch eine Normalverteilung angenommen, also \(t0\sim \mathcal{N}(t0,s_{to})\), letztlich scheint die exakte Verteilung aber keine allzu gravierende Rolle zu spielen (siehe dazu Ratcliff, 2013).

Mathematisch folgt die PDF der RTs dann aus einer Faltung zweier (unabhängiger) Zufallsvariablen bzw. deren PDFs, in diesem Fall die PDFs der First-Passage-Time und der Residualzeit. Für jeden möglichen Gesamtwert (also für jede mögliche RT) summiert die Faltung quasi alle Kombinationen auf, bei denen die beiden Zufallsvariablen zusammen genau diesen Wert ergeben. Die resultierende PDF gibt an, wie wahrscheinlich jede dieser möglichen Summen ist (genauer: sie beschreibt die Wahrscheinlichkeit, dass Werte in einem bestimmten Intervall realisiert werden; also die übliche Interpretation einer Dichtefunktion). Abbildung 4.5 veranschaulicht dies für den Fall, dass \(t0\) einer Normalverteilung folgt. Offensichtlich ist die PDF der RTs im Vergleich zur PDF der First-Passage-Times nach rechts verschoben, außerdem ist sie etwas breiter und steigt glatter an.

Abbildung 4.5: Illustration von PDFs der First-Passage Times und einer normalverteilten Residualzeit \(t0\) (links) sowie deren Faltung (rechts). FPT = First-Passage Time.

4.4 Praktische Umsetzung eines Diffusionsmodells

Nun wollen wir uns mehr mit der praktischen Verwendung eines Diffusionsmodells beschäftigen und sein Verhalten kennenlernen. Aus Übungsgründen beginnen wir dabei mit einer simulationsbasierten Implementation, ganz ähnlich wie beim Random Walk. Im Anschluss führen wir dann das R-Paket dRiftDM ein, mit welchem wir im Folgenden dann weiter arbeiten werden.

4.4.1 Manuelle Implementation per Simulation

Um das Diffusionsmodell zu simulieren, müssen wir vom eigentlich kontinuierlichen Modell wieder Diskretisierungen der Zeit (und auch des States) einführen. Mit der stochastischen Euler Methode wird ausgehend von einem initialen Status \(X_0=\mathbf{X}(0)\) (bei uns zunächst mal auf \(0\) gesetzt) der Wert von \(\mathbf{X}\) “vorwärts” entwickelt, wobei wir dazu die stochastische Differenzgleichung von oben benutzen können.

Programmiertechnisch ergibt sich der (approximierte) Wiener-Prozess dann als Summe standard-normalverteilter und skalierter Zufallsvariablen. In anderen Worten: Im Wesentlichen ziehen wir sehr häufig aus einer Standardnormalverteilung und berechnen dann die kumulierte Summe über die betrachteten Zeitschritte. Statt einer for-Schleife verwenden wir hier hier eine Matrixschreibweise, mit der alle n_trials-vielen Diffusionsprozesse gleichzeitig erzeugt werden.

Der Einfachheit halber betrachten wir den Fall, dass \(\Delta t = 1\) ist. Dies ist allerdings in unseren Kontexten durchaus auch üblich und dieser Zeitschritt würde konzeptuell dann einer Veränderung in der Zeit von 1 ms entsprechen. Die Formel vereinfacht sich dadurch zu \[ \mathbf{X}(t+\Delta t)=\mathbf{X}(t)+\mu+\sigma\cdot{\bf Z}\,. \]

Auch der folgende Code nimmt dies an (für andere Zeitschritte müssten einige Modifikationen eingebaut werden):

# Funktion zur Simulation vieler Diffusionsprozess
simulate_DM <- function(t_max = 200, # Wieviele Zeitschritte?
                        n_trials = 10, # Wieviele Trials / Diff-Prozesse?
                        mu = 1, # Drift Rate
                        boundary = 10, # Boundary
                        t0 = 300, # t0
                        st0 = 30, # "Standardabweichung" t0
                        sigma = 1) { # sigma

  # Matrix aus der Drift Rate erstellen: systematischer Anteil
  D <- matrix(
    data = mu,
    nrow = n_trials,
    ncol = t_max,
    byrow = TRUE
  )
  # Matrix aus standard-normalverteilter Zufallsvariable erstellen: Wiener-Prozess
  W <- matrix(
    rnorm(
      n = n_trials * t_max,
      mean = 0,
      sd = 1
    ) * sigma,
    nrow = n_trials,
    ncol = t_max,
    byrow = TRUE
  )
  # für beide Matrizen zeilenweise die kumulative Summe bilden
  D <- t(apply(
    X = D,
    MARGIN = 1,
    FUN = cumsum
  ))
  W <- t(apply(
    X = W,
    MARGIN = 1,
    FUN = cumsum
  ))
  # Evidenz X = systematischer Anteil + Wiener-Prozess
  X <- D + W

  # zeilenweise die First-Passage Times bestimmen
  firstPassageTimes <- apply(
    X = X,
    MARGIN = 1,
    FUN = function(X) which(abs(X) > boundary)[1]
  )
  # zeilenweise bestimmen, welche Boundary über-/unterschritten wurde
  response <- apply(
    X = X,
    MARGIN = 1,
    FUN = function(X) sign(X[which(abs(X) > boundary)[1]])
  )

  # zur Bestimmung der RTs noch eine gleichverteilte Zufallsvariable zu den
  # First-Passage Times addieren
  RTs <- firstPassageTimes + runif(
    n = n_trials,
    min = t0 - st0,
    max = t0 + st0
  )

  # dann werden Reaktionen, First-Passage Times und RTs zu einem DataFrame
  # zusammengefasst, wobei noch alle NAs entfernt werden, d.h. die Instanzen,
  # wo der Diffusionsprozess bis zum Zeitpunkt t_max keine boundary erreicht hat
  DM <- data.frame(
    response,
    firstPassageTimes,
    RTs
  )
  DM <- na.omit(DM)

  # Schließlich wird der DataFrame an eine aufrufende Funktion zurückgegeben
  return(DM)
}

Nun können wir eine Vielzahl an Diffusionsprozessen gleichzeitig simulieren und die Eigenschaften des Modells explorieren. Zunächst betrachten wir nur Mittelwerte und die Anteile der Reaktionen. Da wir \(\mu>0\) wählen, betrachten wir “untere Reaktionen” als fehlerhafte Reaktionen:

set.seed(3)
myDM <- simulate_DM(
  t_max = 2000, # Wieviele Zeitschritte?
  n_trials = 5000, # Wieviele Trials / Diff-Prozesse?
  mu = 0.3, # Drift Rate
  boundary = 70, # Boundary
  t0 = 300, # t0
  st0 = 30, # "Standardabweichung" t0
  sigma = 1
)

Übung: Schreiben Sie die Funktion analyze_sim() aus Kapitel 3 so um, dass nun die RTs (statt die First-Passage Times) analyiert werden. Dann können Sie bereits das Verhalten des Modells explorieren.

4.4.2 Diffusionsmodell mit dRiftDM

Das R-Paket dRiftDM wurde von uns zur flexiblen Arbeit mit komplexeren Diffusionsmodellen entwickelt (siehe Richter et al., 2023, für die mathematischen Grundlagen). Wir beginnen hier natürlich mit der einfachsten Variante, die wir in diesem Kapitel kennengelernt haben und befassen uns erst im nächsten Kapitel mit einer komplexeren Variante.

Das Modell mit den Parametern \(\mu\), \(b\) und \(t0\) (als Konstante) wird in dRiftDM als ratcliff_dm() bezeichnet und kann sehr leicht erstellt werden:

library(dRiftDM) # Paket laden

## 
##  _________________________________________________________ 
## / Welcome to dRiftDM 0.2.2.9000 Please report any bugs or \
## \ unexpected behavior                                     /
##  --------------------------------------------------------- 
##       \
##        \
## 
##         ^__^ 
##         (oo)\ ________ 
##         (__)\         )\ /\ 
##              ||------w|
##              ||      ||

# Installation der dRiftDM Version, mit der dieses Skript geschrieben wurde: 
# remotes::install_github("bucky2177/dRiftDM", ref = "1d6de63237f688aa46aec47cbf6cc7080c28f871")

myDDM <- ratcliff_dm() # Modell erstellen
print(myDDM) # Modell anzeigen

## Class(es): ratcliff_dm, drift_dm
## 
## Current Parameter Matrix:
##      muc   b non_dec
## null   3 0.6     0.3
## 
## Unique Parameters:
##      muc b non_dec
## null   1 2       3
## 
## Deriving PDFs:
##   solver: kfe
##   values: sigma=1, t_max=3, dt=0.001, dx=0.001, nt=3000, nx=2000
## 
## Observed Data: NULL

Unter Current Parameter Matrix sehen wir, dass es derzeit drei Parameter gibt:

• muc meint die Drift Rate \(\mu\)
• b meint die Boundary, bei welcher eine Entscheidung getroffen wird
• non_dec ist die Residualzeit \(t0\), hier ersteinmal als Konstante

Da für alle drei Parameter unter Unique Parameters eine verschiedene Zahl vergeben ist, sind alle drei Parameter frei veränderbar und würden auch geschätzt werden können. In dRiftDM muss jedes Modell mindestens aus einer Bedingung bestehen. Da dieses einfache Modell allerdings keine zwei Bedingungen unterscheidet, wird die eine Bedingung hier als Platzhalter mit null bezeichnet. Im nächsten Kapitel lernen wir dann auch ein Modell mit mehreren Bedingungen kennen. Unter Deriving PDFs finden wir noch einige weitere Informationen, u.a. begegnen uns hier wieder sigma, dt, dx und t_max; solver: kfe meint, dass die PDFs der (korrekten und fehlerhaften) RTs durch das Lösen der sog. Kolmogorov-Forward Equation vorhergesagt werden, statt durch Simulationen (siehe Richter et al., 2023). Bevor wir uns mit diesem Modell nun im Detail befassen, weisen wir noch (einmal) auf drei Aspekte hin:

1. Standardmäßig verwendet dRiftDM ein accuracy coding, d.h., bei positiver Drift Rate \(\mu\) werden “obere” Reaktionen als korrekt und “untere” Reaktionen als Fehler gewertet.
2. Die Boundaries spannen den symmetrischen Raum \([-b;b]\) auf und der Startpunkt liegt daher bei 0. Andere Softwarelösungen spannen den Raum \([0,a]\) auf und der Startpunkt wäre dann bei \(a/2\). Diese beiden Varianten betreffen aber nur die Notation und verändern nicht das Modellverhalten an sich.
3. Bisher haben wir Zeit in Millisekunden betrachtet, was oft intuitiver ist, wenn Random Walks eingeführt werden. Bei Diffusionsmodellen ist es aber nicht unüblich in der Einheit Sekunden zu rechnen und dRiftDM arbeitet daher auch entsprechend so. Allerdings ist dies auch nur eine Frage der Skalierung und die Parameter wären alle ineinander umrechenbar, wenn eine andere Maßeinheit verwendet werden soll (siehe z.B. hier).

Als erstes können wir uns sehr leicht den erwarteten Verlauf der Evidenzakkumulation sowie auch einzelne Diffusionsprozesse visualisieren. Dazu nutzen wir die Funktion simulate_traces() in der wir mit sigma = 0 die Brown’sche Bewegung quasi unterdrücken:

# Erwarteter Verlauf der Evidenzakkumulation
expected <- simulate_traces(object = myDDM, k = 1, sigma = 0)
# k = 5 Diffusionsprozesse
trajectories <- simulate_traces(object = myDDM, k = 5, sigma = 1)

# Die Ergebnisse können mit der generischen Funktion plot() visualisiert werden
par(mfrow = c(1,2))
plot(expected)
plot(trajectories)

Wir können das Modellverhalten auch weiter explorieren, indem wir einzelne Parameterwerte ändern und schauen, wie sich die Modellvorhersagen dadurch ändern:

coef(myDDM) # Parameter ausgeben

##     muc       b non_dec 
##     3.0     0.6     0.3

calc_stats(myDDM, type = "basic_stats") # Grundlegende Statistiken ausgeben

# Nun verändern wir die Drift Rate und schauen uns die Statistiken erneut an
coef(myDDM)["muc"] <- 2.5
coef(myDDM) # Parameter ausgeben

##     muc       b non_dec 
##     2.5     0.6     0.3

calc_stats(myDDM, type = "basic_stats") # Grundlegende Statistiken ausgeben

Erwartungsgemäß werden bei positiver Drift Rate, also \(\mu>0\) deutlich mehr korrekte Reaktionen (P_corr) als Fehler erzeugt. Wir sehen aber auch, ganz ähnlich wie bei Random Walks, dass die mittleren RTs korrekter und fehlerhafter Reaktionen identisch sind (i.Ü. sind auch die Standardabweichungen identisch und, bis auf die Skalierung, sind sogar die gesamten PDFs identisch) und wenn die Drift Rate kleiner wird, werden die RTs länger und die Fehler nehmen leicht zu. Dass korrekte und fehlerhafte RTs im Mittel gleich lang sind, ist empirisch allerdings eher selten der Fall; i.d.R. sind Fehler entweder langsamer oder schneller (durch Raten). Wir werden darauf später zurückkommen. Zunächst explorieren wir das Modell und sein Verhalten in Abhängigkeit der drei Parameter \(\mu\), \(b\) und \(t0\) etwas systematischer. Im weiteren Verlauf werden wir weitere Anwendungen von dRiftDM kennenlernen, wenngleich die eigentlichen Stärken des Pakets diese Einführung bei Weitem überschreiten (in unserem Online-Buch finden sich dazu viele weitere Informationen).

4.5 Modellverhalten und Validierung

Ein gutes Verständnis des Modellverhaltens ist sowohl für die Anwendung als auch vor allem für die korrekte Interpretation von Modellparametern enorm wichtig. Wir schauen uns daher zunächst das Verhalten des bisher eingeführten Diffusionsmodells an und im Anschluss, was die Bedeutung der einzelnen Parameter ist und wie diese experimentell validiert werden können.

4.5.1 Modellverhalten

Um das Modellverhalten systematischer zu verstehen, kann man einen Parameter variieren und die anderen zwei konstant halten und dann die Auswirkungen auf die RTs und die relative Häufigkeit korrekter Reaktionen (und damit die Fehlerraten) betrachten. Abbildung 4.6 illustriert die Effekte der einzelnen Parameter auf die vorhergesagten mittleren korrekten RTs (schwarze Linie) und die Accuracy (blaue Linie), also die relative Häufigkeit korrekter Reaktionen, wenn systematisch ein Parameter variiert wird. Als Defaultwerte für die anderen Parameter wurde für die Abbildungen angenommen, dass

• \(\mu = 3\),
• \(b = 0.6\) und
• \(t0 = 0.3\)

ist. Dann sind für den linken, mittleren und rechten Teil der Abbildung folgende Parameterbereiche variiert worden, während die jeweils anderen zwei Parameter auf den Defaultwert fixiert waren:

• Drift rate: \(\mu \in \{1.5, 2.0, \ldots, 4.5, 5.0\}\)
• Boundary: \(b \in \{0.2, 0.3, \ldots, 0.8, 0.9\}\)
• Residualzeit: \(t0 \in \{0.1, 0.15, \ldots, 0.40, 0.45\}\)

Abbildung 4.6: Illustration der Wirkungen von \(\mu\) (links), \(b\) (Mitte) und \(t0\) (rechts) auf die Modellvorhersagen (RTs; linke y-Achse) relative Häufigkeit korrekter Reaktionen (proportion correct = PC; rechte y-Achse).

Offenkundig wirkt sich zunächst jeder Parameter auf die RTs aus, aber in jedem Fall sind die Wirkungen auf die Fehlerrate verschieden (eine gut verständliche Einführung in diese (und zusätzliche) Parameter sowie ein Programm zur interaktiven Exploration des Diffusionsmodells findet sich auch in Alexandrowicz, 2020):

• Wenn \(\mu\) größer wird, werden (1) die RTs kürzer, aber (2) auch die Fehlerrate nimmt ab (die Accuracy wird höher).
• Wenn \(b\) zunimmt, werden (1) die RTs länger und (2) die Fehlerrate nimmt ab (die Accuracy wird höher). Dies ist die Situation, die als Speed-Accuracy Tradeoff (SAT) bezeichnet wird (vgl. Liesefeld & Janczyk, 2019, 2023). Hiermit ist das grundsätzliche Problem gemeint, dass Menschen nicht gleichzeitig schnell und korrekt reagieren können, sondern immer eine Balance zwischen beiden Polen finden müssen. Wird ein höherer Wert für \(b\) gewählt, wird dies auch als “vorsichtiger Reaktionstyp” interpretiert.
• Wenn \(t0\) größer wird, werden (1) die RTs länger, während (2) die Fehlerrate gleich bleibt.

4.5.2 Interpretation und Validierung der Parameter

Manche Interpretationen haben wir weiter oben bereits angedeutet. Da wir nun auch wissen, wie das Modell sich in Bezug auf die vorhergesagten RTs und Fehlerraten in Abhängigkeit der Parameter verhält, stellen wir hier die wichtigsten Interpretationen noch einmal vor.

Die Drift Rate wird klassischerweise mit Stimulusqualität in Verbindung gebracht (z.B. Ratcliff & McKoon, 2008). Dies äußert sich auf zweierlei Art: Werden z.B. separate Drift Rates für die beiden Stimuli geschätzt, die jeweils verschiedene Reaktionen erfordern, dann wird einmal \(\mu>0\) und einmal \(\mu<0\) resultieren (dann müssen auch \(-b\) und \(b\) jeweils für die beiden Reaktionen statt für Fehler und korrekte Reaktionen stehen). Allerdings spielt für den absoluten Wert der Drift Rate auch die Qualität des Stimulus insofern eine Rolle, als dass z.B. gut diskriminierbare Stimuli mit einer höheren Drift Rate einhergehen als schlecht diskriminierbare Stimuli. Ein weiterer Aspekt ist, dass die Drift Rate oft mit Geschwindigkeit der Informationsverarbeitung und Intelligenz in Zusammenhang gebracht wird (z.B. Lerche et al., 2020; Schmiedek et al., 2007; Schubert et al., 2016). Insgesamt wird die Drift Rate oft mit der Effizienz der zielgerichteten, kontrollierten Übersetzung eines Stimulus in eine Reaktion, also der Reaktionswahl, in Verbindung gebracht.

Es gibt aber noch eine andere, eher mathematische Interpretation, der Drift Rate, die wir bereits angesprochen haben. Bei einer Standard “2-Alternative Forced-Choice Aufgabe” mag es bspw. sinnvoll sein von einem linearen Anstieg der Evidenz mit der Zeit auszugehen. Diese Funktion würde dann beschrieben werden durch \[\mathbb{E}[\mathbf{X}(t)] = \mu\cdot t\,.\] Die Drift Rate ist dann mathematisch einfach die erste Ableitung dieser Funktion nach \(t\), also \[\frac{d\mathbb{E}[\mathbf{X}(t)]}{dt}=\mu\,,\] und sie gibt damit auch die Steigung der erwarteten Funktion für jeden Zeitpunkt \(t\) an (siehe Cox & Miller, 1965; Schwarz, 2022). Da \(\mu\) eben eine Konstante ist und selber nicht von \(t\) abhängt, hatten wir diese Variante als zeit-unabhängig bezeichnet.

Die Boundary hatten wir oben schon mit SATs in Verbindung gebracht. Höhere Werte stehen eher für einen konservativeren und vorsichtigeren Reaktionsstil, bei dem mehr Wert auf korrekte Reaktionen gelegt wird, während die längeren RTs dabei in Kauf genommen werden. Veränderungen dieses Modellparameters werden manchmal auch als strategische Anpassung des Reaktionskriteriums interpretiert. Werden in einem Block viele Catch Trials untergebracht, also welche, die keine Reaktion erfordern, dann ist i.d.R. die Boundary auch höher als bei wenig/keinen Catch Trials. Dies ist ein Hinweis darauf, dass Personen das Kriterium anpassen, ab dem sie eine Reaktion auswählen, wenn eine vorschnelle Reaktion der Aufgabenstellung eher nicht zuträglich ist.

Die Residualzeit schließlich umfasst alle Prozesse, die nicht mit der Reaktionsauswahl im engeren Sinne zu tun haben. Klassischerweise sind dies sehr frühe perzeptuelle Prozesse, aber auch Prozesse die mit der Programmierung der motorischen Efferenz und der Innervierung des motorischen Systems einhergehen. Allerdings sind auch eine ganze Reihe anderer Prozesse mit der Residualzeit in Verbindung gebracht worden. Werden z.B. in einem Experiment zwei verschiedene Aufgaben bearbeitet, kommt es zu längeren RTs, wenn vom Trial vorher zum aktuellen Trial die Aufgabe gewechselt werden muss, verglichen mit, wenn sie wiederholt wird. Diese sog. Aufgabenwechselkosten schlagen sich (auch) in der Residualzeit nieder (Schmitz & Voss, 2012).

Zunächst einmal sind solche Interpretationen oft eher theoretisch und sie werden “angenommen”. Validierungsstudien testen dann, ob die Parameter sich so verhalten, wenn experimentelle Manipulationen vorgenommen werden, die–idealerweise selektiv–einen bestimmten (aber keine anderen) Parameter betreffen sollten. Eine solche Studie im Bereich Diffusionsmodell stammt von Voss et al. (2004). In den Experimenten dieser Studie wurde den Versuchspersonen ein rechteckiger Stimulus präsentiert, bei dem zufällig ein bestimmter Anteil an Pixeln orange und der andere Anteil blau gefärbt wurde. Die (manuelle) Reaktion sollte auf den jeweils dominanten Farbanteil abgegeben werden. In der Baseline-Bedingung war der Anteil 53% vs. 47%.

In weiteren Blöcken wurden dann folgende Variationen durchgeführt (Experiment 1):

• Aufgabenschwierigkeit: Erhöhung der Aufgabenschwierigkeit durch Veränderung der Farbanteile zu nur noch 51.5% vs. 48.5%. Dadurch sollte die Diskrimination nun schwerer werden und die Drift Rate sollte kleiner ausfallen.
• Accuracy-Instruktion: Induktion eines konservativeren, vorsichtigeren Reaktionsstils durch Verlängerung der Reaktions-Deadline (maximale Zeit bis zur Reaktion) und Instruktion, die das Vermeiden von Fehlern betont. Dadurch sollte die Boundary einen höheren Wert annehmen.
• Reaktion komplexer: Erhöhung der motorischen Schwierigkeit dadurch, dass die Reaktionsfinger auf einem Home-Button liegen müssen und von dort aus zu den Reaktionstasten gehen.

Die folgende Tabelle ist eine konzeptuelle Zusammenfassung der Ergebnisse von Voss et al. (2004, Tab. 2), die die Parameterveränderungen immer in Bezug zur Baseline-Bedingung zeigt. Die betroffenen Parameter werden also tatsächlich durch die experimentellen Manipulationen wie vorhergesagt verändert. Allerdings scheint sich z.B. die Accuracy-Instruktion auch auf die Residualzeit auszuwirken. (Das vollständige Befundmuster ist etwas komplexer, siehe dafür die Originalarbeit.)

4.6 Das vollständige Diffusionsmodell von Ratcliff (1978) und Ratcliff & Rouder (1998)

Wir haben bisher die “Hauptparameter” eines Diffusionsmodells betrachtet, in etwa so, wie es von Ratcliff (1978) im Rahmen einer Studie zu Gedächtnisphänomenen (in die Psychologie) eingeführt wurde. Mit diesen Parametern kann das Modell tatsächlich schon eine ganze Reihe empirischer Phänomene reproduzieren; es gibt aber auch Grenzen. Zum einen kann es sein, dass eine Versuchsperson einen Bias für eine bestimmte Reaktion aufweist. Dies kann dadurch zustande kommen, dass diese Reaktion besonders häufig erfordert ist oder auch, indem diese Reaktion bei einer (korrekten) schnellen Antwort zusätzlich belohnt wird. Im bisherigen Modell hat die Evidenzakkumulation immer bei \(\mathbf{X}(0)=0\) angefangen. Wir können den Startwert aber auch allgemeiner als Bias-Parameter \(z\) bezeichnen (mit \(z\in [-b,+b]\)). Mit \(z=0\) liegt (wie bisher) kein Bias vor, \(z>0\) bzw. \(z<0\) wären dann jeweils ein Bias für die obere bzw. die untere Reaktion.

Zum anderem haben wir oben festgestellt, dass die mittleren RTs für korrekte und fehlerhafte Reaktionen gleich waren. Fassen wir die unteren Reaktionen als “falsch” auf, ist dieses Verhalten aber komisch, da empirische RTs fehlerhafter Reaktionen oft schneller (z.B. durch geratene Reaktionen unter Zeitdruck) oder langsamer als die der richtigen Reaktionen sind. Derartige Beobachtungen kann das Diffusionsmodell vorhersagen, wenn angenommen wird, dass der Startpunkt \(z\) und die Drift Rate \(\mu\) eines Diffusionsprozesses (bei uns bisher immer bei \(\mathbf{X}(0)=z=0\)) über die Trials hinweg nicht konstant, sondern variabel sind (Ratcliff & Rouder, 1998). Eine gewisse Variabilität dieser beiden Parameter erscheint auch psychologisch plausibel, da z.B. die Evidenzaufnahme auch von Fluktuationen der Aufmerksamkeit o.ä. abhängen kann.

Eine Erweiterung des Modells ist es also, Parameter für die Variabilitäten der Parameter \(\mu\) und \(z\) einzuführen, die wir hier kurz \(S_{\mu}\) und \(S_z\) nennen. Während oft angenommen wird, dass \(\mu\) normalverteilt und \(z\) (sowie ja auch \(t0\)) gleichverteilt seien, kommt Ratcliff (2013) zu dem Schluss, dass–mit wenigen Ausnahmen–die Parameterschätzungen eher unabhängig von den tatsächlichen Verteilungen der Parameter sind Lassen wir die Drift Rate variieren, führt dies im Modell zu langsamen Fehlern; lassen wir hingegen den Startpunkt variieren, führt dies zu schnellen Fehlern.

Insgesamt haben wir nun ein Modell mit sieben Parametern:

• Drift Rate \(\mu\): Qualität des Stimulus und Effektivität der Informationsaufnahme
• Boundary \(b\): Reaktions- bzw. Auswahlkriterium (liberal vs. konservativ)
• Residualzeit \(t0\): Prozesse außerhalb der Entscheidung/Auswahl (Wahrnehmung, Motorik)
• Startwert \(z\): ermöglicht Bias in Richtung einer Reaktion
• \(S_z\): Variabilität des Startpunktes \(z\) ermöglicht schnelle Fehler
• \(S_{\mu}\): Variabilität der Drift Rate \(\mu\) ermöglicht langsame Fehler
• \(S_{t0}\): Variabilität der Residualzeit \(t0\)

In vielen Anwendungen werden allerdings nicht alle Parameter benötigt bzw. geschätzt und die wichtigsten Parameter dürften tatsächlich \(\mu\), \(b\) und \(t0\) (ggf. mit \(S_{t0}\)) sein. Auch zeigen manche Arbeiten, dass Diffusionsmodelle mit weniger Parametern durch ihre Sparsamkeit häufig komplexeren Modellen überlegen sind (Lerche & Voss, 2016). Weitere Einführungen und mehr Informationen finden sich auch in Voss et al. (2013), Smith und Ratcliff (2015) oder Ratcliff et al. (2016).

Wollen wir in dRiftDM z.B. Variabilität der Residualzeit (in Form einer Gleichverteilung) und des Startwerts hinzufügen und damit schnelle Fehler ermöglichen, ist dies relativ einfach:

myDDM_2 <- ratcliff_dm(var_non_dec = TRUE, var_start = TRUE)
print(myDDM_2)

## Class(es): ratcliff_dm, drift_dm
## 
## Current Parameter Matrix:
##      muc   b non_dec range_non_dec range_start
## null   3 0.6     0.3          0.05         0.5
## 
## Unique Parameters:
##      muc b non_dec range_non_dec range_start
## null   1 2       3             4           5
## 
## Deriving PDFs:
##   solver: kfe
##   values: sigma=1, t_max=3, dt=0.001, dx=0.001, nt=3000, nx=2000
## 
## Observed Data: NULL

calc_stats(myDDM_2, type = "basic_stats") # Grundlegende Statistiken ausgeben

Wir sehen nun, dass die mittlere RTs der fehlerhaften Reaktionen kürzer ist als die der korrekten Reaktionen.

4.7 Anwendungen des Modells

Wir haben jetzt zugegebenermaßen viel Vorarbeit für ein grundlegendes Verständnis von Diffusionsmodellen geleistet (und dabei dennoch nur einen Teil der möglichen Materie behandelt). Klassisch wurden (und werden) viele Experimente auf Basis mittlerer RTs (und Fehlerraten) ausgewertet. Dies mag in manchen Fällen durchaus OK sein, hat aber einen entscheidenden Nachteil: (Mittlere) RTs geben immer nur einen groben Einblick in das Prozessgeschehen bei der Bearbeitung einer Aufgabe. Tatsächlich werden Diffusionsmodelle, die ja auf der gesamten Verteilung von RTs beruhen, in vielen Feldern der psychologischen Forschung verwendet.

Es mag z.B. sein, dass ältere Versuchspersonen langsamer reagieren als jüngere Versuchspersonen. Damit wissen wir aber noch nicht, warum dies so ist. Die Anwendung eines Diffusionsmodells, d.h. die Parameterschätzung für jede Versuchsperson, zerlegt die RT in ihre Bestandteile und die Parameter können anschließend als abhängige Variablen behandelt und ausgewertet werden. Dies liefert uns Informationen darüber, in welchem Parameter sich z.B. Gruppen von Versuchspersonen unterscheiden und, insofern es psychologische Interpretationen der Parameter gibt, damit auch über die Gründe des eher globalen RT Unterschieds.

Ein anderer Aspekt ist, dass wir gesehen haben, wie mehrere Parameter und damit auch Prozesse Einfluss auf die RT haben. Da es naheliegend ist, dass Versuchspersonen auf allen Parametern mehr oder weniger viel Variation zeigen, steigt die Sensitivität für die Entdeckung eines Unterschieds, wenn dieser in einem Modellparameter ausgedrückt wird, dadurch, dass die Wirkung der anderen Parameter quasi “herausberechnet” wird. Modellparameter wie die Drift Rate können also z.B. ein reineres Maß für Unterschiede in der Bearbeitungsgeschwindigkeit sein, die bei Verwendung von RTs als solche mitunter nicht entdeckt werden.

4.7.1 Praktische Bestimmung der Modellparameter

Ein Punkt, den wir in diesem Kapitel übergangen haben, ist die praktische Bestimmung der Modellparameter für vorliegende empirische Daten. Das prinzipielle Vorgehen entspricht allerdings exakt dem, was wir in Kapitel 3 am Beispiel des Random Walk Modells bereits ausgeführt haben. Allerdings werden in der Regel solche Parameterbestimmungen nicht individuell programmiert, da es eine Reihe gut etablierter Sofwarelösungen gibt, die genau dies effizient tun. Für einfache Modelle gibt es sogar geschlossene Lösungen, auf denen das EZ-Diffusion Model basiert (Groulx et al., 2020; Wagenmakers et al., 2007). Ein weiteres Programm ist fast-dm (Voss & Voss, 2007; Voss et al., 2015). Dieses Programm ist eigenständig und bietet eine effiziente Bestimmung aller Parameter zeit-unabhängiger Diffusionsmodelle. Dabei ist auch eine gewisse Flexibilität gegeben, Parameter zwischen Bedingungen zu fixieren oder frei zu schätzen.

Beide Varianten haben den Nachteil, dass sie nur mit Situationen umgehen können, in denen die Parameter innerhalb eines Prozess konstant, also zeitunabhängig, sind. Im nächsten Kapitel lernen wir auch ein Modell kennen, in welchem die Drift Rate sich zeitabhängig ändert. Zur Bestimmung der Parameter solcher generellen Diffusionsmodelle (die bisher verwendete Ratcliff-Variante ist ja eher ein Spezialfall) kann wiederum das R-Paket dRiftDM (Koob et al., 2024) genutzt werden. Anhand echter Datensätze aus Konfliktaufgaben (Eriksen Flanker-Aufgabe und Simon-Aufgabe) werden wir uns im nächsten Kapitel mit der Parameterschätzung befassen.

4.7.2 Aging-Forschung

Aus vielen Studien ist bekannt, dass ältere Versuchspersonen i.d.R. längere RTs bei der Bearbeitung einer Vielzahl von Aufgaben zeigen, als jüngere Versuchspersonen. Diese Beobachtung wird oft mit einem generalized slowing erklärt, also der Idee, kognitive Prozesse würden mit zunehmendem Alter langsamer, was eben zu längeren RTs führt. Während diese Standarderklärung aber auf (mittlere) RTs fokussiert, erlaubt ein Diffusionsmodell die detailliertere Exploration der Ursachen für die RT-Verlängerung. Eine Verschlechterung in der Qualität der Informationsverarbeitung würde naheliegend eine verringerte Drift Rate vorhersagen.

Ratcliff et al. (2006) haben daher Versuchspersonen aus drei verschiedenen Altersgruppen (Studierende vs. 60-74 Jahre vs. 75-85 Jahre; jeweils \(n=10\)) vier verschiedene Aufgaben bearbeiten lassen:

• Signal Detection: Es wurde eine Menge von 1-100 Asterisken erzeugt, wobei die Menge entweder aus einer Normalverteilung mit Erwartungswert 57.5 (Signal-Verteilung) oder 39.5 (Noise-Verteilung) stammte. Die Asterisken wurden dann an zufälligen Positionen in einem \(10\times 10\)-Gitter präsentiert und die Versuchspersonen sollten angeben, ob die Menge der Asterisken groß oder klein ist. Bei sehr großen oder sehr kleinen Mengen wurde das entsprechende Feedback (“groß” vs. “klein”) als richtige Antwort präsentiert, bei mittelgroßen Mengen wurde das Feedback zufällig bestimmt.
• Letter Discrimination: Hier wurde den Versuchspersonen ein Buchstabe präsentiert und es sollte mit einer Tastenreaktion die Identität des Buchstabens angegeben werden. Der Buchstabe war dabei jeweils nur kurz zu sehen und er wurde nach 10, 20, 30, 40, 50 oder 60 ms maskiert (was die weitere Wahrnehmung unmöglich macht).
• Brightness Discrimination: In einem Gitter von \(64\times 64\) Pixeln wurde die Helligkeit des Stimulus dadurch manipuliert, dass die Wahrscheinlichkeit dafür, dass ein Pixel weiß war von \(p=0.35\) bis \(p=0.65\) in sechs Stufen manipuliert wurde. Die Versuchspersonen sollten angeben, ob der Stimulus hell oder dunkel ist.
• Recognition Memory: Verwendet wurden Wörter mit hoher, geringer und sehr geringer sprachlicher Vorkommensfrequenz. In jedem Lernblock wurden 9 Wörter einmal und 9 Wörter dreimal präsentiert (jeweils 3 Wörter jeder Frequenz) und im Anschluss wurden diese 18 und 18 neue Wörter (jeweils 6 jeder Frequenz) in einer Recognition-Aufgabe präsentiert.

Zusätzlich wurden alle Aufgaben unter “Speed-Instruktionen” und unter “Accuracy-Instruktionen” bearbeitet. Diese Manipulation zeigte auf behavioraler Ebene die beabsichtige Wirkung und generell wurden die RTs länger, je älter die Versuchspersonen waren (während die Fehlerraten eher wenig Unterschiede zeigten). Im Wesentlichen repliziert dies also die typischen Ergebnisse, die als generalized slowing interpretiert werden können.

Die Autor:innen haben nun die Parameter eines Diffusionsmodells für jede Versuchsperson und für jede Aufgabe geschätzt und analysiert. Der wesentliche Befund dieser Studien ist, dass ein Großteil der RT-Unterschiede auf „conservative response criteria, accompanied by relatively small slowdowns in nondecision components of processing“ (S. 634) zurückgeht. In anderen Worten: Ältere Versuchspersonen legen möglicherweise ein vorsichtigeres Reaktionskriterium an; Einbußen in der Qualität der Verarbeitung müssen den längeren RTs aber nicht notwendigerweise zugrunde liegen.

4.7.3 Aufmerksamkeitsbias durch Ängstlichkeit

Es gibt viele Befunde, dass Personen mit höherer Ängstlichkeit einen Threat Bias zeigen, eine präferentielle Verarbeitung bedrohlicher Stimuli im Vergleich zu Personen mit geringer Ängstlichkeit. Dies trifft sowohl auf klinische Stichproben zu, aber auch innerhalb von nicht-klinischen Stichproben auf Personen mit hoher Trait-Ängstlichkeit (siehe Bar-Haim et al., 2007, für eine Metaanalyse). In diesem Zuge wurde vorgeschlagen, dass der Threat Bias nicht auf eine höhere Effizienz der Verarbeitung bedrohlicher Stimuli zurückgeht, sondern auf eine präferierte Zuteilung von Aufmerksamkeit in Situationen mit mehreren Stimuli.

Evidenz dafür stammt aus nicht wenigen Studien mit verschiedenen Aufgaben. In Probe Detection Aufgaben werden z.B. zunächst ein bedrohliches und ein neutrales Wort (oder auch Gesicht) an verschiedenen Lokationen kurz präsentiert, gefolgt von einem Probe-Stimulus an einer der Positionen. Die Versuchspersonen sollen dann die Anwesenheit oder auch Art/Lokation des Probe-Stimulus berichten. Hoch-ängstliche Versuchspersonen zeigen dabei kürzere RTs, wenn der Probe-Stimulus das bedrohliche Wort/Gesicht ersetzt hat. Während solche Studien generell einen Threat Bias bei Anwesenheit mehrerer Stimuli zeigen, gibt es auch direktere Evidenz für die oben beschriebene Hypothese. MacLeod and Mathews (1991) haben eine lexikalische Entscheidungsaufgabe bei einer klinischen Stichprobe mit generalisierter Angststörung und einer Kontrollgruppe verwendet. In einer solchen Aufgabe muss entschieden werden, ob eine Folge an Buchstaben ein Wort ergibt oder nicht. Wichtig ist nun die Unterscheidung zwischen zwei Bedingungen: In der ersten Bedingung wurde jeweils nur eine Buchstabenfolge präsentiert. In der zweiten Bedingung wurden zwei Buchstabenfolgen präsentiert und die “Wort”-Reaktion sollte dann gegeben werden, wenn mindestens eine der Buchstabenfolgen ein Wort ergibt. Im Ergebnis zeigten sich kürzere RTs in der klinischen Stichprobe verglichen mit der Kontrollgruppe nur in der zweiten, nicht aber in der ersten Bedingung.

White et al. (2010a) diskutieren eine Vielzahl derartiger Studien inkl. welcher, die eher gegen die Hypothese von oben sprechen, und argumentieren, (1) dass ein Problem derartiger Studien immer ist, dass das kritische Argument auf der Beibehaltung einer Nullhypothese basiert, also kein Threat Bias, wenn es nur einen Stimulus gibt. (2) Darüber hinaus sollten Threat Biases auch in solchen Situationen auftreten, aber die Effektstärke könnte dabei einfach kleiner sein und daher wären diese Biases mit traditionellen Methoden wie Unterschiede in mittleren RTs nicht entdeckbar. Dies wird durch eher kleine Stichproben und eine geringe Anzahl kritischer Stimuli noch verstärkt, beides Eigenschaften die zu geringer Power der statistischen Tests führen. In ihrer eigenen Studie fokussieren White et al. dann auf die Verwendung sensitiverer abhängiger Variablen. Dabei machen sie sich die Eigenschaft zunutze, dass–wenn RTs und Fehlerraten mit einem Diffusionsmodell analysiert werden–die verschiedenen Einflüsse auf die RTs in den verschiedenen Parametern abgebildet werden. Damit wird die Drift Rate \(\mu\) zu einem reineren Maß für die Verarbeitungseffizienz als es (mittlere) RTs sind, da andere Einflüsse auf die RT, wie z.B. Veränderungen der Boundary oder der Residualzeit, quasi herausgerechnet werden (siehe dazu auch White et al., 2010b).

In Experiment 1 in der Studie von White et al. (2010a) wurden jeweils zwei Gruppen mit hoher und niedriger Trait-Ängstlichkeit miteinander verglichen (bestimmt mit dem Spielberger Trait Anxiety Inventory [STAI]; Spielberger et al., 1970). Die Versuchspersonen bearbeiteten dann eine lexikalische Entscheidungsaufgabe, bei der immer nur eine Buchstabenfolge zur gleichen Zeit präsentiert wurde. Die dabei verwendeten Wörter bestanden zu einem überwiegenden Anteil aus “Filler Items” und eben den kritischen bedrohlichen Wörtern und den Kontrollwörtern. Für mittlere RTs und Fehlerraten gab es einen leichten deskriptiven Trend in Richtung eines größeren Threat Bias für die Versuchspersonen mit hoher Trait-Ängstlichkeit, allerdings war die (statistisch kritische) Interaktion jeweils nicht signifikant. Diese Befunde replizieren also die oben berichtete Situation. Mit der Drift Rate \(\mu\) als abhängiger Variable war die Situation allerdings anders: Hier zeigte sich ein größerer Threat Bias in der Gruppe mit hoher Trait-Ängstlichkeit, während dieser in der Gruppe mit geringer Trait-Ängstlichkeit nicht vorlag.

Derartige Befunde zeigen also einen Nutzen für die Verwendung mathematischer Modelle, hier des Diffusionsmodells, wenn es darum geht, sensitivere Maße zu erzeugen, die einen interessanten Prozess reiner abbilden können, als es einfache Unterschiede in mittleren RTs tun.

4.8 Zusammenfassung zeitunabhängiger Diffusionsmodelle

Wir haben in diesem Kapitel Random Walks zu einfachen Diffusionsmodellen weiterentwickelt. Damit haben wir ein sehr gutes Modell zur Entstehung korrekter und fehlerhafter RTs, mit welchem eine Reihe empirischer Phänomene detailierter analysiert werden kann, als mit reinen Mittelwerten. Insbesondere die Schätzung von Parametern für einzelne Versuchspersonen und deren weitere statistische Analyse hilft weiter, um Ursachen für Mittelwertunterschiede genauer zu verstehen. Im folgenden Kapitel befassen wir uns schließlich noch mit einem speziellen Diffusionsmodell, welches für Konfliktaufgaben entwickelt wurde. Wir werden auch sehen, dass das bisher behandelte Diffusionsmodell manche empirischen Phänomene von Konfliktaufgaben nicht erklären kann.

5 Ein Diffusionsmodell für Konfliktaufgaben

In diesem letzten Kapitel betrachten wir nun noch allgemeinere Formen von Diffusionsmodellen und fokussieren dabei auf ein Modell, welches insbesondere zur Modellierung von Konfliktaufgaben wie der Simon- oder der Eriksen-Flanker Aufgabe geeignet ist und verbale Dual-Route Theorien in ein mathematisches Modell übersetzt.

Der wesentliche Unterschied zum Modell aus Kapitel 4 ist der, dass wir nun zeitabhängige Parameter betrachten, also Parameter, die innerhalb eines Diffusionsprozesses ihren Wert verändern (siehe z.B. Heath, 1992). Dabei beschränken wir uns hier auf den Fall einer zeitabhängigen Drift Rate.

5.1 Konfliktaufgaben, Dual-Route Modelle und Verteilungsanalysen

5.1.1 Konfliktaufgaben und Dual-Route Modelle

Konfliktaufgaben wie die Flanker- oder die Simon-Aufgabe teilen eine Gemeinsamkeit: es gibt immer ein aufgabenrelevantes Merkmal, auf die eine Versuchsperson reagieren soll, und ein aufgabenirrelevantes Merkmal, welches die Versuchsperson ignorieren soll. Angenommen, die Versuchsperson soll auf einen Buchstaben reagieren, z.B. auf ein H mit einer linken und auf ein S mit einer rechten Reaktion.

• Bei einer Flanker-Aufgabe wird dieser Stimulus (das Target) i.d.R. zentral präsentiert. Links und rechts vom Target erscheinen andere Buchstaben, die Flanker, die entweder die gleiche oder die andere Reaktion als das Target erfordern würden, aber eben irrelevant für die Bearbeitung der Aufgabe sind. Erfordern sie die gleiche Reaktion, ist der Trial kongruent (z.B. bei der Zeichenkette SSS), ansonsten ist er inkongruent (z.B. bei der Zeichenkette HSH).
• Bei einer Simon-Aufgabe wird der Stimulus links oder rechts präsentiert, wobei die Lokation für die Bearbeitung der Aufgabe irrelevant ist. Erscheint der Stimulus auf der Seite, auf der die (korrekte) Reaktion abgegeben werden soll, ist der Trial kongruent, ansonsten ist er inkongruent.

Der Standardbefund dabei ist der Kongruenzeffekt: In kongruenten Trials sind die mittleren RTs kürzer (und es werden weniger Fehler gemacht) als in inkongruenten Trials. Typischerweise werden zur Erklärung von Kongruenzeffekten gern Dual-Route Modelle herangezogen (z.B. Kornblum et al., 1990). Die Grundidee ist in Abbildung 5.1 visualisiert: Ein Stimulus wird auf zweierlei Arten verarbeitet. Zum einen gibt es die kontrollierte, intentionale Reaktionsauswahl; zum anderen aktiviert ein Stimulus auch automatisch bestimmte Merkmale (z.B. aktiviert die Lokation eines Stimulus ein räumlich kongruentes Reaktionsmerkmal). Letztere Aktivierung ist transient, d.h. sie nimmt mit der Zeit wieder ab. Da sie aber mit der kontrollierten Reaktionsauswahl interagiert, wird die Reaktionsauswahl in kongruenten Trials erleichtert und in inkongruenten Trials arbeitet die automatische Aktivierung gegen die kontrollierte Reaktionsauswahl und erschwert diese folglich.

Abbildung 5.1: Illustration eines (einfachen) Dual-Route Modells

5.1.2 Verteilungsanalysen bei Konfliktaufgaben: CAFs und Delta-Funktionen

Die Standardauswertung bezieht sich auf gemittelte Daten und nutzt daher nicht die gesamte Information der Daten. Im Detail unterscheiden sich die Kongruenzeffekte zwischen Simon- und Flanker-Aufgaben nämlich tatsächlich. Dazu müssen wir sog. Delta-Funktionen anschauen, die die Kongruenzeffekte in Abhängigkeit vom RT-Level ausdrücken.

Der Abbildung 5.2a zeigt die CDF (also die Verteilungsfunktionen) mit den Quantilen 10%, 20%, …, 90%, getrennt nach kongruenten und inkongruenten Trials (z.B. einer Flanker-Aufgabe). Die horizontalen gepunkteten Linien, die die grünen und roten Kreise verbinden, sind demnach die Kongruenzeffekte für ein bestimmtes Quantil. Genau diese Differenzen werden in einer Delta-Funktion dann in Abhängigkeit der (gemittelten) RT-Quantile dargestellt (siehe Abb. 5.2b; schwarze Punkte). Da der Kongruenzeffekt in diesem Beispiel mit länger-werdenden RTs größer wird, liegt eine positive Delta-Funktion vor. Abbildung 5.2c zeigt schließlich eine Conditional Accuracy Function (CAF). Um eine CAF zu konstruieren, werden die RTs (separat für jede Kongruenzbedingung) in sog. Bins sortiert. Arbeiten wir mit fünf Bins, dann enthält Bin 1 die 20% kürzesten RTs usw. Pro Bin wird die relative Häufigkeit der korrekten Reaktionen, also die Accuracy, berechnet. Ein typisches Muster ist, dass die Accuracy in kongruenten Bedingungen unabhängig vom Bin (also der RT) ist, während in inkongruenten Bedingungen vermehrt schnelle Fehler vorkommen (also die Accuracy in den ersten Bins geringer ist und dann steigt).

Abbildung 5.2: Illustration der (a) Verteilungsfunktion, der (b) resultierenden Delta-Funktion und der (c) Conditional Accuaracy Function (CAF) einer typischen Eriksen-Flanker Aufgabe. In Panel (b) sind zur Illustration weitere Delta-Funktionen dargestellt, wie sie für andere Konfliktaufgaben auftreten können.

Derartige positive Delta-Funktionen sind tatsächlich eher der Standard (und werden nach vielen RT-Modellen auch erwartet; Wagenmakers & Brown, 2007). Interessanter ist daher, dass Daten von Simon-Aufgaben sehr häufig einen anderen Verlauf der Delta-Funtion aufweisen und der Simon-Kongruenzeffekt mit länger werdenden RTs kleiner wird, eine negative Delta-Funktion also.

Wir können anhand der Originaldaten aus der Arbeit von Ulrich et al. (2015) uns diese Unterschiede einmal im Detail anschauen. Die Daten sind im Pake dRiftDM verfügbar und können wie folgt abgerufen werden:

simon_data <- ulrich_simon_data # Simon Daten 
flanker_data <- ulrich_flanker_data # Flanker Daten

head(simon_data) # Die ersten Einträge anzeigen

Dieser Datensatz illustriert auch das Format, in welchem dRiftDM empirische Daten benötigt, um Diffusionsmodelle zu fitten:

• ID: identifiziert typischerweise die Versuchspersonen
• RT: Reaktionszeit in Sekunden
• Error: 0 = korrekter Trial, 1 = fehlerhafter Trial
• Cond: hier werden nun die verschiedenen Bedingungen unterschieden und gemeint ist mit comp und incomp kongruent bzw. inkongruent (oft werden diese Bedingungen auch kompatibel bzw. inkompatibel genannt)

Abbildung 5.3 illustriert die Verteilungsfunktionen, die Delta-Funktionen und die CAFs für die Flanker-Daten (obere Reihe, a-c) und die Simon-Daten (untere Reihe, d-f). Trotz aller konzeptueller Ähnlichkeit der beiden Aufgaben, scheinen sie sich irgendwie auf eine gewichtige Art und Weise zu unterscheiden (siehe v.a. die Delta-Funktionen). Dieser Unterschied würde aber in einfachen aggregierten Maßen gar nicht klar werden.

Abbildung 5.3: Illustration der CDF, Delta-Funktion und CAF von Flanker- (a-c) und Simon-Daten (d-f).

Heißt dies nun, dass wir von grundsätzlich verschiedenen Prozessen ausgehen sollten, die bei Simon- und Flanker-Aufgaben im Spiel sind? Obwohl naheliegend, wäre dies natürlich für die Theoriebildung nicht besonders sparsam.

Die Frage ist nun: Wie kann ein Diffusionsmodell die verschiedenen Verläufe von Delta-Funktionen erzeugen? Tatsächlich ist es so, dass das Diffusionsmodell wie wir es in Kapitel 4 behandelt haben, gar keine negativen Delta-Funktionen erzeugen kann.

5.2 Das Diffusion Model for Conflict tasks (DMC; Ulrich et al., 2015)

Von Hübner et al. (2010) und White et al. (2011) liegen zwei Diffusionsmodelle für Konfliktaufgaben vor; sie sind aber beide in erster Linie für Flanker-Aufgaben bzw. für empirische Muster mit steigenden Delta Funktionen konzeptualisiert worden.³ Das Diffusion Model for Conflict tasks (DMC; Ulrich et al., 2015) greift einen allgemeineren Ansatz auf und übersetzt quasi Dual-Route Modelle in ein Diffusionsmodell (Janczyk et al., 2025b, geben einen ausführlichen Überblick über das Modell, dessen Anwendungen und potenzielle Probleme).

Ausgangspunkt ist dabei die Überlegung, wie die erwarteten Verläufe der Evidenzakkumulation für die kontrollierte und die automatische Verarbeitung aussehen. (1) Für die kontrollierte Verarbeitung wird angenommen, dass diese mit einem linearen Zuwachs der Aktivierung einhergeht, also ganz wie es ein Standarddiffusionsmodell aus Kapitel 4 annimmt. Die entsprechende Drift Rate nennen wir \(\mu_c\) (wobei das angehängte “c” für “controlled” steht). (2) Für die automatische Verarbeitung hingegen wird angenommen, dass die durch sie erzeugte Aktierung zunächst zunimmt und dann, nach einem Maximum, wieder abnimmt und gegen Null geht. Das Argument ist erst einmal i.W., dass diese Aktivierung ja durch ein irrelevantes Merkmal erzeugt wird und daher nicht beibehalten werden sollte. Ein beispielhafter Befund, der für diese Annahme spricht ist eine Beobachtung aus CAFs. Bei kongruenten Durchgängen ist die Accuracy über alle RT-Bins gleich hoch. Bei inkongruenten Durchgängen hingegen nimmt die Accuracy mit dem RT-Bin zu. Dies deutet darauf hin, dass die Aktivierung durch das irrelevante Merkmal zu Beginn der Verarbeitung besonders stark ist und dann zunehmend geringer wird.

Nun wird eine Funktion benötigt, die eine derartige Form hat und eine Möglichkeit dazu ist die Funktion einer Gamma-Verteilung. Die Gamma-Verteilung ist erst einmal eine Dichtefunktion stetiger Zufallsvariablen und hat daher eine Fläche von 1. Wir multiplizieren die Funktion noch mit \(A\), um das Maximum (bzw. das Minimum, wenn \(A<0\)) der Funktion variieren zu können und reden dann von einer reskalierten Gamma-Verteilungsfunktion. Bezeichnen wir mit \(\mathbb{E}(\mathbf{X}_a(t))\) den erwarteten Verlauf der Aktivierung durch den automatischen Prozess, dann kann dieser geschrieben werden als \[ \mathbb{E}[\mathbf{X}_a(t)]=A\cdot e^{\left(-\frac{t}{\tau}\right)}\cdot \left[ \frac{t\cdot e}{(a-1)\cdot \tau}\right]^{(a-1)}\,. \] Diese Funktion hat drei Parameter. Der erste Parameter \(A\) beschreibt die Amplitude, also das Maximum der Funktion. Dieses Maximum wird erreicht zum Zeitpunkt \(t_{\text{max}}=(a-1)\cdot \tau\) (siehe hier für mehr Informationen, warum dies so ist). Setzt man daher \(a=2\) fest, dann beschreibt \(\tau\) den Zeitpunkt des Maximums \(A\). Abbildung 5.4 visualisiert drei Beispiele mit verschiedenen Werten für \(A\) and \(\tau\) (wobei \(a = 2\) fixiert wurde).

Abbildung 5.4: Beispiele (reskalierter) Gamma Verteilungsfunktionen

Was wäre nun die Drift Rate nach der oben bereits erwähnten Interpretation als 1. Ableitung der erwarteten Aktivierung nach der Zeit? Klar ist, dass für die Beispiele die Drift Rate erst positiv ist und nach dem Maximum dann negativ wird, dort ein Minimum erreicht und dann gegen Null geht. Technisch sieht die 1. Ableitung von \(\mathbb{E}[\mathbf{X}_a(t)]\) nach \(t\) so aus: \[ \mu_a(t)=\frac{d\mathbb{E}[\mathbf{X}_a(t)]}{dt}=A\cdot e^{\left(-\frac{t}{\tau}\right)}\cdot \left[ \frac{t\cdot e}{(a-1)\cdot \tau}\right]^{(a-1)}\cdot\left[\frac{a-1}{t}-\frac{1}{\tau}\right]\,. \] Weil die Drift Rate \(\mu_a\) von der Zeit abhängt, schreiben wir entsprechend \(\mu_a(t)\) für sie. Abbildung 5.5 visualisiert die drei Beispiele von gerade, nur dass jetzt die Drift Rate \(\mu_a(t)\) anstelle von \(\mathbb{E}[\mathbf{X}_a(t)]\) gezeigt wird.

Abbildung 5.5: Beispiele zeitabhängiger Drift Rates im Fall der automatischen Aktivierung von DMC.

Anmerkung: Wir haben hier den Fall, dass die Drift Rate ihr Vorzeichen wechselt. Dies widerspricht der inhaltlichen Interpretation aus Kapitel 4 natürlich, dass das Vorzeichen etwas mit der Stimulusqualität zu tun hat. Nach unserem Standpunkt gilt dies nur für zeitunabhängige Drift Rates (siehe dazu Lee & Sewell, 2024, vs. Janczyk et al., 2025a).

Die gesamte Drift Rate (und damit auch die gesamte erwartete Aktivierung, wenn man die Drift Rate über die Zeit hinweg akkumuliert bzw. integriert) ergibt sich nun durch superponieren (d.h. addieren) des kontrollierten und des automatischen Prozesses: \[ \mu(t)=\mu_c + \mu_a(t)\,. \] Wichtig ist noch, dass (angenommen \(\mu_c>0\)) \(A>0\) heißt, dass die automatische Aktivierung in die gleiche Richtung wie die kontrollierte Aktivierung geht. Dies entspricht einem kongruenten Trial. Inkongruente Trials werden dadurch implementiert, dass \(A<0\) gewählt wird, sodass die automatische Aktivierung entsprechend in die andere Richtung geht (die Drift Rate wäre dann auch erst negativ und wechselt ihr Vorzeichen nach positiv). Abbildung 5.6 illustriert die erwarteten Aktivierungen. In schwarz ist hier der lineare Prozess für die kontrollierten Verarbeitung gezeichnet (mit konstanter Drift Rate \(\mu_c\)). Als gestrichelte Linie in rot und grün ist der automatische Prozess in inkongruenten und kongruenten Trials gezeichnet. Zur Erinnerung: Der kontrollierte Prozess beschreibt die Verarbeitung des aufgaben-relevanten Merkmals (z.B. des zentralen Reizes bei einer Flanker-Aufgabe), wohingegen der automatische Prozess die Verarbeitung des aufgaben-irrelevanten Merkmals beschreibt (z.B. die seitlichen Distraktoren bei einer Flanker-Aufgabe). Die durchgezogene rote und grüne Linie ist der eigentliche Entscheidungsprozess, der sich aus der Addition des kontrollierten und des automatischen Prozesses ergibt. Man kann dabei leicht erkennen, dass kongruente Trials schneller die obere Boundary überschreiten als inkongruente Trials. Dieser Unterschied ist dann der Kongruenzeffekt.

Abbildung 5.6: Illustration der erwarteten Verläufe der Evidenzakkumulation für DMC.

Für das vollständige Modell wird nun noch, genau wie einfachen Diffusionsmodell aus Kapitel 4, eine Residualzeit zu dieser First-Passage Time addiert um RTs zu erhalten.

Neben der generellen Produktion eines Kongruenzeffektes kann DMC auch sehr leicht unterschiedliche Delta-Funktionen produzieren. Der wichtige Parameter ist dafür \(\tau\) (angenommen \(a=2\)), also wann die automatische Aktivierung durch das irrelevante Stimulusmerkmal sein Maximum erreicht: Passiert dies eher früh, dann resultiert eine negative Delta-Funktion; passiert dies eher spät, dann resultiert eine positive Delta-Funktion. Dies ist in Abbildung 5.7 für verschiedene Werte von \(\tau\) illustriert.

Abbildung 5.7: Illustration verschiedener Delta-Funktionen als Funktion von \(\tau\), also dem Maximum der automatischen Aktivierung. Die Amplitude wurde zwischen \(A=0.1\) und \(A=0.15\) variiert (linkes vs. rechtes Panel; \(a=2\) wurde überall angenommen).

Die Idee, dass die automatische Aktivierung schneller in einer Simon- als in einer Flanker-Aufgabe abnimmt kann durchaus theoretisch motiviert werden (siehe dazu Ulrich et al., 2015). Für unsere Zwecke hier ist noch eine entscheidende Frage: Was ist der Mechanismus wegen dem die automatische Aktivierung nach dem Erreichen des Maximums wieder sinkt? Hierzu gibt zwei, sich nicht ausschließende, Vorschläge:

• Passiver Zerfall: Da die Aktivation durch das irrelevante Stimulusmerkmal nicht benötigt wird, zerfällt diese nach einer Weile einfach wieder und wird dadurch geringer (z.B. Hommel, 1994).
• Aktive Inhibition: Von Ridderinkhoff (2002) stammt die Idee, dass die aufkommende Interferenz durch die automatische und irrelevante Aktivation dadurch reduziert wird, dass ein (eher aktiver) Mechanismus der Inhibition eingreift und die bereits aufgebaute Aktivation wieder reduziert.

DMC selber macht erst einmal keine Annahmen darüber, ob die spezifische Form des Zeitverlaufs der Aktivierung von Zerfall oder Inhibition herrührt.

5.3 Anwendungen des Modells

Die Gründe für die Anwendung speziellerer mathematischer Modelle wie DMC sind die gleichen wir bereits vorab ausgeführt. Wird in experimentellen Kontexten mit Konfliktaufgaben gearbeitet, ist DMC vermutlich eben auch eine bessere Approximation an die Realität als das einfache Diffusionsmodell es ist. Allerdings bringt dies auch einige Aspekte mit sich, die in der praktischen Arbeit berücksichtigt werden müssen.

5.3.1 Praktische Bestimmung der Modellparameter

Der wichtigste Aspekt ist, dass wir es mit einem Diffusionsmodell mit einem zeitabhängigen Parameter, der Drift Rate, zu tun haben. Viele der in Kapitel 4 aufgeführten Softwarelösungen funktionieren daher nicht mehr, sondern es werden spezielle Verfahren zur Vorhersage der CDF/PDF der korrekten und fehlerhaften RTs benötigt (siehe dazu auch Richter et al., 2023). Die einfachste Möglichkeit ist dann die Simulation, also das Schreiben von Funktionen, die ähnlich wie die in Kapitel 3 und 4 eingeführten Funktionen, eine große Anzahl an Trials simulieren können. Das R-Paket DMCfun (Mackenzie & Dudschig, 2021) nutzt diesen Ansatz, der zwar sehr flexibel, aber auch wenig effizient im Vergleich zu anderen Möglichkeiten ist. Das für uns wichtigste Paket ist dRiftDM (Koob et al., 2024), welches genau zum Zweck der Parameterschätzung im Fall zeitabhängiger Parameter geschaffen worden ist. Es nutzt schnelle und effiziente Methoden auf Basis partieller Differenzialgleichungen und enthält DMC bereits als eines der vor-implementierten Modelle. Mehr Informationen zur praktischen Arbeit mit DMC im Rahmen von dRiftDM gibt es hier.

5.3.2 Modellierung von Simon-Daten mit einem einfachen Diffusionsmodell versus mit DMC

Zu Demonstrationszwecken zeigen wir hier nun, wie DMC und das einfache Diffusionsmodell auf die oben betrachteten Simon-Daten aus Ulrich et al. (2015) mit dRiftDM fitten kann. Dafür benötigen wir zwei Dinge:

• Die Daten, also in unserem Fall die Simon-Daten mit den beiden Bedingungen “kongruent” und “inkongruent”, die in dem Datensatz allerdings mit „kompatibel“ (comp) und „inkompatibel“ (incomp) bezeichnet werden.
• Das jeweilige Modell, welches die Bedingungen der empirischen Daten umfasst.

Die Simon-Daten hatten wir oben bereits geladen und als simon_data gespeichert:

head(simon_data)

Nun benötigen wir noch die Modelle. In dRiftDM ist DMC bereits implementiert, und das entsprechende Modell kann mit Hilfe der Funktion dmc_dm() erzeugt werden:

dmc_model <- dmc_dm(dx = .005, dt = .005)
print(dmc_model)

## Class(es): dmc_dm, drift_dm
## 
## Current Parameter Matrix:
##        muc   b non_dec sd_non_dec  tau a    A alpha
## comp     4 0.6     0.3       0.02 0.04 2  0.1     4
## incomp   4 0.6     0.3       0.02 0.04 2 -0.1     4
## 
## Unique Parameters:
##        muc b non_dec sd_non_dec tau a A alpha
## comp   1   2 3       4          5   0 6 7    
## incomp 1   2 3       4          5   0 d 7    
## 
## Special Dependencies:
## A ~ incomp == -(A ~ comp)
## 
## Custom Parameters:
##        peak_l
## comp     0.04
## incomp   0.04
## 
## Deriving PDFs:
##   solver: kfe
##   values: sigma=1, t_max=3, dt=0.005, dx=0.005, nt=600, nx=400
## 
## Observed Data: NULL

Auch wenn uns die Details des Outputs hier nicht weiter kümmern sollen, lohnt sich ein Blick auf die Parameter-Labels unter Current Parameter Matrix. Hier finden sich nämlich alle wichtigen Parameter von DMC: die Drift Rate des kontrollierten Prozesses, muc, die Amplitude des automatischen Prozesses, A, sowie der Zeitpunkt des Maximums des automatischen Prozesses, tau, (da hier a = 2). Außerdem ist DMC bereits für die Simon-Daten nutzbar, da es die Bedingungen comp und incomp als Bedingungen spezifiziert hat (siehe die Zeilennamen im ersten Teil des Outputs).

Das einfache Diffusionsmodell kann mit Hilfe der Funktion ratcliff_dm() erzeugt werden:

simple_ddm <- ratcliff_dm(dx = .005, dt = .005)
print(simple_ddm)

## Class(es): ratcliff_dm, drift_dm
## 
## Current Parameter Matrix:
##      muc   b non_dec
## null   3 0.6     0.3
## 
## Unique Parameters:
##      muc b non_dec
## null   1 2       3
## 
## Deriving PDFs:
##   solver: kfe
##   values: sigma=1, t_max=3, dt=0.005, dx=0.005, nt=600, nx=400
## 
## Observed Data: NULL

Damit dieses Modell auch für die Simon-Daten verwendet werden kann, müssen wir es noch leicht modifizieren, sodass es (a) die Bedingungen comp und incomp umfasst und (b) die Parameter zwischen den Bedingungen frei variieren können. Was hier genau passiert, soll uns an dieser Stelle nicht zu sehr kümmern und ist dRiftDM-spezifische Syntax (mehr Informationen hierzu findet sich hier):

conds(simple_ddm) <- c("comp", "incomp")

## resetting parameter specifications

simple_ddm <- modify_flex_prms(
  simple_ddm,
  instr = " ~ "
)
print(simple_ddm)

## Class(es): ratcliff_dm, drift_dm
## 
## Current Parameter Matrix:
##        muc   b non_dec
## comp     3 0.6     0.3
## incomp   3 0.6     0.3
## 
## Unique Parameters:
##        muc b non_dec
## comp     1 3       5
## incomp   2 4       6
## 
## Deriving PDFs:
##   solver: kfe
##   values: sigma=1, t_max=3, dt=0.005, dx=0.005, nt=600, nx=400
## 
## Observed Data: NULL

Der wichtigste Teil dieses Outputs ist, dass es nun auch die beiden relevanten Bedingungen comp und incomp gibt und die drei Paramerer muc, b und non_dec jeweils getrennt für beide Bedingungen geschätzt werden. Dies ist nötig, da es ohne mindestens einen Parameter der sich zwischen kongruenten und inkongruenten Bedingungen unterscheidet, unmöglich wäre, den Kongruenzeffekt zu erzeugen.

Im letzten Schritt müssen wir nun die beiden Modelle auf die Simon-Daten fitten. Hierzu nutzen wir die Funktion estimate_model_ids(). Intern wird hierbei das jeweilige Modell auf die Daten jeder einzelnen Versuchspersonen mittels MLE gefittet. Als Optimierungsalgorithmus kommt hierbei (bounded) Nelder–Mead zum Einsatz (vgl. Kap. 2) :

# DMC fitten
estimate_model_ids(
  drift_dm_obj = dmc_model, # das Modell
  obs_data_ids = simon_data, # die Daten
  lower = c(
    muc = 1, b = .3, non_dec = .1, sd_non_dec = .005, tau = .03,
    A = .01, alpha = 2
  ), # untere Grenze des Suchraums
  upper = c(
    muc = 6, b = .9, non_dec = .5, sd_non_dec = .050, tau = .12,
    A = .15, alpha = 9
  ), # obere Grenze des Suchraums
  fit_procedure_name = "dmc_fit", # eine Bezeichnung für den gesamten Fit
  fit_path = getwd(), # hier werden die Daten gespeichert
  use_de_optim = FALSE, # bitte kein DEoptim nehmen (da langsam)
  use_nmkb = TRUE # sondern lieber für nmkb (für die Demo)
)

# einfaches DDM fitten
estimate_model_ids(
  drift_dm_obj = simple_ddm, # das Modell
  obs_data_ids = simon_data, # die Daten
  lower = c(muc = 1, b = .3, non_dec = .1), # untere Grenze des Suchraums
  upper = c(muc = 6, b = .9, non_dec = .5), # obere Grenze des Suchraums
  fit_procedure_name = "simple_fit", # eine Bezeichnung für den gesamten Fit
  fit_path = getwd(), # hier werden die Daten gespeichert
  use_de_optim = FALSE, # bitte kein DEoptim nehmen (da langsam)
  use_nmkb = TRUE # sondern lieber für nmkb (für die Demo)
)

Nachdem alles durchgelaufen ist (dauert ca. 20 Minuten), können wir das Ergebnis laden. Um den Modellfit zu prüfen, schauen wir die CDFs, CAFs und Delta-Funktionen an (Abb. 5.8 zeigt dabei das Ergebnis von DMC, während Abb. 5.9 das Ergebnis des einfachen Diffusionsmodells zeigt):

# laden von DMC
dmc_fits <- load_fits_ids(fit_procedure_name = "dmc_fit")
stats_dmc <- calc_stats(
  dmc_fits, type = c("quantiles", "cafs", "delta_funs"), minuends = "incomp", 
  subtrahends = "comp"
)

# laden vom einfachen DDM
simple_ddm_fits <- load_fits_ids(fit_procedure_name = "simple_fit")

stats_simple_ddm <- calc_stats(
  simple_ddm_fits, type = c("quantiles", "cafs", "delta_funs"), minuends = "incomp", 
  subtrahends = "comp"
)

## Aggregating across ID
## Aggregating across ID
## Aggregating across ID

Abbildung 5.8: Qualitativer Modellfit von DMC (Simon-Daten)

## Aggregating across ID
## Aggregating across ID
## Aggregating across ID

Abbildung 5.9: Qualitativer Modellfit des einfachen Diffusionsmodells (Simon-Daten)

Zwei Punkte sind hierbei wichtig:

1. Der Fit für DMC ist sehr gut; insbesondere kann DMC die negative Delta-Funktion der Simon-Daten vorhersagen.

2. Der Fit des einfachen Diffusionsmodells ist nicht wirklich gut; insbesondere sieht man, dass das einfache Diffusionsmodell eine positive Delta-Funktion vorhersagt, was klar im Widerspruch zur negativen Delta-Funktion der Daten steht. Tatsächlich kann ein einfaches Diffusionsmodell mit sinnvollen Parameterkonstellationen eine negative Delta-Funktion nicht vorhersagen. Auch sehen wir, dass dieses Modell ohne weitere Parameter die schnellen Fehler nicht vorhersagt und die CAFs entsprechend konstant über alle Bins sind.

5.3.3 Anwendungen in der Kognitionspsychologie

Seit der Einführung in die Literatur erweist sich DMC als häufig genutztes Modell in der kognitionspsychologischen Forschung im Kontext von Konfliktaufgaben. Wir wollen hier nicht allzu sehr ins Detail gehen und nur zwei Arbeiten aus unserem Bereich kurz anführen (siehe Janczyk et al., 2025b, für ein ausführliches Review).

Zum einen haben Koob et al. (2023a) DMC verwendet, um die Ursachen sequenzieller Modulationen bei Konfliktaufgaben besser zu verstehen. Der Ansatz dabei war zu untersuchen, welche Parameter sich in Abhängigkeit von der Kongruenz im vorherigen Trial jeweils verändern. Mit verschiedenen Herangehensweisen und anhand einer Reihe empirischer Datensätze wurde gefolgert, dass eine wahrscheinliche Ursache die Reduktion der irrelevanten Aktivation nach inkongruenten Trials ausschlaggebend ist.

Zum anderen wurde die Idee von DMC auch auf kongruenz-ähnliche Situationen bei Doppelaufgaben angewendet (Koob et al., 2023b). Hier lag der Fokus u.a. auf dem zeitlichen Verlauf der durch den zweiten Stimulus (automatisch) verursachten Aktivation. Interessanterweise zeigte sich hierbei, dass die Aktivation mit fortschreitender Zeit eher monoton ansteigt, statt wie bei DMC wieder geringer zu werden. Dies ergibt aber durchaus Sinn, wenn man bedenkt, dass diese Aktivation im Kontext von Doppelaufgaben eben auch nicht irrelevant ist, sondern im weiteren Verlauf der Bearbeitung der Aufgabe relevant wird, um die zweite Reaktion abzugeben.

5.3.4 Inhibitorisches Defizit bei Parkinson

Wenngleich der Einsatz von Konfliktaufgaben im klinischen Bereich nicht unüblich ist, werden bisher oft nur Standard Diffusionsmodelle zur modellbasierten Auswertung der Daten eingesetzt.⁴ Wir betrachten hier eine Studie im Kontext der Parkinson-Krankheit.

Parkinson wird häufig mit Problemen im Bereich Handlungskontrolle und -initiierung in Verbindung gebracht. Patient:innen mit Parkinson zeigen typischerweise längere RTs in Konfliktaufgaben als gesunde Kontrollpersonen (z.B. Praamsta et al., 1998). Dieser Befund wurde mit Defiziten in inhibitorischer Kontrolle, aber auch auch mit Schwierigkeiten zielgerichteter, kontrollierter Verarbeitung erklärt. Vor diesem Hintergrund eignen sich Simon-Aufgaben und die modell-basierte Auswertung der Daten zur detaillierteren Untersuchung. Behavioral zeigt sich zunächst, dass zwischen Patient:innen und Kontrollpersonen kein Unterschied darin besteht, dass in CAFs der Anteil korrekter Reaktionen für inkongruente Trials mit dem RT-Quantil zunimmt. Beide Gruppen sind also offenbar vergleichbar anfällig für die Aktivation durch die irrelevante Lokation. Allerdings fällt die Delta-Funktion bei Patient:innen weniger steil ab, als bei gesunden Personen. Dies wird vor dem Hintergrund der Theorie von Ridderinkhoff (2002) mit einem inhibitorischen Defizit erklärt. Eine passende Dopamin-Medikation kann darüber hinaus die Delta-Funktion wieder vergleichbar(er) mit der gesunder Personen machen (Wylie et al., 2010).

Servant et al. (2018) haben entsprechend in ihrer Studie DMC zur Re-Analyse von Daten einer Simon-Aufgabe genutzt. In ihrer Analyse gehen sie zunächst davon aus, die Reduktion der automatischen Aktivierung würde durch Inhibition erzeugt (Ridderinkhoff, 2002). Sie definieren dafür eine Variable, die sie Suppression Strength (\(S_{\text{strength}}\)) nennen, und die i.W. die Zeit zwischen Beginn der Inhibition (d.h. dem Zeitpunkt des Maximums der Gamma-Funktion) und dem Zeitpunkt, zu dem 90% der automatischen Aktivierung erfolgt sind, erfasst. Je länger dieser Zeitraum ist, als desto weniger effizient kann die zugrunde liegende Inhibition interpretiert werden. Darüber hinaus kann \(A\), also die Amplitude, als Stärke der lokations-induzierten Aktivation betrachtet werden, sowie \(\mu_c\) als Effizienz kontrollierter Verarbeitung.

Wir beschränken uns hier auf den ersten re-analysierten Datensatz, der von einer Simon-Aufgabe mit einer Speed-Accuracy-Manipulation (durch entsprechende Instruktion) herrührt und aus einer Arbeit von van Wouwe (2014) stammt. Abbildung 2A in Servant et al. (2018) zeigt die CAFs (oben) und Delta-Funktionen (unten) getrennt für gesunde Kontrollpersonen (HC) und Parkinson-Patient:innen (PD) sowie für die Accuracy- (acc) und die Speed-Instruktion (spd). Während die CAFs sehr ähnlich für HC- und PD-Personen aussehen, sind die Delta-Funktionen weniger negativ für die PD-Personen und für die spd-Bedingung. In Bezug auf die, aus der DMC-Auswertung resultierenden, Parameter ergab sich, dass die Gruppen sich nicht in der Stärke der automatischen Aktivierung, \(A\), unterschieden. Allerdings setzt die Inhibition bei den Patient:innen sowohl später ein und sie dauert länger an, ist also weniger effizient. Auch die kontrollierte Drift Rate \(\mu_c\) war bei den Patient:innen kleiner als bei den Kontrollpersonen, auch wenn dieser Unterschied statistisch nicht signifikant war (Bei der Re-Analyse eines zweiten Datensatzes war dieser Effekt auch statistisch signifikant.).⁵

Insgesamt werten die Autor:innen ihre Ergebnisse so, dass sowohl die kontrollierte Verarbeitung als auch die inhibitorische Kontrolle zur Reduktion irrelevanter Aktivation bei Parkinson-Patient:innen defizitär sind. Diese Unterscheidung wäre bei Verwendung eines einfachen Diffusionsmodells nicht möglich gewesen und zeigt exemplarisch die Wichtigkeit der Verwendung passender und gut etablierter Modelle in der praktischen Anwendung, um zu tragfähigen theoretischen Schlussfolgerungen zu gelangen.

6 Literatur

Akaike, H. (1974). A new look at the statistical model identification. IEEE transactions on automatic control, 19(6), 716-723.

Alexandrowicz, R. W. (2020). The diffusion model visualizer: An interactive tool to understand the diffusion model parameters. Psychological Research, 84, 1157-1165.

Bar-Haim, Y., Lamy, D., Pergamin, L., Bakermans-Kranenburg, M. J., & van Ijzendoorn, M. H. (2007). Threat-related attentional bias in anxious and nonanxious individuals: A meta-analytic study. Psychological Bulletin, 133, 1–24.

Blurton, S.P., Kesselmeier, M. & Gondan, M. (2012). Fast and accurate calculations for cumulative first-passage time distributions in Wiener diffusion models. Journal of Mathematical Psychology, 56, 470-475.

Box, G.E.P. (1976). Science and statistics. Journal of the American Statistical Association, 71(356), 791–99.

Cox, D.R. & Miller, H.D. (1965). The theory of stochastic processes. Chapman & Hall.

Diederich, A. & Busemeyer, J.R. (2003). Simple matrix methods for analyzing diffusion models of choice probability, choice response time, and simple response time. Journal of Mathematical Psychology, 47(3), 304-322.

Diederich, A. & Mallahi-Karai, K. (2018). Stochastic methods for modeling decision-making.” In W. H. Batchelder, H. Colonius, & E. N. Dzhafarov (eds.), New handbook of mathematical psychology (p. 1-70). Cambridge University Press.

Feller, W. (1968). An introduction to probability theory and its applications: Vol. 1 (3rd ed.). Wiley.

Groulx, J.T., Harding, B. & Cousineau, D. (2000). The EZ Diffusion Model: An overview with derivation, software, and an application to the same-different task. The Quantitative Methods for Psychology, 16, 154-174.

Heath, R. A. (1992). A general nonstationary diffusion model for two-choice decision-making. Mathematical Social Sciences, 23, 283–309.

Heathcote, A., Brown, S. & Mewhort, D. J. (2000). The power law repealed: The case for an exponential law of practice. Psychonomic Bulletin & Review, 7, 185-207.

Henze, N. (2024). Asymptotic stochastics: An introduction with a view towards statistics. Springer.

Hommel, B. (1994). Spontaneous decay of response-code activation. Psychological Research, 56, 261–268

Hübner, R., Steinhauser, M. & Lehle, C. (2010). A dual-stage two-phase model of selective attention. Psychological Review, 117, 304-322.

Hübner, R. & Töbel, L. (2019). Conflict resolution in the Eriksen flanker task: Similarities and differences to the Simon task. PLoS One, 14, e0214203.

Janczyk, M., Mackenzie, I.G. & Koob, V. (2025a). A comment on the Revised Diffusion Model for Conflict Tasks (RDMC). Psychonomic Bulletin & Review, 32, 690-704.

Janczyk, M., Mackenzie, I.G., Ulrich, R., & Koob, V. (2025b). Ten years Diffusion Model for Conflict (DMC) tasks: Theoretical foundations, applications, practical recommendations, and open challenges. Manuscript submitted for publication. Preprint available at https://osf.io/preprints/psyarxiv/rgepd_v1

Koob, V., Mackenzie, I. G., Ulrich, R., Leuthold, H., & Janczyk, M. (2023a). The role of task-relevant and task-irrelevant information in congruency sequence effects: Applying the diffusion model for conflict tasks. Cognitive Psychology, 140, 101528.

Koob, V., Richter, T., Ulrich, R. & Janczyk, M. (2024). An introduction and tutorial to fitting (time-dependent) diffusion models with the R-package dRiftDM. Manuscript under revision. Preprint available at https://osf.io/preprints/osf/3t2vf.

Koob, V., Ulrich, R., & Janczyk, M. (2023b). Response activation and activation-transmission in response-based backward crosstalk: Analyses and simulations with an extended diffusion model. Psychological Review, 130, 102–136

Kornblum, S., Hasbroucq, T. & Osman, A. (1990). Dimensional overlap: Cognitive basis for stimulus-response compatibility – A model and taxonomy. Psychological Review, 97, 253–270.

Lee, P.-S. & Sewell, K., D. (2024). A revised diffusion model for conflict tasks. Psychonomic Bulletin Review, 31, 1–31.

Lerche, V., von Krause, M., Voss, A., Frischkorn, G., Schubert, A.-L. & Hagemann, D. (2020). Diffusion modeling and intelligence: Drift rates show both domain-general and domain-specific relations with intelligence. Journal of Experimental Psychology: General, 149, 2207-2249.

Lerche, V. & Voss, A. (2016). Model complexity in diffusion modeling: Benefits of making the model more parsimonious. Frontiers in Psychology, 7, 1324.

Lewandowsky, S. & Farrell, S. (2013). Computational modeling of cognition: Principles and practice. Sage Publications Inc.

Liesefeld, H.R. & Janczyk, M. (2019). Combining speed and accuracy to control for speed-accuracy tradeoffs (?). Behavior Research Methods, 51, 40-60.

Liesefeld, H.R. & Janczyk, M. (2023). Same same but different: Subtle but consequential differences between two measures to linearly integrate speed and accuracy (LISAS vs. BIS). Behavior Research Methods, 55, 1175-1192.

Luce, R.D. (1995). Four tensions concerning mathematical modeling in psychology. Annual Review of Psychology, 46, 1-26.

MacCallum, R.C. (2003). Working with imperfect models. Multivariate Behavioral Research, 38(1), 13-39.

Mackenzie, I. G. & Dudschig, C. (2021). DMCfun: An R package for fitting Diffusion Model of Conflict (DMC) to reaction time and error rate data. Methods in Psychology, 5, 100074.

MacLeod, C., & Mathews, A. (1991). Biased cognitive operations in anxiety: Accessibility of information or assignment of processing priorities. Behaviour Research and Therapy, 29, 599–610.

Matzke, D. & Wagenmakers, E.-J. (2009). Psychological interpretation of the ex-Gaussian and shifted Wald parameters: A diffusion model analysis. Psychonomic Bulletin & Review, 16, 798–817.

Mosteller, F., Fienberg, S.E., & Rourke, R.E.K. (1983). Beginning statistics with data analysis. Courier corporation.

Nelder, J. A. & Mead, R. (1965). A simplex method for function minimization. Computer Journal, 7, 308-313.

Oberauer, K. & Lewandowsky, S. (2019). Addressing the theory crisis in psychology. Psychonomic Bulletin & Review, 26, 1596–1618.

Praamsta, P., Stegeman, D.F., Cools, A.R. & Horstink, M.W. (1998). Reliance on external cues for movement initiation in Parkinson’s disdease. Evidence from movement-related potentials. Brain, 121, 167-177.

Ratcliff, R. (1978). A theory of memory retrieval. Psychological Review, 85, 59-108.

Ratcliff, R. (2013). Parameter variability and distributional assumptions in the diffusion model. Psychological Review, 120, 281-292.

Ratcliff, R. & McKoon, G. (2008). The diffusion decision model: Theory and data for two-choice decision tasks. Neural Computation, 20, 873–922.

Ratcliff, R. & Rouder, J. N. (1998) Modeling response times for two-choice decisions. Psychological Science, 9, 347–356.

Ratcliff, R., Smith, P. L. Brown, S. & McCoon, G. (2016). Diffusion decision model: Current issues and history. Trends in Cognitive Sciences, 20, 260-281.

Ratcliff, R., Thapar, A. & McCoon, G. (2006). Aging and individual differences in rapid two-choice decisions. Psychonomic Bulletin and Review, 13, 626-635.

Richter, T., Ulrich, R., & Janczyk, M. (2023). Diffusion models with time-dependent parameters: An analysis of computational effort and accuracy of different numerical methods. Journal of Mathematical Psychology, 114, 102756.

Ridderinkhof, K. R. (2002). Activation and suppression in conflict tasks, empirical clarification through distributional analyses. In W. Prinz & B. Hommel (Eds.), Common mechanisms in perception and action. Attention and performance XIX (pp. 494–519). Oxford University Press.

Rosenbaum, D. A. & Janczyk, M. (2019). Who is or was E. R. F. W. Crossman, the champion of the Power Law of Learning and the developer of an influential model of aiming? Psychonomic Bulletin and Review, 26, 1449-1463.

Schmiedek, F., Oberauer, K., Wilhelm, O., Süß, H.-M. & Wittmann, W. W. (2007). Individual differences in components of reaction time distributions and their relations to working memory and intelligence. Journal of Experimental Psychology: General, 136, 414-429.

Schmitz, F. & Voss, A. (2012). Decomposing task-switching costs with the diffusion model. Journal of Experimental Psychology: Human Perception and Performance, 38, 222-250.

Shiffrin, R.M., Lee, M.D., Kim, W. & Wagenmakers, E.-J. (2008). A survey of model evaluation approaches with a tutorial on hierarchical Bayesian methods. Cognitive Science, 32(8), 1248–84.

Spielberger, C. D., Gorsuch, R. L. & Lushene, R. E. (1970). Manual for the State–Trait Anxiety Inventory. Consulting Psychologists Press.

Schubert, A.-L., Frischkorn, G. T., Hagemann, D. & Voss, A. (2016). Trait characteristics of diffusion model parameters. Intelligence, 4, 7.

Schwarz, G. (1978). Estimating the dimension of a model. The annals of statistics, 6(2), 461-464.

Schwarz, W. (2022). Random walk and diffusion models. An introduction for life and behavioral scientists. Springer.

Servant, M, van Wouwe, N., Wylie, S.A. & Logan, G.D. (2018). A model-based quantification of action control deficits in Parkinson’s disease. Neuropsychologia, 111, 26-35.

Smith, P. L. & Ratcliff, R. (2015). An introduction to the diffusion model of decision making. In B. U. Forstmann & E.-J. Wagenmakers (Eds.), An introduction to model-based cognitive neuroscience (S. 49-70). Berlin: Springer.

Storn, R. & Price, K. (1997). Differential Evolution - A simple and efficient heuristic for global optimization over continuous spaces. Journal of Global Optimization, 11, 341–359.

Ulrich, R. (2009). Uncovering unobservable cognitive mechanisms: The contribution of mathematical models. In F. Rösler, C. Ranganath, B. Röder, and R. H. Kluwe (Eds., Neuroimaging of human memory: Linking cognitive processes to neural systems (S. 25–41). Oxford University Press.

Ulrich, R., Schröter, H., Leuthold, H. & Birngruber, T. (2015). Automatic and controlled stimulus processing in conflict tasks: Superimposed diffusion processes and delta functions. Cognitive Psychology, 78, 148-174.

van Wouwe, N.C., van den Wildenberg, W.P., Classen, D.O., Kanoff, K., Bashore, T.R. & Wylie, S.A. (2014). Speed pressure in conflict situations impedes inhibitory action control in Parkinson’s disease. Biological Psychology, 101, 44-60.

Voss, A., Nagler, M. & Lerche, V. (2013). Diffusion models in experimental pychology. Experimental Psychology, 60, 385-402.

Voss, A., Rothermund, K. & Voss, J. (2004). Interpreting the parameters of the diffusion model: An empirical validation. Memory & Cognition, 32, 1206-1220.

Voss, A., & Voss, J. (2007). Fast-dm: A free program for efficient diffusion model analysis. Behavior Research Methods, 39, 767–775.

Voss, A., Voss, J., & Lerche, V. (2015). Assessing cognitive processes with diffusion model analyses: A tutorial based on fast-dm-30. Frontiers in Psychology, 6, 336.

Wagenmakers, E. J. & Brown, S. (2007). On the linear relation between the mean and the standard deviation of a response time distribution. Psychological Review, 114, 830-841.

Wagenmakers, E.-J., Van Der Maas, H.L.J. & Grasman, R.P.P.P. (2007). An EZ-diffusion model for response time and accuracy. Psychonomic Bulletin & Review 14, 3–22.

White, C. N., Ratcliff, R. & Starns, J. J. (2011). Diffusion models of the flanker task: Discrete versus gradual attentional selection. Cognitive Psychology, 63, 210-238.

White, C. N., Ratcliff, R., Vasey, M. W. & McKoon, G. (2010a). Anxiety enhances threat processing even without competition among multiple inputs: A diffusion model analysis. Emotion, 10, 662-677.

White, C. N., Ratcliff, R., Vasey, M. W. & McKoon, G. (2010b). Using diffusion models to understand clinical disorders. Journal of Mathematical Psychology, 54, 39-52.

Wylie, S.A., Ridderinkhoff, K.R., Bashore, T.R. & van den Wildenberg, W.P. (2010). The effect of Parkinson’s disease on the dynamics of on-line and proactive cognitive control during action selection. Journal of Cognitive Neuroscience, 22, 2058-2073.

---

1. Das ist hier etwas vereinfacht gesagt. Tatsächlich variieren die entsprechenden Algorithmen die Parameter so, dass \(Q\) auch durchaus wieder steigen kann.↩︎

2. Die Simulation stochastischer Prozesse in Modellen ist tatsächlich für viele Probleme sehr verbreitet und oft zunächst das Mittel der Wahl. Mit fortschreitender Verwendung und Analyse des Problems werden dann aber mitunter die zugrunde liegenden mathematischen Probleme gelöst, sodass genauere und effizientere Lösungen möglich werden. Das bedeutet nicht, dass es immer exakte und geschlossene Lösungen geben muss; manchmal werden auch “nur” die Methoden zur Annäherung an die exakte Lösung besser und effizienter (vgl. z.B. Richter et al., 2023).↩︎

3. Es gibt aber auch eine Erweiterung des Modell von Hübner et al. (2010), um negative Delta-Funktionen zu produzieren (siehe Hübner & Töbel, 2019).↩︎

4. Inwiefern die dann nicht mögliche Trennung von kontrollierter und automatischer Verarbeitung zu falschen Schlussfolgerungen führen kann ist derzeit unklar. Dies gilt aber auch für nicht-klinische Studien, bei denen Konfliktaufgaben mit Standard Diffusionsmodellen ausgewertet werden.↩︎

5. Die Speed-Accuracy Manipulation hatte zudem Auswirkungen auf den Bondary-Parameter \(b\), aber auch auf \(\mu_c\), \(A\) und \(t_0\).↩︎