normApprox.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Confidence bounds and tests for the population mean, Part I: \n",
    "# The Normal Approximation and Bootstrap Percentile Intervals\n",
    "\n",
    "This notebook examines confidence bounds and tests for a population mean from a random sample,\n",
    "using the normal approximation or the bootstrap.\n",
    "\n",
    "There are many situations in which an exact or conservative (not approximate) confidence interval or test is required, and there are many populations for which asymptotic/approximate tests and confidence intervals are very inaccurate, especially for modest sample sizes.\n",
    "\n",
    "The context for what follows is \"nonstandard\" problems where the population distribution may be quite skewed, and\n",
    "sample sizes may be small.\n",
    "\n",
    "For instance, in applications to financial audits, populations of overstatement errors are often mostly\n",
    "zero or small (i.e., the account is accurately reported), with possible occasional large values.\n",
    "The same is true for election audits based on comparing manual counts of votes to machine counts of votes.\n",
    "\n",
    "This notebook shows that standard approximate methods can fail dramatically in those situations.\n",
    "Later notebooks develop exact and conservative confidence sets and tests that work in those situations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confidence intervals based on the normal approximation\n",
    "\n",
    "One of the first things most statisticians would consider is to use confidence based on the normal approximation, especially if it were possible to draw a moderately large sample (e.g., 100 or more units).\n",
    "\n",
    "The procedure would be:\n",
    "+ choose the desired coverage probability $1-\\alpha$ and the sample size $n$\n",
    "+ draw a simple random sample or a random sample with replacement of size $n$ from the population\n",
    "+ compute the sample mean $\\bar{X}$ and sample standard deviation $s$\n",
    "+ find $\\hat{SE} = f \\times s/\\sqrt{n}$, the estimated standard error of the sample mean, where $f=1$ for sampling with replacement, and $f = \\sqrt{\\frac{N-n}{N-1}}$ for sampling without replacement ($N$ is the\n",
    "population size)\n",
    "+ find $t$, the appropriate critical value of Student's t distribution with $n-1$ degrees of freedom\n",
    "+ the approximate 2-sided confidence interval is $[\\bar{X} - t \\hat{SE}, \\bar{X} + t \\hat{SE}]$.\n",
    "\n",
    "### _Formally_, this gives an approximate confidence interval, but there's no guarantee that the actual coverage probability is close to the nominal coverage probability\n",
    "\n",
    "+ Normal approximation is common, but can be grossly inaccurate, especially with:\n",
    "  - small sample sizes\n",
    "  - skewed populations\n",
    "  \n",
    "In auditing situations, typical for the distributions to be very skewed: point mass at zero, mixed with something else.\n",
    "\n",
    "+ Examine scenarios by simulation.\n",
    "+ First scenario is a population of zeros and ones, in various proportions\n",
    "+ Second scenario is \"nonstandard\" mixture of pointmass at 0 and a Uniform\n",
    "\n",
    "We are estimating the coverage by simulation. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### General considerations for simulations\n",
    "\n",
    "+ It's important to know which pseudo-random number generator (PRNG) you are using (Python uses the Mersenne Twister by default)\n",
    "+ It's important to assess whether that PRNG is adequate  (the Mersenne twister is pretty good for statistical simulations, but not good enough for cryptographic applications). See, e.g., [Stark and Ottoboni, 2018](https://arxiv.org/abs/1810.10985)\n",
    "+ It's important to remember that your estimates will have sampling error\n",
    "+ If the result matters, you should quantify the sampling error, for instance, by a confidence interval for the estimated result.\n",
    "\n",
    "Here, we will use 1,000 or 10,000 iterations in the estimates. \n",
    "For estimating probabilities, the standard error of the result is no larger than $1/(2\\sqrt{n})$: 1.58% for $n=1000$ and 0.5% for $n=10,000$.\n",
    "\n",
    "For estimating things like expected length, guarantees of accuracy are harder, but the kinds of techniques we will see can be applied to _that_ problem too."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bootstrap Percentile 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",
    "\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. \n",
    "However, this approximation generally is not accurate..\n",
    "\n",
    "Let's code it up and see."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# First code cell: set up the Python environment\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.stats\n",
    "from scipy.stats import binom\n",
    "import pandas as pd\n",
    "# For interactive widgets\n",
    "from ipywidgets import interact, interactive, fixed\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import clear_output, display, HTML"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def boot_percentile_ci(data : list, estimator : callable, cl : float=0.95, reps : int=int(10**4), \\\n",
    "                       random_state=None):\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": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Simulated coverage probability and expected length of Student-t and bootstrap percentile confidence intervals for a {0, 1} population</h3><strong>Nominal coverage probability 95.0%</strong>.<br /><strong>Estimated from 10000 replications.</strong>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>fraction of 1s</th>\n",
       "      <th>sample size</th>\n",
       "      <th>Student-t cov</th>\n",
       "      <th>Student-t len</th>\n",
       "      <th>Bootstrap cov</th>\n",
       "      <th>Bootstrap len</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.001</td>\n",
       "      <td>25</td>\n",
       "      <td>2.50%</td>\n",
       "      <td>0.0042</td>\n",
       "      <td>2.50%</td>\n",
       "      <td>0.0031</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.001</td>\n",
       "      <td>50</td>\n",
       "      <td>4.60%</td>\n",
       "      <td>0.0037</td>\n",
       "      <td>4.60%</td>\n",
       "      <td>0.0028</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.001</td>\n",
       "      <td>100</td>\n",
       "      <td>9.72%</td>\n",
       "      <td>0.0039</td>\n",
       "      <td>9.72%</td>\n",
       "      <td>0.0031</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.001</td>\n",
       "      <td>400</td>\n",
       "      <td>32.05%</td>\n",
       "      <td>0.0034</td>\n",
       "      <td>31.96%</td>\n",
       "      <td>0.0028</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.010</td>\n",
       "      <td>25</td>\n",
       "      <td>22.13%</td>\n",
       "      <td>0.0385</td>\n",
       "      <td>22.11%</td>\n",
       "      <td>0.0292</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.010</td>\n",
       "      <td>50</td>\n",
       "      <td>39.74%</td>\n",
       "      <td>0.0352</td>\n",
       "      <td>39.62%</td>\n",
       "      <td>0.0282</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.010</td>\n",
       "      <td>100</td>\n",
       "      <td>63.13%</td>\n",
       "      <td>0.0303</td>\n",
       "      <td>63.13%</td>\n",
       "      <td>0.0258</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.010</td>\n",
       "      <td>400</td>\n",
       "      <td>90.56%</td>\n",
       "      <td>0.0188</td>\n",
       "      <td>90.90%</td>\n",
       "      <td>0.0180</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.100</td>\n",
       "      <td>25</td>\n",
       "      <td>92.83%</td>\n",
       "      <td>0.2337</td>\n",
       "      <td>92.05%</td>\n",
       "      <td>0.2058</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.100</td>\n",
       "      <td>50</td>\n",
       "      <td>87.60%</td>\n",
       "      <td>0.1668</td>\n",
       "      <td>95.57%</td>\n",
       "      <td>0.1573</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>0.100</td>\n",
       "      <td>100</td>\n",
       "      <td>93.20%</td>\n",
       "      <td>0.1180</td>\n",
       "      <td>95.68%</td>\n",
       "      <td>0.1145</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>0.100</td>\n",
       "      <td>400</td>\n",
       "      <td>95.01%</td>\n",
       "      <td>0.0588</td>\n",
       "      <td>95.38%</td>\n",
       "      <td>0.0582</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    fraction of 1s  sample size Student-t cov Student-t len Bootstrap cov  \\\n",
       "0            0.001           25         2.50%        0.0042         2.50%   \n",
       "1            0.001           50         4.60%        0.0037         4.60%   \n",
       "2            0.001          100         9.72%        0.0039         9.72%   \n",
       "3            0.001          400        32.05%        0.0034        31.96%   \n",
       "4            0.010           25        22.13%        0.0385        22.11%   \n",
       "5            0.010           50        39.74%        0.0352        39.62%   \n",
       "6            0.010          100        63.13%        0.0303        63.13%   \n",
       "7            0.010          400        90.56%        0.0188        90.90%   \n",
       "8            0.100           25        92.83%        0.2337        92.05%   \n",
       "9            0.100           50        87.60%        0.1668        95.57%   \n",
       "10           0.100          100        93.20%        0.1180        95.68%   \n",
       "11           0.100          400        95.01%        0.0588        95.38%   \n",
       "\n",
       "   Bootstrap len  \n",
       "0         0.0031  \n",
       "1         0.0028  \n",
       "2         0.0031  \n",
       "3         0.0028  \n",
       "4         0.0292  \n",
       "5         0.0282  \n",
       "6         0.0258  \n",
       "7         0.0180  \n",
       "8         0.2058  \n",
       "9         0.1573  \n",
       "10        0.1145  \n",
       "11        0.0582  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Population of two values, {0, 1}, in various proportions.  Amounts to Binomial random variable\n",
    "seed = 12345678\n",
    "prng = np.random.RandomState(seed)\n",
    "\n",
    "ns = np.array([25, 50, 100, 400])  # sample sizes\n",
    "ps = np.array([.001, .01, 0.1])    # mixture fractions, proportion of 1s in the population\n",
    "alpha = 0.05  # 1- (confidence level)\n",
    "reps = int(10**4) # just for demonstration\n",
    "reps_boot = int(10**3)\n",
    "vals = [0, 1]\n",
    "\n",
    "simTable = pd.DataFrame(columns=('fraction of 1s', 'sample size', 'Student-t cov', 'Student-t len', \\\n",
    "                                 'Bootstrap cov', 'Bootstrap len'))\n",
    "for p in ps:\n",
    "    popMean = p\n",
    "    for n in ns:\n",
    "        tCrit = sp.stats.t.ppf(q=1-alpha/2, df=n-1)\n",
    "        sam = sp.stats.binom.rvs(n, p, size=reps)\n",
    "        cover_boot = 0\n",
    "        aveLen_boot = 0\n",
    "        for x in sam:\n",
    "            sample = [1]*x + [0]*(n-x)\n",
    "            bounds = boot_percentile_ci(sample, np.mean, cl=1-alpha, reps=reps_boot, random_state=prng)\n",
    "            cover_boot += (1 if (bounds[0] <= p <= bounds[1])             \n",
    "                           else 0)\n",
    "            aveLen_boot += bounds[1]-bounds[0]\n",
    "        cover_boot = cover_boot/reps\n",
    "        aveLen_boot = aveLen_boot/reps\n",
    "        samMean = sam/n\n",
    "        samSD = np.sqrt(samMean*(1-samMean)/(n-1))\n",
    "        cover = ( np.fabs(samMean-popMean) < tCrit*samSD).mean()\n",
    "        aveLen = 2*(tCrit*samSD).mean()\n",
    "        simTable.loc[len(simTable)] =  p, n, f'{100*cover :0.2f}%', f'{aveLen :0.4f}', \\\n",
    "                                       f'{100*cover_boot :0.2f}%', f'{aveLen_boot :0.4f}'\n",
    "#\n",
    "ansStr =  '<h3>Simulated coverage probability and expected length of Student-t and bootstrap percentile confidence intervals for a {0, 1} population</h3>' +\\\n",
    "          '<strong>Nominal coverage probability ' + str(100*(1-alpha)) +\\\n",
    "          '%</strong>.<br /><strong>Estimated from ' + str(reps) + ' replications.</strong>'\n",
    "\n",
    "display(HTML(ansStr))\n",
    "display(simTable)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why/When is the coverage bad?\n",
    "\n",
    "+ Because of the high population mass at 0, the sample will often contain _only_ 0s. Then the SD will be zero, the length of the Student-t interval will be 0, and the interval won't cover the mean&mdash;which is bigger than zero.\n",
    "Similarly, the bootstrap will be resampling from a sample that consists entirely of zeros.\n",
    "+ In general, coverage will be bad when sampling distribution of the sample mean is strongly skewed\n",
    "\n",
    "E.g., Binomial with very small or very large $p$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7c81007302e64033981ce2c95b814c35",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=30, description='n', max=300, min=5, step=5), FloatSlider(value=0.001, d…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plotBinomialPmf(n, p)>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def plotBinomialPmf(n, p):\n",
    "    '''\n",
    "       Plots the sampling distribution of the sample mean of a sample of size n\n",
    "       drawn with replacement from a {0, 1} population; i.e., a scaled binomial.\n",
    "       Plots a vertical red line at the expected value.\n",
    "    '''\n",
    "    fig, ax = plt.subplots(nrows=1, ncols=1)\n",
    "    x = np.arange(n+1)\n",
    "    xFrac = x.astype(float)/ x.max()\n",
    "    width = 1.0/n\n",
    "    ax.bar(xFrac, binom.pmf(x, n, p), width, align=\"center\")\n",
    "    plt.xlim([-width,1+width])\n",
    "    plt.axvline(x=p, color='r')\n",
    "    plt.show()\n",
    "\n",
    "interact(plotBinomialPmf, n=widgets.IntSlider(min=5, max=300, step=5, value=30),\\\n",
    "         p=widgets.FloatSlider(min=0.001, max=1, step=0.001, value=.001))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nonstandard mixtures of distributions\n",
    "\n",
    "+ Common in statistical study of auditing to consider \"nonstandard mixtures\" of a pointmass and continuous distribution \n",
    "+ E.g., NRC Monograph and [1989 Statistical Science article](http://projecteuclid.org/download/pdf_1/euclid.ss/1177012665)\n",
    "+ Let's look at a mixture of $U[0,1]$ and a pointmass at zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Simulated coverage probability and expected length of Student-t confidence intervals for mixture of U[0,1] and pointmass at 0.</h3><h4>Nominal coverage probability 95.0%.</h4> Estimated from 10000 replications."
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mass at 0</th>\n",
       "      <th>sample size</th>\n",
       "      <th>Student-t cov</th>\n",
       "      <th>Student-t len</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.900</td>\n",
       "      <td>25</td>\n",
       "      <td>90.31%</td>\n",
       "      <td>0.6428</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.900</td>\n",
       "      <td>50</td>\n",
       "      <td>98.95%</td>\n",
       "      <td>0.6719</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.900</td>\n",
       "      <td>100</td>\n",
       "      <td>99.96%</td>\n",
       "      <td>0.68</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.900</td>\n",
       "      <td>400</td>\n",
       "      <td>100.0%</td>\n",
       "      <td>0.6869</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.990</td>\n",
       "      <td>25</td>\n",
       "      <td>22.61%</td>\n",
       "      <td>0.1011</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.990</td>\n",
       "      <td>50</td>\n",
       "      <td>39.54%</td>\n",
       "      <td>0.1297</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.990</td>\n",
       "      <td>100</td>\n",
       "      <td>62.67%</td>\n",
       "      <td>0.1627</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.990</td>\n",
       "      <td>400</td>\n",
       "      <td>98.05%</td>\n",
       "      <td>0.2098</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.999</td>\n",
       "      <td>25</td>\n",
       "      <td>2.26%</td>\n",
       "      <td>0.0092</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.999</td>\n",
       "      <td>50</td>\n",
       "      <td>4.95%</td>\n",
       "      <td>0.0145</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>0.999</td>\n",
       "      <td>100</td>\n",
       "      <td>9.43%</td>\n",
       "      <td>0.0191</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>0.999</td>\n",
       "      <td>400</td>\n",
       "      <td>32.97%</td>\n",
       "      <td>0.036</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    mass at 0 sample size Student-t cov Student-t len\n",
       "0       0.900          25        90.31%        0.6428\n",
       "1       0.900          50        98.95%        0.6719\n",
       "2       0.900         100        99.96%          0.68\n",
       "3       0.900         400        100.0%        0.6869\n",
       "4       0.990          25        22.61%        0.1011\n",
       "5       0.990          50        39.54%        0.1297\n",
       "6       0.990         100        62.67%        0.1627\n",
       "7       0.990         400        98.05%        0.2098\n",
       "8       0.999          25         2.26%        0.0092\n",
       "9       0.999          50         4.95%        0.0145\n",
       "10      0.999         100         9.43%        0.0191\n",
       "11      0.999         400        32.97%         0.036"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Nonstandard mixture: a pointmass at zero and a uniform[0,1]\n",
    "ns = np.array([25, 50, 100, 400])  # sample sizes\n",
    "ps = np.array([0.9, 0.99, 0.999])    # mixture fraction, weight of pointmass\n",
    "alpha = 0.05  # 1- (confidence level)\n",
    "reps = int(1.0e4) # just for demonstration\n",
    "\n",
    "simTable = pd.DataFrame(columns=('mass at 0', 'sample size', 'Student-t cov', 'Student-t len'))\n",
    "for p in ps:\n",
    "    popMean = (1-p)*0.5  #  p*0 + (1-p)*.5\n",
    "    for n in ns:\n",
    "        tCrit = sp.stats.t.ppf(q=1-alpha/2, df=n-1)\n",
    "        cover = 0\n",
    "        totLen = 0.0\n",
    "        samMean = np.zeros(reps)\n",
    "        for rep in range(int(reps)):\n",
    "            sam = np.random.uniform(size=n)\n",
    "            ptMass = np.random.uniform(size=n)\n",
    "            sam[ptMass < p] = 0.0\n",
    "            samMean[rep] = np.mean(sam)\n",
    "            samSD = np.std(sam, ddof=1)\n",
    "            cover += ( math.fabs(samMean[rep]-popMean) < tCrit*samSD )\n",
    "            totLen += 2*tCrit*samSD\n",
    "        simTable.loc[len(simTable)] =  p, n, str(100*float(cover)/float(reps)) + '%', str(round(totLen/reps,4))\n",
    "#\n",
    "ansStr =  '<h3>Simulated coverage probability and expected length of Student-t confidence intervals for mixture of U[0,1] and ' +\\\n",
    "          'pointmass at 0.</h3>' +\\\n",
    "          '<h4>Nominal coverage probability ' + str(100*(1-alpha)) +\\\n",
    "          '%.</h4> Estimated from ' + str(reps) + ' replications.'\n",
    "\n",
    "display(HTML(ansStr))\n",
    "display(simTable)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again the distribution of the sample mean is skewed, which hurts the coverage of Student-t intervals. \n",
    "\n",
    "When the mass at 0 gets large enough, there is again a large chance that the sample will contain _only_ 0s, in which case the sample SD will be zero, the length of the interval will be zero, and the interval won't cover."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "* The actual coverage probability of confidence intervals based on the normal approximation or the bootstrap percentile method can be *much* less than the nominal coverage probability, even for large-ish sample sizes, for skewed population distributions.\n",
    "* In many applications, the population _is_ very skewed: \"Nonstandard mixture distributions\"\n",
    "* For example, accounting overstatements might be mostly 0, and occasionally large. Then:\n",
    "   * Can be very likely that the sample variance is zero\n",
    "   * Leads to apparent certainty that the result is correct."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What's next?\n",
    "We will consider some conservative methods for constructing confidence bounds by inverting hypothesis tests.\n",
    "\n",
    "- [Next: Duality between hypothesis tests and confidence sets](duality.ipynb)\n",
    "- [Previous: Canonical examples of real-world problems we will consider](canonical.ipynb)\n",
    "- [Index](index.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}