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 Process and Bayes’ Law](#Likelihood-Ratio-Process-and-Bayes’-Law)  \n",
    "  - [Likelihood Ratio Test](#Likelihood-Ratio-Test)  "
   ]
  },
  {
   "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 by frequentist and Bayesian statisticians.\n",
    "\n",
    "We’ll use the simple statistical setting also used in [this lecture](https://python.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 for a Bayesian  a likelihood ratio process and a prior probability determine a posterior probability  \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  "
   ]
  },
  {
   "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 $ before the beginning of time.\n",
    "\n",
    "A **likelihood ratio process** is a useful tool for this task.\n",
    "\n",
    "To begin, we define the 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 $ w_t $.\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\\{ \\ell \\left(w_{t}\\right)\\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.quantecon.org/zreferences.html#neyman-pearson).\n",
    "\n",
    "To help us appreciate how things work, the following Python code evaluates $ f $ and $ g $ as 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 virtually all probability mass to a\n",
    "very small interval near $ 0 $, a peculiar fact about a likelihood\n",
    "ratio process is that 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) $ is a multiplicative process 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 this possibly be true if probability mass of the likelihood\n",
    "ratio process is piling up near $ 1 $ 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 thin-tailed as\n",
    "sufficient small mass shifts to very large values of $ L_t $ to make\n",
    "the mean of $ L_t $ continue to be zero despite the 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 the 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",
    "A useful property is that 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} \\\\\n",
    "    &=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 Process and Bayes’ Law\n",
    "\n",
    "Let $ \\pi_t $ be a Bayesian posterior defined as\n",
    "\n",
    "$$\n",
    "\\pi_t = {\\rm Prob}(q=f|w^t)\n",
    "$$\n",
    "\n",
    "The likelihood ratio process is a principal actor in the formula that governs the evolution\n",
    "of the posterior probability $ \\pi_t $, an instance of **Bayes’ Law**.\n",
    "\n",
    "Bayes’ law implies that $ \\{\\pi_t\\} $ obeys the recursion\n",
    "\n",
    "\n",
    "<a id='equation-eq-recur1'></a>\n",
    "$$\n",
    "\\pi_t=\\frac{\\pi_{t-1} l_t(w_t)}{\\pi_{t-1} l_t(w_t)+1-\\pi_{t-1}} \\tag{1}\n",
    "$$\n",
    "\n",
    "with $ \\pi_{0} $ be a Bayesian prior probability that $ q = f $,\n",
    "i.e., a belief about $ q $ based on having seen no data.\n",
    "\n",
    "Below we define a Python function that updates belief $ \\pi $ using\n",
    "likelihood ratio $ \\ell $ according to  recursion [(1)](#equation-eq-recur1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@njit\n",
    "def update(π, l):\n",
    "    \"Update π using likelihood l\"\n",
    "\n",
    "    # Update belief\n",
    "    π = π * l / (π * l + 1 - π)\n",
    "\n",
    "    return π"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Formula [(1)](#equation-eq-recur1) can be generalized in a useful way.\n",
    "\n",
    "We do this by iterating on recursion [(1)](#equation-eq-recur1) in order to derive an\n",
    "expression for  the time $ t $ posterior $ \\pi_{t+1} $ as a function\n",
    "of the time $ 0 $ prior $ \\pi_0 $ and the likelihood ratio process\n",
    "$ L(w^{t+1}) $ at time $ t $.\n",
    "\n",
    "To begin, notice that the updating rule\n",
    "\n",
    "$$\n",
    "\\pi_{t+1}\n",
    "=\\frac{\\pi_{t}\\ell \\left(w_{t+1}\\right)}\n",
    "{\\pi_{t}\\ell \\left(w_{t+1}\\right)+\\left(1-\\pi_{t}\\right)}\n",
    "$$\n",
    "\n",
    "implies\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\frac{1}{\\pi_{t+1}}\n",
    "    &=\\frac{\\pi_{t}\\ell \\left(w_{t+1}\\right)\n",
    "        +\\left(1-\\pi_{t}\\right)}{\\pi_{t}\\ell \\left(w_{t+1}\\right)} \\\\\n",
    "    &=1-\\frac{1}{\\ell \\left(w_{t+1}\\right)}\n",
    "        +\\frac{1}{\\ell \\left(w_{t+1}\\right)}\\frac{1}{\\pi_{t}}.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\Rightarrow\n",
    "\\frac{1}{\\pi_{t+1}}-1\n",
    "=\\frac{1}{\\ell \\left(w_{t+1}\\right)}\\left(\\frac{1}{\\pi_{t}}-1\\right).\n",
    "$$\n",
    "\n",
    "Therefore\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\frac{1}{\\pi_{t+1}}-1\n",
    "    =\\frac{1}{\\prod_{i=1}^{t+1}\\ell \\left(w_{i}\\right)}\n",
    "        \\left(\\frac{1}{\\pi_{0}}-1\\right)\n",
    "    =\\frac{1}{L\\left(w^{t+1}\\right)}\\left(\\frac{1}{\\pi_{0}}-1\\right).\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Since $ \\pi_{0}\\in\\left(0,1\\right) $ and\n",
    "$ L\\left(w^{t+1}\\right)>0 $, we can verify that\n",
    "$ \\pi_{t+1}\\in\\left(0,1\\right) $.\n",
    "\n",
    "After rearranging the preceding equation, we can express $ \\pi_{t+1} $ as a\n",
    "function of  $ L\\left(w^{t+1}\\right) $, the  likelihood ratio process at $ t+1 $,\n",
    "and the initial prior $ \\pi_{0} $\n",
    "\n",
    "\n",
    "<a id='equation-eq-bayeslaw103'></a>\n",
    "$$\n",
    "\\pi_{t+1}=\\frac{\\pi_{0}L\\left(w^{t+1}\\right)}{\\pi_{0}L\\left(w^{t+1}\\right)+1-\\pi_{0}} . \\tag{2}\n",
    "$$\n",
    "\n",
    "Formula [(2)](#equation-eq-bayeslaw103) generalizes generalizes formula [(1)](#equation-eq-recur1).\n",
    "\n",
    "Formula [(2)](#equation-eq-bayeslaw103)  can be regarded as a one step  revision of prior probability $ \\pi_0 $ after seeing\n",
    "the batch of data $ \\left\\{ w_{i}\\right\\} _{i=1}^{t+1} $.\n",
    "\n",
    "Formula [(2)](#equation-eq-bayeslaw103) shows the key role that the likelihood ratio process  $ L\\left(w^{t+1}\\right) $ plays in determining\n",
    "the posterior probability $ \\pi_{t+1} $.\n",
    "\n",
    "Formula [(2)](#equation-eq-bayeslaw103) is the foundation for the insight that, because of the way the likelihood ratio process behaves\n",
    "as $ t \\rightarrow + \\infty $, the likelihood ratio process dominates the initial prior $ \\pi_0 $ in determining the\n",
    "limiting behavior of $ \\pi_t $.\n",
    "\n",
    "To illustrate this insight, below we will plot  graphs showing **one** simulated\n",
    "path of the  likelihood ratio process $ L_t $ along with two paths of\n",
    "$ \\pi_t $ that are associated with the *same* realization of the likelihood ratio process but *different* initial prior probabilities\n",
    "probabilities $ \\pi_{0} $.\n",
    "\n",
    "First, we specify the two values of $ \\pi_0 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "π1, π2 = 0.2, 0.8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we generate paths of the likelihood ratio process $ L_t $ and the posteior $ \\pi_t $ for a history drawn as IID\n",
    "draws from density $ f $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "T = l_arr_f.shape[1]\n",
    "π_seq_f = np.empty((2, T+1))\n",
    "π_seq_f[:, 0] = π1, π2\n",
    "\n",
    "for t in range(T):\n",
    "    for i in range(2):\n",
    "        π_seq_f[i, t+1] = update(π_seq_f[i, t], l_arr_f[0, t])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots()\n",
    "\n",
    "for i in range(2):\n",
    "    ax1.plot(range(T+1), π_seq_f[i, :], label=f\"$\\pi_0$={π_seq_f[i, 0]}\")\n",
    "\n",
    "ax1.set_ylabel(\"$\\pi_t$\")\n",
    "ax1.set_xlabel(\"t\")\n",
    "ax1.legend()\n",
    "ax1.set_title(\"when f governs data\")\n",
    "\n",
    "ax2 = ax1.twinx()\n",
    "ax2.plot(range(1, T+1), np.log(l_seq_f[0, :]), '--', color='b')\n",
    "ax2.set_ylabel(\"$log(L(w^{t}))$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dotted line in the graph above records the logarithm of the  likelihood ratio process $ \\log L(w^t) $.\n",
    "\n",
    "Please note that there are two different scales on the $ y $ axis.\n",
    "\n",
    "Now let’s study what happens when the history consists of IID draws from density $ g $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "T = l_arr_g.shape[1]\n",
    "π_seq_g = np.empty((2, T+1))\n",
    "π_seq_g[:, 0] = π1, π2\n",
    "\n",
    "for t in range(T):\n",
    "    for i in range(2):\n",
    "        π_seq_g[i, t+1] = update(π_seq_g[i, t], l_arr_g[0, t])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots()\n",
    "\n",
    "for i in range(2):\n",
    "    ax1.plot(range(T+1), π_seq_g[i, :], label=f\"$\\pi_0$={π_seq_g[i, 0]}\")\n",
    "\n",
    "ax1.set_ylabel(\"$\\pi_t$\")\n",
    "ax1.set_xlabel(\"t\")\n",
    "ax1.legend()\n",
    "ax1.set_title(\"when g governs data\")\n",
    "\n",
    "ax2 = ax1.twinx()\n",
    "ax2.plot(range(1, T+1), np.log(l_seq_g[0, :]), '--', color='b')\n",
    "ax2.set_ylabel(\"$log(L(w^{t}))$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we offer Python code that verifies this in a setting in which\n",
    "nature chose permanently to draw from density $ f $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "π_seq = np.empty((2, T+1))\n",
    "π_seq[:, 0] = π1, π2\n",
    "\n",
    "for i in range(2):\n",
    "    πL = π_seq[i, 0] * l_seq_f[0, :]\n",
    "    π_seq[i, 1:] = πL / (πL + 1 - π_seq[i, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "np.abs(π_seq - π_seq_f).max() < 1e-10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Likelihood Ratio Test\n",
    "\n",
    "Having seen how the likelihood ratio process is a key ingredient of the formula [(2)](#equation-eq-bayeslaw103) for\n",
    "a Bayesian’s posteior probabilty that nature has drawn history $ w^t $ as repeated draws from density\n",
    "$ g $, we now turn to how a frequentist statistician would employ the hypothesis testing theory\n",
    "of Neyman and Pearson [[NP33]](https://python.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 figure out\n",
    "which one is the data generating process by executing 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 hospital is to use a **likelihood ratio test**:\n",
    "\n",
    "- reject $ H_0 $ if $ L(W^t) < c $,  \n",
    "- accept $ H_0 $ otherwise.  \n",
    "\n",
    "\n",
    "where the discrimination threshold $ c $ is arbitrarily given.\n",
    "\n",
    "This test is *best* in the sense that it is the uniformly most powerful test.\n",
    "\n",
    "To understand what this means, we have to look at probabilities of two important events.\n",
    "\n",
    "These two probabilities allow us to characterize a test associated with particular\n",
    "threshold $ c $.\n",
    "\n",
    "The two probabities are:\n",
    "\n",
    "1. Probability of detection (= *statistical power* = 1 minus probability\n",
    "  of Type II error):  \n",
    "\n",
    "\n",
    "$$\n",
    "1-\\beta \\equiv \\Pr\\left\\{ L\\left(w^{t}\\right)<c\\mid q=g\\right\\}\n",
    "$$\n",
    "\n",
    "1. Probability of false alarm (= *significance level* = probability of\n",
    "  Type I error):  \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, the likelihood ratio test\n",
    "maximizes the probability of detection for a given probability of false\n",
    "alarm.\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",
    "moving $ 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 the penalties on 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",
    "What we do know that increasing sample size $ t $ improves\n",
    "statistical inference.\n",
    "\n",
    "Below we plot some informative figures that display showing this.\n",
    "\n",
    "We also present a classical frequentist method for choosing a sample\n",
    "size $ t $ that provides confident inferences.\n",
    "\n",
    "Let’s start with a case where 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 in\n",
    "$ t $ and that the probability of a false alarm moves in the opposite\n",
    "direction."
   ]
  },
  {
   "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": [
    "Notice that the threshold $ c $ uniquely pins down  probabilities\n",
    "of both types of error.\n",
    "\n",
    "If we now free up $ c $ and move it, 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).\n",
    "\n",
    "Below, we plot receiver operating characteristic curves for 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",
    "It is up to the test user to decide how to trade off probabilities of\n",
    "making two types of errors.\n",
    "\n",
    "But we know how to choose the smallest sample size given any 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 a case 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 $."
   ]
  },
  {
   "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()"
   ]
  }
 ],
 "metadata": {
  "date": 1584334742.4463534,
  "filename": "likelihood_ratio_process.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Likelihood Ratio Processes"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}