<jupyter_start><jupyter_text>Sheet 5: Decision Making I - <NAME> & <NAME><jupyter_code>import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('png', 'pdf')
plt.rcParams.update({'font.size': 12})
plt.rcParams["figure.figsize"] = (16,8)<jupyter_output><empty_output><jupyter_text>Extended drift-diffusion model  1. Use the Euler-Maruyama method to implement the stochastic differential equation $dX = µdt + σdW$, which is a simple form of the drift diffusion model. Start with $X = 0$. Choose $σ = 0.5$, $µ = 0.2$.  Run five simulations of 10 seconds length with $∆t = 0.01s$ and plot the results.<jupyter_code>def euler_maruyama(ffun, gfun, x0, t, **kwargs):
    time_points = len(t)  # size of the time range
    dt = np.max(t)/time_points  # step size
    x = np.zeros((time_points, len(x0)))  # empty output array
    x[0, :] = x0
    for i in range(time_points-1):  # perform the integration at each time point
        W = np.random.normal(loc=0.0, scale=1.0, size=len(x0))  # sample from a normal dist.
        x[i+1, :] = x[i, :] + (np.array(ffun(x[i, :], t, **kwargs)) * dt) + (np.array(gfun(x[i, :], t, **kwargs)) * np.sqrt(dt) * W) 
    return x
np.random.seed(111)
parameters = {'sigma': 0.5, 'mu': 0.2}
T = 10  # [s]
dt = 0.01  # [s]
time_range = np.arange(0, T, dt)  # [s]
for i in range(5):  # run five simulations
    drift_diffusion = euler_maruyama(ffun=lambda x, t: parameters['mu'], gfun=lambda x, t: parameters['sigma'], x0=[0], t=time_range)
    plt.plot(time_range, drift_diffusion, label='%i' %(i+1))
plt.title(r'Simple Drift Diffusion - $\sigma = 0.5$, $\mu = 0.2$, $x_0 = 0$')
plt.xlabel('time [s]')
plt.legend()
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>Above are five simulations of the simple drift diffusion equation solved using Euler-Maruyama integration.   The only change between simulations are that a new random variable is sampled from the distribution.  2. Plot the percentage of ”up” ($X=+a$) responses as a function of the drift $µ$ by simulating 200 trials for $a = 1$ and the following $µ$ values: −0.1, 0, 0.1, 0.2, 0.5. Resimulate trials in which no decision boundary is reached during the simulation. Plot the results. What happens to this curve if you double, what if you halve the value of $a$?<jupyter_code>def percentage_up_plot(mu_values, trial_num, a, sim_type, color, marker):
    total_up_count = np.zeros(len(mu_values))  # total that reached the boundary
    total_down_count = np.zeros(len(mu_values))  # total that reached the boundary
    resim_index = np.zeros(len(mu_values))  # trials to be resimulated
    for i in range(len(mu_values)):  # loop over mu values
        up_count = np.zeros(int(trial_num[i]))  # count for the mu value
        down_count = np.zeros(int(trial_num[i]))  # count for the mu value
        for j in range(int(trial_num[i])):  # loop over trials
            drift_diffusion = euler_maruyama(ffun=lambda x, t: mu_values[i], gfun=lambda x, t: parameters['sigma'], x0=[0], t=time_range)
            if np.max(drift_diffusion) >= a:  # "up" boundary
                up_count[j] = 1
            elif np.max(drift_diffusion) <= -a:  # "down" boundary
                down_count[j] = 1
            else:  # neither
                resim_index[i] += 1
        total_up_count[i] = np.sum(up_count)
        total_down_count[i] = np.sum(down_count)
        plt.scatter(mu_values[i], total_up_count[i]/trial_num[i], color=color, label=r'%s $a=%%.1f$' %sim_type %a, marker=marker, s=2**7)
    plt.title(r'Simulation')
    plt.ylabel(r'percentage of "up" responses')
    plt.xlabel(r'$\mu$')
    handles, labels = plt.gca().get_legend_handles_labels()
    by_label = dict(zip(labels, handles))
    plt.legend(by_label.values(), by_label.keys())
    return total_up_count, total_down_count, resim_index, up_count, down_count
