{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## A Problem that stumped Milton Friedman\n", "\n", "(and that Abraham Wald solved by inventing sequential analysis)\n", "\n", "#### By [Chase Coleman](https://github.com/cc7768) and [Thomas J. Sargent](http://www.tomsargent.com/), edited by [Alberto Polo](https://github.com/albep) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We begin by importing some packages called by the code that we will be using in this notebook." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "using Distributions\n", "using QuantEcon.compute_fixed_point, QuantEcon.DiscreteRV, QuantEcon.draw, QuantEcon.LinInterp\n", "using Plots\n", "using LaTeXStrings" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sequential analysis\n", "\n", "Key ideas in play are:\n", "\n", " * Bayes' Law\n", " \n", " * Dynamic programming\n", "\n", " * type I and type II statistical errors\n", " \n", " * a type I error occurs when you reject a null hypothesis that is true\n", " \n", " * a type II error is when you accept a null hypothesis that is false \n", " \n", "\n", " * Abraham Wald's **sequential probability ratio test**\n", " \n", " * The **power** of a statistical test\n", " \n", " * The **critical region** of a statistical test\n", " \n", " * A **uniformly most powerful test**\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On pages 137-139 of his book **Two Lucky People** with Rose Friedman, Milton Friedman described a problem presented to him and Allen Wallis during World War II when they worked at the U.S. government's Statistical Research Group at Columbia University. \n", "\n", "Let's listen to Milton Friedman tell what happened.\n", "\n", " \"In order to understand the story, it is necessary to have an idea of a simple statistical problem, and of the\n", " standard procedure for dealing with it. The actual problem out of which sequential analysis grew will serve.\n", " The Navy has two alternative designs (say A and B) for a projectile. It wants to determine which is superior. \n", " To do so it undertakes a series of paired firings. On each round it assigns the value 1 or 0 to A accordingly as\n", " its performance is superior or inferior to that of B and conversely 0 or 1 to B. The Navy asks the statistician \n", " how to conduct the test and how to analyze the results. \n", " \n", " \"The standard statistical answer was to specify a number of firings (say 1,000) and a pair of percentages\n", " (e.g., 53% and 47%) and tell the client that if A receives a 1 in more than 53% of the firings, it can be regarded\n", " as superior; if it receives a 1 in fewer than 47%, B can be regarded as superior; if the percentage is between\n", " 47% and 53%, neither can be so regarded.\n", " \n", " \"When Allen Wallis was discussing such a problem with (Navy) Captain Garret L. Schyler, the captain objected that such a test, to quote from Allen's account, may prove wasteful. If a wise and seasoned ordnance officer like Schyler were on the premises, he would see after the first few thousand or even few hundred [rounds] that the experiment need not be completed either because the new method is obviously inferior or because it is obviously superior beyond what was hoped for $\\ldots$ \"\n", " \n", " Friedman and Wallis struggled with the problem but after realizing that they were not able to solve it themselves told Abraham Wald it. That started Wald on the path that led *Sequential Analysis*. We'll formulate the problem using dynamic programming.\n", "\n", "This started Wald on the path that led him to _Sequential Analysis_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Dynamic programming formulation\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following presentation of the problem closely follows Dmitri Berskekas's treatment in **Dynamic Programming and Stochastic Control**. \n", "\n", "An i.i.d. random variable $z$ can take on values \n", "\n", " * $z \\in [ v_1, v_2, \\ldots, v_n]$ when $z$ is a discrete-valued random variable\n", " \n", " * $ z \\in V$ when $z$ is a continuous random variable. \n", "\n", "A decision maker wants to know which of two probability distributions governs $z$. To formalize this idea,\n", "let $x \\in [x_0, x_1]$ be a hidden state that indexes the two distributions:\n", "\n", "$$ P(v_k \\mid x) = \\begin{cases} f_0(v_k) & \\mbox{if } x = x_0, \\\\\n", " f_1(v_k) & \\mbox{if } x = x_1. \\end{cases} $$\n", " \n", "when $z$ is a discrete random variable and a density\n", "\n", "\n", "$$ P(v \\mid x) = \\begin{cases} f_0(v) & \\mbox{if } x = x_0, \\\\\n", " f_1(v) & \\mbox{if } x = x_1. \\end{cases} $$\n", " \n", "when $v$ is continuously distributed. \n", "\n", " \n", "\n", "Before observing any outcomes, a decision maker believes that the probability that $x = x_0$ is $p_{-1}\\in (0,1)$: \n", "\n", "$$ p_{-1} = \\textrm{Prob}(x=x_0 \\mid \\textrm{ no observations}) $$\n", "\n", "After observing $k+1$ observations $z_k, z_{k-1}, \\ldots, z_0$ he believes that the probability that the distribution is $f_0$ is\n", "\n", "$$ p_k = {\\rm Prob} ( x = x_0 \\mid z_k, z_{k-1}, \\ldots, z_0) $$\n", "\n", "We can compute this $p_k$ recursively by applying Bayes' law:\n", "\n", "$$ p_0 = \\frac{ p_{-1} f_0(z_0)}{ p_{-1} f_0(z_0) + (1-p_{-1}) f_1(z_0) } $$\n", "\n", "and then\n", "\n", "$$ p_{k+1} = \\frac{ p_k f_0(z_{k+1})}{ p_k f_0(z_{k+1}) + (1-p_k) f_1 (z_{k+1}) }. $$\n", "\n", "\n", "After observing $z_k, z_{k-1}, \\ldots, z_0$, the decision maker believes that $z_{k+1}$ \n", "has probability distribution\n", "\n", "$$ p(z_{k+1}) = p_k f_0(z_{k+1}) + (1-p_k) f_1 (z_{k+1}). $$\n", "\n", "This is evidently a mixture of distributions $f_0$ and $f_1$, with the weight on $f_0$ being the posterior probability $f_0$ that the distribution is $f_0$. \n", "\n", "**Remark:** *Because the decision maker believes that $z_{k+1}$ is drawn from a mixture of two i.i.d. distributions, he does *not* believe that the sequence $[z_{k+1}, z_{k+2}, \\ldots] $ is i.i.d. Instead, he believes that it is *exchangeable*. See David Kreps\n", "*Notes on the Theory of Choice*, chapter 11, for a discussion.*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's look at some examples of two distributions. Here we'll display two beta distributions. First, we'll show the two distributions, then we'll show mixtures of these same two distributions with various mixing probabilities $p_k$." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Create two distributions over 50 values for k\n", "# We are using a discretized beta distribution\n", "\n", "p_m1 = linspace(0, 1, 50)\n", "f0 = clamp(pdf(Beta(1, 1), p_m1), 1e-8, Inf)\n", "f0 = f0 / sum(f0)\n", "f1 = clamp(pdf(Beta(9, 9), p_m1), 1e-8, Inf)\n", "f1 = f1 / sum(f1);" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = plot([f0 f1], \n", " xlabel=L\"$k$ Values\",\n", " ylabel=L\"Probability of $z_k$\",\n", " labels=reshape([L\"$f_0$\"; L\"$f_1$\"], 1, 2),\n", " linewidth=2,\n", " ylims=[0.;0.07],\n", " title=\"Original Distributions\")\n", "\n", "mix = Array(Float64, 50, 3)\n", "p_k = [0.25; 0.5; 0.75]\n", "labels = Array(String, 3)\n", "for i in 1:3\n", " mix[:, i] = p_k[i] * f0 + (1 - p_k[i]) * f1\n", " labels[i] = string(L\"$p_k$ = \", p_k[i])\n", "end\n", " \n", "b = plot(mix, \n", " xlabel=L\"$k$ Values\",\n", " ylabel=L\"Probability of $z_k$\",\n", " labels=reshape(labels, 1, 3),\n", " linewidth=2,\n", " ylims=[0.;0.06],\n", " title=\"Mixture of Original Distributions\")\n", "\n", "plot(a, b, layout=(2, 1), size=(600, 800))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Losses and costs\n", "\n", "\n", "After observing $z_k, z_{k-1}, \\ldots, z_0$, the decision maker chooses among three distinct actions:\n", "\n", "* He decides that $x = x_1$ and draws no more $z$'s\n", "\n", "* He decides that $x = x_0$ and draws no more $z$'s\n", "\n", "* He postpones deciding now and instead chooses to draw a $z_{k+1}$\n", "\n", "Associated with these three actions, the decision maker suffers three kinds of losses:\n", "\n", " \n", "* A loss $L_0$ if he decides $x = x_0$ when actually $x=x_1$\n", "\n", "* A loss $L_1$ if he decides $x = x_1$ when actually $x=x_0$\n", "\n", "* A cost $c$ if he postpones deciding and chooses instead to draw another $z$ \n", "\n", "For example, suppose that we regard $x=x_0$ as a null hypothesis. Then \n", "\n", "* We can think of $L_1$ as the loss associated with a type I error\n", "\n", "* We can think of $L_0$ as the loss associated with a type II error" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A Bellman equation\n", "\n", "Let $J_k(p_k)$ be the total loss for a decision maker who with posterior probability $p_k$ who chooses optimally.\n", "\n", "The loss functions $\\{J_k(p_k)\\}_k$ satisfy the Bellman equations\n", "\n", "$$ J_k(p_k) = \\min \\left[ (1-p_k) L_0, p_k L_1, c + E_{z_{k+1}} \\left\\{ J_{k+1} (p_{k+1} \\right\\} \\right] $$\n", "\n", "where $E_{z_{k+1}}$ denotes a mathematical expectation over the distribution of $z_{k+1}$ and the minimization is over the three actions, accept $x_0$, accept $x_1$, and postpone deciding and draw \n", "a $z_{k+1}$. \n", "\n", "Let \n", "\n", "$$ A_k(p_k) = E_{z_{k+1}} \\left\\{ J_{k+1} \\left[\\frac{ p_k f_0(z_{k+1})}{ p_k f_0(z_{k+1}) + (1-p_k) f_1 (z_{k+1}) } \\right] \\right\\} $$\n", "\n", "Then we can write out Bellman equation as\n", "\n", "$$ J_k(p_k) = \\min \\left[ (1-p_k) L_0, p_k L_1, c + A_k(p_k) \\right] $$\n", "\n", "where $p_k \\in [0,1]$. \n", "\n", "Evidently,the optimal decision rule is characterized by two numbers $\\alpha_k, \\beta_k \\in (0,1) \\times (0,1)$\n", "that satisfy\n", "\n", "$$ (1- p_k) L_0 < \\min p_k L_1, c + A_k(p_k) \\textrm { if } p_k \\geq \\alpha_k $$\n", "\n", "and \n", "\n", "$$ p_k L_1 < \\min (1-p_k) L_0, c + A_k(p_k) \\textrm { if } p_k \\leq \\beta_k $$\n", "\n", "The optimal decision rule is then\n", "\n", "$$ \\textrm { accept } x=x_0 \\textrm{ if } p_k \\geq \\alpha_k \\\\\n", " \\textrm { accept } x=x_1 \\textrm{ if } p_k \\leq \\beta_k \\\\\n", " \\textrm { draw another } z \\textrm{ if } \\beta_k \\leq p_k \\leq \\alpha_k $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Infinite horizon version\n", "\n", "An infinite horizon version of this problem is associated with the limiting Bellman equation \n", "\n", "$$ J(p_k) = \\min \\left[ (1-p_k) L_0, p_k L_1, c + A(p_k) \\right] \\quad (*) $$\n", "\n", "where\n", "\n", "$$ A(p_k) = E_{z_{k+1}} \\left\\{ J \\left[\\frac{ p_k f_0(z_{k+1})}{ p_k f_0(z_{k+1}) + (1-p_k) f_1 (z_{k+1}) } \\right] \\right\\} $$\n", "\n", "and again the minimization is over the three actions, accept $x_1$, accept $x_0$, and postpone deciding and draw \n", "a $z_{k+1}$.\n", "\n", "Here\n", "\n", " * $ (1-p_k) L_0$ is the expected loss associated with accepting $x_1$ (i.e., the cost of making a type I error)\n", " \n", " * $p_k L_1$ is the expected loss associated with accepting $x_0$ (i.e., the cost of making a type II error)\n", " \n", " * $ c + A(p_k)$ is the expected cost associated with drawing one more $z$\n", "\n", "\n", "Now the optimal decision rule is characterized by two probabilities $0 < \\beta < \\alpha < 1$ and \n", "\n", "\n", "$$ \\textrm { accept } x=x_0 \\textrm{ if } p_k \\geq \\alpha \\\\\n", " \\textrm { accept } x=x_1 \\textrm{ if } p_k \\leq \\beta \\\\\n", " \\textrm { draw another } z \\textrm{ if } \\beta \\leq p_k \\leq \\alpha $$\n", "\n", "\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One sensible approach is to write the three components of the value function that appears on the rights side of the Bellman equation as separate functions. Later, doing this will help us obey the don't repeat yourself (DRY) rule of coding. Here goes:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "\"\"\"\n", "For a given probability return expected loss of choosing model 0\n", "\"\"\"\n", "function expect_loss_choose_0(p, L0)\n", " return (1-p)*L0\n", "end\n", "\n", "\"\"\"\n", "For a given probability return expected loss of choosing model 1\n", "\"\"\"\n", "function expect_loss_choose_1(p, L1)\n", " return p*L1\n", "end\n", "\n", "\"\"\"\n", "We will need to be able to evaluate the expectation of our Bellman\n", "equation J. In order to do this, we need the current probability\n", "that model 0 is correct (p), the distributions (f0, f1), and a\n", "function that can evaluate the Bellman equation\n", "\"\"\"\n", "function EJ(p, f0, f1, J)\n", " # Get the current distribution we believe (p*f0 + (1-p)*f1)\n", " curr_dist = p*f0 + (1-p)*f1\n", " \n", " # Get tomorrow's expected distribution through Bayes law\n", " tp1_dist = clamp((p*f0) ./ (p*f0 + (1-p)*f1), 0, 1)\n", " \n", " # Evaluate the expectation\n", " EJ = dot(curr_dist, J.(tp1_dist))\n", " \n", " return EJ\n", "end\n", " \n", "expect_loss_cont(p, c, f0, f1, J) = c + EJ(p, f0, f1, J);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To approximate the solution of the Bellman equation (\\*) above, we can deploy a method known as value function iteration (iterating on the Bellman equation) on a grid of points. Because we are iterating on a grid, the current probability, $p_k$, is restricted to a set number of points. However, in order to evaluate the expectation of the Bellman equation for tomorrow, $A(p_{k})$, we must be able to evaluate at various $p_{k+1}$ which may or may not correspond with points on our grid. The way that we resolve this issue is by using *linear interpolation*. This means to evaluate $J(p)$ where $p$ is not a grid point, we must use two points: first, we use the largest of all the grid points smaller than $p$, and call it $p_i$, and, second, we use the grid point immediately after $p$, named $p_{i+1}$, to approximate the function value in the following manner:\n", "\n", "$$ J(p) = J(p_i) + (p - p_i) \\frac{J(p_{i+1}) - J(p_i)}{p_{i+1} - p_{i}}$$\n", "\n", "In one dimension, you can think of this as simply drawing a line between each pair of points on the grid.\n", "\n", "For more information on both linear interpolation and value function iteration methods, see the Quant-Econ [lecture](http://quant-econ.net/py/ifp.html) on the income fluctuation problem." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Compute iterate 5 with error 0.08552607733051265\n", "Compute iterate 10 with error 0.00038782894418165625\n", "Compute iterate 15 with error 1.6097835344730527e-6\n", "Converged in 16 steps\n" ] } ], "source": [ "\"\"\"\n", "Evaluates the value function for a given continuation value\n", "function; that is, evaluates\n", "\n", " J(p) = min(pL0, (1-p)L1, c + E[J(p')])\n", "\n", "Uses linear interpolation between points\n", "\"\"\"\n", "function bellman_operator(pgrid, c, f0, f1, L0, L1, J)\n", " m = length(pgrid)\n", " @assert m == length(J)\n", " \n", " J_out = zeros(m)\n", " J_interp = LinInterp(pgrid, J)\n", "\n", " for (p_ind, p) in enumerate(pgrid)\n", " # Payoff of choosing model 0\n", " p_c_0 = expect_loss_choose_0(p, L0)\n", " p_c_1 = expect_loss_choose_1(p, L1)\n", " p_con = expect_loss_cont(p, c, f0, f1, J_interp)\n", " \n", " J_out[p_ind] = min(p_c_0, p_c_1, p_con)\n", " end\n", " \n", " return J_out\n", "end\n", "\n", "# To solve\n", "pg = linspace(0, 1, 251)\n", "J1 = compute_fixed_point(x -> bellman_operator(pg, 0.5, f0, f1, 5.0, 5.0, x),\n", " zeros(length(pg)), err_tol=1e-6, print_skip=5);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now for some gentle criticisms of the preceding code. Although it works fine, by writing the code in terms of functions, we have to pass around some things that are constant throughout the problem, i.e., $c$, $f_0$, $f_1$, $L_0$, and $L_1$. \n", "\n", "Now that we have a working script, let's turn it into types and type-specific methods for each function. This will allow us to simplify the function calls and make the code more reusable.\n", "\n", "So to illustrate a good alternative approach, we write a type that stores all of our parameters for us internally, as well as the model solution, and is used by many of the same functions that we used above. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "stopping_dist" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"\n", "This type is used to store the solution to the problem presented \n", "in the \"Wald Friedman\" notebook presented on the QuantEcon website.\n", "\n", "Solution\n", "----------\n", "J : vector(Float64)\n", " Discretized value function that solves the Bellman equation\n", "lb : scalar(Real)\n", " Lower cutoff for continuation decision\n", "ub : scalar(Real)\n", " Upper cutoff for continuation decision\n", "\"\"\"\n", "type WFSolution\n", " J::Vector{Float64}\n", " lb::Real\n", " ub::Real\n", "end\n", "\n", "\"\"\"\n", "This type is used to solve the problem presented in the \"Wald Friedman\"\n", "notebook presented on the QuantEcon website.\n", "\n", "Parameters\n", "----------\n", "c : scalar(Real)\n", " Cost of postponing decision\n", "L0 : scalar(Real)\n", " Cost of choosing model 0 when the truth is model 1\n", "L1 : scalar(Real)\n", " Cost of choosing model 1 when the truth is model 0\n", "f0 : vector(Float64)\n", " A finite state probability distribution\n", "f1 : vector(Float64)\n", " A finite state probability distribution\n", "m : scalar(Int64)\n", " Number of points to use in function approximation\n", "\"\"\"\n", "immutable WaldFriedman\n", " c::Real\n", " L0::Real\n", " L1::Real\n", " f0::Vector{Float64}\n", " f1::Vector{Float64}\n", " m::Int64\n", " pgrid::LinSpace{Float64}\n", " sol::WFSolution\n", "end\n", "\n", "function WaldFriedman(c, L0, L1, f0, f1; m=25)\n", " pgrid = linspace(0.0, 1.0, m)\n", "\n", " # Renormalize distributions so nothing is \"too\" small\n", " f0 = clamp(f0, 1e-8, 1-1e-8)\n", " f1 = clamp(f1, 1e-8, 1-1e-8)\n", " f0 = f0 / sum(f0)\n", " f1 = f1 / sum(f1)\n", " J = zeros(m)\n", " lb = 0.\n", " ub = 0.\n", " \n", " WaldFriedman(c, L0, L1, f0, f1, m, pgrid, WFSolution(J, lb, ub))\n", "end\n", "\n", "\"\"\"\n", "This function takes a value for the probability with which\n", "the correct model is model 0 and returns the mixed\n", "distribution that corresponds with that belief.\n", "\"\"\"\n", "function current_distribution(wf::WaldFriedman, p)\n", " return p*wf.f0 + (1-p)*wf.f1\n", "end\n", "\n", "\"\"\"\n", "This function takes a value for p, and a realization of the\n", "random variable and calculates the value for p tomorrow.\n", "\"\"\"\n", "function bayes_update_k(wf::WaldFriedman, p, k)\n", " f0_k = wf.f0[k]\n", " f1_k = wf.f1[k]\n", "\n", " p_tp1 = p*f0_k / (p*f0_k + (1-p)*f1_k)\n", "\n", " return clamp(p_tp1, 0, 1)\n", "end\n", "\n", "\"\"\"\n", "This is similar to `bayes_update_k` except it returns a\n", "new value for p for each realization of the random variable\n", "\"\"\"\n", "function bayes_update_all(wf::WaldFriedman, p)\n", " return clamp(p*wf.f0 ./ (p*wf.f0 + (1-p)*wf.f1), 0, 1)\n", "end\n", "\n", "\"\"\"\n", "For a given probability specify the cost of accepting model 0\n", "\"\"\"\n", "function payoff_choose_f0(wf::WaldFriedman, p)\n", " return (1-p)*wf.L0\n", "end\n", "\n", "\"\"\"\n", "For a given probability specify the cost of accepting model 1\n", "\"\"\"\n", "function payoff_choose_f1(wf::WaldFriedman, p)\n", " return p*wf.L1\n", "end\n", "\n", "\"\"\"\n", "This function evaluates the expectation of the value function\n", "at period t+1. It does so by taking the current probability\n", "distribution over outcomes:\n", "\n", " p(z_{k+1}) = p_k f_0(z_{k+1}) + (1-p_k) f_1(z_{k+1})\n", "\n", "and evaluating the value function at the possible states\n", "tomorrow J(p_{t+1}) where\n", "\n", " p_{t+1} = p f0 / ( p f0 + (1-p) f1)\n", "\n", "Parameters\n", "----------\n", "p : scalar\n", " The current believed probability that model 0 is the true\n", " model.\n", "J : function (interpolant)\n", " The current value function for a decision to continue\n", "\n", "Returns\n", "-------\n", "EJ : scalar\n", " The expected value of the value function tomorrow\n", "\"\"\"\n", "function EJ(wf::WaldFriedman, p, J)\n", " # Pull out information\n", " f0, f1 = wf.f0, wf.f1\n", "\n", " # Get the current believed distribution and tomorrows possible dists\n", " # Need to clip to make sure things don't blow up (go to infinity)\n", " curr_dist = current_distribution(wf, p)\n", " tp1_dist = bayes_update_all(wf, p)\n", " \n", " # Evaluate the expectation\n", " EJ = dot(curr_dist, J.(tp1_dist))\n", "\n", " return EJ\n", "end\n", "\n", "\"\"\"\n", "For a given probability distribution and value function give\n", "cost of continuing the search for correct model\n", "\"\"\"\n", "function payoff_continue(wf::WaldFriedman, p, J)\n", " return wf.c + EJ(wf, p, J)\n", "end\n", "\n", "\"\"\"\n", "Evaluates the value function for a given continuation value\n", "function; that is, evaluates\n", "\n", " J(p) = min( (1-p)L0, pL1, c + E[J(p')])\n", "\n", "Uses linear interpolation between points\n", "\"\"\"\n", "function bellman_operator(wf::WaldFriedman, J)\n", " c, L0, L1, f0, f1 = wf.c, wf.L0, wf.L1, wf.f0, wf.f1\n", " m, pgrid = wf.m, wf.pgrid\n", "\n", " J_out = zeros(m)\n", " J_interp = LinInterp(pgrid, J)\n", " \n", " for (p_ind, p) in enumerate(pgrid)\n", " # Payoff of choosing model 0\n", " p_c_0 = payoff_choose_f0(wf, p)\n", " p_c_1 = payoff_choose_f1(wf, p)\n", " p_con = payoff_continue(wf, p, J_interp)\n", "\n", " J_out[p_ind] = min(p_c_0, p_c_1, p_con)\n", " end\n", "\n", " return J_out\n", "end\n", "\n", "\"\"\"\n", "This function takes a value function and returns the corresponding\n", "cutoffs of where you transition between continue and choosing a\n", "specific model\n", "\"\"\"\n", "function find_cutoff_rule(wf::WaldFriedman, J)\n", " m, pgrid = wf.m, wf.pgrid\n", "\n", " # Evaluate cost at all points on grid for choosing a model\n", " p_c_0 = payoff_choose_f0(wf, pgrid)\n", " p_c_1 = payoff_choose_f1(wf, pgrid)\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", " lb = pgrid[searchsortedlast(p_c_1 - J, 1e-10)]\n", " ub = pgrid[searchsortedlast(J - p_c_0, -1e-10)]\n", "\n", " return lb, ub\n", "end\n", "\n", "function solve_model(wf; tol=1e-7)\n", " bell_op(x) = bellman_operator(wf, x)\n", " J = compute_fixed_point(bell_op, zeros(wf.m), err_tol=tol, print_skip=5)\n", "\n", " wf.sol.J = J\n", " wf.sol.lb, wf.sol.ub = find_cutoff_rule(wf, J)\n", " return J\n", "end\n", "\n", "\"\"\"\n", "This function takes an initial condition and simulates until it\n", "stops (when a decision is made).\n", "\"\"\"\n", "function simulate(wf::WaldFriedman, f; p0=0.5)\n", " # Check whether vf is computed\n", " if sumabs(wf.sol.J) < 1e-8\n", " solve_model(wf)\n", " end\n", " \n", " # Unpack useful info\n", " lb, ub = wf.sol.lb, wf.sol.ub\n", " drv = DiscreteRV(f)\n", "\n", " # Initialize a couple useful variables\n", " decision_made = false\n", " decision = 0\n", " p = p0\n", " t = 0\n", "\n", " while !decision_made\n", " # Maybe should specify which distribution is correct one so that\n", " # the draws come from the \"right\" distribution\n", " k = draw(drv)[1]\n", " t = t+1\n", " p = bayes_update_k(wf, p, k)\n", " if p < lb\n", " decision_made = true\n", " decision = 1\n", " elseif p > ub\n", " decision_made = true\n", " decision = 0\n", " end\n", " end\n", " \n", " return decision, p, t\n", "end\n", "\n", "\"\"\"\n", "Uses the distribution f0 as the true data generating\n", "process\n", "\"\"\"\n", "function simulate_tdgp_f0(wf::WaldFriedman; p0=0.5)\n", " decision, p, t = simulate(wf, wf.f0; p0=p0)\n", "\n", " if decision == 0\n", " correct = true\n", " else\n", " correct = false\n", " end\n", " \n", " return correct, p, t\n", "end\n", "\n", "\"\"\"\n", "Uses the distribution f1 as the true data generating\n", "process\n", "\"\"\"\n", "function simulate_tdgp_f1(wf::WaldFriedman; p0=0.5)\n", " decision, p, t = simulate(wf, wf.f1; p0=p0)\n", "\n", " if decision == 1\n", " correct = true\n", " else\n", " correct = false\n", " end\n", " \n", " return correct, p, t\n", "end\n", "\n", "\"\"\"\n", "Simulates repeatedly to get distributions of time needed to make a\n", "decision and how often they are correct.\n", "\"\"\"\n", "function stopping_dist(wf::WaldFriedman; ndraws=250, tdgp=\"f0\")\n", " if tdgp==\"f0\"\n", " simfunc = simulate_tdgp_f0\n", " else\n", " simfunc = simulate_tdgp_f1\n", " end\n", " \n", " # Allocate space\n", " tdist = Array(Int64, ndraws)\n", " cdist = Array(Bool, ndraws)\n", "\n", " for i in 1:ndraws\n", " correct, p, t = simfunc(wf)\n", " tdist[i] = t\n", " cdist[i] = correct\n", " end\n", " \n", " return cdist, tdist\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's use our types to solve the Bellman equation (*) and check whether it gives the same answer attained above." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Compute iterate 5 with error 0.0855260926408965\n", "Compute iterate 10 with error 0.0003878288254588469\n", "Compute iterate 15 with error 1.6097831208039537e-6\n", "Converged in 16 steps\n", "If this is true then both approaches gave same answer:\n", "true" ] } ], "source": [ "wf = WaldFriedman(0.5, 5.0, 5.0, f0, f1; m=251)\n", "J2 = compute_fixed_point(x -> bellman_operator(wf, x), zeros(wf.m), err_tol=1e-6, print_skip=5)\n", "\n", "@printf(\"If this is true then both approaches gave same answer:\\n\")\n", "print(isapprox(J1, J2; atol=1e-5))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Numerical Example\n", "\n", "Now let's specify the two probability distibutions (the ones that we plotted earlier)\n", "\n", "* for $f_0$ we'll assume a beta distribution with parameters $a=1, b=1$\n", "\n", "* for $f_1$ we'll assume a beta distribution with parameters $a=9, b=9$\n", "\n", "The density of a beta probability distribution with parameters $a$ and $b$ is\n", "\n", "$$ f(z; a, b) = \\frac{\\Gamma(a+b) z^{a-1} (1-z)^{b-1}}{\\Gamma(a) \\Gamma(b)}$$\n", "\n", "where $\\Gamma$ is the gamma function \n", "\n", "$$\\Gamma(t) = \\int_{0}^{\\infty} x^{t-1} e^{-x} dx$$\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Choose parameters\n", "c = 1.25\n", "L0 = 27.0\n", "L1 = 27.0\n", "\n", "# Choose n points and distributions\n", "m = 251\n", "f0 = pdf(Beta(2.5, 3), linspace(0, 1, m))\n", "f0 = f0 / sum(f0)\n", "f1 = pdf(Beta(3, 2.5), linspace(0, 1, m))\n", "f1 = f1 / sum(f1) # Make sure sums to 1\n", "\n", "# Create an instance of our WaldFriedman class\n", "wf = WaldFriedman(c, L0, L1, f0, f1; m=m);" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Compute iterate 5 with error 1.2384971736003685\n", "Compute iterate 10 with error 0.6235689198598084\n", "Compute iterate 15 with error 0.03178165128978527\n", "Compute iterate 20 with error 0.0005980373085616719\n", "Compute iterate 25 with error 1.1172556856564597e-5\n", "Compute iterate 30 with error 2.0872615458245036e-7\n", "Converged in 31 steps\n" ] }, { "data": { "text/plain": [ "(Bool[true,true,true,false,true,false,false,true,false,true … false,false,false,true,false,true,true,true,true,true],[1,1,2,1,1,3,3,1,1,1 … 1,1,1,1,1,1,1,1,1,1])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Solve and simulate the solution\n", "cdist, tdist = stopping_dist(wf; ndraws=5000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The value function equals $ p L_1$ for $p \\leq \\alpha$, and $(1-p )L_0$ for $ p \\geq \\beta$.\n", "Thus, the slopes of the two linear pieces of the value function are determined by $L_1$ and \n", "$- L_0$. \n", "\n", "The value function is smooth in the interior region in which the probability assigned to distribution $f_0$ is in the indecisive region $p \\in (\\alpha, \\beta)$.\n", "\n", "The decision maker continues to sample until the probability that he attaches to model $f_0$ falls below $\\alpha$ or above $\\beta$.\n", "\n", "The value function is smooth in the interior region in which the probability assigned to distribution $f_0$ is in the indecisive region $p \\in (\\alpha, \\beta)$.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now to have some fun, you can change the cost parameters $L_0, L_1, c$, the parameters of two beta distributions $f_0$ and $f_1$, and the number of points and linear functions $m$ to use in our piece-wise continuous approximation to the value function. You can see the effects on the smoothness of the value function in the middle range as you increase the numbers of functions in the piecewise linear approximation. \n", "\n", "The function `stopping_dist` draws a number of simulations from $f_0$, computes a distribution of waiting times to making a decision, and displays a histogram of correct and incorrect decisions. (Here the correct decision occurs when $p_k$ eventually exceeds $\\beta$)." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = plot([f0 f1], \n", " xlabel=L\"$k$ Values\",\n", " ylabel=L\"Probability of $z_k$\",\n", " labels=reshape([L\"$f_0$\"; L\"$f_1$\"], 1, 2),\n", " linewidth=2,\n", " title=\"Distributions over Outcomes\")\n", "\n", "b = plot(wf.pgrid, wf.sol.J, \n", " xlabel=L\"$p_k$\",\n", " ylabel=\"Value of Bellman\",\n", " linewidth=2,\n", " title=\"Bellman Equation\")\n", " plot!(fill(wf.sol.lb, 2), [minimum(wf.sol.J); maximum(wf.sol.J)],\n", " linewidth=2, color=:black, linestyle=:dash, label=\"\", ann=(wf.sol.lb-0.05, 5., L\"\\beta\"))\n", " plot!(fill(wf.sol.ub, 2), [minimum(wf.sol.J); maximum(wf.sol.J)],\n", " linewidth=2, color=:black, linestyle=:dash, label=\"\", ann=(wf.sol.ub+0.02, 5., L\"\\alpha\"),\n", " legend=:none)\n", "\n", "counts = Array(Int64, maximum(tdist))\n", "for i in 1:maximum(tdist)\n", " counts[i] = sum(tdist .== i)\n", "end\n", "c = bar(counts,\n", " xticks=0:1:maximum(tdist),\n", " xlabel=\"Time\",\n", " ylabel=\"Frequency\",\n", " title=\"Stopping Times\",\n", " legend=:none)\n", "\n", "counts = Array(Int64, 2)\n", "for i in 1:2\n", " counts[i] = sum(cdist .== i-1)\n", "end\n", "d = bar([0; 1],\n", " counts, \n", " xticks=[0; 1],\n", " title=\"Correct Decisions\", \n", " ann=(-.4, 0.6 * sum(cdist), \"Percent Correct = $(sum(cdist)/length(cdist))\"),\n", " legend=:none)\n", "\n", "plot(a, b, c, d, layout=(2, 2), size=(1200, 800))" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "### Comparison with Neyman-Pearson formulation\n", "\n", "For several reasons, it is useful to describe the theory underlying test that Navy Captain G. S. Schuyler had been told to use and that led him to approach Milton Friedman and Allan Wallis to convey his conjecture that superior practical procedures existed. Evidently,\n", "the Navy had told Captail Schuyler to use what it knew to be a state-of-the-art Neyman-Pearson test. \n", "\n", "We'll rely on Abraham Wald's elegant summary of Neyman-Pearson theory. For our purposes, watch for there features of the setup:\n", "\n", " * the assumption of a *fixed *sample size $n$\n", " \n", " * the application of laws of large numbers, conditioned on alternative probability models, to interpret the probabilities $\\alpha$ and $\\beta$ defined in the Neyman-Pearson theory.\n", " \n", "Recall that in the sequential analytic formulation above, that\n", "\n", " * The sample size $n$ is not fixed but rather an object to be chosen; technically $n$ is a random variable. \n", " \n", " * The parameters $\\alpha$ and $\\beta$ characterize cut-off rules used to determine $n$ as a random variable.\n", " \n", " * Laws of large numbers make no appearances in the sequential construction.\n", "\n", "In chapter 1 of **Sequential Analysis** Abraham Wald summarizes the Neyman-Pearson approach to hypothesis testing.\n", "\n", "Wald frames the problem as making a decision about a probability distribution that is partially known. (You have to assume that *something* is already known in order to state a well posed problem. Usually, *something* means *a lot*.)\n", "\n", "By restricting what is unknown, Wald uses the following simple structure to illustrate the main ideas.\n", "\n", " * a decision maker wants to decide which of two distributions $f_0$, $f_1$ govern an i.i.d. random variable $z$\n", " \n", " * The null hypothesis $H_0$ is the statement that $f_0$ governs the data.\n", " \n", " * The alternative hypothesis $H_1$ is the statement that $f_1$ governs the data. \n", " \n", " * The problem is to devise and analyze a test of hypthothesis $H_0$ against the alternative hypothesis $H_1$ on the basis of a sample of a fixed number $n$ independent observations $z_1, z_2, \\ldots, z_n$ of the random variable $z$. \n", " \n", "To quote Abraham Wald,\n", "\n", " * *A test procedure leading to the acceptance or rejection of the hypothesis in question is simply a rule specifying, for each possible sample of size $n$, whether the hypothesis should be accepted or rejected on the basis of the sample. This may also be expressed as follows: A test procedure is simply a subdivision of the totality of all possibsle samples of size $n$ into two mutually exclusive parts, say part 1 and part 2, together with the application of the rule that the hypothesis be accepted if the observed sample is contained in part 2. Part 1 is also called the critical region. Since part 2 is the totality of all samples of size 2 which are not included in part 1, part 2 is uniquely determined by part 1. Thus, choosing a test procedure is equivalent to determining a critical region.* \n", " \n", "Let's listen to Wald longer:\n", "\n", " * *As a basis for choosing among critical regions the following considerations have been advanced by Neyman and Pearson: In accepting or rejecting $H_0$ we may commit errors of two kinds. We commit an error of the first kind if we reject $H_0$ when it is true; we commit an error of the second kind if we accept $H_0$ whe $H_1$ is true. After a particular critical region $W$ has been chosen, the probability of committing an error of the first kind, as well as the probability of committing an error of the second kind is uniquely determined. The probability of committing an error of the first kind is equal to the probability, determined by the assumption that $H_0$ is true, that the observed sample will be included in the critical region $W$. The probability of committing an error of the second kind is equal to the probability, determined on the assumption that $H_1$ is true, that the probability will fall outside the critical region $W$. For any given critical region $W$ we shall denote the probability of an error of the first kind by $\\alpha$ and the probability of an error of the second kind by $\\beta$.*\n", "\n", "\n", "Let's listen carefully to how Wald applies a law of large numbers to interpret $\\alpha$ and $\\beta$:\n", "\n", " * *The probabilities $\\alpha$ and $\\beta$ have the following important practical interpretation: Suppose that we draw a large number of samples of size $n$. Let $M$ be the number of such samples drawn. Suppose that for each of these $M$ samples we reject $H_0$ if the sample is included in $W$ and accept $H_0$ if the sample lies outside $W$. In this way we make $M$ statements of rejection or acceptance. Some of these statements will in general be wrong. If $H_0$ is true and if $M$ is large, the probability is nearly $1$ (i.e., it is practically certain) that the proportion of wrong statements (i.e., the number of wrong statements divided by $M$) will be approximately $\\alpha$. If $H_1$ is true, the probability is nearly $1$ that the proportion of wrong statements will be approximately $\\beta$. Thus, we can say that in the long run [ here Wald applies a law of large numbers by driving $M \\rightarrow \\infty$ (our comment, not Wald's) ] the proportion of wrong statements will be $\\alpha$ if $H_0$is true and $\\beta $ if $H_1$ is true.*\n", "\n", "The quantity $\\alpha$ is called the *size* of the critical region, and the quantity $1-\\beta$ is called the *power* of the critical region.\n", "\n", "Wald notes that \n", "\n", " * *one critical region $W$ is more desirable than another if it has smaller values of $\\alpha$ and $\\beta$. Although either $\\alpha$ or $\\beta$ can be made arbitrarily small by a proper choice of the critical region $W$, it is possible to make both $\\alpha$ and $\\beta$ arbitrarily small for a fixed value of $n$, i.e., a fixed sample size.*\n", "\n", "\n", "Wald summarizes Neyman and Pearson's setup as follows:\n", "\n", " * *Neyman and Pearson show that a region consisting of all samples $(z_1, z_2, \\ldots, z_n)$ which satisfy the inequality \n", " $$\\frac{ f_1(z_1) \\cdots f_1(z_n)}{f_0(z_1) \\cdots f_1(z_n)} \\geq k $$ is a most powerful critical region for testing the hypothesis $H_0$ against the alternative hypothesis $H_1$. The term $k$ on the right side is a constant chosen so that the region will have the required size $\\alpha$.*\n", " \n", "Wald goes on to discuss Neyman and Pearson's concept of *uniformly most powerful* test. \n", "\n", "\n", "Here is how Wald introduces the notion of a sequential test\n", "•A rule is given for making one of the following three decisions at any stage of the experiment (at the m th trial for each integral value of m ): (1) to accept the hypothesis H , (2) to reject the hypothesis H , (3) to continue the experiment by making an additional observation. Thus, such a test procedure is carried out sequentially. On the basis of the first observation one of the aforementioned decisions is made. If the first or second decision is made, the process is terminated. If the third decision is made, a second trial is performed. Again, on the basis of the first two observations one of the three decisions is made. If the third decision is made, a third trial is performed, and so on. The process is continued until either the first or the second decisions is made. The number n of observations required by such a test procedure is a random variable, since the value of n depends on the outcome of the observations.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.5.1-pre", "language": "julia", "name": "julia-0.5" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.5.1" } }, "nbformat": 4, "nbformat_minor": 1 }