{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### Lifetime value example \n", "\n", "\n", "Suppose we have a subscription business that has monthly churn, and we'd like to know an estimate of LTV (lifetime value) and build confidence intervals for it. \n", "\n", "\n", "Subscription businesses have a very predictable churn profile (it looks piecewise) but the rates are unknown. \n", "\n", "\n", "\n", "\n", "\n", "We'll use a piecewise-constant hazard model with known breakpoints, $\\tau$.\n", "$$ \n", "h(t) = \\begin{cases}\n", " \\lambda_0 & \\text{if $t \\le \\tau_0$} \\\\\n", " \\lambda_1 & \\text{if $\\tau_0 < t \\le \\tau_1$} \\\\\n", " \\lambda_2 & \\text{if $\\tau_1 < t \\le \\tau_2$} \\\\\n", " ...\n", " \\end{cases}\n", "$$\n", "\n" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "%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", "df = pd.read_csv(\"../churn_data.csv\")\n", "T = df['T'].values\n", "E = df['E'].values" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " fun: 6820.877953596707\n", " hess_inv: <10x10 LbfgsInvHessProduct with dtype=float64>\n", " jac: array([ 0.57181706, -0.16753806, -0.06621519, -0.19042741, -0.09172935,\n", " 0.04749549, 0.1030577 , -0.21740178, -0.08800897, -0.21890358])\n", " message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n", " nfev: 67\n", " nit: 46\n", " status: 0\n", " success: True\n", " x: array([0.00103406, 0.03326294, 0.00104654, 0.01516408, 0.0005534 ,\n", " 0.01317414, 0.00099483, 0.01216659, 0.0012531 , 0.00689811])\n" ] } ], "source": [ "breakpoints = np.array([28, 33, 58, 63, 88, 93, 117, 122, 148, 153])\n", "\n", "def cumulative_hazard(params, times):\n", " # this is NumPy black magic to get piecewise hazards, let's chat after. \n", " times = np.atleast_1d(times)\n", " n = times.shape[0]\n", " times = times.reshape((n, 1))\n", " M = np.minimum(np.tile(breakpoints, (n, 1)), times)\n", " M = np.hstack([M[:, tuple([0])], np.diff(M, axis=1)])\n", " return np.dot(M, params)\n", "\n", "hazard = elementwise_grad(cumulative_hazard, argnum=1)\n", "\n", "def survival_function(params, t):\n", " return np.exp(-cumulative_hazard(params, t))\n", "\n", "def log_hazard(params, t):\n", " return np.log(np.clip(hazard(params, t), 1e-25, np.inf))\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.ones(len(breakpoints)),\n", " method=None, \n", " args=(T, E),\n", " jac=True,\n", " bounds=[(0.0001, None)] * len(breakpoints)\n", ")\n", "\n", "print(results)\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": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEXCAYAAAC9A7+nAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxddZ3/8dcna7MnzdIlSdeElrJ0IYCALK6AIkVEp9VxZFxwfj9BZXBG0Pk5ijPOjPpzR/2BKIpsiogFUURUQESktGXrvu9b2jRbm/Xz++OcpLdp0ty0Se+5N+/n43Efufecc+/53HOT9/3me875HnN3REQk+aUlugARERkeCnQRkRShQBcRSREKdBGRFKFAFxFJEQp0EZEUoUBPEWZ2oZmtSnQd/TGzS8xsawTqeM3MLhmG19loZm8eYF6OmT1iZgfM7Ocnuq4h1jUs72+4mNlvzOwDia5jNMlIdAGjnZltBMYBXTGT73L36wd5ngO17r4WwN2fAWaMUI13AVvd/d9G4vVPFnc/7SSs5hqCz7PU3TtHaiX9fSYn6f3Fzd0vT3QNo40CPRre4e6/T3QRyczMMkYyQIdgMrA6IrXIaOPuuiXwBmwE3jzAvBrgKeAAsBd4IJz+NOBAC9AM/B1wCUGLLfZ1/wV4OVzuToKW42+AJuD3QEnM8j8Hdobreho4LZx+HdABtIfreiScPhH4BbAH2AB8POa1coC7gP3A8rCOrQO8RwO+DuwGGoFXgNPDeX8CPhyz7LXAn2MeO/AxYE1Yw/eAr/Z5/V8B/xy7rcPaDwJjY5abG27jTGA68AegPpx2D1A82GcGfCHcTh3htvoQ8HngpzHLTAnrzoh5j18Eng0/l98BZTHLvx74C9AAbAm3wUCfSW9dQDbwDWB7ePsGkB3OuwTYCtwUbvcdwD/G+zsa+56AMcBPw23VALwAjOv7+fV8dsBXw9+LDcDlMa85leD3rud387bY7aZbfDf1oUfbFwn+wEuAKuDbAO5+UTh/trvnu/sDAzz/XcBbgFOAdxCE+WeAcoL9Jx+PWfY3QC1QASwhCDHc/fbw/pfDdb3DzNKAR4CXgErgTcAnzezS8LX+nSAUpwOXAsfqR30rcFFYYxHwHoJwiNdVwLnALOA+4O/MzADMrCR8/ftjn+Du24HnCLZPj/cCD7p7B8GXzH8RBP+pQDVBiB2Tu/878CWCL958d78zzvfwXuAfCbZ9FvCpsP7JBJ/Ltwk+sznAsv4+k35e87PA68LnzAbOAWK7zMYTbO9Kgi+e28LtNVQfCF+nGigF/ongy7I/5wKrgDLgy8CdPZ8VcC/wt/A1Pg+8/zhqGfUU6NHwsJk1xNw+Ek7vIPgXfqK7H3L3Pw/xdb/t7rvcfRvwDPC8uy9190PALwlapQC4+w/dvcnd2wj+oGabWdEAr3s2UO7ut7p7u7uvB+4AFoTz3wP8p7vvc/ctwLeOUWMHUADMBMzdV7j7jiG8x/8K13MwfI8OXBjOuwZ4Lgzwvu4FFgKEobIgnIa7r3X3J9y9zd33AF8DLh5CTUP1I3dfHb6HnxGEMARB/3t3v8/dO9y93t2Xxfma7wNudffd4Xv4AkeGZEc4v8PdHyNo6R/PPpgOghCucfcud3/R3RsHWHaTu9/h7l3Aj4EJwDgzm0TwO/W58Pfpz8Ci46hl1FOgR8NV7l4cc7sjnP6vBK3Fv4VHMHxwiK+7K+b+wX4e5wOYWbqZ/beZrTOzRoJ/sSFoSfVnMjAx9kuIoOU/Lpw/kaB7oMemgQp09z8A3yH4F3u3md1uZoVxvbtA73rc3Qla4wvDSe8l/E+jH78AzjOzCQT/IXQTfCFgZuPM7H4z2xZuj58y8LYYDjtj7rcSfi4Erd51x/maEzlyu28Kp/Wo9yP7+WPXOxR3A48D95vZdjP7spllDrBs7/t099bwbn5Y176YaXDk74/ESYEeYe6+090/4u4TgY8C3zWzmhFY1XuB+QT9y0UE/bwQfJlA0OqNtQXY0OdLqMDd3xbO30EQRj0mHWvl7v4tdz+LoNvkFII+dwj6/nNjFh3f39P7PL4PuCbsrjiXILj7W+d+gu6svyN4//eHXwgQdJs4cIa7FwJ/z+FtMVTxvIeBbCHoturPYMOkbif44u0xKZx2PAZ8D2EL/wvuPgs4H7gC+Ichvv4OYKyZxa6jeqCFZWAK9Agzs3ebWVX4cD/BH3F3+HgXMG2YVlUAtBH0XecSBFqsvuv6G9BkZp8Oj7tON7PTzezscP7PgFvMrCSs/4aBVmxmZ5vZuWGrrgU4xOH3uAy42sxywy+yDw32Rtx9KcGOzB8Aj7t7wzEWv5cgfK4J7/coIOiCOGBmlRz+gjkey4CLzGxS2IV1yxCeew/wZjN7j5llmFmpmfV0xwz2+d8H/JuZlZtZGfA5gv80jscyYIGZZZpZHcH2AsDM3mBmZ5hZOsFO7Q4Of35xcfdNwGLg82aWZWbnEezzkSFSoEfDI2bWHHP7ZTj9bOB5M2sm6FP8RNhfDUE/94/DLo/3nOD6f0LwL/k2gqNS/tpn/p3ArHBdD4d9oFcQ9PVu4HCA9vS5fyF8vQ0EreC7j7HuQoL+9/3hc+qBr4Tzvk5wJMcugj7XgbpP+rqX4L+NewdZbhHBjuCd7v5SzPQvAPMIjvj5NfBQnOs9irs/ATxAcLTRi8CjQ3juZuBtBEej7CMI1tnh7CM+k36e/h8EIfkywZFDS8Jpx+P/EPynsJ9g28Ru1/HAgwRhvoLgqKxjfd4DeR9wHsHn/x8E26ztOOsdtezwf5kiItFgZg8AK8MjhyROaqGLSMKFXW/TzSzNzC4j2KfT338ecgw6U1REomA8QddWKcFJT/8r3B8iQ6AuFxGRFKEuFxGRFKFAlxERpeF8R2IYVzP7k5l9eDhfcziY2V1mdrxHs0iSUx+6jAgfweF8h8o1jKuMEmqhiyQhM1NjTI6iQJfjFl655xYzW25m+83sR2Y2Jpx3xFWKzGyimf3CzPaY2QYz+3jMvHQz+0w4lkyTmb1oZtXhvJlm9oSZ7TOzVT0nUZnZ1PCkmrTw8R1mtjvmNe82s0+G93u7R8ysxsyesuCKQnvD45051rqOYbKZPRvW/LvwjMye1/q5me0M1/O0mZ0Wsx1iTyJrteBiJYSH7f3BzOrD2u4xs+I+2/vTZvYy0BKePTrXzJaENTxAMJxtz/JlZvZouJ32mdkzPdtLUpM+XDlR7yMYInc6wTgsR13VyAYfbvefCQbUehvBmaMfBFrNLA94guDMxAqCERG/a2az3H0DwdmJPSNGXgQ0m9mp4eOLCc5a7KvfIYmPta5jvPd+h70NDTQc8fZwyNt8d88nGPWyZ3jfeIbtXQi8HSgm+Pt9mODMzLEEY9rHDgl8E8EhgOUEA6d9hsHHgJEkpkCXE/Udd9/i7vuA/+TwSIexBhtu98PAv7n7Kg+85O71BMMLbHT3H7l7Z3hc8i+Ad4fPewq42Mx6Bot6MHw8leCLIfZ0/h4DDUk82Lr6M9Cwt3ENR2xmnyYYNviD4XPiGbb3W+H2Pkgw3nkm8I1wkKwHCS4wEfteJwCTw/nPxAxAJilIgS4nqu8wuRP7WWaw4XYHGiZ2MnBun+e9j8Oj/T1FcPWdiwiudvMnggC8GHjG3fsbJGqgIYkHW1d/+h321uIYjtjMLgc+QTB08sFwWjzD9sZu74nAtj4hHTtk7leAtcDvzGy9md18jPciKUA7VuRE9R0mt78hWnuG260d4DV6hol9tZ/pT7n7WwZ43lMEobU1vP9n4PsEIzb2192Cu+8EPgJgZq8Hfm9mT8exrqGIHY54I8GgZfsJh+A1sxkEg41dHV4ApEfssL37zOwqgrHij3gLMfd3AJVmZjGhPonwy9Hdmwi6XW4ys9OBP5jZC+7+5DC8R4kgtdDlRH3MzKrMbCzBZc/6uxzeYMPt/gD4opnVWuBMMyslGJnwFDN7vwVDt2ZaMObHqQDuvobgQh1/TxDGjQQjM76LAQLdBh6S+JjrGqIBhyO24OIdvwI+60dfgWqow/Y+B3QCHw/rvZrgUnM967oi3AlsBCNHdjHEoW0luSjQ5UTdS7CTcT1By/Cok1riGG73awR90L8j2NF5J5ATtjDfStDXvp2gi+N/CC6A3OMpgqvvbIl5bAQ7IvvT75DEca4rXscajngewfH5X4892iWcN6Rhe929Hbia4ALM+wgu1hH7nFqCCy43E4T/d939j8fxfiRJaCwXOW5mtpHgqu6/T3QtIqIWuohIyhg00M3sh2a228z67rDqmW9m9i0zW2tmL5vZvOEvU0REBhNPC/0u4LJjzL+coK+uFrgO+N6JlyXJwN2nqLtFJDoGDXR3f5pgh8tA5gM/CU8I+StQbGYThqtAERGJz3Ach17JkSc7bA2n7ei7oJldR9CKJy8v76yZM2cOeWW7Gg+xuym4dmx6mpGblU5uVkb4M500s+N4CyIiyeHFF1/c6+7l/c07qScWufvtwO0AdXV1vnjx4iG/Rle3s3Z3M0s272fJpv0s2byfdXtaOAS0G5wyroB5k0uYW13MvMklTC3NIy1NIS8iqcHMNg00bzgCfRtHni1YFU4bEelpxozxBcwYX8DCcyYBcKC1g6Vb9rNkcwNLN+/nkZe2c+/zmwEoyslkTnUxcycVM29SCbOriynKyRyp8kREEmY4An0RcL2Z3Q+cCxxw96O6W0ZSUW4ml8yo4JIZFQB0dzvr9jSzdHNDEPSbGvjmk2voOeS+piKfeZOKmTuphLmTiqmtKCBdrXgRSXKDnlhkZvcRDIBURnBa9b8TjPCGu38/PK34OwRHwrQC/+jug/alHG+Xy/FqOtTBS1sOsGTzfpZtCVry+1s7AMjPzmB2dRFzq4OAn1NdTGn+8ZwgKCIysszsRXev63deos4UPdmB3pe7s7G+laWb97N0cwNLNu9n5c4murqD7TG5NJe51Ydb8adOKCQzXedhiUhiHSvQR+1oi2bG1LI8ppblcfW8YKymg+1dvLy1gSWbG1i2ZT/Prqvn4WXB4IHZGWmcUVnE3JiumglFOYl8CyIiRxi1LfR4uDvbGg6GXTRBN82r2xtp7wwGrBtXmM3c6hLmTCpmbnUxZ1QVkZs1ar8jReQkUAv9OJkZVSW5VJXkcsWZwXUb2ju7WbGjkaU9ffFbGvjta8F1DtLTjBnjCnpb8XOqi5lWpsMmReTkUKAPUVZGGrOri5ld3XvtXuqb21i2paH3tmjZdu4JD5ssHJPB7Ori3v74OdXFlORlJap8EUlh6nIZAb2HTYZdNcu2NLBqZyPh/lamlOaGx8YHAX/qhEKyMrTDVUQGp6NcIqClrZOXtx4IW/HBkTU9QxhkZaRx2sRC5lQHJz9detp4BbyI9EuBHkHuzo4Dh3qPiV+2pYFXth3gUEc3X5x/Gu8/b0qiSxSRCNJO0QgyMyYW5zCxOIe3nREMTtnR1c2b/u9TPLV6jwJdRIZM/9dHSGZ6GhfWlvHcuno6unQtXxEZGgV6xFxYW05LexdLNzckuhQRSTIK9Ig5b3op6WnGM2v2JLoUEUkyCvSI6Rnu9+k1exNdiogkGQV6BF1YW8bLWxtoaG1PdCkikkQU6BF0YW057vDs2vpElyIiSUSBHkGzq4ooGJOhfnQRGRIFegRlpKdxwfQynlmzl0Sd+CUiyUeBHlEXnlLGtoaDrN/bkuhSRCRJKNAj6sKacgCeXaujXUQkPgr0iKoem0PhmAxW72pKdCkikiQU6BFlZkwrz2f9HnW5iEh8FOgRNq0sT4EuInFToEfYtPI8djYeoqWtM9GliEgSUKBH2LTyfAA26EgXEYmDAj3CppXnAbBuT3OCKxGRZKBAj7AppXmYoX50EYmLAj3CxmSmU1mcoy4XEYmLAj3ippXns36vulxEZHAK9IibVpbHhj0tGtNFRAalQI+46eV5tLR3sauxLdGliEjEKdAjbmpZcOjieh3pIiKDUKBHXO+hi9oxKiKDUKBH3PjCMeRkpquFLiKDUqBHXFqaMVVjuohIHBToSWBaeZ4OXRSRQSnQk8C08ny27j/IoY6uRJciIhEWV6Cb2WVmtsrM1prZzf3Mn2xmT5rZy2b2JzOrGv5SR6/p5Xm4w6b61kSXIiIRNmigm1k6cBtwOTALWGhms/os9lXgJ+5+JnAr8F/DXehoNk2HLopIHOJpoZ8DrHX39e7eDtwPzO+zzCzgD+H9P/YzX07A1PDQRV0wWkSOJZ5ArwS2xDzeGk6L9RJwdXj/nUCBmZX2fSEzu87MFpvZ4j179hxPvaNSfnYG4wqzWbdbLXQRGdhw7RT9FHCxmS0FLga2AUftwXP32929zt3rysvLh2nVo0NNRT5r1eUiIscQT6BvA6pjHleF03q5+3Z3v9rd5wKfDac1DFuVQk15Put2N2uQLhEZUDyB/gJQa2ZTzSwLWAAsil3AzMrMrOe1bgF+OLxlSk1FPi3tXew4cCjRpYhIRA0a6O7eCVwPPA6sAH7m7q+Z2a1mdmW42CXAKjNbDYwD/nOE6h21plcER7rocnQiMpCMeBZy98eAx/pM+1zM/QeBB4e3NIlVEwb62t3NXFir/Q8icjSdKZokyvOzKRyTwVod6SIiA1CgJwkzC450UaCLyAAU6EmkpiJffegiMiAFehKpqchnb3M7Da3tiS5FRCJIgZ5EYneMioj0pUBPIjXlBYACXUT6p0BPIpUlOWRnpCnQRaRfCvQkkh5ejk5juohIfxToSUaHLorIQBToSaamIp9tDQc52K7L0YnIkRToSaamIh93dNFoETmKAj3J6NBFERmIAj3JTC3LI83Q1YtE5CgK9CSTnZHOlLI8lu9oTHQpIhIxCvQkNKeqmGVbDujqRSJyBAV6EpozqZi9zW1s19WLRCSGAj0Jza4qBuClLbpsq4gcpkBPQqdOKCQrPY1lCnQRiaFAT0JZGWnMmlioQBeRIyjQk9Sc6mJe2XqAzq7uRJciIhGhQE9Sc6qLOdjRxRodjy4iIQV6kppdHewYVbeLiPRQoCepKaW5FOVk6kgXEemlQE9SZsbs6mK10EWklwI9ic2pLmb1riZa2joTXYqIRIACPYnNqS6i2+HVbQcSXYqIRIACPYn1nDGqbhcRAchIdAFy/Erzs5lcmst//3YlD7ywhVMnFHLqhILwZyETisZgZokuU0ROEgV6krvtvfN4csVuVuxo5JVtB/j1Kzt65xXlZB4R8LMmFFJTkc+YzPQEViwiI0WBnuROryzi9Mqi3sdNhzpYtbOJFTsaWb4j+Hnf3zZzqCM4ozQ9zZhentcb8j2t+vL8bLXmRZKcAj3FFIzJpG7KWOqmjO2d1tXtbKxvYWUY8Ct2NPLChn38atn23mVK87I4dUIhM8cfbtHXVOSTlaHdLCLJQoE+CgSt8nyml+fz9jMn9E5vaG1nxY4mVu5sDIO+iZ/8dRPtnUFrPjM9eN6sCYXMDLtuZo4vpLwgO1FvRUSOQYE+ihXnZnHe9FLOm17aO62zq5sNe1tYvqORlTubWLmjkb+sq+ehpdt6lynLz+7tm+9p0U8vV2teJNEU6HKEjPQ0ascVUDuugPkx0/e1tLNyR+PhoN/ZyF1/2XhUa74n5Geqb17kpIsr0M3sMuCbQDrwA3f/7z7zJwE/BorDZW5298eGuVZJoLF5WZxfU8b5NWW903pa8yvCnbArdzTy3Lp6fhnTmi/Ny2LmhAJmjj/cmteRNiIjwwa70LCZpQOrgbcAW4EXgIXuvjxmmduBpe7+PTObBTzm7lOO9bp1dXW+ePHiEyxfomh/S3tvK35F2KJftbOJts7DR9pMK8tjZtianzWhkHmTSyjKyUxw5SLRZ2Yvuntdf/PiaaGfA6x19/Xhi90PzAeWxyzjQGF4vwjYjoxaJXlH9813dTsb9rawKibol2zazyMvBb8q1WNzeOLGi9VyFzkB8QR6JbAl5vFW4Nw+y3we+J2Z3QDkAW/u74XM7DrgOoBJkyYNtVZJYulpRk1FPjUVRx5p03iogydX7OLGB17ip3/dxIcvnJbAKkWS23AdlrAQuMvdq4C3AXeb2VGv7e63u3udu9eVl5cP06olmRWOyeSdc6u4sLaM7/5pHc0aOVLkuMUT6NuA6pjHVeG0WB8Cfgbg7s8BY4AyROL0qbfOYF9LOz/884ZElyKStOIJ9BeAWjObamZZwAJgUZ9lNgNvAjCzUwkCfc9wFiqpbXZ1MZeeNo47nl7P/pb2RJcjkpQGDXR37wSuBx4HVgA/c/fXzOxWM7syXOwm4CNm9hJwH3CtD3b4jEgfN711Bs3tnXz/6XWJLkUkKcV1HHp4TPljfaZ9Lub+cuCC4S1NRptTxhVw1ZxKfvKXTXzsDTUUjtFhjCJDoXO1JVKuPX8KBzu6jhg4TETio0CXSDmzqohTJxRy/982J7oUkaSjQJdIMTMWnlPNa9sbeWWrrpUqMhQKdImc+XMqGZOZxn0vqJUuMhQKdImcopxM3n7GRBYt206LTjQSiZsCXSJp4TnVNLd18uuXdwy+sIgACnSJqLMml1BTkc+92jkqEjcFukSSmbHg7GqWbWlgxY7GRJcjkhQU6BJZ15xVRVZGGvc+r1a6SDwU6BJZxblZXHHGBH65dJt2jorEQYEukfbecyfR3NbZeyEMERmYAl0i7azJJcwYV8A96nYRGZQCXSLNzHjvuZN4ZdsBXt7akOhyRCJNgS6R9855leRkpmvnqMggFOgSeYVjMnnH7An8atl2djcdSnQ5IpGlQJek8NGLp9Ptzi2/eAVdO0Wkfwp0SQrTy/P59GUzeXLlbh54YUuiyxGJJAW6JI1rz5/CBTWlfPHR5Wyub010OSKRo0CXpJGWZnzlmtmkpRk33L+UXy7dyjNr9rB8eyO7Gw/R2dWd6BJFEiqua4qKRMXE4hy+9M4zuPGBZdz4wEtHzS/JzaQsP5vS/CzK8rPDWxal+dmU5mVRVpBNeTg9Jys9Ae9AZOQo0CXpvGP2RN4ws4I9TW3sbW5jb1Mbe1vag5/NbdQ3t7O3uY1Xtx2gvrmdpgGGDcjNSu8N/9K8nuDPCqcFj3u+FIpzMklLs5P8TkWGRoEuSSk/O4P87AymluUNuuyhji7qW9qpD8N+T0zo93wBbN3fyrItDexvbaer++ijaNLTjLF5WUErP6bVf/i/gfBLoSD4T2BMplr/cvIp0CXljclMp7I4h8rinEGX7e52Gg529Ib93ubgi2Bvcxv7WtrZG34RbN7cyt7mNlrbu/p9nfzsjCO6ekrzsynveRy2/CcW5TCpNHe4366MYgp0kRhpYUt8bF4Wp4wrGHT51vbOmNZ+O/ta2npDv2f6xvoWlmzez76Wdvo2/r/67tlcc1bVCL0bGW0U6CInIDcrg9yxGVSPHbyl3dXt7G9tpz5s9f/P46v4yuMrueLMCeqikWGhwxZFTpL0NKMsP5sZ4ws4v6aMz1w+k12Nbfz4LxsTXZqkCAW6SIKcO62Ui08p53tPraPxUEeiy5EUoEAXSaB/uXQGDa0d3PH0+kSXIilAgS6SQKdXFnHFmRO4888b2NPUluhyJMkp0EUS7Ka3zqCts5vvP7Uu0aVIklOgiyTY1LI83jm3knue36Tx3uWEKNBFIuD6N9TQ0eXc/pT60uX4KdBFImBKWR7z50zkp89vUl+6HDcFukhE3PDGWto7u7njGbXS5fgo0EUiYmpZHlfNqeQnz21kb7Na6TJ0cQW6mV1mZqvMbK2Z3dzP/K+b2bLwttrMGoa/VJHUd/0ba2jv7OZ2HZcux2HQQDezdOA24HJgFrDQzGbFLuPuN7r7HHefA3wbeGgkihVJddPK83tb6epLl6GKp4V+DrDW3de7eztwPzD/GMsvBO4bjuJERqMb3lRLR5fruHQZsngCvRKIvcz61nDaUcxsMjAV+MMA868zs8VmtnjPnj1DrVVkVJhalsfVcyv56V83satRx6VL/IZ7p+gC4EF373fUf3e/3d3r3L2uvLx8mFctkjpueGMtXd3O9/6kVrrEL55A3wZUxzyuCqf1ZwHqbhE5YZNKc7nmrCrufX4zOw4cTHQ5kiTiCfQXgFozm2pmWQShvajvQmY2EygBnhveEkVGp4+9oQbH+eT9y2hobU90OZIEBg10d+8ErgceB1YAP3P318zsVjO7MmbRBcD97n70FXZFZMiqx+bylWtms3RzA1fd9ixrdzcluiSJOEtU/tbV1fnixYsTsm6RZPLipn189O4Xaevo5p8umU5NRT5Ty/KYNDZXl64bhczsRXev63eeAl0k+rY1HOR/37OEl7YcPmfPDCYW5TCtPI8ppXlMKctjalkuU0rzqB6bS2a6TgRPRccKdF0kWiQJVBbn8KuPXcCBgx1s3NvCxvoWNuw9fHt42TaaDnX2Lp+eZlSX5DClLAj7qWVB4E8ry2NicQ7paZbAdyMjRYEukkSKcjKZXV3M7OriI6a7O/ta2sOgb2VjTNj/bcM+WtsPH0mclZ5G9dicIOR7W/bBzwmFY0hT2CctBbpICjAzSvOzKc3P5qzJY4+Y5+7saWpjQ2/LvpUNe5vZuLeVZ9bspa2zu3fZrIw0Jo/NPRzypXlMKctlalke4woU9lGnQBdJcWZGReEYKgrHcO600iPmdXc7OxsPBS36+pawZR+08J9avYf2mLAfk5kWBHxMf/3ksDunoiAbM4V9oinQRUaxtDRjYnEOE4tzOL+m7Ih5Xd3OjgMH2bi3lQ31LWwKu3BW727iyZW76Og6fEBFTmY6k0tzj9o5O0Vhf1Ip0EWkX+lpRlVJLlUluby+9uiw395wsLcbZ+PeVjbWt7B619Fhn5uVzrvmVXHr/NMU7CNMgS4iQ5aeZlSPzaV6bC4XceS4TJ1d3WxvOBQEfX2wU/buv25i1sRCFp4zKUEVjw4KdBEZVhnpaUwqzWVSaRD2f3/uZBpaO/jCI69x9pQSaioKEl1iytKZByIyotLSjK+9Zza5WRnccN8yDnX0OxirDAMFuoiMuIrCMXzlmjNZsaORL/92VaLLSVkKdBE5Kd506jje/7rJ/OgvG3h124FEl5OSFOgictJ86tIZlORm8flFr6GBWYefAl1ETpqinEz+9dIZLN60n2w8O4kAAArXSURBVEUvbU90OSlHgS4iJ9W766o5o7KILz22gpa2zsGfIHFToIvISZWeZnz+ytPY1djGd/64NtHlpBQFuoicdGdNLuHqeZXc+cwGNuxtSXQ5KUOBLiIJcfPlM8nKSNMO0mGkQBeRhKgoGMMn31zLU6v38MTyXYkuJyUo0EUkYT5w/hROGZfPrY8u1xmkw0CBLiIJk5mexq3zT2fr/oN890/rEl1O0lOgi0hCvW5aKfPnTOS7f1zLoy/r2PQTodEWRSThvnjV6exoOMQN9y2l+VAnCzTM7nFRC11EEq5wTCY//uA5XFRbzs0PvcI3f7+G9Xua6erW0S9DYYk6XKiurs4XL16ckHWLSDS1d3Zz4wPL+PUrO4DgotXTy/M5ZVw+p4wroLYinxnjC6guyR21F6w2sxfdva7feQp0EYkSd+fVbY2s2tXE6vC2Zlcz2xoO9i4zJjONmop8TqkooHZcQW/gVxbnpHzQHyvQ1YcuIpFiZpxRVcQZVUVHTG861MGa3c2s2dXE6l3NrN7VxF/W1fPQ0m29y+RmpVNTkU9tRQG144KWfW3F6Ah6UKCLSJIoGJPJvEklzJtUcsT0Awc7WLv7cMiv2dXMM2v28IslW3uXyc1Kp7Yin5qKoDVfm6JBr0AXkaRWlJPJWZPHctbksUdMP9DawZrYoN/d1G/Qp1KLXoEuIimpKDeTuiljqZtyZNA3tLaHXTcDB31OZhj0YcD3BH1VSbSDXoEuIqNKcW4WZ08Zy9lTjm7Rr97dxNrdQdCv3d3Ms2v38tCSw330PTtjaysKwp/51I4rYNLYXNIjEPQKdBERghZ9v0Ef00ffE/bPravnlzE7Y3sOr6ztDfl8ZowvZEppLmYnL+gV6CIixzBQH33ToQ7W7m5mze4g6NfsauLFPpfWO21iIQvPmcT8ORMpGJM54rXqOHQRkWHU0tbJuj3NLNm0n/tf2MLKnU3kZqXzpXeewVVzK0/49Y91HHpcp/6b2WVmtsrM1prZzQMs8x4zW25mr5nZvSdSsIhIssrLzuDMqmKuvWAqv/nEhTz8sQs4vbKITz6wjNufXjeiF/MYNNDNLB24DbgcmAUsNLNZfZapBW4BLnD304BPjkCtIiJJxcyYU13M3R86h7efOYEvPbaSLz66gu4RGqMmnj70c4C17r4+LPB+YD6wPGaZjwC3uft+AHffPdyFiogkq+yMdL69YC7jCsbww2c3UFGYzT9dPH3Y1xNPoFcCW2IebwXO7bPMKQBm9iyQDnze3X/b94XM7DrgOoBJkzQ8poiMHmlpxv+54lRmTSzk8tPHj8w6hul1MoBa4BJgIXCHmRX3Xcjdb3f3OnevKy8vH6ZVi4gkBzPjmrOqyMsemQMM4wn0bUB1zOOqcFqsrcAid+9w9w3AaoKAFxGRkySeQH8BqDWzqWaWBSwAFvVZ5mGC1jlmVkbQBbN+GOsUEZFBDBro7t4JXA88DqwAfubur5nZrWZ2ZbjY40C9mS0H/gj8i7vXj1TRIiJyNJ1YJCKSRE74xCIREYk+BbqISIpQoIuIpAgFuohIilCgi4ikCAW6iEiKUKCLiKQIBbqISIpQoIuIpAgFuohIilCgi4ikCAW6iEiKUKCLiKQIBbqISIpQoIuIpAgFuohIilCgi4ikCAW6iEiKUKCLiKQIBbqISIpQoIuIpAgFuohIilCgi4ikCAW6iEiKUKCLiKQIBbqISIpQoIuIpAgFuohIilCgi4ikCAW6iEiKUKCLiKQIBbqISIpQoIuIpAgFuohIilCgi4ikiLgC3cwuM7NVZrbWzG7uZ/61ZrbHzJaFtw8Pf6kiInIsGYMtYGbpwG3AW4CtwAtmtsjdl/dZ9AF3v34EahQRkTjE00I/B1jr7uvdvR24H5g/smWJiMhQDdpCByqBLTGPtwLn9rPcu8zsImA1cKO7b+m7gJldB1wXPmw2s1VDrLdHGbD3OJ97skS9xqjXB6pxOES9Poh+jVGrb/JAM+IJ9Hg8Atzn7m1m9lHgx8Ab+y7k7rcDt5/oysxssbvXnejrjKSo1xj1+kA1Doeo1wfRrzHq9cWKp8tlG1Ad87gqnNbL3evdvS18+APgrOEpT0RE4hVPoL8A1JrZVDPLAhYAi2IXMLMJMQ+vBFYMX4kiIhKPQbtc3L3TzK4HHgfSgR+6+2tmdiuw2N0XAR83syuBTmAfcO0I1gzD0G1zEkS9xqjXB6pxOES9Poh+jVGvr5e5e6JrEBGRYaAzRUVEUoQCXUQkRSRdoA82DEEC6qk2sz+a2XIze83MPhFOH2tmT5jZmvBnSQRqTTezpWb2aPh4qpk9H27LB8Kd3omqrdjMHjSzlWa2wszOi9o2NLMbw8/4VTO7z8zGJHobmtkPzWy3mb0aM63f7WaBb4W1vmxm8xJU31fCz/llM/ulmRXHzLslrG+VmV060vUNVGPMvJvMzM2sLHx80rfhUCRVoMcMQ3A5MAtYaGazElsVncBN7j4LeB3wsbCmm4En3b0WeDJ8nGif4MgjkP4H+Lq71wD7gQ8lpKrAN4HfuvtMYDZBnZHZhmZWCXwcqHP30wkOEFhA4rfhXcBlfaYNtN0uB2rD23XA9xJU3xPA6e5+JsGJiLcAhH83C4DTwud8N/ybT0SNmFk18FZgc8zkRGzD+Ll70tyA84DHYx7fAtyS6Lr61PgrgnFvVgETwmkTgFUJrquK4I/7jcCjgBGc/ZbR37Y9ybUVARsId9LHTI/MNuTwGdNjCY4OexS4NArbEJgCvDrYdgP+H7Cwv+VOZn195r0TuCe8f8TfM8GRdeclYhuG0x4kaFxsBMoSuQ3jvSVVC53+hyGoTFAtRzGzKcBc4HlgnLvvCGftBMYlqKwe3wD+FegOH5cCDe7eGT5O5LacCuwBfhR2Cf3AzPKI0DZ0923AVwlaazuAA8CLRGcbxhpou0Xx7+eDwG/C+5Gpz8zmA9vc/aU+syJTY3+SLdAjy8zygV8An3T3xth5HnyVJ+z4UDO7Atjt7i8mqoZBZADzgO+5+1yghT7dKxHYhiUEg9JNBSYCefTzb3rUJHq7HYuZfZagy/KeRNcSy8xygc8An0t0LUOVbIE+6DAEiWBmmQRhfo+7PxRO3tVzBm34c3ei6gMuAK40s40Eo2W+kaDPutjMek4uS+S23Apsdffnw8cPEgR8lLbhm4EN7r7H3TuAhwi2a1S2YayBtltk/n7M7FrgCuB94ZcORKe+6QRf3C+FfzNVwBIzG090auxXsgX6oMMQnGxmZsCdwAp3/1rMrEXAB8L7HyDoW08Id7/F3avcfQrBNvuDu78P+CNwTbhYwmp0953AFjObEU56E7CcCG1Dgq6W15lZbviZ99QYiW3Yx0DbbRHwD+GRGq8DDsR0zZw0ZnYZQfffle7eGjNrEbDAzLLNbCrBjse/nez63P0Vd69w9ynh38xWYF74exqJbTigRHfiH8fOi7cR7BlfB3w2AvW8nuBf2peBZeHtbQR91E8Ca4DfA2MTXWtY7yXAo+H9aQR/MGuBnwPZCaxrDrA43I4PAyVR24bAF4CVwKvA3UB2orchcB9Bn34HQfB8aKDtRrAj/Lbwb+cVgiN2ElHfWoJ+6J6/l+/HLP/ZsL5VwOWJ2oZ95m/k8E7Rk74Nh3LTqf8iIiki2bpcRERkAAp0EZEUoUAXEUkRCnQRkRShQBcRSREKdBGRFKFAFxFJEf8fYYAKUNCRRKMAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "t = np.linspace(.001, 150, 100)\n", "plt.plot(t, survival_function(estimates_, t))\n", "plt.ylim(0.5, 1)\n", "plt.title(\"\"\"Estimated survival function using \\npiecewise hazards\"\"\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On day 30, we charge users \\\\$10, and on every 30 days after that, we charge \\\\$20. What's the LTV, and CIs, at the end of day 120?" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def LTV_120(params):\n", " # think about how complicated the gradient of this function is. Now imagine an even more\n", " # complicated function.\n", " ltv = 0\n", " ltv += 10 * survival_function(params, 30)\n", " for t in [60, 90, 120]:\n", " ltv += 20 * survival_function(params, t)\n", " return ltv" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "LTV estimate: [51.82507427]\n" ] } ], "source": [ "ltv_ = LTV_120(estimates_)\n", "print(\"LTV estimate: \", ltv_)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Variance LTV estimate: 1.605494065103963\n" ] } ], "source": [ "from autograd import grad\n", "var_ltv_= grad(LTV_120)(estimates_) @ variance_matrix_ @ grad(LTV_120)(estimates_)\n", "print(\"Variance LTV estimate:\", var_ltv_)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Estimated LTV at day 120: ($49.34, $54.31)\n" ] } ], "source": [ "std_ltv = np.sqrt(var_ltv_)\n", "print(\"Estimated LTV at day 120: ($%.2f, $%.2f)\" % (ltv_ - 1.96 * std_ltv, ltv_ + 1.96 * std_ltv))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From here, we can compute p-values, scenario analysis, sensitvity analysis, etc. \n", "\n", "Let's continue this analytical train to Part 8. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Bonus, if time permits and Cameron is talking too fast. \n", "\n", "In the above model, we are not suggesting to the model much apriori information about this \"predictable\" process. For example, suppose we want a model that gives \"inbetween\" period rates to be close to each other, and likewise for the \"jump\" rates. " ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " fun: 6839.426915809169\n", " hess_inv: <10x10 LbfgsInvHessProduct with dtype=float64>\n", " jac: array([-46.78891729, 23.72599913, -12.19979075, 40.04503677,\n", " 63.39941388, -7.47615848, 8.16437802, -27.76700067,\n", " -11.54179562, -28.46781927])\n", " message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n", " nfev: 45\n", " nit: 32\n", " status: 0\n", " success: True\n", " x: array([0.000976 , 0.01823817, 0.000976 , 0.01823816, 0.00097599,\n", " 0.01823816, 0.000976 , 0.01823816, 0.00097601, 0.01823816])\n" ] } ], "source": [ "\n", "def negative_log_likelihood(params, t, e):\n", " return -log_likelihood(params, t, e) + 1e12 * (np.var(params[::2]) + np.var(params[1::2]))\n", "\n", "from autograd import value_and_grad\n", "\n", "results = minimize(\n", " value_and_grad(negative_log_likelihood), \n", " x0 = np.ones(len(breakpoints)),\n", " method=None, \n", " args=(T, E),\n", " jac=True,\n", " bounds=[(0.0001, None)] * len(breakpoints)\n", ")\n", "\n", "print(results)\n", "estimates_ = results.x" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEmCAYAAACJXlw1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2dfbwcVZnnv7/bl6DIqyG8BkmUqBNfQI2oO+o4ogi+BR3UMKzCDA6zKqvu6o6gI4Oou+o64jqisygog/Imvkx08YU3UWZGICCCoJHIi0kIEAnvEMLtfvaPczq30re7q6rT996uU8/387mf2111quqcqq7znOec33mOzAzHcRzHyTI22xlwHMdxRg83Do7jOM4U3Dg4juM4U3Dj4DiO40zBjYPjOI4zBTcOjuM4zhTcODh9kfQySStnOx/dkPQKSWv67DdJ+81knoog6TZJr+qx74mSvi/pfknfmuF83SjpFTN5zX5I+qGko2Y7H3VlfLYz4EwPkm4Ddgeamc1fN7Pjco4zYJGZrQIws58Dz5imPH4dWGNmfz8d568ohxOe21wzm5iui3S792b2rOm63iCY2aGznYc648Yhbd5gZhfPdibqiqTxASr4fYHfTadhcJwieLdSDZG0n6TLY9fFHyWdF7f/LCb5laSHJL2ts+smdon8D0nXS3pY0umSdo9dAA9KuljSLpn035J0Z7zWzyQ9K24/FjgS+Lt4re/H7XtJ+rak9ZJulfTezLmeKOnrku6VdBPwwgLFfZWkmyXdJ+lUSYrnepqkSyXdE+/BNyXtHPe9Leap/feYpJ/Gfa+T9EtJD0haLemkTP4WxK6sYyT9Abg0bn+7pNvjtT7S57l8DDgRaF//GEknSfpGl2uMx+8/lfRxSf8W7/9PJO2aSf9SSf8ey79a0tF97v3m7i5J20r6vKQ74t/nJW0b971C0hpJH5B0t6R1kv6qT7m26EbLlknSEyR9I96b+yRdLWn3TNneGT8fLekKSZ+Nz/9WSYdmzrkw/r7av8FTs/fNKY8bh3ryceAnwC7AfOCfAMzs5XH//ma2vZmd1+P4vwBeDTwdeAPwQ+DDwDzCb+q9mbQ/BBYBuwHXAt+M1zotfv5MvNYbJI0B3wd+BewNHAS8X9Jr4rn+AXha/HsNUKQ/+vUEI/Jc4K3xOAAB/wvYC/gTYB/gpJi382Keto/7bwHOicc9DLwD2Bl4HfAuSYd1XPPP4jlfI2kx8GXg7fFccwn3fApm9g/A/wTa1z+9QPkA/hL4K8I9ngN8EEDSvoT7/0+EZ3MAcF23e9/lnB8BXhyP2R84EMh2/+0B7ER4TscAp2YbBSU4Kp5nH8K9+S/Aoz3SvghYCewKfAY4vW3sgbOBq+I5TiLcb2crcOOQNt+LrbH239/E7Y8Tui/2MrONZnZFyfP+k5ndZWZrgZ8DV5rZL81sI/Bd4HnthGZ2hpk9aGaPEV7a/SXt1OO8LwTmmdnJZrbJzG4BvgIsi/vfCnzSzDaY2WrgCwXy+ikzu8/M/gBcRqjsMLNVZnaRmT1mZuuBzxEq9c1EY3U28FMz+7/xuJ+a2Q1m1jKz6wlGY4vjgJPM7GEze5QwhvADM/tZvAcfBVoF8l2Gr5nZ7+L1zm+XkWA0Ljazc8zscTO7x8yuK3jOI4GTzezueH8+xpYV7uNx/+NmdiHwEIONTT1OqND3M7OmmV1jZg/0SHu7mX3FzJrAmcCewO6SnkL47ZwYfzdXAMsHyIuTwY1D2hxmZjtn/r4St/8doeV8lYJC5a9LnveuzOdHu3zfHkBSQ9KnJP1e0gPAbTHNrnRnX2CvrEEjeCS7x/17Aasz6W8vkNc7M58fyeRtd0nnSlob8/aNLvn6JLADGU9I0oskXabQ7XU/oaXbeVw2j1vk2cweBu4pkO8ydC0joTX++wHPuRdb3t/b47Y293SMi2SvW4azgB8D58buq89I2qZH2s3lNLNH4se2d7chsw22fAbOALhxqCFmdqeZ/Y2Z7QX8LfAlTY/k8y+BpcCrCF0HC+L2dldAZ0jg1cCtHQZtBzN7bdy/jlDhtXnKVuTtf8brP8fMdgT+cyZfSFoGHAEcbmaPZ447m9Aq3cfMdgL+OXtcl3JtkWdJ2xFaykV5GNgu832PEseuJnTBdSMvHPMdBGPd5ilx2yD0LEP0PD5mZouB/0ToBnxHyfOvA54c722bfXoldorhxqGGSHqLpHa/972EiqLd1XEX8NQhXWoH4DFCS3k7QoWcpfNaVwEPSvqQwuBzQ9KzJbUHns8HTpC0S8z/f93KvD0E3C9pb+B/tHdIeh6hn/6w2KXSedwGM9so6UCCAezHBcDr48DwHOBkyr131wEvl/SU2B13Qoljv0kYkH+rpHFJcyW1u5zynvM5wN9LmhcHuE8keFeDcB2wTNI2kpYQutoAkPTnkp4jqQE8QOhmKtXtZma3AyuAkyTNkfQSwliYsxW4cUib72tL1c134/YXAldKeojQCn5f7N+HMC5wZuzWeetWXv9fCN0Ra4GbgF907D8dWByv9b3Yl/x6Qp/5rcAfga8SvA4I/d63x30/IXRJDMrHgOcD9wP/D/hOZt9SwmD9FZl798O4793AyZIeJFSY5/e7iJndCLyH4HGsIxjjnhP3uhx/EXAecD1wDfCDEsf+AXgt8AFgA6GS3j/u3uLedzn8E4QK93rgBoKY4BNFr93BRwkezL2E+352Zt8eBAP6APAb4HIGe65HAi8hNEQ+Qbhnjw2YXweQL/bjOE5qKMizfxsVYM4AuOfgOE7lkfRChbkrY5IOIXh/3TwipyA+Q9pxnBTYg9A1OJfQbfcuM/vl7Gap2ni3kuM4jjMF71ZyHMdxpuDGwRlZNELhwjUN4aOzsYNGCYX4VYMqk5xE8DEHZ2SZznDhZfHw0U7dcM/BcWqKYmRXx+mGGwdnVonhnE+QdFMMxfw1SU+I+zrDhfcL592Q9OEYx+lBSddI2ifue6akiyRtkLSyPblPIczzfTHAHpK+IunuzDnPkvT++DkbPrpryPN+1+rDvuodbrtXuPO9OiY3PqKwSFPfUOSZ+/0hSdcDD8eZ08+TdG3Mw3nAEzLpd5X0g3ifNkj6eft+OWnjD9kZBY4khNJ+GiEM+JSV4ZQfzvu/E2IhvRbYEfhr4BFJTwIuIszK3Y0Q4fVLkhab2a2EmbntKLIvBx6S9Cfx+58RZux20jXkeb9r9Sl713DbkV7hzu9ohxSPYcW/C5zbvlX0CEWe4QhCuPGdCXXA9wizkp8MfIsQkr3NBwjS0HmEAIgfJj8uk5MAbhycUeCLZrbazDYQIqEe0SVNXjjvdwJ/b2YrLfArM7uHEI7jNjP7mplNRO37t4G3xOMuB/5MUjsY3AXx+0KCkflVl7z0Cnmed61u9Aq3XSjcuaQPAc8kGMNCociBL8T7/ShhzYZtgM/HIHgXAFd3lHVPYN+4/+fm+vda4MbBGQU6w3Dv1SVNXjjvXuGp9wVe1HHckUxGBr0ceAXBa/gZ8FNCZfpnwM/NrFsQuF4hz/Ou1Y1eIcVzw50rrIT2PkKAwEfjtiKhyDtDiq/tqPCzobr/N7AK+ImkWyQd36csTkL4gJQzCnSG4e4WGrodzntRj3O0w1P/usv2y83s1T2Ou5xQAa6Jn68ghOHeSPcuJczsTuBvICzDCVyssMRq3rXKkA13fhsh+OC9xPDgkp5BWPDmzXHhozbZUOQbFFap+2JnETKf1wF7S1LGQDyFaGjN7EFC19IHJD0buFTS1WZ2yRDK6Iww7jk4o8B7JM2X9GTC8pTdlifNC+f9VeDjkhYp8FxJcwlRTJ+usI7zNvHvhe1xBTO7mbBA0X8mVOwPEMJZ/wU9jIN6hzzve62S9Ax3LmlH4F+Bj9jUVfx6hiLvwX8AE8B7Y37fTFgStH2t18cBeBEi2DYZ/kp2zgjixsEZBc4mDPDeQmixTpmAVSCc9+cIffY/IQwynw48MbZ8DyaMTdxB6Mb5NLBt5vSXE1Y2W535LsIgcDe6hjwveK2i9At3/nzC/I9TsqqluK9fKPIpmNkm4M3A0YSw3m/rOGYRcDHB4PwH8CUzu2yA8jgVw2MrObOKpNuAd5rZxbOdF8dxJnHPwXEcx5mCGwfHcRxnCt6t5DiO40zBPQfHcRxnCm4cnJFC0kmSvjHb+cgiaYEkUwxUpwLhuyU9JaqIGjOTy8GJ+XzqbOfDGS18EpzjlKRI+G4z+wNxtvOoE+MzOc4WuOfgODVFHrLb6YMbB2erkXS8JkNl3yTpTZl9R0u6QtJnFUJy3xpjArX3L1QIf/2gpIuYGgcoe51XSFqjEJr7jzH89JGZ/dvG6/xB0l2S/lnSEzuO/YCkuyWtk/RXmWNfJ+mXkh6QtFrSSX3ykQ3f/SttGT7b4rU6u6J+Kunj6h2e+x2SblcItf3RWLZXdbn2ixTCeDcy296kEIIbSQdK+g+F2E7rJH1R0pxMWpP0Hkk3Azdntu2Xdx8yZToq3uM/SvpIZn/psOnO6OLGwRkGvwdeRpit/DHgG5L2zOx/EbCSUPF/Bjg9hmOAMDv6mrjv40DeUpx7xLR7x7SnKcQZAvgUIeT3AcB+Mc2JHcfuFLcfA5wqaZe472HgHYQw1q8D3qUQl6gvZrZ/JnT2f4/l7DWzumt4boWQ3l8iBOnbM5PHbte7Mub1lR3nPTt+bgL/jXCPXkIIbf7ujtMcRngm3UKJF7kPLyXM0D4IOFGT4UFKh03vVkZnRDAz//O/of4B1wFL4+ejgVWZfdsRYhHtQQjwNgE8KbP/bOAbPc77ii7pzwc+Sgh38TDwtMy+lxCC9bWPfRQYz+y/G3hxj2t9Hjglfl4Q8zwev/+UMKs7m/6l8XxP73PM32fSvxv4Ufx8InBOxz3aBLyqR94+AZwRP+8Qy71vj7TvB76b+W7AKzvSGLBfifswP7P/KmBZ/Lyy/dw7zvE2QoTb7Lb/C/zDbP9W/a/3n/c5OluNpHcQWo0L4qbt2bJ7aHNYajN7JDoN7TT3mtnDmbS3s2WU1k66pd+LsBjNdsA1k04JArJqoXvMbCLzPRsi+0UEz+PZhFb9toSFb3KJXSfnA0eZ2e/6JO0anjvmf3Nk1XiP7ulznrOBf5f0LkJcpGvN7PaYl6cT4kwtIdyPcYJnlmU1PSh4H3qVIzdsembbOGGBIWdE8W4lZ6uQtC9h0Z3jgLlmtjMhbLb6HhhYB+wSux3aPCXnmG7p7yAE4nsUeJaZ7Rz/drLiSpyzCUH09jGznQhhu3PLEMc0vkdYLOeHBa/VyTrCinLZc87tldjMbiIYxUPZsksJ4MvAb4FFZrYjYc2LznL0m/k60H2ItMOmd9t+eea57GyhK+5dBc/rzAJuHJyt5UmEymY9QBzkfXaRA2NrdwXwMUlzFNZGeEOBQ9vpX0aI1PotC4vyfIUQqXS3mJe9NbmMaB47ABvMbKOkAwmVbhHOAH5rZp8pmL4bFwBvkPSf4uDxSeRXyGcTFvp5OVu27HcgRKV9SNIzgbIV8KD3AQYMm+6MJm4cnK0itmL/kRDO+S7gOcC/lTjFXxIGRzcA/0AIVd2POwlrKNxBWFP5v5jZb+O+DxFWLfuFwipoFxMGTovwbuBkSQ8SxgDOL3jcMuBNHYqllxU8FgAzuxH4r4R1oNcRwmPfTVjPoRfnEFaru9TM/pjZ/kHCPX2QYCy7rY3Rj0HvA2xd2HRnxPDYSk5lkPQKwmD1/Ly0VUbS9sB9hK6hW2c7P049cc/BcUYASW+QtF0cT/kscAOT60Y7zozjxsFxRoOlhC6XOwirry0zd+udWcS7lRzHcZwpuOfgOI7jTCGJSXC77rqrLViwYLaz4TiOUymuueaaP5rZvG77kjAOCxYsYMWKFbOdDcdxnEoh6fZe+7xbyXEcx5mCGwfHcRxnCm4cHMdxnCm4cXAcx3Gm4MbBcRzHmYIbB8dxHGcKbhwcx3GcKSQxz8HJ5/wVq1mz4ZGpOyTeuP+e7LfbDjOfqWnmtj8+zHd+uRa6hIjZfacncOSL9p2FXE0/5171B+6479GpOyQOO2Avnjqv6PpHo81NdzzAj369ruu+p+22PUsP6LoMd6VptYyv/ftt3P/Ips3bDvqT3dl/n52Hfi03DjXgkU0T/N0F1wOgjiVkzGD9gxv5X29+7izkbHr5xi9u56tX3Nq1zAAHL96DeTuktaTA/Y8+zvHfuQHo/qzve2QTJy8ttBbTyPPly3/P9391R9dyjgneuP9eqHNnxfn9+of4+A9uAiaf7247PsGNgzMYj0+E2vDE1y/mr1+6cIt9f/qpS9k0kWbwxU3NFrtstw2/PPHgLbafd/Uf+NC3b+DxZmuWcjZ9bJoIZfr4Yc/m7S/e0jM68JMXJ1XmTRNNnrnHDvzo/S/fYvsXLrmZz130O1oGjbRsA4/F53va21/Awc/aY1qv5WMONWCiFX5QjbGpb0pjTDRb6VQYWSZa1rXMY7HJ1WylZxTbZRrv8awnmumUudmyzc8yS/uZTyT4u24/326/62HjxqEG9PtBjY+JiQQrSYBms7txGG+0K4/0yr25IdCj0kzJIE60bPOzzNJ+5imVtc2EGwdnmEzktCZTfIkgVh5jU3/ijbgtRY+pTg2BZg/PcHwsXeM/6RlOf9XtxqEG9KswUjYOLetfeSTU/b6ZzQ2BHi3qZkKLezVb1rPBA0HZkxoj160k6RBJKyWtknR8l/3bSjov7r9S0oK4fa6kyyQ9JOmLmfQ7SLou8/dHSZ+P+46WtD6z753DKWp9afapMMYb6RqHiZzKI8U+6VZfz2GMZkJjDr3GlGrhOczASHuuWklSAzgVeDWwBrha0nIzuymT7BjgXjPbT9Iy4NPA24CNwEeBZ8c/AMzsQeCAzDWuAb6TOd95ZnbcwKVytmCyn7J7F0uKLxGEbqP+nkN65c7rQkzpWTdbxpxtGlO2T3YbplPWNv3EJcOmiOdwILDKzG4xs03AuYTF0LMsBc6Mny8ADpIkM3vYzK4gGImuSHo6sBvw89K5dwrRT8EynnC30kSPAelGDVqW3RoCwUtMx1uqtecwIsZhb2B15vuauK1rGjObAO4H5hbMwzKCp5B9kn8h6XpJF0jap9tBko6VtELSivXr1xe8VD3Jk7Km2L0CsU+6W1da0i3LOnkOrb7dhil1obWpm1ppGXBO5vv3gQVm9lzgIiY9ki0ws9PMbImZLZk3r+sSqE6ktp5Dy3p0pcWWZYKVR7NPQyC1Z93LM5yUKqfX6Bk1tdJaINt6nx+3dU0jaRzYCbgn78SS9gfGzeya9jYzu8fMHotfvwq8oEAenT70a22k1prM0kvN0q48Uqoo27QNXj08B5/nMJ0UMQ5XA4skLZQ0h9DSX96RZjlwVPx8OHBpRzdRL45gS68BSXtmvr4R+E2B8zh96NfaSK01mWWix4B0ymql/vMcxpJ61s0enmHaYw7hNzsTYw65aiUzm5B0HPBjoAGcYWY3SjoZWGFmy4HTgbMkrQI2EAwIAJJuA3YE5kg6DDg4o3R6K/Dajku+V9IbgYl4rqO3onwOk63J7hXlWJLdKxAqj20avSuPlCrKNnnzHFKqMHtLlRMeU+rzLg+bQoH3zOxC4MKObSdmPm8E3tLj2AV9zvvULttOAE4oki+nGH3nOSTtORhP2MbVSm3GE4ujVesZ0jMwz2EUBqSdaaavWqmRuFqpR/dKe39q5KqVEvISJ/LUSgn+rkdtzMGpOLVVKzVz1EoJlruvWimx2fC5nkNChrDNqKmVnIrjaqUtGa9By7JXX3xKxiEvPEpKZW3jnoMzVGqtVuojdUy5ZVmLqKw9PMOUQ7LPpFrJjUMN6O85pBxbqc7zHLp3p6VU5t7rOaQ/puSegzMU+rU20vYcahxbqYcyLSXxQS3VSn0mOQ4bNw41oP88BzGR4sIGuFqpk/Q8B1crTSduHGpAnec5uFppkpTGHFoto2W9ywmpPt/gLanLMrDDxo1DDeg75tBIp8LoxNVKW9JIaLGf9op2dVQrzYTXAG4cakFt1UrNvNhK6ZW7r1opoYZA/5ngYVuaarTuXWnTgRuHGlBErVQsTmK1yPUcEqw8Jvo0BFIac+jrIaWsRnPPwRkmeWolgATfo/Ai9ZvnkGCh8+c5pNGV1uwjskh9zME9B2do5M2QDmnSqDSy9HqRJCXVis6St55Dy8JgbtVp/177r+eQ3m+6l8hiOnDjUAP6aaNTDV9tZn1fpFTDhjRbLSQY6/esE+hCzPOQIFHPoemegzNEinkOab1I7eL0epFSC1/dple8IUhr5nDefA5Io5yd+JiDM1T6aaNTHZztF6a8vT01gwi9Zw1DWi3qQmqlBMrZSbPVmpG1HMCNQy3o19poNNJ8kfqFKW9vT7Vl2SuccyOhhoB7DtOPG4ca0E8bneqYQ16YgVQDDvb1HDZHK61+d1reTHBIdZ6Djzk4Q6Sv55CoWikvQNn4mJJoQXfSK94QpNWi7uc5jI0JydVKW0uhq0g6RNJKSaskHd9l/7aSzov7r5S0IG6fK+kySQ9J+mLHMT+N57wu/u3W71zO4PRrbSTvOTTqplaqx5hDv2CSkFYcqSwj5TlIagCnAocCi4EjJC3uSHYMcK+Z7QecAnw6bt8IfBT4YI/TH2lmB8S/u3PO5QxInqSznSYlcsccGomqlfpIHVNSK/ULJglpzQbPMmpjDgcCq8zsFjPbBJwLLO1IsxQ4M36+ADhIkszsYTO7gmAkitL1XCWOdzrop41ONXx1rdVKPSrMpDyHPmolCL/rFMrZyajFVtobWJ35viZu65rGzCaA+4G5Bc79tdil9NGMASh0LknHSlohacX69esLXKq+FBpzSKz/3dVKU0lp5nDe803Wc2iOlucwXRxpZs8BXhb/3l7mYDM7zcyWmNmSefPmTUsGU6GfNjr5MQdXK20mLc+hv2eYUhypLM0eS6NOB0WMw1pgn8z3+XFb1zSSxoGdgHv6ndTM1sb/DwJnE7qvBjqX05/+8xwSVSv1iU4atifasiygVkrBS6yt5zBiaqWrgUWSFkqaAywDlnekWQ4cFT8fDlxqfWJASxqXtGv8vA3weuDXg5zLyaeWaqUcNUvSYw458xxSeNZ5nuH4mJIwgp3MpFppPC+BmU1IOg74MdAAzjCzGyWdDKwws+XA6cBZklYBGwgGBABJtwE7AnMkHQYcDNwO/DgahgZwMfCVeEjPczmD4WqlqdQ5tlIKz3pyHkuP33UjZc9hRIwDgJldCFzYse3EzOeNwFt6HLugx2lf0CN9z3M5g9Hfc0hcrdRH6phqyzJvzCGFZ53vOaQ6pjRaaiWn4hSbIZ3Wi1RsnkNaZYb2PIc8L7H6HpPPc5h+3DjUgGKxlapfYWRxtdJU0vIcaqxWcuPgDIt+2uiUFCxZaq1W6tOaDmmqX+7aqpWao6VWcipOP210SgqWLPmeQw3VSu3xpQQaAoXUSok+X/ccnKHRT62U0sSoLO1uMlcrTZKm59B7fCW1Bg/Ed3mEJsE5FadfayOlYGxZfJ7DVFLyEguplRLwkDpxtZIzVPopHNL1HPqrWdIdc6iJWqnZ3zNM2nNw4+AMi36tjZSCsWXptxgMRLVSki3LuqiV2ut19PaSUjCCnfiYgzNU6uw59BtrSaGS7KRQbKUEyl1btdKIxVZyKk7/MYd0WpNZcj2HRqJjDn1kyynNhne10vTjxqEG9NNGtyuM1LpY+i1AD4mrlWo1z6E+aiUz69ttOGzcONSAvp5DQgqWLPljDum2LHPHHJrVN4rtZ9ernkwxtlJeV9qwceNQAyZaxlivSlLptCaztF+kfuVOzSBCf7XSWFKeQ4vGmOi1gvBYgp7DRM5veti4cagBtVQrNWs65lAjtVK/7pUUYyu55+AMnWIVxkzmaPppWf6AZSuBSrKTImqlZgJrZ7VyBmYbYyIx27D5ufmYgzM0+o05jI0JKUHPIXfAMvRJp7bIYLExh+qXuZaeQ443PGzcONSAvHgsKcr+Juc59K8oEyt2rWIr5XkOKXSfZZmc+OfzHJwhUcsXKW/MIaFQEm1aLcOs98Q/Sck867zJYCk3eNxzcIaCmRV4kVKU/bWQeis7UhqcbbO5K62Pl5iKhLfZzGvwjCXRfZYlb4GjYePGIXHa9UDtPIcC3lI7XSrkdaVBOpP/csccElSjjaTnIOkQSSslrZJ0fJf920o6L+6/UtKCuH2upMskPSTpi5n020n6f5J+K+lGSZ/K7Dta0npJ18W/d259MetLkdZGkoN3BQYsIY3B2TYTOWtYQEKeQ58V7yDdBg+MkOcgqQGcChwKLAaOkLS4I9kxwL1mth9wCvDpuH0j8FHgg11O/VkzeybwPOBPJR2a2XeemR0Q/75aqkTOFhRpbaT6IvVSKsHkoF4KFWWb4p5D9ctcS7VSjgJv2BS5yoHAKjO7xcw2AecCSzvSLAXOjJ8vAA6SJDN72MyuIBiJzZjZI2Z2Wfy8CbgWmL8V5XB6UKS1MT6mBGMrFfQcEqgo2+SFDIFJCW/VKSKyaBlJzWXJW8Bq2BQxDnsDqzPf18RtXdOY2QRwPzC3SAYk7Qy8Abgks/kvJF0v6QJJ+/Q47lhJKyStWL9+fZFL1ZIi2uhGI43WZJZ+k8EgTbVSXphyiJ5DAg2BImolSGPCX5uRHHOYLiSNA+cAXzCzW+Lm7wMLzOy5wEVMeiRbYGanmdkSM1syb968mclwBSmijU5TreSeQzfSGXPIVyu106XC5vHDEVpDei2Qbb3Pj9u6pokV/k7APQXOfRpws5l9vr3BzO4xs8fi168CLyhwHqcHtR1zyJU6JqhWKtDtMN6oiVopxec7gp7D1cAiSQslzQGWAcs70iwHjoqfDwcutZy4BJI+QTAi7+/Yvmfm6xuB3xTIo9ODWquV+s4KT7dlWYt5DgW7DVPoQmsz02ql8bwEZjYh6Tjgx0ADOMPMbpR0MrDCzJYDpwNnSVoFbCAYEAAk3QbsCMyRdBhwMPAA8BHgt8C1MezuF6My6b2S3ghMxHMdPaSy1pLaeg55aqV2yzKhyqNWaqU+K97BpIFMqdEz02qlXOMAYGYXAhd2bDsx83kj8JYexy7ocdquT9bMTgBOKJIvJ5/CaqUEKowsPjxaXP8AAB0dSURBVObQnZTUSttuk2/8U3y+o6RWcipMkdZGmp5DTrdDwi3LXLVSAs+6qFopBUPYpllgkuMwceOQOEW00eNjY0l1r4B7Dr1IZ8yhhmqlEZzn4FSYeo851EytFFuW/ZaRrE1spSSfb35gxWHixiFximijQ5Cy6lcYWfI9h3RblrmeQwJeYmG1UkK/6yKe4TBx45A4tfUcmgXVSgmVu5BaKZHZ8HX2HPqNtQwTNw6J42ql7ozXtGWZklqpULdhAl5SG/ccnKFSa7VSzmQwSKvyKLIAfTJqpWaOWqmRnuCg6Yv9OMOkmOeQRmsyS67nkGLl0SzWEEjhWRdVK6VQ1jbuOThDpYg2Ok3PoX/lkWKfdNEuxBS60iZyw6MkaPx9EpwzTIrNc6ifWilFHXwRqWM6nkP9QrJPFPAMh4kbh8QpWmGkFKAM8mMrpek5FAuymIJBLKpWSqGsbTZ7Dj7PwRkGRfopU12Mvb/nkJ5aqZhsOY3Z8IXVSgn9rn3MwRkqRbTRaY459O92SNNzqJFaKTe2Uuw2TMAQtnG1kjNUCnkOKaqVckI6pxi1s5BsOREvsc6eQ0NuHJwhUKS1kabnYH3HWdoVaApdLG3qolYys3pKlVvGmPrHzhombhwSp5jnUEO1UoqVR7OYbLnqremiIWEgMbVSjshi2LhxSJwi2uhkPYfaqZXy1SwpjDkULSckZvxzGjzDxo1D4hTRRtcxtlKt1UoVf9blPIdqlzVLCCbpxsEZEkW00Y2xMcygldKLlDdJSglWHjVRK02Ws4BaqeJlzdJstWZsjgMUNA6SDpG0UtIqScd32b+tpPPi/islLYjb50q6TNJDkr7YccwLJN0Qj/mCFN5WSU+WdJGkm+P/Xba+mPWl6DyHbNoUyPMcxsbEmFKrPIoHWTSrbrlr6znkKLSGTa5xkNQATgUOBRYDR0ha3JHsGOBeM9sPOAX4dNy+Efgo8MEup/4y8DfAovh3SNx+PHCJmS0CLonfnQEpqlYKaev1IqUm4W2XpV+xU+iLLzoTHCYH6VNgFMccDgRWmdktZrYJOBdY2pFmKXBm/HwBcJAkmdnDZnYFwUhsRtKewI5m9gsLTZh/AQ7rcq4zM9udASiijZ4cnE3jRWq1DLP8RVFSG4hvxxtSn2fdSMBLLOQ5JFDOTkZRrbQ3sDrzfU3c1jWNmU0A9wNzc865psc5dzezdfHzncDuBfLo9KCINjo1z2FzV1pO/+x4IktmtsmLNwSJeA4Fg0lCtcvZySh6DrNG9Cq6Pl1Jx0paIWnF+vXrZzhn1aFIayM1WWfR0MaNRvUnhGVpFlCzpLDOQdFgklDtcnYycmMOwFpgn8z3+XFb1zSSxoGdgHtyzjm/xznvit1O7e6nu7udwMxOM7MlZrZk3rx5BYpRT4q0NlILXz1RYA2L9v7UKo9aeA51ViuNmHG4GlgkaaGkOcAyYHlHmuXAUfHz4cCl1kcOEbuNHpD04qhSegfwr13OdVRmuzMARbTRtfUckhtzMMYb+eMsUO3xpSJjDu1dqfymob006swZh/G8BGY2Iek44MdAAzjDzG6UdDKwwsyWA6cDZ0laBWwgGBAAJN0G7AjMkXQYcLCZ3QS8G/g68ETgh/EP4FPA+ZKOAW4H3jqMgtaVItrozWMOifS/Fw1tnKJaqR6eQ75aSVIScaSyNHPihQ2bXOMAYGYXAhd2bDsx83kj8JYexy7osX0F8Owu2+8BDiqSLyefQpLORvVbk1mKhCkP+1PzHPpP/IOM51DhhkARzwHSiCOVJS9M+bAZ6QFpZ+spNuZQ/dZkluKeQ4qVR7GGQJWfdZGZ4BBng1fYCHaSF6Z82LhxSJxaqpUKSB3b+5PrdqiTWqmAZ1jlcnYyMYID0k6FqbVaqcBYS5W7VzqpzZhDQeM/3hirdDk7cc/BGSrFwkgk5jkU7XZoJDbm0Mz3EpNSKxUx/gk93yLGf5i4cUicItro1MJXFx1zSCF8dZbaeA4F11JOUq3kxsEZFkW00Zs9h0S6WIqqlVIIX52l2WoVak1Dtb3E2qqVmq5WcoZIEW10XdVKjcSWRy3mOVR/fKmUWqnC5ezEPQdnqEy0rG9EVkijNZml3ZWQtxB7Q/WrPNoNzyp7iUXHlMZS8xxcreQMk1rOc2gW8xzGG6lVHvXyHIoILVKb5+DGwRkaYbnMvL736mvfs9Q6tlIt1ErtAem8sqYnOPBuJWdotFrFKkmodmsyS9NKtCwTKTOUUyu1Kr1MaPhf5PlWuZydtNxzcIbJRAEFSwohFbIUHbBMz3OoS2ylYlLW5NRKMxx4z41D4pQZc6hyV0OW5uYxh/zutKQqjyKy5QQaAqXGHBL5TYOPOThDpswM6SpXGFnq6znktyxTmA1fZkypyh5SJ6O4hrRTYcp5Dmm8SEXDK4wnNs+hWSCkcwpxtCaKeoaphUdxz8EZJsWisla/wshSNLxCIzGpY13iaG32HHIndybWbVhgTGmYuHFInFp7DjWb51BuTkt1PaZyYw71er7DxI1D4hRpbWwec2hWt8LIUtcxhzLPuspGsdZqJTcOzrBoFlCwNBrVrzCyFF0MJjW1Ul1mw282/jlhYVJSK7Vahln+xL9h4sYhcYpoo12tlAbFxhzibPgKj7U0W8aYCsTOSshz2NyVNmrzHCQdImmlpFWSju+yf1tJ58X9V0pakNl3Qty+UtJr4rZnSLou8/eApPfHfSdJWpvZ99rhFLWe1HLMIXaPFVtDOo2WJbS9xBy1UiLzHIpIOlMacygq3x0m43kJJDWAU4FXA2uAqyUtN7ObMsmOAe41s/0kLQM+DbxN0mJgGfAsYC/gYklPN7OVwAGZ868Fvps53ylm9tmtL55TT7VSUTVLOpUHlPMSq9wQKDow2xgbq7SHlGXz0rcjNuZwILDKzG4xs03AucDSjjRLgTPj5wuAgyQpbj/XzB4zs1uBVfF8WQ4Cfm9mtw9aCKc3RV6k9u4qVxhZCquVEup2gBqplZrFBmbdc9g6ihiHvYHVme9r4rauacxsArgfmFvw2GXAOR3bjpN0vaQzJO3SLVOSjpW0QtKK9evXFyhGPSmiYJGU1OBd8TGHMczCYF8KFHnW7UHcKhvFZquV6xVC8ByrXM4sReW7w2RWB6QlzQHeCHwrs/nLwNMI3U7rgH/sdqyZnWZmS8xsybx586Y9r1WluAuezotUWK2UkEqr1TJaVmwBnDFVuwuxqKQzpQZP0aVvh0mRK60F9sl8nx+3dU0jaRzYCbinwLGHAtea2V3tDWZ2l5k1zawFfIWp3VBOCUq9SMn0z4Zy5BU7BVlnm6JhykOaakt469jgGVXP4WpgkaSFsaW/DFjekWY5cFT8fDhwqZlZ3L4sqpkWAouAqzLHHUFHl5KkPTNf3wT8umhhnC0po41O6UVqh65WAR08pBGNtkzLsuoD8bVUKzVHUK1kZhOSjgN+DDSAM8zsRkknAyvMbDlwOnCWpFXABoIBIaY7H7gJmADeY2ZNAElPIiig/rbjkp+RdABgwG1d9jsFKaONHm+MJfMiFVn0BtLyHMq0LMcrHq20lFopgWcLGbXSDM5zyDUOAGZ2IXBhx7YTM583Am/pcewngU922f4wYdC6c/vbi+TJyaeMwiEpz6GEmgXSGHMo07JsNKrdF19uzKH6zxZGV63kVJQy2uiUBu+Kew7pzO8o07KsuoS32WoVbvA0W4YlsFToqI45OBWltp5DyxhvFOuThkQ8h5LPusoGsciKd5BWWJhRVSs5FaVsP3QKLxEMMOZQ4f73NuWedbX74ouseAdpBZR0z8EZKmUVLCm8RDCpVspjcp5D9bvT6qZWKlLOtDyHYmHKh4kbh4Qp25pMoQUNrlbKo/pjDsUGpNsGpMplbTO5NKobB2cIlFKwVLzCyFK08khrzKF4y7JRcfHBRMEB6bQ8B1crOUOklIKl4vLGLPVUKxVvWTYSmOdQtJyQRrfhyK7n4FST2qqVmsVn0EIansNECS8xNASqW+aixj9Nz8HVSs4QcLVSf1IIX92mWaJlWfWZw6U9hwp7SW1creQMlVqrlQpOBoO0Ko+iKp4qNwQmCqx4B5OGssplbeNqJWeolFYrJfASQT3VSkUXOIJ2Q6C63lIt1UruOTjDpKyCJYWXCEqolZKaJFX8WVfecyi42E+aYw5uHJwhUEYb7bGVqk15z6G6Za6lWmnzu+wD0s4QKK1WSqDvHdqVR83USiWedfU9hxqrlVzK6gyDcus5VLvCyFJLtVKJlmVjbKzSDYHynkN1y9rGxxycoVJOrZTOgHTh2EoJVh718RyKeIYpdRu6WskZInWKt5OlaEjnJNVKBaOVVrkf3uc5zAxuHBKmfLyd6r9EUDykc7tlmUblEZ71WM662ZCA59AsGFspqXkOrlZyhkh5z6G6rckszYLdDo0EKw9XK02SlFqpVXxMaVi4cUiYOq0OlqXMGsPt9FWndmMOPs9h2ilkHCQdImmlpFWSju+yf1tJ58X9V0pakNl3Qty+UtJrMttvk3SDpOskrchsf7KkiyTdHP/vsnVFrC9ltNEpjTk066hW8thKU0hKrTSK6zlIagCnAocCi4EjJC3uSHYMcK+Z7QecAnw6HrsYWAY8CzgE+FI8X5s/N7MDzGxJZtvxwCVmtgi4JH53BqCMNrqR1GI/9VMrlWlZVtlzMLPaqpUkGBsl4wAcCKwys1vMbBNwLrC0I81S4Mz4+QLgIEmK2881s8fM7FZgVTxfP7LnOhM4rEAenS6UGnNo1NlzqH65myX6pMOEx2p6S+1HVTvPoaC3NEyKGIe9gdWZ72vitq5pzGwCuB+Ym3OsAT+RdI2kYzNpdjezdfHzncDu3TIl6VhJKyStWL9+fYFi1I+6qpWKjzmkF5gtdc+hbAwpSKfbcCbHG2B2B6RfambPJ3RXvUfSyzsTmJkRjMgUzOw0M1tiZkvmzZs3zVmtJrVVKxUM6ZyW5xBX/SvSEKiwl1hWlQWpSJWLhYQZJkWuthbYJ/N9ftzWNY2kcWAn4J5+x5pZ+//dwHeZ7G66S9Ke8Vx7AncXL46TpaxaqWXQqmilkWWi8DyHtCoPqIPnUG7FO0jF+I+m53A1sEjSQklzCAPMyzvSLAeOip8PBy6Nrf7lwLKoZloILAKukvQkSTsASHoScDDw6y7nOgr418GK5pTRRm92wa0+L9LYmJAS6XYooWZpq5Wsgs+6XDlTGnMoJrIYJuN5CcxsQtJxwI+BBnCGmd0o6WRghZktB04HzpK0CthAMCDEdOcDNwETwHvMrClpd+C7YcyaceBsM/tRvOSngPMlHQPcDrx1iOWtFeU8h0llxzaNnMQjTpkXKRUJb1nPAcLg7gwG+RwKm8vZqJtaaeY9h1zjAGBmFwIXdmw7MfN5I/CWHsd+Evhkx7ZbgP17pL8HOKhIvpz+lF3PAarfymq1LFR6BV+kVAbi25WHCoTPyM4cboxVqyUw0JhDAs93ojmaaiWnopTRRm8enK14/3u7W6y451DtCWFtioYph2rPHK61WmmG3Tw3DglTRhs9uWRmtV+kMmHKQ7pUPIfiXWlVblHX1nMYUbWSU1GaLSsUpRMmo3lWvaKc7Hsvlr6RiIS3jOdQZS+xzNhKlcvZSXiXZ/aabhwSppTnkEgrq10R1M9zqMezLiWyUHXL2UkQWbjn4AyJMgqHVCaElR9zSMM4FI03BJOGs1VFKWuJbqWxMTGmapazk2ZrZiOyghuHpAmL3hR7xKlMGCozYNlOl0LLsllCzZKG51Dwd52I4KDZahWa2DlM3DgkTLl+6DTiDJVpWbbTVd0gQv3GHMoMvtft+Q4LNw4JU0bBUmV5Y5aJZvE+6Xa6qhtEKNeyrLIyrUwwSYiTHCtoBDspM6Y0LNw4JMwgrckqVhhZyix6A6HboYot6E4G8hwqaBTLLnrTaCiJeQ7uOThDZRAFSxUrjCwTJfuk0/Ec6jbmUK/wKE2f5+AMk8E8h2q/SKXHHJJqWZZTK1WxITBR0jP0MYfBceOQMEHBUlzVAdWsMLLUVq1UO8+hhmolNw7OsBjIc6h4/7urlfKZHHOonsdUW7VS0z0HZ4gMomCp+otUJrxCO13dWpZVXuRoILVSEs+32AJWw8SNQ8LUWq1Uotuh6gYRyrUsK61WGshzqPZvGtrRDnxA2hkStVQr1XaeQ/GW5eQ8h+qVu6xaqZHIPIcycdKGhRuHhKm1WqnwPIc0Wpa1USs1S3qGjTTGHEZ1DWmnopTRRtdarZRAy7J2aqXCUtY01EqzsYa0G4eEqbXnUGqeQ7XLDK5W6kUqarSR9RwkHSJppaRVko7vsn9bSefF/VdKWpDZd0LcvlLSa+K2fSRdJukmSTdKel8m/UmS1kq6Lv69duuLWU8Gi61UvQojS3m1UhoD0gOplSpY7rJqpZQWc5ppz2E8L4GkBnAq8GpgDXC1pOVmdlMm2THAvWa2n6RlwKeBt0laDCwDngXsBVws6enABPABM7tW0g7ANZIuypzzFDP77LAKWVcGUbBUvYulvFopjQHp2sRWGsBzeLxZfePQbI6mWulAYJWZ3WJmm4BzgaUdaZYCZ8bPFwAHSVLcfq6ZPWZmtwKrgAPNbJ2ZXQtgZg8CvwH23vriOFkGUbBUscLIMsg8h6qXGcqOOcTw7BVsCAykVkrg+U6M6DyHvYHVme9rmFqRb05jZhPA/cDcIsfGLqjnAVdmNh8n6XpJZ0japVumJB0raYWkFevXry9QjPpRRhudzphDaCWWaVkm0e1QomXZqHBDYGIAz7CK5exkZMccpgtJ2wPfBt5vZg/EzV8GngYcAKwD/rHbsWZ2mpktMbMl8+bNm5H8Vo1ya0gnolYaYJ5D1csMNVQrlRhTqqKH1MmoqpXWAvtkvs+P27qmkTQO7ATc0+9YSdsQDMM3zew77QRmdpeZNc2sBXyF0K3lDMAga0hXscLIMsg8h6qXGeKYQ4lIpVBN8UHZ9RxS8BxaLaNlo7mG9NXAIkkLJc0hDDAv70izHDgqfj4cuNTMLG5fFtVMC4FFwFVxPOJ04Ddm9rnsiSTtmfn6JuDXZQvlBMq0NmqtVkqgZVkntZIEYyUW+6l6t2HTyhnEYZGrVjKzCUnHAT8GGsAZZnajpJOBFWa2nFDRnyVpFbCBYECI6c4HbiIolN5jZk1JLwXeDtwg6bp4qQ+b2YXAZyQdABhwG/C3Qyxvrai151BiBm3Vywz1UiuVqSRT8BzKhikfFrnGASBW2hd2bDsx83kj8JYex34S+GTHtiuArk/YzN5eJE9OPuXGHGKFUfFWtKuV8tmsVqpgucsOzKagVior3x0WPkM6Ycpoo9PxHGqqVhpAmVZFozhRcrnMJDyHkiKLYeHGIWHKaKMlJdGKHsRzaFkY9KsyA6mVKugllvccqh9bqd14GcV5Dk5FqaML3hxAzQKTg35VxMxKPeuxMSFVU3xQVtKZhOdQssEzLNw4JMxgL1L1Kowsg6iVoJpdLG3KBhtsp61iQ2CgBk/Fw2f4mIMzVAbRRifhOcTKI6il86myrLPNZoNYotuhql2IE01XK80UbhwSZRBtdAovUhlJJ2QGZyvY/95mMM+hmn3xzRKT/aA9z6F65czinoMzVAZpbaQweFdmMhhkl8ysbtfDxEDPupoNgVqqlUqGKR8WbhwSZZDWxviYKt2Chq3wHCpcgQw+5lA9gzioWskqLDhwz8EZKoNoo1MZcyhbSULVxxzKtyyr6zmUF1kAVLComykbTHJYuHFIlEG00WHJzOq1JrOUmQwGNVcrVdBLHEStBNXuNiwbTHJYuHFIlEG00Ul4DgOoWaDinsMgXmJF184eJLYSVNv4DzKmNAwKxVZKlfOvXs1Xfn7LbGdjWtj8gyoo6YTwIv105Xpe/bnLpytb086dD2xkxydsUzh9u0I9+mtXMadRzbZSexnMMsZhfGyMi266q3LPeu19j/KMPXYonL59T17/T1eUehdGiUcfbwLl3uVhUGvjsPN227Bo9+1nOxvTxv7zd+Kli3YtnP6Yly7k8t9Ve1W9Rbtvz4ufOrdw+hctfDJvft7ebJxoTmOupp/nPWUXXvK04uV+58sW8m+r/jiNOZoeFu2+Pa951h6F07/ymbtx/Zr7K92tBPCSp87lOfN3mtFrqsqj+G2WLFliK1asmO1sOI7jVApJ15jZkm77qulHO47jONOKGwfHcRxnCm4cHMdxnCm4cXAcx3Gm4MbBcRzHmYIbB8dxHGcKbhwcx3GcKbhxcBzHcaaQxCQ4SeuB2wc8fFdg1KeKjnoeRz1/4HkcBqOePxj9PI5a/vY1s3nddiRhHLYGSSt6zRAcFUY9j6OeP/A8DoNRzx+Mfh5HPX9ZvFvJcRzHmYIbB8dxHGcKbhzgtNnOQAFGPY+jnj/wPA6DUc8fjH4eRz1/m6n9mIPjOI4zFfccHMdxnCm4cXAcx3GmUGvjIOkQSSslrZJ0/AjkZx9Jl0m6SdKNkt4Xtz9Z0kWSbo7/dxmBvDYk/VLSD+L3hZKujPfyPElzZjFvO0u6QNJvJf1G0ktG7R5K+m/xGf9a0jmSnjDb91DSGZLulvTrzLau902BL8S8Xi/p+bOUv/8dn/P1kr4raefMvhNi/lZKes10569XHjP7PiDJJO0av8/4PSxDbY2DpAZwKnAosBg4QtLi2c0VE8AHzGwx8GLgPTFPxwOXmNki4JL4fbZ5H/CbzPdPA6eY2X7AvcAxs5KrwP8BfmRmzwT2J+RzZO6hpL2B9wJLzOzZQANYxuzfw68Dh3Rs63XfDgUWxb9jgS/PUv4uAp5tZs8FfgecABDfm2XAs+IxX4rv/GzkEUn7AAcDf8hsno17WJjaGgfgQGCVmd1iZpuAc4Gls5khM1tnZtfGzw8SKrW9Y77OjMnOBA6bnRwGJM0HXgd8NX4X8Erggphk1vIoaSfg5cDpAGa2yczuY8TuIWH99idKGge2A9Yxy/fQzH4GbOjY3Ou+LQX+xQK/AHaWtOdM58/MfmJmE/HrL4D5mfyda2aPmdmtwCrCOz+t9LiHAKcAfwdkFUAzfg/LUGfjsDewOvN9Tdw2EkhaADwPuBLY3czWxV13ArvPUrbafJ7wQ2+v2j4XuC/zks7mvVwIrAe+Fru9virpSYzQPTSztcBnCa3IdcD9wDWMzj3M0uu+jeL789fAD+PnkcmfpKXAWjP7VceukcljN+psHEYWSdsD3wbeb2YPZPdZ0B7Pmv5Y0uuBu83smtnKQw7jwPOBL5vZ84CH6ehCGoF7uAuh1bgQ2At4El26IkaN2b5v/ZD0EUK37DdnOy9ZJG0HfBg4cbbzUpY6G4e1wD6Z7/PjtllF0jYEw/BNM/tO3HxX292M/++erfwBfwq8UdJthK64VxL6+HeOXSQwu/dyDbDGzK6M3y8gGItRuoevAm41s/Vm9jjwHcJ9HZV7mKXXfRuZ90fS0cDrgSNtcuLWqOTvaYRGwK/iOzMfuFbSHoxOHrtSZ+NwNbAoKkTmEAavls9mhmLf/enAb8zsc5ldy4Gj4uejgH+d6by1MbMTzGy+mS0g3LNLzexI4DLg8Jhs1vJoZncCqyU9I246CLiJEbqHhO6kF0vaLj7zdh5H4h520Ou+LQfeERU3Lwbuz3Q/zRiSDiF0cb7RzB7J7FoOLJO0raSFhEHfq2Y6f2Z2g5ntZmYL4juzBnh+/J2OxD3siZnV9g94LUHh8HvgIyOQn5cS3Pbrgevi32sJffqXADcDFwNPnu28xvy+AvhB/PxUwsu3CvgWsO0s5usAYEW8j98Ddhm1ewh8DPgt8GvgLGDb2b6HwDmEMZDHCZXYMb3uGyCC2u/3wA0E5dVs5G8Vod++/b78cyb9R2L+VgKHztY97Nh/G7DrbN3DMn8ePsNxHMeZQp27lRzHcZweuHFwHMdxpuDGwXEcx5mCGwfHcRxnCm4cHMdxnCm4cXAcx3Gm4MbBcRzHmcL/B5KVSFk9v1vWAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "t = np.linspace(.001, 150, 100)\n", "plt.plot(t, hazard(estimates_, t))\n", "plt.title(\"\"\"Estimated hazard function using \\npiecewise hazards \\nand penalizing variance\"\"\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Why do this?\n", "\n", "1) From a Bayesian perspective, this is a way to add prior information into a model. Note that if we have _lots_ of observations, the prior becomes less relevant (just like a traditional Bayesian model). That's good. *We are again using the likelihood to \"add\" information to our system.*\n", "\n", "2) When we have low data sizes, we can \"borrow\" information between periods. That is, deaths in the earlier \"inbetween\" periods can inform the rate in the later \"inbetween\" periods - thus we can do better inference. \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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 }