sshojiro/rudimentary-mcmc-with-pymc.ipynb

## rudimentary-mcmc-with-pymc.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian modeling with `pymc` \n",
    "\n",
    "Pymc example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('data/txtdata.csv', <http.client.HTTPMessage at 0x20642e00dd8>)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from os import makedirs\n",
    "makedirs(\"data\",exist_ok=True)\n",
    "from urllib.request import urlretrieve\n",
    "urlretrieve(\"https://git.io/vXTVC\",\"data/txtdata.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.core.pylabtools import figsize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "figsize(12.5,3.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "count_data=np.loadtxt(\"data/txtdata.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 74)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xm8X/Odx/HXO0ERS5DLxBJRUqVKEFuZjrVFtbRFGTVpizDTVk1X2jFDp6a0qtWplpRW1L5TXTBBS1WIraSxi4pcSSyp2CM+88f3+5Of6y7nJvf8lnvez8fjPn6/s39+557f+fy+3+8536OIwMzMqmtIswMwM7PmciIwM6s4JwIzs4pzIjAzqzgnAjOzinMiMDOrOCeCAiSdLunYXqaHpA36sb5+zW/FSHpR0rtLWO9xks5dzGU/I+mWXqb/TtL4xY+uUAwzJO1a5jYGG0nTJO3Y7DgapfKJIH9JXpE0X9I8SbdKOkLSW/smIo6IiP9ezPXfJOnQgYu4OfLn2DG/P07ScYuxjgFLgN3t14hYISIeG4j1N0pE7BERk6DvpGHlkHS2pO/Uj4uI90XETSVsa19Jd0t6QdI9krYe6G0sjsonguyjEbEisC5wIvAN4KzmhtQalPg4sUFB0lJNDmEkcBiwCnA+cHZTo6mJiEr/ATOAXbuM2xp4E9gkD58NfKdu+teATmAW8DkggA26WfcJwELgVeBF4Cd5fABHAA8DzwOnAapb7nPA9DztWmDdHmLfEZjZ0+fJn2Mq8AIwGzilbr5tgVuBecC9wI51027Ksf8JeAXYII/bMU8/Djguvx8BXJPX8xxwMzCkm1j/mD/3S3lffCqP3wu4Jy9/K7BpHr9+Xt8WeXhN4Jn8mXvbrxvU/c9OA34DzAemAOvXxfMh4EHg78BPgT8Ah/awn48DLgbOyeuaBoyrm3408Gie9lfg43XTPpP34//mbT0A7NJlXx8KbJQ/z8L8mebl6Xvmdc4HngK+2suxfFg+bmpx1PbdDOCrwF9yDBcBy+Zpq+T/31zS8XYNsHaX+P47f4b5wHXAiLrp/wI8ATwLHMvbj78hdfvm2bwPV+0j/kfy//1qYM08/nTg5C7zXgV8ue7YuCx/hseBI7v87y4FziV9Dw7tsp4JwALg9bzff93N9+g44JK8jvnAfcB7gGOAOcCTwIfq1rky6YdkZ/6ffQcY2s3n3Q54ttnnwIhwIqCbRJDH/w341/z+bHIiAHYnnVQ3AYaRsnq3iSDPf1M3B1/kL9xwYFQ+gHfP0/bJX4aNgKWA/wBu7WHdO9J7IvgzcHB+vwKwbX6/Vv5i7pm/rLvl4Y66mP8GvC/HsHQv+++7+Yu6dP77R+qSWjefe4O64S3yF2kbYCgwPsf/rjy9dmJbnpQQTy6wX+sTwXOkZLgUcB5wYZ42gnRS+ESe9iXSyaC3RPBq3l9D82e+rW76fqST0RDgU6RkNzJP+wzwBvDvef98inQyXrXr58jz3tJl253AP+b3q5BP7t3EuB/ppLMVIFLyXrfumLg9x7hq3qdH5GmrAZ/M+3hF0gnvyi77+VHSiW+5PHxinrYx6eS5A7AMcHLej7Xj7yjgNmBt4F3AGcAFPcS/MynRb5Hn/V/gj3naB0knW9Xth1fq9vmdwH/mGN4NPAZ8uO5/t4D0vRoCLNfNts+m7odeN9+j2v//w6Tj5RxSwvlW/p8eBjxet+yV+bMOA1bP+/7wLusfRvp+/qDZ58CIcNVQL2aRvjRd7Q/8MiLuj4iXSAfJ4jgxIuZFxN+AG4GxefzhwHcjYnpEvAH8DzBW0rqLsY0FwAaSRkTEixFxWx7/aeC3EfHbiHgzIq4nlRz2rFv27IiYFhFvRMSCPrYxknTSWRARN0c+0gs4DDgjIqZExMJIdeWvkUorRMTPSaWmKXkb3yq43prLI+L2vB/PY9E+3hOYFhGX52k/Bp7uY1235P21EPgVsFltQkRcEhGz8r68KMdcX/c7B/hR3j8XkUoiHyn4GRYAG0taKSKej4i7epjvUOB7EXFHJI9ExBN103+cY3wO+DV5X0TEsxFxWUS8HBHzSaWtf+qy7l9GxEMR8QrpV31tP+5L+gV9S0S8TjoZ1//vDwe+FREzI+I10ndl3x6qZw4CfhERd+V5jwG2kzSaVMoM0o+M2nb/HBGzSImvIyK+HRGvR2oj+jlwQN26/xwRV+b/zys97L++3BwR1+bj5RKgg/QdXgBcCIyWNFzSGsAewFER8VJEzAF+2CUeSPtxFqkauumcCHq2FukXZVdrkn6d1DzRzTxF1J94Xib9YofUTnFqbriuVbcox9Nfh5B+yT0g6Q5Je9VtY7/aNvJ2diCdbGuepJjvk0ow10l6TNLR/YhvXeArXeJYh7SPa35OKn39bz5B9EdP+/ht/8OcuGb2c13L1k5okv4lN/zVPsMmpFJHzVNdkuMTvP0z9uaTpMT1hKQ/SNquh/nWIf1yLxr/Cjn25SWdIekJSS+QqvCGSxra17K8cz++TCpZ1qwLXFG3X6aTqr7W6Ca+Nan7LkXEi3lda+V9dyFwYJ78z6TEXtvGml2OoW922UbRY7k3s+vevwI8k38U1IYh7Zd1SaWEzrp4ziCVDACQ9B5Ssj04J5ama3bDSUuStBXpxNvdFRydpC9dzag+Vtff7l2fBE6IiPP6nDNVQSxfG8hf3o63NhzxMHBgbuz9BHCppNXyNn4VEYctadz5V+RXSCf09wE3SrojIiYXWLz2WU/obqKkFYAfkepbj5N0Wf5FWzi+HnSSqitq21H9cH/kktrPgV1IvzwXSrqHlLxr1pKkumQwilQH3tU7PlNE3AHsLWlp4AukX5LrdJ2PtC/XX4yP8BVgQ2CbiHha0ljg7i7x96QzLwuApOVIVU31MX0uIv5UYF2zSCfR2rqG5XU9lUddQPqxcSKpKvHjddt4PCLG9LLuvo6VgeyC+UlSqXZELyf5kaQ2oJcHcLtLxCWCOpJWyr+aLwTOjYj7upntYuAzkjaWtDzwX32sdjap3rKo04Fj8kkVSStL2q+HeR8i/TL9SD5R/AepfrX2eT4tqSMi3iQ1xkL6RXYu8FFJH5Y0VNKyknaU1O+ToaS9JG2QT6Yv5PUv7GH2rvvi58ARkrbJVycNy59lxTz9VODOiDiU1Oh7ei/r6o/fAO+XtE/+Vf954B8Wc13DSCeSuQCSPksqEdRbHThS0tL5f7kR8Ntu1jUbWFvSMnldy0g6SNLKuQqitn+7cybwVUlb5n25QcHqxBVJv2jnSVqVvo/nepeSjqMP5JiP5+0J5HTghFockjok7d3Dus4HPitprKR3kapEp0TEDICIuJu0j88Ero2I2vF8O/CCpG9IWi4fz5vkH3NFLcmx9DYR0UlqUP9BPp8MkbS+pPrqtj8Dmw/E9gaKE0Hya0nzSdn8W8ApwGe7mzEifkf6lXoDqUrkhj7WfSqpXvR5ST/uK5CIuAI4CbgwF9XvJ9U5djfv34F/I305niKVEOqrOHYHpkl6McdxQES8GhFPAnuTitBz8+f+Got3PIwB/o/UaPhn4KfR8/XXxwGTcpF5/4iYSmon+AnpipVHSA2m5BPG7qSrqwC+DGwh6aA83K/9Wi8iniE1rn6PVP2wMamNpL9VT0TEX4EfkD77bOD9pCts6k0h7adnSHXw+0bEs7zTDaQrkp6W9EwedzAwIx8LR5Dad7qL45K87vNJV7ZcSfdtXF39iNQI/AypYff3BZapbXMa8EXSD6fOvN05LNqPp5JKPtfl79dtpF/z3a1rMumqo8vyutbnnfXqFwC7kj5jbbmFwEdJ7RaP589xJunKnaLOIrXDzJN0ZT+W68m/kBqu/0o6ri/l7dWu29B9bUPTKAq365kNTrnqbCZwUETc2Ox42lWuypsHjImIx5sdjxXnEoFVUq4WG56rIb5JqtK4rY/FrAtJH80NzsNIl4/eR7r00tqIE4FV1Xakq2yeIVUt7LMElxZW2d6kht5ZpOqvA/px+bC1CFcNmZlVnEsEZmYV1xb3EYwYMSJGjx7d7DDMzNrKnXfe+UxEdPQ1X1skgtGjRzN16tRmh2Fm1lYkFer5wFVDZmYV50RgZlZxTgRmZhXnRGBmVnFOBGZmFedEYGZWcU4EZmYV50RgZlZxTgRmZhXXFncWW2va6vtvfybPHV/buUmRmNmScInAzKziSksEkjaUdE/d3wuSjpK0qqTrJT2cX1cpKwYzM+tbaYkgIh6MiLERMRbYEngZuAI4GpgcEWOAyXnYzMyapFFVQ7sAj0bEE6QnGk3K4ycB+zQoBjMz60ajGosPAC7I79eIiE6AiOiUtHp3C0iaAEwAGDVqVEOCNGtn9Y33bri3/ii9RCBpGeBjwCX9WS4iJkbEuIgY19HR53MVzMxsMTWiamgP4K6ImJ2HZ0saCZBf5zQgBjMz60EjEsGBLKoWArgaGJ/fjweuakAMZmbWg1ITgaTlgd2Ay+tGnwjsJunhPO3EMmMwM7PeldpYHBEvA6t1Gfcs6SoiMzNrAb6z2Mys4pwIzMwqzonAzKzinAjMzCrOicDMrOKcCMzMKs6JwMys4pwIzMwqzonAzKzinAjMzCrOicDMrOKcCMzMKs6JwMys4pwIzMwqzonAzKzinAjMzCrOicDMrOKcCMzMKs6JwMys4pwIzMwqrtREIGm4pEslPSBpuqTtJK0q6XpJD+fXVcqMwczMeld2ieBU4PcR8V5gM2A6cDQwOSLGAJPzsJmZNUlpiUDSSsAHgbMAIuL1iJgH7A1MyrNNAvYpKwYzM+tbmSWCdwNzgV9KulvSmZKGAWtERCdAfl29u4UlTZA0VdLUuXPnlhimmVm1lZkIlgK2AH4WEZsDL9GPaqCImBgR4yJiXEdHR1kxmplVXpmJYCYwMyKm5OFLSYlhtqSRAPl1TokxmJlZH0pLBBHxNPCkpA3zqF2AvwJXA+PzuPHAVWXFYGZmfVuq5PV/EThP0jLAY8BnScnnYkmHAH8D9is5BjMz60WpiSAi7gHGdTNplzK3a2ZmxfWYCCQd2duCEfHjgQ/HzMwarbcSQe1SnTHA1sCv8/BewB/KDMrMzBqnx0QQEccCSLoWGBsRL+ThY4GLGhOemZmVrUgbwbrAq3XDrwHrlRNOdW31/Rveen/H13ZuYiRmVjVFEsH5wBRJlwEBfAI4r9SozMysYfpMBBHxbUm/I/UbBHBERNxRblhmZtYoRW8oGwrMjYgfAI9JGlViTGZm1kB9lggk/QewPbA+cA6wLKm6aIdyQzMzs0YoUiLYF9iT1GkcEfEUsFKZQZmZWeMUSQSvRUSQGoqRtHy5IZmZWSMVSQSXSzoNWFnSZ4HrgF+UG5aZmTVKkauGTpK0B/A66XGTJ0TE70qPzMzMGqJIY/EXgUt88jczG5yKVA2tDtwk6UZJh0saUXZQZmbWOEWqho4FjpW0BfAp4FZJj0XE7qVHZ2alqu/aBNy9SVX15wllTwIzgE7AN5SZmQ0SfSYCSYdJ+j/gZmBt4IsRsXHpkZmZWUMU6XRuQ+DoiJhadjBmZtZ4vT2hbFhEvAR8Ow+/7W7i2vMJzMysvfVWIrgU2AOYRrqrWHXTArcTmJkNCr09oWyP/LrO4q5c0gxgPrAQeCMixklalfSEs9Gkxuf9I+L5xd2GmZktmUJXDUk6QNI38/u1JW3Zj23sFBFjI2JcHj4amBwRY4DJedjMzJqkyFVDPwF2Ag7Oo14GTl+Cbe4NTMrvJwH7LMG6zMxsCRUpEXwgIg4nP7c4Ip4Dlim4/gCuk3SnpAl53BoR0ZnX1Um6c/kdJE2QNFXS1Llz5xbcnJmZ9VeRy0cXSBrCom6oVwPeLLj+7SNilqTVgeslPVA0sIiYCEwEGDduXBRdzszM+qdIieA04DKgQ9LxwC3ASUVWHhGz8usc4Apga2C2pJEA+XXOYsRtZmYDpEhfQ+dIuhPYlXQJ6X4RcX9fy0kaBgyJiPn5/YdI9yRcDYwHTsyvVy1B/GbWA/cjZEUV6YZ6K2B6RJyah1eUNK7AncZrAFdIqm3n/Ij4vaQ7gIslHQL8DdhviT6BmZktkSJtBBOB+stFXwLO6DLuHSLiMdKDbLqOfxbYpR8xmplZiYq0EQyJiLcah/P7pcsLyczMGqlIInhc0r9KGippiKTPk+4INjOzQaBIIjicVJUzO//9E3BYmUGZmVnjFLlqaDawbwNiMTOzJijSxcQGkq6VdG8e3lTSMeWHZmZmjVCkauhM4HgW3U18H/Dp0iIyM7OGKpIIhkXErbWBiAhgQXkhmZlZIxVJBM9KWo9FfQ3tAzxdalRmZtYwRW4o+wJwFvBeSU8AncABpUZlZmYN02sikDQU2Cwidpa0MqCImNeY0MzMrBF6rRqKiIXAUfn9350EzMwGnyJtBNdKOkrSSEkr1f5Kj8zMzBqiSBvB4fn1K6QGY+XXUWUFZWZmjVPkzuJ1GhGImZk1R5GqITMzG8ScCMzMKs6JwMys4op0OretpOXz+wMlfU+S2w3MzAaJIiWCicArkjYFvkl6JsG5pUZlZmYNUyQRvJE7mtsbODUifgCsWG5YZmbWKEUSwUuSvgYcDPxG0hD68czi/IjLuyVdk4fXkzRF0sOSLpK0zOKFbmZmA6FIIvgU6SaywyOiE1gbOKUf2/gSML1u+CTghxExBngeOKQf6zIzswHWZyKIiFnA+XWj5gAXF1m5pLWBj5AeboMkATsDl+ZZJgH79CNeMzMbYEWuGvoccDX5ZE7qWuKqguv/EfB1Fj3dbDVgXkS8kYdnAmv1sN0JkqZKmjp37tyCmzMzs/4qUjV0JLAt8AJARDwErNHXQpL2AuZExJ31o7uZNbpbPiImRsS4iBjX0dFRIEwzM1scRTqdezUiXk+1Om89o6CI7YGPSdoTWBZYiVRCGC5pqVwqWBuY1f+wzcxsoBRJBH+S9HVgWUk7AZ8HrulroYg4BjgGQNKOwFcj4iBJlwD7AhcC4ylezWRm1ra2+v4Nb72/42s7NzGSdypSNfR1YD7wAOkKoMnAt5Zgm98AvizpEVKbwVlLsC4zM1tCRbqhXgj8LP8tloi4Cbgpv38M2Hpx12VmZgOrz0Qg6W7e2aD7d2Aq8N2IeK6MwMzMrDGKtBFcT7rap3YvwQHAQuBF4GzgY6VEZmZmDVEkEXwgInaoG75b0i0RsYOk+8oKzMxsMKtvPIbmNiAXaSxeUdKWtQFJW5AuBQV4o/tFzMysXRR9eP2vJC1NqiJ6HThE0jDge2UGZ2Zm5Sty1dBtwMaSVgMUEc/UTb6gtMjMzKwhipQIkPRh4H2km8oAiIj/KTEuMzNrkCKXj/4UGA58EPgl8EngtpLjMjOzBilSItghIjaVdG9EHCvpe8BlZQdWtla+3dvMFmmlq2sGqyJXDb2SX1+V9A/Aq8Do0iIyM7OGKlIi+J2k4cDJwD2km8nOKTUqMzNrmCJXDR2X316Snzu8nLuVMDMbPIo8oewTklbMg18CTpe0WblhmZlZoxRpIzguIuZL+gDwUeAi4IxywzIzs0YpkggW5te9gJ9GxGXAu8oLyczMGqlIY3GnpNOA3YFxkpahWAIxM7M2UOSEvj/wB+AjEfE8MAI4utSozMysYYqUCEYAV0XEa5J2ADYFzi03LGtXvlHPrP0UKRFcCbwpaX3S/QMbseghNWZm1uaKJII3I2IB8AngRxHxRWCtvhaStKyk2yXdK2mapOPz+PUkTZH0sKSLcpuDmZk1SZGqoTck7QccDOyTxy1dYLnXgJ0j4sX8LINbJP0O+DLww4i4UNLpwCHAzxYjdjNrksFaBVjVfo2KlAg+B+wEfC8iHpO0HgWeQxDJi3lw6fwXwM7ApXn8JBYlFzMza4I+E0FE3A8cRe56OiIej4gTiqxc0lBJ9wBzgOuBR4F5EVF7xOVMClQzmZlZeYo8j+AjwCnAMsB6ksYC/xURH+9r2YhYCIzNndZdQWpofsdsPWx3AjABYNSoUX1tysysFIO1GqxekaqhbwPbAPMAIuIeYIP+bCQi5gE3AdsCwyXVEtDawKwelpkYEeMiYlxHR0d/NmdmZv1QpLF4QUTMqz2iMuv2V3w9SR11yy4H7AqcBNwI7AtcCIwHrup31AZUt2HLzAZWkUQwXdL+wJDcUPwlij2qciQwSdJQUsnj4oi4RtJfgQslfQe4GzhrMWM3M7MBUCQRfAH4T+BNUj3/tcA3+1ooIv4CbN7N+MeArfsXppmZlaXIg2leAr6R/8zMbJApctXQFqRO5kbXzx8RW5QXlpmZNUqRqqHzSVVB95Gqh8zMbBApkgiejYjLS4/ErJ+qcH23tYd2PxaLJILjJZ0B/B+p/yAAIuLq0qIyM7OGKZIIDiI9g2AFFlUNBeBEYGY2CBRJBFtGxCalR2JmZk1RpIuJKZI2LD0SMzNriiIlgq2Bv0h6hNRGIFIv07581MxsECiSCPy8ADOzQazIncWPNiIQMzNrjiJtBGZmNogVqRpqKUVu3Gj3mzvMzBqpzxKBpP8pMs7MzNpTkaqh3bsZ95GBDsTMzJqjx6ohSYcDRwAbSrqrbtKKwJ1lB2ZWz09js974+FgyvbURXAxMBr5L6oa6Zn5EzCk1KjMza5geE0FEPA88L+m0rpeQSjooIs4rPTozMytdkauGTshVQ18ndTw3kXR3sRNBN7ororbDVUztEKM1ho+F9jIQ1WJFGov/EXiK9KD5W4HLI8J3G5uZDRJFSgQrAZsBM4GRwBqSFBHR20KS1gHOAf6B1H31xIg4VdKqwEWkR1/OAPbP1VBm1oLcEDv4FSkR3A7cGBG7kjqgezdwc4Hl3gC+EhEbAdsCn5e0ManheXJEjCE1Rh/dyzrMzKxkRUoEH46IxwEi4iXg3yT1+ZMgIjqBzvx+vqTpwFrA3sCOebZJwE3AN/oduZmZDYginc49LukAYP2IOCFX+fy9PxuRNBrYHJgCrJGTBBHRKWn1HpaZAEwAGDVqVH82ZxXlKoz244bp1lCki4mfADsBn86jXgJOL7oBSSsAlwFHRcQLRZeLiIkRMS4ixnV0dBRdzMzM+qlIG8EHIuJw4FWAiHgOWKbIyiUtTUoC50XE5Xn0bEkj8/SRgG9OMzNroiJtBAskDSE9sB5Jq7HoIfY9kiTgLGB6RJxSN+lqYDxwYn69qr9BV9XiFKPbtbrEVQa9a7X/a6vFY/3TY4lAUi1JnEb6Vd8h6XjgFuCkAuveHjgY2FnSPflvT1IC2E3Sw8BuedjMzJqktxLB7cAWEXGOpDuBXUl3FO8XEff3teKIuCXP351d+h2pmZmVordE8NZJPCKmAdPKD8damatrbEm4+qh19ZYIOiR9uaeJXer9zcysTfWWCIaSOpnrqXrHzMwGgd4SQWdEfLthkZg1gKsn2l+R/6H/z/3T230ELgmYmVVAb4nAV/aYmVVAb08oe66RgZi1O1dHWLsq0sWEmZkNYk4EZmYVV6SvIctc9Lf+WNzjpV2fe23tyyUCM7OKcyIwM6s4Vw1lrvax3vj4sGZo1HHnEoGZWcW1RYlg+tPz2er7N/hX2CDV7F/bboi1/hiMx6tLBGZmFedEYGZWcW1RNWQDy703WqvycdccLhGYmVWcE4GZWcWVVjUk6RfAXsCciNgkj1sVuAgYDcwA9o+I58uKYUkVaZ1fnBZ8V81YFfhqrPZRZongbGD3LuOOBiZHxBhgch42M7MmKi0RRMQfga7PNNgbmJTfTwL2KWv7ZmZWTKOvGlojIjoBIqJT0uo9zShpAjABYJnhPc5mNuBcLWdV07KNxRExMSLGRcS4pYYNb3Y4ZmaDVqMTwWxJIwHy65wGb9/MzLpodNXQ1cB44MT8elWDt28N5moWs9ZXWolA0gXAn4ENJc2UdAgpAewm6WFgtzxsZmZNVFqJICIO7GHSLmVt08zM+q/t+xpy1YOZ2ZJp2auGzMysMZwIzMwqru2rhopw9VFrqXIfNM0+Fqu8761nLhGYmVWcE4GZWcVVomqoSrpWPVRdWVUhza7iMRtILhGYmVWcSwRmZi2iWY35LhGYmVWcE4GZWcW5asjMKsEN/D1zicDMrOKcCMzMKs6JwMys4pwIzMwqzonAzKzifNWQWZO5R9Dm8b5PXCIwM6s4JwIzs4prSiKQtLukByU9IunoZsRgZmZJwxOBpKHAacAewMbAgZI2bnQcZmaWNKNEsDXwSEQ8FhGvAxcCezchDjMzAxQRjd2gtC+we0QcmocPBraJiC90mW8CMCEPbgLc39BAl9wI4JlmB9FP7RgztGfcjrkx2jFmGLi4142Ijr5masblo+pm3DuyUURMBCYCSJoaEePKDmwgOebGace4HXNjtGPM0Pi4m1E1NBNYp254bWBWE+IwMzOakwjuAMZIWk/SMsABwNVNiMPMzGhC1VBEvCHpC8C1wFDgFxExrY/FJpYf2YBzzI3TjnE75sZox5ihwXE3vLHYzMxai+8sNjOrOCcCM7OKa+lE0C5dUUj6haQ5ku6vG7eqpOslPZxfV2lmjF1JWkfSjZKmS5om6Ut5fMvGLWlZSbdLujfHfHwev56kKTnmi/JFCC1F0lBJd0u6Jg+3Q8wzJN0n6R5JU/O4lj0+ACQNl3SppAfysb1dK8csacO8f2t/L0g6qtExt2wiaLOuKM4Gdu8y7mhgckSMASbn4VbyBvCViNgI2Bb4fN6/rRz3a8DOEbEZMBbYXdK2wEnAD3PMzwOHNDHGnnwJmF433A4xA+wUEWPrrmlv5eMD4FTg9xHxXmAz0j5v2Zgj4sG8f8cCWwIvA1fQ6JgjoiX/gO2Aa+uGjwGOaXZcvcQ7Gri/bvhBYGR+PxJ4sNkx9hH/VcBu7RI3sDxwF7AN6Q7Mpbo7blrhj3SvzGRgZ+Aa0k2VLR1zjmsGMKLLuJY9PoCVgMfJF8G0Q8xd4vwQ8KdmxNyyJQJgLeDJuuGZeVy7WCMiOgHy6+pNjqdHkkYDmwNTaPG4cxXLPcAc4HrgUWBeRLyRZ2nF4+RHwNeBN/PwarR+zJDu+L9O0p25yxdo7ePj3cBc4Je5Gu5MScNo7ZjrHQBckN83NOZWTgSFuqKwJSNpBeAy4KiIeKHZ8fQlIhZGKkavTerAcKPuZmtsVD2TtBdwELXHAAAEbklEQVQwJyLurB/dzawtE3Od7SNiC1L17OclfbDZAfVhKWAL4GcRsTnwEi1UDdSb3Eb0MeCSZmy/lRNBu3dFMVvSSID8OqfJ8byDpKVJSeC8iLg8j275uAEiYh5wE6l9Y7ik2s2RrXacbA98TNIMUk+7O5NKCK0cMwARMSu/ziHVW29Nax8fM4GZETElD19KSgytHHPNHsBdETE7Dzc05lZOBO3eFcXVwPj8fjypDr5lSBJwFjA9Ik6pm9SycUvqkDQ8v18O2JXUGHgjsG+eraVijohjImLtiBhNOoZviIiDaOGYASQNk7Ri7T2p/vp+Wvj4iIingSclbZhH7QL8lRaOuc6BLKoWgkbH3OwGkj4aT/YEHiLVA3+r2fH0EucFQCewgPSr5BBSPfBk4OH8umqz4+wS8w6k6oi/APfkvz1bOW5gU+DuHPP9wH/m8e8GbgceIRWt39XsWHuIf0fgmnaIOcd3b/6bVvv+tfLxkeMbC0zNx8iVwCptEPPywLPAynXjGhqzu5gwM6u4Vq4aMjOzBnAiMDOrOCcCM7OKcyIwM6s4JwIzs4pzIrBBQ9Jqdb04Pi3pqbrhW0va5uaSzuxh2gxJIwZwWxdKGjNQ6zOr8eWjNihJOg54MSJOLnk7lwDfiYh7u5k2AxgXEc8M0Lb+Cfh0RBw2EOszq3GJwCpB0ov5dUdJf5B0saSHJJ0o6aD8nIP7JK2f5+uQdJmkO/Lf9t2sc0Vg01oSyCWS63KHZ2dQ16eQpCtz523Tah24STpE0g/r5jlM0in5rt7fKD134X5Jn8qz3AzsWtc1hdmAcCKwKtqM9HyA9wMHA++JiK2BM4Ev5nlOJT0vYCvgk3laV+NIdzjX/BdwS6QOz64GRtVN+1xEbJmXOVLSaqS+hz6W+3wC+CzwS9KzLWZFxGYRsQnwe4CIeJN0J/JmS/LhzbryLwurojsid/Er6VHgujz+PmCn/H5XYOPUJRMAK0laMSLm161nJKnb45oPAp8AiIjfSHq+btqRkj6e368DjImI2yTdAOwlaTqwdETcJ+k14GRJJ5G6pLi5bj1zgDWB+t5MzZaIE4FV0Wt179+sG36TRd+JIcB2EfFKL+t5BVi2y7h3NLpJ2pGUWLaLiJcl3VS33JnAN4EHSKUBIuIhSVuS+n76rqTrIuLbef5l83bNBoyrhsy6dx3whdqApLHdzDMd2KBu+I/AQXn+PUgdngGsDDyfk8B7SV1nAxCpy+R1gH8m9z4paU3g5Yg4FziZ1JVyzXtIncCZDRiXCMy6dyRwmqS/kL4nfwSOqJ8hIh6QtHJdldHxwAWS7gL+APwtz/p74Ii8rgeB27ps62JgbETUqpLeD3xf0pukHm3/FUDSGsArtWots4Hiy0fNloCkfwfmR0S39xIUXMc1pIbpyQW29UJEnLW42zLrjquGzJbMz3h7m0NhkoZLeoj0K7/XJJDNAyYtzrbMeuMSgZlZxblEYGZWcU4EZmYV50RgZlZxTgRmZhXnRGBmVnH/D7GLAMGrTFY6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_count_data=len(count_data)\n",
    "plt.bar(np.arange(n_count_data), count_data,color=\"#348ABD\")\n",
    "plt.xlabel('Time (days)')\n",
    "plt.ylabel('Text messages received')\n",
    "plt.title(\"Did the user''s texting habits change over time?\")\n",
    "plt.xlim(0,n_count_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymc as pm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 1.0 / count_data.mean()\n",
    "lambda_1=pm.Exponential(\"lambda_1\",alpha)\n",
    "lambda_2=pm.Exponential(\"lambda_2\",alpha)\n",
    "tau=pm.DiscreteUniform(\"tau\",lower=0,upper=n_count_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Random output: 16 10 70\n"
     ]
    }
   ],
   "source": [
    "print(\"Random output:\",tau.random(),tau.random(),tau.random())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "@pm.deterministic\n",
    "def lambda_(tau=tau,lambda_1=lambda_1, lambda_2=lambda_2):\n",
    "    out=np.zeros(n_count_data)\n",
    "    out[:tau] = lambda_1\n",
    "    out[tau:] = lambda_2\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "observation = pm.Poisson(\"obs\", lambda_, value=count_data,\n",
    "                        observed=True)\n",
    "model = pm.Model([observation, lambda_1, lambda_2, tau])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " [-----------------100%-----------------] 40000 of 40000 complete in 10.0 sec"
     ]
    }
   ],
   "source": [
    "mcmc = pm.MCMC(model)\n",
    "mcmc.sample(40000, 10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "lambda_1_samples = mcmc.trace(\"lambda_1\")[:]\n",
    "lambda_2_samples = mcmc.trace(\"lambda_2\")[:]\n",
    "tau_samples = mcmc.trace('tau')[:]\n",
    "\n",
    "figsize(14.5, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\users\\shojiro_shibayama\\miniconda3\\envs\\bayesian\\lib\\site-packages\\matplotlib\\axes\\_axes.py:6521: MatplotlibDeprecationWarning: \n",
      "The 'normed' kwarg was deprecated in Matplotlib 2.1 and will be removed in 3.1. Use 'density' instead.\n",
      "  alternative=\"'density'\", removal=\"3.1\")\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Probability')"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2IAAAJgCAYAAADoENAtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XuYHVWZsP37ITRGDiKGIEJCEh0YwiQaYgiMUTwzIJmgDjPCmBdQGQ6BF8XhEMQLkeH1U5kXlG8IyMhRgURRmYgBREQHEJAEoyTEfGIMpI0jMZzBcHy+P3YFN53d3buTrtq7u+/fdfWVXVVrVT21V+9OPXutWhWZiSRJkiSpOpu1OgBJkiRJGmpMxCRJkiSpYiZikiRJklQxEzFJkiRJqpiJmCRJkiRVzERMkiRJkipmIiZJkiRJFTMRkyRJkqSKmYhJUj+JiKUR8a4S9395RJxdxrHq9xcRKyPifWXsux1ExF9HxC8i4smIOKGJ8v36fkiSBCZikgao4uL4zxHxVET8MSIui4itN3F/m3SxnZl/k5k/2ZR99Pexmj2v/oq90fGqfF+adArwk8zcJjPP77rRxGtDrXhPImJZRHRGxN+0w34kqb+ZiEkayP4+M7cGJgN7AZ9tRRARsXkr6w/UY7fQGGBpq4PYFAOt3TYy3gnA/wf8wyYevr/2I0n9ykRM0oCXmb8HbqB2wUVEjI+In0TEY8WwuBnry0bEqRHx+2JY2vKIeG9EfAPYBfh+0cN2SlF2p4j4TkSsiYjf1Q9jK3oITo2IXwFPR8TmXXsNeoljg/pdzysi9oyIe4tY5wHDu9SvP1ZfzqvX2IG9IuL+iHi06G2sP3ZGxF/VLV8eEWf3crz39fae1JU9KSJ+FRGPR8S89cdudI6Nfh96ed9/DLwb+I8ixt261G14DsCkRjEVdbr9PWkQ28qIOK2H93Z2RPy2OMf7I+JDXep2bbfeyp9cxP10RFwSEa+PiBuK8j+KiO16O4+N+Xx0E29T7bdeZr4I3A68padyvemv/UStZ+2piHiu+Hmq+Bm/KfuVNIRlpj/++OPPgPsBVgLvK16PptbD8W9AB/AA8BlgC+A9wJPAXxc/q4CdinpjgTd13V+xvBmwCDij2M8bgRXA39WVX1wc+9UNYuo2ju7qdzm/LYAHgROLfR0MPA+c3eBYTZ9Xk7GvBJYU218H3LH+uMX2BP6qbvnyRnF1bave3pO6sj8HdiqOvQw4pqdz7HKsZo7xE+DIZn63eoqpmd+Tbvbd03v7j8VxNgM+AjwNvKGHduut/F3A64GdgYeBe4E9gVcBPwY+14ff96Y/Hw3ifUsz7dflvXo1tZ6s32zi34p+2U/d/i4BTi/z75s//vgzNH7sEZM0kF0XEY9R+7b7p8AXgH2ArYEvZuZzmflj4HrgUOBFahege0RER2auzMzfdrPvvYCRmXlWsZ8VwH8Ch9SVOT8zV2XmnxvU7ymOZut3AF/JzOcz81rgnm5i7ct5NXNsgP8otj8C/J8ucW+sZt6T9bGtLo79fWASzZ9js8foq0YxQXO/J111+95m5reL47yUmfOA3wBTu8Txcrs1Uf7/zcw/Zq3X+Dbg7sz8RWY+C3yPWlK2MefRbPnzM3MVtQSxr7+j/wf4PfCmKO7/jIhtI+LnRU/UhF7qd7ufYl9/GxF3RsRPI+KaiOhocn9vppZMS9ImMRGTNJB9MDNfm5ljMnNWcXG6E7AqM1+qK/cgsHNmPgB8CjgTeDgi5kbETt3sewywUzG87bEi4fsMtd6F9Vb1EFu3cfSh/u8zM7vU30Afz6uZY3fd/mARz6Zq5j0B+J+6188AW/fhHJs9Rl9tEFPxupnfk666fW8j4rCIWFy3rwnA9t3Ubab8H+te/7nB8saeR7PlV0Hff0cj4m+Bf6J2X9fjxXlB7b0/ELi2u7pN7gdq7/17MvOd1HrzDmpif5sBe2AiJqkfmIhJGmxWA6OLC6b1dqH2jTiZeXVmvp3ahWQCXyrKJK+0Cvhdkeit/9kmMz9QV6ZrnabjaKL+H4CdIyK61G+oD+fVzLGhNpys/rir65afAbasW96xyf028550q4dz7LdjrD9UH8o283vSVcP3NiLGUOtVOh4YkZmvpXbBX/878HJsTZbvr/PYmM/HK+o12X4U98xdSm345yPALynu7yp6h9c0c0I97afY1+q6HuEXgJc23MsGdqF27bSimRgkqScmYpIGm7upDYM6JSI6ovb8qr8H5kbt+VHviYhXAeuo9Qi8WNT7I7X7XNb7OfBEMcHAqyNiWERMiIi9NjWOJuvfSe3i8IRiUoYP88ohZy/r43k167iIGBURr6PW0zGvbtti4J+L92R/4J1123o63ka/J72cY78co8lz6Gpjfk+6e2+3opagrAGIiI/xyh6crvpaflPOY5M+H31oP4CzgDsz8/pieTG14YDditqEMZdvzH4iYhxwALUhrN3ta73XUPv92qKneCSpGSZikgaVzHwOmEHtwupPwBzgsMz8NbV7VL5YrP8fYAdqF8IA/w/w2WKY1UlZm2nt76ndC/S7os7XgW37IY5m638YOAJ4lNpEDN/tpnjT59XMsQtXAz+k9s3/CuDsum2fpPbePAZ8FLiublu3x9vE96Snc+yvY/R6Dg2OtzG/Jw3f28y8H/i/1JLwPwITqU3m0d2x+1R+E89jUz8fTbVfREylNgHJiXWrF9P7jIejqTv3ZvcTEa8BrgD+V/G7s8G+ulhGrWft0YjYvZeYJKlH8crbDyRJUlkiYiW1GRt/1OpYBrKix+rfM3NJRGxBLTl6c2Y+34d9bA78F/B/i0ld2Nh9SdLGsEdMkiQNGBGxANgP+M+IOKKYtXH8RiROhwJ7A2dE7blzH9mEfUlSn5XWIxYRlwLTgYczc4Mx68UN6F8FPkDtxu8jMvPeUoKRJKkN2CMmSVqvzB6xy4H9e9h+ALBr8XMUcGGJsUiS1HKZOdYkTJIEJSZimfnfwCM9FDkIuDJr7gJeGxFvKCseSZIkSWoXm7fw2DvzygdTdhbr/tC1YEQcRa3XjK222uqtu+/uREWSJEmS2sOiRYv+lJkj+1KnlYlYowdONrxhLTMvBi4GmDJlSi5cuLDMuCRJkiSpaRHxYF/rtHLWxE5qz+pYbxSwukWxSJIkSVJlWpmIzQcOi5p9gMczc4NhiZIkSZI02JQ2NDEirgHeBWwfEZ3A54AOgMy8CFhAber6B6hNX/+xsmKRJEmSpHZSWiKWmYf2sj2B4/rjWM8//zydnZ2sW7euP3ankg0fPpxRo0bR0dHR6lAkSZKklmjlZB39prOzk2222YaxY8dSe0602lVmsnbtWjo7Oxk3blyrw5EkSZJaopX3iPWbdevWMWLECJOwASAiGDFihL2XkiRJGtIGRSIGmIQNILaVJEmShrpBk4hJkiRJ0kBhIiZJkiRJFTMRkyRJkqSKDYpZE7u6fsIB/bq/6Utu6Nf9deexxx7j6quvZtasWX2u+7a3vY2f/exnJUQF559/PhdeeCGTJ0/mqquuenn9fffdx/Tp05k9ezbHHntsKceWJEmSBiN7xNrIY489xpw5c/pUJzN56aWX+pSEra/TrDlz5rBgwYJXJGEAEydOZO7cuVx55ZVN70uSJEmSiVi/WblyJbvvvjuHH344b37zmzn44IN55plnADj33HOZMGECEyZM4Ctf+QoATz/9NAceeCBvectbmDBhAvPmzWP27Nn89re/ZdKkSZx88skAfPOb32Tq1KlMmjSJo48+mhdffJGVK1cyfvx4Zs2axeTJk1m1ahVbb731y7E0Ol6jOl01qnfMMcewYsUKZsyYwXnnnbdBnR122IGlS5f275spSZIkDXKDcmhiqyxfvpxLLrmEadOm8fGPf5w5c+bw7ne/m8suu4y7776bzGTvvffmne98JytWrGCnnXbiBz/4AQCPP/44e++9N0uWLGHx4sUALFu2jHnz5nHHHXfQ0dHBrFmzuOqqq9h3331Zvnw5l1122QY9aIsWLWp4vO22267bOj3Vu+iii7jxxhu59dZb2X777TeoN3v2bJ599lkefPBBxowZU8K7KkmSJA0+9oj1o9GjRzNt2jQAZs6cye23387tt9/Ohz70Ibbaaiu23nprPvzhD3PbbbcxceJEfvSjH3Hqqady2223se22226wv1tuuYVFixax1157MWnSJG655RZWrFgBwJgxY9hnn302qNPd8Xqq01u97tx4440v9+yt7xVbsWIFn/jEJzj44IObf+MkSZKkIcZErB91fVBxRJCZDcvutttuLFq0iIkTJ3Laaadx1llnbVAmMzn88MNZvHgxixcvZvny5Zx55pkAbLXVVg33293xeqrTW71G1q1bxymnnMKcOXOYOHEiS5YsAeCNb3wjl1xySZ/2JUmSJA01JmL96KGHHuLOO+8E4JprruHtb387++67L9dddx3PPPMMTz/9NN/73vd4xzvewerVq9lyyy2ZOXMmJ510Evfeey/bbLMNTz755Mv7e+9738u1117Lww8/DMAjjzzCgw8+2GMM3R2vN32td/bZZ3PYYYcxduzYVyRikiRJkno3KO8Rq2q6+a7Gjx/PFVdcwdFHH82uu+7Ksccey5ZbbskRRxzB1KlTATjyyCPZc889uemmmzj55JPZbLPN6Ojo4MILL2TEiBFMmzaNCRMmcMABB3DOOedw9tlns99++/HSSy/R0dHBBRdcwI477thtDJMnT254vJUrV/YYe3f1Glm+fDk333wzd9xxB1CbPfELX/hCX98uSZIkaciKvg5Ja7UpU6bkwoULX7Fu2bJljB8/vkUR1axcuZLp06cP+Z6htWvXcvrpp3PzzTdz5JFHctpppzUs1w5tJkmSJPWHiFiUmVP6UmdQ9oipdUaMGMFFF13U6jAkSZKktuY9Yv1k7NixQ743TJIkSVJzTMQkSZIkqWImYpIkSZJUsUGTiA20SUeGMttKkiRJQ92gSMSGDx/O2rVrvcAfADKTtWvXMnz48FaHIkmSJLXMoJg1cdSoUXR2drJmzZpWh6ImDB8+nFGjRrU6DEmSJKllBkUi1tHRwbhx41odhiRJkiQ1pdShiRGxf0Qsj4gHImJ2g+27RMStEfGLiPhVRHygzHgkSZIkqR2UlohFxDDgAuAAYA/g0IjYo0uxzwLfysw9gUOAOWXFI0mSJEntoswesanAA5m5IjOfA+YCB3Upk8BritfbAqtLjEeSJEmS2kKZidjOwKq65c5iXb0zgZkR0QksAP53ox1FxFERsTAiFjohhyRJkqSBrsxELBqs6zq//KHA5Zk5CvgA8I2I2CCmzLw4M6dk5pSRI0eWEKokSZIkVafMRKwTGF23PIoNhx5+AvgWQGbeCQwHti8xJkmSJElquTITsXuAXSNiXERsQW0yjvldyjwEvBcgIsZTS8QceyhJkiRpUCstEcvMF4DjgZuAZdRmR1waEWdFxIyi2L8C/xIRvwSuAY7IzK7DFyVJkiRpUCn1gc6ZuYDaJBz1686oe30/MK3MGCRJkiSp3ZT6QGdJkiRJ0oZMxCRJkiSpYiZikiRJklQxEzFJkiRJqpiJmCRJkiRVzERMkiRJkipmIiZJkiRJFTMRkyRJkqSKmYhJkiRJUsVMxCRJkiSpYiZikiRJklQxEzFJkiRJqpiJmCRJkiRVzERMkiRJkipmIiZJkiRJFTMRkyRJkqSKmYhJkiRJUsU2b3UAUpWun3BAj9unL7mhokgkSZI0lNkjJkmSJEkVMxGTJEmSpIqZiEmSJElSxUzEJEmSJKliTtahQaO3iTgkSZKkdlFqj1hE7B8RyyPigYiY3U2Zf4qI+yNiaURcXWY8kiRJktQOSusRi4hhwAXA+4FO4J6ImJ+Z99eV2RU4DZiWmY9GxA5lxSNJkiRJ7aLMHrGpwAOZuSIznwPmAgd1KfMvwAWZ+ShAZj5cYjySJEmS1BaaSsQi4jsRcWBE9CVx2xlYVbfcWayrtxuwW0TcERF3RcT+3Rz/qIhYGBEL16xZ04cQJEmSJKn9NJtYXQj8M/CbiPhiROzeRJ1osC67LG8O7Aq8CzgU+HpEvHaDSpkXZ+aUzJwycuTIJkOWJEmSpPbUVCKWmT/KzI8Ck4GVwM0R8bOI+FhEdHRTrRMYXbc8CljdoMx/Zebzmfk7YDm1xEySJEmSBq2mJ+uIiBHATOB/Ab8ArgLeDhxOrUerq3uAXSNiHPB74BBqvWr1rqPWE3Z5RGxPbajiir6dgtR/mpkCf/qSGyqIRJIkSYNZU4lYRHwX2B34BvD3mfmHYtO8iFjYqE5mvhARxwM3AcOASzNzaUScBSzMzPnFtv0i4n7gReDkzFy7aackSZIkSe2t2R6xr2fmgvoVEfGqzHw2M6d0V6mos6DLujPqXifw6eJHkiRJkoaEZifrOLvBujv7MxBJkiRJGip67BGLiB2pTTn/6ojYk7/MhPgaYMuSY5MkSZKkQam3oYl/BxxBbcbDc+vWPwl8pqSYpLbW24QeTuYhSZKk3vSYiGXmFcAVEfEPmfmdimKSJEmSpEGtt6GJMzPzm8DYiNhgQo3MPLdBNUmSJElSD3obmrhV8e/WZQci9aaZZ3xJkiRJA0FvQxO/Vvz7+WrCkSRJkqTBr6np6yPiyxHxmojoiIhbIuJPETGz7OAkSZIkaTBq9jli+2XmE8B0oBPYDTi5tKgkSZIkaRBrNhHrKP79AHBNZj5SUjySJEmSNOj1NlnHet+PiF8DfwZmRcRIYF15YUmSJEnS4NVUj1hmzgb+FpiSmc8DTwMHlRmYJEmSJA1WzfaIAYyn9jyx+jpX9nM8kiRJkjToNZWIRcQ3gDcBi4EXi9WJiZj6yWB6Rlgz5zJ9yQ0VRCJJkqR21WyP2BRgj8zMMoORJEmSpKGg2VkTlwA7lhmIJEmSJA0VzfaIbQ/cHxE/B55dvzIzZ5QSlSRJkiQNYs0mYmeWGYQkSZIkDSVNJWKZ+dOIGAPsmpk/iogtgWHlhiZJkiRJg1NT94hFxL8A1wJfK1btDFxXVlCSJEmSNJg1O1nHccA04AmAzPwNsENZQUmSJEnSYNZsIvZsZj63fqF4qLNT2UuSJEnSRmg2EftpRHwGeHVEvB/4NvD98sKSJEmSpMGr2URsNrAGuA84GlgAfLa3ShGxf0Qsj4gHImJ2D+UOjoiMiClNxiNJkiRJA1azsya+FBHXAddl5ppm6kTEMOAC4P1AJ3BPRMzPzPu7lNsGOAG4u0+RS5IkSdIA1WOPWNScGRF/An4NLI+INRFxRhP7ngo8kJkrivvL5gIHNSj3b8CXgXV9jF2SJEmSBqTeesQ+RW22xL0y83cAEfFG4MKIODEzz+uh7s7AqrrlTmDv+gIRsScwOjOvj4iTuttRRBwFHAWwyy679BKy1P6un3BAj9unL7mhokgkSZLUCr3dI3YYcOj6JAwgM1cAM4ttPYkG616eaTEiNgPOA/61tyAz8+LMnJKZU0aOHNlbcUmSJElqa70lYh2Z+aeuK4v7xDp6qdsJjK5bHgWsrlveBpgA/CQiVgL7APOdsEOSJEnSYNdbIvbcRm4DuAfYNSLGRcQWwCHA/PUbM/PxzNw+M8dm5ljgLmBGZi5sIm5JkiRJGrB6u0fsLRHxRIP1AQzvqWJmvhARxwM3AcOASzNzaUScBSzMzPk91ZckSZKkwarHRCwzh23KzjNzAbVnjtWvazjjYma+a1OOJUmSJEkDRbMPdJYkSZIk9RMTMUmSJEmqWG/3iEn9orfnZkmSJElDiT1ikiRJklQxEzFJkiRJqpiJmCRJkiRVzERMkiRJkipmIiZJkiRJFTMRkyRJkqSKmYhJkiRJUsVMxCRJkiSpYj7QWZvMhzX3v2be0+lLbqggEkmSJJXBHjFJkiRJqpiJmCRJkiRVzERMkiRJkipmIiZJkiRJFTMRkyRJkqSKmYhJkiRJUsVMxCRJkiSpYiZikiRJklQxEzFJkiRJqtjmrQ5A0sa5fsIBPW6fvuSGiiKRJElSX9kjJkmSJEkVKzURi4j9I2J5RDwQEbMbbP90RNwfEb+KiFsiYkyZ8UiSJElSOyhtaGJEDAMuAN4PdAL3RMT8zLy/rtgvgCmZ+UxEHAt8GfhIWTFp4/Q2BE6SJElS35TZIzYVeCAzV2Tmc8Bc4KD6Apl5a2Y+UyzeBYwqMR5JkiRJagtlJmI7A6vqljuLdd35BNBwdoGIOCoiFkbEwjVr1vRjiJIkSZJUvTITsWiwLhsWjJgJTAHOabQ9My/OzCmZOWXkyJH9GKIkSZIkVa/M6es7gdF1y6OA1V0LRcT7gNOBd2bmsyXGI0mSJEltocwesXuAXSNiXERsARwCzK8vEBF7Al8DZmTmwyXGIkmSJElto7RELDNfAI4HbgKWAd/KzKURcVZEzCiKnQNsDXw7IhZHxPxudidJkiRJg0aZQxPJzAXAgi7rzqh7/b4yjy9JkiRJ7ajUREwDg88JkyRJkqpV5j1ikiRJkqQGTMQkSZIkqWIOTZQGqd6GnE5f0vD56ZIkSaqAPWKSJEmSVDETMUmSJEmqmImYJEmSJFXMREySJEmSKmYiJkmSJEkVMxGTJEmSpIo5fb00RPU2vT04xb0kSVJZ7BGTJEmSpIqZiEmSJElSxUzEJEmSJKliJmKSJEmSVDEn6xjkmpmQQZIkSVK17BGTJEmSpIrZIyapW731qDq9vSRJ0saxR0ySJEmSKmaP2ADnPWCSJEnSwGOPmCRJkiRVzB6xNmZvlyRJkjQ4mYhJ2mjNfFnghB6SJEkbKjURi4j9ga8Cw4CvZ+YXu2x/FXAl8FZgLfCRzFxZZkySquXMi5IkSRsq7R6xiBgGXAAcAOwBHBoRe3Qp9gng0cz8K+A84EtlxSNJkiRJ7aLMHrGpwAOZuQIgIuYCBwH315U5CDizeH0t8B8REZmZJcbVNrwHTOqfz4G9apIkaaApMxHbGVhVt9wJ7N1dmcx8ISIeB0YAf6ovFBFHAUcVi09FxPJSIm6t7ely3mprtlc7ieithO01sNheA4vtNXDYVgOL7TWw/HVfK5SZiDW6Mura09VMGTLzYuDi/giqXUXEwsyc0uo41Bzba2CxvQYW22tgsb0GDttqYLG9BpaIWNjXOmU+R6wTGF23PApY3V2ZiNgc2BZ4pMSYJEmSJKnlykzE7gF2jYhxEbEFcAgwv0uZ+cDhxeuDgR8PlfvDJEmSJA1dpQ1NLO75Oh64idr09Zdm5tKIOAtYmJnzgUuAb0TEA9R6wg4pK54BYFAPvRyEbK+BxfYaWGyvgcX2Gjhsq4HF9hpY+txeYQeUJEmSJFWrzKGJkiRJkqQGTMQkSZIkqWImYi0QEZdGxMMRsaTL+v8dEcsjYmlEfLlV8emVGrVXREyKiLsiYnFELIyIqa2MUTURMToibo2IZcXn6JPF+tdFxM0R8Zvi3+1aHat6bK9zIuLXEfGriPheRLy21bGq+/aq235SRGREbN+qGPUXPbWX1xvtp4e/h15vtKGIGB4RP4+IXxbt9fli/biIuLu43phXTFjY/X68R6x6EbEv8BRwZWZOKNa9GzgdODAzn42IHTLz4VbGqZpu2uuHwHmZeUNEfAA4JTPf1cIwBUTEG4A3ZOa9EbENsAj4IHAE8EhmfjEiZgPbZeapLQxV9Nheo6jNovtCRHwJwPZqve7aKzPvj4jRwNeB3YG3ZqYPoW2xHj5fr8frjbbTQ3t9Ba832k5EBLBVZj4VER3A7cAngU8D383MuRFxEfDLzLywu/3YI9YCmfnfbPi8tGOBL2bms0UZ/yi2iW7aK4HXFK+3ZcNn5KkFMvMPmXlv8fpJYBmwM3AQcEVR7Apq/7mpxbprr8z8YWa+UBS7i1piphbr4fMFcB5wCrW/jWoDPbSX1xttqIf28nqjDWXNU8ViR/GTwHuAa4v1vV5vmIi1j92AdxTdmT+NiL1aHZB69CngnIhYBfw7cFqL41EXETEW2BO4G3h9Zv4Bav/ZATu0LjI10qW96n0cuKHqeNSz+vaKiBnA7zPzly0NSt3q8vnyeqPNdWkvrzfaVEQMi4jFwMPAzcBvgcfqvkjs5C9fVjVkItY+Nge2A/YBTga+VXR7qj0dC5yYmaOBE6k9E09tIiK2Br4DfCozn2h1POpZd+0VEacDLwBXtSo2bai+vai1z+nAGS0NSt1q8PnyeqONNWgvrzfaVGa+mJmTqI3amAqMb1Ssp32YiLWPTmpjSjMzfw68BHjDc/s6HPhu8frb1D6AagPFWO3vAFdl5vo2+mMx/n79OHyH4rSJbtqLiDgcmA58NL2ZuW00aK83AeOAX0bESmoXJPdGxI6ti1LrdfP58nqjTXXTXl5vtLnMfAz4CbUvN14bEZsXm0bRy1BSE7H2cR21caVExG7AFoA3O7ev1cA7i9fvAX7TwlhUKL7VvQRYlpnn1m2aT+0/M4p//6vq2LSh7torIvYHTgVmZOYzrYpPr9SovTLzvszcITPHZuZYahf5kzPzf1oYqujx76HXG22oh/byeqMNRcTI9TP6RsSrgfdRu6/vVuDgoliv1xvOmtgCEXEN8C5q30D9Efgc8A3gUmAS8BxwUmb+uFUx6i+6aa/lwFepDfFYB8zKzEWtilE1EfF24DbgPmrf8gJ8hto4+28BuwAPAf+YmV0nYFHFemiv84FXAWuLdXdl5jHVR6h63bVXZi6oK7MSmOKsia3Xw+frR3i90XZ6aK8n8Hqj7UTEm6lNxjGMWsfWtzLzrIh4IzAXeB3wC2Dm+olxGu7HREySJEmSquXQREmSJEmqmImYJEmSJFXMREySJEmSKmYiJkmSJEkVMxGTJEmSpIqZiEmSJElSxUzEJEmSJKliJmKSpAEjIiZGxIMRcWzJx3mqzP1LkmQiJkkaMDLzPuAQ4LBWxyJJ0qYwEZMkDTQPA3/TbOGI+FJEzKpbPjMi/rV4fV1ELIqIpRFxVIO6YyNiSd3ySRFxZvF6ZkT8PCIWR8TXImLYppyUJGloMRGTJA00XwReFRFjmiw/F/hI3fI/Ad8uXn88M98KTAFOiIgRzewwIsYX+5yWmZOAF4GPNhmPJEls3uoAJElqVkTsD2wF/IBar9iDEfFG4HRQwp3lAAAgAElEQVRg28w8uGudzPxFROwQETsBI4FHM/OhYvMJEfGh4vVoYFdgbROhvBd4K3BPRAC8mlpPnSRJTTERkyQNCBExHPgyMAP4GDABWJCZK4BPRMS1PVS/FjgY2JFaDxkR8S7gfcDfZuYzEfETYHiXei/wytEj67cHcEVmnrYp5yRJGrocmihJGig+C1yZmSuB+6glYs2aS22Sj4OpJWUA21LrHXsmInYH9mlQ74/ADhExIiJeBUwv1t8CHBwROwBExOv6MFRSkiQTMUlS+4uIvwbeD3ylWNWnRCwzlwLbAL/PzD8Uq28ENo+IXwH/BtzVoN7zwFnA3cD1wK+L9fdTSwx/WNS/GXhD389MkjRURWa2OoY+2X777XPs2LGtDkOSJEmSAFi0aNGfMnNkX+oMuHvExo4dy8KFC1sdhiRJkiQBEBEP9rWOQxMlSZIkqWImYpIkSZJUMRMxSZIkSarYgLtHrJHnn3+ezs5O1q1b1+pQ1IThw4czatQoOjo6Wh2KJEmS1BKDIhHr7Oxkm222YezYsUREq8NRDzKTtWvX0tnZybhx41odjiRJktQSg2Jo4rp16xgxYoRJ2AAQEYwYMcLeS0mSJA1pg6JHDDAJG0BsK0nt4rzPXV/q/k/8/PRS9y9JGrgGRY+YJEmSJA0kJmKSJEmSVDETMUmSJEmq2KC5R6xef4/5r2qM/2OPPcbVV1/NrFmz+lz3bW97Gz/72c9KiArOP/98LrzwQiZPnsxVV1318vr77ruP6dOnM3v2bI499thSji1JkiQNRvaItZHHHnuMOXPm9KlOZvLSSy/1KQlbX6dZc+bMYcGCBa9IwgAmTpzI3LlzufLKK5velyRJkiQTsX6zcuVKdt99dw4//HDe/OY3c/DBB/PMM88AcO655zJhwgQmTJjAV77yFQCefvppDjzwQN7ylrcwYcIE5s2bx+zZs/ntb3/LpEmTOPnkkwH45je/ydSpU5k0aRJHH300L774IitXrmT8+PHMmjWLyZMns2rVKrbeeuuXY2l0vEZ1umpU75hjjmHFihXMmDGD8847b4M6O+ywA0uXLu3fN1OSJEka5Abl0MRWWb58OZdccgnTpk3j4x//OHPmzOHd7343l112GXfffTeZyd5778073/lOVqxYwU477cQPfvADAB5//HH23ntvlixZwuLFiwFYtmwZ8+bN44477qCjo4NZs2Zx1VVXse+++7J8+XIuu+yyDXrQFi1a1PB42223Xbd1eqp30UUXceONN3Lrrbey/fbbb1Bv9uzZPPvsszz44IOMGTOmhHdVkiRJGnxMxPrR6NGjmTZtGgAzZ87k/PPPp6Ojgw996ENstdVWAHz4wx/mtttuY//99+ekk07i1FNPZfr06bzjHe/g0UcffcX+brnlFhYtWsRee+0FwJ///Gd22GEH9t13X8aMGcM+++yzQQy33357w+PNmDGj2zo91dtzzz27Pd8bb7zx5Z69pUuXMmbMGK677jp+8IMf8PDDD3Pcccex33779fFdlKSavt7v6zO7JEkDiUMT+1HXBxVHBJnZsOxuu+3GokWLmDhxIqeddhpnnXXWBmUyk8MPP5zFixezePFili9fzplnngnwcsLUqE53uqvTW71G1q1bxymnnMKcOXOYOHEiS5YsAeCDH/wg//mf/8nll1/OvHnz+rRPSZIkaagwEetHDz30EHfeeScA11xzDW9/+9vZd999ue6663jmmWd4+umn+d73vsc73vEOVq9ezZZbbsnMmTM56aSTuPfee9lmm2148sknX97fe9/7Xq699loefvhhAB555BEefPDBHmPo7ni96Wu9s88+m8MOO4yxY8e+IhGr337cccf1elxJkiRpKBqUQxNbNTxl/PjxXHHFFRx99NHsuuuuHHvssWy55ZYcccQRTJ06FYAjjzySPffck5tuuomTTz6ZzTbbjI6ODi688EJGjBjBtGnTmDBhAgcccADnnHMOZ599Nvvttx8vvfQSHR0dXHDBBey4447dxjB58uSGx1u5cmWPsXdXr5Hly5dz8803c8cddwC12RO/8IUvALWetdmzZ3PAAQcwefLkPr1/kiRJ0lARfR2S1qedR+wPfBUYBnw9M7/YTbmDgW8De2Xmwp72OWXKlFy48JVFli1bxvjx4/sn6I20cuVKpk+fvkHP0FBz/vnnc8UVV7DXXnsxadIkjjnmmIbl2qHNJLW3Ku4R6+/nTnblfWuSNDRExKLMnNKXOqX1iEXEMOAC4P1AJ3BPRMzPzPu7lNsGOAG4u6xYVJ0TTjiBE044odVhSJIkSW2tzHvEpgIPZOaKzHwOmAsc1KDcvwFfBtaVGEvpxo4dO+R7wyRJkiQ1p8xEbGeg/qnBncW6l0XEnsDozOxxbEhEHBURCyNi4Zo1a/o/UkmSJEmqUJmJWDRY9/INaRGxGXAe8K+97SgzL87MKZk5ZeTIkf0YoiRJkiRVr8xErBMYXbc8Clhdt7wNMAH4SUSsBPYB5kdEn25yW6/MSUfUv2wrSZIkDXVlJmL3ALtGxLiI2AI4BJi/fmNmPp6Z22fm2MwcC9wFzOht1sRGhg8fztq1a73AHwAyk7Vr1zJ8+PBWhyJJkiS1TGmzJmbmCxFxPHATtenrL83MpRFxFrAwM+f3vIfmjRo1is7OTrx/bGAYPnw4o0aNanUYkgaZsqeilySpP5X6QOfMXAAs6LLujG7Kvmtjj9PR0cG4ceM2trokSZIkVarUREySpKFsY3rpfAi0JA0NZd4jJkmSJElqwERMkiRJkipmIiZJkiRJFTMRkyRJkqSKmYhJkiRJUsVMxCRJkiSpYiZikiRJklQxEzFJkiRJqpgPdJYkqY309SHQPgBakgYme8QkSZIkqWImYpIkSZJUMRMxSZIkSaqYiZgkSZIkVcxETJIkSZIqZiImSZIkSRUzEZMkSZKkipmISZIkSVLFTMQkSZIkqWImYpIkSZJUMRMxSZIkSaqYiZgkSZIkVcxETJIkSZIqZiImSZIkSRUzEZMkSZKkipWaiEXE/hGxPCIeiIjZDbYfExH3RcTiiLg9IvYoMx5JkiRJagelJWIRMQy4ADgA2AM4tEGidXVmTszMScCXgXPLikeSJEmS2kWZPWJTgQcyc0VmPgfMBQ6qL5CZT9QtbgVkifFIkiRJUlvYvMR97wysqlvuBPbuWigijgM+DWwBvKfRjiLiKOAogF122aXfA5UkaaA673PX96n8iZ+fXlIkkqS+KLNHLBqs26DHKzMvyMw3AacCn220o8y8ODOnZOaUkSNH9nOYkiRJklStMhOxTmB03fIoYHUP5ecCHywxHkmSJElqC2UmYvcAu0bEuIjYAjgEmF9fICJ2rVs8EPhNifFIkiRJUlso7R6xzHwhIo4HbgKGAZdm5tKIOAtYmJnzgeMj4n3A88CjwOFlxSNJkiRJ7aLMyTrIzAXAgi7rzqh7/ckyjy9JkiRJ7ajUBzpLkiRJkjbUVCIWEd+JiAMjwsRNkiRJkjZRs0MTLwQ+BpwfEd8GLs/MX5cXliRpsOnr864kSRrMmurhyswfZeZHgcnASuDmiPhZRHwsIjrKDFCSJEmSBpumhxpGxAjgCOBI4BfAV6klZjeXEpkkSZIkDVJNDU2MiO8CuwPfAP4+M/9QbJoXEQvLCk6SJEmSBqNm7xH7ejEV/csi4lWZ+WxmTikhLkmSJEkatJodmnh2g3V39mcgkiRJkjRU9NgjFhE7AjsDr46IPYEoNr0G2LLk2CRJkiRpUOptaOLfUZugYxRwbt36J4HPlBSTJEmSJA1qPSZimXkFcEVE/ENmfqeimCRJkiRpUOttaOLMzPwmMDYiPt11e2ae26CaJEmSJKkHvQ1N3Kr4d+uyA5EkSZKkoaK3oYlfK/79fDXhSJIkSdLg1+wDnb9MbQr7PwM3Am8BPlUMW5QkDUHnfe76VocgSdKA1exzxPbLzCeA6UAnsBtwcmlRSZIkSdIg1mwi1lH8+wHgmsx8pKR4JEmSJGnQa2poIvD9iPg1taGJsyJiJLCuvLAkSZIkafBqqkcsM2cDfwtMyczngaeBg8oMTJIkSZIGq2Z7xADGU3ueWH2dK/s5HkmSJEka9JqdNfEbwJuAxcCLxerEREySpAGlr7Ndnvj56SVFIklDW7M9YlOAPTIzywxGkiRJkoaCZmdNXALsWGYgkiRJkjRUNNsjtj1wf0T8HHh2/crMnFFKVJIkSZI0iDWbiJ1ZZhCSJEmSNJQ0O339T4GVQEfx+h7g3t7qRcT+EbE8Ih6IiNkNtn86Iu6PiF9FxC0RMaaP8UuSJEnSgNNUIhYR/wJcC3ytWLUzcF0vdYYBFwAHAHsAh0bEHl2K/YLas8neXOz/y82HLkmSJEkDU7OTdRwHTAOeAMjM3wA79FJnKvBAZq7IzOeAuXR5CHRm3pqZzxSLdwGjmg1ckiRJkgaqZhOxZ4tkCoDioc69TWW/M7CqbrmzWNedTwA3NNoQEUdFxMKIWLhmzZomQ5YkSZKk9tRsIvbTiPgM8OqIeD/wbeD7vdSJBusaJm8RMZPas8rOabQ9My/OzCmZOWXkyJFNhixJkiRJ7anZRGw2sAa4DzgaWAB8tpc6ncDouuVRwOquhSLifcDpwIzMfLbrdkmSJEkabJqavj4zX4qI64DrMrPZsYH3ALtGxDjg98AhwD/XF4iIPalNALJ/Zj7cfNiSJEmSNHD12CMWNWdGxJ+AXwPLI2JNRJzR244z8wXgeOAmYBnwrcxcGhFnRcT6B0GfA2wNfDsiFkfE/E06G0mSJEkaAHrrEfsUtdkS98rM3wFExBuBCyPixMw8r6fKmbmA2jDG+nVn1L1+30ZFLUmSKnHe567vc50TPz+9hEgkaXDp7R6xw4BD1ydhAJm5AphZbJMkSZIk9VFviVhHZv6p68riPrGOckKSJEmSpMGtt0TsuY3cJkmSJEnqRm/3iL0lIp5osD6A4SXEI0mSJEmDXo+JWGYOqyoQSZIkSRoqmn2gsyRJkiSpnzT1QGdJ0uC2MVOUS5KkjWePmCRJkiRVzERMkiRJkipmIiZJkiRJFTMRkyRJkqSKmYhJkiRJUsVMxCRJkiSpYiZikiRJklQxEzFJkiRJqpiJmCRJkiRVzERMkiRJkipmIiZJkiRJFdu81QFIkvrfeZ+7vtUhSJKkHtgjJkmSJEkVMxGTJEmSpIqZiEmSJElSxUzEJEmSJKliJmKSJEmSVDETMUmSJEmqmImYJEmSJFWs1OeIRcT+wFeBYcDXM/OLXbbvC3wFeDNwSGZeW2Y8kiSpfH19jt2Jn59eUiSS1L5K6xGLiGHABcABwB7AoRGxR5diDwFHAFeXFYckSZIktZsye8SmAg9k5gqAiJgLHATcv75AZq4str1UYhySJEmS1FbKvEdsZ2BV3XJnsa7PIuKoiFgYEQvXrFnTL8FJkiRJUquUmYhFg3W5MTvKzIszc0pmThk5cuQmhiVJkiRJrVVmItYJjK5bHgWsLvF4kiRJkjQglJmI3QPsGhHjImIL4BBgfonHkyRJkqQBobRELDNfAI4HbgKWAd/KzKURcVZEzACIiL0iohP4R+BrEbG0rHgkSZIkqV2U+hyxzFwALOiy7oy61/dQG7IoSZIkSUNGqYmYJKl/9PUBuZIkqb2VeY+YJEmSJKkBEzFJkiRJqpiJmCRJkiRVzHvEJElSS/X1HsgTPz+9pEgkqTr2iEmSJElSxUzEJEmSJKliJmKSJEmSVDETMUmSJEmqmImYJEmSJFXMREySJEmSKmYiJkmSJEkV8zlikiRpQOnrc8fAZ49Jaj8mYpLUAhtzISlJkgYPhyZKkiRJUsVMxCRJkiSpYiZikiRJklQxEzFJkiRJqpiJmCRJkiRVzERMkiRJkipmIiZJkiRJFfM5YpK0iXwmmDT49PVz7QOjJfWViZgkSRr0/MJEUrsxEZOkLrxgkyRJZfMeMUmSJEmqWKmJWETsHxHLI+KBiJjdYPurImJesf3uiBhbZjySJEmS1A5KG5oYEcOAC4D3A53APRExPzPvryv2CeDRzPyriDgE+BLwkbJikjQ0OdRQUtmc3ENSX5V5j9hU4IHMXAEQEXOBg4D6ROwg4Mzi9bXAf0REZGaWGJckSVJLVfEFkcme1N7KTMR2BlbVLXcCe3dXJjNfiIjHgRHAn+oLRcRRwFHF4lMRsbyUiFtre7qct9qa7TWw2F4Di+01sNheberTZ22wyrYaWGyvgeWv+1qhzEQsGqzr2tPVTBky82Lg4v4Iql1FxMLMnNLqONQc22tgsb0GFttrYLG9Bg7bamCxvQaWiFjY1zplTtbRCYyuWx4FrO6uTERsDmwLPFJiTJIkSZLUcmUmYvcAu0bEuIjYAjgEmN+lzHzg8OL1wcCPvT9MkiRJ0mBX2tDE4p6v44GbgGHApZm5NCLOAhZm5nzgEuAbEfEAtZ6wQ8qKZwAY1EMvByHba2CxvQYW22tgsb0GDttqYLG9BpY+t1fYASVJkiRJ1Sr1gc6SJEmSpA2ZiEmSJElSxUzEWiAiLo2IhyNiSZf1/zsilkfE0oj4cqvi0ys1aq+ImBQRd0XE4ohYGBFTWxmjaiJidETcGhHLis/RJ4v1r4uImyPiN8W/27U6VvXYXudExK8j4lcR8b2IeG2rY1X37VW3/aSIyIjYvlUx6i96ai+vN9pPD38Pvd5oQxExPCJ+HhG/LNrr88X6cRFxd3G9Ma+YsLD7/XiPWPUiYl/gKeDKzJxQrHs3cDpwYGY+GxE7ZObDrYxTNd201w+B8zLzhoj4AHBKZr6rhWEKiIg3AG/IzHsjYhtgEfBB4Ajgkcz8YkTMBrbLzFNbGKrosb1GUZtF94WI+BKA7dV63bVXZt4fEaOBrwO7A2/NTB9C22I9fL5ej9cbbaeH9voKXm+0nYgIYKvMfCoiOoDbgU8Cnwa+m5lzI+Ii4JeZeWF3+7FHrAUy87/Z8HlpxwJfzMxnizL+UWwT3bRXAq8pXm/Lhs/IUwtk5h8y897i9ZPAMmBn4CDgiqLYFdT+c1OLdddemfnDzHyhKHYXtcRMLdbD5wvgPOAUan8b1QZ6aC+vN9pQD+3l9UYbypqnisWO4ieB9wDXFut7vd4wEWsfuwHvKLozfxoRe7U6IPXoU8A5EbEK+HfgtBbHoy4iYiywJ3A38PrM/APU/rMDdmhdZGqkS3vV+zhwQ9XxqGf17RURM4DfZ+YvWxqUutXl8+X1Rpvr0l5eb7SpiBgWEYuBh4Gbgd8Cj9V9kdjJX76sashErH1sDmwH7AOcDHyr6PZUezoWODEzRwMnUnsmntpERGwNfAf4VGY+0ep41LPu2isiTgdeAK5qVWzaUH17UWuf04EzWhqUutXg8+X1Rhtr0F5eb7SpzHwxMydRG7UxFRjfqFhP+zARax+d1MaUZmb+HHgJ8Ibn9nU48N3i9bepfQDVBoqx2t8BrsrM9W30x2L8/fpx+A7FaRPdtBcRcTgwHfhoejNz22jQXm8CxgG/jIiV1C5I7o2IHVsXpdbr5vPl9Uab6qa9vN5oc5n5GPATal9uvDYiNi82jaKXoaQmYu3jOmrjSomI3YAtAG92bl+rgXcWr98D/KaFsahQfKt7CbAsM8+t2zSf2n9mFP/+V9WxaUPdtVdE7A+cCszIzGdaFZ9eqVF7ZeZ9mblDZo7NzLHULvInZ+b/tDBU0ePfQ6832lAP7eX1RhuKiJHrZ/SNiFcD76N2X9+twMFFsV6vN5w1sQUi4hrgXdS+gfoj8DngG8ClwCTgOeCkzPxxq2LUX3TTXsuBr1Ib4rEOmJWZi1oVo2oi4u3AbcB91L7lBfgMtXH23wJ2AR4C/jEzu07Aoor10F7nA68C1hbr/v/27j3Kzrq+9/j7A6QEQkQIIGI4JLRR0AQCxpAai3K8AEKh2pRipQelHsTUZdeiciSnrUW0y1urHE4BZUkRFQUaK82BSETEHpBrIpGrWcY0QsAa5OKFLCyX7/lj73CGyZ7bYvYzs2fer7VmzXP5/Z79nfmtJ8988vz2s2+pqtOar1B9DTReVbWyT5uNwAKfmjj2Bjm/vo1/b4w7g4zXL/HvjXEnyUG0HsaxPa0bW1dU1dlJ9gcuA3YH7gBO2vpgnI7HMYhJkiRJUrOcmihJkiRJDTOISZIkSVLDDGKSJEmS1DCDmCRJkiQ1zCAmSZIkSQ0ziEmSJElSwwxikiRJktQwg5gkqWckmZfkJ0ne1+XX+XU3jy9JkkFMktQzquou4ETgv411LZIkvRAGMUlSr9kMvGq4jZN8MsnSPutnJfnL9vKVSdYkuSfJqR36zkpyd5/1DyY5q718UpLbkqxN8vkk27+QH0qSNLkYxCRJveYTwI5J9htm+8uAP+6zfgLwz+3lU6rq1cAC4ANJZgzngEkObB9zcVXNB54B3jnMeiRJYoexLkCSpOFKchQwDbia1l2xnyT5A+AYYC/gvKr6Vt8+VXVHkr2S7APsCTxWVfe3d38gydvay/sCc4BHhlHKG4FXA7cnAdiJ1p06SZKGxSAmSeoJSaYCnwKOA94NzAVWVtWVwJVJdgP+HvhWh+7LgSXA3rTukJHkDcCbgN+tqi1JvgtM7dfvaZ4/e2Tr/gCXVNWyF/6TSZImI6cmSpJ6xV8DX6qqjcBdtIJY//3nDdD3MloP+VhCK5QB7Err7tiWJAcAizr0+xmwV5IZSXYEjm1vvw5YkmQvgCS7j2CqpCRJBjFJ0viX5BXAm4Fz2pueC2Jp+STwzar6fqf+VXUPMB14sKp+2t58DbBDkjuBjwK3dOj3FHA2cCtwFfDD9vZ7aQW/b7X7Xwu8dBR+VEnSJJGqGusaRmSPPfaoWbNmjXUZktTzfvTgXc9bn/OyeWNUycD61wjjr85eqFGS1F1r1qz5eVXtOZI+PfcesVmzZrF69eqxLkOSet6Ry/Z/3vqqj4+/f1v71wjjr85eqFGS1F1JfjLSPk5NlCRJkqSGGcQkSZIkqWEGMUmSJElqWM+9R6yTp556ik2bNvHkk0+OdSkagalTpzJz5kymTJky1qVIkiRJjZoQQWzTpk1Mnz6dWbNmkWSsy9EwVBWPPPIImzZtYvbs2WNdjiRJktSoCTE18cknn2TGjBmGsB6ShBkzZngXU5IkSZPShAhigCGsBzlmkiRJmqwmTBCTJEmSpF5hEJMkSZKkhhnEJEmSJKlhE+Kpif0duWz/UT3eqo9vGNXjDeTxxx/nq1/9KkuXLh1x39e+9rXcdNNNXagKzj33XC644AIOPfRQLr300q68hiRJkjSZeEdsHHn88cc5//zzR9Snqnj22WdHFMK29hmu888/n5UrVxrCJEmSpFFiEBslGzdu5IADDuDkk0/moIMOYsmSJWzZsgWAz3zmM8ydO5e5c+dyzjnnAPDEE09wzDHHcPDBBzN37lwuv/xyzjzzTH784x8zf/58zjjjDAC+8pWvsHDhQubPn8973/tennnmGTZu3MiBBx7I0qVLOfTQQ3nggQfYZZddnqul0+t16tNfp36nnXYaGzZs4LjjjuOzn/3s89ovX76cRYsWcfDBB/O6172Ohx9+GIBFixaxceNGAB588EEWLFgwir9pSZIkqfdNyKmJY2XdunVcdNFFLF68mFNOOYXzzz+fI444gosvvphbb72VquKwww7j9a9/PRs2bGCfffbh6quvBuAXv/gFhx12GHfffTdr164F4L777uPyyy/ne9/7HlOmTGHp0qVceumlHH744axbt46LL754mztoa9as6fh6u+2224B9Buv3uc99jmuuuYbrr7+ePfbY43l9jjjiCJYsWQLARz7yEa644gqWLl3K/fffz3777QfAnXfeybx580b9dy1JkiT1Mu+IjaJ9992XxYsXA3DSSSdx4403cuONN/K2t72NadOmscsuu/D2t7+dG264gXnz5vHtb3+bD33oQ9xwww3suuuu2xzvuuuuY82aNbzmNa9h/vz5XHfddWzY0Hq/2n777ceiRYu26TPQ6w3WZ6h+A/niF7/IwoULOfjggzn//POZOnUq69evZ/bs2c99RphBTJIkSdqWQWwU9f+A4iRUVce2L3/5y1mzZg3z5s1j2bJlnH322du0qSpOPvlk1q5dy9q1a1m3bh1nnXUWANOmTet43IFeb7A+Q/Xr5Etf+hK33XYb3/nOd/jBD37AK17xCl71qldx1113PS94rV69moMOOmhEx5YkSZImOoPYKLr//vu5+eabAfja177G6173Og4//HCuvPJKtmzZwhNPPME3vvENfu/3fo+HHnqInXfemZNOOokPfvCDfP/732f69On86le/eu54b3zjG1m+fDmbN28G4NFHH+UnP/nJoDUM9HpDGWm/u+66i9e+9rXssssufP3rX+emm25i3rx5PProo+y0005Aa2rl1Vdf7R0xSZIkqZ8J+R6xph4339+BBx7IJZdcwnvf+17mzJnD+973PnbeeWfe9a53sXDhQgDe8573cMghh7Bq1SrOOOMMtttuO6ZMmcIFF1zAjBkzWLx4MXPnzuXoo4/m05/+NB/72Md4y1vewrPPPsuUKVM477zz2HvvvQes4dBDD+34elsfnjHSfgM5+eSTOf7441m+fDlvfetb2X///Zk2bRpHHnkk5557LieccAJz585lxowZvOQlLxnhb1KSJEma2DLSKWljbcGCBbV69ernbbvvvvs48MADx6iilo0bN3Lsscdy9913j2kdvWY8jJ00WfX/zMWx+k+swXT6XMjxVmcv1ChJ6q4ka6pqRI8Kd2qiJEmSJDXMIDZKZs2a5d0wSZIkScNiEJMkSZKkhhnEJEmSJKlhEyaI9dpDR+SYSZIkafKaEEFs6tSpPPLII/5h30OqikceeYSpU6eOdSmSJElS4ybE54jNnDmTTZs28fDDD491KRqBqVOnMnPmzLEuQ5IkSWrchAhiU6ZMYfbs2VvYB9QAABLISURBVGNdhiRJkiQNS1enJiY5Ksm6JOuTnDlAmxOS3JvkniRf7WY9kiRJkjQedO2OWJLtgfOANwObgNuTrKiqe/u0mQMsAxZX1WNJ9upWPZIkSZI0XnTzjthCYH1Vbaiq/wQuA47v1+a/A+dV1WMAVbW5i/VIkiRJ0rjQzSD2MuCBPuub2tv6ejnw8iTfS3JLkqM6HSjJqUlWJ1ntAzkkSZIk9bpuBrF02Nb/+fI7AHOANwDvAL6Q5MXbdKq6sKoWVNWCPffcc9QLlSRJkqQmdTOIbQL27bM+E3ioQ5t/raqnqurfgXW0gpkkSZIkTVjdDGK3A3OSzE7yW8CJwIp+ba4EjgBIsgetqYobuliTJEmSJI25rgWxqnoaeD+wCrgPuKKq7klydpLj2s1WAY8kuRe4Hjijqh7pVk2SJEmSNB509QOdq2olsLLftg/3WS7g9PaXJEmSJE0KXf1AZ0mSJEnStgxikiRJktQwg5gkSZIkNcwgJkmSJEkNM4hJkiRJUsMMYpIkSZLUMIOYJEmSJDXMICZJkiRJDTOISZIkSVLDDGKSJEmS1DCDmCRJkiQ1zCAmSZIkSQ0ziEmSJElSwwxikiRJktQwg5gkSZIkNcwgJkmSJEkNM4hJkiRJUsMMYpIkSZLUsGEFsSRfT3JMEoObJEmSJL1Aww1WFwB/AvwoySeSHNDFmiRJkiRpQhtWEKuqb1fVO4FDgY3AtUluSvLuJFO6WaAkSZIkTTTDnmqYZAbwLuA9wB3A/6IVzK7tSmWSJEmSNEHtMJxGSf4FOAD4MvD7VfXT9q7Lk6zuVnGSJEmSNBENK4gBX6iqlX03JNmxqn5TVQu6UJckSZIkTVjDnZr4sQ7bbh7NQiRJkiRpshj0jliSvYGXATslOQRIe9eLgJ27XJskSZIkTUhDTU08ktYDOmYCn+mz/VfA/+xSTZIkSZI0oQ0axKrqEuCSJH9YVV9vqCZJkiRJmtCGmpp4UlV9BZiV5PT++6vqMx26SZIkSZIGMdTUxGnt77t0uxBJkiRJmiyGmpr4+fb3jzRTjiRJkiRNfENNTTx3sP1V9YHRLUeSJEmSJr6hpiauaaQKSZIkSZpEhvPUREmSJEnSKNpusJ1Jzml//z9JVvT/GurgSY5Ksi7J+iRnDtJuSZJKsmDkP4IkSZIk9ZahpiZ+uf3970d64CTbA+cBbwY2AbcnWVFV9/ZrNx34AHDrSF9DkiRJknrRoHfEqmpN+/u/ATcDjwGPAje3tw1mIbC+qjZU1X8ClwHHd2j3UeBTwJMjrF2SJEmSetKgQWyrJMcAPwbOBf4RWJ/k6CG6vQx4oM/6pva2vsc9BNi3qq4a4vVPTbI6yeqHH354OCVLkiRJ0rg11NTErf4BOKKq1gMk+W3gauCbg/RJh2313M5kO+CzwLuGevGquhC4EGDBggU1RHNJkiRJGteGdUcM2Lw1hLVtADYP0WcTsG+f9ZnAQ33WpwNzge8m2QgsAlb4wA5JkiRJE91QH+j89vbiPUlWAlfQuqv1R8DtQxz7dmBOktnAg8CJwJ9s3VlVvwD26PNa3wU+WFWrR/gzSJIkSVJPGWpq4u/3Wf4Z8Pr28sPAboN1rKqnk7wfWAVsD/xTVd2T5GxgdVUN+fh7SZIkSZqIhvpA53e/kINX1UpgZb9tHx6g7RteyGtJkiRJUq8Y1sM6kkwF/gx4FTB16/aqOqVLdUmSJEnShDXch3V8GdgbOBL4N1oP3vhVt4qSJEmSpIlsuEHsd6rqb4AnquoS4BhgXvfKkiRJkqSJa7hB7Kn298eTzAV2BWZ1pSJJkiRJmuCG+4HOFybZDfgbYAWwS3tZkiRJkjRCwwpiVfWF9uK/Aft3rxxJkiRJmviGNTUxyYwk/zvJ95OsSXJOkhndLk6SJEmSJqLhvkfsMmAz8IfAEuDnwOXdKkqSJEmSJrLhvkds96r6aJ/1jyX5g24UJEmSJEkT3XDviF2f5MQk27W/TgCu7mZhkiRJkjRRDXpHLMmvgAICnA58pb1rO+DXwN92tTpJkiRJmoAGDWJVNb2pQiRJkiRpshjue8RIchxweHv1u1V1VXdKkiRJkqSJbbiPr/8E8BfAve2vv2hvkyRJkiSN0HDviL0VmF9VzwIkuQS4AzizW4VJkiRJ0kQ13KcmAry4z/Kuo12IJEmSJE0Ww70j9nHgjiTX03qC4uHAsq5VJUmSJEkT2JBBLEmAG4FFwGtoBbEPVdV/dLk2SZIkSZqQhgxiVVVJrqyqVwMrGqhJknrakcv232bbqo9vGINKJEnSeDXc94jdkuQ1Xa1EkiRJkiaJ4b5H7AjgtCQbgSdoTU+sqjqoW4VJkiRJ0kQ13CB2dFerkCRJkqRJZNAglmQqcBrwO8BdwEVV9XQThUmSJEnSRDXUe8QuARbQCmFHA//Q9YokSZIkaYIbamriK6tqHkCSi4Dbul+SJEmSJE1sQ90Re2rrglMSJUmSJGl0DHVH7OAkv2wvB9ipvb71qYkv6mp1kiRJkjQBDRrEqmr7pgqRJEmSpMliuB/oLEmSJEkaJQYxSZIkSWqYQUySJEmSGmYQkyRJkqSGGcQkSZIkqWFdDWJJjkqyLsn6JGd22H96knuT3JnkuiT7dbMeSZIkSRoPuhbEkmwPnAccDbwSeEeSV/ZrdgewoKoOApYDn+pWPZIkSZI0XnTzjthCYH1Vbaiq/wQuA47v26Cqrq+qLe3VW4CZXaxHkiRJksaFbgaxlwEP9Fnf1N42kD8DvtlpR5JTk6xOsvrhhx8exRIlSZIkqXndDGLpsK06NkxOAhYAn+60v6ourKoFVbVgzz33HMUSJUmSJKl5O3Tx2JuAffuszwQe6t8oyZuAvwJeX1W/6WI9kiRJkjQudPOO2O3AnCSzk/wWcCKwom+DJIcAnweOq6rNXaxFkiRJksaNrgWxqnoaeD+wCrgPuKKq7klydpLj2s0+DewC/HOStUlWDHA4SZIkSZowujk1kapaCazst+3DfZbf1M3XlyRJkqTxqKsf6CxJkiRJ2pZBTJIkSZIaZhCTJEmSpIYZxCRJkiSpYQYxSZIkSWqYQUySJEmSGmYQkyRJkqSGGcQkSZIkqWEGMUmSJElqmEFMkiRJkhpmEJMkSZKkhhnEJEmSJKlhBjFJkiRJaphBTJIkSZIaZhCTJEmSpIYZxCRJkiSpYQYxSZIkSWqYQUySJEmSGmYQkyRJkqSGGcQkSZIkqWEGMUmSJElqmEFMkiRJkhpmEJMkSZKkhhnEJEmSJKlhBjFJkiRJaphBTJIkSZIaZhCTJEmSpIYZxCRJkiSpYQYxSZIkSWqYQUySJEmSGmYQkyRJkqSGGcQkSZIkqWEGMUmSJElqmEFMkiRJkhrW1SCW5Kgk65KsT3Jmh/07Jrm8vf/WJLO6WY8kSZIkjQc7dOvASbYHzgPeDGwCbk+yoqru7dPsz4DHqup3kpwIfBL4427VJKm3HLls/222rfr4hjGoRJIkaXR1LYgBC4H1VbUBIMllwPFA3yB2PHBWe3k58I9JUlXVxboksW3IMeBIkiQ1J93KPEmWAEdV1Xva638KHFZV7+/T5u52m03t9R+32/y837FOBU5tr84F7u5K0RptewA/H7KVxprj1Bscp97hWPUGx6l3OFa9YbKP035VtedIOnTzjlg6bOuf+obThqq6ELgQIMnqqlrwwstTtzlWvcFx6g2OU+9wrHqD49Q7HKve4DiNXDcf1rEJ2LfP+kzgoYHaJNkB2BV4tIs1SZIkSdKY62YQux2Yk2R2kt8CTgRW9GuzAji5vbwE+I7vD5MkSZI00XVtamJVPZ3k/cAqYHvgn6rqniRnA6uragVwEfDlJOtp3Qk7cRiHvrBbNWvUOVa9wXHqDY5T73CseoPj1Dscq97gOI1Q1x7WIUmSJEnqrKsf6CxJkiRJ2pZBTJIkSZIaNm6CWJKpSW5L8oMk9yT5SHv7F5P8e5K17a/5A/Q/OcmP2l8nd2qjF26Qcbqhzxg9lOTKAfo/06dd/4e3aJQl2T7JHUmuaq/PTnJr+zy5vP0gnU79liVZn2RdkiObrXpy6jBWl7Z//3cn+ackUwbo5znVoA7j5DVqnOowVl6nxqEkG5Pc1f59r25v2z3Jte3z5dokuw3Q1/OqIQOM06eT/DDJnUm+keTFw+2rlnHzHrEkAaZV1a/bf3DcCPwFcBpwVVUtH6Tv7sBqYAGtzyFbA7y6qh7rfuWTy0DjVFW39GnzdeBfq+pLHfr/uqp2aa7iyS3J6bTOixdV1bFJrgD+paouS/I54AdVdUG/Pq8EvgYsBPYBvg28vKqeabj8SaXDWL0V+GZ791eB/9t/rNr9PKca1GGcvojXqHGp/1j12+d1apxIshFYUFU/77PtU8CjVfWJJGcCu1XVh/r187xq0ADj9BZaTzx/OsknAfqP00B91TJu7ohVy6/bq1PaX8NNiUcC11bVo+0T8FrgqC6UOekNNU5JpgP/Fej4P41qTpKZwDHAF9rroTU2W/9gvAT4gw5djwcuq6rfVNW/A+tphTJ1Sf+xAqiqle3zrYDbaH0Wo8ZQp3EaJq9RDRtsrLxO9YTjaV2jYOBrlefVGKuqb1XV0+3VW/A6NWLjJojBc9MI1gKbaZ1ct7Z3/V37tudnk+zYoevLgAf6rG9qb1MXDDJOAG8DrquqXw7QfWqS1UluSdLpH1aNnnOA/wE8216fATze5x/Ngc4Tz6fm9R+r57TvPP8pcM0AfT2nmjPQOHmNGn8GPKfwOjXeFPCtJGuSnNre9pKq+ilA+/teHfp5XjWr0zj1dQr/fxbHSPtOWuMqiFXVM1U1n1aiXphkLrAMOAB4DbA7sM0tTyCdDte1Qie5AcZpq3fQmtY2kP9SVQuAPwHOSfLbXSx10kpyLLC5qtb03dyhaafzxPOpQQOMVV/n05qWeMMA+z2nGjDIOHmNGmeGcU55nRpfFlfVocDRwJ8nOXyY/TyvmjXgOCX5K+Bp4NKR9p3sxlUQ26qqHge+CxxVVT9tz875DXAxnadIbQL27bM+E3io64VOcn3HCSDJDFrjc/UgfR5qf9/Q7ntIt+ucpBYDx7XnZV9GaxrOOcCLk2z9IPeBzhPPp2ZtM1ZJvgKQ5G+BPYHTB+rsOdWYjuPkNWpcGuyc8jo1zvT5fW8GvkFrfH6W5KUA7e+bO3T1vGrQAONE+yEpxwLvrAEePDFQX42jIJZkz61PW0myE/Am4Id9TsTQmiN8d4fuq4C3JNmt/WSdt7S3aZQNNE7t3X9E603rTw7Qd7et03aS7EHrYnlv96uefKpqWVXNrKpZwIm03kz7TuB6YEm72cnAv3bovgI4McmOSWYDc2i9R0ldMMBYnZTkPbTeA/GOquo0vcpzqkGDjJPXqHFmoLFq7/Y6NY4kmdZ+zx5JptE6N+6mdR3a+hTEga5VnlcNGWickhxFaxbAcVW1ZSR9m6l8/Bs3QQx4KXB9kjuB22m99+gq4NIkdwF3AXsAHwNIsiDJFwCq6lHgo+1+twNnt7dp9A00TtC64D1vukffcQIOBFYn+QGtQPCJqvIC16wPAacnWU/rPWMXASQ5LsnZAFV1D3AFrT8+rgH+3CcmjonPAS8Bbk7rkb8fBs+pcchrVG/xOjW+vAS4sf37vg24uqquAT4BvDnJj4A3t9c9r8bOQOP0j8B04Nr2depzAEn2SbJyiL5iHD2+XpIkSZImi/F0R0ySJEmSJgWDmCRJkiQ1zCAmSZIkSQ0ziEmSJElSwwxikiRJktQwg5gkSZIkNcwgJkmSJEkNM4hJktSW5NdjXYMkaXIwiEmSJElSwwxikqSek+RFSe5Ick+SLUnWJrklyXZ92nwyydI+62cl+cv28pVJ1rT7n9rh+LOS3N1n/YNJzmovn5TktvZrfj7J9l39YSVJE5JBTJLUc6rql1V1CPBu4Nqqml9Vi6rq2T7NLgP+uM/6CcA/t5dPqapXAwuADySZMZzXTXJg+5iLq2o+8Azwzhf440iSJqEdxroASZJegLnAPZ12VNUdSfZKsg+wJ/BYVd3f3v2BJG9rL+8LzAEeGcbrvRF4NXB7EoCdgM0voH5J0iRlEJMk9bJXAt8fZP9yYAmwN607ZCR5A/Am4HerakuS7wJT+/V7mufPGtm6P8AlVbXsBVcuSZrUnJooSepl+wD/Mcj+y4ATaYWx5e1tu9K6O7YlyQHAog79fgbslWRGkh2BY9vbrwOWJNkLIMnuSfYbhZ9DkjTJGMQkSb1sFXBRktd32llV9wDTgQer6qftzdcAOyS5E/gocEuHfk8BZwO3AlcBP2xvvxf4a+Bb7f7XAi8d1Z9IkjQppKrGugZJkiRJmlS8IyZJkiRJDTOISZIkSVLDDGKSJEmS1DCDmCRJkiQ1zCAmSZIkSQ0ziEmSJElSwwxikiRJktSw/wfsUUTaYwjVZgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1044x720 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# histgram of only samples\n",
    "ax = plt.subplot(311)\n",
    "ax.set_autoscaley_on(False)\n",
    "plt.hist(lambda_1_samples,histtype='stepfilled',\n",
    "        bins=30,alpha=.85,color=\"#A60628\",normed=True,\n",
    "        label='posterior of $\\lambda_1$')\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.title(\"Posterior distributions of the parameters \"\n",
    "         r\"$\\lambda_1, \\lambda_2, \\tau$\")\n",
    "plt.xlim([15,30])\n",
    "plt.xlabel('$\\lambda_1$ value')\n",
    "plt.ylabel('Density')\n",
    "\n",
    "ax=plt.subplot(312)\n",
    "plt.hist(lambda_2_samples,histtype='stepfilled',\n",
    "        bins=30,alpha=.85, color=\"#7A68A6\",normed=True,\n",
    "        label='posterior of $\\lambda_2$')\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.xlim([15,30])\n",
    "plt.xlabel('$\\lambda_2$ value')\n",
    "plt.ylabel('Density')\n",
    "\n",
    "ax=plt.subplot(313)\n",
    "w = 1.0 / tau_samples.shape[0] * np.ones_like(tau_samples)\n",
    "plt.hist(tau_samples, bins=n_count_data,\n",
    "         alpha=1, color=\"#467821\",normed=True,\n",
    "        label='posterior of $\\tau $',\n",
    "        weights=w,rwidth=2.)\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.ylim([0, .75])\n",
    "plt.xlim(([35, len(count_data)-20]))\n",
    "plt.xlabel(r'$\\tau$ value')\n",
    "plt.ylabel('Probability')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}