example_notebooks/CensoredSurvivalData.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Survival function with censored data\n",
    "\n",
    "In this notebook we will explore how not dealing with censored data correctly can lead to bias (in this example, a systematic underestimation of the expected lifetime).\n",
    "\n",
    "In order to do this, consider the ages of 100 vampires who were \"born\" (i.e. turned into a vampire) between the years 1400 and 2018.\n",
    "\n",
    "## Terminology\n",
    "\n",
    "- _born_: in this context we will mean \"turned into a vampire\"\n",
    "- _undead_: the state of being \"alive\" as a vampire, so after they were turned, but before they were permanently killed\n",
    "- _died_: for a vampire, we mean that they were killed permanently, e.g. by exposure to sun rays, burning, stake to the heart, or any other way a vampire can die\n",
    "- _age_: the number of years since a vampire was turned\n",
    "- _lifetime_ or _age at time of death_: the age a vampire had when they died\n",
    "- _censored lifetime_: this can be their age at time of death if they have already died, but if they are still undead, it is their current age"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import lifelines\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "% matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_year = 1400\n",
    "current_year = 2018\n",
    "number_of_vampires = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## We have more data from vampires born recently\n",
    "\n",
    "This can be because of a population increase, or because data collection has improved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "birth_dates_increasing = map(int, np.random.triangular(\n",
    "    start_year, current_year, current_year, size=number_of_vampires))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'frequency')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFARJREFUeJzt3X20ZXV93/H3hyc1SsMgF0MQHLSg4lp1wAmSqlRQEMlqwcRUMZqJ0E60ouLSrJLYpSS2jRofstK6NGOgDCpojBBJlgVHysMixdEZHHnolPAgJSgLhggKmmCBb//Y+5bLcO/cMzD7nLnn936tddbd53f2Pvv7mz0zn7uffjtVhSSpXbtMugBJ0mQZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTG7TbpAkaxzz771PLlyyddhiQtKRs3brynqmYWm29JBMHy5cvZsGHDpMuQpCUlyf8ZZT4PDUlS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuOWxJ3FkrQ9zlt/+6RL2GHe9NIDB1+HewSS1DiDQJIaN1gQJHlqkm8l+W6SG5L8ft9+UJL1SW5K8qUkewxVgyRpcUPuETwIHFNVLwZWAMcnORL4CPDJqjoYuBc4dcAaJEmLGCwIqvNA/3b3/lXAMcBf9O1rgZOGqkGStLhBrxpKsiuwEfinwKeAW4D7quqhfpY7gP0XWHY1sBrgwAOHP2sutW6arrTR9hn0ZHFVPVxVK4BnA0cAL5xvtgWWXVNVK6tq5czMog/YkSQ9QWO5aqiq7gMuB44E9koyuyfybOAH46hBkjS/Ia8amkmyVz/9NODVwGbgMuD1/WyrgK8OVYMkaXFDniPYD1jbnyfYBfjzqvrrJP8L+GKS/wh8BzhrwBokSYsYLAiq6lrgsHnab6U7XyBJ2gl4Z7EkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxgwVBkgOSXJZkc5Ibkry7bz8zyfeTbOpfJwxVgyRpcbsN+N0PAe+tqmuS7AlsTLKu/+yTVfWxAdctSRrRYEFQVXcCd/bT9yfZDOw/1PokSU/MWM4RJFkOHAas75tOS3JtkrOTLBtHDZKk+Q0eBEmeAXwFOL2qfgx8GngesIJuj+HjCyy3OsmGJBu2bNkydJmS1KxBgyDJ7nQh8IWqugCgqu6qqoer6hHgs8AR8y1bVWuqamVVrZyZmRmyTElq2pBXDQU4C9hcVZ+Y077fnNleB1w/VA2SpMUNedXQy4C3ANcl2dS3/R5wcpIVQAG3Ab89YA2SpEUMedXQVUDm+ehrQ61TkrT9vLNYkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkho3WBAkOSDJZUk2J7khybv79r2TrEtyU/9z2VA1SJIWN+QewUPAe6vqhcCRwDuSHAqcAVxaVQcDl/bvJUkTMlgQVNWdVXVNP30/sBnYHzgRWNvPthY4aagaJEmLG8s5giTLgcOA9cCzqupO6MIC2HeBZVYn2ZBkw5YtW8ZRpiQ1afAgSPIM4CvA6VX141GXq6o1VbWyqlbOzMwMV6AkNW7RIEiy9xP98iS704XAF6rqgr75riT79Z/vB9z9RL9fkvTkjbJHsD7Jl5OckCSjfnE/71nA5qr6xJyPLgJW9dOrgK+OXK0kaYcbJQgOAdYAbwFuTvKfkxwywnIv65c5Jsmm/nUC8GHg2CQ3Acf27yVJE7LbYjNUVQHrgHVJjgY+D/y7JN8FzqiqqxdY7ipgoT2IVz3BeiVJO9iiQZDkmcCb6X67vwt4J93hnRXAl4GDhixQkjSsRYMAuBr4HHBSVd0xp31Dks8MU5YkaVxGCYLn94eHHqeqPrKD65EkjdkoJ4u/nmSv2TdJliW5ZMCaJEljNEoQzFTVfbNvqupeFrgbWJK09IwSBA8nOXD2TZLnAPMeKpIkLT2jnCN4P3BVkiv690cBq4crSZI0TqPcR3BxksPphpIO8J6qumfwyiRJYzHKHgHAU4Af9vMfmoSqunK4siRJ4zLKDWUfAd4A3AA80jcXYBBI0hQYZY/gJLp7CR4cuhhJ0viNctXQrcDuQxciSZqMUfYIfgpsSnIp8P/3CqrqXYNVJUkam1GC4KL+JUmaQqNcPro2ydOAA6vqxjHUpAact/72SZewQ7zppQcuPpO0kxvlUZX/EtgEXNy/X5HEPQRJmhKjnCw+EzgCuA+gqjbhMwgkaWqMEgQPVdWPtmpzrCFJmhKjnCy+PsmbgF2THAy8C/ifw5YlSRqXUfYI3gm8iO7S0fOBHwOnD1mUJGl8Rrlq6Kd0I5C+f/hyJEnjNspYQ5cxzzmBqjpmkIokSWM1yjmC982Zfirwa8BDw5QjSRq3UQ4Nbdyq6W/mPKRGkrTEjXJD2d5zXvskeQ3wCyMsd3aSu5NcP6ftzCTfT7Kpf53wJOuXJD1Joxwa2kh3jiB0h4S+B5w6wnLnAP8VOHer9k9W1ce2o0ZJ0oBGOTT0hO4irqorkyx/IstKksZnlKuGfnVbn1fVBdu5ztOS/CawAXhvVd27nctLknagUQ4NnQr8c+B/9O+PBi4HfkR3yGh7guDTwIf65T4EfBw4Zb4Zk6wGVgMceKAjPErSUEYJggIOrao7AZLsB3yqqt66vSurqrtmp5N8Fvjrbcy7BlgDsHLlSsc2kqSBjDLExPLZEOjdBRzyRFbWh8is1wHXLzSvJGk8RtkjuDzJJXTjDBXwRuCyxRZKcj7wSmCfJHcAHwRemWRF/z23Ab/9xMqWJO0oo1w1dFqS1wFH9U1rqurCEZY7eZ7ms7azPknSwEbZIwC4Bri/qr6R5OeS7FlV9w9ZmCRpPEa5s/jfAn8B/GnftD/wl0MWJUkan1FOFr8DeBndcwioqpuAfYcsSpI0PqMEwYNV9bPZN0l2w0dVStLUGCUIrkjye8DTkhwLfBn4q2HLkiSNyyhBcAawBbiO7nLPrwH/YciiJEnjs82rhpLsCqytqjcDnx1PSZKkcdrmHkFVPQzMJNljTPVIksZslPsIbqN7KtlFwE9mG6vqE0MVJUkanwX3CJJ8rp98A93gcLsAe855SZKmwLb2CF6S5DnA7cB/GVM9kqQx21YQfAa4GDiI7iEys0J3H8FzB6xLkjQmCx4aqqo/qaoXAv+tqp4753VQVRkCkjQlFr2PoKrePo5CJEmTMcoNZZKkKWYQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMaN+vB67QTOW3/7pEuQNIXcI5CkxhkEktS4wYIgydlJ7k5y/Zy2vZOsS3JT/3PZUOuXJI1myD2Cc4Djt2o7A7i0qg4GLu3fS5ImaLAgqKorgR9u1XwisLafXgucNNT6JUmjGfc5gmdV1Z0A/c99x7x+SdJWdtqTxUlWJ9mQZMOWLVsmXY4kTa1xB8FdSfYD6H/evdCMVbWmqlZW1cqZmZmxFShJrRl3EFwErOqnVwFfHfP6JUlbGfLy0fOBq4HnJ7kjyanAh4Fjk9wEHNu/lyRN0GBDTFTVyQt89Kqh1ilJ2n477cliSdJ4GASS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOZxZLT4LPkdY0cI9AkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktS4iTyYJsltwP3Aw8BDVbVyEnVIkib7hLKjq+qeCa5fkoSHhiSpeZMKggK+nmRjktXzzZBkdZINSTZs2bJlzOVJUjsmFQQvq6rDgdcC70hy1NYzVNWaqlpZVStnZmbGX6EkNWIiQVBVP+h/3g1cCBwxiTokSRMIgiRPT7Ln7DRwHHD9uOuQJHUmcdXQs4ALk8yu/7yqungCdUiSmEAQVNWtwIvHvV5J0vy8fFSSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUuEk+j2Aszlt/+6RLkKSdmnsEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNW4iQZDk+CQ3Jrk5yRmTqEGS1Bl7ECTZFfgU8FrgUODkJIeOuw5JUmcSewRHADdX1a1V9TPgi8CJE6hDksRkgmB/4O/mvL+jb5MkTcAknlmcedrqcTMlq4HV/dsHktw4aFWPtw9wz5jXOW7T3sdp7x9Mfx+nvX+wSB9/48l993NGmWkSQXAHcMCc988GfrD1TFW1BlgzrqK2lmRDVa2c1PrHYdr7OO39g+nv47T3D3aOPk7i0NC3gYOTHJRkD+CNwEUTqEOSxAT2CKrqoSSnAZcAuwJnV9UN465DktSZxKEhquprwNcmse7tMLHDUmM07X2c9v7B9Pdx2vsHO0EfU/W487SSpIY4xIQkNa6pIEhydpK7k1w/p+3MJN9Psql/nTDns9/th8G4Mclr5rTvlENkzNe/vv2dfb03JPnonPYl1T9YcBt+ac72uy3JpjmfLak+LtC/FUm+2fdvQ5Ij+vYk+ZO+D9cmOXzOMquS3NS/Vk2iLwtZoI8vTnJ1kuuS/FWSfzLns6W2DQ9IclmSzf2/uXf37XsnWddvk3VJlvXtk9+OVdXMCzgKOBy4fk7bmcD75pn3UOC7wFOAg4Bb6E5u79pPPxfYo5/n0En3bRv9Oxr4BvCU/v2+S7V/C/Vxq88/DnxgqfZxgW34deC1/fQJwOVzpv873b05RwLr+/a9gVv7n8v66WWT7tsiffw28C/66VOADy3hbbgfcHg/vSfwt30/Pgqc0befAXxkZ9mOTe0RVNWVwA9HnP1E4ItV9WBVfQ+4mW54jJ12iIwF+vd24MNV9WA/z919+5LrH2x7GyYJ8K+B8/umJdfHBfpXwOxvyD/Po/fdnAicW51vAnsl2Q94DbCuqn5YVfcC64Djh69+NAv08fnAlf30OuDX+umluA3vrKpr+un7gc10oyecCKztZ1sLnNRPT3w7NhUE23Bav0t29uzuGgsPhbHUhsg4BHhFkvVJrkjyS337tPRvrlcAd1XVTf37aenj6cAfJfk74GPA7/bt09I/gOuBf9VP/zqP3nS6pPuYZDlwGLAeeFZV3QldWAD79rNNvI8GAXwaeB6wAriT7tACLDwUxkhDZOxEdqPbrTwS+B3gz/vfnKelf3OdzKN7AzA9fXw78J6qOgB4D3BW3z4t/YPucNA7kmykO5zys759yfYxyTOArwCnV9WPtzXrPG1j7WPzQVBVd1XVw1X1CPBZul1OWHgojJGGyNiJ3AFc0O92fgt4hG5sk2npHwBJdgN+FfjSnOZp6eMq4IJ++stM399Rqup/V9VxVfUSujC/pf9oSfYxye50IfCFqprddnf1h3zof84epp18Hyd9YmXcL2A5jz1Jtd+c6ffQHY8EeBGPPUl1K90Jqt366YN49CTViybdr230723AH/TTh9Dtamap9m++PvZtxwNXbNW2JPs4zzbcDLyyn34VsLGf/hUee5LxW3373sD36PYEl/XTe0+6X4v0cfYihl2Ac4FTluo27LfHucAfb9X+Rzz2ZPFHd5btOPE/tDFvoPPpDv/8X7q0PRX4HHAdcC3dmEdzg+H9dL+Z3Eh/1UbffgLdlQC3AO+fdL8W6d8ewOfpjsFeAxyzVPu3UB/79nOAt80z/5Lq4wLb8OXAxv4/u/XAS/p5Q/eQp1v6v8Mr53zPKXQnVm8G3jrpfo3Qx3f32+NvgQ/T3+y6RLfhy+kO4VwLbOpfJwDPBC4Fbup/7r2zbEfvLJakxjV/jkCSWmcQSFLjDAJJapxBIEmNMwgkqXEGgbQDJHlBPzrod5I8b6vPHlhgmT9I8uoFPvutJL845/1tSfbZsVVLHYNA2oYku44460nAV6vqsKq6ZdG5gar6QFV9Y4F1/hbwi49bSBqAQaCpkORDs+O+9+//U5J39dO/k+Tb/cCCvz9nnr9MsrEfM371nPYH+t/W1wO/vNV6Zp8NcG2SC5MsS/cMi9OBf5PksgXq+3iSa5JcmmSmbzsnyev76duSfCDJVXRjJq0EvtDvZTyt/5p39t9xXZIX7IA/NgkwCDQ9zqIbk4ckuwBvpPuP9DjgYLrxeVYAL0lyVL/MKdWNbbMSeFeSZ/btT6cb/uClVXXVVus5F/j3VfXP6O4C/WB1z+D+DPDJqjp6ntqeDlxTVYcDVwAfXKAP/1hVL6+qzwMbgN+oqhVV9Q/95/f03/Fp4H2j/sFIizEINBWq6jbg75McBhwHfKeq/r6fPg74Dt0QGy+gCwbo/vP/LvBNusG9Ztsfphsw7DGS/DywV1Vd0TetpXvIymIe4dHB8D5PNwTBfL60QPus2cHLNtKN1SPtELtNugBpB/ozumPrvwCc3bcF+MOq+tO5MyZ5JfBq4Jer6qdJLgee2n/8j1X18IB1LjSuy08WWe7B/ufD+G9XO5B7BJomF9KNQvpLwCV92yXAKf3Y8CTZP8m+dE/6urcPgRfQjfq4TVX1I+DeJK/om95Cd6hnMbsAr++n3wRsfbhpPvfTjcsvDc7fKjQ1qupn/cna+2Z/o6+qryd5IXB19zweHgDeDFwMvC3JtXSjWn5zxNWsAj6T5OfohkF+6wjL/AR4Uf/QlR8BbxhhmXP69fwDW52wlnY0Rx/V1OhPEl8D/Ho9+rhKSYvw0JCmQpJD6cZsv9QQkLaPewSS1Dj3CCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLj/h/z05xwZ0m6gQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(birth_dates_increasing, kde=False, label='histogram of vampires born in each year')\n",
    "plt.xlabel('year of birth')\n",
    "plt.ylabel('frequency')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## We know all their ages at time of death - even if they haven't died yet (asked a clairvoyant)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Their theoretical average age is 500 years\n",
    "# this is the mean of the distribution we will draw their final ages from\n",
    "theoretical_average_age = 500\n",
    "\n",
    "true_ages_at_time_of_death = map(int, np.random.exponential(\n",
    "    theoretical_average_age, size=number_of_vampires))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'frequency')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFbdJREFUeJzt3X20XXV95/H3xxCEKlNArjYFQtCFWtq1DHgHQVqHoiKlD4DaUVAbLTNpZxVRp86IOmsNXc50iatKx7FjJxZq6oj4iFCKQIo8yFijNxhCMGAQ0EGySFw81ykV+M4fe9/pNb0PJzd3n0Pufr/WOuvu8zt7n/39nZybz91Pv52qQpLUX88YdQGSpNEyCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknttr1AUM4qCDDqoVK1aMugxJ2qNs2LDhR1U1Ntd8e0QQrFixgomJiVGXIUl7lCTfH2Q+dw1JUs8ZBJLUcwaBJPWcQSBJPWcQSFLPdR4ESZYk+XaSK9rnhydZn2Rrks8m2bvrGiRJMxvGFsE7gC1Tnp8PXFBVRwAPAmcNoQZJ0gw6DYIkhwC/DvxF+zzAicAX2lnWAqd1WYMkaXZdbxH8KfAfgafa588BHqqqJ9rn9wIHd1yDJGkWnV1ZnOQ3gO1VtSHJCZPN08xaMyy/GlgNsHz58nnXcfH6H8x72d115svmX7ckDUuXWwTHA7+V5B7gEppdQn8K7J9kMoAOAe6bbuGqWlNV41U1PjY251AZkqR56iwIquq9VXVIVa0A3gh8tareBFwHvL6dbRVwWVc1SJLmNorrCN4D/Pskd9IcM7hwBDVIklpDGX20qq4Hrm+n7wKOGcZ6JUlz88piSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknquc6CIMk+Sb6Z5JYktyX5o7b9k0nuTrKxfazsqgZJ0ty6vFXl48CJVfVYkqXATUm+0r72H6rqCx2uW5I0oM6CoKoKeKx9urR9VFfrkyTNT6fHCJIsSbIR2A6sq6r17Uv/NcmmJBckeWaXNUiSZtdpEFTVk1W1EjgEOCbJLwHvBV4M/EvgQOA90y2bZHWSiSQTO3bs6LJMSeq1oZw1VFUPAdcDJ1fVtmo8DvwlcMwMy6ypqvGqGh8bGxtGmZLUS12eNTSWZP92el/gVcDtSZa1bQFOAzZ3VYMkaW5dnjW0DFibZAlN4Hyuqq5I8tUkY0CAjcDvd1iDJGkOXZ41tAk4apr2E7tapyRp13llsST1nEEgST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzBoEk9VyX9yzeJ8k3k9yS5LYkf9S2H55kfZKtST6bZO+uapAkza3LLYLHgROr6iXASuDkJMcC5wMXVNURwIPAWR3WIEmaQ2dBUI3H2qdL20cBJwJfaNvXAqd1VYMkaW6dHiNIsiTJRmA7sA74HvBQVT3RznIvcPAMy65OMpFkYseOHV2WKUm91mkQVNWTVbUSOAQ4BviF6WabYdk1VTVeVeNjY2NdlilJvTaUs4aq6iHgeuBYYP8ke7UvHQLcN4waJEnT6/KsobEk+7fT+wKvArYA1wGvb2dbBVzWVQ2SpLntNfcs87YMWJtkCU3gfK6qrkjyHeCSJP8F+DZwYYc1SJLm0FkQVNUm4Khp2u+iOV4gSXoa8MpiSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknquS7vWXxokuuSbElyW5J3tO3nJflhko3t45SuapAkzW3OW1UmObCqHpjHez8B/GFV3ZxkP2BDknXtaxdU1Z/M4z0lSQtskC2C9Uk+n+SUJBn0jatqW1Xd3E4/CmwBDp5nnZKkjgwSBC8E1gBvAe5M8sdJXrgrK0myguZG9uvbprOTbEpyUZIDZlhmdZKJJBM7duzYldVJknbBnEFQjXVVdQbwb4BVwDeT3JDkuLmWT/Js4IvAO6vqEeDjwAuAlcA24MMzrHdNVY1X1fjY2NjgPZIk7ZJBjhE8B3gzzRbB/cDbgctp/iP/PHD4LMsupQmBT1fVlwCq6v4pr38CuGI36pck7aY5gwD4O+BTwGlVde+U9okkfz7TQu3xhAuBLVX1kSnty6pqW/v0dGDzrpctSVoogwTBi6qqpnuhqs6fZbnjabYibk2ysW17H3BGkpVAAfcAvzd4uZKkhTZIEFyT5Ler6iGA9uDuJVX1mtkWqqqbgOnOMrpy18uUJHVlkLOGxiZDAKCqHgSe211JkqRhGiQInkyyfPJJksNodutIkhaBQXYNvR+4KckN7fNXAKu7K0mSNExzBkFVXZXkaOBYmn3+76qqH3VemSRpKAbZIgB4JvBAO/+RSaiqG7srS5I0LINcUHY+8AbgNuCptrkAg0CSFoFBtghOo7mW4PGui5EkDd8gZw3dBSztuhBJ0mgMskXwY2BjkmuB/79VUFXndFaVJGloBgmCy9uHJGkRGuT00bVJ9gWWV9UdQ6hJkjREcx4jSPKbwEbgqvb5yiRuIUjSIjHIweLzgGOAhwCqaiOz3INAkrRnGSQInqiqh3dqc6whSVokBjlYvDnJmcCSJEcA5wBf77YsSdKwDLJF8HbgF2lOHf0M8Ajwzi6LkiQNzyBnDf2YZgTS93dfjiRp2AYZa+g6pjkmUFUnzrHcocBfAT9HM0bRmqr6b0kOBD4LrKC5VeW/bm92I0kagUGOEbx7yvQ+wOuAJwZY7gngD6vq5iT7ARuSrAPeClxbVR9Mci5wLvCeXStbkrRQBtk1tGGnpv895SY1sy23DdjWTj+aZAtwMHAqcEI721rgegwCSRqZQXYNHTjl6TOAl9Ls7hlYkhXAUcB64HltSFBV25JMe//jJKtp74S2fPny6WaRJC2AQXYNbaA5RhCa3T13A2cNuoIkzwa+CLyzqh5JMtByVbUGWAMwPj7udQuS1JFBdg3N+yriJEtpQuDTVfWltvn+JMvarYFlwPb5vr8kafcNsmvotbO9PuU/+J2XC3AhsKWqPjLlpcuBVcAH25+XDVytJGnBDbJr6Czg5cBX2+e/SnOA92GaXUbTBgFwPPAW4NYkG9u299EEwOeSnAX8APjteVUuSVoQgwRBAUdOHuBtd+f8WVW9bdaFqm6iOa4wnVfuUpWSpM4MMsTEiskQaN0PvLCjeiRJQzbIFsH1Sa6mGWeogDcC13ValSRpaAY5a+jsJKcDr2ib1lTVpd2WJUkalkG2CABuBh6tqr9N8jNJ9quqR7ssTJI0HIPcqvLfAl8A/mfbdDDw5S6LkiQNzyAHi/+A5lTQRwCqaisw7bAQkqQ9zyBB8HhV/ePkkyR74a0qJWnRGCQIbkjyPmDfJK8GPg/8dbdlSZKGZZAgOBfYAdwK/B5wJfCfuixKkjQ8s541lGQJsLaq3gx8YjglSZKGadYgqKonk4wl2XvqcQIN5uL1Pxh1CUN15su8b4S0JxrkOoJ7aO5Kdjnw95ONO40oKknaQ814jCDJp9rJNwBXtPPuN+UhSVoEZtsieGmSw2iGiv7vQ6pHkjRkswXBnwNXAYcDE1PaQ3MdwfM7rEuSNCQz7hqqqo9W1S8Af1lVz5/yOLyqDAFJWiTmvI6gqv7dfN44yUVJtifZPKXtvCQ/TLKxfZwyn/eWJC2cQS4om69PAidP035BVa1sH1d2uH5J0gA6C4KquhF4oKv3lyQtjC63CGZydpJN7a6jA0awfknSFMMOgo8DLwBWAtuAD880Y5LVSSaSTOzYsWNY9UlS7ww1CKrq/qp6sqqeohm76JhZ5l1TVeNVNT42Nja8IiWpZ4YaBEmWTXl6OrB5pnklScMx6D2Ld1mSzwAnAAcluRf4z8AJSVbSXJB2D82w1pKkEeosCKrqjGmaL+xqfZKk+RnFWUOSpKcRg0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknquc6CIMlFSbYn2Tyl7cAk65JsbX8e0NX6JUmD6XKL4JPAyTu1nQtcW1VHANe2zyVJI9RZEFTVjcADOzWfCqxtp9cCp3W1fknSYIZ9jOB5VbUNoP353CGvX5K0k6ftweIkq5NMJJnYsWPHqMuRpEVr2EFwf5JlAO3P7TPNWFVrqmq8qsbHxsaGVqAk9c2wg+ByYFU7vQq4bMjrlyTtpMvTRz8D/B3woiT3JjkL+CDw6iRbgVe3zyVJI7RXV29cVWfM8NIru1qnJGnXdRYE0jBdvP4HI1nvmS9bPpL1SgvpaXvWkCRpOAwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes77EWjBjOqeAJJ2z0iCIMk9wKPAk8ATVTU+ijokSaPdIvjVqvrRCNcvScJjBJLUe6MKggKuSbIhyeoR1SBJYnS7ho6vqvuSPBdYl+T2qrpx6gxtQKwGWL7cG4RLUldGskVQVfe1P7cDlwLHTDPPmqoar6rxsbGxYZcoSb0x9CBI8qwk+01OAycBm4ddhySpMYpdQ88DLk0yuf6Lq+qqEdQhSWIEQVBVdwEvGfZ6JUnT88piaTeM8mrqM1/mSRRaGF5HIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST3nlcWSdon3ph6uYVxB7haBJPWcQSBJPWcQSFLPGQSS1HMGgST1nGcNSXsoz97RQhnJFkGSk5PckeTOJOeOogZJUmMUN69fAvwZ8GvAkcAZSY4cdh2SpMYotgiOAe6sqruq6h+BS4BTR1CHJInRBMHBwP+Z8vzetk2SNAKjOFicadrqn82UrAZWt08fS3LHPNd3EPCjeS67GPS9/+BnAH4Ge2z/37R7ix82yEyjCIJ7gUOnPD8EuG/nmapqDbBmd1eWZKKqxnf3ffZUfe8/+BmAn0Hf+z+XUewa+hZwRJLDk+wNvBG4fAR1SJIYwRZBVT2R5GzgamAJcFFV3TbsOiRJjZFcUFZVVwJXDml1u717aQ/X9/6DnwH4GfS9/7NK1T87TitJ6hHHGpKknlu0QdCnYSyS3JPk1iQbk0y0bQcmWZdka/vzgLY9ST7afi6bkhw92urnJ8lFSbYn2TylbZf7nGRVO//WJKtG0Zf5mKH/5yX5Yfs92JjklCmvvbft/x1JXjOlfY/8PUlyaJLrkmxJcluSd7TtvfkOLKiqWnQPmoPQ3wOeD+wN3AIcOeq6OuzvPcBBO7V9CDi3nT4XOL+dPgX4Cs31HMcC60dd/zz7/ArgaGDzfPsMHAjc1f48oJ0+YNR9243+nwe8e5p5j2x/B54JHN7+bizZk39PgGXA0e30fsB323725juwkI/FukXgMBZNf9e202uB06a0/1U1vgHsn2TZKArcHVV1I/DATs272ufXAOuq6oGqehBYB5zcffW7b4b+z+RU4JKqeryq7gbupPkd2WN/T6pqW1Xd3E4/CmyhGaGgN9+BhbRYg6Bvw1gUcE2SDe0V2QDPq6pt0PzSAM9t2xfzZ7OrfV6Mn8XZ7a6PiyZ3i7DI+59kBXAUsB6/A/OyWINgoGEsFpHjq+pomhFd/yDJK2aZt2+fDczc58X2WXwceAGwEtgGfLhtX7T9T/Js4IvAO6vqkdlmnaZtUXwGC2GxBsFAw1gsFlV1X/tzO3ApzSb//ZO7fNqf29vZF/Nns6t9XlSfRVXdX1VPVtVTwCdovgewSPufZClNCHy6qr7UNvf6OzBfizUIejOMRZJnJdlvcho4CdhM09/JMyBWAZe105cDv9OeRXEs8PDkpvQisKt9vho4KckB7W6Uk9q2PdJOx3pOp/keQNP/NyZ5ZpLDgSOAb7IH/54kCXAhsKWqPjLlpV5/B+Zt1Eeru3rQnCXwXZqzIt4/6no67Ofzac72uAW4bbKvwHOAa4Gt7c8D2/bQ3Bjoe8CtwPio+zDPfn+GZvfHT2j+qjtrPn0Gfpfm4OmdwNtG3a/d7P+n2v5tovmPb9mU+d/f9v8O4NemtO+RvyfAL9PswtkEbGwfp/TpO7CQD68slqSeW6y7hiRJAzIIJKnnDAJJ6jmDQJJ6ziCQpJ4zCLToJVmR5MxBXksynuSjQ6xtLMn6JN9O8iuzzHdCkivmuY6d+/jWJB+bz3tpcTII1AcrgGmDYOfXqmqiqs4ZQk2TXgncXlVHVdXXOlrHCmbuv2QQaDSSfLkdJO+2KQPlkeSsJN9Ncn2ST0z+5dr+5fzFJN9qH8dP854rknwtyc3t4+XtSx8EfqUdo/9dOy32U69N/cu7Hd9/bZJr0tzz4bVJPpTm3g9XtUMckOSlSW5o+3P1dKO5JjksybXtgHDXJlmeZCXNsMmntOvfd6dlTk5ye5KbgNdOaX9WO6jct9otiVPn0f+fb/uwNcmHBvtX06I16ivafPTzwT9d8bkvzVAIzwF+nubeCgcCS4GvAR9r57sY+OV2ejnN0AI7v+fPAPu000cAE+30CcAVM9TxU69NfU4zvv9NbS0vAX5Me1UuzZhOp7WvfR0Ya9vfAFw0zXr+GljVTv8u8OV2+q2Tfdxp/n1oRsU8guaq2M9NqeuPgTe30/vTXBn8rEH7367zLuBn2/V8Hzh01N8JH6N7jOTm9RJwTpLT2+lDaf7j+jnghqp6ACDJ54EXtvO8CjiyGWIGgH+RZL9qxqKftBT4WPuX9pNTlt0dX6mqnyS5leZGLle17bfS7HJ5EfBLwLq2tiU0Qz/s7Dj+6a/6T9FsCczmxcDdVbUVIMn/Aia3nE4CfivJu9vn+9CE430M3v9rq+rh9r2/AxzGTw/HrB4xCDR0SU6g+Y/9uKr6cZLraf4zm25I4EnPaOf/v7PM8y7gfpq/3p8B/MMClPs4QFU9leQnVTU5JstTNL8/AW6rquN28X0HGdtlpnkCvK6q7vipxuQ8Bu//41Omn8T/C3rNYwQahZ8FHmxD4MU0tw6EZkTMf9WOBLkX8Lopy1wDnD35pP2rd7r33VbNMMxvofnrHOBRmtsZTme21wZxBzCW5Li2rqVJfnGa+b5OM7onwJtodjnN5nbg8CQvaJ+fMeW1q4G3tyNwkuSotn0+/ZcMAo3EVcBeSTYBHwC+AVBVP6TZ/70e+FvgO8DD7TLnAOPtwdbvAL8/zfv+D2BVkm/Q7Bb5+7Z9E/BEklumOVg822tzquYWj68Hzk9yC80omC+fZtZzgLe1fX4L8I453vcfaHYF/U17sPj7U17+AM1usE1pbl7/gbZ9Pv2XHH1UTy9Jnl1Vj7VbBJfSHHi9dNR1SYuZWwR6ujkvyUaaM4nuBr484nqkRc8tAknqObcIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeq5/wfrH01+zQWOawAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(true_ages_at_time_of_death, kde=False)\n",
    "plt.xlabel('age at time of death')\n",
    "plt.ylabel('frequency')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Censored data has a different distribution from the actual lifetimes\n",
    "\n",
    "If the year a vampire will die lies in the future (currently meaning after 2018), we could not know without asking a clairvoyant how old the vampire will be when they die. We only know they will be _at least_ the age they are now. We will therefore use their current age, and record the fact that their age at time of death is _censored_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def year_of_death(year_of_birth, age_at_time_of_death):\n",
    "    return year_of_birth + age_at_time_of_death\n",
    "\n",
    "def censored_age(birth_year, age_at_time_of_death):\n",
    "    death_year = year_of_death(birth_year, age_at_time_of_death)\n",
    "    \n",
    "    return (age_at_time_of_death if (death_year <= 2018)\n",
    "            else (2018 - birth_year))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Having a set of data from each vampire: birth year and age at time of death, we know what the\n",
    "# year is in which they die. If this year lies in the future, we know we could not have observed it\n",
    "# without a clairvoyant, so we replace the age with their current age.\n",
    "ages_with_censoring = map(lambda x: censored_age(x[0], x[1]),\n",
    "                    zip(birth_dates_increasing, true_ages_at_time_of_death))\n",
    "\n",
    "# We mark an age as 'observed' if the age with censoring is equal to the true age at time of death\n",
    "# 'observed' is exactly the opposite of 'censored'\n",
    "# (if a data point is observed, it is not censored, and vice versa)\n",
    "observed = map(lambda x: x[0] == x[1], zip(ages_with_censoring, true_ages_at_time_of_death))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHxpJREFUeJzt3XucVWXd9/HPF5wYFboVHL1RREZekqHooAOKeNZAhUctNc+CWnjIPHTjI+pzG5mVlpWaJA+VaYqpaBK3kaEIHsAwMFBwMETJByEhCIUSbg6/54+1Zu4R57DnsPeGWd/367Vfs9a1Tr9rzd77t691uJYiAjMzy652xQ7AzMyKy4nAzCzjnAjMzDLOicDMLOOcCMzMMs6JwMws45wIzMwyzonAzCzjnAjMzDJuh2IHkIvddtstevToUewwzMy2K3PmzPl7RJQ1Nt92kQh69OjB7Nmzix2Gmdl2RdJfc5nPh4bMzDLOicDMLOOcCMzMMm67OEdgZo3buHEjS5cuZf369cUOxQqstLSUbt26UVJS0qzlnQjM2oilS5fSqVMnevTogaRih2MFEhGsWrWKpUuXUl5e3qx15P3QkKT2kv4s6el0vFzSLEmLJD0m6TP5jsEsC9avX0+XLl2cBDJGEl26dGlRS7AQ5wiuAapqjd8B/Dgi9gP+AVxagBjMMsFJIJta+n/PayKQ1A0YAvw8HRdwPPBEOsuDwOn5jMHMzBqW73MEdwH/G+iUjncB1kTEpnR8KbBXnmMwy6RHZr3Xqus777Du9U5bs2YNjzzyCFdeeWWrbrPaAw88wOzZs7n33nsZO3YsO+20ExdddBELFy7knHPOQRJPPPEEPXv2rHP57373u9x0000140cccQQzZ87MS6zbo7wlAklDgRURMUfSsdXFdcwa9Sw/AhgB0L17/W/AxrT2h6EpGvrgmLUla9as4ac//WmdiWDz5s20b9++1bZ1+eWX1wxPnDiR0047jW9961sNLrN1InAS+KR8HhoaCJwqaQnwKMkhobuAXSRVJ6BuwLK6Fo6IcRFRGRGVZWWNdpVhZkU0atQoFi9eTEVFBddffz3Tp0/nuOOO47zzzqNPnz4sWbKEAw88sGb+O++8k9GjRwOwePFiTjrpJA499FCOOuooFi5c2OC2Ro8ezZ133snkyZO56667+PnPf85xxx0HwMMPP0z//v2pqKjgsssuY/PmzYwaNYqPP/6YiooKzj//fAA6duwIwPTp0znmmGP48pe/TK9evRg1ahTjx4+nf//+9OnTh8WLFwOwcuVKzjjjDPr160e/fv2YMWMGAC+88AIVFRVUVFTQt29f1q5d26r7tVDy1iKIiBuBGwHSFsHIiDhf0gTgTJLkMAz4bb5iMLPCuP3225k/fz5z584Fki/YV199lfnz51NeXs6SJUvqXXbEiBGMHTuW/fbbj1mzZnHllVfy/PPPN7rNU045hcsvv5yOHTsycuRIqqqqeOyxx5gxYwYlJSVceeWVjB8/nttvv5177723JratzZs3j6qqKjp37sy+++7LV77yFV599VXuvvtufvKTn3DXXXdxzTXXcN1113HkkUfy3nvvMXjwYKqqqrjzzjsZM2YMAwcOZN26dZSWljZr/xVbMe4juAF4VNJtwJ+BXxQhBjPLs/79+zd6Xfu6deuYOXMmZ511Vk3Zhg0bmrW9qVOnMmfOHPr16wfAxx9/zO67797ocv369aNr164A9OzZk0GDBgHQp08fpk2bBsBzzz3Hm2++WbPMRx99xNq1axk4cCDf+MY3OP/88/nSl75Et27dmhV7sRUkEUTEdGB6OvwO0L8Q2zWz4tl5551rhnfYYQe2bNlSM159zfuWLVvYZZdd6v213hQRwbBhw/je977XpOU6dOhQM9yuXbua8Xbt2rFp06aaOF955RV23HHHTyw7atQohgwZwuTJkzn88MN57rnn2H///VtYk8JzX0Nm1mKdOnVq8Pj4HnvswYoVK1i1ahUbNmzg6aefBuCzn/0s5eXlTJgwAUi+zOfNm9esGE444QSeeOIJVqxYAcDq1av561+TXphLSkrYuHFjs9YLMGjQIO69996a8erEtXjxYvr06cMNN9xAZWVlo+c3tlXuYsKsjSrkVWtdunRh4MCBHHjggZx88skMGTLkE9NLSkq45ZZbOOywwygvL//Er+bx48dzxRVXcNttt7Fx40bOOeccDj744CbH0Lt3b2677TYGDRrEli1bKCkpYcyYMeyzzz6MGDGCgw46iEMOOYTx48c3ed333HMPX/va1zjooIPYtGkTRx99NGPHjuWuu+5i2rRptG/fnt69e3PyySc3ed3bAkXUefXmNqWysjKa+2AaXz5qWVFVVcXnP//5YodhRVLX/1/SnIiobGxZHxoyM8s4JwIzs4xzIjAzyzgnAjOzjHMiMDPLOCcCM7OM830EZm3V7F+27voqL27d9RVRx44dWbduXb3ly5Yt4+qrr+aJJ5JHp5x77rksWLCAiy++mOuuu67OdU6cOJFevXrRu3dvAG655RaOPvpoTjzxxPxVpJU4EZhZm9Ca3V3vueeeNUngb3/7GzNnzqy5S7k+EydOZOjQoTWJ4NZbb22VWArBh4bMrFX86le/4qCDDuLggw/mwgsvBOrvvnn06NFccsklHHvssey7777cc889APzzn/9kyJAhHHzwwRx44IE89thjQNKhXN++fenTpw+XXHJJTcd0PXr04NZbb+XII49kwoQJ9XZp/e677zJgwAD69evHf/7nfzZal9rdZg8aNIgVK1ZQUVHBSy+9VOc2Zs6cyaRJk7j++uupqKhg8eLFDB8+vCaZ9OjRg5tuuokBAwZQWVnJa6+9xuDBg+nZsydjx46t2e4PfvAD+vXrx0EHHcQ3v/nNBvdJa3KLwMxabMGCBXznO99hxowZ7LbbbqxevRqg3u6bARYuXMi0adNYu3Ytn/vc57jiiit45pln2HPPPfnd734HwIcffsj69esZPnw4U6dOpVevXlx00UXcd999XHvttQCUlpby8ssvA0l/Q3V1aX3NNddwxRVXcNFFFzFmzJgm1W3SpEkMHTq0pn+h+rZx6qmnMnToUM4888w617P33nvzyiuvcN111zF8+HBmzJjB+vXrOeCAA7j88suZMmUKixYt4tVXXyUiOPXUU3nxxRdZuXLlp/ZJa3MiMLMWe/755znzzDPZbbfdAOjcuTNQf/fNAEOGDKFDhw506NCB3XffnQ8++IA+ffowcuRIbrjhBoYOHcpRRx3FvHnzKC8vp1evXgAMGzaMMWPG1CSCs88+G2i4S+sZM2bw5JNPAnDhhRdyww03NKueLek2+9RTTwWS7q3XrVtHp06d6NSpE6WlpaxZs4YpU6YwZcoU+vbtW7OtRYsWcdRRR31qn7Q2JwIza7GIQPr0k2jr674ZPtn9c/v27dm0aRO9evVizpw5TJ48mRtvvJFBgwbVfIHWp7q768a6tK4rvqZqSbfZtbu33rrr602bNhER3HjjjVx22WWfWnbrfXLLLbc0vxJ18DkCM2uxE044gccff5xVq1YB1Bwaqq/75vosW7aMnXbaiQsuuICRI0fy2muvsf/++7NkyRLefvttAB566CGOOeaYTy3bUJfWAwcO5NFHHwVoVu+juWyjsa64GzN48GDuv//+mquZ3n//fVasWFHnPmltbhGYtVUFvNzzgAMO4Oabb+aYY46hffv29O3blwceeKDe7pvr88Ybb3D99dfTrl07SkpKuO+++ygtLeWXv/wlZ511Fps2baJfv36feIB9bfV1aX333Xdz3nnncffdd3PGGWe0qK71beOcc87hq1/9Kvfcc0/NSeKmGDRoEFVVVQwYMABILmV9+OGHefvttz+1T1pb3rqhllQKvAh0IEk4T0TENyU9ABwDVJ/xGB4RDf5McDfUZo1zN9TZ1pJuqPPZItgAHB8R6ySVAC9L+n067fqIaHrKNDOzVpe3RBBJU6P61r2S9LXtPwXHzCxj8nqyWFJ7SXOBFcCzETErnfQdSa9L+rGkDg2swsyaYHt44qC1vpb+3/OaCCJic0RUAN2A/pIOBG4E9gf6AZ2BOi/olTRC0mxJs1euXJnPMM3ahNLSUlatWuVkkDERwapVqygtLW32Ogpy1VBErJE0HTgpIu5MizdI+iUwsp5lxgHjIDlZXIg4zbZn3bp1Y+nSpfiHU/aUlpbSrVu3Zi+ft0QgqQzYmCaBHYETgTskdY2I5Uru7jgdmJ+vGMyypKSkhPLy8mKHYduhfLYIugIPSmpPcgjq8Yh4WtLzaZIQMBeo+4JgMzMriHxeNfQ60LeO8uPztU0zM2s6dzFhZpZxTgRmZhnnRGBmlnFOBGZmGedEYGaWcW2+G+qe701odJ7F3c9qdB4zs7bKLQIzs4xzIjAzyzgnAjOzjHMiMDPLOCcCM7OMcyIwM8s4JwIzs4xzIjAzyzgnAjOzjHMiMDPLOCcCM7OMy1sikFQq6VVJ8yQtkPSttLxc0ixJiyQ9Jukz+YrBzMwal88WwQbg+Ig4GKgATpJ0OHAH8OOI2A/4B3BpHmMwM7NG5C0RRGJdOlqSvgI4HngiLX8QOD1fMZiZWePyeo5AUntJc4EVwLPAYmBNRGxKZ1kK7FXPsiMkzZY0e+XKlfkM08ws0/KaCCJic0RUAN2A/sDn65qtnmXHRURlRFSWlZXlM0wzs0wryFVDEbEGmA4cDuwiqfqBON2AZYWIwczM6pbPq4bKJO2SDu8InAhUAdOAM9PZhgG/zVcMZmbWuHw+qrIr8KCk9iQJ5/GIeFrSm8Cjkm4D/gz8Io8xmJlZI/KWCCLidaBvHeXvkJwvMDOzbYDvLDYzyzgnAjOzjHMiMDPLOCcCM7OMcyIwM8s4JwIzs4xzIjAzyzgnAjOzjHMiMDPLOCcCM7OMcyIwM8s4JwIzs4xzIjAzyzgnAjOzjHMiMDPLuHw+mGa70fO9CY3Os7j7WQWIxMys8NwiMDPLuHw+s3hvSdMkVUlaIOmatHy0pPclzU1fp+QrBjMza1yjh4YkdY6I1c1Y9ybgPyLiNUmdgDmSnk2n/Tgi7mzGOs3MrJXl0iKYJWmCpFMkKdcVR8TyiHgtHV4LVAF7NTNOMzPLk1wSQS9gHHAh8Lak70rq1ZSNSOpB8iD7WWnRVZJel3S/pF3rWWaEpNmSZq9cubIpmzMzsyZoNBFE4tmIOBf4CjAMeFXSC5IGNLa8pI7Ak8C1EfERcB/QE6gAlgM/rGe74yKiMiIqy8rKcq+RmZk1SS7nCLoAF5C0CD4Avg5MIvkinwCUN7BsCUkSGB8RvwGIiA9qTf8Z8HQL4jczsxbK5T6CV4CHgNMjYmmt8tmSxta3UHo+4RdAVUT8qFZ514hYno5+EZjf9LDNzKy15JIIPhcRUdeEiLijgeUGkrQi3pA0Ny27CThXUgUQwBLgstzDNTOz1pZLIpgi6ayIWAOQntx9NCIGN7RQRLwM1HWV0eSmh2lmZvmSy1VDZdVJACAi/gHsnr+QzMyskHJJBJslda8ekbQPyWEdMzNrA3I5NHQz8LKkF9Lxo4ER+QvJzMwKqdFEEBHPSDoEOJzkmP91EfH3vEdmZmYFkWs31B2A1en8vSURES/mLywzMyuUXG4ouwM4G1gAbEmLA3AiMDNrA3JpEZxOci/BhnwHY2ZmhZfLVUPvACX5DsTMzIojlxbBv4C5kqYCNa2CiLg6b1GZmVnB5JIIJqUvMzNrg3K5fPRBSTsC3SPirQLEZGZmBdToOQJJ/wuYCzyTjldIcgvBzKyNyOVk8WigP7AGICLm0sAzCMzMbPuSSyLYFBEfblXmvobMzNqIXE4Wz5d0HtBe0n7A1cDM/IZlZmaFkkuL4OvAASSXjv4a+Ai4Np9BmZlZ4eRy1dC/SHogvTn/4ZiZWaHl0tfQNOo4JxARxzey3N7Ar4B/J+mjaFxE3C2pM/AY0IPkUZVfTh92Y2ZmRZDLOYKRtYZLgTOATTkstwn4j4h4TVInYI6kZ4HhwNSIuF3SKGAUcEPTwjYzs9aSy6GhOVsVzaj1kJqGllsOLE+H10qqAvYCTgOOTWd7EJiOE4GZWdHkcmioc63RdsChJId7ciapB9AXmAXskSYJImK5pDqffyxpBOmT0Lp3717XLGZm1gpyOTQ0h+QcgUgO97wLXJrrBiR1BJ4Ero2IjyTltFxEjAPGAVRWVvq+BTOzPMnl0FCz7yKWVEKSBMZHxG/S4g8kdU1bA12BFc1dv5mZtVwuh4a+1ND0Wl/wWy8n4BdAVUT8qNakScAw4Pb0729zjtbMzFpdLoeGLgWOAJ5Px48jOcH7IckhozoTATAQuBB4Q9LctOwmkgTwuKRLgfeAs5oVuZmZtYpcEkEAvatP8KaHc8ZExMUNLhTxMsl5hbqc0KQozcwsb3LpYqJHdRJIfQD0ylM8ZmZWYLm0CKZL+gNJP0MBnANMy2tUZmZWMLlcNXSVpC8CR6dF4yLiqfyGZWZmhZJLiwDgNWBtRDwnaSdJnSJibT4DMzOzwsjlUZVfBZ4A/m9atBcwMZ9BmZlZ4eRysvhrJJeCfgQQEYuAOruFMDOz7U8uiWBDRPx39YikHfCjKs3M2oxcEsELkm4CdpT0BWAC8F/5DcvMzAoll0QwClgJvAFcBkwG/k8+gzIzs8Jp8KohSe2BByPiAuBnhQnJzMwKqcFEEBGbJZVJ+kzt8wSWm0dmvVfsEArqvMP83Aiz7VEu9xEsIXkq2STgn9WFW/UoamZm26l6zxFIeigdPBt4Op23U62XmZm1AQ21CA6VtA9JV9E/KVA8ZmZWYA0lgrHAM0A5MLtWuUjuI9g3j3GZmVmB1HtoKCLuiYjPA7+MiH1rvcojwknAzKyNaPQ+goi4ojkrlnS/pBWS5tcqGy3pfUlz09cpzVm3mZm1nlxuKGuuB4CT6ij/cURUpK/Jedy+mZnlIG+JICJeBFbna/1mZtY68tkiqM9Vkl5PDx3tWoTtm5lZLYVOBPcBPYEKYDnww/pmlDRC0mxJs1euXFmo+MzMMqegiSAiPoiIzRGxhaTvov4NzDsuIiojorKsrKxwQZqZZUxBE4GkrrVGvwjMr29eMzMrjFyfWdxkkn4NHAvsJmkp8E3gWEkVJDekLSHp1trMzIoob4kgIs6to/gX+dqemZk1TzGuGjIzs22IE4GZWcY5EZiZZZwTgZlZxjkRmJllXN6uGmprer43ocHpi7ufVaBIzMxal1sEZmYZ50RgZpZxTgRmZhnnRGBmlnFOBGZmGedEYGaWcU4EZmYZ50RgZpZxTgRmZhnnRGBmlnFOBGZmGZe3RCDpfkkrJM2vVdZZ0rOSFqV/d83X9s3MLDf5bBE8AJy0VdkoYGpE7AdMTcfNzKyI8pYIIuJFYPVWxacBD6bDDwKn52v7ZmaWm0KfI9gjIpYDpH93L/D2zcxsK9vs8wgkjQBGAHTv3r3I0TSusecVgJ9ZYGbbpkK3CD6Q1BUg/buivhkjYlxEVEZEZVlZWcECNDPLmkIngknAsHR4GPDbAm/fzMy2ks/LR38NvAJ8TtJSSZcCtwNfkLQI+EI6bmZmRZS3cwQRcW49k07I1zbNzKzpttmTxWZN8cis94qy3fMO2/YvZDBrjLuYMDPLOCcCM7OMcyIwM8s4JwIzs4zzyeIC8t3HZrYtcovAzCzjnAjMzDLOicDMLOOcCMzMMs6JwMws43zV0DamsSuLfFWRmbU2twjMzDLOicDMLOOcCMzMMs6JwMws43yy2FpNsZ4JYGYtU5REIGkJsBbYDGyKiMpixGFmZsVtERwXEX8v4vbNzAyfIzAzy7xiJYIApkiaI2lEkWIwMzOKd2hoYEQsk7Q78KykhRHxYu0Z0gQxAqB7dz8g3MwsX4rSIoiIZenfFcBTQP865hkXEZURUVlWVlboEM3MMqPgiUDSzpI6VQ8Dg4D5hY7DzMwSxTg0tAfwlKTq7T8SEc8UIQ4zM6MIiSAi3gEOLvR2zcysbr6z2KwFink39XmH+SIKax2+j8DMLOOcCMzMMs6JwMws45wIzMwyzonAzCzjnAjMzDLOicDMLOOcCMzMMs6JwMws43xncRvU870Jjc6zuPtZBYjE2iI/m7qwCnEHuVsEZmYZ50RgZpZxTgRmZhnnRGBmlnFOBGZmGeerhrYzuVwRZNngq3estRSlRSDpJElvSXpb0qhixGBmZoliPLy+PTAGOBnoDZwrqXeh4zAzs0QxWgT9gbcj4p2I+G/gUeC0IsRhZmYUJxHsBfy/WuNL0zIzMyuCYpwsVh1l8amZpBHAiHR0naS3mrm93YC/N3PZtqCe+o8seCBFlPX3AHgfbLf1P79li++Ty0zFSARLgb1rjXcDlm09U0SMA8a1dGOSZkdEZUvXs73Kev3B+wC8D7Je/8YU49DQn4D9JJVL+gxwDjCpCHGYmRlFaBFExCZJVwF/ANoD90fEgkLHYWZmiaLcUBYRk4HJBdpciw8vbeeyXn/wPgDvg6zXv0GK+NR5WjMzyxD3NWRmlnFtNhFkqRsLSUskvSFprqTZaVlnSc9KWpT+3TUtl6R70v3yuqRDiht980i6X9IKSfNrlTW5zpKGpfMvkjSsGHVpjnrqP1rS++n7YK6kU2pNuzGt/1uSBtcq3y4/J5L2ljRNUpWkBZKuScsz8x5oVRHR5l4kJ6EXA/sCnwHmAb2LHVce67sE2G2rsu8Do9LhUcAd6fApwO9J7uc4HJhV7PibWeejgUOA+c2tM9AZeCf9u2s6vGux69aC+o8GRtYxb+/0M9ABKE8/G+23588J0BU4JB3uBPwlrWdm3gOt+WqrLQJ3Y5HU98F0+EHg9Frlv4rEH4FdJHUtRoAtEREvAqu3Km5qnQcDz0bE6oj4B/AscFL+o2+5eupfn9OARyNiQ0S8C7xN8hnZbj8nEbE8Il5Lh9cCVSQ9FGTmPdCa2moiyFo3FgFMkTQnvSMbYI+IWA7JhwbYPS1vy/umqXVui/viqvTQx/3Vh0Vo4/WX1APoC8zC74FmaauJIKduLNqQgRFxCEmPrl+TdHQD82Zt30D9dW5r++I+oCdQASwHfpiWt9n6S+oIPAlcGxEfNTRrHWVtYh+0hraaCHLqxqKtiIhl6d8VwFMkTf4Pqg/5pH9XpLO35X3T1Dq3qX0RER9ExOaI2AL8jOR9AG20/pJKSJLA+Ij4TVqc6fdAc7XVRJCZbiwk7SypU/UwMAiYT1Lf6isghgG/TYcnARelV1EcDnxY3ZRuA5pa5z8AgyTtmh5GGZSWbZe2OtfzRZL3AST1P0dSB0nlwH7Aq2zHnxNJAn4BVEXEj2pNyvR7oNmKfbY6Xy+SqwT+QnJVxM3FjieP9dyX5GqPecCC6roCXYCpwKL0b+e0XCQPBloMvAFUFrsOzaz3r0kOf2wk+VV3aXPqDFxCcvL0beDiYterhfV/KK3f6yRffF1rzX9zWv+3gJNrlW+XnxPgSJJDOK8Dc9PXKVl6D7Tmy3cWm5llXFs9NGRmZjlyIjAzyzgnAjOzjHMiMDPLOCcCM7OMcyKwNk9SD0nn5TJNUqWkewoYW5mkWZL+LOmoBuY7VtLTzdzG1nUcLune5qzL2iYnAsuCHkCdiWDraRExOyKuLkBM1U4AFkZE34h4KU/b6EH99TdzIrDikDQx7SRvQa2O8pB0qaS/SJou6WfVv1zTX85PSvpT+hpYxzp7SHpJ0mvp64h00u3AUWkf/ddttdgnptX+5Z327/+gpClKnvnwJUnfV/Lsh2fSLg6QdKikF9L6/KGu3lwl7SNpatoh3FRJ3SVVkHSbfEq6/R23WuYkSQslvQx8qVb5zmmncn9KWxKnNaP+e6Z1WCTp+7n916zNKvYdbX5l88X/3PG5I0lXCF2APUmerdAZKAFeAu5N53sEODId7k7StcDW69wJKE2H9wNmp8PHAk/XE8cnptUeJ+nf/+U0loOBf5HelUvSp9Pp6bSZQFlafjZwfx3b+S9gWDp8CTAxHR5eXcet5i8l6RVzP5K7Yh+vFdd3gQvS4V1I7gzeOdf6p9t8B/i3dDt/BfYu9nvCr+K9ivLwejPgaklfTIf3Jvni+nfghYhYDSBpAtArnedEoHfSxQwAn5XUKZK+6KuVAPemv7Q311q2JX4fERslvUHyIJdn0vI3SA65fA44EHg2ja09SdcPWxvA//yqf4ikJdCQ/YF3I2IRgKSHgeqW0yDgVEkj0/FSkuS4jNzrPzUiPkzX/SawD5/sjtkyxInACk7SsSRf7AMi4l+SppN8mdXVJXC1dun8Hzcwz3XAByS/3tsB61sh3A0AEbFF0saIqO6TZQvJ50fAgogY0MT15tK3S33zCDgjIt76RKE0mtzrv6HW8Gb8XZBpPkdgxfBvwD/SJLA/yaMDIekR85i0J8gdgDNqLTMFuKp6JP3VW9d6l0fSDfOFJL/OAdaSPM6wLg1Ny8VbQJmkAWlcJZIOqGO+mSS9ewKcT3LIqSELgXJJPdPxc2tN+wPw9bQHTiT1TcubU38zJwIrimeAHSS9Dnwb+CNARLxPcvx7FvAc8CbwYbrM1UBlerL1TeDyOtb7U2CYpD+SHBb5Z1r+OrBJ0rw6ThY3NK1RkTzi8UzgDknzSHrBPKKOWa8GLk7rfCFwTSPrXU9yKOh36cniv9aa/G2Sw2CvK3l4/bfT8ubU38y9j9q2RVLHiFiXtgieIjnx+lSx4zJry9wisG3NaElzSa4keheYWOR4zNo8twjMzDLOLQIzs4xzIjAzyzgnAjOzjHMiMDPLOCcCM7OMcyIwM8u4/w/1iSy8LERoRAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(true_ages_at_time_of_death, kde=False, label='true lifetimes')\n",
    "sns.distplot(ages_with_censoring, kde=False, label='censored lifetimes')\n",
    "plt.xlabel('age at time of death')\n",
    "plt.ylabel('frequency')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the distribution of censored lifetimes differs significantly from the true distribution (which one would only know if they could look into the future). We can compare the two distributions by drawing survival curves.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "km = lifelines.KaplanMeierFitter()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,0,'age (years)')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEKCAYAAAAb7IIBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X98VOWV+PHPIaABDYgQUgsRCF/QQBIBEbAoJEUQEbFaf9H2VdG6rFpbratb2t0ipd2WpVRR0Vpr/VHrCv5oKcvSalXQiopiQQpEKaBohEoIFBMhQuB8/7hzx2EyP+5M7mR+nffrlddkZu7c+9xJMifnPs9zHlFVjDHG5J8O6W6AMcaY9LAAYIwxecoCgDHG5CkLAMYYk6csABhjTJ6yAGCMMXnKAoAxxuQpCwDGGJOnLAAYY0ye6piuA/fs2VP79euXrsMbY0xWevPNN3erarEf+0pbAOjXrx9r1qxJ1+GNMSYrich2v/Zll4CMMSZPWQAwxpg8ZQHAGGPyVNr6AEzuOHToEHV1dTQ3N6e7KcbkjMLCQvr06UOnTp1SdgwLAKbN6urqKCoqol+/fohIuptjTNZTVRoaGqirq6N///4pO07cS0Ai8qCI7BKRDVGeFxG5S0S2iMh6ERnufzNNJmtubqZHjx724W+MT0SEHj16pDyr9tIH8DAwKcbz5wEDA18zgF+0vVkm29iHvzH+ao+/qbgBQFVfAvbE2ORC4DfqeA04QURO8quB4V6791/Y+JOzWP3kz1N1CGOMyQt+jALqDXwQcr8u8FgrIjJDRNaIyJr6+vqkDzjk4N84/u+/T/r1Jnf885//5N57703Z/h9++GFuuOEGAO677z5+85vfAPD2228zdOhQhg0bxtatW6O+/ic/+clR97/whS+krK3GJMqPABApT4m40ryq3q+qI1R1RHFxcjOZR1//KzYeU8mRlpakXm9yS6wAcPjwYV+Pde211/L1r38dgCVLlnDhhReydu1aBgwYEPU14QHglVde8bVNxrSFHwGgDigNud8H2OHDfmNSYMU7u1J9GJPhZs6cydatWxk6dCi33norK1eupKamhq985StUVlby3nvvUVFREdx+/vz5zJ49G4CtW7cyadIkTj/9dM4++2zefvvtmMeaPXs28+fPZ/ny5SxYsIAHHniAmpoaAH77298ycuRIhg4dyr/+679y+PBhZs6cyYEDBxg6dChf/epXATj++OMBWLlyJePGjeOyyy5j0KBBzJw5k8cee4yRI0dSWVkZzCrq6+v58pe/zBlnnMEZZ5zBqlWrAHjxxRcZOnRoMAtpbGz09X01+cGPYaBLgRtEZBEwCtinqjt92G9sqjQeOJTyw5jE/PB/N7Jpx8e+7nPw57ty2wVDIj43d+5cNmzYwLp16wDng/X1119nw4YN9O/fn/feey/qfmfMmMF9993HwIEDWb16Nddffz0vvPBC3PZMnjyZa6+9luOPP55bbrmF2tpaFi9ezKpVq+jUqRPXX389jz32GHPnzmXhwoXBtoV76623qK2t5cQTT6SsrIxrrrmG119/nTvvvJO7776bBQsWcOONN/Kd73yHs846i/fff59zzz2X2tpa5s+fzz333MOYMWNoamqisLAw/htpTJi4AUBEHgeqgZ4iUgfcBnQCUNX7gOXAZGALsB+4KlWNDdfccpgV7+yi5pRe7XVIkwVGjhwZd+x0U1MTr7zyCpdeemnwsU8//TSp4z3//PO8+eabnHHGGQAcOHCAXr3i/06eccYZnHSSM15iwIABTJw4EYDKykpWrFgBwHPPPcemTZuCr/n4449pbGxkzJgx3HzzzXz1q1/l4osvpk+fPkm13eS3uAFAVafFeV6Bb/rWIo/2HxY2/6OJwo4F7X1oE0O0/9Tb03HHHRf8vmPHjhw5ciR43x1XfeTIEU444YSo/50nQlW58sor+elPf5rQ64499tjg9x06dAje79ChAy2BPq4jR47w6quv0rlz56NeO3PmTM4//3yWL1/O6NGjee655zj11FPbeCYm32RlLaCexzt/KKu27g5mASY/FRUVxbz+XVJSwq5du2hoaODTTz9l2bJlAHTt2pX+/fvz5JNPAs6H+FtvvZVUG8aPH89TTz3Frl3O7+GePXvYvt2p2NupUycOHUr+UuXEiRNZuHBh8L4bsLZu3UplZSXf/e53GTFiRNz+C2MiycoAUFJUSJcCZ6BRaffjqG9stiCQp3r06MGYMWOoqKjg1ltvbfV8p06dmDVrFqNGjWLKlClH/Zf82GOP8etf/5rTTjuNIUOG8Ic//CGpNgwePJgf//jHTJw4kaqqKiZMmMDOnU432IwZM6iqqgp2AifqrrvuYs2aNVRVVTF48GDuu+8+ABYsWEBFRQWnnXYanTt35rzzzktq/ya/iXMFp/2NGDFCk14Q5qHzWf/+br57ws+YNcW55PDB3k+Cl4OKOneyfoF2VFtbS3l5ebqbYUzOifS3JSJvquoIP/aflRmAq3ZnI8/XfgQ4mUBxUSHFRYU2OsgYYzzI2gDQ7RinY2/V1t1pbokxxmSnrA0AJ3ZqoeqEg0dlAcYYY7zL2gAAUPM5Z9y2ZQHGGJO4rA4A5/dppvykolZZgA0NNcaY+LI6AACMGdATODoLKO1+nHUEG2NMHFkfAMaXl1B+UhHbG/YzZ9lG5izbaH0CplURuExRXV1NpOHPoY9PnjyZf/7zn4AzD6C8vDzmPIJ169axfPny4P2lS5cyd+5cn1se3ZIlS5gzZ067HS/bhZYYX7hwIQ899FDa2pITawI7WYCTAWxv2A/spqpPt7S2yeSelpYWOnZM/Z9M6If5vffeyx//+MeYtY3WrVvHmjVrmDx5MgBTp05l6tSpKW+na968eSxdurTdjpeJDh8+TEFB4mVprr76asaMGcNVV7VbCbWjZG0GUNS0jW7/cGqrjy8vYdaUIcyaMoS+PboA1g+QT26//XYqKiqoqKhgwYIFwcdbWlq48sorqaqq4pJLLmH//v2AU0dn8ODBVFVVccsttwDRyy7Pnj2bGTNmMHHiRL7+9a8zatQoNm7cGDxGdXU1b775Jp988glXX301Z5xxBsOGDQvOKj5w4ABXXHEFVVVVXH755Rw4cCDu+fTr14/du3dz7bXXsm3bNqZOncodd9wR8RgHDx5k1qxZLF68mKFDh7J48eKj/sOcPn061113HTU1NZSVlfHiiy9y9dVXU15ezvTp04PHfPbZZznzzDMZPnw4l156KU1NTVHfq1CbN2/m2GOPpWdP51Lsk08+GZyhPHbsWMCpv3TVVVdRWVnJsGHDgoXuHn74Yb70pS9xwQUX0L9/fxYuXMjtt9/OsGHDGD16NHv2OAsReinb3dTUFDxGVVUVTz/9dMzz6tevH7fddhvDhw+nsrIyuM9IZbZVlVtvvZWKigoqKytZvHgxQKvS4xC5LDjAQw89xKBBgxg3blzwdwugS5cu9OvXj9dffz3u70VKqGpavk4//XRN2hsPqt7WVf/5s6H6ypbdR31NWvCiTlrwor6yZbf+YW1d8scwnm3atOmzO8u/q/rgZH+/ln836rHXrFmjFRUV2tTUpI2NjTp48GD961//qu+++64C+vLLL6uq6lVXXaU/+9nPtKGhQQcNGqRHjhxRVdW9e/eqquq0adP0L3/5i6qqbt++XU899VRVVb3tttt0+PDhun//flVVvf3223XWrFmqqrpjxw4dOHCgqqp+73vf00cffTS4z4EDB2pTU5P+/Oc/16uuukpVVd966y0tKCjQN954o9V5jBs3Lvh43759tb6+vtX30Y7x0EMP6Te/+c3gvkLvX3nllXr55ZfrkSNHdMmSJVpUVKTr16/Xw4cP6/Dhw3Xt2rVaX1+vZ599tjY1Namq6ty5c/WHP/xh1Pcq1IMPPqg333xz8H5FRYXW1dUdtf38+fN1+vTpqqpaW1urpaWleuDAAX3ooYd0wIAB+vHHH+uuXbu0a9eu+otf/EJVVW+66Sa94447VFX1i1/8om7evFlVVV977TWtqalp1Y5///d/1xtvvDF4f8+ePVHPy31f77rrLlVVveeee/Qb3/iGqqpOmTIl+DvT2Niohw4d0qeeekrPOeccbWlp0X/84x9aWlqqO3bs0BUrVmiXLl1027Ztqur8HUyZMkUPHjyoqqrXXXedPvLII7pjxw4tLS3VXbt26aeffqpf+MIXjvp5/fjHP9b58+e3Oid3n+GANerT53B2XgIacRW89gs6HPiUbh++yL7e49LdIpMmL7/8MhdddFGwAujFF1/MX/7yF6ZOnUppaSljxowB4Gtf+xp33XUXN910E4WFhVxzzTWcf/75TJkyBYhedhmcSypuNc7LLruMCRMm8MMf/pAnnngiWE762WefZenSpcyfPx9w/ut9//33eemll/j2t78NQFVVFVVVVUmfa7RjxHPBBRcgIlRWVlJSUhL8b3XIkCG899571NXVsWnTpuB7dfDgQc4880y6du0a8b0KtXPnTkJX9xszZgzTp0/nsssu4+KLLwacn9G3vvUtAE499VT69u3L5s2bAaipqaGoqIiioiK6devGBRdcADglsdevX++5bPdzzz3HokWLgve7d+/OsmXLIp6Xy23f6aefzu9+97tg+8PLbL/88stMmzaNgoICSkpKGDduHG+88QZdu3Y9qvR4tLLgq1evprq6Ovg+XX755cHzB+jVq1faivllZwAIOFJwLB0P+rv4iGmj89qv8xGcDDYaEWl1v2PHjrz++us8//zzLFq0iIULF/LCCy9ELbsMR5eX7t27Nz169GD9+vUsXryYX/7yl8F2PP3005xyyilx25GsaMdYvXp1zNeFlpkOL0Hd0tJCQUEBEyZM4PHHH2/12kjvVajOnTuzb9++4P377ruP1atX83//938MHTqUdevWxfwZxSuJ7bVst6q2ep9VNep5hR67oKAgWH47UpntWO0P/d3QKGXBlyxZEvN3oLm5OeLvXXvI2j6AWGx2cP4YO3YsS5YsYf/+/XzyySf8/ve/5+yzzwbg/fff59VXXwXg8ccf56yzzqKpqYl9+/YxefJkFixYEPxgiVZ2OZIrrriCefPmsW/fvuB/0+eeey5333138MNi7dq1wfY99thjAGzYsIH169cnfa7RjhGvJHY8o0ePZtWqVWzZsgWA/fv3s3nz5qjvVajy8vLg68C5Xj9q1CjmzJlDz549+eCDD456DzZv3sz7778fMVBG4rVsd/jPb+/evVHPK5ZIZbbHjh3L4sWLOXz4MPX19bz00kuMHDmy1WujlQUfNWoUK1eupKGhgUOHDgXPxbV58+a0jVjLuQAQOi/AOoJz3/Dhw5k+fTojR45k1KhRXHPNNQwbNgxwPpweeeQRqqqq2LNnD9dddx2NjY1MmTKFqqoqxo0bxx133AFEL7scySWXXMKiRYu47LLLgo/94Ac/4NChQ1RVVVFRUcEPfvADAK677jqampqoqqpi3rx5ET84vIp2jJqaGjZt2hTsBE5UcXExDz/8MNOmTaOqqorRo0fz9ttvR32vQo0dO5a1a9cGg9Ktt95KZWUlFRUVjB07ltNOO43rr7+ew4cPU1lZyeWXX87DDz981H/+8Xgp2/2f//mf7N27N9gBvWLFiqjnFUukMtsXXXQRVVVVnHbaaXzxi19k3rx5fO5zn2v12mhlwU866SRmz57NmWeeyTnnnMPw4cOPet2qVas455xzPL8ffsrOctAAC0fS2CJs6H8VVFx61FNzljmjNGZNGUJ9YzNTh/ZuS1NNHFYOOr/deOONXHDBBWn7EMtma9eu5fbbb+fRRx+N+LyVg46h6NiOdOpgS0Iak07f//73g0NsTWJ2797Nj370o7QdP6sDAEDB4Wa6ffhiupthTN4qKSlp14lnuWTChAn069cvbcfP+gDAif3Yt9eu86dbui4lGpOr2uNvKusDwLDSEyjUQ62yAHckUHPLYZau+9A6g1OosLCQhoYGCwLG+ERVaWhooLCwMKXHyep5AEEn9mPf3g8h0Nc7ZkBPanc2smrrbsaXlwBQ39icxgbmtj59+lBXV0d9fX26m2JMzigsLKRPnz4pPUZ2B4CPNsDmPzFs0CTW7PsHbgHo8eUltkhMO+rUqVPMYmXGmMyUvZeASkc5t9tWetrc5gQYY8zRsjcDKBsHO/4avNv5mAJ2BC7zNLccbrV5affj7DKQMcaEyN4MIMyQ4mOYetwmpg7tTWFHmxtgjDHxZG8AKOwGLSFVAbv3hcaP4O9/Dj5kawUbY0x02RsABk4ACWt+977Q7FQmtLWCjTEmtuwNAHG4awWHs3kBxhjj8BQARGSSiLwjIltEZGaE508WkRUislZE1ovIZP+b6o/S7sdRXFRomYAxJu/FDQAiUgDcA5wHDAamicjgsM3+E3hCVYcBVwD3+t3QqAJzARJlmYAxJt95yQBGAltUdZuqHgQWAReGbaNA18D33YAd/jUxhihzAYo6d+KDvZ/EfqllAsaYPOclAPQGPgi5X0ew6ELQbOBrIlIHLAe+5Uvr4ikbByUVR2cBh5qp6fBWcCjo9ob9zFm20VYIM8aYMF4CQKTFLMOrfk0DHlbVPsBk4FGR8CE6ICIzRGSNiKzxrW5MWbVz62YBYSOB+vbowvaG/VFLQ9jQUGNMvvISAOqA0pD7fWh9iecbwBMAqvoqUAj0DN+Rqt6vqiNUdURxcXFyLQ43aJKTBUQwvryEWVOG0LdHl6gvt6Ghxph85SUAvAEMFJH+InIMTifv0rBt3gfGA4hIOU4AsNKQxhiTweIGAFVtAW4AngFqcUb7bBSROSLiLgP0b8C/iMhbwOPAdM2w4vDhs4KNMSbfeZoHoKrLVXWQqg5Q1f8KPDZLVZcGvt+kqmNU9TRVHaqqz6ay0XEdaqb4Hy8F70aaFRzK+gGMMfkod2YCh44E6t6XrnIgOBQ02qxgl/UDGGPyUW4EgLJq5zZkPsCQz3dNqCqoZQHGmHyT3QGgsBvs3f7ZSKA4s4Jj9QO46wUsXfehzRA2xuSF7A4AAydAp8CiyWXVzm2UFcLi9QPAZ7ODbYawMSYfZHcACBVjPgB81g/gdWaw1QoyxuS67F0SMglOFrCb7Q37gd2MLy+Jum1p9+MAbBlJY0zOyp0MINyhZvrufeWoonChM4OtRpAxJt/lbgDo3pfTekrEkUBeagQZY0yuy70A4GF9AC81glw2PNQYk6tyKwCUVTu3UUYCJcMmiRljclVuBYA4I4Ei8VIjyLIAY0wuyv4A4E4GS4KXuQFgWYAxJjdlfwAInQyWoHg1gowxJpdlfwCAmFlAUedO1Dc2x10jOB67DGSMyTW5EQDCs4CQkUA1p/Ri6tDeMQvDeekHsMtAxphckxsBIFRZtXO7bSUcaoa//znm5l77AYwxJtfkXgAIHQkUskB8NG4/gK0YZozJN7kXAKIo6twpaj+AmwU88PK7zFm2MWqJCOsHMMbkkrwJADWn9IraDzC+vIRrzuofHBEUrUSE9QMYY3JJ3gSAeNzyEPGKxVkWYIzJFbkTANowISxcrGJxlgUYY3JF7gSAGENBE5VIsThjjMlWuRMA4LMsoKzauR82FDRWR7AxxuSb3AoAbhYQZShorI7gaGx4qDEmV+VWAPAgkdIQ0SaJuesF25rBxphslncBwEtpCFe0SWKl3Y+juKiQ4qJC6xA2xmStvAsALq+ZgJWKMMbkqrwNAF4zASsZbYzJVbkXAHycD2CMMbks9wJAGxaIicVGAxljck3HdDcg3dy5AaXdj4u6zZgBPand2ciqrbsZX15y1HPuiKDwfdac0isl7TXGGL94ygBEZJKIvCMiW0RkZpRtLhORTSKyUUT+x99mJsnDbGAvcwNi9QOEjgiykUHGmGwSNwCISAFwD3AeMBiYJiKDw7YZCHwPGKOqQ4CbUtDWxJRVO7ceF4bxkxWMM8ZkAy8ZwEhgi6puU9WDwCLgwrBt/gW4R1X3Aqhq+j/9ElwYxguv/QBWMM4Ykw289AH0Bj4IuV8HjArbZhCAiKwCCoDZqtrq2ouIzABmAJx88snJtNebBEcCtbUfIBLrGzDGZDovGYBEeEzD7ncEBgLVwDTgARE5odWLVO9X1RGqOqK4uDjRtnqX4EigRPoBEskCrG/AGJPJvASAOqA05H4fYEeEbf6gqodU9V3gHZyAkDni9AN4qRRqs4KNMbnESwB4AxgoIv1F5BjgCmBp2DZLgBoAEemJc0lom58NbbM4/QBtHQ1kjDHZJm4AUNUW4AbgGaAWeEJVN4rIHBGZGtjsGaBBRDYBK4BbVbUhVY1OFbc+kNdqoYmy0UHGmEziaSKYqi4Hloc9NivkewVuDnxlFncuwKBJcTcN7aAN78D1Q2n346hvbPZ9v8YYk4zcnglcVu0EgG0rnQDg9gMMnJC2JkUaHQQ2QsgY0/5yOwAMmuR8+Lu694Umb/V83MtBzS2HWw0P3d6wnznLNgJOx7CXYaGuaENNP9j7CSve2WVBwBjTbnI7ALSB+0Ec/t+6MxLIGQW0vWE/4G1eQDx2ecgY095yrxqoK9pksDaWhRhfXsKsKUOYNWUIfXt0aUMDW7NOYmNMe8rdABBtMphPZSFSwUpIGGPaU+4GAJ94mSBmjDHZyAJAHF4miBljTDayTmAPYhWLCx0R5Ep0ZFAoKyJnjGkvFgA8qDmlV8Sx+6EjglxtHRkUKcjY6CBjTCrkRwBIYDZwIsaXl7T6oJ+zbGOwYqgfw0Phs9FBlgUYY/yU+30AZdXObeiEsBRKRcVQGx1kjEmF3M8AwmcDQ1IlIbwsGgNOVrBq6+6U9A1YFmCM8VPuZwCRJDEXIJHRQGMG9Gw1SWx7w/42ZQWWBRhj/JbbGUDobODwfoAUFoaL1jfQVpYFGGP8lNsZgDsbuKzauR9eGC5DZwRHY1mAMcZPuR0AXIMmQUlF68cTrAvU1lnBXtcTjsXqBRlj/JIfASCaBLOAtswK9mt0kGUBxhi/5HYfQLhI8wES7AvwOhooXKzRQZDYCCHrCzDG+CF/MoCyauc2fEhoO2cBkUpIJzpCyLIAY4wf8icDiDQfoJ1FGh0EqZk9bIwx8eRPBpDBUjF72Bhj4sn9ABC+MpjbD5BBxpeXUH5SUUKjhGw0kDGmrXI/AISuDFZW7dym+VJQJIlmAdYPYIxpq9wPAKFizQf421Oe5wQUde5EfWMz9Y3Nvq0W5mYBiXDXDrBMwBiTjPzpBI6le1/ntsnb5ZfQ4ZeR1gloL+5QVFsvwBiTDAsAGcbvKqLGGBONBYAMkooVxowxJhoLAKFSWCHUi2RXGAtdR9jWDzbGeJVfncDxJFEhNLRD2O+OYfA2Oqi0+3EUFxVSXFRoI4OMMZ7lZwbg4xrBkf7b9rNj2K0h5HWmsGUDxhivPGUAIjJJRN4RkS0iMjPGdpeIiIrICP+a6IPQyWBl1c5tCucCtLVsdLhE5ghYNmCM8SpuABCRAuAe4DxgMDBNRAZH2K4I+Daw2u9GtlnoZDB3LkC0GcEJrhEQSVsKxkWSzExh+CwbcL9svoAxJpSXDGAksEVVt6nqQWARcGGE7X4EzAMyf1B6WbVzGykL8GmlML/7BpKpFxSaDVhGYIwJ56UPoDfwQcj9OmBU6AYiMgwoVdVlInKLj+1LDbcyqI99AeH87htItC/AGGPi8ZIBSITHNPikSAfgDuDf4u5IZIaIrBGRNfX19d5b6YfwonBl1c5tBtYFiqatVUNDLwnZ5SBjjJcAUAeUhtzvA+wIuV8EVAArReQ9YDSwNFJHsKrer6ojVHVEcXFx8q1ORmg/AMTvC8hAyfYFuKyD2BgTyksAeAMYKCL9ReQY4Apgqfukqu5T1Z6q2k9V+wGvAVNVdU1KWuynsmrnNjwL8KEjOBI/Rgf5tXaAZQPGmLgBQFVbgBuAZ4Ba4AlV3Sgic0RkaqobmFLRqoP61BEczo/RQclUDY3EsgFjjKeJYKq6HFge9tisKNtWt71ZKeL2A7jVP9PAHR0UrrnlcMILzfvFzQZs4pgx+SW/ZgIPnODU/U+jaB+wiY4Q8nM0kJWVNiY/WS2gaBJcJKatEukfSNUawrbMpDH5xQIAwJ534U8zjx4N1L0vHF+Skr6ASBLpH3D7Ady1A5IZERSJLTNpTH7JvwAQaT7Aif2dIBBpTkCKRgRFkmgW0LdHF7Y37Pc1E7BlJo3JH/kXACLNB5g01wkCkaRoRFAkiWYBs6YMoW+PLr62wR0dZJmAMbkv/wKA8cQyAWNyX36NAoonhbWBvIo0TDTeENHQdYT9Wj/YRgYZk/ssALjKqp0AsG1l6wDQjktFJlpELnQd4VSsHxy6wAzYIjPG5BILAC63Qmgk3ftCkz8jbZLhdg5HygJC1xGes2yj79lA+DEtIzAmd1gAyAI1p/TyNFEs1dkAWEZgTC7JzwAQqyREBvQDRBIrC3DFygbAMgJjzNHycxRQ+FBQV1m1c5vm+QCRJFpIzp0n4PJ7voDLZg8bk73yMwOAyFlABvcDgLcswBWaDUBq+gfAyQgsCzAmO+VnBgDRs4AM1pZy0qEZQaqyAWNMdsnfAACty0K4oq0S1s4F4iJJdlEZd+awO3vYzQb8qCVkl4GMyU75HQAiZQFl1c5tpEtB7VwgLhI/FpXxOxuwInLGZKf87QNwhfcFxOoHyBCJ9AVEEm20UFv6BcKHh3phQ0iNSS8LABmwSEyivM4L8MKdO9DWeQPJBCPrPDYmvSwAROOuEVBWndbSEKnmZgNzlm30dZUxL+JlDZYhGJNaFgAiKat2bve869yGB4AsGxLqxZgBPand2ciqrf7PHo4mXtstQzAmtfK7Ezia0DUCMnREkB+dwaFStcpYW7gZgvtlI42M8ZcFAHA6gps+aj0ktKzauc3QEUFu6ehkhoVGkqpVxpLlLk7jftlII2P8ZQEAnGv5lZe0HhI6aBKUVKSnTR7UnNKLqUN7+5YJhK4ylkmZgCs0I7BswJi2sz4ALzK0QFyquCODanc2BvsFQp9rrz6CcKF9BtY/YEzbWQYQT1m1c5vhcwP85GYC15zVn/KTioKPZ8qlIWOMPywDCJVogTjIiCGhoctIxls+MhHxCspB+jJBElGgAAAQeElEQVQCW5fAmLazABAqmUlhGTAkNPSDb8U7u1ISDODoBWcgdYvOeGHrEhjTdhYAvIrVD5ABWYArNBj4NVvYFSkjaO/JY9G4BeksCzDGO+sDCBepQmhZtXMba62ANA4Hjca9NOR++TVc1OVkBGREv4AVpDMmcRYAwkWqEOoOB3XLQ0SbGJbGMtGRuMNE3S8/J47BZ5PH3Cwg3WyYqDGJsQDgVVm1MzN4z7vRJ4ZlYBYQKjQj8HPyGGROFmCTxozxzlMfgIhMAu4ECoAHVHVu2PM3A9cALUA9cLWqRlhpJYsNmuR8/Wlm9G0yqC8gklT0D4wvL2HV1t2tRgd5kcoRRKGjhGyEkDGRxc0ARKQAuAc4DxgMTBORwWGbrQVGqGoV8BQwz++GtqtoK4XFkwVZgMvPbCB8AXovUj2nwLIBY+LzkgGMBLao6jYAEVkEXAhscjdQ1RUh278GfM3PRra7eMNBc2BmcKSho8kOGw0fHeRFpDkF4dI569iYfOAlAPQGPgi5XweMirH9N4A/RnpCRGYAMwBOPvlkj03MMGXVTgDYtjKrA0AoNxj4PWw0lvA5BeH8nGMQbd0BuzRk8p2XACARHtOIG4p8DRgBjIv0vKreD9wPMGLEiIj7yBhuhdBDzZFnBsdaMCZL+b3GQCzxsgY/Zx1HOx+bPGbynZcAUAeUhtzvA+wI30hEzgH+Axinqp/607w0cjtyI10KKqt2biMtGOOuE1DYLWM7g6Pxc6nJtmqPWcdWTsLkOy8B4A1goIj0Bz4ErgC+ErqBiAwDfglMUtXcGoAdrT5QtBFB7nZpLg+RLLdz2O8yEomKV4fIj/4BKydh8l3cUUCq2gLcADwD1AJPqOpGEZkjIlMDm/0MOB54UkTWicjSlLW4vUWaGOZFBk4M88LvNQb8EjrSKFUjiNyMwCaRmXzhaR6Aqi4Hloc9Nivk+3N8blf2iDYiKAOKxLVFpmQCrtCMIFVVSd3ztEzA5AubCexFtHkBZdXObQ6uFeBmAsVFhSmtJ5SM8HkHtk6BMckR1fQMxhkxYoSuWbMmLcdOyt+ectYADvenmU4WcOYNrbOAvdudy0dZ2CEczdJ1H1JclMQlsRRyq5Jec1Z/XzqJP9j7SatLYNZBbDKFiLypqiP82JdlAG1VVu3cZujC8X5LdYXRZPhdjyh8MXqbTWxyla0H0FbuvIAcmB3sRfh/wZkwbDRSPSK/ZxFHm0yWDMsmTKawAOBVpOGgrrLqnJsd7FV7Th6LJXTeQCrmDPh5ftbJbDKFXQLyauAEKCpxRvaEdwi76wW4WUAeqTmlF8VFhWm/FOQuZD9ryhD69ugSzAbmLNuYEWsVGJOJLAAkYuAEqLwk8ryAsmrnNlJfgDs7OAvnBXhRc0qvjJo30B5zBtrC5huYTGGXgPzi9gVEkuWzg71wO4cTkao5BvHmDLjSVW3U5huYTGEBwPgimU7N9uhAjlZ1NBX9BIlKtmPZOpGNXywAJCNWh3AOVglNlXhZgx8ZQrSqo6maTZyIZM/NMgfjFwsAyRg4wbmeH14uuqzauY1UJdS0Eu+/WHehGpefl4zao9qoMZnOZgK3VaQZwtFmB7szg105NEO4PaRyFrKbEfTt0SXjVyKLNFPZD3ZpKTv4ORPYMoBUKKuOPC8g/JJRDncKp0LoJSO/O5DdjCAbMoFUzbmwS0v5xzKAtvr7n6Hxo9Yf7n+a6VwKOrF/9P6AvduduQWWBSRsxTu7jirP4FdACM0EQmV6VuAHPzILyyJSzzKATBJtAfmyauc2Vn9AlpeMTqdUlaSINGooG7ICP/gRQC2LyC4WAPyQ6Kphxnd+laSINGoo1lyCZORDNmGygwUAP7ijgqINDY3FXTnMLgO1Sc0pvY4aNZTKEUNtkevZhF9F8+xSUvuwAOCXaJeCIPbcALsM5JvQDww/J5lFm0uQjNBsIhczAb+Crl1Kah8WAPwU6VJQWbVza3MD2lUqRwy1RTaNNjK5z0YB+S3WymEAk+a2fi4HVw7LJOEjhlzpDAzRRht5kYuZQ7gP9n5CcVGhXQaKwEYB5Zo8KBaXTtE+RNK5mE2y/Qr5kjmUdj/OLgO1AwsAfivs1rpEhMvtC3CVVR99ScgtG53IsSxjSFp4LaL2zAiS7Vfwe0RSqEzLLPxchc0PudgxbQHAb+4HcviooLLqo7eL1CeQ6AgiyxjaJPyP2R1FlEl9BuH8HJEUKhMzi0z7GeRiRmIBIFXCRwW58wJcbr2gPFhHOFu4ASGT/usM5+eIpFB+ZRaZlkWY2CwApFKsy0Fl1W1fRzj8kpFdEvJFpDLVmZwV+MGPzCITswg/tfWSVCZeQrJRQO0h1sigSFVDk9X0kbNkpfFdtJFEycjVYNKWkU3hcjGTqG9sZurQ3m3ej40CyhVl1W3PAky78PM/t0y+xNQWfvVP5HomkUksALSHaJeC3HWEbRWxvOJeYsq1TMCv/ok5yzZSu7OR52s/siCQYhYA2oN7XT7ZqqEmp7jZRPiKZ17lWuAIN2ZAT2p3NrJqq2UBqWYBoD25mYDrUPPRVUPD5wnEU1ZtASOLJXtZKdnAEU+mBJbx5SWs2ro7ZfMdkpWL/RKeAoCITALuBAqAB1R1btjzxwK/AU4HGoDLVfU9f5uaA8JH6ITOFSirTmxfljHkrVSNJPEjsPgVRFI13yFZudovEXcUkIgUAJuBCUAd8AYwTVU3hWxzPVClqteKyBXARap6eaz95tUooFiirSgWT+iKY67PD4cL7/a3fcYkwK/RUpmSjbj8GOFU1bsb/33JaW1uS3uPAhoJbFHVbYGDLwIuBDaFbHMhMDvw/VPAQhERTdcY02ziriUQbb5ANGXVR9/f8y4c9meYojHJ8is7WfHOLl8W+PFLWzOS7Q37OdRyxL8G+cRLAOgNfBByvw4YFW0bVW0RkX1ADzIph8tkoeUjvJZ3+PwwKC7/LGD8aaYFAJMzwhf4SbeqPt0Y9Lnjkw5Ic5ZtzNoAIBEeC//P3ss2iMgMYAbAySef7OHQeSbRWbyhAaPoc1BwrP9tMiZNMm3WbFsCUq+iY+lU0IFXtzYEHztzQA+/mpY0LwGgDigNud8H2BFlmzoR6Qh0A/aE70hV7wfuB6cPIJkGmxChAcNmABuTUm0JSH7MAE6FDh62eQMYKCL9ReQY4Apgadg2S4ErA99fArxg1/+NMSazxc0AAtf0bwCewRkG+qCqbhSROcAaVV0K/Bp4VES24Pznf0UqG22MMabtPM0DUNXlwPKwx2aFfN8MXOpv04wxxqSSl0tAxhhjcpAFAGOMyVMWAIwxJk9ZADDGmDxlAcAYY/JU2paEFJF6YHuSL8+sUoHtz87fzt/OP3+doqpFfuwobesBqGpxsq8VkTV+VcPLRnb+dv52/vl9/n7tyy4BGWNMnrIAYIwxeSpbA8D96W5Amtn55zc7//zm2/mnrRPYGGNMemVrBmCMMaaNsi4AiMgkEXlHRLaIyMx0tydVROQ9EfmbiKxze/1F5EQR+bOI/D1w2z3wuIjIXYH3ZL2IDE9v6xMnIg+KyC4R2RDyWMLnKyJXBrb/u4hcGelYmSbKuc8WkQ8DP/91IjI55LnvBc79HRE5N+TxrPzbEJFSEVkhIrUislFEbgw8ni8//2jnn/rfAVXNmi+cctRbgTLgGOAtYHC625Wic30P6Bn22DxgZuD7mcB/B76fDPwRZ2W20cDqdLc/ifMdCwwHNiR7vsCJwLbAbffA993TfW5Jnvts4JYI2w4O/N4fC/QP/D0UZPPfBnASMDzwfRGwOXCe+fLzj3b+Kf8dyLYMILhAvaoeBNwF6vPFhcAjge8fAb4U8vhv1PEacIKInJSOBiZLVV+i9SpyiZ7vucCfVXWPqu4F/gxMSn3r2ybKuUdzIbBIVT9V1XeBLTh/F1n7t6GqO1X1r4HvG4FanHXG8+XnH+38o/HtdyDbAkCkBeozc621tlPgWRF5M7CWMkCJqu4E55cGcNeoy9X3JdHzzbX34YbAJY4H3csf5Pi5i0g/YBiwmjz8+YedP6T4dyDbAoCnxedzxBhVHQ6cB3xTRMbG2Daf3heIfr659D78AhgADAV2Aj8PPJ6z5y4ixwNPAzep6sexNo3wWNa/BxHOP+W/A9kWALwsUJ8TVHVH4HYX8Huc9O4j99JO4HZXYPNcfV8SPd+ceR9U9SNVPayqR4Bf4fz8IUfPXUQ64Xz4Paaqvws8nDc//0jn3x6/A9kWALwsUJ/1ROQ4ESlyvwcmAhtwztUd2XAl8IfA90uBrwdGR4wG9rmpc5ZL9HyfASaKSPdAujwx8FjWCevDuQjn5w/OuV8hIseKSH9gIPA6Wfy3ISKCs654rareHvJUXvz8o51/u/wOpLsHPIke88k4veRbgf9Id3tSdI5lOD34bwEb3fMEegDPA38P3J4YeFyAewLvyd+AEek+hyTO+XGcNPcQzn8y30jmfIGrcTrFtgBXpfu82nDujwbObX3gj/ikkO3/I3Du7wDnhTyelX8bwFk4lyrWA+sCX5Pz6Ocf7fxT/jtgM4GNMSZPZdslIGOMMT6xAGCMMXnKAoAxxuQpCwDGGJOnLAAYY0yesgBg8pqIDBORB9J07OdCpvcb0+4sAJh8933g7lTtXEQ6xnj6UeD6VB3bmHgsAJisICJLAoXxNoYUx0NEviEim0VkpYj8SkQWBh4vFpGnReSNwNeYCPssAqpU9S0R6RCoIV8ceK5DoKZ6z2j7EpGRIvKKiKwN3J4SeHy6iDwpIv+LU9DvJBF5KVDTfYOInB1owlJgWmrfOWOii/XfiTGZ5GpV3SMinYE3RORpnHroP8Cppd8IvIAzexrgTuAOVX1ZRE7GKQlQHrbPEQSm16vqERH5LfBVYAFwDvCWqu4Wkf+Jsq+3gbGq2iIi5wA/Ab4c2PeZOMFlj4j8G/CMqv6XiBQAXQLH3BuYzt9DVRv8fbuMic8CgMkW3xaRiwLfl+LUP/kc8KKq7gEQkSeBQYFtzgEGO2VWAOgqIkXq1Ft3nQTUh9x/EKfezAKckgIPxdoX0A14REQG4kzl7xSyrz+77cKp0fJgoODXElVdF7LdLuDzgAUA0+4sAJiMJyLVOB/CZ6rqfhFZCRQSufytq0Ng+wMxtjkQ2A8AqvqBiHwkIl8ERuFkA1H3JSJ3AytU9aJAHfeVIU9/ErLflwLlvM8HHhWRn6nqbwJPFwbaYUy7sz4Akw26AXsDH/6n4iwDCE4FxHGB6o8d+ezyC8CzwA3uHREZGmG/tcD/C3vsAeC3wBOqejjOvroBHwa+nx6t8SLSF9ilqr/Cqfo4PPC44GQx70V7rTGpZAHAZIM/AR1FZD3wI+A1AFX9EOe6+2rgOWATsC/wmm8DI8RZTWkTcG34TlX1baCbW3o7YClwPJ9d/om1r3nAT0VkFc56rNFUA+tEZC1OkLoz8PjpwGuq2hL/LTDGf1YN1GQ1ETleVZsCGcDvgQdV9fcJvP47QKOqPhC4PwKnw/fs2K9sOxG5E1iqqs+n+ljGRGIZgMl2s0VkHc5onneBJQm+/hfApwAiMhNnVabv+drC6DbYh79JJ8sAjDEmT1kGYIwxecoCgDHG5CkLAMYYk6csABhjTJ6yAGCMMXnKAoAxxuSp/w9GI5TPC0gKHgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "km.fit(true_ages_at_time_of_death, np.ones(number_of_vampires), label='true lifetimes')\n",
    "ax = km.plot()\n",
    "\n",
    "km.fit(ages_with_censoring, np.ones(number_of_vampires), label='observed lifetimes (some censored)')\n",
    "plot = km.plot(ax=ax)\n",
    "\n",
    "plot.set_xlabel('age (years)')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Removing all censored data doesn't help!\n",
    "\n",
    "What if we simply remove all data points of vampires that are still undead? Then the distribution of all lifetimes should match the actual distribution obtained using Alice's predictions much better, right? Wrong!\n",
    "\n",
    "We plot the two survival curves we saw before, but now we add a third: only including the final ages for those vampires that have died already. Can you think of a reason why this doesn't make the problem of skewing the distribution towards lower lifetimes any better?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "ages_with_censoring_observed = [age for i, age in enumerate(ages_with_censoring) if observed[i]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,0,'age (years)')"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEKCAYAAAAb7IIBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8lNW5wPHfIQTCMmEJASmJBCk7CSEGAkUgCCgiYrUqUm0F61Wxti5Xr3jbIqK1uFRxodpqXepSENsicmldQQUVAUEUEATZAighgZhAyPrcP2ZhMpktk3cymZnn+/nkk8zMO/Oed2YyzzzvOec5RkRQSikVf1pEugFKKaUiQwOAUkrFKQ0ASikVpzQAKKVUnNIAoJRScUoDgFJKxSkNAEopFac0ACilVJzSAKCUUnGqZaR23KVLF8nIyIjU7pVSKipt2LDhiIikWvFYEQsAGRkZrF+/PlK7V0qpqGSM2WvVY+kpIKWUilMaAJRSKk5pAFBKqTgVsT6AeFdVVUVBQQEnT56MdFOUUs1QUlISaWlpJCYmhm0fGgAipKCgAJvNRkZGBsaYSDdHKdWMiAhFRUUUFBTQq1evsO0n4CkgY8yzxpjDxpgvfdxujDGPGWN2GmM2G2NyrG9m7Dl58iQpKSn64a+UqscYQ0pKStjPEATTB/A8MMnP7ecBfRw/1wJPNr5Z8UE//JVSvjTF50PAU0Ai8oExJsPPJhcCfxP72pKfGGM6GmO6i8ghi9pYx83PTORgbSG5tuH8z/S/hGMXSikVF6wYBdQD2O92ucBxXT3GmGuNMeuNMesLCwtD3uG21jWsL/005PsrOHbsGH/605/C9vjPP/88N954IwBPPfUUf/vb3wD46quvyM7OZujQoezatcvn/e+77746l3/0ox+Fra1KxSsrAoC3PMXrSvMi8hcRyRWR3NTU0GYyL7jmbQZUJKCL2TeOvwBQU1Nj6b6uv/56fv7znwOwdOlSLrzwQjZu3Ejv3r193sczAHz00UeWtkkpZU0AKADS3S6nAQcteNwAhJXbD4d/NzFq9uzZ7Nq1i+zsbG6//XZWrVrFuHHj+OlPf0pmZiZ79uxh8ODBru0feugh5s6dC8CuXbuYNGkSZ555JqNHj+arr77yu6+5c+fy0EMPsWLFChYsWMAzzzzDuHHjAHjppZcYPnw42dnZXHfdddTU1DB79mzKy8vJzs7miiuuAKB9+/YArFq1irFjx3LZZZfRt29fZs+ezcsvv8zw4cPJzMx0ZRWFhYX85Cc/YdiwYQwbNow1a9YA8P7775Odne3KQkpLSy19XpWKJlYMA10G3GiMWQTkASXhOv/vqbS8qil2E3Z3v7GFrQe/t/QxB/4gmbsuGOTz9vnz5/Pll1+yadMmwP7B+umnn/Lll1/Sq1cv9uzZ4/O+1157LU899RR9+vRh7dq13HDDDbz33nsB2zR58mSuv/562rdvz2233ca2bdtYvHgxa9asITExkRtuuIGXX36Z+fPn88QTT7ja5unzzz9n27ZtdO7cmTPOOINrrrmGTz/9lEcffZTHH3+cBQsWcNNNN3HLLbdw1llnsW/fPs4991y2bdvGQw89xMKFCxk1ahRlZWUkJSUFbLdSsSpgADDG/B3IB7oYYwqAu4BEABF5ClgBTAZ2AieAmeFqrKeT1TWs3H6Ycf26NtUuY9rw4cMDjjkuKyvjo48+4tJLL3VdV1FREdL+3n33XTZs2MCwYcMAKC8vp2vXwK/lsGHD6N69OwC9e/fmnHPOASAzM5OVK1cC8M4777B161bXfb7//ntKS0sZNWoUt956K1dccQUXX3wxaWlpIbVdqVgQzCig6QFuF+CXlrUoSLVi2PFtGUktE5p615bz9029KbVr1871d8uWLamtrXVddo5Hrq2tpWPHjj6/nTeEiHDVVVfxhz/8oUH3a926tevvFi1auC63aNGC6upqVzs//vhj2rRpU+e+s2fP5vzzz2fFihWMGDGCd955h/79+zfySJSKTlFZCygxwd7st77a7coCVMPYbDa/57+7devG4cOHKSoqoqKiguXLlwOQnJxMr169WLJkCWD/EP/8889DasP48eN57bXXOHzY/voVFxezd6+90m1iYiJVVaGf4jvnnHN44oknXJedAWvXrl1kZmZyxx13kJubG7D/QqlYFpUBoFVCC1oYoVZqSO/UjsLSkxoEGiglJYVRo0YxePBgbr/99nq3JyYmMmfOHPLy8pgyZUqdb8kvv/wyf/3rXxkyZAiDBg3i9ddfD6kNAwcO5N577+Wcc84hKyuLiRMncuiQvfvo2muvJSsry9UJ3FCPPfYY69evJysri4EDB/LUU08BsGDBAgYPHsyQIUNo06YN5513XkiPr1QsMJEaTpmbmyuhLggz8/lcyioqOfT9vTx48SgA9h897jodZGuT2Oz7BbZt28aAAQMi3QylVDPm7XPCGLNBRHKtePyozACcCooSeHfbdwCkd2pHqi2JVFtSzIwOUkqpcIraANDC2DOXNbuORLglSikVnaI2ACSYWnomH2PboVJXFqCUUip4URsAAAZ3sXf8ahaglFINF9UBILPbAU7vQr0sQIeGKqVUYFEdALq27ECfHvZZqO5ZQHqndtoRrJRSAUR1AADIyqhmQHcbe4tOMG/5FuYt36J9Ao3gWQSuucjPz8fbsGH36ydPnsyxY8cA+zyAAQMG+J1HsGnTJlasWOG6vGzZMubPn29xy31bunQp8+bNa7L9RTv3EuNPPPEEzz33XIRbFP1iYk3gUb27APYMYG/RCeAIWWkdItomdUp1dTUtW4b/reb+Yf6nP/2Jf//7335rG23atIn169czefJkAKZOncrUqVPD3k6nBx54gGXLljXZ/pqjmpoaEhIaXs7l6quvZtSoUcyc2WSlx2JS1GYAX7Rqwe7D71NVW0G3bruZM2UQc6YMomdKW0D7AYLx8MMPM3jwYAYPHsyCBQtc11dXV3PVVVeRlZXFJZdcwokTJwB7HZ2BAweSlZXFbbfdBvguuzx37lyuvfZazjnnHH7+85+Tl5fHli1bXPvIz89nw4YNHD9+nKuvvpphw4YxdOhQ16zi8vJyLr/8crKyspg2bRrl5eUBjycjI4MjR45w/fXX88033zB16lQeeeQRr/uorKxkzpw5LF68mOzsbBYvXlznG+aMGTOYNWsW48aN44wzzuD999/n6quvZsCAAcyYMcO1z7feeouRI0eSk5PDpZdeSllZmc/nyt2OHTto3bo1Xbp0AWDJkiWuGcpjxowB7PWXZs6cSWZmJkOHDnUVunv++ef58Y9/zAUXXECvXr144oknePjhhxk6dCgjRoyguLgYCK5sd1lZmWsfWVlZ/OMf//B7XBkZGdx1113k5OSQmZnpekxvZbZFhNtvv53BgweTmZnJ4sWLAeqVHgfvZcEBnnvuOfr27cvYsWNd7y2Atm3bkpGRwaef6sJQjSIiEfk588wzJVSvvnmzDH5+sPz0L4PkufVvy4KPFstHO4/IRzuPyKQF78ukBe/LRzuPyOsbC0LeR7ht3br11IUVd4g8O9nanxV3+N3/+vXrZfDgwVJWVialpaUycOBA+eyzz2T37t0CyOrVq0VEZObMmfLggw9KUVGR9O3bV2pra0VE5OjRoyIiMn36dPnwww9FRGTv3r3Sv39/ERG56667JCcnR06cOCEiIg8//LDMmTNHREQOHjwoffr0ERGRO++8U1588UXXY/bp00fKysrkj3/8o8ycOVNERD7//HNJSEiQdevW1TuOsWPHuq7v2bOnFBYW1vvb1z6ee+45+eUvf+l6LPfLV111lUybNk1qa2tl6dKlYrPZZPPmzVJTUyM5OTmyceNGKSwslNGjR0tZWZmIiMyfP1/uvvtun8+Vu2effVZuvfVW1+XBgwdLQUFBne0feughmTFjhoiIbNu2TdLT06W8vFyee+456d27t3z//fdy+PBhSU5OlieffFJERG6++WZ55JFHRETk7LPPlh07doiIyCeffCLjxo2r147/+Z//kZtuusl1ubi42OdxOZ/Xxx57TEREFi5cKL/4xS9ERGTKlCmu90xpaalUVVXJa6+9JhMmTJDq6mr59ttvJT09XQ4ePCgrV66Utm3byjfffCMi9v+FKVOmSGVlpYiIzJo1S1544QU5ePCgpKeny+HDh6WiokJ+9KMf1Xm97r33XnnooYfqHVMsqfM54QCsF4s+h6PyFNCl5zzC688MRoB2Rz6nJDk94H1UXatXr+aiiy5yVQC9+OKL+fDDD5k6dSrp6emMGmUvsXHllVfy2GOPcfPNN5OUlMQ111zD+eefz5QpUwDfZZfBfkrFWY3zsssuY+LEidx99928+uqrrnLSb731FsuWLeOhhx4C7N969+3bxwcffMCvf/1rALKyssjKygr5WH3tI5ALLrgAYwyZmZl069bN9W110KBB7Nmzh4KCArZu3ep6riorKxk5ciTJyclenyt3hw4dwn1VvFGjRjFjxgwuu+wyLr74YsD+Gv3qV78CoH///vTs2ZMdO3YAMG7cOGw2GzabjQ4dOnDBBRcA9pLYmzdvDrps9zvvvMOiRYtclzt16sTy5cu9HpeTs31nnnkm//znP13t9yyzvXr1aqZPn05CQgLdunVj7NixrFu3juTk5Dqlx32VBV+7di35+fmu52natGmu4wfo2rWrFvNrpKgMAKcYEqpPRLoRjXde03U8OomfGlDGmHqXW7Zsyaeffsq7777LokWLeOKJJ3jvvfd8ll2GuuWle/ToQUpKCps3b2bx4sX8+c9/drXjH//4B/369QvYjlD52sfatWv93s+9zLRnCerq6moSEhKYOHEif//73+vd19tz5a5NmzaUlJS4Lj/11FOsXbuW//u//yM7O5tNmzb5fY0ClcQOtmy3iNR7nkXE53G57zshIcFVfttbmW1/7Xd/b4iPsuBLly71+x44efKk1/edCl7U9gG4q6qt4MviT1yXdXZwYGPGjGHp0qWcOHGC48eP869//YvRo0cDsG/fPj7++GMA/v73v3PWWWdRVlZGSUkJkydPZsGCBa4PFl9ll725/PLLeeCBBygpKXF9mz733HN5/PHHXR8WGzdudLXv5ZdfBuDLL79k8+bNIR+rr30EKokdyIgRI1izZg07d+4E4MSJE+zYscPnc+VuwIABrvuB/Xx9Xl4e8+bNo0uXLuzfv7/Oc7Bjxw727dvnNVB6E2zZbs/X7+jRoz6Pyx9vZbbHjBnD4sWLqampobCwkA8++IDhw4fXu6+vsuB5eXmsWrWKoqIiqqqqXMfitGPHjmY5Yi2axEQA6JLUnRPVxwHniCD7vADtCPYtJyeHGTNmMHz4cPLy8rjmmmsYOnQoYP9weuGFF8jKyqK4uJhZs2ZRWlrKlClTyMrKYuzYsTzyyCOA77LL3lxyySUsWrSIyy67zHXd7373O6qqqsjKymLw4MH87ne/A2DWrFmUlZWRlZXFAw884PWDI1i+9jFu3Di2bt3q6gRuqNTUVJ5//nmmT59OVlYWI0aM4KuvvvL5XLkbM2YMGzdudAWl22+/nczMTAYPHsyYMWMYMmQIN9xwAzU1NWRmZjJt2jSef/75Ot/8AwmmbPdvf/tbjh496uqAXrlypc/j8sdbme2LLrqIrKwshgwZwtlnn80DDzzAaaedVu++vsqCd+/enblz5zJy5EgmTJhATk5OnfutWbOGCRMmBP18qPqishw0wJXPDEYwTOz4Y+g1lpLKYoZ3HQ/AvOX20SZzpgyisPQkU7N7WNJmK2k5aHXTTTdxwQUX6IdYCDZu3MjDDz/Miy++GOmmhJWWg/YjwRhatojqQ1Bx7H//939dQ2xVwxw5coR77rkn0s2IelHeCQwtait1JJCKSt26dWvSiWexZOLEiZFuQkyI+gCA7TSOlx6hqn1XPj38LgDHq9rRLtEW4YYppVTzFvXnT/p0bU8rqulZdpgOrTrToVVnaqXGNRLoZHUNyzYd0M5gpZTyEPUBAHBkAcdcF/ul2ccmr9l1xLVUpFYHVUqpuqI6AGxsKWzY8Tp9urav0xnsrBCqlFLKt6gNALlJ3QD4qGxPUNvrnIDgaDnopp+V7a6iooIJEyYEnJvQXF8nT3PnznWV4IikOXPm8M4770S6GT7fx5EStZ3AfdNHM3TnP12XW7VswZHySgCKT9aveZLeqR2FpSebrH3qFC0HHbyNGzdSVVUVsISDahhf6y6EWo46VkRtBuCpV8dERiXtZdQPU0lMsKaGTKzTctCRKwddXFzMj3/8Y9dM282bN3P48GGuvPJKNm3aRHZ2Nrt27apznw0bNjBkyBBGjhzJwoULXdfX1NRw++23M2zYMLKyslx1lsC+5kBmZiZDhgxh9uzZgO8y0W+88QZ5eXkMHTqUCRMm8N139nIq3ko9Azz44IOufd51112uff7+97+nX79+TJgwge3bt9c79pKSEjIyMqitrQXspSbS09Opqqri6aefZtiwYQwZMoSf/OQnrvdesK9H+/bt+e///m9ycnIYP348hYWFrvu/9tprrvfJvHnzOOuss1iyZInP58NbiW534qfcdX5+Ppdccgn9+/fniiuuqFcX6a9//Su33HKL6/LTTz/NrbfeWm8fYWdVWdGG/jSmHLSIyP+tmS8/e3qQ/OyZTPn0s2fsP+//Xj5d9yeZ996LMmnB+9LzjuVyzxtbXKWiF6/bK+999V2j9msV9zKv89fOlxn/nmHpz/y18/3uX8tBR7Yc9I033ihz584VEZF3331XhgwZIiIiK1eulPPPP9/ra5aZmSmrVq0SEZHbbrtNBg0aJCIif/7zn+Wee+4REZGTJ0/KmWeeKd98842sWLFCRo4cKcePHxcRkaKiIhHxXSa6uLjY1eann37aVa7aW6nnN998U/7rv/5LamtrpaamRs4//3x5//33Xe+r48ePS0lJifTu3VsefPDBescydepUee+990REZNGiRa6y0keOHHFt85vf/MZVejqY10NEBJCXXnpJRETuvvvuOq/nkiVLRMT+3rj//vtd+/H1fHgr0e3OX7nr5ORk2b9/v9TU1MiIESNc/yPO92tZWZmcccYZrhLYI0eOlM2bN9fbh5aD9qV7Jmzz+KZv6wblRwEbo3p3YduhUtbsOsL4Afb+Aj0NdIqWg45sOejVq1e7Fl85++yzKSoqqlMd1FNJSQnHjh1j7NixAPzsZz/j3//+t+v4Nm/e7PqGW1JSwtdff80777zDzJkzadvWvkhS586d/ZaJLigoYNq0aRw6dIjKykrX6TNvpZ7feust3nrrLVf9qLKyMr7++mtKS0u56KKLXPv0dUpt2rRpLF68mHHjxrFo0SJuuOEGwF7477e//S3Hjh2jrKyMc889N+jXIzs7mxYtWjBt2jTA/t51lq72tn9nu309H95KdLsLVO46LS0NgOzsbPbs2cNZZ53lum+7du04++yzWb58OQMGDKCqqsp1PE0pegOAw8aWwqqiL8hPOfXkVUkl6d12ex0J5JwXYGuTyLh+XZuyqT7dMfyOJt+naDnoiJaD9vb8+zte8VK22f22xx9/vM6HJcB//vOfevfxVyb6V7/6FbfeeitTp05l1apVzJ07F/Bd6vnOO+/kuuuuq/MYCxYsCOp1mzp1KnfeeSfFxcVs2LCBs88+G7Cfqlm6dClDhgzh+eefZ9WqVa77BHo9vPHVFud709/z4a1Ed0pKiut2f/9D7u1zL5vt7pprruG+++6jf//+EVvaMqg+AGPMJGPMdmPMTmPMbC+3n26MWWmM2WiM2WyMmWx9U+ubaO/zZe2xupUKO7c6zVUd1JPOC7DTctCRLQftfnyrVq2iS5cuJCcn+9xXx44d6dChA6tXrwZw3dd5fE8++SRVVfb39I4dOzh+/DjnnHMOzz77rOs8enFxsd8y0SUlJfToYS+c+MILL7ge31up53PPPZdnn33W1edx4MABDh8+zJgxY/jXv/5FeXk5paWlvPHGG16Pp3379gwfPpybbrqJKVOmuDpiS0tL6d69O1VVVXWOMVi1tbWuTOiVV16p863bG3/Ph7cS3e6CLXftS15eHvv37+eVV15h+vTpDTlMywTMAIwxCcBCYCJQAKwzxiwTka1um/0WeFVEnjTGDARWABlhaG8dF1Ya3jPlnEz4vsH3bY6ZQFNyLwcNuMpB79mzx1UO+rrrrqNPnz7MmjWLkpISLrzwQk6ePImI1CkH/ctf/pKsrCyqq6sZM2aMz5LQl1xyCTfddJOrHDPYSzXffPPNZGVlISJkZGSwfPlyZs2axcyZM8nKyiI7O7vR5aC97WPcuHHMnz+f7Oxs7rzzzgY/rnvZZOdpg3vvvRebzeb1uXI3d+5c1/G1bdu2zgeuL8899xxXX301bdu2rfNt/5prrmHPnj3k5OQgIqSmprJ06VImTZrEpk2byM3NpVWrVkyePJn77ruPl19+mVmzZnHvvfdSVVXF5ZdfzpAhQ5g7dy6XXnopPXr0YMSIEezevRuwf6tfuXIlCQkJDBw4kPPOO4/WrVuzbds210ph7du356WXXiInJ4dp06aRnZ1Nz549XV8qvJk2bRqXXnppnW/599xzD3l5efTs2ZPMzMwGB+h27dqxZcsWzjzzTDp06BBUmW9fz8ftt9/O119/jYgwfvx4hgwZUud+F110ER9//DFDhgzBGOMqd92QVcouu+wyNm3aRKdOnRp0nFYJWA7aGDMSmCsi5zou3wkgIn9w2+bPwDcicr9j+z+KyI/8PW5jy0Gv2L2Cgete4u7iddQktuFXOTfZbyg/ysaWWRR8/x3rtvQE7GWhfYlUuWgtB62U9dq3b+/KSqLBlClTuOWWWxg/frzX25tDOegegHvuU+C4zt1c4EpjTAH2b/+/sqJxgRSm5VCT2IaEqnJS99hPWVBdydAWu1xDQfcWnWDe8i26QphSqtk4duwYffv2pU2bNj4//JtCMAHAWy+KZ9owHXheRNKAycCLxph6j22MudYYs94Ys945Prexqlq3B6BzwWf2K2zdoPLU6mA9U9qyt+gEa3Yd8Xp/nSGsVOyIlm//HTt2ZMeOHfWWuWxqwQSAAsC92H4acNBjm18ArwKIyMdAEtDF84FE5C8ikisiuampqaG12ENVUjLr2yTxemJtvdvGD+jGnCmD6JnS1uf90zu1i1iHcKDTb0qp+NUUnw/BBIB1QB9jTC9jTCvgcmCZxzb7gPEAxpgB2AOANV/xA8jr2B+At70EgOYsKSmJoqIiDQJKqXpEhKKiIpKSksK6n4CjgESk2hhzI/AmkAA8KyJbjDHzsM9IWwb8N/C0MeYW7KeHZkgTfbLlp2Tyxd73Am7nXB/AOSks0tLS0igoKMCqU2FKqdiSlJTkmkwWLkFNBBORFdg7d92vm+P291ZglLVNa4TqSpKPfQWnnwfgdVawO2c/QFMOB01MTPRbqEwppcItZorB1RkJZOtGGyrZf9TeGTx+QDe/6wNEsh9AKaUiJSYCQL2RQEBa57YktQy+zKuOBlJKxZvorgXUqh2UfkdVUjJ7qku5iVJObvkrOacNZ1Bix3qb++sHSO/Ujv1Hj7Ns0wGAuJ0hrJSKH9GdAaTlQstW5HXsT0ZL+ymePdWl9WoDgb0fAPA5HwBO1QnSWkFKqXgQ3QHAIT8lk1sH/YJHq2z0qfFe/c/ZDxDszGBnrSA9LaSUilUxEQC8qaitYvv36+pcF8zMYCetGqqUinVRGwBsrWwUlBb4vD0toR227ze6RgJB3ZnBWiNIKRXvojYAjO4xmtYtW/vewNaN023G60ighmQCSikVq6I2APiSUFVO4kn/6wMEUyPISYeHKqViVUwFgOK0HAASK6yrCKiTxJRSsSqmAkBhxkhqEuuvTeuPc26AP5oFKKViUfQHAMdksFAEMzcANAtQSsWm6A8AjslgoQhUI0gppWJZ9AcA8JsF2NokUlh6ss5w0FDoaSClVKyJjQDgkQVsbCls2PE6AOP6dWVqdg+/heGC6QfQ00BKqVgTGwHAzY/aZwDwUdkeqKmEr9/2u32w/QBKKRVrYi4AnNn3QoZWO+oBtUuFkyV+t3f2AwSTBSilVCyJuQDgi61Nos9+AGcW8Mzq3cxbvsVniQjtB1BKxZK4CQDj+nX12Q8wfkA3rjmrl2tEkK8SEdoPoJSKJXETAAJxlocIVCxOswClVKyI6gBQpyJoIyaEefJXLE6zAKVUrIjqAFCnIqjHUNCEqnK6713LhyU7G/y4DSkWp5RS0SqqA0A9zjWCHYvE9/5uB6WV37uGgvrrCFZKqXgTWwHAkQVUJSWfKgrnNhTUX0ewLzo8VCkVq2IrAAShIaUhfE0Sc64XrGsGK6WiWdwFgGBKQzj5miTmXC9Y1wxWSkWzuAsATsFmAloqQikVq6I+AARaHN6XYDMBLRmtlIpVUR8A6i0O36od1OhpGaWUCiTqA0A9ablgjOUPq6OBlFKxpmWkGxBpzrkB6Z3a+dxmVO8ubDtUyppdRxg/oFud25wjgjwfc1y/rmFpr1JKWSWoDMAYM8kYs90Ys9MYM9vHNpcZY7YaY7YYY16xtpmhSS7aRfq+9X63CWZugL9+APcRQToySCkVTQIGAGNMArAQOA8YCEw3xgz02KYPcCcwSkQGATeHoa0NsrGlsMTWjtMObmZF4Wd8+OljTbZvLRinlIoGwWQAw4GdIvKNiFQCi4ALPbb5L2ChiBwFEJEm/fTzHAmUZzsDgDc6dqF1i0RSO/WmtLy40fsJth9AC8YppaJBMAGgB7Df7XKB4zp3fYG+xpg1xphPjDGTvD2QMeZaY8x6Y8z6wsLC0FrshedIoPzUofRLCv4cfDA1gho6H8B9trDOGlZKNUfBBABvQ2rE43JLoA+QD0wHnjHGdKx3J5G/iEiuiOSmpqY2tK3Ba+BIoIb0AzQkC9C+AaVUcxZMACgA0t0upwEHvWzzuohUichuYDv2gNB8BFggPhxZgFJKNWfBBIB1QB9jTC9jTCvgcmCZxzZLgXEAxpgu2E8JfWNlQxstwALxjR0NpJRS0SZgABCRauBG4E1gG/CqiGwxxswzxkx1bPYmUGSM2QqsBG4XkaI8PvySAAAdeElEQVRwNToUFbWVAReHcdYHCrZaaEPp6CClVHMS1EQwEVkBrPC4bo7b3wLc6vhpVpKLdpG652PIGEnh9/v9bus+ectzcpcV0ju1o7D0pOWPq5RSoYiZmcC2VjYKy+0jiyqqKwBcK4N1LviMwoyRp/oB+kyMWDu9zRwGnT2slGp6MRMARvcY7fp7xW57slKVlMz3Kb1PbRSgH8Cd83TQyeqaemUi9hadYN7yLYC9Y9izPIQ/vkpO7D96nJXbD2sQUEo1mZgJAFZzfhB7flu3jwSyjwLaW3QCqF8fKBR6ekgp1dRiNwCYhHploStqK1lR+Bl8tAu6ZwL2U0fu2UMg4wd0c33gO7MAqzg7iTULUEo1hdgNAK1tcPJYnavSkrpAUhcoPwpt7BPRnP0GzYFmAUqpphR76wFYLJgJYkopFY00AAQQzAQxpZSKRrF7CshC/haNcR8R5NTQkUHudIEZpVRT0QAQhHH9unodu+8+IsipsSODvAUZ7RdQSoVDXAQA52zgwoyR9W6rqK7gwwMfNmgkkJP7iCCnecu3uCqGWjE8FHR0kFIqPGK6D2B7+Xe80jUNsM8G9ibNlkZpZall+wxHxVBdYEYpFQ4xmwHkdc9j+9HtvJNQyU/dZwMDVFdCwXr7ugFBCmbxeLBnBWt2HQlL34BmAUopK8VkBmBrZeOHHX9IvzY+Pmxt3aCyYUM7GzIaaFTvLvRMaVvnur1FJxqVFWgWoJSyWkxmAKN7jLbXA3KbDVyvHyCELCBYvvoGGkuzAKWUlWIyA3BpbQNjKE7LATz6AULIAiJNswCllJViOwA4FGaMrFsV1MmZBQSpsbOCg11P2B9dVEYpZZW4CAA+ObKAiuoKVuxewYcHPvS7eWNmBVs1OkizAKWUVWKyD8AXr/MBqitJK/kW0nKDKgwX7GggT/5GB0HDRghpX4BSygpxkwF47QeABvcFNDYL8BwdBA0fIaRZgFLKCnGTARRmjPQ5GaypeBsdBOGZPayUUoHETQbQnIVj9rBSSgUS+wHAJEDpqZE3zn6A5mT8gG4M6G5r0CghHQ2klGqsmA8A+yqPcv+377Gq6Avf/QAOztFAwYwIslpDswDtB1BKNVZMB4C87nmcbjudfRXFrD32lf/5AN+8T1rJt6S2SSW1TarfAnG2NokUlp6ksPSkZauFObOAhnCuHaCZgFIqFDHdCZyfnk9+ej73f3Cn/w1tjo7X8qNBPa778Etv6wQ0FedQVF0vQCkVipgOANHI6iqiSinliwaAZiQcK4wppZQvGgDchbFCaDBCXWHMfR1hXT9YKRWsmO4Edrf9+AFWFX3hfyO3WcHB1gdy7xC2umMYghsdlN6pHam2JFJtSToySCkVtLjIAPJsZ7C9/DvWHvuK6/C/RrBTms2+lGSg+kDevm1b2THsrCEU7ExhzQaUUsEKKgMwxkwyxmw3xuw0xsz2s90lxhgxxkTmHIobWysbBaUFAOSnDqVfkv2DMNBcAEv23ciy0Z4aMkdAswGlVLACBgBjTAKwEDgPGAhMN8YM9LKdDfg1sNbqRoZidI/RtG7Z2n4hLReMAU6tDeBzRnAD1wjwpjEF47wJZaYwnMoGnD86X0Ap5S6YDGA4sFNEvhGRSmARcKGX7e4BHgCa/aB0v1mARSuFWd03EEq9IPdsQDMCpZSnYPoAegD73S4XAHnuGxhjhgLpIrLcGHObhe0LC2dl0GD6AkJldd9AQ/sClFIqkGACgPFynbhuNKYF8AgwI+ADGXMtcC3A6aefHlwLrWIS2Fd+mPt3LQEgsVMSSYldOfu79fQPQwAIh1G9u7DtUClrdoU2L0A7iJVS7oI5BVQApLtdTgMOul22AYOBVcaYPcAIYJm3jmAR+YuI5IpIbmpqauitDkFexgROb93ZdbkqKZmvkpJ4j+N+q4NWVFc0eWE4X0LtC3DSDmKllLtgMoB1QB9jTC/gAHA58FPnjSJSAnRxXjbGrAJuE5HG9aRaLD89n/wqA206ua57eMtfgZN0Lvis3jKRzglhaba0oJaKDEaoy0m6a2wW4KTZgFIqYAYgItXAjcCbwDbgVRHZYoyZZ4yZGu4GhlNVUjI1iW3q32BRR7AnK0YHhVI11BvNBpRSQU0EE5EVwAqP6+b42Da/8c0Kk1bt7IvD2CLXgeocHeTpZHVNozKDxnBmA5oJKBVf4mImsEtaLnzzfkSb4OsDtqEjhKwcDaRlpZWKT3FTC8iXjS2F1xNr69/gWCSmsZPCgtWQ2cPhWkNYl5lUKr7EdQDI69gfgHdNOf1WL6w7GsjWzd5hXHk86MJwjdGQ/gFnP4Bz7YBQRgR5o8tMKhVf4i8AOPsBgPyUTAYn2Kht2Zq2JQe8zwyurnQtFelvmUgrNDQL6JnSlr1FJyzNBHSZSaXiR/wFgLRcaNnKdbEqKZkTHX7AiQ49vG8fphFB3jQ0C5gzZRA9U9pa2gbn6CDNBJSKffEXAFRQNBNQKvbF1yigAMJZGyhY3oaJBhoi6r6OsFXrB+vIIKVinwYAYF95Idd1SuLiynac6zkrGE7NDE7pGfa2NLSInPs6wuFYP9h9xjDorGGlYklMBwDnojDO1b28cY4E2lVeyBsdu3DuCW8P1A3Kj4aplYH5KyHhvo7wvOVbLM8GPPepGYFSsSOm+wDqLArjQ35KJnf0vpTT2wQuTtcUw0G9CbZz2DkyCLB8dJCTLjKjVOyI6QzAJz8lIfz1AwS7TnA4BFNIzl82AJoRKKXqiukMwCePoaBOVa3bAz5WCrNgqcjGaGghOfdsAMKbEWgWoFR0is8MALxmAVVJyXyf0tv79hHuB4CGlZN2zwYgPP0DYM8INAtQKjrFZwYAXrOAfeWF3Ni2mhvbVrOq6IsINcy3xpSTbor+AaVUdInfAAB1ykLkdezv6gjeaWr47NtP62/vLBB3KHLBoSHlItw5Zw47Zw87swEragnpaSClolN8BwC3LMA5Gmhe5+H0q6wksaKs/vbOAnFV5U3c0FOsWFTG6mxAi8gpFZ3itw/AyaMvoDBjJDXFXr79NyONXVrS12ihxvQLeE4YC4ZOKlMqsjQANINFYhpqXL+uDf6w9cU5k7ixs4hDCUbaeaxUZMV8ALC1slFYXkhFdYXfGcGeWlRX0G/1QorTcurPCaiphK/fhj4TLW5t03NmA/OWb7F0lbFgBMoaNENQKrxivg9gdI/RTO41OeCMYHdVrduzvXUrft3qOB9+52Xsf7tUOFliYSsbLtTOYF/CtcqYP+4L03v70X4FpcIr5gNAKHJOG06P9j3Y3roV73G87kphTlUn4YvX7JlABFjRGewuXKuMNYaWnVAqvDQAgL0juPxonZXC7uh9KRktbYCPmcGdekL7bhHNBJylo63KBMK1ylioPDMEzQiUspYGALB3BJ8xtt7EsKqkZGoS20SoUYGN69eVqdk9LMsE3FcZa06ZgJN7RqDZgFKNF/OdwE7BlIb2xbNAXEVtJSsO2/sGbFUVjLa0pZHnHBm07VAp2w6V1skGrCohEQr3kUY6gkipxoubDCCY0tDeeCsQl5bUhdRWyaS2Sqa0JvY+iJyZwDVn9WJAd5vr+uZyakgpZY24yQCC4qVA3Jc1pbx4Wi8u9HX6uRkMCXVfRjLQ8pENEaigHEQuI9CVypRqPA0A7jwmheV17M/24wd4O7GWC6t8JEvNYEio+wffyu2HwxIMoO7ykxCeJSiDpesSKNV4GgD8yE/JZO2xr6DkoP8F46tORjwLcHIPBlbNFnbylhE09eQxX5wF6TQLUCp4cdMHAKdmBbv/FJQW1N3IrUKok9+FYsA+JDTCWYA3zlNDzh8rJ45BZCaP+aIF6ZRquLjKAEb3qD9eZ8XuFXWv8FIbyLlQTNuSA77LQzSjLMDJ89twODKCNbuONKsswHmM2iegVGBxlQE0RnFaDic69KBtyQHfE8OaYRbgzj0jsHLyGDSfLEAnjSkVvKAyAGPMJOBRIAF4RkTme9x+K3ANUA0UAleLyF6L2xox+8oLua0N5A0Yw3XbPqhzm/ucAI4Xwke7oHsmtlY2rxlHJIWjf8CZBXiODgpGOEcQaTagVGABA4AxJgFYCEwECoB1xphlIrLVbbONQK6InDDGzAIeAKaFo8FNwm04aF7H/oA9CABc57FpWlIXt/sl20tKtEml0LF9c2Xl0FHP0UHBCPcIIp00plRgwWQAw4GdIvINgDFmEXAh4AoAIrLSbftPgCutbGQ4eZ0h7NYPkJ+SSX5KJvfvWuK62e+IoCjhbehoqIHAc3RQMLzNKfAUyVnHSsWDYAJAD2C/2+UCIM/P9r8A/u3tBmPMtcC1AKeffnqQTQyv0T1G1+8I9qM4LYfkol10LvgsqgOAO2cwsLqT2J9AWYOVGYKvdQf01JCKd8EEAOPlOvG6oTFXArnAWG+3i8hfgL8A5Obmen2MZsNZIbS6ss7M4MKMkXQu+Mz/iKAo1dilJhsiUNZg5axjX8ejp4ZUvAsmABQA6W6X04CDnhsZYyYAvwHGikiFNc2LoLRc+28vy0UWp+UA0LbE/q2yTgCornTcpwZ6TQ53Ky1l5VKTjdUUs461nISKd8EEgHVAH2NML+AAcDnwU/cNjDFDgT8Dk0Qktur0uk0M21deyP27lpDXsT/5Z/2SfqsX1t/emS18v7/+bVHA2TlsdRmJhgpUh8iK/gEtJ6HiXcB5ACJSDdwIvAlsA14VkS3GmHnGmKmOzR4E2gNLjDGbjDHLwtbiMHDOEK43KxjsmUDLVuR17M/pbVLZV15oLw8RiLNIXJSxeo0BqzgXq4HwVSV1ZgS61oCKF0HNAxCRFcAKj+vmuP09weJ2NSnneH1/ncHeRgOB7xFBFW06sOLgGmh5akJSc5wb4EtzyQSc3DOCcFUldR6nZgIqXsRVKYiQeSkTDf5HBKUldQGxzwlwau5zA9w5z4W7VxcF6yuMhqI5VSVVKpppAAiGl/pAcGpEkM95Ac4O4VbtTnUqR5lw1xMKRbirknobNqodxCoWaS2gRnKOCPJaH8jWDdp0gkprq3BGUrgrjIbC6npEnovRa20hFas0A3Dj7AwGqKiuCGr94IBZQIxprhmBZz0iq2cR+5pMFgrNJlRzoQHAjXsHbb0OYS/rBDjF4uzgYDXl5DF/3PsFwtEnYOXxaSezai70FFCw0nKhbSeorYaauqcDCjNG8n1Kb1cWEE/G9etKqi0p4qeCnAvZz5kyiJ4pbV3ZwLzlW3h3m/fArVS80wDQEGm5kNQRTP3qGH77ApydwYe+CHcLI2Jcv67Nat5AU8wZaAydb6CaCz0F5IOzP6ChfQHeHyy6ZwcHw728dLDCNaQ00JwBp0hVG9X5Bqq50ADgQzCTwxqqorbK7+NF00QxT6F0ajZFB7KvqqPNYe5AqB3L2omsrKIBIAT7Ko5y/45XICHRXhcoJdN1m78qoWmtO9WZGOYpmiaKWSFQ1mBFhuCr6mi4ZhM3RKjHppmDsooGgAbK6+5YCqGilH3lh6GmyhUA/FYJVfUE+hYbzlnIOptYKTAikSnLn5ubK+vXr4/IvhviwwMfcuTEEa/9APd/ej+cPMYdfabXub7f6oUkF+1iz5BL6waB0u+gZatTlz1mCBeWFzI5ykpIN6Vlmw6QaksKy2M7M4KeKW2b/Upk+48eD0unu55aig7GmA0iYklpAc0AAmjoimHgZ16ARy0hyo9a0ML4YeU6xp6cGUE0ZALhmnOhp5bijwaAxjIJ9QrFBb1qWHUlFKyP2jpBTc3bOsZOjQ0Izr4CX6OGmntWYAUrZjtrFhFdNAA0Vmtb3dM6DkH1B9i6aRYQonCVpPA2aigasgIrWJFZaBYRXTQABMHfnIB9pfu4v6LUPkM4IdF1vd9Vw5TlrCpJ4W3UkL+5BKGIh2xCRQcNAEHwNSfANSIIoKLUXiIiIZF9juGc7sNDfdLTQJYY169rndNC4Rwx1Bixnk1YVTRPTyU1DQ0AjZCfnk9+ev6pK755H9p0qrdqmN++AD0NZBn3DwwrJ5n5mksQCvdsIhYzAauCrp5KahoaAKzkpWKozg2IjHCOGGqMaBptpGKfBoAGCFgfyMvKYYUZIynMGOm/L8BZLI4a0HkAlvA3YsgpEoEh0GijYMRi5uDpZHUNK7cf1tNAYaYBoAHCUR8IiIticZHk60MkkovZhNqvEC+ZQ3qndnoaqAloALBaq3ZQW82+iuI6fQFt21Zz/vclnOeWCdTrE6iphC9eC35fSR2gz0QrWh2XPGsRNWVGEGq/gtUjktw1t8zCylXYrBCLHdMaAKyWlkuelMGed1yjggB2tGwByR0474R9M699Au1SoX0D/gHLdKGTxvD8Z3aeKmpOfQaerByR5K45ZhbN7TWIxYxEA0AYuEYHOUYFAdy/awkngO05lwKn6gXFwzrC0cIZEJrTt05PVo5IcmdVZtHcsgjlnwaAcGrVzj7Es7qy3k2WrCNcdbLuKSM9JWQJb2Wqm3NWYAUrMovmmEVYqbGnpJrjKSQNAOHknNzlMTIITtULalQW0Kln3ct6SsgS3v5JfY0kCkVzDCZWZBZW9k80x0yisa9ZczyFpAEgBM7hoEDQS0Z6Y0kWoJqEld/cmvMppsawqn8i1jOJ5kQDQAjcl20Makioj5FBAG3Te9Ki+ji1nz1KTtsfMLmrloSIdc5TTM0xE2gMq/on5i3fwrZDpby77TsNAmGmAaApOEcG7X6z3k1VrduTCHzdoobaE7H5zVDV5cwmQj2tFGuBw9Oo3l3YdqiUNbs0Cwg3DQCN5H46yJ3nqaH89HzyTXuoPH5qo+pK1ySwxz97lBbVFfCf2cHv/Ix86DspxJarSAv1tJKV/RHumktgGT+gG2t2HQnbfIdQNcd+icYKKgAYYyYBjwIJwDMiMt/j9tbA34AzgSJgmojssbapzZP76SB3Xk8NeVb8LFjvWkymqnV7EhCoP2DIu+Ld9t8aAOJOuEaSWBFYrAoi4ZrvEKpY7ZcIGACMMQnAQmAiUACsM8YsE5Gtbpv9AjgqIj80xlwO3A9MC0eDY0parisIVCUlc6JVOxg2J7j7/me2PQi4Zww/yIHMS8LTVhXzrAgsVmUnfU9r36w+bK0Y4ZTVowNTs3tY2KrGCyYDGA7sFJFvAIwxi4ALAfcAcCEw1/H3a8ATxhgjkVpxvhnwdmrI64ghZxCorbb/BOuM/LqXi3fbZx4rFUFWZScrtx+2ZIEfqzQ2I9lbdIKq6lrrGmSRYAJAD8C9SlkBkOdrGxGpNsaUACk0pxyuiXk7NfThgQ+99heQ0pOqgna0qDwR/Fj+HwyF1AGn5gL8Z7YGABUzPBf4ibSstA70Pa19yAFp3vItURsAjJfrPL/ZB7MNxphrgWsBTj/99CB2HVt89RcAfFH4Bd8e/5Z1qfYP9GGnDQv8gF+/fSpg2E6DhNZWNFOpZqG5zZptTEDqamtNYkILPt5V5LpuZO8Uq5oWMhPoLI0xZiQwV0TOdVy+E0BE/uC2zZuObT42xrQEvgVS/Z0Cys3NlfXr11twCEopFT+MMRtExJIJQy2C2GYd0McY08sY0wq4HFjmsc0y4CrH35cA78Xz+X+llIoGAU8BOc7p3wi8iX0Y6LMissUYMw9YLyLLgL8CLxpjdgLF2IOEUkqpZiyoeQAisgJY4XHdHLe/TwKXWts0pZRS4RTMKSCllFIxSAOAUkrFKQ0ASikVpzQAKKVUnNIAoJRScSrgRLCw7diYQmBviHdvXqUCm54evx6/Hn/86iciNiseKGLrAYhIaqj3Ncast2omXDTS49fj1+OP7+O36rH0FJBSSsUpDQBKKRWnojUA/CXSDYgwPf74pscf3yw7/oh1AiullIqsaM0AlFJKNVLUBQBjzCRjzHZjzE5jzOzA94hOxpg9xpgvjDGbnL3+xpjOxpi3jTFfO353clxvjDGPOZ6TzcaYnMi2vuGMMc8aYw4bY750u67Bx2uMucqx/dfGmKu87au58XHsc40xBxyv/yZjzGS32+50HPt2Y8y5btdH5f+GMSbdGLPSGLPNGLPFGHOT4/p4ef19HX/43wMiEjU/2MtR7wLOAFoBnwMDI92uMB3rHqCLx3UPALMdf88G7nf8PRn4N/aV2UYAayPd/hCOdwyQA3wZ6vECnYFvHL87Of7uFOljC/HY5wK3edl2oON93xro5fh/SIjm/w2gO5Dj+NsG7HAcZ7y8/r6OP+zvgWjLAFwL1ItIJeBcoD5eXAi84Pj7BeDHbtf/Tew+AToaY7pHooGhEpEPsK8l4a6hx3su8LaIFIvIUeBtYFL4W984Po7dlwuBRSJSISK7gZ3Y/y+i9n9DRA6JyGeOv0uBbdjXGY+X19/X8fti2Xsg2gKAtwXq/T1R0UyAt4wxGxxrKQN0E5FDYH/TAM5FU2P1eWno8cba83Cj4xTHs87TH8T4sRtjMoChwFri8PX3OH4I83sg2gJAUIvPx4hRIpIDnAf80hgzxs+28fS8gO/jjaXn4UmgN5ANHAL+6Lg+Zo/dGNMe+Adws4h8729TL9dF/XPg5fjD/h6ItgBQAKS7XU4DDkaoLWElIgcdvw8D/8Ke3n3nPLXj+H3YsXmsPi8NPd6YeR5E5DsRqRGRWuBp7K8/xOixG2MSsX/4vSwi/3RcHTevv7fjb4r3QLQFgGAWqI96xph2xhib82/gHOBL7MfqHNlwFfC64+9lwM8doyNGACXO1DnKNfR43wTOMcZ0cqTL5ziuizoefTgXYX/9wX7slxtjWhtjegF9gE+J4v8NY4zBvq74NhF52O2muHj9fR1/k7wHIt0DHkKP+WTsveS7gN9Euj1hOsYzsPfgfw5scR4nkAK8C3zt+N3Zcb0BFjqeky+A3EgfQwjH/HfsaW4V9m8yvwjleIGrsXeK7QRmRvq4GnHsLzqObbPjn7i72/a/cRz7duA8t+uj8n8DOAv7qYrNwCbHz+Q4ev19HX/Y3wM6E1gppeJUtJ0CUkopZRENAEopFac0ACilVJzSAKCUUnFKA4BSSsUpDQAqrhljhhpjnonQvt9xm96vVJPTAKDi3f8Cj4frwY0xLf3c/CJwQ7j2rVQgGgBUVDDGLHUUxtviVhwPY8wvjDE7jDGrjDFPG2OecFyfaoz5hzFmneNnlJfHtAFZIvK5MaaFo4Z8quO2Fo6a6l18PZYxZrgx5iNjzEbH736O62cYY5YYY97AXtCvuzHmA0dN9y+NMaMdTVgGTA/vM6eUb/6+nSjVnFwtIsXGmDbAOmPMP7DXQ/8d9lr6pcB72GdPAzwKPCIiq40xp2MvCTDA4zFzcUyvF5FaY8xLwBXAAmAC8LmIHDHGvOLjsb4CxohItTFmAnAf8BPHY4/EHlyKjTH/DbwpIr83xiQAbR37POqYzp8iIkXWPl1KBaYBQEWLXxtjLnL8nY69/slpwPsiUgxgjFkC9HVsMwEYaC+zAkCyMcYm9nrrTt2BQrfLz2KvN7MAe0mB5/w9FtABeMEY0wf7VP5Et8d629ku7DVannUU/FoqIpvctjsM/ADQAKCanAYA1ewZY/KxfwiPFJETxphVQBLey986tXBsX+5nm3LH4wAgIvuNMd8ZY84G8rBnAz4fyxjzOLBSRC5y1HFf5XbzcbfH/cBRzvt84EVjzIMi8jfHzUmOdijV5LQPQEWDDsBRx4d/f+zLAIK9AuJYR/XHlpw6/QLwFnCj84IxJtvL424Dfuhx3TPAS8CrIlIT4LE6AAccf8/w1XhjTE/gsIg8jb3qY47jeoM9i9nj675KhZMGABUN/gO0NMZsBu4BPgEQkQPYz7uvBd4BtgIljvv8Gsg19tWUtgLXez6oiHwFdHCW3nZYBrTn1Okff4/1APAHY8wa7Oux+pIPbDLGbMQepB51XH8m8ImIVAd+CpSynlYDVVHNGNNeRMocGcC/gGdF5F8NuP8tQKmIPOO4nIu9w3e0/3s2njHmUWCZiLwb7n0p5Y1mACrazTXGbMI+mmc3sLSB938SqAAwxszGvirTnZa20Lcv9cNfRZJmAEopFac0A1BKqTilAUAppeKUBgCllIpTGgCUUipOaQBQSqk4pQFAKaXi1P8D4ar666LKXx4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "km.fit(true_ages_at_time_of_death, np.ones(number_of_vampires), label='true lifetimes')\n",
    "ax = km.plot()\n",
    "\n",
    "km.fit(ages_with_censoring, np.ones(number_of_vampires), label='observed lifetimes (some censored)')\n",
    "ax2 = km.plot(ax=ax)\n",
    "\n",
    "km.fit(ages_with_censoring_observed, np.ones_like(ages_with_censoring_observed), label='observed lifetimes of deceased vampires only')\n",
    "plot = km.plot(ax=ax2)\n",
    "\n",
    "plot.set_xlabel('age (years)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kaplan Meier estimate uses information of censored data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,0,'age (years)')"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEKCAYAAAAb7IIBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXl8VNX5/98nG9mGAElYEwlQVpMQAiRgWBJlUUSoVsWtdanFYm1Rq9+iv6pIrVWLiijWqrjUasWlIiJ1QUEElD2i7IsRw5aQkJiQbSY5vz9mYZLMTCbJTJKZed6vV17M3HvuvefMDPe5z3me53OU1hpBEAQh8Ahq7w4IgiAI7YMYAEEQhABFDIAgCEKAIgZAEAQhQBEDIAiCEKCIARAEQQhQxAAIgiAEKGIABEEQAhQxAIIgCAFKSHtdOC4uTiclJbXX5QVBEHySbdu2ndJax3viXO1mAJKSkti6dWt7XV4QBMEnUUr94KlzyRSQIAhCgCIGQBAEIUARAyAIghCgtFsMINAxGo3k5+dTVVXV3l0RBKEDEh4eTkJCAqGhoV67hhiAdiI/Px+DwUBSUhJKqfbujiAIHQitNUVFReTn59OvXz+vXafJKSCl1EtKqQKl1HdO9iul1GKl1EGl1E6lVLrnu+l/VFVVERsbKzd/QRAaoZQiNjbW6zME7sQAXgEudLH/ImCg5W828I/WdyswkJu/IAjOaIv7Q5MGQGu9Dih20WQm8C9t5mugi1Kql6c62JC3rhnNqknDePOuS7x1CUEQhIDAE1lAfYAf7d7nW7Y1Qik1Wym1VSm1tbCwsMUX7Jev6Zx7qMXHC1BSUsKzzz7rtfO/8sor3HbbbQA899xz/Otf/wJg7969pKWlMWLECA4dcv4dPvzww/Xen3feeV7rqyAEKp4wAI78FIcrzWutn9daj9Jaj4qPb1kl85VvbOH7BOXkCoK7uDIAtbW1Hr3Wb3/7W371q18BsHz5cmbOnMmOHTsYMGCA02MaGoCNGzd6tE+CIHjGAOQDiXbvE4BjHjhvk2x4Y0VbXMYvmTdvHocOHSItLY27776btWvXkpOTwzXXXENKSgp5eXkkJyfb2i9cuJD58+cDcOjQIS688EJGjhzJ+PHj2bt3r8trzZ8/n4ULF7Jq1SoWLVrEiy++SE5ODgD//ve/ycjIIC0tjVtuuYXa2lrmzZtHZWUlaWlpXHvttQBER0cDsHbtWiZOnMiVV17JoEGDmDdvHq+//joZGRmkpKTYvIrCwkJ+8YtfMHr0aEaPHs2GDRsA+OKLL0hLS7N5IWVlZR79XAXBl/BEGugK4Dal1JtAJlCqtT7ugfM2SfXp0ra4jNd58INd7D72k0fPOax3Zx645Fyn+x955BG+++47cnNzAfONdfPmzXz33Xf069ePvLw8p8fOnj2b5557joEDB7Jp0yZuvfVWPv/88yb7NG3aNH77298SHR3NXXfdxZ49e1i2bBkbNmwgNDSUW2+9lddff51HHnmEZ555xta3hnzzzTfs2bOHbt260b9/f26++WY2b97MU089xdNPP82iRYuYO3cud9xxB+PGjePIkSNMnTqVPXv2sHDhQpYsWUJWVhbl5eWEh4c32W9B8FeaNABKqf8A2UCcUiofeAAIBdBaPwesAqYBB4EK4EZvdbYhxro61uwrIGdw97a6pF+TkZHRZM5xeXk5Gzdu5IorrrBtq66ubtH1PvvsM7Zt28bo0aMBqKyspHv3pr/L0aNH06uXOc9gwIABTJkyBYCUlBTWrFkDwOrVq9m9e7ftmJ9++omysjKysrK48847ufbaa7nssstISEhoUd8FwR9o0gBora9uYr8GfuexHrlJv6OaHTu+o7uhMwy+tK0v71FcPam3JVFRUbbXISEh1NXV2d5b85Hr6uro0qWL06fz5qC15vrrr+dvf/tbs47r1KmT7XVQUJDtfVBQECaTydbPr776ioiIiHrHzps3j4svvphVq1YxZswYVq9ezZAhQ1o5EkHwTXxSC2j3MLPdijr4HXVFpyQW0AIMBoPL+e8ePXpQUFBAUVER1dXVrFy5EoDOnTvTr18/3n77bcB8E//mm29a1IcLLriAd955h4KCAgCKi4v54Qez0m1oaChGo7FF5wWYMmUKzzzzjO291WAdOnSIlJQU/vSnPzFq1Kgm4xeC4M/4pAHYlR7G933MyUdd+ydx5vhJMQLNJDY2lqysLJKTk7n77rsb7Q8NDeX+++8nMzOT6dOn13tKfv3111m6dCnDhw/n3HPP5f33329RH4YNG8ZDDz3ElClTSE1NZfLkyRw/bg4fzZ49m9TUVFsQuLksXryYrVu3kpqayrBhw3juuecAWLRoEcnJyQwfPpyIiAguuuiiFp1fEPwBZZ7BaXtGjRqlW7ogzI2vjOKK1yowGvvS+447ASjOP0b1+EkAGCJCO3xcYM+ePQwdOrS9uyEIQgfG0X1CKbVNaz3KE+f3aTG4qtogPi2oZXL3YOJCFHr7egBKgsJ8Pi4gCILgbXxyCghAoUktOoxpozm/u65nL3TXbuiu3Qg6U97OvRMEQej4+KwB2G2Zkh58cDufFni2clUQBCEQ8FkDsCNFcbxH93pegBVlrKF83bp26pkgCIJv4LMGAMA0yJy33tALqOjWnW/35bdXtwRBEHwCnzYAdcn9OdoroZEX0N0QTmWNTAsJgiC4wqcNAIAx1axDN/TQdh7YbeSB3UaJCbSChiJwHYXs7GwcpQ3bb582bRolJSWAuQ5g6NChLusIcnNzWbVqle39ihUreOSRRzzcc+csX76cBQsWtNn1fB17ifFnnnmGl19+uZ175Pv4dBooQGHaYKJ2FTDg5FHmfPIsVbWaPQPSIWNQe3dNsGAymQgJ8f5Pzf5m/uyzz/K///3PpbZRbm4uW7duZdq0aQDMmDGDGTNmeL2fVh577DFWrAjsAsba2lqCg4ObfdxNN91EVlYWN97YZtJjfonPegDfRASzLygXU52R6pQ+hCUk0C9SMaD0GCPydtiE4gTnPPHEEyQnJ5OcnMyiRYts200mE9dffz2pqalcfvnlVFRUAGYdnWHDhpGamspdd90FOJddnj9/PrNnz2bKlCn86le/IjMzk127dtmukZ2dzbZt2zhz5gw33XQTo0ePZsSIEbaq4srKSq666ipSU1OZNWsWlZWVTY4nKSmJU6dO8dvf/pbDhw8zY8YMnnzySYfXqKmp4f7772fZsmWkpaWxbNmyek+YN9xwA3PmzCEnJ4f+/fvzxRdfcNNNNzF06FBuuOEG2zU/+eQTxo4dS3p6OldccQXl5eVOPyt79u/fT6dOnYiLiwPg7bfftlUoT5gwATDrL914442kpKQwYsQIm9DdK6+8ws9//nMuueQS+vXrxzPPPMMTTzzBiBEjGDNmDMXF5gX83JHtLi8vt10jNTWVd9991+W4kpKSeOCBB0hPTyclJcV2Tkcy21pr7r77bpKTk0lJSWHZsmUAjaTHwbEsOMDLL7/MoEGDmDhxou23BRAZGUlSUhKbN29u8nchuEBr3S5/I0eO1C3lrY9v18mvJOtr/zFUv/OvZ/WrLz2iN//3E735v5/oFRMu1ismXKy3Pvuq/vTl/7b4Gt5m9+7dZ9+s+pPWL03z7N+qP7m8/tatW3VycrIuLy/XZWVletiwYXr79u36+++/14Bev3691lrrG2+8Uf/973/XRUVFetCgQbqurk5rrfXp06e11lpfffXV+ssvv9Raa/3DDz/oIUOGaK21fuCBB3R6erquqKjQWmv9xBNP6Pvvv19rrfWxY8f0wIEDtdZa33PPPfq1116znXPgwIG6vLxcP/744/rGG2/UWmv9zTff6ODgYL1ly5ZG45g4caJte9++fXVhYWGj186u8fLLL+vf/e53tnPZv7/++uv1rFmzdF1dnV6+fLk2GAx6586dura2Vqenp+sdO3bowsJCPX78eF1eXq611vqRRx7RDz74oNPPyp6XXnpJ33nnnbb3ycnJOj8/v177hQsX6htuuEFrrfWePXt0YmKirqys1C+//LIeMGCA/umnn3RBQYHu3Lmz/sc//qG11vr222/XTz75pNZa6/PPP1/v379fa631119/rXNychr14//+7//03Llzbe+Li4udjsv6uS5evFhrrfWSJUv0r3/9a6211tOnT7f9ZsrKyrTRaNTvvPOOnjRpkjaZTPrEiRM6MTFRHzt2TK9Zs0ZHRkbqw4cPa63N/xemT5+ua2pqtNZaz5kzR7/66qv62LFjOjExURcUFOjq6mp93nnn1fu+HnroIb1w4cJGY/In6t0nLABbtYfuwz7pAVwx5UmGG80SFhGlBx22qevZSwrCXLB+/XouvfRSoqKiiI6O5rLLLuPLL78EIDExkaysLACuu+461q9fT+fOnQkPD+fmm2/mv//9L5GRkYBZdvm2224jLS2NGTNm2GSXwTylYlXjvPLKK20Ccm+99ZZNTvqTTz7hkUceIS0tjezsbKqqqjhy5Ajr1q3juuuuAyA1NZXU1NQWj9XZNZrikksuQSlFSkoKPXr0ICUlhaCgIM4991zy8vL4+uuv2b17N1lZWaSlpfHqq6/yww8/OP2s7Dl+/Dj2q+JlZWVxww038MILL9ieftevX88vf/lLAIYMGULfvn3Zv38/ADk5ORgMBuLj44mJieGSS8xrZFsX87GX7bY+VVt1luxZvXo1v/vdWTHfrl27Oh2XlcsuuwyAkSNH2taNsMpsL168mJKSEkJCQli/fj1XX301wcHB9OjRg4kTJ7JlyxagvvS4vSx4Wloan332GYcPH2bTpk1kZ2cTHx9PWFgYs2bNqtf37t27c+xYm6w95bf4eAxAEVTbMi36DsVFbRd4tKJdaEAppRq9DwkJYfPmzXz22We8+eabPPPMM3z++edOZZehvrx0nz59iI2NZefOnSxbtox//vOftn68++67DB48uMl+tBRn19i0aZPL4+xlphtKUJtMJoKDg5k8eTL/+c9/Gh3r6LOyJyIigtLSswsaPffcc2zatIkPP/yQtLQ0cnNzXX5HTUliuyvbrbVu9DlrrZ2Oy/7awcHBNvltRzLbrvpv/9vQTmTBly9f7vI3UFVV5fB3J7iPT3oADTHVGfm+7OziHz87eYgDq79sxx51fCZMmMDy5cupqKjgzJkzvPfee4wfPx6AI0eO8NVXXwHwn//8h3HjxlFeXk5paSnTpk1j0aJFthuLM9llR1x11VU89thjlJaW2uZ+p06dytNPP227WezYscPWv9dffx2A7777jp07d7Z4rM6u0ZQkdlOMGTOGDRs2cPCg2QutqKhg//79Tj8re4YOHWo7Dszz9ZmZmSxYsIC4uDh+/PHHep/B/v37OXLkiEND6Qh3Zbsbfn+nT592Oi5XOJLZnjBhAsuWLaO2tpbCwkLWrVtHRkZGo2OdyYJnZmaydu1aioqKMBqNtrFY2b9/f4fMWPMl/MIAxITFUlVrXrCkLDkdAMN32yUQ7IL09HRuuOEGMjIyyMzM5Oabb2bEiBGA+eb06quvkpqaSnFxMXPmzKGsrIzp06eTmprKxIkTefLJJwHnssuOuPzyy3nzzTe58sorbdvuu+8+jEYjqampJCcnc9999wEwZ84cysvLSU1N5bHHHnN443AXZ9fIyclh9+7dtiBwc4mPj+eVV17h6quvJjU1lTFjxrB3716nn5U9EyZMYMeOHTajdPfdd5OSkkJycjITJkxg+PDh3HrrrdTW1pKSksKsWbN45ZVX6j35N4U7st1//vOfOX36tC0AvWbNGqfjcoUjme1LL72U1NRUhg8fzvnnn89jjz1Gz549Gx3rTBa8V69ezJ8/n7FjxzJp0iTS09PrHbdhwwYmTZrk9uchNMYn5aABrnsxGWo1l5aOg15pnDH9xNAu5h/IiUVPAdDr+l9yLH0cM9L6eKTPnkTkoIW5c+dyySWXyE2sBezYsYMnnniC1157rb274lW8LQft0x5AkFJElFQ43a+qq4nMlTQxoWNy77332lJsheZx6tQp/vKXv7R3N3we3w4CBwdBiCai9CBnohovAFPXsxdBx060Q8cEoWl69OjRpoVn/sTkyZPbuwt+gW8bAIDwLlRWlGOK6Mqeku0AxOhaglTzqwsFQRACCZ+eAgLo0SuGqNM/0bPyNFEhnYkK6QxAhUnzaUEtxro6VuQelWCwIAhCA3zeANSl9afOEEtlxdmir6hgc27y+lN1xIUoem9fT+1XG5ydQhAEISDxaQOwI0Szbf/79IwJJzjo7FAMIUZSiw4zZv9XtqUipSpYEAShPj5rAEaF9wBgY3keOrITYcVnb/DFw/oDMCJvh22b1AS4h8hBt31Vtj3V1dVMmjSpydqEjvo9NWT+/PksXLiwvbvB/fffz+rVq9u7G05/x+2FzwaBByWOZ8TB/wLmaaDgI8corzYCcHxwf0y59W/2cSGK019tgMGXtnlfAx2Rg3afHTt2YDQam5RwEJqHs3UXWipH7S/4rAfQkPioEFLCTpLSpwshwY31Q0QcrjEiB91+ctDFxcX8/Oc/t1Xa7ty5k4KCAq677jpyc3NJS0vj0KFD9Y7Ztm0bw4cPZ+zYsSxZssS2vba2lrvvvpvRo0eTmppq01kC85oDKSkpDB8+nHnz5gHOZaI/+OADMjMzGTFiBJMmTeLkyZOAY6lngL///e+2az7wwAO2a/71r39l8ODBTJo0iX379jUae2lpKUlJSdTV1QFmqYnExESMRiMvvPACo0ePZvjw4fziF7+w/fbc/T6io6P54x//SHp6OhdccAGFhYW249955x3b72TBggWMGzeOt99+2+nn4Uii2x7tQu46Ozubyy+/nCFDhnDttdc20kVaunQpd9xxh+39Cy+8wJ133tnoGt7GZz0AQiNA14Ey2zAd2x1VUAClZ4j8yTwsaybQ5O5mC2+qquKrF95gWK8Y22mCOxuIdvDltiWPbn6UvcWuS+2by5BuQ/hTxp+c7t+2bRsvv/wymzZtQmtNZmYmEydOpGvXruzbt4+lS5eSlZXFTTfdxLPPPstNN93Ee++9x969e1FK2aZa5s6dyx133MG4ceM4cuQIU6dOZc+ePbZrrF+/noiICJ588kneeustHnzwQY4fP86xY8cYOXIk9957L+effz4vvfQSJSUlZGRkMGnSJP75z38SGRnJzp072blzZyMZAFc899xzfPTRR6xZs4a4uDin11iwYAFbt261aeG88sor9c5z+vRpPv/8c1asWMEll1zChg0bePHFFxk9ejS5ubkkJCTw0EMPsXr1aqKionj00Ud54oknuO222xx+VvY88MADjBgxguXLl/P555/zq1/9itzcXF588UUWLlzIypUrGx1z44038vTTTzNx4kTuvvtu2/alS5cSExPDli1bqK6uJisriylTprB3716WL1/Opk2biIyMtK0TMHv2bJ577jkGDhzIpk2buPXWW/n8888ZN24cX3/9NUopXnzxRR577DEef/xxFi5cyJIlS8jKyqK8vJzw8HA++eQTDhw4wObNm9FaM2PGDNatW0dUVBRvvvkmO3bswGQykZ6ezsiRI+uNIyYmhuHDh/PFF1+Qk5PDBx98wNSpUwkNDeWyyy7jN7/5DWCWqVi6dCm///3v3fo+0tLSOHPmDOnp6Tz++OMsWLCABx98sJ7WkZXw8HDWr18PmLWIHH0eCxYs4OOPP6ZPnz4Ov8P//ve/5Obm8s0333Dq1ClGjx5tMxQ7duxg165d9O7dm6ysLDZs2MC4ceNsx1ofbh577DFCQ0N5+eWX6xnutsJ3DUCvFNhz9km/Lq0/0B8qT1O7voiYUKgwmTOBrAaga/8kSipqCLGT4TVZnhACDXs5aMAmBz1jxoxGctCLFy/m9ttvt0kcX3zxxUyfPh0wywnv3n1WiM+VHPTkyZN58MEHG8lBr1ixwjZPbC8H/Yc//AHwjBy0o2s0hSM5aMAmB52fn2+TTQaoqalh7Nix9eSg7T8re9avX29bfOX888+nqKionjpoQ0pLSykpKWHixIkA/PKXv+R///ufbXw7d+60PeGWlpZy4MABVq9ezY033miTo+7WrVs9mWgr1dVmRd38/HxmzZrF8ePHqampsU2fWaWer732Wi677DISEhL45JNP+OSTT2z6UeXl5Rw4cICysjIuvfRS2zWdTanNmjWLZcuWkZOTw5tvvsmtt94KmIX//vznP1NSUkJ5eTlTp051+/tIS0sjKCjIJht93XXX2aSrHV3f2m9nn4dVovvKK690eB5nctedO3cmIyODhIQEANLS0sjLy6tnAKKiojj//PNZuXIlQ4cOxWg02sbTlviuAbCwI0SztuhbsmPPfni1mIgIqiAypLEOe01tHRsPFhLZKYS0xK5t2VWnuHpS9xYiB92+ctCOPn9X43Uk22y/7+mnn653swT46KOPGh3jSib697//PXfeeSczZsxg7dq1zJ8/H3Au9XzPPfdwyy231DvHokWL3PreZsyYwT333ENxcTHbtm3j/PPPB8xTNcuXL2f48OG88sorrF271nZMU9+HI5z1xfrbdPV5OJLojo2Nte13V67bXjbbnptvvpmHH36YIUOGtNvSlm7FAJRSFyql9imlDiql5jnYf45Sao1SaodSaqdSaprnu9qYyTXmfzeV1J8+iQruQp2uc3hMd0M4XSLDqKg2fyG6uprydeu82s+OiMhBt68ctP341q5dS1xcHJ07d3Z6rS5duhATE2ObtrAeax3fP/7xD4xGcxLE/v37OXPmDFOmTOGll16yzaMXFxe7lIkuLS2lTx+zcOKrr75qO78jqeepU6fy0ksv2WIeR48epaCggAkTJvDee+9RWVlJWVkZH3zwgcPxREdHk5GRwdy5c5k+fbotEFtWVkavXr0wGo31xugudXV1Nk/ojTfeqPfU7QhXn4cjiW573JW7dkZmZiY//vgjb7zxBldffXVzhukxmvQAlFLBwBJgMpAPbFFKrdBa77Zr9mfgLa31P5RSw4BVQJIX+luPmTWKz1UlVcE/Odz/s5OHGLP/KxjWeI7/rCcQRfJPLb8J+Cr2ctCATQ46Ly/PJgd9yy23MHDgQObMmUNpaSkzZ86kqqoKrXU9Oejf/e53pKamYjKZmDBhglNJ6Msvv5y5c+fa5JjBLNV8++23k5qaitaapKQkVq5cyZw5c7jxxhtJTU0lLS2t1XLQjq6Rk5NjWynsnnvuafZ57WWTrdMGDz30EAaDweFnZc/8+fNt44uMjKx3w3XGyy+/zE033URkZGS9p/2bb76ZvLw80tPT0VoTHx/P8uXLufDCC8nNzWXUqFGEhYUxbdo0Hn74YV5//XXmzJnDQw89hNFo5KqrrmL48OHMnz+fK664gj59+jBmzBi+//57wPxUv2bNGoKDgxk2bBgXXXQRnTp1Ys+ePYwdOxYw39D//e9/k56ezqxZs0hLS6Nv3762hwpHzJo1iyuuuKLeU/5f/vIXMjMz6du3LykpKc020FFRUezatYuRI0cSExPjlsy3s8/j7rvv5sCBA2itueCCCxg+fHi94y699FK++uorhg8fjlLKJnfdlHS2PVdeeSW5ubl07do+sxFNykErpcYC87XWUy3v7wHQWv/Nrs0/gcNa60ct7R/XWp/n6rytlYNe9f0qhm35Nw8Wb6E2NILfp88176g8zY9bjLBvL0M/38zBHgPoeftcp+cpqaghIwZipl/c4r60BJGDFgTPEx0dbfNKfIHp06dzxx13cMEFFzjc3xHkoPsA9r5PvmWbPfOB65RS+Zif/n/vic41RWFCOrWhEQQbK4nPM09ZYKohKegkBcMHcbDHAKpqNQ/sNvJpQa3T8+jqakpXfhiQU0GCILQ9JSUlDBo0iIiICKc3/7bAHQPgKIrS0G24GnhFa50ATANeU0o1OrdSarZSaqtSamuhh7JvjJ2iAeiWb1YCxdADHaaJPHWamFAID1bkVWjWn3IcE6iprWOXjiIkPp7aAJwKEgR/wlee/rt06cL+/fsbLXPZ1rhjAPKBRLv3CcCxBm1+DbwFoLX+CggH4hqeSGv9vNZ6lNZ6VLxdKmZrMIZ3ZmtEOO+Hnr3B16UkUhsaSrdQRb9IRVKk86yE7oZwW0BYEAQhkHDHAGwBBiql+imlwoCrgBUN2hwBLgBQSg3FbADaJME+s8sQAD4NdfyEH3T8KHM+edYcDG4C61SQTAcJghAINGkAtNYm4DbgY2AP5myfXUqpBUopa5XHH4HfKKW+Af4D3KDbaLHh7NgURpgcP+Eb00ZS16sPfYqPMvjgdpdxAIDQhARC4uNlOkgQhIDArUIwrfUqzMFd+233273eDWR5tmutwFRDRMVRTBm/xJSRReWSp6BS16sKtqemto7cH093mMIwQRCEtsBvxODqZQIZehBGLQVlVQB0C1W29QEcEahxgOjoaNvrVatWMXDgQLckEhriSWnivLw8lFL1agVOnTpFaGioTajNGW0t57xo0aJ6i7rby1G3hoYy1YLgLfzCADTKBAJio8IICzYPz5hmFqOyXx+gIVYvwEogpYZ+9tln/P73v+ejjz7inHPOae/u0L9//3piaG+//Tbnnntuk8fNmDHDpnjpDs7kA9yloQFYtWoVXbp0adU5QQyA0Hb4tgEIi4KykxjDO7M3PJy5oWU8sWspa4u+rdfMlJHFwR4DbOqgjuhuCKekooaNBwvZeLAwYFJDv/zyS37zm9/w4YcfMmDAAMC5LPD8+fP55S9/yfnnn8/AgQN54YUXGp0vLy+P8ePHk56eTnp6Ohs3bgTck8i1EhERwdChQ20LZyxbtowrr7zStt+ZBLW9nLO7MtUNcSRxfObMGS6++GKGDx9OcnIyy5YtY/HixRw7doycnBxycnKAs3LUeXl5DBkyhJtvvpnk5GSuvfZaVq9eTVZWFgMHDmTz5s2AWS/ovPPOY8SIEZx33nns27fPoUy1M8lsQWgtvi0GlzAKDn9BZpchbK/6Cagiz1SGsWQvGfSq1zQmFHoWHSPyjWeIsEsLNaaNxJRhDl90N4TbtpdU1LTJEABOPPww1Xs8KwfdaegQet57r8s21dXVzJw5k7Vr1zKo5Td2AAAgAElEQVRkyBDbdmeywAA7d+7k66+/5syZM4wYMYKLL65fQd29e3c+/fRTwsPDOXDgAFdffbXtRt6URK49V111FW+++SY9e/YkODiY3r17c+yYOfvYlQS1FXdlqu1xJnFcWFhI7969+fDDDwGzZk5MTAxPPPGETXK6IQcPHuTtt9/m+eefZ/To0bzxxhusX7+eFStW8PDDD7N8+XKGDBnCunXrCAkJYfXq1dx77728++67jWSqnclZ24vtCUJL8G0DYCE7NoXs2BQGr1/CbaEmKoC6yFBCCk9AX/OUhmH0KPI3QlWt5vsKs0GIO3WMULAZAHusWkHRJYWkrFvX7msGeIPQ0FDOO+88li5dylNPPWXb7kwWGGDmzJlEREQQERFBTk4OmzdvJi0tzbbfaDRy2223kZubS3BwMPv377fta0oi154LL7yQ++67jx49etike624kqB2p429TLU9ziSOx48fz1133cWf/vQnpk+f7lLfxkq/fv3qyRVfcMEFNinjvLw8wGxIrr/+eg4cOIBSyibm5qhfjuSsRUpEaC1+YQAcUTIsnjPb8hmM2QCYMrL4NmkM60/VkVehSYpUPLb+H06Pt3oDJXh/GqipJ3VvERQUxFtvvcWkSZN4+OGHudfSD2eywOBYKtqeJ598kh49evDNN99QV1dHePhZr8qRRO6mTZtsksILFiyw6f6HhYUxcuRIHn/8cXbt2lVPVdKVBLU7bZw9OTuTOAaz17Bq1SruuecepkyZwv333+/gDGdpKFdsL2VsjT3cd9995OTk8N5775GXl0d2drbTfjmTzBaE1uCzMQBDmIH8snyn+xOCo4iozLNlAgFM7h7Mg8NCSYo0y0N8X6EpNrZJuUKHJTIykpUrV/L666+zdOlSwLksMMD7779PVVUVRUVFrF27ltGjR9fbX1paSq9evQgKCuK1116jttZ17UVmZia5ubnk5uY2Wjzkj3/8I48++mg9DXZwT4K6OTLVVpxJHB87dozIyEiuu+467rrrLrZvNycbtFZO2v5ztl+NrOF5nclZC0Jr8VkDML7PeDqFdHLewNCDuHBsmUD2jIsLIilSUVWric8/RMjmDV7sacenW7dufPTRRzz00EO8//77Nlng8ePHN5rfzsjI4OKLL2bMmDHcd9999O7du97+W2+9lVdffZUxY8awf//+Vs1Tn3vuuVx//fWNti9evJitW7eSmprKsGHDHMpPu9OmIVOmTOGaa65h7NixpKSkcPnll1NWVsa3335LRkYGaWlp/PWvf+XPf/4zYF5a8aKLLrIFgZvL//3f/3HPPfeQlZVVz1Dm5OSwe/duWxD4vvvuw2g0kpqaSnJycr0UWUFoDU3KQXuL1spBg1kSOv74bogwF3ANXr+EuaFlVEXHc+e5v6Zi9Xpqel9Ol8gwh8f/b/k6rtj0DrX9fkbl7D84bFNQVkW3siKG9Y336PrBvigHPX/+fKKjox0uci4IgufpCHLQPkNxgnnh8NBq9xQBvx40loM9Brhs090QTnkXkYcQBMH/8KsgcGHSWGqLNzfrmAqTptgIzsOJZ4vEksNdNAoA7IPBgiD4Pr7vAViKwVrCuDjz8KNPHiXi+cVOYwHekopor+k3QRA6Pm1xf/B9A5AwCkIcz/E3xeTuwez7WTpHu/Uh6PhRQnO3uWzvSXmI8PBwioqKxAgIgtAIrTVFRUX10qi9gX9MAVm9AEOPepvrIkOJLimkhHhqauvqVfpa+XrQWL4eNNZlTQCcXTksLb4rJg+sZpaQkEB+fj6eWhlNEAT/Ijw83FY46S38wwBYJCGs7AjRbNv/Puckp5NeZg7gbjzo/Ea7u8xcD9At1PXKYZ6UhwgNDa1XYSsIgtDW+P4UUAPOi04CYGN5HtTWwIlvXba3xgFKjebVwyKeX+wyHiAIguAv+J0BGDlo5tkVwqLiwVjhsv3k7sEMMyhW9RrBqThzUZM78QBBEARfx+8MgDMiO4XUk4WwZ1xcEP9LGsOvRs3h/8bN4VBMb4cSEQ3XDBAEQfBlAsYApCV2dSgLAWYvYHY/sycAZsXQUgfCjIG6cpggCP6JfwSBXRAcFY6psBBdXQ041w6a3D3Ytl7wiU8UvYqPcnzRU5QlpzNw0ln535raOnYfL2OstzsuCILgZXzaA6inCOqkICx6xCBipl+M6uRCOK4BZcnpHO/Wh17FRzF8t73evu6GcCprXCtcCoIg+AI+bQDqKYI2KAgLNlbS64dNfFl6sNnnHThpPD1vn8vxbn0c7lfGGkpXfhgwawYLguCf+LQBaIR1jWDLIvEDTu6nrOYnOPAp4DoQ3BxM8T0JiReBOEEQfBv/MgAWL8AY3pnaUIu8W1Q8VJUS3NnAueqM00CwM1wtJC8IguDL+JcBcEH0hAmoTp2I7BRCSUWNW55ATCikFh3GtLF+UZh1veCNBwvZfbzUW10WBEHwKgFjAKykJXblvJ/Fu+UJGEab11wYfHB7PS+guyGcLpFhdIkMk4CwIAg+i9+ngdoT3NlgE3KLriilrKQWU1UVXfsnOWxvysiicMtWqNSsP1VnSxO1xxoQdnQtT60eJgiC4A183gBYU0ETDE2r5tnfkK15/J8vec3lMd1ClcOiMCvWgHCj7aLyKQhCB8fnp4AaLQ4fFgW1Lu7YgiAIAuAHBqARCaNAOZd1bgl9io9y/aolHFj9pUfPKwiC0J74/BRQa4kIC6agrMrhYjEAxrSRlBs1/U8c5cR3CuxkIeBsRpA9kZ1CAn79YEEQOj5ueQBKqQuVUvuUUgeVUvOctLlSKbVbKbVLKfWGZ7vZMjoXHSLxyFaXbVIGJxBZXOB0vykji4jfzeVErOOqYPuMIOufCMYJguALNOkBKKWCgSXAZCAf2KKUWqG13m3XZiBwD5CltT6tlOrurQ67y44QzduGKLKP7WRVbF8Mm48xPuMPjdpFT5iA/vYHj167praO3T8UMaxBdpBkBgmC0JFwxwPIAA5qrQ9rrWuAN4GZDdr8BliitT4NoLV2/kjtBeqJwgGZhv4AfNAljk5BocR3HUBZZXGrr+NuVXB3QzjlXeJtchEiGyEIQkfEHQPQB/jR7n2+ZZs9g4BBSqkNSqmvlVIXOjqRUmq2UmqrUmqrJxdDb5gJlB0/gsHh7jsh1jiAK2JCoX/pMVLeeMat5SLtq4WlalgQhI6IOwbAUUpNw+WyQoCBQDZwNfCiUqpLo4O0fl5rPUprPSreQe68x2hmJtCwXjFNVgYbRo/iRGwfehYdpWyL67gCOI4NSNWwIAgdCXcMQD6QaPc+ATjmoM37Wmuj1vp7YB9mg9BxqK2xqYI2JLizgeiKUk4fznN6uCkji2+vuY3DMb1dFoYJgiD4Cu4YgC3AQKVUP6VUGHAVsKJBm+VADoBSKg7zlNBhT3a01VhUQR0RPWECY39zDSHhrnM3J3cPJjLEszUGgiAI7UWTBkBrbQJuAz4G9gBvaa13KaUWKKVmWJp9DBQppXYDa4C7tdZF3up0S6iuq2lycZiIsGBKKmrcVgttLsa6Otbsa9P4uCAIglPcKgTTWq8CVjXYdr/daw3cafnrUHQuOkR83leQNJbCn3502XZYrxibrk/D4i5PEBsVzrFKmT8SBKFj4DeVwIYwA4WV5pt2takawLYyWLf87RQmjT0bBxg42eE57NVCQwpPQN9zPNpHVV1Npy9X8/mGxo5XRFgwKYMTpE5AEIQ2w28MwPg+ZyUaVn1vdlaM4Z35KXbA2UYu4gBQXy00/IU3KKmooaa2rpFMRFWt5oHd5if5cXFBDmWiHVHXsxfdnOwrKKvi2335jJX7vyAIbYT/icF5iGG9YhwuHBMTCuHB5kBwXoV5nQBP0N0QLmmigiC0KX7jATRCBTeSha6uq2FV4XbYeAh6pQDmqSN776EpuoUquoXCg8NCbV6Ap7AGiXMGt7uShiAIAYD/GoBOBqgqqbcpITwOwuOg8jREmIO91rhBR0CCxIIgtCX+awA8RGSnkEZy0UHHjxLx/GLmVGh2JI2AYTJxLwiC7yEGwAnWjKBzVTVbgs/qDBnTRhJqed2n+KjllRgAQRB8DzEATrBmBJWu/JBIfdYLMGVkYcrIAuDEoqfqZQRZaU5mUEOqTLWsyD1ab5shIlTiAoIgeBwxAG6QltjVYWFYTCg01MrLq9Bwqq7FBiCxa1SjbYVeqEoWBEEIiDRQWzWwA6pN1Xx5tGVr/XYLVfSLVDw4LNT2lxSp2F3m3roB7lJlqhUJCUEQPI5fewD7Kk/yRvd+/Lbo0Nlq4AYkGBI8mgk0Li6I3WW1rG+BF6Cqqwn9YnWj7f2BnyprKD3grIzMjKw4JghCc/BbA5DZK5N9p/exOriGa+yrgQFMNZC/1bxuQBMEdzZgzM8nslOUw8XjrRlBVmYAKRWa1b1H8AD1DU5TsYG6nr2c7qsMqeK7qhDSErs6bWPy4CI7giD4P345BWQIM/CzLj9jcEQPJw16QM0Zt84VPWECqlMn0hK7NqoKNqaNpK5X48Xizyk5yqRjO+pta23VcHdDuCw2LwiCR/FLD2B8n/FmPSC7amBrHMA2DdQML8AZ9hlB9kQ8v5h+mKuFrXiiarimto7cH0+79AIEQRDcxS89ABudDKAUxQnpgFkV1EYzvICOgngBgiB4Ev82ABYKk8bWVwW1YvUC3MRaFdxSPJEdZPUCBEEQWktAGACnWLyAalM1q75f1WQ6qKM4gLuMizMf11r1UPECBEHwFH4ZA3BGozgAgKmGhNITkDDKaTqoVRZCV1c7zQZqisndg1l/qo68isaVw9C86mGJBQiC4AkCxgNwGAcAt2IB0RMmEDP9YqfZQO4yLi6IpMjGi8o3N0NIvABBEDxBwHgAhUljG9/8vUjD+gCAi9JGMtlB1tADu422+EBLJSQEQRCaS8B4AG2Jo/qAoONHCc3d5rC9p+IDgiAIzcH/PQAVDGUnzVM9OIkDeBhH9QENvQF7rPGB5ngBEgcQBKG1+L0HcKTmNI+e+Jy1Rd86jwNYsGYDuZMR5Gma6wVIHEAQhNbi1wYgs1cm5xjO4Uh1MZtK9rquBzj8BQmlJ4iPiCc+Ip6ymrJGzazZQNElhZRU1FBSUdOqugB7JncPZpihcYDYFTW1dWw8WCh1AYIgtAi/ngLKTswmOzGbR9fd47qhZXqIStc3UqvS5rCVHxISb15T2NE6AW2FNRW1pKKm3fogCILv4tcGwBdxVCfQmhXGBEEQnOHXU0C+hqM6gdaqiAqCIDhDPAB7PKAQ2homdw9u9KTvTo2ANRYQXFKM2lcg6wcLguAWAeMB7DtzlLVF37puZFcV7K4+UGSnEFtA2NOBYXAvO6i7IZwukWF0jgijrLL1stOCIAQGAeEBZBr6s6/yJJtK9nIL7tUCJBgSAJpcLtJRHr4nA8PNrRGoMtWyIvcoAIaIUPEGBEFwilsegFLqQqXUPqXUQaXUPBftLldKaaVU+8yh2GEIM5Bflg9AdvwIBoebb4RN1QJ4gtbKRjekOTUCiV2jiDeEE28IF29AEASXNGkAlFLBwBLgImAYcLVSapiDdgbgD8AmT3eyJYzvM55OIZ3MbxJGgTIHV621AFYvoBHNXCPAEa0RjHOEtUaguesJWL0B69+afQUe65MgCL6PO1NAGcBBrfVhAKXUm8BMYHeDdn8BHgPu8mgPvUBxQjqdiw7RLX9742kgQ48m6wHcwRobsBJaV0fIiWMuJSGMaSMdLjEJZi9gd1kt60/VuZ0Smtg1qt77Qg96JYIg+D7uPKb2AX60e59v2WZDKTUCSNRar/Rg37xGk16AB0hL7Mp5P4u3/fXKmYCpZ2+n7V2JxUHLvQBBEARnuOMBONIn0LadSgUBTwI3NHkipWYDswHOOecc93roKVQwRyoLePTQ2wCEdg0nPLQ755/cyhAvCsNZic7JoTAxmS6RYQ73u/IMrLTEC7BHAsSCINjjjgeQDyTavU8Ajtm9NwDJwFqlVB4wBljhKBCstX5eaz1Kaz0q3iKl0FZkJk3inE7dbO+N4Z3ZGx7O55xx6QVUm6obpYJaNYFMhYUY8/O91ueGtNYLkACxIAj2uOMBbAEGKqX6AUeBq4BrrDu11qVAnPW9UmotcJfWunWRVA+TnZhNtlFBxNm0zSd2LQWqGscC7ArCEgwJjVJBrZpAAKUrP3S7D9bsoOYuJ2lPa70AK+INCILQpAegtTYBtwEfA3uAt7TWu5RSC5RSM7zdQW9iDO9MbWhE4x1uLBPZEjyRHdQS1VBHiDcgCIJbhWBa61XAqgbb7nfSNrv13fISYVH1FodpDxpmB1kJrq3zaOpoc7B6A+IJCEJgERCVwDYSRsHhL9q1C85W8DrQzId6T64hbE0XlTRRQQgsAkYLyBk7QjTvhzqosLUsEtPaojB3CVKKmlr3VD+9tYZwlalWisUEIYAIaAOQ2WUIAJ+pSgavX1I/G8jQwxwwrjnjUhjOPiPI/q+52UHR4aEEuekFWOMA1rUDPFUXkNg1SuIBghBABNYUENSLA2THprD9xGbqQsqJLDRnxDSqDDbVkFB6AhJGORSGs88Isqc52UFWgpRyO0toXFwQnKojr0JDKzOC7JF4gCAEDoHnASSMgpCzxVjG8M5UxPSmIqaP4/ZeyghyRHR4qNuB4Mndg3lwWGijBWRaizU7SDwBQfB/As8ACG5h9QQkJiAI/kvgTQG5wJ11AryNozTRmto6l9NCtnWETRGct+ckFwxtfZqrZAYJgv8jBgA4UlnILV3DuawmiqmOFEKtlcGxfb3el+YuMGONBQDkVQehD53yiAGwYl8xDFI1LAj+hF8bAOuiMNbVvRxhzQQ6VFnIB13imFrh6ESekYhuKa4kJOzXEX5gZyV5RRUsWLkLgKwBca02BiIpLQj+i1/HAOotCuOE7NgU/jTgCs6JaFqczt11gj2NuxIS4w1G+sZGAvBDUQUbDp3yeF9kkRlB8B/82gNwigtJCFdxAHfXCfYG7gjJTe5iInviuQAsWLmLH+y8ARCPQBCE+vi1B+CUBqmgVoydogEn6wV7YKnI1tBcIbmsAXE2bwC86xGIFyAIvklgegDg0Aswhnfmp9gBjts3Mw5grRC2R1dXE5rgPB7RFM2Rk75gaI96T/sNPQJPeANg9gjECxAE3yRwDYADYbgjlYXcFmkCIKXoW7JjU1p8ekcVwi2pDrYnLbGry4wgV2QNiAPMHsAPRRWAZ7OFBEHwPQLXAEA9L8CaDUTNMQ6qWownNjc2AFaBOGqh37Q27y60fFEZe4/A096AdRpI0kMFwbcIzBiAFbtYgDUbaEG3DAbX1BBaXd64vVUgzljZxh09iycWlbGPD3giNiAicoLgmwS2BwCNYgGFSWOpLd7czp1yTWuXlnTmDbTGE2hYMOYOUlQmCO2LGIAOsEhMc2lNLKAh1thAa+MCDdND3UGCx4LQvvi9ATCEGSisLKTaVO2yIrghQaZqBq9fQnFCeuOagNoaOPApDJzs4d62PVZvYMHKXew5XsZnHtIScoemvAbxEATBu/h9DGB8n/FM6zetyYpge4ydotnXKYw/hJ3hy5MOcv+j4qGq1IO9bD7WaSBPYfYE8EqtgDPsF6Z39CdxBUHwLn7vAbSE9J4ZGEv2so+jqKozjHdUGWysgm/fgfAYj3kCNUeOcPJvf3OrbS8gqrKGkKCzNryuWyzGiZNadO0LhvZgw6FTHokHeAoRohME7+L3HoBbhEWZi7zKTgJnM4KSQgyAk8rgrn0huofHPIGoMWMIO+ecZh0TpBSmujq31xJuCmt2kLeqhptLQw9BPAJB8CziAYA5EAyNgsHG8M7mdNA2uO9E5+QQnZPTrGOsz+cbDxbSJTIMdbq4VX2wjwd0JE/Air1HIN6AILSegDEA7khDO6OhQFx1XQ2rCsyxAYOxmvEe7Wn7Y80M2nO8jD3Hy+p5A+1pEOwzjSSDSBBaT8BMAbkjDe0IRwJxCeFxxId1Jj6sM2W1/ncjumBoD+6ffi43j+vH0F4G2/aOMjUkCIJnCBgPwC0cCMR9V1vGaz37MdPZNFAHSAm1LiMZXFnDqdNnWpST74imBOWg/TwCCRALQusRA2BPg6KwzC5D2HfmKJ+G1jHT6MRZ6gApodZlJE2FsN0QbpseqTLVeswYQH1BOWhfUTlZl0AQWo8YABdkx6awqWQvlB5zvWC8sardvQAr9k/BzZVmaApHHkFbF485QwTpBKH5BEwMAM5WBdv/5Zfl129knQayw+VCMWBOCW1nL8ARhohQCsuqbH8/nj7j0fO3R/GYM0SQThCaT0B5AOP7NM7XWfX9qvobHGgDWReKiSw96lwewg0vwNEiMfa0dsGYhjR8GvaGR7Dh0KkO5QVImqgguE9AGYDWUJyQDkBkqfkG08gAdO0L5ScbHlYPR4vE2NPaBWOawuoRgOfiA1kD4mypou1tACRNVBCah1sGQCl1IfAUEAy8qLV+pMH+O4GbARNQCNyktf7Bw31tN45UFnJXBGQOncAte9bV22dfE8CZQth4CHqlYAgzOPQ42hNvxAcaSkg0B29mEIk3IAhN06QBUEoFA0uAyUA+sEUptUJrvduu2Q5glNa6Qik1B3gMmOWNDrcJDlYKO1Jpnrq5pUHThPA4u+M6myUlIuIprPSMXLO38KQ30DA7yB28nUEk3oAgNI07HkAGcFBrfRhAKfUmMBOwGQCt9Rq79l8D13myk97EYYWwXRwgOzaF7NgUHj30tm23y4wgH8H+iXjNvgIKy6pabAgaZge5g6OagoZ0JBkKQfBH3DEAfYAf7d7nA5ku2v8a+J+jHUqp2cBsgHOaKXzmLcb3Gd84EOyC4oR0Ohcdolv+dp82APZYjYGng8SuaMpr8KSH4GzdAZkaEgIddwyAcrBNO2yo1HXAKGCio/1a6+eB5wFGjRrl8BwdBqtCqKmmXmVwYdJYuuVvd50R5KMYIkL50YOVxK5oymvwZNWxs/HI1JAQ6LhjAPKBRLv3CcCxho2UUpOA/wdM1FpXe6Z77YgThVBoIiPIVGM5phb6TfN2Lz1KzuDubeoFuKItqo5FTkIIdNwxAFuAgUqpfsBR4CrgGvsGSqkRwD+BC7XWBR7vZXtiVxh2pLKQRw+9TWaXIWSP+x2D1y9p3N7qLfz0Y+N9PoA1OOxpGYnm0pQOkSfiAyInIQQ6TVYCa61NwG3Ax8Ae4C2t9S6l1AKl1AxLs78D0cDbSqlcpdQKr/XYC1grhBtVBYPZEwgJI7PLEM6JiOdIZaFZHqIprCJxPkbO4O7MSOtDeEhwe3elHtbFasB7qqRWj2DNPv96hhEEZ7hVB6C1XgWsarDtfrvXLVuHsINgzdd3FQx2lA0EzjOCqiNiWHVsA4SclSfoiLUBzugonoAVe4/AW6qk1nGKJyAEClIJ7A4OZKLBdUZQQngcaHNNgJWOXhtgj3Uu3JoiaqUjGISOpEoqCL6MGAB3cKAPBGczgpzWBVgDwmFRZ4PKPoa39YRagrdVSR2ljUqAWPBHAkoN1BtYM4IcKoUaekBEV6jxrApne+JthdGW4GlV0oaL0cuC9IK/Ih6AHdZgMEC1qdqt9YOb9AL8jI7qETTUI/J0FbGzYrKWIN6E0FEQA2CHfYC2UUDYwToBVvyxOthd2rJ4zBX2cQFvxAQ8OT4JMgsdBZkCcpeEURDZFepMUFt/OqAwaSw/xQ6weQGBRM7g7sQbwtt9Ksi6kP3908+lb2ykzRtYsHIXn+1xLdMtCIGKGIDmkDAKwruAaqyO4TIWYA0GH//W2z1sF3IGd+9QdQNtUTPQGqTeQOgoyBSQE6zxgObGAhyfzLerg93BXl7aXbyVUtpUzYCV9lIblXoDoaMgBsAJ7hSHNZfqOqPL88We3ktyfLzT/R2ZlgQ12yKA7Ex1tCPUDrQ0sCxBZMFTiAFoAUeqT/Po/jcgONSsCxSbYtvnSiU0oVPXeoVhDakw7fBanzsiTXkNnvAQnKmOequauDm0dGziOQieQgxAM8nsZVkKobqMI5UFUGu0GYAm1w0W6tHUU6w3q5ClmlgQxAA0ScMVw7ITs8lOzAbg0c2PQlWJrW1h0lgKk8YyeP0Sx3UBNqloCz5cIdwWeLPmwJXaaEdficyTNQn2yNRS4CEGoAmau2IYuKgLaKAlROVpD/QwcPDkOsYNsXoEvuAJeKvmQqaWAg8xAK1FBTcSinN71TBTDeRvFS/ATRytY2yltQbB6hE4yxrq6F6BJ/CEZyFehG8hBqC1dDJASFijzW7FAww9xAtoId6aHnKUNeQLXoEn8IRnIV6EbyEGwA1c1QQcKTvCo9Vl5grh4FDbdperhgkex1OSFI6yhlzVErSEQPAmBN9ADIAbOKsJsGUEAVSXmSUigkM5YhGUs08PdYpMA3mEnMHd600LeTNjqDX4uzfhqQC1TCW1DWIAWoF9RhBgzvCJ6Npo1TCXsQCZBvIY9jcMb2YMtQZfyjZqCZ4yujKV1DaIAfAkDhRDpTagffBmxlBr8KVsI8H/EQPQDJrUB3Kwcph9bYBTrPUBJUc83OPAxVXGkJX2MAxNZRu5gz96Dg2pMtWyZl+BTAN5GTEAzcAb+kCALYW0LkxjKmz5usHBnQ2e6pFf4ewm0p6L2bQ0rhAonkNi1yiZBmoDxAB4mrAoqDNxpLq4XiwgMtLExT+VcpGdJ9AwJlA5NI6Y+Er3rxUeAwMne6TbgUhDLaK29AhaGlfwdEaSPR3Ns/BWxXNL8cfAtBgAT5MwikxdDnmrbVlBAPtDgqBzDBdVmJs5jAlExUN0M/4DlstCJ62h4X9m61RRR4oZNMSTGUn2dETPoqN9B/7okYgB8AK27CBLVhDAo4fepgLYl34FgHO9IKHdsBqEjvTU2RBPZiTZ4ynPoqN5EYJrxAB4k7Aoc4qnqabRLo+sI2ysgm/fOftepoQ8giOZ6o7sFXgCT3gWHdGL8CStnRYh8icAAAqwSURBVJLqiFNIYgC8ibW4q0FmEJzVC2qVF9C1b/33MiXkERz9J3WWSdQSOqIx8YRn4cn4REf0JFr7nXXEKSQxAC3Amg4KuL1kpCM84gUIbYInn9w68hRTa/BUfMLfPYmOhBiAFmBNBwU3U0KdZAYBRCb2Jch0hrrtT5Ee2Ztp3UUSwt+xTjF1RE+gNXgqPrFg5S72HC/jsz0nxQh4GTEAbYE1M+j7jxvtMnaKJhQ4EFRLXYV/PhkK9bF6Ey2dVvI3w9GQrAFx7DlexoZD4gV4GzEArcR+OsiehlND2YnZZKtoqDlztpGpxlYE9vT2pwgyVcNH89y/eP9sGHRhC3sutDctnVbyZDzCno5iWC4Y2oMNh055rd6hpXTEuERrccsAKKUuBJ4CgoEXtdaPNNjfCfgXMBIoAmZprfM829WOif10kD0Op4YaKn7mb7UtJmPsFE0wGhonDDmm+Hvzv2IAAg5vZZJ4wrB4yoh4q96hpfhrXKJJA6CUCgaWAJOBfGCLUmqF1nq3XbNfA6e11j9TSl0FPArM8kaH/YqEUTYjYAzvTEVYFIy+371jP5pnNgL2HkPvdEi53Dt9FfweTxgWT3kng3pGd6ibrScynFL7xDAjrY8He9V63PEAMoCDWuvDAEqpN4GZgL0BmAnMt7x+B3hGKaW01tqDffUpHE0NOcwYshqBOpP5z136Z9d/X/y9ufJYENoRT3kna/YVeGSBH0/RWo/kh6IKjKY6z3XIQ7hjAPoAP9q9zwcynbXRWpuUUqVALB3Jh2tjHE0NfXn0S4fxAmL7YsyPIqimwv1c/t4jIH7o2VqAj+aJARD8hoYL/LQ3qQkxDOoZ3WKDtGDlLp81AMrBtoZP9u60QSk1G5gNcM4557hxaf/CWbwA4NvCbzlx5gRb4s039NE9Rzd9wgOfnjUYhp4Q3MkT3RSEDkFHq5ptjUHqbuhEaHAQXx0qsm0bOyDWU11rMaqpWRql1FhgvtZ6quX9PQBa67/ZtfnY0uYrpVQIcAKIdzUFNGrUKL1161YPDEEQBCFwUEpt01p7pGAoyI02W4CBSql+Sqkw4CpgRYM2K4DrLa8vBz4P5Pl/QRAEX6DJKSDLnP5twMeY00Bf0lrvUkotALZqrVcAS4HXlFIHgWLMRkIQBEHowLhVB6C1XgWsarDtfrvXVcAVnu2aIAiC4E3cmQISBEEQ/BAxAIIgCAGKGABBEIQARQyAIAhCgCIGQBAEIUBpshDMaxdWqhD4oYWHdyypwLZHxi/jl/EHLoO11gZPnKjd1gPQWse39Fil1FZPVcL5IjJ+Gb+MP7DH76lzyRSQIAhCgCIGQBAEIUDxVQPwfHt3oJ2R8Qc2Mv7AxmPjb7cgsCAIgtC++KoHIAiCILQSnzMASqkLlVL7lFIHlVLzmj7CN1FK5SmlvlVK5Vqj/kqpbkqpT5VSByz/drVsV0qpxZbPZKdSKr19e998lFIvKaUKlFLf2W1r9niVUtdb2h9QSl3v6FodDSdjn6+UOmr5/nOVUtPs9t1jGfs+pdRUu+0++X9DKZWolFqjlNqjlNqllJpr2R4o37+z8Xv/N6C19pk/zHLUh4D+QBjwDTCsvfvlpbHmAXENtj0GzLO8ngc8ank9Dfgf5pXZxgCb2rv/LRjvBCAd+K6l4wW6AYct/3a1vO7a3mNr4djnA3c5aDvM8rvvBPSz/H8I9uX/G0AvIN3y2gDst4wzUL5/Z+P3+m/A1zwA2wL1WusawLpAfaAwE3jV8vpV4Od22/+lzXwNdFFK9WqPDrYUrfU6zGtJ2NPc8U4FPtVaF2utTwOfAhd6v/etw8nYnTETeFNrXa21/h44iPn/hc/+39BaH9dab7e8LgP2YF5nPFC+f2fjd4bHfgO+ZgAcLVDv6oPyZTTwiVJqm2UtZYAeWuvjYP7RANZFU/31c2nueP3tc7jNMsXxknX6Az8fu1IqCRgBbCIAv/8G4wcv/wZ8zQC4tfi8n5CltU4HLgJ+p5Sa4KJtIH0u4Hy8/vQ5/AMYAKQBx4HHLdv9duxKqWjgXeB2rfVPrpo62Obzn4GD8Xv9N+BrBiAfSLR7nwAca6e+eBWt9THLvwXAe5jdu5PWqR3LvwWW5v76uTR3vH7zOWitT2qta7XWdcALmL9/8NOxK6VCMd/8Xtda/9eyOWC+f0fjb4vfgK8ZAHcWqPd5lFJRSimD9TUwBfgO81itmQ3XA+9bXq8AfmXJjhgDlFpdZx+nueP9GJiilOpqcZenWLb5HA1iOJdi/v7BPParlFKdlFL9gIHAZnz4/4ZSSmFeV3yP1voJu10B8f07G3+b/AbaOwLegoj5NMxR8kPA/2vv/nhpjP0xR/C/AXZZxwnEAp8BByz/drNsV8ASy2fyLTCqvcfQgjH/B7Oba8T8JPPrlowXuAlzUOwgcGN7j6sVY3/NMradlv/Eveza/z/L2PcBF9lt98n/G8A4zFMVO4Fcy9+0APr+nY3f678BqQQWBEEIUHxtCkgQBEHwEGIABEEQAhQxAIIgCAGKGABBEIQARQyAIAhCgCIGQAholFIjlFIvttO1V9uV9wtCmyMGQAh07gWe9tbJlVIhLna/BtzqrWsLQlOIARB8AqXUcosw3i47cTyUUr9WSu1XSq1VSr2glHrGsj1eKfWuUmqL5S/LwTkNQKrW+hulVJBFQz7esi/Ioqke5+xcSqkMpdRGpdQOy7+DLdtvUEq9rZT6ALOgXy+l1DqLpvt3Sqnxli6sAK727icnCM5x9XQiCB2Jm7TWxUqpCGCLUupdzHro92HW0i8DPsdcPQ3wFPCk1nq9UuoczJIAQxuccxSW8nqtdZ1S6t/AtcAiYBLwjdb6lFLqDSfn2gtM0FqblFKTgIeBX1jOPRazcSlWSv0R+Fhr/VelVDAQabnmaUs5f6zWusizH5cgNI0YAMFX+INS6lLL60TM+ic9gS+01sUASqm3gUGWNpOAYWaZFQA6K6UM2qy3bqUXUGj3/iXMejOLMEsKvOzqXEAM8KpSaiDmUv5Qu3N9au0XZo2WlyyCX8u11rl27QqA3oAYAKHNEQMgdHiUUtmYb8JjtdYVSqm1QDiO5W+tBFnaV7poU2k5DwBa6x+VUieVUucDmZi9AafnUko9DazRWl9q0XFfa7f7jN1511nkvC8GXlNK/V1r/S/L7nBLPwShzZEYgOALxACnLTf/IZiXAQSzAuJEi/pjCGenXwA+AW6zvlFKpTk47x7gZw22vQj8G3hLa13bxLligKOW1zc467xSqi9QoLV+AbPqY7plu8LsxeQ5O1YQvIkYAMEX+AgIUUrt5P+3d8cqDQRRGIXP7QM+jGmFPI6VYKedYBFIl8pGfIVIai0DsTJFCHZWPoBoIwjXYkYEMZGQiIQ5X7cLd3bZYn/u7DID58AUIDOfKPPud8ANsACea80R0I2ym9ICOPw+aGY+AHufS29XY6DD1/TPqrEGQD8iJpT9WJfpAbOIuKeE1LCe3wemmfn++yOQts/VQLXTIqKTma+1AxgBV5k5WqP+GHjJzMt63KV88D1YXbm5iBgC48y8/etrST+xA9CuO4uIGeVvnkfges36C+ANICJOKLsynW71Dpeb+/LXf7IDkKRG2QFIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRn0Afm67ysLkOcsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "km.fit(true_ages_at_time_of_death, np.ones(number_of_vampires), label='true lifetimes')\n",
    "ax = km.plot()\n",
    "\n",
    "km.fit(ages_with_censoring, np.ones(number_of_vampires), label='observed lifetimes (some censored)')\n",
    "ax2 = km.plot(ax=ax)\n",
    "\n",
    "km.fit(ages_with_censoring_observed, np.ones_like(ages_with_censoring_observed), label='observed lifetimes of deceased vampires only')\n",
    "ax3 = km.plot(ax=ax2)\n",
    "\n",
    "km.fit(ages_with_censoring, observed, label='Kaplan-Meier estimate')\n",
    "plot = km.plot(ax=ax3)\n",
    "\n",
    "plot.set_xlabel('age (years)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}