likelihood_ratio_process.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='likelihood-ratio-process'></a>\n",
    "<div id=\"qe-notebook-header\" align=\"right\" style=\"text-align:right;\">\n",
    "        <a href=\"https://quantecon.org/\" title=\"quantecon.org\">\n",
    "                <img style=\"width:250px;display:inline;\" width=\"250px\" src=\"https://assets.quantecon.org/img/qe-menubar-logo.svg\" alt=\"QuantEcon\">\n",
    "        </a>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Likelihood Ratio Processes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Likelihood Ratio Processes](#Likelihood-Ratio-Processes)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [Likelihood Ratio Process](#Likelihood-Ratio-Process)  \n",
    "  - [Nature Permanently Draws from Density g](#Nature-Permanently-Draws-from-Density-g)  \n",
    "  - [Nature Permanently Draws from Density f](#Nature-Permanently-Draws-from-Density-f)  \n",
    "  - [Likelihood Ratio Test](#Likelihood-Ratio-Test)  \n",
    "  - [Sequels](#Sequels)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numba import vectorize, njit\n",
    "from math import gamma\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "This lecture describes likelihood ratio processes and some of their uses.\n",
    "\n",
    "We’ll use the simple statistical setting also used in [this lecture](https://python-programming.quantecon.org/exchangeable.html).\n",
    "\n",
    "Among the things that we’ll learn about are\n",
    "\n",
    "> - A peculiar property of likelihood ratio processes  \n",
    "- How a likelihood ratio process is the key ingredient in frequentist hypothesis testing  \n",
    "- How a **receiver operator characteristic curve** summarizes information about a false alarm probability and power in frequentist hypothesis testing  \n",
    "- How during World War II the United States Navy devised a decision rule that Captain Garret L. Schyler challenged and asked Milton Friedman to justify to him, a topic to be studied in  [this lecture](https://python-programming.quantecon.org/wald_friedman.html)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Likelihood Ratio Process\n",
    "\n",
    "A nonnegative random variable $ W $ has one of two probability density functions, either\n",
    "$ f $ or $ g $.\n",
    "\n",
    "Before the beginning of time, nature once and for all decides whether she will draw a sequence of IID draws from either\n",
    "$ f $ or $ g $.\n",
    "\n",
    "We will sometimes let $ q $ be the density that nature chose once and for all, so\n",
    "that $ q $ is either $ f $ or $ g $, permanently.\n",
    "\n",
    "Nature knows which density it permanently draws from, but we the observers do not.\n",
    "\n",
    "We do know both $ f $ and $ g $ but we don’t know which density nature\n",
    "chose.\n",
    "\n",
    "But we want to know.\n",
    "\n",
    "To do that, we use observations.\n",
    "\n",
    "We observe a sequence $ \\{w_t\\}_{t=1}^T $ of $ T $ IID draws\n",
    "from either $ f $ or $ g $.\n",
    "\n",
    "We want to use these observations to infer whether nature chose $ f $ or\n",
    "$ g $.\n",
    "\n",
    "A **likelihood ratio process** is a useful tool for this task.\n",
    "\n",
    "To begin, we define key component of a likelihood ratio process, namely, the time $ t $ likelihood ratio  as the random variable\n",
    "\n",
    "$$\n",
    "\\ell (w_t)=\\frac{f\\left(w_t\\right)}{g\\left(w_t\\right)},\\quad t\\geq1.\n",
    "$$\n",
    "\n",
    "We assume that $ f $ and $ g $ both put positive probabilities on the\n",
    "same intervals of possible realizations of the random variable $ W $.\n",
    "\n",
    "That means that under the $ g $ density,  $ \\ell (w_t)=\n",
    "\\frac{f\\left(w_{t}\\right)}{g\\left(w_{t}\\right)} $\n",
    "is evidently a nonnegative  random variable with mean $ 1 $.\n",
    "\n",
    "A **likelihood ratio process** for sequence\n",
    "$ \\left\\{ w_{t}\\right\\} _{t=1}^{\\infty} $ is defined as\n",
    "\n",
    "$$\n",
    "L\\left(w^{t}\\right)=\\prod_{i=1}^{t} \\ell (w_i),\n",
    "$$\n",
    "\n",
    "where $ w^t=\\{ w_1,\\dots,w_t\\} $ is a history of\n",
    "observations up to and including time $ t $.\n",
    "\n",
    "Sometimes for shorthand we’ll write $ L_t =  L(w^t) $.\n",
    "\n",
    "Notice that the likelihood process satisfies the *recursion* or\n",
    "*multiplicative decomposition*\n",
    "\n",
    "$$\n",
    "L(w^t) = \\ell (w_t) L (w^{t-1}) .\n",
    "$$\n",
    "\n",
    "The likelihood ratio and its logarithm are key tools for making\n",
    "inferences using a classic frequentist approach due to Neyman and\n",
    "Pearson [[NP33]](https://python-programming.quantecon.org/zreferences.html#neyman-pearson).\n",
    "\n",
    "To help us appreciate how things work, the following Python code evaluates $ f $ and $ g $ as two different\n",
    "beta distributions, then computes and simulates an associated likelihood\n",
    "ratio process by generating a sequence $ w^t $ from *some*\n",
    "probability distribution, for example, a sequence of  IID draws from $ g $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Parameters in the two beta distributions.\n",
    "F_a, F_b = 1, 1\n",
    "G_a, G_b = 3, 1.2\n",
    "\n",
    "@vectorize\n",
    "def p(x, a, b):\n",
    "    r = gamma(a + b) / (gamma(a) * gamma(b))\n",
    "    return r * x** (a-1) * (1 - x) ** (b-1)\n",
    "\n",
    "# The two density functions.\n",
    "f = njit(lambda x: p(x, F_a, F_b))\n",
    "g = njit(lambda x: p(x, G_a, G_b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@njit\n",
    "def simulate(a, b, T=50, N=500):\n",
    "    '''\n",
    "    Generate N sets of T observations of the likelihood ratio,\n",
    "    return as N x T matrix.\n",
    "\n",
    "    '''\n",
    "\n",
    "    l_arr = np.empty((N, T))\n",
    "\n",
    "    for i in range(N):\n",
    "\n",
    "        for j in range(T):\n",
    "            w = np.random.beta(a, b)\n",
    "            l_arr[i, j] = f(w) / g(w)\n",
    "\n",
    "    return l_arr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nature Permanently Draws from Density g\n",
    "\n",
    "We first simulate the likelihood ratio process when nature permanently\n",
    "draws from $ g $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "l_arr_g = simulate(G_a, G_b)\n",
    "l_seq_g = np.cumprod(l_arr_g, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "N, T = l_arr_g.shape\n",
    "\n",
    "for i in range(N):\n",
    "\n",
    "    plt.plot(range(T), l_seq_g[i, :], color='b', lw=0.8, alpha=0.5)\n",
    "\n",
    "plt.ylim([0, 3])\n",
    "plt.title(\"$L(w^{t})$ paths\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evidently, as sample length $ T $ grows, most probability mass\n",
    "shifts toward zero\n",
    "\n",
    "To see it this more clearly clearly, we plot over time the fraction of\n",
    "paths $ L\\left(w^{t}\\right) $ that fall in the interval\n",
    "$ \\left[0, 0.01\\right] $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plt.plot(range(T), np.sum(l_seq_g <= 0.01, axis=0) / N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Despite the evident convergence of most probability mass to a\n",
    "very small interval near $ 0 $,  the unconditional mean of\n",
    "$ L\\left(w^t\\right) $ under probability density $ g $ is\n",
    "identically $ 1 $ for all $ t $.\n",
    "\n",
    "To verify this assertion, first notice that as mentioned earlier the unconditional mean\n",
    "$ E_{0}\\left[\\ell \\left(w_{t}\\right)\\bigm|q=g\\right] $ is $ 1 $ for\n",
    "all $ t $:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "E_{0}\\left[\\ell \\left(w_{t}\\right)\\bigm|q=g\\right]  &=\\int\\frac{f\\left(w_{t}\\right)}{g\\left(w_{t}\\right)}g\\left(w_{t}\\right)dw_{t} \\\\\n",
    "    &=\\int f\\left(w_{t}\\right)dw_{t} \\\\\n",
    "    &=1,\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "which immediately implies\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "E_{0}\\left[L\\left(w^{1}\\right)\\bigm|q=g\\right]  &=E_{0}\\left[\\ell \\left(w_{1}\\right)\\bigm|q=g\\right]\\\\\n",
    "    &=1.\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Because $ L(w^t) = \\ell(w_t) L(w^{t-1}) $ and\n",
    "$ \\{w_t\\}_{t=1}^t $ is an IID sequence, we have\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "E_{0}\\left[L\\left(w^{t}\\right)\\bigm|q=g\\right]  &=E_{0}\\left[L\\left(w^{t-1}\\right)\\ell \\left(w_{t}\\right)\\bigm|q=g\\right] \\\\\n",
    "    &=E_{0}\\left[L\\left(w^{t-1}\\right)E\\left[\\ell \\left(w_{t}\\right)\\bigm|q=g,w^{t-1}\\right]\\bigm|q=g\\right] \\\\\n",
    "    &=E_{0}\\left[L\\left(w^{t-1}\\right)E\\left[\\ell \\left(w_{t}\\right)\\bigm|q=g\\right]\\bigm|q=g\\right] \\\\\n",
    "    &=E_{0}\\left[L\\left(w^{t-1}\\right)\\bigm|q=g\\right] \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "for any $ t \\geq 1 $.\n",
    "\n",
    "Mathematical induction implies\n",
    "$ E_{0}\\left[L\\left(w^{t}\\right)\\bigm|q=g\\right]=1 $ for all\n",
    "$ t \\geq 1 $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Peculiar Property of Likelihood Ratio Process\n",
    "\n",
    "How can $ E_{0}\\left[L\\left(w^{t}\\right)\\bigm|q=g\\right]=1 $ possibly be true when most  probability mass of the likelihood\n",
    "ratio process is piling up near $ 0 $ as\n",
    "$ t \\rightarrow + \\infty $?\n",
    "\n",
    "The answer has to be that as $ t \\rightarrow + \\infty $, the\n",
    "distribution of $ L_t $ becomes more and more fat-tailed:\n",
    "enough  mass shifts to larger and larger values of $ L_t $ to make\n",
    "the mean of $ L_t $ continue to be one despite most of the probability mass piling up\n",
    "near $ 0 $.\n",
    "\n",
    "To illustrate this peculiar property, we simulate many paths and\n",
    "calculate the unconditional mean of $ L\\left(w^t\\right) $ by\n",
    "averaging across these many paths at each $ t $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "l_arr_g = simulate(G_a, G_b, N=50000)\n",
    "l_seq_g = np.cumprod(l_arr_g, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following Python code approximates  unconditional means\n",
    "$ E_{0}\\left[L\\left(w^{t}\\right)\\right] $ by averaging across sample\n",
    "paths.\n",
    "\n",
    "Please notice that while  sample averages  hover around their population means of $ 1 $, there is quite a bit\n",
    "of variability, a consequence of the *fat tail* of the distribution of  $ L\\left(w^{t}\\right) $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "N, T = l_arr_g.shape\n",
    "plt.plot(range(T), np.mean(l_arr_g, axis=0))\n",
    "plt.hlines(1, 0, T, linestyle='--')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nature Permanently Draws from Density f\n",
    "\n",
    "Now suppose that before time $ 0 $ nature permanently decided to draw repeatedly from density $ f $.\n",
    "\n",
    "While the mean of the likelihood ratio $ \\ell \\left(w_{t}\\right) $ under density\n",
    "$ g $ is $ 1 $, its mean under the density $ f $ exceeds one.\n",
    "\n",
    "To see this, we compute\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "E_{0}\\left[\\ell \\left(w_{t}\\right)\\bigm|q=f\\right]  &=\\int\\frac{f\\left(w_{t}\\right)}{g\\left(w_{t}\\right)}f\\left(w_{t}\\right)dw_{t} \\\\\n",
    "    &=\\int\\frac{f\\left(w_{t}\\right)}{g\\left(w_{t}\\right)}\\frac{f\\left(w_{t}\\right)}{g\\left(w_{t}\\right)}g\\left(w_{t}\\right)dw_{t} \\\\\n",
    "    &=\\int \\ell \\left(w_{t}\\right)^{2}g\\left(w_{t}\\right)dw_{t} \\\\\n",
    "    &=E_{0}\\left[\\ell \\left(w_{t}\\right)^{2}\\mid q=g\\right] \\\\\n",
    "    &=E_{0}\\left[\\ell \\left(w_{t}\\right)\\mid q=g\\right]^{2}+Var\\left(\\ell \\left(w_{t}\\right)\\mid q=g\\right) \\\\\n",
    "    &>E_{0}\\left[\\ell \\left(w_{t}\\right)\\mid q=g\\right]^{2} = 1 \\\\\n",
    "       \\end{aligned}\n",
    "$$\n",
    "\n",
    "This in turn implies that the unconditional mean of the likelihood ratio process $ L(w^t) $\n",
    "diverges toward $ + \\infty $.\n",
    "\n",
    "Simulations below confirm this conclusion.\n",
    "\n",
    "Please note the scale of the $ y $ axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "l_arr_f = simulate(F_a, F_b, N=50000)\n",
    "l_seq_f = np.cumprod(l_arr_f, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "N, T = l_arr_f.shape\n",
    "plt.plot(range(T), np.mean(l_seq_f, axis=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also plot the probability that $ L\\left(w^t\\right) $ falls into\n",
    "the interval $ [10000, \\infty) $ as a function of time and watch how\n",
    "fast probability mass diverges  to $ +\\infty $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plt.plot(range(T), np.sum(l_seq_f > 10000, axis=0) / N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Likelihood Ratio Test\n",
    "\n",
    "We now describe how to employ the machinery\n",
    "of Neyman and Pearson [[NP33]](https://python-programming.quantecon.org/zreferences.html#neyman-pearson) to test the hypothesis that  history $ w^t $ is generated by repeated\n",
    "IID draws from density $ g $.\n",
    "\n",
    "Denote $ q $ as the data generating process, so that\n",
    "$ q=f \\text{ or } g $.\n",
    "\n",
    "Upon observing a sample $ \\{W_i\\}_{i=1}^t $, we want to decide\n",
    "which one is the data generating process by performing  a (frequentist)\n",
    "hypothesis test.\n",
    "\n",
    "We specify\n",
    "\n",
    "- Null hypothesis $ H_0 $: $ q=f $,  \n",
    "- Alternative hypothesis $ H_1 $: $ q=g $.  \n",
    "\n",
    "\n",
    "Neyman and Pearson proved that the best way to test this hypothesis is to use a **likelihood ratio test** that takes the\n",
    "form:\n",
    "\n",
    "- reject $ H_0 $ if $ L(W^t) < c $,  \n",
    "- accept $ H_0 $ otherwise.  \n",
    "\n",
    "\n",
    "where $ c $ is a given  discrimination threshold, to be chosen in a way we’ll soon describe.\n",
    "\n",
    "This test is *best* in the sense that it is a **uniformly most powerful** test.\n",
    "\n",
    "To understand what this means, we have to define probabilities of two important events that\n",
    "allow us to characterize a test associated with given\n",
    "threshold $ c $.\n",
    "\n",
    "The two probabities are:\n",
    "\n",
    "- \n",
    "  <dl style='margin: 20px 0;'>\n",
    "  <dt>Probability of detection (= power = 1 minus probability</dt>\n",
    "  <dd>\n",
    "  of Type II error):  \n",
    "  </dd>\n",
    "  \n",
    "  </dl>\n",
    "  \n",
    "\n",
    "\n",
    "$$\n",
    "1-\\beta \\equiv \\Pr\\left\\{ L\\left(w^{t}\\right)<c\\mid q=g\\right\\}\n",
    "$$\n",
    "\n",
    "- \n",
    "  <dl style='margin: 20px 0;'>\n",
    "  <dt>Probability of false alarm (= significance level = probability of</dt>\n",
    "  <dd>\n",
    "  Type I error):  \n",
    "  </dd>\n",
    "  \n",
    "  </dl>\n",
    "  \n",
    "\n",
    "\n",
    "$$\n",
    "\\alpha \\equiv  \\Pr\\left\\{ L\\left(w^{t}\\right)<c\\mid q=f\\right\\}\n",
    "$$\n",
    "\n",
    "The [Neyman-Pearson\n",
    "Lemma](https://en.wikipedia.org/wiki/Neyman–Pearson_lemma)\n",
    "states that among all possible tests, a likelihood ratio test\n",
    "maximizes the probability of detection for a given probability of false\n",
    "alarm.\n",
    "\n",
    "Another way to say the same thing is that  among all possible tests, a likelihood ratio test\n",
    "maximizes **power** for a given **significance level**.\n",
    "\n",
    "To have made a confident inference, we want a small probability of\n",
    "false alarm and a large probability of detection.\n",
    "\n",
    "With sample size $ t $ fixed, we can change our two probabilities by\n",
    "adjusting $ c $.\n",
    "\n",
    "A troublesome “that’s life” fact is that these two probabilities  move in the same direction as we vary the critical value\n",
    "$ c $.\n",
    "\n",
    "Without specifying quantitative losses from making Type I and Type II errors, there is little that we can say\n",
    "about how we *should*  trade off probabilities of the two types of mistakes.\n",
    "\n",
    "We do know that increasing sample size $ t $ improves\n",
    "statistical inference.\n",
    "\n",
    "Below we plot some informative figures that illustrate this.\n",
    "\n",
    "We also present a classical frequentist method for choosing a sample\n",
    "size $ t $.\n",
    "\n",
    "Let’s start with a case in which we fix the threshold $ c $ at\n",
    "$ 1 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "c = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we plot empirical distributions of logarithms of the cumulative\n",
    "likelihood ratios simulated above, which are generated by either\n",
    "$ f $ or $ g $.\n",
    "\n",
    "Taking logarithms has no effect on calculating the probabilities because\n",
    "the log  is a monotonic transformation.\n",
    "\n",
    "As $ t $ increases, the probabilities of making Type I and Type II\n",
    "errors both decrease, which is good.\n",
    "\n",
    "This is because most of the probability mass of log$ (L(w^t)) $\n",
    "moves toward $ -\\infty $ when $ g $ is the data generating\n",
    "process, ; while log$ (L(w^t)) $ goes to\n",
    "$ \\infty $ when data are generated by $ f $.\n",
    "\n",
    "This diverse behavior is what makes it possible to distinguish\n",
    "$ q=f $ from $ q=g $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(2, 2, figsize=(12, 8))\n",
    "fig.suptitle('distribution of $log(L(w^t))$ under f or g', fontsize=15)\n",
    "\n",
    "for i, t in enumerate([1, 7, 14, 21]):\n",
    "    nr = i // 2\n",
    "    nc = i % 2\n",
    "\n",
    "    axs[nr, nc].axvline(np.log(c), color=\"k\", ls=\"--\")\n",
    "\n",
    "    hist_f, x_f = np.histogram(np.log(l_seq_f[:, t]), 200, density=True)\n",
    "    hist_g, x_g = np.histogram(np.log(l_seq_g[:, t]), 200, density=True)\n",
    "\n",
    "    axs[nr, nc].plot(x_f[1:], hist_f, label=\"dist under f\")\n",
    "    axs[nr, nc].plot(x_g[1:], hist_g, label=\"dist under g\")\n",
    "\n",
    "    for i, (x, hist, label) in enumerate(zip([x_f, x_g], [hist_f, hist_g], [\"Type I error\", \"Type II error\"])):\n",
    "        ind = x[1:] <= np.log(c) if i == 0 else x[1:] > np.log(c)\n",
    "        axs[nr, nc].fill_between(x[1:][ind], hist[ind], alpha=0.5, label=label)\n",
    "\n",
    "    axs[nr, nc].legend()\n",
    "    axs[nr, nc].set_title(f\"t={t}\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The graph below shows more clearly that, when we hold the threshold\n",
    "$ c $ fixed, the probability of detection monotonically increases with increases in\n",
    "$ t $ and that the probability of a false alarm monotonically decreases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "PD = np.empty(T)\n",
    "PFA = np.empty(T)\n",
    "\n",
    "for t in range(T):\n",
    "    PD[t] = np.sum(l_seq_g[:, t] < c) / N\n",
    "    PFA[t] = np.sum(l_seq_f[:, t] < c) / N\n",
    "\n",
    "plt.plot(range(T), PD, label=\"Probability of detection\")\n",
    "plt.plot(range(T), PFA, label=\"Probability of false alarm\")\n",
    "plt.xlabel(\"t\")\n",
    "plt.title(\"$c=1$\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a given sample size $ t $,  the threshold $ c $ uniquely pins down  probabilities\n",
    "of both types of error.\n",
    "\n",
    "If for a fixed $ t $ we now free up and move $ c $, we will sweep out the probability\n",
    "of detection as a function of the probability of false alarm.\n",
    "\n",
    "This produces what is called a [receiver operating characteristic\n",
    "curve](https://en.wikipedia.org/wiki/Receiver_operating_characteristic) for a given discrimination threshold $ c $.\n",
    "\n",
    "Below, we plot receiver operating characteristic curves for a given discrimination threshold $ c $ but different\n",
    "sample sizes $ t $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "PFA = np.arange(0, 100, 1)\n",
    "\n",
    "for t in range(1, 15, 4):\n",
    "    percentile = np.percentile(l_seq_f[:, t], PFA)\n",
    "    PD = [np.sum(l_seq_g[:, t] < p) / N for p in percentile]\n",
    "\n",
    "    plt.plot(PFA / 100, PD, label=f\"t={t}\")\n",
    "\n",
    "plt.scatter(0, 1, label=\"perfect detection\")\n",
    "plt.plot([0, 1], [0, 1], color='k', ls='--', label=\"random detection\")\n",
    "\n",
    "plt.arrow(0.5, 0.5, -0.15, 0.15, head_width=0.03)\n",
    "plt.text(0.35, 0.7, \"better\")\n",
    "plt.xlabel(\"Probability of false alarm\")\n",
    "plt.ylabel(\"Probability of detection\")\n",
    "plt.legend()\n",
    "plt.title(\"Receiver Operating Characteristic Curve\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that as $ t $ increases, we are assured a larger probability\n",
    "of detection and a smaller probability of false alarm associated with\n",
    "a given discrimination threshold $ c $.\n",
    "\n",
    "As $ t \\rightarrow + \\infty $, we approach the the perfect detection\n",
    "curve that is indicated by a right angle hinging on the green dot.\n",
    "\n",
    "For a given sample size $ t $, a value discrimination threshold $ c $ determines a point on the receiver operating\n",
    "characteristic curve.\n",
    "\n",
    "It is up to the test designer to trade off probabilities of\n",
    "making the two types of errors.\n",
    "\n",
    "But we know how to choose the smallest sample size to achieve given targets for\n",
    "the probabilities.\n",
    "\n",
    "Typically, frequentists aim for a high probability of detection that\n",
    "respects an upper bound on the probability of false alarm.\n",
    "\n",
    "Below we show an example in which we fix the probability of false alarm at\n",
    "$ 0.05 $.\n",
    "\n",
    "The required sample size for making a decision is then determined by a\n",
    "target probability of detection, for example, $ 0.9 $, as depicted in the following graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "PFA = 0.05\n",
    "PD = np.empty(T)\n",
    "\n",
    "for t in range(T):\n",
    "\n",
    "    c = np.percentile(l_seq_f[:, t], PFA * 100)\n",
    "    PD[t] = np.sum(l_seq_g[:, t] < c) / N\n",
    "\n",
    "plt.plot(range(T), PD)\n",
    "plt.axhline(0.9, color=\"k\", ls=\"--\")\n",
    "\n",
    "plt.xlabel(\"t\")\n",
    "plt.ylabel(\"Probability of detection\")\n",
    "plt.title(f\"Probability of false alarm={PFA}\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The United States Navy evidently used a procedure like this to select a sample size $ t $ for doing quality\n",
    "control tests during World War II.\n",
    "\n",
    "A Navy Captain who had been ordered to perform tests of this kind had second thoughts about it that he\n",
    "presented to Milton Friedman, as we describe in  [this lecture](https://python-programming.quantecon.org/wald_friedman.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sequels\n",
    "\n",
    "Likelihood processes play an important role in Bayesian learning, as described in [this lecture](https://python-programming.quantecon.org/likelihood_bayes.html)\n",
    "and as applied in [this lecture](https://python-programming.quantecon.org/odu.html).\n",
    "\n",
    "Likelihood ratio processes appear again in [this lecture](https://python-advanced.quantecon.org/additive_functionals.html), which contains another illustration\n",
    "of the **peculiar property** of likelihood ratio processes described above."
   ]
  }
 ],
 "metadata": {
  "date": 1591697896.577528,
  "filename": "likelihood_ratio_process.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Likelihood Ratio Processes"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}