ppope/survival_analysis_bug.ipynb

## survival_analysis_bug.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Survival Analysis\n",
    "\n",
    "Original Author: Austin Rochford\n",
    "\n",
    "* This notebook is adapted from [here](https://gist.github.com/AustinRochford/afe6862e622c31494b2f) to depict an issue obtained while running the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import pymc3 as pm\n",
    "from pymc3.distributions.timeseries import GaussianRandomWalk\n",
    "import seaborn as sns\n",
    "from statsmodels import datasets\n",
    "from theano import tensor as T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df = datasets.get_rdataset('mastectomy', 'HSAUR', cache=True).data\n",
    "df.event = df.event.astype(np.int64)\n",
    "df.metastized = (df.metastized == 'yes').astype(np.int64)\n",
    "n_patients = df.shape[0]\n",
    "patients = np.arange(n_patients)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "interval_length = 3\n",
    "interval_bounds = np.arange(0, df.time.max() + interval_length + 1, interval_length)\n",
    "n_intervals = interval_bounds.size - 1\n",
    "intervals = np.arange(n_intervals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/phil/anaconda3/envs/python2/lib/python2.7/site-packages/IPython/kernel/__main__.py:4: VisibleDeprecationWarning: non integer (and non boolean) array-likes will not be accepted as indices in the future\n"
     ]
    }
   ],
   "source": [
    "last_period = np.floor((df.time - 0.01) / interval_length)\n",
    "\n",
    "death = np.zeros((n_patients, n_intervals))\n",
    "death[patients, last_period] = df.event"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/phil/anaconda3/envs/python2/lib/python2.7/site-packages/IPython/kernel/__main__.py:2: VisibleDeprecationWarning: non integer (and non boolean) array-likes will not be accepted as indices in the future\n",
      "  from ipykernel import kernelapp as app\n"
     ]
    }
   ],
   "source": [
    "exposure = np.greater_equal.outer(df.time, interval_bounds[:-1]) * interval_length\n",
    "exposure[patients, last_period] = df.time - interval_bounds[last_period]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "SEED = 5078864 # from random.org"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Applied log-transform to lambda0 and added transformed lambda0_log to model.\n",
      "Applied interval-transform to sigma and added transformed sigma_interval to model.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model:\n",
    "    lambda0 = pm.Gamma('lambda0', 0.01, 0.01, shape=n_intervals)\n",
    "    \n",
    "    sigma = pm.Uniform('sigma', 0., 10.)\n",
    "    tau = pm.Deterministic('tau', sigma**-2)\n",
    "    mu_beta = pm.Normal('mu_beta', 0., 10**-2)\n",
    "    beta = pm.Normal('beta', mu_beta, tau)\n",
    "    \n",
    "    lambda_ = pm.Deterministic('lambda_', T.outer(T.exp(beta * df.metastized), lambda0))\n",
    "    mu = pm.Deterministic('mu', exposure * lambda_)\n",
    "    \n",
    "    obs = pm.Poisson('obs', mu, observed=death)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_samples = 40000\n",
    "burn = 20000\n",
    "thin = 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " [-----------------100%-----------------] 40000 of 40000 complete in 28.2 sec"
     ]
    }
   ],
   "source": [
    "with model:\n",
    "    step = pm.Metropolis()\n",
    "    trace_ = pm.sample(n_samples, step, random_seed=SEED)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "trace = trace_[burn::thin]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the hazard rate for subjects whose cancer has metastized is about one and a half times the rate of those whose cancer has not metastized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.exp(trace['beta'].mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAACKCAYAAAC3gesRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8XHW9//HXpKVpmqVLmq50T/mwKDsUW8pW1iuiIriw\niCCi0Hsv4HJd4aKCokgv6EWvIOD2ExRBUAEpm9CWQilVFqEf6L5C12QmadIkzfz+OCfpJE3mTJNM\nJmnez8ejj8ycbT7z6WROPuf7Pd9vLJlMIiIiIiIiIp2Xl+sARERERERE9hUqsERERERERLqICiwR\nEREREZEuogJLRERERESki6jAEhERERER6SIqsERERERERLqICiyRLDCzE83snb3c53Qz2z9bMYmI\niKSjc5dI11CBJZI9ezvJ3LXAhGwEIiIikiGdu0Q6qX+uAxDZh8XM7FbgQ0ADcBmwBPgRcAawH3Cn\nu99sZt8BZgEHmtl/AX8FfgkcFm73kLt/pfvfgoiI9DE6d4l0klqwRLJnIvCiux8AzAF+CvwXcCBw\nSPjvfDP7N3e/HlgPXODuDwBXAYPd/UDgSOAzZjY9B+9BRET6lono3CXSKSqwRLKnJjzhAPwBOBw4\nG/ipuze4ew3wa+DclH1iAO5+K/CR8HEl8C9gcncFLiIifZbOXSKdpC6CItmztemBu8fNDGAo8D9m\n9j2CE9IA4KXWO5pZOTDHgp0agf2Be7ojaBER6dN07hLpJBVYItkztOmBmQ0JH24BbnL3xyL2vQNY\n7O7nhPvPz06IIiIiLejcJdJJ6iIokj2FZvbh8PH5wCKC7hafM7M8M4uZ2TfN7PRwm3qg6WQ2AvgH\ngJmdBpQDRd0XuoiI9FE6d4l0kgoskex5C/iAmb0FXE1w8+9PgdUE/dLfJLhpuOkK3x+B+83sGuC7\nBN0sXgNmAt8Gvm1mx3XvWxARkT5G5y6RToolk3s73cHeMbNDgYeAOe7+01brTgVuIhgG9HF3vzGr\nwYiIiITMbA5wHMG9Ite4++KUdW2en8zsQuArBFftr3f3x7s9cBER6dGy2oJlZoOAW4G57WxyO/BR\n4HjgdDM7MJvxiIiIAJjZCUC5u08HLgd+3GqTPc5PZjYMuB6YTjCq2ocRERFpJdtdBGuBDwLvtV5h\nZpOAre6+wd2TwGMEk9WJiIhk2yzgYQB3XwoMMbMiSHt+OhV40t13uPt77v6FHMUuIiI9WFYLLHdv\ndPe6dlaPAjanPN8EjM5mPCIiIqHW56At4bK21m0CxhBMwFpoZo+Y2XNmdkp3BCoiIr1LTxqmPRa1\nQTKZTMZikZuJiEju9bYv63TxxoBk+HMYwUSqk4BngQnpDqrzlohIr9IlX9i5LLA20LLFamy4rF2x\nWIzNmxNZDaq3KysrVo4iKEfRlKNoylF6ZWXFuQ4hygZ2t1hB0EK1MWVdW+enauCFsNvgCjNLmNlw\nd9/S3ovovBVNv0vRlKNoylE05ShaV527unOY9hYVobuvBorNbLyZ9Se4Ybi9wTBERES60lzgPAAz\nOxJY7+7VkPb89CRwSjgPUClQmK64EhGRvimrLVhmNg34BVAGNJjZ54F7gRXu/ghwJXA/QdeL+9x9\nWTbjERERAXD3hWb2ipktAHYBs83sEqAi3fnJzP4IvBgu//fcRC8iIj1Z1ufB6mJJNW2mp+bfaMpR\nNOUomnKUXllZsW48Cui8FUG/S9GUo2jKUTTlKFpXnbu6s4ugiIiIiIjIPk0FloiIiIiISBdRgSUi\nIiIiItJFVGCJiIiIiIh0ERVYIiIiIiIiXUQFloiIiIiISBdRgSUiIiIiItJFsjrRsIiISE9lZnOA\n44BG4Bp3X5yy7lTgJqABeNzdb0xZNxB4A/iOu/+6e6MWEZGeTi1YIiLS55jZCUC5u08HLgd+3GqT\n24GPAscDp5vZgSnrrgO2dkugIiLS66jAEhGRvmgW8DCAuy8FhphZEYCZTQK2uvsGd08Cj4XbExZa\nBjyak6hFRKTHy3oXwYguGLOBCwm6YCx29y9mOx4REemdzOx9BK1OD5vZEHev6MThRgGLU55vCZct\nC39uTlm3CZgcPr4FmA1c2onXFhGRfVhWC6zULhjhVb97gOnhuhLgy8Bkd0+a2RNmdqy7L8pmTCIi\n0vuY2bXAp4B8gpan68xse+q9UZ0Ui1pnZhcDz7n7GjOL2qdZWVlx56PbxylH0ZSjaMpRNOWoe2S7\nBatFFwwzG2JmRe5eBewEaoESM6sGCoBtWY5HRER6p08R9IZ4Onz+FeAFoKMF1gaClqomY4CNKetG\np6wbGy77N2CymX0M2B+oNbO17v5MuhfavDnRwRD7hrKyYuUognIUTTmKphxF66oCNNv3YLXuZtHU\nBQN33wncACwHVgIL3H1ZluMREZHeKeHujU1PwseNabaPMhc4D8DMjgTWu3t1eOzVQLGZjTez/sDZ\nwFx3/5S7T3P3DwC/AL4bVVyJiEjf093DtDd3pzCzYoKRmKYCCeAZM3ufu7+R7gBq2oymHEVTjqIp\nR9GUo2613Mz+GxhqZucCnwDe7OjB3H2hmb1iZguAXcBsM7sEqHD3R4ArgfuBJHCfLgCKiEimsl1g\npeuCcRCw3N23A5jZfOBogrlF2qWmzfTU/BtNOYqmHEVTjtLLQvE5G7gaWA9cBMwH7ujMAd39G60W\nvZ6ybj7hPcPt7Pvtzry2iIjsu7LdRbDdLhjAKuAgM8sPnx9NMHqTiIhIa7uAOcCHCM4rtwH1OY1I\nRESkDVltwYrqgmFmtwB/N7N64IXwiqGIiEhrDQTd9ZokgUqgNDfhiIiItC2jAsvMYuFki3stogvG\nXcBdHTmuiIj0He7e3OPCzAYQjFJ7WO4iEhERaVumXQRXm9mNZjY5elMREZHscfc6d38cOC3XsYiI\niLSWaRfBaQQjNv3WzHYQTBj8R3evy1pkIiIiITO7rNWicQTzU4mIiPQoGRVY7r6R4Ibi28zsMOC3\nwE/M7GfAje5em8UYRUREZqY8TgJx4OM5ikVERKRdGQ9yYWYnA5cBxwO/B34FfBB4gGBUJxERkaxw\n90u7+phmNgc4jmDC4mvcfXHKulOBmwgG13jc3W8Ml/+Q4DzYD7jZ3f/U1XGJiEjvlukgF8uBlcCd\nwGXu3jQ07ltm9pFsBSciIn2bma2l5eiBLbj7+A4e9wSg3N2nm9mBBF3fU+e9up3gHq+NwHNm9keC\neR0PCfcZBvwDUIElIiItZNqCdQYQc/d3AMzsCHf/R7huZvu7iYiIdMrxadYN7cRxZwEPA7j7UjMb\nYmZF7l5lZpOAre6+AcDMHgu3/xmwKNy/AhjUmVF2RURk35RpgfUZYAxBF0GAr5vZCnf/mk4sIiKS\nLe6+uumxmR0MDA+f5gM/Bg7q4KFHAYtTnm8Jly0Lf25OWbcJmOzujcCOcNnlwGM6B4qISGuZFlgn\nu/uMpifu/vFw8mDZhzQ2NlJRsT3XYWRdXl4d27Ylch1Gj9YbcjRkyFDy8jKdaUJ6OzO7HTid3UVQ\nOfCjLnyJWKbrzOzDwKVhPCIiIi1kWmANMLMBTcOym1nRXuwrvURFxXZ23nIzQwYOzHUo2TVoAAU7\nNMNAWj08RxW1tVR85WsMG1aa61Ck+xzj7geZ2bPufrKZHQWc34njbSAo1pqMIbjfqmnd6JR1Y8Nl\nmNkZwNeBM9w9o6sQZWXFnQizb1COoilH0ZSjaMpR98i0SPo/ggEtFhOMnHQMcEO2gpLcGTJwIKUF\nBbkOI6sKB+UzMNkv12H0aL0hRzW5DkC6W0P4Mz+87+mVcBTAjppLcB67y8yOBNa7ezUE3RLNrNjM\nxhMUVmcDF5hZCfBDYJa7V2b6Qps39+zW4FwrKytWjiIoR9GUo2jKUbSuKkAznQfrbjN7kqCwSgLX\nuvvaTPaNGAZ3f+A+YD9gibtftZfxi4hI3/CWmf078DzwpJk5UNLRg7n7QjN7JezuvguYbWaXABXu\n/ghwJXA/wTnvPndfZmafA0qBP5hZLFz3aXdf17m3JiIi+5JMh2kfCBxBcDKLAaeZGe5+T8R+UcPg\n3grc4u5/NrOfmNn+OlGJiEhr7v55MxsKVAKfAkYA3+/kMb/RatHrKevm0/J8hbvfBdzVmdcUEZF9\nX6ZdBJ8guMK3OmVZkqBgSifdMLgxguF3Pxmu/4+9ilxERPoMM3sR+DVwv7v/v1zHIyIi0p5MC6z9\n3P3EDhw/3TC4ZUAVcFvY/31eG1cTRUREAL4EfAL4h5n9E/gN8OemwZdERER6ikwLrH+ZWam7b+3k\n68VaPR4L/A+wBnjUzM5y98fTHUCjn0TraI7y8upg0AAKB+V3cUQ9T2Hhvv8eO6sn56g2tovC4cWU\nlub2+0DfR93H3RcAC8zsauBE4CKCiX/LchqYiIhIK5kWWPsDy8zsLXaP5IS7nxCxX7phcLcAq9x9\nFYCZPQ0cAqQtsDT6SXqdGSFm27YEBTvqevzocZ1VWJhPdfXOXIfRo/X0HFXX1FGzJUFj44CcxaDR\nmNLLRvFpZkOAjxAMzz4Z+HmXv4iIiEgnZVpg3dzB46cbBneXma0wsynuvhw4CvhdB19HRET2YWb2\nBMFFuD8BN7n7CzkOSUREpE15mWzk7s8BRcD7w8frCIbKjdpvIdA0DO5thMPgmtmHw02uBX5pZvMJ\nhsb9S0fehIiI7PNuBya4+3+ouBIRkZ4s02HafwBMBSYA/wtcQDBEbuTIfxHD4C4HZmYarIiI9E3u\n/liuYxAREclERi1YwInufi4QB3D37wJHZi0q6fG21NZyzlNPND/fsGMHx/7lT+xKJgGob2zkQ0/9\njR0NDeH6aj7w14dZ8N67zfs8u3ED//bk41zxwjw+t+B5Lpn3dxZt3gRAor6Oo/78EP/YuqV5+2c2\nrueOt/4FwOyF83mzYnvkdpV1dXzv1X/wuQXPc/mC57lq4XyWx+Np39stb7zKo2vXAFDT0MAXXpjH\nn1avBOCUv/21Od6vvPwiG3ZUN+937aKFfHres1zxwjyuWPA8V7wwj6c3rM8on5fO+zubamraXf/3\njRv2WLZtZy2ffv7Z5vfaWkNjI9ctWZzR8VP9bOmb/HnVqhbLViTiXPDcM9y/cnlGx2jLG9u3cdbc\nx/neq/9ga20t//niAi6b/xxff2URyWSSuevXce4zc7nT3+rwa4jsDTObY2YvmNl8Mzu61bpTzewl\nM1tgZt/KZB8RERHI/B6spr/MkgBm1m8v9pV90ND8fBL19c3PH1q9kglFxVTW7WRY/kDmrl/HjBGj\nGNQ/+Jjc6Us5fez+vBOvZMbIYNyT5Yk4F0wu56IpUwHYuGMHly94nkdPO5Pl8QSl+QN5eM0qjigd\nDsCyeJwpxSUArKxKMKmoGK+sbHe7nbt2MfvF+VxSfgDfOOwIAOa9t5H/XrSIn3+g/YbTFfEE54yb\nQHVDPde+tJBzxk/g7HET2FRTw6iCAu6cHuz70uZNXLvoRX534in0i8VYkYjzm5knUzJg7wde2Lyz\nlhEFBW2u25VM8pvl73DS6DEtlj+7cSMnjhrNZw84sM39+ufl8d0jj448fmtvxyv50JRJLZY9um4N\nF0+Zyln7j8voGO0d98Ipwf/3Fxct5JLyAzhqeBm3vP4q8ze9y+lj92d5Is7UksEdfg2RTJnZCUC5\nu083swMJ5nVMnVj4duA0goGZnjOzPxL03Ei3j4iISMZF0gtmdi8wxsy+CJwL/D1rUUmP1y8Wo18s\nxq5kkmQyyWvbtnLEsFIq6uoYlj+QB1ev5FthUbOqKsGKRJyvHXo4v13+TvMxViYSnLn//s3PRw8a\nxI6GBpLJJCuq4nxi0mT+uGolVfX1FO23H8vicU4aNYaq+nr6xWIU9O+fdrsHVq3g0KGlnDZm92vM\nHDma0yaOp7amnte2bePlLZv2KFBWVycoGziQ/3zxBc6fOJkzw6JiRSLOpKKS5u2mlY1gZEEBL2/e\nxBGlw6moq4ssri547hlmjBjJP7dtJRaLccdxM9hcW8vogkEA/GXNah5cvZIkcPyIkXzODuL6JYtZ\nt6Oany19kysPPBiAf27dys+WvsnIggIOG1bKq9u2Mv+9d9nZ2Mi5EyZy3sTJPLR6JZV1dZwxdlzz\n8R9du4YHw9a4Y4aXceWBB1PdUM/1SxYTr6/n0GHDWJVIMKWkhJ01QQH99Ib1PLByBeUlg7HBg3lo\n9SrerNhOQ2MjV9hBHD9yFD956w0q6+rYuGMHd3zgeL788ot8/6hj2S9vdyP5O/FKTho1hg07qqms\nq+Oo4cHo2lNLBrMqUcXMkUERds64CVEfP+mDzOww4G6gyN0PNLPrgLnu/lIHDzkLeBjA3Zea2RAz\nK3L3KjObBGx19w3haz8KnEowJHyb+3Ty7YmIyD4kowLL3b9pZucBOwiGbJ/j7g9lNTLp8Ybl51NZ\nt5MlW7dywqjRVNbVUVFXxzvxSvbLy2Ny2Nr0f0vf4uLyqYwrLGRZSve85Yk4Ewp3D+X8xPp1TB08\nmFgsxvJEnJkjRrF9TB2PrVvDxydNYUVVnEnFxbxZsb352Om2u/Vfr3FVWJCk6hf+0X/osGEcOmxY\ni3WJ+joak0m+tOhFThw1urm4Cl4rwaTilkNPjykYxObaWlZVJRg5MH0LUTKZZF11FSeNOpzZBx3C\nt5a8zIub3yOZDIqM5fE4D61eyd3Hn0i/WIxP/f1pPjR+AtNHjGRcYSFfSHkvh5eWMr6oiJuPOpbS\ngQNZkYhz78yTqG9s5Nxn5nLexMksj8eZVjaCZfFKppYMZn11Nc9sXM89xwdzhn963rOcO2Eif1y1\nkiNKh3PRlKk8tHol++Xl0T8vj6ZB2meNGcvtb77BL2acwENhcXbP8SeyfedOLp3/HMePHMXKRIIj\nS4fzrcOCnsM/Oua4Pd7/O/E4VxxwEAs2vcvRw3dPXVTd0MCAfsH/yfrqasYWFqbNo/RZ/wtcRtCy\nBPB74F5gRgePNwpYnPJ8S7hsWfhzc8q6zcAUoDTNPm166Y2NVMYz657bVw1+r0o5iqAcRVOOoilH\n6Q0rHthlU4xkOsjFZGBJ+K95mbuv6JIopFcamp/P9p11/HnNKr575DE8smYVFXU7eWL9Oj4+cTIA\nb1dW8Fbldr5/1DHEYjEq6+uob2ykXyzG2uoqvvvqEmLAmuoqpo8YyQ+OOhYIWrc+U34AwwcO5Lol\ni/nw+IkA7JeXx4pESlfBNNttqqlpbrmJ19Xx5ZdfZFcySXH+AG475gNtvqcViQS7kklOGT2WuRvW\ncUn5AcRiwfzYK6riTC8b2WL7d2t2MHPUKFYkEmzduZMrXpgHySTEYswaPYZPTJrSvO3a6mrKSwZz\nyNCgqCvo15+6XY2srq6ivKSEv61fy0cnTKJf+HpjBg1iU00tyxNxDhkydI9Y36upYURBAeuqq1my\ndStPb9zAzl27GBy2oi1LxLlwSjl/W7+O8pISHl+/lrXV1c0x1jbsojEJz767gbumB1PajR1UuEcX\nvUR9PYX9+5MXi/HYurXcGHY7HJqfT3VD0Mq1qirR3B2xPRV1Oxman8/KqgTlxbtfY3VVglljxpII\nWyBF2lHv7q+ZGQDu/raZNUTsszdiHViXbp9mg0sy657blylH0ZSjaMpRNOWofUOHdF1uMu0i+DTh\n/VdAPkE/9DeAI7osEul1huXn81ZlBYMH5DN4wABKBgxgc20tS7Zu4avvPwyAO5a+SX1jIxc+/wzJ\nJOxoaGBlIs7Afv2ZUFTMXTOCP+y//PKLnDNuIkPz8wHYVFtD2cACygYWkN+vH39du4aJRcFVheWJ\nBAcOHhK53dD8fCrr6xhRUEDJgAHcOeMEXt++jQfXrmr3PS1PxDln3AQuLp/K2/FKfr9yBZ+cHBRJ\nK+JxLpxc3rxtRd1OllZWcMSw4dzzjnPB5Cnt3g8FQRe58uLdXQxXViW4cEo5T29cz6cml7N4y2aO\nGzEyZX0V44sKWZaI8+HxLbvNbaqpYWR4T9V3Xl3CJeVTmTFiFM+9u5FnwwExNu7YwZhBhSyLV/Kp\nyeU8tHol1x9+JO8burvVrjGZpGJnXXPe34lXUl5S0uK1lsfjzS2G63dUM2ZQ0MK0ubaG0vyB1O7a\nRb9YHoX92y+ONu7Y0VzsJpMwqP/uyaxf276Naw95P29WbGdqq9cWSdEQdt1ruhf4LDIscNqxgaD1\nqckYgvutmtaNTlk3FlgP7EyzT5umvW+0JqSOoEm7oylH0ZSjaMpR98l0HqxJ7j45/DcWOBx4Nruh\nSU83bEA+D6xawbkTJgJQst9+PLZuLWeO3Z+8WIzXtm1jbXU1fz31TH534izuO2kWZ48bz9vxOCsS\n8RZ/TH9o3AQeWBU0iCbq6xm83+57mT4yYSJ3vf1Wc6tV0IJVHLndrNFj+eU7bzePbFi7axcPr1nF\nYaWl7b6nFeH9RwBXH/w+fr38bbbW1gKwurqK8YVFQNAidv2SxXy6/AAG9e/PykScg9poZUq1LBHn\n3dqgaX5NVRXV9fVMLCpmWTxOeXEJYwYVsjxeCcCLm95jQlERQwbksyqRYP/wdZssT2nFe69mBwcN\nHkLdrl08uHolU4pL2L5zJ0Pyw5as8PhDB+SztLICgJe3bOaJ9WvJC1vLGhobqWlo4JE1q/cssBLx\n5mVjBxWyIhF08/zDyhWcPW48KxNxJhWlb1J/O17JAYODVqsJRUXNXUUfXbuGI0pLKejfn7fDrowi\n7fgy8Agww8wqgZuB/+zE8eYC5wGY2ZHAenevBnD31UCxmY03s/7A2eH2T7a3j4iISJMOjQTo7v8y\ns6O6OhjpXYbl57OjoaF59L6S/QbwdmUFtx0bdL/76dJ/cc3B72v+Ix6gvHgw78QrGTxgQIs/po8f\nOYpb3niVirqdrK6qai4eAM4Ysz9z3ni9edmqqgSTikvwyoo2t2tqbfnk5Cnc6W9x2fznyM/LoyGZ\n5Kyx47jogANYvbWSny19k+sObznbwIpEnHPGjQdg+MCBXDRlKj/612t88ZD3s6OhgateXAAEl83P\nnTCJM8YGA2gsTyT4P3+TXy57u7mL4PkTJ7UYYOOdeCWHDh3GZ+c/R32ykesOP5L6xkYakkkK+vfn\nwsnlfHPJyzy9cQMF/fo13880qbiYb7yyiB8cPa35WMsT8eZC8KLJU7ly4QImFBZxwshR3L9yOZNL\niikvLmlx/E9NnsI3lyxm7vp1DOzfn28fHvwKf2ziJC5+/lkmF5dQ0K8fk4taFlgrquJMGz4CgGsP\neT83/PMVBub1Y1JxMVcdeDCPrlvbHAsEw9x/ctIUxqUUhe/EKzmgJGh1PGPsOL66+CWuWPA8ZQML\nmv8P3o5X8pGwi6dIa+7+GnComZUBO909/XwL0cdbaGavmNkCYBcw28wuIZj0/hHgSuB+ghaz+9x9\nGbCs9T6diUFERPZNsWQyGbmRmX2n1aJxwJHuflhWompfUk2b6XWm+Xfbtq0U/OQ2SjMczru3KizM\np7p6Z/SGXezcZ+Zy34mzyO/XL3rjHOvqHH118Ut8wQ7eY5CQVBc9/wx3Tj+heWj/dLbW1FDzH9cw\nbFj7rZHZpq4W6ZWVFXem+14zM/sNu7uo78HdP90Vr5NFOm9F0O9SNOUomnIUTTmK1lXnrkxbsHal\nPE4CrwLfamfbFsxsDnAc0Ahc4+6L29jm+8Bx7n5yhvGI9Co1DQ30j+X1iuIqG96JV3LTa//gR8dM\nY8iA/BbrNtfW8M1XXqamYVdGxZX0OU/lOgAREZG9kelfMzfSxhVEM8sDcPfGtnbKYCJHzOwgYCZQ\ntxdxi/QqBf3784eTT811GDnz0Cmnt7uubGABd4aDnYi05u6/anpsZu8DDiY4H73m7p6zwERERNqR\n0SAXQA1Q38a/hvBne1pM5AgMMbOiVtv8CPj6XsQsIiJ9jJndAvwJ+AjwMeAxM/tubqMSERHZU6Yt\nWN8G3iQYRSkJfAg4yN1viNgv3USOhDcUPw2syTxkERHpg04BDnb3egAzywdeAK7LaVQiIiKtZFpg\nneLuN6U8/72ZPdOB12u+cczMhgIXA6cD48lwPpOummF5X9bRHOXl1cGgARQOyo/euJcrLNz332Nn\n9eQc1cZ2UTi8mNLS3H4f6PuoW71L0GuiSR2wMkexiIiItCvTAqvUzP4NeD58PhMYnsF+6SZyPAUY\nCcwHBgKTzexWd/9SugNq9JP0OjeKYIKCHXUMTO7bAzHkahTB3qSn56i6po6aLQkaGwdEb5wlGo0p\nvSwUn1uAl8OLe3nACcDyplFu3f36rn5BERGRjsi0wLoCuJVgThCAN4CrMthvLnADcFcbEzk+CDwI\nYGYTgHujiisREemzVoT/mjza0QOFkwf/EphA0Cp2qbuvarXNhcDVBKPo3uXu95hZP+BuYArQD/iy\nu7/Q0ThERGTflFGB5e6LgJlmFnP36Imzdu8XNZGjiIhIJHf/dhce7gJgu7tfZGanATcDn2xaaWaD\nCO7tOpqgAHvZzB4iGGBjh7vPNLODgXuBaXscXURE+rSMCiwzO4zgql0RcKCZfQt40t1fitrX3b/R\natHrbWyzmqDLoIiIyB7M7OvAfwEl4aIYkHT3jvRpngU0Df/+FMEUIqmmAYvcvSp87fnADOC37O7J\nsRkY1oHXFhGRfVymXQT/F7gMuD18/geCK3czshGU5E5FbW2uQ8i62tguqms07Vo6PT1HFbW19Nwh\nOCRLPg0cDqzrgmONIiiQcPekmTWaWX93b2i9PrQZGB2ub9rmGuB3XRCLiIjsYzItsOrd/TUzA8Dd\n3zazhoh9pJcZMmQoFV/5GjW5DiTLCocXU7NFgxOk09NzlE/weZU+5V/AOnfftTc7mdlngcsJphiB\noOXr2FabRc0J2WKUWzObDRxBMGWJiIhIC5kWWA1mNonwBGVmZ5HhsOrSe+Tl5TFsWGmuw8i60tLi\nnI4+1xsoR9ID/QZ43cwWkzJcu7tflm4nd7+boIt7MzO7h6CV6vVwwAtSWq8gGAF3dMrzscDCcN/P\nAh8EPpxpsafh/KMpR9GUo2jKUTTlqHtkWmB9CXgEMDOrBFYRdNcQERHpDrcSFFld0UXwSeD88Oc5\nwLOt1r/TrInmAAAOXElEQVREMPptCdAITAeuNrPJwOeBE5omPM6EhvNPT1MeRFOOoilH0ZSjaF1V\ngGZaYG1190PNrAzY6e7xLnl1ERGRzCzrwpEEfw+cZmbzgFrgMwBm9lXg7+7+kpl9jWCqkUbgBndP\nhMuGAY+ZWYygV8fprVq/RESkj8u0wPp/wMnuvjlySxERka73kpl9G1hAyy6Cz+ztgdy9kWDgptbL\nf5Dy+CHgoVbrvwl8c29fT0RE+pZMCyw3s18DLwDNQ4u5e+uhbUVERLLhhFY/IWhB2usCS0REJJvS\nFlhmdqi7v0YwaNcught7t4Srk+w5d4iIiEiXc/eTWy8zs4/lIhYREZF0olqwbgNOcfdLAczsGXfX\nsLQiItKtzGw88O/A8HBRPsEE9Q/mLCgREZE27NXcHyIiIjnya2Ab8AHgFWAEGs1WRER6oKgWrGSr\n53tdcJnZHOA4gpGYrnH3xSnrTga+R3DDsrv75Xt7fBER6RMa3P1mMzvT3e8ws7uBBwiGWhcREekx\nolqwWmtdcKVlZicA5e4+Hbgc+HGrTX4OfMzdZwIlZnbmXsYjIiJ9Q6GZTQAaw/mo6oH9cxyTiIjI\nHqJasKab2ZqU5yPC5zEg6e7jI/afBTwM4O5LzWyImRW5e1W4/uiUObU2A6V7Gb+IiPQNPyAYQfAW\n4J8EAy/9LqcRiYiItCGqwLJOHn8UsDjl+ZZw2TKApuLKzEYDpwHf6uTriYjIPsjdH256bGbDgGJ3\n396RY5lZf+CXwASCLuqXuvuqVttcCFxNUMjdlTotiZmNBN4CPuLuz3ckBhER2XelLbDcfXUXv94e\n93CZ2Qjgz8CVmZwsy8qKuzikfY9yFE05iqYcRVOOss/MSoDL3X1O+PzzwJXAMjOb7e7vdeCwFwDb\n3f0iMzsNuBn4ZMprDgKuA44mKMBeNrOH3L0i3OSHwPIOvykREdmnZTrRcEdtIGixajIG2Nj0xMyK\ngceAr7v705kccPPmRJcGuK8pKytWjiIoR9GUo2jKUXpdWHz+HFgDYGYHAN8HPg5MBm4npTDaC7OA\nX4WPn2LPOR2nAYuaurOb2XxgBvBoODhTJfB6B15XRET6gL0d5GJvzQXOAzCzI4H17l6dsn4OMMfd\nNQqUiIi0ZbK7fzV8fB7wgLs/5e530vIC3t4YRXDfL+6eJBg4o39b60ObgdFmth9BV/ZvomlMRESk\nHVltwXL3hWb2ipktIOjHPtvMLgEqCIqvi4ApZvY5ghEKf+fuv8hmTCIi0qtUpTw+Cbg75Xlj1M5m\n9lmCUWybRsGNAce22izTOSG/BvzM3RNmlro8LXUljaYcRVOOoilH0ZSj7pHtLoK4+zdaLUrtVlGQ\n7dcXEZFerX94r24xwSTDn4Dme7OKonZ297tpWZRhZvcQtFK93tRy5e4NKZtsAEanPB8LLAQuAc4y\nsy8BU4BjzOx8d38rXQzqSpqeuttGU46iKUfRlKNoXVWAZr3AEhER6YSbgTeBQcAN7r7dzAqAecCd\nHTzmk8D54c9zgGdbrX8JuCss4hqB6cDV7v5Y0wZmdi9wb1RxJSIifU+278ESERHpMHd/nKA1aZS7\n/zBcVgN8xd3v6OBhf0/QMjaPYETCrwOY2VfNbJq71xJ0B5wb/rvB3Vtf9k0iIiLShlgy2avOEUk1\nbaan5t9oylE05SiacpReWVmxBoEI6LwVQb9L0ZSjaMpRNOUoWledu9SCJSIiIiIi0kVUYImIiIiI\niHQRFVgiIiIiIiJdRAWWiIiIiIhIF1GBJSIiIiIi0kVUYImIiIiIiHQRFVgiIiIiIiJdpH+2X8DM\n5gDHAY3ANe6+OGXdqcBNQAPwuLvfmO14RESkbzOz/sAvgQkE559L3X1Vq20uBK4GdgF3ufs94fIv\nAxcCdcBV7v5K90UuIiK9QVZbsMzsBKDc3acDlwM/brXJ7cBHgeOB083swGzGIyIiAlwAbHf3mcD3\ngJtTV5rZIOA64BTgZOBaMxtiZgcDHweOBD4PnN2tUYuISK+Q7RasWcDDAO6+NDxBFbl7lZlNAra6\n+wYAM3ss3H5plmMSEZG+bRbwq/DxU8A9rdZPAxa5exWAmc0nuBB4MPAHd08C/wz/iYiItJDte7BG\nAZtTnm8Jl7W1bhMwOsvxiIiINJ9/wmKpMew2uMf60GaC89NEYIKZPW5mT5rZod0Ur4iI9CJZvwer\nlVgH1zVvU1ZW3FWx7LOUo2jKUTTlKJpy1POZ2WcJuqgnw0Ux4NhWm0VdbIyF+8eAPHc/y8xmAL9o\n41h77KvPSTTlKJpyFE05iqYcdY9sF1gb2N1iBTAG2JiyLrXFamy4TEREpEu4+93A3anLzOwegnPT\n600tV+7ekLJJW+enheHPpeH2C8xsQhZDFxGRXirbXQTnAucBmNmRwHp3rwZw99VAsZmND09wZ4fb\ni4iIZNOTwPnh43OAZ1utfwk42sxKzKwImA7MA/4GnAkQDsq0tnvCFRGR3iSrLVjuvtDMXjGzBQRD\n3c42s0uACnd/BLgSuJ+g68V97r4sm/GIiIgAvwdOM7N5QC3wGQAz+yrwd3d/ycy+RnDRrxG4wd0T\nwEtmdpaZvUBw3pqdk+hFRKRHiyWTyeitREREREREJFK2uwiKiIiIiIj0GSqwREREREREuogKLBER\nERERkS7S3fNgZczM5gDHEdxgfI27L05ZdypwE9AAPO7uN+YmytyJyM/JwPcI8uPufnluosytdDlK\n2eb7wHHufnJ3x9cTRHyO9gfuA/YDlrj7VbmJMrcicjQbuJDgd22xu38xN1HmVjjh7kPAHHf/aat1\nffb7OpPvoL7CzH4IHA/0A24GXgZ+Q3ChdyNwsbvXm9mFwNUEA2Pd5e735CjknDCzgcAbwHeAZ1CO\nWgjf+1eAeuB64HWUo2ZmVgj8GhgKDCD4HL2JcrTHeSr8GyejvISjnf8SmEBwLrvU3Vele70e2YJl\nZicA5e4+nWCCyB+32uR24KMEX9anh8Pl9hkZ5OfnwMfcfSZQYmZndneMuZZBjjCzg4CZ7J6AtE/J\nIEe3Are4+3HArvDLqE9JlyMzKwG+DMxw9xOAQ8wsatLZfY6ZDSL4rLQ3zUaf/L7O5DuorzCzk4BD\nwlycBdxG8Iff/7r7icBy4LLws3QdcApwMnCtmQ3JTdQ5cx2wNXz8HeAnylHAzIYRFFXTCab2+QjK\nUWufAZa6+ykEU1Hcjn7X2jtP7c1n5wJge/h39fcILhKl1SMLLGAW8DCAuy8FhoRzkWBmk4Ct7r7B\n3ZPAY+H2fUm7+Qkd7e5NkzZvBkq7Ob6eICpHAD8Cvt7dgfUg6X7PYgR/EP8lXP8f7r4uV4HmULrP\n0U6CIb5LwqtbBcC2nESZW7XAB4H3Wq/o49/XmXwH9RXPs3vesQqgEDgR+HO47C/AacA0YJG7V7l7\nLTAfmNHNseaMmRlgwKNAjCBHfwlXK0dwKvCku+9w9/fc/fPASShHqTax+2++YQR/A+p3re3z1Elk\n9tk5nuD7/E/htk+RQa56aoE1iuBD0WRLuKytdZuA0d0UV0+RLj+4exzAzEYTfGAe69boeoa0OQrn\nY3saWNPNcfUk6XJUBlQBt5nZPDP7XncH10O0myN33wncQHDlayWwoC/O5efuje5e187qvvx9nfY7\nqC8JPyM7wqefJSggCt29PlzW9LkYScucbabvfF4guOj3RYLiCpSj1iYChWb2iJk9Z2anAIOUo93c\n/QFgnJm9QzCB+pfQ56i989Te5KV5eXixsDG8sNqunlpgtRbr4Lq+Yo8cmNkIgisWV7r79u4Pqcdp\nzpGZDQUuJuimEkOfoSaxVo/HAv9DcPXrCDM7KydR9Sypn6Nigq4EU4FJwAwze1+uAusl+vLvWl9+\n7wCY2YeBy4B/Z8/vm7b0mZyZ2cXAc+7e3kW/Pp8jgvc6jKDL8aXAvehz1EJ4/9Bad59K0OpyR6tN\n+nyO2rG3eYmsn3pqgbWBllf6xhDcgNa0LrXKHhsu60vS5afpD7/HgG+4+9PdHFtPkS5HpxBcjZhP\ncMPjEWZ2a/eG1yOky9EWYJW7r3L3RoLWvkO6Ob6eIF2ODgKWu/t2d28g+Dwd3c3x9XR9+fs67fd0\nX2NmZxB0yT7T3RNAwszyw9VjgfX07c/LB4HzzWwhQSvfdUCVctTCe8ALYWvECkCfoz3NAJ4AcPfX\nCd57tXLUpkw/O03LRwE0tVyF5/129dQCay5wHoCZHQmsd/dqAHdfDRSb2fjwTZ5N+zdX76vazU9o\nDsEoKU/mIrgeIt1n6EF3f394w/VHCUbI+1LuQs2ZdDnaBawwsynhtkcBnpMocyvd79oq4KCUL+ij\ngT7XRbCVFlf7+vj3ddT3dJ8RDgjzQ+Bsd68MFz8FfCx8/DHgb8Ai4GgzKwnvV5sOzOvueHPB3T/p\n7tPc/QPALwhuwH+K8DOEcgTB79QpZhYzs1KgCOWotWUEI5diZhMIuvo/iXLUlr35DnqS3feRnkPQ\n/TKtWDLZMwdQC+/5OJFgmMTZwJFAhbs/YmbHE3xZJ4E/uvv/5C7S3GgvPwRfQNuAhQR/7CSB37n7\nL3IUas6k+wylbDMBuDcccafPifg9m0IwLGkMeN3dr8xZoDkUkaPPEXR5qie4svq13EWaG2Y2jeAP\nwjKC4Wu3EXTdWdHXv69bf3bCK8p9Tvh78t/A2+w+L10C3A3kA6sJhj3eZWbnAv9FMLT9j939/txE\nnTtm9t8E93U+QTCMtHIUCj9LlxN8hr4LLEY5ahYO034PQS+dfsC3CC6O/po+nKN2zlNnAL8ig7yY\nWV64/1SCATM+4+7r071mjy2wREREREREepue2kVQRERERESk11GBJSIiIiIi0kVUYImIiIiIiHQR\nFVgiIiIiIiJdRAWWiIiIiIhIF1GBJSIiIiIi0kVUYImIiIiIiHSR/w8g9pl1irneUgAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb74b427ed0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#pm.traceplot(trace, vars=['beta']);\n",
    "pm.traceplot(trace, varnames=['beta']);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuMAAACrCAYAAADM6ZKMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFrJJREFUeJzt3X+0ZWV93/H3hUGygEsRcyMDRvxB+mHV/KjUCA4IyABq\nRPEXaRZoEcWqYAua+CNRFBVclEQUYuxKkBG1UfwRFVshYWgQQRDEaotRvhENaBkqFyp2AInMzOkf\ne189c+bM3DN37r175p73a627zj77efbZ31nry+F7nv3sZ0/0ej0kSZIkLb6dug5AkiRJGlcW45Ik\nSVJHLMYlSZKkjliMS5IkSR2xGJckSZI6YjEuSZIkdcRiXJLGQJIjknxvK485NsnjFiomSZLFuCSN\nk619sMQbgP0XIhBJUmNZ1wFIkhbNRJL3Ac8H1gGvBP4H8GfAs4FdgL+qqvOSvBtYCRyY5M3AfwMu\nBX6n7fe5qnrT4v8TJGlpcWRcksbHE4CvVdW/BC4APgS8GTgQeEr7d0KS36uqdwB3ASdW1WeA04B/\nUVUHAgcBr0iyooN/gyQtKRbjkjQ+ftYW1gCfBv41cBzwoapaV1U/Az4GvLjvmAmAqnof8MJ2+6fA\nPwBPWqzAJWmpcpqKJI2P+2Y2qur/JQF4NPD+JO+lKbwfBdw0eGCSA4AL0hy0AXgcsGoxgpakpcxi\nXJLGx6NnNpLs1W7eC5xbVVfMcuxfALdU1Qva469fmBAlabw4TUWSxsfuSY5vt08AbqaZrvLqJDsl\nmUjytiTHtn0eAWaK9l8DvgmQ5BjgAGCPxQtdkpYmi3FJGh/fBZ6R5LvAGTQ3ZX4IuJNmDvh3aG7m\nnBn1/ixwWZIzgffQTFP5X8AzgXcB70pyyOL+EyRpaZno9bZ22dltl+QC4BCaeYdnVtUtfW1HA+fS\nLLt1ZVWd0+4/CXgTzUjNO6rqykUPXJIkSZpHiz4ynuRw4ICqWgGcClw00OVC4EXAYcCxSQ5Msjfw\nDmAFzZ3/xyNJkiTt4Lq4gXMl8AWAqrotyV5J9qiqB5I8EbivqtYAJLmi7T8NrK6qh4CHgNd2ELck\nSZI0r7qYM74PTXE9495237C2e4B9aR5UsXuSy5Ncm+SoxQhUkiRJWkjbw9KGE7O09drXvWkeOPFE\n4Bpg/9k+uNfr9SYmtvTxkiRJ0jabc8HZRTG+hl+OhEMz8n13X9vyvrb92n0PAjdUVQ/4QZK1SX61\nqu7d0okmJiaYnl47f5FrSZiamjQvtAnzQoPMCQ1jXmiYqanJOR/bxTSVq4CXAiQ5CLirqh4EqKo7\ngckkj0+yjOZmzauA1cBR7Rq4jwF2n60QlyRJkrZ3iz4yXlU3JvlGkq8C64HTk5wM3F9VlwOvAy6j\nmZ7yyaq6HSDJZ4Gvtftfv9hxS5IkSfOtk3XGF1HPS0ka5CVGDWNeaJA5oWHMCw0zNTU55znjPoFT\nkiRJ6ojFuCRJktQRi3FJkiSpIxbjkiRJUkcsxiVJkqSOWIxLkiRJHbEYlyRJkjpiMS5JkiR1xGJc\nkiRJ6ojFuCRJktSRZV2cNMkFwCHABuDMqrqlr+1o4FxgHXBlVZ3T1/YrwLeBd1fVxxY3akmSJGl+\nLfrIeJLDgQOqagVwKnDRQJcLgRcBhwHHJjmwr+0s4L5FCVSSJElaYF1MU1kJfAGgqm4D9kqyB0CS\nJwL3VdWaquoBV7T9aYvyAF/qIGZJkiRp3nVRjO8DTPe9v7fdN6ztHmB5u/2nwBuBiYUOUJIkSVoM\nncwZH7Cl4noCIMnLgWur6odJZjtmI1NTk9sWnZYk80LDmBcaZE5oGPNC86mLYnwNvxwJB9gXuLuv\nbXlf237tvt8DnpTkJcDjgIeT/Kiq/n62k01Pr52XoLV0TE1NmhfahHmhQeaEhjEvNMy2/EDrohi/\nCjgbuDjJQcBdVfUgQFXdmWQyyeNpivDjgBOr6kMzByd5J/BPoxTikiRJ0vZs0YvxqroxyTeSfBVY\nD5ye5GTg/qq6HHgdcBnQAz5ZVbcvdoySJEnSYpjo9Xpdx7CQel5K0iAvMWoY80KDzAkNY15omKmp\nyTkvMOITOCVJkqSOWIxLkiRJHbEYlyRJkjpiMS5JkiR1xGJckiRJ6ojFuCRJktQRi3FJkiSpIxbj\nkiRJUkcsxiVJkqSOLNuazkkmgF88YaiqNszlpEkuAA4BNgBnVtUtfW1HA+cC64Arq+qcdv/5wGHA\nzsB5VfX5uZxbkiRJ2l6MNDKe5E1J7qcpkB/pe91qSQ4HDqiqFcCpwEUDXS4EXkRTeB+b5MAkRwJP\naY95LvCBuZxbkiRJ2p6MOk3llcBvV9XO7d9OVbXzHM+5EvgCQFXdBuyVZA+AJE8E7quqNVXVA65o\n+38FOKE9/n5gt3aUXpIkSdphjTpN5XtV9cN5Ouc+wC197+9t993evk73td0DPKmdDvNQu+9U4Iq2\nWJckSZJ2WKMW47cm+QTwZZopKgBU1ap5iGFLI9wbtSU5HjgFOHbUD5+ampxjWFrKzAsNY15okDmh\nYcwLzadRi/F9gX8GntG3rwfMpRhfQzMC3v/Zd/e1Le9r26/dR5JnA38MPLuq1o56sunpkbtqTExN\nTZoX2oR5oUHmhIYxLzTMtvxAG6kYr6pTAJLsDfSq6idzPiNcBZwNXJzkIOCuqnqwPc+dSSaTPJ6m\nCD8OODHJnsD5wMqq+uk2nFuSJEnaboxUjCdZAXwcmAQmktwHvKx/ScJRVdWNSb6R5KvAeuD0JCcD\n91fV5cDrgMtoRt4/WVW3J3k18Bjg0+2Nmz3g31XV/97a80uSJEnbi4leb/b7IJN8BTitqr7dvn8q\ncGFVHb7A8W2rnpeSNMhLjBrGvNAgc0LDmBcaZmpqcs6r/I26tOH6mUIcoKq+Sd+NnJIkSZK23qg3\ncG5I8mLg6vb9c2immEiSJEmao1FHxl8L/HvgTuAO4OR2nyRJkqQ5GnU1le/RjIZLkiRJmidbLMaT\nXFhVZyS5jmYFk43sADdwSpIkSdut2UbGZx7q8/aFDkSSJEkaN1ssxqvqf7abp1TVK/rbkvwdcO0C\nxSVJkiQtebNNUzmJ5kbN32zXGp+xCxs/0l6SJEnSVpptZPyvk3wZ+GvgnX1NG4B/WMC4JEmSpCVv\n1tVUquou4Mj+fUl2AT4BnLAwYUmSJElL30hLGyZ5GfB+YO921wbgv8/1pEkuAA5pP+fMqrqlr+1o\n4FyaJ3xeWVXnzHaMJEmStCMa9aE/ZwC/BVwH7Am8Hvj4XE6Y5HDggKpaAZwKXDTQ5ULgRcBhwLFJ\nDhzhGEmSJGmHM2ox/tOq+j/AzlX1YFX9JXDKHM+5EvgCQFXdBuyVZA+AJE8E7quqNVXVA74EHL2l\nYyRJkqQd1UjTVIANSV4A/CjJ2TQ3bz5ujufcB+ifYnJvu+/29nW6r20aeDLwmC0cs1lPeAJs2LD7\nHMPUUrXTTuaFNmVeaJA5oWHMCw3zwx/O/dhRi/GXAcuBM4FzgKcC/2Hup93IxBzatnTMRnbaadTB\nf40T80LDmBcaZE5oGPNC82m2dcZnsu3e9g+adce3xRo2XqN8X+DuvrblfW37AXcB/7yFYzbrjjtg\nenrttsSqJWhqatK80CbMCw0yJzSMeaHhJud85Gw/7dYBj7R/6/rez7zOxVXASwGSHATcVVUPAlTV\nncBkkscnWQYc1/ZfvbljJEmSpB3VbA/9mffrMFV1Y5JvJPkqsB44PcnJwP1VdTnwOuAyoAd8sqpu\nB24fPGa+45IkSZIW20Sv15u1U5JHA38C7FNVL0/yfOBrVTU9y6Fd63kpSYO8xKhhzAsNMic0jHmh\nYaamJke+n3HQqCPfHwZ+BDypfb8r8NG5nlSSJEnS6MX4VFVdBPwcoKo+C+y2YFFJkiRJY2DkOeFJ\ndqGZx02SxwIusilJkiRtg1HXGf8g8HVgeZIvAk8HzliwqCRJkqQxMFIxXlWfTnID8AyaNb9fU1Wz\nrvMtSZIkafNGKsaTfLaqXgp8ZoHjkSRJksbGqNNUbk/ySuAG2ps4AarqBwsSlSRJkjQGRi3G/+2Q\nfT1+udShJEmSpK00ajF+WFXdtaCRSJIkSWNm1KUN/8uCRiFJkiSNoVFHxivJx9h0zviqrT1hkmXA\npcD+wDrglKq6Y6DPSTRLJ64HLq6qVUl2Bi4BngzsDPxRVd2wteeXJEmSthejjozvSlMYHww8s/07\nbI7nPBH4SVU9E3gvcF5/Y5LdgLOAo4BnAW9IshfwcuCh9rhTgffP8fySJEnSdmHUdcZPAUiyN9Cr\nqp9swzlXAh9tt68GBkfXDwZurqoH2nNeDxxKM1XmsrbPNLD3NsQgSZIkdW6kkfEkK5J8H7gN+Mck\ntyV52hzPuQ9NMU1V9YAN7dSVTdpb08DyqlpXVQ+3+84EPjHH80uSJEnbhVHnjJ8HHF9V3wZI8lTg\nQuDwLR2U5FU0U0p67a4J4OkD3Wb7QTAx8JmnA08Fnj9K4FNTk6N005gxLzSMeaFB5oSGMS80n0Yt\nxtfPFOIAVfXNJOtmO6iqLqG56fIXkqyiGf2+dWZEvKr6P2sNsLzv/X7Aje2xrwKeR/PDYP0ogU9P\nrx2lm8bI1NSkeaFNmBcaZE5oGPNCw2zLD7RRi/ENSV4CrG7fP4fmhs65WA2c0L6+ALhmoP0m4OIk\newIbgBXAGUmeBLwGOLyqHpnjuSVJkqTtxqjF+GuBPwc+TFMgfwt49RzP+SngmCTXAQ8DrwBI8hbg\ny1V1U5K3Ale15zq7qta2+/YGrkgyQTP15diBUXVJkiRphzHR6/Vm7dTO0z62qo5v318D/E1VfXCB\n49tWPS8laZCXGDWMeaFB5oSGMS80zNTU5MTsvYYbdZ3xlwEv7nt/LHDSXE8qSZIkafRifOeBGyY3\nLEQwkiRJ0jgZdc74F5PcAFxHU8CvBP5mwaKSJEmSxsBII+NVdQ7wZuAe4G7gtKo6dyEDkyRJkpa6\nUUfGqarrgesXMBZJkiRprIw6Z1ySJEnSPLMYlyRJkjpiMS5JkiR1xGJckiRJ6sjIN3DOlyTLgEuB\n/YF1wClVdcdAn5OAM4D1wMVVtaqv7bHAd4EXVtVXFilsSZIkad51MTJ+IvCTqnom8F7gvP7GJLsB\nZwFHAc8C3pBkr74u5wPfX6RYJUmSpAXTRTG+Evh8u301cOhA+8HAzVX1QFU9TLOc4qEASZ4F/BS4\ndZFilSRJkhZMF8X4PsA0QFX1gA3t1JVN2lvTwPIkuwBvB94GTCxSrJIkSdKCWdA540leBZwK9Npd\nE8DTB7rN9oNgpvB+K/Cfq2ptkv79WzQ1NTlasBor5oWGMS80yJzQMOaF5tOCFuNVdQlwSf++JKto\nRr9vnRkRr6p1fV3WAMv73u8H3AicDDw3yR8CTwZ+N8kJVfXdLcUwPb12m/8dWlqmpibNC23CvNAg\nc0LDmBcaZlt+oC36airAauCE9vUFwDUD7TcBFyfZE9gArADOqKorZjok+QjwkdkKcUmSJGl71kUx\n/ingmCTXAQ8DrwBI8hbgy1V1U5K3AlfRFONnV9XgT9AekiRJ0g5uotdb0nVtz0tJGuQlRg1jXmiQ\nOaFhzAsNMzU1OefFRXwCpyRJktQRi3FJkiSpIxbjkiRJUkcsxiVJkqSOWIxLkiRJHbEYlyRJkjpi\nMS5JkiR1xGJckiRJ6ojFuCRJktQRi3FJkiSpI8sW+4RJlgGXAvsD64BTquqOgT4nAWcA64GLq2pV\nu/+PgJOAnwOnVdU3Fi9ySZIkaX4tejEOnAj8pKpeluQY4DzgD2Yak+wGnAU8jaZY/3qSzwH7Ar8P\nHAT8DnA8YDEuSZKkHVYXxfhK4KPt9tXAqoH2g4Gbq+oBgCTXA4cB/wr4dFX1gG+1f5IkSdIOq4s5\n4/sA0wBtYb2hnbqySXtrGlgOPAHYP8mVSVYn+e1FileSJElaEAs6Mp7kVcCpQK/dNQE8faDbbD8I\nJtrjJ4Cdquq5SQ4FPjzkszY5dmpqcuuC1lgwLzSMeaFB5oSGMS80nxa0GK+qS4BL+vclWUUz+n3r\nzIh4Va3r67KGZiR8xn7Aje3rbW3/rybZfwFDlyRJkhZcF9NUVgMntNsvAK4ZaL8JeFqSPZPsAawA\nrgP+FngOQJIDgR8tTriSJEnSwujiBs5PAcckuQ54GHgFQJK3AF+uqpuSvBW4CtgAnF1Va4Gbkjw3\nyQ0001ZO7yB2SZIkad5M9Hq92XtJkiRJmnc+gVOSJEnqiMW4JEmS1BGLcUmSJKkjXdzAueCSXAAc\nQnMD6JlVdUvHIalDSc6neYrrzsB5wNeBj9P8GL0beHlVPdJdhOpCkl8Bvg28G/h7zImxl+Qk4E3A\nI8A7gFsxL8Zakt2BjwGPBh5F833xHcyLsdU+dPJzwAVV9aEkj2NIPrTfJ2cA64GLq2rwifO/sORG\nxpMcDhxQVStoHjh0UcchqUNJjgSe0ubDc4EP0HyZfrCqjgC+D7yyuwjVobOA+9rtdwN/bk6MryR7\n0xTgK4DjgBdiXqhZ8e22qjqKZlnmC/H/IWMryW7A+2hW/JuxyfdE2+8s4CjgWcAbkuy1uc9dcsU4\nsBL4AkBV3Qbs1a5XrvH0FX65rv39wO7AEcAX233/FTi6g7jUoSQBAnyJ5um+R9DkApgT4+poYHVV\nPVRVP66q1wBHYl6Mu3uAx7TbewPT+P+QcfYw8Dzgx337jmTj74ljgIOBm6vqgap6GLgeOHRzH7oU\ni/F9aP5jmXFvu09jqKo2VNVD7dtX0RRfu/ddUryHjZ/4qvHwZ8AbaQpxMCcETwB2T3J5kmuTHAXs\nZl6Mt6r6DPDrSb5H85DCP8Tvi7HV1hQ/H9g9LB8ey8a16DRbyJOlWIwPmpi9i5a6JMfTXEp8PRvn\nhPkxZpK8HLi2qn64mS7mxHiaoBn5fBFwCvAR/K4Ye+283x9V1W/QXHn/i4Eu5oX6bS4ftpgnS7EY\nX8PGI+H70kyo15hK8mzgj4HntE9zXZtk17Z5P5qc0fh4HnBCkhtprpacBTxgToy9HwM3tCNfPwD8\nrhA0Uwv+DqCqbqXJgwfNC/UZ/J64iyYn+kfCt5gnS7EYvwp4KUCSg4C7qurBbkNSV5LsCZwPHFdV\nP213Xw28pN1+CfC3XcSmblTVH1TVwVX1DODDNDffXE37vYE5Ma6uAo5KMpHkMcAemBeC22lWZyPJ\n/sADwGrMC/3SsJriZuBpSfZs71tcAVy3uQ+Y6PV6Cx7lYkvyXpobLNYDp7e/ZjWGkrwaeCfwjzSX\niXrAycAlwK7AncApVbW+syDVmSTvBP6JZuTr45gTY639vjiV5nviPcAtmBdjrV3acBXNHOCdgbcD\nRbPcoXkxZpIcTDOIMwWsA/4v8GzgowzkQ5IXA2+mWWb7oqq6bHOfuySLcUmSJGlHsBSnqUiSJEk7\nBItxSZIkqSMW45IkSVJHLMYlSZKkjliMS5IkSR2xGJckSZI6YjEuSWMoyRFJNvsQCknS4rAYl6Tx\n5YMmJKljy7oOQJLUnSSHAn8KPATsBpxWVd9K8hs0T598GPg08IGqelR3kUrS0uTIuCSNrwlgb+B1\nVXU0cBHwJ23bu4BLq+pI4Oc0jwKXJM0zi3FJGl894B7g/CTXAm8FfrVt+y3g+nb78x3EJkljwWJc\nksbXBM1UlPdW1RHA2/radgI29PWTJC0Ai3FJGm+/Bnwnyc7A7wO7tvu/C/xuu318F4FJ0jiwGJek\n8dUD/hNwDfAl4FLg15P8R+A9wBuTrAb2AtZ3FaQkLWUTvZ4rW0mSNpbk3wDLquqmJE8HVlXVb3Yd\nlyQtNS5tKEka5mfAJUnWAbsAp3UcjyQtSY6MS5IkSR1xzrgkSZLUEYtxSZIkqSMW45IkSVJHLMYl\nSZKkjliMS5IkSR35/+1o9JaKGfXVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb6c7f15c90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.autocorrplot(trace, varnames=['beta']);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}