mu_values = [-0.1, 0, 0.1, 0.2, 0.5]
trial_num = np.full(len(mu_values), 200)  # number of trials for each mu
initial_sim = percentage_up_plot(mu_values, trial_num, a=1, sim_type='Initial', color='k', marker='o')
print(initial_sim[0], ', number that reached the boundary')
print(initial_sim[2], ', number that had to be resimulated')
second_sim = percentage_up_plot(mu_values, trial_num=initial_sim[2], a=1, sim_type='Second', color='b', marker='o')
plt.ylim(-0.1, 1.1)
plt.grid()
plt.show()<jupyter_output>[ 55.  98. 144. 175. 199.] , number that reached the boundary
[145. 102.  56.  25.   1.] , number that had to be resimulated<jupyter_text>These y axis on these plots should be scaled to (0, 100), code would break when tried to fix. We realized after the fact that the above code is note the best implementation of the idea, but it gets the point across.  As the variable $\mu$ is increased, the probability of the results being in the "up" state approach a value of 1.  For each resimulation of the states that did not reach a boundary, we achieved lower percentages of "up" state results.<jupyter_code>a_values = [0.5, 1.0, 2.0]
markers = ['o', 'x', '^']
for i in range(len(a_values)):
    initial_sim = percentage_up_plot(mu_values, trial_num, a=a_values[i], sim_type='Initial', color='k', marker=markers[i])
    second_sim = percentage_up_plot(mu_values, trial_num=initial_sim[2], a=a_values[i], sim_type='Second', color='b', marker=markers[i])
    plt.ylim(-0.1, 1.1)
    plt.grid()
plt.show()<jupyter_output>C:\Users\lcube\AppData\Roaming\Python\Python37\site-packages\ipykernel_launcher.py:18: RuntimeWarning: invalid value encountered in double_scalars<jupyter_text>Adding different values of $a$ and doing the same thing as before, we see that as $a$ is increased, the percentage of "up" states goes down.  Therefore, $a$ and percentage of "up" states are inversely proportional.  Although, for a sufficiently high enough value of $\mu$, they both approach a value of 100%.  3. Use the results from the last exercise to pick a value of $µ$ (with $σ = 0.5$) for which about 70% of the trials result in an ”up” response. Run 2000 trials and record the reaction times, i.e. the moment when the decision variable $X$ crosses one of the two decision thresholds at $X = ±a$ (for $a = 1$), as well as the choice, i.e. which of the boundaries was hit. Got kinda confused by what I was doing above, and didn't make it very reaction time friendly, thus new function: For a value of about $\mu = 0.1$, about 70% of the trials result in an "up" response.<jupyter_code>def reaction_times(mu, a, trial_num):
    parameters = {'sigma': 0.5, 'mu': mu}
    T = 10  # [s]
    dt = 0.01  # [s]
    time_range = np.arange(0, T, dt)  # [s]
    up_reaction = np.zeros(trial_num) 
    down_reaction = np.zeros(trial_num)
    for i in range(trial_num):  # loop over trials
        drift_diffusion = euler_maruyama(ffun=lambda x, t: parameters['mu'], gfun=lambda x, t: parameters['sigma'], x0=[0], t=time_range)
        for j in range(len(drift_diffusion)):  # loop over the time range
            if drift_diffusion[j] >= a:  # reaches "up" boundary
                up_reaction[i] = j
                break
            elif drift_diffusion[j] <= -a:  # reaches "down" boundary
                down_reaction[i] = j
                break
            else:
                pass
        # otherwise it did not reach a boundary
        if up_reaction[i] == 0:
            up_reaction[i] = np.nan
        if down_reaction[i] == 0:
            down_reaction[i] = np.nan
    return up_reaction, down_reaction
rtime = reaction_times(mu=0.1, a=1, trial_num=2000)
print(np.shape(rtime), 'shape of reaction times array')  # boundary, trials<jupyter_output>(2, 2000) shape of reaction times array<jupyter_text>4. Plot a histogram of the reaction times for each of the choices separately, and compare the distributions.  Calculate their individual mean reaction times and discuss whether (and why) this is what you would have expected.<jupyter_code>dt = 0.01
bins = 100

up_reaction = rtime[0]
down_reaction = rtime[1]
up_reaction = up_reaction[~np.isnan(up_reaction)]
down_reaction = down_reaction[~np.isnan(down_reaction)]

up_mean = np.mean(up_reaction) * dt
down_mean = np.mean(down_reaction) * dt
print('mean "up" reaction time: %.2f s' %up_mean)
print('mean "down" reaction time: %.2f s' %down_mean)

