{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chapter 5 - Bayes' Rule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.style.use('seaborn-white')\n",
    "color = '#87ceeb'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# helper function to return an optionally truncated normal distribution\n",
    "def normal(mu, sigma=.75, nvals=11, lbound=0, ubound=1):\n",
    "    dist = []\n",
    "    for x in np.linspace(lbound,ubound,nvals):\n",
    "        dist.append(np.exp(-np.power(x - mu, 2.) / (2 * np.power(sigma, 2.))))\n",
    "    return dist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up\n",
    "Here, we create a vector of values that $\\theta$ can take on.  This is the \"grid\" laid down on the parameter space.  It's a simple grid (only 1 dimension), but that's why what will follow is appropriately called a grid approximation.  By tweaking `n_theta_vals`, this grid can be course-grained (small values of `n_theta_vals`) or it can be very fine-grained (large values of `n_theta_vals`).  Start with small values to get a bit of an understading of how the approximation works and then crank it up to approach an exact solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_theta_vals = 11\n",
    "#n_theta_vals = 1001\n",
    "\n",
    "theta = np.linspace(0, 1, n_theta_vals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prior\n",
    "`p_theta` is the vector of prior probabilities, each associated with the corresponding value of $\\theta$ that we just created.  There are several alternative priors specified here:\n",
    "\n",
    "- Uniform: specifies a belief that all values of $\\theta$ are equally credible\n",
    "- Triangle: values of $\\theta$ near 0.5, are more credible with credibility steadily declining as $\\theta$ moves away from 0.5\n",
    "- (Truncated) Normal: a conventional normal distribution, truncated to the permissible range of $\\theta$ values (i.e., 0-1)\n",
    "\n",
    "Try out each of them.  Tweak the mean and standard deviation of the normal distribution to see how they influence the posterior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# uniform prior\n",
    "p_theta = np.ones_like(theta)\n",
    "\n",
    "# triangle prior\n",
    "#p_theta = np.minimum(theta, 1-theta)\n",
    "\n",
    "# (truncated) normal priors\n",
    "# wide prior\n",
    "#p_theta = normal(.5, 10, n_theta_vals)\n",
    "# narrow prior\n",
    "#p_theta = normal(.5, .04, n_theta_vals)\n",
    "\n",
    "# normalize the vector so that priors sum to 1 (making them proper probabilities)\n",
    "p_theta = p_theta / np.sum(p_theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data\n",
    "This constructs a set of flip outcomes.  Specify the number of heads (i.e., `n_heads`) and the number of tails (i.e., `n_tails`).  There are three scenarios prepared:\n",
    "\n",
    "1. 1 flip that comes up heads\n",
    "2. 4 flips, 1 of which comes up heads (25% heads)\n",
    "3. 40 flips, 10 of which come up heads (25% heads)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# example 1\n",
    "n_heads = 1\n",
    "n_tails = 0\n",
    "\n",
    "# example 2\n",
    "#n_heads = 1\n",
    "#n_tails = 3\n",
    "\n",
    "# example 3\n",
    "#n_heads = 10\n",
    "#n_tails = 30\n",
    "\n",
    "data = np.repeat([1, 0], [n_heads, n_tails])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Likelihood\n",
    "Here, we calculate the likelihood of the data we just constructed (i.e., $p(data|\\theta)$).  We will use a Bernoulli distribution and apply it to each of the values of $\\theta$ we are entertaining.  We then multiply each likelihood by the prior associated with the corresponding value of $\\theta$ and sum all of these values to arrive as what is sometimes called the _evidence_ (i.e., $p(data)$, or the denominator of Bayes' rule)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_data_given_theta = (theta**n_heads) * ( (1-theta)**n_tails )\n",
    "\n",
    "# calculate the evidence (P(D), the prior probability of the data)\n",
    "p_data = np.sum(p_data_given_theta * p_theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inference\n",
    "Here, we put it all together, by applying Bayes' rule to each prior/likelihood pair to generate a corresponding posterior for each value of $\\theta$ we are entertaining."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_theta_given_data = p_data_given_theta * p_theta / p_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize\n",
    "Plot the prior, the likelihood, and the posterior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgEAAAHtCAYAAACAiK0FAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABGc0lEQVR4nO3de1xVdb7/8dcWBFTURNmAgmmmqDimjpdhMBkNvFXH8qSA4O2YTV4wb6PJmNCMaZaXSbPJyrEZtYn0UFPHJi3HOichrCZ10FT0pwSmsFHEiLwA6/eHx30iFWELe7tZ7+fj4ePB2mt91/qs70Ndb9Za3++2GIZhICIiIqbTwNUFiIiIiGsoBIiIiJiUQoCIiIhJKQSIiIiYlEKAiIiISSkEiIiImJRCgIhcV2ZmJqGhofY/3bp149/+7d/YuXPndbdfs2YNoaGh7N2717mFiojDFAJEpEpjxowhNTWVdevW4e3tzRNPPMHx48crbVNRUcGoUaNITU2lY8eONdp/RUVFbZYrIjWgECAiVQoKCqJHjx5EREQwdepULl++THp6OqGhoSxfvpxf/epX/OUvf2HLli3ExMSQnZ0NQG5uLo899hh9+vQhMjKSZ599losXLwLw5JNP0rlzZ1577TV69OjB+fPnXXmKIqalECAi1dawYUMA+8X83XffZeHChURHR1farqysjF//+td8/fXXLF26lDFjxrBhwwZeeukl+zaGYbB//37WrVtH48aNnXcSImKnECAiVaqoqKCsrIzz58/z17/+FQCLxQJAVFQUUVFRtGnTplKbffv2cezYMRISEoiKiuLXv/41d911F9u2bau03eTJkwkPD8fT09M5JyMilehfnohUadWqVaxatQq4cidg9uzZdO3aFQCr1XrdNvn5+QAEBATYP/P392ffvn2VtrtRexFxDoUAEalSQkICDz30EJ6engQHB9O0aVMyMzMBaNDg+jcTr178CwoK7J/l5+cTGBhYabsbtRcR51AIEJEqBQQE8LOf/axGbe655x46dOjA5s2bufvuuzl48CAnTpxg9uzZdVSliDhCMVxEap2npyfr1q2jU6dO/OY3v+E///M/mTJlCv/xH//h6tJE5EcshmEYri5CREREnE93AkRERExKIUBERMSkFAJERERMqt6NDrhw4QJZWVn4+/vj4eHh6nJERETqVHl5OTabjW7duuHj41OjtvUuBGRlZREfH+/qMkRERJxq8+bN9O7du0Zt6l0I8Pf3B650xk8nJhEREalvTp8+TXx8vP36VxP1LgRcfQQQGBhIcHCwi6sRERFxDkcegevFQBEREZOqd3cCrlr/6XGSY93jTsCqD4/Yf54V3cmFlVSPu9ULqtkZ3K1ecL+a3a1eUM3OsP7T4w63rbchYMPuEyTH3uvqMqrlhZ3Z9p/d4S+cu9ULqtkZ3K1ecL+a3a1eUM3OsGH3CbwdbKvHASIiIialECAiImJSCgEiIiImpRAgIiJiUgoBIiIiJqUQICIiYlIKASIiIialECAiImJSTpssaMmSJezbtw+LxUJSUhLdu3e3r0tPT2flypV4eHgwYMAApk2bxpYtW3j33Xft22RlZfHVV18xduxYSktLady4MQDz58+nW7duzjoNERGResMpIWDPnj3k5OSQmprKsWPHSEpKIjU11b5+8eLFrF+/noCAABISEhgyZAijRo1i1KhR9vZ///vf7dsvXbqUTp1u/1mcREREbmdOeRyQkZFBVFQUAB06dKC4uJiSkhIAcnNzad68OUFBQTRo0IDIyEgyMjIqtV+7di1Tp051RqkiIiKm4ZQQUFhYSIsWLezLfn5+2Gw2AGw2G35+ftddB7B//36CgoIqfU/y6tWriY+PZ9GiRVy4cMEJZyAiIlL/uOTFQMMwqr3t1q1befjhh+3L48aNY968eWzevBmLxcLmzZvrokQREZF6zykhwGq1UlhYaF8uKCiw/2b/03X5+flYrVb7cmZmJj179rQvR0dH07ZtWwAGDRrEkSP/95WPIiIiUn1OCQERERFs374dgAMHDmC1WvH19QUgODiYkpIS8vLyKCsrY9euXURERABXAkGTJk3w8vICrtxBmDBhAufPnweuBISOHTs64xRERETqHaeMDujVqxdhYWHExsZisVhITk4mLS2Npk2bEh0dTUpKCnPmzAFg+PDhtG/fHrj2fQGLxcLo0aOZMGECjRo1IiAggMTERGecgoiISL3jtHkC5s6dW2m5c+fO9p/79OlTacjgVd26deO1116r9Nnw4cMZPnx43RQpIiJiIpoxUERExKQUAkRERExKIUBERMSkFAJERERMSiFARETEpBQCRERETEohQERExKQUAkRERExKIUBERMSkFAJERERMSiFARETEpGocAkpLSykvL6+LWkRERMSJbvoFQhUVFWzbto333nuPf/3rX3h5eXHp0iVatGhBZGQksbGx3HnnnTc90JIlS9i3bx8Wi4WkpCS6d+9uX5eens7KlSvx8PBgwIABTJs2jczMTJ544gn7VwV36tSJp556ilOnTjFv3jzKy8vx9/fn+eeft3/VsIiIiFTfTUPAuHHjCA8PZ/bs2XTq1IkGDa7cPDh37hyZmZksX76cqKgoRowYccN97Nmzh5ycHFJTUzl27BhJSUmVvjVw8eLFrF+/noCAABISEhgyZAgAffv2ZfXq1ZX2tXr1asaMGcOwYcNYuXIlW7duZcyYMQ6dvIiIiJndNARs2LCBhg0bXvP5HXfcwZAhQxgyZAiXL1+uch8ZGRlERUUB0KFDB4qLiykpKcHX15fc3FyaN29OUFAQAJGRkWRkZNCpU6fr7iszM5Onn34agIEDB/KnP/1JIUBERMQBN30n4HoBoKbbFBYW0qJFC/uyn58fNpsNAJvNhp+f33XXHT16lMcff5y4uDh2794NwA8//GC//d+yZUv7tiIiIlIzN70TkJeXx+bNm+2/sXfp0oWBAwfSpk0bhw9qGMZNt2nXrh3Tp09n2LBh5ObmMm7cOHbs2FHj/YiIiMj13fROwNSpU7nrrruIj48nPT2dQ4cOkZCQwNNPP82lS5eqdRCr1UphYaF9uaCgAH9//+uuy8/Px2q1EhAQwPDhw7FYLLRt25ZWrVqRn59P48aNuXDhQqVtRUREpOZuGgIqKioYNWoU4eHhNG/enMWLF/Phhx/Spk0bnnrqqWodJCIigu3btwNw4MABrFYrvr6+AAQHB1NSUkJeXh5lZWXs2rWLiIgI3n33XdavXw9ceWRw5swZAgIC+OUvf2nf144dO7j33nsdOnERERGzu+njgPDwcDZt2kRCQgIWi+VKI09PHn30Uftb/DfTq1cvwsLCiI2NxWKxkJycTFpaGk2bNiU6OpqUlBTmzJkDwPDhw2nfvj3+/v7MnTuXnTt3cvnyZVJSUvDy8iIxMZH58+eTmppK69ateeihhxw/exERERO7aQhYsGAB69atY+TIkRQUFJCamoqPjw979+7ljjvuqPaB5s6dW2m5c+fO9p/79OlTacgggK+vLy+//PI1+7FarWzYsKHaxxUREZHru2kIsFgsTJkyhQkTJpCens7XX3/N+fPn6dixI7NmzQKuvKB39S6BiIiIuIdqTRY0ePBg7rvvPvsfgEuXLvHll1/yzjvv0K9fP0aOHFnnxYqIiEjtuWkIeO2119i6dSuzZ88mLy+PZs2acfHiRSoqKoiIiGD8+PF07drVGbWKiIhILbppCPD29iY+Pp74+HguX75MUVERPj4+NGvWzBn1iYiISB256RDBt99+m379+tG3b18WLlxI48aNFQBERETqgZuGgJdeeokNGzbw97//naCgIFatWuWMukRERKSO3TQE+Pr60rVrV1q2bMnMmTPZv3+/M+oSERGROnbTdwJsNhupqancdddddOjQgbKyMmfUJSIiInXspiEgMTGRI0eO8N5773HkyBFKS0uZPHkynTt3JjQ0lAceeMAZdYqIiEgtu2kIiImJqbR8+vRpDh8+zOHDh/nkk08UAkRERNzUTUPATwUGBhIYGEhkZGRd1CMiIiJOctMXA0VERKR+qvGdAEctWbKEffv2YbFYSEpKonv37vZ16enprFy5Eg8PDwYMGMC0adMAeO655/jyyy8pKyvj17/+NYMHD+bJJ5/kwIED9i8vmjRpEr/61a+cdRoiIiL1hlNCwJ49e8jJySE1NZVjx46RlJRU6VsDFy9ezPr16wkICCAhIYEhQ4ZQWFhIdnY2qampFBUV8fDDDzN48GAAZs+ezcCBA51RuoiISL3llBCQkZFBVFQUAB06dKC4uJiSkhJ8fX3Jzc2lefPmBAUFARAZGUlGRgZjxoyx3y1o1qwZP/zwA+Xl5c4oV0RExBSc8k5AYWEhLVq0sC/7+flhs9mAK/MQ+Pn5XbPOw8ODxo0bA7B161YGDBiAh4cHAJs2bWLcuHHMmjWLs2fPOuMURERE6h2XvBhoGEa1t/3oo4/YunUrixYtAmDEiBHMnTuXv/zlL3Tp0oUXX3yxrsoUERGp15wSAqxWK4WFhfblgoIC/P39r7suPz8fq9UKwP/8z//w8ssv8+qrr9K0aVMAwsPD6dKlCwCDBg3iyJEjzjgFERGRescpISAiIoLt27cDcODAAaxWK76+vgAEBwdTUlJCXl4eZWVl7Nq1i4iICL777juee+451q1bZx8JAFdmMMzNzQUgMzOTjh07OuMURERE6h2nvBjYq1cvwsLCiI2NxWKxkJycTFpaGk2bNiU6OpqUlBTmzJkDwPDhw2nfvr19VMDMmTPt+1m2bBnx8fHMnDmTRo0a0bhxY5YuXeqMUxAREal3nDZPwNy5cystd+7c2f5znz59Kg0ZhCvTFf90ymKA1q1b85//+Z91U6SIiIiJaMZAERERk1IIEBERMSmFABEREZNSCBARETEphQARERGTUggQERExKYUAERERk1IIEBERMSmFABEREZNSCBARETEphQARERGTUggQERExKad9gdCSJUvYt28fFouFpKQkunfvbl+Xnp7OypUr8fDwYMCAAUybNu2GbU6dOsW8efMoLy/H39+f559/Hi8vL2edhoiISL3hlDsBe/bsIScnh9TUVJ555hmeeeaZSusXL17MmjVr+Otf/8ru3bs5evToDdusXr2aMWPG8MYbb3DnnXeydetWZ5yCiIhIveOUEJCRkUFUVBQAHTp0oLi4mJKSEgByc3Np3rw5QUFBNGjQgMjISDIyMm7YJjMzk/vuuw+AgQMHkpGR4YxTEBERqXec8jigsLCQsLAw+7Kfnx82mw1fX19sNht+fn6V1uXm5lJUVHTdNj/88IP99n/Lli2x2WzXPebEiHZ1czJ14In7Orq6hBpxt3pBNTuDu9UL7lezu9ULqtkZJka0440djrV12jsBP2YYRq20qWo/k/q3r/ExXGVWdCdXl1Aj7lYvqGZncLd6wf1qdrd6QTU7w6T+7XnDwbZOCQFWq5XCwkL7ckFBAf7+/tddl5+fj9VqpWHDhtdt07hxYy5cuICPj499WxEREak5p4SAiIgI1qxZQ2xsLAcOHMBqteLr6wtAcHAwJSUl5OXlERgYyK5du1i+fDlFRUXXbfPLX/6S7du3M2LECHbs2MG9995b6Vjl5eUAnD592hmnJiIi4lJXr3dXr381YTEcuTfvgOXLl/PFF19gsVhITk7m4MGDNG3alOjoaD7//HOWL18OwODBg5k0adJ123Tu3JmCggLmz5/PxYsXad26NUuXLqVhw4b243zxxRfEx8c745RERERuG5s3b6Z37941auO0EOAsFy5cICsrC39/fzw8PFxdjoiISJ0qLy/HZrPRrVs3fHx8atS23oUAERERqR5NGywiImJSbh0ClixZQkxMDLGxsezfv7/SuvT0dB555BFiYmJYu3atiyqsH6rq588++4zRo0cTGxvLggULqKiocFGV7q2qPr5qxYoVjB071smV1R9V9fGpU6eIi4vjkUceYdGiRS6qsH6oqp83b95MTEwMcXFx18wcK9V35MgRoqKi2LRp0zXranztM9xUZmam8dhjjxmGYRhHjx41Ro8eXWn9sGHDjG+//dYoLy834uLijOzsbFeU6fZu1s/R0dHGqVOnDMMwjMTEROPjjz92eo3u7mZ9bBiGkZ2dbcTExBgJCQnOLq9euFkfz5gxw9ixY4dhGIaRkpJinDx50uk11gdV9fN3331nDBw40Lh8+bJhGIYxceJE46uvvnJFmW7t+++/NxISEoyFCxcaGzduvGZ9Ta99bnsnwJGpiKXmqupngLS0NAIDA4ErszoWFRW5pE53drM+Bnj22WeZNWuWK8qrF6rq44qKCr788ksGDRoEQHJyMq1bt3ZZre6sqn5u2LAhDRs2pLS0lLKyMn744QeaN2/uynLdkpeXF6+++up158hx5NrntiGgsLCQFi1a2JevTisMXHcq4htNLyxVq6qfAft8DwUFBezevZvIyEin1+jubtbHaWlp9O3blzZt2riivHqhqj4+e/YsTZo0YenSpcTFxbFixQpXlen2qupnb29vpk2bRlRUFAMHDuSee+6hfXv3mdn1duHp6XnDEQCOXPvcNgT8lKFBDk5xvX4+c+YMjz/+OMnJyZX+AxDH/LiPz507R1paGhMnTnRhRfXPj/vYMAzy8/MZN24cmzZt4uDBg3z88ceuK64e+XE/l5SUsG7dOj744AN27tzJvn37OHTokAurE3DjEODIVMRSc1X1M1z5hz158mRmzpxJ//79XVGi26uqjz/77DPOnj1LfHw806dP58CBAyxZssRVpbqtqvq4RYsWtG7dmrZt2+Lh4UF4eDjZ2dmuKtWtVdXPx44dIyQkBD8/P7y8vOjduzdZWVmuKrVecuTa57YhICIigu3btwNUORVxWVkZu3btIiIiwpXluq2q+hmuPKseP348AwYMcFWJbq+qPh46dCjvv/8+b731Fi+++CJhYWEkJSW5sly3VFUfe3p6EhISwokTJ+zrdZvaMVX1c5s2bTh27BgXLlwAICsri3bt2rmq1HrJkWufW08W5MhUxFJzN+rn/v3706dPH3r27Gnf9oEHHiAmJsaF1bqnqv4uX5WXl8eCBQvYuHGjCyt1X1X1cU5ODk8++SSGYdCpUydSUlJo0MBtf0dyqar6+c033yQtLQ0PDw969uzJvHnzXF2u28nKymLZsmWcPHkST09PAgICGDRoEMHBwQ5d+9w6BIiIiIjjFHVFRERMSiFARETEpBQCRERETEohQERExKQUAkRERExKIUBERMSkFAJERERMSiFARG5ZeXk5ixcv5v777+fBBx8kNzfX1SWJSDUoBIjILVu3bh0hISFs27aNsWPH8sYbb7i6JBGpBoUAEbklpaWlfPTRR4wfPx64Mn95Tk6Oi6sSkerwdHUBIuLe0tPTOXXqFCNGjACguLiY8PBwF1clItWhECAit+TQoUPMmDGDuLg4AH77298SGhrq4qpEpDr0OEBEbklxcTGNGjUCoKysjN27dzNw4EAXVyUi1aEQICK3pF27duzduxeA119/ncjISEJCQlxblIhUi75KWERuSXFxMZMnT6aoqIgePXrw+9//Hh8fH1eXJSLVoBAgIiJiUnocICIiYlIKASIiIialECAiImJSCgEiIiImpRAgIiJiUgoBIiIiJqUQICIiYlIKASIiIialECAiImJSCgEiIiImpRAgIiJiUgoBIiIiJqUQICIiYlIKASIiIialECAiImJSCgEiIiImpRAgIiJiUgoBIiIiJqUQICIiYlIKASIiIialECAiImJSCgEiIiImpRAgIiJiUgoBIiIiJqUQICIiYlIKASIiIialECAiImJSCgEiIiImpRAgIiJiUgoBIiIiJqUQICIiYlIKASIiIialECAiImJSCgEiIiImpRAgIiJiUgoBIiIiJqUQICIiYlIKASIiIialECAiImJSCgEiIiImpRAgIiJiUgoBIiIiJqUQICIiYlIKASIiIialECAiImJSCgEiIiImpRAgIiJiUgoBIiIiJqUQICIiYlIKASL1VGZmJqGhobzyyivXrBs0aBBDhw4FIC0tjdDQULZt21aj/YeGhjJp0iQA1qxZQ2hoKHv37r1mXV368XmISM15uroAEXG+F198kQYNau93gFGjRnHvvffSsWPHWtuniNQ93QkQMaHp06czc+bMaz6/fPkyY8aM4Ve/+hVnzpyhtLSURYsWER4eTr9+/Xj11Vevu78tW7YQExNDdna2/TNPT0+WLFlCz549eeSRRzh16hQAZWVlvPDCCwwaNIif//znTJgwgWPHjtnbZWRkMHLkSHr27MnQoUPZsmWLfd3+/fu5//776dmzJ6tWraql3hAxL4UAEbFbtmwZBw4cYO3atbRs2ZI1a9bw1ltvMXv2bCZOnMjy5cvZs2dPtfa1b98+rFYrM2bM4F//+hfr168H4NVXX+Wll15i5MiRPP/88xw/fpxf//rXXLp0idzcXB577DGaNGnCmjVr6Nq1KwsXLiQ9PR3DMJg7dy6FhYU8//zzlJSUcPr06brsDpF6T48DRASAbdu2sXPnThYvXkxYWBgA27dvp0uXLjz88MMA/OlPf2L79u307dv3pvtr1aoVjz76KAB//OMfOX78OABvv/027dq1Y/r06QAcPXqUFStWsH//fj7//HMuXbrE7Nmz6dmzJ507d2bbtm3813/9F23atCEnJ4f4+HiioqK49957K90lEJGaUwgQEQB27dqFp6cnBw8etH9WWFjIyZMn7aEA4JtvvqnW/gICAuw/+/r6cunSJQBOnz5Nz5497eusVisABQUF9t/sAwMDgStBokGDBthsNs6ePQtAy5YtAfD29qZZs2Y1Pk8R+T8KASICwJNPPsmJEyd46623SEhIoEOHDvj7++Pr68vixYvt2/n6+t7ScQIDA8nPz7cv//jCf/Xin5+fT1BQEPn5+VRUVBAYGEiLFi0AOHPmDAAXLlzg3Llzt1yPiJkpBIjUc1988QUeHh725UaNGl13u1atWjF8+HDeeecdnnvuOdatW0d0dDSbNm2yv7j35z//mUcffZT27ds7XM/DDz/MH/7wB15++WXuvvtuNm/eTMeOHenevTv+/v689NJLrFq1ismTJ/PWW2/Z27Rt25Y2bdqwbds2fvnLX/Lpp586XIOIXKEQIFLPffLJJ3zyySf25VatWuHt7X3dbf39/Rk3bhwvv/wyGRkZJCYmcv78eZYuXUpFRQXDhg0jKirqlup59NFHuXDhAm+++SYlJSX06tWLp556Ck9PT0JCQnjllVd4/vnnSUxMpE2bNqxatYpevXoBV15cTEpKYv78+UycOJE777yT8vLyW6pHxMwshmEYri5CREREnE9DBEVERExKIUBERMSkFAJERERMqt69GHjhwgWysrLw9/ev9Ea0iIhIfVReXo7NZqNbt274+PjUqG29CwFZWVnEx8e7ugwRERGn2rx5M717965Rm9syBBw5coSpU6cyYcIEEhISKq1LT09n5cqVeHh4MGDAAKZNm1Zpvb+/P3ClM65OPCIiIlJfnT59mvj4ePv1ryZuuxBQWlrK73//e8LDw6+7fvHixaxfv56AgAASEhIYMmQId999t3391UcAgYGBBAcHO6VmERERV3PkEfht92Kgl5cXr776qn0+8R/Lzc2lefPmBAUF0aBBAyIjI8nIyHBBlSIiIu7vtgsBnp6eN3yxwWaz4efnZ1/28/PDZrM5qzQREZHbzvpPjzvc9rYLASIiIlJ9G3afcLitW4UAq9VKYWGhfTk/P/+6jw1ERETk5twqBAQHB1NSUkJeXh5lZWXs2rWLiIgIV5clIiLilm670QFZWVksW7aMkydP4unpyfbt2xk0aBDBwcFER0eTkpLCnDlzABg+fPgtfaWpiIiImd12IaBbt25s3Ljxhuv79OlDamqqEysSERGpn9zqcYCIiIjUntvuToCIiIgrrfrwiP3nWdGdXFhJ3VMIEBER+ZEXdmbbf67vIUCPA0RERExKIUBERMSkFAJERERMSiFARETEpBQCRERETEohQERExKQ0RFBEROqMmcbcuyOFABERqTNmGnPvjvQ4QERExKQUAkRERExKIUBERMSkFAJERERMSiFARETEpBQCRERETEpDBEVE3ITG3Ettq7MQUFpaire3Nx4eHnV1CBERU9GYe6lttRYCKioq2LZtG++99x7/+te/8PLy4tKlS7Ro0YLIyEhiY2O58847a+twIiIicotqLQSMGzeO8PBwZs+eTadOnWjQ4MrrBufOnSMzM5Ply5cTFRXFiBEjauuQIiIicgtqLQRs2LCBhg0bXvP5HXfcwZAhQxgyZAiXL1+urcOJiIjILaq1EJCfn88bb7zBN998Q/PmzenSpQsDBw6kTZs29m2uFxJERETENWptiODUqVNp37498fHxpKenc+jQIRISEnj66ae5dOlSbR1GREREakmtvhg4atQoAJo3b87ixYspKyvj9ddf56mnnmLZsmW1dSgRkVqhIXdidrV2JyA8PJxNmzYBYLFYAPD09OTRRx9l7969tXWYeisvL4+ePXsyduxYEhISGD9+PBkZGQ7v7y9/+QthYWF8//339s8GDRpUaTkzM5MZM2Zc0/bdd9/l3//93xk1ahRbtmwBYM6cOYwdO5ZBgwYxePBgxo4dS0pKisP1idwOXtiZbf8jYka1didgwYIFrFu3jpEjR1JQUEBqaio+Pj7s3buXO+64o7YOU6+1b9+ejRs3AvDNN9/w+OOPs3LlSjp37lyj/bzzzjucOXMGq9Va4xpKS0tZu3YtW7dupWHDhjzyyCNER0ezYsUKANasWUOLFi1ISEio8b5FROT2UmshoEGDBkyZMoUJEyaQnp7O119/zfnz5+nYsSOzZs2qrcOYRtu2bXn88cd54403+N3vfmf/fMuWLbz77ruVtp06dSrh4eH25aioKHx9fXnvvfdqfNx9+/bxs5/9jKZNmwLQq1cv/vnPfzJo0CAHz0RERG5XtT5jYKNGjbjvvvu47777anvXptOtWzfefPPNSp+NGjXK/u7Fjfj6+t5w3eTJk+2zOJ4/f/6aCZwKCwvx8/OzL/v5+WGz2WpauoiIuIFaCwHffvtttbZr1qxZlRcpgCVLlrBv3z4sFgtJSUl0797dvm7QoEEEBgbaL2TLly8nICDA8cJvY99//32tT7v86quv0qRJE+DKOwGbN2+ucnvDMGr1+CIicvuotRAwf/58LBZLlRcNi8XCyJEjeeihh264zZ49e8jJySE1NZVjx46RlJREampqpW1+fCGrz7KysujSpUulz6rzOOBWWK1WCgsL7csFBQX06NGjVvYtIiK3l1oLAVdfaLtVGRkZREVFAdChQweKi4spKSm56d2D+uabb77h9ddfZ8OGDZU+r87jgFtxzz33sHDhQs6fP4+Hhwf//Oc/SUpKqrPjiYiI69T6OwGnT5+muLiYkJAQGjduXOP2hYWFhIWF2ZevPpP+cQhITk7m5MmT/PznP2fOnDn2IYnu7vjx44wdO5ZLly5RXl7OokWLaN26dY3388c//pH09HRsNhuTJ0+mR48ezJs3r8o2r7zyCn369KFnz57MmTOHSZMmYbFYmDZtmv0lQZGqaMy9iPuptRCQl5dHYmIiNpsNHx8fbDYb4eHhzJs3j7vuusvh/f708cKMGTO49957ad68OdOmTWP79u0MHTr0Vst3ueDgYL766qta2deUKVOYMmXKNZ//4x//qLTcr18/+vXrB0BoaKh9ZsehQ4fesE8TExNrpUapf/Q1tyLup9YmC1q+fDkxMTF8+umnfPTRR3z55ZcMHDiQxx57jBMnTlR7P9d7Ju3v729ffuihh2jZsiWenp4MGDCAI0eOXG83UkMNGzbUs38REZOptRBw4sQJYmNj7cuenp7ExMSQkpLC2rVrq72fiIgItm/fDsCBAwewWq32RwHfffcdkyZNsv/G+vnnn9OxY8faOgVT++Uvf4m3t7eryxARESeqtccBN3ou379/f1auXFnt/fTq1YuwsDBiY2OxWCwkJyeTlpZG06ZNiY6OZsCAAcTExODt7U3Xrl3rxaMAERERV6i1EGCz2diyZQudOnXi7rvvrjSEr6Yv7s2dO7fS8o+nzR0/fjzjx4+/tWJFRESk9kJAYmIihw4d4m9/+xvZ2dk0adKEjh070rFjx0rP+EVEROT2UGshYPTo0ZV+4z99+jSHDx/m8OHD9OnTB7jypn99Gc4nUtc05E5E6lqthYBx48YxePBg7rvvPlq3bk1gYCCBgYGEh4fz5ZdfMn/+fPr168fIkSNr65Ai9ZqG3IlIXau1EPDaa6+xdetWZs+eTV5eHs2aNePixYtUVFQQERHB+PHj6dq1a20dTkRERG5RrYUAb29v4uPjiY+P5/LlyxQVFeHj40OzZs1q6xAiIiJSi2ptnoC3336bfv360bdvXxYuXEjjxo0VAERERG5jtRYCXnrpJTZs2MDf//53goKCWLVqVW3tWkREROpArYUAX19funbtSsuWLZk5cyb79++vrV2LiIhIHajVyYJSU1O566676NChA2VlZbW1axEREakDtTpZ0JEjR3jvvfc4cuQIpaWlTJ48mc6dOxMaGsoDDzxQW4cSqTGNuRcRuVathYCYmJhKyz+eLOiTTz5RCBCX0ph7EZFr1VoI+KmrkwVFRkbW1SFERETkFtTai4EiIiLiXhQCRERETEohQERExKQUAkRERExKIUBERMSk6mx0gNRfGnMvIlI/KARIjWnMvYhI/aDHASIiIialECAiImJSCgEiIiImpRAgIiJiUgoBIiIiJqXRAbcBDbkTERFXUAi4DWjInYiIuIIeB4iIiJiUQoCIiIhJKQSIiIiY1G0ZApYsWUJMTAyxsbHs37+/0rr09HQeeeQRYmJiWLt2rYsqFBERcX+3XQjYs2cPOTk5pKam8swzz/DMM89UWr948WLWrFnDX//6V3bv3s3Ro0ddVKmIiIh7u+1CQEZGBlFRUQB06NCB4uJiSkpKAMjNzaV58+YEBQXRoEEDIiMjycjIcGW5IiIibuu2GyJYWFhIWFiYfdnPzw+bzYavry82mw0/P79K63Jzc11RZq164r6Ori6hRtytXlDNzuBu9YL71exu9YJqdoaJEe14Y4djbW+7EPBThmG4uoQ6525zA7hbvaCancHd6gX3q9nd6gXV7AyT+rfnDQfb3naPA6xWK4WFhfblgoIC/P39r7suPz8fq9Xq9BpFRETqg9vuTkBERARr1qwhNjaWAwcOYLVa8fX1BSA4OJiSkhLy8vIIDAxk165dLF++vFL78vJyAE6fPu302kVERJzt6vXu6vWvJm67ENCrVy/CwsKIjY3FYrGQnJxMWloaTZs2JTo6mpSUFObMmQPA8OHDad++faX2NpsNgPj4eKfXLiIi4io2m40777yzRm0sRj176H7hwgWysrLw9/fHw8PD1eWIiIjUqfLycmw2G926dcPHx6dGbetdCBAREZHque1eDBQRERHncOsQoOmFnaOqfv7ss88YPXo0sbGxLFiwgIqKChdV6d6q6uOrVqxYwdixY51cWf1RVR+fOnWKuLg4HnnkERYtWuSiCuuHqvp58+bNxMTEEBcXd81ssFJ9R44cISoqik2bNl2zrsbXPsNNZWZmGo899phhGIZx9OhRY/To0ZXWDxs2zPj222+N8vJyIy4uzsjOznZFmW7vZv0cHR1tnDp1yjAMw0hMTDQ+/vhjp9fo7m7Wx4ZhGNnZ2UZMTIyRkJDg7PLqhZv18YwZM4wdO3YYhmEYKSkpxsmTJ51eY31QVT9/9913xsCBA43Lly8bhmEYEydONL766itXlOnWvv/+eyMhIcFYuHChsXHjxmvW1/Ta57Z3AjS9sHNU1c8AaWlpBAYGAldmcCwqKnJJne7sZn0M8OyzzzJr1ixXlFcvVNXHFRUVfPnllwwaNAiA5ORkWrdu7bJa3VlV/dywYUMaNmxIaWkpZWVl/PDDDzRv3tyV5bolLy8vXn311evOkePItc9tQ0BhYSEtWrSwL1+dXhi47vTCV9dJzVTVz4B9DoeCggJ2795NZGSk02t0dzfr47S0NPr27UubNm1cUV69UFUfnz17liZNmrB06VLi4uJYsWKFq8p0e1X1s7e3N9OmTSMqKoqBAwdyzz33XDPEW27O09PzhiMAHLn2uW0I+ClDgxyc4nr9fObMGR5//HGSk5Mr/QcgjvlxH587d460tDQmTpzoworqnx/3sWEY5OfnM27cODZt2sTBgwf5+OOPXVdcPfLjfi4pKWHdunV88MEH7Ny5k3379nHo0CEXVifgxiFA0ws7R1X9DFf+YU+ePJmZM2fSv39/V5To9qrq488++4yzZ88SHx/P9OnTOXDgAEuWLHFVqW6rqj5u0aIFrVu3pm3btnh4eBAeHk52drarSnVrVfXzsWPHCAkJwc/PDy8vL3r37k1WVparSq2XHLn2uW0IiIiIYPv27QBVTi9cVlbGrl27iIiIcGW5bquqfoYrz6rHjx/PgAEDXFWi26uqj4cOHcr777/PW2+9xYsvvkhYWBhJSUmuLNctVdXHnp6ehISEcOLECft63aZ2TFX93KZNG44dO8aFCxcAyMrKol27dq4qtV5y5Nrn1pMFLV++nC+++MI+vfDBgwft0wt//vnn9u8VGDx4MJMmTXJxte7rRv3cv39/+vTpQ8+ePe3bPvDAA8TExLiwWvdU1d/lq/Ly8liwYAEbN250YaXuq6o+zsnJ4cknn8QwDDp16kRKSgoNGrjt70guVVU/v/nmm6SlpeHh4UHPnj2ZN2+eq8t1O1lZWSxbtoyTJ0/i6elJQEAAgwYNIjg42KFrn1uHABEREXGcoq6IiIhJKQSIiIiYlEKAiIiISSkEiIiImJRCgIiIiEkpBIiIiJiUQoCIiIhJKQSIyC0rLy9n8eLF3H///Tz44IPk5ua6uiQRqQaFABG5ZevWrSMkJIRt27YxduxY3njjDVeXJCLVoBAgIrektLSUjz76iPHjxwNX5i/PyclxcVUiUh2eri5ARNxbeno6p06dYsSIEQAUFxcTHh7u4qpEpDoUAkTklhw6dIgZM2YQFxcHwG9/+1tCQ0NdXJWIVIceB4jILSkuLqZRo0YAlJWVsXv3bgYOHOjiqkSkOhQCROSWtGvXjr179wLw+uuvExkZSUhIiGuLEpFq0VcJi8gtKS4uZvLkyRQVFdGjRw9+//vf4+Pj4+qyRKQaFAJERERMSo8DRERETEohQERExKQUAkRERExKIUBERMSkFAJERERMSiFARETEpBQCRERETEohQERExKQUAkRERExKIUBERMSkFAJERERMSiFARETEpBQCRERETEohQERExKQUAkRERExKIUBERMSkFAJERERMSiFARETEpBQCRERETEohQERExKQUAkRERExKIUBERMSkFAJERERMSiFARETEpBQCRERETEohQERExKQUAkRERExKIUBERMSkFAJERERMSiFARETEpBQCRERETEohQERExKQUAkRERExKIUBERMSkFAJERERMSiFARETEpBQCRERETEohQERExKQUAkRERExKIUBERMSkFAJERERMSiFARETEpBQCRERETEohQERExKQUAkRERExKIUBERMSkFAJERERMSiFARETEpBQCREwgMzOT0NBQ+59u3brxb//2b+zcufOW9vv222/z9ddfO9R20KBBDB069JaOLyK3RiFAxETGjBlDamoq69atw9vbmyeeeILjx487tK/S0lKefvrpGoeAiooKAF588UX+8Ic/1Pi45eXlNW4jItenECBiIkFBQfTo0YOIiAimTp3K5cuXSU9PJyMjg5EjR9KzZ0+GDh3Kli1b7G0+/PBDhg8fTvfu3Rk8eDDvvfceAD179uSHH35gwYIFPPnkkwD813/9F0OHDqV79+5MnjyZc+fOAfDkk0/SuXNnXnvtNXr06MH58+eZPn06M2fOtB/n/fff54EHHqBnz56MGDGCf/zjH/Z1oaGhzJo1i7i4OGbMmFH3HSViEgoBIibVsGFDAC5evMhjjz1GkyZNWLNmDV27dmXhwoWkp6dz8eJFZs+eTZcuXfjTn/5Enz59WLx4MUVFRSxatAiAKVOmMHXqVA4dOsTcuXPp1KkTq1ev5uuvv2bZsmX24xmGwf79+1m3bh2NGzeuVMsXX3zBrFmzuOuuu1izZg0tWrQgMTGRY8eO2bf5xz/+wbBhw0hMTHRC74iYg0KAiIlUVFRQVlbG+fPn+etf/wpcCQGXLl1i9uzZ9O/fn6SkJODKb/WXLl2ioqKCkydPcvbsWX7zm9+QmZlJixYtuPvuuwFo27Ytbdu25cMPP8QwDKZMmUL//v0ZNmwY27dvxzAM+/EnT55MeHg4np6elep65513AFi4cCH9+/dn1qxZlJWV8cEHH9i3CQkJYdy4cXTu3Lkuu0jEVDxvvomI1BerVq1i1apVwJU7AbNnz+bbb78FIDAwEIBWrVrRoEEDbDYbTZs25fnnn2f58uUkJibi6enJiBEj+N3vfnfNvgsLCwF46KGHKn1+9uxZ+89Wq/W6dZ0+fZoGDRrY1/v7+wNgs9lu2lZEHKcQIGIiCQkJPPTQQ3h6ehIcHEzTpk354x//CEB+fj5BQUHk5+dTUVFhDwVDhgxh2LBh5OTksGXLFl577TXuv//+a36bv3rhfvHFF+1tAZo2bWr/uUGD6998DAwMpKKiApvNhr+/P6dPn7Z/fpXFYqmFHhCRH9PjABETCQgI4Gc/+xldunSxX5wfeOABvLy8WLVqFZ9++inPPPMMAA8//DCHDx+mR48evPDCCxQWFtKkSRMAvL298fHxAeCTTz7h4MGDREVFYbFY+Mc//sH58+dZu3Ytf/7zn/Hy8rppXQ8//DAAS5Ys4dNPP2XlypV4e3tz//3310U3iMj/UggQMbmQkBBeeeUVvvvuOxITE/l//+//sWrVKnr16kVoaCgzZ87kb3/7GxMmTGDr1q385je/oXfv3nTp0oVevXrx3//937z11lt07tyZZ599ln/+859MmTKFoqIiHn300WrV8POf/5wVK1Zw+PBhEhMTKS0tZd26dYSEhNTx2YuYm8X48Vs7IiIiYhq6EyAiImJSCgEiIiImpRAgIiJiUgoBIiIiJlXv5gm4cOECWVlZ+Pv74+Hh4epyRERE6lR5eTk2m41u3brZh+5WV70LAVlZWcTHx7u6DBEREafavHkzvXv3rlGbehcCrs5atnnz5kqzjYmIiNRHp0+fJj4+3n79q4l6FwKuPgIIDAwkODjYxdWIiIg4hyOPwPVioIiIiBtb/+lxh9sqBIiIiLixDbtPONxWIUBERMSkFAJERERMSiFARETEpBQCRERETEohQERExKQUAkREREyq3k0WJCIicitWfXjE/vOs6E4urKTuKQSIiIj8yAs7s+0/1/cQoMcBIiIiJuXSOwFLlixh3759WCwWkpKS6N69u33dxYsXWbRoEdnZ2aSlpQGQmZnJE088QceOHQHo1KkTTz31lEtqFxERcXcuCwF79uwhJyeH1NRUjh07RlJSEqmpqfb1zz33HF26dCE7O7tSu759+7J69WpnlysiIlLvuOxxQEZGBlFRUQB06NCB4uJiSkpK7OtnzZplXy8iIiK1z2UhoLCwkBYtWtiX/fz8sNls9mVfX9/rtjt69CiPP/44cXFx7N69u87rFBERqa9um9EBhmHcdJt27doxffp0hg0bRm5uLuPGjWPHjh14eXk5oUIREZH6xWUhwGq1UlhYaF8uKCjA39+/yjYBAQEMHz4cgLZt29KqVSvy8/MJCQmp01pFRMQxZhpz745c9jggIiKC7du3A3DgwAGsVusNHwFc9e6777J+/XoAbDYbZ86cISAgoM5rFRERx7ywM9v+R24/LrsT0KtXL8LCwoiNjcVisZCcnExaWhpNmzYlOjqaGTNmcPr0aY4fP87YsWMZPXo0gwYNYu7cuezcuZPLly+TkpKiRwEiIiIOcuk7AXPnzq203LlzZ/vPNxoG+PLLL9dpTSIiImahGQNFRERMSiFARETEpBQCRERETOq2mSdARESqpuF2UttqJQSUlpbi7e2Nh4dHbexORESuw0xfcSvO4VAIqKioYNu2bbz33nv861//wsvLi0uXLtGiRQsiIyOJjY3lzjvvrO1aRUREpBY5FALGjRtHeHg4s2fPplOnTjRocOXVgnPnzpGZmcny5cuJiopixIgRtVqsiIiI1B6HQsCGDRto2LDhNZ/fcccdDBkyhCFDhnD58uVbLk5ERETqjkOjA64GgNOnT3P48GFKS0tvuI2IiIjcnhy6E5CXl0diYiKFhYV4e3tTWFjIL37xC+bPn0/79u1ru0YRERGpAw7dCVi+fDkxMTH8z//8Dx999BFffPEFAwcOZPLkyZw4caKWSxQREZG64NCdgBMnTvCHP/zh/3bi6UlMTAxt2rRh7dq1PP/889Xaz5IlS9i3bx8Wi4WkpCS6d+9uX3fx4kUWLVpEdnY2aWlp1WrjzvLy8njwwQfp1q0bhmFw6dIlJk+eTHR0NACHDh1i7dq1rFmzhrCwMHr16oVhGBiGQXx8PMOHD+fAgQO88sorvPDCCw7V8Nprr/HBBx9gsViYPn06kZGRldaPHTuW0tJSGjduDMD8+fPp1q3brZ24iAtp3L2YnUMhwGKxXPfz/v37s3LlymrtY8+ePeTk5JCamsqxY8dISkoiNTXVvv65556jS5cuZGdnV7uNu2vfvj0bN24Eroy0ePjhh7n33nvx8fEhOTmZVatWAeDr62vfrrCwkKlTp+Lr68uAAQPw9/fngw8+YOjQoTU6dm5uLu+//z5vvvkmJSUljBkzhv79+18z98PSpUvp1En/WUr9oHH3YnYOPQ6w2Wxs2bKFffv28f3331dad6OA8FMZGRlERUUB0KFDB4qLiykpKbGvnzVrln19ddvUJ3fccQf+/v7YbDa++OILWrZsSevWra/ZrlWrVsyfP58///nPwJXf1q/+fFV5eTljx46t9GfevHmVtsnMzOTee+/Fy8sLPz8/2rRpw9GjR+vuBEVExOUcuhOQmJjIoUOH+Nvf/kZ2djZNmjShY8eOdOzYkcLCwmrto7CwkLCwMPuyn58fNpsNX19f4Mpvu+fOnatRm/okLy+Pc+fOERQUxN/+9jf69Olzw21/9rOf2S/Yd955J6dOneKHH36gUaNGAHh4eNjvHNxIYWEhfn5+9uWrfRsaGlppu9WrV1NUVESHDh1ISkrCx8fH0VMUEREXcygExMTEVFq+OlTw8OHDVV6sqmIYhlPa3M6OHz/O2LFjMQwDb29vli1bhqenJwUFBfziF7+4YbuSkpJKt+1btWpFYWEhISEhDtdyvb4dN24coaGhtG3bluTkZDZv3sykSZMcPoaIiLiWQyHg22+/veazq3cCfry+WbNmN/wt3Wq1VrprUFBQgL+/f5XHdaSNO/nxOwE/VdVjlqysLLp06XLD9eXl5UyYMKHSZ0FBQTz33HP2ZavVyvHjx+3L+fn5WK3WSm2uvqQIMGjQIN5///0bHlNERG5/DoWA+fPnX/fzqxcqwzCwWCyMHDmShx566LrbRkREsGbNGmJjYzlw4ABWq/Wmt/UdaVMfWK1W8vPzr7vuzJkzrFy5kt/97neVPmvVqpV9uTqPA37xi1+wYcMGEhMTKSoqoqCggLvvvtu+3jAMJk6cyOrVq2nWrBmZmZn20CciIu7JoRBwswtKdfTq1YuwsDBiY2OxWCwkJyeTlpZG06ZNiY6OZsaMGZw+fdp+i3z06NE8+OCD17Qxg1/84he8/vrr9t/mS0pKGDt2LJcvX+bChQv8x3/8h32o5DfffENAQID9fYDqat26NaNHjyYhIQGLxUJKSgoNGjTgv//7v8nLy2PMmDGMHj2aCRMm0KhRIwICAkhMTKztUxU3puF2Iu6nVr5K2FFz586ttNy5c2f7z6tXr65Wm/oiODi40nwIP9a7d2+ee+45Tp06RVBQEAcOHLjhfjZt2sS4ceMcquHqyIEfGzBggP3n4cOHM3z4cIf2LfWfhtuJuB+Hpw3evHkzubm5NG/enC5dujBw4EDatGlT2/XJ/3r66adZunTpDcMRwNdff83p06d1oRYRkWpxaJ6AqVOnctdddxEfH096ejqHDh0iISGBp59+mkuXLtV2jQJ06dKlygBQ3W1ERESucigEVFRUMGrUKMLDw2nevDmLFy/mww8/pE2bNjz11FO1XaOIiIjUAYdCQHh4OJs2bQL+b0SAp6cnjz76KHv37q214kRERKTuOPROwIIFC1i3bh0jR46koKCA1NRUfHx82Lt3L3fccUctlygiIiJ1weEvEJoyZQoTJkwgPT2dr7/+mvPnz9OxY0dmzZoF/N9cASIiInJ7cigEjBs3jsGDB3PffffZ/wBcunSJL7/8knfeeYd+/foxcuTIWi1WxEw07l5E6ppDIeC1115j69atzJ49m7y8PJo1a8bFixepqKggIiKC8ePH07Vr19quVcRUNO5eROqaQyHA29ub+Ph44uPjuXz5MkVFRfj4+NCsWbPark9ERETqiEOjA95++2369etH3759WbhwIY0bN1YAEBERcTMOhYCXXnqJDRs28Pe//52goCBWrVpV23WJiIhIHXMoBPj6+tK1a1datmzJzJkz2b9/f23XJSIiInXMoXcCbDYbqamp3HXXXXTo0IGysjKHDr5kyRL27duHxWIhKSnJ/k14AOnp6axcuRIPDw8GDBjAtGnTyMzM5IknnrB/hW2nTp00Q6GIiIiDHAoBiYmJHDlyhPfee48jR45QWlrK5MmT6dy5M6GhoTzwwAM33ceePXvIyckhNTWVY8eOkZSURGpqqn394sWLWb9+PQEBASQkJDBkyBAA+vbtq/nxRUREaoFDISAmJqbS8unTpzl8+DCHDx/mk08+qVYIyMjIICoqCoAOHTpQXFxMSUkJvr6+9m8nDAoKAiAyMpKMjAw6ddIwKXGMxtyLiFzLoRDwU4GBgQQGBhIZGVntNoWFhYSFhdmX/fz8sNls+Pr6YrPZ8PPzq7QuNzeXTp06cfToUR5//HGKi4uZPn06ERERtXEKUs9pzL2IyLVqJQTUBsMwbrpNu3btmD59OsOGDSM3N5dx48axY8cOvLy8nFChiIhI/eLQ6IDaYLVaKSwstC8XFBTg7+9/3XX5+flYrVYCAgIYPnw4FouFtm3b0qpVK/Lz851eu4iISH3gshAQERHB9u3bAThw4ABWqxVfX18AgoODKSkpIS8vj7KyMnbt2kVERATvvvsu69evB66MUDhz5gwBAQGuOgURERG35rLHAb169SIsLIzY2FgsFgvJycmkpaXRtGlToqOjSUlJYc6cOQAMHz6c9u3b4+/vz9y5c9m5cyeXL18mJSVFjwJEREQc5NJ3AubOnVtpuXPnzvaf+/TpU2nIIFyZpOjll192Sm0iIiL13W3zYqC4Dw23ExGpHxQCpMY03E5EpH5w2YuBIiIi4loKASIiIialECAiImJSCgEiIiImpRAgIiJiUgoBIiIiJqUhgrcBjbsXERFXUAi4DWjcvYiIuIIeB4iIiJiUQoCIiIhJufRxwJIlS9i3bx8Wi4WkpCS6d+9uX5eens7KlSvx8PBgwIABTJs27aZtREREpPpcFgL27NlDTk4OqampHDt2jKSkpErfGrh48WLWr19PQEAACQkJDBkyhLNnz1bZRkRERKrPZSEgIyODqKgoADp06EBxcTElJSX4+vqSm5tL8+bNCQoKAiAyMpKMjAzOnj17wzYiIiJSMy4LAYWFhYSFhdmX/fz8sNls+Pr6YrPZ8PPzq7QuNzeXoqKiG7ZxZ0/c19HVJdSIu9ULqtkZ3K1ecL+a3a1eUM3OMDGiHW/scKztbTNE0DAMp7S5HbnbsEB3qxdUszO4W73gfjW7W72gmp1hUv/2vOFgW5eFAKvVSmFhoX25oKAAf3//667Lz8/HarXSsGHDG7YRERGRmnFZCIiIiGDNmjXExsZy4MABrFar/bZ+cHAwJSUl5OXlERgYyK5du1i+fDlFRUU3bHNVeXk5AKdPn3b6OYmIiDjb1evd1etfTbgsBPTq1YuwsDBiY2OxWCwkJyeTlpZG06ZNiY6OJiUlhTlz5gAwfPhw2rdvT/v27a9p81M2mw2A+Ph4p56PiIiIK9lsNu68884atbEY9eXB+v+6cOECWVlZ+Pv74+Hh4epyRERE6lR5eTk2m41u3brh4+NTo7b1LgSIiIhI9WjaYBEREZNSCBARETEptw4BS5YsISYmhtjYWPbv319pXXp6Oo888ggxMTGsXbvWRRXWD1X182effcbo0aOJjY1lwYIFVFRUuKhK91ZVH1+1YsUKxo4d6+TK6o+q+vjUqVPExcXxyCOPsGjRIhdVWD9U1c+bN28mJiaGuLg4nnnmGRdV6P6OHDlCVFQUmzZtumZdja99hpvKzMw0HnvsMcMwDOPo0aPG6NGjK60fNmyY8e233xrl5eVGXFyckZ2d7Yoy3d7N+jk6Oto4deqUYRiGkZiYaHz88cdOr9Hd3ayPDcMwsrOzjZiYGCMhIcHZ5dULN+vjGTNmGDt27DAMwzBSUlKMkydPOr3G+qCqfv7uu++MgQMHGpcvXzYMwzAmTpxofPXVV64o0619//33RkJCgrFw4UJj48aN16yv6bXPbe8E3Oi7B4BK3z3QoEED+3cPSM1V1c8AaWlpBAYGAlemcS4qKnJJne7sZn0M8OyzzzJr1ixXlFcvVNXHFRUVfPnllwwaNAiA5ORkWrdu7bJa3VlV/dywYUMaNmxIaWkpZWVl/PDDDzRv3tyV5bolLy8vXn31VaxW6zXrHLn2uW0IKCwspEWLFvblq98jAFz3uweurpOaqaqfAftkTQUFBezevZvIyEin1+jubtbHaWlp9O3blzZt2riivHqhqj4+e/YsTZo0YenSpcTFxbFixQpXlen2qupnb29vpk2bRlRUFAMHDuSee+6hffv2rirVbXl6et5wGKAj1z63DQE/ZWiko1Ncr5/PnDnD448/TnJycqX/AMQxP+7jc+fOkZaWxsSJE11YUf3z4z42DIP8/HzGjRvHpk2bOHjwIB9//LHriqtHftzPJSUlrFu3jg8++ICdO3eyb98+Dh065MLqBNw4BDjy3QNSc1X1M1z5hz158mRmzpxJ//79XVGi26uqjz/77DPOnj1LfHw806dP58CBAyxZssRVpbqtqvq4RYsWtG7dmrZt2+Lh4UF4eDjZ2dmuKtWtVdXPx44dIyQkBD8/P7y8vOjduzdZWVmuKrVecuTa57YhICIigu3btwNU+d0DZWVl7Nq1i4iICFeW67aq6me48qx6/PjxDBgwwFUlur2q+njo0KG8//77vPXWW7z44ouEhYWRlJTkynLdUlV97OnpSUhICCdOnLCv121qx1TVz23atOHYsWNcuHABgKysLNq1a+eqUuslR659bj1j4PLly/niiy/s3yNw8OBB+3cPfP755yxfvhyAwYMHM2nSJBdX675u1M/9+/enT58+9OzZ077tAw88QExMjAurdU9V/V2+Ki8vjwULFrBx40YXVuq+qurjnJwcnnzySQzDoFOnTqSkpNCggdv+juRSVfXzm2++SVpaGh4eHvTs2ZN58+a5uly3k5WVxbJlyzh58iSenp4EBAQwaNAggoODHbr2uXUIEBEREccp6oqIiJiUQoCIiIhJKQSIiIiYlEKAiIiISSkEiIiImJRCgIiIiEkpBIiIiJiUQoCI3LLy8nIWL17M/fffz4MPPkhubq6rSxKRalAIEJFbtm7dOkJCQti2bRtjx47ljTfecHVJIlINCgEicktKS0v56KOPGD9+PHBl/vKcnBwXVyUi1eHp6gJExL2lp6dz6tQpRowYAUBxcTHh4eEurkpEqkMhQERuyaFDh5gxYwZxcXEA/Pa3vyU0NNTFVYlIdehxgIjckuLiYho1agRAWVkZu3fvZuDAgS6uSkSqQyFARG5Ju3bt2Lt3LwCvv/46kZGRhISEuLYoEakWfZWwiNyS4uJiJk+eTFFRET169OD3v/89Pj4+ri5LRKpBIUBERMSk9DhARETEpBQCRERETEohQERExKQUAkRERExKIUBERMSkFAJERERMSiFARETEpP4/gy7C6bu+N6AAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 576x576 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(8, 8))\n",
    "plt.subplots_adjust(hspace = 1.1)\n",
    "\n",
    "# Plot the prior, the likelihood, and the posterior:\n",
    "for i,dist in enumerate([p_theta, p_data_given_theta, p_theta_given_data]):\n",
    "    plt.subplot(3, 1, i+1)\n",
    "    markerline, stemlines, baseline = plt.stem(theta, dist, markerfmt=' ', basefmt=' ', use_line_collection=True)\n",
    "    plt.setp(stemlines, 'linewidth', 3)\n",
    "    plt.xlim(0, 1)\n",
    "    plt.xlabel('$\\\\theta$')\n",
    "\n",
    "# prior\n",
    "plt.axes(fig.axes[0])\n",
    "plt.xlabel(r'$\\theta$')\n",
    "plt.ylabel('$P(\\\\theta)$')\n",
    "plt.title('Prior', weight='bold')\n",
    "\n",
    "# likelihood\n",
    "plt.axes(fig.axes[1])\n",
    "plt.xlabel(r'$\\theta$')\n",
    "plt.ylabel('$P(D|\\\\theta)$')\n",
    "plt.title('Likelihood', weight='bold')\n",
    "plt.text(0.1, np.max(p_data_given_theta)/2, 'D = %sH,%sT' % (n_heads, n_tails))\n",
    "\n",
    "# posterior\n",
    "plt.axes(fig.axes[2])\n",
    "plt.xlabel(r'$\\theta$')\n",
    "plt.ylabel('$P(\\\\theta|D)$')\n",
    "plt.title('Posterior', weight='bold')\n",
    "plt.text(0.1, np.max(p_theta_given_data)/2, 'P(D) = %g' % p_data)\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}