t-flora/cs146-10.1-pre-class-work.ipynb

## cs146-10.1-pre-class-work.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pre-class work\n",
    "\n",
    "These notes work through the motivation behind testing a model we have seen before. Use the notes to refresh your memory of this particular model, make sure you understand the choice of test statistic used below, and write Python code to compute the p-value for the test statistic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model and dataset\n",
    "\n",
    "In a previous class session, we modeled the dataset shown below using a **normal likelihood with unknown mean and variance** and a **conjugate normal-inverse-gamma prior** over the parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsUAAAF1CAYAAAAA6ZfwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAY8UlEQVR4nO3df9DlV10f8PeHXWLlZzpmKZgfJMVkaHCg4jZAbSu27DT8mkCl0yDC0JHGaIM1lWrqIJVxLFIdVtoE11QziDpGCpSudGm6U4uogN2FIrKBrEsAswbIEgkhioSFT/+4d+vl4dl97m7uvU825/WaeWbv9/s9957PPXtm573nOfd+q7sDAAAje9BmFwAAAJtNKAYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMPSFV1oKqevtl1bKaqen5V3VZV91TVt83R/l1V9bJV1AZwfyMUA6edqvpEVT1jzbmXVtXvHTvu7id097s2eJ3zq6qrauuSSt1sP5fkqu5+WHf/30W+8Hp/B8uwqn4AhGKAJbkfhO3HJjmwyTUAnBaEYuABaXaFsaouqar9VXV3VX2mql43bfbu6Z93TbcYPK2qHlRVr6yqT1bVHVX1pqp65MzrvmR67c6q+ok1/fxkVb2lqn6tqu5O8tJp3++tqruq6lNVdW1VnTHzel1VP1hVf1xVX6iqn6qqx02fc3dVvXm2/Zr3uG6tVfUNVXVPki1J/rCqPnac5++oqo9W1eer6tokNXPtcVX129P3+dmq+vWqOnN67VeTnJfkt6bj9qPT8/+lqj49fb13V9UTZl7vWVV18/Q9/mlVvWLm2nOq6oPTMXpPVT3xRP0ALINQDIzg9Ule392PSPK4JG+env8H0z/PnG4xeG+Sl05/vivJ30zysCTXJklVXZzkDUlelOQxSR6Z5Ow1fV2W5C1Jzkzy60m+kuTqJGcleVqSf5TkB9c859Ik357kqUl+NMn10z7OTfKtSV54nPe1bq3d/aXufti0zZO6+3Frn1hVZyV5a5JXTmv7WJLvmG2S5DVJvjnJ35rW8pNJ0t0vTvInSZ47Hbf/MH3OO5NcmORRST4wff/H/HKS7+/uh0/f029P63hykhuSfH+Sb0ryi0l2V9U3nKAfgIUTioHT1dunK4t3VdVdmYTV4/lykm+pqrO6+57uft8J2r4oyeu6+9buvifJv01y+XQrxAuS/FZ3/15335vkVUl6zfPf291v7+6vdvcXu/v93f2+7j7a3Z/IJPR955rnvLa77+7uA0k+nOR/Tvv/fCZB83gfkjtRrRt5VpKbu/st3f3lJD+f5NPHLnb3oe7eOw3YR5K8bp26v0Z339DdX+juL2USoJ80s8r+5SQXV9Ujuvtz3f2B6fl/keQXu/sPuvsr3f0rSb6UyX8QAFZGKAZOV8/r7jOP/eTrV19nfV+Si5J8tKr2VdVzTtD2m5N8cub4k0m2Jvkb02u3HbvQ3X+R5M41z79t9qCqLqqqd0y3Fdyd5N9nsjI76zMzj7+4zvHDsr4T1bqRte+lZ4+r6lFVdeN0q8PdSX5tnboz035LVf1MVX1s2v4T00vHnvPdmQTxT1bV71TV06bnH5vkR9b8B+fcaX0AKyMUAw943f3H3f3CTH6t/9okb6mqh+brV3mT5PZMgtox5yU5mklQ/VSSc45dqKpvzORX/l/T3ZrjX0jy0SQXTrdv/Hhm9u7eRyeqdSOfyiR8JkmqqmaPM9k60UmeOK37e/O1da99n9+TydaRZ2SyreT8Yy+dJN29r7svy+Tv4O35qy0styX56dn/4HT3Q7r7N47TD8BSCMXAA15VfW9Vbevurya5a3r6K0mOJPlqJvtxj/mNJFdX1QVV9bBMVnZ/s7uPZrJX+LlV9XenH357dTYOuA9PcneSe6rq8Ul+YFHva4NaN/Lfkzyhqv7JdLvFDyV59Jq678nkQ4hnJ/k3a57/mXztuD08k20PdyZ5yLSWJElVnVFVL6qqR063atydyfgnyX9OcmVVPaUmHlpVz66qhx+nH4ClEIqBEVya5MD0Gxlen+Ty7v7L6faHn07y+9Nf3T81kw99/Wom30zx8SR/meTlSTLd8/vyJDdmstL6hSR3ZBIGj+cVmayifiGTAPibC3xfx611I9392ST/NMnPZBJkL0zy+zNNXp3kyUk+n0mAftual3hNkldOx+0VSd6UyfaNP01yc5K1+7ZfnOQT060VV2ay8pzu3p/JvuJrk3wuyaFMPjx4vH4AlqIm28gAOFnT1dm7Mtka8fFNLgeA+8BKMcBJqKrnVtVDpnuSfy7JH+WvPlQGwGlKKAY4OZdl8gG32zPZcnB5+5UbwGnP9gkAAIZnpRgAgOEJxQAADG+eW4EuxVlnndXnn3/+ZnUPAMAg3v/+93+2u7edqM2mheLzzz8/+/fv36zuAQAYRFV9cqM2tk8AADA8oRgAgOEJxQAADE8oBgBgeEIxAADDE4oBABieUAwAwPCEYgAAhicUAwAwPKEYAIDhCcUAAAxPKAYAYHhCMQAAw9u62QUAnKqdew9uWt9X77ho0/oGYPGsFAMAMDyhGACA4QnFAAAMTygGAGB4QjEAAMMTigEAGJ5QDADA8IRiAACGJxQDADA8oRgAgOEJxQAADE8oBgBgeEIxAADDE4oBABieUAwAwPCEYgAAhicUAwAwPKEYAIDhCcUAAAxPKAYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMDyhGACA4c0Viqvq0qq6paoOVdU1x2nz9Kr6YFUdqKrfWWyZAACwPFs3alBVW5Jcl2RHksNJ9lXV7u6+eabNmUnekOTS7v6TqnrUkuoFAJZk596Dm9b31Tsu2rS+IZlvpfiSJIe6+9buvjfJjUkuW9Pme5K8rbv/JEm6+47FlgkAAMszTyg+O8ltM8eHp+dmXZTkr1fVu6rq/VX1kvVeqKquqKr9VbX/yJEjp1YxAAAs2DyhuNY512uOtyb59iTPTvKPk/xEVX3d70G6+/ru3t7d27dt23bSxQIAwDJsuKc4k5Xhc2eOz0ly+zptPtvdf57kz6vq3UmelGTzNicBAMCc5lkp3pfkwqq6oKrOSHJ5kt1r2vy3JH+/qrZW1UOSPCXJRxZbKgAALMeGK8XdfbSqrkpyU5ItSW7o7gNVdeX0+q7u/khV/Y8kH0ry1SS/1N0fXmbhAACwKPNsn0h370myZ825XWuOfzbJzy6uNAAAWA13tAMAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMDyhGACA4QnFAAAMTygGAGB4QjEAAMMTigEAGJ5QDADA8IRiAACGJxQDADA8oRgAgOEJxQAADE8oBgBgeEIxAADDE4oBABieUAwAwPCEYgAAhrd1swsAgPujnXsPblrfV++4aNP6hlFZKQYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMDyhGACA4QnFAAAMTygGAGB4QjEAAMMTigEAGJ5QDADA8IRiAACGJxQDADA8oRgAgOHNFYqr6tKquqWqDlXVNetcf3pVfb6qPjj9edXiSwUAgOXYulGDqtqS5LokO5IcTrKvqnZ3981rmv5udz9nCTUCwFB27j242SXAcOZZKb4kyaHuvrW7701yY5LLllsWAACszjyh+Owkt80cH56eW+tpVfWHVfXOqnrCQqoDAIAV2HD7RJJa51yvOf5Aksd29z1V9awkb09y4de9UNUVSa5IkvPOO+/kKgUAgCWZZ6X4cJJzZ47PSXL7bIPuvru775k+3pPkwVV11toX6u7ru3t7d2/ftm3bfSgbAAAWZ55QvC/JhVV1QVWdkeTyJLtnG1TVo6uqpo8vmb7unYsuFgAAlmHD7RPdfbSqrkpyU5ItSW7o7gNVdeX0+q4kL0jyA1V1NMkXk1ze3Wu3WAAAwP3SPHuKj22J2LPm3K6Zx9cmuXaxpQEAwGq4ox0AAMMTigEAGJ5QDADA8ObaUwxwIm5JC8DpzkoxAADDE4oBABieUAwAwPCEYgAAhicUAwAwPKEYAIDhCcUAAAxPKAYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMDyhGACA4QnFAAAMTygGAGB4QjEAAMMTigEAGJ5QDADA8IRiAACGJxQDADA8oRgAgOEJxQAADE8oBgBgeEIxAADDE4oBABieUAwAwPC2bnYBAHAiO/ce3OwSgAFYKQYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMLy5QnFVXVpVt1TVoaq65gTt/k5VfaWqXrC4EgEAYLk2DMVVtSXJdUmemeTiJC+sqouP0+61SW5adJEAALBM86wUX5LkUHff2t33JrkxyWXrtHt5krcmuWOB9QEAwNLNE4rPTnLbzPHh6bn/r6rOTvL8JLsWVxoAAKzGPKG41jnXa45/PsmPdfdXTvhCVVdU1f6q2n/kyJE5SwQAgOXaOkebw0nOnTk+J8nta9psT3JjVSXJWUmeVVVHu/vts426+/ok1yfJ9u3b1wZrAADYFPOE4n1JLqyqC5L8aZLLk3zPbIPuvuDY46p6Y5J3rA3EAABwf7VhKO7uo1V1VSbfKrElyQ3dfaCqrpxet48YAIDT2jwrxenuPUn2rDm3bhju7pfe97IAAGB13NEOAIDhzbVSDMDX2rn34Kb0e/WOizal3816vwCrYqUYAIDhCcUAAAxPKAYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMDyhGACA4QnFAAAMTygGAGB4QjEAAMMTigEAGJ5QDADA8IRiAACGJxQDADA8oRgAgOEJxQAADE8oBgBgeEIxAADDE4oBABieUAwAwPC2bnYBwGLs3Htws0tgBfw9AyyHlWIAAIYnFAMAMDyhGACA4QnFAAAMTygGAGB4QjEAAMMTigEAGJ5QDADA8IRiAACGJxQDADA8oRgAgOEJxQAADE8oBgBgeHOF4qq6tKpuqapDVXXNOtcvq6oPVdUHq2p/Vf29xZcKAADLsXWjBlW1Jcl1SXYkOZxkX1Xt7u6bZ5r9ryS7u7ur6olJ3pzk8csoGAAAFm2eleJLkhzq7lu7+94kNya5bLZBd9/T3T09fGiSDgAAnCY2XClOcnaS22aODyd5ytpGVfX8JK9J8qgkz17vharqiiRXJMl55513srUCACzUzr0HN6Xfq3dctCn9cnzzrBTXOue+biW4u/9rdz8+yfOS/NR6L9Td13f39u7evm3btpMqFAAAlmWeUHw4ybkzx+ckuf14jbv73UkeV1Vn3cfaAABgJeYJxfuSXFhVF1TVGUkuT7J7tkFVfUtV1fTxk5OckeTORRcLAADLsOGe4u4+WlVXJbkpyZYkN3T3gaq6cnp9V5LvTvKSqvpyki8m+WczH7wDAID7tXk+aJfu3pNkz5pzu2YevzbJaxdbGgAArIY72gEAMDyhGACA4QnFAAAMTygGAGB4QjEAAMMTigEAGJ5QDADA8IRiAACGJxQDADA8oRgAgOEJxQAADE8oBgBgeEIxAADDE4oBABieUAwAwPCEYgAAhicUAwAwPKEYAIDhCcUAAAxPKAYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMDyhGACA4QnFAAAMb+tmFwDLsHPvwU3r++odF21a3wDAqbFSDADA8IRiAACGJxQDADA8oRgAgOEJxQAADE8oBgBgeEIxAADDE4oBABieUAwAwPCEYgAAhjdXKK6qS6vqlqo6VFXXrHP9RVX1oenPe6rqSYsvFQAAlmPrRg2qakuS65LsSHI4yb6q2t3dN880+3iS7+zuz1XVM5Ncn+QpyygY7u927j242SUAnHb828lmm2el+JIkh7r71u6+N8mNSS6bbdDd7+nuz00P35fknMWWCQAAyzNPKD47yW0zx4en547n+5K8c70LVXVFVe2vqv1HjhyZv0oAAFiieUJxrXOu121Y9V2ZhOIfW+96d1/f3du7e/u2bdvmrxIAAJZowz3FmawMnztzfE6S29c2qqonJvmlJM/s7jsXUx4AACzfPCvF+5JcWFUXVNUZSS5Psnu2QVWdl+RtSV7c3XbKAwBwWtlwpbi7j1bVVUluSrIlyQ3dfaCqrpxe35XkVUm+KckbqipJjnb39uWVDQAAizPP9ol0954ke9ac2zXz+GVJXrbY0gAAYDXc0Q4AgOEJxQAADE8oBgBgeEIxAADDE4oBABieUAwAwPCEYgAAhicUAwAwPKEYAIDhCcUAAAxPKAYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMDyhGACA4QnFAAAMTygGAGB4QjEAAMMTigEAGJ5QDADA8IRiAACGJxQDADA8oRgAgOEJxQAADG/rZhcAADCanXsPblrfV++4aNP6vj+zUgwAwPCEYgAAhicUAwAwPKEYAIDhCcUAAAxPKAYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMLy5QnFVXVpVt1TVoaq6Zp3rj6+q91bVl6rqFYsvEwAAlmfrRg2qakuS65LsSHI4yb6q2t3dN880+7MkP5TkecsoEgAAlmmeleJLkhzq7lu7+94kNya5bLZBd9/R3fuSfHkJNQIAwFLNE4rPTnLbzPHh6bmTVlVXVNX+qtp/5MiRU3kJAABYuHlCca1zrk+ls+6+vru3d/f2bdu2ncpLAADAws0Tig8nOXfm+Jwkty+nHAAAWL15QvG+JBdW1QVVdUaSy5PsXm5ZAACwOht++0R3H62qq5LclGRLkhu6+0BVXTm9vquqHp1kf5JHJPlqVf1wkou7++7llQ4AAIuxYShOku7ek2TPmnO7Zh5/OpNtFQAAcNpxRzsAAIY310oxp7+dew9udgkAAPdbVooBABieUAwAwPCEYgAAhicUAwAwPKEYAIDhCcUAAAxPKAYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMDyhGACA4QnFAAAMTygGAGB4QjEAAMPbutkFAACwOjv3HtyUfq/ecdGm9DsvK8UAAAxPKAYAYHhCMQAAwxOKAQAYnlAMAMDwhvz2ic361CUAAPdPVooBABieUAwAwPCEYgAAhicUAwAwPKEYAIDhCcUAAAxPKAYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMLy5QnFVXVpVt1TVoaq6Zp3rVVX/cXr9Q1X15MWXCgAAy7FhKK6qLUmuS/LMJBcneWFVXbym2TOTXDj9uSLJLyy4TgAAWJp5VoovSXKou2/t7nuT3JjksjVtLkvypp54X5Izq+oxC64VAACWYp5QfHaS22aOD0/PnWwbAAC4X9o6R5ta51yfQptU1RWZbK9Iknuq6pY5+j8ZZyX57IJfk+Mz3qtlvFfLeK+W8V4dY71axnvqX6+mm+ON92M3euI8ofhwknNnjs9JcvsptEl3X5/k+jn6PCVVtb+7ty/r9flaxnu1jPdqGe/VMt6rY6xXy3iv1n0Z73m2T+xLcmFVXVBVZyS5PMnuNW12J3nJ9Fsonprk8939qVMpCAAAVm3DleLuPlpVVyW5KcmWJDd094GqunJ6fVeSPUmeleRQkr9I8s+XVzIAACzWPNsn0t17Mgm+s+d2zTzuJP9ysaWdkqVtzWBdxnu1jPdqGe/VMt6rY6xXy3iv1imPd03yLAAAjMttngEAGN5pGYrddnq15hjvp1fV56vqg9OfV21GnQ8EVXVDVd1RVR8+znVze4HmGG9ze0Gq6tyq+t9V9ZGqOlBV/2qdNub3gsw53ub3glTVX6uq/1NVfzgd71ev08b8XpA5x/uk5/dce4rvT2ZuO70jk6+C21dVu7v75plms7edfkomt51+yqprfSCYc7yT5He7+zkrL/CB541Jrk3ypuNcN7cX64058Xgn5vaiHE3yI939gap6eJL3V9Ve/3YvzTzjnZjfi/KlJP+wu++pqgcn+b2qeuf0Lr/HmN+LM894Jyc5v0/HlWK3nV6tecabBenudyf5sxM0MbcXaI7xZkG6+1Pd/YHp4y8k+Ui+/s6n5veCzDneLMh0zt4zPXzw9Gfth7bM7wWZc7xP2ukYit12erXmHcunTX+N8c6qesJqShuSub165vaCVdX5Sb4tyR+suWR+L8EJxjsxvxemqrZU1QeT3JFkb3eb30s0x3gnJzm/T8dQvLDbTjOXecbyA0ke291PSvKfkrx92UUNzNxeLXN7warqYUnemuSHu/vutZfXeYr5fR9sMN7m9wJ191e6+29nclffS6rqW9c0Mb8XaI7xPun5fTqG4oXddpq5bDiW3X33sV9jTL/T+sFVddbqShyKub1C5vZiTff+vTXJr3f329ZpYn4v0EbjbX4vR3ffleRdSS5dc8n8XoLjjfepzO/TMRS77fRqbTjeVfXoqqrp40symVd3rrzSMZjbK2RuL850HH85yUe6+3XHaWZ+L8g8421+L05VbauqM6ePvzHJM5J8dE0z83tB5hnvU5nfp923T7jt9GrNOd4vSPIDVXU0yReTXN7uCnNKquo3kjw9yVlVdTjJv8vkAwTm9hLMMd7m9uJ8R5IXJ/mj6T7AJPnxJOcl5vcSzDPe5vfiPCbJr0y/selBSd7c3e+QTZZmnvE+6fntjnYAAAzvdNw+AQAACyUUAwAwPKEYAIDhCcUAAAxPKAYAYHhCMQAAwxOKAQAYnlAMAMDw/h97Uk9nU4xCcQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Load data: read the particle sizes (in nanometers) from a CSV file.\n",
    "# Log-transform the data so we can model it using a normal likelihood.\n",
    "data = np.log(np.loadtxt('hrtem.csv'))\n",
    "\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.hist(data, bins=20, density=True, alpha=0.5)\n",
    "plt.title('Histogram of dataset')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model was as follows.\n",
    "\n",
    "* Data: $\\{y_i\\}$\n",
    "* Parameters: mean $x$, variance $\\sigma^2$\n",
    "* Likelihood: $y_i \\sim \\text{Normal}(x, \\sigma^2)$\n",
    "* Prior: $(\\mu,\\sigma^2) \\sim \\text{Normal-Inverse-Gamma}(\\mu_0,\\nu_0,\\alpha_0,\\beta_0)$\n",
    "\n",
    "The prior hyperparameter values are given below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posterior hyperparameters:\n",
      "  μ₀ = 2.3\n",
      "  ν₀ = 0.1\n",
      "  α₀ = 2\n",
      "  β₀ = 5\n"
     ]
    }
   ],
   "source": [
    "mu_0 = 2.3\n",
    "nu_0 = 0.1\n",
    "alpha_0 = 2\n",
    "beta_0 = 5\n",
    "\n",
    "print('Posterior hyperparameters:')\n",
    "print('  μ₀ =', mu_0)\n",
    "print('  ν₀ =', nu_0)\n",
    "print('  α₀ =', alpha_0)\n",
    "print('  β₀ =', beta_0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the prior is conjugate to the likelihood, the posterior is also a Normal-Inverse-Gamma distribution. The posterior hyperparameters are calculated below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posterior hyperparameters:\n",
      "  μ₁ = 1.892401480510831\n",
      "  ν₁ = 500.1\n",
      "  α₁ = 252.0\n",
      "  β₁ = 124.45079772182757\n"
     ]
    }
   ],
   "source": [
    "# Sufficient statistics of the data\n",
    "s0 = len(data)\n",
    "s1 = sum(data)\n",
    "s2 = sum(data ** 2)\n",
    "\n",
    "# Posterior parameters\n",
    "mu_1 = (nu_0 * mu_0 + s1) / (nu_0 + s0)\n",
    "nu_1 = nu_0 + s0\n",
    "alpha_1 = alpha_0 + s0 / 2\n",
    "beta_1 = beta_0 + s2/2 - s1**2 / (2*s0) + s0*nu_0/(nu_1 + s0) * (s1/s0 - mu_0)**2/2\n",
    "\n",
    "print('Posterior hyperparameters:')\n",
    "print('  μ₁ =', mu_1)\n",
    "print('  ν₁ =', nu_1)\n",
    "print('  α₁ =', alpha_1)\n",
    "print('  β₁ =', beta_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**One criticism against this model** is that it looks like the data distribution might be bimodal rather than unimodal, which would make a normal likelihood inappropriate. Let’s design a test statistic to check whether this bimodal appearance is really statistically significant or not."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test statistic\n",
    "\n",
    "Our test statistic is the proportion of data values that lie within one standard deviation of the mean of the dataset.\n",
    "\n",
    "So, given a dataset $\\{y_i\\}$, we compute the sample mean $\\bar{\\mu}$ and sample standard deviation $\\bar{\\sigma}$. We then count the proportion of data values that lie in the range $[\\bar{\\mu}-\\bar{\\sigma}, \\bar{\\mu}+\\bar{\\sigma}]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_statistic(data):\n",
    "    mu = np.mean(data)  # sample mean\n",
    "    sigma = np.std(data, ddof=1)  # sample standard deviation\n",
    "    return np.mean((data > mu - sigma) & (data < mu + sigma))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a normally distributed dataset, we expect approximately 68% of the data to lie within one standard deviation of the mean. For a bimodal dataset with equal probability mass in the two modes, we expect far fewer values to lie in this range since there is a gap between the modes.\n",
    "\n",
    "Here is an example to demonstrate why."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs8AAAF1CAYAAAAXywc5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjlklEQVR4nO3df7RdZX3n8fenCegoWEADIglCnaCmrSIrRVvtVEtjgVqDndpCLUaKjcyYUahOm2p/2Fpba9FYl5QM1ihaFXUUTW0UM4y0oxVXAkUUaUJM+RESSEQQKRaIfOePva8eLie5z03uzUnC+7XWWffsZz+/ds4NfPKc5+yTqkKSJEnSxH5k1BOQJEmS9hWGZ0mSJKmR4VmSJElqZHiWJEmSGhmeJUmSpEaGZ0mSJKmR4VnSHpfkZUk+P4Jxn5vkhiT3JDmtof77k/zZHpgaSZYn+cM9MdaelOTGJN9L8sFRz2U6JXlU/3v1wJ76nZE0GoZnaT8xEFLuSXJ7kvclOWgvmNcxSSrJzLGyqvpQVb1wBNP5U+DdVXVQVX1qKjtOckWSV+5q+6o6p6rePJVz2ov8clWdOepJTKequq+qDgI+NJl2Sc5LcluS7yRZkeRRO6k7I8mfJdmc5LtJ/iXJIS19JVmSZG2S+5K8fxcuUVLP8CztX365/x/4CcBPAX8wvsJgiJ1ue3KsRk8Grhv1JKZTOv63fRcleX6SK/bQWL8ILAVOAo4Bfgz4k500+RPgZ4CfBh4HnAn8R2Nfm4E/A1ZM4SVIj0j+B1baD1XVrcBngZ8A6Fd+X53kBuCGvuy3k2xI8u0kK5M8aax9X/81STYm+VaSvxoLZEl+JMkfJLkpydYkH0jyo/25sVXms5PcDPxf4J/6bu/qV8V/OskrknxxYLyfSbKmXzFbk+RnBs5dkeTNSb7Ur7Z9PskTdnTtO7quJN+kCxR/38/jYSt8SZ6V5Op+nI8Cjx44d2iSzyTZluTO/vns/txbgJ8F3t33/e6+/K+T3JLk7iRXJfnZncz7B1tE+gC3Kcnr+j/jLUnO2knbK5K8JcmXgHuBH0tyVpLr+2vZmORVA/V32n+Sxyf5+37ea/rVzsHX62lJVvd/xuuS/NqO5jZkrg97JyINq/bptkVc0M/1+30fleRTDWM+5PetL6sk/7l13n2bOUlW9dddA49zJ9PPgEXAe6vquqq6E3gz8IodjH0ocC7w21V1U3W+XlX/0dJXVX2yf7fljl2cq6Se4VnaDyWZA5wK/MtA8WnAs4F5SX4e+Avg14AjgZuAS8Z18xJgPt0q9kLgt/ryV/SPF9CF0YOAd49r+3PA04FfBP5LX3ZIv13iy+PmehjwD8C7gMcD7wD+IcnjB6r9BnAWcDhwIPD6HVz3Dq+rqp4C3Ey/Ol9V941reyDwKeCDwGHAx4H/OlDlR4D30a1eHw18b+y6q+qNwP8DlvR9L+nbrAGO7/v7MPDxJI+mzROBHwWOAs4GLugD1I6cCSwGDu6veyvwIroVyrOAZUlOaOz/AuDf+zqL+gcASR4LrO6v53DgDOBvkvx443XtqiXAL/DDVdePA1cA503zuIP+hu51P7p//CtwPjvYqpHk6CR3JTl6B/39OPDVgeOvAkeM+90f85PAduBX023NWJ/k1bvYl6TdYHiW9i+fSnIX8EXgH4E/Hzj3F1X17ar6HvAyYEVVXd2HyN8HfjrJMQP1/7KvfzPwTrqQRN/2HVW1saru6duenodu0XhTVf17P9ZEfgm4oao+WFXbq+ojdKHklwfqvK+q1vf9fYwukA7Tcl078hzgAOCdVfVAVf1vuvALQFXdUVWfqKp7q+q7wFvo/pGwQ1X1d3277VX1duBRwFMb5gLwAPCn/VxWAfdM0Pb9/arj9r7NP1TVN/sVyn8EPk+3Or7T/pPMoPtHwx/31/oN4OKBdi8Cbqyq9/VjXQ18AvjVxuvaVS+me21urKp/p9uS9Dzg1mkeF/jBFqRTgD+oqnuq6hZgGXBiVW0b1qaqbq6qQ/q/Q8McBHxn4Hjs+cFD6s6m+8fOccCxdH/eb0qyYBf6krQbDM/S/uW0/n/WT66q/z4uvN4y8PxJdKuTAPQh+A66Vchh9W/q2zysbf98JnDEDtpOZHx/Y30OzuW2gef30gWFCfvawXXtbB63VlWNmwcASR6T5H+l265yN912lEP6sDlUvy3i+nTbUe6iCz873HIyzh1VtX3geGfXDeP+zJOckuTKfovBXXTvRAyOvaP+Z9G9noP9DT5/MvDsfkX1rr7vl9GtUu+2JN/st76MPW7sTx1O987BmLHfu1lTMObSgWv5DPC8cdcH3bsiM4bM4UnsunvoVtHHjD3/7pC6Y3+X/7SqvldV19K9q3LqLvQlaTcYnqVHjsFQuJkuBAE/eCv+8Tx0FW/OwPOj+zYPa9uf2w7cvoOxBp8PM76/sT53ZUWx5bp2ZAtwVJKMm8eY19Gt/D67qh7HD7ejjNV/yHX2+5t/j24LyaFVdQjdauBg/1PpB+On28/9CbotBUf0Y69qHHsb3es5e6Bs8HfhFuAf+3+kjT0Oqqr/Nsn5Dm5fOeQHF1H1lL6/sccx/alNdCuuY47l4b93TeOl36M/MOZbx66FbmX9i4PX11fbBtw3ZA6bGscf5jrgmQPHzwRur6ph+5KvHZvuFPQlaTcYnqVHpg8DZyU5vg9afw58papuHKjzP9N9SG4O8Frgo335R4Dzkhyb7lZ4fw58dNwq5qBtwIN0+6OHWQUcl+Q3ksxM8uvAPLoVwOm4rh35Ml0Ye00/j18BThw4fzDd6t9d/T7tPx7X/nYeeo0H9/1tA2Ym+SMeujI4nQ6k2yKyDdie5BSg6daAVfV94JN0WwIek+RpwMsHqnyG7vU6M8kB/eOnkjx9knM8O92HTxfQ7ec9OMkBO6n/QeDcJD827Pcu3Qcu37+T9s9IcmL/TsHYPunDWidbVQ/S/X69JcnB/Vag32GSt6Yb5wN0fw7z+v3mfwC8fwfjf5NuX/0b+w9PPh34dX7492SnffW/04+mWz2fkeTR2fvuhiPtEwzP0iNQVV0O/CHd6uQW4CnA6eOqfRq4CriG7gN97+3LV9AFmX8C/o3uVln/Yydj3Uu3P/hL/dvgzxl3/g661b7X0W2x+F3gRVX1rWm6rh21vR/4FboPQ95JF0w+OVDlncB/Ar4FXAl8blwXf033Ya47k7wLuIzujifr6d7e/w8mt51ll/V7sl9Dtz/8TroPXK6cRBdL6LaY3Eb3Wn+EbtV1rO8X0v25bu7r/CVdWJ+Mp9CF+9+ju6XaH9Ht592RD9L97n2B7vfu3n6eY+YAX9pJ+2/Q3aptK927ER/nh/8gbHUucBewke73/+/Yya3f+g8M3rOjDwxW1eeAt9Fd0039448H2n82yRsGmpxB987KHXR/J/+w/52fsC+6MP09utvZ/Wb//GG3spQ0sTx0e58kdbfxAuZW1YZRz0Wjl+QvgSdW1aIJKz+87Tq6O59cWlWL+hXbfwMO2Mm7FZMd40C6u0s8o6oeGHL+FcArq+p5UzHeDubwKLp3Hw4A3lZVO7tfs6R9mG/ZSJIeot+qcSDwNbov2zkb2KVvT6yq1ruL7LL+XYPJbhuZ6jncx8DebUn7L8OzJGm8g+m2ajyJbpvD2+m28UjSI57bNiRJkqRGfmBQkiRJamR4liRJkhrtU3uen/CEJ9Qxxxwz6mlI+6R169YB8NSnTvvntyTtg/xvhPRQV1111beq6mHfYrpPhedjjjmGtWvXjnoa0j7p+c9/PgBXXHHFSOchae/kfyOkh0py07Byt21IkiRJjQzPkiRJUiPDsyRJktTI8CxJkiQ1MjxLkiRJjQzPkiRJUiPDsyRJktTI8CxJkiQ1MjxLkiRJjQzPkiRJUiPDsyRJktTI8CxJkiQ1MjxLkiRJjWaOegKSJO0vlq1eP5Jxz1tw3EjGlR6JXHmWJEmSGhmeJUmSpEaGZ0mSJKmR4VmSJElqZHiWJEmSGhmeJUmSpEaGZ0mSJKmR4VmSJElqZHiWJEmSGhmeJUmSpEaGZ0mSJKmR4VmSJElqZHiWJEmSGjWF5yQnJ1mXZEOSpUPOvyzJtf3jn5M8c6K2SQ5LsjrJDf3PQ6fmkiRJkqTpMWF4TjIDuAA4BZgHnJFk3rhq/wb8XFU9A3gzcFFD26XA5VU1F7i8P5YkSZL2Wi0rzycCG6pqY1XdD1wCLBysUFX/XFV39odXArMb2i4ELu6fXwyctstXIUmSJO0BLeH5KOCWgeNNfdmOnA18tqHtEVW1BaD/eXjLhCVJkqRRmdlQJ0PKamjF5AV04fl5k227w8GTxcBigKOPPnoyTSVJkqQp1bLyvAmYM3A8G9g8vlKSZwB/Cyysqjsa2t6e5Mi+7ZHA1mGDV9VFVTW/qubPmjWrYbqSJEnS9GgJz2uAuUmOTXIgcDqwcrBCkqOBTwJnVtX6xrYrgUX980XAp3f9MiRJkqTpN+G2jaranmQJcBkwA1hRVdclOac/vxz4I+DxwN8kAdjerxYPbdt3/VbgY0nOBm4GXjrF1yZJkiRNqZY9z1TVKmDVuLLlA89fCbyytW1ffgdw0mQmK0mSJI2S3zAoSZIkNWpaeZa0f1q2ev3ElabBeQuOG8m4kiTtLleeJUmSpEaGZ0mSJKmR4VmSJElqZHiWJEmSGhmeJUmSpEaGZ0mSJKmR4VmSJElqZHiWJEmSGhmeJUmSpEaGZ0mSJKmR4VmSJElqZHiWJEmSGhmeJUmSpEaGZ0mSJKmR4VmSJElqZHiWJEmSGhmeJUmSpEaGZ0mSJKmR4VmSJElqZHiWJEmSGhmeJUmSpEaGZ0mSJKmR4VmSJElqNHPUE5AkSbtn2er1u93Hpju/t0t9nbfguN0eW9qXuPIsSZIkNTI8S5IkSY2awnOSk5OsS7IhydIh55+W5MtJ7kvy+oHypya5ZuBxd5Jz+3NvSnLrwLlTp+yqJEmSpGkw4Z7nJDOAC4AFwCZgTZKVVfWNgWrfBl4DnDbYtqrWAccP9HMrcOlAlWVVdf5uzF+SJEnaY1pWnk8ENlTVxqq6H7gEWDhYoaq2VtUa4IGd9HMS8M2qummXZytJkiSNUEt4Pgq4ZeB4U182WacDHxlXtiTJtUlWJDl0WKMki5OsTbJ227ZtuzCsJEmSNDVablWXIWU1mUGSHAi8GPj9geILgTf3fb0ZeDvwWw8bqOoi4CKA+fPnT2pcaV8xFbeZmsiu3oZKkiT9UMvK8yZgzsDxbGDzJMc5Bbi6qm4fK6iq26vq+1X1IPAeuu0hkiRJ0l6rJTyvAeYmObZfQT4dWDnJcc5g3JaNJEcOHL4E+Pok+5QkSZL2qAm3bVTV9iRLgMuAGcCKqrouyTn9+eVJngisBR4HPNjfjm5eVd2d5DF0d+p41biu35bkeLptGzcOOS9JkiTtVZq+nruqVgGrxpUtH3h+G912jmFt7wUeP6T8zEnNVJIkSRoxv2FQkiRJamR4liRJkhoZniVJkqRGhmdJkiSpkeFZkiRJamR4liRJkhoZniVJkqRGhmdJkiSpkeFZkiRJamR4liRJkhoZniVJkqRGhmdJkiSpkeFZkiRJamR4liRJkhoZniVJkqRGM0c9AUmStO9atnr9SMY9b8FxIxlXcuVZkiRJamR4liRJkhoZniVJkqRGhmdJkiSpkeFZkiRJamR4liRJkhoZniVJkqRGhmdJkiSpkeFZkiRJamR4liRJkhoZniVJkqRGhmdJkiSp0cyWSklOBv4amAH8bVW9ddz5pwHvA04A3lhV5w+cuxH4LvB9YHtVze/LDwM+ChwD3Aj8WlXduXuXI2lfsGz1+pGNfd6C40Y2tiRp3zfhynOSGcAFwCnAPOCMJPPGVfs28BrgfIZ7QVUdPxace0uBy6tqLnB5fyxJkiTttVq2bZwIbKiqjVV1P3AJsHCwQlVtrao1wAOTGHshcHH//GLgtEm0lSRJkva4lvB8FHDLwPGmvqxVAZ9PclWSxQPlR1TVFoD+5+HDGidZnGRtkrXbtm2bxLCSJEnS1GoJzxlSVpMY47lVdQLdto9XJ/kvk2hLVV1UVfOrav6sWbMm01SSJEmaUi3heRMwZ+B4NrC5dYCq2tz/3ApcSrcNBOD2JEcC9D+3tvYpSZIkjUJLeF4DzE1ybJIDgdOBlS2dJ3lskoPHngMvBL7en14JLOqfLwI+PZmJS5IkSXvahLeqq6rtSZYAl9Hdqm5FVV2X5Jz+/PIkTwTWAo8DHkxyLt2dOZ4AXJpkbKwPV9Xn+q7fCnwsydnAzcBLp/TKJEmSpCnWdJ/nqloFrBpXtnzg+W102znGuxt45g76vAM4qXmmkiRJ0oj5DYOSJElSI8OzJEmS1MjwLEmSJDUyPEuSJEmNDM+SJElSI8OzJEmS1MjwLEmSJDUyPEuSJEmNDM+SJElSI8OzJEmS1MjwLEmSJDWaOeoJSNKetGz1+pGMe96C40YyriRparnyLEmSJDUyPEuSJEmNDM+SJElSI8OzJEmS1MjwLEmSJDUyPEuSJEmNDM+SJElSI+/zLEna74zqft6S9n+uPEuSJEmNDM+SJElSI8OzJEmS1MjwLEmSJDUyPEuSJEmNDM+SJElSI8OzJEmS1KgpPCc5Ocm6JBuSLB1y/mlJvpzkviSvHyifk+QLSa5Pcl2S1w6ce1OSW5Nc0z9OnZpLkiRJkqbHhF+SkmQGcAGwANgErEmysqq+MVDt28BrgNPGNd8OvK6qrk5yMHBVktUDbZdV1fm7exGSJEnSntCy8nwisKGqNlbV/cAlwMLBClW1tarWAA+MK99SVVf3z78LXA8cNSUzlyRJkvawlvB8FHDLwPEmdiEAJzkGeBbwlYHiJUmuTbIiyaGT7VOSJEnak1rCc4aU1WQGSXIQ8Ang3Kq6uy++EHgKcDywBXj7DtouTrI2ydpt27ZNZlhJkiRpSrWE503AnIHj2cDm1gGSHEAXnD9UVZ8cK6+q26vq+1X1IPAeuu0hD1NVF1XV/KqaP2vWrNZhJUmSpCnXEp7XAHOTHJvkQOB0YGVL50kCvBe4vqreMe7ckQOHLwG+3jZlSZIkaTQmvNtGVW1PsgS4DJgBrKiq65Kc059fnuSJwFrgccCDSc4F5gHPAM4Evpbkmr7LN1TVKuBtSY6n2wJyI/CqKbwuSZIkacpNGJ4B+rC7alzZ8oHnt9Ft5xjviwzfM01Vndk+TUmSJGn0/IZBSZIkqZHhWZIkSWpkeJYkSZIaGZ4lSZKkRoZnSZIkqZHhWZIkSWpkeJYkSZIaGZ4lSZKkRoZnSZIkqVHTNwxKjwTLVq8f9RQkSdJezpVnSZIkqZHhWZIkSWpkeJYkSZIaGZ4lSZKkRoZnSZIkqZHhWZIkSWpkeJYkSZIaGZ4lSZKkRoZnSZIkqZHhWZIkSWpkeJYkSZIaGZ4lSZKkRoZnSZIkqZHhWZIkSWpkeJYkSZIaGZ4lSZKkRoZnSZIkqZHhWZIkSWrUFJ6TnJxkXZINSZYOOf+0JF9Ocl+S17e0TXJYktVJbuh/Hrr7lyNJkiRNnwnDc5IZwAXAKcA84Iwk88ZV+zbwGuD8SbRdClxeVXOBy/tjSZIkaa/VsvJ8IrChqjZW1f3AJcDCwQpVtbWq1gAPTKLtQuDi/vnFwGm7dgmSJEnSntESno8Cbhk43tSXtdhZ2yOqagtA//PwYR0kWZxkbZK127ZtaxxWkiRJmnot4TlDyqqx/91p21Wuuqiq5lfV/FmzZk2mqSRJkjSlWsLzJmDOwPFsYHNj/ztre3uSIwH6n1sb+5QkSZJGoiU8rwHmJjk2yYHA6cDKxv531nYlsKh/vgj4dPu0JUmSpD1v5kQVqmp7kiXAZcAMYEVVXZfknP788iRPBNYCjwMeTHIuMK+q7h7Wtu/6rcDHkpwN3Ay8dIqvTZIkSZpSE4ZngKpaBawaV7Z84PltdFsymtr25XcAJ01mspIkSdIo+Q2DkiRJUiPDsyRJktTI8CxJkiQ1MjxLkiRJjQzPkiRJUiPDsyRJktTI8CxJkiQ1MjxLkiRJjQzPkiRJUiPDsyRJktTI8CxJkiQ1MjxLkiRJjQzPkiRJUiPDsyRJktTI8CxJkiQ1MjxLkiRJjQzPkiRJUqOZo56AJGn/tGz1+lFPQZKmnCvPkiRJUiPDsyRJktTI8CxJkiQ1MjxLkiRJjQzPkiRJUiPvtiFJkvY5o7yby3kLjhvZ2Bo9V54lSZKkRoZnSZIkqZHhWZIkSWrUFJ6TnJxkXZINSZYOOZ8k7+rPX5vkhL78qUmuGXjcneTc/tybktw6cO7UKb0ySZIkaYpN+IHBJDOAC4AFwCZgTZKVVfWNgWqnAHP7x7OBC4FnV9U64PiBfm4FLh1ot6yqzp+C65AkSZKmXcvK84nAhqraWFX3A5cAC8fVWQh8oDpXAockOXJcnZOAb1bVTbs9a0mSJGkEWsLzUcAtA8eb+rLJ1jkd+Mi4siX9No8VSQ5tmIskSZI0Mi3hOUPKajJ1khwIvBj4+MD5C4Gn0G3r2AK8fejgyeIka5Os3bZtW8N0JUmSpOnREp43AXMGjmcDmydZ5xTg6qq6faygqm6vqu9X1YPAe+i2hzxMVV1UVfOrav6sWbMapitJkiRNj5bwvAaYm+TYfgX5dGDluDorgZf3d914DvCdqtoycP4Mxm3ZGLcn+iXA1yc9e0mSJGkPmvBuG1W1PckS4DJgBrCiqq5Lck5/fjmwCjgV2ADcC5w11j7JY+ju1PGqcV2/LcnxdNs7bhxyXpIkSdqrTBieAapqFV1AHixbPvC8gFfvoO29wOOHlJ85qZlKkiRJI+Y3DEqSJEmNmlaepT1p2er1o56CJEnSUK48S5IkSY0Mz5IkSVIjw7MkSZLUyPAsSZIkNTI8S5IkSY0Mz5IkSVIjw7MkSZLUyPAsSZIkNTI8S5IkSY0Mz5IkSVIjw7MkSZLUaOaoJyBJjwTLVq8f2djnLThuZGNL0v7GlWdJkiSpkeFZkiRJamR4liRJkhoZniVJkqRGhmdJkiSpkeFZkiRJamR4liRJkhoZniVJkqRGhmdJkiSpkeFZkiRJauTXc0vSfm6UXw0uSfsbw7MkSdIkjOofpOctOG4k4+qh3LYhSZIkNTI8S5IkSY2awnOSk5OsS7IhydIh55PkXf35a5OcMHDuxiRfS3JNkrUD5YclWZ3khv7noVNzSZIkSdL0mDA8J5kBXACcAswDzkgyb1y1U4C5/WMxcOG48y+oquOrav5A2VLg8qqaC1zeH0uSJEl7rZaV5xOBDVW1saruBy4BFo6rsxD4QHWuBA5JcuQE/S4ELu6fXwyc1j5tSZIkac9rCc9HAbcMHG/qy1rrFPD5JFclWTxQ54iq2gLQ/zx8MhOXJEmS9rSWW9VlSFlNos5zq2pzksOB1Un+tar+qXWCfeBeDHD00Ue3NpMkSZKmXMvK8yZgzsDxbGBza52qGvu5FbiUbhsIwO1jWzv6n1uHDV5VF1XV/KqaP2vWrIbpSpIkSdOjJTyvAeYmOTbJgcDpwMpxdVYCL+/vuvEc4DtVtSXJY5McDJDkscALga8PtFnUP18EfHo3r0WSJEmaVhNu26iq7UmWAJcBM4AVVXVdknP688uBVcCpwAbgXuCsvvkRwKVJxsb6cFV9rj/3VuBjSc4GbgZeOmVXJUmSJE2Dpq/nrqpVdAF5sGz5wPMCXj2k3UbgmTvo8w7gpMlMVpIkSRolv2FQkiRJamR4liRJkhoZniVJkqRGhmdJkiSpkeFZkiRJamR4liRJkhoZniVJkqRGhmdJkiSpkeFZkiRJamR4liRJkhoZniVJkqRGhmdJkiSpkeFZkiRJamR4liRJkhoZniVJkqRGhmdJkiSpkeFZkiRJajRz1BPQ3mnZ6vWjnoIkSdJex5VnSZIkqZHhWZIkSWpkeJYkSZIaGZ4lSZKkRoZnSZIkqZHhWZIkSWpkeJYkSZIaGZ4lSZKkRoZnSZIkqZHhWZIkSWrUFJ6TnJxkXZINSZYOOZ8k7+rPX5vkhL58TpIvJLk+yXVJXjvQ5k1Jbk1yTf84deouS5IkSZp6MyeqkGQGcAGwANgErEmysqq+MVDtFGBu/3g2cGH/czvwuqq6OsnBwFVJVg+0XVZV50/d5UiSJEnTp2Xl+URgQ1VtrKr7gUuAhePqLAQ+UJ0rgUOSHFlVW6rqaoCq+i5wPXDUFM5fkiRJ2mNawvNRwC0Dx5t4eACesE6SY4BnAV8ZKF7Sb/NYkeTQYYMnWZxkbZK127Zta5iuJEmSND1awnOGlNVk6iQ5CPgEcG5V3d0XXwg8BTge2AK8fdjgVXVRVc2vqvmzZs1qmK4kSZI0PVrC8yZgzsDxbGBza50kB9AF5w9V1SfHKlTV7VX1/ap6EHgP3fYQSZIkaa/VEp7XAHOTHJvkQOB0YOW4OiuBl/d33XgO8J2q2pIkwHuB66vqHYMNkhw5cPgS4Ou7fBWSJEnSHjDh3TaqanuSJcBlwAxgRVVdl+Sc/vxyYBVwKrABuBc4q2/+XOBM4GtJrunL3lBVq4C3JTmebnvHjcCrpuiaJEmSpGkxYXgG6MPuqnFlyweeF/DqIe2+yPD90FTVmZOaqSRJkjRifsOgJEmS1MjwLEmSJDUyPEuSJEmNDM+SJElSI8OzJEmS1MjwLEmSJDUyPEuSJEmNDM+SJElSI8OzJEmS1MjwLEmSJDUyPEuSJEmNDM+SJElSI8OzJEmS1MjwLEmSJDUyPEuSJEmNZo56ApIkSZrYstXrRzb2eQuOG9nYextXniVJkqRGhmdJkiSpkds29nKjfItGkiQJRpdH9sbtIq48S5IkSY0Mz5IkSVIjw7MkSZLUyD3PDdx3LEmSJHDlWZIkSWpmeJYkSZIaGZ4lSZKkRoZnSZIkqVFTeE5ycpJ1STYkWTrkfJK8qz9/bZITJmqb5LAkq5Pc0P88dGouSZIkSZoeE4bnJDOAC4BTgHnAGUnmjat2CjC3fywGLmxouxS4vKrmApf3x5IkSdJeq2Xl+URgQ1VtrKr7gUuAhePqLAQ+UJ0rgUOSHDlB24XAxf3zi4HTdu9SJEmSpOnVEp6PAm4ZON7Ul7XU2VnbI6pqC0D/8/D2aUuSJEl7XsuXpGRIWTXWaWm788GTxXRbQQDuSbJuMu01Ek8AvjXqSWi433nhU3e1qa/r/snXdf+0y6/rbvw3QtPvEff39XdGO/yThxW2hOdNwJyB49nA5sY6B+6k7e1JjqyqLf0Wj63DBq+qi4CLGuapvUSStVU1f9Tz0NTydd0/+brun3xd90++rnuHlm0ba4C5SY5NciBwOrByXJ2VwMv7u248B/hOvxVjZ21XAov654uAT+/mtUiSJEnTasKV56ranmQJcBkwA1hRVdclOac/vxxYBZwKbADuBc7aWdu+67cCH0tyNnAz8NIpvTJJkiRpiqVqUluQpQklWdxvt9F+xNd1/+Trun/ydd0/+bruHQzPkiRJUiO/nluSJElqZHjWtEry+iSV5Amjnot2X5K/SvKvSa5NcmmSQ0Y9J+2aJCcnWZdkQxK/4XU/kGROki8kuT7JdUleO+o5aeokmZHkX5J8ZtRzeaQzPGvaJJkDLKD7QKj2D6uBn6iqZwDrgd8f8Xy0C5LMAC4ATgHmAWckmTfaWWkKbAdeV1VPB54DvNrXdb/yWuD6UU9ChmdNr2XA7zLJL8bR3quqPl9V2/vDK+nu3a59z4nAhqraWFX3A5cAC0c8J+2mqtpSVVf3z79LF7TGfyOw9kFJZgO/BPztqOciw7OmSZIXA7dW1VdHPRdNm98CPjvqSWiXHAXcMnC8CUPWfiXJMcCzgK+MeCqaGu+kW4x6cMTzEG3fMCgNleT/AE8ccuqNwBuAF+7ZGWkq7Ox1rapP93XeSPcW8Yf25Nw0ZTKkzHeI9hNJDgI+AZxbVXePej7aPUleBGytqquSPH/E0xGGZ+2GqvqFYeVJfhI4FvhqEuje2r86yYlVddsenKJ2wY5e1zFJFgEvAk4q73W5r9oEzBk4ng1sHtFcNIWSHEAXnD9UVZ8c9Xw0JZ4LvDjJqcCjgccl+buq+s0Rz+sRy/s8a9oluRGYX1XfGvVctHuSnAy8A/i5qto26vlo1ySZSfeBz5OAW4E1wG8MfAOs9kHpVisuBr5dVeeOeDqaBv3K8+ur6kUjnsojmnueJU3Gu4GDgdVJrkmyfNQT0uT1H/pcAlxG96Gyjxmc9wvPBc4Efr7/+3lNv1opaQq58ixJkiQ1cuVZkiRJamR4liRJkhoZniVJkqRGhmdJkiSpkeFZkiRJamR4liRJkhoZniVJkqRGhmdJkiSp0f8HYNlgystTa20AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Generate samples from a mixture of two normal distributions\n",
    "N = 1000\n",
    "x = np.concatenate((stats.norm.rvs(loc=-2, scale=1, size=1000), stats.norm.rvs(loc=2, scale=1, size=1000)))\n",
    "mu = np.mean(x)  # sample mean\n",
    "sigma = np.std(x, ddof=1)  # sample standard deviation\n",
    "\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.hist(x, bins=20, alpha=0.5, density=True)\n",
    "plt.axvline(mu - sigma, color='black')\n",
    "plt.axvline(mu + sigma, color='black')\n",
    "plt.title('Proportion of data in range [μ–σ, μ+σ]: %.3f' % test_statistic(x))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The black lines show the mean plus and minus one standard deviation.\n",
    "\n",
    "In this case, 60.2% of the data lie between the black lines. \n",
    "\n",
    "If we apply this test statistic to the original dataset, we get the following."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsUAAAF1CAYAAAAA6ZfwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAeVklEQVR4nO3dfbSlV10f8O/PmUSroBEyiOaFRJqIaMWXMYgvlbZGA2KDLatGWCgUTbFGJWo1VQHfLVqNWoJpVg1Uag1akE5xEGepEfGtExCpATOOEcgYkCEmYCASBn/94zwXDzf3zpyZnHPvTPbns9ZZ65zn2c+z97PvvnO/s88+56nuDgAAjOyjtrsBAACw3YRiAACGJxQDADA8oRgAgOEJxQAADE8oBgBgeEIxcL9U1dOq6je2od4vqqo/r6q7q+rJC5R/SVX98BY0LVV1bVU9dyvq2kpV9daquqeqXrrdbVmlqvroaVx9cKvGDLD9hGI4Cc2Fj7ur6q+r6sVV9aCToF3nVVVX1c61bd39i9395dvQnB9M8sLuflB3v3KZJ66qG6vqG070+O5+dnf/0DLbdBL5qu5++nY3YpW6+wPd/aAkv3g8x1XVlVX1zqp6T1VdX1UffZSyO6rqh6vq9qr626r646o6Y4Nyv7X+d66qHlJVv1pV76uqt1XVU4+nncDGhGI4eX3V9If5c5N8fpLvW19g/g/lqm1lXQt6RJKbt7sRq1Qz/p0+QVX1+Kq6cYvq+ookVyX5F0nOS/KpSX7gKIf8QJIvTPK4JB+f5OlJ/m7dOZ+WZKPfu2uS3Jvkk5I8LcnPVdVn3L8rAPxjCye57v6rJK9O8plJMs0afXNV/XmSP5+2fWNVHayqv6mqPVX1KWvHT+W/tapurap3V9VPrAWtqvqoqvq+abbpXVX1C1X1CdO+tVnhZ1XV25P8VpLXTqe9a5rFflxVPaOqXjdX3xdW1f5ptmx/VX3h3L4bq+qHqur3ptmx36iqMze79s2uq6r+IrPQ8X+mdtxnRq6qPqeq3jDV87IkHzO37xOr6lVVdbiq7pyenz3t+5EkX5LkhdO5Xzht/5mquq2q3ltVr6+qLzlKuz+8VGMKZoeq6jumPn5HVT3zKMfeWFU/UlW/l+T9ST61qp5ZVW+ZruXWqvp3c+WPev6qemhV/Z+p3fun2cn5n9ejqmrf1Me3VNW/2axtG7T1Pu8c1AKz7DVbnnDN1NYPTefoqnrlAnV+xHibtnVV/eNF2z0dc05V7Z2uu+cezzme88z5+iQ/3903d/edSX4oyTM2qfsTkzwnyTd299t65k+7++/mynxCkucn+a51x35ckn+d5LndfXd3vy7JnsxCNXA/CMVwkquqc5I8Mckfz21+cpLHJnl0Vf3zJD+W5N8k+eQkb0tyw7rTfHWS3ZnNOl+a5N9O258xPf5ZZiHzQUleuO7YL03y6Um+Isk/nbadMS1b+IN1bX1Ikl9L8rNJHprkp5L8WlU9dK7YU5M8M8nDkpye5Ds3ue5Nr6u7H5nk7Zlm07v7A+uOPT3JK5O8NMlDkvxKZkFizUcleXFms83nJrln7bq7+3uT/G6SK6ZzXzEdsz/JZ0/n+59JfqWqPiaLeXiST0hyVpJnJblmCkabeXqSy5M8eLrudyV5UmYzis9McnVVfe6C578myfumMl8/PZJ8OGDtm67nYUm+NsmLavWzjlck+bL8wyzpryS5McmVK6533osy+7mfOz3+LMl/ziZLJqrq3Kq6q6rO3eR8n5HkT+Ze/0mST1o39tf8kyRHkjylZsstDlTVN68r86NJfi7JO9dtvzDJh7r7wLq6zBTD/SQUw8nrlVV1V5LXJfmdzP5Irvmx7v6b7r4ns7dPr+/uN0zh8D8meVxVnTdX/gVT+bcn+enMwk+mY3+qu2/t7runYy+rj1wq8f3d/b6prmP5yiR/3t0v7e4j3f1LmYWNr5or8+LuPjCd75czC5obWeS6NvMFSU5L8tPd/cHu/l+ZhdokSXff0d0v7+73d/ffJvmRzML/prr7f0zHHenun0zy0Uk+bYG2JMkHk/zg1Ja9Se4+xrEvmWYcj0zH/Fp3/8U0o/g7SX4js9nso56/qnZk9p+B50/X+uYk/33uuCcleWt3v3iq6w1JXp7kKQte14n6l5n9bN7a3e/LbGnQFyf5qxXXm+TDS4GekOT7ptnW25JcneSi7j680THd/fbuPmP6HdrIg5K8Z+712vMHb1D27Mz+E3NhkvMz6+/vr6qLp/btTvJFSf7LAvWs1bVRPcBxEIrh5PXk6Y/wI7r7368LpbfNPf+UzGYTkyRTuL0js1nDjcq/bTrmPsdOz3dmtlZxo2OPZf351s4535b5ma/3Z/ZH/pjn2uS6jtaOv+ruXteOJElVfWxV/deaLRt5b2bLQs6YQuSGpuUJb6nZspC7Mgs1my79WOeO7j4y9/po152s6/OqekJV/eH0Vv9dmb1zMF/3ZuffldnPc/58888fkeSx0wzoXdO5n5bZrPL9VlV/MS1BWXu8ddr1sMxm+tesjbtdS6jzqrlreVWSL153fcnsXYwdG7ThU3Li7s5s1nvN2vO/3aDs2u/yD3b3Pd39pszeBXlizZY2vSjJt637mW5Wz1pdG9UDHAehGE5N82Hv9szCTZIPvyX+0HzkrNs5c8/PnY65z7HTviNJ/nqTuuafb2T9+dbOeSIzgItc12bekeSsqqp17VjzHZnN1D62uz8+/7AsZK38R1zntH74uzNbyvGJ3X1GZrNz8+dfpg/XX7P10i/P7K39T5rq3rtg3Ycz+3mePbdtfizcluR3pv98rT0e1N3fdJztnV9GcsaHL6L7kdP51h7nTbsOZTZDuub83HfcLVTftPb2w7r7P61dS2Yz4a+bv76p2OEkH9igDYcWrH8jNyd5zNzrxyT56+6+Y4Oyb1pr7gb7Pj6zpU4vq6p35h/e4Tg0jcMDSXZW1QXr6npAf+gUtoJQDKe+/5nkmVX12VOA+tEkf9Tdb50r8x9q9uGyc5J8W5KXTdt/KcmVVXV+zb7y7UeTvGyTGapkFib+PrP1xxvZm+TCqnpqVe2sqq9J8ujMZuxWcV2b+YPMQta3Tu34V0kumtv/4Mxm6+6a1kE/f93xf52PvMYHT+c7nFkgeV7uO1u3KqdntlTjcJIjVfWEJAt9BV53fyjJKzJ7a/5jq+pRSb5ursirMvt5Pb2qTpsen19Vn36cbXxWzT60eXFm62UfXFWnHaX8S5M8p6o+daNxV7MPKr7kKMd/VlVdNM3sr61Dfsiije3uv89sfP1IVT14WpLz7TnOr2Bb5xcy64dHT+u5vy/JSzap/y8yW7f+vTX70OGnJ/mazH4e78lsxvqzp8cTp8M+L7Px/77MfqY/WFUfV1VflNnnBB7Q3x0NW0EohlNcd/9mkudmNpv4jiSPTHLZumL/O8nrk7wxsw/C/fy0/frM/pi+NslfZvaVUN9ylLren9n629+b3o7+gnX778hsdu47Mlvq8F1JntTd717RdW127L1J/lVmHyK8M7PA8Yq5Ij+d5B8leXeSP0zy6+tO8TOZfQjqzqr62SSvyewbQA5k9jb73+X4lpWcsGnN87dmtv76zsw+qLjnOE5xRWZLPd6Z2c/6lzKbJV0795dn1q+3T2VekFkIPx6PzCy0f3dmXzX2vMzWy27mpZmNvd/ObNy9f2rnmnOS/N5Rjn9zkh/O7AOID83sg3ovO0r5jTwnyV1Jbs1s/P+PqU0bmj5od/dmH7Tr7l9P8uOZXdPbpsfz545/dVV9z9whX5vZOyF3ZPY7+dzu/s1p3fg71x6Z9Wsym3W+d3r+7zMbv+/K7Of5Td1tphjup/rIJXfAA01VdZILuvvgdreF7VdVL0jy8O7++mMWvu+xt2T2TSC/2t1fP82w/mWS047y7sLx1nF6Zt+m8Fnd/cEN9j8jyTd09xcvo75N2vDRmb1bcFqSH+/uo33fMPAAcbJ9GT8ASzQtmTg9yf/L7CYwz0pyQnfr6+5Fv23jhE2zoce7fGPZbfhA5tZGA2MQigEe2B6c2Vvsn5LZ2+0/mdlyGgDmWD4BAMDwfNAOAIDhCcUAAAxv29YUn3nmmX3eeedtV/XA/XTLLbckST7t01b+2SvgJOH3nlPV61//+nd391Hvmrltofi8887LTTfdtF3VA/fT4x//+CTJjTfeuK3tALaO33tOVVX1tmOVsXwCAIDhCcUAAAxPKAYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMDyhGACA4QnFAAAMTygGAGB4QjEAAMMTigEAGN7O7W4AwIm6et+Bbav7yosv3La6AVg+M8UAAAxPKAYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMDyhGACA4QnFAAAMTygGAGB4QjEAAMMTigEAGJ5QDADA8IRiAACGJxQDADA8oRgAgOEJxQAADE8oBgBgeEIxAADDE4oBABieUAwAwPCEYgAAhicUAwAwPKEYAIDhCcUAAAxPKAYAYHgLheKquqSqbqmqg1V11SZlHl9Vb6yqm6vqd5bbTAAAWJ2dxypQVTuSXJPk4iSHkuyvqj3d/ea5MmckeVGSS7r77VX1sBW1FwBYkav3HTjq/kN33rNQuRNx5cUXLv2ccDwWmSm+KMnB7r61u+9NckOSS9eVeWqSV3T325Oku9+13GYCAMDqLBKKz0py29zrQ9O2eRcm+cSqurGqXl9VX7fRiarq8qq6qapuOnz48Im1GAAAlmyRUFwbbOt1r3cm+bwkX5nkK5I8t6ru8z5Id1/X3bu7e/euXbuOu7EAALAKx1xTnNnM8Dlzr89OcvsGZd7d3e9L8r6qem2SxyRZ/qIjAABYskVmivcnuaCqzq+q05NclmTPujL/O8mXVNXOqvrYJI9N8pblNhUAAFbjmDPF3X2kqq5I8pokO5Jc3903V9Wzp/3XdvdbqurXk7wpyd8n+W/d/aerbDgAACzLIssn0t17k+xdt+3ada9/IslPLK9pAACwNdzRDgCA4QnFAAAMTygGAGB4QjEAAMMTigEAGJ5QDADA8IRiAACGJxQDADA8oRgAgOEJxQAADE8oBgBgeEIxAADDE4oBABieUAwAwPCEYgAAhicUAwAwPKEYAIDhCcUAAAxPKAYAYHhCMQAAwxOKAQAY3s7tbgAAnIyu3ndg2+q+8uILt61uGJWZYgAAhicUAwAwPKEYAIDhCcUAAAxPKAYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMDyhGACA4QnFAAAMTygGAGB4QjEAAMMTigEAGN5CobiqLqmqW6rqYFVdtcH+x1fVe6rqjdPjectvKgAArMbOYxWoqh1JrklycZJDSfZX1Z7ufvO6or/b3U9aQRsBYChX7zuw3U2A4SwyU3xRkoPdfWt335vkhiSXrrZZAACwdRYJxWcluW3u9aFp23qPq6o/qapXV9VnLKV1AACwBY65fCJJbbCt171+Q5JHdPfdVfXEJK9McsF9TlR1eZLLk+Tcc889vpYCAMCKLDJTfCjJOXOvz05y+3yB7n5vd989Pd+b5LSqOnP9ibr7uu7e3d27d+3adT+aDQAAy7NIKN6f5IKqOr+qTk9yWZI98wWq6uFVVdPzi6bz3rHsxgIAwCocc/lEdx+pqiuSvCbJjiTXd/fNVfXsaf+1SZ6S5Juq6kiSe5Jc1t3rl1gAAMBJaZE1xWtLIvau23bt3PMXJnnhcpsGAABbwx3tAAAYnlAMAMDwhGIAAIa30JpigKNxS1oATnVmigEAGJ5QDADA8IRiAACGJxQDADA8oRgAgOEJxQAADE8oBgBgeEIxAADDE4oBABieUAwAwPCEYgAAhicUAwAwPKEYAIDhCcUAAAxPKAYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMDyhGACA4QnFAAAMTygGAGB4QjEAAMMTigEAGJ5QDADA8IRiAACGt3O7GwAAR3P1vgPb3QRgAGaKAQAYnlAMAMDwhGIAAIYnFAMAMDyhGACA4QnFAAAMb6FQXFWXVNUtVXWwqq46SrnPr6oPVdVTltdEAABYrWOG4qrakeSaJE9I8ugkX1tVj96k3AuSvGbZjQQAgFVaZKb4oiQHu/vW7r43yQ1JLt2g3LckeXmSdy2xfQAAsHKLhOKzktw29/rQtO3DquqsJF+d5NrlNQ0AALbGIqG4NtjW617/dJLv7u4PHfVEVZdX1U1VddPhw4cXbCIAAKzWzgXKHEpyztzrs5Pcvq7M7iQ3VFWSnJnkiVV1pLtfOV+ou69Lcl2S7N69e32wBgCAbbFIKN6f5IKqOj/JXyW5LMlT5wt09/lrz6vqJUletT4QAwDAyeqYobi7j1TVFZl9q8SOJNd3981V9expv3XEAACc0haZKU53702yd922DcNwdz/j/jcLAAC2jjvaAQAwvIVmigH4SFfvO7At9V558YXbUu92XS/AVjFTDADA8IRiAACGJxQDADA8oRgAgOEJxQAADE8oBgBgeEIxAADDE4oBABieUAwAwPCEYgAAhicUAwAwPKEYAIDhCcUAAAxPKAYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMDyhGACA4QnFAAAMTygGAGB4O7e7AcByXL3vwJbWd+jOe7al3tHpb4DVMFMMAMDwhGIAAIYnFAMAMDyhGACA4QnFAAAMTygGAGB4QjEAAMMTigEAGJ5QDADA8IRiAACGJxQDADA8oRgAgOEJxQAADG+hUFxVl1TVLVV1sKqu2mD/pVX1pqp6Y1XdVFVfvPymAgDAauw8VoGq2pHkmiQXJzmUZH9V7enuN88V+80ke7q7q+qzkvxykketosEAALBsi8wUX5TkYHff2t33JrkhyaXzBbr77u7u6eXHJekAAMAp4pgzxUnOSnLb3OtDSR67vlBVfXWSH0vysCRfudGJquryJJcnybnnnnu8bQUAWKqr9x3YlnqvvPjCbamXzS0yU1wbbLvPTHB3/2p3PyrJk5P80EYn6u7runt3d+/etWvXcTUUAABWZZFQfCjJOXOvz05y+2aFu/u1SR5ZVWfez7YBAMCWWCQU709yQVWdX1WnJ7ksyZ75AlX1j6uqpuefm+T0JHcsu7EAALAKx1xT3N1HquqKJK9JsiPJ9d19c1U9e9p/bZJ/neTrquqDSe5J8jVzH7wDAICT2iIftEt3702yd922a+eevyDJC5bbNAAA2BruaAcAwPCEYgAAhicUAwAwPKEYAIDhCcUAAAxPKAYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMDyhGACA4QnFAAAMTygGAGB4QjEAAMMTigEAGJ5QDADA8IRiAACGJxQDADA8oRgAgOEJxQAADE8oBgBgeEIxAADDE4oBABieUAwAwPCEYgAAhicUAwAwvJ3b3QBYhav3Hdi2uq+8+MJtqxsAODFmigEAGJ5QDADA8IRiAACGJxQDADA8oRgAgOEJxQAADE8oBgBgeEIxAADDE4oBABieUAwAwPAWCsVVdUlV3VJVB6vqqg32P62q3jQ9fr+qHrP8pgIAwGrsPFaBqtqR5JokFyc5lGR/Ve3p7jfPFfvLJF/a3XdW1ROSXJfksatoMJzsrt53YLubAHDK8W8n222RmeKLkhzs7lu7+94kNyS5dL5Ad/9+d985vfzDJGcvt5kAALA6i4Tis5LcNvf60LRtM89K8uqNdlTV5VV1U1XddPjw4cVbCQAAK7RIKK4NtvWGBav+WWah+Ls32t/d13X37u7evWvXrsVbCQAAK3TMNcWZzQyfM/f67CS3ry9UVZ+V5L8leUJ337Gc5gEAwOotMlO8P8kFVXV+VZ2e5LIke+YLVNW5SV6R5OndbaU8AACnlGPOFHf3kaq6IslrkuxIcn1331xVz572X5vkeUkemuRFVZUkR7p79+qaDQAAy7PI8ol0994ke9dtu3bu+Tck+YblNg0AALaGO9oBADA8oRgAgOEJxQAADE8oBgBgeEIxAADDE4oBABieUAwAwPCEYgAAhicUAwAwPKEYAIDhCcUAAAxPKAYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMDyhGACA4QnFAAAMTygGAGB4QjEAAMMTigEAGJ5QDADA8IRiAACGJxQDADA8oRgAgOHt3O4GAACM5up9B7at7isvvnDb6j6ZmSkGAGB4QjEAAMMTigEAGJ5QDADA8IRiAACGJxQDADA8oRgAgOEJxQAADE8oBgBgeEIxAADDE4oBABjeQqG4qi6pqluq6mBVXbXB/kdV1R9U1Qeq6juX30wAAFidnccqUFU7klyT5OIkh5Lsr6o93f3muWJ/k+Rbkzx5FY0EAIBVWmSm+KIkB7v71u6+N8kNSS6dL9Dd7+ru/Uk+uII2AgDASi0Sis9Kctvc60PTtuNWVZdX1U1VddPhw4dP5BQAALB0i4Ti2mBbn0hl3X1dd+/u7t27du06kVMAAMDSLRKKDyU5Z+712UluX01zAABg6y0SivcnuaCqzq+q05NclmTPapsFAABb55jfPtHdR6rqiiSvSbIjyfXdfXNVPXvaf21VPTzJTUk+PsnfV9Vzkjy6u9+7uqYDAMByHDMUJ0l3702yd922a+eevzOzZRUAAHDKcUc7AACGt9BMMae+q/cd2O4mAACctMwUAwAwPKEYAIDhCcUAAAxPKAYAYHhCMQAAwxOKAQAYnlAMAMDwhGIAAIYnFAMAMDyhGACA4QnFAAAMTygGAGB4QjEAAMMTigEAGJ5QDADA8IRiAACGt3O7GwAAwNa5et+Bban3yosv3JZ6F2WmGACA4QnFAAAMTygGAGB4QjEAAMMTigEAGN6Q3z6xXZ+6BADg5GSmGACA4QnFAAAMTygGAGB4QjEAAMMTigEAGJ5QDADA8IRiAACGJxQDADA8oRgAgOEJxQAADE8oBgBgeEIxAADDWygUV9UlVXVLVR2sqqs22F9V9bPT/jdV1ecuv6kAALAaxwzFVbUjyTVJnpDk0Um+tqoeva7YE5JcMD0uT/JzS24nAACszCIzxRclOdjdt3b3vUluSHLpujKXJvmFnvnDJGdU1Scvua0AALASi4Tis5LcNvf60LTteMsAAMBJaecCZWqDbX0CZVJVl2e2vCJJ7q6qWxao/3icmeTdSz4nm9PfW+uk7O9v//JP2+4mrMpJ2d8PYPp769zvvn4A/96vgrE9+fatqWaz/n7EsQ5cJBQfSnLO3Ouzk9x+AmXS3dcluW6BOk9IVd3U3btXdX4+kv7eWvp7a+nvraW/t46+3lr6e2vdn/5eZPnE/iQXVNX5VXV6ksuS7FlXZk+Sr5u+heILkrynu99xIg0CAICtdsyZ4u4+UlVXJHlNkh1Jru/um6vq2dP+a5PsTfLEJAeTvD/JM1fXZAAAWK5Flk+ku/dmFnznt10797yTfPNym3ZCVrY0gw3p762lv7eW/t5a+nvr6Outpb+31gn3d83yLAAAjMttngEAGN4pGYrddnprLdDfj6+q91TVG6fH87ajnQ8EVXV9Vb2rqv50k/3G9hIt0N/G9pJU1TlV9dtV9Zaqurmqvm2DMsb3kizY38b3klTVx1TV/62qP5n6+wc2KGN8L8mC/X3c43uhNcUnk7nbTl+c2VfB7a+qPd395rli87edfmxmt51+7Fa39YFgwf5Okt/t7idteQMfeF6S5IVJfmGT/cb2cr0kR+/vxNheliNJvqO731BVD07y+qra59/ulVmkvxPje1k+kOSfd/fdVXVaktdV1aunu/yuMb6XZ5H+To5zfJ+KM8VuO721FulvlqS7X5vkb45SxNheogX6myXp7nd09xum53+b5C25751Pje8lWbC/WZJpzN49vTxteqz/0JbxvSQL9vdxOxVDsdtOb61F+/Jx09sYr66qz9iapg3J2N56xvaSVdV5ST4nyR+t22V8r8BR+jsxvpemqnZU1RuTvCvJvu42vldogf5OjnN8n4qheGm3nWYhi/TlG5I8orsfk+S/JHnlqhs1MGN7axnbS1ZVD0ry8iTP6e73rt+9wSHG9/1wjP42vpeouz/U3Z+d2V19L6qqz1xXxPheogX6+7jH96kYipd222kWcsy+7O73rr2NMX2n9WlVdebWNXEoxvYWMraXa1r79/Ikv9jdr9igiPG9RMfqb+N7Nbr7riQ3Jrlk3S7jewU26+8TGd+nYih22+mtdcz+rqqHV1VNzy/KbFzdseUtHYOxvYWM7eWZ+vHnk7ylu39qk2LG95Is0t/G9/JU1a6qOmN6/o+SfFmSP1tXzPhekkX6+0TG9yn37RNuO721FuzvpyT5pqo6kuSeJJe1u8KckKr6pSSPT3JmVR1K8vzMPkBgbK/AAv1tbC/PFyV5epL/N60DTJLvSXJuYnyvwCL9bXwvzycn+e/TNzZ9VJJf7u5XySYrs0h/H/f4dkc7AACGdyounwAAgKUSigEAGJ5QDADA8IRiAACGJxQDADA8oRgAgOEJxQAADE8oBgBgeP8fPFAyVV/zmFoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mu = np.mean(data)\n",
    "sigma = np.std(data)\n",
    "\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.hist(data, bins=20, alpha=0.5, density=True)\n",
    "plt.axvline(mu - sigma, color='black')\n",
    "plt.axvline(mu + sigma, color='black')\n",
    "plt.title('Proportion of data in range [μ–σ, μ+σ]: %.3f' % test_statistic(data))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the test statistic evaluated on the real dataset is 0.640. But is this statistically significantly different from the expected value, which is 0.680 if we assume normally distributed data in our model?\n",
    "\n",
    "We can only answer this question by comparing the test statistic for the real dataset (0.640) to the test statistic on replicated data from the posterior predictive distribution of our model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task: Calculate the p-value\n",
    "\n",
    "* Generate samples from the posterior predictive distribution and compute the test statistic for each replicated dataset.\n",
    "* Plot a histogram of the test statistic on the replicated datasets.\n",
    "* Plot a vertical line on the histogram to show value of the test statistic on the real dataset (0.640).\n",
    "* Compute the p-value as the proportion of replicated test statistic values that are greater than the real test statistic value.\n",
    "\n",
    "### How to generate samples\n",
    "\n",
    "* Generate 1000 (or more) samples of $(x,\\sigma^2)$ from the posterior Normal-Inverse-Gamma distribution.\n",
    "* For each sample from the posterior, generate a replicated dataset $\\{y^{\\text{(rep)}}_i\\}$ with the same size (the same number of data points) as the real dataset. It is important the each replicated dataset has the same size as the real dataset, to make sure their statistical behavior is the same.\n",
    "* Compute the test statistic for each replicated dataset to get the samples from the replicated test statistic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def norminvgamma_rvs(mu, nu, alpha, beta, size=1):\n",
    "    '''\n",
    "    Generate n samples from the normal-inverse-gamma distribution. This function\n",
    "    returns a (size x 2) matrix where each row contains a sample, (x, sigma2).\n",
    "    '''\n",
    "    sigma2 = stats.invgamma.rvs(a=alpha, scale=beta, size=size)  # Sample sigma^2 from the inverse-gamma\n",
    "    x = stats.norm.rvs(loc=mu, scale=np.sqrt(sigma2 / nu), size=size)  # Sample x from the normal\n",
    "    return np.vstack((x, sigma2)).transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9986\n"
     ]
    }
   ],
   "source": [
    "posterior_samples = norminvgamma_rvs(mu_1, nu_1, alpha_1, beta_1, size = 10000)\n",
    "\n",
    "test_values = []\n",
    "\n",
    "# Generate a replicated dataset with each parameter value\n",
    "for posterior_sample in posterior_samples:\n",
    "    replicated_dataset = stats.norm.rvs(loc = posterior_sample[0], scale = np.sqrt(posterior_sample[1]), size = s0)\n",
    "    test_values.append(test_statistic(replicated_dataset))\n",
    "    \n",
    "p_value = sum(test_values > test_statistic(data))/len(test_values)\n",
    "print(p_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEICAYAAABGaK+TAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVbklEQVR4nO3de7hldX3f8fcnXKLcBGRQUMpEZCaAT7g44gW1NEgCJgZ4KibEIKh0NImJiIml2DzStEltvJCLjRQRQx8VQ5FES4w6JYBQrTpcVHDkoqIMTpmD3JFGgW//WOvI5nDOnH32Ppf5Me/X86znrPv6/vae/dlr/9bae1JVSJLa8zNLXYAkaTQGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwAZDk5CRXLcFxVya5Nsn9SX5/sY+/WJK8LMmNI277j0lOmu+a1D4DvAFJbk3yUJIHktyR5CNJdljquubJO4DLq2rHqvrLqQuTXJ7klHEPkuTwJOsXa9skleS5k9NVdWVVrRxiuzOTfHRwXlUdXVXnz+X4M+x7eZLLkvwoybeSvGIT6+6c5PwkG/vhzCnLX5LkK/0b79eTvHRgWZK8M8n3k9yX5BNJdhq3fj2RAd6OV1XVDsAhwAuAf7/E9cyXvYEblrqILcQFwLXA04F3AhclWTbDumcB2wHLgUOBE5O8HiDJrsCngfcAOwN/BvzPJLv0274OOBE4DNgTeCrwV/PfHFFVDpv5ANwKvGJg+j3AJdOsdzbw3inzPgWc1o+fDnwbuB/4JnDcwHonA1f148uBArYeWH45cMrA9BuAdcDdwOeAvTdR/6/RhfQ9/X726+f/E/AI8P+AB4AVU7b7kynLP9DP/3lgDXAXcCPwmoFtXtm37X7gduAPgO2Bh4BH+/08AOw5TZ1Db0sXal/q27QB+ACwbb+fL/SP34P9+r8OHA6sHzjWv+2PcX/fhiOAo4AfAz/pt/vaDI/9v+kf+8nn8ZAh/g2tAP4Z2HFg3pXAm2dY/07gBQPTZwBX9uO/CtwwZf2bgDf24xcBfziw7CX9c7jdUr+WnmzDkhfgMMSTNBDgwF59GP7HadZ7OXAbkH56lz589uynj+/D52f6UHkQ2KNfdjJDBjhwLHALsB+wNd2ngS/OUPuK/jhHAtvQdZncMhB2jwunabafGl7b9218fX/sQ/qwOaBfvgF42UD7D+nHHxegMxxr6G2B5wMv6mtY3gfqqQPLC3juwPRP9wGs7Nsw+bwsB/bpx88EPjrTY9A/h7fTfQoL8Fz6N0/gr4G/nqFtxwHrpsz7APBXM6x/J3DowPQ7gbv78VcB35yy/s3AWf34J4F3DCw7rH88Dlzq19KTbbALpR1/n+Qe4CrgCuBPp1nnSroXysv66VcDX6qqHwBU1f+oqh9U1aNV9bd0L7pDR6jlTcB/rqp1VfVwX8tBSfaeZt1fB/6hqtZU1U+A99J9pH7JCMeF7uzv1qr6SFU9XFXX0AXGq/vlPwH2T7JTVd3dLx/W0NtW1dVV9X/6Gm4F/hvwL4c8ziPAz/bH2qaqbq2qbw+57SnAn1XVV6tzS1V9r6/pd6rqd2bYbgfg3inz7gV2nGH9zwKnJ9mx78t/A12XCsAXgT2TnJBkm/4C6z4Dy/8ROKXvc38a3acNBpZrnhjg7Ti2qnauqr37F+pDSc7oL2w+kOTs6k53PgGc0G/zm8DHJneQ5HVJrktyT/9m8DxgtxFq2Rv4i4H93EV3NvisadbdE/je5ERVPUp39jndusMe+4WTx+6P/1rgmf3yf03XFfK9JFckefEc9j30tklWJLkkyf9Nch/dm9hQj2VV3QKcSne2vbG/yLfnkDXuRdcNNlcPAFMvJO5E1w0znd+n+/R2M1033AXAeoCq+iFwDHAacAdd18//mlwOnNevfzndp8XL+vkjXUTWzAzwhlXVn1bVDv3w5n72BcCr+7PhF9KdndJPfwh4C/D0qtoZuJ4ueKd6sP87eMb0zIHx24A39W8ok8NTq+qL0+zrB3ShS19H6ELo9mGbOWX6NuCKKcfeoap+G6A/Mz0G2B34e+DCGfbzxAPNbdsPAt8C9q2qnej6iKd7LGc61ser6qV0j00B/2XIOm+jO9udqxuA5yQZPOM+kBkuIFfVXVX12qp6ZlUdQJcVXxlYfkVVvaCqdqW7YLlycnn/Ce9dVbW8qp7dH+N2hn/ONSQD/Emmqq4FJoBzgc9V1T39ou3pwmECoL+j4Hkz7GOC7sX2W0m2SvIGHh8aZwP/LskB/b6eluT4GUq6EPiVJEck2QZ4O93FtOnCfjp3AM8ZmL4EWJHkxP7j+zZJXpBkvyTbJnltkqf13TX30XVXTO7n6f1H+icYYdsd+3UeSPLzwG/PUvfgsVYm+cUkP0t3ce+hKcdanmSm1+a5wB8keX5/u95zZ+i6epyqugm4DnhXkqckOQ74Bfo3+Glq3CfJ0/vn/2hgNfCfBpYf3D/2O9F1i62vqs/1y3btt0+S/YH3A3/cf/rSfFrqTniH2Qem3IUyxPp/RBfWx0+Z/yd03R130r2oruCxi2Mn01/E7KePBr5Ld5fF+wbX7ZefCHyDLsRuA87bRD3H0d0tcW+/nwMGll3Opi9ivpjuDoe7gb/s560E/oHuzeiHdHezHARsS9d3e3df11eBlw7s67x+/XuYchfKXLelu2D8LbquiSuBP57y+L2Z7qLoPcBrePxFzF+gO1u9v38+LuGxC5pPp7vOcTdwzXSPUb/vG/tjXw8c3M8/Gzh7E4/l8n5fD/XbD97Z9DLggYHp19B9evoRXfD/8pR9XdA/n/cCfwvsPrBsRb//H9F1n5221K+hJ+swebeCJKkxdqFIUqMMcElqlAEuSY0ywCWpUVsv5sF22223Wr58+WIecmQ33tj98ufKlbP+gJwkLairr776zqp6wg+PLWqAL1++nLVr1y7mIUd2+OGHA3D55ZcvaR2SlOR70823C0WSGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhq1qN/ElDZXZ625aazt33bkinmqRBqeZ+CS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDVq1gBP8pQkX0nytSQ3JPkP/fxdk6xJcnP/d5eFL1eSNGmYM/B/Bn6xqg4EDgKOSvIi4HTg0qraF7i0n5YkLZJZA7w6D/ST2/RDAccA5/fzzweOXYgCJUnTG6oPPMlWSa4DNgJrqurLwDOqagNA/3f3GbZdnWRtkrUTExPzVLYkaagAr6pHquog4NnAoUmeN+wBquqcqlpVVauWLVs2YpmSpKnmdBdKVd0DXA4cBdyRZA+A/u/G+S5OkjSzYe5CWZZk5378qcArgG8BnwZO6lc7CfjUAtUoSZrGMD8nuwdwfpKt6AL/wqq6JMmXgAuTvBH4PnD8AtYpSZpi1gCvqq8DB08z/4fAEQtRlCRpdn4TU5IaZYBLUqMMcElqlP8npjQPxvk/Nf3/NDUqz8AlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjfKLPNqs+IUYaXiegUtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY2aNcCT7JXksiTrktyQ5K39/DOT3J7kun545cKXK0maNMyPWT0MvL2qrkmyI3B1kjX9srOq6r0LV54kaSazBnhVbQA29OP3J1kHPGuhC5Mkbdqc+sCTLAcOBr7cz3pLkq8nOS/JLjNsszrJ2iRrJyYmxqtWkvRTQwd4kh2ATwKnVtV9wAeBfYCD6M7Q3zfddlV1TlWtqqpVy5YtG79iSRIwZIAn2YYuvD9WVRcDVNUdVfVIVT0KfAg4dOHKlCRNNcxdKAE+DKyrqvcPzN9jYLXjgOvnvzxJ0kyGuQvlMOBE4BtJruvnnQGckOQgoIBbgTctQH2SpBkMcxfKVUCmWfSZ+S9HkjQsv4kpSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVHD/BaK1ISz1ty01CVIi8ozcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1atYAT7JXksuSrEtyQ5K39vN3TbImyc39310WvlxJ0qRhzsAfBt5eVfsBLwJ+N8n+wOnApVW1L3BpPy1JWiSzBnhVbaiqa/rx+4F1wLOAY4Dz+9XOB45doBolSdOYUx94kuXAwcCXgWdU1QboQh7YfYZtVidZm2TtxMTEmOVKkiYNHeBJdgA+CZxaVfcNu11VnVNVq6pq1bJly0apUZI0jaECPMk2dOH9saq6uJ99R5I9+uV7ABsXpkRJ0nSGuQslwIeBdVX1/oFFnwZO6sdPAj41/+VJkmYyzP+JeRhwIvCNJNf1884A3g1cmOSNwPeB4xekQknStGYN8Kq6CsgMi4+Y33IkScPym5iS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGzRrgSc5LsjHJ9QPzzkxye5Lr+uGVC1umJGmqYc7A/wY4apr5Z1XVQf3wmfktS5I0m1kDvKq+ANy1CLVIkuZgnD7wtyT5et/FsstMKyVZnWRtkrUTExNjHE6SNGjUAP8gsA9wELABeN9MK1bVOVW1qqpWLVu2bMTDSZKmGinAq+qOqnqkqh4FPgQcOr9lSZJms/UoGyXZo6o29JPHAddvan1JC+OsNTeNvO3bjlwxj5VoKcwa4EkuAA4HdkuyHngXcHiSg4ACbgXetHAlSpKmM2uAV9UJ08z+8ALUIkmaA7+JKUmNMsAlqVEjXcTUk9s4F8bAi2PSYvEMXJIa5Rm45t24Z/BbGh8vjcozcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjZo1wJOcl2RjkusH5u2aZE2Sm/u/uyxsmZKkqYY5A/8b4Kgp804HLq2qfYFL+2lJ0iKaNcCr6gvAXVNmHwOc34+fDxw7v2VJkmYzah/4M6pqA0D/d/f5K0mSNIwFv4iZZHWStUnWTkxMLPThJGmLMWqA35FkD4D+78aZVqyqc6pqVVWtWrZs2YiHkyRNNWqAfxo4qR8/CfjU/JQjSRrWMLcRXgB8CViZZH2SNwLvBo5McjNwZD8tSVpEW8+2QlWdMMOiI+a5FknSHPhNTElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVGz/haK2nTWmpuWugRJC8wzcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQobyOUtlDj3Gr6tiNXzGMlGpVn4JLUKANckhplgEtSowxwSWqUAS5JjRrrLpQktwL3A48AD1fVqvkoSpI0u/m4jfBfVdWd87AfSdIceB+4pDkb9+eKvY98fozbB17A55NcnWT1dCskWZ1kbZK1ExMTYx5OkjRp3AA/rKoOAY4GfjfJy6euUFXnVNWqqlq1bNmyMQ8nSZo0VoBX1Q/6vxuBvwMOnY+iJEmzGznAk2yfZMfJceCXgOvnqzBJ0qaNcxHzGcDfJZncz8er6rPzUpUkaVYjB3hVfQc4cB5rkSTNgd/ElKRGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQo/0OHzdi4P5ov6cnNM3BJapQBLkmNsgtF0qIbp3vQ/0/zMQb4ArMfW9JCsQtFkhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaNdY3MZMcBfwFsBVwblW9e16q2sz4bUpp8+HX8B8z8hl4kq2A/wocDewPnJBk//kqTJK0aeOcgR8K3FJV3wFI8gngGOCb81HYVIt9Frz+7oeW5LiSFs5Svp4X4ux/nAB/FnDbwPR64IVTV0qyGljdTz6Q5MYxjrnYdjvtl1beudRFLJLdgC2hrVtKO8G2blZOG2/zvaebOU6AZ5p59YQZVecA54xxnCWTZG1VrVrqOhbDltLWLaWdYFu3BOPchbIe2Gtg+tnAD8YrR5I0rHEC/KvAvkl+Lsm2wG8An56fsiRJsxm5C6WqHk7yFuBzdLcRnldVN8xbZZuHJrt+RrSltHVLaSfY1ie9VD2h21qS1AC/iSlJjTLAJalRW2yAJzkqyY1Jbkly+gzrHJ7kuiQ3JLliyrKtklyb5JLFqXg047Qzyc5JLkryrSTrkrx48SqfuzHb+rZ+3vVJLkjylMWrfG5ma2eSP+zbeF3fnkeS7DrMtpubUduaZK8kl/X/bm9I8talqH/BVdUWN9BddP028BxgW+BrwP5T1tmZ7lul/6Kf3n3K8tOAjwOXLHV7FqqdwPnAKf34tsDOS92mhWgr3ZfSvgs8tZ++EDh5qds0ajunrP8q4J9G2XaphzHbugdwSD++I3DT5tzWUYct9Qz8pz8DUFU/BiZ/BmDQbwIXV9X3Aapq4+SCJM8GfgU4d5HqHdXI7UyyE/By4MP9/B9X1T2LVfgIxnpO6e7IemqSrYHt2Hy/0zBMOwedAFww4rZLbeS2VtWGqrqmH78fWEf3Rv2ksqUG+HQ/AzD1yV0B7JLk8iRXJ3ndwLI/B94BPLqgVY5vnHY+B5gAPtJ3FZ2bZPuFL3lkI7e1qm4H3gt8H9gA3FtVn1+EmkcxTDsBSLIdcBTwybluu5kYp62Dy5YDBwNfnv8Sl9aWGuDD/AzA1sDz6c60fxn4oyQrkvwqsLGqrl7gGufDyO3s5x8CfLCqDgYeBDbnPtNxntNd6M7sfg7YE9g+yW8tZLFjGOonLHqvAv53Vd01wrabg3Ha2u0g2YEu1E+tqvvmub4lN9bvgTdsmJ8BWA/cWVUPAg8m+QJwIF2o/VqSVwJPAXZK8tGq2hxf8OO080pgfVVNnrVcxOYd4OO0FeC7VTUBkORi4CXARxe25JHM5ScsfoPHuk/muu3mYJy2kmQbuvD+WFVdvCAVLrWl7oRfioHujes7dGdckxdHDpiyzn7Apf262wHXA8+bss7hbN4XMcdqJ12Ir+zHzwTes9RtWoi20v2K5g39vNBdvP29pW7TqO3s13sacBew/Vy33VyGMdsa4L8Df77U7VjIYYs8A68ZfgYgyZv75WdX1boknwW+TtfXfW5VXb90Vc/dPLTz94CP9b918x3g9YvfiuGM29YkFwHXAA8D17KZfjV7mHb2qx4HfL66Txub3HZxWzC8cdoKHAacCHwjyXX9vDOq6jOLU/3i8Kv0ktSoLfUipiQ1zwCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5Jjfr/fUumErxll1oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(test_values, bins=20, alpha=0.5, density=True)\n",
    "plt.axvline(test_statistic(data), color='black', label = 'Test statistic for real data')\n",
    "plt.title('P-value of test statistic: %.3f' % p_value)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}