{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hessians to compute variances\n",
    "\n",
    " - We have our point estimates, but now we want variances of the estimates\n",
    " - We can approximate them (under some asymptotic conditions) using the _Hessian at the (negative) MLE minimum_\n",
    " - Think of the the Hessian as a matrix representing the curvature at the minimum point. If we move slightly away from that minimum, how much change do we see in the likelihood?\n",
    "   - high curvature (think of a tight peak) => low variance of estimate.\n",
    "   - low curvature (think of a broad hill) => high variance of estimate.\n",
    "\n",
    "To extract the parameter estimates, we invert the Hessian at the MLE, $\\hat{\\theta}$, and take the diagonal:\n",
    "\n",
    "$$\\text{Var}(\\hat{\\theta}_i) = \\text{diag}(H(\\hat{\\theta})^{-1})_i $$\n",
    "\n",
    "The Hessian is the matrix of all second derivatives. `autograd` has this built in:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from autograd import numpy as np\n",
    "from autograd import elementwise_grad, value_and_grad, hessian\n",
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = (np.random.exponential(size=1000)/1.5) ** 2.3\n",
    "E = np.random.binomial(1, 0.95, size=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      fun: 249.99413038558657\n",
      " hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>\n",
      "      jac: array([-0.00046244,  0.00022406])\n",
      "  message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n",
      "     nfev: 11\n",
      "      nit: 7\n",
      "   status: 0\n",
      "  success: True\n",
      "        x: array([0.46459167, 0.44953271])\n"
     ]
    }
   ],
   "source": [
    "# seen all this...\n",
    "def cumulative_hazard(params, t):\n",
    "    lambda_, rho_ = params\n",
    "    return (t / lambda_) ** rho_\n",
    "\n",
    "hazard = elementwise_grad(cumulative_hazard, argnum=1)\n",
    "\n",
    "def log_hazard(params, t):\n",
    "    return np.log(hazard(params, t))\n",
    "\n",
    "def log_likelihood(params, t, e):\n",
    "    return np.sum(e * log_hazard(params, t)) - np.sum(cumulative_hazard(params, t))\n",
    "\n",
    "def negative_log_likelihood(params, t, e):\n",
    "    return -log_likelihood(params, t, e)\n",
    "\n",
    "from autograd import value_and_grad\n",
    "\n",
    "results = minimize(\n",
    "        value_and_grad(negative_log_likelihood), \n",
    "        x0 = np.array([1.0, 1.0]),\n",
    "        method=None, \n",
    "        args=(T, E),\n",
    "        jac=True,\n",
    "        bounds=((0.00001, None), (0.00001, None)))\n",
    "\n",
    "print(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 887.54177234, -758.67183105],\n",
       "       [-758.67183105, 8295.5360413 ]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "estimates_ = results.x\n",
    "\n",
    "# Note: this will produce a matrix\n",
    "hessian(negative_log_likelihood)(estimates_, T, E)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.00122226 0.00013077]\n"
     ]
    }
   ],
   "source": [
    "H = hessian(negative_log_likelihood)(estimates_, T, E)\n",
    "variance_ = np.diag(np.linalg.inv(H))\n",
    "print(variance_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.03496082 0.01143546]\n"
     ]
    }
   ],
   "source": [
    "std_error = np.sqrt(variance_)\n",
    "print(std_error)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lambda_  0.465 (±0.03)\n",
      "rho_     0.450 (±0.01)\n"
     ]
    }
   ],
   "source": [
    "print(\"lambda_  %.3f (±%.2f)\" % (estimates_[0], std_error[0]))\n",
    "print(\"rho_     %.3f (±%.2f)\" % (estimates_[1], std_error[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From here you can construct confidence intervals (CIs), and such. Let's move to Part 6, which is where you connect these abstract parameters to business logic. "
   ]
  }
 ],
 "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
}