confidence-sets.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Confidence Sets\n",
    "\n",
    "## This is a rough work in progress!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Types of Parameters\n",
    "\n",
    "Many different kinds of things are called \"parameters.\" Here are several categories.\n",
    "\n",
    "### Population parameters\n",
    "\n",
    "Any property of a population may be called a _parameter_.\n",
    "Examples include the population mean, percentiles, number of modes, etc.\n",
    "\n",
    "If the population has more than one \"value\" per item, a parameter could involve more than one of them. E.g., if the population is a group of people each of whom has a height and a weight, then the population correlation between height and weight is a parameter.\n",
    "\n",
    "Similarly, consider a group of individuals and the values of some quantity (the \"response\") for each of those individuals without and with some intervention (notionally, a \"treatment\").\n",
    "The difference between the average response without the intervention and the average response with the intervention is a parameter (the _average treatment effect_).\n",
    "\n",
    "If we are sampling at random from a population, the probability distribution of the sample\n",
    "depends on the values in the population, and thus, in general, on the \n",
    "values of population parameters.\n",
    "(It will also depend on the sampling design.)\n",
    "\n",
    "### Functional parameters of probability distributions\n",
    "\n",
    "Suppose $X \\sim \\mathbb{P}$, where $\\mathbb{P}$ is a probability distribution on some space $\\mathcal{X}$ of possible outcomes.\n",
    "We assume that $\\mathbb{P} \\in \\mathcal{P}$, some known set of possible distributions.\n",
    "\n",
    "A _functional parameter_ $\\theta(\\mathbb{P})$ is a function of $\\mathbb{P}$.\n",
    "For instance the (population) mean is a functional parameter:\n",
    "\n",
    "\\begin{equation}\n",
    "\\theta(\\mathbb{P}) = \\mathbb{E}X \\equiv \\int_\\mathcal{X} x d\\mathbb{P}(x).\n",
    "\\end{equation}\n",
    "\n",
    "So are other moments of the probability distribution:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\theta(\\mathbb{P}) = \\mathbb{E}X^n \\equiv \\int_\\mathcal{X} x^n d\\mathbb{P}(x), \\;\\; n=1, 2, \\ldots .\n",
    "\\end{equation*}\n",
    "\n",
    "Other properties of $\\mathbb{P}$, such as percentiles of a univariate\n",
    "distribution, are also functional parameters.\n",
    "For instance, if $X$ is a real-valued random variable,\n",
    "then the $\\alpha$ percentile of $\\mathbb{P}$,\n",
    "\\begin{equation}\n",
    "\\theta(\\mathbb{P}) = \\inf \\left \\{x: \\int_{-\\infty}^x d\\mathbb{P}(x) \\ge \\alpha \\right \\},\n",
    "\\end{equation}\n",
    "is a functional parameter.\n",
    "\n",
    "For multivariate distributions, correlations among the components of $X$\n",
    "are functional parameters.\n",
    "\n",
    "In general, there can be distinct distributions $\\mathbb{P}$ and $\\mathbb{Q}$ such that \n",
    "$\\mathbb{P} \\ne \\mathbb{Q}$ but $\\theta(\\mathbb{P}) = \\theta(\\mathbb{Q})$. \n",
    "For instance there are infinitely many normal distributions with the\n",
    "same mean (but different variances).\n",
    "\n",
    "### Parameters as indices for sets of distributions\n",
    "\n",
    "Another use of the term \"parameter\" is as an abstract\n",
    "index that points to a particular distribution in\n",
    "a family of distributions.\n",
    "For instance, we might have a multiset of distributions\n",
    "$\\mathcal{P} = \\{\\mathbb{P}_\\eta\\}_{\\eta \\in \\Theta}$.\n",
    "In that case, $\\eta$ is an index parameter.\n",
    "For index parameters, if for all parameters $\\eta$, $\\nu \\in \\Theta$ such that\n",
    "$\\eta \\ne \\nu$, $\\mathbb{P}_\\eta \\ne \\mathbb{P}_\\nu$, the parameter is said to be _identifiable_.\n",
    "That is, $X$ contains enough information to\n",
    "identify the value of the parameter with arbitrarily high accuracy, given enough\n",
    "observations.\n",
    "Otherwise, the parameter is _non-identifiable_ or _unidentifiable_: the data\n",
    "do not contain enough information to distinguish among different values\n",
    "of the parameter, no matter how many observations are made.\n",
    "\n",
    "\n",
    "#### Special case: location-scale families\n",
    "\n",
    "Many indexed families of distributions are related through the\n",
    "value of their parameter in a particular way. \n",
    "For instance, suppose that the outcome space $\\mathcal{X}$ is a real vector space,\n",
    "so it makes sense to add elements of $\\mathcal{X}$ and to multiply them by scalars.\n",
    "\n",
    "If $X \\sim \\mathbb{P}$, then for any $x \\in \\mathcal{X}$ and $a \\in \\Re \\backslash \\{0\\}$,\n",
    "we could define $\\mathbb{P}_{x,a}$ to be the distribution of $aX+x$.\n",
    "Then $\\{P_{x, a} \\}_{x \\in \\mathcal{X}, a \\in \\Re \\backslash \\{0\\}}$ is a \n",
    "_location-scale family_\n",
    "with parameter $\\theta = (x, a)$\n",
    "As $x$ varies, the probability distribution \"shifts\" its location.\n",
    "As $a$ varies, the probability distribution $P$ is \"stretched\" or re-scaled.\n",
    "The family of univariate normal distributions is a location-scale family over the two-dimensional parameter $\\theta = (\\mu, \\sigma)$ with $\\mu \\in \\Re$, $\\sigma \\in \\Re \\backslash \\{0\\}$.\n",
    "\n",
    "#### Notation for index parameters\n",
    "To keep the notation for index parameters\n",
    "parallel with the notation for functional parameters, we will define $\\theta(\\mathbb{P}) \\equiv \n",
    "\\{ \\eta: \\mathbb{P} = \\mathbb{P}_\\eta\\}$.\n",
    "If $\\theta$ is identifiable, $\\theta(\\mathbb{P})$ is a singleton set; otherwise,\n",
    "it may contain more than one element. \n",
    "\n",
    "### Parametric families of distributions\n",
    "\n",
    "A _parametric family of distributions_ is an indexed collection of probability distributions\n",
    "that depends on the index parameter (which might be multidimensional) in a \n",
    "fixed functional way.\n",
    "(We can think of things like the mean and standard deviation of a normal distribution as either a multidimensional parameter or as a collection of parameters.)\n",
    "\n",
    "Most distributions that have names are parametric families, e.g., Bernoulli (the parameter $p$), Binomial (the two-dimensional parameter $(n, p)$), Geometric ($p$), Hypergeometric (the three-dimensional parameter $(N, G, n)$), Negative Binomial $(p, k)$, Normal $(\\mu, \\sigma)$, Student's $T$ $(\\mu, \\sigma, \\nu)$, continuous uniform (the endpoints of the interval of support, the two-dimensional parameter $(a, b)$), and so on. \n",
    "\n",
    "### Nuisance parameters\n",
    "\n",
    "When the probability distribution of the data depends on a multi-dimensional parameter\n",
    "but only some components of that parameter are of interest, the other components are called _nuisance parameters_. \n",
    "For instance, in estimating the mean of a normal distribution, the\n",
    "variance of the distribution is a nuisance parameter: we don't care what it is,\n",
    "but it affects the probability distribution of the data.\n",
    "\n",
    "Similarly, in estimating a population mean from a stratified sample, the means\n",
    "within the different strata are nuisance parameters.\n",
    "\n",
    "\n",
    "### Abstract parameters\n",
    "\n",
    "For most of the theory in this chapter, $\\theta$ will be an abstract parameter:\n",
    "the development applies to functional parameters, index parameters, parameters of\n",
    "parametric families, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confidence sets\n",
    "\n",
    "What can we learn about the value of $\\theta(\\mathbb{P})$ from observations?\n",
    "[The chapter on testing](./tests.ipynb) discusses testing hypotheses, including\n",
    "hypotheses about parameters.\n",
    "Here we explore a different approach to quantifying what a sample tells us about\n",
    "$\\theta$: confidence sets.\n",
    "The treatment will be abstract but informal. \n",
    "(For instance, we shall ignore measurability issues.)\n",
    "\n",
    "In an abuse of notation, we will let $\\theta$ denote both the value of a parameter, and\n",
    "the mapping from a distribution to the value of the parameter for that distribution,\n",
    "as if $\\theta$ were a functional parameter even if it is an index parameter (or some other\n",
    "kind of parameter).\n",
    "Thus, $\\theta: \\mathcal{P} \\rightarrow \\Theta$, $\\mathbb{P} \\mapsto \\theta(\\mathbb{P})$.\n",
    "If $\\mathbb{P} = \\mathbb{P}_\\eta$, then $\\theta(\\mathbb{P}) = \\eta$.\n",
    "The set $\\Theta$ will denote the possible values of $\\theta$. \n",
    "Lowercase Greek letters such as $\\eta$ will denote\n",
    "generic elements of $\\Theta$.\n",
    "\n",
    "We shall observe $X \\sim \\mathbb{P}$, where $X$ takes values in the outcome space $\\mathcal{X}$.\n",
    "We do not know $\\mathbb{P}$, but we know that $\\mathbb{P} \\in \\mathcal{P}$, a known set of distributions.\n",
    "Let $\\mathcal{I}(\\cdot)$ be a set-valued function that assigns a subset of $\\Theta$ to each possible observation $x \\in \\mathcal{X}$.\n",
    "For instance, we might observe $X \\sim N(\\theta, 1)$, and $\\mathcal{I}(x)$ might be $[x-c, x+c]$.\n",
    "\n",
    "Fix $\\alpha \\in (0, 1)$.\n",
    "Suppose that for all $\\eta \\in \\Theta$, if $\\theta(\\mathbb{P}) = \\eta$ then\n",
    "\\begin{equation}\n",
    "\\mathbb{P} \\{\\mathcal{I}(X) \\ni \\eta \\} \\ge 1-\\alpha.\n",
    "\\end{equation}\n",
    "Then $\\mathcal{I}(\\cdot)$ is _a $1-\\alpha$ confidence set procedure_ for $\\theta(\\mathbb{P})$.\n",
    "It maps outcomes to sets in such a way that the chance is at least\n",
    "$1-\\alpha$ that the resulting set will contain the true value of the\n",
    "parameter $\\theta(\\mathbb{P})$.\n",
    "\n",
    "If we observe $X=x$, $\\mathcal{I}(x)$ is _a $1-\\alpha$ confidence set for $\\theta$_.\n",
    "The _confidence level_ of the set is $1-\\alpha$.\n",
    "\n",
    "When $\\mathcal{I}(x) \\ni \\theta$, we say that the confidence set _covers_ $\\theta$.\n",
    "The _coverage probability_ of the confidence set procedure $\\mathcal{I}$\n",
    "is \n",
    "\\begin{equation*}\n",
    "\\inf_{\\eta \\in \\Theta} \\inf_{\\mathbb{P} \\in \\mathcal{P} : \\theta(\\mathbb{P}) = \\eta} \\mathbb{P} \\{\\mathcal{I}(X) \\ni \\eta \\}.\n",
    "\\end{equation*}\n",
    "\n",
    "Before the data $X$ are observed, the chance that $\\mathcal{I}(X)$ will contain $\\theta$\n",
    "is the coverage probability, $1-\\alpha$. \n",
    "After the data $X=x$ are observed, the set $\\mathcal{I}(x)$ either\n",
    "does or does not contain $\\theta$: there is nothing random anymore.\n",
    "\n",
    "A _confidence interval_ is a special case of a confidence set, when the set is an interval of real numbers.\n",
    "_One-sided confidence intervals_ are the special case that the confidence set is a semi-infinite interval, i.e., a set of real numbers of the form $(-\\infty, c]$ or $[c, \\infty)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: confidence interval for a Normal mean\n",
    "\n",
    "Suppose $X \\sim N(\\theta, 1)$: $\\theta$ is an index parameter for the (parametric) family of unit variance normal distributions ($\\mathcal{P} \\equiv \\{N(\\eta, 1)\\}_{\\eta \\in \\Re}$) and also a functional parameter, since $\\mathbb{E} X = \\theta$.\n",
    "\n",
    "Define $\\mathcal{I}(x) \\equiv [x - z_{1-\\alpha/2}, x + z_{1-\\alpha/2}]$, where\n",
    "$z_{1-\\alpha/2}$ is the $1-\\alpha/2$ percentile of the standard normal distribution.\n",
    "Then \n",
    "\\begin{equation}\n",
    "   \\mathbb{P}_\\theta \\{ \\mathcal{I}(X) \\ni \\theta \\} = 1-\\alpha\n",
    "\\end{equation}\n",
    "whatever be $\\theta \\in \\Re$.\n",
    "Thus $[x - z_{1-\\alpha/2}, x + z_{1-\\alpha/2}]$ is a $1-\\alpha$\n",
    "confidence interval for $\\theta$.\n",
    "\n",
    "Why is the coverage probability of $[X - z_{1-\\alpha/2}, X + z_{1-\\alpha/2}]$ \n",
    "equal to $1-\\alpha$?\n",
    "\n",
    "The distribution of $X-\\theta$ is a standard normal ($N(0,1)$), so\n",
    "\\begin{equation*}\n",
    "\\mathbb{P}_\\theta  \\{|X-\\theta| \\le z_{1-\\alpha/2} \\} = 1-\\alpha.\n",
    "\\end{equation*}\n",
    "But whenever $|X-\\theta| \\le z_{1-\\alpha/2}$,\n",
    "the interval $[X - z_{1-\\alpha/2}, X + z_{1-\\alpha/2}]$ contains $\\theta$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Duality between hypothesis tests and confidence sets\n",
    "\n",
    "One of the most versatile ways of constructing confidence sets is to _invert_\n",
    "hypothesis tests.\n",
    "\n",
    "Suppose we have a (possibly randomized) family of significance level $\\alpha$ hypothesis tests for all possible values of a parameter $\\theta \\in \\Theta$. \n",
    "That is, for each $\\eta \\in \\Theta$, we have a _test function_ (aka _critical function_)\n",
    "$\\phi_\\eta : \\mathcal{X} \\rightarrow [0, 1]$ such that if $\\theta(\\mathbb{P}) = \\eta$,\n",
    "\\begin{equation*}\n",
    "\\mathbb{E}_{\\mathbb{P}} \\phi_\\eta(X) \\ge 1-\\alpha.\n",
    "\\end{equation*}\n",
    "The test function\n",
    "$\\phi_\\eta(x)$ is the probability of not rejecting the hypothesis $\\theta(\\mathbb{P}) = \\eta$ \n",
    "when $X=x$.\n",
    "When $\\phi_\\eta(X) = 0$, we certainly reject the null; when \n",
    "$\\phi_\\eta(X)=1$, we certainly do not reject the null; values between 0 and 1 correspond to\n",
    "rejecting the null hypothesis with probability $1-\\phi(X)$.\n",
    "The test involves both $X$ and a uniformly distributed random variable $U \\sim U[0,1]$\n",
    "independent of $X$. The test rejects if $U \\ge \\phi(X)$.\n",
    "See [the chapter on hypothesis tests](./tests.ipynb).\n",
    "\n",
    "Consider the set \n",
    "\\begin{equation*}\n",
    "\\mathcal{I}(X,U) \\equiv \\{ \\eta \\in \\Theta : \\phi_\\eta(X) > U \\}.\n",
    "\\end{equation*}\n",
    "That is, $\\mathcal{I}$ is the set of possible parameters $\\eta \\in \\Theta$ for which the corresponding test $\\phi_\\eta$ does not reject the hypothesis that $\\theta(\\mathbb{P}) = \\eta$, for the observed values of $X$ and $U$.\n",
    "(If $\\phi$ can only take the values 0 and 1, i.e., if the test is not randomized, then $\\mathcal{I}$ does not depend on $U$.)\n",
    "\n",
    "**Claim:** $\\mathcal{I}(X,U)$ is a $1-\\alpha$ confidence procedure. That is,\n",
    "whatever the true value of $\\theta(\\mathbb{P})$ happens to be,\n",
    "\\begin{equation}\n",
    "\\mathbb{P} \\{ \\mathcal{I}(X,U) \\ni \\theta(\\mathbb{P}) \\} \\ge 1-\\alpha.\n",
    "\\end{equation}\n",
    "\n",
    "**Proof:** The set \n",
    "$\\mathcal{I}(X,U) \\ni \\theta(\\mathbb{P})$ whenever the test of the null hypothesis that $\\theta(\\mathbb{P}) = \\theta$\n",
    "does not reject, that is, when $\\phi_\\theta(X) \\ge U$.\n",
    "But when $\\theta(\\mathbb{P}) = \\theta$, \n",
    "\\begin{equation*}\n",
    "\\mathbb{P} \\{\\phi_\\theta(X) \\ge U\\} = \\mathbb{E}_{\\mathbb{P}} \\phi(X) \\ge 1-\\alpha.\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, there are many ways to construct a set of tests $\\{\\phi_\\eta\\}_{\\eta \\in \\Theta}$ with significance level $\\alpha$.\n",
    "Inverting different tests will lead to confidence sets with different properties.\n",
    "As a simple example, inverting 1-sided tests for a real parameter will give a 1-seded confidence interval for the parameter, while inverting 2-sided tests will yield a 2-sided confidence interval.\n",
    "\n",
    "It is often possible to design confidence sets that have desirable characteristics--such as avoiding including zero to the extent possible (so that they determine the sign of their parameter), or being on average as small as possible--by inverting suitably chosen tests. \n",
    "See, e.g., [Benjamini, Hochberg, and Stark (1996)](https://amstat.tandfonline.com/doi/abs/10.1080/01621459.1998.10474112#.YGk8ERRKg-Q); [Benjamini and Stark (1996)](https://www.tandfonline.com/doi/abs/10.1080/01621459.1996.10476692); \n",
    "[Evans, Hansen, and Stark (2005)](https://projecteuclid.org/journals/bernoulli/volume-11/issue-4/Minimax-expected-measure-confidence-sets-for-restricted-location-parameters/10.3150/bj/1126126761.full); and [Benjamini, Hechtlinger, and Stark (2019)](https://arxiv.org/abs/1906.00505)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: confidence interval for Binomial $p$ (known $n$)\n",
    "\n",
    "Suppose we will observe $X \\sim \\mbox{Binom}(n, p)$, where $n$ is known but $p$ is not.\n",
    "We seek a one-sided lower confidence interval for $p$, that is, a set of the form\n",
    "$[f(X,U), \\infty)$ such that for all $q \\in [0, 1]$, if $X \\sim \\mbox{Binom}(n, q)$\n",
    "and $U$ is an independent uniform random variable, then\n",
    "\\begin{equation*}\n",
    "\\mathbb{P} \\{ [f(X,U), \\infty) \\ni q \\} \\ge 1-\\alpha.\n",
    "\\end{equation*}\n",
    "Since it is certain that $p \\le 1$, the upper endpoint of the interval can be reduced from $\\infty$ to 1 without sacrificing coverage probability. That is, the same $f$ will satisfy\n",
    "\\begin{equation*}\n",
    "\\mathbb{P} \\{ [f(X,U), 1] \\ni q \\} \\ge 1-\\alpha.\n",
    "\\end{equation*}\n",
    "\n",
    "To make such a lower confidence bound, we can invert one-sided hypothesis tests that\n",
    "reject when $X$ is \"too big.\"\n",
    "That is, we want a family of tests of the hypotheses $p = q$ for all $q \\in [0, 1]$\n",
    "that reject for large values of $X$. Such tests give evidence that $p$ is _at least_ a given size.\n",
    "\n",
    "To keep things simple, we will use conservative non-randomized tests rather than exact randomized tests. \n",
    "Because we are basing the confidence intervals on\n",
    "_conservative_ tests, we expect the coverage probability to be greater than $1-\\alpha$.\n",
    "We reject the hypothesis $p = q$ if, on the assumption that $p=q$, the chance that\n",
    "$X$ would be greater than or equal to its observed value is not greater than $\\alpha$.\n",
    "That is, \n",
    "\\begin{equation*}\n",
    "\\phi_q(x) = \\left \\{ \\begin{array}{ll}\n",
    "                      1, & \\sum_{k=x}^n \\binom{n}{k}q^k(1-q)^{n-k} \\ge \\alpha \\\\\n",
    "                      0, & \\mbox{otherwise.}\n",
    "                      \\end{array}\n",
    "            \\right .\n",
    "\\end{equation*}\n",
    "The lower endpoint of the one-sided confidence interval is the smallest value of \n",
    "$q$ for which the\n",
    "corresponding test does not reject:\n",
    "\\begin{equation*}\n",
    "f(x) \\equiv \\min \\left \\{q \\in [0, 1]: \\sum_{k=x}^n \\binom{n}{k}q^k(1-q)^{n-k} \\ge \\alpha \\right \\}.\n",
    "\\end{equation*}\n",
    "Note that the upper tail probability, $\\sum_{k=x}^n \\binom{n}{k}q^k(1-q)^{n-k}$,\n",
    "increases continuously and \n",
    "monotonically as $q$ increases, so finding where it crosses $\\alpha$ is a straightforward\n",
    "root-finding problem.\n",
    "\n",
    "Let's code this up in python. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "from scipy.optimize import brentq  # Brent's root-finding algorithm\n",
    "from scipy.stats import binom      # the Binomial distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def binom_lower_ci(n: int, x: int, cl: float=0.95) -> float:\n",
    "    '''\n",
    "    lower confidence bound for a binomial p\n",
    "    \n",
    "    Assumes x is a draw from a binomial distribution with parameters\n",
    "    n (known) and p (unknown). Finds a lower confidence bound for p \n",
    "    at confidence level cl.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    n : int\n",
    "        number of trials, nonnegative integer\n",
    "    x : int\n",
    "        observed number of successes, nonnegative integer not larger than n\n",
    "    cl : float\n",
    "        confidence level, between 0 and 1\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    lb : float\n",
    "        lower confidence bound\n",
    "    '''\n",
    "    if x < 0 or x > n:\n",
    "        raise ValueError('impossible arguments')\n",
    "    if cl <= 0 or cl >= 1:\n",
    "        raise ValueError('silly confidence level')\n",
    "    lb = 0\n",
    "    if x > 0:\n",
    "        lb = brentq(lambda q: binom.sf(x-1,n,q)-(1-cl), 0, 1)\n",
    "    return lb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "n: 10\n",
      "p: 0.01 covered:  99.51%\n",
      "p: 0.05 covered:  98.82%\n",
      "p: 0.10 covered:  98.90%\n",
      "p: 0.20 covered:  96.73%\n",
      "p: 0.40 covered:  98.69%\n",
      "p: 0.50 covered:  99.07%\n",
      "p: 0.70 covered:  97.24%\n",
      "p: 0.90 covered:  100.00%\n",
      "p: 0.95 covered:  100.00%\n",
      "p: 0.99 covered:  100.00%\n",
      "\n",
      "n: 50\n",
      "p: 0.01 covered:  98.43%\n",
      "p: 0.05 covered:  96.22%\n",
      "p: 0.10 covered:  97.58%\n",
      "p: 0.20 covered:  96.80%\n",
      "p: 0.40 covered:  96.91%\n",
      "p: 0.50 covered:  96.76%\n",
      "p: 0.70 covered:  95.82%\n",
      "p: 0.90 covered:  96.87%\n",
      "p: 0.95 covered:  100.00%\n",
      "p: 0.99 covered:  100.00%\n",
      "\n",
      "n: 100\n",
      "p: 0.01 covered:  98.27%\n",
      "p: 0.05 covered:  97.33%\n",
      "p: 0.10 covered:  95.57%\n",
      "p: 0.20 covered:  96.56%\n",
      "p: 0.40 covered:  95.70%\n",
      "p: 0.50 covered:  95.96%\n",
      "p: 0.70 covered:  95.21%\n",
      "p: 0.90 covered:  97.75%\n",
      "p: 0.95 covered:  96.29%\n",
      "p: 0.99 covered:  100.00%\n"
     ]
    }
   ],
   "source": [
    "# some simulations to test the confidence intervals\n",
    "reps = int(10**4)\n",
    "for n in [10, 50, 100]:\n",
    "    print(f'\\nn: {n}')\n",
    "    for p in [0.01, 0.05, 0.1, 0.2, 0.4, 0.5, 0.7, 0.9, 0.95, 0.99]:\n",
    "        cover = 0\n",
    "        xs = binom.rvs(n, p, size=reps)\n",
    "        for x in xs:\n",
    "            cover += (1 if binom_lower_ci(n,x) <= p\n",
    "                 else 0)\n",
    "        print(f'p: {p:.2f} covered: {100*cover/reps : .2f}%')   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Think about how you might make a 2-sided confidence interval for $p$ instead of a lower 1-sided confidence interval. \n",
    "There are countless ways of constructing acceptance regions for the underlying tests, as mentioned in the [testing](./tests.ipynb) chapter. \n",
    "For instance, we could trim the same probability $\\alpha/2$ from both tails, or use the acceptance region that contains the fewest outcomes. \n",
    "The latter choice in general will lead to shorter confidence intervals.\n",
    "This approach is due to [Sterne (1954)](https://www.jstor.org/stable/2333026).\n",
    "\n",
    "\n",
    "Let's implement Sterne's approach now. The development is analogous to the randomized hypergeometric test in [the chapter on testing](./tests/ipynb), repeated below in this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import lru_cache\n",
    "\n",
    "@lru_cache(maxsize=None)  # decorate the function to cache the results \n",
    "                          # of calls to the function\n",
    "def binom_accept(n: int, p: float, alpha: float=0.05, randomized: bool=False) -> tuple:\n",
    "    '''\n",
    "    Acceptance region for a randomized binomial test\n",
    "    \n",
    "    If randomized==True, find the acceptance region for a randomized, exact \n",
    "    level-alpha test of the null hypothesis X~Binomial(n,p). \n",
    "    The acceptance region is the smallest possible. (And not, for instance, symmetric.)\n",
    "\n",
    "    If randomized==False, find the smallest conservative acceptance region.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    n : integer\n",
    "        number of independent trials\n",
    "    p : float\n",
    "        probability of success in each trial\n",
    "    alpha : float\n",
    "        desired significance level  \n",
    "    ramndomized : Boolean\n",
    "        return randomized exact test or conservative non-randomized test?\n",
    "  \n",
    "    Returns\n",
    "    --------\n",
    "    If randomized:\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",
    "    If not randomized:\n",
    "    I : list\n",
    "        values for which the test does not reject\n",
    "    \n",
    "    '''\n",
    "    if alpha <= 0 or alpha >= 1:\n",
    "        raise ValueError('bad significance level')\n",
    "    x = np.arange(0, n+1)\n",
    "    I = list(x)                    # start with all possible outcomes (then remove some)\n",
    "    pmf = binom.pmf(x,n,p)         # \"frozen\" binomial 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 outcomes for which test is randomized\n",
    "    p_tail = 0                     # probability of outcomes excluded from I\n",
    "    while p_tail < alpha:          # 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:                # remove outcomes from acceptance region\n",
    "            I.remove(j)\n",
    "    return_val = None\n",
    "    if randomized:\n",
    "        gamma = (p_tail-alpha)/p_J     # probability of accepting H_0 when X in J \n",
    "                                       # to get exact level alpha\n",
    "        return_val = I, J, gamma\n",
    "    else:\n",
    "        while p_tail > alpha:\n",
    "            j = J.pop()            # move the outcome into the acceptance region\n",
    "            p_tail -= pmf[j]\n",
    "            I.append(j)\n",
    "        return_val = I\n",
    "    return return_val \n",
    "\n",
    "\n",
    "def binom_ci(n: int, x: int, cl: float=0.95, eps: float=10**-3) -> tuple:\n",
    "    '''\n",
    "    two-sided confidence bound for a binomial p\n",
    "    \n",
    "    Assumes x is a draw from a binomial distribution with parameters\n",
    "    n (known) and p (unknown). Finds a confidence interval for p \n",
    "    at confidence level cl by inverting conservative tests\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    n : int\n",
    "        number of trials, nonnegative integer\n",
    "    x : int\n",
    "        observed number of successes, nonnegative integer not larger than n\n",
    "    cl : float\n",
    "        confidence level, between 1/2 and 1\n",
    "    eps : float in (0, 1)\n",
    "        resolution of the grid search\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    lb : float\n",
    "        lower confidence bound\n",
    "    ub : float\n",
    "        upper confidence bound\n",
    "    '''\n",
    "    if x < 0 or x > n:\n",
    "        raise ValueError('impossible arguments')\n",
    "    if cl <= 0 or cl >= 1:\n",
    "        raise ValueError('silly confidence level')\n",
    "    lb = 0\n",
    "    ub = 1\n",
    "    alpha = 1-cl\n",
    "    if x > 0:\n",
    "        while x not in binom_accept(n, lb, alpha, randomized=False):\n",
    "            lb += eps\n",
    "        lb -= eps\n",
    "    if x < n:\n",
    "        while x not in binom_accept(n, ub, alpha, randomized=False):\n",
    "            ub -= eps\n",
    "        ub += eps\n",
    "    return lb, ub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "n: 10\n",
      "p: 0.01 covered: 99.57% mean width: 0.31\n",
      "p: 0.05 covered: 98.89% mean width: 0.36\n",
      "p: 0.10 covered: 98.77% mean width: 0.41\n",
      "p: 0.20 covered: 96.87% mean width: 0.48\n",
      "p: 0.40 covered: 97.98% mean width: 0.54\n",
      "p: 0.50 covered: 97.90% mean width: 0.55\n",
      "p: 0.70 covered: 96.32% mean width: 0.52\n",
      "p: 0.90 covered: 98.71% mean width: 0.41\n",
      "p: 0.95 covered: 99.03% mean width: 0.36\n",
      "p: 0.99 covered: 99.66% mean width: 0.31\n",
      "\n",
      "n: 50\n",
      "p: 0.01 covered: 98.67% mean width: 0.09\n",
      "p: 0.05 covered: 95.91% mean width: 0.14\n",
      "p: 0.10 covered: 97.09% mean width: 0.17\n",
      "p: 0.20 covered: 95.18% mean width: 0.22\n",
      "p: 0.40 covered: 95.51% mean width: 0.27\n",
      "p: 0.50 covered: 96.75% mean width: 0.27\n",
      "p: 0.70 covered: 95.60% mean width: 0.25\n",
      "p: 0.90 covered: 96.70% mean width: 0.17\n",
      "p: 0.95 covered: 96.18% mean width: 0.13\n",
      "p: 0.99 covered: 98.50% mean width: 0.09\n",
      "\n",
      "n: 100\n",
      "p: 0.01 covered: 98.18% mean width: 0.05\n",
      "p: 0.05 covered: 96.50% mean width: 0.09\n",
      "p: 0.10 covered: 95.19% mean width: 0.12\n",
      "p: 0.20 covered: 95.30% mean width: 0.16\n",
      "p: 0.40 covered: 95.64% mean width: 0.19\n",
      "p: 0.50 covered: 96.60% mean width: 0.20\n",
      "p: 0.70 covered: 94.96% mean width: 0.18\n",
      "p: 0.90 covered: 95.29% mean width: 0.12\n",
      "p: 0.95 covered: 96.60% mean width: 0.09\n",
      "p: 0.99 covered: 98.19% mean width: 0.05\n"
     ]
    }
   ],
   "source": [
    "cl = 0.95\n",
    "eps = 10**-4\n",
    "plist = [0.01, 0.05, 0.1, 0.2, 0.4, 0.5, 0.7, 0.9, 0.95, 0.99]\n",
    "for n in [10, 50, 100]:\n",
    "    print(f'\\nn: {n}')\n",
    "    for p in plist:\n",
    "        cover = 0\n",
    "        mean_width = 0\n",
    "        xs = binom.rvs(n, p, size=reps)\n",
    "        for x in xs:\n",
    "            bounds = binom_ci(n, x, cl=cl, eps=eps)\n",
    "            cover += (1 if (bounds[0] <= p <= bounds[1])             \n",
    "                      else 0)\n",
    "            mean_width += bounds[1]-bounds[0]\n",
    "        mean_width /= reps\n",
    "        print(f'p: {p:.2f} covered: {100*cover/reps:.2f}% mean width: {mean_width:.2f}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bootstrap approximate confidence intervals\n",
    "\n",
    "The bootstrap approximates sampling from a population by sampling with replacement from a sample from the population. That is, it approximates the population distribution by the empirical distribution of the observed sample, which is\n",
    "the nonparametric maximum likelihood estimate of the population distribution.\n",
    "(See [notes on the bootstrap](./bootstrap.ipynb).)\n",
    "\n",
    "The bootstrap is commonly used to estimate the variability of an estimator of a functional parameter. It generally does a good job of estimating things like the variance of an estimator.\n",
    "\n",
    "It can also be used to construct approximate confidence intervals in a variety of ways, including the _percentile method_, which approximates percentiles of an estimator from percentiles of the resampling distribution of the estimator. However, this generally is not a good approximation.\n",
    "\n",
    "Let's code it up and see."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def boot_ci(data: list, estimator: callable, cl: float=0.95, reps: int=int(10**4), \n",
    "            random_state: object=None) -> tuple:\n",
    "    \"\"\"\n",
    "    Bootstrap percentile approximate confidence interval for a functional parameter\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    data : array-like\n",
    "        original sample\n",
    "    estimator : callable\n",
    "        function defined on the empirical distribution of a sample.\n",
    "        Applying it to the population distribution would give the\n",
    "        true value of the parameter of interest. Applying it to a \n",
    "        sample yields an estimator of the parameter of interest\n",
    "    cl : float in (0,1)\n",
    "        confidence level\n",
    "    reps : int, nonnegative\n",
    "        number of bootstrap samples to use\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    lb : float\n",
    "        estimated lower confidence bound\n",
    "    ub : float\n",
    "        estimated upper confidence bound\n",
    "    \"\"\"\n",
    "    if random_state is None:\n",
    "        prng = np.random.randomstate()\n",
    "    else:\n",
    "        prng = random_state\n",
    "    n = len(data)\n",
    "    estimates = []\n",
    "    for j in range(reps):\n",
    "        estimates.append(estimator(prng.choice(data, size=n, replace=True)))\n",
    "    estimates = np.array(estimates)\n",
    "    tail = (1-cl)/2\n",
    "    lower = np.quantile(estimates, tail)\n",
    "    upper = np.quantile(estimates, 1-tail)\n",
    "    return (lower, upper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "n: 10\n",
      "p: 0.01 covered:  9.20% mean width: 0.03\n",
      "p: 0.05 covered:  38.60% mean width: 0.13\n",
      "p: 0.10 covered:  62.20% mean width: 0.24\n",
      "p: 0.20 covered:  89.20% mean width: 0.42\n",
      "p: 0.40 covered:  93.80% mean width: 0.56\n",
      "p: 0.50 covered:  97.40% mean width: 0.58\n",
      "p: 0.70 covered:  96.40% mean width: 0.52\n",
      "p: 0.90 covered:  64.00% mean width: 0.24\n",
      "p: 0.95 covered:  40.50% mean width: 0.14\n",
      "p: 0.99 covered:  9.60% mean width: 0.03\n",
      "\n",
      "n: 50\n",
      "p: 0.01 covered:  36.60% mean width: 0.03\n",
      "p: 0.05 covered:  93.00% mean width: 0.11\n",
      "p: 0.10 covered:  96.70% mean width: 0.16\n",
      "p: 0.20 covered:  96.00% mean width: 0.22\n",
      "p: 0.40 covered:  95.90% mean width: 0.27\n",
      "p: 0.50 covered:  96.50% mean width: 0.27\n",
      "p: 0.70 covered:  95.50% mean width: 0.25\n",
      "p: 0.90 covered:  96.20% mean width: 0.16\n",
      "p: 0.95 covered:  90.90% mean width: 0.11\n",
      "p: 0.99 covered:  37.30% mean width: 0.03\n",
      "\n",
      "n: 100\n",
      "p: 0.01 covered:  63.80% mean width: 0.03\n",
      "p: 0.05 covered:  96.80% mean width: 0.08\n",
      "p: 0.10 covered:  94.90% mean width: 0.11\n",
      "p: 0.20 covered:  95.60% mean width: 0.15\n",
      "p: 0.40 covered:  96.00% mean width: 0.19\n",
      "p: 0.50 covered:  96.50% mean width: 0.19\n",
      "p: 0.70 covered:  95.30% mean width: 0.18\n",
      "p: 0.90 covered:  96.20% mean width: 0.11\n",
      "p: 0.95 covered:  96.10% mean width: 0.08\n",
      "p: 0.99 covered:  62.90% mean width: 0.03\n"
     ]
    }
   ],
   "source": [
    "reps = int(10**3)\n",
    "reps_boot = int(10**3)\n",
    "seed = 12345678\n",
    "prng = np.random.RandomState(seed)\n",
    "for n in [10, 50, 100]:\n",
    "    print(f'\\nn: {n}')\n",
    "    for p in [0.01, 0.05, 0.1, 0.2, 0.4, 0.5, 0.7, 0.9, 0.95, 0.99]:\n",
    "        cover = 0\n",
    "        mean_width = 0\n",
    "        xs = binom.rvs(n, p, size=reps)\n",
    "        for x in xs:\n",
    "            sample = [1]*x + [0]*(n-x)\n",
    "            bounds = boot_ci(sample, np.mean, reps=reps_boot, random_state=prng)\n",
    "            cover += (1 if (bounds[0] <= p <= bounds[1])             \n",
    "                      else 0)\n",
    "            mean_width += bounds[1]-bounds[0]\n",
    "        mean_width /= reps\n",
    "        print(f'p: {p:.2f} covered: {100*cover/reps : .2f}% mean width: {mean_width:.2f}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the coverage probability of bootstrap percentile confidence intervals can be \n",
    "much lower than their nominal level (here, as low as about 9% when they should be 95%). Where their coverage is about right, their average width is comparable to the average width of the conservative intervals derived by inverting two-sided tests.\n",
    "\n",
    "This example is a fairly simple situation: estimating the mean of a population of zeros and ones from a random sample with replacement. In more complicated examples, bootstrap confidence intervals generally do not attain their nominal level.\n",
    "There is a substantial body of work on how to improve the coverage of bootstrap CIs.\n",
    "_Pre-pivoting_ can help substantially. See the [chapter on the bootstrap](./bootstrap.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other approximate confidence intervals for Binomial($p$)\n",
    "\n",
    "There are many approximate methods for binomial confidence intervals, for instance, based on the normal approximation or the normal approximation to a transformation of the data.\n",
    "They generally do not attain their nominal confidence level, even for modest sample sizes,\n",
    "and their coverage probability varies erratically as $n$ and $p$ vary. See [Brown, Cai, and DasGupta (2001)](https://www.jstor.org/stable/2676784?seq=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 2: confidence interval for Hypergeometric $G$ (known $N$, $n$)\n",
    "\n",
    "Suppose we will observe $X \\sim \\mbox{Hyper}(N, G, n)$, where $N$ and $n$ are known but $G$ is not.\n",
    "We seek a one-sided lower confidence interval for $G$, that is, a set of the form\n",
    "$[f(X,U), \\infty)$ such that for all $G \\in \\{0, 1, \\ldots, N\\}$, if $X \\sim \\mbox{Hyper}(N, G, n)$\n",
    "and $U$ is an independent uniform random variable, then\n",
    "\\begin{equation*}\n",
    "\\mathbb{P} \\{ [f(X,U), \\infty) \\ni G \\} \\ge 1-\\alpha.\n",
    "\\end{equation*}\n",
    "Since it is certain that $G \\le N$, the upper endpoint of the interval can be reduced from $\\infty$ to $N$ without sacrificing coverage probability. That is, the same $f$ will satisfy\n",
    "\\begin{equation*}\n",
    "\\mathbb{P} \\{ [f(X,U), N] \\ni q \\} \\ge 1-\\alpha.\n",
    "\\end{equation*}\n",
    "\n",
    "To make such a lower confidence bound, we can invert one-sided hypothesis tests that\n",
    "reject when $X$ is \"too big.\"\n",
    "That is, we want a family of tests of the hypotheses $G = H$ for all $H \\in \\{0, 1, \\ldots, N\\}$\n",
    "that reject for large values of $X$. \n",
    "Such tests give evidence that $G$ is _at least_ a given size.\n",
    "\n",
    "To keep things simple, we will use conservative non-randomized tests rather than exact randomized tests. \n",
    "Because we are basing the confidence intervals on\n",
    "_conservative_ tests, we expect the coverage probability to be greater than $1-\\alpha$.\n",
    "We reject the hypothesis $G = I$ if, on the assumption that $G=I$, the chance that\n",
    "$X$ would be greater than or equal to its observed value is not greater than $\\alpha$.\n",
    "That is, \n",
    "\\begin{equation*}\n",
    "\\phi_q(x) = \\left \\{ \\begin{array}{ll}\n",
    "                      1, & \\sum_{k=x}^n \\frac{\\binom{I}{k}\\binom{N-I}{n-k}}{\\binom{N}{n}} \\ge \\alpha \\\\\n",
    "                      0, & \\mbox{otherwise.}\n",
    "                      \\end{array}\n",
    "            \\right .\n",
    "\\end{equation*}\n",
    "(Of course, we know $G \\ge X$.)\n",
    "The lower endpoint of the one-sided confidence interval is the smallest value of \n",
    "$I$ for which the\n",
    "corresponding test does not reject:\n",
    "\\begin{equation*}\n",
    "f(x) \\equiv \\min \\left \\{I \\in \\{x, x+1, \\ldots, N\\} : \\frac{\\binom{I}{k}\\binom{N-I}{n-k}}{\\binom{N}{n}} \\ge \\alpha \\right \\}.\n",
    "\\end{equation*}\n",
    "Note that the upper tail probability does not increase monotonically as $I$ increases;\n",
    "moreover, because $I$ must be an integer, the optimization problem is discrete,\n",
    "not continuous.\n",
    "Thus standard root-finding methods will _not_ let us find where the tail probability\n",
    "crosses $\\alpha$.\n",
    "Instead, we will use a search.\n",
    "\n",
    "Let's code this up in python. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import hypergeom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hypergeom_lower_ci(N: int, n: int, x: int, cl: float=0.95) -> float:\n",
    "    '''\n",
    "    lower confidence bound for a hypergeometric G\n",
    "    \n",
    "    Assumes x is a draw from a hypergeometric distribution with parameters\n",
    "    N (known), n (known), and G (unknown). Finds a lower confidence bound for G \n",
    "    at confidence level cl.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    N : int\n",
    "        population size, nonnegative integer\n",
    "    n : int\n",
    "        number of trials, nonnegative integer <= N\n",
    "    x : int\n",
    "        observed number of successes, nonnegative integer <= n\n",
    "    cl : float\n",
    "        confidence level, between 0 and 1\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    lb : float\n",
    "        lower confidence bound\n",
    "    '''\n",
    "    if x < 0 or x > n:\n",
    "        raise ValueError('impossible arguments')\n",
    "    if n > N or n < 1:\n",
    "        raise ValueError('impossible sample size')\n",
    "    if cl <= 0 or cl >= 1:\n",
    "        raise ValueError('silly confidence level')\n",
    "    lb = x\n",
    "    tail = hypergeom.sf(x-1, N, lb, n)\n",
    "    while tail < (1-cl):\n",
    "        lb += 1\n",
    "        tail = hypergeom.sf(x-1, N, lb, n)\n",
    "    return lb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N=20, G=10, n=10: covered  98.71%\n",
      "N=50, G=30, n=20: covered  98.22%\n",
      "N=100, G=25, n=30: covered  97.41%\n"
     ]
    }
   ],
   "source": [
    "NN = [20, 50, 100]\n",
    "nn = [10, 20, 30]\n",
    "GG = [10, 30, 25]\n",
    "reps = int(10**4)\n",
    "\n",
    "for j in range(len(NN)):\n",
    "    cover = 0\n",
    "    xs = hypergeom.rvs(NN[j], GG[j], nn[j], size=reps)\n",
    "    for x in xs:\n",
    "        cover += (1 if hypergeom_lower_ci(NN[j], nn[j], x, cl=0.95) <= GG[j]\n",
    "                 else 0)\n",
    "    print(f'N={NN[j]}, G={GG[j]}, n={nn[j]}: covered {100*cover/reps : .2f}%')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can make two-sided confidence bounds for $G$ by inverting two-sided tests, such as those given in the [chapter on testing](./tests.ipynb).\n",
    "Recall that we discussed two different constructions of two-sided tests there, one based on trimming the same probability from each tail, and one based on minimizing the number of outcomes in the acceptance region."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confidence intervals from permutation tests\n",
    "\n",
    "### The two-sample problem\n",
    "\n",
    "Recall from the chapter on [hypothesis tests](./tests.ipynb) that the two-sample problem asks whether two groups plausibly resulted from allocating their union randomly into two groups of their observed sizes.\n",
    "\n",
    "That is, we have two groups of data, $\\{x_j\\}_{j=1}^n$ and $\\{y_j\\}_{j=1}^m$, and hypothesize that they arose by taking the multiset of $n+m$ values $\\{x_1, \\ldots, x_n, y_1, \\ldots, y_m\\}$ and randomly selecting $n$ items to comprise the first group, with the remaining $m$ comprising the second group.\n",
    "\n",
    "This problem arises in many contexts, including randomized controlled trials:\n",
    "We have a group of $n+m$ subjects of whom $n$ are selected at random to be the control/placebo group\n",
    "and the other $m$ receive the active treatment.\n",
    "The _strong null hypothesis_ of no treatment effect is that each subject would have had the same response, no matter which treatment the subject was assigned.\n",
    "It is as if the response were determined before the assignment occurred. \n",
    "The assignment just reveals the response corresponding to the assigned treatment.\n",
    "\n",
    "Thus, testing the strong null hypothesis is an instance of the two-sample problem: did the two sets of values plausibly arise from partitioning a single set of values into two groups by simple random sampling?"
   ]
  },
  {
   "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
}