prob_meaning.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a5934ae7",
   "metadata": {},
   "source": [
    "# Two Meanings of Probability"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "317a3ac2",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "This lecture  illustrates two distinct interpretations of a  **probability distribution**\n",
    "\n",
    "- A frequentist interpretation as **relative frequencies** anticipated to occur in a large i.i.d. sample  \n",
    "- A Bayesian interpretation as a **personal opinion** (about a parameter or list of parameters) after seeing a collection of observations  \n",
    "\n",
    "\n",
    "We recommend watching this video about **hypothesis testing** within  the frequentist approach\n",
    "\n",
    "After you watch that video, please watch the following video on the Bayesian approach to constructing **coverage intervals**\n",
    "\n",
    "After you are familiar with the material in these videos, this lecture uses the Socratic method to  to help consolidate your understanding of the different questions that are answered by\n",
    "\n",
    "- a frequentist confidence interval  \n",
    "- a Bayesian coverage interval  \n",
    "\n",
    "\n",
    "We do this  by inviting you to  write some  Python code.\n",
    "\n",
    "It would be especially useful if you tried doing this after each question that we pose for you,  before\n",
    "proceeding to read the rest of the lecture.\n",
    "\n",
    "We provide our own answers as the lecture unfolds, but you’ll learn more if you try writing your own code before reading and running ours.\n",
    "\n",
    "**Code for answering questions:**\n",
    "\n",
    "In addition to what’s in Anaconda, this lecture will deploy the following library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0d2740a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "pip install prettytable"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed67648a",
   "metadata": {},
   "source": [
    "To answer our coding questions, we’ll start with some imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d381103",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import prettytable as pt\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import binom\n",
    "import scipy.stats as st"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1d577c0",
   "metadata": {},
   "source": [
    "Empowered with these Python tools, we’ll now  explore the two meanings described above."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac54041f",
   "metadata": {},
   "source": [
    "## Frequentist Interpretation\n",
    "\n",
    "Consider the following classic example.\n",
    "\n",
    "The random variable  $ X $ takes on possible values $ k = 0, 1, 2, \\ldots, n $  with probabilties\n",
    "\n",
    "$$\n",
    "\\textrm{Prob}(X =  k | \\theta) =\n",
    "\\left(\\frac{n!}{k! (n-k)!} \\right) \\theta^k (1-\\theta)^{n-k}\n",
    "$$\n",
    "\n",
    "where the fixed parameter $ \\theta \\in (0,1) $.\n",
    "\n",
    "This is called   the **binomial distribution**.\n",
    "\n",
    "Here\n",
    "\n",
    "- $ \\theta $ is the probability that one toss of a coin will be a head, an outcome that we encode as  $ Y = 1 $.  \n",
    "- $ 1 -\\theta $ is the probability that one toss of the coin will be a tail, an outcome that we denote $ Y = 0 $.  \n",
    "- $ X $ is the total number of heads that came up after flipping the coin $ n $ times.  \n",
    "\n",
    "\n",
    "Consider the following experiment:\n",
    "\n",
    "Take $ I $ **independent** sequences of $ n $  **independent** flips of the coin\n",
    "\n",
    "Notice the repeated use of the adjective **independent**:\n",
    "\n",
    "- we use it once to describe that we are drawing $ n $ independent times from a **Bernoulli** distribution with parameter $ \\theta $ to arrive at one draw from a **Binomial** distribution with parameters\n",
    "  $ \\theta,n $.  \n",
    "- we use it again to describe that we are then drawing $ I $  sequences of $ n $ coin draws.  \n",
    "\n",
    "\n",
    "Let $ y_h^i \\in \\{0, 1\\} $ be the realized value of $ Y $ on the $ h $th flip during the $ i $th sequence of flips.\n",
    "\n",
    "Let $ \\sum_{h=1}^n y_h^i $ denote the total number of times  heads come up during the $ i $th sequence of $ n $ independent coin flips.\n",
    "\n",
    "Let $ f_k $ record the fraction of samples of length $ n $ for which $ \\sum_{h=1}^n y_h^i = k $:\n",
    "\n",
    "$$\n",
    "f_k^I = \\frac{\\textrm{number of samples of length n for which } \\sum_{h=1}^n y_h^i = k}{\n",
    "    I}\n",
    "$$\n",
    "\n",
    "The probability  $ \\textrm{Prob}(X =  k | \\theta) $ answers the following question:\n",
    "\n",
    "- As $ I $ becomes large, in what   fraction of  $ I $ independent  draws of  $ n $ coin flips should we anticipate  $ k $ heads to occur?  \n",
    "\n",
    "\n",
    "As usual, a law of large numbers justifies this answer."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bf0a709",
   "metadata": {},
   "source": [
    "## Exercise 10.1\n",
    "\n",
    "1. Please write a Python class to compute $ f_k^I $  \n",
    "1. Please use your code to compute $ f_k^I, k = 0, \\ldots , n $ and compare them to\n",
    "  $ \\textrm{Prob}(X =  k | \\theta) $ for various values of $ \\theta, n $ and $ I $  \n",
    "1. With the Law of Large numbers in mind, use your code to say something  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9067eb84",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 10.1](https://python.quantecon.org/#pm_ex1)\n",
    "\n",
    "Here is one solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13ef5198",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "class frequentist:\n",
    "\n",
    "    def __init__(self, θ, n, I):\n",
    "\n",
    "        '''\n",
    "        initialization\n",
    "        -----------------\n",
    "        parameters:\n",
    "        θ : probability that one toss of a coin will be a head with Y = 1\n",
    "        n : number of independent flips in each independent sequence of draws\n",
    "        I : number of independent sequence of draws\n",
    "\n",
    "        '''\n",
    "\n",
    "        self.θ, self.n, self.I = θ, n, I\n",
    "\n",
    "    def binomial(self, k):\n",
    "\n",
    "        '''compute the theoretical probability for specific input k'''\n",
    "\n",
    "        θ, n = self.θ, self.n\n",
    "        self.k = k\n",
    "        self.P = binom.pmf(k, n, θ)\n",
    "\n",
    "    def draw(self):\n",
    "\n",
    "        '''draw n independent flips for I independent sequences'''\n",
    "\n",
    "        θ, n, I = self.θ, self.n, self.I\n",
    "        sample = np.random.rand(I, n)\n",
    "        Y = (sample <= θ) * 1\n",
    "        self.Y = Y\n",
    "\n",
    "    def compute_fk(self, kk):\n",
    "\n",
    "        '''compute f_{k}^I for specific input k'''\n",
    "\n",
    "        Y, I = self.Y, self.I\n",
    "        K = np.sum(Y, 1)\n",
    "        f_kI = np.sum(K == kk) / I\n",
    "        self.f_kI = f_kI\n",
    "        self.kk = kk\n",
    "\n",
    "    def compare(self):\n",
    "\n",
    "        '''compute and print the comparison'''\n",
    "\n",
    "        n = self.n\n",
    "        comp = pt.PrettyTable()\n",
    "        comp.field_names = ['k', 'Theoretical', 'Frequentist']\n",
    "        self.draw()\n",
    "        for i in range(n):\n",
    "            self.binomial(i+1)\n",
    "            self.compute_fk(i+1)\n",
    "            comp.add_row([i+1, self.P, self.f_kI])\n",
    "        print(comp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21c52395",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "θ, n, k, I = 0.7, 20, 10, 1_000_000\n",
    "\n",
    "freq = frequentist(θ, n, I)\n",
    "\n",
    "freq.compare()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85281fa4",
   "metadata": {},
   "source": [
    "From the table above, can you see the law of large numbers at work?\n",
    "\n",
    "Let’s do some more calculations.\n",
    "\n",
    "**Comparison with different $ \\theta $**\n",
    "\n",
    "Now we fix\n",
    "\n",
    "$$\n",
    "n=20, k=10, I=1,000,000\n",
    "$$\n",
    "\n",
    "We’ll vary $ \\theta $ from $ 0.01 $ to $ 0.99 $ and plot outcomes against $ \\theta $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69020e24",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "θ_low, θ_high, npt = 0.01, 0.99, 50\n",
    "thetas = np.linspace(θ_low, θ_high, npt)\n",
    "P = []\n",
    "f_kI = []\n",
    "for i in range(npt):\n",
    "    freq = frequentist(thetas[i], n, I)\n",
    "    freq.binomial(k)\n",
    "    freq.draw()\n",
    "    freq.compute_fk(k)\n",
    "    P.append(freq.P)\n",
    "    f_kI.append(freq.f_kI)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74efe4bc",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "ax.grid()\n",
    "ax.plot(thetas, P, 'k-.', label='Theoretical')\n",
    "ax.plot(thetas, f_kI, 'r--', label='Fraction')\n",
    "plt.title(r'Comparison with different $\\theta$', fontsize=16)\n",
    "plt.xlabel(r'$\\theta$', fontsize=15)\n",
    "plt.ylabel('Fraction', fontsize=15)\n",
    "plt.tick_params(labelsize=13)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46cc621d",
   "metadata": {},
   "source": [
    "**Comparison with different $ n $**\n",
    "\n",
    "Now we fix $ \\theta=0.7, k=10, I=1,000,000 $ and vary $ n $ from $ 1 $ to $ 100 $.\n",
    "\n",
    "Then we’ll plot outcomes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0822ff09",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "n_low, n_high, nn = 1, 100, 50\n",
    "ns = np.linspace(n_low, n_high, nn, dtype='int')\n",
    "P = []\n",
    "f_kI = []\n",
    "for i in range(nn):\n",
    "    freq = frequentist(θ, ns[i], I)\n",
    "    freq.binomial(k)\n",
    "    freq.draw()\n",
    "    freq.compute_fk(k)\n",
    "    P.append(freq.P)\n",
    "    f_kI.append(freq.f_kI)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57cd9f6a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "ax.grid()\n",
    "ax.plot(ns, P, 'k-.', label='Theoretical')\n",
    "ax.plot(ns, f_kI, 'r--', label='Frequentist')\n",
    "plt.title(r'Comparison with different $n$', fontsize=16)\n",
    "plt.xlabel(r'$n$', fontsize=15)\n",
    "plt.ylabel('Fraction', fontsize=15)\n",
    "plt.tick_params(labelsize=13)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aba6ade2",
   "metadata": {},
   "source": [
    "**Comparison with different $ I $**\n",
    "\n",
    "Now we fix $ \\theta=0.7, n=20, k=10 $ and vary $ \\log(I) $ from $ 2 $ to $ 7 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8893adc5",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "I_log_low, I_log_high, nI = 2, 6, 200\n",
    "log_Is = np.linspace(I_log_low, I_log_high, nI)\n",
    "Is = np.power(10, log_Is).astype(int)\n",
    "P = []\n",
    "f_kI = []\n",
    "for i in range(nI):\n",
    "    freq = frequentist(θ, n, Is[i])\n",
    "    freq.binomial(k)\n",
    "    freq.draw()\n",
    "    freq.compute_fk(k)\n",
    "    P.append(freq.P)\n",
    "    f_kI.append(freq.f_kI)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbff1e33",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "ax.grid()\n",
    "ax.plot(Is, P, 'k-.', label='Theoretical')\n",
    "ax.plot(Is, f_kI, 'r--', label='Fraction')\n",
    "plt.title(r'Comparison with different $I$', fontsize=16)\n",
    "plt.xlabel(r'$I$', fontsize=15)\n",
    "plt.ylabel('Fraction', fontsize=15)\n",
    "plt.tick_params(labelsize=13)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "014a6178",
   "metadata": {},
   "source": [
    "From the above graphs, we can see that **$ I $, the number of independent sequences,** plays an important role.\n",
    "\n",
    "When $ I $ becomes larger, the difference between theoretical probability and frequentist estimate becomes smaller.\n",
    "\n",
    "Also, as long as $ I $ is large enough, changing $ \\theta $ or $ n $ does not substantially change the accuracy of the observed fraction\n",
    "as an approximation of $ \\theta $.\n",
    "\n",
    "The Law of Large Numbers is at work here.\n",
    "\n",
    "For each draw of an independent sequence, $ \\textrm{Prob}(X_i =  k | \\theta) $  is the same, so aggregating all draws forms an i.i.d sequence of a binary random variable $ \\rho_{k,i},i=1,2,...I $, with a mean of $ \\textrm{Prob}(X =  k | \\theta) $ and a variance of\n",
    "\n",
    "$$\n",
    "n \\cdot \\textrm{Prob}(X =  k | \\theta) \\cdot (1-\\textrm{Prob}(X =  k | \\theta)).\n",
    "$$\n",
    "\n",
    "So, by the LLN, the average of $ P_{k,i} $ converges to:\n",
    "\n",
    "$$\n",
    "E[\\rho_{k,i}] = \\textrm{Prob}(X =  k | \\theta) = \\left(\\frac{n!}{k! (n-k)!} \\right) \\theta^k (1-\\theta)^{n-k}\n",
    "$$\n",
    "\n",
    "as $ I $ goes to infinity."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40539184",
   "metadata": {},
   "source": [
    "## Bayesian Interpretation\n",
    "\n",
    "We again use a binomial distribution.\n",
    "\n",
    "But now we don’t regard  $ \\theta $ as being a fixed number.\n",
    "\n",
    "Instead, we think of it as a **random variable**.\n",
    "\n",
    "$ \\theta $ is described by a probability distribution.\n",
    "\n",
    "But now this probability distribution means something different than a relative frequency that we can anticipate to occur in a large i.i.d. sample.\n",
    "\n",
    "Instead, the probability distribution of $ \\theta $ is now a summary of our views about  likely values of $ \\theta $ either\n",
    "\n",
    "- **before** we have seen **any** data at all, or  \n",
    "- **before** we have seen **more** data, after we have seen **some** data  \n",
    "\n",
    "\n",
    "Thus, suppose that, before seeing any data, you have a personal prior probability distribution saying that\n",
    "\n",
    "$$\n",
    "P(\\theta) = \\frac{\\theta^{\\alpha-1}(1-\\theta)^{\\beta -1}}{B(\\alpha, \\beta)}\n",
    "$$\n",
    "\n",
    "where $ B(\\alpha, \\beta) $ is a  **beta function** , so that $ P(\\theta) $ is\n",
    "a **beta distribution** with parameters $ \\alpha, \\beta $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dc6dd8c",
   "metadata": {},
   "source": [
    "## Exercise 10.2\n",
    "\n",
    "**a)**  Please write down the **likelihood function** for a sample of length $ n $ from a binomial distribution with parameter $ \\theta $.\n",
    "\n",
    "**b)** Please write down the **posterior** distribution for $ \\theta $ after observing  one flip of the coin.\n",
    "\n",
    "**c)** Now pretend that the true value of $ \\theta = .4 $ and that someone who doesn’t know this has a beta prior distribution with parameters  with $ \\beta = \\alpha = .5 $. Please write a Python class to simulate this person’s personal posterior distribution for $ \\theta $  for a *single* sequence of $ n $ draws.\n",
    "\n",
    "**d)** Please plot the posterior distribution for $ \\theta $ as a function of $ \\theta $ as $ n $ grows as $ 1, 2, \\ldots $.\n",
    "\n",
    "**e)** For various $ n $’s, please describe and compute  a Bayesian coverage interval for the interval $ [.45, .55] $.\n",
    "\n",
    "**f)** Please tell what question a Bayesian coverage interval answers.\n",
    "\n",
    "**g)** Please compute the Posterior probabililty that $ \\theta \\in [.45, .55] $ for various values of sample size $ n $.\n",
    "\n",
    "**h)** Please use your Python class to study what happens to the posterior distribution as $ n \\rightarrow + \\infty $, again assuming that the true value of $ \\theta = .4 $, though it is unknown to the person doing the updating via Bayes’ Law."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d4378a5",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 10.2](https://python.quantecon.org/#pm_ex2)\n",
    "\n",
    "**a)** Please write down the **likelihood function** and the **posterior** distribution for $ \\theta $ after observing  one flip of our coin.\n",
    "\n",
    "Suppose the outcome is **Y**.\n",
    "\n",
    "The likelihood function is:\n",
    "\n",
    "$$\n",
    "L(Y|\\theta)= \\textrm{Prob}(X =  Y | \\theta) =\n",
    "\\theta^Y (1-\\theta)^{1-Y}\n",
    "$$\n",
    "\n",
    "**b)** Please write the **posterior** distribution for $ \\theta $ after observing  one flip of our coin.\n",
    "\n",
    "The prior distribution is\n",
    "\n",
    "$$\n",
    "\\textrm{Prob}(\\theta) = \\frac{\\theta^{\\alpha - 1} (1 - \\theta)^{\\beta - 1}}{B(\\alpha, \\beta)}\n",
    "$$\n",
    "\n",
    "We can derive the posterior distribution for $ \\theta $ via\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "  \\textrm{Prob}(\\theta | Y) &= \\frac{\\textrm{Prob}(Y | \\theta) \\textrm{Prob}(\\theta)}{\\textrm{Prob}(Y)} \\\\\n",
    "  &=\\frac{\\textrm{Prob}(Y | \\theta) \\textrm{Prob}(\\theta)}{\\int_{0}^{1} \\textrm{Prob}(Y | \\theta) \\textrm{Prob}(\\theta) d \\theta }\\\\\n",
    "  &= \\frac{\\theta^Y (1-\\theta)^{1-Y}\\frac{\\theta^{\\alpha - 1} (1 - \\theta)^{\\beta - 1}}{B(\\alpha, \\beta)}}{\\int_{0}^{1}\\theta^Y (1-\\theta)^{1-Y}\\frac{\\theta^{\\alpha - 1} (1 - \\theta)^{\\beta - 1}}{B(\\alpha, \\beta)} d \\theta } \\\\\n",
    "  &= \\frac{ \\theta^{Y+\\alpha - 1} (1 - \\theta)^{1-Y+\\beta - 1}}{\\int_{0}^{1}\\theta^{Y+\\alpha - 1} (1 - \\theta)^{1-Y+\\beta - 1} d \\theta}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "which means that\n",
    "\n",
    "$$\n",
    "\\textrm{Prob}(\\theta | Y) \\sim \\textrm{Beta}(\\alpha + Y, \\beta + (1-Y))\n",
    "$$\n",
    "\n",
    "Now please pretend that the true value of $ \\theta = .4 $ and that someone who doesn’t know this has a beta prior with $ \\beta = \\alpha = .5 $.\n",
    "\n",
    "**c)** Now pretend that the true value of $ \\theta = .4 $ and that someone who doesn’t know this has a beta prior distribution with parameters  with $ \\beta = \\alpha = .5 $. Please write a Python class to simulate this person’s personal posterior distribution for $ \\theta $  for a *single* sequence of $ n $ draws."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6edb8fc",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "class Bayesian:\n",
    "\n",
    "    def __init__(self, θ=0.4, n=1_000_000, α=0.5, β=0.5):\n",
    "        \"\"\"\n",
    "        Parameters:\n",
    "        ----------\n",
    "        θ : float, ranging from [0,1].\n",
    "           probability that one toss of a coin will be a head with Y = 1\n",
    "\n",
    "        n : int.\n",
    "           number of independent flips in an independent sequence of draws\n",
    "\n",
    "        α&β : int or float.\n",
    "             parameters of the prior distribution on θ\n",
    "\n",
    "        \"\"\"\n",
    "        self.θ, self.n, self.α, self.β = θ, n, α, β\n",
    "        self.prior = st.beta(α, β)\n",
    "\n",
    "    def draw(self):\n",
    "        \"\"\"\n",
    "        simulate a single sequence of draws of length n, given probability θ\n",
    "\n",
    "        \"\"\"\n",
    "        array = np.random.rand(self.n)\n",
    "        self.draws = (array < self.θ).astype(int)\n",
    "\n",
    "    def form_single_posterior(self, step_num):\n",
    "        \"\"\"\n",
    "        form a posterior distribution after observing the first step_num elements of the draws\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        step_num: int.\n",
    "               number of steps observed to form a posterior distribution\n",
    "\n",
    "        Returns\n",
    "        ------\n",
    "        the posterior distribution for sake of plotting in the subsequent steps\n",
    "\n",
    "        \"\"\"\n",
    "        heads_num = self.draws[:step_num].sum()\n",
    "        tails_num = step_num - heads_num\n",
    "\n",
    "        return st.beta(self.α+heads_num, self.β+tails_num)\n",
    "\n",
    "    def form_posterior_series(self,num_obs_list):\n",
    "        \"\"\"\n",
    "        form a series of posterior distributions that form after observing different number of draws.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        num_obs_list: a list of int.\n",
    "               a list of the number of observations used to form a series of posterior distributions.\n",
    "\n",
    "        \"\"\"\n",
    "        self.posterior_list = []\n",
    "        for num in num_obs_list:\n",
    "            self.posterior_list.append(self.form_single_posterior(num))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d4f29fd",
   "metadata": {},
   "source": [
    "**d)** Please plot the posterior distribution for $ \\theta $ as a function of $ \\theta $ as $ n $ grows from $ 1, 2, \\ldots $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f99f562",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "Bay_stat = Bayesian()\n",
    "Bay_stat.draw()\n",
    "\n",
    "num_list = [1, 2, 3, 4, 5, 10, 20, 30, 50, 70, 100, 300, 500, 1000, # this line for finite n\n",
    "            5000, 10_000, 50_000, 100_000, 200_000, 300_000]  # this line for approximately infinite n\n",
    "\n",
    "Bay_stat.form_posterior_series(num_list)\n",
    "\n",
    "θ_values = np.linspace(0.01, 1, 100)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 6))\n",
    "\n",
    "ax.plot(θ_values, Bay_stat.prior.pdf(θ_values), label='Prior Distribution', color='k', linestyle='--')\n",
    "\n",
    "for ii, num in enumerate(num_list[:14]):\n",
    "    ax.plot(θ_values, Bay_stat.posterior_list[ii].pdf(θ_values), label='Posterior with n = %d' % num)\n",
    "\n",
    "ax.set_title('P.D.F of Posterior Distributions', fontsize=15)\n",
    "ax.set_xlabel(r\"$\\theta$\", fontsize=15)\n",
    "\n",
    "ax.legend(fontsize=11)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fed502c",
   "metadata": {},
   "source": [
    "**e)** For various $ n $’s, please describe and compute  $ .05 $ and $ .95 $ quantiles for  posterior probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f638760a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "upper_bound = [ii.ppf(0.05) for ii in Bay_stat.posterior_list[:14]]\n",
    "lower_bound = [ii.ppf(0.95) for ii in Bay_stat.posterior_list[:14]]\n",
    "\n",
    "interval_df = pd.DataFrame()\n",
    "interval_df['upper'] = upper_bound\n",
    "interval_df['lower'] = lower_bound\n",
    "interval_df.index = num_list[:14]\n",
    "interval_df = interval_df.T\n",
    "interval_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d0dfa9f",
   "metadata": {},
   "source": [
    "As $ n $ increases, we can see that Bayesian coverage intervals narrow and move toward $ 0.4 $.\n",
    "\n",
    "**f)** Please tell what question a Bayesian coverage interval answers.\n",
    "\n",
    "The Bayesian coverage interval tells the range of $ \\theta $ that corresponds to the [$ p_1 $, $ p_2 $] quantiles of the cumulative probability distribution (CDF)  of the posterior distribution.\n",
    "\n",
    "To construct the coverage interval we first compute a posterior distribution of the unknown parameter $ \\theta $.\n",
    "\n",
    "If the CDF is $ F(\\theta) $, then the Bayesian coverage interval $ [a,b] $ for the interval $ [p_1,p_2] $ is described by\n",
    "\n",
    "$$\n",
    "F(a)=p_1,F(b)=p_2\n",
    "$$\n",
    "\n",
    "**g)** Please compute the Posterior probabililty that $ \\theta \\in [.45, .55] $ for various values of sample size $ n $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa34b444",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "left_value, right_value = 0.45, 0.55\n",
    "\n",
    "posterior_prob_list=[ii.cdf(right_value)-ii.cdf(left_value) for ii in Bay_stat.posterior_list]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "ax.plot(posterior_prob_list)\n",
    "ax.set_title('Posterior Probabililty that '+ r\"$\\theta$\" +' Ranges from %.2f to %.2f'%(left_value, right_value),\n",
    "             fontsize=13)\n",
    "ax.set_xticks(np.arange(0, len(posterior_prob_list), 3))\n",
    "ax.set_xticklabels(num_list[::3])\n",
    "ax.set_xlabel('Number of Observations', fontsize=11)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45e2bdeb",
   "metadata": {},
   "source": [
    "Notice that in the graph above the posterior probabililty that $ \\theta \\in [.45, .55] $ typically exhibits a hump shape as $ n $ increases.\n",
    "\n",
    "Two opposing forces are at work.\n",
    "\n",
    "The first force is that the individual  adjusts his belief as he observes new outcomes, so his posterior probability distribution  becomes more and more realistic, which explains the rise of the posterior probabililty.\n",
    "\n",
    "However, $ [.45, .55] $ actually excludes the true $ \\theta =.4 $ that generates the data.\n",
    "\n",
    "As a result, the posterior probabililty drops as larger and larger samples refine his  posterior probability distribution of $ \\theta $.\n",
    "\n",
    "The descent seems precipitous only because of the scale of the graph  that has the number of observations increasing disproportionately.\n",
    "\n",
    "When the number of observations becomes large enough, our Bayesian becomes so confident about $ \\theta $ that he considers $ \\theta \\in [.45, .55] $ very unlikely.\n",
    "\n",
    "That is why we see a nearly horizontal line when the number of observations exceeds 500.\n",
    "\n",
    "**h)** Please use your Python class to study what happens to the posterior distribution as $ n \\rightarrow + \\infty $, again assuming that the true value of $ \\theta = .4 $, though it is unknown to the person doing the updating via Bayes’ Law.\n",
    "\n",
    "Using the Python class we made above, we can see the evolution of posterior distributions as $ n $ approaches infinity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d139b15e",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 6))\n",
    "\n",
    "for ii, num in enumerate(num_list[14:]):\n",
    "    ii += 14\n",
    "    ax.plot(θ_values, Bay_stat.posterior_list[ii].pdf(θ_values),\n",
    "            label='Posterior with n=%d thousand' % (num/1000))\n",
    "\n",
    "ax.set_title('P.D.F of Posterior Distributions', fontsize=15)\n",
    "ax.set_xlabel(r\"$\\theta$\", fontsize=15)\n",
    "ax.set_xlim(0.3, 0.5)\n",
    "\n",
    "ax.legend(fontsize=11)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14207d1c",
   "metadata": {},
   "source": [
    "As $ n $ increases, we can see that the probability density functions *concentrate* on $ 0.4 $, the true value of $ \\theta $.\n",
    "\n",
    "Here the  posterior means  converges to $ 0.4 $ while the posterior standard deviations converges to $ 0 $ from above.\n",
    "\n",
    "To show this, we compute the means and variances statistics of the posterior distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90d9117e",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "mean_list = [ii.mean() for ii in Bay_stat.posterior_list]\n",
    "std_list = [ii.std() for ii in Bay_stat.posterior_list]\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, figsize=(14, 5))\n",
    "\n",
    "ax[0].plot(mean_list)\n",
    "ax[0].set_title('Mean Values of Posterior Distribution', fontsize=13)\n",
    "ax[0].set_xticks(np.arange(0, len(mean_list), 3))\n",
    "ax[0].set_xticklabels(num_list[::3])\n",
    "ax[0].set_xlabel('Number of Observations', fontsize=11)\n",
    "\n",
    "ax[1].plot(std_list)\n",
    "ax[1].set_title('Standard Deviations of Posterior Distribution', fontsize=13)\n",
    "ax[1].set_xticks(np.arange(0, len(std_list), 3))\n",
    "ax[1].set_xticklabels(num_list[::3])\n",
    "ax[1].set_xlabel('Number of Observations', fontsize=11)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90a9960e",
   "metadata": {},
   "source": [
    "How shall we interpret the patterns above?\n",
    "\n",
    "The answer is encoded in the  Bayesian updating formulas.\n",
    "\n",
    "It is natural to extend the one-step Bayesian update to an $ n $-step Bayesian update.\n",
    "\n",
    "$$\n",
    "\\textrm{Prob}(\\theta|k) = \\frac{\\textrm{Prob}(\\theta,k)}{\\textrm{Prob}(k)}=\\frac{\\textrm{Prob}(k|\\theta)*\\textrm{Prob}(\\theta)}{\\textrm{Prob}(k)}=\\frac{\\textrm{Prob}(k|\\theta)*\\textrm{Prob}(\\theta)}{\\int_0^1 \\textrm{Prob}(k|\\theta)*\\textrm{Prob}(\\theta) d\\theta}\n",
    "$$\n",
    "\n",
    "$$\n",
    "=\\frac{{N \\choose k} (1 - \\theta)^{N-k} \\theta^k*\\frac{\\theta^{\\alpha - 1} (1 - \\theta)^{\\beta - 1}}{B(\\alpha, \\beta)}}{\\int_0^1 {N \\choose k} (1 - \\theta)^{N-k} \\theta^k*\\frac{\\theta^{\\alpha - 1} (1 - \\theta)^{\\beta - 1}}{B(\\alpha, \\beta)} d\\theta}\n",
    "$$\n",
    "\n",
    "$$\n",
    "=\\frac{(1 -\\theta)^{\\beta+N-k-1}* \\theta^{\\alpha+k-1}}{\\int_0^1 (1 - \\theta)^{\\beta+N-k-1}* \\theta^{\\alpha+k-1} d\\theta}\n",
    "$$\n",
    "\n",
    "$$\n",
    "={Beta}(\\alpha + k, \\beta+N-k)\n",
    "$$\n",
    "\n",
    "A beta distribution with $ \\alpha $ and $ \\beta $ has the following mean and variance.\n",
    "\n",
    "The mean is $ \\frac{\\alpha}{\\alpha + \\beta} $\n",
    "\n",
    "The variance is $ \\frac{\\alpha \\beta}{(\\alpha + \\beta)^2 (\\alpha + \\beta + 1)} $\n",
    "\n",
    "- $ \\alpha $ can be viewed as the number of successes  \n",
    "- $ \\beta $ can be viewed as the number of failures  \n",
    "\n",
    "\n",
    "The random variables $ k $ and $ N-k $ are governed by Binomial Distribution with $ \\theta=0.4 $.\n",
    "\n",
    "Call this the true data generating process.\n",
    "\n",
    "According to the Law of Large Numbers, for a large number of observations, observed frequencies of $ k $ and $ N-k $ will be described by the true data generating process, i.e., the population probability distribution that we assumed when generating the observations on the computer. (See Exercise 10.1).\n",
    "\n",
    "Consequently, the  mean of the posterior distribution converges to $ 0.4 $ and the variance withers to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5af98365",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "upper_bound = [ii.ppf(0.95) for ii in Bay_stat.posterior_list]\n",
    "lower_bound = [ii.ppf(0.05) for ii in Bay_stat.posterior_list]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 6))\n",
    "ax.scatter(np.arange(len(upper_bound)), upper_bound, label='95 th Quantile')\n",
    "ax.scatter(np.arange(len(lower_bound)), lower_bound, label='05 th Quantile')\n",
    "\n",
    "ax.set_xticks(np.arange(0, len(upper_bound), 2))\n",
    "ax.set_xticklabels(num_list[::2])\n",
    "ax.set_xlabel('Number of Observations', fontsize=12)\n",
    "ax.set_title('Bayesian Coverage Intervals of Posterior Distributions', fontsize=15)\n",
    "\n",
    "ax.legend(fontsize=11)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aec07683",
   "metadata": {},
   "source": [
    "After observing a large number of outcomes, the  posterior distribution collapses around $ 0.4 $.\n",
    "\n",
    "Thus, the Bayesian statististian  comes to believe that $ \\theta $ is near $ .4 $.\n",
    "\n",
    "As shown in the figure above, as the number of observations grows, the Bayesian coverage intervals (BCIs) become narrower and narrower   around  $ 0.4 $.\n",
    "\n",
    "However, if you take a closer look, you will find that the centers of  the BCIs are not exactly $ 0.4 $, due to the persistent influence of the prior distribution and the randomness of the simulation path."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e15d09fe",
   "metadata": {},
   "source": [
    "## Role of a Conjugate Prior\n",
    "\n",
    "We have made  assumptions that link functional forms of  our likelihood function and our prior in a way that has eased our calculations considerably.\n",
    "\n",
    "In particular, our assumptions that the likelihood function is **binomial** and that the prior distribution is a **beta distribution** have the consequence that the posterior distribution implied by Bayes’ Law is also a **beta distribution**.\n",
    "\n",
    "So posterior and prior are both beta distributions, albeit ones with different parameters.\n",
    "\n",
    "When a likelihood function and prior fit together like hand and glove in this way, we can  say that the  prior and posterior are **conjugate distributions**.\n",
    "\n",
    "In this situation, we also sometimes  say that we have **conjugate prior** for the likelihood function $ \\textrm{Prob}(X | \\theta) $.\n",
    "\n",
    "Typically, the functional form of the likelihood function determines the functional form of a **conjugate prior**.\n",
    "\n",
    "A natural question to ask is why should a person’s personal prior about a parameter $ \\theta $ be restricted to be described by a conjugate prior?\n",
    "\n",
    "Why not some other functional form that more sincerely describes the person’s beliefs?\n",
    "\n",
    "To be argumentative, one could ask, why should the form of the likelihood function have *anything* to say about my personal beliefs about $ \\theta $?\n",
    "\n",
    "A dignified response to that question is, well, it shouldn’t, but if you want to compute a posterior easily you’ll just be happier if your prior is conjugate to your likelihood.\n",
    "\n",
    "Otherwise, your posterior won’t have a convenient analytical form and you’ll be in the situation of wanting to apply the Markov chain Monte Carlo techniques deployed in [this quantecon lecture](https://python.quantecon.org/bayes_nonconj.html).\n",
    "\n",
    "We also apply these powerful methods to approximating Bayesian posteriors for non-conjugate priors in\n",
    "[this quantecon lecture](https://python.quantecon.org/ar1_bayes.html) and [this quantecon lecture](https://python.quantecon.org/ar1_turningpts.html)"
   ]
  }
 ],
 "metadata": {
  "date": 1733369122.9903927,
  "filename": "prob_meaning.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Two Meanings of Probability"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}