import numpy as np
import matplotlib.pyplot as plt
from statsmodels.graphics import gofplots

LATER model

def later_model(theta, mu, sigma, rt_max, N):
    r = np.random.normal(mu, sigma, size=N)
    T = theta / r
    T[(T < 0) | (T > rt_max)] = rt_max
    return T

rt = later_model(theta=0.7, mu=1, sigma=0.25, rt_max=1.5, N=10000)

def plot_distributions(rt, bins=50):
    fig, axes = plt.subplots(1,2,figsize=(10,4))
    axes[0].hist(rt, bins=bins);
    axes[0].set_xlabel("Time (s)");
    axes[0].set_ylabel("Count");
    axes[0].set_title("Decision times");
    axes[1].hist(1./rt, bins=bins)
    axes[1].set_xlabel("Rate (1/s)");
    axes[1].set_ylabel("Count");
    axes[1].set_title("Inverse decision times");

plot_distributions(rt)

gofplots.qqplot(1./rt, ms=1);

Drift-diffusion model

def drift_diffusion_model(theta, mu, sigma, N):
    t = np.zeros(N)
    y = np.zeros(N)
    for n in range(N):
        while y[n] < theta:
            t[n] += 1
            y[n] += np.random.normal(mu, sigma)
    return t

rt = drift_diffusion_model(theta=2, mu=0.05, sigma=0.25, N=10000)

plot_distributions(rt)

gofplots.qqplot(1./rt, ms=1);

Race model

def race_model(theta, mu1, mu2, sigma1, sigma2, N):
    t = np.zeros(N)
    y1 = np.zeros(N)
    y2 = np.zeros(N)
    code = np.zeros(N)
    for n in range(N):
        while True:
            if y1[n] > theta:
                code[n] = 1
                break
            if y2[n] > theta:
                code[n] = 2
                break
            t[n] += 1
            y1[n] += np.random.normal(mu1, sigma1)
            y2[n] += np.random.normal(mu2, sigma2) 
    return t, code

rt, codes = race_model(theta=2, mu1=0.05, mu2=0.06, sigma1=0.25,sigma2=0.25, N=1000)

plot_distributions(rt)

gofplots.qqplot(1./rt, ms=1);

fig, ax = plt.subplots(1,1,figsize=(6,4));
ax.bar(range(2),np.bincount(codes.astype(int))[1:]/1000.);
ax.set_xticks([0,1]);
ax.set_xticklabels([1,2]);
ax.set_xlabel("Code");
ax.set_ylabel("Percentage");
ax.set_title("Outcome");