COGMOD-HWI/HW03/hw3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## HW3 Problem 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArMAAAIjCAYAAAAQgZNYAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAxkFJREFUeJzs3XlcVOX3wPHPMOwgOwqICiLggqC5peauuaSZZlrulUuZlUuLZlnaamVqX5dWNZfUsrTcc993xQ0FRNwQxIV9Z+b+/tCZnygqAwPDct6v17yKO/c+99wZlTPPnHselaIoCkIIIYQQQpRBZqYOQAghhBBCiMKSZFYIIYQQQpRZkswKIYQQQogyS5JZIYQQQghRZkkyK4QQQgghyixJZoUQQgghRJklyawQQgghhCizJJkVQgghhBBlliSzQgghhBCizJJkVogKJDU1lWHDhuHh4YFKpWLMmDEAXL9+nT59+uDq6opKpWLmzJkmjdMQD7um4nbx4kVUKhULFy7Ms33jxo00aNAAa2trVCoViYmJJRKPEEJUVJLMClHGLVy4EJVK9dDHgQMH9Pt+8cUXLFy4kNdff53FixczaNAgAMaOHcumTZuYOHEiixcvpkuXLkaP84svvmD16tXFMm5+15QfHx8f/etiZmaGk5MT9evXZ8SIERw8eLDIsdy6dYu+fftiY2PDnDlzWLx4MXZ2dgZduy5JvjdOFxcXunbtyv79+wsd29y5cx9IvMsS3XvXsWPHfJ//+eef9a/ZkSNHSji6grv/76u5uTlVq1Zl6NChxMTEADB06NBH/p3WPYYOHWraixGilFApiqKYOgghROEtXLiQl19+malTp+Lr6/vA8126dMHNzQ2AJ598EnNzc/bs2ZNnHw8PDzp27MiSJUuKLU57e3v69Olj9ITqYdeUHx8fH5ydnRk/fjwAKSkpnD17lj///JO4uDjGjh3Ld999V6DzKopCVlYWFhYWqNVq4M6sbNeuXdm8eXOepMuQa7948SK+vr689NJLdOvWDY1GQ0REBHPnziUjI4PDhw9Tv379AsV4r6CgINzc3NixY4fBx5YGPj4+XL9+nezsbGJiYvDw8MjzfNu2bTl48CCZmZkcPnyYxo0bmyjSR7v/72tmZiYHDhxg4cKF+Pj4cPr0aY4fP05UVJT+mOjoaCZPnsyIESNo1aqVfrufnx/Nmzc3xWUIUaqYmzoAIYRxdO3a9bG/wOPj46lbt26+252cnIopsuL1sGt6mKpVqzJw4MA826ZNm0b//v2ZMWMG/v7+vP766w89Pjc3F61Wi6WlJdbW1g/EAhjltXziiSfyxNmqVSu6du3KvHnzmDt3bpHHL4tatmzJ4cOHWbFiBW+//bZ++9WrV9m9eze9evXir7/+MmGEBXfv39dhw4bh5ubGtGnT+Pfff+nbt2+eJPXIkSNMnjyZ5s2bP/BnVwghZQZCVAg7duxApVIRHR3NunXr9F9T6r7yVBSFOXPm6LfrJCYmMmbMGKpVq4aVlRW1atVi2rRpaLXaPONrtVpmzZpF/fr1sba2xt3dnS5duui/7lWpVKSlpfHbb78V+CvS+Ph4Xn31VapUqYK1tTUhISH89ttvj72mixcvGvz62NjYsHjxYlxcXPj888/RfWGl+8r/22+/ZebMmfj5+WFlZUVYWNgDNbNt27ZlyJAhADRp0kR/jYW59vzoZuTunbEDWLBgAe3bt6dy5cpYWVlRt25d5s2bl2cfHx8fzpw5w86dO/UxtG3bVv98Qd/nh5k7dy716tXDysoKLy8v3njjjQdqhdu2bUtQUBBhYWG0a9cOW1tbqlatytdff13g18Da2prevXvz+++/59m+bNkynJ2d6dy5c77HnTt3jj59+uDi4oK1tTWNGzfm33//zbPP7du3eeedd6hfvz729vY4ODjQtWtXTpw4kWc/3Z+7P/74g88//xxvb2+sra3p0KED58+fL/C13O9h7+/DxMXF8fLLL+Pt7Y2VlRWenp707NmzUH/+hSjrZGZWiHIiKSmJmzdv5tmmUqlwdXWlTp06LF68mLFjx+Lt7a3/mr1hw4b6OtNOnToxePBg/bHp6em0adOGmJgYRo4cSfXq1dm3bx8TJ04kNjY2z01ir776KgsXLqRr164MGzaM3Nxcdu/ezYEDB2jcuDGLFy9m2LBhNG3alBEjRgB3viJ9mIyMDNq2bcv58+cZPXo0vr6+/PnnnwwdOpTExETefvvth16Tu7t7oV4/e3t7evXqxa+//kpYWBj16tXTP7dgwQIyMzMZMWIEVlZWuLi4PJDoTZo0icDAQH766Sf9V8h+fn507NjRoGt/GF2S4uzsnGf7vHnzqFevHs8++yzm5uasWbOGUaNGodVqeeONNwCYOXMmb775Jvb29kyaNAmAKlWqAIa9z/n55JNPmDJlCh07duT1118nPDycefPmcfjwYfbu3YuFhYV+34SEBLp06ULv3r3p27cvK1eu5P3336d+/fp07dq1QK9D//79efrpp4mKitK/jr///jt9+vTJcy6dM2fO0LJlS6pWrcqECROws7Pjjz/+4LnnnuOvv/6iV69eAFy4cIHVq1fzwgsv4Ovry/Xr1/nxxx9p06YNYWFheHl55Rn3q6++wszMjHfeeYekpCS+/vprBgwYUOja64e9vw/z/PPPc+bMGd588018fHyIj49n8+bNXL58GR8fn0LFIESZpQghyrQFCxYoQL4PKyurPPvWqFFDeeaZZx4YA1DeeOONPNs+/fRTxc7OTomIiMizfcKECYparVYuX76sKIqibNu2TQGUt95664FxtVqt/v/t7OyUIUOGFOiaZs6cqQDKkiVL9Nuys7OV5s2bK/b29kpycvJjryk/j9t3xowZCqD8888/iqIoSnR0tAIoDg4OSnx8fJ59dc8tWLBAv033Xhw+fDjPvoZcu27cKVOmKDdu3FDi4uKU3bt3K02aNFEA5c8//8yzf3p6+gNjdO7cWalZs2aebfXq1VPatGnzwL4FfZ/zEx8fr1haWipPP/20otFo9Ntnz56tAMr8+fP129q0aaMAyqJFi/TbsrKyFA8PD+X5559/6Dl0dO9dbm6u4uHhoXz66aeKoihKWFiYAig7d+7M9/Xv0KGDUr9+fSUzM1O/TavVKi1atFD8/f312zIzM/Ncg6LceS+srKyUqVOn6rdt375dAZQ6deooWVlZ+u2zZs1SAOXUqVOPvA5djFu2bFFu3LihXLlyRVm5cqXi7u6uWFlZKVeuXHngmMOHD+f5s5aQkKAAyjfffPPY102IikDKDIQoJ+bMmcPmzZvzPDZs2FDo8f78809atWqFs7MzN2/e1D86duyIRqNh165dAPz111+oVCo+/vjjB8a4t2TBEOvXr8fDw4OXXnpJv83CwoK33nqL1NRUdu7cWbiLegx7e3vgzo1h93r++ecLPeNbWB9//DHu7u54eHjQqlUrzp49y/Tp0+nTp0+e/WxsbPT/r5udb9OmDRcuXCApKemx5yno+5yfLVu2kJ2dzZgxYzAz+/9fJ8OHD8fBwYF169bl2d/e3j5PzaelpSVNmzblwoULj41TR61W07dvX5YtWwbA0qVLqVatWp4bo3Ru377Ntm3b6Nu3LykpKfpru3XrFp07dyYyMlLfQcDKykp/DRqNhlu3bmFvb09gYCDHjh17YOyXX34ZS0tL/c+68xf0Wjp27Ii7uzvVqlWjT58+2NnZ8e+//+Lt7f3YY21sbLC0tGTHjh0kJCQU6HxClGdSZiBEOdG0aVOj3sEdGRnJyZMnH5rE6W52ioqKwsvLCxcXF6Od+9KlS/j7++dJkADq1Kmjf744pKamAlCpUqU82/PrElHcRowYwQsvvEBmZibbtm3j+++/R6PRPLDf3r17+fjjj9m/fz/p6el5nktKSsLR0fGR5yno+5wf3fsQGBiYZ7ulpSU1a9Z84H3y9vZ+4AOOs7MzJ0+efGSM9+vfvz/ff/89J06c4Pfff+fFF1/M94PT+fPnURSFjz76iI8++ijfseLj46lataq+7nvu3LlER0fnea1dXV0fOK569eoPXAdQ4OR
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## HW3 Problem 3\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "def simulate_diffusion(v, a, beta, tau, dt=1e-3, scale=1.0, max_time=10., rng=None):\n",
    "    \"\"\"\n",
    "    Simulates one realization of the diffusion process given\n",
    "    a set of parameters and a step size `dt`.\n",
    "\n",
    "    Parameters:\n",
    "    -----------\n",
    "    v        : float\n",
    "        The drift rate (rate of information uptake)\n",
    "    a        : float\n",
    "        The boundary separation (decision threshold).\n",
    "    beta     : float in [0, 1]\n",
    "        Relative starting point (prior option preferences)\n",
    "    tau      : float\n",
    "        Non-decision time (additive constant)\n",
    "    dt       : float, optional (default: 1e-3 = 0.001)\n",
    "        The step size for the Euler algorithm.\n",
    "    scale    : float, optional (default: 1.0)\n",
    "        The scale (sqrt(var)) of the Wiener process. Not considered\n",
    "        a parameter and typically fixed to either 1.0 or 0.1.\n",
    "    max_time : float, optional (default: .10)\n",
    "        The maximum number of seconds before forced termination.\n",
    "    rng      : np.random.Generator or None, optional (default: None)\n",
    "        A random number generator with locally set seed or None\n",
    "        If None provided, a new generator will be spawned within the function.\n",
    "        \n",
    "    Returns:\n",
    "    --------\n",
    "    (x, c) - a tuple of response time (y - float) and a \n",
    "        binary decision (c - int) \n",
    "    \"\"\"\n",
    "\n",
    "    # Inits (process starts at relative starting point)\n",
    "    y = beta * a\n",
    "    num_steps = tau\n",
    "    const = scale * np.sqrt(dt)\n",
    "    if rng is None:\n",
    "        rng = np.random.default_rng()\n",
    "\n",
    "    # Loop through process and check boundary conditions\n",
    "    while (y <= a and y >= 0) and num_steps <= max_time:\n",
    "        # Perform diffusion equation\n",
    "        z = rng.normal()\n",
    "        y += v * dt + const * z\n",
    "\n",
    "        # Increment step counter\n",
    "        num_steps += dt\n",
    "\n",
    "    if y >= a:\n",
    "        c = 1.\n",
    "    else:\n",
    "        c = 0.\n",
    "    return (round(num_steps, 4), c)\n",
    "\n",
    "def simulate_diffusion_n(num_sims, v, a, beta, tau, dt=1e-3, scale=1.0, max_time=10., rng=None):\n",
    "    \"\"\"Simulate multiple realizations of the diffusion process.\"\"\"\n",
    "    \n",
    "    # Inits\n",
    "    data = np.zeros((num_sims, 2))\n",
    "    if rng is None:\n",
    "        rng = np.random.default_rng()\n",
    "    \n",
    "    # Create data set\n",
    "    for n in range(num_sims):\n",
    "        data[n, :] = simulate_diffusion(v, a, beta, tau, dt, scale, max_time, rng)\n",
    "    return data\n",
    "\n",
    "def visualize_data(data, figsize=(10, 5)):\n",
    "    \"\"\"Helper function to visualize a simple response time data set.\"\"\"\n",
    "\n",
    "    f, axarr = plt.subplots(1, 2, figsize=figsize)\n",
    "    \n",
    "    # Histogram of response times\n",
    "    sns.histplot(\n",
    "        data[:, 0][data[:, 1] == 1], ax=axarr[0], color='#AA0000', alpha=0.8, lw=2, label=f'Response 1')\n",
    "    sns.histplot(\n",
    "        data[:, 0][data[:, 1] == 0], ax=axarr[0], color='#0000AA', alpha=0.8, lw=2, label=f'Response 0')\n",
    "\n",
    "    # Barplot of categorical responses\n",
    "    response, frequency = np.unique(data[:, 1], return_counts=True)\n",
    "    sns.barplot(x=response.astype(np.int32), y=frequency, ax=axarr[1], alpha=0.8, color='#00AA00')\n",
    "\n",
    "    # Labels and embelishments\n",
    "    axarr[0].set_xlabel('Response time (s)', fontsize=16)\n",
    "    axarr[0].legend(fontsize=16)\n",
    "    axarr[0].set_ylabel('Count', fontsize=16)\n",
    "    axarr[1].set_xlabel('Response', fontsize=16)\n",
    "    axarr[1].set_ylabel('Frequency', fontsize=16)\n",
    "    for ax in axarr:\n",
    "        sns.despine(ax=ax)\n",
    "        ax.grid(alpha=0.1, color='black')\n",
    "\n",
    "    f.suptitle('Data Summary', fontsize=18)\n",
    "\n",
    "    f.tight_layout()\n",
    "\n",
    "# Baseline parameters\n",
    "a = 2.0\n",
    "beta = 0.5\n",
    "tau = 0.5\n",
    "scale = 1.0\n",
    "max_time = 10.0\n",
    "dt = 1e-3\n",
    "num_sims = 2000\n",
    "\n",
    "# Vary drift rate\n",
    "v_values = np.linspace(0.5, 1.5, 25)\n",
    "mean_rt_upper = []\n",
    "mean_rt_lower = []\n",
    "\n",
    "for v in v_values:\n",
    "    data = simulate_diffusion_n(num_sims, v, a, beta, tau, dt, scale, max_time)\n",
    "    mean_rt_upper.append(data[data[:, 1] == 1, 0].mean())\n",
    "    mean_rt_lower.append(data[data[:, 1] == 0, 0].mean())\n",
    "\n",
    "# Plot results\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.plot(v_values, mean_rt_upper, label='Upper Boundary (Correct)', color='maroon')\n",
    "plt.plot(v_values, mean_rt_lower, label='Lower Boundary (Incorrect)', color='gray')\n",
    "plt.xlabel('Drift Rate (v)')\n",
    "plt.ylabel('Mean Response Time (s)')\n",
    "plt.legend()\n",
    "plt.title('Effect of Drift Rate on Mean RTs')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}