{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Suppose we generate data from a mixture of exponentials with no censoring (this could be extended to censoring, but for our application we don't observe censoring). Furthermore, the data is not observed directly, but our instrument has binned the data into integer buckets. For example, the actual value is 4.3, but our instrument's resolution is too poor to measure the decimal point. So 4.3 goes into the \"4.x\" bucket. \n", "\n", "This model is easy to handle in *lifelines*. Instead of worrying about the binning, we can treat the system as an interval-censoring problem. That is, if an observation landed in the $i$th bin, then we know the _true_ data point occured somewhere between the $i$th and $i+1$th bin. This is precisely interval censoring. \n", "\n", "We can use *lifelines* custom model semantics to create a mixture model as well. The true model is:\n", "\n", "$$S(t) = p \\exp\\left(-\\frac{t}{\\lambda_1}\\right) + (1-p)\\exp\\left(-\\frac{t}{\\lambda_2}\\right)$$\n", "\n", "Therefore the cumulative hazard is:\n", "\n", "$$H(t) = -\\log(S(t)) = -\\log\\left(p \\exp\\left(-\\frac{t}{\\lambda_1}\\right) + (1-p)\\exp\\left(-\\frac{t}{\\lambda_2}\\right)\\right) $$" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "from numpy.random import exponential\n", "from lifelines.fitters import ParametricUnivariateFitter\n", "import autograd.numpy as np\n", "np.random.seed(10)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "T = [exponential(20) for _ in range(10000)] + [exponential(40) for _ in range(500)]\n", "counts_obs = np.bincount(T)\n", "T_obs = np.arange(np.amax(T))" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0, 100)" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwoAAAE/CAYAAADmEDQWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X+0XWV95/H3pzcEG4GIQlWSaAKGtkErtCloO1VbUWOxxGlRQ7Vih05KK9VWXTa2Fm2sHaQO1Y60JS2ZWixGRIfeahxqRe1YBXMRKiaWRYhIElGUnyIIBL7zx96ph7PvzT25P3Jv4P1a666c/fzY59n37rVzPmc/zzmpKiRJkiSp1w/N9AAkSZIkzT4GBUmSJEkdBgVJkiRJHQYFSZIkSR0GBUmSJEkdBgVJkiRJHQYFSY84STYned5Mj2MmJfmvSbYnuTvJcXvZ9zNJfmO6xjYdkvxdkj8ZsO3iJJVkznSPq+c5H/XnpKT9j0FB0n4lyY1JTuwre02Sz+3erqpjquoz4+xnn79Y3MfeDZxZVQdV1dUzPZhHk9FCyyDnpCTNNgYFSZoGsyCAPBXYPMNjGNMs+P1MyP46bkmaCIOCpEec3rsOSY5PMpLkriTfSnJu2+xf23/vaKfnPDvJDyV5a5KvJ7klyd8nmd+z31e3dbcm+aO+53l7kkuSfCDJXcBr2uf+QpI7ktyc5H1J5vbsr5L8dpLrk3w3yTuSHJXk8+14L+5t33eMo441yYFJ7gaGgH9PcsMY/X8myaYkd7b//kxfk6OSfLEdxz8meXzb7zHtMd7aHtemJE9s6+YnuaA91p1J/iTJUFv3miT/luTPk9wKvKPt//SeMR2e5N4kP9JuvyTJNW27zyf5iZ62xyX5Uvt7+xDwmD2cD0NJ3p3kO0m2ASf11T/sLlX7t/xA+3j3nafTk9wEXN6WfzjJN9vf378mOaYtXw28Enhze179U/9ztH+j9yT5RvvzniQHtnXPS7IjyRvbv+vNSX69Z2y/mGRLe9w7k7xprOOWpMkyKEh6pHsv8N6qOgQ4Cri4LX9O++/j2uk5XwBe0/78PHAkcBDwPoAky4C/pHkR+GRgPrCg77lWApcAjwP+AXgQ+D3gMODZwPOB3+7r8yLgp4BnAW8G1gGvAhYBTwdOHeO4Rh1rVd1XVQe1bZ5ZVUf1d2xf9H8c+AvgCcC5wMeTPKGn2auB/9Ye6662LcBp7bEvavueAdzb1v1d2/ZpwHHAC4HetQ4nANuAJwJrgY/2Hd/Lgc9W1S1p1lWsB36zfZ7zgeH2RfZc4FLgQuDxwIeBXxnj9wTw34GXtGNaDpyyh7ZjeS7w4zR/L4BPAEuBHwG+RPP3pqrWtY/Pac+rXxplX39I8/c+FngmcDzw1p76J/GD8+t04Lwkh7Z1FwC/WVUH05wfl0/gWCRpIAYFSfujS9t3me9IcgfNC/ixPAA8LclhVXV3VV2xh7avBM6tqm1VdTfwFmBVmukmpwD/VFWfq6r7gbOA6uv/haq6tKoeqqp7q+qqqrqiqnZV1Y00L3af29fnnKq6q6o2A18B/rl9/jtpXoyOtRB5T2Mdz0nA9VV1YTu2DwL/AfS+qL2wqr5SVd8D/gh4eXt34AGaF+5Pq6oH22O8q72r8IvA71bV96rqFuDPgVU9+/xGVf2v9jnvBS7qq//VtgxgNXB+VV3ZPs/7gftoXmA/CzgAeE9VPVBVlwCb9nC8L2/bbq+q24D/McDvqN/b2+O6F6Cq1lfVd6vqPuDtwDPTc/dpHK8E1lbVLVX1beCPgV/rqX+grX+gqjYCdwM/2lO3LMkhVXV7VX1pAsciSQMxKEjaH720qh63+4fuu/S9TgeOBv6jnSbzkj20PQL4es/214E5NO+AHwFs311RVfcAt/b13967keToJB9rp6jcBfwpzd2FXt/qeXzvKNsHMbo9jXU8/X139++9Q7K9r+4AmrFfCFwGbGinzZyT5ACaNREHADf3BLjzad5xH22fAJ8G5iU5IclimnfY/09b91TgjX2BcFE79iOAnVXVG9T6j6f/ePuPZ2/9Z/92KtPZSW5o/643tlX9f9s9jaf/b3dEz/atVbWrZ/sefnAe/ApNIPt6ks8mefZeHIMk7RWDgqRHtKq6vqpOpXnB+i7gkiSPpXs3AOAbNC9Qd3sKzVSabwE3Awt3VyT5YZp31h/2dH3bf0XzTv3SdurTHwCZ+NEMPNa97bu7/86e7UV9dQ8A32nf5f7jqloG/AzNlJ5X07yQvg84rCfEHVJVx/Ts52G/n6p6kGYq2Kntz8eq6rtt9Xbgnb2BsKrmtXc/bgYWJOn9XT5lD8d78yjH0+t7wLye7SeNso/esf8qzTSzE2mmCC1uyzNK29GM9rf7xjh9mh1XbaqqlTTn86X8YCqdJE05g4KkR7Qkr0pyeFU9BNzRFj8EfLv998ie5h8Efi/JkiQH0dwB+FD77u4lwC+1i4Dn0kw3Ge9F/8HAXcDdSX4M+K2pOq5xxjqejcDRSX41yZwkrwCWAR/rafOqJMuSzKNZT3BJVT2Y5OeTPKOdhnQXTYB4qKpuBv4Z+J9JDkmz2PqoJP1TrfpdBLyCZjrORT3lfwOc0d5tSJLHJjkpycHAF2hC0euSHJDkl2nm+Y/l4rbtwnau/5q++mtopm0dkGSQNQwH04SiW2kCxp/21X+Lh59X/T4IvDXN4u3DaKaxfWCc5yTJ3CSvTDK/qh6g+f0/NF4/SZoog4KkR7oVwOY0nwT0XmBVu37gHuCdwL+1U1ueRbN49kKaT0T6GvB94HcA2jUEvwNsoHmH+m7gFpoXjGN5E827z9+leeH7oSk8rjHHOp6qupXmTsAbaV7svhl4SVV9p6fZhTSLk79J84lCr2vLn0QTmu4Cvgp8tm0LzZ2FucAW4Pa23ZPHGcuVNO/oH0GzJmN3+QjNIuT3tfvaSrN4m3aNyC+327fRBI2P7uFp/oZmutS/0yw87m/7RzQL3W+nWS9wEXv29zTThXbSHGv/upcLaNYR3JHk0lH6/wkwAnwZuLYd00BfFkezluHGdsrTGTQBS5KmRR4+xVOSNIj2Xfw7aKYVfW2mxyNJ0lTzjoIkDSjJLyWZ165xeDfNu8E3zuyoJEmaHgYFSRrcSppFp9+g+Qz9VeVtWUnSI5RTjyRJkiR1eEdBkiRJUodBQZIkSVLHnJkeQL/DDjusFi9ePNPDkCRJkvZrV1111Xeq6vCJ9p91QWHx4sWMjIzM9DAkSZKk/VqSr0+mv1OPJEmSJHUYFCRJkiR1GBQkSZIkdRgUJEmSJHUYFCRJkiR1GBQkSZIkdRgUJEmSJHUYFCRJkiR1GBQkSZIkdRgUJEmSJHUYFCRJkiR1zJnpAfS7duedLF7z8THrbzz7pH04GkmSJOnRyTsKkiRJkjoGCgpJViS5LsnWJGv20O5XklSS5T1lb2n7XZfkRVMxaEmSJEnTa9ypR0mGgPOAFwA7gE1JhqtqS1+7g4HXA1f2lC0DVgHHAEcA/5Lk6Kp6cOoOQZIkSdJUG+SOwvHA1qraVlX3AxuAlaO0ewfwLuD7PWUrgQ1VdV9VfQ3Y2u5PkiRJ0iw2SFBYAGzv2d7Rlv2nJD8JLKqq/lXI4/Zt+69OMpJk5MF77hxo4JIkSZKmz6QXMyf5IeBc4I0T3UdVrauq5VW1fGje/MkOSZIkSdIkDfLxqDuBRT3bC9uy3Q4Gng58JgnAk4DhJCcP0FeSJEnSLDTIHYVNwNIkS5LMpVmcPLy7sqrurKrDqmpxVS0GrgBOrqqRtt2qJAcmWQIsBb445UchSZIkaUqNe0ehqnYlORO4DBgC1lfV5iRrgZGqGt5D381JLga2ALuA1/qJR5IkSdLsN9A3M1fVRmBjX9lZY7R9Xt/2O4F3TnB8kiRJkmaA38wsSZIkqcOgIEmSJKnDoCBJkiSpw6AgSZIkqcOgIEmSJKnDoCBJkiSpw6AgSZIkqcOgIEmSJKnDoCBJkiSpw6AgSZIkqcOgIEmSJKnDoCBJkiSpw6AgSZIkqcOgIEmSJKnDoCBJkiSpw6AgSZIkqcOgIEmSJKnDoCBJkiSpw6AgSZIkqcOgIEmSJKnDoCBJkiSpw6AgSZIkqWOgoJBkRZLrkmxNsmaU+jOSXJvkmiSfS7KsLV+c5N62/Jokfz3VByBJkiRp6s0Zr0GSIeA84AXADmBTkuGq2tLT7KKq+uu2/cnAucCKtu6Gqjp2aoctSZIkaToNckfheGBrVW2rqvuBDcDK3gZVdVfP5mOBmrohSpIkSdrXBgkKC4DtPds72rKHSfLaJDcA5wCv66lakuTqJJ9N8nOTGq0kSZKkfWLKFjNX1XlVdRTw+8Bb2+KbgadU1XHAG4CLkhzS3zfJ6iQjSUYevOfOqRqSJEmSpAkaJCjsBBb1bC9sy8ayAXgpQFXdV1W3to+vAm4Aju7vUFXrqmp5VS0fmjd/0LFLkiRJmiaDBIVNwNIkS5LMBVYBw70Nkizt2TwJuL4tP7xdDE2SI4GlwLapGLgkSZKk6TPupx5V1a4kZwKXAUPA+qranGQtMFJVw8CZSU4EHgBuB05ruz8HWJvkAeAh4Iyqum06DkSSJEnS1Bk3KABU1UZgY1/ZWT2PXz9Gv48AH5nMACVJkiTte34zsyRJkqQOg4IkSZKkDoOCJEmSpA6DgiRJkqQOg4IkSZKkDoOCJEmSpA6DgiRJkqQOg4IkSZKkDoOCJEmSpA6DgiRJkqQOg4IkSZKkDoOCJEmSpA6DgiRJkqQOg4IkSZKkDoOCJEmSpA6DgiRJkqQOg4IkSZKkDoOCJEmSpA6DgiRJkqQOg4IkSZKkDoOCJEmSpA6DgiRJkqQOg4IkSZKkjoGCQpIVSa5LsjXJmlHqz0hybZJrknwuybKeure0/a5L8qKpHLwkSZKk6TFuUEgyBJwHvBhYBpzaGwRaF1XVM6rqWOAc4Ny27zJgFXAMsAL4y3Z/kiRJkmaxQe4oHA9sraptVXU/sAFY2dugqu7q2XwsUO3jlcCGqrqvqr4GbG33J0mSJGkWmzNAmwXA9p7tHcAJ/Y2SvBZ4AzAX+IWevlf09V0woZFKkiRJ2membDFzVZ1XVUcBvw+8dW/6JlmdZCTJyIP33DlVQ5IkSZI0QYMEhZ3Aop7thW3ZWDYAL92bvlW1rqqWV9XyoXnzBxiSJEmSpOk0SFDYBCxNsiTJXJrFycO9DZIs7dk8Cbi+fTwMrEpyYJIlwFLgi5MftiRJkqTpNO4aharaleRM4DJgCFhfVZuTrAVGqmoYODPJicADwO3AaW3fzUkuBrYAu4DXVtWD03QskiRJkqZIqmr8VvvQgU9eWk8+7T1j1t949kn7cDSSJEnS/inJVVW1fKL9/WZmSZIkSR0GBUmSJEkdBgVJkiRJHQYFSZIkSR0GBUmSJEkdBgVJkiRJHQYFSZIkSR0GBUmSJEkdBgVJkiRJHQYFSZIkSR0GBUmSJEkdBgVJkiRJHQYFSZIkSR0GBUmSJEkdBgVJkiRJHQYFSZIkSR0GBUmSJEkdBgVJkiRJHQYFSZIkSR0GBUmSJEkdBgVJkiRJHQYFSZIkSR0GBUmSJEkdAwWFJCuSXJdka5I1o9S/IcmWJF9O8qkkT+2pezDJNe3P8FQOXpIkSdL0mDNegyRDwHnAC4AdwKYkw1W1pafZ1cDyqronyW8B5wCvaOvurapjp3jckiRJkqbRIHcUjge2VtW2qrof2ACs7G1QVZ+uqnvazSuAhVM7TEmSJEn70iBBYQGwvWd7R1s2ltOBT/RsPybJSJIrkrx0tA5JVrdtRh68584BhiRJkiRpOo079WhvJHkVsBx4bk/xU6tqZ5IjgcuTXFtVN/T2q6p1wDqAA5+8tKZyTJIkSZL23iB3FHYCi3q2F7ZlD5PkROAPgZOr6r7d5VW1s/13G/AZ4LhJjFeSJEnSPjBIUNgELE2yJMlcYBXwsE8vSnIccD5NSLilp/zQJAe2jw8DfhboXQQtSZIkaRYad+pRVe1KciZwGTAErK+qzUnWAiNVNQz8GXAQ8OEkADdV1cnAjwPnJ3mIJpSc3fdpSZIkSZJmoYHWKFTVRmBjX9lZPY9PHKPf54FnTGaAkiRJkvY9v5lZkiRJUodBQZIkSVKHQUGSJElSh0FBkiRJUodBQZIkSVKHQUGSJElSh0FBkiRJUodBQZIkSVKHQUGSJElSh0FBkiRJUodBQZIkSVKHQUGSJElSh0FBkiRJUodBQZIkSVKHQUGSJElSh0FBkiRJUodBQZIkSVKHQUGSJElSh0FBkiRJUodBQZIkSVKHQUGSJElSh0FBkiRJUsdAQSHJiiTXJdmaZM0o9W9IsiXJl5N8KslTe+pOS3J9+3PaVA5ekiRJ0vQYNygkGQLOA14MLANOTbKsr9nVwPKq+gngEuCctu/jgbcBJwDHA29LcujUDV+SJEnSdBjkjsLxwNaq2lZV9wMbgJW9Darq01V1T7t5BbCwffwi4JNVdVtV3Q58ElgxNUOXJEmSNF0GCQoLgO092zvasrGcDnxign0lSZIkzQJzpnJnSV4FLAeeu5f9VgOrAYYOOXwqhyRJkiRpAga5o7ATWNSzvbAte5gkJwJ/CJxcVfftTd+qWldVy6tq+dC8+YOOXZIkSdI0GSQobAKWJlmSZC6wChjubZDkOOB8mpBwS0/VZcALkxzaLmJ+YVsmSZIkaRYbd+pRVe1KcibNC/whYH1VbU6yFhipqmHgz4CDgA8nAbipqk6uqtuSvIMmbACsrarbpuVIJEmSJE2ZgdYoVNVGYGNf2Vk9j0/cQ9/1wPqJDlCSJEnSvuc3M0uSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6hgoKCRZkeS6JFuTrBml/jlJvpRkV5JT+uoeTHJN+zM8VQOXJEmSNH3mjNcgyRBwHvACYAewKclwVW3paXYT8BrgTaPs4t6qOnYKxipJkiRpHxk3KADHA1urahtAkg3ASuA/g0JV3djWPTQNY5QkSZK0jw0y9WgBsL1ne0dbNqjHJBlJckWSl+7V6CRJkiTNiEHuKEzWU6tqZ5IjgcuTXFtVN/Q2SLIaWA0wdMjh+2BIkiRJkvZkkDsKO4FFPdsL27KBVNXO9t9twGeA40Zps66qllfV8qF58wfdtSRJkqRpMkhQ2AQsTbIkyVxgFTDQpxclOTTJge3jw4CfpWdtgyRJkqTZadygUFW7gDOBy4CvAhdX1eYka5OcDJDkp5PsAF4GnJ9kc9v9x4GRJP8OfBo4u+/TkiRJkiTNQgOtUaiqjcDGvrKzeh5vopmS1N/v88AzJjlGSZIkSfuY38wsSZIkqcOgIEmSJKnDoCBJkiSpw6AgSZIkqcOgIEmSJKnDoCBJkiSpw6AgSZIkqcOgIEmSJKnDoCBJkiSpw6AgSZIkqcOgIEmSJKnDoCBJkiSpw6AgSZIkqcOgIEmSJKnDoCBJkiSpw6AgSZIkqcOgIEmSJKljzkwPYG8tXvPxMetuPPukfTgSSZIk6ZHLOwqSJEmSOgwKkiRJkjoMCpIkSZI6DAqSJEmSOgwKkiRJkjoGCgpJViS5LsnWJGtGqX9Oki8l2ZXklL6605Jc3/6cNlUDlyRJkjR9xg0KSYaA84AXA8uAU5Ms62t2E/Aa4KK+vo8H3gacABwPvC3JoZMftiRJkqTpNMgdheOBrVW1raruBzYAK3sbVNWNVfVl4KG+vi8CPllVt1XV7cAngRVTMG5JkiRJ02iQoLAA2N6zvaMtG8RAfZOsTjKSZOTBe+4ccNeSJEmSpsusWMxcVeuqanlVLR+aN3+mhyNJkiQ96g0SFHYCi3q2F7Zlg5hMX0mSJEkzZJCgsAlYmmRJkrnAKmB4wP1fBrwwyaHtIuYXtmWSJEmSZrFxg0JV7QLOpHmB/1Xg4qranGRtkpMBkvx0kh3Ay4Dzk2xu+94GvIMmbGwC1rZlkiRJkmaxOYM0qqqNwMa+srN6Hm+imVY0Wt/1wPpJjFGSJEnSPjYrFjNLkiRJml0MCpIkSZI6DAqSJEmSOgwKkiRJkjoMCpIkSZI6DAqSJEmSOgwKkiRJkjoMCpIkSZI6DAqSJEmSOgwKkiRJkjoMCpIkSZI6DAqSJEmSOgwKkiRJkjoMCpIkSZI6DAqSJEmSOgwKkiRJkjoMCpIkSZI6DAqSJEmSOgwKkiRJkjoMCpIkSZI6DAqSJEmSOgwKkiRJkjoGCgpJViS5LsnWJGtGqT8wyYfa+iuTLG7LFye5N8k17c9fT+3wJUmSJE2HOeM1SDIEnAe8ANgBbEoyXFVbepqdDtxeVU9Lsgp4F/CKtu6Gqjp2isctSZIkaRoNckfheGBrVW2rqvuBDcDKvjYrgfe3jy8Bnp8kUzdMSZIkSfvSIEFhAbC9Z3tHWzZqm6raBdwJPKGtW5Lk6iSfTfJzkxyvJEmSpH1g3KlHk3Qz8JSqujXJTwGXJjmmqu7qbZRkNbAaYOiQw6d5SJIkSZLGM8gdhZ3Aop7thW3ZqG2SzAHmA7dW1X1VdStAVV0F3AAc3f8EVbWuqpZX1fKhefP3/igkSZIkTalBgsImYGmSJUnmAquA4b42w8Bp7eNTgMurqpIc3i6GJsmRwFJg29QMXZIkSdJ0GXfqUVXtSnImcBkwBKyvqs1J1gIjVTUMXABcmGQrcBtNmAB4DrA2yQPAQ8AZVXXbdByIJEmSpKkz0BqFqtoIbOwrO6vn8feBl43S7yPARyY5RkmSJEn7mN/MLEmSJKnDoCBJkiSpw6AgSZIkqWO6v0dhn1q85uMT6nfj2SdNeL/j9ZUkSZL2R95RkCRJktRhUJAkSZLUYVCQJEmS1GFQkCRJktTxiFrMPFETXQQtSZIkPVJ5R0GSJElSh0FBkiRJUodBQZIkSVKHQUGSJElSh0FBkiRJUodBQZIkSVKHQUGSJElSh0FBkiRJUodBQZIkSVKHQUGSJElSx5yZHsD+bvGaj49Zd+PZJ02o33gmut899RvPdO1XkiRJs5N3FCRJkiR1eEdhGk3mrsFM7Hd/Ml2/A++OzE6Pljta+9Nx7k9jlSRNjHcUJEmSJHUMFBSSrEhyXZKtSdaMUn9gkg+19VcmWdxT95a2/LokL5q6oUuSJEmaLuNOPUoyBJwHvADYAWxKMlxVW3qanQ7cXlVPS7IKeBfwiiTLgFXAMcARwL8kObqqHpzqA9H4xpuuMxPTBWbbNKqZGM90LTKf6HNO10L7PZmJaXqz7fc+meecbVN9ZuIcGs90fPDEo+Xvub/xdzs5M/H7eyRdqydqNl43YbA7CscDW6tqW1XdD2wAVva1WQm8v318CfD8JGnLN1TVfVX1NWBruz9JkiRJs9ggQWEBsL1ne0dbNmqbqtoF3Ak8YcC+kiRJkmaZVNWeGySnACuq6jfa7V8DTqiqM3vafKVts6PdvgE4AXg7cEVVfaAtvwD4RFVd0vccq4HV7ebTga9M/tD0KHcY8J2ZHoT2a55DmizPIU2W55Am60er6uCJdh7k41F3Aot6the2ZaO12ZFkDjAfuHXAvlTVOmAdQJKRqlo+6AFIo/E80mR5DmmyPIc0WZ5DmqwkI5PpP8jUo03A0iRLksylWZw83NdmGDitfXwKcHk1tyqGgVXtpyItAZYCX5zMgCVJkiRNv3HvKFTVriRnApcBQ8D6qtqcZC0wUlXDwAXAhUm2ArfRhAnadhcDW4BdwGv9xCNJkiRp9hvom5mraiOwsa/srJ7H3wdeNkbfdwLv3IsxrduLttJYPI80WZ5DmizPIU2W55Ama1Ln0LiLmSVJkiQ9+gz0zcySJEmSHl1mVVBIsiLJdUm2Jlkz0+PR7JdkUZJPJ9mSZHOS17flj0/yySTXt/8eOtNj1eyWZCjJ1Uk+1m4vSXJlez36UPthDtKokjwuySVJ/iPJV5M82+uQ9kaS32v/H/tKkg8meYzXIY0nyfokt7RfVbC7bNRrTxp/0Z5PX07yk+Ptf9YEhSRDwHnAi4FlwKlJls3sqLQf2AW8saqWAc8CXtueN2uAT1XVUuBT7ba0J68Hvtqz/S7gz6vqacDtwOkzMirtL94L/N+q+jHgmTTnktchDSTJAuB1wPKqejrNh8eswuuQxvd3wIq+srGuPS+m+QTSpTTfX/ZX4+181gQF4Hhga1Vtq6r7gQ3Ayhkek2a5qrq5qr7UPv4uzX/OC2jOnfe3zd4PvHRmRqj9QZKFwEnA37bbAX4B2P3lkJ5DGlOS+cBzaD4BkKq6v6ruwOuQ9s4c4Ifb76OaB9yM1yGNo6r+leYTR3uNde1ZCfx9Na4AHpfkyXva/2wKCguA7T3bO9oyaSBJFgPHAVcCT6yqm9uqbwJPnKFhaf/wHuDNwEPt9hOAO6pqV7vt9Uh7sgT4NvC/2+lrf5vksXgd0oCqaifwbuAmmoBwJ3AVXoc0MWNde/b6tfZsCgrShCU5CPgI8LtVdVdvXfvlf368l0aV5CXALVV11UyPRfutOcBPAn9VVccB36NvmpHXIe1JO4d8JU3oPAJ4LN3pJNJem+y1ZzYFhZ3Aop7thW2ZtEdJDqAJCf9QVR9ti7+1+3Za++8tMzU+zXo/C5yc5EaaKY+/QDPf/HHtFADweqQ92wHsqKor2+1LaIKD1yEN6kTga1X17ap6APgozbXJ65AmYqxrz16/1p5NQWETsLRd4T+XZhHP8AyPSbNcO5f8AuCrVXVuT9UwcFr7+DTgH/f12LR/qKq3VNXCqlpMc925vKpeCXwaOKVt5jmkMVXVN4HtSX60LXo+sAWvQxrcTcCzksxr/1/bfQ55HdJEjHXtGQZe3X760bOAO3umKI1qVn3hWpJfpJkrPASsb7/VWRpTkv8C/D/gWn4wv/wPaNYpXAw8Bfg68PKq6l/sIz1MkucBb6qqlyQ5kuYOw+OBq4FXVdV9Mzk+zV5JjqVZDD8X2Ab8Os2bcV6HNJAkfwy8gube0FnRAAAAdklEQVTT/K4GfoNm/rjXIY0pyQeB5wGHAd8C3gZcyijXnjaEvo9mWts9wK9X1cge9z+bgoIkSZKk2WE2TT2SJEmSNEsYFCRJkiR1GBQkSZIkdRgUJEmSJHUYFCRJkiR1GBQkSZIkdRgUJEmSJHUYFCRJkiR1/H/y/IQ0c48thgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(13,5))\n", "plt.hist(counts_obs, bins=T_obs, density=True)\n", "plt.title(\"Histogram of observed durations\")\n", "plt.xlim(0, 100)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "20.015428571428572\n" ] } ], "source": [ "# To help the model avoid the \"non-identibility\" problem, we can set the *upper* bound of the first lambda to be \n", "# the average of the observed data, and the *lower* bound of the second lambda to the same value. Why? \n", "# We'd like to partition the postive reals into two halves, each containing one of the lambdas. The sample\n", "# mean of the data is v = p * lambda1 + (1-p) * lambda2, which has the property lambda1 < v < lambda2, therefore \n", "# it will partition the space correctly. \n", "mean_obs = np.average(T_obs, weights=counts_obs)\n", "print(mean_obs)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
modellifelines.MixtureExponential
number of observations10500
number of events observed0
log-likelihood-42208.35
hypothesislambda1 != 10.0077, lambda2 != 40.0309, p != 0.5
\n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
coefse(coef)coef lower 95%coef upper 95%zp-log2(p)
lambda119.001.0117.0120.988.88<0.00560.36
lambda236.4610.3216.2456.68-0.350.730.46
p0.910.100.711.114.06<0.00514.30
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "class MixtureExponential(ParametricUnivariateFitter):\n", "\n", " _fitted_parameter_names = ['lambda1', 'lambda2', 'p']\n", " _bounds = [(0, mean_obs), (mean_obs, None), (0, 1)]\n", "\n", " def _cumulative_hazard(self, params, t):\n", " l1, l2, p = params\n", " return -np.log(p * np.exp(-t / l1) + (1-p) * np.exp(-t / l2))\n", "\n", "\n", "model = MixtureExponential()\n", "model.fit_interval_censoring(\n", " lower_bound=T_obs, \n", " upper_bound=T_obs + 1, \n", " weights=counts_obs, \n", " initial_point=np.array([mean_obs / 2, mean_obs * 2, 0.5])\n", ")\n", "\n", "model.print_summary()" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAEyCAYAAADJFbiyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xmc3FWd7//Xt7aupfd9X7IvJCQQCBACBGQRGRFFRUFRr+I43tG5Vx3AEVHuzPycGZ07o+M4w8wo1xEVNwQFkcWEEEAggZCls3Rn7/TenV6qurZvfc/vj+o0CSSk0/vyfj4eeZCuVFed/qbTb875ns/nWMYYREREZPK4JnsAIiIis53CWEREZJIpjEVERCaZwlhERGSSKYxFREQmmcJYRERkkimMRUREJpnCWEREZJIpjEVERCaZZyLfrLCw0NTW1k7kW4qIiEyaLVu2dBpjis70vAkN49raWjZv3jyRbykiIjJpLMs6NJznaZlaRERkkimMRUREJpnCWEREZJJN6D1jERGAZDJJU1MTsVhssociMib8fj+VlZV4vd4Rfb7CWEQmXFNTE1lZWdTW1mJZ1mQPR2RUjDF0dXXR1NREXV3diF5Dy9QiMuFisRgFBQUKYpkRLMuioKBgVCs9CmMRmRQKYplJRvv9rDAWERGZZApjERGRSaYwFpFZybIsbrvttqGPbdumqKiIG264AYBHH32Ub3zjG2/7Gg888ADNzc1jNqaDBw8SCARYsWLF0K8f/vCHY/b6Y+1v//ZvT/r4kksuOePnZGZmjtn7/9M//RMDAwNDH19//fX09PSM+nW3bt3K448/PurXORsKYxGZlUKhEDt27CAajQLw1FNPUVFRMfTn7373u7nrrrve9jVGEsa2bb/tn8+dO5etW7cO/froRz96Vq8/kd4cxi+88MKEvv+bw/jxxx8nNzd31K87VmGcsB1we33Dea5Km0RkUn39Nzupb+4b09dcUp7NvX+y9IzPu/7663nssce4+eab+clPfsKHPvQhnnvuOSAdtJs3b+Zf/uVfuPHGG3nf+97HRz/6Uf793/+djRs3ctNNN7F582ZuvfVWAoEAL774IosXL2bz5s0UFhayefNmvvjFL7Jhwwa+9rWvsW/fPvbv3091dTU/+tGPuOuuu9iwYQPxeJzPfvazfPrTnz7tOA8dOsQ73vEOXnzxRfLz87n88su55557WLBgAddddx3nn38+r776KkuXLuWHP/whwWCQZ555hi9+8YvYts0FF1zA9773PTIyMqitreX222/nN7/5Dclkkp///OcsWrSISCTCn//5n7Njxw6SySRf+9rXuPHGG3nggQd49NFHGRgYYN++fdx00038/d//PXfddRfRaJQVK1awdOlSHnzwQTIzMwmHw4TDYW688UaOHTtGMpnkr//6r7nxxhuH9Xf3D//wD/zsZz8jHo9z00038fWvf51IJMIHPvABmpqaSKVS3HPPPbS1tdHc3My6desoLCxk/fr1Q+cfhMNhrrvuOi666CJeeOEFLrjgAj7+8Y9z77330t7ezoMPPsiFF17Iyy+/zOc//3lisRiBQIAf/OAH1NXV8dWvfpVoNMqmTZu4++67ueGGG055bU7HTjk0HYtS39KLKyM4rKUAzYxFZNa65ZZb+OlPf0osFmPbtm2sXr36lM+7//77ue+++3juuef41re+xXe+8x1uvvlmVq1axYMPPsjWrVsJBAJv+1719fU8/fTT/OQnP+G//uu/yMnJ4ZVXXuGVV17hP/7jPzhw4AAA+/btO2mZ+rnnnqOmpoY777yTz3zmM3zrW99iyZIlXHPNNQDs2bOHP/uzP2PXrl1kZ2fzr//6r8RiMT72sY/x0EMPsX37dmzb5nvf+97QWAoLC3n11Vf5zGc+wze/+U0A/uZv/oYrr7ySl19+mfXr1/OlL32JSCQCpGeKx1/roYce4siRI3zjG98gEAiwdetWHnzwwZO+Vr/fz8MPP8yrr77K+vXr+cIXvoAx5ox/H08++SQNDQ28/PLLbN26lS1btrBx40aeeOIJysvLef3119mxYwfXXXcdn/vc5ygvL2f9+vWsX7/+La/V2NjIF77wBXbv3s3u3bv58Y9/zKZNm/jmN785NKNftGgRzz33HK+99hr33XcfX/7yl/H5fNx333188IMfZOvWrXzwgx9822tzImOguWeA9Xvaeb2ph5DPy7C+cDQzFpFJNpwZ7HhZvnw5Bw8e5Cc/+QnXX3/9aZ9XUlLCfffdx7p163j44YfJz88/6/d697vfPRTYTz75JNu2beMXv/gFAL29vTQ0NLBgwYKhZeo3++QnP8nPf/5z/u3f/u2kP6+qqmLNmjUA3HbbbXz729/m6quvpq6ujgULFgBw++23893vfpe/+Iu/AOC9730vAOeffz6/+tWvhsb06KOPDoVzLBbj8OHDAFx11VXk5OQAsGTJEg4dOkRVVdVpv1ZjDF/+8pfZuHEjLpeLo0eP0tbWRmlp6dteoyeffJInn3ySlStXAhAOh2loaGDt2rV84Qtf4M477+SGG25g7dq1b/s6AHV1dSxbtgyApUuXctVVV2FZFsuWLePgwYNA+rrffvvtNDQ0YFkWyWTytOM61bVZvHjx0NdrO4a4neLlg8fI8Xspzjq7TlwKYxGZ1d797ncPLSd3dXWd9nnbt2+noKDgbe8RezweHMcBeEsDiFAoNPR7Ywzf+c53uPbaa096zvGQOJWBgQGampqAdEhlZWUBb61vHU69a0ZGBgBut3voHrYxhl/+8pcsXLjwpOe+9NJLQ89/8+eczoMPPkhHRwdbtmzB6/VSW1s7rIYYxhjuvvvuUy7Zv/rqqzz++ON85Stf4aqrruKrX/3qsL5GAJfLNfSxy+UaGv8999wz9D9YBw8e5IorrjjtuE51bSC9JB1Lpkg66QlwSZb/jF/nqWiZWkRmtU984hPce++9Q7OoU3n55Zf53e9+x2uvvcY3v/nNoSXlrKws+vv7h55XW1vLli1bAPjlL3952te79tpr+d73vjc0E9u7d+8plz1PdOedd3Lrrbdy33338alPfWro8cOHD/Piiy8C8OMf/5hLL72UhQsXcvDgQRobGwH47//+by6//PK3ff1rr72W73znO0Orqq+99trbPh/A6/WecjbZ29tLcXExXq+X9evXc+jQsI705dprr+X73/8+4XAYgKNHj9Le3k5zczPBYJDbbruNL33pS7z66qvAW6//2ert7R3atPfAAw8MPf7m1z3VtUk5hkjcpj9ukzLgdbkYTdsPhbGIzGqVlZV87nOfO+2fx+NxPvWpT/H973+f8vJyvvWtb/GJT3wCYwwf+9jH+NM//VNWrFhBNBrl3nvv5fOf/zyrVq3C7Xaf9jU/+clPsmTJEs477zzOOeccPv3pTw/N1t58z/jb3/42zz77LK+88spQIPt8Pn7wgx8AsHDhQr773e+yePFijh07xmc+8xn8fj8/+MEPeP/738+yZctwuVz86Z/+6dteh3vuuYdkMsny5ctZunQp99xzzxmv3R133MHy5cu59dZbT3r81ltvZfPmzSxbtowf/vCHLFq06IyvBXDNNdfw4Q9/mIsvvphly5Zx880309/fz/bt27nwwgtZsWIFX//61/nKV74y9P7XXXcd69atG9brv9lf/uVfcvfdd7Ny5cqTZvvr1q2jvr6eFStW8NBDD73l2nz5r75CfyyJnXLwuizcY9BNzhrmveUxsWrVKrN58+YJez8RmZp27do1dL9NRu7gwYPccMMN7NixY7KHMuM5xpCw00vSAB6XBW+aCzfu3U00VHbSY5cum7M/NdA790yvr3vGIiIip2FOCGEDuF0W1qgWpE9NYSwiMk3V1tZOy1nx9u3b+chHPnLSYxkZGbz00kuTNKK3MsaQTDnEkg6OMeMWwscpjEVkUhhjdHLTLLVs2bJTlm9NBcfLlGKJFLYxuC0Lj+vM26uMMTCK72dt4BKRCef3++nq6hpuPwSRCWGnHCJxm3DcxpDeIe0aRsAaY+jt6cZxjXx+q5mxiEy4yspKmpqa6OjomOyhiOAYg50y2I6DhYXrbCe4loXj8hD35Y14DApjEZlwXq+Xurq6yR6GzHKxZIrG9jD7O8J43C5yQ95Ju3WiMBYRkVklYTsc6oqwpzXd2KMglIHrrKfDY0thLCIis8Lx05R2tfaRTDnkBX3D2pw1Ec4YxpZlfR+4AWg3xpwz+Fg+8BBQCxwEPmCMOTZ+wxQRERkZxzG09kWpb+4jkkiRG/DhC0yNED5uOKN5ALjuTY/dBTxjjJkPPDP4sYiIyJRhjKEzHGfj3g5ePnAMt8tFcZYfn2dqBTEMY2ZsjNloWVbtmx6+Ebhi8Pf/D9gA3DmG4xIRERmxnoEE9S19tPfFyczwUJI9stOUJspI7xmXGGNaBn/fCpSc7omWZd0B3AFQXV09wrcTERE5s3DcZm9bH0e6o/g97ikfwseNegOXMcZYlnXayn1jzP3A/ZA+KGK07yciIvJmby5TKsrMmFYd3kYaxm2WZZUZY1osyyoD2sdyUCIiIsMxFcuURmKkYfwocDvwjcH/PjJmIxIRETkDO+VwtCfKzuZebMdMqTKlkRhOadNPSG/WKrQsqwm4l3QI/8yyrP8BHAI+MJ6DFBERgTfKlHY29zGQsMkNZEzJ3dFnazi7qT90mj+6aozHIiIickrpMqUEO4/20hNLkuP3UpwVmOxhjRl14BIRkSntLWVKWdNjh/TZUBiLiMiUFI7b7G7po+lYlIB3+pQpjYTCWEREppRoIkVjez8HOiN4PS6Ks6ZXmdJIKIxFRGRKSNgOBwfLlFzW9C1TGgmFsYiITKrjpynVt/SSciA36J3WZUojoTAWEZFJ4TiGlt50mVI0OXia0gwoUxoJhbGIiEwoYwwd4Tg7j/bRG02QE/CR5fdO9rAmlcJYREQmzLFIgvrmPjrCcbL8HkqyZ06t8GgojEVEZNz1xZLsae3n6LEoAd/MLlMaCYWxiIiMm4GETWN7mAMdEXyzpExpJBTGIiIy5uJ2igMdEfa2h3FZUJiVgUshfFoKYxERGTPJlMOR7gF2tfThGMgP+nDPklrhE3WG47zQ2Dns5yuMRURk1FKOobknys6WXhJJh9ygD697dpUpheM2L+3vYlNjJ7sHz1ceLoWxiIiMmDGGtr4YO5v76I/b5Aa85Ph9kz2sCZOwHV47fIxNjZ28dqSHlGMoz/Hz/vMrWTOvkPd9e3ivozAWEZER6QrHqW/uoysSJ8vvnZGnKZ2K4xjqW/rY1NjJywe6BxuWeLl2SQlr5hVSVxg6601qCmMRETkrvdEku1v6aOmNEvJ5Z0WtsDGGg10DPN/YyQv7Ojk2kCTgdXNBbR6Xzi9iaVn2qPpoK4xFRGRYInGbhrYwh7oiZHjdFGf5Z3yZUntfjOf3dfF8YydHe6K4XRYrqnJZM7eQ82vyxqx9p8JYRETeViyZYl9HmH3tYTxu14wvU+qLJXlpfxfPN3axpy29EWthSRafWFPHRXPyx6V1p8JYREROKWE7HBo80hBm9pGGcTvFlkPHeL6xk9eP9JIyhsq8AB+8oIo1cwsoGuf74QpjERE5yfEjDXe19pFMOeQFfTPySMOUY9jZ3Mumxk5eOdhNLOmQH/LxzmWlXDqvkOr84IQtwyuMRUQESO8Sbu1LH2k4kBg80jAws0LYGMP+zgibGjt5cV8XvdEkQZ+bi+cUsGZeIYtLR7cRa6QUxiIis9ypjjQszppZRxq29cXY1NjJ842dtPTG8LgsVlbnsmZeISurxm4j1kgpjEVEZrHuSIJdM/RIw95okj8OdsRqbA9jAYvLsrlheTkX1uWTmTF1InDqjERERCbMiUcaBmfQkYaxZIrNgxuxtjX14Biozg/y4QuruWRuAQWZGZM9xFNSGIuIzCIDiXSt8MGuCD73zDjSMOUYth/tYVNjF5sPdhO3HQozfdywvJxL5xVSlR+c7CGekcJYRGQWiCVT7O8I09gexu2yKMyc3rXCxhj2dRzfiNVJX8wmlOHm0nmFrJlXyMLSrGn19SmMRURmsITtcHiwVtjBTPta4ba+GM83drJpcCOW121xXnUel84r5Nyq3Gl7UpTCWERkBrJTDkd7otS3DNYKB3x4pmlQhWM2L+5Pt6Q83hFrSVk2f3JuOavr8gn6pn+UTf+vQEREhpxcK2yTG8iYlrXCyZTDa4d7eK6hY+howorcALdcUMWaeYUUTtGNWCOlMBYRmQFOXSs8vcqUHGPY29bPpoZO/ri/i0giRc7g0YSXzi+itmDiOmJNNIWxiMg01x1JUN/cS2c4MS1rhZt7omxq7GRTQycd4TgZHhcX1OZz6bxCzqnIwT2N73EPl8JYRGSaOvFc4aDPM61qhXujSV7cl96Ita8jgmXBsvIc3r+qkgtq8/F73ZM9xAmlMBYRmWbCcZvGwXOFfR7XtDlX+PjJSM81vNGQo6YgyG2ra7hkXgF5Qd9kD3HSKIxFRKaJWDJFY3uYfR1hvNPkXGHHMdS39LGpsZOXD3QTTabID02vhhwTQWEsIjLFxe0Uh7oG2NPajwUUToNa4cPdA2xq6OD5fV10RxIEvG5W1+Vz6fxCFpdlT/n/iZhoCmMRkSkqOXiucH1zL46B3KB3Sp8r3B1J8MK+9EasQ90DuC2Lc6tyuG11NefX5E/6yUhT2ajC2LKs/wV8EjDAduDjxpjYWAxMRGS2SjmG5mNRdrb2kkg65AZ9U7azVDSR4pWD3TzX2MnOo70YYG5RiI9dUsvFcwrIDsysoxjHy4jD2LKsCuBzwBJjTNSyrJ8BtwAPjNHYRERmFccxtPXF2NnSSySerrHN8U+9TU1DBzM0dPLKwWMkUg7FWRnctLKCS+cVUpY7vUqrpoLRLlN7gIBlWUkgCDSPfkgiIrOLMYbOcIKdR3vpiSbJ9nspzpp6M8pDXRE2NnTyfGMnvdEkoQw3a+cXsnZ+EQtKMqfFju6pasRhbIw5alnWN4HDQBR40hjz5JiNTERkFuiOJNjV3Ed7f3ywYcfUqhXuGUjwwr4uNu7tGLoPvLI6l7Xzi1hZPX0PZphqRrNMnQfcCNQBPcDPLcu6zRjzozc97w7gDoDq6upRDFVEZOZ4o2FHjKDPTWnO1AnhhO2w5dAxNjZ0DNUDD90HnltAtn/qzdqnu9EsU78DOGCM6QCwLOtXwCXASWFsjLkfuB9g1apVZhTvJyIy7YXjNg1t/RzuGhhs2JExJZZ3jTHsbQuzsaGDP+7vYiDxRj3w2vmFVOapHng8jSaMDwMXWZYVJL1MfRWweUxGJSIyw0QTKRrb+znQGcHjdlE0RUK4vS/GxoZONjV20NaX7gt9YW0+axcUsbQse8rXM88Uo7ln/JJlWb8AXgVs4DUGZ8AiIpIWS6Y42Blhb3sYtwUFU6Bhx0DC5qX93Wxs6GD3YCORJeXZ3LSyktV1s68v9FQwqt3Uxph7gXvHaCwiIjNGwnY43B1hT0s/BsgP+ib19KF0OVIvGxs62Hywm2TKUJbj54Orqrh0/sw7H3i6UQcuEZExZA92zdrV2kfSdsgL+vBM4o7jw90DbNzbwfONnfQMliNdsbCYy+YXMrdI5UhThcJYRGQMpBxDc0+U+pY+YskUuQEfvsDkhPBQOVJDB4e63ihHumx+EStUjjQlKYxFREbhxK5Z4XiK3IB3Ukp/jpcjPdfQweuD5UhzCkPcfnEtl8xTOdJUpzAWERkBYwwd4Tj1R/vojSbJ8nspmeCuWcYYGtvDPLu3gxdVjjTlOGb41bwKYxGRs9QVjlPf0kdXOEFmhofiCe6a1R1J8FxDBxv3dtDcG8PndnFhXT6XqRxpUhhjSKQc4kmHeCoFxuL4rXjjpOzhvIbCWERkmHoGEuxu7ae1N0rQN7GtK0/sivV6Uw/GwMKSLO5YXs7qOfkEffpxPt6MMSRThridIm47GGPAAow12Mo0g7ygj2CGh6DPTcDrxsQj/cN5bf3tiYicQV8syd7WfpqORQl43RRn+SdkF7Ixhv2dEZ7d28EL+zqJxNPL0DeeW8FlCwopy9HpSOPFTjnEbId4MoVjwLLAYAj6PBRkZpAX9JLp9xD0eQh43aMuW1MYi4icRjhu09gW5lB3BJ974lpX9gwk2NTYybN7O2g6FsXrtrigNp/LFxRxTnmOlqHHUMpJz3RjSYeU4ww97ve6yQ16ySsIkh3wEvC5CXrd41ampjAWEXmT460r93dG8LpdFGZm4BrnELZTDq8e7uHZvR1sPXIMx8D84kz+x6V1XDyngFCGflyPhmMMCdshlkyRTL0Ruh63i9yAl/JcP7lBH0Gvh4DPjc8zseVf+tsVERkUS6bY3xGhsSPdurJwAlpXHuyK8OyeDjY1dhKO2+QGvdywvJzL5hdRkadl6JFI2M7J93VJb6jK9nuoyAuQF/QRGryvm+FxTYnGJwpjEZn14naKQ10D7G3tx2DID2aMa+vKvmiS5/d18uye9BnBHpfF+TV5XLGwiGUVuZPaNnM6cRxDbHCJ2XEcjhcSBX1u8kM+8kM+Mv0eQoP3dafy8r7CWERmrWTK4VBXhD2t/aQcyAt58bjGZ3nSdhy2Hulh494OXj3cQ8oxzCkM8fFLarlkbiGZfv04Pp3ju5hjydRg6VD6cbfbRV7AS1muP72L2eshmOGelh3G9LcvIrPOif2jbdshdxz7Rx/pHmDD3vQydF80SXbAy3VLS7lsQRHV+WrK8WYnzXaNgzEWWIZQhofi7Azygz5Cg7Ndv3dqLDGPBYWxiMwaKcfQfCzKztZeEkmHnHHqHz2QsHlhXxcb9rSzryOC27I4ryaXyxcUc25VzrjNvqebhO0Qs1PE7bfOdityA+QGfQR8bkK+8dvFPFUojEVkxnMcQ0tv+hCHSDxFbtBLjt83pu9hjGF3az/r97Tz0v5uEimHqrwAH7mohkvnFZIdmL29oU/ayew4Q8Eb9LkpzMyg4IR7uzNptns2FMYiMmO9cYhDH+FYkpyAj5LssQ3FnoEEG/d2sH5PB619MQJeN2vnF3LFwmLmFoVmXbCknPS93VgyRcoYLAssLHICHqryg+SHfAR9boI+z4SXD01lCmMRmXGMMXT0x9nZnD7EIdvvpSR77MqEUo7htSPH2LCng9cOp2uCF5Vm8Z6VFayuy8fvdY/Ze01ldsohNtiPOd0vw+B2uygI+qguCJIT8BL0uQn5PFN6J/NUoDAWkRnDGENnOEF9Sx89A+lDHMayf3RLT5QNe9MHNPREk+QEvLxrWRlXLCymPHdm1wQnUw7R5Mn3d30eN/khL4WZIbL8XkIZ6X7Ms201YCwojEVk2jPG0B1JnHySUtbYhHAsmeKlA91s2NPO7tZ+XBasqMpj3cIiVlTnzsjNWMfv7ybsFCZ9EgJ+b/r+buHx+7sZnlmzAjARFMYiMq11RxLsbu2jvS9OaIxOUjLGsK8jwoY97bywr4toMkVptp9bLqjisgVF5AXHdvPXZDo+400k3wjeoM9NcfbxjVXpGW+GR8E7nhTGIjItnXicYWCMQrg/lmRTYyfr93RwpHsAn9vF6jn5rFtYzKLSrGm//GoPBm/shKVmv9dNUdbxGa+Cd7IojEVkWukdSLK3vZ+jY3ScoWMM9c19/GF3O68c7MYe7Iz1iTV1rJlXMG3PCU45Zuge7/HNVT6Pm8JMH4WZGWQHFLxTyfT8LhORWeeNM4UHyPC4R32c4fGSpD/saaetL04ow807FpdwxcIiagpCYzjy8ecYQzyZnvXajoNlgdvloiDkY05WiJyAl0zd453SFMYiMqX1x5I0tPdzpDs6eKbwyGfCjjHsONrLH3a3s/ngMVLGsKg0i5vPr+LC2vxpU/easAfv86ZSYNInEuUGvZTnhsgL+cgcPJFoui+rzyYKYxGZksJxm8a2MIe602cKF2WOfCbcM5Bgw94O1u9up70/TmaGh2vPKeXKRcVUTPGSJGdwuTlmp3Cc9I3eoM9NWU4GBZkZZPnTs16d9DS9KYxFZEqJxG0a28Mc6orgcbsozMzANYIQdoxhe1N6FrzlUHoWvKQsmw+squKCKTwLTtgOAwkbO+VgLHBbFgWhDGoLg+QGfVpunqEUxiIyJUTiNvs6whzoTM+EC0YYwt2RBBv2tLN+Tzud4QRZfg/vXFbKlQuLKZtis2DHpFtHRhMp0pPedFlRRV6AwswMsgb7Nat71cynMBaRSTWQsNnXng5hl8sa0UzYcQyvN/Xwh93tvDrYnvKc8mw+fGENq2rzpsz5tinHMJCwiSVTGMBlWeQFvVTnp2e9WX7NemcrhbGITIpoIsW+jjD7O8K4XOml2LOdAXaF40P3grsiCbIDXm5YXs66hcWU5oxdG8yRSqYcBhLHO1kxeO/bR0l2FtmDO5xn+tGAMjwKYxGZUNFEigOdYRrbw7hHEMKOMWxr6uHpXelZsDGwvCKHj1xUw/k1eZMabgnbIZKwSaYcLNINNcpyMijK8pPl95CZ4dEOZzklhbGITIihEO6I4IKzDuG+aJINezt4Zlcb7f1xsgNe3n1uehY8lodBnI24nSIST5FKd9UgmOGhOj9IUVb6fu90bRgiE0/fKSIyrmLJFAc6IzS0h3EB+UHfsMtwjDHsbQvz1K42Xtrfhe0YFpdlccsF6R3REz0LTtgO4bidbqwBhDI81BUGhzpa6X6vjJTCWETGxfEQbmwPA4b8YMawQziaSLGpsYOndrVzpHuAgNfNVYtLeMfiYirzguM78BMcLzNKpt6Y+dYWpGe+Cl8ZSwpjERlTbw7hvLMI4UNdEZ7e1camxk5iSYfagiCfXFvHmrmFExJ8tuMQiac3XEH6nm/V4LJzjsJXxpHCWETGRCyZ4uDgcvTZhHDCdnjpQBdP72pjb1sYr9vi4jkFXL2khLlFmeO64ckxhoFEilgyhWMMPreL0hw/pTl+sv1etZSUCaMwFpFROR7Cje1hzFmEcFtfjKd3tbFhTwfhuE1ptp/bVtdw+YIiMv3j96MplkwRSdg4jsFlWRRmZbCwJJPckI8s7XaWSaIwFpEROTGEAXKHsTHLcQyvHenhqfpWXm/qxWXBqpp83rGkhKXl2SPquHUmx5c54tydAAAgAElEQVSe43YKC8gOeFlQkpXedOVXna9MDaMKY8uycoH/BM4hfVT1J4wxL47FwERkajpxOdoCcoJePK63D7RwzGbD3naeqk+XJeUFvbzvvEquXFRMfsg3puMzJn2wwkDCxhgLj9uiPNdPaXaA3KDu+8rUNNqZ8T8DTxhjbrYsywdM3DZHEZlQJ90TNobckO+MIXyoK8Lvd7bxfGMniZTD4rIsPnxhNatq88f0lKE3z37zQj7qCnMoCKXrfdXbWaa6EYexZVk5wGXAxwCMMQkgMTbDEpGp4qTd0cMIYdtxeOXAMZ6sb2V3az8+t4tL5xdyzZISagpCYzquSMIm5Rg8bhcVuX7KcgLa9SzT0mhmxnVAB/ADy7LOBbYAnzfGRE58kmVZdwB3AFRXV4/i7URkIr05hPNCb78xq2cgwR92t/P0rjaODSQpzsrg1tXVXLGwmMyM0W9POb7zOZq0MQay/F4WlmSla379Xs1+ZVqzjDEj+0TLWgX8EVhjjHnJsqx/BvqMMfec7nNWrVplNm/ePLKRisiEOJs6YWMM+zrCPLGzjT/u7yLlGJZX5nDtklJWVOWOOiBTjiEct4eWn4uz/FTkBcgP+QiNQcCLjDfLsrYYY1ad6Xmj+W5uApqMMS8NfvwL4K5RvJ6ITKKzCeGE7fDH/V38fmcr+zsjBLxurl5cwtVLSigf5ZnByZRDf8wm5Ti4XRYVeQHKcwLkBn34PNr5LDPTiMPYGNNqWdYRy7IWGmP2AFcB9WM3NBGZCCce4GCdIYS7wnGe3tXGM7vb6Y/ZVOQG+PiaWtbOKyLgG/l92oTt0B9PknIMGR43dYVBSrL9wyqXEpkJRrvO8+fAg4M7qfcDHx/9kERkIrz5FKXTHeBw/LCG3+1o4ZWD3Rjg/Oo8rllayjnl2SNukhFLpgjHkzgGQj43C0qyKB5sO6nGGzLbjCqMjTFbgTOuhYvI1DHcELZTDi8d6Obx7S3s74wQynBz/bIyrllSQlHWyI4sjCVThGNJHCAzw8PismyKs/3qfCWznnZAiMwS0USK/Z1h9rWHcVnWaUO4P5bkmV3tPFnfyrGBJOU5fj6xppa184tGVDI0FMAm3f1qaUXO4Hm/3rH4skRmBIWxyAw3kLA50BlhX3sYt8uiIJRxyl3OTccGeGJHKxsbOkimDMsrcrjjslKWV+aedZvKuJ2ifzCAMzM8CmCRM1AYi8xQAwmb/R0R9necPoSNMbze1MvvdrSwrakXr9ti7fwirltaSlX+2TXUO3ETVijDzeKybEqy/WRqCVrkjBTGIjPMcEI4bqd4rqGTJ3a0crQnSm7QywdWVXHVomKyA8OfvdqOQ1/UxnYc/F4384oyKcsNkO1XAIucDYWxyAxxYgi7ThPC3ZEET9a38syudsJxm7rCEJ9dN4+L6vKHfXqRYwz9MZuYncLrtqgpCFGeGyA3oC5YIiOlMBaZ5gYSNvvawxzojJw2hA90RnhsWzN/3N+Ng+GC2nzeeU4pC0uyhjWDNYOtKAcSNgDleQGq84MUnKFFpogMj8JYZJo6Uwin7wf38NttLexs7iPgdXPt0hKuXVpKcfbwSpMStkNfNIEDFIR8LC7LozArgwyPDmIQGUsKY5Fp5kwhnEw5vLCvk99ua6HpWJT8kI9bV1dz5aJigr4z/5NPOYa+WJJkyiHkc7O4PJuynIB6QYuMI/3rEpkmzhTC4bjNM7vaeGJnKz0DSarzg/zZFXO5eE7BGe8HG2OIDC5De1wW1flBKvOC5AbVDUtkIiiMRaa4M23M6uiP8fiOVtbvbiduOyyryOEzl5exrCLnjEGasB36YgkcA4WZPpaW51GU5cc7zM1cIjI2FMYiU9SJHbNOVaK0vyPMb7e38NL+LiwsLplbwLuWl1FTEHrb13WMITy4G9rvdbOoNJvyXC1Di0wm/esTmWKGekefIoQdY9h6pIfHtrVQ35LelHX9sjKuW1pKQWbG275u3E7RF7UBQ3lugJqCEAUhn8qRRKYAhbHIFBFLptjf8cYBDieGcDLlsKmxk8e2tXC0Z/ibshxj6IsmSQxuxlpWmUNZjn9EPaZFZPwojEUmWSyZ4kBnhMb2MBYnn6IUTaT4w+52Ht/RQnckQc0wN2XF7RR9sSRgUZkboKYgSH7Ip81YIlOUwlhkksSSKQ4OhrDBkBd8o4FGXzTJ73e28vv6ViLxFEvKsrlj7RyWV55+U9bxzljxZIpQhptlFbmaBYtMEwpjkQkWt1Mc6hqgobUfA+QEvXhc6VluZzjOY9ta+MPudhIph1U1ebz73HLml2Sd9vUStjN4QlL6XnBtYYj8oO4Fi0wnCmORCZKwHQ53RdjT1o9jIPeEEG46NsBvXm/m+cYuANbMK+BPzi2nMu/UJycZY4jEUwwkbTK8LhaVZVGRGyTg0yxYZDpSGIuMs2TK4Uj3ALta+0g5hryAb+h+b0NbP4++3szmQ8fI8Li4emkJ71pWRuFpdkanHENvNIGdMhRlZbC8KofCTPWHFpnuFMYi48ROOTQdi7KrtY+k7ZAb9OF1u9I9o4/08OjrzdS39BHKcPPe8yq4dmkp2f5TH18YS6bojyVxWRa1hSGqC4Knfa6ITD8KY5ExlnIMzcei7GztJZ50yA348AVcOI7hj/u7eGTrUQ52DZAf8vGRi2q4clHxKTdZGWPoj9vEEilCfjcrqnIpzQng86g7lshMozAWGSOOY2jpjVLf0kckniI36CXH7yPlGJ5r6ODXW4/S3BOjLMfPHZfN4dJ5hadsO5lyDD3RBCnHUJLt5/yaPG3IEpnhFMYio+Q4hvb+GDub++iPJckJ+CjJ9mKnHP6wu51Hth6lvT9OVX6Qz105n9V1+acM1nSHrPRS9JyiEFX5QbK0FC0yKyiMRUbIGENnOMHO5l56B5Jk+j2UZAdI2A5P7mzl0deb6YokmFMY4iMX1XBeTR6uU9QIh2M2AwmbUIaWokVmK4WxyAh0RxLsau6jIxwnM8NDcbafWDLFY9ta+O22ZnqiSRaUZPLJtXM49xSNOhxj6BlIYjsORZkZnFudQ+GbDoIQkdlDYSxyFnoHkuxq7aO1N0rQ56Ek289AwubXrx3lse0thOM255Rn8+dXzmNxWfZbQthOOfQMJDBATUGI2oIQOUEtRYvMdgpjkWHojyVpaO/ncNcAfq+H4iw/kXiKn285wu93tBJJpFhRlctNKytYcIpuWbFkule0x22xsCybqjw16BCRNyiMRd5GNJGiob2fg50RPG4XxVl++mM2P33lCE/WtxJLOlxYm897VlZQV/jWc4SH7gf73ZxXlUdprv+UO6hFZHZTGIucwvGTlBraw7gsKMjMIBxPh/Dvd7aSsB0umlvATSsqqMo/uWXl8WML47ZDYaaPc6sLdD9YRN6WwljkBMlUun/07tZ+jEkfZxhJ2Dx0QghfMreAm86rpCI3cNLnHq8PdhxDRV6AuUWZ5AZ9k/SViMh0ojAW4YSuWS29JGyHvKCPgWSKn21+I4QvnlvAe1dWUpF3cggnUw49A0lcFswpClFTECKUoX9aIjJ8+okhs5oxhra+GDuae9NdswI+XFaKX7zaxO93thJPnj6Ej2/K8rpdLCnPojIvqLODRWREFMYya3WF49Q399EVSZDt9xLwwsOvNfHEYAhfNLeA966seMsxhgMJm/64TaZPm7JEZGwojGXW6Ysl2dXcR0tvlJDPS8jn4dHXj74RwnMKeO95J4fw8UMbookUeSEfF88poChTm7JEZGwojGXWGEjYNLSFOdgVwed2EfR5eGx7C7/f2UosmeKiOQXctPLk3dHGGHoHd0aXZPtZVZNHfsj3lmYeIiKjoTCWGS9hOxzoDLOnLV2mFPJ5+P3OVh7b3kI0kWL1nHzeu7LypBB2BndG2ylDZX6AeUVZ6pQlIuNGYSwzVsoxHOkeoL6lFzsFoQz34ClKzYTjNqtq8rj5/EpqCkInfU7PQALHGGoLQ9QVhnRykoiMO4WxzDhDO6SP9hJJpAj5PLxyoJOHtx6lZyDJ8socPrCqirlFmUOfY6cceqJJLGBucSa1BSG1qxSRCTPqMLYsyw1sBo4aY24Y/ZBERu5YJMGO5l66wglCPg87m/v41atNdIYTLCrN4nNXzmdxWfbQ8xO2Q180gctlsag0i6p8lSeJyMQbi5nx54FdQPaZnigyXsJxm90tfTQdi+L3umhsD/OLLU209sWYWxTiU2vnsKzijaMM43aK3miSDI+LcypzqcjVGcIiMnlGFcaWZVUC7wL+BvjfYzIikbMQt1Psaw/T2B7G7bI41BXhF1uaOHIsSlV+kC9cvYDza/KGQvh4ow6/N10jXJbrx6MaYRGZZKOdGf8T8JfAW8+MG2RZ1h3AHQDV1dWjfDuRtOObs3Y29+I4hqPHovx8SxP7OyOU5fj58yvncdGcAlyDITyQsOmPJcn0e1hVk0dpTgC3aoRFZIoYcRhblnUD0G6M2WJZ1hWne54x5n7gfoBVq1aZkb6fCKQ3Z7X3xdg+2L6ysz/Oz7c0sbu1n8JMH5++bA5r5xcNBe1AwiYcs8kKeFhdV0BJtl+NOkRkyhnNzHgN8G7Lsq4H/EC2ZVk/MsbcNjZDEzlZ70CSHc29tPfHCEdT/Pr1o2w5dIycgJePX1LLlYuKh5acI3GbcDxJbtDHxXMLKFS3LBGZwkYcxsaYu4G7AQZnxl9UEMt4iCZS7G3r52BnhEjc5omdrTzX0Inf6+aDq6q47pzSoR3Q4ZhNJGGTH/JxSVUhRZkZ6pYlIlOe6oxlykqmHA51RtjV2k8kbrNhbwdP17dhWXD9sjJuXFE+1JCjP5ZkIJGiINPHiupCCjPVslJEpo8xCWNjzAZgw1i8logxhpbeKNuP9tE3kODFA938bnsrMTvF5fOLuPn8SgoyM4A3QrgoM4Pz1TdaRKYpzYxlSjnetKOtN8ZrR3r4zbYW+qJJLqjN44OrqofOFD4ewsVZCmERmf4UxjIlDCRs9raG2d8ZZufRXn6zrYX2/jhLyrL50DULmFecrp7rjyUZSKYo1kxYRGYQhbFMquP3hetb+tjV0sdj21s53D1ATUGQu65bxPLKdNesE2fCq2rzyQt6FcIiMmMojGVSGGNo7Y2x/Wgve9v6eWx7C7ta+inJzjipYYeWo0VkNlAYy4TrHUiy42gvu1v7eLK+jZcOdJPl9/CxS2q5arBW+HiJkjZmichsoDCWCRNLpuuFdxztZf3uDtbvacey4D0ryvmTc8sJ+jyE4zYDAwkKQj5W1hRSoBAWkVlAYSzj7ngf6a1HjvHs3g6erG9jIJ7isgVFvH+wTCkSt2nri5IfymBFleqERWR2URjLuOroj7P1yDGea+jk8e0tdIYTLK/M4cMXVlNTEGIgYdPWFyM36OWSeeqYJSKzk8JYxkU4brOrpY8Ne9p5bFsLB7vSO6Q/tXYOyytziSZStPXHyPZ7uGRuAUVZCmERmb0UxjKmkimH/R1hNuzp4Lfbmtl+tI/8kI/PXD6XS+cXkrAd2vpjZGa4WV2br1OURERQGMsYOV6qtLGhg0e2NvPi/i78Hje3XFDFO88pA6AzHCfgdXOBzhMWETmJwlhGrTeaZPPBbn62+Qgb9nRgpwzXLinlpvMqCHjd9Awk8HlcrKzKpTw3MHTMoYiIpCmMZcTidoo9rf389OXDPLa9ld5oktV1+dxyQTWFWT56BpKkHIdzKnKoyg/iVQiLiJySwljOmuMYjh6L8svXjvDLLUc5cizK3KIQ//vqBcwtyuTYQJzegSQLS7OoKQiS4XFP9pBFRKY0hbGclWORBE/vauO/XzzEtqO9FIR8/M9181g9J5/eaJLuSJx5xZnMKcrE71UIi4gMh8JYhiWWTLHl0DH+87n9bGzoxOOyeP/5lVy/rJRo0qE7nKCuKMS84kyCPn1biYicDf3UlLeVcgwHOsP8x8b9PLa9lUjc5rIFRXxgVRUel0Vv1KYyL8DC0iyy/N7JHq6IyLSkMJbT6o4k+NGLh/jxy4dp7YuxuCyLj1xUS3FWBuG4TV4og9Vz8skN+iZ7qCIi05rCWN4ilkzx1M42vruhkd2t6WMN//fVCzinPJu+mI3X42KN+keLiIwZhbEMSTmG7U09/ONTe9nU2Inf6+a21TWsW1hEOGHjGMPqOnXNEhEZawpjAaC1N8r/faqBR14/SsJ2eMfiEt6zshzbMSQch5VVuVTkBdU1S0RkHCiMZ7lowuZHfzzMv2/cR2c4wYqqXD50QRWhDA+OA4tLs6kpCOHzqGGHiMh4URjPUinHsKmhg797Yjf1Lf2U5/i589qF1BSGsFMONQVB5pdkqVZYRGQCKIxnoUNdEf728V08Vd9GhsfNRy6q4eK5BSRsh+KsDJUpiYhMMIXxLDIQt/nXDft44IWDROI2Vyws5sZzy8GCLL+HpeU55IdUpiQiMtEUxrOA4xh+t6OFv3tiN4e7o8wvzuTW1TXkhbyEfG7OqcihJNuvMiURkUmiMJ7h9rT28fXf1PPCvi7ygl7+7Iq5LCrNwuO2WFKWPk1JO6RFRCaXwniG6osm+Pvf7+FnrzThGMN7VpRz2YIiPC6L+SVZzCkK6TQlEZEpQmE8w6RSDj995Qj/+PReusIJVtXkcdPKCkIZHqoLAiwsySaUob92EZGpRD+VZ5Ath47x1Ud2sLO5j4rcAF+6ZiEVeQEKMn2cU55DnjZniYhMSQrjGaCjP87/+e1OfrutBb/XzUcuqua86jwy/R6WaXOWiMiUpzCexuyUw388t59/Wd/IQCLFuoXFXLOkhKwMD0sqcqjKC+Bxq3OWiMhUpzCepp5v7OCeX+9kf2eEhSVZvH9VJUVZGcwpDKlzlojINKMwnmba+qJ89ZGd/H5nGzkBL3esrWNRaTbleQEWl2WTE1DnLBGR6UZhPE2kHMO/PdvId9fvI5ZMcd3SUtYtKqIgM4PllTkUZWbovrCIyDSlMJ4GNjV08JVf7+Bg1wCLSrN433kVlOcGWFquph0iIjPBiMPYsqwq4IdACWCA+40x/zxWAxNo641xzyM7eLK+jdygl0+tncPS8izmFmXqvrCIyAwympmxDXzBGPOqZVlZwBbLsp4yxtSP0dhmLTvl8G/P7ue76xtJpBzeubSUyxcWUl0QYml5ju4Li4jMMCMOY2NMC9Ay+Pt+y7J2ARWAwngUnmvo4CsP7+BQ9wBLyrJ5z8py6gpDqhcWEZnBxuSesWVZtcBK4KWxeL3ZqLU3yj2/3slTu9rIC3r55No6llfksKQ8m9qCkOqFRURmsFGHsWVZmcAvgb8wxvSd4s/vAO4AqK6uHu3bzTjJlMO/PbuP765vJJky6V3SC4uZXxpiUWk2QZ/22ImIzHSj+klvWZaXdBA/aIz51ameY4y5H7gfYNWqVWY07zfTbNzbwV/9ejtHuqMsLc/mxhXlLCjJYnllLvnqIy0iMmuMZje1BfwXsMsY849jN6SZL70kvYOndrVTEPLxqbV1rKjKYXlFHhV5AVwqVRIRmVVGMzNeA3wE2G5Z1tbBx75sjHl89MOameyUw39tOsA/Pd1AMuVw3dISrlpcwpLybOYVZ+p8YRGRWWo0u6k3AZrCDdPmg93c9avtNLaHBxt3VLKkLJtzKlWqJCIy22l30Dg7Fknwf35bz69eO0pOwMvtl9Rw0Zx8zq3MpSwnoFIlERFRGI8XxzH89JXDfON3uwnHbS5fUMS7lpWyrDKXuUWZ+DwqVRIRkTSF8Tiob+7lrl9tZ1tTL7UFQT592RzOr81naXk2WX4tSYuIyMkUxmMoHLf5+yd286M/HsLvdfOBVZWsW1TMyqpcdc8SEZHTUhiPAWMMv93WzNd/U09nOMGFdfm8Z0U559fkMacoE6+6Z4mIyNtQGI/Swc4Id/9qOy/u76I8x8/nrpzLJfOKtCQtIiLDpjAeoVgyxb/8oZF/37gPl2XxJ+eW8c5zyjivWkvSIiJydhTGI7BhTzt/9fAOjvZEWV6Rw/vOr2R1XT5zi7UkLSIiZ09hfBZae2Pc88gOnqpvozDTx6curePKxcWcU5GjJWkRERkxhfEw2CmHHzx/gH98qgHbcXjHkmJuXFHOqup8SnO0JC0iIqOjMD6DbU09fOkX29jT2s+Ckkzef34la+YVMq84S407RERkTCiMT6M/luTvntjNj186TCjDw4curOZd55SyrCpXvaRFRGRMKYzfxBjD73a08NVHdtIVTnBBXT4fWFXJ6roCKvPUS1pERMaewvgER7oH+KuHt7OxoZPSHD+fXTeXa5aWsrA0S8cbiojIuFEYA8mUw38+d4B/fnovDnDt0hI+sKqK86rzyAv5Jnt4IiIyw836MH718DHu/MU2GtrDLCzN4pZVVaxbVEx1fhCXS0vSIiIy/mZtGPdGk3zjd7v46ctHyAp4uO2iat57XgWLS3MI+LQkLSIiE2fWhfHxQx2+9mg93QMJLpqTz4dXV3Px3EIKMzMme3giIjILzaowPtw1wJcf3s6mxk4qcv38xVXzedeyMmoLQ3jUxlJERCbJrAjjZMrh/o37+ednGrCAG5aV8eHV1SyvyiUzY1ZcAhERmcJmfBJtPtjNnb/cxr6OCEvLsvnQ6iquWlSiNpYiIjJlzNgw7hlI8P89vpuHNh8hN+jl45fUcvOqSuarjaWIiEwxMy6MjTE8srWZr/9mJ73RJGvnFXLrRdVcPKeQnKDaWIqIyNQzo8L4YGeELz+8nRf2dVGVF+BTa+u4YXk5lXmqGRYRkalrRoTx0AatpxtwueA9K8q5dXUNyypz8HtVMywiIlPbtA/j1w4f485fbmNvW5hzKrK5dXU1Vy8pVc2wiIhMG9M2jMNxm394Yjc/fPEQ2QEvn1hTyy0XVjNHNcMiIjLNTMswfrq+ja/8egdtfTEumVvAbRfVsHZBkWqGRURkWppW6dXeF+Nrj+7k8R2tlGb7+V9XL+Dm8ysoy9E5wyIiMn1NizB2HMNDm4/wt4/tIppMcf05pdx+SS0rqnN1zrCIiEx7Uz6MG9vD3P2rbbxy8BjzijO5/eIablhernOGRURkxpiyYZywHb63YR//sr4Br9vFhy+s5uNraphTlIVbNcMiIjKDTMkw3nywm7t+tZ3G9jDnVefy0YtruWZpCUHflByuiIjIqEypdOuLJfm73+3mwZcOkx/y8dl1c/nIRbWU5vgne2giIiLjZsqE8RM7WvnqIzvoCMe5clExH1tTy8VzCvCqZlhERGa4SQ/j1t4YX31kB0/Wt1GZF+Ar1y/m5vOrdKiDiIjMGpMWxo5jePClQ3zjid0kbIf3nVfBJy6pZXF5jg51EBGRWWVUYWxZ1nXAPwNu4D+NMd8Yzuftbevnrl9u49XDPSwpy+bja2q4YXkFAZ9qhkVEZPYZcRhbluUGvgtcDTQBr1iW9agxpv50n2MMfOvJPXzv2X34PW7uWFvHbRfVUF0QGukwREREpr3RzIwvBBqNMfsBLMv6KXAjcNow3tvez3f+0MiawX7SVy8p0aEOIiIy640mjCuAIyd83ASsfvOTLMu6A7gDwF86l7uuW8jNq6p0xKGIiMigcZ+WGmPuN8asMsasmleSzacvn6sgFhEROcFowvgoUHXCx5WDj51Whsel05VERETeZDRh/Aow37KsOsuyfMAtwKNjMywREZHZY8T3jI0xtmVZ/xP4PenSpu8bY3aO2chERERmiVHVGRtjHgceH6OxiIiIzEqqKxIREZlkCmMREZFJpjAWERGZZApjERGRSaYwFhERmWQKYxERkUmmMBYREZlkljFm4t7MsvqBPRP2hrNXIdA52YOYBXSdJ4au88TQdR4fNcaYojM9aVRNP0ZgjzFm1QS/56xjWdZmXefxp+s8MXSdJ4au8+TSMrWIiMgkUxiLiIhMsokO4/sn+P1mK13niaHrPDF0nSeGrvMkmtANXCIiIvJWWqYWERGZZApjERGRSTYhYWxZ1nWWZe2xLKvRsqy7JuI9ZwvLsg5alrXdsqytlmVtHnws37KspyzLahj8b95kj3M6sizr+5ZltVuWteOEx055ba20bw9+j2+zLOu8yRv59HGaa/w1y7KODn5Pb7Us6/oT/uzuwWu8x7Ksaydn1NOPZVlVlmWttyyr3rKsnZZlfX7wcX0/TxHjHsaWZbmB7wLvBJYAH7Isa8l4v+8ss84Ys+KEGsG7gGeMMfOBZwY/lrP3AHDdmx473bV9JzB/8NcdwPcmaIzT3QO89RoD/N/B7+kVxpjHAQZ/btwCLB38nH8d/PkiZ2YDXzDGLAEuAj47eD31/TxFTMTM+EKg0Riz3xiTAH4K3DgB7zub3Qj8v8Hf/z/gPZM4lmnLGLMR6H7Tw6e7tjcCPzRpfwRyLcsqm5iRTl+nucancyPwU2NM3BhzAGgk/fNFzsAY02KMeXXw9/3ALqACfT9PGRMRxhXAkRM+bhp8TMaGAZ60LGuLZVl3DD5WYoxpGfx9K1AyOUObkU53bfV9Prb+5+Dy6PdPuM2iazwGLMuqBVYC/3979+8aRRRFcfx78FehdoE0WgSxD1aCElJZpLMRGw1iYRELaxtbG21TiHYqBGIwheifoGkEjbYKBkk6LayMJ8Wb4CBuIcnO21nPp5ndNwtzudzlwts7s69JPY+MDHD133nbZyjbSguSZtonXe5dy/1rQ5DcDs0icAqYBr4C9+qGMz4kHQOWgVu2v7fPpZ7r6qIZbwAnW+9PNGuxD2xvNMctYIWybbe5u6XUHLfqRTh2BuU2db5PbG/a3rb9C3jA763o5HgPJB2iNOLHtp81y6nnEdFFM14DTkuaknSYMoCx2sF1x56ko5KO774GLgDvKfmdbz42DzyvE+FYGpTbVeBqM4V6FvjW2v6Lf/DHb5MXKTUNJceXJR2RNEUZLnrTdXx9JEnAQ+Cj7futU6nnETH0f22y/VPSTeAVcAB4ZHt92Nf9T0wCK+V7xktFmNAAAACRSURBVEHgie2XktaAJUnXgc/ApYox9pakp8AsMCHpC3AHuMvfc/sCmKMMFf0ArnUecA8NyPGspGnKlukn4AaA7XVJS8AHynTwgu3tGnH30DngCvBO0ttm7Tap55GRx2FGRERUlgGuiIiIytKMIyIiKkszjoiIqCzNOCIiorI044iIiMrSjCMiIipLM46IiKhsB4eda9jDHrdwAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "model.plot_cumulative_hazard(figsize=(8,5))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 2 }