tests.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hypothesis Testing, Combinations of Tests, and Stratified Tests"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hypotheses\n",
    "\n",
    "A _scientific hypothesis_ is an assertion about the world, for instance \"this drug increases longevity in patients with Alzheimer's Disease,\" \"the outcome of the 2020 presidential election would have been different, but for fraud,\" \"to 6 significant digits, the gravitational constant is $6.67430e-11 Nm^2kg^{-2}$,\" \"the Moderna COVID-19 vaccine is at least 90% effective for at least 6 months for all ages,\" \"students give female instructorss lower ratings than equally effective male instructors,\" \"the universe is expanding at an increasing rate,\" \"there is a black hole at the center of the universe,\" or \"female faculty job applicants are interrupted more frequently during their job talks than male applicants are.\"\n",
    "\n",
    "A _statistical hypothesis_ is an assertion about the probability distribution $\\mathbb{P}$ of data $X$.\n",
    "\n",
    "Generally, statistical hypotheses are interesting only insofar as they help us evaluate scientific \n",
    "hypotheses.\n",
    "Whenever a statistical hypothesis is being tested, pay attention to the connection between that hypothesis and the scientific hypothesis it purports to represent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Types of Statistical Hypotheses\n",
    "\n",
    "A _simple_ statistical hypothesis is an assertion that completely specifies the probability distribution of the data, e.g., $\\mathbb{P} = \\mathbb{P}_0$ for some distribution $\\mathbb{P}_0$. Here are some simple hypotheses:\n",
    "\n",
    "* $X \\sim N(0,1)$ The data have a standard normal distribution.\n",
    "* $X \\sim U[0,1]$ The data have a standard uniform distribution.\n",
    "* $X \\sim N(\\mu, \\Sigma)$ for a given $\\mu \\in \\Re^n$ and positive semi-definite $n \\times n$ matrix $\\Sigma$. The data are jointly normally distributed with mean $\\mu$ and covariance matrix $\\Sigma$.\n",
    "* $\\{X_j\\}\\; \\mbox{IID} \\; U[0, 17]$. The components of the data are independent and identically distributed with a uniform distribution on $[0, 17]$.\n",
    "* $X \\sim \\mbox{Binom}(n, p)$ for given values of $n$ and $p$. The data have a binomial distribution with parameters $n$ and $p$.\n",
    "\n",
    "A _composite_ hypothesis is an assertion that does not completely specify the distribution; it only says the\n",
    "distribution is in some specified set of distributions. \n",
    "E.g., $\\mathbb{P} \\in \\mathcal{P}_0$ for some set\n",
    "$\\mathcal{P}_0$ of distributions. \n",
    "Here are some composite null hypotheses:\n",
    "\n",
    "* $X \\sim N(\\mu,1)$, $\\mu \\in [0, 7]$ The data have a normal distribution with variance 1 and mean between 0 and 7. \n",
    "* $X \\sim \\mbox{Binom}(n, p)$, $p > 1/2$ ($n$ given). The data have a binomial distribution with known parameter $n$ and unknown parameter $p > 1/2$.\n",
    "* $\\mathbb{E}X=0$ The expected value of the data is zero.\n",
    "* $\\mathbb{E}X \\le 1/2$ The expected value of the data is at most 1/2.\n",
    "* The distribution $\\mathbb{P}$ of $X$ has at most 3 modes.\n",
    "* $\\{X_j\\}_{j=1}^n$ are _exchangeable_ (A collection of random variables is _exchangeable_ if their joint distribution is invariant under permuting their labels. If a collection of random variables is IID, it is also exchangeable.)\n",
    "* The distribution of $X$ is spherically symmetric\n",
    "* The $q$th quantile of the probability distribution of the data is $x$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A common statistical question is to use data to assess whether a statistical hypothesis is true: to _test_ the hypothesis.\n",
    "There is always more than one explanation for any particular set of data, so in general, it is not\n",
    "possible to prove that a (simple) statistical hypothesis is _true_.\n",
    "But the data might provide evidence that a hypothesis is _false_.\n",
    "This is tied to Popper's notion that only hypotheses that can be falsified by evidence (potentially shown to be false) are scientific hypotheses.\n",
    "\n",
    "We speak of \"rejecting\" a hypothesis when there is sufficiently strong evidence that the hypothesis is false.\n",
    "But if data do not cast doubt on a hypothesis, that is not evidence that the hypothesis is true--it is only absense of evidence that the hypothesis is false.\n",
    "Additional evidence might cast doubt on the hypothesis.\n",
    "When we speak (informally) of \"accepting\" a hypothesis, it means only that we did not reject it on the basis of a particular set of data, not that we\n",
    "have affirmative evidence that it is true."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To test a statistical hypothesis, one specifies a set $A$ of possible data values (before examining the data). \n",
    "If the data fall inside that set, i.e., if $X \\in A$, the null hypothesis is _not rejected_; otherwise, the null hypothesis is _rejected_.\n",
    "Each (measurable) set $A$ implicitly defines a hypothesis test.\n",
    "The set $A$ is called the _acceptance region_ for the test.\n",
    "\n",
    "Sometimes the set $A$ is defined explicitly, but\n",
    "more often it is defined implicitly in terms of a _test statistic_, a function of the data that does not depend on any unknown parameters.\n",
    "\n",
    "For instance an acceptance region for testing the hypothesis that $\\{X_j\\}_{j=1}^n$ are IID $N(0,1)$ might be\n",
    "\\begin{eqnarray*}\n",
    "A := \\left \\{x = \\{x_1, \\ldots, x_n\\}: \\frac{1}{n} \\sum_{j=1}^n x_j \\le \\frac{c}{\\sqrt{n}} \\right \\}.\n",
    "\\end{eqnarray*}\n",
    "In this example, the test statistic is the sample mean, and the region $A$ is the set of all data for which the\n",
    "sample mean does not exceed $c/\\sqrt{n}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test functions versus acceptance regions\n",
    "\n",
    "Instead of working with a set $A$, we can work with the indicator function of the set $A$, \n",
    "\\begin{eqnarray*}\n",
    "1_A &:& \\mathcal{X} \\rightarrow \\{0, 1\\} \\\\\n",
    "       && x  \\mapsto \\left \\{ \\begin{array}{ll}\n",
    "                    1, & x \\in A \\\\\n",
    "                    0, & \\mbox{ otherwise. } \n",
    "                    \\end{array} \\right .\n",
    "\\end{eqnarray*}\n",
    "Then we reject the null hypothesis if $1_A(X) = 0$ and do not reject it if $1_A(X) = 1$.\n",
    "\n",
    "This can be generalized by using a _test function_ or _critical function_ $\\phi(\\cdot)$ that take \n",
    "values in $[0,1]$ rather than just $\\{0, 1\\}$, as described below.\n",
    "In that case, when $X=x$, the test rejects the null with probability $\\phi(x)$, i.e.,\n",
    "if $\\phi(X) = 1$, the test certainly rejects the null;\n",
    "if $\\phi(X) = 0$, the test certainly does not reject the null; and if $\\phi(X) = q \\in (0, 1)$,\n",
    "it rejects the null with probability $q$.\n",
    "This is called a _randomized test_, discussed further below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aside on notation\n",
    "\n",
    "The following all mean the same thing, namely, the probability that the random variable $X$ takes a value in the set $A$ if the distribution of $X$ is $\\mathbb{P}_0$:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\mathbb{P}_0 \\{X \\in A\\},\n",
    "\\end{equation*}\n",
    "\\begin{equation*}\n",
    "\\Pr \\{X \\in A || X \\sim \\mathbb{P}_0 \\},\n",
    "\\end{equation*}\n",
    "and\n",
    "\\begin{equation*}\n",
    "\\Pr_{X \\sim \\mathbb{P}_0} \\{X \\in A\\}.\n",
    "\\end{equation*}\n",
    "\n",
    "\n",
    "For expectations, we use similar notation:\n",
    "\\begin{equation*}\n",
    "\\mathbb{E}_{\\mathbb{P}_0} X,\n",
    "\\end{equation*}\n",
    "\\begin{equation*}\n",
    "\\mathbb{E} (X || X \\sim \\mathbb{P}_0),\n",
    "\\end{equation*}\n",
    "\\begin{equation*}\n",
    "\\mathbb{E}_{X \\sim \\mathbb{P}_0} X,\n",
    "\\end{equation*}\n",
    "and sometimes\n",
    "\\begin{equation*}\n",
    "\\mathbb{E}_0 X\n",
    "\\end{equation*}\n",
    "or even\n",
    "\\begin{equation*}\n",
    "\\mathbb{P}_0 X\n",
    "\\end{equation*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Significance level\n",
    "\n",
    "The _significance level_ of the test $A$ of the simple null hypothesis $X \\sim \\mathbb{P}_0$ is the probability that the test rejects the null hypothesis when the null hypothesis is true:\n",
    "\\begin{equation*}\n",
    "   \\alpha := \\mathbb{P}_0 \\{X \\notin A \\}.\n",
    "\\end{equation*}\n",
    "Because the expected value of an indicator function is the probability of the set, this can be written\n",
    "\\begin{equation*}\n",
    "   \\alpha := 1-\\mathbb{E}_{\\mathbb{P}_0} 1_A(X).\n",
    "\\end{equation*}\n",
    "For a test based on the test function $\\phi: \\mathcal{X} \\rightarrow [0,1]$,\n",
    "the significance level of the simple null hypothesis $X \\sim \\mathbb{P}_0$ is\n",
    "\\begin{equation*}\n",
    "   \\alpha := \\mathbb{E}_{\\mathbb{P}_0} \\phi(X).\n",
    "\\end{equation*}\n",
    "\n",
    "Why is this expression the probability of rejecting the null when the null is true?\n",
    "The test rejects when $U \\le \\phi(X)$.\n",
    "By the law of total expectation (see [Inequalities and Identities](./math-inequalities.ipynb)) and \n",
    "using the fact that $X$ and $U$ are independent,\n",
    "\\begin{equation*}\n",
    "\\mathbb{P}_{0,U} \\{ U \\le \\phi(X)\\} = \n",
    "\\mathbb{E}_{\\mathbb{P}_0, U} 1_{U \\le \\phi(X)} = \n",
    "\\mathbb{E}_{\\mathbb{P}_0} \\mathbb{E}_U (1_{U \\le \\phi(X)} | \\phi(X)) =\n",
    "\\mathbb{E}_{\\mathbb{P}_0} \\phi(X).\n",
    "\\end{equation*}\n",
    "\n",
    "The _significance level_ of the test $A$ of the composite null hypothesis $H_0: X \\sim \\mathbb{P} \\in \\mathcal{P}_0$ is\n",
    "\\begin{equation*}\n",
    "   \\alpha := \\sup_{\\mathbb{P} \\in \\mathcal{P}_0} \\mathbb{P} \\{X \\notin A \\}.\n",
    "\\end{equation*}\n",
    "That is, it is the largest probability that the test rejects the null hypothesis when the null hypothesis is true.\n",
    "For a test based on a a test function $\\phi: \\mathcal{X} \\rightarrow [0,1]$,\n",
    "the significance level of the composite null hypothesis $X \\sim \\mathcal{P}_0$ is\n",
    "\\begin{equation*}\n",
    "   \\alpha := \\sup_{\\mathbb{P_0} \\in \\mathcal{P}_0} \\mathbb{E}_{\\mathbb{P}_0} \\phi(X).\n",
    "\\end{equation*}\n",
    "\n",
    "When $\\mathcal{P}_0$ contains only one distribution, the two definitions of significance level coincide."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Power\n",
    "\n",
    "The _power_ of the test $A$ against the simple hypothesis $H_1: X \\sim \\mathbb{P}_1$ is \n",
    "\n",
    "\\begin{equation*}\n",
    "  \\beta := \\mathbb{P}_1 \\{X \\notin A \\}\n",
    "\\end{equation*}\n",
    "That is, the power is the chance that the test rejects the null hypothesis $H_0$ when the alternative hypothesis $H_1$ is true.\n",
    "\n",
    "For a test based on a test function $\\phi$, the power is\n",
    "\\begin{equation*}\n",
    "  \\beta := \\mathbb{E}_{\\mathbb{P}_1} \\phi(X).\n",
    "\\end{equation*}\n",
    "\n",
    "The _power_ of the test $A$ against the composite alternative hypothesis $H_1: X \\sim \\mathbb{P} \\in \\mathcal{P}_1$ is\n",
    "\\begin{equation*}\n",
    "   \\beta := \\inf_{\\mathbb{P_1} \\in \\mathcal{P}_1} \\mathbb{P}_1 \\{X \\notin A \\}.\n",
    "\\end{equation*}\n",
    "That is, it is the smallest probability that that the test rejects the null hypothesis when the alternative hypothesis is true.\n",
    "\n",
    "For a test based on a test function $\\phi$, the power is\n",
    "\\begin{equation*}\n",
    "   \\beta := \\inf_{\\mathbb{P_1} \\in \\mathcal{P}_1}\\mathbb{E}_{\\mathbb{P}_1} \\phi(X).\n",
    "\\end{equation*}\n",
    "\n",
    "When the alternative $\\mathcal{P}_1$ contains only one distribution, the two definitions coincide.\n",
    "\n",
    "Power and significance level are the (extremal) probability of the same event--namely, rejecting the null hypothesis--but computed under different assumptions.\n",
    "The power is computed under the assumption that the alternative hypothesis is true; the significance level is computed under the assumption that the null hypothesis is true."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Seminal work by Jerzy Neyman (the founder of Berkeley's Department of Statistics) and Egon Pearson showed how to find the most powerful test of a simple null hypothesis against a simple alternative hypothesis among all tests with a given significance level. They showed that the most powerful test had an acceptance region characterized by the likelihood ratio; see below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Type I and Type II errors\n",
    "\n",
    "A type I error occurs when a test rejects the null hypothesis but the null hypothesis is true.\n",
    "The chance of a type I error is the significance level of the test.\n",
    "\n",
    "\n",
    "A type II error occurs when a test does not reject the null hypothesis the null hypothesis is false.\n",
    "The chance of a type II error when a particular alternative is true is 100% minus the power of the test\n",
    "against that alternative, i.e., $\\beta$.\n",
    "\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Type III errors\n",
    "\n",
    "There are a number of informal definitions of type III errors. \n",
    "\n",
    "One is that a Type III error occurs when a test correctly rejects the null hypothesis \"for the wrong reason.\" For instance, suppose that we test the null hypothesis $X \\sim N(0, 1)$\n",
    "at significance level $\\alpha = 0.05$ by defining $A := \\{x: |x| \\le 1.96\\}$. \n",
    "Suppose that in reality $X \\sim N(-1, 1)$ and we observe $X = 2$. \n",
    "Then we would _correctly_ reject $H_0$, but because $X$ was \"too big,\" while it was much more likely to be \"too small\" since in reality the distribution of $X$ has mean $-1$ instead of $0$. (If, after rejecting $H_0$, we concluded that $\\mathbb{E}X > 0$, that directional conclusion would be wrong.)\n",
    "\n",
    "Similarly, suppose we want to test whether the mean of some finite population is zero from a simple random sample $X = (X_1, \\ldots, X_n)$.\n",
    "We set $A = \\{x \\in \\Re^n : |\\bar{x}| \\le z_\\alpha/\\sqrt{n} \\}$, where $z_\\alpha$ is the $1-\\alpha$ percentage\n",
    "point of the standard normal distribution.\n",
    "This is a significance level $\\alpha$ test of the null hypothesis that $X \\sim N(0, 1)$ from an IID sample, not \n",
    "of the null hypothesis that the population mean is zero from a random sample without replacement.\n",
    "Nonetheless, this incorrect test might correctly reject the null hypothesis.\n",
    "\n",
    "Another informal definition is that a Type III error occurs when one gets the right answer to the wrong question.\n",
    "One example is testing the hypothesis $\\mathbb{E}X=0$ by testing the hypothesis $X \\sim N(0, \\sigma^2)$ when\n",
    "there is no reason to think that $X$ has a normal distribution. Even if the test is performed correctly,\n",
    "it is testing the wrong null.\n",
    "\n",
    "In my experience, the most frequent and pernicious Type III errors involve testing a _statistical hypothesis_ that has little or nothing to do with the _scientific hypothesis_ (aside, perhaps, from having some words in common).\n",
    "*Many* hypothesis tests in applications have this sort of Type III error baked in. \n",
    "Indeed, most analyses of clinical trial data I've seen test a statistical null hypothesis that involves selecting subjects at random from a superpopulation--which did not occur--rather than test a statistical hypothesis based on the randomization of subjects into treatment or control--which actually did occur.\n",
    "\n",
    "It is quite common to base $\\mathbb{P}_0$ on a statistical model that has no connection to the scientific hypothesis and how the data were generated, then to claim that rejecting $\\mathbb{P}_0$ says something about the world. \n",
    "That logic is flawed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Randomized tests\n",
    "\n",
    "Sometimes it is useful for a hypothesis test to depend not only on the data but also on \"auxilliary\" randomness, typically a $U[0, 1]$ variable $U$ that is independent of the data $X$. \n",
    "Then the acceptance region $A$ is a subset of the Cartesian product of the data space and $[0, 1]$. \n",
    "The null hypothesis is not rejected if $(X, U) \\in A$; otherwise it is rejected.\n",
    "\n",
    "Randomized tests arise in a number of situations, including tests involving discrete distributions.\n",
    "Randomized tests are also useful for proving theorems about tests. \n",
    "For instance, _the Neyman-Pearson lemma_ shows that the most powerful test of a simple null hypothesis $X \\sim \\mathbb{P}_0$ against a simple alternative hypothesis $X \\sim \\mathbb{P}_1$ at significance level $\\alpha$ is\n",
    "a randomized test of the form:\n",
    "\n",
    "\\begin{equation*}\n",
    "A = \\{ (x, u): \\mathcal{L}_0(x)/\\mathcal{L}_1(x) > c \\} \\cup \\{ (x, u): (\\mathcal{L}_0(x)/\\mathcal{L}_1(x) = c) \\mbox{ and }\n",
    "  (u \\le d) \\}\n",
    "\\end{equation*}\n",
    "for suitable choices of $c$ and $d$,\n",
    "where $\\mathcal{L}_j(x)$ is the likelihood of hypothesis $j$ for data $x$. \n",
    "(This assumes $\\mathbb{P}_0$ and $\\mathbb{P}_1$ are _absolutely continuous_ with respect to each other.)\n",
    "For a definition of the likelihood function, see [Bayesian and Frequentist Estimation and Inference](./bayes.ipynb).\n",
    "\n",
    "A test based on a test function $\\phi: \\mathcal{X} \\rightarrow [0, 1]$ is implicitly\n",
    "a randomized test if $\\phi$ can have a value strictly between 0 and 1.\n",
    "Then, we also use an auxilliary randomness, $U \\sim U[0,1]$ and reject the null\n",
    "if $U \\le \\phi(X)$.\n",
    "This is completely equivalent to the description above, but instead of thinking\n",
    "in terms of a _set_ $A \\subset \\mathcal{X} \\times [0, 1]$, we think of a _function_\n",
    "$\\phi: \\mathcal{X} \\rightarrow [0, 1]$.\n",
    "In both formulations, we use $X$ and $U$ to decide whether to reject $H_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $P$-values\n",
    "\n",
    "Here are two approaches to defining $P$-values, in terms of families of hypothesis tests and in terms of a test statistic.\n",
    "\n",
    "Family of tests:\n",
    "+ Suppose you have a set of nested (monotone) hypothesis tests:\n",
    "\n",
    "    - $\\{A_\\alpha : \\alpha \\in (0, 1] \\}$\n",
    "\n",
    "    - $\\mathbb{P}_0 \\{ X \\notin A_\\alpha \\} \\le \\alpha$ (or more generally, $\\mathbb{P} \\{ X \\notin A_\\alpha \\} \\le \\alpha, \\; \\forall \\mathbb{P} \\in \\mathcal{P}_0$)\n",
    "\n",
    "    - $A_\\alpha \\subset A_\\beta$ if $\\beta < \\alpha$ (Can always re-define $A_\\alpha \\leftarrow \\cup_{\\beta \\ge \\alpha } A_\\beta$)\n",
    "    \n",
    "+ If we observe $X = x$, the $P$-value is $\\sup \\{ \\alpha: x \\in A_\\alpha \\}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test statistic definition of a $P$-value:\n",
    "\n",
    "If $P = P(X)$ is a random variable whose probability distribution is _dominated_ by the uniform distribution on $[0, 1]$ when\n",
    "the null hypothesis is true, then $P$ is a $P$-value.\n",
    "\n",
    "That is, $P$ is a $P$-value if\n",
    "\\begin{equation*}  \n",
    " \\mathbb{P}_0 \\{ P(X) \\le x \\} \\le x \\;\\; \\forall x \\in [0, 1].\n",
    "\\end{equation*}\n",
    "\n",
    "You may hear someone say that a $P$-value is the chance that a test statistic is \"as extreme or more extreme than\n",
    "observed.\"\n",
    "That is not really a precise mathematical definition; moreover, not every $P$-value can be expressed naturally\n",
    "in that form and not everything of that form is a $P$-value.\n",
    "\n",
    "For randomized tests expressed in terms of acceptance regions, $P(X, U)$ is a $P$-value if\n",
    "\\begin{equation*}  \n",
    " \\mathbb{P}_0 \\{ P(X,U) \\le x \\} \\le x \\;\\; \\forall x \\in [0, 1].\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The $Z$-test\n",
    "\n",
    "A $Z$-test tests the hypothesis that some function of the data has a standard Normal distribution, that is,\n",
    "the hypothesis that for some given test statistic $f(x)$, $f(X) \\sim N(0,1)$.\n",
    "\n",
    "What kinds of things have a standard normal distribution?\n",
    "\n",
    "+ A random draw from a population that has a standard normal distribution.\n",
    "+ The sample mean of $n$ IID draws from a population that has a standard normal distribution, after multiplying the sample mean by $\\sqrt{n}$.\n",
    "+ The sample mean of $n$ IID draws from a population that has a $N(\\mu, \\sigma^2)$ distribution, after subtracting $\\mu$ and multiplying by $\\sqrt{n}/\\sigma$.\n",
    "\n",
    "The $Z$-test is often used as an _approximate_ test rather than an _exact_ test. \n",
    "An _approximate_ test is one that has approximately its nominal significance level.  \n",
    "What kinds of things have distributions that are approximately standard normal?\n",
    "\n",
    "+ Things that have a Binomial$(n, p)$ distribution, after subtracting $np$ and dividing by $\\sqrt{np(1-p)}$, provided $np$ and $n(1-p)$ are not small. For instance, consider tossing a fair coin 100 times, independently. The number $X$ of heads has a Binomial$(100, 1/2)$ distribution. The distribution of $(X-50)/5$ is approximately a standard normal. ($np(1-p) = 100\\times 1/2 \\times 1/2$, so $\\sqrt{np(1-p)} = 5$.)\n",
    "+ The sample mean of $n$ IID draws with replacement from a finite population of numbers, after subtracting the population mean $\\mu$ and dividing by $\\sigma/\\sqrt{n}$, where $\\sigma$ is the population standard deviation, provided $n$ is sufficiently large and the population distribution is sufficiently \"bell shaped.\" If nothing is known about the population, it is in general impossible to know how accurate the normal approximation to the sample mean is.\n",
    "(The Binomial distribution is a special case where the population values are known to be 0 and 1.)\n",
    "\n",
    "When the test statistic is only approximately normally distributed under the null hypothesis, the accuracy\n",
    "of the approximation matters--but is rarely addressed.\n",
    "\n",
    "### Example of an approximate $Z$-test with a built-in Type III error\n",
    "\n",
    "There are two binary populations, $\\{x_j\\}_{j=1}^n$ and $\\{y_j\\}_{j=1}^m$. (A population is _binary_ if the only possible values in the population are 0 and 1.)\n",
    "We are interested in whether the populations are \"surprisingly different.\"\n",
    "The null hypothesis is that the two populations were formed by selecting $n$ items at random from the overall group of $n+m$ items to form the first population, with the remaining $m$ items comprising the second population.\n",
    "\n",
    "Let $\\bar{x} := \\frac{1}{n} \\sum_j x_j$ and $\\bar{y} := \\frac{1}{m} y_j$.\n",
    "If $x$ and $y$ were independent random samples with replacement from the same binary \"super-population\" that had a fraction $p$ of\n",
    "1s and a fraction $(1-p)$ of zeros,\n",
    "$n\\bar{x}$ would be a random variable with a binomial distribution with parameters $n$ and $p$, and $m\\bar{y}$ would \n",
    "be a random variable with\n",
    "a binomial distribution with parameters $m$ and $p$, and would be independent of $n\\bar{x}$.\n",
    "\n",
    "The expected value of $\\bar{x}$ would be $p$ and its variance would be $p(1-p)/n$;\n",
    "the expected value of $\\bar{y}$ would be $p$ and its variance would be $p(1-p)/m$.\n",
    "The expected value of $\\bar{x} - \\bar{y}$ would be 0, and its variance would be $p(1-p)(1/n+1/m)$.\n",
    "Moreover, if $x$ and $y$ were random samples with replacement from the same binary super-population,\n",
    "$\\hat{p} = (n\\bar{x} + m\\bar{y})/(n+m)$ would be an unbiased estimate of $p$,\n",
    "and for sufficiently large $m$ and $n$, the distribution of \n",
    "\\begin{equation*}\n",
    "f(x,y) := \\frac{\\bar{x} - \\bar{y}}{\\sqrt{\\hat{p}(1-\\hat{p})(1/n+1/m)}}\n",
    "\\end{equation*}\n",
    "would be approximately $N(0,1)$. \n",
    "(The accuracy of the approximation would depend on $n$, $m$, and $p$.)\n",
    "\n",
    "If we define the acceptance region $A := \\{ x, y: |f(x, y)| \\le z_{\\alpha/2}$, where $z_\\alpha$ is the $1-\\alpha$ \n",
    "percentage point of the standard normal distribution, we get a (two-sided) $Z$ test at nominal significance\n",
    "level $\\alpha$.\n",
    "\n",
    "But this is an approximate test of a different hypothesis: an approximate answer to a different question,\n",
    "a Type III error.\n",
    "The hypothesis test has little to do with the original hypothesis.\n",
    "The true significance level of the test for the original null hypothesis that the two groups are a random partition\n",
    "of the $n+m$ items could be quite different from $\\alpha$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numerical example\n",
    "\n",
    "There are two binary populations, $\\{x_j\\}_{j=1}^N$ and $\\{ y_j\\}_{j=1}^M$. \n",
    "Let $p_x$ denote the mean of the first population and $p_y$ the mean of the second.\n",
    "We wish to know whether $p_x = p_y$.\n",
    "\n",
    "A random sample with replacement of size $n = 100$ will be drawn from the first list and a random sample with replacement of size $m = 300$ will be drawn from the second list, independent of the other sample.\n",
    "Let $X$ denote the sum of the numbers in the sample from the first population and let $Y$ denote the sum of the numbers in the sample from the second population.\n",
    "Then $X \\sim \\mbox{Binom}(n, p_x)$, $Y \\sim \\mbox{Binom}(m, p_y)$, and $X$ and $Y$ are independent.\n",
    "\n",
    "(An equivalent problem is that there are two coins, one with chance $p_x$ of landing heads and one with chance $p_y$ of landing heads. The first coin is tossed $n$ times and the second coin is tossed $m$ times.\n",
    "All $m+n$ tosses are independent. Let $X$ be the number of times the first coin lands heads and let $Y$\n",
    "denote the number of times the second coin lands heads.)\n",
    "\n",
    "If $p_x = p_y = p$, $X+Y \\sim \\mbox{Binom}(n+m, p)$ and $\\hat{p} := (X+Y)/(n+m)$ is an unbiased estimate of $p$ with standard deviation $\\sqrt{p(1-p)/(m+n)} \\le 1/(2 \\sqrt{m+n}) = 0.025$.\n",
    "\n",
    "Let $\\bar{X} := X/n$ and $\\bar{Y} := Y/m$.\n",
    "If $p_x = p_y = p$ and $p$ is not close to 0 or 1, then the distribution of $\\bar{X}$ is approximately Gaussian\n",
    "with mean $p$ and standard deviation $\\sqrt{p(1-p)/n}$ and the distribution of $\\bar{Y}$ is approximately Gaussian\n",
    "with mean $p$ and standard deviation $\\sqrt{p(1-p)/m}$, and $\\bar{X}$ and $\\bar{Y}$ are independent.\n",
    "\n",
    "It follows that the distribution of $\\bar{X}-\\bar{Y}$ is approximately Gaussian with mean $p-p = 0$ and standard deviation \n",
    "\\begin{eqnarray*}\n",
    "\\sigma := \\sqrt{p(1-p)/n + p(1-p)/m} = \\sqrt{p(1-p)}\\sqrt{1/n + 1/m} \\le 0.0577.\n",
    "\\end{eqnarray*}\n",
    "Thus the distribution of \n",
    "\\begin{equation*}\n",
    "\\frac{\\bar{X} - \\bar{Y}}{\\sqrt{p(1-p)}\\sqrt{1/n + 1/m}}\n",
    "\\end{equation*}\n",
    "is approximately the standard normal.\n",
    "If we use the \"plug-in\" estimator of $p$, $\\hat{p}$, the distribution is also approximately normal (but the accuracy of the approximation is, in general, worse). \n",
    "Thus the distribution of\n",
    "\\begin{equation*}\n",
    "Z := \\frac{\\bar{X} - \\bar{Y}}{\\sqrt{\\hat{p}(1-\\hat{p})}\\sqrt{1/n + 1/m}}\n",
    "\\end{equation*}\n",
    "is approximately $N(0,1)$ if $p_x = p_y$.\n",
    "\n",
    "Note that if $p$ is close to 0 or to 1, $\\hat{p}$ could end up equal to 0 or 1, in which case the denominator will\n",
    "vanish (clearly underestimating the standard deviation of $\\bar{X}-\\bar{Y}$).\n",
    "\n",
    "To perform a $Z$ test, we pretend that the distribution of $Z$ is exactly the standard normal.\n",
    "If $p_x > p_y$, we expect $\\bar{X} > \\bar{Y}$ and hence $Z > 0$. If $p_x < p_y$, we expect $\\bar{X} < \\bar{Y}$ and hence $Z < 0$.\n",
    "If we want the test to have power against both the possibility that $p_x > p_y$ and the possibility that $p_x < p_y$, we would design the test to reject when $|Z|$ is large.\n",
    "To test at (approximate) level $\\alpha$, we would reject with $|Z| > z_{1-\\alpha/2}$, the $1-\\alpha/2$ quantile of the standard normal distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max SD of alpha-hat: 0.0005\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:52: RuntimeWarning: invalid value encountered in double_scalars\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p=0.01 estimated alpha=0.033563\n",
      "p=0.05 estimated alpha=0.039816\n",
      "p=0.1 estimated alpha=0.049005\n",
      "p=0.3 estimated alpha=0.050683\n",
      "p=0.4 estimated alpha=0.052481\n",
      "p=0.5 estimated alpha=0.057022\n",
      "p=0.6 estimated alpha=0.052372\n",
      "p=0.7 estimated alpha=0.050798\n",
      "p=0.9 estimated alpha=0.049088\n",
      "p=0.95 estimated alpha=0.039562\n",
      "p=0.99 estimated alpha=0.03344\n"
     ]
    }
   ],
   "source": [
    "# simulating the significance level of the Z-test\n",
    "import math\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy.stats import norm # the normal distribution\n",
    "from numpy.random import default_rng, binomial  # this is the Mersenne Twister; there's much to consider in picking a PRNG.\n",
    "rng = default_rng(1592592021)   # set the seed for reproducibility\n",
    "\n",
    "alpha = 0.05 # significance level 5%\n",
    "z_c = norm.ppf(1-alpha/2)  # normal quantile, critical value for the test\n",
    "\n",
    "p = [.01, .05, .1, .3, .4, .5, .6, .7, .9, .95, .99]  # assortment of values for p\n",
    "\n",
    "reps = int(10**6)\n",
    "n = 50  \n",
    "m = 100 \n",
    "\n",
    "def absZ(n: int, m: int, x: float, y: float) -> float:\n",
    "    '''\n",
    "    absolute value of the Z statistic for the difference in sample sums from two binary populations\n",
    "    \n",
    "    Used in a Z test of the hypothesis that x and y are the sample sums of IID draws from binary\n",
    "    populations with the same (unspecified) population fraction of 1s against the alternative\n",
    "    hypothesis that the two population fractions differ.\n",
    "    \n",
    "    When the population percentages are close to 0 or 1, the sample percentages can be 0 or 1,\n",
    "    in which case the estimated variance will be zero. In that situation, return -inf or inf, \n",
    "    as appropriate.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    n : int\n",
    "        sample size from first population\n",
    "    m : int\n",
    "        sample size from second population\n",
    "    x : float\n",
    "        sample sum of the draws from the first population\n",
    "    y : float\n",
    "        sample sum of the draws from the second population\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    z : float\n",
    "        absolute value of the Z statistic\n",
    "        \n",
    "    '''\n",
    "    pHat = (x+y)/(n+m)\n",
    "    return (abs((x/n - y/m)/math.sqrt(pHat*(1-pHat)*(1/n + 1/m)))\n",
    "           if pHat*(1-pHat) > 0\n",
    "           else np.inf*np.sign(x/n - y/m))\n",
    "\n",
    "# how precise should we expect the result to be?\n",
    "print('max SD of alpha-hat: {}'.format((1/(2*math.sqrt(reps)))))  # max SD of a binomial is 1/2\n",
    "\n",
    "for pj in p:\n",
    "    reject = 0\n",
    "    for i in range(reps):\n",
    "        X = rng.binomial(n, pj) # simulate X\n",
    "        Y = rng.binomial(m, pj) # simulate Y. X and Y could be simulated together but the code would be less clear\n",
    "        reject += (1 \n",
    "                  if absZ(n,m,X,Y) > z_c\n",
    "                  else 0)\n",
    "    print('p={} estimated alpha={}'.format(pj, reject/reps))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the simulation estimate of the actual significance level differs from the nominal significance level, 5%, by more than simulation variability accounts for: the test is only approximate, and the accuracy of the approximation depends on the true value of $p$. The actual significance level is below 5% when the true $p$ is near 0 or 1, and above 5% when the true $p$ is near 1/2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's look at the distribution of approximate $P$-values for (a one-sided version of) this test. \n",
    "If $Z$ really had a standard normal distribution, then for $z \\ge 0$,\n",
    "\\begin{equation*}\n",
    "\\Pr \\{|Z| \\le z \\} = 1 - 2\\Phi(-z) =  2\\Phi(z) - 1,\n",
    "\\end{equation*}\n",
    "where $\\Phi$ is the standard normal CDF.\n",
    "Let $x := 2\\Phi(z) - 1$.\n",
    "Then \n",
    "\\begin{equation*}\n",
    "x = \\Pr \\{ |Z| \\le \\Phi^{-1}((x+1)/2) \\}\n",
    "\\end{equation*}\n",
    "\\begin{equation*}\n",
    " = \\Pr \\{ \\Phi(|Z|) \\le (x+1)/2 \\}\n",
    "\\end{equation*}\n",
    "\\begin{equation*}\n",
    " = \\Pr \\{ 2\\Phi(|Z|)-1 \\le x \\}.\n",
    "\\end{equation*}\n",
    "Thus $2\\Phi(|Z|)-1$ would have a uniform distribution: it would be a $P$-value.\n",
    "\n",
    "Let's check whether that's true."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([136995.,  89683.,  86794.,  82063.,  78132., 139366.,  60161.,\n",
       "        136581.,  95286.,  94939.]),\n",
       " array([0.        , 0.09999988, 0.19999976, 0.29999964, 0.39999952,\n",
       "        0.49999939, 0.59999927, 0.69999915, 0.79999903, 0.89999891,\n",
       "        0.99999879]),\n",
       " <BarContainer object of 10 artists>)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD4CAYAAADy46FuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWCUlEQVR4nO3df6zd9X3f8edrdsPIMogBw5hNZ1q8toAaNXiO125VNm/gJFXMJNictcXKPFllrM2mTQ200pASWQJtKh3doLKCh2FZwKLZ8JbS1IJl2VR+5JImIUApdyEDFzfcxB5l6aA1ee+P87F2fHP8udf33B/c3OdDOjrf8/5+Pp/z+cjovO73xzmkqpAk6VT+zFJPQJL01mZQSJK6DApJUpdBIUnqMigkSV2rl3oC8+28886rDRs2LPU0JGlZefLJJ79ZVWtH7fueC4oNGzYwMTGx1NOQpGUlyf861T5PPUmSugwKSVKXQSFJ6jIoJEldBoUkqWvGoEiyL8krSb46Yt8/T1JJzhuq3ZRkMslzSa4aql+R5Km27/YkafUzktzf6o8n2TDUZ2eS59tj59irlSSdttkcUdwNbJteTHIR8LeBF4dqlwI7gMtanzuSrGq77wR2Axvb48SYu4BjVXUJcBtwaxvrHOBm4D3AZuDmJGtOb3mSpHHNGBRV9Xng6IhdtwG/CAz/Tvl24L6qeqOqXgAmgc1JLgTOqqpHa/C75vcAVw/12d+2HwC2tqONq4BDVXW0qo4BhxgRWJKkhTWnaxRJPgj8QVV9edqudcBLQ68Pt9q6tj29flKfqjoOvAqc2xlr1Hx2J5lIMjE1NTWXJUmSTuG0v5md5O3ALwNXjto9olad+lz7nFys2gvsBdi0aZP/Jya9ZW248TNL8r5fv+UDS/K++t4wlyOKHwQuBr6c5OvAeuCLSf4Cg7/6Lxpqux54udXXj6gz3CfJauBsBqe6TjWWJGkRnXZQVNVTVXV+VW2oqg0MPtDfXVV/CBwEdrQ7mS5mcNH6iao6AryWZEu7/nAd8GAb8iBw4o6ma4BH2nWMzwJXJlnTLmJf2WqSpEU046mnJJ8C3gucl+QwcHNV3TWqbVU9neQA8AxwHLihqt5su69ncAfVmcBD7QFwF3BvkkkGRxI72lhHk3wc+EJr97GqGnVRfV55akCSTjZjUFTVh2bYv2Ha6z3AnhHtJoDLR9RfB649xdj7gH0zzVHSW5d/fC1/fjNbktRlUEiSugwKSVKXQSFJ6jIoJEldBoUkqcugkCR1GRSSpC6DQpLUZVBIkroMCklSl0EhSeoyKCRJXQaFJKnLoJAkdRkUkqQug0KS1GVQSJK6DApJUpdBIUnqmjEokuxL8kqSrw7V/mWS30vylST/Mck7h/bdlGQyyXNJrhqqX5Hkqbbv9iRp9TOS3N/qjyfZMNRnZ5Ln22PnfC1akjR7szmiuBvYNq12CLi8qn4U+H3gJoAklwI7gMtanzuSrGp97gR2Axvb48SYu4BjVXUJcBtwaxvrHOBm4D3AZuDmJGtOf4mSpHHMGBRV9Xng6LTab1fV8fbyMWB9294O3FdVb1TVC8AksDnJhcBZVfVoVRVwD3D1UJ/9bfsBYGs72rgKOFRVR6vqGINwmh5YkqQFNh/XKP4B8FDbXge8NLTvcKuta9vT6yf1aeHzKnBuZ6zvkmR3kokkE1NTU2MtRpJ0srGCIskvA8eBT54ojWhWnfpc+5xcrNpbVZuqatPatWv7k5YknZY5B0W7uPxTwE+300kw+Kv/oqFm64GXW339iPpJfZKsBs5mcKrrVGNJkhbRnIIiyTbgo8AHq+qPh3YdBHa0O5kuZnDR+omqOgK8lmRLu/5wHfDgUJ8TdzRdAzzSguezwJVJ1rSL2Fe2miRpEa2eqUGSTwHvBc5LcpjBnUg3AWcAh9pdro9V1c9V1dNJDgDPMDgldUNVvdmGup7BHVRnMrimceK6xl3AvUkmGRxJ7ACoqqNJPg58obX7WFWddFFdkrTwZgyKqvrQiPJdnfZ7gD0j6hPA5SPqrwPXnmKsfcC+meYoSVo4fjNbktRlUEiSugwKSVKXQSFJ6jIoJEldBoUkqcugkCR1GRSSpC6DQpLUZVBIkroMCklSl0EhSeoyKCRJXQaFJKnLoJAkdRkUkqQug0KS1GVQSJK6DApJUpdBIUnqMigkSV0zBkWSfUleSfLVodo5SQ4leb49rxnad1OSySTPJblqqH5FkqfavtuTpNXPSHJ/qz+eZMNQn53tPZ5PsnPeVi1JmrXZHFHcDWybVrsReLiqNgIPt9ckuRTYAVzW+tyRZFXrcyewG9jYHifG3AUcq6pLgNuAW9tY5wA3A+8BNgM3DweSJGlxzBgUVfV54Oi08nZgf9veD1w9VL+vqt6oqheASWBzkguBs6rq0aoq4J5pfU6M9QCwtR1tXAUcqqqjVXUMOMR3B5YkaYHN9RrFBVV1BKA9n9/q64CXhtodbrV1bXt6/aQ+VXUceBU4tzPWd0myO8lEkompqak5LkmSNMp8X8zOiFp16nPtc3Kxam9VbaqqTWvXrp3VRCVJszPXoPhGO51Ee36l1Q8DFw21Ww+83OrrR9RP6pNkNXA2g1NdpxpLkrSI5hoUB4ETdyHtBB4cqu9odzJdzOCi9RPt9NRrSba06w/XTetzYqxrgEfadYzPAlcmWdMuYl/ZapKkRbR6pgZJPgW8FzgvyWEGdyLdAhxIsgt4EbgWoKqeTnIAeAY4DtxQVW+2oa5ncAfVmcBD7QFwF3BvkkkGRxI72lhHk3wc+EJr97Gqmn5RXZK0wGYMiqr60Cl2bT1F+z3AnhH1CeDyEfXXaUEzYt8+YN9Mc5QkLRy/mS1J6prxiEKSlqMNN35myd7767d8YMneeyEYFJI0z5YqpBYqoDz1JEnqMigkSV0GhSSpy6CQJHUZFJKkLoNCktRlUEiSuvwexVuEXw6S9FblEYUkqcugkCR1GRSSpC6vUeh77ndpJM0vjygkSV0GhSSpy6CQJHUZFJKkLoNCktRlUEiSusa6PTbJPwX+IVDAU8CHgbcD9wMbgK8Df7eqjrX2NwG7gDeBX6iqz7b6FcDdwJnAbwIfqapKcgZwD3AF8C3g71XV18eZs946vC1XWh7mfESRZB3wC8CmqrocWAXsAG4EHq6qjcDD7TVJLm37LwO2AXckWdWGuxPYDWxsj22tvgs4VlWXALcBt851vpKkuRn31NNq4MwkqxkcSbwMbAf2t/37gavb9nbgvqp6o6peACaBzUkuBM6qqkerqhgcQQz3OTHWA8DWJBlzzpKk0zDnoKiqPwD+FfAicAR4tap+G7igqo60NkeA81uXdcBLQ0McbrV1bXt6/aQ+VXUceBU4d/pckuxOMpFkYmpqaq5LkiSNMM6ppzUM/uK/GPiLwJ9L8jO9LiNq1an3+pxcqNpbVZuqatPatWv7E5cknZZxTj39LeCFqpqqqj8FPg38OPCNdjqJ9vxKa38YuGio/3oGp6oOt+3p9ZP6tNNbZwNHx5izJOk0jXPX04vAliRvB/4vsBWYAL4N7ARuac8PtvYHgf+Q5FcYHIFsBJ6oqjeTvJZkC/A4cB3wa0N9dgKPAtcAj7TrGNKcLeX/JEpajuYcFFX1eJIHgC8Cx4HfBfYC7wAOJNnFIEyube2fTnIAeKa1v6Gq3mzDXc//vz32ofYAuAu4N8kkgyOJHXOdryRpbsb6HkVV3QzcPK38BoOji1Ht9wB7RtQngMtH1F+nBY0kaWn4zWxJUpdBIUnqMigkSV0GhSSpy6CQJHUZFJKkLoNCktRlUEiSugwKSVKXQSFJ6jIoJEldBoUkqcugkCR1GRSSpC6DQpLUZVBIkroMCklSl0EhSeoyKCRJXQaFJKlrrKBI8s4kDyT5vSTPJvmrSc5JcijJ8+15zVD7m5JMJnkuyVVD9SuSPNX23Z4krX5Gkvtb/fEkG8aZryTp9I17RPGvgd+qqh8G3gU8C9wIPFxVG4GH22uSXArsAC4DtgF3JFnVxrkT2A1sbI9trb4LOFZVlwC3AbeOOV9J0mmac1AkOQv4SeAugKr6k6r638B2YH9rth+4um1vB+6rqjeq6gVgEtic5ELgrKp6tKoKuGdanxNjPQBsPXG0IUlaHKvH6PsDwBTw75K8C3gS+AhwQVUdAaiqI0nOb+3XAY8N9T/can/atqfXT/R5qY11PMmrwLnAN8eYt7TibLjxM0s9BS1j45x6Wg28G7izqn4M+DbtNNMpjDoSqE691+fkgZPdSSaSTExNTfVnLUk6LeMExWHgcFU93l4/wCA4vtFOJ9GeXxlqf9FQ//XAy62+fkT9pD5JVgNnA0enT6Sq9lbVpqratHbt2jGWJEmabs5BUVV/CLyU5IdaaSvwDHAQ2NlqO4EH2/ZBYEe7k+liBhetn2inqV5LsqVdf7huWp8TY10DPNKuY0iSFsk41ygAfh74ZJK3AV8DPswgfA4k2QW8CFwLUFVPJznAIEyOAzdU1ZttnOuBu4EzgYfaAwYXyu9NMsngSGLHmPOVJJ2msYKiqr4EbBqxa+sp2u8B9oyoTwCXj6i/TgsaSdLS8JvZkqQug0KS1GVQSJK6DApJUpdBIUnqMigkSV0GhSSpy6CQJHUZFJKkLoNCktRlUEiSugwKSVKXQSFJ6jIoJEldBoUkqcugkCR1GRSSpC6DQpLUZVBIkroMCklSl0EhSeoaOyiSrEryu0n+S3t9TpJDSZ5vz2uG2t6UZDLJc0muGqpfkeSptu/2JGn1M5Lc3+qPJ9kw7nwlSadnPo4oPgI8O/T6RuDhqtoIPNxek+RSYAdwGbANuCPJqtbnTmA3sLE9trX6LuBYVV0C3AbcOg/zlSSdhrGCIsl64APAJ4bK24H9bXs/cPVQ/b6qeqOqXgAmgc1JLgTOqqpHq6qAe6b1OTHWA8DWE0cbkqTFMe4Rxa8Cvwh8Z6h2QVUdAWjP57f6OuCloXaHW21d255eP6lPVR0HXgXOnT6JJLuTTCSZmJqaGnNJkqRhcw6KJD8FvFJVT862y4hadeq9PicXqvZW1aaq2rR27dpZTkeSNBurx+j7E8AHk7wf+LPAWUn+PfCNJBdW1ZF2WumV1v4wcNFQ//XAy62+fkR9uM/hJKuBs4GjY8xZknSa5nxEUVU3VdX6qtrA4CL1I1X1M8BBYGdrthN4sG0fBHa0O5kuZnDR+ol2euq1JFva9YfrpvU5MdY17T2+64hCkrRwxjmiOJVbgANJdgEvAtcCVNXTSQ4AzwDHgRuq6s3W53rgbuBM4KH2ALgLuDfJJIMjiR0LMF9JUse8BEVVfQ74XNv+FrD1FO32AHtG1CeAy0fUX6cFjSRpafjNbElSl0EhSeoyKCRJXQaFJKnLoJAkdRkUkqQug0KS1GVQSJK6DApJUpdBIUnqMigkSV0GhSSpy6CQJHUZFJKkLoNCktRlUEiSugwKSVKXQSFJ6jIoJEldBoUkqcugkCR1zTkoklyU5L8meTbJ00k+0urnJDmU5Pn2vGaoz01JJpM8l+SqofoVSZ5q+25PklY/I8n9rf54kg1jrFWSNAfjHFEcB/5ZVf0IsAW4IcmlwI3Aw1W1EXi4vabt2wFcBmwD7kiyqo11J7Ab2Nge21p9F3Csqi4BbgNuHWO+kqQ5mHNQVNWRqvpi234NeBZYB2wH9rdm+4Gr2/Z24L6qeqOqXgAmgc1JLgTOqqpHq6qAe6b1OTHWA8DWE0cbkqTFMS/XKNopoR8DHgcuqKojMAgT4PzWbB3w0lC3w622rm1Pr5/Up6qOA68C5454/91JJpJMTE1NzceSJEnN2EGR5B3AbwD/pKr+qNd0RK069V6fkwtVe6tqU1VtWrt27UxTliSdhrGCIsn3MQiJT1bVp1v5G+10Eu35lVY/DFw01H098HKrrx9RP6lPktXA2cDRceYsSTo949z1FOAu4Nmq+pWhXQeBnW17J/DgUH1Hu5PpYgYXrZ9op6deS7KljXndtD4nxroGeKRdx5AkLZLVY/T9CeBngaeSfKnVfgm4BTiQZBfwInAtQFU9neQA8AyDO6ZuqKo3W7/rgbuBM4GH2gMGQXRvkkkGRxI7xpivJGkO5hwUVfU/GH0NAWDrKfrsAfaMqE8Al4+ov04LGknS0vCb2ZKkLoNCktRlUEiSugwKSVKXQSFJ6jIoJEldBoUkqcugkCR1GRSSpC6DQpLUZVBIkroMCklSl0EhSeoyKCRJXQaFJKnLoJAkdRkUkqQug0KS1GVQSJK6DApJUpdBIUnqWhZBkWRbkueSTCa5cannI0kryVs+KJKsAv4t8D7gUuBDSS5d2llJ0srxlg8KYDMwWVVfq6o/Ae4Dti/xnCRpxVi91BOYhXXAS0OvDwPvGW6QZDewu738P0meG+P9zgO+OUb/5WilrXmlrRdc84qQW8da81861Y7lEBQZUauTXlTtBfbOy5slE1W1aT7GWi5W2ppX2nrBNa8UC7Xm5XDq6TBw0dDr9cDLSzQXSVpxlkNQfAHYmOTiJG8DdgAHl3hOkrRivOVPPVXV8ST/GPgssArYV1VPL+BbzssprGVmpa15pa0XXPNKsSBrTlXN3EqStGIth1NPkqQlZFBIkrpWZFDM9JMgGbi97f9KkncvxTzn0yzW/NNtrV9J8jtJ3rUU85xPs/3plyR/JcmbSa5ZzPkthNmsOcl7k3wpydNJ/ttiz3G+zeK/7bOT/OckX25r/vBSzHO+JNmX5JUkXz3F/vn//KqqFfVgcEH8fwI/ALwN+DJw6bQ27wceYvAdji3A40s970VY848Da9r2+1bCmofaPQL8JnDNUs97Ef6d3wk8A3x/e33+Us97Edb8S8CtbXstcBR421LPfYw1/yTwbuCrp9g/759fK/GIYjY/CbIduKcGHgPemeTCxZ7oPJpxzVX1O1V1rL18jMH3VZaz2f70y88DvwG8spiTWyCzWfPfBz5dVS8CVNVyX/ds1lzAn08S4B0MguL44k5z/lTV5xms4VTm/fNrJQbFqJ8EWTeHNsvJ6a5nF4O/SJazGdecZB3wd4BfX8R5LaTZ/Dv/ZWBNks8leTLJdYs2u4UxmzX/G+BHGHxR9yngI1X1ncWZ3pKY98+vt/z3KBbAjD8JMss2y8ms15PkbzAIir+2oDNaeLNZ868CH62qNwd/bC57s1nzauAKYCtwJvBokseq6vcXenILZDZrvgr4EvA3gR8EDiX571X1Rws8t6Uy759fKzEoZvOTIN9rPxsyq/Uk+VHgE8D7qupbizS3hTKbNW8C7mshcR7w/iTHq+o/LcoM599s/9v+ZlV9G/h2ks8D7wKWa1DMZs0fBm6pwQn8ySQvAD8MPLE4U1x08/75tRJPPc3mJ0EOAte1uwe2AK9W1ZHFnug8mnHNSb4f+DTws8v4r8thM665qi6uqg1VtQF4APhHyzgkYHb/bT8I/PUkq5O8ncEvMT+7yPOcT7NZ84sMjqBIcgHwQ8DXFnWWi2veP79W3BFFneInQZL8XNv/6wzugHk/MAn8MYO/SJatWa75XwDnAne0v7CP1zL+5c1Zrvl7ymzWXFXPJvkt4CvAd4BPVNXI2yyXg1n+O38cuDvJUwxOy3y0qpbtz48n+RTwXuC8JIeBm4Hvg4X7/PInPCRJXSvx1JMk6TQYFJKkLoNCktRlUEiSugwKSVKXQSFJ6jIoJEld/w/Q5hcum1FQogAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "reps = int(10**6)\n",
    "n = 50\n",
    "m = 100\n",
    "p = 0.5  # one of the worst offenders\n",
    "\n",
    "x = np.zeros(reps)\n",
    "for i in range(reps):\n",
    "    X = rng.binomial(n, p) # simulate X\n",
    "    Y = rng.binomial(m, p) # simulate Y\n",
    "    x[i] = 2*norm.cdf(absZ(n,m,X,Y)) - 1\n",
    "    \n",
    "plt.hist(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEJCAYAAACE39xMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAt00lEQVR4nO3dd3hUZdrH8e+dTm+hho6ACgoiAmIBQaQoYBdxxXVdWXftZdW1ri77WtaCCBZU1t5RqSqgi4gC0nsx1IRmKElIz8zc7x9nwCGkTGCSyczcn+viIjPzzDn3c+bMLyenPEdUFWOMMaEvKtgFGGOMCQwLdGOMCRMW6MYYEyYs0I0xJkxYoBtjTJiwQDfGmDBhgW6MMWEi7AJdRNaKSN8KmvY2EbmwIqZdzLzeFpExx9M2kMvAd1qB7n9FflalzLOjiCwXkUMickdlzruYWgLS/8pcLyu7hqr0eYWCkAx0ETlXRH4WkQwROSAiP4nIWQCq2klV5wa5vqB+wfxZBv7WGKjlWdz8gvRZ3Q/MVdVaqjqukud9lMrqv4jUExEVkSwRyRGRXSJyV0XPN0AC/nkFYnmISH0R+VJEskVku4iMLKXtXBHJ884vS0Q2nnAnShBygS4itYHpwMtAfSAJeALID2Zd4UZEYoJdQwVpBawNdhGVrCuQpqo1VbU68FfgRRFpHtyy/HLcn1cp63BXTnx5TAAKgMbAdcCrItKplPa3eedXU1U7lmM+5RJygQ50AFDVj1TVraq5qjpLVVfB0VuC3p//LiKrvL9J3xKRxiLytfdPuDkiUu/whL2/tU/yeVzibg8ReVBENnuns05ELvM+/x7QEpjm/W18v/f5ZiIyWUTSRGRr0T8fReQMEVnmnd4nQEJJC6CstkWWwQMistPbdqOI9C+lxm3e9quAbBGJKWbL+ixvfw+KyH9FJMFnvsUuvzLmd7jOU7xbMuneXRHDiunTfd7PMkNEPvGdd5G2xU5LRL4HLgDGe+voUML7S5xXaXUex/p21LItY77Frm9+6gos9nm8yPt/XJF+Pyginxd57iURGVeeGsr6HpX1XfBpd8zn5ed6ctQ6fLzLoyQiUgO4AnhUVbNUdT4wFbjen/dXKFUNqX9AbWA/8A4wGKhX5PVtwIU+Py/E+S2aBPwGLAPOAOKB74HHfd6rwEk+j98GxpQw7auAZji/FK8BsoGmRdt5H0cBS4HHcFaatsAWYKD39ThgO3A3EAtcCRT6zttnWmW2PTx/oCOQAjTzPt8aaFdcjT7PrQBaANVKWJ5rvK/XB34qMt8Sl18p87vQ249k4CFv//oBh4CORdr+4l3m9YH1wC3FLJ9SpwXMBf5cxjpW7Lz8mPY2yre+HbVMSusj5VjfiunPu4fnC9QF3gKWAFKkXSsgB6jtfRwN7AZ6laeGMtaDUr8LxdR+5PMqx3qyAp91uDzLA+ev//QS/k33vucMILfINO8DppXShzRgH853pm9F5WPIbaGraiZwLs5K8waQJiJTRaRxCW95WVX3qupO4EdgkaouV9V84EucD+d46vhMVXepqkdVPwF+BXqU0PwsoKGqPqmqBaq6xVv7CO/rvXBW1rGqWqiqn3P0FoSv8rR14wTJqSISq6rbVHVzGV0bp6opqppbwuvjva8fAP4NXFvG9PzRC6gJPO1dPt/jfLGKTnucd5kfAKbhbGkd77TKUty8/Jn2ia5vxfaxnOtbUV2Bv4vIAZwwVWCoetPmMFXdjvML6FLvU/2AHFVdGIAaDivru1Ca8qwnpa3DXSlheajqJapat4R/l3jfXxPIKDLNDKBWCfN7AOcXVxIwEecv1XZ+9LfcQnI/qaquB/4IICInA+8DYyn+S7vX5+fcYh7XPJ4aRGQUcA/OVi/e6SSW0LwV0ExE0n2ei8b5woOz1bOzyBdsewnT8rutqiaLc7Dnn0AnEfkWuEdVd5UwbXC26Evj+/p2bz0nqhmQoqqeItNOKtJuj8/POSXM299plaW4efkz7RNd34rtYznXtyNEJB44BWijqqlltQc+xPkevQuM9D4+PK3jqqGIsr4LpfH3sy1xHT6O5VGcLJw9Bb5q4/y1cAxVXeTz8B0RuRYYgnMcMKBCbgu9KFXdgPMnXecATC4HqO7zuElxjUSkFc5WxW1AA1Wti7MrQg6XVeQtKcDWIr/ta6nqEO/ru4EkERGf97QsocbytEVVP1TVc3G+SAo8U0KNlPH8YS2KzNf3l0Npy6+06e4CWoiI7/rYEthZRi0VPa3KnHaJ/FjfStMZyC5HeH0G9BXnAOFleAO9nDWUth6U9V0ojb/Lv7R1rdTl4T3ekVXCv6+9zTYBMSLS3uetXfD/4K3i32dXbiEX6CJysojc613hEJEWOFsUCwMw+RXASBGJFpFBQJ8S2tXA+VDSvDXcyNG/UPbi/Il12C9ApvdgTTXv9DuL91RLYAHgAu4Q50Dk5ZT8p6zfbcU5h7efd6skD2cL0V1Cjf66VUSai0h9nH2Zn/i8toKSl19p81uEsz/2fhGJFefc7KHAx8dRXyCnVZnTLk1Z61tpzqAcZ4moahrOPt//4gTv+uOoYQUlrwdlfRdKE4jlX+ryUNXB+vvZKEX/Dfa2yQa+AJ4UkRoicg4wHHiv6PREpK6IDBSRBO/39TrgfODbctTst5ALdJw/a3oCi0QkGyfI1wD3BmDad+KsIOk4pyJ9VVwjVV0HPI8TrnuB03AOdhz2FPCI90j8farq9k63K7AV5+DIm0Ad7/QKgMtxdiMdxDng9EUJ8/a7Lc7+86e989sDNMIJ4WNqLOH9xfkQmIVzIGsL4HsWUGnLr8T5efs0DOcg9z7gFWCU96+vcgnktCpz2mXMt6z1rTRdcb4f5fEhzsHqI7tbyllDietBWd+F0gRo+Xel/MujOH8DquEc+P4I+KuqroUjW/mHv2exON+RwwdFbwcuVdUKORddihwXMcYYE6JCcQvdGGNMMSzQjTEmTFigG2NMmLBAN8aYMBG0C4sSExO1devWwZq9McaEpKVLl+5T1YbFvRa0QG/dujVLliwJ1uyNMSYkiUhJV5HbLhdjjAkXFujGGBMmLNCNMSZMWKAbY0yYsEA3xpgwUWagi8gkEflNRIod0EYc40QkWZxbZ3ULfJnGGGPK4s8W+tvAoFJeHwy09/4bDbx64mUZY4wprzIDXVXnAQdKaTIceNd7+6aFQF0RaRqoAo0xJly4Cwv4+e2HSV4+r0KmH4h96EkcfcunVEq43ZeIjBaRJSKyJC0tLQCzNsaY0LA/eTE7nulF723j2bf4swqZRyACvbhbKRU7yLqqTlTV7qravWHDYq9cNcaYsFKQl8NPr99BnfcvombhPn4+80V63jyuQuYViEv/Uzn6PpPNOfo+k8YYE5G2LJ2DTLudc9jF99UG0OraF+ndskXZbzxOgQj0qcBtIvIxzq3hMlR1dwCma4wxISnn0EHmv34nFx6aym4SWXjOm/S76KoKn2+ZgS4iHwF9gUQRSQUex7lPHqr6GjATGAIk49zt+8aKKtYYY6q6OdM+oMvyf3KhO40FiVdQf9gYerWqnPNEygx0Vb22jNcVuDVgFRljTAhyZe1n1yf3cGHKV2yTJGZ3nsjIq66u1BqCNnyuMcaEA1Vl3PjnGblvHM04xCfVr+HSO19kZEKNSq/FAt0YY47XoT0cmnwXd+7/mpSEDnzTbSKDzutPfEJsUMqxQDfGmHJauu0AP342lj9lv0m8FvC0awQ33/MCl9Su/K1yXxboxhhTHge3E/fRH7krfxlba3RhRuuHaNumEw2CHOZggW6MMWUqdHtYk3KAGqv+S9uVL9DW5WFmq/sY8seHuS2q6gxaW3UqMcaYKsjtUe58+WPckwbTYdkY5he0Z0D+s9Q+769QhcIcbAvdGGNK5i4kefIYXjw4Hk9cDbb1fpHEDpfzXnwMbRODv4ulKAt0Y4wpzq7lMOV2Ou5dzXfR53DBnW/TulajYFdVKgt0Y4zx2rjnEE9NW8YFuydxnXsqB6jDI4V30+bca+hfxcMcLNCNMQYAj0d54Y1JPOF+lVbsZlniUOa2voMOcXUY1btVsMvziwW6McbkZXJo+iO87n6HrOrN4copdGvbl1C7n2bVOkRrjDGV7dfZ8MrZ1F7zLm+5BrPjmjnQtm+wqzoutoVujIkoBS4Pz8/ayOYdO7h63ytc5JrLVmnBo/yLRdqOtUmNg13icbNAN8ZElNfmJpM6/0Oej3+HmprFt4k3MCfxDzSOiueR5nWIiwndHRcW6MaYyJG5m56/3M4dcQuh6RkwbDwDm3RmYLDrCpDQ/VVkjDH+UoVl76ITetClYCnfNLsVbpoDTToHu7KAsi10Y0x4O7AVpt0BW+dR2Lw3AzdfyY2d+kF0+MVf+PXIGBPRClweXvpuEz9u3MOg7KncmP8+HqKYEHMLn+/pz29aSNuGNYNdZoWwQDfGhJXHpqxh6ZIFvFJzEu0LN7C2Ri8mN7uPzLhGnA/USoihR5v6wS6zQligG2PCxra9B2m0/CW+SfiK6Ng6MOwtOnW+gk4iwS6tUligG2PCw86l1P/kFu6J2URGu+HUuewFqJEY7KoqlZ3lYowJbQU5MOsRePNCJC+dPxfeS/w1/424MAfbQjfGhLKtPzpnsBzYAmf+kSczrmTDzkISYqODXVlQWKAbY0JCToGLp2ZuYPG2A2Rl7OcOfZ+rmcMObczj+igLf+lMniuT89o3DHapQWOBbowJCU9MXcenS1O4tVkyN8e8TK3C/fzceCRzm91Mu+gEThIQEQZ3bhLsUoPGAt0YU+UVuDx8t3QtUxpP5vT9s6HRqTDsE3o3P5PewS6uCrFAN8ZUbaocWvwR38Y9QN2MXOj7EJx7N8TEBbuyKscC3RhTdWXshBn30GDTN6zQdmwd/DLde5wT7KqqLAt0Y0zV4/HAsndg9mPgLmR15we4fMlpzGp9erArq9Is0I0xVYbbo7jSkomZeRfR2+fjbn0+riEvMne14mETTetUC3aJVZoFujEmaLbvz+b/Zq5n2Y50snLz+IPO4N6Yz8gmhn+7buaTDX1hw0YA6lSLpUa8RVZp/Fo6IjIIeAmIBt5U1aeLvF4HeB9o6Z3mc6r63wDXaowJI26PcvO7S9idkceotllct/c5mmWvY2uDPvzY4R+0jm/EAz7tT21WO2i1hooyA11EooEJwAAgFVgsIlNVdZ1Ps1uBdao6VEQaAhtF5ANVLaiQqo0xIW/epjS27T3ItK6L6LhpIiTUhSsn0abT5bSJkMG0As2fLfQeQLKqbgEQkY+B4YBvoCtQS0QEqAkcAFwBrtUYE0Y2LPmeGQljaL8hBU6/BgY+BTUaBLuskOZPoCcBKT6PU4GeRdqMB6YCu4BawDWq6ik6IREZDYwGaNmy5fHUa4wJdQXZ8P2/+UvyK+yPagAjPoUO4XJXz+DyZ7TF4v720SKPBwIrgGZAV2C8iByzw0tVJ6pqd1Xt3rBh5I63YEzE2vIDvNobFk5gRtwgnmgxycI8gPwJ9FSghc/j5jhb4r5uBL5QRzKwFTg5MCUaY0JebjpMvR3eHQYSDX+cyaOuP1Gvnu1iCSR/An0x0F5E2ohIHDACZ/eKrx1AfwARaQx0BLYEslBjTIjaMAMm9ITl78M5d8JffyK3WS/ScwppUich2NWFlTL3oauqS0RuA77FOW1xkqquFZFbvK+/BvwLeFtEVuPsonlAVfdVYN3GmCpm/e5MPl2SwvrdmWTkuqiWv5/bCt6gn2s+m6Na82zCs2xa2R5WLsLlcQ6xNaltgR5Ifp2HrqozgZlFnnvN5+ddwEWBLc0YEyp2Z+RyzesLyHd56NysNlfEzOfajFeI01ymJ/6JOfVHEi8xnAYcPiOxd9tE+nS0Y2mBZJddGWNO2Ns/byO7wM3cm0+ixc8Pw6+zoPlZMGw8lzQ6mUuCXWCEsEA3xpwQVWX68lQeb7KAFh/dDOqBQc9Aj5shKjJvBRcsFujGmBOSvXsDL+Q9Qs+CDdC2Lwx9Ceq1DnZZEcmfs1yMMeZYbhfMH0v1t/pysuxgadcxcP1XFuZBZFvoxpjy27MaptwKu1dysMVFDPp1OGM7D/r9iKcJCgt0Y4z/XPkw7z8w/0WoVg+ueof5Bd1J+3UljWrFB7u6iGeBbozxz45FztWe+zZCl5Ew8N9QvT5pPzrXEDaqZeeUB5sFujHmiLxCN1NW7CT1YC6qoCgxrhzO2/EqZ+79jMy4xszs8CKbY3uh3+/Bo7tZtv0gcTFR1K5mcRJs9gkYYwD4bv1exsxYz9Z92QBECZwbtZp/x7xJC0njPc9AxuaOIGd9dUR2ECXijNwn0KdDQ8T2nwedBboxBpfbw2NT1rIzPZfnr+rC5afUQGY9AivehwbtYdh7XN/qbK4PdqGmVBboxhiWp6SzMz2Xl689g6FxS+GVeyF7H5x7D/R5AGJt/3gosEA3xrArPZeGpNN35X2wZQY0OQ1GfgrNuga7NFMOFujGRDpV6m76jNnxz1Jzuwv6Pwa974Do2GBXZsrJAt2YSJa+A6bdRZ/N37GUjnS75T1o2DHYVZnjZIFuTCTyeGDxmzDnnwBMbnwnL6afz3wL85BmgW5MpNn3K0y5DVIWQrv+MHQsX3y+i4Z13MGuzJwgC3RjIoW7EH4eB3OfgdhqcOmr0OVaEGF3xhbaN6oZ7ArNCbJANyYS7F7pDKa1ZzWcOhwG/wdqNQZgT0YeW9KyuerMFmVMxFR1FujGhAFV5cvlO5m8LJVd6XkUuJx7dsZpATe6PmGk+yvSqc1zMffzv+Se8PJqVFejcKRtv5MbBbEHJhAs0I0JcR6P8vBXq/nolxRaNahOp2a1SYiNpm3Oaq7e9QyN3Cn8UncIM5reikbXpi/OKLdy+Lp9oHm9anRobLtcQp0FujEh7vsNv/HRLync0qcdfx/YkejCLJjzBKx7A+q2hGu+pEe7fvQIdqGmwlmgGxPi1u3OBODO/u2J3vIdTLsLMlKh5y3Q71GIty3vSGGBbkyI25yWxcl1XFSbcSus/AgSO8CfvoWWPYNdmqlkFujGhDJVGqd+w5OFr8HqLDjvPjj/7zaYVoSyQDcmVB3aAzPu5aGs6eyq1pE6o6ZB09ODXZUJoqhgF2CMKSdVWP4+TOiBJs/hqcJrmXrWexbmxgLdmJBycBu8d6lzkVCjTuy8Zjavu4fSsE6NYFdmqgDb5WJMKPC44Zc34LsnQKLg4ufhzD+xe0c6kELDWvHBrtBUARboxlR1v22AqbdD6i9w0gC45EWo61ymn3YoH4BGtS3QjQW6MVWGqqIK6v0ZdyHy80tE/fgfiKuJe/hr6GlXO5d4up3L9fdk5AHQsKYFuvEz0EVkEPASEA28qapPF9OmLzAWiAX2qWqfgFVpTBjamZ7L6z9sZt6mNPZlFZCV7zryWmfZwn9iJ3JK1A6muXvxz+wb2P9Jbfjkm2OmExcdRb3qcZVZuqmiygx0EYkGJgADgFRgsYhMVdV1Pm3qAq8Ag1R1h4jYKD/GlEJV+dv7S1m/5xDnt2/IBSdXo1ZCLLGePM7e8Qbddr5PTlx9pp/0HNsT+3JjKdM6qVEtoqKk0mo3VZc/W+g9gGRV3QIgIh8Dw4F1Pm1GAl+o6g4AVf0t0IUaE05mrdvLytQMnr3idK4+yzts7bafnH3lBzZDt1HUHPAvLqlWN6h1mtDiz2mLSUCKz+NU73O+OgD1RGSuiCwVkVHFTUhERovIEhFZkpaWdnwVGxMGXv9hM20Ta3B5tyTIy4Tp98DbQ8DjglFTYNjLYGFuysmfLfTi/pbTYqZzJtAfqAYsEJGFqrrpqDepTgQmAnTv3r3oNIyJGDsO5DDg1CbEbJ4D0++CzF3Q61bo9zDE2Tnl5vj4E+ipgO+tTJoDu4pps09Vs4FsEZkHdAE2YYw5isejkLOf63e9Cau+gYYnw02zocVZwS7NhDh/drksBtqLSBsRiQNGAFOLtJkCnCciMSJSHegJrA9sqcaEAVXyln/GN7F/5+T9c6DPg/CXeRbmJiDK3EJXVZeI3AZ8i3Pa4iRVXSsit3hff01V14vIN8AqwINzauOaiizcmJCTuRtm3EP1jTP5Vduysu/z9O/TL9hVmTDi13noqjoTmFnkudeKPP4P8J/AlWZMmFCFZe/CrEfBnU/KWQ9x+Y+n8FaTzsGuzIQZu1LUmIp0YAtMuxO2zoNW58KwcazfUwP3j0tpUMOu7jSBZYFuTEXwuGHhq/D9GIiKgUvGQrcbICqKg1t2AFCvRmxwazRhxwLdmEDbuw6m3gY7l0KHQXDxC1Dn90s3Ug7kEh0lJNr4KybALNCNCRRXAcx/AeY9Bwm14Yq3oPMVzmBaPhZvO0CnZrVJiI0OUqEmXFmgG+OHA9kFPP31ejbtzXJGRQQ83tERPQrtCzdyd8442ni2811MHybITaTPqgPf/uC082m/Kz2XG89pE+wumTBkgW5MGRZs3s89n64g7VA+vdo2IDpKiBIQERI0j+EH32FA1mQyouszockY1tbqTVOEpgJRIggcaS8CsVFRXNezZbC7ZcKQBboxZRg7ZxO7M/J4+dozGNql2e8vbJ0HU++AzK1w5o3UG/AEtybUCV6hJuJZoBtTioPZBSzedoA7+p30e5jnZcDsx2Dp21CvDdwwDdqcH9Q6jQELdGNKtW1/Nh6Fri3rOk9s/Bqm3w1Ze6H37dD3IYirHtQajTnMAt2YUhy+i1A9MuHzh2DN59CoE4z4AJLODHJ1xhzNAt2YUmTlFjIs6idO++pWKMhytsjPvRti7JZvpuqxQDemJBk76TL/rwyO+4H8Ot2IueIVaHRKsKsypkT+DJ9rTGTxeGDJJJjQk0b7FvGvwj+Qe/1MC3NT5VmgG+Nr/2Z4Z6hz4DPpDN7v9jFvuYdQo5pdpm+qPgt0YwDcLvhpHLzaG/asdu7pOWoqu6OakBAbRWy0fVVM1Wf70I3Zs8YZTGvXcuh4MVz8PNRuCsChfBc14+1rYkKDrakmcrny4cfnnX8JdeHK/0Kny44aTCsrzwLdhA5bU01kSlnsbJWnbYDTr4FBT0P1+sc0y8p3UTPBviYmNNiaasLKxj2HWLztAC63B8W5+xuA9z+iXTmcuXkCnVM+JDu+EfNPH8eOxHNhyUFUDx5pd/h9yb9l0axuQiX3wpjjY4FuwkJOgYu/fbCMuRvTSmzTO2oNT8e8QcuoNN51DeDZvGvI+qU6sKHUaV/QsWGAqzWmYligm7AweWkqczemcfeFHbiye3Oqe28eIQKSl0G1/z1O3OoPcNdrS/bgt7i0ZW8u875XvEPcHm4PcPgZEYiPsTNcTGiwQDdhYdqq3ZzatDZ39D8J8b1D0IYZMP0eyE6Dc+4iuu+D1IitFrxCjalAFugmLKQcyKF3u8TfwzzrN/j6flj7JTQ+DUZ+DM3OCG6RxlQwC3QT8grdHvZm5pFUr5pzNHPVJ/DNg1CQDf0egXPugujYYJdpTIWzQDchb09GHh6FDvHp8MFVkDwbmveA4eOhYcdgl2dMpbFANyEv9UA2f4iezaB5n4IoDHoGetwMUdHBLs2YSmWBbkLbvmTazxzNmNil5DQ9n5jLx0O9VsGuypigsPOxTGhyu2D+i/Bqb2plbuK+wr9QOGKyhbmJaBboJvTsXgVv9oM5/4T2A5jU5RO+8PShVjU78GkimwW6CR2FefDdkzCxL2TuhqvfhREfkOqqTd3qcURFSZmTMCac2T50Exp2LHIG09q3CbqMhIH/PjKYVnpOIXWr29a5MX5toYvIIBHZKCLJIvJgKe3OEhG3iFwZuBJNRMvPgpn3w6SBUJgLf5gMl7161MiIB3MKqFfdbtpsTJlb6CISDUwABgCpwGIRmaqq64pp9wzwbUUUaiJQ8ncw7S7ISHFOQ+z/GMTXOqbZwZxCkmxERGP82uXSA0hW1S0AIvIxMBxYV6Td7cBk4KyAVmgiT+5B+PZhWPEB2qA93DgTbXG2M7StR1Hv2LaHh7pNzymgc7PawarWmCrDn0BPAlJ8HqcCPX0biEgScBnQj1ICXURGA6MBWrZsWd5aTQjLzCvkgc9XMXdjGi6PB1V8xh7/PaAHRv3CkzFvU59MXncPY9zOy8l/9SAws9Tp169pu1yM8SfQizt1QIs8Hgs8oKruo0a6K/om1YnARIDu3bsXnYYJY898vYHZ6/YyokcLaic4BzB9h6qtUbiP/lueo8OB7/mtRgc+Omk8+TU7cgvFD2nrPHZERQmXnpFUeZ0xporyJ9BTgRY+j5sDu4q06Q587A3zRGCIiLhU9atAFGlCW1a+iy+X7+TybkmMufS0o19UhRUfwrcPOQc9+z9Oo963c70NpmVMufkT6IuB9iLSBtgJjABG+jZQ1TaHfxaRt4HpFubmsOU7DpJT4GZol2ZHv3BwO0y/CzZ/Dy3PhmEvQ2L7oNRoTDgoM9BV1SUit+GcvRINTFLVtSJyi/f11yq4RhPiVqVmAHB6Ul3nCY8HFr8Bc55w9p8MeQ663wRRdp2bMSfCrwuLVHUmRY5KlRTkqvrHEy/LhJPVqRm0blCdOtVjIW0TTL0dUhZCu/4wdCzUtQPkxgSCXSlqKtyujFza1o+Hec/BD89AbHW49DXoMuL3I5zGmBNmgW4qXJPsjTx66BVISYZThzu7WGo2CnZZxoQdC3RTcQpz4YdneCX3JXJi6sE178MpQ4NdlTFhywLdVIztC5zBtPYn85WnLztOf4h7TulZ9vuMMcfNTiswgZV/CGbcB/8dBO4CXNd9yX0Fo4mpUb/s9xpjTohtoZvA+XW2M5hW5k7o+Vfo9wiHXHHAbGrG26pmTEWzb5k5cTkH4Jt/wKqPIbEj3DQLWvQAICs7B4BaCbaqGVPR7Ftmjp8qrPsKZv7dGSHx/L87/2LijzQ5lOcCLNCNqQz2LTPH59AemHEvbJgOTbvC9V9Ck9OOaZaV7wR6zXgbm8WYimaBbspHFZa/74xX7s6HAU9Cr1shuvhV6VBeIQA1bQvdmApn37II9+veQ0z4XzL7sgrILnDhcpc8qnFj9x5uzX6ZMwpXsCamMy/XvINdy5Jg2cIS35OeWwBgB0WNqQT2LYtw9322klU7MzijRV1qxscQG33smayibvofmsIVGZPwEMW79e/kh1oXoxJFwzKm37BWPL3aNKBVg+oV0wFjzBEW6BGs0O1h/e5DjD6vLf8YckrxjX7b4FwgdGAxnDQAho5lVJ3mjKrcUo0xfrBAj2Cb07IocHs4tbj7cboK4KexMO8/EFcTLn8DTrvKBtMypgqzQI9gG/ccAuDkJkUCfecyZ4jbvWug8xUw6BmoWdbOFWNMsFmgR7D0HOcMlMTDN1guzIX//R8sGA81G8OIj+DkIUGs0BhTHhboESynwA1A9bgY2Dbf2So/sAW63eCcjlitbnALNMaUiwV6BMspcFFLckj49j5YOgnqtYZRU6Ftn2CXZow5DhboEazZb/OYHfcUsuwgnH0bXPAQxNUIdlnGmONkgR6JsvfDNw9ybfKnbJYWcNOn0Lx7sKsyxpwgC/RIogprJsPX90NeJl83uIHncy9hjoW5MWHBAj1SZO5yBtPaOBOadYPh45n8TS6xntxgV2aMCRAL9HCnCsvegVmPgrsQLhoDvf4GUdHkFCykRlx0sCs0xgSIBXo4O7AFpt4B236E1ufB0JegQbsjL2cXuKlTzYa1NSZcWKCHI48bFr4K34+B6Fi4ZKxzbnnU0QNv5Ra4aFo7ITg1GmMCzgI93Oxd5wymtXMpdBgEF78AdZKKbZqd76Z6vO1yMSZcWKCHC1cBzH8B5j0HCbXhireccVhKGUwrp8BFdduHbkzYsECvgvZm5vHugm3sO1RAgduDask3nQBombueEXuepVn+VpbUvpAvG91G9tq6sHZFqe/LzHNRI85WAWPChX2bq6A7PlrOoq0HaFw7nriYKKJK2MqO13xuKviAK13T2S/1eDD+YRa4z4LdAOllzqdl/er0bFs/oLUbY4LHAr2KSf4ti0VbD/DQkJMZfX67khtunQdTH4CcbXDmjTQa8ARPJ9SptDqNMVXPsfcbK4aIDBKRjSKSLCIPFvP6dSKyyvvvZxHpEvhSI0PaoXwATkuqW3yDvAznVMR3hgICN0yHoWPBwtyYiFfmFrqIRAMTgAFAKrBYRKaq6jqfZluBPqp6UEQGAxOBnhVRcLjLK3SGtE2ILeZ37cavYfrdkLUXet8Bff8BcXavTmOMw59dLj2AZFXdAiAiHwPDgSOBrqo/+7RfCDQPZJGR5HCgV/M9+yR7nzP+yprJ0KgTjPgQkroFqUJjTFXlT6AnASk+j1Mpfev7JuDr4l4QkdHAaICWLVv6WWJkyT0c6LHRzmX7qz+Drx+A/ENwwcNwzl0QExfcIo0xVZI/gV7cKRbFnkcnIhfgBPq5xb2uqhNxdsfQvXv30s/Fi1CHA7163h748B/w67eQ1B2Gj4dGpwS5OmNMVeZPoKcCLXweNwd2FW0kIqcDbwKDVXV/YMqLPHkFLq6LnkPiO6NBPTDwKej5F4iyC4CMMaXzJ9AXA+1FpA2wExgBjPRtICItgS+A61V1U8CrjBT7NzNoyc0kxS7F06wPMuwlqN8m2FUZY0JEmactqqoLuA34FlgPfKqqa0XkFhG5xdvsMaAB8IqIrBCRJRVWcThyu+Cnl+DV3jTI2sgDhaORUV9ZmBtjysWvC4tUdSYws8hzr/n8/Gfgz4EtLULsWQ1TboPdK6Djxbwa/xemr8jlmSi/LhEwxpgjLDWCxZXvDG87sS9k7oSr3oYRH7BP6h99yqIxxvjJLv0PhpRfnK3yfRvh9BEw6Cmo7oypklvoJj7GAt0YU34W6JWpIBu++xcseg1qJ8F1n0P7AUc1ySt02xa6Mea4WKBXls3/g2l3QPoOOOvP0P9xZ9zyInIL3M5FRcYYU04W6BUtNx1mPQzL34f67eDGr6FV7xKb5xV6ih/HxRhjymCBXpHWT4cZ90J2Gpx7N/R5AGKrlfqW3EI3tRLsYzHGlJ8lhx8+XZLCuO9+JbfA7Vf7+prO/foWA3QBG2nNk1H/x4YF7WDB/DLfm55bSP+TG51oycaYCGSB7oepK3aRejCXP/QqY0AxVboe/IYhO8cRp3nMbjqaHxuNpI3EUJ5LhIZ1Kf6mzsYYUxoLdD/sycxjYKfGjLn0tJIbpafA9Ltgxxxo3gOGj2dAw44MKPkdxhgTUBboftibkce5JyUW/6LHA0vegjn/dIa7HfyscxaLDaZljKlkFuhlyM53cSjfRePaCce+uO9XmHo77FgAbS+AoS9BvVaVX6QxxmCBXqY9mXkANKkT//uT7kL4+WWY+zTEJsDwV6DrSJDiho43xpjKYYFehr3eQG9cy7uFvnulc9n+nlVwylAY8jzUahzECo0xxmGBXobMXBcAdePc8N2TMH8sVG8AV78Lpw4PbnHGGOPDAr0M2fkuzpSNtP/iUUjfDF1GwsB/HxlMyxhjqgoL9NLkZ3HqyjF8Fvcx6mkOf5gMJ10Y7KqMMaZYNmhISZLnwCtnc/KOj3nHfRGFf/nJwtwYU6VZoBeVcwC+/Cu8fwXExPNRp9cZ4/kj8dWPHRnRGGOqEgt0X+umwISesOoTOO9euGU+m+I7UyMuGrFTEo0xVZztQwc4tAdm3gfrp0GT05195U1PByAr30XNeFtMxpiqL7KTShVWfAjf/gMK8+DCf8LZt0F07JEm2fkualigG2NCQOQm1cHtMO1O2PI/aHk2DHsZEtsf0ywr30V1C3RjTAiIvKTyuGHxmzDnCedS/SHPQfebIKr4wwnZ+S5qxttAW8aYqi+yAj1tozOYVsoi5xTES16EuqWPcZ5T4CaxZnypbYwxpioIuUCftXYPD36xulzviVEXN+gUbtbPyCGB5+V2Zmw5H17eCGws9b3pOQWc2sxOWTTGVH0hF+jN6lbj4tOa+t2+ac5GLtvxfzTNS2ZN3X5MT7obYutzcTnmeeWZzctfqDHGVLKQC/TOSXXonFSn7IaFuc7wtitfhhqJcM37dD5lKJ0rvkRjjAmKkAt0v2z/2dlXvj8ZzrgeLvoXVKsX7KqMMaZChVeg52XCd084Z7HUbQnXfwXtLgh2VcYYUynCJ9B/nQ3T7oLMndDrb9DvEYirEeyqjDGm0oR+oOccgG/+Aas+hsSOcNMsaNEj2FUZY0yl82twLhEZJCIbRSRZRB4s5nURkXHe11eJSLfAl1qEKqz5AsafBWs+h/Pvh1t+tDA3xkSsMrfQRSQamAAMAFKBxSIyVVXX+TQbDLT3/usJvOr9v2Jk7nYG09owHZp2hVFToImdv2KMiWz+7HLpASSr6hYAEfkYGA74Bvpw4F1VVWChiNQVkaaqujvgFW+aBZP/DO58GPAk9LoVokN/z5Exxpwof5IwCUjxeZzKsVvfxbVJAo4KdBEZDYwGaNmy9EvuS9SgHbQ4CwY/6/xsjDEG8G8fenF3dtDjaIOqTlTV7qravWHDhv7Ud6wG7Zzxyi3MjTHmKP4EeirQwudxc2DXcbQxxhhTgfwJ9MVAexFpIyJxwAhgapE2U4FR3rNdegEZFbL/3BhjTInK3Ieuqi4RuQ34FogGJqnqWhG5xfv6a8BMYAiQDOQAN1ZcycYYY4rj1+khqjoTJ7R9n3vN52cFbg1sacYYY8rDrwuLjDHGVH0W6MYYEyYs0I0xJkxYoBtjTJgQ53hmEGYskgZsP863JwL7AlhOKLA+Rwbrc2Q4kT63UtVir8wMWqCfCBFZoqrdg11HZbI+Rwbrc2SoqD7bLhdjjAkTFujGGBMmQjXQJwa7gCCwPkcG63NkqJA+h+Q+dGOMMccK1S10Y4wxRVigG2NMmKjSgV4lb05dwfzo83Xevq4SkZ9FpEsw6gyksvrs0+4sEXGLyJWVWV9F8KfPItJXRFaIyFoR+aGyaww0P9btOiIyTURWevsc0qO2isgkEflNRNaU8Hrg80tVq+Q/nKF6NwNtgThgJXBqkTZDgK9x7pjUC1gU7Loroc+9gXrenwdHQp992n2PM+rnlcGuuxI+57o49+1t6X3cKNh1V0KfHwKe8f7cEDgAxAW79hPo8/lAN2BNCa8HPL+q8hb6kZtTq2oBcPjm1L6O3JxaVRcCdUWkaWUXGkBl9llVf1bVg96HC3HuDhXK/PmcAW4HJgO/VWZxFcSfPo8EvlDVHQCqGur99qfPCtQSEQFq4gS6q3LLDBxVnYfTh5IEPL+qcqCXdOPp8rYJJeXtz004v+FDWZl9FpEk4DLgNcKDP59zB6CeiMwVkaUiMqrSqqsY/vR5PHAKzu0rVwN3qqqncsoLioDnl183uAiSgN2cOoT43R8RuQAn0M+t0Ioqnj99Hgs8oKpuZ+Mt5PnT5xjgTKA/UA1YICILVXVTRRdXQfzp80BgBdAPaAfMFpEfVTWzgmsLloDnV1UO9Ei8ObVf/RGR04E3gcGqur+Saqso/vS5O/CxN8wTgSEi4lLVryqlwsDzd93ep6rZQLaIzAO6AKEa6P70+UbgaXV2MCeLyFbgZOCXyimx0gU8v6ryLpdIvDl1mX0WkZbAF8D1Iby15qvMPqtqG1Vtraqtgc+Bv4VwmIN/6/YU4DwRiRGR6kBPYH0l1xlI/vR5B85fJIhIY6AjsKVSq6xcAc+vKruFrhF4c2o/+/wY0AB4xbvF6tIQHqnOzz6HFX/6rKrrReQbYBXgAd5U1WJPfwsFfn7O/wLeFpHVOLsjHlDVkB1WV0Q+AvoCiSKSCjwOxELF5Zdd+m+MMWGiKu9yMcYYUw4W6MYYEyYs0I0xJkxYoBtjTJiwQDfGmDBhgW6MMWHCAt0YY8LE/wPCWVojhfuXUAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "y = np.linspace(0,1,num=500)\n",
    "cumul = np.array([np.sum(x <= yi) for yi in y])/reps  # simulated CDF\n",
    "\n",
    "plt.plot(y, cumul)\n",
    "plt.plot([0,1],[0,1])\n",
    "plt.title('Simulated distribution of nominal $P$ value for $P$={}'.format(p))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The (simulated) distribution is not dominated by the uniform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:52: RuntimeWarning: invalid value encountered in double_scalars\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([5.78410e+04, 1.00000e+00, 2.34340e+04, 1.13167e+05, 4.68700e+03,\n",
       "        2.37444e+05, 1.13710e+05, 6.48320e+04, 1.24849e+05, 3.86530e+04]),\n",
       " array([0.        , 0.0999593 , 0.19991861, 0.29987791, 0.39983722,\n",
       "        0.49979652, 0.59975583, 0.69971513, 0.79967444, 0.89963374,\n",
       "        0.99959305]),\n",
       " <BarContainer object of 10 artists>)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD4CAYAAADy46FuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQn0lEQVR4nO3df6zd9V3H8edLuhF0gxUohLTM4kAdEIejFuLUMJtQxv6AJZB0mtEsTarIzEz8Y7A/xIw0gT8UQxQWHA0/ogPCNsFsDBtQ0YwBF8P4KXIdCBUC3UoYzoApe/vH+TQ7vdx+7un9dbi9z0dycs55f7+fz/l80ua8zvfnTVUhSdL+/My4ByBJenczKCRJXQaFJKnLoJAkdRkUkqSuFeMewHw7+uija+3ateMehiQtKY888sgPqmrVdMsOuqBYu3YtExMT4x6GJC0pSf5rf8vc9SRJ6jIoJEldBoUkqcugkCR1GRSSpC6DQpLUZVBIkroMCklSl0EhSeo66K7Mlt7N1l76zbF87vNXfnIsn6uDg1sUkqQug0KS1GVQSJK6DApJUpdBIUnqMigkSV0GhSSpy6CQJHUZFJKkLoNCktRlUEiSugwKSVKXQSFJ6jIoJEldBoUkqcugkCR1GRSSpC6DQpLUZVBIkroMCklSl0EhSeoyKCRJXQaFJKnLoJAkdRkUkqQug0KS1GVQSJK6DApJUteMQZHk+CT/mOTpJE8m+XyrH5lkR5Jn2/PKoTaXJZlM8kySjUP105M83pZdkyStfmiS21r9wSRrh9psbp/xbJLN8zp7SdKMRtmi2AP8cVV9GDgTuCTJycClwL1VdRJwb3tPW7YJOAU4B7g2ySGtr+uArcBJ7XFOq28BXquqE4GrgataX0cClwNnAOuBy4cDSZK08GYMiqp6uar+rb1+A3gaWA2cB9zUVrsJOL+9Pg+4tareqqrngElgfZLjgMOr6oGqKuDmKW329nUHsKFtbWwEdlTV7qp6DdjBT8NFkrQIDugYRdsl9KvAg8CxVfUyDMIEOKatthp4cajZzlZb3V5Pre/Tpqr2AK8DR3X6mjqurUkmkkzs2rXrQKYkSZrByEGR5H3A14A/qqof9Vadplad+mzb/LRQdX1VrauqdatWreoMTZJ0oEYKiiTvYRASf1NVX2/lV9ruJNrzq62+Ezh+qPka4KVWXzNNfZ82SVYARwC7O31JkhbJKGc9BbgBeLqq/nxo0V3A3rOQNgN3DtU3tTOZTmBw0PqhtnvqjSRntj4vmtJmb18XAPe14xj3AGcnWdkOYp/dapKkRbJihHU+BnwGeDzJo632ReBK4PYkW4AXgAsBqurJJLcDTzE4Y+qSqnq7tbsYuBE4DLi7PWAQRLckmWSwJbGp9bU7yRXAw229L1XV7tlNVZI0GzMGRVX9K9MfKwDYsJ8224Bt09QngFOnqb9JC5pplm0Hts80TknSwvDKbElSl0EhSeoyKCRJXQaFJKnLoJAkdRkUkqQug0KS1GVQSJK6DApJUpdBIUnqMigkSV0GhSSpy6CQJHUZFJKkLoNCktRlUEiSugwKSVKXQSFJ6jIoJEldBoUkqcugkCR1GRSSpC6DQpLUZVBIkroMCklSl0EhSeoyKCRJXQaFJKnLoJAkdRkUkqQug0KS1GVQSJK6DApJUpdBIUnqMigkSV0zBkWS7UleTfLEUO1Pk/x3kkfb49yhZZclmUzyTJKNQ/XTkzzell2TJK1+aJLbWv3BJGuH2mxO8mx7bJ63WUuSRjbKFsWNwDnT1K+uqtPa41sASU4GNgGntDbXJjmkrX8dsBU4qT329rkFeK2qTgSuBq5qfR0JXA6cAawHLk+y8oBnKEmakxmDoqruB3aP2N95wK1V9VZVPQdMAuuTHAccXlUPVFUBNwPnD7W5qb2+A9jQtjY2AjuqandVvQbsYPrAkiQtoLkco/hcksfarqm9v/RXAy8OrbOz1Va311Pr+7Spqj3A68BRnb7eIcnWJBNJJnbt2jWHKUmSppptUFwHfAg4DXgZ+LNWzzTrVqc+2zb7Fquur6p1VbVu1apVnWFLkg7UrIKiql6pqrer6ifAXzM4hgCDX/3HD626Bnip1ddMU9+nTZIVwBEMdnXtry9J0iKaVVC0Yw57fQrYe0bUXcCmdibTCQwOWj9UVS8DbyQ5sx1/uAi4c6jN3jOaLgDua8cx7gHOTrKy7do6u9UkSYtoxUwrJPkqcBZwdJKdDM5EOivJaQx2BT0P/B5AVT2Z5HbgKWAPcElVvd26upjBGVSHAXe3B8ANwC1JJhlsSWxqfe1OcgXwcFvvS1U16kF1SdI8mTEoqurT05Rv6Ky/Ddg2TX0COHWa+pvAhfvpazuwfaYxSpIWjldmS5K6DApJUpdBIUnqMigkSV0GhSSpy6CQJHUZFJKkLoNCktRlUEiSugwKSVKXQSFJ6jIoJEldBoUkqcugkCR1GRSSpC6DQpLUNeMfLpKkpWjtpd8c22c/f+Unx/bZC8EtCklSl0EhSeoyKCRJXQaFJKnLoJAkdRkUkqQug0KS1GVQSJK6DApJUpdXZmtsxnXl7MF21ewovEpZc+EWhSSpy6CQJHUZFJKkLoNCktRlUEiSugwKSVKXQSFJ6poxKJJsT/JqkieGakcm2ZHk2fa8cmjZZUkmkzyTZONQ/fQkj7dl1yRJqx+a5LZWfzDJ2qE2m9tnPJtk87zNWpI0slG2KG4EzplSuxS4t6pOAu5t70lyMrAJOKW1uTbJIa3NdcBW4KT22NvnFuC1qjoRuBq4qvV1JHA5cAawHrh8OJAkSYtjxqCoqvuB3VPK5wE3tdc3AecP1W+tqreq6jlgElif5Djg8Kp6oKoKuHlKm7193QFsaFsbG4EdVbW7ql4DdvDOwJIkLbDZHqM4tqpeBmjPx7T6auDFofV2ttrq9npqfZ82VbUHeB04qtOXJGkRzffB7ExTq059tm32/dBka5KJJBO7du0aaaCSpNHMNiheabuTaM+vtvpO4Pih9dYAL7X6mmnq+7RJsgI4gsGurv319Q5VdX1VrauqdatWrZrllCRJ05ltUNwF7D0LaTNw51B9UzuT6QQGB60farun3khyZjv+cNGUNnv7ugC4rx3HuAc4O8nKdhD77FaTJC2iGW8znuSrwFnA0Ul2MjgT6Urg9iRbgBeACwGq6skktwNPAXuAS6rq7dbVxQzOoDoMuLs9AG4AbkkyyWBLYlPra3eSK4CH23pfqqqpB9UlSQtsxqCoqk/vZ9GG/ay/Ddg2TX0COHWa+pu0oJlm2XZg+0xjlCQtHK/MliR1GRSSpC6DQpLUZVBIkroMCklSl0EhSeoyKCRJXQaFJKnLoJAkdRkUkqQug0KS1GVQSJK6DApJUpdBIUnqMigkSV0GhSSpy6CQJHUZFJKkLoNCktRlUEiSugwKSVKXQSFJ6jIoJEldBoUkqcugkCR1rRj3ACQd3NZe+s1xD0Fz5BaFJKnLoJAkdRkUkqQug0KS1GVQSJK6POtpinGdofH8lZ8cy+dK0kzcopAkdRkUkqQug0KS1DWnoEjyfJLHkzyaZKLVjkyyI8mz7Xnl0PqXJZlM8kySjUP101s/k0muSZJWPzTJba3+YJK1cxmvJOnAzccWxcer6rSqWtfeXwrcW1UnAfe29yQ5GdgEnAKcA1yb5JDW5jpgK3BSe5zT6luA16rqROBq4Kp5GK8k6QAsxK6n84Cb2uubgPOH6rdW1VtV9RwwCaxPchxweFU9UFUF3Dylzd6+7gA27N3akCQtjrkGRQH/kOSRJFtb7diqehmgPR/T6quBF4fa7my11e311Po+bapqD/A6cNTUQSTZmmQiycSuXbvmOCVJ0rC5Xkfxsap6KckxwI4k/95Zd7otgerUe232LVRdD1wPsG7duncslyTN3py2KKrqpfb8KvANYD3wStudRHt+ta2+Ezh+qPka4KVWXzNNfZ82SVYARwC75zJmSdKBmXVQJPm5JO/f+xo4G3gCuAvY3FbbDNzZXt8FbGpnMp3A4KD1Q2331BtJzmzHHy6a0mZvXxcA97XjGJKkRTKXXU/HAt9ox5ZXAH9bVd9O8jBwe5ItwAvAhQBV9WSS24GngD3AJVX1duvrYuBG4DDg7vYAuAG4Jckkgy2JTXMYryRpFmYdFFX1feAj09R/CGzYT5ttwLZp6hPAqdPU36QFjSRpPLwyW5LUZVBIkroMCklSl0EhSeryDxdJ0jw72P4AmlsUkqQug0KS1GVQSJK6DApJUpdBIUnqMigkSV0GhSSpy6CQJHUZFJKkLoNCktRlUEiSugwKSVKXNwXU2G5gJmlpcItCktRlUEiSugwKSVKXQSFJ6jIoJEldBoUkqcugkCR1GRSSpC6DQpLUZVBIkroMCklSl0EhSeoyKCRJXQaFJKnLoJAkdRkUkqQu/3CRlh3/UJN0YJbEFkWSc5I8k2QyyaXjHo8kLSfv+qBIcgjwV8AngJOBTyc5ebyjkqTl410fFMB6YLKqvl9V/wfcCpw35jFJ0rKxFI5RrAZeHHq/EzhjeIUkW4Gt7e3/JHlmDp93NPCDObSflVy12J+4j7HMeYyW23zBOS8LuWpOc/75/S1YCkGRaWq1z5uq64Hr5+XDkomqWjcffS0Vy23Oy22+4JyXi4Wa81LY9bQTOH7o/RrgpTGNRZKWnaUQFA8DJyU5Icl7gU3AXWMekyQtG+/6XU9VtSfJ54B7gEOA7VX15AJ+5Lzswlpiltucl9t8wTkvFwsy51TVzGtJkpatpbDrSZI0RgaFJKlrWQbFTLcEycA1bfljST46jnHOpxHm/Lttro8l+U6Sj4xjnPNp1Fu/JPm1JG8nuWAxx7cQRplzkrOSPJrkyST/vNhjnG8j/N8+IsnfJ/lem/NnxzHO+ZJke5JXkzyxn+Xz//1VVcvqweCA+H8CvwC8F/gecPKUdc4F7mZwDceZwIPjHvcizPnXgZXt9SeWw5yH1rsP+BZwwbjHvQj/zh8AngI+2N4fM+5xL8Kcvwhc1V6vAnYD7x332Ocw598CPgo8sZ/l8/79tRy3KEa5Jch5wM018F3gA0mOW+yBzqMZ51xV36mq19rb7zK4XmUpG/XWL38IfA14dTEHt0BGmfPvAF+vqhcAqmqpz3uUORfw/iQB3scgKPYs7jDnT1Xdz2AO+zPv31/LMSimuyXI6lmss5Qc6Hy2MPhFspTNOOckq4FPAV9exHEtpFH+nX8RWJnkn5I8kuSiRRvdwhhlzn8JfJjBhbqPA5+vqp8szvDGYt6/v97111EsgBlvCTLiOkvJyPNJ8nEGQfEbCzqihTfKnP8C+EJVvT34sbnkjTLnFcDpwAbgMOCBJN+tqv9Y6MEtkFHmvBF4FPht4EPAjiT/UlU/WuCxjcu8f38tx6AY5ZYgB9ttQ0aaT5JfAb4CfKKqfrhIY1soo8x5HXBrC4mjgXOT7Kmqv1uUEc6/Uf9v/6Cqfgz8OMn9wEeApRoUo8z5s8CVNdiBP5nkOeCXgYcWZ4iLbt6/v5bjrqdRbglyF3BRO3vgTOD1qnp5sQc6j2acc5IPAl8HPrOEf10Om3HOVXVCVa2tqrXAHcAfLOGQgNH+b98J/GaSFUl+lsGdmJ9e5HHOp1Hm/AKDLSiSHAv8EvD9RR3l4pr3769lt0VR+7klSJLfb8u/zOAMmHOBSeB/GfwiWbJGnPOfAEcB17Zf2HtqCd95c8Q5H1RGmXNVPZ3k28BjwE+Ar1TVtKdZLgUj/jtfAdyY5HEGu2W+UFVL9vbjSb4KnAUcnWQncDnwHli47y9v4SFJ6lqOu54kSQfAoJAkdRkUkqQug0KS1GVQSJK6DApJUpdBIUnq+n93wuus3oJ6iQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "reps = int(10**6)\n",
    "n = 50\n",
    "m = 100\n",
    "p = 0.01  # one of the worst offenders\n",
    "\n",
    "x = np.zeros(reps)\n",
    "for i in range(reps):\n",
    "    X = rng.binomial(n, p) # simulate X\n",
    "    Y = rng.binomial(m, p) # simulate Y\n",
    "    x[i] = 2*norm.cdf(absZ(n,m,X,Y)) - 1\n",
    "    \n",
    "plt.hist(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEJCAYAAACE39xMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAr4UlEQVR4nO3dd5wU9f3H8dfnGr13jyoKdkCPYkeQADZiYqIhamLMz5+JaOwtxm5iiQmx/7CEmMSQRE0oUsRCUBEFlI7oiZSj93Icd7e7n98f3z2yHLt3e3e7N7uzn+fjcQ9ud2ZnPrM3++a735n5jqgqxhhj0l+W1wUYY4xJDAt0Y4zxCQt0Y4zxCQt0Y4zxCQt0Y4zxCQt0Y4zxCQt0Y4zxibQPdBFZJiKDk7Ts1SJybjKWHWVd40Xk4drMm8j3IHJZid7+ZP6tqlhnbxH5XET2isgN9bnuKLUkZPvrc7+s7xpS6e+VjtIi0EXkDBGZIyK7RWSHiHwkIv0BVPV4VZ3lcX2efsDieQ/irTFR72e09Xn0t7odmKWqzVT1qXpe9yHqa/tFpJWIqIjsE5H9IrJBRG5M9noTJOF/r0S8HyLSWkT+JSLFIrJGREbXdl4RGSMi80WkVETG126rokv5QBeR5sAU4GmgNZAPPACUelmX34hIjtc1JEk3YJnXRdSzvsBWVW2qqo2BnwG/F5HO3pYVl1r/varYh/tS9/fjWaAM6AD8EHheRI6v5bwbgIeBV2qw/vioakr/AAXAriqmrwbOjfj9NmAxUAy8HH5TpwF7gXeAVhGvVeCoiMfjgYdjLPtO4OvwcpYDF4ef/zMQAkqAfcDt4eePAN4AtgLfADdUqrsf8Fl4eX8HJkSuuybzVqrzDmB9eN6VwNAqalwdnn8x7j/InCjv513h7d0J/BFoWN37V836KpZ9LDAL2IX7AF8U5e96a7i23eHtbhjj/Ym6LOA9IAgcCNfRq4p9KOq6qqqTmu9vB7c/jvVG3d+iLSfK9twEvBXxuGP4b3VkpfnuBF6v9NwfgKdqUkOs/SDicZWfhYj5Dvt7xbmfHLIP1/b9qOL9bIIL6F4Rz/0ZeLSO8z4MjK9pJlZZayIXlowfoDmwHfgTMDLyAxJlx1oNzA1/qPKBLbgg7Ac0CO8w90W8tiaB/r3wjpkFXIr7AHeqPF/4cRawALgXyAOOBFYBw8PT84A14R0tF7gEKCdKoMczb8X6gd7AOuCI8PPdgZ7Raox4biHQBWgU4/1cGp7eGvio0npjvn9VrO/c8HYUAneHt28ILjR6V5r30/B73hpYAVwb5f2pclm4MPhpNftY1HXFsezV1Gx/O+Q9qWobqcH+FmV7Xq1YL9AS9x/NfEAqzdcN2A80Dz/OBjYCg2pSQzX7QZWfhSi1H/x71WA/WUjEPlyT9wP37X9XjJ8p4df0A0oqLfNWYHKUddVk3oQHesp3uajqHuAM3E7zIrBVRCaJSIcYL3laVTer6nrgA+ATVf1cVUuBf+He8NrU8U9V3aCqIVX9O/AVMCDG7P2Bdqr6oKqWqeqqcO2XhacPwu2sY1W1XFVfB+bFWFZN5g3iguQ4EclV1dWq+nU1m/aUqq5T1ZIY058JT98BPAL8oJrlxWMQ0BTXailT1fdwH6zKy34q/J7vACbjvjrXdlnVibaueJZd1/0t6jbWcH+rrC9wm4jswIWpAhdqOEUqqOoa3H9A3w4/NQTYr6pzE1BDheo+C1WpyX5S1T7clxjvh6peoKotY/xcEH59U9w3qEi7gWZR1lWTeRMuLfpNVXUF8GMAETkG+Aswlugf2s0Rv5dEedy0NjWIyJXAzbhWL+HltI0xezfgCBHZFfFcNu4DD67Vs77SB2xNjGXFPa+qFoYP9twPHC8iM4CbVXVDjGWDa9FXJXL6mnA9dXUEsE5VQ5WWnV9pvk0Rv++Pse54l1WdaOuKZ9l13d+ibmMN97eDRKQBrpuih6oWVTc/8Bruc/QqMDr8uGJZtaqhkuo+C1WJ928bcx+uxfsRzT5cT0Gk5rhvC3WZN+FSvoVemap+gftKd0ICFrcfaBzxuGO0mUSkG65VMQZoo6otcV0RUlFWpZesA76p9L99M1U9Lzx9I5AvIhLxmq4xaqzJvKjqa6p6Bu6DpMBjMWqkmucrdKm03sj/HKp6/6pa7gagi4hE7n9dcX3/NZXIZdXnsmOKY3+ryglAcQ3C65/A4PABwosJB3oNa6hqP6jus1CVeN//qva1Kt8PEZkWPvsl2s+08GxfAjkicnTES/sQ/eBtTeZNuJQPdBE5RkRuqTgiLSJdcC2KuQlY/EJgtIhki8gI4OwY8zXB7TRbwzVcxaH/oWzG9Q1W+BTYIyJ3iEij8PJPqDjVEvgYCAA3iEiOiHyH2F9l4543fA7vkHCr5ACuhRiMUWO8rhORziLSGteX+feIaQuJ/f5Vtb5PcP2xt4tIbvjc7AtxB3trKpHLqs9lV6W6/a0q/ahBeKjqVly/9R9xwbuiFjUsJPZ+UN1noSqJeP+rfD9UdaS6s1+i/YwMz1MMvAk8KCJNROR0YBTuYGfl5VU7b/hz3BD3TSVbRBom6iyzlA903FeVgcAnIlKMC/KlwC0JWPYvcDvILtzpRf+ONpOqLgeexIXrZuBE3AHCCr8B7hGRXSJyq6oGw8vtizuqvw14CWgRXl4Z8B1cN9JO3AGnN2OsO+55cf3nj4bXtwlojwvhw2qM8fpoXgPexh3IWoU7kFOhqvcv5vrC23QR7iD3NuA54Mrwt68aSeSy6nPZ1ay3uv2tKn1xn4+aeA13sPpgd0sNa4i5H1T3WahKgt7/vtT8/Yjm50Aj3IHvvwE/U9VlcLCVf3c884bdg2ts3QlcHv79ngTU6I56G2OMSX/p0EI3xhgTBwt0Y4zxCQt0Y4zxCQt0Y4zxCc8uLGrbtq12797dq9UbY0xaWrBgwTZVbRdtmmeB3r17d+bPn+/V6o0xJi2JSKyryq3LxRhj/MIC3RhjfMIC3RhjfMIC3RhjfMIC3RhjfKLaQBeRV0Rki4hEHeBGnKdEpFBEFovIyYkv0xhjTHXiaaGPB0ZUMX0kcHT45xrg+bqXZYwxpqaqDXRVnQ3sqGKWUcCr4ds5zQVaikinRBVojDG+ESyHD34H6xckZfGJ6EPP59BbQBUR4/ZfInKNiMwXkflbt25NwKqNMSZNbFwELw6Bdx+A5ZOSsopEXCka7ZZUUQdZV9VxwDiAgoICG4jdGON/5Qdg9uPw4Vho3Aa+/yocNyopq0pEoBdx6H0nO3PofSeNMSYzrZ0LE8fA9q+g7+Uw/GFo1Cppq0tEl8sk4Mrw2S6DgN2qujEByzXGmPRUuhem3gavjIBAKVz+Jnz72aSGOcTRQheRvwGDgbYiUgTcB+QCqOoLwFTgPKAQd/fvq5JVrDHGpLzCd2DyjbC7CAb+Lwz5FTRoWi+rrjbQVfUH1UxX4LqEVWSMMelo/w6Y8UtY9Bq07QU/mQ5dB9VrCZ4Nn2uMMb6xfCK8dSvs3w5n3gpn3Qa5Deu9DAt0Y4yprb2bYOqtsGIydOoDl78BnU7yrBwLdGOMqSlVWPhXmHG3Oy3x3Pvh1Osh29tItUA3xpia2LkGJv8CVr0PXU+Di56Gtkd5XRVggW6MMfEJBeHTF+HdB0EEzvstFFwNWakzaK0FujHGVGfrSph0Paz7BI46Fy4YCy27VPuy+maBbowxsQTL4aOx8J/HIa8JXDwOTvq+a6GnIAt0Y4yJZsPnMPF62LwEjr8YRj4BTdt5XVWVLNCNMSZSeQnMehTmPA1N2sGlf4VjL/C6qrhYoBtjTIXVH7m+8h1fQ78r4FsPQ6OWXlcVNwt0Y4w5sMeNUz7vJWjZDa6cCEcO9rqqGrNAN8Zktq9musG09qyHQT+HIfe4A6BpyALdGJOZ9u+A6XfB4gnQ7hi4eiZ06e91VXVigW6MySyqsOxfbrzyA7vg7DvgzFsgp4HXldWZBboxJnPs2Qhv3QIr34Ij+sFFE6HjCV5XlTAW6MYY/1OFz/8MM+6BYCkMe8j1l3s8mFai+WtrjDGmsh3fwOQb4JvZ0O0MuOgpaNPT66qSwgLdGONPoSB88n/w3kMg2XDB7+HkH6fUYFqJZoFujPGfLStg4hhYPx+OHu7CvEW+11UlnQW6McY/AmXw4e9h9hPQsDl892U44bspO5hWolmgG2P8Yf0CN5jWlmVwwiUw8jFo0tbrquqVBboxJr2V7YdZv4aPn4WmHeEHE6D3SK+r8oQFujEmfX3zgTuDZccqOOXHMOxBaNjC66o8Y4FujEk/B3bDzPtgwR+hVQ/40WTocZbXVXnOAt0Yk15WTocpN8G+TXDa9TD4bshr7HVVKcEC3RiTHoq3wbQ7YOnr0P44uPQv0PkUr6tKKRboxpjUpgpL34Bpt7txywffDWfcBDl5XleWcizQjTGpa/d6eOtm+HI65J8CFz0DHY7zuqqUZYFujEk9oRB89ieYeS8Ey2H4r2HgtZCV7XVlKc0C3RiTWrZ/DZN/Aas/cGeuXPgHaH2k11WlBQt0Y0xqCAZg7nPw/iOQnQcXPgUnX5kxl+0nQlzDjonICBFZKSKFInJnlOktRGSyiCwSkWUiclXiSzXG+NbmZfDyMJj5K+g5BK77BE75kYV5DVXbQheRbOBZYBhQBMwTkUmqujxituuA5ap6oYi0A1aKyF9VtSwpVRtj/CFQCh886X4atoRLXoHjv2NBXkvxdLkMAApVdRWAiEwARgGRga5AMxERoCmwAwgkuFZjjJ8UzXdD3G5dASddCsN/A03aeF1VWosn0POBdRGPi4CBleZ5BpgEbACaAZeqaqjygkTkGuAagK5du9amXmNMuisrhvcecf3lzY+A0f+AXsO9rsoX4gn0aN99tNLj4cBCYAjQE5gpIh+o6p5DXqQ6DhgHUFBQUHkZxhi/W/UfN5jWztVQcDWce78bt9wkRDyBXgR0iXjcGdcSj3QV8KiqKlAoIt8AxwCfJqRKY0x6K9nlDnh+9iq07gk/ngrdT/e6Kt+JJ9DnAUeLSA9gPXAZMLrSPGuBocAHItIB6A2sSmShxpg09cVbMOVmKN4Cp/8CBt8FuY28rqrefPz1dv4+by0l5UFKykOUBYJc1Cef0QMT3+1cbaCrakBExgAzgGzgFVVdJiLXhqe/ADwEjBeRJbgumjtUdVvCqzXGpI99W934K8vehA4nwA/+Bvkne11VQk1etIF3V2wmqBAKKYFQiGAIQqoEQ+7n87U7yc3JokOzhjTMy6ZBdhZ6WK91YsR1YZGqTgWmVnruhYjfNwDfSmxpxpi0pAqL/wHT73AHQM+5B864EbJzva4s4Z5+7yvW7yyhQ/OGZGUJ2SLu3yzIzsoiW2DgkW144KLj6dI6+UP82pWixpjE2V3kxir/6m3o3N8NptX+GK+rSpqNuw9wySmdeWDUCV6XAligG2MSIRSCBa+4uwhpCEY8BgP+x9eDae0rDbD3QIBOLVPneIAFujGmbrYVwqTrYe0cOHKwG0yrVXevq0q6TbtLAOjUoqHHlfyXBboxpnaCAfj4GZj1G8hpAKOehb4/9M1l+yVlQbbtKz34WCsdx1xctBuAjs0t0I0x6WzTEph4HWxcBMdcAOc/Cc06el1VXLbvK6WkPEio4mwUVUKhin/dcyXlQa5/7XM27TlQ7fI618PBznhZoBtj4hcohdlPwIe/h0at4Ht/guNGedIq/9una5m0sPI1jlUrLgscbFlXJy87i4dGHU+jvENjMnJL2zZrQL71oRtj0s7aT1xf+baV0Gc0DH8EGrf2pJRAMMQTM1bSICeLLq3ibyE3zMnmZ4N70qNNk4OnF2aJkCVCdpaEf4fsLKFnu6Z0b9skiVuReBboxpiqle6D9x6CT/4PWnSGy9+Ao871tKQLnv6QHcVlvHD5KYw4IT26euqDBboxJrav33O3g9u1FgZcA0PvhQbNvK6KVVuLaZKXzdBj23tdSkqxQDfGHK5kJ8y4Bxb+BdocDVdNh26nel3VQUFV/uf0HuRmx3XTtYxhgW6MOdSKyfDWLVC8Dc64Gc6+A3JT59S8ijFS8rL9e9FSbVmgG2OcvZth2m2wfCJ0PNHdeOKIvl5XdZiygLt3Tl6Otc4rs0A3JtOpwqK/wfS7oLzE9ZOfdkPKDqZlgR6bBboxmWzXWph8I3z9LnQZBBc9De16eV1VlUoDQcACPRoLdGMyUSgE816Cd+53j0c+Af1/ClmpH5Kl4RZ6AzsgehgLdGMyzbavYOIYWDcXeg6FC8dCy/S5aXtZ0LpcYrFANyZTBMthzlMw6zF3C7hvPw99fpB2g2lZH3psFujGZIKNi9xgWpuWuLFXRj4BzTp4XVWtVAR6Awv0w1igG+Nn5QfgP4/CR09Bk7bw/T/DcRd5XVWdWJdLbBboxvjVmo9h0hjYXgh9L4fhD7sREtPcwS4XOyh6GAt0Y/ymdC+88wDMe9Ed7LziX9BziNdVJYz1ocdmgW6MnxS+484r310EA6+FIb+CBk29riqhSi3QY7JAN8YP9u+AGXe7Kz7b9oKfzICuA72uKikq+tDtoOjhLNCNSWeqbuyVqbe6ERLPvBXOui2lBtNKtP/2odvgXJVZoBuTrvZucqMifjEFOvWBy9+ETid5XVXSWR96bBboxqQbVVj4V9fFEiiFcx+AU8dAdmZ8nMtsLJeYMmMPMMYvdq52dxBaNQu6nuYG02p7lNdV1Ss7Dz02C3Rj0kEoCJ++CO8+AJIF5z8Jp/wkLQbTSrTScjsPPRYLdGNS3ZYvYNL1UPQpHDUMLvg9tOzidVWeqWih52an1xg09cEC3ZhUFSyHD8fC7MchrylcPA5O+n7aDaZVIRhS3vysiH2lgYPPqcKCNTtZvH7XwceR/0bS8JN7DgRokJOFpOn7kExxBbqIjAD+AGQDL6nqo1HmGQyMBXKBbap6dsKqNCbTbPjcDXG7eSkc/x0Y+Tg0bed1VXWycN1Obnt98WHP52VnMeSY9jTOC5+GKBX//DewRQ6ZxLGdmiex0vRVbaCLSDbwLDAMKALmicgkVV0eMU9L4DlghKquFZH2SarXGH8rL4FZv4E5T0OT9nDZa3DM+V5XlRB7SlzL/M9XD+DE/BYHn2+Qk02jPDunPBHiaaEPAApVdRWAiEwARgHLI+YZDbypqmsBVHVLogs1xvdWf+T6ynd8DSdfCcMegkYtva4qYSq6Wjo2b0jLxnkeV+NP8RwmzgfWRTwuCj8XqRfQSkRmicgCEbky2oJE5BoRmS8i87du3Vq7io3xmwN7YMrNMP48CAXgyonudEQfhTnA/jIX6I0b2KG7ZInnnY125KHyIYsc4BRgKNAI+FhE5qrql4e8SHUcMA6goKAgymEPYzLMl2/DlBthzwYYdB0M+SXkNfG6qqQoLnUXBDWx7pWkiSfQi4DIc6Q6AxuizLNNVYuBYhGZDfQBvsQYc7ji7TD9TljyD2h3DFw9E7r097qqpCoOd7k0zrMWerLE0+UyDzhaRHqISB5wGTCp0jwTgTNFJEdEGgMDgRWJLdUYH1CFpW/AswNg2Ztw9p3wv7N9H+YAxWVB8rKz7ArPJKr2v0pVDYjIGGAG7rTFV1R1mYhcG57+gqquEJHpwGIghDu1cWkyCzcm7ezZCG/dDCunwhH9YNQk6HC811XVm/1lARo3sO6WZIrru4+qTgWmVnruhUqPnwCeSFxpxviEKnz2Krz9KwiWwrcehoE/y5jBtCoUlwZpYt0tSWXvrjHJtGOVG0zrm9nQ7Qy46Clo09PrqjxRXBqgibXQk8oC3ZhkCAVh7vPw3sOQlQMXjIWTf5SRg2lVKC4L2AHRJLN315hE27wcJo2B9Qug1wg4/3fQovKlG5lnf1nQWuhJZoFuTKIEyuDD38Hs30LD5vDdl+GE76bcYFqhkPK7mV+ydW/pwec04tKSygNjRT48fFrsmStfaFK4ZR8DerSuecEmbhboxiTC+gVuMK0ty+HE78GIR6FJW6+riqpoZwnPvF9I84Y5h3SBVP5/p/J/Q1WNbhg56fDluCdaNMrlrF7pPcBYqrNAN6YuyvbD+4/A3OegaUf4wQToPdLrqqpUGr6F28MXn8hFfY7wuBqTSBboxtTWN7Nh0g2w8xs45SoY9gA0bFH96zxWGrA7/viVBboxNXVgN8y8FxaMh1Y94EeTocdZXlcVt4o7/jSwKzZ9xwLdmJpYOQ2m3AT7NsNp18PguyGvsddV1UhZwG6y7FcW6MbEo3gbTLsDlr4O7Y+Hy/4K+ad4XVWtVAS6tdD9xwLdmKqowpLXYdrtULrXtcjPuAly0vcGDdZC9y8LdGNi2b3eDab15XTIL4BRz0D7Y72uqs4q+tAt0P3HAt2YykIh+Gw8vH2vu4PQ8F/DwGshyx9XOZbZWS6+ZYFuTKTtX7tTEdd86M5cufApaN3D66oSyrpc/MsC3RiAYMBdHPT+I5DdwN3Ts98VKXfZfiKUWpeLb1mgG7NpqRtMa8Pn0Pt8OP9JaN7J66qS5uBZLtn+6EIy/2WBbjJXoBQ+eNL9NGwJl/wRjr/Yl63ySNbl4l8W6CYzrZvnWuVbv4CTLnWDaTXOjJEALdD9ywLdZJayYnfTibnPQ/MjYPQ/ode3vK6qXpUFg2RnCdlZ/v4mkoks0E3mWDXLncGyaw30/ykMvc+NW55hygIhO2XRpyzQjf+V7IK374HP/wyte8KPp0L3072uyjNlgZB1t/iUBbrxty/egik3Q/FWOP1GGHwn5DbyuipPlQUt0P3KAt34074tbvyVZf+CDifC6AlwRD+vq0oJpeXW5eJXFujGX1Rh8d9h+p3uAOiQe1zLPDvX68pSRmkwZCMt+pQFuvGPXevcWOWFM6HzADeYVrveXleVcqwP3b8s0E36C4Vg/svwzv2gIRjxGAz4n6QNplVcGmD7vrKkLDuWdTv3s2T97oQs6+ut+2jawD76fmR/VZPethXCpOth7Rw48hy48A/QqlvSVqeqjPzDB6zdsT9p66gPF/fL97oEkwQW6CY9BQPw8dPw/m8gtyGMeg76jk76ZftFO0tYu2M/l/XvQv/u9XdlaYPcLM44qi0NchLzraNhrnW5+JEFukk/Gxe7y/Y3LoJjLnCDaTXrWC+rfvnDbwC4fFA3TshvUS/rNCZeFugmfZQfgNmPw4djoXEb+P6rcNyoelv9pt0HGD9nNQC9Ozart/UaEy8LdJMe1n7iWuXbvoQ+o2H4I/U+mNa+0gAAv7rgOHLtPG6TguLaK0VkhIisFJFCEbmzivn6i0hQRC5JXIkmo5Xug6m3wyvDobwELn8DLn7ek5ERQ6oAdGzesN7XbUw8qm2hi0g28CwwDCgC5onIJFVdHmW+x4AZySjUZKDCd2HyjbB7nTsNcei90MC7ro5A0AW6Nc5Nqoqny2UAUKiqqwBEZAIwClheab7rgTeA/gmt0GSekp0w45ew8K/Q5mi4ahp0O9Xrqg620LOzLNFNaoon0POBdRGPi4CBkTOISD5wMTCEKgJdRK4BrgHo2rVrTWs1mWD5JJh6KxRvgzNuhrPvcKclpoBAyFroJrXFE+jRTuzVSo/HAneoalCqOA9YVccB4wAKCgoqL8Nksr2bXZCvmAQdT4Qf/hM69fG6qkMEQ9ZCN6ktnkAvArpEPO4MbKg0TwEwIRzmbYHzRCSgqv9ORJHGx1Rh4Wsw42530HPofXDa9Sk5mNbBQPf5PUdN+oon0OcBR4tID2A9cBkwOnIGVe1R8buIjAemWJibau1cA1NuhK/fg66nwkVPQ9ujva4qpv+20C3QTWqqNtBVNSAiY3Bnr2QDr6jqMhG5Njz9hSTXaPwmFIJ5L8I7D7hL9c/7LRRcDSnelWGBblJdXBcWqepUYGql56IGuar+uO5lGd/a+qUbTGvdXOg5FC4cCy3T4wB5UC3QTWqzK0VN/QiWw0d/gP88BrmN4dsvQJ/Lkj6YViIFQyHAAt2kLgt0k3wbFrrL9jctcWOvnPdbaNre66pqLOjynBwLdJOiLNBN8pSXuBb5R09Bk7Zw6V/g2Au9rqrWKlroWWn0rcJkFgt0kxxrPnat8u2F0O9y+NbD0KiV11XVycEWerYFuklNFugmsUr3urNX5r3oDnZe8W/oeY7XVSVEwPrQTYqzQDeJ89VMN5jWnvUw8Gcw5B5o0NTrqhLm4Fgu1uViUpQFuqm7/Ttg+l2weAK07Q1Xvw1dBnhdVcL9d7RFC3STmizQTe2pwvJ/w9Tb3AiJZ93mfnIaeF1ZUoTsPHST4izQTe3s3QRv3QJfTIFOfeGKf7lBtXysYrRFO23RpCoLdFMzqvD5X9x45cFSGPYgDLoOsv2/K1Vc+p9lgW5SlP8/hSZxdq6Gyb+AVbOg2+lw4VPQ9iivq6o3QWuhmxRngW6qFwrCp+Pg3QdBsuH838EpV6X8YFqJZi10k+os0E3VtnzhLhAqmgdHDXODabXo7HVVnrAWukl1FugmukAZfDQWZj8BeU3hOy/Cid9Lq8G0Eq3ioKhd+m9SlQW6Odz6z9wQt5uXwgnfhRGPQdN2XlfluZC10E2Ks0A3/1VeAu//Gj5+Bpp2gMv+Bsec53VVKSNgN7gwKc4C3TirP3St8h2r4OQfudMRG7X0uqqUElIlS6CqG6Eb4yUL9Ex3YA+8cx/MfwVadYcrJ8GRZ3tdVUoKhNRa5yalWaBnsi9nwJSbYO9GOHUMnHM35DXxuqqUFbJANynOAj0TFW+H6XfCkn9Au2Ph+69C5wKvq0p5gZCSk2Hn3pv0YoGeSVRh6Rsw7XbX1XL2nXDmLZCT53VlaSEYcn3oxqQqC/RMsWeDG0xr5VQ44mQY9Qx0ON7rqtJKMKTkZFsL3aQuC3S/U4XP/gRv/wqC5e5WcIN+DlnZXleWdoKqdlGRSWkW6H62YxVMugFWfwDdz4QL/wBtenpdVdoKBtUuKjIpzQLdj0JBmPs8vPcwZOfCBWPdueV2QK9OgmpnuZjUZoHuN5uXu8G01i+AXiPcyIgt8r2uyheCdtqiSXEW6H4RKIMPfwezfwsNm8N3X3bjsFifb8LYhUUm1Vmg+0HRAtcq37LcjYg44jFo0sbrqnzHLiwyqc4CPZ2V7Yf3H4G5z0HTjvCDv0PvEV5X5VuBUIhs+8ZjUpgFerr6ZrYbTGvnanf3oGEPQMMWXlfla8GQjbRoUltcpz2IyAgRWSkihSJyZ5TpPxSRxeGfOSLSJ/GlGgAO7HanIv7pQkDgR1PcXYQszJMuGApZoJuUVm0LXUSygWeBYUARME9EJqnq8ojZvgHOVtWdIjISGAcMTEbBGW3lNDeY1r7NcNoNMPguyGvsdVUZI6jWQjepLZ4ulwFAoaquAhCRCcAo4GCgq+qciPnnApl508lkKd7mxl9Z+ga0Px4uew3yT/a6qrS3o7iMDbtK4p5/d0k5uRboJoXFE+j5wLqIx0VU3fq+GpgWbYKIXANcA9C1a9c4S8xgqrDknzDtDijdC+f8Ek6/0ZeDae0rDbBo3S5UQdGDz2v4V42YV8NPRj7HwfmivFYPm43SQJC73lzC3gOBGtV5di+7FZ9JXfEEerQmiUZ5DhE5BxfoZ0SbrqrjcN0xFBQURF2GCdtdBFNuhq9mQH6BG0yr/bFeV5U0v52xkvFzVtfrOps3zOGZ0f1okBP/uDYn5DdPYkXG1E08gV4EdIl43BnYUHkmETkJeAkYqarbE1NeBgqFYMEfYeZ9oEEY/hsY+L++H0xr5/4yOjZvyDOj+x187tAzBOWw5yMnR94WTirN5547/PX5LRvRqon/vu2YzBVPoM8DjhaRHsB64DJgdOQMItIVeBO4QlW/THiVmWL71+4MljUfQo+z3WBarXt4XVW9KAuEaNYwh4Lurb0uxZi0VW2gq2pARMYAM4Bs4BVVXSYi14anvwDcC7QBngu3lAKqarfAiVcwAHOfhfd/DdkN4KJnoN/lGXXZflkgRF6ODR5mTF3EdWGRqk4FplZ67oWI338K/DSxpWWITUtg4hjYuBB6nw/nPwnNO3ldVb0rC1qgG1NXdqWoVwKlMPsJ+PD30KgVfG88HPftjGqVRyoNhMizuwEZUycW6F5Y96lrlW9bCSddBiN+A40zu++4og/dGFN79gmqT2XF8O5D8MkL0Dwffvg6HD3M66pSQlkgRAPrcjGmTizQ68vX78PkG2DXWuj/Uxh6nxu33ADWh25MIligJ1vJLnj7l/D5X6B1T7hqGnQ7zeuqUk6Z9aEbU2cW6Mm0Ygq8dQsUb4UzboKz74DcRl5XlZLstEVj6s4CPRn2bYGpt8Hyf0OHE2H0BDiiX7Uvy2TW5WJM3VmgJ5IqLJoA0++E8v0w5Fdw+i8gO9frylKe63Lx9/AGxiSbBXqi7FoHU26Ewneg8wA3mFa73l5XlTasy8WYurNAr6tQCOa/DO/c71roIx93Z7H4fDCtRFJV63IxJgEs0Oti21fuvp5rP4Yjz3GDabXq5nVVaacsGAKw89CNqSML9NoIlsOcp2HWo5DbEEY9B31HZ+xl+3VVFnCBbqctGlM3Fug1tXGRu2x/02I49kI470lo1sHrqtLawUC3FroxdWKBHq/yAzD7cfhwLDRuA99/FY4b5XVVvlDR5WKBbkzdWKDHY+1c1yrf/hX0GQ3DH8n4wbQSqaKFnmtdLsbUiQV6VUr3wbsPwqfjoEUXuPwNOOpcr6vyHetyMSYxLNBjKXwHJt8Eu9fBgGtg6L3QoKnXVflSqR0UNSYhLNAr278DZvwSFr0GbY6Gn0yHroO8rirprh4/jyXrd3uy7vKK0xZzLdCNqQsL9EjLJ8Jbt8L+7XDmLXDW7e60RJ/bXVLOu19s4ZRurejVwZtvIY1ycyjo1sqTdRvjFxboAHs3wdRbYcVk6HiS6yvvdJLXVdWbwi17Afj54J4MPdZOwTQmXWV2oKvCwtdgxl3utMRz74dTx/huMK1gSFlUtItAUKNOn7VyCwC9OjSrz7KMMQmWdoG+dW8pX2zac9jzGj2riPE0Dfauo9e8e2i96SN2tStg5YBHKGl+JFq4s0YL0igTYtZSwxo1xgtizx/9+cmLNvDWko0xXuW0aJRLfksbq92YdJZ2gf7pNzu47rXPav36LEJckT2T23MmoAj3BK7ir+uGouu2A9sTV2iK+fFp3Rl2XOzulPyWjcjKsqELjElnaRfop/Zsw+vXnhp1WuyhVNyEhrsK6fbRHTTdsoDd+Wez7vRH+E7Tznwn6tzRlh99Suz5Y1Vz+ISaDgNTk2Xn5WRxVHs75dIYv0u7QG/dJI/WTWp4lWawHD4aC/95HPKawMX/R4uTLqWFDaZljPGRtAv0Gtuw0F22v3kJHPdtOO8JaNre66qMMSbh/Bvo5SVueNs5T0OTtnDpX9zoiMYY41P+DPQ1c9yNJ7YXQr8r4FsPQSO7aMUY42/+CvQDe+DdB2DeS9CyK1zxb+h5jtdVGWNMvfBPoH81EybfCHvWw6Cfw5B73AFQY4zJEOkf6Pt3wPS7YPEEaNsbrn4bugzwuipjjKl3cQ1vJyIjRGSliBSKyJ1RpouIPBWevlhETk58qZWowtI34Zn+sPR1N5DWtR9YmBtjMla1LXQRyQaeBYYBRcA8EZmkqssjZhsJHB3+GQg8H/43OfZsdINpfTEFOvWFKydCxxOStjpjjEkH8XS5DAAKVXUVgIhMAEYBkYE+CnhV3eAjc0WkpYh0UtWqBxCpjS/fhjd+CsFSGPYgDLoOstO/58gYY+oqniTMB9ZFPC7i8NZ3tHnygUMCXUSuAa4B6Nq1a01rddr0hC79YeTj7ndjjDFAfH3o0a6PrzyuXzzzoKrjVLVAVQvatWsXT32Ha9PTjVduYW6MMYeIJ9CLgC4RjzsDG2oxjzHGmCSKJ9DnAUeLSA8RyQMuAyZVmmcScGX4bJdBwO6k9J8bY4yJqdo+dFUNiMgYYAaQDbyiqstE5Nrw9BeAqcB5QCGwH7gqeSUbY4yJJq7TQ1R1Ki60I597IeJ3Ba5LbGnGGGNqIq4Li4wxxqQ+C3RjjPEJC3RjjPEJC3RjjPEJccczPVixyFZgTS1f3hbYlsBy0oFtc2awbc4Mddnmbqoa9cpMzwK9LkRkvqoWeF1HfbJtzgy2zZkhWdtsXS7GGOMTFujGGOMT6Rro47wuwAO2zZnBtjkzJGWb07IP3RhjzOHStYVujDGmEgt0Y4zxiZQO9JS8OXWSxbHNPwxv62IRmSMifbyoM5Gq2+aI+fqLSFBELqnP+pIhnm0WkcEislBElonIf+q7xkSLY99uISKTRWRReJvTetRWEXlFRLaIyNIY0xOfX6qakj+4oXq/Bo4E8oBFwHGV5jkPmIa7Y9Ig4BOv666HbT4NaBX+fWQmbHPEfO/hRv28xOu66+Hv3BJ3396u4cftva67Hrb5buCx8O/tgB1Ante112GbzwJOBpbGmJ7w/ErlFvrBm1OrahlQcXPqSAdvTq2qc4GWItKpvgtNoGq3WVXnqOrO8MO5uLtDpbN4/s4A1wNvAFvqs7gkiWebRwNvqupaAFVN9+2OZ5sVaCYiAjTFBXqgfstMHFWdjduGWBKeX6kc6LFuPF3TedJJTbfnatz/8Oms2m0WkXzgYuAF/CGev3MvoJWIzBKRBSJyZb1VlxzxbPMzwLG421cuAX6hqqH6Kc8TCc+vuG5w4ZGE3Zw6jcS9PSJyDi7Qz0hqRckXzzaPBe5Q1aBrvKW9eLY5BzgFGAo0Aj4Wkbmq+mWyi0uSeLZ5OLAQGAL0BGaKyAequifJtXkl4fmVyoGeiTenjmt7ROQk4CVgpKpur6fakiWebS4AJoTDvC1wnogEVPXf9VJh4sW7b29T1WKgWERmA32AdA30eLb5KuBRdR3MhSLyDXAM8Gn9lFjvEp5fqdzlkok3p652m0WkK/AmcEUat9YiVbvNqtpDVburanfgdeDnaRzmEN++PRE4U0RyRKQxMBBYUc91JlI827wW940EEekA9AZW1WuV9Svh+ZWyLXTNwJtTx7nN9wJtgOfCLdaApvFIdXFus6/Es82qukJEpgOLgRDwkqpGPf0tHcT5d34IGC8iS3DdEXeoatoOqysifwMGA21FpAi4D8iF5OWXXfpvjDE+kcpdLsYYY2rAAt0YY3zCAt0YY3zCAt0YY3zCAt0YY3zCAt0YY3zCAt0YY3zi/wHYIleL3utCZAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "y = np.linspace(0,1,num=500)\n",
    "cumul = np.array([np.sum(x <= yi) for yi in y])/reps  # simulated CDF\n",
    "\n",
    "plt.plot(y, cumul)\n",
    "plt.plot([0,1],[0,1])\n",
    "plt.title('Simulated distribution of nominal $P$ value for $P$={}'.format(p))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, the (simulated) distribution is not dominated by the uniform. The _nominal_ $P$-value based on treating the $Z$-statistic as if it had a standard normal distribution is not an _actual_ $P$-value for this problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### An exact conditional test based on invariance: permutation methods\n",
    "\n",
    "We now explore a different approach that yields an exact test rather than approximate test.\n",
    "\n",
    "Let $N := n+m$.\n",
    "Let $I_1, \\ldots, I_n$ be the values of the draws from the first population and $I_{n+1}, \\ldots, I_N$ be\n",
    "the values of the draws from the second population. \n",
    "Because the draws are all independent, $\\{I_j\\}$ are independent random variables with a Bernoulli distribution. \n",
    "If $p_x = p_y$, they are identically distributed also.\n",
    "\n",
    "Define the random vector $I := (I_j)_{j=1}^N$.\n",
    "Let $\\pi$ be a permutation of $1, 2, \\ldots, N$.\n",
    "(That is, $\\pi = (\\pi_1, \\ldots, \\pi_N)$ is a vector of length $N$ in which every number between $1$ and $N$ appears exactly once.)\n",
    "Define the vector $I_\\pi := (I_{\\pi_j})_{j=1}^N$. \n",
    "This vector has the same components as $I$, but in a different order (unless $\\pi$ is the identity permutation).\n",
    "Because the draws are IID, the joint probability distribution of this permutation of the draws is\n",
    "the same as the joint probability distribution of the original draws: the components of $I$ are _exchangeable_\n",
    "random variables if the null hypothesis is true.\n",
    "\n",
    "Equivalently, let $z$ be a vector of length $N$.\n",
    "Then $\\mathbb{P}_0 \\{I = z \\} = \\mathbb{P}_0 \\{I = z_\\pi \\}$ for all $N!$ permutations $\\pi$.\n",
    "Whatever vector of values $z$ was actually observed, if the null hypothesis is true, all permutations of those values\n",
    "were equally likely to have been observed instead.\n",
    "In turn, that implies that all $\\binom{N}{n}$ multisubsets of size $n$ of the $n+m$ components of $z$ \n",
    "were equally likely to be the sample from the first population (with the other $m$ values \n",
    "comprising the sample from the second population).\n",
    "\n",
    "Let $\\{z\\}$ denote the _multiset_ of elements of $x$. \n",
    "(See, e.g., [Mathematical Foundations](./math-foundations.ipynb)\n",
    "for the distinction between a set and a multiset.)\n",
    "Consider the event that $\\{I\\} = \\{z\\}$, that is, that the multiset of observed data is equal to the multiset of\n",
    "elements of $z$.\n",
    "Then\n",
    "\\begin{equation*}\n",
    "  \\mathbb{P}_0(I = z_\\pi | \\{I\\} = \\{z\\}) = \\mathbb{P}_0(I = z | \\{I\\} = \\{z\\})\n",
    "\\end{equation*}\n",
    "for all permutations $\\pi$ of $\\{1, \\ldots, n+m\\}$.\n",
    "\n",
    "It follows that the sample from the first population is equally likely to be any multisubset of $n$ of the \n",
    "$n+m$ elements of the pooled sample, given the elements of the pooled sample.\n",
    "That is, the conditional probability distribution of the sample from the first population is like that of a simple random sample of size $n$ from the $N$ elements of the pooled sample, given the elements of the pooled sample.\n",
    "\n",
    "The pooled sample in this example consists of $G = \\sum_{j=1}^N I_j$ 1s and $N-G$ 0s.\n",
    "Given $G$, the conditional distribution of the number of 1s in the sample from the first population is hypergeometric with parameters $N=N$ (population size), $G=G$ (number of \"good\" items in the population), and $n=n$ (sample size from the population).\n",
    "Recall that $X$ is the sum of the draws from the first population and $Y$ is the sum of the draws from the second\n",
    "population.\n",
    "We have\n",
    "\\begin{equation*}\n",
    "\\mathbb{P}_0 \\{ X = k | G=g \\} = \\left \\{ \\begin{array}{ll}\n",
    "   \\frac{\\binom{g}{k}\\binom{N-g}{n-k}}{\\binom{N}{n}}, & \\max(0, g-m) \\le k \\le \\min(n,g) \\\\\n",
    "   0, & \\mbox{otherwise.}\n",
    "   \\end{array}\n",
    "   \\right .\n",
    "\\end{equation*}\n",
    "Moreover, given $G$, $X$ and $Y$ are dependent, since $X+Y = G$.\n",
    "\n",
    "We can base a (conditional) test of the hypothesis $p_0 = p_1$ on the conditional\n",
    "hypergeometric distribution of $X$.\n",
    "Because the hypergeometric is discrete, for some values of $n$, $m$, $G$, and $\\alpha$ we will need a randomized\n",
    "test $A=A_{\\alpha;g}(\\cdot, \\cdot)$ to attain exactly level $\\alpha$.\n",
    "\n",
    "What shall we use as the acceptance region?\n",
    "To have power against the alternative that $p_x > p_y$ and the alternative $p_x < p_y$, we should\n",
    "reject if $X$ is \"too small\" or if $X$ is \"too big.\" \n",
    "There are any number of ways we could trade off small and big.\n",
    "For instance, we could make the acceptance region $A = A_{\\alpha;g}(x,u)$ nearly symmetric around \n",
    "$\\mathbb{E}(X|G) = nG/N$.\n",
    "Or we could make the chance of rejecting when $X$ is too small equal to the chance of rejecting when $X$ is too big.\n",
    "\n",
    "We will do something else: pick the acceptance region to include\n",
    "as few values in $\\{0, \\ldots, n\\}$ as possible.\n",
    "That means omitting from $A$ the possible values of $X$ that have the lowest (conditional) probabilities.\n",
    "Because the hypergeometric distribution is unimodal, this yelds a test that rejects when $X$ is \"in the tails\" of the hypergeometric distribution, as desired.\n",
    "Let $\\mathcal{K}$ denote the largest subset of $\\{0, \\ldots, n\\}$ such that \n",
    "\\begin{equation*}\n",
    "\\mathbb{P}_0 \\{X \\in \\mathcal{K} | G=g \\} < \\alpha,\n",
    "\\end{equation*}\n",
    "\n",
    "Define\n",
    "\\begin{equation*}\n",
    "\\mathcal{J} := \\arg \\min_{j \\notin \\mathcal{K}} \\mathbb{P}_0 \\{ X = j | G=g \\}\n",
    "\\end{equation*}\n",
    "be the outcome or outcomes with smallest (conditional) probability not included in $\\mathcal{K}$,\n",
    "and let \n",
    "\\begin{equation*}\n",
    "\\mathcal{I} := \\{0, \\ldots, n\\} \\setminus (\\mathcal{K} \\cup \\mathcal{J}).\n",
    "\\end{equation*}\n",
    "Let \n",
    "\\begin{equation*}\n",
    "\\gamma := \\frac{(1-\\alpha) - \\mathbb{P}_0 \\{X \\in \\mathcal{I} | G=g \\}}{\\mathbb{P}_0 \\{X \\in \\mathcal{J} |G=g\\}}\n",
    "\\end{equation*}\n",
    "be the fraction of the conditional probability of $\\mathcal{J}$ that needs to be added to the \n",
    "conditional probability of $\\mathcal{I}$ to make\n",
    "the conditional probability of their union equal $1-\\alpha$.\n",
    "Define the acceptance region \n",
    "\\begin{equation*}\n",
    "A_{\\alpha; g} = \\{ (x,u) \\in \\mathcal{I} \\times [0,1] \\} \\cup \\{ (x, u): x \\in \\mathcal{J} \\mbox{ and } u \\le \\gamma \\}.\n",
    "\\end{equation*}\n",
    "Let $U$ be a $U[0,1]$ random variable independent of $X$.\n",
    "Then\n",
    "\\begin{equation*}\n",
    "\\mathbb{P}_0 \\{ (X, U) \\in A_{\\alpha; g} | G=g \\} = 1-\\alpha.\n",
    "\\end{equation*}\n",
    "That is, $A$ is an acceptance region for a test of the null hypothesis with conditional significance level $\\alpha$. \n",
    "\n",
    "Conditioning on $G$ eliminates the nuisance parameter $p$, the fraction of ones in the two populations\n",
    "from which the sample was drawn: whatever $p$ might be, and whatever the fraction of ones in the sample, every random sample of $n$ of the $n+m$ observations is equally likely to have been the sample from the first population, if the null hypothesis is true.\n",
    "See [Introduction to Permutation Tests](./permute-intro.ipynb).\n",
    "\n",
    "This test, in a slightly different form, is called _Fisher's Exact Test_.\n",
    "It is useful for a wide range of problems, including clinical trials and A/B testing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Illustration\n",
    "\n",
    "To visualize what's going on, suppose $n=m=10$ and $G=5$. Under the null, the conditional distribution of $X$ given $G=5$ is hypergeometric with parameters $N=n+m=20$, $G=5$, and $n=10$.\n",
    "The pmf of that distribution is plotted below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAY3klEQVR4nO3df5RfdX3n8eeLiVFMoalmsJofTbSxbI4liLPh59KiguHHGrS1BKiICmmUrLhdT0m7Pfaop7t4Wl21JxKzmGJbIYuUaJRAsBRlLUIzASSEH3VOADMEN+GHEhEDgdf+ce/UL5M7M3eGufOdzLwe53zP93s/9/O53/cXQl7cX58r20RERPR3ULsLiIiI8SkBERERlRIQERFRKQERERGVEhAREVFpSrsLGE0zZszw3Llz211GRMQBY8uWLY/Z7qxaN6ECYu7cuXR3d7e7jIiIA4akhwdal0NMERFRKQERERGVEhAREVEpAREREZUaDQhJiyU9IKlH0sqK9Usk3S3pLkndkk5oWfeQpK1965qsMyIi9tfYVUySOoBVwMlAL7BZ0gbb97Z0uwnYYNuSjgCuBg5vWX+S7ceaqjEiIgbW5B7EIqDH9nbbzwLrgCWtHWz/zL+cTnYakKllIyLGiSYDYiawo2W5t2x7EUnvknQ/cB3wgZZVBm6UtEXSsoG+RNKy8vBU9+7du0ep9IiIaPJGOVW07beHYHs9sF7SicCngLeXq463vVPSYcC3Jd1v+5aK8WuANQBdXV3ZA6lh7srrXrT80KWnT+jvjYiRaXIPoheY3bI8C9g5UOfyL/83SJpRLu8s33cB6ykOWUVExBhpMiA2A/MlzZM0FVgKbGjtIOk3Jan8fBQwFXhc0jRJh5Tt04BTgHsarDUiIvpp7BCT7X2SVgCbgA5gre1tkpaX61cDvwecJ+k54BngrPKKptdQHHbqq/FK2zc0VWtEROyv0cn6bG8ENvZrW93y+dPApyvGbQcWNllbREQMLndSR0REpQRERERUSkBERESlBERERFRKQERERKUEREREVEpAREREpQRERERUSkBERESlBERERFRKQERERKUEREREVEpAREREpQRERERUSkBERESlBERERFRKQERERKUEREREVEpAREREpUYDQtJiSQ9I6pG0smL9Ekl3S7pLUrekE+qOjYiIZjUWEJI6gFXAqcAC4GxJC/p1uwlYaPtI4APA5cMYGxERDWpyD2IR0GN7u+1ngXXAktYOtn9m2+XiNMB1x0ZERLOaDIiZwI6W5d6y7UUkvUvS/cB1FHsRtcdGRERzmgwIVbR5vwZ7ve3DgTOBTw1nLICkZeX5i+7du3ePtNaIiOinyYDoBWa3LM8Cdg7U2fYtwBskzRjOWNtrbHfZ7urs7HzpVUdEBNBsQGwG5kuaJ2kqsBTY0NpB0m9KUvn5KGAq8HidsRER0awpTW3Y9j5JK4BNQAew1vY2ScvL9auB3wPOk/Qc8AxwVnnSunJsU7VGRMT+GgsIANsbgY392la3fP408Om6YyMiYuzkTuqIiKiUgIiIiEoJiIiIqJSAiIiISgmIiIiolICIiIhKCYiIiKiUgIiIiEoJiIiIqJSAiIiISgmIiIiolICIiIhKCYiIiKiUgIiIiEoJiIiIqJSAiIiISgMGhKSXj2UhERExvgy2B/F9AEl/P0a1RETEODLYI0enSnofcJykd/dfafva5sqKiIh2GywglgPnAtOB/9xvnYEERETEBDZgQNj+HvA9Sd22vzySjUtaDHwe6AAut31pv/XnApeUiz8DPmT7B+W6h4A9wPPAPttdI6khIiJGZrA9iD5XSHonMLe1v+3PDjZIUgewCjgZ6AU2S9pg+96Wbg8Cv2P7SUmnAmuAo1vWn2T7sVq/JCIiRlWdgPgm8AtgK/DCMLa9COixvR1A0jpgCfDvAWH71pb+twGzhrH9iIhoUJ2AmGX7iBFseyawo2W5lxfvHfT3QeD6lmUDN0oy8CXba6oGSVoGLAOYM2fOCMqMiIgqdW6Uu17SKSPYtiraXNlROokiIC5paT7e9lHAqcBFkk6sGmt7je0u212dnZ0jKDMiIqrUCYjbgPWSnpH0lKQ9kp6qMa4XmN2yPAvY2b+TpCOAy4Elth/va7e9s3zfBaynOGQVERFjpE5AfAY4Fnil7UNtH2L70BrjNgPzJc2TNBVYCmxo7SBpDsXlsu+1/W8t7dMkHdL3GTgFuKfWL4qIiFFR5xzED4F7bFceHhqI7X2SVgCbKC5zXWt7m6Tl5frVwMeBVwNflAS/vJz1NRR7LX01Xmn7huF8f0REvDR1AuJR4DuSrgf29jUOdZlr2WcjsLFf2+qWzxcAF1SM2w4srFFbREQ0pE5APFi+ppaviIiYBIYMCNufAJB0aLHoPY1XFRERbTfkSWpJXZK2AncDWyX9QNJbmi8tIiLaqc4hprXAh23/XwBJJwB/C4zk5rmIiDhA1LnMdU9fOMC/T+KXw0wRERNcnT2If5X0JeAqijuhz6K4qukoANt3NFhfRES0SZ2AOLJ8/4t+7cdRBMZbR7OgiIgYH+pcxXTSWBQSERHjS51zEBERMQklICIiolICIiIiKtW5Ue4iSdNbln9N0ocbrSoiItquzh7EhbZ/0rdg+0ngwsYqioiIcaFOQBykct5tAEkdZNK+iIgJr859EJuAqyWtprjvYTmQZzNERExwdQLiEuCPgA9RPGf6RopHhEZExARW50a5F4DLyldEREwSAwaEpKtt/0E51fd+jxu1ndlcIyImsMH2IC4u388Yi0IiImJ8GfAqJtuPllcsfdn2w/1fdTYuabGkByT1SFpZsf5cSXeXr1slLaw7NiIimjXoZa62nwd+LulXh7vhMlxWAacCC4CzJS3o1+1B4HfKw1WfAtYMY2xERDSozlVMv6B41Oi3gaf7Gm1/ZIhxi4Ae29sBJK0DlgD3tmzj1pb+twGz6o6NiIhm1QmI68pXq/1OWleYCexoWe4Fjh6k/weB64c7VtIyYBnAnDlzapQVERF11AmI6bY/39og6eKBOrd2q2irDBZJJ1EExAnDHWt7DeWhqa6urjrBFRERNdSZauN9FW3n1xjXC8xuWZ4F7OzfSdIRFDfeLbH9+HDGRkREcwa7D+Js4BxgnqQNLasOAR6vHvUim4H5kuYBjwBLy+21fscc4Frgvbb/bThjIyKiWYMdYroVeBSYAXympX0PcPdQG7a9T9IKirmcOoC1trdJWl6uXw18HHg18MVyPsB9trsGGjvsXxcRESM2YECU9zo8DBwr6TeA+bb/SdLBwMEUQTEo2xuBjf3aVrd8vgC4oO7YiIgYO3UeGHQhcA3wpbJpFvD1BmuKiIhxoM5J6ouA44GnAGz/EDisyaIiIqL96gTEXtvP9i1ImkK9+yAiIuIAVicgvivpz4CDJZ0MfA34ZrNlRUREu9UJiJXAbmArxYODNgJ/3mRRERHRfnUfGPS/y1dEREwSda5iOkPSnZKekPSUpD2SnhqL4iIion3qzMX0OeDdwFbbOTkdETFJ1DkHsQO4J+EQETG51NmD+BNgo6TvAnv7Gm1/trGqIiKi7eoExF8CPwNeAUxttpyIiBgv6gTEq2yf0nglERExrtQ5B/FPkhIQERGTTN25mG6Q9IvyEtdc5hoRMQnUuVHukLEoJCIixpc65yCQ9E7gxHLxO7a/1VxJERExHtS5k/pS4GLg3vJ1cdkWERETWJ09iNOAI8s5mZD0FeBOikn8IiJigqpzkhpgesvnX22gjoiIGGfq7EH8T+BOSTcDojgX8aeNVhUREW035B6E7auAY4Bry9exttfV2bikxZIekNQjab9DUpIOl/R9SXslfazfuockbZV0l6Tuej8nIiJGS52T1ALeRnEe4hvAVEmLaozrAFYBpwILgLMlLejX7QngI8BfD7CZk2wfabtrqO+LiIjRVeccxBeBY4Gzy+U9FH/xD2UR0GN7e/lM63XAktYOtnfZ3gw8V7/kiIgYC3UC4mjbFwG/ALD9JPUm7ZtJMVV4n96yrS4DN0raImnZQJ0kLZPULal79+7dw9h8REQMpk5APFceLjKApE7ghRrjVNE2nGdKHG/7KIpDVBdJOrGqk+01trtsd3V2dg5j8xERMZg6AfEFYD1wmKS/BL4H/I8a43qB2S3Ls4CddQuzvbN831V+/5DnPSIiYvTUmYvpq5K2UJyoFnCm7ftqbHszMF/SPOARYClwTp2iJE0DDrK9p/x8CvDJOmMjImJ0DBkQkl4F7AKuaml7me1BTyzb3idpBbAJ6ADW2t4maXm5frWkXwe6gUOBFyR9lOKKpxnA+uICKqYAV9q+YQS/LyIiRqjOjXJ3UBwqepJiD2I68KikXcCFtrcMNND2RmBjv7bVLZ9/THHoqb+ngIU1aouIiIbUOQdxA3Ca7Rm2X01x0vhq4MMUl8BGRMQEVCcgumxv6luwfSNwou3bgJc3VllERLRVnUNMT0i6hOJGN4CzgCfLS1/rXO4aFeauvO5Fyw9denqbKpkc8s87Yvjq7EGcQ3Ge4Ovla3bZ1gH8QVOFRUREew26B1HuJXzO9h8O0KVn9EuKiIjxYNA9CNvPA52S6kytERERE0idcxAPAf8iaQPwdF+j7c82VVRERLRfnYDYWb4OAg5ptpyIiBgv6ky18Qkopr+w/fRQ/SMiYmKo88CgYyXdC9xXLi+UlBvkIiImuDqXuX4OeAfwOIDtH1A8lzoiIiawOgGB7R39mp5voJaIiBhH6pyk3iHpOMDl5a4foTzcFBERE1edPYjlwEUUjwt9BDiyXI6IiAmszlVMjwHnjkEtERExjtS5iun1kr4pabekXZK+Ien1Y1FcRES0T51DTFdSPP/htcDrgK/R8nS5iIiYmOoEhGz/ve195esfADddWEREtFedq5hulrSS4nkQpngexHXls6qx/USD9UVERJvU2YM4C/gj4GbgO8CHgA8AW4DuwQZKWizpAUk9Zcj0X3+4pO9L2ivpY8MZGxERzapzFdO8kWy4fJbEKuBkoBfYLGmD7Xtbuj1BcV/FmSMYGxERDapzFVO3pA9Lmj7MbS8Cemxvt/0sxSGqJa0dbO+yvRl4brhjIyKiWXUOMS2luEmuW9I6Se+QpBrjZgKtU3T0lm111B4raVkZYt27d++uufmIiBjKkAFhu8f2fwfeSHHJ61rgR5I+0XeiegBVIVL36qfaY22vsd1lu6uzs7Pm5iMiYii1JuuTdATwGeCvgH8Efh94CvjnQYb1ArNblmdRPHiojpcyNiIiRsGQJ6klbQF+AnwZWGl7b7nqdknHDzJ0MzBf0jyKOZyWAufUrOuljI2IiFFQ5z6I99jeXrXC9rsHGmR7n6QVwCagA1hre5uk5eX61ZJ+neJS2UOBFyR9FFhg+6mqscP5YRER8dLUCYhHJJ0DzG3tb/uTQw20vRHY2K9tdcvnH1McPqo1NiIixk6dgPgG8FOKG+P2DtE3IiImiDoBMcv24sYriYiIcaXOVUy3SvrtxiuJiIhxZcA9CElbKe49mAK8X9J2ikNMAmz7iLEpMSIi2mGwQ0xnjFkVEREx7gwYELYfHstCIiJifKl1J3VEREw+CYiIiKiUgIiIiEoJiIiIqJSAiIiISgmIiIiolICIiIhKCYiIiKiUgIiIiEoJiIiIqJSAiIiISgmIiIiolICIiIhKjQaEpMWSHpDUI2llxXpJ+kK5/m5JR7Wse0jSVkl3Sepuss6IiNhfnUeOjoikDmAVcDLQC2yWtMH2vS3dTgXml6+jgcvK9z4n2X6sqRojImJgTe5BLAJ6bG+3/SywDljSr88S4O9cuA2YLum1DdYUERE1NRkQM4EdLcu9ZVvdPgZulLRF0rLGqoyIiEqNHWKieHZ1fx5Gn+Nt75R0GPBtSffbvmW/LynCYxnAnDlzXkq9ERHRosk9iF5gdsvyLGBn3T62+953AespDlntx/Ya2122uzo7O0ep9IiIaDIgNgPzJc2TNBVYCmzo12cDcF55NdMxwE9tPyppmqRDACRNA04B7mmw1oiI6KexQ0y290laAWwCOoC1trdJWl6uXw1sBE4DeoCfA+8vh78GWC+pr8Yrbd/QVK0REbG/Js9BYHsjRQi0tq1u+Wzgoopx24GFTdYWERGDy53UERFRKQERERGVEhAREVEpAREREZUSEBERUSkBERERlRIQERFRKQERERGVEhAREVEpAREREZUSEBERUSkBERERlRIQERFRKQERERGVGp3u+0Awd+V1L1p+6NLT21RJTET58xUHsuxBREREpQRERERUSkBERESlBERERFRKQERERKVGA0LSYkkPSOqRtLJivSR9oVx/t6Sj6o6NiIhmNRYQkjqAVcCpwALgbEkL+nU7FZhfvpYBlw1jbERENKjJPYhFQI/t7bafBdYBS/r1WQL8nQu3AdMlvbbm2IiIaJBsN7Nh6feBxbYvKJffCxxte0VLn28Bl9r+Xrl8E3AJMHeosS3bWEax98GcOXPe8vDDDzfyeyIiJiJJW2x3Va1rcg9CFW3902igPnXGFo32Gttdtrs6OzuHWWJERAykyak2eoHZLcuzgJ01+0ytMTYiIhrU5B7EZmC+pHmSpgJLgQ39+mwAziuvZjoG+KntR2uOjYiIBjW2B2F7n6QVwCagA1hre5uk5eX61cBG4DSgB/g58P7BxjZVa0RE7K+xk9Tt0NXV5e7u7naXERFxwGjXSeqIiDiAJSAiIqJSAiIiIiolICIiotKEOkktaTcw0lupZwCPjWI5B4L85olvsv1eyG8ert+wXXmX8YQKiJdCUvdAZ/InqvzmiW+y/V7Ibx5NOcQUERGVEhAREVEpAfFLa9pdQBvkN098k+33Qn7zqMk5iIiIqJQ9iIiIqJSAiIiISpM+ICQtlvSApB5JK9tdT9MkzZZ0s6T7JG2TdHG7axorkjok3Vk+yXDCkzRd0jWS7i//fR/b7pqaJum/ln+u75F0laRXtLum0SZpraRdku5paXuVpG9L+mH5/muj8V2TOiAkdQCrgFOBBcDZkha0t6rG7QP+m+3/ABwDXDQJfnOfi4H72l3EGPo8cIPtw4GFTPDfLmkm8BGgy/abKB4VsLS9VTXiCmBxv7aVwE225wM3lcsv2aQOCGAR0GN7u+1ngXXAkjbX1Cjbj9q+o/y8h+IvjZntrap5kmYBpwOXt7uWsSDpUOBE4MsAtp+1/ZO2FjU2pgAHS5oCvJIJ+CRK27cAT/RrXgJ8pfz8FeDM0fiuyR4QM4EdLcu9TIK/LPtImgu8Gbi9zaWMhc8BfwK80OY6xsrrgd3A35aH1S6XNK3dRTXJ9iPAXwM/Ah6leELlje2tasy8pnwaJ+X7YaOx0ckeEKpomxTX/Ur6FeAfgY/afqrd9TRJ0hnALttb2l3LGJoCHAVcZvvNwNOM0mGH8ao87r4EmAe8Dpgm6Q/bW9WBbbIHRC8wu2V5FhNwl7Q/SS+jCIev2r623fWMgeOBd0p6iOIw4lsl/UN7S2pcL9Bru2/v8BqKwJjI3g48aHu37eeAa4Hj2lzTWPl/kl4LUL7vGo2NTvaA2AzMlzRP0lSKE1ob2lxToySJ4rj0fbY/2+56xoLtP7U9y/Zcin/H/2x7Qv+fpe0fAzsk/VbZ9Dbg3jaWNBZ+BBwj6ZXln/O3McFPzLfYALyv/Pw+4BujsdEpo7GRA5XtfZJWAJsornhYa3tbm8tq2vHAe4Gtku4q2/7M9sb2lRQN+S/AV8v/+dkOvL/N9TTK9u2SrgHuoLha704m4LQbkq4CfheYIakX+AvgUuBqSR+kCMr3jMp3ZaqNiIioMtkPMUVExAASEBERUSkBERERlRIQERFRKQERERGVEhARIyDpfEmva3cdEU1KQESMzPkU0zlETFi5DyKiJOmPgQ+Ui5cDXwe+VU4djaSPAb8C3EMx5fIjwDPAscCbKKbXngbspbiL9zngMqCL4satP7Z9s6TzKWbb7CjHfQaYSnED417gNNtPSHoDxXT0ncDPgQtt3y/pPRQ3Rz1PMSHdiY38A4lJb1LfSR3RR9JbKO40PppiEsfbge9W9bV9TXkH/sdsd5d3Kv8f4Czbm8uptp+heP4Etn9b0uHAjZLeWG7mTRQz6b4C6AEusf1mSf8LOI9i9tk1wHLbP5R0NPBF4K3Ax4F32H5E0vTR/mcR0ScBEVE4AVhv+2kASdcC/6nm2N8CHrW9GaBvdlxJJwB/U7bdL+lhoC8gbi6fx7FH0k+Bb5btW4Ejytl2jwO+VkwrBMDLy/d/Aa6QdDXFhHQRjUhARBSqpn6fzovP0w30+EpRPU181Tb77G35/ELL8gsU/10eBPzE9pH9B9peXu5RnA7cJelI248P8l0RI5KT1BGFW4Azy5lApwHvAq4HDpP0akkvB85o6b8HOKT8fD/wOkn/EUDSIeUTzW4Bzi3b3gjMAR6oU0y5F/Jgeb4BFRaWn99g+3bbHwce48VT1keMmuxBRAC275B0BfCvZdPl5fmET1Kcj3iQIgj6XAGsltR3kvos4G8kHUxx/uHtFOcMVkvaSnGS+nzbe1sOGQ3lXOAySX8OvIziWRY/AP5K0nyKPZSbyraIUZermCIiolIOMUVERKUEREREVEpAREREpQRERERUSkBERESlBERERFRKQERERKX/D1kg5f2QZu3/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pmf: [(0, 0.016253869969040196), (1, 0.13544891640866824), (2, 0.34829721362229055), (3, 0.34829721362229055), (4, 0.13544891640866824), (5, 0.016253869969040196), (6, 0.0), (7, 0.0), (8, 0.0), (9, 0.0), (10, 0.0)]\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import hypergeom # the hypergeometric distribution\n",
    "n=10\n",
    "m=10\n",
    "N=n+m\n",
    "G=5\n",
    "x = np.arange(0,n+1)\n",
    "\n",
    "pmf = hypergeom.pmf(x,N,G,n)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.vlines(x, 0, pmf, lw=4)\n",
    "ax.set_xlabel('outcomes')\n",
    "ax.set_ylabel('hypergeometric pmf')\n",
    "plt.show()\n",
    "print('pmf: {}'.format(list(zip(x,pmf))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0645714285714292"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# probability of rejecting when X \\in {1, 5} to get overall level 0.05\n",
    "alpha = 0.05\n",
    "(alpha-(pmf[0]+pmf[5]))/(pmf[1]+pmf[4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To construct the acceptance region for a level $\\alpha = 0.05$ test, first notice that $\\mathbb{P}_0\\{X > 5 | G=5\\} = 0$, so the outcomes 6, 7, ..., 10 should be outside $\\mathcal{I}$.\n",
    "\n",
    "$\\mathbb{P}_0\\{X = 0 | G=5\\} = \\mathbb{P}_0\\{X = 5 | G=5\\} = 0.01625$, so if we exclude those outcomes from $\\mathcal{I}$, the chance of rejecting the null is $2 \\times 0.01625387 = 0.032508$.\n",
    "For the test to have level $\\alpha$, we need to reject more often: an additional $0.017492$ of the time.\n",
    "If we always rejected when $X=1$ or $X=4$, we would reject too often, because each of those possibilities has probability $0.13545$. \n",
    "If when $X \\in \\{1, 4\\}$ we reject with probability $\\gamma = .01749/(2\\times 0.13545) = 0.06457$, the overall chance of erroneously rejecting the null will be 5%, as desired."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unconditional tests from conditional tests\n",
    "\n",
    "If you always test conditionally at level not greater than $\\alpha$, that yields a test that has unconditional level not greater than $\\alpha$, as we shall see.\n",
    "\n",
    "Suppose we condition on the value of a discrete random variable, such as $G$ in the previous example, then test at level not greater than $\\alpha$, conditional on the value of $G$.\n",
    "The unconditional significance level of the test is\n",
    "\\begin{equation*}\n",
    "\\mathbb{P}_0 \\{\\mbox{reject null}\\} =\n",
    "\\sum_g \\mathbb{P}_0 \\{ \\mbox{reject null} | G=g \\} \\mathbb{P}_0 \\{G=g\\}\n",
    "\\le \\sum_g \\alpha \\mathbb{P}_0 \\{G=g\\} = \\alpha \\sum_g \\mathbb{P}_0 \\{G=g\\} = \\alpha,\n",
    "\\end{equation*}\n",
    "where the sums are over all values $g$ that $G$ can take.\n",
    "\n",
    "A similar proof establishes the result when $G$ does not have a discrete distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numerical comparison\n",
    "\n",
    "We saw that the $Z$-test does not necessarily have its nominal level in this problem. \n",
    "We now implement the permutation test (a randomized version of Fisher's exact test) for comparison."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Algorithmic considerations\n",
    "\n",
    "#### Finding $\\mathcal{I}$ and $\\mathcal{J}$\n",
    "\n",
    "In principle, we could sort the $n+1$ values of the pmf to find the outcomes with the largest probabilities to include in $\\mathcal{I}$, and the largest-probability outcomes not included in $\\mathcal{I}$ to be $\\mathcal{J}$ (if the test needs to be randomized to attain the desired significance level).\n",
    "\n",
    "But the best sorting algorithms still take $n\\log n$ operations, and we don't need a complete sort--we only need to divide the probabilities into two groups, the biggest and the rest.\n",
    "\n",
    "However, because the hypergeomtric distribution is unimodal, the smallest probability will either be $\\mathbb{P}_0 \\{X = 0\\}$ or $\\mathbb{P}_0 \\{X = n\\}$ (or they will be equal). \n",
    "After removing that outcome from consideration, the second-smallest probability will be either that of the largest remaining outcome or the smallest remaining outcome, etc. This means we can identify the outcomes to include in $\\mathcal{I}$ with a number of operations that is linear in $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simulate the significance level of the permutation test\n",
    "def fisher_accept(N: int, G: int, n: int, alpha: float=0.05) -> tuple:\n",
    "    '''\n",
    "    Acceptance region for randomized hypergeometric test\n",
    "    \n",
    "    Find the acceptance region for a randomized, exact level alpha test of \n",
    "    the null hypothesis X~Hypergeometric(N, G, n). The acceptance region is\n",
    "    the smallest possible. (And not, for instance, symmetric.)\n",
    "\n",
    "    If a non-randomized, conservative test is desired, use the union of I and J \n",
    "    as the acceptance region.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    N:  integer\n",
    "        population size\n",
    "    G:  integer\n",
    "        number of \"good\" items in the population\n",
    "    n:  integer\n",
    "        sample size\n",
    "    alpha : float\n",
    "        desired significance level    \n",
    "  \n",
    "    Returns\n",
    "    --------\n",
    "    I:  list\n",
    "        observed values for which the test never rejects\n",
    "    J:  list \n",
    "        observed values for which the test sometimes rejects\n",
    "    gamma : float\n",
    "        probability the test does not reject when the observed value is in J\n",
    "    '''\n",
    "    x = np.arange(0, n+1)\n",
    "    I = list(x)                    # start with all possible outcomes, then remove some\n",
    "    pmf = hypergeom.pmf(x,N,G,n)   # hypergeometric pmf\n",
    "    bottom = 0                     # smallest outcome still in I\n",
    "    top = n                        # largest outcome still in I\n",
    "    J = []                         # outcomes for which the test is randomized\n",
    "    p_J = 0                        # probability of the outcomes for which test is randomized\n",
    "    p_tail = 0                     # probability of outcomes not in I\n",
    "    while p_tail < alpha:          # still need to remove outcomes from the acceptance region\n",
    "        pb = pmf[bottom]\n",
    "        pt = pmf[top]\n",
    "        if pb < pt:                # the smaller possibility has smaller probability\n",
    "            J = [bottom]\n",
    "            p_J = pb\n",
    "            bottom += 1\n",
    "        elif pb > pt:              # the larger possibility has smaller probability\n",
    "            J = [top]\n",
    "            p_J = pt\n",
    "            top -= 1\n",
    "        else:                      \n",
    "            if bottom < top:       # the two possibilities have equal probability\n",
    "                J = [bottom, top]\n",
    "                p_J = pb+pt\n",
    "                bottom += 1\n",
    "                top -= 1\n",
    "            else:                  # there is only one possibility left\n",
    "                J = [bottom]\n",
    "                p_J = pb\n",
    "                bottom +=1\n",
    "        p_tail += p_J\n",
    "        for j in J:\n",
    "            I.remove(j)\n",
    "    gamma = (p_tail-alpha)/p_J     # probability of accepting H_0 when X in J to get \n",
    "                                   # exact level alpha\n",
    "    return I, J, gamma\n",
    "        \n",
    "    \n",
    "def apply_test(A, x, U):\n",
    "    '''\n",
    "    Apply a randomized test A to data x using auxiliary uniform randomness U\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    A:  triple\n",
    "        first element is a list or set, the values of x for which the test never rejects\n",
    "        second element is a list or set, the values of x for which the test sometimes rejects\n",
    "        third element is a float in [0, 1], the probability of rejecting when x is in the second element\n",
    "    x:  number-like\n",
    "        observed data\n",
    "    U:  float in [0, 1]\n",
    "        observed value of an independent U[0,1] random variable\n",
    "    Returns\n",
    "    -------\n",
    "    reject: Boolean\n",
    "        True if the test rejects\n",
    "    '''\n",
    "    assert 0 <= A[2] < 1, 'probability out of range:{}'.format(A[2])\n",
    "    assert 0 <= U <= 1, 'uniform variable out of range:{}'.format(U)\n",
    "    return not(x in A[0] or (x in A[1] and U <= A[2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise: write unit tests of the two functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max SD of p-hat: 0.0005\n"
     ]
    }
   ],
   "source": [
    "# set up the simulations\n",
    "alpha = 0.05 # significance level 5%\n",
    "p = [.01, .05, .1, .3, .4, .5, .6, .7, .9, .95, .99]  # assortment of values for p\n",
    "\n",
    "reps = int(10**6)\n",
    "n = 50\n",
    "m = 100\n",
    "N = n+m\n",
    "\n",
    "# how precise should we expect the result to be?\n",
    "print('max SD of p-hat: {}'.format((1/(2*math.sqrt(reps)))))  # max SD of a binomial is 1/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p=0.01 estimated alpha=0.05001\n",
      "p=0.05 estimated alpha=0.049823\n",
      "p=0.1 estimated alpha=0.050227\n",
      "p=0.3 estimated alpha=0.049927\n",
      "p=0.4 estimated alpha=0.050298\n",
      "p=0.5 estimated alpha=0.049955\n",
      "p=0.6 estimated alpha=0.049959\n",
      "p=0.7 estimated alpha=0.049856\n",
      "p=0.9 estimated alpha=0.049854\n",
      "p=0.95 estimated alpha=0.050085\n",
      "p=0.99 estimated alpha=0.049817\n",
      "CPU times: user 47min 31s, sys: 15.4 s, total: 47min 46s\n",
      "Wall time: 47min 56s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# brute force: re-compute the test for each replication\n",
    "for pj in p:\n",
    "    reject = 0\n",
    "    for i in range(reps):\n",
    "        X = rng.binomial(n, pj) # simulate X\n",
    "        Y = rng.binomial(m, pj) # simulate Y\n",
    "        A = fisher_accept(N, X+Y, n, alpha=alpha)\n",
    "        reject += (1 if apply_test(A, X, rng.uniform())\n",
    "                   else 0)\n",
    "    print('p={} estimated alpha={}'.format(pj, reject/reps))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This test has true significance level equal to its nominal significance level: it is an _exact_ test rather than an _approximate_ test.\n",
    "\n",
    "Randomized tests have some drawbacks in practice.\n",
    "For instance, it would be hard to explain to a consulting client or a judge that for a given set of data, the test you\n",
    "propose to use sometimes\n",
    "rejects the null hypothesis and sometimes does not, depending on a random factor ($U$) that has nothing to do \n",
    "with the data or the experiment.\n",
    "\n",
    "In such situations, it might make sense to use a _conservative_ test rather than an exact randomized test.\n",
    "A conservative test is one for which the chance of a Type I error is not greater than $\\alpha$.\n",
    "In the previous example, if the test rejects only when the test "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Algorithmic considerations\n",
    "\n",
    "The previous simulation is **slow**, in part because it finds the conditional acceptance region \n",
    "`reps` times for each `pj in p`.\n",
    "\n",
    "But there are relatively few possible acceptance regions, one for each possible value of $G$, i.e., $n+m+1$.\n",
    "\n",
    "What happens if we pre-compute the acceptance regions for all possible values of $G$, and look them up as needed?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p=0.01 estimated alpha=0.050358\n",
      "p=0.05 estimated alpha=0.050045\n",
      "p=0.1 estimated alpha=0.049983\n",
      "p=0.3 estimated alpha=0.050171\n",
      "p=0.4 estimated alpha=0.05002\n",
      "p=0.5 estimated alpha=0.049699\n",
      "p=0.6 estimated alpha=0.049793\n",
      "p=0.7 estimated alpha=0.049458\n",
      "p=0.9 estimated alpha=0.049673\n",
      "p=0.95 estimated alpha=0.050049\n",
      "p=0.99 estimated alpha=0.050015\n",
      "CPU times: user 1min 17s, sys: 295 ms, total: 1min 18s\n",
      "Wall time: 1min 18s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# Smarter approach: pre-compute the tests. There are only n+m+1 possible tests, far fewer than reps\n",
    "# Include the \"cost\" of precomputing the tests in the comparison.\n",
    "\n",
    "AA = []\n",
    "for j in range(N+1):\n",
    "    AA.append(fisher_accept(N, j, n, alpha=alpha))\n",
    "\n",
    "for pj in p:\n",
    "    reject = 0\n",
    "    for i in range(reps):\n",
    "        X = rng.binomial(n, pj) # simulate X\n",
    "        Y = rng.binomial(m, pj) # simulate Y\n",
    "        reject += (1 \n",
    "                   if apply_test(AA[X+Y], X, rng.uniform()) \n",
    "                   else 0)\n",
    "    print('p={} estimated alpha={}'.format(pj, reject/reps))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This easy change speeds up the simulation by a factor of more than 40.\n",
    "\n",
    "Alternatively, we can use the `functools` library to automatically cache the acceptance regions as they are computed.\n",
    "\n",
    "Here's the acceptance region function with the `@lru_cache` decorator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import lru_cache\n",
    "\n",
    "@lru_cache(maxsize=None)  # decorate the function to cache the results of calls to the function\n",
    "def fisher_accept(N, G, n, alpha=0.05):\n",
    "    '''\n",
    "    Acceptance region for randomized hypergeometric test\n",
    "    \n",
    "    Find the acceptance region for a randomized, exact level alpha test of \n",
    "    the null hypothesis X~Hypergeometric(N, G, n). The acceptance region is\n",
    "    the smallest possible. (And not, for instance, symmetric.)\n",
    "\n",
    "    If a non-randomized, conservative test is desired, use the union of I and J as \n",
    "    the acceptance region.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    N:  integer\n",
    "        population size\n",
    "    G:  integer\n",
    "        number of \"good\" items in the population\n",
    "    n:  integer\n",
    "        sample size\n",
    "    alpha : float\n",
    "        desired significance level    \n",
    "  \n",
    "    Returns\n",
    "    --------\n",
    "    I:  list\n",
    "        values for which the test never rejects\n",
    "    J:  list \n",
    "        values for which the test sometimes rejects\n",
    "    gamma : float\n",
    "        probability the test does not reject when the value is in J\n",
    "    '''\n",
    "    x = np.arange(0, n+1)          # all possible values of X\n",
    "    I = list(x)                    # start with all possible outcomes, then remove some\n",
    "    pmf = hypergeom.pmf(x,N,G,n)   # hypergeometric pmf\n",
    "    bottom = 0                     # smallest outcome still in I\n",
    "    top = n                        # largest outcome still in I\n",
    "    J = []                         # outcome for which the test is randomized\n",
    "    p_J = 0                        # probability of the randomized outcome\n",
    "    p_tail = 0                     # probability of outcomes excluded from I\n",
    "    while p_tail < alpha:          # still need to remove outcomes from the acceptance region\n",
    "        pb = pmf[bottom]\n",
    "        pt = pmf[top]\n",
    "        if pb < pt:                # the lower possibility has smaller probability\n",
    "            J = [bottom]\n",
    "            p_J = pb\n",
    "            bottom += 1\n",
    "        elif pb > pt:              # the upper possibility has smaller probability\n",
    "            J = [top]\n",
    "            p_J = pt\n",
    "            top -= 1\n",
    "        else:                      \n",
    "            if bottom < top:       # the two possibilities have equal probability\n",
    "                J = [bottom, top]\n",
    "                p_J = pb+pt\n",
    "                bottom += 1\n",
    "                top -= 1\n",
    "            else:                  # there is only one possibility left\n",
    "                J = [bottom]\n",
    "                p_J = pb\n",
    "                bottom +=1\n",
    "        p_tail += p_J\n",
    "        for j in J:\n",
    "            I.remove(j)\n",
    "    gamma = (p_tail-alpha)/p_J     # probability of accepting H_0 when X in J to get exact level alpha\n",
    "    return I, J, gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p=0.01 estimated alpha=0.050092\n",
      "p=0.05 estimated alpha=0.050045\n",
      "p=0.1 estimated alpha=0.049978\n",
      "p=0.3 estimated alpha=0.050282\n",
      "p=0.4 estimated alpha=0.050424\n",
      "p=0.5 estimated alpha=0.049985\n",
      "p=0.6 estimated alpha=0.050117\n",
      "p=0.7 estimated alpha=0.049918\n",
      "p=0.9 estimated alpha=0.049929\n",
      "p=0.95 estimated alpha=0.050044\n",
      "p=0.99 estimated alpha=0.0501\n",
      "CPU times: user 1min 21s, sys: 329 ms, total: 1min 21s\n",
      "Wall time: 1min 22s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# re-run the original \"brute-force\" code but with function caching\n",
    "for pj in p:\n",
    "    reject = 0\n",
    "    for i in range(reps):\n",
    "        X = rng.binomial(n, pj) # simulate X\n",
    "        Y = rng.binomial(m, pj) # simulate Y\n",
    "        A = fisher_accept(N, X+Y, n, alpha=alpha)\n",
    "        reject += (1 if apply_test(A, X, rng.uniform())\n",
    "                   else 0)\n",
    "    print('p={} estimated alpha={}'.format(pj, reject/reps))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, this speeds up the simulation by a factor of more than 40, with even less coding effort. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise: How would you find a $P$-value for this test?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}