{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "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": "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. \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": 80,
   "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": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 100)"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0IAAAE/CAYAAABrfXNCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dfbRddX3n8ffHGwOCgii0I0kkUSMatINtDHa01qUgsbTEWaMSWitYO1l2ZGq1WmNrUdPaou1o6yqtMhW1WgyIrb1TYxkUH9oqmKAoJpgxRCSJCAqIDygQ+M4fe996ON6be3LvuQ/Jfr/WOou9fw9n//a5253zcf/2PqkqJEmSJKlLHjDXA5AkSZKk2WYQkiRJktQ5BiFJkiRJnWMQkiRJktQ5BiFJkiRJnWMQkiRJktQ5BiFJ2ockNyQ5eYK6X0iyfbbHNJ+k8e4ktyf53H72XZqkkiyYqfHNhHbMjxmw7RuSvH+mx9Szvc4fk5I0KIOQJE1RVf1rVR0/WbvZ/jI8y54GnAIsrqpVcz2YrukPZYMek5Ikg5AkHdDmwdWU44AbquoHczyOcbVXrA7If+vmwd9Wkg5qB+Q/DpI0y05M8qUkdyS5OMmhAEmekWT3WKMkr0myJ8n3kmxP8qwkq4HfB85I8v0kX2zbHptkNMltSXYk+e897/OgJO9tp5tdl+T3+rZzQ7utLwE/SLIgyfok17fb3pbkv/a0PzvJvyd5W5LvJNmZ5L+05buS3JLkrIl2fqKxJnkJ8LfAz7f79sZx+j4gyeuSfL3dzt8lObKv2W8k+UaSm5K8qqfvqiRbknw3yc1J3tpT95Qkn2n354tJntFT98kkb0ry78CdwKuTbOkb1yuSjLbLhyT58yQ3ttt5R5IH9bR9dTu2byT5jYk+p7btsiSfav8OlwNH99Td73hpy/5j6mV75fDSJO9P8l3g7PYz+Gy7nzcl+askC9v2n27f5ovt53/GOMfk49vP4ztJtiY5vafuPUnOT/KRdrxXJXl0W5f2eLml/fyvTfKEfe27JB1wqsqXL1++fE3wAm4APgccCzwMuA54aVv3DGB3u3w8sAs4tl1fCjy6XX4D8P6+9/008NfAocCJwLeAZ7Z15wGfAo4CFgNfGttOz5iuAZYAD2rLnt+O8QHAGcAPgEe0dWcDe4EXAyPAHwM3AucDhwDPBr4HPHiCz2BfYz0b+Ld9fH6/AewAHgU8GPgH4H09n1EBHwAOB57YvvfJbf1ngV9vlx8MPKVdXgTcCvxSu7+ntOvHtPWfbPfvBGABcGS7f8t7xrUZWNsuvw0Ybf++DwH+D/Cnbd1q4GbgCe0YL2rH/JgJ9vezwFvbz/Xp7Xbf33+89P0tx/b3DcA9wHPb/XoQ8HPAU9r9WEpz/P1OT//7jYX7H5MPbD/73wcWAs9sx3N8W/+e9nNb1b7/3wMb27pTgauBhwIBHk97PPny5cvXwfLyipAkTe7tVfWNqrqN5kvyieO0uZfmy++KJA+sqhuq6vrx3izJEuCpwGuq6kdVdQ3NlZUXtU1eAPxJVd1eVbuBt08wpl1V9UOAqvpgO8b7qupi4Ks0X3DHfK2q3l1V9wIX04SoDVV1V1X9X+Bu4CceADDAWCfza8Bbq2pnVX0feC2wNvef9vXGqvpBVV0LvBs4sy2/B3hMkqOr6vtVdWVb/kJgU1Vtavf3cmALTTAa856q2lpVe6vqDuCfxt43yXLgccBokgDrgFdU1W1V9T3gT4C17fu8AHh3VX25mul/b5hoR5M8Engy8Ift5/ppmuNlf3y2qj7c7tcPq+rqqrqy3Y8bgHcCvzjgez2FJkCeV1V3V9UVwD/z488X4B+r6nNVtZcmCI0d2/fQhMLHAamq66rqpv3cF0ma1wxCkjS5b/Ys30nz5fJ+qmoH8Ds0X5RvSbIxybETvN+xwNiX7jFfp7nSMVa/q6eud3ncsiQvSnJNOwXqOzRXMI7uaXJzz/JYeOov+4n9GmCskzm2bd/bdwHw0xPsy9fbPgAvAR4LfCXJ5iS/3JYfBzx/bF/b/X0a8IgJ3hOaKzljAeBXgQ9X1Z3AMcBhwNU97/UvbfnY+PvHt699vb3uf7/UvtqPp//v+tgk/5zkm+10uT/h/n/XfTkW2FVV9/WNp/dvN+6x3Yamv6K5anhLkguSHLF/uyJJ85tBSJKGpKouqqqn0XxRL+DNY1V9Tb8BPCzJQ3rKHgnsaZdvopkSN2bJeJsbW0hyHPC/gXOAh1fVQ4Ev00xpmq7JxjpI/+P6+u7l/sFsSV/9NwCq6qtVdSbwUzSf5aVJDqcJC++rqof2vA6vqvN63qf/M78cOCbJiTSB6KK2/Ns0IfCEnvc6sqrGQuFN44xvIjcBR7VjHK/9D2hCFwBJRvhx4Jpo3H8DfIVmWt8RNNPcBv27fgNYkvs/LGLgv11Vvb2qfg5YQRNIXz3gdiXpgGAQkqQhSHJ8kmcmOQT4Ec2X67H/J/5mYOnYF9Kq2gV8BvjTJIcm+Rmaqx9jj9i+BHhtkqOSLKIJOPtyOM0X6G+1Y3kxzRWhaRtgrJP5APCK9iECD6a5onFxOxVrzB8mOSzJCTT3MV3c7scLkxzTXtH4Ttv2vnbbv5Lk1CQj7biekaQ3PPbvxz3AB4E/o7kX6PK2/D6aEPm2JD/VbndRklPbrpfQPLRgRZLDgNfvYxtfp5mi98YkC5M8DfiVnib/Dzg0yWlJHgi8jmY65b48BPgu8P0kjwN+q6/+Zpr7r8ZzFc1Vnt9L8sA0D5T4FWDjJNskyZOTnNSO8wc0x/R9k3STpAOKQUiShuMQmoccfJtmutFP0dwPA80XcIBbk3y+XT6T5ub3bwD/CLy+qj7W1m0AdgNfAz4GXArcNdGGq2ob8L9obtS/meahA/8+jJ0aYKyTuRB4H80DF75G84X6f/a1+RTNTf0fB/68vWcJmgcVbE3yfeAvaR5u8MM2nK2huTryLZorRK9m8n/TLgJOBj7YF8Re027/ynb62cdoHn5BVX0U+AvgirbNFZNs41eBk4DbaELT341VtPcq/Q+ae6z20ASM3eO8R69Xte/5PZrAdnFf/RuA97bT+l7QW1FVd9MEn+fQHJd/Dbyoqr4yyTYBjmi3dzvNdLpbaUKkJB00UtV/FV6SNJ8k+S2aEDDoTfKSJGkSXhGSpHkmySOSPDXNb/AcD/wuzZUYSZI0JP5qtSTNPwtpHpO8jObemI0005okSdKQODVOkiRJUuc4NU6SJElS5xiEJEmSJHXOvLtH6Oijj66lS5fO9TAkSZIkzWNXX331t6uq/4epBzbvgtDSpUvZsmXLXA9DkiRJ0jyW5OvT6e/UOEmSJEmdYxCSJEmS1DkGIUmSJEmdYxCSJEmS1DkGIUmSJEmdYxCSJEmS1DkGIUmSJEmdYxCSJEmS1DkDBaEkq5NsT7Ijyfp9tPtvSSrJyp6y17b9tic5dRiDliRJkqTpWDBZgyQjwPnAKcBuYHOS0ara1tfuIcDLgat6ylYAa4ETgGOBjyV5bFXdO7xdkCRJkqT9M8gVoVXAjqraWVV3AxuBNeO0+yPgzcCPesrWABur6q6q+hqwo30/SZIkSZozk14RAhYBu3rWdwMn9TZI8rPAkqr6SJJX9/W9sq/von1t7No9d7B0/UcmrL/hvNMGGLIkSZIkTWyQILRPSR4AvBU4exrvsQ5YBzByxDHTHZIkSZIk7dMgQWgPsKRnfXFbNuYhwBOATyYB+E/AaJLTB+gLQFVdAFwAcMgjltd+jF+SJEmS9tsg9whtBpYnWZZkIc3DD0bHKqvqjqo6uqqWVtVSmqlwp1fVlrbd2iSHJFkGLAc+N/S9kCRJkqT9MOkVoaram+Qc4DJgBLiwqrYm2QBsqarRffTdmuQSYBuwF3iZT4yTJEmSNNcGukeoqjYBm/rKzp2g7TP61t8EvGmK45MkSZKkoRvoB1UlSZIk6WBiEJIkSZLUOQYhSZIkSZ1jEJIkSZLUOQYhSZIkSZ1jEJIkSZLUOQYhSZIkSZ1jEJIkSZLUOQYhSZIkSZ1jEJIkSZLUOQYhSZIkSZ1jEJIkSZLUOQYhSZIkSZ1jEJIkSZLUOQYhSZIkSZ1jEJIkSZLUOQYhSZIkSZ1jEJIkSZLUOQYhSZIkSZ1jEJIkSZLUOQYhSZIkSZ1jEJIkSZLUOQYhSZIkSZ1jEJIkSZLUOQMFoSSrk2xPsiPJ+nHqX5rk2iTXJPm3JCva8qVJftiWX5PkHcPeAUmSJEnaXwsma5BkBDgfOAXYDWxOMlpV23qaXVRV72jbnw68FVjd1l1fVScOd9iSJEmSNHWDXBFaBeyoqp1VdTewEVjT26CqvtuzejhQwxuiJEmSJA3XIEFoEbCrZ313W3Y/SV6W5HrgLcBv91QtS/KFJJ9K8gvjbSDJuiRbkmy598479mP4kiRJkrT/hvawhKo6v6oeDbwGeF1bfBPwyKp6EvBK4KIkR4zT94KqWllVK0cOO3JYQ5IkSZKkcQ0ShPYAS3rWF7dlE9kIPBegqu6qqlvb5auB64HHTm2okiRJkjQcgwShzcDyJMuSLATWAqO9DZIs71k9DfhqW35M+7AFkjwKWA7sHMbAJUmSJGmqJn1qXFXtTXIOcBkwAlxYVVuTbAC2VNUocE6Sk4F7gNuBs9ruTwc2JLkHuA94aVXdNhM7IkmSJEmDmjQIAVTVJmBTX9m5Pcsvn6Dfh4APTWeAkiRJkjRsQ3tYgiRJkiQdKAxCkiRJkjrHICRJkiSpcwxCkiRJkjrHICRJkiSpcwxCkiRJkjrHICRJkiSpcwxCkiRJkjrHICRJkiSpcwxCkiRJkjrHICRJkiSpcwxCkiRJkjrHICRJkiSpcwxCkiRJkjrHICRJkiSpcwxCkiRJkjrHICRJkiSpcwxCkiRJkjrHICRJkiSpcwxCkiRJkjrHICRJkiSpcwxCkiRJkjrHICRJkiSpcwYKQklWJ9meZEeS9ePUvzTJtUmuSfJvSVb01L227bc9yanDHLwkSZIkTcWkQSjJCHA+8BxgBXBmb9BpXVRVT6yqE4G3AG9t+64A1gInAKuBv27fT5IkSZLmzCBXhFYBO6pqZ1XdDWwE1vQ2qKrv9qweDlS7vAbYWFV3VdXXgB3t+0mSJEnSnFkwQJtFwK6e9d3ASf2NkrwMeCWwEHhmT98r+/oumtJIJUmSJGlIhvawhKo6v6oeDbwGeN3+9E2yLsmWJFvuvfOOYQ1JkiRJksY1SBDaAyzpWV/clk1kI/Dc/elbVRdU1cqqWjly2JEDDEmSJEmSpm6QILQZWJ5kWZKFNA8/GO1tkGR5z+ppwFfb5VFgbZJDkiwDlgOfm/6wJUmSJGnqJr1HqKr2JjkHuAwYAS6sqq1JNgBbqmoUOCfJycA9wO3AWW3frUkuAbYBe4GXVdW9M7QvkiRJkjSQVNXkrWbRIY9YXo846y8mrL/hvNNmcTSSJEmS5qMkV1fVyqn2H9rDEiRJkiTpQGEQkiRJktQ5BiFJkiRJnWMQkiRJktQ5BiFJkiRJnWMQkiRJktQ5BiFJkiRJnWMQkiRJktQ5BiFJkiRJnWMQkiRJktQ5BiFJkiRJnWMQkiRJktQ5BiFJkiRJnWMQkiRJktQ5BiFJkiRJnWMQkiRJktQ5BiFJkiRJnWMQkiRJktQ5BiFJkiRJnWMQkiRJktQ5BiFJkiRJnWMQkiRJktQ5BiFJkiRJnTNQEEqyOsn2JDuSrB+n/pVJtiX5UpKPJzmup+7eJNe0r9FhDl6SJEmSpmLBZA2SjADnA6cAu4HNSUaraltPsy8AK6vqziS/BbwFOKOt+2FVnTjkcUuSJEnSlA1yRWgVsKOqdlbV3cBGYE1vg6r6RFXd2a5eCSwe7jAlSZIkaXgGCUKLgF0967vbsom8BPhoz/qhSbYkuTLJc6cwRkmSJEkaqkmnxu2PJC8EVgK/2FN8XFXtSfIo4Iok11bV9X391gHrAEaOOGaYQ5IkSZKknzDIFaE9wJKe9cVt2f0kORn4A+D0qrprrLyq9rT/3Ql8EnhSf9+quqCqVlbVypHDjtyvHZAkSZKk/TVIENoMLE+yLMlCYC1wv6e/JXkS8E6aEHRLT/lRSQ5pl48Gngr0PmRBkiRJkmbdpFPjqmpvknOAy4AR4MKq2ppkA7ClqkaBPwMeDHwwCcCNVXU68HjgnUnuowld5/U9bU6SJEmSZt1A9whV1SZgU1/ZuT3LJ0/Q7zPAE6czQEmSJEkatoF+UFWSJEmSDiYGIUmSJEmdYxCSJEmS1DkGIUmSJEmdYxCSJEmS1DkGIUmSJEmdYxCSJEmS1DkGIUmSJEmdYxCSJEmS1DkGIUmSJEmdYxCSJEmS1DkGIUmSJEmdYxCSJEmS1DkGIUmSJEmdYxCSJEmS1DkGIUmSJEmdYxCSJEmS1DkGIUmSJEmdYxCSJEmS1DkGIUmSJEmdYxCSJEmS1DkGIUmSJEmdYxCSJEmS1DkDBaEkq5NsT7Ijyfpx6l+ZZFuSLyX5eJLjeurOSvLV9nXWMAcvSZIkSVMxaRBKMgKcDzwHWAGcmWRFX7MvACur6meAS4G3tH0fBrweOAlYBbw+yVHDG74kSZIk7b9BrgitAnZU1c6quhvYCKzpbVBVn6iqO9vVK4HF7fKpwOVVdVtV3Q5cDqweztAlSZIkaWoGCUKLgF0967vbsom8BPjoFPtKkiRJ0oxbMMw3S/JCYCXwi/vZbx2wDmDkiGOGOSRJkiRJ+gmDXBHaAyzpWV/clt1PkpOBPwBOr6q79qdvVV1QVSurauXIYUcOOnZJkiRJmpJBgtBmYHmSZUkWAmuB0d4GSZ4EvJMmBN3SU3UZ8OwkR7UPSXh2WyZJkiRJc2bSqXFVtTfJOTQBZgS4sKq2JtkAbKmqUeDPgAcDH0wCcGNVnV5VtyX5I5owBbChqm6bkT2RJEmSpAENdI9QVW0CNvWVnduzfPI++l4IXDjVAUqSJEnSsA30g6qSJEmSdDAxCEmSJEnqHIOQJEmSpM4xCEmSJEnqHIOQJEmSpM4xCEmSJEnqHIOQJEmSpM4xCEmSJEnqHIOQJEmSpM4xCEmSJEnqHIOQJEmSpM4xCEmSJEnqHIOQJEmSpM4xCEmSJEnqHIOQJEmSpM4xCEmSJEnqHIOQJEmSpM4xCEmSJEnqHIOQJEmSpM4xCEmSJEnqHIOQJEmSpM4xCEmSJEnqHIOQJEmSpM4ZKAglWZ1ke5IdSdaPU//0JJ9PsjfJ8/rq7k1yTfsaHdbAJUmSJGmqFkzWIMkIcD5wCrAb2JxktKq29TS7ETgbeNU4b/HDqjpxCGOVJEmSpKGYNAgBq4AdVbUTIMlGYA3wH0Goqm5o6+6bgTFKkiRJ0lANMjVuEbCrZ313WzaoQ5NsSXJlkufu1+gkSZIkaQYMckVouo6rqj1JHgVckeTaqrq+t0GSdcA6gJEjjpmFIUmSJEnqskGuCO0BlvSsL27LBlJVe9r/7gQ+CTxpnDYXVNXKqlo5ctiRg761JEmSJE3JIEFoM7A8ybIkC4G1wEBPf0tyVJJD2uWjgafSc2+RJEmSJM2FSYNQVe0FzgEuA64DLqmqrUk2JDkdIMmTk+wGng+8M8nWtvvjgS1Jvgh8Ajiv72lzkiRJkjTrBrpHqKo2AZv6ys7tWd5MM2Wuv99ngCdOc4ySJEmSNFQD/aCqJEmSJB1MDEKSJEmSOscgJEmSJKlzDEKSJEmSOscgJEmSJKlzDEKSJEmSOscgJEmSJKlzDEKSJEmSOscgJEmSJKlzDEKSJEmSOscgJEmSJKlzDEKSJEmSOscgJEmSJKlzDEKSJEmSOscgJEmSJKlzDEKSJEmSOscgJEmSJKlzDEKSJEmSOscgJEmSJKlzDEKSJEmSOmfBXA9gfy1d/5F91t9w3mmzNBJJkiRJByqvCEmSJEnqHIOQJEmSpM4xCEmSJEnqnIGCUJLVSbYn2ZFk/Tj1T0/y+SR7kzyvr+6sJF9tX2cNa+CSJEmSNFWTBqEkI8D5wHOAFcCZSVb0NbsROBu4qK/vw4DXAycBq4DXJzlq+sOWJEmSpKkb5IrQKmBHVe2sqruBjcCa3gZVdUNVfQm4r6/vqcDlVXVbVd0OXA6sHsK4JUmSJGnKBglCi4BdPeu727JBDNQ3ybokW5JsuffOOwZ8a0mSJEmamnnxsISquqCqVlbVypHDjpzr4UiSJEk6yA0ShPYAS3rWF7dlg5hOX0mSJEmaEYMEoc3A8iTLkiwE1gKjA77/ZcCzkxzVPiTh2W2ZJEmSJM2ZSYNQVe0FzqEJMNcBl1TV1iQbkpwOkOTJSXYDzwfemWRr2/c24I9owtRmYENbJkmSJElzZsEgjapqE7Cpr+zcnuXNNNPexut7IXDhNMYoSZIkSUM1Lx6WIEmSJEmzySAkSZIkqXMMQpIkSZI6xyAkSZIkqXMMQpIkSZI6xyAkSZIkqXMMQpIkSZI6xyAkSZIkqXMMQpIkSZI6xyAkSZIkqXMMQpIkSZI6xyAkSZIkqXMMQpIkSZI6xyAkSZIkqXMMQpIkSZI6xyAkSZIkqXMMQpIkSZI6xyAkSZIkqXMMQpIkSZI6xyAkSZIkqXMMQpIkSZI6xyAkSZIkqXMMQpIkSZI6Z6AglGR1ku1JdiRZP079IUkubuuvSrK0LV+a5IdJrmlf7xju8CVJkiRp/y2YrEGSEeB84BRgN7A5yWhVbetp9hLg9qp6TJK1wJuBM9q666vqxCGPe0JL139kn/U3nHfaLI1EkiRJ0nw1yBWhVcCOqtpZVXcDG4E1fW3WAO9tly8FnpUkwxumJEmSJA3PIEFoEbCrZ313WzZum6raC9wBPLytW5bkC0k+leQXpjleSZIkSZq2SafGTdNNwCOr6tYkPwd8OMkJVfXd3kZJ1gHrAEaOOGaGhyRJkiSp6wa5IrQHWNKzvrgtG7dNkgXAkcCtVXVXVd0KUFVXA9cDj+3fQFVdUFUrq2rlyGFH7v9eSJIkSdJ+GCQIbQaWJ1mWZCGwFhjtazMKnNUuPw+4oqoqyTHtwxZI8ihgObBzOEOXJEmSpKmZdGpcVe1Ncg5wGTACXFhVW5NsALZU1SjwLuB9SXYAt9GEJYCnAxuS3APcB7y0qm6biR2RJEmSpEENdI9QVW0CNvWVnduz/CPg+eP0+xDwoWmOUZIkSZKGaqAfVJUkSZKkg4lBSJIkSVLnGIQkSZIkdY5BSJIkSVLnGIQkSZIkdY5BSJIkSVLnGIQkSZIkdY5BSJIkSVLnGIQkSZIkdY5BSJIkSVLnLJjrAcy2pes/ss/6G847bcr9J+srSZIkaX7wipAkSZKkzuncFaH5bLpXqw7UbUuSJEmzzStCkiRJkjrHICRJkiSpc5wa12eyKWLT6ev0Mh1IPJ4lSdLBzCtCkiRJkjrHICRJkiSpc5wadwCZzrS9mZ7G5DSqmeHnqmHyeJIk6ce8IiRJkiSpcwxCkiRJkjrHqXGaFfuakuN0HEmShsMpsNLgvCIkSZIkqXMGCkJJVifZnmRHkvXj1B+S5OK2/qokS3vqXtuWb09y6vCGLkmSJElTM+nUuCQjwPnAKcBuYHOS0ara1tPsJcDtVfWYJGuBNwNnJFkBrAVOAI4FPpbksVV177B35EAwnae+Hcjbnq6ZHvtk0wRmcvsH8xSFmZwOOZNTP+Z6Wsl0tj/T/1uZz0+unMzBOvb5fDzONadkj+9A/ptO5mD9m0/3bzbT/adjLrcNg10RWgXsqKqdVXU3sBFY09dmDfDedvlS4FlJ0pZvrKq7quprwI72/SRJkiRpzgwShBYBu3rWd7dl47apqr3AHcDDB+wrSZIkSbMqVbXvBsnzgNVV9Zvt+q8DJ1XVOT1tvty22d2uXw+cBLwBuLKq3t+Wvwv4aFVd2reNdcC6dvUJwJenv2vSwI4Gvj3Xg1BneLxpNnm8aTZ5vGm2HV9VD5lq50Een70HWNKzvrgtG6/N7iQLgCOBWwfsS1VdAFwAkGRLVa0cdAek6fKY02zyeNNs8njTbPJ402xLsmU6/QeZGrcZWJ5kWZKFNA8/GO1rMwqc1S4/D7iimktNo8Da9qlyy4DlwOemM2BJkiRJmq5JrwhV1d4k5wCXASPAhVW1NckGYEtVjQLvAt6XZAdwG01Yom13CbAN2Au8rKtPjJMkSZI0fwwyNY6q2gRs6is7t2f5R8DzJ+j7JuBN+zGmC/ajrTQMHnOaTR5vmk0eb5pNHm+abdM65iZ9WIIkSZIkHWwGuUdIkiRJkg4q8yoIJVmdZHuSHUnWz/V4dHBJsiTJJ5JsS7I1ycvb8ocluTzJV9v/HjXXY9XBI8lIki8k+ed2fVmSq9rz3MXtQ2ikoUjy0CSXJvlKkuuS/LznOM2UJK9o/z39cpIPJDnUc5yGJcmFSW5pf6ZnrGzc81kab2+Puy8l+dlBtjFvglCSEeB84DnACuDMJCvmdlQ6yOwFfreqVgBPAV7WHmPrgY9X1XLg4+26NCwvB67rWX8z8LaqegxwO/CSORmVDlZ/CfxLVT0O+M80x57nOA1dkkXAbwMrq+oJNA/UWovnOA3Pe4DVfWUTnc+eQ/N06uU0v036N4NsYN4EIWAVsKOqdlbV3cBGYM0cj0kHkaq6qao+3y5/j+YLwiKa4+y9bbP3As+dmxHqYJNkMXAa8LfteoBnAmM/Ku3xpqFJciTwdJonuVJVd1fVd/Acp5mzAHhQ+xuShwE34TlOQ1JVn6Z5GnWvic5na4C/q8aVwEOTPGKybcynILQI2NWzvrstk4YuyVLgScBVwE9X1U1t1TeBn56jYeng8xfA7wH3tesPB75TVXvbdc9zGqZlwLeAd7fTMf82yeF4jtMMqKo9wJ8DN9IEoDuAq/Ecp5k10flsSjliPgUhaVYkeTDwIeB3quq7vXXtDwH7KEVNW5JfBm6pqqvneizqjAXAzwJ/U1VPAn5A3zQ4z3EalvbejDU0AfxY4HB+chqTNGOGcT6bT0FoD7CkZ31xWyYNTZIH0oSgv6+qf2iLbx67fNr+95a5Gp8OKk8FTk9yA81U32fS3L/x0HYaCXie03DtBnZX1VXt+qU0wchznGbCycDXqupbVXUP8A805z3PcZpJE53PppQj5lMQ2gwsb582spDmhrvROZo99BgAAAFISURBVB6TDiLt/RnvAq6rqrf2VI0CZ7XLZwH/NNtj08Gnql5bVYurainN+eyKqvo14BPA89pmHm8amqr6JrAryfFt0bOAbXiO08y4EXhKksPaf1/HjjfPcZpJE53PRoEXtU+PewpwR88UugnNqx9UTfJLNHPqR4ALq+pNczwkHUSSPA34V+BafnzPxu/T3Cd0CfBI4OvAC6qq/+Y8acqSPAN4VVX9cpJH0VwhehjwBeCFVXXXXI5PB48kJ9I8nGMhsBN4Mc3/6ek5TkOX5I3AGTRPZf0C8Js092V4jtO0JfkA8AzgaOBm4PXAhxnnfNaG8b+imZ55J/Diqtoy6TbmUxCSJEmSpNkwn6bGSZIkSdKsMAhJkiRJ6hyDkCRJkqTOMQhJkiRJ6hyDkCRJkqTOMQhJkiRJ6hyDkCRJkqTOMQhJkiRJ6pz/D6Jvg31fs63mAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1008x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(14,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": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20.75914285714286\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 wish 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, thefore  \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": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<lifelines.MixtureExponential: fitted with 10500 total observations, 10500 interval-censored observations>\n",
      "number of observations = 10500\n",
      "number of events observed = 0\n",
      "    log-likelihood = -42593.81\n",
      "        hypothesis = lambda1 != 10.3796, lambda2 != 41.5183, p != 0.5\n",
      "\n",
      "---\n",
      "         coef  se(coef)  lower 0.95  upper 0.95      p  -log2(p)\n",
      "lambda1 20.76      0.63       19.53       21.98 <0.005    203.27\n",
      "lambda2 38.17     17.50        3.86       72.47   0.85      0.24\n",
      "p        0.97      0.06        0.86        1.09 <0.005     49.37\n"
     ]
    }
   ],
   "source": [
    "\n",
    "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()"
   ]
  }
 ],
 "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
}