Created
August 31, 2022 17:18
-
-
Save ljbelenky/db9e0eac06552c3d1fb1eaa46cfe4912 to your computer and use it in GitHub Desktop.
Bayesian Approach to Determining P of a Bernoulli Trial.ipynb
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "id": "8f432b21", | |
| "metadata": {}, | |
| "source": [ | |
| "# Bayesian Approach to Determining $P$ of a Bernoulli Trial\n", | |
| "\n", | |
| "\n", | |
| "Consider the following experiment: We have a coin that is not necessarily fair. We flip it 20 times and obtain 13 heads and 7 tails. We wish to determine the value $p$ of this coin, that is, the long-term probability of heads.\n", | |
| "\n", | |
| "Considering the limited data available, the most obvious conclusion is:\n", | |
| "$$ p= \\frac{Heads}{Trials} = \\frac{13}{20} = 0.65 $$\n", | |
| "\n", | |
| "This suggests that the coin has a small bias towards heads.\n", | |
| "\n", | |
| "While this value provides the best possible (point) estimate of $p$, it is not the only possible value. Certainly, if the true value for $p$ was 0.64, the result obtained would not be surprising. Similarly, these same results could be obtained from a coin with $p=0.8$ that had a bit of bad luck, or from a fair coin ($p=.5$) that had a bit of good luck, or even from a coin biased towards tails ($p=0.48$) that had a run of good luck.\n", | |
| "\n", | |
| "Therefore, although the coin does have a fixed value for $p$, this true value is hidden from us. The data available provides us only the abilty to describe $p$ as a random variable within a given range. I.e, we can say $p$ has a maximum likelihood at 0.65, but may be as low as $\\large p_{lower}$ or as high as $\\large p_{upper}$\n", | |
| "\n", | |
| "To asses this distribution, we start by considering the binomial distribution which describe the probability of an observered outcome given a value $p$.\n", | |
| "\n", | |
| "$$P_x = \\mathrm{C}_{n}^{x}p^x(1-p)^{n-x}$$\n", | |
| "\n", | |
| "Where $n$ is the number of trials (20), $x$ is the number of heads (13), $p$ is the fixed probability of a positive outcome, and $C$ is the combinations function:\n", | |
| "\n", | |
| "$$ \\mathrm{C}_{n}^{x} = \\frac{n!}{x!(n-x)!}$$" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 1, | |
| "id": "5f97a8d7", | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "from math import factorial as fac\n", | |
| "import numpy as np\n", | |
| "from scipy.stats import beta\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "\n", | |
| "C = lambda n,r: fac(n)/(fac(r)*fac(n-r))\n", | |
| "P = lambda n,x,p: C(n,x)*p**(x)*(1-p)**(n-x)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "id": "521d6ccd", | |
| "metadata": {}, | |
| "source": [ | |
| "Using this equation, we can determine the probability of obtaining exactly 13 heads in 20 flips from a coin with $p=0.65$" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 2, | |
| "id": "aca39f64", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "0.18440118638393127" | |
| ] | |
| }, | |
| "execution_count": 2, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "P(20,13, .65)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "id": "70cb81be", | |
| "metadata": {}, | |
| "source": [ | |
| "The first result, which may be a bit surprising, is that it's actually fairly unlikely (less than 1 out of five) that a coin performs exactly as expected. While 13 heads is the most probable outcome from this experiment, it is actually more probable that the result would be something other than 13 (e.g., 11, 12, 14 or 15)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 3, | |
| "id": "33be8492", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "0.5756216677999606" | |
| ] | |
| }, | |
| "execution_count": 3, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "P(20,11, .65) + P(20,12, .65) + P(20,14, .65) + P(20,15, .65) " | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "id": "490c5038", | |
| "metadata": {}, | |
| "source": [ | |
| "Similarly, we can assess the likelihood that the probabilty $p$ is any value in the range of 0 to 1 given the observed results." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 4, | |
| "id": "5f1e8f76", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAs40lEQVR4nO3deXxU93nv8c8zM1oQQmgFBAgkQGwOi0FgbGOHJV7wEuLYybWT2K6bXkpj13HTpHHT3ja9N01dN6tT146TOnHaNK7jJSExAdvg3dggMJvYJMQmkNACaEXrPPePGbXTiUBHYqQzy/N+vfTSzFlmnvNCzHfO7/x+5yeqijHGmMTjcbsAY4wx7rAAMMaYBGUBYIwxCcoCwBhjEpQFgDHGJCif2wUMRG5urhYWFrpdhjFx4WDDQQBm5MxwuRIz1LZv316vqnnhy2MqAAoLCyktLXW7DGPiwrKfLgPgjT94w9U6zNATkWN9LbcmIGOMSVAWAMYYk6AsAIwxJkFZABhjTIKyADDGmARlAWCMMQnKAsAYYxKUBYAxCarHr9S3dFB69IzbpRiXWAAYk2D2nmzkay/tYcfxs1TUtvCpH27hO68eosdvc4MkmpgaCWyMGZyWjm7W7TzFL7YeZ8/JRlKTPGRnJJM3KoUVORN5bFM5O46d5Xt3zic3PcXtcs0wsTMAY+JcTWM7y/7pDb720h46u/383ccv44OvfYypeelkpCbxrU/N49Hb57Lt6BlueewdaxJKIHYGYEwcU1W++sJuWjq6eHbNEq4oykZEfm+7Ty8q4LIJGXzh5zu486n3eXjVTD6/tKjPbU38sDMAY+LYL7ae4M1DdfzlqlksmZJz0Q/0y8aP5jd/upSVs8bwjZf38+fP7RrGSo0bLACMiVPHG9r4xsv7uHpaDncvmexon4zUJJ783ELWfnQqL354kvcq6oe4SuMmCwBj4pDfr3z5+V14RXj0jnl4PM6bckSEhz5WTP7oVP7plYOoWu+geOUoAETkRhE5KCIVIvJwH+tnisgWEekQkS+HLJ8hIjtDfppE5KHguq+LyMmQdTdF7KiMSXBPv3uErUfO8De3zmZC5ogB75+a5OXBlcV8ePwcm/bXDkGFJhr0GwAi4gUeB1YBs4G7RGR22GZngAeBb4UuVNWDqjpfVecDC4E24KWQTb7bu15V1w/+MIwxvSpqm3l040E+NmsMdyycOOjXuWPhRCbnpPGtVw7itzECccnJGcBioEJVK1W1E3gWWB26garWquo2oOsir7MSOKyqfc5MY4y5dN09fr703C5GJnv55ifnXFIvniSvhy9dN50DNc38dk91BKs00cJJAEwAToQ8rwouG6g7gV+ELXtARHaLyNMiktXXTiKyRkRKRaS0rq5uEG9rTOL4lzcOs7uqkW98Yg5jRqVe8uvdOnc8M8eN4ruvHqK7xx+BCk00cRIAfX2FGND5oIgkAx8Hfhmy+AlgKjAfqAa+3de+qvqUqpaoakle3u/NaWyMCaqobeGxTeV8fN54bp6bH5HX9HiEL103nSP1rbywoyoir2mih5MAqAIKQp5PBE4N8H1WATtU9XTvAlU9rao9quoHfkSgqckYM0g/efcIHo/wN7eGX6K7NNfNHsu8gky+/1o5Hd09EX1t4y4nAbANKBaRouA3+TuBdQN8n7sIa/4RkdCvKLcBewf4msaYoMbzXby44yQfnzc+4vfyERG+cv0MTjW28x8fHI/oaxt39RsAqtoNPABsBPYDz6lqmYisFZG1ACIyTkSqgC8Bfy0iVSKSEVyXBlwHvBj20o+KyB4R2Q0sB/4sYkdlTIL5ZekJznf18AdXFQ7J6189LYcrp+Tw+OsVtHV2D8l7mOHn6F5AwS6a68OWPRnyuIZA01Bf+7YBOX0sv3tAlRpj+uT3K//2/jEWTs7iIxNGD8l7iAhfvmEGtz/xHj959yj3L582JO9jhpeNBDYmxr15qI5jDW3cO0Tf/nstnJzFyplj+OGbh2npsLOAeGABYEyM++l7RxkzKoVVHxk35O+1dtlUmtq7eW3f6f43NlHPAsCYGFZZ18Kbh+r47BWTSfIO/X/nhZOyyB+dym92DbQjoIlGFgDGxLCfbTlGkle464qC/jeOAI9HuHlOPm+V19HYdrGB/yYWWAAYE6NaOrp5fnsVN8/Jj8ioX6dumTeerh5l476aYXtPMzQsAIyJUS/uqKKlo3vIL/6GmzdxNAXZI/jtbrs/UKyzADAmBqkqz7x3lHkTR3P5pD5vozVkRIRb5o7n3Yp6Glo6hvW9TWRZABgTg96pqOdwXeuwf/vvdcvcfHr8yoYyawaKZRYAxsSgZ947Sm56csRu+jZQs/MzmJI7kt/usmagWGYBYEyMOd7QxqYDtdy1eBIpPq8rNQSagfJ5/0gDtU3trtRgLp0FgDEx5pfbTyDAZ69wNtH7ULl13nhUYb1NFhOzLACMiSGqysu7q7lqai7jRg9f18++FI8dxYyxo6w3UAyzADAmhuyvbqayvpWb5rjT9h/ulrn5lB47y6lz590uxQyCBYAxMeTlPafweoQbLhvrdilAYFAYWDNQrLIAMCZGqCrr99Rw5ZQcciI86ctgFeWO5CMTMuzeQDHKAsCYGLGvuokj9a2udf28kFvmjmdXVSPHG9rcLsUMkAWAMTHi5d3Vweafob/t80DcHLwe8ds9dhYQaywAjIkBgeafaq6amkP2yGS3y/kfCrLTmF+QaYPCYpCjABCRG0XkoIhUiMjDfayfKSJbRKRDRL4ctu5ocO7fnSJSGrI8W0ReFZHy4O/hvaGJMTGk7FQTRxva/uvbdrS5dd549lU3cbiuxe1SzAD0GwAi4gUeB1YBs4G7RGR22GZngAeBb13gZZar6nxVLQlZ9jCwSVWLgU3B58aYPqzfE2j+uT7Kmn969c5Gtnl/rcuVmIFwcgawGKhQ1UpV7QSeBVaHbqCqtaq6DRjIDBGrgWeCj58BPjGAfY1JGKrKy1Ha/NNrfOYIpuaN5O2KerdLMQPgJAAmACdCnlcFlzmlwCsisl1E1oQsH6uq1QDB32P62llE1ohIqYiU1tXVDeBtjYkPZaeaONbQxi1R1vsn3DXFeWw90kB7V4/bpRiHnASA9LFMB/AeV6vqAgJNSPeLyLUD2BdVfUpVS1S1JC8vbyC7GhMXXu5t/pkdnc0/vZZOy6W9y8+OY2fdLsU45CQAqoDQCUcnAo77e6nqqeDvWuAlAk1KAKdFJB8g+NsaD40J03vvn6un5ZIVpc0/vZZMzcHnEWsGiiFOAmAbUCwiRSKSDNwJrHPy4iIyUkRG9T4Grgf2BlevA+4NPr4X+PVACjcmEZSdauL4mTZunhPd3/4B0lN8XD4pk3fKLQBiha+/DVS1W0QeADYCXuBpVS0TkbXB9U+KyDigFMgA/CLyEIEeQ7nASyLS+17/oaobgi/9CPCciHweOA58KqJHZkwc+O3uanwx0PzTa+m0PL636RBnWzuj/ozFOAgAAFVdD6wPW/ZkyOMaAk1D4ZqAeRd4zQZgpeNKjUkwgd4/p2Ki+afX0uJcvvvaId49XM8tc8e7XY7ph40ENiZK7T3ZxIkz56N28Fdf5k0czahUnzUDxQgLAGOi1G/3nAo0/0TJrZ+d8Hk9XDklh7fL61EdSGdB4wYLAGOi1Gv7TnPl1Bwy02Kj+afXNcW5nDx3nqN2d9CoZwFgTBQ6Wt/K4bpWVs7sc3xkVFtaHBiv8065DdyMdhYAxkShzQcCw2JWzIyd5p9ehTlpTMgcwdt2HSDqWQAYE4U2H6ileEw6k3LS3C5lwESEa6fnsuVwA909frfLMRdhAWBMlGlu7+KDIw2smBV7zT+9lk7Lo7mjm11VjW6XYi7CAsCYKPNOeT1dPcrKGGz+6XXV1BxEsO6gUc4CwJgos+lALaNHJLFgUqbbpQxa1shk5kwYzTsVdiE4mlkAGBNF/H7l9QO1LJuRh88b2/89l07L5cPj52huH8g0IWY4xfZfmDFxZlfVORpaO1kRg90/wy0tzqXbr7xfecbtUswFWAAYE0U2H6jF6xE+Oj32575YODmLEUleGw8QxSwAjIkim/bXsnByVsyN/u1Lis/L4qJsmx8gilkAGBMlqhvPs6+6KSZH/17INcW5VNa1curcebdLMX2wADAmSvSO/l0Zw/3/wy0tzgWsO2i0sgAwJkps3l/LpOw0pualu11KxMwYO4rc9GTer2xwuxTTBwsAY6LA+c4e3qmoZ8XMMQRn0IsLIsLiomw+OGI9gaKRBYAxUWBLZT0d3f64av7ptagwm5PnznPSrgNEHUcBICI3ishBEakQkYf7WD9TRLaISIeIfDlkeYGIvC4i+0WkTES+GLLu6yJyUkR2Bn9uiswhGRN7Nu2vZWRyoNdMvOk9pm12FhB1+g0AEfECjwOrCEz0fpeIzA7b7AzwIPCtsOXdwJ+r6ixgCXB/2L7fVdX5wZ/1GJOAVJXNB2q5pjiPFJ/X7XIibua4DEal+qwZKAo5OQNYDFSoaqWqdgLPAqtDN1DVWlXdBnSFLa9W1R3Bx83AfmBCRCo3Jk7sr26murE9pu/+eTFej1AyOYutR+xCcLRxEgATgBMhz6sYxIe4iBQClwMfhCx+QER2i8jTIpJ1gf3WiEipiJTW1dmIQhN/Nh84DcDyGfEZAACLi3I4XNdKfUuH26WYEE4CoK8uCQOa7VlE0oEXgIdUtSm4+AlgKjAfqAa+3de+qvqUqpaoakleXuwPjzcm3KYDtcwryCRvVIrbpQyZ3usApUetGSiaOAmAKqAg5PlE4JTTNxCRJAIf/j9X1Rd7l6vqaVXtUVU/8CMCTU3GJJQzrZ3sPHGO5TPi+8vNnAmjSU3ysPXIWbdLMSGcBMA2oFhEikQkGbgTWOfkxSXQoflfgf2q+p2wdfkhT28D9jor2Zj48dahOlTju/kHINnn4fKCLLYetesA0cTX3waq2i0iDwAbAS/wtKqWicja4PonRWQcUApkAH4ReYhAj6G5wN3AHhHZGXzJrwV7/DwqIvMJNCcdBf44gsdlTEx442AtOcHJU+Ld4qJsfrC5nOb2LkalJrldjsFBAAAEP7DXhy17MuRxDYGmoXDv0Pc1BFT1budlGhN/evzKm4fqWD5jDB5P/Iz+vZDFRdn4FbYfO8uyOD/jiRU2EtgYl+yuOsfZti4+Guft/70un5SJzyNstfEAUcMCwBiXvHGwDo/AtcWJEQBpyT7mTBxtARBFLACMcckbB2uZX5BJ1sjYn/zFqcWF2eyuaqS9q8ftUgwWAMa4oqGlg90nGxOuLXxxUTadPX52njjndikGCwBjXPFWeaD757IEaf/vVTI5GxG7MVy0sAAwxgWvH6gjNz2Zj4yP/+6foUanJTFj7Ci22ojgqGABYMww6/Erb5XXce30vITo/hnuiqJsth87S3eP3+1SEp4FgDHDbOeJc5xr60q49v9ei4qyaevsoexUU/8bmyFlAWDMMHvzYG2w+2eu26W4YnFh4MZw1h3UfRYAxgyzNw7VcfmkLDLTEqf7Z6gxGakU5Y60CWKigAWAMcOorrmD3VWNLJueWL1/wi0qzKL02Bn8/gHdWd5EmAWAMcPorUOBSY2Wz0zM9v9ei4tyONfWRXlti9ulJDQLAGOG0RuH6shNT2F2fobbpbjqiuAEMdYd1F0WAMYMk+4eP28dquOjCdr9M9TErBGMy0jlg0qbH8BNFgDGDJNdVedoPN/F8pmJ3f4PICIsKsqm9OhZVO06gFssAIwZJq8fCNz985ppFgAAiwuzqGlqp+rsebdLSVgWAMYMkzcO1bJgUhaj02w2LICS4HiAbXYdwDUWAMYMg9rmdvaebEq4m79dzIyxoxiV6mPbUZso3i2OAkBEbhSRgyJSISIP97F+pohsEZEOEfmyk31FJFtEXhWR8uDvrEs/HGOi0xsHAt0/V8wc63Il0cPjEUomZ9kZgIv6DQAR8QKPA6sITPR+l4jMDtvsDPAg8K0B7PswsElVi4FNwefGxKXNB2rJH53KrPxRbpcSVRYVZVNR28KZ1k63S0lITs4AFgMVqlqpqp3As8Dq0A1UtVZVtwFdA9h3NfBM8PEzwCcGdwjGRLeO7h7eLq9j+cwxiCR2989wi4LXAUrtLMAVTgJgAnAi5HlVcJkTF9t3rKpWAwR/9zk0UkTWiEipiJTW1dU5fFtjose2I2dp7exhRYLe/fNi5k4cTbLPQ+kxuw7gBicB0NdXFqcddy9l38DGqk+paomqluTl2QU0E3s2HThNss/DVdNy3C4l6qT4vMyzieJd4yQAqoCCkOcTgVMOX/9i+54WkXyA4O9ah69pTMxQVTYfqOWqqTmkJfvcLicqlRRms/dkI+c7baL44eYkALYBxSJSJCLJwJ3AOoevf7F91wH3Bh/fC/zaednGxIbK+laONbSxIsFv/nYxiwuz6fYrH56wZqDh1u9XElXtFpEHgI2AF3haVctEZG1w/ZMiMg4oBTIAv4g8BMxW1aa+9g2+9CPAcyLyeeA48KkIH5sxrnv9QODEdrm1/1/QgslZiEDp0bNcNTUxJ8lxi6NzUlVdD6wPW/ZkyOMaAs07jvYNLm8AVg6kWGNizab9tUwfm05BdprbpUSt0SMCE8XbeIDhZyOBjRkiTe1dbDt6xgZ/ObCoMJsdNlH8sLMAMGaIvFNeT7dfrf3fgUVF2bR29nCgptntUhKKBYAxQ2TT/lpGj0hiwaRMt0uJeosKA3eCse6gw8sCwJgh4Pcrbxys5aPT8/B57b9Zf/JHj2Bi1ghKj1kADCf7yzRmCOyqOkdDaycrZ1nzj1OLCrPZesQmiBlOFgDGDIHXD9TiEfjodBu97lRJYRb1LR0ca2hzu5SEYQFgzBDYdKCWhZOzyExLdruUmLG40CaKH24WAMZEWE1jO2WnmlhuvX8GZGpeOplpSXZn0GFkAWBMhL1+MDD6d6X1/x+QwAQx2TZD2DCyADAmwjYfqGVC5gimj013u5SYs6gwiyP1rdQ1d7hdSkKwADAmgtq7eninvJ4VNvnLoCwqsglihpMFgDER9H5lA+e7emz07yB9ZPxoUpM81gw0TCwAjImgjWWnGZns5cqpNvnLYCT7PMwvyLQbww0TCwBjIqTHr7y6r4blM8eQmuR1u5yYtbgwm7JTjTS3h08xbiLNAsCYCNl+7Cz1LZ3c+JFxbpcS05ZMycGvgfkBzNCyADAmQjbsrSHZ52GZTf5ySRZMziLZ62FLZYPbpcQ9CwBjIkBV2VhWw7XFuaSn2Ny/lyI1ycv8SZlsOWwBMNQcBYCI3CgiB0WkQkQe7mO9iMhjwfW7RWRBcPkMEdkZ8tMUnC4SEfm6iJwMWXdTRI/MmGG092QTJ8+d54bLrPknEq6ckkPZqUYaz9t1gKHUbwCIiBd4HFgFzAbuEpHZYZutAoqDP2uAJwBU9aCqzlfV+cBCoA14KWS/7/auD04daUxM2lBWjdcjfGyWjf6NhN7rANtsfoAh5eQMYDFQoaqVqtoJPAusDttmNfAzDXgfyBSR/LBtVgKHVfXYJVdtTJTZsLeGJVOyyRppN3+LhMsnZZLss+sAQ81JAEwAToQ8rwouG+g2dwK/CFv2QLDJ6GkRyerrzUVkjYiUikhpXV2dg3KNGV4Vtc0crmvlRmv+iZjUJC8LJmXyvgXAkHISAH2NZw+fseGi24hIMvBx4Jch658ApgLzgWrg2329uao+paolqlqSl2f3VjfRZ8PeGgCutwCIqCun5LKvuolzbZ1ulxK3nARAFVAQ8nwicGqA26wCdqjq6d4FqnpaVXtU1Q/8iEBTkzExZ0NZDQsmZTI2I9XtUuLKlVNzUIUP7DrAkHESANuAYhEpCn6TvxNYF7bNOuCeYG+gJUCjqlaHrL+LsOafsGsEtwF7B1y9MS47caaNvSebbPDXEJhXMJoUn8eagYZQvx2WVbVbRB4ANgJe4GlVLRORtcH1TwLrgZuACgI9fe7r3V9E0oDrgD8Oe+lHRWQ+gaaio32sNybqbSwLNP9Y98/IS/F5KSnMsvEAQ8jRiJVgF831YcueDHmswP0X2LcN+L07Y6nq3QOq1JgotLGshln5GUzOGel2KXHpyik5fOuVQ5xt7bQeVkPARgIbM0i1ze2UHjtrvX+G0JIpge+OHxyxs4ChYAFgzCC9uu80qlj7/xCaOzGTEUleawYaIhYAxgzShr01FOWOtKkfh1Cyz0NJYRbvV1pPoKFgAWDMIDS2dbHlcAM3XDbOpn4cYkum5HDwdDMNLTZPcKRZABgzCK/tP023X635Zxj0zq5mZwGRZwFgzCCs23WKCZkjmDthtNulxL05E0aTluy18QBDwALAmAGqbWrn7fI6brt8Ah6PNf8MtSSvh0WF2XZjuCFgAWDMAK3bdQq/wm0Lwu93aIbKkik5VNS2UNds1wEiyQLAmAF6YcdJ5hVkMjXPev8Ml/++DmBnAZFkAWDMAOyvbmJ/dROfvNy+/Q+nj4zPID3FZ81AEWYBYMwAvPThSXwe4dZ5490uJaH4vB4WFWbZGUCEWQAY41CPX/nVhydZPnMM2XZfmmF35dQcKutaOd3U7nYpccMCwBiH3q2op7a5w5p/XLJ0WmBCqDcP2syAkWIBYIxDL+6oIiPVx4pZY9wuJSHNyh9F/uhUNh043f/GxhELAGMcaOnoZkNZDbfMG0+Kz+t2OQlJRFg+cwzvlNfT0d3jdjlxwQLAGAc27K2hvcvP7db331UrZoyhtbOHbUfOul1KXLAAMMaBF3dUMTknjQWTstwuJaFdNS2HZJ+HzQdq3S4lLlgAGNOPU+fOs6Wygdsun2B3/nRZWrKPK6fksNmuA0SEowAQkRtF5KCIVIjIw32sFxF5LLh+t4gsCFl3VET2iMhOESkNWZ4tIq+KSHnwt321MlHpVztPogq3We+fqLBy1hiONrRRWdfidikxr98AEBEv8DiwCpgN3CUis8M2WwUUB3/WAE+ErV+uqvNVtSRk2cPAJlUtBjYFnxsTVVSVF3ecpGRyls37GyWWzwj0wrJmoEvn5AxgMVChqpWq2gk8C6wO22Y18DMNeB/IFJH8fl53NfBM8PEzwCecl23M8Nh7somK2hY+uWCi26WYoILsNIrHpPP6QQuAS+UkACYAJ0KeVwWXOd1GgVdEZLuIrAnZZqyqVgMEf/fZuVpE1ohIqYiU1tXZABAzvF7YUUWy18PNc/r7PmOG04qZY9h65AzN7V1ulxLTnARAX1e9dADbXK2qCwg0E90vItcOoD5U9SlVLVHVkry8vIHsaswlae3o5oUdVVx/2VhGpyW5XY4JsWLmGLp6lHfK690uJaY5CYAqoCDk+UTglNNtVLX3dy3wEoEmJYDTvc1Ewd92Pmeiygs7qmhu7+a+q4vcLsWEWTg5i4xUn10HuEROAmAbUCwiRSKSDNwJrAvbZh1wT7A30BKgUVWrRWSkiIwCEJGRwPXA3pB97g0+vhf49SUeizER4/crP3n3KPMKMlkwKdPtckwYn9fDtdPzeP1gHX5/eIOEcarfAFDVbuABYCOwH3hOVctEZK2IrA1uth6oBCqAHwFfCC4fC7wjIruArcDLqrohuO4R4DoRKQeuCz43Jiq8eaiOI/Wt/OHVhdb3P0qtmDmG+pYO9p5qdLuUmOVzspGqrifwIR+67MmQxwrc38d+lcC8C7xmA7ByIMUaM1yefvcIYzNSuMku/katZTPGIAKb9tcyd2Km2+XEJBsJbEyYQ6ebebu8nnuuLCTJa/9FolX2yGQuL8i07qCXwP66jQnzk3ePkuLzcNfiSW6XYvqxYuYYdlc1Uttsk8QMhgWAMSHOtnby4o4qPrlggs36FQOWzwwMH3rDJokZFAsAY0L8x9bjdHT7retnjJidn8G4jFRet+6gg2IBYExQV4+ff9tyjGuKc5k+dpTb5RgHeieJebu8ns5uv9vlxBwLAGOCfre3hpqmdu67utDtUswArJg5hpaObrYeOeN2KTHHAsCYoKffOUJR7kiWTbc5f2PJ0mm5jEz2sm7XSbdLiTkWAMYAO46fZeeJc9x3dSEejw38iiUjkr3cPDefl3dX09bZ7XY5McUCwBgCXT9Hpfq43W77HJPuWFhAa2cPG/bWuF1KTLEAMAmv/HQzL+8+xWcWT2JkiqPB8SbKLCrMYlJ2Gs9vr3K7lJhiAWAS3j9uOMjIZB9//NGpbpdiBklEuGPhRLZUNlB1ts3tcmKGBYBJaFuPnOG1/adZu2yqDfyKcZ9cMAFVeGmHXQx2ygLAJCxV5R9+t59xGan8oQ38inkTs9K4ckoOz++oInB/StMfCwCTsDbsreHD4+f40nXTGZHsdbscEwF3LJzIsYY2So+ddbuUmGABYBJSV4+fRzceZPrYdG5faD1/4sWqOeMYmezl+VK7GOyEBYBJSM9uPc6R+lYeXjUTr/X7jxtpyT5umpPPy3tsTIATFgAm4bR0dPO918q5oiib5TNs1G+8uWPhRFo6utlYZmMC+mMBYBLOU29V0tDayV/eNMume4xDiwqzbUyAQ44CQERuFJGDIlIhIg/3sV5E5LHg+t0isiC4vEBEXheR/SJSJiJfDNnn6yJyUkR2Bn9uitxhGdO32qZ2fvx2JTfPzWd+Qabb5Zgh4PEIty+YyHuHGzh57rzb5US1fgNARLzA48AqYDZwl4jMDttsFVAc/FkDPBFc3g38uarOApYA94ft+11VnR/8+R9zDhszFL63qZzObj9fuX6G26WYIfTfYwLsLOBinJwBLAYqVLVSVTuBZ4HVYdusBn6mAe8DmSKSr6rVqroDQFWbgf3AhAjWb4xje0828p/bTvDZKyZRmDvS7XLMECrITmPJlGye325jAi7GSQBMAE6EPK/i9z/E+91GRAqBy4EPQhY/EGwyelpEsvp6cxFZIyKlIlJaV2fTvpnBOd/Zw4PPfkhuejIPfWy62+WYYXDHwgKONrSx3cYEXJCTAOjrKll4pF50GxFJB14AHlLVpuDiJ4CpwHygGvh2X2+uqk+paomqluTl5Tko15jf9/fr91FZ18p3Pj2fLLvlQ0JY9ZHAmICfbTnmdilRy0kAVAEFIc8nAqecbiMiSQQ+/H+uqi/2bqCqp1W1R1X9wI8INDUZE3Gv7TvNv79/nP99TRFXT8t1uxwzTEam+LjnqkJ+s/sU+6ub+t8hATkJgG1AsYgUiUgycCewLmybdcA9wd5AS4BGVa2WQB+7fwX2q+p3QncQkfyQp7cBewd9FMZcQG1zO199YTez8jP48g124TfRrL12KqNSfHxr40G3S4lK/QaAqnYDDwAbCVzEfU5Vy0RkrYisDW62HqgEKgh8m/9CcPnVwN3Aij66ez4qIntEZDewHPiziB2VMQRu9vaVX+6mpaObx+6cT4rP7veTaEanJbF22VQ2Hahl21GbMzico9kvgl0014ctezLksQL397HfO/R9fQBVvXtAlRozQD/bcow3D9Xxf1dfRvHYUW6XY1xy31VF/PTdozy64QDP/fGVNvgvhI0ENnHp0Olmvrl+P8tn5HH3kslul2NcNCLZy5+uLGbb0bO8cdB6EoayADBxp6O7hwd/8SHpKT4evWOefeMz3LmogMk5afzjhgP4/TYuoJcFgIkr3T1+/uw/d3KgpplH75hL3qgUt0syUSDJ6+FL103nQE0zv9kd3okxcVkAmLjR41e+8vxu1u+p4a9vnsXKWWPdLslEkVvnjmdWfgbffuUQnd1+t8uJChYAJi74/cpfvbSHlz48yVdumMEfXTPF7ZJMlPF4hL+4YQbHz7Txn6Un+t8hAVgAmJinqvzdb8p4dtsJ/nTFNO5fPs3tkkyUWjYjj0WFWTy2qdwmjMECwMQ4VeWR3x3gmS3H+KOlRXzpOrvPj7kwEeEvbpxJXXMHP3n3qNvluM4CwMS0775Wzg/fquTuJZP5q5ttghfTv0WF2Vw3eyzff62c7ccSe3CYBYCJSV09fv5h/X4e21TOpxZO5O8+fpl9+BvHHr19LhOyRrDmZ9s5cabN7XJcYwFgYs6JM218+odb+OFblXzmikk8cvtcPDaxuxmArJHJ/Ou9JXT7lT/86Taa2rvcLskVFgAmpvx29ylu+v7bVNS28M+fuZxv3jYHr334m0GYkpfOE59bwJH6Vu7/+Q66ehKva6gFgIkJbZ3dPPzCbh74jw+ZNjad9Q9ewy1zx7tdlolxV03N5Zu3zeHt8nq+vq4s4WYPc3QzOGPctPdkI1989kMq61u5f/lUHvrYdJK89t3FRManFxVwuL6FH75ZyZS8dD6/tMjtkoaNBYCJWntPNvKDzeVsLDtN3qgU/v3zV9iELmZIfPWGmRytb+UbL++jMCctYUaRWwCYqPPh8bP8YHMFmw/UMirVx4Mri/n81UWMTktyuzQTpzwe4bv/az6f/uEW/uTnO/jiymLWXDsl7s80LQBMVOjxKx9UNvDEm4d5u7yezLQkvnz9dO65qpCMVPvgN0MvLdnHT+9bzP/51V7+aeNB1u08xTc/OYeFk7PcLm3IWAAY17R39fBuRT2vlJ3mtf2naWjtJDc9mb9cNZPPLZnMyBT78zTDKzc9hSc+t5BX953mb369lzuefI/PXTGZr9w4Iy6/iNj/MDNsevzKkfoWdp5o5LV9p3nzUB3nu3oYleJj+cwxXDd7LB+bNZYRyTZ1o3HXdbPHcuXUHL79ykF++t5RXtlXw9/eehnXzx6LL46ahRwFgIjcCHwf8AI/VtVHwtZLcP1NQBvwB6q642L7ikg28J9AIXAU+LSqnr30QzJu8/uVM22dnDjTxr7qJspONbHvVBMHappo7wr0tR6bkcLtCydw/exxLJmSQ7Ivfv5TmfiQnuLjb2+9jE/Mn8DDL+7hCz/fQUaqj6un5XJNcR7XFOdSkJ3mdpmXpN8AEBEv8DhwHVAFbBORdaq6L2SzVUBx8OcK4Angin72fRjYpKqPiMjDwedfjdyhmcFSVbr9Ske3n87gT0d3D60dPTS3d9HS0U1zezfN7V00tXdztrWTmqZ2ahrbqWlqp7apg86QQTUZqT5mj8/gM4snM3t8BpeNz2DG2FE2etfEhHkFmax74GpeKTvNW4fqeKu8jt/trQGgKHckS6flMiVvJDnpKeSmJ5ObnkJuegqZI5Ki/m/cyRnAYqBCVSsBRORZYDUQGgCrgZ8FJ4d/X0QyRSSfwLf7C+27GlgW3P8Z4A2GKAB+sKmcdbtiZxagCw1FudggFQ170PtcVfErKIoqwZ/Ash5V/P7Ah73fr/So0uNXOnv8DGQ8TGqSh3EZqYwbnUrJ5CzGjR7BuIwUxmeOYFZ+BhOzRth9ekxMS/J6uHluPjfPzUdVOVzXytvldbxdXs/z26s439Xze/t4PUJakhevV/B5PPg8gtcjJHkFjwgE/0v0/s/o/T9yof8p3/zkHBYVZkf0uJwEwAQgdPaEKgLf8vvbZkI/+45V1WoAVa0WkTF9vbmIrAHWAEyaNMlBub8vb1QKxWPTB7WvW+RCfwYX+Ry90B+SRwLLAn9zEnwe+AP1egSvCJ7gb69HSPZ5SPZ6SEkK/E72eUn2eRiZ7GVUahKjUn2kp/oYleojIzWJFJ/HPuBNwhARpo1JZ9qYdO67ugi/X2k830V9Swf1LZ3Ut3TQEHzc1tlDj99Pl1/p6VG6/H56/IEvWtDXF7cLf/MakRT5a2NOAqCv/9nhVV5oGyf7XpSqPgU8BVBSUjKocdp3Lp7EnYsHFx7GGHMxHo+QNTKZrJHJFMfY+DEnV96qgIKQ5xOB8PaUC21zsX1PB5uJCP6udV62McaYS+UkALYBxSJSJCLJwJ3AurBt1gH3SMASoDHYvHOxfdcB9wYf3wv8+hKPxRhjzAD02wSkqt0i8gCwkUBXzqdVtUxE1gbXPwmsJ9AFtIJAN9D7LrZv8KUfAZ4Tkc8Dx4FPRfTIjDHGXJSjcQCqup7Ah3zosidDHitwv9N9g8sbgJUDKdYYY0zk2OgbY4xJUBYAxhiToCwAjDEmQVkAGGNMgpJYmgNTROqAY4PcPReoj2A5scCOOTHYMSeGSznmyaqaF74wpgLgUohIqaqWuF3HcLJjTgx2zIlhKI7ZmoCMMSZBWQAYY0yCSqQAeMrtAlxgx5wY7JgTQ8SPOWGuARhjjPmfEukMwBhjTAgLAGOMSVBxFwAicqOIHBSRiuBcw+HrRUQeC67fLSIL3Kgzkhwc82eDx7pbRN4TkXlu1BlJ/R1zyHaLRKRHRO4YzvoizcnxisgyEdkpImUi8uZw1xhpDv6uR4vIb0RkV/CY73OjzkgSkadFpFZE9l5gfWQ/v1Q1bn4I3HL6MDAFSAZ2AbPDtrkJ+B2B2cqWAB+4XfcwHPNVQFbw8apEOOaQ7TYTuBvtHW7XPcT/xpkE5tqeFHw+xu26h+GYvwb8Y/BxHnAGSHa79ks87muBBcDeC6yP6OdXvJ0B/NcE9qraCfROQh/qvyawV9X3gd4J7GNVv8esqu+p6tng0/cJzMwWy5z8OwP8KfACsT/bnJPj/QzwoqoeB1DVRDhmBUZJYELqdAIB0D28ZUaWqr5F4DguJKKfX/EWABeanH6g28SSgR7P5wl8g4hl/R6ziEwAbgOeJPY5+TeeDmSJyBsisl1E7hm26oaGk2P+Z2AWgWlm9wBfVFX/8JTnmoh+fjmaECaGXMoE9rHK8fGIyHICAbB0SCsaek6O+XvAV1W1J/AFMaY5OV4fsJDAJEsjgC0i8r6qHhrq4oaIk2O+AdgJrACmAq+KyNuq2jTEtbkpop9f8RYAlzKBfaxydDwiMhf4MbBKA7OxxTInx1wCPBv88M8FbhKRblX91bBUGFlO/67rVbUVaBWRt4B5QKwGgJNjvg94RAON4xUicgSYCWwdnhJdEdHPr3hrArqUCexjVb/HLCKTgBeBu2P4G2Gofo9ZVYtUtVBVC4HngS/E6Ic/OPu7/jVwjYj4RCQNuALYP8x1RpKTYz5OcFpZERkLzAAqh7XK4RfRz6+4OgPQS5jAPlY5POa/AXKAfwl+I+7WGL6TosNjjhtOjldV94vIBmA34Ad+rKp9diWMBQ7/jf8f8FMR2UOgaeSrqhrTt4gWkV8Ay4BcEakC/hZIgqH5/LJbQRhjTIKKtyYgY4wxDlkAGGNMgrIAMMaYBGUBYIwxCcoCwBhjEpQFgDERErwnDSLy9dDnxkQrCwBjIme+iDwGZIvIJ4C/d7keYy4qrgaCGTNcRKQQ2AB8AFxO4JYL9wD/AmwBklT1T1wr0BgHbCCYMYMQDIAjwFJVfVdEniZ4P34Co1I3ASWq+tfuVWnMxVkTkDGDd0JV3w0+/nfgGlV9EGgI3nfo/7hWmTEOWAAYM3jhp89+AFX9evC3nV6bqGYBYMzgTRKRK4OP7wLecbMYYwbKAsCYwdsP3Csiu4Fs4AmX6zFmQKwXkDGD51fVtW4XYcxg2RmAMcYkKOsGaowxCcrOAIwxJkFZABhjTIKyADDGmARlAWCMMQnKAsAYYxLU/wedDVeFnGGMTgAAAABJRU5ErkJggg==\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "x = np.linspace(0,1,50)\n", | |
| "likelihoods = np.array([P(20,13,p) for p in x])\n", | |
| "\n", | |
| "plt.plot(x, likelihoods)\n", | |
| "plt.axvline(13/20, color = 'green')\n", | |
| "plt.xlabel('p*');" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "id": "34f8e4bf", | |
| "metadata": {}, | |
| "source": [ | |
| "From this chart, we can see that the most likely value of $p$ is indeed 0.65, however reasonable estimates of $p$ could fall anywhere between approximately 0.45 and 0.8.\n", | |
| "\n", | |
| "This distibution is the Beta distribution which takes two parameters:\n", | |
| "* α = 1 + successes = (1 + heads)\n", | |
| "* β = 1 + failures = (1 + tails)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 5, | |
| "id": "71a83498", | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "β = beta(1+13,1+7)\n", | |
| "β_distribution = β.pdf(x)/β.pdf(x).sum()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "id": "242d18e4", | |
| "metadata": {}, | |
| "source": [ | |
| "We can demonstrate the equivalnce of these two distributions with an overlay plot of their normalzied values." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 6, | |
| "id": "f79bf29e", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA+TElEQVR4nO3dd3hUZfbA8e9Jo7eEgCSABESwIAgJ1cKAICCCupBF10XdVda2uj9dwcbqrq4lllVXRUFRsYABdUVBsDCIICVBAUFAQoISAkLoLZByfn/MBIeQMoFMJjNzPs8zTzL3vnfuuSlz5r5VVBVjjDGmpDB/B2CMMaZmsgRhjDGmVJYgjDHGlMoShDHGmFJZgjDGGFOqCH8HUJWaNm2qbdq08XcYxhgTMJYvX56rqrGl7QuqBNGmTRvS09P9HYYxxgQMEfm5rH1WxWSMMaZUliCMMcaUyhKEMcaYUvm0DUJEBgHPA+HAa6r6RIn94t4/BDgEXK+q37n3/R9wI6DAD8ANqprny3iNMTVffn4+2dnZ5OXZ20Fl1K5dm5YtWxIZGen1MT5LECISDrwEDACygTQRmamqP3oUGwy0dz96ABOAHiISD9wBnK2qh0UkFRgFvOmreI0xgSE7O5sGDRrQpk0bXJ8xTUVUlZ07d5KdnU1CQoLXx/myiqk7kKGqmap6FJgGDC9RZjgwRV2WAI1FpIV7XwRQR0QigLpAjg9jNcZ4SFmUgjPLedw2Z5aTlEUpforoN3l5ecTExFhyqAQRISYmptJ3Xb5MEPHAZo/n2e5tFZZR1S3A08AvwFZgr6p+XtpJRGSMiKSLSPqOHTuqLHhjQllSXBLJM5KPJQlnlpPkGckkxSX5OTIXSw6VdzI/M18miNKiKTm3eKllRKQJrruLBCAOqCci15Z2ElWdqKqJqpoYG1vqWA9jTCU5Ehw8N+Athk8dwd/n3k/yjGRSR6TiSHD4OzRTjXyZILKBVh7PW3JiNVFZZS4BslR1h6rmAx8CvX0YqzEhp7RqpLkZX3L9jAcZNXExD0xTwg4N5JkljzP8jOstOXgIDw+nS5cudO7cma5du/Ltt98CkJOTw4gRI3x67vT0dO64445yy8yfP5+hQ4ee8rl8mSDSgPYikiAiUbgamWeWKDMTGC0uPXFVJW3FVbXUU0Tquns69QfW+jBWY0KOZzXSxh0H+PPUN7js3RHM+a4eW/YcZniP3Uj9Lzg98o+8sXIid3/8DrbAmEudOnVYsWIFK1eu5PHHH+e+++4DIC4ujhkzZvj03ImJibzwwgs+PUcxnyUIVS0Abgfm4npzT1XVNSJys4jc7C42G8gEMoBJwK3uY5cCM4DvcHVxDQMm+ipWY0KRI8FB6ohUhk8bQbfnb+LN9X9jSNyTpF5/Iw+PgPcy/o8Pfz+dVXe9zmXxKTz3/W2c/ezv+GTdF8e9Tk1pvPaXffv20aRJEwA2bdrEueeeC8Cbb77JVVddxaBBg2jfvj1jx449dszUqVPp1KkT5557LuPGjTu2vX79+owbN45u3bpxySWXsGzZMvr27Uvbtm2ZOdP1+drz7mDZsmX07t2b888/n969e7N+/foqvTafjoNQ1dm4koDntlc8vlfgtjKOfQh4yJfxGRPq2jbsTlTepeyMnMpdPe7jmUF/ASBlUfpxbQ4f3ziGsZ804NX0t7nq/ZG8OuRt/pR0+bHG69QRqX6J/5+frOHHnH1V+ppnxzXkocvPKbfM4cOH6dKlC3l5eWzdupV58+aVWm7FihV8//331KpViw4dOvDXv/6V8PBwxo0bx/Lly2nSpAkDBw7kf//7H1dccQUHDx6kb9++PPnkk1x55ZU8+OCDfPHFF/z4449cd911DBs27LjX79ixIwsWLCAiIoIvv/yS+++/nw8++KDKfhZBNVmfMcZ7hUXKH9+ZzB6Zzd963MuUHyYxtMMAHAkOxvYZe1xZEeGpYdcwqvMgrn37DcbMupb5m27ks01TQrLxuriKCWDx4sWMHj2a1atXn1Cuf//+NGrUCICzzz6bn3/+mZ07d9K3b1+KO9X84Q9/YMGCBVxxxRVERUUxaNAgADp16kStWrWIjIykU6dObNq06YTX37t3L9dddx0bNmxARMjPz6/S67QEYUyIGvvJuyza/SD/6D2JhwYmM6zDwAp7K3U7PZqFf7udpP+u4u0fn+X+Cx70a3Ko6JN+dejVqxe5ubmU1s2+Vq1ax74PDw+noKCg3HacyMjIY91Rw8LCjh0fFhZGQUHBCeXHjx+Pw+Hgo48+YtOmTfTt2/cUr+Z4NheTMSHop1/381b6PAa3eJJ/DBgJ/NYmkZaTVu6xq3Z8yy75lEb5o/jv0pdP6AkVatatW0dhYSExMTFele/Rowdff/01ubm5FBYWMnXqVC6++OKTOvfevXuJj3cNL3vzzTdP6jXKY3cQxoSY/MIi7kpdQeuoq3nzDxcdN4DKkeAo946guM3hw+TpvDS3Nj/u7Eby9GRSR4ZWNVNxGwS4prF46623CA8P9+rYFi1a8Pjjj+NwOFBVhgwZwvDhJSeZ8M7YsWO57rrrePbZZ+nXr99JvUZ5JJi6rSUmJqotGGRM+f7zxU88/9UGXrm2K4PObVHxAR5SFqWQFJeEI8HB4o07uXrSEkb03kOzmC0ntFv4ytq1aznrrLOq5VzBprSfnYgsV9XE0spbFZMxQc5zQNwP2Xt50ZlB4pk5rNr7dqVfa2yfscfuFHq1i+HC9k2Zt7IZtybeVaUxm5rBEoQxQa54QNzcDV9yV+oKouquZX7u/VUyr9LfB3Zg18GjTF6YVQWRmprGEoQxQa648fl3qcmk7XqFrRGPM72K2gw6t2rMwLObM2lBJrsPHq2CaE1NYgnCmBDQM/4i6hUMZm/kNO7ocWuVNijfPbADB44W8MqCjVX2mqZmsARhTAh47Kvp5Oqn3NDpbiakT6jSrqkdTmvA8M5xvPXtJrbvs1XegoklCGOC3LzMeTyZdjNd6/+T1698itQRqcet9VAVakfPYn/RCl50ZhzbFupzNAUDSxDGBLkP13xNdN447rroKkTE6wFxlXFp+z7srvMUr6d9wuZdh2rcAkNVrazpvsuyZ88eXn755WqKrupYgjAmyOm+4TSv1Y0ruvy2oGNp8y2dCkeCgylXvMevEU/w+2n/V2MWGPLV0qllTfddFksQxpgaZ+vew8xZs43fJ7WiTpR3I31P1ohzB9G92SgW75jITef/xe/JAapn6VTP6b4BnnrqKZKSkjjvvPN46CHXhNT33nsvGzdupEuXLtxzzz0cOHCA/v3707VrVzp16sTHH39cZfFUJZtqw5gg9u6SXyhS5Y89T/f5uZxZTn7c9wGN8kfxUtoEBrTr7/ckUVydljwjmVsSb2FC+oQqubMpa7rvzz//nA0bNrBs2TJUlWHDhrFgwQKeeOIJVq9efWwG2IKCAj766CMaNmxIbm4uPXv2ZNiwYTVurW27gzAmSOXlFzJ12S/079icVtF1fXqu4k/mM0a+T7taf8LR9LEqbwg/WY4EB7ck3sIjCx7hlsRbqiRpFVcxrVu3jjlz5jB69GhUlc8//5zPP/+c888/n65du7Ju3To2bNhwwvGqyv333895553HJZdcwpYtW/j1119POa6q5tMEISKDRGS9iGSIyL2l7BcRecG9f5WIdHVv7yAiKzwe+0Tkb76M1ZhgM2vVVnYePMr1vdv4/FxpOWmkjkjlknb9GdzpNNb90pq3hk+t0obwk+XMcjIhfQLjLxpf5V184fjpvlWV++67jxUrVrBixQoyMjL485//fMIx7777Ljt27GD58uWsWLGC5s2bk5dX87oI+yxBiEg48BIwGDgbuFpEzi5RbDDQ3v0YA0wAUNX1qtpFVbsA3YBDwEe+itWYYKOqvLV4E2c0q0+fM7ybhvpUeM7RdPl5ceTlF1Fw+Oxqm8CvLJ4r3v3L8S+fdPH1nO770ksvZfLkyRw4cACALVu2sH37dho0aMD+/fuPHbN3716aNWtGZGQkTqeTn3/+ucriqUq+bIPoDmSoaiaAiEwDhgM/epQZDkxxLz26REQai0gLVd3qUaY/sFFVa+ZP0Jga6PvNe1iVvZdHhp9T7fXaSW2iad6wFp+szGFY57hqPXdJxXc2xcnLs4vvqVQ1lTXd98CBA1m7di29evUCXGtMv/POO7Rr144+ffpw7rnnMnjwYMaNG8fll19OYmIiXbp0oWPHjqd8rb7gywQRD2z2eJ4N9PCiTDzgmSBGAVPLOomIjMF190Hr1q1PIVxjApvnVNxvfbuJBrUiiInJIGVRarV+kg8LE4Z0asG7S35hX14+DWtHVtu5Syrtuita88IbhYWFZe678847ufPOO0/Y/t577x33fPHixacUQ3XwZRtEaR9bSi4+UW4ZEYkChgHTyzqJqk5U1URVTSxe49WYUFTcpfPDNXOZ/cNWunXI4bqPr/bLYLWh58VxtLCIL9bUvIZX4z1fJohsoJXH85ZATiXLDAa+U1X7KzOmAsXVJ6P/dzU75G0+2XyP3wardW3dmPjGdfh0Vcl/eRNIfJkg0oD2IpLgvhMYBcwsUWYmMNrdm6knsLdE+8PVlFO9ZIw53sWn96Vx0WXsjZzGbd2rdtbWyhARhp7Xgm825PpkGvBgWgmzupzMz8xnCUJVC4DbgbnAWiBVVdeIyM0icrO72GwgE8gAJgG3Fh8vInWBAcCHvorRmGAzadknbC2cycgz7/RJl87KGHpeHAVFytw126r0dWvXrs3OnTstSVSCqrJz505q165dqeNsTWpjgoQzy8nQ935Hk7yxrH3gbtK3LvTrnEiqSt+n59OqSV3eubFk/5STl5+fT3Z2do0cN1CT1a5dm5YtWxIZeXyngfLWpLapNowJEsu2LKO1PsAF7S+mQe3IKuvSebKKq5kmzN9I7oEjNK1fq0peNzIykoSEhCp5LVM+m2rDmCBxSaubOHygI5d1anFsW1XP2lpZl3eOo0jhsx+2VlzY1DiWIIwJEp+u2kpURBj9z2rm71CO6dC8AWc0q88nqyxBBCJLEMYEgaIiZfYPW+l7ZiwN/DgwraSnvn2KDq1/Jm3TLn51L0dqK80FDksQxgSB737Zza/7jnDZeS0qLlyNkuKSmLrxLg7LKmat2hr0K80FG2ukNiYI/Fa91NzfoRzHkeDgg+TpXPr2laQs/omdfFojVpoz3rE7CGMCXFGR8tlqV/VS/Vo17zOfI8FBv5bXsv7Qm1x33k2WHAKIJQhjAtzyGlq9VMyZ5WTp9mk0yh/FpO9erRGLCBnvWIIwJsDNWrWVWjWwegl+W49henIqrSOuZ3DckzVmpTlTMUsQxgSwY72XOtTM6qVjK8217Ufvdk35JSeB93/3fo1Yac5UrOb9RRljvJb+82627z/CZef5d2GesngO0rugfVPmrNnG6Q0upl/bfn6MynjL7iCMCWCzf3BXL3WsOYPjynJh+6YALNyQ6+dIjLcsQRgToArd1UuODs2oVwOrl0o6PaYeraLr8I0liIBhCcKYAJOyKAVnlpP0Tbvc1UstAmZ08gVnxLIkcyf5hUX+DsV4wRKEMQGmeGnRl779mFoRYYTX+TFgRidf2L4pB44UsHLzHn+HYrxgCcKYAONIcDDtd+8zdeNdNGj2Idd9fHXAjE7u3S4GEayaKUD4NEGIyCARWS8iGSJybyn7RURecO9fJSJdPfY1FpEZIrJORNaKSC9fxmpMIImOOJ96+YNZvvs1bkm8JSCSA0DjulGcF9+IhRmWIAJBmS1bnm/WpVHV78rbLyLhwEu4lg3NBtJEZKaq/uhRbDDQ3v3oAUxwfwV4HpijqiPca1rXreBajAkZE5fOZH/EbO7pdT8T0ifgaOMImCRxQfumvPJ1Jvvy8mlYg2aeNScq7w7iGffjJWApMBHXutFLgRe8eO3uQIaqZqrqUWAaMLxEmeHAFHVZAjQWkRYi0hC4CHgdQFWPquoe7y/LmODlzHLy+to7uTD636QM/DepI1IDanTyBWfEUlikLNm409+hmAqUmSBU1aGqDuBnoKuqJqpqN+B8IMOL144HNns8z3Zv86ZMW2AH8IaIfC8ir4lIvdJOIiJjRCRdRNJ37NjhRVjGBLavMr8lOm8cV3ceDHDc0qKBoOvpjakbFW7VTAHAmzaIjqr6Q/ETVV0NdPHiOCllm3pZJgLoCkxQ1fOBg8AJbRjueCa6k1dibGysF2EZE9jOqn8ttYvOO27lOH8vLVoZtSLC6ZEQbQPmAoA3CWKt+xN8XxG5WEQmAWu9OC4baOXxvCWQ42WZbCBbVZe6t8/AlTCMCXnz1m6ndXRd2sXW93coJ+2C9rFk5h5ky57D/g7FlMObBHEDsAa4E/gb8KN7W0XSgPYikuBuZB4FzCxRZiYw2t2bqSewV1W3quo2YLOIdHCX6+8+rzEh7fDRQhZm5NKvYzNESrsBDwy/Tbth1cI1WYXj81U1T0ReAr7EVf2zXlXzvTiuQERuB+YC4cBkVV0jIje7978CzAaG4GrTOMTxieevwLvu5JKJd0nJmKC2ODOXIwVFx1UvBaL2zerTvGEtvtmQy++TWvs7HFOGChOEiPQF3gI24WozaCUi16nqgoqOVdXZuJKA57ZXPL5X4LYyjl0BJFZ0DmNCyVdrt1MvKpzuCdH+DuWUiAh9zmiKc912ioqUsLDAvRsKZt5UMT0DDFTVi1X1IuBS4D++DcsYU5KqMm/ddi5sH0utiHB/h3PKLmzflN2H8lmTs8/foZgyeJMgIlV1ffETVf0JsNEtxlSztVv3s3VvHv0CvHoJXBMOatQaAL7JcLVDBMqEg6HEmwSRLiKvu3sx9XX3Ylru68CMMcebt+5XABwdAj9BJMUl8ZfZ1xLbdAMLN+QeW5o0ECYcDCXeTCJ/C652gjtwtUEsAF72ZVDGmBN9tW47nVs1JrZBLX+HcsqKB/cNfe93ZG4exBfTvyB1ZGBMOBhKKryDUNUjwIvAQ8B44EX3NmNMNck9cIQVm/cExMpx3nIkOPhdhxvYFT6Voe2us+RQA1WYINy9mDbgShIvAz+JyEW+DcsY42n++h2oQr8gShDOLCezM6fQuGAU09dNDpi5pEKJ9WIyJgDMW/crzRvW4py4hv4OpUoUtzlMH5lK79hb6d7wkYCacDBUWC8mY2q4owVFLPgpl34dmwf06GlPaTlpxxY56p4QTc72trx75bSAmXAwVHjTSJ0uIq8Db7uf/wHrxWRMtUnbtIsDRwqCqv3Bc2LBHgnRvLFoEzGRXRnbp78fozIleXMHcQuuuZjuwDUf04/Azb4MyhjjGivgzHLy1drt1IoIc408DsKxAoltXKPCl2Xt8nMkpiSvejGp6rOqepWqXqmq/7FeTMb4XlJcEskzkvlgzRx6t4thyZYFQTlWoGn9WrSLrUfaJksQNY03vZj6iMgXIvKTiGQWP6ojOGNCmSPBwX8GvMXqvH+yO+JdkmckH6u3DzbdE2JI27SLwqKSS8YYf/Kmiul14FngAiDJ42GM8bG8/R1pUDCEjzNf4JbEW4IyOQB0T2jC/rwC1m/b7+9QjAdvEsReVf1MVber6s7ih88jM8YwbdUcDkV9xviLxjMhfULQdgPtnhADwLIse2upScpMECLSVUS6Ak4ReUpEehVvc283xvjQ7J++ZP6O+/hj+//wL8e/SB2RGrRjBeIb1yG+cR2WWTtEjVJeN9dnSjz3XJtBgX5VH44xpthHq7+m6dF7+XPS5cBv8xel5aQFZVVT94RovtmQi6oGzXiPQFdmglDVU/4LFJFBwPO4VpR7TVWfKLFf3PuH4FpR7npV/c69bxOwHygEClTVFg8yIaWZJBMTseVYN1BwJYlgTA7gShAffb+FrNyDtA3g9baDSZkJQkSuVdV3ROSu0var6rPlvbCIhAMvAQOAbCBNRGaqqufa0oOB9u5HD2CC+2sxh6rmenUlxgQRVWX++h30PqMpURHeNBUGvuJV8pZl7bIEUUOU95dXz/21QRmPinQHMlQ1U1WPAtOA4SXKDAemqMsSoLGItKjMBRgTjDK2H2DLnsNBsfaDt9o2rUfT+lHWDlGDlFfF9Kr76z9P8rXjgc0ez7M5/u6grDLxwFZc7Ryfi4gCr6rqxNJOIiJjgDEArVvb4ucmOMxf71plrW+HWD9HUn1EhKQ20TaiugYpr4rphfIOVNU7Knjt0lqZSo6CKa9MH1XNEZFmwBcisk5VF5QSx0RgIkBiYqKNsjFBYf5P2zmzeX3iGtfxdyjVqntCNJ+t3kbOnsMhd+01UXm9mE51Qr5soJXH85ZAjrdlVLX463YR+QhXldUJCcKYYHPgSAHLsnbxpz4J/g6l2iW5G+TTNu1ieJd4P0djyqtiesvzuYjUU9WDlXjtNKC9iCQAW4BRwDUlyswEbheRabiqn/aq6lYRqQeEqep+9/cDgX9V4tzGBKxvM3LJL1QuDqHqpWJntWhIg1oRLM2yBFETeDMXUy8R+RFY637eWUQqXJNaVQuA24G57mNTVXWNiNwsIsWzwc4GMoEMYBJwq3t7c2ChiKwElgGzVHVO5S7NmMDkXL+DelHhJJ4eXXHhIBMeJnRr04Q0a4eoEbxZD+I5XKvIzQRQ1ZXeLjmqqrNxJQHPba94fK/AbaUclwl09uYcxgQTVeXr9du5oH3odG8tqXtCNCnr17PzwBFi6tfydzghzau/QFXdXGJToQ9iMSbkbdh+gJy9efQNoe6tJfVIKG6H2O3nSIw3CWKziPQGVESiROTvuKubjDFVy7luOxBa3VtL6hTfmFoRYbY+RA3gTYK4GVc1UDyuXkdd+K2twBhTheav30HH0xrQolHodvF8bunTxDXfeNx4iGBcSS8QeJMgklT1D6raXFWbqeq1QLKvAzMm1OzPyyf9510h2XvJU1JcEmn7/kH6tm84cKQAZ5YzKFfSCwTeJIjxInJs5lYRGcuJU2YYY07Rooyd5Bcqfc8M3fYHcE1I+NjFk9ke+QS3fjIuqFfSq+m8SRDDgMdE5EIR+Teu8QrDfBuWMaHn65+2U79WBIltmvg7FL+7sftQGhYO4e01zwb1Sno1XYUJwj2b6jBcM7PGASNUNd/XgRkTClIWpeDMcqKqONft4IIzmrLwl69Dvr596ZZvOBT5GWfXuyGoV9Kr6cpbUW6/iOwTkf24BrKdCYwE9onIvuoK0JhglhSXRPKMZKZ8N4tt+/KIbboh5Ovbi9scrmn/Hwr2jOSdK6cG7Up6NV15U214M6W3MeYUFK8SN2zq7wiPuJSXf/iSGSNDu749LSeN1BGpHDl4Fs6Vy2ka2S2oV9KrycqbzbWjqq4ra/3p4pXfjDGnxpHgoFXUFazNf4PxSeND/k1wbJ+xAOQeOAK4Ju67+eLgXUmvJitvqo27cK2zUHJtarA1qY2pMp+s+5z1Bz7gkpa3MiF9Ao429mYI0LR+LdrG1iN90y64uJ2/wwlJ5VUxjXF/PeEvVUR6+jIoY0KFM8vJtR9dTezRe3ly4M3sLRxh3To9JJ0ezZw12ygqUsLCSls+xvjSyc4GllqlURgTotJy0rik+ePE1Umkc8vGx9ok0nLS/B1ajZCUEM3ew/lk7Djg71BC0skmCEvlxlSBu3vdw8bsNvQ9M5Zw9ydkR4LjWD18qEtyjwmxZUj942QThC3taUwV+P6X3ew5lE+/s0J79HRZWkfXJbZBLVc7hKl25fVi+oTSE4EAMT6LyJgQ8tW67YSHCRe2D+35l8oiInRvE21Tf/tJeb2Ynj7JfceIyCDgeSAceE1VnyixX9z7hwCHgOs9u8+KSDiQDmxR1aHenNOYQOJct52kNk1oVCfS36HUWIltmjDrh61s2XOY+MahO8utP5TXi+nrU3lh95v7S8AAXNOEp4nITFX90aPYYKC9+9EDmOD+WuxOXGtPNDyVWIypibbsOcy6bfu5f0hHf4dSoyW1cS0glL5pF/G2TnW18uWaht2BDFXNVNWjwDROnAV2ODBFXZYAjUWkBYCItAQuA17zYYzG+M089+JA/To293MkNdtZLRpSv1aELSDkB75MEPGA51Kl2e5t3pZ5DhgLFJV3EhEZIyLpIpK+Y8eOUwrYmOrkXLed1tF1aRdbz9+h1GjhYULX05uQbu0Q1c6XCaK0rrAlG71LLSMiQ4Htqrq8opOo6kRVTVTVxNhYa+gzgeHw0UIWZeTSr2MzXE1xpjxJpzdh/a/72XvIJpKuTifTiwkAVa1oTYhsoJXH85ZAjpdlRgDDRGQIUBtoKCLvuFezMybgLc7M5UhBEf06WvdWbyQlRKMKy3/ZZVVy1ai8O4incc3DlAUcBia5HweA1V68dhrQXkQSRCQKGAXMLFFmJjBaXHoCe1V1q6rep6otVbWN+7h5lhxMMPlq7XbqRoXTo220v0MJCJ1bNiYyXFiWZdVM1anCXkwi8oiqXuSx6xMRWVDRC6tqgYjcDszF1c11sqquEZGb3ftfAWbj6uKagaub6w0nfSXGBAjX4kDbueCMptSKCPd3OAGhTlQ458Y3sgFz1ay8cRDFYkWkrapmAohIAuBVZb+qzsaVBDy3veLxvQK3VfAa84H53pzPmECwbtt+cvbmcecl7f0dSkDp3iaaNxZtIi+/kNqRllirgzeN1P8HzBeR+SIyH3ACf/NlUMYEs+LurY4O1v5QGYltojlaWMSq7L3+DiVkVHgHoapzRKQ9UDyaZ52qHvFtWMYEr3nrttMpvhHNGtb2dygBJfF018R9aZt20T3B2m6qQ4V3ECJSF7gHuF1VVwKt3d1QjTFeSlmUgjPLya6DR/n+l904OjbDmeUkZVGKv0MLGE3qRdG+WX0bMFeNvKliegM4CvRyP88GHvVZRMYEoaS4JJJnJPPit/+jSKFho/Ukz0gmKS7J36EFjJRFKTSLzWD5pt0UFrl64FuS9S1vEkQ7VU0B8gFU9TC2HoQxlVK8ENDjS//C0bpTue/rP9mqcZWUFJfEp9lj2ZH/Heu37ceZ5bQk62Pe9GI6KiJ1cA+aE5F2gLVBGFNJvVpeRP38wWwNe5fxieMtOVSSI8HBpMve4eoPR3Hfl1v59teplmR9zJs7iIeBOUArEXkX+ArXHEnGmEp46dv/sUtm8cdz7mJC+gScWU5/hxRwRnYaRIvwYczM+i+3JN5iycHHKkwQqvo5cBVwPTAVSHSPTTDGeMmZ5WT8NzfRWh/gtSueInVEKskzki1JVNL8TfPZJZ9yGn+wJFsNvOnF9BXQQ1VnqeqnqporIhOrITZjgsbS7GXEFd7L8LMGEBURdqxNIi0nzd+hBYziNoe/J75CrcNX80z/Ny3J+pg3VUwJwDgRechjW6KP4jEmKF0U92cKDp/DoHNPO7bNkeBgbB+rrfVWWk4aqSNSuan75QBo3jmWZH3Mm0bqPUB/4AX3DK82aZ4xlTR39TZqR4Zx0Zk2Jf3J8kym8Y3rsCRzJ9f1dlg7hA95cwchqlqgqrcCHwALAZsjwBgvFRUpc9f8ysVnxlI3ypvPZKYiPdvGsCRzJ0VFZa5IYKqANwnCc3K9N3E1Vn/uo3iMCTors/ewbV/ecdVL5tT0ahfD7kP5rP91v79DCWplJggRaej+drqIRBc/cK0P8fdqic6YIDBnzTYiwsQWuqlCPd3raCzJ3OnnSIJbeXcQ77m/LgfS3V+Xezw3xlRAVZm7ehu9z2hKozqR/g4naLRsUpdW0XVYvNEShC+Vt2DQUPfXhOoLx5jgsv7X/WzaeYgxF7XzdyhBp1fbGOau+ZWiIiUszGb/8YXyqpi6lvfw5sVFZJCIrBeRDBG5t5T9IiIvuPevKn5dEaktIstEZKWIrBGRf578JRrjP3NWb0MEBpxt1UtVrWfbGPYezmfttn3+DiVoldel4ply9inQr7wXFpFw4CVgAK4ZYNNEZKaq/uhRbDDQ3v3oAUxwfz0C9FPVAyISCSwUkc9UdUlFF2RMTTJn9TaSTo8mtkEtf4cSdHq1iwFg8cadnBPXyM/RBKfyqphOtXNxdyDDY6nSacBwwDNBDAemuJceXSIijUWkhapuBQ64y0S6H9afzQSUTbkHWbdtP+OHnu3vUIJSi0Z1aBNTlyWZu7jxwrb+DicoedUpW0TOBc4Gji2BpapTKjgsHtjs8Twb191BRWXiga3uO5DlwBnAS6q6tIzYxgBjAFq3bl3htRhTXeau2QbApedY9ZKv9Gwbw6wftlJYpIRbO0SV82YupoeA/7ofDiAFGObFa5f22yp5F1BmGVUtVNUuQEuguztJnVhYdaKqJqpqYmysjVI1NcecNdvoFN+Ilk3q+juUoNWrXQz78wr4McfaIXzBm4FyI3BNtbFNVW8AOgPeVKhmA608nrcEcipbRlX3APOBQV6c0xi/Kl5adNvePL7/ZQ+Dzj3NVj3zoZ5tXe0QNh7CN7xJEIdVtQgocA+e2w54U+GXBrQXkQQRiQJGATNLlJkJjHb3ZuoJ7FXVrSISKyKNAdyLFV0CrPPukozxn+KlRZ9f+BEAjRr/ZKue+VDzhrVp27Qeiy1B+IQ3bRDp7jfrSbjaBA4Ayyo6SFULROR2YC4QDkxW1TUicrN7/yvAbGAIkAEcAm5wH94CeMvdDhEGpKrqp5W5MGP8oXga78HvXEWLRsP4vy9n26pnPtazXQwzV+RQUFhERLg3n3mNtypMEO5J+gBeEZE5QENVXeXNi6vqbFxJwHOb59xOCtxWynGrgPO9OYcxNc1Z0b2oc3QQm4qmML6nLS3qaz3bxvDe0l9Yk7OPzq0a+zucoOJVuhWR80RkGNAVOENErvJtWMYErhTnDPaFz+avieNs1bNqUDwvk1UzVT1vejFNBiYDvwMudz+G+jguYwLSvMx5vLjqdno1fpQXLnvClhatBs0a1KZdbD1rqPYBb9ogeqqqjfQxxguz1i8kOm8cfxno6gnuubSoVTX5Tq92MXz03RbyC4uItHaIKuPNT3KxiFiCMMYLjQtH0EA6M7RTi2PbbGlR30pZlEK9Bus5eLSQH7bsBbCuxVXEmzuIt3AliW245kgSXO3L5/k0MmMCTEFhER+vyKFfx2Y0qRfl73BCRlJcEiOnJxMRdhdLMjuwt3AFyTOSSR2R6u/QAp43CWIy8EfgB6DIt+EYE7i+ycgl98ARrura0t+hhBRHgoPpI1O59O0reeW7DLamzbSuxVXEmwTxi6qWHOBmjCnho++20LhuJI4OtmR7dXMkOOjd/Gq+3vYK9/Z5wJJDFfGmDWKdiLwnIleLyFXFD59HZkwA2Z+Xz9w127j8vDiiIqyRtLo5s5ys2D2dRvmjmJBmXYurijd/yXVwtT0MxLq5GlOqz37YxpGCIq7sGu/vUEKOM8vpbnN4n/iw6xgSn2Jdi6tIuVVM7qkuclX1nmqKx5iA9OH32SQ0rcf5NpK32qXlpB1rc7iw/XJWZkfy/u/ft67FVaDcBKGqhd4uL2pMqMrefYglmbu4a8CZiNiaBNXNswtxv47NmLNmG6fVvpCxfcpd9NJ4wZtG6hUiMhOYDhws3qiqH/osKmMCyMcrXDPUX3m+VS/5W9+OrjVh5q37lbPjGvo5msDnTRtENLAT1xrU1gZhjAdV5YPvsuneJppW0bYwkL81a1Cb81o2Yt667f4OJSh4M5vrDRWVMSbUpCxKISkuiSYR55O54yBjLmyLM8tJWk6ajZr2M0eHZrwwbwO7Dh4l2gYsnhJvJutrKSIfich2EflVRD4QERsJZEJa8cJAz379IVERYdRvtN4WBqoh+nVship8/ZPdRZwqb6qY3sC18lscEA984t5mTMhyJDh498ppvLfxLho1/5AbZl5jo3driE7xjWhavxbz1u3wdygBz5sEEauqb6hqgfvxJhDrzYuLyCARWS8iGSJybyn7RURecO9fVdxjSkRaiYhTRNaKyBoRubNSV2VMNTi0vyP18weTvus1bkm8xZJDDREWJjg6xPL1+u0UFNrsQKfCmwSRKyLXiki4+3EtrkbrcrnHULwEDAbOBq4uZVbYwUB792MMMMG9vQC4W1XPAnoCt9mMsqYmUVVS5n/AocjPePDCB21hoBqmX8dm7MsrYPnPu/0dSkDzJkH8CUgGtgFbgRHubRXpDmSoaqaqHgWmAcNLlBkOTFGXJUBjEWmhqltV9TsAVd0PrMVVvWVMjfDqkpks3TueuxMn8Ei/R2xhoBrmgvZNiQwX5q23dohTUWGCUNVfVHWYqsaqajNVvUJVf/biteOBzR7PsznxTb7CMiLSBtf61EtLO4mIjBGRdBFJ37HD6hxN9ZiyfB6n8wDjB4wEjl8YyPhfg9qRJLWJxmndXU9Jmd1cReQf5RynqvpIBa9d2pBSrUwZEakPfAD8TVX3lRHIRGAiQGJiYsnXN6bKZe8+xLYtA7nporbUjfrtX8iR4LB2iBqkX8dmPDprLZt3HbIxKiepvDuIg6U8AP4MjPPitbOBVh7PWwI53pYRkUhcyeFdG7VtapK3F/+MiDC6Vxt/h2LK0a+ja9p1p1UznbQyE4SqPlP8wPUJvQ5wA662hLZevHYa0F5EEkQkChiFq7usp5nAaHdvpp7AXlXdKq4JbV4H1qrqs5W/LGN849DRAqYu+4VB55xGfOM6/g7HlKNtbH3axNS1UdWnoNw2CBGJFpFHgVW4qqO6quo4Va3wJ66qBcDtwFxcjcypqrpGRG4WkZvdxWYDmUAGMAm41b29D65V7PqJyAr3Y8hJXJ8xVeqD77awL6+AP13Qxt+hGC84OjZj8cadHD5a6O9QAlJ5bRBPAVfhunvopKoHKvviqjobVxLw3PaKx/cK3FbKcQspvX3CGL8pKlLeWJRF55aN6Nq6ib/DMRVIWZRCkyZncqQgkm835tL/rOY2HUollXcHcTeu0dMPAjkiss/92C8ipTYYGxPMFmzYQeaOg/zpggSb1jsAJMUl8cjimyBqNV+t235sYSGbDsV7Zd5BqKqtm2iMh8mLNtGsQS0Gn9vC36EYLxR3PR7yzlVMWb2aSRlzbDqUSrIkYEwZUhalHBv4tuHX/Sz4aQc9z9rKc0uf9nNkxluOBAdD213HlqJ3GHbG9ZYcKskShDFlKJ6x1Znl5I1vN1EY+QPvZdxlVRQBxJnlZH72uzQtupqpa163ke6V5M2KcsaEpOIqipHTk+HgAA7VmsOskR/Yp9AAUdzmkDoylc+WR5O6qjPJ013P7XfoHbuDMKYcjgQHXaJHsjNsKjd0vsneWAJIWk7asTaHEd1awdFzub3LSzYdSiVYgjCmHDNWz2F+9rt0j7mJ1HWTrYoigIztM/ZYQk9q04TW0XX56ZfW1sW1EixBGFMGZ5aT0f+7htMK7mPGNc/ZjK0BTEQY0a0l327cSfbuQ/4OJ2BYgjCmDLPXL6TR4Xu4uecwWkXXtRlbA9yV58ejCh99t8XfoQQMSxDGlGHXtkE0jezKbY4zjm1zJDisiiJAtYquS6+2Mcz4LhvXJA6mIpYgjCnFksydfLVuO7f2PYMm9aL8HY6pIiO6teTnnYdIt5XmvGIJwpgSVJXHZ6+lRaPa3NCnjb/DMVVocKfTqBcVzoz0bH+HEhAsQRhTwqwftrIyey93DTiT2pHh/g7HVKG6UREM6dSCWT9s5dDRAn+HU+NZgjDGw9GCIp6au56OpzXgqq4t/R2O8YER3Vpy4EgBc9ds83coNZ4lCBPyPOdcmrrsF37eeYiBXXN5ZvFTfo7M+EJSm2haR9flg+XWm6kiliBMyCuec2nW+i94/qsNtI3fxONL/2JzLgWpsDDhqq7xLNqYy5Y9h/0dTo3m0wQhIoNEZL2IZIjIvaXsFxF5wb1/lYh09dg3WUS2i8hqX8ZoTPH4huQZvyfzyGRWHHrIpoUOYimLUmjRbKN7TISrsdqZ5SRlUYqfI6t5fJYgRCQceAkYDJwNXC0iZ5coNhho736MASZ47HsTGOSr+IzxFFXYicjDl7I3chq3d7/VkkMQS4pL4va5o0mIz2LG8mzmZc6zhYTK4Ms7iO5AhqpmqupRYBowvESZ4cAUdVkCNBaRFgCqugDY5cP4jAFg98Gj3DjtLQ5Ffsa9fR5gQvoEm04jiBXfMS4/8BAr9k7kqtRku2Msgy8TRDyw2eN5tntbZcuUS0TGiEi6iKTv2LHjpAI1oUtVue69N/ip4F+8MuRtHr/kUZtzKQQ4Ehzc3v0W9kZOo1nYUPq26evvkGokXyaI0hbtLTm+3Zsy5VLViaqaqKqJsbGxlTnUGKYvz+abn5dwR5eX+FPS5QA251IIcGY5mfTdq/yu/Z1sPPQRTzo/8HdINZIvFwzKBlp5PG8J5JxEGWN8YlPuQR6euYZLW4/h6WE9j9vnSHBYlUOQOraQ0IhULjq9L92fbsM/vrmRxNObcEm7/v4Or0bx5R1EGtBeRBJEJAoYBcwsUWYmMNrdm6knsFdVt/owJhPCPMc75BcW8bf3V5AXtooz2s4jPKy0m1kTjDwXEgoPEx67bBTRR8bxRtpX/g6txvHZHYSqFojI7cBcIByYrKprRORm9/5XgNnAECADOATcUHy8iEwF+gJNRSQbeEhVX/dVvCb4FY93SB2RyvcZ8SzZsoDDDZ6mf7sZ/g7NVKOSs/EOPLs5PeMvIjMrj7z8QptexYME07S3iYmJmp6e7u8wTA3mzHJy1fsj4eAAjtSaw6w/fGhVSYZvN+ZyzaSlPHjZWdx4YVt/h1OtRGS5qiaWts9GUpuQcnr9JOocHcSeiGn8tYeNdzAuvds15cL2TXnJmcH+vHx/h1NjWIIwIeOXnYe47NVX2aGfclu3cUxeMdG6sppjxl7akd2H8pn0TZa/Q6kxLEGYkJCz5zBDXp3AhsJHeG3oO7w49Akb72CO06llI+Jafs7z33xE7oEjx7aH8jQcliBM0Nu+L49rJi0h98haXh38Ntd1GwrYeAdzolsvuJTssMe4++P3gN+6xIbqNBy+HAdhjF+kLEohKS4JR4KDnQeO8IfXlpK1fxnJ3Vvy5+6XH1fWxjsYT9d2GcKslc/w7oa/0WT2T7y75rWQnobD7iBM0CnuzvrJus+59vVlrNu9mH11n+LKcy72d2gmADx7xR9oopfxfNoT3HT+X0I2OYAlCBOEHAkOXrx0CiNSk1m2cwL76z3Nh7+fHtL/6MZ763YtIb/2XBoVjOI/i1/iq8x5/g7JbyxBmKDzycocHv0gjEZFl7ErfCp3WHdW46XiNoePRs0gZcCjNDo8luFTR4RsRwZLECZgeU6dAXDoaAGj3nqN0dMfoHGTnyis8znjLxpv03cbr3lOw3HThW25vttl1D94D5OXfenv0PzCGqlNwPKcOqN5rW5c+/brrDz8MEPa3cKS3H8xI9n1j+5o4zhWzu4kTHk8p+EQER654lx+2XWIJat2sTRpJz3axvgxuupnU22YgPZl5jyunDaSqLxL2Rc+m2f7v8lh+elYL6ZiziwnaTlpJ8zDY0xF9h7Kp+vztyNH2/HlbbfRpmk9IHj+psqbasPuIEyN5tlltZgzy8mS7GWcUecaXp4fRvihgeyKnMrdPe/jrxdcUerrWHdWc7Ia1Y3k0cHD+ePH13Dl65HM/+vtfPfrwmN3pcHM2iBMjVZcjVTchjA340uGTx3Bu99EcM+MVRxiJVrvcx688EHeWjXJ2hqMT1zdZQgvDHiTVYcfpvNzf+aq90eGRJWlJQhTI5RscIbfbuFTR6Qycnoyw6fcwdB3R1D3wN9p16gnfx2cR6Y+yv9GzeCRfo/Y1BnGp27tfQWjO93Ez/lvowcGMPe7GPYeCu6J/SxBmBqh5J2CM8vJiOnJ7N7Tisnz6lK0/xJmZv2XM+peSer1N/K/W3tTGLHxuE9xNnWG8SVnlpNPM9/i3j4PkF9nLlO+m8U5T9/Ko19Ox7MtN5jmbrIEYXyirDuCIe8OKfNO4a3hU7nq/ZEMnnw7g96+ksg9dzH1m4aszv2Wo7XnclPnv5PLpxRGrkZEGNtn7Am3+I4ER8A3Gpqax3OZ0scveZRPr/mAvIbP0KhObR5aeCODJkwgc8eBY+U27tpY6t95oCUOnyYIERkkIutFJENE7i1lv4jIC+79q0Skq7fHVoXKvolVdnvKohSfn8Of5y7vHBt3bTzhjiB5RjKXtL2E5BnJvL/qM5zrt3PnR+8w5J2reG9hJLe9cQQ9MIA5m1+iTe0r+cfAkTz6eyE77DE+ueYDJl7xlFUjGb/wHB8Brg8iHyRP54+9WnJ3twl8tf0+uj5/E5dOuZIhcSmc2XAgI6ef+PefFJdULf+XVUZVffLAtczoRqAtEAWsBM4uUWYI8BkgQE9gqbfHlvbo1q2bVsa8zHnaNKWpzsucd9zzZ759pkq2z8uc5/Nz+OPcX238Sps+2VQ/Wfe5frh6jkY/2VTf/n6Wrty8W1/45iNt+Fi0jv9smt7z8bta79Em2uflv2jtfzXWHk8/p90e+Vyb3/eYhj3UUBs9MErDHmqo5z72tP5lSrrePmOKNnosWu/67L5j53ty4ZPHzuv5e3ty4ZOV+l0b40v/99l9ysPoec/+STs8OFtPH/epxj3wmNb6Z2Md8Nqt2uDf0frkvBmavmmXTl0xW2Oe9O3/fmUA6VrGe6rPxkGISC/gYVW91P38PndCetyjzKvAfFWd6n6+Htc61G0qOrY0JzMOwpnlZNDbVxIfOZzs/I/pVPthoiPOZ1fB9/yQ9zAtT3E7UGWvVbw9PnI4W/I/5lyP7avzHiY+Yhhb8mdybu2HaOI+9+6C71md90/iIoeRkz+Tc2o/ROPw89ld+D0/5v2TFhHDyCmYSYfI8TQKO5+9Rd+zPv8RmoUNZXvRp5zO/dTVzuzXlWyWf9Oo6DL2hM0i9ui91C46D4C8sFXsiHqCBgVD2B8x+7h9eyLeYW/kNNrVuo4Lm9/KaY1qc1rDOiz49UU+zHiBe3rdT8rAfx93C+9IcJzw3Jiaqvhv9ZbEW5iQPoF3rphKfenCgg07eH3lk2w88haN8kfRuODaY8fkha0iN+oJohnKbplF27AHiYk4n/2sYP3RR479vxb/L+8uXMHqct4Tzqh7FTv0k5P6f/HXOIh4YLPH82yghxdl4r08FgARGQOMAWjdunWlg3QkOOgSncyynZPo1uRGusdc6N5zIbJzJMt3v3aK26vytX7bntjkRrrHXOTefhHhO0eSvvs1EqNvpHuM56ylF7v27XLt69m0r+vnRl+icteStmsSPWJuondsfxAIk0tosP0nluROpE+zv+BoMYTwMCFMWuHcmsX8ra8woOWtDG87kloRYdSKCCMqohMfZPzK9J+e50+d7uaunrfSsHYkK7Yv5E+ffMn4RNd0F6P7XocjoTvOLCePfv/esWkwBre/pNRb+OIGZ0sQpqYq+UHGc9R+z47w7OpZ3HfBA7y6/FX+4biGhIY92HngKLkHOpL602acORO4oPlfuKDZpRQUFlFY1J8629eTtmvSb//LCnARYRW8J4y/aHzV/6+UdWtxqg9gJPCax/M/Av8tUWYWcIHH86+Abt4cW9qjslVMqr/dlo2fN77U27VT3V4d5/DnucvaV9nqrcreFhtTE5RVBTpm5phy/86r4//SW5RTxeTLBNELmOvx/D7gvhJlXgWu9ni+HmjhzbGlPawNonrPXd45xswcU+o/zuB3Blubggl65bWdVcf/ZWX4K0FEAJlAAr81NJ9TosxlHN9IvczbY0t7VDZBlPVLLOtNrLLbn1z4pM/P4c9zl3cOe8M3pnTV8X9ZGeUlCJ9O1iciQ4DncPVKmqyq/xaRmwFU9RUREeBFYBBwCLhBVdPLOrai89lkfcYYUznlNVLbbK7GGBPCyksQNpLaGGNMqSxBGGOMKZUlCGOMMaWyBGGMMaZUQdVILSI7gJ9P8vCmQG4VhhMI7JqDX6hdL9g1V9bpqhpb2o6gShCnQkTSy2rJD1Z2zcEv1K4X7JqrklUxGWOMKZUlCGOMMaWyBPGbif4OwA/smoNfqF0v2DVXGWuDMMYYUyq7gzDGGFMqSxDGGGNKFVIJQkQGich6EckQkXtL2S8i8oJ7/yoR6eqPOKuSF9f8B/e1rhKRb0Wksz/irEoVXbNHuSQRKRSREdUZny94c80i0ldEVojIGhH5urpjrGpe/G03EpFPRGSl+5pv8EecVUVEJovIdhFZXcb+qn//Kmse8GB74Jo2fCPQlt/WmDi7RJkhHL8+xVJ/x10N19wbaOL+fnAoXLNHuXnAbGCEv+Ouht9zY+BHoLX7eTN/x10N13w/8KT7+1hgFxDl79hP4ZovAroCq8vYX+XvX6F0B9EdyFDVTFU9CkwDhpcoMxyYoi5LgMYi0qK6A61CFV6zqn6rqrvdT5cALas5xqrmze8Z4K/AB8D26gzOR7y55muAD1X1FwBVDfTr9uaaFWjgXnemPq4EUVC9YVYdVV2A6xrKUuXvX6GUIOKBzR7Ps93bKlsmkFT2ev6M6xNIIKvwmkUkHrgSeKUa4/Ilb37PZwJNRGS+iCwXkdHVFp1veHPNLwJnATnAD8CdqlpUPeH5RZW/f0WcUjiBRUrZVrKPrzdlAonX1yMiDlwJ4gKfRuR73lzzc8A4VS10fbgMeN5ccwTQDegP1AEWi8gSVf3J18H5iDfXfCmwAugHtAO+EJFvVHWfj2Pzlyp//wqlBJENtPJ43hLXJ4vKlgkkXl2PiJwHvAYMVtWd1RSbr3hzzYnANHdyaAoMEZECVf1ftURY9bz9285V1YPAQRFZAHQGAjVBeHPNNwBPqKuCPkNEsoCOwLLqCbHaVfn7VyhVMaUB7UUkQUSigFHAzBJlZgKj3b0BegJ7VXVrdQdahSq8ZhFpDXwI/DGAP016qvCaVTVBVduoahtgBnBrACcH8O5v+2PgQhGJEJG6QA9gbTXHWZW8ueZfcN0xISLNgQ5AZrVGWb2q/P0rZO4gVLVARG4H5uLqATFZVdeIyM3u/a/g6tEyBMgADuH6BBKwvLzmfwAxwMvuT9QFGsAzYXp5zUHFm2tW1bUiMgdYBRQBr6lqqd0lA4GXv+dHgDdF5Adc1S/jVDVgpwEXkalAX6CpiGQDDwGR4Lv3L5tqwxhjTKlCqYrJGGNMJViCMMYYUypLEMYYY0plCcIYY0ypLEEYY4wplSUIY6qBez4gRORhz+fG1GSWIIypHl1E5AUgWkSuAP7t53iMqVDIDJQzprqISBtgDrAUOB/XdBajgZeBxUCkqt7itwCN8ZINlDOmirkTRBZwgaouEpHJuNdiwDWi9ysgUVUf9F+UxlTMqpiM8Y3NqrrI/f07wIWqegew0z3v03i/RWaMlyxBGOMbJW/NiwBU9WH3V7t1NzWeJQhjfKO1iPRyf381sNCfwRhzMixBGOMba4HrRGQVEA1M8HM8xlSa9WIyxjeKVPVmfwdhzKmwOwhjjDGlsm6uxhhjSmV3EMYYY0plCcIYY0ypLEEYY4wplSUIY4wxpbIEYYwxplT/D764MH2cw7ZlAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "plt.plot(x, likelihoods/likelihoods.sum(), label = 'Binomial')\n", | |
| "plt.plot(x, β_distribution, marker = 'x', color = 'green', linestyle='', label = 'Beta')\n", | |
| "plt.xlabel('p*')\n", | |
| "plt.ylabel('Normalized Likelihood')\n", | |
| "plt.legend();" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "id": "96a09f54", | |
| "metadata": {}, | |
| "source": [ | |
| "We can use this Beta distribution to determine the expectation value (mean) of $p$.\n", | |
| "Note that because the Beta distribution is slighly asymmetric, the expectation value will be slightly closer to 0.5 than the mode (peak) value." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 7, | |
| "id": "49d47cf8", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "0.6363636363636361" | |
| ] | |
| }, | |
| "execution_count": 7, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "β.expect()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "id": "904ef7b8", | |
| "metadata": {}, | |
| "source": [ | |
| "We can also calculate the 90% credible interval of $p$" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 8, | |
| "id": "1b81dd37", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "(0.4640640923907075, 0.7942501012493469)" | |
| ] | |
| }, | |
| "execution_count": 8, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "β.interval(.90)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 9, | |
| "id": "ffca90aa", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA1MklEQVR4nO3deXzU9Z348dc7d0JIQg6uBEgCEUS5wi14gIpgVdS6XWy7VtuuZdWt3bZW7WN3bbe7bXdt+6uuB4vWa626ttUWKwothyc3cpNACAECgRxAyEHOef/+mImOYZJMIJNvZub9fDzmkfl+v5/vzPtLwrzn+zlFVTHGGGPai3A6AGOMMX2TJQhjjDE+WYIwxhjjkyUIY4wxPlmCMMYY41OU0wH0pPT0dM3OznY6DBMGCqsKARidNtrhSEwb+52cny1btlSqaoavYyGVILKzs9m8ebPTYZgwcNULVwGw9s61jsZhPmO/k/MjIoc6OmZVTMYYY3yyBGGMMcYnSxDGGGN8CmgbhIjMBx4DIoFnVfXn7Y6L5/j1QD1wp6pu9Rz7J+CbgAI7gbtUtSGQ8Rpjwk9zczOlpaU0NIT2x0tcXBxZWVlER0f7fU7AEoSIRAJPAtcCpcAmEVmmqnu8ii0A8jyP6cDTwHQRyQS+DYxV1bMi8jqwCHghUPEaY8JTaWkp/fv3Jzs7G/d31tCjqlRVVVFaWkpOTo7f5wWyimkaUKSqxaraBLwGLGxXZiHwkrqtB1JEZIjnWBQQLyJRQAJwLICxGmPCVENDA2lpaSGbHABEhLS0tG7fJQUyQWQCR7y2Sz37uiyjqkeBXwCHgTKgWlVX+noTEblbRDaLyOaKiooeC94YEz5COTm0OZ9rDGSC8BVN+7nFfZYRkQG47y5ygKFAPxH5qq83UdWlqjpFVadkZPgc62GMOQ/FFbW8vP4QlbWNTodiHBLIBFEKDPPazuLcaqKOylwDHFTVClVtBt4ALgtgrMYYoLGllWXbj7Fo6Trm/vI9/vmPu/jC4x+wqeSk06GFtMjISCZOnMiECRPIz8/n448/7rT86dOneeqppwIeVyATxCYgT0RyRCQGdyPzsnZllgF3iNsM3FVJZbirlmaISIKnp9PVwN4AxmpMWDtQUct/vL2HGT9dxbdf/YSjp8/ywHWjeeXvp5MQE8WipetZ+v4BbIGxwIiPj2fbtm1s376dn/3sZzz88MOdlg/6BKGqLcB9wArcH+6vq+puEVksIos9xZYDxUAR8Axwj+fcDcDvga24u7hGAEsDFasx4ez/Nh3m6l++x/MflTAjN42Xvj6N974/h3vnjOKyken86b5ZzBs7iJ8uL+Bb/7uF6rPNTocc0s6cOcOAAQM+3X700UeZOnUq48eP55FHHgHgoYce4sCBA0ycOJEHHniA2tparr76avLz8xk3bhx/+tOfeiSWgI6DUNXluJOA974lXs8VuLeDcx8BHglkfMaEu0NVdfxo2R4uG5nGrxdNZGD/uHPKJMVF89RX8nnuoxJ+tnwvN/73hzz1lXwuzUx2IOLA+vFbu9lz7EyPvubYoUk8cuMlnZY5e/YsEydOpKGhgbKyMlavXg3AypUr2b9/Pxs3bkRVuemmm3j//ff5+c9/zq5du9i2bRsALS0tvPnmmyQlJVFZWcmMGTO46aabLrjx3UZSGxOmWl3K917fTlSk8MsvTfCZHNqICN+YncP/fWsGTS0ubn36Y9YUlPditKGtrYqpoKCAd999lzvuuANVZeXKlaxcuZJJkyaRn59PQUEB+/fvP+d8VeWHP/wh48eP55prruHo0aOcOHHiguMKqdlcjTH++82HxWw+dIpffWkCQ5Lj/Tpn8ohU3v72bG5/Zj2PLNvN7Lx0oiND53tmV9/0e8PMmTOprKykoqICVeXhhx/mW9/61ufKlJSUfG77t7/9LRUVFWzZsoXo6Giys7N7ZGR46PxmjTF+23eihl+s2Me8sYO4ZVL74UmdS0uM5aEFYzh8sp7XNx/p+gTTLQUFBbS2tpKWlsZ1113Hc889R21tLQBHjx6lvLyc/v37U1NT8+k51dXVDBw4kOjoaNasWcOhQx3O4N0tdgdhTJhpbnXx3de3kRgXxU9vHXde9dRzRg9k8ogBPL5qP1/MzyIuOjIAkYaPtjYIcFcXvfjii0RGRjJv3jz27t3LzJkzAUhMTOTll19m5MiRzJo1i0svvZQFCxbw4IMPcuONNzJlyhQmTpzImDFjeiQuSxDGhJknVhex6+gZlnw1n/TE2PN6DRHh+/NGc/sz63l5/SG+eXluD0cZXlpbWzs8dv/993P//fefs/+VV1753Pa6det6PC6rYjImjOwsreaJNUXcMimT+ZcO6fqETswcmcbleek8tfYAtY0tPRSh6UssQRgTJhqaW/nu69vISIzlRz3UGPv9eaM5WdfEcx8e7JHXM32LJQhjwsRTa4rYX17Lf942nuQE/9cE6MyEYSnMGzuIZ94v5lRdU4+8phPCYYT4+VyjJQhjwsDZplZeXHeI+ZcM5sqLenZSy+/NG01tUwtL3j/Qo6/bW+Li4qiqqgrpJNG2HkRcXMdjXXyxRmpjwsCfth2l+mwzd87K7vHXHj24PwsnDOXFj0v4xqwcBiZ170PIaVlZWZSWlhLqywW0rSjXHZYgjAlxqsoLH5cwZnB/puekBuQ9vnPNRfx5RxlPrCni3xZeGpD3CJTo6OhurbIWTqyKyZgQt/HgSQqO1/C1ywK3pGZ2ej++NHUYr248zJGT9QF5D9P7LEEYE+JeWneI5Phobp7YvRHT3fWPc0chIjyxuiig72N6jyUIY0JYWfVZ3t19nL+dOoz4mMCOdh6SHM8tEzP5845jNDR3PPDLBA9LEMaEsN+uP4xLlb+bMaJX3u/GCUOpa2q1mV5DhCUIY0JUQ3Mrr248zNVjBjEsNaFX3nNGbipp/WL4846yXnk/E1gBTRAiMl9ECkWkSEQe8nFcRORxz/EdIpLv2T9aRLZ5Pc6IyHcCGasxoebtHWVU1TVx52XZvfaeUZERLBg3mFUFJ6iz6TeCXsAShIhEAk8CC4CxwO0iMrZdsQVAnudxN/A0gKoWqupEVZ0ITAbqgTcDFasxoUZVeXFdCaMGJjJrVFqvvveN44fS0OxilVUzBb1A3kFMA4pUtVhVm4DXgIXtyiwEXlK39UCKiLSfQexq4ICq9swE58aEgU+OnGZHaTVfmzkiYF1bOzI1O5VBSbG8tf1Yr76v6XmBTBCZgPdqIqWefd0tswh4taM3EZG7RWSziGwO9ZGQxvjrxY9L6B8bxa353Rs52xMiIoTrxw3hvcIKzjQ09/r7m54TyATh62tL+8lOOi0jIjHATcDvOnoTVV2qqlNUdUpGRs/OMWNMMCqvaWD5zjJum5JFv1hnJku4YfxQmlpd/GX3ha+LbJwTyARRCgzz2s4C2t9zdlVmAbBVVe2vzBg/vbLhMM2tyh0zsx2LIX94Cpkp8fx5h1UzBbNAJohNQJ6I5HjuBBYBy9qVWQbc4enNNAOoVlXv/nG300n1kjHm81wu5fVNR7jyogxy0vs5FoeIcMP4IXywvzKopwEPdwFLEKraAtwHrAD2Aq+r6m4RWSwiiz3FlgPFQBHwDHBP2/kikgBcC7wRqBiNCTWfHDnNseoGbpkU2Gk1/HHD+KG0uJQVu487HYo5TwGtoFTV5biTgPe+JV7PFbi3g3Prgd7tn2dMkHt7RxkxURFcffFAp0Ph0swkRqQl8OcdZSyaNtzpcMx5sJHUxoQIl0tZvrOMKy/KoH9cz6wYdyHaqpk+PlBJZW2j0+GY82AJwpgQ8cmRUxw/08AXxrUfSuScGycMxaXwzk6beiMYWYIwJkT8uQ9VL7UZPag/owYm8pbNzRSULEEYEwLaqpeu6iPVS23aqpk2lZzkxJkGp8Mx3WQJwpgQsPXwKU6caeQL4/tO9VKbG8YPRdXdgG6CiyUIY0LAZ9VLg5wO5RyjBiZy8ZAk3rJBc0HHEoQxQc7lUt7Z5a5eSnRoao2uLLh0MJ8cPk2V9WYKKpYgjAlyW/pw9VKby/PSAfjoQJXDkZjusARhTJB7e0cZsX20eqnN+KwUkuKi+HC/zbgcTCxBGBPEPu29NLrvVi8BREYIl41M58P9lbgnUDDBwBKEMUFs86FTlNc08oXxQ50OpUuz89I5Vt1AcWWd06EYP1mCMCaILd/pqV4a03cGx3WkrR3iw/2VDkdi/GUJwpgg1eqpXpozeqBjCwN1x4i0fgxLjecDSxBBwxKEMUFqc8lJT/VS3+291N7sURmsL66iudXldCjGD5YgjAlSbdVLc4OgeqnN5Xnp1Da2sP3IaadDMX6wBGFMEHK5lOW7jjN3THBUL7W5bGQaIlg1U5AIaIIQkfkiUigiRSLykI/jIiKPe47vEJF8r2MpIvJ7ESkQkb0iMjOQsRoTTLaXnqaippH5lw52OpRuSUmIYXxmMh8WWYIIBh1+9fD+sPZFVbd2dlxEIoEncS8bWgpsEpFlqrrHq9gCIM/zmA487fkJ8Bjwrqre5lnTOqGLazEmbKwuKCcyQrjyogynQ+m22XnpLHmvmDMNzST1oZlnzbk6uzf9pednHDAF2A4IMB7YAMzu4rWnAUWqWgwgIq8BCwHvBLEQeMmz9Oh6z13DEKAOuAK4E0BVmwBb+dwYj1V7y5k8YgApCTFOh9Jts0dl8OSaA6w/UMW8S4LrDijcdFjFpKpzVHUOcAjIV9UpqjoZmAQU+fHamcARr+1Szz5/yuQCFcDzIvKJiDwrIv18vYmI3C0im0Vkc0WFDeM3oa+s+ix7ys4ExdgHX/JHpJAQE2nVTEHAnzaIMaq6s21DVXcBE/04T3zsaz/GvqMyUUA+8LSqTsJ9R3FOG4YnnqWe5DUlIyP4breN6a7VBeUAfWrluO6IjYpkek6qDZgLAv4kiL2eb/BXiciVIvIMsNeP80qBYV7bWUD7CeE7KlMKlKrqBs/+3+NOGMaEvdV7yxmemsDIjESnQzlvs/MyKK6s4+jps06HYjrhT4K4C9gN3A98B3cbwl1+nLcJyBORHE8j8yJgWbsyy4A7PL2ZZgDVqlqmqseBIyIy2lPuaj7fdmFMWDrb1MqHRZXMHTMQEV834MHhs2k3rFq4L+uyA7WqNojIk8BfcVf/FKpqsx/ntYjIfcAKIBJ4TlV3i8hiz/ElwHLgetxtGvV8PvH8I/BbT3Ipxr+kZExIW1dcSWOLK2irl9rkDUxkUFIsH+yv5G+nDnc6HNOBLhOEiFwFvAiU4G4zGCYiX1PV97s6V1WX404C3vuWeD1X4N4Ozt2Gu/eUMcZj1d5y+sVEMi0n1elQLoiIMGtUOmsKynG5lIiI4L0bCmX+VDH9Epinqleq6hXAdcD/C2xYxpj2VJXVBeVcnpdBbFSk0+FcsMvz0jlV38zuY2ecDsV0wJ8EEa2qhW0bqroPsNEtxvSyvWU1lFU3MDfIq5fazBrlbof4oMjaIfoqfxLEZhH5jacX01WeXkxbAh2YMebzVhecAGDO6NBIEAP7xzFmcH/r7tqH+ZMg/gF3L6Zv4+7JtAdYHMigjDHnWlVQzoRhKWT0j3U6lB4ze1Q6m0tOcbap1elQjA9dJghVbQSeAB4B/gV4wrPPGNNLKmsb2XbkdNCOnu7IrLx0mlpdfHL4lNOhGB+6TBCeXkz7cSeJp4B9InJFYMMyxnhbW1iBKkG19oM/Jo8YQITAhoMnnQ7F+ODPRPJtvZgKAUTkIuBVYHIgAzPGfGZ1wQkGJcVyydAkp0PpUUlx0Vw8JImNliD6JOvFZEwf19Ti4v19lcwdMyioR093ZFpOKp8cOUVTiy1D2tdYLyZj+rhNJSepbWwJufaHNtNzUmlodrHzaLXToZh2rBeTMX3cqr3lxEZFfDpuINRMyXaPCrdqpr7Hn7mYGoFfeR7GmF6kqqwqOMFlI9OIjwn+0dO+pCfGMjKjH5tKTvIPjHQ6HOPFn15Ms0TkLyKyT0SK2x69EZwx4a64so5DVfXMvXiQ06EE1LScNDaVnKTV1X7JGOMkf6qYfoP77mE2MNXrYYwJsDWexYHmjA7txbCm5QygpqGFwuM1TodivPjTzbVaVd8JeCTGmHOsLawgb2AiWQMSnA4loKblpAGw8WAVY0OsK28w6/AOQkTyRSQfWCMij4rIzLZ9nv3GmACqa2xh48GTXBXidw8AmSnxZKbEs7HEGqr7ks7uIH7Zbtt7bQYF5vZ8OMaYNh8fqKKp1RUyk/N1ZVpOKh/sr0RVQ3K8RzDqMEGo6pwLfXERmQ88hntFuWdV9eftjovn+PW4V5S7U1W3eo6VADVAK9CiqrZ4kAkrawvdiwO1dQMNddNyUnnzk6McrKwjN4jX2w4lHSYIEfmqqr4sIt/1dVxVO+32KiKRwJPAtUApsElElqmq99rSC4A8z2M68LTnZ5s5qmpzAZuwo6qsLazgslHpxET505ck+LWtkrfx4ElLEH1EZ395/Tw/+3fw6Mo0oEhVi1W1CXgNWNiuzELgJXVbD6SIyJDuXIAxoaiovJajp8+GTfUSQG56P9ITY6wdog/prIrpfzw/f3yer50JHPHaLuXzdwcdlckEynC3c6wUEQX+R1WX+noTEbkbuBtg+HBb/NyEhrWF7lXWwqGBuo2IMDU71UZU9yGdVTE93tmJqvrtLl7bVytT+1EwnZWZparHRGQg8BcRKVDV933EsRRYCjBlyhQbZWNCwtp95Vw0KJGhKfFOh9KrpuWk8s6u4xw7fTbsrr0v6qwX04VOyFcKDPPazgKO+VtGVdt+lovIm7irrM5JEMaEmlpP99avz8pxOpReN9XTIL+p5CQLJ2Y6HI3prIrpRe9tEemnqnXdeO1NQJ6I5ABHgUXAl9uVWQbcJyKv4a5+qlbVMhHpB0Soao3n+Tzg37rx3sYErY+LKmluVa4Mo+qlNhcPSaJ/bBQbDlqC6Av8mYtppojsAfZ6tieIyFNdnaeqLcB9wArPua+r6m4RWSwibbPBLgeKgSLgGeAez/5BwIcish3YCLytqu9279KMCU5rCivc3VtHhEf3Vm+REcLk7AFssnaIPsGfqTZ+DVyH+9s+qrrd3yVHVXU57iTgvW+J13MF7vVxXjEwwZ/3MCaUqCrvFZYzOy98ure2Ny0nlf8qLKSqtpG0xFinwwlrfv0FquqRdrtaAxCLMWFvf3ktx6obuCqMure2Nz2nrR3ilMORGH8SxBERuQxQEYkRke/jqW4yxvSsttlbw6l7a3vjMlOIjYpgk42HcJw/CWIx7mqgTNy9jibyWVuBMaYHrS2sYMzg/gxJDt8unjFREUwanmLjIfoAfxLEVFX9iqoOUtWBqvpV4EuBDsyYcFPT0MzmQyfDsvdSe9Ny0th9rJraxhanQwlr/iSIfxGRT2duFZEfcO6UGcaYC/RRURXNrcpVF4Vv+0ObadmpuBQ2WzWTo/xJEDcBPxWRy0XkP3CPV7gpsGEZE37e21dOYmwUU7IHOB2K4/JHpBAZIWy2hmpHddnNVVUrReQm4K+4R1ff5umeaozpIarKmoIKZo9KJzoyPLu3ekuIieLSoUnWUO2wzlaUqxGRMyJSg3sg20XA3wBnRORMbwVoTDgoPFHD8TMNYd17qb0p2alsO3KaxhbrVe+UDhOEqvZX1SSvn3Gqmti23ZtBGhPq1hS4Z2+1BurPTM1OpbHFxa6j9n3UKZ3N5jpGVQs6Wn+6beU3Y8yFW1NQztghSWHdvbW9traYTSUnmTzC2mWc0FkbxHdxr7PQfm1qsDWpjekxp+ub2HL4FP9w5UinQ+lT0hNjyc3o5+7JZP82juhsNte7PT/PWZtaRGYEMihjwsl7+ypodSlzL7bure1NHZHKu7uP43IpERG+lo8xgXS+3SVe79EojAljawrKSe0Xw4SsFKdD6XOm5qRSfbaZoopap0MJS+ebICyVG9MDWl3K2n0VXHVRBpH2DfkcUz3tEDbthjPON0HYOAhjesAnh09xur7Zqpc6MDw1gYz+sTai2iGd9WJ6C9+JQIC0gEVkTBhZVVBOZIRweZ51b/VFRJiWnWpTfzuks15MvzjPY58SkfnAY0Ak8Kyq/rzdcfEcvx6oB+707j4rIpHAZuCoqt7gz3saE0zWFJQzNXsAyfHRjsXwx0+O8uiKQo6dPsvQlHgeuG40N0/qO8t9TskewNs7yzh6+iyZKdYNuDd11ovpvQt5Yc+H+5PAtbinCd8kIstUdY9XsQVAnucxHXja87PN/bjXnrCBeSbkHD19loLjNfzw+jGOxfDHT47y8Bs7Odvc+mlMD7+xE6DPJImp2e4FhDaXnCTT1qnuVYGc9GUaUKSqxaraBLzGubPALgReUrf1QIqIDAEQkSzgC8CzAYzRGMes9iwONHfMIMdieHRF4afJoc3Z5lYeXVHoUETnunhIEomxUTYvkwMCmSAyAe+lSks9+/wt82vgB4CrszcRkbtFZLOIbK6oqLiggI3pTWsKyhmemsDIjH6OxXDs9Nlu7XdCZISQP2KAzezqgEAmCF999to3evssIyI3AOWquqWrN1HVpao6RVWnZGRYQ58JDmebWvmoqJK5YwbibopzxtAO6vQ72u+UqSMGUHiihur6ZqdDCSvn04sJAFXtak2IUmCY13YWcMzPMrcBN4nI9UAckCQiL3tWszMm6K0rrqSxxcXcMc52b33gutGfa4MAiI+O5IHrRjsY1bmm5qSiClsOn3S0Si7cdHYH8Qvc8zAdBM4Cz3getcAuP157E5AnIjkiEgMsApa1K7MMuEPcZgDVqlqmqg+rapaqZnvOW23JwYSSVXvLSYiJZHpuqqNx3Dwpk5/dOo7MlHgEyEyJ52e3juszDdRtJmSlEB0pbDxo1Uy9qcteTCLyE1W9wuvQWyLyflcvrKotInIfsAJ3N9fnVHW3iCz2HF8CLMfdxbUIdzfXu877SowJEu7FgcqZPSqd2KhIp8Ph5kmZfS4htBcfE8mlmck2YK6XdbmiHJAhIrmqWgwgIjmAX5X9qrocdxLw3rfE67kC93bxGmuBtf68nzHBoOB4DceqG7j/mjynQwkq07JTef6jEhqaW4mLdj6xhgN/Gqn/CVgrImtFZC2wBvhOIIMyJpS1dW+dM9qm1+iOKdmpNLW62FFa7XQoYcOfNanfFZE8oG00T4GqNgY2LGNC1+qCcsZlJjMwKc7pUILKlBGfLSA0LcfZtptw0eUdhIgkAA8A96nqdmC4pxuqMaabTtY18cnhU8xxuPdSMBrQL4a8gYk2YK4X+VPF9DzQBMz0bJcC/x6wiIwJYe/tK8elcLUliPMyNSeVLSWnaHXZhNK9wZ8EMVJV/wtoBlDVs9h6EMacl1V7y0lPjGFcZrLToQSlqdkDqGlsofB4jdOhhAV/EkSTiMTjGTQnIiMBa4MwppsamltZU1DOtWMH2fKZ56lt4j6rZuod/iSIHwHvAsNE5LfAKtxzJBljuuHjA5XUNbVy3SWDnQ4laGWmxDM0Oc5WmOsl/vRiWikiW4AZuKuW7lfVyoBHZkyIeXfXcfrHRnHZyHSnQwlaIsKM3DTe21eBqjo6j1U48KcX0ypguqq+rap/VtVKEVnaC7EZEzJaWl38Zc8Jrr54IDFRgZwjM/TNyE2jqq6J/eW1TocS8vz5S80BHhSRR7z2TQlQPMaEpI0lJzlV38z8S6166ULNHOle8XjdgSqHIwl9/iSI08DVwCAReUtErPuFMd20Ytdx4qIjuOIim5L+Qg1LTSAzJZ71xZYgAs2fBCGq2qKq9wB/AD4ErBO3MX5yuZQVu09w5UUZJMT4M/2Z6cqM3DTWF1fhsvEQAeVPgvCeXO8F4E5gZYDiMSbkbC89zfEzDVa91INmjkzjVH0zhSdsPEQgdZggRCTJ8/R3IpLa9sC9PsT3eyU6Y0LAu7uPExUhttBND5rhWUfDqpkCq7M7iFc8P7cAmz0/t3htG2O6oKqs2HWcy0alkxwf7XQ4ISNrQALDUuOtoTrAOlsw6AbPz5zeC8eY0FJ4ooaSqnruvmKk06GEnJm5aazYfQKXS21keoB0VsWU39nDnxcXkfkiUigiRSLykI/jIiKPe47vaHtdEYkTkY0isl1EdovIj8//Eo1xzru7jiMC14616qWeNiM3jeqzzew9fsbpUEJWZ10qftnJMQXmdvbCIhIJPAlci3sG2E0iskxV93gVWwDkeR7Tgac9PxuBuapaKyLRwIci8o6qru/qgozpS97ddZypI1LJ6B/rdCghx3s8xCVDrfd9IHRWxTTnAl97GlDktVTpa8BCwDtBLARe8iw9ul5EUkRkiKqWAW3DJKM9D+vPZoJKSWUdBcdr+JcbxjodSkgakhxPdloC64tP8s3Lc50OJyT51SlbRC4FxgKfLoGlqi91cVomcMRruxT33UFXZTKBMs8dyBZgFPCkqm7oILa7gbsBhg8f3uW1GNNbVuw+DsB1l1j1UqDMyE3j7Z1ltj5EgPgzF9MjwH97HnOA/wJu8uO1fbUatf8tdlhGVVtVdSKQBUzzJKlzC6suVdUpqjolI8NGqZq+493dxxmXmUzWgASnQwlZM0emUdPQwp5j1g4RCP4MlLsN91Qbx1X1LmAC4E+FaikwzGs7CzjW3TKqehpYC8z34z2N6ROOVzfwyeHTNjguwGbkutshbDxEYPiTIM6qqgto8QyeKwf8qfDbBOSJSI6IxACLgGXtyiwD7vD0ZpoBVKtqmYhkiEgKgGexomuAAv8uyRjnrdzTVr1kCSKQBiXFkZvej3WWIALCnzaIzZ4P62dwtwnUAhu7OklVW0TkPmAFEAk8p6q7RWSx5/gSYDlwPVAE1AN3eU4fArzoaYeIAF5X1T9358KMcdI7O48zamAiowYmOh1KyJsxMo1l246Rlg22PETP8mfBoHs8T5eIyLtAkqru8OfFVXU57iTgvc97bicF7vVx3g5gkj/vYUxfc7y6gQ0Hq7hvziinQwkLM3LTeGXDYeKaWkiMtckQe5K/vZjGA9lt5UVklKq+EcC4jAlaf9p2FJfCLflZTocSFtrmZTpzttkSRA/r8l9TRJ4DxgO7AZdntwKWIIxpR1V5Y+tRJg1PISe9n9PhhIWB/eMYmdGP3Q0tDHU6mBDjT7qdoao20scYP+wpO0PhiRp+crPPXtkmQGaOTGP9tmbUhkP0KH96Ma0TEUsQxvjhja1HiY4Ubhg3xOlQwsrM3HRaXUpdU4vToYQUf+4gXsSdJI7jniNJcLcvjw9oZMYEmZZWF3/adoy5YwYyoF+M0+GElele7RCm5/iTIJ4D/g7YyWdtEMaYdj4oqqSytpFbrXG616UnxpIQE0m1JYge5U+COKyq7Qe4GWPaeXPrUVISopkz2pZsd0JyfAwnzjRQ39Ria3/3EH/aIApE5BURuV1Ebm17BDwyY4JITUMzK3Yf58bxQ4mJ8ue/lelpKQnRuFT5qMhGVfcUf9JsPO62h3le+6ybqzFe3tl5nMYWF7fkZzodSthKiosmMkJYXVBuCzT1kE4ThGeqi0pVfaCX4jEmKL3xSSk56f2YNCzF6VDClggkx0eztrAcVUVs3o0L1um9sKq2An4tL2pMuCo9Vc/64pPcMinTPpQcNiAhhrLqBvaW1TgdSkjwp4ppm4gsA34H1LXttKk2jHH70zb3DPW3TLLqJaelJERTB6wuOMHYoUlOhxP0/GlNSwWqcK9BfaPncUMggzImWKgqf9hayrTsVIal2sJATouOjGB8VjKrC8qdDiUk+DOb611dlTEmXG0vraa4oo67bU3kPmPO6IE8vno/J+uaSLUBixfEnyVHs0TkTREpF5ETIvIHEbGRQMYAb24tJSYqguvH29QafcXcMQNRhff22V3EhfKniul53Cu/DQUygbc8+4wJa6qwbPsxrh07iKS4aKfDMR7jMpNJT4xldUGF06EEPX8SRIaqPq+qLZ7HC0CGPy8uIvNFpFBEikTkIR/HRUQe9xzfISL5nv3DRGSNiOwVkd0icn+3rsqYXlBV18ip+mb+ZrLdUPclERHCnNEZvFdYTkurzQ50IfxJEJUi8lURifQ8voq70bpTnjEUTwILgLHA7T5mhV0A5HkedwNPe/a3AN9T1YuBGcC9NqOs6WuOVzeQm9GPK/L8+r5ketHcMQM509DClkOnnA4lqPmTIL4OfAk4DpQBt3n2dWUaUKSqxaraBLwGLGxXZiHwkrqtB1JEZIiqlqnqVgBVrQH24q7eMqZPqGloobaxhbtm5RARYWMf+prZeelERwqrC60d4kJ0mSBU9bCq3qSqGao6UFVvVtVDfrx2JnDEa7uUcz/kuywjItm416fe4OtNRORuEdksIpsrKqzO0fSO49UNREUIX7SpNfqk/nHRTM1OZY11d70gHXZzFZF/7eQ8VdWfdPHavr5WtV/vqdMyIpII/AH4jqqe6SCQpcBSgClTpth6UibgSk/Vc7K+kSHJ8TZraB82d8xA/v3tvRw5WW9jVM5TZ3cQdT4eAN8AHvTjtUuBYV7bWcAxf8uISDTu5PBbG7Vt+pL/XXcIEAYlxTkdiunE3DHuadfXWDXTeeswQajqL9seuL+hxwN34W5L8GdU0CYgT0RyRCQGWIS7u6y3ZcAdnt5MM4BqVS0T94Q2vwH2quqvun9ZxgRGfVMLr248TGpCDLE2rXeflpuRSHZago2qvgCd/oWLSKqI/DuwA3d1VL6qPqiqXf6Lq2oLcB+wAncj8+uqultEFovIYk+x5UAxUAQ8A9zj2T8L9yp2c0Vkm+dx/XlcnzE96g9bj3KmoYXByXb3EAzmjBnIugNVnG1qdTqUoNRZG8SjwK247x7GqWptd19cVZfjTgLe+5Z4PVfgXh/nfYjv9gljHONyKc9/dJAJWclUx1nbQzC4eswgnv+ohI8PVHL1xbZGRHd1dgfxPdyjp/8ZOCYiZzyPGhHx2WBsTCh7f38FxRV1fH12jtOhGD9Ny0mlX0wkq6ya6bx01gYRoarxqtpfVZO8Hv1V1ebRNWHnuY9KGNg/lgWX2rxLwSImKoLZeemsKXAvImS6x1rZjPHD/hM1vL+vgjtmjrA1p4PMdZcMpqy6gU0lNqq6u+wv3Rg/PP9xCbFREdw+bbjToZhumn/pYPrFRPL7LUe6Lmw+xxKEMV04Xd/EG1tLuWVSJmmJsU6HY7opISaKL4wfwts7yqhvanE6nKBiCcKYLix9v5iGZhd3zbLG6WB12+Rh1DW18u6u406HElQsQRjTibLqs/zmw4PcPHEoowf3dzocc56mZg9geGoCf9ha6nQoQcUShDGd+H9/2YcqfG/eaKdDMRdARLhtchYfH6ii9FS90+EEDUsQxnRg34kafr+llDtmjrDJ3kLALZMyUYU3tx51OpSgYQnCmA785zsF9IuN4t45o5wOxfSAYakJzMxN4/dbS21MhJ8sQRjjw/riKlYVlHPPVaMY0C/G6XBMD7ltchaHqurZbCvN+cUShDHtqCo/W76XIclx3DUr2+lwTA9aMM4zJmKzNVb7wxKEMe28vbOM7aXVfPfai4iLjnQ6HNODEmKiuH7cEN7eaWMi/GEJwhgvTS0uHl1RyJjB/bk1P8vpcEwA3DY5i9rGFlbstjERXbEEYYyXVzce5lBVPQ8uGENkhM04H4qmZqe6x0Rssd5MXbEEYYxHTUMzj63az8zcNK66KMPpcEyAREQIt+Zn8tGBSo6ePut0OH1aQBOEiMwXkUIRKRKRh3wcFxF53HN8h4jkex17TkTKRWRXIGM0ps1jf93PybomHr5+DO5Vb02o+mJ+lmdMhDVWdyZgCUJEIoEngQXAWOB2ERnbrtgCIM/zuBt42uvYC8D8QMVnjLePiip59sODfGX6cMZnpTgdjgmwYakJzMhN5fdbbExEZwJ5BzENKFLVYlVtAl4DFrYrsxB4Sd3WAykiMgRAVd8HTgYwPmMAOFXXxPde305uRj/++Qvtv8OYUHXb5GGUVNXzwf5Kp0PpswKZIDIB7wnYSz37ulumUyJyt4hsFpHNFRUV5xWoCV+qyg/f3ElVXSOPL5pEfIx1aw0XN04YQmZKPL9YWWh3ER0IZILwVYnb/rfgT5lOqepSVZ2iqlMyMqxh0XTP77aU8s6u43xv3mguzUx2OhzTi2KjIvmnay9iR2m1TQPegUAmiFJgmNd2FnDsPMoYExAllXX8aNluZuSm8veX5zodjnHALZMyyRuYyKMrC2lpdTkdTp8TyASxCcgTkRwRiQEWAcvalVkG3OHpzTQDqFbVsgDGZAwAza0uvvN/24iKEH71pYk25iFMRUYI379uNMUVdbZWhA8BSxCq2gLcB6wA9gKvq+puEVksIos9xZYDxUAR8AxwT9v5IvIqsA4YLSKlIvKNQMVqws9/ry5i25HT/PTWcQxNiXc6HOOgeWMHMXFYCr/+634amludDqdPiQrki6vqctxJwHvfEq/nCtzbwbm3BzI2E742l5zkidX7uTU/kxvGD3U6HOMwEeEH80fz5Wc28PL6Q3zTqhs/ZSOpTVgprqhl8ctbyRwQz49vusTpcEwfcdnIdC7PS+fJNUXUNDQ7HU6fYQnChI3DVfV8+ZkNgPL8ndPoHxftdEimD/nBdWM4Vd/MMx8cdDqUPsMShAkLx06f5cvPrqehpZWXvzmdUQMTnQ7J9DHjspL5wrghPPtBMZW1jU6H0ydYgjAhr/xMA19+Zj3V9c3879enM2ZwktMhmT7qu/MuorHFxROri5wOpU+wBGFCWlVtI195dgPlNY288PWpjMuywXCmYyMzEvmbyVm8suEwR07WOx2O4yxBmJB1ur6Jr/5mI4dP1vObr01l8ohUp0MyQeD+a/KIihTue/WTsO/2agnChKQjJ+v5yrMbOFBey9I7pjBzZJrTIZkgMSQ5nl//7UR2lJ7me69vx+UK33maLEGYkPPW9mNc/9gHHD5Zz//83WSutMV/TDfNu2QwDy8Yw9s7y/jVX/Y5HY5jAjpQzpjeVN/Uwo+X7eH/Nh8hf3gKjy2axLDUBKfDMkHq7y/PpbiijifWFJGT3o8vTg6/NcotQZiQsOfYGf7x1a0UV9Zx75yRfOeai4iOtBtkc/5EhJ/cfCmHT9bz0Bs7yBoQz/Tc8KqqtP9BJqi1upQXPjrIzU99RE1DC7/9xnQeuG6MJQfTI6IjI3j6K5MZlprAt17eQkllndMh9Sr7X2SCUnOri99tPsI1v3qPH721h1kj03jn/su5bFS606GZEJOcEM1zX5sKwNdf2ER1ffhMxWEJwgSVxpZWXtlwmDm/WMsDv99BfHQkS76az3N3TiUtMdbp8EyIyk7vx/98dTJHTtVzy9MfsaG4yumQeoW1QZigUH22mT9+cpQl7x2grLqBCcNS+PFNlzB3zEBEbC0HE3jTc9N44a5pPPiHHfzt0vUsmjqMhxdcTHJC6M7pZQnC9Fll1Wf5654TrNxzgnUHqmhxKVNGDOA/vziey/PSLTGYXjdrVDor/+kKHvvrfp798CB/3XuCf73xEm4cPyQk/x4tQZg+o66xhb1lZ1hfXMXKPSfYUVoNQG56P755eS7XXeJe2CUU/yOa4JEQE8XD11/MTROH8vAbO/n2q5/why2lPHLjWHIzQmsSyIAmCBGZDzwGRALPqurP2x0Xz/HrgXrgTlXd6s+5Jng1tbg4caaBoopa9hw7436UnaGkqg71DFqdOCyFH8wfzbyxg23mVdMnXTI0mTfvmcVL60p4dEUhc3/5HrkZ/bgiL4PL89KZkZtGv9jg/g4esOhFJBJ4ErgWKAU2icgyVd3jVWwBkOd5TAeeBqb7ea5xiKrS2OKiqdVFY7P7Z0NzK3WNLdQ0tFDT0MyZhhZqG9zb5TUNnDjTQFm1+2dlbdPnXm9YajyXDEnmlkmZjB2SxPhhyQzsH+fQ1Rnjv8gI4a5ZOSy4dAhv7yzjg/0VvLbpMC98XEJ0pJA/fADTc9MYnBRHemIMaYmxZCTGkpYYExTJI5ARTgOKVLUYQEReAxYC3h/yC4GXPEuPrheRFBEZAmT7cW6PufG/PwyqSbk6mhlG1fcRbfekbVtVUUAVXKqffnt3qeJSpdUFrS4XrS7Fpe4xB60upanV1a14ByREMzg5nsFJsYzPSmZwUjyDk2MZkdaPi4ckkRwfuo18JjwMTo7jG7Nz+MbsHBqaW9ly6BTv76/gg32VPL5qv89zYqMiiImKIDoygsgIITpCiIwUoiIiaKtFbatMbatW7ahydUBCDK8vntmzF0VgE0QmcMRruxT3XUJXZTL9PBcAEbkbuBtg+PDh5xXoyIx+3f7Qc5p09KfSxe72f2giECGeV/N6HiHuP9ZIESIjxL0dAZEREZ/+YcdGeT+PJDE2iv5xUSTGRZEUF03/uCj6xUbZoDUTVuKiI5k1Kp1Zo9J5eIG7a/bJuiaqapuoqG2kqraJytpGTtU10dTqoqVVaXEpLa3uL2PNLv30yxvg9cWu40kDkwK0OmIgE4Svj6r2V9hRGX/Ode9UXQosBZgyZcp5Tbv460WTzuc0Y4zpUmxUJEOS4xmSHO90KN0WyARRCgzz2s4CjvlZJsaPc40xxgRQIO/9NwF5IpIjIjHAImBZuzLLgDvEbQZQraplfp5rjDEmgAJ2B6GqLSJyH7ACd1fV51R1t4gs9hxfAizH3cW1CHc317s6OzdQsRpjjDlXQPtZqepy3EnAe98Sr+cK3OvvucYYY3qPdS8xxhjjkyUIY4wxPlmCMMYY45MlCGOMMT5JR9MzBCMRqQAOnefp6UBlD4YTDOyaQ1+4XS/YNXfXCFXN8HUgpBLEhRCRzao6xek4epNdc+gLt+sFu+aeZFVMxhhjfLIEYYwxxidLEJ9Z6nQADrBrDn3hdr1g19xjrA3CGGOMT3YHYYwxxidLEMYYY3wKqwQhIvNFpFBEikTkIR/HRUQe9xzfISL5TsTZk/y45q94rnWHiHwsIhOciLMndXXNXuWmikiriNzWm/EFgj/XLCJXicg2EdktIu/1dow9zY+/7WQReUtEtnuu+S4n4uwpIvKciJSLyK4Ojvf855eqhsUD97ThB4Bc3AsSbQfGtitzPfAO7hXtZgAbnI67F675MmCA5/mCcLhmr3Krcc8YfJvTcffC7zkF95ruwz3bA52Ouxeu+YfAf3qeZwAngRinY7+Aa74CyAd2dXC8xz+/wukOYhpQpKrFqtoEvAYsbFdmIfCSuq0HUkRkSG8H2oO6vGZV/VhVT3k21+NevS+Y+fN7BvhH4A9AeW8GFyD+XPOXgTdU9TCAqgb7dftzzQr0F/dC7Im4E0RL74bZc1T1fdzX0JEe//wKpwSRCRzx2i717OtumWDS3ev5Bu5vIMGsy2sWkUzgFmAJocGf3/NFwAARWSsiW0Tkjl6LLjD8ueYngItxL1e8E7hfVV29E54jevzzK6ALBvUx4mNf+z6+/pQJJn5fj4jMwZ0gZgc0osDz55p/DTyoqq3uL5dBz59rjgImA1cD8cA6EVmvqvsCHVyA+HPN1wHbgLnASOAvIvKBqp4JcGxO6fHPr3BKEKXAMK/tLNzfLLpbJpj4dT0iMh54FligqlW9FFug+HPNU4DXPMkhHbheRFpU9Y+9EmHP8/dvu1JV64A6EXkfmAAEa4Lw55rvAn6u7gr6IhE5CIwBNvZOiL2uxz+/wqmKaROQJyI5IhIDLAKWtSuzDLjD0xtgBlCtqmW9HWgP6vKaRWQ48Abwd0H8bdJbl9esqjmqmq2q2cDvgXuCODmAf3/bfwIuF5EoEUkApgN7eznOnuTPNR/GfceEiAwCRgPFvRpl7+rxz6+wuYNQ1RYRuQ9YgbsHxHOqultEFnuOL8Hdo+V6oAiox/0NJGj5ec3/CqQBT3m+UbdoEM+E6ec1hxR/rllV94rIu8AOwAU8q6o+u0sGAz9/zz8BXhCRnbirXx5U1aCdBlxEXgWuAtJFpBR4BIiGwH1+2VQbxhhjfAqnKiZjjDHdYAnCGGOMT5YgjDHG+GQJwhhjjE+WIIwxxvhkCcKYXuCZDwgR+ZH3tjF9mSUIY3rHRBF5HEgVkZuB/3A4HmO6FDYD5YzpLSKSDbwLbAAm4Z7O4g7gKWAdEK2q/+BYgMb4yQbKGdPDPAniIDBbVT8SkefwrMWAe0TvKmCKqv6zc1Ea0zWrYjImMI6o6kee5y8Dl6vqt4Eqz7xP/+JYZMb4yRKEMYHR/tbcBaCqP/L8tFt30+dZgjAmMIaLyEzP89uBD50MxpjzYQnCmMDYC3xNRHYAqcDTDsdjTLdZLyZjAsOlqoudDsKYC2F3EMYYY3yybq7GGGN8sjsIY4wxPlmCMMYY45MlCGOMMT5ZgjDGGOOTJQhjjDE+/X8W/dS4CddI4QAAAABJRU5ErkJggg==\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "plt.plot(x, β_distribution, label = 'Beta')\n", | |
| "plt.xlabel('p*')\n", | |
| "plt.ylabel('Normalized Likelihood')\n", | |
| "plt.scatter(β.expect(), .04)\n", | |
| "[plt.axvline(i, color = 'green') for i in β.interval(.9)]\n", | |
| "plt.legend();" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "id": "f41f6a60", | |
| "metadata": {}, | |
| "source": [ | |
| "## Conclusion:\n", | |
| "\n", | |
| "After flipping a coin 20 times and obtaining 13 heads and 7 tails, we may conclude that the true value is $p=0.636$ with a credible range of 0.464 to 0.794" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "id": "95472284", | |
| "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.9.7" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 5 | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment