{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Jeffreys interval\n",
    "\n",
    "Copyright 2020 Allen Downey\n",
    "\n",
    "[MIT License](https://opensource.org/licenses/MIT)\n",
    "\n",
    "\n",
    "Suppose you have run 174 trials and 173 were successful.  You want to report an estimate of the probability of success and a confidence interval for the estimate.\n",
    "\n",
    "According to our friends at [Wikipedia](https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval) there are several ways to compute it, based on different assumptions and requirements.\n",
    "\n",
    "In my opinion, the clear best option is the [Jeffreys interval](https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval#Jeffreys_interval\n",
    ")."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">The Jeffreys interval has a Bayesian derivation, but it has good frequentist properties. In particular, it has coverage properties that are similar to those of the Wilson interval, but it is one of the few intervals with the advantage of being equal-tailed (e.g., for a 95% confidence interval, the probabilities of the interval lying above or below the true value are both close to 2.5%). In contrast, the Wilson interval has a systematic bias such that it is centred too close to p = 0.5.\n",
    ">\n",
    ">The Jeffreys interval is the Bayesian credible interval obtained when using the non-informative Jeffreys prior for the binomial proportion p. The Jeffreys prior for this problem is a Beta distribution with parameters (1/2, 1/2), it is a conjugate prior. After observing x successes in n trials, the posterior distribution for p is a Beta distribution with parameters (x + 1/2, n – x + 1/2).\n",
    ">\n",
    ">When x ≠0 and x ≠ n, the Jeffreys interval is taken to be the 100(1 – α)% equal-tailed posterior probability interval, i.e., the α / 2 and 1 – α / 2 quantiles of a Beta distribution with parameters (x + 1/2, n – x + 1/2). These quantiles need to be computed numerically, although this is reasonably simple with modern statistical software.\n",
    ">\n",
    ">In order to avoid the coverage probability tending to zero when p → 0 or 1, when x = 0 the upper limit is calculated as before but the lower limit is set to 0, and when x = n the lower limit is calculated as before but the upper limit is set to 1.\n",
    "\n",
    "In my opinion, that sentence is an unnecessary hack.\n",
    "\n",
    "Here's how to compute a Jeffrey's interval for the example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trials = 174\n",
    "successes = 173\n",
    "failures = trials-successes\n",
    "failures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a beta distribution that represents the posterior distribution for the proportion, assuming a Jeffrey's prior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import beta\n",
    "\n",
    "dist = beta(successes+1/2, failures+1/2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I think the best point estimate is the posterior mean:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "99.14285714285714"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "estimate = dist.mean() * 100\n",
    "estimate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the confidence interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.050000000000000044, 0.95)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = 0.95\n",
    "a = (1-p)\n",
    "a, 1-a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([97.34566219, 99.9379198 ])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ci = dist.ppf([a/2, 1-a/2]) * 100\n",
    "ci"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So you could report an estimate of 99.1% with 95% CI (97.3, 99.9).\n",
    "\n",
    "If anyone asks how you computed that, tell them it's a Jeffreys interval.  If they ask why, send them [this paper](https://repository.upenn.edu/cgi/viewcontent.cgi?article=1440&context=statistics_papers)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}