2-Data-and-Spike-Stats.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 2: Data Input/Output and Point Process Statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Representing neural data in digital computers\n",
    "\n",
    "Computational neuroscience often involves a lot of data. In this lesson, you'll learn about:\n",
    "\n",
    "- different kinds of data: how they are represented digitally\n",
    "- how to read data from the disk into your program\n",
    "- how to store data from your program to the disk\n",
    "\n",
    "We'll consider three kinds of data:\n",
    "\n",
    "- [time series](#Time-series-data)\n",
    "- [point processes](#Point-processes)\n",
    "- [structured records](#Structured-records)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Time series data\n",
    "\n",
    "Typically represented as an **array** of measurements: \n",
    "\n",
    "$$\\mathbf{x} = \\{x_0, x_1, \\ldots, x_N\\}$$\n",
    "\n",
    "Can be **multichannel** if more than one measurement taken at a time. Each time point is now a **vector** ($\\vec{x}$):\n",
    "\n",
    "$$\\mathbf{X} = \\{\\vec{x}_0, \\vec{x}_1, \\ldots, \\vec{x}_N\\}$$\n",
    "\n",
    "Multichannel time series are represented as two-dimensional arrays. One dimension correponds to time and the other to the component of the measurement vector.\n",
    "\n",
    "Note that the \"channels\" can be repeated **trials** rather than simultaneous measurements."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Time series in Python\n",
    "\n",
    "We use `numpy` arrays to store single- and multichannel time series in Python. Let's look at an example using some Gaussian white noise.\n",
    "\n",
    "Gaussian noise is drawn from a normal distribution, and it's called white noise because it has equal power at all frequencies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load matplotlib inline mode\n",
    "%matplotlib inline\n",
    "\n",
    "# import some useful libraries\n",
    "import sys\n",
    "from pathlib import Path\n",
    "import numpy as np                # numerical analysis linear algebra\n",
    "import matplotlib.pyplot as plt   # plotting\n",
    "sys.path.insert(0,\"/project/psyc5270-cdm8j/comp-neurosci\")\n",
    "from comp_neurosci_uva import signal, data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)                  # set random seed so we all get the same results\n",
    "x = np.random.randn(1000)          # generate 100 random WN samples\n",
    "\n",
    "plt.plot(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can determine the number of elements in a 1D array using the `len()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"The number of points in x is:\", len(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Q**: Do you recall from the last exercise how to access subsets of a numpy array? In the code cell below, write an expression to evaluate the mean of the first 100 samples of `x`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numpy (and Python in general) supports **negative indexing**, which means that negative indices are interpreted as referencing elements from the **end** of the array. The following expression gives us the mean of the last 100 elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(x[-100:1000])\n",
    "# you can leave out the second index in the slice if it refers to the end of the array, so this is equivalent:\n",
    "np.mean(x[-100:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Multichannel time series\n",
    "\n",
    "For multichannel data, the array has two dimensions. There is a (weak) convention that the first dimension of the array represents time. That means each column represents a separate channel.\n",
    "\n",
    "Here is an example of a 3-channel array, again using Gaussian white noise. I've added some correlations between the channels to make things interesting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.random.multivariate_normal(mean=[0, 0, 0], cov=[[1.0, 0.2, 0.0],[0.2, 1.0, 0.1], [0.0, 0.1, 1.0]], size=1000)\n",
    "plt.plot(y);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how the plot now has three different-colored traces? Matplotlib assumes that time is the first dimension when you give it an array to plot.\n",
    "\n",
    "The size along each dimension of the array is called its **shape**. You can get the shape (and therefore the dimension of an array) using the `.shape` attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"The shape of y is:\", y.shape)\n",
    "# note that len returns the number of elements along the first dimension\n",
    "print(\"The number of time points in y is:\", len(y))\n",
    "print(\"The total size of y is:\", y.size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For multichannel data, we need two indices or slices to access values in the array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the first time point in the first channel. Note the comma.\n",
    "y[0, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use `:` to indicate all the values along one dimension. This gives all the values for the first channel\n",
    "plt.plot(y[:, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to get all the channels at a time point:\n",
    "print(\"y_0 =\", y[0, :])\n",
    "# you can leave out the trailing indices\n",
    "print(\"y_0 =\", y[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Q**: Calculate and plot the time series defined by the mean of all three channels (i.e., at each point in time, compute the average value across all channels)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Enter code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Input/Output \n",
    "\n",
    "Presumably you'd like to look at more than just randomly generated noise. So how do you get data from a recording into your program?\n",
    "\n",
    "Usually, data are stored long-term on your computer's drive or in the cloud. There are advantage to both approaches, which we'll discuss later. For now, we're going to retrieve some files stored in a shared directory. Once we do this, we'll see how to load the data from these files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### I/O for Time Series\n",
    "\n",
    "Unfortunately, there is no agreed-upon standard for storing time series data, so you'll often have to do some sleuthing.\n",
    "\n",
    "There are three major kinds of storage formats: text, binary, and custom."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Text\n",
    "\n",
    "One way of storing numbers is how you would write them (i.e., as **text**). \n",
    "\n",
    "When reading a text file, the main thing you need to know is how the elements are separated. \n",
    "\n",
    "For single-channel data, usually each number goes on its own line. \n",
    "\n",
    "For multi-channel data, there will be multiple numbers per line, typically separated by white space (tabs and/or spaces) or by commas.\n",
    "\n",
    "When you're storing data in text format, you also need to be mindful of the precision. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Reading text files\n",
    "\n",
    "One advantage of storing data as text is that it's human-readable. However, this isn't as much of an advantage as you might think. \n",
    "\n",
    "To see an example of an extracellular recording in text format, switch to the main Juptyer tab in your browser, navigate to the `data/io-examples` folder, and click on one of the files that ends in `.txt`\n",
    "\n",
    "Numpy can easily load single- and multi-channel data from text files using the `loadtxt` command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = data.data_path / \"io-examples\" / \"st11_1_2_A8.txt\"\n",
    "d = np.loadtxt(path)\n",
    "plt.plot(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Binary\n",
    "\n",
    "Storing numbers as text is very inefficient. Let's see why."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## text is read into Python first as a string:\n",
    "s = open(path, \"r\").readline().strip()\n",
    "print(\"the number as text:\", s)\n",
    "print(\"size of the text (in bytes):\", len(s))\n",
    "\n",
    "## to use it as a number, python has to parse the text\n",
    "f = float(s)\n",
    "print(\"the number as a float:\", f)\n",
    "print(\"size of a float (in bytes):\", d.dtype.itemsize)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Reading binary data\n",
    "\n",
    "Not only does text-formatted data take up a lot more space, it also requires additional work for Python to translate into a numerical representation that it can do math on (i.e., floats and ints).\n",
    "\n",
    "This inefficiency becomes a consideration for large datasets. Thus, we often want to store the data on disk in a binary format, i.e., the same format as it would be in memory. When you load binary data from disk, it's stored directly in an array, without any need for parsing.\n",
    "\n",
    "Another reason to store data in a binary format is that it allows your program to directly map the file into memory as a **memory map**. Instead of reading the entire file into memory, Python only reads the parts of the file that you actually use, when you want to use them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = data.data_path / \"io-examples\" / \"st11_1_2_A8.dat\"\n",
    "d = np.memmap(path, mode=\"r\", dtype='d')\n",
    "plt.plot(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Structured formats\n",
    "\n",
    "Text and raw binary formats both have shortcomings and tradeoffs.\n",
    "\n",
    "A shortcoming they both have in common is that it can be difficult to store metadata.\n",
    "\n",
    "Without metadata, it may be hard to know how to interpret the contents of the file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Some critical metadata we need for time series include:\n",
    "\n",
    "- sampling rate\n",
    "- dimensions of the array\n",
    "- ordering of the array (i.e., time first or last) and what's in each channel\n",
    "- measurement units\n",
    "\n",
    "**Q**: What other metadata do you think are important for time series data? Write a couple of ideas in the cell below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Structured data formats can be text- or binary-based.\n",
    "\n",
    "Some formats are used widely and are well-documented, like [Javascript Object Notation](http://json.org) or [HDF5](https://support.hdfgroup.org/HDF5/). These formats are likely to have well-supported Python packages for I/O.\n",
    "    \n",
    "Other formats are more obscure or proprietary, like Axon Binary Format (ABF) or [Elan](http://elan.lyon.inserm.fr/). It may be difficult to find Python support to read these files, though the situation is improving thanks to projects like [Neo IO](https://neo.readthedocs.io/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Point processes\n",
    "\n",
    "Typically represented as an ordered sequence of times in some interval from 0 to $T$:\n",
    "\n",
    "$$\\{0 \\leq t_0 < t_1 < \\ldots < t_N \\leq T\\}$$\n",
    "\n",
    "In contrast to time series, there is not a fixed relationship between the number of events and the duration of the analysis interval."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Point process data in Python\n",
    "\n",
    "Point processes are also typically stored in `numpy` arrays, but the elements of the array are event times, not measurements.\n",
    "\n",
    "Because point processes vary in the number of events, multi-channel point-processes are represented by **lists of arrays**, not by 2D arrays.\n",
    "\n",
    "Let's look at some examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from comp_neurosci_uva import pprox\n",
    "resp = data.load_pprox(Path(\"starling\", \"pprox\", \"st11_1_2_1\"))\n",
    "resp_A8 = pprox.select_stimulus(resp, \"A8\")\n",
    "resp_A8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variable `resp_A8` refers to a Python **list**. Lists are like arrays, but they can store heterogeneous data types. The syntax for accessing elements and slices is the same.\n",
    "\n",
    "**Q:** Using what you know from previous exercises, complete the following code cell to print out some information about the data: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"The number of trials is:\", ???)\n",
    "print(\"The number of events in trial 0 is:\", ???)\n",
    "print(\"The time of the first event in trial 2 is:\", ???)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## I/O for Point Processes\n",
    "\n",
    "Just as there is no agreed-upon standard for storing time series data, there is also no standard format for point-process data.\n",
    "\n",
    "Because point process data tend to be smaller than time series, text formats are more common than binary.\n",
    "\n",
    "A very simple text format is to put each trial (or channel) on a separate line and separate the events on each line with a space. Take a look at `data/io-examples/st_11_2_1_A8.txt` for an example.\n",
    "\n",
    "The [PySpike](http://mariomulansky.github.io/PySpike/) library has a function for loading data from such files, but we're going to write our own so that we can learn a bit about basic I/O in Python and looping."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a list where we will store our trials\n",
    "trials = []\n",
    "# open the file for reading\n",
    "fp = open(data.data_path / \"io-examples\" / \"st11_1_2_1_A8.txt\", mode=\"r\")\n",
    "# loop through the lines of the file with a for statement\n",
    "for line in fp:\n",
    "    # read the line into an array\n",
    "    arr = np.fromstring(line, sep=\" \")\n",
    "    # append the array to our list\n",
    "    trials.append(arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you learned how to program in Java or C or another low-level programming language, take a moment to appreciate how simple this task is in Python."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Checking our work\n",
    "\n",
    "**Q:** The `trials` list we created in the last code cell should be the same as `resp_A8`. In the cell below, complete **three assert** statements to check that this is true. I've provided you with one to get started. If you complete your task correctly, the cell will not emit any errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert len(trials) == len(resp_A8), \"The number of trials is not the same\"\n",
    "assert True == False, \"The total number of events is not the same\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Point process metadata\n",
    "\n",
    "As with time series data, it's important to keep track of metadata. Here are some important metadata that need to be associated with point process files:\n",
    "\n",
    "- type of event (e.g., spike, behavioral action, stimulus start/stop)\n",
    "- number of channels\n",
    "- unit scaling (e.g., milliseconds or seconds?)\n",
    "- start time\n",
    "- other experimental variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Structured Records\n",
    "\n",
    "In both point process and time series data, the elements of the arrays have been **homogeneous** (i.e., all the same type). What if that's not the case?\n",
    "\n",
    "The third (and final) kind of data we'll consider today consists of **records**. Each record in turn comprises **fields**, which may have different types.\n",
    "\n",
    "This kind of data is also called **tabular data**. If you're coming from the R world, you might think of this kind of data as a `data.frame`.\n",
    "\n",
    "It's common to encounter structured records when you have independent observations; for example, from different neurons or animals or populations. The fields in each record might include:\n",
    "\n",
    "- a unique identifier for the observation\n",
    "- group identifiers (e.g., cell, animal, population)\n",
    "- independent variables (e.g., sex, treatment, age)\n",
    "- dependent variable(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Structured records in Python\n",
    "\n",
    "Python used to lag pretty badly behind R for handling this kind of data, but we now have [pandas](http://pandas.pydata.org/), which is beginning to approach `numpy` in popularity and maturity.\n",
    "\n",
    "As with numpy, there is a convention for importing pandas: \n",
    "\n",
    "```python\n",
    "import pandas as pd\n",
    "```\n",
    "\n",
    "For a detailed introduction to pandas, take a look at [Chapter 3](https://jakevdp.github.io/PythonDataScienceHandbook/03.00-introduction-to-pandas.html) of the Python Data Science Handbook.\n",
    "\n",
    "There are two main concepts to understand in using pandas: `Series` and `DataFrames`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "A `Series` is essentially a column of a table. Like a numpy array, all the elements of a series are the same type. Unlike a numpy array, the indices of a `Series` do not have to be sequential integers, but can be any label you like.\n",
    "\n",
    "For example, here's a `Series` that might represent the ages of several subjects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "ages = pd.Series([391, 442, 183], index=['st11', 'st22', 'st231'])\n",
    "ages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can access the elements of a `Series` using the standard Python bracket syntax:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ages['st11']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "A `DataFrame` is a collection of `Series`, i.e. a table of columns. Here's how we might represent the ages and sexes of a set of subjects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sex = pd.Series(['M', 'F', 'M'], index=['st11', 'st22', 'st231'])\n",
    "subjects = pd.DataFrame({'age': ages, 'sex': sex})\n",
    "subjects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how we now have a table, which Jupyter renders nicely for us with the row and column indices are indicated in bold.\n",
    "\n",
    "The bracket syntax for `DataFrames` accesses **columns**. It's important to remember that this is different from numpy arrays, where a single index gives you a row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subjects['age']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "To index by row and column, you have to use the `loc` attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subjects.loc['st11', 'age']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equivalently, you can use `iloc` and the numerical indices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subjects.iloc[0, 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## I/O for Structured Records\n",
    "\n",
    "Structured record data is usually stored on disk in text-based formats. This is because human readability is often quite important. There are two very common formats:\n",
    "\n",
    "- In _comma-separated-value_ files, each record is on a separate line, and fields are separated by commas.\n",
    "- In _whitespace-delimited-value_ files, each record is on a separate line, and fields are separated by white space (tabs or spaces)\n",
    "\n",
    "In both kinds of files, it's common that the first line of the file is a header giving the name for each column.\n",
    "\n",
    "Take a look at `data/stimuli/motifs.csv` for an example of a comma-delimited file.\n",
    "\n",
    "One really good reason to use pandas is that it provides some nice I/O functions for these kinds of files. It's trivial to load tabular data into Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "motifs = pd.read_csv(data.data_path / \"starling\" / \"stimuli\" / \"motifs.csv\")\n",
    "motifs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "You can tell pandas that certain columns should be used as indices.\n",
    "\n",
    "This allows you to select a subset of rows using the `loc` syntax. For example, to see all the rows where `song` is equal to `A8`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "motifs = motifs.set_index(['song'])\n",
    "motifs.loc['A8']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exercise\n",
    "\n",
    "Write code in the cell below to compute the following:\n",
    "\n",
    "- the number of different songs\n",
    "- the number of motifs for each song\n",
    "- the average motif duration in each song"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spike Train Statistics\n",
    "\n",
    "As we've seen from previous examples, much of the data generated by the brain consists of spikes (i.e., action potentials). We all know that neurons spike when they are excited, but what does this mean quantitatively?\n",
    "\n",
    "In this part of the lesson, we'll dive into a fairly simple but foundational model that attempts to quantify spiking as the function of an underlying **rate**.\n",
    "\n",
    "Our goals are to:\n",
    "\n",
    "- understand homogeneous and inhomogeneous Poisson process models\n",
    "- be able to estimate the latent rate variable in Poisson models\n",
    "\n",
    "Follow along in the notebook during the lecture, and then work on the cells marked **Q** with help from your instructor. Submit the completed notebook to Collab."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### More point process math\n",
    "\n",
    "Recall that a point process is an (ordered) sequence of event times:\n",
    "\n",
    "$$X = \\{t_0, t_1, \\ldots, t_{N-1}\\}$$\n",
    "\n",
    "We can represent this as a function by making each spike a [Dirac delta function](https://en.wikipedia.org/wiki/Dirac_delta_function):\n",
    "\n",
    "$$\\rho(t) = \\sum_{i=0}^{N-1} \\delta(t - t_i)$$\n",
    "\n",
    "![spike_train](https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l5_spike_train.png \"spike train delta function\")\n",
    "\n",
    "Because the area under each delta function is 1, this allows us to count spikes or calculate any continuous function of a spike train through integration.\n",
    "\n",
    "For example, the **rate** is defined as the number of spikes $N$ that occurred in some interval divided by the duration of the interval, $T$:\n",
    "\n",
    "$$R = \\frac{N}{T} = \\frac{1}{T} \\int_{0}^{T} d\\tau\\; \\rho(\\tau)$$\n",
    "\n",
    "![spike_train_rate](https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l5_rate_integral.png \"spike train rate integral\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Spiking as a random variable\n",
    "\n",
    "Neural responses are **stochastic**. Even under \"identical\" conditions, spike trains will vary from trial to trial.\n",
    "\n",
    "In other words, the response $\\rho(t)$ is a **random variable**. The probability of observing a particular response is given by a **distribution**, $p(\\rho(t))$.\n",
    "\n",
    "We can also represent the stimulus $\\vec{s}(t)$ as a random variable with a distribution $p(\\vec{s}(t))$, even if it's under experimental control."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Review of probability\n",
    "\n",
    "Leaving out the notation indicating that both $\\rho$ and $\\vec{s}$ are functions of time, the following distributions are of interest:\n",
    "\n",
    "The joint distribution $p(\\rho, \\vec{s})$ is the probability of $\\rho$ and $\\vec{s}$ occurring together. This is likely to be a very small number because the space of $\\vec{s}$ and $\\rho$ are potentially quite large.\n",
    "\n",
    "Thus, we are often more interested in the conditional probability $p(\\rho|\\vec{s})$, which is the probability of observing $\\rho$ when $\\vec{s}$ is presented.\n",
    "\n",
    "The marginal probability $p(\\rho)$ indicates the probability of observing $\\rho$ irrespective of which stimulus was presented.\n",
    "\n",
    "There are two mathematical identities that allow us to convert between these distributions. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Conditional Probability\n",
    "\n",
    "$$p(\\rho,\\vec{s}) = p(\\rho|\\vec{s})p(\\vec{s}) = p(\\vec{s}|\\rho)p(\\rho)$$\n",
    "\n",
    "If $p(\\rho|\\vec{s}) = p(\\rho)$, then $\\rho$ and $\\vec{s}$ are **independent**. Independence means $p(\\rho)$ is the same regardless of what $\\vec{s}$ is. This means that for independent variables, the joint distribution is simply the product of the marginal distributions.\n",
    "\n",
    "$$p(\\rho,\\vec{s}) = p(\\rho)p(\\vec{s})$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Marginal Probability\n",
    "\n",
    "We can convert a joint probability distribution to a marginal probability distribution by summing (integrating) the probabilities of one of the variables.\n",
    "\n",
    "$$p(\\rho) = \\int_S d\\vec{s}\\; p(\\rho, \\vec{s}) = \\int_S d\\vec{s}\\; p(\\rho|\\vec{s}) p(\\vec{s})$$\n",
    "\n",
    "The $S$ under the integral means that we are summing up $p(\\rho,\\vec{s})$ over every possible stimulus. Remember that every probability density function has to integrate to 1.0:\n",
    "\n",
    "$$\\int_S d\\vec{s}\\; p(\\vec{s}) = 1$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Spike-train statistics\n",
    "\n",
    "Let's apply these concepts to spike trains:\n",
    "\n",
    "The probability of of a sequence of $N$ spikes $X = \\{t_0,\\ldots,t_{N-1}\\}$ is the joint probability density of all the individual spikes: \n",
    "\n",
    "$$p(t_0, t_1, \\ldots, t_{N-1})$$\n",
    "\n",
    "If the spikes are independent, then this joint distribution is simply the product of the distributions for each spike:\n",
    "\n",
    "$$p(t_0, \\ldots, t_{N-1}) = \\prod_{i=0}^{N-1}p(t_i)$$\n",
    "\n",
    "When each spike is independent of every other spike, we have a **Poisson process**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Homogeneous Poisson Processes\n",
    "\n",
    "If $t_i$ is independent of all the other spikes, what does it depend on?\n",
    "\n",
    "In the simplest case, $p(t_i)$ is a constant: the probability of observing a spike at any given time is a single number, which corresponds to the **rate** of spiking, $R$.\n",
    "\n",
    "If the rate is constant, the Poisson process is **homogeneous**. In this case, in an interval $(t_i, t_i + \\Delta)$, we would expect to observe $\\lambda = R\\Delta$ events. The distribution of the number of events we actually observe, $n$, is given by the Poisson distribution:\n",
    "\n",
    "$$p(n|\\lambda) = \\frac{\\lambda^n}{n!}\\exp(-\\lambda)$$\n",
    "\n",
    "The parameter $\\lambda$ is often called the **intensity** of the process."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's explore some properties of the Poisson distribution in Python. Here's a graph of the distribution from wikipedia:\n",
    "\n",
    "![poisson_distro](https://upload.wikimedia.org/wikipedia/commons/1/16/Poisson_pmf.svg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# load matplotlib inline mode\n",
    "%matplotlib widget\n",
    "\n",
    "# import some useful libraries\n",
    "import numpy as np                # numerical analysis linear algebra\n",
    "# notebook modules\n",
    "import IPython\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "from IPython.display import display\n",
    "import ipywidgets as widgets\n",
    "\n",
    "# import the poisson distro\n",
    "from comp_neurosci_uva import dists"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# first, let's define our *support*: the values over which we want to evalute p(n):\n",
    "supp = np.arange(0, 20)\n",
    "\n",
    "# next, we *instantiate* the distribution object with our parameter lambda\n",
    "dist = dists.poisson(1.0)\n",
    "\n",
    "# you can get the probability of any value in the distribution with .pmf. Note that we have to use\n",
    "# pmf (probability mass function) rather than pdf.\n",
    "print(\"p(5|lambda=1) =\", dist.pmf(5))\n",
    "\n",
    "# we can also evaluate the distribution over a vector of numbers\n",
    "prob = dist.pmf(supp)\n",
    "\n",
    "# and plot the distribution with plt.plot\n",
    "plt.plot(supp, prob, lw=3);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Q**: In the cell below, try to evaluate the `prob` distribution for negative or non-integral numbers. Given the definition of the Poisson distribution, why is this the case?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Q** Recall that $\\lambda = R\\Delta$. If $R = 1$ Hz and you steadily reduce $\\Delta$ from 1.0 s to 1.0 microseconds, what is the probability of observing one spike in that interval? Write a *for loop* to evaluate and plot this. Does the result make sense? What is the probability of a spike occurring at some *exact* time?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Delta = [1e0, 5e-1, 1e-1, 5e-2, 1e-2, 5e-3, 1e-3, 5e-4, 1e-4, 5e-5, 1e-5, 5e-6, 1e-6]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In addition to plotting the probability distribution, Python can generate random samples (i.e, **draw**) from the Poisson distribution, too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# note: 'lambda' is a reserved symbol in python\n",
    "lam = 1.0\n",
    "dist = dists.poisson(lam)\n",
    "# rvs stands for random value(s)\n",
    "n = dist.rvs(1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Q** In the cell below, calculate the sample mean and standard deviation for n, then try with a few other values of lambda."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"The sample mean when lambda=\", lam, \"is\", np.mean(n))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q**: Based on your exploratory analysis, fill out the following identities:\n",
    "\n",
    "- Mean: $\\mu =$\n",
    "- Standard deviation: $\\sigma =$\n",
    "- Fano factor $\\sigma^2/\\mu =$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### From Poisson distribution to Poisson process\n",
    "\n",
    "How can we generate a series of spike times from the Poisson distribution? The trick is to divide your response interval up into a set of smaller intervals (or **bins**) such that the probability of observing more than one spike in a single bin is very small, then draw from $p(n|\\lambda)$ for each bin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "T     = 100     # s\n",
    "rate  = 4.0     # Hz\n",
    "Delta = 0.005   # s\n",
    "dist   = dists.poisson(rate * Delta)\n",
    "hom_spikes = dist.rvs(int(T / Delta))\n",
    "bins   = np.arange(0, T, Delta)\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(6, 3))\n",
    "axes.plot(bins, hom_spikes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "The plot of the `spikes` array should only show `0` and `1` values. If it doesn't, try adjusting the `Delta` variable in the code cell above. What direction does it need to change to fix the problem? Why?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Spike arrays and spike times\n",
    "\n",
    "You can think of the `spikes` array as a sort of time series representation of the point process. \n",
    "\n",
    "To get the actual spike times, we need to find the bins where there is a spike and then look up the times in the `bins` array. This allows us to generate a *raster plot*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hom_spike_i = np.nonzero(hom_spikes)[0]\n",
    "hom_spike_t = bins[hom_spike_i]\n",
    "# here's one way to plot a raster of spike times\n",
    "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(6, 3))\n",
    "axes.plot(hom_spike_t, np.zeros_like(hom_spike_t), \"|\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### More spike train statistics\n",
    "\n",
    "**Q** Here's a slightly harder question about the properties of Poisson processes. Calculate the interspike *intervals* from `spike_t` (hint: look at the documentation for `np.diff`), then plot a histogram. What function does this look like? Calculate the sample mean and variance. How do these relate to the rate of the process? What is the coefficient of variation of the interspike distribution (CV = $\\mu / \\sigma$)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Inhomogeneous Poisson Processes\n",
    "\n",
    "It's also possible for the rate of a Poisson process to vary in time; that is, for $\\lambda$ to be a function of $t$.\n",
    "\n",
    "$$p(n|\\lambda(t)) = \\frac{\\lambda(t)^n}{n!}\\exp(-\\lambda(t))$$\n",
    "\n",
    "As before, we need to discretize time and determine the probability that there is a spike in some interval $(t, t + \\Delta)$; the only difference is that some intervals are more likely to have spikes than others.\n",
    "\n",
    "We could simulate an inhomogeneous Poisson process in much the same way as we did above, but we need to vary $\\lambda$ in each bin.\n",
    "\n",
    "Let's look at a simple example where the intensity linearly ramps up from zero and then back down."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "T     = 100     # s\n",
    "Delta = 0.001   # s\n",
    "N     = int(T / Delta)\n",
    "bins  = np.arange(0, T, Delta)\n",
    "# rate is now a function of time\n",
    "inh_rate  = np.concatenate([np.linspace(0.0, 4.0, N//2),\n",
    "                            np.linspace(4.0, 0.0, N//2)])\n",
    "\n",
    "# generate N values from a uniform distribution\n",
    "rand = dists.uniform().rvs(N)\n",
    "# this is an alternative method of simulating spiking based on the Bernoulli distribution\n",
    "# compare each value to lambda = rate * Delta; if it's greater, then the bin gets a spike\n",
    "lam  = inh_rate * Delta\n",
    "inh_spikes = (inh_rate * Delta) > rand\n",
    "inh_spike_i = np.nonzero(inh_spikes)[0]\n",
    "inh_spike_t = bins[inh_spike_i]\n",
    "# here's one way to plot a raster of spike times\n",
    "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(9, 3))\n",
    "axes.plot(bins, inh_rate)\n",
    "axes.plot(inh_spike_t, np.zeros_like(inh_spike_t), \"k|\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Calculate the interspike intervals for this spike train (`inh_spike_t`) and plot the distribution. What are the mean, standard deviation, and CV? Is the relationship with the rate the same as you discovered for the homogeneous process?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Estimating spike rates\n",
    "\n",
    "Let's think about how we can estimate $\\lambda$ for Poisson processes.\n",
    "\n",
    "If we assume that the process is homogeneous over each trial, then we have a simple observational model where the number of spikes is a random sample from the Poisson distribution.\n",
    "\n",
    "$$p(y_i|\\lambda) = \\frac{\\lambda^n}{n!}\\exp(-\\lambda)$$\n",
    "\n",
    "Given a set of trials, we can estimate $\\lambda$ from the sample mean of the spike count:\n",
    "\n",
    "$$\\hat{\\lambda} = \\frac{T}{N}\\sum_i y_i$$\n",
    "\n",
    "where $N$ is the number of time bins and $T$ is the duration of the observation interval.\n",
    "\n",
    "**Q**: What is the estimated intensity of the homogeneous spike train we generated above (`hom_spike_t`)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The problem is a lot trickier if the process is inhomogeneous, because now we're trying to estimate a continuous function of time, $\\lambda(t)$.\n",
    "\n",
    "In this setting, people usually talk about rate rather than intensity ($\\lambda$), so we'll use $r(t)$ from here on.\n",
    "\n",
    "The issue we confront is that $r(t)$ is a continuous function. We can discretize it into small intervals of $(t, t + \\Delta)$ and count the number of spike in each interval, but as we make $\\Delta$ smaller to get higher temporal resolution, we reach the point at which each bin has either one or zero spikes, which doesn't tell us much about the rate. We can address this problem by averaging across multiple trials. If we use $\\langle \\rangle$ to denote averaging across trials, this looks like:\n",
    "\n",
    "$$r(t) = \\frac{1}{\\Delta} \\int_t^{t+\\Delta} d\\tau\\; \\langle \\rho(t) \\rangle$$\n",
    "\n",
    "You hopefully can see that as $\\Delta$ gets smaller, the number of trials you need to average to get a smooth function gets larger. So part of our problem is to determine what $\\Delta$ should be. More practically, at what time scale do we think the rate is changing?\n",
    "\n",
    "There are a number of different ways of approximating $r(t)$. We'll look at a couple."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Spike time histogram\n",
    "\n",
    "For historical reasons, this is also called a peri-stimulus spike time histogram (PSTH), even when there isn't a stimulus.\n",
    "\n",
    "The simplest way of approximating the rate is to divide the interval up into a fixed number of bins of duration $\\Delta$ and count how many spikes occurred in each bin. The rate is simply the number of spikes divided by $\\Delta$.\n",
    "\n",
    "The main problem with histograms is that setting the bin size is largely subjective. Try adjusting the bin count (which is just the inverse of the bin size) variable and see what gives you the best tradeoff between variability and temporal resolution.\n",
    "\n",
    "There is still active development of new methods for adaptively setting bin sizes in timing histograms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(9, 3))\n",
    "bin_size = 5.0\n",
    "bin_count = widgets.IntSlider(value=int(T / bin_size), min=1, max=100, step=1, \n",
    "                              description=\"bin count:\", continuous_update=True)\n",
    "display(bin_count)\n",
    "\n",
    "ax.plot(bins, inh_rate)\n",
    "ax.plot(inh_spike_t, np.zeros_like(inh_spike_t), \"k|\")\n",
    "r_est, edges  = np.histogram(inh_spike_t, bins=np.arange(0, T + bin_size, bin_size))\n",
    "p = ax.step(edges[1:], r_est / bin_size)\n",
    "\n",
    "def update(bin_count):\n",
    "    bin_size = T / bin_count\n",
    "    r_est, edges  = np.histogram(inh_spike_t, bins=np.arange(0, T + bin_size, bin_size))\n",
    "    print(\"bin size:\", bin_size)\n",
    "    p[0].set_data(edges[1:], r_est / bin_size)\n",
    "    \n",
    "widgets.interactive_output(update, {\"bin_count\": bin_count})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Smoothing\n",
    "\n",
    "Another problem with the PSTH is that the count in each bin can depend very strongly on where the edges of the bins are.\n",
    "\n",
    "One solution to this problem is to use a **sliding window**. The simplest window is simply a square function with a defined width and a total area equal to 1.0.\n",
    "\n",
    "For example, here's a 10 ms window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(3, 3))\n",
    "window_size = 2.0\n",
    "w_t = np.arange(-10.0, 10.0, Delta)\n",
    "w = np.zeros_like(w_t)\n",
    "w[(-window_size/2 < w_t) & (w_t <window_size/2)] = 1. / window_size\n",
    "ax.plot(w_t, w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Convolution\n",
    "\n",
    "We can express the sliding window operation as a sum of the window function for the values of the spike times.\n",
    "\n",
    "$$r(t) \\approx \\sum_{i=0}^{N-1} w(t - t_i)$$\n",
    "\n",
    "This is equivalent to doing an integral over the response function:\n",
    "\n",
    "$$r(t) \\approx \\int_{-\\infty}^{\\infty} d\\tau\\; w(\\tau) \\rho(t - \\tau)$$\n",
    "\n",
    "This integral is also called a linear **convolution** or filter, and we'll be seeing a lot of them.\n",
    "\n",
    "Numpy has a function that can calculate this convolution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(9, 3))\n",
    "r_est = np.convolve(inh_spikes, w, mode='same')\n",
    "ax.plot(bins, inh_rate)\n",
    "ax.plot(inh_spike_t, np.zeros_like(inh_spike_t), \"k|\")\n",
    "ax.plot(bins, r_est)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "You can use any function as a window as long as it goes to zero outside $\\tau = 0$ and its integral is 1.0.\n",
    "\n",
    "A popular choice is to use a Gaussian, which smooths the function by downweighting points further away from $\\tau = 0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(3, 3))\n",
    "sigma = 2.0\n",
    "w = 1 / np.sqrt(2 * np.pi) / sigma * np.exp(-w_t**2 / 2 / sigma**2)\n",
    "plt.plot(w_t, w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(9, 3))\n",
    "sigma_w = widgets.FloatSlider(value=sigma, min=0.1, max=5.0, step=0.1,\n",
    "                              description=\"sigma:\", continuous_update=False)\n",
    "display(sigma_w)\n",
    "\n",
    "ax.plot(bins, inh_rate)\n",
    "ax.plot(inh_spike_t, np.zeros_like(inh_spike_t), \"k|\")\n",
    "r_est = np.convolve(inh_spikes, w, mode='same')\n",
    "p = ax.plot(bins, r_est)\n",
    "\n",
    "def update(sigma):\n",
    "    w_T = 10\n",
    "    w_t = np.arange(-w_T, w_T, Delta)\n",
    "    w = 1 / np.sqrt(2 * np.pi) / sigma * np.exp(-w_t**2 / 2 / sigma**2)\n",
    "    r_est = np.convolve(inh_spikes, w, mode='same')\n",
    "    p[0].set_data(bins, r_est)\n",
    "    \n",
    "widgets.interactive_output(update, {\"sigma\": sigma_w})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Averaging trials\n",
    "\n",
    "Hopefully, these exercises have illustrated the fundamental tradeoff between variance and temporal resolution. As you increase the bin (or window) size ($\\Delta$), the estimated rate becomes less variable, but the temporal resolution decreases. Thus, smoothing can interfere with detecting rapid changes in the underlying rate function.\n",
    "\n",
    "As noted above, the solution to this problem is to average across trials. In essence, this gives you multiple independent estimates of the rate at any given instant, thereby reducing the amount of smoothing you need.\n",
    "\n",
    "We can represent spiking responses over multiple trials using a two-dimensional array. Here's a simulation of 10 trials to the above rate function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we're keeping all the parameters the same as above, but resetting the random seed\n",
    "np.random.seed(1)\n",
    "n_trials = 10\n",
    "\n",
    "# it's often a good idea to initialize the array first, then fill it up\n",
    "# recall that N is the number of bins\n",
    "mt_spikes = np.zeros((N, n_trials))\n",
    "# we're going to store the spike times in a ragged array (i.e. a list of arrays)\n",
    "mt_spikes_t = []\n",
    "for trial in range(n_trials):\n",
    "    # generate N values from a uniform distribution\n",
    "    rand = dists.uniform().rvs(N)\n",
    "    lam  = inh_rate * Delta\n",
    "    spikes = (inh_rate * Delta) > rand\n",
    "    idx = np.nonzero(spikes)[0]\n",
    "    mt_spikes[:, trial] = spikes\n",
    "    mt_spikes_t.append(bins[idx])\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(9, 3))\n",
    "axes.plot(bins, inh_rate)    \n",
    "for i, trial in enumerate(mt_spikes_t):\n",
    "    axes.plot(trial, np.ones_like(trial) * i / n_trials, \"k|\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The rate function is calculated by averaging smoothed estimates across trials. Compare the plot below to the smoothed estimate based on 1 trial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma = 2.0\n",
    "w_T = 10\n",
    "w_t = np.arange(-w_T, w_T, Delta)\n",
    "w = 1 / np.sqrt(2 * np.pi) / sigma * np.exp(-w_t**2 / 2 / sigma**2)\n",
    "r_est = np.convolve(mt_spikes.mean(1), w, mode='same')\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(9, 3))\n",
    "axes.plot(bins, inh_rate)\n",
    "axes.plot(bins, r_est)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Putting it together\n",
    "\n",
    "For the rest of the class period, work on the following exercises.\n",
    "\n",
    "#### Exercise 1\n",
    "\n",
    "Consider an inhomogeneous Poisson process with a time-varying rate specified by:\n",
    "\n",
    "$$r(t) = r_\\mathrm{max}[\\sin(2\\pi\\omega t) + 1]$$\n",
    "\n",
    "Use $\\omega$ = 2.1 Hz, $r_\\mathrm{max} = 50$ Hz, and a response interval from 0 to 2 s."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Plot r(t) with a time step of 1 ms (0.001 s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Generate 20 independently simulated spike trains and plot them as rasters. There is code in previous notebooks you can use to make the raster plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q**. Using a bin size of 10 ms, calculate the PSTHs averaged from the first 10 trials and the last 10 trials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Now simulate 1000 trials and calculate a PSTH from the first and second half. How do these PSTHs relate to $\\lambda(t)$? Do more trials give you a more precise estimate of the rate?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Calculate a PSTH *for each trial* of the 1000-trial simulation. This will yield a 1000 x 200 array, with trials along one axis and time along the other. For each time bin, calculate and plot the mean, variance, and Fano factor. Each of these will be a 200-element array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q** Repeat the analysis in the previous question but with a bin size of 150 ms. What happens to the Fano factor? Why is it important to choose a bin size such that the rate is not changing much within each bin?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.4"
  },
  "rise": {
   "scroll": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}