fig, axs = plt.subplots(2, 1, sharex=True, sharey=True)
axs[0].hist(up_reaction*dt, bins=bins, histtype='stepfilled', color='k', label='up')
axs[1].hist(down_reaction*dt, bins=bins, histtype='stepfilled', color='b', label='down')
axs[0].axvline(x=up_mean, color='r'), axs[1].axvline(x=down_mean, color='r')
axs[0].legend(), axs[1].legend()
axs[0].set_ylabel('frequency'), axs[1].set_ylabel('frequency')
axs[0].set_title('Reaction Times')
axs[1].set_xlabel('time [s]')
plt.show()<jupyter_output>mean "up" reaction time: 3.38 s
mean "down" reaction time: 3.16 s<jupyter_text>Above are the reaction time histograms for the up and down boundary conditions, where the red line represents the average reaction time for that condition.  The reaction times are roughly the same, which is what one would expect from this model as it is taken randomly with a start of $X=0$, and thus unbiased.  To make these match behavioral experiments, we will change this next:  5. In behavioral experiments error trials typically show longer reaction times than correct trials.  To reproduce these results, implement an extended version of the drift-diffusion model, in which the drift $µ$ for each trial is drawn from a Gaussian distribution with mean $m_µ$ and standard deviation $s_µ$. Find values for $m_µ$ and $s_µ$ that keep 70% ”up” rate but have reaction time distributions that (qualitatively) fit the experimental results better.  Plot your resulting distributions and mean reaction times for ”up” and ”down” responses from 2000 trials. Why does the extended verion of the drift-diffusion model produce this assymetry in the reaction time distributions?<jupyter_code>parameters = {'m_mu': 0.3, 's_mu':0.5}
T = 10  # [s]
dt = 0.01  # [s]
time_range = np.arange(0, T, dt)  # [s]
for i in range(5):  # run five simulations
    mu = np.random.normal(parameters['m_mu'], parameters['s_mu'])
    extended_drift = euler_maruyama(ffun=lambda x, t: mu, gfun=lambda x, t: 0.5, x0=[0], t=time_range)
    plt.plot(time_range, extended_drift, label='%i' %(i+1))
plt.title(r'Extended Drift Diffusion - $\sigma = 0.5$, $x_0 = 0$, $m_{\mu} = %.1f$, $s_{\mu} = %%.1f$' %parameters['m_mu'] %parameters['s_mu'])
plt.xlabel('time [s]')
plt.legend()
plt.grid()
plt.show()
def reaction_times_extended(m_mu, s_mu, a, trial_num):
    sigma = 0.5
    T = 10  # [s]
    dt = 0.01  # [s]
    time_range = np.arange(0, T, dt)  # [s]
    up_reaction = np.zeros(trial_num) 
    down_reaction = np.zeros(trial_num)
    for i in range(trial_num):  # loop over trials
        mu = np.random.normal(m_mu, s_mu)
        extended_drift = euler_maruyama(ffun=lambda x, t: mu, gfun=lambda x, t: sigma, x0=[0], t=time_range)
        for j in range(len(extended_drift)):  # loop over the time range
            if extended_drift[j] >= a:  # reaches "up" boundary
                up_reaction[i] = j
                break
            elif extended_drift[j] <= -a:  # reaches "down" boundary
                down_reaction[i] = j
                break
            else:
                pass
        # otherwise it did not reach a boundary
        if up_reaction[i] == 0:
            up_reaction[i] = np.nan
        if down_reaction[i] == 0:
            down_reaction[i] = np.nan
    return up_reaction, down_reaction
dt = 0.01
bins = 100

rtime = reaction_times_extended(m_mu=parameters['m_mu'], s_mu=parameters['s_mu'], a=1, trial_num=2000)
up_reaction = rtime[0]
down_reaction = rtime[1]
up_reaction = up_reaction[~np.isnan(up_reaction)]
down_reaction = down_reaction[~np.isnan(down_reaction)]

up_mean = np.mean(up_reaction) * dt
down_mean = np.mean(down_reaction) * dt
print('mean "up" reaction time: %.2f s' %up_mean)
print('mean "down" reaction time: %.2f s' %down_mean)
up_percentage = 100 * np.shape(up_reaction)[0] /(np.shape(up_reaction)[0] + np.shape(down_reaction)[0])
print(r'%.3f%% , percentage of reactions that are "up"' %up_percentage)
fig, axs = plt.subplots(2, 1, sharex=True, sharey=True)
axs[0].hist(up_reaction*dt, bins=bins, histtype='stepfilled', color='k', label='up')
axs[1].hist(down_reaction*dt, bins=bins, histtype='stepfilled', color='b', label='down')
axs[0].axvline(x=up_mean, color='r', label='mean'), axs[1].axvline(x=down_mean, color='r', label='mean')
axs[0].legend(), axs[1].legend()
axs[0].set_ylabel('frequency'), axs[1].set_ylabel('frequency')
axs[0].set_title(r'Reaction Times - $m_{\mu} = %.2f$, $s_{\mu} = %%.2f$' %parameters['m_mu'] %parameters['s_mu'])
axs[1].set_title(r'percentage "up" = %.3f %%' %up_percentage)
axs[1].set_xlabel('time [s]')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13540614069179946, 'fraction_numerical': 0.022386319471434123, 'mean_word_length': 4.049450549450549, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 30, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '39622', 'license': 'MIT', 'n_tokens_mistral': 4709, 'n_tokens_neox': 4419, 'n_words': 1356, 'path': 'MHBF/.ipynb_checkpoints/MHBF_CC5-checkpoint.ipynb', 'repo_name': 'lcubelongren/BCCN_classwork', 'size': '12882'}