{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3Q2KuaRamHfr"
   },
   "source": [
    "# Survival Estimates that Vary with Time\n",
    "\n",
    "Welcome to the third assignment of Course 2. In this assignment, we'll use Python to build some of the statistical models we learned this past week to analyze surivival estimates for a dataset of lymphoma patients. We'll also evaluate these models and interpret their outputs. Along the way, you will be learning about the following: \n",
    "\n",
    "- Censored Data\n",
    "- Kaplan-Meier Estimates\n",
    "- Subgroup Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outline\n",
    "\n",
    "- [1. Import Packages](#1)\n",
    "- [2. Load the Dataset](#2)\n",
    "- [3. Censored Data](#)\n",
    "    - [Exercise 1](#Ex-1)\n",
    "- [4. Survival Estimates](#4)\n",
    "    - [Exercise 2](#Ex-2)\n",
    "    - [Exercise 3](#Ex-3)\n",
    "- [5. Subgroup Analysis](#5)\n",
    "    - [5.1 Bonus: Log Rank Test](#5-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "UopnLTeLkViX"
   },
   "source": [
    "<a name='1'></a>\n",
    "## 1. Import Packages\n",
    "\n",
    "We'll first import all the packages that we need for this assignment. \n",
    "\n",
    "- `lifelines` is an open-source library for data analysis.\n",
    "- `numpy` is the fundamental package for scientific computing in python.\n",
    "- `pandas` is what we'll use to manipulate our data.\n",
    "- `matplotlib` is a plotting library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "TZyXoADQmYlt"
   },
   "outputs": [],
   "source": [
    "import lifelines\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from util import load_data\n",
    "\n",
    "from lifelines import KaplanMeierFitter as KM\n",
    "from lifelines.statistics import logrank_test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "5rp2TD1qnGmp"
   },
   "source": [
    "<a name='2'></a>\n",
    "## 2. Load the Dataset\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "WEbu3MtrVsnU"
   },
   "source": [
    "Run the next cell to load the lymphoma data set. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "e3wHdLrEnSNa"
   },
   "outputs": [],
   "source": [
    "data = load_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3hrHa0dPqU08"
   },
   "source": [
    "As always, you first look over your data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 221
    },
    "colab_type": "code",
    "id": "QEd504pKqWuc",
    "outputId": "7297830a-d316-4623-bb6a-77f8f96b8805"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data shape: (80, 3)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Stage_group</th>\n",
       "      <th>Time</th>\n",
       "      <th>Event</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>6</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>19</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>32</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>42</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>42</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Stage_group  Time  Event\n",
       "0            1     6      1\n",
       "1            1    19      1\n",
       "2            1    32      1\n",
       "3            1    42      1\n",
       "4            1    42      1"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(\"data shape: {}\".format(data.shape))\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "dblUOLQS4UU0"
   },
   "source": [
    "The column `Time` states how long the patient lived before they died or were censored.\n",
    "\n",
    "The column `Event` says whether a death was observed or not. `Event` is 1 if the event is observed (i.e. the patient died) and 0 if data was censored.\n",
    "\n",
    "Censorship here means that the observation has ended without any observed event.\n",
    "For example, let a patient be in a hospital for 100 days at most. If a patient dies after only 44 days, their event will be recorded as `Time = 44` and `Event = 1`. If a patient walks out after 100 days and dies 3 days later (103 days total), this event is not observed in our process and the corresponding row has `Time = 100` and `Event = 0`. If a patient survives for 25 years after being admitted, their data for are still `Time = 100` and `Event = 0`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "L0d2s2wtn2Pf"
   },
   "source": [
    "<a name='3'></a>\n",
    "## 3. Censored Data\n",
    "\n",
    "We can plot a histogram of the survival times to see in general how long cases survived before censorship or events."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 285
    },
    "colab_type": "code",
    "id": "cFrvXrODZklx",
    "outputId": "a260b523-e792-47dd-9833-7677cd4b5741"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEJCAYAAACT/UyFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de5gcZZ328e9NQAmJcjDsbDgZ4EV8kSCQUVFAExRXYQV8F0UWENQ1uoqg4q7ouhJRV1wFPCAiCnJwJYh4iKIvAmZA8QAJh0wgIogRROSkIuFo4Ld/PE+TYtLdUzPp6ulK35/r6muqqqur7q6p+U3101VPKSIwM7P+sc5EBzAzs+5y4Tcz6zMu/GZmfcaF38ysz7jwm5n1GRd+M7M+U1nhl7SlpIWSbpR0g6Sj8/R5ku6QdF1+7FNVBjMzW52qOo9f0nRgekRcI+kZwGLgAOD1wIqI+HQlKzYzs7bWrWrBEXEncGcefkDSMmDz8Sxr2rRpMWPGjNLzP/jgg0yZMmU8q5oQdcsL9ctct7xQv8x1ywv1yzzWvIsXL743IjZd7YmIqPwBzABuA54JzAOWA0uAM4GNR3v9rFmzYiwWLlw4pvknWt3yRtQvc93yRtQvc93yRtQv81jzAouiSU2trKmnQdJU4HLg4xHxLUkDwL1AAB8lNQe9ucnr5gJzAQYGBmbNnz+/9DpXrFjB1KlTOxG/K+qWF+qXuW55oX6Z65YX6pd5rHnnzJmzOCIGV3ui2X+DTj2A9YCLgfe2eH4GsHS05fiIv/fULXPd8kbUL3Pd8kbUL3OnjvirPKtHwBnAsog4qTB9emG21wJLq8pgZmarq+zLXWB34DBgWNJ1edoHgYMl7Uxq6lkOvK3CDGZmNkKVZ/X8FFCTp35Q1TrNzGx0vnLXzKzPuPCbmfUZF34zsz7jwm9m1meqPKunJ8w49qIJW/fyE/adsHWbmbXiI34zsz7jwm9m1mdc+M3M+owLv5lZn3HhNzPrMy78ZmZ9xoXfzKzPuPCbmfUZF34zsz7jwm9m1mdc+M3M+owLv5lZn3HhNzPrMy78ZmZ9xoXfzKzPuPCbmfUZF34zsz7jwm9m1mdc+M3M+kype+5K2hjYDHgYWB4RT1SayszMKtOy8EvaEHgncDDwNOAeYH1gQNIvgFMjYmFXUpqZWce0O+L/JnAOsGdE/KX4hKRZwGGStomIM6oMaGZmndWy8EfE3m2eWwwsriSRmZlVatQvdyXtLmlKHj5U0kmSnl19NDMzq0KZs3q+CDwk6fnAMcBvSE1AZmZWQ2UK/8qICGB/4JSI+ALwjGpjmZlZVcqczvmApA8AhwIvlbQOsF61sczMrCpljvgPAh4F3hIRfwS2AD5VaSozM6tMmSP+90TE+xsjEXGbpOdVmMnMzCpU5oi/2Wmdrx7tRZK2lLRQ0o2SbpB0dJ6+iaRLJN2cf2481tBmZjZ+LQu/pH+VNAxsL2lJ4fFbYLjEslcCx0TEDsBuwDsl7QAcC1wWEdsBl+VxMzPrknZNPV8Hfgh8gqcW5wci4k+jLTgi7gTuzMMPSFoGbE46O2h2nu1sYAh4f5NFmJlZBdpduXs/cD9wsKRJwECef6qkqRFxW9mVSJoB7AL8EhjI/xQA/piXa2ZmXaJ0in6bGaQjgXnAXUCjV86IiJ1KrUCaClwOfDwiviXpLxGxUeH5P0fEau38kuYCcwEGBgZmzZ8/v8zqAFixYgVTp04FYPiO+0u/rtNmbr5hqfmKeeuibpnrlhfql7lueaF+mcead86cOYsjYnDk9DKF/xbgRRFx31hDSloP+D5wcUSclKfdBMyOiDslTQeGImL7dssZHByMRYsWlV7v0NAQs2fPBmDGsReNNXbHLD9h31LzFfPWRd0y1y0v1C9z3fJC/TKPNa+kpoW/zFk9t5OafMZEkoAzgGWNop8tAA7Pw4cD3x3rss3MbPzKnMd/KzAk6SLShVwAjCjmzewOHAYMS7ouT/sgcALwDUlvAX4HvH7Mqc3MbNzKFP7b8uNp+VFKRPwUUIunX152OWZm1lmjFv6I+AiApA0i4qHqI5mZWZXK9Mf/Ykk3Ar/K48+XdGrlyczMrBJlvtz9DPAPwH0AEXE98NIqQ5mZWXXKFH4i4vYRkx6vIIuZmXVBmS93b5f0EiDyeflHA8uqjWVmZlUpc8T/duCdpH527gB2zuNmZlZDZc7quRc4pAtZzMysC1oWfkn/HhH/LenzwGr9OkTEUZUmMzOzSrQ74m+045fvJMfMzHpeu26Zv5cHH4qIC4rPSXpdpanMzKwyZb7c/UDJaWZmVgPt2vhfDewDbC7pc4Wnnkm6raKZmdVQuzb+P5Da9/cDFhemPwC8p8pQZmZWnXZt/NcD10v6ekT8rYuZzMysQmWu3J0h6RPADsD6jYkRsU1lqczMrDJlvtz9KvBFUrv+HOAc4GtVhjIzs+qUKfyTI+Iy0v15fxcR84ByN5M1M7OeU6ap51FJ6wA3SzqS1F9PfW5Lb2ZmT1HmiP9oYAPgKGAW6T66h7d9hZmZ9awynbRdDZCP+o+KiAcqT2VmZpUpc+vFQUnDwBJgWNL1kmZVH83MzKpQpo3/TOAdEfETAEl7kM702anKYGZmVo0ybfyPN4o+QET8FHfZYGZWW2WO+C+X9CXgPFK//AcBQ5J2BYiIayrMZ2ZmHVam8D8//zxuxPRdSP8I9upoIjMzq1SZs3rmdCOImZl1R5k2fjMzW4u48JuZ9ZmWhb9xe0VJW3cvjpmZVa3dEX/j9ooXdiOImZl1R7svd++T9CNga0kLRj4ZEftVF8vMzKrSrvDvC+wKnAuc2J04ZmZWtXa3XnwM+IWkl0TEPZKm5ukrupbOzMw6rsxZPQOSrgVuAG6UtFjSjhXnMjOzipQp/KcD742IZ0fEVsAxeVpbks6UdLekpYVp8yTdIem6/Nhn/NHNzGw8yhT+KRGxsDESEUPAlBKvOwt4VZPpJ0fEzvnxg1IpzcysY8oU/lsl/aekGfnxIeDW0V4UEVcAf1rjhGZm1lFlCv+bgU2Bb5HO6Z+Wp43XkZKW5KagjddgOWZmNg6KiOoWLs0Avh8RO+bxAeBeUq+eHwWmR0TTfyKS5gJzAQYGBmbNnz+/9HpXrFjB1KnpfvDDd9w//jfQJQOT4a6HO7e8mZtv2LmFtVDcxnVQt7xQv8x1ywv1yzzWvHPmzFkcEYMjp3e18Jd9bqTBwcFYtGhR6fUODQ0xe/ZsAGYce1Hp102UY2au5MThMj1kl7P8hH07tqxWitu4DuqWF+qXuW55oX6Zx5pXUtPC39VO2iRNL4y+Fljaal4zM6tG28NMSZOAoyLi5LEuWNJ5wGxgmqTfk27kMlvSzqSmnuXA28a6XDMzWzNtC39EPC7pYGDMhT8iDm4y+YyxLsfMzDqrTMPylZJOAc4HHmxM9L12zczqqUzh3zn/PL4wzffaNTOrKd9z18ysz4x6Vo+kAUlnSPphHt9B0luqj2ZmZlUoczrnWcDFwGZ5/NfAu6sKZGZm1SpT+KdFxDeAJwAiYiXweKWpzMysMmUK/4OSnkX6QhdJuwG93w+CmZk1VeasnvcCC4BtJV1J6rDtwEpTmZlZZcqc1XONpJcB2wMCboqIv1WezMzMKjFq4Ze0PvAOYA9Sc89PJJ0WEY9UHc7MzDqvTFPPOcADwOfz+D8D5wKvqyqUmZlVp0zh3zEidiiML5R0Y1WBzMysWmXO6rkmn8kDgKQXAeU7xzczs57S8ohf0jCpTX894GeSbstPbQX8qgvZzMysAu2aev6xaynMzKxrWhb+iPhdYzjfFH3LEfP/brUXmZlZzytzOudHgSOA35Cv3sXdMpuZ1VaZs3peD2wbEY9VHcbMzKpX5qyepcBGVQcxM7PuKHPE/wngWklLgUcbEyNiv8pSmZlZZcoU/rOBTwLD5K6ZzcysvsoU/oci4nOVJzEzs64oU/h/IukTpK6Zi00911SWyszMKlOm8O+Sf+5WmObTOc3MaqpMf/xzuhHEzMy6o8wFXB9uNj0iju98HDMzq1qZpp4HC8Prk/rwWVZNHDMzq1qZpp4Ti+OSPg1cXFkiMzOrVJkrd0faANii00HMzKw7yrTxN/rlB5gEbAq4fd/MrKbKtPEX++VfCdwVESsrymNmZhUbtakn98v/e+BvpCP+zSRtVXUwMzOrRpmmnncBxwF3saqvngB2qjCXmZlVpExTz9HA9hFxX9VhzMysemXO6rkduL/qIGZm1h1ljvhvBYYkXcRTO2k7qd2LJJ1J+mL47ojYMU/bBDgfmAEsB14fEX8eV3IzMxuXMkf8twGXAE8DnlF4jOYs4FUjph0LXBYR2wGX5XEzM+uiMlfufmQ8C46IKyTNGDF5f2B2Hj4bGALeP57lm5nZ+LQ84pf0ZUkzWzw3RdKbJR0yxvUNRMSdefiPwMAYX29mZmtIEdH8CWln4IPATNIN1+8hddK2HfBM4EzgtIh4tOkC0jJmAN8vtPH/JSI2Kjz/54jYuMVr5wJzAQYGBmbNnz+/9JtasWIFU6dOBWD4jt7/XnpgMtz1cOeWN3PzDTu3sBaK27gO6pYX6pe5bnmhfpnHmnfOnDmLI2Jw5PSWTT0RcR3weklTgUFgOvAwsCwibhp7ZADukjQ9Iu6UNB24u836TwdOBxgcHIzZs2eXXsnQ0BCN+Y849qJxRu2eY2au5MThMt+zl7P8kNkdW1YrxW1cB3XLC/XLXLe8UL/Mncpbpo1/BaktvhMWAIcDJ+Sf3+3Qcs3MrKTx9M5ZiqTzgJ8D20v6vaS3kAr+3pJuBl6Rx83MrIs6174wQkQc3OKpl1e1TjMzG92oR/ytzuwxM7N6KtPUc6qkqyS9Q1L1p4uYmVmlynTLvCdwCLAlsFjS1yXtXXkyMzOrRKkvdyPiZuBDpKtsXwZ8TtKvJP2/KsOZmVnnlWnj30nSycAyYC/gNRHxf/PwyRXnMzOzDitzVs/nga8AH4yIJ68vjYg/SPpQZcnMzKwSZQr/vsDDEfE4gKR1gPUj4qGIOLfSdGZm1nFl2vgvBSYXxjfI08zMrIbKFP71c7cNwJNdOGxQXSQzM6tSmcL/oKRdGyOSZpE6azMzsxoq08b/buACSX8ABPw9cFClqczMrDJleue8WtJzge3zpJsi4m/VxjIzs6qU7aTtBaQbpK8L7CqJiDinslS2RmZ04R4Ex8xcudq9DpafsG/l6zWzNTdq4Zd0LrAtcB3weJ4cgAu/mVkNlTniHwR2iFb3aDQzs1opc1bPUtIXumZmthYoc8Q/DbhR0lXAkzdWj4j9KktlZmaVKVP451UdwszMuqfM6ZyXS3o2sF1EXCppA2BS9dHMzKwKZbplfivwTeBLedLmwHeqDGVmZtUp8+XuO4Hdgb/Ckzdl+bsqQ5mZWXXKFP5HI+KxxoikdUnn8ZuZWQ2VKfyXS/ogMDnfa/cC4HvVxjIzs6qUKfzHAvcAw8DbgB+Q7r9rZmY1VOasnieAL+eHmZnVXJm+en5Lkzb9iNimkkRmZlapsn31NKwPvA7YpJo4ZmZWtVHb+CPivsLjjoj4DOkG7GZmVkNlmnp2LYyuQ/oEULYffzMz6zFlCviJheGVwHLg9ZWkMTNro9M3GWp2Q6Fm1rabDJU5q2dON4KYmVl3lGnqeW+75yPipM7FMTOzqpU9q+cFwII8/hrgKuDmqkKZmVl1yhT+LYBdI+IBAEnzgIsi4tAqg5mZWTXKFP4B4LHC+GN52rhJWg48QLp5+8qIGGz/CjMz65Qyhf8c4CpJ387jBwBnd2DdcyLi3g4sx8zMxqDMWT0fl/RDYM886U0RcW21sczMrCpleucE2AD4a0R8Fvi9pK3XcL0B/EjSYklz13BZZmY2Bopof08VSceRzuzZPiKeI2kz4IKI2H3cK5U2j4g7JP0dcAnwroi4YsQ8c4G5AAMDA7Pmz59fevkrVqxg6tSpAAzfcf94Y3bNwGS46+GJTjE2zTLP3HzDiQnD6L/nqrZxle+5uB/3klbbem3Zj3tNcR8b6z4xZ86cxc2+Qy1T+K8DdgGuiYhd8rQlEbFT6bW3X/48YEVEfLrVPIODg7Fo0aLSyxwaGmL27NlA56/0q8IxM1dy4nC9esFolnkir24c7fdc1Tau8j0X9+Ne0mpbry37ca8p7mNj3SckNS38ZZp6Hov03yHygqaUXmvzIFMkPaOwrFcCS9dkmWZmVl6Zf3XfkPQlYCNJbwXezJrdlGUA+Lakxvq/HhH/fw2WZ2ZmY1DmrJ5P53vt/hXYHvhwRFwy3hVGxK3A88f7ejMzWzNtC7+kScCluaO2cRd7MzPrHW3b+CPiceAJSRN3uoaZmXVUmTb+FcCwpEuABxsTI+KoylKZmVllyhT+b+WHWVt1OHW206p8z2VvEmI2Vi0Lv6StIuK2iOhEvzxmZtYj2rXxf6cxIOnCLmQxM7MuaFf4VRjepuogZmbWHe0Kf7QYNjOzGmv35e7zJf2VdOQ/OQ+TxyMinll5OjMz67iWhT8iJnUziJmZdUfZ/vjNzGwt4cJvZtZnXPjNzPqMC7+ZWZ9x4Tcz6zMu/GZmfcaF38ysz7jwm5n1GRd+M7M+48JvZtZnXPjNzPqMC7+ZWZ9x4Tcz6zMu/GZmfcaF38ysz7jwm5n1GRd+M7M+48JvZtZnXPjNzPqMC7+ZWZ9x4Tcz6zMu/GZmfcaF38ysz7jwm5n1mQkp/JJeJekmSbdIOnYiMpiZ9auuF35Jk4AvAK8GdgAOlrRDt3OYmfWriTjifyFwS0TcGhGPAfOB/Scgh5lZX5qIwr85cHth/Pd5mpmZdYEiorsrlA4EXhUR/5LHDwNeFBFHjphvLjA3j24P3DSG1UwD7u1A3G6pW16oX+a65YX6Za5bXqhf5rHmfXZEbDpy4rqdy1PaHcCWhfEt8rSniIjTgdPHswJJiyJicHzxuq9ueaF+meuWF+qXuW55oX6ZO5V3Ipp6rga2k7S1pKcBbwAWTEAOM7O+1PUj/ohYKelI4GJgEnBmRNzQ7RxmZv1qIpp6iIgfAD+ocBXjaiKaQHXLC/XLXLe8UL/MdcsL9cvckbxd/3LXzMwmlrtsMDPrM2tV4a9LVxCSlksalnSdpEV52iaSLpF0c/658QTmO1PS3ZKWFqY1zafkc3mbL5G0aw9lnifpjrydr5O0T+G5D+TMN0n6hwnIu6WkhZJulHSDpKPz9J7czm3y9vI2Xl/SVZKuz5k/kqdvLemXOdv5+SQTJD09j9+Sn5/RI3nPkvTbwjbeOU8f/z4REWvFg/RF8W+AbYCnAdcDO0x0rhZZlwPTRkz7b+DYPHws8MkJzPdSYFdg6Wj5gH2AHwICdgN+2UOZ5wHvazLvDnn/eDqwdd5vJnU573Rg1zz8DODXOVdPbuc2eXt5GwuYmofXA36Zt903gDfk6acB/5qH3wGcloffAJzfI3nPAg5sMv+494m16Yi/7l1B7A+cnYfPBg6YqCARcQXwpxGTW+XbHzgnkl8AG0ma3p2kq7TI3Mr+wPyIeDQifgvcQtp/uiYi7oyIa/LwA8Ay0hXsPbmd2+RtpRe2cUTEijy6Xn4EsBfwzTx95DZubPtvAi+XpC7FbZe3lXHvE2tT4a9TVxAB/EjS4nyFMsBARNyZh/8IDExMtJZa5ev17X5k/hh8ZqH5rKcy5yaFXUhHeD2/nUfkhR7expImSboOuBu4hPTJ4y8RsbJJricz5+fvB541kXkjorGNP5638cmSnj4yb1Z6G69Nhb9O9oiIXUk9lL5T0kuLT0b6HNezp1v1er6CLwLbAjsDdwInTmyc1UmaClwIvDsi/lp8rhe3c5O8Pb2NI+LxiNiZ1EPAC4HnTnCktkbmlbQj8AFS7hcAmwDvX9P1rE2Fv1RXEL0gIu7IP+8Gvk3aIe9qfEzLP++euIRNtcrXs9s9Iu7Kf0hPAF9mVVNDT2SWtB6piP5PRHwrT+7Z7dwsb69v44aI+AuwEHgxqUmkcQ1TMdeTmfPzGwL3dTkq8JS8r8rNbBERjwJfpQPbeG0q/LXoCkLSFEnPaAwDrwSWkrIenmc7HPjuxCRsqVW+BcAb8xkGuwH3F5oqJtSI9s7XkrYzpMxvyGdxbA1sB1zV5WwCzgCWRcRJhad6cju3ytvj23hTSRvl4cnA3qTvJhYCB+bZRm7jxrY/EPhx/tQ1kXl/VTgQEOn7iOI2Ht8+0c1vrat+kL7l/jWpHe8/JjpPi4zbkM52uB64oZGT1JZ4GXAzcCmwyQRmPI/0sf1vpHbDt7TKRzqj4At5mw8Dgz2U+dycaUn+I5lemP8/cuabgFdPQN49SM04S4Dr8mOfXt3ObfL28jbeCbg2Z1sKfDhP34b0T+gW4ALg6Xn6+nn8lvz8Nj2S98d5Gy8FvsaqM3/GvU/4yl0zsz6zNjX1mJlZCS78ZmZ9xoXfzKzPuPCbmfUZF34zsz7jwr+GJG0h6btKvSn+RtJnC739HSHplB7IeICkHQrjx0t6RQeWu/OI3hj3UwW9okoakjSm+4xKOi9f4v6eTucZsZ4xZyu89ghJmxXGl0ua1rl0vU+pd8/3tXjuZ+NY3jclbdNkekf/FiXNlHRWp5bXbS78ayBfUPEt4DsRsR3wHGAq8PEK1zmeu6YdQOotEYCI+HBEXNqBODuTzuVuLHdBRJzQgeWuEUl/D7wgInaKiJNLvmYi7kZ3BLDZaDON1QS9l7YkTRrrayLiJWNcx/NIPYDeOtZ1jVVEDANbSNqq6nVVwYV/zewFPBIRX4XUzwbwHuDNkjbI82yZjwpvlnQcPHn17kVK/W4vlXRQnj5L0uW587aLC1fsDUn6jFLf/f8h6XeS1iks63ZJ60l6q6Sr83IvlLSBpJcA+wGfUurLe1ul/r0PzK9/uaRrle4PcKZyB1D56PMjkq7Jzz2lj5P8qeZ44KC83IOKR1V5HV+U9AtJt0qanZe/rHikJOmVkn6e13OBUl8wzRyW17NU0gsL7/1MpT7Mr5XU6I31R8Dmef498yeTX+RPAN/Wqj7ui9v1aKUrJy/M2/BqSbuPDCFpsqT5+X18G5g82nuR9OG8vKWSTldyIDAI/E/O2VjOu1pt87ys9SV9NT9/raQ5efoRkhZI+jHpArCRr3tjfv/XSzo3T2v6fpWOws/M2+dWSUcVtnez/bbdPvRJSdcAr5N0lFJ//kskzS/E22HkuvLrV+SfsyVdkdd9k6TTlPf/EQ6hcMW7pDdJ+rWkq4DdC9Nfo9Tf/rWSLpU0IGkdpb/RTfM86yj1c7+ppNfl93u9pCsK6/seqYeA+un21XRr0wM4Cji5yfRrSVfhHUG6mvRZpAKxlPTH/k/Alwvzb0jqgvVnwKZ52kGkG9EDDAGnFub/LjCnMN9X8vCzCvN8DHhXHj6LQn/ejXHSlYq3A8/J088hdb4F6Z4Bjde/o7GOEe/zCOCUZuN5HfNJVxfuD/wVmEk62FhM+rQwDbgCmJJf837y1Yoj1jPU2F6kfveX5uH/Ag7NwxuRrtqeAszgqf3yLwFeloePBz7TYrt+ndSBHsBWpO4JRmZ5b+H3shOwMv9OW74XCldhk650fU1h/YOF58ps82MK638ucFv+PR5BumJ5tSu+geflbTOtmKfV+yX1sf8zUl/600j91axH8/12tH3o3wvz/4FVV8lu1G5d+bkV+eds4BHSFbeTSL1sNuuf/nJgZh6enrfNpqT7c1zJqn1zY1bddvZfgBPz8HGF7K8ELszDw8Dmxdx5eHfgexNdh8bz8BF/9S6JiPsi4mFSs9AepB1p73w0tGdE3A9sD+wIXKLULeuHSJ0uNZw/YvigPPyGwnM7SvqJpGHS0c/zRsm2PfDbiPh1Hj+bVFgbGh2HLSYV07H6XqS/kGHgrogYjtSZ1w15ebuRmqCuzO/5cODZLZZ1HjzZ7/4zlfo0eSVwbH7tEKkIPeWjt6QNSX+sl7d4j8Xt+grglLy8BXk9Iz+BvJR02TwRsYT0T4VR3sucfIQ5TPqU2O73Mto236Ow/l8BvyM1MULa15rdk2Av4IKIuDe/rjFPu/d7UaS+9O8ldRQ3QOv9tt0+VNy+S0ifcA4l/cNsaLauka6KdK+Nx0n7wh5N5pkO3JOHXwQMRcQ9ke7PUcyxBXBx/n38G6t+H2cCb8zDbyZ1iAbpn8ZZkt5K+sfTcDcVNNV1Q8+1BdbMjazq7AkASc8kFZ9bSHeEGtknRkTEr5Vuk7YP8DFJl5F66bwhIl7cYl0PFoYXAP8laRNgFqkvD0hH2QdExPWSjiAdKa2JR/PPxxnfvtJ4/ROF4cb4unm5l0TEwSWWtdp2JH2a+KeIuKn4hMZ2y7zidl0H2C0iHhnD659cLU3ei6T1gVNJR/a3S5pH+gfVypps8wdHn+Upmr5fpXuPFH9fjwPrtthvR+tMsJhpX9I/hdeQmixn5umrravJcpr9/kd6mPbbtuHzwEkRsUDSbNKnDvLv5y5Je5F6wDwkT3+7pBfl/IslzYqI+/K6Hi6xvp7jI/41cxmwgaQ3wpNfYJ0InBURD+V59la6j+pk0pesVyqdyfFQRHwN+BTpH8RNwKaSXpyXtZ7Sl1WriXSXnquBzwLfz0dBkG6Jd6dS97mHFF7yQH5upJuAGZL+Tx4/jPRxuaxWyy3rF8DujfXnNuTntJi30Z68B6kXwvuBi0lt4srP7TLyRXm+P0vaM09q9x5/BLyrMaJ8b9MRrgD+OT+/I6m5p917aRSie/PRdPFAYTzb7yfk321e/lak32M7Pya1sT8rv26TPL3M+6XwfKv9dtR9KLfJbxkRC0nNYBuSToQo64VKPe+uQ9oXftpknmVAI8cvgZdJelb+e3hdYb4NWdV98eE81VdIn6guaPxdSdo2In4ZER8mfaJodIX8HFb1lFkrLvxrIDdjvJb0R3UzqR31EeCDhdmuIvVhvoTUZriI1NZ9Vf6IfRzwsfxx9EDgk5KuJ/V+2O6shvOBQ3nqR9j/JO3wVwK/KkyfD/xb/jJr20L+R4A3ARfkj71PkO5BWtZC0hdz1zW+6BuLiLiH1DZ9nqQlwM9pfaOMRyRdm/O9JU/7KKnteYmkG/J4M4eTvtxeQvpu4fgW8x0FDOYvH28E3t5kni8CUyUty7feefQAAAD/SURBVMtZ3O69ROpX/cukAnEx6R92w1nAaXrql7ujORVYJ/++zgeOiNRPe0sRcQPpTLPL877V6Fa5zPstarbflt2HJgFfy/NcC3wub5uyrgZOIRX335I+IY90EflTbqTuieeRfg9X5tc1zMt5FwP3jljGAtI/pK8Wpn1K6YvrpaTvI67P0+fkddaOe+c0s56Wm2PeFxH/OMp8k0kHI7sXPgWPdV2DpBM29hxlvqeTPtnsEatu41gbPuI3s7VCPoHiOMZ5b1+liw8vJN3qcDRbAcfWseiDj/jNzPqOj/jNzPqMC7+ZWZ9x4Tcz6zMu/GZmfcaF38ysz7jwm5n1mf8FITNv65RX+QIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data.Time.hist();\n",
    "plt.xlabel(\"Observation time before death or censorship (days)\");\n",
    "plt.ylabel(\"Frequency (number of patients)\");\n",
    "# Note that the semicolon at the end of the plotting line\n",
    "# silences unnecessary textual output - try removing it\n",
    "# to observe its effect"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ohXd70ZBaRWJ"
   },
   "source": [
    "<a name='Ex-1'></a>\n",
    "### Exercise 1\n",
    "\n",
    "In the next cell, write a function to compute the fraction ($\\in [0, 1]$) of observations which were censored. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>    \n",
    "<summary>\n",
    "    <font size=\"3\" color=\"darkgreen\"><b>Hints</b></font>\n",
    "</summary>\n",
    "<p>\n",
    "<ul>\n",
    "    <li>Summing up the <code>'Event'</code> column will give you the number of observations where censorship has NOT occurred.</li>\n",
    "    \n",
    "</ul>\n",
    "</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "9fkHnfJ6bD0f",
    "outputId": "21ec0301-5884-48f2-8ed1-10795d567c31"
   },
   "outputs": [],
   "source": [
    "# UNQ_C1 (UNIQUE CELL IDENTIFIER, DO NOT EDIT)\n",
    "def frac_censored(df):\n",
    "    \"\"\"\n",
    "    Return percent of observations which were censored.\n",
    "    \n",
    "    Args:\n",
    "        df (dataframe): dataframe which contains column 'Event' which is \n",
    "                        1 if an event occurred (death)\n",
    "                        0 if the event did not occur (censored)\n",
    "    Returns:\n",
    "        frac_censored (float): fraction of cases which were censored. \n",
    "    \"\"\"\n",
    "    result = 0.0\n",
    "    \n",
    "    ### START CODE HERE ###\n",
    "    censored_count = sum(df['Event'] == 0)\n",
    "    result = censored_count/len(df)\n",
    "    ### END CODE HERE ###\n",
    "    \n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.325\n"
     ]
    }
   ],
   "source": [
    "print(frac_censored(data))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Expected Output:\n",
    "```CPP\n",
    "0.325\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "BpzYWXUpbk6x"
   },
   "source": [
    "Run the next cell to see the distributions of survival times for censored and uncensored examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 545
    },
    "colab_type": "code",
    "id": "1k3qlTQLbulW",
    "outputId": "de331041-6612-4df7-953d-41efbc741e49"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEWCAYAAABsY4yMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAXm0lEQVR4nO3de5hkdX3n8fdHIDIyCFGwVxAdNd5Yxws0RlejMxojgtesGzHquokra+IlrmQf8fIo+sSsJuIta6KIrrfECaIYleiCl5G4BnFGgQERRZmIiCAqlwaW63f/qNNStN3T1T11+vLj/XqefqbqnFO/3/fXp+ozp3916lSqCklSe+6w3AVIkvphwEtSowx4SWqUAS9JjTLgJalRBrwkNcqAl1aIJBuS/Hi561A7DHitakn+MMmWJFNJLkny+SSPWe66pJXAgNeqleSVwDuBvwQmgHsCfws8fTnrAkiy63LXIBnwWpWS7AW8CXhJVX2qqq6pqhur6rNV9T+S3CHJ0Ul+kOTnSU5IcpfuseuSVJIXJPlRksuTvHao7Ud0fxVcleTSJG8fWve0JOcmuSLJ5iQPGlq3PcmrkpwNXJNk1yT7Jflkkp8luTDJy4e2X5PkQ0l+meQ7wCFL8bvT7YcBr9XqUcDuwElzrH8Z8AzgccB+wC+B98zY5jHAA4AnAK8fCut3Ae+qqjsD9wVOAEhyf+DjwCuAfYF/Bj6b5DeG2nwOcDiwN3AL8FngLGD/rp9XJHlSt+0buvbvCzwJeMGCfgPSPAx4rVZ3BS6vqpvmWP9i4LVV9eOquh44BnjWjKmTN1bVdVV1FoMQfmi3/Ebgt5LsU1VTVXV6t/zZwMlVdWpV3Qi8DVgD/IehNt9dVRdV1XUMjsj3rao3VdUNVfVD4P3AEd22fwC8uap+UVUXAe9e/K9D+nUGvFarnwP77GCu+17ASd1UyhXAecDNDObqp/106Pa1wNru9guB+wPfTfLNJE/plu8H/Nv0A6rqFuAiBkfn0y6aUcN+0zV0dbxmqIb9Zmz/b0hjZMBrtfpX4HoG0zCzuQh4clXtPfSze1VdPF/DVfX9qnoOcDfgrcCJSfYAfsIgtAFIEuAAYLjN4cuzXgRcOKOGPavqsG79Jd3jp91zvtqkhTDgtSpV1ZXA64H3JHlGkjsl2S3Jk5P8FfBe4M1J7gWQZN8kI51dk+R5SfbtjtCv6BbfwmAu/vAkT0iyG3AUg/9kvj5HU2cAV3dvvK5JskuSByeZfjP1BODVSX4zyT0YvG8gjY0Br1Wrqo4FXgm8DvgZgyPmlwKfZvBG6WeAU5JcDZwO/PaITR8KnJtkqmvniG6u/nzgecDfAJcDTwWeWlU3zFHfzcBTgIcBF3aPOR7Yq9vkjQymZS4ETgE+OvLgpRHEL/yQpDZ5BC9JjTLgJalRBrwkNcqAl6RGragLIu2zzz61bt263tq/5ppr2GOPPXprf6k4jpWjhTGA41hpFjKOrVu3Xl5V+862bkUF/Lp169iyZUtv7W/evJkNGzb01v5ScRwrRwtjAMex0ixkHEnm/AS0UzSS1CgDXpIaZcBLUqMMeElqlAEvSY0y4CWpUb2eJplkO3A1gy9auKmqJvvsT5J0q6U4D35jVV2+BP1IkoY4RSNJjer1evBJLmTwbfYFvK+qjptlmyOBIwEmJiYO3rRp06L62nbxlfNuM7EGLr1uUc3Paf3+e82/0ZhNTU2xdu3a+Tdc4VoYRwtjgJU1jlFey3Pp4zW+FGbmyEL2x8aNG7fONf3dd8DvX1UXJ7kbcCrwsqo6ba7tJycna7GXKlh39MnzbnPU+ps4dtt4Z6W2v+XwsbY3itvjx7FXqhbGACtrHKO8lufSx2t8KczMkQVeqmDOgO91imb6C46r6jLgJOARffYnSbpVbwGfZI8ke07fBn4POKev/iRJt9Xn3zITwElJpvv5h6r6Qo/9SZKG9BbwVfVD4KF9tS9J2jFPk5SkRhnwktQoA16SGmXAS1KjDHhJapQBL0mNMuAlqVEGvCQ1yoCXpEYZ8JLUKANekhplwEtSowx4SWqUAS9JjTLgJalRBrwkNcqAl6RGGfCS1CgDXpIaZcBLUqMMeElqlAEvSY0y4CWpUQa8JDXKgJekRhnwktQoA16SGmXAS1KjDHhJapQBL0mNMuAlqVEGvCQ1qveAT7JLkm8n+VzffUmSbrUUR/B/Bpy3BP1Ikob0GvBJ7gEcDhzfZz+SpF+Xquqv8eRE4H8CewJ/XlVPmWWbI4EjASYmJg7etGnTovradvGV824zsQYuvW5Rzc9p/f57jbfBEUxNTbF27dol73fcWhhHC2OA2ccxymtqpenjNb4UZubIQp5XGzdu3FpVk7Ot23XnS5tdkqcAl1XV1iQb5tquqo4DjgOYnJysDRvm3HSH/svRJ8+7zVHrb+LYbeMd8vbnbhhre6PYvHkzi/09rSQtjKOFMcDs4xjlNbXS9PEaXwozc2Rcz6s+p2geDTwtyXZgE/D4JB/rsT9J0pDeAr6qXl1V96iqdcARwJer6nl99SdJui3Pg5ekRi3JZFVVbQY2L0VfkqQBj+AlqVEGvCQ1yoCXpEYZ8JLUKANekhplwEtSowx4SWqUAS9JjTLgJalRBrwkNcqAl6RGGfCS1CgDXpIaZcBLUqMMeElqlAEvSY0y4CWpUQa8JDXKgJekRhnwktQoA16SGmXAS1KjDHhJapQBL0mNMuAlqVEGvCQ1yoCXpEaNFPBJ1vddiCRpvEY9gv/bJGck+dMke/VakSRpLEYK+Kr6HeC5wAHA1iT/kOSJvVYmSdopI8/BV9X3gdcBrwIeB7w7yXeT/H5fxUmSFm/UOfiHJHkHcB7weOCpVfWg7vY7eqxPkrRIu4643d8AxwOvqarrphdW1U+SvG62ByTZHTgNuGPXz4lV9YadrFeSNKJRA/5w4LqquhkgyR2A3avq2qr66ByPuR54fFVNJdkN+FqSz1fV6TtftiRpPqPOwX8RWDN0/07dsjnVwFR3d7fupxZcoSRpUUYN+N2Hwpru9p3me1CSXZKcCVwGnFpV31hcmZKkhUrV/AfVSf4v8LKq+lZ3/2Dgf1XVo0bqJNkbOKlr45wZ644EjgSYmJg4eNOmTQsbQWfbxVfOu83EGrj0unk3W5D1+y/9xwKmpqZYu3btkvc7bi2Mo+8xjPK8Hoc+XhvLYbWOY2aOLOR5tXHjxq1VNTnbulED/hBgE/ATIMC/A55dVVtHqmDQxuuBa6vqbXNtMzk5WVu2bBm1ydtYd/TJ825z1PqbOHbbqG87jGb7Ww4fa3uj2Lx5Mxs2bFjyfsethXH0PYZRntfj0MdrYzms1nHMzJGFPK+SzBnwI/0mquqbSR4IPKBbdH5V3ThPp/sCN1bVFUnWAE8E3jpSxZKknbaQ/+oOAdZ1jzkoCVX1kR1sf3fgw0l2YTDXf0JVfW7RlUqSFmSkgE/yUeC+wJnAzd3iAuYM+Ko6G3j4zhYoSVqcUY/gJ4EDa5QJe0nSijDqaZLnMHhjVZK0Sox6BL8P8J0kZzD4hCoAVfW0XqqSJO20UQP+mD6LkCSN36inSX41yb2A+1XVF5PcCdil39IkSTtj1MsFvwg4EXhft2h/4NN9FSVJ2nmjvsn6EuDRwFXwqy//uFtfRUmSdt6oAX99Vd0wfSfJrnhlSEla0UYN+K8meQ2wpvsu1k8An+2vLEnSzho14I8GfgZsA/4b8M8Mvp9VkrRCjXoWzS3A+7sfSdIqMOq1aC5kljn3qrrP2CuSJI3FQq5FM2134D8Bdxl/OZKkcRlpDr6qfj70c3FVvZPBF3FLklaoUadoDhq6ewcGR/Sr72tTJOl2ZNSQPnbo9k3AduAPxl6NJGlsRj2LZmPfhUiSxmvUKZpX7mh9Vb19POVIksZlIWfRHAJ8prv/VOAM4Pt9FCVJ2nmjBvw9gIOq6mqAJMcAJ1fV8/oqTJK0c0a9VMEEcMPQ/Ru6ZZKkFWrUI/iPAGckOam7/wzgw/2UJEkah1HPonlzks8Dv9Mt+qOq+nZ/ZUmSdtaoUzQAdwKuqqp3AT9Ocu+eapIkjcGoX9n3BuBVwKu7RbsBH+urKEnSzhv1CP6ZwNOAawCq6ifAnn0VJUnaeaMG/A1VVXSXDE6yR38lSZLGYdSAPyHJ+4C9k7wI+CJ++YckrWijnkXztu67WK8CHgC8vqpO7bUySdJOmTfgk+wCfLG74JihLkmrxLxTNFV1M3BLkr2WoB5J0piM+knWKWBbklPpzqQBqKqX91KVJGmnjRrwn+p+JEmrxA4DPsk9q+pHVbXg684kOYDBNWwmGJxeeVz3KVhJ0hKYbw7+09M3knxygW3fBBxVVQcCjwRekuTABbYhSVqk+QI+Q7fvs5CGq+qSqvpWd/tq4Dxg/4WVJ0larAw+oDrHyuRbVXXQzNsL7iRZB5wGPLiqrpqx7kjgSICJiYmDN23atJgu2HbxlfNuM7EGLr1uUc3Paf3+S39y0dTUFGvXrl3yfsethXH0PYZRntfj0MdrYzms1nHMzJGFPK82bty4taomZ1s3X8DfzOCsmQBrgGunVwFVVXeer/Mka4GvAm+uqh2+UTs5OVlbtmyZr8lZrTv65Hm3OWr9TRy7bdT3lUez/S2Hj7W9UWzevJkNGzYseb/j1sI4+h7DKM/rcejjtbEcVus4ZubIQp5XSeYM+B3+JqpqlxHrm6vj3YBPAn8/X7hLksZrIdeDX5AkAT4AnFdVb++rH0nS7HoLeODRwPOBxyc5s/s5rMf+JElDepusqqqvcduzcCRJS6jPI3hJ0jIy4CWpUQa8JDXKgJekRhnwktQoA16SGmXAS1KjDHhJapQBL0mNMuAlqVEGvCQ1yoCXpEYZ8JLUKANekhplwEtSowx4SWqUAS9JjTLgJalRBrwkNcqAl6RGGfCS1CgDXpIaZcBLUqMMeElqlAEvSY0y4CWpUQa8JDXKgJekRhnwktQoA16SGmXAS1KjDHhJalRvAZ/kg0kuS3JOX31IkubW5xH8h4BDe2xfkrQDvQV8VZ0G/KKv9iVJO5aq6q/xZB3wuap68A62ORI4EmBiYuLgTZs2LaqvbRdfOe82E2vg0usW1fyc1u+/13gbHMHU1BRr164dacwrWR/7Y6m1MAZwHMttZo5Mv8ZHsXHjxq1VNTnbumUP+GGTk5O1ZcuWRfW17uiT593mqPU3cey2XRfV/ly2v+XwsbY3is2bN7Nhw4aRxryS9bE/lloLYwDHsdxm5sj0a3wUSeYMeM+ikaRGGfCS1Kg+T5P8OPCvwAOS/DjJC/vqS5L063qbrKqq5/TVtiRpfk7RSFKjDHhJapQBL0mNMuAlqVEGvCQ1yoCXpEYZ8JLUKANekhplwEtSowx4SWqUAS9JjTLgJalRBrwkNcqAl6RGGfCS1CgDXpIaZcBLUqMMeElqlAEvSY0y4CWpUQa8JDXKgJekRhnwktQoA16SGmXAS1KjDHhJapQBL0mNMuAlqVEGvCQ1yoCXpEYZ8JLUKANekhplwEtSo3oN+CSHJjk/yQVJju6zL0nSbfUW8El2Ad4DPBk4EHhOkgP76k+SdFt9HsE/Arigqn5YVTcAm4Cn99ifJGlIqqqfhpNnAYdW1X/t7j8f+O2qeumM7Y4EjuzuPgA4v5eCBvYBLu+x/aXiOFaOFsYAjmOlWcg47lVV+862Ytfx1bM4VXUccNxS9JVkS1VNLkVffXIcK0cLYwDHsdKMaxx9TtFcDBwwdP8e3TJJ0hLoM+C/Cdwvyb2T/AZwBPCZHvuTJA3pbYqmqm5K8lLg/wC7AB+sqnP76m9ESzIVtAQcx8rRwhjAcaw0YxlHb2+ySpKWl59klaRGGfCS1KhmAz7J9iTbkpyZZEu37C5JTk3y/e7f31zuOmdK8sEklyU5Z2jZrHVn4N3dpSDOTnLQ8lV+W3OM45gkF3f75Mwkhw2te3U3jvOTPGl5qv51SQ5I8pUk30lybpI/65avqn2yg3Gsqn2SZPckZyQ5qxvHG7vl907yja7ef+xO7CDJHbv7F3Tr1y1n/V1Nc43hQ0kuHNoXD+uWL/45VVVN/gDbgX1mLPsr4Oju9tHAW5e7zlnqfixwEHDOfHUDhwGfBwI8EvjGctc/zziOAf58lm0PBM4C7gjcG/gBsMtyj6Gr7e7AQd3tPYHvdfWuqn2yg3Gsqn3S/V7Xdrd3A77R/Z5PAI7olr8X+JPu9p8C7+1uHwH84woew4eAZ82y/aKfU80ewc/h6cCHu9sfBp6xjLXMqqpOA34xY/FcdT8d+EgNnA7sneTuS1Ppjs0xjrk8HdhUVddX1YXABQwudbHsquqSqvpWd/tq4Dxgf1bZPtnBOOayIvdJ93ud6u7u1v0U8HjgxG75zP0xvZ9OBJ6QJEtU7qx2MIa5LPo51XLAF3BKkq3d5RAAJqrqku72T4GJ5Sltweaqe3/goqHtfsyOX7QrwUu7PzM/ODRFtirG0f15/3AGR1yrdp/MGAessn2SZJckZwKXAacy+Oviiqq6qdtkuNZfjaNbfyVw16Wt+NfNHENVTe+LN3f74h1J7tgtW/S+aDngH1NVBzG4muVLkjx2eGUN/vZZdeeIrta6O38H3Bd4GHAJcOzyljO6JGuBTwKvqKqrhtetpn0yyzhW3T6pqpur6mEMPh3/COCBy1zSgs0cQ5IHA69mMJZDgLsAr9rZfpoN+Kq6uPv3MuAkBk+ES6f/tOn+vWz5KlyQuepeVZeDqKpLuyf2LcD7ufVP/hU9jiS7MQjFv6+qT3WLV90+mW0cq3WfAFTVFcBXgEcxmLaY/uDmcK2/Gke3fi/g50tc6pyGxnBoN41WVXU98L8Zw75oMuCT7JFkz+nbwO8B5zC4VMILus1eAPzT8lS4YHPV/RngP3fvsj8SuHJo2mDFmTFv+EwG+wQG4ziiO+Ph3sD9gDOWur7ZdPO1HwDOq6q3D61aVftkrnGstn2SZN8ke3e31wBPZPB+wleAZ3Wbzdwf0/vpWcCXu7+4ls0cY/ju0AFDGLyHMLwvFvecWu53lPv4Ae7D4AyAs4Bzgdd2y+8KfAn4PvBF4C7LXesstX+cwZ/KNzKYa3vhXHUzeFf9PQzmILcBk8td/zzj+GhX59ndk/buQ9u/thvH+cCTl7v+oboew2D65WzgzO7nsNW2T3YwjlW1T4CHAN/u6j0HeH23/D4M/gO6APgEcMdu+e7d/Qu69fdZwWP4crcvzgE+xq1n2iz6OeWlCiSpUU1O0UiSDHhJapYBL0mNMuAlqVEGvCQ1yoDXqpLkrkNX2/vpjCshfr2nPh+e5ANzrNueZJ8x9rUpyf3G1Z5u3zxNUqtWkmOAqap6W8/9fAL4i6o6a5Z12xmcl3z5mPp6HPC8qnrRONrT7ZtH8GpGkqnu3w1Jvprkn5L8MMlbkjy3uwb3tiT37bbbN8knk3yz+3n0LG3uCTxkOty7vyBO6a7jfTyDD6FMb/vp7uJ2505f4C7JHyd559A2L+ouJLVHkpMzuCb4OUme3W3yL8DvDn3sXlo0A16teijwYuBBwPOB+1fVI4DjgZd127wLeEdVHQL8x27dTJPc+pFxgDcAX6uqf8/gGkf3HFr3x1V1cPeYlye5K4PrlD+1uw4MwB8BHwQOBX5SVQ+tqgcDXwCowTVhLujql3aKRwlq1Teru15Hkh8Ap3TLtwEbu9u/Cxw4dHnwOydZW7deqxsGX5Txs6H7jwV+H6CqTk7yy6F1L0/yzO72AcD9qur0JF8GnpLkPGC3qtqW5Hrg2CRvBT5XVf8y1M5lwH7A1kWPXsKAV7uuH7p9y9D9W7j1eX8H4JFV9f920M51DK5nskNJNjD4D+NRVXVtks1DjzseeA3wXQZXCaSqvpfBV68dBvxFki9V1Zu67Xfv+pV2ilM0uj07hVuna0j3HZgznAf81tD904A/7LZ/MjD9BRl7Ab/swv2BDL5aDYAafJnDAd3jPt49dj/g2qr6GPDXDL7ecNr9ue20kLQoHsHr9uzlwHuSnM3gtXAag3n7X6mq7ybZK8meNfiquzcCH09yLvB14Efdpl8AXtxNw5wPnD6jrxOAh1XV9JTOeuCvk9zC4IqbfwKQZAK4rqp+Ouax6nbI0ySleST578DVVTXbm7CjtvE5Bm/ofmmEvq6qqlnPu5cWwikaaX5/x23n9EeWZO8k32NwVL7DcO9cwa1fEi3tFI/gJalRHsFLUqMMeElqlAEvSY0y4CWpUQa8JDXq/wNLEYKHrgH38wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAaK0lEQVR4nO3de5xdZX3v8c+XABIZGoTgFgI6oEhLM0LJ1tJidUapRS6i59BXRaCAlNQeFdoTjo3SHqhHT/ESLT2ltimmIiBT5CYl1nIpA7VcEwSGi9w0hQRI5GLCQLiE/M4fa0X27O7J7JnZa63Mfr7v12tes9ez1l7r+c2a+c7aa6/9LEUEZmaWjq2q7oCZmZXLwW9mlhgHv5lZYhz8ZmaJcfCbmSXGwW9mlhgHv1kXkLRC0sFV98OmBwe/TUuSQtLbmtrOlHRBVX0ymy4c/GZbGElbV90H624OfutKkvolrZS0QNIaSU9IOrFh/kxJiyT9p6S1kn4oaWY+70BJN0n6uaS7JPU3PG9I0v+R9B+SnpN0taTZ+bztJF0g6en8ubdLquXzdpN0paRnJD0s6eSGdZ4p6ZL8ueuAEyRtJWmhpEfy9V0saaeG5xyX9/1pSacX/xO1buLgt272JmAWMAc4CThH0hvyeV8F5gG/CewEfAbYKGkOsBT4Qt5+GnCppF0a1vsx4ETgjcC2+TIAx+fb2wPYGfgEsD6fNwisBHYDjgL+r6T3NazzSOASYEfgQuDTwIeB9+bPeRY4B0DSvsA3gOPyeTsDu0/uR2QpcvBbN3sF+HxEvBIR3wdGgH0kbQV8HDg1IlZFxKsRcVNEvAQcC3w/Ir4fERsj4hpgGXBow3r/MSIejIj1wMXA/g3b2xl4W77O5RGxTtIewEHAn0bEixFxJ3Au8PsN67w5Iq7It7me7J/G6RGxMu/XmcBR+Wmgo4CrIuLGfN6fAxs7/tOzruVziTZdvQps09S2DVn4bvJ0RGxomH4B6AFmA9sBj7RY71uA35V0RNN6r2+YfrLFOgHOJzvaH5S0I3ABcDrZUfkzEfFcw/P+E6g3TD/Woh+XS2oM9FeBWr6+XywfEc9LerpFLWYt+YjfpqtHgd6mtj3JAnU8TwEvAm9tMe8x4PyI2LHha/uIOGu8leavLP4iIvYlO4V0ONlR/ePATpJ2aFj8zcCqxqe36McHm/qxXUSsAp4g+wcDgKTXk73SMGuLg9+mq38C/kzS7vkboQcDR5CdJ9+siNgILAG+lr/pOkPSb0h6HdlR+hGSfidv3y5/o3jcc+iSBiT1SZoBrCN79bExIh4DbgL+Ml/fO8jec9jcpad/B3xR0lvyde8i6ch83iXA4ZLeLWlb4PP4b9kmwL8sNl19nixMf0j2xueXgWMi4p42n38aMAzcDjwDfAnYKg/pI4HPAT8jO/L+X7T3t/ImslBeB9wP3EB2+gfgaLJXKI8DlwNnRMS1m1nX2cCVwNWSngNuAX4dICLuBT4JfIfs6P9ZsjeOzdoi34jFzCwtPuI3M0uMg9/MLDEOfjOzxDj4zcwSMy0+wDV79uzo7e0dd7nnn3+e7bffvvgObUFSrBnSrNs1p6GTNS9fvvypiNiluX1aBH9vby/Lli0bd7mhoSH6+/uL79AWJMWaIc26XXMaOlmzpJYfaPSpHjOzxDj4zcwS4+A3M0uMg9/MLDEOfjOzxDj4zcwSU1jwS1qS3+v0nqb2T0v6saR7JX25qO2bmVlrRR7xfws4pLFB0gDZkLf7RcSvkt331MzMSlRY8EfEjWTjnDf6I+Cs/D6hRMSaorZvZmatFToev6ResptCz82n7wS+R/ZK4EXgtIi4fYznzgfmA9RqtXmDg4Pjbm9kZISenp5RbcOr1k6+gCnqmzOr8G20qjkFKdbtmtPQyZoHBgaWR0S9ub3sIRu2BnYCDgTeCVwsaa9o8d8nIhYDiwHq9Xq08xHmVh91PmHh0il3erJWHNNf+DZS/Eg7pFm3a05DGTWXfVXPSuCyyNwGbARml9wHM7OklR38VwADAJLeDmwLPFVyH8zMklbYqR5JFwH9wGxJK4EzgCXAkvwSz5eB41ud5jEzs+IUFvwRcfQYs44taptmZjY+f3LXzCwxDn4zs8Q4+M3MEuPgNzNLjIPfzCwxDn4zs8Q4+M3MElP2WD1J6S1hnKAFfRv+y3hEK846rPDtmtn05SN+M7PEOPjNzBLj4DczS4yD38wsMQ5+M7PEOPjNzBLj4DczS4yD38wsMYUFv6Qlktbkd9tqnrdAUkjy/XbNzEpW5BH/t4BDmhsl7QF8AHi0wG2bmdkYCgv+iLgReKbFrK8DnwF8r10zswqUeo5f0pHAqoi4q8ztmpnZaxRR3IG3pF7gqoiYK+n1wPXAByJiraQVQD0inhrjufOB+QC1Wm3e4ODguNsbGRmhp6dnVNvwqrVTKWGLV5sJq9ePbuubM6uazpSo1b7udq45DZ2seWBgYHlE1Jvbywz+PuA64IV89u7A48C7IuLJza2nXq/HsmXLxt3e0NAQ/f39o9rKGCGzSgv6NrBoePQgqymMztlqX3c715yGTtYsqWXwlzYsc0QMA29s6NAKNnPEb2ZmxSjycs6LgJuBfSStlHRSUdsyM7P2FXbEHxFHjzO/t6htm5nZ2PzJXTOzxDj4zcwS4+A3M0uMg9/MLDEOfjOzxDj4zcwS4+A3M0uMg9/MLDEOfjOzxDj4zcwS4+A3M0uMg9/MLDEOfjOzxDj4zcwS4+A3M0uMg9/MLDEOfjOzxBR568UlktZIuqeh7SuSfizpbkmXS9qxqO2bmVlrRR7xfws4pKntGmBuRLwDeBD4bIHbNzOzFgoL/oi4EXimqe3qiNiQT94C7F7U9s3MrDVFRHErl3qBqyJibot5/wz8U0RcMMZz5wPzAWq12rzBwcFxtzcyMkJPT8+otuFVayfc7+mkNhNWrx/d1jdnVjWdKVGrfd3tXHMaOlnzwMDA8oioN7dv3ZG1T5Ck04ENwIVjLRMRi4HFAPV6Pfr7+8dd79DQEM3LnbBw6RR6uuVb0LeBRcOjd+OKY/qr6UyJWu3rbuea01BGzaUHv6QTgMOB90eRLzfMzKylUoNf0iHAZ4D3RsQLZW7bzMwyRV7OeRFwM7CPpJWSTgL+BtgBuEbSnZL+rqjtm5lZa4Ud8UfE0S2av1nU9szMrD3+5K6ZWWIc/GZmiXHwm5klxsFvZpYYB7+ZWWIc/GZmiXHwm5klxsFvZpYYB7+ZWWIc/GZmiXHwm5klxsFvZpYYB7+ZWWIc/GZmiXHwm5klxsFvZpYYB7+ZWWLaCn5JfRNdsaQlktZIuqehbSdJ10h6KP/+homu18zMpqbdI/6/lXSbpP8haVabz/kWcEhT20LguojYG7gunzYzsxK1FfwR8VvAMcAewHJJ35H02+M850bgmabmI4Hz8sfnAR+eWHfNzGyqFBHtLyzNIAvrvwbWAQI+FxGXjbF8L3BVRMzNp38eETvmjwU8u2m6xXPnA/MBarXavMHBwXH7NzIyQk9Pz6i24VVr2ylt2qrNhNXrR7f1zWn3Rdn01WpfdzvXnIZO1jwwMLA8IurN7Vu382RJ7wBOBA4DrgGOiIg7JO0G3Ay0DP7NiYiQNOZ/nYhYDCwGqNfr0d/fP+46h4aGaF7uhIVLJ9q1aWVB3wYWDY/ejSuO6a+mMyVqta+7nWtOQxk1t3uO//8BdwD7RcQnI+IOgIh4HPizCWxvtaRdAfLvaybSWTMzm7p2g/8w4DsRsR5A0laSXg8QEedPYHtXAsfnj48HvjeB55qZWQe0G/zXAjMbpl+ft41J0kVkp4H2kbRS0knAWcBvS3oIODifNjOzErV1jh/YLiJGNk1ExMimI/6xRMTRY8x6f7udMzOzzmv3iP95SQdsmpA0D1i/meXNzGwL1e4R/x8D35X0ONklnG8Cfq+wXpmZWWHaCv6IuF3SLwP75E0PRMQrxXXLzMyK0u4RP8A7gd78OQdIIiK+XUivzMysMO1+gOt84K3AncCreXMADn4zs2mm3SP+OrBvTGR8BzMz2yK1e1XPPWRv6JqZ2TTX7hH/bOA+SbcBL21qjIgPFdIrMzMrTLvBf2aRnTAzs/K0eznnDZLeAuwdEdfmn9qdUWzXzMysCO3eevFk4BLg7/OmOcAVRXXKzMyK0+6bu58EDiK7+QoR8RDwxqI6ZWZmxWk3+F+KiJc3TUjamuw6fjMzm2baDf4bJH0OmJnfa/e7wD8X1y0zMytKu8G/EPgZMAz8IfB9JnbnLTMz20K0e1XPRuAf8i8zM5vG2h2r56e0OKcfEXt1vEdmZlaoiYzVs8l2wO8CO012o5L+BPgDsn8mw8CJEfHiZNdnZmbta+scf0Q83fC1KiL+iuwG7BMmaQ5wClCPiLlkHwT76GTWZWZmE9fuqZ4DGia3InsFMJGx/Fttd6akV8hu3P74FNZlZmYToHZGWpZ0fcPkBmAF8NWIeGBSG5VOBb5Idt/eqyPimBbLzAfmA9RqtXmDg4PjrndkZISenp5RbcOr1k6mi9NGbSasbrr7cd+cWdV0pkSt9nW3c81p6GTNAwMDyyOi3tzeVvB3kqQ3AJeS3bP352SfCbgkIi4Y6zn1ej2WLVs27rqHhobo7+8f1da7cOlUurvFW9C3gUXDo198rThrUmfhppVW+7rbueY0dLJmSS2Dv91TPf9zc/Mj4msT6MvBwE8j4mf5ui8DfhMYM/jNzKxzJnJVzzuBK/PpI4DbgIcmsc1HgQPzET7XA+8Hxj+cNzOzjmg3+HcHDoiI5wAknQksjYhjJ7rBiLhV0iXAHWTvF/wIWDzR9ZiZ2eS0G/w14OWG6ZfztkmJiDOAMyb7fDMzm7x2g//bwG2SLs+nPwycV0yXzMysSO2O1fNFSf8C/FbedGJE/Ki4bpmZWVHaHZ0Tsg9arYuIs4GVkvYsqE9mZlagdm+9eAbwp8Bn86Zt8OWXZmbTUrtH/B8BPgQ8DxARjwM7FNUpMzMrTrvB/3JkH/ENAEnbF9clMzMrUrvBf7Gkvwd2lHQycC2+KYuZ2bTU7lU9X83vtbsO2Af43xFxTaE9s2mnzHGRFvRt4ISG7aUwPpFZp4wb/JJmANdGxADgsDczm+bGPdUTEa8CGyV1/1i/ZmYJaPeTuyPAsKRryK/sAYiIUwrplZmZFabd4L8s/zIzs2lus8Ev6c0R8WhEeFweM7MuMd45/is2PZB0acF9MTOzEowX/Gp4vFeRHTEzs3KMF/wxxmMzM5umxntzdz9J68iO/Gfmj8mnIyJ+qdDemZlZx202+CNiRhEblbQjcC4wl+yVxMcj4uYitmVmZqO1ezlnp50N/CAijpK0LdlY/2ZmVoLSgz//BPB7gBMAIuJlRt/P18zMCqRstOUSNyjtDywG7gP2A5YDp0bE803LzQfmA9RqtXmDg4PjrntkZISenp5RbcOr1nam41uo2kxYvX50W9+cakbXKPNn3Vx3VTWXqdXvd7dzzVMzMDCwPCLqze1VBH8duAU4KCJulXQ22S0d/3ys59Tr9Vi2bNm46x4aGqK/v39UW5kjRlZhQd8GFg2PfuFW1UiVZY/O2Vh3CqNztvr97naueWoktQz+idxzt1NWAisj4tZ8+hLggAr6YWaWpNKDPyKeBB6TtE/e9H6y0z5mZlaCqq7q+TRwYX5Fz0+AEyvqh5lZcioJ/oi4E/gv553MzKx4VZzjNzOzCjn4zcwS4+A3M0uMg9/MLDEOfjOzxDj4zcwS4+A3M0uMg9/MLDEOfjOzxDj4zcwS4+A3M0uMg9/MLDEOfjOzxDj4zcwS4+A3M0uMg9/MLDEOfjOzxFQW/JJmSPqRpKuq6oOZWYqqPOI/Fbi/wu2bmSWpkuCXtDtwGHBuFds3M0uZIqL8jUqXAH8J7ACcFhGHt1hmPjAfoFarzRscHBx3vSMjI/T09IxqG161thNd3mLVZsLq9aPb+ubMqqQvZf6sm+uuquYyNf9+V/W7XebPutXfdLfrZM0DAwPLI6Le3L51R9Y+AZIOB9ZExHJJ/WMtFxGLgcUA9Xo9+vvHXPQXhoaGaF7uhIVLp9DbLd+Cvg0sGh69G1cc019JX8r8WTfXXVXNZWr+/a7qd7vMn3Wrv+luV0bNVZzqOQj4kKQVwCDwPkkXVNAPM7MklR78EfHZiNg9InqBjwL/FhHHlt0PM7NU+Tp+M7PElH6Ov1FEDAFDVfbBzCw1PuI3M0uMg9/MLDEOfjOzxDj4zcwS4+A3M0uMg9/MLDEOfjOzxFR6Hb8Vo7fLxycys6nxEb+ZWWIc/GZmiXHwm5klxsFvZpYYB7+ZWWIc/GZmiXHwm5klxsFvZpYYB7+ZWWJKD35Je0i6XtJ9ku6VdGrZfTAzS1kVQzZsABZExB2SdgCWS7omIu6roC9mZskp/Yg/Ip6IiDvyx88B9wNzyu6HmVmqFBHVbVzqBW4E5kbEuqZ584H5ALVabd7g4OC46xsZGaGnp2dU2/CqtR3q7ZapNhNWr6+6F+VLsW7XXJ2+ObNK21arHJusgYGB5RFRb26vLPgl9QA3AF+MiMs2t2y9Xo9ly5aNu86hoSH6+/tHtXX7SJUL+jawaDi9QVZTrNs1V2fFWYeVtq1WOTZZkloGfyVX9UjaBrgUuHC80Dczs86q4qoeAd8E7o+Ir5W9fTOz1FVxxH8QcBzwPkl35l+HVtAPM7MklX7yLCJ+CKjs7ZqZWcaf3DUzS4yD38wsMQ5+M7PEOPjNzBLj4DczS4yD38wsMQ5+M7PEVD8IhpnZFq7MMb8W9G3ghIbtFTFOkI/4zcwS4+A3M0uMg9/MLDEOfjOzxDj4zcwS4+A3M0uMg9/MLDEOfjOzxDj4zcwSU9XN1g+R9ICkhyUtrKIPZmapquJm6zOAc4APAvsCR0vat+x+mJmlqooj/ncBD0fETyLiZWAQOLKCfpiZJUkRUe4GpaOAQyLiD/Lp44Bfj4hPNS03H5ifT+4DPNDG6mcDT3Wwu9NBijVDmnW75jR0sua3RMQuzY1b7OicEbEYWDyR50haFhH1grq0RUqxZkizbtechjJqruJUzypgj4bp3fM2MzMrQRXBfzuwt6Q9JW0LfBS4soJ+mJklqfRTPRGxQdKngH8FZgBLIuLeDq1+QqeGukSKNUOadbvmNBRec+lv7pqZWbX8yV0zs8Q4+M3MEtM1wZ/KMBCSVkgalnSnpGV5206SrpH0UP79DVX3cyokLZG0RtI9DW0ta1Tmr/P9frekA6rr+eSNUfOZklbl+/pOSYc2zPtsXvMDkn6nml5PjaQ9JF0v6T5J90o6NW/v2n29mZrL3dcRMe2/yN4kfgTYC9gWuAvYt+p+FVTrCmB2U9uXgYX544XAl6ru5xRrfA9wAHDPeDUChwL/Agg4ELi16v53sOYzgdNaLLtv/jv+OmDP/Hd/RtU1TKLmXYED8sc7AA/mtXXtvt5MzaXu62454k99GIgjgfPyx+cBH66wL1MWETcCzzQ1j1XjkcC3I3MLsKOkXcvpaeeMUfNYjgQGI+KliPgp8DDZ38C0EhFPRMQd+ePngPuBOXTxvt5MzWMpZF93S/DPAR5rmF7J5n+Y01kAV0tang9rAVCLiCfyx08CtWq6Vqixauz2ff+p/LTGkoZTeF1Xs6Re4NeAW0lkXzfVDCXu624J/pS8OyIOIBvd9JOS3tM4M7LXh119jW4KNea+AbwV2B94AlhUbXeKIakHuBT444hY1zivW/d1i5pL3dfdEvzJDAMREavy72uAy8le9q3e9JI3/76muh4WZqwau3bfR8TqiHg1IjYC/8BrL/G7pmZJ25AF4IURcVne3NX7ulXNZe/rbgn+JIaBkLS9pB02PQY+ANxDVuvx+WLHA9+rpoeFGqvGK4Hfz6/4OBBY23CaYFprOn/9EbJ9DVnNH5X0Okl7AnsDt5Xdv6mSJOCbwP0R8bWGWV27r8equfR9XfW73B18t/xQsnfIHwFOr7o/BdW4F9k7/HcB926qE9gZuA54CLgW2Knqvk6xzovIXu6+QnZO86SxaiS7wuOcfL8PA/Wq+9/Bms/Pa7o7D4BdG5Y/Pa/5AeCDVfd/kjW/m+w0zt3AnfnXod28rzdTc6n72kM2mJklpltO9ZiZWZsc/GZmiXHwm5klxsFvZpYYB7+ZWWIc/NY1JO3cMLrhk02jHd5U0DZ/TdI3x5i3QtLsDm5rUNLenVqfpcuXc1pXknQmMBIRXy14O98FvhARd7WYt4LsWvOnOrSt9wLHRsTJnVifpctH/JYESSP5935JN0j6nqSfSDpL0jGSblN2n4O35svtIulSSbfnXwe1WOcOwDs2hX7+iuPqfJz1c8k+cLRp2SvygfXu3TS4nqSPS/qrhmVOlvT1/BPaSyXdJekeSb+XL/LvwMGSSr9XtnUXB7+laD/gE8CvAMcBb4+IdwHnAp/Olzkb+HpEvBP47/m8ZnVe+2g9wBnADyPiV8nGUXpzw7yPR8S8/DmnSNoZuBg4Ih+7BeBEYAlwCPB4ROwXEXOBHwBENo7Lw3n/zSbNRw6WotsjH+NF0iPA1Xn7MDCQPz4Y2DcbWgWAX5LUExEjDevZFfhZw/R7gP8GEBFLJT3bMO8USR/JH+8B7B0Rt0j6N+BwSfcD20TEsKSXgEWSvgRcFRH/3rCeNcBuwPJJV2/Jc/Bbil5qeLyxYXojr/1NbAUcGBEvbmY964HtxtuYpH6yfyS/EREvSBpqeN65wOeAHwP/CBARDyq7reChwBckXRcRn8+X3y7frtmk+VSPWWtX89ppHyTt32KZ+4G3NUzfCHwsX/6DwKabacwCns1D/5fJbhsIQETcSvYK4GNkA7UhaTfghYi4APgK2S0ZN3k7o08vmU2Yj/jNWjsFOEfS3WR/JzeSvS/wCxHxY0mzJO0Q2W30/gK4SNK9wE3Ao/miPwA+kZ/OeQC4pWlbFwP7R8SmU0N9wFckbSQbrfOPACTVgPUR8WSHa7XE+HJOsymQ9CfAcxHR6s3fdtdxFdkbyde1sa11EdHycwNm7fKpHrOp+Qaj3zNom6QdJT1IdhS/2dDP/ZzXbkJuNmk+4jczS4yP+M3MEuPgNzNLjIPfzCwxDn4zs8Q4+M3MEvP/AbpoFAauso24AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_censored = data[data.Event == 0]\n",
    "df_uncensored = data[data.Event == 1]\n",
    "\n",
    "df_censored.Time.hist()\n",
    "plt.title(\"Censored\")\n",
    "plt.xlabel(\"Time (days)\")\n",
    "plt.ylabel(\"Frequency\")\n",
    "plt.show()\n",
    "\n",
    "df_uncensored.Time.hist()\n",
    "plt.title(\"Uncensored\")\n",
    "plt.xlabel(\"Time (days)\")\n",
    "plt.ylabel(\"Frequency\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "WmFDrzYrn-JA"
   },
   "source": [
    "<a name='4'></a>\n",
    "## 4. Survival Estimates\n",
    "\n",
    "We'll now try to estimate the survival function:\n",
    "\n",
    "$$\n",
    "S(t) = P(T > t)\n",
    "$$\n",
    "\n",
    "To illustrate the strengths of Kaplan Meier, we'll start with a naive estimator of the above survival function. To estimate this quantity, we'll divide the number of people who we know lived past time $t$ by the number of people who were not censored before $t$.\n",
    "\n",
    "Formally, let $i$ = 1, ..., $n$ be the cases, and let $t_i$ be the time when $i$ was censored or an event happened. Let $e_i= 1$ if an event was observed for $i$ and 0 otherwise. Then let $X_t = \\{i : T_i > t\\}$, and let $M_t = \\{i : e_i = 1 \\text{ or } T_i > t\\}$. The estimator you will compute will be:\n",
    "\n",
    "$$\n",
    "\\hat{S}(t) = \\frac{|X_t|}{|M_t|}\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name='Ex-2'></a>\n",
    "### Exercise 2\n",
    "Write a function to compute this estimate for arbitrary $t$ in the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 306
    },
    "colab_type": "code",
    "id": "qUoKpBHJjZM-",
    "outputId": "9477925d-ee39-4489-e8be-9af3780f9fc4"
   },
   "outputs": [],
   "source": [
    "# UNQ_C2 (UNIQUE CELL IDENTIFIER, DO NOT EDIT)\n",
    "def naive_estimator(t, df):\n",
    "    \"\"\"\n",
    "    Return naive estimate for S(t), the probability\n",
    "    of surviving past time t. Given by number\n",
    "    of cases who survived past time t divided by the\n",
    "    number of cases who weren't censored before time t.\n",
    "    \n",
    "    Args:\n",
    "        t (int): query time\n",
    "        df (dataframe): survival data. Has a Time column,\n",
    "                        which says how long until that case\n",
    "                        experienced an event or was censored,\n",
    "                        and an Event column, which is 1 if an event\n",
    "                        was observed and 0 otherwise.\n",
    "    Returns:\n",
    "        S_t (float): estimator for survival function evaluated at t.\n",
    "    \"\"\"\n",
    "    S_t = 0.0\n",
    "    \n",
    "    ### START CODE HERE ###\n",
    "    X = sum(df['Time'] > t)\n",
    "    M = sum( (df['Time'] > t) | (df['Event'] == 1) )\n",
    "    S_t = X / M\n",
    "    ### END CODE HERE ###\n",
    "    \n",
    "    return S_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test Cases\n",
      "Sample dataframe for testing code:\n",
      "   Time  Event\n",
      "0     5      0\n",
      "1    10      1\n",
      "2    15      0\n",
      "\n",
      "\n",
      "Test Case 1: S(3)\n",
      "Output: 1.0, Expected: 1.0\n",
      "\n",
      "Test Case 2: S(12)\n",
      "Output: 0.5, Expected: 0.5\n",
      "\n",
      "Test Case 3: S(20)\n",
      "Output: 0.0, Expected: 0.0\n",
      "\n",
      "Test case 4: S(5)\n",
      "Output: 0.5, Expected: 0.5\n"
     ]
    }
   ],
   "source": [
    "print(\"Test Cases\")\n",
    "\n",
    "sample_df = pd.DataFrame(columns = [\"Time\", \"Event\"])\n",
    "sample_df.Time = [5, 10, 15]\n",
    "sample_df.Event = [0, 1, 0]\n",
    "print(\"Sample dataframe for testing code:\")\n",
    "print(sample_df)\n",
    "print(\"\\n\")\n",
    "\n",
    "print(\"Test Case 1: S(3)\")\n",
    "print(\"Output: {}, Expected: {}\\n\".format(naive_estimator(3, sample_df), 1.0))\n",
    "\n",
    "print(\"Test Case 2: S(12)\")\n",
    "print(\"Output: {}, Expected: {}\\n\".format(naive_estimator(12, sample_df), 0.5))\n",
    "\n",
    "print(\"Test Case 3: S(20)\")\n",
    "print(\"Output: {}, Expected: {}\\n\".format(naive_estimator(20, sample_df), 0.0))\n",
    "\n",
    "# Test case 4\n",
    "sample_df = pd.DataFrame({'Time': [5,5,10],\n",
    "                          'Event': [0,1,0]\n",
    "                         })\n",
    "print(\"Test case 4: S(5)\")\n",
    "print(f\"Output: {naive_estimator(5, sample_df)}, Expected: 0.5\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "gKMsSOzfmGwD"
   },
   "source": [
    "In the next cell, we will plot the naive estimator using the real data up to the maximum time in the dataset. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 281
    },
    "colab_type": "code",
    "id": "nauOrVd9mNDs",
    "outputId": "fc884506-6758-4d22-ab04-6310145f9ed3"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd5xU9bnH8c+zHbawwNKL9KoiuIhYQUTFhprYk6g3kTQ1MRpLikFzE6OJydUbr5EYSzR2TURFFI0lFhBQuoBIkaUL0jv73D/mgMOy5eyys2d25/t+vea1c8qc+c5h2WfO75zz+5m7IyIiqSst6gAiIhItFQIRkRSnQiAikuJUCEREUpwKgYhIilMhEBFJcSoEUqfM7Gdm9kDUOQ5GbX0GMxttZo/VRqYq3ucvZvbLRL+P1F8qBFItZrbYzFabWW7cvO+Y2VthXu/uv3X37yQgV18ze83M1pnZejObaman1/b7QOI+QzwzG2JmpWa2ucxjcBWvu9zM3i2T93vu/usE5XQz65aIbUvdUSGQmkgHfhR1iDJeBCYArYGWwDXAxppsyMzSazHXwVju7nllHh9EHUoaHhUCqYnfA9ebWWF5C83sbjNbamYbg2/mx8ct29ccYmavmNlVZV473czOC573MrMJwbf8eWZ2QQXvVwR0Bv7q7juDx3vu/m6w/IBvyfHfZM3sYTO7z8zGmdmW4LOtjC8IZnaumc2owWeocF8cjOAzLTSzTWa2yMwuNbPewF+AwcHRw/q4z/ffwfMhZlZiZjcER3YrzOwcMzvdzOYH+/pnce9zlJl9EBxlrTCzP5tZVrDsnWC16cH7XRjMP9PMpgWved/MDq+NzyyJo0IgNTEFeAu4voLlk4EjgGbA48AzZpZTznpPABfvnTCzPsAhwMtB09OE4PUtgYuA/wvWKWstsAB4LPij1qoGn+kS4DdAPnA3sAU4qczyx6vzGYJZYfdFaMG+uQcY4e75wDHANHf/BPge8EFw9FBuoSZ21JQDtANuAf4KfAM4Ejge+KWZdQ7W3QNcCxQBg4FhwA8A3P2EYJ1+wfs9ZWb9gQeB7wLNgfuBsWaWfTCfWRJLhUBq6hbgajNrUXaBuz/m7mvdfbe73wVkAz3L2cY/gSPM7JBg+lLgeXffAZwJLHb3h4LtfAw8B5xfzvs5MBRYDNwFrDCzd8ysezU+zwvBUUSpu28n7g+8meUDpwfzqvMZqrMvytM2+FYd/9h7bqYUONTMGrn7CnefXY3Pugv4jbvvAp4k9kf+bnffFGxnDtAvyD/V3ScG+RcT+8N+YiXbHgXc7+6T3H2Puz8C7ACOrkY+qWMqBFIj7j4LeAm4qewyM7vezD4xsw1B80QTYn9sym5jE7FvzhcFsy4G/hE8PwQYFP9HkNgf2dYV5Clx96vcvWvw2i3A36vxkZaWmX4cOC/4Jnse8JG7L6nmZwi9Lyqw3N0Lyzy2uPsW4EJi3/5XmNnLZtarGp91rbvvCZ5vC36uilu+DcgL8vcws5eCprKNwG+ryH8IcF2Zf7cOQNtq5JM6pkIgB+NXwJXEmhgACNrAbwAuAJoGzRMbAKtgG08AFwdXw+QAbwbzlwJvl/kjmOfu368qlLsvBe4FDg1mbQEax2Usr5js1w2vu88BlgAjqLhZqNLPUIN9EZq7v+ruw4E2wFxizTsHfI5acF+w/e7uXgD8jMrzLyV2tBH/79bY3cs7mpIkoUIgNebuC4CniF2hs1c+sBtYA2SY2S1AQSWbGUfsW+RtwFPuXhrMfwnoYWbfNLPM4DEwOCG6HzNrama3mlk3M0sLTh7/FzAxWGU60NfMjgja50eH/IiPE7s66gTgmRp8hurui1DMrJWZjQyaiXYAm4k1FUHsm337vSd0a0E+sauvNgdHHWUL8SqgS9z0X4Hvmdkgi8k1szOC5jVJUioEcrBuA3Ljpl8FxgPziX2j3s6BzS77BG3pzwMnE/etO2hyOYVYk8tyYCVwB7E29rJ2Ap2A14n90ZpF7A/k5cG25gc5Xwc+Bd4tZxvleYJYe/i/3f2L6n4GqrkvytHWDryP4GvE/t/+hNh+WRdk3PsH+t/AbGClmVWYuRquJ3ZEtInYH/mnyiwfDTwSNANd4O5TiB0l/hn4kthJ/MtrIYckkGlgGhGR1KYjAhGRFKdCICKS4lQIRERSnAqBiEiKy4g6QHUVFRV5p06doo4hIlKvTJ069Qt3P6AnAKiHhaBTp05MmTIl6hgiIvWKmR1wZ/xeahoSEUlxKgQiIilOhUBEJMWpEIiIpDgVAhGRFJewQmBmDwZD4c2qYLmZ2T1mtsDMZpjZgERlERGRiiXyiOBh4LRKlo8AugePUcT6PRcRkTqWsPsI3P0dM+tUySojgb8HwwxONLNCM2vj7isSkWfy4nX8Z/6afdMn9W7FER0qGtJVRCR1RHlDWTv275u9JJh3QCEws1HEjhro2LFjjd7soyVf8r9vLgDAHV6bs4rxPz6hileJiDR89eJksbuPcfdidy9u0aLcO6Sr9N0Tu7Lo9jNYdPsZ3Hp2X+au3MSMkvW1nFREpP6JshAsIzao9V7tg3kJd+bhbTCDs//8Hi/NWF4XbykikrSiLARjgW8FVw8dDWxI1PmBsprnZTPmm8UAvP/Z2rp4SxGRpJWwcwRm9gQwBCgysxLgV0AmgLv/hdiA36cTG9N0K3BForKUZ3ifVhzVuRnzVm6qy7cVEUk6ibxq6OIqljvww0S9fxi9Wufz/EfLcHfMLMooIiKRqRcnixOlV+sCNu/YzcxlG1i9aTvbd+2JOpKISJ2rd+MR1KY+bQuA2EljgFYF2fznhpPIykjp+igiKSalC0G/9k245+L+bNy2i6VfbuX+txfy1rzVnNK3ddTRRETqTEoXAjPj7H5tAdi1p5RnppTwwrTlKgQiklLUBhLITE9jSI8WTF3yZdRRRETqlApBnPbNGrNq03Z27SmNOoqISJ1RIYjTrjAHd1i5YXvUUURE6kyoQmBmx5nZFcHzFmbWObGxotGusDEAy9ZviziJiEjdqbIQmNmvgBuBm4NZmcBjiQwVlbaFOQAsVyEQkRQS5ojgXOBsYAuAuy8H8hMZKiptCxsBKgQiklrCFIKdQXcQDmBmuYmNFJ2czHSa52YxYc4qtu3UXcYikhrCFIKnzex+oNDMrgReBx5IbKzodC7KZXrJBn41ttyhlkVEGpwqC4G7/wF4FngO6Anc4u73JDpYVP73kv4c3r4JL89YoaMCEUkJVd5ZbGZ3uPuNwIRy5jU4bZo04uYRvbn4rxO54P4PyM/JoFluFn84vx85melRxxMRqXVhmoaGlzNvRG0HSSaDOjfjvP7tyMlMY/OO3bw0YwUT5qyKOpaISEJUeERgZt8HfgB0MbMZcYvygfcSHSxKaWnGHy88AoA9pc4xv3uDP02YzxEdCunQrHHE6UREaldlRwSPA2cRG1LyrLjHke7+jTrIlhTS04xz+7dn4RdbOO++99X9hIg0OBUWAnff4O6L3f1id18CbCN2CWmemXWss4RJ4Nrh3blqaDfWbNrBuwu+iDqOiEitCnNn8Vlm9imwCHgbWAy8kuBcSSU7I51rhnWnSaNMxk5bHnUcEZFaFeZk8X8DRwPz3b0zMAyYmNBUSSgrI43juhcxZcm6qKOIiNSqMIVgl7uvBdLMLM3d3wSKE5wrKfVunc/SddvYvGN31FFERGpNmEKw3szygHeAf5jZ3QT9DqWanq1jYxzPW7kp4iQiIrUnTCEYCWwFrgXGA58Ru3oo5fRqHetrT4VARBqSSu8sNrN04CV3HwqUAo/USaok1b5pI/KyM5izYkPUUUREak2lRwTuvgcoNbMmdZQnqZkZR3dpxutzVrOn1KOOIyJSK6rsawjYDMw0swnEnRtw92sSliqJnX1EO17/ZDWTFq3lmK5FUccRETloYQrB88FDgOG9W5Gblc4LHy9XIRCRBqHKQuDuKX1eoKxGWemc2rc142at4NaRfdUjqYjUe6EGr5f9jezfjk3bd3PY6FdZuGZz1HFERA6KCkENHNetiIsGdmDXHuf9z9ZGHUdE5KCoENRAeppx+3mHkZ+dwdyVG6OOIyJyUCobj+BFggHry+PuZ1e1cTM7DbgbSAcecPfflVnekdi9CYXBOje5+7hw0aNlZvRsna+by0Sk3qvsZPEfDmbDwc1o9xIb4awEmGxmY919TtxqvwCedvf7zKwPMA7odDDvW5d6tcnnhWnLcXfMLOo4IiI1UmEhcPe3D3LbRwEL3H0hgJk9Say7ivhC4EBB8LwJUK/6eO7VuoDHtn9Oj1+8ghErBFcc24mbT+8dcTIRkfDCDF7fHbgd6APk7J3v7l2qeGk7YGncdAkwqMw6o4HXzOxqIBc4uYIMo4BRAB07Js+YOGcd3pY1m3awMxi1bOLCtfxj0udcO7yHLisVkXojzMnih4D7gN3AUODvwGO19P4XAw+7e3vgdOBRMzsgk7uPcfdidy9u0aJFLb31wWvSOJNrh/fgxtN6ceNpvbhueE8279jN659ooHsRqT/CFIJG7v4GYO6+xN1HA2eEeN0yoEPcdPtgXrxvA08DuPsHxI446u3tuoO7NqewcSbvzF8TdRQRkdDCFIIdwbf0T83sKjM7F8gL8brJQHcz62xmWcBFwNgy63xObMQzzKw3sUJQb/+KpqcZfdoU6EoiEalXwhSCHwGNgWuAI4FvAJdV9SJ33w1cBbwKfELs6qDZZnabme299PQ64Eozmw48AVzu7vW6W8+erfOZv2qzeicVkXojTKdze9x9M7FeSK+ozsaDewLGlZl3S9zzOcCx1dlmsuvduoBtu/bw+bqtdC7KjTqOiEiVwhwR3GVmn5jZr83s0IQnqud6tdk7ipnuOBaR+qHKQhCMTjaUWNv9/WY208x+kfBk9VSPVvnkZKbxn0+/iDqKiEgoofoacveV7n4P8D1gGnBLFS9JWTmZ6ZzSpzUvz1zBzt2lUccREalSlYXAzHqb2Wgzmwn8L/A+sUtBpQLn9G/L+q27dBmpiNQLYY4IHgTWA6e6+xB3v8/dVyc4V712fPcWNMvN4l/Tyt42ISKSfMKMUDa4LoI0JJnpaZxxWBuembqUzTt2k5cd5uIsEZFoVHhEYGZPBz9nmtmMuMdMM5tRdxHrp3P6t2X7rlJenbUy6igiIpWq7Kvqj4KfZ9ZFkIZmQMemdGjWiH9NW8bXjtQpFRFJXhUeEbj7iuDp14BdQT9D+x51E6/+MjNG9mvHewu+YPWm7VHHERGpUJiTxfnABDP7T9DXUKtEh2oozunfllKHF6evqHplEZGIhLmh7FZ37wv8EGgDvG1mryc8WQPQrWU+h7Yr4AVdPSQiSaw6g9evBlYCa4GWiYnT8JxzRDtmlGxg4ZrNUUcRESlXmBvKfmBmbwFvAM2BK9398EQHayjO6tcWM3j4/cXMKFnPrGUb2L1HdxyLSPIIc4F7e+DH7j4t0WEaolYFORzXrYi/f7CEv38QO8d+3fAeXD2se8TJRERiKj0iMLN04DwVgYPzpwuP4G+XFfO3y4o5okMhz31UQj0fdkFEGpBKC4G77wHmmVnyjBhfDxXlZTOsdyuG9W7FJUd1ZPHarUwv2RB1LBERINzJ4qbAbDN7w8zG7n0kOlhDddphrcnKSONfH+tKIhFJDmHOEfwy4SlSSEFOJsN6teSlGcv5xRm9yUivzoVbIiK1L8x9BG+X96iLcA3VeQPa88Xmndz64pyoo4iIhLp8dJOZbQwe281sj5lpHMaDcHLvlnxr8CE8OnEJn+n+AhGJWJgjgnx3L3D3AqARsb6H/i/hyRowM+OC4g4AzF2xKeI0IpLqqtVA7TH/Ak5NUJ6U0a1lHmmmQe5FJHpVniw2s/PiJtOAYkDdaR6knMx0OhXlMneljghEJFphrho6K+75bmAxMDIhaVJMr9b5zF6uIwIRiVaYoSqvqIsgqahf+0LGzVzJ3JUb6dW6IOo4IpKiwlw1dKeZFZhZZnBT2Roz+0ZdhGvoLhzYgdysdC796yTGzdSYBSISjTAni09x943EhqxcDHQDfprIUKmisHEW1w7vwdotO7nlhVnqlVREIhGmEOxtPjoDeMbd1UlOLfrO8V2479IBfLF5J5MWrYs6joikoDCF4CUzmwscCbxhZi3QVUO1amivluRmpfPSjOVRRxGRFBTmhrKbgGOAYnffBWxFVw3VqpzMdIb3acUrs1ayc7eah0SkboW6oczd1wVdUuPuW9x9ZZjXmdlpZjbPzBaY2U0VrHOBmc0xs9lm9nj46A3LmYe3Zf3WXfxq7Gzufv1THnl/MaWlGrNARBIvzH0ENRIManMvMBwoASab2Vh3nxO3TnfgZuBYd//SzFJ2LOTjexTRsVljnvjw833zOhflckKPFhGmEpFUkLBCABwFLHD3hQBm9iSxJqX4LjevBO519y8B3H11AvMkteyMdN7+6RDcYcfuUgb+5nVemrFchUBEEi5MFxMGXAp0cffbgtHKWrv7h1W8tB2wNG66BBhUZp0ewXu8B6QDo919fDkZRgGjADp2bLiDpZkZZtAoK51T+rTi2aklvDZn1b7lGWnG78/vx9CeKXvgJCIJEOaI4P+AUuAk4DZgE/AcMLCW3r87MARoD7xjZoe5+/r4ldx9DDAGoLi4OCUazq8e1p0mjTP3O0/w8swV/GPiEhUCEalVYQrBIHcfYGYfAwRt+VkhXrcM6BA33T6YF68EmBRcjbTIzOYTKwyTQ2y/QetclMuvzuq737ysjDQefn8xz04tISPNAOjdpoCerfOjiCgiDUSYQrArOPHrAMF9BGGucZwMdDezzsQKwEXAJWXW+RdwMfCQmRURaypaGDJ7yjmnfzseeHcR1z8zfd+8orxsPrj5JDI15KWI1FCYQnAP8E+gpZn9Bvg68IuqXuTuu83sKuBVYu3/D7r7bDO7DZji7mODZaeY2RxgD/BTd19bw8/S4PVt24RJNw9jy849AExd8iXXPzOdMe8sZNQJXVQMRKRGzL3qJncz6wUMAwx4w90/SXSwihQXF/uUKVOievuksqfUOeHON1m2fhuXH9OJ0Wf3rfpFIpKSzGyquxeXtyxM76P3AM3c/V53/3OURUD2l55mPPrtozikeWNenL5cndaJSI2EaRqaCvzCzHoSayJ60t31lTxJdGmRx80jevG9xz7iZ/+cSbPc7Nj8olwuGNihileLiIQbmOYR4BEza0Zs4Po7zKyju3dPeDoJZUjPlnQuyuVf02Kd1pWWOrtLnWO6Nad908YRpxORZFedO4u7Ab2AQwA1DyWRnMx03rx+yL7pJWu3cOLv3+KVmSu58oQu0QUTkXohzJ3FdwLnAp8BTwG/LnvDlySXQ5rncmi7Av7w2jweeDfc1bjZGek8cFkxPVrpngSRVBPmiOAzYLC7f5HoMFJ7fn56H16YVvb+vYo991EJT09eyi/O7JPAVCKSjCosBGbWy93nErsxrGPQx9A+7v5RosNJzQ3u2pzBXZuHXn/Nph2Mm7mCk3q1pF3TRhzSPDeB6UQkmVR2RPATYh293VXOMifW95A0EGcf0ZY35q7mkgcmUZCTwYc/P5mczPSoY4lIHaiwELj7qODpCHffb2hKM8tJaCqpc2cd3paOzRozc9kGbnlhNm/PX8OpfVtHHUtE6kCYPgneDzlP6rG0NKN/x6ZcfFRHmjbO5K7X5nH9M9OZuuTLqKOJSIJVWAjMrLWZHQk0MrP+ZjYgeAwBdHF6A5WZnsZ3ju/Clh17eHH6cu56bV7UkUQkwSo7R3AqcDmx7qP/GDd/E/CzBGaSiP1waDd+OLQbf3xtHn9+cwFfbN5BUV521LFEJEEqPCJw90fcfShwubsPjXuc7e7P12FGicjph7eh1GH8rJVRRxGRBArTxcRzZnYG0BfIiZt/WyKDSfR6tsqnS4tcxs1cwTeOPiTqOCKSIGF6H/0LcCFwNbFuqM8n1s2ENHBmxpmHtWHiwrV8sXlH1HFEJEHCXDV0jLt/C/jS3W8FBhMMOi8Nn5qHRBq+MIVgW/Bzq5m1BXYBbRIXSZLJ3uahl2esiDqKiCRImELwkpkVAr8HPgIWA08kMpQkj73NQ5MWreVbD37Ir1+aE3UkEallYU4W/zp4+pyZvQTkuPuGxMaSZHJ+cQcmLVpHybqtvDN/DRcO7KBeSkUakMo6nTuvkmXoEtLU0aFZY5767mBWb9rOoN++wcszVtBjuAqBSENR2RHBWZUsc0CFIMW0zM9hUOdmvDxzBT8+uTtmFnUkEakFlXU6d0VdBpH64YzD2vDLF2Yzf9VmerbWUYFIQxBmhLJbypuvG8pS06mHtuaWsbP5w2vzGNS5WYXrDenZkm4t8+owmYjUVJgRyrbEPc8BzkRjFqeslvk5DOvVkglzVjFhzqoK13ttziqe/u7gOkwmIjUV5qqh/QamMbM/AK8mLJEkvTHfLGbLzt0VLr//7YXc+9YCVm/cTssCDV0hkuzCHBGU1ZhYj6SSotLSjPyczAqXjzyiLX9+cwFf+8v75Gd/td6h7Qq48+v96iKiiFRDmHMEM4ldJQSQDrQAdH5AKtS9VT7fO7ErC1Zv3jdvzabtPD2lhFEndNW5A5EkE+aI4My457uBVe5ecbuACHDTiF77Ta/csJ2jb3+DJz78nPOLq3dA2So/h6a5WbUZT0TihDlHsMTMmgIdgvVbBTeUfZTwdNJgtG6Sw1GdmvG3dxfxt3cXVeu1+TkZvPrjE2hb2ChB6URSW5imoV8TG6nsM75qInLgpMTFkobo7ouPYNrn66v1mh27S7nxuRl855Ep9Gh1YJNS28JG/PTUnrq5TeQghGkaugDo6u47Ex1GGrY2TRrR5rDqf6vftmsP97/9GR8v3b+IbN+1h1Ubd3Bq39b061BYWzFFUk6YQjALKARWV3fjZnYacDexk8wPuPvvKljva8CzwEB3n1Ld95GG7eKjOnLxUR0PmL9+606K//t1xs1aoUIgchDCFILbgY/NbBawb5gqdz+7sheZWTpwLzAcKAEmm9lYd59TZr184EfApGpmlxRX2DiLY7sVMeadhTz07mKuOLYTN5/eO+pYIvVOmELwCHAHMBMorca2jwIWuPtCADN7EhgJlO3Q/tfB9n9ajW2LAPCz03vTd1oB73+2lic+/JzrTulJVkaYYTZEZK8whWCru99Tg223A5bGTZcAg+JXMLMBQAd3f9nMKiwEZjYKGAXQseOBTQSSunq2zueG03rxxier+PYjU/jbu4voXkv3KfRsnU+HZo1rZVsiySxMIfiPmd0OjGX/pqGDunzUzNKAPxK7IqlS7j4GGANQXFzsVawuKei47kU0bZzJHePn1to287IzeOunQyjKy661bYokozCFoH/w8+i4eWEuH11G7N6DvdoH8/bKBw4F3gou/WsNjDWzs3XCWKorOyOdV398Aqs27qh65RDWbtnBtx+Zwn89PJlDmueWu86xXZtzUTknsUXqmzA3lA2t4bYnA93NrDOxAnARcEncdjcARXunzewt4HoVAamplgU5tdrJ3U2n9eKJDz9n9rIDR2Zdu2Un734aG7ZT9zBIfZew8QjcfbeZXUWsp9J04EF3n21mtwFT3H1sTQKL1JUrT+jClSd0KXfZox8s5pcvzGb5hu200x3PUs8ldDwCdx8HjCszr6LCMiTMNkWSwaHtmgAws2SDCoHUexqPQKQGercpID3NeGry5wzp2YKczPSoI4nUWE0uuNZ4BJLycjLT6d0mnzfnreG2l8reGiNSv1RZCMxsppnNCB6zgXnA/yQ+mkhyG/PNYnq0ymP8rJXs3lOdey1FkkuYI4IzgbOCxylAW3f/c0JTidQDbQsbce3JPVi3ZSc/eXo6m3domA6pn8IUgjbAOndf4u7LgEZmNqiqF4mkghN7tqBZbhZjpy/n8UlLoo4jUiNhCsF9wOa46S3BPJGU1zgrg8k/P5lerfN5dfaqqOOI1EiYQmDuvq9bB3cvpWaD3os0SOlpxpmHt2Hqki/57qNT2LF7T9SRRKolTCFYaGbXmFlm8PgRsDDRwUTqk/MGtKddYSNenb2K9xZ8EXUckWoJUwi+BxxDrJuIvT2IjkpkKJH6pm1hI968fgj52RmMn7Uy6jgi1RLmhrLVxPoJEpFKZGWkMax3S56eUsLzHy0rd50OzRoz/sfHk52hG9AkeaitX6QWXTu8Bx2aNabUD+wtffXGHTwztYT3F6xlaK+WEaQTKZ8KgUgtOqR5Lted0rPcZTt272H8rJU8NXkpeTnV/6/XuiBHA+VIQqgQiNSR7Ix0hvdpxfMfL2P87OqfR8jNSufDn59Mbrb+20rtqvA3ysx+UtkL3f2PtR9HpGG7dWRfzhtQ/a66Fn6xmVtemM3b89dw+mFtEpBMUlllXy3yg589gYHEhqqEWFcTHyYylEhDlZ+TyXHdi6pesYzBXZvzP69/yph3FrJg9Vf3dzbKTOdbxxyik89yUCosBO5+K4CZvQMMcPdNwfRo4OU6SSciQOymtfOL23P/2wuZtnT9fsua5mbx9SPVIbDUXJjGxlbAzrjpncE8EalDN4/ozQ2n9to37e4cf+ebvDp7pQqBHJQwheDvwIdm9s9g+hzgkcRFEpGKpKfFj49snNq3NY9P+pyLx0wEYFjvlnzn+PKH1xSpSJV3Frv7b4ArgC+DxxXu/ttEBxORql06qCMDOzdlT6mzfMM27hw/jy8274g6ltQz5uXc+HLASmbHAd3d/SEzawHkufuihKcrR3FxsU+ZMiWKtxZJap+t2cywu96mZ6t8WuRnV7heWppx/Sk9OLx9YR2mk6iZ2VR3Ly5vWZgRyn4F3AjcHMzKBB6rvXgiUhu6tsjjB0O6kpeTwbZdeyp8fLhoLQ+9tzjquJJEwpwjOBfoD3wE4O7LzSy/8peISBRuOK1Xletc9/R0JsxZya49pWSm12TYcmlowhSCne7uZuYAZpab4EwikkCn9m3Fcx+VcPjo19h77nlw1yIeuKzcVgNJAWEKwdNmdj9QaGZXAv8FPJDYWCKSKEN7teQnw3uwcdsuABav3crrn6xiRsl6nTdIUWFPFg8nNnC9Aa+6+4REB6uIThaL1K6N284vXOYAAA+eSURBVHdx9G/foFFmOoWNMw9YPqBjU35/fr8IkkltquxkcZVHBGZ2h7vfCEwoZ56I1HMFOZncenZf3pq/5oBlqzZs55mpJVxxbGf6tC2IIJ3UhSqPCMzsI3cfUGbeDHc/PKHJKqAjApG6s37rTgb99g06F+XSuajy04NpacY1J3WnZ2tdS5KManREYGbfB34AdDGzGXGL8oH3ajeiiCSjwsZZXDW0Gy/OWM5nazZXuu6StVsBuPeSAZWuJ8mnwiMCM2sCNAVuB26KW7TJ3dfVQbZy6YhAJDnd9uIcHp24mBtO7UVaXFcYx3RtTu82alaKWo2OCNx9A7ABuDjYSEsgB8gzszx3/zwRYUWkfrpkUEcem7iE34z7ZL/5fdoUMO5Hx0eUSsIIc7L4LOCPQFtgNXAI8AnQN8RrTwPuBtKBB9z9d2WW/wT4DrAbWAP8l7svqeZnEJEk0K1lHjNGn8KO3aX75v39/cXcNWE+y9Zvo11howjTSWXC3Fb438DRwHx37wwMAyZW9SIzSwfuBUYAfYCLzaxPmdU+BoqDE8/PAndWI7uIJJmczHSaNMrc9xgRjKb27Ycns3Td1ojTSUXCFIJd7r4WSDOzNHd/EwhzC+JRwAJ3X+juO4EngZHxK7j7m+6+97djIqBO1UUakK4tchnUuRlzV27iTxPmRx1HKhCmEKw3szzgHeAfZnY3sCXE69oBS+OmS4J5Ffk28Ep5C8xslJlNMbMpa9YceK2ziCQnM+Op7w7mkkEdGTdrBdOWrmf28g3MXr6BdVt2Vr0BqRNhupgYCWwHrgUuBZoAt9VmCDP7BrGjjBPLW+7uY4AxELtqqDbfW0QS7/wj2/P4pM85596vrjxvV9iId24YWmawHYlClYXA3bcAmFkB8GI1tr0M6BA33T6Ytx8zOxn4OXCiu2tEDZEGqH/Hpjw16mjWB/0bzV6+kXve+JQPPlvLcd2LIk4nYa4a+i5wK7GjglJi/Q05UNV4eJOB7mbWmVgBuAi4pMy2+wP3A6e5++pqpxeRemNQl+b7np/YowUPvbeInz47nTZNcip8zYUDO3DhwI51ES+lhTlHcD1wqLt3cvcu7t7Z3ascFNXddwNXAa8Su9z0aXefbWa3mdnZwWq/B/KAZ8xsmpmNreHnEJF6JCcznZtG9KJbyzxyszPKfazcsJ0/TpjPnlK1BidamL6GxgPnxV3dEyndWSySGl6cvpyrn/iYW87sQ+cW4YZBadukkfo6qsBB9T5KbIjK981sErCvDd/dr6mlfCIiBzi5dysKcjK47aU5oV+TlZHG5J+dTJNyutOWioUpBPcD/wZmEjtHICKScI2y0nn12hNYuWF7qPWXfrmNa574mFdnr+SCgR2qfoHsE6Zp6GN3719HeaqkpiERKY+7M+QPb7FjVyk94pqHcjLS+M25h9EiPzvCdNGrrGkozMniV4IbutqYWbO9j1rOKCJyUMyMH5/cndZNcti4bRcbt+1iw7ZdvDZnFc9/VBJ1vKQW5ohgUTmzPcyVQ4mgIwIRqY6z/vdd0tKMF354bNRRInVQRwTB5aJlH5EUARGR6hpxWGumL11PyZdJceFjUqqwEJjZScHP88p71F1EEZGaG3ForAfU8bNWRpwkeVV21dCJxK4WOqucZQ48n5BEIiK1qHNRLr3bFPDcR8toVVDxXcwAvdvk061l6t2HUNkIZb8Knt7m7vudJwi6jRARqRdGHtGW370yl6uf+LjS9doVNuLdG4dillod4YW5j+A5oOxo1M8CR9Z+HBGR2jfq+C6c0qcVpZVcHPPGJ6u5/ZW5zF6+kUPbNanDdNGrsBCYWS9iw1E2KXNOoIDY2MUiIvVCWprRpUVepes0bZzFHePn8oN/fMQDlxXTo1XqNBFVdtVQT+BMoJDYeYK9jwHAlYmPJiJSd5rnZTPi0DZ8vm4rd7wyN+o4daqycwQvAC+Y2WB3/6AOM4mIROLeSwdQ9MIsnpy8lG0799AoKz3qSHUizDmCc81sNrANGA8cDlzr7o8lNJmISASG9W7FIx8s4ZonP6YoL2vf/AEdm3J+ccPswyhMITjF3W8ws3OBxcB5xMYvViEQkQZnUJdm9O9YyPSl6/fN27pzD//8eBln9WtLTmbDO0oIUwj29ud6BvCMu29ItUurRCR1ZGek888f7N8dxVvzVnP5Q5P5YOFahvZsGVGyxAlTCF40s7nEmoa+b2YtiA1bKSKSEo7u0pyczDT+Z8J83vhk1b75TRtn8aNh3clID9N/Z/IKM3j9TWZ2J7DB3feY2VZgZOKjiYgkh5zMdC4a2JEXpy+n5MttAOwudTZs28XRXZpzbLeiiBMenMruI7jB3e8MJoe5+zMA7r7FzH4O/KwuAoqIJIPRZ/dl9Nl9901v2bGb/rdN4O35axpuIQAuAvYWgpuBZ+KWnYYKgYiksNzsDAZ2bsrLM1bQtHFW1S+oQL/2TTgm4kJSWSGwCp6XNy0iknLOOrwtNz0/kzvG1/wGtMLGmUz9xXDS06L7s1pZIfAKnpc3LSKSci46qiPn9G9X49e/PGMF1z0znVnLNtCvQ2EtJqueygpBPzPbSOzbf6PgOcG0+hoSEYGDuq9gSM8WADw7tYQdu0v3zS9snFmnfR1V1sVEw7trQkQkiTTPy6Zfh0IenbiERycu2TffDN68bgidinLrJEeY+whERCRB/vrNI/l09eZ90+u27OTqJz7mP5+uUSEQEUkFLQtyaBk3cpq787tX5vLegrV8c3CnOsmgQiAikkTMjGO6Nuf5j5cx6LevA3Dl8V34zvFdEvae9fu+aBGRBmjUCV24oLgDQ3u2JCcznaenLE3o++mIQEQkyXRvlc/t5x0GwP+9tYA7x8/ji807KMrLTsj76YhARCSJHd2lOQAfLlqXsPdQIRARSWKHtWtC46x0Ji5cm7D3SGghMLPTzGyemS0ws5vKWZ5tZk8FyyeZWadE5hERqW8y09Mo7tSsfhYCM0sH7gVGAH2Ai82sT5nVvg186e7dgD8BdyQqj4hIfTWoczPmr9rM2s07ErL9RB4RHAUscPeF7r4TeJIDxzEYCTwSPH8WGGYa/kxEZD+JPk+QyELQDoi/5qkkmFfuOu6+G9gANC+7ITMbZWZTzGzKmjVrEhRXRCQ5Hd6+CSf1aknj7MRc6FkvLh919zHAGIDi4mL1fCoiKSUzPY0HLx+YsO0n8ohgGdAhbrp9MK/cdcwsA2gCJO6MiIiIHCCRhWAy0N3MOptZFrERz8aWWWcscFnw/OvAv91d3/hFROpQwpqG3H23mV0FvAqkAw+6+2wzuw2Y4u5jgb8Bj5rZAmAdsWIhIiJ1KKHnCNx9HDCuzLxb4p5vB85PZAYREamc7iwWEUlxKgQiIilOhUBEJMWpEIiIpDirb1drmtkaYEmVK5avCPiiFuMkmvImlvImlvImVnXzHuLuLcpbUO8KwcEwsynuXhx1jrCUN7GUN7GUN7FqM6+ahkREUpwKgYhIiku1QjAm6gDVpLyJpbyJpbyJVWt5U+ocgYiIHCjVjghERKQMFQIRkRSXMoXAzE4zs3lmtsDMboo6T3nMbLGZzTSzaWY2JZjXzMwmmNmnwc+mEeZ70MxWm9msuHnl5rOYe4L9PcPMBiRJ3tFmtizYx9PM7PS4ZTcHeeeZ2akR5O1gZm+a2Rwzm21mPwrmJ90+riRrMu/fHDP70MymB5lvDeZ3NrNJQbangm7zMbPsYHpBsLxTEmR92MwWxe3fI4L5B/e74O4N/kGsG+zPgC5AFjAd6BN1rnJyLgaKysy7E7gpeH4TcEeE+U4ABgCzqsoHnA68AhhwNDApSfKOBq4vZ90+we9FNtA5+H1Jr+O8bYABwfN8YH6QK+n2cSVZk3n/GpAXPM8EJgX77WngomD+X4DvB89/APwleH4R8FQSZH0Y+Ho56x/U70KqHBEcBSxw94XuvhN4EhgZcaawRgKPBM8fAc6JKoi7v0Ns3Ih4FeUbCfzdYyYChWbWpm6SxlSQtyIjgSfdfYe7LwIWEPu9qTPuvsLdPwqebwI+ITaud9Lt40qyViQZ9q+7++ZgMjN4OHAS8Gwwv+z+3bvfnwWGmZlFnLUiB/W7kCqFoB2wNG66hMp/aaPiwGtmNtXMRgXzWrn7iuD5SqBVNNEqVFG+ZN7nVwWHzw/GNbUlVd6gGaI/sW+CSb2Py2SFJN6/ZpZuZtOA1cAEYkcm6919dzm59mUOlm8AmkeV1d337t/fBPv3T2aWXTZroFr7N1UKQX1xnLsPAEYAPzSzE+IXeuwYMGmv9032fIH7gK7AEcAK4K5o4xzIzPKA54Afu/vG+GXJto/LyZrU+9fd97j7EcTGUD8K6BVxpAqVzWpmhwI3E8s8EGgG3Fgb75UqhWAZ0CFuun0wL6m4+7Lg52rgn8R+UVftPcQLfq6OLmG5KsqXlPvc3VcF/8FKgb/yVfNEUuQ1s0xif1j/4e7PB7OTch+XlzXZ9+9e7r4eeBMYTKwZZe9ojfG59mUOljcB1tZx1PispwVNcu7uO4CHqKX9myqFYDLQPbg6IIvYiZ+xEWfaj5nlmln+3ufAKcAsYjkvC1a7DHghmoQVqijfWOBbwdUMRwMb4po3IlOm3fRcYvsYYnkvCq4U6Qx0Bz6s42xGbBzvT9z9j3GLkm4fV5Q1yfdvCzMrDJ43AoYTO7fxJvD1YLWy+3fvfv868O/giCyqrHPjvhAYsXMZ8fu35r8LdXUWPOoHsbPq84m1Cf486jzl5OtC7KqK6cDsvRmJtUm+AXwKvA40izDjE8QO93cRa4P8dkX5iF29cG+wv2cCxUmS99Egz4zgP0+buPV/HuSdB4yIIO9xxJp9ZgDTgsfpybiPK8mazPv3cODjINss4JZgfhdiRWkB8AyQHczPCaYXBMu7JEHWfwf7dxbwGF9dWXRQvwvqYkJEJMWlStOQiIhUQIVARCTFqRCIiKQ4FQIRkRSnQiAikuIyql5FJDWZ2d7LNgFaA3uANcH0Vnc/JpJgIrVMl4+KhGBmo4HN7v6HqLOI1DY1DYnUgJltDn4OMbO3zewFM1toZr8zs0uDvuRnmlnXYL0WZvacmU0OHsdG+wlEvqJCIHLw+gHfA3oD3wR6uPtRwAPA1cE6dwN/cveBwNeCZSJJQecIRA7eZA/6dTGzz4DXgvkzgaHB85OBPnHd2ReYWZ5/1ee8SGRUCEQO3o6456Vx06V89X8sDTja3bfXZTCRMNQ0JFI3XuOrZiL2jjUrkgxUCETqxjVAcTCy1Bxi5xREkoIuHxURSXE6IhARSXEqBCIiKU6FQEQkxakQiIikOBUCEZEUp0IgIpLiVAhERFLc/wMz0dQ4y6h0mAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "max_time = data.Time.max()\n",
    "x = range(0, max_time+1)\n",
    "y = np.zeros(len(x))\n",
    "for i, t in enumerate(x):\n",
    "    y[i] = naive_estimator(t, data)\n",
    "    \n",
    "plt.plot(x, y)\n",
    "plt.title(\"Naive Survival Estimate\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Estimated cumulative survival rate\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "jg4VTizxqgKM"
   },
   "source": [
    "<a name='Ex-3'></a>\n",
    "### Exercise 3\n",
    "\n",
    "Next let's compare this with the Kaplan Meier estimate. In the cell below, write a function that computes the Kaplan Meier estimate of $S(t)$ at every distinct time in the dataset. \n",
    "\n",
    "Recall the Kaplan-Meier estimate:\n",
    "\n",
    "$$\n",
    "S(t) = \\prod_{t_i \\leq t} (1 - \\frac{d_i}{n_i})\n",
    "$$\n",
    "\n",
    "where $t_i$ are the events observed in the dataset and $d_i$ is the number of deaths at time $t_i$ and $n_i$ is the number of people who we know have survived up to time $t_i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>    \n",
    "<summary>\n",
    "    <font size=\"3\" color=\"darkgreen\"><b>Hints</b></font>\n",
    "</summary>\n",
    "<p>\n",
    "<ul>\n",
    "    <li>Try sorting by Time.</li>\n",
    "    <li>Use <a href=\"https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.Series.unique.html\">pandas.Series.unique<a> </li>\n",
    "    <li>If you get a division by zero error, please double-check how you calculated `n_t`</li>\n",
    "</ul>\n",
    "</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 527
    },
    "colab_type": "code",
    "id": "jnwysrz7CzNG",
    "outputId": "a26f0a84-98c2-4fe2-c40e-b18402f8bf20"
   },
   "outputs": [],
   "source": [
    "# UNQ_C3 (UNIQUE CELL IDENTIFIER, DO NOT EDIT)\n",
    "def HomemadeKM(df):\n",
    "    \"\"\"\n",
    "    Return KM estimate evaluated at every distinct\n",
    "    time (event or censored) recorded in the dataset.\n",
    "    Event times and probabilities should begin with\n",
    "    time 0 and probability 1.\n",
    "    \n",
    "    Example:\n",
    "    \n",
    "    input: \n",
    "    \n",
    "         Time  Censor\n",
    "    0     5       0\n",
    "    1    10       1\n",
    "    2    15       0\n",
    "    \n",
    "    correct output: \n",
    "    \n",
    "    event_times: [0, 5, 10, 15]\n",
    "    S: [1.0, 1.0, 0.5, 0.5]\n",
    "    \n",
    "    Args:\n",
    "        df (dataframe): dataframe which has columns for Time\n",
    "                          and Event, defined as usual.\n",
    "                          \n",
    "    Returns:\n",
    "        event_times (list of ints): array of unique event times\n",
    "                                      (begins with 0).\n",
    "        S (list of floats): array of survival probabilites, so that\n",
    "                            S[i] = P(T > event_times[i]). This \n",
    "                            begins with 1.0 (since no one dies at time\n",
    "                            0).\n",
    "    \"\"\"\n",
    "    # individuals are considered to have survival probability 1\n",
    "    # at time 0\n",
    "    event_times = [0]\n",
    "    p = 1.0\n",
    "    S = [p]\n",
    "    \n",
    "    ### START CODE HERE (REPLACE INSTANCES OF 'None' with your code) ###\n",
    "    \n",
    "    # get collection of unique observed event times\n",
    "    observed_event_times = df.Time.unique()\n",
    "  \n",
    "    # sort event times\n",
    "    observed_event_times = sorted(observed_event_times)\n",
    "    \n",
    "    # iterate through event times\n",
    "    for t in observed_event_times:\n",
    "  \n",
    "        # compute n_t, number of people who survive to time t\n",
    "        n_t = len(df[df.Time >= t])\n",
    "  \n",
    "        # compute d_t, number of people who die at time t\n",
    "        d_t = len(df[(df.Time == t) & (df.Event == 1)])\n",
    "        \n",
    "        # update p\n",
    "        p = p*(1 - float(d_t)/n_t)\n",
    "  \n",
    "        # update S and event_times (ADD code below)\n",
    "        # hint: use append\n",
    "        event_times.append(t)\n",
    "        S.append(p)\n",
    "    \n",
    "    ### END CODE HERE ###\n",
    "  \n",
    "    return event_times, S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TEST CASES:\n",
      "\n",
      "Test Case 1\n",
      "\n",
      "Test DataFrame:\n",
      "   Time  Event\n",
      "0     5      0\n",
      "1    10      1\n",
      "2    15      0\n",
      "\n",
      "Output:\n",
      "Event times: [0, 5, 10, 15], Survival Probabilities: [1.0, 1.0, 0.5, 0.5]\n",
      "\n",
      "Expected:\n",
      "Event times: [0, 5, 10, 15], Survival Probabilities: [1.0, 1.0, 0.5, 0.5]\n",
      "\n",
      "Test Case 2\n",
      "\n",
      "Test DataFrame:\n",
      "   Time  Event\n",
      "0     2      0\n",
      "1    15      0\n",
      "2    12      1\n",
      "3    10      1\n",
      "4    20      1\n",
      "\n",
      "Output:\n",
      "Event times: [0, 2, 10, 12, 15, 20], Survival Probabilities: [1.0, 1.0, 0.75, 0.5, 0.5, 0.0]\n",
      "\n",
      "Expected:\n",
      "Event times: [0, 2, 10, 12, 15, 20], Survival Probabilities: [1.0, 1.0, 0.75, 0.5, 0.5, 0.0]\n"
     ]
    }
   ],
   "source": [
    "print(\"TEST CASES:\\n\")\n",
    "\n",
    "\n",
    "print(\"Test Case 1\\n\")\n",
    "\n",
    "print(\"Test DataFrame:\")\n",
    "sample_df = pd.DataFrame(columns = [\"Time\", \"Event\"])\n",
    "sample_df.Time = [5, 10, 15]\n",
    "sample_df.Event = [0, 1, 0]\n",
    "print(sample_df.head())\n",
    "print(\"\\nOutput:\")\n",
    "x, y = HomemadeKM(sample_df)\n",
    "print(\"Event times: {}, Survival Probabilities: {}\".format(x, y))\n",
    "print(\"\\nExpected:\")\n",
    "print(\"Event times: [0, 5, 10, 15], Survival Probabilities: [1.0, 1.0, 0.5, 0.5]\")\n",
    "\n",
    "print(\"\\nTest Case 2\\n\")\n",
    "\n",
    "print(\"Test DataFrame:\")\n",
    "\n",
    "sample_df = pd.DataFrame(columns = [\"Time\", \"Event\"])\n",
    "sample_df.loc[:, \"Time\"] = [2, 15, 12, 10, 20]\n",
    "sample_df.loc[:, \"Event\"] = [0, 0, 1, 1, 1]\n",
    "print(sample_df.head())\n",
    "print(\"\\nOutput:\")\n",
    "x, y = HomemadeKM(sample_df)\n",
    "print(\"Event times: {}, Survival Probabilities: {}\".format(x, y))\n",
    "print(\"\\nExpected:\")\n",
    "print(\"Event times: [0, 2, 10, 12, 15, 20], Survival Probabilities: [1.0, 1.0, 0.75, 0.5, 0.5, 0.0]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "G7OAiWjS7hLA"
   },
   "source": [
    "Now let's plot the two against each other on the data to see the difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 265
    },
    "colab_type": "code",
    "id": "JbPlC5717gM_",
    "outputId": "06553ed7-9396-4f16-eab3-43a1c6ecddfd"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3yV9fn/8deVk00SVtgzIMieCaCAgjjAAVZRHF8HteKviqNiLVpt0dbZql+1fqu4q60gaisq4oYWKSPIkiGyCRCm7JVx/f44N+EQck7ujJP7JOd6Ph555NzrnHduI1fuz+dzf25RVYwxxkSvGK8DGGOM8ZYVAmOMiXJWCIwxJspZITDGmChnhcAYY6JcrNcByio9PV1bt27tdQxjjKlWFixYsFNVG5S0rdoVgtatW5Odne11DGOMqVZEZEOwbdY0ZIwxUc4KgTHGRDkrBMYYE+WqXR+BMSZy5OXlkZOTw5EjR7yOYhyJiYk0b96cuLg418dYITDGlFtOTg6pqam0bt0aEfE6TtRTVXbt2kVOTg4ZGRmujwtb05CIvCYi20Xk+yDbRUSeE5HVIrJERHqFK4sxJjyOHDlC/fr1rQhECBGhfv36Zb5CC2cfwRvA0BDbhwHtnK8xwF/DmMUYEyZWBCJLef57hK1pSFX/LSKtQ+wyAvib+ufBniMidUSkiapuDUee3Ml3sWXPYWZk3A3AOR0b0aNFnXB8lDHGVCtejhpqBmwKWM5x1p1CRMaISLaIZO/YsaNcH1a4dQnHNi/m+W9W89zXqxn//pJyvY8xJrKICOPGjSta/vOf/8yECRNCHjN16lQef/zxMCerPqrF8FFVnaiqmaqa2aBBiXdIl6pp7ST6ZdRn3WMX8dDwzqzM3c+SnD2VnNQYU9USEhL44IMP2Llzp+tjhg8fzvjx48OYqnrxshBsBloELDd31oVP7lJ4/SKuWf5LJsX/gXf++jAfL9kS1o80xoRXbGwsY8aM4Zlnnjll20cffUTfvn3p2bMn5557Ltu2bQPgjTfeYOzYsezdu5dWrVpRWFgIwMGDB2nRogV5eXmsWbOGoUOH0rt3bwYOHMjKlSur9OeqSl4OH50KjBWRSUBfYG+4+gcA6Dqy6GWcL4beCTlwdDZT1+zi4m5Nw/axxkSLhz5axvIt+yr1PTs1TeP3l3Qudb/bbruNbt26ce+99560fsCAAcyZMwcR4ZVXXuHJJ5/kqaeeKtpeu3ZtevTowcyZMxk8eDAff/wxF1xwAXFxcYwZM4YXX3yRdu3aMXfuXG699Va+/vrrSv35IkXYCoGIvAMMAtJFJAf4PRAHoKovAtOAC4HVwCFgdLiyAJA52v/liHv9IlK37uWH3P1h/VhjTPilpaVx/fXX89xzz5GUlFS0Picnh1GjRrF161aOHTtW4tj6UaNGMXnyZAYPHsykSZO49dZbOXDgALNnz+aKK64o2u/o0aNV8rN4IZyjhq4uZbsCt4Xr891Ijo/lh9z9qKoNgTOmgtz85R5Od911F7169WL06BN/8N1+++3cfffdDB8+nBkzZpTYiTx8+HDuv/9+du/ezYIFCzjnnHM4ePAgderUYdGiRVX4E3inWnQWh0tyvI8DR/NZunkv2/cf4UhegdeRjDHlVK9ePa688kpeffXVonV79+6lWTP/YMQ333yzxONSUlLIysrizjvv5OKLL8bn85GWlkZGRgZTpkwB/HfsLl68OPw/hEeiuhCkH1hV1Gnc55GvOPtP33Asv9DrWMaYcho3btxJo4cmTJjAFVdcQe/evUlPTw963KhRo3j77bcZNWpU0bq///3vvPrqq3Tv3p3OnTvz4YcfhjW7l8TfQlN9ZGZmaqU8mCb7dXTpFPK3LGF3yum81v4vvDRzLROv6835nRtX/P2NiQIrVqygY8eOXscwxZT030VEFqhqZkn7R+8VQeZoZPQ04pp2p9GhH/lN7jjeS/wjO2e+5HUyY4ypUtFbCI7rOhIadyUGoXPMRjrs+MzrRMYYU6VsGuqAYaU7nx7MsT2HySsoJM5nNdIYEx3sX7sACbH+05G71x6yYYyJHlYIAiTE+gDYvOewx0mMMabqWCEIEO9cEWyxQmCMiSJWCAIkWCEwptpJSUkpej1t2jTat2/Phg0byvw+69evp0uXLpWSaf369YgIDzzwQNG6nTt3EhcXx9ixY0Me68UU2VYIAsSIEBsjfLF8G4eP2V3GxlQnX331FXfccQeffvoprVq18joOGRkZfPLJJ0XLU6ZMoXPn0qfhKOsU2fn5+eXKF8gKQTGZLKfT1g/4/dQSH7VsjIlA//73v7n55pv5+OOPadu2LRB8CuoJEyZw3XXXccYZZ9CuXTtefvnlU95v/fr1DBw4kF69etGrVy9mz54NwIwZMxg0aBAjR46kQ4cOXHvttQS7KTc5OZmOHTty/AbYyZMnc+WVVxZt37FjB5dffjlZWVlkZWXx7bffAiemyA61z/GfoX///lx33XUVPn82fDRQ15GwYRbXJM3jqiXn89DwLiTF+7xOZUz18Ol4/zM/KlPjrjAsdDPJ0aNHufTSS5kxYwYdOnQoWh9qCuolS5YwZ84cDh48SM+ePbnoootOes+GDRvyxRdfkJiYyI8//sjVV19d9A/6woULWbZsGU2bNqV///58++23DBgwoMRsV111FZMmTaJRo0b4fD6aNm3Kli3+Z6Dceeed/OpXv2LAgAFs3LiRCy64gBUrVpx0fKh9li9fzqxZs06abbW8rBAEyhwNS9+j45bFvMoEXn32P8yucwn1asXz5yu6kxhnRcGYSBMXF8eZZ57Jq6++yrPPPlu0PtQU1CNGjCApKYmkpCQGDx7MvHnz6NGjR9H2vLw8xo4dy6JFi/D5fKxatapoW58+fWjevDkAPXr0YP369UELwdChQ3nwwQdp1KjRSfMYAXz55ZcsX768aHnfvn0cOHDA9T7Dhw+vlCIAVghO1XUkPpTumxaTemwmnx8dxuw1u7igc2Mu6W4PsDEmqFL+cg+XmJgY3n33XYYMGcKjjz7K/fffD4Segrr4tPPFl5955hkaNWrE4sWLKSwsJDExsWhbQkJC0Wufz0d+fj5z587llltuAeDhhx+mW7duAMTHx9O7d2+eeuopli9fztSpU4uOLSwsZM6cOSe9d3Gh9qlVq1Zpp8Y16yMozpmDKKlFDzrHbOTDWo/yftIjrJr2PJt2H/I6nTGmBMnJyXzyySdFM4ZC6CmoP/zwQ44cOcKuXbuYMWMGWVlZJ23fu3cvTZo0ISYmhrfeeouCgtCDR/r27cuiRYtYtGgRw4cPP2nbuHHjeOKJJ6hXr95J688//3yef/75ouWSnn3gZp/KYIUgGGcOIkHoHLOB/oe/4bK/ziavwKapNiYS1atXj+nTp/PHP/6RqVOnhpyCulu3bgwePJh+/frx4IMP0rTpyVf7t956K2+++Sbdu3dn5cqVFfrru3Pnztxwww2nrH/uuefIzs6mW7dudOrUiRdffLFc+1SG6J2GugwKX7uQzXsOM3D7Pbw+OovBpzes0s83JlJVx2moJ0yYQEpKCvfcc4/XUcLGpqEOgxgRmh9dzZTEP7Ljm/BUZGOM8Yp1FrvRdSQCdN60kNjt04HfeZ3IGFNOJT23ONrZFYEbmaNh9Cf8lHo6R/MLOXC04nfyGVNTVLfm5ZquPP89rBCUQVK8/wLqh9z9HicxJjIkJiaya9cuKwYRQlXZtWtXyCGpJbGmoTJIdu4y/iF3P71b1fU4jTHea968OTk5OezYscPrKMaRmJhYdMObW64KgYi0Atqp6pcikgTEqmrU/VmcEBtDv5gV5C79G/R9oPQDjKnh4uLiTrpj11RPpTYNicjNwHvA8ae6Nwf+Fc5QkUq6XgFAy83TKCi0S2FjTM3gpo/gNqA/sA9AVX8EonMgfeZodqZn0UuXsf6zv3idxhhjKoWbQnBUVY8dXxCRWCBq/xxOy7ra/2Lpe94GMcaYSuKmEMwUkfuBJBE5D5gCfBTeWJErvu9NrE7uwe5DxziSZw+vMcZUf24KwXhgB7AUuAWYpqq/DWuqCJeeEk8Wy/njQ/eydseB0g8wxpgI5qYQ3K6qL6vqFao6UlVfFpE7w54sgqVlXQPAxTHfMnvNLo/TGGNMxbgpBKdOmwc3VnKOaiUmazTaqj8+EVbm7vM6jjHGVEjQ+whE5GrgGiBDRKYGbEoFdrt5cxEZCjwL+IBXVPXxYttbAm8CdZx9xqvqtDL9BB4RhOR4n91lbIyp9kLdUDYb2AqkA08FrN8PLCntjUXEB7wAnAfkAPNFZKqqLg/Y7QHgXVX9q4h0AqYBrcv0E3goOd7Hytz9qOopTzgyxpjqImghUNUNwAbgjHK+dx9gtaquBRCRScAIILAQKJDmvK4NbCnnZ3ki4+AiLs77jPYPFCD4C8Ho/q2578LqNT+7MSa6ubmzuJ+IzBeRAyJyTEQKRMRNw3gzYFPAco6zLtAE4H9EJAf/1cDtQTKMEZFsEcmOmDlNuo4E4Nb6C/nFwDbcNDCDzs3S+PvcjTas1BhTrbiZa+gvwFX47x/IBK4H2lfS518NvKGqT4nIGcBbItJFVU96HqSqTgQmgv8JZZX02RWTORqWvkeL3KX8Zus4APaSxwrdx6ppK+g24i6PAxpjjDuupqFW1dWAT1ULVPV1YKiLwzYDLQKWmzvrAt0EvOt8xn+BRPx9EtWD81zj49KSYukcs4FaP/zTw1DGGFM2bq4IDolIPLBIRJ7E34HspoDMB9qJSAb+AnAV/lFIgTYCQ4A3RKQj/kIQIW0/LmSO9n85BNj06EAKjlnTkDGm+nDzD/p1+Id2jgUO4v8r//LSDlLVfOeYz4AV+EcHLRORh0VkuLPbOOBmEVkMvAPcqNX8CRfJ8T4O5xXY7KTGmGqj1CsCZ/QQwGHgobK8uXNPwLRi634X8Ho5/plNa4zkeB+FqmzcfYiM9FpexzHGmFK5GTV0sYgsFJHdIrJPRPa7HDUUlZLjY+kXs4KDs1/2Oooxxrjipmnof/FPM1FfVdNUNVVV00o7KFol9LwSgGTrMDbGVBNuCsEm4Pvq3nZfVeL63sSPSd3ZffAYx/ILSz/AGGM85mbU0L3ANBGZCRw9vlJVnw5bqmouPSWBdocXs+LT5+l4SVRP1GqMqQbcXBE8AhzCP7QzNeDLBJHax/8Us9jl73ucxBhjSufmiqCpqnYJe5IaJDbr56yb8Td+OnCEA0fzSUlwc5qNMcYbbq4IponI+WFPUsOkp8ZTqPDZ97leRzHGmJDcFIJfAtNF5LANH3UvJSGWhNgY/rWo+KwaxhgTWUotBM5w0RhVTbLho+4JQs/CZbRYO5nt+494HccYY4IKWghEpIPzvVdJX1UXsZpypqke7pvNR4u3ehzGGGOCC9WLeTcwhpOfTnacAueEJVFN4UxT3W/DLBbMeQUGPOJ1ImOMKVGoJ5SNcV4OU9WT2jZEJDGsqWqKriNhwyx67/uKtTvuo02DFK8TGWPMKdx0Fs92uc4UlzmaY83PBOCN2etZkrOH7zfvJb/A7jg2xkSOoFcEItIY/6Mlk0SkJ3D86expQHIVZKsR4n0x9ItZwYfzXmP4f4cAMO689tw+pJ3HyYwxxi/UFcEFwJ/xP1nsqYCvu4H7wx+thnA6je9tupRXb8ikR4s6vP9dDjZ1kzEmUkhp/yCJyOWqGjFzJWRmZmp2drbXMcrm9Ysgdyk07sr2/UdYu/MgzQZeT4vzbvU6mTEmSojIAlXNLGmbmz6C5iKSJn6viMh3dqdxGQU827herXg6yQYKFr/rcShjjPFzUwh+rqr7gPOB+vgfXfl4WFPVNJmjYfQnMPoTYm/6lNykdrQ+sJCC+a95ncwYY1wVguOdxBcCf1PVZQHrTDnkd/Y/8nnjjL95nMQYY9wVggUi8jn+QvCZiKQCNv6xAjpefAfrUnqybf8R1uw44HUcY0yUc1MIbgLGA1mqegiIB0aHNVUNJyI0TE0AYOXW/R6nMcZEOzeFQIFOwB3Oci38D6kxFZAU5wPgh1ybyNUY4y03heD/gDOAq53l/cALYUsUJWJESIzzsTLXrgiMMd5yUwj6quptwBEAVf0Jf/OQqaD2uo5frr8Dsl/3OooxJoq5KQR5IuLD30SEiDTAOosrrutI9qR1oG3BOg4tmOR1GmNMFHNTCJ4D/gk0FJFHgFnAo2FNFQ0yR5M8ZjoraU3B1iXs+su5dmVgjPGEmyeU/R24F3gM2ApcqqpTwh0sGtRJjudox8tYVtiK+J3L0KV2Wo0xVS/Ug2mKqOpKYGWYs0SlgVf9mk+X/g/LpvyMjofzqe11IGNM1HHTNGTCbHCHhsSIUHv7XGseMsZUOSsEESAxzse6JsMAKFxizUPGmKpVaiEQkdtFpG553lxEhorIDyKyWkTGB9nnShFZLiLLROQf5fmcmqD+Wbcwp7AjR3MWk/P0YHKfHULhfLs6MMaEn5srgkbAfBF51/mH3dWEc86Q0xeAYfjvTL5aRDoV26cdcB/QX1U7A3eVKX0NMrB9Ov9JHMTi/Bbk7DlM8u7l7JsftXXRGFOFSu0sVtUHRORB/NNQjwb+IiLvAq+q6poQh/YBVqvqWgARmQSMAJYH7HMz8IJzkxqqur18P0b1lxDr4577HkMVjuYX8v0jA6h34Bh1vA5mjKnx3I4aUhHJBXKBfKAu8J6IfKGq9wY5rBmwKWA5B+hbbJ/2ACLyLeADJqjq9OJvJCJjgDEALVu2dBO5WhIRRCAp3kfdWvE0OLiK7AlnnNgOtGmYQt2kYjd2dx3pf+aBMcaUg5s+gjtFZAHwJPAt0FVVfwn0Bi6v4OfHAu2AQfjnMnpZRE75I1hVJ6pqpqpmNmjQoIIfWT3U6XMNu1NPp35KQtEXwPZ9R07eMXcpLH3Pg4TGmJrCzRVBPeAyVd0QuFJVC0Xk4hDHbQZaBCw3d9YFygHmqmoesE5EVuEvDPNd5KrR0gfdQvqgW05a949PlvPG7PU8dl43YmP8XTXnz/s5yV4ENMbUGG46i9sULwIi8haAqq4Icdx8oJ2IZIhIPHAVMLXYPv/CfzWAiKTjbypa6y569Lm0ZzPyC5V7pizmrsmLuGvyIlZs3Y/mLoHXL7J7EIwx5eLmiqBz4IIzGqh3aQepar6IjAU+w9/+/5qqLhORh4FsVZ3qbDtfRJYDBcCvVXVXWX+IaNG5aW3m3jeEg8cKAFiw4Sfe+6AfjWsn0jR3if/5odZXYIwpI1HVkjeI3AfcDyQBh46vBo4BE1X1vipJWExmZqZmZ2d78dERp6BQOevJb9i85zAz0v9E67y10LjriR2sE9kY4xCRBaqaWdK2oE1DqvqYqqYCf1LVNOcrVVXre1UEzMl8McJbN/WhVf1k3j7YB23c5cRG60Q2xrgUtGlIRDo4k81NEZFexber6ndhTWZcadMghfuGdeD/vX2IfanXUq+Wf3TRz/ffRkOPsxljqodQfQTj8N/w9VQJ2xQ4JyyJTJkNOr0hGem1+NeiLQAUFiqDYg+SfmAVMa9f5P6NrCnJmKgUtBCo6s3O98FVF8eUR2Kcj2/uGVS0vGHXQV58+kxaJi+iqds3yV3q/26FwJioE6pp6LJQB6rqB5Ufx1SGVvVrsbTxzxi87TzqHIpzdcyL+b/j9Lx8uyfBmCgUqmnokhDbFLBCEMF+e2EnPlxU/P694I4tLiRm2/f++xHAmomMiSKhmobsX4Fq7Iy29TmjbX3X+7+55QLidn9B28N5pO5Z4R9OZoXAmKgQqmnof1T1bRG5u6Ttqvp0+GKZqlZn4Bgum9QHNsKUxD/SW9WeWmRMlAjVNFTL+Z5aFUGMty7p1pSW9ZJZunkvBZ8qBVuXnBhxZM1ExtRooZqGXnK+P1R1cYxXYmKEni3r0qVZbR7/fCBJhbNJ2XGAVvlr/L8kVgiMqbHcTEPdRkQ+EpEdIrJdRD4UkTZVEc5UvThfDPXOuoVb4/7AhfvGs4oM2DDLJrQzpgZz0wz8D+BdoAnQFJgCvBPOUMZbtw0+jW/Hn8MtZ7Xh7YNZ/pU2XYUxNZabQpCsqm+par7z9TaQGO5gxnsXdmvCPwqGsK1uifNUGWNqiKCFQETqiUg94FMRGS8irUWklYjcC0yruojGK6c3SqVNg1rsOnjM6yjGmDAKNWpoAf4bx8RZDnxclgI2A2kNJyJc3LUJ+2blkVdQiLt7lI0x1U2oUUMZVRnERKYLuzVhzyzYffAYjbwOY4wJCzdPKENEugCdCOgbUNW/hSuUiRynN0plcZyPlD0rbPoJY2qoUguBiPwe/3OFO+HvGxgGzAKsEEQBEWFHq0s4svpf1M7dR0b+Wv9fA1YIjKkx3IwaGgkMAXKd+Ye6A7XDmspElA4X38EzzZ7htrg/sCivBYfy8r2OZIypRG4KwWFVLQTyRSQN2A60CG8sE0la1Etm8i1nMOmWfgAkb5ljN5gZU4O4KQTZIlIHeBn/SKLvgP+GNZWJSA1TE1la9zwAdOkUj9MYYypLqX0Eqnqr8/JFEZkOpKnqkvDGMpEqsd9NzPn0C7K2LsVX2mMwrVPZmGrB7aihy4AB+O8fmAVYIYhSF3RpzP9+fCZ1Y7JJ23s46H4ND/2ID6wQGFMNuBk19H/AaZyYX+gWETlXVW8LazITkRqmJrK9/dVcsGJIyP0+SXuMzlWUyRhTMW6uCM4BOqqqAojIm8CysKYyEW3idZkcPBZ85NBLM9ey/9t8CnMDnmlQUdbMZEzYuCkEq4GWwAZnuYWzzkSpmBghNTH4hBMjejTltX+fSULeHHyb9xStrxUfS9sGKWX/wNyl/u9WCIwJi1CPqvwIf59AKrBCROY5m/oA84IdZ0y7RqnUHjCGF7ZfU7Rux/4jLM7Zy5c3ns1pDctYDCrrqsIYU6JQVwR/rrIUpsYZP6zDScu5e4/Q77GveGfeRq7IbF6m92pbWEjc9mXBC4I1GxlTIaEmnZt5/LWINAKcJ5QwT1W3hzuYqVka106kT+t6vDprHa/OWlemY0cndmF8MyWhpI3WbGRMhbkZNXQl8CdgBv4pqZ8XkV+rqj2yypTJs1f3YNHGPaXvGOBofiG/eT+GuYeG077RqU1Kd8ffTQu0aK50Y0zZueks/i2QdfwqQEQaAF8CVghMmTSpnUSTrkllPu5wXgEvzVzDwk0nF5EjeQVsOXKYZkdXl35zW3lZs5OJAm4KQUyxpqBduJuaAhEZCjwL+IBXVPXxIPtdjr+wZKlqtpv3NtHj6j4tubpPy1PW7zl0jD892p8mCQtoFY4PtmYnEyXcFILpIvIZJ24oG4WLR1WKiA94ATgPyAHmi8hUVV1ebL9U4E5gblmCG1MnOZ6cNqMY9OM5xO2OYXT/1tx3YcfK+wAbrWSihJu5hn4dMMUEwERV/aeL9+4DrFbVtQAiMgkYASwvtt8fgCeAX7tObYzj/gs70nlRGrPX7OKdeRsZd/7pxMe6umB1J3dp6IJgTUemBghZCJy/6r9U1cHAB2V872bApoDlHKBvsffvBbRQ1U9EJGghEJExwBiAli1PbSIw0ev0xqncO7QDX63Yxk1vZvPqrHW0K+t9CkH0ajuceqF2sKYjU0OELASqWiAihSJSW1X3VuYHi0gM8DRwY2n7qupEYCJAZmamVmYOUzMMaJdO3eQ4npi+stLeMyWhNTN+/QHpKSUOXLWmI1NjuOkjOAAsFZEvgIPHV6rqHaUct5mTH2DT3Fl3XCrQBZghIgCNgakiMtw6jE1ZJcT6+Oyus9i272ilvN+ug0e56c1sfv7GfFrVr1XiPg/uP0LDgz/6C4I1EZlqzE0h+ICyNwsBzAfaiUgG/gJwFVA054BzhZF+fFlEZgD3WBEw5dUwLZGGaYmV9n7jh3bgnXkbWbb51IvhXQeP8RK9eKBlAmJNRKaac9NZ/KaIxAMd8M899IOqHnNxXL6IjAU+wz989DVVXSYiDwPZqjq1gtmNCaubz2rDzWe1KXHbW/9dz4Mf5vHznz1Ms3+OrNpgxlQyN3cWXwi8BKzBf2dxhojcoqqflnasqk6j2FBTVf1dkH0HuQlsTCTo0qw2AEtz9tIMYMMs/3Oc7arAVENuxtk9DQxW1UGqejYwGHgmvLGMiWwdm6ThixEmz99IXqfL/CuX2s32pnpyUwj2q2rg8wfWAvvDlMeYaiExzkfHJql888MOfr+lD7QaUPpBxkQoN4UgW0SmiciNInID8BH+u4Qvc240MyYqTbwuk/aNUpj+fS6Knrj5LPt1r6MZUyZuCkEisA04GxgE7ACSgEuAi8OWzJgI17ROEr86tz27Dx5j8pF+FDTq4i8G1kRkqhk3o4as98uYIM4+vQH1asUzfkNv9l14LWOktNtrjIk8bu4jMMYEkRwfy/zfnstFz/2Hz5ZtY0wiNj+RqXYqcXYuY6KTL0a4uFsTFmz4ibcPZVHYqEvwna3pyEQguyIwphJc1qs578zbxAObsmh64y85p0Ojkne0+YlMBApaCETk7lAHqurTlR/HmOqpaZ0kvrlnEL3/8AXTv88NXgig9KajUKxZyYRBqCuC1CpLYUwNEB8bw5CODXk3O4cPvttc4j6/TO3Grxpp+dpkbU4jEyZBC4GqPlSVQYypCX51Xnta1EumUE+dLX37vqM8v2AAvX72KwZ3aFj2N7dmJRMmbuYaSgRuAjrjv6cAAFX9eRhzGVMttapfi3Hnn17itqP5BUz/PpfJ8zeRklj27rmu+QUk7lxeekGw5iNTRm5+G98CVgIXAA8D1wIrwhnKmJooIdbHeZ0a8cHCzUxfllvm42+M78KDrWPwhdrJmo9MObgpBKep6hUiMsKZkvofwH/CHcyYmuihEZ25rFfzMh+3ducBfvch9Mkax4VdmwTf0ZqPTDm4KQR5zvc9ItIFyAXK0cBpjElNjGNAu/TSdyzmjLb1+d8vf2Tiv9eyevuBovVJcT6uP7MVCbEB1wkVGZVUGaxpqtpxUwgmikhd4EFgKpDivDbGVBFfjHBFZnNemrmWRZv2nLStbq14RvZ2rjK6evyQHGuaqpZESxjdcNIOIj5VLaiiPKXKzMzU7Gx7mqWJTgWFJ/5/VVUGPvkNXZrV5uXrMz1MFazvdyAAABBqSURBVOD4lcjoT7zNYU4hIgtUtcRfFDdXBOtEZDowGfhaS6scxpiw8cVIwJJwQefG/GPuRq6eOAeAIR0b8ouBJT9e05hg3BSCDvinm74NeE1EPgImqeqssCYzxpTq2r4t+XH7fvLylW37j/Dk9B+4tGcz0lMSvAvldR9FTda4Kwx7vNLf1s001IeAd4F3nb6CZ4GZEHoUmzEm/No1SuXvv+gHwJodBxjy1EyufXkuDVKDF4KYGOGe89vTrXmdyg/kdR+FKRdXd7WIyNnAKGAokA1cGc5Qxpiya9sghVsHtWXuut0czgverbdsy15e/3Y9z4zqUfkhMkdbR3E15ObO4vXAQvxXBb9W1YPhDmWMKZ97h3YodZ9x7y7mi+W55BUUEuezmeiNuyuCbqq6L+xJjDFV4oLOjXj/uxy6Tfic433PZ7RN55UbImTkkalyoaahvldVnwQeEZFTRgqpqj2Tz5hqaHCHhtx9Xnv2HfbfK7p+1yG+XLGNJTl7wtNvYCJeqCuC4/MJ2aB9Y2qQOF8MdwxpV7S870ge/R79itGvz6dOctwp+/dqWZc/XdG9KiOaKhZqGuqPnJdLVfW7KspjjKliaYlxPDS8MzNW7Thl27a9R5iyIIfR/TPo1DTNg3SmKri5s/gboDHwHjBZVb+vimDB2J3FxlSdPYeO0ffRr8hIr0VGeq2Q+8bECHec047TG9szrSJRhe4sVtXBItIY/5DRl0QkDX9B+GMl5zTGRJg6yfGMHXwaHy3ZwpodB0Luu2HXIQBeuKZXVUQzlajUK4KTdhbpCtwLjFLV+LClCsGuCIyJTA9/tJy35qzn3gs6EBMwFcaZbevTsYk1K3mtQlcEItIR/81klwO78M85NK5SExpjqr1r+rbk7TkbeGTayc+t6tQkjWl3DvQolXHDzX0ErwGTgAtUdUtZ3lxEhuKfksIHvKKqjxfbfjfwCyAf2AH8XFU3lOUzjDGR4bSGKSyZcD5H8wuL1v1t9nqe+mIVm/ccplmdJA/TmVBC3lYoIj5gnao+W44i4ANeAIYBnYCrRaRTsd0WApmq2g1/Z/STZfkMY0xkSYzzUTspruhrmPM0tZvemM+m3Yc8TmeCCVkInOcQtBCR8vQH9AFWq+paVT2G/6piRLH3/8aZ1A5gDlD2Z/gZYyJW2wa16JtRj5W5+3nmi1VexzFBuHoeAfCtiEwFiuYZUtWnSzmuGbApYDkH6Bti/5uAT0vaICJjgDEALVu2dBHZGBMJRITJt5zB/f9cygff5XD9ma2J8/k7kpvUTqJeLU/GnJhi3BSCNc5XDBCWAcIi8j9AJnB2SdtVdSIwEfyjhsKRwRgTPlf0bs4/5m7k0he+LVrXrE4S/753cLGH7RgvuLmP4KFyvvdmoEXAcnNn3UlE5Fzgt8DZqnq0nJ9ljIlgPVvWZfKYfuxx5jdatmUfz331I/9ds4sB7dI9TmfcDB/9Bihp0rlzSjl0PtBORDLwF4CrgGuKvXdP4CVgqKpudxvaGFP99G1Tv+j12e0b8Pq36/j1e4tpUjsx6DGjslowKsuag8PNTdPQPQGvE/HfT5Bf2kGqmi8iY4HP8A8ffU1Vl4nIw0C2qk4F/gSkAFNEBGCjqg4v489gjKlmEuN8jB/Wgenf5wbdZ832Azz9xSpG9m5hzUdhVqY7i4sOEpmnqn3CkKdUdmexMdHho8VbuP2dhfzu4k5kNAg9z9FxTWsn2VxHQVT0zuJ6AYsxQG+gdiVlM8aYEp3bsRFpibE8/PFy18fEx8Yw//5zqV3CdNomODdNQwvw9xEI/iahdfiHehpjTNgkxfv47Fdnkbv3iKv9N/10mDveWchny3K5MqtF6QeYIuVqGvKSNQ0ZY0qiqgz68wyO5hXSPqB5KDE2hkd+1pUGqQkepvNeqKahoHcWi0iWM/308eXrReRDEXmuWHORMcZ4TkS469x2NK6dyL7Deew7nMfew3l8vnwbH3yX43W8iBb0ikBEvgPOVdXdInIW/ikibgd6AB1VdWTVxTzBrgiMMWVxyfOziIkRPrytv9dRPFWuKwLAp6q7ndejgImq+r6qPgicVtkhjTEmHIZ1bcziTXvI+ckmvQsmZCEQkeOdyUOArwO2uelkNsYYzw3r4p8BNdQ9C9Eu1D/o7wAzRWQncBj4D4CInAbsrYJsxhhTYRnptejYJI33v9tMo7TgdzEDdGySymkNo+8+hKCFQFUfEZGvgCbA53qiMyEGf1+BMcZUCyN6NOXxT1dy+zsLQ+7XrE4Ss34zGGemg6gRsolHVeeUsM4mFTfGVCtjBrbh/E6NKAwxXP6rFdt57NOVLNuyjy7NouueWWvrN8bUeDExQpsGKSH3qZsczxPTV3Lr37/jlRsyad8oepqIQj6hzBhjokX9lASGdWnCxt2HeOLTlV7HqVJWCIwxxvHCtb244YxWzFq9k8PHCryOU2WsacgYYwIM6diIN/+7gTsmLSQ95cSjNHu1rMsVmTVzDiMrBMYYE6Bvm3r0bFmHxZv2FK07dKyAfy7czCXdm5IY5/MwXXhYITDGmAAJsT7+eevJ01HM+GE7N74+n/+u3cXg0xt6lCx8rBAYY0wp+rWpT2JcDP/7xSq+WrGtaH3d5HjuHNKOWF/17m61QmCMMaVIjPNxVVZLPlq8hZyfDgOQX6jsPZxHvzb16X9auscJK8YKgTHGuDBheGcmDO9ctHzwaD49H/6Cmat2WCEwxphoVCshlqyMunyyZCt1k+NLPyCI7s1rc6bHhcQKgTHGlNMl3Zoy/oOlPDG9/Deg1UmOY8ED5+GL8W5+IysExhhTTlf1acmlPZuV+/hPlmxl3JTFfL95L91b1KnEZGVjhcAYYyqgIvcVDDq9AQDvLcjhaH5h0fo6yXFVOteRFQJjjPFI/ZQEureow1tzNvDWnA1F60Xgm3GDaJ1eq0pyWCEwxhgPvXxdb37cfqBoeffBY9z+zkL+8+MOKwTGGBMNGqYl0jDgyWmqyuOfruTb1bu47ozWVZLBCoExxkQQEeHMtvX5YOFm+j76JQA3D2zDLwa2CdtnVu/7oo0xpgYac1YbrsxsweDTG5IY5+Pd7E1h/Ty7IjDGmAjTrlEqj13WFYD/m7GaJ6f/wM4DR0lPSQjL59kVgTHGRLB+beoDMG/d7rB9hhUCY4yJYF2b1SY53sectbvC9hlhLQQiMlREfhCR1SIyvoTtCSIy2dk+V0RahzOPMcZUN3G+GDJb16uehUBEfMALwDCgE3C1iHQqtttNwE+qehrwDPBEuPIYY0x11TejHqu2HWDXgaNhef9wXhH0AVar6lpVPQZMAkYU22cE8Kbz+j1giIh4N/OSMcZEoHD3E4SzEDQDAsc85TjrStxHVfOBvUD94m8kImNEJFtEsnfs2BGmuMYYE5m6Na/NOR0akpwQnoGe1WL4qKpOBCYCZGZmqsdxjDGmSsX5YnjtxqywvX84rwg2Ay0Clps760rcR0RigdpA+HpEjDHGnCKchWA+0E5EMkQkHrgKmFpsn6nADc7rkcDXqmp/8RtjTBUKW9OQquaLyFjgM8AHvKaqy0TkYSBbVacCrwJvichqYDf+YmGMMaYKhbWPQFWnAdOKrftdwOsjwBXhzGCMMSY0u7PYGGOinBUCY4yJclYIjDEmylkhMMaYKCfVbbSmiOwANpS6Y8nSgZ2VGCfcLG94Wd7wsrzhVda8rVS1QUkbql0hqAgRyVbVTK9zuGV5w8vyhpflDa/KzGtNQ8YYE+WsEBhjTJSLtkIw0esAZWR5w8vyhpflDa9KyxtVfQTGGGNOFW1XBMYYY4qxQmCMMVEuagqBiAwVkR9EZLWIjPc6T0lEZL2ILBWRRSKS7ayrJyJfiMiPzve6HuZ7TUS2i8j3AetKzCd+zznne4mI9IqQvBNEZLNzjheJyIUB2+5z8v4gIhd4kLeFiHwjIstFZJmI3Omsj7hzHCJrJJ/fRBGZJyKLncwPOeszRGSuk22yM20+IpLgLK92treOgKxviMi6gPPbw1lfsd8FVa3xX/inwV4DtAHigcVAJ69zlZBzPZBebN2TwHjn9XjgCQ/znQX0Ar4vLR9wIfApIEA/YG6E5J0A3FPCvp2c34sEIMP5ffFVcd4mQC/ndSqwyskVcec4RNZIPr8CpDiv44C5znl7F7jKWf8i8Evn9a3Ai87rq4DJEZD1DWBkCftX6HchWq4I+gCrVXWtqh4DJgEjPM7k1gjgTef1m8ClXgVR1X/jf25EoGD5RgB/U785QB0RaVI1Sf2C5A1mBDBJVY+q6jpgNf7fmyqjqltV9Tvn9X5gBf7nekfcOQ6RNZhIOL+qqgecxTjnS4FzgPec9cXP7/Hz/h4wRETE46zBVOh3IVoKQTNgU8ByDqF/ab2iwOciskBExjjrGqnqVud1LtDIm2hBBcsXyed8rHP5/FpAU1tE5XWaIXri/0swos9xsawQwedXRHwisgjYDnyB/8pkj6rml5CrKLOzfS9Q36usqnr8/D7inN9nRCSheFZHmc5vtBSC6mKAqvYChgG3ichZgRvVfw0YseN9Iz2f469AW6AHsBV4yts4pxKRFOB94C5V3Re4LdLOcQlZI/r8qmqBqvbA/wz1PkAHjyMFVTyriHQB7sOfOQuoB/ymMj4rWgrBZqBFwHJzZ11EUdXNzvftwD/x/6JuO36J53zf7l3CEgXLF5HnXFW3Of+DFQIvc6J5IiLyikgc/n9Y/66qHzirI/Icl5Q10s/vcaq6B/gGOAN/M8rxpzUG5irK7GyvDeyq4qiBWYc6TXKqqkeB16mk8xsthWA+0M4ZHRCPv+NnqseZTiIitUQk9fhr4Hzge/w5b3B2uwH40JuEQQXLNxW43hnN0A/YG9C84Zli7aY/w3+OwZ/3KmekSAbQDphXxdkE/3O8V6jq0wGbIu4cB8sa4ee3gYjUcV4nAefh79v4Bhjp7Fb8/B4/7yOBr50rMq+yrgz4g0Dw92UEnt/y/y5UVS+411/4e9VX4W8T/K3XeUrI1wb/qIrFwLLjGfG3SX4F/Ah8CdTzMOM7+C/38/C3Qd4ULB/+0QsvOOd7KZAZIXnfcvIscf7naRKw/2+dvD8AwzzIOwB/s88SYJHzdWEknuMQWSP5/HYDFjrZvgd+56xvg78orQamAAnO+kRnebWzvU0EZP3aOb/fA29zYmRRhX4XbIoJY4yJctHSNGSMMSYIKwTGGBPlrBAYY0yUs0JgjDFRzgqBMcZEudjSdzEmOonI8WGbAI2BAmCHs3xIVc/0JJgxlcyGjxrjgohMAA6o6p+9zmJMZbOmIWPKQUQOON8HichMEflQRNaKyOMicq0zl/xSEWnr7NdARN4XkfnOV39vfwJjTrBCYEzFdQf+H9ARuA5or6p9gFeA2519ngWeUdUs4HJnmzERwfoIjKm4+erM6yIia4DPnfVLgcHO63OBTgHT2aeJSIqemHPeGM9YITCm4o4GvC4MWC7kxP9jMUA/VT1SlcGMccOahoypGp9zopmI48+aNSYSWCEwpmrcAWQ6T5Zajr9PwZiIYMNHjTEmytkVgTHGRDkrBMYYE+WsEBhjTJSzQmCMMVHOCoExxkQ5KwTGGBPlrBAYY0yU+/+hUNo3ws1h5QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "max_time = data.Time.max()\n",
    "x = range(0, max_time+1)\n",
    "y = np.zeros(len(x))\n",
    "for i, t in enumerate(x):\n",
    "    y[i] = naive_estimator(t, data)\n",
    "    \n",
    "plt.plot(x, y, label=\"Naive\")\n",
    "\n",
    "x, y = HomemadeKM(data)\n",
    "plt.step(x, y, label=\"Kaplan-Meier\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Survival probability estimate\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "iY__6ufG3sDk"
   },
   "source": [
    "### Question\n",
    "\n",
    "What differences do you observe between the naive estimator and Kaplan-Meier estimator? Do any of our earlier explorations of the dataset help to explain these differences?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "i7tElIKVoQ4R"
   },
   "source": [
    "<a name='5'></a>\n",
    "## 5. Subgroup Analysis\n",
    "\n",
    "We see that along with Time and Censor, we have a column called `Stage_group`. \n",
    "- A value of 1 in this column denotes a patient with stage III cancer\n",
    "- A value of 2 denotes stage IV. \n",
    "\n",
    "We want to compare the survival functions of these two groups.\n",
    "\n",
    "This time we'll use the `KaplanMeierFitter` class from `lifelines`. Run the next cell to fit and plot the Kaplan Meier curves for each group. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 279
    },
    "colab_type": "code",
    "id": "Ge6P3fgVrZLS",
    "outputId": "efbf8e54-7623-4d96-e24a-2fbf07ae2aac"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de7yVZZ338c9XPKCIpsCkchAsFEEcsCU2j74SDwlaSWOU4NOUZGElWjrZozVa2eNkJ6csnEIHtWYmRMpHUkyzKIfGDFACwQOEmhsxFUdQ0Tz9nj/ue+Fiuffa9957ne69vu/Xa7/2ug/rXr/FgvXjun73dV2KCMzMrHXt0OgAzMyssZwIzMxanBOBmVmLcyIwM2txTgRmZi1ux0YH0FUDBw6M4cOHNzoMM7NcWb58+dMRMai9Y7lLBMOHD2fZsmWNDsPMLFckPdrRMXcNmZm1OCcCM7MW50RgZtbiclcjMLPW8corr9DW1sZLL73U6FByo2/fvgwZMoSddtop83OcCMysabW1tdG/f3+GDx+OpEaH0/Qigk2bNtHW1saIESMyP69mXUOS5kp6UtJ9HRyXpCskrZO0UtJhtYrFzPLppZdeYsCAAU4CGUliwIABXW5B1bJGcC0wucLxE4GR6c9M4F9rGIuZ5ZSTQNd058+rZl1DEXGnpOEVTpkC/CiSebB/L+ktkvaNiI2VrvvixgdY/c9Hbdv+3a7H8KvdTmLKuMGcdsSwaoRuZtZSGnnX0GDgsZLttnTfm0iaKWmZpGWl6ycMf2U9R764mDUbt3DTig21jdbMWtKll17KmDFjOPTQQxk3bhx33303AN/5znfYunVrTV/7kUce4ZBDDgHgN7/5De9973sBuPbaa5k1a1bVXicXxeKImAPMASgUCjHmC0uSA9e8hzHA6Jf3aFxwZtZr3XXXXdx8883cc8897LLLLjz99NO8/PLLQJIIPvzhD7Pbbrs1OMqea2SLYAMwtGR7SLrPzKwpbNy4kYEDB7LLLrsAMHDgQPbbbz+uuOIKHn/8cY455hiOOeYYAD71qU9RKBQYM2YMX/rSl7ZdY9GiRYwaNYp3vOMdnHPOOdv+V//CCy/wsY99jAkTJjB+/Hhuuumm+r/BVCNbBAuBWZLmAUcAmzurD7TriVVcHOfzu12PAf6uyiGaWbP4ys9Xs+bxLVW95uj99uBL7xvT4fETTjiBSy65hAMPPJDjjz+eU089laOPPppzzjmHyy+/nMWLFzNw4EAg6ULae++9ee211zjuuONYuXIlBx54IGeeeSZ33nknI0aMYPr06duufemll3Lssccyd+5cnn32WSZMmMDxxx9Pv379qvoes6jl7aM/Ae4CDpLUJukMSZ+U9Mn0lEXAemAdcBXw6S6/yNipsM/YbbUCM7Nq2n333Vm+fDlz5sxh0KBBnHrqqVx77bXtnjt//nwOO+wwxo8fz+rVq1mzZg0PPPAABxxwwLZ7+ksTwe23385ll13GuHHjmDhxIi+99BJ//vOf6/G23qSWdw1N7+R4AGf16EUKM6Awg0dK7iIys96p0v/ca6lPnz5MnDiRiRMnMnbsWK677jpOP/307c55+OGH+da3vsXSpUvZa6+9OP300zu9lz8i+OlPf8pBBx1Uw+iz8VxDZmYdePDBB1m7du227RUrVrD//vsD0L9/f5577jkAtmzZQr9+/dhzzz35y1/+wq233grAQQcdxPr163nkkUcAuP7667dda9KkSXzve9+jeCfkvffeW4+31K5c3DWUxfBX1sM173ljx9ipSYvBzKybnn/+ec4++2yeffZZdtxxR97+9rczZ84cAGbOnMnkyZPZb7/9WLx4MePHj2fUqFEMHTqUI488EoBdd92VK6+8ksmTJ9OvXz8OP/zwbde+6KKL+OxnP8uhhx7K66+/zogRI7j55psb8j5Vel9+HhQKhShfmGbOv1zEkS8uZsy+eyY7nlgF+4yFGbc0IEIzq5b777+fgw8+uNFh9Mjzzz/P7rvvTkRw1llnMXLkSM4999yavmZ7f26SlkdEob3ze0XX0K92O4lLBnwz+eKfcUuSBMzMmsBVV13FuHHjGDNmDJs3b+bMM89sdEhv0mu6hszMmtG5555b8xZAT/WKFkG7Hl0Cy65pdBRmZk2v17QI1mzcwqk/vAuA47YexkyWsPr2f+OS5aPaPd+T1JmZJXpFIpgybvu56n6120kVB5it2ZiMTnQiMDPrJYngtCOGvflL/ZrkDqLrZ7x52oliy8HMzHpzjcDMrAqaYRrqrVu3MmDAALZs2X6upfe///3bDVLrLicCM7MOlE5DvXLlSu644w6GDk0mTa5HIijabbfdmDRpEjfeeOO2fZs3b2bJkiW8733v6/H1nQjMzDrQTNNQT58+nXnz5m3bvvHGG5k0aVJV1kPoFTUCM2sBt16QzBpQTfuMhRMv6/BwM01DPWnSJD7+8Y+zadMmBgwYwLx586q2SplbBGZmHWimaah33nlnTj75ZBYsWMDTTz/Nvffey6RJk6ryPt0iMLN8qPA/91pqpmmop0+fzle/+lUigilTprDTTjt15y29iVsEZmYdaLZpqCdOnMjatWuZPXv2dq2Lnqppi0DSZOC7QB/g6oi4rOz4/sBcYBDwDPDhiGirZUxmZlk12zTUO+ywA1OnTmX+/PkcffTRVXufNZuGWlIf4CHg3UAbsBSYHhFrSs65Abg5Iq6TdCwwIyL+odJ125uGul3XvOeN6ajL1iYoDii7/kyvcWzWzDwNdfc00zTUE4B1EbE+Il4G5gFTys4ZDfw6fby4nePdl65nzBOrYNWCql3WzKwrWn0a6sHAYyXbbcARZef8ETiFpPvo74H+kgZExKbSkyTNBGYCDBuWcX6gdD3j7VYtK1E6SV1WnqjOzLrK01B37nPA0ZLuBY4GNgCvlZ8UEXMiohARhUGDBvX4RaeMG8zofffo0nPWbNzCTSs29Pi1zaxr8raKYqN158+rli2CDcDQku0h6b5tIuJxkhYBknYHPhARz1Y9kidWbdcyOG3sVE47s2vrGXuiOrP669u377YBVJIaHU7Tiwg2bdpE3759u/S8WiaCpcBISSNIEsA04LTSEyQNBJ6JiNeBC0nuIKqusVO33y6OTPTC9mZNb8iQIbS1tfHUU081OpTc6Nu3L0OGDOnSczIlgvQ2z5ERcYekXYEdI+K5Ss+JiFclzQJuI7l9dG5ErJZ0CbAsIhYCE4GvSQrgTuCsLkWfRbFWUNRBzcDMms9OO+20bVSu1U6niUDSJ0gKtXsDbyPp4vkBcFxnz42IRcCisn0XlzxeAPiWHjOzBspSLD4LOBLYAhARa4G/qWVQNVesGXhNYzOzTF1Df42Il4uFGkk7Avkt4xdrBq4VmJkB2VoEv5X0BWBXSe8GbgB+XtuwaqgwA2bckgw2MzOzTIngAuApYBVwJrAoIr5Y06jMzKxusnQNnR0R3wWuKu6Q9Jl0n5mZ5VyWFsFH29l3epXjMDOzBumwRSBpOskAsBGSFpYc6k8yZbSZmfUClbqG/hvYCAwEvl2y/zlgZS2DalZZJ6rz5HRmlicdJoKIeBR4FPCk/SRf7lms2bgFwInAzHIjy8jidwLfAw4GdiaZLuKFiOja9J3N6NElyaCyDGMJTjtiWKYvd09OZ2Z5k6VY/H1gOrAW2BX4ODC7lkHVRXFgmRetMbMWl2k9gohYB/SJiNci4hpgcm3DqoPCDNj/qEZHYWbWcFnGEWyVtDOwQtI3SArIjV7QxszMqiTLF/o/kNQFZgEvkCw284FaBlVXxTqBmVmL6jQRRMSjEfFiRGyJiK9ExHlpV1H+uU5gZtZ5IpD0Xkn3SnpG0hZJz0nakuXikiZLelDSOkkXtHN8mKTF6fVXSjqpO2+i21wnMDPL1DX0HZJpJgZExB4R0T/LraOS+pDcXXQiMBqYLml02Wn/BMyPiPEkS1le2aXozcysx7IkgseA+yKiq2sQTADWRcT6iHgZmAdMKTsngGJS2RN4vIuvUR1eqMbMWliWu4Y+DyyS9Fvgr8WdEXF5J88bTJJEitqAI8rO+TJwu6SzgX7A8e1dSNJMkuUyGTasyiN2vVCNmbW4LC2CS4GtQF+SCeeKP9UwHbg2IoYAJwE/lvSmmCJiTkQUIqIwaNCgKr10ygvVmFmLy9Ii2C8iDunGtTeQ3GpaNCTdV+oM0sFpEXGXpL4kk9w92Y3XMzOzbsjSIlgk6YRuXHspMFLSiHRA2jRgYdk5fwaOA5B0MEmr46luvJaZmXVTlkTwKeAXkl7syu2jEfEqySC024D7Se4OWi3pEkknp6f9I/AJSX8EfgKc3o2itJmZ9UCnXUMR0e16QEQsAhaV7bu45PEa4MjuXt/MzHqu0gployLiAUmHtXc8Iu6pXVhmZlYvlVoE55Hcsvntdo4FcGxNImqk4niCUmOn+pZSM+vVKq1QNjN9eGJEvFR6LL27p3cpjico5bEFZtYCstw++t9AefdQe/vyrTDjzV/45a0DM7NeqFKNYB+S0cG7ShoPKD20B7BbHWIzM7M6qNQimAScTjIQ7Nu8kQieA75Q27CaSBfWNTYzy6NKNYLrgOskfSAiflrHmJrH2KlJIli1wInAzHqtLAPKhkjaQ4mrJd3TzZHG+eP1CsysBWRJBB+LiC3ACcAAkqUrL6tpVGZmVjdZ7hoq1gZOAn6UThOhSk9odXc//Ayn/vCuRofRY1PGDea0I6o87beZNZ0siWC5pNuBEcCFkvoDr9c2rPyaMm5wo0OoijUbk+mknAjMer8sieAMYBywPiK2ShoAuHLagdOOGNYrvjx7Q4vGzLLJUiMIkjWHz0m3+5FMF21mZr1AlkRwJfB3JKuJQTKOYHbNImpGxbEEZma9UJZEcEREnAW8BBAR/wPsXNOomklxDqJVCxobh5lZjWRJBK9I6kPSRYSkQbRSsdhjCcysl8uSCK4AbgT+RtKlwBLgn7NcXNJkSQ9KWifpgnaO/4ukFenPQ5Ke7VL0ZmbWY1lWKPsPSctJ1hYW8P6IuL+z56WtiNnAu4E2YKmkhemqZMVrn1ty/tnA+K6/hToprlXg9QnMrJfJcvsoEfEA8EAXrz0BWBcR6wEkzQOmAGs6OH868KUuvkZ9FOsEXp/AzHqhLF1D3TUYeKxkuy3d9yaS9icZsPbrDo7PlLRM0rKnnnqq6oF2qjADZtwC+4yt/2ubmdVYLRNBV0wDFkTEa+0djIg5EVGIiMKgQYPqHJqZWe/WaSKQdLakvbpx7Q3A0JLtIem+9kwDftKN16g/jykws14mS4vgrSSF3vnpXUBZJ5xbCoyUNELSziRf9gvLT5I0CtgLaP45DTymwMx6oSx3Df2TpItIpqGeAXxf0nzg3yLiTxWe96qkWcBtQB9gbjpz6SXAsogoJoVpwLyIiJ6+mZorzGipJLBm4xbPOZQznjHWuiPrXUMh6QngCeBVkv/BL5D0y4j4fIXnLQIWle27uGz7y10N2mqvt8yi2ko8Y6x1V6eJQNJngI8ATwNXA+dHxCuSdgDWAh0mgl6rBdYx7i2zqLYSt96su7LUCPYGTomISRFxQ0S8AhARrwPvrWl0zch1AjPrZbIkggMi4tHSHZJ+DJBlhHGv47mHzKyXyZIIxpRupFNHvKM24ZiZWb11WCOQdCHwBWBXSVuKu4GXgTl1iK25FeceKvIcRGaWUx0mgoj4GvA1SV+LiAvrGFPzK9YJijwHkZnlWKUWwah0srkbJB1Wfjwi7qlpZM2sMGP7L/3SloGZWc5Uun30H4FPAN9u51gAx9YkIjMzq6tKXUOfSH8fU79wcqwFxhaYWe9UqWvolEpPjIifVT+cnBo7NUkEqxY4EZhZ7lTqGnpfhWMBOBEUtdgcRGbWu1TqGvJ/bc3MWkClrqEPR8S/SzqvveMRcXntwsqpR5e8+Q4ijy8wsyZXqWuoX/q7fz0Cyb3ysQXg8QVmlguVuoZ+mP7+Sv3CybHysQXg8QVmlgtZpqE+APgu8E6SIvFdwLkRsT7Dcyenz+0DXB0Rl7VzzoeAL6fX/mNEnNaVN2Bmb8jjYkJeTKfxskw695/AfGBfYD/gBjKsL5xOTjcbOBEYDUyXNLrsnJHAhcCRETEG+GyXos8Dr3FsdTJl3GBG77tHo8PokjUbt3DTio6WMrd6ybJC2W4R8eOS7X+XdH6G500A1hVbDpLmAVOANSXnfAKYHRH/AxART2YLOyc8vsDqKI+LCeWt9dJbddgikLS3pL2BWyVdIGm4pP0lfZ6y5Sc7MBh4rGS7Ld1X6kDgQEm/k/T7tCup9/DaBWaWA5VaBMtJ+u2Vbp9ZcixIunSq8fojgYnAEOBOSWMj4tnSkyTNBGYCDBuWr//xmJk1u0p3DY3o4bU3AENLtoek+0q1AXeny18+LOkhksSwtCyWOaRrIBQKhehhXPVXunaBxxWYWZPJUiNA0iEkBd++xX0R8aNOnrYUGClpBEkCmAaU3xH0/4DpwDWSBpJ0FXV6N1KulI4v8LgCM2tCWW4f/RJJ181oktrAicASoGIiiIhXJc0CbiO5fXRuRKyWdAmwLCIWpsdOkLQGeA04PyI29eD9NJ/S8QUeV2BmTShLi2Aq8LfAvRExQ9JbgX/PcvGIWERZYTkiLi55HMB56Y+ZmTVAlnEEL0bE68CrkvYAnmT7vn8zM8uxLC2CZZLeAlxFcifR8ySji607yhe9L3IR2cwapNNEEBGfTh/+QNIvgD0iYmVtw+ql2puYDlxENrOGynrX0CnAUSTjB5YATgTd0d7EdOAispk1VKc1AklXAp8EVgH3AWdKml3rwMzMrD6ytAiOBQ5O7/BB0nXA6ppG1Yo6qh2A6wdmVlNZEsE6YBjwaLo9NN1n1dJR7QBcPzCzmqu0VOXPSWoC/YH7Jf0hPTQB+ENHz7Nu6Kh2AK4fmFnNVWoRfKtuUZiZWcNUmnTut8XH6Wjiw9PNP/S6dQOaXXv1A9cNzKxKstw19CGSrqAPAh8C7pZUoVPbqmrsVNhn7Pb7nliVLHZjZlYFWYrFXwQOL7YCJA0C7gD8TVQP7dUPXDcwsyrKkgh2KOsK2kS2OYrMzDq1ZuMWL1nZYFkSwS8k3cYbC9afSralKs3MKpoyrnz1WmsEpePEKp/0xhQTAP8VETfWNKoKCoVCLFu2rFEv3xyueU9SJyivHXTGBWazliVpeUQU2jtWsUUgqQ9wR0QcA/ysFsFZN1QagNYRD0wzsw5UTAQR8Zqk1yXtGRGbu3pxSZOB75KsUHZ1RFxWdvx04Ju8sZbx9yPi6q6+TsupNACtIy4wm1kHstQIngdWSfol8EJxZ0ScU+lJaWtiNvBukkXql0paGBFryk69PiJmdS1sMzOrliyJ4Gd0r1toArAuItYDSJoHTAHKE4HVS+nANNcLzCyVZWGa6yTtDIwimXvowYh4OcO1BwOPlWy3AUe0c94HJL0LeAg4NyIeKz9B0kxgJsCwYcMyvLS9SWldwfUCMyuRZWTxScCfgCuA7wPrJJ1Ypdf/OTA8Ig4Ffglc195JETEnIgoRURg0aFCVXrrFFGbAjFuSn67ebWRmvVqWrqHLgWMiYh2ApLcBtwC3dvK8DWy/yP0Q3igKAxARm0o2rwa+kSEeMzOroiwjhJ8rJoHUeuC5DM9bCoyUNCLtWpoGLCw9QdK+JZsnA/dnuK6ZmVVRlhbBMkmLgPkkNYIPktwBdApARLRbSI6IVyXNAm4juX10bkSslnQJsCwiFgLnSDoZeBV4Bji9p2/IMnLh2MxSWRJBX+AvwNHp9lPArsD7SBJDh3cURcQiyqajiIiLSx5fCFzYtZCtx1w4NrMSWe4a8jdEb1M6IM0DzcxanmcRNTNrcU4EZmYtzonA4NElsOyaRkdhZg3SYY1A0nmVnhgRl1c/HKu7sVOTRLBqgQvGZi2qUrG4f92isMYpzPD6x2YtrsNEEBFfqWcgZmbWGJ3ePiqpL3AGMIZkTAEAEfGxGsZl9VasE7h7yKzlZCkW/xjYB5gE/JZkzqAsU0xYXhQHmLmLyKwlZUkEb4+Ii4AXIuI64D20P5205VVhBux/VOfnmVmvlCURvJL+flbSIcCewN/ULiQzM6unLHMNzZG0F3ARyeyhu6ePzcysF8iSCK6JiNdI6gMH1DgeayTPSGrWkrJ0DT0saY6k4ySp5hFZY4yd+sbKZU+scuHYrIVkSQSjgDuAs4BHJH1fkiuLvY2XsjRrWZ0mgojYGhHzI+IUYBywB0k3UackTZb0oKR1ki6ocN4HJIWkQubIzcysKrLUCJB0NHAqMBlYBnwow3P6ALOBdwNtJKuaLYyINWXn9Qc+A9zdtdCtpkrrBZ1xPcEs17KMLH4EuJdkqcrzI+KFjNeeAKyLiPXpdeYBU4A1Zed9Ffg6cH7G61qtla5g1hmvcGaWe1laBIdGxJZuXHsw8FjJdhtlA9EkHQYMjYhbJDkRNIvSFcw64xXOzHKv0jTUn4+IbwCXSory4xFxTk9eWNIOwOVkWLBe0kxgJsCwYcN68rJmZlamUovg/vT3sm5eewMwtGR7SLqvqD9wCPCb9K7UfYCFkk6OiO1eMyLmAHMACoXCm5KSmZl1X6VpqH+ePlwVEfd049pLgZGSRpAkgGnAaSXX3wwMLG5L+g3wufIkYDnQUWHZRWSzXMgyjuDbku6X9NV0rqFMIuJVYBZwG0nrYn5ErJZ0iaSTuxmvNZvSgWilPCjNLDcU0XlPi6R9SG4ZPZVkHMH1EfF/axxbuwqFQixb5kZD0yu2EGbc0tg4zAwAScsjot2xWpkWr4+IJyLiCuCTwArg4irGZ2ZmDZRlHMHBJC2BDwCbgOuBf6xxXNYbdGVQWle5/mBWNVnGEcwF5gGTIuLxGsdjvUVXBqV1lQexmVVVxUSQThPxcER8t07xWG/RlUFpXeVBbGZVVbFGkK5DMFTSznWKx8zM6ixL19DDwO8kLQS2zTMUEZfXLCozM6ubLIngT+nPDiSjgc0az6upmVVNp4kgIr5Sj0DMMistRLtwbNZjWW4fXQy0N+ncsTWJyKwzpYVoF47NeixL19DnSh73JRlP8GptwjEzs3rL0jW0vGzX7yT9oUbxmHVdLQeudcb1CesFsnQN7V2yuQPwDmDPmkVk1hW1HLjWGdcnrJfI0jW0nKRGIJIuoYeBM2oZlFlmtRy41hnXJ6yXyNI1NKIegZiZWWN0OLJY0uHp9NPF7Y9IuknSFWXdRWZmlmOVppj4IfAygKR3AZcBPwI2ky4badbyioXqZdc0OhKzbquUCPpExDPp41OBORHx04i4CHh7lotLmizpQUnrJF3QzvFPSlolaYWkJZJGd/0tmDVIcXU2r8ZmOVcxEUgq1hCOA35dcizL3UZ9gNnAicBoYHo7X/T/GRFjI2Ic8A3A8xdZfhRmJCuwtbdUp1mOVPpC/wnwW0lPAy8C/wUg6e0k3UOdmQCsi4j16fPmAVOANcUTImJLyfn9aGcEs5mZ1VaHiSAiLpX0K2Bf4PZ4Y3HjHYCzM1x7MPBYyXYbcET5SZLOAs4DdgbanbZC0kxgJsCwYcMyvLSZmWXV2XoEv4+IGyOidPrphyLinmoFEBGzI+JtwP8B/qmDc+ZERCEiCoMGDarWS5tVz6NLXDC23Mq0eH03bQCGlmwPSfd1ZB7w/hrGY1YbxdHNLhhbTtUyESwFRkoaka5wNg1YWHqCpJElm+8B1tYwHrPaKMyA/Y9qdBRm3ZZlioluiYhXJc0CbgP6AHMjYrWkS4BlEbEQmCXpeOAV4H+Aj9YqHjMza1/NEgFARCwCFpXtu7jk8Wdq+fpmZta5WnYNmbUWjzK2nKppi8CsZRQLxp6a2nLILQKzavAoY8sxJwIzsxbnRGBm1uJcIzCrtixrKHutY2siTgRm1ZRlDWUXlK3JOBGYVVOWNZS91rE1GdcIzMxanFsEZo2QpY7QGdcZrEqcCMzqLUsdoTOuM1gVORGY1VuWOkJnXGewKnKNwMysxTkRmJm1OHcNmeVVNQrOpVx8bllOBGZ5VI2CcykXn1taTROBpMnAd0lWKLs6Ii4rO34e8HHgVeAp4GMR8WgtYzLrFapRcC7l4nNLq1mNQFIfYDZwIjAamC5pdNlp9wKFiDgUWAB8o1bxmJlZ+2pZLJ4ArIuI9RHxMjAPmFJ6QkQsjoit6ebvgSE1jMfMzNpRy66hwcBjJdttwBEVzj8DuLW9A5JmAjMBhg0bVq34zKxUtYvP1jw6WTCpKYrFkj4MFICj2zseEXOAOQCFQiHqGJpZa6h28dlypZaJYAMwtGR7SLpvO5KOB74IHB0Rf61hPGbWkWoXn60Jfb3DI7WsESwFRkoaIWlnYBqwsPQESeOBHwInR8STNYzFzMw6ULNEEBGvArOA24D7gfkRsVrSJZJOTk/7JrA7cIOkFZIWdnA5MzOrkZrWCCJiEbCobN/FJY+Pr+Xrm5lZ5zzXkJlZi3MiMDNrcU4EZmYtzonAzKzFORGYmbU4ReRroK6k54AHGx1HDwwEnm50ED2U9/eQ9/gh/+/B8dff/hExqL0DTTHFRBc9GBGFRgfRXZKW5Tl+yP97yHv8kP/34Pibi7uGzMxanBOBmVmLy2MimNPoAHoo7/FD/t9D3uOH/L8Hx99EclcsNjOz6spji8DMzKrIicDMrMXlKhFImizpQUnrJF3Q6HiykPSIpFXpNNvL0n17S/qlpLXp770aHWcpSXMlPSnpvpJ97casxBXpZ7JS0mGNi3xbrO3F/2VJG9LPYYWkk0qOXZjG/6CkSY2J+g2ShkpaLGmNpNWSPpPuz8VnUCH+PH0GfSX9QdIf0/fwlXT/CEl3p7Fen661gqRd0u116fHhjYy/yyIiFz9AH+BPwAHAzsAfgdGNjitD3I8AA8v2fQO4IH18AfD1RsdZFt+7gMOA+zqLGTiJZK1pAe8E7m7S+L8MfK6dc0enf5d2AUakf8f6NDj+fYHD0sf9gYfSOHPxGVSIP0+fgYDd08c7AXenf7bzgWnp/h8An0offxr4Qfp4GnB9I+Pv6k+eWgQTgHURsT4iXgbmAVMaHFN3TQGuSx9fB7y/gbG8SUTcCTxTtrujmKcAP4rE74G3SNq3PpG2rz6MvbEAAAPQSURBVIP4OzIFmBcRf42Ih4F1JH/XGiYiNkbEPenj50gWdhpMTj6DCvF3pBk/g4iI59PNndKfAI4FFqT7yz+D4mezADhOkuoUbo/lKREMBh4r2W6j8l+uZhHA7ZKWS5qZ7ntrRGxMHz8BvLUxoXVJRzHn6XOZlXadzC3pjmvq+NMuhvEk/yPN3WdQFj/k6DOQ1EfSCuBJ4JckLZVnI1l9EbaPc9t7SI9vBgbUN+Luy1MiyKujIuIw4ETgLEnvKj0YSVsyV/fw5jFm4F+BtwHjgI3AtxsbTuck7Q78FPhsRGwpPZaHz6Cd+HP1GUTEaxExDhhC0kIZ1eCQaiZPiWADMLRke0i6r6lFxIb095PAjSR/of5SbLqnv59sXISZdRRzLj6XiPhL+g/7deAq3uh6aMr4Je1E8iX6HxHxs3R3bj6D9uLP22dQFBHPAouBvyPpdivO0VYa57b3kB7fE9hU51C7LU+JYCkwMq3a70xSkGnqxe4l9ZPUv/gYOAG4jyTuj6anfRS4qTERdklHMS8EPpLeufJOYHNJ90XTKOsz/3uSzwGS+Keld32MAEYCf6h3fKXSvuV/A+6PiMtLDuXiM+go/px9BoMkvSV9vCvwbpJax2Jganpa+WdQ/GymAr9OW2350OhqdVd+SO6OeIikr+6LjY4nQ7wHkNwN8UdgdTFmkr7DXwFrgTuAvRsda1ncPyFpur9C0g96Rkcxk9xdMTv9TFYBhSaN/8dpfCtJ/tHuW3L+F9P4HwRObIL4jyLp9lkJrEh/TsrLZ1Ah/jx9BocC96ax3gdcnO4/gCRJrQNuAHZJ9/dNt9elxw9o9Hvoyo+nmDAza3F56hoyM7MacCIwM2txTgRmZi3OicDMrMU5EZiZtTgnArMKJL1F0qfTx/tJWtDZc8zyxrePmlWQzpVzc0Qc0uBQzGpmx85PMWtplwFvSycfWwscHBGHSDqdZObJfiQjYb9FMj36PwB/BU6KiGckvY1ksNcgYCvwiYh4oP5vw6xj7hoyq+wC4E+RTD52ftmxQ4BTgMOBS4GtETEeuAv4SHrOHODsiHgH8DngyrpEbdYFbhGYdd/iSObbf07SZuDn6f5VwKHp7Jv/C7ihZGr6XeofplllTgRm3ffXksevl2y/TvJvaweS+evH1Tsws65w15BZZc+RLLfYZZHMwf+wpA/CtrWF/7aawZlVgxOBWQURsQn4naT7gG924xL/GzhDUnEG2rwur2q9mG8fNTNrcW4RmJm1OCcCM7MW50RgZtbinAjMzFqcE4GZWYtzIjAza3FOBGZmLe7/A/P+czNnjJbfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "S1 = data[data.Stage_group == 1]\n",
    "km1 = KM()\n",
    "km1.fit(S1.loc[:, 'Time'], event_observed = S1.loc[:, 'Event'], label = 'Stage III')\n",
    "\n",
    "S2 = data[data.Stage_group == 2]\n",
    "km2 = KM()\n",
    "km2.fit(S2.loc[:, \"Time\"], event_observed = S2.loc[:, 'Event'], label = 'Stage IV')\n",
    "\n",
    "ax = km1.plot(ci_show=False)\n",
    "km2.plot(ax = ax, ci_show=False)\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('Survival probability estimate')\n",
    "plt.savefig('two_km_curves', dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "M4DwaOVEs19Q"
   },
   "source": [
    "Let's compare the survival functions at 90, 180, 270, and 360 days"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "11dhdsUOtEqe"
   },
   "outputs": [],
   "source": [
    "survivals = pd.DataFrame([90, 180, 270, 360], columns = ['time'])\n",
    "survivals.loc[:, 'Group 1'] = km1.survival_function_at_times(survivals['time']).values\n",
    "survivals.loc[:, 'Group 2'] = km2.survival_function_at_times(survivals['time']).values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 173
    },
    "colab_type": "code",
    "id": "-zRlM1SAtYdl",
    "outputId": "3642dc3e-01b0-4e96-e91f-8e39c6c0e3e5"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>time</th>\n",
       "      <th>Group 1</th>\n",
       "      <th>Group 2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>90</td>\n",
       "      <td>0.736842</td>\n",
       "      <td>0.424529</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>180</td>\n",
       "      <td>0.680162</td>\n",
       "      <td>0.254066</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>270</td>\n",
       "      <td>0.524696</td>\n",
       "      <td>0.195436</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>360</td>\n",
       "      <td>0.524696</td>\n",
       "      <td>0.195436</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   time   Group 1   Group 2\n",
       "0    90  0.736842  0.424529\n",
       "1   180  0.680162  0.254066\n",
       "2   270  0.524696  0.195436\n",
       "3   360  0.524696  0.195436"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "survivals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "RA3amMk__J6e"
   },
   "source": [
    "This makes clear the difference in survival between the Stage III and IV cancer groups in the dataset. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3VoOQREQoXny"
   },
   "source": [
    "<a name='5-1'></a>\n",
    "## 5.1 Bonus: Log-Rank Test\n",
    "\n",
    "To say whether there is a statistical difference between the survival curves we can run the log-rank test. This test tells us the probability that we could observe this data if the two curves were the same. The derivation of the log-rank test is somewhat complicated, but luckily `lifelines` has a simple function to compute it. \n",
    "\n",
    "Run the next cell to compute a p-value using `lifelines.statistics.logrank_test`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "_7-QIy8ovsgC",
    "outputId": "c7582d94-4c42-4cae-d83d-72a8873fe985"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.009588929834755544"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def logrank_p_value(group_1_data, group_2_data):\n",
    "    result = logrank_test(group_1_data.Time, group_2_data.Time,\n",
    "                          group_1_data.Event, group_2_data.Event)\n",
    "    return result.p_value\n",
    "\n",
    "logrank_p_value(S1, S2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "nUbv_csdJRSw"
   },
   "source": [
    "If everything is correct, you should see a p value of less than `0.05`, which indicates that the difference in the curves is indeed statistically significant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Congratulations!\n",
    "\n",
    "You've completed the third assignment of Course 2. You've learned about the Kaplan Meier estimator, a fundamental non-parametric estimator in survival analysis. Next week we'll learn how to take into account patient covariates in our survival estimates!"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [
    "3VoOQREQoXny"
   ],
   "include_colab_link": true,
   "name": "C2M3_Assignment.ipynb",
   "provenance": [],
   "toc_visible": true
  },
  "coursera": {
   "schema_names": [
    "AI4MC2-3"
   ]
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}