{ "cells": [ { "cell_type": "markdown", "id": "06215005", "metadata": {}, "source": [ "$$\n", "\\newcommand{\\argmax}{arg\\,max}\n", "\\newcommand{\\argmin}{arg\\,min}\n", "$$" ] }, { "cell_type": "markdown", "id": "f6ea9565", "metadata": {}, "source": [ "\n", "\n", "
\n", " \n", " \"QuantEcon\"\n", " \n", "
" ] }, { "cell_type": "markdown", "id": "518eb90c", "metadata": {}, "source": [ "# Bayesian versus Frequentist Decision Rules" ] }, { "cell_type": "markdown", "id": "b0db23c3", "metadata": {}, "source": [ "## Contents\n", "\n", "- [Bayesian versus Frequentist Decision Rules](#Bayesian-versus-Frequentist-Decision-Rules) \n", " - [Overview](#Overview) \n", " - [Setup](#Setup) \n", " - [Frequentist Decision Rule](#Frequentist-Decision-Rule) \n", " - [Bayesian Decision Rule](#Bayesian-Decision-Rule) \n", " - [Was the Navy Captain’s Hunch Correct?](#Was-the-Navy-Captain’s-Hunch-Correct?) \n", " - [More Details](#More-Details) \n", " - [Distribution of Bayesian Decision Rule’s Time to Decide](#Distribution-of-Bayesian-Decision-Rule’s-Time-to-Decide) \n", " - [Probability of Making Correct Decision](#Probability-of-Making-Correct-Decision) \n", " - [Distribution of Likelihood Ratios at Neyman-Pearson’s $ t $](#Distribution-of-Likelihood-Ratios-at-Neyman-Pearson’s-$-t-$) " ] }, { "cell_type": "code", "execution_count": null, "id": "02141214", "metadata": { "hide-output": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from numba import jit, prange, float64, int64\n", "from numba.experimental import jitclass\n", "from math import gamma\n", "from scipy.optimize import minimize" ] }, { "cell_type": "markdown", "id": "09a708d4", "metadata": {}, "source": [ "## Overview\n", "\n", "This lecture follows up on ideas presented in the following lectures:\n", "\n", "- [A Problem that Stumped Milton Friedman](https://python.quantecon.org/wald_friedman.html) \n", "- [A Bayesian Formulation of Friedman and Wald’s Problem](https://python.quantecon.org/wald_friedman_2.html) \n", "- [Exchangeability and Bayesian Updating](https://python.quantecon.org/exchangeable.html) \n", "- [Likelihood Ratio Processes](https://python.quantecon.org/likelihood_ratio_process.html) \n", "\n", "\n", "[A Problem that Stumped Milton Friedman](https://python.quantecon.org/wald_friedman.html) described a problem\n", "that a Navy Captain presented to Milton Friedman during World War II.\n", "\n", "The Navy had told the Captain to use a decision rule for quality control.\n", "\n", "In particular, the Navy had ordered the Captain to use an instance of a **frequentist decision rule**.\n", "\n", "The Captain doubted that that rule was a good one.\n", "\n", "Milton Friedman recognized the Captain’s conjecture as posing a challenging statistical problem that he and other members of the US Government’s Statistical Research Group at Columbia University proceeded to try to solve.\n", "\n", "A member of the group, the great mathematician and economist Abraham Wald, soon solved the problem.\n", "\n", "A good way to formulate the problem is to use some ideas from Bayesian statistics that we describe in\n", "this lecture [Exchangeability and Bayesian Updating](https://python.quantecon.org/exchangeable.html) and in this lecture\n", "[Likelihood Ratio Processes](https://python.quantecon.org/likelihood_ratio_process.html), which describes the link between Bayesian\n", "updating and likelihood ratio processes.\n", "\n", "The present lecture uses Python to generate simulations that evaluate expected losses under the Neyman-Pearson **frequentist** procedure that the Navy captain questioned and the **Bayesian** decision rule described in [A Bayesian Formulation of Friedman and Wald’s Problem](https://python.quantecon.org/wald_friedman_2.html).\n", "\n", "The simulations confirm the Navy Captain’s hunch that there is a better rule than the Neyman-Pearson likelihood ratio test that the Navy had told him to use." ] }, { "cell_type": "markdown", "id": "01eb1711", "metadata": {}, "source": [ "## Setup\n", "\n", "To formalize the problem that had confronted the Navy Captain, we consider a setting with the following parts.\n", "\n", "- Each period a decision maker draws a non-negative random variable\n", " $ Z $. He knows that two probability distributions are possible,\n", " $ f_{0} $ and $ f_{1} $, and that which ever distribution it\n", " is remains fixed over time. The decision maker believes that before\n", " the beginning of time, nature once and for all had selected either\n", " $ f_{0} $ or $ f_1 $ and that the probability that it\n", " selected $ f_0 $ is probability $ \\pi^{*} $. \n", "- The decision maker observes a sample\n", " $ \\left\\{ z_{i}\\right\\} _{i=0}^{t} $ from the distribution\n", " chosen by nature. \n", "\n", "\n", "The decision maker wants to decide which distribution actually governs\n", "$ Z $.\n", "\n", "He is worried about two types of errors and the losses that they will\n", "impose on him.\n", "\n", "- a loss $ \\bar L_{1} $ from a **type I error** that occurs if he decides that\n", " $ f=f_{1} $ when actually $ f=f_{0} $ \n", "- a loss $ \\bar L_{0} $ from a **type II error** that occurs if he decides that\n", " $ f=f_{0} $ when actually $ f=f_{1} $ \n", "\n", "\n", "The decision maker pays a cost $ c $ for drawing\n", "another $ z $.\n", "\n", "We mainly borrow parameters from the quantecon lecture\n", "[A Bayesian Formulation of Friedman and Wald’s Problem](https://python.quantecon.org/wald_friedman_2.html) except that we increase both $ \\bar L_{0} $\n", "and $ \\bar L_{1} $ from $ 25 $ to $ 100 $ to encourage the\n", "Bayesian decision rule to take more draws before deciding.\n", "\n", "We set the cost $ c $ of taking one more draw at $ 1.25 $.\n", "\n", "We set the probability distributions $ f_{0} $ and $ f_{1} $ to\n", "be beta distributions with $ a_{0}=b_{0}=1 $, $ a_{1}=3 $, and\n", "$ b_{1}=1.2 $, respectively.\n", "\n", "Below is some Python code that sets up these objects." ] }, { "cell_type": "code", "execution_count": null, "id": "9ea3cf82", "metadata": { "hide-output": false }, "outputs": [], "source": [ "@jit\n", "def p(x, a, b):\n", " \"Beta distribution.\"\n", "\n", " r = gamma(a + b) / (gamma(a) * gamma(b))\n", "\n", " return r * x**(a-1) * (1 - x)**(b-1)" ] }, { "cell_type": "markdown", "id": "b8e4d35b", "metadata": {}, "source": [ "We start with defining a `jitclass` that stores parameters and\n", "functions we need to solve problems for both the Bayesian and\n", "frequentist Navy Captains." ] }, { "cell_type": "code", "execution_count": null, "id": "cc38633f", "metadata": { "hide-output": false }, "outputs": [], "source": [ "wf_data = [\n", " ('c', float64), # unemployment compensation\n", " ('a0', float64), # parameters of beta distribution\n", " ('b0', float64),\n", " ('a1', float64),\n", " ('b1', float64),\n", " ('L0', float64), # cost of selecting f0 when f1 is true\n", " ('L1', float64), # cost of selecting f1 when f0 is true\n", " ('π_grid', float64[:]), # grid of beliefs π\n", " ('π_grid_size', int64),\n", " ('mc_size', int64), # size of Monto Carlo simulation\n", " ('z0', float64[:]), # sequence of random values\n", " ('z1', float64[:]) # sequence of random values\n", "]" ] }, { "cell_type": "code", "execution_count": null, "id": "8e1553d6", "metadata": { "hide-output": false }, "outputs": [], "source": [ "@jitclass(wf_data)\n", "class WaldFriedman:\n", "\n", " def __init__(self,\n", " c=1.25,\n", " a0=1,\n", " b0=1,\n", " a1=3,\n", " b1=1.2,\n", " L0=100,\n", " L1=100,\n", " π_grid_size=200,\n", " mc_size=1000):\n", "\n", " self.c, self.π_grid_size = c, π_grid_size\n", " self.a0, self.b0, self.a1, self.b1 = a0, b0, a1, b1\n", " self.L0, self.L1 = L0, L1\n", " self.π_grid = np.linspace(0, 1, π_grid_size)\n", " self.mc_size = mc_size\n", "\n", " self.z0 = np.random.beta(a0, b0, mc_size)\n", " self.z1 = np.random.beta(a1, b1, mc_size)\n", "\n", " def f0(self, x):\n", "\n", " return p(x, self.a0, self.b0)\n", "\n", " def f1(self, x):\n", "\n", " return p(x, self.a1, self.b1)\n", "\n", " def κ(self, z, π):\n", " \"\"\"\n", " Updates π using Bayes' rule and the current observation z\n", " \"\"\"\n", "\n", " a0, b0, a1, b1 = self.a0, self.b0, self.a1, self.b1\n", "\n", " π_f0, π_f1 = π * p(z, a0, b0), (1 - π) * p(z, a1, b1)\n", " π_new = π_f0 / (π_f0 + π_f1)\n", "\n", " return π_new" ] }, { "cell_type": "code", "execution_count": null, "id": "adc3457d", "metadata": { "hide-output": false }, "outputs": [], "source": [ "wf = WaldFriedman()\n", "\n", "grid = np.linspace(0, 1, 50)\n", "\n", "plt.figure()\n", "\n", "plt.title(\"Two Distributions\")\n", "plt.plot(grid, wf.f0(grid), lw=2, label=\"$f_0$\")\n", "plt.plot(grid, wf.f1(grid), lw=2, label=\"$f_1$\")\n", "\n", "plt.legend()\n", "plt.xlabel(\"$z$ values\")\n", "plt.ylabel(\"density of $z_k$\")\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "5b87d775", "metadata": {}, "source": [ "Above, we plot the two possible probability densities $ f_0 $ and\n", "$ f_1 $" ] }, { "cell_type": "markdown", "id": "8a6f183a", "metadata": {}, "source": [ "## Frequentist Decision Rule\n", "\n", "The Navy told the Captain to use a Neyman-Pearson likelihood ratio decision rule.\n", "\n", "That decision rule is characterized by\n", "\n", "- a sample size $ t $, and \n", "- a cutoff value $ d $ of a likelihood ratio \n", "\n", "\n", "Let\n", "$ L\\left(z^{t}\\right)=\\prod_{i=0}^{t}\\frac{f_{0}\\left(z_{i}\\right)}{f_{1}\\left(z_{i}\\right)} $\n", "be the likelihood ratio associated with observing the sequence\n", "$ \\left\\{ z_{i}\\right\\} _{i=0}^{t} $.\n", "\n", "The decision rule associated with a sample size $ t $ is:\n", "\n", "- decide that $ f_0 $ is the distribution if the likelihood ratio\n", " is greater than $ d $ \n", "- decide that $ f_1 $ is the distribution if the likelihood ratio is less than $ d $ \n", "\n", "\n", "For our purposes here, we want to compute an expected loss from using this rule, where we borrow\\\\\n", "\n", "\n", "loss parameters $ \\bar L_1 $ and $ \\bar L_2 $ from [A Bayesian Formulation of Friedman and Wald’s Problem](https://python.quantecon.org/wald_friedman_2.html).\n", "\n", "Let null and alternative\n", "hypotheses be\n", "\n", "- null: $ H_{0} $: $ f=f_{0} $, \n", "- alternative $ H_{1} $: $ f=f_{1} $. \n", "\n", "\n", "Given sample size $ t $ and cutoff $ d $, under the model\n", "described above, the mathematical expectation of total loss is\n", "\n", "\n", "\n", "$$\n", "\\begin{aligned}\n", "\\bar{V}_{fre}\\left(t,d\\right)=ct+\\pi^{*}PFA\\times \\bar L_{1}+\\left(1-\\pi^{*}\\right)\\left(1-PD\\right)\\times \\bar L_{0}\n", "\\end{aligned} \\tag{31.1}\n", "$$\n", "\n", "$$\n", "\\begin{aligned}\n", "\\textrm{where} \\quad PFA & =\\Pr\\left\\{ L\\left(z^{t}\\right) tol:\n", " h_new = Q(h, wf)\n", " error = np.max(np.abs(h - h_new))\n", " i += 1\n", " h = h_new\n", "\n", " if error > tol:\n", " print(\"Failed to converge!\")\n", "\n", " return h_new" ] }, { "cell_type": "code", "execution_count": null, "id": "0d8ea611", "metadata": { "hide-output": false }, "outputs": [], "source": [ "h_star = solve_model(wf)" ] }, { "cell_type": "code", "execution_count": null, "id": "427b8f8c", "metadata": { "hide-output": false }, "outputs": [], "source": [ "@jit\n", "def find_cutoff_rule(wf, h):\n", "\n", " \"\"\"\n", " This function takes a continuation value function and returns the\n", " corresponding cutoffs of where you transition between continuing and\n", " choosing a specific model\n", " \"\"\"\n", "\n", " π_grid = wf.π_grid\n", " L0, L1 = wf.L0, wf.L1\n", "\n", " # Evaluate cost at all points on grid for choosing a model\n", " payoff_f0 = (1 - π_grid) * L0\n", " payoff_f1 = π_grid * L1\n", "\n", " # The cutoff points can be found by differencing these costs with\n", " # The Bellman equation (J is always less than or equal to p_c_i)\n", " B = π_grid[np.searchsorted(\n", " payoff_f1 - np.minimum(h, payoff_f0),\n", " 1e-10)\n", " - 1]\n", " A = π_grid[np.searchsorted(\n", " np.minimum(h, payoff_f1) - payoff_f0,\n", " 1e-10)\n", " - 1]\n", "\n", " return (B, A)\n", "\n", "B, A = find_cutoff_rule(wf, h_star)\n", "cost_L0 = (1 - wf.π_grid) * wf.L0\n", "cost_L1 = wf.π_grid * wf.L1\n", "\n", "fig, ax = plt.subplots(figsize=(10, 6))\n", "\n", "ax.plot(wf.π_grid, h_star, label='continuation value')\n", "ax.plot(wf.π_grid, cost_L1, label='choose f1')\n", "ax.plot(wf.π_grid, cost_L0, label='choose f0')\n", "ax.plot(wf.π_grid,\n", " np.amin(np.column_stack([h_star, cost_L0, cost_L1]),axis=1),\n", " lw=15, alpha=0.1, color='b', label='minimum cost')\n", "\n", "ax.annotate(r\"$B$\", xy=(B + 0.01, 0.5), fontsize=14)\n", "ax.annotate(r\"$A$\", xy=(A + 0.01, 0.5), fontsize=14)\n", "\n", "plt.vlines(B, 0, B * wf.L0, linestyle=\"--\")\n", "plt.vlines(A, 0, (1 - A) * wf.L1, linestyle=\"--\")\n", "\n", "ax.set(xlim=(0, 1), ylim=(0, 0.5 * max(wf.L0, wf.L1)), ylabel=\"cost\",\n", " xlabel=r\"$\\pi$\", title=\"Value function\")\n", "\n", "plt.legend(borderpad=1.1)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "032757a6", "metadata": {}, "source": [ "The above figure portrays the value function plotted against the decision\n", "maker’s Bayesian posterior.\n", "\n", "It also shows the cutoff probabilities $ A $ and $ B $.\n", "\n", "The Bayesian decision rule is:\n", "\n", "- accept $ H_0 $ if $ \\pi \\geq A $ \n", "- accept $ H_1 $ if $ \\pi \\leq B $ \n", "- delay deciding and draw another $ z $ if\n", " $ B \\leq \\pi \\leq A $ \n", "\n", "\n", "We can calculate two “objective” loss functions under this situation\n", "conditioning on knowing for sure that nature has selected $ f_{0} $,\n", "in the first case, or $ f_{1} $, in the second case.\n", "\n", "1. under $ f_{0} $, \n", " $$\n", " V^{0}\\left(\\pi\\right)=\\begin{cases}\n", " 0 & \\text{if} A \\leq\\pi,\\\\\n", " c+EV^{0}\\left(\\pi^{\\prime}\\right) & \\text{if }B\\leq\\pi< A,\\\\\n", " \\bar L_{1} & \\text{if }\\pi= A] = wf.L0\n", "\n", " V_old = np.empty_like(V)\n", "\n", " while True:\n", " V_old[:] = V[:]\n", " V[(B <= wf.π_grid) & (wf.π_grid < A)] = 0\n", "\n", " for i in prange(len(wf.π_grid)):\n", " π = wf.π_grid[i]\n", "\n", " if π >= A or π < B:\n", " continue\n", "\n", " for j in prange(len(z_arr)):\n", " π_next = wf.κ(z_arr[j], π)\n", " V[i] += wf.c + np.interp(π_next, wf.π_grid, V_old)\n", "\n", " V[i] /= wf.mc_size\n", "\n", " if np.abs(V - V_old).max() < 1e-5:\n", " break\n", "\n", " return V" ] }, { "cell_type": "code", "execution_count": null, "id": "f4fd7564", "metadata": { "hide-output": false }, "outputs": [], "source": [ "V0 = V_q(wf, 0)\n", "V1 = V_q(wf, 1)\n", "\n", "plt.plot(wf.π_grid, V0, label='$V^0$')\n", "plt.plot(wf.π_grid, V1, label='$V^1$')\n", "plt.vlines(B, 0, wf.L0, linestyle='--')\n", "plt.text(B+0.01, wf.L0/2, 'B')\n", "plt.vlines(A, 0, wf.L0, linestyle='--')\n", "plt.text(A+0.01, wf.L0/2, 'A')\n", "plt.xlabel(r'$\\pi$')\n", "plt.title(r'Objective value function $V(\\pi)$')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "02dde8fe", "metadata": {}, "source": [ "Given an assumed value for\n", "$ \\pi^{*}=\\Pr\\left\\{ \\text{nature selects }f_{0}\\right\\} $, we can\n", "then compute $ \\bar{V}_{Bayes}\\left(\\pi_{0}\\right) $.\n", "\n", "We can then determine an initial Bayesian prior $ \\pi_{0}^{*} $ that\n", "minimizes this objective concept of expected loss.\n", "\n", "The figure below plots four cases corresponding to\n", "$ \\pi^{*}=0.25,0.3,0.5,0.7 $.\n", "\n", "We observe that in each case $ \\pi_{0}^{*} $ equals $ \\pi^{*} $." ] }, { "cell_type": "code", "execution_count": null, "id": "6c3a74fe", "metadata": { "hide-output": false }, "outputs": [], "source": [ "def compute_V_baye_bar(π_star, V0, V1, wf):\n", "\n", " V_baye = π_star * V0 + (1 - π_star) * V1\n", " π_idx = np.argmin(V_baye)\n", " π_optimal = wf.π_grid[π_idx]\n", " V_baye_bar = V_baye[π_idx]\n", "\n", " return V_baye, π_optimal, V_baye_bar" ] }, { "cell_type": "code", "execution_count": null, "id": "e0b10492", "metadata": { "hide-output": false }, "outputs": [], "source": [ "π_star_arr = [0.25, 0.3, 0.5, 0.7]\n", "\n", "fig, axs = plt.subplots(2, 2, figsize=(15, 10))\n", "\n", "for i, π_star in enumerate(π_star_arr):\n", " row_i = i // 2\n", " col_i = i % 2\n", "\n", " V_baye, π_optimal, V_baye_bar = compute_V_baye_bar(π_star, V0, V1, wf)\n", "\n", " axs[row_i, col_i].plot(wf.π_grid, V_baye)\n", " axs[row_i, col_i].hlines(V_baye_bar, 0, 1, linestyle='--')\n", " axs[row_i, col_i].vlines(π_optimal, V_baye_bar, V_baye.max(), linestyle='--')\n", " axs[row_i, col_i].text(π_optimal+0.05, (V_baye_bar + V_baye.max()) / 2,\n", " r'${\\pi_0^*}=$'+f'{π_optimal:0.2f}')\n", " axs[row_i, col_i].set_xlabel(r'$\\pi$')\n", " axs[row_i, col_i].set_ylabel(r'$\\overline{V}_{baye}(\\pi)$')\n", " axs[row_i, col_i].set_title(r'$\\pi^*=$' + f'{π_star}')\n", "\n", "fig.suptitle(r'$\\overline{V}_{baye}(\\pi)=\\pi^*V^0(\\pi) + (1-\\pi^*)V^1(\\pi)$', fontsize=16)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "d80831e2", "metadata": {}, "source": [ "This pattern of outcomes holds more generally.\n", "\n", "Thus, the following Python code generates the associated graph that\n", "verifies the equality of $ \\pi_{0}^{*} $ to $ \\pi^{*} $ holds\n", "for all $ \\pi^{*} $." ] }, { "cell_type": "code", "execution_count": null, "id": "430fb7ea", "metadata": { "hide-output": false }, "outputs": [], "source": [ "π_star_arr = np.linspace(0.1, 0.9, n_π)\n", "V_baye_bar_arr = np.empty_like(π_star_arr)\n", "π_optimal_arr = np.empty_like(π_star_arr)\n", "\n", "for i, π_star in enumerate(π_star_arr):\n", "\n", " V_baye, π_optimal, V_baye_bar = compute_V_baye_bar(π_star, V0, V1, wf)\n", "\n", " V_baye_bar_arr[i] = V_baye_bar\n", " π_optimal_arr[i] = π_optimal\n", "\n", "fig, axs = plt.subplots(1, 2, figsize=(14, 5))\n", "\n", "axs[0].plot(π_star_arr, V_baye_bar_arr)\n", "axs[0].set_xlabel(r'$\\pi^*$')\n", "axs[0].set_title(r'$\\overline{V}_{baye}$')\n", "\n", "axs[1].plot(π_star_arr, π_optimal_arr, label='optimal prior')\n", "axs[1].plot([π_star_arr.min(), π_star_arr.max()],\n", " [π_star_arr.min(), π_star_arr.max()],\n", " c='k', linestyle='--', label='45 degree line')\n", "axs[1].set_xlabel(r'$\\pi^*$')\n", "axs[1].set_title(r'optimal prior given $\\pi^*$')\n", "axs[1].legend()\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "8132e84a", "metadata": {}, "source": [ "## Was the Navy Captain’s Hunch Correct?\n", "\n", "We now compare average losses obtained by our frequentist Neyman-Pearson\n", "and Bayesian decision rules.\n", "\n", "As a starting point, let’s compare average loss functions when\n", "$ \\pi^{*}=0.5 $." ] }, { "cell_type": "code", "execution_count": null, "id": "6436bb0c", "metadata": { "hide-output": false }, "outputs": [], "source": [ "π_star = 0.5" ] }, { "cell_type": "code", "execution_count": null, "id": "c7207596", "metadata": { "hide-output": false }, "outputs": [], "source": [ "# frequentist\n", "V_fre_arr, PFA_arr, PD_arr = compute_V_fre(L0_arr, L1_arr, π_star, wf)\n", "\n", "# bayesian\n", "V_baye = π_star * V0 + (1 - π_star) * V1\n", "V_baye_bar = V_baye.min()" ] }, { "cell_type": "code", "execution_count": null, "id": "447e3938", "metadata": { "hide-output": false }, "outputs": [], "source": [ "plt.plot(range(T), V_fre_arr, label=r'$\\min_{d} \\overline{V}_{fre}(t,d)$')\n", "plt.plot([0, T], [V_baye_bar, V_baye_bar], label=r'$\\overline{V}_{baye}$')\n", "plt.xlabel('t')\n", "plt.title(r'$\\pi^*=0.5$')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "4c0d1cae", "metadata": {}, "source": [ "Evidently, there is no sample size $ t $ at which the Neyman-Pearson\n", "decision rule attains a lower loss function than does the Bayesian rule.\n", "\n", "Furthermore, the following graph indicates that the Bayesian decision\n", "rule does better on average for all values of $ \\pi^{*} $." ] }, { "cell_type": "code", "execution_count": null, "id": "256fee36", "metadata": { "hide-output": false }, "outputs": [], "source": [ "fig, axs = plt.subplots(1, 2, figsize=(14, 5))\n", "\n", "axs[0].plot(π_star_arr, V_fre_bar_arr, label=r'$\\overline{V}_{fre}$')\n", "axs[0].plot(π_star_arr, V_baye_bar_arr, label=r'$\\overline{V}_{baye}$')\n", "axs[0].legend()\n", "axs[0].set_xlabel(r'$\\pi^*$')\n", "\n", "axs[1].plot(π_star_arr, V_fre_bar_arr - V_baye_bar_arr, label='$diff$')\n", "axs[1].legend()\n", "axs[1].set_xlabel(r'$\\pi^*$')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "c23796db", "metadata": {}, "source": [ "The right panel of the above graph plots the difference\n", "$ \\bar{V}_{fre}-\\bar{V}_{Bayes} $.\n", "\n", "It is always positive." ] }, { "cell_type": "markdown", "id": "1fcce2c7", "metadata": {}, "source": [ "## More Details\n", "\n", "We can provide more insights by focusing on the case in which\n", "$ \\pi^{*}=0.5=\\pi_{0} $." ] }, { "cell_type": "code", "execution_count": null, "id": "e700c5e4", "metadata": { "hide-output": false }, "outputs": [], "source": [ "π_star = 0.5" ] }, { "cell_type": "markdown", "id": "87c74d81", "metadata": {}, "source": [ "Recall that when $ \\pi^*=0.5 $, the frequentist Neyman-Pearson decision rule sets a\n", "sample size `t_optimal` **ex ante**.\n", "\n", "For our parameter settings, we can compute its value:" ] }, { "cell_type": "code", "execution_count": null, "id": "e7d2716a", "metadata": { "hide-output": false }, "outputs": [], "source": [ "t_optimal" ] }, { "cell_type": "markdown", "id": "d3005e02", "metadata": {}, "source": [ "For convenience, let’s define `t_idx` as the Python array index\n", "corresponding to `t_optimal` sample size." ] }, { "cell_type": "code", "execution_count": null, "id": "3bcd640f", "metadata": { "hide-output": false }, "outputs": [], "source": [ "t_idx = t_optimal - 1" ] }, { "cell_type": "markdown", "id": "90a99703", "metadata": {}, "source": [ "## Distribution of Bayesian Decision Rule’s Time to Decide\n", "\n", "We use simulations to compute the frequency distribution of the time to\n", "decide for the Bayesian decision rule and compare that time to the\n", "frequentist rule’s fixed $ t $.\n", "\n", "The following Python code creates a graph that shows the frequency\n", "distribution of Bayesian times to decide of Bayesian decision maker,\n", "conditional on distribution $ q=f_{0} $ or $ q= f_{1} $\n", "generating the data.\n", "\n", "The blue and red dotted lines show averages for the Bayesian decision\n", "rule, while the black dotted line shows the frequentist optimal sample\n", "size $ t $.\n", "\n", "On average the Bayesian rule decides **earlier** than the frequentist\n", "rule when $ q= f_0 $ and **later** when $ q = f_1 $." ] }, { "cell_type": "code", "execution_count": null, "id": "786ace32", "metadata": { "hide-output": false }, "outputs": [], "source": [ "@jit(parallel=True)\n", "def check_results(L_arr, A, B, flag, π0):\n", "\n", " N, T = L_arr.shape\n", "\n", " time_arr = np.empty(N)\n", " correctness = np.empty(N)\n", "\n", " π_arr = π0 * L_arr / (π0 * L_arr + 1 - π0)\n", "\n", " for i in prange(N):\n", " for t in range(T):\n", " if (π_arr[i, t] < B) or (π_arr[i, t] > A):\n", " time_arr[i] = t + 1\n", " correctness[i] = (flag == 0 and π_arr[i, t] > A) or (flag == 1 and π_arr[i, t] < B)\n", " break\n", "\n", " return time_arr, correctness" ] }, { "cell_type": "code", "execution_count": null, "id": "4eccf07f", "metadata": { "hide-output": false }, "outputs": [], "source": [ "time_arr0, correctness0 = check_results(L0_arr, A, B, 0, π_star)\n", "time_arr1, correctness1 = check_results(L1_arr, A, B, 1, π_star)\n", "\n", "# unconditional distribution\n", "time_arr_u = np.concatenate((time_arr0, time_arr1))\n", "correctness_u = np.concatenate((correctness0, correctness1))" ] }, { "cell_type": "code", "execution_count": null, "id": "16df43c8", "metadata": { "hide-output": false }, "outputs": [], "source": [ "n1 = plt.hist(time_arr0, bins=range(1, 30), alpha=0.4, label='f0 generates')[0]\n", "n2 = plt.hist(time_arr1, bins=range(1, 30), alpha=0.4, label='f1 generates')[0]\n", "plt.vlines(t_optimal, 0, max(n1.max(), n2.max()), linestyle='--', label='frequentist')\n", "plt.vlines(np.mean(time_arr0), 0, max(n1.max(), n2.max()),\n", " linestyle='--', color='b', label='E(t) under f0')\n", "plt.vlines(np.mean(time_arr1), 0, max(n1.max(), n2.max()),\n", " linestyle='--', color='r', label='E(t) under f1')\n", "plt.legend();\n", "\n", "plt.xlabel('t')\n", "plt.ylabel('n')\n", "plt.title('Conditional frequency distribution of times')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "029bef1a", "metadata": {}, "source": [ "Later we’ll figure out how these distributions ultimately affect\n", "objective expected values under the Neyman-Pearson and Bayesian decision rules.\n", "\n", "To begin, let’s look at simulations of the Bayesian’s beliefs over time.\n", "\n", "We can compute updated beliefs at any time $ t $ using\n", "the one-to-one mapping from $ L_{t} $ to $ \\pi_{t} $ given\n", "$ \\pi_0 $ described in this lecture [Likelihood Ratio Processes](https://python.quantecon.org/likelihood_ratio_process.html)." ] }, { "cell_type": "code", "execution_count": null, "id": "5b4ef88f", "metadata": { "hide-output": false }, "outputs": [], "source": [ "π0_arr = π_star * L0_arr / (π_star * L0_arr + 1 - π_star)\n", "π1_arr = π_star * L1_arr / (π_star * L1_arr + 1 - π_star)" ] }, { "cell_type": "code", "execution_count": null, "id": "0d99d200", "metadata": { "hide-output": false }, "outputs": [], "source": [ "fig, axs = plt.subplots(1, 2, figsize=(14, 4))\n", "\n", "axs[0].plot(np.arange(1, π0_arr.shape[1]+1), np.mean(π0_arr, 0), label='f0 generates')\n", "axs[0].plot(np.arange(1, π1_arr.shape[1]+1), 1 - np.mean(π1_arr, 0), label='f1 generates')\n", "axs[0].set_xlabel('t')\n", "axs[0].set_ylabel(r'$E(\\pi_t)$ or ($1 - E(\\pi_t)$)')\n", "axs[0].set_title('Expectation of beliefs after drawing t observations')\n", "axs[0].legend()\n", "\n", "axs[1].plot(np.arange(1, π0_arr.shape[1]+1), np.var(π0_arr, 0), label='f0 generates')\n", "axs[1].plot(np.arange(1, π1_arr.shape[1]+1), np.var(π1_arr, 0), label='f1 generates')\n", "axs[1].set_xlabel('t')\n", "axs[1].set_ylabel(r'var($\\pi_t$)')\n", "axs[1].set_title('Variance of beliefs after drawing t observations')\n", "axs[1].legend()\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "83df04fb", "metadata": {}, "source": [ "The above figures compare averages and variances of updated Bayesian\n", "posteriors after $ t $ draws.\n", "\n", "The left graph compares $ E\\left(\\pi_{t}\\right) $ under\n", "$ f_{0} $ to $ 1-E\\left(\\pi_{t}\\right) $ under $ f_{1} $:\n", "they lie on top of each other.\n", "\n", "However, as the right hand side graph shows, there is significant\n", "difference in variances when $ t $ is small: the variance is lower\n", "under $ f_{1} $.\n", "\n", "The difference in variances is the reason that the Bayesian decision\n", "maker waits longer to decide when $ f_{1} $ generates the data.\n", "\n", "The code below plots outcomes of constructing an unconditional\n", "distribution by simply pooling the simulated data across the two\n", "possible distributions $ f_0 $ and $ f_1 $.\n", "\n", "The pooled distribution describes a sense in which on average the\n", "Bayesian decides earlier, an outcome that seems at least partly to\n", "confirm the Navy Captain’s hunch." ] }, { "cell_type": "code", "execution_count": null, "id": "94f28b1a", "metadata": { "hide-output": false }, "outputs": [], "source": [ "n = plt.hist(time_arr_u, bins=range(1, 30), alpha=0.4, label='bayesian')[0]\n", "plt.vlines(np.mean(time_arr_u), 0, n.max(), linestyle='--',\n", " color='b', label='bayesian E(t)')\n", "plt.vlines(t_optimal, 0, n.max(), linestyle='--', label='frequentist')\n", "plt.legend()\n", "\n", "plt.xlabel('t')\n", "plt.ylabel('n')\n", "plt.title('Unconditional distribution of times')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "120b41ef", "metadata": {}, "source": [ "## Probability of Making Correct Decision\n", "\n", "Now we use simulations to compute the fractions of samples in which the\n", "Bayesian and the frequentist Neyman-Pearson decision rules decide correctly.\n", "\n", "For the frequentist Neyman-Pearson rule, the probability of making the correct decision\n", "under $ f_{1} $ is the optimal probability of detection given\n", "$ t $ that we defined earlier, and similarly it equals $ 1 $\n", "minus the optimal probability of a false alarm under $ f_{0} $.\n", "\n", "Below we plot these two probabilities for the frequentist rule, along\n", "with the conditional probabilities that the Bayesian rule decides before\n", "$ t $ *and* that the decision is correct." ] }, { "cell_type": "code", "execution_count": null, "id": "af40b62c", "metadata": { "hide-output": false }, "outputs": [], "source": [ "# optimal PFA and PD of frequentist with optimal sample size\n", "V, PFA, PD = V_fre_t(t_optimal, L0_arr, L1_arr, π_star, wf)" ] }, { "cell_type": "code", "execution_count": null, "id": "51c77894", "metadata": { "hide-output": false }, "outputs": [], "source": [ "plt.plot([1, 20], [PD, PD], linestyle='--', label='PD: fre. chooses f1 correctly')\n", "plt.plot([1, 20], [1-PFA, 1-PFA], linestyle='--', label='1-PFA: fre. chooses f0 correctly')\n", "plt.vlines(t_optimal, 0, 1, linestyle='--', label='frequentist optimal sample size')\n", "\n", "N = time_arr0.size\n", "T_arr = np.arange(1, 21)\n", "plt.plot(T_arr, [np.sum(correctness0[time_arr0 <= t] == 1) / N for t in T_arr],\n", " label='q=f0 and baye. choose f0')\n", "plt.plot(T_arr, [np.sum(correctness1[time_arr1 <= t] == 1) / N for t in T_arr],\n", " label='q=f1 and baye. choose f1')\n", "plt.legend(loc=4)\n", "\n", "plt.xlabel('t')\n", "plt.ylabel('Probability')\n", "plt.title('Cond. probability of making correct decisions before t')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "e995611c", "metadata": {}, "source": [ "By averaging using $ \\pi^{*} $, we also plot the unconditional\n", "distribution." ] }, { "cell_type": "code", "execution_count": null, "id": "dcfb11bc", "metadata": { "hide-output": false }, "outputs": [], "source": [ "plt.plot([1, 20], [(PD + 1 - PFA) / 2, (PD + 1 - PFA) / 2],\n", " linestyle='--', label='fre. makes correct decision')\n", "plt.vlines(t_optimal, 0, 1, linestyle='--', label='frequentist optimal sample size')\n", "\n", "N = time_arr_u.size\n", "plt.plot(T_arr, [np.sum(correctness_u[time_arr_u <= t] == 1) / N for t in T_arr],\n", " label=\"bayesian makes correct decision\")\n", "plt.legend()\n", "\n", "plt.xlabel('t')\n", "plt.ylabel('Probability')\n", "plt.title('Uncond. probability of making correct decisions before t')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "35f95326", "metadata": {}, "source": [ "## Distribution of Likelihood Ratios at Neyman-Pearson’s $ t $\n", "\n", "Next we use simulations to construct distributions of likelihood ratios\n", "after $ t $ draws.\n", "\n", "To serve as useful reference points, we also show likelihood ratios that\n", "correspond to the Bayesian cutoffs $ A $ and $ B $.\n", "\n", "In order to exhibit the distribution more clearly, we report logarithms\n", "of likelihood ratios.\n", "\n", "The graphs below reports two distributions, one conditional on\n", "$ f_0 $ generating the data, the other conditional on $ f_1 $\n", "generating the data." ] }, { "cell_type": "code", "execution_count": null, "id": "6c447137", "metadata": { "hide-output": false }, "outputs": [], "source": [ "LA = (1 - π_star) * A / (π_star - π_star * A)\n", "LB = (1 - π_star) * B / (π_star - π_star * B)" ] }, { "cell_type": "code", "execution_count": null, "id": "1d0b142c", "metadata": { "hide-output": false }, "outputs": [], "source": [ "L_min = min(L0_arr[:, t_idx].min(), L1_arr[:, t_idx].min())\n", "L_max = max(L0_arr[:, t_idx].max(), L1_arr[:, t_idx].max())\n", "bin_range = np.linspace(np.log(L_min), np.log(L_max), 50)\n", "n0 = plt.hist(np.log(L0_arr[:, t_idx]), bins=bin_range, alpha=0.4, label='f0 generates')[0]\n", "n1 = plt.hist(np.log(L1_arr[:, t_idx]), bins=bin_range, alpha=0.4, label='f1 generates')[0]\n", "\n", "plt.vlines(np.log(LB), 0, max(n0.max(), n1.max()), linestyle='--', color='r', label='log($L_B$)')\n", "plt.vlines(np.log(LA), 0, max(n0.max(), n1.max()), linestyle='--', color='b', label='log($L_A$)')\n", "plt.legend()\n", "\n", "plt.xlabel('log(L)')\n", "plt.ylabel('n')\n", "plt.title('Cond. distribution of log likelihood ratio at frequentist t')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "2fdd3efa", "metadata": {}, "source": [ "The next graph plots the unconditional distribution of Bayesian times to\n", "decide, constructed as earlier by pooling the two conditional\n", "distributions." ] }, { "cell_type": "code", "execution_count": null, "id": "5395dd31", "metadata": { "hide-output": false }, "outputs": [], "source": [ "plt.hist(np.log(np.concatenate([L0_arr[:, t_idx], L1_arr[:, t_idx]])),\n", " bins=50, alpha=0.4, label='unconditional dist. of log(L)')\n", "plt.vlines(np.log(LB), 0, max(n0.max(), n1.max()), linestyle='--', color='r', label='log($L_B$)')\n", "plt.vlines(np.log(LA), 0, max(n0.max(), n1.max()), linestyle='--', color='b', label='log($L_A$)')\n", "plt.legend()\n", "\n", "plt.xlabel('log(L)')\n", "plt.ylabel('n')\n", "plt.title('Uncond. distribution of log likelihood ratio at frequentist t')\n", "\n", "plt.show()" ] } ], "metadata": { "date": 1770028421.9181511, "filename": "navy_captain.md", "kernelspec": { "display_name": "Python", "language": "python3", "name": "python3" }, "title": "Bayesian versus Frequentist Decision Rules" }, "nbformat": 4, "nbformat_minor": 5 }