empirical_power_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Empirical power analysis\n",
    "\n",
    "tdsmith@mozilla.com, 2018-07-17\n",
    "\n",
    "This notebook demonstrates how to run an empirical power analysis on data from a Poisson distribution, which is intended to approximate the distribution of indices of search bar results chosen by Firefox users."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from IPython.display import display\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import scipy.stats as stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate a parent distribution that we'll sample from iteratively.\n",
    "\n",
    "For a real analysis, the parent distribution should be a large sample of the actual observed telemetry data.\n",
    "\n",
    "$\\mu$ is the location parameter and arithmetic mean of the Poisson distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 3 1 2 0 0 2 1 0 4]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAD8CAYAAABZ/vJZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFWRJREFUeJzt3X+s3XWd5/Hna1pRRlcpcrdh25qSsRlTSSzahc662biwAwU3WyZRA8lIYxo7G2FXN2bX6j/MqCSY7MguiZIwQ4fiuiJBJzRap9MgiZk/+HHRDlDQcJcfQ7uFdig/dI244Hv/OJ/uHDqn995+bm/PpX0+km/O97y/n8/38/kC6Yvz/X7OaaoKSZKO1W+NewKSpDcmA0SS1MUAkSR1MUAkSV0MEElSFwNEktTFAJEkdTFAJEldDBBJUpfF457A8XbWWWfVypUrxz0NSXpDefDBB/++qiaOpc9JFyArV65kcnJy3NOQpDeUJE8fax9vYUmSuhggkqQuBogkqYsBIknqYoBIkroYIJKkLgaIJKmLASJJ6mKASJK6nHTfRJ+LlVu+P+4pnHBPXf/hcU9B0huUn0AkSV0MEElSFwNEktTFAJEkdTFAJEldZgyQJG9Jcn+Sv02yJ8mftPqtSZ5Msrtta1o9SW5MMpXkoSTvHzrXxiSPt23jUP0DSR5ufW5MklY/M8mu1n5XkiXH/x+BJKnHbD6BvAJcWFXvA9YA65Osa8f+c1WtadvuVrsUWNW2zcBNMAgD4FrgAuB84NqhQLgJ+ORQv/WtvgW4u6pWAXe395KkBWDGAKmBX7S3b2pbTdNlA3Bb63cvcEaSs4FLgF1VdaiqXgB2MQijs4G3V9W9VVXAbcDlQ+fa1va3DdUlSWM2q2cgSRYl2Q0cYBAC97VD17XbVDckeXOrLQOeGeq+t9Wmq+8dUQdYWlX72/6zwNLZXZYkab7NKkCq6rWqWgMsB85Pci7weeA9wD8HzgQ+N2+zHMyhOMonnySbk0wmmTx48OB8TkOS1BzTKqyqehG4B1hfVfvbbapXgL9g8FwDYB+wYqjb8labrr58RB3guXaLi/Z64Cjzurmq1lbV2omJiWO5JElSp9mswppIckbbPx34feCnQ3+wh8GziUdal+3AVW011jrgpXYbaidwcZIl7eH5xcDOduzlJOvaua4C7ho61+HVWhuH6pKkMZvNjymeDWxLsohB4NxRVd9L8sMkE0CA3cC/b+13AJcBU8AvgU8AVNWhJF8CHmjtvlhVh9r+p4BbgdOBH7QN4HrgjiSbgKeBj/VeqCTp+JoxQKrqIeC8EfULj9K+gKuPcmwrsHVEfRI4d0T9eeCimeYoSTrx/Ca6JKmLASJJ6mKASJK6GCCSpC4GiCSpiwEiSepigEiSuhggkqQuBogkqYsBIknqYoBIkroYIJKkLgaIJKmLASJJ6mKASJK6GCCSpC4GiCSpiwEiSepigEiSuswYIEnekuT+JH+bZE+SP2n1c5Lcl2QqybeTnNbqb27vp9rxlUPn+nyr/yzJJUP19a02lWTLUH3kGJKk8ZvNJ5BXgAur6n3AGmB9knXAV4AbqurdwAvAptZ+E/BCq9/Q2pFkNXAF8F5gPfD1JIuSLAK+BlwKrAaubG2ZZgxJ0pjNGCA18Iv29k1tK+BC4M5W3wZc3vY3tPe04xclSavfXlWvVNWTwBRwftumquqJqvo1cDuwofU52hiSpDGb1TOQ9klhN3AA2AX8L+DFqnq1NdkLLGv7y4BnANrxl4B3DteP6HO0+junGUOSNGazCpCqeq2q1gDLGXxieM+8zuoYJdmcZDLJ5MGDB8c9HUk6JRzTKqyqehG4B/g94Iwki9uh5cC+tr8PWAHQjr8DeH64fkSfo9Wfn2aMI+d1c1Wtraq1ExMTx3JJkqROs1mFNZHkjLZ/OvD7wGMMguQjrdlG4K62v729px3/YVVVq1/RVmmdA6wC7gceAFa1FVenMXjQvr31OdoYkqQxWzxzE84GtrXVUr8F3FFV30vyKHB7ki8DPwFuae1vAb6RZAo4xCAQqKo9Se4AHgVeBa6uqtcAklwD7AQWAVurak871+eOMoYkacxmDJCqegg4b0T9CQbPQ46s/wr46FHOdR1w3Yj6DmDHbMeQJI2f30SXJHUxQCRJXQwQSVIXA0SS1MUAkSR1MUAkSV0MEElSFwNEktTFAJEkdTFAJEldDBBJUhcDRJLUxQCRJHUxQCRJXQwQSVIXA0SS1MUAkSR1MUAkSV0MEElSlxkDJMmKJPckeTTJniSfbvU/TrIvye62XTbU5/NJppL8LMklQ/X1rTaVZMtQ/Zwk97X6t5Oc1upvbu+n2vGVx/PiJUn9ZvMJ5FXgs1W1GlgHXJ1kdTt2Q1WtadsOgHbsCuC9wHrg60kWJVkEfA24FFgNXDl0nq+0c70beAHY1OqbgBda/YbWTpK0AMwYIFW1v6p+3PZ/DjwGLJumywbg9qp6paqeBKaA89s2VVVPVNWvgduBDUkCXAjc2fpvAy4fOte2tn8ncFFrL0kas2N6BtJuIZ0H3NdK1yR5KMnWJEtabRnwzFC3va12tPo7gRer6tUj6q87Vzv+UmsvSRqzWQdIkrcB3wE+U1UvAzcBvwOsAfYDfzovM5zd3DYnmUwyefDgwXFNQ5JOKbMKkCRvYhAe36yq7wJU1XNV9VpV/Qb4Mwa3qAD2ASuGui9vtaPVnwfOSLL4iPrrztWOv6O1f52qurmq1lbV2omJidlckiRpjmazCivALcBjVfXVofrZQ83+AHik7W8HrmgrqM4BVgH3Aw8Aq9qKq9MYPGjfXlUF3AN8pPXfCNw1dK6Nbf8jwA9be0nSmC2euQkfBD4OPJxkd6t9gcEqqjVAAU8BfwRQVXuS3AE8ymAF19VV9RpAkmuAncAiYGtV7Wnn+xxwe5IvAz9hEFi0128kmQIOMQgdSdICMGOAVNXfAKNWPu2Yps91wHUj6jtG9auqJ/iHW2DD9V8BH51pjpKkE89vokuSuhggkqQuBogkqYsBIknqYoBIkroYIJKkLgaIJKmLASJJ6mKASJK6GCCSpC4GiCSpiwEiSepigEiSuhggkqQuBogkqYsBIknqYoBIkroYIJKkLgaIJKmLASJJ6jJjgCRZkeSeJI8m2ZPk061+ZpJdSR5vr0taPUluTDKV5KEk7x8618bW/vEkG4fqH0jycOtzY5JMN4Ykafxm8wnkVeCzVbUaWAdcnWQ1sAW4u6pWAXe39wCXAqvathm4CQZhAFwLXACcD1w7FAg3AZ8c6re+1Y82hiRpzGYMkKraX1U/bvs/Bx4DlgEbgG2t2Tbg8ra/AbitBu4FzkhyNnAJsKuqDlXVC8AuYH079vaqureqCrjtiHONGkOSNGbH9AwkyUrgPOA+YGlV7W+HngWWtv1lwDND3fa22nT1vSPqTDPGkfPanGQyyeTBgweP5ZIkSZ1mHSBJ3gZ8B/hMVb08fKx9cqjjPLfXmW6Mqrq5qtZW1dqJiYn5nIYkqZlVgCR5E4Pw+GZVfbeVn2u3n2ivB1p9H7BiqPvyVpuuvnxEfboxJEljNptVWAFuAR6rqq8OHdoOHF5JtRG4a6h+VVuNtQ54qd2G2glcnGRJe3h+MbCzHXs5ybo21lVHnGvUGJKkMVs8izYfBD4OPJxkd6t9AbgeuCPJJuBp4GPt2A7gMmAK+CXwCYCqOpTkS8ADrd0Xq+pQ2/8UcCtwOvCDtjHNGJKkMZsxQKrqb4Ac5fBFI9oXcPVRzrUV2DqiPgmcO6L+/KgxJEnj5zfRJUldDBBJUhcDRJLUZTYP0XUSW7nl+2MZ96nrPzyWcSUdP34CkSR1MUAkSV0MEElSFwNEktTFAJEkdTFAJEldDBBJUhcDRJLUxQCRJHUxQCRJXQwQSVIXA0SS1MUAkSR1MUAkSV0MEElSlxkDJMnWJAeSPDJU++Mk+5LsbttlQ8c+n2Qqyc+SXDJUX99qU0m2DNXPSXJfq387yWmt/ub2fqodX3m8LlqSNHez+QRyK7B+RP2GqlrTth0ASVYDVwDvbX2+nmRRkkXA14BLgdXAla0twFfaud4NvABsavVNwAutfkNrJ0laIGYMkKr6EXBolufbANxeVa9U1ZPAFHB+26aq6omq+jVwO7AhSYALgTtb/23A5UPn2tb27wQuau0lSQvAXJ6BXJPkoXaLa0mrLQOeGWqzt9WOVn8n8GJVvXpE/XXnasdfau3/kSSbk0wmmTx48OAcLkmSNFu9AXIT8DvAGmA/8KfHbUYdqurmqlpbVWsnJibGORVJOmV0BUhVPVdVr1XVb4A/Y3CLCmAfsGKo6fJWO1r9eeCMJIuPqL/uXO34O1p7SdIC0BUgSc4eevsHwOEVWtuBK9oKqnOAVcD9wAPAqrbi6jQGD9q3V1UB9wAfaf03AncNnWtj2/8I8MPWXpK0ACyeqUGSbwEfAs5Kshe4FvhQkjVAAU8BfwRQVXuS3AE8CrwKXF1Vr7XzXAPsBBYBW6tqTxvic8DtSb4M/AS4pdVvAb6RZIrBQ/wr5ny1kqTjZsYAqaorR5RvGVE73P464LoR9R3AjhH1J/iHW2DD9V8BH51pfpKk8fCb6JKkLgaIJKmLASJJ6mKASJK6GCCSpC4GiCSpiwEiSepigEiSuhggkqQuBogkqYsBIknqYoBIkroYIJKkLgaIJKmLASJJ6mKASJK6GCCSpC4GiCSpiwEiSeoyY4Ak2ZrkQJJHhmpnJtmV5PH2uqTVk+TGJFNJHkry/qE+G1v7x5NsHKp/IMnDrc+NSTLdGJKkhWE2n0BuBdYfUdsC3F1Vq4C723uAS4FVbdsM3ASDMACuBS4AzgeuHQqEm4BPDvVbP8MYkqQFYMYAqaofAYeOKG8AtrX9bcDlQ/XbauBe4IwkZwOXALuq6lBVvQDsAta3Y2+vqnurqoDbjjjXqDEkSQtA7zOQpVW1v+0/Cyxt+8uAZ4ba7W216ep7R9SnG+MfSbI5yWSSyYMHD3ZcjiTpWM35IXr75FDHYS7dY1TVzVW1tqrWTkxMzOdUJElNb4A8124/0V4PtPo+YMVQu+WtNl19+Yj6dGNIkhaA3gDZDhxeSbURuGuoflVbjbUOeKndhtoJXJxkSXt4fjGwsx17Ocm6tvrqqiPONWoMSdICsHimBkm+BXwIOCvJXgarqa4H7kiyCXga+FhrvgO4DJgCfgl8AqCqDiX5EvBAa/fFqjr8YP5TDFZ6nQ78oG1MM4YkaQGYMUCq6sqjHLpoRNsCrj7KebYCW0fUJ4FzR9SfHzWGJGlh8JvokqQuBogkqYsBIknqYoBIkroYIJKkLgaIJKmLASJJ6mKASJK6GCCSpC4GiCSpiwEiSepigEiSuhggkqQuBogkqYsBIknqYoBIkrrM+BdKSfNh5Zbvj23sp67/8NjGlk4mfgKRJHUxQCRJXeYUIEmeSvJwkt1JJlvtzCS7kjzeXpe0epLcmGQqyUNJ3j90no2t/eNJNg7VP9DOP9X6Zi7zlSQdP8fjE8i/rqo1VbW2vd8C3F1Vq4C723uAS4FVbdsM3ASDwAGuBS4AzgeuPRw6rc0nh/qtPw7zlSQdB/NxC2sDsK3tbwMuH6rfVgP3AmckORu4BNhVVYeq6gVgF7C+HXt7Vd1bVQXcNnQuSdKYzTVACvjrJA8m2dxqS6tqf9t/Flja9pcBzwz13dtq09X3jqhLkhaAuS7j/ZdVtS/JPwV2Jfnp8MGqqiQ1xzFm1MJrM8C73vWu+R5OksQcP4FU1b72egD4SwbPMJ5rt59orwda833AiqHuy1ttuvryEfVR87i5qtZW1dqJiYm5XJIkaZa6AyTJW5P8k8P7wMXAI8B24PBKqo3AXW1/O3BVW421Dnip3eraCVycZEl7eH4xsLMdeznJurb66qqhc0mSxmwut7CWAn/ZVtYuBv5nVf1VkgeAO5JsAp4GPtba7wAuA6aAXwKfAKiqQ0m+BDzQ2n2xqg61/U8BtwKnAz9omyRpAegOkKp6AnjfiPrzwEUj6gVcfZRzbQW2jqhPAuf2zlGSNH/8JrokqYsBIknqYoBIkroYIJKkLgaIJKmLASJJ6mKASJK6GCCSpC4GiCSpiwEiSepigEiSuhggkqQuBogkqYsBIknqMte/0lZ6w1m55ftjGfep6z88lnGl+eInEElSFwNEktTFAJEkdTFAJEldFnyAJFmf5GdJppJsGfd8JEkDCzpAkiwCvgZcCqwGrkyyeryzkiTBwl/Gez4wVVVPACS5HdgAPDrWWUkdxrV8GFxCrPmx0ANkGfDM0Pu9wAVjmov0huV3XzQfFnqAzEqSzcDm9vYXSX42hmmcBfz9GMYdJ6/51NB9zfnKcZ7JiXMq/nv+3WPtsNADZB+wYuj98lZ7naq6Gbj5RE1qlCSTVbV2nHM40bzmU4PXfGpIMnmsfRb0Q3TgAWBVknOSnAZcAWwf85wkSSzwTyBV9WqSa4CdwCJga1XtGfO0JEks8AABqKodwI5xz2MWxnoLbUy85lOD13xqOOZrTlXNx0QkSSe5hf4MRJK0QBkgx8Gp9nMrSVYkuSfJo0n2JPn0uOd0IiRZlOQnSb437rmcKEnOSHJnkp8meSzJ7417TvMpyX9q/00/kuRbSd4y7jnNhyRbkxxI8shQ7cwku5I83l6XzHQeA2SOTtGfW3kV+GxVrQbWAVefAtcM8GngsXFP4gT778BfVdV7gPdxEl9/kmXAfwTWVtW5DBbuXDHeWc2bW4H1R9S2AHdX1Srg7vZ+WgbI3P3/n1upql8Dh39u5aRVVfur6sdt/+cM/lBZNt5Zza8ky4EPA38+7rmcKEneAfwr4BaAqvp1Vb043lnNu8XA6UkWA78N/O8xz2deVNWPgENHlDcA29r+NuDymc5jgMzdqJ9bOan/MB2WZCVwHnDfeGcy7/4b8F+A34x7IifQOcBB4C/arbs/T/LWcU9qvlTVPuC/An8H7Adeqqq/Hu+sTqilVbW/7T8LLJ2pgwGibkneBnwH+ExVvTzu+cyXJP8WOFBVD457LifYYuD9wE1VdR7wf5jFbY03qnbPfwOD4PxnwFuT/OF4ZzUeNVieO+MSXQNk7mb1cysnmyRvYhAe36yq7457PvPsg8C/S/IUg1uUFyb5H+Od0gmxF9hbVYc/Xd7JIFBOVv8GeLKqDlbV/wW+C/yLMc/pRHouydkA7fXATB0MkLk75X5uJUkY3Bd/rKq+Ou75zLeq+nxVLa+qlQz+/f6wqk76/zOtqmeBZ5Ic/pG9izi5/yqFvwPWJfnt9t/4RZzEiwZG2A5sbPsbgbtm6rDgv4m+0J2iP7fyQeDjwMNJdrfaF9qvBujk8h+Ab7b/OXoC+MSY5zNvquq+JHcCP2aw0vAnnKTfSE/yLeBDwFlJ9gLXAtcDdyTZBDwNfGzG8/hNdElSD29hSZK6GCCSpC4GiCSpiwEiSepigEiSuhggkqQuBogkqYsBIknq8v8ALD7RcNreGysAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "parent = stats.poisson.rvs(mu=1, size=int(1e6))\n",
    "\n",
    "# show that the Poisson distribution yields integers\n",
    "print(parent[:10])\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.hist(parent, bins=np.arange(parent.max() + 2), align=\"left\")\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sanity check: if we take two independent samples from our parent distribution, we should call a difference about 5% of the time at p=0.05."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0506"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_iterations = 5000\n",
    "n_samples = 1000\n",
    "n_positive = 0\n",
    "p = 0.05\n",
    "\n",
    "for _ in range(n_iterations):\n",
    "    # Sampling with replacement is the default, but let's be explicit.\n",
    "    a = np.random.choice(parent, n_samples, replace=True)\n",
    "    b = np.random.choice(parent, n_samples, replace=True)\n",
    "    _u, test_p = stats.mannwhitneyu(a, b)\n",
    "    # multiply by 2 for a 2-sided alternative hypothesis\n",
    "    if test_p*2 < p:\n",
    "        n_positive += 1\n",
    "\n",
    "n_positive / n_iterations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's close enough for government work.\n",
    "\n",
    "## Power analysis\n",
    "\n",
    "Our empirical power analysis will tell us how often we call a statistically significant difference between two samples, which we'll call $a$ and $b$, each of size $N$, when there is a true difference $\\frac{\\textrm{mean}(\\textrm{parent}_b)}{\\textrm{mean}(\\textrm{parent}_a)} = 1+\\Delta$ in the means of the respective parent distributions of the samples.\n",
    "\n",
    "Customarily, we want to choose $N$ so that we call a difference at least $(1-\\beta) = 80\\%$ of the time.\n",
    "\n",
    "To do this, we will sweep over a range of values for $N$. At each $N$, we will repeat many times the process of drawing samples $a$ and $b$ and performing a U test, counting the number of times our test reports $p < \\alpha$. The fraction of iterations that yield a significant value will be our power.\n",
    "\n",
    "---\n",
    "\n",
    "Concretely, we'll let `a` be a sample from `parent` at each iteration. How do we generate a sample `b`?\n",
    "\n",
    "There are at least a couple of options; they both seem reasonable to me:\n",
    "\n",
    "* If the parent distribution is adequately modeled by a parametric distribution, choose `b` from a parametric distribution with parameters reflecting your model of the intervention.\n",
    "* Take `b` as an independent sample from `parent` and manipulate it under a model of the effect you expect (or fear) your intervention will have.\n",
    "\n",
    "I'll pursue the latter here since it's a little more interesting to write.\n",
    "\n",
    "I'll use a model suggesting that the effect of a bad frecency guess is generally small and distributed over the sample.\n",
    "\n",
    "A simple model is to choose $\\Delta \\cdot \\textrm{mean(parent)} \\cdot N$ values from the sample with replacement and increment them by 1. This will increase the mean of the sample by a factor of about $1+\\Delta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "parent_mean = parent.mean()\n",
    "delta = 0.05\n",
    "\n",
    "def test_sample(n_samples):\n",
    "    a = np.random.choice(parent, n_samples, replace=True)\n",
    "    b = np.random.choice(parent, n_samples, replace=True)\n",
    "    n_to_frob = int(np.round(parent_mean*delta*n_samples))\n",
    "    to_frob = np.random.randint(low=0, high=n_samples, size=n_to_frob)\n",
    "    for i in to_frob:\n",
    "        b[i] += 1\n",
    "    return {\n",
    "        \"a_mean\": a.mean(),\n",
    "        \"b_mean\": b.mean(),\n",
    "        \"p\": stats.mannwhitneyu(a, b)[1]\n",
    "    }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try the function out on a small sample and run it for a bunch of iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>a_mean</th>\n",
       "      <th>b_mean</th>\n",
       "      <th>p</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.023</td>\n",
       "      <td>1.109</td>\n",
       "      <td>0.034640</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.004</td>\n",
       "      <td>1.010</td>\n",
       "      <td>0.382725</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.036</td>\n",
       "      <td>1.043</td>\n",
       "      <td>0.246348</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.983</td>\n",
       "      <td>1.026</td>\n",
       "      <td>0.240245</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.983</td>\n",
       "      <td>1.054</td>\n",
       "      <td>0.083785</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   a_mean  b_mean         p\n",
       "0   1.023   1.109  0.034640\n",
       "1   1.004   1.010  0.382725\n",
       "2   1.036   1.043  0.246348\n",
       "3   0.983   1.026  0.240245\n",
       "4   0.983   1.054  0.083785"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "0.282"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_samples = 1000\n",
    "n_iterations = 1000\n",
    "alpha = 0.05\n",
    "\n",
    "rows = [test_sample(n_samples) for _ in range(n_iterations)]\n",
    "results = pd.DataFrame(rows)\n",
    "display(results.head())\n",
    "\n",
    "power = (results[\"p\"] < alpha).sum() / len(results)\n",
    "power"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like our observed power's in the neighborhood of 30%, which is less than the 80% we'd like to obtain.\n",
    "\n",
    "Let's try sweeping over a range of values for N and plotting the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_grid = [int(i) for i in [1e3, 3e3, 1e4, 3e4]]\n",
    "partial_results = []\n",
    "\n",
    "for n_samples in N_grid:\n",
    "    rows = [test_sample(n_samples) for _ in range(n_iterations)]\n",
    "    partial_results.append(pd.DataFrame(rows).assign(N=n_samples))\n",
    "results = pd.concat(partial_results)\n",
    "results[\"significant\"] = results[\"p\"] < alpha\n",
    "\n",
    "power = (\n",
    "    results\n",
    "    .groupby(\"N\")\n",
    "    [\"significant\"]\n",
    "    .sum()\n",
    "    / n_iterations\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'Power for delta=0.05, p=0.05')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEOCAYAAABmVAtTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmYFOW1x/HvcViGHWWVfRsQVLjoAC5xX3Ah6k001y1G4w3ojdGYQBSjhrjEJMSoSbwmxAhqTFwRiaCoKOKNooDINqyyybCD7LNz7h9V0zTDTE/PMN09y+/zPPN0d3V11amGrlPve6reMndHREQE4KhUByAiItWHkoKIiEQoKYiISISSgoiIRCgpiIhIhJKCiIhEKCmIiEiEkoKIiEQoKYiISISSgoiIRNRLdQAV1bp1a+/WrVuqwxARqVHmzp27zd3blDdfjUsK3bp1Y86cOakOQ0SkRjGztfHMp+4jERGJUFIQEZEIJQUREYlQUhARkQglBRERiVBSEBGRCCUFERGJUFIQEZEIJQUREYlQUhARkQglBRERiVBSEBGRCCUFERGJSFhSMLNnzGyLmS0q430zsz+Y2UozW2BmJyUqFhERiU8ih86eAPwJeK6M9y8GMsK/IcBT4aOIiESZNC+bsdOWsWFnDh1aNmLU0D5cMbBjQtaVsJaCu88EdsSY5XLgOQ/MAlqa2bGJikdEpCaaNC+b0RMXkr0zBweyd+YweuJCJs3LTsj6UllT6Ah8FfV6fThNRKTOKzrgbN6dy8NTl5BTUHTIezkFRYydtiwh660Rd14zs+HAcIAuXbqkOBoRkSOzJ7eAzbtz2bw7j027ctm0O5fNu3PZtCuXzXvy2Lwrl6178yg64GUuY8POnITElsqkkA10jnrdKZx2GHcfB4wDyMzMLPtbEhFJoYKiA2zdkxfs5HeFO/rdeVE7/GD6vvyiwz7bPL0e7Vuk0655OhltW9O+eTrtWqTz2LvL2bEv/7D5O7RslJBtSGVSmAzcZmYvEhSYd7n7xhTGIyJSKndnV05BcGQf7vCLj+6DHX8um3blsX1fHl7isLV+mtG2WTrtW6TTt31zzurdhvbN0yMJIPhrSOMGpe+OmzWsx+iJCw/pQmpUP41RQ/skZFsTlhTM7J/A2UBrM1sP/AKoD+DufwamApcAK4H9wE2JikVEpCx5hUVsKT6aLz6iL+7aidrx5xYcOOyzRzeuT7twB39Chxa0bZ4e7vAbRnb4xzRuwFFHWaXjKz7LKFlnH5mXTGvVXGZmps+ZMyfVYYhINefu7NiXH7Vjz4vs8DeFrzfvzi21a6ZhvaOCnX3YhdO++cGdfPsWwfQ2zRqSXj8tBVtWOWY2190zy5uvRhSaRUSi5RYUHVKg3Rx230Qf7W/dk0d+0aFH92bQqklD2rdoSIcW6Qzs0jI4sm+eTtvmDSM7/BaN6mNW+aP7mkxJQUSqjaIDzva9eZGum0MLtgfP0NmdW3jYZxs3SAuO7JunM7j7McFOPupov13zdNo2a0j9NI3uE4uSgogkxd68QjbtymVL8dF8ZId/sO9+y57DT8M8yqBNs2AH361VE07p0epgV07Yf9+2eTrNGtars0f3VUlJQURKFe/QCoVFB9i6N+/wAu0h3Tt57M07/Oi+WXq9yNF9z56tad+iYdiVkx45Q6d104akHUGhVipGSUFEDlM8tELxaZDZO3MY9ep83lm8iWOaNjik/37b3sNPw6x3lEVOtezdrhlnZLSJ9NcXT2/fIr3M0zAldfQvIiKHGTtt2WFDKxQUOVMXbaJl4/qRnXu/Y5vTrnnD8Aydg+fdt2pyZKdhSuooKYjIIdyd7DKGUDDgi/svTG5AklQqw4tIxMZdOXxv/Owy30/U0ApSfSgpiAjuzitzvuLCx2Yye/UOvn1SRxrVP3T3kMihFaT6UPeRSB23eXcuoycu5P2lWxjc/RjGXtmfrq2acEZGm6QNrSDVh5KCSB3l7rw+L5sxkxeTX3SAX3yzH987tVukQHzFwI5KAnVQ3EnBzLoDA4Esd1+auJBEJNG27MnlnomLeG/JZk7uejS/u2oA3Vs3SXVYUg2UWVMws0lRzy8H3ge+CbxhZjcmPjQRqWruzhtfZHPhYzOZuWIr917al5dHnKqEIBGxWgpdo57fBZzr7qvNrDUwHZiQyMBEpGpt25vHva8v4u3Fm/iPzi353VUD6NW2aarDkmomVlKIvkaxnruvBnD3bWZ2+MDiIlJtTVmwkfveWMTe3ELuvvg4fnBGDw0dIaWKlRQGmNlugutVGprZse6+0cwaADVnEHGROmzHvnzue2MRUxZspH+nFjx61QAy2jVLdVhSjZWZFNy9rB1/Y2BEYsIRkary9qJN3DtpIbtyChg1tA8jzuxBPQ0bLeUo9+wjM2sHFJ+Xlu3um4FPEhqViFTa1/vyGfOvxbzxxQZO6Nicv//3EI5r3zzVYUkNUWZSMLOBwFNACyA7nNzJzHYCt7r7vCTEJyIV8G7WZu55fSFf78vnJxf05taze+qmMlIhsVoK44ER7v5p9EQzO4XgzKMBCYxLRCpg1/4CfvmvxUycl03fY5sz4aZBHN+hRarDkhooVlJoUjIhALj7LDPTSc0i1cT7SzczeuJCtu3N5/bzMrjtnF40qKfWgVROrKTwlplNAZ4DvgqndQZuAN5OdGAiEtvu3AIe/FcWr8xdT592zXj6hkGc2EmtAzkysc4+ut3MLgYuJ6rQDDzp7lOTEZyIlO7D5Vu5+7UFbN6dyw/P6cnt52XQsJ7OFJcjF/PsI3d/C3grSbGISDn25Bbwq6lL+OdnX9GrbVNe/5/TGdC5ZarDklqkUqOkmtlwdx9X1cGISNn+b8U27nptARt35TDirB7ceX5v0uurdSBVq7JDZ+v6eJEk2ZdXyCNvLeHvs9bRo3UTXrnlNE7uenSqw5JaqlJJwd3/UtWBiMjhPvlyO6NenU/2zhx+cEZ3fnphH7UOJKFiJgUzGwpcwaGF5jfcXWcfiSTQ/vxCfvPWUp79ZC3dWjXmlRGnktntmFSHJXVArCuaHwd6E5ySuj6c3Am43cwudvc7khCfSJ3z2eodjHxlPut27Oem07vxs6HH0aiBWgeSHLFaCpe4e++SE83sJWA5oKQgUoVy8osYO20Z4z9eTeejG/Pi8FM4pUerVIcldUyspJBrZoPcfXaJ6YOA3ATGJFLnzF27g5GvLGD1tn3ccGpX7rroOJo01C3UJfli/a+7EXjKzJpxsPuoM7ArfK9cZnYR8ATB/Reedvdfl3i/C/As0DKc525dGCd1SW5BEb9/dzl//WgVHVs24h8/GMJpPVunOiypw2Jd0fw5MMTM2nPo0Nmb4lmwmaUBTwIXECSV2WY22d2zoma7F3jZ3Z8ys37AVKBbxTdDpOaZt+5rfvrKfFZt3cd1Q7ow+pK+NFXrQFKs3P+B7r7JzPa7+24zq8ig7IOBle6+CsDMXiQYMiM6KThQvMwWwIYKLF+kRsotKOLx91YwbuaXtG+ezvM3D+aMjDapDksEiP86hRnASVGP8ejIwYH0IGgtDCkxzxjgHTP7EdAEOD/OZYvUSAvW7+SnL89nxZa9XD2oM/dc2pfm6fVTHZZIREXbqlV9JfM1wAR3f9TMTgWeN7MT3P3AISs1Gw4MB+jSpUsVhyCSeHmFRfxx+kqe+vBL2jRtyPibBnFOn7apDkvkMInswMwmKEwX68TBO7gVuxm4CMDdPzGzdKA1sCV6pnCcpXEAmZmZnqiARRJhUfYuRr4yn6Wb9nDlyZ24b1g/WjRS60Cqp0QmhdlAhpl1J0gGVwPXlphnHXAeMMHM+gLpwNYExiSSNPmFB3jyg5U8+cFKjmnSgL99L5Pz+rZLdVgiMVU0KcR9lO7uhWZ2GzCN4HTTZ9x9sZk9AMxx98nAT4G/mtmd4bJvdHe1BKTGy9qwm5GvzCdr427+c2BHfvHNfrRs3CDVYYmUK96kYCUe4xJeczC1xLT7o55nAadXZJki1VlB0QH+PONL/vD+Clo0asC4757Mhce3T3VYInGLNyn8V4lHESlh2aY9jHxlPguzd3HZgA788rLjObqJWgdSs8SbFNaYWX/AzayBu+cnMiiRmqSw6AB/mbmKJ95bQbP0ejx13UlcfOKxqQ5LpFLKTQpmdinwZ+BLgu6j7mY2IrxVp0idtnLLHn768nzmr9/FJSe258HLT6BV04apDkuk0uJpKTwKnOPuKwHMrCcwBd27WeqwogPO0x+t4tF3l9OkQRp/unYgw/p3SHVYIkcsnqSwpzghhFYBexIUj0i19+XWvYx6ZT6fr9vJ0OPb8dAVJ9KmmVoHUjvEkxTmmNlU4GWC00avIhjc7lsA7j4xgfGJVBtFB5zx/17N2GnLSK+fxhNX/weXDeiAmW5ZLrVHPEkhHdgMnBW+3go0Ar5JkCSUFKTWW7NtH6Nenc/sNV9zft+2/Oo/T6Rt8/RUhyVS5eIZJfWmZAQiUh0dOOA898kafv32UuqnHcWjVw3gWyd1VOtAai0N3i5ShnXb9zPq1fl8unoH5/RpwyPf6k/7FmodSO2mpCBSwoEDzgufruWRt5aSZsZvr+zPVSd3UutA6gQlBZEoX+3Yz12vLeDjL7dzRkZrfvPt/nRo2SjVYYkkTYWTgpldDmxy908TEI9ISrg7//zsKx6eEtwY8JFvncjVgzqrdSB1TmVaCkOAE82snrtfXNUBiSTbhp053PXaAj5asY3Te7XiN9/uT6ejG6c6LJGUqHBScPd7EhGISLK5Oy/P+YqH3lxCkTsPXXEC1w3potaB1GlxJ4XwZjkDgSx3X5q4kEQSb+OuHO5+bSEfLt/KKT2OYeyVA+h8jFoHIkeV9YaZTYp6fjnwPsEFa5PN7MbEhyZS9dydV+eu58LHZvLZ6h388rLj+cd/n6KEIBKK1VLoGvX8LuBcd19tZq2B6cCERAYmUtW27M5l9MSFTF+6hUHdjmbslQPo1rpJqsMSqVZiJYXo22LWc/fVAO6+zcwOJDYskarj7rzxxQZ+MXkxuQVF3DesHzee1o20o1Q7ECkpVlIYYGa7Ce6h0NDMjnX3jWbWgOCeyyLV3tY9efz89YW8k7WZk7q05HdXDaBHm6apDkuk2iozKbh7WTv+xsCIxIQjUjXcnTcXbOT+NxaxL7+Iey45jpu/0UOtA5FyxHPntXZAx/BltrtvBj5JaFQiR2D73jzue2MRUxduYkDnljx6VX96tW2W6rBEaoQyk4KZDQSeAloA2eHkTma2E7jV3eclIT6RCnlr4UbunbSIPbmF/OyiPgw/owf10so8yU5ESojVUhgPjCg5nIWZnUJw5tGABMYlUiE79uVz/xuLeHPBRk7s2IJHvzOA3u3UOhCpqFhJoUlp4xu5+ywz03l8Um1MW7yJn7++kF05BYy8sDcjzupJfbUORColVlJ4y8ymAM8BX4XTOgM3AG8nOjCR8uzcn8+YyYuZ9MUG+h3bnOdvHkLfY5unOiyRGi3W2Ue3m9nFwOVEFZqBJ919ajKCEynL9CWbuXviQr7el8+Pz8/gh+f0UutApArEPPvI3d8C3kpSLCLl2pVTwAP/yuK1z9dzXPtmjL9xECd0bJHqsERqjUrdZMfMhrv7uKoORiSWD5Zt4e7XFrBtbz4/OrcXPzo3gwb11DoQqUqVvfOargCSpNmdW8DDby7hpTlfkdG2KX+9IZP+nVqmOiyRWqlSScHd/1LVgYiU5qMVW7nr1QVs2p3LrWf35MfnZ9CwnkZZEUmUmEnBzIYCV3BoofkNd9fZR5JQe/MKeXjKEv752Tp6tmnCa7eexsAuR6c6LJFaL9YVzY8DvQlOSV0fTu4E3G5mF7v7HeUt3MwuAp4gGEDvaXf/dSnzfAcYQzAq63x3v7aiGyG1y79XbuNnry5gw64cRpzZgzsv6E16fbUORJIhVkvhEnfvXXKimb0ELAdiJgUzSwOeBC4gSCqzzWyyu2dFzZMBjAZOd/evzaxtJbZBaol9eYX8+q2lPD9rLd1bN+HVW07l5K7HpDoskTolVlLINbNB7j67xPRBQG4cyx4MrHT3VQBm9iLBNQ9ZUfP8gOC6h68B3H1L3JFLrTJr1XZGvTqf9V/ncPM3ujPywj40aqDWgUiyxUoKNwJPmVkzDnYfdQZ2he+VpyMHr4QmXMaQEvP0BjCzfxN0MY0prV5hZsOB4QBdunSJY9VSU+zPL+S3by9jwsdr6NqqMS8NP5XB3dU6EEmVWFc0fw4MMbP2HDp09qYqXn8GcDZBvWKmmZ3o7jtLxDIOGAeQmZnpJRciNcekedmMnbaMDTtzaNW0AbizbV8BN57WjZ9d1IfGDSp7lrSIVIVyf4HuvsnM9rv7bjOryMAy2QQti2KdODgEd7H1wKfuXgCsNrPlBEmiZJeV1AKT5mUzeuJCcgqKANi2Nx8DfnhOT0YNPS61wYkIAPFeDjqjxGM8ZgMZZtY9vIXn1cDkEvNMImglYGatCbqTVlVgHVKDjJ22LJIQijkwad6G1AQkIoep6BgBcV/J7O6FwG3ANGAJ8LK7LzazB8zssnC2acB2M8sCPgBGufv2CsYkNcDKLXvJ3plT6nsbypguIsmX0A7ccDTVqSWm3R/13IGfhH9SC+3cn8/j763g+VlrMYKWQUkdWjZKdlgiUgZV9SQhCooO8MKstTz23gr25BZwzeAu9D22GQ9PWXpIF1Kj+mmMGtonhZGKSLSKJgWd+SPl+mDZFh56M4svt+7j9F6tuG9YP45rH5yj0LRh/cjZRx1aNmLU0D5cMbBjOUsUkWSJNylYiUeRw6zcsoeHpixhxrKtdG/dhKdvyOS8vm0xO/jf5oqBHZUERKqxeJPCf5V4FIn4el8+T0wP6gaNG6Rx76V9ueHUbrrXgUgNFFdScPfl0Y8iENQN/j5rLY9H1Q1+ckFvWjVtmOrQRKSSyk0K4aB1jwD9gPTi6e7eI4FxSTUXXTf4Rq/W3Dusb6RuICI1VzwthfHAL4DHgHOAm6j49Q1SS6zYHNQNPlxedt1ARGqueJJCI3efbmbm7muBMWY2F7i/vA9K7fH1vnwef285f/90neoGIrVYPEkhz8yOAlaY2W0E4xc1TWxYUl0UFB3g+U/W8sT0oG5w3ZCu/Pj8DNUNRGqpeJLCHUBj4HbgQYIupBsSGZSknrszY9lWHpySxaqwbnDfsH70ad8s1aGJSALFkxS6hTfa2UtQT8DMrgI+TWRgkjorNu/hwSlLmBnWDf72vUzOPU51A5G6IJ6kMBp4JY5pUsPtCOsGL3y6jiYN0rhvWD++e0pX1Q1E6pAyk4KZXQxcAnQ0sz9EvdUcKEx0YJI8xXWDx99bzt68Qq4b0pU7L+jNMU0apDo0EUmyWC2FDcBc4LLwsdge4M5EBiXJ4e7B9QZTlrBq6z7OyGjNvZeqbiBSl8W6Hed8YL6Z/T28N4LUIss37+HBN7P4aMU2erRuwjM3ZnJOH9UNROq6WN1HCwlHRS1tR+Hu/RMXliTKjn35PPbucv7xmeoGInK4WN1Hw5IWhSRcfuEBnp+1lifeW86+/CKuG9KFH5+vuoGIHCpW99Ha4udm1hXIcPf3zKxRrM9J9eLuvL90Cw9PWcKqbUHd4L5h/ejdTnUDETlcPAPi/QAYDhwD9AQ6AX8GzktsaHKkVDcQkYqK54j/h8BgwovV3H2FmbVNaFRyRIrrBi98upamDetx/7B+fPfUrtRPU91ARGKLa+wjd88vPro0s3rotpzVUn7hAZ77ZA1PTF/B/vwirj+lK3ee35ujVTcQkTjFkxQ+NLN7gEZmdgHwP8C/EhuWVIS7M33JFh6euoTVqhuIyBGIJyncDdwMLARGAFOBpxMZlMRv2aY9PDQlrBu0acL4Gwdxdp82qhuISKWUmxTc/QDw1/BPqonte/N47L3l/OPTdaobiEiVievitdLo4rXUKFk3+O4pXfmx6gYiUkXiuXjth+Hj8+Hj9ajQnHQl6wZn9m7DfZf2JUN1AxGpQuVevGZmF7j7wKi37jKzzwlqDZIESzft5qE3l/B/K7fRs00Txt80iHP66KxgEal68RSazcxOd/d/hy9OA9RxnQTb9+bx+3eX88/P1tEsvT6/+GY/rj9FdQMRSZx4ksLNwDNm1iJ8vRP4fuJCkpJ1gxtO7cYd52WobiAiCRfP2UdzgQHFScHddyU8qjrK3XlvyRYenpLFmu37Oat3G+4b1pdebVU3EJHkiHtgu8okAzO7CHgCSAOedvdflzHft4FXgUHuPqei66kNlm7azYNvZvHvldtVNxCRlEnYaKdmlgY8CVwArAdmm9lkd88qMV8z4A7CsZXqmpJ1gzHf7Md1qhuISIrETApmdhRwirt/XIllDwZWuvuqcFkvApcDWSXmexD4DTCqEuuosfILD/Dsx2v4w/QV7C8I6gY/Pj+Dlo1VNxCR1ImZFNz9gJk9CQyMNV8ZOgJfRb1eDwyJnsHMTgI6u/sUM6sTScHdeTdrM7+auoQ12/dzdp823Hup6gYiUj3E0300Pezzn+juVXbRWtgK+T1wYxzzDie4pwNdunSpqhCSbsnG3Tw0Jagb9GrbVHUDEal24kkKI4CfAEVmlgMY4O7evJzPZQOdo153CqcVawacAMwIB29rD0w2s8tKFpvdfRwwDiAzM7PGXU29LawbvPjZOpo3qs8vLzuea4d0Ud1ARKqdeE5JrWy/xmwgw8y6EySDq4Fro5a7C2hd/NrMZgAja9PZR3mFRTz78Rr+OH2l6gYiUiPEdfaRmV0GnBm+nOHub5b3GXcvNLPbgGkEp6Q+4+6LzewBYI67T65s0NVdcd3g4alLWLt9P+f0acPPVTcQkRognns0/xoYBLwQTrojHPZidHmfdfepBPdfiJ52fxnznl1utDXAko3B9QYffxnUDSbcNIizVTcQkRoinpbCJcB/hPdVwMyeBeYB5SaFumTb3jwefWc5L81W3UBEaq54L15rCewIn7eINWNdk1dYxIR/r+FP768kp6CI750WjFOkuoGI1ETxJIVHgHlm9gHBmUdnomGzcXfeCa83WLt9P+ce15Z7LulLr7ZNUx2aiEilxbrzWvFw2ROBGQR1BYC73H1TEmKrtrI2BHWDT1YFdYNnvz+Ys3q3SXVYIiJHLFZL4Q/AycAn7n4SUGvPFopXUDdYxouzv6JFo/o8cPnxXDu4C/VUNxCRWiJWUigws3FAJzP7Q8k33f32xIVVvRTXDf74/kpyC4q46bTu3HFeBi0a1091aCIiVaq8ezSfDwwF5iYnnOrF3Zm2OKgbrNsR1A1+fmlferZR3UBEaqdY92jeBrxoZkvcfX4SY6oWFm/YxYNvZjFr1Q4yVDcQkToinmEu6lRC2Lonj9+/G9QNWjaqz4OXH881qhuISB2RsJvs1DR5hUWMD683UN1AROqqeG6yc6W7v5ykeJIuqBts4ldTl7Jux37OO64t96huICJ1VDw32fkZUCuTwqLsoG7w6eod9G7XlOe+P5gzVTcQkTosnu6j98xsJPASsK94orvvKPsj1cukedmMnbaMDTtz6NCyESPO6kHWht28NCesG1xxAtcM6qy6gYjUeVbezdTMbHUpk93deyQmpNgyMzN9zpz4b7kwaV42oycuJKeg6JDpBnz/G925/VzVDUSk9jOzue6eWd588Zx91L1qQkqNsdOWHZYQANo0a8h9w/qlICIRkeqr3P4SM2tsZveGVzdjZhlmNizxoVWNDTtzSp2+dU9ekiMREan+4ulEHw/kA6eFr7OBhxIWURXr0LJRhaaLiNRl8SSFnu7+W6AAwN33E3TJ1wijhvahUf20Q6Y1qp/GqKF9UhSRiEj1Fc/ZR/lm1ghwADPrCdSYvpcrBnYEOOTso1FD+0Smi4jIQfEkhTHA20BnM3sBOB24MYExVbkrBnZUEhARiUM8Zx+9Y2ZzgVMIuo3uCAfLExGRWqbcpGBmfwc+BD5y96WJD0lERFIlnkLz34BjgT+a2Soze83M7khwXCIikgLxdB99YGYzCe7RfA5wC3A88ESCYxMRkSSLp/toOtAE+AT4CBjk7lsSHZiIiCRfPN1HCwguXjsB6A+cEJ6iKiIitUw83Ud3AphZM4JTUccD7YGGCY1MRESSLp7uo9uAM4CTgTXAMwTdSCIiUsvEc/FaOvB7YK67FyY4HhERSaF4uo9+Z2YDgFvMDILrFeYnPDIREUm6eIbOvh14AWgb/v3dzH6U6MBERCT54jn76L+BIe5+v7vfTzDcxQ/iWbiZXWRmy8xspZndXcr7PzGzLDNbYGbTzaxrxcIXEZGqFE9SMCD61mVFxDF0tpmlAU8CFwP9gGvMrOStzuYBme7eH3gV+G08QYuISGLEU2geD3xqZq+Hr68gGPqiPIOBle6+CsDMXgQuB7KKZ3D3D6LmnwVcH0/QIiKSGPEUmn9vZjOAb4STbnL3eXEsuyPwVdTr9cCQGPPfDLxV2htmNhwYDtClS5c4Vi0iIpVRZlIws3SCcY56AQuB/03UKalmdj2QCZxV2vvuPg4YB5CZmemJiEFERGK3FJ4luAXnRwR1gb7Ajyuw7Gygc9TrTuG0Q5jZ+cDPgbPcvcbc0U1EpDaKlRT6ufuJAGb2N+CzCi57NpBhZt0JksHVwLXRM5jZQOAvwEUaZE9EJPVinX1UUPykMt1G4WduA6YBS4CX3X2xmT1gZpeFs40FmgKvmNkXZja5ousREZGqE6ulMMDMdofPDWgUvjbA3b15eQt396nA1BLT7o96fn7FQxYRkUQpMym4e1oyAxERkdSL5+I1ERGpI5QUREQkQklBREQilBRERCRCSUFERCKUFEREJEJJQUREIpQUREQkQklBREQilBRERCRCSUFERCKUFEREJEJJQUREIpQUREQkQklBREQilBRERCRCSUFERCKUFEREJEJJQUREIpQUREQkQklBREQilBRERCRCSUFERCKUFEREJEJJQUREIpQUREQkQklBREQilBRERCRCSUFERCISmhTM7CIzW2ZmK83s7lLeb2hmL4Xvf2pm3RIZj4iIxJawpGBmacCTwMVAP+AaM+tXYrabga/dvRfwGPCbRMUjIiLlS2RLYTCw0t1XuXs+8CJmK7LDAAAIPUlEQVRweYl5LgeeDZ+/CpxnZpbAmEREJIZEJoWOwFdRr9eH00qdx90LgV1AqwTGJCIiMdRLdQDxMLPhwPDw5V4zW1bKbC0IkkosrYFtVRlbDRHPd5NsyYipqtdRFcurzDIq+pmKzK/fTdlq2++ma1xzuXtC/oBTgWlRr0cDo0vMMw04NXxej+A/nlVyfePimGdOora3Ov/F893Uxpiqeh1VsbzKLKOin6nI/PrdJPbfuybGlMjuo9lAhpl1N7MGwNXA5BLzTAa+Fz6/Enjfwy2vhH9V8nN1QXX8bpIRU1WvoyqWV5llVPQzFZm/Ov7fqC6q43eT8Jis8vvgOBZudgnwOJAGPOPuD5vZAwRHHpPNLB14HhgI7ACudvdVCYxnjrtnJmr5IrWRfjd1S0KTQnVjZsPdfVyq4xCpSfS7qVvqVFIQEZHYNMyFiIhEKCmIiEiEkoKIiETU6aRgZn3N7M9m9qqZ3ZrqeERqCjNrYmZzzGxYqmORqlXrkoKZPWNmW8xsUYnph43Y6u5L3P0W4DvA6amIV6Q6qMjvJnQX8HJyo5RkqHVJAZgAXBQ9IdaIrWZ2GTAFmJrcMEWqlQnE+bsxswuALGBLsoOUxKsRYx9VhLvPLOW+DJERWwHMrHjE1ix3nwxMNrMpwD+SGatIdVHB301ToAlBosgxs6nufiCJ4UoC1bqkUIbSRmwdYmZnA98CGqKWgkhJpf5u3P02ADO7EdimhFC71JWkUCp3nwHMSHEYIjWSu09IdQxS9WpjTaE02UDnqNedwmkiUjb9buqgupIU4hmxVUQOpd9NHVTrkoKZ/RP4BOhjZuvN7GYP7up2G8H9G5YAL7v74lTGKVKd6HcjxTQgnoiIRNS6loKIiFSekoKIiEQoKYiISISSgoiIRCgpiIhIhJKCiIhEKClImczMzezRqNcjzWxMFS17gpldWRXLKmc9V5nZEjP7INHrKieONWbW+giXcYuZ3VBVMVVgvWPMbGQc8x3xNkrqKSlILHnAt6rbD93MKjJm183AD9z9nETFkyzu/md3fy7VcUjtpqQgsRQC44A7S75R8kjfzPaGj2eb2Ydm9oaZrTKzX5vZdWb2mZktNLOeUYs5P7x71/LiO3iZWZqZjTWz2Wa2wMxGRC33IzObTDCWf8l4rgmXv8jMfhNOux/4BvA3MxtbYv5jzWymmX0RfuaMcPpTYUyLzeyXUfOvMbNHwvnnmNlJZjbNzL40s1uiYpxpZlPCG9P82cwO+42Z2fXh9/GFmf0l3Oa08DtdFG5Had955IjdzGaY2W/C5Swvjr+Uz4yK+i6jt2eSmc0Nt3N41PSLzOxzM5tvZtOjFtUvXOcqM7u9tHWVt43h9JvDeD8zs7+a2Z/KW5YkmbvrT3+l/gF7gebAGqAFMBIYE743Abgyet7w8WxgJ3AswZDk2cAvw/fuAB6P+vzbBAcmGQTDMqcDw4F7w3kaAnOA7uFy9wHdS4mzA7AOaEMw8u/7wBXhezOAzFI+81Pg5+HzNKBZ+PyYqGkzgP7h6zXAreHzx4AFQLNwnZujtj0X6BF+/t3i7yj8fGugL/AvoH44/X+BG4CTgXej4mtZSsxjgJFR2/Vo+PwS4L1S5r+QIKlb+D2/CZxZYjsbAYuAVuG2fFX8HUfNMwb4OPz3aA1sL46/xPrK28YO4TzHAPWBj4A/pfr/uf4O/avTQ2dL+dx9t5k9B9wO5MT5sdnuvhHAzL4E3gmnLwSiu3Fe9mAs/hVmtgo4jmBH1j+qFdKCIGnkA5+5++pS1jcImOHuW8N1vgCcCUyKFSPwjJnVBya5+xfh9O+ER871CBJbP4IEAAcHg1sINHX3PcAeM8szs5bhe5/5wZvS/JOgpfJq1HrPI0gAs80Mgp3yFoKdaA8z+yPBnQDfoXwTw8e5QLdS3r8w/JsXvm5K8F3OBG43s/8Mp3cOp7cBZhZ/x+6+I2pZU9w9D8gzsy1AO4JEXpqytnEw8GHxcs3sFaB3HNspSaSkIPF4HPgcGB81rZCw+zHsImkQ9V5e1PMDUa8PcOj/uZIDbznBUe2P3H1a9BsW3BBpX+XCP5wHdxo7E7gUmGBmvyc4ch0JDHL3r81sAkHrpVj0dpTcxuLtKm2bohnwrLuPLhmTmQ0AhgLF9w3/fjmbURxDEaX/lg14xN3/UmI9ZwPnA6e6+34zm8Gh2xlrXbHWF73ew7bRzK4oZx1SDaimIOUKj+xeJijaFltDcDQIcBlBd0BFXWVmR4V1hh7AMoIROW8Nj+Axs95m1qSc5XwGnGVmrcO+62uAD2N9wMy6EnT7/BV4GjiJoKtsH7DLzNoR3Ju4ogZbMNT0UcB/Af9X4v3pwJVm1jaM4xgz62pBMf8od38NuDeM50hNA75vZk3DdXUM19sC+DpMCMcBp4TzzwLONLPuxbFVcr2lbiNB6+wsMzvagpMFvl3pLZOEUUtB4vUowTDKxf4KvGFm8wlqA5U5il9HsENvDtzi7rlm9jRBV8jnFvQ9bAViHmG6+0Yzuxv4gOAodYq7v1HOus8GRplZAUHt5AZ3X21m84ClBH3r/67ENs0G/gT0CuN5vUSsWWZ2L/BOmDgKgB8SdM2NjypMH9aSqCh3f8fM+gKfhN04e4HrCf69bjGzJQSJeFY4/9aw62xiGMcW4IJKrLfUbXT3WWb2K4J/8x0E3/OuI91OqVoaOlukioTdMiPdfViqY6muzKypu+8NWwqvA8+4++vlfU6SR91HIpJMY8zsC4IznlYT+2QASQG1FEREJEItBRERiVBSEBGRCCUFERGJUFIQEZEIJQUREYlQUhARkYj/B0LwJ6vheHWrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.semilogx(power.index, power.values, \"o-\")\n",
    "ax.set_ylim(0, 1.1)\n",
    "ax.set_xlabel(\"Number of samples in each leg\")\n",
    "ax.set_ylabel(f\"Power for delta={delta}, p={p}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this model, around 10,000 samples in each leg looks like enough to detect a one-sided difference of 5% in the mean with high confidence.\n",
    "\n",
    "But we should make sure we're running enough iterations at each N to get a reliable answer. To build confidence, you can repeat the analysis for a given number of iterations at each N many times, and plot the average and some measure of spread of the observed powers in order to convince yourself that you're getting reliable results."
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

