{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "hide_input" ] }, "outputs": [], "source": [ "from IPython.core.display import HTML, Markdown, display\n", "\n", "import numpy.random as npr\n", "import numpy as np\n", "import pandas as pd\n", "import seaborn as sns\n", "import matplotlib.pyplot as plt\n", "import scipy.stats as stats\n", "import statsmodels.formula.api as smf\n", "\n", "import ipywidgets as widgets\n", "\n", "import requests\n", "import zipfile\n", "import os\n", "import shutil\n", "import nibabel\n", "from scipy.stats import gamma\n", "\n", "# Enable plots inside the Jupyter Notebook\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Signal detection simulation -- Part 1" ] }, { "cell_type": "markdown", "metadata": { "tags": [ "popout" ] }, "source": [ "Authored by *Liang Zhli*, *Clay Curtis*, *Brenden Lake* and *Todd Gureckis*.\n", "Justin Gardner's [Signal Detection Tutorial](http://gru.stanford.edu/doku.php/tutorials/sdt) but translated from MATLAB to Python." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We're going to simulate a signal detection experiment and an “ideal observer” (an observer who behaves exactly according to signal detection theory). This is always a useful thing to do when trying to understand a decision model, experimental method or an analysis tool. You get to control and play around with the simulation to see what effect it has on the analysis that comes out.\n", "\n", "\n", "On each trial, our observer sees an element sampled from either the signal present gaussian distribution or the signal absent distribution, which is also gaussian with the same standard deviation. The observer chooses to say “signal present” when the signal they see on that trial is above criterion and “signal absent” otherwise. The picture you \n", "should have in your head is this:\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Question 1
\n", " To check your understanding we will first have you create SDT type figures matching a few example situations. Using the code provided below, please create figures matching the following situations:\n", "
\n", "
• Draw two bell-shaped curves on the same plot (Gaussian distributions) with the same mean, but different standard deviations.
• \n", "
• Next draw two bell-shaped curves (Gaussian distributions) with the same standard deviations, but different means.
• \n", "
• Next draw one signal-detection-theory plot for a subject who has POOR discrimination, and another signal-detection-theory plot for a different subject is has GOOD discrimination.
• \n", "
• On the SDT plots that you just completed, draw a liberal criterion for one subject, and a conservative criterion for the other. Label each of these clearly.
• \n", "
• Describe in a markdown cell which parts of the graphs reflect the criterion, misses, false alarms, hits, and correct rejections.
• \n", "
\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Enter the code for each graph below.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = np.linspace(-10,30,100)\n", "noise=stats.norm.pdf(x,10,3.5)\n", "signalnoise=stats.norm.pdf(x,15,3.5)\n", "\n", "# plot the \n", "plt.plot(x,noise,color='steelblue')\n", "# this fills in the color of the plot\n", "plt.fill_between(x,noise, interpolate=True,facecolor='lightblue',alpha=0.2)\n", "\n", "\n", "plt.plot(x,signalnoise,color='red')\n", "# this fills in the color of the plot\n", "plt.fill_between(x,signalnoise,interpolate=True,facecolor='pink',alpha=0.2)\n", "\n", "plt.ylim(0,0.15)\n", "plt.ylabel(\"Probability\")\n", "plt.xlabel(\"Internal response\")\n", "plt.title(\"Your title here\")\n", "# this line will add a line for the criterion.\n", "#plt.plot([10,10],[0,0.14],color='k')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Make signal_present array" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First let's make a sequence of n = 1000 trials in which 50% randomly are going to have “signal” and 50% are randomly going to have no signal. To do that we want to create an array called signal_present that is True 50% of the time and False 50% of the time." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "signal_present = np.random.rand(1000) > .5" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can make a heatmap to show the binary values:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sns.heatmap(signal_present.reshape(25, 40), linewidths=0, xticklabels=False, yticklabels=False);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Make signal array " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ok, now we want to simulate the **signal** that the observer actually gets to base the decision on. Remember that in applying signal detection to cognitive neuroscience you should be thinking of the signal as a neuron or population of neurons response and that the magnitude of this response (e.g. spikes per second) is monotonically related to the actual signal strength. The signal is corrupted by random noise. In particular, signal trials should come from some gaussian distribution and noise trials should come from another gaussian distribution that differ only in the means. This is an assumption about the process that is termed $iid$ - the signal and noise come from independent identical distributions.\n", "\n", "Ok, let's make a new array from our signal_present array such that on signal present trials (i.e. when signal_present == 1), values are picked from a gaussian distribution with standard deviation of 1 and whose mean is 2. On signal absent trials (i.e. when signal_present == 0), values are picked from a gaussian distribution with standard deviation of 1 (that's the identical part of the $iid$ assumption) and whose mean is 0." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "signal = np.random.normal(0, 1, size=signal_present.size)\n", "signal[signal_present] = np.random.normal(2, 1, size=signal_present.sum())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can look at a heatmap of the signal values:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sns.heatmap(signal.reshape(25, 40), linewidths=0, xticklabels=False, yticklabels=False);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is somewhat helpful, but it doesn't really show us the distribution of values. To do that, let's first put all the data into a pandas data frame." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "df=pd.DataFrame({\"trial\": range(len(signal)), \"signal_present\": signal_present, \"signal\": signal})\n", "df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then let's make a histogram and density plot." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Question 2
\n", " Using seaborn's histplot() create a histogram of the signal column from our data frame.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Question 3
\n", " Using histplot() create a histogram of the signal present data separately from the signal absent data. Hint, you can select only the signal present data by df.signal[df.signal_present]. Second hint: the ~ character means \"not\" or \"opposite\".\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Question 4
\n", " Now confirm numerically that the means and standard deviations of the signal present and signal absent data are approximately what they should be (hint: use the .mean() function in pandas).\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ok, we're good to go. We have created the stimulus that our ideal observer will get to see (signal) and we know which of the trials come from the signal and which from the absent distributions (signal_present)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Make ideal observer responses" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we are going to simulate an *ideal observer* which will behave just as signal detection says. They will choose signal present (response = 1) when the signal they get to see (signal array from above) is greater than their internal criterion and they will choose signal absent (response = 0) when the signal is below their internal criterion.\n", "\n", "Let's start by making the criterion right in between the signal present and absent distributions that we created above. That is, let's set criterion to 0.5 and make an array of responses:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "df['response'] = df.signal > .5" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's make a heatmap like before" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sns.heatmap(df.response.values.reshape(25, 40), linewidths=0, xticklabels=False, yticklabels=False);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, not too informative, but it would help us catch strange cases like response always being True, which might indicate a problem in the code." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculate hits, misses, correct-rejects and false-alarms" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ok, so now we have our experiment dataframe (signal_present) a simulation of the signal it generates in the observer (signal) and the ideal observers responses (response).\n", "\n", "From these you should be able to calculate hits, misses, correct rejects and false alarms." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hit = df.response[df.signal_present]\n", "miss = ~df.response[df.signal_present]\n", "fa = df.response[~df.signal_present]\n", "cr = ~df.response[~df.signal_present]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"Hit rate: {:.2f}\".format(hit.mean()))\n", "print(\"Miss rate: {:.2f}\".format(miss.mean()))\n", "print(\"False alarm rate: {:.2f}\".format(fa.mean()))\n", "print(\"Correct rejection rate: {:.2f}\".format(cr.mean()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculate $d'$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's calculate the exact $d'$, based on how we made the signal and noise distributions. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Question 5
\n", "What should $d'$ be, given how we made the signal and noise distributions? Refer to the chapter reading for the formula for d'.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The picture you should have in mind is:\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ok, now that we know what $d'$ should be in an ideal world, we're ready to think about the approach for computing $d'$ from behavioral data. Here is how we estimate $d'$ from a set of behavioral data,\n", "\n", "$$d' = z(\\textrm{hits}) - z(\\textrm{false alarms}).$$\n", "\n", "So, how do we calculate those $z$'s?\n", "\n", "$z$ means the distance in units of standard deviation from the center of a gaussian distribution such that the area under the curve to the right of that is the proportion hits or false alarms that we measured in the experiment. By convention, center of the distribution is 0, and criterion to the left are positive. The picture you should have for example if $z$ is 1 (giving an area under the curve to the right of that of 0.84 is):\n", "\n", "\n", "\n", "Wee don't have the z-score for the hits directly. Instead, we have the empirical proportion of hits or false alarms. Interestingly those values relate directly to the area under the curve of the signal strength/internal response distributions. Thus we can use the **inverse cumulative density** function or inverse CDF to look up the z-score associated with a particular area under the curve for a standard normal (mean=0, sd=1).\n", "\n", "To get this area you can use the stats.norm.ppf() function which gives the inverse of the cumulative density function. It's the cumulative density function since you are interested in the area under the gaussian probability density function, and it's the inverse since you are going from area back to the $z$ (units of std of the gaussian)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy import stats\n", "dprime = stats.norm.ppf(hit.mean()) - stats.norm.ppf(fa.mean())\n", "print(\"d prime: {:.2f}\".format(dprime))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Stop and think
\n", " How close is the d' you got to the expected value? Why is it different?\n", "
\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Different criterion\n", "\n", "Now let's simulate an observer that doesn't want to miss as much and see if we get a similar $d'$ from the data." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "df['response'] = signal > 0.2" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hit = df.response[df.signal_present]\n", "fa = df.response[~df.signal_present]\n", "dprime = stats.norm.ppf(hit.mean()) - stats.norm.ppf(fa.mean())\n", "print(\"d prime: {:.2f}\".format(dprime))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Stop and think
\n", " Did d' change? Try a few other threshold rules. Is it always the same value? Is that what you expected? If not why do you think it changes?\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Bias or criterion\n", "\n", "Another value of interest in a signal detection analysis is the bias or criterion. Response bias in a yes/no task is often quantified with $\\beta$. Use of this measure assumes that responses are based on a likelihood ratio. However, a somewhat more popular alternative is the measure/statistic known as the criterion, $c$. $c$ is defined via the distance from a neutral point, where neither response is favored. If the criterion is located at this neutral point, c has a value of 0. Deviations from the neutral point are measured in standard deviation units. Negative values of c signify a bias toward responding yes (the criterion lies to the left of the neutral point), whereas positive values signify a bias toward the no response (the criterion lies to the right of the neutral point)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "c = -(stats.norm.ppf(hit.mean()) + stats.norm.ppf(fa.mean()))/2.0\n", "print(\"c: {:.2f}\".format(c))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Stop and think
\n", " What value of the response threshold above will make c closer to or equal to zero? Try adjusting it until you find the value. Think about why that is.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## ROC calculation\n", "\n", "Now, let's make an ROC curve. In some ways a ROC curve is a more theoretical object because it can be difficult to construct it for a simple detection experiment. In fact, for a single run of simulated data we have been dealing with so far, it is not possible to plot the entire ROC curve because we don't know exactly how the subject might perform using a threshold different than the one we saw them adopt.\n", "\n", "However, if we can draw simulated samples from the signal and noise distributions, we can draw the expected curve by applying different values of the threshold ourselves. In values of the spike counts for the neurons preferred direction (signal) and the neurons anti-preferred direction (noise, what they called the anti-neuron). Will just do it by pulling actual values from a signal distribution that is gaussian with mean of 2 and standard deviation of 1, just like we did above. For the noise, let's do mean of 0 and standard deviation of 1. Make two arrays of length n = 1000, one for signal and noise using the techniques from above:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n = 1000\n", "signal = np.random.normal(2, 1, n)\n", "noise = np.random.normal(0, 1, n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check that the procedure worked by computing the means and standard deviations of signal and noise and making sure that they are what you expected (always good to double check everything when programming to make sure that things are working the way you expected)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"Signal mean +/- std: {:.2f} ({:.2f})\".format(signal.mean(), signal.std()))\n", "print(\"Noise mean +/- std: {:.2f} ({:.2f})\".format(noise.mean(), noise.std()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It's also helpful to make plots like we did above:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sns.histplot(signal, color=\".2\", label=\"Signal\")\n", "sns.histplot(noise, color=\"darkred\", label=\"Noise\")\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ok, how do we compute the ROC curve? What's on those two axes. Remember?\n", "\n", "\n", "\n", "Ok, so that is computed for various criterions, right? So, if you think you know how, try to compute a curve from the signal and noise distributions.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First we need to create a large number of different criterion values. How about 100 different values between 10 and -10 using np.linspace()." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "criterions = np.linspace(10, -10, 100)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now for each setting of the criterion we can compute the hit and false alarm rate." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hit_rates = [(signal > c).mean() for c in criterions]\n", "fa_rates = [(noise > c).mean() for c in criterions]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ok, now plot it:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f, ax = plt.subplots(figsize=(6, 6))\n", "ax.plot(fa_rates, hit_rates)\n", "ax.plot([0, 1], [0, 1], ls=\"--\", c=\".5\")\n", "ax.set(xlabel=\"FA Rate\", ylabel=\"Hit Rate\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, one interesting fact is that the area under the ROC is equal to performance on a two alternative forced choice task (2AFC). So, calculate the area under the curve and see what you get." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.integrate import trapz\n", "auc = trapz(hit_rates, fa_rates) # this does the integral under the ROC curve\n", "print(\"Area under ROC curve: {:.2f}\".format(auc))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What did you get? Is it right? Well, one way to check is to see if a simulation of 2AFC will give you the same (or numerically similar) values. The experimental logic says that if signal is greater than noise then the ideal observer will get the 2AFC right, otherwise the observer will choose the noise interval and get it wrong. So, that's pretty easy to simulate. Try it!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "afc_perf = (signal > noise).mean()\n", "print(\"2AFC performance: {:.2f}\".format(afc_perf))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Question 6
\n", " In what we just described, I suggested that you can't compute the ROC curve unless you have samples from the underlying internal response distributions (e.g., neural spike counts). Is this correct? Think about what the meaning of the ROC curve is. Can you think of a way to assess it using an experiment where the primary data are just the false alarms and hits for different conditions? \n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " In the example reported above we computed the hits, false alarms, misses, and correct rejections, along with the $d'$ and $c$ values for a case where the signal absent distribution had mean 0 and standard deviation 1, while the signal present distribution had a mean 2 and standard deviation 1. Recompute these statistics for a case where the mean of the signal distirbution is 1 and the criterion is 0.4. What effect does this have on d'? \n", "