{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Bayes Factor for Two Different Coin Models\n", "This is taken from chapter 10 of Kruschke's book. We hypothesize two types of coins. One type of coin is hypothesized to be tail-biased and the other type of coin is head-biased. These two possibilities will form our two hypotheses and we will calculate Bayes factors to evaluate their relative credibility.\n", "\n", "**Disclaimer**: The models implemented in the notebook are intended for illustrative, teaching purposes only. Do not use these models for \"real\" analyses." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import arviz as az\n", "import numpy as np\n", "import pymc as pm" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "with pm.Model() as eqPrior:\n", " pm1 = pm.Categorical('pm1', [.5, .5])\n", " \n", " omega_0 = .25\n", " kappa_0 = 12\n", " theta_0 = pm.Beta('theta_0', mu=.25, sigma=.25)\n", " \n", " omega_1 = .75\n", " kappa_1 = 12\n", " theta_1 = pm.Beta('theta_1', mu=.75, sigma=.25)\n", " \n", " theta = pm.math.switch(pm.math.eq(pm1, 0), theta_0, theta_1)\n", " \n", " y2 = pm.Bernoulli('y2', theta, observed=[1,1,0,0,0,0,0,0])" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Multiprocess sampling (2 chains in 2 jobs)\n", "CompoundStep\n", ">BinaryGibbsMetropolis: [pm1]\n", ">NUTS: [theta_0, theta_1]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [40000/40000 00:37<00:00 Sampling 2 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 2 chains for 10_000 tune and 10_000 draw iterations (20_000 + 20_000 draws total) took 38 seconds.\n", "The acceptance probability does not match the target. It is 0.5074, but should be close to 0.8. Try to increase the number of tuning steps.\n" ] } ], "source": [ "with eqPrior:\n", " idata1 = pm.sample(10000, tune=10000)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "pm1 = idata1.posterior['pm1'].mean().item() # mean value of model indicator variable" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The posterior is provided by the estimated value of the model indicator variable, `pm1`." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Posterior: p(model 1|data) = 0.29\n" ] } ], "source": [ "print(f'Posterior: p(model 1|data) = {pm1:.2f}')" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Posterior: p(model 2|data) = 0.71\n" ] } ], "source": [ "print(f'Posterior: p(model 2|data) = {(1-pm1):.2f}')" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Posterior odds: p(model 2|data)/p(model 1|data) = 2.50\n" ] } ], "source": [ "print(f'Posterior odds: p(model 2|data)/p(model 1|data) = {(1-pm1)/pm1:.2f}')" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Bayes factor: p(model 2|data)/p(model 1|data) * p(model 1)/p(model 2) = 2.50\n" ] } ], "source": [ "print(f'Bayes factor: p(model 2|data)/p(model 1|data) * p(model 1)/p(model 2) = {(1-pm1)/pm1 * (.5/.5):.2f}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So our posterior odds are identical to our Bayes factor. This is because our prior on the model indicator variable gave equal credibility to each model." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "with pm.Model() as uneqPrior:\n", " pm1 = pm.Categorical('pm1', [.25, .75])\n", " \n", " omega_0 = .25\n", " kappa_0 = 12\n", " theta_0 = pm.Beta('theta_0', mu=.25, sigma=.25)\n", " \n", " omega_1 = .75\n", " kappa_1 = 12\n", " theta_1 = pm.Beta('theta_1', mu=.75, sigma=.25)\n", " \n", " theta = pm.math.switch(pm.math.eq(pm1, 0), theta_0, theta_1)\n", " \n", " y2 = pm.Bernoulli('y2', theta, observed=[1,1,0,0,0,0,0,0])" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Multiprocess sampling (2 chains in 2 jobs)\n", "CompoundStep\n", ">BinaryGibbsMetropolis: [pm1]\n", ">NUTS: [theta_0, theta_1]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [40000/40000 00:42<00:00 Sampling 2 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 2 chains for 10_000 tune and 10_000 draw iterations (20_000 + 20_000 draws total) took 43 seconds.\n" ] } ], "source": [ "with uneqPrior:\n", " idata2 = pm.sample(10000, tune=10000)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "pm1 = idata2.posterior['pm1'].mean().item() # mean value of model indicator variable" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Posterior: p(model 1|data) = 0.52\n" ] } ], "source": [ "print(f'Posterior: p(model 1|data) = {pm1:.2f}')" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Posterior: p(model 2|data) = 0.48\n" ] } ], "source": [ "print(f'Posterior: p(model 2|data) = {(1-pm1):.2f}')" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Posterior odds: p(model 1|data)/p(model 2|data) = 1.08\n" ] } ], "source": [ "print(f'Posterior odds: p(model 1|data)/p(model 2|data) = {pm1/(1-pm1):.2f}')" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Bayes factor: p(model 1|data)/p(model 2|data) * p(model 2)/p(model 1) = 3.25\n" ] } ], "source": [ "print(f'Bayes factor: p(model 1|data)/p(model 2|data) * p(model 2)/p(model 1) = {pm1/(1-pm1) * (.75/.25):.2f}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, the posterior odds and the Bayes factor are different because we gave more (prior) credibility to model 2. So the posterior probabilies of the two model are nearly identical, but that reflects the our priors (favoring model 2) and the likelihoods (favoring model 1) more or less cancelling each other out." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Bayes Factor with a \"Null\" Hypothesis\n", "Let's test a more traditional \"null hypothesis\". Here, we will posit two types of coins. One type is characterized by a value of theta that is exactly 0.5. We have absolute confidence that such a coin's value of theta is not .4999999999 nor .5000000001, etc. The other type of coin could have any value of theta (0-1) and all values are equally credible a priori. We then observe some data and ask whether such data should convince us that the coin is \"fair\" (H_0) or not (H_1)." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "n_heads = 2\n", "n_tails = 8\n", "data3 = np.repeat([1, 0], [n_heads, n_tails])" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "with pm.Model() as model3:\n", " pm1 = pm.Categorical('pm1', [.5, .5])\n", " \n", " theta_0 = 0.5\n", " \n", " theta_1 = pm.Beta('theta_1', 1, 1)\n", " \n", " theta = pm.math.switch(pm.math.eq(pm1, 0), theta_0, theta_1)\n", " \n", " y2 = pm.Bernoulli('y2', theta, observed=data3)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Multiprocess sampling (2 chains in 2 jobs)\n", "CompoundStep\n", ">BinaryGibbsMetropolis: [pm1]\n", ">NUTS: [theta_1]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [40000/40000 00:31<00:00 Sampling 2 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 2 chains for 10_000 tune and 10_000 draw iterations (20_000 + 20_000 draws total) took 32 seconds.\n" ] } ], "source": [ "with model3:\n", " idata3 = pm.sample(10000, tune=10000)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "pm1 = idata3.posterior['pm1'].mean().item() # mean value of model indicator variable" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.67195" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pm1" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Posterior: p(model 1|data) = 0.67\n" ] } ], "source": [ "print(f'Posterior: p(model 1|data) = {pm1:.2f}')" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Posterior: p(model 2|data) = 0.33\n" ] } ], "source": [ "print(f'Posterior: p(model 2|data) = {(1-pm1):.2f}')" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Bayes factor: p(model 1|data)/p(model 2|data) * p(model 2)/p(model 1) = 2.05\n" ] } ], "source": [ "print(f'Bayes factor: p(model 1|data)/p(model 2|data) * p(model 2)/p(model 1) = {pm1/(1-pm1) * (.5/.5):.2f}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So we have no good evidence that would allow us to choose between our 2 hypotheses. The data isn't particularly consistent with our \"null hypothesis\". A priori, the alternative hypothesis entails many credible values of theta that are much more consistent with the observed data (e.g., theta = .2). However, this alternative hypothesis also entails many values of theta that are **highly** inconsistent with the observed data (e.g., theta = .9999). So the \"null\" suffers because there is poor agreement with the data (i.e., likelihood) whereas the alternative hypothesis suffers because it is too agnostic about the possible values of theta.\n", "\n", "Let's compare this to a traditional, frequentist procedure to compare these models. We will first find the most likelihood of theta permitted each model, find the likelihood that each of these values yields, and then take the ratio of these likelihoods.\n", "\n", "To get a quick approximation of the maximum likelihood associated with our alternative hypothesis, we can plot the posterior and request a mode from the kernel density estimate." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAE5CAYAAACkkJZkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA7IklEQVR4nO3dd3xUVf7/8ddJ75WEVBIIkAACAQQEKQIKIoiKBXRVWL4uX921rVh21wKy+7Og7td1195w1VV2QUVsKEonICtNIKGXkIT03mfm/P6YEElIQhImc2eSz/PxmAfMnXvvfO4h3HfOvefeq7TWCCGEEB3NxegChBBCdA0SOEIIIexCAkcIIYRdSOAIIYSwCwkcIYQQdiGBI4QQwi4kcIRoglLqWqXUA42mXaaU0kqpy230HfFKqUVKqV62WF8z33G1UupfSqmDSimLUmpdR32XEOcjgSNE064FHjjfTBcoHlgIdFjgYN2OZGArcKoDv0eI83IzugAhRIf6jdbaAqCU2mR0MaJrkx6OEI0opZYCc4DoukNoWil1/KxZfJRS/1BK5SmlcpVSHyilghqtw00p9UelVJpSqloplamUekEp5VX3+WXA2rrZvzvrey6r+3y2UuqHuvWXKaV2KqXmtHVbzoSNEI5AejhCnOvPQBgwHJhRN60aCKz7+9+AL4BbgERgCWDGGlJnfABcDTwLbAH61a03Hrge2AH8DngZuBfYXrfc/ro/ewHLgWcACzAOeEsp5a21fs1mWyqEHUngCNGI1vqIUioXqNFabz0z/UzvA9igtb6n7u/fKqUSgTuUUnO11lopNRaYBczRWv+zbr41SqkC4AOlVLLWepdS6ky4pJ79PXU1PHXW97oA64BI4C5AAkc4JTmkJkTbfdno/c+AJ9C97v2VQA2wou7QmptSyg34tu7zcef7AqVUH6XUR0qpDKC27nUH1h6VEE5JejhCtF1Bo/fVdX961f0ZDngAZc0sH9rSypVSfsB3QAXwB+AI1gC7C5jXjnqFcAgSOELYXj5QBYxt5vPM8yw/CogDxmqt60eW1fWShHBa8gMsRNOqAe92LvsN8AgQqLX+/jzfQRPf41P3Z+2ZCUqpYOCadtYjhEOQwBGiafuBEKXUXcB/sfZYWkVrvU4p9RGwXCn1V+BHrCPN4oGrgEe01geBg4AJmFc3oKAaOIB1VFsJ8LJSaiHgCzwG5PHLSLlWUUrFYR1tB9ZDeRal1A1177drrU+0ZX1CXAgJHCGa9hZwCfAUEAScAOa2YflbgXuwnnN5FGuYHAdWA9kAWut8pdTdWHtD6wFXYEJdYF0HvIB1aHQm1qHYIVjvTNAWE4B3G037T92fvwaWtnF9QrSbkkdMCyGEsAcZFi2EEMIu5JCaEE6o7mLQln5h1Fprs73qEaI1pIcjhHN6h18uCG3q1dLoOCEMIedwhHBCSql4oFsLs5RqrQ/YqRwhWqUzBI7Tb4AQQnQiqrkP5JCaEEIIu5DAEUIIYRcSOEIIIexCAkcIIYRdSOAIIYSwCwkcIYQQdiGBI4QQwi4kcIQQQtiFBI4QQgi7kMARQghhF3K3aNGizKJKvtl7mrJqE8PighmdEIpSzd65QgghmiWBI5pkMlt4bf0RXvr+MDVmS/30IT2CeP7GwSSE+RlYnRDCGcnNO8U5as0W7l+2iy/3ZDFtUCSPTEkiPMCTz3ZmsGT1AUxmC2/cfjGX9Ao1ulQhhONp9hCIBI5ooMZk4Z6PdrB6XzZ/uiqJ+eMSGnyeXlDBr5duJ7Ookn/OG8HF8SEGVSqEcFASOF2FxaJ5d8tx3tl0jLyyaqKCvLl+aDRzL+2Jn2fLR1CrTWZ+9+EO1qTmsPDq/vz60p5NzpdTWsWs17eSV1rNh78ZyaCYoA7YEiGEk5LA6Qq01jz62V7+te0kl/YO5aKoQH7OKGbLkXy6B3jy2LT+TB8U2eRJ/4oaE3d9sIP1B3P587UXcdslcS1+V2ZRJTe+lkJ5jYll80eRGOHfUZslhHAuEjhdwYqfTrHgP7u5c3wCj1yZWB8sO08W8vjKvezNKGF0QihPXN2fpIiA+uUO55Rx38c7Sc0q4emZA5k1vEervu9kfgU3vr4FswX+c+coenbz7ZDtEkI4FQmczq6kqpaJz68jNsSHFXeOxsWl4b+52aL517YTPP/tQYoraxnZM4R+kQGkF1Sw7mAuPh6uvHTzECYkhrfpew/nlDHr9RQ83VxY9r+jiA3xseVmCSGcjwROZ/fM12m8vuEIn/9uDANjApudr7C8hg+3neCLPVmkF1QQHuDFpKRw7rosgVA/z3Z99/7MEma/kUKQjwcfz7+EqCDv9m6GEML5SeB0ZhU1Ji556nvG9g3j5VuGGlLDrvQibntrG6F+Hnw8fxQRgV6G1CGEMFyzgSO3tukEPtuZSUmViV+PjjeshuTYIJbOG0FuaTW3vLmVgvIaw2oRQjgmCZxO4P2tJ+gfGcCwuGBD6xgWF8zSeSM4VVTJXR/8RI3Jcv6FhBBdhgSOkzuUXUpqVgmzhsc6xD3OhseH8NwNg9h2rIDnVqcZXY4QwoFI4Di5L3/OQimYelGE0aXUuyY5mlsv6cFbm46x5Uie0eUIIRyEBI6T+3JPFiPiQwgPcKyT9H+6qh/xob788ZOfqTaZjS5HCOEAJHCc2OGcMg7llDFtUKTRpZzDx8ONRTMGcCK/gqWbjxtdjhDCAUjgOLF1B3IAmJjUtos17WV83zAmJoXzjx8OU1xRa3Q5QgiDSeA4sfUHc+kd7kdMsONe3f/QlERKq028u+WY0aUIIQwmgeOkKmpMbDtawGV9w4wupUX9IgO4on933tl0jNIq6eUI0ZVJ4DiprUfzqTFbGJ/o2IED8LsJvSmpMrHip1NGlyKEMJAEjpNafyAXb3dXhjvBA9CSY4MYHBvE+1tP0AlupSSEaCcJHCe17mAuoxJC8XJ3NbqUVrntkjiO5JaTcjTf6FKEEAaRwHFC6QUVnMivYFyfbkaX0mrTB0US5OPOB1tPGF2KEMIgEjhO6EwvYXRv5wkcL3dXbro4ltX7sskpqTK6HCGEASRwnNDWo/mE+HrQJ9zP6FLaZNbwWMwWzee7M40uRQhhAAkcJ6O1ZtvRAi7pFeIQN+tsi4QwPwbHBPLZrgyjSxFCGEACx8mcKqwko6iSUb1CjS6lXa4dEs3ejBIOZZcaXYoQws4kcJxMyhHr+ZtLnDRwpg+KwtVFSS9HiC5IAsfJbD2aT6ivB73bcf7miy++QCnF8ePHbV9YI9XV1SxYsIDw8HB8fX2ZNm0ax48fJ8zfkzG9u/HZzkwslobX5JSUlLBw4UJGjBhBYGAgERERXHfddRw8eLDBfOvXr2fChAmEh4fj6elJr169WLBgASUlJR2+XUKI9pPAcSJaa7YezeeSXqEOf/7m3nvvZenSpTz//PMsX76cvLw8rrjiCqqqqrhuSDQZRZX890Rhg2VOnjzJm2++yZQpU1i+fDmvv/46WVlZjBw5kvT09Pr5CgoKGDJkCC+//DKrV69mwYIFvPfee9xyyy323kwhRFtorZ391WWcyCvXcY98of+55Vi7ll+1apUG9LFj7Vu+tdLT07Wrq6t+77336qedOnVKu7u76zfffFOXV9fqfo9/rf+wYk+D5crKynRFRUWDafn5+drX11cvWrSoxe984403NKDz8/NttyFCiPZodn8tPRwDzZ07l4svvpgvv/yS/v374+Pjw7Rp0ygoKODw4cNMmDABX19fLr74Yvbs2cPWuutvkiO9uffee4mIiMDLy4vhw4fz7bffNli31ppFixYRHh6Ov78/t99+e5OHnKqqqnj44YeJjY3F09OTwYMH89VXX13Qdp2pZebMmfXToqOjGTNmDF9//TU+Hm5MGRDBVz9nUWOy1M/j6+uLt7d3g3WFhIQQFxdHTk5Oi98ZGmo9p1VTU3NBtQshOo4EjsFOnjzJE088wV/+8hfeeOMNtmzZwvz585k9ezazZ89m+fLlmEwmZs+eTcqRPLr5efDMow/w7rvv8uijj/Lpp58SGxvLtGnT2LRpU/16X3rpJRYvXsz8+fNZvnw53t7ePPzww+d8/w033MDSpUv505/+xKpVqxg+fDgzZsxg165d9fNYLBZMJlOLL7P5l6d6pqWlERMTg59fw/NM/fr1Iy0tDYAZyVEUV9ay/mBui+2Tm5vL4cOH6d+//zmfmc1mqqur2bVrF3/5y1+YOXMmERGO86htIUQjLXV/nOTltObMmaNdXV314cOH66c99NBDGmhwOOrLL7/UgB58/9t69pLlWimlly5dWv+52WzWAwYM0JMnT9Zaa20ymXRkZKS+8847G3zf5Zdf3uCQ2po1azSg161b12C+sWPH6htuuKFBnUCLr/Hjx9fPf8cdd+jBgwefs72PPvqojoyM1FprXWMy66GLv9W//fCnFtvotttu0yEhITovL++czxITE+u/f8qUKbq8vLzFdQkh7KLZ/bWbESEnfhEfH09CQkL9+969ewMwceLEc6ZlZWUx0r8IrTU33nhj/ecuLi7ceOONLFmyBID09HSysrK45pprGnzXzJkzWbNmTf37NWvWEBERwaWXXorJZKqfPmnSJJYuXVr/ftGiRdx9990tboe/v3+D900NatBa1093d3Vh2qBIlm1Pp6zahJ/nuT+Kr776Kh988AErVqyoP2R2thUrVlBcXMzPP//M4sWLufHGG+tH4gkhHI8EjsGCgoIavPfw8Dhn+plp2lRDgDbj5+eHj0/Dp3x2796diooKqqurOX36NADh4Q0fPd34fV5eHqdPn8bd3f2culxdf7kLdY8ePYiJiWlxO87eyQcHB1NUVHTOPEVFRQ2265rkaP6ZcoLVe09z/bCG6//888+55557ePbZZ7nuuuua/M4BAwYAMHr0aPr168f48eNZu3Ztg7AWQjgOOYfjRAK83biodxxlZWVUVFQ0+Cw7OxsfHx88PT3rz2M0PtHe+H1ISAjR0dFs3779nNfWrVvr55s3bx7u7u4tviZNmlQ/f1JSEunp6ZSXlzf4vrS0NJKSkurfD+0RRGyI9zkXgW7ZsoXZs2dz55138tBDD7WqbYYOHQrA0aNHWzW/EML+pIfjBHTdQ8sSu/szYsQIlFIsX76c22+/vf7z5cuXM2bMGABiY2OJiIhg5cqVXHnllfXr+eSTTxqsd9KkSbzwwgv4+fk1CILG2npIbfLkyQB8+umn3HrrrQBkZmayceNGXnnllfr5lFJcMziaV9YdJre0mjB/T/bt28f06dO58soreemll87bNmds3rwZgJ49e7Z6GSGEfUngOIGMokoAEiMC6NevHzfffDN33303JSUl9O7dmzfffJO0tDReffVVwHo47OGHH+bBBx+kW7dujB07lhUrVpCamtpgvVdccQVTpkzhiiuu4JFHHmHAgAGUlJSwa9cuqqqqePrppwHreab4+PhW1xsTE8P//M//cP/996O1JiwsjEWLFhEXF1cfQACLFy9m8eLFxDy4ki/2ZDKtjy9XXnklfn5+3Hvvvfz444/18wYEBNSPVLvtttvo27cvycnJ+Pj4sGPHDpYsWcKoUaOYMGFCu9pYCGEHLY0ocJKX05ozZ44eNmxYg2nvvvuuBnRpaWn9tJc+26QB/fr7y7TWWpeXl+u7775bh4eHaw8PDz1s2DD9zTffNFiPxWLRjz32mO7WrZv28/PTt9xyi/7www/PufCzqqpKP/HEEzohIUG7u7vr7t276ylTpugvvvjigratqqpK//73v9fdunXTPj4+eurUqfro0aMN5lm4cKEG9NQXN+gZ/9ik165d26pRcC+99JIeOnSoDggI0L6+vvqiiy7SixcvbtBmQgjDNLu/Vlo7/TPmnX4Dzuf3y3ax8VAu2x+9vFOOwHpjwxGe+iqNdQ9eRnw3X6PLEUJcmGZ3UjJowMFpbb1/2kgnuH9ae109OAqlYOUueTCbEJ2ZBI6DO1lQQVZxldM+jqA1IgO9GdkzhJW7MugEPW4hRDMkcBzcmfunjeoVYnAlHeva5GiO5pWzN0MeMSBEZyWB4+BSjuTTzc+DhLC2P//GmUy9KBIPVxd5MJsQnZgEjgPTWrPlSOc+f3NGoI87lyWGsWp3JmaLHFYTojOSwHFgB7PLyCmtZmzvbkaXYhfXDokmp7S6/jCiEKJzkcBxYBsPWW/dP6ZP1wiciUnh+Hm68dlOOawmRGckgePANh3Oo2c3X2KCfc4/cyfg5e7KlRdF8M3e01TVms+/gBDCqUjgOKhqk5ltRwsY00UOp51xTXIUpdUm1qa1/IRPIYTzkcBxUDtOFFFZa+4yh9POGJ3QjW5+njJaTYhOSALHQW06nIuri2JUQue94LMpri6KqwdHsjYtl+LKWqPLEULYkASOg9p0KI/BMYEEeJ37cLTO7trkaGrMFr76OcvoUoQQNiSB44CKKmrYk1HMmD5hRpdiiEExgSSE+fLJjlNGlyKEsCEJHAe0+XA+WsPYLnb+5gylFDOHxrD9eCEn8svPv4AQwilI4Dig79OyCfR2Z0hskNGlGOa6IdEoBZ/skMEDQnQWEjgOxmzRrDuQy4TEMNxcu+4/T1SQN6MTQvlk5ym5g7QQnUTX3aM5qF3phRSU1zCxX3ejSzHc9UNjSC+oZPvxQqNLEULYgASOg1mTmoOri2J83645YOBsUwZE4OPhyoqfZPCAEJ2BBI6D+SE1h+HxwQR6d73h0I35erox9aJIvvw5S251I0QnIIHjQNILKjiQXcrlcjit3vVDoymrNrF632mjSxFCXCAJHAfyQ939wyYmhRtcieO4pFcoUYFeMlpNiE5AAseBrEnNplc3X3p18qd7toWLi+K6odFsPJRLdkmV0eUIIS6ABI6DKKs2se1oAZP6Se+msZlDY7BoWCk39BTCqUngOIhNh3KpMVuYmCTnbxpLCPMjOTaIFT9lyDU5QjgxCRwHsSY1hwAvNy6ODza6FId0/bAYDmSXsi+zxOhShBDtJIHjACwWzdq0HMYnhuPehe8u0JKrB0Xi4erCCrmhpxBOS/ZuDmDXqSLyy2u4XM7fNCvIx4NJ/cL5fFcmtWaL0eUIIdpBAscB/CB3F2iV64fGkF9ew/oDuUaXIoRoBwkcB7AmNZthccEE+XgYXYpDG58YRqivhxxWE8JJSeAYLKOokrTTpXI4rRXcXV2YkRzF96k5FFXUGF2OEKKNJHAM9kNqNoAMh26l64fGUGO28MUeefy0EM5GAsdga1JziA/1ISHM1+hSnMKAqAASwnz5fHem0aUIIdpIAsdA5dUmUo7kM6lfd5RSRpfjFJRSXJMczY/HCsgoqjS6HCFEG0jgGGjT4TxqzBYmyc062+Sa5CgAVkkvRwinIoFjoB9Sc/D3dGN4zxCjS3EqcaG+JMcGsXKXBI4QzkQCxyAWi+b7tBzGJYbJ3QXa4ZrkKFKzSjiYXWp0KUKIVpI9nUF+zigmr6xahkO307RBkbgo+Fx6OUI4DQkcg3yfloOLgvF9JXDaI9zfi0t7d2PlbrmDtBDOQgLHID+kZTO0RzAhvnJ3gfa6Jjma9IJKdpwsMroUIUQrSOAY4HRxFXszSpgoh9MuyJQB3fF0c+FzeTCbEE5BAscAaw/kADBJ7i5wQfy93Lm8X3e+2JOFSe4gLYTDk8AxwPepOUQHedO3u5/RpTi9GclR5JfXsPlIvtGlCCHOQwLHzqpqzWw+nMekfuFydwEbuCwxDH8vN1bulMNqQjg6CRw7SzmaT2WtmYlydwGb8HRz5aqLIlm97zSVNWajyxFCtEACx85+SM3B292VS3qFGl1Kp3HNkCjKa8ysqbvzthDCMUng2JHWmrUHcri0dze83F2NLqfTGNkzlO4BnnKrGyEcnASOHR3Pr+BUYSXjE+VR0rbk6qK4elAU6w/Kg9mEcGQSOHa06VAuAGN7dzO4ks7n2iHR1Jo1X/182uhShBDNkMCxow2H8ogN8SYu1MfoUjqdMw9mWykXgQrhsCRw7MRktrD1SD5jeofJcOgOcObBbNuOFZApD2YTwiFJ4NjJ7lNFlFabGNdHDqd1FHkwmxCOTQLHTjYczMNFwegECZyOcubBbJ/JaDUhHJIEjp1sPpzHwJggAn3cjS6lU7tWHswmhMOSwLGDyhozu08VMUou9uxw0wZF4eqiZPCAEA5IAscOdqYXUmvWjOwZYnQpnV6Yv6f1wWy7MuXBbEI4GAkcO9h2tAAXBcPig40upUu4ZnAUpwor+e+JQqNLEUKcRQLHDn48VkD/qAACvOT8jT1ceVEEvh6u/Ht7utGlCCHOIoHTwapNZnacLGREvJy/sRdfTzeuHhzFlz9nUVZtMrocIUQdCZwO9vOpYqpNFkb2kvM39nTT8Fgqasx8IdfkCOEwJHA62LZjBQAMj5fAsachsUH0Dvdj2X/lsJoQjkICp4NtP15An3A/Qnw9jC6lS1FKMeviWHaeLOKQXJMjhEOQwOlAWmt2nixiWJyMTjPCdUOjcXNRLJPBA0I4BAmcDnQsr5ziylqG9AgyupQuqZufJ5MHdOc/P52iokYGDwhhNAmcDrTzZBEAQ3pID8cov760J8WVtXyyQ+48IITRJHA60K70Ivw93egd5md0KV3WxXHBDIoJ5J3Nx7BY5M4DQhhJAqcD7UwvZFBsIC4u8vwboyilmHdpT47mlrO+7omrQghjSOB0kMoaM6lZpQyJlcNpRrtqYCTdAzx5Z9Mxo0sRokuTwOkgP2cUY7ZoGTDgADzcXLh9VDwbD+WRdrrE6HKE6LIkcDrIzpPWG0cmxwYZW4gA4Fcje+Dr4crLa48YXYoQXZYETgfZebKIuFAfQv08jS5FAEE+Htw+Op4v9mRyOEcuBBXCCBI4HUBrzY6ThQyR3o1DuWNMT7zcXPn7D4eNLkWILkkCpwNkFVeRU1ot1984mFA/T24fFceq3ZkcyS0zuhwhuhwJnA6w51QxAINiAg2uRDR2x9heeLi58LL0coSwOwmcDrA/sxhXF0W/yACjSxGNhPl7cuvIOD7blSG9HCHsTAKnA+zNLKF3mB9e7q5GlyKacOdlCXi6ufK3NYeMLkWILkUCpwPszShmQJT0bhxVNz9P5oyOZ9WeTA7KowuEsBsJHBvLKbUOGBgQLedvHNn8cb3wcXflxTUHjS5FiC5DAsfG9mVar2S/SHo4Di3E14N5Y3ry1c+n2ZdZbHQ5QnQJEjg2ti/DuvPqL4Hj8O4Y0wt/LzdelHM5QtiFBI6N7cssIT7UB38vd6NLEecR6OPOHWN68d3+bPacKjK6HCE6PQkcG9ubWSznb5zIvDHxBHq783/fybkcITqaBI4NFVfUkl5QKSPUnIi/lzvzx/Vi7YFcfjpRaHQ5QnRqEjg2tC/Lev7moijp4TiTuaPjCfX1kBFrQnQwCRwb2pdhHaEmPRzn4uvpxp3jE9h4KI8fjxUYXY4QnZYEjg3tzSwmKtBLHknghG69JI4wf09e+PYAWmujyxGiU5LAsaG9GcX0l8NpTsnbw5XfXpbAtmMFpBzNN7ocITolCRwbqagxcTSvnIui5XCas7p5RA/C/T35+/dyJ2khOoIEjo2kZpWgtQwYcGZe7q787/gEUo7my7kcITqABI6N7D0zYEB6OE7tlhE96ObnwUvfy90HhLA1CRwb2ZdZTKivBxEBXkaXIi6At4cr88f1YtPhPH46Ib0cIWxJAsdG9maUMCA6EKWU0aWIC3TrJXGE+HrwkpzLEcKmJHBsoNpk5mB2qVx/00n4eLjxm7G9WH8wl13pRUaXI0SnIYFjA4eyyzBZtAwY6ERuGxVHkI87f5dzOULYjASODeyteySBDInuPPw83bhjTE++T8up//cVQlwYCRwb2JtZjL+nG7HBPkaXImzo9tHxBHi58cK3B4wuRYhOQQLHBvZlltA/KgAXFxkw0JkEeLnzuwm9WXsgl82H84wuRwinJ4FzgUxmC6lZJVwkz8DplOaMjicm2Ju/fJmK2SL3WBPiQkjgXKCjeeVU1VpkhFon5eXuysNXJpGaVcKKHaeMLkcIpyaBc4H2ZVpPKA+QEWqd1tWDIkmODeL51QcorzYZXY4QTksC5wLtzSjBy92FhDBfo0sRHUQpxePT+5FTWs3fZJi0EO0mgXOB9mYU0y8yADdXacrObFhcCDePiOXtTcfqe7VCiLaRveQFsFg0+zNL5ILPLuKRK5MI9nHnT5/8LAMIhGgHCZwLcLKggtJqk1zw2UUE+Xjw+PT+7D5VzOsbjhhdjhBORwLnAuyVAQNdzozBUVw1MIL/++4gqVklRpcjhFORwLkAezNKcHdV9O3ub3Qpwk6UUvzl2oEEenvw+2W7qDaZjS5JCKchgXMB9mUWkxjhj4ebNGNXEuLrwbPXDyTtdCkvrpFRa0K0luwp20lrzd6MYhkw0EVN6ted2cNjeX39Ef57XB7UJkRrSOC0U2ZxFYUVtQyQW9p0WY9N709UkDcL/rObihq5IFSI85HAaaf6RxLILW26LD9PN56/cTAnCyp45us0o8sRwuFJ4LTTvoxiXF0U/SIlcLqyS3qFMu/Snvwz5QSbDskdpYVoiQROO+3NLKF3mB9e7q5GlyIM9tCURBLCfHlo+W5KqmqNLkcIhyWB0057M4oZIBd8Cqx3lP7rTcnklFbz5Of7jS5HCIclgdMOOSVV5JRWywi1DvbZZ58xaNAgPD096dmzJ3/9619bnP/+++9HKcWDDz7YYHpaWhojR44kMDCQ2bNnU1ZW1uDzDRs2EB0dfc70pixduhSl1DnzDo4Nonf6V7w4dxzrDuQAcPz4cZRS9S9fX18SEhL41a9+xcaNG89Z99y5c7n44ovPW4MQzkoCpx32ZVqvMJeHrnWczZs3M3PmTEaMGMGqVauYN28ejzzyCC+++GKT8+/fv5933nmHgIBze51z586ld+/e/Pvf/2b//v089dRT9Z9ZLBbuv/9+nn76afz8/C6o5uE9Q3BV8NhnexuMWnv++edJSUnhq6++4vHHHyc/P59x48bx5JNPXtD3CeFs3IwuwBntOVWMUtBfRqh1mMWLFzNmzBjeeustACZPnkxhYSGLFy/mt7/9LR4eHg3mv/fee7nvvvt4//33G0wvKytj27ZtrFq1irCwMIqKinj++efrQ+ftt9/G3d2d22677YJrdnNxwcfTjVOFlfxtzSFu7u8NQGJiIpdccgkA48ePZ+7cuTzxxBMsWrSI8ePHc9lll13wdwvhDKSH0w670gvpE+6Hn6fkdUfZtWsXl19+eYNpZ0InJSWlwfTly5eTmprKH/7wh3PWU1NTA4C3t3Xn7+PjUz+tpKSExx9/nL/97W8opWxSt5uLYvbwWN7adIyD2c3fa23hwoVERUXx2muv2eR7hXAGEjhtpLVmV3oRg2OCjC6lU6uqqjqnF+Pp6QlAampq/bTKykoWLFjAM888g6/vuQ/BCwkJoWfPnvz973+noKCAN954o/48yZ///Gcuv/zy+t5HW5jNZkwmU4OXxWIB4I9T+xHs486Sbw40u7yrqysTJ05k69atbf5uIZyV/IreRicLKiisqCW5R5DRpXRqvXv3Zvv27Q2m/fjjjwAUFPxyK5mnn36ayMhIbr311mbX9fLLL3PjjTfypz/9iT59+vDyyy9z+PBh3n77bfbs2dOu+oKCgpqcHhoaSqCPO49P78/vXl/d4jpiYmLIzs5u1/cL4Yykh9NGu9KLAEiODTK0js7uzjvvZOXKlbz55psUFhayevVqXnjhBcDaOwA4duwYzz//PC+++GKLh8SmTp1KTk4OBw4cIDU1lR49evDAAw/w+9//npiYGF5++WV69OhBjx49eOWVV1pV34YNG9i+fXuD129+85v6z2cMjmJI3c9IWXXT1+ZoLQ9xE12L9HDaaFd6Ed7uriTKIwk61Lx589i9ezd33XUX8+fPx8fHh2effZZ77rmH7t27A/CHP/yBqVOnkpSURFFREWAddVZdXU1RURGBgYH1QeTj40Pfvn0BWLNmDbt372bZsmXs3r2bxx9/nC1btgAwatQoxowZw6BBg1qsb8iQIeeMavviiy/q/66U4v4r+vDFo/DZzkxmX3/uOjIyMuq3RYiuQHo4bbQrvYiB0YG4uUrTdSRXV1f+8Y9/kJuby549e8jOzq4/13LmzwMHDvDJJ58QHBxc/0pPT+cf//gHwcHBZGRknLNek8nE/fffz5IlS/D29mbdunVMnDiRpKQkkpKSmDRpEuvXr7fJNiSEWX8pWXcgh7TTDQcQmEwmfvjhB0aNGmWT7xLCGUgPpw1qTBb2ZZYwZ1Sc0aV0GWeCBOCVV15h9OjRJCUlAfDWW2+dcwHm7NmzGT9+PHfddRdhYWHnrO+1114jODiYWbNm1U+rqKio/3t5ebnND3X5eLiy6PN9fPSbS+p7XIsXLyYzM5M777zTpt8lhCOTwGmD1KwSakwWkmODjS6l09u6dSubNm0iOTmZkpISPvroI1avXs2mTZvq52nqqnwvLy9iY2ObvLalsLCQJ598ktWrfzmZP27cOB5++GHeeecdAH744QeeeeYZm27LJaE1rN24hSXup4lwKeHjjz/mm2++qb8OR4iuQgKnDeoHDMgItQ7n7u7OsmXLWLRoES4uLowdO5bNmzczcODAdq9z4cKFzJgxg6FDh9ZPGzJkCEuWLOHRRx8FrHcFGDx48AXXf7aP/v7/APjjxx7ExUQzevQoNmzYwNixY236PUI4OtUJRsrYbQN+v2wXmw7n8eOfJtnsQkHRNWw7ms+sN7bywBV9uXdSH6PLEaIjNbtzlDPfbbArvYjk2CAJG9FmI3uFMvWiCF5dd4TTxVVGlyOEISRwWqmgvIZjeeVy/Y1otz9O7YfZolmyWp4OKromCZxW2n7cenX7iJ4hBlcinFWPUB/mjenJJzsy2F13PlCIrkQCp5W2HyvAw82FQTHySALRfr+bkEA3Pw8Wf7Ff7jQguhwJnFb68XgBybFBeLrJI6VF+/l7ufPg5ER+OlHIF3uyjC5HCLuSwGmF8moT+zJLGBEvh9PEhbvx4lj6RQbwzNdpVNWajS5HCLuRwGmFHScLMVs0w+X8jbABVxfF49P7kVFUydubjhldjhB2I4HTCtuPFeCiYFic3GFA2MbohG5MGdCdl9ceJqdEhkmLrkECpxW2HStgQFSgPOFT2NSfrupHrdnCc6ubf1CbEJ2JBM55VNaY2ZlexEg5nCZsLC7Ul3mX9mT5jlP8fKrY6HKE6HASOOex9Vg+NSYLY/uee+dhIS7U7yb2JsTHgz/LMGnRBUjgnMfGg3l4urlID0d0iAAvdxZMTuTH4wV8vfe00eUI0aEkcM5jw6FcRvQMwctdrr8RHWPW8FiSIvx56qtUGSYtOjUJnBZkFlVyOKeM8XI4TXQg6zDp/pwqrOTltYeNLkeIDiOB04INB3MBGCeBIzrYpb27MXNoNK+sO8LeDBlAIDonCZwWbDiUS0SAF33C/YwuRXQBT0zvT7CPBw8v30Ot2WJ0OULYnAROM6pqzaw/kMtliWHy/BthF0E+Hvzl2ovYn1XCa+uOGF2OEDYngdOMTYfyKK8xM3VgpNGliC7kyosimD4okr99f4ifThQYXY4QNiWB04yv9mYR4OXGqF6hRpciupj/d91AooK8uftfOykorzG6HCFsRgKnCTUmC2v2Z3NF/wg83KSJhH0Fervzyq+Gkl9Wwz0f7aDGJOdzROcge9MmfJ+aTUmViemD5XCaMMZF0YE8PXMgmw/ns+A/u7FY5C4EwvnJ3Sib8PH2dCIDvRjXR4ZDC+NcPyyG3LJqnvk6DV8PV/7fdQNxdZEBLMJ5SeA0kllUyYZDudw9obf85xaG+99xvSivNvH3Hw6TV1bDX2cNJsDL3eiyhGgXOaTWyHspx1HATRfHGl2KECilWDA5kSdnDGDtgRymvbSRH9Ky5UafwilJ4JylsLyGD1JOMH1QFLEhPkaXI0S9OaPj+ff/XoK7iwvzlv6X297+kbVpOZjl3I5wIqoT/KZksw1Y8k0ar6w7wur7x5EY4W+r1QphMzUmC+9vPcGr646QV1ZNVKAXM5KjuSY5iqQIf7lIWTiCZn8IJXDq7M8sYcY/NnH14Cj+b1ayLVYpRIepMVn4bn82y39KZ8OhPMwWTd/uflyTHM2MwdJDF4aSwGlJSVUtN76aQn55Dd/9fhzBvh62qEsIu8gvq+arvadZuTOD/54oxEXBjMFR3DOpDwlhch9AYXcSOE05mlvG5sN5vLvlOOkFFbwzdzhjZSi0cGLpBRW8v/UE76ecoNpk5toh0TxwRV9igqXHI+xGAqcp76cc5/GV+4gL9eHpmQMZndDNlnUJYZi8smpeX3+E91JOgIbbR8Xxuwm9pfcu7EECpylFFTWU15iJDPDCRa65EZ1QZlElL645yPKfTuHr4cadlyXw60vj8fGQS/BEh5HAEaIrO5hdypJvDrAmNRt/TzemD45kyoAIhseH4Osp4SNsSgJHCAE/nSjkX9tO8tXPWVTWmnFzUfTs5kvvcD8SwvxIjPBnQFQA8aG+0usX7SWBI4T4RUWNiZ9OFLLtaAFpp0s5mlvGiYKK+gtJ/T3dGNc3jMkDujNlQARe7q4GVyyciASOEKJl1SYzh3PK2JdRwo6ThaxJzSGvrJpgH3dmj+jBHWN6EurnaXSZwvFJ4Agh2sZi0Ww5ks/7W4/z3f5svNxdmTM6nvlje8loN9ESCRwhRPsdzinjpe8PsWpPJj7urtw6Ko7/GdOTcH8vo0sTjkcCRwhx4Q5ml/L3Hw7z5Z5M3FxduOniGG4e0YP+kQHtuo+bxaIxa43ZotEaPNxcOu1jQfLLqjmUU8aJ/HKO51eQWVRJebWJGrMm0NudiABPkmODGd4z2NmDXAJHCGE7x/PKeX3DEZb/dIpas6ZnN18ujgumX2QAoX4eeLu7UllrpqLGTFFFLQXl1eSX1ZBXXlP/9/zymiYfn+3j4Yqvpxv+Xm7EBPsQH+pDjxAf4kN9ie/mQ0ywj0MPYtBaczSvnN3pRaSdLiU1q4S006XkllbXz+PuqogM9Mbfyw03F0VJlYnMokqqTRaUgpE9Q5g+KIppAyOd8fClBI4QwvYKymv4Zu9pvtt/mj2niskvr2lyPi93F0J9PQn18yDU14NQP09CfT3w9nDFVSlcXBQuSlFVa6a82kRZtYniylrSCys4kV9BaZWpfl1KQWSAFz1CfQjy9sDDzQVPNxdclEJj7SlpwFL3Fw83F6KDvIkJ8SY22Ife4X4E+dh2J34yv4L1B3NIOZrPj8cKyCuztoOHmwt9u/uRFBFAUoQ/fbv707ObL5GBXri5Nnw6TK3ZQmpWCWvTcvl8dwZHcstxd1VMTArn+qExXJYYjoebUzxRRgJHCNGxtNbkl9dQVFFLVa0ZL3dXfDxcCfJxv6A7G2itKaqo5Xh+OScLKjieV1F3WKqcsmoTNSYL1SYLFq1RKJSy7vHOHOKrNpnrA+CMcH9PEiP86RPuT9/ufvSN8KdPuB/+rXyaalFFDTtOFrLhYB7rD+ZyLK8cgOggb0b2DGFEzxCGxgXTq5vvOcHS2m3en1XCpzsy+GxXZv1owRmDo5iQFM6wuOBmaz3TXsfyyzmZX0F6QQUnCypIL6ygoLyG8mozFTXWdlPK2l4uZ/05c0g0j03v3+aazyKBI4TouqpqzWQUVXIiv5xD2WUczC7jYHYph3JKqar95bBedJA3vcJ8CfX1IMjHg0Bvd1xdFCazhdyyGnJLq+rOw1QA1p7bqF6hjO8bxvjEcOJDfWz+TCKT2cLGQ3ms2HGKb/dn1wWFtZcXHeyNr6cbbi4ulFbVUlxZS2ZRJSVn9QjBGrCxIT508/PA19MNXw83PNxc6nqD1l6hRWssWjMsLpjrhsRcSMkSOEII0ZjFojlVWMmB7FJrAGWXciy/gsLyGgrLayittu64lYJQXw+6+XkSH+rLoNhAkmOCGBoXbNfzSeXVJnaeLGLHyUKO55eTUVhJZa2ZWrPG38uNIG93wgOsNfbs5ktcqCHnvCRwhBCirc7ceUGB3Oqn9ZptqFYdXFRK9VdKfa+UqlBKZSqlFiulzhuZSqlApdS7SqlCpVSxUupDpVRoE/Ndo5T6WSlVpZTar5Sa1ehzD6XUc0qpjUqpSqVUkyEzd+7cumOSDV9paWmt2UwhhGjA1UXh6qLOGzYrV65k4MCBeHl50b9/f5YtW9aq9X/88ccMHToUPz8/oqOjuf3228nMzKz/PCsri4ceeojBgwfj5+dHbGwsc+bMaTCPMzlv4CilgoE1WHsS1wCLgQXAk61Y/zLgMuAOYC4wHPis0frHACuAtcBU4EvgI6XU5LNm86lbRwWwpaUvTEpKIiUlpcErPj6+FaUKIUTbbdq0ieuvv54JEybw9ddfM23aNG6++Wa+/fbbFpf7/PPPufnmmxk9ejQrV67k2WefZcOGDUyfPh2LxXpe6aeffuLTTz/l5ptvZtWqVTz33HNs27aN0aNHU1ZWZo/Nsy2tdYsv4I9AIRBw1rSHse78A1pYbhTWkBp31rQRddMuP2vaauCHRst+BWxqNO3M4b+7rWXXf1Zvzpw5etiwYVoIIexl8uTJesKECQ2mTZ06VV966aUtLjdr1iw9dOjQBtNWrlypAb1//36ttdaFhYW6tra2wTwHDhzQgF66dKkNqu8QzeZJaw6pTQVWa61Lzpr2MeANjD/Pctla6w1nhduPwLG6z1BKeQITgH83WvZjYJRSKvCsZeVcjRDCoVRXV7N27VpuuummBtNnz55NSkoKxcXFzS5bW1tLYGBgg2lBQUEAZ37JJigoCDe3hkPK+/bti4+PDzk5OTbYAvtqTeAkAQ1OgmitT2Lt4SS1Zbk6qWctlwC4NzFfal1tfVtRXwP79+8nICAAT09PxowZw/r169u6CiGEaJUjR45QW1tLUlLDXWG/fv2wWCwcPHiw2WXnzZvHxo0b+ec//0lJSQkHDx7kscceY8KECfTv3/x1MHv27KGioqLFeRxVawInGChqYnph3WcXstyZPxvPV9jo81YZMmQIL7zwAqtWreLDDz/EbDZzxRVX8OOPP7ZlNUII0SqFhdZd1ZmeyRnBwcENPm/KtGnTWLp0KfPnzycwMJDExETMZjOffPJJs8tYLBbuu+8++vTpw+TJk5udz1G19vLfpg5nqWamt2e5xu9VM9NbdN999zV4P23aNPr3789TTz3FZ5991pZVCSFEqzW+2PPMIbGWLgJdu3Ytd955J/fddx9Tp04lOzubRYsWcd1117FmzRpcXc8dCPzHP/6RlJQU1q9fj7t76+6K4EhaEziFQFAT0wNpugdz9nJhTUwPOmu5wrOmNZ6H86z/vLy9vbnqqqtYtWrVhaxGCCGadKYnU1RU1GD6mfeNez5nW7BgATNmzODZZ5+tn5acnExSUhIrV65k5syZDeZ/5ZVXeO655/joo48YOXKkTeq3t9YcUkuj0bkapVQs4EvT52iaXa7O2ed2jgC1TcyXBFiA5g+AtoGtbzUhhBAACQkJuLu7n3OtX1paGi4uLvTt2/xp6LS0NJKTkxtMS0xMxNvbmyNHjjSYvmLFCu655x6WLFnCrFkNLlN0Kq0JnK+BKUop/7OmzQIqgZbOyH8NRNRdZwOAUupioFfdZ2itq7Fef3Njo2VnASla6+aHeLRCZWUlX3/9NcOGDbuQ1QghRJM8PT2ZMGEC//nPfxpMX7ZsGaNGjTpnFNrZ4uLi2LFjR4NpqampVFZWNrh2cN26dfzqV7/i7rvv5sEHH7Rp/XbX0pjpuuOQwUAW8B1wOTAfKAP+0mi+w8DbjaZ9AxwFZgLXAgeAjY3mGQOYgBexXiS6BGvvZnKj+aYCNwBvYT23cwNww/Hjx7XWWhcVFekxY8bo1157Ta9Zs0Z//PHHeuTIkdrDw0Nv3769A4ecCyG6so0bN2pXV1d933336bVr1+qHHnpIK6X06tWr6+c5fvy4dnV11e+99179tBdffFErpfQDDzygv/vuO/3BBx/ovn376vj4eF1WVqa11nr//v06MDBQDx48WG/evFmnpKTUvw4fPmz3bW2l5vOkpQ/1Lzv7/sAPWHs1WcCfAddG8xwHljaaFgS8i/VcTAnwL6BbE+u/FtgLVGM93Da7iXmO1wVNg9e7776rtda6srJSX3fddTomJkZ7eHjogIAAPWXKFJ2SktLhrSuE6No+/fRTPWDAAO3h4aETExP1Rx991ODzY8eONdhfaa21xWLRr7zyih44cKD28fHRUVFR+qabbtJHjhypn+fdd989Z5935jVnzhw7bV2bNZslcvNOIYQQtnRhN+8UQgghLpQEjhBCCLuQwBFCCGEXEjhCCCHsQgJHCCGEXUjgCCGEsAsJHCGEEHYhgSOEEMIuJHCEEELYRWufh9NpLVq0iCeffNLoMoQQDm7hwoUsWrTI6DKcmvRwhBBC2IUEjhBCCLuQm3cKIYSwJbl5pxBCCGNJ4AghhLALCRwhhBB2IYEjhBDCLiRwhBBC2IUEjhBCCLuQwBFCCGEXEjhCCCHsQgJHCCGEXUjgCCGEsAsJHCGEEHYhgSOEEMIuJHCEEELYhQSOEEIIu3D6xxM8+eST3wDdLmAVUUCmjcpxVtIG0gYgbQDSBmdcSDvkLVy48MomP9Fad+nXokWLtNE1GP2SNpA2kDaQNrBHO8ghNSGEEHYhgQNPGl2AA5A2kDYAaQOQNjijQ9rB6c/hCCGEcA7SwxFCCGEXEjhCCCHsQgJHCCGEXUjgCCGEsItOHThKqd8qpY4ppaqUUj8ppcaeZ/6BSqn1SqlKpVSGUuoJpZSyV70dpS3toJS6TCm1UimVpZSqUErtUUrNs2e9HaGtPwtnLddHKVWqlCrr6Bo7Wjv+Pyil1P1KqTSlVHXdz8Qz9qq3I7SjDaYopVLqfgby6v5v9LVXvbamlBqnlPq8bv+mlVJzW7GMzfaLnTZwlFKzgL8BTwFDgC3A10qpHs3MHwB8B2QDw4F7gYeAB+xScAdpazsAo4GfgRuAi4BXgTeUUrfYodwO0Y42OLOcB/AxsKHDi+xg7WyDF4DfAo8A/YCrcOK2aMc+oSewEthYN//lgDfwlV0K7hh+wF7gPqDyfDPbfL9o9BWtHfUCtgFvNpp2CHi6mfnvAkoA77OmPQZkUDd83BlfbW2HZtbxb2CF0dti7zYA/g94F5gLlBm9HfZsAyARqAX6GV27gW1wA2AGXM+aNgHQQDejt8cG7VEGzD3PPDbdL3bKHk7db6bDgG8bffQt1t/gmzIK2Ki1Pjv1V2O9p1C8rWu0h3a2Q1MCgEJb1WVP7W0DpdQ0YDrW3+icWjvb4BrgKHClUuqoUuq4Uuo9pVR4B5baYdrZBv/FGrp3KKVclVL+wBxgu9Y6r8OKdSw23S92ysDBejNPV6zdwLNlAxHNLBPRzPxnPnNG7WmHBpRS04FJwBu2Lc1u2twGSqlI4E3gNq11aceWZxft+TnoBcQBs7H28G4DkoBVSiln3G+0uQ201seBK7BedV8NFAMDsf4i0lXYdL/ojD84bdH4NgqqiWnnm7+p6c6mre1gnUmpS4F/AfdqrX/siMLsqC1t8AHwqtZ6a8eWZHdtaQMXwBNr6G7QWm/EGjojsB7Ld1atbgOlVATwNvBPrNt8GVAK/NtJQ7e9bLZf7KyNlof12GvjBA7n3LQ+43Qz89PCMo6uPe0AgFJqDPA18ITW+tWOKc8u2tMGE4GFSimTUsqEdafjW/d+fseV2mHa0wZZgElrffCsaYcAE9DiYAsH1Z42+B1QrrV+WGu9U2u9AbgVGE/bDkk7M5vuFztl4Gita4CfsHaHz3YF1pEpTUkBxiqlvBrNnwkct3WN9tDOdkApNQ5r2DyptX6xwwq0g3a2wUAg+azXE1hH9CQD/7F9lR2rnW2wGXBTSiWcNa0X4AacsHmRHaydbeCDNaTOduZ9p9x3NsG2+0WjR0p04AiMWUANcAfWIZ1/wzoqI67u86eB78+aPxBrmn+MdTjwTKyjMxYYvS12bofLgHLgOay/2Zx5hRm9LfZqgyaWn4vzj1Jr68+BC9Yd9HqsQ4KH1P19K+Bi9PbYqQ0mAhZgIdAHGAp8A5wEfI3enna2gR+//CJVgfWXqWSgRzNtYNP9ouEN0MGN+1usKVxd959n3FmfLQWON5p/INbrDKqwHlJYiBMPiW5PO9S91028jtu7biN/Fhot6/SB0542ACKx9uhKgRzgQ6C70dth5zaYDeyoC6ZcYBXQ3+jtuIDtv6yZ/99LW2gDm+0X5fEEQggh7KKrHIcUQghhMAkcIYQQdiGBI4QQwi4kcIQQQtiFBI4QQgi7kMARQghhFxI4Qggh7EICRwghhF38f7jSqv9raiZlAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "az.plot_posterior(idata3, var_names=['theta_1'], point_estimate='mode');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So the value of theta that gives us the maximum likelihood is 0.2 (which makes sense because we observed 2 heads in our 10 flips). So we can use that. Of course our null hypothesis has theta fixed at 0.5." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "mle_h0 = .5\n", "mle_h1 = .2 # 20% of flips were heads" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These are the values of theta that maximize the likelihood of the observed data. Now we need to know what the likelihood of our observed data under each of these values of theta. We know how to calculate the likelihood of a set of flips from earlier." ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "def likelihood(theta, n_flips, n_heads):\n", " return (theta**n_heads) * ( (1-theta)**(n_flips - n_heads) )" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "likelihood_h0 = 0.0010\n", "likelihood_h1 = 0.0067\n" ] } ], "source": [ "likelihood_h0 = likelihood(mle_h0, n_heads+n_tails, n_heads)\n", "print(f'likelihood_h0 = {likelihood_h0:.4f}')\n", "likelihood_h1 = likelihood(mle_h1, n_heads+n_tails, n_heads)\n", "print(f'likelihood_h1 = {likelihood_h1:.4f}')" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Likelihood Ratio = 6.87\n" ] } ], "source": [ "print(f'Likelihood Ratio = {likelihood_h1 / likelihood_h0:.2f}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the limit of large data, likelihood ratios (multiplied by 2) are chi-squared distributed, with a degree of freedom equal to the difference in the number of parameters in the two models. Here, our alternative hypothesis has 1 parameter (theta) and our null hypothesis doesn't have any. So the df=1. Let's calculate a p-value." ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "p = 0.0002\n" ] } ], "source": [ "from scipy.stats import chi2\n", "\n", "print(f'p = {1 - chi2.cdf(2 * (likelihood_h1 / likelihood_h0), 1):.4f}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So our likelihood-ratio test is suggesting that we should be extremely skeptical of our null hypothesis whereas the Bayes factor was basically ambivalent. What is going on?\n", "\n", "The key difference between the Bayes factor and the likelihood-ratio test is that the Bayes factor treats our alternative hypothesis as embodying the full prior (i.e., theta~U(0,1)), whereas the likelihood ratio test, being a frequestist test, doesn't know anything about our priors. As a result, the likelihood-ratio test permits the alternative hypothesis to reflect whatever value of theta is most consistent with the observed data (i.e., the maximium likelihood estimate). But that's an extrodinary degree of flexibility. Our alternative hypothesis gets to adapt itself to the data it is seeking to explain, no matter how credible the final estimate was before we observed the data. This means that we should construct our hypotheses so as to be as open-minded and agnostic as possible, because we are only penalized when we observe data that are inconsistent with any configuration of our hypothesis (e.g., combination of parameter values). We are penalized for being unparsimonious, but only coarsely (i.e., the alternative hypothesis is penalized for having 1 more parameter than our \"null\").\n", "\n", "In the Bayes factor, our agnosticism about the cedible values of theta represents a substantial tradeoff. Being uncertain is good because an uncertain hypothesis will be somewhat consistent with many different patterns of data that **might** be observed. However, an uncertain hypotheiss will also be consistent with many different patterns of data that **were not** observed. The former is good, but the latter is bad. The Bayes factor (and all Bayesian approaches) appropriately balance both of these facets and does so thoroughly (incorporating the prior credibility of each parameter value and how the likleihood of the data in light of each parameter value). This is the sense in which people say that Bayesian approaches naturally ensure parsimony. The more agnostic you are (regardless of how many parameters your model has), the less parsimonious your hypotheses are, and the lower the likelihood of the overall model will be.\n", "\n", "To see this in action, let's consider the same 2 hypothesis but evaluate them on a data set that is highly likely under the \"null\". In this data set, 50% of flips come up heads. In a frequentist context, our two hypotheses are indistinguisable. In a Bayesian context, the parsimony of the \"null\" should cause it to win out over the more agnostic alterntive hypothesis." ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "n_heads = 10\n", "n_tails = 10\n", "data3 = np.repeat([1, 0], [n_heads, n_tails])" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "with pm.Model() as model4:\n", " pm1 = pm.Categorical('pm1', [.5, .5])\n", " \n", " theta_0 = 0.5\n", " \n", " theta_1 = pm.Beta('theta_1', 1, 1)\n", " \n", " theta = pm.math.switch(pm.math.eq(pm1, 0), theta_0, theta_1)\n", " \n", " y2 = pm.Bernoulli('y2', theta, observed=data3)" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Multiprocess sampling (2 chains in 2 jobs)\n", "CompoundStep\n", ">BinaryGibbsMetropolis: [pm1]\n", ">NUTS: [theta_1]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [40000/40000 00:31<00:00 Sampling 2 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 2 chains for 10_000 tune and 10_000 draw iterations (20_000 + 20_000 draws total) took 32 seconds.\n" ] } ], "source": [ "with model4:\n", " idata4 = pm.sample(10000, tune=10000)" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [], "source": [ "pm1 = idata4.posterior['pm1'].mean().item() # mean value of model indicator variable" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Posterior: p(model 1|data) = 0.21\n" ] } ], "source": [ "print(f'Posterior: p(model 1|data) = {pm1:.2f}')" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Posterior: p(model 2|data) = 0.79\n" ] } ], "source": [ "print(f'Posterior: p(model 2|data) = {(1-pm1):.2f}')" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Bayes factor: p(model 2|data)/p(model 1|data) * p(model 1)/p(model 2) = 3.65\n" ] } ], "source": [ "print(f'Bayes factor: p(model 2|data)/p(model 1|data) * p(model 1)/p(model 2) = {(1-pm1)/pm1 * (.5/.5):.2f}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using our t-shirt guide to interpreting Bayes factor, we have \"substantial evidence\" in favor of our \"null\" hypothesis. Why? Because our alternative hypothesis was agnostic and implied that many theta values were credible. The \"null\", in contrast, committed to exactly 1. So our null hypothesis is far more parsimonious than the alternative.\n", "\n", "What does the likelihood-ratio test have to say about this?" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEOCAYAAACKDawAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAx8UlEQVR4nO3deVyU5f7/8dfFvgoIiiCIiihq7vu+U9nJ+nosrSw9ph49J83So5WZWv6qY1pWZh7LpeykllqmVpqamlsuqVQC7jsuyCayDnP9/hjkgCwCMszAfJ6Px/3Auea67/tzCbzn5p7rvkdprRFCCFH12Vm6ACGEEBVDAl8IIWyEBL4QQtgICXwhhLAREvhCCGEjJPCFEMJGSOCLSk8p9ahS6sU72noqpbRSqm857aOuUmqGUqp+eWyviH08rJT6Uil1XCllVEptN9e+hG2SwBdVwaPAi3frdI/qAtMBswU+pnG0BPYBF824H2GjHCxdgBAi1yittRFAKbXL0sWIqkeO8EWlppRaBgwDauecwtFKqbN5urgppeYrpeKUUteVUl8opbzv2IaDUuplpVS0UipDKXVZKTVXKeWS83xP4Oec7j/l2U/PnOeHKKW25Ww/RSl1WCk1rLRjuR32QpiLHOGLyu4NoAbQDhiQ05YBeOX8+31gA/Ak0AiYDWRjepG47QvgYeDfwB6gcc526wJ/BX4D/gl8BIwHDuSsdyzna31gNfA2YAS6A58qpVy11gvLbaRC3CMJfFGpaa1PKaWuA5la6323228ffQM7tdbjcv69WSnVCBiplBqutdZKqW7AYGCY1vrznH5blFLxwBdKqZZa6yNKqdvhHpV3Pzk1vJlnv3bAdiAAGAtI4AurIad0RFW38Y7HvwPOgH/O4weATGBNzqkdB6WUA7A55/nud9uBUipMKbVCKXUJyMpZRmL6i0IIqyFH+KKqi7/jcUbOV5ecrzUBJyCliPV9i9u4UsoD+AlIBV4CTmF6ARkLjChDvUKYjQS+sHU3gHSgWxHPX77L+p2AEKCb1jp3Zk3OXwlCWBX5oRRVQQbgWsZ1fwSmAF5a66132QeF7Mct52vW7QallA/wSBnrEcJsJPBFVXAMqK6UGgscxHTEXiJa6+1KqRXAaqXUu8B+TDNt6gL9gSla6+PAccAAjMh5QzcDiME0qycZ+EgpNR1wB14F4vjfTKESUUqFYJptBKZTSUal1KCcxwe01udKsz0h7iSBL6qCT4GOwJuAN3AOGF6K9YcC4zCdc5+KKczPApuAqwBa6xtKqecw/TWwA7AHeuW8YPwfMBfT1MzLmKaCVsd0ZW5p9AKW3tH2dc7XvwHLSrk9IfJR8hGHQghhG2RaphBC2Ag5pSOEmeVcjFXcwZXWWmdXVD3CdskRvhDmt4T/XZBV2FLc7CAhyo2cwxfCzJRSdQG/Yrrc1FrHVFA5woZJ4AshhI2w9nP48mokhBClpwprlHP4QghhIyTwhRDCRkjgCyGEjZDAF0IIGyGBL4QQNkICXwghbIQEvhBC2AgJfCGEsBHWfuGVEBZzNu4Wqw9d5JcT1zlxLYW0rGx83JxoGliNfk38ebRVbaq5OFq6TCFKzNpvrWDVxYmqKTYpjbmbj7P2t4sAtA2pTtPa1XBzsudGSiaHziVw4loKns4OTOjXkGc6heBoL38sC6tS6JW2EvhC5LHuyCWmffsHGQYjQzuG8Pfu9alZzaVAv98vJvHO5hh2Hr9Ok4BqzH+yFfVreFigYiEKJYEvRFGyjZo3Nhxj2Z6ztK7jzXuDWxLi617sOlprNv15hZfX/k6Gwci//9qch1sEVlDFQhRLAl+IwqRkGBi/4jDboq/xbNd6vPxgOA6lOEVzJSmdcSt+48DZBJ7vE8aEvmEoVejvmxAVRQJfiDslp2fx9OL9/HEpiZkDmjK0Y0iZtpNhyOaVtX+w5reLPNwikHcGNcfF0b6cqxWixCTwhcgrOT2LZxbv58/LSSx4qg39mvjf0/a01ny84xSzf4yhXV0fPnmmLd5uTuVUrRClIoEvxG2pmQaGfvorkReTWPBUayKa1iq3bX939DKTvjpKUHVXlg1vTx1ft3LbthAlJIEvBJjeoP378kNsi77Kgqda88B9AeW+j19P32D08kM42isWD2tHi2Dvct+HEMWQD0ARAmDWxmNsibrK9IebmiXsATrU92XN2M64ONozZNE+fjp21Sz7EaI0JPCFTVm+7xxLd5/l2a71GNa5rln31aCmB9/8owth/h78fflBPvr5JEaj/NEqLEdO6Qibcfh8Ao//Zy9dGvixeFg77O0qZupkaqaByasj2RAZS69GNXj38Zb4uMubucKs5By+sF03UjL4y4e7sLdTbBjXtcJnz2itWb7vHG9sOEZNTxfmP9mKVnV8KrQGYVPkHL6wTdlGzfMrj3DjViYLh7axyFRJpRTPdKrL6jGdUQoeW7iXT385jZUfcIkqRgJfVHnzt51k18k4Zj1yH/fV9rJoLS2Cvdk4rhu9w2sya2MUoz4/SGJqpkVrErZDTumIKu3QOdN5+wEtAnlvcEtLl5NLa82yPWd58/so/Ku5sGR4Oxr6e1q6LFF1yDl8YVtSMgz0f/8XjFrz/fPdrPLe9UcuJDLq84OkZ2azYGhruoXVsHRJomqQc/jCtsz47k8uJqQyb3BLqwx7gJbB3qz7Zxdq+7jy7LKDbI2S+frCfCTwRZX04x9XWH3oIs/1akDbutUtXU6xAr1dWTm6I+EBnoz54hA/x1yzdEmiipJTOqLKSUzNpO+7O/Gv5sy3/+xSaT6NKjk9iycW7TN9tOLYzjQOqGbpkkTlJad0hG2YtTGKxNRMZg9qXmnCHqCaiyOLh7XD08WRZ5cdIC4lw9IliSqm8vw2CFECO45fZ/Whi4zpEUrTwIJTMDds2IBSirNnz5q9loyMDCZOnEjNmjVxd3fnoYceuut+a3m58OmwtsTdymTiV0d57733UEoxaNCgAn137dpFp06dcHFxITAwkKlTp2IwGMw0GlEVSOCLKuNWhoFX1v5OaA13nuvdwNLlMH78eJYtW8acOXNYvXo1cXFx9OvXj/T09GLXu6+2F9Measy2wyd4dfpMatQoOHPnzJkz9OvXD39/f7755htefvll3n//fSZNmmSu4YgqwMHSBQhRXt7feoJLiWmsHtPJ4p82dfHiRRYvXsySJUt45plnAGjevDn16tXjiy++YOTIkcWuP7RjCP/vpedJCmlNPc+sAs+//fbbBAQEsHr1ahwcTL/GWmsmTpzIlClTCAgwz11AReUmR/jCooYPH07btm3ZuHEjTZo0wc3NjYceeoj4+HhOnjxJr169cHd3p23btkRGRuaul5qayvjx46lVqxYuLi40b9WG+Z+v4fG2QbmzcrTWzJgxg5o1a+Lp6ckzzzxDcnJygRrS09OZPHkywcHBODs706JFC77//vt7GtfmzZsBGDhwYG5b7dq16dq1Kz/88MNd1z948CAXf9tG8P0jOR2Xwp1zK44cOULPnj1zwx4gIiICg8GQu28h7iSBLyzu/PnzvPbaa8yaNYtFixaxZ88eRo8ezZAhQxgyZAirV6/GYDAwZMiQ3HvPjBo1iqVLlzJ16lTWrl1LovIk9usZdK92I3e7H3zwAa+//jqjR49m9erVuLq6Mnny5AL7HzRoEMuWLeOVV15h/fr1tGvXjgEDBnDkyJHcPkajEYPBUOySnZ2d2z86OpqgoCA8PDzy7atx48ZER0cX+/+htea5555jyuTJvPFkd26mGzgffytfn/T0dJyc8t8TyNnZGYCoqKhity9smNbamhdRxQ0bNkzb29vrkydP5rb961//0oD+7LPPcts2btyoAX3s2DF97NgxrZTSy5Yt01pr/eMfsbrO5O907XoNdUREhNZaa4PBoAMCAvSYMWPy7a9v374a0GfOnNFaa71lyxYN6O3bt+fr161bNz1o0KB8dWKaJlzk0qNHj9z+I0eO1C1atCgw3qlTp+qAgIBi/08WL16s69Spo1NTU7XRaNQ1G7bS1Rp31TdSMnL7DBw4ULdu3TrfeitXrtSAHjVqVLHbFzah0EyVc/jC4urWrUtoaGju4wYNTG+49u7du0DbpUuXuHz5MlprHnvsMdKzsnljwzHCa3nR8eknmDPnHQAuXLhAbGwsjzzySL59DRw4kC1btuQ+3rJlC7Vq1aJLly75Zrj06dOHZcuW5T6eMWMGzz33XLHj8PTMfy8cpQpOhdZaF9p+W1JSEq+88goffPABrq6uANT1defPeM17Px3njUfvA2Ds2LH069ePN954g7Fjx3Ly5Eleeukl7O3tsbe37PsXwnpJ4AuL8/b2zvf49qmKvO2329LT04mNjcXDwwM3Nzc++vkkFxPS+HJUByJ/OkZqaioZGRlcuXIFgJo1a+bb9p2P4+LiuHLlCo6OBW+9kDc469SpQ1BQULHjyBvkPj4+JCYmFuiTmJhYYLx5vfnmmwQHBxMREZG7vqOdpraXI5/v+JMn2wfRONCbvn37MmvWLN544w1ee+01HB0dee211/jggw/w9/cvtk5huyTwRaUTEBBASkoKF64l8PH2U/Rr4k/nUD9++uIqbm5uODs7U6tWLQCuXct/m4I7H1evXp3atWvz7bffFrvPESNG8NlnnxXbp0ePHmzfvh2A8PBwLly4wK1bt3B3d8/tEx0dTXh4eJHbiImJ4eDBg/j4FPLhKPu3M8F5IT++NRqlFFOnTuX555/nzJkzBAUFkZ2dzbRp0+jYsWOxdQrbJYEvKp127dqhlGL8mx+T6tqCKQ80QmvN6tWr6dq1KwDBwcHUqlWLdevW8cADD+Suu3bt2nzb6tOnD3PnzsXDw6PYIC7tKZ2IiAgAvvnmG4YOHQrA5cuX+eWXX1iwYEGR25g1axYTJkzI1zZhwgS8vLxo/9e/s+qMA7tOxuXeVdPDw4NmzZoBMHPmTEJCQujbt2+xdQrbJYEvKp3GjRszYOBjrPt4Fn2fnsDJ3+yZ8sknREdH8/HHHwOm0zGTJ09m0qRJ+Pn50a1bN9asWVNgBku/fv24//776devH1OmTKFp06YkJydz5MgR0tPTeeuttwDT+wx169YtcY1BQUE8++yzTJgwAa01NWrUYMaMGYSEhOS+AAC8/vrrvP7667nvH9x3330FtuXt7Y2fnx+z/jGY3e9sZ87m49QikRUrVtC+fXsMBgMbNmxgyZIlbNy4Md9UTSHykp8MUSn5PzQer3NpHP5uMY8sn0uzZs3YsGFD7hE+mI6M4+PjWbhwIfPmzWPAgAHMnj2bp556KrePUoq1a9fy5ptvMm/ePM6fP0/16tVp2bIl48aNu6caP/jgA9zd3XnxxRdJTU2lR48erFixAhcXl9w+RqMx33TO4jg72DO+Txgvrf2dA+dc2LJlC3PnzsVgMNCuXTu2bt1Kt27d7qlmUbXJ3TJFpXPkQiKPfrSb8b0b8GJEI0uXU6Gyso30e3cHrk4ObBzXFTu7omf8CJsmd8sUlZ/Wmrd/iMLX3YnRPULvvkIV42hvx4S+DYmKTeaHP65YuhxRyUjgi0plz6kb7Dsdz7jeDfBwts0zkg+3CCSspgfv/hRDtlH+CBYlJ4EvKg2tNXM3xxDg5cKQ9nUsXY7F2NspJkY05NT1W3x7+JKlyxGViAS+qDR2HL/Ob+cTea53A4vfDdPS7m9ai6aB1Xh/6wmyso2WLkdUEhL4olLQWvPuT8cJ8nHlsTbBli7H4pQyHeWfj09l9aGLli5HVBIS+KJS2Bp1jciLSYzvHYaTg/zYAvRqVJNWdbz5YOsJ0rNKNrVT2Db5zRFWz2g0Hd3X9XVjYOvali7HaiilmBTRiNikdFbuP2/pckQlIIEvrN7mY1c4FpvM833DcKhEH0peETqH+tKxfnXm/3yKtEw5yhfFk98eYdWMRs17P52gfg13BrSQo/s7mc7lNyIuJYPP9561dDnCykngC6v2wx9XiLl6k+f7hGEvV5UWql3d6nRvWIOFO06RkmG4+wrCZkngC6tlNGre33qcBjU9+EvzQEuXY9Um9mtIQmoWS3edsXQpwopJ4AurtfH3WI5fTZGj+xJoEexNvyb+LPrlNEmpWZYuR1gpCXxhlbKNmve3nqChvwcPNQuwdDmVwov9GnIz3cAnv5y2dCnCSkngC6u0IfIyJ6+l8HyfhnJHyBJqHFCNvzQPYMnuM9xIybB0OcIKSeALq5Nt1Hyw9QThtTx58L5ali6nUpnQtyHpWdks3HHK0qUIKySBL6zO+qOXOXX9Fs/3CZOj+1JqUNODR1vV5vO957ianG7pcoSVkcAXVsWQbcw9ur+/qRzdl8XzfcLINmrmbztp6VKElZHAF1blu6OXOR13iwl95dx9WYX4uvN4u2BW7D/Pqespli5HWBEJfGE1bh/dNwmoxv1N/S1dTqX2Qt+GuDja8/82Rt29s7AZEvjCanxz+BJnb6QyoW8YSsnR/b2o4enMuN4N2BZ9jR3Hr1u6HGElJPCFVUjPyua9n47TPMiLfk3k6L48DO9SlxBfN97YcEw+JEUAEvjCSnyx7xyXk9J56YFwObovJ84O9kzt35iT11JYultuuSAk8IUVSE7PYv7PJ+kW5kfnBn6WLqdK6dfEn76Na/LeTye4EJ9q6XKEhUngC4v7z45TJKZmMeWBcEuXUuUopZj5yH0oBdPW/YHW2tIlCQuSwBcWdS05ncW7zjCgRSD31faydDlVUm1vVyZFNGJ7zHU2RMZauhxhQRL4wqLe2RRDtlEzMaKhpUup0oZ1rkvzIC9mrj9GYmqmpcsRFiKBLyzm8PkEvj50kRFd6xHi627pcqo0ezvFWwObkZiaybR1f1q6HGEhEvjCIoxGzYzv/qSmpzPjeodZuhyb0DTQiwl9w1h/9DLrjlyydDnCAiTwhUWsPnSRoxeTeLl/OB7ODpYux2aM6RFK6zreTPv2Dy4nplm6HFHBJPBFhYu/lcnbP0bTJsSHR1vKB5NXJAd7O959vCUGo2bS10cxGmXWji2RwBcV7vX1f3IzPYu3BjaTi6wsoK6fO9P+0oQ9p27wsdw336ZI4IsK9XPMNb49cpl/9GxAQ39PS5djs4a0C+bhFoHM3RzDnlNxli5HVBBl5RdiWHVxonRupmdx/3s7cXd2YMP4rjg72Fu6JJt2K8PAgPm7SErLYuP4bvhXc7F0SaL8FPqnsxzhiwqhtWbqN39w9WYG/x7UXMLeCrg7O7BwaBtSM7MZvfwQ6VnZli5JmJkEvqgQXx+6yHdHL/NC3zBa1/GxdDkiR5i/J+8+3pLIi4nyJq4NkMAXZnfyWgrT1/1Jp/q+jO3ZwNLliDs8cF8tpjwQzobIWOZtOW7pcoQZyQRoYVZJaVmMXn4QVyd75g1pib18bKFV+nv3+py+nsIH205Sx9edQW2CLF2SMAMJfGE22UbN+BWHOX8jlS9GdpA3Ba2YUopZjzbjUmIaU9ZE4uHswAP3yYfIVzVySkeYzVvfR7Hj+HVef+Q+Otb3tXQ54i6cHOxY9HRbmgd5MX7FYXbKRyNWORL4wiw+/eU0n+46w7BOITzZoY6lyxEl5O7swLLh7alfw53Ryw9y8Gy8pUsS5UgCX5S7NYcuMmtjFP2b1eK1h5tauhxRSl5ujix/tgOBXq78bekBjl5ItHRJopxI4Ity9cPvsUxeE0mXBr68N1jepK2sang688XIDni7OzJ08a8S+lWEBL4oN+uOXOK5FYdpEeTFf55uKxdXVXKB3q6sHN0JbzdT6EdeTLR0SeIeSeCLcvH1wQtMWHWEtiE+fP5sB7nlcRVR29uVFaM64uXqyNBPf+X3i0mWLkncAwl8cc++/PU8/1odSdcGfiz7W3sJ+yomyMeNFaM64uliOtL/45KEfmUlgS/uybLdZ3jlm9/p1agGnzzTFlcnOY1TFQVXd2Pl6I54ODvw1KcS+pWV3C1TlNnH20/x7x+jiWjiz/wnW+PkIMcPVd2F+FSGLNpHSoaB/47swH21vSxdkihcobMlJPBFqWmtmbv5OPN/PsmAFoHMfbwFjvYS9rbidujfyjSFftNACX0rJIEv7p3Wmtc3HGPp7rMMaRfM//u/ZjL10gadv5HKkEV7Sc3KltC3ThL44t5kGzWvrP2dVQcv8GzXerz6UGP5iEIbdv5GKoMX7SUtK5svR3akSWA1S5ck/kcCX5RdVraRF1YdYUNkLOP7hPFC3zAJe8G5G7cYsmgf6VnZ/FdC35pI4IuyyTZqxq88zMbIWF5+MJy/9wi1dEnCipyNM4V+hiGbL0d1pHGAhL4VkI84FKWntebVb39nY2Qsr/SXsBcF1fVzZ+Xojjg72PPUp78SfSXZ0iWJIkjgiyJprXn7h2hW7L/AP3uFMrq7hL0oXF0/d1aM7oijveLJTyT0rZUEvijSgu2n+M/O0zzdMYRJEY0sXY6wcvX83Fk5uhOO9oonFu2TG65ZoSoZ+OlZ2RiyjZYuo1Jbvu8c72yK4dGWgcwc0FTeoBUlUs/PnVWjO+Hu7MCTn+xjz6k4S5ck8qiSb9pOXn2Urw5exMnBDjcne6q7ORFc3Y061d1oUNODFsHeNAmoJleGFmHdkUtMWHWEPuE1+XhoG7moSpTalaR0nlnyK2dvpPLhE624v6l8XGIFs51ZOtuir/LHpWRuZRpIy8wmLiWDC/FpnLtxi+R0A2D6OLeWwd70alSTPo1rElbTQ45iga1RVxm9/BBtQ3z4bER7XBzl3jiibBJTMxm+9ACRFxOZ/nBTnukUIr9jFcd2Ar84sUlpHD6fyOHzCew6eYOoWNObS3V93fhr6yAGtgmitrdree+2Uth76gbDl+6nUS1P/juyA54ujpYuSVRytzIMjF9xmK3R1xjYujZv/l8zOYioGBL4hbmSlM626GusP3qZvadvoBR0DvXlsTbB9G8WYDOnfSIvJvLEon0EeLvy1d87Ud3dydIliSrCaNR8uO0k87Yep3Gtanz4ZCtCa3hYuqyqTgL/bi7Ep7L2t0us+e0i5+NTqeHpzNMdTR/C7efhXJGlVKgTV2/y+H/24u7swOoxnanl5WLpkkQVtC36Ki+sOkpaZjb/6BXK2J6h8qlo5iOBX+Kdas0vJ+JYsvsM22Ou4+RgxyMtAvlbl3pV7tLxczdu8fh/9pJthNVjOlHXz93SJYkq7NrNdGZtiOK7o5cJru7KuF5h/F/r2jIxoPxJ4JfFyWspLNtzhjWHLpGWlU3PRjUY1zuMNiE+li7tnl1KTOPxhXu5lWlg5eiOhNeqWi9mwnrtPH6ddzbF8PulJGp6OjO4XTB/aR5IQ3+ZPFFOJPDvRWJqJv/99TyLd50h/lYmnUN9Gdc7jI71q1fKH9CryekM/s9ebtzK5MuRHWkWJLe3FRVLa832mOt8vvcs249fR2sI8nGlawM/7qvtRZPAagR5u+Lr4VziW3BnZRtJy8omLTNnycrOfexgp/D1cKaGhzPVXB0q5e9tKUjgl4fUTAP/3Xee/+w8TVxKBu3q+jCudxjdwvwqzQ9QbFIaQz/9lStJ6Swf2YHWdSr/XyuicruanM7WqGtsjbrKofMJJKZm5T5nCmonXBztcXaww8nBDq0hw2Ak02AkPU+oG4wli4xqLg40CaxG00AvWgZ70znUF9+q9T6dBH55Ss/KZuX+8yzccZoryem0CPbm+T4N6NWoplUH/6nrKTyzeD/JaVksHt6O9vWqW7okIfLRWnMxIY3oKze5kpRGbFI6cSkZuQGfYTCiAGdHO5zs7XBxtMfVyR5Xx5zF6X+P3ZzsTc872pOVrYlLySAuJYPTcbf483Iy0bHJZBhMV+U3DqhGl1BfujTwo3296rg7O1j2P+LeSOCbQ4YhmzWHLrFg+0kuJqTRPMiLCX3DrDL4D51LYPTnBwH4bER7+TxSYfMM2UZ+v5TEnlM32H0yjoPnEsg0GHGwU7Su40OXBn50aeBLi2DvyvbGsgS+OWVlG1n720U+3GZ9wa+15r+/nmfm+j8J8HJl2d/aUV/mQQtRQHpWNgfPJrDrZBx7TsXx+6UktAZ3J3s61velTV0fWgZ50yzI654vTNRac/1mBsevpnD86k1OXU/hanI6129mkJxu4P0hLWke5F3WzUvgVwRrC/5rN9N5Ze0fbIm6Ss9GNXh/cCu83OQKWiFKIjE1k72nbrD7VBx7Tt7gdNwtAJSC0BoeNKrlSV1fN+r6uhPi6051d0e83ZzwdHHATim0hpvpWcTfyuRyUjoncoL9xNUUTlxLISntf+9VeLk6EuDlQs1qLni5OjK+dwPC/D3LWroEfkW6M/gbB1RjRJe6PNwisEIuLc80GPl871k+2HqCDIORf93fiL91qScfOC7EPUhMzSTyYhJHLyRy9GIiJ6+lcCEhjewSvlkM4OPmSFhNT0JretDI34OG/p6E+Xvi5+FUngeFEviWkJVt5JvDl1j8yxlirt7Ez8OJQW2CebRVoFnmvcffymTF/vN8se8csUnp9GhYg9cebiKXsgthJlnZRi4lpHE+PpWE1EyS0rJIznPk7uHsQHUPZ2p6OhNW06OiZgNJ4FuS1po9p26wdPcZfo65TrZRE17Lk4gm/nQM9aV1HZ8yHfkbso1Exd7k4Ll4tkVfY9/pG2Rla7qF+TGqW326N6xhhtEIIaycBL61iEvJYGNkLN8dvczh8wkYNdjbKUJruNOoVjUCvFzwdXfCw8UBRzs7HOwVhmzNzQwDKekGktKyuJCQyvkbqZyLv0V6lmlaWT0/dyKa+DOoTdC9nPsTQlR+EvjWKDk9iwNn4vntfAIxV24Sc/Um15IzcucGF8bV0Z4gH1dCfN0I8XWnZbA3bUJ8CLTR2zoLIQqQwK8stNakZBi4lZFNVrYRg1HjaK/wdHbE3dkeh8o1H1gIUfEk8IUQwkYUGvhyqCiEEDZCAl8IIWyEBL4QQtgICXwhhLAREvhCCGEjJPCFEMJGSOALIYSNkMAXQggbIYEvhBA2QgJfCCFshAS+EELYCAl8IYSwERL4QghhIyTwhRDCRkjgC6vy7bff0rx5c5ydnalXrx7vvvtusf0nTJiAUopJkybla4+OjqZDhw54eXkxZMgQUlJS8j2/c+dOateuXaC9MMuWLUMpVWjfGTNm4Ofnl/v47NmzKKVyF3d3d0JDQ3nqqaf45ZdfCqw/fPhw2rZte9cahCgPEvjCauzevZuBAwfSvn171q9fz4gRI5gyZQrz5s0rtP+xY8dYsmQJ1aoV/DD44cOH06BBA7766iuOHTvGm2++mfuc0WhkwoQJvPXWW3h4mOfD3efMmcPevXv5/vvvmTZtGjdu3KB79+7MnDnTLPsTokS01ta8CBsSERGhu3Xrlq/thRde0D4+PjojI6NA/z59+uhXX31Vh4SE6IkTJ+a237x5UwP62rVrWmutV65cqdu2bZv7/KJFi3T79u210WgsUV1Lly7VgL5582aB56ZPn659fX1zH585c0YDev369QX6Tps2TQP6559/zm0bNmyYbtOmTYnqEKIUCs1UOcIXVuPIkSP07ds3X1tERAQJCQns3bs3X/vq1auJioripZdeKrCdzMxMAFxdTZ/x6+bmltuWnJzMtGnTeP/991Gq0A8FMpvp06cTGBjIwoULK3S/QtwmgS+sRnp6Ok5OTvnanJ2dAYiKisptS0tLY+LEibz99tu4u7sX2E716tWpV68eH374IfHx8SxatCj3PPkbb7xB37596dixY6nry87OxmAw5FuMxqI/bP5O9vb29O7dm3379pV630KUBwdLFyDEbQ0aNODAgQP52vbv3w9AfHx8bttbb71FQEAAQ4cOLXJbH330EY899hivvPIKYWFhfPTRR5w8eZLFixcTGRlZpvq8vb0Lbff19S3xNoKCgrh69WqZ9i/EvZIjfGE1xowZw7p16/jkk09ISEhg06ZNzJ07FzAdHQOcOXOGOXPmMG/evGJPyTz44INcu3aNmJgYoqKiqFOnDi+++CIvvPACQUFBfPTRR9SpU4c6deqwYMGCEtW3c+dODhw4kG8ZNWpUqcaotS5VfyHKkxzhC6sxYsQIjh49ytixYxk9ejRubm78+9//Zty4cfj7+wPw0ksv8eCDDxIeHk5iYiJgmnWTkZFBYmIiXl5euS8Ebm5uNGzYEIAtW7Zw9OhRVq1axdGjR5k2bRp79uwBoFOnTnTt2pXmzZsXW1+rVq0KzOrZsGFDqcZ46dKl3LEIUdHkCF9YDXt7e+bPn8/169eJjIzk6tWruefab3+NiYlh7dq1+Pj45C4XLlxg/vz5+Pj4cOnSpQLbNRgMTJgwgdmzZ+Pq6sr27dvp3bs34eHhhIeH06dPH3bs2GH28RkMBrZt20anTp3Mvi8hCiNH+MLq3A5ygAULFtC5c2fCw8MB+PTTTwtcADVkyBB69OjB2LFjqVGjRoHtLVy4EB8fHwYPHpzblpqamvvvW7duVcipltdff53Lly8zZswYs+9LiMJI4AursW/fPnbt2kXLli1JTk5mxYoVbNq0iV27duX2KeyqVBcXF4KDg+nZs2eB5xISEpg5cyabNm3KbevevTuTJ09myZIlAGzbto233367XMcSExODn58fmZmZnDlzhpUrV/Ljjz8yY8YMevToUa77EqKkJPCF1XB0dGTVqlXMmDEDOzs7unXrxu7du2nWrFmZtzl9+nQGDBhA69atc9tatWrF7NmzmTp1KmC6KrZFixb3XH9et2/14OLiQkBAAJ06dWLnzp1069atXPcjRGkoK581YNXFCSGElSp0Cpu8aSuEEDZCAl8IIWyEBL4QQtgICXwhhLAREvhCCGEjJPCFEMJGSOALIYSNkMAXQggbIYEvhBA2QgJfCCFshAS+EELYCAl8IYSwERL4QghhIyTwhRDCRkjgCyGEjZDAF0IIGyGBL4QQNkICXwghbIQEvhBC2AgJfCGEsBES+EIIYSMk8IUQwkZI4AshhI2QwBdCCBshgS+EEDZCAl8IIWyEBL4QQtgICXwhhLAREvhCCGEjJPCFEMJGSOALIYSNkMAXQggbIYEvhBA2QgJfCCFshAS+EELYiBIFvlKqiVJqq1IqVSl1WSn1ulLKvgTreSmlliqlEpRSSUqp/yqlfAvp94hS6nelVLpS6phSanDe5//8808eeOABAgMDcXZ2pk6dOowcOZLY2NiSj1QIISrIunXraNasGS4uLjRp0oRVq1aVaL1vv/2W5s2b4+zsTL169Xj33XfzPR8bG8u//vUvWrRogYeHB8HBwQwbNozLly+XaPt3DXyllA+wBdDAI8DrwERgZgm2vwroCYwEhgPtgG/v2H5XYA3wM/AgsBFYoZSKuN0nKSmJevXqMWfOHDZt2sTMmTPZsmUL/fv3x2AwlKAMIYSoGLt27eKvf/0rvXr14ocffuChhx7iiSeeYPPmzcWut3v3bgYOHEj79u1Zv349I0aMYMqUKcybNy+3z6FDh/jmm2944oknWL9+Pe+88w6//vornTt3JiUl5e7Faa2LXYCXgQSgWp62yUBq3rZC1uuE6UWie5629jltffO0bQK23bHu98AuXYzNmzdrQB86dKi4bkIIUaEiIiJ0r1698rU9+OCDukuXLnddr1u3bvnaXnjhBe3j46MzMjK01lonJCTorKysfH1iYmI0oJctW5a3udBcLskpnQeBTVrr5DxtKwFXoMdd1ruqtd6Z58VlP3Am5zmUUs5AL+CrO9ZdCXRKSkoqcuO+vqYzQ5mZmSUYghBCmF9GRgY///wzjz/+eL72IUOGsHfvXorLtCNHjtC3b998bRERESQkJLB3714AvL29cXBwyNenYcOGuLm5ce3atbvWV5LADwei8zZorc9jOsIPL816OaLyrBcKOBbSLwqwO378eL5Go9FIZmYmMTExvPTSS7Rr14727duXYAhCCGF+p06dIisri/Dw/NHYuHFjjEYjd2ZaXunp6Tg5OeVrc3Z2BiAqKqrI9SIjI0lNTaVJkyZ3ra8kge8DJBbSnpDz3L2sd/vrnf0SABISEvI19u/fH2dnZ8LDw4mPj2fDhg3Y2clEIyGEdbidWd7e3vnafXx88j1fmAYNGnDgwIF8bfv37wcgPj6+0HWMRiPPP/88YWFhREREFNonr5KmpS6kTRXRXpb17nysAJRS+Ro//PBD9u3bx/Lly0lJSeHBBx8kPT39LiUIIUTFujO7tOm9yQLteY0ZM4Z169bxySefkJCQwKZNm5g7dy4A9vaFT4p8+eWX2bt3L8uXL8fR0fGudZUk8BMA70LavSj8CP5u63nnWS8hT9udfQq8SoaFhdGhQweGDh3Kpk2bOHz4MF9++WUxJQghRMW5fSSfmJiYr/324zszLa8RI0YwZswYxo4dS/Xq1Rk4cCCvvfYaAP7+/gX6L1iwgHfeeYfPPvuMDh06lKi+kgR+NHecq1dKBQPuFH6Ovsj1cuQ9t38KyCqkXzhgbNiwYZEbDwkJoXr16pw+fbrY4oUQoqKEhobi6OhIdHT+aIyOjsbOzo7iMs3e3p758+dz/fp1IiMjuXr1Kh07dgTI/XrbmjVrGDduHLNnz2bw4MGFba5QJQn8H4D7lVKeedoGA2nAjrusVytnnj0ASqm2QP2c59BaZ2Caf//YHesOBvZ6eXkVufGYmBhu3LhBvXr1SjAEIYQwP2dnZ3r16sXXX3+dr33VqlV06tSJ4jLtNh8fH5o1a4aHhwcLFiygc+fO+d4E3r59O0899RTPPfcckyZNKlV96va5pSI7mC68Ogb8AfwbU2C/C8zTWr+ap99JYIfW+tk8bT8CDYFJgDFn/Wta6255+nQFtgPzMV2U1T+n/wNa600AkyZNwsHBgQ4dOuDt7U1UVBSzZ8/GwcGBo0eP4u7uXqpBCyGEuezatYuePXvy3HPP8eijj/L9998zZ84cfvzxx9w3Vs+dO0doaChLlizhmWeeAWDfvn3s2rWLli1bkpyczIoVK9i0aRO7du2iefPmgGm2TqdOnahbty4LFizIN2mlRo0ahIaG3n5Y+JsFRU3Q1/kvhGoCbMN0VB8LvAHY39HnLLDsjjZvYCmmc/bJwJeAXyHbfxTTC0oGptM9Q3Ke01prvWLFCt25c2ft4+OjXV1ddaNGjfSLL76or1+/XuyFDEIIYQnffPONbtq0qXZyctKNGjXSK1asyPf8mTNnNKCXLl2a23bw4EHdtm1b7e7urj09PXX//v11ZGRkvvWWLl2qMU1yKbAMGzYsb9dCs/yuR/gWZtXFCSGElSr0CF8msQshhI2QwBdCCBshgS+EEDZCAl8IIWyEBL4QQtgICXwhhLAREvhCCGEjJPCFEMJWFHVFlpUs92T69OlFXpUmiyyyyHIvy/Tp0+81osypzB9xKIQQogqQwBdCCBsh99IRQoiqR+6lI4QQtkwCXwghbIQEvhBC2AgJfCGEsBES+EIIYSMk8IUQwkZI4AshhI2QwBdCCBshgS+EEDZCAl8IIWyEBL4QQtgIq76XzsyZM38E/Mq4eiBwuRzLqQxkzLZBxmwb7mXMcdOnT3+gQGtR902u7MuMGTOKvCd0VV1kzLaxyJhtYzHHmOWUjhBC2IiqHPgzLV2ABciYbYOM2TaU+5it+hy+EEKI8lOVj/CFEELkIYEvhBA2QgJfCCFsRKUNfKXUP5RSZ5RS6UqpQ0qpbnfp30wptUMplaaUuqSUek0pVejnPlqr0oxZKdVTKbVOKRWrlEpVSkUqpUZUZL3lobTf5zzrhSmlbiqlUsxdY3krw8+2UkpNUEpFK6Uycr7nb1dUveWhDGO+Xym1N+d7HJfzs96wouq9F0qp7kqp73JySCulhpdgnfLJL0vPNS3LAgwGsoBRQGPgQyAFqFNE/2rAFeAr4D7gr8BNYKKlx2LGMb8CzAK6APWBsYABeNLSYzHXmPOs5wQcAjYCKZYeh7nHDLwLHAceyfletwL6W3os5hozUA9IB2YDDYCWwGbgpKXHUsLx9gfeBAYBqcDwu/Qvt/yy+ODL+B/2K/DJHW0ngLeK6D8WSAZc87S9ClwiZ6aStS+lHXMR2/gKWGPpsZh7zMB7wFJgeCUM/NL+bDfKCcvGlq69Asc8CMgG7PO09QI04Gfp8ZRy7CklCPxyy69Kd0pHKeUEtMH0ip7XZqBzEat1An7RWqfladuE6dLluuVdY3kr45gLUw1IKK+6zKmsY1ZKPQT8BRhvvurMo4xjfgQ4DTyglDqtlDqrlPpMKVXTjKWWmzKO+SCmF7mRSil7pZQnMAw4oLWOM1uxllNu+VXpAh/TvXXsgat3tF8FahWxTq0i+t9+ztqVZcz5KKX+AvQBFpVvaWZT6jErpQKAT4CntdY3zVueWZTl+1wfCAGGYPqL5mkgHFivlKoMv9+lHrPW+izQD9OFSRlAEtAM0wt9VVRu+VUZfiCKcucVY6qQtrv1L6zdmpV2zKZOSnUBvgTGa633m6MwMyrNmL8APtZa7zNvSWZXmjHbAc6YXuR2aq1/wRT67YF25iux3JV4zEqpWsBi4HNMY+yJ6Zz2V5XkRa4syiW/KuN/Thym83d3vrLVpOCr4G1XiuhPMetYk7KMGQClVFfgB+A1rfXH5inPLMoy5t7AdKWUQSllwBQK7jmPR5uv1HJTljHHAgat9fE8bScwvUFfp9wrLH9lGfM/gVta68la68Na653AUKAHpTvFWVmUW35VusDXWmdimoHR746n+gF7ilhtL9BNKeVyR//LwNnyrrG8lXHMKKW6Ywr7mVrreWYr0AzKOOZmmGZs3F5eA9Jy/v11+VdZvso45t2Ag1IqNE9bfcABOFfuRZazMo7ZDdOLRF63H1e6TCuB8ssvS79LXcZ3tgcDmcBITNO43sf0bndIzvNvAVvz9PfC9Cq5EtO0poGY3vWubNMySzPmnsAt4B1MRwe3lxqWHou5xlzI+sOpfLN0Svt9tsMUmDswTcdslfPvfYCdpcdjpjH3BozAdCAMaA38CJwH3C09nhKM14P/HZSkYjowaUnONFRz5pfFB38P/2n/wPTqlpHzA989z3PLgLN39G8G7MQ0fzc254elUkzJLMuYcx7rQpazFV13RX6f71i30gV+WcYMBGD6C+YmcA34L+Bv6XGYecxDgN9yXhiuA+uBJpYeRwnH2rOI381lxYy3XPJL7pYphBA2oiqe7xJCCFEICXwhhLAREvhCCGEjJPCFEMJGSOALIYSNkMAXQggbIYEvhBA2QgJfCCFshAS+EELYiP8P40/tI0YOp60AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "az.plot_posterior(idata4, var_names=['theta_1'], point_estimate='mode');" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [], "source": [ "mle_h0 = .5\n", "mle_h1 = .5 # 50% of flips were heads" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "likelihood_h0 = 0.00000095\n", "likelihood_h1 = 0.00000095\n" ] } ], "source": [ "likelihood_h0 = likelihood(mle_h0, n_heads+n_tails, n_heads)\n", "print(f'likelihood_h0 = {likelihood_h0:.8f}')\n", "likelihood_h1 = likelihood(mle_h1, n_heads+n_tails, n_heads)\n", "print(f'likelihood_h1 = {likelihood_h1:.8f}')" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Likelihood Ratio = 1.00\n" ] } ], "source": [ "print(f'Likelihood Ratio = {likelihood_h1 / likelihood_h0:.2f}')" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "p = 1.0000\n" ] } ], "source": [ "print(f'p = {1 - chi2.cdf(-2 * (likelihood_h1 / likelihood_h0), 1):.4f}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result of this likelihood-ratio test is pretty trivial, but confirms the expectation described above. By treating each hypothesis as synomymous with the corresponding maximum likelihood estimate of theta, the two hypotheses end up being identical when we observe heads on 50% of our flips.\n", "\n", "### Take home message\n", "Bayes factors are fine, but I would almost never recommend them. They are useful for performing NHST-style \"tests\" in a Bayesian framework, but the idea that you have to **choose** between two or more hypotheses is something that (I think) researchers never should have been doing in the first place.\n", "\n", "This is particularly true in the case of these point-estimate, \"null\"-style models. We either believe that the value of $\\theta$ is *exactly* 0.50000000 or whether have no clue whatsoever what the value of $\\theta$ is? It seems plausible that we don't actually believe one of these (e.g., as is the case in most NHST settings).\n", "\n", "On top of that, Bayes factors do not reflect our prior beliefs about the credibility of the models we are comparing. Bayes factors speak to the \"evidential value\" of our data. But, as we saw above, the data can strongly imply one model and our priors can strongly favor the other. In such cases, Bayes factors only provide part of the relevant story.\n", "\n", "What do I recommend? If you have data and a model and you would like to answer questions about the values of model parameters in light of the data, then **estimate** the credible values of those parameters and use the posterior to answer your questions. If you are not sure whether a coin is fair or not, build a model in which theta can take on many values of theta and ask how credible the values far from 0.5 are (for one or more definitions of \"far\"). Ask how credible values of theta close to 0.5 are (for one or more definitions of \"close\"). But dichotomizing the world seems unwise. Also the sampler doesn't like it." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.5" } }, "nbformat": 4, "nbformat_minor": 4 }