PDM4.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Perceptual decision making (4)\n",
    "\n",
    "In this section, we introduce how the standard drift-diffusion model(DDM) can be used to fit real data."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the standard DDM, sensory evidence in each frame follows a Gaussian distribution\n",
    "$$\n",
    "e = N(u, \\sigma^2)\n",
    "$$\n",
    "The visibility of stimuli becomes higher as coherence increases. In theory, a higher coherence provides stronger sensory evidence in each frame. Assuming $\\sigma=0.1$, we further assume a linear relationship between the mean $u$ of the distribution and cohernece level $coh$\n",
    "$$\n",
    "u = k*coh\n",
    "$$\n",
    "\n",
    "Let's simulate a drift-diffusion process. We have two settings:\n",
    "* The answer is correct if the accumulated evidence reaches the upper bound, and vice verse for the lower bound\n",
    "* Sensory evidence starts at the middle point between the upper and lower bounds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 3.0)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAEYCAYAAADmugmLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA1wElEQVR4nO3de1hU1f4/8PeeYWC4CYrgBRXEC4oXFC+ZIiqiiJ7jPU94PGgeK+2ifo+XsjqJxyzSsuxkdrAMKUt/qZUoGmpapJaKiPcLgoqYKHhBcgaGmfX7gxglQDbIMBt4v55nHmDPnr0/Qzlv1tprryUJIQSIiIgUSGXtAoiIiMrDkCIiIsViSBERkWIxpIiISLEYUkREpFgMKSIiUqwaDan//ve/aNOmDRo0aIARI0bgypUrpfaJi4uDr68vXFxcMG7cOOTk5NRkiUREpCA1FlI//fQTZs6ciaCgILz33nvYtWsXXnzxxRL7ZGVlYcKECfD29sYHH3yA+Ph4zJs3r6ZKJCIihbGpqRP169cPFy5cQNOmTXH69Gmo1WpoNJoS+yQkJECv12P27NkICwvDl19+iS1bttRUiUREpDA11pJSq9Xw8fHBjh070LNnTzRo0ABRUVEl9snIyAAAuLu7m7/m5ORAp9PVVJlERKQgNT5womfPntixYwcaNGiAv/zlLzAYDKX2kSQJAFA8Y1Pxz8UiIyMhSVKpBxER1S01FlJpaWlYv349XF1dERoailGjRuH06dO4fPmyeR9PT08AwI0bNwAA2dnZcHNzg1arLXGsyMhICCFKPIiIqO6psWtSZ8+eRXh4OCIiIvCXv/wF69evh7e3NwBg165d6NOnD0JCQmBra4sVK1bg+vXrSExMRHh4eE2VSEREClNjLamwsDC8//772Lt3L5566in4+voiPj4en3/+OYYMGYLU1FR4enpi06ZNSE9PxwsvvIBhw4Zh6dKlNVUiEREpjFRXluqQJIndfkREdQxnnCAiIsViSBERkWIxpIiISLEYUkREpFgMKSIiUiyGFBERKRZDioiIFIshRUREisWQIiIixWJIERGRYjGkiIhIsRhSRESkWAwpIiJSLIYUEREpFkOKiIgUiyFFRESKxZAiIiLFYkgREZFiMaSIiEixGFJERKRYDCkiIlIshhQRESkWQ4qIiBSLIUVERIrFkCIiIsViSBERkWIxpIiISLEYUkREpFgMKSIiUqwaC6nly5fDy8sLjo6OCAwMREpKSql95s+fD0mSzA9XV9eaKo+IiBRIEkIIS58kISEBoaGhmDZtGoYNG4ZZs2ZBq9UiNTW1xH6BgYHQ6XSIioqCJEnQaDQYMGCArHNIkoQaeCtERFSDKh1Ser0eer2+Uq2cQ4cOYfv27Xj++efh5uaGMWPGYOvWrdDpdLCxsQEAFBQUoEGDBtBqtcjLy0Pbtm0RHR2NoKAgeW+EIUVEVOfI7u774osv4O/vDycnJ0ydOhUxMTGYM2cOTCZTha/t1asXXn/9dbi5uSExMRHx8fEYOnSoOaAA4MqVK/Dz88OUKVMQFxcHIQTGjx+Pu3fvljpeZGRkiW5BSZLkvg0iIqpFZLWkoqOjMWPGDNjZ2UGv12P06NFo1qwZPv74Y7z00kt48803ZZ1s27ZtmDBhAho0aIADBw7A29u73H1XrFiB2bNn48cff5TVmmJLioio7pHVklq2bBmaNWuGzMxM87Z///vfaNKkCWJjY2WdaMOGDRg9ejRatGiBffv2lQqo06dPY8mSJebrVIWFhQAAW1tbWccnIqK6R1ZIZWRkoFOnTmjYsKF5W9OmTdGhQwfk5ORU+PqUlBRERERAo9Fg0aJFSEtLw65du3Du3Dns2rULeXl5EEJg4cKFmDFjBjZu3IiVK1eiffv26NWrV9XfHRER1WqyQsrX1xeJiYnYvHkzACA3Nxdr1qxBYmIiOnbsWOHr33//fRQUFECn0yE8PBxDhgzBkCFD8NFHH2HIkCFITU2Fn58f1q1bh0uXLmHSpEnw9PTEli1boFarH+0dEhFRrSXrmtSOHTswcuRIGI1GCCHMAxUkScLmzZsxcuRIixdaEV6TIiKqe2QPQT99+jTefvttnDx5EgaDAX5+fpg5cyb69Olj6RplYUgREdU9lbpPSqfTwd7eHgCQmZkJT09PixVWWQwpIqK6R9Y1qZs3b2LYsGGYOXOmeVvPnj0REhIia+AEERFRVcgKqdmzZyMhIcE8BF2v18PJyQl79uzBvHnzLFogERHVX7K6+zw8PNC4cWOkpKRAo9EAKLqPyd/fHzk5Obh27ZrFC60Iu/uIiOoeWS2pvLw8NGzY0BxQAGBjYwNnZ2fk5uZarDgiIqrfbCreBejduzcSExPx/PPPY8iQISgsLMT333+PQ4cOITAw0NI1EhFRPSWru+/QoUMICQnB3bt3zfdICSHg7OyMnTt3onfv3hYvtCLs7iMiqntkD0G/fv06/vvf/+LUqVMAgE6dOuH5559HkyZNLFqgXAwpIqK6p0YWPawJDCkiorpH1jUpg8GA2NhYJCUlIT8/v0QYSJKETz/91GIFEhFR/SWrJfXUU08hNja2zJaKJEkwGo0WKa4y2JIiIqp7ZLWktm7dCrVajeeeew4dO3YssaIuERGRpchqSXl5eaFNmzb44YcfaqKmKmFLioio7pF1M+/ChQtx9OhR/PLLL5auh4iIyExWSyo4OBjJycnIzc2Fk5MTHBwcSqwp9eCy8tbClhQRUd0jK6RUqvIbXBw4QUREliJrBER6erql6yAiIiqlUjfzZmdnIykpCY6OjvD394ezs7Mla6sUtqSIiOoeWQMnCgoKMG3aNDRv3hzDhw/H8uXLER0djV69eiErK8vSNRIRUT0lK6ReeuklfPbZZwgICDC3VrKyspCUlIRZs2ZZtEAiIqq/ZHX3NW3aFB4eHjh27BhUKhVGjx6NTZs2oXPnzrh27ZoilpBndx8RUd0jqyWVm5sLV1fXEtskSYKDgwMMBoMl6iIiIpIXUkFBQdi3bx9mzpwJAEhLS0NERASSkpLQt29fixZIRET1l6zuvtTUVISGhpYaiu7u7o49e/bAz8/PYgXKxe4+IqK6R/YQ9IKCAnz55Zc4efIkDAYD/Pz8MHHiRDg5OVm6RlkYUkREdQ8XPSQiIsUqd8YJtVot6wCSJKGwsLDaCiIiotrh2+RMLPv+LK7e1qG5qz3mhfpidHfPaj1HuQMnhBAVPuzs7NCyZUtZJ1q+fDm8vLzg6OiIwMBApKSklNonLi4Ovr6+cHFxwbhx4xQxtJ2IiEr7NjkTCzYfQ+ZtHQSAzNs6LNh8HN8mV++E4+WGlMFgMD/Wr18PrVaLmJgYZGdn4+bNm/jkk08ghMA777xT4UkSEhIwZ84cDB06FLGxsbh48SLGjRtXYp+srCxMmDAB3t7e+OCDDxAfH4958+Y9+jskIqJqt3jrKegMphLbdAYjln1/tlrPI+uaVOvWreHj44Pdu3eX2B4cHIzLly8jNTX1oa8/dOgQtm/fjueffx5ubm4YM2YMtm7dCp1OZ17l9/PPP0dERATi4+MRFhaG0NBQJCUlITs7W94b4TUpIiKLO5F5B2/vOIPE82V/NksA0qNGVNv5ZM2Cnp2djby8PKSnp6N169YAioalp6SkyLqZt1evXujVqxcAIDExEfHx8Rg6dGiJZegzMjIAFA1rL/6ak5MDnU4He3t7WW+meI0rIiKyHuntyu3/sAaGrJt5//rXvyInJwcdO3aEv78/unbtCj8/P9y6dQsTJkyQXci2bdswbNgwNGrUCCtXrixzn+KgKS66rOCJjIyEJEklHkREVPfIakl9/PHH0Gg0WLduHY4fPw4A0Gg0eOGFF7B06VJZJ9qwYQMmTZoEHx8fbN++Hd7e3iWe9/QsGhFy48YNAEWtNzc3N2i12lLHioyMRGRkZIlt7O4jIqo+eoMRa/dfxMo9qcjVF43gHt2tOeYM9UXLRg4AamZ0X6Xuk8rNzUV6ejokSUKbNm3g6Ogo63UpKSno3bs31Go11qxZg8aNGwMAWrVqhcuXL6NPnz64c+cOfHx8EBwcjPDwcEyfPh3h4eH49NNP5b0RhhQR0SMzmgS+Tc7EuwlncfWOHgAQ2LYxXg7rgM6eLjVeT7khde7cOTg7O6NZs2Y4d+7cQw/Svn37hz7/1FNPISYmptT2WbNmYcWKFUhOTka3bt2wdetWzJ07F1evXkVISAhWr14NNzc3eW+EIUVEVGVCCPx47gaitp/BmWt3AQAdmzXAgrAOCGrvbrW6yg0plUqFMWPGYNOmTVCpVOVe91HKzbwMKSKiqjmReQdvbT+NfalF96Z6utpjztD2GN3NEyqVda/5l3tNqlWrVuaRdq1ateLgBCKiOibj5j28k3AW3x29CgBooLXBC8FtEfG4N7QaebMOWRrn7iMiqgceHOTQ1EWL9h5OOJB2EwVGE2xtVJjS1xvPDWwDVwdba5dagqyQGjBgAKZMmYLx48fD2dm5JuqqNIYUEVHZiqYwOg6dwVjqubHdPfGvoe3RoqGDFSqrmKyQKr4mZW9vjzFjxiAiIgIhISGK6gJkSBERla1f1A/IvK0rtd3D2Q4HXw2xQkXyybqZ94cffsC0adNgb2+PdevWYdiwYWjVqhVeeeUVnD592tI1EhHRI7haRkABwI27+TVcSeVV6ppUYWEhEhISsHHjRmzevBl3797l6D4iIgU7duU2Rq3ch7I+Hj1d7bHv5eCaL6oSZLWkimVkZCAlJQXHjh1Dbm4uhBCKWZmXiIhKOnL5Fv6++lcIAfx5JLm9Ro15ob7WKawSZE2L9N5772H9+vU4fPiweVtwcDCmTJlSaskNIiKyvsMXb2LKZ4eQl1+I4V2aYnAHDyzfed6iUxhZguyBEwDQtm1bREREYPLkybIXO6wp7O4jIiryS1oOpsYcwr0CI0b6N8fyCf6wUVeq40wxZLWkpk6diilTpiAwMNDS9RAR0SPYl5qNf649BL3BhLHdPbHsCX+orTxrxKOQPXCisLAQcXFxOHjwILy8vBAcHAxHR0fz7OXWxpYUEdV3e89ex7OfJyG/0IQJPVvgrbFda3VAATJbUleuXEFoaCjOnDkDABg1ahSysrLw3nvvYefOneYFDYmIyDp2n87CjC+OoMBowsTHWuGNUZ2tPu9edZDVSTlz5kycPXsWc+fONbdWmjRpgrt372Lu3LkWLZCIiB5ux4lrmP5FEgqMJkzp640lo+tGQAEyu/ucnZ3Ro0cP7N27FyqVCqNHj8bmzZsRGBiIo0ePIi8vryZqfSh29xFRfbTt2G+YtT4ZhSaBaYGt8eqIjoqaDehRyerus7W1xY0bN0qEgMFgwNWrV+HgoMz5noiI6rrvjmbi/zYchUkAMwa2wfxQ3zoVUIDM7r5x48bh9OnTCAgIgCRJOHjwIDp16oRLly5h1KhRlq6RiIj+ZGPSFXNAzRzcrk4GFCCzu+/evXt45pln8NVXX5VoTQ0bNgxfffUVXFxqfknhP2N3HxHVF+sPXsaCb45DCGDOkPZ4cXA7a5dkMZWau+/ixYs4efIkDAYD/Pz8Klw2viYxpIioPvj8l0v497cnAAAvDeuAGQPbWLkiy+Kih0REtcRn+9KxKO4UAOC1ER0xrb+PlSuyPFkDJ4iIyLpW/5SGJfFFSyP9Z1QnRDzubd2CaghDiohI4VbuScWy788CAN4c0wUTH2tl5YpqDkOKiEihhBBYsfs83t91HpIEvD2uKyb0VNbk3pZWbkj99NNPsg8SFBRULcUQEdU23yZnYtn3Z6t1CYyiY55B5m29edvyCf4Y073Fo5Zb65QbUgMHDpQ95t5oNFZbQUREtcW3yZlYsPk4dIaiz8DM2zos2HwcAPCXrs1QYDShoNCE/ML7X/MLjSgoLL29wFi0/WD6TWxJuQqD8f5AMI1agoS6dw+UHOWO7gsMDCwRUgcOHIBKpUKHDh2gUqlw6tQp2NnZYfz48fjss89qrODycHQfEdW0flE/IPO2rkbOVRuWereEcltSP//8s/n7hQsX4syZMzhw4ADatm0LADh79iz69u2ruMUPiYhqwvmsuw8NKJUE2NqoYGej/uOrCrY2KtiqVbDTqGGnVsFOU/Tzg8//v8NXyjze1RoKQ6WRdZ9Uo0aN0L17d+zevbvE9kGDBuH48ePIzs62WIFysSVFRDUhK1eP93aew/87nAFTOR85zV212P/y4Codv7zWWX1tScmau0+j0WDfvn344osvcOvWLdy8eRMxMTHYv38/NBqNpWskIrK6u3oD3vn+LAYs24P1hzIgSRL6tWkErU3Jj1F7jRrzQztU+TzzQn1hr1GXOua8UN8qH7M2kzUEfcaMGfjPf/6DyZMnl9guhMD8+fMtUhgRkRIUFJqw7tdL+O8Pqbj5ewEAIKxzU8wL9YWPu1O1j+4rfm11jxisrWRPi/TZZ59h5cqVuHDhAiRJQseOHTFnzhyMHTtW9snu3buHnj17wtbWFkePHi1zn/nz52PZsmXmn11cXHD79u0Kj83uPiKqTkIIbD32G5Z9fxaXb94DAPTyboiXwzqih1dDK1dXf9TY3H2HDx/G9OnTkZSUBH9//3JDKjAwEDqdDlFRUZAkCRqNBgMGDKjw+AwpIqouBy7kIGr7aaRcuQMAaOPuiJeGdcAQvyZ1cjkMJZM948T58+fx4Ycf4uDBg+jRowfGjh2LvLw8jBw5Utbre/XqhTFjxiAjI6PcfQoKCnD48GFotVqEhYWhbdu2iI6OllsiEdEjOXvtLqK2n8aeszcAAO7Odvi/kPaY0LMFbNSyLuFTNZMVUj///DNCQ0Oh0+kgSRKaNWuGHTt24N1338XatWsxadKkCo+RlJSEgIAAeHt7l7vPlStX4Ofnh6CgIISGhmL27NkYP348Lly4AGdnZ9lvioioMn67o8PyhHPYdOQKTAJwtFXj2QFtMK1/azjYcvY4a5L1p8HcuXOhUqmwbds2c5daWFgY7O3t8dZbb8k6UUBAQIX7+Pj44MiRI3j//fcRFhaG5557Djdu3EBycnKJ/SIjIyFJUokHEVFl3dEZ8PaOMxi4bC++TroClSRh8uNe+HH+IMwc3I4BpQCy/gukpKRgwIABCAsLM28bNGgQHnvsMezfv7/aijl9+jQ2b96Mv/3tb2jbti0KCwsBALa2tiX2i4yMRGRkZIltDCoieSw311ztOmYDexsUFJqgM5gAACO6NMPcUF+0buz4SOeg6iUrpNzd3XHs2DFcv37dvC01NRXJyclo1qzZIxWQlpaGtLQ09OnTB0IILFy4EHv37sWzzz6LlStXon379ujVq9cjnYOIijxsrrmKAsBgNEFnMEJvMCLfYILeYITOYMTOU1n4349pKDCazMecv/EYUq/fxeCOTe7PpqAumnnB9oGZFzTqsntC5NZpMBbVof+jnvzC+9+bvxYaoSsw4tf0m/juaKZ5Trw7uqI/gn0aO+LdCf7o3ooj9pRI1ui+d999F/PmzYONjQ2MRiM0Gg0KCwthMpmwZMkSLFiwQPYJvb294erqah7dFxkZiUWLFiE5ORndunXDhg0b8O9//xuXL19Gr1698Mknn8DXt+Kb2Di6j6hi5c1moLVRoad3I3MIlfiQNxihLzTBWN70Co9AkgCNWgU79f3wsrNR4cotHQrLOJ9akuDqoKnWmjxdtdhXxdkhyPJkD0H/5JNP8MYbb+Dy5csAgKZNm2L27NmKuZmXIUVUsdYvb0NV/5WoVRK0NirY26phZ6OGVlP0/YnM3HJf49/C5Y8Zvk3mmb8f/L6sIKoMlQRoNeqih40KWo0adpo/aiverlFBa6PG5uTMMo8hAUiPGvFIdZDlyL4qOG3aNEybNg2///47DAYDXF1dkZ+fj+zsbDRu3NiSNRLRIxJCYPORTEgSUNbfco0cbfH+37qZP9SLP+DtNEUf/PYaNTTlDMF+2Fxz370Q+NC6jCZxf9kK4/0lLJ6M/gXX7+aX2r9JAzvEvRj4Ryipy+0uLMuv6TfLrLO5q72s15N1yBrdp1arMW7cOACAo6MjXF1dAQBDhw7l9SIihbuY/Tsmffor5nydApMoan08yF6jxut/8UNQe3f0bt0IXVu4ol0TZ7Rs5AAPZy0aaDXlBhTwaHPNqVUS7G3VcHHQwMNZixYNHeDj7oRXhncs85gLwjqaa7K1UVVqwBTnxKudym1JxcbGYseOHQCK/go7dOgQJk6caH7eZDIhOTkZJpPJ8lUSUaUZjCasTkzDil3nkV9ogquDBq+N8INaAt5JOKfoueZqyzHJ8sq9JpWRkYGOHTvi3r175r9Wytp19OjR2Lx5s2WrlIHXpIjuS758Cws2H8eZa3cBAGO6e+K1ER3h5mRn5cqIKuehAyd+/PFHpKenY+rUqejevTtmzpxpfk6lUsHDwwODBg0qdR+TNTCkiIC8/EK88/1ZrD1wEUIALRvZY8noLghq727t0oiqRNbovh9//BGNGzdGp06daqKmKmFIUX2381QWXv/uBH67o4daJWFa/9aYPbg97G3VFb+YSKFkhVR+fj5WrVqF8+fPQ6/Xm8MgLy8PBw8exMWLFy1dZ4UYUlRfXc/VY+GWk9h+4hqAomHfb47tgk7NXaxcGdGjkzUE/cUXX8Snn34KIUSpMHB05BQiRNZgMgl8degyorafwV19IRxs1Zg71BeT+3pD/echfES1lKwh6N999x3c3d2xceNG2NjY4JNPPsGbb74JSZKwcOFCS9dIRH9yPusuJvzvAF795gTu6gsxuIMHdv5rAKYGtmZAUZ0iq7tPq9Vi0KBB2L59O/r374+nnnoKU6dOxcCBA5GZmYnz58/XRK0Pxe4+qg/0BiM+2nsBq/amwmAUaOxkh0UjO2F4l6acZJnqJFndfc2bN8ehQ4dw5swZ9O3bF5999hk6dOiA8+fP486dO5aukWoZS8xeXRtYetZuNydbSBJw424BACC8dyu8PKwDXBw01VE+kSLJaklFRUXhlVdeweLFixEcHIzAwKKpToQQGDBgAPbs2WPxQivClpQy/Hn2aqDorv63xnZ5pA9spQdfZd+3EAIGo0CB0YR8g/GPr6YHvhqx+/R1fJKYbp5dvJiHsy0+nNgDvVs3svj7IrI22RPMbtiwAW3atEHPnj0RHR2NDz74AK1bt8aKFSvg4+Nj6TorxJCqmqp8+JtMArl6A27fM+C2zoBb9wpw554Bt+8V4N2d53BXX1jqNU52Nng2yKdoctI/JgO1ty2af02rUcPeVvXHpKXqP7arzBOHxqVctUjwVafy5q+zUUlo1cgB+YWmPx5G8ySrVf3ftbmLFvsXcNZuqh9kh5TSMaQqr6y//jVqCX/t2gyt3ByLQuheAW7rDCW+v6MzVPkDtjp5utpj38vBVq3h9r0CxB+/hle+OV7p19qopPtrLdkUhfSDPydfvl3m6zhrN9Un5V6TenCevoeRJAnr1q2rtoLI8oQQOP3bXfz7uxMlAgoADEaBzclXKzyGs9YGrg4aNHSwhYu9Bq4OtnC11+Db5EzczS/dknLW2uAffbygNxQtnJf/x2J0eoMJuoL73+f/sZCeeT2jQmO5gZh5W4eYfekY0qkpPGtwJut7BYXYdfo6thzNxI/nbpgX0SuLh7Mdvny6D+zKCKOKRuGV1zrjrN1Un5TbklKpZI1OhyRJMBqNFe9oYWxJPZzBaMKh9JtIOJWFnaeyyvzwe9DM4LZw+SN4Gjpq4GJvaw6lBlob2JQzK3Z1X5MSQqBf1A+4ekf/0P06ezbAkI5NMbRTE3Ro6lztI90MRhMSz9/Ad0evYuepLNwrKHp/Kgno17YxWjS0xzfJmdAb7l8/etQuSUtd3yOqTcoNqbVr18o+yOTJk6utoKpiSJWWl1+IH8/ewM5T1/DDmevIfeBaUWMnO9wrKDR/2D7oUbvRqnuQQ1kf1lqNCuN6eOJmngE/nrtR4n20bGRvDqyeXg3LDdSKmEwChy/dwndHMxF//DfcumcwPxfQyhWjunlieJdmcHe2M9dpydF9ShwwQmRpvCZVx1zP1WPn6aLW0v7UnBIjw9p6OGGIXxMM8WuCbi1csaUWDEgo9rAPa73BiP0XspFwMgu7TmchO6/A/LqGDhoM7lj0noPauVc4j11xV+h3KZmIO3q1RAuunYcTRnf3xF+7NkcrNwfLvFEiKkFWSE2dOrX8A0gSPv3002otqirqQ0iV9UE9qltznL+eh52nspBwKgspGbfN+0sS0KNVQwzt1ARD/JqidePSU1jVtb/UjSaBoxm3kHAyC9+fvIaLOffMz2k1KvRv544hfk0wuIMHEs9nm9+7RwM79GjVEOev5+H89Tzzazxd7fFX/+YY1a25RboRiejhZIVUWdenikOB16RqRlldXmqVBFd7G+T8fr8bys6m6IN4qF8TBHf0QON6vH6QEAKp1/OQUEaAA0XXk0zlLKU+okszjOzWHD1aNYSK0wwRWY2skHqwpSSEQH5+Pn755Rds2bIFH3/8McLDwy1apBxKC6lHaaEYTQI3fy/A9bt63Libjxt387F466kS15Qe1MjRFsEdPDDUrwn6y+jSqq+u3bnfFfrTuRtl7uPmaItfXhn80OXSiajmPNI1qbFjx8JkMuHbb7+txpKqRkkhVd6orMi/+uExHzfcyMvH9dx83Lirv/99XlEYXb+bj5y8/DL/wi+LBCD1zeGcVLSSWr+8DWX9inkPEpGyyJq7ryw6nQ6ZmZk4ceJEddZTJyz7/myp+490BiNe2iz/hs9GjrZwd7KDu7MdPJztkHAqC3ll3H/U3NWeAVUFzV3teQ8SUS0ge4LZB5lMJty5cwcFBQXw9fW1SGG1UaHRVOE9SC0b2ZvDpyiAtEXfO9nBo0HRNjdHO9jalOxuKq91Ni+Uv/+qmBfqy98nUS0gK6SuXbtW5vYWLVrg448/rtaCaqPsvHysP3gZ6369jN8ectOpp6sWifOrdv9R8fWsujQSz5r4+ySqHWRdk7p06VKpbQ4ODnB3d7dIUVVhjWtSRzNuI3b/RWw99pv5fiSfxo7o1tIF8SeuVevsA0RE9RFv5q0kvcGIbcd+Q+yBi0i5cuePcwODOzTB5L5e6NemMVQqqc7df0REZA2yQuro0aN45plncPz4cRQUFJR4TpIkFBaWPTS6Jlk6pK7e1mHdr5ew/mAGcn4v+h24Omjwt14tMekxL7RsxBkIiIiqm6yQ6tatG44dO1bu8yaTqdznaoolQkoIgQNpOYjdfwkJp66Zh4X7NWuAKX29MbJbc2g1vCeJiMhSZIWUg4MDmjdvjl27dsHT07PUDBRqtbwP6nv37qFnz56wtbXF0aNHy9wnLi4Oc+fOxbVr1xASEoLo6Gi4ublV/EYeIaT+3DU3M7gtDCaB2AMXcS6raIocG5WE4V2aYXJfLwS0asjpcYiIaoCskAoJCUFOTg6Sk5OrfKLDhw9j+vTpSEpKgr+/f5khlZWVBW9vbwQFBWHixImYPn06wsPDsWbNmgqPX9WQKmto94Pcne3w98daYWLvVvBooK308YmIqOpkDUFfs2YN/P39ERgYiMDAQDg4lLz+8vrrr1d4jF69emHMmDHIyMgod5+EhATo9XrMnj0bYWFh+PLLL7FlyxY5JVZZWTfeAoCtWoV3J/gjtFPTUvcsERFRzZAVUu+99x7u3LmD/fv3Y//+/eauruIJZuWEVFJSEgICAuDt7V3uPsUBVjy03d3dHTk5OdDpdLC3vz8TQGRkJBYtWiSn9ApdLefGW4PRhL/6Ny/zOSIiqhmyQuqTTz6BjY0NJk6cWOY1KTkCAgJk7/tgCD74c7HIyEhERkaW+ZrK4vQ4RETKJSuk3Nzc4OXlhZiYGIsW4+lZdB/RjRtFM1RnZ2fDzc0NWq3lrgVxehwiIuWS1SR66623kJycjLi4OOTm5qKgoKDE41GkpaVh165dyMvLQ0hICGxtbbFixQrExsYiMTERo0aNeqTjV2R0d0+8NbYLPF3tIaFokTvODEFEpAyyRve1b98eV65cQX5+fukDVPJmXm9vb7i6uppH9xVfX0pOTka3bt2wdetWzJ07F1evXkVISAhWr15t8SHoRESkTFVemfdBdfVmXiIisi5Z16SUEEJERFT/yAqpy5cvP/T5Vq1aVUsxRERED5Ld3VfeEO/6MsEsERHVPNkr8z5471J+fj5u3rwJR0dH9OzZ06IFEhFR/SUrpK5cuVJq2/nz5zFkyBBMnDix2osiIiICHnHRw4iICCQmJiI9Pb06a6oSdvcREdU9slpS0dHRJX42mUy4evUqvvnmGy5ZQUREFiMrpKZPn14qjIpbLf/85z+rvyoiIiLIDKmIiIhSIeXg4ICAgABMnjzZIoURERE90jUpJeE1KSKiukf2mhtRUVFYsmSJ+ee+ffvijTfesEhRREREgMyQWrRoEV555RUkJCQAAPR6PQ4ePIiFCxdi8eLFFi2QiIjqL1ndfd7e3tDr9fjpp5/Qvn17AEB6ejr69u0Le3t7pKWlWbzQirC7j4io7pHVkrp27Rq6dOliDigAaN26NTp16oTffvvNYsUREVH9JiukfHx88PPPP+Prr7/G7du3kZ2djXXr1iExMRGtW7e2dI1ERFRPyeru++KLL8ochi6EQExMDCIiIixWoFzs7iMiqntkD0HfsWMHlixZglOnTgEAOnXqhAULFiAsLMyiBcrFkCIiqnt4nxQRESlWpe6TevPNN80/8z4pIiKytErdJ/X9998D4H1SRERUM3ifFBERKRbvkyIiIsXifVJERKRYvE+KiIgUi/dJERGRYj3yfVLnzp0rca3KWhhSRER1j+z7pB6Um5uL1atXo1+/fvDz86vumoiIiADIXD6+2M6dOxETE4Nvv/0Wer0eQohS16mIiIiqS4UtqfPnz+PVV19Fq1atMGzYMKxfvx46nQ729vaYM2cOzp8/L/tkcXFx8PX1hYuLC8aNG4ecnJxS+8yfPx+SJJkfrq6ulXpDRERUd5R7TWr16tWIiYnBL7/8AqBoJJ+joyOGDx+Or7/+Gp07d8axY8dknygrKwve3t4ICgrCxIkTMX36dISHh2PNmjUl9gsMDIROp0NUVBQkSYJGo8GAAQMqfiO8JkVEVOeUG1IqlQqSJMHOzg5jx47FuHHjEBYWBq1WC5VKVemQ+vzzzxEREYH4+HiEhYUhNDQUSUlJyM7ONu9TUFCABg0aQKvVIi8vD23btkV0dDSCgoIqfiMMKSKiOueh3X1CCOTn5yMlJQUpKSlITU2t8okyMjIAAO7u7uavOTk50Ol05n2uXLkCPz8/TJkyBXFxcRBCYPz48bh7926Vz0tERLVXuSF14sQJzJ8/H56enjh58iQWL14Mf39/tGvXDgCq3GopHmhR/PoHB174+PjgyJEjeP/99xEWFobnnnsON27cQHJycoljREZGlrhuxcEbRER1U7kh5efnh6ioKFy6dAm7du3CpEmT4OjoiAsXLgAATp06hT59+uDTTz+VdSJPT08AwI0bNwAA2dnZcHNzg1arNe9z+vRpLFmyxNxiKywsBADY2tqWOFZkZCSEECUeRERU91TqZl6dTodNmzbhiy++wO7du2E0GiFJEoxGY4WvzczMhI+PD4KDgxEeHm4eOPHqq68iLS0Nffr0weXLl9G1a1cMGjQIzz77LObPnw+NRoNTp05BrVY//I3wmhQRUZ1TqZt57e3tMWnSJOzYsQMZGRlYunQpunTpIuu1np6e2LRpE9LT0/HCCy9g2LBhWLp0KWJjYzFkyBCkpqbCz88P69atw6VLlzBp0iR4enpiy5YtFQYUERHVTVw+noiIFKtK0yIRERHVBIYUEREpFkOKiIgUiyFFRESKxZAiIiLFYkgREZFiMaSIiEixGFJERKRYDCkiIlIshhQRESkWQ4qIiBSLIUVERIrFkCIiIsViSBERkWIxpIiISLEYUkREpFgMKSIiUiyGFBERKRZDioiIFIshRUREisWQIiIixWJIERGRYjGkiIhIsRhSRESkWAwpIiJSLIYUEREpFkOKiIgUiyFFRESKxZAiIiLFqtGQiouLg6+vL1xcXDBu3Djk5ORUaR8iIqofaiyksrKyMGHCBHh7e+ODDz5AfHw85s2bV+l9iIio/qixkEpISIBer8fs2bMxefJkBAUFYcuWLZXeh4iI6g+bmjpRRkYGAMDd3d38NScnBzqdDvb29rL3AYDIyEgsWrSo1DkkSbLoeyAiIssQQpS5vcZCqlhxkBQXVFawVLRPZGQkIiMjLVJbeb8oJWGd1as21FkbagRYZ3VjnTXY3efp6QkAuHHjBgAgOzsbbm5u0Gq1ldqHiIjqjxprSYWEhMDW1hYrVqzA9evXkZiYiPDwcKSlpSEtLQ19+vQpdx8iIqqfarQltWnTJqSnp+OFF17AsGHDsHTpUsTGxmLIkCFITU0tdx8iIqqnBJktXLjQ2iXIwjqrV22oszbUKATrrG6sUwhJiFpwVY6IiOolTotERESKxZAiIiLFYkih9swXuHz5cnh5ecHR0RGBgYFISUmxdkkP9fHHH0OSJMTExFi7lDIdP34cAwYMgLOzM7p3744jR45Yu6QyrVq1Cl5eXnBwcEDfvn1x4sQJa5dkdu/ePfj5+aFbt27mbQaDATNnzkTjxo3RsmVLvP/++1arDyi7xuzsbIwfPx6urq5wd3fHs88+i/z8fOsVibLrLFZYWIj+/fsrYsKC8ur86KOP0LZtWzg7O+Mf//gHCgoKqueEFrvaVUtcu3ZNaLVaMXToUBETEyO0Wq146qmnrF1WKd9//70AIKZNmyY2btwoPD09RZs2baxdVrmOHj0qtFqtACA+++wza5dTyp07d4Sbm5vo06eP+Oqrr0Tr1q1F586drV1WKRcuXBCSJInw8HCxbt064e3tLQIDA61dlhBCiEOHDokePXoIAMLf39+8fcWKFQKAePvtt8XTTz8tAIgff/xRUTWGh4cLrVYrYmJixH/+8x8BQCxevNgqNT6szmIvvfSSACCs/ZFdXp3r1q0TAMTLL78s3n77bQFAvPPOO9VyznofUrGxsQKAiI+PF0IIMXToUOHm5mblqko7ePCgWLRokcjOzhZCCDF69GhhY2MjDAaDlSsrLTc3V7Rr10507dpVsSG1YcMGAUBs27ZNmEwmkZ6eLi5dumTtsko5d+6cACDeeustcfnyZdG7d28RHBxs7bKEEEIAEGPGjBEeHh4lPrCCg4OFh4eHEEKIq1evCgBizpw5iqpxzZo1Ijo6WgghxO3btwUAMWXKFKvUKET5dQohxPbt24VGoxF+fn5WD6ny6gwLCzP/NzeZTCIlJUXcvn27Ws5Z77v7HjZfoJL06tULr7/+Otzc3JCYmIj4+HgMHToUNjY1PrNVhZ599lnY2dkhKirK2qWU68KFCwCA1atXw8nJCcHBwYrqRivWrl07LF68GAsWLECrVq1w5swZq3efFUtKSsLmzZtLzKsJFP2bevDfU/E2ayivxqeeegpPP/00AGDBggUAgOHDh9d4fcXKqzMzMxP/+Mc/8Prrr6NXr15Wqu6+8uq8cOEC7Ozs0L9/f2i1Wrz00kvV1n1a70OqmJw5BZVg27ZtGDZsGBo1aoSVK1dau5xSoqOj8c0332Dt2rVQq9UAiq5RGI1GK1dWUvF/ZyEEvvvuO7i6umLChAm4ffu2dQv7k7Nnz+Ltt9/GpEmTsGPHDjRp0gR/+9vfUFhYaO3SEBAQUO5zf/73Y61/Tw+r0Wg04p///CdWrVqFJ598Ek888UQNVlZSWXUajUaEh4ejY8eO+Ne//mX+N6TX62u6PLPyfp9CCGRkZGDixIlYvXo1du3ahX/961/Vcs56H1K1ab7ADRs2YPTo0WjRogX27dsHb29va5dUypdffgm9Xo8ePXogNDQUAPDMM8/g888/t3JlJXl5eQEAIiIiEBISgieeeAK///67uYWlFHFxccjLy8OMGTMQGhqK0aNH4/Tp07h48aK1SyuXp6en+d9T8dcWLVpYs6RSCgsL8cQTT2DNmjV49tln8cUXX1i7pFIyMjKQmJiIxMREODo6mmv8cytGCby8vODs7IwZM2YgIiICPj4+1TawS3l9RTWstswXmJKSgoiICGg0GixatMg85+HAgQMV1eX30UcfITc3FwBw5MgRPP/883jttdcwYsQIK1dWUlhYGFxcXLB8+XKo1Wps3LgRrq6uaN++vbVLK6F4BNXSpUsRERGBr7/+Gh4eHuaQVaIRI0Zg7969eOedd5CamgoAGDlypJWrKum1117DN998g379+mH8+PHYs2cPPDw80LVrV2uXZtasWTMcOHDA/PPixYsRHx9fYptSPPnkk9i1axdee+01tG7dGmlpaZg8eXL1HLxarmzVcnFxccLX11c4OzuLMWPGmAcnKMmUKVPMo3sefNy6dcvapZVrz549ih04IUTRYJTHHntMODg4iO7du4u9e/dau6Qyffjhh6J169bC3t5e9O7dW/zyyy/WLqkELy+vEhfRCwoKxKxZs0Tjxo1F8+bNxXvvvWe12oo9WKNOpzOPPH3wMWrUKKvWKETp3+WDJk+ebPWBE8X+XKfJZBJLly4VLVu2FC4uLuLJJ58UN2/erJZzcVokIiJSrHp/TYqIiJSLIUVERIrFkCIiIsViSBERkWIxpIiISLEYUkQyrF69GpIkISUlBbt27YIkSVi7dm2Vj7du3Tq0bNkSWq0WY8eOLXOfOXPmwNXVFU5OTnj33XerfC6i2owhRSTD0aNHYWtrCz8/Pxw7dgwAylxSQa5Zs2bhypUr6Nq1Kzp37lzq+RMnTmD58uX4/fff8dhjjylydhGimqCcqQqIFCwlJQV+fn7QaDRISUmBRqOBn59flY5lMpnMa5bt37+/zBlDrl+/DgB47LHHsHv37qoXTlTLsSVF9BDe3t6QJAn79u3D0aNHIUkSYmNjYTAY0K5duzJfYzKZEBUVhbZt28LOzg6+vr5YtWoVAODixYvmiXcBQKPRlFoUMiYmBoMHDwYA7Nu3D5Ik4eLFixg4cCAkSUJUVBSaNGmCtm3bIj8/HwcOHEDfvn3h5OQEJycn9OvXD0lJSQCAvXv3QpIkTJ06FS+//DIaNmyIZs2aYe3atdizZw/8/Pxgb2+PESNG4ObNm+YaduzYAX9/f9jZ2aFt27b43//+V52/ViL5qmXeCqI66umnnxZBQUECgAgICBDDhw8XAES7du3E008/XeZrnn/+eQFAuLq6iuDgYOHo6CgAiDfeeENkZWWJkSNHlpiKZ/fu3SVev3v3btG3b18BQLi5uYlRo0aJrKwsMWDAAAFA2NjYiH79+om///3vIi8vTzRq1EioVCrRr18/0a1bNwFAdOnSRQhxf2oqe3t74eHhYV6wTqPRCEdHRxEUFCRcXV1LrPl06tQpYWNjI+zs7ERISIho1qyZACDWrl1r2V82URkYUkQVWL16tQAgDh8+LH766ScBQMTGxpa5b0ZGhlCpVMLBwUGkp6cLIYQ4cuSIUKvVwt7eXuTm5gqDwVDhKqs7d+4UAMSAAQPM24pDaubMmeZtOTk54sMPPxQxMTFCCCEMBoNo2LCh0Gq1Qoj7IeXg4CCuXLkiTCaT8Pb2FgDE/PnzhRBFCwACEEOHDhVC3J8j7uuvvzafQ6vVKnLlYqr72N1H9BDPPPMM3nrrLQDAokWL8OKLLwIA1qxZg2eeeabU/r/++itMJhP69+9vHuzQvXt3dO7cGTqdrlqWL3j88cfN3zdq1Ajjx4/HzZs3MWbMGDRv3hy3bt0qteZQu3bt4OnpCUmS0KxZMwDAgAEDANxfrqZ4kbrixR+feOIJSJIENzc36PV6nDx5Evfu3Xvk+okqgwMniB4iISEBly5dAlC0tlOxvXv3Ij09vdT+KlXZf/eJalxM09nZ2fz9xYsX0a9fP+j1esyaNQuzZ8/G3//+d2RmZpZ4zYProxXX6ODgUGbNBoMBADBo0CA0aNCgxHN6vd78OqKawJYU0UOcO3cOGo0G06dPhxACvr6+GDhwIIQQZS48GBAQAEmSkJiYaH4+OTkZJ0+ehKOjI/z9/R+5pgcHXmzcuBFXr17FmDFj8Prrr8Pb2xvZ2dmPdPziUYsTJ07Et99+i1WrVqF58+YYOXIkGjVq9EjHJqoshhTRQ5w8eRIGgwHdunWDTqdDamrqQ++P8vLywpQpU3Dv3j10794dgwcPRv/+/WE0GrFw4UI4OTlVa33FXYoxMTEICQlBQECAudsuLy+vSsecOXMmVCoVpk+fjkGDBiEgIACrVq3C8ePHq6tsItkYUkQPUXwNqVu3bjhx4gSMRmOFN/FGR0djyZIlcHNzw88//4wWLVogOjoa8+bNq/b6xo0bh//7v/+Di4sLkpOTMXjwYIwePRoAqryC6+OPP45Nmzahc+fO2L9/PyRJwvz587Fs2bJqrJxIHi56SEREisWWFBERKRZDioiIFIshRUREisWQIiIixWJIERGRYjGkiIhIsRhSRESkWAwpIiJSLIYUEREp1v8HGQghPl2zyjYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "B = 2 # upper bounds, assuming lower bounds is 0\n",
    "k = 1 # linear coeffient between drift rate and coherence\n",
    "coh = 0.064 # choose one coefficient\n",
    "driftRate = k * coh # drift rate\n",
    "sigma = 0.1\n",
    "\n",
    "evidence = B/2\n",
    "accumEvidence = []\n",
    "accumEvidence.append(evidence)\n",
    "\n",
    "i = 0\n",
    "while np.abs(evidence) < B:\n",
    "    # update frame number\n",
    "    i = i+1\n",
    "    # update evidence\n",
    "    # for every frame, we sample a \n",
    "    evidence = evidence + norm.rvs(loc=driftRate, scale=sigma, size=1)[0]\n",
    "    \n",
    "    accumEvidence.append(evidence)\n",
    "\n",
    "plt.plot(np.arange(i+1), accumEvidence, '-o')\n",
    "plt.axhline(B, color='k', label='Decision boundary')\n",
    "plt.xlabel('# of frame')\n",
    "plt.ylabel('Accumulated evidence')\n",
    "plt.ylim([0, B+1])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The standard DDM has four free paremters:\n",
    "* drift coefficient: $k$, $u = k * coh$ \n",
    "* decision boundary: $B$\n",
    "* initial bias: $alpha$ (default: 0.5). In some scenarios, the starting point could be biased towards either the upper or the lower bounds, indicating the prior bias in decision-making processes.\n",
    "* non-decision time(ndt): the whole decision process includes both drift-diffusion process and the duration of other non-decision processes (e.g., the time for a button press)\n",
    "\n",
    "In a typical experiment, we know the three varibles in each trial:\n",
    "* coherence \n",
    "* Correct/wrong($c_i$=1/0)\n",
    "* Reaction time (ms), $rt_i$\n",
    "\n",
    "We use the data to infere the four free parameters of each subject. Let $\\theta=(k, B, a, ndt)$, where $k$ is drift coefficient, $B$ is decision boundary，$a$ is initial bias, $ndt$ is non-decision time.\n",
    "\n",
    "We used the method of maximum likelihood estimate to obtain the parameters\n",
    "$$\n",
    "p(rt_i|\\theta, coh_i, c_i)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize\n",
    "from ddm import ddmpdf # ddm probility density function"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's load the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(300, 3)\n"
     ]
    }
   ],
   "source": [
    "data = np.loadtxt('exampledata.txt')\n",
    "print(data.shape)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The example data is a 300 (trials) x 3 array.\n",
    "* 1st column: coherence\n",
    "* 2nd column: reaction time\n",
    "* 3rd column: correct/wrong (1/0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define a negative log-likelihood objective functions\n",
    "def negloglikeli(params):\n",
    "    '''\n",
    "    <params>:(4,) array, drift coefficient, decision boundary, initial bias, non-decision time\n",
    "    '''\n",
    "    k = params[0] # drift coefficient\n",
    "    B = params[1] # decision boundary\n",
    "    a = params[2] # initial bias (0, 1)\n",
    "    ndt = params[3] # non-decision time\n",
    "\n",
    "    nTrial = data.shape[0] # how many trials\n",
    "    pp = np.empty(nTrial)\n",
    "    for i in range(nTrial): # loop trial\n",
    "        pp[i] = ddmpdf(k, a, B, ndt, coh=data[i, 0], rt=data[i, 1], correct=data[i, 2])\n",
    "\n",
    "    pp=0.999*pp + np.finfo(np.float32).eps # to avoid p=0\n",
    "    return -np.log(pp).sum() # take log, sum，add negative"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's optimize this objective function to obtain the fitted parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'minimize' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m/Users/ruyuan/Dropbox/backup/RuyuanZhang/Presentation/20221210 Neuroimagingworkshop/20221210_CITUWORKSHOP/PDM4.ipynb 单元格 11\u001b[0m in \u001b[0;36m<cell line: 1>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> <a href='vscode-notebook-cell:/Users/ruyuan/Dropbox/backup/RuyuanZhang/Presentation/20221210%20Neuroimagingworkshop/20221210_CITUWORKSHOP/PDM4.ipynb#X13sZmlsZQ%3D%3D?line=0'>1</a>\u001b[0m res \u001b[39m=\u001b[39m minimize(fun\u001b[39m=\u001b[39mnegloglikeli, x0\u001b[39m=\u001b[39m(\u001b[39m1.6\u001b[39m, \u001b[39m2.3\u001b[39m, \u001b[39m0.5\u001b[39m, \u001b[39m0.01\u001b[39m), bounds\u001b[39m=\u001b[39m((\u001b[39m0\u001b[39m, \u001b[39m20\u001b[39m), (\u001b[39m0\u001b[39m, \u001b[39m5\u001b[39m), (\u001b[39m0\u001b[39m, \u001b[39m1\u001b[39m), (\u001b[39m0\u001b[39m, \u001b[39m1\u001b[39m))) \u001b[39m# use bounds to constrain the range of the parameters\u001b[39;00m\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/ruyuan/Dropbox/backup/RuyuanZhang/Presentation/20221210%20Neuroimagingworkshop/20221210_CITUWORKSHOP/PDM4.ipynb#X13sZmlsZQ%3D%3D?line=2'>3</a>\u001b[0m \u001b[39mprint\u001b[39m(\u001b[39m'\u001b[39m\u001b[39m\\n\u001b[39;00m\u001b[39mfitted drift coefficient is \u001b[39m\u001b[39m'\u001b[39m, res\u001b[39m.\u001b[39mx[\u001b[39m0\u001b[39m])\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/ruyuan/Dropbox/backup/RuyuanZhang/Presentation/20221210%20Neuroimagingworkshop/20221210_CITUWORKSHOP/PDM4.ipynb#X13sZmlsZQ%3D%3D?line=3'>4</a>\u001b[0m \u001b[39mprint\u001b[39m(\u001b[39m'\u001b[39m\u001b[39mfitted decision boundary is \u001b[39m\u001b[39m'\u001b[39m, res\u001b[39m.\u001b[39mx[\u001b[39m1\u001b[39m])\n",
      "\u001b[0;31mNameError\u001b[0m: name 'minimize' is not defined"
     ]
    }
   ],
   "source": [
    "res = minimize(fun=negloglikeli, x0=(1.6, 2.3, 0.5, 0.01), bounds=((0, 20), (0, 5), (0, 1), (0, 1))) # use bounds to constrain the range of the parameters\n",
    "\n",
    "print('\\nfitted drift coefficient is ', res.x[0])\n",
    "print('fitted decision boundary is ', res.x[1])\n",
    "print('fitted initial bias is ', res.x[2])\n",
    "print('fitted nondecision time is ', res.x[3])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Summary\n",
    "* The standard DDM assumes a linear relationship between drift rate and coherence\n",
    "* We can use a predefined probility density function to calculate log likelihood of data given parameter\n",
    "* We use the method of maximum likelihood estimation to infer free paremeters\n",
    " "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mne",
   "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.10.4"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "6c5a876fa049f083f2e571ef64f1d25be9eec14fbc0627aabafac0fa7f44fa64"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}