examples/sampling/nested-ellipsoidal-sampling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ellipsoidal nested sampling\n",
    "\n",
    "This example demonstrates how to use ellipsoidal nested rejection sampling [1] to sample from the posterior distribution for a logistic model fitted to model-simulated data.\n",
    "\n",
    "[1] \"A nested sampling algorithm for cosmological model selection\", Pia Mukherjee, David Parkinson and Andrew R. Liddle, [arXiv:astro-ph/0508461v2](https://arxiv.org/abs/astro-ph/0508461)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First create fake data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xd4VFX6wPHvO6mQ0BIC0kOHCKEH\nEAWkSRNRQUBUFIS19766Lrq71rX/RFBUVFZAFEEF6V0EE0B6r4FAQoCQhPSc3x/3AgEHCJDJTSbv\n53nmydwzZ2beOxfmnVPuuWKMQSmllDqXy+kAlFJKFU2aIJRSSrmlCUIppZRbmiCUUkq5pQlCKaWU\nW5oglFJKuaUJQimllFuaIJRSSrmlCUIppZRbvk4HcCUqVqxowsPDnQ5DKaWKlZiYmCPGmLCL1SvW\nCSI8PJzo6Ginw1BKqWJFRPbmp552MSmllHJLE4RSSim3NEEopZRyq1iPQbiTlZVFbGws6enpTodS\nbAUGBlK9enX8/PycDkUp5SCvSxCxsbGUKVOG8PBwRMTpcIodYwyJiYnExsZSu3Ztp8NRSjnI67qY\n0tPTCQ0N1eRwmUSE0NBQbYEppbwvQQCaHK6Qfn5KKfDSBKGUOo+dCyDuT6ejUMWEJohCcu+997Jp\n06YCf93g4OACf03lpXJz4Lu7YcowyMl2OhpVDHjdIHVR9dlnnzkdgirp4tZCepJ1W/8dNB/idESq\niNMWhAekpqbSp08fmjVrRpMmTZg8eTKdO3c+vSzI+PHjadCgAZ07d2bkyJE89NBDANx999088sgj\nXHPNNdSpU4epU6cCkJKSQteuXWnZsiVNmzZl+vTpju2bKsZ2Lbb+htSBJW9pK6IgTbsfxnaEzT+D\nMU5HU2C8ugUx+qeNbDp4okBfM6JqWV6+8eoL1vn111+pWrUqv/zyCwBJSUmMGTMGgIMHD/Lqq6+y\nevVqypQpQ5cuXWjWrNnp58bFxbFs2TK2bNlCv379GDBgAIGBgUybNo2yZcty5MgR2rVrR79+/XQw\nWV2aXYug0tXQ+TmYcids/AEibyvY90g7BinxENawYF+3KDsRB+smgU8ATB4KVZpDz9eg1jX5e/7x\nfbDxR2g2GIIrnb+eMVDI/+c92oIQkT0isl5E1opItF0WIiJzRWS7/beCXS4i8oGI7BCRdSLS0pOx\neVLTpk2ZN28ezz77LEuXLqVcuXKnH1u1ahWdOnUiJCQEPz8/Bg4ceNZz+/fvj8vlIiIigsOHDwPW\nuQkvvPACkZGRdOvWjQMHDpx+TKl8yUqDfb9Dnc7QqC9UioDFb1rjElfq+D6YOBDebghvhMP/RcG2\n2WfXSYmH95vD9rlX/n5FzfrvwOTCqEVw0//ByaMw8TZIPnTx58b9CZ91g7kvwfvNYN4/refnlZsD\nyz+A12vB0v8WagulMFoQ1xtjjuTZfg6Yb4x5XUSes7efBXoB9e1bW2CM/feyXeyXvqc0aNCAmJgY\nZs6cyfPPP0+PHj1OP2YucnADAgL+UnfixIkkJCQQExODn58f4eHhep6CujT7V0JOBtTpBC4XdHrG\nGrDeOA2aDriy1577D9izDK6+GcIawcqx8NuH0OCGM3VWfQrHdsOqcVC/++W/lzEw+wVY8w34lbJu\nlZvAjR9AUOiZeukn4OguqNr88t8rv/H8+S1UbwOVGlm3mu3h43bW53LLuLPrpydBQFmrJbBzAUy+\nEwLLw+3fwbrJsOw967Nq2Asa94OQ2vDzExC7yuoanP8KHN8Pvd8GH89/fTvRxXQT0Nm+PwFYhJUg\nbgK+Mta34u8iUl5Eqhhj4hyI8YocPHiQkJAQ7rjjDoKDg/nyyy9PPxYVFcXjjz/OsWPHKFOmDN9/\n/z1Nmza94OslJSVRqVIl/Pz8WLhwIXv35mulXqXO2LUIXL5nuj0a3wQVG8LvY64sQRxabyWZjk9D\nlxftQmN9OR5aD1c1hcyT8Men1vvvmA+pRyCo4oVf1xj46VHwDbS6a1w+Vvmab+D3j6Fhb6s7JiMF\nNv8E47vB0KkQWhf2LIdp90HSPrh3AVRvdfn7l5/9j98Eff57piy0LlzzsPVrv9XdmJrtMZknke9H\nINtmYVx+EFwZUg5BxYbkDJmCKVsV6naDax7DtWosrq2/IOu/sz6KwPLk3DSWnCa34rPwX/j+9h45\nSQdgwOf4BJbx3L7h+QRhgDkiYoCxxphxQOVTX/rGmDgROdXpVg3Yn+e5sXZZsUsQ69ev5+mnn8bl\ncuHn58eYMWN46qmnAKhWrRovvPACbdu2pWrVqkRERJzVBeXO0KFDufHGG2ndujXNmzenUaNGhbEb\nypvsWmz9yg2wv1BcLmv8YcGrcOIglK3q/nnGWDfXeXqjF/4HAstB+4fOlLW8Cxa9Tu6KMSTf8D7m\njy8pn3aMA9f8i2q/vcj2hV+xM/x2MrJzyMjKtf5m55KRnUtmdi6ZObk0iZtGn70TAFizdSdfX/Uc\noel7eXrfE+wIbMFbaY+QfdJFdo6hXrn2PH1sNHzUieX+HbghYzaHXFdRWsqy84vHearUK+QaY91y\nOXPfnNo9Q1BuMo3NTqJpTCZ+p8vNqY8AY/8F8my/4PMVd7h8iPqhDEk//HL6Iwg0VzM3oCLJ40cx\nNPPvjPV/h1aynXE5fcjFRaVjx0k1jXlr3yCS3/gTyHtuSi986EGUawuRsotp6dcSP7kMTJ4DRDHU\nZzivbP+CP394m5a3j768fw/5JBfr8riiFxepaow5aCeBucDDwAxjTPk8dY4ZYyqIyC/Aa8aYZXb5\nfOAZY0zMOa85ChgFULNmzVbn/prevHkzjRs39tg+FYSUlBSCg4PJzs7m5ptvZvjw4dx8881Oh3WW\n4vA5qnxKOwZv1IZOz8L1z58pT9hqjRf0fhuiRp4pX/01/PaB1ReedgxTrjrxt88nPtOfIykZJKZm\nciw1k4D4Ndy1YTgzQoczudQgTqRlk5yeRXJ6Nk9kjWWALOTajPf5zv8VjlGGmzNH86v/c6QSyK2Z\n5/9iq+2TwM++z7JB6rHK1ZyHc79hkU97apv9BJtUHi77ASl+Ifi6BF+XC5cLquXG8dSRl6iSHcuS\n4N5MDXuA65JnMvDIx3xY/R12BLfERwQRwSUQlJNMjYyt1Dy5kQYpq6iRuhEXufxR8WZm1noaQRAB\nweoNEhEEQMBlDxT75GZz/5q+HCjTnJ8avcG5w8f1ji6k35ZnSfcti29OGrMbvsqOit3OGmfO+5yz\nyi8yGF35xHoatexIkxqhF6x3PiISY4xpfbF6Hm1BGGMO2n/jRWQaEAUcPtV1JCJVgHi7eixQI8/T\nqwMH3bzmOGAcQOvWrYvlfLJ//vOfzJs3j/T0dHr06EH//v2dDkl5sz3LAGMNUOcV1hAqNiBn0ww2\nV7uN/UdPcjAxiduW/YNUKc0aV2uO4suQY78w84MHGZ097Kynf+P/LsdcZRmb0YMAnxxCg/0JrxhE\nmUBfUnJHELB+Lj9XHkvlpMMkRD3PZ+FtCNw2mEZr32LO3TXwDa1DgJ8P/j4uAvxc+Pu48HeB66sb\n4ZA/be+fRNvyNWBFBJ1nvwDigrum87/aHd3vZ3pXOLqLjlVb0BEgqw18OJ2Hzf9g0HDrG3jbbJj9\nd0jcfuZ5VVtAqyfhRBxt1n5Dmx5Dzh4/OZ9ts2HVMRr0GMmTjdzM2jINYOI8AvetgNun0ufcz/+K\n1CvA1zo/jyUIEQkCXMaYZPt+D+AVYAYwDHjd/ntqUv8M4CERmYQ1OJ1UHMcf8uPtt992OgRVkuxa\nBP7BpIY1Z/Oeo2w5lMy2w9bthqNNuDPhR+74cBbHKUN/1zJG+B/ltaBHOVjxGqqUC2RdfCnuPjSV\nxjeMxL9Wa0KD/KkUO5dS09ZDj3/zyzW93LxpU0jrTuUdc6F8Ldr0vMsaR6h6F6x9iwaHZ0Ojp89U\nTz8BR3bDhu9h73K46WMob/9ebP8glLnK6us5X3IAq6uraosz236BVqvpp0escYq4tda4QKWroevL\nVt2qzaFUBat+doZVZ/qDcP8KCL7AJZuNgZgJUCoE6p1n0F0EBk+0xkmCLu+XvtM82YKoDEyzm0q+\nwP+MMb+KyB/AFBEZAewDTs3znAn0BnYAJ4F7PBibUl4tOyeXLYeSWb03kV5rZ7IztyFDXl1weoZk\ncIAv9SsHkxTeE99dPzDh2qP4NO9Jo5/fhOwG/OfBR8/0eaT/Fz5aSrtNr8I1C2DtN9bMmqsioc2I\n8wfR/gHYMdcanzg1yFy+JtS8BtZPgeuegHVTYOnbkLjjzPMa3wjNbz/7tZrcenkfRPPbYfn71owt\nkwMth0GvN63kcS7fALjlUxjXGWY8DB0ehW2/Wi2wet2sgXgfXys5zHkRtv4CnZ4DX//zv79vgHUr\npjyWIIwxu4BmbsoTga5uyg3woKfiUapYMQZOJp5/tk92pjW9cuMPcPXNmBZ3sflQCst3HOG3nUdY\ntfsoqZk59HMt5y7/g/xY6R4ea9+AJtXK0rhKWaqUC7T6uY2Bd1+hWfIyMG3g0FprTCJvH3hgOWsm\n0dR7YMKNsO8361fzwC+saabnU7cLjFwAVVqcXR45EH5+HD5sCcf2WDOdur5szf6pUBsqX11wJ4T5\n+EGPf1mtghv+/dfEc67KEdDtnzD7edg2y5p5FdYIFr9utcRu/RRW/B+s/ASiRlknHXoxrz6TWqli\na/GbsOg1GPS19Yv6lNwciP7cmi9/IpbMgFD8dy3it18m8PjJEcRTgTphQdzcshpta5Wh16LnMaWa\nMnLUU+5nIolA474Q86W1HVDWOqP3XFffDGsnwo550Oqe/M/Dr+ZmimlEf5j7Mvj4w8AvrSm355sl\nVRAa9YaGu/KfdNreZ413lKlsJbnAcrDuOzuptYKcTGhnJxwvX81AE4RSRc2+ldYvVh9/az5/aH3r\nBKycbMyP9yPrp7CndFPeYBi/JjVmpP88npL/sbTM86T2/j9CWnS2XmfVp5C0F/p+f+Ev4EZ9rV/E\nW3+Bdg+cmQqbl4jV/RIbbZ3odiVfjKVD4PEN4B98puvJ0y4lXpcL2t13dlnkQKjWEn55Amq0hc7P\ne31yAF2sr9Ds2bOHJk2aOB3GX+RdRFAVAekn4IeRUK4G3LfU6sKZdDuH42LZPmYQsn4Kb2XdRt/U\nFynVuAfj747iyX+8g/+DvxEQGk7IjGHWmcwZKbD4DQi/Dur9pUf3bDXbQ+lQQKDNveevVzoEGvQo\nmC/GwHKFlxwKSmhduGs6XP9CiUgOoC2IYi07OxtfXz2EXmXWs5C0H+6ZBWEN2XH9GGr/PIjAT6Ko\nLKlMLDeSRl2eIDqiMoF+eb5gK9aH4b/C9/fCrGesGTapCTD424t/mfn4WgOyyYetL0GlbNqC8JB3\n3nmHJk2a0KRJE9577z3A+kIfNmwYkZGRDBgwgJMnTwLw3HPPERERQWRk5OkzrhMSErj11ltp06YN\nbdq0Yfny5YB1DsWoUaPo0aMHd911F23btmXjxo2n37dz587ExMSQmprK8OHDadOmDS1atDi9RHha\nWhqDBw8mMjKSQYMGkZaWVpgfi7qQbbPhz//BdU+xPLMeg8auoNvUTP5thlNG0jna8V8MffxtbmxW\n9ezkcIp/EAz6xuofj99odR3VaJO/9+7wKPT8T8Hujyr2vPvn56znrLVSCtJVTaHX6xesEhMTwxdf\nfMHKlSsxxtC2bVs6derE1q1bGT9+PB06dGD48OF8/PHHDB8+nGnTprFlyxZEhOPHjwPw6KOP8vjj\nj3Pttdeyb98+brjhBjZv3nz69ZctW0apUqV49913mTJlCqNHjyYuLo6DBw/SqlUrXnjhBbp06cLn\nn3/O8ePHiYqKolu3bowdO5bSpUuzbt061q1bR8uWxXbRXO/z5ySySlXknu2dWDZnJVeVDeTFPo0Z\n1KYHLtdoQvyDLv4aLh/ri75hL+vfqlJXwLsThEOWLVvGzTffTFCQ9R/6lltuYenSpdSoUYMOHToA\ncMcdd/DBBx/w2GOPERgYyL333kufPn3o27cvAPPmzTvrEqUnTpwgOTkZgH79+lGqlDW98LbbbqN7\n9+6MHj2aKVOmnF4+fM6cOcyYMeP0SXnp6ens27ePJUuW8MgjjwAQGRlJZGRkIXwi6izJhyE57qyV\nRmOPJBGyeTYzMtuwJSuNl2+MYEhUzTwtBb9Le4/a1xVcvKrE8u4EcZFf+p5yvvWtzl1fRUTw9fVl\n1apVzJ8/n0mTJvHRRx+xYMECcnNzWbFixelEkNepxAPW4n+hoaGsW7eOyZMnM3bs2NMxfP/99zRs\n+NclAPRCQw6b9YzVnfTkFtJ9yzBuyS5iFk1ngk8qZZr2ZcnN11Pa37v/a6riQccgPKBjx478+OOP\nnDx5ktTUVKZNm8Z1113Hvn37WLFiBQDffvst1157LSkpKSQlJdG7d2/ee+891q5dC0CPHj346KOP\nTr/mqXJ3Bg8ezJtvvklSUtLppcNvuOEGPvzww9PJas2aNadjmzhxIgAbNmxg3bp1Bf8BqPPLSrcu\nmpOdxv4lX9Hng6W8M3cbd4VuxvgE0Kf/EE0OqsjQBOEBLVu25O677yYqKoq2bdty7733UqFCBRo3\nbsyECROIjIzk6NGj3H///SQnJ9O3b18iIyPp1KkT7777LgAffPAB0dHRREZGEhERwSeffHLe9xsw\nYACTJk3ittvOXD7ypZdeIisri8jISJo0acJLL70EwP33309KSgqRkZG8+eabREVFefbDUGfbvQSy\nUsl0lSJp+XjSMnOYMDyKrq41SO3rICDY6QiVOs2jy317WuvWrc25c/h1meqCoZ+jZyR/9wA+G3/g\nv1m38pLfN6TcPZ/g4HLwUeu/LrutlIfkd7lvbUEoVZBm/91ayM3ND6/ZGw6QufFnltCCawc+Br6l\nCN4w0VoQDvK3xLRShUg7O5UqKLk5EP0FZKVa291fBRFycw1vzdnKysWzuCEgidY976Ri8/qwu7+1\nxk9YA2sJ6vI1nY1fqXN4ZQuiOHebFQX6+V2mI9us5BDWGH77EJb+l/SsHB76djVjFu3k0erbMC4/\nKjbvY9VvOQwyk+FADDTs6WzsSrnhdQkiMDCQxMRE/ZK7TMYYEhMTCQx0s16+urAD9tVxB34JkYNg\nwatMfu9Jft1wkBd7N6JjzkprIDrQvgZ5zXZQsYF1v4EmCFX0eF0XU/Xq1YmNjSUhIcHpUIqtwMBA\nqlev7nQYxc+B1dZy2RUbcLDT22zfuI9hqV/Qp8YWKlZ5ChbstC6ic4oIXPektW6Su2WxlXKY1yUI\nPz8/ateu7XQYqiQ6EANVm7P/eDpDPosmKfsxfu7Ql1rRr8FE+4poDXuf/Zxmg91ff0GpIsDrEoRS\njshKh8MbOdb8b9w2dgUnM3OYOLIdtar3hDY3WtcR8CsNZas6HalS+aYJQqmCcHgD5Gbx+rpSZJpc\nJo1qR+MqZa3HKtSCO753Nj6lLoPXDVIr5YTkXasAWJ1dh4kj255JDkoVY5oglLpCJ9KzWLVsDgmm\nPK/f3ZNGV2lyUN5BE4RSVyAzO5dRX0VTK2MrUr0VrcJDnA5JqQKjCUKpy2SM4aUfN7BxVyz15CAV\nG7R3OiSlCpQmCKUu0+fL9zA5ej8vtsywCqrp1fmUd9EEodRlWLg1nn//sokbrq7MwKvircKqLZwN\nSqkCpglCqUu0L/Ekj3y7hoZXleWdgc1wHYiGkDpQWscflHfR8yCUugQZ2Tk8ODGGUfzAKP9tBLyz\n3VpwL3KQ06EpVeA0QSh1CV6buYWyh5bzsP8k8G0FzYdYC+417ud0aEoVOI8nCBHxAaKBA8aYviJS\nG5gEhACrgTuNMZkiEgB8BbQCEoFBxpg9no5PqfyatT6OL3/bzfKQ6eBXA+6ZBb4BToellMcUxhjE\no8DmPNtvAO8aY+oDx4ARdvkI4Jgxph7wrl1PqSLhwPE0nvl+HSMrbaXayU3Q6RlNDsrreTRBiEh1\noA/wmb0tQBdgql1lAtDfvn+TvY39eFe7vlKOys01PDP1T3Jzc3jabyqE1IVmtzsdllIe5+kWxHvA\nM0CuvR0KHDfGZNvbsUA1+341YD+A/XiSXV8pR339+16W70hkXIv9+CdugutfAB8dvlPez2MJQkT6\nAvHGmJi8xW6qmnw8lvd1R4lItIhE60WBlKftik9m0qx5/KvKcq7Z/QFUioCrb3E6LKUKhSd/BnUA\n+olIbyAQKIvVoigvIr52K6E6cNCuHwvUAGJFxBcoBxw990WNMeOAcQCtW7fW64oqj8k5eRzX2E7M\n8tljjZZVqA193gGXnj6kSgaP/Us3xjxvjKlujAkHBgMLjDFDgYXAALvaMGC6fX+GvY39+AKjF5ZW\nDvp9xqeE5+xhw9VPwSNr4dG1UEvXW1IlhxMdqc8Ck0TkX8AaYLxdPh74WkR2YLUc9DqMyjEHj6cR\ntHky+/3DufrWv2urQZVIhZIgjDGLgEX2/V1AlJs66cDAwohHqYsZM3UWr8p2jrd9GdHkoEoo/Zev\n1DlmbzxElT3TyBUfyrcd6nQ4SjlGE4RSubmwdRZkpZOSkc3oH9cx0G851OsKZSo7HZ1SjtHJ3Ept\nng7f3Q21OvBZ2D+pmxpDmH8iNNfWgyrZNEEotWM++JbC7F9F7z330KtCFcgpDw17OR2ZUo7SLiZV\nshkDOxdA/e68U/k1qpJIw9Q/oOkAXWtJlXiaIFTJdmQbnDjA9rJRfLi7Kr+0+Rzq94B2DzgdmVKO\n0y4mVbLtXADA6I2VqRVamv49O4Jvb4eDUqpo0BaEKtl2LuBE6VosOxLEC70bE+Dr43REShUZmiBU\nyZWdgdmzjFlpEUSFh9AjQqe0KpWXdjGpkmvf70jWSeZkRvBsr0bo5UeUOpsmCFVipW2Zi6/xIbjh\n9bSqVcHpcJQqcjRBqBLr+PrZ7DUNeKR3C6dDUapI0jEIVSLF7t9HlbRtJFW9jrphwU6Ho1SRpAlC\nlUhLZk0CoHWXWx2ORKmiSxOEKnF2JaRQYf9ckv0qElrvLyvPK6VsmiBUifPx3A10dP2Jb+PeeiEg\npS5AB6lVibIjPpnEjfMI8suApv2cDkepIk1/PqkS5b152+nluxrjFwS1OzodjlJFmiYIVWJsPZTM\nzPUH6O2/FqnfTVdrVeoiNEGoEuPDBdtp67+X4KxEaNTX6XCUKvI0QagSYWdCCr+sj+ORattAfKB+\nd6dDUqrI0wShSoSPF+4kwNdFm/QVEN4BSunSGkpdjCYI5fX2Hz3Jj2sP8FCkC9+j26BhH6dDUqpY\n0AShvN6YxTvxERiRNREQaKQXBFIqP/Q8COXV4pLSmBody8c15lNq23ToNhrK13Q6LKWKBW1BKK/2\n6ZLd9GAF3Q59Bs2GQIdHnQ5JqWJDWxDKax0/mcnqP5Yyxf8TqN4Wbnwf9KJASuWbtiCU1/p6xV4G\n5c7Cx88fBn2jJ8YpdYk8liBEJFBEVonInyKyUURG2+W1RWSliGwXkcki4m+XB9jbO+zHwz0Vm/J+\n6Vk5fLl8NzcErMenXhcIruR0SEoVO55sQWQAXYwxzYDmQE8RaQe8AbxrjKkPHANG2PVHAMeMMfWA\nd+16Sl2W76L3E5q2i5CcI1Cvm9PhKFUseSxBGEuKveln3wzQBZhql08A+tv3b7K3sR/vKnoVeXUZ\nsnNyGbd0F0NDtloFmiCUuiweHYMQER8RWQvEA3OBncBxY0y2XSUWqGbfrwbsB7AfTwJCPRmf8iLJ\nh+DwJgBmbjjE/qNp9A3aBJWuhrJVHQ5OqeLJownCGJNjjGkOVAeigMbuqtl/3bUWzLkFIjJKRKJF\nJDohIaHgglXF2+y/w/jumJR4xi3ZSURFFyFHYqBeV6cjU6rYKpRZTMaY48AioB1QXkROTa+tDhy0\n78cCNQDsx8sBR9281jhjTGtjTOuwsDBPh66Ki8MbIDOFwz+9woYDJ3i2QTySm6WL8il1BS4pQYiI\nS0TK5rNumIiUt++XAroBm4GFwAC72jBgun1/hr2N/fgCY8xfWhBK/UV2JiTuAL/SVNz6P5qXTqQD\na8EvCGq0czo6pYqtiyYIEfmfiJQVkSBgE7BVRJ7Ox2tXARaKyDrgD2CuMeZn4FngCRHZgTXGMN6u\nPx4ItcufAJ679N1RJVLiDsjN5kibp8gwvrwZMh3fXfOgTifw9Xc6OqWKrfycSR1hjDkhIkOBmVhf\n8DHAWxd6kjFmHdDCTfkurPGIc8vTgYH5CVqpsyRsBmBiQm18TF8eOvK9Vd7hMQeDUqr4y08Xk5+I\n+GFNR51ujMnCzeCxUo6J34wRF59u9uFI5N8gyB6b0gFqpa5IfloQY4E9wJ/AEhGpBZzwZFBKXZL4\nzRwLrEFKmi93droaIt6BXYugQrjTkSlVrF00QRhjPgA+yFO0V0Su91xISl2a3PgtrEm7im6NK1E3\nLBjC+kFEP6fDUqrYy88gdWURGS8is+ztCM7MNlLKWVnpcHQXG7KrMuLaOk5Ho5RXyc8YxJfAbODU\n6ajbAB39U0WCObIVF7mcLNeAdnVCnA5HKa+SnwRR0RgzBciF08tg5Hg0KqXyacu6PwBo1eYadOku\npQpWfhJEqoiEYs9csldkTfJoVErl0/aNf5CND53b6wlxShW0/MxiegLrLOe6IrIcCOPMmdBKOWbb\n4WRKHdtGUplahAYEOh2OUl4nP7OYVotIJ6Ah1oJ6W+1zIZRy1OfLdvOAK5bgGu2dDkUpr3TRBCEi\nd51T1FJEMMZ85aGYlLqoxJQMZq3ZyX98E3BVudrpcJTySvnpYmqT534g0BVYDWiCUI6ZuHIfNXNi\ncfkaqNTI6XCU8kr56WJ6OO+2iJQDvvZYREpdREZ2Dl//vpeRVZKsBeErRTgdklJe6XKuB3ESqF/Q\ngSiVXz//GUdCcga9Kh8HH3+oUNvpkJTySvkZg/iJM4vzuYAIYIong1LqfIwxjF+2myZhflRPWAxh\njcAnPz2lSqlLlZ//WW/nuZ8N7DXGxHooHqUuaOXuo2yKO8GCelOQ2O0wdKrTISnltfIzBrG4MAJR\nKj/GL9vNPaWWUif2R+j4DNTv5nRISnmt8yYIEUnG/XUfBDDGmHxdelSpgrLnSCoHtqzi48DPIbwT\ndNaLDirlSedNEMaYMoUZiFJwGTmpAAAWZUlEQVQX8+Vve/i37+e4SpWHWz8Dl4/TISnl1fI9i0lE\nKolIzVM3Twal1LmS0rJYGr2GFq7t+LS/H4IrOR2SUl4vP9eD6Cci24HdwGKsq8vN8nBcSp1l0qp9\ndMr53dporBcDUqow5KcF8SrQDthmjKmNdSb1co9GpVQe2Tm5TPhtD7cFrYbKTSC0rtMhKVUi5CdB\nZBljEgGXiLiMMQuB5h6OS6nTZm04RFbSIRpmboLGNzodjlIlRn7OgzguIsHAUmCiiMRjnQ+hlMed\nOjFuSNl1SKbR7iWlCtF5WxAi8pGIdABuwlpe4zHgV2AnoD/jVKGI2XuMtfuPM7jMWgitB5UaOx2S\nUiXGhVoQ27HOoq4CTAa+NcZMKJSolFr1KZQO5dPVNagZmE6Vo39Ah0dBLyuqVKG50HkQ7wPvi0gt\nYDDwhYgEAv8DJhtjthVSjKqkOXkUZj4NGG7LaYF/eBSyPwcitHtJqcJ00UFqY8xeY8wbxpgWwO3A\nLcBmj0emSq69vwGG6NAbae/axHX7x0K5mlBF50YoVZjycx6En4jcKCITsc5/2Abc6vHIVMm1ZynG\ntxTD4wfzXoMJENEfrntCu5eUKmQXWoupOzAE6AOsAiYBo4wxqfl5YRGpgXXVuauAXGCcMeZ9EQnB\nGtMIxzrp7jZjzDEREeB9oDfWoPjdxpjVl7lfqjjbvYT9ZSI5kSLc2vVauKqX0xEpVSJdqAXxArAC\naGyMudEYMzG/ycGWDTxpjGmMdaLdgyISATwHzDfG1Afm29sAvbAuRFQfGAWMubRdUV4h9QjEb+Kn\npHp0bBBGw6t0STClnHKhQerrr+SFjTFxQJx9P1lENgPVsKbNdrarTQAWAc/a5V8ZYwzwu4iUF5Eq\n9uuokmLPUgDmpjXk6Y51HA5GqZLtci45eslEJBxoAawEKp/60rf/nlp1rRqwP8/TYu2yc19rlIhE\ni0h0QkKCJ8NWDjC7l5JKKXKrNOOauqFOh6NUiebxBGGfhf098Jgx5sSFqrop+8v1KIwx44wxrY0x\nrcPCwgoqTFVEpG5dwMqchozs1ADRQWmlHOXRBCEifljJYaIx5ge7+LCIVLEfrwLE2+WxQI08T68O\nHPRkfKpoMScOEpy8m82BzenV5Cqnw1GqxPNYgrBnJY0HNhtj3snz0AxgmH1/GDA9T/ldYmkHJOn4\nQ8myK3o2AOGte+LrUyi9n0qpC8jPYn2XqwNwJ7BeRNbaZS8ArwNTRGQEsA8YaD82E2uK6w6saa73\neDA2VQQdWD2bMILo0qmr06EopfBggjDGLMP9uAJY15Q4t74BHvRUPKpo2xJ3nPAT0cRXbEO9QH+n\nw1FKUUizmJS6mN1T/0FNVwJV2g9yOhSllE0ThHLc4VXf0ytxAusq9iGo9RCnw1FK2TRBKGfFb6H8\nrw+yztSlytAxut6SUkWIJgjlnKx0sr+9neQcP+Y0eZuwCuWcjkgplYcmCOWctd/ge2wnT2c/wJDu\n7Z2ORil1Dk9Oc1Xq/HKyyVn2PutMfSo260W18qWcjkgpdQ5tQShnbPwBn6R9jMnuxwNd6jsdjVLK\nDW1BqMKXm0v24v+yy1SnbOSN1K4Y5HRESik3tAWhCt/22fgmbuGTnH483K2B09Eopc5DWxCqcBlD\n5qK3OGzC8IscQK1QbT0oVVRpC0IVroOr8Y+L4bOcvjzYtbHT0SilLkAThCpUKTGTyTQ+0HQgNUNL\nOx2OUuoCNEGowpObS866H1himnNv9xZOR6OUughNEKrQxK5fSLnsBJLq9qNGiLYelCrqNEGoQrN9\n/gTSjD/X97vL6VCUUvmgCUIVitV7EmiStJDYsI6EVAhxOhylVD5oglAeZ4xh5owphMkJanS8w+lw\nlFL5pAlCedzCrfHUi59Dlk9pAhv3dDocpVQ+aYJQHpWZncsbP/1Jb99ofCL6gp8uyqdUcaEJQnlO\nbi7Lpo/jk+SHKEsKruZ6tTilihNdakN5RlIs2d8OpcuhtewPCMcMmIzU7eJ0VEqpS6AtCFXwcnNh\n2n3kxG/lqez7SR++GGmoYw9KFTeaIFTBWzUW9izl5Yw7CGpzB/WrlHc6IqXUZdAuJlWwErZh5v2T\n1QFRzJbuLNDlvJUqtrQFoQpOTjZM+xuZEsh9SXfzfJ8IKgT5Ox2VUuoyaYJQV84Y2DQDxl4HB1fz\nYvZwaofXYUDL6k5HppS6AtrFpK7MkR0w9R44tA5C6/NV9dFM29WAmTc3weUSp6NTSl0BbUGoKzPv\nZTi2F/p/woqev/CPHfUZ2bEODSqXcToypdQV8liCEJHPRSReRDbkKQsRkbkist3+W8EuFxH5QER2\niMg6EWnpqbhUATq6C7b8AlEjSW08kGembaRWaGke6VLf6ciUUgXAky2IL4FzJ78/B8w3xtQH5tvb\nAL2A+vZtFDDGg3GpgrJyLLh8IWokr8/aQuyxNN4a0IxS/j5OR6aUKgAeSxDGmCXA0XOKbwIm2Pcn\nAP3zlH9lLL8D5UWkiqdiUwUg7Tis/hqaDmD5YV++/n0vwzvUJqq2LuWtlLco7DGIysaYOAD7byW7\nvBqwP0+9WLvsL0RklIhEi0h0QkKCR4NVF7D6K8hKJbXlKJ6Zuo46FYN4+oaGTkellCpARWWQ2t10\nF+OuojFmnDGmtTGmdVhYmIfDUm7lZFndS+HX8dLvLuKS0nhrYDMC/bRrSSlvUtgJ4vCpriP7b7xd\nHgvUyFOvOnCwkGNT+RXzJZyIZXnYIH5Yc4CHu9SnVa0KTkellCpghZ0gZgDD7PvDgOl5yu+yZzO1\nA5JOdUWpIiQ7A2Y+AzOfIq1qO0b+Hkrb2iE80lVnLSnljTx2opyIfAt0BiqKSCzwMvA6MEVERgD7\ngIF29ZlAb2AHcBK4x1Nxqct0fD9MvgPi1pLd9gEGbulGoH8u7w9ugY+eEKeUV/JYgjDGnO/qMF3d\n1DXAg56KRRWA2S9A4g7MoIn8fWNNNhzezxd3t+GqcoFOR6aU8pCiMkitirLMk7B9LjQbwlfHmjA5\nej8PXl+X6xtVuvhzlVLFliYIdXE750N2GhvKdeSVnzfRrXElnuyuU1qV8na6WJ+6uM0/kRNYgWHz\n/ahTsTTvDmquC/EpVQJoC0JdWHYmuVtnMSe7JVnGxad3taZMoJ/TUSmlCoEmCHVB6TsW4co4wbSM\nlnw2rA3hFYOcDkkpVUg0QajzysrJZcXPX5JiAhk48E5dZ0mpEkYThHIrOyeXJyevpknyUhKrdqZ7\nZC2nQ1JKFTIdpFZny84kO+UIj/9ykEMbFhMWcAI6DHI6KqWUAzRBqDOO7cF81R/fY7t52/hCUCkw\nAVC/u9ORKaUcoAlCWeK3YL7qz8mTKbyXdTt96/nTrEwyVG0JAXr5UKVKIk0QJV1WGuxbQe7UERzP\nMAxJf5GhN/akWftwpyNTSjlME0RJteRtWDcFEreDyeWwVOLOrBd46vZe9GxyldPRKaWKAE0QJdHh\nTbDgVajRlj0R9/PhpkBW0ZT3RnbW6zoopU7TBFESrRyD8Q1kbNV/8caSBBpdVZZv7mhJrVA9CU4p\ndYYmiJIm9Qjmz8ksLNWN1xcncEvLavy7f1NK+evlQpVSZ9MEUcJs+eV9GuVk8Pbx6/lX/yYMbVsT\nEV14Tyn1V5ogSoj45HRem/Enz2/7itX+rfjwwSHUDQt2OiylVBGmCcLLZefk8vXve3lnzjZ65iyi\nku9xQgc+j48mB6XURehaTMVJ2nH44zPIyT67PDcXts+zzmmwGWNYtDWevh8uY/RPm7i2mot/V5oP\nFRviU79bIQeulCqOtAVRnCx41UoQPv7Q8q4z5esmwY/3Q/h1cPtkYuIyefPXLazcfZQaIaWY1COL\ntmueQlITYOCXoGMOSql80ARRXBzdDTFfAgKL34LIweDrD9kZsPA1TNlqsHc5m97uyR0nHiMouBxv\n3VCRWzJ/xmfJhxBaF4bMg6rNnd4TpVQxoQmiuFj4H3D5wY3vw/QHYc1X0OZeMlZ9QUDSPl4MGk1y\nRjzvmjEsqvQeFUNC8Fm8GDDQ4g7o+QYE6LiDUir/NEEUB4c2wPrvoMOj0HwoZvXXZC18kzf3NOZv\nG15jl2lEjG8LRtxSB+PXgso//g1c1aDTMxA5yGo9KKXUJdIEURwseBUTWJbNdYczc8429sb34cOM\nlxiw4T7C5DiJfT5lVpuO9vkMNaB+VwgsDy6dg6CUunyaIIqK9CQ4thcqNzn9xZ6cmsqBeZ/QaNuv\nfOIzlDfGbcAl0KFeW+LT2tLoyEqo151GUT3Ofq3SemlQpdSV0wRxJXJzIXo81OsKIXXy95yEbXDy\nCNRsf2Y20d4VMHU4JB8krXRV1pXrwubkUvRI/oFGksga04BNNYfw1tXhdGlUidDgADjwGky+E7q9\n7Ln9U0qVaGKMcTqGy9a6dWsTHR1dcC+Ykw3JB6F8zYvXNQZmPQOrxkH1NjBi7sWnj6YkwJhrIDWe\nrEqR7Kw/gqRDO2m98yPipDKfZPaki2sN17nW4yc5xJZpRnKbx6jTvh8BfprLlVIFQ0RijDGtL1ZP\nv3UAjuyAtd/A2m8h5RCMXAjVWl74OYvftJJDtdYQ+wds+B6aDjiryon0LPYfPcn+o2nsPZJCp5gH\nqZN6jLFyO70PLaRR/KMAzHddw/Saz9GgZlVK1QohOzQHv4wjVK8c4ak9VkqpiypSCUJEegLvAz7A\nZ8aY1z3+pn9Ogml/A/GBet1gZ6L1ZZ83QRgDW2dBTgY5fsFkxK6l9JL/cKTuLcQ0e4U2c2/F96e/\n8+aWWsSmGOKOpxOXlMaJ9DNnPN/pM4dGfr8zvsx9xNYaypLKT5KZHU218v50bX4zXf/S+qjk8V1X\nSqkLKTJdTCLiA2wDugOxwB/AEGPMpvM953K7mA4eT2P3kVSyT8TTflZPTgTXZn7Tt0h0hdBz7SOE\nnNzJv+pNJiUjhxPpWdQ7sYpXkl866zXm5bTgvqzHycaX9q6NfOv/bz50DWVOhdupVtaXjqymSkAG\npcNqERbsT+25I5Dw62Dod3oms1LKUcWxiykK2GGM2QUgIpOAm4DzJojL9dOfB3lt1hb+6zcGXCkM\nPnQ7O+ISgAT2+DbhTd/lHNuxiv2lGlI20I/+Ob+S4lOOr+q9TwX/HMr6Q3bV1nxRNoiQIH/CgruS\n+/MaHt4zg4eb1IDVX1ljGQBb7DcNCoP+H2tyUEoVG0UpQVQD9ufZjgXanltJREYBowBq1szHYLIb\nfZtVpaPfFhrPWcrRlg/zxbV3EhTgS2l/HwKy2sPb4xnf5gB0HwVJB+C93+Gah3ig+83nf9Eer8LH\nbWHx61C3C/R9F8IaWM8/cRBqREGwdhsppYqPopQg3P20/kv/lzFmHDAOrC6my3mjasEuqsX8AyqE\nE9Lr74T4lTrzoF8I1O4Im2dAt3/Cmq/B5ECruy/8ohXrwfDZEFAGwhqeKc/v9FellCpiitKptrFA\njTzb1YGDHnmn3z6AxB3Q57+QNzmcEnETHN0FcWshZgLUzed5DtVbn50clFKqGCtKCeIPoL6I1BYR\nf2AwMMMj79R8KPR83Zq15E6jviAu+OlRayyh9XCPhKGUUkVZkUkQxphs4CFgNrAZmGKM2eiRNytb\nFdrdf/7HgypC+LUQ9yeUqQoNenokDKWUKsqK0hgExpiZwEyn4wCsbqbdS6DVMPApUh+TUkoVCv3m\nO5+mt1lnWEeNcjoSpZRyhCaI8wksC708fyK3UkoVVUVmDEIppVTRoglCKaWUW5oglFJKuaUJQiml\nlFuaIJRSSrmlCUIppZRbmiCUUkq5pQlCKaWUW0XminKXQ0QSgL2X+fSKwJECDKc40H0uGXSfS4Yr\n2edaxpiwi1Uq1gniSohIdH4uuedNdJ9LBt3nkqEw9lm7mJRSSrmlCUIppZRbJTlBjHM6AAfoPpcM\nus8lg8f3ucSOQSillLqwktyCUEopdQElMkGISE8R2SoiO0TkOafjKQgiUkNEForIZhHZKCKP2uUh\nIjJXRLbbfyvY5SIiH9ifwToRaensHlw+EfERkTUi8rO9XVtEVtr7PNm+xjkiEmBv77AfD3cy7ssl\nIuVFZKqIbLGPd3tvP84i8rj973qDiHwrIoHedpxF5HMRiReRDXnKLvm4isgwu/52ERl2JTGVuAQh\nIj7A/wG9gAhgiIhEOBtVgcgGnjTGNAbaAQ/a+/UcMN8YUx+Yb2+Dtf/17dsoYEzhh1xgHsW6jvkp\nbwDv2vt8DBhhl48Ajhlj6gHv2vWKo/eBX40xjYBmWPvutcdZRKoBjwCtjTFNAB9gMN53nL8Eep5T\ndknHVURCgJeBtkAU8PKppHJZjDEl6ga0B2bn2X4eeN7puDywn9OB7sBWoIpdVgXYat8fCwzJU/90\nveJ0A6rb/3G6AD8DgnXykO+5xxuYDbS37/va9cTpfbjE/S0L7D43bm8+zkA1YD8QYh+3n4EbvPE4\nA+HAhss9rsAQYGye8rPqXeqtxLUgOPOP7ZRYu8xr2E3qFsBKoLIxJg7A/lvJruYtn8N7wDNArr0d\nChw3xmTb23n36/Q+248n2fWLkzpAAvCF3a32mYgE4cXH2RhzAHgb2AfEYR23GLz7OJ9yqce1QI93\nSUwQ4qbMa6ZyiUgw8D3wmDHmxIWquikrVp+DiPQF4o0xMXmL3VQ1+XisuPAFWgJjjDEtgFTOdDu4\nU+z32e4iuQmoDVQFgrC6WM7lTcf5Ys63jwW67yUxQcQCNfJsVwcOOhRLgRIRP6zkMNEY84NdfFhE\nqtiPVwHi7XJv+Bw6AP1EZA8wCaub6T2gvIj42nXy7tfpfbYfLwccLcyAC0AsEGuMWWlvT8VKGN58\nnLsBu40xCcaYLOAH4Bq8+zifcqnHtUCPd0lMEH8A9e0ZEP5Yg10zHI7piomIAOOBzcaYd/I8NAM4\nNZNhGNbYxKnyu+zZEO2ApFNN2eLCGPO8Maa6MSYc6zguMMYMBRYCA+xq5+7zqc9igF2/WP2yNMYc\nAvaLSEO7qCuwCS8+zlhdS+1EpLT97/zUPnvtcc7jUo/rbKCHiFSwW1497LLL4/SgjEMDQb2BbcBO\n4O9Ox1NA+3QtVlNyHbDWvvXG6nudD2y3/4bY9QVrNtdOYD3WDBHH9+MK9r8z8LN9vw6wCtgBfAcE\n2OWB9vYO+/E6Tsd9mfvaHIi2j/WPQAVvP87AaGALsAH4GgjwtuMMfIs1xpKF1RIYcTnHFRhu7/sO\n4J4riUnPpFZKKeVWSexiUkoplQ+aIJRSSrmlCUIppZRbmiCUUkq5pQlCKaWUW74Xr6KUEpFT0w0B\nrgJysJa8ADhpjLnGkcCU8iCd5qrUJRKRfwIpxpi3nY5FKU/SLialrpCIpNh/O4vIYhGZIiLbROR1\nERkqIqtEZL2I1LXrhYnI9yLyh33r4OweKOWeJgilClYzrOtTNAXuBBoYY6KAz4CH7TrvY13HoA1w\nq/2YUkWOjkEoVbD+MPZaRyKyE5hjl68HrrfvdwMirGWFACgrImWMMcmFGqlSF6EJQqmClZHnfm6e\n7VzO/H9zYV3QJq0wA1PqUmkXk1KFbw7w0KkNEWnuYCxKnZcmCKUK3yNAa/ti85uA+5wOSCl3dJqr\nUkopt7QFoZRSyi1NEEoppdzSBKGUUsotTRBKKaXc0gShlFLKLU0QSiml3NIEoZRSyi1NEEoppdz6\nf9ajMeVaQKDZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x26de8c2abe0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pints\n",
    "import pints.toy as toy\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Load a forward model\n",
    "model = toy.LogisticModel()\n",
    "\n",
    "# Create some toy data\n",
    "r = 0.015\n",
    "k = 500\n",
    "real_parameters = [r, k]\n",
    "times = np.linspace(0, 1000, 100)\n",
    "signal_values = model.simulate(real_parameters, times)\n",
    "\n",
    "# Add independent Gaussian noise\n",
    "sigma = 10\n",
    "observed_values = signal_values + pints.noise.independent(sigma, signal_values.shape)\n",
    "\n",
    "# Plot\n",
    "plt.plot(times,signal_values,label = 'signal')\n",
    "plt.plot(times,observed_values,label = 'observed')\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Values')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the nested sampler that will be used to sample from the posterior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create an object with links to the model and time series\n",
    "problem = pints.SingleOutputProblem(model, times, observed_values)\n",
    "\n",
    "# Create a log-likelihood function (adds an extra parameter!)\n",
    "log_likelihood = pints.GaussianLogLikelihood(problem)\n",
    "\n",
    "# Create a uniform prior over both the parameters and the new noise variable\n",
    "log_prior = pints.UniformLogPrior(\n",
    "    [0.01, 400, sigma * 0.5],\n",
    "    [0.02, 600, sigma * 1.5])\n",
    "\n",
    "# Create a nested ellipsoidal rejectection sampler\n",
    "sampler = pints.NestedController(log_likelihood, log_prior, method=pints.NestedEllipsoidSampler)\n",
    "\n",
    "# Set number of iterations\n",
    "sampler.set_iterations(8000)\n",
    "\n",
    "# Set the number of posterior samples to generate\n",
    "sampler.set_n_posterior_samples(1600)\n",
    "\n",
    "# Do proposals in parallel\n",
    "sampler.set_parallel(True)\n",
    "\n",
    "# Use dynamic enlargement factor\n",
    "sampler._sampler.set_dynamic_enlargement_factor(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the sampler!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running Nested ellipsoidal sampler\n",
      "Number of active points: 400\n",
      "Total number of iterations: 8000\n",
      "Total number of posterior samples: 1600\n",
      "Iter. Eval. Time m:s Delta_log(z) Acceptance rate\n",
      "0     1       0:05.6 -inf          1             \n",
      "0     2       0:05.6 -inf          1             \n",
      "0     21      0:05.6 -inf          1             \n",
      "0     41      0:05.7 -inf          1             \n",
      "0     61      0:05.7 -inf          1             \n",
      "0     81      0:05.8 -inf          1             \n",
      "0     101     0:05.8 -inf          1             \n",
      "0     121     0:05.8 -inf          1             \n",
      "0     141     0:05.9 -inf          1             \n",
      "0     161     0:05.9 -inf          1             \n",
      "0     181     0:06.0 -inf          1             \n",
      "0     201     0:06.0 -inf          1             \n",
      "0     221     0:06.1 -inf          1             \n",
      "0     241     0:06.1 -inf          1             \n",
      "0     261     0:06.2 -inf          1             \n",
      "0     281     0:06.2 -inf          1             \n",
      "0     301     0:06.3 -inf          1             \n",
      "0     321     0:06.3 -inf          1             \n",
      "0     341     0:06.4 -inf          1             \n",
      "0     361     0:06.5 -inf          1             \n",
      "0     381     0:06.5 -inf          1             \n",
      "1     404     0:06.6 -inf          0.25          \n",
      "20    424     0:06.6 -11961.96419  0.833333333   \n",
      "40    444     0:06.6 -8352.604183  0.909090909   \n",
      "60    464     0:06.7 -6359.778919  0.9375        \n",
      "80    488     0:06.7 -5275.768444  0.909090909   \n",
      "100   520     0:06.8 -4590.47119   0.833333333   \n",
      "120   544     0:06.8 -4326.757781  0.833333333   \n",
      "140   576     0:06.9 -3719.789341  0.795454545   \n",
      "160   600     0:06.9 -3392.855978  0.8           \n",
      "180   636     0:07.0 -3108.993827  0.762711864   \n",
      "200   660     0:07.1 -2754.993979  0.769230769   \n",
      "220   688     0:07.1 -2529.775272  0.763888889   \n",
      "240   720     0:07.2 -2351.601106  0.75          \n",
      "260   756     0:07.2 -2175.081064  0.730337079   \n",
      "280   792     0:07.3 -2044.625718  0.714285714   \n",
      "300   820     0:07.3 -1918.872831  0.714285714   \n",
      "320   856     0:07.4 -1805.151712  0.701754386   \n",
      "340   888     0:07.4 -1689.317568  0.696721311   \n",
      "360   936     0:07.5 -1586.3349    0.671641791   \n",
      "380   980     0:08.3 -1489.927673  0.655172414   \n",
      "400   1032    0:08.4 -1401.254961  0.632911392   \n",
      "420   1064    0:08.4 -1310.834638  0.63253012    \n",
      "440   1096    0:08.5 -1216.232591  0.632183908046\n",
      "460   1132    0:08.6 -1149.102977  0.628415301   \n",
      "480   1172    0:08.7 -1086.938189  0.621761658   \n",
      "500   1216    0:09.9 -1032.899956  0.612745098   \n",
      "520   1252    0:09.9 -971.9526094  0.610328638   \n",
      "540   1284    0:10.0 -919.9015005  0.610859729   \n",
      "560   1336    0:10.1 -865.2953872  0.598290598   \n",
      "580   1380    0:10.2 -818.8343722  0.591836735   \n",
      "600   1428    0:10.3 -770.6446078  0.583657588   \n",
      "620   1464    0:10.3 -743.0924604  0.582706767   \n",
      "640   1508    0:10.4 -699.0575667  0.577617329   \n",
      "660   1560    0:10.5 -665.7779445  0.568965517   \n",
      "680   1616    0:10.7 -642.617847   0.559210526   \n",
      "700   1664    0:10.8 -610.8239256  0.553797468   \n",
      "720   1704    0:10.9 -582.4655629  0.552147239   \n",
      "740   1728    0:10.9 -553.9587075  0.557228916   \n",
      "760   1768    0:11.0 -529.239175   0.555555556   \n",
      "780   1816    0:11.1 -503.4733143  0.550847458   \n",
      "800   1884    0:11.2 -477.9142571  0.539083558   \n",
      "820   1912    0:11.3 -443.8719624  0.542328042328\n",
      "840   1964    0:11.4 -420.2479985  0.537084398977\n",
      "860   2008    0:17.5 -396.8890682  0.534825871   \n",
      "880   2048    0:17.6 -377.7338834  0.533980583   \n",
      "900   2092    0:17.6 -352.209204   0.531914893617\n",
      "920   2140    0:17.7 -324.9226736  0.528735632   \n",
      "940   2200    0:17.8 -307.4450264  0.522222222   \n",
      "960   2248    0:17.8 -290.8736334  0.519480519   \n",
      "980   2292    0:17.9 -277.4842442  0.517970402   \n",
      "1000  2360    0:18.0 -259.9878235  0.510204082   \n",
      "1020  2408    0:18.1 -247.6515773  0.50796812749 \n",
      "1040  2460    0:18.1 -237.333319   0.504854368932\n",
      "1060  2504    0:19.0 -218.9319717  0.503802281   \n",
      "1080  2552    0:19.1 -207.3628985  0.501858736   \n",
      "1100  2596    0:19.1 -195.246171   0.500910747   \n",
      "1120  2640    0:19.2 -186.031126   0.5           \n",
      "1140  2684    0:19.2 -174.9244619  0.499124343   \n",
      "1160  2728    0:19.3 -165.854887   0.498281787   \n",
      "1180  2768    0:19.3 -158.9586757  0.498310811   \n",
      "1200  2816    0:19.4 -153.9111851  0.496688742   \n",
      "1220  2864    0:19.5 -147.3558753  0.49512987    \n",
      "1240  2916    0:19.5 -145.8282542  0.492845787   \n",
      "1260  2968    0:19.6 -140.7466328  0.490654206   \n",
      "1280  3012    0:20.5 -137.1815534  0.490045941807\n",
      "1300  3056    0:20.5 -133.8850999  0.489457831   \n",
      "1320  3104    0:20.6 -129.4724271  0.48816568    \n",
      "1340  3160    0:20.6 -125.7228883  0.485507246   \n",
      "1360  3196    0:20.7 -120.279964   0.486409156   \n",
      "1380  3268    0:20.8 -116.3331368  0.481171548   \n",
      "1400  3340    0:20.8 -113.3900068  0.476190476   \n",
      "1420  3424    0:20.9 -110.5637194  0.46957672    \n",
      "1440  3540    0:21.0 -108.0741926  0.458598726   \n",
      "1460  3612    0:21.1 -105.543546   0.454545455   \n",
      "1480  3664    0:21.2 -103.2284003  0.453431372549\n",
      "1500  3732    0:21.2 -101.2561546  0.450180072   \n",
      "1520  3828    0:21.4 -99.38890983  0.443407235   \n",
      "1540  3940    0:21.9 -97.54519601  0.435028249   \n",
      "1560  4056    0:23.8 -95.80301875  0.426695842   \n",
      "1580  4164    0:24.0 -93.82367642  0.419766206   \n",
      "1600  4312    0:24.2 -91.55523354  0.408997955   \n",
      "1620  4336    0:24.2 -89.43104957  0.411585366   \n",
      "1640  4368    0:24.3 -87.64257243  0.413306452   \n",
      "1660  4396    0:24.3 -86.06338007  0.415415415   \n",
      "1680  4428    0:24.4 -84.44832862  0.417080437   \n",
      "1700  4464    0:24.4 -83.18414238  0.418307087   \n",
      "1720  4492    0:24.5 -81.95574752  0.420332356   \n",
      "1740  4516    0:26.0 -80.79838234  0.422740525   \n",
      "1760  4548    0:26.0 -79.43406283  0.424300868   \n",
      "1780  4580    0:26.1 -78.0530443   0.425837321   \n",
      "1800  4612    0:26.1 -76.75921356  0.427350427   \n",
      "1820  4648    0:26.2 -75.52914362  0.428436911   \n",
      "1840  4684    0:26.2 -74.49560163  0.429505135   \n",
      "1860  4720    0:26.3 -73.44291573  0.430555556   \n",
      "1880  4752    0:26.3 -72.35234981  0.431985294   \n",
      "1900  4784    0:26.4 -71.32262181  0.433394161   \n",
      "1920  4824    0:26.4 -70.0908834   0.433996383   \n",
      "1940  4860    0:26.5 -68.69185073  0.434977578   \n",
      "1960  4912    0:26.5 -67.42729521  0.434397163   \n",
      "1980  4956    0:26.6 -66.25324711  0.434591747   \n",
      "2000  4996    0:26.6 -65.14478357  0.43516101    \n",
      "2020  5052    0:27.4 -63.97010481  0.43422184    \n",
      "2040  5096    0:27.4 -62.82926098  0.434412266   \n",
      "2060  5140    0:27.5 -61.77683486  0.434599156   \n",
      "2080  5212    0:27.6 -60.48984022  0.432252702   \n",
      "2100  5280    0:27.7 -59.24114129  0.430327869   \n",
      "2120  5308    0:27.7 -58.11618716  0.43194784    \n",
      "2140  5336    0:27.8 -57.16139475  0.433549433   \n",
      "2160  5360    0:27.9 -56.13369893  0.435483871   \n",
      "2180  5404    0:27.9 -54.82890499  0.435651479   \n",
      "2200  5436    0:28.0 -53.53591305  0.436854647   \n",
      "2220  5456    0:28.0 -52.47235975  0.439082278481\n",
      "2240  5484    0:28.0 -51.42696157  0.440597954   \n",
      "2260  5508    0:28.1 -50.50509087  0.442443226   \n",
      "2280  5540    0:28.1 -49.45177092  0.443579766537\n",
      "2300  5572    0:28.2 -48.52654179  0.444702243   \n",
      "2320  5604    0:28.2 -47.47265284  0.445810914681\n",
      "2340  5632    0:28.3 -46.63052003  0.447247706422\n",
      "2360  5672    0:28.3 -45.65977965  0.447647951   \n",
      "2380  5716    0:28.4 -44.75842374  0.447705041   \n",
      "2400  5752    0:28.5 -43.96567828  0.448430493   \n",
      "2420  5776    0:28.5 -43.30026266  0.45014881    \n",
      "2440  5796    0:28.6 -42.5866327   0.452186805   \n",
      "2460  5820    0:28.6 -41.81378689  0.453874539   \n",
      "2480  5852    0:28.6 -40.91586948  0.454878943507\n",
      "2500  5884    0:28.7 -40.09987659  0.45587162655 \n",
      "2520  5916    0:29.2 -39.31693858  0.456852792   \n",
      "2540  5944    0:29.3 -38.45608851  0.458152958153\n",
      "2560  5972    0:29.3 -37.58393358  0.45944005743 \n",
      "2580  6016    0:29.4 -36.68273763  0.459401709   \n",
      "2600  6056    0:29.4 -35.91479435  0.459688826   \n",
      "2620  6076    0:29.5 -35.20875886  0.461592670895\n",
      "2640  6104    0:29.6 -34.42035725  0.4628331     \n",
      "2660  6132    0:29.9 -33.52393033  0.46406141    \n",
      "2680  6156    0:30.0 -32.69261296  0.465601112   \n",
      "2700  6180    0:30.0 -31.95660541  0.467128028   \n",
      "2720  6216    0:30.0 -31.49113766  0.467675378   \n",
      "2740  6244    0:30.1 -30.78795045  0.468856947   \n",
      "2760  6288    0:30.2 -30.07688502  0.46875       \n",
      "2780  6332    0:30.2 -29.42305105  0.468644639   \n",
      "2800  6368    0:30.3 -28.84538666  0.469168901   \n",
      "2820  6408    0:30.3 -28.24930599  0.469374168   \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2840  6448    0:30.4 -27.47903096  0.46957672    \n",
      "2860  6496    0:30.5 -26.70709881  0.469160105   \n",
      "2880  6532    0:31.0 -25.93762189  0.469667319   \n",
      "2900  6576    0:31.0 -25.27702144  0.469559585   \n",
      "2920  6596    0:31.1 -24.61540232  0.471271788   \n",
      "2940  6620    0:31.1 -24.0467441   0.47266881    \n",
      "2960  6644    0:31.2 -23.48313358  0.474055093   \n",
      "2980  6668    0:31.2 -22.93188037  0.475430759   \n",
      "3000  6692    0:31.3 -22.39463485  0.476795931   \n",
      "3020  6720    0:31.3 -21.89199034  0.477848101   \n",
      "3040  6748    0:31.3 -21.38515542  0.478890989288\n",
      "3060  6780    0:31.4 -20.8674129   0.479623824   \n",
      "3080  6808    0:31.4 -20.37271826  0.480649189   \n",
      "3100  6836    0:31.5 -19.84836779  0.481665631   \n",
      "3120  6876    0:31.6 -19.29009278  0.481778876   \n",
      "3140  6904    0:31.6 -18.73268854  0.482779828   \n",
      "3160  6948    0:31.7 -18.1433521   0.482590104   \n",
      "3180  6988    0:31.7 -17.54468311  0.482695811   \n",
      "3200  7036    0:32.6 -16.94597411  0.482218204   \n",
      "3220  7060    0:32.6 -16.43274182  0.483483483   \n",
      "3240  7084    0:32.7 -15.98274995  0.484739677   \n",
      "3260  7108    0:32.7 -15.56086098  0.485986881   \n",
      "3280  7140    0:32.8 -15.26496322  0.486646884273\n",
      "3300  7180    0:32.8 -14.90929632  0.486725664   \n",
      "3320  7220    0:32.9 -14.56597025  0.486803519   \n",
      "3340  7252    0:32.9 -14.2088416   0.48744892    \n",
      "3360  7292    0:33.0 -13.83524883  0.487521764   \n",
      "3380  7332    0:33.1 -13.45993124  0.487593768   \n",
      "3400  7372    0:33.1 -13.08426307  0.487664945   \n",
      "3420  7400    0:33.2 -12.67106451  0.488571429   \n",
      "3440  7428    0:33.2 -12.27595072  0.489470689   \n",
      "3460  7452    0:33.3 -11.91206828  0.490640953   \n",
      "3480  7484    0:33.3 -11.55856633  0.491247883   \n",
      "3500  7516    0:33.4 -11.19961339  0.491849354   \n",
      "3520  7548    0:33.4 -10.83760494  0.492445439   \n",
      "3540  7584    0:34.0 -10.4824231   0.492761693   \n",
      "3560  7624    0:34.0 -10.36336527  0.492801772   \n",
      "3580  7656    0:34.1 -10.0575758   0.493384785   \n",
      "3600  7692    0:34.1 -9.776437933  0.493691717   \n",
      "3620  7716    0:34.2 -9.517409977  0.494805904866\n",
      "3640  7740    0:34.2 -9.276152848  0.495912807   \n",
      "3660  7768    0:34.3 -9.039111206  0.496742671   \n",
      "3680  7800    0:34.4 -8.798967429  0.497297297   \n",
      "3700  7836    0:34.4 -8.552023806  0.497579344   \n",
      "3720  7864    0:34.5 -8.29645501   0.498392283   \n",
      "3740  7892    0:34.5 -8.04755083   0.499199146   \n",
      "3760  7920    0:34.6 -7.807561001  0.5           \n",
      "3780  7964    0:34.6 -7.56765474   0.49973559    \n",
      "3800  8004    0:35.4 -7.335283629  0.499736981   \n",
      "3820  8028    0:35.5 -7.11477005   0.500786576   \n",
      "3840  8056    0:35.5 -6.900830111  0.501567398   \n",
      "3860  8080    0:35.5 -6.689921705  0.502604167   \n",
      "3880  8108    0:35.6 -6.487720734  0.503373119   \n",
      "3900  8144    0:35.7 -6.29099825   0.503615702   \n",
      "3920  8172    0:35.7 -6.094812697  0.504374678   \n",
      "3940  8208    0:35.8 -5.905216015  0.504610656   \n",
      "3960  8236    0:35.8 -5.72046549   0.505359877   \n",
      "3980  8260    0:35.9 -5.538291728  0.506361323   \n",
      "4000  8292    0:35.9 -5.361538973  0.506842372   \n",
      "4020  8324    0:36.0 -5.183169443  0.507319536   \n",
      "4040  8356    0:36.0 -5.005499661  0.507792860734\n",
      "4060  8384    0:36.1 -4.8302518    0.508517034   \n",
      "4080  8432    0:36.1 -4.661826     0.50796812749 \n",
      "4100  8468    0:36.2 -4.499569064  0.508180466   \n",
      "4120  8512    0:37.2 -4.341571823  0.507889546   \n",
      "4140  8552    0:37.2 -4.18961608   0.507850834   \n",
      "4160  8600    0:37.3 -4.041397337  0.507317073   \n",
      "4180  8652    0:37.4 -3.897877916  0.506543868   \n",
      "4200  8700    0:37.4 -3.759561273  0.506024096   \n",
      "4220  8752    0:37.5 -3.62655309   0.505268199   \n",
      "4240  8836    0:37.6 -3.495091471  0.502607871   \n",
      "4260  8892    0:37.6 -3.367122903  0.50164861    \n",
      "4280  8960    0:37.7 -3.243861381  0.5           \n",
      "4300  9012    0:38.6 -3.123581891  0.499303298   \n",
      "4320  9036    0:38.6 -3.005269354  0.500231589   \n",
      "4340  9064    0:38.6 -2.891499668  0.500923361   \n",
      "4360  9088    0:38.7 -2.782204296  0.501841621   \n",
      "4380  9116    0:38.7 -2.675514353  0.502524094   \n",
      "4400  9140    0:38.8 -2.570142208  0.503432494   \n",
      "4420  9176    0:38.8 -2.467872788  0.503646308113\n",
      "4440  9212    0:38.9 -2.369097611  0.503858375   \n",
      "4460  9248    0:38.9 -2.274235917  0.504068716094\n",
      "4480  9284    0:39.0 -2.182427352  0.504277353   \n",
      "4500  9328    0:39.1 -2.092884827  0.504032258   \n",
      "4520  9348    0:39.1 -2.005826719  0.505140814   \n",
      "4540  9376    0:39.1 -1.921487903  0.505793226   \n",
      "4560  9404    0:39.2 -1.839416046  0.506441582   \n",
      "4580  9428    0:39.2 -1.766283029  0.507310589   \n",
      "4600  9464    0:39.3 -1.68907854   0.507502207   \n",
      "4620  9484    0:39.3 -1.625758745  0.508586526   \n",
      "4640  9512    0:39.4 -1.552340948  0.509218613   \n",
      "4660  9536    0:39.9 -1.481584005  0.510070053   \n",
      "4680  9564    0:39.9 -1.413641665  0.51069402    \n",
      "4700  9596    0:40.0 -1.348443856  0.511091779   \n",
      "4720  9620    0:40.0 -1.285759203  0.511930586   \n",
      "4740  9652    0:40.1 -1.225563026  0.51232166    \n",
      "4760  9692    0:40.1 -1.167871931  0.512268618   \n",
      "4780  9732    0:40.2 -1.11271321   0.512216031   \n",
      "4800  9772    0:40.3 -1.059999624  0.512163892   \n",
      "4820  9796    0:40.3 -1.009452042  0.512984249   \n",
      "4840  9824    0:40.4 -0.960922     0.513582343   \n",
      "4860  9844    0:40.4 -0.914457     0.514612452   \n",
      "4880  9872    0:40.5 -0.870034     0.515202703   \n",
      "4900  9900    0:40.5 -0.827619     0.515789474   \n",
      "4920  9920    0:40.6 -0.787142     0.516806723   \n",
      "4940  9952    0:40.6 -0.748524     0.517169179   \n",
      "4960  9984    0:40.6 -0.711621     0.517529215   \n",
      "4980  10016   0:40.7 -0.676434     0.517886855   \n",
      "5000  10044   0:40.7 -0.642811     0.518457072   \n",
      "5020  10080   0:40.8 -0.610705     0.518595041   \n",
      "5040  10128   0:40.8 -0.580168     0.518092105   \n",
      "5060  10176   0:40.9 -0.551077     0.517594108   \n",
      "5080  10212   0:40.9 -0.524139     0.517733388   \n",
      "Convergence obtained with Delta_z = -0.4989821806669852\n",
      "Done!\n"
     ]
    }
   ],
   "source": [
    "samples = sampler.run()\n",
    "print('Done!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot posterior samples versus true parameter values (dashed lines)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAGoCAYAAAATsnHAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xu0HGWd7//3h6AkQQVCCJOBhCAT\nEOFIgA2IeNlchRwBccJlRCZcJMjloMj5SUQGGC+zOIoT8TAiQYHAcBGGW9SICfGEMIOYGxESIJAB\nJCFZhCDKJUgkfH9/1LNJk+xL9d7dXdW7P6+19urq6uqqTxVhf3dVPfU8igjMzMzKZpOiA5iZmXXG\nBcrMzErJBcrMzErJBcrMzErJBcrMzErJBcrMzErJBcrMzErJBcrMzErJBcrMzEpp03qtWNK1wGeA\nVRGxe5p3KXA68GJa7MKImJY++zpwGrAOODcift3TNoYOHRqjRo2qfXgzM6vK/PnzV0fENrVcZ90K\nFHA9cCVwwwbzJ0XE5ZUzJH0YOAHYDfhb4D5JO0fEuu42MGrUKObNm1e7xGZm1iuS/lDrddbtEl9E\nzAb+mHPxo4FbI+LNiHgGWArsW69sZmZWfkXcgzpH0iOSrpW0VZq3HbCsYpnlad5GJE2QNE/SvBdf\nfLGzRczMrB9odIG6CtgJGAOsBL6f5quTZTvtZj0iJkdEW0S0bbNNTS93mpVCe3s77e3tRccwK1xD\nC1REvBAR6yLibeAa1l/GWw6MqFh0e2BFI7OZmVm5NLRASRpe8fYYYFGangqcIGkzSTsCo4E5jcxm\nZmblUs9m5rcA7cBQScuBS4B2SWPILt89C5wBEBGLJd0GPAa8BZzdUws+s/7quOOOKzqCWSmomUfU\nbWtrCzczNzMrnqT5EdFWy3XW8zkoM+vEpBlPdvv5GQdsD8DgwYMbEcestFygzEpm7NixAMyaNavY\nIGYFc198ZmZWSi5QZmZWSi5QZmZWSi5QZmZWSm4kYVYyJ598ctERzErBBcqsZFygzDIuUGY11NMz\nTnmsXr0agKFDh/Z5XWbNzAXKrGTGjRsH+DkoMzeSMDOzUnKBMjOzUvIlPrMq1OIek5nl4wJl1mTy\nFMnzDt25AUnM6ssFyqxkzjzzzKIjmJVCrgIlafeIWNTzkmbWV8cff3zREcxKIe8Z1I8lvRe4Hrg5\nIv5Uv0hm9dPT5bEyXBpbtmwZACNGjCg4iVmxchWoiPi4pNHAqcA8SXOA6yJiRl3TmbWgk046CfBz\nUGa5m5lHxFPARcAFwKeAH0p6QtLn6hXOzMxaV64CJekjkiYBjwMHAUdGxK5pelId85mZWYvKew/q\nSuAa4MKIeKNjZkSskHRRXZKZmVlLy1ugxgJvRMQ6AEmbAAMjYk1E3Fi3dGZm1rLy3oO6DxhU8X5w\nmtclSddKWiVpUcW8IZJmSHoqvW6V5kvSDyUtlfSIpL2q3RGz/uL888/n/PPPLzqGWeHyFqiBEfFa\nx5s0PbiH71wPHL7BvInAzIgYDcxM7wGOAEannwnAVTlzmfU7Rx55JEceeWTRMcwKl7dAvV55ViNp\nb+CNbpYnImYDf9xg9tHAlDQ9BfhsxfwbIvMQsKWk4TmzmfUrS5YsYcmSJUXHMCtc3ntQXwFul7Qi\nvR8O9OZx920jYiVARKyUNCzN3w5YVrHc8jRv5YYrkDSB7CyLkSNH9iKCWbmdccYZgJ+DMsv7oO5c\nSR8CdgEEPBERf61hDnW22S6yTAYmA7S1tXW6jJmZNb9qOovdBxiVvrOnJCLihiq394Kk4ensaTiw\nKs1fDlT267I9sGKjb5uZWcvI21nsjcBOwEJgXZodQLUFaiowHrgsvd5TMf8cSbcC+wF/7rgUaGbl\n4yE/rBHynkG1AR+OiNyX1CTdArQDQyUtBy4hK0y3SToNeA44Ni0+jexZq6XAGuCUvNsxM7P+KW+B\nWgT8DZ00WuhKRPxDFx8d3MmyAZydd91m9VKGEXMvusids5hB/gI1FHgs9WL+ZsfMiDiqLqnMWtgh\nhxzS53X0tdD68pyVQd4CdWk9Q5jZegsXLgRgzJgxBScxK1beZub3S9oBGB0R90kaDAyobzSz1vSV\nr3wF8HNQZnmH2zgd+A/g6jRrO+DueoUyMzPL29XR2cABwCvwzuCFw7r9hpmZWR/kvQf1ZkSslbIO\nHyRtShc9PZhZ3yx/OevmsgwtCs2KlPcM6n5JFwKDJB0K3A78vH6xzMys1eU9g5oInAY8CpxB9mDt\nT+oVyqyVjT3lvKIjmJVC3lZ8b5MN+X5NfeOY2Y67ebxOM8jfF98zdHLPKSI+WPNEZi3umcULABcq\ns2r64uswkKwPvSG1j2PWe/2lUcG06yYBcPblNxacxKxYuRpJRMRLFT/PR8QPgIPqnM3MzFpY3kt8\nldcaNiE7o3p/XRKZmZmR/xLf9yum3wKeBY6reRozM7Mkbyu+A+sdxMzMrFLeS3xf7e7ziPjX2sQx\ns8+eeWHREcxKoZpWfPuQDc0OcCQwG1hWj1BmrWy7nXYtOoJZKVQzYOFeEfEqgKRLgdsj4ov1CmbW\nqp5c8CAAO+/1sYKTmBUrb4EaCayteL8WGFXzNGbGjJuvAlygzPIWqBuBOZLuIutR4hjghrqlMjOz\nlpe3Fd93JP0K+ESadUpEPFy/WGZm1urynkEBDAZeiYjrJG0jaceIeKY3G5X0LPAqsA54KyLaJA0B\nfkZ26fBZ4LiIeLk367fm1F+6KjKz2sg75PslwAXA19Os9wD/3sdtHxgRYyKio5+/icDMiBgNzEzv\nzcysReU9gzoG2BNYABARKyTVuqujo4H2ND0FmEVWFM1ayrFf/ueiI/hs1kohb4FaGxEhKQAkbd7H\n7QYwPa3v6oiYDGwbESsBImKlpGGdfVHSBGACwMiRI/sYw6x8ho3wKDZmkL9A3SbpamBLSacDp9K3\nwQsPSGdhw4AZkp7I+8VUzCYDtLW1bTRGlVmzW/zb3wCw2/4eMMBaW95WfJdLOhR4BdgFuDgiZvR2\noxGxIr2uSk3X9wVekDQ8nT0NB1b1dv1mzWzWHdcB/b9A9XQZ8bxDd25QEiurHguUpAHAryPiEKDX\nRalifZsDm0TEq2n6MOCbZN0ojQcuS6/39HVbZmbWvHosUBGxTtIaSVtExJ9rsM1tgbskdWz/5oi4\nV9JcskuJpwHPkY3aa2ZmLSrvPai/AI9KmgG83jEzIs6tdoMR8TSwRyfzXwIOrnZ91hzcKqz1+L+5\n9VXeAvXL9GNmZtYQ3RYoSSMj4rmImNKoQGat7vNf+27REcxKoaeeJO7umJB0R52zmBmw1bDhbDVs\neNExzArXU4FSxbSfHjRrgIdnTePhWdOKjmFWuJ7uQUUX02ZWJw/+4hYA9mwfW3ASs2L1VKD2kPQK\n2ZnUoDRNeh8R8YG6pjMzs5bVbYGKiAGNCmJmZlYp13AbZmZmjVbNgIXWwtxvmpk1mguU1YR7Daid\n8f90RdERzErBBcqsZN63xZCiI5iVgu9BmZXMnOl3Mmf6nUXHMCucC5RZycydfhdzp99VdAyzwvkS\nn5k1LTfe6d9coMyslNzwxnyJz8zMSskFyszMSsmX+Azw5ZQyOf3bk4uOYFYKLlBmJfPegYOKjtBv\nuBFFc3OBagE+O2ou/zX1JgAOOOrEgpOYFcsFyqxkFs6+F3CBagSfYZWbC1Q/4DMkM+uPSlegJB0O\nXAEMAH4SEZcVHKlwLkBm1opKVaAkDQD+DTgUWA7MlTQ1Ih4rNpmZtaI8fxz2dBnQlxF7r1QFCtgX\nWBoRTwNIuhU4GqhLgWrEP76+rt/Myq2vvwNq8Xuor9so6+8hRUTRGd4haRxweER8Mb0/CdgvIs6p\nWGYCMCG93QVY0odNDgVW9+H7jdAMGcE5a60ZcjZDRnDOWusq5w4RsU0tN1S2Myh1Mu9dFTQiJgM1\neZJR0ryIaKvFuuqlGTKCc9ZaM+RshozgnLXWyJxl6+poOTCi4v32wIqCspiZWYHKVqDmAqMl7Sjp\nvcAJwNSCM5mZWQFKdYkvIt6SdA7wa7Jm5tdGxOI6brIZOj1rhozgnLXWDDmbISM4Z601LGepGkmY\nmZl1KNslPjMzM8AFyszMyioimuoHOJzs2aelwMROPt8M+Fn6/HfAqDR/a+D/Aa8BV3ax7qnAoor3\nY4CHgIXAPGDfNF/AD9M2HgH2KmHGduDPaf5C4OKCj+UewG+BR4GfAx+o+OzraRtLgE+XMScwCnij\n4nj+uBE5gVlpnR3bHdbduno6nmXIWMJj+UlgAfAWMG6D74wHnko/40ucc13F8lMLzvlVss4VHgFm\nkj0flet4bpSrpwXK9EPWcOK/gQ8C7wV+D3x4g2XO6vgHT9YK8GdpenPg48CXNjzQ6fPPATfz7l9W\n04Ej0vRYYFbF9K/ICtVHgd+VMGM78IsSHcu5wKfS9KnAt9L0h9O2NwN2TJkGlDDnqMrlGnU8yX4J\ntHWyva7W1eXxLFHGsh3LUcBHgBuo+MUPDAGeTq9bpemtypYzffZaiY7ngcDgNH1mxbq6PZ6d/TTb\nJb53ukKKiLVAR1dIlY4GpqTp/wAOlqSIeD0i/hP4y4YrlfQ+sqr/7Q0+CuADaXoL1j+TdTRwQ2Qe\nAraUNLxkGXvS6Jy7ALPT9Azg7yu2cWtEvBkRz5D9BbdvCXP2pC45u9Hpuuj+eJYlY08amjMino2I\nR4C3N/jo08CMiPhjRLxM9u/h8BLm7Emjc/6/iFiT3j5E9jwr9Hw8N9JsBWo7YFnF++VpXqfLRMRb\nZJe5tu5hvd8Cvg+s2WD+V4DvSVoGXE526aSnHGXJCLC/pN9L+pWk3brK0KCci4Cj0vSxrH8gu6cc\nZckJsKOkhyXdL+kTXWWocU6A6yQtlPRPFb/gu1pXEf82q80I5TqWXSnq32a1OQEGSpon6SFJn61y\nP+qZ8zSyq015c7xLsxWoHrtCyrnM+oWlMcDfRcRdnXx8JnBeRIwAzgN+mmMbZcm4gOza7x7A/wXu\n7kWGWuY8FThb0nzg/cDanNsoS86VwMiI2JPszOtmSR+o+F7NcyYnRsT/AD6Rfk7qYV0N/bfZy4xl\nO5Zdafi/zaTanJAdzzbg88APJO1UdE5JXwDagO/1dhvNVqDydIX0zjKSNiW77PXHbta5P7C3pGeB\n/wR2ljQrfTYeuDNN3876SyXd5ShFxoh4JSJeS9PTgPdIGtpZhkbkjIgnIuKwiNgbuIXsmnieHKXI\nmS6ZvZSm56f5lV1A1yMnEfF8en2V7H7ZRv8GN1hXo/9tVp2xhMeyK0X82+xNTiJiRXp9muze0J5F\n5pR0CPAN4KiIeLOKHO/SbAUqT1dIU8l+aQOMA34T6Q5dZyLiqoj424gYRXYz8MmIaE8frwA+laYP\nImt50rGNf1Tmo8CfI2JlmTJK+puOU25J+5L9t36pYrUNzSlpWHrdBLgI+HHFNk6QtJmkHYHRwJyy\n5ZS0jbLxypD0wZTz6XrmlLRpxx8Vkt4DfIbsEmR36+rueJYiYwmPZVd+DRwmaStJWwGHpXmlypny\nbZamhwIH8O4hihqaU9KewNVkxWlVxdd6Op4bi25aUJTxh6yl2pNkf3V9I837ZjoYAAPJziSWkv2P\n+cGK7z5L9lfBa2TVfMOWLKN4d4uujwPzyVq9/A7YO80X2cCK/03WHLmthBnPARan+Q8BHyv4WH45\nbetJ4DJSLybps2+kDEtILRLLlpOssUTH8VwAHFnvnGQtqOaTNdddTBppOse6ujyeZchYwmO5T1ru\ndbI/4hZXrOvUtI2lwCllzAl8jOz30O/T62kF57wPeIFOmr33dDw3/HFXR2ZmVkrNdonPzMxahAuU\nmZmVkguUmZmVkguUmZmVkguUmZmVkguUWSJpXeq2ZZGk2yUNLjoTgKQLa7COYyUtlvS2pLZa5DKr\nNxcos/XeiIgxEbE7WRdHX8r7xY4HT+uk6gLVSZ5FZL23z+5kcbNScoEy69wDwN8BSLpb0vx0BjKh\nYwFJr0n6pqTfkXXMe7GkuekMbHJFTx6zJE2SNFvS45L2kXSnpKckfbtifV+QNCedxV0taYCky4BB\nad5NXS3XWZ7KnYmIxyNiSb0PmlktuUCZbSD1RXYE2VP5AKdG1jdfG3CupI5enjcn64Viv8iGJLgy\nIvZJZ2CDyLp/6bA2Ij5J1nXSPcDZwO7AyZK2lrQrcDxwQESMIRuA7sSImMj6M7sTu1quizxmTW3T\nogOYlcggSQvT9AOs7xn+XEnHpOkRZH3HvURWHO6o+P6Bkr4GDCYblG0x2ai8sL7vs0fJuqhZCSDp\n6bTOjwN7A3PTidcgoLIfsw4Hd7PchnnMmpoLlNl6b6SzkndIagcOAfaPiDXKej0fmD7+S0SsS8sN\nBH5E1i/jMkmXViwH0NGj89sV0x3vNyXr33FKRFSO59WZ7pZ7J49Zf+BLfGbd2wJ4ORWnDwEf7WK5\njmK0WtloveOq3M5MYFxFb+pDJO2QPvtr6jG6p+XM+hUXKLPu3QtsKukRslF4H+psoYj4E3AN2SW8\nu8mGOMgtIh4jG95jetrWDGB4+ngy8Iikm3pYrkuSjpG0nKzxxC8ldT/MgVkJuDdzMzMrJZ9BmZlZ\nKblAmZlZKblAmZlZKblAmZlZKblAmZlZKblAmZlZKblAmZlZKblAmZlZKblAmZlZKblAmZlZKblA\nmZlZKblAmZlZKblAmZlZKblAmZlZKTX1iLpDhw6NUaNGFR3DzKzlzZ8/f3VEbFPLdTZ1gRo1ahTz\n5s0rOoaZWcuT9Idar9OX+MzMrJRcoMzMrJTqVqAkXStplaRFFfMulfS8pIXpZ2zFZ1+XtFTSEkmf\nrlcuMzNrDvU8g7oeOLyT+ZMiYkz6mQYg6cPACcBu6Ts/kjSgjtnMWkJ7ezvt7e1FxzDrlboVqIiY\nDfwx5+JHA7dGxJsR8QywFNi3XtnMzKz8imjFd46kfwTmAedHxMvAdsBDFcssT/M2ImkCMAFg5MiR\ndY5q1tyOO+64oiOY9VqjG0lcBewEjAFWAt9P89XJstHZCiJickS0RUTbNtvUtMm9Wb9z1llncdZZ\nZxUdw6xXGnoGFREvdExLugb4RXq7HBhRsej2wIoGRjOriUkznuxxmfMO3bkBSTJr1qwBYPDgwQ3b\nplmtNPQMStLwirfHAB0t/KYCJ0jaTNKOwGhgTiOzmfVHY8eOZezYsT0vaFZCdTuDknQL0A4MlbQc\nuARolzSG7PLds8AZABGxWNJtwGPAW8DZEbGuXtnMzKz86lagIuIfOpn9026W/w7wnXrlMTOz5uKe\nJMzMrJRcoMzMrJSaujdzM+veySefXHQEs15zgTLrx1ygrJn5Ep9ZP7Z69WpWr15ddAyzXvEZlFnJ\n9PSwb54HfTvW8W//+yQAzr78xqrXYVY0n0GZmVkpuUCZmVkpuUCZmVkpuUCZmVkpuZGEWT/2sc90\n1uOYWXPIVaAk7R4Ri3pe0szKZM9292RuzSvvJb4fS5oj6SxJW9Y1kZnVzMurVvLyqpVFxzDrlVwF\nKiI+DpxINqjgPEk3Szq0rsnMrM9u/u7XuPm7Xys6hlmv5G4kERFPARcBFwCfAn4o6QlJn6tXODMz\na125CpSkj0iaBDwOHAQcGRG7pulJdcxnZmYtKm8rviuBa4ALI+KNjpkRsULSRXVJZmZmLS1vgRoL\nvNExDLukTYCBEbEmIm7s/qtmZmbVy1ug7gMOAV5L7wcD04GP1SOUmdVG+9+fUnQEs17LW6AGRkRH\ncSIiXpM0uE6ZzKxGdtv/oKIjmPVa3gL1uqS9ImIBgKS9gTd6+I6Z1UFPw3FUWrXsaQCGjfhgveKY\n1U3eZuZfAW6X9ICkB4CfAed09wVJ10paJWlRxbwhkmZIeiq9bpXmS9IPJS2V9IikvXq7Q2a23u1X\nXMLtV1xSdAyzXsn7oO5c4EPAmcBZwK4RMb+Hr10PHL7BvInAzIgYDcxM7wGOAEannwnAVXlymZlZ\n/1VNZ7H7AKPSd/aURETc0NXCETFb0qgNZh8NtKfpKcAssgd/jwZuiIgAHpK0paThEeE+WqzfqeYS\nnVkry9tZ7I3ATsBCYF2aHUCXBaoL23YUnYhYKWlYmr8dsKxiueVp3kYFStIEsrMsRo4cWeXmzfqm\nVYpLLYadN+urvGdQbcCH0xlOPaiTeZ1uKyImA5MB2tra6pXHWlSrFKC+ynOcXMSsr/IWqEXA39DJ\nGU2VXui4dCdpOLAqzV9O1hFth+2BFX3cllnLO/TzZxYdwazX8haoocBjkuYAb3bMjIijqtzeVGA8\ncFl6vadi/jmSbgX2A/7s+09mfbfzXn6W3ppX3gJ1abUrlnQLWYOIoZKWA5eQFabbJJ0GPAccmxaf\nRtad0lJgDeDH381q4Pn/fhyA7XbateAkZtXLVaAi4n5JOwCjI+K+1IvEgB6+09VY0wd3smwAZ+fJ\nYmb53X3VvwBw9uXuMtOaT97hNk4H/gO4Os3aDri7XqHMzMzy9iRxNnAA8Aq8M3jhsG6/YWZm1gd5\nC9SbEbG2442kTemiGbiZmVkt5C1Q90u6EBgk6VDgduDn9YtlZmatLm8rvonAacCjwBlkre5+Uq9Q\nZlYbY085r+gIZr2WtxXf22RDvl9T3zhmVks77uaBAax55e2L7xk6uecUER5kxqzEnlm8AHChsuZU\nTV98HQaSPWA7pPZxzKyWpl03Cdj4OahG9DnoDmetr/KOB/VSxc/zEfEDwGNJm5lZ3eS9xFd5fWAT\nsjOq99clkZmZGfkv8X2/Yvot4FnguJqnMTMzS/K24juw3kHMzMwq5b3E99XuPo+If61NHDOrpc+e\neWHREcx6rZpWfPuQjdsEcCQwm3cP025mJVPmYTbcys96Us2AhXtFxKsAki4Fbo+IL9YrmFk9tNqQ\n7k8ueBDwwIXWnPIWqJHA2or3a4FRNU9jZjU14+arABcoa055C9SNwBxJd5H1KHEMcEPdUpmZWcvL\n24rvO5J+BXwizTolIh6uXywzM2t1ec+gAAYDr0TEdZK2kbRjRDxTr2Bm1Wq1+0tm/V3eId8vAS4A\nvp5mvQf493qFMjMzy3sGdQywJ7AAICJWSOp1V0eSngVeBdYBb0VEm6QhwM/IGl88CxwXES/3dhtm\nBsd++Z+LjmDWa3lH1F0bEUEackPS5jXY9oERMSYiOnpKnwjMjIjRwMz03sz6YNiIDzJshEfFseaU\nt0DdJulqYEtJpwP3UfvBC48GpqTpKcBna7x+s5az+Le/YfFvf1N0DLNeyduK73JJhwKvALsAF0fE\njD5sN4DpkgK4OiImA9tGxMq0vZWShnX2RUkTgAkAI0eO7EMEs/5v1h3XAbDb/h4dx5pPjwVK0gDg\n1xFxCNCXolTpgHQfaxgwQ9ITeb+YitlkgLa2to1G+TUzs/6hx0t8EbEOWCNpi1ptNCJWpNdVwF3A\nvsALkoYDpNdVtdqemZk1n7yt+P4CPCppBvB6x8yIOLfaDaYGFptExKtp+jDgm2Qd0Y4HLkuv91S7\nbjMz6z/yFqhfpp9a2Ba4S1LH9m+OiHslzSVrjHEa8BxwbI22Z03CvVubWaVuC5SkkRHxXERM6W65\nakTE08Aencx/CTi4VtsxM/j8175bdASzXuvpHtTdHROS7qhzFjOrsa2GDWerYcOLjmHWKz0VKFVM\n+2k/sybz8KxpPDxrWtExzHqlp3tQ0cW0mTWBB39xCwB7to8tOIlZ9XoqUHtIeoXsTGpQmia9j4j4\nQF3TmZlZy+q2QEXEgEYFMTMzq5S3Lz4zM7OGcoEyM7NSqmZEXTNrMuP/6YqiI5j1mguUNQ0P6V69\n920xpOgIZr3mAmXWj82ZficA+x72uYKT1Ie7x+rfXKDM+rG50+8CmrNA+YzZ3EjCzMxKyQXKzMxK\nyZf4rCZ8OcbMas0FysxaVi3+sHJDjPpxgTLrx07/9uSiI5j1mguUWT/23oGDio5g1mtuJGHWj/3X\n1Jv4r6k3FR3DrFd8BmXWjy2cfS8ABxx1YsFJiuHGO83NBcoA/49s1lvuzaJ+SlegJB0OXAEMAH4S\nEZcVHMnMrNdcwHqvVAVK0gDg34BDgeXAXElTI+KxYpOZmdWHm7p3rVQFCtgXWBoRTwNIuhU4Gui3\nBaoRl9b66z9eM8v017M0RUTRGd4haRxweER8Mb0/CdgvIs6pWGYCMCG93QVYUsUmhgKraxS32Xjf\nW1cr77/3vXF2iIhtarnCsp1BqZN576qgETEZ6NXTh5LmRURbb77b7Lzvrbnv0Nr7731v7n0v23NQ\ny4ERFe+3B1YUlMXMzApUtgI1FxgtaUdJ7wVOAKYWnMnMzApQqkt8EfGWpHOAX5M1M782IhbXcBOt\n3DGZ9711tfL+e9+bWKkaSZiZmXUo2yU+MzMzwAXKzMxKqt8VKEkDJD0s6Rfp/UGSFkhaJGmKpE03\nWH4fSevSM1hNrZp9l9QuaaGkxZLuLy51beTdd0lbSPq5pN+nfT+l2OR9J+lZSY+m/57z0rwhkmZI\neiq9bpXmS9IPJS2V9IikvYpN3zdV7vuJaZ8fkfSgpD2KTd931ex/xXea5ndevytQwJeBxwEkbQJM\nAU6IiN2BPwDjOxZMXSv9H7JGGf1Brn2XtCXwI+CoiNgNOLaYuDWV97/72cBjEbEH0A58P7UYbXYH\nRsSYiudeJgIzI2I0MDO9BzgCGJ1+JgBXNTxp7eXd92eAT0XER4Bv0Q8aESR597/pfuf1qwIlaXvg\nfwI/SbO2Bt6MiI5+QGYAf1/xlf8F3AGsaljIOqly3z8P3BkRzwFERFPvf5X7HsD7JQl4H/BH4K0G\nxm2Uo8mKNOn1sxXzb4jMQ8CWkoYXEbCOOt33iHgwIl5O8x8ie86yP+rqvz002e+8flWggB8AXwPe\nTu9XA++R1PGXxTjSg8CStgOOAX7c6JB1knvfgZ2BrSTNkjRf0j82NmrNVbPvVwK7kj0A/ijw5Yh4\nm+YWwPT037KjG7BtI2IlQHodluZvByyr+O7yNK9ZVbPvlU4DftWgjPWUe/+b8XdeqZ6D6gtJnwFW\nRcR8Se0AERGSTgAmSdoMmM76v5Z/AFwQEeuyP6abVy/2fVNgb+BgYBDwW0kPVZxxNI1e7PungYXA\nQcBOwAxJD0TEK41PXzMHRMT5G2bWAAATCUlEQVQKScPI9ueJbpbtsTuxJlPNvgMg6UCyAvXxuqer\nv2r2v+l+5/WbAgUcABwlaSwwEPiApH+PiC8AnwCQdBjZ2QNAG3Br+g81FBgr6a2IuLvx0fus2n1f\nDqyOiNeB1yXNBvYAmq5AUf2+nwJcFtkDgEslPQN8CJjT+Oi1EREr0usqSXeRjQrwgqThEbEyXcLr\nuKTTr7oTq3LfkfQRskvBR0TES4WErqEq97/5fudFRL/7Ibv5/Ys0PSy9bkZ2w/CgTpa/HhhXdO5G\n7TvZJa6ZZH+gDAYWAbsXnb1B+34VcGma3hZ4HhhadPY+7PPmwPsrph8EDge+B0xM8ycC303T/5Ps\n0paAjwJzit6HBu77SGAp8LGisxex/xt8tyl+5/WnM6iu/H/pMtAmwFUR8ZuiAzVQp/seEY9Luhd4\nhOy+zU8iYlGBOeuhq//u3wKul/Qo2S/pCyKimYdj2Ba4K/1VvClwc0TcK2kucJuk04DnWN9Scxow\nluwX9RqyM8pmVe2+X0zWgOZH6TtvRXP39l3t/jcdd3VkZmal1N9a8ZmZWT/hAmVmZqXkAmVmZqXk\nAmVmZqXkAmVmZqXkAmWWpB6eFyrrAf12SYOLzgQg6cIarON7kp5IPXnflToMNis1Fyiz9d6IrFfo\n3YG1wJfyfjH1El0vVReoTvLMIHsY+yNkPYZ8vRbBzOrJBcqscw8Afwcg6e7UGefiig45kfSapG9K\n+h2wv6SLJc1NZ2CTU4/ppE55J0maLenxNB7PnWm8nm9XrO8Lkuaks7irlY1xdRkwKM27qavlOstT\nuTMRMT0iOvoj7M89eVs/4gJltgFlgxseQdbbOcCpEbE3WV9m50raOs3fHFgUEftFxH8CV0bEPukM\nbBDwmYrVro2IT5L1JH0P2bhUuwMnS9pa0q7A8WSdf44B1gEnRsRE1p/ZndjVcl3k6cqp9I+evK2f\na4WujszyGiRpYZp+APhpmj5X0jFpegTZYH8vkRWHOyq+f6Ckr5H1bzgEWAz8PH02Nb0+CiyONByC\npKfTOj9O1sP83HTiNYjOx+w5uJvlNsyzEUnfIOvZ/abuljMrAxcos/XeSGcl70hDeBwC7B8RayTN\nIus1HeAvEbEuLTeQbJTitohYJunSiuUA3kyvb1dMd7zflKxfwCkR0dO9oe6WeydPp1+UxpOd1R0c\n7uPMmoAv8Zl1bwvg5VScPkTWA3hnOorRaknvIxsksRozgXFpXB8kDZG0Q/rsr5Lek2O5Lkk6HLgA\nOCoi1lSZzawQPoMy6969wJckPQIsIWtgsJGI+JOka8gu4T0LzK1mIxHxmKSLyEZH3QT4K9l9qj8A\nk4FHJC1I96G6Wq47V5INPTIjXRp8KCJyt1I0K4J7Mzczs1LyJT4zMyslFygzMyslFygzMyslFygz\nMyslFygzMyslFygzMyslFygzMyslFygzMyslFygzMyslFygzMyslFygzMyslFygzMyslFygzMysl\nFygzMyslFygzMyulph6wcOjQoTFq1KiiY5iZtbz58+evjohtarnOpi5Qo0aNYt68eUXHMDNreZJ6\nGtW5ar7EZ2ZmpeQCZWZmpeQCZdYC2tvbaW9vLzqGWVVcoMzMrJRcoMzMrJSauhWfmeVz3HHHFR3B\nrGouUGYt4Kyzzio6glnVXKCsaUya8WS3n5936M4NStJ81qxZA8DgwYMLTmKWnwuUWQsYO3YsALNm\nzSo2iFkVXKDMquCzOLPGcSs+MzMrJRcoMzMrJV/iM6uhni4Bgi8DmuXlAmUtpVXvIZ188slFRzCr\nmguUWYU8Z0DNyAXKmlHd7kFJulbSKkmLKuYNkTRD0lPpdas0X5J+KGmppEck7VWvXGataPXq1axe\nvbroGGZVqWcjieuBwzeYNxGYGRGjgZnpPcARwOj0MwG4qo65zFrOuHHjGDduXNExzKpSt0t8ETFb\n0qgNZh8NtKfpKcAs4II0/4aICOAhSVtKGh4RK+uVz6worXofzKxajW5mvm1H0Umvw9L87YBlFcst\nT/M2ImmCpHmS5r344ot1DWtmZsUpy3NQ6mRedLZgREyOiLaIaNtmm23qHMvMzIrS6AL1gqThAOl1\nVZq/HBhRsdz2wIoGZzMzsxLJdQ9K0u4RsajnJXs0FRgPXJZe76mYf46kW4H9gD/7/pNZ7Zx55plF\nRzCrWt5GEj+W9F6ylnk3R8SfevqCpFvIGkQMlbQcuISsMN0m6TTgOeDYtPg0YCywFFgDnFLFPphZ\nD44//viiI5hVLVeBioiPSxoNnArMkzQHuC4iZnTznX/o4qODO1k2gLPzZDGz6i1blrVBGjFiRA9L\nmpVH7mbmEfGUpIuAecAPgT0lCbgwIu6sV0Az67uTTjoJ8HhQ1lxyNZKQ9BFJk4DHgYOAIyNi1zQ9\nqY75zMysReU9g7oSuIbsbOmNjpkRsSKdVZmZmdVU3gI1FngjItYBSNoEGBgRayLixrqlM6tCf+3o\n1axV5X0O6j5gUMX7wWmemZlZXeQ9gxoYEa91vImI1yQNrlMmM6ux888/v+gIZlXLW6Bel7RXRCwA\nkLQ38EYP37EW4g5Qy+3II48sOoJZ1fIWqK8At0vq6H5oOOAn/8yaxJIlSwDYZZddCk5ill/eB3Xn\nSvoQsAtZx65PRMRf65rMzGrmjDPOAPwclDWXasaD2gcYlb6zpyQi4oa6pDIzs5aXt7PYG4GdgIXA\nujQ7ABcoswbL05ze9/ysP8h7BtUGfDj1mWdmZlZ3eZ+DWgT8TT2DmJmZVcp7BjUUeCz1Yv5mx8yI\nOKouqcyspi66yD2SWfPJW6AurWcIM3BXRbW08bEcCcCjab7vUVkzyNvM/H5JOwCjI+K+1IvEgPpG\ns/7ExadYz//34wBst9OuBScxyy/vcBunA/8BXJ1mbQfcXa9QZlZbd1/1L9x91b8UHcOsKnkbSZwN\nHAC8AtnghcCweoUyMzPLW6DejIi1HW8kbUr2HJSZmVld5G0kcb+kC4FBkg4FzgJ+Xr9YZq3L9+vM\nMnnPoCYCLwKPAmcA0wC3WzUzs7rJ24rvbbIh36+pbxwzq4exp5xXdASzquXti+8ZOrnnFBEfrHki\nM6u5HXfbq+gIZlWrpi++DgOBY4EhtY9jZvXwzOIFgAuVNZdc96Ai4qWKn+cj4gfAQXXOZmY1Mu26\nSUy7blLRMcyqkvcSX+WfXZuQnVG9vy6JzMzMyH+J7/sV028BzwLH1TyNmZlZkrcV34H1DmJmZlYp\n7yW+r3b3eUT8azUblfQs8CrZ6LxvRUSbpCHAz8iGlX8WOC4iXq5mvWZm1n9U04pvH2Bqen8kMBtY\n1odtHxgRqyveTwRmRsRlkiam9xf0Yf1mlnz2zAuLjmBWtWoGLNwrIl4FkHQpcHtEfLGGWY4G2tP0\nFGAWLlBmNeFhNqwZ5e3qaCSwtuL9WrJLcb0VwHRJ8yVNSPO2jYiVAOnVvaWb1ciTCx7kyQUPFh3D\nrCp5z6BuBOZIuousuBwD3NCH7R4QESskDQNmSHoi7xdTQZsAMHLkyD5EsEruoLR/m3HzVQDsvNfH\nCk5ill/eVnzfkfQr4BNp1ikR8XBvNxoRK9LrqlT09gVekDQ8IlZKGg6s6uK7k4HJAG1tbR7yw6wO\nevqDxUPGWyPkvcQHMBh4JSKuAJZL2rE3G5S0uaT3d0wDhwGLyBpgjE+LjQfu6c36zcysf8jbzPwS\nspZ8uwDXAe8B/p1slN1qbQvcJalj+zdHxL2S5gK3SToNeI6svz8zM2tRee9BHQPsCSyA7BJdx1lQ\ntSLiaWCPTua/BBzcm3WamVn/k7dArY2IkBTwzqU5M2sSx375n9/13o1irBnkLVC3Sboa2FLS6cCp\nePBCs6YxbISHbrPmk7cV3+WSDgVeIbsPdXFEzKhrMjOrmcW//Q0Au+3vUXKsefRYoCQNAH4dEYcA\nLkpmTWjWHdcBtStQeS4Ruim69VWPzcwjYh2wRtIWDchjZmYG5L8H9RfgUUkzgNc7ZkbEuXVJZWZm\nLS9vgfpl+jEzM2uIbguUpJER8VxETGlUIDMzM+j5DOpuYC8ASXdExN/XP5KZ1drnv/bdoiOYVa2n\nAqWKaT9IYdaktho2vOgIZlXrqRVfdDFtZk3k4VnTeHjWtKJjmFWlpzOoPSS9QnYmNShNk95HRHyg\nrunMrCYe/MUtAOzZPrZh2/SQHdZX3RaoiBjQqCBmZmaVqhkPyszMrGFcoMzMrJTyPqhrJeZr/WbW\nH7lAtQCP/WPj/+mKoiOYVc0FyqwFvG+LIUVHMKua70GZtYA50+9kzvQ7i45hVhWfQZm1gLnT7wJg\n38M+V3CS9Rpx6dn3X5ubz6DMzKyUXKDMzKyUXKDMzKyUXKDMzKyU3EiiCfg5Juur0789uegIZlVz\ngTJrAe8dOKjoCGZVc4EyawH/NfUmAA446sSCk5SLuwkrNxcosxawcPa9QOsVqEZcHneRqx8XqD7y\nP06z/sv3f4tVugIl6XDgCmAA8JOIuKzgSH3iAmZm1julKlCSBgD/BhwKLAfmSpoaEY/VY3v+68jM\nmkGr/qFbqgIF7AssjYinASTdChwN1KVAlYGLpFn/Vob/x5u1wCkiis7wDknjgMMj4ovp/UnAfhFx\nTsUyE4AJ6e0uwJKGBy2HocDqokOUlI9N53xcOufj0rVqjs0OEbFNLTdetjModTLvXRU0IiYDLf/U\noaR5EdFWdI4y8rHpnI9L53xculb0sSlbV0fLgREV77cHVhSUxczMClS2AjUXGC1pR0nvBU4Aphac\nyczMClCqS3wR8Zakc4BfkzUzvzYiFhccq6xa/jJnN3xsOufj0jkfl64VemxK1UjCzMysQ9ku8ZmZ\nmQEuUGZmVlIuUE1I0nmSFktaJOkWSQOLzlQGkr6cjsliSV8pOk+RJF0raZWkRRXzhkiaIemp9LpV\nkRmL0MVxOTb9m3lbUss2N+/i2HxP0hOSHpF0l6QtG5nJBarJSNoOOBdoi4jdyRqTnFBsquJJ2h04\nnaw3kj2Az0gaXWyqQl0PHL7BvInAzIgYDcxM71vN9Wx8XBYBnwNmNzxNuVzPxsdmBrB7RHwEeBL4\neiMDuUA1p02BQZI2BQbjZ8UAdgUeiog1EfEWcD9wTMGZChMRs4E/bjD7aGBKmp4CfLahoUqgs+MS\nEY9HRKv2SPOOLo7N9PT/E8BDZM+mNowLVJOJiOeBy4HngJXAnyNierGpSmER8ElJW0saDIzl3Q99\nG2wbESsB0uuwgvNYczkV+FUjN+gC1WTSfYOjgR2BvwU2l/SFYlMVLyIeB/4P2SWJe4HfA291+yUz\ny0XSN8j+f7qpkdt1gWo+hwDPRMSLEfFX4E7gYwVnKoWI+GlE7BURnyS7VPFU0ZlK5gVJwwHS66qC\n81gTkDQe+AxwYjT4wVkXqObzHPBRSYMlCTgYeLzgTKUgaVh6HUl20/uWYhOVzlRgfJoeD9xTYBZr\nAmkA2QuAoyJiTcO3754kmo+kfwaOJzvlfhj4YkS8WWyq4kl6ANga+Cvw1YiYWXCkwki6BWgnGy7h\nBeAS4G7gNmAk2R86x0bEhg0p+rUujssfgf8LbAP8CVgYEZ8uKmNRujg2Xwc2A15Kiz0UEV9qWCYX\nKDMzKyNf4jMzs1JygTIzs1JygTIzs1JygTIzs1JygTIzs1JygTJLJK2TtDD1iH576jKpcJIurME6\nvpV6pF4oabqkv61FNrN6cjNzs0TSaxHxvjR9EzA/Iv4153cHRMS6eueq4jvvyiPpAxHxSpo+F/hw\nI59nMesNn0GZde4B4O8AJN0taX4aM2hCxwKSXpP0TUm/A/aXdLGkuekMbHLq6QNJsyRNkjRb0uOS\n9pF0ZxqX6dsV6/uCpDnpLOdqSQMkXUbWc/3CVDQ7Xa6zPJU701Gcks0B/2VqpecCZbaBNIzJEcCj\nadapEbE30AacK2nrNH9zYFFE7BcR/wlcGRH7pHG6BpH1X9Zhbeoj8MdkXQydDewOnJx6YN+VrHeQ\nAyJiDLCOrO+zicAbETEmIk7sarku8my4X9+RtCwtf3ENDpVZXblAma03SNJCYB5ZV0A/TfPPlfR7\nsvFwRgAdAyGuA+6o+P6Bkn4n6VHgIGC3is+mptdHgcURsTJ1T/V0WufBwN7A3JThYOCDnWTsbrkN\n87xLRHwjIkaQ9Uh9TrdHwqwENi06gFmJvJHOSt4hqZ2sB/n9I2KNpFnAwPTxXzru80gaCPyIbKTj\nZZIurVgOoKOvxLcrpjvebwoImBIRPY1Y2t1yf8l5H+xm4Jdkfa2ZlZbPoMy6twXwcipOHwI+2sVy\nHcVotaT3AeOq3M5MYFxFj+xDJO2QPvurpPfkWK5LkkZXvD0KeKLKfGYN5zMos+7dC3xJ0iPAErLL\nfBuJiD9JuobsEt6zwNxqNhIRj0m6CJguaROyHtnPBv4ATAYekbQg3YfqarnuXCZpF7Iztj8AbsFn\npedm5mZmVkq+xGdmZqXkAmVmZqXkAmVmZqXkAmVmZqXkAmVmZqXkAmVmZqXkAmVmZqX0/wN0XEH/\nnrLOiwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x26df2a39e80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot output\n",
    "import pints.plot\n",
    "\n",
    "pints.plot.histogram([samples], ref_parameters=[r, k, sigma])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot posterior predictive simulations versus the observed data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAADQCAYAAAAK/RswAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3Xd8k+X6+PHPnaRpupsu6ABappAC\nRcoSBNEDCILinrjluI56/MJxK84jP/QcPW6POFAEERdHHDhAoDjY0LCxpbQFWiAt3W2S+/dH0tDS\nAGWUQrnerxevJnfu58mVtOS5ck+ltUYIIYQQoiUxNHcAQgghhBDHmyQ4QgghhGhxJMERQgghRIsj\nCY4QQgghWhxJcIQQQgjR4kiCI4QQQogWRxIcIYQQQrQ4kuAIIYQQosWRBEcIIYQQLY6puQM4FjEx\nMTo5Obm5wxBCCCHECbJ8+fLdWuvYw9U7pROc5ORkli1b1txhCCGEEOIEUUpta0w96aISQgghRIsj\nCY4QQgghWhxJcIQQ4iSTk5ODw+GoV+ZwOMjJyWmmiIQ49ZzSY3D8qampITc3l8rKyuYORbQwFouF\npKQkAgICmjsU0cKFhYVht9ux2WxYrVYcDofvvhCicVpcgpObm0tYWBjJyckopZo7HNFCaK3Zs2cP\nubm5pKSkNHc4ooWzWq3YbDbsdjuJiYnk5eX5kh1x+snJySEsLKze79/hcFBSUkLbtm2bMbKTW4vr\noqqsrCQ6OlqSG3FcKaWIjo6WlkFxwlitVhITE8nOziYxMfGkSG4a03Um3WvHX22LXu37arfbWbZs\nGWFhYb46R/seH+3v61T4Pbe4BAeQ5EY0Cfm7EieSw+EgLy+P5ORk8vLyGlxMmpq/C5jb7SYjI8NX\nXtt1VvdCe+DF2F+dpo7Tbrdjt9vrlZ1sF98jUbdFLysri/z8/HqPH8t77O/3lZGRgdvtbvAcdd+/\nE/17PhpNmuAopbKVUmuVUquUUsu8ZVFKqR+UUpu9P63ecqWU+o9SaotSao1S6symjE0IIU5Wdcfc\npKSkYLPZyMjIIDs7u0G9prpo+7uA5ebmkpqa6rvQ1h0nVOvAi7G/Oo3VmFYCf3Hm5+eTn59/yIvv\nqdACUVfdFr2OHTuSnp5+2Pf4wNeYk5NDdnZ2g9dotVrrnSs1NZXc3NxDthj5O87f77g539MT0YIz\nVGudprVO995/EPhJa90J+Ml7H2Ak0Mn7bzzwRlMH1lR/4Lm5uVx00UV06tSJDh06cO+991JdXe23\nbn5+Ppdddtlhzzlq1CiKioqOKp5JkybxwgsvHLZeaGjoIR8vKiri9ddfP6oYhBCNV1JSUu+CZbVa\nSU1NJTMz84i/MR/t59zBEpXk5OTDdp0dr+61xrQS+IszPT29QQJwYAxhYWEsW7bM19JzLC0QJyJZ\nOrBFDyAxMZGsrCwSEhKIjIxEa43WGrfbjdvtJiQkhLVr17J7zx5cLjc1Tie/LFxMjdOJ0+mmcM8e\nVq9dS1zr1rSKj2fLn38S1zqexKS2dD6jK6vXrmXTlq3k5ObidLupcbmpdrrZVeg5LiauNbGt4tm8\n9U9iW8UTHdeaVWvWsrNwDzVOd7O36iitddOdXKlsIF1rvbtO2UbgHK31DqVUPLBAa91FKfWW9/aM\nA+sd7Pzp6en6wJWM169fT9euXRsVX91vSQfOVDja/5Baa/r168cdd9zBTTfdhMvlYvz48URFRTFl\nypR6dZ1OJyZT04/znjRpEqGhoUyYMOGQ9UJDQyktLT3o49nZ2YwePZrMzMzjHeIp40j+voQ43mo/\now428NjfYNTs7GwyMzMZOHDgEX/Oud2aDVu2sjVrG/FJbWidkEThHgfrN20mKiaOHTsLSGjTFpMl\nhMoaF5XVLiqdLvbuK2V73g4sIWHsLdpHuDUKtyEAp9NNlctNtctNdY2baqeLaqemyu3C6dI43Rq3\nC2rcbtxujcutqXY5qaxyYTCA0+UNTHn+ud2gAbfef7v2YQ24Aa3doDUohVKe7/Ta7UI7q0EZMAQE\nesq0G3el5/PPGBTuew9cFSVoVw1GSxjK5JlB6a4qx11VijIHY7R4vhhqVzXOkr0ogxFT+P5dBJzF\nu9AuJ6aIOJTRc7yrzIG7shRDcCTGIM/F311dibN4FwazBVNEK9/x1YXbAE1ATFtf/M7iAtxVZZjC\nYzF4n99VWYqruACDJRRTRJzvddYUZoNSmOPa+85ZszcPd3UFAdYEDIHBnuNLHThL92AMjsQUHuOJ\nqaaKmt05KFMA5tjkOjFlo501mGPboUxmT0z7CnCVFWEKi8UY6vm7+mh0RJMMjldKLa/TaHJQTX11\n1cA8pZQG3tJavw20qk1avElOnLduIrC9zrG53rKDJjjHqilmKvz8889YLBZuuukmAIxGI//+979J\nSUnhySefZNasWcydO5fKykrKysp49913fUlDeXk5N954Ixs2bKBr165kZ2fz2muvkZ6e7tuWorS0\nlJEjRzJo0CCWLFlCYmIiX331FUFBQfz3v//l7bffprq6mo4dO/Lhhx8SHBx80FizsrK45pprcDqd\nnH/++b7y0tJSLrroIhwOBzU1NTzzzDNcdNFFPPjgg2zdupW0tDSGDRvGE0884beeEKJp1G0ZSU5O\n9n1Waa2pcrpwGQP5+Y81xLdNwWUKYnvBXtZt2Q6WGD7/+A+cJgsFjlLcJgvli5dTVumk0umixglO\nvBfu6nIMgaEYzBbAc+GrLszGGLLOd5HUzmpK134KShGWNtIXX8nKb6hx5BN25gUERMYDeyjf8gdl\n6xYQ1KEPobahgOcCvefblzGGRRNzwf2+4wtmP4mzuIC4yyf5koSijBmUrplHxIArCUvzfE5Vbs+k\n8PNnCEzqRtylj/uO3/7SlbhrKmlz7ye++HfPmUL5hkXEjJlISLchAJStX8ier18kuNsQYsdM9Lym\n6gpy/3MNyhxE279/6jtn4efPUJVrp9XV/8TStrvnda76lqIF7xHe9xKsQ2/2vKaiXeS/cwemqCQS\nb3vTd/yumY/iLNpBwm1vERCVCEDxb7MpWfYV1nNvJbzPWACqdmykYOYjBLbtQeurn/Mdv/PD+9E1\nVbT5+2yU9zU5FrzX4DVVbF3a8DXVVLHj/XsbvKY9371C1fbMeq+p1P5zg9fk2lfAzml/b/CaCr94\nFqej/mvat/Qrz2saegvhfS8GaPbB8U2d4AzUWud7k5gflFIbDlHX3wjOBs1LSqnxeLqwjsv0uIN9\nYBwtu91O796965WFh4fTtm1btmzZAsCvv/7KmjVriIqKqten/vrrr2O1WlmzZg2ZmZmkpaX5fY7N\nmzczY8YM/vvf/3LFFVfw2Wefcd1113HJJZdw2223AfDoo48ydepU/va3vx001nvvvZc77riD66+/\nntdee81XbrFY+OKLLwgPD2f37t3079+fCy+8kOeff57MzExWrVoFeFqg/NWTwbjidHQsU3kPPLbG\n5Wb7zt1sK3DgDook31FBzt5ysncVsTnfQYlLUfL9OqpqVoLRhDIYAaguyKK6MBtzXIrvG3f17hz2\n/T4bU2Q8kQOv9j5jJfnv3IGrooTE26f6WjH2fPsfyjcuJubCfxDSdTAAFdtWeS6cXYcQe6H3wuly\nsnfe6yhzUL0Ep2z9Qqq2ZxLcsa83wYGaPbmUr1+IMTQKvAmOdlZTuW01JmtCvfehZm8eTke+p3XF\ny11ZimtfIbq6fH9Ft6e1xV1dUe947XaB20W9S4cyeP5Rt8iICrD4WlS8pRgsoShzUL26hqAwDCGR\nYNh/uTSYgzCGxaAC63yBNJgwRbTCFBZd73hTRJzn6ub9HQEYQyIxRSX5Wl8ADAEWAqLbYgqPq3d8\nQExbdE11vSukKTyWgNhkDIEh+4+3hBIQl1Kv9QilCIhrjyHAUv+c1gR0dUW912oMicTcqgPG0Nr4\nNcpkxty6o5+Y2mEIDPG1aNXGZG7dyfNeeeXl5REZGdlsSU6TJjha63zvzwKl1BdAX2CXUiq+ThdV\ngbd6LtCmzuFJQP2h4p5zvQ28DZ4uqmON8cB+zWP9ZWit/V7g65YPGzaMqKioBnUWL17MvffeC0Bq\naio9evTw+xwpKSm+5Kd3796+JCkzM5NHH32UoqIiSktLGTFixCFjzcjI4LPPPgNg3LhxPPDAA75Y\nH374YRYuXIjBYCAvL49du3b5fU3+6rVu3fqQzyvEyexoE5XGLs5X5XSxu6SKLYUlbNlZypbCUtbl\nOvizoIx9VdU4S/eiDCZfN4GrrIiijBkoo4mo827znkWTP/VOanbnkHDrmwREJwFQZp/Pvj8+J3LI\njb4Ex11ZQlnmzwQmnAG+BMdzXndlCdpZBd4ExxAcjjHEWi8hMIVFY2nXE3NsO1+ZMpkJTRuJ8h7n\new/SRhLUoS+myP2fAcEd+mAKiyIgev/HuzE8hrgrn/ElVrViL3kU3K56F9SIs64kPP1CDJb94zgC\nk7qRdM/HKEP9S1ibe2aAUr5kwgC0unAiypuYed45iOw6mPCugzEAgSYwGICAYELvmwlAkAWMxgAU\nLqxXPYrFbMDtduOqAbMZTOeNxX32KGpqINAMZnMg2tqWlPumYjKD2wkBgQZqqt10u+dFzAEBOJ1O\nykorUApCL76ewMtvoaS8FFe1m9DwYMzt+mAYOJ2aaifV1ZU4a1wEB1uouOffhIWFEmD2vFaT0QCd\nJ9RrEVBKQZfzUaPPp7y0nLKyCkJDgwkIiKB44n+IiIjAbA6gsqKcXYW7aXXHg1iCgnE5ayguKiYi\nMoIA21Vw5VUooKKsHFOAiQBzHGrgxzirq9mzuwhLsJmIyEhcnV+kqNhBZIQVs9mT5Kju46mqvomi\nIgeRkVZsbazYbK2PedjHsWiyBEcpFQIYtNYl3tvDgaeAOcANwPPen195D5kD3K2Umgn0A4oPNf7m\neDiwLzoyMvKYfxk2m82XNNTat28f27dvp0OHDixfvpyQkBC/xzZ2PFRg4P4PBaPRSEWF51vMjTfe\nyJdffknPnj15//33WbBgwWHP5S8Zmz59OoWFhSxfvpyAgACSk5P9rv/S2HpCnEoam6gcmAhZrVaS\nkpJYvDiDDl26smrLdgqqzXz13RbWFZSxeec+qvX+b/Gla3+kKm89YWde4Ov6KVn9PY4f3yI0bSTR\nI+4CPGNDSlfOxRAcWSfBARVgAaOpXiuGuXVHgrsO8SU8AAFRSUSPug9jWEy9+ONvfAllNPtaEcxA\n+/PvxHLBnbSKMJEYF0VsqJnYv6QQdc9VhAYGEGYxERJowhJgxPT3czAaFEalyM/PY0deHm2vvoW2\nbdv4ypVSGBSU7Ctm48aNdOt6BlFWK8XFRWxYH0ebpETy8/JI8v5MTb2KoqIiX4u60Wj0m2zWzpA6\nkqEFB37e2+128vPzSU9P9x17sHMD9YYyJCUlkZub2+D+gX8zB9az2XrVO1dtq356etoBf2v96rwP\n3Rq9uOj+MVop5OXlYbVaSUjY/97k5OTgdqdgMBh8yXpjknfP33q3wyb9nnpd69Wz2WyUlJS0rAQH\naAV84b2AmoCPtdbfKaWWArOUUrcAOcDl3vrfAKOALUA5cFMTxgb4n6lwrL+M8847jwcffJBp06Zx\n/fXX43K5+L//+z9uvPHGQ46HARg0aBCzZs1i6NChrFu3jrVr1x7x64mPj6empobp06eTmJh4yPoD\nBw5k5syZXHfddUyfPt1XXlxcTFxcHAEBAcyfP59t2zw704eFhVFSUnLYekKcyho7Ni8sLIyVazIx\nx7Zj4x4nv2zIY+mfRZQUFeL+PAtzq/ZAGVq7yX/nTpxFO2n791m+QZkVWSspX/8LgUk2X4JjCovB\nGBaLoW7XQVA41vPGe1pW6mh97WSCjSZSYgOwJUbTJd5K9MWdKMg5h6T41nQ7ozMVZSWsWr6M7nfd\nTYeUdgQYDZgMiqKiogYX8oYX446N+hx0OBzsKttNP1t78vLyCCEOa0T946JD44gICsBut+P0nr9X\nj1SsVismg/IlNEC9FvWDJS4Oh+OIv5ge+Hlvs9lISEio93lvtVqxWq1kZWXVi6nu+ZVSLF68mEGD\nBpGcnExkZCQZGRmkpqb6vZb4GwJRd7p3bfwHJlRH2rNwsC/sCQn7uwH9JTG1r/lQGnvc0Z6/qTRZ\ngqO1/hPo6ad8D3Cen3IN3NVU8fjTFL8MpRRffPEFd955J08//TRut5tRo0bx3HPPHfbYO++8kxtu\nuIEePXrQq1cvevToQURERKOf++mnn6Zfv360a9eO7t2710tG/Hn55Ze55pprePnll7n00kt95dde\ney1jxowhPT2dtLQ0zjjjDACio6MZOHAgqampjBw5kgceeMBvPSFOJf66pAACAgJ8F6aSkhLKq51s\nK4FFm3bz47odbMgtpHrHZozBEd5kBipzN7Jr+kQC4tqTcNN/ADwzX9xOcDtxluwmwDvuJDT1XCxt\nbAQmdPE9Z1TnAZwzeBAp4Zpz0zrTPjGOILORYPNIggKMBJoMGAz7W133X9Q6+FolQiNNpNuSsVpD\nICaE+PCBlJSUEGze/3F/4IU8Ojq6XgtEY5OGI2kF9zfese4Qgf2tGemHPNfRfjFt7Of9gcMWysvL\n6z2fwWBg0KBBGAwG3zkGDhzY4PO2tr7dbq+XqEDD5KXu+1J7zLEmcMfjC/uprkmniTe1Y50mfrJx\nuVzU1NRgsVjYunUr5513Hps2bcJsNjd3aMLrVP77Ev75Wy5i6dKl7Klws8MQw7eZu9iws4TqHZs8\ngzhDPePn9i39EsfP7xCadj7RI+4GPANi8965ncD4zvVm9zj3FWIICqd7q0Au7NWOoIpd9OlpI6lV\nNNXlJaxYvhzYf3E/kqnch5s63pjjahd3q73A1j7euK6Lxo1XOjDOA7t1DtZldCL3WzqeS4f4O1ft\n9aru77m2rGPHjnW6lRJk36lDOFmmiYsjUF5eztChQ6mpqUFrzRtvvCHJjRBHqbEXX6vVSrdu3Zj/\nx2r+rArm02UF5O7agSmyNUp5Btfvnfc65RsWEXX+PYT1HA6AOb4TgYldCbDWdgVrDJZQLp38OdcP\n7ED3pGhUTTnFe3ZRXBRRJwFJBDxJhcVdQV5eHgkJCfUuakfy7ftoZoIerOUlIiKiQXfNoRxJq8iB\nz3dgt86huoxOlOPZCuLvXLXdRQeeKyEhgZSUFL/dSrX1T9dWmGMhCc5JpHZlTSHEsTvcYOHyKie/\n/bmH9zK2smhL7Sq0ZeyYNoHqHRuJv/lV30wkS9vuuEr2+BZFA0juZGPUX6Zw2ZBedGjTivKSfQ2+\n7Tsc1eQUFfntbqiblPgbRNrYi9rRzAQ90d0Z/p7vYN06zXkhP57DFvyd68CB6iUlJfVarKRb6fiS\nBEcI0SL5GyycmNKZr9c5+Ne839hbXk3RwmlU5a2n1TXP+9aSCYhKwFm0A+e+Ql+Ck9R3FHfcfytt\njPvo2r4N+fn5DboSAv1cnA6WSNQO8D3W5SmOdiboiR4MerINPj1ZyPvStCTBEUKc1I5lAT2r1UqI\nNY43flrH5+srKSv4ncBEz2B4ZTRRvvk3nI58qvI3YknqBkDUsDtICLUw+XIbZ9mSCTAa6iQSPXwX\noMZ0JRwsvqOZBeSPDCwV4uAkwRFCnNQauy5NXZU1Lr5fm8cjn6+l1OkZ/Jv7+o2gNUl/m+5bxt86\n9BYMgUHEtevEa9f34qxO8RgNypdABRg9M2WaemzG0Z5LWgCEODhJcIQQJzV/XU0Ha+3IzC3igdkr\nWb5sKdU7thDex7M3msESijk2Ga01rtI9GKISuTA1hmcnPUaYJaDBeRrTEnM8x2ZIUiLE8ScJThM4\n3K7cAIsWLeL2228nICCAX3/9laCgoEPWPxZNtQv4448/zuDBg/nLX/5yXM8rxIEOnClUOzjVarVS\nUe3ks2U5PDpnPeDZVXnXzEfA7Sb4jIGYvCv4trr6OcZ0Dua2/gn07NHdd26ZgitEyyQJTjOZPn06\nEyZM8O06fjhaa7TWvsWlwLNujtFoPMRRTcflcvHUU081y3OL08+BM4WSkpJYtGwNMze5mP/bCir+\nXEHk2dcCYLSEEtpjBAazBaWMDG0XyH9vG4rJZPR1bzkcjkZ3dwkhTk0tOsF58n921uXvO67n7JYQ\nzhNjGvdhuGDBAiZNmkRMTAyZmZn07t2bjz76iKlTpzJr1iy+//57fvzxR6ZPn86UKVOYNWsWVVVV\nXHzxxTz55JNkZ2czcuRIhg4dyq+//sqXX36JzWbj/vvv5/vvv+fFF18kKCiI+++/n9LSUmJiYnj/\n/feJj49n+fLl3HzzzQQHBzNo0CC/8e3YsYMrr7ySffv24XQ6eeONNzj77LOZN28eTzzxBFVVVXTo\n0IH33nuP0NBQkpOTufnmm5k3bx5333033333HaNHj+ayyy5j+fLlfuP4z3/+w5tvvonJZKJbt27M\nnDnzeP46xGmgwYJpLjOXvbOK3U5w11RR8Okk3FVlBHXqR2DrjgB0v/hO3h/Xg4LcbGw2GyaT54vA\nkXR3CSFObS06wTkZrFy50jfbYuDAgWRkZHDrrbeyePFiX3Iwb948Nm/ezB9//IHWmgsvvJCFCxfS\ntm1bNm7cyHvvvcfrr78OQFlZGampqTz11FPU1NQwZMgQvvrqK2JjY/nkk0945JFHePfdd7npppt4\n5ZVXGDJkCBMnTvQb28cff8yIESN45JFHcLlclJeXs3v3bp555hl+/PFHQkJCmDx5Mv/61794/HHP\nqqwWi4XFixcD8N133wFQU1PD3/72N79xPP/882RlZREYGEhRUdEJeMdFS1M7KDenDPo9+zl7N2QQ\nYjsHpQwYAgIJ6zMWXVOFMcTK3We1pn9UhSd58SY3jdkyQAjR8rToBKexLS1NqW/fviQleXb2TUtL\nIzs7u0GLyrx585g3bx69enl2mi0tLWXz5s20bduWdu3a0b9/f19do9Ho2zdq48aNZGZmMmzYMMDT\nbRQfH09xcTFFRUUMGTIEgHHjxvHtt982iK1Pnz7cfPPN1NTUMHbsWNLS0vjll19Yt24dAwcOBKC6\nupoBAwb4jrnyyisbnOdgcQD06NGDa6+9lrFjxzJ27NijeAfF6a48IJLhkxdRoTU7P/oH1bu2YrCE\nEtyxLwCxA6/mqbOMnH92P6zW+hsl+ktejmZhPCHEqadFJzgng8DAQN9to9GI0+lsUEdrzUMPPcRf\n//rXeuXZ2dmEhITUK7NYLL5xN1prbDYbv/76a706RUVFeHdxP6TBgwezcOFC5s6dy7hx45g4cSJW\nq5Vhw4YxY8YMv8ccGM+h4gCYO3cuCxcuZM6cOTz99NPY7XZMJvmzE4f35+5Sxr74I8UuA8pgRClF\nSLdzQCkMARasRphzdzpt4lv5urGSkpIOmbwc7cJ4QohTj+HwVURTGzFiBO+++65v5lVeXh4FBQWH\nPa5Lly4UFhb6EouamhrsdjuRkZFERET4upKmT5/u9/ht27YRFxfHbbfdxi233MKKFSvo378/GRkZ\nvp19y8vL2bRp01HF4Xa72b59O0OHDuX//b//R1FR0WFnlwlRVFbNkOe+of9t/2Td6+Mps8/3PRaW\nfiEj/vE2i5+9mp8nnEWb+FaAp9spKSmJzMxMbDYbKSkpvrE2DofDd/yh1qARQrQs8lX6JDB8+HDW\nr1/v6woKDQ3lo48+OuwMKbPZzOzZs7nnnnsoLi7G6XRy3333YbPZeO+993yDjEeMGOH3+AULFjBl\nyhQCAgIIDQ1l2rRpxMbG8v7773P11VdTVVUFwDPPPEPnzp2POI7OnTtz3XXXUVxcjNaav//970RG\nRh7luyRaOpfLzRNfreKjP3YA4K6pxLWvkPKNGYR2/wu2uCC+vGcIASb//y8MBgMDBw485AJ6sgaN\nEKcPpbVu7hiOWnp6uj5wc8r169fTtWvXZopItHTy93V81W7D8EduOTe/9h1Oxw6CUjxj0bTbRcWW\n30lJG8Dnf+2Hq6pC1qoRQqCUWq61Tj9cPWnBEUI0G6fBTN/JSygr3MaOD+7DYA4icfzbGCyhmAxG\n1r3/GEZXhaxVI4Q4YpLgCCFOOK01U75Zw+uLcgEIiGlLYEIXTOGxaLeLf1/UnljXHkr27JC1aoQQ\nR0USHCHECZVVWMqQZ+dSlDGDiLOuxBgUjlKKVlc8zV86hfHf287FYFCHne4thBCH0uQJjlLKCCwD\n8rTWo5VSKcBMIApYAYzTWlcrpQKBaUBvYA9wpdY6u6njE0KcOA/NWsqMFQXs+e5VyjcuRldXED3y\nHozA57fYKNmzi+Jiz4KQslaNEOJYnIgWnHuB9UC49/5k4N9a65lKqTeBW4A3vD8dWuuOSqmrvPUa\nrionhDjp1Q4erk1K9pRU0vvZn3yPRw4eh7uyhLA+Y7koUfPEtWcRFRWFwxFD7cSB9PR0WatGCHHU\nmnQdHKVUEnAB8I73vgLOBWZ7q3wA1C5ve5H3Pt7Hz1ONWa1OCHHSCQsL861BM/WXjaTe+SrFSz7x\nPR4QlcgZ1z/LSxe04e9jUomKigI8U7YTEhJISEiQtWqEEMekqRf6ewn4B+D23o8GirTWtcv55gKJ\n3tuJwHYA7+PF3vr1KKXGK6WWKaWWFRYWNmXsR00pxbhx43z3nU4nsbGxjB49ulni2bBhA2lpafTq\n1YutW7fWe2zUqFGH3SPqUMc3hdDQ0ON+zjfffJNp06Yd9/MK/6xWK53P6MrAyUuY9Nlyds18lKJF\nH1K5bQ0A9w2I4pXhEQw7ZyDJycn1jrXZbA1mTFmtVpkiLoQ4Ik2W4CilRgMFWuvldYv9VNWNeGx/\ngdZva63TtdbpsbGxxyHS4y8kJITMzEwqKioA+OGHH0hMTDzMUU3nyy+/5KKLLmLlypV06NCh3mPf\nfPPNYRffO9TxB+NyuQ55/0RyOp3cfvvtXH/99c0Ww+kmq6CEvi/8RjlgDI4gYtC1hPe9hKDErnxx\nfSd6hleSmJgoXU5CiCbTlC04A4ELlVLZeAYVn4unRSdSKVU79icJyPfezgXaAHgfjwD2HmsQSqkG\n+zKNGTMGpRT/+9//fGVvv/02SinGjx/vK8vPz0cpRUJCwhE/78iRI5k7dy4AM2bM4Oqrr/Y9VlZW\nxs0330yfPn3o1asXX331FeDvsI+kAAAgAElEQVTZe+rss8/mzDPP5Mwzz2TJkiWAZ8Xhc845h8su\nu4wzzjiDa6+9Fn8LNK5atYr+/fvTo0cPLr74YhwOB9988w0vvfQS77zzDkOHDm1wTHJyMrt37yY7\nO5uuXbty2223YbPZGD58OBUVFX6P/+ijj+jbty9paWn89a9/9SUvoaGhPP744/Tr149ff/2V5ORk\nnnrqKQYNGsSnn37K1q1bOf/88+nduzdnn302GzZsACArK4sBAwbQp08fHnvsMb/vZ1lZGRdccAE9\ne/YkNTWVTz7xdHcsX76cIUOG0Lt3b0aMGMGOHZ5VcM855xwefvhhhgwZwssvv8ykSZN44YUXAA4a\nx6effkpqaio9e/Zk8ODBR/LrPq3k5OTU2/4APHs85eTkAPD6j3YG3PsqVTu3+B6P6HcJ1/z1PpY9\nPJiyvQW+wcMHnkcIIY4brXWT/wPOAb723v4UuMp7+03gTu/tu4A3vbevAmYd7ry9e/fWB1q3bl29\n+3hageqVjR49WgN6zpw5vrK33npLA/q2227zleXl5WlAx8fHN3ieQwkJCdGrV6/Wl156qa6oqNA9\ne/bU8+fP1xdccIHWWuuHHnpIf/jhh1prrR0Oh+7UqZMuLS3VZWVluqKiQmut9aZNm3Tt65s/f74O\nDw/X27dv1y6XS/fv318vWrSowfN2795dL1iwQGut9WOPPabvvfderbXWTzzxhJ4yZYrfWNu1a6cL\nCwt1VlaWNhqNeuXKlVprrS+//HJfjHWPX7dunR49erSurq7WWmt9xx136A8++EBr7XmvP/nkk3rn\nnjx5su/+ueeeqzdt2qS11vq3337TQ4cO1VprPWbMGN85Xn31VR0SEtIgztmzZ+tbb73Vd7+oqEhX\nV1frAQMG6IKCAq211jNnztQ33XST1lrrIUOG6DvuuMNXv+5rOFgcqampOjc3V2vt+b34c+Df1+lo\n7969etGiRXrv3r317hfu3q0HPD1Xx13xlEYZtDE8TifdO1O3e+Br/fOqLXrevHl63rx5DY6rvS+E\nEI0BLNONyD2aYx2cB4CZSqlngJXAVG/5VOBDpdQWPC03Vx2PJ9N+WjrqttzUGj9+fL3WG4CEhAS/\nxzdGjx49yM7OZsaMGYwaNareY/PmzWPOnDm+FoXKykpycnJISEjg7rvvZtWqVRiNxnqbXPbt25ek\npCQA0tLSyM7OZtCgQb7Hi4uLKSoqYsiQIQDccMMNXH755UcUc0pKCmlpaQD07t2b7OzsBnV++ukn\nli9fTp8+fQCoqKggLi4O8OyWfumll9arf+WVnolwpaWlLFmypF5MtXtdZWRk8NlnnwEwbtw4Hnjg\ngQbP2717dyZMmMADDzzA6NGjOfvss8nMzCQzM5Nhw4YBnm6w+Pj4Bs9d16HiGDhwIDfeeCNXXHEF\nl1xyyaHeqtNa7aBfu91OYmIieXl5mIJCSZ/yGwCWNt0xt+6IJbkXweYglk0aTpglALu90nd83fPU\n3StKCCGOlxOS4GitFwALvLf/BPr6qVMJHNkV+SR34YUXMmHCBBYsWMCePXt85VprPvvsM7p06VKv\n/qRJk2jVqhWrV6/G7XZjsVh8jwUGBvpuG41GnE4nx9uBz1E7hqgurTU33HAD//znPxs8ZrFYGmwQ\nGhISAoDb7SYyMpJVq1b5fe7DTZjr3Lkzy5cv55tvvuGhhx5i+PDhXHzxxdhsNt8u5geqfe66DhXH\nm2++ye+//87cuXNJS0tj1apVREc3GOcu8CQniYmJZGdnUxMcy83vLMUYFo0yGFGmAFpfO5n02ADe\nvnUAYZYAAL9bLchGl0KIptLUs6hOazfffDOPP/443bt3r1c+YsQIXnnlFV/r0MqVKwFPK0x8fDwG\ng4EPP/zwiAbmRkREYLVaWbRoEQAffvihrzXneDrvvPOYPXs2BQUFAOzdu5dt27Yd9rjw8HBSUlL4\n9NNPAU+itHr1asDTcjJz5kwApk+f7vf4/Px8goODue6665gwYQIrVqygS5cuFBYW+hKcmpoa7Hb7\nUcexdetW+vXrx1NPPUVMTAzbt28/7Os6XTkcDvLy8vhtt4nrX/4fO96/B8fP7/gev+tMC+/cdhbR\n3unfQghxokmC04SSkpK49957G5Q/9thj1NTU0KNHD1JTU30Da++8804++OAD+vfvz6ZNm/y2QBzK\nBx98wMSJE+nRowerVq3i8ccfPy6vo65u3brxzDPPMHz4cHr06MGwYcN8A3sPZ/r06UydOpWePXti\ns9l8g6tffvllXnvtNfr06UNxcbHfY9euXesb2Pzss8/y6KOPYjabmT17Ng888AA9e/YkLS3NNzD7\naOKYOHEi3bt3JzU1lcGDB9OzZ89GviunF4fDQWZmJv9aVs5by/eBwYS7uhJn8S60y8nkQWYu799F\nWmaEEM1KHe0Yk5NBenq6rl31tNb69evp2rVrM0UkWjr5+4Ks7GzGvmunuHp/WdWOTYS16sBroyLp\n2K6NbJAphGgySqnlWuv0w9WTFhwhRKNVVNUw5JUVbP3sBap2bPaV9+5p4+0LIknv2Z2UlBTfIGSZ\nBi6EaC6ym7gQolEK91XS57mfKFn2FWXrFlC1cwsJt7zGbYNSGNcznPDwcJkhJYQ4abTIBEdrfdhZ\nOUIcqVO5O/do1W6aubNcM/IVz2Du8L4XU7M3n/C+F/PylWmM7e1/CwWZISWEaE4tLsGxWCzs2bOH\n6OhoSXLEcaO1Zs+ePfWm7p8OwsLCmL1gKZO+zcYU0QplNKGMAcRccB+f3JpOv46tmjtEIYTwq8Ul\nOElJSeTm5nKybsQpTl0Wi8W32OLpYmleGY9/vpaCT58gqPNZRI+8B6UUP/x9MJ1ahTV3eEIIcVAt\nLsEJCAggJSWlucMQ4pQ3fckWHpmzEa3daLcTd1UpuF38/uhwWkUENXd4QghxSC0uwRFCHLuXvrfz\n0vxsACxJ3Wh93RQCY9rx1vmRmN2VgCQ4QoiTm0wTF0LU8+QXK3j29WlU5a33lUW2ao/96VH06dVd\npn8LIU4JkuAIIXwe/3Q5r8/6jt1zplDw6SScpXuJCoQVT55PcKCp3vRvIYQ4mUkXlRACgIc/WcbH\nK3dhSbIR3HkA5lYdaJ8Yw6e39qFgZx5t23qmg8v0byHEqUASHCEE/5i5lE9W7kIphTKaiBn7IGe0\nCuLD63uxacN6vzuBCyHEyUy6qIQ4zd0/43emfvARe+b+C+327GDfJtjAqxem+JIbabERQpxqpAVH\niNNUTk4OL/2ynU+XbsPx8zvo6nKCuw6mZ+8+vH7ZGeRuzyE5OVmSGyHEKUkSHCFOU28tyedzexHG\n4AjiLn+C6h1b6GzrwxuXd2XnjnySk5PJy8sjMjJSkhwhxClHEhwhTkMvfrOGD3/fjiEwGABLko2O\nHW28P86T3NR2S0VGRmK326WbSghxymmyMThKKYtS6g+l1GqllF0p9aS3PEUp9btSarNS6hOllNlb\nHui9v8X7eHJTxSbE6eytn9cz+Z1Z5L11K1V5GwCIN8N7152BOcBUL5mRaeFCiFNVUw4yrgLO1Vr3\nBNKA85VS/YHJwL+11p0AB3CLt/4tgENr3RH4t7eeEOIo5OTkNFiMz+Fw8J+vf+ef8/6kbMNi3BX7\nqMhaQYwJpl1/Brt27iAsLKxBS43VavVNERdCiFNFkyU42qPUezfA+08D5wKzveUfAGO9ty/y3sf7\n+HlKtgMX4qiEhYXVW3HY4XDwzne/8a/FuwGIGX0/0aP/j7YDr+ab+/rTqWMHbDabrFIshGgxmnSa\nuFLKqJRaBRQAPwBbgSKttdNbJRdI9N5OBLYDeB8vBqL9nHO8UmqZUmqZ7BguhH+1XUt2u52srCw+\n+XkZLy3I9U0DVwYj8bahzJvQn7iY6HrHSHeUEKIlaNIER2vt0lqnAUlAX6Crv2ren/5aa3SDAq3f\n1lqna63TY2Njj1+wQrQwVquVxMREfl6xkafmbmTntPvZ/fWLaJeTYAW/PzWC1rExDY6R7ighREtw\n2ARHKdVKKTVVKfWt9343pdQthzuuLq11EbAA6A9EKqVqZ28lAfne27lAG+9zmIAIYO+RPI8QYj+H\nw8Ef6/7k6aUu3FVlaLcL7azGqBR/PDGcqrIScnJymjtMIYRoEo1pwXkf+B5I8N7fBNx3uIOUUrFK\nqUjv7SDgL8B6YD5wmbfaDcBX3ttzvPfxPv6z1rpBC44Q4vAcDgeLlq7mgUWVAFjapNL6uinEXvgP\nlj8xgpqKUux2O2FhYc0cqRBCNI3GJDgxWutZgBt842NcjTguHpivlFoDLAV+0Fp/DTwA3K+U2oJn\njM1Ub/2pQLS3/H7gwSN6JUIIn7xCB3d9nU/Vzi2+MnNsMs/2D6CoIF/WthFCtHiNWeivTCkVjXc8\njHeqd/HhDtJarwF6+Sn/E894nAPLK4HLGxGPEOIQKqqcjHpzBQWfTqJmby6trniKwMSuLJx4Dq59\nBWRnZ8sWDEKIFq8xCc79eLqPOiilMoBY9ncxCSFOIk6ni9QnvkcZTZis8bgrSzBGtOL7+84mzFCN\nPS9PtmAQQpwWDpvgaK1XKKWGAF3wzHTaqLWuafLIhBBHxO12c+ZT3+EClDGAmDETcJcVM2fiSOIC\nXfW6pWQLBiFES3fYBEcpdf0BRWcqpdBaT2uimIQQR+GcyT+Sv+JnQroNQRmMKGXggzuG0qtdNDk5\nOQfdgkESHCFES9SYLqo+dW5bgPOAFYAkOEKcJK5+YyFrvn6f4iUzqPhzObEXTuT/XdKNc22edTT9\nrW1jtVoluRFCtFiN6aL6W937SqkI4MMmi0gIcUQe/GQpv24rIaj9mZSs+paQbkOYMKIDV/RNae7Q\nhBCi2RzNSsblQKfjHYgQ4si9+eM6Zq4sACAwsSuJf/0vt189lmvTWskifkKI01pjxuD8j/1bJhiA\nbsCspgxKCHF4X6/czhNvz8YQGExgomcXlBHdE7lvaDvfAGIhhDhdNWYMzgt1bjuBbVrr3CaKRwjR\nCEuzdjP+1TkUfvEcKIi/4SU6JbXhoSGtWLduncyOEkKc9hozBueXExGIEOLwcnJy2F2luHzqGgKi\n2xDc9WxA0Touif9eewbbtm2TRfyEEIJDJDhKqRL87OaNZy0crbUOb7KohBB+VWHmkqnLAVAGI9Ej\n78Wk3cy4qRu7dubLIn5CCOF10EHGWuswrXW4n39hktwIceKVVFRz7gs/45j/LtrlWWtTKcX/bulG\nwc58bDYbKSkp2Gw27HY7DoejmSMWQojm05gxOAAopeLwrIMDgNZapmgIcYJUVbtInTSPwq8mU5Wz\nBu1yEvWX8XwxrhOhQWbiZRE/IYSo57DTxJVSFyqlNgNZwC9ANvBtE8clhPByuVykPv4dSimsQ2/G\nnNCFiP6XM+PqjpQ5CggLC2uQyFitVr+L+wkhxOmiMevgPA30BzZprVPwrGSc0aRRCXEaysnJadCt\ntHfvXs584jtqN38LbN2R1te9wBf3nMOAnl2kO0oIIQ6iMQlOjdZ6D2BQShm01vOBtCaOS4jTTlhY\nWL1kxeFwMPqlJWR99zaVOWt89T64NpX0zklA/e4oIYQQ+zUmwSlSSoUCi4DpSqmX8ayHI4Q4jmqT\nFbvdTlZWFuPeXMLmVb9QsnwOhV/8E3dVOe+OS6Nnm4h6qxRLd5QQQjR0qGnirwIzgIuACuA+4Fog\nAnjqhEQnxGnGarWSmJjIY7PXkVkMwWecTVjeBoLap/PCNen0SgiWVYqFEKIRDjWLajOeVYzjgU+A\nGVrrD05IVEKcphwOB6/OW8ei3RqlFMpgJGrY7VzRCfrE4EtuZHaUEEIc2qHWwXlZaz0AGALsBd5T\nSq1XSj2mlOp8uBMrpdoopeZ7j7Erpe71lkcppX5QSm32/rR6y5VS6j9KqS1KqTVKqTOP02sU4pTg\ncDh4Y85vTJu7gN1f/hPtrAbg+l7R3HFuN7Kzs0lMTJTkRgghGuGwY3C01tu01pO11r2Aa4BLgPWN\nOLcT+D+tdVc8s7DuUkp1Ax4EftJadwJ+8t4HGIlnl/JOwHjgjSN9MUKcyr5Zkc2bq8vY+9M7lG9a\nQvmmXxndJZwrU8PIy8vzrVIsM6aEEOLwGrMOToBSaoxSajqe9W82AZce7jit9Q6t9Qrv7RI8SVEi\nnjE9tV1dHwBjvbcvAqZpj9+ASKVU/JG+ICFORUs27uSRH/IxmINodeXTRA6+njFjxvD02FQcDoes\nUiyEEEfoUIOMhwFXAxcAfwAzgfFa67IjfRKlVDLQC/gdaKW13gGeJMi7QjJ4kp/tdQ7L9ZbtOOBc\n4/G08MjMEdEirMndy1VvZ2AI8CwUbo5LoU/Pbky9dRDbt2+vN+ZGVikWQojGOVQLzsPAr0BXrfUY\nrfX0o0xuQoHPgPu01vsOVdVPWYPNPrXWb2ut07XW6bGxsUcajhAnhdpF/f4s2MeIh98n761bqcha\nAUAHq5Ev/nYOSinatm0rqxQLIcRROGgLjtZ66LGeXCkVgCe5ma61/txbvEspFe9tvYkHCrzluUCb\nOocnAfnHGoMQJ6OwsDB++m0lE36poGz9QtxlRZStX0hK6pn8MHE4BkNjlqgSQghxMI3ebPNIKaUU\nMBVYr7X+V52H5gA3AM97f35Vp/xupdRMoB9QXNuVJURLU+4OYMIvFQBYz70Fc2w72vQ8j4yHR0py\nI4QQx0GTJTjAQGAcsFYptcpb9jCexGaWUuoWIAe43PvYN8AoYAtQDtzUhLEJccLk5OTU2xAz31FG\n34c+xhQeizKZUcpATI/hLHtmJEajJDdCCHE8NFmCo7VejP9xNeDZsPPA+hq4q6niEaK51O4xZbPZ\nqFJm0ie8z64ZD2OO70TcxY9gMZl5a1QY+/YVy8BhIYQ4TpqyBUcIwf6ZT4v/WMnf5leA3j923gLY\nnx1Jyb5iWaVYCCGOI0lwhDgBKjB7khs808BbXzeF0PBY1k2+CKPRINO/hRDiOJMER4gmlre3jPT7\n38FdVkRQh3QAIqMSWfNs/TE3VqtVkhshhDhOZESjEMdZ7Ro3ANmF++j78Ex2zXiYgi+epWrnFkKB\n171jboQQQjQNacER4jirHVRsaZXMZVNXY4poRUi3IbgrSkhq047fnxxJcXGxdEcJIUQTkgRHiOPM\narXiCovn0ndWoZRCKUXU8DvoEqH4/qELUEpJd5QQQjQx6aIS4jj7btU2LnvsTQo/ewrtrAbg7JQw\n5j08Gs/6l0IIIZqaJDhCHEf/nb+B8dOW4Zj/LhVbl1K2YRH9ouGVK3s2d2hCCHFakS4qIY6Tf3z8\nO7PW7MZgthB3+SQqczKZ+NcbGT+ojaxxI4QQJ5gkOEIcg9ptGG5493f+WJmJJakrAOa49rxx1xgu\n6p0MIGvcCCHECSYJjhDHICQkhD7PLSRv9pNU5qwl7vJJBCWnMe2aMxjcI9lXTwYVCyHEiSVjcIRo\npLrr2wBUVrvoPeU3nEYTATFtMVhCMASG8N0dvRnco0MzRiqEEEJacIRopLqbZlZipv/zP6OU5zuC\ndejNhPe9mO/uOYsz2rVu5kiFEEJIC44QjVS7X9QnP/xO9+seY9fHD+KuqQLAbDCy4P5BFBXuqtfK\nI4QQonlIgiPEQRzYJQXw4Ocree73MooyPqYqdx3lGxfTIwz+eOAsOnbogM1mw263S5IjhBDNTLqo\nhGD/bKi6A4HdbjcZGRkMHDiQiIgIOj/8LU7AEGAh7tLHqMpbz6Sbx3Dt0DTfcbIruBBCnBykBUcI\n9o+vqW15cTgc5Obmkpqayq/LVpJ4w79wrPrOV98cm8z81x4kLS6gwbmsVitt27Y9YbELIYRoSFpw\nhGB/y4vdbicxMZG8vDxsNhsLtzq464tt7JzxELjdmONSCEzowoanRmAxm3BEhUhrjRBCnISaLMFR\nSr0LjAYKtNap3rIo4BMgGcgGrtBaO5Rng56XgVFAOXCj1npFU8UmhD9Wq5XExESys7Np164dw19e\nQmElmMLjCO8zFu2sYWS/M3j3rlH1jpHkRgghTj5N2UX1PnD+AWUPAj9prTsBP3nvA4wEOnn/jQfe\naMK4hPDL4XCQl5dHeGxrznrkE/Lztvseixx8A1+/+wI39wyXAcRCCHEKaLIER2u9ENh7QPFFwAfe\n2x8AY+uUT9MevwGRSqn4popNnN78zY7Kzs4mIyODH7drzn/sQ3bNfITCOVPQzhoA1j81gsHdO/oG\nEAshhDi5nehBxq201jsAvD/jvOWJwPY69XK9ZQ0opcYrpZYppZYVFhY2abCiZfI3oHjZqkzuXqx4\nf/U+gpJ7YYpKILhTPy7tYib7+QsIDvQMJpYBxEIIcWo4WQYZKz9l2l9FrfXbwNsA6enpfusIUdeB\nU8CtVitJSUlkZGRgs9l47vPVfD7vB0J7jkAZjBgCg0m46RVWPjmKqJDAZo5eCCHE0TjRCc4upVS8\n1nqHtwuqwFueC7SpUy8JyD/BsYkWqu4WC1ar1TcFPLZ1AkPfWseuWc9QmbUcXVNFeN+LSYuFL+4f\ni2fsuxBCiFPRiU5w5gA3AM97f35Vp/xupdRMoB9QXNuVJcSx8jcFfMqSYlbvLQYgrPdonEX5mOM7\n8cdDQ4mLCG7miIUQQhyrppwmPgM4B4hRSuUCT+BJbGYppW4BcoDLvdW/wTNFfAueaeI3NVVc4vRU\nOwV87uJMJs34Ce12EpY2EoDgDn0YktaLt+8cjFWSGyGEaBGaLMHRWl99kIfO81NXA3c1VSxC7Cwo\nZOhb66jMXc+eb19GmYMI7tQfY4iVdU8Op6q8VBbsE0KIFuRkGWQsRJNwu910/r/ZOANDALAk2QhN\nO5/AJBsf3tSHc87sAkBwoCzYJ4QQLYnsRSVarLEvziWy71j+fP1GnPsKfOW33vcQm6c9QvuYoGaM\nTgghRFOSFhxxyvK3A7jdbueJr7NZ5l3Hz1W+D11dSUXWKtqmDWf5syMxGDx5vbTYCCFEyyUJjjhl\nHTj9e8jE91k5/wtCu/8Fc6v2AFiH3EDEgMtZNflK4mJimjliIYQQJ4okOOKkdLDWGQCbzQZ4WmAS\nEhJIn5yBC8XeBV9SsnwOrjIHsRc9AMC0K9vTq2cPaa0RQojTjCQ44qTkb3G+/HzP2o8JCQm4lIlu\nt0zBEBSGJakbAOF9LsZVXkREv0uZd2sX8rZvJ7ldW0luhBDiNCSDjMVJqe7ifFlZWdjtdtLT0/l1\nr4Fek5fQ5coHKfz8aYoXf+w7JjQ8hpyFn2F/8ToKdu4kOTmZvLw82f1bCCFOQ9KCI05atYvzLV+1\nhjs/XAqGeQR37AdASNfBlKz8hqAO6dx1ZhATrzjXs2nmsmUApKenY7VaiYyMrNcSJIQQ4vQgCY44\nKblcLjo88h0A5Rv/oPDL5wiIS/ElOAZLKB+9+Qp9z0wlOTkZ2D8mp/Z27U+bzSaL+AkhxGlGEhzR\n7GoHFIeHh9P2jqmUrPgaY4iViAGXA4qgDn0IbJNKcMe+DDG7eOnvZ/vG5djtdiIiInzJS+0A5Lqs\nVlnETwghTjeS4IhG8TeryeFwUFJSQtu2bRtd50Br1m9i1CsZmMLjAHBX7KNk+f8whsUQ3v8ylFKE\nmQLIyl7N+vXrAWmdEUIIcXiS4IhG8Terqfb+kdTJysri/LfWUQWUrfuF3V+/SHDXs4kdMxEAc0IX\nIgZeQ1T73mx6diQm0/4/UWmdEUII0ViS4IhGqTurKTExkby8vAYDd/3VCQ8Pp9/kJRRvXca+ZV8R\n3GUgYWnnA2Bu1cFzoLMGgD6hcOfAUJJvf4aUlJQT/hqFEEK0HJLgiEarndWUnZ1NcnKy35aTRWv/\n5LbX5lOx5Q9Cuv+FwNYdAXCVOajMXokKCPQlOKaoRKZ9NIPRI4cB1EuMIiMjpWVGCCHEUZME5zRz\nLGNp8vPzcTgcJCcnY7fbufGtddgzf6Z6x0Yizx6HwRIKQMWW3ylZMRdjWLQvwbGk9CJmzEQeG9mD\ncVedRVFRkTdR8nQ71e3akqndQgghjpUkOM3oaAblHuu5DhwnY7fbyc/PJz093VfH7XaTkZHBWWed\nRWhoKJ0f/wF3dQVFv3yAq7LEO15GAZqSlXOpzt9IcOeBWNr1ACCo0wCMIVGEtu/NByND6dGjhzeh\naeub0p2Xl+dbiK+8vLxeMiODh4UQQhwrSXAOozF7IsHRJSaNGZR7LOfKyMggNTW1Qd3axCYxMZG1\na9fyxJwsSj/LISCyNQCV2zPZ881LmOPaE3vxwwAok5mS1d+By4l7+F0YAoMBRVjPEbg6n4UpshVz\nrmlDmzZtsFovqPe+1G352bJlC9BwIb7a9WvqxijJjRBCiKMlCc5hHG5PpIMlJo3dLDIpKYmMjAxs\nNpvfgbuNVXeAb1xcHPPnz2fKwiKcbRTgeV7Hgvep3rWVqBF3eZOZYhwLZrPv98+IHHw9EQOuAEAF\nWHAW7UQFWHznVwYj0cPvxBAcwZIJA0lK9CQkWVndfGNyDjYw2OFw+F5XeXm57/2rG7e01gghhDie\nTqoERyl1PvAyYATe0Vo/38wh+Z0ZVNudc6gZRY1NjHJzc2nbtm2Dgbt2u51Vq1bx6uowdhiMAJRv\n+YPqnVsI7jwAc5wnmajYuhTHLx8QlNwL67m3AOAqzyH3ldsxBIbQ5r5PfDFV5a2nKteOq7jA11oT\nENOOwKRuGIIjffWCYtqxaNEiqqqqaN++fZ3Xt79lBjyJS92uJn8Dg0tKSuq9NzabjYSEhHoJjbTW\nCCGEON5OmgRHKWUEXgOGAbnAUqXUHK31uuaIZ8uWLezZs4fOnTvz2WdLmL4FNhdvonrXn5giVvim\nOLurdlC69keUyeybHQRQ/NunOIt3ETHgKkzhMQCUrJ5HWeZPhPYYTmj38wCoyt/Irk8exRzXntbX\nTgY8L3f7y1fhriwl6Z4ZGIPCACjfmEFZ5k+YwmN9CY7+/+3dfWydZRnH8e9vK4WVvXTDDtmK0tFN\n3HhtZoGVGOR1IDojmI9o4TEAAAkySURBVGAIomKIibxoTAzGRPA/SVDESEgIoGgMKIMo8IdAJtEY\nI6yggRYQOjahG7IZDrNsk2708o9zn7PT07OtZT3PaZ/9PsmTc577udve5zp3zq7dL+d5fze7t22i\nad7C8t+eMWsOhy1cwsxZc4kYQSreU7X1rCvQnmF+9JkOLrpoVUqyVtHbW2x3Z2cnmzdvpr29ncHB\nQbq6uva56Ldy1Gp/C4NrTdk5oTEzs3qbMgkO0A0MRMRrAJIeANZQ+hc/A6VppdNu+StbH7yZXa/1\n0nbp92np7AZg18bnePvxO5h96mqOuvBaAEaGd1JYdxczZy8YleDsfPkvDL+1gdmnrC4nOO8PbeO9\nwX6O+MhJ5Xqa2UQM72JkeOeotjQffTwjw/+DkT3lslmd3TTNbeOwhXungpZ95GTuffRRduzYwYIF\nzeW1LYUbz6K/v7+crBRHmuaV79dUnWAsWrSIjo4OWltby2t39rfot3pkxlNNZmY2lUylBGcx8EbF\n+SBwenUlSdcA10Dt0YGDUZpWguJ3tDTv2s6M5lnl603zjmZW5+k0tx1XLpvR3MKcrkvKW6RL5nZ/\nnpH3dtA056hy2TEnnsvlJ57M+ecvpatrOf39/ZxwwsUsvP3LtLS00NTUtHdR7iU/GzP9VWu9y8aN\nG9m0aRNtbW0sW7ZsVMLR3t5OX18fPT09+1zQOzQ0VE6KSj/X09PD0NDQqNdTPerikRkzM5vKFBGN\nbgMAkr4AXBgRX0vnVwLdEXHdvn5m5cqV0dvbO6ntKE29lJKLWiMe49lFVT2FUygUKLW1PMpSVadS\nKXEpJTPV7Sr97f2tA5rMbehmZmZTgaRnI2LlgepNpRGcQeDYivN2YMs+6tZN9bf11toZVL2Nu5RI\nVCYNW7ZsGTWiUUqUSs9Lj7WmdaoX70picHBw1HqX6mSp1hoYj7KYmdmhakajG1BhPbBUUoekZuBy\n4JGsG1GdXBQKhQP+TGlqq1S3UChQKBTGfLfLihUrxiRH8+fPH5WIVI7qdHR0sGLFCvr6+mhvbx+T\nLFWOLFUmS2ZmZoe6KTOCExF7JF0LPE5xm/i9EdGfZRvGuzOo2nhuRDletRbv1loT4ztrm5mZ7duU\nWYPzQUz2GpyDXbNSvW7GzMzMJtd0XIPTcAezZmU8X3pnZmZm2ZhKa3CmrVrrZirX5JiZmVm2nOBM\ngv196Z2ZmZllz1NUk8Dbsc3MzKYWj+CYmZlZ7jjBMTMzs9yZ1tvEJW0D/lWnX/8h4D91+t1Wm2Oe\nPcc8e4559hzz7NUz5h+NiLYDVZrWCU49Seodzz57mzyOefYc8+w55tlzzLM3FWLuKSozMzPLHSc4\nZmZmljtOcPbtrkY34BDkmGfPMc+eY549xzx7DY+51+CYmZlZ7ngEx8zMzHLHCY6ZmZnljhOcGiSt\nlvRPSQOSbmx0e/JC0rGSnpL0kqR+STek8gWSnpT0anqcn8ol6afpfXheUldjX8H0JGmmpL9Leiyd\nd0h6OsX7N5KaU/nh6XwgXT+uke2eriS1Slor6eXU1890H68vSd9Knyl9ku6XdIT7+eSSdK+krZL6\nKsom3K8lXZXqvyrpqnq22QlOFUkzgTuAi4DlwBclLW9sq3JjD/DtiPg4cAbwjRTbG4F1EbEUWJfO\nofgeLE3HNcCd2Tc5F24AXqo4vwW4LcW7AFydyq8GChHRCdyW6tnE3Q78ISJOAE6hGHv38TqRtBi4\nHlgZEScCM4HLcT+fbL8AVleVTahfS1oA3AScDnQDN5WSonpwgjNWNzAQEa9FxDDwALCmwW3KhYh4\nMyKeS8+HKH7wL6YY3/tStfuAz6Xna4BfRtHfgFZJx2Tc7GlNUjvwaeDudC7gHGBtqlId79L7sBY4\nN9W3cZI0F/gkcA9ARAxHxDu4j9dbEzBLUhPQAryJ+/mkiog/A29XFU+0X18IPBkRb0dEAXiSsUnT\npHGCM9Zi4I2K88FUZpMoDQufBjwNHB0Rb0IxCQIWpmp+Lw7eT4DvACPp/CjgnYjYk84rY1qOd7q+\nPdW38VsCbAN+nqYF75Z0JO7jdRMRm4FbgdcpJjbbgWdxP8/CRPt1pv3dCc5YtTJ576WfRJJmAw8B\n34yI/+6vao0yvxfjJOkSYGtEPFtZXKNqjOOajU8T0AXcGRGnATvYO2xfi2N+kNIUxxqgA1gEHElx\niqSa+3l29hXjTGPvBGesQeDYivN2YEuD2pI7kg6jmNz8OiIeTsVvlYbl0+PWVO734uD0AJ+VtIni\nVOs5FEd0WtNQPoyOaTne6fo8xg5J2/4NAoMR8XQ6X0sx4XEfr5/zgI0RsS0idgMPA6twP8/CRPt1\npv3dCc5Y64GlaQV+M8XFao80uE25kOa57wFeiogfV1x6BCitpr8K+H1F+ZfSivwzgO2l4VA7sIj4\nbkS0R8RxFPvxHyPiCuAp4LJUrTrepffhslTf/7OdgIj4N/CGpI+lonOBF3Efr6fXgTMktaTPmFLM\n3c/rb6L9+nHgAknz08jbBamsPiLCR9UBXAy8AmwAvtfo9uTlAM6iOBz5PPCPdFxMcf57HfBqelyQ\n6ovijrYNwAsUd0k0/HVMxwM4G3gsPV8CPAMMAA8Ch6fyI9L5QLq+pNHtno4HcCrQm/r574D57uN1\nj/kPgJeBPuBXwOHu55Me4/sprnHaTXEk5uoP0q+Br6bYDwBfqWebfasGMzMzyx1PUZmZmVnuOMEx\nMzOz3HGCY2ZmZrnjBMfMzMxyxwmOmZmZ5U7TgauYmdWXpNJ2U4APA+9TvOUBwM6IWNWQhpnZtOVt\n4mY2pUi6GXg3Im5tdFvMbPryFJWZTWmS3k2PZ0v6k6TfSnpF0g8lXSHpGUkvSDo+1WuT9JCk9eno\naewrMLNGcIJjZtPJKcANwEnAlcCyiOgG7gauS3VuB26LiE8Al6ZrZnaI8RocM5tO1ke6V5OkDcAT\nqfwF4FPp+XnA8uJtiQCYK2lORAxl2lIzaygnOGY2nbxX8Xyk4nyEvZ9nM4AzI2JXlg0zs6nFU1Rm\nljdPANeWTiSd2sC2mFmDOMExs7y5Hlgp6XlJLwJfb3SDzCx73iZuZmZmueMRHDMzM8sdJzhmZmaW\nO05wzMzMLHec4JiZmVnuOMExMzOz3HGCY2ZmZrnjBMfMzMxy5/+uXyhTiAsMwAAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x26df3d27160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pints.plot.series(samples[:100], problem)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Marginal likelihood estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "marginal log-likelihood = -375.348895962891 ± 0.07908966826933428\n"
     ]
    }
   ],
   "source": [
    "print('marginal log-likelihood = ' + str(sampler.marginal_log_likelihood())\n",
    "      + ' ± ' + str(sampler.marginal_log_likelihood_standard_deviation()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Effective sample size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "effective sample size = 1597.4248390822784\n"
     ]
    }
   ],
   "source": [
    "print('effective sample size = ' + str(sampler.effective_sample_size()))"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}