ab-test.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# AB Test\n",
    "\n",
    "\n",
    "PyMC and Bernoulli distribution\n",
    "\n",
    "[Pythonで体験するベイズ推論](https://www.amazon.co.jp/Python%E3%81%A7%E4%BD%93%E9%A8%93%E3%81%99%E3%82%8B%E3%83%99%E3%82%A4%E3%82%BA%E6%8E%A8%E8%AB%96-PyMC%E3%81%AB%E3%82%88%E3%82%8BMCMC%E5%85%A5%E9%96%80-%E3%82%AD%E3%83%A3%E3%83%A1%E3%83%AD%E3%83%B3-%E3%83%87%E3%83%93%E3%83%83%E3%83%89%E3%82%BD%E3%83%B3-%E3%83%94%E3%83%AD%E3%83%B3/dp/4627077912)より"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymc as pm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = [12, 8]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# conversion rate\n",
    "p = pm.Uniform('p', lower=0, upper=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "74\n"
     ]
    }
   ],
   "source": [
    "p_true = 0.05\n",
    "N = 1500\n",
    "\n",
    "occurences = pm.rbernoulli(p_true, size=N)\n",
    "print(occurences.sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " [-----------------100%-----------------] 20000 of 20000 complete in 4.1 sec"
     ]
    }
   ],
   "source": [
    "obs = pm.Bernoulli('obs', p, value=occurences, observed=True)\n",
    "mcmc = pm.MCMC([obs, p])\n",
    "mcmc.sample(20000, 10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEcCAYAAADQqlM0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmYHWWZ/vHvbXbCEkgCAglJ0CiBCNnYJwEBEQFNQBlQkCAqMjI6oowsM2NgfjhmBEEZ3KKsw74GRMeFDEgACXQCEUJgCDEhzZZO2LJvPL8/qhpOOqf7nOrT51Qv9+e6+uqut6reeuo9VefpemtTRGBmZlauD+QdgJmZdSxOHGZmlokTh5mZZeLEYWZmmThxmJlZJk4cZmaWiROHmZll4sRhZmaZOHG0gqR5kg6tYv3XSrq4GssqrE/SIklHVKPuWpH0UUlPSloh6Zu1XHaTONq0LduTpm1czTbPYxuqtda2X3tqm3afONIdco2klZJel3SNpK0rrK+iHTwi9oqIByupo62XVe56tVXsxZZXy3Yp8F3gwYjYJiKuqNVCq5ko2mESatrGbdLm7WgbqrVWtV9h27TVNiLpQUlvSuqVZb52nzhSn46IrYExwL7Av+YRhKTuec7fUZddZUOAeXkHUUs5fJZN27jLtXkbaxftJ2koMB4I4DOZZo6Idv0DLAKOKBi+BLgv/XsE8CDwFskH8ZmC6c4FXgZWAM8DhwP/DbwLrAFWAt9Np90FuBNoAP4GfLPJ8s8F/gqsA7oXiamlOLaYv8g6jgbmpLHeCtwCXNzM+mdZrxZjT/8+H3gWeBO4BuhdsKwAPlwwfC1wcYnlHVGqTQqmPSeN7e10vXs3bZuW6gL+F9gErE3j+EiReU8GHk3rfw1YAnyqwm1yi/UvtT60sI2VqruFz7Lo55NleaXiK9LGW7R5qWUBg4G70vHLgStLbUPAecAdTer5CXBFOW1ayedBkX2spfJyt9cM22xzy29sm8zfY83E+T3gEeAy0u/UsveBSnagWvyw+ZfR4PSD+H9AD2ABcAHQEzgsbeiPpj9LgF3S+YYCH2paXzr8AWB22og9gd2BhcAnC6Z/Kl12nyIxNRtHc/M3Wb+ewGLg7LSuzwEbKJI4sqxXmbEvAp5Jx++QbkQXF8zf0hdTc8s7olSbFEz7OMnGvgMwHzizSPuUat8Hga+0sP1MJdnB/j6t6xxgcVtul6XWhxLbWKm6W/gsm0vsWZdXah/YrI0Lh8uYtxswF7gc6Av0Bv6ujG1oCLAa2LagnleBAzLst5k/D5rZx5orz7q9ltpmW1oOW+67ZX+PNbOsBcDXgbEk3zk7lbv9d5SuqumS3gIeBv4M/AdwALA1MDUi1kfE/wL3AZ8nyei9gD0l9YiIRRHxYjN17wsMjIh/T+tZCPwKOKlgmisiYklErCkyf0txlDt/D+DHEbEhIu4Anmgm1izrVc6yAa5Mx78BfL9J3K1VTps0xvZKuuzfAKMqqKs5HwMuj4jbImIDcD2wm6TejRNI2kvSJkmDsqxkEc2tTznbWLn1t/RZNsq6vEriKzXvfiRf3v8cEasiYm1EPFyq0ohYTHIUPiktOgxYHRGPZYi5NZ9Hc/tYuftepdtra/bxctvjPZL+jiQ53xYRs4EXgS+UGSMdpd97UkTcX1ggaRdgSUS8W1C8GNg1IhZI+hZwIbCXpD8A346IV4rUPQTYJU1MjboBMwuGl7QQW7NxZJj/5Uj/BSiYfwsZ16ucZTcdvziNp1LltAkkXUeNVjez7HLras7HgH8rGN4RWBkRawvKziU5/B8B1JdZbzHNrU8521g5Sn2WjbIur5L4Ss07mOQIb2NZkW/uJpIv3OtJvtRuyhhz5s+jhX2s3H2vou21lft4i+vUzPSTgT9GxLJ0+Ka07PJy4uwoRxzFvAIMllS4DruR9A0SETdFRGNWDeA/02mCzS0B/hYR/Qp+tomIowumaTpP2XGUMf+rwK6S1GT+ojKsVznLhmTHLlxu4Qa6GtiqYPiDZdZbTpuUq9V1SepHsn4NBcWfA/6nYJq9ST6DP5AkjnKVatdC5Wxj5dTdtLy5zyfr8rJOn2XeJSRHeMX+SS3VhrcDh6ZHgsexeeKoWszN7WMt7HuFKt72y1xOa77HAJDUh6Tr9hBJr0l6jaSrfB9J+5QTY0dOHLOAVcB3JfVIr2/+NHBLep30YeklZmtJ+rg3pfO9TtL/1+hx4B1J50rqI6mbpJGS9q00jjLn/wuwEfimpO6Sjic5vN9CxvUq11mSBknagaRf9taCcU8BX0jb5CjgkIJxLS2v0jZpq7o+RtI+X0jb9hiSPt0LC6Y5m2THfJaCxKHkXpprW6g7S3tn3cbKrbu5zyfr8irZB0rN+zhJYp4qqa+k3pIOLmc9I6KB5HzANSRfivOrHXNz+1iJfa9QRdt+huVU8j02Ka1zT5Luu1Ek2/5M4NRy4uywiSMi1pNcQvYpYBnwM+DUiHiOpI9walr+Gkn3xAXprD8A/lXSW5LOiYhNJB/sKJIrEZYBvwa2a4M4yp3/eOA0kiubTiS5AqWYsternGWnbgL+SHIibSHJydVG/0TSNm+RXJ00vWBcs8urtE3asK6PATcCB5K07UUk3Z7PAkgaBRxM0gaXkpyYbDSY5GKB5pTd3q3Yxsqtu+jnk3V5lewDpeYtGP9h4CWSrsATM6znTSQnywuPNqoZc3P7WEv7XmHdlW77ZS2Hyr7HJgPXRMRLEfFa4w/J1W4nN3N0uBlt3rVu1nlI+jnwfxFRtN9W0i3AP0TEm+nw4xGxn6SeJFcC7R3JCXUzK9BhjzjMyvAxksswtyBpLLCmMWmk1krqn16VMsJJw6w4H3FYp5VeYTIqIhblHYtZZ+LEYWZmmbiryszMMukoNwAWNWDAgBg6dGjeYZiZdSizZ89eFhEDWzt/h04cQ4cOpa6uLu8wzMw6FElFn05RLndVmZlZJk4cZmaWiROHmZll0qHPcZhZyzZs2EB9fT1r164tPbF1Or1792bQoEH06NGjTet14jDrxOrr69lmm20YOnQomz+A2Tq7iGD58uXU19czbNiwNq3bXVVmndjatWvp37+/k0YXJIn+/ftX5WjTicOsk3PS6Lqq9dk7cZiZWSZOHGZmlokTh1kNPfroozz66KN5h2FWEScOsxo66KCDOOigg/IOo6beeustfvazn+UdRlFr1qzhkEMOYdOmYm9nbdmiRYsYOXJkFaIq3/r165kwYQIbN26s6XKdOMxqqCsecTSXOCKCd999N4eI3nf11Vdz/PHH061bt1zjaK2ePXty+OGHc+utt9Z0uU4cZjV0wQUXcMEFxV4h3Xmdd955vPjii4waNYoTTjiBESNG8PWvf50xY8awZMmSLf5zv/TSS7nwwgsBuOGGG9hvv/0YNWoUX/va14oeGZx00kmceOKJ7L///gwZMoTf/va3Zcd24403MnHiRGDLI4jGOBYtWsSIESP46le/yl577cWRRx7JmjVrNqtn4cKFjB49mieeeKLk9JdddhkjR45k5MiR/PjHPwbghz/8IVdccQUAZ599NocddhgAM2bM4JRTTmmxzkmTJnHjjTeWvc5twYnDrAs59NBDt/hpPBpYvXp10fHXXnstAMuWLdtiXDmmTp3Khz70IZ566ikuueQSnn/+eU499VSefPJJhgwZ0ux88+fP59Zbb+WRRx7hqaeeolu3bkW/IOfOncvuu+/OrFmzuPHGG7nooovKimv9+vUsXLiQcl7N8MILL3DWWWcxb948+vXrx5133vneuOeff57PfvazXHPNNey7774tTj979myuueYaZs2axWOPPcavfvUrnnzySSZMmMDMmTMBqKurY+XKlWzYsIGHH36Y8ePHt1jnyJEjeeKJJ8pa57bixGFmNTVkyBAOOOCAktPNmDGD2bNns++++zJq1ChmzJjBwoULN5tmzZo1LFu2jClTpgCw55578uabb242zYknnsiPfvSjLepftmwZ/fr1KyvmYcOGMWrUKADGjh3LokWLAGhoaGDixInccMMN741vafqHH36Y4447jr59+7L11ltz/PHHM3PmTMaOHcvs2bNZsWIFvXr14sADD6Suro6ZM2e+lziaq7Nbt2707NmTFStWlLUubcGPHDHrQh588MFmx2211VYtjh8wYECL48vVt2/fzYa7d+++2bmOxjudI4LJkyfzgx/8oNm6nnnmGYYPH07v3r0BmDNnDvvss8974++55x6OPfZY7r///i3m7dOnz2Z3VTcXB0CvXr3e+7tbt27vdRNtt912DB48mEceeYS99tqr5PTNvaq7R48eDB06lGuuuYaDDjqIvffemwceeIAXX3yRESNGsHjx4mbrBFi3bt17bVALPuIws6raZpttWvxveKeddmLp0qUsX76cdevWcd999wFw+OGHc8cdd7B06VIA3njjDRYv3vz9Q3PnzuWll15i7dq1rFq1iilTpnD22WcDyRf/7bffzhe/+EXefvvtLZa7/fbbs2nTpvcSRHNxtKRnz55Mnz6d66+/nptuuqnk9BMmTGD69OmsXr2aVatWcffdd793RDFhwgQuvfRSJkyYwPjx4/nFL37BqFGjSt79vXz5cgYOHNjmDzJsiY84zGqo8WRoV9K/f38OPvhgRo4cyYgRI7YY36NHD773ve+x//77M2zYMPbYYw8g6Xa6+OKLOfLII3n33Xfp0aMHP/3pTzc7LzJ37lxOPvlkDj30UN555x0uuOACDj74YAAuueQSVq5cyZlnnsm8efNYs2YNffr02WzZRx55JA8//DBHHHFEs3GU0rdvX+677z4+8YlP0Ldv382OeJoaM2YMp512Gvvttx8AX/nKVxg9ejQA48eP5/vf/z4HHnggffv2pXfv3u8llZY88MADHH300WXF2mYioio/wNXAUuCZJuXfAJ4H5gE/LCg/H1iQjvtkOcsYO3ZsmFnznn322bxDqKrx48fHc889t0X54sWL4/TTT39v+MILL4zHHntsi+nmzJkTp5xySlVjrLbjjjuuaBs0KrYNAHVRwfd7NY84rgWuBK5vLJD0cWAisHdErJO0Y1q+J3ASsBewC3C/pI9ERPa7cszasca+9iOOOCLnSDqHF198keHDh29Rvttuu3HVVVe9N9x48ryp0aNH8/GPf5xNmzZ1yHs51q9fz6RJk/joRz9a0+VWLXFExEOShjYp/gdgakSsS6dZmpZPBG5Jy/8maQGwH/CXasVnloeLL74YcOJoKy+//HLFdZx++ultEEk+evbsyamnnlrz5db65PhHgPGSZkn6s6R90/JdgSUF09WnZVuQdIakOkl1DQ0NVQ7XzMyaqnXi6A5sDxwA/DNwm5JLBopdNlD0urWImBYR4yJi3MCBA6sXqZmZFVXrxFEP3JWen3kceBcYkJYPLphuEPBKjWMz65SimXsHrPOr1mdf68QxHTgMQNJHgJ7AMuBe4CRJvSQNA4YDj9c4NrNOp3fv3ixfvtzJowuK9J3j1bgxsGonxyXdDBwKDJBUD0whuUT3aknPAOuByemlYfMk3QY8C2wEzvIVVdYZ/fKXv6zp8gYNGkR9fT0+H9g19e7dm0GDBrV5verI/4mMGzcu6urq8g7DzKxDkTQ7Isa1dn7fOW5WI0PP+y2rF8wCYKsP7190mkVTj6llSGat4sRhVkPvPH430HziMOsI/JBDMzPLxInDzMwyceIwM7NMnDjMzCwTnxw3q6EBx34n7xDMKubEYVZD3bf189Ws43NXlVkNrZr/EKvmP5R3GGYV8RGHWQ2tePJ3APQdMSHnSMxaz0ccZmaWiY84rFMbet5vS07jx3yYZeMjDjMzy8SJw8zMMnFXlVkNDZx0ft4hmFWsakcckq6WtDR9aVPTcedICkkD0mFJukLSAkl/lTSmWnGZ5anbVtvRbavt8g7DrCLV7Kq6FjiqaaGkwcAngJcKij9F8rrY4cAZwM+rGJdZblY+fT8rn74/7zDMKlK1xBERDwFvFBl1OfBdoPDVgxOB6yPxGNBP0s7Vis0sL04c1hnU9OS4pM8AL0fE3CajdgWWFAzXp2XF6jhDUp2kOr9H2cys9mqWOCRtBfwL8L1io4uUFX0ZekRMi4hxETFu4EA/98fMrNZqeVXVh4BhwFxJAIOAOZL2IznCGFww7SDglRrGZmZmZarZEUdEPB0RO0bE0IgYSpIsxkTEa8C9wKnp1VUHAG9HxKu1is3MzMpXtSMOSTcDhwIDJNUDUyLiqmYm/x1wNLAAWA18qVpxmeVpxxMuzDsEs4pVLXFExOdLjB9a8HcAZ1UrFrP24gM9eucdglnF/MgRsxpaMee3rJhT+sGLZu2ZE4dZDa16biarnpuZdxhmFXHiMDOzTPyQQ7MylHqvh9/pYV2JjzjMzCwTJw4zM8vEXVVmNfTBL0ytaH6/CtfaAx9xmJlZJk4cZjX09qy7eHvWXXmHYVYRJw6zGlrz4uOsefHxvMMwq4gTh5mZZeLEYWZmmThxmJlZJr4c17q8ci5xbSvq3qtmyzKrFicOsxra6e8vyjsEs4pVratK0tWSlkp6pqDsEknPSfqrpLsl9SsYd76kBZKel/TJasVlZmaVqeY5jmuBo5qU/QkYGRF7A/8HnA8gaU/gJGCvdJ6fSepWxdjMcvHWIzfz1iM35x2GWUWqljgi4iHgjSZlf4yIjengY8Cg9O+JwC0RsS4i/kbyCtn9qhWbWV7WLp7L2sVz8w7DrCJ5XlV1OvA/6d+7AksKxtWnZWZm1s7kkjgk/QuwEbixsajIZNHMvGdIqpNU19DQUK0QzcysGTVPHJImA8cCJ0dEY3KoBwYXTDYIeKXY/BExLSLGRcS4gQMHVjdYMzPbQk0vx5V0FHAucEhErC4YdS9wk6TLgF2A4YAf6GOdTrc+27Y4vpb3lJi1VtUSh6SbgUOBAZLqgSkkV1H1Av4kCeCxiDgzIuZJug14lqQL66yI2FSt2MzyMvC4C/IOwaxiVUscEfH5IsVXtTD994HvVyseMzNrG75z3Dqsjtit8+afrwVg+0NOyzUOs0o4cZjV0LqXn8s7BLOK+em4ZmaWiROHmZll4sRhZmaZ+ByHWQ1132ZA3iGYVcyJw6yGBnz6nLxDMKuYu6rMzCwTH3GYtYFy7yl54/5pAOxwxBnVDMesqpw4zGpo/dKFeYdgVjF3VZmZWSZOHGZmlokTh5mZZeJzHGY11GMHvxHZOj4nDrMa6n/UN/IOwaxi7qoyM7NMykocku6UdIykshONpKslLZX0TEHZDpL+JOmF9Pf2abkkXSFpgaS/ShqTfVXM2r/lv/8vlv/+v/IOw6wi5XZV/Rz4EnCFpNuBayOi1IsFrgWuBK4vKDsPmBERUyWdlw6fC3yK5D3jw4H90+XtX+5KWOfUEV/UVMqGN17OOwSzipV1BBER90fEycAYYBHJO8MflfQlST2amech4I0mxROB69K/rwMmFZRfH4nHgH6Sds62KmZmVgtZup76A6cBXwGeBH5Ckkj+lGF5O0XEqwDp7x3T8l2BJQXT1adlxeI4Q1KdpLqGhoYMizYzs7ZQ7jmOu4CZwFbApyPiMxFxa0R8A9i6DeJQkbIoNmFETIuIcRExbuDAgW2waDMzy6Lccxy/jojfFRZI6hUR6yJiXIblvS5p54h4Ne2KWpqW1wODC6YbBLySoV6zDqHnjrvnHYJZxcrtqrq4SNlfWrG8e4HJ6d+TgXsKyk9Nr646AHi7sUvLrDPZ4Ygz/GRc6/BaPOKQ9EGScw19JI3m/S6lbUm6rVqa92bgUGCApHpgCjAVuE3Sl4GXgBPSyX8HHA0sAFaTXMFlZmbtUKmuqk+SnBAfBFxWUL4CuKClGSPi882MOrzItAGcVSIWsw5v2W8uBfwmQOvYWkwcEXEdcJ2kz0bEnTWKyazT2rhiWd4hmFWsVFfVKRFxAzBU0rebjo+Iy4rMZmZmnViprqq+6e+2uOTWzMw6gVJdVb9Mf19Um3DMzKy9K/cGwB9K2lZSD0kzJC2TdEq1gzPrbHrtuge9dt0j7zDMKlLufRxHRsQ7wLEkN+t9BPjnqkVl1kltf8hpbH/IaXmHYVaRchNH44MMjwZujoimDy80M7MuotxHjvxG0nPAGuDrkgYCa6sXllnn1HD3fwAw8LgWb4Mya9fKfaz6ecCBwLiI2ACsInkUupllsGnNO2xa807eYZhVJMs7x0eQ3M9ROM/1zU1sZvko5wVYi6YeU4NIrLMqK3FI+m/gQ8BTwKa0OHDiMDPrcso94hgH7Jk+U8rMzLqwchPHM8AHAT/q3KwCvYfsk3cIZhUrN3EMAJ6V9DiwrrEwIj5TlajMOql+Bzf30GizjqPcxHFhNYMwM7OOo9zLcf8MLAJ6pH8/Acxp7UIlnS1pnqRnJN0sqbekYZJmSXpB0q2Sera2frP26vXbpvD6bVPyDsOsIuU+q+qrwB3AL9OiXYHprVmgpF2Bb5LcEzIS6AacBPwncHlEDAfeBL7cmvrN2rPYuI7YuK70hGbtWLmPHDkLOBh4ByAiXgB2rGC53UleR9ud5BW0rwKHkSQngOuASRXUb2ZmVVJu4lgXEesbB9Iv/FZdmhsRLwOXkrxz/FXgbWA28FZEbEwnqyc5qtmCpDMk1Umqa2hoaE0IZmZWgXITx58lXUBylPAJ4HbgN61ZoKTtSR5XMgzYheRlUZ8qMmnRxBQR0yJiXESMGzhwYGtCMDOzCpR7VdV5JOccnga+BvwO+HUrl3kE8LeIaACQdBdwENBPUvf0qGMQ8Eor6zdrt/p8aL+8QzCrWFmJIyLelTQdmN74hV+Bl4ADJG1F8rTdw4E64AHgc8AtwGTgngqXY9bubLf/8XmHYFaxFruqlLhQ0jLgOeB5SQ2SvtfaBUbELJKT4HNIjmA+AEwDzgW+LWkB0B+4qrXLMDOz6il1xPEtkqup9o2IvwFI2h34uaSzI+Ly1iw0IqYATS9mXwj4ON46tdduOg+AD35has6RmLVeqZPjpwKfb0waABGxEDglHWdmZl1MqcTRIyKWNS1Mz3P0KDK9mZl1cqUSx/pWjjMzs06q1DmOfSQVe8+lgN5ViMfMzNq5FhNHRHSrVSBmXUHfPcbnHYJZxbK8c9zMKrTNmPbxru9S7yX3O8mtJeU+csTM2sC7G9by7oa1eYdhVhEfcZjV0NLbLwR8H4d1bD7iMDOzTJw4zMwsEycOMzPLxInDzMwy8clxy0Wpy0E7q60/dkTeIZhVzInDrIacOKwzcFeVWQ1tWv02m1a/nXcYZhVx4jCroYbpP6Bh+g/yDsOsIrkkDkn9JN0h6TlJ8yUdKGkHSX+S9EL6e/s8YjMzs5bldcTxE+D3EbEHsA8wHzgPmBERw4EZ6bCZmbUzNU8ckrYFJpC+Uzwi1kfEW8BE4Lp0suuASbWOzczMSsvjiGN3oAG4RtKTkn4tqS+wU0S8CpD+3rHYzJLOkFQnqa6hoaF2UZuZGZBP4ugOjAF+HhGjgVVk6JaKiGkRMS4ixg0cOLBaMZpVxTajj2ab0UfnHYZZRfK4j6MeqI+IWenwHSSJ43VJO0fEq5J2BpbmEJtZVfUdMSHvEMwqVvMjjoh4DVgi6aNp0eHAs8C9wOS0bDJwT61jM6u2je80sPEdd7Fax5bXnePfAG6U1BNYCHyJJIndJunLwEvACTnFZlY1y+77EeD3cVjHlkviiIingHFFRh1e61jMzCwb3zluZmaZOHGYmVkmThxmZpaJH6tuVkPb7ndc3iGYVcyJw6yGtvrw/nmHYFYxJw5rc1317X7l2LC8HoAe/QflHIlZ6/kch1kNLf/DlSz/w5V5h2FWER9xmNkWyjlqXDT1mBpEYu2RjzjMzCwTJw4zM8vEicPMzDLxOQ6zGtruoJPyDsGsYk4cZjXUZ+iovEMwq5i7qsxqaP3rC1n/+sK8wzCriBOHWQ29MWMab8yYlncYZhXJLXFI6ibpSUn3pcPDJM2S9IKkW9OXPJmZWTuT5xHHPwHzC4b/E7g8IoYDbwJfziUqMzNrUS6JQ9Ig4Bjg1+mwgMOAO9JJrgMm5RGbmZm1LK8jjh8D3wXeTYf7A29FxMZ0uB7YtdiMks6QVCeprqGhofqRmpnZZmp+Oa6kY4GlETFb0qGNxUUmjWLzR8Q0YBrAuHHjik5j1l71mzA57xDMKpbHfRwHA5+RdDTQG9iW5Aikn6Tu6VHHIOCVHGIzq6reg0bkHUKbaavH5/thiR1PzbuqIuL8iBgUEUOBk4D/jYiTgQeAz6WTTQbuqXVsZtW2tn4+a+vnl57QrB1rT/dxnAt8W9ICknMeV+Ucj1mbe+uh63jroevyDsOsIrk+ciQiHgQeTP9eCOyXZzxmZlZaezriMDOzDsCJw8zMMnHiMDOzTPxYdbMa2uHwM/IOwaxiThxmNdRzp93zDsGsYk4clllb3fjVFa1Z9BTgFzpZx+bEYVZDbz96C+DEYR2bT46bmVkmThxmZpaJE4eZmWXixGFmZpn45LhZDfX/5D/mHYJZxZw4zGqoR/9BeYdgVjEnDtuM79GortULZgGw1Yf3zzkSs9Zz4jCroXcevxtw4rCOreYnxyUNlvSApPmS5kn6p7R8B0l/kvRC+nv7WsdmZmal5XFV1UbgOxExAjgAOEvSnsB5wIyIGA7MSIfNzKydyeOd469GxJz07xXAfGBXYCLQ+E7N64BJtY7NzMxKy/U+DklDgdHALGCniHgVkuQC7NjMPGdIqpNU19DQUKtQzcwsldvJcUlbA3cC34qIdySVNV9ETAOmAYwbNy6qF6FZ2xtw7HfyDqHdKXUl36Kpx9QoEitXLolDUg+SpHFjRNyVFr8uaeeIeFXSzsDSPGIzq6bu2w7MOwSziuVxVZWAq4D5EXFZwah7gcnp35OBe2odm1m1rZr/EKvmP5R3GGYVyeOI42Dgi8DTkp5Kyy4ApgK3Sfoy8BJwQg6xmVXViid/B0DfERNyjsSs9WqeOCLiYaC5ExqH1zIWMzPLzneOm1m7Vs5jcHwCvbacOLoQP4fKOisnl9ry+zjMzCwTH3GY1dDASefnHYJZxZw4OhF3RbV/3bbaLu8QzCrmriqzGlr59P2sfPr+vMMwq4iYK99aAAAHoUlEQVQTh1kNOXFYZ+DEYWZmmThxmJlZJk4cZmaWiROHmZll4stxzWpoxxMuzDsEs4o5cZjV0Ad69M47BLOKuavKrIZWzPktK+b4Rk3r2HzEYVZDq56bCcA2Y/zAvc6qK7wK14nDzLqErvCFXivtLnFIOgr4CdAN+HVETM05pKrzM6bMrCNpV4lDUjfgp8AngHrgCUn3RsSz+UZWGScGs/bP7/QoX7tKHMB+wIKIWAgg6RZgItDmiaOtNhInBbOuoy3297b6zsgziSkiclt4U5I+BxwVEV9Jh78I7B8R/1gwzRnAGengR4HnaxjiAGBZDZfXEbmNSnMbleY2Kq2SNhoSEQNbu+D2dsShImWbZbaImAZMq004m5NUFxHj8lh2R+E2Ks1tVJrbqLQ826i93cdRDwwuGB4EvJJTLGZmVkR7SxxPAMMlDZPUEzgJuDfnmMzMrEC76qqKiI2S/hH4A8nluFdHxLycwyqUSxdZB+M2Ks1tVJrbqLTc2qhdnRw3M7P2r711VZmZWTvnxGFmZpl02cQh6ShJz0taIOm8IuN7Sbo1HT9L0tAm43eTtFLSOeXW2dFUqY0WSXpa0lOS6qq/FtXV2jaSNFTSmrQdnpL0i4J5xqZttEDSFZKKXabeYVSpjR5M62wct2Pt1qg6KtnfJO0t6S+S5qXbTu+0vDrbUkR0uR+SE+8vArsDPYG5wJ5Npvk68Iv075OAW5uMvxO4HTin3Do70k812igtWwQMyHv98m4jYCjwTDP1Pg4cSHJf0/8An8p7XdthGz0IjMt7/dpJO3UH/grskw73B7pVc1vqqkcc7z3aJCLWA42PNik0Ebgu/fsO4PDGbC1pErAQKLziq5w6O5JqtFFnU1EbFSNpZ2DbiPhLJHv+9cCktg+9Ztq8jTqpStrpSOCvETEXICKWR8Smam5LXTVx7AosKRiuT8uKThMRG4G3gf6S+gLnAhe1os6OpBptBMmTAP4oaXb6+JiOrNVtlI4bJulJSX+WNL5g+voSdXYk1WijRtek3VT/1gkSTSXt9BEgJP1B0hxJ3y2YvirbUru6j6OGSj7apIVpLgIuj4iVTbbVcursSKrRRgAHR8QraZ/0nyQ9FxEPVR5uLippo1eB3SJiuaSxwHRJe5VZZ0fS5m0UEe8AJ0fEy5K2IekS/SLJf9QdVSXt1B34O2BfYDUwQ9Js4J0y6myVrpo4ynm0SeM09ZK6A9sBbwD7A5+T9EOgH/CupLXA7DLq7EjavI0i4sqIeAUgIpZKupvkEL2jJo5Wt1HadbAOICJmS3qR5D/H+rSelursSKrRRnUR8XJavkLSTSTbUUdOHJXsb/XAnyNiGYCk3wFjgBuo1raU90mhPH5IEuZCYBjvn4jaq8k0Z7H5iajbitRzIe+fHC9ZZ0f6qVIb9QW2Kfj7UZKnIee+vrVuI2Ag75/A3B14GdghHX4COID3T2genfe6tqc2SusckJb3IOnvPzPvdc2xnbYH5gBbpfXcDxxTzW0p9wbL8YM6Gvg/kisZ/iUt+3fgM+nfvUmuCFpAcmXC7kXqeO9Lsbk6O/JPW7dRuvPPTX/mdeU2Aj6btsHcdKf/dEGd44Bn0jqvJH3CQ0f9aes2IvmnYzbJlUTzSN8Ymvd65tVO6bhT0rZ4BvhhtbclP3LEzMwy6apXVZmZWSs5cZiZWSZOHGZmlokTh5mZZeLEYWZmmThxmJlZJk4cZmaWiROHdRnpOxw+2aTsW5J+1sI8K2sQ1zclzZd0Y7WXZdYWnDisK7mZ5FENhU5Ky/P0dZJHQZyccxxmZXHisK7kDuBYSb0gecMcsAvwsKTp6aPe5xV73Hv6NrpnCobPkXRhwfApkh5PH/P9S0nditTxbUnPpD/fSst+QfIolnslnV1knlvSt77NkrRY0jGVNoJZpZw4rMuIiOUkz/g5Ki1qfItaAKdHxFiSZ/t8U1L/ZqrZgqQRwIkkj4wfBWwCTm4yzVjgSyRPDj4A+Kqk0RFxJskTSz8eEZcXqX4fYGFE7J/WOaXsFTarEicO62oKu6sKu6m+KWku8BjJo6uHZ6jzcGAs8ISkp9Lh3ZtM83fA3RGxKiJWAncBTV9MtBlJfYABvP9CrGdJnoTaOP5WSd/JEKdZm+iq7+Owrms6cJmkMUCfiJgj6VDgCODAiFgt6UGSJ5EW2sjm/2gVjhdwXUSc38JyW/OGupHACxGxNh0eQ/KkWCRNBO5L4zarKR9xWJeS/rf/IHA17x9tbAe8mSaNPUi6kpp6HdhRUv/0HMmxBeNmkLy4akcASTtIGtJk/oeASZK2Sl+texwws0S4+wC7SeqdznMRcLmk3sAJEfHfaexmNeUjDuuKbibpKmrssvo9cKakvwLPk3RXbSYiNkj6d2AW8DfguYJxz0r6V5J3qX8A2EDy0p3FBdPMkXQtyTkWgF9HxJMl4twHuJEk0W0L/EdEPCLp34Ct0xPre0nqExFrsjSAWSX8Pg6zdkrSQ8BXI+L5grLdgCkR8eV0eArw+4iYlVOY1gU5cZi1U5JeBgZHxLt5x2JWyInDzMwy8clxMzPLxInDzMwyceIwM7NMnDjMzCwTJw4zM8vEicPMzDJx4jAzs0z+P9FuxmqQGtpDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.vlines(p_true, 0, 170, linestyles=\"--\",\n",
    "          label='true $p_A$ (unknown)')\n",
    "plt.hist(mcmc.trace('p')[:],\n",
    "        bins=35, histtype=\"stepfilled\", normed=True)\n",
    "plt.title('Posterior distribution of $p_A$, '\n",
    "         'the true effectiveness of site A')\n",
    "plt.xlabel('Value of $p_A$')\n",
    "plt.ylabel('Density')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## AB Test 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_true_a = 0.05\n",
    "p_true_b = 0.04\n",
    "\n",
    "n_a = 1500\n",
    "n_b = 750\n",
    "\n",
    "observations_a = pm.rbernoulli(p_true_a, n_a)\n",
    "observations_b = pm.rbernoulli(p_true_b, n_b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " [-----------------100%-----------------] 40000 of 40000 complete in 6.8 sec"
     ]
    }
   ],
   "source": [
    "p_a = pm.Uniform('p_a', lower=0, upper=1)\n",
    "p_b = pm.Uniform('p_b', lower=0, upper=1)\n",
    "\n",
    "rate_conv_a = pm.Bernoulli('obs_a', p_a, value=observations_a, observed=True)\n",
    "rate_conv_b = pm.Bernoulli('obs_b', p_b, value=observations_b, observed=True)\n",
    "# delta = pm.Bernoulli('delta', p_a - p_b)\n",
    "@pm.deterministic\n",
    "def delta(p_a=p_a, p_b=p_b):\n",
    "    return p_a - p_b\n",
    "\n",
    "mcmc = pm.MCMC([p_a, p_b, rate_conv_a, rate_conv_b, delta])\n",
    "mcmc.sample(40000, 10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize']=[12.5, 10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuYAAAJCCAYAAACMFAKrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xt0lPW97/HPtyEhJFjBQG011sRTqgE2BAjBGzEVpaJuiWgLVCoe2lKrPXsX27WrXatVN3ZjXRTZnK1lYYV6QcFCRQ9aWTWFFthKSbhUKVo0O2DUlhAscgn37/kjYxaQ20xmJs8zw/u11qzMzHP7PPhz8s1vfs/vMXcXAAAAgGB9KugAAAAAACjMAQAAgFCgMAcAAABCgMIcAAAACAEKcwAAACAEKMwBAACAEKAwBwAAAEKAwhwAAAAIAQpzAAAAIAS6deXB+vTp4wUFBV15SAAAACCpqqurd7l733j306WFeUFBgaqqqrrykAAAAEBSmdn2ROyHoSwAAABACFCYAwAAACEQVWFuZtPMbIuZvWlmz5pZtpkVmtk6M9tmZovNLCvZYQEAAIB01eEYczM7V9K/SOrv7o1m9pykCZKulfSwuy8ys7mSviHpF0lNCwAAEIMjR46orq5OBw8eDDoK0kB2drby8/OVmZmZlP1He/FnN0k9zOyIpBxJH0q6UtLXIsufkHSfKMwBAECI1NXV6YwzzlBBQYHMLOg4SGHuroaGBtXV1amwsDApx+hwKIu7vy9ppqQdairI90iqlvQPdz8aWa1O0rlJSQgAANBJBw8eVF5eHkU54mZmysvLS+q3Lx0W5mbWW9JYSYWSzpGUK2lMK6t6G9tPNbMqM6uqr6+PJysAIEpjxozRmDGtfVQDpx+KciRKsttSNENZrpL0P+5eHwn0G0mXSuplZt0iveb5kj5obWN3nydpniSVlJS0WrwDABLrt7/9bdARAAAximZWlh2SLjazHGv6M2GUpL9IWinp5sg6kyW9kJyIAAAAQPqLZoz5OklLJG2Q9EZkm3mSfijpLjN7R1KepMeTmBMAEIPp06dr+vTpQccAAMQgqnnM3f1ed7/I3Qe6+9fd/ZC717h7qbt/wd2/4u6Hkh0WABCdyspKVVZWBh0DgKTy8nKtWLHipPdmz56tO+64o81tevbsmexYmjNnjoqKinTLLbckbJ+NjY264oordOzYsU5tX1tbq4EDByYsT2ccPnxYZWVlOnr0aMcrJxh3/gQAAEiiiRMnatGiRSe9t2jRIk2cODGgRE0effRRvfzyy1q4cGHC9jl//nyNGzdOGRkZCdtnV8vKytKoUaO0ePHiLj82hTkAAEAS3XzzzVq+fLkOHWoaXFBbW6sPPvhAl19+uSoqKjRs2DANGDBA8+bNa7HtqT3IM2fO1H333df8+umnn1ZpaamKi4v17W9/u9We6lmzZmngwIEaOHCgZs+eLUm6/fbbVVNToxtuuEEPP/xwi20mTJig8ePHa8SIETr//PP10ksvRXWuCxcu1NixYzvMXltbq6KiIn3rW9/SgAEDNHr0aDU2Np60r5qaGg0ZMkTr169vd/3Wzu+hhx7SnDlzJEnTpk3TlVdeKanp28RJkyZ1ePyKioqE/sESLQpzAABw2igvL2/xePTRRyVJBw4caHX5r371K0nSrl27WiyLRl5enkpLS/XKK69IauotHz9+vMxM8+fPV3V1taqqqjRnzhw1NDREfS5bt27V4sWLtXbtWm3atEkZGRktisnq6motWLBA69at0+uvv67HHntMGzdu1Ny5c3XOOedo5cqVmjZtWot9b968WRdccIHWrVunhQsX6v777+8wz+HDh1VTU6OCgoKo8m/btk133nmntmzZol69emnp0qXNy95++23ddNNNWrBggYYPH97m+m2dX1lZmVavXi1Jqqqq0r59+3TkyBGtWbNGI0eO7PD4AwcO1Pr166M6j0SK9s6fAIAQKri79V6s+g+PdHESAO35ZDjL2LFjtWjRIs2fP19S0zjv559/XpL03nvvadu2bcrLy4tqn5WVlaqurm4uXBsbG/WZz3zmpHXWrFmjG2+8Ubm5uZKkcePGafXq1RoyZEib+21sbNSuXbt07733SpL69++vjz76qHn5+PHjVVpaqu9///snbbdr1y716tUrquySVFhYqOLiYknSsGHDVFtbq8svv1z19fUaO3asli5dqgEDBrS7fkNDQ6vn953vfEfV1dXau3evunfvrqFDh6qqqkqrV69u7klvbX+fyMjIUFZWlvbu3aszzjgj6nOKF4U5AKShvjf+KOgIQCitWrWqzWU5OTntLu/Tp0+7y9tTUVGhu+66Sxs2bFBjY6OGDh2qVatW6dVXX9Vrr72mnJwclZeXt7irZLdu3XT8+PHm1ycud3dNnjxZM2bMaPO47rHfQubNN99Uv379lJ2dLUnasGGDBg8eLEl64YUXdP311+vVV19tsV2PHj1Oytdedknq3r178/OMjIzmoSRnnnmmzjvvPK1du/akwry19ds6v8zMTBUUFGjBggW69NJLNWjQIK1cuVLvvvuuioqKtH379jaP/4lDhw41/xt0FYayAAAAJFnPnj1VXl6uKVOmNF/0uWfPHvXu3Vs5OTl666239Prrr7fY7uyzz9bOnTvV0NCgQ4cOafny5c3LRo0apSVLlmjnzp2SpN27d2v79u0nbV9WVqZly5bpwIED2r9/v55//vnmoRxt2bx5s3bs2KGDBw9q//79uvfeezVt2jQdPHhQv/71r/X1r39de/bsabFd7969dezYseYCvL3s7cnKytKyZcv05JNP6plnnml33fbOr6ysTDNnzlRZWZlGjhypuXPnqri4OKq7dzY0NKhv377KzMyMKnOi0GMOAGnooz/8KvLsuiBjADjBxIkTNW7cuOYZWq655hrNnTtXgwYN0oUXXqiLL764xTaZmZn6yU9+ohEjRqiwsFAXXXRR87L+/fvrgQce0OjRo3X8+HFlZmbqkUce0fnnn9+8ztChQ3XbbbeptLRUkvTNb36z3WEsUlNhfsstt6i8vFwff/yxfvSjH+myyy7T9OnTtW/fPt1+++3asmWLGhsb1aNHj5O2HT16tNasWaOrrrqq3ewdyc3N1fLly3X11VcrNze3ucf+VO2d38iRI/XTn/5Ul1xyiXJzc5Wdnd3hHyWfWLlypa699tqo8yaKdeYrjs4qKSnxqqqqLjseAKS7tsaY/+2ZuyVJB3e80ZVxgNDZunWrioqKgo6RUsrKyvTYY4/pwgsvbH5vx44duv/++/X44033k7z//vt1zTXXaMSIESdtu3HjRs2aNUtPPfVUl2ZOtHHjxmnGjBkn/Rt8orU2ZWbV7l4S73HpMQcAAECzd999V/369Tvpvc9//vPNRbmk5gtDTzVkyBB96Utf0rFjx1J2LvPDhw+roqKi1aI82SjMAQAA0Oz999+Pa/spU6YkKEkwsrKydOuttwZybC7+BAAAAEKAHnMASEPdzugjqe0x6CeqfZALRAEgDCjMASAN9fnnHwQdAQgNd49qijygI8meNIWhLAAAIG1lZ2eroaEh6QUV0p+7q6GhIak3HaLHHADS0O5X50mSzrpqasBJgGDl5+errq5O9fX1QUdBGsjOzlZ+fn7S9k9hDgBp6PDOmqAjAKGQmZmpwsLCoGMAUWEoCwAAABACFOYAAABACFCYAwAAACHAGHMACEgy5xjPPOvcTm0HAAgOhTkApKG8a/5P1Ot29AcCNyACgK7BUBYAAAAgBCjMASANNbzyf9Xwyv8NOgYAIAYMZQGANHRk9/tBRwAAxIgecwAAACAEKMwBAACAEKAwBwAAAEKAMeYAkIayPnNB0BEAADGiMAeANHTWVVODjgAAiFFUhbmZ9ZL0S0kDJbmkKZLelrRYUoGkWklfdfePkpISAE5T0dwdFACQHqIdY/6fkl5x94skDZa0VdLdkirdvZ+kyshrAEAI7Pp/M7Xr/80MOgYAIAYdFuZm9mlJZZIelyR3P+zu/5A0VtITkdWekFSRrJAAgNgc3btLR/fuCjoGACAG0fSYXyCpXtICM9toZr80s1xJZ7v7h5IU+fmZJOYEAAAA0lo0hXk3SUMl/cLdh0jarxiGrZjZVDOrMrOq+vr6TsYEAAAA0ls0hXmdpDp3Xxd5vURNhfrfzexzkhT5ubO1jd19nruXuHtJ3759E5EZAAAASDsdFubu/jdJ75nZhZG3Rkn6i6QXJU2OvDdZ0gtJSQgAiFn3cy9S93MvCjoGACAG0c5j/n8kLTSzLEk1kv63mor658zsG5J2SPpKciICAGLV+4rbgo4AAIhRVIW5u2+SVNLKolGJjQMAAACcnqKdxxwAkELqn/8P1T//H0HHAADEINqhLACAFHKs8eOE7auju4/WPnhdwo4FAKczeswBAACAEKDHHACSpKOeZgAATkSPOQAAABAC9JgDQBrKPn9w0BEAADGiMAeANNTrsolBRwAAxIihLAAAAEAIUJgDQBr6+3P36u/P3Rt0DABADBjKAgBpyI8eCjoCACBG9JgDAAAAIUBhDgAAAIQAhTkAAAAQAowxB4A01ON/lXbZsaK5w2ntg9d1QRIASG0U5gDQCdEUo0E6c8S4oCMAAGLEUBYAAAAgBCjMASAN/e2Zu/W3Z+4OOgYAIAYU5gAAAEAIUJgDAAAAIUBhDgAAAIQAhTkAAAAQAkyXCACtCPt0iB3JvWhk0BEAADGiMAeANHTGUG7oAwCphqEsAJCGjh85qONHDgYdAwAQA3rMASAN7fz1fZKkz37twWCDAACiRo85AAAAEAIU5gAAAEAIUJgDAAAAIUBhDgAAAIRA1Bd/mlmGpCpJ77v79WZWKGmRpLMkbZD0dXc/nJyYAIBY9Pynq4KOAACIUSw95v8qaesJr38m6WF37yfpI0nfSGQwAEDn9fynqyjOASDFRFWYm1m+pOsk/TLy2iRdKWlJZJUnJFUkIyAAIHbHDuzRsQN7go4BAIhBtD3msyX9m6Tjkdd5kv7h7kcjr+sknZvgbACATqpfNkP1y2YEHQMAEIMOC3Mzu17STnevPvHtVlb1NrafamZVZlZVX1/fyZgAAABAeoumx/wySTeYWa2aLva8Uk096L3M7JOLR/MlfdDaxu4+z91L3L2kb9++CYgMAAAApJ8OC3N3v8fd8929QNIESb9391skrZR0c2S1yZJeSFpKAAAAIM1FPV1iK34oaZGZPSBpo6THExMJAJBuCu5+qd3ltQ9e10VJACC8YirM3X2VpFWR5zWSShMfCQAQrzOGXBt0BABAjOLpMQcAhFRuUVnQEQAAMYrlBkMAgBRx9ON6Hf2YmbAAIJVQmANAGtq1/OfatfznQccAAMSAwhwAAAAIAQpzAAAAIAQozAEAAIAQoDAHAAAAQoDpEgEgDX269MagIwAAYkRhDgBpKOcLI4KOAACIEUNZACANHWmo05GGuqBjAABiQI85AKShhhX/JUn67NceDDhJdArufqnd5bUPXtdFSQAgOBTmAE5LHRWCAAB0NYayAAAAACFAYQ4AAACEAIU5AAAAEAKMMQeANHTmpROCjgAAiBGFOQCkoR4FxUFHAADEiKEsAJCGDv+9Rof/XhN0DABADCjMASAN7a6cp92V84KOAQCIAYU5AAAAEAIU5gAAAEAIUJgDAAAAIUBhDgAAAIQA0yUCQBrqVTY56AgAgBhRmANAGsrOLwo6AgAgRhTmAJCGDtZtlZQ+BXrB3S91uE7tg9d1QRIASB7GmANAGvrHH5/QP/74RNAxAAAxoMccQNqJpncVAICwocccAAAACAEKcwAAACAEOizMzew8M1tpZlvNbIuZ/Wvk/bPM7Hdmti3ys3fy4wIAAADpKZox5kclfd/dN5jZGZKqzex3km6TVOnuD5rZ3ZLulvTD5EUFAETrrFFTg44AAIhRh4W5u38o6cPI871mtlXSuZLGSiqPrPaEpFWiMAeAUMg6+4KgIwAAYhTTGHMzK5A0RNI6SWdHivZPivfPtLHNVDOrMrOq+vr6+NICAKLSWLtJjbWbgo4BAIhB1IW5mfWUtFTS99z942i3c/d57l7i7iV9+/btTEYAQIz2/Pci7fnvRUHHAADEIKrC3Mwy1VSUL3T330Te/ruZfS6y/HOSdiYnIgAAAJD+opmVxSQ9Lmmru886YdGLkiZHnk+W9ELi4wEAAACnh2hmZblM0tclvWFmnwxY/JGkByU9Z2bfkLRD0leSExEAAABIf9HMyrJGkrWxeFRi4wAAAACnp2h6zAEAKSbvy98NOgIAIEYU5gCQhjLz8oOOAACIEYU5gJRTcPdLQUcIvQPvrJMk5XxhRMBJAADRojAHgDT08Z+el0RhDgCpJKY7fwIAAABIDgpzAAAAIAQYygIASAsdXXtQ++B1XZQEADqHHnMAAAAgBOgxB4A01Of67wcdIXToUQcQdhTmAJCGun26b9ARAAAxYigLAKSh/Vv/qP1b/xh0DABADOgxB4A0tHfjy5Kk3KKygJMAAKJFjzkAAAAQAhTmAAAAQAhQmAMAAAAhQGEOAAAAhAAXfwJAGupbcU/QEVJOR/OcS8x1DiC5KMwBIA1l5JwZdAQAQIwYygIAaWjfG69q3xuvBh0DABADCnMASEMU5gCQehjKAiB0ohnrCwBAuqHHHAAAAAgBeswBdCl6w5HKOmq/zNoCIB4U5gAAJAiFO4B4UJgDQBr6zFfuCzoCACBGFOYAkIY+lZkddAQAQIwozAEkFGPIw2Hvhqb/DmcMZehEmHB3UQDtYVYWAEhD+99arf1vrQ46BgAgBhTmAAAAQAjENZTFzK6R9J+SMiT90t0fTEgqAKHFUBUAAJKj04W5mWVIekTS1ZLqJK03sxfd/S+JCgcAwOmGKReB01c8Pealkt5x9xpJMrNFksZKojAHUhS94UD4JeL/U4p7IJziKczPlfTeCa/rJI04dSUzmyppauTlITN7M45j4vTRR9KuoEMgZdBe2rD9Z9cHHSFsaCuS7GdBJ0gZtBdE68JE7CSewtxaec9bvOE+T9I8STKzKncvieOYOE3QVhAL2guiRVtBLGgviJaZVSViP/HMylIn6bwTXudL+iC+OAAAAMDpKZ7CfL2kfmZWaGZZkiZIejExsQAAAIDTS6eHsrj7UTP7rqQVapoucb67b+lgs3mdPR5OO7QVxIL2gmjRVhAL2guilZC2Yu4thoUDAAAA6GLc+RMAAAAIAQpzAAAAIAQSUpib2TVm9raZvWNmd7eyvLuZLY4sX2dmBScsuyfy/ttm9uVE5EG4dba9mNnVZlZtZm9Efl7Z1dnRteL5bIks/7yZ7TOzH3RVZgQnzt9Fg8zsNTPbEvmMye7K7OhacfweyjSzJyJtZKuZ3dPV2dH1omgvZWa2wcyOmtnNpyybbGbbIo/JHR7M3eN6qOnCz3clXSApS9JmSf1PWecOSXMjzydIWhx53j+yfndJhZH9ZMSbiUd4H3G2lyGSzok8Hyjp/aDPh0c428oJy5dK+rWkHwR9PjzC217UNBHCnyUNjrzO43dR+j7ibCtfk7Qo8jxHUq2kgqDPiUfg7aVA0iBJT0q6+YT3z5JUE/nZO/K8d3vHS0SPeamkd9y9xt0PS1okaewp64yV9ETk+RJJo8zMIu8vcvdD7v4/kt6J7A/pq9Ptxd03uvsnc+VvkZRtZt27JDWCEM9ni8ysQk0fgh3NFoX0EE97GS3pz+6+WZLcvcHdj3VRbnS9eNqKS8o1s26Sekg6LOnjromNgHTYXty91t3/LOn4Kdt+WdLv3H23u38k6XeSrmnvYIkozM+V9N4Jr+si77W6jrsflbRHTT0S0WyL9BJPeznRTZI2uvuhJOVE8DrdVswsV9IPJd3fBTkRDvF8tnxRkpvZisjX0f/WBXkRnHjayhJJ+yV9KGmHpJnuvjvZgRGoeGrVmLft9DzmJ7BW3jt1Dsa21olmW6SXeNpL00KzAZJ+pqZeLqSveNrK/ZIedvd9kQ50pL942ks3SZdLGi7pgKRKM6t298rERkRIxNNWSiUdk3SOmoYmrDazV929JrERESLx1Koxb5uIHvM6Seed8Dpf0gdtrRP5+udMSbuj3BbpJZ72IjPLl/S8pFvd/d2kp0WQ4mkrIyQ9ZGa1kr4n6UeRG6IhfcX7u+gP7r7L3Q9IelnS0KQnRlDiaStfk/SKux9x952S1koqSXpiBCmeWjXmbRNRmK+X1M/MCs0sS00XSbx4yjovSvrkStSbJf3em0bFvyhpQuTq50JJ/ST9KQGZEF6dbi9m1kvSS5Lucfe1XZYYQel0W3H3ke5e4O4FkmZL+g93/6+uCo5AxPO7aIWkQWaWEynCrpD0ly7Kja4XT1vZIelKa5Ir6WJJb3VRbgQjmvbSlhWSRptZbzPrraZv+le0t0HcQ1nc/WikJ2qFmq5cne/uW8zs3yVVufuLkh6X9JSZvaOmvzgnRLbdYmbPqekD8KikO7ngJr3F014kfVfSFyT92Mx+HHlvdKTXAmkmzraC00ycv4s+MrNZavoF7JJedveXAjkRJF2cny2PSFog6U01DVNYELnoD2kqmvZiZsPV9G1+b0n/bGb3u/sAd99tZtPV9NkiSf/e0TUJFpnOBQAAAECAuPMnAAAAEAIU5gAAAEAIUJgDAAAAIUBhDgAAAIQAhTkAAAAQAhTmAAAAQAhQmAMAAAAhQGEOAAAAhACFOQAAABACFOYAAABACFCYAwAAACFAYQ4AAACEAIU5AAAAEAIU5gAAAEAIUJgDAAAAIdCtKw/Wp08fLygo6MpDAgAAAElVXV29y937xrufLi3MCwoKVFVV1ZWHBAAAAJLKzLYnYj8MZQEAAABCgMIcAAAACIGoCnMzm2ZmW8zsTTN71syyzazQzNaZ2TYzW2xmWckOCwAAAKSrDseYm9m5kv5FUn93bzSz5yRNkHStpIfdfZGZzZX0DUm/SGpaAACAGBw5ckR1dXU6ePBg0FGQBrKzs5Wfn6/MzMyk7D/aiz+7SephZkck5Uj6UNKVkr4WWf6EpPtEYQ4AAEKkrq5OZ5xxhgoKCmRmQcdBCnN3NTQ0qK6uToWFhUk5RodDWdz9fUkzJe1QU0G+R1K1pH+4+9HIanWSzk1KQgAAgE46ePCg8vLyKMoRNzNTXl5eUr996bAwN7PeksZKKpR0jqRcSWNaWdXb2H6qmVWZWVV9fX08WQGghTFjxmjMmNY+kgCgCUU5EiXZbSmaoSxXSfofd6+PBPqNpEsl9TKzbpFe83xJH7S2sbvPkzRPkkpKSlot3gGgs377298GHQEAgISIZlaWHZIuNrMca/ozYZSkv0haKenmyDqTJb2QnIgAAABA+otmjPk6SUskbZD0RmSbeZJ+KOkuM3tHUp6kx5OYEwBaNX36dE2fPj3oGAAAxC2qeczd/V53v8jdB7r71939kLvXuHupu3/B3b/i7oeSHRYATlVZWanKysqgYwBAm8rLy7VixYqT3ps9e7buuOOONrfp2bNnsmNpzpw5Kioq0i233JKwfTY2NuqKK67QsWPHOrV9bW2tBg4cmLA8nXH48GGVlZXp6NGjHa+cYNz5EwAAIIkmTpyoRYsWnfTeokWLNHHixIASNXn00Uf18ssva+HChQnb5/z58zVu3DhlZGQkbJ9dLSsrS6NGjdLixYu7/NgU5gAAAEl08803a/ny5Tp0qGlwQW1trT744ANdfvnlqqio0LBhwzRgwADNmzevxban9iDPnDlT9913X/Prp59+WqWlpSouLta3v/3tVnuqZ82apYEDB2rgwIGaPXu2JOn2229XTU2NbrjhBj388MMttpkwYYLGjx+vESNG6Pzzz9dLL70U1bkuXLhQY8eO7TB7bW2tioqK9K1vfUsDBgzQ6NGj1djYeNK+ampqNGTIEK1fv77d9Vs7v4ceekhz5syRJE2bNk1XXnmlpKZvWSdNmtTh8SsqKhL6B0u0KMwBAMBpo7y8vMXj0UcflSQdOHCg1eW/+tWvJEm7du1qsSwaeXl5Ki0t1SuvvCKpqbd8/PjxMjPNnz9f1dXVqqqq0pw5c9TQ0BD1uWzdulWLFy/W2rVrtWnTJmVkZLQoJqurq7VgwQKtW7dOr7/+uh577DFt3LhRc+fO1TnnnKOVK1dq2rRpLfa9efNmXXDBBVq3bp0WLlyo+++/v8M8hw8fVk1NjQoKCqLKv23bNt15553asmWLevXqpaVLlzYve/vtt3XTTTdpwYIFGj58eJvrt3V+ZWVlWr16tSSpqqpK+/bt05EjR7RmzRqNHDmyw+MPHDhQ69evj+o8EonCHEBKy8vLU15eXtAxAKBdJw5nOXEYy5w5czR48GBdfPHFeu+997Rt27ao91lZWanq6moNHz5cxcXFqqysVE1NzUnrrFmzRjfeeKNyc3PVs2dPjRs3rrlgbUtjY6N27dqle++9V5LUv39/ffTRR5KkkpIS3Xnnnbriiiu0ZcuWk7bbtWuXevXqFXX+wsJCFRcXS5KGDRum2tpaSVJ9fb3Gjh2rp59+unl5W+u3dX7Dhg1TdXW19u7dq+7du+uSSy5RVVWVVq9e3VyYt3V8ScrIyFBWVpb27t0b9fkkQjTzmANAaJ3YwwEAHVm1alWby3Jyctpd3qdPn3aXt6eiokJ33XWXNmzYoMbGRg0dOlSrVq3Sq6++qtdee005OTkqLy9vcVfJbt266fjx482vT1zu7po8ebJmzJjR5nHdY7+FzJtvvql+/fopOztbkrRhwwYNHjxY7733nkpLS/XII49o1qxZqqur04ABA5q369Gjx0n52ssuSd27d29+npGR0TyU5Mwzz9R5552ntWvXnrT/1tZv6/wyMzNVUFCgBQsW6NJLL9WgQYO0cuVKvfvuuyoqKtL27dvbPP4nDh061Pxv0FXoMQcAAEiynj17qry8XFOmTGnuLd+zZ4969+6tnJwcvfXWW3r99ddbbHf22Wdr586damho0KFDh7R8+fLmZaNGjdKSJUu0c+dOSdLu3bu1ffv2k7YvKyvTsmXLdODAAe3fv1/PP/98c49xWzZv3qwdO3bo4MGD2r9/v+69915NmzZN1dXV+utf/6opU6bo979gCZ1gAAAgAElEQVT/vb785S+ftF3v3r117Nix5gK8veztycrK0rJly/Tkk0/qmWeeaXfd9s6vrKxMM2fOVFlZmUaOHKm5c+equLg4qrt3NjQ0qG/fvsrMzIwqc6JQmANIaffcc4/uueeeoGMAQIcmTpyozZs3a8KECZKka665RkePHtWgQYP04x//WBdffHGLbTIzM/WTn/xEI0aM0PXXX6+LLrqoeVn//v31wAMPaPTo0Ro0aJCuvvpqffjhhydtP3ToUN12220qLS3ViBEj9M1vflNDhgxpN+fmzZt1yy23qLy8XMOHD9d3vvMdXXbZZaqurtbPf/5zzZ8/X927d9f+/ftbbDt69GitWbOmw+wdyc3N1fLly/Xwww/rhRfavodle+c3cuRIffjhh7rkkkt09tlnKzs7u8M/Sj6xcuVKXXvttVHnTRTrzFccnVVSUuJVVVVddjwA6e+Ti686+/UygPS2detWFRUVBR0jpZSVlemxxx7ThRdeeNL71157rc4//3x96lOfUu/evfXAAw+02Hbjxo2aNWuWnnrqqa6KmxTjxo3TjBkzWvwbSK23KTOrdveSeI/LGHMAAAA0e/fdd9WvX78W77/88ssdbjtkyBB96Utf0rFjx1J2LvPDhw+roqKi1aI82SjMAQAA0Oz999+Pa/spU6YkKEkwsrKydOuttwZybMaYAwAAACFAjzmAlJafnx90BAAAEoLCHEBKe/rpp4OOACDk3D2qKfKAjiR70hSGsgAAgLSVnZ2thoaGpBdUSH/uroaGhqTedIgecwAp7Xvf+54kafbs2QEnARBG+fn5qqurU319fdBRkAays7OTOoSSwhxAStu0aVPQEQCEWGZmpgoLC4OOAUSFoSwAAABACFCYAwAAACFAYQ4AAACEAGPMAaS0L37xi0FHAAAgISjMAaS0efPmBR0BAICEYCgLAAAAEAIU5gBS2tSpUzV16tSgYwAAEDeGsgBIaX/961+DjgAAQELQYw4AAACEAIU5AAAAEAIU5gAAAEAIMMYcQEorLi4OOgIAAAlBYQ4gpc2ePTvoCAAAJERUQ1nMrJeZLTGzt8xsq5ldYmZnmdnvzGxb5GfvZIcFAAAA0lW0Y8z/U9Ir7n6RpMGStkq6W1Klu/eTVBl5DQBdatKkSZo0aVLQMQAAiFuHQ1nM7NOSyiTdJknufljSYTMbK6k8stoTklZJ+mEyQgJAW+rq6oKOAABAQkTTY36BpHpJC8xso5n90sxyJZ3t7h9KUuTnZ5KYEwAAAEhr0RTm3SQNlfQLdx8iab9iGLZiZlPNrMrMqurr6zsZEwAAAEhv0RTmdZLq3H1d5PUSNRXqfzezz0lS5OfO1jZ293nuXuLuJX379k1EZgAAACDtdDjG3N3/ZmbvmdmF7v62pFGS/hJ5TJb0YOTnC0lNCgCtuOSSS4KOAABAQpi7d7ySWbGkX0rKklQj6X+rqbf9OUmfl7RD0lfcfXd7+ykpKfGqqqp4MwMAAAChYWbV7l4S736iusGQu2+S1NrBRsUbAAAAAAB3/gQQQgV3vxTVerUPXqebbrpJkrR06dJkRgIAIOkozAGktIaGhqAjAACQENHe+RMAAABAElGYAwAAACFAYQ4AAACEAGPMAaS0UaOYHAoAkB4ozAGktB//+MdBRwAAICEYygIAAACEAIU5gJQ2ZswYjRkzJugYAADEjaEsAFJaY2Nj0BEAAEgIeswBAACAEKDHHECXKLj7paTs8281DQnff+2D1yVsXwAARIsecwAAACAE6DEHkNJ6/K/SoCMAAJAQFOYAUtqZI8YFHQEAgIRgKAsAAAAQAhTmAFLa3565W3975u6gYwAAEDcKcwAAACAEKMwBAACAEODiTwA4RbRzojPfOQAgkegxBwAAAEKAHnMAKS33opFBRwAAICEozAGktDOGBjecJNohLxLDXgAAHWMoC4CUdvzIQR0/cjDoGAAAxI0ecwApbeev75MkffZrDwYbBACAONFjDgAAAIQAhTkAAAAQAhTmAAAAQAhQmAMAAAAhEPXFn2aWIalK0vvufr2ZFUpaJOksSRskfd3dDycnJgC0ruc/XRV0BAAAEiKWWVn+VdJWSZ+OvP6ZpIfdfZGZzZX0DUm/SHA+AGhXqhTmzHkOAOhIVENZzCxf0nWSfhl5bZKulLQkssoTkiqSERAA2nPswB4dO7An6BgAAMQt2jHmsyX9m6Tjkdd5kv7h7kcjr+sknZvgbADQofplM1S/bEbQMQAAiFuHhbmZXS9pp7tXn/h2K6t6G9tPNbMqM6uqr6/vZEwAAAAgvUXTY36ZpBvMrFZNF3teqaYe9F5m9skY9XxJH7S2sbvPc/cSdy/p27dvAiIDAAAA6afDiz/d/R5J90iSmZVL+oG732Jmv5Z0s5qK9cmSXkhiTgAhFMsFjQAAoH3xzGP+Q0l3mdk7ahpz/nhiIgEAAACnn1imS5S7r5K0KvK8RlJp4iMBQPTOGHJt0BESLtpvIphWEQDSS0yFOQCETW5RWdARAABIiHiGsgBA4I5+XK+jHzPjEwAg9VGYA0hpu5b/XLuW/zzoGAAAxI3CHAAAAAgBCnMAAAAgBCjMAQAAgBBgVhYALXDjIAAAuh6FOYCU9unSG4OOAABAQlCYA0hpOV8YEXQEAAASgjHmAFLakYY6HWmoCzoGAABxozAHkNIaVvyXGlb8V9AxAACIG4U5AAAAEAIU5gAAAEAIcPEnAKSoWKa1rH3wuiQmAQAkAj3mAAAAQAjQYw4gpZ156YSgIwAAkBAU5gBSWo+C4qAjAACQEAxlAZDSDv+9Rof/XhN0DAAA4kZhDiCl7a6cp92V84KOAQBA3CjMAQAAgBCgMAcAAABCgIs/AeA0wJznABB+FObAaSKWwgwAAHQ9CnMAKa1X2eSgIwAAkBAU5gBSWnZ+UdARAABICC7+BJDSDtZt1cG6rUHHAAAgbhTmAFLaP/74hP7xxyeCjgEAQNwozAEAAIAQoDAHAAAAQqDDwtzMzjOzlWa21cy2mNm/Rt4/y8x+Z2bbIj97Jz8uAAAAkJ6i6TE/Kun77l4k6WJJd5pZf0l3S6p0936SKiOvAQAAAHRCh9MluvuHkj6MPN9rZlslnStprKTyyGpPSFol6YdJSQkAbThr1NSgIwAAkBAxzWNuZgWShkhaJ+nsSNEud//QzD6T8HQA2sXdPKWssy8IOgIAAAkRdWFuZj0lLZX0PXf/2Myi3W6qpKmS9PnPf74zGQGgTY21myRJPQqKA06SPqL9g6/2weuSnAQATi9RzcpiZplqKsoXuvtvIm//3cw+F1n+OUk7W9vW3ee5e4m7l/Tt2zcRmQGg2Z7/XqQ9/70o6BgAAMQtmllZTNLjkra6+6wTFr0oaXLk+WRJLyQ+HgAAAHB6iGYoy2WSvi7pDTPbFHnvR5IelPScmX1D0g5JX0lORAAAACD9RTMryxpJbQ0oH5XYOAAAAMDpiTt/AgAAACEQ03SJABA2eV/+btARAABICApzACktMy8/6AgAACQEQ1kApLQD76zTgXfWBR0DAIC40WMOIKV9/KfnJUk5XxgRcJLTTyx3nuVmRADQMXrMAQAAgBCgMAcAAABCgMIcAAAACAEKcwAAACAEuPgTQErrc/33g44AAEBCUJgDSGndPt036AgAACQEhTkQQrFMQ3e627/1j5Kk3KKygJMAABAfCnMAKW3vxpclUZgDAFIfhTkAIOm4GREAdIxZWQAAAIAQoDAHAAAAQoChLACAUIl22AtDXgCkGwpzACmtb8U9QUcAACAhKMwBpLSMnDODjgAAQEJQmANdhLnJk2PfG69Kknr+01UBJwEAID5c/Akgpe1749Xm4hwAgFRGYQ4AAACEAIU5AAAAEAIU5gAAAEAIUJgDAAAAIWDu3mUHKykp8aqqqi47HpBszLQSvONHDkqSPpWZHXAShBk3IwKQTGZW7e4l8e6H6RIBpDQKckQjlj+iKeIBBIWhLABS2t4NL2nvBr65AACkPnrMAaS0/W+tliSdMZReTiRGtL3r9KwDSLS4CnMzu0bSf0rKkPRLd38wIamAADFuHAAABKHThbmZZUh6RNLVkuokrTezF939L4kKBwBAWDFuHUCixdNjXirpHXevkSQzWyRprCQKcwAATkARDyAa8RTm50p674TXdZJGxBcHiA3DTgCkm2R8rlHsA6khnsLcWnmvxaToZjZV0tTIy0Nm9mYcx8Tpo4+kXUGHQMros/1n19NeEI3T8rPFfhZ0gpR1WrYXdMqFidhJPIV5naTzTnidL+mDU1dy93mS5kmSmVUlYvJ1pD/aCmJBe0G0aCuIBe0F0TKzhNxBM555zNdL6mdmhWaWJWmCpBcTEQoAAAA43XS6x9zdj5rZdyWtUNN0ifPdfUvCkgEAAACnkbjmMXf3lyW9HMMm8+I5Hk4rtBXEgvaCaNFWEAvaC6KVkLZi7i2u1wQAAADQxeIZYw4AAAAgQRJSmJvZNWb2tpm9Y2Z3t7K8u5ktjixfZ2YFJyy7J/L+22b25UTkQbh1tr2Y2dVmVm1mb0R+XtnV2dG14vlsiSz/vJntM7MfdFVmBCfO30WDzOw1M9sS+YzJ7srs6Fpx/B7KNLMnIm1kq5nd09XZ0fWiaC9lZrbBzI6a2c2nLJtsZtsij8kdHszd43qo6cLPdyVdIClL0mZJ/U9Z5w5JcyPPJ0haHHneP7J+d0mFkf1kxJuJR3gfcbaXIZLOiTwfKOn9oM+HRzjbygnLl0r6taQfBH0+PMLbXtR0vdWfJQ2OvM7jd1H6PuJsK1+TtCjyPEdSraSCoM+JR+DtpUDSIElPSrr5hPfPklQT+dk78rx3e8dLRI95qaR33L3G3Q9LWiRp7CnrjJX0ROT5EkmjzMwi7y9y90Pu/j+S3onsD+mr0+3F3Te6+ydz5W+RlG1m3bskNYIQz2eLzKxCTR+CzBZ1eoinvYyW9Gd33yxJ7t7g7se6KDe6XjxtxSXlmlk3ST0kHZb0cdfERkA6bC/uXuvuf5Z0/JRtvyzpd+6+290/kvQ7Sde0d7BEFObnSnrvhNd1kfdaXcfdj0rao6YeiWi2RXqJp72c6CZJG939UJJyInidbitmlivph5Lu74KcCId4Plu+KMnNbEXk6+h/64K8CE48bWWJpP2SPpS0Q9JMd9+d7MAIVDy1aszbxjVdYoS18t6pU720tU402yK9xNNemhaaDZD0MzX1ciF9xdNW7pf0sLvvi3SgI/3F0166Sbpc0nBJByRVmlm1u1cmNiJCIp62UirpmKRz1DQ0YbWZveruNYmNiBCJp1aNedtE9JjXSTrvhNf5kj5oa53I1z9nStod5bZIL/G0F5lZvqTnJd3q7u8mPS2CFE9bGSHpITOrlfQ9ST+K3BAN6Sve30V/cPdd7n5ATffnGJr0xAhKPG3la5Jecfcj7r5T0lpJJUlPjCDFU6vGvG0iCvP1kvqZWaGZZanpIokXT1nnRUmfXIl6s6Tfe9Oo+BclTYhc/VwoqZ+kPyUgE8Kr0+3FzHpJeknSPe6+tssSIyidbivuPtLdC9y9QNJsSf/h7v/VVcERiHh+F62QNMjMciJF2BWS/tJFudH14mkrOyRdaU1yJV0s6a0uyo1gRNNe2rJC0mgz621mvdX0Tf+K9jaIeyiLux+N9EStUNOVq/PdfYuZ/bukKnd/UdLjkp4ys3fU9BfnhMi2W8zsOTV9AB6VdCcX3KS3eNqLpO9K+oKkH5vZjyPvjY70WiDNxNlWcJqJ83fRR2Y2S02/gF3Sy+7+UiAngqSL87PlEUkLJL2ppmEKCyIX/SFNRdNezGy4mr7N7y3pn83sfncf4O67zWy6mj5bJOnfO7omgTt/AgAAACHAnT8BAACAEKAwBwAAAEKAwhwAAAAIAQpzAAAAIAQozAEAAIAQoDAHAAAAQoDCHAAAAAgBCnMAAAAgBCjMAQAAgBCgMAcAAABCgMIcAAAACAEKcwAAACAEKMwBAACAEKAwBwAAAEKAwhwAAAAIAQpzAAAAIAS6deXB+vTp4wUFBV15SAAAACCpqqurd7l733j306WFeUFBgaqqqrrykAAAAEBSmdn2ROyHoSwAAABACFCYAwAAACEQVWFuZtPMbIuZvWlmz5pZtpkVmtk6M9tmZovNLCvZYQEAAIB01eEYczM7V9K/SOrv7o1m9pykCZKulfSwuy8ys7mSviHpF0lNCwAAEKcjR46orq5OBw8eDDoKUkx2drby8/OVmZmZlP1He/FnN0k9zOyIpBxJH0q6UtLXIsufkHSfKMwBAEDI1dXV6YwzzlBBQYHMLOg4SBHuroaGBtXV1amwsDApx+hwKIu7vy9ppqQdairI90iqlvQPdz8aWa1O0rlJSQgAAJBABw8eVF5eHkU5YmJmysvLS+o3LR0W5mbWW9JYSYWSzpGUK2lMK6t6G9tPNbMqM6uqr6+PJysAADEZM2aMxoxp7VcWTncU5eiMZLebaIayXCXpf9y9PhLoN5IuldTLzLpFes3zJX3Q2sbuPk/SPEkqKSlptXgHACAZfvvb3wYdAQCiFs2sLDskXWxmOdb0Z8IoSX+RtFLSzZF1Jkt6ITkRAQAAgPQXzRjzdZKWSNog6Y3INvMk/VDSXWb2jqQ8SY8nMScAADGbPn26pk+fHnQMoIXy8nKtWLHipPdmz56tO+64o81tevbsmexYmjNnjoqKinTLLbe0uvz555+Xmemtt95qcx+NjY264oordOzYsZiPX1tbq4EDB8a8XSIdPnxYZWVlOnr0aMcrJ1hU85i7+73ufpG7D3T3r7v7IXevcfdSd/+Cu3/F3Q8lOywAALGorKxUZWVl0DGAFiZOnKhFixad9N6iRYs0ceLEgBI1efTRR/Xyyy9r4cKFrS5/9tlnVVJS0iL7iebPn69x48YpIyMjWTGTKisrS6NGjdLixYu7/Njc+RMAAKCL3XzzzVq+fLkOHWrq16ytrdUHH3ygyy+/XBUVFRo2bJgGDBigefPmtdj21F7lmTNn6r777mt+/fTTT6u0tFTFxcX69re/3WrP9axZszRw4EANHDhQs2fPliTdfvvtqqmp0Q033KCHH364xTb79u3TH/7wBz3++ON69tln2zy3hQsXauzYse1mra2tVVFRkb71rW9pwIABGj16tBobG0/aT01NjYYMGaL169e3u35r5yJJDz30kObMmSNJmjZtmq688kpJTX+wT5o0qd19VlRUtPnHSTJRmAMAgNNaeXl5i8ejjz4qSTpw4ECry3/1q19Jknbt2tViWTTy8vJUWlqqV155RVJTb/n48eNlZpo/f76qq6tVVVWlOXPmqKGhIepz2bp1qxYvXqy1a9dq06ZNysjIaFFgVldXa8GCBVq3bp1ef/11PfbYY9q4caPmzp2rc845RytXrtS0adNa7HvZsmW66qqrNGjQIOXm5mrDhg0t1jl8+LBqampUUFDQYdZt27bpzjvv1JYtW9SrVy8tXbq0ednbb7+tm266SQsWLNDw4cPbXL+tc5GksrIyrV69WpJUVVWlffv26ciRI1qzZo1GjhzZboaBAwdq/fr1UfyLJxaFOQAAQABOHM5y4jCWOXPmaPDgwbr44ov13nvvadu2bVHvs7KyUtXV1Ro+fLiKi4tVWVmpmpqak9ZZs2aNbrzxRuXm5qpnz54aN25ccwHbnmeffVZf/epXJUlf/epXW+0137Vrl3r16hVV1sLCQhUXF0uShg0bptraWklSfX29xo4dq6effrp5eVvrt3cuw4YNU3V1tfbu3avu3bvrkksuUVVVlVavXt1cmLeVISMjQ1lZWdq7d29U55Io0d75EwCAlJOXlxd0BKSAVatWtbksJyen3eV9+vRpd3l7KioqdNddd2nDhg1qbGzU0KFDtWrVKr366qt67bXXlJOTo/Ly8hY3tOnWrZuOHz/e/PrE5e6uyZMna8aMGW0e1z322asbGhr0pz/9Sb/5zW8kSePHj9cVV1yhhx566KS5vXv06HFSnvaydu/evfl5RkZG8zCSM888U+edd57Wrl2rAQMGtLt+e+eSmZmpgoICLViwQJdeeqkGDRqklStX6t1331VRUZG2b9/eZgZJOnTokLKzs6P7B0oQeswBAGlr6dKlJ309DoRJz549VV5erilTpjT3lu/Zs0e9e/dWTk6O3nrrLb3++usttjv77LO1c+dONTQ06NChQ1q+fHnzslGjRmnJkiXauXOnJGn37t3avn37SduXlZVp2bJlOnDggPbv36/nn3++uQe5LUuWLNG1117bXMgWFhbqs5/9rNasWXPSer1799axY8eaC/D2srYlKytLy5Yt05NPPqlnnnmm3XU7OpeysjLNnDlTZWVlGjlypObOnavi4uIObxTU0NCgvn37KjMzs8O8iUSPOQAAQEAmTpyocePGNQ9pueaaazR37lwNGjRIF154oS6++OIW22RmZuonP/mJRowYocLCQl100UXNy/r3768HHnhAo0eP1vHjx5WZmalHHnlE559/fvM6Q4cO1W233abS0lJJ0je/+U0NGTKk3ZzPPvus/vznP580dryhoUHPPPNMi6J+9OjRWrNmja666qp2s7YnNzdXy5cv19VXX63c3FwNHjy41fU6OpeRI0fqpz/9qS655BLl5uYqOzu7wz9CJGnlypW69tpro8qaSNaZrzM6q6SkxKuqqrrseACA09s999wjSe1+rY/Tz9atW1VUVBR0jLS1ceNGzZo1S0899VTQUTpt3LhxmjFjhi688MIWy1prP2ZW7e4l8R6XHnMAQNp67bXXgo4AnHaGDBmiL33pSzp27FhKzmV++PBhVVRUtFqUJxuFOQAAABJqypQpQUfotKysLN16662BHJuLPwEAAIAQoDAHAAAAQoChLACAtJWfnx90BISUu3c4ZR5wqmRPmkJhDgBIW08//XTQERBC2dnZamhoUF5eHsU5oubuamhoSOpNhyjMAQDAaSU/P191dXWqr68POgpSTHZ2dlK/iaMwBwCkre9973uSpNmzZwecBGGSmZmpwsLCoGMALVCYAwDS1qZNm4KOAABRY1YWAAAAIAQozAEAAIAQoDAHAAAAQoAx5gCAtPXFL34x6AgAEDUKcwBA2po3b17QEQAgagxlAQAAAEKAwhwAkLamTp2qqVOnBh0DAKLCUBYAQNr661//GnQEAIgaPeYAAABACFCYAwAAACFAYQ4AAACEAGPMAQBpq7i4OOgIABA1CnMAQNqaPXt20BEAIGpRDWUxs15mtsTM3jKzrWZ2iZmdZWa/M7NtkZ+9kx0WAAAASFfRjjH/T0mvuPtFkgZL2irpbkmV7t5PUmXkNQAAoTFp0iRNmjQp6BgAEJUOh7KY/f/27jZWjvK8w/h1l9eadxubEkxiECjC6QdILQhNi1zcAgY3kNogqoqaNJUV0UhJG1SgFIm0tDJINLSlAlmklT+kdRJohBVoEVAsSKM6OYApIYbguG44gYBtyottwLFy98OOwzE5jsdnZndmZ6+ftNqdndk9/721x3v7Oc88G0cC5wBXAmTmTmBnRFwMzC8OWwmsAa7pR0hJkqZifHy86QiSVFqZEfOTgc3AP0XEkxFxV0QcBhyXmS8BFNez+phTkiRJ6rQyjfmBwIeBOzLzDGA7+zFtJSKWRcRYRIxt3rx5ijElSZKkbivTmI8D45m5tti+m16j/nJEHA9QXL8y2YMzc0VmzsvMeTNnzqwjsyRJktQ5+5xjnpk/iogXIuKDmfkcsAD4bnFZCiwvru/ta1JJkvbT2Wef3XQESSotMnPfB0WcDtwFHAxsBD5Bb7T9K8D7gR8Al2bmqz/veebNm5djY2NVM0uSJEmtERGPZ+a8qs9T6guGMnMdMNkPW1A1gCRJkqTy65hLkjR0Fi9ezOLFi5uOIUmllBoxlyRpGG3durXpCJJUmiPmkiRJUgvYmEuSJEktYGMuSZIktYBzzCVJnbVggYuHSRoeNuaSpM664YYbmo4gSaU5lUWSJElqARtzSVJnLVy4kIULFzYdQ5JKcSqLJKmz3nrrraYjSFJpjphLkiRJLWBjLkmSJLWAjbkkSZLUAs4xlyR11qJFi5qOIEml2ZhLkjrr6quvbjqCJJXmVBZJkiSpBWzMJUmdNX/+fObPn990DEkqxcZckiRJagEbc0mSJKkFbMwlSZKkFrAxlyRJklrA5RIlSZ112WWXNR1BkkqzMZckddZVV13VdARJKs2pLJKkztqxYwc7duxoOoYkleKIuSSpsy688EIA1qxZ02wQSSrBEXNJkiSpBWzMJUmSpBZwKoskjbA5195Xy/NsWn5RLc/TtjySNEiOmEuSJEktUHrEPCIOAMaAH2bmoog4CVgFTAeeAK7IzJ39iSlJarO6RrrrduWVVzYdQZJK25+pLJ8B1gNHFts3A1/IzFURcSfwSeCOmvNJkjRlNuaShkmpqSwRMRu4CLir2A7gXODu4pCVwCX9CChJ0lRt2bKFLVu2NB1DkkopO2J+G/CnwBHF9gzgtczcVWyPAyfUnE2SpEqWLFkCuI65pOGwzxHziFgEvJKZj0+8e5JDcy+PXxYRYxExtnnz5inGlCRJkrqtzFSWjwIfi4hN9E72PJfeCPrREbF7xH028OJkD87MFZk5LzPnzZw5s4bIkiRJUvfscypLZl4HXAcQEfOBqzPz9yLiq8ASes36UuDePuaUJKm03avE/Gjj1j2295froUsapCpfMHQNsCoibgKeBL5YTyRJ0r60dXnCrvELjyQN0n415pm5BlhT3N4InFl/JEmS6nHEGRc2HUGSSqsyYi5JUqsddto5TUeQpNJKrWMuSdIw2vXGZna94YpgkoaDjbkkqbO2fP1Wtnz91qZjSFIpNuaSJElSC9iYS5IkSS1gYy5JkiS1gI25JEmS1AIulyhJ6qwjz/x40xEkqTQbc0lSZ0075aymI0hSaU5lkdw7fSwAAAqqSURBVCR11o+3jvPjreNNx5CkUmzMJUmdtfWB29n6wO1Nx5CkUmzMJUmSpBawMZckSZJawJM/JUnqsznX3lfL82xaflEtzyOpnRwxlyRJklrAEXNJUmcd9auXNx1BkkqzMZckddYvzjm96QiSVJpTWSRJnbXz5Y3sfHlj0zEkqRQbc0lSZ7368ApefXhF0zEkqRSnskjSgNS1MockqZscMZckSZJawMZckiRJagEbc0mSJKkFnGMuSeqso89Z2nSEWvkNolK32ZhLkjrr0NmnNR1BkkpzKoskqbPeHl/P2+Prm44hSaXYmEuSOuu1R1fy2qMrm44hSaXYmEuSJEkt4BxzSdoHvxhIkjQI+xwxj4gTI+KRiFgfEc9ExGeK+6dHxIMR8XxxfUz/40qSJEndVGYqyy7gc5l5GvAR4I8iYi5wLfBwZp4KPFxsS5IkSZqCfU5lycyXgJeK229GxHrgBOBiYH5x2EpgDXBNX1JKkjQF0xcsazqCJJW2X3PMI2IOcAawFjiuaNrJzJciYlbt6SRJquDg405uOoIklVZ6VZaIOBy4B/hsZr6xH49bFhFjETG2efPmqWSUJGlK3tq0jrc2rWs6hiSVUqoxj4iD6DXlX8rMfy3ufjkiji/2Hw+8MtljM3NFZs7LzHkzZ86sI7MkSaW8/s1VvP7NVU3HkKRSyqzKEsAXgfWZ+TcTdq0Glha3lwL31h9PkiRJGg1l5ph/FLgCeDoidv898M+A5cBXIuKTwA+AS/sTUZIk1anOtfk3Lb+otueSRl2ZVVm+AcRedi+oN44kSZI0mvzmT0md5Td2SpKGiY25JKmzZpz/6aYjSFJpNuaSpM46aMbspiNIUmml1zGXJGnY7Niwlh0b1jYdQ5JKccRcktRZb3zrawBMO+WshpNI0r45Yi5JkiS1gCPmklrH1VSk4VHX76vroUuOmEuSJEmtYGMuSZIktYBTWSRJnXXsos81HUGSSrMxlyR11oFHzmw6giSV5lQWSVJnbV//KNvXP9p0DEkqxRFzSVJnvfnk/QAcdto5DSeRpH2zMZdUG5c5lCRp6pzKIkmSJLWAI+aSJKlxflGR5Ii5JEmS1AqOmEuSOmvmJdc1HUGSSrMxlyR11gHTjmo6giSV5lQWSVJnbXv6IbY9/VDTMSSpFBtzSVJn2ZhLGiY25pIkSVILOMdckiR1hssuapjZmEvyGzslSWoBp7JIkiRJLeCIuTTEHOmWfr5Zl97YdARJKs3GXJLUWb9w0KFNR5Ck0mzMJUmd9eYTvb8qHfFhT+TT/vEkUjXBOeaSpM7a/uxjbH/2saZjSFIplUbMI+IC4G+BA4C7MnN5LamkjnNuuCSNBkfetT+mPGIeEQcA/wAsBOYCvxsRc+sKJkmSJI2SKiPmZwIbMnMjQESsAi4GvltHMKkujk5LkoadI++joUpjfgLwwoTtceCsanHUBf7jIUlSO7VxsMrP+3dVacxjkvvyZw6KWAYsKzbfiYjvVPiZetexwJamQ/RT3DywH9X5Wg6QtayPtazPsf978yJrWQ/fl/WxloUaPu/bUMsP1PEkVRrzceDECduzgRffe1BmrgBWAETEWGbOq/AzVbCW9bGW9bGW9bGW9bGW9bGW9bGW9elSLassl/ht4NSIOCkiDgYuB1bXE0uSJEkaLVMeMc/MXRHxaeABessl/mNmPlNbMkmSJGmEVFrHPDPvB+7fj4esqPLztAdrWR9rWR9rWR9rWR9rWR9rWR9rWZ/O1DIyf+Z8TUmSJEkDVmWOuSRJkqSa1N6YR8T0iHgwIp4vro/Zy3FLi2Oej4ilk+xfPepLK1atZUT8e0Q8FRHPRMSdxbe1jqQqtYyIaRFxX0Q8W9Ry+WDTt0sN78u/iogXImLb4FK3S0RcEBHPRcSGiLh2kv2HRMSXi/1rI2LOhH3XFfc/FxHnDzJ320y1jhExIyIeiYhtEXH7oHO3UYVa/lZEPB4RTxfX5w46e9tUqOWZEbGuuDwVER8fdPa2qfJvZbH//cXv+dWDylxZZtZ6AW4Bri1uXwvcPMkx04GNxfUxxe1jJuz/HeCfge/UnW+YLlVrCRxZXAdwD3B5069pGGsJTAN+ozjmYOAxYGHTr2kYa1ns+whwPLCt6dfSUP0OAL4PnFy8n54C5r7nmKuAO4vblwNfLm7PLY4/BDipeJ4Dmn5NQ1jHw4BfAz4F3N70a2n6UrGWZwDvK27/MvDDpl/PENdyGnBgcft44JXd26N4qVLLCfvvAb4KXN306yl76cdUlouBlcXtlcAlkxxzPvBgZr6amf8HPAhcABARhwN/AtzUh2zDplItM/ON4pgD6b2pR/mEginXMjN3ZOYjAJm5E3iC3rr9o6rq+/K/MvOlgSRtpzOBDZm5sXg/raJX04km1vhuYEFERHH/qsx8JzP/B9hQPN8omnIdM3N7Zn4DeHtwcVutSi2fzMzd32HyDHBoRBwykNTtVKWWOzJzV3H/oYz2ZzZU+7eSiLiE3qDQUK0Y2I/G/LjdH7rF9axJjjkBeGHC9nhxH8BfArcCO/qQbdhUrSUR8QC9/3W/Se9NO6oq1xIgIo4Gfht4uE85h0EttRxhZWrz02OKD+rXgRklHzsqqtRRe6qrlouBJzPznT7lHAaVahkRZ0XEM8DTwKcmNOqjaMq1jIjDgGuAzw8gZ62mtFxiRDwE/NIku64v+xST3JcRcTpwSmb+8XvnCXVVv2r50xuZ50fEocCXgHPpjVx2Ur9rGREHAv8C/F1mbtz/hMOj37UccWVqs7djrOu7qtRRe6pcy4j4EHAzcF6NuYZRpVpm5lrgQxFxGrAyIv4tM0f1LztVavl54AuZua0YQB8aU2rMM/M397YvIl6OiOMz86WI2D1H6r3GgfkTtmcDa4CzgV+JiE1FtlkRsSYz59NRfazlxJ/xdkSspvcnn8425gOo5Qrg+cy8rYa4rTaI9+UIGwdOnLA9G3hxL8eMF/8hPAp4teRjR0WVOmpPlWoZEbOBrwG/n5nf73/cVqvlfZmZ6yNiO715+2P9i9tqVWp5FrAkIm4BjgZ+EhFvZ2brT/bux1SW1cDuFRiWAvdOcswDwHkRcUz0VnQ4D3ggM+/IzPdl5hx6J+Z8r8tNeQlTrmVEHF40TbtHei8Enh1A5raaci0BIuImer/wnx1A1rarVEvxbeDUiDgpIg6md8LS6vccM7HGS4D/yN6ZTKuBy4uVCE4CTgW+NaDcbVOljtrTlGtZTO+7D7guM/9zYInbq0otTyo+r4mIDwAfBDYNJnYrTbmWmfnrmTmn6CdvA/56GJpyoC+rssygN//2+eJ6enH/POCuCcf9Ab0TlzYAn5jkeebgqixTriVwHL039X/TO/Hh7xnts7ur1HI2vT+NrQfWFZc/bPo1DWMti/tvoTfK8ZPi+samX1MDNbwQ+B69FQeuL+77C+Bjxe1D6a0ksIFe433yhMdeXzzuOUZ4daAa6riJ3sjatuJ9OHfQ+dt0mWotgT8Htk/4t3EdMKvp1zOktbyi+LxeR2+RgUuafi1NX6r8jk94jhsZolVZ/OZPSZIkqQX85k9JkiSpBWzMJUmSpBawMZckSZJawMZckiRJagEbc0mSJKkFbMwlSZKkFrAxlyRJklrAxlySJElqgf8H0QS3rlx/ZgAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 900x720 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = plt.subplot(311)\n",
    "plt.vlines(p_true_a, 0, 80, linestyles='--', label='Value of $p_A$ (unknown)')\n",
    "plt.hist(mcmc.trace('p_a')[:], bins=35,\n",
    "        histtype='stepfilled', normed=True)\n",
    "plt.xlim([0, .1])\n",
    "plt.legend(loc='upper right')\n",
    "# plt.subtitle('Posterior distribution of $p_A$')\n",
    "\n",
    "ax = plt.subplot(312)\n",
    "plt.vlines(p_true_b, 0, 80, linestyles='--', label='Value of $p_B$ (unknown)')\n",
    "plt.hist(mcmc.trace('p_b')[:], bins=35,\n",
    "        histtype='stepfilled', normed=True)\n",
    "plt.xlim([0, .1])\n",
    "plt.legend(loc='upper right')\n",
    "# plt.subtitle('Posterior distribution of $p_A$')\n",
    "\n",
    "ax = plt.subplot(313)\n",
    "plt.vlines(0.01, 0, 80, linestyles='--', label='Value of $\\Delta$ (unknown)')\n",
    "plt.hist(mcmc.trace('delta')[:], bins=35,\n",
    "        histtype='stepfilled', normed=True)\n",
    "# plt.subtitle('Posterior distribution of $p_A$')\n",
    "plt.legend(loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "No handles with labels found to put in legend.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvQAAAJeCAYAAADbW5GeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xu0pGddJ/rvz3QgEgO5TIOhG0jUcA800DB4GBWJIBilUaOiiJHDWTkzgiPjcXFxzSzE0RE8DohHh1kZETMMTIAW7GhUZALoDOMCOklLCJcJhAS6CUmTi9whgd/5oyrjptPdu7p3V1U/uz+ftfaqeu+/qv2uvb/11PM+b3V3AACAMX3LsgsAAAAOn0APAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABibQAwDAwAR6gP2oqqur6olz3P8fV9VvzONYK/dXVddV1Q/MY9+LUlUPqqorq+rzVfUvF3nsw7Xy9wswbwI9sHTT0PnlqvpCVd1YVa+rqm9b4/7WFGK7+2Hd/e617ONIH2vW13Wkat/f8Rb5vqzwwiTv7u6Tuvv3FnxsgKOeQA8cLX6ku78tyaOTPDbJv15GEVW1YZnbj3rsOXtAkquXXQTA0UqgB44q3b0nyV8meXiSVNVDqurdVXXbtLvH0+9ct6peVFV7pl0xPlpV51TV65PcP8mfTVv8Xzhd975V9SdVtbeqPrGy68a0JfpFVfWBJF+sqg37tk6vUsddtt/3dVXVo6rqimmtb0pywj7brzzWobyuVWtP8tiq+lBV3Tr99mPlsbuqvmvF9B9X1W+scrwfWO09WbHur1TVB6rqH6rqTSuPvc+6+91XVb0zyfcn+f1pHQ/cz7bPqqr/Od3/Z6rqU1X1tP0d52Cq6sVV9fHp+/6hqvrRWV/LwX6/86oX4E4CPXBUqar7JfmhJFdW1fFJ/izJXye5d5JfTPKGmvSpflCS5yd5bHeflOQHk1zX3c9O8slMW/y7+7er6lum+/n7JJuSnJPkBVX1gysO/dNJzk1ycnffsU9NB6xjxu3vluRPk7w+yalJ3pLkxw/w+md+XbMce+pZ0/18Z5IHZoZvP1Y53qzvSZL8ZJKnJjkzySOS/Px+XvMB99XdT0ry35M8f1rH/9pPuWcneVSSP0lyvySvTvIfV3uN+/HxJN+T5F5JXpbkv1TV6au9lkP5/R7hegGSCPTA0eNPq+q2JP8jyd8k+XdJHp/k25K8vLu/1t3vTPLnmQTYrye5e5KHVtXx3X1dd3/8APt+bJKN3f3r0/1cm+Q/JXnminV+r7s/1d1f3s/2B6tj1u2PT/K73X17d29P8v4D1Hoor2uWYyfJ70+X35LkN/ep+3DN8p7cWdunp8f+syRb1rCvAzk7yau6+83dfXuS/5zk/vu0oD+sqr5eVZsPtJPufsu01m9095uSXJPkcTO8lkP5/a5ab1W9pareM/1554G+1QC4k0APHC2e0d0nd/cDuvsXpuH0vkk+1d3fWLHe9Uk2dffHkrwgya8luamqLq6q+x5g3w9Ict9pd47bph8cfjXJfVas86mD1HbAOg5h+z3d3ftsfxeH+LpmOfa+y6+f1rNWs7wnSfKZFc+/lElwP9x9HcjZSbavmL53ki9091dWzHtRJi3oDznQTqrq56pq14pz5OFJ/smKVQ70Wmb+/c5Y7wOTfF93PyHJ5zL5NgDggAR64Gj26ST3m3aZudP9k+xJku5+Y3f/s0wCeyd5xXSdzjf7VJJPTD8w3PlzUnf/0Ip19t1m5jpm2P6GJJuqqvbZfr8O4XXNcuxk0q1j5XE/vWL6S0nusWL622fc7yzvyawOe19VdXImr2/vitnnZXIdxp3rPCKT38Hbc4BAX1UPyORbm+cnOa27T07ywSS1v/X3MfPvd7V6p913vqW776jJSE+nJ9lfNyOA/02gB45m703yxSQvrKrjazL++Y8kuXjaj/5JVXX3JF9J8uVMuqskyY1JvmPFft6X5HPTi0e/taqOq6qHV9Vj11rHjNv/XZI7kvzL6UWrP5Zv7srxvx3i65rV86pqc1Wdmsk3E29asWxXkp+ZvidPTfJ9K5Yd7HhrfU+O1L7OzuT9+Znpe3tukl/I5BuOO/2rTD4UfSgHbqE/MZMPMHuTpKqek+mF2TOY+fc7Q70PSfLtVfXuTLr8/EF33zZjHcAxSqAHjlrd/bUkT0/ytCSfTfIfkvxcd38kk37mL5/O/0wm3RZ+dbrpbyX519OuE7/S3V/PJCBuSfKJ6TZ/mMnFj2utY9btfyyTiyhvTfJTSd56gNVnfl2zHHvqjZlccHrt9GflDY9+KZP35rZMLp790xXLDni8tb4nR3BfZyd5Q5LvzuS9fVkm3bc+lCRVtSXJEzJ5D34nyb4X7d5Zw4eS/PtMwvmN0/2+5xDqn/X3e9B6M/kQ8ZrufmKSByd5ySw1AMe2+uYufwAwjqp6TZL/1d2vOsDyi5P8i+6+dTr9vu4+UOv53M1Q728leX93v7WqzkhycXc/foElAgNarzchAeDYcHaSHftbUFWPSfLlO8P81Feq6rTuvnkh1d3VAeudeliSp1TV85J8I8k/X0hVwNDm2kJfVf8qyf+VSb/Eq5I8J5MLfC7OZKzeK5I8e/p1JQAckuloNFu6+7pl1zKL0eoFxjC3QF9VmzIZT/qh3f3lqnpzkr/I5IYxb+3ui6vqPyb5++5+zVyKAACAdW7eF8VuSPKtNbkN+j0yGdrrSfnH8XcvSvKMOdcAAADr1tz60Hf3nqr6nUxuHf7lTEZYuDzJbStuTb47B7hxSFVdkOSCJDnxxBMf8+AHP3hepQIAwMJdfvnln+3ujWvdz9wCfVWdkmRbkjMzGQ7tLZkMSbav/fb56e4Lk1yYJFu3bu2dO3fOqVIAAFi8qjrYXaVnNs8uNz+QyZ0Z93b37ZmMyft/JDl52gUnSTbnm+9YCAAAHIJ5BvpPJnl8Vd1jejvsczK5S9+7MrnNdZKcn4MP3wUAABzE3AJ9d783k4tfr8hkyMpvyaQLzYuS/HJVfSzJaUleO68aAABgvZvrjaW6+6VJXrrP7GuTLO0ufQAAcKTdfvvt2b17d77yla/cZdkJJ5yQzZs35/jjj5/Lsd0pFgAA1mj37t056aSTcsYZZ2TS23yiu3PzzTdn9+7dOfPMM+dy7HmPQw8AAOveV77ylZx22mnfFOaTpKpy2mmn7bfl/kgR6AEA4AjYN8yvNv9IEegBAGBgAj0AAAxMoAcAgCOguw9p/pEi0AMAwBqdcMIJufnmm+8S3u8c5eaEE06Y27ENWwkAAGu0efPm7N69O3v37r3LsjvHoZ8XgR4AANbo+OOPn9s486vR5QYAAAYm0AMAwMAEegAAGJhADwAAAxPoAQBgYAI9AAAMTKAHAICBCfQAADAwgR4AAAYm0AMAwMAEegAAGJhADwAAAxPoAQBgYAI9AAAMTKAHAICBbVh2AQAs3o5dew66fNuWTQuqBIC10kIPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABibQAwDAwAR6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABibQAwDAwAR6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABibQAwDAwAR6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIFtWHYBABx9duzas+o627ZsWkAlAKxGCz0AAAxMoAcAgIEJ9AAAMDCBHgAABibQAwDAwOYW6KvqQVW1a8XP56rqBVV1alW9o6qumT6eMq8aAABgvZtboO/uj3b3lu7ekuQxSb6U5G1JXpzksu4+K8ll02kAAOAwLKrLzTlJPt7d1yfZluSi6fyLkjxjQTUAAMC6s6hA/8wk/3X6/D7dfUOSTB/vvb8NquqCqtpZVTv37t27oDIBAGAscw/0VXW3JE9P8pZD2a67L+zurd29dePGjfMpDgAABreIFvqnJbmiu2+cTt9YVacnyfTxpgXUAAAA69IiAv1P5x+72yTJJUnOnz4/P8mOBdQAAADr0lwDfVXdI8mTk7x1xeyXJ3lyVV0zXfbyedYAAADr2YZ57ry7v5TktH3m3ZzJqDcAAMAauVMsAAAMTKAHAICBCfQAADAwgR4AAAYm0AMAwMAEegAAGJhADwAAAxPoAQBgYAI9AAAMTKAHAICBCfQAADAwgR4AAAa2YdkFAHBoduzas+o627ZsWkAlABwNtNADAMDABHoAABiYLjcAHJbVuv7o9gOwGFroAQBgYAI9AAAMTKAHAICBCfQAADAwgR4AAAYm0AMAwMAEegAAGJhADwAAAxPoAQBgYAI9AAAMTKAHAICBCfQAADAwgR4AAAYm0AMAwMAEegAAGJhADwAAAxPoAQBgYAI9AAAMTKAHAICBCfQAADAwgR4AAAYm0AMAwMAEegAAGJhADwAAAxPoAQBgYAI9AAAMTKAHAICBCfQAADAwgR4AAAYm0AMAwMAEegAAGJhADwAAAxPoAQBgYBuWXQAAR96OXXuWXQIAC6KFHgAABibQAwDAwAR6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABjbXQF9VJ1fV9qr6SFV9uKq+u6pOrap3VNU108dT5lkDAACsZ/NuoX91kr/q7gcneWSSDyd5cZLLuvusJJdNpwEAgMMwt0BfVfdM8r1JXpsk3f217r4tybYkF01XuyjJM+ZVAwAArHfzbKH/jiR7k7yuqq6sqj+sqhOT3Ke7b0iS6eO997dxVV1QVTuraufevXvnWCYAAIxrnoF+Q5JHJ3lNdz8qyRdzCN1ruvvC7t7a3Vs3btw4rxoBAGBo8wz0u5Ps7u73Tqe3ZxLwb6yq05Nk+njTHGsAAIB1bW6Bvrs/k+RTVfWg6axzknwoySVJzp/OOz/JjnnVAAAA692GOe//F5O8oaruluTaJM/J5EPEm6vquUk+meQn5lwDAACsW3MN9N29K8nW/Sw6Z57HBQCAY4U7xQIAwMAEegAAGJhADwAAAxPoAQBgYAI9AAAMTKAHAICBCfQAADAwgR4AAAYm0AMAwMAEegAAGJhADwAAAxPoAQBgYAI9AAAMTKAHAICBbVh2AQCsTzt27Tno8m1bNi2oEoD1TaAH1o+rts+23tnnzbcOAFggXW4AAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABibQAwDAwAR6AAAYmHHoAY4yq92QCQBW0kIPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABibQAwDAwAR6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABibQAwDAwAR6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABibQAwDAwAR6AAAYmEAPAAADE+gBAGBgAj0AAAxsw7ILAFi4q7bPvu7Z582vDgA4AgR6gIMR/gE4ygn0wNHtUAI1AByD9KEHAICBCfQAADAwgR4AAAYm0AMAwMBcFAvAUuzYtWfVdbZt2bSASgDGNtdAX1XXJfl8kq8nuaO7t1bVqUnelOSMJNcl+cnuvnWedQAAwHq1iC4339/dW7p763T6xUku6+6zklw2nQYAAA7DMrrcbEvyxOnzi5K8O8mLllAHwMLN0s0EAA7FvFvoO8lfV9XlVXXBdN59uvuGJJk+3nt/G1bVBVW1s6p27t27d85lAgDAmObdQv+E7v50Vd07yTuq6iOzbtjdFya5MEm2bt3a8yoQAABGNtcW+u7+9PTxpiRvS/K4JDdW1elJMn28aZ41AADAeja3QF9VJ1bVSXc+T/KUJB9MckmS86ernZ9kx7xqAACA9W6eXW7uk+RtVXXncd7Y3X9VVe9P8uaqem6STyb5iTnWAAAA69rcAn13X5vkkfuZf3OSc+Z1XAAAOJYsYhx6AABgTgR6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABibQAwDAwAR6AAAYmEAPAAADE+gBAGBgAj0AAAxsw7ILAFg3rtq+6iqbdt+SPZvPXUAxABwrtNADAMDABHoAABiYQA8AAAMT6AEAYGACPQAADMwoNwBH0M7rb1l2CQAcY7TQAwDAwAR6AAAYmEAPAAADE+gBAGBgAj0AAAzMKDcAC7Zp96Uzr7tn87lzrASA9UALPQAADEygBwCAgQn0AAAwMIEeAAAGJtADAMDABHoAABiYQA8AAAMT6AEAYGACPQAADEygBwCAgQn0AAAwMIEeAAAGJtADAMDABHoAABiYQA8AAAMT6AEAYGACPQAADEygBwCAgQn0AAAwMIEeAAAGNlOgr6ofrirhHwAAjjKzhvRnJrmmqn67qh4yz4IAAIDZzRTou/tnkzwqyceTvK6q/q6qLqiqk+ZaHQAAcFAzd6Pp7s8l+ZMkFyc5PcmPJrmiqn5xTrUBAACrmLUP/dOr6m1J3pnk+CSP6+6nJXlkkl+ZY30AAMBBbJhxvfOSvKq7/3blzO7+UlX9n0e+LAAAYBazdrm5Yd8wX1WvSJLuvuyIVwUAAMxk1hb6Jyd50T7znrafeQAcQZt2Xzrzuns2nzvHSgA4Wh000FfVv0jyC0m+s6o+sGLRSUneM8/CAACA1a3WQv/GJH+Z5LeSvHjF/M939y1zqwoAAJjJaoG+u/u6qnrevguq6lShHgAAlmuWFvofTnJ5kk5SK5Z1ku+YU10AAMAMDhrou/uHp49nLqYcAADgUMx6Y6knVNWJ0+c/W1WvrKr7z7jtcVV1ZVX9+XT6zKp6b1VdU1Vvqqq7HX75AABwbJt1HPrXJPlSVT0yyQuTXJ/k9TNu+0tJPrxi+hWZ3KTqrCS3JnnujPsBAAD2MWugv6O7O8m2JK/u7ldnMnTlQVXV5iTnJvnD6XQleVKS7dNVLkryjEMtGgAAmJg10H++ql6S5GeTXFpVxyU5fobtfjeTFv1vTKdPS3Jbd98xnd6dZNP+NqyqC6pqZ1Xt3Lt374xlAgDAsWXWQP9TSb6a5Lnd/ZlMQvj/e7ANquqHk9zU3ZevnL2fVXt/23f3hd29tbu3bty4ccYyAQDg2LLasJVJkmmIf+WK6U8m+c+rbPaEJE+vqh9KckKSe2bSYn9yVW2YttJvTvLpwykcAACYfZSbH5uOSvMPVfW5qvp8VX3uYNt090u6e3N3n5HkmUne2d3PSvKuJOdNVzs/yY411A8AAMe0Wbvc/HaSp3f3vbr7nt19Unff8zCP+aIkv1xVH8ukT/1rD3M/AABwzJupy02SG7v7w6uvtn/d/e4k754+vzbJ4w53XwAAwD+aNdDvrKo3JfnTTC6OTZJ091vnUhUAADCTWQP9PZN8KclTVszrJAI9AAAs0ayj3Dxn3oUAAACHbqZAX1UPTPKaJPfp7odX1SMyuUj2N+ZaHbA+XbV99XUgyY5dew66fNuW/d6bEOCYMusoN/8pyUuS3J4k3f2BTIaiBAAAlmjWQH+P7n7fPvPuONLFAAAAh2bWQP/ZqvrOTC6ETVWdl+SGuVUFAADMZNZRbp6X5MIkD66qPUk+keRZc6sKgEO2afelM623Z/O5c64EgEU6aKCvql9eMfkXSd6VSav+F5P8eJJXzq80AABgNau10J80fXxQkscm2ZGkkjw7yd/OsS4AAGAGBw303f2yJKmqv07y6O7+/HT615K8Ze7VAQAABzXrRbH3T/K1FdNfS3LGEa8GAAA4JLNeFPv6JO+rqrdlMtLNjya5aG5VAQAAM5kp0Hf3b1bVXyb5nums53T3lfMrCwAAmMWsLfTp7iuSXDHHWgAAgEM0ax96AADgKCTQAwDAwAR6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABibQAwDAwAR6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABibQAwDAwAR6AAAY2IZlFwCsI1dtX3YFAHDM0UIPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABibQAwDAwAR6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMLANyy4AAA7Xjl17Drp825ZNC6oEYHm00AMAwMAEegAAGJhADwAAAxPoAQBgYC6KBTgEO6+/ZdklAMA30UIPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMLC5jXJTVSck+dskd58eZ3t3v7SqzkxycZJTk1yR5Nnd/bV51QHAN9u0+9KZ192z+dw5VgLAkTDPFvqvJnlSdz8yyZYkT62qxyd5RZJXdfdZSW5N8tw51gAAAOva3AJ9T3xhOnn89KeTPCnJ9un8i5I8Y141AADAejfXPvRVdVxV7UpyU5J3JPl4ktu6+47pKruTbDrAthdU1c6q2rl37955lgkAAMOaa6Dv7q9395Ykm5M8LslD9rfaAba9sLu3dvfWjRs3zrNMAAAY1kJGuenu25K8O8njk5xcVXdejLs5yacXUQMAAKxHcwv0VbWxqk6ePv/WJD+Q5MNJ3pXkvOlq5yfZMa8aAABgvZvbsJVJTk9yUVUdl8kHhzd3959X1YeSXFxVv5HkyiSvnWMNAACwrs0t0Hf3B5I8aj/zr82kPz0AALBG7hQLAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABibQAwDAwAR6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABibQAwDAwAR6AAAYmEAPAAADE+gBAGBgG5ZdAABHr027L5153T2bz51jJQAciBZ6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDDDVgIHd9X2ZVcAAByEFnoAABiYQA8AAAMT6AEAYGACPQAADEygBwCAgQn0AAAwMIEeAAAGJtADAMDA3FgKgHVrx649q66zbcumBVQCMD9a6AEAYGACPQAADEygBwCAgQn0AAAwMIEeAAAGJtADAMDADFsJMLXz+luWXQIAHDIt9AAAMDCBHgAABibQAwDAwAR6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABibQAwDAwAR6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIFtWHYBwBJctX3ZFQAAR8jcWuir6n5V9a6q+nBVXV1VvzSdf2pVvaOqrpk+njKvGgAAYL2bZ5ebO5L8P939kCSPT/K8qnpokhcnuay7z0py2XQaAAA4DHML9N19Q3dfMX3++SQfTrIpybYkF01XuyjJM+ZVAwAArHcLuSi2qs5I8qgk701yn+6+IZmE/iT3PsA2F1TVzqrauXfv3kWUCQAAw5n7RbFV9W1J/iTJC7r7c1U103bdfWGSC5Nk69atPb8KATgSNu2+dKb19mw+d86VABxb5tpCX1XHZxLm39Ddb53OvrGqTp8uPz3JTfOsAQAA1rN5jnJTSV6b5MPd/coViy5Jcv70+flJdsyrBgAAWO/m2eXmCUmeneSqqto1nferSV6e5M1V9dwkn0zyE3OsAQAA1rW5Bfru/h9JDtRh/px5HRcAAI4lCxnlBgAAmA+BHgAABibQAwDAwAR6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABibQAwDAwAR6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABibQAwDAwAR6AAAYmEAPAAAD27DsAgAWZef1tyy7BAA44rTQAwDAwAR6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIEZthKAhdq0+9KZ192z+dw5VgKwPmihBwCAgQn0AAAwMIEeAAAGJtADAMDABHoAABiYQA8AAAMT6AEAYGACPQAADMyNpQA4pu3Yteegy7dt2bSgSgAOjxZ6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDCBHgAABmYcegA4COPUA0c7LfQAADAwgR4AAAYm0AMAwMAEegAAGJhADwAAAxPoAQBgYIathPXkqu3LrgAAWDAt9AAAMDCBHgAABibQAwDAwAR6AAAYmEAPAAADM8oNAEetTbsvnXndPZvPnWMlAEevubXQV9UfVdVNVfXBFfNOrap3VNU108dT5nV8AAA4Fsyzy80fJ3nqPvNenOSy7j4ryWXTaQAA4DDNLdB3998muWWf2duSXDR9flGSZ8zr+AAAcCxY9EWx9+nuG5Jk+njvA61YVRdU1c6q2rl3796FFQgAACM5ake56e4Lu3trd2/duHHjsssBAICj0qID/Y1VdXqSTB9vWvDxAQBgXVl0oL8kyfnT5+cn2bHg4wMAwLoyz2Er/2uSv0vyoKraXVXPTfLyJE+uqmuSPHk6DQAAHKa53Viqu3/6AIvOmdcxAQDgWHPUXhQLAACsbm4t9ACwSJt2XzrTens2nzvnSgAWSws9AAAMTKAHAICBCfQAADAwfegBYA127Nqz6jrbtmxaQCXAsUoLPQAADEygBwCAgQn0AAAwMIEeAAAGJtADAMDABHoAABiYQA8AAAMzDj0c7a7avuwKAICjmBZ6AAAYmEAPAAADE+gBAGBgAj0AAAxMoAcAgIEJ9AAAMDDDVgLrxs7rb1l2CQCwcFroAQBgYAI9AAAMTKAHAICBCfQAADAwF8UCwJzt2LXnoMu3bdm0oEqA9UigB4ZgBBsA2D9dbgAAYGACPQAADEygBwCAgelDDwBL5qJZYC200AMAwMC00MMyXLV92RUAAOuEQA/AMWXT7ktnXnfP5nPnWAnAkaHLDQAADEygBwCAgQn0AAAwMIEeAAAGJtADAMDAjHIDAAdgRBxgBAI9cFTYef0tyy4Bjlqr3Uk2cTdZOJbpcgMAAAMT6AEAYGACPQAADEwfegA4Ama9gNbFs8CRpoUeAAAGJtADAMDABHoAABiYQA8AAANzUSwcKVdtX3YFAMAxSAs9AAAMTKAHAICB6XIDLMTO629ZdglwVJh1vPrEmPXAbLTQAwDAwLTQA8A6sGPXnoMu37Zl04IqARZNCz0AAAxMCz0AHKX0twdmIdADAKt22Ul024GjlUAPq3EcccXwAAAHb0lEQVTDqFUZwQYAlmcpgb6qnprk1UmOS/KH3f3yZdTBMUxIB45VB/j7t2n3XT+Y68YDY1h4oK+q45L8QZInJ9md5P1VdUl3f2jRtQCz0QIPR7/V+tvv3D3/GmYeaWfWRpWzz1tjRXBsWEYL/eOSfKy7r02Sqro4ybYkAj0AHEX2/ZBwsA8Fs7Tm3xn49/dtwH73+fW7fkDQjx/uahmBflOST62Y3p3kn+67UlVdkOSC6eRXq+qDC6iN8f2TJJ9ddhEMw/nCrJwrHArnC7N60JHYyTICfe1nXt9lRveFSS5Mkqra2d1b510Y43OucCicL8zKucKhcL4wq6raeST2s4wbS+1Ocr8V05uTfHoJdQAAwPCWEejfn+Ssqjqzqu6W5JlJLllCHQAAMLyFd7np7juq6vlJ3p7JsJV/1N1Xr7LZhfOvjHXCucKhcL4wK+cKh8L5wqyOyLlS3Xfpvg4AAAxiGV1uAACAI0SgBwCAgS010FfVU6vqo1X1sap68X6W372q3jRd/t6qOmPFspdM53+0qn5wkXWzHId7vlTVk6vq8qq6avr4pEXXzmKt5W/LdPn9q+oLVfUri6qZ5Vnj/6JHVNXfVdXV078xJyyydhZrDf+Hjq+qi6bnyIer6iWLrp3Fm+F8+d6quqKq7qiq8/ZZdn5VXTP9OX+1Yy0t0FfVcUn+IMnTkjw0yU9X1UP3We25SW7t7u9K8qokr5hu+9BMRsd5WJKnJvkP0/2xTq3lfMnk5h4/0t1nJzk/yesXUzXLsMZz5U6vSvKX866V5Vvj/6INSf5Lkn/e3Q9L8sQkty+odBZsjX9bfiLJ3af/hx6T5P/etyGB9WXG8+WTSX4+yRv32fbUJC/N5Marj0vy0qo65WDHW2YL/eOSfKy7r+3uryW5OMm2fdbZluSi6fPtSc6pqprOv7i7v9rdn0jysen+WL8O+3zp7iu7+857HVyd5ISquvtCqmYZ1vK3JVX1jCTXZnKusP6t5Xx5SpIPdPffJ0l339zdX19Q3SzeWs6VTnLi9EPgtyb5WpLPLaZslmTV86W7r+vuDyT5xj7b/mCSd3T3Ld19a5J3ZNKAfUDLDPSbknxqxfTu6bz9rtPddyT5hySnzbgt68tazpeVfjzJld391TnVyfId9rlSVScmeVGSly2gTo4Oa/nb8sAkXVVvn35t/sIF1MvyrOVc2Z7ki0luyKRV9ne6+5Z5F8xSrSWrHvK2Cx+HfoXaz7x9x9A80DqzbMv6spbzZbKw6mGZfP35lCNYF0eftZwrL0vyqu7+wrTBnvVvLefLhiT/LMljk3wpyWVVdXl3X3ZkS+QosZZz5XFJvp7kvklOSfLfq+q/dfe1R7ZEjiJryaqHvO0yW+h3J7nfiunNST59oHWmX1PdK8ktM27L+rKW8yVVtTnJ25L8XHd/fO7VskxrOVf+aZLfrqrrkrwgya9Ob4TH+rXW/0V/092f7e4vJfmLJI+ee8Usy1rOlZ9J8lfdfXt335TkPUm2zr1ilmktWfWQt11moH9/krOq6syqulsmF7less86l2RyEWOSnJfknT25E9YlSZ45vZr8zCRnJXnfgupmOQ77fKmqk5NcmuQl3f2ehVXMshz2udLd39PdZ3T3GUl+N8m/6+7fX1ThLMVa/he9Pckjquoe0/D2fUk+tKC6Wby1nCufTPKkmjgxyeOTfGRBdbMcs5wvB/L2JE+pqlOmF8M+ZTrvgJbW5aa775i2fL09yXFJ/qi7r66qX0+ys7svSfLaJK+vqo9l8gn3mdNtr66qN2fyh/OOJM9zIdL6tpbzJcnzk3xXkn9TVf9mOu8p01YS1pk1niscY9b4v+jWqnplJv+4O8lfdPelS3khzN0a/7b8QZLXJflgJt0pXje9GJJ1apbzpaoem0nvgVOS/EhVvay7H9bdt1TVv83kb0uS/Ppq11zU5IMjAAAwIneKBQCAgQn0AAAwMIEeAAAGJtADAMDABHoAABiYQA9AkqSqvrDsGgA4dAI9AAAMTKAHWKeq6hVV9Qsrpn+tql5aVZdV1RVVdVVVbdvPdk+sqj9fMf37VfXz0+ePqaq/qarLq+rtVXX6Ql4MAAck0AOsXxcn+akV0z+Zyd0qf7S7H53k+5P8+6qqWXZWVccn+f+SnNfdj0nyR0l+88iWDMCh2rDsAgCYj+6+sqruXVX3TbIxya1Jbkjyqqr63iTfSLIpyX2SfGaGXT4oycOTvGP6GeC46f4AWCKBHmB9257kvCTfnkmL/bMyCfeP6e7bq+q6JCfss80d+eZvcO9cXkmu7u7vnmvFABwSXW4A1reLkzwzk1C/Pcm9ktw0DfPfn+QB+9nm+iQPraq7V9W9kpwznf/RJBur6ruTSRecqnrY3F8BAAelhR5gHevuq6vqpCR7uvuGqnpDkj+rqp1JdiX5yH62+VRVvTnJB5Jck+TK6fyvVdV5SX5vGvQ3JPndJFcv6OUAsB/V3cuuAQAAOEy63AAAwMAEegAAGJhADwAAAxPoAQBgYAI9AAAMTKAHAICBCfQAADCw/x9J41jDwKjc3gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 900x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(mcmc.trace('p_a')[:], bins=35,\n",
    "        histtype='stepfilled', normed=True,alpha=.35)\n",
    "plt.hist(mcmc.trace('p_b')[:], bins=35,\n",
    "        histtype='stepfilled', normed=True,alpha=.35)\n",
    "plt.xlim([0, .1])\n",
    "plt.ylim([0, 80])\n",
    "plt.xlabel('value')\n",
    "plt.ylabel('density')\n",
    "plt.legend(loc='upper right')\n",
    "plt.title('Posterior distribution of $p_A$ and $p_B$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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
}

bayesian-gaussian-mixture-model.ipynb

bayesian-least-squares.ipynb

challenger_causal_analysis.ipynb

recieving-messages.ipynb