likelihood_bayes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e5b92567",
   "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",
   "id": "51261eb1",
   "metadata": {},
   "source": [
    "# Likelihood Ratio Processes and Bayesian Learning"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "325de6c4",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "This lecture describes the role that **likelihood ratio processes** play in  **Bayesian learning**.\n",
    "\n",
    "As in [this lecture](https://python.quantecon.org/likelihood_ratio_process.html), we’ll use a simple statistical setting from [this lecture](https://python.quantecon.org/exchangeable.html).\n",
    "\n",
    "We’ll focus on how a likelihood ratio process and a **prior** probability determine a **posterior** probability.\n",
    "\n",
    "We’ll derive a convenient recursion for today’s posterior as a function of yesterday’s posterior and\n",
    "today’s multiplicative increment to a likelihood process.\n",
    "\n",
    "We’ll also present a useful generalization of that formula that represents today’s posterior in terms of an initial prior and\n",
    "today’s realization of the likelihood ratio process.\n",
    "\n",
    "We’ll study how, at least  in our setting, a Bayesian eventually learns the probability distribution that generates the data, an outcome that\n",
    "rests on the asymptotic behavior of likelihood ratio processes studied in [this lecture](https://python.quantecon.org/likelihood_ratio_process.html).\n",
    "\n",
    "We’ll also drill down into the psychology of our Bayesian learner and study dynamics  under his subjective beliefs.\n",
    "\n",
    "This lecture provides technical results that underly outcomes to be studied in [this lecture](https://python.quantecon.org/odu.html)\n",
    "and [this lecture](https://python.quantecon.org/wald_friedman.html) and [this lecture](https://python.quantecon.org/navy_captain.html).\n",
    "\n",
    "We’ll begin by loading some Python modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4db3a96",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams[\"figure.figsize\"] = (11, 5)  #set default figure size\n",
    "import numpy as np\n",
    "from numba import vectorize, jit, prange\n",
    "from math import gamma\n",
    "import pandas as pd\n",
    "\n",
    "import seaborn as sns\n",
    "colors = sns.color_palette()\n",
    "\n",
    "@jit\n",
    "def set_seed():\n",
    "    np.random.seed(142857)\n",
    "set_seed()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7d062f3",
   "metadata": {},
   "source": [
    "## The Setting\n",
    "\n",
    "We begin by reviewing the setting in [this lecture](https://python.quantecon.org/likelihood_ratio_process.html), which we adopt here too.\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 $ f $ or from $ g $.\n",
    "\n",
    "We will sometimes let $ q $ be the density that nature chose once and for all, so\n",
    "that $ q $ is either $ f $ or $ g $, permanently.\n",
    "\n",
    "Nature knows which density it permanently draws from, but we the observers do not.\n",
    "\n",
    "We do know both $ f $ and $ g $, but we don’t know which density nature\n",
    "chose.\n",
    "\n",
    "But we want to know.\n",
    "\n",
    "To do that, we use observations.\n",
    "\n",
    "We observe a sequence $ \\{w_t\\}_{t=1}^T $ of $ T $ IID draws\n",
    "from either $ f $ or $ g $.\n",
    "\n",
    "We want to use these observations to infer whether nature chose $ f $ or\n",
    "$ g $.\n",
    "\n",
    "A **likelihood ratio process** is a useful tool for this task.\n",
    "\n",
    "To begin, we define 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 the random variable $ W $.\n",
    "\n",
    "That means that under the $ g $ density,  $ \\ell (w_t)=\n",
    "\\frac{f\\left(w_{t}\\right)}{g\\left(w_{t}\\right)} $\n",
    "is evidently a nonnegative  random variable with mean $ 1 $.\n",
    "\n",
    "A **likelihood ratio process** for sequence\n",
    "$ \\left\\{ w_{t}\\right\\} _{t=1}^{\\infty} $ is defined as\n",
    "\n",
    "$$\n",
    "L\\left(w^{t}\\right)=\\prod_{i=1}^{t} \\ell (w_i),\n",
    "$$\n",
    "\n",
    "where $ w^t=\\{ w_1,\\dots,w_t\\} $ is a history of\n",
    "observations up to and including time $ t $.\n",
    "\n",
    "Sometimes for shorthand we’ll write $ L_t =  L(w^t) $.\n",
    "\n",
    "Notice that the likelihood process satisfies the *recursion* or\n",
    "*multiplicative decomposition*\n",
    "\n",
    "$$\n",
    "L(w^t) = \\ell (w_t) L (w^{t-1}) .\n",
    "$$\n",
    "\n",
    "The likelihood ratio and its logarithm are key tools for making\n",
    "inferences using a classic frequentist approach due to Neyman and\n",
    "Pearson [[Neyman and Pearson, 1933](https://python.quantecon.org/zreferences.html#id255)].\n",
    "\n",
    "We’ll again deploy the following Python code from [this lecture](https://python.quantecon.org/likelihood_ratio_process.html) that\n",
    "evaluates $ f $ and $ g $ as two different\n",
    "beta distributions, then computes and simulates an associated likelihood\n",
    "ratio process by generating a sequence $ w^t $ from *some*\n",
    "probability distribution, for example, a sequence of  IID draws from $ g $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "081ec7bf",
   "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 = jit(lambda x: p(x, F_a, F_b))\n",
    "g = jit(lambda x: p(x, G_a, G_b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2644ffeb",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@jit\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",
   "id": "ecaaf5b4",
   "metadata": {},
   "source": [
    "We’ll also use the following Python code to prepare some informative simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22cfadcd",
   "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": "code",
   "execution_count": null,
   "id": "5b87b7a2",
   "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": "markdown",
   "id": "126c3a86",
   "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{54.1}\n",
    "$$\n",
    "\n",
    "with $ \\pi_{0} $ being a Bayesian prior probability that $ q = f $,\n",
    "i.e., a personal or subjective belief about $ q $ based on our having seen no data.\n",
    "\n",
    "Below we define a Python function that updates belief $ \\pi $ using\n",
    "likelihood ratio $ \\ell $ according to  recursion [(54.1)](#equation-eq-recur1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ce92cf8",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@jit\n",
    "def update(π, l):\n",
    "    \"Update π using likelihood l\"\n",
    "\n",
    "    # Update belief\n",
    "    π = π * l / (π * l + 1 - π)\n",
    "\n",
    "    return π"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f022f504",
   "metadata": {},
   "source": [
    "Formula [(54.1)](#equation-eq-recur1) can be generalized  by iterating on it and thereby deriving 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{54.2}\n",
    "$$\n",
    "\n",
    "Formula [(54.2)](#equation-eq-bayeslaw103) generalizes formula [(54.1)](#equation-eq-recur1).\n",
    "\n",
    "Formula [(54.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 [(54.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 [(54.2)](#equation-eq-bayeslaw103) is the foundation for the insight that, because of how 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 $ \\pi_{0} $.\n",
    "\n",
    "First, we tell Python two values of $ \\pi_0 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd3d90ea",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "π1, π2 = 0.2, 0.8"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80d88730",
   "metadata": {},
   "source": [
    "Next we generate paths of the likelihood ratio process $ L_t $ and the posterior $ \\pi_t $ for a\n",
    "history of IID draws from density $ f $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7402307",
   "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,
   "id": "979b8754",
   "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",
   "id": "6367d2f1",
   "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,
   "id": "c91913fc",
   "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,
   "id": "477a32c3",
   "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",
   "id": "56a43373",
   "metadata": {},
   "source": [
    "Below we offer Python code that verifies that nature chose permanently to draw from density $ f $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e762b9d",
   "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,
   "id": "7b29d201",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "np.abs(π_seq - π_seq_f).max() < 1e-10"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b26c45d8",
   "metadata": {},
   "source": [
    "We thus conclude that  the likelihood ratio process is a key ingredient of the formula [(54.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 $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b54bcbd",
   "metadata": {},
   "source": [
    "## Behavior of  posterior probability $ \\{\\pi_t\\} $  under the subjective probability distribution\n",
    "\n",
    "We’ll end this lecture by briefly studying what our Baysian learner expects to learn under the\n",
    "subjective beliefs $ \\pi_t $ cranked out by Bayes’ law.\n",
    "\n",
    "This will provide us with some perspective  on our application of  Bayes’s law as a theory of learning.\n",
    "\n",
    "As we shall see, at each time $ t $, the Bayesian learner knows that he will be surprised.\n",
    "\n",
    "But he expects that new information will not lead him  to change his beliefs.\n",
    "\n",
    "And it won’t on average under his subjective beliefs.\n",
    "\n",
    "We’ll continue with our setting in which a McCall worker  knows that successive\n",
    "draws of his wage are drawn from either $ F $ or $ G $, but  does not know which of these two  distributions\n",
    "nature has drawn once-and-for-all before time $ 0 $.\n",
    "\n",
    "We’ll review and reiterate and rearrange some formulas that we have encountered above and in associated lectures.\n",
    "\n",
    "The worker’s initial beliefs induce a joint probability distribution\n",
    "over a potentially infinite sequence of draws $ w_0, w_1, \\ldots $.\n",
    "\n",
    "Bayes’ law is simply an application of  laws of\n",
    "probability to compute the conditional distribution of the $ t $th draw $ w_t $ conditional on $ [w_0, \\ldots, w_{t-1}] $.\n",
    "\n",
    "After our worker puts a subjective probability $ \\pi_{-1} $ on nature having selected distribution $ F $, we have in effect assumes from the start that the   decision maker **knows** the joint distribution  for the process $ \\{w_t\\}_{t=0} $.\n",
    "\n",
    "We assume that the worker also knows the laws of probability theory.\n",
    "\n",
    "A respectable view is that Bayes’ law is less a theory of learning than a statement  about the consequences of information inflows for a decision maker who thinks he knows the truth (i.e., a joint probability distribution) from the beginning."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54db5bad",
   "metadata": {},
   "source": [
    "### Mechanical details again\n",
    "\n",
    "At time $ 0 $ **before** drawing a wage offer, the worker attaches probability $ \\pi_{-1} \\in (0,1) $ to the distribution being $ F $.\n",
    "\n",
    "Before drawing a wage at time $ 0 $, the  worker thus believes that the density of $ w_0 $\n",
    "is\n",
    "\n",
    "$$\n",
    "h(w_0;\\pi_{-1}) = \\pi_{-1} f(w_0) + (1-\\pi_{-1}) g(w_0).\n",
    "$$\n",
    "\n",
    "Let $ a \\in \\{ f, g\\} $ be an index that indicates whether  nature chose permanently to draw from distribution $ f $ or from distribution $ g $.\n",
    "\n",
    "After drawing $ w_0 $, the worker uses Bayes’ law to deduce that\n",
    "the posterior  probability $ \\pi_0 = {\\rm Prob}{a = f | w_0} $\n",
    "that the density is $ f(w) $ is\n",
    "\n",
    "$$\n",
    "\\pi_0 = { \\pi_{-1} f(w_0) \\over \\pi_{-1} f(w_0) + (1-\\pi_{-1}) g(w_0)} .\n",
    "$$\n",
    "\n",
    "More generally,  after making the $ t $th draw and having   observed   $ w_t, w_{t-1}, \\ldots, w_0 $, the worker believes that\n",
    "the probability that $ w_{t+1} $ is  being drawn from  distribution  $ F $ is\n",
    "\n",
    "\n",
    "<a id='equation-eq-like44'></a>\n",
    "$$\n",
    "\\pi_t = \\pi_t(w_t | \\pi_{t-1}) \\equiv { \\pi_{t-1} f(w_t)/g(w_t) \\over \\pi_{t-1} f(w_t)/g(w_t) + (1-\\pi_{t-1})} \\tag{54.3}\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\n",
    "\\pi_t=\\frac{\\pi_{t-1} l_t(w_t)}{\\pi_{t-1} l_t(w_t)+1-\\pi_{t-1}}\n",
    "$$\n",
    "\n",
    "and that the density of $ w_{t+1} $ conditional on $ w_t, w_{t-1}, \\ldots, w_0 $ is\n",
    "\n",
    "$$\n",
    "h(w_{t+1};\\pi_{t}) = \\pi_{t} f(w_{t+1}) + (1-\\pi_{t}) g(w_{t+1}) .\n",
    "$$\n",
    "\n",
    "Notice that\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "E(\\pi_t | \\pi_{t-1}) & = \\int \\Bigl[  { \\pi_{t-1} f(w) \\over \\pi_{t-1} f(w) + (1-\\pi_{t-1})g(w)  } \\Bigr]\n",
    " \\Bigl[ \\pi_{t-1} f(w) + (1-\\pi_{t-1})g(w) \\Bigr]  d w \\cr\n",
    "& = \\pi_{t-1} \\int  f(w) dw  \\cr\n",
    "              & = \\pi_{t-1}, \\cr\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "so that the process $ \\pi_t $ is a **martingale**.\n",
    "\n",
    "Indeed, it is a **bounded martingale** because each $ \\pi_t $, being a probability,\n",
    "is between $ 0 $ and $ 1 $.\n",
    "\n",
    "In the first line in the above string of equalities, the term in the first set of brackets\n",
    "is just $ \\pi_t $ as a function of $ w_{t} $, while the term in the second set of brackets is the density of $ w_{t} $ conditional\n",
    "on $ w_{t-1}, \\ldots , w_0 $ or equivalently conditional on the *sufficient statistic* $ \\pi_{t-1} $ for $ w_{t-1}, \\ldots , w_0 $.\n",
    "\n",
    "Notice that here we are computing $ E(\\pi_t | \\pi_{t-1}) $ under the **subjective** density described in the second\n",
    "term in brackets.\n",
    "\n",
    "Because $ \\{\\pi_t\\} $ is a bounded martingale sequence, it follows from the **martingale convergence theorem** that $ \\pi_t $ converges almost surely to a random variable in $ [0,1] $.\n",
    "\n",
    "Practically, this means that  probability one is  attached to   sample paths\n",
    "$ \\{\\pi_t\\}_{t=0}^\\infty $ that  converge.\n",
    "\n",
    "According to the theorem,  it  different sample  paths  can converge to different limiting values.\n",
    "\n",
    "Thus, let $ \\{\\pi_t(\\omega)\\}_{t=0}^\\infty $ denote a particular sample path indexed by a particular $ \\omega\n",
    "\\in \\Omega $.\n",
    "\n",
    "We can think of nature as drawing an $ \\omega \\in \\Omega $ from a probability distribution\n",
    "$ {\\textrm{Prob}} \\Omega $ and then generating a single realization (or *simulation*) $ \\{\\pi_t(\\omega)\\}_{t=0}^\\infty $ of the process.\n",
    "\n",
    "The limit points of  $ \\{\\pi_t(\\omega)\\}_{t=0}^\\infty $ as $ t \\rightarrow +\\infty $ are realizations of a random variable that  is swept out as we sample $ \\omega $ from $ \\Omega $ and construct repeated draws of $ \\{\\pi_t(\\omega)\\}_{t=0}^\\infty $.\n",
    "\n",
    "By staring at law of motion [(54.1)](#equation-eq-recur1) or [(54.3)](#equation-eq-like44) , we can figure out some things about the probability distribution of the limit points\n",
    "\n",
    "$$\n",
    "\\pi_\\infty(\\omega) = \\lim_{t \\rightarrow + \\infty} \\pi_t(\\omega).\n",
    "$$\n",
    "\n",
    "Evidently, since the likelihood ratio $ \\ell(w_t) $ differs from $ 1 $ when $ f \\neq g $,\n",
    "as we have assumed, the only possible fixed points of [(54.3)](#equation-eq-like44) are\n",
    "\n",
    "$$\n",
    "\\pi_\\infty(\\omega) =1\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "\\pi_\\infty(\\omega) =0\n",
    "$$\n",
    "\n",
    "Thus, for some realizations, $ \\lim_{\\rightarrow + \\infty} \\pi_t(\\omega) =1 $\n",
    "while for other realizations,  $ \\lim_{\\rightarrow + \\infty} \\pi_t(\\omega) =0 $.\n",
    "\n",
    "Now let’s remember that $ \\{\\pi_t\\}_{t=0}^\\infty $ is a martingale and apply the law of iterated expectations.\n",
    "\n",
    "The law of iterated expectations implies\n",
    "\n",
    "$$\n",
    "E_t \\pi_{t+j}  = \\pi_t\n",
    "$$\n",
    "\n",
    "and in particular\n",
    "\n",
    "$$\n",
    "E_{-1} \\pi_{t+j} = \\pi_{-1}.\n",
    "$$\n",
    "\n",
    "Applying the above formula to $ \\pi_\\infty $, we obtain\n",
    "\n",
    "$$\n",
    "E_{-1} \\pi_\\infty(\\omega) = \\pi_{-1}\n",
    "$$\n",
    "\n",
    "where the mathematical expectation $ E_{-1} $ here is taken with respect to the probability\n",
    "measure $ {\\textrm{Prob}(\\Omega)} $.\n",
    "\n",
    "Since the only two values that $ \\pi_\\infty(\\omega) $ can take are $ 1 $ and $ 0 $, we know that for some $ \\lambda \\in [0,1] $\n",
    "\n",
    "$$\n",
    "{\\textrm{Prob}}\\Bigl(\\pi_\\infty(\\omega) = 1\\Bigr) = \\lambda, \\quad {\\textrm{Prob}}\\Bigl(\\pi_\\infty(\\omega) = 0\\Bigr) = 1- \\lambda\n",
    "$$\n",
    "\n",
    "and consequently that\n",
    "\n",
    "$$\n",
    "E_{-1} \\pi_\\infty(\\omega) = \\lambda \\cdot 1 + (1-\\lambda) \\cdot 0 = \\lambda\n",
    "$$\n",
    "\n",
    "Combining this equation with equation (20), we deduce that\n",
    "the probability that $ {\\textrm{Prob}(\\Omega)} $ attaches to\n",
    "$ \\pi_\\infty(\\omega) $ being $ 1 $ must be $ \\pi_{-1} $.\n",
    "\n",
    "Thus, under the worker’s subjective distribution, $ \\pi_{-1} $ of the sample paths\n",
    "of $ \\{\\pi_t\\} $ will converge pointwise to $ 1 $ and $ 1 - \\pi_{-1} $ of the sample paths will\n",
    "converge pointwise to $ 0 $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72ddb3ae",
   "metadata": {},
   "source": [
    "### Some simulations\n",
    "\n",
    "Let’s watch the martingale convergence theorem at work in some simulations of our learning model under the worker’s subjective distribution.\n",
    "\n",
    "Let us simulate $ \\left\\{ \\pi_{t}\\right\\} _{t=0}^{T} $, $ \\left\\{ w_{t}\\right\\} _{t=0}^{T} $ paths where for each $ t\\geq0 $, $ w_t $ is drawn from the subjective distribution\n",
    "\n",
    "$$\n",
    "\\pi_{t-1}f\\left(w_{t}\\right)+\\left(1-\\pi_{t-1}\\right)g\\left(w_{t}\\right)\n",
    "$$\n",
    "\n",
    "We’ll plot a large sample of paths."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8eac68af",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@jit\n",
    "def martingale_simulate(π0, N=5000, T=200):\n",
    "\n",
    "    π_path = np.empty((N,T+1))\n",
    "    w_path = np.empty((N,T))\n",
    "    π_path[:,0] = π0\n",
    "\n",
    "    for n in range(N):\n",
    "        π = π0\n",
    "        for t in range(T):\n",
    "            # draw w\n",
    "            if np.random.rand() <= π:\n",
    "                w = np.random.beta(F_a, F_b)\n",
    "            else:\n",
    "                w = np.random.beta(G_a, G_b)\n",
    "            π = π*f(w)/g(w)/(π*f(w)/g(w) + 1 - π)\n",
    "            π_path[n,t+1] = π\n",
    "            w_path[n,t] = w\n",
    "\n",
    "    return π_path, w_path\n",
    "\n",
    "def fraction_0_1(π0, N, T, decimals):\n",
    "\n",
    "    π_path, w_path = martingale_simulate(π0, N=N, T=T)\n",
    "    values, counts = np.unique(np.round(π_path[:,-1], decimals=decimals), return_counts=True)\n",
    "    return values, counts\n",
    "\n",
    "def create_table(π0s, N=10000, T=500, decimals=2):\n",
    "\n",
    "    outcomes = []\n",
    "    for π0 in π0s:\n",
    "        values, counts = fraction_0_1(π0, N=N, T=T, decimals=decimals)\n",
    "        freq = counts/N\n",
    "        outcomes.append(dict(zip(values, freq)))\n",
    "    table = pd.DataFrame(outcomes).sort_index(axis=1).fillna(0)\n",
    "    table.index = π0s\n",
    "    return table\n",
    "\n",
    "\n",
    "\n",
    "# simulate\n",
    "T = 200\n",
    "π0 = .5\n",
    "\n",
    "π_path, w_path = martingale_simulate(π0=π0, T=T, N=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d718dbf6",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "for i in range(100):\n",
    "    ax.plot(range(T+1), π_path[i, :])\n",
    "\n",
    "ax.set_xlabel('$t$')\n",
    "ax.set_ylabel('$\\pi_t$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "203c57a6",
   "metadata": {},
   "source": [
    "The above graph indicates that\n",
    "\n",
    "- each of paths converges  \n",
    "- some of the paths converge to $ 1 $  \n",
    "- some of the paths converge to $ 0 $  \n",
    "- none of the paths converge to a limit point not equal to $ 0 $ or $ 1 $  \n",
    "\n",
    "\n",
    "Convergence actually occurs pretty fast, as the following graph of the cross-ensemble distribution of $ \\pi_t $ for various small $ t $’s indicates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ee0176e",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "for t in [1, 10, T-1]:\n",
    "    ax.hist(π_path[:,t], bins=20, alpha=0.4, label=f'T={t}')\n",
    "\n",
    "ax.set_ylabel('count')\n",
    "ax.set_xlabel('$\\pi_T$')\n",
    "ax.legend(loc='lower right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63ffc614",
   "metadata": {},
   "source": [
    "Evidently, by $ t = 199 $, $ \\pi_t $ has converged to either $ 0 $ or $ 1 $.\n",
    "\n",
    "The fraction of paths that have converged to $ 1 $ is $ .5 $\n",
    "\n",
    "The fractions of paths that have converged to $ 0 $ is also $ .5 $.\n",
    "\n",
    "Does the fraction $ .5 $ ring a bell?\n",
    "\n",
    "Yes, it does: it equals the value of $ \\pi_0 = .5 $ that we used to generate each sequence\n",
    "in the ensemble.\n",
    "\n",
    "So let’s change $ \\pi_0 $ to $ .3 $ and watch what happens to the distribution of the ensemble of\n",
    "$ \\pi_t $’s for various $ t $’s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8508db56",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# simulate\n",
    "T = 200\n",
    "π0 = .3\n",
    "\n",
    "π_path3, w_path3 = martingale_simulate(π0=π0, T=T, N=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b334b57",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "for t in [1, 10, T-1]:\n",
    "    ax.hist(π_path3[:,t], bins=20, alpha=0.4, label=f'T={t}')\n",
    "\n",
    "ax.set_ylabel('count')\n",
    "ax.set_xlabel('$\\pi_T$')\n",
    "ax.legend(loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bea8d747",
   "metadata": {},
   "source": [
    "For the preceding ensemble that assumed $ \\pi_0 = .5 $, the following graph shows two  paths of\n",
    "$ w_t $’s and the $ \\pi_t $ sequences that gave rise to them.\n",
    "\n",
    "Notice that one of the paths involves systematically higher $ w_t $’s, outcomes that push $ \\pi_t $ upward.\n",
    "\n",
    "The luck of the draw early in a simulation push the subjective distribution to draw from\n",
    "$ F $ more frequently along a sample path, and this pushes $ \\pi_t $ toward $ 0 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93c2caa3",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "for i, j in enumerate([10, 100]):\n",
    "    ax.plot(range(T+1), π_path[j,:], color=colors[i], label=f'$\\pi$_path, {j}-th simulation')\n",
    "    ax.plot(range(1,T+1), w_path[j,:], color=colors[i], label=f'$w$_path, {j}-th simulation', alpha=0.3)\n",
    "\n",
    "ax.legend(loc='upper right')\n",
    "ax.set_xlabel('$t$')\n",
    "ax.set_ylabel('$\\pi_t$')\n",
    "ax2 = ax.twinx()\n",
    "ax2.set_ylabel(\"$w_t$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ebfb6c2",
   "metadata": {},
   "source": [
    "## Initial Prior is Verified by Paths Drawn from Subjective Conditional Densities\n",
    "\n",
    "Now let’s use our Python code to generate a table that checks out our earlier claims about the\n",
    "probability distribution of the pointwise limits $ \\pi_{\\infty}(\\omega) $.\n",
    "\n",
    "We’ll use our simulations to generate a histogram of this distribution.\n",
    "\n",
    "In the following table, the left column in bold face reports an assumed value of $ \\pi_{-1} $.\n",
    "\n",
    "The second column reports the fraction of $ N = 10000 $ simulations for which $ \\pi_{t} $  had converged to $ 0 $  at the terminal date $ T=500 $ for each simulation.\n",
    "\n",
    "The third column reports the fraction of $ N = 10000 $ simulations for which $ \\pi_{t} $  had converged to $ 1 $ as the terminal date $ T=500 $ for each simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eac1fb49",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# create table\n",
    "table = create_table(list(np.linspace(0,1,11)), N=10000, T=500)\n",
    "table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48bb05d2",
   "metadata": {},
   "source": [
    "The fraction of simulations for which $ \\pi_{t} $  had converged to $ 1 $ is indeed always  close  to $ \\pi_{-1} $, as anticipated."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23acb89b",
   "metadata": {},
   "source": [
    "## Drilling Down a Little Bit\n",
    "\n",
    "To understand how the local dynamics of $ \\pi_t $ behaves, it is enlightening to consult the  variance of $ \\pi_{t} $ conditional on $ \\pi_{t-1} $.\n",
    "\n",
    "Under the subjective distribution this conditional variance is defined as\n",
    "\n",
    "$$\n",
    "\\sigma^2(\\pi_t | \\pi_{t-1})  = \\int \\Bigl[  { \\pi_{t-1} f(w) \\over \\pi_{t-1} f(w) + (1-\\pi_{t-1})g(w)  } - \\pi_{t-1} \\Bigr]^2\n",
    " \\Bigl[ \\pi_{t-1} f(w) + (1-\\pi_{t-1})g(w) \\Bigr]  d w\n",
    "$$\n",
    "\n",
    "We can use  a Monte Carlo simulation to approximate this conditional variance.\n",
    "\n",
    "We approximate it for  a grid of points $ \\pi_{t-1} \\in [0,1] $.\n",
    "\n",
    "Then we’ll plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bb441e6",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@jit\n",
    "def compute_cond_var(pi, mc_size=int(1e6)):\n",
    "    # create monte carlo draws\n",
    "    mc_draws = np.zeros(mc_size)\n",
    "\n",
    "    for i in prange(mc_size):\n",
    "        if np.random.rand() <= pi:\n",
    "            mc_draws[i] = np.random.beta(F_a, F_b)\n",
    "        else:\n",
    "            mc_draws[i] = np.random.beta(G_a, G_b)\n",
    "\n",
    "    dev = pi*f(mc_draws)/(pi*f(mc_draws) + (1-pi)*g(mc_draws)) - pi\n",
    "    return np.mean(dev**2)\n",
    "\n",
    "pi_array = np.linspace(0, 1, 40)\n",
    "cond_var_array = []\n",
    "\n",
    "for pi in pi_array:\n",
    "    cond_var_array.append(compute_cond_var(pi))\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(pi_array, cond_var_array)\n",
    "ax.set_xlabel('$\\pi_{t-1}$')\n",
    "ax.set_ylabel('$\\sigma^{2}(\\pi_{t}\\\\vert \\pi_{t-1})$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c10303b",
   "metadata": {},
   "source": [
    "The shape of the the conditional variance as a function of $ \\pi_{t-1} $ is informative about the behavior of sample paths of $ \\{\\pi_t\\} $.\n",
    "\n",
    "Notice how the conditional variance approaches $ 0 $ for $ \\pi_{t-1} $ near  either $ 0 $ or $ 1 $.\n",
    "\n",
    "The conditional variance is nearly zero only when the agent  is almost sure that $ w_t $ is drawn from $ F $,  or is almost sure it is drawn from $ G $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c2ffd8e",
   "metadata": {},
   "source": [
    "## Sequels\n",
    "\n",
    "This lecture has been devoted to building some useful infrastructure that will help us understand inferences that are the foundations of\n",
    "results described  in [this lecture](https://python.quantecon.org/odu.html) and [this lecture](https://python.quantecon.org/wald_friedman.html) and [this lecture](https://python.quantecon.org/navy_captain.html)."
   ]
  }
 ],
 "metadata": {
  "date": 1735960782.2025597,
  "filename": "likelihood_bayes.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Likelihood Ratio Processes and Bayesian Learning"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}