Day 6 - Permutation tests, multiple comparisons.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple Comparisons in Statistical Tests\n",
    "\n",
    "While our goal with this course is not to teach you statistics, we do need to take a slight detour to address some common statistical issues that come up in ECoG analyses. First, we'll touch on the problem of **multiple comparisons** ( https://en.wikipedia.org/wiki/Multiple_comparisons_problem ). This refers to the idea that, if you run many statistical tests, you're more likely to get a false positive with increasing numbers of tests. For example, if a subject has 100 electrodes, and we'd like to test whether there's an increase in power at each electrode, some will -- purely by chance -- show an increase in power even if it's not meaningfully correlated with memory. \n",
    "\n",
    "The easiest way to correct for multiple comparisons is to set a more stringent threshold for determining significance. Typically, these thresholds are set as <i>p</i> values, which you're likely quite familiar with. The <i>p</i> value is the proportion of test statistic values that, under the null hypothesis, exceed the true test statistic. If too many null values exceed the true value, it's hard to say that the true value is \"significant.\" But if you set the <i>p</i> value to something very small -- as in, almost no null values exceed the true value -- it becomes quite difficult to find a \"significant\" result. \n",
    "\n",
    "<center>\n",
    "<img src=\"https://med.stanford.edu/news/all-news/2016/03/misleading-p-values-showing-up-more-often-in-journals/_jcr_content/main/image.img.620.high.jpg\" width=500>\n",
    "</center>\n",
    "\n",
    "Indeed, the **Bonferonni correction** is a common and simple way to account for the multiple comparisons problem. The basic idea is to take the <i>P</i> value you'd ordinarily use, and divide it by the number of tests you're running. So if usually you'd want a P<0.05 to declare significance, but you run 100 tests, now you'd need a P<0.0005 for significance. This is a very conservative way to control for multiple comparisons (strictly, Bonferonni correction controls for the \"familywise error\" rate). \n",
    "\n",
    "If a test survives Bonferonni correction (and assuming other assumptions about the test were also met), that's a strong indication that there's something real there (i.e. you're not looking at a false positive), and you can make some conclusions about it. Bonferonni falls under the category of **family-wise error rate (FWER) control**, or the probably that *any* test leads to a Type I (false positive) error. \n",
    "\n",
    "You'll also see \"False Discovery Rate\" controlling procedures, or **FDR correction**, in the ECoG literature. These are less conservative tests that let you screen for effects in large amounts of data, but do not offer the same kind of strict control that Bonferonni does (therefore, be more cautious in your interpretation of the \"significant\" effects.). A common FDR procedure is called \"Benjamini-Hochberg,\" though there are several algorithms with unique properties. FDR correction limits the total *proportion* of Type I errors across all tests.\n",
    "\n",
    "Many types of correction can be applied using the stats toolbox within MNE. Let's explore a few of these now, and use them on real data. \n",
    "\n",
    "<center>\n",
    "<img src=\"https://github.com/esolomon/PythonBootcamp2019/blob/master/figures/fwer.png?raw=true\" width=300>\n",
    "<center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, let's load up some data.\n",
    "import numpy as np\n",
    "import os\n",
    "# With dimensions [frequencies, electrodes, time bins].\n",
    "pvals = np.load(os.path.join('npy_files','R1001P_pvals.npy'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Count')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVoAAADUCAYAAADdnOtqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAARY0lEQVR4nO3debBkZXnH8e9PFESJuDAgAc3FYBmMiSjjnjIajXsQLUWNCyBqaSxFTVQmVEWtiqVW3JdoUBQ1KIMTiCgIuOFSARQQRcUFYVBcANFyj4o++eOcG3quPTPdM/123+75fqq6bp/3nNPnuaemfnPu6e73SVUhSWrnBrMuQJIWnUErSY0ZtJLUmEErSY0ZtJLUmEErSY3dcNYFTNsee+xRS0tLsy5D0oK54IILflhVa4at2+GCdmlpifPPP3/WZUhaMEmu2Nw6bx1IUmMGrSQ1ZtBKUmMGrSQ1ZtBKUmPNgjbJO5NcneTLQ9b9U5JKssfA2Loklyb5epIHD4wflOTift0bk6Qf3yXJ+n78vCRLrX4XSdoeLa9ojwcesnIwyW2AvwW+PTB2R+DxwJ/3+/x7kp361W8FngHcvn8sv+aRwI+ran/gdcCrmvwWWmjrX3Y061929KzLGNkpr7mQU15z4UyOfcWTn8IVT37K2PsdccYRHHHGEQ0qmh/NgraqPg38aMiq1wEvAgYnwn0kcGJV/bqqLgcuBe6eZG/gZlV1TnUT574HOGRgn3f3zzcAD1i+2pWk1WSq92iTHAx8t6q+uGLVPsB3Bpav7Mf26Z+vHN9kn6q6DvgJcKsGZUvSdpnaN8OS3AQ4BnjQsNVDxmoL41vaZ9ixn0F3+4Hb3va2W61VkiZpmle0fwrsB3wxyUZgX+DCJLemu1K9zcC2+wLf68f3HTLO4D5JbgjszvBbFVTVsVW1tqrWrlkz9KvIktTM1IK2qi6uqj2raqmqluiC8q5V9QPgVODx/ScJ9qN70+tzVfV94GdJ7tnff30K8MH+JU8FDuufPwb4RNkATdIq1PLjXe8HzgHukOTKJEdubtuq+gpwEvBV4Azg2VX1u371s4B30L1B9i3gI/34ccCtklwKvACYn7eOJe1Qmt2jraonbGX90orllwMvH7Ld+cCdhoz/L/DY7atS0qJYOvq0ib3Wxlc+fGKvBX4zTJKaM2glqTGDVpIaM2glqTGDVpIaM2glqTGDVpIaM2glqTGDVpIaM2glqTGDVpIaM2glqTGDVpIam2oX3CT/luRrSb6U5JQkNx9YZxdcSQtp2l1wPwrcqar+EvgGsA7sgitpsU21C25VndU3UgQ4l+vb1NgFV9LCmuU92qdyfbcEu+BKWlgzCdokxwDXAScsDw3ZbKJdcJOcn+T8a665ZtxyJWm7TD1okxwGPAJ44kAzRbvgSlpYUw3aJA8BXgwcXFW/HFhlF1xJC6tZc8a+C+79gD2SXAm8hO5TBrsAH+3ftzq3qp5ZVV9JstwF9zr+sAvu8cCudPd0B7vgvrfvgvsjuk8tSNKqM+0uuMdtYXu74EpaSH4zTJIaM2glqTGDVpIaM2glqTGDVpIaM2glqTGDVpIaM2glqTGDVpIaM2glqTGDVpIaM2glqTGDVpIam3YX3Fsm+WiSb/Y/bzGwzi64khbStLvgHg18vKpuD3y8X7YLrqSFNtUuuGzaufbdbNrR1i64khbStO/R7tW3p6H/uWc/3rQLrs0ZJc3SankzrGkXXJszSpqlaQftVf3tAPqfV/fjTbvgStIsTTtoBzvXHsamHW3tgitpIU27C+4rgZOSHAl8m765ol1wJS2yaXfBBXjAZra3C66khbRa3gyTpIVl0EpSYwatJDVm0EpSYwatJDVm0EpSYwatJDVm0EpSYyMFbZL7jDImSfpDo17RvmnEMUnSClv8Cm6SewH3BtYkecHAqpsBOw3fS5I0aGtzHewM7NZv90cD4z+lmzFLkrQVWwzaqvoU8Kkkx1fVFZM6aJLnA0+jm6j7YuAI4CbAemAJ2AgcWlU/7rdfR9cj7HfAc6vqzH78IK6f2et04CinSpS02ox6j3aXJMcmOSvJJ5Yf23LAJPsAzwXWVtWd6G5BPJ7JNm6UpFVj1GkSPwC8DXgH3VXlJI67a5Lf0l3Jfg9YRzd/LXRNF88GXsxA40bg8n7+2bsn2UjfuBEgyXLjxo8gSavIqEF7XVW9dRIHrKrvJnk13cTfvwLOqqqzkmzSuDHJYOPGcwdeYrlB42/ZfONGSVo1Rr118KEk/5Bk7yS3XH5sywGT3ILuKnU/4I+BmyZ50pZ2GTK2tcaNK49pF1xJMzPqFe1yb64XDowVcLttOOYDgcur6hqAJCfTfYTsqiR791ez29u4cRNVdSxwLMDatWt9s0zSVI10RVtV+w15bEvIQnfL4J5JbtI3XHwAcAmTbdwoSavGSFe0SZ4ybLyq3jPuAavqvCQbgAvpGjF+ge5qczcm17hRklaNUW8d3G3g+Y3prkIvBMYOWoCqegldV9xBv2ZCjRslaTUZKWir6jmDy0l2B97bpCJJWjDbOk3iL+nulUqStmLUe7Qf4vqPTu0EHACc1KooSVoko96jffXA8+uAK6rqys1tLEm63qgf7/oU8DW6GbxuAfymZVGStEhG7bBwKPA5uo9cHQqcl8RpEiVpBKPeOjgGuFtVXQ2QZA3wMWBDq8IkaVGM+qmDGyyHbO/aMfaVpB3aqFe0ZyQ5E3h/v/w4uom2JUlbsbWeYfsDe1XVC5M8GvgrulmzzgFOmEJ9kjT3tvbn/+uBnwFU1clV9YKqej7d1ezrWxcnSYtga0G7VFVfWjnYzzGw1KQiSVowWwvaG29h3a6TLESSFtXWgvbzSZ6+crCfyvCCbT1okpsn2ZDka0kuSXKvvmvDR5N8s/95i4Ht1yW5NMnXkzx4YPygJBf3697Yz0srSavK1oL2ecARSc5O8pr+8Sm6VuFHbcdx3wCcUVV/BtyZbuJvu+BKWkhb/NRBVV0F3DvJ/bl+3tfTqmqbWo0DJLkZcF/g8P4YvwF+k+SR2AVX0gIadT7aTwKfnNAxbwdcA7wryZ3pbkEcRfcxMrvgSlo4s/h21w2BuwJvraq7AL+gv02wGXbBlTTXZhG0VwJXVtV5/fIGuuC9qu9+S4suuFW1tqrWrlmzZmK/iCSNYupBW1U/AL6T5A790APoGi/aBVfSQhp1roNJew5wQpKdgcuAI+hC3y64khbOTIK2qi4C1g5ZZRdcSQtnVle0knZwS0efNusSpsY5ZSWpMYNWkhozaCWpMYNWkhozaCWpMYNWkhozaCWpMYNWkhozaCWpMYNWkhozaCWpMYNWkhqbWdAm2SnJF5J8uF+2C66khTTLK9qj6LrfLrMLrqSFNJOgTbIv8HDgHQPDj6Trfkv/85CB8ROr6tdVdTmw3AV3b/ouuFVVwHsG9pGkVWNWV7SvB14E/H5gbJMuuMBgF9zvDGy33O12H+yCK2kOTD1okzwCuLqqLhh1lyFjdsGVNDdmcUV7H+DgJBuBE4G/SfKf2AVX0oKaRRfcdVW1b1Ut0b3J9YmqehJ2wZW0oFZTz7BXYhdcSQtopkFbVWcDZ/fPr8UuuJIWkN8Mk6TGDFpJasyglaTGDFpJasyglaTGDFpJasyglaTGDFpJasyglaTGDFpJasyglaTGDFpJasyglaTGpj57V5Lb0PX3ujVdK5tjq+oNSW4JrAeWgI3AoVX1436fdcCRwO+A51bVmf34QVw/TeLpwFF9/zBJDZx72bX89dGnjbXPrre9FoClMfdbJLO4or0O+MeqOgC4J/DsvtOtXXAlLaRZdFj4flVd2D//GV3L8X2wC66kBTXTe7RJloC7AOfRsAuuzRklzdLMgjbJbsB/Ac+rqp9uadMhY2N1wbU5o6RZmknQJrkRXcieUFUn98PNuuBK0ixNPWj7jrXHAZdU1WsHVtkFV9JCmkVzxvsATwYuTnJRP/bP2AVX0oKaetBW1WcZfn8V7IIraQH5zTBJasyglaTGDFpJasyglaTGDFpJasyglaTGZvE5WklTMsmpCV912bUTe60djVe0ktSYQStJjRm0ktSYQStJjflmmLTKbOkNrMf9bGcAnr8D99+aR3N/RZvkIUm+nuTSJEfPuh5JWmmug7Zv0vgW4KHAHYEn9M0cJWnVmPdbB3cHLq2qywCSnEjXzPGrM61KO5QduY22RjPvQTusceM9ZlSL5shyOD7q+z8C4MWGpRpK16l7PiV5LPDgqnpav/xk4O5V9ZwV2z0DeEa/eAfg62MeanfgJ2OuXzm2peVhz5d/7gH8cAr1bkuNk6p3W2tetHM8rLZJ1butNS/yOR5W5/bU+ydVNbz7a1XN7QO4F3DmwPI6YF2D4xw77vqVY1taHvZ84Of506h3W2qcVL2e4z9c9hy3P8ebqXO7zvHmHnP9ZhjweeD2SfZLsjPweLpmjpP2oW1Yv3JsS8vDnm/tmOPWM8r6cWscfL499Y6y/45wjgeXPcftz/Hg80md46Hm+tYBQJKHAa8HdgLeWV1/sYWR5PyqWjvrOkY1b/XC/NU8b/XC/NU86Xrn/c0wqup04PRZ19HQsbMuYEzzVi/MX83zVi/MX80TrXfur2glabWb93u0krTqGbSS1JhBK0mNGbRzLMn9knwmyduS3G/W9YwiyU2TXJDkEbOuZRRJDujP74Ykz5p1PVuT5JAkb0/ywSQPmnU9o0hyuyTHJdkw61o2p/93++7+3D5x3P0N2hlJ8s4kVyf58orxcWYjK+DnwI3pvn7czITqBXgxcFKbKjc1iZqr6pKqeiZwKND040kTqve/q+rpwOHA4xqWu1zbJGq+rKqObFvpHxqz9kcDG/pze/DYB5vktx98jPWtl/sCdwW+PDC2E/At4HbAzsAX6WYl+wvgwyseewI36PfbCzhhDup9IN2XSg4HHjEP57jf52Dgf4C/n4d6+/1eA9x1Xs5xv9+G1vVuR+3rgAP7bd437rHm/nO086qqPp1kacXw0NnIquoVwJb+1P4xsEuLOpdNot4k9wduSvcP91dJTq+q36/mmvvXORU4NclpwPtWc71JArwS+EhVXdiq1mUT/nc8VePUTvcX477ARWzDnQCDdnUZazayJI8GHgzcHHhz29KGGqveqjoGIMnhwA9bhuwWjHuO70f3Z+MuzOaLMePOUPccur8cdk+yf1W9rWVxmzHuOb4V8HLgLknW9YE8K5ur/Y3Am5M8nG34mq5Bu7pkyNhmv1FSVScDJ7crZ6vGqvf/N6g6fvKljGzcc3w2cHarYkYwbr1vpAuFWRq35muBZ7YrZyxDa6+qXwBHbOuL+mbY6nIlcJuB5X2B782ollHMW70wfzXPW70wnzUva1K7Qbu6TGs2skmZt3ph/mqet3phPmte1qb2ab7L52OTdzzfD3wf+C3d/6JH9uMPA75B987nMbOuc17rncea563eea15FrU7qYwkNeatA0lqzKCVpMYMWklqzKCVpMYMWklqzKCVpMYMWi2kJL9LclGSLyf5QJKbTOA1D08yizklNOcMWi2qX1XVgVV1J+A3rJ7v0msHZNBqR/AZYP/BgSQ3SLIxyc0Hxi5NsleSv0tyXpIvJPlYkr1WvmCS45M8ZmD55wPPX5jk80m+lORljX4nzRGDVgstyQ2BhwIXD45XN0XjB4FH9dvdA9hYVVcBnwXuWVV3AU4EXjTG8R4E3J5uXtMDgYOS3HcCv4rmmNMkalHtmuSi/vlngOOGbLMe+BfgXXSTh6zvx/cF1ifZm26W/cvHOO6D+scX+uXd6IL302NVr4Vi0GpR/aqqDhwcSPJs4On94sOAc4D9k6wBDgH+tV/3JuC1VXVqP/H3S4e8/nX0fxH2XQ12Xj4M8Iqq+o/J/Sqad9460A6jqt7Sv0F2YFV9r7oZlU4BXgtcUt0E1AC7A9/tnx+2mZfbCBzUP38kcKP++ZnAU5PsBpBknyR7TvhX0ZzxilY7uvV0c5AePjD2UuADSb4LnAvsN2S/twMfTPI54OPALwCq6qwkBwDndBe6/Bx4EnB1o/o1B5wmUZIa89aBJDVm0EpSYwatJDVm0EpSYwatJDVm0EpSYwatJDVm0EpSY/8Hy8ApZSelmI4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's take a look at the distribution of these p-values.\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Plot a histogram of pvalues on a log scale.\n",
    "plt.figure(figsize=(5, 3))\n",
    "ax = plt.subplot(111)\n",
    "_, bins, _ = plt.hist(pvals.ravel(), bins=13, alpha=0)\n",
    "logbins = np.logspace(np.log10(bins[0]),np.log10(bins[-1]),len(bins))\n",
    "plt.hist(pvals.ravel(), bins=logbins, color='C0')\n",
    "plt.vlines([0.1, 0.05, 0.01, 0.001], ymin=ax.get_ylim()[0], ymax=ax.get_ylim()[1], \n",
    "       colors=('C2', 'C3', 'C4', 'C5'),)\n",
    "plt.xscale('log')\n",
    "plt.xlabel('P-value')\n",
    "plt.ylabel('Count')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's clear that there aren't many P-values less than 0.001. Would any of these survive Bonferonni correction?\n",
    "\n",
    "(Perhaps testing individual P-values across all timepoints, frequencies, and electrodes isn't the best approach. What if we restricted ourselves to running tests at the electrode level, for a given time window of interest?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reload the raw power values so we can re-run some tests.\n",
    "rec = np.load(os.path.join('npy_files', 'R1001P_recalled.npy'))\n",
    "pows = np.load(os.path.join('npy_files', '/home1/esolo/bootcamp2019/R1001P_pows.npy'))\n",
    "\n",
    "# pows are dimensions [freqs, events, electrodes, time bins].\n",
    "\n",
    "# Let's average powers into across the 400-700 ms window.\n",
    "# This begins after ERP effects die out,\n",
    "# but not so long as to include irrelevant activity.\n",
    "# Furthermore, let's only analyze the theta frequencies,\n",
    "# which were the first 5 (trust me).\n",
    "avg_pow = np.mean(np.mean(pows[:5, :, :, 4:7], 0), -1)\n",
    "\n",
    "# avg_pow should be an events X electrodes array. Excellent!\n",
    "# Let's do our 2-sample T-test on these averaged values\n",
    "from scipy.stats import ttest_ind\n",
    "ts, ps = ttest_ind(avg_pow[rec==1, :], avg_pow[rec==0, :], axis=0, equal_var=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[False False False False False False  True  True False False False False\n",
      "  True False  True False False  True  True False  True False False False\n",
      " False False False False False False False False False False False False\n",
      "  True False False False  True False False False False False False False\n",
      " False False False False False False  True False False False False  True\n",
      "  True False  True False False False False False False False False False]\n"
     ]
    }
   ],
   "source": [
    "# If we did an uncorrected test, what would we find?\n",
    "sig_elecs = ps<0.05\n",
    "\n",
    "print(ps<0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[False False False False False False False False False False False False\n",
      " False False False False False False False False False False False False\n",
      " False False False False False False False False False False False False\n",
      " False False False False False False False False False False False False\n",
      " False False False False False False False False False False False False\n",
      " False False False False False False False False False False False False]\n"
     ]
    }
   ],
   "source": [
    "# What about Bonferonni correction?\n",
    "\n",
    "from mne.stats import bonferroni_correction\n",
    "\n",
    "reject, pval_cor = bonferroni_correction(ps, 0.05)\n",
    "\n",
    "print(reject)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No electrodes survive Bonferroni correction! Is it too conservative? Remember that Bonferonni assume **no statistical dependency between tests.** What about FDR correction (Benjamini-Hochberg)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[False False False False False False  True False False False False False\n",
      "  True False  True False False False  True False False False False False\n",
      " False False False False False False False False False False False False\n",
      "  True False False False False False False False False False False False\n",
      " False False False False False False False False False False False  True\n",
      " False False False False False False False False False False False False]\n"
     ]
    }
   ],
   "source": [
    "from mne.stats import fdr_correction\n",
    "\n",
    "reject, pval_cor = fdr_correction(ps, 0.05)  #by default, this implements Benjamini-Hochberg with a 0.05 threshold. \n",
    "\n",
    "print(reject)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A few electrodes are back! With this threshold, we would expect, on average, 5% of the significant tests to result in a Type I error. \n",
    "\n",
    "The significant electrodes should not be over-interpreted, since we used a less-conservative test. Some still might be false positives. But they may carry interesting signal, and combined with other information, these signals could teach us something about the brain. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Permutation testing\n",
    "\n",
    "The statistical tests we have discussed so far are **parametric**, meaning that they rely on models of data established by a finite number of parameters. For instance, the t-test models our data as normal distributions, with fixed means and standard deviations. These tests are often powerful and extremely fast, but they do have their drawbacks: Namely, they rely on a number of assumptions about the underlying data, including their distributional forms and **independence**. Many parametric tests assume that datapoints or distributions are uncorrelated or unrelated to each other; if this isn't actually true, you can end up with a false positive!\n",
    "\n",
    "There are several **nonparametric** tests that do not make assumptions about the distributional form of your data. Today, we're going to implement a general approach to nonparametric testing, called a **permutation test**. Permutation tests (also called randomization tests) make no assumptions about the data; instead, they literally create the null distribution against which to compare your test statistic, instead of modeling it based on parameters derived from your data. In fact, there are many ways in which the general idea of resampling statistics are applied throughout data science -- learning about permutation tests is a good way to get some initial exposure to these methods. \n",
    "\n",
    "How to do a permutation test: \n",
    "\n",
    "1. **Decide what you want to test**: Are you looking for the difference in means between two distributions? A difference in their variance? Whether the sum of one distribution is significantly greater than another? A permutation test can, in theory, handle anything you throw at it. \n",
    "2. **Compute your \"true\" statistic.** For whatever you decided in step (1), compute that number. It could be a difference score, a t-statistic, an F-statistic, or anything else. \n",
    "3. **Randomly shuffle your data.** For example, if you're testing something between two distributions, randomly swap data between those distributions, but <i>maintain the original size of the distributions</i>. Remember, a permnutation test is about keeping everything exactly the same, except for the data values themselves. \n",
    "4. For each random shuffle in step (3), **recompute your test statistic**. \n",
    "5. **Repeat steps 3-4 many times.** There's no strict rule on exactly how many times, and the more shuffles you do, the precise P-value you can get. But more shuffles means more compute time, and it may quickly become intractable. 1000 shuffles is common and a good rule-of-thumb, but typically you'll need at least several hundred. \n",
    "6. **Compare true statistic to null distribution.** Typically this takes the form of computing a p-value or a z-score.\n",
    "\n",
    "**WARNING: ONLY SHUFFLE OVER THE DIMENSION THAT REPRESENTS YOUR VARIABLE OF INTEREST. DO NOT SHUFFLE ALL OF YOUR DATA -- PERMUTATION TESTS MUST PRESERVE THE ORIGINAL STRUCTURE OF THE DATA AS MUCH AS POSSIBLE.**\n",
    "\n",
    "#### Example\n",
    "\n",
    "Let's try this out on our power data from R1001P. Yesterday, we ran a 2-sample (or \"independent\") t-test to ask whether there was a significant difference in power between remembered and not-remembered conditions. But in doing so, we violated an assumption of the test: The two distributions are not actually independent! Recalled and not-recalled trials can occur right after one another, and EEG data is autocorrelated in time. While, in practice, this violation might not be a huge deal, it was technically inappropriate to use this test for statistical inference. A permutation test allows us to keep everything above-board. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reload the raw power values and recalled array.\n",
    "rec = np.load(os.path.join('npy_files', 'R1001P_recalled.npy'))\n",
    "pows = np.load(os.path.join('npy_files', 'R1001P_pows.npy'))\n",
    "\n",
    "# pows are dimensions [freqs, events, electrodes, time bins].\n",
    "\n",
    "# Let's average powers into across the 400-700 ms window.\n",
    "# This begins after ERP effects die out,\n",
    "# but not so long as to include irrelevant activity.\n",
    "# Furthermore, let's only analyze the theta frequencies,\n",
    "# which were the first 5 (trust me).\n",
    "avg_pow = np.mean(np.mean(pows[:5, :, :, 4:11], 0), -1)\n",
    "\n",
    "# Step 1: We'd like to test the difference in mean theta power\n",
    "# between rec and not-rec distributions.\n",
    "# Step 2: Compute the true difference in means.\n",
    "\n",
    "diff_mu = np.mean(avg_pow[rec==1, :], 0)-np.mean(avg_pow[rec==0, :], 0)\n",
    "\n",
    "# Step 3, 4, 5: Randomly shuffle data between rec and not-rec groups,\n",
    "# and recompute test statistic (x 1000)\n",
    "\n",
    "# We need to know how many data points fall in each group.\n",
    "num_rec = np.sum(rec)\n",
    "# Vector of indices for later shuffling.\n",
    "idx_vec = np.arange(len(rec))\n",
    "\n",
    "# Note that this function shuffles in-place.\n",
    "from numpy.random import shuffle\n",
    "\n",
    "# Number of random shuffles.\n",
    "numshuf = 1000\n",
    "null_mus = []\n",
    "# Iterating over shuffles.\n",
    "for i in range(numshuf):\n",
    "    shuffle(idx_vec)\n",
    "    # Creating null groups of 'rec' and 'not-rec'.\n",
    "    null_rec = idx_vec[:num_rec]\n",
    "    null_nrec = idx_vec[num_rec:]\n",
    "    diff_mu_null = np.mean(avg_pow[null_rec, :], 0)-np.mean(avg_pow[null_nrec, :], 0)  #recomputing the test statistic\n",
    "    null_mus.append(diff_mu_null)   \n",
    "null_mus = np.array(null_mus)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Count')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQAAAADQCAYAAADoO7/9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAPIElEQVR4nO3df5BdZX3H8feH0CAQKKFJFAPThTEVKbVAg4PSUlSgDDAGrVIcrEnBUqbFqlRKKG2RsUxjazvUtkJTBKOmUExhjICS8ENpnRCzCZAQAgYhwEqaDWrFH51g8Ns/zrPNZdkf594959y7+3xeMzv33OeeH9+9yf3uc859zvNVRGBmedqr2wGYWfc4AZhlzAnALGNOAGYZcwIwy5gTgFnG9u52ABMxa9as6Ovr63YYZj1h/fr1z0fE7Ha2mdQJoK+vj/7+/m6HYdYTJD3d7jY+BTDLmBOAWcacAMwy5gRgljEnALOMTepvAax+fYvvGHedbUvObCASq4N7AGYZcwIwy5gTgFnGnADMMuYEYJYxJwCzjDkBmGXMCcAsY04AZhmrLQFIukHSoKRHWtr+VtJjkjZKuk3SQS2vXS7pCUmPS/qtuuIysz3q7AF8Fjh9WNtq4OiIeCPwLeByAElHAecCv5y2+bSkaTXGZmbUeC9ARNwvqW9Y26qWpw8A707LC4CbI2IX8JSkJ4A3AWvqis/KjfO3qa2b1wDOB76SlucCz7a8NpDazKxGXUkAkq4AdgPLh5pGWG3EooWSLpTUL6l/586ddYVoloXGE4CkhcBZwHmxpzLpAHBYy2qHAs+NtH1ELI2I+RExf/bstiZANbNhGk0Akk4HLgPeERE/aXlpJXCupH0kHQ7MA77ZZGxmOartIqCkm4CTgVmSBoArKa767wOslgTwQERcFBGbJd0CPEpxavBHEfFSXbFZtTxpyORV57cA7x2h+TNjrH81cHVd8ZjZK3kkoFnGnADMMuYEYJYxJwCzjDkBmGXMCcAsY04AZhlzAjDLmBOAWcacAMwy5gRgljEnALOMuTz4FOXpvqwM9wDMMuYEYJaxpusCHCxptaSt6XFmy2uuC2DWsKbrAiwG7omIecA96bnrAph1SW0JICLuB743rHkBsCwtLwPObmm/OSJ2RcRTwFBdADOrUdPXAF4dEdsB0uOc1F66LoCnBTerTq9cBCxdF8DTgptVp+kEsEPSIQDpcTC1l64LYGbVaToBrAQWpuWFwJda2l0XwKxhTdcFWALcIukC4BngPQCuC2DWHU3XBQB4+yjruy6AWcN65SKgmXWBE4BZxpwAzDLmBGCWMScAs4w5AZhlzAnALGNOAGYZcwIwy5gTgFnGnADMMuZpwSchT/ltVXEPwCxjTgBmGevKKYCkjwAfoJj2axPwe8B+wL8DfcA24JyI+H434rPqlTlt2bbkzAYisVaN9wAkzQX+GJgfEUcD0yimBB9xynAzq0+3TgH2BvaVtDfFX/7nGH3KcDOrSeMJICK+A3ySYkqw7cAPImIVo08ZbmY1KZUAJJ1Ypq3kvmZS/LU/HHgtsL+k97WxvesCmFWkbA/gH0u2lXEK8FRE7IyInwK3Am9h9CnDX8Z1AcyqM+a3AJLeTPHhnC3pkpaXDqS4eNeJZ4ATJO0H/C/FJKH9wI8ppgpfwsunDDezmoz3NeB0YEZa74CW9heAd3dywIhYK2kFsIFiCvAHgaXpOK+YMtzM6jNmAoiIrwNfl/TZiHi6qoNGxJUUdQJa7WKUKcPNrB5lBwLtI2kpxSCd/98mIt5WR1Bm1oyyCeCLwHXA9YAr9phNEWUTwO6IuLbWSMyscWW/BvyypD+UdIikg4d+ao3MzGpXtgcwVNH30pa2AI6oNhwza1KpBBARh9cdiJk1r1QCkPT+kdoj4nPVhmNmTSp7CnB8y/KrKL6v3wA4AZhNYmVPAT7Y+lzSzwOfryUiy5YnDWlep7cD/wSYV2UgZta8stcAvkxx1R+Km4DeANxSV1Bm1oyy1wA+2bK8G3g6IgZqiMfMGlTqFCDdFPQYxR2BM4EX6wzKzJpRdkagc4BvUtyiew6wVlJHtwObWe8oewpwBXB8RAwCSJoN3A2sqCuwnLnyjzWl7LcAew19+JPvtrHtK0g6SNIKSY9J2iLpzen+gtWStqbHmZ3u38zKKfsh/qqkuyQtkrQIuAO4cwLH/QfgqxFxJPCrwBZcF8CscePNCfg6ium6L5X0LuDXAQFrgOWdHFDSgcBJwCKAiHgReFHSAuDktNoy4GvAZZ0cw8zKGa8HcA3wQ4CIuDUiLomIj1D89b+mw2MeAewEbpT0oKTrJe2P6wKYNW68BNAXERuHN0ZEP8X0YJ3YGzgOuDYijqWYDbh0d991AcyqM14CeNUYr+3b4TEHgIGIWJuer6BICK4LYNaw8RLAOkm/P7wxTd29vpMDRsR/A89Ken1qejvwKLCSPROPuC6AWQPGGwfwYeA2Seex5wM/n6JewDsncNwPAsslTQeepCgPvheuC2DWqPHqAuwA3iLprcDRqfmOiLh3IgeNiIcoEslwrgtg1qCy8wHcB9xXcyxm1rDGy4ObWe9wAjDLmBOAWcbK3g1oFfGdftZLnABsUhkvgXrS0Pb4FMAsY04AZhlzAjDLmBOAWcacAMwy5gRgljEnALOMOQGYZcwJwCxjXUsAkqalSUFvT89dF8CsYd3sAXyIoh7AENcFMGtYVxKApEOBM4HrW5oXUNQDID2e3XRcZrnpVg/gGuBPgZ+1tJWqC+Bpwc2q03gCkHQWMBgRnc4q7GnBzSrSjduBTwTeIekMiroDB0r6AqkuQERsH6sugJlVp/EeQERcHhGHRkQfcC5wb0S8D9cFMGtcL40DWAKcKmkrcGp6bmY16uqMQBHxNYoqwETEd3FdALNG9VIPwMwa5gRgljEnALOMOQGYZcwJwCxjTgBmGXMCMMuYE4BZxpwAzDLmBGCWMScAs4y5OnCFXPrbJhsnAJtSyiRhlxDfw6cAZhlzAjDLWDfmBDxM0n2StkjaLOlDqd11Acwa1o1rALuBP4mIDZIOANZLWg0soqgLsETSYoq6AJd1IT6b4nydYI9uzAm4PSI2pOUfUhQHmYvrApg1rqvXACT1AccCa3FdALPGdbM24AzgP4APR8QLZbdzXQCz6nSrNNjPUXz4l0fEral5R6oHgOsCmDWjG98CCPgMsCUi/r7lJdcFMGtYtyoD/S6wSdJDqe3PKOoA3CLpAuAZ4D1diM0sK40ngIj4L0CjvOy6AGYN8r0AZiPIZayAhwKbZcw9gJJ8q69NRe4BmGXMCcAsY04AZhlzAjDLmBOAWcacAMwy5gRgljGPA8Df8Vu+3AMwy5h7AGYdmgr3C0z5BODuvdnoeu4UQNLpkh6X9ESaHdjMatJTPQBJ04B/Bk4FBoB1klZGxKPdjcysM1X1QOs6lei1HsCbgCci4smIeBG4mWK6cDOrQa8lgLnAsy3PB1KbmdWgp04BGHmqsHjZCtKFwIXp6Y8kPV5TLLOA52vadzt6IY5eiAF6I46uxKBPlIrjF9vdb68lgAHgsJbnhwLPta4QEUuBpXUHIqk/IubXfZzJEEcvxNArcfRCDFXG0WunAOuAeZIOlzQdOJdiunAzq0FP9QAiYreki4G7gGnADRGxucthmU1ZPZUAACLiTuDObsdBA6cZJfVCHL0QA/RGHL0QA1QUhyJi/LXMbErqtWsAZtagrBOApIMlrZa0NT3OHGW9GyQNSnqkk+0rimHEIdKSPibpO5IeSj9ntHn8MYdeq/Cp9PpGSceV3bahGLZJ2pR+9/5OYygZx5GS1kjaJemj7WzbUAztvxcRke0P8DfA4rS8GPjEKOudBBwHPNLJ9hONgeKC6LeBI4DpwMPAUem1jwEf7fD3H3W/LeucAXyFYozGCcDastvWHUN6bRswq4L/C2XimAMcD1zd+p43/F6MGEOn70XWPQCKYcbL0vIy4OyRVoqI+4Hvdbp9BTHUNUS6zH4XAJ+LwgPAQal8e1UxTSSGKo0bR0QMRsQ64Kcd/A51x9CR3BPAqyNiO0B6nNPw9mX3Md4Q6YtT1/iGNk9Dygy9Hm2dqoZtTyQGKEaKrpK0Po0S7dREfp8m34uxtP1e9NzXgFWTdDfwmhFeumISxTDWEOlrgY+n5x8H/g44v4L9jrdOmW3rjgHgxIh4TtIcYLWkx1KPrY446ti2yv20/V5M+QQQEaeM9pqkHZIOiYjtqUs52ObuS21fQQyjDpGOiB0t+/pX4PY24h936PUY60wvsW3dMRARQ4+Dkm6j6EZ3kgDKxFHHtpXtp5P3IvdTgJXAwrS8EPhSw9uX3ceoQ6SHnQu/E3hkhO1HU2bo9Urg/elK/AnAD9KpSlXDtjuOQdL+kg4AkLQ/cBrt/f7txlHHtpXsp+P3YqJXTyfzD/ALwD3A1vR4cGp/LXBny3o3AdspLrwMABeMtX1NMZwBfIviKvEVLe2fBzYBG9N/lkPaPP4r9gtcBFyUlkUxScu303HmjxdTB+9BRzFQXC1/OP1snkgMJeN4Tfr3fwH4n7R8YMPvxYgxdPpeeCSgWcZyPwUwy5oTgFnGnADMMuYEYJYxJwCzjDkBTHGSXkp3h22W9LCkSyTtlV6bL+lTaXkfSXendX9H0m+kbR6StG93fwuri78GnOIk/SgiZqTlOcC/Ad+IiCuHrXcCxZ2Iv5meX0dx192NJY8jiv9PP6v0F7BauQeQkYgYpJhS/eI0qu5kSbenxPAF4Jj0F/8PgHOAv5S0HEDSpZLWpZuOrkptfZK2SPo0sAE4TNJp6X71DZK+KGko+WyTdFVq3yTpyNQ+Q9KNqW2jpN9O7SPux6rlBJCZiHiS4t99TkvbIPAB4D8j4piI+BeKUYWXRsR5kk4D5lGMLT8G+DVJJ6XNX09xq+6xwI+BPwdOiYjjgH7gkpbDP5/arwWGJrP4C4qhvb8SEW8E7pU0a5z9WEWm/M1ANqKR7joby2np58H0fAZFQngGeDqKe/ShmKzjKOAbxRkB04E1Lfu5NT2uB96Vlk+hGPMOQER8X9JZ4+zHKuIEkBlJRwAvUdx1+IaymwF/nXoGrfvqo/ir37re6oh47yj72ZUeX2LP/z0x8u2/Y+3HKuJTgIxImg1cB/xTtHf19y7g/Jbz+bnpusFwDwAnSnpdWm8/Sb80zr5XARe3xDizw/1YB5wApr59h74GBO6m+MBd1c4OImIVxbcHayRtAlYAB4yw3k5gEXCTpI0UH+Qjx9n9XwEzJT0i6WHgrR3uxzrgrwHNMuYegFnGnADMMuYEYJYxJwCzjDkBmGXMCcAsY04AZhlzAjDL2P8BhsBkieAbSYQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 252x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Inspect the null distribution to see that everything worked as intended.\n",
    "plt.figure(figsize=(3.5, 3))\n",
    "ax = plt.subplot(111)\n",
    "plt.hist(null_mus[:, 11], 20)\n",
    "plt.xlabel('Difference')\n",
    "plt.ylabel('Count')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks pretty good! Null distributions don't *have* to look any particular way, but if you see very strange-looking distributions, or an unreasonable range of values, you may want to go back and check your code. (Or at least understand why the distribution should look that way.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "One-sided: 0.19999999999999996\n",
      "Two-sided: 0.3999999999999999\n",
      "-0.8397919451081097\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQAAAADQCAYAAADoO7/9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAQJElEQVR4nO3df5QdZX3H8fcHaBAMCDQJIJAuVCqmVoEuHpCWWvlRChTUKsLRigWbcloUS0GCtIUc5TTUH7W2FZoKiIpQoHCMFCUBEVoPIJtAEn4alACBlASRUn+cWPDbP+ZZclnv7p17d2bu3X0+r3P23Jlnfn3vTfa7z8x95juKCMwsT1v0OwAz6x8nALOMOQGYZcwJwCxjTgBmGXMCMMvYVv0OYDJmzZoVQ0ND/Q7DbCAsX778mYiY3c02tSUASZcCxwAbIuL1Y5adCXwCmB0Rz6S2c4BTgBeBD0XETZ2OMTQ0xMjISOWxm01Fkh7rdps6TwG+ABw5tlHSHsDhwOMtbfOAE4BfT9t8TtKWNcZmZtSYACLiduDZNov+HvgI0DoE8TjgqojYFBGPAo8Ab6orNjMrNHoRUNKxwJMRsXLMot2AJ1rm16U2M6tRYxcBJW0LnAsc0W5xm7a2NylImg/MB5g7d25l8ZnlqMkewK8CewIrJa0FdgdWSNqF4i/+Hi3r7g481W4nEbE4IoYjYnj27K4ueJrZGI31ACJiNTBndD4lgeGIeEbSEuArkj4NvBrYG/hOU7FNeeef337arIM6vwa8EngLMEvSOuC8iLik3boRcb+kq4EHgBeAP4+IF+uKbdpZuHDztBOAdaG2BBARJ3ZYPjRm/gLggrriMbNf5KHAZhlzAjDLmBOAWcacAMwy5gRglrEpfTuw1W9owX90XGftoqMbiMTq4B6AWcacAMwy5gRgljEnALOMOQGYZcwJwCxjTgBmGXMCMMtYbQlA0qWSNki6r6XtE5IekrRK0vWSdmhZdo6kRyQ9LOn36orLzDZruiz4MuD1EfEG4LvAOeCy4Gb9UmdBkNslDY1pW9oyeyfwzjT9Ullw4FFJo2XB76grPis3zNemt35eAzgZ+HqaLl0WXNJ8SSOSRjZu3FhziGbTW18SgKRzKWr/XTHa1Ga1tmXBXRXYrDqN3w0o6SSKZwYeGhGjv+Sly4KbWXWafjLQkcDZwLER8ZOWRUuAEyRtLWlPXBbcrBGNlgWnuOq/NbBMEsCdEXGqy4JPba4ZMHU1XRa87XMB0vouC27WMI8ENMuYE4BZxpwAzDLmBGCWMScAs4w5AZhlzAnALGNOAGYZcwIwy5gTgFnGnADMMuaHg05TrvZjZbgHYJaxpqsC7yRpmaQ16XXHlmWuCmzWsKarAi8AbomIvYFb0ryrApv1SW0JICJuB54d03wccHmavhx4W0v7VRGxKSIeBUarAptZjZq+BrBzRKwHSK9zUnvpqsBmVp1BuQhYuiqwy4KbVafpBPC0pF0B0uuG1F66KrDLgptVp+kEsAQ4KU2fBHy1pd1Vgc0a1nRV4EXA1ZJOAR4H3gXgqsBm/dF0VWCAQ8dZ31WBzRo2KBcBzawPnADMMuYEYJYxJwCzjDkBmGXMCcAsY04AZhlzAjDLmBOAWcacAMwy5gRgljFXBZ6Cxlb8XTvBMrOJuAdglrFSCUDSwWXazGxqKXsK8I/A/iXaSpH0F8AHKMp+rQb+GNgW+DdgiKJXe3xE/LCX/dvgKXNqsnbR0Q1EYq0mTACSDgLeDMyWdEbLou2Bnsp2S9oN+BAwLyJ+mgqBnADMoygZvkjSAoqS4Wf3cgwzK6fTKcAMYCZFotiu5ed54J2TOO5WwDaStqL4y/8U45cMN7OaTNgDiIjbgNskfSEiHqvigBHxpKRPUpQE+ymwNCKWSnpZyXBJc9ptL2k+MB9g7ty5VYRklq2y1wC2lrSY4vz8pW0i4q3dHjA9Duw4YE/gOeAaSe8tu31ELAYWAwwPD7ctHW5m5ZRNANcAFwOfByZbrPMw4NGI2Agg6TqK6wxPS9o1/fVvLRluZjUpmwBeiIiLKjrm48CBkralOAU4FBgBfkxRKnwRLy8ZbmY1KZsAvibpz4DrgU2jjREx9tl/HUXEXZKuBVZQlAC/h6JLP5M2JcPNrD5lE8DowzzOamkLYK9eDhoR51E8J6DVJsYpGW5m9SiVACJiz7oDMbPmlUoAkt7Xrj0ivlhtOGbWpLKnAAe0TL+Coqu+AnACMJvCyp4CfLB1XtKrgC/VEpGZNabX24F/QvEEXzObwspeA/gaxVV/KG4Ceh1wdV1BmVkzyl4D+GTL9AvAYxGxroZ4zKxBZa8B3CZpZzZfDFxTX0iWK9cMaF7ZikDHA9+hGJ13PHCXpMncDmxmA6DsKcC5wAERsQFA0mzgZuDaugIzs/qV/RZgi9Ff/uQHXWxrZgOqbA/gG5JuAq5M8+8GbqwnJDNrSqeagK8Bdo6IsyS9A/gtQMAdwBUNxGdmNerUA/gM8FGAiLgOuA5A0nBa9ge1RpcpP9zDmtLpPH4oIlaNbYyIEYryYD2RtIOkayU9JOlBSQdJ2knSMklr0uuOve7fzMrplABeMcGybSZx3H8AvhER+wBvBB6kKAN+S0TsDdyS5s2sRp0SwN2S/mRsY6ras7yXA0raHjgEuAQgIn4WEc/hsuBmjet0DeDDwPWS3sPmX/hhiucFvL3HY+4FbAQuk/TGtN/TKS42uiy4WYMm7AFExNMR8WZgIcXjutYCCyPioIj47x6PuRXFI8Uuioj9KIqBlu7uR8TiiBiOiOHZs2f3GIKZQfl7AW4Fbq3omOuAdRFxV5q/liIBuCy4WcMaH82Xeg5PSHptajoUeABYwubioy4LbtaAsiMBq/ZB4ApJM4DvUzwdeAtcFtysUX1JABFxL8XFxLFcFtysQb6hxyxjTgBmGXMCMMtYvy4CZss3+tggcQKwKaVTAnXNwO74FMAsY04AZhlzAjDLmBOAWcacAMwy5gRgljEnALOMOQGYZaxvCUDSlpLukXRDmndVYLOG9bMHcDpFNeBRrgps1rC+JABJuwNHA59vaXZVYLOG9asH8BngI8DPW9peVhUYaFsV2Myq03gCkHQMsCEien2uwHxJI5JGNm7cWHF0ZnnpRw/gYOBYSWuBq4C3SvoyqSowwERVgV0W3Kw6/agKfE5E7B4RQ8AJwDcj4r24KrBZ4wZpHMAi4HBJa4DD07yZ1aivBUEi4lvAt9L0D3BVYLNGDVIPwMwa5gRgljEnALOMOQGYZcwJwCxjTgBmGXMCMMuYE4BZxpwAzDLmBGCWMT8bsEJ+8KdNNe4BmGXMPQCbVsr0wvwE4c3cAzDLWD9Kgu0h6VZJD0q6X9Lpqd1lwc0a1o9TgBeAv4yIFZK2A5ZLWga8n6Is+CJJCyjKgp/dh/hsmvNpwmb9KAm2PiJWpOn/pXg2wG64LLhZ4/p6DUDSELAfcBcly4K7KrBZdfr5aLCZwL8DH46I58tu56rAZtXp15OBfonil/+KiLguNZcqC25m1enHtwACLgEejIhPtyxyWXCzhvXjW4CDgT8CVku6N7V9lKIM+NWSTgEeB97Vh9jMstJ4AoiI/wI0zmKXBTdrkIcCm7WRy1gBDwU2y5h7ACX5Vl+bjtwDMMuYE4BZxpwAzDLmBGCWMScAs4w5AZhlzF8D4q/4LF/uAZhlzD0Asx5Nh+HC0z4BuHtvNr6BOwWQdKSkhyU9koqDmllNBioBSNoS+Gfg94F5wImS5vU3KrPpa9BOAd4EPBIR3weQdBVFteAH+hqVWY+qOgWt61rCQPUAKMqDP9Eyvy61mVkNBq0H0K5SULxsBWk+MD/N/kjSwzXFMgt4pqZ9d6NjHC/70C48pi8xNGQQ4uhLDLqwVBy/0u1+By0BrAP2aJnfHXiqdYWIWAwsrjsQSSMRMVz3caZCHIMQw6DEMQgxVBnHoJ0C3A3sLWlPSTOAEyiqBZtZDQaqBxARL0g6DbgJ2BK4NCLu73NYZtPWQCUAgIi4Ebix33HQwGlGSYMQxyDEAIMRxyDEABXFoYjovJaZTUuDdg3AzBqUdQKQtJOkZZLWpNcdx1nvUkkbJN3Xy/YVxdB2iLSk8yU9Kene9HNUl8efcOi1Cp9Ny1dJ2r/stg3FsFbS6vTeR3qNoWQc+0i6Q9ImSWd2s21DMXT/WUREtj/A3wEL0vQC4MJx1jsE2B+4r5ftJxsDxQXR7wF7ATOAlcC8tOx84Mwe3/+4+21Z5yjg6xTDDQ4E7iq7bd0xpGVrgVkV/F8oE8cc4ADggtbPvOHPom0MvX4WWfcAKIYZX56mLwfe1m6liLgdeLbX7SuI4aUh0hHxM2B0iPRkldnvccAXo3AnsEN6enNVMU0mhip1jCMiNkTE3cD/9fAe6o6hJ7kngJ0jYj1Aep3T8PZl99FpiPRpqWt8aZenIWWGXo+3TlXDticTAxQjRZdKWp5GifZqMu+nyc9iIl1/FgP3NWDVJN0M7NJm0blTKIaJhkhfBHwszX8M+BRwcgX77bROmW3rjgHg4Ih4StIcYJmkh1KPrY446ti2yv10/VlM+wQQEYeNt0zS05J2jYj1qUu5ocvdl9q+ghjGHSIdEU+37OtfgRu6iL/j0OsJ1plRYtu6YyAiRl83SLqeohvdSwIoE0cd21a2n14+i9xPAZYAJ6Xpk4CvNrx92X2MO0R6zLnw24H72mw/njJDr5cA70tX4g8E/iedqlQ1bLvnGCS9UtJ2AJJeCRxBd++/2zjq2LaS/fT8WUz26ulU/gF+GbgFWJNed0rtrwZubFnvSmA9xYWXdcApE21fUwxHAd+luEp8bkv7l4DVwKr0n2XXLo//C/sFTgVOTdOiKNLyvXSc4U4x9fAZ9BQDxdXylenn/snEUDKOXdK///PAc2l6+4Y/i7Yx9PpZeCSgWcZyPwUwy5oTgFnGnADMMuYEYJYxJwCzjDkBTHOSXkx3h90vaaWkMyRtkZYNS/psmt5a0s1p3XdL+u20zb2Stunvu7C6+GvAaU7SjyJiZpqeA3wF+HZEnDdmvQMp7kT8nTR/McVdd5eVPI4o/j/9vNI3YLVyDyAjEbGBoqT6aWlU3Vsk3ZASw5eBfdNf/D8Fjgf+RtIVAJLOknR3uuloYWobkvSgpM8BK4A9JB2R7ldfIekaSaPJZ62khal9taR9UvtMSZeltlWS/jC1t92PVcsJIDNRPHVpC1ruOkyJ4QPAf0bEvhHxLxSjCs+KiPdIOgLYm2Js+b7Ab0o6JG3+WopbdfcDfgz8FXBYROwPjABntBz+mdR+ETBazOKvKYb2/kZEvAH4pqRZHfZjFZn2NwNZW+3uOpvIEennnjQ/kyIhPA48FsU9+lAU65gHfLs4I2AGcEfLfq5Lr8uBd6TpwyjGvAMQET+UdEyH/VhFnAAyI2kv4EWKuw5fV3Yz4G9Tz6B1X0MUf/Vb11sWESeOs59N6fVFNv/fE+1v/51oP1YRnwJkRNJs4GLgn6K7q783ASe3nM/vlq4bjHUncLCk16T1tpX0ax32vRQ4rSXGHXvcj/XACWD622b0a0DgZopfuIXd7CAillJ8e3CHpNXAtcB2bdbbCLwfuFLSKopf5H067P7jwI6S7pO0EvjdHvdjPfDXgGYZcw/ALGNOAGYZcwIwy5gTgFnGnADMMuYEYJYxJwCzjDkBmGXs/wFdbMibR1j5cQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 252x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Step 6: Get a p-vaue or z-score.\n",
    "\n",
    "# Plot where the true value falls relative to the null distribution.\n",
    "plt.figure(figsize=(3.5, 3))\n",
    "ax = plt.subplot(111)\n",
    "plt.hist(null_mus[:, 11], 20)\n",
    "plt.xlabel('Difference')\n",
    "plt.ylabel('Count')\n",
    "plt.vlines([diff_mu[11]], ymin=ax.get_ylim()[0], ymax=ax.get_ylim()[1], color='r', linewidths=3)\n",
    "\n",
    "# Quantify with a p-value.\n",
    "# If we're doing a 2-sided test, multiply the resulting p-value by 2. \n",
    "pval = np.sum(null_mus[:, 11]>diff_mu[11])/float(null_mus.shape[0])\n",
    "print('One-sided: '+str(1.-pval))\n",
    "print('Two-sided: '+str((1.-pval)*2))\n",
    "\n",
    "# Quantify with a z-score.\n",
    "Z = (diff_mu[11]-np.mean(null_mus[:, 11]))/np.std(null_mus[:, 11])\n",
    "print(Z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The true value falls pretty close to the middle, which is why it has a nonsignificant p-value and a z-score of -0.744. Nothing significant here!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'T-stat')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAADTCAYAAACsut7+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2deXxU5fX/34cQICwahLCDqCi0qICyKbhgKSiaCCjWDWu10qq1tWoqauuC+pUKtfbl+vNX5SuFqlVoFKE/qgJuQDAYEVlCLbIlgGxhM4Qs5/fHzOAkmcksuTN3ZnLer9e8nLlz57nnRj73Oc95znMeUVUMw0hNmrhtgGEYscMEbhgpjAncMFIYE7hhpDAmcMNIYUzghpHCNHXbgGhp37699uzZ020zDMN1Vq5cuVtVswJ9l7QC79mzJwUFBW6bYRiuIyKbg31nLrphpDAmcMNIYUzghpHCmMANI4VJ2iCbYbhNXmEx0xYWUVJaRpfMDHJH92bsgK6Otb9u3Tp69OhBq1atom7DenDDiIK8wmLum7ua4tIyFCguLeO+uavJKyx27BpFRUXk5uY2qA0TuGFEwbSFRZRVVNU4VlZRxbSFRQ1qV1VZtWoVAGPHjmX69OkNas8EbhhRUFJaFtHxcDh69Cg333wzZ599Nl988QUALVu2jLo9MIEbRlR0ycwIeFyBYVMXReyq7969m5EjRzJjxgzuv/9+zjzzTAesNIEbRlTkju5NRnpawO8iHY+vWbOGwYMHs2LFCv7+978zZcoUmjRxRpoWRTeMENQXLZ+2sIjiAG65bzweTlR93rx5lJWV8eGHHzJkyBBHbZdkrck2cOBAtVx0wwnqE7AvWu4fUMtIT+OJ8WccO+ekyfMJpCIBvpl6acBrqirFxcV069YNVWX37t1kZQVcLxISEVmpqgMDfZcQLrqItBCRFSKySkTWiMgjbttkNA5CTXeFEy0PNh4Pdvzo0aNMmjSJfv36UVxcjIhELe5QJITAgXLgIlXtB/QHLhaRoS7bZDQCQgk4nGh5oPF4RnoauaN71/ndnj17GDVqFH/961+59dZb6dy5c0NvoV4SYgyunnHCIe/HdO8rOccORlIRSsBdMjMCjrH9e2f/8Xh9WW3r1q0jOzubbdu2MXv2bK699lqnbiMoCSFwABFJA1YCvYDnVDXfZZOMRkAoAeeO7h1wDF67dx47oGvIgNr06dM5dOgQS5YsYejQ+DioieKio6pVqtof6AYMFpHTa58jIpNEpEBECnbt2hV/I42UI5R7PXZAV54YfwZdMzMQoGtmRo0AWyhUlQMHDgDwzDPP8Nlnn8VN3JCgUXQReQg4rKpB8/Qsim44RawWjVRUVHDHHXfw6aefsmzZMlq3bu2AtXWpL4qeEC66iGQBFapaKiIZwEjgjy6bZTQSwnGvI2Xv3r1MmDCBRYsWMXny5AannEZLQggc6Ay86h2HNwH+oarvumyTYURFUVERl112GVu2bGHmzJlMnDjRNVsSQuCq+iUwwG07DMMJbr31Vvbv38+iRYsYNmyYq7YkhMANIxWoqKggPT2dV199laqqKhKhrLcJ3DAaSEVFBXfeeSfFxcXMnTuX7t27u23SMRJmmswwkpF9+/ZxySWX8Pzzz9O7d28SbVbKenDDiJINGzaQnZ3NN998w4wZM7jxxhvdNqkOJnDDCECoufGqqipycnLYu3cvixYtYvjw4S5aGxwTuGHUovYSUd8KM/DMmasqaWlpzJw5k6ysLE466SQ3za2XhMxkCwfLZDNixbCpiwLnp7dpxtk73iEzM5NHH33UBcsCk/CZbIYRDrGuQ+4j0Aqz6iOHKHzjjyzbVMg999yDqiIijl/baUzgRlIQym12ktorzCr2lfDtW1Oo2r+Dl19+mZtuusnR68USmyYzkoJY1SEPhP8Ks+qjR9g5+1607ACPvvh6UokbrAc3koRY1CEPRs0CDtBr7K+58yej+GX2uY5fK9aYwI2kIJzKKv40ZLxeWVnJhzOn85tzzuGqq64CAhdOTAYSwkUXke4islhE1nmLLv7GbZuMxCKSumcN2Tds//79ZGdn8/TTT1NYWOiU+a6REAIHKoG7VfUHwFDgdhH5ocs2GQlEJJVVoh2v//e//+Wcc87h/fff56WXXuKJJ55w8hZcISFcdFXdDmz3vj8oIuuArsBaVw0zHKWh01zhFmaIZry+Y8cOhgwZgqry3nvvceGFF4ZtVyKTKD34MUSkJ5614VZ0MYWIx3a7PiKtUw7QqVMn7r33XvLz81NG3JBgAheR1sAc4E5VPRDgeyu6mKS4Nc3lI9B4vaqqismTJx8ba+fm5tKrVy/H7XGThHDRAUQkHY+4Z6vq3EDnqOpLwEvgSVWNo3lGA4nHNJdvCFBcWoZ/klnbluk8lN23hnt/4MABrrnmGhYsWECrVq0YMCA1CwolhMDFk/P3MrBOVZ9y2x7DeSKd5oqU2plu/kssjlRU1zh348aN5OTkUFRUxIsvvsgvfvELR2xIRBLFRR8GTAQuEpEvvK8xbhtlOEck01zREGgI4MN/KFBUVMSQIUMoKSlh4cKFKS1uSJAeXFU/wbMZo5GihLu9T7SEcvV935988slcccUV3H333Zx66qmOXDuRSQiBG42DWNQf9xFsCACg1VUcWP4WPe7YR/euXcj9xR849dTY2JFoJIqLbhgNItAQAKC6/Dt2zX2MfR//jUPrPorp9FwiYj24kRL4DwGKS8tIE6G8dAffzplCxe6tnPDjW2lzlien3Dcmb6g3Ea/16Q3BBG6kDP5DgMLCQkaPvpGqA9/R4aopZPTsX+Pchk7PxXN9ekMwF91IePIKixk2dREnTZ7PsKmLwnKve/TowcCBA+l3+3N1xA2ezefDbSsQ8UzcaQgmcCOhiSTFtaqqihdeeIHy8nLatWvHggULePD6kQHH5oRoKxTxXJ/eEEzgRkITbk956NAhxo8fz2233cabb7557Lj/KrRARNvrRpPv7gYmcCOhCaen3Lx5M8OGDePdd9/lmWee4brrrqtx7tgBXfl08kVBEy2i6XVjnbjjFBZkMxKaYPPbx2ekA5Cfn09OTg7l5eX861//YtSoURG3FU2vG+vEHacwgRsJTe7o3uS+uYqK6ppriw4frSSvsJhTW7emW7duzJ49mz59+hz7PtAUVu7o3jUi39CwXjeWiTtOYS66kdCMHdCV1i1q9kOq1exf9ylP/r/19O3bl4KCgjrizn1rVY3AXO5bqwDCrgqTKlgPbiQ8pd9VHHtffbSM3e/+ibL/LKdJs5bAj+psQPDIvDVUVNXs8SuqlEfmraHwwVEpLejamMCNiIl3Bpdv7Fx54Fu+nfMoFbs20/ZHt3Byv6EBz9/n90AI53gqkzAuuoi8IiLfishXbttiBCeepZd85I7uDTs3sH3mXVSW7qTDlQ/R8Zzx/O7iPqF/3MhJGIED/wtc7LYRRv24kcE1dkBXru2XSbMWreg8cTq9zhpe79g50xthD/d4KpMwLrqqfuQtuGgkMLHI4Arm8ldXV7NixQqGDh3K/9z5Mx6+7TqaNWsWsr2Hc/rWibynNxEezukbtY3JSiL14CGxoovu43QGVzCX//Wl/+Gqq65i+PDhfPWVZ9QWjrjB0+NPm9CvRrR82oR+jSq45iNhevBwsKKL7uP0XHIgl//gnp38fMIdfLf9v0yfPp2+fev2vKECfckwRx0Pkkrghvs4ncFV27UvLyli19zHqK44wrx587j00rr7giXLUs1EwARuRIyTvWPt9NEjW1YjTZvR75bpAcUN9Qf6TOA1SZgxuIi8BiwDeovINhG52W2bjMiIZt127ujetGgqVOzbDsBxQ67gpFue46EbRgf9TbIs1UwEwhK4iHwW5PhypwxR1WtUtbOqpqtqN1V92am2jdgT7fz4qN5tabfiBb79211UHdpLt7YtefLaofX2xMmyVDMRCLcHDxZBOc0pQ4zkJpr58eLiYs4//3yWvT+fJ6Y8yJZnrufTyReFdLOTZalmIlDvGFxEXvG+be733kdPYF0sjDKSj0jd5oKCAnJycjh48CBvv/022dnZYV8rWZZqJgKhgmzFQd4rsBJ4w3GLjKQk0rXWzz33HM2aNWPp0qWcccYZEV/PpsHCo16Bq+ofwDPWVtX58THJSEbCmR9XVfbu3Uu7du14/vnnOXjwIB06dHDD3EZDWNNkqjrfu/tnL6A9ftsMqepHMbLNcJlIVo2FcpvLysq46aab+PLLL1mxYgWtWrUiI8OCYrEmLIGLyDnAm8BxQEvgsPe/24EeMbPOcI1okkmCuc0lJSWMHTuWgoICpk6dSsuWLWNnuFGDcBNdngb+AkwH9qpqWxF5BCiNmWWNjETbJcOpZJLPP/+cnJwcSktL+ec//8nll18e8Lxw7z/R/k6JTrgC7w38SVXVr3rG48BG4M+xMKwxkYipl04kk6gqE2+5jV2HK2h/1RM8ua4N2qO4zj2Fe/+J+HdKdMKdBz8ItPG+3yEifYBMv2NGA0jEXTIakkyiqhw5coS3vyih/Lw76HD9n0jvcHLQ5Jdw7z8R/06JTrgCzwMu876fASwGCoC5sTCqseF06mU0KaO1iTaZpKysjOuuu44JEybw5L/WUtkik7RWbb//PoAgw71/S1GNnLAErqp3qOps7/sngWuA2wHLF3cAJ1MvnSqp5L8jSLgVSHfs2MGIESN47bXXOPfccynZXx7wvNqCDPf+LUU1csLNRX/K/7OqLlHVecCTMbGqkeFk6mUwN/bON77glPsW0DOCXt23I8g3Uy8NmUJaWFjIoEGDWL16NXPmzOG+++6ja9vA0fLaggz3/i1FNXLCddF/HuT4TU4ZIiIXi0iRiHwtIpOdajcZiKa3DEZ97mqVempkOFUo0TcU6Pm7dzh35GWUVVTxySefMH78eCC0IH2//+0bX9AivQmZGen13r+Tf6fGQqhc9Bt854nIRKixvdPJwG4njBCRNOA54MfANuAzEXlHVdc60X4y4FTqZbCU0do0dP10XmExk+d8SVlFJdIkjbY599Lq+PZspgMDvOfUl/xSOyK+77sKMtLT+PNP+tdrk6WoRkaoabJbvP9tBkzyO67ATuBnDtkxGPhaVTcCiMjrwOVAoxF4tNSeFx7RJ4s5K4vruOmBaEhw6o/vrmZb3jTS2rSn7QU/pVnHU6iAOg+NYIK0og3xIVQu+nkAIjJVVWPpNncFtvp93gYMqX2SiEzC+6Dp0cMS6ALNC89ZWcwVZ3dl8fpdIXvyUMGpYEklO3fupPDFOykvKeL4866v8ZtwHxoWEY8P4Y7BnxSRlgAi0kREJorINVJ7z5joCdROnaKKqvqSqg5U1YFZWVkOXTp5CdYLLl6/i08nX8TTP+lfZwzsI1RwKlg0/s9vvMegQYOo2LWJ9mPvI/Pcq2v8LtyItkXE40O4Al8A+LaReBx4ALgfT+qqE2wDuvt97gaUONR2SpJXWBy0h/b1gmMHdOWKs7uSVus5HE5wKtDD4/Chg/zuZ1dQXV3Nk/+bR/vTz6/xfSQRbYuIx4dIUlULve+vB4YDh4DVwN0O2PEZcKqInIRn3fnVwLUOtJuS+HrXYPh6wbzCYuasLD4WPYfvRRRqnBvIVW7SvCUnjLmTz565lc6dO3NKA/LCAwXgRvTJYtrCIn77xhfH2qt9juWeR0a4Aq8C0kXkNOCgqm72uuetnTBCVStF5FfAQiANeEVV1zjRdjIQ6QKKQL2rD/9esCGBLF80Xisr2LPwOTJOGkCrH15Ar4EX0rlzZ6DhEW3/3weKJ+S+tQqUYzuUWO555ITroi8EXgee9/4X4Id4los6gqouUNXTVPUUVX3cqXYTnWgyz+oLnvm73sECVsWlZSHTWHNH9yb96EF2vv4Ah796n8r9O2PqQgd6GFVUaY3th8ByzyMl3B7853imxCrwbBII0AGYEgObGhWR9rJ5hcUIASKQeMbW/r+pb07c/2EC1OhJpy0sYtN/1rFr7qNUHtpHVs69nHrOqJi6x5FEzy3SHj7h5qKXqerzqvp/gXbeY4t9+elG9EQ6XTRtYVFAcQuEldpZG/8e0edNbN6yle2zcqmurKDjtVPpctZFcdkDPBbnNnai2fhgg+NWNGIinS4KJnyl7ri0dmpnMHxt+ryJpse1J/O8iXS64Smadz6N0rKKuOwBXvthlJ4mpDepablF2iMjYXY2aaxEOl0UTPhdMzMCLhP1XzDStmXg/bG7ZGZw9OhRVv/jT5Rv/w8Axw3MoWmb9sfOicce4LXzzKdd2a/OLqGWex4Z0exN5lRyi0HkNb6DVS8d0Ser3moneYXFHDpSWae99DRh0uD2jBw5koOFH9M0szPNO58a8NqxHvsGi8qboKMn1GKTa1T1tVqHj4+hPY2SSKabgj0QQgXrpi0sqhORBmiybxsP33w7O3bs4K4nnmPB4ZODTsHZ2Df5CNWD/x+ghsBVtTp25jRuwp0PD/RA+O0bXwRs09frBup9j+7axJZZuXQ84Xg+/PBDBg8ezHmFxTwybw37vquoca6NfZOTUAI3dzxONLSgYKidRQJ9n96uO50GjWHF60/TrVu3Y9fyufQNySCz6qeJQSiBp4nICOoRuqouctakxskj89Y0aPlkqJ1FRvTJYvbyLVRXVVD68SzaDLycNm2zeOKZv9CtW3heQrhY9dPEIZTAmwMvE1zgiqfwg9EA8gqL67jEPsINbIUqrjBnZTGV3+1nV94TlG/9ivTMztz4y0l1yhLXzg1fvH5XwF64vh7a1nonDqEEflhVTcAxpr7pp0gCW/UVV9i//Rt2zZlC5cE9tM++h1Y/vJDF63cdO+f3eauZvXzLsSSa4tIyZi3fcux7/14YqLeHtrXeiUM002SGw9T3D9+JwNbGNYXs/MeDSHozOl07leZdete4bl5hcQ1xB8N/Lry+HjrSnUaN2BEq0SXmQTYRmSAia0SkWkQGxvp6iUiwf/iZGekNdmnzCotp3r4bGSedRecbnjombv/rBkt/DURJaVnIHtrWeicO9QpcVeOxc8lXwHig0e5SGkwQD+f0jbrNiooKbvztH7j3HyuheWuyxk6m6XHfb9XrL7hIXOcumRkh02ut+mni4HqqqqquU9VGvf7PaUHs2bOH0aNH8+rTj7F3/fI636eJ1Gg/XNfZ91AIpxyyTZElBqIarnMWW0RkCXCPqhbUc45/0cWzN2/eHCfrkof169eTnZ3Nli1bOO7Hv6LV6RcFPE+gRtWU2lNsApx7ygls2lMWURS99hQZeMRvPXjsEJGVqhpweBsXgYvI+0CnAF89oKpve89ZQgiB+zNw4EAtKAjr1EbDkiVLGDt2LM2bNycvL497PiwLWVnVJz5wpjTSsKmLAl6za2YGn04O/LAxGkZ9Ao9LFF1VR8bjOolKvFzWTp060a9fP2bOnMmJJ55Ibou6vWltfNHvUFsThYtNkSUWro/BUx2nNgMMRkVFBbNmzUJV6dOnD0uWLOHEE08EIl8P7gRWDjmxcF3gIjJORLYB5wDzRWSh2zY5SSz3tN67dy+XXHIJEydO5OOPPwagdql6//XgXcMUX0O2H7YpssTCdYGr6j9VtZuqNlfVjqo62m2bnCRWLuuGDRsYOnQoH330ETNmzOD8888P+ZtwxNdQj8OmyBILy2SLMbHI6vrggw+48soradq0KYsWLWL48OFh/S6c4hJO5JHbBoGJgwk8xoRa5VWbcAJy5eXlnHjiieTl5dGzZ8+I7AklPguSpRauu+ipTiQua33ucWVlJYsXLwZgzJgxrFy5MmJxh4MFyVKLhEl0iZRUnAcPNofcsXklxy97jg8++IC1a9fSu3f4AatIp+gsUSX5cH0evDHRkDnvQOKu2FvM53OmUFm6k1PG3cW671oTrryjKbwQaRFII7ExgTtIQyqZBIpSl21exe68J0Ca0PHqx6jsfnpElVGiDZhZkCx1sDG4g0Q75x1st9CjO74mrdUJdLrhKVp0Pz3s9nxYwMywHtxBohWU/4NBq6uo2FtMs/Y9OG7weNoMuJQmzVpE1J4PK7xgWA/uINFGoH2CrS4/zLdvPcKOWblUHS5FROqIO5z2fFhWmWE9uAP4AmvFpWV1dv4MR1BdMjPY9M1Gds2ZQsW+Ek4YdRtprTIDnhuJQC1gZpjAG0jtwJrCMZF3DVNQl7Tfy0NT7gKFjj95lBY9zgSgbct0Lj2zc9DKpuFgAbPGjQm8gQSqZ+4Td7jrn//76Xy6dOpEhyseZG9aW+tpDcdwXeAiMg3IBo4C/wV+pqql7loVHg2pZ15ZWcnu3bvp1KkTL7zwAuXl5Rx/vG37ZjhLIgTZ3gNOV9Uz8ew9fp/L9oRNtPXM9+/fT05ODhdeeCFlZWW0aNHCxG3EBNd7cFX9t9/H5cCVbtkSKdHUM9+4cSPZ2dls2LCBZ599lowMm7IyYofrAq/FTcAbbhsRLsHmmYPVM//oo48YP3481dXVLFy4kIsushplRmyJi4suIu+LyFcBXpf7nfMAUAnMrqedSSJSICIFu3btCnZa3Iiknrmqcv/999O+fXvy8/NN3EZcSIiiiyLyU+Ay4Edaz/I2VX0JeAk8q8kcNTIKwplnrqqq4s38jTzz0Va29P8lXdq1Yc2hlpwaJxutRnnjxnUXXUQuBu4FLlDV79y2J1Lqm2c+cOAAF106nvXb99PuigdJa92WneXEbStd28bXSIQo+rNAG+A9EflCRF502yAn2LRpE8OGDWPl0sU0P2UwIt//qZ0quhiKWBZ8NJID13twVe3ltg1O88knnzBu3DgqKyvpOGEKLXr2r3NOPFZ02WoyIxF68JTi6NGjTJw4kbZt25Kfn88p/c8JeF48VnRZ+SXDBO4QVVVVVFVV0axZM+bNm0d+fj6nnXaaqyu6km01WUPqsRuBcd1FTwUOHjzIddddR58+fXjyySc5/fTTj33n5oquZFpNZgHB2GACbyCbN28mOzubtWvXMmrUqIDnuLmiK1lWkzlRj92oiwm8ASxdupRx48ZRXl7OggULggrcCI0FBGNDoxV4QxNASktLueSSS8jKyuLDDz+kT58+MbQ29bHyUrGhUQbZGrL/li/RLjMzkzfffJP8/HwTtwMkW0AwWWiUAo82AeTQoUOMHz+eWbNmATBq1CjatWsXMzsbE7ZpYWxolC56NOO9LVu2kJOTw+rVqxk5st7UeiNKkiUgmEw0SoFHOt5btmwZ48aNo6ysjPnz53PxxRfH2kTDcIRG6aJHMt7btGkTI0aMoFWrVixfvtzEbSQVjbIHjyQBpGfPnjz77LOMGzfOxttG0uH67qIi8ihwOVANfAvcqKoloX4Xy91FDx8+zKRJk/j1r3/NkCFDYnINw3CK+nYXTQQXfZqqnqmq/YF3gQfdNGbr1q2cd955vP7666xatSpu17U8bCMWuO6iq+oBv4+tqLkxSFxZsWIFl19+OYcPH2bevHmMGTMmLte1PGwjViRCD46IPC4iW4HrcKkHX7lyJRdccAEZGRksW7YsbuIGK8xgxI6EKLqoqg+oanc8BRd/VU87MSu62L9/f+666y7y8/Pp27du0cRYYnnYRqyIi8BVdaSqnh7g9XatU/8OXFFPOy+p6kBVHZiVldVgu7777jtuv/12SkpKSEtL4/HHH8eJdiPFCjMYscJ1F11E/AuM5gDr43Hd4uJizj//fF544QUWL14cj0sGxfKwjVjhepANmCoivfFMk20GfhnrCxYUFJCTk8PBgwd55513uOyyy2J9yXpJpsIMRnLhusBVNahLHgsWL17MmDFj6NSpEwsXLuSMM86I5+WDYnnYRixw3UWPNwMGDODqq68mPz8/YcRtGLGiUQi8rKyMhx56iCNHjpCZmcmMGTPo0KGD22ZZcosRc1x30WNNSUkJY8eOpaCggEGDBrk+3vZhyS1GPEjpHvzzzz9n8ODBrF27lry8vIQRN1hyixEfUlbgCxYsYPjw4aSlpbF06VJycnLcNqkGltxixIOUFXivXr0YMWIEK1as4Mwzz3TbnDpYcosRD1JW4Keddhrz58+nY8eObpsSEEtuMeJBygfZEhVLbjHigQncRSy5xYg1KeuiG4ZhAjeMlMYEbhgpjAncMFIY16uqRouI7MKzvDQZaQ/sdtsIh0iVe0nm+zhRVQNWKklagSczIlIQrMxtspEq95Iq91Ebc9ENI4UxgRtGCmMCd4eX3DbAQVLlXlLlPmpgY3DDSGGsBzeMFMYE7hIiMk1E1ovIlyLyTxHJdNumSBCRi0WkSES+FpHJbtsTLSLSXUQWi8g6EVkjIr9x2yYnMRfdJURkFLBIVStF5I8Aqnqvy2aFhYikARuAHwPbgM+Aa1R1rauGRYGIdAY6q+rnItIGWAmMTcZ7CYT14C6hqv9W1Urvx+VANzftiZDBwNequlFVjwKv49kCOulQ1e2q+rn3/UFgHZAyS/xM4InBTcC/3DYiAroCW/0+byMFRCEiPYEBQL67ljiHrQePISLyPtApwFcP+PZlE5EHgEo8Gy8mCxLgWFKP9USkNTAHuLPWltZJjQk8hqjqyPq+F5GfApcBP9LkCoZsA7r7fe4GlLhkS4MRkXQ84p6tqnPdtsdJLMjmEiJyMfAUcIGqOrsXcowRkaZ4gmw/AorxBNmuVdU1rhoWBSIiwKvAXlW90217nMYE7hIi8jXQHNjjPbRcVWO+8aJTiMgY4GkgDXhFVR932aSoEJHhwMfAajwbYALcr6oL3LPKOUzghpHCWBTdMFIYE7hhpDAmcMNIYUzghpHCmMANI4UxgRtGCmMCTwJEZJOIlInIIRHZKSIzvKmV8bbjPK8NtV/VIvJKvO0xQmMCTx6yVbU1cBYwCPh9pA14M9CiRlU/VtXW/i9gPHAIT1Ze3GnoPaU6JvAkQ1WL8aw8Ox1ARI4XkZdFZLuIFIvIY9712ojIjSLyqYj8WUT2Ag/XOlYqIhtF5Fzv8a0i8q03Rz4kItIdzyKZ21T1qyDntBeRd73X2isiH4tIE9/vRWSuiOwSkT0i8qz3eBMR+b2IbPbaM1NEjvd+11NEVERuFpEtwCLv8aEistR7nVUicmED/swpgwk8yfCKagxQ6D30Kp7VaL3wLHR25B8AAAKJSURBVHUcBfzc7ydDgI1AB+Bxv2NfAu2Av+NZzz3I28b1wLOhhgDeBRr/AN5S1Vn1nHo3nsUpWUBH4H5AvQ+hd/FsXtETz3LT172/udH7GgGcDLQGnq3V7gXAD4DRItIVmA88BpwA3APMEZGAmwE0KlTVXgn+AjbhcYNL8QjieSADj2DKgQy/c68BFnvf3whsqdXWjcB//D6fgWepZ0e/Y3uA/iFsegbPIpPmIc6bArwN9Kp1/BxgF9A0wG8+wOMV+D73BirwrH7s6bX3ZL/v7wX+VquNhcBP3f5/5/bLxi/Jw1hVfd//gIicAaQD2z2LogCPV+ZfjMH/vY+dfu/LAFS19rGgPbiIXA1cC5ylquV+x3sAx0odqWeMPg14GPi318aXVHUqnuWmm/X7qjb+dKHmtlSb8Yi7Y5D7OhGYICLZfsfSgcXB7qGxYAJPbrbi6cHbBxEKOFyIQUR+gKeG+NWqWmNvOFXdQq0Hg3rKIN0N3C0ifYHFIvKZ1/YeItI0gO0leETroweeYchOvi9t5X9fW/H04Lc06OZSEBuDJzGquh34N/AnETnOG5w6RUQuiMX1RKQVnsIIf9Ewl1OKyGUi0su77voAUOV9rQC2A1NFpJWItBCRYd6fvQb8VkRO8sYC/gd4o56H2CwgW0RGi0iat60LRSSZ6tzFBBN48nMD0AyPa7wPeAvoHKNrXYEnsHVXgLnwYDXlTgXexxNDWAY8r6pLVLUKyMYT2NuCJxD3E+9vXgH+BnwEfAMcAe4IZpSqbsVT9PF+POP6rUAu9u/b1oMbRirT6J9whpHKmMANI4UxgRtGCmMCN4wUxgRuGCmMCdwwUhgTuGGkMCZww0hhTOCGkcL8f9CRu1jVslhPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 252x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# How do the permutation results compare to our original t-test?\n",
    "# An easy way to check would be to compare the permutation z-score to the t-statistics \n",
    "# (with large enough samples, these statistics are essentially the same).\n",
    "\n",
    "# Compute z-score for each electrode.\n",
    "zs = []\n",
    "for i in range(len(diff_mu)):\n",
    "    zs.append((diff_mu[i]-np.mean(null_mus[:, i]))/np.std(null_mus[:, i]))\n",
    "\n",
    "# Plot results.\n",
    "plt.figure(figsize=(3.5, 3))\n",
    "ax = plt.subplot(111)\n",
    "plt.plot(zs, ts, 'o')\n",
    "plt.plot([-3.5, 3.5], [-3.5, 3.5], 'k--')\n",
    "plt.xlabel('Perm Z-score', fontsize=12)\n",
    "plt.ylabel('T-stat', fontsize=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The t-stats and z-scores look very similar, meaning that the independence violation probably wasn't a big deal. There is a hint that the t-stats might have been slightly overstating our effects at the low end of the range. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise: Generalize the approach above to run a permutation t-test. How do the results compare to the original parametric t-test? How do they compare to the difference in means permutation test?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise: How does the correlation of between the t-test and the permutation z-scores change as a function of the number of random shuffles?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise: Can you imagine a way to extend permutation tests to help correct for multiple comparisons?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cluster testing "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mne.stats import permutation_cluster_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "environmentname",
   "language": "python",
   "name": "environmentname"
  },
  "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}