Chapter 10 - Parametric Inference.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 10. Parametric Inference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Parametric models** are of the form\n",
    "\n",
    "$$ \\mathfrak{F} = \\bigg\\{ f(x; \\theta) : \\; \\theta \\in \\Theta \\bigg\\} $$\n",
    "\n",
    "where $\\Theta \\subset \\mathbb{R}^k$ is the parameter space and $\\theta = (\\theta_1, \\dots, \\theta_k)$ is the parameter.  The problem of inference then reduces to the problem of estimating parameter $\\theta$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10.1 Parameter of interest\n",
    "\n",
    "Often we are only interested in some function $T(\\theta)$.  For example, if $X \\sim N(\\mu, \\sigma^2)$ then the parameter is $\\theta = (\\mu, \\sigma)$.  If our goal is to estimate $\\mu$ then $\\mu = T(\\theta)$ is called the **parameter of interest** and $\\sigma$ is called  a **nuisance parameter**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10.2 The Method of Moments\n",
    "\n",
    "Suppose that the parameter $\\theta = (\\theta_1, \\dots, \\theta_n)$ has $k$ components.  For $1 \\leq j \\leq k$ define the $j$-th **moment**\n",
    "\n",
    "$$ \\alpha_j \\equiv \\alpha_j(\\theta) = \\mathbb{E}_\\theta(X^j) = \\int x^j dF_\\theta(x)$$\n",
    "\n",
    "and the $j$-th **sample moment**\n",
    "\n",
    "$$ \\hat{\\alpha}_j = \\frac{1}{n} \\sum_{i=1}^n X_i^j $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **method of moments estimator** $\\hat{\\theta}_n$ is defined to be the value of $\\theta$ such that\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\alpha_1(\\hat{\\theta}_n) &= \\hat{\\alpha_1} \\\\\n",
    "\\alpha_2(\\hat{\\theta}_n) &= \\hat{\\alpha_2} \\\\\n",
    "\\vdots \\\\\n",
    "\\alpha_k(\\hat{\\theta}_n) &= \\hat{\\alpha_k}              \n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "This defines a system of $k$ equations with $k$ unknowns."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 10.6**.  Let $\\hat{\\theta}_n$ denote the method of moments estimator.  Under the conditions given in the appendix, the following statements hold:\n",
    "\n",
    "(1) The estimate $\\hat{\\theta}_n$ exists with probability tending to 1.\n",
    "\n",
    "(2) The estimate is consistent: $\\hat{\\theta}_n \\xrightarrow{\\text{P}} \\theta$.\n",
    "\n",
    "(3) The estimate is asymptotically Normal:\n",
    "\n",
    "$$\\sqrt(n)(\\hat{\\theta}_n - \\theta) \\leadsto N(0, \\Sigma) $$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\Sigma = g \\mathbb{E}_\\theta (Y Y^T) g^T \\\\\n",
    "Y = (X, X^2, \\dots, X^k)^T, \\quad g = (g_1, \\dots, g_k) \\quad \\text{and} \\quad g_j = \\partial \\alpha_j^{-1}(\\theta)/\\partial\\theta\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last statement in Theorem 10.6 can be used to find standard errors and confidence intervals.  However, there is an easier way: the bootstrap."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10.3 Maximum Likelihood"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $X_1, \\dots, X_n$ be iid with PDF $f(x; \\theta)$.\n",
    "\n",
    "The **likelihood function** is defined by \n",
    "\n",
    "$$ \\mathcal{L}_n(\\theta) = \\prod_{i=1}^n f(X_i; \\theta) $$\n",
    "\n",
    "The **log-likelihood function** is defined by $\\ell_n(\\theta) = \\log \\mathcal{L}_n(\\theta)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The likelihood function is just the joint density of the data, except we treat is as a function of parameter $\\theta$.  Thus $\\mathcal{L}_n : \\Theta \\rightarrow [0, \\infty)$.  The likelihood function is not a density function; in general it is not true that $\\mathcal{L}_n$ integrates to 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **maximum likelihood estimator** MLE, denoted by $\\hat{\\theta}_n$, is the value of $\\theta$ that maximizes $\\mathcal{L}_n(\\theta)$.\n",
    "\n",
    "The maximum of $\\ell_n(\\theta)$ occurs at the same place as the maximum of $\\mathcal{L}_n(\\theta)$, so maximizing either leads to the same answer.  Often it's easier to maximize the log-likelihood."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10.4 Properties of Maximum Likelihood Estimators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Under certain conditions on the model, the MLE $\\hat{\\theta}_n$ possesses many properties that make it an appealing choice of estimator.\n",
    "\n",
    "The main properties of the MLE are:\n",
    "\n",
    "- It is **consistent**: $\\hat{\\theta}_n \\xrightarrow{\\text{P}} \\theta_*$, where $\\theta_*$ denotes the true value of parameter $\\theta$.\n",
    "- It is **equivariant**: if $\\hat{\\theta}_n$ is the MLE of $\\theta$ then $g(\\hat{\\theta}_n)$ is the MLE of $g(\\theta)$.\n",
    "- If is **asymptotically Normal**: $\\sqrt{n}(\\hat{\\theta} - \\theta_*) / \\hat{\\text{se}} \\leadsto N(0, 1)$ where $\\hat{\\text{se}}$ can be computed analytically.\n",
    "- It is **asymptotically optimal** or **efficient**: roughly, this means that among all well behaved estimators, the MLE has the smallest variance, at least for large samples.\n",
    "- The MLE is approximately the Bayes estimator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10.5 Consistency of Maximum Likelihood Estimator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $f$ and $g$ are PDFs, define the **Kullback-Leibler distance** between $f$ and $g$ to be:\n",
    "\n",
    "$$ D(f, g) = \\int f(x) \\log \\left( \\frac{f(x)}{g(x)} \\right) dx $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be shown that $D(f, g) \\geq 0$ and $D(f, f) = 0$.  For any $\\theta, \\psi \\in \\Theta$ write $D(\\theta, \\psi)$ to mean $D(f(x; \\theta), f(x; \\psi))$.  We will assume that $\\theta \\neq \\psi$ implies $D(\\theta, \\psi) > 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $\\theta_*$ denote the true value of $\\theta$.  Maximizing $\\ell_n(\\theta)$ is equivalent to maximizing\n",
    "\n",
    "$$M_n(\\theta) = \\frac{1}{n} \\sum_i \\log \\frac{f(X_i; \\theta)}{f(X_i; \\theta_*)}$$\n",
    "\n",
    "By the law of large numbers, $M_n(\\theta)$ converges to:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{E}_{\\theta_*} \\left( \\log \\frac{f(X_i; \\theta)}{f(X_i; \\theta_*)} \\right)\n",
    "& = \\int \\log \\left( \\frac{f(x; \\theta)}{f(x; \\theta_*)} \\right) f(x; \\theta_*) dx \\\\\n",
    "& = - \\int \\log \\left( \\frac{f(x; \\theta_*)}{f(x; \\theta)} \\right) f(x; \\theta_*) dx \\\\\n",
    "&= -D(\\theta_*, \\theta)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Hence $M_n(\\theta) \\approx -D(\\theta_*, \\theta)$ which is maximized at $\\theta_*$, since the KL distance is 0 when $\\theta_* = \\theta$ and positive otherwise.  Hence, we expect that the maximizer will tend to $\\theta_*$.\n",
    "\n",
    "To prove this formally, we need more than $M_n(\\theta) \\xrightarrow{\\text{P}} -D(\\theta_*, \\theta)$.  We need this convergence to be uniform over $\\theta$.  We also have to make sure that the KL distance is well-behaved.  Here are the formal details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 10.13**.  Let $\\theta_*$ denote the true value of $\\theta$.  Define\n",
    "\n",
    "$$M_n(\\theta) = \\frac{1}{n} \\sum_i \\log \\frac{f(X_i; \\theta)}{f(X_i; \\theta_*)}$$\n",
    "\n",
    "and $M(\\theta) = -D(\\theta_*, \\theta)$.  Suppose that\n",
    "\n",
    "$$ \\sup _{\\theta \\in \\Theta} |M_n(\\theta) - M(\\theta)| \\xrightarrow{\\text{P}} 0 $$\n",
    "\n",
    "and that, for every $\\epsilon > 0$,\n",
    "\n",
    "$$ \\sup _{\\theta : |\\theta - \\theta_*| \\geq \\epsilon} M(\\theta) < M(\\theta_*)$$\n",
    "\n",
    "Let $\\hat{\\theta}_n$ denote the mle.  Then $\\hat{\\theta}_n \\xrightarrow{\\text{P}} \\theta_*$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10.6 Equivalence of the MLE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 10.14**.  Let $\\tau = g(\\theta)$ be a one-to-one function of $\\theta$.  Let $\\hat{\\theta}_n$ be the MLE of $\\theta$.  Then $\\hat{\\tau}_n = g(\\hat{\\theta}_n)$ is the MLE of $\\tau$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proof**.  Let $h = g^{-1}$ denote the inverse of $g$.  Then $\\hat{\\theta}_n = h(\\hat{\\tau}_n)$.  For any $\\tau$, $L(\\tau) = \\prod_i f(x_i; h(\\tau)) = \\prod_i f(x_i; \\theta) = \\mathcal{L}(\\theta)$ where $\\theta = h(\\tau)$.  Hence, for any $\\tau$, $\\mathcal{L}_n(\\tau) = \\mathcal{L}(\\theta) \\leq \\mathcal{L}(\\hat{\\theta}) = \\mathcal{L}_n(\\hat{\\tau})$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10.7 Asymptotic Normality"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **score function** is defined to be\n",
    "\n",
    "$$ s(X; \\theta) = \\frac{\\partial \\log f(X; \\theta)}{\\partial \\theta} $$\n",
    "\n",
    "The **Fisher information** is defined to be\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "I_n(\\theta) &= \\mathbb{V}_\\theta \\left( \\sum_{i=1}^n s(X_i; \\theta) \\right) \\\\\n",
    "&= \\sum_{i=1}^n \\mathbb{V}_\\theta(s(X_i; \\theta))\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For $n = 1$ we sometimes write $I(\\theta)$ instead of $I_1(\\theta)$.\n",
    "\n",
    "It can be shown that $\\mathbb{E}_\\theta(s(X; \\theta)) = 0$.  It then follows that $\\mathbb{V}_\\theta(s(X; \\theta)) = \\mathbb{E}_\\theta((s(X; \\theta))^2)$.  A further simplification of $I_n(\\theta)$ is given in the next result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 10.17**.\n",
    "\n",
    "$$ I_n(\\theta) = n I(\\theta)$$\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "I(\\theta) & = -\\mathbb{E}_\\theta \\left( \\frac{\\partial^2 \\log f(X; \\theta)}{\\partial^2 \\theta^2} \\right) \\\\\n",
    "&= - \\int \\left( \\frac{\\partial^2 \\log f(x; \\theta)}{\\partial^2 \\theta^2} \\right) f(x; \\theta) dx\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 10.18 (Asymptotic Normality of the MLE)**.  Under appropriate regularity conditions, the following hold:\n",
    "\n",
    "(1)  Let $\\text{se} = \\sqrt{1 / I_n(\\theta)}$.  Then,\n",
    "\n",
    "$$ \\frac{\\hat{\\theta}_n - \\theta}{\\text{se}} \\leadsto N(0, 1) $$\n",
    "\n",
    "(2) Let $\\hat{\\text{se}} = \\sqrt{1 / I_n(\\hat{\\theta}_n)}$.  Then,\n",
    "\n",
    "$$ \\frac{\\hat{\\theta}_n - \\theta}{\\hat{\\text{se}}} \\leadsto N(0, 1) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first statement says that $\\hat{\\theta}_n \\approx N(\\theta, \\text{se})$.  The second statement says that this is still true if we replace the standard error $\\text{se}$ by its estimated standard error $\\hat{\\text{se}}$.\n",
    "\n",
    "Informally this says that the distribution of the MLE can be approximated with $N(\\theta, \\hat{\\text{se}})$.  From this fact we can construct an asymptotic confidence interval."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 10.19**.  Let\n",
    "\n",
    "$$ C_n = \\left( \\hat{\\theta_n} - z_{\\alpha/2} \\hat{\\text{se}}, \\; \\hat{\\theta_n} + z_{\\alpha/2} \\hat{\\text{se}} \\right) $$\n",
    "\n",
    "Then, $\\mathbb{P}_\\theta(\\theta \\in C_n) \\rightarrow 1 - \\alpha$ as $n \\rightarrow \\infty$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proof**  Let $Z$ denote a standard random variable.  Then,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{P}_\\theta(\\theta \\in C_n) \n",
    "&= \\mathbb{P}_\\theta(\\hat{\\theta}_n - z_{\\alpha/2} \\hat{\\text{se}} \\leq \\theta \\leq \\hat{\\theta}_n + z_{\\alpha/2} \\hat{\\text{se}}) \\\\\n",
    "&= \\mathbb{P}_\\theta(-z_{\\alpha/2} \\leq \\frac{\\hat{\\theta}_n - \\theta}{\\hat{\\text{se}}} \\leq z_{\\alpha/2}) \\\\\n",
    "&\\rightarrow \\mathbb{P}(-z_{\\alpha/2} \\leq Z \\leq z_{\\alpha/2}) = 1 - \\alpha\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10.8 Optimality"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose that $X_1, \\dots, X_n \\sim N(0, \\sigma^2)$.  The MLE is $\\hat{\\theta}_n = \\overline{X}_n$.  Another reasonable estimator is the sample median $\\overline{\\theta}_n$.  The MLE satisfies\n",
    "\n",
    "$$ \\sqrt{n}(\\hat{\\theta}_n - \\theta) \\leadsto N(0, \\sigma^2) $$\n",
    "\n",
    "It can be proved that the median satisfies\n",
    "\n",
    "$$ \\sqrt{n}(\\overline{\\theta}_n - \\theta) \\leadsto N\\left(0, \\sigma^2 \\frac{\\pi}{2} \\right) $$\n",
    "\n",
    "This means that the median converges to the right value but has a larger variance than the MLE."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More generally, consider two estimators $T_n$ and $U_n$ and suppose that\n",
    "\n",
    "$$\n",
    "\\sqrt{n}(T_n - \\theta) \\leadsto N(0, t^2) \n",
    "\\quad \\text{and} \\quad \n",
    "\\sqrt{n}(U_n - \\theta) \\leadsto N(0, u^2)\n",
    "$$\n",
    "\n",
    "We define the **asymptotic relative efficiency** of U to T by $ARE(U, T) = t^2/u^2$.  In the Normal example, $ARE(\\overline{\\theta}_n, \\hat{\\theta}_n) = 2 / \\pi = 0.63$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 10.23**.  If $\\hat{\\theta}_n$ is the MLE and $\\overline{\\theta}_n$ is any other estimator then\n",
    "\n",
    "$$ ARE(\\overline{\\theta}_n, \\hat{\\theta}_n) \\leq 1 $$\n",
    "\n",
    "Thus, MLE has the smallest (asymptotic) variance and we say that MLE is **efficient** or **asymptotically optimal**.\n",
    "\n",
    "The result is predicated over the model being correct -- otherwise the MLE may no longer be optimal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10.9 The Delta Method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $\\tau = g(\\theta)$ where $g$ is a smooth function.  The maximum likelihood estimator of $\\tau$ is $\\hat{\\tau} = g(\\hat{\\theta})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 10.24 (The Delta Method)**.  If $\\tau = g(\\theta)$ where $g$ is differentiable and $g'(\\theta) \\neq 0$ then\n",
    "\n",
    "$$ \\frac{\\sqrt{n}(\\hat{\\tau}_n - \\tau)}{\\hat{\\text{se}}(\\hat{\\tau})} \\leadsto N(0, 1) $$\n",
    "\n",
    "where $\\hat{\\tau}_n = g(\\hat{\\theta})$ and\n",
    "\n",
    "$$ \\hat{\\text{se}}(\\hat{\\tau}_n) = |g'(\\hat{\\theta})| \\hat{\\text{se}} (\\hat{\\theta}_n) $$\n",
    "\n",
    "Hence, if\n",
    "\n",
    "$$ C_n = \\left( \\hat{\\tau}_n - z_{\\alpha/2} \\hat{\\text{se}}(\\hat{\\tau}_n), \\; \\hat{\\tau}_n + z_{\\alpha/2} \\hat{\\text{se}}(\\hat{\\tau}_n) \\right) $$\n",
    "\n",
    "then $\\mathbb{P}_\\theta(\\tau \\in C_n) \\rightarrow 1 - \\alpha$ as $n \\rightarrow \\infty$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10.10 Multiparameter Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can extend these ideas to models with several parameters.\n",
    "\n",
    "Let $\\theta = (\\theta_1, \\dots, \\theta_n)$ and let $\\hat{\\theta} = (\\hat{\\theta}_1, \\dots, \\hat{\\theta}_n)$ be the MLE.  Let $\\ell_n = \\sum_{i=1}^n \\log f(X_i; \\theta)$,\n",
    "\n",
    "$$\n",
    "H_jj = \\frac{\\partial^2 \\ell_n}{\\partial \\theta_j^2}\n",
    "\\quad \\text{and} \\quad\n",
    "H_jk = \\frac{\\partial^2 \\ell_n}{\\partial \\theta_j \\partial \\theta_k}\n",
    "$$\n",
    "\n",
    "Define the **Fisher Information Matrix** by\n",
    "\n",
    "$$\n",
    "I_n(\\theta) = -\n",
    "\\begin{bmatrix}\n",
    "\\mathbb{E}_\\theta(H_{11}) & \\mathbb{E}_\\theta(H_{12}) & \\cdots & \\mathbb{E}_\\theta(H_{1k}) \\\\\n",
    "\\mathbb{E}_\\theta(H_{21}) & \\mathbb{E}_\\theta(H_{22}) & \\cdots & \\mathbb{E}_\\theta(H_{2k}) \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "\\mathbb{E}_\\theta(H_{k1}) & \\mathbb{E}_\\theta(H_{k2}) & \\cdots & \\mathbb{E}_\\theta(H_{kk})\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Let $J_n(\\theta) = I_n^{-1}(\\theta)$ be the inverse of $I_n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 10.27**. Under appropriate regularity conditions,\n",
    "\n",
    "$$ \\sqrt(n)(\\hat{\\theta} - \\theta) \\approx N(0, J_n(\\theta))$$\n",
    "\n",
    "Also, if $\\hat{\\theta}_j$ is the $j$-th component of $\\hat{\\theta}$, then\n",
    "\n",
    "$$ \\frac{\\sqrt{n}(\\hat{\\theta_j} - \\theta_j)}{\\hat{\\text{se}}_j} \\approx N(0, 1) $$\n",
    "\n",
    "where $\\hat{\\text{se}}_j^2$ is the $j$-th diagonal element of $J_n$.  The approximate covariance of $\\hat{\\theta}_j$ and $\\hat{\\theta}_k$ is $\\text{Cov}(\\hat{\\theta}_j, \\hat{\\theta}_k) \\approx J_n(j, k)$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is also a multiparameter delta method.  Let $\\tau = g(\\theta_1, \\dots, \\theta_k)$ be a function and let\n",
    "\n",
    "$$ \\nabla g = \\begin{pmatrix}\n",
    "\\frac{\\partial g}{\\partial \\theta_1} \\\\\n",
    "\\vdots \\\\\n",
    "\\frac{\\partial g}{\\partial \\theta_k}\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "be the gradient of $g$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 10.28 (Multiparameter delta method)**.  Suppose that $\\nabla g$ evaluated at $\\hat{\\theta}$ is not 0. Let $\\hat{\\tau} = g(\\hat{\\theta})$.  Then\n",
    "\n",
    "$$ \\frac{\\sqrt{n}(\\hat{\\tau} - \\tau)}{\\hat{\\text{se}}(\\hat{\\tau})} \\leadsto N(0, 1) $$\n",
    "\n",
    "where\n",
    "\n",
    "$$ \\hat{\\text{se}}(\\hat{\\tau}) = \\sqrt{\\left(\\hat{\\nabla} g \\right)^T \\hat{J}_n \\left(\\hat{\\nabla} g \\right)} ,$$\n",
    "\n",
    "$\\hat{J}_n = J_n(\\hat{\\theta}_n)$ and $\\hat{\\nabla}g$ is $\\nabla g$ evaluated at $\\theta = \\hat{\\theta}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10.11 The Parametric Bootstrap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For parametric models, standard errors and confidence intervals may also be estimated using the bootstrap.  There is only one change.  In nonparametric bootstrap, we sampled $X_1^*, \\dots, X_n*$ from the empirical distribution $\\hat{F}_n$. In the parametric bootstrap we sample instead from $f(x; \\hat{\\theta}_n)$.  Here, $\\hat{\\theta}_n$ could be the MLE or the method of moments estimator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10.12 Technical Appendix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 10.12.1 Proofs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proof of Theorem 10.13**.  Since $\\hat{\\theta}_n$ maximizes $M_n(\\theta)$, we have $M_n(\\hat{\\theta}) \\geq M_n(\\theta_*)$.  Hence,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "M(\\theta_*) - M(\\hat{\\theta}_n) \n",
    "&= M_n(\\theta_*) - M(\\hat{\\theta}_n) + M(\\hat{\\theta}_*) - M_n(\\theta_*) \\\\\n",
    "&\\leq M_n(\\hat{\\theta}) - M(\\hat{\\theta}_n) + M(\\theta_*) - M_n(\\theta_*) \\\\\n",
    "&\\leq \\sup_\\theta | M_n(\\theta) - M(\\theta) |  + M(\\theta_*)  - M_n(\\theta_*) \\\\\n",
    "&\\xrightarrow{\\text{P}} 0\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "It follows that, for any $\\delta > 0$, \n",
    "\n",
    "$$\\mathbb{P}(M(\\hat{\\theta}_n) < M(\\theta_*) - \\delta) \\rightarrow 0$$\n",
    "\n",
    "Pick any $\\epsilon > 0$.  There exists $\\delta > 0$ such that $|\\theta - \\theta_*| \\geq \\epsilon$ implies that $M(\\theta) < M(\\theta_*) - \\delta$.  Hence,\n",
    "\n",
    "$$\\mathbb{P}(|\\hat{\\theta}_n - \\theta_*| > \\epsilon) \\leq \n",
    "\\mathbb{P}\\left( M(\\hat{\\theta}_n) < M(\\theta_*) - \\delta \\right) \\rightarrow 0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Lemma 10.31**.  The score function satisfies\n",
    "\n",
    "$$\\mathbb{E}[s(X; \\theta)] = 0$$\n",
    "\n",
    "**Proof**.  Note that $1 = \\int f(x; \\theta) dx$.  Differentiate both sides of this equation to get\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "0 &= \\frac{\\partial}{\\partial \\theta} \\int f(x; \\theta)dx = \\int \\frac{\\partial}{\\partial \\theta} f(x; \\theta) dx \\\\\n",
    "&= \\int \\frac{\\frac{\\partial f(x; \\theta)}{\\partial \\theta}}{f(x; \\theta)} f(x; \\theta) dx\n",
    "= \\int \\frac{\\partial \\log f(x; \\theta)}{\\partial \\theta} f(x; \\theta) dx \\\\\n",
    "&= \\int s(x; \\theta) f(x; \\theta) dx = \\mathbb{E}[s(X; \\theta)]\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proof of Theorem 10.18**.  Let $\\ell(\\theta) = \\log \\mathcal{L}(\\theta)$.  Then\n",
    "\n",
    "$$0 = \\ell'(\\hat{\\theta}) \\approx \\ell'(\\theta) + (\\hat{\\theta} - \\theta) \\ell''(\\theta)$$\n",
    "\n",
    "Rearrange the above equation to get $\\hat{\\theta} - \\theta = -\\ell'(\\theta) / \\ell''(\\theta)$, or\n",
    "\n",
    "$$ \\sqrt{n}(\\hat{\\theta} - \\theta) = \\frac{\\frac{1}{\\sqrt{n}}\\ell'(\\theta)}{-\\frac{1}{n}\\ell''(\\theta)} = \\frac{\\text{TOP}}{\\text{BOTTOM}}$$\n",
    "\n",
    "Let $Y_i = \\partial \\log f(X_i, \\theta) / \\partial \\theta$.  From the previous lemma $\\mathbb{E}(Y_i) = 0$ and also $\\mathbb{V}(Y_i) = I(\\theta)$.  Hence,\n",
    "\n",
    "$$\\text{TOP} = n^{-1/2} \\sum_i Y_i = \\sqrt{n} \\overline{Y} = \\sqrt{n} (\\overline{Y} - 0) \\leadsto W \\sim N(0, I)$$\n",
    "\n",
    "by the central limit theorem.  Let $A_i = -\\partial^2 \\log f(X_i; \\theta) / \\partial theta^2$.  Then $\\mathbb{E}(A_i) = I(\\theta)$ and\n",
    "\n",
    "$$\\text{BOTTOM} = \\overline{A} \\xrightarrow{\\text{P}} I(\\theta)$$\n",
    "\n",
    "by the law of large numbers.  Apply Theorem 6.5 part (e) to conclude that\n",
    "\n",
    "$$\\sqrt{n}(\\hat{\\theta} - \\theta) \\leadsto \\frac{W}{I(\\theta)} \\sim N \\left(0, \\frac{1}{I(\\theta)} \\right)$$\n",
    "\n",
    "Assuming that $I(\\theta)$ is a continuous function of $\\theta$, it follows that $I(\\hat{\\theta}_n) \\xrightarrow{\\text{P}} I(\\theta)$.  Now\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\frac{\\hat{\\theta}_n - \\theta}{\\hat{\\text{se}}}&= \\sqrt{n} I^{1/2}(\\hat{\\theta_n})(\\hat{\\theta_n} - \\theta) \\\\\n",
    "&= \\left\\{ \\sqrt{n} I^{1/2}(\\theta)(\\hat{\\theta}_n - \\theta)\\right\\} \\left\\{ \\frac{I(\\hat{\\theta}_n)}{I(\\theta)} \\right\\}^{1/2}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The first term tends in distribution to $N(0, 1)$.  The second term tends in probability to 1.  The result follows from Theorem 6.5 part (e)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Outline of proof of Theorem 10.24**.  Write,\n",
    "\n",
    "$$\\hat{\\tau} = g(\\hat{\\theta}) \\approx g(\\theta) + (\\hat{\\theta} - \\theta)g'(\\theta) = \\tau + (\\hat{\\theta} - \\theta)g'(\\theta)$$\n",
    "\n",
    "Thus,\n",
    "\n",
    "$$\\sqrt{n}(\\hat{\\tau} - \\tau) \\approx \\sqrt{n}(\\hat{\\theta} - \\theta)g'(\\theta)$$\n",
    "\n",
    "and hence\n",
    "\n",
    "$$\\frac{\\sqrt{n}I^{1/2}(\\theta)(\\hat{\\theta} - \\theta)}{g'(\\theta)} \\approx \\sqrt{n}I^{1/2}(\\theta)(\\hat{\\theta} - \\theta)$$\n",
    "\n",
    "Theorem 10.18 tells us that the right hand side tends in distribution to $N(0, 1)$, hence\n",
    "\n",
    "$$\\frac{\\sqrt{n}I^{1/2}(\\theta)(\\hat{\\theta} - \\theta)}{g'(\\theta)} \\leadsto N(0, 1)$$\n",
    "\n",
    "or, in other words,\n",
    "\n",
    "$$\\hat{\\tau} \\approx N(\\tau, \\text{se}^2(\\hat{\\tau}_n))$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\\text{se}^2(\\hat{\\tau}_n) = \\frac{(g'(\\theta))^2}{nI(\\theta)}$$\n",
    "\n",
    "The result remains true if we substitute $\\hat{\\theta}$ for $\\theta$ by Theorem 6.5 part (e)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 10.12.2 Sufficiency"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A **statistic** is a function $T(X^n)$ of the data.  A sufficient statistic is a statistic that contains all of the information in the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write $x^n \\leftrightarrow y^n$ if $f(x^n; \\theta) = c f(y^n; \\theta)$ for some constant $c$ that might depend on $x^n$ and $y^n$ but not $\\theta$.  A statistic is **sufficient** if $T(x^n) \\leftrightarrow T(y^n)$ implies that $x^n \\leftrightarrow y^n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that if $x^n \\leftrightarrow y^n$ then the likelihood functions based on $x^n$ and $y^n$ have the same shape.  Roughly speaking, a statistic is sufficient if we can calculate the likelihood function knowing only $T(X^n)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A statistic $T$ is **minimally sufficient** if it is sufficient and it is a function of every other sufficient statistic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 10.36**. $T$ is minimally sufficient if $T(x^n) = T(y^n)$ if and only if $x^n \\leftrightarrow y^n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The usual definition of sufficiency is this: $T$ is sufficient if the distribution of $X^n$ given $T(X^n) = t$ does not depend on $\\theta$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 10.40 (Factorization Theorem)**.  $T$ is sufficient if and only if there are functions $g(t, \\theta)$ and $h(x)$ such that $f(x^n; \\theta) = g(t(x^n); \\theta)h(x^n)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 10.42 (Rao-Blackwell)**.  Let $\\hat{\\theta}$ be an estimator and let $T$ be a sufficient statistic.  Define a new estimator by\n",
    "\n",
    "$$\\overline{\\theta} = \\mathbb{E}(\\hat{\\theta} | T)$$\n",
    "\n",
    "Then, for every $\\theta$, \n",
    "\n",
    "$$R(\\theta, \\overline{\\theta}) \\leq R(\\theta, \\hat{\\theta})$$\n",
    "\n",
    "where $R(\\theta, \\hat{\\theta}) = \\mathbb{E}_\\theta[(\\theta - \\hat{\\theta})^2]$ denote the MSE of an estimator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 10.12.3 Exponential Families"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We say that $\\{f(x; \\theta) : \\theta \\in \\Theta\\}$ is a **one-parameter exponential family** if there are functions $\\eta(\\theta)$, $B(\\theta)$, $T(x)$ and $h(x)$ such that\n",
    "\n",
    "$$f(x; \\theta) = h(x) e^{\\eta(\\theta)T(x) - B(\\theta)}$$\n",
    "\n",
    "It is easy to see that $T(X)$ is sufficient.  We call $T$ the **natural sufficient statistic**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can rewrite an exponential family as\n",
    "\n",
    "$$f(x; \\eta) = h(x) e^{\\eta T(x) - A(\\eta)}$$\n",
    "\n",
    "where $\\eta = \\eta(\\theta)$ is called the **natural parameter** and\n",
    "\n",
    "$$A(\\eta) = \\log \\int h(x) e^{\\eta T(x)} dx$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $X_1, \\dots, X_n$ be iid from an exponential family. Then $f(x^n; \\theta)$ is an exponential family:\n",
    "\n",
    "$$f(x^n; \\theta) = h_n(x^n) e^{\\eta(\\theta) T_n(x^n) - B_n(\\theta)}$$\n",
    "\n",
    "where $h_n(x^n) = \\prod_i h(x_i)$, $T_n(x^n) = \\sum_i T(x_i)$ and $B_n(\\theta) = nB(\\theta)$.  This implies that $\\sum_i T(X_i)$ is sufficient."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 10.47**. Let $X$ have an exponential family.  Then,\n",
    "\n",
    "$$\n",
    "\\mathbb{E}(T(X)) = A'(\\eta),\n",
    "\\quad\n",
    "\\mathbb{V}(T(X)) = A''(\\eta)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $\\theta = (\\theta_1, \\dots, \\theta_n)$ is a vector, then we say that $f(x; \\theta)$ has exponential family form if\n",
    "\n",
    "$$ f(x; \\theta) = h(x) \\exp \\left\\{ \\sum_{j=1}^k \\eta_j(\\theta) T_j(x) - B(\\theta) \\right\\}$$\n",
    "\n",
    "Again, $T = (T_1, \\dots, T_k)$ is sufficient and $n$ iid samples also has exponential form with sufficient statistic $\\left(\\sum_i T_1(X_i), \\dots, \\sum_i T_k(X_i)\\right)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 10.13 Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 10.13.1**.  Let $X_1, \\dots, X_n \\sim \\text{Gamma}(\\alpha, \\beta)$.  Find the method of moments estimator for $\\alpha$ and $\\beta$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**. \n",
    "\n",
    "The first two moments are:\n",
    "\n",
    "$$\n",
    "\\alpha_1 = \\mathbb{E}(X) = \\frac{\\alpha}{\\beta}\\\\\n",
    "\\alpha_2 = \\mathbb{E}(X^2) = \\mathbb{V}(X) + \\mathbb{E}(X)^2 = \\frac{\\alpha}{\\beta^2} + \\frac{\\alpha^2}{\\beta^2} = \\frac{\\alpha(\\alpha + 1)}{\\beta^2} \n",
    "$$\n",
    "\n",
    "We have the sample moments:\n",
    "\n",
    "$$\\hat{\\alpha}_1 = \\frac{1}{n}\\sum_{i=1}^n X_i\n",
    "\\quad \\quad\n",
    "\\hat{\\alpha}_2 = \\frac{1}{n}\\sum_{i=1}^n X_i^2$$\n",
    "\n",
    "Equating these we get:\n",
    "\n",
    "$$\n",
    "\\hat{\\alpha}_1 = \\frac{\\hat{\\alpha}_n}{\\hat{\\beta}_n}\n",
    "\\quad \\quad\n",
    "\\hat{\\alpha}_2 = \\frac{\\hat{\\alpha}_n(\\hat{\\alpha}_n + 1)}{\\hat{\\beta}_n^2}\n",
    "$$\n",
    "\n",
    "Solving these we get the method of moments estimators for $\\alpha$ and $\\beta$:\n",
    "\n",
    "$$\n",
    "\\hat{\\alpha}_n = \\frac{\\hat{\\alpha}_1^2}{\\hat{\\alpha}_2 - \\hat{\\alpha}_1^2}\n",
    "\\quad \\quad\n",
    "\\hat{\\beta}_n = \\frac{\\hat{\\alpha}_1}{\\hat{\\alpha}_2 - \\hat{\\alpha}_1^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 10.13.2**. Let $X_1, \\dots, X_n \\sim \\text{Uniform}(a, b)$ where $a, b$ are unknown parameters and $a < b$.\n",
    "\n",
    "**(a)** Find the method of moments estimators for $a$ and $b$.\n",
    "\n",
    "**(b)** Find the MLE $\\hat{a}$ and $\\hat{b}$.\n",
    "\n",
    "**(c)** Let $\\tau = \\int x dF(x)$.  Find the MLE of $\\tau$.\n",
    "\n",
    "**(d)** Let $\\hat{\\tau}$ be the MLE from the previous item.  Let $\\tilde{\\tau}$ be the nonparametric plug-in estimator of $\\tau = \\int x dF(x)$.  Suppose that $a = 1$, $b = 3$ and $n = 10$.  Find the MSE of $\\hat{\\tau}$ by simulation.  Find the MSE of $\\tilde{\\tau}$ analytically.  Compare."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "**(a)**\n",
    "\n",
    "The first two moments are:\n",
    "\n",
    "$$ \n",
    "\\alpha_1 = \\mathbb{E}(X) = \\frac{a + b}{2} \\\\\n",
    "\\alpha_2 = \\mathbb{E}(X^2) = \\mathbb{V}(X) + \\mathbb{E}(X)^2 = \\frac{(b - a)^2}{12} + \\frac{(a + b)^2}{4}\n",
    "= \\frac{a^2 + ab + b^2}{3}\n",
    "$$\n",
    "\n",
    "We have the sample moments:\n",
    "\n",
    "$$\\hat{\\alpha}_1 = \\frac{1}{n}\\sum_{i=1}^n X_i\n",
    "\\quad \\quad\n",
    "\\hat{\\alpha}_2 = \\frac{1}{n}\\sum_{i=1}^n X_i^2$$\n",
    "\n",
    "Equating these we get:\n",
    "\n",
    "$$\n",
    "\\hat{\\alpha}_1 = \\frac{\\hat{a} + \\hat{b}}{2}\n",
    "\\quad \\quad\n",
    "\\hat{\\alpha}_2 = \\frac{(\\hat{b} - \\hat{a})^2}{12} + \\frac{(\\hat{a} + \\hat{b})^2}{4}\n",
    "$$\n",
    "\n",
    "Solving these we get the method of moment estimators for $a$ and $b$:\n",
    "\n",
    "$$\n",
    "\\hat{a} = \\hat{\\alpha}_1 - \\sqrt{3}(\\hat{\\alpha}_1^2 - \\hat{\\alpha}_2)\n",
    "\\quad \\quad\n",
    "\\hat{b} = \\hat{\\alpha}_1 + \\sqrt{3}(\\hat{\\alpha}_1^2 - \\hat{\\alpha}_2)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(b)**\n",
    "\n",
    "The probability density function for each $X_i$ is\n",
    "\n",
    "$$f(x; (a, b)) = \\begin{cases}\n",
    "(b - a)^{-1} & \\text{if } a \\leq x \\leq b \\\\\n",
    "0 & \\text{otherwise}\n",
    "\\end{cases}$$\n",
    "\n",
    "The likelihood function is\n",
    "\n",
    "$$\\mathcal{L}_n(a, b) = \\prod_{i=1}^n f(X_i; (a, b)) = \\begin{cases}\n",
    "(b-a)^{-n} & \\text{if } a \\leq X_i \\leq b \\text{ for all } X_i\\\\\n",
    "0 & \\text{otherwise}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "The parameters that maximize the likelihood function make the $b - a$ as small as possible -- that is, we should pick the maximum $a$ and the minimum $b$ for which the likelihood function is non-zero.  So the MLEs are:\n",
    "\n",
    "$$\\hat{a} = \\min \\{X_1, \\dots, X_n \\}\n",
    "\\quad \\quad\n",
    "\\hat{b} = \\max \\{X_1, \\dots, X_n \\}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(c)**\n",
    "\n",
    "$\\tau = \\int x dF(x) = \\mathbb{E}(x) = (a + b)/2$, so since the MLE is equivariant, the MLE of $\\tau$ is \n",
    "\n",
    "$$\\hat{\\tau} = \\frac{\\hat{a} + \\hat{b}}{2} = \\frac{\\min \\{X_1, \\dots, X_n\\} + \\max\\{X_1, \\dots, X_n\\}}{2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(d)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "a = 1\n",
    "b = 3\n",
    "n = 10\n",
    "\n",
    "X = np.random.uniform(low=a, high=b, size=n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE for tau_hat: \t 0.150\n"
     ]
    }
   ],
   "source": [
    "tau_hat = (X.min() + X.max()) / 2\n",
    "\n",
    "# Nonparametric bootstrap to find MSE of tau_hat\n",
    "B = 10000\n",
    "t_boot = np.empty(B)\n",
    "for i in range(B):\n",
    "    xx = np.random.choice(X, n, replace=True)\n",
    "    t_boot[i] = (xx.min() + xx.max()) / 2\n",
    "    \n",
    "se = t_boot.std()\n",
    "print(\"MSE for tau_hat: \\t %.3f\" % se)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analytically, we have:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{V}(\\tilde{\\tau}) \n",
    "&= \\mathbb{E}(\\tilde{\\tau}^2) - (\\mathbb{E}(\\tilde{\\tau}))^2 \\\\\n",
    "&= \\frac{1}{n^2}\\left(\\mathbb{E}\\left[ \\left(\\sum_{i=1}^n X_i\\right)^2\\right] - \\left(\\mathbb{E}\\left[\\sum_{i=1}^n X_i \\right]\\right)^2\\right) \\\\\n",
    "&= \\frac{1}{n^2}\\left( \\mathbb{E}\\left[ \\sum_{i=1}^n X_i^2 + \\sum_{i=1}^n \\sum_{j=1, j \\neq i}^n X_i X_j \\right] - \\left(n \\frac{a + b}{2}\\right)^2\\right) \\\\\n",
    "&= \\frac{1}{n^2}\\left( \\sum_{i=1}^n \\mathbb{E}[X_i^2] + \\sum_{i=1}^n \\sum_{j=1, j \\neq i}^n \\mathbb{E}[X_i]\\mathbb{E}[X_j]  - \\left(n \\frac{a + b}{2}\\right)^2\\right) \\\\\n",
    "&= \\frac{1}{n^2}\\left( n \\frac{a^2 + ab + b^2}{3} + n(n-1) \\left(\\frac{a+b}{2}\\right)^2  - n^2\\left(\\frac{a + b}{2}\\right)^2\\right) \\\\\n",
    "&= \\frac{1}{n^2}\\left( n \\frac{a^2 + ab + b^2}{3} - n \\left(\\frac{a+b}{2}\\right)^2 \\right) \\\\\n",
    "&= \\frac{1}{n} \\left( \\frac{a^2 + ab + b^2}{3} - \\frac{a^2 + 2ab + b^2}{4}\\right) \\\\\n",
    "&= \\frac{1}{n} \\frac{(b - a)^2}{12}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Therefore,\n",
    "\n",
    "$$\\text{se}(\\tilde{\\tau}) = \\sqrt{\\frac{1}{n} \\frac{(b - a)^2}{12}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE for tau_tilde: \t 0.183\n"
     ]
    }
   ],
   "source": [
    "se_tau_tilde = np.sqrt((1/n) * ((b - a)**2 / 12))\n",
    "\n",
    "print(\"MSE for tau_tilde: \\t %.3f\" % se_tau_tilde)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 10.13.3**.  Let $X_1, \\dots, X_n \\sim N(\\mu, \\sigma^2)$.  Let $\\tau$ be the 0.95 percentile, i.e. $\\mathbb{P}(X < \\tau) = 0.95$.\n",
    "\n",
    "**(a)** Find the MLE of $\\tau$.\n",
    "\n",
    "**(b)** Find an expression for an approximate $1 - \\alpha$ confidence interval for $\\tau$.\n",
    "\n",
    "**(c)** Suppose the data are:\n",
    "\n",
    "$$\n",
    "\\begin{matrix}\n",
    "3.23 & -2.50 &  1.88 & -0.68 &  4.43 & 0.17 \\\\ \n",
    "1.03 & -0.07 & -0.01 &  0.76 &  1.76 & 3.18 \\\\\n",
    "0.33 & -0.31 &  0.30 & -0.61 &  1.52 & 5.43 \\\\\n",
    "1.54 &  2.28 &  0.42 &  2.33 & -1.03 & 4.00 \\\\\n",
    "0.39 \n",
    "\\end{matrix}\n",
    "$$\n",
    "\n",
    "Find the MLE $\\hat{\\tau}$.  Find the standard error using the delta method.  Find the standard error using the parametric bootstrap."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(a)**  \n",
    "\n",
    "Let $Z \\sim N(0, 1)$, so $(X - \\mu) / \\sigma \\sim Z$.  We have:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{P}(X < \\tau) &= 0.95 \\\\\n",
    "\\mathbb{P}\\left(\\frac{X - \\mu}{\\sigma} < \\frac{\\tau - \\mu}{\\sigma}\\right) &= 0.95 \\\\\n",
    "\\mathbb{P}\\left(Z < \\frac{\\tau - \\mu}{\\sigma}\\right) &= 0.95 \\\\\n",
    "\\frac{\\tau - \\mu}{\\sigma} &= z_{5\\%} \\\\\n",
    "\\tau &= \\mu + z_{5\\%} \\sigma \n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Since the MLE is equivariant, $\\hat{\\tau} = \\hat{\\mu} + z_{5\\%} \\hat{\\sigma}$, where $\\hat{\\mu}, \\hat{\\sigma}$ are the MLEs for the Normal distribution parameters:\n",
    "\n",
    "$$ \\hat{\\mu} = n^{-1} \\sum_{i=1}^n X_i\n",
    "\\quad \\quad\n",
    "\\hat{\\sigma} = \\sqrt{n^{-1} \\sum_{i=1}^n (X_i - \\overline{X})^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(b)**\n",
    "\n",
    "Let's use the multiparameter delta method.  \n",
    "\n",
    "We have $\\tau = g(\\mu, \\sigma) = \\mu + z_{5\\%} \\sigma$, so\n",
    "\n",
    "$$ \\nabla g = \\begin{bmatrix}\n",
    "\\partial g / \\partial \\mu \\\\\n",
    "\\partial g / \\partial \\sigma\n",
    "\\end{bmatrix}\n",
    "= \\begin{bmatrix}\n",
    "1 \\\\\n",
    "z_{5\\%}\n",
    "\\end{bmatrix}$$.\n",
    "\n",
    "The Fisher Information Matrix for the Normal process is\n",
    "\n",
    "$$ I_n(\\mu, \\sigma) = \\begin{bmatrix}\n",
    "n / \\sigma^2 & 0 \\\\\n",
    "0 & 2n / \\sigma^2\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "then its inverse is\n",
    "\n",
    "$$ J_n = I_n^{-1}(\\mu, \\sigma) = \\frac{1}{n} \\begin{bmatrix}\n",
    "\\sigma^2 & 0 \\\\\n",
    "0 & \\sigma^2/2\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "and the standard error estimate for our new parameter variable is\n",
    "\n",
    "$$\\hat{\\text{se}}(\\hat{\\tau}) = \\sqrt{(\\hat{\\nabla} g)^T \\hat{J}_n (\\hat{\\nabla} g)} = \\hat{\\sigma} \\sqrt{n^{-1}(1 + z_{5\\%}^2 / 2)}$$\n",
    "\n",
    "A $1 - \\alpha$ confidence interval for $\\hat{\\tau}$, then, is\n",
    "\n",
    "$$ C_n = \\left(\n",
    "\\hat{\\mu} + \\hat{\\sigma}\\left( z_{5\\%} - z_{\\alpha / 2} \\sqrt{n^{-1}(1 + z_{5\\%}^2 / 2)} \\right), \\;\n",
    "\\hat{\\mu} + \\hat{\\sigma}\\left( z_{5\\%} + z_{\\alpha / 2} \\sqrt{n^{-1}(1 + z_{5\\%}^2 / 2)} \\right) \\right) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(c)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "\n",
    "z_05 = norm.ppf(0.95)\n",
    "z_025 = norm.ppf(0.975)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.array([\n",
    "    3.23, -2.50,  1.88, -0.68,  4.43, 0.17,\n",
    "    1.03, -0.07, -0.01,  0.76,  1.76, 3.18,\n",
    "    0.33, -0.31,  0.30, -0.61,  1.52, 5.43,\n",
    "    1.54,  2.28,  0.42,  2.33, -1.03, 4.00,\n",
    "    0.39   \n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated tau: 4.180\n"
     ]
    }
   ],
   "source": [
    "# Estimate the MLE tau_hat\n",
    "\n",
    "n = len(X)\n",
    "mu_hat = X.mean()\n",
    "sigma_hat = X.std()\n",
    "tau_hat = mu_hat + z_05 * sigma_hat\n",
    "\n",
    "print(\"Estimated tau: %.3f\" % tau_hat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated tau (delta method, 95% confidence interval): \t (3.088, 5.273)\n"
     ]
    }
   ],
   "source": [
    "# Confidence interval using delta method\n",
    "\n",
    "se_tau_hat = sigma_hat * np.sqrt((1/n) * (1 + z_05 * z_05 / 2))\n",
    "confidence_interval = (tau_hat - z_025 * se_tau_hat, tau_hat + z_025 * se_tau_hat)\n",
    "\n",
    "print(\"Estimated tau (delta method, 95%% confidence interval): \\t (%.3f, %.3f)\" % confidence_interval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated tau (parametric bootstrap, 95% confidence interval): \t (2.887, 5.474)\n"
     ]
    }
   ],
   "source": [
    "# Confidence interval using parametric bootstrap\n",
    "\n",
    "n = len(X)\n",
    "mu_hat = X.mean()\n",
    "sigma_hat = X.std()\n",
    "tau_hat = mu_hat + z_05 * sigma_hat\n",
    "\n",
    "B = 10000\n",
    "t_boot = np.empty(B)\n",
    "for i in range(B):\n",
    "    xx = norm.rvs(loc=mu_hat, scale=sigma_hat, size=n)\n",
    "    t_boot[i] = np.quantile(xx, 0.95)\n",
    "    \n",
    "se_tau_hat_bootstrap = t_boot.std()\n",
    "confidence_interval = (tau_hat - z_025 * se_tau_hat_bootstrap, tau_hat + z_025 * se_tau_hat_bootstrap)\n",
    "\n",
    "print(\"Estimated tau (parametric bootstrap, 95%% confidence interval): \\t (%.3f, %.3f)\" % confidence_interval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 10.13.4**  Let $X_1, \\dots, X_n \\sim \\text{Uniform}(0, \\theta)$.  Show that the MLE is consistent.\n",
    "\n",
    "Hint: Let $Y = \\max \\{ X_1, \\dots, X_n \\}$.  For any c, $\\mathbb{P}(Y < c) = \\mathbb{P}(X_1 < c, X_2 < c, \\dots, X_n < c) = \\mathbb{P}(X_1 < c)\\mathbb{P}(X_2 < c)\\dots\\mathbb{P}(X_n < c)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "The probability density function is\n",
    "\n",
    "$$ f(x, \\theta) = \\mathbb{P}(Y < x) = \\prod_{i = 1}^n \\mathbb{P}(X_i < x) = f_{\\text{Uniform}(0, \\theta)}(x)^n $$\n",
    "\n",
    "The probability density function for the original distribution is\n",
    "\n",
    "$$ f_{\\text{Uniform}(0, \\theta)}(x) = \\begin{cases}\n",
    "\\theta^{-1} & \\text{if } 0 \\leq x \\leq \\theta \\\\\n",
    "0 & \\text{otherwise}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "so\n",
    "\n",
    "$$ f(x, \\theta) = \\begin{cases}\n",
    "\\theta^{-n} & \\text{if } 0 \\leq x \\leq \\theta \\\\\n",
    "0 & \\text{otherwise}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "The likelihood is maximized when $\\theta$ is as small as possible while keeping all samples within the first case, so $\\hat{\\theta}_n = \\max \\{X_1, \\dots, X_n \\}$.\n",
    "\n",
    "For a given $\\epsilon > 0$, we have\n",
    "\n",
    "$$\\mathbb{P}(\\hat{\\theta}_n < \\theta - \\epsilon) = \\prod_{i=1}^n \\mathbb{P}(X_i < \\theta - \\epsilon) = \\left(1 - \\frac{\\epsilon}{\\theta} \\right)^n$$\n",
    "\n",
    "which goes to 0 as $n \\rightarrow \\infty$, so $\\lim _{n \\rightarrow \\infty} \\hat{\\theta}_n = \\theta$, and thus the MLE is consistent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 10.13.5**.  Let $X_1, \\dots, X_n \\sim \\text{Poisson}(\\lambda)$.  Find the method of moments estimator, the maximum likelihood estimator, and the Fisher information $I(\\lambda)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "The first moment is:\n",
    "\n",
    "$$\\mathbb{E}(X) = \\lambda$$\n",
    "\n",
    "We have the sample moment:\n",
    "\n",
    "$$\\hat{\\alpha}_1 = \\frac{1}{n} \\sum_{i=1}^n X_i$$\n",
    "\n",
    "Equating those, the method of moments estimator for $\\hat{\\lambda}$ is:\n",
    "\n",
    "$$\\hat{\\lambda} = \\hat{\\alpha_1} = \\frac{1}{n} \\sum_{i=1}^n X_i$$\n",
    "\n",
    "The likelihood function is\n",
    "\n",
    "$$\\mathcal{L}_n(\\lambda) = \\prod_{i=1}^n f(X_i; \\lambda) = \\prod_{i=1}^n \\frac{\\lambda^{X_i}e^{-\\lambda}}{(X_i)!}$$\n",
    "\n",
    "so the log likelihood function is\n",
    "\n",
    "$$\\ell_n(\\lambda) = \\log \\mathcal{L}_n(\\lambda) = \\sum_{i=1}^n (\\log(\\lambda^{X_i}e^{-\\lambda}) - \\log X_i!)\n",
    "= \\sum_{i=1}^n (X_i \\log \\lambda - \\lambda - \\log X_i!)\n",
    "= -n \\lambda + (\\log \\lambda) \\sum_{i=1}^n X_i - \\sum_{i=1}^n \\log X_i!\n",
    "$$\n",
    "\n",
    "To find the MLE, we differentiate this equation with respect to 0 and equate it to 0:\n",
    "\n",
    "$$ \\frac{\\partial \\ell_n(\\lambda)}{\\partial \\lambda} = 0 \\\\\n",
    "-n + \\frac{\\sum_{i=1}^n X_i}{\\hat{\\lambda}} = 0 \\\\\n",
    "\\hat{\\lambda} = \\frac{1}{n} \\sum_{i=1}^n X_i\n",
    "$$\n",
    "\n",
    "The score function is:\n",
    "\n",
    "$$ s(X; \\lambda) = \\frac{\\partial \\log f(X; \\lambda)}{\\partial \\lambda} = \\frac{X}{\\lambda} - 1$$\n",
    "\n",
    "and the Fisher information is:\n",
    "\n",
    "$$ I_n(\\lambda) = \\sum_{i=1}^n \\mathbb{V}\\left( s(X_i; \\lambda) \\right) \n",
    "= \\sum_{i=1}^n \\mathbb{V} \\left( \\frac{X_i}{\\lambda} - 1 \\right)\n",
    "= \\frac{1}{\\lambda^2}  \\sum_{i=1}^n \\mathbb{V}(X_i) = \\frac{n}{\\lambda}$$\n",
    "\n",
    "In particular, $I(\\lambda) = I_1(\\lambda) = 1 / \\lambda$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 10.13.6**.  Let $X_1, \\dots, X_n \\sim N(\\theta, 1)$.  Define\n",
    "\n",
    "$$Y_i = \\begin{cases}\n",
    "1 & \\text{if } X_i > 0 \\\\\n",
    "0 & \\text{if } X_i \\leq 0\n",
    "\\end{cases}$$\n",
    "\n",
    "Let $\\psi = \\mathbb{P}(Y_1 = 1)$.\n",
    "\n",
    "**(a)**  Find the maximum likelihood estimate $\\hat{\\psi}$ of $\\psi$.\n",
    "\n",
    "**(b)**  Find an approximate 95% confidence interval for $\\psi$.\n",
    "\n",
    "**(c)**  Define $\\overline{\\psi} = (1 / n) \\sum_i Y_i$.  Show that $\\overline{\\psi}$ is a consistent estimator of $\\psi$.\n",
    "\n",
    "**(d)**  Compute the asymptotic relative efficiency of $\\overline{\\psi}$ to $\\hat{\\psi}$.  Hint:  Use the delta method to get the standard error of the MLE.  Then compute the standard error (i.e. the standard deviation) of $\\overline{\\psi}$.\n",
    "\n",
    "**(e)**  Suppose that the data are not really normal.  Show that $\\psi$ is not consistent.  What, if anything, does $\\hat{\\psi}$ converge to?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "Note that, from the definition, $Y_1, \\dots, Y_n \\sim \\text{Bernoulli}(\\Phi(\\theta))$, where $\\Phi$ is the CDF for the normal distribution.  Let $p = \\Phi(\\theta)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(a)**  We have $\\psi = \\mathbb{P}(Y_1 = 1) = p$, so the MLE is $\\hat{\\psi} = \\hat{p} = \\Phi(\\hat{\\theta})\n",
    "= \\Phi(\\overline{X})$, where $\\overline{X} = n^{-1} \\sum_{i=1}^n X_i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(b)**  Let $g(\\theta) = \\Phi(\\theta)$.  Then $g'(\\theta) = \\phi(\\theta)$, where $\\phi$ is the standard normal PDF.  By the delta method, $\\hat{\\text{se}}(\\hat{\\psi}) = |g'(\\hat{\\theta})| \\hat{\\text{se}}(\\hat{\\theta}) = \\phi(\\overline{X}) n^{-1/2}$.\n",
    "\n",
    "Then, an approximate 95% confidence interval is\n",
    "\n",
    "$$ C_n = \\left(\\Phi(\\overline{X}) \\left(1 - \\frac{z_{2.5\\%}}{\\sqrt{n}}\\right), \\; \n",
    "\\Phi(\\overline{X}) \\left(1 + \\frac{z_{2.5\\%}}{\\sqrt{n}}\\right) \\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(c)**  $\\overline{\\psi}$ has mean $p$, so consistency follows from the law of large numbers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(d)** We have $\\mathbb{V}(Y_1) = \\psi (1 - \\psi)$, since $Y_1$ follows a Bernoulli distribution, so $\\mathbb{V}(\\overline{\\psi}) = \\mathbb{V}(Y_1) / n = \\psi (1 - \\psi) / n$.\n",
    "\n",
    "From (b), $\\mathbb{V}{\\hat{\\psi}} = \\phi(\\theta) / n$.\n",
    "\n",
    "Therefore, the asymptotic relative efficiency is\n",
    "\n",
    "$$\\frac{\\psi(1 - \\psi)}{\\phi(\\theta)} = \\frac{\\Phi(\\theta)(1 - \\Phi(\\theta))}{\\phi(\\theta)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(e)** By the law of large numbers, we still have that $\\overline{X}$ converges to $\\mathbb{E}(Y_1) = \\mathbb{P}(Y_1 = 1) \\cdot 1 + \\mathbb{P}(Y_1 = 0)\\cdot 0 = \\mathbb{P}(Y_1 = 1) = 1 - F_X(0) = \\mu_Y$.  Then $\\hat{\\psi} = \\Phi(\\overline{X})$ converges to $\\Phi(\\mu_Y)$.  But the true value of $\\psi$ is $\\mathbb{P}(Y_1 = 1) = 1 - F_X(0)$.\n",
    "\n",
    "But for an arbitrary distribution $1 - F_X(0) \\neq \\Phi(1 - F_X(0))$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 10.13.7**. (Comparing two treatments).  $n_1$ people are given treatment 1 and $n_2$ people are given treatment 2.  Let $X_1$ be the number of people on treatment 1 who respond favorably to the treatment and let $X_2$ be the number of people on treatment 2 who respond favorably.  Assume that $X_1 \\sim \\text{Binomial}(n_1, p_1)$, $X_2 \\sim \\text{Binomial}(n_2, p_2)$.  Let $\\psi = p_1 - p_2$.\n",
    "\n",
    "**(a)** Find the MLE of $\\psi$.\n",
    "\n",
    "**(b)** Find the Fisher Information Matrix $I(p_1, p_2)$.\n",
    "\n",
    "**(c)** Use the multiparameter delta method to find the asymptotic standard error of $\\hat{\\psi}$.\n",
    "\n",
    "**(d)** Suppose that $n_1 = n_2 = 200$, $X_1 = 160$ and $X_2 = 148$.  Find $\\hat{\\psi}$.  Find an approximate 90% confidence interval for $\\psi$ using (i) the delta method and (ii) the parametric bootstrap."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "**(a)** The MLE is equivariant, so\n",
    "\n",
    "$$\\hat{\\psi} = \\hat{p_1} - \\hat{p_2} = \\frac{X_1}{n_1} - \\frac{X_2}{n_2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(b)**\n",
    "\n",
    "The probability mass function is\n",
    "\n",
    "$$f((x_1, x_2); \\psi) = f_1(x_1; p_1) f_2(x_2; p_2) = \n",
    "\\binom{n_1}{x_1} p_1^{x_1} (1 - p_1)^{n_1 - x_1}\n",
    "\\binom{n_2}{x_2} p_2^{x_2} (1 - p_2)^{n_2 - x_2}\n",
    "$$\n",
    "\n",
    "The log likelihood is\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\ell_n &= \\log f((x_1, x_2); \\psi) \\\\\n",
    "&= \\sum_{i=1}^2 \\log \\binom{n_i}{x_i} + x_i \\log p_i + (n_i - x_i) \\log (1 - p_i)\n",
    "\\end{align}$$\n",
    "\n",
    "Calculating the partial derivatives and their expectations,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "H_{11} & = \\frac{\\partial^2 \\ell_n}{\\partial p_1^2}\n",
    "= \\frac{\\partial}{\\partial p_1} \\left( \\frac{x_1}{p_1} - \\frac{n_1 - x_1}{1 - p_1}\\right)\n",
    "= -\\frac{x_1}{p_1^2} - \\frac{n_1 - x_1}{(1 - p_1)^2} \\\\\n",
    "\\mathbb{E}[H_{11}] &= -\\frac{\\mathbb{E}[x_1]}{p_1^2} - \\frac{\\mathbb{E}[n - x_1]}{(1 - p_1)^2}\n",
    "= -\\frac{n_1 / p_1}{p_1^2} - \\frac{n_1/(1 - p_1)}{(1 - p_1)^2}\n",
    "= -\\frac{n_1}{p_1} - \\frac{n_1}{1 - p_1} = -\\frac{n_1}{p_1(1 - p_1)}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "$$ \n",
    "\\begin{align}\n",
    "H_{22} &= -\\frac{x_2}{p_2^2} - \\frac{n_2 - x_2}{(1 - p_2)^2} \\\\\n",
    "\\mathbb{E}[H_{22}] &= -\\frac{n_2}{p_2(1 - p_2)}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "$$ H_{12} = \\frac{\\partial^2 \\ell_n}{\\partial p_1 \\partial p_2} = 0$$\n",
    "$$ H_{21} = 0$$\n",
    "\n",
    "So the Fisher Information Matrix is:\n",
    "\n",
    "$$ I(p_1, p_2) = \\begin{bmatrix}\n",
    "\\frac{n_1}{p_1(1 - p_1)} & 0\\\\\n",
    "0 & \\frac{n_2}{p_2(1 - p_2)}\n",
    "\\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(c)** Using the multiparameter delta method, $g(\\psi) = p_1 - p_2$, so\n",
    "\n",
    "$$ \\nabla g = \\begin{bmatrix}\n",
    "\\partial g / \\partial p_1 \\\\ \\partial g / \\partial p_2\n",
    "\\end{bmatrix}\n",
    "= \\begin{bmatrix}\n",
    "1 \\\\ -1\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "The inverse of the Fisher Information Matrix is\n",
    "\n",
    "$$J(p_1, p_2) = I(p_1, p_2)^{-1} = \\begin{bmatrix}\n",
    "\\frac{p_1(1 - p_1)}{n_1} & 0 \\\\\n",
    "0 & \\frac{p_2(1 - p_2)}{n_2}\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "Then the asymptotic standard error of $\\hat{\\psi}$ is:\n",
    "\n",
    "$$\\hat{\\text{se}}(\\hat{\\psi}) = \\sqrt{(\\hat{\\nabla} g)^T \\hat{J}_n (\\hat{\\nabla} g)}\n",
    "= \\sqrt{\\frac{p_1(1 - p_1)}{n_1} + \\frac{p_2(1 - p_2)}{n_2}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(d)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import norm, binom\n",
    "\n",
    "n = 200\n",
    "X1 = 160\n",
    "X2 = 148"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated psi: \t 0.060\n"
     ]
    }
   ],
   "source": [
    "p1_hat = X1 / n\n",
    "p2_hat = X2 / n\n",
    "psi_hat = p1_hat - p2_hat\n",
    "\n",
    "print(\"Estimated psi: \\t %.3f\" % psi_hat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "90% confidence interval (delta method): \t -0.009, 0.129\n"
     ]
    }
   ],
   "source": [
    "# Confidence using delta method\n",
    "\n",
    "z = norm.ppf(.95)\n",
    "\n",
    "se_delta = np.sqrt(p1_hat * (1 - p1_hat)/n + p2_hat * (1 - p2_hat) / n)\n",
    "confidence_delta = (psi_hat - z * se_delta, psi_hat + z * se_delta)\n",
    "\n",
    "print(\"90%% confidence interval (delta method): \\t %.3f, %.3f\" % confidence_delta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "90% confidence interval (parametric bootstrap): \t -0.010, 0.130\n"
     ]
    }
   ],
   "source": [
    "# Confidence using parametric bootstrap\n",
    "\n",
    "B = 1000\n",
    "xx1 = binom.rvs(n, p1_hat, size=B)\n",
    "xx2 = binom.rvs(n, p2_hat, size=B)\n",
    "t_boot = xx1 / n - xx2 / n\n",
    "\n",
    "se_bootstrap = t_boot.std()\n",
    "confidence_delta = (psi_hat - z * se_bootstrap, psi_hat + z * se_bootstrap)\n",
    "\n",
    "print(\"90%% confidence interval (parametric bootstrap): \\t %.3f, %.3f\" % confidence_delta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 10.13.8**.  Find the Fisher information matrix for Example 10.29:\n",
    "\n",
    "Let $X_1, \\dots, X_n \\sim N(\\mu, \\sigma^2)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**  The log likelihood is:\n",
    "\n",
    "$$\\ell_n = \\sum_i \\log f(x; (\\mu, \\sigma))\n",
    "= n \\left[ \\log \\left( \\frac{1}{\\sigma \\sqrt{2 \\pi}} \\right) + \\left( -\\frac{1}{2} \\left(\\frac{x - \\mu}{\\sigma} \\right)^2\\right) \\right]$$\n",
    "\n",
    "From this,\n",
    "\n",
    "$$H_{11} = \\frac{\\partial^2 \\ell_n}{\\partial \\mu^2} = -\\frac{n}{\\sigma^2}$$\n",
    "$$H_{22} = \\frac{\\partial^2 \\ell_n}{\\partial \\sigma^2} = -\\frac{n}{\\sigma^2} - \\frac{n}{\\sigma^2} = -\\frac{2n}{\\sigma^2}$$\n",
    "$$H_{12} = H_{21} = \\frac{\\partial^2 \\ell_n}{\\partial \\mu \\partial \\sigma} = 0$$\n",
    "\n",
    "So the Fisher Information Matrix is\n",
    "\n",
    "$$I(\\mu, \\sigma) = -\\begin{bmatrix}\n",
    "\\mathbb{E}[H_{11}] & \\mathbb{E}[H_{12}] \\\\\n",
    "\\mathbb{E}[H_{21}] & \\mathbb{E}[H_{22}]\n",
    "\\end{bmatrix} = \\begin{bmatrix}\n",
    "\\frac{n}{\\sigma^2} & 0 \\\\\n",
    "0 & \\frac{2n}{\\sigma^2}\n",
    "\\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercises 10.13.9 and 10.13.10**.  See final exercises from chapter 9."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}