SDA-Recitation 3.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.signal as sig\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "import logging\n",
    "logging.basicConfig(level=logging.INFO, format='%(levelname)s: %(asctime)s: %(message)s')\n",
    "\n",
    "%matplotlib widget\n",
    "# %matplotlib inline\n",
    "\n",
    "import ipywidgets as widgets\n",
    "\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "# plt.style.use(\"https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Single spike train analysis:\n",
    "### ISI - Time difference between adjacent spikes\n",
    "### TIH - Histogram of ISI\n",
    "### Survivor function - The probability of surviving at least to time t (before an event) \n",
    "### $ Survivor(t) = 1 - CDF(ISI(t)) = 1 - \\sum_{i=1}^{t}ISI(i) $\n",
    "### Hazard function - The probability of an event at a time n after the previous event\n",
    "### $ Hazard(t) = \\frac{ISI(t)}{Survivor(t)} $\n",
    "### ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### properties of the hazard function:\n",
    "### $ h(t) \\geq 0 $\n",
    "### $ h(t) $ is not probability so $ \\int_{}h(t) \\longrightarrow \\infty $\n",
    "### very noisy for long ISIs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "prob = 0.3\n",
    "prob_vec = np.array([prob]*10)\n",
    "ISI_vals = [prob]\n",
    "rest = 1\n",
    "for i in range(9):\n",
    "    rest -= ISI_vals[i]\n",
    "    ISI_vals.append(prob*rest)\n",
    "\n",
    "CDF = np.hstack((0,np.cumsum(ISI_vals)[:-1]))\n",
    "survivor = 1-CDF\n",
    "hazard = ISI_vals/survivor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "41b5a0e3e68144cc967253c5b06d90f8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(nrows=5, figsize=(10,8))\n",
    "ax[0].bar(np.arange(10), prob_vec)\n",
    "ax[0].set_ylim([0, prob+0.1])\n",
    "ax[0].set_title('probability of a spike')\n",
    "ax[1].bar(np.arange(10), np.array(ISI_vals))\n",
    "ax[1].set_ylim([0, prob+0.1])\n",
    "ax[1].set_title('TIH')\n",
    "ax[2].bar(np.arange(10), CDF)\n",
    "ax[2].set_title('CDF(ISI)')\n",
    "ax[3].bar(np.arange(10), survivor)\n",
    "ax[3].set_title('Survivor function')\n",
    "ax[4].bar(np.arange(10), hazard)\n",
    "ax[4].set_ylim([0, prob+0.1])\n",
    "ax[4].set_title('Hazard function')\n",
    "fig.tight_layout()\n",
    "for x in ax:\n",
    "    x.set_xlabel('bins')\n",
    "    x.set_ylabel('p')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d1f7dd262e934f7fb9d44b327ebe9040",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Generate a poisson neuron with refractory period\n",
    "samp = 1000\n",
    "rate = 20 / samp \n",
    "duration = 3600    \n",
    "ref_period = 10 #in ms\n",
    "spk_vec=np.zeros(duration*samp)\n",
    "rand_vec=np.random.uniform(low=0, high=1, size=duration*samp)\n",
    "idx=0\n",
    "while idx<duration*samp:\n",
    "    if rand_vec[idx] <= rate:\n",
    "        spk_vec[idx] = 1\n",
    "        idx = idx+ref_period+1\n",
    "    else:\n",
    "        idx=idx+1\n",
    "\n",
    "# first we calculate TIH and TIH probabilities\n",
    "spk_times = np.where(spk_vec==1)[0]\n",
    "isi_vec = np.diff(spk_times)\n",
    "tih, _ = np.histogram(isi_vec, bins=np.arange(1,400,1), density=True)\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (10,3))\n",
    "ax.bar(np.arange(1,399,1),tih)\n",
    "__ = ax.set_title('TIH')\n",
    "ax.set_xlabel('Interval [ms]')\n",
    "ax.set_ylabel('ISI(t)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a72d6e1df46849f3857dc3aae52b3d87",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "survivor = 1 - np.cumsum(tih)\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (10,3))\n",
    "ax.bar(np.arange(1,399,1),survivor)\n",
    "__ = ax.set_title('Survivor')\n",
    "ax.set_xlabel('Interval [ms]')\n",
    "ax.set_ylabel('Survivor(t-1)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "18dce4708c204db9bd4c5b039a3a0e20",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# calculatin hazard function \n",
    "hazard = tih[1:] / survivor[0:-1]\n",
    "\n",
    "#plotting hazard function\n",
    "fig, ax = plt.subplots(figsize = (10,3))\n",
    "ax.plot(hazard)\n",
    "ax.set_title('Hazard function')\n",
    "ax.set_xlabel('Interval [ms]')\n",
    "ax.set_ylabel('ISI(t)/Survivor(t-1)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Autocorrelation definition:\n",
    "### $ Q(x,y)[n] = \\sum_{i=-\\infty}^{\\infty} x(i)y(i+n)   $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = 1/samp\n",
    "vec_length = len(spk_vec)\n",
    "autocorr = sig.correlate(spk_vec, spk_vec, mode='full', method='auto')\n",
    "autocorr[(vec_length-1):(vec_length+1)] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f1d81d05d90f48a0bfebdd4522980dc2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "max_offset = 200\n",
    "fig, ax = plt.subplots(figsize=(10,3))\n",
    "ax.plot(np.arange( (-max_offset)/samp,(max_offset)/samp,dt),autocorr[vec_length-max_offset:vec_length+max_offset])\n",
    "ax.set_title('autocorrelation poisson neuron with refractory period')\n",
    "__ = ax.set_ylabel('count')\n",
    "__ = ax.set_xlabel('offset(sec)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Autocorrelation is symmetric"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}