{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8537f44f",
   "metadata": {},
   "source": [
    "### STK 4021 2021 Oblig"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "322e9d47",
   "metadata": {},
   "source": [
    "#### Problem 4 - Rats"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f19c2f48",
   "metadata": {},
   "source": [
    "#### Import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "22d6d573",
   "metadata": {},
   "outputs": [],
   "source": [
    "import autograd.numpy as np\n",
    "from numpy.linalg import det, inv\n",
    "import matplotlib.pyplot as plt\n",
    "from autograd import grad, hessian\n",
    "from scipy.optimize import minimize\n",
    "import arviz as az\n",
    "from scipy.stats import gaussian_kde as kde\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5896e58c",
   "metadata": {},
   "source": [
    "#### Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "853c9404",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 6\n",
    "x = np.array([0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0])\n",
    "y = np.array([2, 1, 2, 2, 3, 4, 1, 2, 5, 3])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "749a2bd3",
   "metadata": {},
   "source": [
    "#### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "975b3dd5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def H(arg):\n",
    "    return np.exp(arg) / (1 + np.exp(arg))\n",
    "\n",
    "\n",
    "def p(x, a, b):\n",
    "    arg = a + b * x\n",
    "    return H(arg)\n",
    "\n",
    "\n",
    "def ell(a, b, x, y):\n",
    "    # omitting the binomial coefficient term\n",
    "    p_vec = p(x, a, b)\n",
    "    return np.sum(y * np.log(p_vec) + (m - y) * np.log(1 - p_vec))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77f039b9",
   "metadata": {},
   "source": [
    "#### (a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "dc9f4cda",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "  message: Optimization terminated successfully.\n",
       "  success: True\n",
       "   status: 0\n",
       "      fun: 39.46178286453997\n",
       "        x: [-1.173e+00  7.471e-01]\n",
       "      nit: 6\n",
       "      jac: [ 0.000e+00  4.768e-07]\n",
       " hess_inv: [[ 3.636e-01 -2.554e-01]\n",
       "            [-2.554e-01  2.226e-01]]\n",
       "     nfev: 24\n",
       "     njev: 8"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res = minimize(lambda params: -ell(params[0], params[1], x, y), x0=np.zeros(2))\n",
    "res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "cb962d87",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum likelihood estimate of a: -1.1735\n",
      "Maximum likelihood estimate of b: 0.7471\n"
     ]
    }
   ],
   "source": [
    "mle = res.x\n",
    "print(f'Maximum likelihood estimate of a: {np.round(mle[0], 4)}')\n",
    "print(f'Maximum likelihood estimate of b: {np.round(mle[1], 4)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23b9394c",
   "metadata": {},
   "source": [
    "#### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "c037904c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAALEwAACxMBAJqcGAAAJ2ZJREFUeJzt3Xt8VPWd//HXhxEIIBcpICioaAFz0VWJooLcpVyECCRtcdtVm6o/W8Fdq0Ubt93aUrrb/myVtVpb/K1YiTUhJJGbARLAWLCAWqqgFtGKihXlIsgtiZ/fHzOwIeYySCaTzHk/H495ZOackzPvkwPnM+f7nfM95u6IiEhwtYp3ABERiS8VAhGRgFMhEBEJOBUCEZGAUyEQEQk4FQIRkYCLWSEws8fM7EMze6WO+WZmD5rZVjPbZGaXxCqLiIjULZZnBP8DjK1n/jigX+RxM/BwDLOIiEgdYlYI3H0NsKueRTKAeR62DuhiZr1ilUdERGp3Shzf+0xge7XX70am7ai5oJndTPisgQ4dOgw8//zzmySgiEii2Lhx40fu3r22efEsBFFz90eBRwHS09N9w4YNcU4kItKymNnf65oXz28NvQf0qfa6d2SaiIg0oXgWgmLgXyLfHroc2Ovun2sWEhGR2IpZ05CZ5QLDgW5m9i7wI6A1gLs/AiwBxgNbgQPAjbHKIiIidYtZIXD3aQ3Md+C7sXp/ERGJjq4sFhEJOBUCEZGAUyEQEQk4FQIRkYBTIRARCTgVAhGRgFMhEBEJOBUCEZGAUyEQEQk4FQIRkYBTIRARCTgVAhGRgFMhEBEJOBUCEZGAUyEQEQk4FQIRkYBTIRARCTgVAhGRgFMhEBEJOBUCEZGAUyEQEQk4FQIRkYBTIRARCTgVAhGRgFMhEBEJOBUCEZGAUyEQEQk4FQIRkYBTIRCRRpObm0taWhqhUIi0tDRyc3PjHUmicEq8A4hIYsjNzSUnJ4e5c+cyZMgQysvLyc7OBmDatGlxTif1MXePd4YTkp6e7hs2bIh3DBGpIS0tjTlz5jBixIhj08rKypg+fTqvvPJKHJMJgJltdPf0WuepEIhIYwiFQhw6dIjWrVsfm1ZRUUFSUhJVVVVxTCZQfyFQH4GINIrk5GTKy8uPm1ZeXk5ycnKcEkm0VAhEpFHk5OSQnZ1NWVkZFRUVlJWVkZ2dTU5OTryjSQPUWSwijeJoh/D06dPZsmULycnJzJo1Sx3FLUBM+wjMbCzwABACfu/uP68x/yzgcaBLZJm73X1JfetUH4GIyImLSx+BmYWAh4BxQAowzcxSaix2L/C0u18MfB34TazyiIhI7WLZR3AZsNXdt7n7EeApIKPGMg50ijzvDLwfwzwiIlKLWBaCM4Ht1V6/G5lW3X8A3zCzd4ElwPTaVmRmN5vZBjPbsHPnzlhkFREJrHh/a2ga8D/u3hsYDzxhZp/L5O6Punu6u6d37969yUOKiCSyWBaC94A+1V73jkyrLht4GsDd1wJJQLcYZhIRkRpiWQjWA/3MrK+ZtSHcGVxcY5l3gFEAZpZMuBCo7UdEpAnFrBC4eyVwG/AssIXwt4NeNbP7zGxSZLHvATeZ2V+AXOAGb2ljXoiItHAxvaAsck3AkhrTfljt+WZgcCwziIhI/eLdWSwiInGmQiAiEnAqBCIiAadCICIScCoEIiIBp0IgIhJwKgQiIgH3hQqBmZ3e2EFERCQ+oi4EZtbFzLLNbCXwUgwziYhIE6r3ymIza0f4HgLXARcDHYFrgTUxTyYiIk2izjMCM5sPvAFcDcwBzgF2u/sqd/+saeKJiEis1dc0lALsJjxg3BZ3ryJ8RzEREUkgdRYCd78I+Crh5qAVZlYOdFRHsYhIYqm3s9jdX3P3H7n7+cDtwOPAejP7U5OkExGRmIt6GGp33whsNLO7gKtiF0lERJrSCd+PIHLjGH1rSEQkQejKYhGRgGuwEJhZ32imiYhIyxTNGcGCWqblN3YQERGJjzr7CMzsfCAV6GxmU6rN6gQkxTqYiIg0jfrOCAYA1wBdgInVHpcAN8U8mUgzk5ubS1paGqFQiLS0NHJzc+MdSaRR1HlG4O5FQJGZXeHua5swk0izk5ubS05ODnPnzmXIkCGUl5eTnZ0NwLRp0+KcTuTkWPjboPUsYNad8BnAOVQrHO7+rZgmq0N6erpv2LAhHm8tAZaWlsacOXMYMWLEsWllZWVMnz6dV155JY7JJJHt3r2bLl26YGZs3ryZlJSUL7wuM9vo7um1zYums7gI6AysABZXe4gExpYtWxgyZMhx04YMGcKWLVvilEgS1euvv84vf/lLhg0bRvfu3dm0aRMAp512WszeM5oLytq7+8yYJRBpAZKTkykvLz/ujKC8vJzk5OQ4ppJE8uqrrzJ58mT+9re/AXDBBRcwc+bMYwWgV69eMXvvaM4IFpnZ+JglEGkBcnJyyM7OpqysjIqKCsrKysjOziYnJyfe0aQF2r17N7m5uVx33XU88MADAJx99tn069ePOXPm8Pbbb7Np0yZmzZrFWWedFfM80ZwR3A78wMyOAEcAIzzSRKeYJhNpRo52CE+fPp0tW7aQnJzMrFmz1FEsJ+SRRx7h6aefZs2aNVRVVdG9e3fS0tIAOPXUU1m8OD6t7g12Fjc36iwWkZbA3Xn55Zf585//zC233ALAuHHj2L59O5MmTWLixIlcdtllhEKhJslTX2dxg2cEZmbAPwN93f0nZtYH6OXuf27knCIiLVpFRQVr1qyhsLCQ4uJi3nnnHUKhEFlZWXTt2pUFCxbQvn37eMf8nGj6CH4DXEH4vsUA+4GHYpZIRKQF2bdvH59++ikAc+fOZfTo0cydO5eLLrqIxx57jB07dtC1a1eAZlkEILo+gkHufomZvQTg7rvNrE2Mc4mINFs7duyguLiYwsJCSktLefjhh/nWt77F5MmT6dWrF1dffXWzPejXJppCUGFmISL3K45cYKab14tI4Bw8eJARI0bwwgsvAHDeeedx2223cemllwJw+umnk5GREc+IX0g0heBBYCHQw8xmAZnAvTFNJSISZ1VVVaxdu5aioiIqKir49a9/Tbt27ejfvz8TJ07k2muvJSUlhXA3asvWYCFw9yfNbCMwivBXR691d11OKSIJ6fnnn2fevHkUFhby4Ycf0rp1a6655hrcHTNj3rx58Y7Y6KK5Mc2DQFd3f8jd//tEioCZjTWz181sq5ndXccyXzWzzWb2qpnNP4HsIiIn7cCBAxQWFnL48GEAli1bxpNPPsnw4cN56qmn2LlzJwUFBQnxyb8u0Qw6dz3wNcLDUi8EnnL3Br/IH+lXeAO4GngXWA9Mc/fN1ZbpBzwNjIx0Qvdw9w/rW6+uIxCRk7V3714WL15MQUEBS5cu5cCBAyxevJjx48ezZ88ekpKSSEpKrNuunNR1BO7+OPC4mXUFpgL/aWZnuXu/Bn71MmCru2+LhHgKyAA2V1vmJuAhd98dea96i4CIyBd1tGnntdde48ILL6SiooJevXpx4403MmXKFIYOHQpAly5d4hs0DqLpLD7qy8D5wNlANM1DZwLbq71+FxhUY5n+AGb2PBAC/sPdl9VckZndDNwMNMm4GyKSGLZv387ChQspKCjgwgsv5MEHH6R///7cfffdjBs3jkGDBtGqVTSXUyW2aK4s/i9gMvAm8BTwE3ff04jv3w8YDvQG1pjZBTXX7+6PAo9CuGmokd5bRBLUI488wmOPPcb69euB8P0kzjvvPABatWrFfffdF894zU40ZwRvAle4+0cnuO73gD7VXveOTKvuXeAFd68A3jKzNwgXhvUn+F4iElDuzl/+8heeffZZvv/972NmvPjii5gZP//5z5k8eTL9+/ePd8xmLZrO4laEh5c4193vM7OzgJ4NjTVkZqcQ7iweRbgArAeuc/dXqy0zlnAH8vVm1g14CbjI3T+ua73qLBaRyspKysvLj3X4btu2jVatWrF582YGDBhAVVVVkw3m1lKc7B3KHiI81tDR8Xb3EcVYQ+5eCdwGPEu4T+Fpd3/VzO4zs0mRxZ4FPjazzUAZcFd9RUBEgquqqoo1a9Ywffp0evTowYgRI3jggQcYMGAAv/vd79ixYwcDBgwAUBE4QTEda8jdlwBLakz7YbXnDtwReYiIHKeqqory8nLy8vJYsGABH3zwAUlJSYwZM4bk5GTuueceOnfuHO+YLZ7GGhKRZuXoJ/+8vDwKCgr4xz/+Qbt27Rg/fjxZWVlMmDCBU089Nd4xE4rGGhKRuKusrDzu4P/hhx/Svn17JkyYQGZmJhMmTKBDhw7xjpmwNNaQiMRFZWUlq1atIj8/n4KCAnbu3En79u255ppryMrKYty4cTr4N5GoLihz99eA12KcRUQSXGVlJWVlZeTl5bFw4UI++ugjOnTowMSJE8nMzGTcuHEtahz/RHEiVxaLiJywiooKSktLyc/PZ+HChXz88ceceuqpTJw4kaysLMaOHUu7du3iHTPQ6iwEZtbW3Q83ZRgRSQwVFRWsXLmSvLw8CgsL2bVrFx07dmTSpElkZmbyla98RQf/ZqS+M4K1wCVm9oS7f7OpAolIy3TkyBFWrFhBfn4+hYWF7N69m06dOjFp0iSysrIYM2ZMwo3omSjqKwRtzOw64Eozm1JzprsXxC6WiLQER44cYfny5eTl5VFUVMSePXvo1KkTGRkZxw7+bdu2jXdMaUB9heD/AP8MdAEm1pjngAqBSAAdPnyYkpIS8vPzKSoqYu/evXTu3Jlrr72WrKwsRo8erYN/C1NnIXD3cqDczDa4+9wmzCQizcyhQ4coKSkhLy+P4uJiPvnkE7p06cLkyZOPHfzbtIlqwAFphqL51tATZjYDGBp5vRp4JDJiqIgkqEOHDrFs2TLy8/MpLi5m3759nHbaaWRmZpKZmcmoUaN08E8Q0RSC3wCtIz8Bvgk8DHw7VqFEJD4OHjzIsmXLyMvL45lnnmH//v107dqVr371q2RlZTFy5Ehat24d75jSyKIpBJe6+z9Ve11qZn+JVSARaVoHDhw4dvBftGgR+/fv50tf+hLTpk0jMzOTESNG6OCf4KIZhrrKzM47+sLMzgWqYhdJ5IvJzc0lLS2NUChEWloaubm58Y7UbB04cID8/Hy+9rWv0aNHD6ZOncqKFSu47rrrWL58OR988AGPPvooY8aMaXZFQPs5Bty93gfhMYbeAVYR7h94GxjR0O/F6jFw4EAXqWn+/Pnet29fLy0t9SNHjnhpaan37dvX58+fH+9ozcb+/fv96aef9qysLG/fvr0D3r17d7/lllt8xYoVXlFREe+IDdJ+/uKADV7Xcb6uGcctBG2BCyOPttH8TqweKgRSm9TUVC8tLT1uWmlpqaempsYpUfOwb98+f+qpp3zq1Knerl07B7xHjx5+6623emlpqVdWVsY74gnRfv7i6isEDd6qsrnRrSqlNqFQiEOHDh3XjFFRUUFSUhJVVcFqydy3bx+LFy8mLy+PpUuXcvDgQU4//XSmTp1KVlYWV111VYu9g5f28xdX360qNeicJITk5GTKy8sZMWLEsWnl5eUkJyfHMVXT2bdvH8888wx5eXksW7aMQ4cO0atXL7Kzs8nKymLw4MEt9uBfXdD3c8zUdarQXB9qGpLaBLHteO/evf6HP/zBMzIyvG3btg74GWec4dOnT/c1a9Z4VVVVvCM2uiDu58bCyfQREB5KYgLQqqFlm+KhQiB1mT9/vqempnqrVq08NTU1IQ8Oe/bs8Xnz5vnEiRO9TZs2DviZZ57pM2bM8Oeeey4hD/41BWE/x0J9haDBPgIzGw3cCFwO5AH/z91fj+FJSr3URyBBs2fPHoqLi8nLy6OkpIQjR47Qu3dvMjMzycrK4vLLL6dVq2i+CS5BdlJ9BO6+AlhhZp2BaZHn24HfAX9wDTUh0uh2795NUVER+fn5lJSUUFFRQZ8+ffjud79LVlYWgwYN0sFfGk1UncVm9iXgG4SHl3gJeBIYAlwPDI9VOJEg2bVrF0VFReTl5bFixQoqKio466yzmDFjBllZWVx22WWYWbxjSgJqsBCY2UJgAPAEMNHdd0Rm/dHM1EYjchI+/vhjCgsLyc/PZ8WKFVRWVnLOOedw++23k5WVxaWXXqqDv8RcNGcEv3P3JdUnHL2NZV3tTSJSt48++ojCwkLy8vIoLS2lsrKSvn37cscdd5CVlcXAgQN18JcmFU0h+CmwpMa0tcAljR9HJDH94x//OPbJv6ysjKqqKs4991y+973vkZWVxSWXXKKDv8RNfTev7wmcCbQzs4uBo/9KOwHtmyCbSIv2/vvvU1BQwIIFC1izZg2fffYZ/fr1Y+bMmUydOpWLL75YB39pFuo7I/gKcAPQG7i/2vR9wA9imEmkxdq+fTsLFiwgPz+fP/3pT7g7KSkp3HvvvWRmZpKWlqaDvzQ79d2q8nHgcTOb6u4LmjCTSIvy1ltvHTv4v/DCCwBceOGF/PjHP2bq1KmkpKTEOaFI/eprGvqGu/8BOMfM7qg5393vr+XXRALhb3/7G/n5+SxYsICNGzcCMHDgQGbPns3UqVPp169fnBOKRK++pqEOkZ+nNkUQkeZuy5Yt5Ofnk5+fz6ZNmwAYNGgQv/jFL5g6dSp9+/aNc0KRL6a+pqHfRn7+uOniiDQf7s4rr7xy7OC/efNmAAYPHsyvfvUrpkyZwllnnRXnlCInr76moQfr+0V3n9H4cUTiy915+eWXjx3833jjDcyMoUOHMmfOHCZPnsyZZ54Z75gijaq+pqGNTZZCJI7cnfXr1x9r89+2bRuhUIjhw4fzb//2b1x77bX07Nkz3jFFYqahbw2JJKSqqirKy8tZuHAhCxcu5J133uGUU05h9OjR3HPPPWRkZNC9e/d4xxRpEvU1Df3a3f/VzJ4BPjdWtbtPamjlZjYWeAAIAb9395/XsdxUIB+41N01fpHExOHDhyktLaWgoICioiJ27txJ27ZtGTNmDD/+8Y/JyMjgtNNOi3dMkSZXX9PQE5Gfv/wiKzazEPAQcDXwLrDezIrdfXON5ToCtwMvfJH3EanP/v37WbZsGQUFBSxevJhPPvmEjh07MmHCBKZMmcLYsWPp2LFjvGOKxFV9TUMbIz9Xm1kb4HzCZwavu/uRKNZ9GbDV3bcBmNlTQAawucZyPwH+E7jrxOOLfN6uXbt45plnKCgooKSkhEOHDtGtWzeysrKYMmUKo0aNom3btvGOKdJsRDMM9QTgEeBNwuMN9TWzW9x9aQO/eiawvdrrd4FBNdZ9CdDH3RebWZ2FwMxuBm4G9HU9qdX7779PYWEhCxcuPDaoW+/evbn55puZMmUKgwcP5pRTorr9hkjgRPM/4/8CI9x9K4CZnQcsBhoqBPUys1aExzC6oaFl3f1R4FEI36ryZN5XEsebb77JwoULKSgoYO3atQD079+f73//+0yePJn09HSN6yMShWgKwb6jRSBiG+GB5xryHtCn2uvekWlHdQTSgFWR/6w9gWIzm6QOY6nN0Qu8CgoKKCgoOHZ17yWXXMJPf/pTJk+eTHJysg7+Iieovm8NTYk83WBmS4CnCfcRZAHro1j3eqCfmfUlXAC+Dlx3dKa77wW6VXu/VcCdKgJSXVVVFevWraOoqIiFCxeydetWzIwhQ4Zw//33M3nyZM4555x4xxRp0eo7I5hY7fk/gGGR5zuBdg2t2N0rzew24FnCXx99zN1fNbP7gA3uXvwFM0uCO3jwICtWrKCoqIhnnnmGDz/8kNatWzNy5EjuuusuMjIyOP300+MdUyRhmHvLanJPT0/3DRt00pBoPvroIxYtWkRRURElJSUcOHCATp06MX78eDIyMhg3bhydO3eOd0yRFsvMNtZ1e+FovjWUBGQDqUDS0enu/q1GSyiB9Oabb1JUVERhYSHPP/88n332Gb179+bGG28kIyODYcOG0aZNm3jHFEl40XQWPwG8RviOZfcB/wxsiWUoSUyfffYZGzZsoKioiKKiIl599VUgfBOXnJwcMjIydO9ekTiIphB82d2zzCzD3R83s/nAc7EOJonh8OHDlJWVUVRURHFxMe+//z6hUIihQ4dy0003MWnSJI3jLxJn0RSCisjPPWaWBnwA9IhdJGnpdu/ezZIlSygqKmLZsmXs27ePDh06MHbsWDIyMpgwYQJdu3aNd0wRiYimEDxqZqcB/w4UE75j2b/HNJW0OK+//jqLFi1i0aJFPPfcc1RVVdGzZ0+mTZtGRkYGI0eOJCkpqeEViUiTa9XQAu7+e3ff7e6r3f1cd+9x9O5l8r9yc3NJS0sjFAqRlpZGbm5uvCOdkBPNf+TIEUpLS7njjjvo378/559/PnfeeSe7du1i5syZrF27lvfee4/f/va3jB8/nqSkpBb/N0oE2gdSK3ev9wF8CZgDvEj4ZjW/Br7U0O/F6jFw4EBvbubPn+99+/b10tJSP3LkiJeWlnrfvn19/vz58Y4WlWjz79y50+fNm+dZWVneqVMnB7xt27Y+btw4f+ihh/zvf//7Sb+HxI72QbARvn6r9uN8XTOOLQDLCTcF9Y087gVWNPR7sXo0x0KQmprqpaWlx00rLS311NTUOCU6MXXlT0lJ8b/+9a/+s5/9zK+88ko3Mwe8Z8+e/u1vf9uLiop8//79J/UeLeVvlAi0D4KtvkLQ4AVlZvaKu6fVmPZXd7+g8c5LotccLygLhUIcOnSI1q1bH5tWUVFBUlISVVVVcUwWner5Dx06xKpVqyguLubhhx8+tszAgQOZOHEi11xzDRdffDGtWjXYqljnexzVkv5GiUD7INhO6oIyoMTMvk54rCGATMLDRkhEcnIy5eXljBgx4ti08vJykpOT45gqel/+8pf5wQ9+wNatW1m+fDmffvopbdu2pWPHjtx///2MHz+eM84446Teo6X/jRKB9oHUqa5TBcIjjH4S+fkZUBl5fAZ8UtfvxfrRHJuGWlrba0VFha9Zs8bvuecev+iii5zwYILeo0cPv+WWW3z27Nl+9tlnN2r+lvY3SkTaB8HGyfQRNLdHcywE7uH/ZKmpqd6qVStPTU1tdv+53nvvPZ87d65nZmZ6586dHfBQKORDhw712bNn++zZsz0lJSWm+Zv73ygItA+Cq75CENWgc2Y2CRgaebnK3Rc18olJ1JpjH0FzVFlZydq1a1m6dClLly7l5ZdfBuCMM85g7NixjB8/ntGjR2sgN5GAONlB534OXAo8GZl0u5kNdvd7GjGjNIIdO3awbNkyli5dyvLly9mzZw+hUIjBgwcze/Zsxo0bx4UXXqixfETkONF0Fo8HLnL3zwDM7HHgJUCFIM4qKytZt24dS5YsOe5Tf69evZgyZQrjxo1j9OjRdOnSJa45RaR5i/Zu3l2AXZHnakuIo7feeouSkhJKSkpYuXIle/fuJRQKceWVV/Kzn/2M8ePH61O/iJyQaArBz4CXzKwMMMJ9BXfHNJUc88knn1BWVnbs4L91a/j20X369CErK4uvfOUr+tQvIiel3kJgZq0If130csL9BAAz3f2DWAcLqqqqKjZs2EBJSQnPPvss69ato6qqig4dOjBixAhmzJjBmDFj6N+/vz71i0ijqLcQuPtnZvZ9d3+a8MijEgNvv/02y5cvp6SkhBUrVrBnzx7MjIEDBzJz5kzGjBnDFVdcobt1iUhMRNM0tMLM7gT+CHx6dKK776r7V6Q++/btO9bcs3z5ct544w0AevfuzZQpUxgzZgyjRo2iW7ducU4qIkEQTSH4WuTnd6tNc+Dcxo+TmA4fPsy6detYuXIlK1eu5M9//jOVlZW0b9+e4cOH853vfIcxY8Zw/vnnq7lHRJpcg4XA3XUfwRNUVVXFiy++SGlpKStXrqS8vJyDBw/SqlUr0tPTueuuu7j66qu58soradu2bbzjikjARXNBWRLwHWAI4TOB54BH3P1QjLO1GO7Oli1bjn3iX7VqFXv37gUgNTWVm266iVGjRjF06FB9u0dEmp1omobmER54bk7k9XXAE0BWrEK1BG+//faxT/ylpaV88EH4i1TnnHMOmZmZjBo1ipEjR3L66afHOamISP2iKQRp7p5S7XWZmW2OVaDm6r333mP16tWsWrWKlStXsm3bNgB69OjByJEjGTVqFKNGjaJvX7WkiUjLEk0heNHMLnf3dQBmNghI+FHftm/fzqpVq1i9ejWrV68+diFXp06dGD58ODNmzGDUqFGkpqaqg1dEWrRoCsFA4E9m9k7k9VnA62b2V8Dd/cKYpWtCb7/99rFP/KtXr+att94CoEuXLlx11VXceuutDBs2jIsuuohQKBTntCIijSeaQjA25imamLvz1ltvHTvor1q1infeCde5rl27MnToUGbMmMHw4cO54IILdOAXkYQWzddH/94UQWLtnXfeoaSk5NjB/9133wWgW7duDBs2jDvvvJPhw4eTmpp6wvfjFRFpyaIdfbTFy83N5e6776ZHjx4MHz6cYcOGMWzYMFJSUtTGLyKBFphCcP3115ORkcGAAQN04BcRqSYwhaBnz5707Nkz3jFERJodNYaLiAScCoGISMDFtBCY2Vgze93MtprZ5+5qZmZ3mNlmM9tkZivN7OxY5hERkc+LWSEwsxDwEDAOSAGmmVlKjcVeAtIjF6XlA/8VqzwiIlK7WJ4RXAZsdfdt7n4EeArIqL6Au5e5+4HIy3VA7xjmERGRWsSyEJwJbK/2+t3ItLpkA0trm2FmN5vZBjPbsHPnzkaMKCIizaKz2My+AaQDv6htvrs/6u7p7p7evXv3pg0nIpLgYnkdwXtAn2qve0emHcfMRgM5wDB3PxzDPCIiUotYnhGsB/qZWV8zawN8HSiuvoCZXQz8Fpjk7h/GMIuIiNQhZoXA3SuB24BngS3A0+7+qpndZ2aTIov9AjgVyDOzl82suI7ViYhIjMS0j8Ddl7h7f3c/z91nRab90N2LI89Hu/vp7n5R5DGp/jUGV25uLmlpaYRCIdLS0sjNzY13JIkB7WeJh8CMNdSS5ebmkpOTw9y5cxkyZAjl5eVkZ2cDMG3atDink8ai/Sxx4+4t6jFw4EAPmtTUVC8tLT1uWmlpqaempsYpkcSC9rPEErDB6ziuWnh+y5Genu4bNiT8LZOPEwqFOHToEK1btz42raKigqSkJKqqquKYTBqT9rPEkpltdPf02uY1i+sIpH7JycmUl5cfN628vJzk5OQ4JZJY0H6WeFEhaAFycnLIzs6mrKyMiooKysrKyM7OJicnJ97RpBFpP0u8qLO4BTjaUTh9+nS2bNlCcnIys2bNUgdigtF+lnhRH4GISACoj0BEROqkQiAiEnAqBCIiAadCICIScCoEIiIBp0IgIhJwKgQiIgGnQiAiEnAqBCIiAadCICIScCoEIiIBp0IgIhJwKgQiIgGnQiAiEnAqBCIiAadCICIScCoEIiIBp0IgIhJwKgQiIgGnQiAiEnAqBCIiAadCICIScCoEIiIBp0IgIhJwKgQiIgGnQiAiEnAqBCIiAadCICIScDEtBGY21sxeN7OtZnZ3LfPbmtkfI/NfMLNzYplHREQ+L2aFwMxCwEPAOCAFmGZmKTUWywZ2u/uXgV8B/xmrPCIiUrtYnhFcBmx1923ufgR4CsiosUwG8HjkeT4wyswshplERKSGU2K47jOB7dVevwsMqmsZd680s73Al4CPqi9kZjcDN0de7jez179gpm411x0A2uZg0DYHw8ls89l1zYhlIWg07v4o8OjJrsfMNrh7eiNEajG0zcGgbQ6GWG1zLJuG3gP6VHvdOzKt1mXM7BSgM/BxDDOJiEgNsSwE64F+ZtbXzNoAXweKayxTDFwfeZ4JlLq7xzCTiIjUELOmoUib/23As0AIeMzdXzWz+4AN7l4MzAWeMLOtwC7CxSKWTrp5qQXSNgeDtjkYYrLNpg/gIiLBpiuLRUQCToVARCTgErIQBHFoiyi2+QYz22lmL0ce345HzsZiZo+Z2Ydm9kod883MHoz8PTaZ2SVNnbGxRbHNw81sb7V9/MOmztjYzKyPmZWZ2WYze9XMbq9lmYTZ11Fub+PvZ3dPqAfhjuk3gXOBNsBfgJQay3wHeCTy/OvAH+Oduwm2+Qbgv+OdtRG3eShwCfBKHfPHA0sBAy4HXoh35ibY5uHAonjnbORt7gVcEnneEXijln/bCbOvo9zeRt/PiXhGEMShLaLZ5oTi7msIf9OsLhnAPA9bB3Qxs15Nky42otjmhOPuO9z9xcjzfcAWwiMSVJcw+zrK7W10iVgIahvaouYf8rihLYCjQ1u0VNFsM8DUyKlzvpn1qWV+Ion2b5JorjCzv5jZUjNLjXeYxhRpwr0YeKHGrITc1/VsLzTyfk7EQiC1ewY4x90vBJbzv2dEkjheBM52938C5gCF8Y3TeMzsVGAB8K/u/km888RaA9vb6Ps5EQtBEIe2aHCb3f1jdz8cefl7YGATZYuXaP4dJBR3/8Td90eeLwFam1m3OMc6aWbWmvBB8Ul3L6hlkYTa1w1tbyz2cyIWgiAObdHgNtdoM51EuO0xkRUD/xL5RsnlwF533xHvULFkZj2P9nWZ2WWE/3+35A84RLZnLrDF3e+vY7GE2dfRbG8s9nOLGH30RHjzHNoipqLc5hlmNgmoJLzNN8QtcCMws1zC357oZmbvAj8CWgO4+yPAEsLfJtkKHABujE/SxhPFNmcCt5pZJXAQ+HoL/4ADMBj4JvBXM3s5Mu0HwFmQkPs6mu1t9P2sISZERAIuEZuGRETkBKgQiIgEnAqBiEjAqRCIiAScCoGISMCpEEjgmdl/mNmd8c4hEi8qBCIiAadCIIFkZjlm9oaZlQMDqk2/yMzWRQbnW2hmp0Wmz4iMEb/JzJ6KTLvMzNaa2Utm9iczGxCZ3t7Mno4sv9DC97xIj8wbE/mdF80sLzKmTPVcp5jZejMbHnk928xmNckfRQJLF5RJ4JjZQOB/gEGEr65/kfD9KX5pZpuA6e6+OnJldid3/1czex/o6+6HzayLu+8xs07AgciV3aOBW919aqSZqZ+732JmacDLhMfJfxsoAMa5+6dmNhNo6+731ciXSnh49OnAL4BBkeHFRWIi4YaYEInCVcBCdz8AYGbFkZ+dgS7uvjqy3ONAXuT5JuBJMyvkf0d77Aw8bmb9ACcy3AMwBHgAwN1fiRQXCBeDFOD5yFAxbYC1NcNFhgd5AlgEXKEiILGmQiASnQmE7xA2EcgxswuAnwBl7j45Mnb8qgbWYcByd58WxftdAOwBenzRwCLRUh+BBNEa4Foza2dmHQkf3HH3vcBuM7sqstw3gdVm1gro4+5lwEzCZwKnRn4eHe74hmrrfx74KoCZpRA+qAOsAwab2Zcj8zqYWf+a4cxsCtCVcOGZY2ZdGmOjReqiMwIJHHd/0cz+SPjezh8SHsb7qOuBR8ysPbCN8EiWIeAPkaYjAx6M9BH8F+GmoXuBxdXW8ZvI9M3Aa8CrhIdG3mlmNwC5ZtY2suy9hO9LC0BkXPmfA6PcfbuZ/TfhZqbrEYkRdRaLNDIzCwGt3f2QmZ0HrAAGqK1fmiudEYg0vvZAWeROUwZ8R0VAmjOdEYiIBJw6i0VEAk6FQEQk4FQIREQCToVARCTgVAhERALu/wPN+7CCZkqt+QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mplot = 200\n",
    "space1 = np.linspace(0, 2.0, mplot)\n",
    "space2 = np.linspace(2.0, 2.5, mplot)\n",
    "\n",
    "plt.plot(space1, H(mle[0] + space1 * mle[1]), color='black')\n",
    "plt.plot(space2, H(mle[0] + space2 * mle[1]), color='black', linestyle='dashed')\n",
    "plt.scatter(x, y / m, facecolor='white', edgecolor='black')\n",
    "plt.ylim(0, 1)\n",
    "plt.xlabel('dosage x')\n",
    "plt.ylabel('probability of event A')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35667026",
   "metadata": {},
   "source": [
    "#### (b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "fd5cd204",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = 10_000\n",
    "sample = np.zeros((sim, 2))\n",
    "\n",
    "for s in range(sim):\n",
    "    if s == 0:\n",
    "        sample[s] = 16 * np.random.random(2) - 8\n",
    "    else:\n",
    "        current = sample[s - 1]\n",
    "        proposal = current + (2 * np.random.random(2) - 1)\n",
    "        if np.max(proposal) > 8 or np.min(proposal < -8):\n",
    "            sample[s] = current\n",
    "        else:\n",
    "            log_acceptance_prob = ell(proposal[0], proposal[1], x, y) - ell(current[0], current[1], x, y)\n",
    "            U = np.random.random()\n",
    "            if np.log(U) < log_acceptance_prob:\n",
    "                sample[s] = proposal\n",
    "            else:\n",
    "                sample[s] = current\n",
    "                \n",
    "### burn-in\n",
    "sample = sample[int(.05 * sim):]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "c74a1c65",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAALEwAACxMBAJqcGAAARiZJREFUeJztnXmYFMX5x7+19y4LuyywcspyCYIoyCF4AYLibbwiamLwCGrUJGr0hzEajVExHokmRkM0osYD44UHqIAgKgisHMqxnCI3u1y7XHvX74/pmunp6aO6u7p7eqY+z8PDzkx3VXV11dtvv/XW+xJKKSQSiUQSXjKCboBEIpFI3CEFuUQikYQcKcglEokk5EhBLpFIJCFHCnKJRCIJOVlBVNq2bVtaVlYWRNUSiUQSWr799tvdlNJ22u8DEeRlZWUoLy8PomqJRCIJLYSQH/W+l6YViUQiCTlSkEskEknIkYJcIpFIQo4U5BKJRBJypCCXSCSSkCMFuUQikYQcKcglEokk5EhBLpGkONWHG/DKgk1obGoOuikSjxAiyAkh/yGEVBJCVogoTyKRiOMfc9bh/mkrMX/DnqCbIvEIURr5FABnCypLIpEIZF3lQQBAY7PUyFMVIYKcUjoPwF4RZUkkErHsPlgHAMjLygy4JRKvkDZyiSTF2XuwHgCQmy2ne6ri250lhEwghJQTQsqrqqr8qlYiSXvyc5gmTgJth8Q7fBPklNLJlNLBlNLB7dolRGGUSCQeQQgT4DLReqoi37UkkjSBSjmesohyP3wDwAIAvQkhWwkh14soVyKRuEfq46mPkMQSlNIrRZQjkUi8Q2rkqYs0rUgkKQ4zkVMpyVMWKcglkhSHKMaVZinHUxYpyCWSFCeqkUsrecoiBblEki5IOZ6ySEEukaQJqSLHZRTHRKQgl0jShFRY62xqphj40Ew8Mn110E1JKqQgl0hSHLazMxVs5EcamnCgthGT520MuilJhRTkEkmakAoa+QtfSgGuhxTkEkmKk0o7O99dsi3oJiQlUpCHkBXbqjHl6x+CboYkJDD3w+YUUMkLc4VsRk85ZK+EkPP//hUAYPwp3QJuiSQMpFLwQyIj8eoiNXKJJE1IhcVOiT5SkAfA5j2Hsf9wfdDNkKQJbIt+ClhWpEZugBTkATDyiTm4+b9Lgm6GJM1IBUEu0UcK8gBopsCCjXuCboYkTYjFWgk/JInT1VUfbsAdby3Dim3VvtctBblEkuJE1zqlSu4p323bj3eXbMPTs9f5XrcU5BJJqkNSJ4xtMtvIWf/WNjT5XrcU5BJJ2pC8knxXTS227jscdDNcwfz0SQBPGynIJZIUJ2ZaCbQZppw86XOc+tgcy8iGnVvnAwB6lhb60Sx7KP2bEcBbgxTkkrSDUoo3F23G4fpG0+N2H6wLZOHKK5JYjqNJsUs0WTxt2GJn7/YtPW+TXZhGniE1conEe+Zv2IOJ736PP324yvS4s/82L7qLNszEcnYG2w4erNoY3dSUhNfCbORBmPGlIPeZQ3XmWqDEe9g92H3QfFOW1e9hIww7O4N82Exbtg0X/uMrx4kraNRGLrJVfMhYKz7zeUVl0E1Ie5JfnIklDDZyhtXDhlK+45xw73srcLCuEYfqm1CUb1/HZS2Si50SiY8ksyubSGKJJZKP+sZm3Pf+iuhn3oeNFw+lqNugjbK37D2MBkWDp1EbueiWWSMFuc/kZkW6PDOIuy1JS5J5Q9CMFTvw6jc/Rj9btdDLS7Ab7nf3wTqc9pc50bWWmI1cauQpDxsi7Vvl+VZnUzPF9O934Ei9/xsVrNhVUyuT6fpEEsrxhDbxClEvriUaXIzz+OojDQCAr9bvjmtTEG96UpD7TBCTacGGPfjVa0vwz7nr/a/chO37j+CkR2bj/g9W+lpvMgo0L2GC5YBqob3SowfoNf9ZhAmvlDs+n9trxQbvL92Gd5dstdEGZwNEuh8K5t/zNuLHPYeCboYB/kuR+qaIJv59kvlE7z8c0WhmrtoVSP3pYtxiWXVWbI3c//2H6zH0kdn488fiM9HPW1uFz9zcT14buY159Nupy3DHW8utD4yaVriLVp8We5uQGrl79h+ux8PTV+Oqfy8MuilJQ4Djy5Tg3OHSSyUvaZELAGihCHT2APXag+rzil2oazQ352mVV26vFU9MK3xt0KI9WmrkAmhUHqdHAghcw0OQr/Vqt6g9B+sMj2tupnh1waZAgv+kEzO+34GL//k19h7yx1+dCaiYm5x3dS3ZvA/XTSnHo9MrbJ0XZGAvuxuntN0XjbUirkncpJwgZzjpzPJNe1FT2yC0HW9/uxX97v8kao9MBl1w1qpdGPTnWZivLNJomb5iB+6bthJPzVzraTvSzVat5cmZa7F08378sNtbM2BUgGv620uBs095OK3aUYNmG9KZ1z7txdCxm0lJe1izsuQgFzsF4FQ4VB6oxWXPL8CfPzLftq3l5y8uxIzvdxj+/uAHK3GovgmHFe02GYTX4h/3AgCWb9W3mbOdj/t80hSDImg/cibg/Mpuv2zLfuyqqfXFDZFVseiHvfiTjTll6X7ovEmW2HU/1MLOCq1phRByNiFkDSFkPSFkoogyncK0jz02hdD2/bUAgCWb9/PXRSm+XLcbN78WTNo2ryakX36wQdnuk+FhCgAZyl6CxiaPG6QUv2zLfpzxxNzo117uQFQLwynzN3GfF+SGoJiN3N7xjFgYW1Et4se1ICeEZAJ4FsA5APoCuJIQ0tdtuY5xeIOZ6SM7k79LGmxMQC+3FovGrzYG3RdBpw1je8KafDQMH6pvEtbrG6sO4pg/zMDGqoN44IOV2LI3Fk/c6SVZKSexn73rMzumIDUfLt8OIJhxJSLWylAA6ymlGwGAEPImgIsA2LNRCMLxAGL/WwykusYmHKlvQnFBDhqbOfxwtavyAscfpcGbB9yQLJpxULBXcKvQrW7Rli7qTejdJdtQ39iMhz9ejdkVlVi2Zb+qDmfXtLOmFq1b5NhSqETh9A2FrXF8uW63Uo6wJnEjorc6Adii+rxV+S4OQsgEQkg5IaS8qqpKQLX6ONXyeMfdb95YhlHK6+nhJNwpqYfdOeWlRlHb0ISyiR/jvaVbQ/BuIoaGpmZc9tx8fLZyZ9z3UUHOoxB4gcvbzOZazNsjdkedKlQX/uNr3PPu96a1RupyVr4ZooKLpXSsFUrpZErpYErp4Hbt2nlYj9PzIidW7DyAL9cZP2g+WbkT+xQ/XCeLIm7H3/rKA7bLirqb2azrw++22zzDmt2K2+MTn64NLOyn3XuwYls1np3jfFds9ZEGlP+4D394f0Xca3uGMvv8j1AgRgrGtqQn3kA3ZrO3v7XehWlUemNTM57/YoMz11mXi52MsC52bgPQRfW5s/KdJ7yyYBN+/uJCw1c3tyvOAPDzFxc5KoOrHpeD5H8cg1wUtQ3+SJigTCy88+38v3+Fxz9d47gedn2VB+pw4bOxRBWZAjTyu99ejqmLN1vUr9/B6stncUMYB+saUTbxYzz/xQbjcjXlqGvh1cjtmlCsxsq7S7Zh0owKPOMgk73txU6DARRW08piAL0IId0IITkAxgH4QEC5utw/bSW+XLfbcKAcdJC4obmZ4tsf99lvjB0BJEhYqT0c7D4UuAeYDwOR0qCXOt1DKbVcGGtuppi2LKbXrNhWE/2bCYJGF4udb5Vvxf+9k2iKOFTXiBnf7wClFJ+ujN8yr9Wkv99ajRMe/CyunWzD2OsLjR8SP1RFbMN6ey+0Y/N/5VsSjgGcRwE1GvssfZ8TOcD6w61GHsp45JTSRgC3AvgUwGoAb1FKPY+CZHQjGxrt34TJX250pHG5ud3q9j8yfTWe+ixS/95D9abCQT1Ewi8Ig26BO6bM34Tuv59uuDOztqEJ3X8/XTemycrt1dHFQS+8Vv4xZz1ufm0Jlm3Zn7DLWatJr9oR2U/wtcEGMSM+UWz+32zcm/CbVhje9fZ3tso2wqqn3AhRuzs7ed50GE3NFEs3O1AWORFiI6eUTqeUHkMp7UEpfVhEmZZ1Gn5vf1Ks3XXA+iBB6N37yfM24pnP12Pb/iM48aGZeHi6cTAjvXFKKcWcisqE12ORHKxrxPrKg8LKi0y4gGzkguTmW+URM9f2/Ud0f1/0Q6KAY9ypCuLkhSBnSaNrao010yMNTXhq5lrU67jRuu0jr9dvvdnZqZTNefFGt03vvj89ex0u/ud8z4R5qHZ2qgPwGPW1owHoeIHU+hitjIqt9CdKL7aT8pMVOxN+i5anOo/Vv3J7Da6dshiPqDS/5maKmat2RUwYAiTX+P8swpinvnBdDiPSLmHFGTJ/w258sVZ/8drtA8Rqsdbo8qoO1KFiZ0x58NKPXK9prN+37juCZ2avw2tKYgc9byU7faS+n15dEetzL8eOjd0hut+u01F4KnZETGq7aoxjHLkhVIL8jqkxLYYJxFXba+Kecn6+rbux8joVrnrnsQdcherNYmr5FvzylXJMXbyFu5Vrdx0wbFe5kzUEC/y4V1f9eyF+8Z/4xWvR1nkjd02jvtSaMLwQ5GbDS3v9zMODhZF47JMKnPfMl67qd2tnDiK0cTQlHrdphb9sr8d6qAS5+uayTjz3mS9x8T/nq7637rLp3++I+vTOXVMpLFJibUMTPv5uB8omfpyw2LJ0yz7c/N9vcbvqYaSFNX3b/iN4dcEm3WNYKAEgNiGZHFCvG+2ojhy3kzO2xtLN+3DWX+fhxa9+iPvey+w9YbeRM63arkauPd7NYidj2COz40yEWh/vuHZpqtuujJUPl2/H/A278dzcDThU7y42kNu3wFXba3S/p5r/RZIQV9wCJ23wyowYKkGu1iQMTSsc5fzqtSWY8Oq3WF95EONfWowZOqaMBg4Bpm5DTW0D+tz3CW55PRJ3Zeu+w3HH3vX2d7r1xLc9VuB902Lrxc3NFM/N3YADtQ34WCdAF1scjRsjVK+vjEfRZmV79fKt1aZuaW6JNw2xdjsf3X94/3tcP2Wx63a5we3kFBE0a2dNbdyGo9juTYKTe7SJO1ZraqpvjI117cJllkOvErfPJlFvTVUH+E0ZvIuddY1N2H+4Hrtqas0P9BERW/R9I94Gp9/bdubEAZOQtf/39nd46ooB5u1R/T3PwA7LsDOgtHxeUYnHPqnAxqp421ssfos5trwkPXYL1NPU3NT432/MfajN8DwmBudlMZfSusYmPD93I24a2R25WZnWxWv6Mm5+RF0ME+fEpBkmMcI1B7cpzLFsR/RU1QW7fThZrYHxavz2FBE+98OfvbAQizd554HihFBp5GqM+1qMGProO+PQtLE2xOq69fWlusfYcYcyuqY6RWMy8o2NRV1T1aWzKJpQn0Ff1ak0NC99Yr18YMxdU2n4ViXapGNoI+e8woqdETPClK834a+z1iaYt3ihcX/H3tLcPCjbFubqfn/WX80Xvt1r5ObUNTTH7XJmMPMSW3ewk7SDd6g7EeJemxHDK8gFlGEqpDQ/7TtUH/X1BiJPepYyS49DdU1c5hk1vDbVhPOir9E6XyJxIlsNqq17Dxv+5tb2yfp8e3Vt3Ou/SBZu3IPxLy3GE585342ph3rDjBpDGzlnVzE3RraTtpYzho+2fD2NHDoauQjW7jJ3RXXtKWVw/o7qiKvnok17MeapeQmbkdgGo65tCgAAD9hI7G3XRs5DbUMTvlq3my/AngtCZVpRY9TZXj35HvhwJaYti8UeOeHBz0yPv/S5+RhzbKmtOpx7skT+5xUok+dtjByvEqDqh1qXkgLTutwo6dmZqjcFj3RyFov+x93GDyQn/ObNZbhoQEI8uGhyYy1Gt/OFL801bvVpO6tr0aZQPxqgtnhDc6NpbeY4dT90Ggo2WpbOd7NX70p4gOyqrkWrvOyENjDTlFUrNu85jNYtstFSpwzGim3VKG2Zi9JWebzNj/Lfb370JMm1llBp5Or+rTpQhw1ViVqBV28wdQ7ijsxaLTa5rZEGFssVaGBaSfzKFPUCV6IfvECUwnbW1Ea3VgspVufB9tF32+MW9SIHiKnPKEiSUV99v00/M5O2mKZmimGPzsYjBhvETG3krEwQf31yFVybVnTOX64Kk8swSt7MvrcK5XL643Nw2XML4s7RVn3+37/C6Y/PsWhxPDsVTyBtghuvDJWhEuRqRj/5BUY/mWin80ojzxDcU6yZ6ldDQ9MKZ1nx+/fVXiv8nUJN2sFLbUNTwsJsXAWJf2L5Fn3h5gQ917tbX1+Kj5Rojn7JNTsPpxXbYt5C7IHD7Lwvfb1J9xzT61Avdrra78Dv5745LrGE+F7WL9EgcJXyP08kwjWKXZ2YLHbaDSA37NHZAIDcLH9EbGgFuRGi0p9pb79TO67VqrmB7NUlwdatfP5m4x4Axi6TujZ0h+14+OPVlh44d/5vOc548gs89ol5BvV4m25ipZPnbcDgP8907AKpvWfadjvVju6YugxlEz+2PM5OONOJ734XTf22dV/EDmwlgBPe0OL+ptE22JkSbP8BY8aKnbhIFbHRjAMm4QDswvvwSdDIldPum7YSTc2U+x5QShPcD1fvqMEHy92Fcs6RgjwRHiHtmbbl0TvREdXC1v7DzpIdPzc3Emp0k0E2dp5JoScEgEg8DvUD4j9f/4AZK8w9er5SMqWwdjGamyleVm10smrXM7PXY/dBPn/dfYfqLbdvi3LAeXdp/KLnml0HdDdw2RmLK7bVoK3i6leQo9h3dQp459ut0TgqCf2nOoEp0cSmZWV7dWLcGHXERl68cD/U+87sllYdqOOOrvi3WesSYq1cN2Uxfv2GvjcaL35teguVIOdBrRUYCTYG72Bbt+sAPuZwR3TC5f9aEP37+pfLbZ2rbb6eFw4Bwby15lHtHv+0Im7Aqss9edLnuPal+A03dpIFq/M4frZqF56dox/fmpW4s7o2eg5bGOVJEjDwoZl4VYkZomtqQuS6mputQ8/a5ZbXluC+aSsTXN1yMu09OcwWrce/tAjrdh3Anf9bjvP//lXc8dHz48pi6yb23lJFeRB57X5odNyciti6FCGxmO9WPD17XWyLvvKd9u3ECdqx5pU7b6gEOc/Nvevt2BZ4K/cz3vFtFsTKLT/usfasYPc+Ia609jjV398q8WcoKN5bap7nQytcteV+ZTO8qZor//1N9G+tzVjd/3WNTWhWFvdO+0tkYSkvO6KZ8saW/mJNlVJuTIhpmfBqOX47dRln6/lg7dO6mNnNFBO3QIn4/pm7pgp/eH+F+fl6GqtNjVyUnHGrkT83dwMqD8QLUqM3uNmrd2Hl9shbynat8HVwPSIf9F7nY2WESpDzoPbt3mdpqjDuZC8Wa7yGTcJV22vw9fo9Cb9z2esoXIfaVGuAzN4b+T6hqijXTSnHr9+Mf41lglAdXwYA5q/fjeunLMY2TfjYnKwMUErxXxbNTyOVHp1REedJ5OQWPz3LJPOMoCET854w9koBEgNLTVeZvOK1czHt2n+4Pu4Ny4x1Fn7mPLw8f1PcZyNzy/Uvl+O8ZxLt+ATOLKJmoX/tIvrtz4iUE+RqrAawOpSoFi9Di9rHOLqenkvWzpqYgFO/KudnW2/7XrZlP95d4lmmvji0r/za3bQlLSI247zs+GH6/LyNmF1RmXDtOVkZ2F5dG915ZzWJ9eLWWLF8637D30SPmMQHX/wXt2nstxurYqZEdXwdexq5fq/VNjRhzFPzom9LVoiYPzzJza02/Cw0iQmvhbU5y8AktqumFo9/ar6An1Cm1MgTsdsnVsff+57xq2ozBSa+E8lq4sWtsHMtZq+7P3thoeq4yIENBjZsnnyfWi3XCbyXdv8080lo9CptFJGxRW6W5xqQmb2ZZ83CvGzNZ83vew7yLYbvqD6iCtpGbQ22OoP1iCP1TdHE2Tzwxtoya1qeRvHQ85ZSm/2048LqQa29ly1yzRWdxz9dY7jGo8e2/UcSEmtLP3IHuN05+OZi/TyDfmN085dpNFJ2nN3QAE54ZcEmlE382Fm2coXNJq/pz83dEOctMaeiMiGipJbubVvECQYeOWpVph1ce2po/N+1gmajxeI9Y/ijn2O3IvQptaeIGGmwIh6PB2obMHdN/CY5s63rWh/sJRYmvzcWxQdQqzkSbyKpa2zCh8u3W3o3GV3stn3GSs5Pn1+Q8N0pkz7H5Hn8gt8NoRLkd43tbet4UW81Xrwd2VHWjDQ7bZ5RdtjqHTEBWP4j/6ulHZhroXbnWlOzTclhgNYH/dopi3HO3yLJDtj90HvQqh/ePF38zrcxM1JzM8UrCzYlJKLgxWzx2R7x3hNuoBAzfjft4XuImHHbG0sx/qXF3OFftYvFORYRIQ9oFsW10+apmWtx2xtLMVeJVKo1e8S8Vux32KJN+vPMLwttaGOt8KDXh3Y3DEXCuoq/G7ZMKwbfH9IOXOVI9S7JL9c59zgxg/nnqk0Zc9ZU4lf/XSIsUQeD9ZV2olqFDuYxbTCN8IUvN7qOiaEdW3Y9QMyCYDklopG7L6jG5qYsvRpZzlfecBfa7rN609QKfu35lUqaNZZS0Sxek575zu/csnYIlUZuG5379JnNFFJBOq/c8PJifF6xy3QAqYVby7zIc9mNuyAvTJCrs9tMXbRFuBC3Q2Oz/Tygf/98PQAICWykrnv+ht2YovG64CV6vx2Mvfkb4u+9qNyodhcvtXUW5GRaBnfToj3OSpBb+YwzQc/GbJ0m9g47+8nP1qLnvTP4GmkTmSEIDjQcnZlgdxvxY59WJNxwv5i1uhL/98733Nc9oEux7mDv3rZF9G9RIQyYIFdPcMst5UJqNq5HmzDByn/eWd0mvyk/PjN7Ha7690Lbb0Pasp1o0lf9e2Hc59kVlVhpkDbNDry+/Aye9YIV26qx1sRzTIvVngvtQr22BSyAFnuLNFrfWbVDv7/mb0h06U0WUtu0IkBy/OuLjVxue15RfbiBe7cdIfoTSB0CVZS2zrQfJsjrG5sNAwt1v+djDC4rwU8HdxFSt9l9DdJplAnep2audXT+D4r7YGyruPs2sZDFfmNm1mSKCduhaoTaNLaPI0GE9g1IOxeiygeleG/pVtP8uWEjVBq5XfQGk5M3G6+DwptBQQ0brT9ZDMpQUMfCdmM71Wrk17+8OCEXJKOZAots+PO6wSosg5e49Raqb4rXEM21f38fWXar04a0IFDtXHVgX3CSoDrRHBTbMfvxd4m7tc2apd3LkGykuEaucd+qOhiNFJgMjOSIcUyN5XjCTCfQj3Sn/s5I2NpFK8iNzAhXDj06wS3MLnamcJD3941FWxJjnttAHegK8F9YmyFigw9PFE7Dcx0oHVqN3I19WlQMGq9s5CktyLVj7wyd+OVBQUGxiSPOSjOlhhqM3tDWM60YusvSeFdFO2RFFzvNBVdxQbbp704w1VSF1xaP2TysOdLAtenKiFiMGOvZ7reMd+sjTwhR5ZYFZnE4HcQNe0EePGYFep6M20OS+31Bg92OTh59JpFdNXy75JpNNHL9jPTxZGUQU81OHdTKCVaxyUUIHDuhREVpsXWN+gthZqXzhkw1gj0TWZjgIB9YousjACqVsUJAcMMrNiN9OqjTbCis3sG/yArw71QNilAJctsImtRGW979ws7rmFZzys7MMJwEdQ1NpgmkzchXYmZbCS9P8nKaFPnygh+FVHHb6/wPD4bb1+ap5ZENTrsP1mPow7PwvonXjf828mDngJPqExKJqO6P3VAUXoWfFUVKC/Jk1sjtYPQmojXNEAJQPUuHQUc8o/hQO6Egh9Mqp6pblDDYf8TYg8GNjVqN3f0GEcRN9soDdabZafwe2yJ3KHJn/1H1pxOFwMjM5eRSRN3Z77e6dwXVI7UFeapIchvopYPjjdGhpV/HVgLaIx4nHgz+ILZdpm6WIbORqyWhk6JE3HI3wli7q9gpG3e7D++rR6gEuYgNQWHEyHrBUoKp0Q54u0lj1Tx39SDd79VVWE1KtRZeZSN6nnmhYooRjWgvVTPh6ffYFquR80HihL/IBuiXleTWE1NCJcjtkjIaucEA69w6P/4wItbnPdfAd5ZSGm2Slaam/vmDZe4S2bIQBIFu+jHdjJS6GrleTlKnOBHKqTKXtREdReGqVELI5YSQlYSQZkLIYFGNEkWq3Hw9G/k5x7XXvb7rpixO/NIhWRxL9VZdrP5dvcPUCQdqG1FT2yB84c1ODHMzP3zR482sOJEZ63lY6zLjT5wnIWc/pcr8VdMqT7w7LuBeI18B4BIA8wS0xRV6W3hFhN5MBvRe+YxyMTrJeG5EFofgtZpsL34V20kqIn3enW8tF66RazPtOEW06d7sgfWvL/yJc+0F3IJcdadFCPVYCr3g8KpuV4KcUrqaUmqe4dgnBj40E5v3HEaXkpi5oX2rvABb5B2dW+cruzi9HZJGr4GUAvXKdnQ7bRDR3C17DwvV1DoU5TlK+aaHn6YVv1KI6WGUoYkXJ/0ksm/LlVSAarjy2SYxvrWeEDKBEFJOCCmvqnK2TdzqRX/TnkNx7mep8mamvW5KI19qr0/03DbKAk9Bo1vyb35tCfZbJrlm54lB5KTuUCTwYS/ctGKy2Bng4B7051m2z1EnNHZiWhGikSszSS+xRd8OrXzxFffqvlk6AxNCZgFor/PTvZTSabwVUUonA5gMAIMHD/bkcuatrYrbMRn0JgZRaAdYdLFRc3miL9doXGvreZVzE46o+5Gst1V0s5L1OqttJpnQ4uSyRLqcti7ISfhu2Zb9rnfm8uCVt5GlIKeUjvGkZg/QxkxO0nlgG61ApYgId6+vzzA0gOaz327dIgWcSAEhYg1AjekW/WSV8hw4afuvXvtWWP1G3lgiAoMFRbgNQxryNX7VTibWsR3cb4IRjVagNisauXZCiH7aGwbr0lTDa6/dyxFTmgeRQkyk8BXutWJSYIhlDvcopZRGQwO79ZqJL1dYUUlTt1v3w4sJIVsBDAfwMSHkUzHN0mevhS02T5MAwkmn5Sdh3OHLNBm6KdX3WvFrgCbsHuWsuNIiwFYQ7D4g5uECiH8D3FBl7HUV5s1uvOP0m4170eveGcJi2TO9RPSbUzLgKowtpfQ9AO8JaoslnYvzTX/X+gOr75dVlL4wwSIiasejmxCqehhHXdS2x7+JUbHzgNAFyp2cGd158LMfwi2L+BrPslmJjjEf6q4zIPnUTxOsVpXNbFxDHra/0p68UMVG7u2Q5F3EF6nV8pC0QszHdiVrF/Bg9/6Jut9sOGuzF/lJvUvXTSNCJciN3OEYWlttmBeEzKAGGrlojG3k8RVr1ya8Jllfjf00d7y+0F3WpSBZkERZuvzm9YWbud117RAqQW6lIWo1cifTame1uFdtr2hWJHlQ8kxUzG+nJKcYdx80K8xBm+xw/7SVto43SvIRVvY5zAFgRqgEuZWbZ4IgdzDjt4dAkEf2AwU367XrDX5ryEmqkDsKWNa2MDf6d5sWuSZHpi/PhzgcgR5eJHIOlSC3spHv1oRJDfPKvhkXHN8x4rWSJBLNf1NHcly3llmrK22fM6BLkepTcl5X0IhytUyWLD9eKGGhEuRWNvJ8jfvhrpo6/LD7kK3odsnOyN7t8Ifzj43YyINujMKPHEmkRZIkzy/hpNAwTUqSRfHxQsEMmSA3/12ve0Y9MRevLQzWpiuSsjYtkJuVqWjkQbcmAou74hdJctnCSdZF3FQhWXrXi9scMkFuLsmNNJr7bC6uJDOHlDAEBN67HyYryaJZiSaV3hyTkWQZNl48sEMlyK1MXOmg0bAF3WTSyP0mlS5bdIQ/vzm+cxGeuXJg0M3g4t0l+hvmWhd4k+zBCKmRW0nyEE4EpxgllmCM7lPqqvwWPvuG2yGMAs8IQgjeunE4AG8VES/jbV94QkfT3/0WlEYcqtd3Y7SUKyEgpQR5epkaiKlAe/bqE4XU4kNkT9ukkmnlkYuPiybw8PKqsjMISlokhm+V+O+/L00rVgp56sxvS4heQHIVedmZQhK9zv3dqOjfvz+3j+vyRJBKt7m0VZ6rYE592rfkOs4r17tUmHN+uyVK04q0kUfxeos+G9xHFcU2qUw4vYd3FdqAJyl0mGB+xU7uJ69ZQC/ssV8k+6z0ezh50R+hEuR+PDmTPXcfGwRWNnJRBLmD1IiepYVBN0EobFjXNdrfGZqRAbx143C896uTLSpx0DAOUsC8LNRGfv/5fS2PSXvTitX1f7pyl+s6Jl3S33UZfuBH8uVkZbFO8tx0hYBgaLcSlLVpYXEcvyJ0as+2AloWX3cyMnXCMGyadJ5QQd61TYHlMWlvWvGDS07sbPuc7Ez/hioT3mYa+V1je0ePcUsqaFzJjps+ZmYBqzIyMgi3Jphu91zk9fZox/O2mOYaebLy8nVDfa/TzEaul3hB7QLWtpDfeyHN5nQguDFfZdgw8DY2SRu5HiI18g7F1klPvNj3FTJBnpxDQptizg8IIYbZzJmAVwuI5342CH+7YgA729vGSWyR4WIWMiFk9TAgSC9nADsU5dv3c591xwjH9UnTiobxJ5fhnZuHB90MVKrShW2adJ6rskQoB3rjJDcrA0O7lQAAzj++g3U7ou2RQt9r3Gjk7Eyrh4Gd+yjinndvF7PZ85Y28Zw++PLuUdYHCsaJ10qxi01OaR80S0tpq1wM6loSdDPQrmV8HOke7fQXnjZNOi8hQqOW9q3ycPox7Szr3GWSa1JvEZQQgo7F+Vh+/1n44wXWK+vR87iPlDjFlY1ckUIt88wFi50qRNxzJ2WM7dceXUoKcIbLXcl2EaWstMrjS4HsNgGJHqES5Fr55NWb4jFH2XNv66p4DLCdc7PvHGmoma9+6GwM627+8DHTENglz99gnC5Lr1tYmUUF2VLLTiF4tUmvY/O0VITY2f3aK/XZH2MFSliIF38xWFzDTFC78to+V9OXZ/QpxZd3n8H1dpX2GrlfFj6thm2F3XFgtbjCUx6bMHowjVydMcnu67ubQS6xh5su5l2oO/2YdtwCxMk9Zy6LzMSjLkKdBcm03mj9/g46J4ud2r48uqQARZzmFmkj16BN7cZ4c8IwjO13lC9t6N7W3H9XD0tBzjGwHriwn+FvbKC0VL3qOZ0bUnt3B8/rtjv3Q76T7xrbO06AmGnybu64nsLwwIX9MG5IF56TA8FptXeeeUz0b737YBTSQgpyTo7rVGS5QUIEV590NGaqVq95N+iYzb3crAyugWWWuZ61Ijsz+W5vqc23nbBzH8dOPzcSjPchkKk50GwdxsnDm52i/p89xPq0b4mHfnKc7TL9QiuE7z67t/VJFLhtdK/oxyydvSSFufoaujStaK7/htO6AQBO6xW/E83tg93OEzMzg9ge+GbHT7nWvU86czPLVKldyaJYp0LIUDvwhHxw0yU8OwkZ6mFtVuWq7TWO28PG9tpdB5GVGYvqyHOJQYWD0Pb/iUe3tl1Gpg3XF6mRayjIiTzxB3YpDrYhCrz3x+iWn9arLcratnBtzji9V0TbihPkNieJV1PK7qWJCr1q5S1kxPGdiyyPGX9ymeFvXrtuXzPcuO44COIGqNkY22niEcVTDYM9tJubKTIzCPp1bGV+rqBBd/2p3eI+3zLKPNibVrngaQbrSvZmw954eK5hv8H+DzeEWpAzMjVOtF4HlNImXWiRG/l87cnd9A5PwOjhHYvq5671XUoKNOWFVyM/+zjjRV179To7j2fL9Yld7WtwaqyadrRyP5+8/ISE37i9VnzSdtW3l1n2migFIQT/u8l8z4edFpo9YLu0zo/7bGTiiNZLtJ/5WzLYwb2vOlBn+xwrQiXIkyVxxMRzjwUQE7e5WZnYNOk8/GZML+OTVBgJM/ZAMtPi7Gh4LVU71rwwaQwpsz+IjTaudNZMPkbQb1tue02v27UhEqwEBxPkBjU4aofb63pA2YtQo9Eu1eWyN2bmlCByDP5mtPFc05o5LOPQaDVyjmbGdk/bY2i3Ek/85EMlyLu3jWlHY4417gynr9F2aXIYu8Jo4rLgW05ebfXSbeVkequRO9HyjCYzIcCYY917GhmFuLVyKTVKbedWddB78HYoin9oWfWiqHunVoTcmu9KW0ViirBricZUV34fc+xReGn8ENx2Rk90KlaO0VSpHbN22mSk0Nw1tjc6t+ZbN4gKYyd+5MqVmo0PbbmPX3Y83rpxuCeZmkIlyPt2bIVHLmZhZo2FlJcuc7eO6om6hkjuv0P1jY7KMGoe0yRGcOzs1GLloZIklhXDdhAQ3X5R30smEJzA5z3iDLt9q4154ma46p1boPNQIogXfm6nCLuGlhr3SuYRnEGAsrYtcOdZvaP30Inma1W/mmOOKsQto3raviEi3hR4ivByoT9UghwA2iivpV7afM3MF78b2zv6RLXKVDNuSBe0zE30IzY6jQXfuuqko/kaagO7/VXmwD+erx3GGrneL+rvjHyR/3LZ8brHM1rmZmG0hbZ/m8Gr+r3nHYsF95xheq5Z3+r9ph1fvG82ug86neNuOK27zrmE22uFBxoV2PElNZuYUbTfOVlkjNajM0nvOLO3bjl233i4FjsNdpmbyQ43wdGscFU0IeRxQkgFIeQ7Qsh7hJBiQe2y1w6f62Nb8vt0MF+Fn3Tp8Xhx/JCE77UT97FL+yvfK7/blLonddPf8h8/qPjK/PDWU3HjiO6Ycm1iu3k5x2SBkl3aE5qFu3FDnD+8Tu7RJqF8NTyhXm8akejZ0DIvC20LcxNMIXbQm9h2oxDaHQ/Vh+ttHQ8APx1sLw4/u4KocFKayK5N9z4ItNOr9wKe0rMNbjy9O0Yr5lZ1f53V19pc58hGztfM+HI9lFRunxEzARxHKT0ewFoA97hvUvIzqGtrfHTbqZigo/lo0Rcs2mPi7Yt2ue/8vpYLwbyyoH/nItxzzrFow7Gt+n6d4Fvn9e+AX43saXgOc9Pq3ynidXDBCR2x9L4zcdOI7rptHFIWe0g1cQhAvcnidSYluxN098F4QWt1b8x+1hPy3XU8bbRHaf3P7fpORxOcaEo20tQB6weSneeV+mFYnJ+De849NmpeVBfzr58PsixLp6WW52jHFGu7dh6qd357aUVwJcgppZ9RSpmh+BsA9tPrCMbIZ/SPF/TFJSd2in5u7SIMJRDZPcqj6emaC1R3dMWDY6PH8Ghqekfo7SrjaYcbWuVl6S4s8j5QKChWPDgWf/3pCWjdIgeEEPyw+1DC8UerBE7vo/gyxie2KTj0Jm9etr1p10vpZ56HK6BvuosEzYr1xITTe+AfVw3Escpbpbqdo/uUWvp8Gy0UNplo5Fr6dSqKCwerfiicZxFqmTc5AyH6ay/aY+xilINA+wD7/Hcjo9fi5dqdSKvNdQBmGP1ICJlACCknhJRXVVU5roQNoByDxb0TOhfr1w/EzeiZd4zA9F+f5rgdblDf7MLcrNgN5nBp0g1RyyGmRS+0fPfAWMPfzKoa1Tvy+lvSIgeFuVnR3X9AZDegGV1Nwi5MnTAMf7qon37dHktys+tV365nrzoRQGKMIKtbc8Np3fH8zwbhNJ1cmrqnKgWqc28SxNvIC3OzcP7xHdHQ1Bz9nZGZQTD1xuGmma+Y0hFLbBGBar43Y0CXYiy978zYOpLqFNZXRozoZRZiwLLqOJzsMTAKBJadmYGOmgxducru3hwPU0JaRvQhhMwCoGf0vJdSOk055l4AjQBeMyqHUjoZwGQAGDx4sOOpNbxHGwztVhJZnVbgiTqmfRq2LcxNuBndPFjg47EVZsTLcSF1JNYpZhDN/d3I6M40noU8LXef3QfXntINpS2tU2JpKdRZOGZ1ntS9DU7q3gavL9yc+Lvtmsxp3yrPkYso87tvaNIKcvN7U1yQbbgxymx8RcJH6N8Tdh5LiqJ+G8zKJCjMzUIxR+Yc7VhWe63ose7hcwAAh+uaovP2mSsH4u+fr9N1DFBTmJuFg3URA0BRQTb+cdVA3Pr60oTjtIqNkaLD3h7P6d8en62KJW43ux13nHlMnNskS6uoXkc5rVc7TC3fEq31/vP7okNRHs7sK2Zzmx6WgpxSOsbsd0LIeADnAxhNfUjrXpSfjbdujN8hdvVJXXH/tJUAYpNW2xSrp+7c341EieIRozUPLLp3NIY+PNthi60tcGzgOE3FZSWiR/cpNdxw071tC2zUMWkYofZmMfMJNyIzg6C9Tk5RHo4qchZwS/SwLMjV7uzlSyjA3j4abWQWuHxQZ9upBNl90a516wv0yLHdOO6rmphpReO1EjWt6JfB7NhFBbE3sVF9SjGKY5OMUaucbBS88ISOGKwkpfnJgE64fepyy3oA4Nca76bLBnVG28JcjOxt/IZQXJCDu8bqR0IUhVuvlbMB3A3gQkrpYTFNsk9mBsGZVqvTGvcrLWVtW6CVQZYVJ9qjGYluV8rEc7FBQQ92vTeN7GFoz1cvet1zjvFg09tdqGvF8PBRru03vY0mTIBcNsj+cg2v7ZpVV5ibhUmX9MfpvRJNHnq0LsjGFYO7JARGM7vdVvkk9bTN6Bue6mYk7rWIP0YteK3cagHjePWix7ATeOp+5sqB0YBm6mu/5MROtmzZhBCM6lPqqf2bB7c28n8AaAlgJiFkGSHkeQFtEoK2Y+10s0g3IdaM+EiE+u5OPDJQb4eiG68HtX/1jToueEDET3vyNdar/7H6vBnUakF+19je2HMoErOirjGm4RLN/wC/aWXRvWMwVOUlYyTQWDvyczIxbujRupNYf8GR4LHLjscATdgB0TJArz3ab7J04hMxtLGL9NDayK2+F4JBkYmmFM3vNprid5o5Ubj1WulJKe1CKR2g/LtJVMOct0n/ezsDS2RMF1aremIbCQgrE8A5x7WPWxuIq8VhkwdzxEv56eAu6NM+0YtBT2B4GQ9H3W23jIpt/W5uTtQ81U2jFnZbRqu8bOSpdkX+39n6byjMRHKCSeAmO+PNVfJlXRu5SmkwOC8hHkncb9b1xtwP44kJcusyrLhc81bFbVqxqFsb60ZLYa4/IT5EErqdnU6x81Q+zWRF3CnqycVsrEwQaP3IjSb2+JPLdFfLXW3xdiFEjGIwe/WWqX1wRMOk6tl+da4ri0NCMQE15dohGDdUf5NS28IcPHn5CXg4Gi4iEVuC3ORQJ4/FmGlFXYfRgl+sDcw910gjv0C1yKf1F2flM/O/CI183ND4nbw87r6A8aI4Y87vRpr+3rO0JSb/fBCuGd6Vqz4tQQT3SzlBbjR+CPgXvW4e0QPf3DNaUHsiDdIb2BcO6KT8pnxh0Tyjn91MGbeak3bnHKXiBble7BAgZt9X27Z1NXKl57I5LlYrmIy4dFBnHNXKZO3ERh+Ifu7pNd26DmK5KUhdBguVfILGTGS2s9M+zsyjZqGHx/Y7Ci111sJO6dkm7vNZ/dpbPhCs8NNs7q6lIcKOhpDhwrNCC6vVLIMI0x6tvFaMfnaz0OJ2kUZ3B58A0TSgSzGuUGKrvHr9SdhQmehj/tQVAzB//R5d/3I900rrFjk4VH/EtN5Jl/bHs3PW45QebQyP4ekzUXNYW86y+8/EgD/NjH7WGxMxr5VIHHDzuMjsHMOfEuhVWoiTe7TBZ7efHhWaUV3EwJvFCYkLtARHlxREQzIYjTOzuWYUXK44P9HcEvD6pS1STpAbjtmAboqehsigGu3F6oUhJ8vAlAH9SWf2BjK6TylmV1QK1zqozndOeP+WU6J/D+raGoN0Avi3ykv0r45N7lgjWDTJN345DG+Vb8HfP19vWG+Honz8+SfGJhMjltx3Jk58KCZgbZkWbBxaXGAdApUJUXUCA0KAvh1aYdWO+DRuVHOOFR/ddio6t84HIQTHqHbadlTWK9gag9GGPTskLFoC+Oz206NlG5kwnJh1kiXXgVNSxrRyvBK7o4OBJq0Wdn+94gTdY7zELO1abEOQ8WvplUOPxoAu+q++VgNX7+fnfjYI5X8Yg0KOLO92yw4S7YNzaFkJnrlyIICIOeCGU2PxcZyECzZCG2PanhwXvdgZ+V9tYiAgeHrcAJM2mJfJOK5Tke7D5OaRPfDXK07AIxcfh1N6tsHPHdqX49qk492Vl51paStPTCwR+2xHXLt9s/R+V02MlBHkt4zqiY9uOzXBZscIKukvGwzm9esv2qlPufqkow1fGTsZbPYxIycrA20Lc12/qOgFTWLf9CwtxIJ7zkCf9s5ipIigS0lB/IYaVXPNtqCbwdNnfix2XjO8q+5Gr6hpJW6xE+hlEqvG7fQoys/GxQM7o3PrArx2wzCuFHlW6Bnt4j/pN9rJuo+ZP76IsrwmZUwrGRkEx3WKuYNpzQpGu+nenDDM025nE8RsYhuZVoxMJloyiPPdi24nsKkQohQdivIxpKwEFTsPcJd5gyZ5rq32KP9r33Kiv/s0x0RZVg6bJC+5ZVRPU59xCppQ9qw7RmDb/tg6gVEUQyCSfi87k2B4j7aYt9Z5fCSnGG1iYhiZQ7R9csxR1g8V3bKS7XXThJTRyK34+Pud0b/Vg3aYEqPDa9QmQ6ZBs3RZMSEf75vrx24xt9pDYhup5QS0ootpjkpzGjWJDTZWxYcf8GtqOq1HGxe9ptZYkBvVYdbfPUsLdU1KhCQqDW0Lc7H6T2fjkoGdEo73A6uNPryc1qsdfmuRT1ek22QQ9vaU0cit8CEMjCnqQXntyWXoWlIQC4SvfK81rWQQoMnjduVkZaB72xY4q5+zgD4JYpzGvqUGx1jh5l6xyIIsDnR9Y/ybmBfeFHrw+jxHyosde/uZvXDrGT0xc9VO3D51edxmJ8aw7iX4ZuPe6BZzI3hSuxlttWfw+N57Rb7G7VTbxvaKInS8QcRTNVFTj8HQYuF01Yu0TkfKjSN6oHzTPsdzyglpI8gJIb4uPsTqjfyvtm9nZBCMUflfs3ji2ZlaDSRmXLESHuzS1FvMeS43M4Pgc4sNEmboySsWtpPZxv2MQ8E0cubjHdTj2+qKLx/UOeofrz42NysTuVlAfnZ8Bno1z4wbiDW7Dph4sMQ/SCPf6Lcolg0++cwIPUsL0bl1Prbui5iCjtTHqzWDurbGG78chhO7FluWZTUEbxvdC7lZGThVFTfH6bDt0a7Q1ZxyQtoIcjVOb9AQju3sTjitVzvcOKJ7zKMi0YPOFLWg1M/36d0k1Qppioht9ZGL+6tSb3lWfQIlinA7pn1L3DyyB64YHL87UNuUywd1xqzVu2DFL0/rhn9/+QN3O6xe0R9XpbrTO5Q9+PUEeWmrvKhZTg91ebx9H3cO3ym+cMEJHfHc3A0AgAcu7Bf3GyEEw038/e3QqTgfD150XHz5SdUT5qSVIHejnVU8dDZXVDgtbDCYTajMDIJ7zjk24fuMuMllXjcLQ6oOqsXiSWs1fZHo7iIkJO6B4ueE+OOFffHfbzaja0mBbqwUbXsf1+QONeKec45F346tcPvU5VxXo67npWuHYPp3O4yP1QsnwAS5i9dIO6cSEtlb8GAGwdUnxbsOmrnFeg27hrvG9sYlJ7pPQGbHfi0iXoxfpJUgdwILI2o3JjTDTeZsOwsvt47qiSFlJThFlRXmyZ8OwPtLt0XzY3pBovth4kTxUwD061iERy8x3tDj9KGSkUEsQ8qq6VJSgN0H61GYm4VRvUujmZEMGpVAb8Us9ZMB/AuN155ShkFdW6uKo6q/9AVYVEiDoGNxPtY/ci53fX4SIgeSQEhZQa4XP9su79w8XDcugx0yHYzAqBav+q5HqXn2oqzMjDghDkQ2qVznwpWPB57LS6Y56JdAuGlED+yqqTVNOGBGx+J8bJp0nq1z/nhBxPQwU8l2Qynwwi+G4D9f/WC50zJZBWXYd1z6Rcq6H/50SBf88+pY3j8n43RQ15K4bchOiEY2tDEeiwuy0aEoD3++OGazy81KztCa2tdPvcu048GRzNi5h0X52bhmeBmX14doIaoubsQx7fDydUONox/SxHOSiT0H6wEIcJN1cH6Yxm3KCvLcrEyc2988E7cfmAXwMSI7MwML7hmN8/p3tD44YHgmSHimAx88Xjj24pF7A89zx26sFb+p2BmJDVPMkZeXh4C9kD0jZQW5HkH4kjsxrTDCoBAk7LbT62Kb1+HlXXIjr+wuHvIfK/ZGOykuWU1krM/drvMk6XNKGClrIzfDT+0jtv3evngKKj6MHfhCuibPdYhoC08Jdh7CnmnkPGMuyU0rzP3S7lT459UnYsted2mERURw9IvwtNQl6i3I1Yfrfa03iHP9Qi+MrdUxRvxCQMQ8K9z0KYtR369jYto7nZqcV+SSqPJg65zkHGzsWWTXRHlu/w66OWhtrXMo5py+HXjud7CkjUZOKZCrPGGrDvopyJ1PkGSdXGp45lf1kQauspL9eo/rVIQPbz0VfTkEuR25k2ux1d4u7K3DlinI5Lcg7crMjz6It1NWY/d2LRLiuCcb6SPIARytpAYLYmSm6BoLlx/5QZPAT37jVhz0N0m4HFePDcEjPJ6JDY08yM0+PLA0fk6cBkSR7AoGkEaCHAjGVstqTNXVcjbG87MzcaRBP8RXMl26X5MyyIVqO1VrkygnG8//bBDe+XZbNAiaU9xcXXL2TDxpYyMf6lGcFCuSdH4IgwkAq0h8ovjDecfiUhdbtf26Hcm0wGsGC+2QrOO0c+sC/GZML2EPYCcbjJK1b9SkvEbeo10LbKg6hLH92mO6Kia536T6DrVepYVYsb0a44YkBu2yOw/M3l5uOK278Y8ceDkpOxTlYUd1ref18MLjtfL0uIH4Ym0VOhXbzzIVJpy5ZCbusE5WUl6Q6+GnSG2tROO7fcwxPtbqH2yCdCzOx1s3DtfdDZdMjzAvTSvTbjkFQx+ZDSBYU4Wda+zbsRXX4q3dcpMN1vaCHPsiLwzXnfKCXH0TWiqJht1mjrdDXnam7ZgZYYIJLELEbWkOwbzRpbRVHo7t0Aqrd9SgOcBFEdHrMqf2aouczAxce0qZmAIDYFi3NrhsUGfcNIL/jS5MwzDlBTmDAtHM3iIyfEsiEM3/IgjzwjBbK2ho0s8R6wcxP3IxHVnaMg9rHz5HSFlBUVSQjSc4QxZrCYNAT3lBrr4J2ZkZnkcDTDeYFm72+hl0mj0/YXsVtCnm/CQsC60ScaSN10oayZJAmLum0nUZzPSVm+39sBQVhEnLzaN6oH2rPPQotc7c7jVyzIshDN2Y8hq5xFtYLIx9h413b/IuFo0/uQz7DzfgIhvJFJzw9LgBGNjFG3fUUb1L8c3vR3tSNi9hXWNINsLUj1KQS1zBM9Z5TSttCnPx0E+Osz7QJV4/KIKGLea7jaWf7rB+bG2Y5Dp5SHnTSpieqmFE9m/ycXznIjx5+Qm4/czUdHn1izP7HoVHL+mPu8/uDUB8TByRuNLICSEPAbgIQDOASgDjKaXbRTRMNKm+IUciYRBCcOkg94mK0x1CCK4cGtngNn/iGch3mLfXD9w+Yh6nlB5PKR0A4CMA97tvkljkCr638Ni/5SNUEnY6FuejdYvkNbG4EuSUUnVsxxZI4jkrV/AlEkmq4nqxkxDyMIBrAFQDGGVy3AQAEwDg6KMT43F4RVZmRGNk3hUSiUSSalhq5ISQWYSQFTr/LgIASum9lNIuAF4DcKtROZTSyZTSwZTSwe3atRN3BRawnXZ1AW7QSGWk4UoiCR5LjZxSOoazrNcATAfwR1ctEsyjl/THpBkVOK5T8qdrSlXKWEIPiUTiCa5s5ISQXqqPFwGocNcc8fRp3wpTrh2K3KzkXXFOdX6bopEfJZJkwa2NfBIhpDci7oc/ArjJfZMkoUKxrQztVmJ4SHaIspFLJGHElSCnlF4qqiGScDOyt/m6x8Rz+iRV7k6JJJWQW/QlruD1079pRA+PWyKRpC/ynTfEFOV7E8FPIpGEC6mRh5iPbjsVK7ZVB90MAHLDlUQSJFKQh5guJQXoUhKsa58MmiWGl8YPQV1jU9DNkIQUKcglkiRgVJ/SoJsgCTHSRi6RSCQhRwpyiSukZUUiCR4pyCVCSKcEyxJJsiEFucQVcrFTIgkeKcglEokk5EhBLhGCtKxIJMEhBbnEFTKVnkQSPFKQS4QgFXKJJDikIJe4olFJoberpjbglkgk6YsU5BJXlLbMBQAcqZfbyyWSoJCCXOKKNoU5AID6JpkTVSIJCinIJa5g2X8am6SVXCIJCinIJa7Iyoh4rTQ2S41cIgkKKcglrsjJYkNIuiFKJEEhw9hKXHFqz7a47pRuuHFE96CbIpGkLVKQS1yRlZmB+y/oG3QzJJK0RppWJBKJJORIQS6RSCQhRwpyiUQiCTlSkEskEknIkYJcIpFIQo70Wklyplw7BIfqZBwTiURijBTkSc7I3qVBN0EikSQ50rQikUgkIUcKcolEIgk5UpBLJBJJyBEiyAkhdxJCKCGkrYjyJBKJRMKPa0FOCOkC4CwAm903RyKRSCR2EaGR/xXA3ZD5dyUSiSQQXAlyQshFALZRSpdzHDuBEFJOCCmvqqpyU61EIpFIVFj6kRNCZgFor/PTvQB+j4hZxRJK6WQAkwFg8ODBUnuXSCQSQRBKnclUQkh/ALMBHFa+6gxgO4ChlNKdFudWAfjRUcVAWwC7HZ6bCsjrl9efztcPpHcfdKWUttN+6ViQJxREyCYAgymlnnYwIaScUjrYyzqSGXn98vrT+foB2Qd6SD9yiUQiCTnCYq1QSstElSWRSCQSfsKokU8OugEBI68/vUn36wdkHyQgzEYukUgkkmAIo0YukUgkEhVSkEskEknICZUgJ4ScTQhZQwhZTwiZGHR7REAI6UIImUMIWUUIWUkI+Y3yfQkhZCYhZJ3yf2vle0IIeUbpg+8IISeqyvqFcvw6QsgvgromJxBCMgkhSwkhHymfuxFCFirXOZUQkqN8n6t8Xq/8XqYq4x7l+zWEkLEBXYojCCHFhJC3CSEVhJDVhJDh6TQGCCG3K+N/BSHkDUJIXrqNAVdQSkPxD0AmgA0AugPIAbAcQN+g2yXgujoAOFH5uyWAtQD6AvgLgInK9xMBPKb8fS6AGQAIgGEAFirflwDYqPzfWvm7ddDXZ6Mf7gDwOoCPlM9vARin/P08gJuVv38F4Hnl73EApip/91XGRC6AbspYyQz6umxc/8sAblD+zgFQnC5jAEAnAD8AyFfd+/HpNgbc/AuTRj4UwHpK6UZKaT2ANwFcFHCbXEMp3UEpXaL8fQDAakQG9kWITG4o//9E+fsiAK/QCN8AKCaEdAAwFsBMSuleSuk+ADMBnO3flTiHENIZwHkAXlA+EwBnAHhbOUR7/axf3gYwWjn+IgBvUkrrKKU/AFiPyJhJegghRQBOB/AiAFBK6yml+5FGYwARV+h8QkgWgAIAO5BGY8AtYRLknQBsUX3eqnyXMiiviAMBLARwFKV0h/LTTgBHKX8b9UOY++dviETQbFY+twGwn1LaqHxWX0v0OpXfq5Xjw3z93QBUAXhJMS+9QAhpgTQZA5TSbQCeQCQU9g5E7um3SK8x4IowCfKUhhBSCOAdAL+llNaof6OR98aU9BMlhJwPoJJS+m3QbQmQLAAnAniOUjoQwCFETClRUnwMtEZEm+4GoCOAFgjPm0RSECZBvg1AF9Xnzsp3oYcQko2IEH+NUvqu8vUu5XUZyv+VyvdG/RDW/jkFwIVKrJ43EXmdfhoRcwHbeay+luh1Kr8XAdiD8F4/ENEct1JKFyqf30ZEsKfLGBgD4AdKaRWltAHAu4iMi3QaA64IkyBfDKCXspKdg8gixwcBt8k1im3vRQCrKaVPqX76AADzOvgFgGmq769RPBeGAahWXr8/BXAWIaS1ouGcpXyX1FBK76GUdqaREA/jAHxOKb0awBwAlymHaa+f9ctlyvFU+X6c4tHQDUAvAIt8ugxX0Ei00C2EkN7KV6MBrEKajAFETCrDCCEFynxg1582Y8A1Qa+22vmHyGr9WkRWo+8Nuj2CrulURF6ZvwOwTPl3LiI2v9kA1gGYBaBEOZ4AeFbpg+8RiTjJyroOkQWe9QCuDfraHPTFSMS8VrojMgnXA/gfgFzl+zzl83rl9+6q8+9V+mUNgHOCvh6b1z4AQLkyDt5HxOskbcYAgAcBVABYAeBVRDxP0moMuPknt+hLJBJJyAmTaUUikUgkOkhBLpFIJCFHCnKJRCIJOVKQSyQSSciRglwikUhCjhTkEolEEnKkIJdIJJKQ8/9MDFs8o+A4EAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.arange(len(sample[:, 0])), sample[:, 0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d896d444",
   "metadata": {},
   "source": [
    "#### Convergence diagnostics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "409d6ff0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ESS for a: 376.35\n",
      "ESS for b: 403.93\n"
     ]
    }
   ],
   "source": [
    "ess1 = az.ess(sample[:, 0])\n",
    "ess2 = az.ess(sample[:, 1])\n",
    "\n",
    "print(f'ESS for a: {np.round(ess1, 2)}')\n",
    "print(f'ESS for b: {np.round(ess2, 2)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "898029c8",
   "metadata": {},
   "source": [
    "#### Posterior means and standard deviations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "add64f42",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posterior mean of a: -1.2233\n",
      "Posterior mean of b: 0.7768\n",
      "\n",
      "Posterior standard deviation of a: 0.6286\n",
      "Posterior standard deviation of b: 0.4921\n"
     ]
    }
   ],
   "source": [
    "post_mean_a = np.mean(sample[:, 0])\n",
    "post_mean_b = np.mean(sample[:, 1])\n",
    "\n",
    "post_sd_a = np.std(sample[:, 0])\n",
    "post_sd_b = np.std(sample[:, 1])\n",
    "\n",
    "print(f'Posterior mean of a: {np.round(post_mean_a, 4)}')\n",
    "print(f'Posterior mean of b: {np.round(post_mean_b, 4)}')\n",
    "print('')\n",
    "print(f'Posterior standard deviation of a: {np.round(post_sd_a, 4)}')\n",
    "print(f'Posterior standard deviation of b: {np.round(post_sd_b, 4)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbc386ca",
   "metadata": {},
   "source": [
    "#### Lazy Bayes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "ac748502",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posterior lazy Bayes standard deviation of a: 0.6038\n",
      "Posterior lazy Bayes standard deviation of b: 0.4746\n"
     ]
    }
   ],
   "source": [
    "hess = hessian(lambda params: ell(params[0], params[1], x, y))(mle)\n",
    "A = -inv(hess)\n",
    "\n",
    "print(f'Posterior lazy Bayes standard deviation of a: {np.round(np.sqrt(A[0][0]), 4)}')\n",
    "print(f'Posterior lazy Bayes standard deviation of b: {np.round(np.sqrt(A[1][1]), 4)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f523d235",
   "metadata": {},
   "source": [
    "We see that lazy Bayes does very well."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a842c44",
   "metadata": {},
   "source": [
    "#### (c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "ff7914c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "space = np.linspace(0, 2.5, mplot).reshape((mplot, 1))\n",
    "a_plus_bx_sample = (sample[:, 0] + space * sample[:, 1]).T\n",
    "lower90 = np.quantile(a_plus_bx_sample, .05, axis=0)\n",
    "upper90 = np.quantile(a_plus_bx_sample, .95, axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f21cc932",
   "metadata": {},
   "source": [
    "#### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "c2619603",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAALEwAACxMBAJqcGAAAOTZJREFUeJzt3Xl4lNX1wPHvSVjCvq+CCrIYEhFlc0EEQSAghCXRYq1oqaJVaOuvrmjditparbu4oFULQRIhQRAMkABG2XFhE0V2QQHZBUISzu+POyEhJclAMnmTzPk8zzyZ+86bmTMZeM+89973XFFVjDHGBK8QrwMwxhjjLUsExhgT5CwRGGNMkLNEYIwxQc4SgTHGBDlLBMYYE+QClghE5B0R2SUiq/N5XETkJRHZICLfiMilgYrFGGNM/gJ5RvAfoF8Bj0cBrX2324HXAxiLMcaYfAQsEajqQmBvAbtEA++rsxioLSJNAhWPMcaY06vg4WufA2zL1d7u27Yz744icjvurIFq1ap1vPDCC0skQGOMKS9WrFixR1UbnO4xLxOB31T1TeBNgE6dOuny5cs9jsgYY8oWEdmS32Nezhr6EWieq93Mt80YY0wJ8jIRTAdu9s0eugw4oKr/0y1kjDEmsALWNSQicUAPoL6IbAceBSoCqOp44BOgP7ABOALcGqhYjDHG5C9giUBVhxfyuAJ3Ber1jTHG+MeuLDbGmCBnicAYY4KcJQJjjAlylgiMMSbIWSIwxpggZ4nAGGOCnCUCY4wJcpYIjDEmyFkiMMaYIGeJwBhjgpwlAmOMCXKWCIwxJshZIjDGmCBnicAYY4KcJQJjjAlylgiMMSbIWSIwxpggZ4nAGGOCnCUCY4wJcpYIjDEmyFkiMMaYIGeJwBhjgpwlAmOMKQuOHQvYU1siMMaY0ir74L9/P9x9d8BexhKBMcaURv36wU03ufu1a8MNNwTspSwRGGNMafD889CrV067f3+49tqcdu77xcwSgTHGeCEtDWJiID3dtWvUgAYNctpjxsCoUSUSiiUCY4wpCXv2wHPPwbZtrn3oEKxYAZs2ufZtt8HkyVC5comHZonAGFNs4uLiiIyMJDQ0lMjISOLi4rwOyVtbtrgbuAHfv/4VUlJcu29f2LgRLrzQs/CyVfA6AGNM+RAXF8fYsWOZMGEC3bp1Iy0tjZEjRwIwfPhwj6MrQZmZUKGCm/ETEQE33wyvvQatWsHWrdC8udsvpPR8DxdV9TqGM9KpUyddvny512EYY/KIjIzk5ZdfpmfPnie3paamMnr0aFavXu1hZCVo5Ej48UeYPdu1p02Diy+Gli29jQsQkRWq2ul0j5WelGSMKdPWrVtHt27dTtnWrVs31q1b51FEJSA11R38s79QX3opXHFFTnvIkFKRBApjicAYUyzCw8NJS0s7ZVtaWhrh4eEeRRQA6emQlASHD7v2li2QnOzOAgDuugv+9jcQ8S7Gs2CJwBhTLMaOHcvIkSNJTU0lIyOD1NRURo4cydixY70OrWhUISPD3V+6FAYPho8/du3f/tYlg2bNPAuvONhgsTGmWGQPCI8ePZp169YRHh7OuHHjyvZA8bFj0LUrDBvmvul36wZz5sDVV7vHK1b0Nr5iEtBEICL9gBeBUOBtVX0mz+PnAu8BtX37PKCqnwQyJmNM4AwfPrxsH/gBZs500zpHj4awMOjRA9q2dY+JQO/enoYXCAHrGhKRUOBVIApoBwwXkXZ5dnsYmKKqlwC/AV4LVDzGGJOv7Iu8wM30eeUVOHHCtV98MaB1fkqDQI4RdAE2qOpGVT0OTAai8+yjQE3f/VrAjgDGY4wx/+vdd+G88+CHH1z7X/+CtWtL1Tz/QAvkOz0HyJVm2e7blttjwE0ish34BBh9uicSkdtFZLmILN+9e3cgYjXGBIsTJ2DKFPj6a9fu2xeefhrq1nXt2rUhNNSz8LzgdcobDvxHVZsB/YEPROR/YlLVN1W1k6p2atCgQYkHaYwpB7K7en79Fe68E954w7WbNoX774c6dbyLzWOBHCz+EWieq93Mty23kUA/AFVdJCJhQH1gVwDjMsYEm8ceg4ULXZ2fGjVc5c82bbyOqtQI5BnBMqC1iLQQkUq4weDpefbZCvQCEJFwIAywvh9jTNGouoN9VpZrn3OOm/mTXeI5PDzoun8KErAzAlXNFJG7gU9xU0PfUdU1IvIEsFxVpwP/B7wlIn/BDRzfomWt+JExpvSZPdst7DJtmrsA7LbbvI6oVLOic8aYsi97ADgszB34MzNh4kS4/nqoUsXr6EoFKzpnjCmfDhxwP0Xcoi8TJrh2hQowYoQlAT9ZIjDGlE2PPeYqex475hLB9OmuIJw5Y1ZryBhTNqSnw4cfQr9+0LAhREW5b/4ZGa5LqEkTryMss+yMwBhTNmze7Lp7Jk927a5d4eGH3XRQUyR2RmCMKb1eeskt+v7EE2765/LlbvEXU6wsERhjShfVnIVdVq92i75kb+vY0dvYyinrGjLGlB4rVkC7dvDtt6792muuLHQZW/GrrDmrRCAijYo7EGNMkDpxAvbudffPPRfq1cuZFlrBOi1Kgt9/ZRGpDQwDbgTCgaYBiskYEyxU3UIvVaq4b/4NGrjSEKZEFZgIRKQKbg2BG4FLgBrAYGBhwCMzxpRfmzfD+ee7Lp+bbnLTP3OPDZgSlW/XkIhMAr4DrgVeBs4H9qnqfFU9UTLhGWPKnTlz4IILIDnZtX//e7jxRksCHipojKAdsA9YB6xT1SxcYThjjDkzq1bBokXufvfubv6/zQAqNfLtGlLVDiJyIW7xmLkisgeoISKNVPXnEovQGFO2qbrib/Xquf7/ypXh8ce9jsrkUuCsIVX9VlUfVdULgT8B7wHLROSLEonOGFM2ffQRdOniqoCKwKRJrhaQKZX8njWkqiuAFSJyL3BV4EIyxpRZ2QO+Vaq4b/6//AKNGsEll3gdmSnAGV9HoI7NGjLG5Dh6FIYNgxdfdO2oKLc0ZCO75KgssKs1jDFFFxYGISE5M39sBlCZUugZgYi08GebMSbILF8Ol18Ou3e7A/+UKfCnP3kdlTkL/nQNfXSabQnFHYgxpozIzHQ/q1WDfftcUTiws4AyLN+uId/U0QiglogMzfVQTSAs0IEZY0qZEyfcesDnnguvvALh4bBmDYSGeh2ZKaKCzgjaAtcBtYGBuW6XArcFPDJjSpm4uDgiIyMJDQ0lMjKSuLg4r0MqGTt3up8hIRAZCa1b5zxmSaBcKOiCsiQgSUQuV9VFJRiTMaVOXFwcY8eOZcKECXTr1o20tDRGjhwJwPDhwz2OLoD++1+49VZ3ZfCFF8JTT3kdkQkAUS24aoSINMCdAZxPrsShqr8PaGT56NSpky5fvtyLlzZBLDIykpdffpmePXue3Jaamsro0aNZvXq1h5EFwJYtbh3gVq3cQPALL8D//R/Uret1ZEFn37591K5dGxFh7dq1tGvX7qyfS0RWqGqn0z7mRyL4AvgMWAFkZW9X1dMNIgecJQLjhdDQUI4dO0bFihVPbsvIyCAsLIysrKwCfrOMychwYwCXXQbTpnkdTVBav349H3/8MR9//DGff/45K1as4OKLL2bnzp00adLkrJ+3oETgz3UEVVX1/rN+dWPKgfDwcNLS0k45I0hLSyM8PNzDqIrJnj0QHw933gkVK8K770JEhNdRBZ01a9YwZMgQvv/+ewAuuugi7r//furUqQNQpCRQGH+mj84Qkf4Bi8CYMmDs2LGMHDmS1NRUMjIySE1NZeTIkYwdO9br0IouLg7uugvWrnXtfv2geXNvYyrn9u3bR1xcHDfeeCMv+q7GPu+882jdujUvv/wymzdv5ptvvmHcuHGce+65AY/HnzOCPwEPichx4DgguEoTNQMamTGlSPaA8OjRo1m3bh3h4eGMGzeubA4UHz0KzzwDV1wBffvCbbdBr15urWATUOPHj2fKlCksXLiQrKwsGjRoQGRkJADVq1dn5syZnsRV6BhBaWNjBMacpfR0VwjuxAlo08ZdE/Cvf3kdVbmlqnz11VcsXbqUUaNGARAVFcW2bdsYNGgQAwcOpEuXLoSW0BTcog4WC/BboIWqPikizYEmqrq0+EMtnCUCY87CuHFuKuiqVW5B+CNHoGpVr6MqdzIyMli4cCGJiYlMnz6drVu3Ehoayq5du6hbty5Hjhyhqkd/94ISgT9jBK8Bl+PWLQY4DLxaTLEZYwLl119dNxBAhw6uIuixY65tSaDYHDp0iF9//RWACRMm0Lt3byZMmECHDh1455132LlzJ3V9U2+9SgKF8ScRdFXVu4BjAKq6D6gU0KiMMUXz449uwPfNN117wAB4/nmoXt3buMqJnTt38sYbbxAVFUX9+vX58MMPARgyZAiJiYns2bOHpKQkbr31Vho0aOBxtIXzZ7A4Q0RC8a1X7LvAzBavN6a0OXzYdf1cfjmccw6MHu0GhE2xOXr0KD179mTJkiUAXHDBBdx999107twZgEaNGhEdHe1liGfFn0TwEjANaCgi44AY4OGARmWMOXN33QVJSbB9u/vmb+sCF0lWVhaLFi0iKSmJjIwMXnjhBapUqUKbNm0YOHAggwcPpl27dkg5qLpaaCJQ1YkisgLohZs6OlhV1wU8MmNMwX74Af72N3j2WWjaFO6/H/74R+v+KaLPP/+c999/n8TERHbt2kXFihW57rrrUFVEhPfff9/rEIudPwvTvATUVdVXVfWVM0kCItJPRNaLyAYReSCffa4XkbUiskZEJp1B7MYEn6wsOHDA3ReBWbPgq69cu1076NrVs9DKqiNHjpCYmEh6ejoAs2fPZuLEifTo0YPJkyeze/dupk6dWi6++efHn+mjI4AbcGWppwGTVbXQ+Zu+cYXvgGuB7cAyYLiqrs21T2tgCnCNqu4TkYaququg57XpoyZoZWW52T9dusCECW5b9rUB5owcOHCAmTNnMnXqVGbNmsWRI0eYOXMm/fv3Z//+/YSFhREWVr6WXSlSrSFVfQ94T0TqAsOAf4jIuaraupBf7QJsUNWNviAmA9HA2lz73Aa86puJRGFJwJigs28fpKS4heFDQ91VwOedl/O4JQG/ZXftfPvtt7Rv356MjAyaNGnCrbfeytChQ+nevTsAtWvX9jZQD5zJ4vWtgAuB8wB/uofOAbblam8H8p63tgEQkc+BUOAxVZ2d94lE5HbgdqBE6m4YU2q89ho8/LAbD2jZEsaM8TqiMmXbtm1MmzaNqVOn0r59e1566SXatGnDAw88QFRUFF27diUkxJ9Z9OVboYlARP4JDAF+ACYDT6rq/mJ8/dZAD6AZsFBELsr7/Kr6JvAmuK6hYnptY0qfI0fg3/+Gnj3d1M9Ro+C661wSMH4bP34877zzDsuWLQPcehIXXHABACEhITzxxBNehlfq+HNG8ANwuaruOcPn/hHIXcKwmW9bbtuBJaqaAWwSke9wiWHZGb6WMeXHK6/A8eMuEdSv724mX6rK119/zaeffsp9992HiLBy5UpEhGeeeYYhQ4bQpk0br8Ms1fwZLA7BlZdoqapPiMi5QOPCag2JSAXcYHEvXAJYBtyoqmty7dMPN4A8QkTqA18CHVT1l/ye1waLTbnz2Wdu8Pedd9y6wHv32mpghcjMzCQtLe3kgO/GjRsJCQlh7dq1tG3blqysrBIr5lZWFLXW0Ku4WkPZ9XYP4UetIVXNBO4GPsWNKUxR1TUi8oSIDPLt9inwi4isBVKBewtKAsaUGwcPutXAwPX/p6XBjh2ubUngtLKysli4cCGjR4+mYcOG9OzZkxdffJG2bdvy1ltvsXPnTtq2bQtgSeAM+XNGsFJVLxWRL1X1Et+2r1X14hKJMA87IzBl3nffQefOrgvod79z3UAAlayEV15ZWVmkpaURHx/PRx99xE8//URYWBh9+vQhPDycBx98kFq1ankdZplQ1KUqrdaQMUX19deuEFz//tC6tZsG2qGDe8wSwCmyv/nHx8czdepUfv75Z6pUqUL//v2JjY1lwIABVLerp4uV1RoyJtBU4e673drAUVHuimBbEOYUmZmZpxz8d+3aRdWqVRkwYAAxMTEMGDCAatWqeR1muWW1howJhLQ0N///o4+gXj2YOBFq1HBJwADu4D9//nwSEhKYOnUqu3fvpmrVqlx33XXExsYSFRVlB/8S4tcFZar6LfBtgGMxpuw6eNDN+unRw3X5VK/uzgC2bXOJwC6EBNzBPzU1lfj4eKZNm8aePXuoVq0aAwcOJCYmhqioqFK7eEt5diZXFhtjcsvKgv373YE+IwPuuw+ee84lgg4d3NoAdgZARkYGKSkpJCQkMG3aNH755ReqV6/OwIEDiY2NpV+/flSpUsXrMINavolARCqranpJBmNMmXLtta67JynJJYNNm9yCMNmCOAlkZGQwb9484uPjSUxMZO/evdSoUYNBgwYRExND37597eBfihR0RrAIuFREPlDV35VUQMaUWkuWwCefwGOPuYN8bCzkXoYwdxIIQsePH2fu3LkkJCSQmJjIvn37qFmzJoMGDSI2NpY+ffqUu4qe5UVBiaCSiNwIXCEiQ/M+qKpTAxeWMaXEV1/BhRdCWBh8842rA3TTTW4K6J13eh2d544fP86cOXOIj48nKSmJ/fv3U7NmTaKjo08e/CtbhdRSr6BEcAfwW6A2MDDPYwpYIjDl26JFrt7PpEkwfLhLACNGBP28//T0dJKTk0lISCApKYkDBw5Qq1YtBg8eTGxsLL1797aDfxmTbyJQ1TQgTUSWq+qEEozJGO+8/z5UrOgO/F27wptvQt++7rEg7tM+duwYycnJxMfHM336dA4ePEjt2rUZMmTIyYN/pSBPkGWZP7OGPhCRMUB3X3sBMN5XMdSYsm3HDli6FAYPdu3s9WiHD3cF4G67zbPQvHbs2DFmz55NQkIC06dP59ChQ9SpU4eYmBhiYmLo1auXHfzLCX8SwWtARd9PgN8BrwN/CFRQxgSUas6Mnsceg7g42LnTzf2fMgXq1PE0PC8dPXqU2bNnEx8fz8cff8zhw4epW7cu119/PbGxsVxzzTVUrFjR6zBNMfOn6Nz/FJizonOmzFq+3PXzJya6Ad9t29w1AEG88MuRI0dOHvxnzJjB4cOHqVevHkOHDiUmJoaePXvawb8cKGoZ6iwRuSDXk7UEsoorOGOKS1xcHJGRkYSGhhIZGUlcXBwcPerq+syZ43Zq3hxq13ZrAWe3gzAJHDlyhISEBG644QYaNmzIsGHDmDt3LjfeeCNz5szhp59+4s0336RPnz6lLgmc9nM2RaOqBd5wNYa2AvNx4wObgZ6F/V6gbh07dlRj8po0aZK2aNFCU1JS9PjOnbp4wgRt0aKFxn3wgWr9+qp/+5vXIXru8OHDOmXKFI2NjdWqVasqoA0aNNBRo0bp3LlzNSMjw+sQC3XK53z8uKakpGiLFi100qRJXodW6gHLNb/jfH4PnLITVAba+26V/fmdQN0sERjdt0917dqc9gcf6Fv162tKSopr9+mj2qmTpqSkaEREhOrevZ6EWRocOnRIJ0+erMOGDdMqVaoooA0bNtQ777xTU1JSNDMz0+sQz0hERETO5+xz8nM2BSooERQ6RlDa2BhBEFq1CubPh9GjXfvOOyE+3hV1A7j7br549VU6Hz/uujGSk6F6dTI6dyYsLIysrODqyTx06BAzZ84kPj6eWbNmcfToURo1asSwYcOIjY3lqquuKrMreIWGhnLs2LFTuqsyMjKC8nM+UwWNEXj2zf5sb3ZGECS2blXN/rb67LOqIqq7d7v20qWq06apnjjh2hkZQf9N8eDBgzpx4kQdPHiwhoWFKaBNmjTRu+++WxcsWFDmvvnnJ9g/56KgqF1DpelmiaAc+vVX1aQk1U2bXHvpUnfgnzXLtffuVf3llwKfIhj7jg8cOKD//e9/NTo6WitXrqyANm3aVEePHq0LFy7UrKwsr0MsdsH4OReXIiUCXCmJAUBIYfuWxM0SQTmwe7dqbKzq5MmuvX27+6f4wguu/csvquPGqW7bdkZPO2nSJI2IiNCQkBCNiIgolweH/fv36/vvv68DBw7USpUqKaDnnHOOjhkzRj/77LNyefDPKxg+50AoKBH4cx1Bb+BW4DIgHnhXVdcXT6/VmbMxgjLq3nuhWjV3AVdWlqvhM2oU/P737vHFi+Hii4O6jEN+9u/fz/Tp04mPjyc5OZnjx4/TrFkzYmJiiI2N5bLLLiMkxJ+Z4CaYFWnxelWdC8wVkVrAcN/9bcBbwH/VSk0Ep59+gu++g+6+yiOPPebakya59rBhsHUrLFvm2hs2uCqeAKGh7sCfu17/ZZeVWOhlwb59+0hKSiIhIYHk5GQyMjJo3rw5d911F7GxsXTt2tUO/qbY+LVCmYjUA27ClZf4EpgIdANGAD0CFZwpRb75xl2Udc897gD+97/De++5C7MqVICaNSEzM2f/wYPd6l3Zpk499cAfxIu25Gfv3r0kJSURHx/P3LlzycjI4Nxzz2XMmDHExsbSpUsXxP5uJgD86RqaBrQFPgD+o6o7cz22PL9TjUCxrqFicvw4bNkCTZu6LpstW+DTTyEmBurWdYuu/+Uv8OWXbvWtl16CP/3JnQk0agTr1rl1ejt3dsXZzFn55ZdfSExMJCEhgblz55KZmcn5559/stunc+fOdvA3xaJIXUPAW6r6SZ4nrKyq6SWdBMwZyMiAH36Axo1dSYVvv4XHH4eHH4aICDfXfuBAV3P/ssvcN/5Ro+DSS10iaNIELr/cPQ+4+jwjRkCtWq4dHu7ZWyvr9uzZQ2JiIvHx8aSkpJCZmUmLFi245557iI2NpWPHjnbwNyUrv1Hk7Buw0p9tJXUrV7OGTpxQzb6sPytL9euvVXfscO30dNX331ddtcq1Dx5Uvece1YULXfvnn1V79VKdMcO1f/hBtV491bg411692s3EiY937VWrVFu2VE1Nde2dO1Xfe091zx7X/vVXN0vn+PGAvuVg9dNPP+n48eO1d+/eGhoaqoC2bNlS77//fl2+fLmeyL4mwpgAoYBZQ/me04tIYxHpCFQRkUtE5FLfrQdQNfApKgDS0+GXX3Lac+e6b8bZXnwRxo/PaQ8c6LpDst16Kzz5ZE77kUfgP//Jab/xxqnPd+edOfXtwS108txz7r6qW+kq+/kyM92smQm+NYCysuDmm90auQAnTrjnX73atStWhGPHcvrl69SB3/wGWrRw7WbN4IMPcgZhIyPdGUKPHq7duLF7/nr1XLtqVfc7pazAWFm2Y8cOXnnlFXr27EnTpk2544472LJlC/fffz8rVqxgw4YNPPPMM3YGYDxXUNdQX+AWoBnwfK7th4CHAhhT4Fx3HRw+7LpDAJ56yh1I+/Rx7Y8/doOed9zh2pGRrlslW3p6TlcJwOzZ7uB+yy2u/eijMGBAzvMtXer607O1bJlz4BVxiaRbN9euVMn1y0dEuHZYmJtp07Cha9eq5WLPVqcOpKWd2n7llZx2rVpuaUVTorZt28ZHH31EQkICX3zxBapKu3btePjhh4mJiSEyMtIO+qbU8WeweJiqflRC8RSqSIPFiYnuYH7DDa69bZub8dKkSfEEd/iwO6Dbqk1BZdOmTScP/kuWLAGgffv2xMTEMGzYMNq1a+dxhMac5WCxiNykqv8FzheRe/I+rqrPn+bXSrfs5QizNW9evM9fvXrxPp8ptb7//nsSEhL46KOPWLFiBQAdO3bk6aefZtiwYbRu3drjCI3xX0FdQ9V8P+3oZgywbt06EhISSEhI4JtvvgGga9euPPvsswwbNowW2eMzxpQx+SYCVX3D9/PxkgvHmNJDVVm9evXJg//atWsBuPLKK/n3v//N0KFDOffccz2O0piiK6hr6KWCflFVxxR/OMZ4S1X56quvTh78v/vuO0SE7t278/LLLzNkyBDOOeccr8M0plgV1DW0osSiMMZDqsqyZctO9vlv3LiR0NBQevTowV/+8hcGDx5M48aNvQ7TmIApqGvovZIMxJiSlJWVRVpaGtOmTWPatGls3bqVChUq0Lt3bx588EGio6Np0KCB12EaUyIK6hp6QVX/LCIfA/8zx1RVBxX25CLSD3gRCAXeVtVn8tlvGJAAdFZVKyRkAiI9PZ2UlBSmTp1KUlISu3fvpnLlyvTp04fHH3+c6Oho6tSp43WYxpS4grqGPvD9/NfZPLGIhAKvAtcC24FlIjJdVdfm2a8G8Cdgydm8jjEFOXz4MLNnz2bq1KnMnDmTgwcPUqNGDQYMGMDQoUPp168fNWrU8DpMYzxVUNfQCt/PBSJSCbgQd2awXlWP+/HcXYANqroRQEQmA9HA2jz7PQn8A7j3zMM35n/t3buXjz/+mKlTp5KcnMyxY8eoX78+sbGxDB06lF69elG5cmWvwzSm1Ci0+qiIDADGAz8AArQQkVGqOquQXz0H2JarvR3omue5LwWaq+pMEck3EYjI7cDtgE3XM6e1Y8cOEhMTmTZtGqmpqWRlZdGsWTNuv/12hg4dypVXXkmFCn4tv2FM0PHnf8ZzQE9V3QAgIhcAM4HCEkGBRCQEV8PolsL2VdU3gTfBlZgoyuua8uOHH35g2rRpTJ06lUW++lFt2rThvvvuY8iQIXTq1Mnq+pizpwrbt+dUIPj0U1ff69FHXWma5GT44gtX2r1CBVfU8eef3TKs4GqFHTgAHTu69s8/w6FD0KqVf6+/fTscPQrZV6l/9x20aVO879HHnxVFDmUnAZ+NuMJzhfkRyF3DoZlvW7YaQCQwX0Q249ZEni4itsaBOS1VZdWqVTz++ONcfPHFtGrVinvvvZf09HT+/ve/s2bNGr799lueeuopW9ClvFJ1lXnBrY6XlOSWRAW3fccOt8/p/PyzqwacXYF46VIYMsQdwMEVkRwwwO0H8PzzcO657nXAVf4dNy5ndb1Fi+Af/3BLrwK8/jpce23O673wQk4BSnAJ48orc9qjR+ckDYDbboOrrz61feONOe2FCwv80xRJfvWpgaG+2+vAJ7hv7iOAGcBr+f1ert+vgEsaLYBKwNdARAH7zwc6Ffa85Wo9AlOozMxMTUtL03vvvVdbtWqlgIqIXnXVVfr888/rpk2bvA7RnInjx1WPHctpb9igunGju3/ihOoNN6hOmJDz+LBhqhMnuvu7d6uGhKi+9JJr//ijW3Pj7bdd+/vvXfvdd1172TLVqlVVZ8927dRU93hysmsnJ6tGRua8/qefqnbqpLpli2uvWqU6frzqgQM58eReN+LECdXMzJz299/nrPehqrp+vercuTntJUtUp0/Pab/9tuqDD+a033hD9bHHctoLFpz6+wcPalFQwHoEBR3I3y3olt/v5XmO/sB3uPGFsb5tTwCDTrOvJQKjqqpHjhzR6dOn68iRI7Vhw4YKaMWKFbVv3776xhtv6E8//eR1iCZbZuapB6jPPlP95JOc9oMPqj73XE77ggtUb7wxp928ueott+S0r7hC9Z//dPePH1ft2FH1lVdcOz1ddexY1S++cO0TJ1RXrFDdtcu1d+92+27Y4Nrbt6v+9a+qa9a49q+/qn77bdAuvlRQIii0DHVpY2sWl0979uxhxowZJCUlkZyczJEjR6hZsyb9+/cnOjqaqKgoamUvk2kCKysrp7sjNRU2boSRI1171CjX9z1vnmsPGOC6Y7780rV794b9+yH7/+jgwdCgAbz1lmu/8YYr+z7IdxlSYqJbp6N9+xJ4Y8GtSGsWi0gYMBKIAMKyt6vq74stQhOUfvjhB5KSkkhMTOTzzz/nxIkTNGvWjFtvvZXo6GiuvvpqKtnaDsXr6FHYswfOOQdCQuD77+GrryAmxvV9P/ggvPNOTj/5lCluwaTsRNC5s1vTOtvIke45s73yinvubImJp77+qFGntvOWhjfeyO9UQXO6bOJxc/1/wI0RJAMvFvZ7gbpZ11DZlZWVpUuWLNGHHnpIIyIiFHddirZv314feeQRW7u3OOzZ49ap3r/ftVNSVAcOVP3uO9f+979dj/Deva79z3+69qFDrj1zpuqjj+b0fe/ff2ofuSmzKKBryJ/po61UNVZEolX1PRGZBHwWmLRkypv09HRSU1NJSkpi+vTp7Nixg9DQULp3785tt93GoEGDrI5/URw8CJMnu9kp55/vvt3Hxrqpjdde61bL27w5Zy3q3r3h7bch+4K6m2+G/v2hShXX7t/f3bJZd1xQ8CcRZC/Su19EIoGfgIaBC8mUdfv27eOTTz4hKSmJ2bNnc+jQIapVq0a/fv2Ijo5mwIAB1M3dvWAKdvSoW2K1dm23xvYTT7hph/36uf74UaPgpZfcdMTLLoNly+Cii9zvXnkl+BbRAdw63JGROe1GjU5dV9sEJX8SwZsiUgd4BJiOW7HskYBGZcqc9evXM2PGDGbMmMFnn31GVlYWjRs3Zvjw4URHR3PNNdcQFhZW+BMZmDjR9cNHRbk58Q0bur74F15wg7gvv+z269fPzXNfvz7noqNq1aCTXYpjzlB+fUal9VZaxwgmTZqkERERGhISohERETpp0iSvQzojZxp/enq6zps3T//yl79o69atT+nvf+ihh3TRokWalZVVpNco13LPpX/7bdUXX8xpR0SoDh6c0x4/XnX+/Jx27rnrZ8g+g+DF2VxHcHIHqAe8DKzELVbzAlCvsN8L1K00JoJJkyZpixYtNCUlRY8fP64pKSnaokWLMvOfzN/4d+/ere+//77GxsZqzZo1FdDKlStrVFSUvvrqq7ol+0KcIrxGuXT0qOo33+S0H3nEXbiULTpatUePnPbOnaoZGcUeRlB/BqbIiWAOriuohe/2MDC3sN8L1K00JoKIiAhNSUk5ZVtKSopGRER4FNGZyS/+du3a6apVq/Spp57SK664QkVEAW3cuLH+4Q9/0KSkJD18+HCRXqOs/I3OyM8/q374Yc5VqPfeqxoW5i5oUlV94QXVxx/P+Waf58wpUILqMzD/o6BEUOgFZSKyWlUj82xbpaoXFVPv1BkpjReUhYaGcuzYMSpmz8wAMjIyCAsLIyu7Lkopljv+Y8eOMX/+fKZPn87rr79+cp+OHTsycOBArrvuOi655BJCQvwpU3X618hWlv5GhVq/3hUnq1rVXTR1xx2wbh1ceKEbrN20yfXpe1j+utx/BqZARbqgDEgWkd8AU3ztGODT4gquPAgPDyctLY2ePXue3JaWlkZ4eLiHUfmvVatWPPTQQ2zYsIE5c+bw66+/UrlyZWrUqMHzzz9P//79adq0aZFeo6z/jU7rxAl3Udbnn0O3bjB1qitidv317sKr7AHc9u1LxZWz5fIzMMUjv1MFXIXRg76fJ4BM3+0EcDC/3wv0rTR2DZW1vteMjAxduHChPvjgg9qhQ4eTA70NGzbUUaNG6dNPP63nnXdescZf1v5Gp5Vdo2b/ftUWLVz3jqrr4nnpJdVSXgOpXHwG5qxRlDGC0nYrjYlAtfTPxvjxxx91woQJGhMTo7Vq1VJAQ0NDtXv37vr000/r008/re3atQto/KX9b/Q/cs/Oueoq1ZiYnPaDD6pOnVryMRVRmfsMTLEpKBH4VXRORAYB3X3N+ao6o5hPTPxWGscISqPMzEwWLVrErFmzmDVrFl999RUATZs2pV+/fvTv35/evXtbIbfcfvkF6tVz9+++G5YscRdngZvDX7s23HKLR8EZUzRFLTr3DNAZmOjb9CcRuVJVHyzGGE0x2LlzJ7Nnz2bWrFnMmTOH/fv3ExoaypVXXsnTTz9NVFQU7du3twVbsh075gZvReCRR9zBft8+t9pU585Qp467oEsE/vxnr6M1JmD8GSzuD3RQ1RMAIvIe8CVgicBjmZmZLF68mE8++eSUb/1NmjRh6NChREVF0bt3b2rXru1pnKXGnj2udk7FihAXBzfd5Eoqt2jhZvTUqQMZGS4RjBjhdbTGlBh/V/OuDez13be+BA9t2rSJ5ORkkpOTmTdvHgcOHCA0NJQrrriCp556iv79+9u3/tNZvBiuusotbdi/P3To4JYOzC57ceWVpy4jaEwQ8ScRPAV8KSKpgODGCh4IaFTmpIMHD5Kamnry4L9hg1s+unnz5sTGxtK3b1/71n866eluLn+3bq5Oz6WXwr335iz+HR4Ojz/ubYzGlBIFJgIRCcFNF70MN04AcL+q/hTowIJVVlYWy5cvJzk5mU8//ZTFixeTlZVFtWrV6NmzJ2PGjKFPnz60adPGvvXnlZwMu3a5Lp/Kld3KWm3buscqVYKnnvI2PmNKKX+uLF6e30izF8rjrKHNmzczZ84ckpOTmTt3Lvv370dE6NixI3369KFPnz5cfvnltlpXXhkZsHIldO3q2jExsGaNu6IXcgZ6jTFFvrJ4roj8FfgQ+DV7o6ruzf9XTEEOHTp0srtnzpw5fPfddwA0a9aMoUOH0qdPH3r16kX9+vU9jrQUyv7iIuLq8j/zjFszt0EDt0xi7nUOLAkY4xd/EsENvp935dqmQMviD6d8Sk9PZ/HixcybN4958+axdOlSMjMzqVq1Kj169OCPf/wjffr04cILL7TunoIsW+bKN3z0kevzv+UWdzaQPT7SuLGX0RlTZhWaCFTV1hE8Q1lZWaxcuZKUlBTmzZtHWloaR48eJSQkhE6dOnHvvfdy7bXXcsUVV1DZwyJkpV5WlluGsVEjt8TiBRdAu3ZuO7j2BRd4G6Mx5YA/F5SFAX8EuuHOBD4DxqvqsQDHVmaoKuvWrTv5jX/+/PkcOHAAgIiICG677TZ69epF9+7dbXbPmcjIgAcegOholwjq1oWZM72Oyphyx5+uofdxhed86+NxI/ABEBuooMqCzZs3n/zGn5KSwk8/uYlU559/PjExMfTq1YtrrrmGRrYe7Jl59113sdenn7o5/gsXwnnneR2VMeWaP4kgUlXb5WqnisjaQAVUWv34448sWLCA+fPnM2/ePDZu3AhAw4YNueaaa+jVqxe9evWiRQvrSTtja9a4g3316q6ef1YW7N3r6v7Y39OYgPMnEawUkctUdTGAiHQFytf8zdPYtm0b8+fPZ8GCBSxYsODkhVw1a9akR48ejBkzhl69ehEREWEDvEXx5Zdu4Pfdd93g7w03uJsxpsT4kwg6Al+IyFZf+1xgvYisAlRVvV9xoxhs3rz55Df+BQsWsGnTJgBq167NVVddxZ133snVV19Nhw4dCA0N9TjaMiwry833v/RSV+itQwd46y0YMMDryIwJWv4kgn4Bj6KEqSqbNm06edCfP38+W7e6PFe3bl26d+/OmDFj6NGjBxdddJEd+Itq2TJYvRpuvRVCQ6FmTahWzT0mAn/4g7fxGRPk/FqPoDQ52yuLt27dSnJy8smD//bt2wGoX78+V199NVdffTU9evQgIiLijNfjNaexcyc0aeLu33GHm/u/Y4er/GmMKXFFvbK4XIiLi+OBBx6gYcOG9OjR4+TBv127dtbHX9wmToTf/c4t6N66NTz2GDz7rCUBY0qpoEkEI0aMIDo6mrZt29qBv7ilp8Mbb8Bll0GXLtCjBzz5pKv9D3bFrzGlXND0gTRu3NhKOBSnfftylnEMCYFHH4WpU137nHNg7Fho2NC7+IwxfguaMwJTzPr2dSt5ffGF6/JZtQqaNfM6KmPMWQiaMwJTRBs2wO23w/Hjrv2vf8H48TmPWxIwpswKaCIQkX4isl5ENojI/6xqJiL3iMhaEflGROaJiNUSKE0OHYKDB939TZtgyhT3zR+ge3doXy4uITEm6AUsEYhIKPAqEAW0A4aLSLs8u30JdPJdlJYA/DNQ8ZgztGePG+R9/XXXvuYa+P576NjR27iMMcUukGcEXYANqrpRVY8Dk4Ho3DuoaqqqHvE1FwPWv+Clxx93Uz0B6td3A8B9+rh2aKhb/MUYU+4EMhGcA2zL1d7u25afkcCs0z0gIreLyHIRWb579+5iDDHIrVkDr76a09661Y0FZLvvPrjkkpKPyxhTokrFYLGI3AR0Ap493eOq+qaqdlLVTg3sW2nRqOYs9/jhh/B//wf797v222/Df//rWWjGGG8EMhH8CDTP1W7m23YKEekNjAUGqWp6AOMx33zjyjp/8YVr33MP/PhjzlKPdo2FMUEpkIlgGdBaRFqISCXgN8D03DuIyCXAG7gksCuAsQSngwfhuutciWeAli1dV09mpmvXru1q/htjglrALihT1UwRuRv4FAgF3lHVNSLyBLBcVafjuoKqA/G+K363quqgQMUUFN57D06ccJU+a9Rw5R+yD/zVq8O0ad7GZ4wpdQI6RqCqn6hqG1W9QFXH+bb9zZcEUNXeqtpIVTv4bpYE8hEXF0dkZCShoaFERkYSFxfnHsjKcou7ZPvwQ1f0DVxXz5w5cNttJR+wOSv5fs7GBJCVmCgD4uLiGDt2LBMmTKBbt24snj2bEWPGADB8/Xr4xz9g2zY35TMuztX7N2VO3s85LS2NkSNHAjB8+HCPozPlmqqWqVvHjh012ERERGhKSoprJCWphobqkgkTNCIiQvWrr1QTElQzM70N0hTZKZ+zT0pKivucjSkiXJf8aY+rQbMwTZm1dy/v1KvHzTNmUGHAALe4y+uvk3HrrYS1bk1WVpbXEZpiEhoayrFjx6iYa92GjIwMwsLC7HM2RVbQwjSl4joCk8fBg7Bunbtfowa9KlXihzlzXLtpU3jySdK2bCE8PNy7GE2xCw8PJy0t7ZRtaWlp9jmbgLNEUBqouvn82QYNghtucNsrVmTRO+8QNX06qampZGRkkJqaysiRIxk7dqx3MZtiN3bsWEaOHGmfsylxNljslawsV78H4N57YcIE2LXL1fZ/9FEICzu5629++1s0JITRo0ezbt06wsPDGTdunA0gljPZn6d9zqak2RhBScku7RASAh9/DDffDKtXu9W8Fi1yV/2OGHFKAjDGmOJiYwReW7PGXdWbmura558P11/vLvYCuPxyGDXKkoAxxhOWCAIhPd1V7vzwQ9du2RI6dMg50F90kVvsvWVLz0I0xphslgiKy8qVMHOmu1+pEsye7bp7AKpUcaUdrrzSu/iMMSYfNlh8tvbtcwf/Xr1c+29/g2+/hf79XWmHlSvd4u7GGFPK2RmBvw4dgpQUN9sH4MknISoKfv3VtV94AZYuzSnlbEnAGFNGWCLIT3o6fPZZzqIt06a5b//ffuvad98NCxbk9Pu3agV163oSqjHGFIUlgmyZmW7Bli1bXPvLL6F7d5g3z7X79nX9/uef79otW7rZPtnXAhhjTBkVvIlA1Q3mrlnj2vv3u8HcDz5w7Y4dITExZwygUSOXDKpV8yJaY4wJmODqyN60Cfbsgc6dXbt3b3dw/+ADV8J59mzo5LveomJFiI72LlZjjCkhwZUIbr7Z9f1nD+pOmQIXXJDzeN++3sVmjDEeCa5E8NxzbvnGbD16eBaKMcaUFsGVCLp08ToCY4wpdYJ3sNgYYwxgicAYY4KeJQJjjAlylgiMMSbIWSIwxpggZ4nAGGOCnCUCY4wJcpYIjDEmyFkiMMaYIGeJwBhjgpwlAmOMCXKWCIwxJshZIjDGmCBnicAYY4KcJQJjjAlylgiMMSbIBTQRiEg/EVkvIhtE5IHTPF5ZRD70Pb5ERM4PZDzGGGP+V8ASgYiEAq8CUUA7YLiItMuz20hgn6q2Av4N/CNQ8RhjjDm9QJ4RdAE2qOpGVT0OTAai8+wTDbznu58A9BIRCWBMxhhj8gjkmsXnANtytbcDXfPbR1UzReQAUA/Yk3snEbkduN3XPCwi688ypvp5nzsI2HsODvaeg0NR3vN5+T1QJhavV9U3gTeL+jwislxVOxVDSGWGvefgYO85OATqPQeya+hHoHmudjPfttPuIyIVgFrALwGMyRhjTB6BTATLgNYi0kJEKgG/Aabn2Wc6MMJ3PwZIUVUNYEzGGGPyCFjXkK/P/27gUyAUeEdV14jIE8ByVZ0OTAA+EJENwF5csgikIncvlUH2noODvefgEJD3LPYF3BhjgptdWWyMMUHOEoExxgS5cpkIgrG0hR/v+RYR2S0iX/luf/AizuIiIu+IyC4RWZ3P4yIiL/n+Ht+IyKUlHWNx8+M99xCRA7k+47+VdIzFTUSai0iqiKwVkTUi8qfT7FNuPms/32/xf86qWq5uuIHpH4CWQCXga6Bdnn3+CIz33f8N8KHXcZfAe74FeMXrWIvxPXcHLgVW5/N4f2AWIMBlwBKvYy6B99wDmOF1nMX8npsAl/ru1wC+O82/7XLzWfv5fov9cy6PZwTBWNrCn/dcrqjqQtxMs/xEA++rsxioLSJNSia6wPDjPZc7qrpTVVf67h8C1uEqEuRWbj5rP99vsSuPieB0pS3y/iFPKW0BZJe2KKv8ec8Aw3ynzgki0vw0j5cn/v5NypvLReRrEZklIhFeB1OcfF24lwBL8jxULj/rAt4vFPPnXB4TgTm9j4HzVbU9MIecMyJTfqwEzlPVi4GXgURvwyk+IlId+Aj4s6oe9DqeQCvk/Rb751weE0EwlrYo9D2r6i+qmu5rvg10LKHYvOLPv4NyRVUPquph3/1PgIoiUt/jsIpMRCriDooTVXXqaXYpV591Ye83EJ9zeUwEwVjaotD3nKfPdBCu77E8mw7c7JtRchlwQFV3eh1UIIlI4+yxLhHpgvv/XZa/4OB7PxOAdar6fD67lZvP2p/3G4jPuUxUHz0TWjpLWwSUn+95jIgMAjJx7/kWzwIuBiISh5s9UV9EtgOPAhUBVHU88AluNskG4AhwqzeRFh8/3nMMcKeIZAJHgd+U8S84AFcCvwNWichXvm0PAedCufys/Xm/xf45W4kJY4wJcuWxa8gYY8wZsERgjDFBzhKBMcYEOUsExhgT5CwRGGNMkLNEYIKeiDwmIn/1Og5jvGKJwBhjgpwlAhOURGSsiHwnImlA21zbO4jIYl9xvmkiUse3fYyvRvw3IjLZt62LiCwSkS9F5AsRaevbXlVEpvj2nyZuzYtOvsf6+H5npYjE+2rK5I6rgogsE5EevvbTIjKuRP4oJmjZBWUm6IhIR+A/QFfc1fUrcetT/EtEvgFGq+oC35XZNVX1zyKyA2ihqukiUltV94tITeCI78ru3sCdqjrM183UWlVHiUgk8BWuTv5mYCoQpaq/isj9QGVVfSJPfBG48uijgWeBrr7y4sYERLkrMWGMH64CpqnqEQARme77WQuoraoLfPu9B8T77n8DTBSRRHKqPdYC3hOR1oDiK/cAdANeBFDV1b7kAi4ZtAM+95WKqQQsyhucrzzIB8AM4HJLAibQLBEY458BuBXCBgJjReQi4EkgVVWH+GrHzy/kOQSYo6rD/Xi9i4D9QMOzDdgYf9kYgQlGC4HBIlJFRGrgDu6o6gFgn4hc5dvvd8ACEQkBmqtqKnA/7kyguu9ndrnjW3I9/+fA9QAi0g53UAdYDFwpIq18j1UTkTZ5gxORoUBdXOJ5WURqF8ebNiY/dkZggo6qrhSRD3FrO+/ClfHONgIYLyJVgY24SpahwH99XUcCvOQbI/gnrmvoYWBmrud4zbd9LfAtsAZXGnm3iNwCxIlIZd++D+PWpQXAV1f+GaCXqm4TkVdw3UwjMCZAbLDYmGImIqFARVU9JiIXAHOBttbXb0orOyMwpvhVBVJ9K00J8EdLAqY0szMCY4wJcjZYbIwxQc4SgTHGBDlLBMYYE+QsERhjTJCzRGCMMUHu/wEIbEdWXC43LwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "space = np.linspace(0, 2.5, mplot)\n",
    "\n",
    "plt.plot(space1, H(mle[0] + space1 * mle[1]), color='black')\n",
    "plt.plot(space2, H(mle[0] + space2 * mle[1]), color='black', linestyle='dashed')\n",
    "plt.plot(space, H(lower90), color='red', linestyle='dotted')\n",
    "plt.plot(space, H(upper90), color='red', linestyle='dotted')\n",
    "plt.scatter(x, y / m, facecolor='white', edgecolor='black')\n",
    "plt.ylim(0, 1)\n",
    "plt.xlabel('dosage x')\n",
    "plt.ylabel('probability of event A')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "202265b5",
   "metadata": {},
   "source": [
    "#### (d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "d2b1e70e",
   "metadata": {},
   "outputs": [],
   "source": [
    "xnew = 2.5\n",
    "pxnew_sample = H(sample[:, 0] + sample[:, 1] * xnew)\n",
    "density = kde(pxnew_sample)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e383a68",
   "metadata": {},
   "source": [
    "#### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "c17615ae",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAALEwAACxMBAJqcGAAALalJREFUeJzt3Xl4VdXVx/HvysAko4RBmQcREmIYLkqgAg4IooIWEVQUASuat1XESqnaSqGKRcVqHRAEEZlSlCgCAopVkTAYgjJPMonKICoQIYSQ9f5xLzSFABeSm32H9Xme+3iHk3N+R5KsnL332VtUFWOMMZErynUAY4wxblkhMMaYCGeFwBhjIpwVAmOMiXBWCIwxJsLFuA5wruLi4rRu3bquYxhjTEhZvnz5j6papaDPAlYIRKQWMBGoBigwRlVfPGmbDsD7wFbfWzNUddiZ9lu3bl0yMjKKPK8xxoQzEdl+us8CeUWQCzyiqpkiUg5YLiIfqerak7ZbqKo3BjCHMcaYMwhYH4Gq/qCqmb7nB4F1QI1AHc8YY8z5KZbOYhGpCzQHlhbwcbKIfC0iH4pIQnHkMcYY818B7ywWkbLAu8BAVT1w0seZQB1VzRKRLsB7wCUF7OM+4D6A2rVrBzawMcZEmIBeEYhILN4iMFlVZ5z8uaoeUNUs3/M5QKyIxBWw3RhV9aiqp0qVAju9jTHGnKeAFQIREWAcsE5VR51mm+q+7RCRy3159gUqkzHGmFMFsmmoLXAXsEpEvvK99xhQG0BVRwO3Ag+ISC5wGOilNh2qMcYUq4AVAlX9ApCzbPMy8HKgMhhjQsPBgwdZunQpmzZt4scff6RSpUpcdtlltG3blujoaNfxwl7I3VlsjAle33//PaNGjeLDDz9k27ZtXHDBBdSqVYvExESSkpJo0qQJlStXJjs7m+3bt7N48WLS09NZuXIleXl5p+yvatWqDBkyhJSUFEqWLOngjCKDhFpLjMfjUbuz2JjgM2PGDHr37k1OTg6dOnWiUaNGHD58mK1bt/L111+ze/fuU76mbNmytG7dmjZt2tC2bVuaNm1KXFwcP/30E1988QVjx45l/vz5xMfHM336dOLj4x2cWXgQkeWq6inwMysExpjCGjduHL/73e+44oormDRpEg0aNDhlm927d7N582Z+/vlnSpUqRfXq1WnSpMlZm37mzJlD3759ycrKYtq0adx0002BOo2wZoXAGBMw6enptG/fnmuuuYYZM2ZQpkyZIj/G999/T7du3VixYgVvvvkmd911V5EfI9ydqRDYNNTGmPO2b98+evToQZ06dZg2bVpAigDAxRdfzCeffEL79u3p06cPqampATlOpLLOYmPMefvTn/7Enj17+PLLL6lYsWJAj1WuXDk++OADOnfuTO/evalUqRLXXXddQI8ZKeyKwBhzXtLT0xk3bhyDBg2iWbNmxXLMMmXK8MEHHxAfH89tt93GunXriuW44c76CIwJQ1lZWTz77LN89NFH5ObmctNNN/GHP/yhyP5qz8vLo1WrVuzZs4d169ZRtmzZItmvv7Zv387ll19O+fLlyczMpFy5csV6/FBkfQTGRJDvv/+edu3aMXz4cABiYmJ48skn8Xg8rFq1qkiOkZqaSmZmJiNGjCj2IgBQp04d3nnnHbZs2UJKSkqxHz/cWCEwJozk5ORw/fXXs3HjRmbNmkV6ejrp6eksXLiQQ4cO0a5du0I3p+Tk5PD444+TlJTEHXfcUUTJz92VV17JX//6VyZNmsS0adOc5QgHVgiMCSPDhw9n5cqVTJ06lS5dupx4v23btqSnp1OyZEluuOEG9uzZc97HGD16NFu3bmXEiBFERbn9FfLEE0/QqlUrBg4cyC+//OI0S0hT1ZB6tGzZUo0xp1q5cqVGR0drnz59TrvNsmXLtFSpUtq5c2c9duzYOR9j//79GhcXp1dffbXm5eUVIm3RWb58uUZFRWlKSorrKEENyNDT/F61KwJjwsTf//53ypQpw6hRBc76DkCrVq0YNWoUc+fO5aWXXjrnY4wcOZIff/yRkSNH4ptB3rkWLVqQkpLC6NGj2bRpk+s4oel0FSJYH3ZFYMypNmzYoCKiQ4YMOeu2eXl52rVrVy1RooSuWLHC72Ps3LlTS5curb169SpE0sDYtWuXli5dWu+++27XUYIWdkVgTHj7xz/+QalSpXj44YfPuq2IMG7cOCpXrsztt9/OoUOH/DrG0KFDyc3N5amnnips3CJXrVo1UlJSmDRpkl0VnAcrBMaEuF9++YUpU6bQp08fqlat6tfXxMXF8fbbb7NhwwYGDRp01u1Xr17N+PHjSUlJoX79+oWNHBCPPvooJUqUOGPTmCmYFQJjQty0adPIzs7m3nvvPaevu+aaaxg8eDCvv/46aWlpp93u6NGj9OnThwsvvJAnnniisHEDplq1avTo0YPJkyeTlZXlOk5IsUJgTIgbP348l112GS1atDjnrx02bBgej4d+/fqd9mazv/3tb2RmZjJmzBji4uIKGzeg7r//fg4ePMjUqVNdRwkpVgiMCWGrV6/myy+/pF+/fuc1iqdEiRJMnz6dMmXK0LFjR9auXfs/n//rX//iqaeeom/fvtxyyy1FFTtgkpOTadq0Ka+//rrrKCHFCoExIWzq1KlER0cX6g7funXr8vHHH5OXl0eLFi0YMmQIkydP5rbbbuPBBx/klltuCZlfrCLCvffey/Lly9mwYYPrOCHDCoExIUpVmT59Oh06dKBKlSqF2leTJk1YuXIlXbt25R//+Ae9e/dm/vz5PPbYY6SmphIbG1tEqQOve/fuALz77ruOk4QOm33UmBC1cuVKkpKSGD16NAMGDCiy/f7888/s3LmTunXrhuysnsnJyeTk5LB8+XLXUYKGzT5qTBh65513iIqKKvK2+0qVKpGYmBiyRQC8VwWZmZls3brVdZSQYIXAmBD17rvv0r59e7/vHYgkx5uHZsyY4ThJaLBCYEwI+uabb1i7di0333yz6yhBqV69eiQkJDBv3jzXUUKCFQJjQtDs2bMBuOGGGxwnCV4dO3Zk4cKFZGdnu44S9KwQGBOCZs2aRePGjWnQoIHrKEGrY8eOZGdn88UXX7iOEvSsEBgTYg4ePMinn37KjTfe6DpKUGvXrh2xsbF8/PHHrqMEPSsExoSYjz76iKNHj1ohOIuyZcuSnJzMRx995DpK0LNCYEyImT17NhUrVqRNmzauowS9jh07smLFCvbt2+c6SlCzQmBMCMnLy2P27Nl07tw5pO72deXKK69EVVmyZInrKEHNCoExIWT58uXs3r3bmoX81KpVK6Kjo0lPT3cdJagFrBCISC0R+Y+IrBWRNSLyUAHbiIi8JCKbRWSliJz7PLrGRJBZs2YRFRVF586dXUcJCWXKlKF58+YsWrTIdZSgFsgrglzgEVWNB1oD/yci8Sdtcz1wie9xH/BaAPMYE/JmzZpFcnIylStXdh0lZLRt25Zly5Zx9OhR11GCVsAKgar+oKqZvucHgXVAjZM26wZM9K2tvASoKCIXBSqTMaFs586dZGZmWrPQOWrTpg2HDx/m66+/dh0laBVLH4GI1AWaA0tP+qgG8G2+1zs5tVggIveJSIaIZOzduzdgOY0JZrNmzQKgW7dujpOEluOjq6x56PQCXghEpCzwLjBQVQ+czz5UdYyqelTVU9h5140JVTNnzqRhw4Y0btzYdZSQUrNmTWrVqmUjh84goIVARGLxFoHJqlrQNIDfAbXyva7pe88Yk09WVhYLFiyga9eu57UkZaRr2bKlrU1wBoEcNSTAOGCdqo46zWYzgbt9o4daA/tV9YdAZTImVM2fP5+cnBy6du3qOkpI8ng8bNq0if3797uOEpQCeUXQFrgLuFpEvvI9uojI/SJyv2+bOcAWYDMwFkgJYB5jQtbMmTOpVKkSbdu2dR0lJLVs2RKAFStWOE4SnGICtWNV/QI44zWsetfJ/L9AZTAmHBw7doxZs2bRpUsXYmIC9iMb1o4XgoyMDDp06OA2TBCyO4uNCXKLFy9m37591ixUCFWqVKF27drWT3AaVgiMCXIzZ84kNjaWTp06uY4S0qzD+PSsEBgTxFSVmTNn0qFDBypUqOA6Tkhr2bKldRifhhUCY4LYqlWr2LBhA7fccovrKCGvefPmgPf/qflfVgiMCWJTp04lOjqaW2+91XWUkJeYmAhYISiIFQJjgpSqMnXqVK677jrsjvrCq1mzJhUqVGDlypWuowQdKwTGBKnFixezfft2br/9dtdRwoKIkJiYaFcEBbBCYEyQeuuttyhdujQ333yz6yhhIzExkdWrV+O9hckcZ4XAmCB04MABJk+ezO233065cuVcxwkbiYmJ7N+/n2+//fbsG0cQKwTGBKG3336bX3/9lQceeMB1lLBy2WWXAdZhfDIrBMYEGVXltddeo2XLlng8HtdxwkrTpk0BKwQns4lLjAkys2bNYs2aNbz55puuo4SdChUqULt2bRs5dBK7IjAmiKgqQ4cOpX79+tx5552u44QlGzl0KisExgSRmTNnkpmZyV/+8hdiY2NdxwlLiYmJrF+/npycHNdRgoYVAmOCxK+//srDDz9Mo0aN6N27t+s4YSsxMZHc3Fw2bNjgOkrQsEJgTJB48skn2bp1K2PGjLF1BwLIppo4lRUCY4LAnDlzeOGFFxgwYADt27d3HSesXXrppcTExFiHcT5WCIxxbMWKFfTs2ZNmzZrx3HPPuY4T9kqUKEGTJk3siiAfKwTGODRv3jw6dOhApUqV+OCDDyhbtqzrSBHBRg79LysExjiQlZXFI488QpcuXahXrx6LFi3i4osvdh0rYiQmJvLtt9/yyy+/uI4SFKwQGFOMsrOzefHFF2nYsCGjRo3id7/7HQsXLqRWrVquo0WU4x3Ga9ascZwkOFghMKYY5OTk8Prrr9OwYUMGDhxIQkIC6enpjB492iaVcyA+Ph6wQnCcjVEzJsA+/vhjBgwYwJYtW0hOTubtt9/mqquuch0rotWpU4cyZcpYIfCxKwJjAiQvL4/BgwfTsWNHYmNjmT17NosWLbIiEASioqKIj49n7dq1rqMEBSsExgTA0aNHueuuu3j22WcZMGAAmZmZdOnSBRFxHc34xMfH2xWBjxUCY4qYqpKSksKUKVMYMWIEr732GmXKlHEdy5wkISGBH374gZ9//tl1FOesEBhTxF566SXeeOMNHnvsMYYMGWJXAUEqISEBwJqHsEJgTJFavXo1jz76KN26dWP48OGu45gzOF4IrHnICoExRebYsWP079+fChUqMHbsWKKi7McrmNWuXdtGDvnY8FFjisiYMWNYtmwZU6ZMoUqVKq7jmLOwkUP/ZX+yGFMEfv31V4YNG8aVV15Jr169XMcxfkpISLArAqwQGFMkXn75ZXbt2sXTTz9tncMhJD4+3kYOEcBCICLjRWSPiKw+zecdRGS/iHzle/w1UFmMCaRDhw4xcuRIunTpwm9+8xvXccw5sA5jr0BeEUwAOp9lm4Wq2sz3GBbALMYEzMSJE/npp5/485//7DqKOUc2hNQrYIVAVT8HfgrU/o0JBnl5ebz44ot4PB7atm3rOo45RzZyyMt1H0GyiHwtIh+KSMLpNhKR+0QkQ0Qy9u7dW5z5jDmj+fPns379egYOHGh9AyHo+MghKwTuZAJ1VDUJ+Bfw3uk2VNUxqupRVY8NyzPBZMyYMVStWpUePXq4jmLOU0JCgjUNuTqwqh5Q1Szf8zlArIjEucpjzLnas2cPH3zwAXfffTclSpRwHcecJxs55LAQiEh18V1Li8jlviz7XOUx5lxNmjSJ3Nxc+vbt6zqKKQQbORTY4aNTgcXApSKyU0T6i8j9InK/b5NbgdUi8jXwEtBLVTVQeYwpSqrK+PHjueKKK06sdmVCk40cCuAUE6p6+1k+fxl4OVDHNyaQ1q1bx5o1a3j5ZfsWDnW1a9fmggsusCsCY8y5SUtLA+Dmm292G8QUWlRUFE2aNLFCcDYicpOIWNEwxmfGjBm0bt2aGjVquI5iikCkzznk7y/3nsAmERkpIo0DGciYYLd9+3YyMzO55ZZbXEcxRSQhIYFdu3bx00+ReQ+sX4VAVXsDzYFvgAkisth3k1e5gKYzJgi99957AFYIwsjxDv9I7TD2u7lHVQ8A7wDTgIuAW4BMEflDgLIZE5TS0tJo2rQpl1xyiesopohE+hBSf/sIuolIGvApEAtcrqrXA0nAI4GLZ0xw2bt3LwsXLrSrgTBzfORQpF4R+Dt89LfAC76J5E5Q1UMi0r/oYxkTnGbOnEleXp4VgjAT6SOH/G0a2nVyERCRfwCo6oIiT2VMkEpLS6Nu3bo0a9bMdRRTxCJ55JC/haBjAe9dX5RBjAl2Bw4c4KOPPuKWW26xmUbDUCSPHDpjIRCRB0RkFdBYRFbme2wFVhZPRGOCw4cffkhOTg6//e1vXUcxARDJI4fO1kcwBfgQGAEMyff+QVWNvLJpIlpaWhpVq1YlOTnZdRQTAPlHDkXakqNnaxpSVd0G/B9wMN8DEbkwsNGMCR7Z2dnMnj2bbt26ER0d7TqOCYBInnPInyuCG4HlgAL5G0YVqB+gXMYElQULFpCVlWXNQmHs+Gpl1jR0ElW90fffesUTx5jglJaWRvny5bn66qtdRzEBFB8fz7x581zHKHb+3lDWVkQu8D3vLSKjRKR2YKMZExxyc3N5//33ueGGG2wlsjAXqSOH/B0++hpwSESO30n8DfB2wFIZE0QWLVrEjz/+aM1CESBSp5rwtxDk+lYP6wa8rKqvADbhnIkIaWlplCxZks6dO7uOYgKsadOmAKxatcpxkuLl7xQTB0Xkz0BvoJ1vbYLYwMUyJjioKmlpaVx33XWULVvWdRwTYLVq1aJixYp8/fXXrqMUq3NZj+AI0F9VdwE1gWcDlsqYIJGZmcmOHTusWShCiAhJSUlWCAqiqrtUdZSqLvS93qGqEwMbzRj30tLSiI6O5qabbnIdxRSTpKQkVq1axbFjx1xHKTb+jhr6rYhsEpH9InJARA6KyIFAhzPGtbS0NNq1a0flypVdRzHFJCkpiUOHDvHNN9+4jlJs/G0aGgl0VdUKqlpeVcupavlABjPGtQ0bNrB27VprFoowSUlJABHVPORvIditqusCmsSYIJOWlgbAzTff7DaIKVYJCQlER0dHVCHwd9RQhoikAu/h7TQGQFVnBCKUMcFg+vTpXH755dSsWdN1FFOMSpUqxaWXXsrKlZEzwbK/haA8cAi4Lt97ClghMGFp8+bNZGZm8vzzz7uOYhxISkpi0aJFrmMUG78Kgar2DXQQY4JJamoqAD169HCcxLiQlJTE1KlT+fnnn6lUqZLrOAHn76ihRiKyQERW+15fJiJPBDaaMe6kpqbStm1batWq5TqKceB4h3GkNA/521k8FvgzcBRAVVcCvQIVyhiX1q1bx6pVq+jZs6frKMaRSBs55G8hKKOqy056L7eowxgTDFJTUxERbr31VtdRjCPVq1cnLi7OCsFJfhSRBng7iBGRW4EfApbKGEdUldTUVNq1a8dFF13kOo5xJNKmmvC3EPwf8DreRey/AwYC9wcqlDGurFq1ivXr11uzkCEpKYnVq1eTmxv+jR9nHDUkIoPyvZwD/Adv8fgV6A6MOsPXjse7zOUeVW1awOcCvAh0wTs09R5VzTzXEzCmKKWmphIVFUX37t1dRzGOJSUlceTIETZu3Eh8fLzrOAF1tiuCcr6HB3gAqARUxHs10OIsXzsBONME7tcDl/ge9+Fd/MYYZ/Ly8pg0aRLXXnstVatWdR3HONasWTMAvvrqK6c5isMZC4Gq/k1V/4Z32ukWqvpHVX0EaAmccalKVf0cONN6b92Aieq1BKgoItYoa5z55JNP2LFjB3372m0zBpo0aULJkiXJzAz/hgp/+wiqATn5Xuf43iuMGsC3+V7v9L1njBPjx4+nYsWKNreQASA2NpZmzZqRkZHhOkrA+VsIJgLLRGSoiAwFluJt+ikWInKfiGSISMbevXuL67Amgvzyyy+kpaVx5513UqpUKddxTJDweDwsX76cvLw811ECyt+FaZ4C+gI/+x59VXVEIY/9HZD/ts2avvcKOv4YVfWoqqdKlSqFPKwxp3rrrbfIzs6mX79+rqOYIOLxeMjKymLjxo2uowSUv5PO4RvRU5SNZTOB34vINOAKYL+q2r0JptipKq+++iqtW7emRYuzjYEwkcTj8QCQkZFB48aNHacJHH+bhs6ZiEwFFgOXishOEekvIveLyPH7D+YAW4DNeKewSAlUFmPOZMGCBWzcuJGUFPsWNP+rcePGlClTJuz7Cfy+IjhXqnr7WT5XvDeqGePUK6+8QlxcnM00ak4RExND8+bNw74QBOyKwJhQsHHjRt5//30GDBhgncSmQB6PhxUrVoT1HcZWCExEe+GFFyhRogR/+MMfXEcxQcrj8XDo0CHWr1/vOkrAWCEwEWvPnj1MmDCBu+++m2rVCntbjAlX+TuMw5UVAhOxnn32WXJycnjkkUdcRzFBrFGjRpQtW9YKgTHhZteuXbzyyivccccdXHrppa7jmCAWFRVFy5YtrRAYE25GjBhBTk4OTz75pOsoJgR4PB6++uorjh496jpKQFghMBFn586djB49mj59+tCwYUPXcUwI8Hg8HDlyhDVr1riOEhBWCEzEeeqpp1BV/vKXv7iOYkLE8Q7jL7/80nGSwLBCYCLKtm3bGDduHP3796du3bqu45gQ0aBBAy688EKWLFniOkpAWCEwEWXIkCFER0fz+OOPu45iQoiIkJycTHp6uusoAWGFwESMhQsXkpqayuDBg6lZs6brOCbEJCcns379en766UzrbYUmKwQmIuTl5TFw4EBq1KjB4MGDXccxIahNmzYALF261HGSomeFwESECRMmkJmZyciRI7ngggtcxzEhqFWrVkRFRYVl85AVAhP2Dhw4wGOPPUZycjK3337GSXGNOa2yZcty2WWXsXjxYtdRipwVAhP2hg4dyu7du/nnP/+JiLiOY0JYmzZtWLp0KceOHXMdpUhZITBhLTMzkxdffJEBAwZw+eWXu45jQlxycjJZWVlhd2OZFQITto4dO8aAAQOoUqUKI0YUdoltY7yFAAi7fgIrBCZsvfrqq2RkZPDPf/6TSpUquY5jwkD9+vWpWrVq2PUTWCEwYem7777j8ccfp1OnTvTs2dN1HBMmjt9YZoXAmCCnqgwYMIDc3FxeffVV6yA2RSo5OZlNmzaxd+9e11GKjBUCE3beeustZs+ezYgRI6hfv77rOCbMHL+xLJzmHbJCYMLKzp07eeihh7jyyittHWITEB6Ph9jYWBYtWuQ6SpGxQmDChqpy7733kpuby5tvvklUlH17m6JXunRpWrVqxWeffeY6SpGxnxQTNsaOHcu8efMYOXIkDRo0cB3HhLH27duTkZFBVlaW6yhFwgqBCQtr165l4MCBXHvttTzwwAOu45gw1759e3Jzc8PmfgIrBCbkHT58mF69elG2bFkmTpxoTUIm4Nq0aUN0dHTYNA/FuA5gTGH98Y9/ZNWqVXz44YdcdNFFruOYCFCuXDlatmwZNoXA/nQyIS0tLY1XX32VRx55hM6dO7uOYyJI+/btWbZsGYcOHXIdpdCsEJiQtWPHDvr374/H4+Hpp592HcdEmA4dOnD06NGwuJ/ACoEJSbm5udx5550cPXqUqVOnUqJECdeRTIT5zW9+Q1RUVFg0D1kfgQlJf//73/niiy94++23adiwoes4JgKVL1+e5s2b8+mnn7qOUmh2RWBCzmeffcbw4cO5++676d27t+s4JoK1b9+epUuXkp2d7TpKoQS0EIhIZxHZICKbRWRIAZ/fIyJ7ReQr3+PeQOYxoW/fvn3ceeedNGjQgFdeecV1HBPh2rdvz5EjR0J+QfuAFQIRiQZeAa4H4oHbRSS+gE1TVbWZ7/FGoPKY0Keq9OvXjz179jBt2jTKli3rOpKJcFdeeSUiEvLNQ4G8Irgc2KyqW1Q1B5gGdAvg8UyYe+WVV5g5cyYjR46kRYsWruMYQ6VKlWjZsiUfffSR6yiFEshCUAP4Nt/rnb73TtZdRFaKyDsiUqugHYnIfSKSISIZ4TQHuPHf119/zR//+Ee6dOnCQw895DqOMSdcd911LFmyhP3797uOct5cdxZ/ANRV1cuAj4C3CtpIVceoqkdVPVWqVCnWgMa9gwcP0rNnTy688EImTJhgC82YoNKpUyeOHTvGggULXEc5b4EsBN8B+f/Cr+l77wRV3aeqR3wv3wBaBjCPCUGqyn333cemTZuYOnUq9oeACTbJycmUK1eO+fPnu45y3gJZCL4ELhGReiJSAugFzMy/gYjknximK7AugHlMCHr99deZNm0aw4cPp3379q7jGHOK2NhYrr76aubNm4equo5zXgJWCFQ1F/g9MA/vL/h/q+oaERkmIl19mz0oImtE5GvgQeCeQOUxoSczM5OHHnqIzp07M2TIKaOPjQkanTp1Ytu2bWzatMl1lPMioVbBPB6PZmRkuI5hAmz//v20aNGCnJwcVqxYQVxcnOtIxpzWli1baNCgAf/617/4/e9/7zpOgURkuap6CvrMdWexMac4fr/Ajh07SE1NtSJggl79+vVp0KAB8+bNcx3lvFghMEHnueeeY8aMGTzzzDO0adPGdRxj/NKpUyf+85//kJOT4zrKObNCYILK3LlzGTJkCD169GDQoEGu4xjjt06dOvHrr7+G5PKVVghM0Ni0aRO9evUiMTGRN9980+4XMCHlqquuIiYmhrlz57qOcs6sEJigcODAAbp160ZMTAzvvfceF1xwgetIxpyTcuXK0a5dOz744APXUc6ZFQLj3NGjR+nRowcbN25k+vTp1K1b13UkY85Lt27dWLt2LZs3b3Yd5ZxYITBOqSoPPPAA8+fPZ8yYMVx11VWuIxlz3rp1886r+f777ztOcm6sEBinnn76acaNG8cTTzxBv379XMcxplDq1KlDUlIS7733nuso58QKgXFmwoQJPPHEE/Tu3Zthw4a5jmNMkejWrRvp6ens2bPHdRS/WSEwTkyfPp3+/ftz7bXX8sYbb9gIIRM2unfvTl5eHjNmzHAdxW9WCEyxmz17NnfccQfJycm89957lCxZ0nUkY4pMYmIijRs35t///rfrKH6zQmCK1bx58+jevTtJSUnMnj3bhomasCMi3HbbbXz22Wfs2rXLdRy/WCEwxWbSpEnceOONNG7cmHnz5lGhQgXXkYwJiNtuu428vDzeeecd11H8YoXABJyq8txzz3HXXXfRrl07PvvsMypXruw6ljEBk5CQQNOmTZk8ebLrKH6xQmAC6vDhwwwYMIBHH32Unj17MmfOHLsSMBGhT58+LFmyhHXrgn+9LSsEJmCWLVvGFVdcwdixYxkyZAhTpkyxjmETMXr37k10dDRvvVXgUuxBxQqBKXJbtmyhX79+tG7dmn379jFnzhxGjBhBVJR9u5nIUb16da6//nomTpxIbm6u6zhnZD+Zpkjs27ePSZMm0aVLFy655BKmTJnCww8/zLp167j++utdxzPGib59+/LDDz8wa9Ys11HOyJaqNOdl+/btfPrpp6Snp5Oens6aNWtQVWrUqME999xDSkoKF198seuYxjiVm5tLvXr1aNSoEQsWLHCa5UxLVcYUdxgTmo4dO8ann37Kv//9bxYsWMA333wDQMWKFUlOTqZnz55cd911eDweawIyxicmJoaUlBQee+wx1qxZQ0JCgutIBbIrAnNGR44cYfz48Tz33HNs2bKFcuXKcdVVV3HNNddw1VVXkZCQYL/4jTmDH3/8kZo1a3L33XczZswYZzls8XpzXmbNmkVCQgIpKSlUqVKFadOmsXv3bt5//30efPBBEhMTrQgYcxZxcXHcc889TJgwgW+//dZ1nALZT7E5xS+//ELv3r256aabKFGiBB9++CGLFy+mZ8+elC5d2nU8Y0LOn//8ZwCeeeYZx0kKZoXA/I+VK1fi8XhITU1l6NChfPXVV3Tu3NlmBzWmEOrUqcM999zDG2+8wY4dO1zHOYUVAnPC5MmTad26NYcPH+bzzz/nySefpESJEq5jGRMWHn/8caKiohg8eLDrKKewQmDIycnhwQcfpHfv3rRq1Yrly5eTnJzsOpYxYaVOnTr86U9/IjU1lU8//dR1nP9ho4Yi3M6dO+nZsyfp6ekMGjSIZ555htjYWNexjAlLhw8fJj4+npIlS7J8+fJinYbdRg2ZAs2dO5dmzZqxcuVKpk2bxvPPP29FwJgAKl26NG+88QYbN25k0KBBruOcYIUgAh05coTHHnuM66+/nosvvpiMjAx69uzpOpYxEeGaa65h8ODBjBkzhgkTJriOA1ghiDgLFiygefPmjBgxgnvvvZelS5dy6aWXuo5lTEQZPnw4HTt25N577w2KeYisEAQ5VWX//v1s3LiR9evXs23bNnbv3s2RI0f83sfRo0dJS0vj6quv5tprr+XQoUPMnj2bsWPH2n0BxjgQGxvLu+++S7Nmzfjtb3/LxIkTneaxuYaCSF5eHmvXrj0xkduSJUvYvn072dnZBW5fqVIlqlevTvXq1alWrdqJ53FxcRw5coSdO3eyatUqPvvsMw4ePEiNGjUYNWoUDzzwAKVKlSrmszPG5FeuXDk+/vhjunfvTp8+ffjiiy949tlnnSzcFNBRQyLSGXgRiAbeUNVnTvq8JDARaAnsA3qq6rYz7TOcRg3t37+fL7/88sQv/sWLF3PgwAEAqlSpQps2bWjUqBHVqlWjWrVqxMTEkJ2dzeHDh/npp5/YtWvXKY+srKwT+4+JiaFhw4a0b9+eG2+8kc6dOxMTY7XfmGCSk5PDE088wfPPP0/lypV56KGH6N+/P9WrVy/S45xp1FDACoGIRAMbgY7ATuBL4HZVXZtvmxTgMlW9X0R6Abeo6hl7LUOlEKgq2dnZ7N+/nwMHDvDDDz+wdetWtm3bxtq1a8nMzDwxg6eIkJiYSJs2bU486tevf15382ZlZbFv3z5Kly5NpUqVbBSQMSFi+fLl/PWvf2XOnDmICG3atKFdu3bEx8dTr1496taty0UXXXTe83u5KgTJwFBV7eR7/WcAVR2Rb5t5vm0Wi0gMsAuoomcIdb6FYO7cuTz88MP4Mpx4nPza3/fOtM2xY8c4cOBAgasSiQh169alRYsWtGjRAo/HwxVXXGHr+BpjAFi7di3Tpk1j7ty5rFix4n9+jwwcOJAXXnjhvPbraj2CGkD+qfZ2AlecbhtVzRWR/UBl4Mf8G4nIfcB9ALVr1z6vMBUqVCAxMfHEX9kicuJx8mt/3zvdNtHR0ZQvX57y5ctToUIFypcvT9WqValXrx61atWydXuNMacVHx/PsGHDGDZsGNnZ2Wzfvv1Ea0LTpk0DcsyQaDBW1THAGPBeEZzPPpKTk23aBGNMSClVqhSXXnppwId4B3L46HdArXyva/reK3AbX9NQBbydxsYYY4pJIAvBl8AlIlJPREoAvYCZJ20zE+jje34r8MmZ+geMMcYUvYA1Dfna/H8PzMM7fHS8qq4RkWFAhqrOBMYBb4vIZuAnvMXCGGNMMQpoH4GqzgHmnPTeX/M9zwZ6BDKDMcaYM7MpJowxJsJZITDGmAhnhcAYYyKcFQJjjIlwIbdUpYjsBbaf55fHcdJdyxHAzjky2DlHhsKccx1VrVLQByFXCApDRDJON9dGuLJzjgx2zpEhUOdsTUPGGBPhrBAYY0yEi7RCMMZ1AAfsnCODnXNkCMg5R1QfgTHGmFNF2hWBMcaYk1ghMMaYCBeWhUBEOovIBhHZLCJDCvi8pIik+j5fKiJ1HcQsUn6c8yARWSsiK0VkgYjUcZGzKJ3tnPNt111EVERCfqihP+csIrf5/q3XiMiU4s5Y1Pz43q4tIv8RkRW+7+8uLnIWFREZLyJ7RGT1aT4XEXnJ9/9jpYi0KPRBT153N9QfeKe8/gaoD5QAvgbiT9omBRjte94LSHWduxjO+SqgjO/5A5Fwzr7tygGfA0sAj+vcxfDvfAmwAqjke13Vde5iOOcxwAO+5/HANte5C3nO7YAWwOrTfN4F+BAQoDWwtLDHDMcrgsuBzaq6RVVzgGlAt5O26Qa85Xv+DnCNHF9wODSd9ZxV9T+qesj3cgneFeNCmT//zgDDgX8A2cUZLkD8OeffAa+o6s8AqrqnmDMWNX/OWYHyvucVgO+LMV+RU9XP8a7PcjrdgInqtQSoKCIXFeaY4VgIagDf5nu90/degduoai6wH6hcLOkCw59zzq8/3r8oQtlZz9l3yVxLVWcXZ7AA8uffuRHQSEQWicgSEelcbOkCw59zHgr0FpGdeNc/+UPxRHPmXH/ezyokFq83RUdEegMeoL3rLIEkIlHAKOAex1GKWwze5qEOeK/6PheRRFX9xWWoALsdmKCqz4tIMt5VD5uqap7rYKEiHK8IvgNq5Xtd0/degduISAzey8l9xZIuMPw5Z0TkWuBxoKuqHimmbIFytnMuBzQFPhWRbXjbUmeGeIexP//OO4GZqnpUVbcCG/EWhlDlzzn3B/4NoKqLgVJ4J2cLV379vJ+LcCwEXwKXiEg9ESmBtzN45knbzAT6+J7fCnyivl6YEHXWcxaR5sDreItAqLcbw1nOWVX3q2qcqtZV1bp4+0W6qmqGm7hFwp/v7ffwXg0gInF4m4q2FGPGoubPOe8ArgEQkSZ4C8HeYk1ZvGYCd/tGD7UG9qvqD4XZYdg1Dalqroj8HpiHd8TBeFVdIyLDgAxVnQmMw3v5uBlvp0wvd4kLz89zfhYoC0z39YvvUNWuzkIXkp/nHFb8POd5wHUishY4BjyqqiF7tevnOT8CjBWRh/F2HN8Tyn/YichUvMU8ztfv8SQQC6Cqo/H2g3QBNgOHgL6FPmYI//8yxhhTBMKxacgYY8w5sEJgjDERzgqBMcZEOCsExhgT4awQGGNMhLNCYIwxEc4KgTHGRDgrBMb4QURKi8hnIhIdgH2XEJHPfdOdGFPsrBAY459+wAxVPVbUO/ZNr7wA6FnU+zbGH1YIjMlHRKb6Vq9bJiLbReQG30d3Au/n2+4/ItLR9/zvIvKvQh76Pd8xjCl2VgiM+V9JwBZVvRzvL+YnfZOd1VfVbfm2exJ4XETuBJoDAwt53NVAq0Luw5jzYm2SxviISCmgCvA331trgUp4pzT+Jf+2qvq5b1W7QUCHwjYZqeoxEckRkXKqerAw+zLmXFkhMOa/mgKbVPX4spYt8K6Rexjv1MYniEgicBGw7/gvbhGpAUzCO01wa6A33qJSBu/V93JgK9AE7xoBLwH9VHWob7clCY8lNU2IsaYhY/4rCagtIqVE5AK8v8Rf8K3/G+27YsC3PuxkvGvHZuVbDjIJmKKqLwC5wH1AabxXExWAn/EumNMIyMG7oMrrvn1WBn5U1aPFcaLG5GeFwJj/SgJmAEvxLojymqou8n02H/iNiJTxbfOIqq4DhuPtLzj+9Qt9zxVv38EQVR2qqn3wFoLbgdl4rxLK51tQ5Crf+8YUO2saMua/koD7VPXBAj57BXhYVT8Gko+/qaqf53vdENjoWxlsF/AZMEFEvgU+wbuMZBNV/cS3sMrAfPu/AxhSxOdjjF9sYRpjfHyrQdU+3aLnItIPeKuo7yU4vgSjqk4syv0a4y8rBMYYE+Gsj8AYYyKcFQJjjIlwVgiMMSbCWSEwxpgIZ4XAGGMinBUCY4yJcFYIjDEmwv0/M7zVR4UFy8gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "yspace = np.linspace(0, 1, mplot)\n",
    "plt.plot(yspace, density(yspace), color='black')\n",
    "plt.xlabel(r'$p(x_{new})$')\n",
    "plt.ylabel('density')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4104611a",
   "metadata": {},
   "source": [
    "Assumptions: The model may be extrapolated to dosages for which we have no data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa6bd9ab",
   "metadata": {},
   "source": [
    "#### (e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "82cf61aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "predictive = np.zeros(m + 1)\n",
    "\n",
    "for k in range(m + 1):\n",
    "    predictive[k] = math.comb(m, k) * np.mean(pxnew_sample**k * (1 - pxnew_sample)**(m - k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "af5ca89f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEICAYAAABF82P+AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAALEwAACxMBAJqcGAAAFndJREFUeJzt3X+wJWV95/H3h0HFH6gYZncNA8yoYwSMEXMDGozG31goWK6uaLkhhi00iOJSWWtcU6BjKgu6xhiWKKyOiz/RaNgaZZR1BVxdgs6AKAKSjJNBZqLLKEjEH+DAd/84PfHMpefchrl9+8zl/ao6dU4/3X3uZyjqfm8/z9NPp6qQJGm2vYYOIEmaThYISVIrC4QkqZUFQpLUygIhSWq199AB5sv+++9fy5cvHzqGJO1Rrrzyyh9W1dK2fYumQCxfvpwNGzYMHUOS9ihJbtzVPruYJEmtLBCSpFYWCElSKwuEJKmVBUKS1MoCIUlqZYGQJLWyQEiSWlkgJEmter2TOsnRwHuBJcAHqurMWftPA/4DsB3YBvxRVd3Y7LsLuKY59HtVdWyfWaXFYvmqi4aOsJPNZx4zdATdR70ViCRLgHOA5wFbgPVJ1lbVdWOHfQOYqaqfJflj4J3AK5p9P6+qJ/eVT5I0WZ9dTEcAG6tqU1XdCVwAHDd+QFVdWlU/azavAJb1mEeSdC/0WSAOAG4a297StO3KicDnx7b3SbIhyRVJXtJ2QpKTmmM2bNu2bbcDS5J+ZSpWc03yamAGeOZY88FVtTXJY4BLklxTVd8dP6+qzgPOA5iZmakFCyxp3jhmMr36vILYChw4tr2sadtJkucCbwWOrao7drRX1dbmfRNwGXB4j1klSbP0WSDWAyuTrEjyQOB4YO34AUkOB85lVBxuHmvfL8mDms/7A0cB44PbkqSe9dbFVFXbk5wCXMxomuuaqro2yWpgQ1WtBd4FPAz4myTwq+mshwDnJrmbURE7c9bsJ0lSz3odg6iqdcC6WW2nj31+7i7Ouxz4zT6zSZIm805qSVIrC4QkqZUFQpLUygIhSWplgZAktbJASJJaWSAkSa0sEJKkVhYISVIrC4QkqZUFQpLUygIhSWplgZAktbJASJJaWSAkSa0sEJKkVhYISVIrC4QkqZUFQpLUygIhSWplgZAktbJASJJaWSAkSa0sEJKkVhYISVIrC4QkqZUFQpLUygIhSWplgZAktbJASJJaWSAkSa16LRBJjk5yQ5KNSVa17D8tyXVJvpXkS0kOHtt3QpJ/aF4n9JlTknRPvRWIJEuAc4AXAocCr0xy6KzDvgHMVNWTgE8D72zOfRRwBnAkcARwRpL9+soqSbqnPq8gjgA2VtWmqroTuAA4bvyAqrq0qn7WbF4BLGs+vwD4YlXdUlW3Al8Eju4xqyRplj4LxAHATWPbW5q2XTkR+Px9PFeSNM/2HjoAQJJXAzPAM+/leScBJwEcdNBBPSSTpPuvPq8gtgIHjm0va9p2kuS5wFuBY6vqjntzblWdV1UzVTWzdOnSeQsuSeq3QKwHViZZkeSBwPHA2vEDkhwOnMuoONw8tuti4PlJ9msGp5/ftEmSFsicBSLJu5Mcdm+/uKq2A6cw+sV+PfCpqro2yeokxzaHvQt4GPA3Sa5OsrY59xbgHYyKzHpgddMmSVogXcYgrgfOS7I38CHgE1V1W5cvr6p1wLpZbaePfX7uhHPXAGu6/BxJ0vyb8wqiqj5QVUcBfwAsB76V5ONJntV3OEnScDqNQTQ3vT2hef0Q+CZwWpILeswmSRrQnF1MSd4DvAi4BPjzqvp6s+usJDf0GU4a2vJVFw0dYSebzzxm6Ai6H+kyBvEt4E+r6qct+46Y5zySpCnRpYvp1bOLQ5IvAXQdrJYk7Xl2eQWRZB/gIcD+zb0IaXY9HJe9kKRFb1IX02uBNwG/Dlw11v7PwH/rMZMkaQrsskBU1XuB9yZ5Q1WdvYCZJElTYFIX07Or6hJga5KXzt5fVX/bazJJ0qAmdTE9k9HU1he37CvAAiFJi9ikLqYzmvfXLFwcSdK0mNTFdNqkE6vqL+Y/jiRpWkzqYtp3wVJIkqbOpC6mty9kEEnSdJnUxfTmqnpnkrMZDUrvpKre2GsySdKgJnUxXd+8b1iIIJKk6TKpi+mzzfv5AEkePtqsnyxQNknSgLo8cnQmyTWMVnX9dpJvJvnt/qNJkobUZbnvNcDJVfUVgCRPZ/To0Sf1GUySNKwuy33ftaM4AFTVV4Ht/UWSJE2DSbOYntJ8/HKSc4FPMJrN9Argsv6jSZKGNKmL6d2zts8Y+3yPaa+SpMVl0iymZy1kEEnSdOkySE2SY4DDgH12tFXV6r5CSZKG12Wa6/sZjTu8gdFjR18OHNxzLknSwLrMYvrdqvoD4NZmfaanAY/vN5YkaWhdCsTPm/efJfl14JfAo/uLJEmaBl3GID6X5JHAu4CrGM1g+u99hpIkDW/OAlFV72g+fibJ54B9quq2fmNJkoY2Z4FIsg9wMvB0RlcPX03yvqr6Rd/hJEnD6dLF9GHgJ8DZzfargI8wms0kSVqkuhSIJ1bVoWPblya5rq9AkqTp0GUW01VJnrpjI8mRdHyIUJKjk9yQZGOSVS37n5HkqiTbk7xs1r67klzdvNZ2+XmSpPkzabG+axiNOTwAuDzJ95pdBwHfmeuLkywBzgGeB2wB1idZW1XjVx/fA/4Q+JOWr/h5VT25w79BktSDSV1ML9rN7z4C2FhVmwCSXAAcB/xLgaiqzc2+u3fzZ0mS5tkuu5iq6sYdL+CRwIub1yObtrkcANw0tr2laetqnyQbklyR5CVtByQ5qTlmw7Zt2+7FV0uS5tJlLaZTgY8B/6p5fTTJG/oOBhxcVTOMZk39ZZLHzj6gqs6rqpmqmlm6dOkCRJKk+48us5hOBI6sqp8CJDkL+Dt+Ne11V7YCB45tL2vaOqmqrc37piSXAYcD3+16viRp93SZxRTgrrHtu5q2uawHViZZkeSBwPFAp9lISfZL8qDm8/7AUYyNXUiS+tflCuJDwNeSXNhsvwT44FwnVdX2JKcAFwNLgDVVdW2S1cCGqlqb5HeAC4H9gBcneXtVHQYcApzbDF7vBZw5a/aTJKlnEwtEkr2AKxg9g/rpTfNrquobXb68qtYB62a1nT72eT2jrqfZ510O/GaXnyFJ6sfEAlFVdyc5p6oOZ7SSqyTpfqLLGMSXkvzbJF3GHSRJi0SXMYjXAqcBdyXZsYJrVdXD+4slSdNr+aqLho6wk81nHtPL93Z5HsS+vfxkSdJU63IFQZKX8qvnQXylqv5nn6EkScPrcif1XwOvA64Bvg28Lsk5fQeTJA2ryxXEs4FDqqoAkpwPXNtrKknS4LrMYtrIaInvHQ5s2iRJi1iXK4h9geuTfJ3RGMQRwIYdD/GpqmN7zCdJGkiXAnH63IdIkhabLtNcv7wQQSRJ06XLGIQk6X7IAiFJatWpQCR5cJLf6DuMJGl6dLlR7sXA1cAXmu0n75jBJElavLpcQbyN0dTWHwNU1dXAit4SSZKmQpcC8cuqum1WW/URRpI0PbrcB3FtklcBS5KsBN4IXN5vLEnS0LpcQbwBOAy4A/g4cBvwph4zSZKmQJcriCdU1VuBt/YdRpI0PbpcQbw7yfVJ3pHkib0nkiRNhTkLRFU9C3gWsA04N8k1Sf6092SSpEF1ulGuqn5QVX/F6MFBV+MCfpK06HW5Ue6QJG9Lcg1wNqMZTMt6TyZJGlSXQeo1wCeBF1TVP/WcR5I0Jbos9/20hQgiSZouuywQST5VVf+u6Voav3M6QFXVk3pPJ0kazKQriFOb9xctRBAtfstXXTR0hJ1sPvOYoSNIU22Xg9RV9f3m48lVdeP4Czh5YeJJkobSZZrr81raXjjfQSRJ02XSGMQfM7pSeGySb43t2hcX65OkRW/SGMTHgc8D/wVYNdb+k6q6pddUkqTBTRqDuK2qNgPvBW4ZG3/YnuTILl+e5OgkNyTZmGRVy/5nJLkqyfYkL5u174Qk/9C8Trh3/yxJ0u7qMgbxPuD2se3bm7aJkiwBzmE0XnEo8Mokh8467HvAHzK6Whk/91HAGcCRjJ5md0aS/TpklSTNky4FIlX1L/dBVNXddLsD+whgY1Vtqqo7gQuA48YPqKrNVfUt4O5Z574A+GJV3VJVtwJfBI7u8DMlSfOkS4HYlOSNSR7QvE4FNnU47wDgprHtLU1bF53OTXJSkg1JNmzbtq3jV0uSuuhSIF4H/C6wldEv6iOBk/oM1VVVnVdVM1U1s3Tp0qHjSNKi0mUtppuB4+/Dd28FDhzbXta0dT3392ede9l9yCBJuo8m3Qfx5qp6Z5Kz2XktJgCq6o1zfPd6YGWSFYx+4R8PvKpjrouBPx8bmH4+8JaO50qS5sGkK4jrm/cN9+WLq2p7klMY/bJfAqypqmuTrAY2VNXaJL8DXAjsB7w4ydur6rCquiXJOxgVGYDV3nshSQtrlwWiqj7bvJ9/X7+8qtYB62a1nT72eT27ePhQVa1h9CwKSdIAJnUxfZaWrqUdqurYXhJJkqbCpC6m/9q8vxT4N8BHm+1XAv+vz1CSpOFN6mL6MkCSd1fVzNiuzya5T+MSkqQ9R5f7IB6a5DE7NppZSQ/tL5IkaRp0WTLjPwKXJdnE6HGjBwOv7TWVJGlwXW6U+0KSlcATmqbvVNUd/caSJA1tzi6mJA8B/hNwSlV9Ezgoic+plqRFrssYxIeAO4GnNdtbgT/rLZEkaSp0KRCPrap3Ar8EqKqfMRqLkCQtYl0KxJ1JHkxz01ySxwKOQUjSItdlFtMZwBeAA5N8DDiK0VPgJEmL2MQCkWQvRgvpvRR4KqOupVOr6ocLkE2SNKCJBaKq7m6W/f4UcNECZZIkTYEuYxD/O8mfJDkwyaN2vHpPJkkaVJcxiFc0768fayvgMS3HSpIWiS53Uq9YiCCSpOkyZ4FIsg9wMvB0RlcOXwHeX1W/6DmbJGlAXbqYPgz8BDi72X4V8BHg5X2FkiQNr0uBeGJVHTq2fWmS6/oKJEmaDl1mMV2V5Kk7NpIcCfjAIEla5LpcQfw2cHmS7zXbBwE3JLkGqKp6Um/pJEmD6VIgju49hSRp6nSZ5nrjQgSRJE2XLmMQkqT7IQuEJKmVBUKS1MoCIUlqZYGQJLWyQEiSWlkgJEmtLBCSpFYWCElSq14LRJKjk9yQZGOSVS37H5Tkk83+ryVZ3rQvT/LzJFc3r/f3mVOSdE9d1mK6T5IsAc4BngdsAdYnWVtV40uFnwjcWlWPS3I8cBa/esTpd6vqyX3lkyRN1ucVxBHAxqraVFV3AhcAx8065jjg/Obzp4HnJEmPmSRJHfVZIA4Abhrb3tK0tR5TVduB24Bfa/atSPKNJF9O8nttPyDJSUk2JNmwbdu2+U0vSfdz0zpI/X3goKo6HDgN+HiSh88+qKrOq6qZqppZunTpgoeUpMWszwKxFThwbHtZ09Z6TJK9gUcAP6qqO6rqRwBVdSXwXeDxPWaVJM3S2yA1sB5YmWQFo0JwPPCqWcesBU4A/g54GXBJVVWSpcAtVXVXkscAK4FNPWbd4yxfddHQEXay+cxjho4gaZ71ViCqanuSU4CLgSXAmqq6NslqYENVrQU+CHwkyUbgFkZFBOAZwOokvwTuBl5XVbf0lVWSdE99XkFQVeuAdbPaTh/7/Avg5S3nfQb4TJ/ZJEmTTesgtSRpYBYISVIrC4QkqZUFQpLUygIhSWplgZAktbJASJJaWSAkSa0sEJKkVhYISVIrC4QkqZUFQpLUygIhSWplgZAktbJASJJaWSAkSa0sEJKkVhYISVIrC4QkqZUFQpLUygIhSWplgZAktbJASJJaWSAkSa0sEJKkVhYISVKrvYcOMC2Wr7po6Aj/YvOZxwwdQZK8gpAktbNASJJaWSAkSa0sEJKkVhYISVKrXgtEkqOT3JBkY5JVLfsflOSTzf6vJVk+tu8tTfsNSV7QZ05J0j31ViCSLAHOAV4IHAq8Msmhsw47Ebi1qh4HvAc4qzn3UOB44DDgaOCvm++TJC2QPq8gjgA2VtWmqroTuAA4btYxxwHnN58/DTwnSZr2C6rqjqr6R2Bj832SpAXS541yBwA3jW1vAY7c1TFVtT3JbcCvNe1XzDr3gNk/IMlJwEnN5u1Jbpif6Ltlf+CHu/MFOWueknSz23nBzB2YuX97Wl6YjswH72rHHn0ndVWdB5w3dI5xSTZU1czQObra0/KCmRfKnpZ5T8sL05+5zy6mrcCBY9vLmrbWY5LsDTwC+FHHcyVJPeqzQKwHViZZkeSBjAad1846Zi1wQvP5ZcAlVVVN+/HNLKcVwErg6z1mlSTN0lsXUzOmcApwMbAEWFNV1yZZDWyoqrXAB4GPJNkI3MKoiNAc9yngOmA78PqququvrPNsqrq8OtjT8oKZF8qelnlPywtTnjmjP9glSdqZd1JLklpZICRJrSwQ82SuZUWmTZI1SW5O8u2hs3SV5MAklya5Lsm1SU4dOtMkSfZJ8vUk32zyvn3oTF0lWZLkG0k+N3SWLpJsTnJNkquTbBg6TxdJHpnk00m+k+T6JE8bOtNsjkHMg2YZkL8Hnsfopr71wCur6rpBg02Q5BnA7cCHq+qJQ+fpIsmjgUdX1VVJ9gWuBF4yrf+dm1UBHlpVtyd5APBV4NSqumKOUweX5DRgBnh4Vb1o6DxzSbIZmKmq3b7pbKEkOR/4SlV9oJnp+ZCq+vHAsXbiFcT86LKsyFSpqv/DaObYHqOqvl9VVzWffwJcT8sd9tOiRm5vNh/QvKb+L7Iky4BjgA8MnWWxSvII4BmMZnJSVXdOW3EAC8R8aVtWZGp/cS0Gzcq/hwNfGzjKRE1XzdXAzcAXq2qq8zb+EngzcPfAOe6NAv5XkiubJXim3QpgG/ChpivvA0keOnSo2SwQ2uMkeRjwGeBNVfXPQ+eZpKruqqonM1oN4IgkU92dl+RFwM1VdeXQWe6lp1fVUxitHv36pgt1mu0NPAV4X1UdDvwUmLqxSwvE/HBpkAXS9OV/BvhYVf3t0Hm6aroPLmW0fP00Owo4tunTvwB4dpKPDhtpblW1tXm/GbiQ6V/9eQuwZeyK8tOMCsZUsUDMjy7Limg3NYO+HwSur6q/GDrPXJIsTfLI5vODGU1i+M6goeZQVW+pqmVVtZzR/8eXVNWrB441UZKHNpMWaLppng9M9ey8qvoBcFOS32iansNo5Yipskev5jotdrWsyMCxJkryCeD3gf2TbAHOqKoPDptqTkcB/x64punXB/jPVbVuuEgTPRo4v5nlthfwqaraI6aN7mH+NXDh6O8H9gY+XlVfGDZSJ28APtb8UbkJeM3Aee7Baa6SpFZ2MUmSWlkgJEmtLBCSpFYWCElSKwuEJKmVBUKS1MoCIUlqZYGQdkOSJya5fGz7KUm+NGQmab54o5y0G5LsBfwTcEBV3ZXkMuC0HcuSS3syl9qQdkNV3Z3kWuCwJCuBGy0OWiwsENLuu4LROlEn06zWmuQA4KOMFm18KvBq4O3AQxh17V4J/CNwCLAS+Cvgj6rqbQucXdolxyCk3XcF8GfAhTuWnQZ+i9Gice8BtgMnAQ8Gfgw8ArgV2Bd4PHAncCJw7sLGlibzCkLafd8B7gDOGmv7LUbPJYDR084OB15fVXcAJPk94LXAGuBYRk8o/f6CJZY6sEBIu+9U4C1V9dOxtscBf59kf+AHwJeB/5HkJuASRg+MOaSqLkmyGnjTAmeW5uQsJuk+SvJY4CLg/1bViUPnkeabBUKS1MpBaklSKwuEJKmVBUKS1MoCIUlqZYGQJLWyQEiSWlkgJEmt/j+DNZqPZblHYQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar(np.arange(m + 1), predictive)\n",
    "plt.xlabel(r'$y_{new}$')\n",
    "plt.ylabel('predictive probability')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "586186cd",
   "metadata": {},
   "source": [
    "#### (f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "f846ca2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = 100_000\n",
    "sample = np.zeros((sim, 2))\n",
    "\n",
    "for s in range(sim - 1):\n",
    "    if s == 0:\n",
    "        sample[s] = 16 * np.random.random(2) - 8\n",
    "    else:\n",
    "        current = sample[s]\n",
    "        proposal = sample[s] + (2 * np.random.random(2) - 1)\n",
    "        if np.max(proposal) > 8 or proposal[0] < -8 or proposal[1] < 0:\n",
    "            sample[s + 1] = current\n",
    "        else:\n",
    "            log_acceptance_prob = ell(proposal[0], proposal[1], x, y) - ell(current[0], current[1], x, y)\n",
    "            U = np.random.random()\n",
    "            if np.log(U) < log_acceptance_prob:\n",
    "                sample[s + 1] = proposal\n",
    "            else:\n",
    "                sample[s + 1] = current\n",
    "                \n",
    "### burn-in\n",
    "sample = sample[int(.05 * sim):]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "89acfa4a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ESS for a: 5015.97\n",
      "ESS for b: 5536.45\n"
     ]
    }
   ],
   "source": [
    "ess1 = az.ess(sample[:, 0])\n",
    "ess2 = az.ess(sample[:, 1])\n",
    "\n",
    "print(f'ESS for a: {np.round(ess1, 2)}')\n",
    "print(f'ESS for b: {np.round(ess2, 2)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88085538",
   "metadata": {},
   "source": [
    "#### LD50"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "fc5dff05",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posterior 0.05 quantile of LD50: 0.9895\n",
      "Posterior 0.50 quantile of LD50: 1.5448\n",
      "Posterior 0.95 quantile of LD50: 3.6891\n"
     ]
    }
   ],
   "source": [
    "LD50_sample = -sample[:, 0] / sample[:, 1]\n",
    "print(f'Posterior 0.05 quantile of LD50: {np.round(np.quantile(LD50_sample, .05), 4)}')\n",
    "print(f'Posterior 0.50 quantile of LD50: {np.round(np.quantile(LD50_sample, .50), 4)}')\n",
    "print(f'Posterior 0.95 quantile of LD50: {np.round(np.quantile(LD50_sample, .95), 4)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "979f4820",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}