requirements.txt

appnope==0.1.0
backcall==0.1.0
bleach==2.1.3
cycler==0.10.0
decorator==4.3.0
entrypoints==0.2.3
html5lib==1.0.1
ipykernel==4.8.2
ipython==6.4.0
ipython-genutils==0.2.0
ipywidgets==7.3.0
jedi==0.12.1
Jinja2==2.10
jsonschema==2.6.0
jupyter==1.0.0
jupyter-client==5.2.3
jupyter-console==5.2.0
jupyter-core==4.4.0
kiwisolver==1.0.1
MarkupSafe==1.0
matplotlib==2.2.2
mistune==0.8.3
nbconvert==5.3.1
nbformat==4.4.0
notebook==5.6.0
numpy==1.14.5
pandas==0.23.3
pandocfilters==1.4.2
parso==0.3.1
pexpect==4.6.0
pickleshare==0.7.4
prometheus-client==0.3.0
prompt-toolkit==1.0.15
ptyprocess==0.6.0
Pygments==2.2.0
pyparsing==2.2.0
python-dateutil==2.7.3
pytz==2018.5
pyzmq==17.1.0
qtconsole==4.3.1
scipy==1.1.0
seaborn==0.9.0
Send2Trash==1.5.0
simplegeneric==0.8.1
six==1.11.0
terminado==0.8.1
testpath==0.3.1
tornado==5.1
traitlets==4.3.2
wcwidth==0.1.7
webencodings==0.5.1
widgetsnbextension==3.3.0