{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Delta Method\n",
    "\n",
    "> CFO: What's our variation in churn this year?\n",
    "\n",
    "> Data scientist: Our $\\lambda$ value has been increasing, but $\\rho$ is staying the same, so we should see-\n",
    "\n",
    "> CFO: Our banana value is increasing?\n",
    "\n",
    "We want to connect our parameters to business logic _and_ carry over variance in estimates. \n",
    "\n",
    "Example: It's silly to present a point estimate without confidence intervals (CIs), since arguably the CIs contains more useful information than a point estimate. \n",
    "\n",
    "We'll start with asking: what is the CI for the survival function, $S(t; \\hat{\\theta})$? \n",
    "\n",
    "\n",
    "We will use the **DELTA METHOD** to do this (bolded because it's awesome):\n",
    "\n",
    "$$\\text{Var}(f(\\hat{\\theta})) \\approx \\text{grad}(f)(\\hat{\\theta}) \\cdot \\text{Var}(\\hat{\\theta}) \\cdot \\text{grad}(f)(\\hat{\\theta}) ^ T $$\n",
    "\n",
    "1. $f$ in our case is the survival function, $S$\n",
    "2. We know $\\text{Var}(\\hat{\\theta})$ (inverse of the Hessian)\n",
    "3. Do we need to compute $\\text{grad}(f)$ by hand? Heck no, use `autograd`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# seen all this...\n",
    "%matplotlib inline\n",
    "from autograd import numpy as np\n",
    "from autograd import elementwise_grad, value_and_grad, hessian\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "# N = 50 for this example\n",
    "T = (np.random.exponential(size=50)/1.5) ** 2.3\n",
    "E = np.random.binomial(1, 0.95, size=50)\n",
    "\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",
    "estimates_ = results.x\n",
    "H = hessian(negative_log_likelihood)(estimates_, T, E)\n",
    "variance_matrix_ = np.linalg.inv(H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.11163842, -0.24322003])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from autograd import grad\n",
    "\n",
    "def survival_function(params, t):\n",
    "    return np.exp(-cumulative_hazard(params, t))\n",
    "\n",
    "grad_sf = grad(survival_function)\n",
    "grad_sf(estimates_, 5.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def variance_at_t(t):\n",
    "    return grad_sf(estimates_, t) @ variance_matrix_ @ grad_sf(estimates_, t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.00036047601896770017"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "variance_at_t(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(.001, 10, 100)\n",
    "\n",
    "std_sf = np.sqrt(np.array([variance_at_t(_) for _ in t]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Estimated survival function with CIs (Delta method)')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZRc51nn8e9Tt9beu9Xa1ZY3eZFN4kVxQowTZwMnJHHYkpgJJJkQwzkYGMgBAjMnOIYZMgMHMgxmcUJWYhuTQBCJIbuzObYlx0ssybJkWVa3pG71vtf+zB/vre7b1dXd1a3qpaqezzl9uureW7feW8vvvvXe975XVBVjjDHVL7TeBTDGGFMZFujGGFMjLNCNMaZGWKAbY0yNsEA3xpgaYYFujDE1ou4CXURuEpGj612OUkTkZhHp2QDlOCQiN1dgPSdF5PULzEuIyL+LyKiI/PP5Ptcyy1WR7Vut5xeRh0TkV9awSAsSke+LyLUVXN+Cn4mNRkTuFJF/rNC63iMi3/Nvx0TkWRHZXIl1B1VNoPsfhGkRmQj8/XUZj1MRubRwX1W/q6qXr1IZPyUif7Ia615LqnqVqj60yk/z88BWYJOq/sJqPUmp92SNtm9BweevRGiISIuIfFRETvnfi+f9+53+/BWFqIi8BRhX1ScCZc2IyLj/95yI/LWIbF9huSsWmOdrLStTqpoCPgF8sNLrrppA971FVZsCf3esd4GqjYiE17sMvt3Ac6qaXe+CVDMRiQLfAK4CbgFagB8HBoEbznP1vwZ8tmjaP6lqM9AB/AywDXh8paFex+4F3i0isYquVVWr4g84Cbx+gXmXAt8GRoEB3IcO4DuAApPABPAO4Gagp2i9vws87S/3D7ia438A48DXgfbA8v8M9PrP9R3gKn/67UAGSPvP9e/+9B3AF4B+4AXgNwPrSgCfAoaBw345ehbYRgH+EjgHjAE/Aq725z0E/Epg2fcA3wvcV+DXgWN+Gf4W+POi9f8b8DvB19ov+zTQEVjuWv81jgCXAN/EhccA8Dmgban3DPiw/zpl/NfqfcCdwD8GlrnQL3c4sI1/DHzff1++CnQGlv8J4GFgBOj2X4OF3pOZcgEx4KPAGf/vo0DMn3cz0AN8wH/dzwLvXeD9eQ3wo8D9rwEHAve/C7yt6PW9peh1eKqcbS163l8B+oCmcr47LPBdKfGYqP/e7wpMm/Me+dM84CkCnyfgzcCT/nvxMPCS4rIssu3vBY74230C+NVFtus9/mv0l/5znQBe6U/v9t+zdweWjwF/DpzyX7O/w30HG/1tzftlmcB99u8EHgA+45fnELAvsL4r/fdqxJ/31sC8TcB+3Hf1Mf/9/F5R+Y8Br65oTlZyZav5x+KBfh/w33G/OOLATwTmKXBp4P7NzA/0R3AhvtP/EPwQF1xxXGD9UWD5/wo0MxsETwbmfQr4k8D9EPA48CH/C3Kx/6H7KX/+R3Bf9A6gC3iGhQP9p/x1teHC/Upgu84GwFKB/jX/eRLAq/wPvPjz2/0P9I4SAfBN4P2Bdf0Z8HeBcHiD/1psxu3gPlrme3YncwO8+P6FzA/054HL/G14CPiIP2837gt3G25Hswm4ptR7UmL77vLf/y3+NjwM/HHgs5L1l4kAbwKmCOzgA+tMAEmg01+2DziN+6wk/Nd3U4nnn7PdS21riee9H/h0ud8dFvmuFD3mKmBysfcsMP0u4FH/9rW479DLcWH/bv/5Y2Vu+0/jKgoCvNp/va9boIzv8d+f9/rP9Se4sL4b95n8Sf9z0eQv/5e4kO3w35d/B/60VC4Eypf033cP+FPgEX9eBDgO/CHuu/1a/7kuD7wvD+B2Flf7n4XiQN9PoIJXib9qa3L5ooiMBP7e70/P4L7UO1Q1qarfW+Z6/5+q9qnqaVzAPqqqT6hqEvhX3IcUAFX9hKqOq2sHuxN4qYi0LrDelwGbVfUuVU2r6gngY8A7/flvB/6nqg6pajfwV4uUMYP7EF6BC+Ijqnp2Gdv4p/7zTPvbqMBN/ryfB36gqmdKPO5eXFAiIuKX/V4AVT2uql9T1ZSq9gN/gfsSrpZPqupz/jY8AFzjT/9F4Ouqep+qZlR1UFWfLHOd/wW4S1XP+dvwYeCXAvMz/vyMqj6Iq73NOwbjl+kAbmd5Pa7W+n3gRuAVwDFVHazAthbbhPvlUK5yvyttuIAqxxlcSIL7VfT3qvqoquZU9dNACvcaLElVv6yqz6vzbdyvk5sWecgLqvpJVc0B/4SrGN3lfya/ivsVcKn/2b0d+G3/ezAO/C9mv4sL+Z6qPuiv/7PAS/3prwCacDvatKp+E/gScJuIeMDPAR9S1UlVfQb4dIl1j+Ne54qptkB/m6q2Bf4+5k//Pdwe/TG/B8F/XeZ6+wK3p0vcbwIQEU9EPuIfdBrD1TbA1cpK2Q3sCO6EcHv0rf78HbiacsGLCxXQ/8D8Na72cU5E7hGRlrK2zpl5HnXVg/vxgxoXiJ9b4HFfAH7cbyN9Fe5n6XcBRGSriNwvIqf91+MfWfi1qITewO0p/PcF9yV+foXr3MHc1/1Ff1rBoM5t5w8+b7Fv42p6r/JvP4Tbwb3av78cC21rsUFgOe3X5X5XhnEViHLsBIb827uBDxR95ruY+5ouSETeKCKPiMiQ/9g3sfhnqvi7iqqW+v5uBhpw7f2Fcv2nP30xxe9D3D8OtQPoVtV8YP6LuNdiMxBm6e92M665pmKqLdBLUtVeVX2/qu4AfhX4m2DPlgr6ReBWXBtgK65ZANwXBFytN6gbV4MI7oSaVfVN/vyzuA97wQWLPbmq/pWqXg/sxf0c/11/1iTuw1qwrdTDi+7fB/y8iOzG/Tz+wgLPOYyrJb0Dt/33+zsEcDUcBX5MVVuAdzH7WixXOduwkG7cz/RSlhpO9AwuhAou8KetRHGgf5ulA/18hzv9OvBTItJYzsLL+K4cx/0o27nY+kQkBLwFfyePey/+Z9FnvkFV7ytVnKJ1xXCfwz8HtqpqG/AgK/9MBQ3gwv2qQLlaVbWwo1zu+3AG6PK3v+ACXNNKP64paKnv9pW4X3IVUxOBLiK/ICK7/LvDuDensOfsw7VdV0Iz7ufjIC58/lfR/OLnegwYF5Hf9/tdeyJytYi8zJ//APAHItLul/83FnpiEXmZiLxcRCK48Esyu41PAj8rIg3+l/N9S22Iuq5oA8DHga+o6mI1hXuBX8Y1zdwbmN6Ma4IY9b/4v1viseV6EniViFzgN2H9wTIe+zng9SLydhEJi8gmESk0USz1/t8H/A8R2ex38/sQ7pfGSjyMa465AXhMVQ/hdhYvxx1fKKUPuLAoGJbjs7gQ/YKIXCEiIX/7/1BE3lS88BLflRmqmsbtLEo2ofmv85W4128brrkNXJPir/mfVRGRRhH5aREpVdsv3vYoru27H8iKyBtx7eDnza9Jfwz4SxHZ4m/DThH5qUBZNi3SfFrsUVyN/fdEJCLuvIK34Co8OeBfgDv97+Re3LGEGf73pQN3/KZiqi3Q/13m9kP/V3/6y4BHRWQCd6Dht/z2anDt3J/2f2a9/Tyf/zO4n06ncb1Sit+MfwD2+s/1Rf+NfTOu/fMFZgO08KH5sL++F3C14OIuYkEtuA/ksP+YQdwBSnAHe9K4D+WnWbj5pNi9uF8b9y6x3H5gD9CrqsEaxYeB63A9Jr6M+xCviKp+DdcG+jTu4O+XlvHYU7if5h/A/fR/ktm2zjnvSYmH/wlw0H/eH+EOiK/oXAJVnfQff8gPRIAfAC+q6rkFHlY4qWpQRH64gudM4d7DZ3EHvgu9KjpxoVNsse9Ksb9n7vEEgHf4jx31Hz8IXF84/qKqB4H345oHh3E1/fcssP452+63a/8mrqIzjPtFuH/BjV++3/fL84jfRPh1/OMhqvosbud0wv+sLNpE5L+/bwHeiPte/w3wy/56AO7ANfX04g7Mf7JoFb+IO5idqsB2zSj0cjDGmHlE5PvAHf4vOlMBftPSU8CrFtnRr2zdFujGGFMblmxyEZFPiMg5EXlmgfkiIn8lIsdF5GkRua7yxTTGGLOUctrQP4U7q2shb8S1r+7B9fP82/MvljHGmOVaMtBV9TvM9jEt5VbgM/6JAI8AbWLjOhhjzJqrxEBNO5nbgb7Hnzbv7DURuR1Xi6exsfH6K664ogJPb4wx9ePxxx8fUNWSJ0St6ch7qnoPcA/Avn379ODBg2v59MYYU/VEZMEzyivRD/00c8+I2uVPWzW5vPXMMcaYYpUI9P3AL/u9XV4BjC5z0Khly+TmndhmjDF1b8kmFxG5Dzc+Rae4K3r8EW7oSFT173BjLbwJdwbWFG4oy1WVtRq6McbMs2Sgq+ptS8wvXDxhTZzon+DAySF+4fouQqFKjNljjDG1odrGcuFrh/v4/S/8iOlMbr2LYowxG0rVBXpjzP2omEzbpSiNMSaoCgPdA2AyZTV0Y4wJqr5Aj/o19JTV0I0xJqj6At1vcplKWw3dGGOCqjbQrYZujDFzVV+gR/02dDsoaowxc1RdoDf4NfTxpAW6McYEVV2gN0ULgZ5Z55IYY8zGUnWB3uB3W7QaujHGzFV1gR7xQkQ8YcIC3Rhj5qi6QAdoiHqMWy8XY4yZo0oDPWzdFo0xpkhVBnoi6lmgG2NMkaoM9MZomEk7U9QYY+aoykBviHpM2YlFxhgzR1UGemMsbGO5GGNMkeoM9KjHtAW6McbMUZ2BHgvbFYuMMaZIVQZ6UyxsNXRjjClStYGezSvpbH69i2KMMRtGVQZ6Y7xwkQvr6WKMMQVVGejN/hC6E3ZykTHGzKjKQG/ya+hj0zaErjHGFFRloBcuQzc2bTV0Y4wpqM5A9y9yMWYXuTDGmBnVGeh2kQtjjJmnOgO9cBm6lNXQjTGmoDoDvdDLxWroxhgzo0oD3TW5WKAbY8ysqgz0RMRDgAk7scgYY2ZUZaCLCPGIx2TKxnMxxpiCqgx0sMvQGWNMseoN9IhnF7kwxpiA6g10uwydMcbMUVagi8gtInJURI6LyAdLzL9ARL4lIk+IyNMi8qbKF3Uud11Rq6EbY0zBkoEuIh5wN/BGYC9wm4jsLVrsfwAPqOq1wDuBv6l0QYs1RD2SdtUiY4yZUU4N/QbguKqeUNU0cD9wa9EyCrT4t1uBM5UrYmmNUbsMnTHGBJUT6DuB7sD9Hn9a0J3Au0SkB3gQ+I1SKxKR20XkoIgc7O/vX0FxZzXEPJKZPPm8ntd6jDGmVlTqoOhtwKdUdRfwJuCzIjJv3ap6j6ruU9V9mzdvPq8nbIyGSWZyZPJ2GTpjjIHyAv000BW4v8ufFvQ+4AEAVf0BEAc6K1HAhTTGPNLZPBm7rqgxxgDlBfoBYI+IXCQiUdxBz/1Fy5wCXgcgIlfiAv382lSW0BQLo8C4nVxkjDFAGYGuqlngDuArwBFcb5ZDInKXiLzVX+wDwPtF5CngPuA9qrqqjduzVy2yIXSNMQYgXM5Cqvog7mBncNqHArcPAzdWtmiLa/IDfdQuQ2eMMUAVnynaZDV0Y4yZo2oDvTkeAWB4Kr3OJTHGmI2hagN9tsnFaujGGANVHOjNcQt0Y4wJqtpAb/IDfTyZYZU71BhjTFWo2kBvjrk29Ol0npSdXGSMMdUb6A3+haJT2ZyNumiMMVRxoEe8EBFPSGbyNuqiMcZQxYEOEAt7pLI5pu1CF8YYU92Bnoh6VkM3xhhfdQd6JETS2tCNMQao8kCPR9xl6OzaosYYU+WB3hANk8rkSdmVi4wxproDPRHxZvqgJ7NWSzfG1LeqDvSGqDfTfm49XYwx9a7qA71QQ7d2dGNMvavuQI+FZ2ro1tPFGFPvqjrQG6Me2bySzVlfdGOMqe5A98dET2bz1oZujKl7VR3ohYtcTKdzVkM3xtS9qg70ba1xAIYm02RzSiZnw+gaY+pXVQf6rrYEAAOTKQCrpRtj6lpVB/r2tgQCDE64C0UnrR3dGFPHqjrQExGPtoYIAxOuhm590Y0x9ayqAz0UEjY1xmZq6BOp7DqXyBhj1k9VB7oXEjY1RRn029DHkpl1LpExxqyf6g50ETqbYgxNpsnllel0jpyNumiMqVNVHeihEGxqipJXGJlKo2rNLsaY+lXVge6J0NkYA2DA2tGNMXWuugPdb0OH2b7o49aOboypU1Ud6CJCZ3N0Tl/0iaTV0I0x9amqAx0gFp7bF308lUXVDowaY+pP1Qe6V9QXPZdTkhkb08UYU3+qP9Blbl90gPGUtaMbY+pPWYEuIreIyFEROS4iH1xgmbeLyGEROSQi91a2mAsLheb2RQcYt3Z0Y0wdCi+1gIh4wN3AG4Ae4ICI7FfVw4Fl9gB/ANyoqsMismW1Clys0NOl0Bd9U1PMDowaY+pSOTX0G4DjqnpCVdPA/cCtRcu8H7hbVYcBVPVcZYu5sFCJvuhWQzfG1KNyAn0n0B243+NPC7oMuExEvi8ij4jILaVWJCK3i8hBETnY39+/shIXKdUXPZnJ2cUujDF1p1IHRcPAHuBm4DbgYyLSVryQqt6jqvtUdd/mzZsr8sSeCB2Nc/uig/VHN8bUn3IC/TTQFbi/y58W1APsV9WMqr4APIcL+FUX9oSIF5rTFx2s2cUYU3/KCfQDwB4RuUhEosA7gf1Fy3wRVztHRDpxTTAnKljOBcUjHsCcvuhgQ+kaY+rPkoGuqlngDuArwBHgAVU9JCJ3ichb/cW+AgyKyGHgW8DvqurgahU6KFEI9KK+6MNT6YUeYowxNWnJbosAqvog8GDRtA8FbivwO/7fmopH3D6psynGgZND5PKKFxJSmTzT6RyJqLfWRTLGmHVR9WeKxgM19LzOrZmPTFst3RhTP6o+0GPhECKwpdn1Re8dTc7MG560dnRjTP2o+kAXEeIRj672BgC6h6dm5lkN3RhTT6o+0MG1ozfGwnQ0Rukemp6ZPpXKkcrm1rFkxhizdmoi0GNh147e1Z6YU0MHGJ2yZhdjTH2oiUAvHBjt6migdyxJOjt72v+wBboxpk7URKAXuiZ2tTegCqdHZptdRqw/ujGmTtREoMfDbjO6OhLA3AOjE6msDdRljKkLtRHofpNLZ1OMWDhE99BsoKvC6LQ1uxhjal9NBXpIhF3tCXqGp+fMt2YXY0w9qIlA90JCxG92uaCjge7hKfKqM/ODg3YZY0ytqolAh9l29F3tDSQz+TkhPp7MksxYf3RjTG2rmUCf6elS4sAoMGesdGOMqUU1E+iFdvSdbQlEmHNgFGavN2qMMbWqdgLdP1s0FvbY2hynu+jA6PBkmlxeSz3UGGNqQu0EenR2U7o6EvQUNbnk8moXvTDG1LTaCfTI7IUsutobGJhIM5Wee11Ra0c3xtSy2gn0cCDQO9xQuqeK29HHrYZujKldNRPo0XAILyQAXNTZiADHzk3MWSaZyTFuF482xtSomgl0gJh/fdGmWJid7Qme6x2ft4z1djHG1KqaCvRgO/plW5t5vn+SbNHAXP3j1o5ujKlNNRXoiUCgX761mXQuz8nBue3oY9OZeQdLjTGmFtRUoDdGwzO3L9vaBMBzffObXc4GLiRtjDG1oqYCvTkeDtyOsKM1XjLQ+yzQjTE1qKYCvSkQ6ODa0Y+dm5h3huhUOmfXGjXG1JyaCvSIF5p3YDSVzc/rjw7QO2a1dGNMbampQIe5tfTLtzUDpdvR+8aSqNrYLsaY2lFzgR5sR29NRNjaEuNoiUBPZ/MMTlqfdGNM7ai9QI/NbUe/fGszx/omyJcYabHXDo4aY2pIzQV68YHRPVubmc7k5l1nFODceJJ0Nj9vujHGVKOaC/RExMPzZOb+FX47+jNnRuctm8/D2dH5QW+MMdWo5gJdROY0u7Q3RNm9qYEnu0dKLt8zPG0HR40xNaHmAh3mN7tc09XGCwOTjE7P73s+nc7ZwVFjTE2ozUCPzQ90BZ5apJZujDHVriYDvTkWmXN/V1uCzqbogs0uA+MpptO5tSiaMcasmrICXURuEZGjInJcRD64yHI/JyIqIvsqV8Tla4qHkdnjoogI13S1cfjsGMlM6eAuvgapMcZUmyUDXUQ84G7gjcBe4DYR2VtiuWbgt4BHK13I5fJCQiLqzZl2TVcb2bxy6MxYycecHpmeN3a6McZUk3Jq6DcAx1X1hKqmgfuBW0ss98fA/wY2xNk6xc0ue7Y00xD1Fmx2yebU2tKNMVWtnEDfCXQH7vf402aIyHVAl6p+ebEVicjtInJQRA729/cvu7DL0VzU08ULCS/Z1crTPSPzRl8sODU0teA8Y4zZ6M77oKiIhIC/AD6w1LKqeo+q7lPVfZs3bz7fp15Ue0N03rRrutqYTOc4XnTx6IJ0Ns+ZEaulG2OqUzmBfhroCtzf5U8raAauBh4SkZPAK4D9631gtDkennPGKMDVO1qJeiEeOzm04ONODk6WHPfFGGM2unIC/QCwR0QuEpEo8E5gf2Gmqo6qaqeqXqiqFwKPAG9V1YOrUuIyhUJCW2JuO3o84nHd7jYee2FowTFcUpk8Z22sdGNMFVoy0FU1C9wBfAU4AjygqodE5C4ReetqF/B8lGp2ufGSTqYzuQUPjgKcHLBaujGm+oSXXgRU9UHgwaJpH1pg2ZvPv1iV0d44P9Av39ZMR2OU7z8/wA0XdZR83HQ6x+mRabo6Gla7iMYYUzE1eaZoQUuJdvSQCK+8ZBOHz44xPLXwGC4nBiatX7oxpqrUdKCLSMlml1desglV+MHzgws+NpPNc3LQzh41xlSPmg50gPaGyLxpW5rj7NnSxMPPDy46dO6pockFhwowxpiNpvYDvUQ7OriDo71jSU4MTC742Hwenu8v3WfdGGM2mpoP9ObY/HZ0gH0XthMLh3jo6OJnrJ4dSZYcR90YYzaamg/0hdrR4xGPGy/t5LGTQ4wscnAU4MjZMbuqkTFmw6v5QAfoKBHoAK+/cgv5vPKtJWrpE8ks3UM2JIAxZmOri0Df3BwrOX1Lc5xrLmjjoaPnSGUXP/j5fP+EHSA1xmxodRHoiag37zqjBW+4ciuT6RyPnFh4fBeAXF452ju+GsUzxpiKqItAB9iyQC19z5YmLtzUwNcO95Ffop28fzxFn43zYozZoOon0FviJaeLCG/Yu5XesSQ/Oj265HqOLHIZO2OMWU91E+hNsTANMa/kvOt3t9PRGOXLT59dsjdLNqcctl4vxpgNqG4CHdxB0FLCoRBvfsl2TgxM8lTP0rX0oYm0Xa7OGLPh1Fegt5RuRwc3vsuW5hhffPL0km3pAMfOjTORylayeMYYc17qKtBb4hHikdLNLuFQiFtfuoOe4WkOnhxecl35PDzdPWIjMhpjNoy6CnRYvJb+sos62NmW4N+eOl3WxaKn0jkOnx2rZPGMMWbF6i7Qty7Q2wXcWOlvu2YHfWMpHn5+oKz1nRtLccqG2TXGbAB1F+iticiCJxkBXNPVxsWdjXzxyTNMp8vrnnjs3DjDk4uPB2OMMaut7gIdYGdbYsF5IsIv3nABY9MZ9j91pqz1qcJTPSNMpe0gqTFm/dRloG9rjeOF5g+pW3BhZyM37enkG8/2cbrM7onZnPLkqRHSWTtIaoxZH3UZ6BEvtOjBUYCfvXYXiYjH5x57seyTiKbSOZ7uGSFfxgFVY4yptLoMdIBdbQ2Lzm+Kh/nZ63bxXN8Ej51cfOCuoJGpDIfO2Jmkxpi1V7eB3toQoTG28MFRgJsu7eTCTQ3cf6CbsWVctahvLMmRszYyozFmbdVtoAPsal/44ChAKCS898aLmE7n+Oyj5Te9AJwZmea5Pgt1Y8zaqetA394aJ1zieqNBO9sSvO2anTxxamTJMdOLnRqc4vg5u8i0MWZt1HWgh70QXR2Lt6UD/OTerVy6uYl7HzvF0DL7m58cmOT4OaupG2NWX10HOkBXe8OiXRih0PRyITlVPvn9F5bdi+XkwJQ1vxhjVl3dB3o0HGLnEm3p4IYMuO1lXRzpHWf/0+WdcBR0anCKIzaOujFmFdV9oANc0NFAqIxX4icu7eTGSzbxpafPlnV1o2Knh6d5ume0rIG/jDFmuSzQgXjEY1vL0rV0EeG/vHw3Xe0JPv7dEwxOpJb9XP3jKZ44NUzGht01xlSYBbrvos7Gsmrp0XCIX3v1JeQV7n7oeVIruL7oyFSGAy8MMWkXyDDGVJAFui8R9djVvnSPF3Dt6e+/6SK6h6f42HeXf5AU3DABB04OraiWb4wxpVigB1zU2UgkXN5L8pJdbdz2sgt4smeEBx7vXtHzZXPKk90jNp66MaYiLNADIl6Iizsby17+tVds4Q1XbuXrR87x9SN9K3pOVXiub5yne0asXd0Yc17KCnQRuUVEjorIcRH5YIn5vyMih0XkaRH5hojsrnxR18bOtgQNsdLXHS3lF67fxbUXtPFPB7r5wfODK37ec2MpDrwwxHiy/DFjjDEmaMlAFxEPuBt4I7AXuE1E9hYt9gSwT1VfAnwe+D+VLuhaCYWEPVual7X87TddzBXbmvnEwy9wYBkjMxYrtKufGpyy/urGmGUrp4Z+A3BcVU+oahq4H7g1uICqfktVCw3BjwC7KlvMtbW5ObbkeOlBES/EHa+5lEs3N/Hx777AE6eGV/zc+bxrgnmie4TkCnrQGGPqVzmBvhMIHvXr8act5H3Af5SaISK3i8hBETnY399ffinXweXbmpccuCsoFvH4zdfu4YJNDfzdd05w8MWV19QBhibSPHJikN7R5HmtxxhTPyp6UFRE3gXsA/6s1HxVvUdV96nqvs2bN1fyqSsuFva4fFv5TS/guj7+9uv3cOGmBv7+Oyf43rGB8ypDNqc8c3qUJ622bowpQzmBfhroCtzf5U+bQ0ReD/x34K2qWhOdq7e3JuhsLr/pBaAhGuZ3Xn8Ze7e18KkfnOSrh3vPuxwD4ykeOTFI95C1rRtjFlZOoB8A9ojIRSISBd4J7A8uICLXAn+PC/NzlS/m+rlimU0v4Jpf7njtpVy/u50HDvZw32Onzvs6o9mccrR3nAMnhxmznjDGmBKWDHRVzQJ3AF8BjgAPqOohEblLRN7qL/ZnQBPwzyLypIjsX2B1VSce8bhqR+uyHxfxQvzqTRfzhiu38o1nz/H/vnWc6fT5N5uMTYusyasAABDUSURBVGd47MQQh8+MkcpaM4wxZpas10/4ffv26cGDB9fluVfi+LlxTg6s7IzOh46e497HTrG9NcGvv+YStjTHK1ImzxMu7mykq72B0BJjuhtjaoOIPK6q+0rNszNFy3RxZxNtDZEVPfbmy7fw3153GcNTaf74S0fOq1tjUC6nHOub4OHnBzkzMm3t68bUOQv0MoVCwtU7W4mWOdZLsb07WvjQm/eypSXG3Q89z+cf7yGbr8yp/slMjsNnxnjkxBDnxq2bozH1ygJ9GeIRj5d2tS15ybqFdDbF+OAtV/Dqyzbzn4d6+T//eZS+scoF8GQqy9PdozxyYtCC3Zg6ZIG+TK2JCFftbFnx4yNeiF96xW5+9VUX0zuW5MNfOsx3nuuvaHPJRHI22HtHk9YUY0ydsEBfgS3NcS7buryTjoq97MIO7nzLVVzS2chnHnmRj379GAMVHht9IpnlmdOjPPy868Nul74zprZZL5fzcD49Xwryqjx0tJ8v/LAHBX7mmp287ootq9JrJewJu9ob2NWeIB4pf0RJY8zGsVgvFwv08/Rc33hFLlAxOJHis4++yDOnx7igo4Hbbuha1qiPyyHifmXsak/Q3hhdlecwxqwOC/RVdrR3nO6h8w91VeXAyWH++fFuhqcy/PjFm/i563bS1rB6odsYC7OzLcG21viKe/AYY9aOBfoaqFRNHSCVyfHlZ87y1UN9hER4w96t3HLVNhLR1WsmCYVcrX17a5yOxigidqKSMRuRBfoaeXFwkmN9ExVbX/94in994jSPnRyiKRbmTT+2jZsv27LqNelYJMT21jjbWhM0xcKr+lzGmOWxQF9DZ0enOXxmjEq+rCcHJvnCD3s40jtOayLCLVdt49WXbV6TJpKmeJhtLXG2tsRX9ReCMaY8FuhrbGgyzY9Oj5LJVvaiz8/1jbP/qTM82ztOczzM667Ywmsu30LjGtWiWxIRtrbE2NJs4W7MerFAXwfT6RxP9YwwkcxWfN1He8f5j2fO8syZMWLhEDft6eR1V2xl8zLHbj8fTfEwW5pjdDbHaImvbIwbY8zyWaCvk1xeOXJ2bNUuI3dqaIqvHOrl4Mlh8qq8dFcbr7liM1dubyG0hgc14xGPzuYomxpjdDRGVzw0gjFmaRbo6+zs6DTP9o6Ty63Oaz08leaho/18+7l+JlJZNjfHeNWeTl55SSetibWtPYdC0NYQpbMxRkdT1A6qGlNhFugbwFQ6y6EzY4xOrd7VhjK5PD98cZhvH+vnub4JPBGu3tnCKy/p5KW7Wgl7a9/PPBYJ0d4QZVNTlPaGqJ2hasx5skDfIFSV7qFpnu+fWPVxVc6OTvP944P84MQgo9MZGqIe+3a3c8NFHVy2pXndLojREPVoa4jS0RilrSFiAW/MMlmgbzDT6RxHescYmkiv+nMV2vF/cGKQJ7tHSGXztCYiXHdBG9dd0M5lW5vXtc07HvFoa4jQ1hChNRGhKRa2k5qMWYQF+gZ1bizJsXMTFbnWaDlS2RxPdY/y+IvD/Oj0KOlcnqZYmB/b2co1XW1ctaNl3WvMnie0JiK0xF3AtyTCxMJWizemwAJ9A8vllVNDU5wcnFy1g6alpDI5njkzxg9PuXCfSucIh4Q9W5u4ekcrP7azle2t8Q1RW45HPFoSYVriEZrjYVoSESLrcDzAmI3AAr0KpLN5XhycpHt4igpdma5subxy/NwET/WM8MyZUc6MuG6W7Q0RrtzewpXbW7hiWzPtqzhI2HIloh5NsTDN8TBNcRf26/3rwpi1YIFeRZKZHKeGpjg9PL1uF6QYnEhx6MwYh8+O8WzvOBMpd3LU1uYYl29rZs+WZi7d0kRn08YaxCvsCU0xF/CNURf2jbGw1eZNTbFAr0LpbJ7u4Sl6hqcrPoTAcuRV6R6a4mjfOM/1TnC0b5zpjGvzb0tEuGRLExd3NnLJ5iZ2b2rYkOEZDYdojIVpjHk0Rl3IN0Q9q9GbqmSBXsXyeaV3LMmpoalVGUZg2eVR5fTINMf7Jjh2boITAxMM+L11PBF2tie4qLOR3Zsa2N3RwM62xLr0fy+H5wmNURfu7i9Mwr+9EXdMxoAFes0YncpwemSavrHkhro+6MhUmhMDk5wcmOSFgUlODk7N1OLDIWFHW4Ku9gRdHe7ydzvbEjRv8PFfIuEQiYg3U5NviHokIh6JqEcsHNpQTU2mvlig15hsLs+58RS9Y8k16cu+XHlV+sdTvDg4xYtDk/QMTXNqeIrxwC+M1kSEHW1xdrYl2NGaYHtbnO0tCZriG3+ogFAI4mGPeNRz/yMhEjO3XeCv14lbpvYtFugb/9tj5gl7IXa0JdjRliCZydE/nqJvLMnIKg4rsBwhEbb6Y6jfcFEH4M6SHZ3O0DM8zemR2b/vHBsgHThG0BRz469va42ztSXG1pY4W5pjbG6ObZj+6Pk8TKVzTC1y/kAsEiLmh30h5Av/Y2ELfbM6rIZeQ1JZF+4DE2mGJlNr3v1xJfKqDE6kOTs6Te9Ykt7RJGdHk5wbTzE6PXcH1ZaIsNkP986mGJ1NUTqbYmxqdOPEVFtARsIhYuEQUf9/IegL06LhEFEvtGGPQZj1YTX0OhELe+xqb2BXewO5vDI8lWZoMs3ARIqp1NqcjbpcIZGZkH5J0bzpdI6+MRfu58bd/4GJFEfOjjEylUHnrIeZQcA6GqN0+OPFtPth396w8YYVyGTzZfVg8kIyE/oRb27YRwPTIp4Q9ax9v55ZoNcoLyR+LTbGZVubSWZyjExlGJpMMzKVXrS5YKNIRD0u7Gzkws7GefMyufzMzmpgIs3gZIrBCbcDO9Y3wchUhlzRr89wSGbGjGlriLr/CXffDTMQoSUepjke2VBjuufyynQ6V/YQEWE/2Ath7/7E/feDPxJyt8Mht2y1/boxpVmg14l4xGNbq8e21jjgmmdGpzOMTmUYS2YYm85uqJ4zS4l4oZl2+lLyeWUsmWFoKs3wZIbhqTQjUxlGpt3/08Pu2q+F3jhBAjTGwrQkwjTH3HADzX7QNwfOTm2Kzf5tpGaRbE7J5hZv4y/mhYSwJ4RDs+EfLvwPzd73Qm5nUJjnhYRwSOxXwQZhgV6nYmGPLc0eW5pdIKoqk+kc48kM48ksY9MZJlJZsms4vkwlhUJCW0OUtoYodC68XCqbY2w663Zu04WdW4axZJaxZIaJZJZuv4fOogdBwyGaYu6kpcaoN3Py0sz/mf7u7n+hv3si6hEOrf/OIJdXcnklxcoOvHieC/bCDsHzb3shmXvfX66wIwguV/gzK2eBbgAQkZna5vbW2enJTI7xZJbJVJbJdJbJVI7JdHZNBxJbTbGwx+Zmr6zrsWbzeSZTbqc3kcq6v6T7P5nKuf9p91oND08zmXY7gaV++UQ91+0xEZntAulu+33fIx4xv7dMPOIRL/SYiYTmdJWMRULrtnPI5ZRcbuU7hAIR5oR7MOzDISEkMme+J4Lnuf+hkDu5LRwKudvB5UXqolnJAt0sqhAixYGXzLg23alMjmk/uKbSOZKZXNXW6pcSDoVoTYSWdVk/VSWVzTPlt4EXQn46nWMqnWW68Dqmc0xn3Os3nXG/Ggrzktkc5XZGCwcOoMbCXuD27P3CtGjRAdbInAOtMmdaoR2+MG21rlmrWmgyKmxw5Y71iLhfbp4f8rM7B1ehCU4v7BxCheX8aYXHiPjzZ5aZ3YEUz1tLFuhmRQpB315iXiaXnwmnZDpPKpsjmZn9n87lqqJLZSWIyMxrxfxju2VRVdLZPMls3r2mmbmvZ/B/KpsnFbydzZPO5plIZRmaTJPOzU5LZ/OsdNcbbFcvhP289vdQaOZ2OCSEA8uFA000wfnhQFv+nKaZUtOK/sKhECFhwfZ8Vf+XxIq3emUKO4LC35YW11FhNZQV6CJyC/B/AQ/4uKp+pGh+DPgMcD0wCLxDVU9WtqimWhR6VrQscnp/OutCJ53Nu5DJ5MkUwibnuvNlclpX4b8QESEW8YhFvIpe9FtVyebdL4hMbjbk07nZ++49KLwfhXlKxl+mcDvrL5fN5cnm3Q5oMpcjm8uTyevM9Gxh+VU8AB9sjgl7s80tc5ppQoXaOPOml6rFB2vghfnFtfiQEFjevx+s4fs7m47GKPGwxwWbGiq+7UsGuoh4wN3AG4Ae4ICI7FfVw4HF3gcMq+qlIvJO4H8D76h4aU3NKPzUL0cuPzdACiGRybpgyOTyM8u4/+4AXzbv7q/TuXMbnojM1KjXmqqSU51pXsnm/cAvhH/OzQ++rzPvrbr2+sL7W1hPYfngX9b/nw/OC9zO+/fT2by/HDPL5Es8Nq/M3FdlXtfYcjXFw7xr0+4Kv6rl1dBvAI6r6gkAEbkfuBUIBvqtwJ3+7c8Dfy0iout1GqqpKa4GtfLhbgvhns8zGwLzvthuXvALW/gyF77kqrNf6kIgFZZX1br/JbEcIkJYhHAI2NjjtC0pXwj94s9HiR1AXpX2hig/sWeRrlfnoZxA3wl0B+73AC9faBlVzYrIKLAJGAguJCK3A7f7dydE5OhKCo3riDaw5FK1xba5Ptg214fz2eYFq/ZrelBUVe8B7jnf9YjIwYXGMqhVts31wba5PqzWNpfTeHYa6Arc3+VPK7mMiISBVtzBUWOMMWuknEA/AOwRkYtEJAq8E9hftMx+4N3+7Z8Hvmnt58YYs7aWbHLx28TvAL6C67b4CVU9JCJ3AQdVdT/wD8BnReQ4MIQL/dV03s02Vci2uT7YNteHVdnmdRsP3RhjTGWt/6hAxhhjKsIC3RhjakTVBbqI3CIiR0XkuIh8cL3Ls9pEpEtEviUih0XkkIj81nqXaS2IiCciT4jIl9a7LGtBRNpE5PMi8qyIHBGRH1/vMq02Eflt/zP9jIjcJyKlB7evYiLyCRE5JyLPBKZ1iMjXROSY/7/UkEgrUlWBHhiG4I3AXuA2Edm7vqVadVngA6q6F3gF8Ot1sM0AvwUcWe9CrKH/C/ynql4BvJQa33YR2Qn8JrBPVa/GdbhY7c4U6+FTwC1F0z4IfENV9wDf8O9XRFUFOoFhCFQ1DRSGIahZqnpWVX/o3x7HfdF3rm+pVpeI7AJ+Gvj4epdlLYhIK/AqXG8xVDWtqiPrW6o1EQYS/rkrDcCZdS5Pxanqd3A9/4JuBT7t3/408LZKPV+1BXqpYQhqOtyCRORC4Frg0fUtyar7KPB7cJ5XS6geFwH9wCf9ZqaPi8gKB9utDqp6Gvhz4BRwFhhV1a+ub6nWzFZVPevf7gW2VmrF1RbodUtEmoAvAP9NVcfWuzyrRUTeDJxT1cfXuyxrKAxcB/ytql4LTFLBn+Ebkd9ufCtuZ7YDaBSRd61vqdaefwJmxfqOV1uglzMMQc0RkQguzD+nqv+y3uVZZTcCbxWRk7gmtdeKyD+ub5FWXQ/Qo6qFX16fxwV8LXs98IKq9qtqBvgX4JXrXKa10ici2wH8/+cqteJqC/RyhiGoKeIuv/IPwBFV/Yv1Ls9qU9U/UNVdqnoh7v39pqrWdM1NVXuBbhG53J/0OuYOT12LTgGvEJEG/zP+Omr8QHBAcKiUdwP/VqkVV9Ul6BYahmCdi7XabgR+CfiRiDzpT/tDVX1wHctkKu83gM/5FZUTwHvXuTyrSlUfFZHPAz/E9eR6ghocAkBE7gNuBjpFpAf4I+AjwAMi8j7gReDtFXs+O/XfGGNqQ7U1uRhjjFmABboxxtQIC3RjjKkRFujGGFMjLNCNMaZGWKAbY0yNsEA3xpga8f8B8FPmG1MTLfoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t, survival_function(estimates_, t))\n",
    "plt.fill_between(t, \n",
    "                 y1=survival_function(estimates_, t) + 1.65 * std_sf, \n",
    "                 y2=survival_function(estimates_, t) - 1.65 * std_sf,\n",
    "                 alpha=0.3\n",
    "                )\n",
    "plt.ylim(0, 1)\n",
    "plt.title(\"Estimated survival function with CIs (Delta method)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we will explore a subscription service LTV example. Move to Part 7! "
   ]
  }
 ],
 "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
}