Case Study 1- Bayesian_Changepoint_Detection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Probabilistic Programming Primer\n",
    "## Bayesian Changepoint detection\n",
    "A very useful kind of approach when modelling is to have the ability to detect changepoints.\n",
    "\n",
    "Some examples \n",
    "- When the server began to get overloaded\n",
    "- When a marketing campaign began to kick in. \n",
    "- Lots of examples in process control/ manufacturing. \n",
    "- When behaviour of users say conversion rate began to change"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Coal mining disasters ## \n",
    "Consider the following time series of recorded coal mining disasters in the UK from 1851 to 1962 (Jarrett, 1979). The number of disasters is thought to have been affected by changes in safety regulations during this period. Unfortunately, we also have pair of years with missing data, identified as missing by a NumPy MaskedArray using -999 as the marker value.\n",
    "\n",
    "Next we will build a model for this series and attempt to estimate when the change occurred. At the same time, we will see how to handle missing data, use multiple samplers and sample from discrete random variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/anaconda/lib/python3.6/site-packages/h5py/__init__.py:34: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEDCAYAAADKhpQUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de5BcZ3nn8W/f5mZd2iaDo9HIgWTx\nkxiHkQ0D2cAmJBAw4ABZyQE7ESknG5LFdoBU5QIbiIGUi7AOixeRWxnsIOyEWIIiCwRCrq6wwVFs\na+xgeENig0casZaxW5Y8l56+7B+nW+qZ6dPn0n26e45+n7LKo+5z3vd53vedRz1vnzmdqdfriIhI\nOmQHHYCIiPSOirqISIqoqIuIpIiKuohIiqioi4ikiIq6iEiK5AfZ+YkTp4b6esotW0Y5fXpl0GEk\nJs35KbfNK8359Sq3ycmtGb/n9Eq9g3w+N+gQEpXm/JTb5pXm/PqRm4q6iEiKqKiLiKSIirqISIqo\nqIuIpEiiV7+YWRG4FbgUqAM/75z7pyT77NZiucqBw/McnFvg5FKF7eN59s5MsW92FxMj4d7k6EUb\nIiJxJP1K/RbgC8657wdmgK8l3F9XFstVrr3zfj5++CilpQp1oLRU4eOHj3LtnfezWK72pQ0RkbgS\nK+pmth34EeCjAM65snOulFR/vXDg8DxHS8uUq7U1j5erNY6WljlweL4vbYiIxJXkK/VnAyeA28zs\nfjO71czOS7C/rh2cW9hQjJvK1RqH5o73pQ0RkbiS3FPPA5cDNzjn7jGzW4DfBN7VPGDLltGh+kWD\nk0uVjs+XllcpFicSb6Nfcrns0MTSa8pt80pzfv3ILcmifhQ46py7p/H3g3hF/Yxh+1Xg7eN5Sh2K\ncnGsQKm0mHgb/VIsTgxNLL2m3DavNOfXq9wmJ7f6PpfY9otz7tvAvJlZ46GXAQ8l1V8v7J2ZYiTX\nfkhGcln2zOzoSxsiInElffXLDcAdZvYAsBu4KeH+urJvdhfTxbENRXkkl2W6OMa+2V19aUNEJK7M\nID94ehjv0ti8xvzQ3HFv/3uswJ6ZHbGuUz80d5zS0irF8eht9IN+zN2c0pwbpDu/Hm6/+N6lUUW9\ngzQvLkh3fspt80pzfv0o6rpNgIhIiqioi4ikiIq6iEiKqKiLiKSIirqISIqoqIuIpIiKuohIiqio\ni4ikiIq6iEiKqKiLiKSIirqISIqoqIuIpIiKuohIiqioi4ikiIq6iEiKqKiLiKSIirqISIqoqIuI\npIiKuohIiqioi4ikiIq6iEiKqKiLiKSIirqISIrkk2zczL4JnAKqQMU594Jetb1YrnLg8DwH5xY4\nuVRh+3ievTNT7JvdxcRIrlfdJMIv9j0zUxyaW4idU2u7paUKo/kMGTKsVGod29rMYykia2Xq9Xpi\njTeK+gucc4+3e/7EiVOxOl8sV7n2zvs5WlqmXK2deXwkl2W6OMZt11zWk2JULE5QKi123U4rv9gL\n2Qx4/1Gunh2WsDn5tdtqfVvF4gQLj53qy1j2WxJzNyzSnBukO79e5TY5uTXj99ym3H45cHi+bfEq\nV2scLS1z4PD8gCIL5hf7aq3OarW+pqBD+Jz82g1qazOPpYhslHRRrwN/ZWb3mtmbe9XowbkF3+JV\nrtY4NHe8V131XKfY/YTJKWy769vazGMpIhsluqcOvMQ5d8zMngl8ycy+7py7u/nkli2j5PPRf7Q/\nuVTp+HxpeZVicSJyu+vlctmetNMqKHY/QTlFabfZVi6X7dtY9lsSczcs0pwbpDu/fuSWaFF3zh1r\n/P8xM/s08ELgTFE/fXolVrvbx/OUOhSj4lihJ/tWSeztBcXuG0tATlHabbZVLE70bSz7Tfuym1ea\n8+vhnrrvc4ltv5jZeWa2tfk18ArgX3vR9t6ZKUZy7UMfyWXZM7OjF90kolPsfsLkFLbd9W1t5rEU\nkY2S3FO/EPhHM5sD/hn4nHPuC71oeN/sLqaLYxuKUfOKjX2zu3rRTSL8Yi9kMxRyGUZya9/UDpuT\nX7tBbW3msRSRjRK9pDFI3Esa4ey11YfmjlNaWqU4XmDPzI6eXlud1I+BfrE3r1OPm1Nru08urTKW\n9wr1SqXWtq1mfv0Yy37Tj/CbV5rz68cljZu2qPdDmhcXpDs/5bZ5pTk/XacuIiKRqKiLiKSIirqI\nSIqoqIuIpIiKuohIiqioi4ikiIq6iEiKqKiLiKSIirqISIqoqIuIpIiKuohIiqioi4ikiIq6iEiK\nqKiLiKSIirqISIqoqIuIpIiKuohIiqioi4ikiIq6iEiKqKiLiKSIirqISIqoqIuIpIiKuohIiuST\n7sDMcsC/AMecc1d209ZiucqBw/McnFvg5FKF7eN59s5MsWdmikNzCxse3ze7i4mRXOzzosTQ6Zw4\nObVrr/XY0lKF0XyGDBlWKrXQcXQbf6/zF2m1WK7yJ3/zDT5xz7cirS+ty7My9Xq94wFmdpVz7q6g\nxzqc/6vAC4Bt64v6iROnOnfeYrFc5do77+doaZlytXbm8UI2A95/lKtnmxvJZZkujvGRvc/juoMP\nRD7vtmsuY+qZWymVFgNjaD0nygKK0p7fsa2C4ug2/l7nn6RicWLN3KVJWnNrrq9jJ5dZqYRfX+fi\nupyc3Jrxey7M9ss7Qj62gZlNA68Bbg1zfCcHDs+3LWirtTqr1fqawgxQrtY4WlrmnZ99KNZ5Bw7P\nh46h0zlxcmrXnt+xUeLoNv5e5y/Sqrm+Wgs6JL+u08a3qJvZq8zsw8BOM/vfLX9uByoh2/8Q8OuA\nfyUK6eDcQseC1k65WuPIsadinXdo7nikGPzO6SRKe2Hz7xRHt/H3On+RVnHXl9blWp321Bfw9sJf\nC9zb8vgp4O1BDZvZlcBjzrl7zeyl7Y7ZsmWUfD7cj0Unl8L+O7JW6P2ddUrLq+RyWYrFidAxlJZX\n1xwfJEp7UfL3i6Pb+Hudf5LWz12apDW3uOtL63It36LunJsD5szsTufcaoy2Xwy81sxeDYwB28zs\nE865n20ecPr0SujGto/nKcUo7BniFfbiWIFqtbZm/ysohuJYIdJ+WZT2ouTvF0e38fc6/ySldd8Z\n0ptb3PV1Lq7Lycmtvs+F2VN/oZl9ycz+zcweNrNHzOzhoJOcc+9wzk07554FvBH429aCHtXemSlG\nctGuwBzJZdm9c1us8/bM7IgUg985nURpL2z+neLoNv5e5y/SKu760rpcK0y1+yjwQeAlwCzelSyz\nSQbVzr7ZXUwXxzZMXiGboZDLMJJb+2Zw853vm668JNZ5+2Z3hY6h0zlxcmrXnt+xUeLoNv5e5y/S\nqrm+RvPR1pfW5VphLmm8xzn3oiQ6j3JJI5y9FvXQ3HFKS6sUxwvsmdlx5nrz9Y+vv0496nntflTy\na6vb69TDtNd67JNLq4w1Fv9KpRY6jjX9La9SHIsWf6/zT0patygg3bktlqvc9eC3ueOeRyOtr3Nt\nXXa6pDFMUX8/kAM+BZzZBHfO3ddtYFGLer+l+ZsH0p2fctu80pxfP4p6mN8obb5Kf0HLY3Xgx7sJ\nSkREei+wqDvnfqwfgYiISPcCi7qZvbvd48659/Y+HBER6UaY7ZenW74eA64EvpZMOCIi0o0w2y+/\n1/p3M7sZ+GJiEYmISGxx7qc+AUz3OhAREelemD31Bzn7m/Y5YBLQfrqIyBAKs6feeg/0CvD/nHPx\n7q4lIiKJCtx+cc59CygCPwn8FHBJ0kGJiEg8gUXdzN4K3AE8s/HnDjO7IenAREQkujDbL78AvMg5\n9zSAmf0u8E/Ah5MMTEREogtz9UsGqLb8vdp4TEREhkyYV+q3AfeY2acbf3893u14RURkyIR5o/SD\nwLXAE40/1zrnPpR0YCIiEl2Y69R/CPhq81a7ZrbNzF7knLsn8ehERCSSMHvqfwCcbvn76cZjIiIy\nZEK9UeqcO/NhFs65GuH24kVEpM/CFOeHzexXOPvq/C1A4AdPi4hI/4V5pf7LwA8Dx4CjeJ+E9OYk\ngxIRkXjC3Hr3MeCNfYhFRES6FOfWuyIiMqRU1EVEUqRjUTezrJn9dL+CERGR7nTcU3fO1czs14E/\nj9qwmY0BdwOjjX4OOud+O1aUPbRYrnLg8DwH5xY4uVRh+3ievTNT7JvdxcRIbs2xT69U+KMvf5OD\ncwuUliqM5jNkyLBSqXU8r11/dx05xsnlKhm8TxzZPpbnqt3tz48SYzc5X//yi7s6f5Cxx20vTh/9\niGGYxVnDMjiZer3e8QAzez/wOPBJWj6E2jn3RMB5GeA859xpMysA/wi81Tn3leYxJ06c6tx5jy2W\nq1x75/0cLS1TrtbOPD6SyzJdHOO2ay47szgXy1X+2yfn+NZ3Ftcc26rdee36m39yidXaxlRHchmm\ni+Mb+g0bY7c5f88zJrj1DTOB/yhFGbN+xR7UXrE4Qam02HUfH9n7PK47+EDPcuqFsLn1Qpw13K1+\n5tdvvcptcnKr700Vw+ypvwG4Du9V972NP/8SdJJzru6ca/4maqHxp69FfL0Dh+c3fHMClKs1jpaW\nOXB4fs2xjz7hX9D9zmvXX7tvBu/8ett+w8YYRqf2Hn1iMbC9qGPWr9jjtBenj3d+9qHEYxhmcdaw\nDFaYG3o9u82f7w3TuJnlzOwI8BjwpUHfL+bg3IJvkS5XaxyaO77m2JWKf0H3Oy9sf536DRtjGJ3a\nW6kEtxd1zPoVe5z24vRx5NhTiccwzOKsYRmsMDf0mgB+FbjIOfdmM3sOYM65zwad65yrArvNrAh8\n2swudc79a/P5LVtGyef796PryaXOH61aWl6lWJwIdazfeVH6i9OvX19+um2vl2PWz9hzuWyovoL6\nCPrRMmpOvRA2t16Is4a71c/8+q0fuYW9n/q9eL9VCt5vlt4FBBb1Judcycz+DrgCOFPUT59eCR9p\nD2wfz1PqsEiLY4Uz+11Bx/qdF6W/OP369eWn2/Z6OWb9jD3s3mVQH803BePEkJR+7jnHWcPd0p56\nsMnJrb7PhdlT/z7n3AeAVQDn3CIhPvnIzCYbr9Axs3HgJ4Cvhwk4KXtnphjJtU95JJdlz8yONceO\n5oOHZ/15Yfvr1G/YGMPo1N5oPri9qGPWr9jjtBenj907tyUewzCLs4ZlsMIU9XKjKNcBzOz7gDAv\nsXcAf2dmDwCH8fbUQ7+6T8K+2V1MF8c2LNLmlQz7ZnetOfaiCyY6Luh257Xrr5Bt/2+gd+XAxn7D\nxhhGp/YuumAisL2oY9av2OO0F6ePm668JPEYhlmcNSyDFeaSxlcA/wO4BPgr4MV4n370d9123u9L\nGuHsNbeH5o5TWlqlOF5gz8yOttfaFsZH2P/X/8ahueM8ubTKWOOV+0ql1vG8dv0dPLJAably5sf5\nYohrvcPE2E3O17/8YlaXyj0ds37FHtRelB9zg/rodU7d6vf2RJw13A1tvwTrdEljYFEHMLNnAD+E\nt+3yFefc411HxWCKehRpXlyQ7vyU2+aV5vz6UdTDXP3yN865lwGfa/OYiIgMEd+i3vg1/wngu8zs\nfM6+OboN2NmH2EREJKJOr9R/CXgbMIV3SWOzqD8F7E84LhERicG3qDvnbgFuMbMbnHMf7mNMIiIS\nU5hLGr9tZlsBzOy3zOxTZnZ5wnGJiEgMYYr6u5xzp8zsJcDLgY9y9kOoRURkiIQp6tXG/18D/LFz\n7nPASHIhiYhIXGGK+jEz+yO8W/B+3sxGQ54nIiJ9FqY4/zTwReCVzrkScAHwa4lGJSIisYS5n/qi\nc+5TwEkzuwjvwy4GemMuERFpL8xvlL4W+D2869UfAy7CK+rPTTY0ERGJKsz2y/vw7vvyb865Z+Nd\nAfOVzqeIiMgghCnqq8657wBZM8s27s74goTjEhGRGMJ88lHJzLbgffD0HWb2GPB0smGJiEgcYV6p\nvw5YBN4OfAH4D+AnkwxKRETiCX29uXOuAvwT8E28m3qJiMiQCVPU7wbGzGwn3icf7QNuTzIoERGJ\nJ0xRzzQ+bPq/Ar/vnLsKXc4oIjKUQhV1M/vPwM9w9tOP+v/BjCIiEihMUX8b8A7g0865r5rZ9wJd\nf+i0iIj0XuAljc65fwD+oeXvDwO/kmRQIiIST6fPKP2Qc+5tZvZ/gPr6551zr000MhERiazTK/UD\njf/f3I9ARESke50+o/Texv//wcwmG1+fCNuwme0CPg5ciPdK/48bn3uaSovlKgcOz3NwboGTSxW2\nj+fZOzPFvtldTIwM3/vKi+Uqf/I33+AT93xrTbx7ZqY4NLfAwbkFSksVRvMZMmRYqdQ25BQl5802\nPn5a8wgan37H08txjdtuL8973aXfDZkMn3nw+KZeM/2Wqdc37KycYWY3AtfjvaGaASrAh51z7w1q\n2Mx2ADucc/c1PuP0XuD1zrmHmsecOHHKv/MhUCxOUCotBh63WK5y7Z33c7S0TLlaO/P4SC7LdHGM\n2665bKgWYTPeYyeXWamcjbeQzYD3H+Vq+6lp5vSRvc/juoMPhMp5EOMTdu6i8MujVT/mvJlbUuMa\nt91enwfeWmxdicP6PRVWr9bl5OTWjN9zvle/mNmvAi8GZp1zFzjnzgdeBLzYzN4e1Klz7rhz7r7G\n16eArwE7owa/GRw4PN92QZarNY6WljlweH5AkbXXjLe1oAOs1uqsVuu+BR3O5vTOzz4UOufNNj5+\n/PJo1c+ckhrXuO32+jzY+GbeZlszg9DpksZ9wNXOuUeaDzSufPlZ4E1ROjGzZwGXAffEiHHoHZxb\n8P1GL1drHJo73ueIOusUbxjlao0jx54KnfNmGx8/YcetXzklNa5x203ivKhtSec3SgvOucfXP+ic\nO2FmhbAdNO7weAh4m3NuzT1jtmwZJZ8f3h+hcrksxeJE4HEnlyodny8tr4Zqp1+C4g0jaN+sNedB\njE/YuYsiyrglOefN3JIa17jtJnVelLaGXRLrcr1ORb0c87kzGsX/EHBH4yPx1jh9eiVMMwMTdv9r\n+3ieUoeFWRwr9Hx/txtB8Yaxfq9zvdacBzE+SeypRxm3JOe8mVtS4xq33aTOi9LWsOvhnrrvc522\nX2bM7Kk2f04BPxjUqZllgI8CX3POfTBy1JvI3pkpRnLth3Ikl2XPzI4+R9RZp3jDGMll2b1zW+ic\nN9v4+Ak7bv3KKalxjdtuEudFbUs6FHXnXM45t63Nn63OuTDbLy/G25f/cTM70vjz6p5FPkT2ze5i\nuji2YWE236nfN7trQJG114x3NL823kI2QyGXYSTn+8b6mZxuuvKS0DlvtvHx45dHq37mlNS4xm23\n1+eB9xNhlLYk4JLGpKXlkkY4e53tobnjlJZWKY4X2DOzY2ivqV0sV7nrwW9zxz2Prom3eZ36obnj\nPLm0ylij8K9UahtyipJzv8cnie0XWJtH0PgkpTW3pMY1brs9OW95leJYgddeeiFkMvzFg9/eFN9T\nYfTjkkYV9Q6SKgzDIs35KbfNK835DfQ6dRER2XxU1EVEUkRFXUQkRVTURURSREVdRCRFVNRFRFJE\nRV1EJEVU1EVEUkRFXUQkRVTURURSREVdRCRFVNRFRFJERV1EJEVU1EVEUkRFXUQkRVTURURSREVd\nRCRFVNRFRFJERV1EJEVU1EVEUkRFXUQkRVTURURSJD/oAJKyWK5y4PA8B+cWOLlUYft4nr0zU+yb\n3cXESG7Q4YUWlMew5BkljtZjS0sVRvMZqMNKtU4GqAOjuQyZTIaVSi10W639/vyPfC8f+/I3uevI\nMU4uV8+0u30sz1W7o7UVZSzb5ZYhw3KldubrZk6vu/S7IZPhMw8eX9PfnpkpDs0tRI4jzhxEGZ8o\nfYUZh3ZjsndmiutffnHksU5qnHqRf7t1meT3a6Zer/ekofXM7GPAlcBjzrlL2x1z4sSpRDpfLFe5\n9s77OVpaplytnXl8JJdlujjGbddcFmoAi8UJSqXFJEIMJSiPj+x9HtcdfCB2nr3KL8p4+x0bJEpb\nhWyGTAbq9TqrbboYyWWYLo6HaivKmomTW7OYtsaO9x/l6tlnWuOYeubWDfMWZw7mn1xitbbxW7Dd\n+ETpy29dhjGSy/I9z5jg1jfMBBbmOPOV5DwH5d9cl+A/t2EL++Tk1ozfc0luv9wOXJFg+74OHJ5v\nu6DK1RpHS8scODw/iLAiC8rjnZ99aCjyjDLefscGidLWaq1Oudq+oHtt1UO3FWUs4+S2vqSu1uqs\nVutrvunDxBFnDtoVdO+cjeMTpS+/dRlGuVrj0ScWA8c77nwlOc9B+TfXZdS5jSqxou6cuxt4Iqn2\nOzk4t+C7oMrVGofmjvc5oniC8jhy7KmhyDPKeHc6Nsig2go7lt3EE0anOHo9B9301WldhrFSCR7v\nuPOV9DzHzb+X368D3VPfsmWUfL73+74nlyodny8tr1IsTgS2k8tlQx2XlKA8gvaugvLsVX5Rxjvo\n2CCDaivMmuk2njBKy6tt5y2JOfDLudt12U3fYWOIG3sv5jlu/mHrUpCBFvXTp1cSaXf7eJ5Sh4Ev\njhVC7SUPek89KI/1+7HrBeXZq/yijHfQsUEG1VaYNdNtPGEUxwpUq7UNsSQxB345d7suwwga77jz\n1Y95jpt/2LoEMDm51fe5VF7SuHdmipFc+9RGcln2zOzoc0TxBOWxe+e2ocgzynh3OjbIoNoKO5bd\nxBNGpzh6PQfd9NVpXYYxmg8e77jzlfQ8x82/l9+vqSzq+2Z3MV0c2zCwzXeZ983uGlBk0QTlcdOV\nlwxFnlHG2+/YIFHaKmQzjOQyFHy68K7uCNdWlLGMk9v6SxgK2QyFnBd/lDjizEEh2/4CinbjE6Uv\nv3UZxkguy0UXTASOd9z5SnKeg/JvrsuocxtV7sYbb+xJQ+uZ2Z8CvwNctH///l/av3//yRtuuOH+\n1mMWF8uJdF7IZXnNJReSAR75ziIrlRrnjxe4+vKdvPsKC33Z0NhYgeXl1SRCDCUoj+3jha7y7FV+\nUcZ7/bHLlRpj+Sy5DFTrZ4vcaD5DIZulWquHbutMv8/fyS1XX0a9UuPhx59muVI7025xPM/Vl0+H\nbyvCmvHLLZ/NUKnVz3zdzOkNl00xM72dR59YWhP7+179A4zksr5xtJu3WHOQIfT4ROlr/br0G4d2\nY3L15Tu5+aoZcgGXWsedryTmOfT3ZWNd1iq1ruoSwHnnjb7H77nErlMPI6nr1Htl0HvqSUtzfspt\n80pzfr3KbVDXqYuISJ+pqIuIpIiKuohIiqioi4ikiIq6iEiKqKiLiKSIirqISIqoqIuIpIiKuohI\niqioi4ikiIq6iEiKqKiLiKSIirqISIqoqIuIpIiKuohIiqioi4ikiIq6iEiKqKiLiKSIirqISIqo\nqIuIpIiKuohIiqioi4ikiIq6iEiK5AcdQK8tlqscODzPwbkFTi5V2D6eZ+/MFPtmdzExkot0zvUv\nv7jnfUl8zfG+68gxTi5XyQB1YPtYnqt2T7FnZopDcwu+89E6X6WlCqP5DBkyrFRqHY8NmtteHNsu\n9tdd+t2QyfCZB4+viXe5Ugsde9zzhm09x5275rHUYaVa77hm2o3V+vNGcxkymbX9NtuIuy57LVOv\n13veaJOZXQHcAuSAW51z7299/sSJUz3tfLFc5do77+doaZlytXbm8ZFcluniGLddc1nbbzK/c77n\nGRPc+oaZtgMfp69hUyxOUCotDjqMUJrjPf/kEqu1jcumkAUyGTJAuXr2+eZ8fGTv87ju4AMb5qtV\n0LHt5jbKOvA7tpDNgPffmtiBMwUiSJQ8u805aevXpd+4tYqTv9+aiaI5d9TrrLbpLmhdxh3Xycmt\nGb/nEtt+MbMc8BHgVcAlwNVmdklS/QEcODzfdjLL1RpHS8scODwf6ZxHn1hse07cviS+5ni3K+gA\nqzVYrdY3fHM25+Odn30o8Bs96Nh2cxtlHfgdu1qrt40dwhX0qHlGOW8Y1rPfuLWKk7/fmomiOXft\nCnqnPpIc1yT31F8I/Ltz7mHnXBn4M+B1CfbHwbkF38ksV2scmjse6ZyVSvtz4vYl8XUa7yDlao0j\nx54KdX7QsevnNso66CaHMKLkGeW8Qa/nsOMWN/9BSWpck9xT3wm0/jN0FHhR6wFbtoySz/fuR7qT\nS5WOz5eWVykWJ7o+p5vzhkkulx36GJuCxjtIlNdiQce2zm2UddBtDmHEfc0ZJeekrV+XUcYtuc3k\nZCQxrgN9o/T06ZWetrd9PE+pwwIojhU27CHHOaeb84bJZtpTDxrvIGH3psMc2zq3UdZBtzmEESXP\nKOf1cz2vX5dRxi1u/oMSd1wnJ7f6Ppfk9ssxYFfL36cbjyVm78wUI7n2KY3ksuyZ2RHpnNF8+3Pi\n9iXxdRrvICO5LLt3bgt1ftCx6+c2yjroJocwouQZ5bxBr+ew4xY3/0FJalyTzP4w8Bwze7aZjQBv\nBP4iwf7YN7uL6eLYhkltvtO8b3ZXpHMuumCi7Tlx+5L4muNdyLZ/07+QhUIuw0hu7fPN+bjpykva\nzleUY9vNbZR14HdsIZtpGzt4rzzDiJJnlPOGYT37jVurOPn7rZkomnNX8OkuaF0mMa65G2+8seeN\nAtxwww21/fv3fwO4A7gB+IRz7lDrMYuL5Z52Xshlec0lF5IBHvnOIiuVGuePF7j68p28+wpre+lQ\np3NuvmqGnM8ln3H6GjZjYwWWl1cHHUYoZ8Y7Aw8//jTLldqZglccz3PN86d536t/gJFc1puP6tr5\n2D5eWDNfy5UaY/ks+WyGaq3e8dhOcxtlHfge+/yda2NvPP6Gy6aYmd7Oo08srYm3UquHjj3ueYNc\nz+vX5fpxizJ3zWNzGajW6bhm2o3V+vNG8xkK2ezZfptzl8/GWpdxx/W880bf4/dcotepB+n1deq9\ntpn2nONIc37KbfNKc369ym0g16mLiEj/qaiLiKSIirqISIqoqIuIpMhA3ygVEZHe0it1EZEUUVEX\nEUkRFXURkRRJ3ScfBTGzjwFXAo855y5tPLYb+ENgDKgAb3HO/XPjuZcCHwIKwOPOuR9tPN7xA0AG\nIUpuZrYd+ARwEd46uNk5d1vjnJ8DfqvR7O845/6kv5ls5JPbDF5uW4BvAj/jnHuq8dw7gF8AqsCv\nOOe+2Hh86OYNouVnZj8BvB4u1MsAAAX7SURBVB8YAcrArznn/rZxzvOB24Fx4PPAW51zA33jLOrc\nNZ6/CHgIuNE5d3PjsaGbuxjr8nnAHwHbgBow65xb7uW8nYuv1G8Hrlj32AeA9zjndgPvbvwdMysC\nvw+81jn3XOCqxuN9/wCQkG4nZG7AdcBDzrkZ4KXA75nZiJldAPw23m2SXwj8tpmd34fYg9zOxtxu\nBX7TOfeDwKeBXwNozMUbgec2zvl9M8sN8bxBhPyAx4GfbDz+c8CBlnP+APhF4DmNP+vbHITbCZ9b\n0weBv2z+ZYjn7nbCr8s83gupX27Uk5cCzfsh9Gzezrmi7py7G3hi3cN1vH85AbYDC42vrwE+5Zx7\ntHHuY43H+/4BIGFEzK0ObDWzDN4riifwXsm/EviSc+4J59yTwJcYgsLgk9vFwN2Nr78E7Gl8/Trg\nz5xzK865R4B/x5uzoZw3iJafc+5+51xzHr8KjJvZqJntALY5577SeJX3ceD1yUffWcS5w8xeDzyC\nl1vTUM5dxNxeATzgnJtrnPsd51y11/N2zhV1H28D/qeZzQM3A+9oPH4xcL6Z/b2Z3Wtmb2o83u4D\nQHb2Ldpo/HLbD/wAXpF/EO/HvRqbK7evcvYb+yrO3urZL4fNlBv459dqD3Cfc24FL5ejLc8Nc35t\nczOzLcBvAOtvWLWZ5s5v3i4G6mb2RTO7z8x+vfF4T+dNRd3z34G3O+d2AW8HPtp4PA88H3gN3ivY\nd5nZxYMJMTa/3F4JHAGmgN3AfjPb1r6JofXzwFvM7F5gK97+cpp0zM/Mngv8LvBLA4itW3653Qj8\nL+fc6UEF1gN+ueWBlwA/0/j/T5nZy3rd+Tn3RqmPnwPe2vj6Lrw9MfD+xfyOc+5p4GkzuxuYaTze\n1w8A6YJfbtcC72/8uPfvZvYI8P14eby05fxp4O/7EmlEzrmv4/1IS+Mf29c0nur0AS2bZd465YeZ\nTePt177JOfcfjYeP4eXUNLT5dcjtRcBeM/sAUARqZrYM3MsmmbsOuR0F7nbOPd547vPA5Xj77D2b\nN71S9ywAP9r4+seBbzS+/gzwEjPLm9kE3oL7GgP4AJAu+OX2KPAyADO7EDDgYeCLwCvM7PzGG6Sv\naDw2dMzsmY3/Z/Gu1vnDxlN/Abyxsc/8bLw3nv6ZzTVvvvk13sD/HN6bcV9uHu+cOw48ZWY/1Hiv\n5E14a3jo+OXmnPsvzrlnOeeehXfV2U3Ouf1sornrsC6/CPygmU003jT9UbyLFXo6b+fcK3Uz+1O8\nV6LfZWZH8a70+EXglsZALwNvBnDOfc3MvgA8gHf50a3OuX9ttHM93iTlgI855766vq9+i5Ib8D7g\ndjN7EO8zAH6j5RXE+/C+iQDe65xb/0ZQ3/nktsXMrmsc8ingNgDn3FfN7M/xLomrANc556qNdoZu\n3iBafsD1wH8C3m1m72489orGG/lv4eylcX9JyxUkgxIxt7acc5VhnLuI6/JJM/sg3vdWHfi8c+5z\njeN6Nm+694uISIpo+0VEJEVU1EVEUkRFXUQkRVTURURSREVdRCRFVNTlnGFmGTP7RzN7VctjVzUu\nWxVJBV3SKOcUM7sU7zdrL8P7PY37gStafiszTpt551ylRyGKdEVFXc45jV9Bfxo4DzjlnHtf4x7y\n1+Hdo/z/Atc752pm9sd4v8o9DnzSOffeRhtH8X69+5V4v/V41wBSEdlA2y9yLnoP3m2VXwV8oPHq\n/aeAH27cdz6P92vo4P0q/gvw7vnzE+vu4f2Yc+4yFXQZJirqcs5p3KDtk8CBxi1rXw7MAv9iZkfw\n7snxfY3Drzaz+4D78G5V3FrUP9m/qEXCOefu/SLSUGv8Ae/eNx9zzr2r9QAzew7eHS5f6Jwrmdkn\n8D4WsOnpvkQqEoFeqYvAXwM/bWbfBWBmz2h8RuY24BTeHfR24O2fiww1FXU55znnHsTbZ/9rM3sA\n+CvgQrwtl4eAr+N9xNiXfRsRGRK6+kVEJEX0Sl1EJEVU1EVEUkRFXUQkRVTURURSREVdRCRFVNRF\nRFJERV1EJEVU1EVEUuT/AwP15KHPmF8yAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pymc3 as pm\n",
    "plt.style.use('seaborn-darkgrid')\n",
    "\n",
    "# Initialize random number generator\n",
    "np.random.seed(123)\n",
    "%matplotlib inline \n",
    "disaster_data = np.ma.masked_values([4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,\n",
    "                            3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,\n",
    "                            2, 2, 3, 4, 2, 1, 3, -999, 2, 1, 1, 1, 1, 3, 0, 0,\n",
    "                            1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,\n",
    "                            0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,\n",
    "                            3, 3, 1, -999, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,\n",
    "                            0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1], value=-999)\n",
    "year = np.arange(1851, 1962)\n",
    "\n",
    "plt.plot(year, disaster_data, 'o', markersize=8);\n",
    "plt.ylabel(\"Disaster count\")\n",
    "plt.xlabel(\"Year\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see here that most of the disasters are earlier in the time series - around mid- to late 19th century. \n",
    "\n",
    "Our question is *when did the safety regulation kick in* and there were probably numerous policy interventions. But we want to know when it all began to kick in. \n",
    "\n",
    "This could be \n",
    "- *when did our advertising campaign begin to kick in*\n",
    "- *when did our new product launch begin to resonate with users*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Occurrences of disasters in the time series is thought to follow a Poisson process with a large rate parameter in the early part of the time series, and from one with a smaller rate in the later part. We are interested in locating the change point in the series, which perhaps is related to changes in mining safety regulations. \n",
    "\n",
    "In our model,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "D_t &\\sim \\text{Pois}(r_t), r_t= \\begin{cases}\n",
    "l, & \\text{if } t \\lt s \\\\ \n",
    "e, & \\text{if } t \\ge s \n",
    "\\end{cases} \\\\ \n",
    "s &\\sim \\text{Unif}(t_l, t_h)\\\\\n",
    "e &\\sim \\text{exp}(1)\\\\\n",
    "l &\\sim \\text{exp}(1)\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "the parameters are defined as follows:\n",
    "* $D_t$: The number of disasters in year $t$,\n",
    "* $r_t$: The rate parameter of the Poisson distribution of disasters in year $t$. \n",
    "* $s$: The year in which the rate parameter changes (the switchpoint).\n",
    "* $e$: The rate parameter before the switchpoint $s$.\n",
    "* $l$: The rate parameter after the switchpoint $s$.\n",
    "* $t_l$, $t_h$: The lower and upper boundaries of year $t$. \n",
    "\n",
    "This model is built much like our previous models. The major differences are the introduction of discrete variables with the Poisson and discrete-uniform priors and the novel form of the deterministic random variable `rate`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![img](http://docs.pymc.io/api/distributions/discrete-5.hires.png)\n",
    "Poisson log-likelihood.\n",
    "\n",
    "Often used to model the number of events occurring in a fixed period of time when the times at which events occur are independent. The pmf ([probability mass function](https://en.wikipedia.org/wiki/Probability_mass_function)) of this distribution is\n",
    "\n",
    "$f(x∣μ)=\\dfrac{e^{-\\mu} μ^{x}}{{x!}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/anaconda/lib/python3.6/site-packages/pymc3/model.py:1331: UserWarning: Data in disasters contains missing values and will be automatically imputed from the sampling distribution.\n",
      "  warnings.warn(impute_message, UserWarning)\n"
     ]
    }
   ],
   "source": [
    "# I have modelled above poisson's distribution with pymc3 lib below\n",
    "\n",
    "with pm.Model() as disaster_model:\n",
    "\n",
    "    switchpoint = pm.DiscreteUniform('switchpoint', lower=year.min(), upper=year.max(), testval=1900)\n",
    "\n",
    "    # Priors for pre- and post-switch rates number of disasters\n",
    "    early_rate = pm.Exponential('early_rate', 1)\n",
    "    late_rate = pm.Exponential('late_rate', 1)\n",
    "\n",
    "    # Allocate appropriate Poisson rates to years before and after current\n",
    "    rate = pm.math.switch(switchpoint >= year, early_rate, late_rate)\n",
    "\n",
    "    disasters = pm.Poisson('disasters', rate, observed=disaster_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The logic for `rate` is implemented here\n",
    "* `rate = switch(switchpoint >= year, early_rate, late_rate)`\n",
    "\n",
    "is implemented using `switch`, a Theano function that works like an if statement. It uses the first argument to switch between the next two arguments."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Missing values are handled transparently by passing a `MaskedArray` or a `pandas.DataFrame` with NaN values to the observed argument when creating an observed stochastic random variable. Behind the scenes, another random variable, `disasters.missing_values` is created to model the missing values. All we need to do to handle the missing values is ensure we sample this random variable as well.\n",
    "\n",
    "Unfortunately because they are discrete variables and thus have no meaningful gradient, we cannot use `NUTS` (No U-Turn Sampler) for sampling switchpoint or the missing disaster observations. Instead, we will sample using a `Metroplis` step method, which implements adaptive `Metropolis-Hastings`, because it is designed to handle discrete values. PyMC3 automatically assigns the correct sampling algorithms. This is one of the best features of `PyMC3` the ability to abstract away these choices from you. \n",
    "\n",
    "Briefly - the NUTS is considered best in class, it's often the most accurate. It doesn't work for discrete problems though and we have to use an older and less-efficient sampler called Metropolis-Hastings. We'll talk about this in more depth, but briefly NUTS is more likely to give less-biased results. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "CompoundStep\n",
      ">CompoundStep\n",
      ">>Metropolis: [disasters_missing]\n",
      ">>Metropolis: [switchpoint]\n",
      ">NUTS: [late_rate, early_rate]\n",
      "Sampling 2 chains: 100%|██████████| 21000/21000 [00:18<00:00, 1117.43draws/s]\n",
      "The number of effective samples is smaller than 10% for some parameters.\n"
     ]
    }
   ],
   "source": [
    "with disaster_model:\n",
    "    trace = pm.sample(10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the trace plot below we can see that there’s about a 10 year span that’s plausible for a significant change in safety, but a 5 year span that contains most of the probability mass. The distribution is jagged because of the jumpy relationship between the year switchpoint and the likelihood and not due to sampling error.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[<matplotlib.axes._subplots.AxesSubplot object at 0x7fd6c67adcf8>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7fd6c5ed84e0>],\n",
       "       [<matplotlib.axes._subplots.AxesSubplot object at 0x7fd6c5f77780>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7fd6c5f3b550>],\n",
       "       [<matplotlib.axes._subplots.AxesSubplot object at 0x7fd6cab90668>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7fd6c5f9ecf8>],\n",
       "       [<matplotlib.axes._subplots.AxesSubplot object at 0x7fd6c980c4e0>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7fd6c678f6a0>],\n",
       "       [<matplotlib.axes._subplots.AxesSubplot object at 0x7fd6c4ef79e8>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7fd6c5f030f0>]],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "ename": "ValueError",
     "evalue": "Image size of 1296x1440000 pixels is too large. It must be less than 2^16 in each direction.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m/usr/local/anaconda/lib/python3.6/site-packages/IPython/core/formatters.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, obj)\u001b[0m\n\u001b[1;32m    330\u001b[0m                 \u001b[0;32mpass\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    331\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 332\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mprinter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    333\u001b[0m             \u001b[0;31m# Finally look for special method names\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    334\u001b[0m             \u001b[0mmethod\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_real_method\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprint_method\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/anaconda/lib/python3.6/site-packages/IPython/core/pylabtools.py\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(fig)\u001b[0m\n\u001b[1;32m    235\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    236\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0;34m'png'\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mformats\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 237\u001b[0;31m         \u001b[0mpng_formatter\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfor_type\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mFigure\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mfig\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mprint_figure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfig\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'png'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    238\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0;34m'retina'\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mformats\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0;34m'png2x'\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mformats\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    239\u001b[0m         \u001b[0mpng_formatter\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfor_type\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mFigure\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mfig\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mretina_figure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfig\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/anaconda/lib/python3.6/site-packages/IPython/core/pylabtools.py\u001b[0m in \u001b[0;36mprint_figure\u001b[0;34m(fig, fmt, bbox_inches, **kwargs)\u001b[0m\n\u001b[1;32m    119\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    120\u001b[0m     \u001b[0mbytes_io\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mBytesIO\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 121\u001b[0;31m     \u001b[0mfig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcanvas\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprint_figure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbytes_io\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkw\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    122\u001b[0m     \u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbytes_io\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgetvalue\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    123\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mfmt\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'svg'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/anaconda/lib/python3.6/site-packages/matplotlib/backend_bases.py\u001b[0m in \u001b[0;36mprint_figure\u001b[0;34m(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, **kwargs)\u001b[0m\n\u001b[1;32m   2054\u001b[0m                         \u001b[0morientation\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0morientation\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2055\u001b[0m                         \u001b[0mdryrun\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2056\u001b[0;31m                         **kwargs)\n\u001b[0m\u001b[1;32m   2057\u001b[0m                     \u001b[0mrenderer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cachedRenderer\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2058\u001b[0m                     \u001b[0mbbox_artists\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"bbox_extra_artists\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/anaconda/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py\u001b[0m in \u001b[0;36mprint_png\u001b[0;34m(self, filename_or_obj, metadata, pil_kwargs, *args, **kwargs)\u001b[0m\n\u001b[1;32m    525\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    526\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 527\u001b[0;31m             \u001b[0mFigureCanvasAgg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdraw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    528\u001b[0m             \u001b[0mrenderer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_renderer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    529\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mcbook\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_setattr_cm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrenderer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdpi\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdpi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/anaconda/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py\u001b[0m in \u001b[0;36mdraw\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    384\u001b[0m         \u001b[0mDraw\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mfigure\u001b[0m \u001b[0musing\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mrenderer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    385\u001b[0m         \"\"\"\n\u001b[0;32m--> 386\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrenderer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_renderer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcleared\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    387\u001b[0m         \u001b[0;32mwith\u001b[0m \u001b[0mRendererAgg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlock\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    388\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdraw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrenderer\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/anaconda/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py\u001b[0m in \u001b[0;36mget_renderer\u001b[0;34m(self, cleared)\u001b[0m\n\u001b[1;32m    397\u001b[0m                           and getattr(self, \"_lastKey\", None) == key)\n\u001b[1;32m    398\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mreuse_renderer\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 399\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrenderer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mRendererAgg\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mw\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mh\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdpi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    400\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_lastKey\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkey\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    401\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0mcleared\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/anaconda/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, width, height, dpi)\u001b[0m\n\u001b[1;32m     84\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwidth\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mwidth\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     85\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mheight\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mheight\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 86\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_renderer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_RendererAgg\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mwidth\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mheight\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdpi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     87\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_filter_renderers\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     88\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: Image size of 1296x1440000 pixels is too large. It must be less than 2^16 in each direction."
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 1296x1.44e+06 with 10 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.traceplot(trace, textsize=9, figsize=(18,len(trace['switchpoint']*2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Results \n",
    "We can see above that there's a 5 year span that's most likely for the safety regulation changes. Between about 1888 and 1893. \n",
    "\n",
    "In your marketing example you may have several days or months - with some uncertainty over them. \n",
    "\n",
    "Hopefully that gives you a nice introduction to **Bayesian Changepoint detection** in **PyMC3**\n",
    "\n",
    "**References**\n",
    "- R.G. Jarrett. A note on the intervals between coal mining disasters. Biometrika, 66:191–193, 1979."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# done"
   ]
  }
 ],
 "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}