{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ddcc3d99-2ca7-47c3-b0ea-fb61220b0c8b",
   "metadata": {},
   "source": [
    "# Lecture 4\n",
    "\n",
    "## RCT for lower back pain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "800bd383-f248-42ca-ab33-1aa07f5c1777",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from matplotlib.pyplot import subplots\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11023792-ab07-4945-a2be-2cf08c238790",
   "metadata": {},
   "outputs": [],
   "source": [
    "bta = pd.read_csv('bta.csv')\n",
    "bta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b17d345-a7f0-4678-bc46-814b516ba3eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "means = bta.groupby('Group')['Result'].mean()\n",
    "means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95e86124-b26a-48a6-babe-b6b3e8962805",
   "metadata": {},
   "outputs": [],
   "source": [
    "obout = pd.read_csv('observed_outcomes.csv')\n",
    "obout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b04ec9d2-8c6d-4653-8bf9-3f512a75fad8",
   "metadata": {},
   "outputs": [],
   "source": [
    "bta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6dcad476-25d8-42c3-9df0-8d561fb15118",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_repetitions = 20000\n",
    "L = []\n",
    "for i in range(num_repetitions):\n",
    "    shuffle = bta.sample(frac=1.0, replace=False)['Group']\n",
    "    shuffled_bta = bta.assign(Shuffled_groups=shuffle.reset_index().iloc[:, 1])\n",
    "    means = shuffled_bta.groupby('Shuffled_groups')['Result'].mean().reset_index()\n",
    "    L.append(means.iloc[0, 1] - means.iloc[1, 1])\n",
    "L = np.array(L)\n",
    "L = np.abs(L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79ed76ca-82b5-4c6f-88d8-daed19d27ecc",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = subplots(figsize = (4, 4))\n",
    "ax.hist(L, bins = np.arange(0, 0.7, 0.1));\n",
    "ax.scatter(0.475, 0, color='red');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d047b74a-f2f2-494a-bd91-b5c8ca852d11",
   "metadata": {},
   "outputs": [],
   "source": [
    "empirical_pvalue = np.sum(L >= 0.475)/num_repetitions\n",
    "empirical_pvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4510f46-f66c-4d89-b1a6-0603167725ba",
   "metadata": {},
   "source": [
    "## Deflategate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e4cc17f-82cb-471a-9bf1-601eb2a8be84",
   "metadata": {},
   "outputs": [],
   "source": [
    "balls = pd.read_csv('deflategate.csv')\n",
    "balls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55ea6dc7-7a02-4227-9038-423a9ac799dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "average_balls = balls.assign(Combined = (balls['Blakeman']+balls['Prioleau'])/2)\n",
    "average_balls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46b87b6e-fe43-4342-a3b8-6f15d1910350",
   "metadata": {},
   "outputs": [],
   "source": [
    "L = np.ones(11)*12.5\n",
    "L = np.append(L, np.ones(4)*13)\n",
    "L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "548b3fc2-2e09-4e29-a038-4d12265e9668",
   "metadata": {},
   "outputs": [],
   "source": [
    "balls2 = average_balls.assign(Drop = L - average_balls['Combined'])\n",
    "balls2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4970be7-3cde-4d20-9909-53ead959c6f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "means = balls2.groupby('Team')['Drop'].mean().reset_index()\n",
    "means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "cf104c07-1a8b-4812-8254-e1ca28708ea6",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_repetitions = 30000\n",
    "L = []\n",
    "for i in range(num_repetitions):\n",
    "    shuffle = balls2.sample(frac=1.0, replace=False)['Team']\n",
    "    shuffled_balls = balls2.assign(Shuffled_teams=shuffle.reset_index().iloc[:, 1])\n",
    "    means = shuffled_balls.groupby('Shuffled_teams')['Drop'].mean().reset_index()\n",
    "    L.append(means.iloc[1, 1] - means.iloc[0, 1])\n",
    "\n",
    "L = np.array(L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "32142fd4-2b03-4355-84fb-7b5ed40ac376",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAFfCAYAAABa/eebAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAkTUlEQVR4nO3dCXQUVb7H8X8CJKxJWBOQxTgoi2wCCozACGRYjDMqOCOCgLINGFB2yAwiokcQRIRhiQoazhkYBI+4gGwGwYU9DgIBMqDhACKgIgQQAoF653/fq37dMSBZSCe3v59z6lS663b17UrnV9W3/l0JchzHEQCANYL93QEAQP4i2AHAMgQ7AFiGYAcAyxDsAGAZgh0ALEOwA4Blioulrl69KseOHZNy5cpJUFCQv7sDAHmmXzs6e/asVKtWTYKDgwMv2DXUa9So4e9uAEC+O3LkiFSvXj3wgl2P1N0NEBYW5u/uAECepaenmwNWN98CLtjd4RcNdYIdgE1+a3iZk6cAYBmCHQAsQ7ADgGUIdgCwDMEOAJYh2AHAMgQ7AFiGYAcAyxDsAGAZgh0ALEOwA4BlCHYAsIy1FwEDCtKt41YW6PMdmhJboM+HooUjdgCwDMEOAJYh2AHAMgQ7AFiGYAcAyxDsAGAZgh0ALEOwA4BlCHYACPRg/+677+Txxx+XihUrSqlSpaRhw4ayY8cOz3LHcWTChAlStWpVszwmJkYOHDjgs45Tp05Jz549JSwsTCIiIqRfv35y7tw5nza7du2SNm3aSMmSJaVGjRoyderUvLxOAAgYOQr2n3/+We69914pUaKErFq1Svbu3SvTp0+X8uXLe9poAM+aNUsSEhJk69atUqZMGenUqZNcvHjR00ZDPSUlRdatWycrVqyQzz77TAYOHOhZnp6eLh07dpRatWpJcnKyTJs2TSZOnChvvPFGfr1uALBWkKOH2Ddo3Lhx8uWXX8rnn3+e7XJdVbVq1WTkyJEyatQoc9+ZM2ckMjJSEhMTpXv37rJv3z6pX7++bN++XZo3b27arF69Wu6//345evSoefy8efPkH//4hxw/flxCQkI8z/3+++/L/v37b6ivunMIDw83z6+fDICbiWvFoCDcaK7l6CJgH374oTn6/stf/iIbN26UW265RZ566ikZMGCAWZ6WlmbCWIdfXNqJFi1ayObNm02w61yHX9xQV9o+ODjYHOE//PDDpk3btm09oa70eV9++WXzqcH7E4IrIyPDTN4bAIGroIMWKLJDMd9++605mr799ttlzZo1MnjwYHn66adl4cKFZrmGutIjdG96212m8ypVqvgsL168uFSoUMGnTXbr8H6OrCZPnmx2Iu6k4/IAEIhyFOxXr16Vpk2byksvvSR33XWXGRfXo3UdT/e3+Ph48/HEnY4cOeLvLgFA4Q92rXTR8XFv9erVk8OHD5ufo6KizPzEiRM+bfS2u0znJ0+e9FmemZlpKmW822S3Du/nyCo0NNSMOXlPABCIchTsWhGTmprqc99///tfU72ioqOjTfAmJSX5jHXr2HmrVq3MbZ2fPn3aVLu41q9fbz4N6Fi820YrZS5fvuxpoxU0derUyXZ8HQCQy2AfPny4bNmyxQzFHDx4UBYvXmxKEOPi4szyoKAgGTZsmLz44ovmROvu3buld+/eptLloYce8hzhd+7c2QzhbNu2zVTZDBkyxJxY1XaqR48e5sSp1rdrWeQ777wjM2fOlBEjRuSkuwAQkHJUFXP33XfL8uXLzXj2pEmTzBH6a6+9ZurSXWPGjJHz58+b8Xc9Mm/durUpZ9QvGrkWLVpkwrxDhw6mGqZbt26m9t2lJz/Xrl1rdhjNmjWTSpUqmS89ede6AwDyoY69KKGOPbDZXu5IHXtgSr/BXONaMQBgGYIdACxDsAOAZQh2ALAMwQ4AliHYAcAyBDsAWIZgBwDLEOwAYBmCHQAsQ7ADgGUIdgCwDMEOAJYh2AHAMgQ7AFiGYAcAyxDsAGAZgh0ALEOwA4BlCHYAsAzBDgCWIdgBwDIEOwBYhmAHAMsQ7ABgGYIdACxDsAOAZQh2ALAMwQ4AliHYAcAyBDsAWIZgBwDLEOwAYBmCHQAsQ7ADgGUIdgAI5GCfOHGiBAUF+Ux169b1LL948aLExcVJxYoVpWzZstKtWzc5ceKEzzoOHz4ssbGxUrp0aalSpYqMHj1aMjMzfdps2LBBmjZtKqGhoVK7dm1JTEzM6+sEgIBRPKcPuPPOO+WTTz75/xUU//9VDB8+XFauXCnLli2T8PBwGTJkiHTt2lW+/PJLs/zKlSsm1KOiomTTpk3y/fffS+/evaVEiRLy0ksvmTZpaWmmzaBBg2TRokWSlJQk/fv3l6pVq0qnTp3y51WjQN06bqW/uwAElBwHuwa5BnNWZ86ckQULFsjixYulffv25r63335b6tWrJ1u2bJGWLVvK2rVrZe/evWbHEBkZKU2aNJEXXnhBxo4daz4NhISESEJCgkRHR8v06dPNOvTxX3zxhcyYMYNgB4CbMcZ+4MABqVatmtx2223Ss2dPM7SikpOT5fLlyxITE+Npq8M0NWvWlM2bN5vbOm/YsKEJdZeGdXp6uqSkpHjaeK/DbeOu41oyMjLMerwnAAhEOQr2Fi1amPHu1atXy7x588ywSZs2beTs2bNy/Phxc8QdERHh8xgNcV2mdO4d6u5yd9n12mhQX7hw4Zp9mzx5shn+cacaNWrk5KUBQGAOxXTp0sXzc6NGjUzQ16pVS5YuXSqlSpUSf4qPj5cRI0Z4buuOgHAHEIjyVO6oR+d33HGHHDx40Iy7X7p0SU6fPu3TRqti3DF5nWetknFv/1absLCw6+48tIJG23hPABCI8hTs586dk2+++cZUrDRr1sxUt2gViys1NdWMwbdq1crc1vnu3bvl5MmTnjbr1q0zIVy/fn1PG+91uG3cdQAA8jHYR40aJRs3bpRDhw6ZcsWHH35YihUrJo899pgZ1+7Xr58ZDvn000/NydQnn3zSBLJWxKiOHTuaAO/Vq5d8/fXXsmbNGhk/frypfdcjbqVljt9++62MGTNG9u/fL3PnzjVDPVpKCQDI5zH2o0ePmhD/6aefpHLlytK6dWtTyqg/Ky1JDA4ONl9M0ioVrWbRYHbpTmDFihUyePBgE/hlypSRPn36yKRJkzxttNRRa+E1yGfOnCnVq1eX+fPnU+oIADcoyHEcRyykJ0/1U4TW1zPe7l98QSn/HZoS6+8uoBDnGteKAQDLEOwAYBmCHQAsQ7ADgGUIdgCwDMEOAJYh2AHAMgQ7AFiGYAcAyxDsAGAZgh0ALEOwA4BlCHYAsAzBDgCWIdgBwDIEOwBYhmAHAMsQ7ABgGYIdACxDsAOAZQh2ALAMwQ4AliHYAcAyBDsAWIZgBwDLFPd3BwDk3K3jVhbo8x2aElugz4e84YgdACxDsAOAZQh2ALAMwQ4AliHYAcAyBDsAWIZgBwDLEOwAYBmCHQAsk6dgnzJligQFBcmwYcM89128eFHi4uKkYsWKUrZsWenWrZucOHHC53GHDx+W2NhYKV26tFSpUkVGjx4tmZmZPm02bNggTZs2ldDQUKldu7YkJibmpasAEDByHezbt2+X119/XRo1auRz//Dhw+Wjjz6SZcuWycaNG+XYsWPStWtXz/IrV66YUL906ZJs2rRJFi5caEJ7woQJnjZpaWmmTbt27WTnzp1mx9G/f39Zs2ZNbrsLAAEjV8F+7tw56dmzp7z55ptSvnx5z/1nzpyRBQsWyKuvvirt27eXZs2aydtvv20CfMuWLabN2rVrZe/evfKvf/1LmjRpIl26dJEXXnhB5syZY8JeJSQkSHR0tEyfPl3q1asnQ4YMkUceeURmzJiRX68bAKyVq2DXoRY9oo6JifG5Pzk5WS5fvuxzf926daVmzZqyefNmc1vnDRs2lMjISE+bTp06SXp6uqSkpHjaZF23tnHXkZ2MjAyzDu8JAAJRjq/uuGTJEvnqq6/MUExWx48fl5CQEImIiPC5X0Ncl7ltvEPdXe4uu14bDesLFy5IqVKlfvXckydPlueffz6nLwcAAvuI/ciRI/LMM8/IokWLpGTJklKYxMfHm6Egd9K+AkAgylGw61DLyZMnTbVK8eLFzaQnSGfNmmV+1qNqHSc/ffq0z+O0KiYqKsr8rPOsVTLu7d9qExYWlu3RutLqGV3uPQFAIMpRsHfo0EF2795tKlXcqXnz5uZEqvtziRIlJCkpyfOY1NRUU97YqlUrc1vnug7dQbjWrVtngrh+/fqeNt7rcNu46wAA5NMYe7ly5aRBgwY+95UpU8bUrLv39+vXT0aMGCEVKlQwYT106FATyC1btjTLO3bsaAK8V69eMnXqVDOePn78eHNCVo+61aBBg2T27NkyZswY6du3r6xfv16WLl0qK1cW7H+NAYCiKN//NZ6WJAYHB5svJmmlilazzJ0717O8WLFismLFChk8eLAJfN0x9OnTRyZNmuRpo6WOGuJaEz9z5kypXr26zJ8/36wLAHB9QY7jOGIhraAJDw83J1IZbw+s/8+J/Mf/PC1auca1YgDAMgQ7AFiGYAcAyxDsAGAZgh0ALEOwA4BlCHYAsAzBDgCWIdgBwDIEOwBYhmAHAMsQ7ABgGYIdACxDsAOAZQh2ALAMwQ4AliHYAcAyBDsAWIZgBwDLEOwAYBmCHQAsQ7ADgGUIdgCwDMEOAJYh2AHAMgQ7AFiGYAcAyxDsAGAZgh0ALEOwA4BlCHYAsAzBDgCWIdgBwDIEOwBYhmAHAMsQ7AAQyME+b948adSokYSFhZmpVatWsmrVKs/yixcvSlxcnFSsWFHKli0r3bp1kxMnTvis4/DhwxIbGyulS5eWKlWqyOjRoyUzM9OnzYYNG6Rp06YSGhoqtWvXlsTExLy+TgAIGDkK9urVq8uUKVMkOTlZduzYIe3bt5cHH3xQUlJSzPLhw4fLRx99JMuWLZONGzfKsWPHpGvXrp7HX7lyxYT6pUuXZNOmTbJw4UIT2hMmTPC0SUtLM23atWsnO3fulGHDhkn//v1lzZo1+fm6AcBaQY7jOHlZQYUKFWTatGnyyCOPSOXKlWXx4sXmZ7V//36pV6+ebN68WVq2bGmO7h944AET+JGRkaZNQkKCjB07Vn744QcJCQkxP69cuVL27NnjeY7u3bvL6dOnZfXq1Tfcr/T0dAkPD5czZ86YTxfwn1vHrfR3F5BHh6bE+rsLkBvPtVyPsevR95IlS+T8+fNmSEaP4i9fviwxMTGeNnXr1pWaNWuaYFc6b9iwoSfUVadOnUxn3aN+beO9DreNu45rycjIMOvxngAgEOU42Hfv3m3Gz3X8e9CgQbJ8+XKpX7++HD9+3BxxR0RE+LTXENdlSufeoe4ud5ddr40G9YULF67Zr8mTJ5s9mTvVqFEjpy8NAAIz2OvUqWPGvrdu3SqDBw+WPn36yN69e8Xf4uPjzccTdzpy5Ii/uwQAflE8pw/Qo3KtVFHNmjWT7du3y8yZM+XRRx81J0V1LNz7qF2rYqKioszPOt+2bZvP+tyqGe82WStp9LaOJ5UqVeqa/dJPEDoBQKDLcx371atXzfi2hnyJEiUkKSnJsyw1NdWUN+oYvNK5DuWcPHnS02bdunUmtHU4x23jvQ63jbsOAEA+HrHrcEeXLl3MCdGzZ8+aChitOddSRB3X7tevn4wYMcJUymhYDx061ASyVsSojh07mgDv1auXTJ061Yynjx8/3tS+u0fbOm4/e/ZsGTNmjPTt21fWr18vS5cuNZUyAIB8DnY90u7du7d8//33Jsj1y0oa6n/84x/N8hkzZkhwcLD5YpIexWs1y9y5cz2PL1asmKxYscKMzWvglylTxozRT5o0ydMmOjrahLjWxOsQj9bOz58/36wL+YPyQ8Buea5jL6yoY782gh05RR17gNSxAwAKJ4IdACxDsAOAZQh2ALAMwQ4AliHYAcAyBDsAWIZgBwDLEOwAYBmCHQAC/bK9AAJPQV+GgksY5A1H7ABgGYIdACxDsAOAZQh2ALAMwQ4AliHYAcAyBDsAWIZgBwDLEOwAYBmCHQAsQ7ADgGUIdgCwDMEOAJYh2AHAMgQ7AFiGYAcAyxDsAGAZgh0ALEOwA4BlCHYAsAzBDgCWIdgBwDIEOwBYhmAHgEAO9smTJ8vdd98t5cqVkypVqshDDz0kqampPm0uXrwocXFxUrFiRSlbtqx069ZNTpw44dPm8OHDEhsbK6VLlzbrGT16tGRmZvq02bBhgzRt2lRCQ0Oldu3akpiYmJfXCQABI0fBvnHjRhPaW7ZskXXr1snly5elY8eOcv78eU+b4cOHy0cffSTLli0z7Y8dOyZdu3b1LL9y5YoJ9UuXLsmmTZtk4cKFJrQnTJjgaZOWlmbatGvXTnbu3CnDhg2T/v37y5o1a/LrdQOAtYIcx3Fy++AffvjBHHFrgLdt21bOnDkjlStXlsWLF8sjjzxi2uzfv1/q1asnmzdvlpYtW8qqVavkgQceMIEfGRlp2iQkJMjYsWPN+kJCQszPK1eulD179nieq3v37nL69GlZvXr1DfUtPT1dwsPDTZ/CwsJy+xKtdOu4lf7uAnBdh6bE+rsLhdKN5lqexth15apChQpmnpycbI7iY2JiPG3q1q0rNWvWNMGudN6wYUNPqKtOnTqZDqekpHjaeK/DbeOuIzsZGRlmHd4TAASiXAf71atXzRDJvffeKw0aNDD3HT9+3BxxR0RE+LTVENdlbhvvUHeXu8uu10bD+sKFC9cc/9c9mTvVqFEjty8NAAIz2HWsXYdKlixZIoVBfHy8+QThTkeOHPF3lwDAL4rn5kFDhgyRFStWyGeffSbVq1f33B8VFWVOiupYuPdRu1bF6DK3zbZt23zW51bNeLfJWkmjt3VMqVSpUtn2SatndAKAQJejI3Y9z6qhvnz5clm/fr1ER0f7LG/WrJmUKFFCkpKSPPdpOaSWN7Zq1crc1vnu3bvl5MmTnjZaYaOhXb9+fU8b73W4bdx1AADy6Yhdh1+04uWDDz4wtezumLiOaeuRtM779esnI0aMMCdUNayHDh1qAlkrYpSWR2qA9+rVS6ZOnWrWMX78eLNu94h70KBBMnv2bBkzZoz07dvX7ESWLl1qKmUAAPl4xD5v3jwzfn3fffdJ1apVPdM777zjaTNjxgxTzqhfTNISSB1Wee+99zzLixUrZoZxdK6B//jjj0vv3r1l0qRJnjb6SUBDXI/SGzduLNOnT5f58+ebyhgAwE2sYy/MqGO/NurYUdhRx+7HOnYAQOFDsAOAZQh2ALAMwQ4AliHYAcAyBDsAWIZgBwDLEOwAYBmCHQAsQ7ADgGUIdgCwDMEOAJYh2AHAMgQ7AFiGYAcAyxDsAGAZgh0ALEOwA4BlCHYAsAzBDgCWIdgBwDIEOwBYhmAHAMsQ7ABgGYIdACxDsAOAZQh2ALAMwQ4AliHYAcAyxf3dAYjcOm6lv7sAwCIcsQOAZQh2ALAMwQ4AliHYAcAyBDsAWIZgB4BAD/bPPvtM/vSnP0m1atUkKChI3n//fZ/ljuPIhAkTpGrVqlKqVCmJiYmRAwcO+LQ5deqU9OzZU8LCwiQiIkL69esn586d82mza9cuadOmjZQsWVJq1KghU6dOze1rBICAkuNgP3/+vDRu3FjmzJmT7XIN4FmzZklCQoJs3bpVypQpI506dZKLFy962miop6SkyLp162TFihVmZzFw4EDP8vT0dOnYsaPUqlVLkpOTZdq0aTJx4kR54403cvs6ASBgBDl6iJ3bBwcFyfLly+Whhx4yt3VVeiQ/cuRIGTVqlLnvzJkzEhkZKYmJidK9e3fZt2+f1K9fX7Zv3y7Nmzc3bVavXi3333+/HD161Dx+3rx58o9//EOOHz8uISEhps24cePMp4P9+/ffUN905xAeHm6eXz8ZFGZ8QQnwdWhKrL+7UCjdaK7l6xh7WlqaCWMdfnFpJ1q0aCGbN282t3Wuwy9uqCttHxwcbI7w3TZt27b1hLrSo/7U1FT5+eefs33ujIwM86K9JwAIRPka7BrqSo/Qveltd5nOq1Sp4rO8ePHiUqFCBZ822a3D+zmymjx5stmJuJOOywNAILKmKiY+Pt58PHGnI0eO+LtLAFD0gz0qKsrMT5w44XO/3naX6fzkyZM+yzMzM02ljHeb7Nbh/RxZhYaGmjEn7wkAAlG+Bnt0dLQJ3qSkJM99OtatY+etWrUyt3V++vRpU+3iWr9+vVy9etWMxbtttFLm8uXLnjZaQVOnTh0pX758fnYZAKyT48v2ar35wYMHfU6Y7ty504yR16xZU4YNGyYvvvii3H777Sbon332WVPp4lbO1KtXTzp37iwDBgwwJZEa3kOGDDEVM9pO9ejRQ55//nlT3z527FjZs2ePzJw5U2bMmJGfrx1AIeWPSrFDFlXi5DjYd+zYIe3atfPcHjFihJn36dPHlDSOGTPG1LprXboembdu3dqUM+oXjVyLFi0yYd6hQwdTDdOtWzdT++7Sk59r166VuLg4adasmVSqVMl86cm71h0AcBPq2Asz6tgB2HbE7pc6dgCA/xHsAGAZgh0ALEOwA4BlCHYAsAzBDgCWIdgBwDIEOwBYhmAHAMsQ7ABgGYIdACxDsAOAZQh2ALAMwQ4AliHYAcAyBDsAWIZgBwDLEOwAYBmCHQAsQ7ADgGUIdgCwTHF/d6AwunXcSn93AQByjSN2ALAMwQ4AliHYAcAyBDsAWIZgBwDLUBUDAFLw1XCHpsTetHVzxA4AliHYAcAyBDsAWIZgBwDLEOwAYBmCHQAsQ7kjAORC8NUrcs/RFKly7mc5Wba8bKt+p1wNLiaFAcEOADnUKXWTPJf0hlQ7+6PnvmPlKsnzHQbKmjq/F38r1EMxc+bMkVtvvVVKliwpLVq0kG3btvm7SwACXKfUTTLv/ZckyivUld7W+3W5vxXaYH/nnXdkxIgR8txzz8lXX30ljRs3lk6dOsnJkyf93TUAATz88lzSG//7c9Zl/zfX5drOnwrtUMyrr74qAwYMkCeffNLcTkhIkJUrV8pbb70l48aN+1X7jIwMM7nOnDlj5unp6Tl+7qsZv+Sp7wDs1OzoHil79kc5d502uvyutGTZXr3BddeVm2xyH+M4zvUbOoVQRkaGU6xYMWf58uU+9/fu3dv585//nO1jnnvuOX2lTExMTI7t05EjR66boYXyiP3HH3+UK1euSGRkpM/9env//v3ZPiY+Pt4M3biuXr0qp06dkooVK0pQUFCu+6J7yBo1asiRI0ckLCxMioKi1uei1t+i2Oei1t+i2Of0AuivHqmfPXtWqlWrdt12hTLYcyM0NNRM3iIiIvJt/fqLKgpvrqLc56LW36LY56LW36LY57Cb3N/w8PCiefK0UqVKUqxYMTlx4oTP/Xo7KirKb/0CgKKgUAZ7SEiINGvWTJKSknyGVvR2q1at/No3ACjsCu1QjI6X9+nTR5o3by733HOPvPbaa3L+/HlPlUxB0eEdLbnMOsxTmBW1Phe1/hbFPhe1/hbFPocWov4G6RlUKaRmz54t06ZNk+PHj0uTJk1k1qxZ5otKAIAiGuwAAEvG2AEAuUewA4BlCHYAsAzBDgCWIdhFzKUHevbsab4tpt9W7devn5w7d+3L/Bw6dMhcpiC7admyZZ522S1fsmRJgfdX3Xfffb/qy6BBg3zaHD58WGJjY6V06dJSpUoVGT16tGRmZua5v7nps7YfOnSo1KlTR0qVKiU1a9aUp59+2nNxt/zexjm9RLT+nuvWrWvaN2zYUD7++GOf5VqTMGHCBKlatarpf0xMjBw4cCDH/cqvPr/55pvSpk0bKV++vJm0P1nbP/HEE7/alp07d/ZLfxMTE3/VF31cYd7G92XzN6aT/k0V1Db2yMdrdxVZnTt3dho3buxs2bLF+fzzz53atWs7jz322DXbZ2ZmOt9//73P9Pzzzztly5Z1zp4962mnm/ftt9/2aXfhwoUC76/6wx/+4AwYMMCnL2fOnPF5TQ0aNHBiYmKc//znP87HH3/sVKpUyYmPj89zf3PT5927dztdu3Z1PvzwQ+fgwYNOUlKSc/vttzvdunXzaZcf23jJkiVOSEiI89ZbbzkpKSlmO0VERDgnTpzItv2XX35pLlI3depUZ+/evc748eOdEiVKmD67pkyZ4oSHhzvvv/++8/XXX5uL10VHR+fL7z83fe7Ro4czZ84c87vdt2+f88QTT5j+HT161NOmT58+5vfkvS1PnTrll/7q7zQsLMynL8ePH/dpU9i28U8//eTT3z179pj3ib6WgtjG3gI+2PUPU8Nh+/btnvtWrVrlBAUFOd99990Nr6dJkyZO3759fe7T9Wa9QqW/+qvB/swzz1xzuQZ5cHCwzx/PvHnzzB+XXm3TH33OaunSpeYP7fLly/m6je+55x4nLi7Oc/vKlStOtWrVnMmTJ2fb/q9//asTGxvrc1+LFi2cv/3tb+bnq1evOlFRUc60adM8y0+fPu2EhoY6//73v/PU19z2OSvdkZcrV85ZuHChT+g8+OCD+dK/vPZXw1BD+1qKwjaeMWOG2cbnzp0rkG3sLeCHYjZv3myGBvQbri79SBccHCxbt269oXUkJyfLzp07zfBCVnFxcebaN/rtWb2WfF6/NpCX/i5atMj0pUGDBuZqmL/88ovPenVIwfuKmvqPTfSKdSkpKX7rszcdhtGhnOLFi+fbNr506ZL5/Wl/XNovva39vtbr8W7vbiu3fVpamvlSnXcbvXCTfpS/1jpzIjd9zkp/95cvX5YKFSr43L9hwwYzDKdDYIMHD5affvrJb/3VobpatWqZKyY++OCDPu/DorCNFyxYIN27d5cyZcrc9G1cZC4pUFD0zaEb2ZsGh77hddmN/gLr1asnv/+97/86nDRpkrRv396MWa9du1aeeuop82bVseKC7m+PHj3MH4le7nPXrl0yduxYSU1Nlffee8+z3uwuk+wu8/c21ks5v/DCCzJw4MB83ca5uUT0tbaV+1rc+fXa5EVu+pyV/v71veAdXDrW27VrV4mOjpZvvvlG/v73v0uXLl1MkOlF+Qqyvxp6upNu1KiR2aG/8sor5u9Lw7169eqFfhtv27ZN9uzZY7LB283axgET7Ppfll5++eXrttm3b1+en+fChQuyePFiefbZZ3+1zPu+u+66y1zrRi+RkF3o3Oz+egeiHpnrCacOHTqYN9fvfve7Qr2N9VODnoCqX7++TJw4MdfbGP9rypQp5gSzHjl6n5DUo0vv94iGqr43tJ2+VwqSXuzP+4J/Gup68PT666+bHXxht2DBArMN9VOkt4LaxtYG+8iRI80Z6Ou57bbbzGWAs/4fVa0E0aqMG7lE8Lvvvms+1vbu3fs32+rHRH1T6r/wy3qhoILqr3df1MGDB80bSx+b9Yy/e9nka623IPqs/1RAj3LKlSsny5cvlxIlSuR6G+fXJaL1/uu1d+d6n+5AvdvoNY/8eVlrPfLVYP/kk09MqPzW706fS98jeQmd/LgMt/7edcetfSns2/j8+fNmx6mfJn9Lfm3jX3ECnHtib8eOHZ771qxZc8Mn9vSkZNZKjWt58cUXnfLly/u1v64vvvjCrEerCbxPnnqf8X/99dfNydOLFy/6pc9atdOyZUuzjc+fP3/TtrGeJBsyZIjPSbJbbrnluidPH3jgAZ/7WrVq9auTp6+88orPa8nvE3s56bN6+eWXze9z8+bNN/Qc+u/X9Hf0wQcf+KW/WU/21qlTxxk+fHih3sbuiV/tx48//ugU5Db2FvDBrrT86K677nK2bt1qAk/L6rxL8bQkTN9UutzbgQMHzC9FKzyy0jK9N99805TAabu5c+c6pUuXdiZMmFDg/dVywUmTJplgTUtLM2+i2267zWnbtu2vyh07duzo7Ny501m9erVTuXLlfC13zEmf9Y9UK00aNmxo+u9dHqZ9zc9trGVt+oeYmJhodkIDBw40ZW1uhVCvXr2ccePG+ZQ7Fi9e3ISKlg7q/9vNrtxR16HbeteuXaYSIr9L8XLSZ+2PVhS9++67PtvSLc/V+ahRo0zo63vkk08+cZo2bWp+T3ndseemv1o+rDv/b775xklOTna6d+/ulCxZ0pQdFtZt7GrdurXz6KOPOlnd7G3sjWD/v/pTDRmtQ9cjmieffNKnHl1/CXrE+emnn/o8TkOvRo0aZk+elYa9lkDqOsuUKWNquBMSErJte7P7e/jwYRPiFSpUMG9UrSEfPXq0Tx27OnTokNOlSxenVKlSpoZ95MiRPqWFBdlnnV/rH/lq2/zexv/85z+dmjVrmvDTIzWtt3fpJwYtU8taennHHXeY9nfeeaezcuVKn+V6RPnss886kZGRZpt36NDBSU1NzXG/8qvPtWrVynZb6k5J/fLLL2anrjtz3Ulpe63bzlo7XlD9HTZsmKetbsP777/f+eqrrwr1Nlb79+8323Xt2rVOVgWxjV1cthcALBPwdewAYBuCHQAsQ7ADgGUIdgCwDMEOAJYh2AHAMgQ7AFiGYAcAyxDsAGAZgh0ALEOwA4DY5X8Amf6YHShLEWQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots(figsize = (4, 4))\n",
    "ax.hist(L);\n",
    "ax.scatter(0.73352, 0, color='red');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "7dbffa74-b431-4929-b6d4-cf7650fbbb4d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.float64(0.0033666666666666667)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "empirical_pvalue = np.sum(L >= 0.73352)/num_repetitions\n",
    "empirical_pvalue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cff66fbf-f517-4652-ae32-c29491b5c087",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}