diff --git a/.all-contributorsrc b/.all-contributorsrc index 8d6f60fa8..e05333be5 100644 --- a/.all-contributorsrc +++ b/.all-contributorsrc @@ -184,6 +184,16 @@ "contributions": [ "doc" ] + }, + { + "login": "meraldoantonio", + "name": "Meraldo Antonio", + "avatar_url": "https://avatars.githubusercontent.com/u/37468543?v=4", + "profile": "https://github.com/meraldoantonio", + "contributions": [ + "code", + "example" + ] } ] } diff --git a/examples/04_bayesian_conjugate.ipynb b/examples/04_bayesian_conjugate.ipynb new file mode 100644 index 000000000..3e21a9574 --- /dev/null +++ b/examples/04_bayesian_conjugate.ipynb @@ -0,0 +1,527 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Bayesian Conjugate" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "\n", + "from skpro.distributions import Beta\n", + "from skpro.regression.bayesian_proportion import BayesianProportionEstimator" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Conjugate Bayesian analysis is a statistical method where the posterior distribution belongs to the same family as the prior distribution. \n", + "\n", + "In this notebook, we will use the `skpro` estimator `BayesianProportionEstimator` to demonstrate the application of Bayesian conjugates in the estimation of proportion." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Based on the data above, we aim to estimate the proportion $p$ of heads. However, due to the small sample size, we might be uncertain about our estimate. This is where Bayesian analysis, and specifically the use of conjugate priors, becomes valuable." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Motivation for Estimating Proportions: The Coin Toss Example\n", + "\n", + "Estimating proportions is a fundamental problem in statistics, often arising in situations where we want to understand the likelihood of a specific outcome. A classic example is the coin toss.\n", + "\n", + "Imagine we have a coin, and we want to determine if it is fair. A fair coin has an equal probability (0.5) of landing heads or tails. To estimate this probability, we can perform a series of coin tosses and observe the outcomes.\n", + "\n", + "Suppose we toss the coin 10 times and observe the following results:\n", + "\n", + "- Heads: 8 times\n", + "- Tails: 2 times\n", + "\n", + "This coin toss sample is represented in the pandas Series `y` shown below." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0 1\n", + "1 1\n", + "2 1\n", + "3 1\n", + "4 1\n", + "5 1\n", + "6 1\n", + "7 0\n", + "8 0\n", + "9 1\n", + "Name: head, dtype: int64" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Define the number of ones and zeros\n", + "ones = [1] * 8\n", + "zeros = [0] * 2\n", + "\n", + "# Combine them into a single list\n", + "data = ones + zeros\n", + "\n", + "# Shuffle the list to randomize the order\n", + "np.random.shuffle(data)\n", + "\n", + "# Create a Series from the shuffled list\n", + "y = pd.Series(data, name=\"head\")\n", + "\n", + "y" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Based on this data, we aim to estimate the proportion `p` of heads. However, due to the small sample size, we might be uncertain about our estimate. This is where Bayesian analysis becomes valuable." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "# Prior\n", + "\n", + "A Bayesian analysis begins with a **prior belief** about the parameter of interest, which in this case is the proportion `p`. This prior belief is represented by a **prior _distribution_**, which encapsulates our initial assumptions about `p` before observing any data. The choice of distribution can be based on previous knowledge, expert opinion, or even non-informative (if we have no prior knowledge).\n", + "\n", + "A **conjugate prior** is a specific type of prior distribution that, when combined with the likelihood of the observed data, yields a posterior distribution that belongs to the same family as the prior distribution. This property simplifies the process of updating our beliefs with new data, as the mathematical form of the distribution remains consistent.\n", + "\n", + "The **Beta distribution** is often chosen as the prior distribution because it is a conjugate prior to the Binomial likelihood, which describes the probability of observing a given number of heads in a series of coin tosses.\n", + "\n", + "\n", + "The distribution, denoted as $\\text{Beta}(\\alpha, \\beta)$, is parameterized by two positive parameters, $\\alpha$ and $\\beta$. These parameters reflect our prior beliefs about the proportion $p$. For example:\n", + "\n", + "- If $\\alpha = 1$ and $\\beta = 1$, we have a uniform prior, which indicates that all values of $p$ between 0 and 1 are equally likely.\n", + "- If $\\alpha > 1$ and $\\beta > 1$, the prior might be more peaked around a particular value, indicating a stronger prior belief about $p$.\n", + "\n", + "Suppose in our case that we have some conviction that the coin is fair, meaning we believe the proportion $p$ is centered around 0.5. To reflect this belief, we will choose a Beta distribution with parameters $\\alpha = 2$ and $\\beta = 2$. This distribution is symmetric and centered around 0.5, indicating our prior belief that the coin is likely fair but allowing for some uncertainty.\n", + "\n", + "To implement this in our analysis, we will instantiate a new `BayesianProportionEstimator` object, `B`. This object encapsulates our prior belief on `p` and provides the framework for updating this belief based on observed data. \n", + "\n", + "Here's how we can set this up in code:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "B = BayesianProportionEstimator(prior=Beta(2, 2))\n", + "# alternative instantiation:\n", + "B = BayesianProportionEstimator(prior_alpha=2, prior_beta=2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can use the plot function to visualize the probability density function (PDF) of our prior distribution. This visualization helps us understand our initial belief about `p`. By plotting the PDF of the $Beta(2, 2)$ distribution, we can see how our belief is symmetrically distributed around the value 0.5, reflecting our prior conviction that the coin is fair with some allowance for variability." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABTNElEQVR4nO3deVxU5eIG8GcWdplRVtkE3BFUEMIts0wpNcyy8mpllt5f3ltpeaubWdp2s7pl3hatTLtZaqapLVpKmYq7IKiIKyA7sskMO8zM+f0xQnHdAIF35szz/Xzmj6Yz8nBczsN53/O+CkmSJBARERHJhFJ0ACIiIqK2xHJDREREssJyQ0RERLLCckNERESywnJDREREssJyQ0RERLLCckNERESyohYdoKOZTCbk5eXB1dUVCoVCdBwiIiJqBkmSUF5eDl9fXyiV1743Y3PlJi8vDwEBAaJjEBERUStkZ2fD39//msfYXLlxdXUFYD45Go1GcBoiIiJqDr1ej4CAgMbr+LXYXLlpGIrSaDQsN0RERFamOVNKOKGYiIiIZIXlhoiIiGSF5YaIiIhkheWGiIiIZIXlhoiIiGSF5YaIiIhkheWGiIiIZIXlhoiIiGSF5YaIiIhkheWGiIiIZIXlhoiIiGSF5YaIiIhkxeY2ziQieTGaJNQajKitN6HOaILRJMFOpYS9Sgk7tQIOahVUyutvtEdE8sFyQ0QWq6rOgIziSmSVVCGztAqZJVW4oK9BSWUdSipqUVpZh6o643V/Ha2THbo426GLiz3cXRwQ4OaEbm7O6ObmjEB3ZwS5u0Ct4o1sIrlguSEii1BTb0RSVhmO5ZThRJ4eJ/J0SC+uhCQ1/9dQKRVQKRSoM5qavK+rroeuuh7nS6qu+Dl7tRJ9vF0R4uOKfj4aRAW5IcRHwzs+RFaK5YaIhKg1GHE44yL2pxfjYHopjuaUod54eZPp4myHbu4uCHRzRpC7M3w7O8G9kwPcXOzh0ckeGkc7ONiZh6Ea7r5IkgSjSUK9UUJlnQFlVXUoraxHaWUdisprkH2xGtmlVcgqrcL54kpU1hlxPFeH47m6xq/byUGNQYFdEB3UBbf09kSYrxZKlh0iq6CQpJb8XGT99Ho9tFotdDodNBqN6DhENqW4ohY7ThVix8lCxJ8tQuX/DCl5axwQGdgFob5a9PPVINRXAy9Xx3bNZDJJyL5YhZP5eqTml+N4ThkSMi+ivMbQ5DhPVwfc2tsTo/p6YWQfTzjb82dDoo7Ukus3yw0RtSt9TT22pRTgh6N52HuuGKY//Yvj6eqAEb08MCTYHYO7u6GbmzMUCvF3R4wmCacLynEoowT700uw52xxkyLmZKfC6H7emDDQF7f09oCDWiUwLZFtYLm5BpYbovZnMknYm1aMbw5lI+7kBdQZ/pgDE+qrwe0h3hgd4mU1Qz21BiMSzl/EjlOFiEu9gKzSP+buaBzVmBDuiynR3RDqqxWYkkjeWG6ugeWGqP2UVNRifWIO1hzMalIAeni6YGK4HyaE+yLQ3UVgwhsnSRKO5ejww9E8/HQsDxf0tY3/b2BAZ0yNDkDsQF8OWxG1MZaba2C5IWp7aUUVWL47HRuP5DY+qeTqoMa9g/xwf1QAQn01FjHc1NZMJgn700uw5lAWtp8oaJwQ3dnZDtOGBOLhoUHwdHUQnJJIHlhuroHlhqjtHMm6iE93pWF76oXGR7YH+Gvx0OBA3DXQx6buXhRX1OK7xBys/tNdK3u1EpMG+eH/bumBYA/rvmNFJBrLzTWw3BDduMTMi3hv+2nsSytpfG90iDdmjeyOqCA3gcnEM5okbD9RgE93pyM5uwwAoFQA9w7yx+xRvdDN3VlsQCIrxXJzDSw3RK2XmqfHe9tP47dThQAAO5UCE8P98PjI7ujp5So4nWWRJAmJmRexbGda4/lSKxW4P8ofT9zWE/5dWHKIWoLl5hpYbohaLru0Cm//cgo/HcsHYF4J+L5B/pg9uhf8OjsJTmf5krPL8H7cGew6UwQAsFcp8ejNQXjytp5wdbQTnI7IOrDcXAPLDVHzVdUZsPT3NHwWn974OHfsQF88M7oXunt2EpzO+iScL8XiuDONw3nuLvb4R0wfTL4pgFs9EF0Hy801sNwQXZ8kSfjhaB4WbT2FAn0NAGBYD3e8NL4f+vny782NkCQJO08X4fUtqUgvqgQA9O3qioWxoRjaw11wOiLLxXJzDSw3RNeWVlSBeRuP41BGKQDAv4sTXhrfD3eEesvycW5R6o0mfH0gE0t+PQtddT0A4P5If7w4LgRdXOwFpyOyPCw318ByQ3Rl9UYTPt2Vhg9+O4c6owlOdir8/dYe+Ost3eFox+0F2svFyjq8u/00Vh/MAgC4udjj5btCMDHcj2WS6E9Ybq6B5Ybockezy/DP747hVEE5AOCW3p7418QwBLjxiZ6OkphZinkbj+PMhQoAwIheHnhr0gBO2Ca6hOXmGlhuiP5QZzBhya9n8MmuNJgkoIuzHRbE9uNdA0HqDCYsj0/Hf347izqDCa4OaiycEIpJg/j7QcRycw0sN0Rm5wrL8fS6ZKTk6gEAd4f7YsFd/eDeidsFiJZRXIl/fJuMI1llAIA7Qr3x5j39+XtDNo3l5hpYbsjWSZKErw5k4l9bTqLWYEJnZzu8eU9/jOvvIzoa/YnBaMKnu9Ox5NczqDdK8Ohkj7cnDcDtId6ioxEJwXJzDSw3ZMsuVtZh7rfJ+P20eTG5Eb088O79A+GtcRScjK7mRJ4Oz6xLbpyL89cRwXj+zr6wUykFJyPqWC25fgv927F7927ExsbC19cXCoUCmzdvbvZn9+7dC7VajfDw8HbLRyQnR7IuYvwH8fj9dBHs1UosjO2HLx+NZrGxcKG+Wvzw5M14dHgQAGB5fAYe+HQ/ci5WiQ1GZMGElpvKykoMHDgQH330UYs+p9PpMG3aNNx+++3tlIxIPiRJwoo9GXjgk/3I09Ug2MMFm/8+HI8OD4aSq+JaBUc7FRbGhuLThyPh6qhGUlYZxn+wB7+mXhAdjcgiWcywlEKhwKZNmzBx4sTrHvuXv/wFvXr1gkqlwubNm5GcnNzsr8NhKbIl5TX1eH7DMfycUgAAGN/fB29N6s/9jKxYdmkVnlxzBEdzdACAJ27rgX+M6cOiSrJnNcNSrfHFF18gLS0NCxcubNbxtbW10Ov1TV5EtuB8cSXuWboPP6cUwE6lwKsTQvHR1AgWGysX4OaM9bOGNQ5Tffx7GmauSoC+pl5sMCILYlXl5uzZs3jhhRewevVqqNXqZn1m0aJF0Gq1ja+AgIB2Tkkk3t5zxbj74704V1gBb40Dvn18KB4ZFsS1UmTCPGcqFO9PHggHtRI7ThVi4sd7kVZUIToakUWwmnJjNBoxdepUvPrqq+jdu3ezPzdv3jzodLrGV3Z2djumJBJLkiR8ue88pq08BF11PcIDOuPHJ29GRLcuoqNRO7gnwh8bZg2Dj9YR6UWVmPjRXuw4xXk4RFYz56asrAxdunSBSvXHHjcmkwmSJEGlUmH79u0YNWrUdb8O59yQXBmMJiz44QTWXNqj6N4IP7x5b3/uC2UDispr8ffViTh8/iKUCuDVCaF4eGiQ6FhEbaol1+/mje1YAI1Gg+PHjzd5b+nSpdixYwc2bNiA4OBgQcmIxKusNeCJNUew83QRlApg3tgQzBwRzGEoG+Hp6oDVM4fgpc3H8W1CDl7+/gSySqswb2wIJxqTTRJabioqKnDu3LnG/87IyEBycjLc3NzQrVs3zJs3D7m5uVi1ahWUSiXCwsKafN7LywuOjo6XvU9kSwrLa/DYfw8jJVcPRzslPpwyCGP6cRVbW2OvVuLtSQMQ6O6Cf287jeXxGci5WI33J4fz7h3ZHKFzbhISEhAREYGIiAgAwNy5cxEREYEFCxYAAPLz85GVlSUyIpFFSyuqwL1L9yElVw83F3us/esQFhsbplAo8MRtPfGfv4TDXqXEzykFmLL8AEor60RHI+pQFjPnpqNwzg3JRcL5UsxclYCyqnoEujvjy0ejEeThIjoWWYiD6SX4v68SoauuR0+vTvhqRjR8tE6iYxG1mqzXuSEiYOfpQjy04iDKquoxMKAzNv5tGIsNNTG4uzu++9tQ+Ggdca6wAvct24+M4krRsYg6BMsNkZX5JSUff12VgJp6E27t44m1fx0M904OomORBerp5Yr1s4Yi2MMFuWXVuP+T/UjN40KmJH8sN0RWZOORHDyxJgn1Rgnj+/vgs4ej4GxvNQ89kgD+XZyxftZQ9PPRoLiiFpM/24+E86WiYxG1K5YbIivx9YFMzP32KIwmCfdH+uODKRGwV/OvMF2fRycHrP2/IbgpqAvKawx4aMVB7D1XLDoWUbvhv4xEVmD57nS8tDkFADB9WBDenjQAKq5fQi2gdbLDqscG49Y+nqipN+Gx/x7GnrMsOCRPLDdEFu7z+HT8a+tJAOYdoBfG9uPCbNQqTvYqfPpwJEaHeKHWYMKMLw8j/myR6FhEbY7lhsiCrdyTgTe2mIvN06N74bk7+nLVYbohDmoVPn5wUGPBmfllAgsOyQ7LDZGFWrX/PF77KRUAMHtUTzw9uvkbxhJdi4NahaUPRmJ0iHdjwdl9hgWH5IPlhsgCfX0gEwu+PwEA+PutPfDMGBYbalv2aiWWPmjeqqPWYMLMVQnYx0nGJBMsN0QW5ptDWY2Thx8f2R3P3dGHQ1HULuzVSnw81Vxw6i4VnCNZF0XHIrphLDdEFuSnY3mYt+k4AGDmzcF44U7OsaH2Za9W4qOpERjRywNVdUZMX3kIJ/O50B9ZN5YbIgux60wRnlmXDEkCHhzcDfPHh7DYUIdwUJufoooM7AJ9jQEPrzjErRrIqrHcEFmAxMxSzPoqEfVGCXcN8MFrd4ex2FCHcrZXY+X0mxpXMn7o84PILasWHYuoVVhuiAQ7VaDHo18cRnW9ESN7e2LxA+FcoI+E0DrZYdWMaHS/tBfVw58fRElFrehYRC3GckMkUFZJFR5ecQj6GgMGdeuMZQ8N4pYKJJRHJwd8PXMw/Do7Ib24Eo99mYCqOoPoWEQtwn9FiQQprazDtJUHUVReiz7erlg5/SZugkkWwbezE1bNiEZnZzsczS7DU2uSYDCaRMciajaWGyIBauqNmPnlYZwvqYJf44XEXnQsokY9PDthxSNRcFAr8dupQrz8/QlIkiQ6FlGzsNwQdTCTScIz65JxJKsMGkc1vnzsJnhrHEXHIrpMZKAb/vOXCCgUwNpDWfhoxznRkYiaheWGqIO9ufUkfk4pgL1Kic+mRaGnl6voSERXdWdYV7w6IRQA8F7cGaxPyBaciOj6WG6IOtB/92bg8z0ZAIB/3z8AQ7q7C05EdH3ThgZh1sgeAIB5G49jz1lu00CWjeWGqINsP1GAVy9thPncHX1wd7if4EREzff8HX1wd7gvDCYJf1+diLSiCtGRiK6K5YaoA6Tm6THnG/Pqw1Oiu+Hvt/YQHYmoRZRKBd6eNACDunWGvsaAmV8moKyqTnQsoitiuSFqZ8UVtfjrqgRU1xtxc08PvH53KFcfJqvkaKfCpw9Hwa+zEzKKK/H31UdQz0fEyQKx3BC1o1qDEX/7OhG5ZdUI9nDBx1MHQa3iXzuyXp6uDvj8kSi42KuwL60EC/iIOFkg/itL1E4kScLLm1Nw+PxFuDqqsXxaFLTOdqJjEd2wEB9Nk0fEv9h7XnQkoiZYbojaycq95/FtQg6UCuDDKRHo6dVJdCSiNjO6nzdeHBsCAHhjSyrizxYJTkT0B5Ybonaw60wR/rXF/GTUi+NCcGsfL8GJiNrezBHBuD/SHyYJeGptErJLq0RHIgLAckPU5rJKqvDUmiMwScADUf6YcXOw6EhE7UKhUOD1iWEY6K9FWVU9Hv8qEdV1RtGxiFhuiNpSTb0Rs75OhL7GgPCAznh9YhifjCJZc7RTYdlDkXB3sUdqvh4vbjrOCcYkHMsNURuRJAnzN6UgNV8Pdxd7LHtoEBzUKtGxiNqdb2cnfDR1EFRKBTYl5eK/+86LjkQ2juWGqI2sOZSF7478MYHYR+skOhJRhxnawx3zxvYFALyx5SQOpJcITkS2jOWGqA0kZV3EKz+cAAA8f2dfDOvpITgRUcebcXMw7g73hdEk4ck1R1CgqxEdiWwUyw3RDSqpqL20UquEO0K98fgt3UVHIhJCoVDgrXsHoG9XVxRX1GH22iQYuIIxCcByQ3QDjCYJs79JQr6uBt09XPDu/QM5gZhsmpO9eYKxi70Kh86XYnHcGdGRyAax3BDdgA93nMXecyVwtlfhk4cj4erIFYiJgj1c8NakAQCApTvTsPN0oeBEZGtYbohaaX9aCT747SwA4M17+qO3t6vgRESWI3agLx4a0g0A8My6ZOTrqgUnIlvCckPUCiUVtZjzTRJMEnB/pD8mRviJjkRkcV4a3w+hvhpcrKrHU2uSuIM4dRiWG6IWMpkk/GP9URSW16KnVye8eneo6EhEFsnRToWlDw6Cq4MaCZkX8e7206IjkY0QWm52796N2NhY+Pr6QqFQYPPmzdc8fuPGjRgzZgw8PT2h0WgwdOhQbNu2rWPCEl2yPD4dO08XwUGtxMdTB8HZXi06EpHFCnR3wTv3mefffLorHTtOXRCciGyB0HJTWVmJgQMH4qOPPmrW8bt378aYMWOwdetWJCYm4rbbbkNsbCySkpLaOSmR2ZGsi/j3NvNPnwtjQ9GnK+fZEF3P2P4+mD4sCADw7PpjKCzn+jfUvhSShWwColAosGnTJkycOLFFnwsNDcXkyZOxYMGCK/7/2tpa1NbWNv63Xq9HQEAAdDodNBrNjUQmG6Orqse4D+KRW1aNuwb44MMpEXzsm6iZag1GTPx4H07m63FLb0/8d/pNUCr594eaT6/XQ6vVNuv6bdVzbkwmE8rLy+Hm5nbVYxYtWgStVtv4CggI6MCEJBeSJOHFzceRW1aNbm7OWHRvfxYbohZwUKvwwV/C4aBWYveZInzB/aeoHVl1uXnvvfdQWVmJBx544KrHzJs3DzqdrvGVnZ3dgQlJLjYl5WLLsXyolQp8OCWC69kQtUIvb1e8dFc/AMDbP59Cap5ecCKSK6stN2vXrsUrr7yCdevWwcvL66rHOTg4QKPRNHkRtUR2aRUWfm/eN+rp0b0wMKCz2EBEVuyhwd0wOsQbdUYT5nyThJp6o+hIJENWWW7WrVuHGTNm4Ntvv8Xo0aNFxyEZM5ok/OPboyivNSAqsAv+dmtP0ZGIrJpCocDbk/rD09UBZwsr8K8tJ0VHIhmyunKzdu1aTJ8+HWvWrMH48eNFxyGZ+3R3Gg6dL0UnBzXenxwOFSdAEt0w904OeO/+gQCArw5k4reTfDyc2pbQclNRUYHk5GQkJycDADIyMpCcnIysrCwA5vky06ZNazx+7dq1mDZtGt577z0MGTIEBQUFKCgogE6nExGfZC4lV4fF282b/i2M7YcAN2fBiYjk45benphxczAA4PkNx1BcUXudTxA1n9Byk5CQgIiICERERAAA5s6di4iIiMbHuvPz8xuLDgB8+umnMBgMeOKJJ+Dj49P4mjNnjpD8JF/VdUbM+SYJBpOEsWFdcV+kv+hIRLLz/J190LerK0oq6/DSphRYyMokJAMWs85NR2nJc/JkuxZ+n4Iv92fCy9UB256+BV1c7EVHIpKl1Dw97v54D+qNEpZMDuc+bXRVNrPODVF72HuuGF/uzwQAvHv/QBYbonbUz1eDObf3AgAs+D4FBTquXkw3juWG6E/Ka+rx/IZjAICHhwTilt6eghMRyd+skT0w0F8LfY0B//zuGIen6Iax3BD9yZtbTyG3rBoBbk54YWxf0XGIbIJapcR7D5hXL951pgjfHOZiq3RjWG6ILok/W4S1h8wT2N+ZNBAuDtztm6ij9PTqhOfu6AMAeOOnVGSXVglORNaM5YYI5uGof14ajnpkaCCG9nAXnIjI9jw2PBjRwW6orDPi2fVHYTJxeIpah+WGCMC/tpxEnq4G3dyc8U8ORxEJoVQq8O59A+Fsr8LBjFJ8uf+86EhkpVhuyOb9eYz/3/cNgLM9h6OIROnm7ox540IAAP/edprDU9QqLDdk0/Q19XjhO/Nw1PRhQRjcncNRRKI9GN0N0cFuqKoz4sVNx/n0FLUYyw3ZtLd+PoV8XQ0C3Z3x/J19RMchIpiHp966tz8c1ErEny3Gd0dyRUciK8NyQzbrUEYp1hw0Px319iQORxFZku6enfDMmN4AgNd/SkVhORf3o+ZjuSGbVFNvxLyN5uGoKdEBGMLhKCKLM/PmYIT5aaCrrscrP5wQHYesCMsN2aSlv59DWlElPF0d8MLYENFxiOgK1Col3p40ACqlAluPF+CXlHzRkchKsNyQzTlzoRzLdqUBAF6dEAqtk53gRER0NaG+Wswa2R0A8PL3J6CrqheciKwByw3ZFJNJwgvfHUO9UcLoEG+MDesqOhIRXcdTo3qhu6cLispr8caWVNFxyAqw3JBN+fpgJo5klaGTgxqvTwyFQqEQHYmIrsPRToV3Jg0AAKxPzMHB9BLBicjSsdyQzcjXVeOdX04DAJ6/sw98tE6CExFRc0UFuWFKdDcAwPzNKagzmAQnIkvGckM2QZIkvLw5BRW1Bgzq1hkPDQ4UHYmIWuiFO/vC3cUe5worsDw+XXQcsmAsN2QTtp0owK8nC2GnUuCtSQOgVHI4isjaaJ3t8NJd5qcbP/jtLDJLKgUnIkvFckOyV1lrwGs/michPn5LD/T2dhWciIhaa2K4H4b3dEetwYQF35/g1gx0RSw3JHsf7DiLPF0NAtyc8OSonqLjENENUCgUeP3uMNirlNh1pghbjnPtG7ocyw3J2tkL5VgRnwEAeCU2FI52KsGJiOhGdffshL/d2gMA8NqPqdDXcO0baorlhmRLkiS8tDkFBpN5TZvbQ7xFRyKiNvK3W3sg2MMFheW1eG/badFxyMKw3JBsfZ+ch4MZpXC0U2JhbD/RcYioDTnaqfDGxDAAwKoDmTieoxOciCwJyw3Jkq66Hm9sOQnAvLppgJuz4ERE1NaG9/TA3eG+kCRgwQ8pMJk4uZjMWG5Ilt6PO4Piilp093DBzBHBouMQUTt5cVwIXOxVSMoqw8akXNFxyEKw3JDspOTqsGr/eQDAa3eHwUHNScREcuWtccTs23sBAN76+SQnFxMAlhuSGZPJPInYJAF3DfDBzb08REcionb26PBgdPd0QXFFHf7z61nRccgCsNyQrHx3JAfJ2WVwsVfhpfGcRExkC+zVSrwSGwoA+O++8zhzoVxwIhKN5YZko7ymHm9f2hhz9u290FXrKDgREXWUW3p74o5QbxhNEhZy5WKbx3JDsvHR7+dQXFGLIHdnTB8eJDoOEXWwl8b3g4Naif3pJdh6vEB0HBKI5YZkIaO4Eiv3mFcifvmufpxETGSDAtycG1cufmNLKqrqDIITkSgsNyQL/9qSinqjhFt6e2JUXy/RcYhIkFkje8C/ixPydTX4+PdzouOQICw3ZPV2nSnCrycLoVYqsOCuECgUCtGRiEgQRzsVXr7L/DDB8t0ZyCypFJyIRGC5IatWbzTh9Z9SAQDThgahp5er4EREJFpMP2+M6OWBOqMJb/18SnQcEoDlhqzaV/szca6wAm4u9pgzupfoOERkARQKBV4a3w9KBfBzSgEOppeIjkQdjOWGrFZJRS3e//UMAODZmD7QOtkJTkRElqJPV1f8JbobAOCNLSe575SNYbkhq7U47gzKawwI8dFg8k0BouMQkYWZO6Y3XB3UOJ6r475TNkZoudm9ezdiY2Ph6+sLhUKBzZs3X/czu3btQmRkJBwdHdG9e3d88skn7R+ULM7pgnKsPZQFAHglth9USk4iJqKmPDo54IlRPQEA/952io+G2xCh5aayshIDBw7ERx991KzjMzIyMG7cOIwYMQJJSUl48cUXMXv2bHz33XftnJQszZtbT8IkAXeGdsXg7u6i4xCRhXp0eBAC3JxwQV+LT3ali45DHUQt8ouPHTsWY8eObfbxn3zyCbp164YlS5YAAEJCQpCQkIB3330XkyZNaqeUZGl2nynCrjNFUCsVeGFsX9FxiMiCOahVmDc2BH9ffQSf7U7DlOgA+GidRMeidmZVc27279+PmJiYJu/dcccdSEhIQH39lbe5r62thV6vb/Ii62U0SXhz60kAwMNDAxHk4SI4ERFZurFhXREd5IaaehPeubT/HMmbVZWbgoICeHt7N3nP29sbBoMBxcXFV/zMokWLoNVqG18BAZx4as2+O5KDUwXl0DiqMXsUH/0moutTKBR46a4QAMCmpFwkZ5eJDUTtzqrKDYDLVp9t2Pn1aqvSzps3DzqdrvGVnZ3d7hmpfVTVGfDedvNPXU+N6oUuLvaCExGRtRjg3xn3RvgBAN74KZW7hsucVZWbrl27oqCg6U6vhYWFUKvVcHe/8qRSBwcHaDSaJi+yTp/HZ+CCvhb+XZwwbVig6DhEZGWev7MvHO2USMi8iO2pF0THoXZkVeVm6NChiIuLa/Le9u3bERUVBTs7LuAmZ4XlNfhkVxoA4J939uWu30TUYl21jphxczAA4J1fTsFgNAlORO1FaLmpqKhAcnIykpOTAZgf9U5OTkZWlnn9knnz5mHatGmNx8+aNQuZmZmYO3cuTp48iZUrV2LFihV49tlnRcSnDvR+3FlU1RkRHtAZdw3wER2HiKzU4yN7oIuzHdKKKvFtQo7oONROhJabhIQEREREICIiAgAwd+5cREREYMGCBQCA/Pz8xqIDAMHBwdi6dSt27tyJ8PBwvP766/jggw/4GLjMnblQjnWHzX8O5o/nrt9E1HoaRzs8delhhCW/nuHCfjKlkGxsVpVer4dWq4VOp+P8Gyvx6BeH8PvpItwZ2hWfPBwpOg4RWblagxGjF+9Cdmk1no3pjSf55KVVaMn126rm3JDtOZBegt9Pmxfs+ycX7COiNuCgVuHZmD4AgE92paOkolZwImprLDdksSRJwju/nAIA/CU6AMFcsI+I2kjsAF+E+WlQUWvAhzvOiY5DbYzlhizWrycLcSSrDI52Si7YR0RtSqlU4IU7zQv7rT6YiaySKsGJqC2x3JBFMpok/Hub+a7NY8OD4aVxFJyIiOTm5l4eGNHLA/VGCe9u57YMcsJyQxZpU1IuzlyogNbJDo+P7CE6DhHJ1Atj+0KhAH44mofjOTrRcaiNsNyQxak1GPF+3BkAwN9v7QGtExdoJKL2EeqrxcRw87YMb/1yUnAaaissN2RxVh/IQm5ZNbw1DnhkWJDoOEQkc3PH9IadSoG950qwL+3KmzCTdWG5IYtSUWvAR7+bn1x4enRvONpxmwUial8Bbs6YGt0NAPDuttPcVFMGWG7IoizfnY7Syjp093DB/ZH+ouMQkY14YlRPONopcSSrDDtOFYqOQzeI5YYsRnFFLT6PTwcA/COmD9Qq/vEkoo7h5erYOAz+7vYzMJl498aa8epBFuPj38+hss6I/n5ajA3rKjoOEdmYWbf0gKuDGifz9diaki86Dt0AlhuyCPm6aqw+YN4c8/k7+0Cp5OaYRNSxurjYY+aI7gCAxXFnYDCaBCei1mK5IYvw8e/nUGc0ITrYDTf39BAdh4hs1GM3B6GLsx3SiyqxMSlXdBxqJZYbEi7nYhXWHc4GYH4kU6HgXRsiEsPV0Q5/u9W8cOh/fj2LWoNRcCJqDZYbEu6jHedQb5QwvKc7hnR3Fx2HiGzctKFB8HJ1QG5ZdeMPXmRdWG5IqMySSqxPzAFgvmtDRCSao50KT43qCQD4cMc5VNfx7o21Ybkhof7z21kYTRJG9vZEZKCb6DhERACAyTd1g38XJxSV1+LL/edFx6EWYrkhYdKKKrD50oQ93rUhIktir1bi6dHmf5c+3ZWGylqD4ETUEiw3JMx/fj0LkwSMDvHCwIDOouMQETUxMdwXwR4uuFhVj1X7M0XHoRZguSEhTheU48djeQCAZ3jXhogskFqlbJx789lu3r2xJiw3JMR/fjsDSQLGhnVFqK9WdBwioiuaMJB3b6wRyw11uBN5Omw9XgCFAo1j2kREloh3b6wTyw11uCW/ngUA3DXAF326ugpOQ0R0bbx7Y31YbqhDpeTqEJd6AUoFMOf2XqLjEBFdF+/eWB91Sz8gSRJ27dqF+Ph4nD9/HlVVVfD09ERERARGjx6NgICA9shJMvHhDvNdm9iBvujp1UlwGiKi5pkw0Bcf7jiHjOJKrNqf2bhFA1mmZt+5qa6uxptvvomAgACMHTsWW7ZsQVlZGVQqFc6dO4eFCxciODgY48aNw4EDB9ozM1mpUwV6bDtxAQoF8ORtPUXHISJqNt69sS7NLje9e/fGkSNH8Mknn0Cv1+PAgQP47rvv8PXXX2Pr1q3IyspCWloaRowYgcmTJ2P58uXtmZus0Ec7zgEAxoX5oJc359oQkXXh3BvroZAkSWrOgSkpKQgLC2vWL1pXV4fMzEz06mV5cyr0ej20Wi10Oh00Go3oODbjXGEFxry/C5IE/DxnBEJ8eO6JyPpsPJKDud8eRRdnO+z55yi4OLR4dge1Ukuu382+c9PcYgMA9vb2FllsSJyPfz8HSQLG9PNmsSEiq/Xnuzfcc8pyteppqZdffhlG4+W7pOp0OkyZMuWGQ5G8nC+uxPfJ5j2kZo9i6SUi66VWKRvnDK6Iz+CO4RaqVeVm1apVGD58ONLS0hrf27lzJ/r374/z58+3VTaSiaU7z8EkAbf18UR/f65GTETWbUK4L/y7OKGksg7fHM4SHYeuoFXl5tixYwgKCkJ4eDiWL1+O5557DjExMZg+fTr27NnT1hnJimWXVmHjEfNdm6e4rg0RyYCdSolZI82Pgn+2Ox21Bt69sTStmgml1WrxzTffYP78+Xj88cehVqvx888/4/bbb2/rfGTllu1Kg8EkYUQvDwzq1kV0HCKiNnFfpD8++O0s8nU12HgkF1Oiu4mORH/S6hWKP/zwQ7z//vuYMmUKunfvjtmzZ+Po0aNtmY2sXL6uGhsScgAAT3GuDRHJiKOdCv93S3cAwLKdaTAYTYIT0Z+1qtyMHTsWr776KlatWoXVq1cjKSkJt9xyC4YMGYJ33nmnrTOSlfp0VzrqjCYMDnZDdLCb6DhERG1q6uBucHOxR1ZpFX46li86Dv1Jq8qNwWDAsWPHcN999wEAnJycsGzZMmzYsAHvv/9+mwYk61RYXoO1h8wT7WZzrg0RyZCzvRqPDQ8CYF7uwmRq1rJx1AFaVW7i4uLg6+t72fvjx4/H8ePHbzgUWb8V8RmoNZgwqFtnDOvhLjoOEVG7eHhoEFwd1DhbWIHtqRdEx6FL2nxXcA8PDwDmDTabY+nSpQgODoajoyMiIyMRHx9/zeNXr16NgQMHwtnZGT4+Pnj00UdRUlJyw7mp7eiq6vH1AfPS5E+O6gmFQiE4ERFR+9A62WHasEAADYuV8u6NJWh2uQkJCcGaNWtQV1d3zePOnj2Lv/3tb3j77bev+2uuW7cOTz/9NObPn4+kpCSMGDECY8eORVbWldcN2LNnD6ZNm4YZM2bgxIkTWL9+PQ4fPoyZM2c299ugDvD1wUxU1hnRt6srbuvjJToOEVG7emx4MJzsVDieq8Pus8Wi4xBasLfUjh078M9//hPnzp1DTEwMoqKi4OvrC0dHR1y8eBGpqanYs2cPUlNT8eSTT+LFF1+87t4PgwcPxqBBg7Bs2bLG90JCQjBx4kQsWrTosuPfffddLFu2rMnigR9++CHeeecdZGdnN+sb5t5S7au6zoib396Bkso6/Ocv4bg73E90JCKidvfaj6lYuTcD0UFu+HbWUNFxZKld9pYaNWoUDh8+jC1btqBr165Ys2YNnnzySTz44IN45ZVXcPbsWUybNg05OTl46623rvuF6+rqkJiYiJiYmCbvx8TEYN++fVf8zLBhw5CTk4OtW7dCkiRcuHABGzZswPjx46/6dWpra6HX65u8qP2sT8xGSWUd/Ls4YXx/H9FxiIg6xP/d0h32KiUOnS/FwXROlRCtxYv4DRs2DMOGDbvhL1xcXAyj0Qhvb+8m73t7e6OgoOCqX3v16tWYPHkyampqYDAYMGHCBHz44YdX/TqLFi3Cq6++esN56frqjSZ8uisdAPD4Ld2hVrX5lC4iIovUVeuISZH+WHsoCx/9fg6Du/NBCpGaffVxc3NDcbF5LPGxxx5DeXl5mwT438mmkiRddQJqamoqZs+ejQULFiAxMRG//PILMjIyMGvWrKv++vPmzYNOp2t8NXf4ilrup2N5yC2rhkcne9wfFSA6DhFRh/rbyB5QKoD4s8VIydWJjmPTml1u6urqGod0vvzyS9TU1NzQF/bw8IBKpbrsLk1hYeFld3MaLFq0CMOHD8dzzz2HAQMG4I477sDSpUuxcuVK5OdfeQElBwcHaDSaJi9qeyaThGU7zXOhHh0eDEc7leBEREQdq5u7M8YPMC+T8tnudMFpbFuzh6WGDh2KiRMnIjIyEpIkYfbs2XBycrrisStXrrzur2dvb4/IyEjExcXhnnvuaXw/Li4Od9999xU/U1VVBbW6aWSVynwR5eN3Yu04VYgzFyrQyUGNh4YEio5DRCTE47d0x49H87DleD6eu6MPAtycRUeySc2+c/P1119j3LhxqKiogEKhgE6nw8WLF6/4aq65c+fi888/x8qVK3Hy5Ek888wzyMrKahxmmjdvHqZNm9Z4fGxsLDZu3Ihly5YhPT0de/fuxezZsxEdHX3FRQWpY0iShKU7zwEAHhoSCK2TneBERERihPlpMaKXB4wmCZ/H8+6NKM2+c+Pt7Y233noLABAcHIyvvvoK7u43NmFq8uTJKCkpwWuvvYb8/HyEhYVh69atCAw0/+Sfn5/fZM2b6dOno7y8HB999BH+8Y9/oHPnzhg1alSz1tSh9nP4/EUcySqDvVrZuBQ5EZGtevyWHog/W4x1CdmYM7o33FzsRUeyOc1e50YuuM5N25v+xSHsPF2EqYO74c17+ouOQ0QklCRJiP1oD1Jy9Xh6dC88Pbq36Eiy0JLrd7Pv3HzwwQfNDjB79uxmH0vWLTVPj52ni6BUmMeaiYhsnUKhwOO39MBTa5Pw5b7z+L9busPZvsUrr9ANaPbZ/t/dvouKilBVVYXOnTsDAMrKyuDs7AwvLy+WGxvyyS7zE1LjB/gi0N1FcBoiIsswNqwrAtyckF1ajfUJOXhkWJDoSDal2ROKMzIyGl//+te/EB4ejpMnT6K0tBSlpaU4efIkBg0ahNdff70985IFyblYhS3HzY/g864NEdEf1Col/m+E+d/F5fHpMBhNghPZllYtIfvyyy/jww8/RJ8+fRrf69OnD95//3289NJLbRaOLNsXe8/DaJIwvKc7wvy0ouMQEVmU+yID4OZij5yL1diacuWV96l9tKrc5Ofno76+/rL3jUYjLly4cMOhyPLpquvxzSHzk2x/HcG7NkRE/8vJXoXpl4ajPtmZxvXYOlCrys3tt9+Ov/71r0hISGj8zUpISMDjjz+O0aNHt2lAskzfHMpCZZ0RfbxdMbK3p+g4REQW6eEhgXCyUyE1X48954pFx7EZrSo3K1euhJ+fH6Kjo+Ho6AgHBwdER0fDx8cHn3/+eVtnJAtTZzDhi73nAQAzRwRfdS8wIiJb18XFHpNvMu+117CxMLW/Vj2b5unpia1bt+Ls2bM4efIkDAYDwsLC0Ls3n+W3BT8dy0OBvgaerg6YEM6VoYmIrmXmiGB8dSATe86ZN9TkHMX216o7NwCwYsUK3HPPPbj//vsxZcoU3HvvvbxrYwMkSWrcEG76sCA4qLlBJhHRtfh3ccZdA3wAACv2ZAhOYxta/bTUnDlzEBsbi/Xr12P9+vWIjY3FM888w6elZG7vuRKcKiiHs70KDw7uJjoOEZFVmHFzMADgx6N5KNDVCE4jf60allq2bBmWL1+OKVOmNL43YcIEDBgwAE899RTeeOONNgtIluWzSxvBPRAVgM7O3C+FiKg5Bvh3RnSwGw5llOLL/efxzzv7io4ka626c2M0GhEVFXXZ+5GRkTAYDDcciizTyXw9dp8xb7XQ8FMIERE1z8xL/26uPpCJylpeK9tTq8rNQw89hGXLll32/meffYYHH3zwhkORZfo83jxWPDbMBwFuzoLTEBFZl9tDvBHk7gx9jQHfHckRHUfWWr2T14oVK7B9+3YMGTIEAHDgwAFkZ2dj2rRpmDt3buNxixcvvvGUJNwFfQ1+OJoLwDzzn4iIWkalVOCxm4Ox4PsTWLknAw8NDoRSyaU02kOryk1KSgoGDRoEAEhLM2+c6OnpCU9PT6SkpDQex/VP5OO/+86j3ighOsgNEd26iI5DRGSV7ov0x3vbz+B8SRV+O1WIMf28RUeSpVaVm99//72tc5AFq6w1YPWBTADAX7lBJhFRqznbqzF1cDcs25mGz+PTWW7aSavXuSHbsfFIDvQ1BgS5O+P2vl6i4xARWbVHhgZBrVTgYEYpjufoRMeRJZYbuiaTSWrcauHR4cEcHyYiukFdtY6IHWhe3X3FHm7J0B5Ybuiadp0tQnpxJVwd1JgU6S86DhGRLDQsp/HTsXzk66oFp5Eflhu6ppWXlgqffFMAOjm0+uE6IiL6kzA/LYZ0d4PBJOHLfZmi48gOyw1d1dkL5Yg/WwylAnhkWJDoOEREsjLzZvMDGmsOclG/tsZyQ1f1xb7zAIAx/by5aB8RURsb1dcLwR4u0NcYsCGRi/q1JZYbuqKyqjpsvLSC5mPDuWgfEVFbUyoVeGx4EADgy33nYTJJYgPJCMsNXdHaQ9moqTehn48G0cFuouMQEcnSvYP84eqoRnpxJXafLRIdRzZYbugy9UYTVu0/DwB47OZgrjRNRNROXBzUeCAqAIB5JXhqGyw3dJltJwqQr6uBRyd7xA70ER2HiEjWpg0NhEIB7DxdhPSiCtFxZIHlhi7T8Pj3g4MD4aBWCU5DRCRvge4ujau/r9rPx8LbAssNNZGcXYYjWWWwVynx4JBuouMQEdmE6cPMD25sSMxBeU294DTWj+WGmvhir/muzV0DfeDl6ig4DRGRbRje0x09vTqhotaA7/hY+A1juaFGF/Q12HIsHwAf/yYi6kgKhaJxsdQv92fysfAbxHJDjVYfyITBJCE6yA1hflrRcYiIbMq9EX5wdVQjo7gSu/hY+A1huSEAQJ3BhDWHsgFwqwUiIhFcHNSY3PBY+N7zYsNYOZYbAgD8nJKP4opaeGscEBPqLToOEZFNmjY0CAoFsOtMEdL4WHirsdwQgD8eP5waHQg7Ff9YEBGJ0M3dGbf3Nf+A+RUfC281XsUIKbk6JGZehFqpwJToANFxiIhs2vRLUwPWJ2TzsfBWYrmhxp8Oxvb3gZeGj38TEYnU8Fh4ZZ2Ru4W3EsuNjSurqsP3R3MBmJcAJyIisRQKRePdG+4W3josNzZufUIOaupNCPHRICqwi+g4REQE4N5BfnB1UON8SRX2phWLjmN1hJebpUuXIjg4GI6OjoiMjER8fPw1j6+trcX8+fMRGBgIBwcH9OjRAytXruygtPJiMkn46oB5SMq8cRt3/yYisgTO9mpMivQHwInFrSG03Kxbtw5PP/005s+fj6SkJIwYMQJjx45FVlbWVT/zwAMP4LfffsOKFStw+vRprF27Fn379u3A1PKx60wRskqroHFU4+5wX9FxiIjoTx66tL/frycvIK+sWnAa66IW+cUXL16MGTNmYObMmQCAJUuWYNu2bVi2bBkWLVp02fG//PILdu3ahfT0dLi5uQEAgoKCrvk1amtrUVtb2/jfer2+7b4BK7dq/3kAwP1RAXC2F/pHgYiI/kdPL1cM7e6O/eklWHsoC/+I6SM6ktUQduemrq4OiYmJiImJafJ+TEwM9u3bd8XP/PDDD4iKisI777wDPz8/9O7dG88++yyqq6/eaBctWgStVtv4Cgjgo84AkFlSiZ1nzMt7PzyEE4mJiCzRw5ce9Fh7KBt1BpPgNNZDWLkpLi6G0WiEt3fT1XC9vb1RUFBwxc+kp6djz549SElJwaZNm7BkyRJs2LABTzzxxFW/zrx586DT6Rpf2dnZbfp9WKuvD2RCkoCRvT0R5OEiOg4REV3BmH7e8HJ1QHFFLbaduPK1kS4nfELx/05ilSTpqhNbTSYTFAoFVq9ejejoaIwbNw6LFy/Gf//736vevXFwcIBGo2nysnXVdUasO9ywjxTv2hARWSo7lRJTos1zbxoeAKHrE1ZuPDw8oFKpLrtLU1hYeNndnAY+Pj7w8/ODVvvHjtUhISGQJAk5OVzoqLm+T86FvsaAADcnjOztJToOERFdw5ToblApFTiUUYrTBeWi41gFYeXG3t4ekZGRiIuLa/J+XFwchg0bdsXPDB8+HHl5eaio+GMzsTNnzkCpVMLf379d88qFJEmN+0g9PCQQKiUf/yYismRdtY6I6Wf+of9r3r1pFqHDUnPnzsXnn3+OlStX4uTJk3jmmWeQlZWFWbNmATDPl5k2bVrj8VOnToW7uzseffRRpKamYvfu3Xjuuefw2GOPwcnJSdS3YVWOZJUhNV8PB7USD0RxcjURkTVoePBj45EcVNQaBKexfEKf/508eTJKSkrw2muvIT8/H2FhYdi6dSsCA82/ifn5+U3WvOnUqRPi4uLw1FNPISoqCu7u7njggQfwxhtviPoWrM6ag+bzGTvQF52d7QWnISKi5hjawx3dPV2QXlSJTUm5fMr1OhSSJNnUphV6vR5arRY6nc7mJheXVdVh8Ju/odZgwsa/D8OgbtxugYjIWnyxNwOv/piKPt6u+OXpETa3qnxLrt/Cn5aijvPdkVzUGsz7SEUEdBYdh4iIWuDeQf5wslPh9IVyHD5/UXQci8ZyYyMkScKag+aJaFMHd7O5xk9EZO20TnaYGGHeKoePhV8by42NOJhRirSiSjjbqzCR+0gREVmlhy7NtfklJR9F5bXXOdp2sdzYiIaJxHeH+8HV0U5wGiIiao1QXy0GdeuMeqOEdYevvsm0rWO5sQHFFbX4OSUfAPDg4G6C0xAR0Y1ouHuz9lA2TCabeiao2VhubMCGxBzUGyUM9NcizE97/Q8QEZHFGtffB1onO+SWVSP+XLHoOBaJ5UbmTCYJaw+Zb10+OJjrIhARWTtHOxXuHeQHAI0PilBTLDcytzetGJklVXB1UOOugT6i4xARURto2Ezz15OFKNTXCE5jeVhuZK5hIvG9g/zgbC90QWoiImojvb1dERXYBUaThPWJ3Dj6f7HcyNgFfQ22p14AAEzlkBQRkaw03L1ZeyiLE4v/B8uNjH17OBtGk4SowC7o09VVdBwiImpD4wf4QOOoRs7FauzhxOImWG5kymiS8M3hbADAg0P4+DcRkdyYJxb7A0DjgyNkxnIjU7vOFCK3rBqdne0wNowTiYmI5KhhaCou9QIKyzmxuAHLjUw1TCS+b5A/HO1UgtMQEVF76NPVFZGBXWAwSdjAicWNWG5kqEBXgx2nCgEAU7giMRGRrDXcvfmGKxY3YrmRoQ2J2TBJQHSwG3p4dhIdh4iI2tH4/j5wdVQjq7QKe9M4sRhguZEdk0nCugTzROK/3BQgOA0REbU3J3sV7o0wr1jMicVmLDcysz+9BNml1XB1VHMiMRGRjWiYgrD9xAUUldcKTiMey43MNDz+PTHcD072nEhMRGQL+nbVIKJbZ04svoTlRkYuVtZhW0oBAGAyh6SIiGzK1IaJxYe5YjHLjYxsSspFndGEMD8Nwvy0ouMQEVEHumuAL1wd1cgsqcK+tBLRcYRiuZEJSZKw7tKQ1OSb+Pg3EZGtcbJXYWK4eWJxw4MltorlRiaSsstw+kI5HO2UmDDQV3QcIiISoGFKwrYTBSirqhOcRhyWG5lYd8jc0sf194HWyU5wGiIiEiHMT4t+PhrUGUz4PjlPdBxhWG5koKLWgB+Pmf8Q/4VDUkRENu2BKPNmmg1TFWwRy40M/HQ0D1V1RnT3cMFNQV1ExyEiIoEmRvjBXqVEar4eKbk60XGEYLmRgW8aJxIHQKFQCE5DREQidXa2R0yoNwDgWxudWMxyY+VOFeiRnF0GtVKBewf5i45DREQWoGFi8eakXNTUGwWn6XgsN1auYUx1dIg3PF0dBKchIiJLMLyHB/w6O0FfY8C2EwWi43Q4lhsrVlNvxKakXADA5GiuSExERGZKpQL3RZrv5tvi0BTLjRXbnnoBZVX18NE64pZenqLjEBGRBbk/yh8KBbD3XAmyS6tEx+lQLDdW7NtLQ1L3R/pDpeREYiIi+oN/F2cM7+EBAFhvY3dvWG6sVG5ZNfamFQMA7o/ikBQREV3ugUsTizck5sBoQ5tpstxYqY2JOZAkYEh3NwS4OYuOQ0REFiimnze0TnbI09Vgz7li0XE6DMuNFZIkCRuO5AAA7o/kXRsiIroyRzsVJoab9xv81oZWLGa5sUKHz19EZkkVXOxVGNu/q+g4RERkwRqGpranFqC00jY202S5sUINE8PGD/CBs71acBoiIrJkob5ahPlpUG+UsPnS8iFyJ7zcLF26FMHBwXB0dERkZCTi4+Ob9bm9e/dCrVYjPDy8fQNamMpaA7YczwfAicRERNQ8D1y6XnybkA1Jkv/EYqHlZt26dXj66acxf/58JCUlYcSIERg7diyysrKu+TmdTodp06bh9ttv76CkluPnlAJU1RkR5O6MqEBukklERNd390A/2KuVOFVQjuM2sJmm0HKzePFizJgxAzNnzkRISAiWLFmCgIAALFu27Jqfe/zxxzF16lQMHTr0ul+jtrYWer2+ycuaNQxJ3Rfpz00yiYioWbTOdrgj1DxH87vEHMFp2p+wclNXV4fExETExMQ0eT8mJgb79u276ue++OILpKWlYeHChc36OosWLYJWq218BQRY71BOVkkVDmaUQqEAN8kkIqIWmTTIDwDww9E81BlMgtO0L2Hlpri4GEajEd7e3k3e9/b2RkHBlTf5Onv2LF544QWsXr0aanXzJtLOmzcPOp2u8ZWdbb2PwjU8/n1zTw/4dnYSnIaIiKzJiF6e8HJ1wMWqeuw4VSg6TrsSPqH4f4dWJEm64nCL0WjE1KlT8eqrr6J3797N/vUdHByg0WiavKyRySQ13kps2AyNiIiouVRKBe65dPdmg8yHpoSVGw8PD6hUqsvu0hQWFl52NwcAysvLkZCQgCeffBJqtRpqtRqvvfYajh49CrVajR07dnRUdCEOpJcgt6waro7qxnFTIiKilrjv0pSGnacLUVJRKzhN+xFWbuzt7REZGYm4uLgm78fFxWHYsGGXHa/RaHD8+HEkJyc3vmbNmoU+ffogOTkZgwcP7qjoQqy/1LJjB/rC0U4lOA0REVmjXt6uGOCvhcEk4fvkPNFx2o3QFeDmzp2Lhx9+GFFRURg6dCg+++wzZGVlYdasWQDM82Vyc3OxatUqKJVKhIWFNfm8l5cXHB0dL3tfbspr6vFzyqW1bTgkRUREN+C+SH8cy9FhQ2IOHrs5WHScdiG03EyePBklJSV47bXXkJ+fj7CwMGzduhWBgYEAgPz8/OuueWMLthzLR029CT08XRAe0Fl0HCIismKxA3zx+k+pSM3X42S+HiE+1jkX9VoUki0sVfgner0eWq0WOp3OaiYXT1q2D4mZF/HC2L6YNbKH6DhERGTlZn2ViF9OFGDmzcF46a5+ouM0S0uu38KflqJrSy+qQGLmRSgVwD0RfqLjEBGRDEy6NMVhc3Iu6o3yW/OG5cbCbTxi3uRsZG9PeGscBachIiI5uLWPJ9xd7FFcUYfdZ4pEx2lzLDcWzGSSsOnSDq5ckZiIiNqKnUqJu8PNowHfHZHfmjcsNxbs0PlS89o2DmqM6Xf52j9EREStNSnSXG5+TS1EWVWd4DRti+XGgm281KbH9ffh2jZERNSmQn21CPHRoM5owo9H5bXmDcuNhaqpN+Ln4+bVm+8dxInERETU9ho209xwaX6nXLDcWKi41AsorzXAr7MTbgpyEx2HiIhk6O5wP6iUChzNLsO5wnLRcdoMy42FahiSuifCD0rl5RuJEhER3ShPVwfc1scTALAhUT53b1huLFBReS12ny0GgMYdXImIiNrDpEtP425KyoHRJI91fVluLNCPR/NgNEkYGNAZPTw7iY5DREQyNirEC1onO1zQ1+JAeonoOG2C5cYCbUwyD0ndyxWJiYionTmoVRg/wAcAGtdWs3YsNxbmzIVypOTqoVYqEDvQV3QcIiKyAQ3b+/ySUoDqOqPgNDeO5cbCNGy3cGsfL7i52AtOQ0REtiCyWxf4d3FCRa0Bv568IDrODWO5sSAmk4Tvk83lZhInEhMRUQdRKhWYeGk7hs0yGJpiubEgB9JLkK+rgcZRjVEhXqLjEBGRDZkYYZ4KsetMEUoqagWnuTEsNxbku0tDUuMH+MJBze0WiIio4/T0ckV/Py0MJgk/HcsXHeeGsNxYiOo6I35JMf9h4pAUERGJ0DCx2NqfmmK5sRDbUwtQWWdENzdnRAZ2ER2HiIhsUOxAX6iUCiRnlyGjuFJ0nFZjubEQDUNS90T4QaHgdgtERNTxPF0dcHNPDwDWPbGY5cYCFOprsOdsEYA/bgkSERGJ0HAd2pycC0myzu0YWG4swPfJeTBJwKBunRHk4SI6DhER2bCYUG8426uQWVKFI1llouO0CsuNBdh8aW2bey5tXkZERCSKs70ad4R2BWC9Q1MsN4KdKyzHiTzzdgvj+/uIjkNERNQ4NPXTsTzUGUyC07Qcy41gPyTnAQBG9vbkdgtERGQRhvVwh6erAy5W1WP3mSLRcVqM5UYgSZLw/VFzuZkQzk0yiYjIMqhVSky4tHnzpmTrG5piuREoObsMmSVVcLJTYUw/b9FxiIiIGjUMTf2aegH6mnrBaVqG5Uag7y8NSZlnpqsFpyEiIvpDqK8GPb06odZgwi/HC0THaRGWG0EMRlPj3h0NO7ESERFZCoVCYbXbMbDcCLI/vQTFFbXo4myHm3t5iI5DRER0mbsvzQc9kFGCfF214DTNx3IjSMOQ1PgBPrBT8beBiIgsj38XZ0QHuUGSgJ+OWs9O4byqClBTb8QvKebxy7s5JEVERBYs9tLdmx8uPd1rDVhuBNhxqhAVtQb4dXZCZDfuAE5ERJZrXFhXqJQKHM/VIb2oQnScZmG5EeD7S2sGTAj3hVLJHcCJiMhyuXf6Y6dwa7l7w3LTwXTV9fj9lHm1x7u5cB8REVmBhgX9fjiaZxU7hbPcdLBtKQWoM5rQx9sVfbtqRMchIiK6rphQbziolUgvqsSJPL3oONfFctPBNv9pSIqIiMgauDraYVRfLwDAj1YwNMVy04Eu6GuwP70EwB+3+IiIiKxBw3Xrx6N5MJkse2hKeLlZunQpgoOD4ejoiMjISMTHx1/12I0bN2LMmDHw9PSERqPB0KFDsW3btg5Me2N+PJoHSQKiArsgwM1ZdBwiIqJmu62vFzo5qJGnq0Fi1kXRca5JaLlZt24dnn76acyfPx9JSUkYMWIExo4di6ysrCsev3v3bowZMwZbt25FYmIibrvtNsTGxiIpKamDk7dOwyxzTiQmIiJr42inwh2hXQEAPyRb9tCUQhI47Xnw4MEYNGgQli1b1vheSEgIJk6ciEWLFjXr1wgNDcXkyZOxYMGCZh2v1+uh1Wqh0+mg0XTchN70ogqMem8XVEoFDr14O9w7OXTY1yYiImoLu84U4ZGVh+DuYo+DL94OdQeusN+S67ewOzd1dXVITExETExMk/djYmKwb9++Zv0aJpMJ5eXlcHNzu+oxtbW10Ov1TV4iNGy3cEsvDxYbIiKySsN7uMPdxR4llXXYm1YiOs5VCSs3xcXFMBqN8Pb2bvK+t7c3Cgqat7X6e++9h8rKSjzwwANXPWbRokXQarWNr4CAgBvK3RqSJDXOLud2C0REZK3UKiXG9fcBYNlDU8InFCsUTVfolSTpsveuZO3atXjllVewbt06eHl5XfW4efPmQafTNb6ys7NvOHNLncjTI724Eo52Sozu5339DxAREVmohqVMtp8oQE29UXCaKxNWbjw8PKBSqS67S1NYWHjZ3Zz/tW7dOsyYMQPffvstRo8efc1jHRwcoNFomrw62k/HzDupjro005yIiMhaRXbrAl+tI8prDdh5ulB0nCsSVm7s7e0RGRmJuLi4Ju/HxcVh2LBhV/3c2rVrMX36dKxZswbjx49v75g3TJIk/HTMfOvurgF8SoqIiKybUqlA7EDL3ilc6LDU3Llz8fnnn2PlypU4efIknnnmGWRlZWHWrFkAzENK06ZNazx+7dq1mDZtGt577z0MGTIEBQUFKCgogE6nE/UtXFdydhlyLlbD2V6F2/pcffiMiIjIWjSUm19PFqK8pl5wmssJLTeTJ0/GkiVL8NprryE8PBy7d+/G1q1bERgYCADIz89vsubNp59+CoPBgCeeeAI+Pj6Nrzlz5oj6Fq6rYUhqTD9vONmrBKchIiK6caG+GnT3dEGdwYTtJy6IjnMZoevciNCR69yYTBKGvbUDBfoaLJ8WhTGcTExERDLxn1/P4v1fz2Bkb098+Vh0u389q1jnxhYkZF5Egb4Gro5q3NLbQ3QcIiKiNtPw1NSec8UoqagVnKYplpt21DCROKZfVzioOSRFRETyEezhgv5+WhhNEn5Oad76dB2F5aadGIwmbD1unm8TO9BHcBoiIqK213B9a/hh3lKw3LSTgxmlKK6oQxdnOwzvySEpIiKSn4bVig9mlKJQXyM4zR9YbtpJQ4u9M6wr7DpwYzEiIqKO4t/FGRHdOkOS0DhaYQl41W0H9UZT4/gjF+4jIiI5a7jObWG5kbc954pRVlUPj04OGNLdXXQcIiKidjP+0tDU4fMXka+rFpzGjOWmHfx01Nxex/XvCpXy+puAEhERWauuWkfcFNQFALDlmGXcvWG5aWO1BiO2n+CQFBER2Y6G691PLDfytOt0EcprDeiqcURUYBfRcYiIiNrd2P5doVSY91PMLq0SHYflpq01tNbxA3yg5JAUERHZAC9XRwwONs8xtYSnplhu2lB1nRG/njRvIHbXAC7cR0REtuOuxgX9WG5kZcepQlTVGeHfxQnhAZ1FxyEiIuowd4aaH6I5nqvD+eJKoVlYbtpQw8J9dw3whULBISkiIrId7p0cMKyHeWhK9Jo3LDdtpKLWgB2nCgFwSIqIiGxTw/Xvx6Ni95piuWkjWSVV8NI4INjDBaG+GtFxiIiIOtwdoV2hViqgUipQXlMvLIdCkiRJ2FcXQK/XQ6vVQqfTQaNp2xIiSRKKymvhpXFs01+XiIjIWhTqa9rlOtiS6zfv3LQhhULBYkNERDbNEq6DLDdEREQkKyw3REREJCssN0RERCQrLDdEREQkKyw3REREJCssN0RERCQrLDdEREQkKyw3REREJCssN0RERCQrLDdEREQkKyw3REREJCssN0RERCQrLDdEREQkK2rRATqaJEkAzFunExERkXVouG43XMevxebKTXl5OQAgICBAcBIiIiJqqfLycmi12mseo5CaU4FkxGQyIS8vD66urlAoFM3+nF6vR0BAALKzs6HRaNoxIf0Zz7sYPO8dj+dcDJ73jtfacy5JEsrLy+Hr6wul8tqzamzuzo1SqYS/v3+rP6/RaPgXQACedzF43jsez7kYPO8drzXn/Hp3bBpwQjERERHJCssNERERyQrLTTM5ODhg4cKFcHBwEB3FpvC8i8Hz3vF4zsXgee94HXHObW5CMREREckb79wQERGRrLDcEBERkayw3BAREZGssNwQERGRrLDc/MnSpUsRHBwMR0dHREZGIj4+/prH79q1C5GRkXB0dET37t3xySefdFBSeWnJed+4cSPGjBkDT09PaDQaDB06FNu2bevAtPLQ0j/rDfbu3Qu1Wo3w8PD2DShTLT3vtbW1mD9/PgIDA+Hg4IAePXpg5cqVHZRWPlp63levXo2BAwfC2dkZPj4+ePTRR1FSUtJBaa3f7t27ERsbC19fXygUCmzevPm6n2nz66lEkiRJ0jfffCPZ2dlJy5cvl1JTU6U5c+ZILi4uUmZm5hWPT09Pl5ydnaU5c+ZIqamp0vLlyyU7Oztpw4YNHZzcurX0vM+ZM0d6++23pUOHDklnzpyR5s2bJ9nZ2UlHjhzp4OTWq6XnvEFZWZnUvXt3KSYmRho4cGDHhJWR1pz3CRMmSIMHD5bi4uKkjIwM6eDBg9LevXs7MLX1a+l5j4+Pl5RKpfSf//xHSk9Pl+Lj46XQ0FBp4sSJHZzcem3dulWaP3++9N1330kApE2bNl3z+Pa4nrLcXBIdHS3NmjWryXt9+/aVXnjhhSse//zzz0t9+/Zt8t7jjz8uDRkypN0yylFLz/uV9OvXT3r11VfbOppstfacT548WXrppZekhQsXsty0QkvP+88//yxptVqppKSkI+LJVkvP+7///W+pe/fuTd774IMPJH9//3bLKGfNKTftcT3lsBSAuro6JCYmIiYmpsn7MTEx2Ldv3xU/s3///suOv+OOO5CQkID6+vp2yyonrTnv/8tkMqG8vBxubm7tEVF2WnvOv/jiC6SlpWHhwoXtHVGWWnPef/jhB0RFReGdd96Bn58fevfujWeffRbV1dUdEVkWWnPehw0bhpycHGzduhWSJOHChQvYsGEDxo8f3xGRbVJ7XE9tbuPMKykuLobRaIS3t3eT9729vVFQUHDFzxQUFFzxeIPBgOLiYvj4+LRbXrlozXn/X++99x4qKyvxwAMPtEdE2WnNOT979ixeeOEFxMfHQ63mPxmt0Zrznp6ejj179sDR0RGbNm1CcXEx/v73v6O0tJTzbpqpNed92LBhWL16NSZPnoyamhoYDAZMmDABH374YUdEtkntcT3lnZs/USgUTf5bkqTL3rve8Vd6n66tpee9wdq1a/HKK69g3bp18PLyaq94stTcc240GjF16lS8+uqr6N27d0fFk62W/Fk3mUxQKBRYvXo1oqOjMW7cOCxevBj//e9/efemhVpy3lNTUzF79mwsWLAAiYmJ+OWXX5CRkYFZs2Z1RFSb1dbXU/4YBsDDwwMqleqyJl9YWHhZm2zQtWvXKx6vVqvh7u7eblnlpDXnvcG6deswY8YMrF+/HqNHj27PmLLS0nNeXl6OhIQEJCUl4cknnwRgvuhKkgS1Wo3t27dj1KhRHZLdmrXmz7qPjw/8/Pyg1Wob3wsJCYEkScjJyUGvXr3aNbMctOa8L1q0CMOHD8dzzz0HABgwYABcXFwwYsQIvPHGG7wr3w7a43rKOzcA7O3tERkZibi4uCbvx8XFYdiwYVf8zNChQy87fvv27YiKioKdnV27ZZWT1px3wHzHZvr06VizZg3HwVuopedco9Hg+PHjSE5ObnzNmjULffr0QXJyMgYPHtxR0a1aa/6sDx8+HHl5eaioqGh878yZM1AqlfD392/XvHLRmvNeVVUFpbLppVGlUgH4424Cta12uZ62eiqyzDQ8LrhixQopNTVVevrppyUXFxfp/PnzkiRJ0gsvvCA9/PDDjcc3PLr2zDPPSKmpqdKKFSv4KHgrtPS8r1mzRlKr1dLHH38s5efnN77KyspEfQtWp6Xn/H/xaanWael5Ly8vl/z9/aX77rtPOnHihLRr1y6pV69e0syZM0V9C1appef9iy++kNRqtbR06VIpLS1N2rNnjxQVFSVFR0eL+hasTnl5uZSUlCQlJSVJAKTFixdLSUlJjY/fd8T1lOXmTz7++GMpMDBQsre3lwYNGiTt2rWr8f898sgj0siRI5scv3PnTikiIkKyt7eXgoKCpGXLlnVwYnloyXkfOXKkBOCy1yOPPNLxwa1YS/+s/xnLTeu19LyfPHlSGj16tOTk5CT5+/tLc+fOlaqqqjo4tfVr6Xn/4IMPpH79+klOTk6Sj4+P9OCDD0o5OTkdnNp6/f7779f8d7ojrqcKSeJ9NiIiIpIPzrkhIiIiWWG5ISIiIllhuSEiIiJZYbkhIiIiWWG5ISIiIllhuSEiIiJZYbkhIiIiWWG5ISIiIllhuSEiIiJZYbkhIiIiWWG5ISIiIllhuSEiq1dUVISuXbvizTffbHzv4MGDsLe3x/bt2wUmIyIRuHEmEcnC1q1bMXHiROzbtw99+/ZFREQExo8fjyVLloiORkQdjOWGiGTjiSeewK+//oqbbroJR48exeHDh+Ho6Cg6FhF1MJYbIpKN6upqhIWFITs7GwkJCRgwYIDoSEQkAOfcEJFspKenIy8vDyaTCZmZmaLjEJEgvHNDRLJQV1eH6OhohIeHo2/fvli8eDGOHz8Ob29v0dGIqIOx3BCRLDz33HPYsGEDjh49ik6dOuG2226Dq6srfvrpJ9HRiKiDcViKiKzezp07sWTJEnz11VfQaDRQKpX46quvsGfPHixbtkx0PCLqYLxzQ0RERLLCOzdEREQkKyw3REREJCssN0RERCQrLDdEREQkKyw3REREJCssN0RERCQrLDdEREQkKyw3REREJCssN0RERCQrLDdEREQkKyw3REREJCv/D3SwoB79y52MAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "B.prior.plot(\"pdf\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Likelihood and Posterior\n", + "\n", + "The likelihood function represents the probability of observing the data given the parameter $p$. For coin tosses, this is described by the Binomial distribution. If we observe $k$ heads out of $n$ tosses, the likelihood is given by:\n", + "\n", + "$$\\text{Binomial}(k | n, p) = \\binom{n}{k} p^k (1 - p)^{n - k}$$\n", + "\n", + "Using Bayes' theorem, we combine the prior distribution with the likelihood of the observed data to obtain the posterior distribution, which represents our updated beliefs about $p$ after observing the data. \n", + "\n", + "\n", + "Because the Beta distribution is a conjugate prior to the Binomial likelihood, the posterior distribution is also a Beta distribution with updated parameters. That is, if the prior is $\\text{Beta}(\\alpha, \\beta)$ and we observe $k$ heads in $n$ tosses, the posterior distribution is:\n", + "\n", + "$$\\text{Beta}(\\alpha + k, \\beta + n - k)$$\n", + "\n", + "This means that we simply add the number of observed heads to $\\alpha$ and the number of observed tails to $\\beta $ to get the parameters of the posterior distribution." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Using the `BayesianProportionEstimator`, the process of updating our prior belief with observed data is straightforward and can be accomplished in a single line of code. We simply use the fit method on our data.\n", + "\n", + "However, to maintain compatibility with the estimator's API, we need to supply an additional argument, `X`, which is expected to be a DataFrame or Series of the same length as `y`. This X parameter will be ignored during the fitting process, so for simplicity, we can simply set `X` to be the same as `y`." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
BayesianProportionEstimator(prior=Beta(alpha=2, beta=2), prior_alpha=2,\n",
+       "                            prior_beta=2)
Please rerun this cell to show the HTML repr or trust the notebook.
" + ], + "text/plain": [ + "BayesianProportionEstimator(prior=Beta(alpha=2, beta=2), prior_alpha=2,\n", + " prior_beta=2)" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "X = y.copy()\n", + "B.fit(X, y)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As expected, after fitting the BayesianProportionEstimator to our observed data, we obtain a posterior distribution represented by the Beta distribution with updated parameters $\\alpha=10$ and $\\beta=4$:" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
Beta(alpha=10, beta=4)
Please rerun this cell to show the HTML repr or trust the notebook.
" + ], + "text/plain": [ + "Beta(alpha=10, beta=4)" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "B._posterior" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A plot of this posterior shows that it skews to the right with a mean of around 0.71, reflecting our updated belief that the coin is more likely to land on heads than tails. This skewness arises from the combination of our prior belief and the observed data, which had a higher number of heads compared to tails." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.7142857142857143" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "B._posterior.mean()" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAG0CAYAAADO5AZFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABW/ElEQVR4nO3deVzUdeI/8NccMIDAIHILCsihggeCB5qmoph2aFnbdtdWv7U8Kte27Nxqd22rb5qVdqiVqdkWVm4eiSl4HyieHKJyiZwqwz3DzHx+f3CUiQrIzHuO1/PxmMcu4wy+5hPOvHh/3p/3WyZJkgQiIiIiGyEXHYCIiIioK7HcEBERkU1huSEiIiKbwnJDRERENoXlhoiIiGwKyw0RERHZFJYbIiIisiksN0RERGRTWG6IiIjIprDcEBERkU0RWm6WLl2KgQMHwt3dHe7u7oiPj8emTZuu+viUlBTIZLIrbllZWWZMTURERJZMKfIvDwwMxNtvv42wsDAAwFdffYWpU6ciPT0dUVFRV31ednY23N3dW7/29vZu999pNBpx/vx5uLm5QSaTdT48ERERmY0kSaiurkZAQADk8uuMzUgWpnv37tKyZcva/LPt27dLAKRLly51+vsXFhZKAHjjjTfeeOONNyu8FRYWXvezXujIze8ZDAZ89913qK2tRXx8/DUfGxMTg4aGBvTv3x+vvPIKxo0bd9XHarVaaLXa1q+l5k3QCwsLLxv9ISIiIstVVVWFoKAguLm5XfexwsvN8ePHER8fj4aGBri6uuKHH35A//7923ysv78/PvvsM8TGxkKr1eLrr79GQkICUlJSMGbMmDafs2DBArzxxhtX3N8yz4eIiIisR3umlMiklqEMQXQ6HQoKClBZWYmkpCQsW7YMqampVy04f3T77bdDJpNh/fr1bf75H0duWpqfRqNhuSEiIrISVVVVUKvV7fr8Fj5y4+jo2DqhOC4uDgcPHsQHH3yATz/9tF3PHzFiBFatWnXVP1epVFCpVF2SlYiIiCyfxa1zI0nSZSMt15Oeng5/f38TJiIiIiJrInTk5qWXXsLkyZMRFBSE6upqrF27FikpKdi8eTMAYP78+SgqKsLKlSsBAIsWLUJwcDCioqKg0+mwatUqJCUlISkpSeTLICIiIgsitNyUlpbioYceQnFxMdRqNQYOHIjNmzdj4sSJAIDi4mIUFBS0Pl6n02HevHkoKiqCs7MzoqKisGHDBkyZMkXUSyAiIiILI3xCsbl1ZEISERERWYaOfH5b3JwbIiIiohvBckNEREQ2heWGiIiIbArLDREREdkUlhsiIiKyKSw3REREZFOEb79ARERi6PRG1Gj10BuMkMlk8HBxgIOCv/OS9WO5ISKyA5q6Ruw8XY70gkocO1eJ/At1KK/R4o8rnXm4OCDcxxX9/d0xNMQTo8O9oXZ2EBOaqJO4iB8RkY2q0+nxv6Pnsf7oeew7exEGY9tv93IZIAFXFB0AUMhlGBHqiT/FBeGWaD+olArThia6CqvaFZyIiLpWWVUDPkk9i+/SClGt1bfeH+HriuEhPTA4yAPhvq4I8HBGdxdHKOQyGIwSKut0KK3SIru0CieKqrDjVDlyymqw+/QF7D59AZ7dHPHE6BA8Eh+Mbip+fJDl4sgNEZGNqKjRYmnKGazalw+t3ggA6N3DBfcODcKUaH8Ee3Xr8PcsvFiH7w+dw3/TClGsaQAAdHdxwNzESNw/rBcUclmXvgaiq+nI5zfLDRGRlTMaJXxzsAD/2ZSFqoamkZrY3t0xe3wYxoR7Q94FBURvMGL90fNY/GsO8i7UAQAGBarx77sGICpAfcPfn+h6WG6ugeWGiGxJTmk1nv/+GI4UVgIA+vu744XJfTEm3AsyWdePqugNRqzeX4D3fslGtVYPB4UMf5/UF4/fFNIlJYroalhuroHlhohsgSRJWLW/AP/8OQNavRGuKiXmTozAw/G9oTTD5dxlVQ149acT+OVkKQBgTIQ3PvxzDNQuvLKKTIPl5hpYbojI2lU1NOJv/z2K5IymYjE63Avv3j0Ifmons+aQJAlrDhTgrZ8z0NBoRIhXNyx7JA59vF3NmoPsQ0c+v7laExGRFTlbXoNpH+9GckYpHBQyvHJrP3z12DCzFxsAkMlkeGB4b6x7ahR6ejgjt6IWd368GwdyL5o9C9HvsdwQEVmJHafKMfXj3ThbXgt/tROSnhqJJ0aHCp/r0j/AHT/OHIUhvTxQ1aDHwyv2IyW7TGgmsm8sN0REVuDH9CL85cuDqG7QI7Z3d/w0axQGBnqIjtXK202FNU+OwLhIbzQ0GvHkyjT8crJEdCyyUyw3REQW7ovduXj22yPQGyVMGxyANU8Oh4+b+U9DXY+TgwKfPhSHWwf6o9EgYfaadOzMKRcdi+wQyw0RkYWSJAkLk0/hjf9lAAAeHRmM9/802KK3QHBUyvHBvYMxOdoPOoMR/2/lIRzKvyQ6FtkZlhsiIgu1cGsOPvg1BwAwd2IEXr+9v/D5Ne2hVMix6M+DMTrcC/WNBjz2xQGcKq0WHYvsCMsNEZEF+nj7aSxuLjav3NoPcxLCTbIon6molAp8+lAs4np3R1WDHo9/dRAXa3WiY5GdYLkhIrIwn+84i3d/yQYAvHBLXzwxOlRwos5xcVTis4fj0MvTBYUX6zHj60PQNe95RWRKLDdERBbkv2mF+NfGTADAcxMi8NTYPoIT3RjPbo5Y/kgc3FRKHMi7iFd+PA47WzuWBGC5ISKyECnZZZi/7jgAYMbNfTAnIUxwoq4R7uuGxffHQC4D/pt2Dt+lnRMdiWwcyw0RkQU4fk6Dp1cfhsEo4a6YnnjhlkirmmNzPeMiffC3xEgAwKs/nUBWSZXgRGTLWG6IiAQ7d6kOj315EHU6A24K88Lb0wfaVLFp8dTNfTAmwhtavREzVx9GrVYvOhLZKJYbIiKB6nR6PLnyECpqtOjr54alDw6Bo9I235rlchkW/mkQfN1VOFNei9d+Oik6Etko2/wXRERkBSRJwvPfH0NmcRW8XB2x4tGhcHNyEB3LpHq4qvDhfUMglwFJh89xiwYyCZYbIiJBlqScwYZjxXBQyLD0wVgEeDiLjmQWw0I88f/GNF0F9tK647hQoxWciGwNyw0RkQDbskrx3pamtWz+cUcUhgZ7Ck5kXs9NDEekrxsu1Orw8g8neHk4dSmWGyIiMyu4UIdn1h6BJAEPDO+FB4b3Fh3J7FRKBf7vT4OglMuw+WQJ1h89LzoS2RCWGyIiM9LqDZj1zWFUN+gxpJcHXr89SnQkYaJ7qjEnIRwA8NbPGais4/YM1DVYboiIzGjBxiwcO6eBh4sDPrrfdq+Maq8ZN/dBmI8rKmp0eHtTlug4ZCPs+18VEZEZbT5Rgi/35AEA/u+eQXYzgfhaHJVyLLhrAABg7cFCHMi9KDgR2QKWGyIiMyi8WIfnvz8KAPjrmFAk9PMVnMhyDA32xH3DggAAL/1wnJtr0g1juSEiMjG9wYhn1qa3zrOZNylSdCSL8+It/eDl6ojTZTVYuTdPdByyciw3REQm9umOszhcUAk3lRKL74uBg4JvvX+kdnHA3yf1BQB8sDUHFVz7hm4A/4UREZnQ8XMaLEw+BQB4Y2oUAru7CE5kue6ODUR0T3dUa/X4vy2nRMchKya03CxduhQDBw6Eu7s73N3dER8fj02bNl3zOampqYiNjYWTkxNCQ0PxySefmCktEVHHNDQa8Oy36dAbJUwZ4Ic7Y3qKjmTR5HIZXrut6dL4tQcLcPK8RnAislZCy01gYCDefvttpKWlIS0tDePHj8fUqVNx8mTbm6nl5uZiypQpGD16NNLT0/HSSy9hzpw5SEpKMnNyIqLre3tTFs6U18LHTYV/TRtgkzt9d7VhIZ64baA/JAl4838ZXLmYOkUmWdhPjqenJ9599108/vjjV/zZCy+8gPXr1yMzM7P1vhkzZuDo0aPYu3dvu75/VVUV1Go1NBoN3N3duyw3EdHv7cwpx0PLDwAAvvrLMNwc4S04kfUoqqzH+PdSoNUb8cmDsbgl2k90JLIAHfn8tpg5NwaDAWvXrkVtbS3i4+PbfMzevXuRmJh42X2TJk1CWloaGhsbzRGTiOi6arR6vJh0HADw0IjeLDYd1NPDGU+ODgUAvLclGwajRf0OTlZAeLk5fvw4XF1doVKpMGPGDPzwww/o379/m48tKSmBr+/la0P4+vpCr9ejoqKizedotVpUVVVddiMiMqX/bMpCUWU9gjydMX9KX9FxrNL/uzkUHi4OOF1Wg3WHz4mOQ1ZGeLmJjIzEkSNHsG/fPjz11FN45JFHkJGRcdXH//GcdctZtaudy16wYAHUanXrLSgoqOvCExH9wYHci/h6Xz4A4O27BsLFUSk4kXVyd3LAUzf3AQAs2poDrd4gOBFZE+HlxtHREWFhYYiLi8OCBQswaNAgfPDBB20+1s/PDyUlJZfdV1ZWBqVSiR49erT5nPnz50Oj0bTeCgsLu/w1EBEBTVdHvZB0DADw56FBGBXmJTiRdXtkZDB83VUoqqzH6n0FouOQFRFebv5IkiRotW0v3hQfH4/k5OTL7tuyZQvi4uLg4ODQ5nNUKlXrpeYtNyIiU1i49RRyK2rh667C/Cn9RMexek4OCjyTEAEA+Hj7adRo9YITkbUQWm5eeukl7Ny5E3l5eTh+/DhefvllpKSk4IEHHgDQNOry8MMPtz5+xowZyM/Px9y5c5GZmYkVK1Zg+fLlmDdvnqiXQEQEADh2rhKf7zgLAPjntAFQO7f9Cxd1zD1xgQjx6oYLtTp8uTtXdByyEkLLTWlpKR566CFERkYiISEB+/fvx+bNmzFx4kQAQHFxMQoKfhuKDAkJwcaNG5GSkoLBgwfjrbfewuLFizF9+nRRL4GICHqDES8mHYdRAm4fFICJ/bkpZldxUMjx7IRwAMCyXbkcvaF2sbh1bkyN69wQUVdbsSsXb/6cAbWzA379283wclWJjmRTDEYJE99PxdmKWrxwS188NbaP6EgkgFWuc0NEZI1KNA34vy3ZAIAXbunLYmMCCrkMM8eFAQA+33kWdTqO3tC1sdwQEd2At37OQK3OgJheHvjzUC41YSpTBwegdw8XXKzV8copui6WGyKiTko9VY4Nx4shlwH/nBYNuZx7R5mKUiFvHb35dMdZ1Ou47g1dHcsNEVEnNDQa8NpPJwAAj44MQVSAWnAi23dnTE8EdndGRY0W3xzg6A1dHcsNEVEnLNl+GvkX6uDn7oS5iRGi49gFB4UcT49tGr1ZtvMsGg1GwYnIUrHcEBF1UG5FLT5JbVrT5rXb+8NVxS0WzOWuIT3h5arCeU0DNhwrFh2HLBTLDRFRB/3z5wzoDEaMifDG5Gg/0XHsipODAo+O7A2gae6Nna1mQu3EckNE1AEp2WX4NasMSrkMr93W/6qb9pLpPDiiN1wcFcgsrsKu0xWi45AFYrkhImonnd6IN3/OAAA8OjIYYT6ughPZJw8XR/wprumy+0+bTw8S/R7LDRFRO63cm4ez5bXo0c0Rc5q3BCAxHr8pBAq5DLtOV+BEkUZ0HLIwLDdERO1QUaPFB1tzAADPT4qEuxM3xhQpyNMFUwb4A2hatZjo91huiIja4b1fslGt1SO6pzvuieNKxJbgr2NCAQA/HytGUWW94DRkSVhuiIiu4/g5Db5NKwQA/OP2KCi4ErFFiO6pxsg+PWAwSvh6b77oOGRBWG6IiK5BkiS88b+TkKSm/Y3igj1FR6LfeXRkMABg7cECNDRySwZqwnJDRHQNG44XIy3/EpwdFHhxcl/RcegPEvr5IrC7MyrrGvHTkSLRcchCsNwQEV2FVm/AO5uzAQD/b0wo/NXOghPRHynkMjw0omlRvy/35HNRPwLAckNEdFVf781HwcU6+Lip8NebQ0XHoau4d2gQnBzkyCyuwsG8S6LjkAVguSEiaoOmrhEfbjsNAJg7MQIujtw/ylJ5uDjizpieAICv9uSJDUMWgeWGiKgNH23Pgaa+EZG+brz02wo80jyxePPJEpznZeF2j+WGiOgPCi7U4as9TZcWz5/Sl5d+W4G+fu4YHuIJg1HC6v28LNzesdwQEf3BO79kQWcw4qYwL9wc4S06DrXTY6OCAQBrDxRCq+dl4faM5YaI6HfSCy7h52PFkMmaRm2467f1mNDPF77uKlyo1SE5o1R0HBKI5YaIqJkkSfj3xkwAwPQhgYgKUAtORB2hVMhxb/P8qDX7CwSnIZFYboiImm3JKMXBvEtwcpDjb4kRouNQJ/xpaBBkMmDPmQvIragVHYcEYbkhIgKgNxjxzuYsAMDjN4VwwT4rFdjdBWOb50mtPcDRG3vFckNEBGDd4SKcKa9FdxcH/PXmPqLj0A24f3jTisXfHTrHicV2iuWGiOxeQ6MBi7aeAgA8PTYM7k4OghPRjRgX6Q0/dydcrNXhl5OcWGyPWG6IyO6t2peP85oG+Kud8FB8b9Fx6AYpFXLcO7RpYvE3nFhsl1huiMiuVTc04uPtTdssPDshHE4OCsGJqCvcOzQIchmw9+wFnC2vER2HzIzlhojs2rKdubhU14hQ726YPiRQdBzqIgEezhgX6QMA+PZgoeA0ZG4sN0RktypqtFi28ywAYF5iJJQKviXakj81n5pal14EvcEoOA2ZE/8lE5Hd+nj7adTqDBjQU43J0X6i41AXG9/XBz26OaK8WosdOeWi45AZsdwQkV06d6kOq/c1TTb9+y2R3GbBBjko5JgW0xMA8F3aOcFpyJxYbojILi3amgOdwYiRfXrgpjAv0XHIRO6Ja5pHtTWzFBdrdYLTkLmw3BCR3ckprca6w02/yf/9Fm6Oacv6+rljQE81Gg0SfjpSJDoOmQnLDRHZnfe2ZMMoAZOifDE4yEN0HDKxltEbnpqyHyw3RGRXjp/T4JeTpZDLmq6QItt3x6AAOCrkyCiuwsnzGtFxyAxYbojIrryfnA0AmDa4J8J93QSnIXPwcHHExP6+ADh6Yy9YbojIbhwuuITt2eVQyGWYkxAuOg6Z0d3Np6Z+OlIEnZ5r3tg6lhsishsLk5s2x5w+pCeCvboJTkPmNCbcG77uKlyqa8SvmdxM09YJLTcLFizA0KFD4ebmBh8fH0ybNg3Z2dnXfE5KSgpkMtkVt6ysLDOlJiJrdCD3InbmVEApl2H2eI7a2BuFXIa7mrfXWJfOq6ZsndByk5qaipkzZ2Lfvn1ITk6GXq9HYmIiamtrr/vc7OxsFBcXt97Cw/lmRURX1zLX5k9DgxDk6SI4DYlwZ/OCfinZZais45o3tkwp8i/fvHnzZV9/8cUX8PHxwaFDhzBmzJhrPtfHxwceHh4mTEdEtmLP6QrsO3sRjgo5Zo0LEx2HBInwdUN/f3dkFFdhw/FiPDC8t+hIZCIWNedGo2m6RM/T0/O6j42JiYG/vz8SEhKwffv2qz5Oq9WiqqrqshsR2Q9JkvB/zXNt7h/eCwEezoITkUjTYgIAAD+lnxechEzJYsqNJEmYO3cubrrpJkRHR1/1cf7+/vjss8+QlJSEdevWITIyEgkJCdixY0ebj1+wYAHUanXrLSgoyFQvgYgs0I6cChzKvwSVUo6nx/YRHYcEu2NQT8hkwIG8izh3qU50HDIRmSRJkugQADBz5kxs2LABu3btQmBgYIeee/vtt0Mmk2H9+vVX/JlWq4VWq239uqqqCkFBQdBoNHB3d7/h3ERkuSRJwrSPd+PoOQ2euCkEr9zWX3QksgD3f74Pe85cwPOTIjGTpymtRlVVFdRqdbs+vy1i5Gb27NlYv349tm/f3uFiAwAjRoxATk5Om3+mUqng7u5+2Y2I7MO2rDIcPaeBs4MCMzhqQ81adgr/Mb0IFvL7PXUxoeVGkiTMmjUL69atw7Zt2xASEtKp75Oeng5/f/8uTkdE1kySJLzfPNfmkZHB8HJVCU5EluKWaD84KuXIKatBRjHnYdoioVdLzZw5E2vWrMFPP/0ENzc3lJSUAADUajWcnZsm/c2fPx9FRUVYuXIlAGDRokUIDg5GVFQUdDodVq1ahaSkJCQlJQl7HURkeX45WYKT56vgqlLir2NCRcchC+Lu5ICJ/Xyx4XgxfkwvQlSAWnQk6mJCR26WLl0KjUaDsWPHwt/fv/X27bfftj6muLgYBQUFrV/rdDrMmzcPAwcOxOjRo7Fr1y5s2LABd911l4iXQEQWyGiUsDC56VT1X0YFo3s3R8GJyNJMHdx01dT6o+dhMPLUlK2xmAnF5tKRCUlEZJ3+d/Q8Zn+TDjcnJXb9fTzULg6iI5GF0emNGPbvraisa8TqJ4ZjVJiX6Eh0HVY3oZiIqKsYjRIW/9o0avPk6FAWG2qTo1KOKQOa5mr+yO0YbA7LDRHZlE0nSpBTVgN3JyUeHRUsOg5ZsJbtGDadKEFDo0FwGupKLDdEZDOMRgkfbmuea3NTCNydOGpDVxfbqzt6ejijRqtHSna56DjUhVhuiMhmbMkoRVZJNdxUSjw2snNLS5D9kMtluHVg06mpn49xOwZbwnJDRDZBkn6ba/PoqGDOtaF2ua253PyaWYY6nV5wGuoqLDdEZBN+zSxDRnEVujkq8JdRHLWh9hnQU41eni6obzRgW1aZ6DjURVhuiMjqSZKED5pHbR4ZyXVtqP1kMlnr6M3PR4sFp6GuwnJDRFYvJbscx4s0cHFU4InRXI2YOqZl3s327DLUaHlqyhaw3BCRVfv9qM1DI3rDk6M21EH9/d0R6tUNWr0RWzNKRcehLsByQ0RWbWdOBY4UVsLJQc5RG+qUy05N8aopm8ByQ0RW6/ejNg8M7w1vN+78TZ1z26CmvaZST5VDU98oOA3dKJYbIrJae85cwKH8S1Ap5dz5m25IhK8bInxd0WiQsOVkieg4dINYbojIarWM2tw3rBd83J0EpyFrd9vAptGbn4/xqilrx3JDRFZp39kLOJB7EY4KOWbc3Ed0HLIBLfNudp+uwKVaneA0dCNYbojIKrWsRnzv0CD4qTlqQzcu1NsV/f3doTdK2MxTU1aN5YaIrM7BvIvYc+YCHBQyzBjLURvqOrcNahq92Xicp6asGcsNEVmdllGbu2OD0NPDWXAasiWTo5vKzd4zF1BZx1NT1orlhoisyqH8S9iZUwGlXIanOWpDXSzEqxv6+rlBb5SQzAX9rBbLDRFZlQ+3NY3aTB8SiCBPF8FpyBbdEu0HANh8gvNurBXLDRFZjaOFlUjJLodCLsPT4zhqQ6bRcmpqZ04Fqhu4oJ81YrkhIqvRMmozbXBP9O7RTXAaslURvq4I9e4GncGIbVllouNQJ7DcEJFVOFGkwdbMMshlwKzxYaLjkA2TyWSYzFNTVo3lhoisQssVUlMH90SIF0dtyLRaTk1tzy5DnU4vOA11FMsNEVm8jPNV2JJRCpkMmDmOozZkelEB7gjs7oyGRiNSs8tFx6EOYrkhIov30famUZvbBgYgzMdVcBqyB78/NbWJp6asDssNEVm07JJqbDze9OEym3NtyIwmD2g6NbUtqwxavUFwGuoIlhsismgtV0hNGeCHCF83wWnIngwO9ICfuxNqtHrsyqkQHYc6gOWGiCzW6bJqbGje42f2+HDBacjeyOWy1gX9eGrKurDcEJHF+mjbaUgSMCnKF/383UXHITvUUm6SM0rRaDAKTkPtxXJDRBbpbHkN1h89D4CjNiTO0GBPeLk6QlPfiL1nLoiOQ+3EckNEFumj7adhlIAJ/XwQ3VMtOg7ZKYVchsSo5gX9TvLUlLVguSEii5NXUYufjjSN2sxJ4KgNiTWpudxszSiF0SgJTkPtwXJDRBZnScppGIwSxkV6Y2Cgh+g4ZOfiQ3vATaVEWbUWR89Vio5D7cByQ0QWpfBiHdYdLgIAzOaoDVkAR6UcY/v6AAC2ZJQKTkPtwXJDRBZlScoZ6I0SRod7YUiv7qLjEAEAEvv7AgC2cN6NVWC5ISKLUVRZj+8PFQIAnuGoDVmQsZHecFDIcKa8FqfLakTHoetguSEii7E05TQaDRJG9umBuGBP0XGIWrk5OWBkHy8ATWvekGVjuSEii1Csqcd/D54DwCukyDIlRjWfmsrgqSlLx3JDRBbh09Sz0BmMGB7iiRGhPUTHIbrCxH5N5Sa9oBJlVQ2C09C1CC03CxYswNChQ+Hm5gYfHx9MmzYN2dnZ131eamoqYmNj4eTkhNDQUHzyySdmSEtEplJa1YA1BwoAcNSGLJePuxNienkAAJIzeWrKkgktN6mpqZg5cyb27duH5ORk6PV6JCYmora29qrPyc3NxZQpUzB69Gikp6fjpZdewpw5c5CUlGTG5ETUlT5NPQud3oi43t0xsg9HbchyJfZvWtBvy0mWG0smkyTJYpZbLC8vh4+PD1JTUzFmzJg2H/PCCy9g/fr1yMzMbL1vxowZOHr0KPbu3Xvdv6OqqgpqtRoajQbu7tyIj0i0suoGjP7Pdmj1Rnz9+DCMDvcWHYnoqk6X1WDC+6lwUMhw+NWJcHNyEB3JbnTk89ui5txoNBoAgKfn1a+S2Lt3LxITEy+7b9KkSUhLS0NjY6NJ8xFR1/t8x1lo9UbE9PLATWFeouMQXVOYjytCvbuh0SAhJbtcdBy6CospN5IkYe7cubjpppsQHR191ceVlJTA19f3svt8fX2h1+tRUVFxxeO1Wi2qqqouuxGRZaio0WLVvt/m2shkMsGJiK6v9dQULwm3WBZTbmbNmoVjx47hm2++ue5j//gG2HJmra03xgULFkCtVrfegoKCuiYwEd2wz3eeRX2jAYMC1RgbwdNRZB1aLgnfnlUGrd4gOA21xSLKzezZs7F+/Xps374dgYGB13ysn58fSkouX2OgrKwMSqUSPXpcORFx/vz50Gg0rbfCwsIuzU5EnXOxVoev9+YD4KgNWZfBgR7wdlOhRqvHvrMXRcehNggtN5IkYdasWVi3bh22bduGkJCQ6z4nPj4eycnJl923ZcsWxMXFwcHhyoldKpUK7u7ul92ISLzlu86iTmdAdE93jG/elJDIGsjlMkzkXlMWTWi5mTlzJlatWoU1a9bAzc0NJSUlKCkpQX19fetj5s+fj4cffrj16xkzZiA/Px9z585FZmYmVqxYgeXLl2PevHkiXgIRdUJlnQ5f7WketRnPURuyPi0baSZnlMJotJiLjqmZ0HKzdOlSaDQajB07Fv7+/q23b7/9tvUxxcXFKCgoaP06JCQEGzduREpKCgYPHoy33noLixcvxvTp00W8BCLqhBW7clGj1aOfv3vrb8BE1iS+Tw+4qpQoq9biyLlK0XHoD5Qi//L2LLHz5ZdfXnHfzTffjMOHD5sgERGZmqauEV/szgMAPJMQxlEbskoqpQJjI73x87FibM0oxZBe3UVHot+xiAnFRGQ/vtiTi2qtHpG+bq2X1BJZownNe039mlkmOAn9EcsNEZlNVUMjVuzKBQDMTgiDXM5RG7JeYyO9oZDLkF1ajcKLdaLj0O+w3BCR2Xy1Ow9VDXqE+7hiSrS/6DhEN8TDxRFxvZtOR23lRpoWheWGiMyiRqvHsuZRm1njOWpDtqFlQjzLjWVhuSEis1i5Nw+a+kaEenfDbQMDRMch6hIJzfNu9p+9iKoG7m9oKVhuiMjkarV6fL7jLABg9vgwKDhqQzYixKsb+nh3g94oIZUbaVoMlhsiMrlV+/Jxqa4RwT1ccDtHbcjG/HbVFE9NWQqWGyIyqTqdHp81j9rMGh8OpYJvO2RbJjTPu9meXQ69wSg4DQGdWMRPkiSkpqZi586dyMvLQ11dHby9vRETE4MJEyZw120iusya/QW4UKtDL08XTB3MURuyPUN6dUd3FwdcqmtEWv4ljAi9chNnMq92/wpVX1+Pf//73wgKCsLkyZOxYcMGVFZWQqFQ4PTp03j99dcREhKCKVOmYN++fabMTERWoqHRgE9Sm0ZtZo7rAweO2pANUshlGNe8+evWDJ6asgTtfqeJiIjA4cOH8cknn6Cqqgr79u1DUlISVq1ahY0bN6KgoABnzpzB6NGjce+99+Lzzz83ZW4isgJr9hegokaLnh7OuGtIoOg4RCbTOu8mi6sVW4J2n5batGkToqOjr/mY3r17Y/78+fjb3/6G/Pz8Gw5HRNaradTmDABg5rgwjtqQTRsT4Q1HhRy5FbU4U16DPt6uoiPZtXa/21yv2Pyeo6MjwsPDOxWIiGzD6v0FKKtuGrW5O5ajNmTbXFVKDA/1BMBTU5agU79KvfrqqzAYDFfcr9FocN99991wKCKybvU6A5amNI3azBofBkclR23I9nEjTcvRqXeclStXYtSoUThz5kzrfSkpKRgwYADy8vK6KhsRWalV+/JRUaNFkCdHbch+JPRrmlScln8Rl2p1gtPYt06Vm2PHjiE4OBiDBw/G559/jueffx6JiYl49NFHsWvXrq7OSERWpE6nb51rM3tcOOfakN0I7O6Cvn5uMErA9myO3ojU4XVuAECtVmPt2rV4+eWX8de//hVKpRKbNm1CQkJCV+cjIivz9d58XKjVoXcPF9w5pKfoOERmNbG/L7JKqrE1s5RXCArU6V+pPvzwQyxcuBD33XcfQkNDMWfOHBw9erQrsxGRlanV6vFp6x5SHLUh+9OykeaOUxXQ6q+cm0rm0al3nsmTJ+ONN97AypUrsXr1aqSnp2PMmDEYMWIE3nnnna7OSERW4qu9ebhYq0OIVzdM42rEZIcG9lTD202FGq0e+89eFB3HbnWq3Oj1ehw7dgx33303AMDZ2RlLly7F999/j4ULF3ZpQCKyDtUNja17SM1JCOMeUmSX5HIZEppXK+ZGmuJ06t0nOTkZAQFX/lZ266234vjx4zccioisz1d78lBZ14hQ7264YxDn2pD9arkkfGtmGSRJEpzGPnX5r1ZeXl4AwP+gRHakqqERn+/MBQA8kxAOhVwmOBGROKPCvKBSylFUWY+skmrRcexSu8tNv379sGbNGuh01752PycnB0899RT+85//3HA4IrIOX+zKg6a+EWE+rrhtIOfakH1zdlRgdHjTL/pcrViMdl8K/vHHH+OFF17AzJkzkZiYiLi4OAQEBMDJyQmXLl1CRkYGdu3ahYyMDMyaNQtPP/20KXMTkYXQ1Ddi2a6muTbPTuCoDRHQdNXU1swybM0qw+wEbkdkbu0uN+PHj8fBgwexZ88efPvtt1izZg3y8vJQX18PLy8vxMTE4OGHH8aDDz4IDw8PE0YmIkuyfFcuqhv0iPR1w5Rof9FxiCxCy6Tio4WVKKtqgI+7k+BE9qXDi/iNHDkSI0eONEUWIrIylXU6fLGraa7NsxPCIeeoDREAwMfdCYMC1Th6ToNtWWX487BeoiPZlXbPufH09ERFRQUA4C9/+QuqqzlJisjeLduZi2qtHn393DApyk90HCKLkvC7q6bIvNpdbnQ6HaqqqgAAX331FRoaGkwWiogs36VaHb7Y3TRq89zECI7aEP1By0aau06Xo6GRqxWbU7tPS8XHx2PatGmIjY2FJEmYM2cOnJ2d23zsihUruiwgEVmmT3ecRa3OgKgAdyT29xUdh8ji9Pd3R4DaCec1DdhzpgLj+/Lfibm0e+Rm1apVmDJlCmpqaiCTyaDRaHDp0qU2b0Rk28qqG/DlnqZRm7kTIyCTcdSG6I9kMhnGN4/e8NSUebV75MbX1xdvv/02ACAkJARff/01evToYbJgRGS5Pt52Gg2NRgzp5YHxzVeFENGVEvr5YtW+AmzLLIM0TeIvAmbSqRWKc3NzWWyI7NS5S3VYc6AAADBvUiTfrImuIT60B1wcFSipasDJ81Wi49iNdo/cLF68uN3fdM6cOZ0KQ0SW74OtOWg0SBgV1gMj+3iJjkNk0ZwcFLgpzAtbMkqxNbMU0T3VoiPZhXaXmz/u9l1eXo66urrWBfsqKyvh4uICHx8flhsiG3WmvAZJh88BAOYlRgpOQ2QdJvTzbS03z06IEB3HLrT7tFRubm7r7V//+hcGDx6MzMxMXLx4ERcvXkRmZiaGDBmCt956y5R5iUig95NPwSg1vVnH9OouOg6RVRjX1wcyGXCiqAolGi6jYg6dmnPz6quv4sMPP0Rk5G+/uUVGRmLhwoV45ZVXuiwcEVmOk+c12HCsGDIZ8LdE/vZJ1F7ebioMCvQAAPyaxY00zaFT5aa4uBiNjY1X3G8wGFBayv9wRLbo/S2nAAC3DwxAP393wWmIrMvE5rWgfuUl4WbRqXKTkJCAJ598EmlpaZAkCQCQlpaGv/71r5gwYUKXBiQi8Q7lX8KvWWVQyGV4biJHbYg6qmW14t2nK1Cv42rFptapcrNixQr07NkTw4YNg5OTE1QqFYYNGwZ/f38sW7asqzMSkUCSJOHdX7IAAPfEBiLEq5vgRETWJ9LXDT09nKHVG7HrdIXoODavw7uCA4C3tzc2btyInJwcZGZmQq/XIzo6GhER/I2OyNbsPn0B+85ehKNCjtkJ4aLjEFklmUyGCf188NXefPyaWdp6mopMo1MjNwCwfPly3Hnnnbjnnntw33334a677urwqM2OHTtw++23IyAgADKZDD/++OM1H5+SkgKZTHbFLSsrq7Mvg4iuQZIkvLslGwBw//Be6OnR9n5yRHR9LbuE/5pVBqNREpzGtnVq5ObVV1/FwoULMXv2bMTHxwMA9u7di+eeew55eXn45z//2a7vU1tbi0GDBuGxxx7D9OnT2/33Z2dnw939twmN3t7eHXsBRNQuyRmlOFpYCWcHBWaOCxMdh8iqDQ/1RDdHBcqrtThepMGgIA/RkWxWp8rN0qVL8fnnn+O+++5rve+OO+7AwIEDMXv27HaXm8mTJ2Py5Mkd/vt9fHxaFw8kItMwGCW81zxq89ioYHi7qQQnIrJuKqUCYyK8selECX7NLGW5MaFOnZYyGAyIi4u74v7Y2Fjo9fobDnU9MTEx8Pf3R0JCArZv337Nx2q1WlRVVV12I6LrSzp8DqdKa6B2dsBfb+4jOg6RTWg5NcVdwk2rU+XmwQcfxNKlS6+4/7PPPsMDDzxww6Guxt/fH5999hmSkpKwbt06REZGIiEhATt27LjqcxYsWAC1Wt16CwoKMlk+IlvR0GjAwuSmdW1mjQuD2tlBcCIi2zAu0hsyGZBRXIXzlfWi49gsmdSyUE0HzJ49GytXrkRQUBBGjBgBANi3bx8KCwvx8MMPw8HhtzfC999/v31BZDL88MMPmDZtWoey3H777ZDJZFi/fn2bf67VaqHValu/rqqqQlBQEDQazWXzdojoN5+mnsGCTVkIUDth27yxcHJQiI5EZDOmL92DQ/mX8NbUKDwUHyw6jtWoqqqCWq1u1+d3p+bcnDhxAkOGDAEAnDlzBkDTpF5vb2+cOHGi9XEymawz375DRowYgVWrVl31z1UqFVQqzhUgai9NXSM+3n4aADA3MZLFhqiLJfTzwaH8S9iaWcZyYyKdKjfXm+diTunp6fD39xcdg8hmLEk9jaoGPSJ93XBnTE/RcYhszoR+vnhnczb2nrmAWq0e3VSd+iimaxB6RGtqanD69OnWr3Nzc3HkyBF4enqiV69emD9/PoqKirBy5UoAwKJFixAcHIyoqCjodDqsWrUKSUlJSEpKEvUSiGzK+cp6fLE7DwDwwuRIKOSmH30lsjfhPq7o5emCgot12JlTgVui/URHsjlCy01aWhrGjRvX+vXcuXMBAI888gi+/PJLFBcXo6CgoPXPdTod5s2bh6KiIjg7OyMqKgobNmzAlClTzJ6dyBYt2noKOr0Rw0I8MS7SR3QcIpskk8mQ0M8HX+zOw6+ZpSw3JtCpCcXWrCMTkojsyanSatyyaAeMErDu6ZEY0qu76EhENmv36Qo8sGw/vFwdceClCZBzlPS6OvL53entF4jItryzOQtGCbglyo/FhsjEhgZ7wk2lREWNDkfOVYqOY3NYbogIB3IvYmtmGRRyGZ6/JVJ0HCKb56iUY0xk09ZBv2aWCk5je1huiOycJEl4e1MmAOBPcUHo4+0qOBGRfZjQr2le269crbjLsdwQ2blfTpbicEElnBzkeHZCuOg4RHZjbIQP5DIgq6Qa5y7ViY5jU1huiOyYTm9sHbV5/KYQ+Lo7CU5EZD+6d3NEXG9PABy96WosN0R2bNW+fORdqIOXqyOeGhsmOg6R3UloPjW1lfNuuhTLDZGdqqzT4YNfcwAAcydGwpWrpBKZXcsu4fvOXkB1Q6PgNLaD5YbITn247TQ09Y2I9HXDn+ICRcchskt9vLshuIcLGg0SduZUiI5jM1huiOxQXkUtVu7NAwC8dGs/KBV8KyASQSaTYULz6A1PTXUdvqMR2aG3N2Wh0SDh5ghv3BzhLToOkV1rOTWVkl0Og9GuNg0wGZYbIjtzIPciNp8sgVwGvHxrP9FxiOxeXHB3uDspcbFWh/SCS6Lj2ASWGyI7YjRK+OeGDADAvUN7IcLXTXAiInJQyDE2suWqKV4S3hVYbojsyPqj53HsnAbdHBWYOzFCdBwiapbQulox5910BZYbIjvR0GjAO5uzAABPjwuDt5tKcCIiajE2wgcKuQw5ZTUouMDVim8Uyw2RnVi+KxfnNQ0IUDvh8ZtCRMchot9RuzhgaHB3ALxqqiuw3BDZgRJNAz7efhoA8Pdb+sLJQSE4ERH9Ucsl4b9msdzcKJYbIjvwn81ZqNMZMKSXB6YODhAdh4ja0HJJ+P6zF1HF1YpvCMsNkY07lH8RP6QXQSYD/nFHFGQymehIRNSGEK9uCPXuBr1RQmp2ueg4Vo3lhsiGGYwS/rG+6dLvP8UGYWCgh9hARHRNraemOO/mhrDcENmw79IKcbxIAzeVEs/fEik6DhFdR0u52Z5dDr3BKDiN9WK5IbJRmvpGvPtLNgDgmQnh8HLlpd9Elm5ILw94uDhAU9+IQ/lcrbizWG6IbNQHW3NwoVaHMB9XPDIyWHQcImoHpUKOcc2rFf+axdWKO4vlhsgG5ZRWt+76/dpt/eHAXb+JrEbLasVc76bz+I5HZGMkScKbP2dAb5Qwsb8vxnDXbyKrMibCG0q5DGfLa5FbUSs6jlViuSGyMckZpdiZUwFHhRyvcNdvIqvj7uSA4aGeAHjVVGex3BDZkHqdAW/8r+nS7ydGh6B3j26CExFRZyT0bbpqiqemOoflhsiGfLgtB0WV9ejp4YxZ48NExyGiTmq5JPxg3iVo6rhacUex3BDZiNNl1fh851kAwOu394eLo1JwIiLqrF49XBDu4wqDUULKKV411VEsN0Q2QJIkvPrjSTQaJCT09cHE/r6iIxHRDWrZa2prJstNR7HcENmA9UfPY+/ZC1Ap5dw/ishGTGi+JDwluwyNXK24Q1huiKxcVUMj3vo5EwAwe3wYgjxdBCcioq4Q06s7PLs5orpBj4N5F0XHsSosN0RW7v0tp1BRo0Wodzc8OSZUdBwi6iIKuey31Yp5aqpDWG6IrNiJIk3rSsRvTY2GSqkQG4iIutTE/k3lZktGCSRJEpzGerDcEFkpo1HCyz+egFEC7hgUgFFhXqIjEVEXGxPhDZVSjsKL9cgsrhYdx2qw3BBZqTUHCnC0sBKuKiVXIiayUS6OytYtVH45WSI4jfVguSGyQiWaBvxnUxYA4G+JEfBxdxKciIhMZVKUHwCWm45guSGyQq/9dALVWj0GB3ng4fhg0XGIyIQm9POBQi5DVkk18i9wI832YLkhsjKbTxRjS0YplHIZ3p4+AAo517QhsmUeLo4YHtK0kSZHb9qH5YbIimjqG/HqTycBADNu7oO+fu6CExGROfx2aoobabaH0HKzY8cO3H777QgICIBMJsOPP/543eekpqYiNjYWTk5OCA0NxSeffGL6oEQW4u1NWSiv1iLUqxs3xiSyI4lRTVsxHC64hLLqBsFpLJ/QclNbW4tBgwbho48+atfjc3NzMWXKFIwePRrp6el46aWXMGfOHCQlJZk4KZF4+85ewDcHCgAAC+4aACcHrmlDZC/81c4YFKiGJAHJGRy9uR6h2wZPnjwZkydPbvfjP/nkE/Tq1QuLFi0CAPTr1w9paWl47733MH36dBOlJBKvodGAl9YdBwDcNywIw0N7CE5EROaWGOWHo+c0+OVkKR4Y3lt0HItmVXNu9u7di8TExMvumzRpEtLS0tDY2CgoFZHpfbTtNM5W1MLbTYUXJ3NNGyJ7dEt007ybPacroKnnZ961WFW5KSkpga+v72X3+fr6Qq/Xo6Kios3naLVaVFVVXXYjsiYnijRYmnoGAPDmHVFQOzsITkREIvTxdkWYjyv0Rgnbs7jX1LVYVbkBAJns8steW/ba+OP9LRYsWAC1Wt16CwoKMnlGoq6i1Rsw77ujMBglTBngh8kD/EVHIiKBJjVPLOYl4ddmVeXGz88PJSWX/wctKyuDUqlEjx5tz0GYP38+NBpN662wsNAcUYm6xIe/nkZWSTV6dHPEW1OjRcchIsFaLglPyS5HQ6NBcBrLZVXlJj4+HsnJyZfdt2XLFsTFxcHBoe2hepVKBXd398tuRNbg2LnK1tNRb02LRg9XleBERCTagJ5qBKidUN9owM6ctqdjkOByU1NTgyNHjuDIkSMAmi71PnLkCAoKmi53nT9/Ph5++OHWx8+YMQP5+fmYO3cuMjMzsWLFCixfvhzz5s0TEZ/IZH5/Ouq2gf6YwtNRRISmKRiJ3GvquoSWm7S0NMTExCAmJgYAMHfuXMTExOC1114DABQXF7cWHQAICQnBxo0bkZKSgsGDB+Ott97C4sWLeRk42ZwPtubgVGkNvFwd8SZPRxHR77Qs6Lc1sxR6g1FwGsskk1pm5NqJqqoqqNVqaDQanqIii3SksBJ3LdkNowR88mBs6+WfREQAoDcYMfRfW3GprhFrnhiOkWFeoiOZRUc+v61qzg2RrWtobDodZZSAqYMDWGyI6ApKhRwT+jWN3mzmqak2sdwQWZB3NmfjdFkNvFxV+MftUaLjEJGFmjyg6RefTSdKYDDa1QmYdmG5IbIQO3PKsWJ3LgDgP9MHoHs3R8GJiMhS3RTmDTcnJcqrtUjLuyg6jsVhuSGyAJdqdZj33VEAwIMjeiGhn+91nkFE9sxRKUdi/6bRm43HiwWnsTwsN0SCSZKE+euOo7RKiz7e3fDylP6iIxGRFbh1IE9NXQ3LDZFg3x06h80nS6CUy/DBn2Pg7KgQHYmIrEDLqakynpq6AssNkUD5F2rxxvqTAIC5iRGI7qkWnIiIrIWjUo6J/ZtOYfPU1OVYbogE0RuMePbbI6jVGTAsxBN/HdNHdCQisjK3Nq9evulECYw8NdWK5YZIkA+3nUZ6QSXcnJRYeO9gKORt72xPRHQ1N4V7wU3VfGoq/5LoOBaD5YZIgL1nLuDDbTkAgH9Oi0ZPD2fBiYjIGqmUCkyM4qmpP2K5ITKzihotnlmbDqME3BMbiKmDe4qORERW7LdTU8U8NdWM5YbIjIxGCc99ewRl1VqE+7jijalchZiIbkzLqanSKi0OFfDUFMByQ2RWn+w4g505FXBykOPjB4bAxVEpOhIRWTmVUtF61dSGYzw1BbDcEJlNWt5F/N+WUwCAN+6IQoSvm+BERGQrpvDU1GVYbojM4FKtDrO/SYfBKGHq4AD8KS5IdCQisiGjI7zg5tR0auoAF/RjuSEyNaNRwtz/HkGxpgEhXt3wrzsHQCbjZd9E1HVUSgUmRzdtx/DTkfOC04jHckNkYou35WB7djlUSjk+uj8GrirOsyGirtdy5eWmE8XQ6Y2C04jFckNkQtuzyvDBr03r2fzrzgGICuD2CkRkGiNCe8DbTYXKukbszCkXHUcolhsiEym4UIdn1qZDkoAHR/TC3bGBoiMRkQ1TyGW4bWDTxGJ7PzXFckNkAvU6A/666hCqGvQYHOSBV2/rLzoSEdmBllNTyRmlqNPpBacRh+WGqItJkoSXfziOzOIqeLk6YumDQ6BSKkTHIiI7MChQjd49XFDfaEByRqnoOMKw3BB1sZV787EuvQgKuQwf3jcE/mruG0VE5iGTyXDHoAAAwHo7PjXFckPUhXafrsCbP2cAAF68pS/i+/QQnIiI7M3UwU3lJvVUOS7V6gSnEYPlhqiL5FbU4unVh2EwSrgrpieeGB0iOhIR2aEwHzf093eH3ihh04kS0XGEYLkh6gKa+kY8/tVBaOobEdPLA/++iwv1EZE4LaM3Px0pEpxEDJYbohukNxgxa81hnC2vRYDaCZ8+FAsnB04gJiJxbm+ed7M/9yKKKusFpzE/lhuiG/SvjZnYmVMBZwcFPns4Dj5uTqIjEZGdC/BwxohQTwDAj+n2N3rDckN0A9bsL8AXu/MAAO//aRCie3IFYiKyDNOHNC0c+v2hc5Ak+9opnOWGqJO2Z5Xh1Z9OAADmTozA5AH+ghMREf1m8gB/ODsokFtRi8MFlaLjmBXLDVEnHDtX2Xpl1N2xgZg9Pkx0JCKiy7iqlJg8oGmn8KTD5wSnMS+WG6IOKrxYh798eRD1jQaMDvfCAl4ZRUQW6u7mU1P/O3oeDY0GwWnMh+WGqAMu1erwyBcHUFGjQz9/dyx5YAgcFPxnRESWaURoD/T0cEZ1g96utmPguzJROzU0GvDkyrTWS76/eHQo3JwcRMciIroquVyGu4Y0bab5/SH7OTXFckPUDnqDEXO+SUda/iW4OSnx5V+GwU/NS76JyPLd1XxqamdOOUqrGgSnMQ+WG6LrMBol/D3pGLZklMJRKcdnD8UhwtdNdCwionYJ8eqGuN7dYZSAH+xkzRuWG6JrkCQJb/6cgXWHm3b5/ui+GG6GSURWZ3ps0+hNkp2secNyQ3QNC7fm4Ms9eQCA9+4ZiMQoP7GBiIg64daB/lAp5cgpq8GRwkrRcUyO5YboKpbtPIvFv+YAAN6cGoU7YwIFJyIi6hx3Jwfc2rzQ6DcHCgSnMT2WG6I2fHOgAP/ckAkAmJcYgYfjg8UGIiK6QfcN7wUA+N/RYlQ1NApOY1osN0R/sPZAAeavOw4AeHJ0CGaO4+rDRGT94np3R7iPK+obDfjJxicWCy83S5YsQUhICJycnBAbG4udO3de9bEpKSmQyWRX3LKyssyYmGzZ2gMFeLG52Dw6MhgvTenH1YeJyCbIZDLcN6xp9Gb1/gKbnlgstNx8++23ePbZZ/Hyyy8jPT0do0ePxuTJk1FQcO3zgdnZ2SguLm69hYeHmykx2bJvD15ebF6/vT+LDRHZlLuG9ISjUo6skmqbnlgstNy8//77ePzxx/HEE0+gX79+WLRoEYKCgrB06dJrPs/Hxwd+fn6tN4VCYabEZKv+e7CQxYaIbJ6HiyNus4OJxcLKjU6nw6FDh5CYmHjZ/YmJidizZ881nxsTEwN/f38kJCRg+/bt13ysVqtFVVXVZTei31uzvwAvrDsGSWKxISLbZw8Ti4WVm4qKChgMBvj6+l52v6+vL0pKStp8jr+/Pz777DMkJSVh3bp1iIyMREJCAnbs2HHVv2fBggVQq9Wtt6CgoC59HWTdPkk9g5d+OM5iQ0R2wx4mFgufUPzHDxJJkq764RIZGYknn3wSQ4YMQXx8PJYsWYJbb70V77333lW///z586HRaFpvhYWFXZqfrJMkSfjP5iy8valpMvrTY/uw2BCRXfj9xOJV+2xzYrGwcuPl5QWFQnHFKE1ZWdkVoznXMmLECOTk5Fz1z1UqFdzd3S+7kX0zGiW88uMJLE05AwB4cXJf/P2Wviw2RGQ3pscGwtlBgezSauw9e0F0nC4nrNw4OjoiNjYWycnJl92fnJyMkSNHtvv7pKenw9/fv6vjkY1qNBjx3H+PYPX+AshkwL/vHIAZN/cRHYuIyKzUzg6YHtsTAPDl7jyxYUxAKfIvnzt3Lh566CHExcUhPj4en332GQoKCjBjxgwATaeUioqKsHLlSgDAokWLEBwcjKioKOh0OqxatQpJSUlISkoS+TLISlQ3NOLp1YexM6cCSrkM7987GHcMChAdi4hIiEfig7FqXwG2Zpai8GIdgjxdREfqMkLLzb333osLFy7gzTffRHFxMaKjo7Fx40b07t0bAFBcXHzZmjc6nQ7z5s1DUVERnJ2dERUVhQ0bNmDKlCmiXgJZiWJNPR774iCySqrh7KDAkgeGYFxfH9GxiIiECfd1w+hwL+zMqcDX+/Lx0pR+oiN1GZlkizOJrqGqqgpqtRoajYbzb+xEZnEVHvviIEqqGuDlqsKKR+MwMNBDdCwiIuF+zSzF41+lwd1JiX0vJcDFUeiYxzV15PNb+NVSRKa041Q57vlkL0qqGhDm44ofnh7JYkNE1GxcpA9693BBVYMe6w7bzmXhLDdks1bvz8dfvjyIGq0ew0M8kTRjpE2dUyYiulFyuQyPxAcDAL7ckwej0TZO5rDckM3R6Y146YfjePmHE9AbJUwbHICVjw+D2sVBdDQiIotzd1wg3FRKnC6rwbasMtFxugTLDdmU8motHli2D2uaL/X++y2RWHjvYKiU3H+MiKgt7k4OuH9E06J+n6SeEZyma7DckM04UaTB1I924WDeJbiplFjxyFA8PTaMi/MREV3H46NC4KiQIy3/EtLyLoqOc8NYbsgmfJdWiOlL9+C8pgGh3t3w46xRvNSbiKidfNydcNeQpkX9bGH0huWGrFqdTo+//fconv/+GLR6I8ZFeuPHmaPQx9tVdDQiIqvy/8aEQiYDtmaW4VRpteg4N4TlhqxWTmk1pn60G0mHz0EuA56fFInljwyFuxMnDhMRdVSotysm9fcDAHyaelZwmhvDckNW6ftD53DHR7uRU1YDHzcVvnlyBGaOC4Nczvk1RESdNWNs0157Px0pQuHFOsFpOo/lhqxKZZ0Os79Jx7zvjqK+0YDR4V7Y+MxoDA/tIToaEZHVGxzkgZvCvKA3Svh4+2nRcTqN5Yasxs6cckxatAP/O3oeCrkM8xIj8NVjw+DlqhIdjYjIZjw3MRwA8N2hcyi4YJ2jNyw3ZPHqdQb8Y/1JPLT8AEqrtAj16oZ1T43ErPHhPA1FRNTFYnt7YkyENwxGCR9uyxEdp1NYbsiiHcq/iFs/3Ikv9+QBAB6J740Nc0ZjUJCH0FxERLbsuQlNozfr0ouQV1ErOE3HsdyQRapuaMRrP53A3Z/sxdnyWvi4qfDVX4bhjanRcHbkasNERKYU06s7xvf1gcEoYfGv1jd6w3JDFmdrRikmvr8DK/fmQ5KAe2IDseW5Mbg5wlt0NCIiu/HchAgAwA9HinDyvEZwmo5huSGLUaJpwMw1h/HEyjSUVDWgl6cLVj8xHO/eMwgeLo6i4xER2ZUBgWrcPigAkgQs2JgFSbKeHcOVogMQafUGLN+Vi4+2nUadzgCFXIYnRofg2YQInoIiIhLo75Mi8cuJEuw6XYHUU+UYG2kd29qw3JBQ27JK8eb/MpDXfLnhkF4eeHNqNKJ7qgUnIyKiIE8XPDKyNz7fmYsFG7MwOtwbCiu4SpXlhoQ4VVqNtzdlYVtWGQDAx02F+VP6YtrgntzFm4jIgswaF47/pp1Ddmk1vj1YiPuH9xId6bpYbsiszlfWY2HyKSQdPgejBDgoZPjLTSGYPT4crir+OBIRWRq1iwPmJITjrZ8z8M4vWbgl2g+e3Sx7HiQ/TcgsKut0WJpyBl/syYNObwQATI72w/OTIhHKHbyJiCzaw/G98V1aIbJKqvHO5iy8PX2g6EjXxHJDJlVZp8OK3Xn4Yncuqhv0AIDhIZ54YXJfDOnVXXA6IiJqDweFHG9Ni8Y9n+zF2oOFuCcuCLG9Lfc9nOWGTOJCjRbLduVi5Z481OoMAIC+fm54YXJfjI3w5rwaIiIrMzTYE/fEBuK7Q+fwyo8n8L9Zo6BUWOaKMiw31KVKNA1YtvMsVu8vQH1jU6np5++O2ePDcEuUH/eCIiKyYi9O7ostGaXILK7CpzvOYua4MNGR2sRyQ13i+DkNlu86i5+PFUNvbFroaWCgGrPHh2NCPx+O1BAR2YAeriq8fnt/zP3vUSzaegrjIn3QP8BddKwrsNxQpxmMEpIzSrFiVy4O5F1svX9YsCeeHtcHN/P0ExGRzbkzpic2nyjBloxSzP3vEayfdRMclZZ1eorlhjqstKoB/z1YiG/TCnHuUj0AQCmX4baB/nj8plAMCOQCfEREtkomk+Ffdw5AWv4lZJVUY9HWU/j7LX1Fx7oMyw21i8EoIfVUGdbsL8T27DIYmk89ebg44P5hvfBwfDD81E6CUxIRkTl4u6nwr2nReGr1YSxNPYP4Pj0wOtxyNjdmuaFryiqpwk9HzuPH9CIUaxpa7x8a3B33DeuFKQP84eTA/Z+IiOzN5AH++PPQIKw9WIhn1h7Bhjk3wV/tLDoWAJYbakPhxTqsP3oe64+cR3Zpdev93V0cMH1IIP48LAhhPm4CExIRkSX4xx1ROHZOg4ziKsxak45vnhxhEfNvZJI17WHeBaqqqqBWq6HRaODubnkzvEXJv1CLLSdLsflkCQ7lX2q931Ehx7i+3pg6uCcS+vlApeQoDRER/Sb/Qi1uW7wL1Vo97o4NxLt3DzTJxSQd+fzmyI2dkiQJx4s0SM4oxZaTpZeN0MhkQHxoD0wdHIBbov2hdnYQmJSIiCxZ7x7dsPj+GDz+5UF8f+gcQry6CV//huXGjlyo0WLX6QrszKnArpwKlFT9NodGIZdhRKgnEvv7YVKUHycHExFRu42L9MEbU6Px6o8n8O4v2QjwcMKdMYHC8rDc2LB6nQHpBZewI6cCu06X40RR1WV/7uKowNhIb0zs74vxkb5Qu3CEhoiIOuehEb2RX1GLZbty8cL3xzEitIewCcYsNzbkUq0OafmXcDDvIg7mXcSJIg0aDZdPqern747R4V4YHe6FocGevNKJiIi6zEtT+qG6QY/REV5Cr5xiubFSDY0GZJVU43iRBsfPVSK9oBI5ZTVXPM7XXYVRYU1lZlSYF3zceLqJiIhMQy6X4T93DxQdg+XGGlQ1NCKntBqZxdU4UaTBsXManCqtbt3D6ffCfFwxNLg7hgZ7YmiwJwK7O3MLBCIisissNxakqqERZ8trcaq0GqdKqnGqrAY5pdWXLZ73e57dHDGgpxoDeqoxMFCN2N7d0cNVZebUREREloXlxowkScKFWh3yL9Qh/0Ltb/97sQ75F+pwsVZ31ef6uTshws8N0QHuGBioxoBADwSonTgqQ0RE9AfCy82SJUvw7rvvori4GFFRUVi0aBFGjx591cenpqZi7ty5OHnyJAICAvD3v/8dM2bMMGPitkmShEt1jSjW1KO4sqHpfzUNzbff/r9Ob7zm9/FyVSHSzxXhPm6I8HVDpJ8rwnzcuNYMERFROwktN99++y2effZZLFmyBKNGjcKnn36KyZMnIyMjA7169bri8bm5uZgyZQqefPJJrFq1Crt378bTTz8Nb29vTJ8+XcAr+M3Zilok/F/qdR8nkwH+7k7o3aMbgr1c0MuzG4J7uKBXDxf07tENrirhfZOIiMiqCd1+Yfjw4RgyZAiWLl3ael+/fv0wbdo0LFiw4IrHv/DCC1i/fj0yMzNb75sxYwaOHj2KvXv3tuvvNNX2C3U6Pfq/9gu8XB3hr3aGn9oJAWon+Hs4w1/tBH910//6ujtZxL4bRERE1sQqtl/Q6XQ4dOgQXnzxxcvuT0xMxJ49e9p8zt69e5GYmHjZfZMmTcLy5cvR2NgIBwdxp25cHJXI/uct3HuJiIhIMGHlpqKiAgaDAb6+vpfd7+vri5KSkjafU1JS0ubj9Xo9Kioq4O/vf8VztFottFpt69dVVVVXPKarsNgQERGJJ/z8yB+v9pEk6ZpXALX1+Lbub7FgwQKo1erWW1BQ0A0mJiIiIksmrNx4eXlBoVBcMUpTVlZ2xehMCz8/vzYfr1Qq0aNHjzafM3/+fGg0mtZbYWFh17wAIiIiskjCyo2joyNiY2ORnJx82f3JyckYOXJkm8+Jj4+/4vFbtmxBXFzcVefbqFQquLu7X3YjIiIi2yX0tNTcuXOxbNkyrFixApmZmXjuuedQUFDQum7N/Pnz8fDDD7c+fsaMGcjPz8fcuXORmZmJFStWYPny5Zg3b56ol0BEREQWRuiiKvfeey8uXLiAN998E8XFxYiOjsbGjRvRu3dvAEBxcTEKCgpaHx8SEoKNGzfiueeew8cff4yAgAAsXrxY+Bo3REREZDmErnMjgqnWuSEiIiLT6cjnt/CrpYiIiIi6EssNERER2RSWGyIiIrIpLDdERERkU1huiIiIyKaw3BAREZFNYbkhIiIimyJ0ET8RWpb1MeXu4ERERNS1Wj6327M8n92Vm+rqagDg7uBERERWqLq6Gmq1+pqPsbsVio1GI86fPw83NzfIZLJ2PaeqqgpBQUEoLCzkqsbNeEzaxuNyJR6TtvG4XInHpG08Lk0kSUJ1dTUCAgIgl197Vo3djdzI5XIEBgZ26rncVfxKPCZt43G5Eo9J23hcrsRj0jYeF1x3xKYFJxQTERGRTWG5ISIiIpvCctMOKpUKr7/+OlQqlegoFoPHpG08LlfiMWkbj8uVeEzaxuPScXY3oZiIiIhsG0duiIiIyKaw3BAREZFNYbkhIiIim8JyQ0RERDaF5QbAkiVLEBISAicnJ8TGxmLnzp1XfeyuXbswatQo9OjRA87Ozujbty8WLlxoxrTm05Hj8nu7d++GUqnE4MGDTRtQkI4cl5SUFMhksituWVlZZkxseh39WdFqtXj55ZfRu3dvqFQq9OnTBytWrDBTWvPpyHF59NFH2/xZiYqKMmNi0+voz8rq1asxaNAguLi4wN/fH4899hguXLhgprTm09Hj8vHHH6Nfv35wdnZGZGQkVq5caaakVkKyc2vXrpUcHBykzz//XMrIyJCeeeYZqVu3blJ+fn6bjz98+LC0Zs0a6cSJE1Jubq709ddfSy4uLtKnn35q5uSm1dHj0qKyslIKDQ2VEhMTpUGDBpknrBl19Lhs375dAiBlZ2dLxcXFrTe9Xm/m5KbTmZ+VO+64Qxo+fLiUnJws5ebmSvv375d2795txtSm19HjUllZednPSGFhoeTp6Sm9/vrr5g1uQh09Jjt37pTkcrn0wQcfSGfPnpV27twpRUVFSdOmTTNzctPq6HFZsmSJ5ObmJq1du1Y6c+aM9M0330iurq7S+vXrzZzcctl9uRk2bJg0Y8aMy+7r27ev9OKLL7b7e9x5553Sgw8+2NXRhOrscbn33nulV155RXr99ddtstx09Li0lJtLly6ZIZ0YHT0mmzZtktRqtXThwgVzxBPmRt9bfvjhB0kmk0l5eXmmiCdER4/Ju+++K4WGhl523+LFi6XAwECTZRSho8clPj5emjdv3mX3PfPMM9KoUaNMltHa2PVpKZ1Oh0OHDiExMfGy+xMTE7Fnz552fY/09HTs2bMHN998sykiCtHZ4/LFF1/gzJkzeP31100dUYgb+XmJiYmBv78/EhISsH37dlPGNKvOHJP169cjLi4O77zzDnr27ImIiAjMmzcP9fX15ohsFl3x3rJ8+XJMmDABvXv3NkVEs+vMMRk5ciTOnTuHjRs3QpIklJaW4vvvv8ett95qjshm0ZnjotVq4eTkdNl9zs7OOHDgABobG02W1ZrYdbmpqKiAwWCAr6/vZff7+vqipKTkms8NDAyESqVCXFwcZs6ciSeeeMKUUc2qM8clJycHL774IlavXg2l0jb3Y+3McfH398dnn32GpKQkrFu3DpGRkUhISMCOHTvMEdnkOnNMzp49i127duHEiRP44YcfsGjRInz//feYOXOmOSKbxY28twBAcXExNm3aZPfvKyNHjsTq1atx7733wtHREX5+fvDw8MCHH35ojshm0ZnjMmnSJCxbtgyHDh2CJElIS0vDihUr0NjYiIqKCnPEtni2+SnUQTKZ7LKvJUm64r4/2rlzJ2pqarBv3z68+OKLCAsLw3333WfKmGbX3uNiMBhw//3344033kBERIS54gnTkZ+XyMhIREZGtn4dHx+PwsJCvPfeexgzZoxJc5pTR46J0WiETCbD6tWrW3f4ff/993H33Xfj448/hrOzs8nzmktn3lsA4Msvv4SHhwemTZtmomTidOSYZGRkYM6cOXjttdcwadIkFBcX4/nnn8eMGTOwfPlyc8Q1m44cl1dffRUlJSUYMWIEJEmCr68vHn30UbzzzjtQKBTmiGvx7HrkxsvLCwqF4op2XFZWdkWL/qOQkBAMGDAATz75JJ577jn84x//MGFS8+rocamurkZaWhpmzZoFpVIJpVKJN998E0ePHoVSqcS2bdvMFd2kbuTn5fdGjBiBnJycro4nRGeOib+/P3r27NlabACgX79+kCQJ586dM2lec7mRnxVJkrBixQo89NBDcHR0NGVMs+rMMVmwYAFGjRqF559/HgMHDsSkSZOwZMkSrFixAsXFxeaIbXKdOS7Ozs5YsWIF6urqkJeXh4KCAgQHB8PNzQ1eXl7miG3x7LrcODo6IjY2FsnJyZfdn5ycjJEjR7b7+0iSBK1W29XxhOnocXF3d8fx48dx5MiR1tuMGTMQGRmJI0eOYPjw4eaKblJd9fOSnp4Of3//ro4nRGeOyahRo3D+/HnU1NS03nfq1CnI5XIEBgaaNK+53MjPSmpqKk6fPo3HH3/clBHNrjPHpK6uDnL55R9TLSMTko1si3gjPysODg4IDAyEQqHA2rVrcdttt11xvOyWiFnMlqTlErzly5dLGRkZ0rPPPit169at9QqFF198UXrooYdaH//RRx9J69evl06dOiWdOnVKWrFiheTu7i69/PLLol6CSXT0uPyRrV4t1dHjsnDhQumHH36QTp06JZ04cUJ68cUXJQBSUlKSqJfQ5Tp6TKqrq6XAwEDp7rvvlk6ePCmlpqZK4eHh0hNPPCHqJZhEZ/8NPfjgg9Lw4cPNHdcsOnpMvvjiC0mpVEpLliyRzpw5I+3atUuKi4uThg0bJuolmERHj0t2drb09ddfS6dOnZL2798v3XvvvZKnp6eUm5sr6BVYHrsvN5IkSR9//LHUu3dvydHRURoyZIiUmpra+mePPPKIdPPNN7d+vXjxYikqKkpycXGR3N3dpZiYGGnJkiWSwWAQkNy0OnJc/shWy40kdey4/Oc//5H69OkjOTk5Sd27d5duuukmacOGDQJSm1ZHf1YyMzOlCRMmSM7OzlJgYKA0d+5cqa6uzsypTa+jx6WyslJydnaWPvvsMzMnNZ+OHpPFixdL/fv3l5ydnSV/f3/pgQcekM6dO2fm1KbXkeOSkZEhDR48WHJ2dpbc3d2lqVOnSllZWQJSWy6ZJNnI2B4RERER7HzODREREdkelhsiIiKyKSw3REREZFNYboiIiMimsNwQERGRTWG5ISIiIpvCckNEREQ2heWGiIiIbArLDREREdkUlhsiIiKyKSw3RGT1ysvL4efnh3//+9+t9+3fvx+Ojo7YsmWLwGREJAL3liIim7Bx40ZMmzYNe/bsQd++fRETE4Nbb70VixYtEh2NiMyM5YaIbMbMmTOxdetWDB06FEePHsXBgwfh5OQkOhYRmRnLDRHZjPr6ekRHR6OwsBBpaWkYOHCg6EhEJADn3BCRzTh79izOnz8Po9GI/Px80XGISBCO3BCRTdDpdBg2bBgGDx6Mvn374v3338fx48fh6+srOhoRmRnLDRHZhOeffx7ff/89jh49CldXV4wbNw5ubm74+eefRUcjIjPjaSkisnopKSlYtGgRvv76a7i7u0Mul+Prr7/Grl27sHTpUtHxiMjMOHJDRERENoUjN0RERGRTWG6IiIjIprDcEBERkU1huSEiIiKbwnJDRERENoXlhoiIiGwKyw0RERHZFJYbIiIisiksN0RERGRTWG6IiIjIprDcEBERkU1huSEiIiKb8v8BXkCRAU8qUO0AAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "B._posterior.plot(\"pdf\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Updating\n", + "\n", + "A significant advantage of the Bayesian approach is its ability to incorporate new data and continuously update our posterior distribution. Suppose we conduct 100 additional coin tosses and observe that, unexpectedly, 80 of these new tosses result in tails:" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0 0\n", + "1 0\n", + "2 0\n", + "3 0\n", + "4 0\n", + " ..\n", + "95 0\n", + "96 0\n", + "97 1\n", + "98 0\n", + "99 0\n", + "Name: head, Length: 100, dtype: int64" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Define the number of ones and zeros\n", + "ones = [1] * 20\n", + "zeros = [0] * 80\n", + "\n", + "# Combine them into a single list\n", + "data = ones + zeros\n", + "\n", + "# Shuffle the list to randomize the order\n", + "np.random.shuffle(data)\n", + "\n", + "# Create a Series from the shuffled list\n", + "y2 = pd.Series(data, name=\"head\")\n", + "\n", + "y2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This new data can be integrated into our existing model to further refine our belief about the coin's fairness. We do so using the `update` method." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
BayesianProportionEstimator(prior=Beta(alpha=2, beta=2), prior_alpha=2,\n",
+       "                            prior_beta=2)
Please rerun this cell to show the HTML repr or trust the notebook.
" + ], + "text/plain": [ + "BayesianProportionEstimator(prior=Beta(alpha=2, beta=2), prior_alpha=2,\n", + " prior_beta=2)" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "X = y2\n", + "B.update(X, y2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We see that after this update, the posterior becomes $\\text{Beta}(\\alpha=30, \\beta=84)$, which is a left-skewing distribution with a mean of 0.263." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
Beta(alpha=30, beta=84)
Please rerun this cell to show the HTML repr or trust the notebook.
" + ], + "text/plain": [ + "Beta(alpha=30, beta=84)" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "B._posterior" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.2631578947368421" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "B._posterior.mean()" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGwCAYAAABcnuQpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABNdElEQVR4nO3dd3hUVeI+8Hf6pEN6Qio1QICQ0FGwYFCKIqICoui6rqy4gu66K6v7sy6ouypfG5ZVRJAiiCsKKCiC9BKSUBIIAVIgvfc2c39/TJIlUkyZmXPvzPt5nnl8mEzCy3Ugb8459xyVJEkSiIiIiBRKLToAERERUVewzBAREZGiscwQERGRorHMEBERkaKxzBAREZGiscwQERGRorHMEBERkaJpRQewNbPZjJycHHh4eEClUomOQ0RERO0gSRIqKysRHBwMtfraYy8OX2ZycnIQGhoqOgYRERF1QnZ2NkJCQq75GocvMx4eHgAsF8PT01NwGiIiImqPiooKhIaGtn4fvxaHLzMtU0uenp4sM0RERArTniUiQhcA//LLL5g6dSqCg4OhUqnw3//+t83HJUnCCy+8gODgYLi4uOCGG27AyZMnxYQlIiIiWRJaZqqrqzFkyBC8++67V/z466+/jjfffBPvvvsuDh8+jMDAQNxyyy2orKy0c1IiIiKSK6HTTLfddhtuu+22K35MkiQsXboUzz77LKZPnw4AWLFiBQICArB69Wo8+uijV/y8+vp61NfXt/66oqLC+sGJiIhINmS7z8z58+eRl5eH+Pj41ucMBgPGjx+Pffv2XfXzlixZAi8vr9YH72QiIiJybLItM3l5eQCAgICANs8HBAS0fuxKFi1ahPLy8tZHdna2TXMSERGRWLK/m+nXq5glSbrmymaDwQCDwWDrWERERCQTsh2ZCQwMBIDLRmEKCgouG60hIiIi5yXbMhMZGYnAwEBs37699bmGhgbs2rULY8aMEZiMiIiI5EToNFNVVRXS09Nbf33+/HkkJSXB29sbYWFhWLhwIRYvXow+ffqgT58+WLx4MVxdXTF79myBqYmIiEhOhJaZI0eO4MYbb2z99VNPPQUAmDt3Lj777DP89a9/RW1tLR577DGUlpZi5MiR2LZtW7u2NiYiIiLnoJIkSRIdwpYqKirg5eWF8vJyHmdARESkEB35/i3bNTNERERE7SH7W7OJSP6aTGZU1TehttEEV50W7kYtNOrfPhyOiMgaWGaIqMPyK+rw86kC7E4vQmpOBTKKq2G+ZMJap1Eh3McN0cGeGNvbFzdF+cPHnfs/EZFtsMwQUbtIkoRdaYVYsS8Du9IK25SXFlq1Ck1mCY0mCekFVUgvqMJ/k3KgVatwY5Q/HhobgdE9fa658SURUUexzBDRbzpwrhj/+uE0EjJLW5+LCe2GG/v5Iy68O3r7u8PHXQ+dRo2GJjMKKutwpqAKCRml2JlWgBMXK7A9JR/bU/IxIsIbf47vi5E9fQT+iYjIkfBuJiK6qvKaRryyOQXrEy4AAIw6Ne4bGY45o8IR6evW7q+Tll+Jlfszse5wNhpMZgDAjLgQPDupP7q76W2SnYiUrSPfv1lmiOiK9p0twsK1SSiorIdKBcweEYYFN/eBv6ex018zr7wO//fTGaw5lAUA8HHT4617YzCur5+1YhORg2CZuQTLDFHHSJKET/acx5Ktp2AyS+jp54Z/zRiMuHBvq/0eCZklWLTxONLyq6BSAY/f2BsLJ/TlHVBE1Ir7zBBRpzSZzHh6wzG8sjkVJrOE6UN7YMsT11u1yABAXLg3Nj1+HWaPDIMkAe/sSMe8VQmobTBZ9fchIufAMkNEAIC6RhPmrTqKDQkXoFGr8OLtA/HGPUNg1Gls8vsZdRosvnMQlt4bA71Wje0p+Zj18QEUV9Xb5PcjIsfFMkNEqGs04eEVh/Fjaj70WjU+mBOHuWMi7HIL9bShPfDF70fCy0WHpOwyzPr4AIpYaIioA1hmiJxco8mM+V8cxd70YrjpNVjx0AjcMiDArhmGR3jjqz+OQYCnAWn5VZjNQkNEHcAyQ+TEzGYJT32ZjJ9OFcCgVePTB4djdC8x+7/09nfH2j+Mbi00c/5zEBV1jUKyEJGysMwQObE3tp/Gt8k50GlU+OD+OOEb2UX6umHtH0bDz8OAU3mVmLcyAQ1NZqGZiEj+WGaInNQ3SRfx3s9nAQCv3TUYN/bzF5zIItLXDcsfHA43vQb7zhbjrxuS4eA7SBBRF7HMEDmhxKxSPL3hGADgjzf0wvTYEMGJ2oru4YVlc+KgVavw36QcvPXjGdGRiEjGWGaInExJdQPmrbJM30zoH4Cn4/uJjnRF4/r6YfH0QQCAt386gx9T8gUnIiK5YpkhciKSJOEv65ORX1GPXn5uWDozBmoZ77p7z7BQPDgmAgDw5LoknC+qFhuIiGSJZYbIiXy6NwM7ThVAr1Xj3dmxcDdoRUf6Tc9O7o/hEd1RWd+ER1ceQU1Dk+hIRCQzLDNETuL4hXK8ujUVAPCPKQPQP0gZZ5XpNGq8NzsWfh6WW7Zf+jZFdCQikhmWGSInUNdowlNfJqHRJOG26EDMGRkmOlKH+Hsa8fbMoVCpgLWHs7H1eK7oSEQkIywzRE7g7Z/O4ExBFXzdDVh85yC7HFNgbaN7+WDe+F4AgGc2Hkduea3gREQkFywzRA7u2IUyfPjLOQDAP++MRnc3veBEnffkhL4YHOKF8tpGPLUuGWYz958hIpYZIodW32TC0+uPwWSWMHVIMCYODBQdqUv0WjWW3hsDF50G+88VY/WhLNGRiEgGWGaIHNh/dp/H6fxK+Ljp8eLtA0XHsYqefu74662WvXFe3XoKOWWcbiJydiwzRA4qu6QG7+yw7Jz7jykD4K3g6aVfe2B0BGLDuqGqvgnP/fcEjzsgcnIsM0QO6qXvUlDXaMbISG/cERMsOo5VadQqvD5jMPQaNXacKsCm5BzRkYhIIJYZIge041Q+tqfkQ6tW4eVp0Yq8e+m39Pb3wBM39wYAvLDpJEqqGwQnIiJRWGaIHEx9kwkvbLJsLPe76yLRN8BDcCLbeXR8L0QFeqC0phH/+uGU6DhEJAjLDJGDWbk/E1klNfD3MOCJm/uIjmNTOo0ar0yLBmDZTC85u0xsICISgmWGyIGU1TTg7Z8si37/HN9XEWcvddWwCG9MH9oDkgT8v29OcO8ZIifEMkPkQN7dkY6Kuib0C/DAjLhQ0XHs5plJUXA3aJF8oRxfHskWHYeI7IxlhshBZBXXYMX+DADAoklR0Kgdb9Hv1fh7GLFwgmVK7bXvT6GshouBiZwJywyRg3j9h1NoNEm4vo8vxvf1Ex3H7uaOiUDfAHeU1jTinR3pouMQkR2xzBA5gOMXyvHdsVyoVMCi2/o75K3Yv0WnUeO5yQMAAJ/vz0BmcbXgRERkLywzRA7gze2nAQDTYnpgQLCn4DTijOvrh+v7+KLRJOH1H06LjkNEdsIyQ6RwR7NK8fPpQmjUKoe/Fbs9/j6pP1QqYPOxXCRkloqOQ0R2wDJDpHBvbU8DAEwf2gORvm6C04jXP8gTM2JDAACLt6Ty3CYiJ8AyQ6Rgh86XYPeZImg5KtPGn+P7wahTIyGzFD+czBMdh4hsjGWGSMFaRmXuHhaKUG9XwWnkI9DLiEeu7wkA+Pe2NJi4kR6RQ2OZIVKo/WeLsf9cMfQaNR6/qbfoOLLz++t7wstFh/SCKmxKvig6DhHZEMsMkUK9v9Oyl8rdw0LQo5uL4DTy4+Wiwx/GWUZnlv54Bo0ms+BERGQrLDNECnT8Qjl2nymCRq3CvPG9RMeRrQfHRMDXXY/M4hp8lXBBdBwishGWGSIFWrbLMiozdXAQ18pcg5tBiz/eYJmCe/unM6hvMglORES2wDJDpDDnCquw9YTlDp15N3BU5rfcNzIMgZ5G5JTXYc3BLNFxiMgGWGaIFObDXecgScDNUf6ICnTe3X7by6jT4E83W0Zn3v35LGobODpD5GhYZogUJLe8FhsTLWs/HruRozLtdXdcKEK6u6Coqh5rD3N0hsjRsMwQKch/dp9Ho0nCiEhvxIV7i46jGHqtGn9snpL7cNc5rp0hcjAsM0QKUV7biDWHLKMKf+RamQ6bEReCAE8D8irqsPEo950hciQsM0QKse5wFmoaTOgb4I4b+vqJjqM4Bq0Gj46zlMD3d6ajifvOEDkMlhkiBWgymbFiXyYA4OHrIqFSqQQnUqZZI8Lg46ZHdkktNiXniI5DRFbCMkOkAD+czMfFslp4u+lxR0wP0XEUy0Wvwe+bz2x67+d0mHlmE5FDYJkhUoBP954HAMwZGQajTiM4jbLNGRUGT6MWZwur8T1P1CZyCCwzRDKXlF2GhMxS6DQqzBkVLjqO4nkYdXhobCQA4J0d6ZAkjs4QKR3LDJHMfbrHMiozdUgw/D2NgtM4hofGRsBNr0FqbgV2phWKjkNEXcQyQyRjueW12HI8FwDwu+bRBOq6bq56zB4ZBgD4aNc5wWmIqKtYZohkbOX+TDSZJYyM9EZ0Dy/RcRzKQ2MjoVWrsP9cMY5fKBcdh4i6gGWGSKbqm0xYdzgbgGVahKwruJsLpg4JBgB8+MtZwWmIqCtkXWaamprw3HPPITIyEi4uLujZsydeeuklmM3c7Ioc3/cn8lBc3YAATwMm9A8QHcchPdJ8m/aW47nILqkRnIaIOkvWZea1117DBx98gHfffRepqal4/fXX8a9//QvvvPOO6GhENrfqgGWTvNkjwqHVyPqvqmINCPbE9X18YZaAT5oXWhOR8sj6X8j9+/fjjjvuwOTJkxEREYEZM2YgPj4eR44cER2NyKZO5VXgcEYpNGoVZo4IFR3Hof1hnGV0Zt3hbJRWNwhOQ0SdIesyc9111+Gnn35CWloaACA5ORl79uzBpEmTrvo59fX1qKioaPMgUpqWUZmJAwMQwNuxbeq63r4YEOSJ2kZT63UnImWRdZn529/+hlmzZiEqKgo6nQ5Dhw7FwoULMWvWrKt+zpIlS+Dl5dX6CA3lT7WkLFX1Tfi6+VRnbpJneyqVqnV0ZsX+DNQ1mgQnIqKOknWZWbduHVatWoXVq1fj6NGjWLFiBf79739jxYoVV/2cRYsWoby8vPWRnZ1tx8REXfd14kVUN5jQy88No3v6iI7jFCYPDkKwlxFFVQ3YlMQDKImURtZl5umnn8YzzzyDmTNnYtCgQbj//vvx5JNPYsmSJVf9HIPBAE9PzzYPIqWQJAmr9lumOuaMCufp2Hai06jxwJgIAJZzsHjEAZGyyLrM1NTUQK1uG1Gj0fDWbHJYhzNKcTq/Ei46DabHhoiO41RmDg+FUafGqbxKHDhXIjoOEXWArMvM1KlT8c9//hObN29GRkYGvv76a7z55pu48847RUcjsonVBy2jMnfEBMPLRSc4jXPp5qpvLZDL9/I2bSIlkXWZeeeddzBjxgw89thj6N+/P/7yl7/g0Ucfxcsvvyw6GpHVldc0YsuJPADArBFhgtM4p4eap5q2p+ZzEz0iBZF1mfHw8MDSpUuRmZmJ2tpanD17Fq+88gr0er3oaERW99+ki2hoMiMq0AODQ3gOkwh9AjxwfR9fSBKwYl+G6DhE1E6yLjNEzkKSJKxtPodp5vBQLvwVqOUcrHVHslFd3yQ2DBG1C8sMkQycuFiB1NwK6LVqTBvaQ3Qcp3ZDX39E+Liisq4JG49eEB2HiNqBZYZIBtYdyQIATBwYiG6unEYVSa1WYW7z2pnl+zJgNvM2bSK5Y5khEqy2wYRvEi0btc0czh2r5WBGXAjcDVqcK6zGL2cKRcchot/AMkMk2NYTuaisb0Kotwt3/JUJD6MOdw9ruU07Q2wYIvpNLDNEgrUs/L0nLhRqNRf+ysWDYyKgUgG70gpxtrBKdBwiugaWGSKBzhVW4dD5EqhVwIxh3PFXTsJ93HBTP38AwBcHsgSnIaJrYZkhEujLI5a7Zcb39UOQl4vgNPRrc0ZbTi1fn5CNmgbepk0kVywzRII0mcz4qvnW33u58FeWxvfxQ6i3CyrrmvBtMk/TJpIrlhkiQX45U4jCynr4uOlxU1SA6Dh0BWq1CnNGWkZnPt+fydO0iWSKZYZIkK8SLgIA7ojpAb2WfxXl6u5hodBr1TiZU4Gk7DLRcYjoCvgvKJEA5TWN2J6aDwCYHssdf+XM202PKYODAAArD2QKTkNEV8IyQyTAd8dzWg+VHBjsKToO/YY5oyxTTd8dy0VpdYPgNET0aywzRAJ8lWBZ+Ds9tgcPlVSAoaHdMDDYEw1NZqxPyBYdh4h+hWWGyM7OF1XjaFYZ1CpgWgynmJRApVLh/ubRmVUHsnheE5HMsMwQ2VnLSczX9/GDv6dRcBpqr9tjguFh1CKrpIbnNRHJDMsMkR2ZzRI2HrXcxXRXHHf8VRJXvRYzmv+freJCYCJZYZkhsqOD50twsawWHgYt4gdwbxmlaVkI/NOpAmSX1AhOQ0QtWGaI7Khlimny4CAYdRrBaaijevm5Y2xvH0gSsOYQz2sikguWGSI7qWlowpbjuQCA6bGcYlKqlh2BvzxyAY0ms+A0RASwzBDZzbaT+ahuMCHU2wXDI7qLjkOdNGFAAHzdDSiqqsdPzRsfEpFYLDNEdtJyqOT0oSHcW0bBdBo17h5mGVlbfYh7zhDJAcsMkR3kV9RhT3oRAOAuTjEp3szmU853nynkQmAiGWCZIbKDb5NzIElAXHh3hPm4io5DXRTu44brevtCkoB1hzk6QyQaywyRHXyTlAMAmBYTLDgJWcusEWEAgC+PZHMhMJFgLDNENna2sArHL5ZDo1Zh0qAg0XHISm4ZEAAfNz0KKuux41SB6DhETo1lhsjGNjWPylzfxxc+7gbBacha9Fo1ZjQvBOaeM0RiscwQ2ZAkSdiUbCkzd3CKyeHMHG6ZatqVVogLpVwITCQKywyRDR27UI7zRdUw6tSIHxAoOg5ZWaSvG8b0suwI/CUXAhMJwzJDZEMtC39vGRAIN4NWcBqyhZaFwOuOZKOJC4GJhGCZIbIRk1nCt8eap5iGcIrJUcUPDIC3mx75FfX4+XSh6DhETollhshGDpwrRmFlPbq56jCur5/oOGQjBq0GM+K4EJhIJJYZIhv5JukiAGDSoCDotfyr5shadgTeeboAOWW1gtMQOR/+C0tkA3WNJmw9ngeAU0zOoKefO0b39IFZsmyiR0T2xTJDZAM7Txegsr4JQV5GDI/wFh2H7GDWyOaFwIezYTJLgtMQOReWGSIbaLmL6fYhwVCreUK2M5g4MADdXXXILa/DztPcEZjInlhmiKysoq4RPzVvb39HTA/BacheDFpN64noPHySyL5YZois7IcTeWhoMqOPvzv6B3mIjkN2dG/zQuCfThWgoLJOcBoi58EyQ2RlLVNMd8QEQ6XiFJMz6RPggaFh3WAyS9h49KLoOEROg2WGyIoKKuqw72wRAOD2IZxickYtt2l/eTgbksSFwET2wDJDZEWbj+fCLAFDw7ohzMdVdBwSYPLgYLjqNThXVI0jmaWi4xA5BZYZIiv67lguAMtdTOSc3A1aTBkcBABYe4gLgYnsgWWGyEpyymqRkFkKlcqy6y85r3uHW/ac2XI8F5V1jYLTEDk+lhkiK9ly3DIqMzzCGwGeRsFpSKTYsG7o7e+O2kYTvk3OFR2HyOGxzBBZybfNU0wtUwzkvFQqFe4dZlkIvO4wD58ksjWWGSIryC6pQXJ2GdQq4NboQNFxSAbujO0BrVqF5AvlOJVXIToOkUNjmSGygs3NU0wjI33g78EpJgJ83Q24ZUAAAO4ITGRrLDNEVvDdMctGeVOGcIqJ/uee5j1nvk68iPomk+A0RI6LZYaoizKKqnHiYgU0ahVuHcgpJvqfcX38EORlRFlNI7adzBcdh8hhscwQdVHLFNOYXj7wcTcITkNyolGrcHec5fDJL49wqonIVlhmiLroO97FRNdwd/NdTXvSi5BdUiM4DZFjYpkh6oKzhVVIza2AVq3CRE4x0RWEertibG8fSBKwPuGC6DhEDollhqgLvmveEO26Pr7o5qoXnIbk6p7m0ZkNR7JhMvPwSSJrY5kh6oLNxy13MU3m8QV0DRMHBsLLRYec8jrsSS8SHYfI4bDMEHVSWn4l0vKroNeoEc8pJroGo06DO4f2AAB8yT1niKyOZYaok1oW/o7r6wsvF53gNCR3LVNN21LyUFxVLzgNkWNhmSHqBEmSWjfKm8y7mKgdBgR7YnCIFxpNEr5OvCg6DpFDYZkh6oTU3EqcK6yGXqvGhP4BouOQQrSMznx5JBuSxIXARNbCMkPUCS0Lf2/o6wcPI6eYqH1ujwmGUadGWn4VErPLRMchchiyLzMXL17EnDlz4OPjA1dXV8TExCAhIUF0LHJilimm5o3yhgQLTkNK4mnUYVK0ZVqSC4GJrEfWZaa0tBRjx46FTqfD1q1bkZKSgjfeeAPdunUTHY2c2MmcCmQW18CoU+PmKH/RcUhh7m0+fPLb5BxU1zcJTkPkGLSiA1zLa6+9htDQUCxfvrz1uYiICHGBiAB827zw96Yof7gZZP1XiGRoRKQ3In3dcL6oGpuP5baerE1EnSfrkZlNmzZh2LBhuPvuu+Hv74+hQ4fi448/vubn1NfXo6Kios2DyFokScLm1rOYOMVEHadSqXD3MMvhk+t4+CSRVci6zJw7dw7Lli1Dnz598MMPP2DevHl44okn8Pnnn1/1c5YsWQIvL6/WR2gof+oh60m+UI4LpbVw1WtwYz9OMVHnzIgNgUatQkJmKdILKkXHIVI8WZcZs9mM2NhYLF68GEOHDsWjjz6KRx55BMuWLbvq5yxatAjl5eWtj+xs/uRD1vNdsmWK6eb+AXDRawSnIaXy9zTipub1Vuu4EJioy2RdZoKCgjBgwIA2z/Xv3x9ZWVlX/RyDwQBPT882DyJrMJslbDlumWLiWUzUVfc27znz1dGLaGgyC05DpGyyLjNjx47F6dOn2zyXlpaG8PBwQYnImSVmlyKnvA5ueg1u6OcnOg4p3A39/ODvYUBJdQN+Ss0XHYdI0WRdZp588kkcOHAAixcvRnp6OlavXo2PPvoI8+fPFx2NnNC3yZZRmVsGBMCo4xQTdY1Wo8aMOMtC4LWcaiLqElmXmeHDh+Prr7/GmjVrEB0djZdffhlLly7FfffdJzoaOZlLp5h4FxNZS8vxBr+cKcTFslrBaYiUS/abZEyZMgVTpkwRHYOc3OGMEhRU1sPDqMX1fX1FxyEHEeHrhtE9fbD/XDE2HLmABRP6iI5EpEiyHpkhkovNzaMy8QMCYdByiomsp2VH4C+PZMNk5uGTRJ3BMkP0G0xmCVuO5wEApgzhXUxkXbdGB8LTqMXFslrsTS8SHYdIkVhmiH7DwXPFKKqqh5eLDmN7cYqJrMuo0+DOoT0AcEdgos5imSH6Dd81TzHdOjAQei3/ypD1tZzPtO1kHkqqGwSnIVIe/stMdA1NJjO+P8EpJrKtgcFeGNTDC40mCRuPXhAdh0hxWGaIrmH/uWKUVDfA202P0T19RMchB3bpQmBJ4kJgoo5gmSG6hu+aN8q7NToQWg3/upDt3B4TDKNOjbT8KiRml4mOQ6Qo/NeZ6Coamsz4/mTzFBPPYiIb8zTqMKn5fbbuEBcCE3UEywzRVew9W4Ty2kb4uusxklNMZAczh4cBAL49loOq+ibBaYiUg2WG6Cpapphuiw6CRq0SnIacwfCI7ujp64aaBhM2H8sRHYdIMVhmiK6gvsmEbSnNU0yDOcVE9qFSqVpv0+bhk0TtxzJDdAW704pQWdcEfw8Dhkd4i45DTmR6bA9o1SokZpUhLb9SdBwiRWCZIbqClrOYJg0KgppTTGRH/h5G3BTlDwBYx9EZonZhmSH6lbpGE7an5APgFBOJMXOEZapp49ELqG8yCU5DJH8sM0S/svN0IarqmxDkZURsWHfRccgJjevjh0BPI0prGvFjSoHoOESyxzJD9CvfNd9FMplTTCSIVqPGjLgQAMDaw1mC0xDJn7ajnyBJEnbt2oXdu3cjIyMDNTU18PPzw9ChQzFhwgSEhobaIieRXdQ2mPBTquUn4SlDggWnIWd2z7BQvPtzOvakF+FCaQ1CuruKjkQkW+0emamtrcXixYsRGhqK2267DZs3b0ZZWRk0Gg3S09Px/PPPIzIyEpMmTcKBAwdsmZnIZnacKkBtowkh3V0wJMRLdBxyYmE+rhjb2weSBKw/wsMnia6l3SMzffv2xciRI/HBBx9g4sSJ0Ol0l70mMzMTq1evxr333ovnnnsOjzzyiFXDEtla6xTT4CCoVJxiIrHuGRaKvenFWH8kG0/c3IebNxJdRbvLzNatWxEdHX3N14SHh2PRokX485//jMzMzC6HI7Kn6vom7DhlmWKaOphTTCTexIGB8HLRIae8DnvSizC+r5/oSESy1O5ppt8qMpfS6/Xo06dPpwIRifJjaj7qm8wI93HFwGBP0XGIYNRpcOfQHgCAdVwITHRVnbqb6R//+AdMpsv3PigvL8esWbO6HIpIhO+OWTbKm8IpJpKRe5uPN9ieko/iqnrBaYjkqVNl5vPPP8fYsWNx9uzZ1ud27tyJQYMGISMjw1rZiOymsq4Ru04XAgCmcIqJZKR/kCeGhHih0STh68SLouMQyVKnysyxY8cQERGBmJgYfPzxx3j66acRHx+PBx98EHv27LF2RiKb256SjwaTGb383BAV6CE6DlEblx4+KUmS4DRE8tPhfWYAwMvLC2vXrsWzzz6LRx99FFqtFlu3bsXNN99s7XxEdtEyxTR5cDCnmEh2bh8SjFe+S0V6QRWOZpUiLpyHnxJdqtM7AL/zzjt46623MGvWLPTs2RNPPPEEkpOTrZmNyC7Kaxqx+4xlimkqz2IiGfIw6jC5+b25+iAPnyT6tU6Vmdtuuw0vvvgiPv/8c3zxxRdITEzEuHHjMGrUKLz++uvWzkhkUz+k5KHRJKFfgAf6BHCKieRp1ogwAJa9kMprGgWnIZKXTpWZpqYmHDt2DDNmzAAAuLi4YNmyZdiwYQPeeustqwYksrX/TTFxVIbkKzasG6ICPVDfZMbGRO4ITHSpTpWZ7du3Izj48js+Jk+ejOPHj3c5FJG9lFQ3YG96EQDLLdlEcqVSqTB7pGV0ZvXBLC4EJrqE1U/N9vX1BQD+RSNF+OFkHkxmCQOCPNHTz110HKJrmja0B1x0GpwpqMKRzFLRcYhko91lpn///li9ejUaGhqu+bozZ87gj3/8I1577bUuhyOytZazmKYM4agMyZ+nUYepQ1oWAnNHYKIW7b41+7333sPf/vY3zJ8/H/Hx8Rg2bBiCg4NhNBpRWlqKlJQU7NmzBykpKXj88cfx2GOP2TI3UZcVVtZj/9liAMCUQdwoj5Rh9shwfHnkAjYfz8X/mzIA3d30oiMRCdfuMnPTTTfh8OHD2LdvH9atW4fVq1cjIyMDtbW18PX1xdChQ/HAAw9gzpw56Natmw0jE1nH9yfzYJaAwSFeCPNxFR2HqF2GhHhhQJAnUnIr8NXRC/j99T1FRyISrsOb5o0ZMwZjxoyxRRYiu/ouuXmKiQt/SUFaFgI/998TWH0oCw9fF8mNHsnptXvNjLe3N4qKLHd9/O53v0NlZaXNQhHZWn5FHQ5llAAAJg1imSFluSMmGK56Dc4VVuPg+RLRcYiEa3eZaWhoQEVFBQBgxYoVqKurs1koIlvbejwXkgQMDeuGkO6cYiJl8TDqcEeMZZ0XFwITdWCaafTo0Zg2bRri4uIgSRKeeOIJuLi4XPG1n376qdUCEtlCy0Z5PCGblGr2iHCsOZSN70/koaS6Ad5cCExOrN0jM6tWrcKkSZNQVVUFlUqF8vJylJaWXvFBJGc5ZbWte3RMGhQoOA1R5wwK8cKgHl5oMJmxIYHnNZFza/fITEBAAF599VUAQGRkJFauXAkfHx+bBSOylS3HLaMywyO6I8jryqOLREowe2QYFm08jjWHsvHI9T25EJicVqd2AD5//jyLDCnWpua7mKYO4RQTKdvtQ4LhbtDifFF1655JRM6o3SMzb7/9dru/6BNPPNGpMES2dr6oGsculEOjVvEuJlI8N4MW04YGY9WBLHxxKAtjevuKjkQkRLvLzK9Pwy4sLERNTU3rBnllZWVwdXWFv78/ywzJ1qYky6jM2N6+8HU3CE5D1HWzR4Rj1YEsbDuZh6Kqer6vySm1e5rp/PnzrY9//vOfiImJQWpqKkpKSlBSUoLU1FTExsbi5ZdftmVeok6TJAmbki8CsAzPEzmCAcGeiAnthkaThPVHLoiOQyREp9bM/OMf/8A777yDfv36tT7Xr18/vPXWW3juueesFo7ImlJyK3C2sBp6rRoTBwaIjkNkNbNHhgEAvjiYCZNZEpyGyP46VWZyc3PR2Nh42fMmkwn5+fldDkVkCy0Lf2+O8oeHUSc4DZH13D4kGN1cdbhQWoudpwtExyGyu06VmZtvvhmPPPIIjhw5Akmy/BRw5MgRPProo5gwYYJVAxJZg9ks4dvm9TKcYiJHY9RpcM+wUADA5/szBachsr9OlZlPP/0UPXr0wIgRI2A0GmEwGDBixAgEBQXhP//5j7UzEnVZQlYpcsrr4G7Q4sYof9FxiKxuzshwqFTArrRCZBRVi45DZFcdPjUbAPz8/LBlyxacOXMGqampaGpqQnR0NPr27WvtfERW0XIXU/zAABh1GsFpiKwvzMcVN/T1w8+nC7HqQCaemzJAdCQiu+nUyAwAfPLJJ7jzzjtx9913Y9asWZg+fTpHZUiWGk3m1l1/74jpITgNke08MDoCAPDlkWzUNpjEhiGyo07fzbRgwQJMnToV69evx/r16zF16lQ8+eSTvJuJZGdvehGKqxvg46bH2F7cuZoc1/i+fgjzdkVFXVPrNgREzqBT00zLli3Dxx9/jFmzZrU+d/vtt2Pw4MH405/+hFdeecVqAYm6quUupkmDgqDVdHowkkj21GoV5owKw+Itp/D5/kzcMyyU5zWRU+jUv+wmkwnDhg277Pm4uDg0NTV1ORSRtdQ1mrDtpGW7gDtieBcTOb57hoXCoFXjZE4FjmaViY5DZBedKjNz5szBsmXLLnv+o48+wn333dflUETW8vOpAlTVN6FHNxfEhnUXHYfI5rq56lu3H1i5P0NsGCI76dQ0E2BZALxt2zaMGjUKAHDgwAFkZ2fjgQcewFNPPdX6ujfffLPrKYk6qWWKacqQIKjVHG4n5/DA6AisT7iALcfz8Ozkevh58LwmcmydKjMnTpxAbGwsAODs2bMALLdr+/n54cSJE62v41wtiVRR14ifTll2Q+VGeeRMBoV4ISa0G5Kyy7DucBYev6mP6EhENtWpMvPzzz9bOweR1W07mY+GJjN6+7tjQJCn6DhEdvXA6HAkZZfhi4NZmDe+Fxe/k0Pju5scVssU0+1DgjlKSE5n0qAgeLvpkVtehx9TeV4TOTaWGXJIBRV12HOmEACnmMg5GXUa3Dvccl7TygMZYsMQ2RjLDDmkTck5MEtAbFg3RPi6iY5DJMR9I8OgVgF704uRll8pOg6RzSiqzCxZsgQqlQoLFy4UHYVkbuNRy+6nd8aGCE5CJE5Id1dMHBgIAFi+N0NsGCIbUkyZOXz4MD766CMMHjxYdBSSudN5lUjJrYBOo8KUQUGi4xAJ9dDYSADAxqMXUFrdIDgNkW0oosxUVVXhvvvuw8cff4zu3bnxGV3b14mWUZkb+vmju5tecBoisYZHdEd0D0/UN5mx+lCW6DhENqGIMjN//nxMnjwZEyZM+M3X1tfXo6Kios2DnIfZLOGbJEuZmT6UJ2QTqVQq/K55dGbl/kw0msyCExFZn+zLzNq1a3H06FEsWbKkXa9fsmQJvLy8Wh+hoaE2TkhycuB8MXLL6+Bp1OLGKH/RcYhkYfLgIPi6G5BXUYetJ/JExyGyOlmXmezsbCxYsACrVq2C0Whs1+csWrQI5eXlrY/s7GwbpyQ5+bp54e/kwUEw6jSC0xDJg0Grwf2jwgEAn+45LzgNkfXJuswkJCSgoKAAcXFx0Gq10Gq12LVrF95++21otVqYTKbLPsdgMMDT07PNg5xDbYOp9afOO4fyLiaiS903Kgx6jRpJ2WU4mlUqOg6RVcm6zNx88804fvw4kpKSWh/Dhg3Dfffdh6SkJGg0/Mmb/md7aj6q6psQ0t0Fw8K5UJzoUr7uBtweY9lAkrdpk6ORdZnx8PBAdHR0m4ebmxt8fHwQHR0tOh7JzH+b72K6c2gPnpBNdAUPjY0AAGw5novc8lqxYYisSNZlhqi9iqrqsSvNcnzBNN7FRHRFA4O9MKqnN0xmCZ/vzxQdh8hqOnVqtkg7d+4UHYFk6NvkHJjMEoaEeKGXn7voOESy9dDYSBw4V4LVB7PwxE194KLndD0pH0dmyCFcOsVERFc3oX8AQr1dUF7biK+OXhAdh8gqWGZI8dILqpB8oRwatQpTeEI20TVp1P/bRO8/u8/BZJYEJyLqOpYZUrwNCZafLm/o6wdfd4PgNETyd8+wUHi56JBRXIPtKdxEj5SPZYYUrclkxsbmofK7h3FvGaL2cDNoWzfR+2DXOUgSR2dI2VhmSNF2nylCQWU9vN30uCkqQHQcIsWYOyYCeq1lE73DGdxEj5SNZYYUbX2C5biKO2KCodfy7UzUXn4eBtwVa1kw/9EvZwWnIeoa/utPilVa3YAfUwoAAHfH8UBRoo76/fU9oVIBP6YWIL2gUnQcok5jmSHF2pScgwaTGQODPTEgmGdwEXVULz93TOhvmZ79+BceQEnKxTJDitUyxTQjjgt/iTrr0XE9AQBfJ15EQUWd4DREncMyQ4qUmluBExcroNOocEcMN8oj6qxhEd6IDeuGBpMZn+3LEB2HqFNYZkiR1h+x3I49oX8AvN30gtMQKdsfxvUCAKw6kImq+ibBaYg6jmWGFKehyYz/JlmOL+DeMkRdd8uAAET6uqGirgnrDmeLjkPUYSwzpDg/ny5ASXUD/DwMGNfHT3QcIsXTqFV45HrL2pn/7D6Hhiaz4EREHcMyQ4rTMsU0fWgPaDV8CxNZw/TYHvD3MCC3vA5fJ/IASlIWficgRSmsrMfO0817y3CKichqjDpN6+jMsp1n0WTi6AwpB8sMKcqGhAtoMksYGtYNvf09RMchciizR4ahm6vlAMrNx3NFxyFqN5YZUgyzWcLaw1kAgFkjwgSnIXI8bgYtfjc2EgDw/s9nYTbzAEpSBpYZUowD54qRWVwDD4MWUwYHiY5D5JDmjo6Au0GL0/mV+DE1X3QconZhmSHFWH3IMipzx9BguOq1gtMQOSYvVx3uHx0OAHjv53RIEkdnSP5YZkgRiqvqse2k5afEmcM5xURkSw9fFwmjTo3kC+XYk14kOg7Rb2KZIUXYePQiGkxmDA7xQnQPL9FxiByar7uh9YeGd3ekC05D9NtYZkj2JEnCmuaFvxyVIbKPR8f3hE6jwsHzJTiSUSI6DtE1scyQ7B06X4JzhdVw1Wtwe0yw6DhETiHIywV3xVr2cnqHozMkcywzJHtrmhf+3j4kGO4GLvwlspfHbugNjVqFXWmFSMgsFR2H6KpYZkjWymoasOVEHgDuLUNkb2E+rpjRPDqz9Mc0wWmIro5lhmRt49GLaGgyo3+QJwaHcOEvkb09flNvaNUq7D5ThEPnuXaG5IllhmRLkiR8cTATADB7RChUKpXgRETOJ9TbFfcMDwUAvLWdozMkTywzJFv7zhbjbGE13PQaTBvaQ3QcIqc1/8be0GvU2H+uGPvOct8Zkh+WGZKtz/dnAACmx4bAw6gTG4bIifXo5oKZIyyjM0u3n+GuwCQ7LDMkSzlltdieYtnxt2VrdSIS57EbekOvVeNQRgn2pheLjkPUBssMydLqg1kwS8Donj7oG+AhOg6R0wv0MuK+kZY7Ct/cfpqjMyQrLDMkO/VNpta9ZR7gqAyRbPzxhl4w6tQ4mlWGXWmFouMQtWKZIdnZejwPxdUNCPQ04pYBAaLjEFEzfw8j7h9l+QHj39tOw2zm6AzJA8sMyU7Lwt/ZI8Og1fAtSiQn88b3gptegxMXK7D5eK7oOEQAWGZIZk5cLMfRrDLoNKrWuyeISD583A34w7heAIA3tp1Go8ksOBERywzJzMr9lk3ybo0Ogr+HUXAaIrqS318fCV93PTKKa7D2cLboOEQsMyQfZTUN+Cb5IgAu/CWSMzeDFn+6qQ8A4P9+PIPq+ibBicjZscyQbKw+lIW6Rss5TMPCu4uOQ0TXMGtEGMK8XVFUVY9P95wXHYecHMsMyUKjyYzP91mmmB6+LpLnMBHJnF6rxp/j+wIAPvzlHEqqGwQnImfGMkOysOV4LvIq6uDrbsDUIUGi4xBRO0wdHIyBwZ6oqm/CuzvSRcchJ8YyQ8JJkoRPmoep544Oh0GrEZyIiNpDrVbhb7dGAQBWHchEdkmN4ETkrFhmSLjDGaU4dqEcBq0a943iwl8iJbm+jy/G9vZBg8mMV78/JToOOSmWGRLukz3nAADTY3vA200vOA0RdYRKpcJzkwdArQI2H8vF4YwS0ZHICbHMkFBZxTXY1nw69u/GRgpOQ0Sd0T/IE/cOtxxC+dK3KTzmgOyOZYaEWr7vPCQJGN/XD314OjaRYj11S1+4G7Q4frEcXydeFB2HnAzLDAlTUdeIL5t3D334Oo7KECmZn4cBj9/UGwDw+g+nuJEe2RXLDAmz6kAmqhtM6Bvgjuv7+IqOQ0Rd9NDYCIR6uyC/oh4f7jorOg45EZYZEqKu0YRP92QAsJzCy03yiJTPoNXg77f1B2DZSO9iWa3gROQsWGZIiA0JF1BUVY8e3VwwdUiw6DhEZCW3RgdiRKQ36pvMeHUrb9Um+2CZIbtrMpnx0S+W27EfuT4SOg3fhkSOQqVS4f9NGQCVCvg2OQf70otERyInwO8iZHdbTuQhq6QG3m761ts5ichxRPfwwpyRlg0w/9+mk2hoMgtORI6OZYbsSpIkLNtpWRj44JgIuOh5dAGRI/pLfD/4uOmRXlCFT/fyVG2yLZYZsqudaYVIza2Am16DuaMjRMchIhvxctVh0STLYuD/+/EMcrgYmGyIZYbsqmVUZvbIMHi56gSnISJbuiu2B4ZHdEdtowkvf5ciOg45MJYZspsjGSU4dL4EOo0KD1/XU3QcIrIxlUqFl6dFQ6NWYeuJPOxKKxQdiRwUywzZzf/9dAYAcFdsCAK9jILTEJE9RAV64sExEQCA5785gbpGk9hA5JBYZsgujmSUYPeZImjVKsy/sbfoOERkRwsn9IG/hwEZxTV47+d00XHIAbHMkF20jMrMiAtBqLer4DREZE8eRh1evH0gAMu6udTcCsGJyNGwzJDNcVSGiG6NDkT8gAA0mSU889UxmMyS6EjkQGRdZpYsWYLhw4fDw8MD/v7+mDZtGk6fPi06FnUQR2WIqGUxsIdRi+QL5VjOvWfIimRdZnbt2oX58+fjwIED2L59O5qamhAfH4/q6mrR0aidOCpDRC0CPI34e/PeM29sS0N2SY3gROQotKIDXMv333/f5tfLly+Hv78/EhISMG7cOEGpqCOW/shRGSL6n5nDQ/FN0kUcOFeCv399HJ//bgRUKpXoWKRwsh6Z+bXy8nIAgLe391VfU19fj4qKijYPEmNfehH2pBdBp+GoDBFZqFQqLJk+GAatGrvPFOHLI9miI5EDUEyZkSQJTz31FK677jpER0df9XVLliyBl5dX6yM0NNSOKamFJEl47QfL+qbZI8I4KkNErSJ93fDn+L4AgJe+TeF0E3WZYsrM448/jmPHjmHNmjXXfN2iRYtQXl7e+sjOZusX4YeT+UjOLoOrXoPHb+ojOg4RyczD1/XE8IjuqG4w4S/rk2Hm3U3UBYooM3/605+wadMm/PzzzwgJCbnmaw0GAzw9Pds8yL6aTGb8e5tlVObh6yLh52EQnIiI5EajVuGNu2Pgqtfg4PkSLN+XIToSKZisy4wkSXj88cexceNG7NixA5GRkaIjUTtsTLyI9IIqdHPV4ZFxPIOJiK4szMcVz0623N302venkF5QKTgRKZWsy8z8+fOxatUqrF69Gh4eHsjLy0NeXh5qa3mUvFzVNZqwdHsaAGD+Db3haeTJ2ER0dbNHhGF8Xz80NJnx1JfJaDSZRUciBZJ1mVm2bBnKy8txww03ICgoqPWxbt060dHoKlbuz0ROeR2CvIy4f3S46DhEJHMqlQqv3TUYXi46HLtQjv9r3s6BqCNkXWYkSbri48EHHxQdja6guKoeb++w/EP05IS+MOo0ghMRkRIEehnxzzstd6m+tzMd+9KLBCcipZF1mSFlWfrjGVTWNWFgsCfuirv2Qm0ioktNGRyMmcNDIUnAwnVJKK6qFx2JFIRlhqwiLb8SXxzMBAD8Y8oAaNTc0ZOIOub5qQPR298dBZX1eHrDMUgSb9em9mGZoS6TJAkvf5cCswTcOjAQo3r6iI5ERArkotfg3dlDodeqseNUAT7dmyE6EikEywx12c7Thdh9pgh6jRqLJkWJjkNEChYV6Il/NN+u/erWVBy7UCY2ECkCywx1SaPJjFc2pwAAHhobgXAfN8GJiEjp5owKx8SBAWg0SfjjqqMorW4QHYlkjmWGuuTTPedxtrAaPm56zL+Jh0kSUdepVCq8PmMIInxccbGsFk+sTYSJxx3QNbDMUKfllNViafOeEIsm9ecGeURkNV4uOnxwfxyMOsvp2kt/TBMdiWSMZYY67cVvT6K20YQREd64K7aH6DhE5GCiAj3x6vTBAIB3dqRje0q+4EQkVywz1Ck7TuXjh5P50KhVeHlaNFQq3opNRNY3bWgPPDgmAgDw1LoknC+qFhuIZIllhjqsrtGE5zedBGA5FbtfoIfgRETkyP4+qT/iwrujsr4JD684jPKaRtGRSGZYZqjD3vs5HdkltQjyMmLBzX1ExyEiB6fXqrHsvlgEeRlxrrAa81cf5YGU1AbLDHXIyZxyLNt5FgDw/NQBcDNoBSciImfg72nExw8Mg4tOgz3pRXjx25PcIZhascxQuzWazPjrhmNoMku4LToQt0YHiY5ERE4kuocXls6MgUoFrDqQhRX7MkRHIplgmaF2+3DXWZzMqUA3Vx1evGOg6DhE5IQmDgzEXydadhp/6bsU/Hy6QHAikgOWGWqXtPxKvP1TOgDL9JK/h1FwIiJyVvPG98SMuBCYJeCxVUeRlF0mOhIJxjJDv6nJZMbTG46hwWTGzVH+mBbDPWWISByVSoXFdw7C9X18Udtowu8+O4yzhVWiY5FALDP0mz7YdRbJ2WXwMGrxzzsHcU8ZIhJOr1Vj2Zw4DOrhhZLqBjzwySHkV9SJjkWCsMzQNSVmleKt5iMLnp86EIFenF4iInlwN2ix/KHhrWc4zf30ECrquAeNM2KZoauqqm/CgrVJMJklTBkcxCMLiEh2fN0NWPnwSPh5GHAqrxIPLT+M6vom0bHIzlhm6Kpe2HQSWSU16NHNhdNLRCRbod6uWPHQCHgatUjILMVDnx1GTQMLjTNhmaEr+u5YDjYkXIBaBbx1bwy8XHgiNhHJ14BgT6x8eCQ8DFocOl+CRz4/grpGk+hYZCcsM3SZjKJqLNp4HADw2A29MSLSW3AiIqLfNiS0Gz773XC46TXYm16MR1cmoL6JhcYZsMxQG7UNJsxblYDKuibEhXfHggk8e4mIlCMu3BufPjgcLjoNdqUV4g+fJ6C2gYXG0bHMUCtJkvDs18dxKq8Svu56vDc7FjoN3yJEpCwje/rgk7nDYNSpsSutEHOXH0Il73JyaPxORa2+OJiFjYkXoVGr8M6sWN6GTUSKNaa3b5s1NHP+cxCl1Q2iY5GNsMwQAMt+Mi99mwIA+OvEfhjdy0dwIiKirhke4Y3Vj4xCd1cdki+UY+ZHB1DAjfUcEssMIaesFn9YmYAGkxm3DgzEH8b1FB2JiMgqBoV44ctHR8Pfw4DT+ZW464N9SC/g0QeOhmXGyVXXN+HhFUdQWFmPqEAP/PueIdxPhogcSp8AD2yYNwbhPq7ILqnFXcv24eC5YtGxyIpYZpyYySxhwdpEpOZWwNfdgP/MHQZ3g1Z0LCIiqwvzccXGP47B0LBuKK9txP2fHMI3SRdFxyIrYZlxUpIk4ZXNKfgxtQB6rRofPxCHkO6uomMREdmMj7sBax4ZhVsHBqLBZMaCtUl456czkCRJdDTqIpYZJ/X+zrNYvjcDAPDG3UMwNKy72EBERHZg1Gnw3n2xePi6SADAG9vT8NgXR3mek8KxzDihNYey8K8fTgMA/jFlAKYOCRaciIjIfjRqFf4xZQCWTB8EnUaFrSfycOf7e5FRVC06GnUSy4yT2Xo8F89+3XJUQa/Wn06IiJzNrBFhWPsHy51OaflVuP3dPfj5dIHoWNQJLDNOZNvJPDyxNhFmCZg1IhRPT+wnOhIRkVBx4d3x7Z+uQ2xYN1TUNeGh5YexZEsqGprMoqNRB7DMOIntKfmYv/ooGk0SpgwOwivTBvEWbCIiAAGeRqz5wyg8MDocAPDhL+dw94f7kVVcIzgZtRfLjBP4MSUfj32R0Fpklt4bA42aRYaIqIVBq8FLd0Tjgzlx8HLRITm7DJPf3o1vki7ybicFYJlxcN8kXcS8VW2LjJaHRxIRXdGt0YHYsuB6DI/ojsr6JixYm4THvjiKoqp60dHoGvhdzYEt33seC9Ymocks4Y6YYBYZIqJ26NHNBWseGYWFE/pAq7bc7RT/1i/YfCxXdDS6Cn5nc0CSJOHfP5zGi80HRz44JgJv3cMiQ0TUXlqNGgsn9MV/549FVKAHSqobMH/1UTz2RQLyynlYpdyoJAefDKyoqICXlxfKy8vh6ekpOo7N1TWa8PSGY/g2OQcA8Jf4vph/Y28u9iUi6qSGJjPe3XEG7+08C5NZgptegydv6Yu5YyKg4w+JNtOR798sMw4kv6IOj3x+BMculEOrVuGVadGYOSJMdCwiIodwMqccz/33BBKzygAAUYEeeHlaNIZHeIsN5qBYZi7hLGUmObsMf1h5BPkV9ejuqsOyOXEY1dNHdCwiIodiNktYn5CNJVtPoaymEQAwdUgwno7vhzAfnm9nTSwzl3D0MiNJEj7bl4HFW1LRaJLQx98dn8wdzr9UREQ2VFLdgNe/P4V1R7IhSYBOo8L9oyLwp5t6o7ubXnQ8h8AycwlHLjNlNQ14esMxbE/JBwBMHBiAf989BB5GneBkRETO4WROOV7degq7zxQBADyMWswb3wtzx0TA3aAVnE7ZWGYu4ahlZt/ZIjy9/hgultVCr1Hj2cn98cDocC70JSIS4Je0QizZegqpuRUAAC8XHX43NhIPjo2Alwt/wOwMlplLOFqZqapvwpItqfjiYBYAIMLHFe/OjkV0Dy/ByYiInJvZLOGb5It456d0nGs+gdvDoMXcMRF4cGwEfN0NghMqC8vMJRypzOxKK8TfNx7HxbJaAMB9I8OwaFJ/DmUSEcmIySxh8/FcvLvjDNLyqwAAeq0a02KC8dDYSPQPUvb3InthmbmEI5SZrOIavLI5Bdua18aEervgtemDMaa3r+BkRER0NWazhG0peVi28yySL5S3Pj+mlw8eHBOBm6L8uZnpNbDMXELJZaamoQkf7DyLD345h4YmMzRqFeaOjsCf4/vCjaMxRESKIEkSjmaV4dO95/H9iTyYzJZvu34eBtwVG4J7h4ci0tdNcEr5YZm5hBLLTF2jCV8czMKynekoqmoAYGnyL9w+EH0DPASnIyKizrpYVovP92dgw5ELKK5uaH1+RKQ37o4LQfzAQC4YbsYycwkllZm6RhM2JFzAuzvSkVdhOfsjzNsVf58UhYkDA3mnEhGRg2hoMmPHqXysO5yNXWmFaB6sgV6jxri+vpgyOBgTBgQ49ZpIlplLKKHMlFY3YNWBTKzYn9E6EhPkZcQTN/fBjLgQnv1BROTAcstr8VXCBWxKzmldMAwABq0a4/r64eYof9wY5Y8AT6PAlPbHMnMJuZYZSZJw7EI51h7OxteJF1DXaAYABHsZ8ci4npg1IgxGnUZwSiIisqe0/Ep8l5yD747ltt7e3SK6hyduigrAjf38MDikGzRqxx6tZ5m5hNzKTFlNA75JysHaw9mtmysBljfpI9f3xKRBQRyJISJycpIkITW3Ej+m5uOnUwU4dqEMl3639jBqMTLSG6N6+mB0Lx/0D/SE2sHKDcvMJeRQZspqGrAtJR+bj+Vib3oRmponR/VaNSZFB2LmiDCMjPTmmhgiIrqiwsp67DxdgB2nCrDnTBEq65vafLybqw5xYd0xNKwbYkK7Y3CoFzwVfrQNy8wlRJQZSZJwpqAKu04XYldaIQ6cK24tMADQP8gTM4eHYlpMD3i5KvvNRkRE9mUySziZU479Z4ux/1wxDp8vQXWDqc1rVCqgl587YkK7ITrYE1FBnugf6Kmo7zksM5ewR5mRJAmZxTVIyCzFofMl+OVMIXLL69q8JirQA5MHBWHS4CD08nO3SQ4iInI+jSYzTlwsR2JWGRKzy5CUXYrsktorvjbYy4ioIE9EBXqgb4AHIn3dEOHrJsvbwVlmLmGrMpOaW4Ff0gqRkFmKo1mlrXchtTBo1RjZ0wfj+/rhxn5+6MkCQ0REdlJUVY/k7DIkZ5chJbcSp/IqcKH0ygUHAHzc9Ij0dbM8/NwQ2t0Vwd1c0KObC/w9DELW47DMXMJWZeb170/h/Z1nW3+t16gxKMQLceHdMba3L0ZGevNuJCIiko2KukaczqvEqdwKpORW4mxhFc4XVaOwsv6an6fTqBDoZUSwl6XcBHdzgb+nAf4eBvh5GJv/a7D697yOfP9WxG4877//Pv71r38hNzcXAwcOxNKlS3H99dcLzXRdb1+cLaxCXHh3xIV3R3QPLxi0LC9ERCRPnkYdhkd4Y3iEd5vnK+sakVlcg3NF1ThfWI2M4mpcLK3FxbJa5FXUodEkIbuk9qpTVwDw4JgIvHD7QFv/Ea5K9mVm3bp1WLhwId5//32MHTsWH374IW677TakpKQgLCxMWK4xvX150CMRESmeh1GH6B5eiO7hddnHmkxmFFTWI6fMUm5yyuqQU1aLgso6FFbWo6D54edhEJD8f2Q/zTRy5EjExsZi2bJlrc/1798f06ZNw5IlSy57fX19Perr/zdkVlFRgdDQUNnsM0NERORIJEmCySxZ/QTwjkwzyXp3toaGBiQkJCA+Pr7N8/Hx8di3b98VP2fJkiXw8vJqfYSGhtojKhERkVNSqVRWLzIdJesyU1RUBJPJhICAgDbPBwQEIC8v74qfs2jRIpSXl7c+srOz7RGViIiIBJH9mhkAl+2MK0nSVXfLNRgMMBjEzt0RERGR/ch6ZMbX1xcajeayUZiCgoLLRmuIiIjIOcm6zOj1esTFxWH79u1tnt++fTvGjBkjKBURERHJieynmZ566incf//9GDZsGEaPHo2PPvoIWVlZmDdvnuhoREREJAOyLzP33nsviouL8dJLLyE3NxfR0dHYsmULwsPDRUcjIiIiGZD9PjNdJeLUbCIiIuoah9lnhoiIiOi3sMwQERGRorHMEBERkaKxzBAREZGiscwQERGRorHMEBERkaLJfp+Zrmq587yiokJwEiIiImqvlu/b7dlBxuHLTGVlJQAgNDRUcBIiIiLqqMrKSnh5eV3zNQ6/aZ7ZbEZOTg48PDyuetJ2Z1VUVCA0NBTZ2dnckM+GeJ3tg9fZPnid7YPX2X5sda0lSUJlZSWCg4OhVl97VYzDj8yo1WqEhITY9Pfw9PTkXxY74HW2D15n++B1tg9eZ/uxxbX+rRGZFlwATERERIrGMkNERESKxjLTBQaDAc8//zwMBoPoKA6N19k+eJ3tg9fZPnid7UcO19rhFwATERGRY+PIDBERESkaywwREREpGssMERERKRrLDBERESkay8yvvP/++4iMjITRaERcXBx279591dfm5uZi9uzZ6NevH9RqNRYuXHjZaz777DOoVKrLHnV1dTb8U8hfR67zxo0bccstt8DPzw+enp4YPXo0fvjhh8te99VXX2HAgAEwGAwYMGAAvv76a1v+ERTB2teZ7+cr68h13rNnD8aOHQsfHx+4uLggKioKb7311mWv4/v5cta+znw/X1lHrvOl9u7dC61Wi5iYmMs+ZvP3s0St1q5dK+l0Ounjjz+WUlJSpAULFkhubm5SZmbmFV9//vx56YknnpBWrFghxcTESAsWLLjsNcuXL5c8PT2l3NzcNg9n1tHrvGDBAum1116TDh06JKWlpUmLFi2SdDqddPTo0dbX7Nu3T9JoNNLixYul1NRUafHixZJWq5UOHDhgrz+W7NjiOvP9fLmOXuejR49Kq1evlk6cOCGdP39eWrlypeTq6ip9+OGHra/h+/lytrjOfD9frqPXuUVZWZnUs2dPKT4+XhoyZEibj9nj/cwyc4kRI0ZI8+bNa/NcVFSU9Mwzz/zm544fP/6qZcbLy8tKCR1DV65ziwEDBkgvvvhi66/vuece6dZbb23zmokTJ0ozZ87sWlgFs8V15vv5cta4znfeeac0Z86c1l/z/Xw5W1xnvp8v19nrfO+990rPPfec9Pzzz19WZuzxfuY0U7OGhgYkJCQgPj6+zfPx8fHYt29fl752VVUVwsPDERISgilTpiAxMbFLX0/JrHGdzWYzKisr4e3t3frc/v37L/uaEydO7PL/O6Wy1XUG+H6+lDWuc2JiIvbt24fx48e3Psf3c1u2us4A38+X6ux1Xr58Oc6ePYvnn3/+ih+3x/uZZaZZUVERTCYTAgIC2jwfEBCAvLy8Tn/dqKgofPbZZ9i0aRPWrFkDo9GIsWPH4syZM12NrEjWuM5vvPEGqqurcc8997Q+l5eXZ/X/d0pmq+vM93NbXbnOISEhMBgMGDZsGObPn4/f//73rR/j+7ktW11nvp/b6sx1PnPmDJ555hl88cUX0GqvfHa1Pd7PDn9qdkepVKo2v5Yk6bLnOmLUqFEYNWpU66/Hjh2L2NhYvPPOO3j77bc7/XWVrrPXec2aNXjhhRfwzTffwN/f3ypf05FZ+zrz/XxlnbnOu3fvRlVVFQ4cOIBnnnkGvXv3xqxZs7r0NR2dta8z389X1t7rbDKZMHv2bLz44ovo27evVb5mZ7HMNPP19YVGo7msKRYUFFzWKLtCrVZj+PDhTtv8u3Kd161bh4cffhjr16/HhAkT2nwsMDDQ5v/vlMRW1/nX+H7u/HWOjIwEAAwaNAj5+fl44YUXWr/J8v3clq2u86/x/dyx61xZWYkjR44gMTERjz/+OADL9LQkSdBqtdi2bRtuuukmu7yfOc3UTK/XIy4uDtu3b2/z/Pbt2zFmzBir/T6SJCEpKQlBQUFW+5pK0tnrvGbNGjz44INYvXo1Jk+efNnHR48efdnX3LZtm1X/3ymJra7zr/H9bJ1/NyRJQn19feuv+X5uy1bX+Uof5/u5/dfZ09MTx48fR1JSUutj3rx56NevH5KSkjBy5EgAdno/W20psQNouSXtk08+kVJSUqSFCxdKbm5uUkZGhiRJkvTMM89I999/f5vPSUxMlBITE6W4uDhp9uzZUmJionTy5MnWj7/wwgvS999/L509e1ZKTEyUHnroIUmr1UoHDx60659NTjp6nVevXi1ptVrpvffea3P7ZFlZWetr9u7dK2k0GunVV1+VUlNTpVdffZW3strgOvP9fLmOXud3331X2rRpk5SWlialpaVJn376qeTp6Sk9++yzra/h+/lytrjOfD9frjPfBy91pbuZ7PF+Zpn5lffee08KDw+X9Hq9FBsbK+3atav1Y3PnzpXGjx/f5vUALnuEh4e3fnzhwoVSWFiYpNfrJT8/Pyk+Pl7at2+fnf408tWR6zx+/PgrXue5c+e2+Zrr16+X+vXrJ+l0OikqKkr66quv7PSnkS9rX2e+n6+sI9f57bfflgYOHCi5urpKnp6e0tChQ6X3339fMplMbb4m38+Xs/Z15vv5yjr6ffBSVyozkmT797NKkiTJeuM8RERERPbFNTNERESkaCwzREREpGgsM0RERKRoLDNERESkaCwzREREpGgsM0RERKRoLDNERESkaCwzREREpGgsM0RERKRoLDNERESkaCwzREREpGgsM0SkKIWFhQgMDMTixYtbnzt48CD0ej22bdsmMBkRicKDJolIcbZs2YJp06Zh3759iIqKwtChQzF58mQsXbpUdDQiEoBlhogUaf78+fjxxx8xfPhwJCcn4/DhwzAajaJjEZEALDNEpEi1tbWIjo5GdnY2jhw5gsGDB4uORESCcM0MESnSuXPnkJOTA7PZjMzMTNFxiEggjswQkeI0NDRgxIgRiImJQVRUFN58800cP34cAQEBoqMRkQAsM0SkOE8//TQ2bNiA5ORkuLu748Ybb4SHhwe+++470dGISABOMxGRouzcuRNLly7FypUr4enpCbVajZUrV2LPnj1YtmyZ6HhEJABHZoiIiEjRODJDREREisYyQ0RERIrGMkNERESKxjJDREREisYyQ0RERIrGMkNERESKxjJDREREisYyQ0RERIrGMkNERESKxjJDREREisYyQ0RERIr2/wEqZHqL1bbDQwAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "B._posterior.plot(\"pdf\")" + ] + } + ], + "metadata": { + "hide_input": false, + "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.12.3" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": false + }, + "vscode": { + "interpreter": { + "hash": "3e631b8a076cc106144e9b132b7d31cae2f1e2660b47e5f9fcb0397caae5fbd5" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/skpro/regression/bayesian_proportion.py b/skpro/regression/bayesian_proportion.py new file mode 100644 index 000000000..e638fa58a --- /dev/null +++ b/skpro/regression/bayesian_proportion.py @@ -0,0 +1,172 @@ +"""Bayesian proportion estimator for probabilistic regression.""" + +__author__ = ["meraldoantonio"] + +from skpro.distributions import Beta +from skpro.regression.base import BaseProbaRegressor + + +class BayesianProportionEstimator(BaseProbaRegressor): + """Bayesian probabilistic estimator for proportions. + + This estimator uses a Beta prior and Beta posterior with a Binomial likelihood + for Bayesian inference of proportions. It provides methods for updating the + posterior, making predictions, and various utilities for analysis and visualization. + + """ + + _tags = { + # packaging info + # -------------- + "authors": ["meraldoantonio"], + "python_dependencies": ["scipy", "matplotlib"], + "capability:multioutput": False, + "capability:missing": True, + # estimator tags + # -------------- + "capability:multioutput": False, # can the estimator handle multi-output data? + "capability:missing": True, # can the estimator handle missing data? + "X_inner_mtype": "pd_DataFrame_Table", # type seen in internal _fit, _predict + "y_inner_mtype": "pd_Series_Table", # type seen in internal _fit + } + + def __init__(self, prior_alpha=None, prior_beta=None, prior=None): + """Initialize the Bayesian inference class with priors. + + Parameters + ---------- + prior_alpha : float, optional + The alpha parameter for the Beta prior distribution. Default is None. + prior_beta : float, optional + The beta parameter for the Beta prior distribution. Default is None. + prior : Beta, optional + An existing Beta distribution prior. Default is None. + + Raises + ------ + ValueError + If neither (prior_alpha and prior_beta) nor prior are provided. + TypeError + If the provided prior is not an instance of Beta. + """ + if prior is None: + if prior_alpha is None or prior_beta is None: + raise ValueError( + "Must provide either (prior_alpha and prior_beta) or prior." + ) + self.prior_alpha = prior_alpha + self.prior_beta = prior_beta + self.prior = Beta(alpha=prior_alpha, beta=prior_beta) + else: + if not isinstance(prior, Beta): + raise TypeError("Prior must be an instance of Beta.") + self.prior = prior + self.prior_alpha = prior.alpha + self.prior_beta = prior.beta + + super().__init__() + + def _fit(self, X, y): + """Fit regressor to the observed data. + + Writes to self: + Sets fitted model attributes ending in "_". + + Parameters + ---------- + X : pandas DataFrame + feature instances to fit regressor to; + will be ignored + y : pandas Series, must be same length as X + represents a series of binary experiments + whose outcome are either True (1) or False (0) + + Returns + ------- + self : reference to self + """ + assert y.apply( + lambda x: isinstance(x, bool) or x in [0, 1] + ).all(), "Values in y must be boolean or convertible to boolean (0 or 1)" + self._posterior = self._perform_bayesian_inference(self.prior, X, y) + return self + + def _predict_proba(self, X=None): + """Predict distribution over labels for data from features. + + Parameters + ---------- + X : pandas DataFrame, must have same columns as X in `fit` + data to predict labels for; + will be ignored + + Returns + ------- + y_pred : skpro BaseDistribution, same length as `X` + labels predicted for `X` + """ + y_pred = Beta(alpha=self._posterior.alpha, beta=self._posterior.beta) + return y_pred + + def _perform_bayesian_inference(self, prior, X, y): + """Perform Bayesian inference using a conjugate prior (Beta distribution). + + This method calculates the posterior Beta distribution parameters + given observed binary outcomes. + + Parameters + ---------- + prior : Beta + The prior Beta distribution from skpro distributions. + X : pandas DataFrame + Feature data corresponding to the observed outcomes `y`. + y : array-like, must be binary (0 or 1) + Observed binary outcomes. + + Returns + ------- + posterior : Beta + The posterior Beta distribution with updated parameters. + """ + n = len(y) + successes = y.sum() + posterior_alpha = prior.alpha + successes + posterior_beta = prior.beta + n - successes + return Beta(alpha=posterior_alpha, beta=posterior_beta) + + def update(self, X, y): + """Update the posterior with new data. + + Parameters + ---------- + X : pandas DataFrame + New feature instances + y : pandas Series + New labels + + Returns + ------- + self : reference to self + """ + self._posterior = self._perform_bayesian_inference(self._posterior, X, y) + return self + + @classmethod + def get_test_params(cls, parameter_set="default"): + """Return testing parameter settings for the estimator. + + Parameters + ---------- + parameter_set : str, default="default" + Name of the set of test parameters to return, for use in tests. If no + special parameters are defined for a value, will return `"default"` set. + + Returns + ------- + params : dict or list of dict, default = {} + Parameters to create testing instances of the class + """ + params1 = {"prior_alpha": 1, "prior_beta": 1} + params2 = {"prior_alpha": 2, "prior_beta": 2} + + return [params1, params2]