Use this Skill to analyze reaction time data: outlier removal, ex-Gaussian distribution fitting, EZdiff DDM parameters, and HDDM-equivalent drift-diffusion models with PyMC.
Use this Skill to analyze reaction time data: outlier removal, ex-Gaussian distribution fitting, EZdiff DDM parameters, and HDDM-equivalent drift-diffusion models with PyMC.
TL;DR — Complete pipeline for reaction time (RT) data quality control,
ex-Gaussian distribution fitting, EZdiff DDM parameter extraction, and
full Bayesian drift-diffusion modeling with PyMC. Covers outlier removal,
Q-Q plots, Vincentile curves, delta plots, and condition comparison.
When to Use
Use this Skill when you need to:
Clean RT data from keyboard-response behavioral experiments
Fit ex-Gaussian distributions to characterize RT distributions per condition
Extract drift-diffusion model (DDM) parameters without full Bayesian fitting
(EZdiff; Wagenmakers et al., 2007)
Build a full Bayesian DDM with PyMC to compare conditions (accuracy, speed)
Visualize RT distributions with Q-Q plots, density plots, and delta plots
Compare conditions using Vincentile (quantile-average) curves
Background
Reaction times are right-skewed because fast responses are bounded by motor
minimum latency (~150 ms) while slow responses have no upper bound. Two
complementary frameworks model this shape:
Ex-Gaussian Distribution
The ex-Gaussian is the convolution of a Gaussian (mean μ, SD σ) and an
exponential (rate λ, mean τ = 1/λ):
μ (mu): the Gaussian component mean — reflects fast, Gaussian-like RTs
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
from typing importOptional, Dict, Tuple, List# ── Absolute and SD-based outlier removal ───────────────────────────────────defclean_rt_data(
df: pd.DataFrame,
rt_col: str = "rt",
accuracy_col: str = "correct",
min_rt: float = 150.0,
max_rt: float = 3000.0,
sd_cutoff: float = 2.5,
group_col: Optional[str] = None,
) -> Tuple[pd.DataFrame, Dict]:
"""
Remove RT outliers using absolute cutoffs and SD-based trimming.
Args:
df: DataFrame with one row per trial.
rt_col: Column name for reaction times (milliseconds).
accuracy_col: Column name for accuracy (1=correct, 0=error).
min_rt: Absolute lower bound (ms); RTs below are fast-guesses.
max_rt: Absolute upper bound (ms); RTs above are likely lapses.
sd_cutoff: Number of SDs for within-condition trimming.
group_col: Optional grouping column (e.g., 'condition') for
SD-based trimming within groups.
Returns:
Tuple of (cleaned_df, stats_dict).
"""
n_original = len(df)
# Absolute cutoffs
mask_abs = (df[rt_col] >= min_rt) & (df[rt_col] <= max_rt)
df_clean = df[mask_abs].copy()
n_after_abs = len(df_clean)
# SD-based trimming (within group if specified)if group_col and group_col in df_clean.columns:
groups = df_clean[group_col].unique()
keep_mask = pd.Series(True, index=df_clean.index)
for g in groups:
g_mask = df_clean[group_col] == g
g_rts = df_clean.loc[g_mask, rt_col]
mean_rt = g_rts.mean()
sd_rt = g_rts.std(ddof=1)
out_mask = (
(df_clean[rt_col] < mean_rt - sd_cutoff * sd_rt) |
(df_clean[rt_col] > mean_rt + sd_cutoff * sd_rt)
)
keep_mask[g_mask & out_mask] = False
df_clean = df_clean[keep_mask].copy()
else:
mean_rt = df_clean[rt_col].mean()
sd_rt = df_clean[rt_col].std(ddof=1)
df_clean = df_clean[
(df_clean[rt_col] >= mean_rt - sd_cutoff * sd_rt) &
(df_clean[rt_col] <= mean_rt + sd_cutoff * sd_rt)
].copy()
n_final = len(df_clean)
stats_dict = {
"n_original": n_original,
"n_removed_absolute": n_original - n_after_abs,
"n_removed_sd": n_after_abs - n_final,
"n_final": n_final,
"pct_removed": round(100 * (1 - n_final / n_original), 2),
"mean_rt_clean": round(df_clean[rt_col].mean(), 2),
"median_rt_clean": round(df_clean[rt_col].median(), 2),
"sd_rt_clean": round(df_clean[rt_col].std(ddof=1), 2),
"mean_accuracy": round(df_clean[accuracy_col].mean(), 4),
}
print(
f"RT cleaning: {n_original} → {n_final} trials "f"({stats_dict['pct_removed']:.1f}% removed)\n"f" Absolute cutoffs: {stats_dict['n_removed_absolute']} removed\n"f" SD trimming (±{sd_cutoff}): {stats_dict['n_removed_sd']} removed\n"f" Mean RT = {stats_dict['mean_rt_clean']} ms, "f"Accuracy = {stats_dict['mean_accuracy']:.1%}"
)
return df_clean, stats_dict
Step 2 — Ex-Gaussian Fitting and Q-Q Plot
from scipy.stats import exponnorm
deffit_exgaussian(
rts: np.ndarray,
label: str = "data",
plot: bool = True,
) -> Dict:
"""
Fit an ex-Gaussian distribution to RT data and produce a Q-Q plot.
scipy.stats.exponnorm parameterization:
K = tau / sigma (shape)
loc = mu (Gaussian mean, milliseconds)
scale = sigma (Gaussian SD)
Therefore:
mu = loc
sigma = scale
tau = K * scale
Args:
rts: 1-D array of reaction times (milliseconds, already cleaned).
label: Label for plot title and print output.
plot: Whether to generate the Q-Q plot.
Returns:
Dict with mu, sigma, tau, K, loc, scale, and fit statistics.
"""# MLE fit
K_fit, loc_fit, scale_fit = exponnorm.fit(rts, floc=None)
mu = loc_fit
sigma = scale_fit
tau = K_fit * scale_fit
mean_theoretical = mu + tau
var_theoretical = sigma ** 2 + tau ** 2# KS goodness-of-fit
ks_stat, ks_p = stats.kstest(rts, lambda x: exponnorm.cdf(x, K_fit, loc_fit, scale_fit))
params = {
"label": label,
"mu": round(mu, 3),
"sigma": round(sigma, 3),
"tau": round(tau, 3),
"K": round(K_fit, 4),
"loc": round(loc_fit, 3),
"scale": round(scale_fit, 3),
"mean_theoretical": round(mean_theoretical, 3),
"sd_theoretical": round(np.sqrt(var_theoretical), 3),
"ks_statistic": round(ks_stat, 4),
"ks_p": round(ks_p, 4),
}
print(
f"Ex-Gaussian fit [{label}]: μ={params['mu']:.1f}, "f"σ={params['sigma']:.1f}, τ={params['tau']:.1f} ms | "f"KS p={params['ks_p']:.3f}"
)
if plot:
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
# Density overlay
ax = axes[0]
rt_range = np.linspace(rts.min(), rts.max(), 500)
ax.hist(rts, bins=40, density=True, alpha=0.4, color="steelblue", label="Observed")
ax.plot(rt_range, exponnorm.pdf(rt_range, K_fit, loc_fit, scale_fit),
color="crimson", linewidth=2, label="Ex-Gaussian fit")
ax.set_xlabel("Reaction Time (ms)")
ax.set_ylabel("Density")
ax.set_title(f"Ex-Gaussian fit — {label}")
ax.legend()
# Q-Q plot vs ex-Gaussian
ax2 = axes[1]
n = len(rts)
probs = (np.arange(1, n + 1) - 0.5) / n
emp_quantiles = np.sort(rts)
theo_quantiles = exponnorm.ppf(probs, K_fit, loc_fit, scale_fit)
ax2.scatter(theo_quantiles, emp_quantiles, s=10, alpha=0.5, color="steelblue")
lims = [min(theo_quantiles.min(), emp_quantiles.min()),
max(theo_quantiles.max(), emp_quantiles.max())]
ax2.plot(lims, lims, "r--", linewidth=1.5, label="Identity line")
ax2.set_xlabel("Theoretical ex-Gaussian quantiles (ms)")
ax2.set_ylabel("Empirical quantiles (ms)")
ax2.set_title(f"Q-Q Plot — {label}")
ax2.legend()
fig.tight_layout()
plt.show()
return params
defvincentile_plot(
rt_dict: Dict[str, np.ndarray],
n_bins: int = 5,
output_path: Optional[str] = None,
) -> plt.Figure:
"""
Create Vincentile (quantile-average) curves for multiple conditions.
Vincentiles average participants' quantiles across the sample, preserving
distributional shape while allowing condition comparison.
Args:
rt_dict: Dict mapping condition label to array of RTs.
n_bins: Number of quantile bins (default 5 = quintiles).
output_path: Optional path to save figure.
Returns:
Matplotlib Figure.
"""
quantile_positions = np.linspace(0, 1, n_bins + 2)[1:-1]
fig, ax = plt.subplots(figsize=(8, 5))
colors = plt.cm.tab10(np.linspace(0, 0.6, len(rt_dict)))
for (label, rts), color inzip(rt_dict.items(), colors):
quantiles = np.quantile(rts, quantile_positions)
ax.plot(quantile_positions * 100, quantiles, marker="o",
color=color, linewidth=2, markersize=6, label=label)
ax.set_xlabel("Percentile")
ax.set_ylabel("Reaction Time (ms)")
ax.set_title("Vincentile Plot")
ax.legend()
ax.grid(alpha=0.3)
fig.tight_layout()
if output_path:
fig.savefig(output_path, dpi=150)
plt.show()
return fig
Step 3 — EZdiff DDM Parameter Extraction
from scipy.special import ndtri # probit functiondefezdiff(
mean_rt: float,
var_rt: float,
accuracy: float,
scaling_factor: float = 0.1,
) -> Dict:
"""
Extract DDM parameters (v, a, Ter) using the EZdiff method.
Reference: 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(1), 3–22.
Assumptions: equal-variance, unbiased starting point (z = a/2),
no inter-trial variability. Best used for quick parameter estimation
when full DDM fitting is not feasible.
Args:
mean_rt: Mean reaction time (seconds or ms — consistent units).
var_rt: Variance of reaction times (same units squared).
accuracy: Proportion correct (0 < Pc < 1; values near 0.5 or 1
require edge correction).
scaling_factor: DDM scaling constant (default 0.1 for standard units).
Returns:
Dict with drift rate v, boundary a, and non-decision time Ter.
"""# Edge correction: avoid Pc = 0.5 (v=0) or Pc = 1.0 (v→∞)if accuracy <= 0.5:
accuracy = max(accuracy, 0.501)
if accuracy >= 1.0:
accuracy = min(accuracy, 0.999)
s = scaling_factor
z_pc = ndtri(accuracy) # probit# Drift rate
v = (
np.sign(accuracy - 0.5) * s *
(z_pc / mean_rt - z_pc ** 3 / var_rt * mean_rt) ** 0.5if (z_pc / mean_rt - z_pc ** 3 / var_rt * mean_rt) >= 0else0.0
)
# Alternative closed-form (Wagenmakers 2007, eq. 3)# More numerically stable form:
L = ndtri(accuracy)
x = L * (L ** 2 * mean_rt - var_rt * L ** 2 / mean_rt - var_rt / mean_rt)
if x > 0:
v = np.sign(accuracy - 0.5) * s * (x ** 0.5) / np.sqrt(var_rt)
else:
v = 0.0# Boundary separation
a = s * L / v if v != 0else np.nan
# Non-decision time
Ter = mean_rt - a / (2 * v) if v != 0else np.nan
result = {
"v": round(float(v), 4),
"a": round(float(a), 4),
"Ter": round(float(Ter), 4),
"mean_rt": mean_rt,
"var_rt": var_rt,
"accuracy": accuracy,
}
print(
f"EZdiff: v={result['v']:.3f}, a={result['a']:.3f}, "f"Ter={result['Ter']:.3f} | Pc={accuracy:.3f}"
)
return result
defezdiff_by_condition(
df: pd.DataFrame,
rt_col: str = "rt",
accuracy_col: str = "correct",
condition_col: str = "condition",
rt_units: str = "ms",
) -> pd.DataFrame:
"""
Apply EZdiff to each condition in a DataFrame.
Args:
df: Trial-level DataFrame.
rt_col: Column of reaction times.
accuracy_col: Column of accuracy (1/0).
condition_col: Column identifying conditions.
rt_units: 'ms' or 's' — if 'ms', divides by 1000 before fitting.
Returns:
DataFrame with one row per condition and DDM parameters.
"""
rows = []
for cond, group in df.groupby(condition_col):
rts = group[rt_col].values
if rt_units == "ms":
rts_s = rts / 1000.0else:
rts_s = rts
mrt = rts_s.mean()
vrt = rts_s.var(ddof=1)
pc = group[accuracy_col].mean()
params = ezdiff(mrt, vrt, pc)
params["condition"] = cond
params["n_trials"] = len(group)
rows.append(params)
result_df = pd.DataFrame(rows).set_index("condition")
print("\nEZdiff parameters by condition:")
print(result_df[["v", "a", "Ter", "n_trials"]].round(3))
return result_df
Advanced Usage
PyMC Drift-Diffusion Model
import pymc as pm
import numpy as np
deffit_ddm_pymc(
rt_correct: np.ndarray,
rt_error: np.ndarray,
draws: int = 2000,
tune: int = 1000,
target_accept: float = 0.90,
seed: int = 42,
) -> pm.backends.base.MultiTrace:
"""
Fit a simple Wiener drift-diffusion model using PyMC.
Models both correct and error RTs via the Wiener first-passage time
distribution. Assumes equal-variance, zero between-trial variability.
Priors (weakly informative):
v ~ Normal(0, 2) — drift rate
a ~ HalfNormal(1) — boundary separation (positive)
Ter ~ Uniform(0.1, 0.5) — non-decision time (seconds)
z ~ Uniform(0.3, 0.7) — relative starting point (0.5 = unbiased)
Args:
rt_correct: Array of RTs for correct responses (seconds).
rt_error: Array of RTs for error responses (seconds).
draws: Number of posterior samples per chain.
tune: Number of tuning steps.
target_accept: NUTS target acceptance rate.
seed: Random seed.
Returns:
PyMC InferenceData object with posterior samples.
Notes:
- Requires PyMC >= 5.0 and the hssm package or manual Wiener
likelihood. Below uses a simplified normal approximation for
demonstration; for production use hssm (pip install hssm).
- Do NOT use this function with hard-coded API keys or credentials.
"""
all_rts = np.concatenate([rt_correct, rt_error])
responses = np.concatenate([
np.ones(len(rt_correct)),
np.zeros(len(rt_error))
])
with pm.Model() as ddm_model:
# Priors
v = pm.Normal("v", mu=0, sigma=2)
a = pm.HalfNormal("a", sigma=1)
Ter = pm.Uniform("Ter", lower=0.05, upper=0.5)
# Likelihood: approximate via normal on decision time# Decision time = RT - Ter (must be positive)
decision_time = pm.Deterministic("decision_time", all_rts - Ter)
# Expected decision time under DDM (Ratcliff 1978, eq for unbiased start)# E[DT | correct] ≈ a/(2v) × tanh(av/2) for large a
expected_dt = pm.Deterministic(
"expected_dt",
(a / (2 * v)) * pm.math.tanh(a * v / 2)
)
# Variance of decision time (approximation)
var_dt = pm.Deterministic(
"var_dt",
(a / (2 * v ** 3)) * (1 - pm.math.tanh(a * v / 2) ** 2)
)
# Normal approximation likelihood
obs = pm.Normal(
"obs",
mu=expected_dt,
sigma=pm.math.sqrt(var_dt),
observed=decision_time
)
# Sample
idata = pm.sample(
draws=draws,
tune=tune,
target_accept=target_accept,
random_seed=seed,
progressbar=True,
)
return idata
defcompare_ddm_conditions(
idata_cond1,
idata_cond2,
param_names: List[str] = ["v", "a", "Ter"],
labels: Tuple[str, str] = ("Condition 1", "Condition 2"),
) -> pd.DataFrame:
"""
Compare posterior DDM parameters across two conditions.
Args:
idata_cond1: InferenceData from condition 1.
idata_cond2: InferenceData from condition 2.
param_names: DDM parameters to compare.
labels: Labels for the two conditions.
Returns:
DataFrame with posterior mean, HDI, and P(cond1 > cond2).
"""import arviz as az
rows = []
for param in param_names:
post1 = idata_cond1.posterior[param].values.flatten()
post2 = idata_cond2.posterior[param].values.flatten()
diff = post1 - post2
p_greater = (diff > 0).mean()
rows.append({
"parameter": param,
f"mean_{labels[0]}": round(post1.mean(), 4),
f"mean_{labels[1]}": round(post2.mean(), 4),
"mean_difference": round(diff.mean(), 4),
"HDI_2.5%": round(np.percentile(diff, 2.5), 4),
"HDI_97.5%": round(np.percentile(diff, 97.5), 4),
f"P({labels[0]}>{labels[1]})": round(p_greater, 3),
})
comparison_df = pd.DataFrame(rows).set_index("parameter")
print(f"\nDDM parameter comparison ({labels[0]} vs {labels[1]}):")
print(comparison_df)
return comparison_df
Delta Plot Analysis
defdelta_plot(
rt_dict: Dict[str, np.ndarray],
n_quantiles: int = 9,
output_path: Optional[str] = None,
) -> plt.Figure:
"""
Generate a delta plot showing condition differences across the RT distribution.
A delta plot plots quantile differences (Condition A - Condition B) against
the mean quantile RT. Positive slopes indicate the effect grows with RT
(common in congruency effects); negative slopes suggest strategic slowing.
Args:
rt_dict: Dict with exactly two conditions (first - second = delta).
n_quantiles: Number of quantile points (default 9 = deciles 0.1–0.9).
output_path: Optional path to save figure.
Returns:
Matplotlib Figure.
"""assertlen(rt_dict) == 2, "Delta plot requires exactly 2 conditions."
labels = list(rt_dict.keys())
rts_a, rts_b = rt_dict[labels[0]], rt_dict[labels[1]]
quantile_ps = np.linspace(0.1, 0.9, n_quantiles)
q_a = np.quantile(rts_a, quantile_ps)
q_b = np.quantile(rts_b, quantile_ps)
mean_q = (q_a + q_b) / 2
delta = q_a - q_b
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# Quantile functions
ax1 = axes[0]
ax1.plot(quantile_ps * 100, q_a, marker="o", label=labels[0], color="steelblue")
ax1.plot(quantile_ps * 100, q_b, marker="s", label=labels[1], color="tomato")
ax1.set_xlabel("Percentile")
ax1.set_ylabel("RT (ms)")
ax1.set_title("Conditional Accuracy Function")
ax1.legend()
ax1.grid(alpha=0.3)
# Delta plot
ax2 = axes[1]
ax2.plot(mean_q, delta, marker="o", color="purple", linewidth=2)
ax2.axhline(0, color="gray", linestyle="--", linewidth=1)
ax2.set_xlabel("Mean RT (ms)")
ax2.set_ylabel(f"Δ RT ({labels[0]} − {labels[1]}) (ms)")
ax2.set_title("Delta Plot")
ax2.grid(alpha=0.3)
fig.tight_layout()
if output_path:
fig.savefig(output_path, dpi=150)
plt.show()
return fig
Troubleshooting
Problem
Likely Cause
Solution
exponnorm.fit returns very large K
Outlier RTs inflate the tail
Apply absolute cutoffs before fitting
EZdiff produces negative Ter
Mean RT too close to a/(2v)
Check accuracy is well above 0.5; use edge correction
EZdiff v = 0
Accuracy exactly 0.5
Edge-correct accuracy to 0.501
PyMC divergences during sampling
Model mis-specification or bad priors
Increase target_accept to 0.95; reparameterize
PyMC Ter samples near 0
Prior too wide
Tighten Uniform(0.1, 0.4) based on task
K-S test rejects ex-Gaussian
Bimodal RT distribution (e.g., two strategies)
Separate trials by strategy or use mixture model
ndtri returns inf
Accuracy = 1.0 exactly
Clamp accuracy to 0.999 before calling EZdiff
External Resources
Wagenmakers, E.-J., et al. (2007). An EZ-diffusion model.
Psychonomic Bulletin & Review, 14(1), 3–22.
https://doi.org/10.3758/BF03194023
Ratcliff, R. (1978). A theory of memory retrieval. Psychological Review, 85(2).
Matzke, D., & Wagenmakers, E.-J. (2009). Psychological interpretation of ex-Gaussian.
Psychonomic Bulletin & Review, 16(5), 798–817.
