5-Decoding-Models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Week 7: Decoding Models and Information Theory\n",
    "\n",
    "Up until now, we have mostly been considering models of how various stimulus features are **encoded** in neural activity.\n",
    "\n",
    "We'll now consider the other side of the coin: given a pattern of neural activity, what can we tell about what stimulus was presented?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# load matplotlib inline mode\n",
    "%matplotlib inline\n",
    "\n",
    "# import some useful libraries\n",
    "import sys\n",
    "import numpy as np                # numerical analysis linear algebra\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt   # plotting\n",
    "from comp_neurosci_uva import dists\n",
    "\n",
    "# set some style options\n",
    "mpl.rcParams['image.origin'] = 'lower'\n",
    "mpl.rcParams['image.aspect'] = 'auto'\n",
    "mpl.rcParams['image.cmap'] = 'jet'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Encoding vs decoding models\n",
    "\n",
    "We can think of the sensory stimulus and the response it evokes as being random variables. In each trial, a stimulus is selected from the stimulus distribution, and this results in some spiking pattern which comes from a large (but finite) set of possible responses.\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l9_distributions.png\" alt=\"V1 tuning\" style=\"width: 500px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Encoding models\n",
    "\n",
    "Goal is to model how the response depends on which stimulus was presented. More formally, a probabilistic model of the response conditional on the stimulus, $p(r|s)$. \n",
    "\n",
    "RF models covered in our last lesson are an example. The response at time $t_i$ is modeled as a weighted sum of the present and past values of the stimulus, plus some random noise, $\\varepsilon_i$.\n",
    "\n",
    "$$\n",
    "r_i = \\mathbf{s}_i \\cdot \\mathbf{h} + \\varepsilon_i\n",
    "$$\n",
    "\n",
    "If the errors are normally distributed with variance $\\sigma^2$, then the conditional probability is simply\n",
    "\n",
    "$$\n",
    "p(r|\\mathbf{s}) = \\mathrm{N}(\\mathbf{s} \\cdot \\mathbf{h}, \\sigma^2)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Decoding models\n",
    "\n",
    "In contrast, a decoding model attempts to reconstruct the stimulus based on the response. More formally, it represents $p(\\mathbf{s}|r)$, the probability that a stimulus was presented conditional on the response that was observed.\n",
    "\n",
    "Decoding models are related to encoding models through Bayes' rule:\n",
    "\n",
    "$$\n",
    "p(\\mathbf{s}|r) = \\frac{p(r|\\mathbf{s})p(\\mathbf{s})}{p(r)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Signal detection\n",
    "\n",
    "Let's consider a simple problem in which we have to detect a signal against a noisy background. We'll start by describing the encoding model and then develop a decoding model.\n",
    "\n",
    "On a given trial, the stimulus $s$ is drawn from a random variable $S$ that can either be present ($s = 1$) or absent ($s = 0$). The distribution $p(s)$ gives the *prior probability* that the signal will be present $p(s = 1)$ or absent $p(s = 0)$. Of course, because there are only two possible outcomes, we only need to know $p(s = 1)$ to determine $p(s = 0)$ so for notational convenience we'll just use $p(s)$ to indicate the probability of there being a signal.\n",
    "\n",
    "In each trial, the neuron that we're monitoring generates a response that we'll summarize by the rate $r$, which is drawn from a random variable $R$ that can be any non-negative real number. In contrast to distributions over $S$ that can be expressed as a single number, distributions over $R$ will be *functions* of $r$ that have to obey the rule that all the probabilities sum to 1.0 ($\\int_R p(r) dr = 1$). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Our encoding model is defined by the distribution of responses conditional on whether or not the stimulus is present. In signal detection theory, we assume that these conditional distributions are Gaussian:\n",
    "\n",
    "$$\n",
    "p(r|s) = \\mathrm{N}(\\mu, \\sigma^2)\n",
    "$$\n",
    "\n",
    "Furthermore, let's say that the average response $\\mu$ is a simple linear function of whether or not the stimulus is present:\n",
    "\n",
    "$$\n",
    "\\mu = \\beta_0 + \\beta_1 s\n",
    "$$\n",
    "\n",
    "Let's look at the distributions in Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "beta  = [10, 5]\n",
    "sigma = np.sqrt(5)\n",
    "\n",
    "# support of the distribution\n",
    "r = np.arange(0, 30, 0.1)\n",
    "# probability density of responses to noise\n",
    "pr_noise = dists.normal(mean=beta[0], std=sigma)\n",
    "# probability density of responses to signal\n",
    "pr_signal = dists.normal(mean=beta[0] + beta[1], std=sigma)\n",
    "plt.plot(r, pr_noise.pdf(r), lw=2, label=\"p(r|s=0)\")\n",
    "plt.plot(r, pr_signal.pdf(r), lw=2, label=\"p(r|s=1)\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"Rate\")\n",
    "plt.ylabel(\"Probability\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now let's use Bayes's Rule to derive the decoding model, $p(s|r)$.\n",
    "\n",
    "First, we need to know the *prior probability* that a stimulus will be present. If you're the experimenter, you get to decide what this is. If you're the receiver, this is your *belief* about how often the signal is present. This belief can be based on other knowledge you have about the world.\n",
    "\n",
    "Next, we need to calculate the marginal response probability $p(r)$, which is done by integrating the conditional response distributions:\n",
    "\n",
    "\\begin{align}\n",
    "p(r) & = \\sum_{s \\in S} p_{S,R}(s, r) \\\\\n",
    "     & = \\sum_{s \\in S} p_R(r|s) p_S(s) \\\\\n",
    "     & = p_R(r|s=0)p_S(s=0) + p_R(r|s=1)p_S(s=1)\n",
    "\\end{align}\n",
    "\n",
    "If the probability of the signal being present $p_S(s=1) = 0.5$, then\n",
    "\n",
    "$$p(r) = \\frac{1}{2} p(r|s=0) + \\frac{1}{2} p(r|s=1)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# probability of a stimulus\n",
    "ps = 0.5\n",
    "# marginal probability density of the response\n",
    "pr = (1 - ps) * pr_noise.pdf(r) + ps * pr_signal.pdf(r)\n",
    "plt.plot(r, pr)\n",
    "plt.xlabel(\"Rate\")\n",
    "plt.ylabel(\"p(r)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Finally, we'll use Bayes's Rule to calculate $p(s|r)$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bayes's Rule: probability of a stimulus conditional on response\n",
    "psr = (pr_signal.pdf(r) * ps) / pr\n",
    "plt.plot(r, psr)\n",
    "plt.xlabel(\"Rate\")\n",
    "plt.ylabel(\"p(s=1|r)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Important Note*: the above plot is not a probability *distribution* like the plots of $p(r|s)$ and $p(r)$ shown earlier. For a given value of $r$, $p(s|r)$ is just a single\n",
    "number, the posterior probability that there was a signal (conditional on the value of $r$). Because it's just a number, we can plot it as a function of $r$ to visualize how this probability\n",
    "would change depending on the value of $r$ we observe in any given trial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Write a **function** to calculate $p(s|r)$ for any value of $p(s)$. Plot $p(s|r)$ when $p(s) = 0.1$ and when $p(s) = 0.9$. How does the decoding model shift its output based on this prior knowledge?\n",
    "\n",
    "Consult the Software Carpentry [module for functions](http://swcarpentry.github.io/python-novice-inflammation/08-func/index.html) if you need a refresher."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## ROC analysis\n",
    "\n",
    "The foregoing approach is called a **naive Bayesian decoder**. It works well if you happen to know or control the stimulus probabilities, and if you can collect enough data to accurately model $p(r|s)$.\n",
    "\n",
    "However, in the real world, it may not be possible to know what the prior probability of a stimulus is. In this case, the only strategy is to set a **criterion**, or threshold, for deciding whether the signal is present or not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(r, pr_noise.pdf(r), lw=2, label=\"p(r|s=0)\")\n",
    "plt.plot(r, pr_signal.pdf(r), lw=2, label=\"p(r|s=1)\")\n",
    "plt.vlines(12, 0, 0.2, label=\"criterion\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"Rate\")\n",
    "plt.ylabel(\"Probability\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the response is above the criterion, the signal is assumed to be present, and if it's below, it's assumed to be absent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exercise 2\n",
    "\n",
    "With $p(s) = 0.5$, what criterion would you choose to minimize the probability of an error? Give your best guess and explain why. We'll come up with a more formal solution later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "This problem is at the heart of **signal detection theory** (SDT). SDT applies not only to neurons but to the receiver at the end of any noisy communication channel. A core concept is the **reciever operating characteristic** (ROC), which describes how the receiver's ability to accurately detect the signal depends on where it sets its criterion.\n",
    "\n",
    "Because the criterion and the signal are both binary, this leads to a 2-by-2 table of possible outcomes. Let's let $y = 0$ represent a response of \"no signal\" (i.e. below criterion) and $y = 1$ represent a response of \"signal\" (above criterion):\n",
    "\n",
    "\n",
    "|  -  | s = 0 | s = 1 |\n",
    "|----|------|------|\n",
    "|y = 0  | correct rejection| miss |\n",
    "|y = 1  | false alarm | hit |\n",
    "\n",
    "It should be clear that false alarms and correct rejections are complementary, and so are misses and hits. Thus, we only need two values from the table, the false alarm rate and the hit rate. If we let $\\gamma$ be the criterion, then these are:\n",
    "\n",
    "\\begin{align}\n",
    "p(FA) & = Pr(r > \\gamma|s = 0) \\\\\n",
    "p(H) & = Pr(r > \\gamma|s = 1)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "As the criterion is changed, the false alarm and hit rate also change. You should be able to see this by looking at the plot of the two conditional response distributions above.\n",
    "\n",
    "Let's illustrate in Python. An important concept we need to cover first is the **cumulative probability distribution** or **cdf**. This is a function that gives the area under the pdf up to the specified value. Given $p(x)$, the cdf is defined as:\n",
    "\n",
    "$$\n",
    "P(y) = Pr(x \\leq y) = \\int_{-\\infty}^y p(x) dx\n",
    "$$\n",
    "\n",
    "The conditional cdf is defined in terms of the conditional pdf:\n",
    "\n",
    "$$\n",
    "P(y|z)  = Pr(x \\leq y|z) = \\int_{-\\infty}^y p(x|z) dx\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exercise 3\n",
    "\n",
    "Edit the equations in the next cell to define $p(FA)$ and $p(H)$ in terms of the **cdfs** for $p(r|s)$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "p(FA) & =  \\\\\n",
    "p(H) & = \n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The Python implementation of the normal distirbution has a method `cdf()` that will evaluate the **cdf** for you.\n",
    "\n",
    "This allows us to generate a plot of $p(H)$ versus $p(FA)$ for different values of the criterion. This is called the ROC curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a range of criteria to test\n",
    "gamma = np.arange(5, 20)\n",
    "# probability of false alarm\n",
    "pfa = 1 - pr_noise.cdf(gamma)\n",
    "# probability of hit\n",
    "phit = 1 - pr_signal.cdf(gamma)\n",
    "plt.plot(pfa, phit, 'o', label=r\"$\\sigma^2 = {}$\".format(5))\n",
    "plt.legend(loc='lower right')\n",
    "plt.xlabel(\"P(FA)\");\n",
    "plt.ylabel(\"P(H)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that there is an inherent tradeoff. You can only increase the hit rate by also increasing the false alarm rate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exercise 4\n",
    "\n",
    "Go back to the original definitions of $p(r|s)$. What do you think will happen if you increase or decrease $\\sigma^2$?\n",
    "\n",
    "Write a function to calculate the ROC as a function of $\\sigma^2$. Plot ROCs (using the same array of $\\gamma$ values as above) for $\\sigma^2 = 2, 10, 20, 100$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Excursion\n",
    "\n",
    "If we assume that the conditional response distributions are normally distributed with equal variance, then we can summarize the relationship between them using a single number, $d'$ (d-prime). This is defined as the difference between the means divided by the standard deviation:\n",
    "\n",
    "$$\n",
    "d' = \\frac{\\beta_1}{\\sigma} = \\frac{\\mu_s - \\mu_n}{\\sigma}\n",
    "$$\n",
    "\n",
    "Furthermore, because the difference between the two means is itself a normal distribution, we can estimate $d'$ from the false alarm and hit rate at a single criterion value:\n",
    "\n",
    "$$\n",
    "d' = z(H) - z(FA)\n",
    "$$\n",
    "\n",
    "where $z(H)$ is the z-score of the hit rate and $z(FA)$ is the z-score of the false alarm rate. The z-score is simply the inverse of the normal cdf. We can calculate it in Python using the `ppf` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pfa = 1 - pr_noise.cdf(10)\n",
    "phit = 1 - pr_signal.cdf(10)\n",
    "\n",
    "std_normal = dists.normal()\n",
    "d_prime = std_normal.ppf(phit) - std_normal.ppf(pfa)\n",
    "print(\"d' =\", d_prime)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 5\n",
    "\n",
    "1. With $\\sigma^2 = 5$, verify that $d'$ is the same for any value of $\\gamma$ in the support.\n",
    "2. Verify that this empirical calculation is correct given the value defined for $\\beta_1$ above.\n",
    "3. What is $d'$ when $\\sigma^2$ is increased to 10?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Information Theory\n",
    "\n",
    "Signal detection theory is powerful, but it's not trival to apply to **discrimination** problems where there are more than two stimuli.\n",
    "\n",
    "Another way of summarizing decoding models is through information theory, which builds on the foundational work of Claude Shannon.\n",
    "\n",
    "We need to understand two key concepts: **entropy** and **information**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Entropy\n",
    "\n",
    "The idea of uncertainty is inherent to probability theory. If some event has a probability that's less than 1.0, that means that we're uncertain about whether it will happen.\n",
    "\n",
    "A probability distribution quantifies our uncertainty about all the values a particular random variable can take on.\n",
    "\n",
    "One way of summarizing the uncertainty of the whole distribution is by its **entropy**, which is defined as follows:\n",
    "\n",
    "$$\n",
    "H(R) = - \\sum_{r \\in R} p(r) \\log_2 p(r)\n",
    "$$\n",
    "\n",
    "For continuous distributions, the entropy is defined by an integral over the support of the pdf:\n",
    "\n",
    "$$\n",
    "H(R) = - \\int_{R} p(r) \\log_2 p(r) dr\n",
    "$$\n",
    "\n",
    "If, as above, a base-2 logarithm is used, then entropy has units of **bits**. You can interpret the entropy as the number of digital bits needed to represent the possible values of the distribution.\n",
    "\n",
    "For example, a fair coin has equal probability of being heads or tails, which gives it an entropy of 1 bit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coin = dists.bernoulli(p=0.5)\n",
    "coin.entropy() / np.log(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "When the natural log is used, the units are called **nats**. The `scipy` distribution functions return entropy in nats, so we have to convert to bits by dividing by $\\log 2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exercise 6\n",
    "\n",
    "In constrast, if we have a coin that always lands heads, then there is no uncertainty about the outcome, and we need 0 bits to store it.\n",
    "\n",
    "Thus, we can infer that the entropy of the Bernoulli distribution depends on its parameter.\n",
    "\n",
    "Plot the entropy of the Bernoulli distribution as a function of its parameter $p$ (the probability that the value will be 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Information\n",
    "\n",
    "Information has many colloquial meanings, but in information theory it is defined as a reduction in uncertainty.\n",
    "\n",
    "**Mutual information** is the reduction in uncertainty about the value of one random variable when you know the value of another.\n",
    "\n",
    "$$\n",
    "I(S;R) = \\sum_{s \\in S, r \\in R} p(s, r) \\log_2 \\frac{p(s,r)}{p(s) p(r)}\n",
    "$$\n",
    "\n",
    "To understand this formula, think about the definition of a joint probability distribution. If the two variables are independent, that means that their joint distribution is just the product of the marginal distributions, \n",
    "$p(s, r) = p(s) p(r)$.\n",
    "\n",
    "In this case, the last term in the sum above is\n",
    "\n",
    "$$\n",
    "\\log_2 \\frac{p(s,r)}{p(s) p(r)} = \\log_2 \\frac{p(s)p(r)}{p(s)p(r)} = \\log_2 1 = 0\n",
    "$$\n",
    "\n",
    "In contrast, if $S$ and $R$ are not mutually independent, then $p(s,r) \\neq p(s)p(r)$ and $I(S; R) > 0$. MI is always non-negative."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 7\n",
    "\n",
    "For the signal detection problem, with $\\beta_0 = 10$, $\\beta_1 = 5$, $\\sigma^2 = 5$, and $p(s) = 0.5$, calculate the mutual information between the signal $S$ and the response $R$.\n",
    "\n",
    "Calculating the joint probability distributions can be a little tricky, so it helps to factor them out:\n",
    "\n",
    "\\begin{align}\n",
    "I(S;R) & = \\sum_{s \\in S, r \\in R} p(s, r) \\log_2 \\frac{p(s,r)}{p(s) p(r)} \\\\\n",
    " & = \\sum_{s \\in S, r \\in R} p(r|s)p(s) \\log_2 \\frac{p(r|s)}{p(r)} \\\\\n",
    " & = \\sum_{s \\in S} p(s) \\sum_{r \\in R} p(r|s) \\log_2 \\frac{p(r|s)}{p(r)}\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta  = [10, 5]\n",
    "sigma = np.sqrt(5)\n",
    "r = np.arange(0, 30, 0.1)\n",
    "pr_noise = dists.normal(mean=beta[0], std=sigma)\n",
    "pr_signal = dists.normal(mean=beta[0] + beta[1], std=sigma)\n",
    "ps = 0.5\n",
    "# insert your code here to calculate I(S;R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happens when you change the discretization of $𝑟$? Why? What would you need to do to remove this dependence?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Information and entropy\n",
    "\n",
    "Mutual information can also be calculated as a reduction in entropy:\n",
    "\n",
    "\\begin{align}\n",
    "I(S;R) & = H(S) - H(S|R) \\\\\n",
    "I(S;R) & = H(R) - H(R|S)\n",
    "\\end{align}\n",
    "\n",
    "The second term in each of the equations above is a **conditional entropy**, which can be calculated much like the standard entropy, but integrated over all the possible values of the conditioned variable:\n",
    "\n",
    "$$\n",
    "H(R|S) = - \\sum_{s \\in S} p(s) \\sum_{r \\in R} p(r|s) \\log_2 p(r|s)\n",
    "$$\n",
    "\n",
    "In the neural setting, the conditional entropy represents how much noise there is in the response. In other words, if the stimulus is known, how much does the response vary from trial to trial?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "One way of interpreting this relationship is as an *information channel* between the stimulus and the neural response. \n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l9_information_channel.png\" alt=\"entropy and mutual information\" style=\"width: 400px;\"/>\n",
    "\n",
    "The amount of information carried by this channel can be thought of as either the total amount of information the neuron *could* represent in its response distribution $H(R)$, reduced by the variability in the response that is independent of the stimulus, $H(R|S)$, or as the total amount of uncertainty about the stimulus $H(S)$ reduced by the remaining uncertainty about the stimulus when the response is known, $H(S|R)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exercise 8\n",
    "\n",
    "Use one of the entropy formulas to calculate mutual information for the SDT problem, and verify that it agrees with the value you calculated previously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary and parting thoughts\n",
    "\n",
    "Decoding models represent the probability of the stimulus conditional on the response.\n",
    "\n",
    "They can only represent a **best-case scenario**. They don't actually tell us how information in the neural response is intepreted by neurons downstream.\n",
    "\n",
    "MI can be used to identify stimulus regimes where a neuron carries more information.\n",
    "\n",
    "MI can also be used to compare neurons and brain regions.\n",
    "\n",
    "Calculating MI with small amounts of data requires a lot of assumptions."
   ]
  }
 ],
 "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
}