poisson_factor_test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymc3 as pm\n",
    "import arviz as az\n",
    "import numpy as np\n",
    "\n",
    "az.style.use('arviz-darkgrid')\n",
    "\n",
    "\n",
    "def make_data(N=1000, factor=.5):\n",
    "    U = np.random.exponential(10, N)\n",
    "\n",
    "    obs_nodamp = np.random.poisson(U)\n",
    "    obs_factor = np.random.poisson(U*factor)\n",
    "    \n",
    "    return obs_nodamp, obs_factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000\n",
    "\n",
    "obs_nodamp, obs_factor = make_data(N, factor=.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [f, U]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 00:10<00:00 Sampling 4 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 11 seconds.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model:    \n",
    "\n",
    "    U = pm.Exponential('U', 10, shape=N)\n",
    "    f = pm.Uniform('f', 0, 1)\n",
    "    \n",
    "    p_nodamp = pm.Poisson('p_nodamp', U,   observed=obs_nodamp)\n",
    "    p_factor = pm.Poisson('p_factor', U*f, observed=obs_factor)\n",
    "    \n",
    "    idata = pm.sample(1_000, return_inferencedata=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['U', 'f']"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(idata.posterior.data_vars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtsAAAHrCAYAAAAe4lGYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHDUlEQVR4nO3dd3xUVcL/8e8kIUASSoDQS2gJvYQWulJUEBvrWtBFxdUVFZSVB919UHdd/ekjrg1kwbY2bIgNREGkGSARMQiE0IuEFkpCAiH9/P64mSFjEkggk2mf9+uV19w5986dM+Ph+s3JuefYjDFGAAAAACpdgLsrAAAAAPgqwjYAAADgIoRtAAAAwEUI2wAAAICLELYBAAAAFyFsAwAAAC5C2AYAAABchLANAAAAuAhhGwAAAHARwjYAeInly5dr3LhxiomJUXR0tKKjo5WQkODuagEAziPI3RUAAFxYcnKyJk+erMLCQsXGxioiIkI2m00NGjRwd9UAAOdB2AYAL7Bs2TLl5eXpvvvu05QpU9xdHQBAOTGMBAC8wKFDhyRJrVq1cnNNAAAVYTPGGHdXAgBQupkzZ2rWrFml7uvbt6/ef//9Kq4RAKAiGEYCAB6sY8eOuuGGG7Rhwwb99ttviomJcfRut2nTxs21AwBcCD3bAOAFHnvsMX3xxRd69tlnNXbsWHdXBwBQTozZBgAAAFyEsA0AAAC4CGEbAAAAcBHCNgAAAOAihG0AAADARQjbAAAAgIsQtgEAAAAXIWwDAAAALkLYBgAAAFyEFSQBAAAAF6FnGwAAAHARwjYAAADgIoRtAAAAwEUI2wAAAICLELYBAAAAFyFsAwAAAC5C2AYAAABcJKi8B6alpbmyHl6lTp06OnXqlLurAR9F+4Kr0cbgarQxuJq721h4eHi5j6Vn+yIEBPC1wXVoX3A12hhcjTYGV/OmNuY9NQUAAAC8DGEbAAAAcBHCNgAAAOAihG0AAADARQjbAAAAgIsQtgEAAAAXIWwDAAAALkLYBgAAAFyEsA0AAAC4CGEbAAAAcBHCNgAAAOAihG0AAADARQjbAAAAgIsQtgEAAAAXIWwDAAAALkLYBgAAgFfIzzfakmSUn2/cXZVyI2wDAADA46WkGE180Oi+B4z++262u6tTbkHurgAAAABQFmOMvlsivfiK0dmzUliY1D+2miTvCNyEbQAAAHgkY4xemWX02QLreY/u0uP/a1OH6CClpbm3buVF2AYAAIBHmr9A+myBZLNJf55g0+3jpMBAm7urVSGEbQAAAHicteuMZs22boS8/z6bbr3Zu0K2HTdIAgAAwKPs3mP05FNGhYXSNWOkW25yd40uHmEbAAAAHuPsWaPH/te6GbJXjPTIwzbZbN7Zqy0RtgEAAOBB/vue0eHDUqNG0r/+YVNQkPcGbYmwDQAAAA+xe4/RJ59a2399yKbatb07aEuEbQAAAHiAwkKjF140KiiQhgyWBg7w/qAtEbYBAADgAb5ZLG3eItWsKT08yTeCtkTYBgAAgJudOmU0e641zd+fJ9jUsCFhGwAAAKgUH39qlJkptW0j/eEGd9emchG2AQAA4DanThl99rm1ffdd3j/7yO8RtgEAAOA2n8y35tRu304aPMjdtal8hG0AAAC4RUbGuV7tO+/w7sVrykLYBgAAgFt8Mt8oK0tq21YaPNDdtXENwjYAAACqXEbmuV7tu8bbFBDge73aEmEbAAAAbvDFl9KZM9YMJEMGu7s2rkPYBgAAQJXKzzf68itrXu1xt/hur7ZE2AYAAEAV+zFOOnZcCg+XLr/M3bVxLcI2AAAAqtSCL6xe7WvHSMHBvturLRG2AQAAUIV27zHa+KsUGCBdd41vB22JsA0AAIAq9HlRr/bgwVLDhoRtAAAAoFJkZhot+d7a/sMNvh+0JcI2AAAAqsi330nZ2VKb1lKP7u6uTdUgbAMAAMDljDH6+htrCMnY631zafbSELYBAADgctu2S/v2ScHB0ojh7q5N1SFsAwAAwOUWf2f1ag8dIoWF+UevtkTYBgAAgIvl5hot+8HaHn2V/wRtibANAAAAF1uzVsrMlCIaSDE93V2bqkXYBgAAgEt9WzSE5KorpcBAerYBAACASnHihFHCT9b2qCv9K2hLhG0AAAC40NJlUkGh1KWz1LIlYRsAAACoNOeGkPhf0JYI2wAAAHCRPXuM9uyVqlWThl3u7tq4B2EbAAAALvHDCqtXu19fqXYterYBAACASmGM0Q8rrO1hl/tn0JYI2wAAAHCBHTullBSpenVp0AB318Z9CNsAAACodPYhJP1jpZAQerYBAACASmGM0fLl1vZwPx5CIhG2AQAAUMmStkpHjko1a1o92/6MsA0AAIBKtbxoCMmggVKNGvRsAwAAAJWisNBoxUpr29+HkEiEbQAAAFSiLUnSseNSWKjUt4+7a+N+hG0AAABUmlU/WkNIBg6UgoPp2Q5ydwUAAADgegUFBVqxYoW2bt2qrVu3atu2bcrOztbo0aP1xBNPVMp7GGO0erW1PXQwQVsibAMAAPiFrKwsTZ8+3aXvsXOXdPiItZANQ0gshG0AAAA/EBQUpFGjRqlDhw7q1KmT9u/fr6effrpS32PVamsISWw/ZiGxI2wDAAD4gZo1a+rJJ590PD9y5Eilv8eqH61HhpCcQ9gGAABe6dChQxo7dqx69uypF198UXPnztXy5ct16tQptWrVSvfcc48GDx4sSfrhhx80b9487dmzRzVr1tSIESP0wAMPqEaNGk7nzM7O1ieffKIffvhBBw4ckCS1adNGY8eO1dVXX12iDhs3btSyZcuUmJio1NRU5ebmqnHjxhoyZIjGjx+vWrVqOR2/YcMGPfDAAxo9erQeeughzZkzR6tXr1ZGRoZatGihW2+9Vddcc42LvjHX2r/faN8+KSiIhWyKI2wDAACvlp+frwcffFCHDh1Sz549lZ6ero0bN+qxxx7TSy+9pN27d2vWrFnq2bOn+vXrp40bN2r+/Pk6deqUnnrqKcd5Tp48qcmTJ2vXrl2qX7++evbsKWOMNm/erH/9619KTk7W1KlTnd575syZ2rVrl9q2bavevXsrNzdX27dv1/vvv681a9bozTffVEhISIk6nz59Wvfcc4/Onj2rHj16OOr8zDPPqLCwUNddd53Lv7fKtjrOeozpKdWqRc+2HWEbAAB4tc2bN6t37976/PPPVbNmTUnSokWL9PTTT+v5559XRkaG3nzzTXXs2FGSdOzYMY0fP15Lly7VX/7yFzVr1kyS9PTTT2vXrl26+eab9cADDyg4OFiSdOLECU2dOlWfffaZBg4cqP79+zve++6771a3bt0UFhbmKMvNzdWLL76oL7/8Uh999JHuvvvuEnVevXq1Ro4cqccff9zxPqtWrdKjjz6q//73vyXC9sSJE5WYmFih72X69OkaM2ZMhV5zKezjtYcOIWgXR9gGAABeLSAgQNOmTXMEbUkaPXq0Zs2apZSUFN11112OoC1JERERuvLKK/Xxxx8rMTFRzZo1044dO7R27Vp16tRJDz30kAICzi1FUr9+fT322GO644479PnnnzuF7QEDBpSoT3BwsB5++GEtXLhQq1evLjVsh4aGaurUqY6gLUlDhw5V27ZttXv3bh06dEhNmzZ17Ovfv7+aNGlSoe+lefPmFTr+Uhw5arRtu2SzSYMHVtnbegXCNgAA8GpNmjRRy5YtncoCAgLUuHFjpaenq1+/fiVeY+/NPnHihCQpISFBkjRkyBCnoG0XHR2tkJAQbd26tcS+1NRUxcXFaf/+/Tpz5owKCwslSdWqVXOM+/69Dh06qE6dOiXKW7Rood27d+vEiRNOYXv8+PGlnsdT/Fg0hKRrF6lePXq2iyNsAwAArxYREVFqub2nu7T99n25ubmSpMOHD0uS5syZozlz5pT5Xjk5OU7PP/zwQ82ePVv5+fkVqnPDhg1LLbeP77bXy1vErbGGkAxhFpISCNsAAMCr2WznD3gX2i9ZKx9KUvfu3R293heyZcsWvfrqqwoLC9OUKVMUExOj+vXrO4aGjBkzRsePH7/oOhX33nvvad++fRV6zbXXXqsePXpU6DUXIzPTaOOv1vbAkqNq/B5hGwAA+D177/eQIUN02223les1K1eulCTdd999JaYFzM7OdgxRqQzr1q2r8A2SMTExVRK2E36SCgqkyFZSi+b0bP8eYRsAAPi9vn376vXXX9eqVavKHbYzMzMllT4kZPny5Y7e8srwn//8p9LOVdni1lqfk17t0pW8AwAAAMDPdOnSRX379tWmTZs0Y8YMnTlzpsQxO3fu1Lp16xzP7Tdlfv31105jtvfu3avXXnvN9ZX2APn5RvHx1vaggfRql4aebQAAAEn/+Mc/9PDDD2vBggVaunSp2rdvrwYNGujMmTPatWuXjh49qptvvtkx9d+YMWP04YcfKi4uTjfddJM6duyojIwMJSYmaujQoUpKSnLJkuiX4vnnn9f27dslSadOnZIkrV271ml6wrfeeqvc5/t1k3T6jFS3rtSp4wUP90uEbQAAAEn16tXTG2+8oa+++krff/+9duzYoc2bN6tevXpq2rSpbrrpJo0cOdJxfJ06dfTf//5Xs2bNUmJiouLi4tSkSRPde++9uu222/SHP/zBjZ+mdHv37lVSUpJTWXp6utLT0y/qfPZZSAb0lwID6dkujc2Uc0BRWlqaq+viNcLDw/k+4DK0L7gabQyuRhvzD8YY3TTO6PBh6f/9y1al0/65u42Fh4eX+1jGbAMAAKDC9u6VDh+WgqtJfXq7uzaei7ANAACACotbaz327iXVrMkQkrIQtgEAAFBha9cVjdceQNA+H8I2AAAAKiQ93Shpq7XdP9a9dfF0hG0AAABUyE8/S8ZIbdtIjRrSs30+hG0AAABUSHy8NYQktp+bK+IFCNsAAAAot4ICo4SfrO3+sfRqXwhhGwAAAOWWvE06lSGFhUpdOru7Np6PsA0AAIByi0+whpD06SMFBdGzfSGEbQAAAJTbunjrsX8/gnZ5ELYBAABQLidOGG3fYW336+veungLwjYAAADKxX5jZHSUVL8+PdvlQdgGAABAuawrGq/NQjblR9gGAADABeXnG61fb20z5V/5EbYBAABwQVuSpNNnpLp1pA7R7q6N9yBsAwAA4ILWFa0a2a+vFBhIz3Z5EbYBAABwQfEJ1mMsQ0gqhLANAACA8zqaarR7jxQQIPXt7e7aeBfCNgAAAM7L3qvduZNUpw492xVB2AYAAMB5xReN145l1cgKI2wDAACgTLm5Rj9vsLaZX7viCNsAAAAo06+bpLPZUv36Uvt27q6N9yFsAwAAoEzxCfYhJJLNxjCSiiJsAwAAoEzr4q3H/ozXviiEbQAAAJTq4EGj3w5IgYFS717uro13ImwDAACgVPYp/7p1lcLC6Nm+GIRtAAAAlGpd0Xjt/qwaedEI2wAAACghO9vol0RrO7afe+vizQjbAAAAKCFxo5SbKzVqJLWOdHdtvBdhGwAAACWsK1o1sj9T/l0SwjYAAACcGGMcU/7FMl77khC2AQAA4OS336TDR6Rq1aRePd1dG+9G2AYAAICTdUVT/vXsIdWsSc/2pSBsAwAAwIl9vHYsq0ZeMsI2AAAAHLKyjH7dZG33j3VvXXwBYRsAAAAOP2+Q8vOl5s2kFs3p2b5UhG0AAAA42FeNZCGbykHYBgAAgCRryr94+5R/jNeuFIRtAAAASJJ275GOHZeqV7dmIsGlI2wDAABAkhwL2fTqKVWvTs92ZSBsAwAAQJKU8JM1XrsfQ0gqDWEbAAAAysw02rzZ2ubmyMpD2AYAAIDWb5AKCqVWLaVmTenZriyEbQAAACjePuUfC9lUKsI2AACAnyssNEpIsLZj+9KrXZkI2wAAAH5u5y7pxEmpZg2pezd318a3ELYBAAD8nH3Kv969pOBgerYrE2EbAADAz9mn/IuNJWhXNsI2AACAHzt1yihpq7XNlH+Vj7ANAADgx376WSoslNq0lho1pGe7shG2AQAA/Fh8PFP+uRJhGwAAwE8VFhol/GRtM+WfaxC2AQAA/NS27VL6KSkkROrW1d218U2EbQAAAD+1rmgISd/eUlAQPduuQNgGAADwU/H2ISRM+ecyhG0AAAA/lJZmtG2btR3b17118WWEbQAAAD+UsF4yRmrfTmrQgJ5tVyFsAwAA+KF1TPlXJQjbAAAAfqagwOin9dZ2/370arsSYRsAAMDPbE2WMjOlsDCpU0d318a3EbYBAAD8jH0ISb8+TPnnaoRtAAAAPxOfYD0y5Z/rEbYBAAD8yPETRjt2WttM+ed6hG0AAAA/klC0kE2HaCk8nJ5tVyNsAwAA+BH7eO3+TPlXJQjbAAAAfiI/3+jnn63tWKb8qxKEbQAAAD+xJUk6fUaqU9saRgLXI2wDAAD4CceUf32lwEB6tqsCYRsAAMBPMOVf1SNsAwAA+IGjqUa790g2m7WYDaoGYRsAAMAPrF1nPXbuJNWpQ892VSFsAwAA+IE1a63x2gMHELSrEmEbAADAx2VlGW34xdoeOMC9dfE3hG0AAAAft36DlJcnNW0qtY50d238C2EbAADAxzmGkPSXbDaGkVQlwjYAAIAPKygwjpsjGa9d9QjbAAAAPix5m5SeLoWFSj26u7s2/oewDQAA4MPsQ0j69ZWCgujZrmqEbQAAAB8Wt8Z6ZAiJexC2AQAAfNShw0Z790mBAVJsP3fXxj8RtgEAAHzUmrXWY9euUu3a9Gy7A2EbAADAR8WtYdVIdyNsAwAA+KDTp402/mpts2qk+xC2AQAAfFDCeqmgQGrZQmrZgp5tdyFsAwAA+KA1jiEkbq6InyNsAwAA+Jj8fKN1CdY247Xdi7ANAADgY7YkSZmZUu3aUpfO7q6NfyNsAwAA+Bj7qpH9+7FqpLsRtgEAAHxMXNH82gMYQuJ2hG0AAAAf8tsBowMHpMBAqV8fd9cGhG0AAAAfYl81smcPKSyMnm13I2wDAAD4kB/jWDXSkxC2AQAAfMSJE0abt1jbQwa7ty6wELYBAAB8xI9xkjFSx45So4b0bHsCwjYAAICPWPWjNYRk6GCCtqcgbAMAAPiAjAyjXxKtbYaQeA7CNgAAgA9Ys04qKJDatJZatqBn21MQtgEAAHzA6tVFQ0iGuLkicELYBgAA8HJZWUYJ663tIYzX9iiEbQAAAC8X/5OUmys1ayq1a+vu2qA4wjYAAICXW1VsCInNRs+2JyFsAwAAeLGcHKO166xthpB4HsI2AACAF9vwi3T2rNSggdSpo7trg98jbAMAAHgx+0I2QwZJAQH0bHsawjYAAICXys83iouztocOIWh7IsI2AACAl9q0WTqVIdWpLXXv5u7aoDSEbQAAAC+1cpU1hGTQICkoiJ5tT0TYBgAA8EKFhUar7UNImIXEYxG2AQAAvNDWZOn4cSkkROoV4+7aoCyEbQAAAC+0umgWkv6xUvXq9Gx7KsI2AACAlzHGaNVqa5tZSDwbYRsAAMDL7NwlHTwkBQdLsX3dXRucD2EbAADAyyxfcW4ISUgIPduejLANAADgRYwx+mGFtT3scoK2pyNsAwAAeJFt26XDh6UaNaQBse6uDS6EsA0AAOBFflhuDSEZ2F+qWZOebU9H2AYAAPAShYVGy1da28OGEbS9AWEbAADASyRtlVJTrYVsmIXEOxC2AQAAvIR9FpLBA1nIxlsQtgEAALxAYaHRilXW9uXMQuI1CNsAAABeYNNm6fhxKSxU6tvb3bVBeRG2AQAAvMDSZdYQkiFDpOBgera9BWEbAADAw+XmGq1YaW1fOZKg7U0I2wAAAB4uPkHKzJQaNJB6dHd3bVARhG0AAAAPt/R7awjJyOFSYCA9296EsA0AAODBMjON1q6ztq8YQdD2NoRtAAAAD7ZqtZSbJ7WOlNq1c3dtUFGEbQAAAA9mn4XkipE22Wz0bHsbwjYAAICHOppqlLjR2h453K1VwUUibAMAAHioZT9Ixkjdu0mNG9Or7Y0I2wAAAB7IGKNvl5wbQgLvRNgGAADwQMnbpH37pOBgafjl7q4NLhZhGwAAwAMt/s7q1R46RAoLo2fbWxG2AQAAPExOjtGyH6zt0VcRtL0ZYRsAAMDDxK2RTp+WGjaUYnq6uza4FIRtAAAAD2MfQjLqSpZn93aEbQAAAA9y7JjR+p+t7VFXErS9HWEbAADAg3y3VCostObWbt6csO3tCNsAAAAewhhzbggJN0b6BMI2AACAh0jcKB04INWsKQ27zN21QWUgbAMAAHiIrxcVrRg5QgoJoWfbFxC2AQAAPEB6utGq1db2tdcQtH0FYRsAAMADfLdUysuToqOk6CjCtq8gbAMAALiZMUZfL7SGkNCr7VsI2wAAAG726ybptwNSzRrSyOHurg0qE2EbAADAzey92iO4MdLnELYBAADcKD3daOUqa/u6MQRtX0PYBgAAcKOF30i59hsjo91dG1Q2wjYAAICb5OcbffGVNYTkxrE22Wz0bPsawjYAAICbrFkrpaZKdetIwy53d23gCoRtAAAAN/nsc6tX+5prpOrV6dX2RYRtAAAAN9i12yhxoxQYIN1wLUHbVxG2AQAA3GDBF1av9uDBUsOGhG1fRdgGAACoYhmZRku/t7ZvHEvQ9mWEbQAAgCr21ddSTo7Utq3UvZu7awNXImwDAABUodxco88WWENIbrmJ6f58HWEbAACgCi39XjpxUopoII0Y5u7awNUI2wAAAFWksNDoo0+sXu2b/mhTtWr0avs6wjYAAEAVWRsv7f9NCg2Vrh3j7tqgKhC2AQAAqsiHH1m92tddK4WG0qvtDwjbAAAAVWBLktGmzVJQkHTTHwja/oKwDQAAUAU++NDq1b5ipNSgAWHbXxC2AQAAXGznTqO4NZLNJt1+K0HbnxC2AQAAXOzd961e7eHDpJYtCdv+hLANAADgQnv2GK1cbfVqj7+doO1vCNsAAAAuZO/VHjpEatOasO1vCNsAAAAusn+/0fKV1vYdfyJo+yPCNgAAgIu8+76RMdLggVL7doRtf0TYBgAAcIE9e4y+/8HavmM8QdtfEbYBAABc4I23rF7ty4ZIHaIJ2/6KsA0AAFDJtiQZ/bhGCgiQ/nw3QdufEbYBAAAqkTFGr79pzUBy1ZVSZCvCtj8jbAMAAFSinzdIvyRK1apJE+4kaPs7wjYAAEAlKd6rff21UuNGhG1/R9gGAACoJN//ICVvk2rWYLVIWAjbAAAAlSA722jO61av9p9utyk8nLANwjYAAECl+PhTKTVVatRIuvmP7q4NPAVhGwAA4BIdP270wYdWr/bEe22qXp1ebVgI2wAAAJfo9TeNsrOlLp2l4cPcXRt4EsI2AADAJUjeZrT4O2t78oM22Wz0auMcwjYAAMBFKigweuFFa/jIlVdInToStOGMsA0AAHCRvvxK2r5DCguTHriPoI2SCNsAAAAX4fgJo9ffsnq1/3KPTfXqEbZREmEbAADgIsx8zejMGaljR+naMe6uDTwVYRsAAFzQ3r179cQTT+jqq6/W4MGDdf311+uFF15Qenp6uc/xzDPPKDY2VrGxsdq4cWOJ/YWFhXr99dc1ZswYDR06VBMnTtTOnTtLPVd+fr5uu+023XPPPTLGVPjz2OtxPosWLVJsbKyeeuqpUsuXLB6g/OwB2pw4QMOHX64xY8Zo4sSJmjVrlvbs2VPh88I3EbYBAMB5/fzzz7rrrru0dOlShYWFaeDAgQoODtZnn32m8ePHKzU19YLn2LBhgxYuXHjemTref/99vf322woNDVWfPn20ZcsWTZ48WWfOnClx7Pz587V3715NnTq1ymf/yM0t2rA1U2Sb0Ro9erSGDBmiNm3aaO/evfrggw80btw4Pfnkk6XWHf4lyN0VAAAAnis7O1tPPPGEsrOzdffdd+uee+6RJBljNGvWLM2bN0/PPPOMXnnllTLPkZOTo+eee05t2rRRaGioNm/eXOKY/Px8ffDBB2rfvr3eeustBQcH67vvvtM//vEPffnll7rtttscx544cUJvvvmmrr/+ekVHR1f+h76A5SusnvSQ0O56+83HFRJyLuwbY7RmzRr9+9//1pIlS5SamqqZM2cqKIjI5a/o2QYAAGVasWKFTp48qVatWunuu+92lNtsNk2cOFFNmjRRQkJCmcM9JOntt99WSkqKpk2bVmboPHTokDIzMzVy5EgFBwdLkq644gpVr15dO3bscDr2tddeU1BQkP7yl79UwiesmJ83GP2SaG136yKnoC1Z38ugQYP01ltvKSIiQomJiVqwYEGV1xOeg7ANAADKtH37dklSjx49FBDgHBuCgoLUrVs3SdLq1atLff2uXbs0b948jRkzRj169CjzfTIzMyVJtWrVcpQFBAQoNDTUsU+SNm3apG+//VYTJ05UnTp1LuozXaysLKPnnj83PrxevbKPrVevnuOvAPPnz3d11eDBCNsAAKBMZ8+eleQcgouzB97SerYLCwv13HPPqVatWnrwwQfP+z6NGzeWJP3222+OsoyMDKWnp6tRo0aO873wwgvq0KGDrr322op/mEv02n+MjhyV6pYz4w8fPlwBAQFKSUkp17h2+CbCNgAAKFPdunUlSUeOHCl1/6FDh8rc/9lnn2nLli2aNGnSBXuh69evr+joaH3zzTfauHGjMjIy9Morr6iwsFADBw6UJH3++efauXOnpk6dWqKX3dV+jDP6aqG1PXp0+W7IDA0NVdOmTSVZs7nAPzFaHwAAlKlnz5569913tXbtWqWnpzvCtySlpqZq/fr1kqSsrCyn16WmpmrOnDmKiYnR6NGjy/VekydP1sMPP6z77rvPUTZgwAANGjRIp06dckwL2LlzZ8f+nJwcVatW7aLD94Wm/7PeQ47hI7fcJEW2LP/569atq5SUFGVkZFxU/eD9CNsAAKBM/fr1U3R0tLZv364pU6Zo6tSpat26tXbv3q3nnntO+fn5klRi+r0ZM2YoLy9P06ZNK/d79erVS++++66+/fZbnT59Wp07d9ZVV10lSZo9e7Yk6YEHHpAkrV+/Xi+++KL27t2r6tWra9SoUZoyZYqqV69eoc93vl8EUlJStGnTJm3aLJ3KkNq3k+79s01Ll5b//PY5wKt6ekJ4DsI2AAAok81m03PPPadHHnlEycnJTjOS1KtXT3/+8581d+5c1a5d21G+ZMkS/fjjj5owYYIiIyMr9H5t2rRxBGq75ORkLVy4UH/9619Vt25dpaamaurUqWrbtq2effZZ7d27V2+99ZZq1Kihhx9+uELv98QTT5S5b9GiRdq0aZOOH5dCaklPPm5TcHDFQvOpU6ckyen7gX8hbAMAgPNq0qSJ3nvvPa1atUqbN29WTk6OWrdurSuvvFIrV66UJLVu3dpx/IoVKyRJP/30kxITE53OZb+R8sUXX1RoaKiuvvpqjRlT9lrnxhjNmDFD7dq10w033CBJWrBggXJzc/X000+radOmuvzyy5WSkqIFCxbovvvuU40aNSrlcx8+fG570v02RbaqWNA+c+aMDh48KMn5+4F/IWwDAIALCgoK0vDhwzV8+HCncvsCNTExMSVes2XLljLPZ587u7TXFbdw4UIlJydrzpw5CgwMlCTt27dPdevWddx8KEmdOnXS4sWLdeDAAbVv3758H+o8MjKNvvjSGgLSqJF03UVMfrJs2TIZY9SyZUtFRERccp3gnQjbAADgopw4cULLly9XnTp1dNlllznKn3vuOT366KOlvmbixIlKTEzUnDlzzjvvtmTNvf2f//xHo0aNUvfu3Z325eTkOD3Pzs6WpEqZpcQYo2efM0q3RoCoa5eKj7k+efKk3njjDUnSTTfddMl1gvdi6j8AAHBeu3fvLhFuU1NT9T//8z/KysrS5MmTK23oRnFz585Vbm5uiTHcbdq0UVZWlmMhnfz8fC1fvlzBwcFq1qzZJb/vJ/OlH9dIgUVdktWqlf+1xhitXbtWd999t44fP67evXvr+uuvv+Q6wXvRsw0AAM5r3rx5WrVqlaKjo9WgQQOdPHlSmzZtUm5uriZMmKCrr7660t9z586d+uKLLzRp0iTVr1/fad+NN96oTz75RNOnT1e/fv2UkpKivXv3avz48Zcc+n/dZPSfudbwkZHDbVr8TdnHbtq0SU899ZQkK/CfOnVK27dvV3p6uiRp1KhRmjp1aplL1MM/8F8fAACc19ChQ3Xy5Ent3LlTmzZtUq1atRQbG6ubb75ZvXr1csl7/vvf/1ZkZKRuvPHGEvvq16+vl19+WTNnzlR8fLzCwsJ022236d57772k9zyaajT9SaOCAmn4MCmmp84btlNSUpSSkiJJql69umrVqqXWrVurc+fOGj16tNq0aXNJ9YFvsBn7BJAXkJaW5uq6eI3w8HC+D7gM7QuuRhuDq3ljG8vJMbp/ktH2HVK7ttJ/ZtlUsyZzY3sqd7ex8PDwch/LmG0AAODXjDF6/gUraNepLT37NEEblYewDQAA/NpHn0hLvpcCA6Sn/mFTkyYEbVQewjYAAPBby1cazZ5jjah98AGbesUQtFG5CNsAAMAvbd5i9PQzVtC+caz1A1Q2wjYAAPA7KSlGj/3dKDdPGjRQmvSArcIL1wDlQdgGAAB+JT3daOqjRqcypA7R0pPTbQoMJGjDNQjbAADAb+TkGP1tulHKQalxI+n//h8zj8C1CNsAAMAvFBYaPf2s0eYtUliY9MLzNtWvT9CGaxG2AQCAX5jzutGKlVJQkPT//mVTZCuCNlyPsA0AAHzehx8bffixtf23aTbF9CRoo2oQtgEAgE9b9M25ubTvu9emK68gaKPqELYBAIDPWrnK6Pl/W0F73K3S7eMI2qhahG0AAOCT1v9s9M+njQoLpWuulibeS9BG1SNsAwAAn7Mlyejv043y8qTLL5Om/pVFa+AehG0AAOBT9uwx+p/HjM5mS337SI//nUVr4D6EbQAA4DMOHjSa8j9GmZlSl87SM0/ZFBxM0Ib7ELYBAIBPOHzYaPJfjU6ckNq2kZ5/jtUh4X6EbQAA4PWOplpB++hRqUUL6cUZNtWuRdCG+xG2AQCAVzt2zGjyw0aHD0vNmkqvvsgy7PAchG0AAOC1TpywerQPHpKaNJFefdmmiAiCNjwHYRsAAHiltDSjhx4xOnBAatTI6tFu1JCgDc9C2AYAAF4nPd3o4UeM9u2TIhpYQbtJE4I2PA9hGwAAeJVTp4ymTDXavUeqX1965SWbmjUjaMMzEbYBAIDXSEu3xmjv3CWFh0uvvGhTyxYEbXguwjYAAPAKaWlGD00x2r1bql9PmvmyTZGtCNrwbEHurgAAAMCFnDhh3Qy5b581dGTmSza1bEnQhucjbAMAAI92/LjR5ClGvx0ouhnyZZtaNCdowzsQtgEAgMc6dsxo0hSjlBSpYUOrR5ubIeFNGLMNAAA80tFUo0kPW0G7cSNp1isEbXgfwjYAAPA4B1KM7p9klHLQWhly5ss2NWUebXghwjYAAPAou/cYPTDJ6OhRqUULq0ebBWvgrQjbAADAYyRvs4aOnEyT2rWVXnuFJdjh3QjbAADAIyRutGYdyciQOneyZh2pV4+gDe9G2AYAAG63Lt7okWlGZ89KvWKkl16wqXYtgja8H2EbAAC41fKVRo/9r1FurjRwgPT8szaFhBC04RuYZxsAALjNF18ZvfiykTHSiOHS9L/ZFBRE0IbvIGwDAIAqZ4zR2+8Y/fdd6/n110pTHrIpMJCgDd9C2AYAAFWqoMDqzf5qofV8wp023XWHZLMRtOF7CNsAAKDK5OQYPfWM0arVks0mPfKwTddfR8iG7yJsAwCAKnH6tHUj5MZfpWrVpCen23TZUII2fBthGwAAuNzx40ZTHzXatVsKCZGee8ammJ4Ebfg+wjYAAHCpnbuMHv2bUeoxqV649MLzNkW1J2jDPzDPNgAAcJm164zun2QF7VYtpf+8RtCGf6FnGwAAuMRnnxu9OsuosNBaFfLpf9pUi1Uh4WcI2wAAoFLl5xu99EqhFnxhPR8zWpr6VxargX8ibAMAgEqTlWX098cztfpH6/n999l0683MoQ3/RdgGAACV4uBBo78/brR7T56qV5ee+F+bhg4hZMO/EbYBAMAlW5dg9M9/GZ0+LTVoYNOzT0sdOxC0AcI2AAC4aIWFRu/Pk95828gYqUtnaebLdVWt2il3Vw3wCIRtAABwUc6cMXr6WaMf46zn118rPTTJpoYNA5SW5t66AZ6CsA0AACps/35rfPb+36yl1x952KYxVzNsBPg9wjYAAKiQ5SuNnnveKCtLimggPfMvmzp1JGgDpSFsAwCAcsnONnr1NaOvF1rPe3SXnnrSpnr1CNpAWQjbAADggvbsNXryn0Z790k2m3T7bdLdd7JQDXAhhG0AAFAmY4wWfSO9PNMoJ0eqFy49/r829elNyAbKg7ANAABKdfq00Yx/G/2wwnret480/W8MGwEqgrANAABK+CXR6Nn/Mzp8RAoMkO75s03jbpECAgjaQEUQtgEAgMPZs0ZzXjda8IX1vElj6cnHberSmZANXAzCNgAAkCRt2mz0/54zSjloPb/uGumBiTaFhBC0gYtF2AYAwM/l5Bi98ZbRJ/MlY6SGEdJj02zq24eQDVwqwjYAAH7s101Gz79grQQpSaOvkiY9YFOtWgRtoDIQtgEA8ENpaUaz5xh9u8R6Xr+eNG2qTQMHELKBykTYBgDAjxQWGn29SJr7hlFmprVAzbVjpL/cY1Pt2gRtoLIRtgEA8BPbdxi98JJRcrL1vH076ZEpzDQCuBJhGwAAH5eWbvT2O0ZffS0VFkohIdI9d9t0w3ViuXXAxQjbAAD4qNxco88+l9593+jMGats+DBp0v02NWhAyAaqAmEbAAAfY4zR8pXSnNeNDh+2ytq3s2YZielJyAaqEmEbAAAfsiXJaNZsoy1J1vMGDaR777bpyiukwECCNlDVCNsAAPiAPXuN3nzbaPWP1vMaNaRxt9h0681SzZqEbMBdCNsAAHixg4esmx+Xfm+t/hgQIF11pdWbzbhswP0I2wAAeKFjx4zefd9o4TdSQYFVdtkQ6e4JNrWOJGQDnoKwDQCAFzlxwujDT4y++FLKzbXK+vaxerI7dCBkA56GsA0AgBc4ftxo3kdGXy08F7K7dZXu/bNNPboTsgFPRdgGAMCDpaZaIXvhIik3zyrr3Em66w6b+vWVbDaCNuDJCNsAAHigQ4eNPvrEaNE3Ul5RyO7aRZpwp029exGyAW9B2AYA+Kxt27bpp59+UlJSkrZu3apjx45JkuLj491cs7Lt2Gn1ZK9YaS2tLkk9uls92TE9CdnwP97477g4wjYAwGe9/fbbWr16tburcUHGGG34RZr3kdH6n8+V9+0j3T6OVR/h37zl33FZCNsAAJ/VpUsXtWvXTh07dlSnTp10ww03KNd+d6EHyM83WrVamvex0Y4dVllggDRsmDTuZpvatydkA57+7/hCCNsAAJ81fvx4d1ehVDk5Rt98K338qdGhQ1ZZ9erSNVdLN//RpiZNCNmAnaf+Oy4vwjYAAFUkNdXoi6+Mvl4oncqwyurUlv4w1qax10t16xKyAV9D2AYAwIWMMdq8Rfrsc6NVq6SCopsemzSWbrnJpqtHSzVqELIBX0XYBgDABXJyjJavkOZ/fm48tiT17CHd+AebBg2QAgMJ2YCvI2wDAFCJjh83+vJroy+/ltLTrbLgYOnKkdZwkXZtCdiAPyFsAwBwiQoLjX7eIC38xmj1j1JBgVXeMEK64Xqbrrma8diAvyJsAwBwkY6mGi3+VvpmsdGRo+fKu3W1hooMGSQFBRGyAX9G2AYAoAJycozWrpMWLTb6ab1kjFUeFmYNFRkzmvmxAZxD2AYA4ALy861hIsuWW8NEsrLO7YvpKY252qahg6Xq1QnZAJwRtgEAKEVhoTVl3/c/GK1cKaWfOrevUSPpihFWL3azZgRsAGUjbAMAUOTsWaP1P0tr1llDRdLSzu2rW1cadpk0YrhNXTpLAQGEbAAXRtgGAPisNWvW6O2333Y8z8vLkyTdfffdjrKxY+9SXv5ArVlntGGDlJt37vVhodKQwVbAjunJzY6AO5Tn3/GECRM0cODAKq9beRC2AQA+Ky0tTUlJSSXKi5clb09TQJBxPG/SRBo0QBo4wKbu3aRq1QjYgDuV599xWvE/Q3kYmzHGXPgwz/4QVS08PJzvAy5D+4Kr+VsbO37CaNs2adNmo59/kXbuPDeDiCQFBEidOlrheuAAqXWkZLMRsC+Fv7UxVD13t7Hw8PByH0vPNgDAZ2RlGW3bLiVvk5KTjbZuk1JTSx4XGSn1jpF6xdjUo7tUqxbhGoBrELYBAF4pP99o9x4pOVnaus0oOVnat9+511qyeq5bR0odO0gxPW2KiZEa1CdcA6gahG0AgMczxujgQWlrUY918jZpx04pN7fksY0aWcG6U0ebOnWUotpLISGEawDuQdgGAHictDSjrclS8jb7o5SZWfK4sDBrvHXHDlLHjjZ1jJbq02sNwIMQtgEAblVQYA0H+XWTtHmL0dat0pGjJY+rVk1q307q1Enq1MGmjh2l5s24mRGAZyNsAwCqVE6OUdJWadNma4aQLUnOy59Lks0mtWp5rse6UwepbVum4QPgfQjbAACXOnXKWvZ802ajXzdJ23dI+fnOx4SGSl06S926WqszRkdJYWEEawDej7ANAKg0xhgdOWL1Wv+62WjTZmnfvpLH1a8vde8mde9qU7euUps2UmAg4RqA7yFsAwAuWmGh0Z69ReF6k9HmzVLqsZLHtWopdetm9Vx36yo1bcJYawD+gbANACi3nBxr0Rj7eOvNm6XTZ5yPCQy0hoF062qF665dpfC6BGsA/omwDQAoU2amdQPjr5usISHbtkm5ec7H1Kx5brx1927WTY01axKuAUAibAMAihQWGu3bL21NlrYmWzOG7NlTckXGeuFFvdbdrCEh7dpKQUGEawAoDWEbAPxUWppRUrKUtNWa23rbdunMmZLHNW9u3cxoH2/N3NYAUH6EbQDwccYYHU2Vdu+Wdu+RfvstU79uKtThIyWPrVFD6hBtLRzTuaNNXbtI9eoRrAHgYhG2AcCHZGVZs4Ps3iPt3m2KHn9/E2OupKKFY1pZy5136mhT505S60iGhABAZSJsA4CXyc83OnpUOpBS9HPAFD2Wvsy5ZM0Q0qqltQpj504hatXyrDp28L+FY44fP+7uKviF/Px8paenu7sa8GHh4eHurkK5EbYBwMNkZVnDPo4eVdGj8/Njx0quwFhcgwZSm9bWjYtt29jUtq0VtO1LnYeH11RaWnYVfRrPEhUV5e4qAKgE5vd3bnswwjYAVIHCQqOMDOnkSenESSktzf5odOKkVX7ypLUgTGbmhc8XXM26cbFFC6lFc6lFC5taNJdatpDqMqc1AHgMwjYAVEBurtGZLCnrjHQmSzp9WsrIkDIypVOnpIwMo4xMKeNUUVmGtX3qlFRQWP73CQuTGjWSGjW0P9oczxs2lBpGSAEBhGoA8HSEbQA+JT/fKDtbysmRsnOsx5wcOcoc27nnyrKzzbnj7fvOSllnranwiofrvLwL1+F86tSW6tX73U+4rehRioiwAnVoKEEaAHwBYRuASxljlJtbLOgWC8DFfxzlRWE3J8c4QrFz+bmwnGMP0EXBOTenYr3Hl6JmTSk0VAoNkWrXtn7q1JZq1Zbq1Lapdi2pdh2pdi2rvE4dKTz83LhpuMeOHTvcXQW/ULduXW6QBIoQtgEfYIxRYaFUUHDuJz/f+bmjvJQy+/G5uVZwzS31x5zbn1P6cfbQnFssQOfmllyBsCrYbFKN6lL16lL1GsW2q1tzSTttBxcvsznKQorCtD1Uh4RKYaFW0A4MJDR7owYNGri7Cn4hPDxcQUFEDEAibAMlGGOUl3euxzQnR049s3l51k9urpSbJ+XZH4vKrEfjOC6nqMx+nD3c2sPx70OyzZau3NxCFRaeC8aOY+3HFXutPVR7g8DAooBbLPgW/6lRXQouFoDtYTg42OYUjO0B2ikkF20HF20HB7PKIQDA/QjbuCjGGBmjSv8pNJKMFSDzi3pb8/OKbeef21dQtJ1XrLygaF9u0TAExzCFUkKzo4c2p9i+ov3unVGo8pKzzSYFBVkh9/c/QaWUBQef+7EH1urVncuDg22lHmc/trQAbd9msRQAgL/xirC9dp3R/AVGBQXOIai0QFTa/rKCU2n7HWVOBzrvCwxKV35e4e93O72o1HOf55wXem/7/sKiUCp7QC0qN4Xntp32FTofp2LncBxbyjl+f6zTcd4zteUls9nO9ZoWD53VqjlvV6tmTcVWzR48qxUrLwqn9uf2kBsQcC7kBtgDb4BUp06YsrJOlzgm8PevC7C27WG6tFDNbBUAALiXV4Tt75Yarf/Z3bUozkv+Zu8FAgIkmyTZnLeDgop+AqWgasW2i8oDf/e8mr2s6Mfe42oPyNWr25xDc/VzIbpEj2yx11WrVvVDEcLDg5WWRkgGAMAXeEXYnvpXmy4fao1PlVSUyIo2i2+XVmb73U7nTcd+pzxlO/9xYWFhOn36dKnnLO29y8pqpe0v7djiZQEB58pK/SmqjM0mBdjObdv32QKKbdvKONZ+/oBi2/b99jqUcVzx+p1vm7G0AADAH3hF2K5dy6bLL3N3Lc6h5xEAAADlEeDuCgAAAAC+irANAAAAuAhhGwAAAHARwjYAAADgIoRtAAAAwEUI2wAAAICLELYBAAAAFyFsAwAAAC5C2AYAAABchLANAAAAuAhhGwAAAHARwjYAAADgIoRtAAAAwEUI2wAAAICLELYBAAAAFyFsAwAAAC5C2AYAAABcxGaMMe6uhDfJzMzUhg0b1KtXL9WqVcvd1YGPoX3B1WhjcDXaGFzN29oYPdsVdPr0aa1atUqnT592d1Xgg2hfcDXaGFyNNgZX87Y2RtgGAAAAXISwDQAAALgIYbuCwsLCNHToUIWFhbm7KvBBtC+4Gm0MrkYbg6t5WxvjBkkAAADARejZBgAAAFyEsA0AAAC4CGEbAAAAcBHCNgAAAOAihG0AAADARYLcXYGLtWnTJs2cOVOJiYnKz89XVFSU7rzzTo0ePbrc59i9e7dmz56t+Ph4nTp1ShERERo+fLgefPBB1a1bt8TxOTk5evPNN/XNN9/owIEDqlGjhrp3766JEyeqV69epb5HUlKS5s6dqy1btuj48eOqX7++OnTooAkTJqhPnz4lji8sLNS8efP06aefav/+/QoJCdGAAQM0ZcoUtWjRotyfDZfOF9tYQkKCxo8fX2Z9n332WY0dO7bcnw8Xz1va1+rVq/X2229ry5Ytys3NVatWrTR27FiNHz9egYGBJY7nGuY5fLGNcQ3zHF999ZU2bNigLVu2aMeOHcrLy7uo7/9irhk//vij5s6dq6SkJNlsNnXu3Fn333+/+vfvX+rxe/fu1csvv6z4+HidPXtWkZGRuuWWW3TrrbfKZrOVOP706dOaOXOmli5dqmPHjqlhw4a68sor9eCDDyo0NLRCn0/y0qn/4uPj9ec//1nBwcG6+uqrFRoaqqVLl+rgwYN69NFHNWHChAueY+PGjbrrrruUnZ2t4cOHq0WLFtq2bZvWrl2ryMhIffzxxwoPD3ccn5OTozvuuEOJiYmKjo5WbGysMjMztWTJEmVnZ+vVV1/ViBEjnN5j2bJlmjRpkoKDgzVixAg1adJEhw8f1rJly5SdnV1qo5w+fbrmz5+v9u3ba+jQoUpNTdW3336r0NBQffLJJ4qMjKyU7xDn56ttzP4/qr59+6pv374l6jxixAh17NjxEr45lIe3tK/33ntPzzzzjMLCwnTFFVeoVq1aWrdunXbs2KErr7xSr776aol6cQ3zDL7axriGeY5hw4bp4MGDCg8PV0hIiA4ePHhRYbui14yvvvpK06ZNU7169Ry/OC5evFhpaWl6+eWXddVVVzkdv2vXLt1yyy3Kzs7WqFGj1LBhQ61atUo7d+7U7bffrscff9zp+KysLI0bN07JyckaNGiQOnbsqOTkZMXFxalr166aN2+eqlevXrEvy3iZvLw8M2LECNOlSxezdetWR3lGRoa54oorTOfOnU1KSsoFzzNmzBgTFRVlli1b5lT+xhtvmKioKPP44487lb/55psmKirKTJ482eTn5zvK9+/fb2JiYkxsbKzJzMx0es2oUaNMdHS0Uz2NMWbLli0mOjraDBs2zKl83bp1Jioqytx2220mJyfHUb5y5UoTFRVlJkyYcMHPhUvny20sPj7eREVFmVdfffWC9YdreEv7OnLkiOnSpYvp06ePU33y8vLMxIkTTVRUlFm0aJHTe3AN8wy+3Ma4hnmONWvWOP67zZ0710RFRZkFCxZU6BwVvWakp6eb3r17m379+pnDhw87yg8fPmz69etn+vXrV+L/k7fddpuJiooyK1eudJTl5OSYcePGmaioKPPLL784Hf/KK6+YqKgoM2PGDKfyGTNmmKioKDNnzpwKfUZjjPG6Mdvx8fH67bffNGbMGKffXmvVqqX77rtPeXl5+uKLL857jt9++007duxQ165dNXz4cKd9EyZMUN26dfX1118rKyvLUf7DDz9IkiZNmuT0Z62WLVvqD3/4g06ePKklS5Y4nevAgQOKiIgo8Vt2586dFRERoZMnTzqVz58/X5L00EMPKTg42FE+dOhQ9e3bV3FxcTp06NB5PxsunS+3Mbift7Sv1atXKzc3VzfeeKOaNWvmKA8KCtLkyZMlSR999JHTe3MN8wy+3MbgOQYMGOD03+1iVPSa8d133ykjI0O33367Gjdu7Chv3Lixbr/9dqWlpWnZsmWO8r1792r9+vXq16+fhg4d6igPDg7WQw89JEn69NNPHeXGGM2fP18hISG6//77nep6//33KyQkxFHnivC6sP3TTz9JkgYNGlRin71s/fr15z3HsWPHJEnNmzcvsS8gIEBNmzbV2bNn9euvvzrKjx8/XuZr7GXx8fFO5e3bt9exY8eUnJzsVJ6UlKRjx46VGFuUkJCgkJAQxcTElHiPwYMHSzr3+eE6vtzG7Pbt26d33nlHc+fO1ZdffqmjR4+e9/Og8nhL+yrP8YmJicrNzXWUcw3zDL7cxuy4hvmGil4zytO2y3t8r169FBIS4vRvYd++fUpNTVVMTIxCQkKcjrfX88CBAzp8+HC5P6PkhTdI7tu3T5LUqlWrEvsiIiIUEhKi/fv3n/cc9jFmKSkpJfYVFhY6fovau3evI6yEh4dr//79SklJUbt27ZxeYz+PvW52f//733Xvvffqlltu0ciRI9WkSRMdOnRIy5YtU9++ffXPf/7TcWxWVpaOHTumqKioUm86sn/eC302XDpfbWPFLVq0SIsWLXI8DwoK0u23365p06aV2v5QebylfZ3vPexl+fn5OnDggNq2bcs1zIP4ahsrjmuY97uYa8b52nZFjw8MDFTz5s21a9cu5efnKygoyPHasu4tiYyMVFxcnPbt26cmTZpc+EMW8bqe7dOnT0uy/hxWmrCwMGVmZp73HK1bt1aLFi20efNmrVy50mnfu+++q/T0dElyOo/9N6zXXntNBQUFjvIDBw7o888/lyRlZGQ4nat379768MMP1bhxYy1cuFCvv/66Fi1apPr162vs2LGKiIhwHGt/r7CwsDI/1+/rBNfw1TYmSfXq1dMjjzyiRYsWKTExUWvXrtVrr72mli1b6p133tGMGTPO+7lw6bylfQ0cOFCBgYFasGCBUy9Ofn6+Zs2a5Xhufw3XMM/hq21M4hrmSy7mmnG+tl3R4yUpNDRUhYWFOnPmTIXqZD9veXldz3ZlsNlsevLJJzVx4kRNnDhRI0aMUIsWLbR9+3bFxcUpKipKO3bscJoO5s4779S3336rxYsXa8+ePYqNjVVGRoaWLl2qZs2aafv27QoIcP7dZeXKlXrkkUd02WWXafbs2WrWrJkOHjyo2bNn69FHH9WOHTs0bdq0qv74qAKe2sbat2+v9u3bO56HhIRoxIgR6t69u6699lq9//77uueee1S/fn3Xf0m4aFXRvlq0aKG//OUvmj17tq655hqNHDnSMVPE4cOH1bRpUx06dKhEm4Rv8NQ2xjUM3sjrrpIX6h05ffp0mb/BFDd48GDNmzdPQ4YMUXx8vN5//32lpaXptddec0wnVPwfa1hYmD766CPdeeedyszM1Lx587RmzRrdcssteuKJJ0ocn5aWpqlTp6pVq1aaMWOG2rZtqxo1aqht27aaMWOGOnfurHfeecfxpzh7ncv6belCv52h8vhqGzsf+9y5+fn5TmMwUfm8pX1J1k1Lzz//vFq3bq3Fixfrs88+U+PGjfXRRx85Pke9evUkcQ3zJL7axs6Ha5j3uZhrxvnadkWPl6QzZ87IZrM55s4ub53K6vkui9f1bNvH0ezfv19dunRx2nfs2DFlZWWpW7du5TpX9+7dNXfu3BLl7777riSVOH/t2rX1t7/9TX/729+cyu1/Hit+fGJiojIzM9W3b98SPT8BAQHq06ePkpKStH37djVt2lQhISGKiIhQSkqKCgoKSoxfso8jKm3cESqXr7axC7GPnzx79mw5Phkulre0L7vrrrtO1113nVNZbm6u9u3bp7p16zoWneAa5jl8tY1dCNcw73Ix14zIyEht2bJF+/fvd5rj/XzHF99XXEFBgVJSUtS8eXMFBQU5vfb390fZ2csrul6A1/Vs21fEi4uLK7HPXlbayozldfDgQW3YsEHt2rVTdHR0uV6zcOFCSXJalSsvL0+Sypx6zV5efKqbvn37KisrS7/88kuJ43/88UdJl/bZUD6+3MbOx94bdKlTOeH8vKV9nc93332n3NzcEsdzDfMMvtzGzodrmPep6DWjPG27+GJH5zt+w4YNysrKcjp/ZGSkGjZsqF9++cVpWktJjno2b968QjdHSl4Ytvv3768WLVpo0aJFTtOdZWZmas6cOapWrZquv/56R3lqaqp2795d4k8IZ86ckfnd4pmZmZmaNm2aCgoK9Ne//rXEe5f2Z4V33nlHa9eu1ciRI516Crp3767AwEAtWbJE27Ztc3pNcnKylixZopo1a6p79+6O8ptuukmS9MorrzhNdbRq1Sr99NNPGjRoEBeRKuDLbWzLli2lfuZ3331XCQkJioyMVNeuXUs9BpXDW9pXWcfv3r1b//d//6ewsDDde++9Tvu4hnkGX25jXMO808mTJ7V79+4SnUMVvWaMGjVKtWrV0gcffKAjR444yo8cOaIPPvhA4eHhTquUtmnTRn369FFCQoJWrVrlKM/NzdUrr7wiSfrjH//oKLfZbPrjH/+orKwszZ4926mus2fPVlZWlqPOFeHzy7U/9thj+uKLL0osIfrVV1/ppZdeUmxsrBo2bKgTJ05o+fLlOnnypB566KESk5lLUs+ePdWvXz9FRkbKZrMpISFBSUlJ6tKli95++23VqVPH6fiXXnrJcWEbOXKkmjZtqoMHD2rZsmXKy8vT9OnT9ac//cnpNb9ftvTYsWNavHixQkND9fHHH6t169aV/G2iNL7axoYNG6agoCB16dJFjRo1csyTu3XrVtWuXVtvvfVWuf+8jIvnLe1r+vTp2rp1q7p06aK6detq3759WrFihQICAvTaa6+VOnct1zDP4KttjGuY55g/f742bNggSdqxY4eSkpIUExPjGIrRq1cvR5CdOXOmZs2apQcffFCTJk1yOk9FrxnnW679pZde0qhRo5yO37lzp2699VZlZ2dr9OjRioiIuOBy7bfeequ2bdumQYMGqVOnTtq6datjufYPPvhANWrUqNB35XVjtiUpNjZWH374oV599VUtXrxY+fn5ioqK0tSpU8v9J6fo6Gh16NBBcXFxSk9PV1hYmHr06KE777xTsbGxpb7m2muvVUJCguLj42Wz2RQZGalp06bpT3/6U6l/qp8yZYqio6P18ccfKy4uTmfOnFGtWrXUr18/jR8/3mk1I7unnnpKUVFR+vTTT/Xee+8pJCREI0eO1JQpU9SyZcuKfVG4aL7axm655RbFxcVp/fr1Sk9PdyxOcccdd2jChAlOK3LBdbylfQ0cOFB79uzRkiVLdObMGTVo0EDXXHON7rvvvjKvR1zDPIOvtjGuYZ5jw4YNJVYi/eWXX5yGhBTvNS5LRa8Z1113ncLDwzV37lynewEmTpyoAQMGlDi+ffv2+vTTT/Xyyy9r1apVysrKUmRkpJ544gmNGzeuxPEhISH64IMPNHPmTC1dulQJCQmKiIjQhAkT9MADD1Q4aEte2rMNAAAAeAOvG7MNAAAAeAvCNgAAAOAihG0AAADARQjbAAAAgIsQtgEAAAAXIWwDAAAALkLYBgAAAFyEsA0AAAC4CGEbAAAAcBHCNgAAAOAihG0AAADARQjbAAAAgIv8f7e+iUwkFVrkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# why don't i recover the .5??\n",
    "\n",
    "az.plot_posterior(idata, var_names=['f']);"
   ]
  },
  {
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}