{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "# Statistical Hypothesis Testing\n", "\n", "![](./img/stats/preface2.png)\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "The statistics of hypothesis can be thought of as observations of random variables from known distributions, \n", "- which allows us to make statements about how likely those assumptions are to hold." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "- H1: data scientists prefer Python to R. \n", "- H2: people tends to jump onto the bandwagon of collective gatekeepers.\n", "- H3: media agenda determines the public agenda. \n", "- H4: the rich people can better use mobile phone compared with the poor. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## The logic of falsification \n", "\n", "证伪的逻辑\n", "\n", "- **Null hypothesis H0** represents some default position \n", "- **Alternative hypothesis H1** \n", " \n", " - We use statistics to decide whether we can **reject H0** as false or not." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "We do not prove (just support) alternative hypothesis;\n", "\n", "\n", "\n", "We just support alternative hypothesis when we reject the null hypothesis." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "## Example: Flipping a Coin\n", "Imagine we have a coin and we want to test whether it’s fair. \n", "\n", "**Null hypothesis H0**: The coin is fair. \n", "- The probability of landing heads is 0.5\n", "- The alternative hypothesis p ≠ 0.5" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Each coin flip is a Bernoulli trial, \n", "- X is a Binomial(n,p) random variable\n", "- we can approximate X using the normal distribution:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "ExecuteTime": { "end_time": "2018-10-27T01:14:27.896305Z", "start_time": "2018-10-27T01:14:27.891124Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "from probability import normal_cdf, inverse_normal_cdf\n", "import math, random" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "ExecuteTime": { "end_time": "2018-10-16T06:23:22.736087Z", "start_time": "2018-10-16T06:23:22.731811Z" }, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "def normal_approximation_to_binomial(n, p):\n", " \"\"\"finds mu and sigma corresponding to a Binomial(n, p)\"\"\"\n", " mu = p * n\n", " sigma = math.sqrt(p * (1 - p) * n)\n", " return mu, sigma" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Whenever a random variable follows a normal distribution, we can use **normal_cdf** to figure out the probability that its realized value lies within (or outside) a particular interval:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "ExecuteTime": { "end_time": "2018-10-27T01:15:44.324145Z", "start_time": "2018-10-27T01:15:44.321242Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "# the normal cdf _is_ the probability the variable is below a threshold\n", "normal_probability_below = normal_cdf" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "ExecuteTime": { "end_time": "2018-10-27T01:15:06.947081Z", "start_time": "2018-10-27T01:15:06.937135Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "# it's above the threshold if it's not below the threshold\n", "def normal_probability_above(lo, mu=0, sigma=1):\n", " return 1 - normal_cdf(lo, mu, sigma)\n", "\n", "# it's between if it's less than hi, but not less than lo\n", "def normal_probability_between(lo, hi, mu=0, sigma=1):\n", " return normal_cdf(hi, mu, sigma) - normal_cdf(lo, mu, sigma)\n", "\n", "# it's outside if it's not between\n", "def normal_probability_outside(lo, hi, mu=0, sigma=1):\n", " return 1 - normal_probability_between(lo, hi, mu, sigma)\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Lower and Upper Bunds 上、下边界\n", "We can also find either the nontail region or the (symmetric) interval around the mean that accounts for a certain level of likelihood. \n", "\n", "For example, if we want to find an interval centered at the mean and containing 60% probability\n", "- we find the cutoffs where the upper and lower tails each contain 20% of the probability (leaving 60%):" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "ExecuteTime": { "end_time": "2018-10-27T01:15:11.216282Z", "start_time": "2018-10-27T01:15:11.210067Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "def normal_upper_bound(probability, mu=0, sigma=1):\n", " \"\"\"returns the z for which P(Z <= z) = probability\"\"\"\n", " return inverse_normal_cdf(probability, mu, sigma)\n", "\n", "def normal_lower_bound(probability, mu=0, sigma=1):\n", " \"\"\"returns the z for which P(Z >= z) = probability\"\"\"\n", " return inverse_normal_cdf(1 - probability, mu, sigma)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "ExecuteTime": { "end_time": "2018-10-27T01:15:12.344145Z", "start_time": "2018-10-27T01:15:12.335066Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "def normal_two_sided_bounds(probability, mu=0, sigma=1):\n", " \"\"\"returns the symmetric (about the mean) bounds\n", " that contain the specified probability\"\"\"\n", " tail_probability = (1 - probability) / 2\n", "\n", " # upper bound should have tail_probability above it\n", " upper_bound = normal_lower_bound(tail_probability, mu, sigma)\n", "\n", " # lower bound should have tail_probability below it\n", " lower_bound = normal_upper_bound(tail_probability, mu, sigma)\n", "\n", " return lower_bound, upper_bound" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "In particular, when we choose to flip the coin n = 1000 times. \n", "If our hypothesis of fairness is true, \n", "- X should be distributed approximately normally \n", " - with mean 500 and\n", " - standard deviation 15.8" ] }, { "cell_type": "code", "execution_count": 66, "metadata": { "ExecuteTime": { "end_time": "2018-11-06T13:17:24.594734Z", "start_time": "2018-11-06T13:17:24.589891Z" }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "15.811388300841896" ] }, "execution_count": 66, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import numpy as np\n", "np.sqrt(1000*0.5*(1-0.5))" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "ExecuteTime": { "end_time": "2018-10-17T05:02:15.761379Z", "start_time": "2018-10-17T05:02:15.756894Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "mu_0 500.0\n", "sigma_0 15.811388300841896\n" ] } ], "source": [ "mu_0, sigma_0 = normal_approximation_to_binomial(1000, 0.5)\n", "print(\"mu_0\", mu_0)\n", "print(\"sigma_0\", sigma_0)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Type 1 Error and Type 2 Error\n", "\n", "对于原假设H0，容易犯的两类错误：Type 1“弃真”或 Type 2“纳伪”？\n", "\n", "We need to make a decision about **significance**:\n", "How willing we are to make a **type 1 error** (“false positive”), in which we reject H0 even though it’s true. " ] }, { "cell_type": "markdown", "metadata": { "ExecuteTime": { "end_time": "2018-10-17T05:21:30.174120Z", "start_time": "2018-10-17T05:21:30.168100Z" }, "slideshow": { "slide_type": "subslide" } }, "source": [ "Assuming **p** really equals 0.5 (i.e., H0 is true), there is just a 5% chance we observe an X that lies outside this interval, which is the exact significance we wanted. \n", "\n", "if H0 is true, approximately 19 times out of 20, this test will give the correct result." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "ExecuteTime": { "end_time": "2018-10-17T05:14:11.538421Z", "start_time": "2018-10-17T05:14:11.532364Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "normal_two_sided_bounds(0.95, mu_0, sigma_0)\n", "power of a test\n", "95% bounds based on assumption p is 0.5\n", "lo 469.01026640487555\n", "hi 530.9897335951244\n" ] } ], "source": [ "print(\"normal_two_sided_bounds(0.95, mu_0, sigma_0)\")\n", "lo, hi = normal_two_sided_bounds(0.95, mu_0, sigma_0)\n", "print(\"power of a test\")\n", "print(\"95% bounds based on assumption p is 0.5\")\n", "print(\"lo\", lo)\n", "print(\"hi\", hi)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Power of a test\n", "- the probability of not making a type 2 error, \n", "- in which we fail to reject H0 even though it’s false. \n", "\n", "Let’s check what happens if p is really 0.55:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "ExecuteTime": { "end_time": "2018-10-17T05:40:39.243326Z", "start_time": "2018-10-17T05:40:39.237541Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "actual mu and sigma based on p = 0.55\n", "mu_1 550.0\n", "sigma_1 15.732132722552274\n", "Now the coin is slightly biased toward heads.\n" ] } ], "source": [ "print(\"actual mu and sigma based on p = 0.55\")\n", "mu_1, sigma_1 = normal_approximation_to_binomial(1000, 0.55)\n", "print(\"mu_1\", mu_1)\n", "print(\"sigma_1\", sigma_1)\n", "print(\"Now, the coin is slightly biased toward heads.\")" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "ExecuteTime": { "end_time": "2018-10-17T05:46:40.365634Z", "start_time": "2018-10-17T05:46:40.359618Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "type 2 probability 0.0636203880618903\n", "power 0.9363796119381097\n" ] } ], "source": [ "# a type 2 error means we fail to reject the null hypothesis\n", "# which will happen when X is still in our original interval\n", "type_2_probability = normal_probability_between(lo, hi, mu_1, sigma_1)\n", "power = 1 - type_2_probability # 0.887\n", "\n", "print(\"type 2 probability\", type_2_probability)\n", "print(\"power\", power)\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### One-sided test\n", "\n", "Assume that:\n", "\n", "**null hypothesis** H0 was that `the coin is not biased toward heads, or that p ≤ 0.5.`\n", "\n", "In that case we want a one-sided test that rejects the null hypothesis when X is much larger than 500 but not when X is smaller than 500. " ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "ExecuteTime": { "end_time": "2018-10-17T05:46:45.668057Z", "start_time": "2018-10-17T05:46:45.660493Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "one-sided test\n", "hi 526.0073585242053\n", "type 2 probability 0.06362051966928273\n", "power 0.9363794803307173\n" ] } ], "source": [ "print(\"one-sided test\")\n", "hi = normal_upper_bound(0.95, mu_0, sigma_0)\n", "\n", "print(\"hi\", hi) # is 526 (< 531, since we need more probability in the upper tail)\n", "type_2_probability = normal_probability_below(hi, mu_1, sigma_1)\n", "power = 1 - type_2_probability # = 0.936\n", "print(\"type 2 probability\", type_2_probability)\n", "print(\"power\", power)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "This is a **more powerful** test, since \n", "- it no longer rejects H0 when X is below 469\n", " - which is very unlikely to happen if H1 is true\n", "- it instead rejects H0 when X is between 526 and 531 \n", " - which is somewhat likely to happen if H1 is true" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "ExecuteTime": { "end_time": "2018-10-27T01:15:20.318371Z", "start_time": "2018-10-27T01:15:20.312611Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "def two_sided_p_value(x, mu=0, sigma=1):\n", " if x >= mu:\n", " # if x is greater than the mean, the tail is above x\n", " return 2 * normal_probability_above(x, mu, sigma)\n", " else:\n", " # if x is less than the mean, the tail is below x\n", " return 2 * normal_probability_below(x, mu, sigma)" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "ExecuteTime": { "end_time": "2018-10-19T01:33:51.658934Z", "start_time": "2018-10-19T01:33:51.655300Z" }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "two_sided_p_value(529.5, mu_0, sigma_0) is 0.06207721579598857\n" ] } ], "source": [ "# If we were to see 530 heads, we would compute:\n", "print(\"two_sided_p_value(529.5, mu_0, sigma_0) is \", \\\n", " two_sided_p_value(529.5, mu_0, sigma_0))\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "**Why did we use 529.5 instead of 530? **\n", "- This is what’s called a [continuity correction](https://en.wikipedia.org/wiki/Continuity_correction#Binomial). 连续性校正\n", "- `normal_probability_between(529.5, 530.5, mu_0, sigma_0)` is a better estimate of the probability of seeing 530 heads than `normal_probabil ity_between(530, 531, mu_0, sigma_0)` is." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "One way to convince yourself that this is a sensible estimate is with a simulation:" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "ExecuteTime": { "end_time": "2018-10-19T01:41:54.080802Z", "start_time": "2018-10-19T01:41:54.070102Z" }, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "def count_extreme_values():\n", " extreme_value_count = 0\n", " for _ in range(100000):\n", " num_heads = sum(1 if random.random() < 0.5 else 0 # count # of heads\n", " for _ in range(1000)) # in 1000 flips\n", " if num_heads >= 530 or num_heads <= 470: # and count how often\n", " extreme_value_count += 1 # the # is 'extreme'\n", "\n", " return extreme_value_count / 100000" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "ExecuteTime": { "end_time": "2018-10-19T01:43:23.343104Z", "start_time": "2018-10-19T01:42:52.258488Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/plain": [ "0.06118" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "count_extreme_values() # 0.062" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Since the p-value is greater than our 5% significance, we don’t reject the null. " ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "ExecuteTime": { "end_time": "2018-10-19T01:40:03.507796Z", "start_time": "2018-10-19T01:40:03.504084Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "two_sided_p_value(531.5, mu_0, sigma_0) is 0.046345287837786575\n" ] } ], "source": [ "print(\"two_sided_p_value(531.5, mu_0, sigma_0) is \",\\\n", " two_sided_p_value(531.5, mu_0, sigma_0))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "It is smaller than the 5% significance, and we would reject the null." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "For our one-sided test, if we saw 525 heads we would compute:" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "ExecuteTime": { "end_time": "2018-10-19T01:45:43.269071Z", "start_time": "2018-10-19T01:45:43.266010Z" }, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "upper_p_value = normal_probability_above\n", "lower_p_value = normal_probability_below" ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "ExecuteTime": { "end_time": "2018-10-19T01:47:12.353848Z", "start_time": "2018-10-19T01:47:12.348686Z" }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "0.06062885772582083" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "upper_p_value(524.5, mu_0, sigma_0) # 0.061\n", "# wouldn’t reject the null." ] }, { "cell_type": "markdown", "metadata": { "ExecuteTime": { "end_time": "2018-10-19T01:47:41.824017Z", "start_time": "2018-10-19T01:47:41.819950Z" }, "slideshow": { "slide_type": "subslide" } }, "source": [ "For our one-sided test, if we saw 527 heads we would compute:" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "ExecuteTime": { "end_time": "2018-10-19T01:46:44.057479Z", "start_time": "2018-10-19T01:46:44.053492Z" }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "0.04686839508859242" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "upper_p_value(526.5, mu_0, sigma_0) # 0.047\n", "# we would reject the null." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "\n", "\n", "\n", "Make sure your **data is roughly normally distributed** before using normal_probability_above to compute p-values.\n", "\n", "- There are various statistical tests for normality, \n", "- but even plotting the data is a good start.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Confidence Intervals\n", "\n", "> to construct a confidence interval around the observed value of the parameter.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "**Example**\n", "\n", "We can estimate the probability of the unfair coin by looking at the average value of the Bernoulli variables corresponding to each flip\n", "- 1 if heads, \n", "- 0 if tails. \n", "\n", "If we observe 525 heads out of 1,000 flips, then we estimate `p equals 0.525`." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "How confident can we be about this estimate? " ] }, { "cell_type": "markdown", "metadata": { "ExecuteTime": { "end_time": "2018-10-24T05:16:14.494393Z", "start_time": "2018-10-24T05:16:14.482809Z" }, "slideshow": { "slide_type": "subslide" } }, "source": [ "\n", "If we knew the exact value of p, the `central limit theorem` tells us \n", "> the average of those Bernoulli variables should be approximately normal, with mean p and standard deviation $\\sigma = \\sqrt{p(1 - p) / 1000}$:" ] }, { "cell_type": "code", "execution_count": 69, "metadata": { "ExecuteTime": { "end_time": "2018-11-06T13:32:45.994391Z", "start_time": "2018-11-06T13:32:45.989718Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.525 0.015791611697353755\n" ] } ], "source": [ "p_hat = 525 / 1000\n", "mu = p_hat\n", "sigma = math.sqrt(p_hat * (1 - p_hat)/ 1000) # 0.0158\n", "\n", "print(mu, sigma)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "we are “95% confident” that the following interval contains the true parameter p:" ] }, { "cell_type": "code", "execution_count": 47, "metadata": { "ExecuteTime": { "end_time": "2018-10-24T05:19:45.441765Z", "start_time": "2018-10-24T05:19:45.431566Z" }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "(0.4940490278129096, 0.5559509721870904)" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "normal_two_sided_bounds(0.95, mu, sigma)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "\n", "if you were to repeat the experiment many times, 95% of the time the “true” parameter (which is the same every time) would lie within the observed confidence interval (which might be different every time).\n", "\n", "we do not conclude that the coin is unfair, since 0.5 falls within our confidence interval." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "If instead we’d seen 540 heads, then we’d have:" ] }, { "cell_type": "code", "execution_count": 70, "metadata": { "ExecuteTime": { "end_time": "2018-11-06T13:33:23.900320Z", "start_time": "2018-11-06T13:33:23.894828Z" }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.54 0.015760710643876435 (0.5091095927295919, 0.5708904072704082)\n" ] } ], "source": [ "p_hat = 540 / 1000\n", "mu = p_hat\n", "sigma = math.sqrt(p_hat * (1 - p_hat) / 1000) \n", "print(mu, sigma, normal_two_sided_bounds(0.95, mu, sigma) )" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Here, “fair coin” doesn’t lie in the confidence interval. \n", "- The “fair coin” hypothesis doesn’t pass a test that you’d expect it to pass 95% of the time if it were true." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## P-hacking P值操纵\n" ] }, { "cell_type": "markdown", "metadata": { "ExecuteTime": { "end_time": "2018-10-24T05:42:33.915702Z", "start_time": "2018-10-24T05:42:33.910679Z" }, "slideshow": { "slide_type": "subslide" } }, "source": [ "- Test enough hypotheses against your data set, and one of them will almost certainly appear significant. \n", "- Remove the right outliers, and you can probably get your p value below 0.05.\n", "\n", "\n", "P-hacking是科研人员不断的尝试统计计算直到p<.05，当然有时这可能是无意识的。\n", "- 最早应该是美国宾夕法尼亚大学的Simmons等人提出来的。\n", "- P-hacking 按照字面的意思来看是「P值黑客]，但是实际上的意思是「P值篡改」或者「P值操纵」。\n", "- P-hacking危害很明显，那就是很容易引起假阳性，导致实验的不可重复性\n", "\n", "> Simmons JP, Nelson LD, Simonshohn U. False-positive psychology: undisclosed flexibility in data collection and analysis allows presenting anything as significant. Psychol Sci. 2011 Nov;22(11):1359-66." ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "ExecuteTime": { "end_time": "2018-10-24T05:41:20.421246Z", "start_time": "2018-10-24T05:41:20.408579Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "def run_experiment():\n", " \"\"\"flip a fair coin 1000 times, True = heads, False = tails\"\"\"\n", " return [random.random() < 0.5 for _ in range(1000)]\n", "\n", "def reject_fairness(experiment):\n", " \"\"\"using the 5% significance levels\"\"\"\n", " num_heads = len([flip for flip in experiment if flip])\n", " return num_heads < 469 or num_heads > 531" ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "ExecuteTime": { "end_time": "2018-10-24T05:41:55.809801Z", "start_time": "2018-10-24T05:41:55.446143Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "P-hacking\n", "46\n" ] } ], "source": [ "print(\"P-hacking\")\n", "\n", "random.seed(0)\n", "experiments = [run_experiment() for _ in range(1000)]\n", "num_rejections = len([experiment\n", " for experiment in experiments\n", " if reject_fairness(experiment)])\n", "print(num_rejections) # 46" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "\n", "- you should determine your hypotheses before looking at the data, \n", "- you should clean your data without the hypotheses in mind, and \n", "- you should keep in mind that p-values are not substitutes for common sense. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "## Example: Running an A/B Test\n", "\n", "To help choosing between advertisement A (“tastes great!”) and advertisement B (“less bias!”)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "If 990 out of 1,000 A-viewers click their ad while only 10 out of 1,000 B-viewers click their ad, you can be pretty confident that A is the better ad. \n", "\n", "But what if the differences are not so stark? \n", "- Here’s where you’d use statistical inference." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Let’s say that $N_A$ people see ad A, and that $n_A$ of them click it.\n", "- We can think of each ad view as a Bernoulli trial \n", " - where $p_A$ is the probability that someone clicks ad A. \n", " \n", "$\\frac{n_A}{N_A}$ is approximately a normal random variable with mean $p_A$ and standard deviation $\\sigma_A$\n", "\n", "$$\\sigma_A = \\sqrt \\frac{P_A(1-P_A)}{N_A}$$ " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Similarly,\n", "\n", "$$\\sigma_B = \\sqrt \\frac{P_B(1-P_B)}{N_B}$$ " ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "ExecuteTime": { "end_time": "2018-10-27T01:14:13.336748Z", "start_time": "2018-10-27T01:14:13.332300Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "# running an A/B test\n", "def estimated_parameters(N, n):\n", " p = n / N\n", " sigma = math.sqrt(p * (1 - p) / N)\n", " return p, sigma" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "If we assume those two normals are independent (which seems reasonable, since the individual Bernoulli trials ought to be), then their difference should also be normal with \n", "- mean $p_B − p_A$ \n", "- standard deviation $\\sqrt {\\sigma_A^2 + \\sigma_B^2}$.\n", "\n", "Null Hypothesis: $p_A$ and $p_B$ are the same, that is $p_A - p_B = 0$" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "ExecuteTime": { "end_time": "2018-10-27T01:14:17.226285Z", "start_time": "2018-10-27T01:14:17.220053Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "def a_b_test_statistic(N_A, n_A, N_B, n_B):\n", " p_A, sigma_A = estimated_parameters(N_A, n_A)\n", " p_B, sigma_B = estimated_parameters(N_B, n_B)\n", " return (p_B - p_A) / math.sqrt(sigma_A ** 2 + sigma_B ** 2)\n", "\n", "# it should approximately be a standard normal." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "For example, if “tastes great” gets 200 clicks out of 1,000 views and “less bias” gets 180 clicks out of 1,000 views, the statistic equals:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "ExecuteTime": { "end_time": "2018-10-27T01:17:48.099418Z", "start_time": "2018-10-27T01:17:48.090969Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1000 rejections out of 1000\n", "A/B testing\n", "a_b_test_statistic(1000, 200, 1000, 180) -1.1403464899034472\n", "p-value 0.254141976542236\n", "which is large enough that you can’t conclude there’s much of a difference. \n" ] } ], "source": [ "num_rejections = 1000\n", "print(num_rejections, \"rejections out of 1000\")\n", "print(\"A/B testing\")\n", "\n", "z = a_b_test_statistic(1000, 200, 1000, 180)\n", "print(\"a_b_test_statistic(1000, 200, 1000, 180)\", z)\n", "print(\"p-value\", two_sided_p_value(z))\n", "print('which is large enough that you can’t conclude there’s much of a difference. ')" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "On the other hand, if “less bias” only got 150 clicks, we’d have:\n" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "ExecuteTime": { "end_time": "2018-10-27T01:18:49.148965Z", "start_time": "2018-10-27T01:18:49.143574Z" }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "a_b_test_statistic(1000, 200, 1000, 150) -2.948839123097944\n", "p-value 0.003189699706216853\n" ] } ], "source": [ "z = a_b_test_statistic(1000, 200, 1000, 150)\n", "print(\"a_b_test_statistic(1000, 200, 1000, 150)\", z)\n", "print(\"p-value\", two_sided_p_value(z))\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "- it means there’s only a 0.003 probability you’d see such a large difference if the ads were equally effective.\n", "- there’s only a .3% chance you’d observe such an extreme statistic if our null hypothesis were true." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "## Bayesian Inference" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "An alternative approach to inference involves treating the unknown parameters themselves as random variables. \n", "\n", "- To start with a `prior distribution` for the parameters \n", "- and then uses the observed data and Bayes’s Theorem to get an updated `posterior distribution` for the parameters. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "For example, when the unknown parameter is a probability (as in our coin-flipping example), we often use a `prior` from the `Beta distribution`, which puts all its probability between 0 and 1:\n" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "ExecuteTime": { "end_time": "2018-10-27T01:19:25.637110Z", "start_time": "2018-10-27T01:19:25.629514Z" }, "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "def B(alpha, beta):\n", " \"\"\"a normalizing constant so that the total probability is 1\"\"\"\n", " return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta)\n", "\n", "def beta_pdf(x, alpha, beta):\n", " if x < 0 or x > 1: # no weight outside of [0, 1]\n", " return 0\n", " return x ** (alpha - 1) * (1 - x) ** (beta - 1) / B(alpha, beta)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "The larger alpha and beta are, the “tighter” the distribution is." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "- if alpha and beta are both 1, it’s just the uniform distribution (centered at 0.5, very dispersed). \n", "- If alpha is much larger than beta, most of the weight is near 1.\n", "- if alpha is much smaller than beta, most of the weight is near zero. " ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "ExecuteTime": { "end_time": "2018-10-27T04:25:59.497947Z", "start_time": "2018-10-27T04:25:59.490577Z" }, "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "x = [i*10/1000 for i in range(101)]\n", "y1 = [beta_pdf(i, 1, 1) for i in x]\n", "y2 = [beta_pdf(i, 10, 10) for i in x]\n", "y3 = [beta_pdf(i, 4, 16) for i in x]\n", "y4 = [beta_pdf(i, 16, 4) for i in x]" ] }, { "cell_type": "code", "execution_count": 56, "metadata": { "ExecuteTime": { "end_time": "2018-10-27T04:28:30.049822Z", "start_time": "2018-10-27T04:28:29.872764Z" }, "slideshow": { "slide_type": "skip" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsnXd8lFX2h593Jr2HMJESCTUhNIEEDEhZCB0EQcEKurq6uhI7Vlx113V15aeuILi6FlZARRQUlYBKExCBhCYkJPROJj0hbTJzf3+8TKSFTJKZeafcx08+Mpl37vslCSd3zj3nexQhBBKJRCJxH3RaC5BIJBJJw5CBWyKRSNwMGbglEonEzZCBWyKRSNwMGbglEonEzZCBWyKRSNwMmwK3oigRiqIsURQlS1GUTEVR+jlamEQikUguj4+N1/0bSBNC3KQoih8Q5EBNEolEIrkCSn0NOIqihAM7gPZCdutIJBKJ5tiy424HGIGPFEW5BkgHHhZCnD3/IkVR7gPuAwgODk7s3LmzvbVKJBKJx5Kenp4nhDDYcq0tO+4kYDNwnRDiV0VR/g2UCCGer+s1SUlJYtu2bQ3RLJFIJF6NoijpQogkW6615XDyOHBcCPHrucdLgN6NFSeRSCSSplFv4BZCnAaOKYoSf+5TKcBeh6qSSCQSSZ3YWlWSCiw8V1FyEPij4yRJJBKJ5ErYFLiFEDsAm3IvEs/FZDJx/PhxKisrtZbitQQEBBATE4Ovr6/WUiQaYuuOWyLh+PHjhIaG0rZtWxRF0VqO1yGEID8/n+PHj9OuXTut5Ug0RLa8S2ymsrKSqKgoGbQ1QlEUoqKi5DseiQzckoYhg7a2yK+/BGTglkgkErdDBm6JRCJxM2TglrgVer2enj17cs0119C7d282bdp0xeuLioqYO3euTWtXVFQwePBgzGYzAKNGjSIiIoJx48bZ9Pr169fTu3dvfHx8WLJkSe3njUYjo0aNsmkNicQWZOCWuBWBgYHs2LGDnTt38s9//pNnnnnmitc3JHB/+OGHTJo0Cb1eD8CMGTP45JNPbNbWpk0bPv74Y2677bYLPm8wGGjZsiUbN260eS2J5ErIckBJo3hp+R72niyx65pdWoXxwvVdbb6+pKSEyMjI2sevv/46ixcvpqqqiokTJ/LSSy/x9NNPc+DAAXr27Mnw4cN54YUXmDBhAoWFhZhMJl5++WUmTJgAwMKFC1m0aFHteikpKaxdu9ZmPW3btgVAp7t0P3TDDTewcOFCrrvuOpvXk0jqQgZuiVtRUVFBz549qays5NSpU6xevRqAVatWkZOTw5YtWxBCMH78eNavX8+rr77Kb7/9xo4dOwCoqalh6dKlhIWFkZeXR3JyMuPHj8dkMnHw4MHa4GtvkpKSmDlzpkPWlngfMnBLGkVDdsb2xJoqAfjll1+YNm0av/32G6tWrWLVqlX06tULgLKyMnJycmjTps0FrxdC8Oyzz7J+/Xp0Oh0nTpzgzJkzWCwWIiIiHKY7OjqakydPOmx9iXchA7fEbenXrx95eXkYjUaEEDzzzDP8+c9/vuCaw4cPX/B44cKFGI1G0tPT8fX1pW3btlRWVhIeHu7QxpbKykoCAwMdtr7Eu5CHkxK3JSsrC7PZTFRUFCNHjuTDDz+krKwMgBMnTpCbm0toaCilpaW1rykuLiY6OhpfX1/WrFnDkSNHAIiMjMRsNtsUvJ955hmWLl3aIK3Z2dl069atQa+RSOpC7rglboU1xw1q2mP+/Pno9XpGjBhBZmYm/fqpc6xDQkJYsGABHTp04LrrrqNbt26MHj2ap556iuuvv57u3buTlJTE+ZOaRowYwYYNGxg2bBgAAwcOJCsri7KyMmJiYvjggw8YOXIku3fvZvz48Zdo27p1KxMnTqSwsJDly5fzwgsvsGfPHgDWrFnD2LFjHf3lkXgJ9U7AaQxyAo5nkpmZSUJCgtYyHEZGRgZvvvlmvSWAI0eOZOXKlQ1ae9CgQXz99dcXVME0Fk//Pngr9p6AI5F4Bb1792bIkCG1DTh10dCgbTQaeeyxx+wStCUSkKkSieQC7r77bruvaTAYuOGGG+y+rsR7kTtuiUQicTNk4JZIJBI3QwZuiUQicTNk4JZIJBI3QwZuiVvhCrauhw4d4tprr6Vjx47cfPPNVFdXX3Hd/Px8hgwZQkhICNOnT7/gufT0dLp3707Hjh156KGHsJbnPvHEE7U+LBLJxcjALXErXMHW9amnnuLRRx9l//79REZG8sEHH1xx3YCAAP7+978za9asS5574IEHeP/998nJySEnJ4e0tDQAUlNTefXVV23SLfE+ZOCWNJ6Pxl76seV99bnq8ss/v32h+vzZ/EufayCXs3Xt06cPPXr04IUXXgC4wNZ1xowZlJWVkZKSQu/evenevTtff/117esXLlxYa/EKqq1raGjoBfcUQrB69WpuuukmAO68806WLVt2RZ3BwcEMGDCAgICACz5/6tQpSkpKSE5ORlEUpk2bVrtWbGws+fn5nD59usFfF4nnI+u4JW6F1rau+fn5RERE4OOj/tOJiYnhxIkTjfq7nDhxgpiYmNrHF6/Vu3dvNm7cyI033tio9SWeiwzcksbzx+/qfs4v6MrPB0dd+fk6cFdb18YgrWAldSEDt8Rt0cLWNSoqiqKiImpqavDx8eH48eO0bt26Ufpbt27N8ePHax9fvJa0gpXUhcxxS9wWLWxdFUVhyJAhtcOA58+fX5sXX7p0ab2HpefTsmVLwsLC2Lx5M0II/ve//12QY5dWsJK6sGnHrSjKYaAUMAM1tjpYSST2xhVsXV977TVuueUWZs6cSa9evbjnnnsAOHDgAGFhYZfV3bZtW0pKSqiurmbZsmWsWrWKLl26MHfuXO666y4qKioYPXo0o0ePBsBkMrF//36SkuQ/Ncml2GTrei5wJwkh8mxZVNq6eiaebidqq61rXdxxxx28+eabGAyGJmtZunQpGRkZ/P3vf7/kOU//PngrDbF1lTluieQc59u6Wmu5G8KCBQvspqWmpobHH3/cbutJPAtbA7cAVimKIoD/CCHec6AmiUQzHGHr2hgmT56stQSJC2Nr4B4ghDihKEo08IOiKFlCiPXnX6Aoyn3AfcAlJVgSiUQisR82VZUIIU6c+38usBToe5lr3hNCJAkhkuyR45NIJBLJ5ak3cCuKEqwoSqj1z8AI4DdHC5NIJBLJ5bElVXIVsFRRFOv1i4QQaQ5VJZFIJJI6qXfHLYQ4KIS45txHVyHEP5whTOIZGMuN3JV2F3kVNlWS1oszbV1BNbKKiYm5xI71cnzxxRd07doVnU7HxeWwu3btol+/fnTt2pXu3bvXNvoMGzaMwsJCm/RJJFZk56TEoby7610yzmQwb+c8u6znTFtXgOeff55BgwbZ9Ppu3brx1VdfXXJ9TU0Nd9xxB++++y579uxh7dq1+Pr6AjB16lSb9UkkVmQdt6RRvLblNbIKsup8Pv1MOoLfm7sW71vM4n2LUVBIvCrxsq/p3KwzT/V9ymYNl7N1Xbx4MVVVVUycOJGXXnrpAlvX4cOH88ILLzBhwgQKCwsxmUy8/PLLtW3mCxcuZNGiRb//HdLTOXPmDKNGjbpkB3056mqKWbVqFT169OCaa64BVL8TK+PHj2fgwIE899xzNv+9JRIZuCUOoXvz7hwvPU5hVSECgYJCZEAkV4dc3aR1nWXrarFYePzxx1mwYAE//vhjkzRnZ2ejKAojR47EaDRyyy238OSTTwKqR0pVVRX5+fkXBHSJ5ErIwC1pFLbsjP/2y99Ykr0EP70fJrOJYbHDeD75+Sbd11m2rnPnzmXMmDEX+GU3lpqaGjZs2MDWrVsJCgoiJSWFxMREUlJSgN/tW2XgltiKDNx2xFhuZMb6GcwaPIvmgc0veextFFQWMCV+CpPjJvNF9hd2O6C04khb119++YWff/6ZuXPnUlZWRnV1NSEhIY0aJxYTE8OgQYNo3lz9GRgzZgwZGRm1gVvat0oaijyctCMXH8TZ+2DO3XhryFvMTJ5JfLN4ZibP5K0hb9l1fUfaui5cuJCjR49y+PBhZs2axbRp02qD9rRp09iyZYvNOkeOHMnu3bspLy+npqaGdevW0aVLF0B9B3D69Ol6J+9IJOcjd9x2IHFBItXm3yd9Ww/iLn6soLB6ymqv3H3bC2fautbFrl27aNWq1SWfX7p0KampqRiNRsaOHUvPnj1ZuXIlkZGRPPbYY/Tp0wdFURgzZgxjx6ozNtPT00lOTq4dheZNnN28meJvluPTLBKfVq3wb9+eoGuv5VzPiOQK2GTr2lC8ydbVWG7kkTWPEBUYxbrj67AIyyXX6BQdLYJacOrsKSbHT25ynlcrPN1O1BZb15KSEu655x6++OILu9zz4YcfZvz48bVpE1tw9+/D6b+/TOHChVe8pvmDD2JIrb923pNoiK2rTJU0kXd3vcvuvN1sOrkJi7CgV/QoKLQPb4+CunOwCAsnz55EIFi8bzHd53cnccHlS+Ik2nG+rWtdhIWF2S1og1r73ZCg7e6UrFhB4eLF6IKCMDz2GPE7d9B57x46rl9H8788AEDwgAE0++MfNVbq2njf+zM7cXF6pMpcBaijrabETWH10dVMiZ/CsDbDeHHTi5w4q07v9tP7MbDVQIyVRvIq8mTaxMVwtq3rvffe69T7aYVx9hzy3nmn9rHFZML4xhuIqioMqdPxjY7G8NBD+LZqxakXXuTItKnEfvQR+vBwDVW7LjJwN5K0SWnM2jaLlYdXYhZm/HR+DI8dzhN9nqB5YHNmJs+svbZ/6/4syV6CQFBtruZA8QGOlBxh3s55bps2kUgaQsSkiRR+9hm6oCDaLv4cn/Map87HdOo0mM1U7c0k+9rk2s97Y+rkSshUSSMxBBmoqKnALMzoFT0mi4lgv+DL7qCtZXE+OvX35OGSwzJtIvEazGVnOXb/A4jqaq5+d16dQRvAkDqdhKxMIu+4A3Q62i39ioSsTBm0L0IG7kZgLDdy54o72XJqC4E+gXw86mOmxE8hvyL/stdby+JW3biK3tG9az8foA9gbLuxrLxxpbOkSyROwzh7DpmdE8hOSqIqJwdLaSkHx47DOHtOva81pE5HHx7O6Zf/gSMKKNwdmSppBO/uepftudsRCF4f/Do9o3vSM7pnva8zBBnoENGBjNwMQM2L17VLl0jcHUPqdEKHD+PQpBuJvO02Wsy03Y9FHx6O4bFHOf38Xyn59lvCr7/egUrdD7njbgCJCxLpPr87i/ctrjVQmrFuRoNSHQWVBUzoMIEAfQAtg1vWuUuXXB5n2rqOGjWKiIgIxo0bd8F1Qgiee+454uLiSEhI4O2337Zp/ctZxHqyrasQgjP/eAV9WFijUh0RN95IQPfu5P7rdcxlZx2g0H2RgbsBpE1KY0y7Mfgo6hsVP71fg1Mdbw15i5cHvMwDPR/g5NmTDI8dble/ak/HmbauM2bMuGxN98cff8yxY8fIysoiMzOTW265xab1L2cR68m2rqUrV1K+dSuGRx5uVHVI3jtzqdy9mxqjkeykJDI7J5DZOcGmVIunI1MlDcAQZCDIJ4gaUYOCgslc94FkfdyRcAdf5XzFy5tf5qzprNtVmJx+5RWqMuu2dW0M/gmdafHsszZf72hb15SUFNauXXvJfefNm8eiRYvQ6dR9T3R0dL1a67KI9VRbV0tFBWde+xf+nTsT0ciJ9YbU6RhSp3Noys1Yys/Sfvly2VV5Dhm4G8j+ov0APJr4KCfKTjR6p9zv036XbZP30/uRfke6XbRqSfWJE9ScPHnJ531atcKvdetGr+ssW9crceDAAT7//HOWLl2KwWDg7bffplOnTnVefyWLWE+1dS34ZAE1p07R+l+voZw3mKIxRE6ZzKmZz1OxfQdBvXvZSaF7IwN3A4kKjKJZQDPuSLgDX71vo9ex1oGvOLQCgSBAH0BKmxSe6POEHdU6jobsjO2Js2xdr0RVVRUBAQFs27aNr776irvvvpuff/65zuvrs4j1JFvXixttjkydBjStDjts9GjOvPJPihYvloH7HDJwNwBjuZG1x9Yyrcu0JgVtUNMuwb7BtY9lhUnDcaSt65WIiYlh0qRJAEycOJE/1tOeXZ9FrCfZuhpSp+NjaM7pF1+izfz5BF/bt8lr6oKDCbv+eoq//pqrnn0GfViYHZS6N/Jw0kaM5UamrpiKWZiZ1GmSXdYsqCxgcvxkYkJiCPMLkxUmDcSRtq5X4oYbbmDNmjUArFu3jri4OAC2bNnCtGnTLrn+ShaxnmbrKsxm8j/8iIAePQjq28du60ZMmYyorKR4+XK7renOyB23jczbOY8TZScwBBpoG97WLmta/amX7V/G8xufr60w8dbBC7bgTFvXgQMHkpWVRVlZGTExMXzwwQeMHDmSp59+mttvv50333yTkJAQ/vvf/wJw9OjRBu+cPc3WtXTVKkxHjxL9xON2PUgM7NqVgC5dKPp8MZG33eb1h5TS1rUeLjaTsmLPQ0ST2cSor0YhhCCvIs9lrV/d3U60Pmyxdb0SM2bMYOrUqfTo0cPm13iSrasQgsM33oTl7Fnaf/9dkw8lL6bws885/eKLtF38OYEN+Bq7C9LW1Y5Ya7f1ivpD6K/3t3ubevKnyeSW52KsMEoPEw2xxdb1Srz++usNCtrgWbau5b/8QuXevTS75267B23j7DmcfvFFAA5Pudnra7pl4K4H6yGiWZjRoaPaXG33Q8S0SWmMiB1R+9iVPUw83Tfi7rvvrm3AcQYNtXV1xa+/1ZPk6N33AHD6ry/YPahazaeC+iXj16EDCVmZXm0+JQO3DRwpUQ+w7r/m/iuaSTUWQ5CBcP/fO8tctcIkICCA/Px8lwwe3oAQgvz8fAICArSWcgGG1Ol0XLcWdDqi7r3XoUE1dMgQqg8coProUbuv7U54xomIg+nTog9bT29lcvxkhwXTgsoCRsaOZOWRlfQ09HTJCpOYmBiOHz+O0WjUWorXEhAQUGc9uJYUffklWCxETL7JofcJ+cMfOPPKPylbu5Zml6ng8RZk4LaBdcfW0b15d4fugK0VJie/O0lxdTHzR8932L0ai6+vL+3atdNahsTFEGYzRV8sIbh/P/wuaniyN35t2uDXoQOla9Z4deC2OVWiKIpeUZTtiqJ860hBroax3Mhv+b8x+OrBTrnf5LjJHCw+yPbc7U65n0TSVM5u2EDNqVNETLnZKfcLHfIHyrduw3xefb630ZAc98NApqOEuCo/n1BbmQfHOCdwj2w7kmDfYBbsXSBdAyVuQeHiL9BHRRE6dIhT7hcyZAjU1HB240an3M8VsSlwK4oSA4wF/utYOa6FsdzIG+lvEB0YTVxknFPuGeQbxLj24/jp6E9knMlg3s55TrmvRNIYTGfOULZ2LRGTJqL4+TnlnoE9e6KPiKDsXPeqN2Jrjvst4EkgtK4LFEW5D7gPuMTYx12Zu2MuxVXFdAzv6LROrYsbfjzNNVDiWRR/9RWYzY22bm0Mil5PyOBBlK1bjzCb7V4z7g7Uu+NWFGUckCuEuGLUEEK8J4RIEkIkGQwGuwnUAuukmyU5SwDYX7zfaQ0x1oYfBfUXhSvXdEu8F2vttvHf6vSfAyNGOrUhJmTIEMxFRVScc4r0NmxJlVwHjFcU5TDwGTBUUZQFDlWlMc7olqwLa8OPdTSaq9Z0S7wbQ+p02n7+GQAt//GyUxtijLPncOKRRwE4cvsdXtlFWW/gFkI8I4SIEUK0BW4BVgsh7nC4Mg05v1tSQXFIt+SVsM6lVFDo3KyzS9Z0SyTF3yxH8fcndMSI+i+2I9YuyoAuXQi69lqv7KKUddx1kFueC8BNcTehU3ROre6w1nTnVeZxsOggn437zGn3lkhsQZhMlHz/PSFDh6APrfPoy6EE9elD4WefYamuRuekg1FXoUEt70KItUKIcfVf6f7cHK/WpI5sO5KZyTNrg6kzub799Zw6e4r0M/JQUuJalP28AXNhIeHjx2umIahPEqKqisrfftNMg1ZIr5I6SD+Tjo/Ohx4G7ewjh7YZSpBPEEuyl8iabolLUbz8G/SRkYQMGKCZhsBEtVigfMtWzTRohQzcdbDtzDa6RXUj0Ee7kVKBPoEMix3GD4d/kDXdEpfBXFpK2U+rCRszBsW3aSP8moJPZCT+nTpR7iHe/w1BBu7LUFFTwZ68PSRepa0fduKCRL458A0mYZI+3RKXoXTlSkR1NeETtEuTWAnqk0RFRgaipkZrKU5FBu7LsNO4kxpRQ1ILm4ZROIy0SWmMbju69rGs6ZZoibV2+9RMdTqTdaCBlmV4QX36YCkvpzLTu9w4ZOC+DOln0tEpOnoaemqqwxBkIMQvpPaxrOmWaIkhdTqdNvys+m4/cL9LlOF5a55bBu7LsO30Njo363xB0NSKgsoCUtqoo636tugra7olmlKyahVYLISNHl3/xU7ANzoav9hYr8tzy8B9ESdKT5B+Jp2uUV21lgKoNd1v/OENWgS3IMAnQJOyRInESun3K/Dr2IGAOOeYrtlCUN8+lKenIxo5K9QdkYH7Iv619V8IBCfLTmotpRadomNk7Eg2ntxIcVWx1nIkXorpTC7l6ekus9u2EpSUhKWkhKqcHK2lOA0ZuM9hNZZafWw1ABtPbnSpCo5R7UZRY6lh9dHVWkuReCmlK9NACMJGj9FaygUE9ekDQPlW70mXyMB9DquxlE5RvySuVsHRNaorMSExLD+wXDbjSDSh5PsV+HfujH971xlfZ5w9h/1D1TOgM//4h9cYTsnAfQ6rsZRFWNApOper4FAUhZFtR7L1zFbZjNMQynJhXxoc3vD757YvgOxVcFYe9NqK6eRJKnbscLk0idVwKmToUPzatXOJShdnIE2mzuNU2SkA7ux6J+Wmcpfa1Z4/YMHajCMHLFyB07th0xz47UuwmCBuFLQdAGYTLH8YLDWg94ceU6DfdIjurLVil6ZkRRoAYWNcK3BbCezRnbLVqzGXlKAPC9NajsORgfs8boq7iQ0nN5DSJoVrDNdoLecC0ialMWvbLFYcWoFAEKAPIKVNCk/0eUJraa7HNw9BxnzwDYY+90DXidD8XBWEzgce3wd52bD7C9ixCLZ/AmP/D/r8SVvdLohx9hzy3nmn9vGB4aqFa/MHH3SpXW1At+4AVP72G8H9+2usxvHIwH0ev+X/ho/iQ+dmrrf7kgMWGkC3GyGijRq0AyMvfE5RILi5+hHbH4bMhB0LoMe5CeVmE+i1899wNQyp04m4cRL7h6ZgePwxmt97r9aSLktgN7V8t2K3dwRumeM+j915u4lrFoe/3l9rKZeloLKAEbHqjifpqiTZjGNFCFj7qvoB0H4wDHri0qB9OYKj4LqHwT8UTBXw/hDYNNuxet2M0h9+ACDMyQMTGoI+IgK/2Fgqdu/SWopTkIH7HBZhYU/eHrpFddNaSp28NeQtZg2eReuQ1vj7+MtmHFCD9g9/hbX/hKKj6uPGYqmBZu1h1UxY97r9NLo5JStX4d+5M36xsVpLuSIB3btTuds7vLll4D7HkZIjlJnK6NbcdQM3qNUlI2JHsPnUZtmMIwT89DfY9DYk3QMT3lFTIY3FPxRu+gh63AJrXoaf37CfVjfFdOYMFRkZhI103d22lcAe3ak5cwbTmTNaS3E4MnCf47c89Te1qwdugOGxw6mx1LDu+DqtpWjL2ldhwxuQeBeMmdW0oG1Fp4cb5kL3yfDTS7Dl/aav6caU/vAjAKEjR2qspH5qDyh379ZYieORgfscu/N2E+gTSPvw9lpLqZduzbvRIrgF3x741rubccJj1J322DdBZ8cfZZ0ebngXBj4Bca4fsBxJ6cqV+HfqiH971/93EdAlAfR6KrwgXSID9zn25O2ha1RX9Dq91lLqRVEUhscOZ/Opzd7ZjGM5ZybUeyqMe8O+QduK3gdSnlerUywWr2zWqcnLo3zbNkJHuMcvL11AAP7xcVR6wQGlDNzAydKT7M7b7Ra7bVCbcT7Z+wni3H9eNRmnogj+Mwj2LHPePb+ZDvOvh+qzzrunC1D6448gBKFukN+2EtitOxW/7UFYLFpLcSgycAOzts1CIDhWekxrKTaRNimN0e28cDKOELD0fjBmQWgL5923242QuxeWP+K8e2qIddLN6RdfAuDQ+Alu4/8R2KM7lpISqo8c0VqKQ/HqBpzz28gBfjn1C93nd3f5NnJDkIEQXy+cjLP9E8heASP/CW2SnXffjikw5FlY8w/oPEbtxPRgDKnTibzjdnIGDCTqnnuIfuxRrSXZTED33w8o/du5jhmWvfHqHbfVEVCvqHltd9q5FlQWMCRmCAD9W/X3/Gac4uOw8jloOxCuvd/59x/wGLTqBd89Dmc9/zC4bPUaMJsJdeGmm4sxzp7DoQk3AHDyyac82inQq3fc1jZyszCjoLjVzvWtIW9RY6lh6OKhhPmH8a9B/9JakmPZt0I9lBw/2zGHkfWh94EJc+GzW9VGn2DX/xlpCqWrVuHbujUBXbtoLcVmDKnTMaRO59CUm9EFBhI7/2OtJTkMrw7cAHkVeSgoXN/hegJ9At2qtM5H58PQNkNJO5xGtbkaP72f1pIcR997ofM4CGupnYarusD0dDWIezDmsjLObtpE5O23o9ijNt7JBHTuTElaGkIIt9RvC16dKgGY3ms6AkH/Vv2ZmTzT7drIU9qkcNZ0ls2nNmstxTGUGVWLVtA2aFvR+0BNNayfpVa4eCBla9chTCZCRwzXWkqjCOiSgKWkhJqTrjN+0N7UG7gVRQlQFGWLoig7FUXZoyjKS84Q5iz2FewDcElHQFtIbplMqG8oPxz5QWspjuGnF+G/w1yrjjovG1a/DOte01qJQyhdtQofg4HAnj21ltIoAjqr/5Yrs7I0VuI4bNlxVwFDhRDXAD2BUYqiOPFI37FkFWThr/cnNsy1DXTqwlfvy+CrB/PT0Z+4c8WdbpXqqZcT6eq0mr73qi5+rkKLbmqb/a//gVzPCg6WigrKfv6Z0OHDULQ4S7AD/nFxoChUZnrW9+Z86v3OCJWycw99z300wYLNtdhXsI+4yDh8dO6btxzWZhil1aVsz93uOV2UFguseAo8+UM6AAAgAElEQVSCo2HQk1qruZShz4N/CKQ91TRHQhejbMMGREUFocPdM00CoAsKwq9dOyqzMrWW4jBs+pWqKIpeUZQdQC7wgxDi18tcc5+iKNsURdlmNBrtrdMhCCHILMgkvlm81lIaTeKCRB5ZqzaGeFQX5a7P4fhWGPYiBLjgKKrgKPjDs3BwLWR9p7Uau1G66gf0ERG1k9PdlYDOnana6+WBWwhhFkL0BGKAvoqiXGKhJ4R4TwiRJIRIMhgM9tbpEE6fPU1JdQmdI90zvw2uP52+0VQUQuwAuOZWrZXUTZ97VH3hMVoraTLWbsmS5csxFxWR1a27W9dA+yd0xnTyJOZiz7Q+blB+QAhRpCjKGmAU4PYWXFkFag7MnXfctSPNhIeNNOv3F7XRxpXzrHpfmPiu1irsgiF1OoE9unPsz/dz9X/eJWTwYK0lNYmAzgkAVGbtI/javhqrsT+2VJUYFEWJOPfnQGA44BFZ/6zCLBQU4iLjtJbSJAoqC5jUaRI+ig+dIjq5dxdleQFkr1Tzxq4ctM+n9DSkPQtVpVoraRIlq1ahCwkhqF8/raU0mYAE9V10lYfmuW3ZcbcE5iuKokcN9IuFEN86VpZz2Fewj9iwWIJ8g7SW0iSstef5FflkFWbx5h/e1FhRE9jwBmyaA6npENVBazW2UXwcNr8DgREw2AUPUm1AmEyU/fgTIUOGoPNz/0Yun+bN8TEYPLayxJaqkl1CiF5CiB5CiG5CiL85Q5gzyCrIctv67csxLHYYp8+erp3m43YUn4Bf31Pzxu4StAFiktSuzo1vu1a9eQMo37oVc3Gx2zbdXA7/hM5UZnrmjttN3ovan4PFBzlRdoKYUPc/WLLyh6v/gI/iww9H3bQZZ91rICzwh6e1VtJwhj4PprPqOwY3pGTVKpTAQEIGDNBait0I6JxA1YEDWKqr67/YzfDawP1WuppeyC7I1liJ/Qj3D6dvy76sPLSSu1a42Uiz/ANqs02feyDSDZuhojur7xS2vA8l7tVqLcxmSn/8iZBBg9AFBmotx24EdEmAmhqq9+/XWord8brAnbggke7zu7Pm2BoA1p9Y7xl1z+cYHjuck2dPkpHrZiPNio9BZFvVPtVdGfwkdL9JaxUNpmL7dsx5eW4xyb0h1La+e2CeWxEO6PpKSkoS27Zts/u69sBYbmTWtlmsPLwSszAToA8gpU0KT/R5wu1L6C4eDGHF1QdD1CKEfSa1SxrE6Vdeoeizz+m0aRP6kGCt5dgF4+w55L3zziWfb/7ggxhSp2ugqH4URUkXQiTZcq3X7bjd2YO7PmqbcXCzZpzDG1XHPU8J2ie3w/aFWquoF2vTTeH/PkFUV5OdlOTWTTfnY0idTkJWJgHX9CDo2mtJyMokISvTZYN2Q3Ffg44mkF+Rjw4dY9qPIdg32L1ywVegthkHN2rGKToK/xsP/R+CYS9orcY+bH4X9n4NnUZAiOt2ERtSpxM84DqO3Hobrf71GuHjx2stye4ExMVT+sMPHufN7XU7boBHEh/BgoXklslu6cF9JQoqCxjfQf0H2CWqi+s342x4ExQd9PmT1krsx6AnwFwFv8zWWkm9lKatRPH1JWTIEK2lOAT/uDjMRUXUuIl/kq145Y47pzAHgE6RnTRWYn+sv4SOlx2nqLLItX8pFZ+AjE+g91QIb621GvvRvJM6GX7Lf+G6RyComdaKLouwWChZtYrg665DHxqqtRyH4B+ndkVXZefgGx2tsRr74ZU77uzCbHSKjvbh7bWW4jBGxI7gQPEBDhQd0FpK3Wyeq9ZtX/eI1krsz8DH1bruLe9praROKnftoubUKUJHjdRaisPwj1M3Z1X79mmsxL54ZeDOKcwhNiyWAJ8AraU4jOGxw1FQWHV4ldZSLo8QcCJDLZ9zx7rt+ohOUOu6XfhnrCRtJfj6Ejp0qNZSHIZPZCQ+0dFUZXtOvwZ4a6qkKIeEZglay3AohiADvaJ78f2h7/n19K/MGjzLtQ4pFQX++D1Un9VaieNwYedAIQQlq1YS0r8/+jAX9Du3I/5xcVTmeFbg9rodd7mpnGOlxzwyv30xI9qO4HDJYTLOuFgzjqlC9dtWFHWKjCdjsUD2KjCbtFZyAZW7dlFz8hSho0ZpLcXh+MfHUb3/AKKmRmspdsPrAvf+IrX91d2tXOsjcUEir255FXDByTgZn8Cb3aDomNZKHM/h9bBoMuxeorWSC6hNk6R4bprESkBcHKK6muojR7SWYje8LnBnF6pvmTx9x21txlFQa1ddphnHbIJNs+GqbhBxtbZanEG7wRDdFTa+pe6+NcbadFPw0UdgMpHd91qPabqpi9rKEg86oPS6wJ1TmEOgTyCtQzyo/OwyWJtxrLhMM86epVB8FAZ4YCXJ5VAUGPAoGLMgR/sOVkPqdGI/XQRAq9de9biOwsvh16ED6PVUetABpfcF7qIcOkV0qp3R6MlYm3EUFBKiErRvxhFC3W03j4NOnluCdgldb4Dwq9UBES5AyfcrUPz8CElJ0VqKU9D5+eHXri1V2TlaS7EbXlVVIoQguzCbYW2GaS3FKVibb06Xn+b02dPaT8Y5vUv9uP5t9xlLZg/0vpD8AGyep45m07AhR5jNlKalETJ4EPoQDz8YPo+AuHgqdu7UWobd8KJ/PerEm+KqYloGt9RailMZ3XY0R0qOkFmg8TSQltfA/Rugx83a6tCCPn+Ch3Zo3kVZvi2dGqORsDFjNNXhbPzj4jCdOIG5rExrKXbBqwL3nO3qW9U9+Xs0VuJchsUOw0fxIe1QmnYirPbBLbqDr+s2pTgMH3/Q+0BNFVQWayaj5PvvUYKC3H6Ke0M5v/XdE/CKwG0dnrD+xHoA1hxb4zqlcU4g3D+c/q37893B77SbjPPtI/DNQ86/rythqoTZibD2VU1uL0wmSleuJHTIEHRB7j0gu6EExFsDt2ccUHpF4LaWxukVPeBCpXFOZFTbUeRW5GozGaf0DOxYpOZ6vRnfAIjtD+nz1QYkJ3N282bMRUWEjfWuNIlx9hz2p6jnWqdffJHMzgluXwLpFYeTnjw8wRbOn4xjbcZZvG+x8ybjbH1frd9O/ovj7+Xq9JsOuz5Xg7eTSyJLvvseXWgowR40ENgWDKnTMaRO5/Ctt6Ho9cQu+ERrSU3GK3bc8PvwhHHtxzElfor2pXFOpHYyjqLBZJzqctj6AXQeC1EdHH8/V6dlD2g3CH79j9Pa4K1NN8XLlmEpLWVfj2vcfsfZGPw7daIqJwdHjGt0Nl6x4wZ4OPFhVh9bzbUtr2VCxwlay3EqtZNxhAaTcXZ9BhUF0O9Bx9/LXeg3HRZNgZxV6i80B2NInY5/p46ceORR2nz0IcH9+jn8nq6If1wcRYsXU5NrxPcq9/bm9podtycPT7CFgsoCboq7iWCfYNqEtnHeO45OI2DEP6CNdwaLy9JxOPxxBcQ7L9dc/M1yfKKjCerb12n3dDX8O53z5vaAA0qv2XHvL9rv8cMTrsT5k3C+Pfgti69f7Jwbh8dAf89tp24UOp16SOkkagoLKVu/nmZ3TkPR6512X1ejdqhCTg4hA907z+9VO+42oW08eniCLYxrP46Kmgp+OvqT42/240tweIPj7+OurH0Nvvqzw29TsmIF1NR45DDghuATGYmPweARO+56A7eiKFcrirJGUZS9iqLsURTlYWcIszf7i/Z7bZrkfHpG96R1SGu+O/idY290ejdseAOOb3PsfdwZc7VaYZLv2PFyJd8sxz8ujoD4eIfexx2wHlC6O7bsuGuAx4UQXYBk4EFFUbo4VpZ9qaip4GjJUTpGdNRaiuboFB1j2o1h08lN3P797Y5rxvnlHfANhsQ7HbO+J9D3XtD5wK+Om5RTffQoFTt2ED7+eofdw53wj4ujav9+hNmstZQmUW/gFkKcEkJknPtzKZAJuJUn6sGigwiE3HGfY1z7cQgEu4y7HNOMU3paHRzQ6w4IjLT/+p5CaAvoPhm2L3RYQ07x8uWgKISNdXz1ijvg36kToqqK6qNHtZbSJBqU41YUpS3QC/j1Ms/dpyjKNkVRthmNRvuosxM5RepbI7njVptxJnz9ezmkQybjbHkfLDWQfL/91vRU+v1FnQafPt/uSwshKPlmOUF9++Lb0ruM1eqi1rPEzdMlNgduRVFCgC+BR4QQJRc/L4R4TwiRJIRIMhgM9tTYZHIKc/DX+9MmtI3WUjTH2ozjo1MLivz0fvZvxglrCUl3QzPvrOBpEC26w7AXodNwuy5rnD2HrIQuVB85Qvmvv3pEm7c98O/YARTF7c2mbCoHVBTFFzVoLxRCfOVYSfZnf9F+2oe3R6/z3lIoK7Xt/xY1x1dtrrZ/M06fP9lvLW9gwKN2X9KQOh3T6VOUrkij08/r0QUH1/8iL0AXGIhvm6s9f8etKIoCfABkCiHecLwk+5NTmCPz2+dRUFnAlPgp9GvZDz+dH8ZyO6W2LBbY+43LTTR3C87sgdUv/25/20QsZ89SsiKN0DGjZdC+iIC4OLcvCbQlVXIdMBUYqijKjnMfbmMvVlRZhLHCSKcIGbitvDXkLWYmz2Rql6lUW6q5voOdKg72/wCLp0LmN/ZZz5s4uhnWvw5HNtlluZK0lYjyciImTbLLep6Ef6dOVB85gqWyUmspjcaWqpINQghFCNFDCNHz3Mf3zhBnD7ae2QpAdJB7exM4gv6t+hMdFM2y/cvss+Cm2RDWGhK8u9GjUVxzKwQ2U8so7UDR0q/wa9uWwF697LKeJ+EfFwcWC9UHD2otpdF4fOfkJ3tVC8cNJ2QH38XodXomdJjAz8d/5vbvmljTfWonHP4Zrv2z9N1uDH5B6tnAvu+b3JBTdegQFdvSCb9xEmqmU3I+Vs8Sd5767rGB2zr1ZnvudgCWH1zuVVNvbGVCxwlqTXdeE2u6N80BvxDoLRtuGk3fe0Hv1+Rdd/HSZaDTET7eu1wwbcUvNhbFz8+tK0s8NnBby94U1B2HN069qY/EBYmMWzqu9nGja7rNJjBmQu9pEBhhZ5VeREi0uutuwkBhUVND8bJlhAwc6PbWpY7AOHsOWd26I6qrKfjwQ7ctk/TYwF3rQY1Ap+i8buqNLVh/ufnq1NSGn66RNd16X/jzzzD0eQeo9DJGvQJDZzbqpdagVJObS9m6dW4blByJIXU6CVmZhN9wAz4GAwlZmSRkZWJIdS8HS4+2dT1RdgKA+7rfR2FVoTZDcl0Y6y+3GksNANWWRtR0V5eDMIN/qJqnlTQdIeDAaoi9Tp1TaSOG1OlU7NxJ1f79dPzxBxQfj/7n3ST84+MpXraMmoICfJo1/h2OVnjsjhtgSvwUAAbFDGJm8swLPKklKtaa7uvbX4+Cwsmykw1bYNuH8EZXKDnlGIHeyNHNsGCS6hzYAKqPHuXshg1ETL5JBu16cPep7x793c0uzEZBoUOEnHVYF9ZfZkdKjrD84HJ6RTegfMxsgs1z1bbtMOmFYTfaJEOLHmp5Za+p6uAFGyhavBj0eiJuusnBAt0f/3MWt1X79hGcnKyxmobj0TvunMIc2oS1IchXvoWvj9iwWJJbJvN51ufcueJO29JKv30JJSfgOre0aHddFEX9mubnQPYKm15iqa6m6MuvCB06FN+rrnKwQPfHJyoKffPmVO5zzx23Rwfu7MJs4iLjtJbhNkyJn0JuRS7bc7fXXxooBGz8N0R3sbtBkgTocgNEtFG/xjZQunIl5sJCIm652cHCPIeAuDiq9u3TWkaj8NjAXW4q52jJUelRYiOJCxJ5bO1jAAhE/aWBx7ZA7l7o/5C6Q5TYF70P9EuFomNQllvv5YWffY5vbBuvneDeGPzj49WhCjU1WktpMB4buA8UHUAg5I7bRqylgXpFdVD01/tfuTSwzbVw7xrodqMTVXoZiXfCwzvV+u46MM6eQ2bnBCrS0zEdOUpWl66yBNBG/OPj3HaogsceTmYXqrkrGbhtw1oaaBEWgCvXvQuh7rJb93aySi/Dx1/9f001VJddtjHHkDod0/FjlPzwI53WrkEfFuZkke5LwHkHlP7t3cs73mN33NmF2QT5BNE6xK2mrGmKtTRwcMxg9Iqe02WnL3/hkj/CyuecK85bMdfAvP6w8tnLPm3KzaX4+xVETJokg3YD8evQAfR6Kt0wz+2xgXtf4T46RXZCp3jsX9HuWO1eU3ulYhZmEltcJr+dmwl7loKf9Hh2Cnof9fB312IoPHzJ04WLFkFNDc2mTXW+NjdH5+eHf/t2VGXJwO0SCCFkRUkTiG8WT98Wfflk7yeXlgb+/IY6vf1aOU/SafRPBUUHG9++4NOWykqKPvuckKFD8Wsjx/I1Bv+4eCqzZeB2Cc6Un6G0ulQG7iYwtctU8iryLiwNLDgIvy2BpD82yQhJ0kDCWkHP22D7Aij9PX1V/PU3mIuKaHbnNA3FuTf+8fHUnDyFueSSMboujUcG7q2n5fCEppC4IJHU1anARaWByyeCzlfdAUqcy4BHwGKCnZ8BICwWCv73P/y7JBDUp4/G4twXd21998jA/WnWpwCsPbZWWyFuyiWT4K2ugSPmw8R5ENpCY4VeSLP2cN86uO5h1QWwS1eqDxygam8mWQldZAlgI7G2vrvbAaVHlQMmLkik2lxd+3jp/qUs3b8UP70f6Xeka6jMvbhkErzVNbBlT2jZU2N1XkzLHgA0f+A+ytavx1xURIcV30tDqSbgc9VV6MLD3e6A0qN23HJ4gv2wlgZOTVCrFQ7uWw5G93o76ZHsWcrZGddQuXs3UX/6kwzaTcA4ew5ZCV2wFBdT9MUXbuVf7lHfdUOQAX+9PwKBXtHL4QlNwOoaWG4q55usz/A5W8hdv/6VWSlz5NdTSwydyd9WhU9EBOETb9BajVtjSJ2OIXU6Z/71OoWffEJ8RjqKr3vMS/WoHTeo9qQAT/Z5kinxU8ivyNdYkXsTVF7AHYUFbA4KICN/b9PmUkqaTPnRs5Qb/YnqVIjOVKq1HI8goEsXhMlE1f79WkuxGY/acQMMix1GRm4GI9qOkDtDO5C4bCzVEaHA7xUmi/ctlucGGpH37n/QR4QREXsKNr0Nw1/SWpLbE9C1CwCVe/cSkJCgsRrb8Lgd9978vRgCDTJo24OiY6QdP8UYv6tqzaf89I2cSylpElYzqbMbNmAuKmHfkpZkpi7G+Ob/aS3N7fGLjUUXHEzlnj1aS7EZj9txZ+Zn0iWqi9YyPIPACAwDniCYfCyHVUP/anMj5lJKmkzz6Q9SvmULVQcP0nHVSnRVuVBZDK1klU9TUXQ6AhISqNyzV2spNuNRO+5yUzmHSg6REOUeb3dcHv9QGPwkBZYqpsRP4db4WwE4WHRQY2Hex9kNGynfupXmDzyALjgYmrX7PWhbLNqK8wACunalct8+t/Hm9qjAnV2YjUVY6NJM7ribzMrnYJ+6y7aaTz2S+AhRAVFU1VRx14q7bBtvJmkywmIh98038G3dmsgpky988tvH4BvZydpUArp2QVRWUnXQPTYl9QZuRVE+VBQlV1GU35whqCnszVff6sgddxM5tgV+mQOnL/yWB/kG8edr/szu/N1k5GbIChMnUZqWRtXeTAwPpaL4+V34pF8w7Fh4yfdK0jACunYF1ANKd8CWHffHwCgH67ALe/P30iygGVcFyWGpjUYI+OEFCI6Gfn+54KnEBYm88usr6mW2jDeTNBlhMmH899v4d+pE2Lhxl14w4FEICIOf/uZ8cR6EX9u2KEFBbpPnrjdwCyHWAwVO0NJkMgsySYhKQJEzEBtP9ko4ugn+8PQlntvWzlRfndqk4KvzlRUmDsQ4ew5Z3XtQfeQIVTk5ZHXtdmlnX1AzGPAY5KyEwxu0E+vmKHo9AZ07u82O26WqSl5avoe9Jxtnr2ihmmz//RQYO3Lzf36xszLvQC9MvG58HEXfmifS4zBnXPp1POVThklfAyiYzCY25BRxKDMHyHG6Xk8n2NSRVL9AjrfswIKJj144lPm8n3Ff0ZO3dM0xLnqKF6NmaaDUMxhtiaTXng3cMm8jQte4478urcJ44fqudlZ2KXY7nFQU5T5FUbYpirLNaDTaa1mbOavsA8WCr0WWqTUWC3qWhdzCh+EPYlYu/zu9Rikl0jyIVqZ7ADiry+Sw3+vUUOxMqV5BysYv8a2pZsUfbr0waF+ESfHn7chneDviKSeq8zxORrfFz1RFVNEZraXUi9123EKI94D3AJKSkkRj1mjKb6q7VvyHY7nQOz6fVwf1a/Q6kuvqef73r+1LvxSxJHsJZgrpec1Wnk9+3rHSvIiKXbs4/OYGmt19N3OfnGTDK859X4QAswl8/K58ueQSKvdFcWjVB7zSzZfw6107hihC1B9jFUVpC3wrhOhmy6JJSUli27ZtTVNmIxdbuVqRLdkNZOVzEBEL195n0+Xy6+44hMXC4VtvxXTyJB1WrEAfEmLbC801sGASXNUNRr3iWJEeiKipYV9iEpG33spVTzv/3YuiKOlCiCRbrrWlHPBT4BcgXlGU44qi3NNUgfZEWrnagZPb4Zd3oOiIzS+5eNiCPKi0D9YhCZU7d2E25pGd1Md2q1G9D0S2hS3/AaN7+UtrjXH2HLK6dUdUVVHw8ccub/Fab6pECHGrM4Q0FkOQAT+9n7RybSwWs9rEEdwcBj9p88vOH7agoGCymNDr9PLr3kTCJ04k/6OPCOrVi6v/+37DK6RS/gp7lsF3j8Ody6+YG5f8zsUWr3HbtqLz99daVp14ROfkoeJDAMxMnimtXBvKlvfgZAaMehUCwhv0UuuwhTeHvIkOHRtObOCuNNlR2ViEEJz+619RgJZ/e6lxZa3BzVXHwMM/q405kgYRlNgbYTK5vOGURwTuvi36olf0jGk3hpnJM2uHAEjqoaIQVr8MHYdDtxsb/HJrK3xKmxQe6v0QBZUFZJyRHZWNpfjLLzm7aRPRM57At3Xrxi/U+05o0x9+fVc9rJTYTGCvXgBUZGRorOTKuFQdd2PZnrudzs06E+QbpLUU9yIwEm77HCLaNOkt9fkHldKzu3GYTp3izKuvEdSnDxE339y0xXQ6uPG/6jsomSppED7NmuHXrh3l6RlE/UlrNXXj9jtuk9nE7rzd9IrupbUU96KqTP1/2wFq4G4C1oNKP71agqZDx5h2Y+RBpY3k/vtt9g8ZiqWsjPKtW8nq0rXpB2PhrcE/BEyV8qCygQQm9qYiIwPhwq6Lbh+49xbspcpcRe+remstxX0oOQVv94Qt79tlOetBpclswkfxwYKFg8UHeWLdEzLfbQvnAkSrWbNIyMqs/TCkTm/62l/9CeaPh3K3cK1wCYJ69cZcXEy1CzsFun3g3pG7A0DuuG3FYoGv/6LuuNv/wW7LWg8qPx37KTEhMWQVZMl8tw2c3bSJ/PfeI/ymGwkfN9b+Nxg0A8rzYflDMt9tI0GJ6iaw3IXz3G4fuDPOZHB16NWyDM1WtrwHB1bDyH9A8052W9Z6UHn7its5XnYckA6C9WHKzeXEk0/h16E9LZ57zjE3aXkNDJ0JmcthxyLH3MPD8I2NRR8VRUW6DNwOQQjB9tztcrdtK7mZ8MNfIW4UJN3tkFvU5rt1ar5bQWF029Ey330e1vmR+wcNxpyXR/X+A+zr1dtxzR79U6HtQFjxJBQccsw9PAhFUQjq3UvuuB1FRm4GhVWFdIqw387Ro8ndq9b5jp/jsGqD2ny3Rc13CwS/5f3GE2tlvttK8wf/QujoUaAoxMx9x7457cuh08MN89Tdt8QmAnsnYjp2DFNurtZSLotbB+53tr8DwJ581y6Wdxm63QjTt0KIwaG3sea7Pxv3GV2iunCs7Bjpueky332OvDnvULoijegnHid06FDn3DTiarjrO3VWpRAy310P1jx3RcZ2jZVcHptMphqKo02mpMFRA/n1PxByFXS9wam3ld+nSyn68ktOPTeT8EmTaPmPl50/9MNUAUvvVw+mk/7o3Hu7EcJkYl+fvkRMmUyLZ591yj3tajLliljzqFaksdQVOLgO0p6GPUudvsuyfp/89b97PiRGJ7JozCKva4235rVPPTcTgOKvviIroYvzTYz0flBdBt/PgGNbnXtvN8E6eUhUVlL4v09c0nDKLQO3IciAyWICwEfxkcZSdZGbCYunQvM4mOC4vHZdWPPd1ebq2sPKjNwM3s542+tKBf07x4NeT1CfPsRvz3B8XrsudHqY9L7aoPPZrVDgurXKWmFInU5CViZXPfM0AB1WrdTme3UF3DJwA+wv2g/AvOHzpLHU5Sg+DgtuBJ9AuP0L8A/VRIY1371o7CJ0ig6BYP2J9V5VKlj644+ceOxxArt1I2bePHSBgdoKCmoGty9RnSE/mQRlrnkApzUh584fSlev1ljJpbitV0mAPoCehp4kt0wmuWWy1nJcjz3LoKoU/riiyS3tTeF8w68fb/qRV359hZ+O/oRA4KvzZVDrQRgrjeRV5HncOybj7DnkvfNO7eOKnTvJTkqi+YMPar97a94JblsMX9wFJScgJFpbPS6I39VX49+pE2U/rSbqrru0lnMBbrnjPlV2isyCTIa0GaK1FNel/3T4yy/QwqahRU7BEGQgMiASoNbDe4dxB7uNuz0ubSKEAJ2amgoeNJD4jHTt0iN1cXUfeCgDWp3rg6i59CDZ2wlJGUp5ejo1hYVaS7kAtwzca46tAWDI1TJwX0B5AXwyEc6cK48Mj9FWz2Wwpk6sk3PyK/MvSJv0mN/D7Q8tLRUVnHr6GfJmzyF8wgSufucddEEu6lzpc+7g+Of/g/nXQ2WJtnpcjNCUFLBYKFu3TmspF+C2gbttWFvahbfTWorrUJYLH4+DwxtVEykXxdoav/LGlYxuOxq9ogfUHXir4FYAbr37Pv2Pf7CvV2+Kv/4agOKvvyarew+Xqki4LM3aw4lt8L8J0pDqPAK6dsUnOpqyn1wrz+12Oe5DxYfYfGozN8c10TlJyUMAAAyCSURBVLPYkzBmw6c3Q+lpuH2xXc2jHIUhyECIXwgWoTrjCQQnz54EqPXzVlBYPWW12+S+S1auovjrb9CFhdH69X8RMniw1pJsp+tE8AmAxdPgozFw66dqs46Xo+h0hAwZQvHy5ViqqlxmnJnb7bhf+VWdXm2sMGqsxEU4tQv+m6IeRE5d5hZB24o1bfL+8PdpGdSy9vN6RU9MiJrmcYfd95nXXyezcwInHn4YS0kJlpISjv35ftffZV9M/Gi12qT0FHwwQv2ZkhCaMhRRXk755s1aS6nFbTonZRdeHdRUq+ZBAx9X25rdlL/98jeWZC9BcPmfR1fcfQshKPn2O868+irmkhIMf3mAqD/9CcXXV2tpTSP/AJxIhx5T1MdCePUkHUt1NTnJ/QgbN46Wf3vJYffxyM7JtElppFydUvvYq7sli47BkrvVXKSPH1z/llsHbbhw990iqEXt5xUUWgSrj+ftnIex3OgSXZcnn5tJVkIXTs6YgTk/H0wmjP9+m7x3/6OpLrsQ1eH3oJ35LSyaoqbhvBDj7Dns63ENlvJyihYvdpkuSrfJcRuCDBwrOwaAn87PO7slhYCdn8KKp0BYoPc0t0qNXInz670Hxgys3X0LBKfPqkHDmvsGNYg/n/y803VW7NlD/rv/ofSHH9AbmhP96GOE3zABRec2e6CGUVEIh9bD3H4wdhZ0neRVu29D6nQMqdOp3LePQxNuoHnqdAwPPqi1LPfZcRdWFrK/aD+xobEsGrvI+7olj2+DD4bDsgcgugvcv8FjgvbFnL/7jgmJQeHSQOHM8kFhsXD2l1/YP2Ikh2+8idIffgDAbMzj1LPPkvfOXIfeX1N6T4U//wyRseq7vI/HwendWqtyOgHx8YQMHUrh/z7BXHZWaznuEbiN5UamfDsFi7Dw76H/Jr5ZPDOTZ16wS/N4Nv4bio7ChLlqN6QHn/hbSwaTWyXTr1U/AHyVC/PGCgqR/mozz8UpFHulU2ry8zly9z1kdenK0T/ejeno0drnmj/4oOs11DgKQxzc8yOM/T/V0906fNjLrGGb3/9nzMXFFH3+mdZS3ONw8q8b/8rS/UuJCYlhxY0r7LauyyIEHNkEm96GlL/CVV3VOm3fQM08R7TikTWP0DywOZPjJvPU+qc4UHzgitdPiVdzs1/s+4LJ8ZO5v8f9zFg/g1mDZ9mUVqvJz6f0p58oWbGC8l+3gMVCYK9eRN5yM6EjR6ILCLDL38ttqSwG/zA1XbL+dTi1E/o/DDFJXpFCOXr3PVTu20fHn360+89CQw4nXTpwe10lSdFR2Pm5mscuOACBzWD8bEgYp7Uyl8AaxIe1GcbfN/+d46XHsWCx6bVT4qdcEMSFEMxYP4PX+/yN6l17SfvqdQYdD6EmM7vONVzCY8SV2Pg2/DxLDebN46HnrdB9iuo86KGUb93KkanTuOq552g29Q67ru0RgdtYbuSRNY+gKAo7jTsBtZIkpU0KT/R5wjMOJU0VUHYGItuqTm3/aqf+I2g7EK65VR184BestUqXxFo+6KPzwWQxoVN0tc08deFnElxthJuURMKPFtLs1xya1ZOu3DCyNRNf/aw20J8f9C/+s0f8TDaUqlL47St1EPGxzeqUpZs+VJ87nq565fi4RtNKU7nYNMyKvX6h270cUFGUUYqi7FMUZb+iKE83TV7dnJ+bfHfXu+zO210btD2ikuT0btj2EXz/JLw3BP4ZAwsnq8/p9HDDu/DwLrjrW+h1uwzaV8B6gPnp2E/pEN4Bi7AQLHwxFAnijgsG7IUJv1i4J83Mc5+aeeedGhbMMvPP+WYSP95Cx9UXBu2lyQpTnvGp/Vjy8W0s+fg2ZvfOZd7Oeby7691aD/G6/lxXnt2WPwNNer1mfzZXcpdxDXm3LcR474/c5Vemfv7EVu769hbyXmuD8b9DuGvRYPI2z8V4eqfbfi1mdNqGYfvPRKZ9SV6LQPDzI/TNfzKj0za7nq3YQr07bkVR9EA2MBw4DmwFbhVC7K3rNY3dcf99899ry70uxkfnw42dbiSvIs95h5JCqGV3ik7N35kq1ekhNVVgrlJ3zKYKaNkT9D5q5cexLVCeB2eNal667Azcu0Z9/dIHYOci8A1WHdmu7gMxfdWONQ/KDwohwGIBi+X3P5vNCOv/z31gNiNqahCmGoTJhKgxIaqrEVXVCFM1lspKRGUVoqoSS3m5+nG2HHNZKZaSUsylJZiLizlzMofAshr8KkwN0vnFAIUvBurt/vcf2XYkOnSkHU5jXPtxKCgsP7icCR0noKCwbP8yJnWahILClzlfclPcTTyf/Dwv//pybW4ecO8/m018sX8Zk4NioczIF0oZk0vLoMNQvjBuY3LLQfD/7Z1bbBRVGMd//267BQk3AUEEuQWIBB8kRDAmqIEoISoPsqYmeIkoAYMP+maICZEHookkmpBoH4iXBITyoDVeSFQISihqAnJLVMrNKpcGEUUf2m0/H2aKpXTplM7OOPD9ksmes+fs7v8/Z/bbmXPmnD2xi7qqFgr50VDVn7q/GylMeJCX56xh9Y6V1B39mMLYeVCRo+74VgoTH4JcNXU/b6EwZVHweT9tScXn53s281r9IIYcP8fGOWL03XP5c/RgNjR9RGFq4apuVY21q0TSXcAqM3sgzL8EYGZrSr2mt4G7oy+79o0i1cXLyytVQT5Xjaw9CJ5dqeoHyoG1lSjvHwTf9mIQdEuWt5YoHxAE1rbWIGB3JT8AELS1BBuEgVjB+1Z2DGJ07Otkg3SkzrDOx8GV0mHeOuc7b/8ThtTUsHZWM4OHjrpkYDOfy18cN+mczilHm7UhVHL2puN0UPi6jcI3lx8nF6rh6RdyIPV6LC7uwL0ImG9mz4T5x4FZZraiS72lwNIwOxX4MargUVWVY4ZVVI6MWt+5PjjbVjx5qlj8LY73yo/MT7I2ay3+VWzOj8hPAmhpbmm8JD08P0FV6o9hF28e702662P017dbuxVVoUpExVV//rWRzty+uPECDO1mWZez7cXTp1qLTZeXlGScmY2IUjG2mZNmVgvU9vV9JH0f9VfnWsE9X/tcb37BPZeTKIOTvwKdF8IYEz7nOI7jpECUwP0dMFnSBEl5oAaoL68sx3EcpxQ9dpWYWVHSCmArkAPWm9nBMmrqc3dLBnHP1z7Xm19wz2WjLBNwHMdxnPKRiUWmHMdxnP/wwO04jpMxUgvcPU2jl1QtaVNYvlvS+ORVxkcEvy9KOiRpn6QvJY1LQ2ecRF0qQdIjkkxS5m8di+JZ0qNhWx+UtCFpjXET4di+VdI2SXvC43tBGjrjQtJ6SWckHShRLklvhvtjn6QZsYsws8Q3gkHORmAikAd+AKZ1qfMc8FaYrgE2paE1Qb/3ATeE6eVZ9hvVc1hvILADaABmpq07gXaeDOwBhob5m9LWnYDnWmB5mJ4GHEtbdx89zwFmAAdKlC8APiOYIj0b2B23hrTOuO8EDpvZETNrAT4AFnapsxB4N0xvAeZKmV3Qo0e/ZrbNzP4Jsw0E98tnmShtDLAaeBXoZq2CzBHF87PAOjM7B2BmZxLWGDdRPBswKEwPBmKZDZsWZrYD+P0KVRYC71lAAzBE0s1xakgrcN8C/NIp3xQ+120dMysC54FhiaiLnyh+O7OE4Bc7y/ToObyEHGtmnyQprIxEaecpwBRJOyU1SJqfmLryEMXzKmCxpCbgU+D5ZKSlRm+/770mM38WfL0gaTEwE7gnbS3lRFIFsBZ4KmUpSVNJ0F1yL8FV1Q5Jt5vZH6mqKi+PAe+Y2evhonXvS5pu1sMC6k5J0jrjjjKN/mIdSZUEl1hZ/XfgSMsGSJoHrAQeNrNuliHMFD15HghMB7ZLOkbQF1if8QHKKO3cBNSbWauZHSVYMnlyQvrKQRTPS4DNAGa2C+gHZHRR/UiUfZmQtAJ3lGn09cCTYXoR8JWFPf8ZpEe/ku4A3iYI2lnv94QePJvZeTMbbmbjzWw8Qb/+w2YW35+VJk+U4/pDgrNtJA0n6Do5kqTImIni+QQwF0DSbQSBuzlRlclSDzwR3l0yGzhvZidj/YQUR2YXEJxtNAIrw+deIfjyQtC4dcBh4FtgYtqjyWX2+wVwGtgbbvVpay635y51t5Pxu0oitrMIuogOAfuBmrQ1J+B5GrCT4I6TvcD9aWvuo9+NwEmgleAKagmwDFjWqY3XhftjfzmOa5/y7jiOkzF85qTjOE7G8MDtOI6TMTxwO47jZAwP3I7jOBnDA7fjOE7G8MDtOI6TMTxwO47jZIx/AVa/T7BICNDvAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt \n", "\n", "plt.plot(x, y1, label = 'Beta(1, 1)')\n", "plt.plot(x, y2, '--', label = 'Beta(10, 10)')\n", "plt.plot(x, y3, '-*', label = 'Beta(4, 16)')\n", "plt.plot(x, y4, '_-', label = 'Beta(16, 4)')\n", "plt.legend(loc=0,numpoints=1,fontsize=10)\n", "plt.ylim(0, 6)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "So let’s say we assume a prior distribution on p. \n", "- Maybe we don’t want to take a stand on whether the coin is fair, and we choose alpha and beta to both equal 1. \n", "- Or maybe we have a strong belief that it lands heads 55% of the time, and we choose alpha equals 55, beta equals 45." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "Then we flip our coin a bunch of times and see h heads and t tails. \n", "- Bayes’s Theorem tells us that the posterior distribution for p is again a Beta distribution \n", " - but with parameters alpha + h and beta + t.\n", " \n", " " ] }, { "cell_type": "markdown", "metadata": { "ExecuteTime": { "end_time": "2018-10-27T04:34:36.368214Z", "start_time": "2018-10-27T04:34:36.363726Z" }, "slideshow": { "slide_type": "skip" } }, "source": [ "Beta is the conjugate prior (共轭先验) to the Binomial distribution. \n", "\n", "Whenever you update a Beta prior using observations from the corresponding binomial, you will get back a Beta posterior." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "Let’s say you flip the coin 10 times and see only 3 heads.\n", "\n", "- If you started with the uniform prior, your posterior distribution would be a Beta(4, 8), centered around 0.33. \n", " - Since you considered all probabilities equally likely, your best guess is something pretty close to the observed probability.\n", " \n", " " ] }, { "cell_type": "code", "execution_count": 61, "metadata": { "ExecuteTime": { "end_time": "2018-10-27T04:57:45.304161Z", "start_time": "2018-10-27T04:57:45.163044Z" }, "slideshow": { "slide_type": "skip" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl4VdXd9vHvL3MCgTAFCAECyAyCGBREBa0CTtA61TqgFUWtWvtoq7U+Vm37vFiHOlSxUrVOONd5qKhMIiAGCPOMgAEJSYCEQBIyrPePxBRkyAk5Jzv75P5cF5dJzs7h3jnhdp+1117bnHOIiIh/RHgdQEREakfFLSLiMypuERGfUXGLiPiMiltExGdU3CIiPhNQcZtZkpm9ZWarzGylmQ0NdTARETm0qAC3ewz4j3PuQjOLARJCmElERI7AaroAx8yaA5lAV6erdUREPBfIEXcXIAf4l5kNABYAtzjn9uy/kZlNACYANGnS5PhevXoFO6uISNhasGBBrnOuTSDbBnLEnQ7MA4Y55742s8eAAufc3Yf7nvT0dJeRkVGbzCIijZqZLXDOpQeybSAnJ7OALOfc11WfvwUMOtpwIiJSNzUWt3NuG/CdmfWs+tJPgBUhTSUiIocV6KySm4EpVTNKNgC/DF0kERE5koCK2zmXCQQ09iIi4aW0tJSsrCyKi4u9jhIW4uLiSE1NJTo6+qifI9AjbhFppLKyskhMTCQtLQ0z8zqOrznnyMvLIysriy5duhz18+iSdxE5ouLiYlq1aqXSDgIzo1WrVnV+96LiFpEaqbSDJxg/SxW3iIjPqLhFRHxGxS0iDV5kZCQDBw5kwIABDBo0iDlz5hxx+127djFp0qSAnruoqIjhw4dTXl5e/bWCggJSU1O56aabavz+zMxMhgwZwsCBA0lPT2f+/PkAfPjhh/zxj38MKENtqbhFpMGLj48nMzOTxYsXM3HiRO68884jbl+b4n7uuec4//zziYyMrP7a3XffzamnnhrQ999+++3cc889ZGZm8qc//Ynbb78dgHPOOYcPPviAvXv3BvQ8taHpgCJSKz9/eu5BXzv32PZcMTSNon3lXPWv+Qc9fuHxqVyU3pEde/Zxw8sLDnjs9etqt7x/QUEBLVq0qP78wQcf5I033qCkpISf/exn3Hffffz+979n/fr1DBw4kDPPPJN77rmHsWPHsnPnTkpLS/nLX/7C2LFjAZgyZQqvvPJK9fMtWLCA7OxsRo8eTSBrLpkZBQUFAOTn55OSklL99REjRvDhhx9y8cUX12ofa6LiFpEGr6ioiIEDB1JcXMz333/PtGnTAJg6dSpr165l/vz5OOcYM2YMs2bN4v7772fZsmVkZmYCUFZWxjvvvEOzZs3Izc1lyJAhjBkzhtLSUjZs2EBaWhoAFRUV3Hbbbbz88st8/vnnAWV79NFHGTVqFL/97W+pqKg4YBgnPT2dL7/8UsUtIt460hFyfEzkER9v2SSm1kfY8N+hEoC5c+cybtw4li1bxtSpU5k6dSrHHXccAIWFhaxdu5ZOnTod8P3OOf7whz8wa9YsIiIi2LJlC9nZ2VRUVJCUlFS93aRJkzj77LNJTU0NONtTTz3FI488wgUXXMAbb7zB+PHjq0s/OTmZrVu31np/a6LiFhFfGTp0KLm5ueTk5OCc48477+S66647YJuNGzce8PmUKVPIyclhwYIFREdHk5aWRnFxMc2bNz/gYpi5c+fy5ZdfMmnSJAoLC9m3bx9Nmzbl/vvvP2yeF154gcceewyAiy66iGuuuab6seLiYuLj44Ow1wfSyUkR8ZVVq1ZRXl5Oq1atGDVqFM899xyFhYUAbNmyhe3bt5OYmMju3burvyc/P5/k5GSio6OZPn06mzZtAqBFixaUl5dXl/eUKVPYvHkzGzdu5KGHHmLcuHHVpT1u3LjqGSP7S0lJYebMmQBMmzaN7t27Vz+2Zs0a+vXrF/SfgY64RaTB+2GMGyqHPV544QUiIyMZOXIkK1euZOjQyuGXpk2b8vLLL9OtWzeGDRtGv379OOuss7jjjjs477zz6N+/P+np6ex/h66RI0cye/ZszjjjjCNmWLJkSfWJx/3985//5JZbbqGsrIy4uDgmT55c/dj06dOZOHFiMH4EB6jxDjhHQ3fAEQkfK1eupHfv3l7HCJmFCxfyyCOP8NJLLx12m4KCAsaPH8+bb74Z8PNmZ2dz6aWX8sUXXxz02KF+psG+A46ISNgaNGgQp5122gEX4PxYs2bNalXaAJs3b+bhhx+ua7xD0lCJiNTIORfWC01dffXVQX/OwYMHH/LrwRjl0BG3iBxRXFwceXl5QSmcxu6H9bjj4uLq9Dw64haRI0pNTSUrK4ucnByvo4SFH+6AUxcqbhE5oujo6DrdrUWCT0MlIiI+o+IWEfEZFbeIiM+ouEVEfEbFLSLiMypuERGfUXGLiPiMiltExGcCugDHzDYCu4FyoCzQFaxERCT4anPl5GnOudyQJRERkYBoqERExGcCLW4HTDWzBWY2IZSBRETkyAIdKjnZObfFzJKBz8xslXNu1v4bVBX6BOCgOyyLiEjwBHTE7ZzbUvXf7cA7wAmH2Gaycy7dOZfepk2b4KYUEZFqNRa3mTUxs8QfPgZGAstCHUxERA4tkKGStsA7VbctigJecc79J6SpRETksGosbufcBmBAPWQREZEAaDqgiIjPqLhFRHxGxS0i4jMqbhERn1Fxi4j4jIpbRMRnVNwiIj6j4hYR8RkVt4iIz6i4RUR8RsUtIuIzKm4REZ9RcYuI+IyKW0TEZ1TcIiI+o+IWEfEZFbeIiM+ouEVEfEbFLSLiMypuERGfUXGLiPiMiltExGdU3CIiPqPiFhHxGRW3iIjPqLhFRHxGxS0i4jMqbhERnwm4uM0s0swWmdmHoQwkIiJHVpsj7luAlaEKIiIigQmouM0sFTgHeCa0cUREpCaBHnE/CtwOVBxuAzObYGYZZpaRk5MTlHAiInKwGovbzM4FtjvnFhxpO+fcZOdcunMuvU2bNkELKCIiBwrkiHsYMMbMNgKvAaeb2cshTSUiIodVY3E75+50zqU659KAS4BpzrnLQ55MREQOSfO4RUR8Jqo2GzvnZgAzQpJEREQCoiNuERGfUXGLiPiMiltExGdU3CIiPqPiFhHxGRW3iIjPqLhFRHxGxS0i4jMqbhERn1Fxi4j4jIpbRMRnVNwiIj6j4hYR8RkVt4iIz6i4RUR8RsUtIuIzKm4REZ9RcYuI+IyKW0TEZ1TcIiI+o+IWEfEZFbeIiM+ouEVEfEbFLSLiMypuERGfUXGLiPiMiltExGdqLG4zizOz+Wa22MyWm9l99RFMREQOLSqAbUqA051zhWYWDcw2s0+cc/NCnE1ERA6hxuJ2zjmgsOrT6Ko/LpShRETk8AIa4zazSDPLBLYDnznnvj7ENhPMLMPMMnJycoKdU0REqgQyVIJzrhwYaGZJwDtm1s85t+xH20wGJgOkp6friLwBKq9w5OwuIW9PCXmF+yivcHRqlUC3Nk29jiYitRBQcf/AObfLzKYDo4FlNW0v3ispKyc2KpKKCsew+6exraD4gMfHDEjh8V8ch3OOa19cQI+2TTn5mNYM6dqKiAjzKLWIHEmNxW1mbYDSqtKOB84E/hryZHLUnHPMXpfL819tJKewhPduHEZEhDHh1K7EREXQumksrZrGEBlhNI2t/BUoKC4jt7CEGau3M2nGelJbxPPz9I5clN6Rds3jPN4jEdlfIEfc7YEXzCySyjHxN5xzH4Y2lhytRZt38sf3lrN0Sz6tm8Zw6YmdKatwREcaV5/c5bDf1zw+mndvHMbefWV8vnI7r3+zmYc/W0OXNk0499iUetwDEamJVU4aCa709HSXkZER9OeVI5u1Jocr/zWf5MRYbhvZk7EDU4iNijzq5/tux146JMUTEWE8O/tbCovLuH5E1zo9p4gcmpktcM6lB7Jtrca4pWEqLi0nLjqSod1a8btRPRk3NK16CKQuOrZMqP549bYC3sjI4qOlW3nwwgEM6JhU5+cXkaOjS959zDnHS3M3csbfZpKzu4ToyAh+NeKYoJT2jz1w4QD+9cvBFBSV8bNJX/HAf1ZRWl4R9L9HRGqm4vapsvIK7vj3Eu5+bzk92iYSExn6l/K0nslMvfVULjw+lX/MXM+SrPyQ/50icjANlfjQvrIKbnltEZ8s28avTz+G35zRo96m7jWLi+aBCwcw4dSuHJOcCMD2gmKSm2nmiUh90RG3Dz32xRo+WbaNP57bh1tH9vRkvvUPpT1nXS4nPzCddxdtqfcMIo2Vjrh96Prh3eib0pyz+7f3Ogp9OzTn+E4t+M3rmWwrKOa6U7tipgt3REJJR9w+UVHheHb2t+zdV0ZiXHSDKG2onP/9/NWDOW9ACvd/sor7PlhBRYVWPBAJJR1x+8RfP13F0zM3kBgbxcWDO3od5wCxUZE89vOBtE2M5ZnZ3zKka0tG92sY/2MRCUcqbh94Yc5Gnp65gcuHdOKi9FSv4xxSRIRx1zm9GdEzmWHHtPI6jkhY01BJAzdnfS73fbCcM3q35b4x/Rr0+LGZcXL31pgZq7ft5t73l1OuYRORoFNxN2AVFY573ltOl9ZNePSSgUT6aLW+2etyeX7ORu56ZymhWFZBpDHTUEkDFhFhPH/1CRSXlofkashQGn9yF3bu2ccT09fRPD6a35/Vq0G/WxDxE3+1QSMye20uJ3VrRYekeK+jHLXbRvagoLiUp2dtoFl8NDeedozXkUTCgoZKGqDpq7dz+bNf8+o3m72OUidmxr3n9eWnA1OYsz5Xa5uIBImOuBuYHXv2cftbS+jZNpELBjXMGSS1ERFhPHjRAJyD6HpYT0WkMdC/pAbEOcf/vruUXXv38cjPBxIXHR7rXkdHRhATFcGOPfu44tmvWbZFi1OJ1IWKuwF5N3MLHy/dxv+c2YM+Kc28jhN0peUVrN9eyNXPf8PWXUVexxHxLRV3A9K+eTxjBqRw3andvI4SEm2bxfH81SdQtK+ca17IYO++Mq8jifiSirsBGdK1FY//4jhfzdeurR5tE/n7pcexalsB//N6ptY1ETkKKu4G4PMV2Uz8ZCXFpeVeR6kXI3omc9c5fVi1bTe5e0q8jiPiOypujxWWlHH3e8uYsSonrI+0f+zqYWl8/OtTSE7UDRhEakvF7bGHPl3NtoJi/t/5/RvVdDkzo0lsFCVl5fzxvWUs36qZJiKBajxN0QAt25LPC3M3cvmJnTm+cwuv43hid3EZU5dnc91LC9i5Z5/XcUR8QcXtob98tIKWCTH8dlRPr6N4pnXTWJ6+4ni27y7hplcXUqarK0VqpOL20J/H9uOhiwfQPD7a6yieGtAxif/7aT++WpfHg5+u9jqOSIOn4vbAD1PgurdN5LSeyR6naRguSu/IFUM6M+XrzWzfXex1HJEGTcXtgYemruamVzQs8GN3n9uHD28+WTNNRGpQY3GbWUczm25mK8xsuZndUh/BwtXmvL088+W3xERFENWIZpEEIiYqgrTWTXDO8daCLApLdGWlyKEE0hxlwG3OuT7AEOBGM+sT2ljh66+friIywrhjdC+vozRYa7ILuf2txdz5tu6eI3IoNRa3c+5759zCqo93AyuBDqEOFo4Wbd7JR0u+59pTutC2mYYDDqdnu0RuG9mTDxZv5aV5m7yOI9Lg1Oq9upmlAccBXx/isQlmlmFmGTk5OcFJF2Ye/2ItrZvGMmF4eC4iFUw3DO/G6b2S+fOHK1i0eafXcUQaFAv0raiZNQVmAv/nnHv7SNump6e7jIyMIMQLL3mFJWzI3cPgtJZeR/GFXXv3ce7fZ+McfHHb8LBZn1zkUMxsgXMuPZBtA7oDjplFA/8GptRU2nKw8gpHhEGrprG0ahrrdRzfSEqIYdJlg8gr3KfSFtlPILNKDHgWWOmc+1voI4Wf177ZzE8nzWHXXl3SXVvHpiZxWq/Kue7ZBZrfLQKBjXEPA64ATjezzKo/Z4c4V9goLi3n8S/WEhVhjf4Kybr4fEU2pzwwnXkb8ryOIuK5QGaVzHbOmXPuWOfcwKo/H9dHuHDwwpyNZBeUcPuonlS+eZGjMaRbK1KT4rn51UW6slIaPV0BEkL5RaVMmrGe4T3acGLXVl7H8bWmsVFMunwQu4tL+c1rmZTrzjnSiKm4Q+jleZvILyrld4149b9g6tWuGX8a24856/N47Iu1XscR8UxAs0rk6FxzShd6t0+kX4fmXkcJGxend2Thpp0kxGiWiTReKu4Qcc4RGxXJ6b3aeh0l7Ew8v7/OF0ijpqGSEMguKOYnf5vJ3PWaAREKP5T2jNXbmfBihlZZlEZHxR0Ck6avY3PeXjokxXsdJazlF5UydUU2D07VzRekcVFxB9mWXUW8Ov87LkpPpVOrBK/jhLWxAztw2YmdeHrmBj5bke11HJF6o+IOsiemrQPgptO7e5ykcbj73D7069CM297I5Lsde72OI1IvVNxBtGVXEW9mfMclJ3TUMEk9iYuOZNKlx+OAtxZkeR1HpF5oVkkQpTSPY9JlgxjQMcnrKI1Kp1YJfPzrU0htof9ZSuOgI+4gMjNG9m2nmyR4oGPLBMyM9TmFGu+WsKfiDpL7PljOE9N0NZ/XJn68kptfXciqbQVeRxEJGRV3EGzK28OLczeRt0fLtnrt/53fn2Zx0Vz/0gIKiku9jiMSEiruIPj7tHVERRg36JZknktOjOPJywaRtbOI295YTIUWo5IwpOKuo425e3hn0RYuO7EzyRrbbhAGp7XkrnN689mKbF7P+M7rOCJBp1kldfTE9Mqj7etHdPU6iuznqpPSiI+O5KfHdfA6ikjQqbjr6MqhaZzQpSXJiTrabkjMjEtO6ARUXhpfWFKmufUSNlTcddQ/tTn9U7Vsa0PlnOOqf82naF85b//qJBJi9Csv/qcx7qO0MXcPv31zsW5g28CZGf9zRg/WZO/md28twTmdrBT/U3EfpSemr+ODxVvRstAN36k92nDH6F58tOR7npq53us4InWm4j4KP8wkuXxIZ41t+8SEU7ty3oAUHvx0NdNW6cpK8TcN+B2FH+ZtXzdcM0n8wsx48MJjaRobSb8UnZMQf1Nx19KGnELeWZTF1cO66GjbZ+KiI5l4/rEAlJVXsLe0nGZx0R6nEqk9DZXUUtO4KK4Y0pnrdJWkbznnuPGVhVzzfAYlZeVexxGpNRV3LSUnxnHf2H60SYz1OoocJTPjnGNTmL9xB7//91LNNBHfUXHXwlMz1rNg0w6vY0gQjBmQwm1n9uCdRVt4/It1XscRqRUVd4DWZu/mgU9X8fnK7V5HkSC56fRjuGBQKo98voZ3F23xOo5IwGo8OWlmzwHnAtudc/1CH6lheuTzNSRER3LtKZpJEi7MjInn96ekrJxubZp6HUckYIEccT8PjA5xjgZt2ZZ8Pl66jfGndKVlkxiv40gQxURF8MSlg6qXLdihNdXFB2osbufcLKBRD+w+PHU1zeOjueaULl5HkRD6x8z1jH50lu4WLw2exrhr4JxjcJeW3Dayh+b8hrnTeyVTXFrOlc/NJ7ewxOs4IocVtOI2swlmlmFmGTk5OcF6Ws+ZGb8acQzjhqZ5HUVCrEfbRJ69ajBb84u48rn5uvWZNFhBK27n3GTnXLpzLr1NmzbBelpPZWzcwXuZW3T7q0ZkcFpL/nH58azJ3s01L2TotZcGSZe8H0ZFheO+D1aQV1jCqL7tiIuI9DqS1JMRPZN59OfHUVZRQUSEln+UhqfGI24zexWYC/Q0sywzGx/6WN77cOn3LN2Sz20jexIXrdJubM45tj1jB1be9mzh5p26NF4alBqPuJ1zv6iPIA3JvrIKHvp0Nb3aJeqehY3c1l1FXDJ5HsN7tOHJSwcRE6Xz+eI9/RYewitfb2Lzjr3ccVYvIvVWuVFLSYrnf6vuGH/zqwvZV1bhdSQRFfehpCTFc9HxqYzoER4nWaVuxg1N497z+vDp8myuf3kBxaUaNhFv6eTkIYzs246Rfdt5HUMakKuGdSEmKpK73l3KmwuyuGJIZ68jSSOm4t7Pdzv28v7irYw/uYtOSMpBLj2xEz3bNWVQpxZeR5FGTkMl+5n4yUqemLaOnXu1XoUc2vGdW2JmbM7by2XPzGPrriKvI0kjpOKuMm9DHh8v3cYNI7rRvnm813GkgcveXcyS7/K54Kk5rM3e7XUcaWRU3EB51cU2HZLimXCqlm2Vmg1Oa8lr1w2hrMJxwVNzmLMu1+tI0oiouIHXvtnMyu8L+MPZvTW2LQHrm9Kct284ibbN4hj33HxmrgmfNXqkYVNxAwNSk7jqpDTO7q+ZJFI7HVsm8O9fncS4oWmkd9ZJS6kfKm6gX4fm3DumL2a62EZqr1lcNH88rw9NYqPYU1LGnW8v1bKwElKNurhnrcnhd28uZreW75QgWZy1i38vzOK8v89m0eadXseRMNVoi7toXzl3vbuUBZt2av0JCZqTurXm7RtOIjLCuPjpuTzz5QYtDStB12gb65HP1/DdjiL+72f9iY3SCUkJnn4dmvPhzSczomcyf/loJQ98utrrSBJmGuWVk/M25PHPLzfwixM6MbRbK6/jSBhKSohh8hXHM+XrzQyvWvOmpKxcBwkSFI2uuJ1z/OmDFXRumcD/ntPb6zgSxsyMy6vWNHHOceOUhTSJjeLe8/rSokmMx+nEzxpdcZsZz1yZTn5RKU1iG93ui0cqXOUQyhPT1vHVujzuPrc3YwakaCaTHJVGNca9IaeQigpHSlI8vds38zqONCKREcZvzujBezcNo33zOG55LZPLn/1aa53IUWk0xb02ezfn/X02f/3PKq+jSCPWN6U57944jD+P7cv3u4ppElP5rs85zTyRwDWK4s4vKuXaFzNIiI3il8O6eB1HGrnICOOKoWl8dutwmidEU1Zewc8nz+OZLzfoJg0SkLAv7vIKxy2vLWLLriKeumwQ7ZrHeR1JBKD6tngFxWXERkXwl49WMuLBGbw0d6NuTixHFPbF/bfPVjNjdQ73julLelpLr+OIHKRlkxhevPoEXh5/Iqkt4rn7veWMeHAGG3P3eB1NGqiwn1Zxeq9kyivgshN1qylpuMyMk7u3ZtgxrfhqXR5vZHxHx5YJAHy1LpdOLROqPxexUJwUSU9PdxkZGUF/3tpYvjWfvinNPc0gUlflFY6hE78gp7CE4T3acNmJnTmtZxuiIsP+zXKjY2YLnHPpgWwblq/+s7O/5ZzHZ/PZimyvo4jUSWSE8c6Nw7j5tGNYsbWAa1/MYMjEL3h30Ravo4mHwmqopLS8gvs+WM7L8zYzum87RvRs43UkkTrrkBTPrSN7cvNPujN91XbezdxCcrNYAFZsLeDfC7MY1bcdx3duUX3CU8Jb2BT3zj37+NWUhczdkMf1w7vxu1E99UssYSU6MoKRfdsxsu9/b/ixfGs+L83dxLOzv6VZXBQnd2/Nqd3b8NPjOuhuTmEsbIo787tdLNi0k79dPIDzB6V6HUekXlyU3pGz+rdn5uocZqzezqy1OUxbtb3638C7i7aQX1TKsanN6d2+mco8TPj65OS3uXtYkrWLsQM7ALAxdw9prZuE/O8Vaaicc2zZVURqi8oZKFfudy/M6Eije3Iip/RozZ1nVS6wtn13Ma2axOrdaQNQm5OTAR1xm9lo4DEgEnjGOXd/HfLVSXmFY8bq7bw0bxMz1+TQIiGGM/u0JSEmSqUtjZ6ZVZc2wPO/HMz3+cUsydrF4qx8VmwtIH/vf+/49LMn55Czu4SOLePp2DKBji0SGHZMK0b3aw/Aprw9JCXE0CwuSgtiNSA1FreZRQJPAmcCWcA3Zva+c25FKIM55ygoLiNndwnJzWJpFhfNp8u3cefbS9mxZx/JibH8+vTuXHZiJxJiwmbERySozIyUpHhSkuKry3h/t5zRnfXbC9mUt5fvdu5l4aadlDvH6H7tKa9wnP7wTMorHDFREbRMiCEpIZpLBnfkqmFdKCkr5/5PVpEYG0VCbBRNYiKJj4ni2NTm9GibSElZOUuz8omJiiAmKoKoiAiiI42WTWJIjKu81H9PSTkREZWzZyLMiIwwIs2I0DuAIwqk8U4A1jnnNgCY2WvAWCDoxb0kaxcXPz2X8gpHWYXjh1GcJy8dxDnHtie1RTyn9UzmJ72TObNPW6I1l1WkTi5O73jQ18rKK4DKd7cPXXQsubv3kbunhLzCfezaW0pC1XLIhcVlvJmRRWFJ2QHf/7tRPenRNpHtBSVc+I+5Bz3/fWP6cuVJaazdXshZj3150OMPXTSAC49PJWPjDi5+ei5mhgFmYBh/v/Q4RvVtx5drc7j2xf8OyVZuBc9emc5Jx7TmP8u2cesbmQc8twFTrh3CwI5JvL0wi7vfXXbQ3//eTcM4JjmRl+Zt4v6PVx70+NRbh9MhKZ6nZ67n8S/WHvDYBzefTNc2TQ/6nmCrcYzbzC4ERjvnrqn6/ArgROfcTT/abgIwoerTnsDR3q+pNZB7lN/rV9rn8NfY9he0z7XV2TkX0BzmoI0xOOcmA5Pr+jxmlhHoAH240D6Hv8a2v6B9DqVAxhq2APu/n0qt+pqIiHggkOL+BuhuZl3MLAa4BHg/tLFERORwahwqcc6VmdlNwKdUTgd8zjm3PISZ6jzc4kPa5/DX2PYXtM8hE5ILcEREJHQ0n05ExGdU3CIiPuNZcZvZaDNbbWbrzOz3h3g81sxer3r8azNLq/+UwRPA/t5qZivMbImZfWFmvr9lT037vN92F5iZMzPfTx0LZJ/N7OKq13q5mb1S3xmDLYDf7U5mNt3MFlX9fp/tRc5gMbPnzGy7mR189U7l42Zmj1f9PJaY2aCgh3DO1fsfKk9yrge6AjHAYqDPj7b5FfCPqo8vAV73Ims97u9pQELVxzf4eX8D3eeq7RKBWcA8IN3r3PXwOncHFgEtqj5P9jp3PezzZOCGqo/7ABu9zl3HfT4VGAQsO8zjZwOfUHmh5hDg62Bn8OqIu/oyeufcPuCHy+j3NxZ4oerjt4CfmH9Xualxf51z051ze6s+nUflfHk/C+Q1Bvgz8FeguD7DhUgg+3wt8KRzbieAc257PWemA9DIAAACL0lEQVQMtkD22QHNqj5uDmytx3xB55ybBew4wiZjgRddpXlAkpkdvFBMHXhV3B2A7/b7PKvqa4fcxjlXBuQDreolXfAFsr/7G0/l/7H9rMZ9rnoL2dE591F9BguhQF7nHkAPM/vKzOZVrbzpZ4Hs873A5WaWBXwM3Fw/0TxT23/vtaZl9RoYM7scSAeGe50llMwsAvgbcJXHUepbFJXDJSOofFc1y8z6O+d2eZoqtH4BPO+ce9jMhgIvmVk/51yF18H8yqsj7kAuo6/exsyiqHyLlVcv6YIvoGUDzOwM4C5gjHOupJ6yhUpN+5wI9ANmmNlGKscC3/f5CcpAXucs4H3nXKlz7ltgDZVF7leB7PN44A0A59xcII7KxZjCVciXCfGquAO5jP594Mqqjy8EprmqkX8fqnF/zew44GkqS9vv455Qwz475/Kdc62dc2nOuTQqx/XHOOdCf+uk0Ank9/pdKo+2MbPWVA6dbKjPkEEWyD5vBn4CYGa9qSzunHpNWb/eB8ZVzS4ZAuQ7574P6t/g4ZnZs6k82lgP3FX1tT9R+Y8XKl/cN4F1wHygq9dnk0O8v58D2UBm1Z/3vc4c6n3+0bYz8PmskgBfZ6NyiGgFsBS4xOvM9bDPfYCvqJxxkgmM9DpzHff3VeB7oJTKd1DjgeuB6/d7jZ+s+nksDcXvtS55FxHxGV05KSLiMypuERGfUXGLiPiMiltExGdU3CIiPqPiFhHxGRW3iIjP/H/SJw49L2k8mwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "y5 = [beta_pdf(i, 4, 8) for i in x]\n", "plt.plot(x, y5, '--', label = 'Beta(4, 8)')\n", "plt.legend(loc=0,numpoints=1,fontsize=10)\n", "plt.ylim(0, 6)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "- If you started with a Beta(20, 20) \n", " - expressing the belief that the coin was roughly fair\n", " - your posterior distribution would be a Beta(23, 27), centered around 0.46, \n", " - indicating a revised belief that maybe the coin is slightly biased toward tails.\n", " " ] }, { "cell_type": "code", "execution_count": 60, "metadata": { "ExecuteTime": { "end_time": "2018-10-27T04:57:21.998295Z", "start_time": "2018-10-27T04:57:21.843852Z" }, "slideshow": { "slide_type": "skip" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd4VFX6wPHvyaQCSaghhAABpLdAooIodlBRsKCyFiwo6q6KrvuzrA3b2ta6LO66iqJgVyzYEGkiKCb03glJSCGQ3mfO7487CS0hk2Rm7tyZ9/M8eabdufPemeSdk3PPeY/SWiOEEMI6gswOQAghRONI4hZCCIuRxC2EEBYjiVsIISxGErcQQliMJG4hhLAYlxK3Uqq1UuozpdQWpdRmpdQITwcmhBCibsEubvca8IPWeoJSKhRo4cGYhBBCnIBqaAKOUioaWAP00DJbRwghTOdKi7s7kAu8o5QaAqQCU7XWJUdupJSaAkwBaNmyZVLfvn3dHasQQvit1NTUA1rrDq5s60qLOxn4DRiptf5dKfUaUKi1frS+5yQnJ+uUlJTGxCyEEAFNKZWqtU52ZVtXTk6mA+la69+dtz8DhjU1OCGEEM3TYOLWWmcB+5RSfZx3nQts8mhUQggh6uXqqJK7gDnOESW7gJs8F5IQQogTcSlxa63XAC71vQghfFdVVRXp6emUl5ebHUrACg8PJz4+npCQkCbvw9UWtxDCD6SnpxMZGUlCQgJKKbPDCThaa/Ly8khPT6d79+5N3o9MeRcigJSXl9OuXTtJ2iZRStGuXbtm/8cjiVuIACNJ21zueP8lcQshhMVI4hZCCIuRxC2sqygL3rkQirLNjkQ0gs1mIzExkSFDhjBs2DCWL19+wu3z8/OZMWOGS/suKyvjzDPPxG63195XWFhIfHw8d955Z4PPX7NmDcOHDycxMZHk5GRWrlwJwLx583jsscdcisEbJHELazkyWS9+HtJ+gyXPmx2VaISIiAjWrFnD2rVrefbZZ3nooYdOuH1jEvfMmTO5/PLLsdlstfc9+uijjBo1yqXn33///Tz++OOsWbOGJ598kvvvvx+AsWPH8s0331BaWurSfjxNhgMKa1nyAuxdDi/1PnxfytvGDwru2wqRHU0Lz2qu/u+K4+67eHAnrh+RQFmlnRvfWXnc4xOS4rkyuQsHSyq5Y3bqUY99fFvjSvUXFhbSpk2b2tsvvvgin3zyCRUVFVx22WU88cQTPPjgg+zcuZPExETOP/98Hn/8ccaPH8+hQ4eoqqri6aefZvz48QDMmTOHDz74oHZ/qampZGdnc8EFF+BK/SSlFIWFhQAUFBQQFxdXe/9ZZ53FvHnzuOqqqxp1jJ4giVtYw9MxUF1Rz4MKWneF/DSj9X3xy14NTTROWVkZiYmJlJeXs3//fhYuXAjA/Pnz2b59OytXrkRrzbhx41i6dCnPPfccGzZsYM2aNQBUV1czd+5coqKiOHDgAMOHD2fcuHFUVVWxa9cuEhISAHA4HNx3333Mnj2bBQsWuBTbq6++ypgxY/jb3/6Gw+E4qhsnOTmZX375RRK3EC6bug4+ugYynC08ZQNd04+pIX+vcbWm9R0cBo/kmBKqlZyohRwRajvh421bhja6hQ2Hu0oAVqxYwaRJk9iwYQPz589n/vz5DB06FIDi4mK2b99O165dj3q+1pq///3vLF26lKCgIDIyMsjOzsbhcNC6deva7WbMmMFFF11EfHy8y7G98cYbvPLKK1xxxRV88sknTJ48uTbpx8TEkJmZ2ejj9QRJ3MIagsMga8Ph69UV0KEvXPAcfPUXKMxwPhYB/S6G0c+YF6tw2YgRIzhw4AC5ublorXnooYe47bbbjtpmz549R92eM2cOubm5pKamEhISQkJCAuXl5URHRx81sWXFihX88ssvzJgxg+LiYiorK2nVqhXPPfdcvfHMmjWL1157DYArr7ySW265pfax8vJyIiIi3HDUzScnJ4XvK8qCGSPAXgH9L4VbfobkydDuJOh5NvQaAzgnNVSXQ1iU9HNbxJYtW7Db7bRr144xY8Ywc+ZMiouLAcjIyCAnJ4fIyEiKiopqn1NQUEBMTAwhISEsWrSIvXuN/7batGmD3W6vTd5z5swhLS2NPXv28M9//pNJkybVJu1JkybVjhg5UlxcHEuWLAFg4cKF9OrVq/axbdu2MXDgQM+8EY0kLW7h+354EIr2Gy3sq2YZ9x3Zj12SA8MmwaavILQlFMvwQF9W08cNRrfHrFmzsNlsjB49ms2bNzNihNH90qpVK2bPnk3Pnj0ZOXIkAwcO5MILL+SBBx7gkksuYdCgQSQnJ3PkalujR49m2bJlnHfeeSeMYd26dbUnHo/0v//9j6lTp1JdXU14eDhvvvlm7WOLFi3i2Wefdcdb0GwNroDTFLICjnCL+k5I1td//cvL8PMTMGUxxA31dHSWtHnzZvr162d2GB6zatUqXnnlFd5///16tyksLGTy5Ml8+umnLu83Ozuba665hp9//tkdYdb5Obh7BRwhzDF1HZx0/uHbwREw6EqYur7u7U+eDGHRsPAfMjEnQA0bNoyzzz77qAk4x4qKimpU0gZIS0vjpZdeam54biOJW/iuyFgozjKu28KMPu4T9V+HR8OpU2DHfNi7QibmBKibb775qAk47nDyySfXdu/4AunjFr5LaziwA6I6wzUfQ8o7J+6/PqprRcvQQOG3pMUtfFd6ClSXwTmPQOwg44TkxDn1bz91HQy80hjjDQ13rQhhUZK4he/a8JnRRdJ3rGvbR8ZCWCRoh3FbhgYKPyWJW/gmhx02zoVe5xt9164qyYFh10NQMHQcIEMDhV+SxC1806avjaTb68TjcY8zcQ6M+xf0PAcqiuDq2Z6JL5C4uXyut8q6rlmzhhEjRjBgwAAGDx7Mxx9/XLvd5MmTGTJkCIMHD2bChAm1k37q89NPP5GUlMSgQYNISkqqra9SVFREYmJi7U/79u255557AJg+fTozZ850Ke5G01q7/ScpKUkL0Sz/Olnrx6O0/uqupj0/5R3j+fvXuzUsq9u0aVPjn/TNvVpPa21cukHLli1rr//www961KhRJ9x+9+7desCAAS7te/r06frVV1/VWmu9detWvW3bNq211hkZGTo2NlYfOnRIa611QUFB7XPuvfde/eyzz55wv6tWrdIZGRlaa63Xr1+v4+Li6txu2LBhesmSJVprrUtKSnRiYmKd29X1OQAp2sUcK6NKhG85dtLNqlnGT2NHhvS+ELgHtnwLsb4xTdnnfP8gZJ3gxG3ar8bInho1o3SUgq4j635O7CC4sP5aIMfyZFnX3r0Pl/6Ni4sjJiaG3NxcWrduTVRUFGA0XMvKyhpcB7Km8BXAgAEDKCsro6KigrCwsNr7t23bRk5ODmeccQYALVq0ICEhgZUrV3LKKae4/J64QhK38C1T18Hc22HXIuN2U4tGRXaE+JONfvLdi2HCu3KSsrHiToZDu6Eszzjhq4KgRTto071Zu/VWWdcjrVy5ksrKSnr27Fl730033cR3331H//79GzW55vPPP2fYsGFHJW2Ajz76iKuvvvqoL4GaUrCSuIV/i4yF8gLjuiuTbk6k71hY8DigpE53XVxpGX9zL6x6F4LDwV4J/cY1+330VlnXGvv37+f6669n1qxZBAUdPq33zjvvYLfbueuuu/j444+56aabGox948aNPPDAA8yfP/+4xz766KPjptrHxMSwZcuWht+URpLELXxP/l6IaAs3fN3wpJv6yGQc9yjJgaSbIPmmpn8WJ+DJsq5gdMWMHTuWZ555huHDhx/3+jabjYkTJ/LCCy80mLjT09O57LLLeO+9945quQOsXbuW6upqkpKSjrrfU6VgZVSJ8C1V5VBRDEP+5Nqkm/rUTMapKfcqk3GaZuIc4zNozmdxAp4s61pZWclll13GpEmTmDBhQu3ztdbs2LGj9vrXX39dW2Fw7ty5da6BmZ+fz9ixY3nuuecYOfL4/v0PP/yQP/3pT8fd76lSsC61uJVSe4AiwA5UaxcrWAnRaOl/GN0j3c9o3n5qJuPgPLnWnC4X4VbeKuv6ySefsHTpUvLy8nj33XcBePfddxk8eDA33HADhYWFaK0ZMmQIb7zxBgA7d+6sPXF5pOnTp7Njxw6efPJJnnzyScDok4+JiQHgk08+4bvvvjvueb/++ivTpk1z23tXy5WhJ8AeoL2rQ1VkOKBosoXPGEPPyvKbv68Pr9F61nhjWOBH1xu3A1yThgNaSGpqqr7uuuua/Pxrr71W5+TkuCWWVatW1RuLDAcU/mX3L9BpSONmS9Zn4hwozISX+0GXU+C0O5u/T+HTjizr2pQKgbNnu2/C1oEDB3jqqafctr8judrHrYH5SqlUpdQUj0QiRGWp0VWS0MxukiNFxRnD1/aeeHZeINEeWDzFl3iirGtTnH/++XUOTXTH++9q4j5daz0MuBD4i1Jq1LEbKKWmKKVSlFIpubm5zQ5MBKB9v4OjCrof9+vVPN1GQtoKcDjcu18LCg8PJy8vz++Tt6/SWpOXl0d4eHiz9uNSV4nWOsN5maOUmgucAiw9Zps3gTfBWLqsWVGJwLTte+OybQ/37rfbCFgzGw5shRj/XbbLFfHx8aSnpyONK/OEh4cTHx/frH00mLiVUi2BIK11kfP6aODJZr2qEHXZMNe4XPFv906W6Xaacbn314BP3CEhIXTv3ryZj8J8rrS4OwJzndM4g4EPtNY/eDQqEViOrU/i7skybbpDZCejn/vkW5q/PyFM1mDi1lrvAoZ4IRYRqKaugy+mwO4lxu2m1iepj1JGq3vvcqNoUgMFhYTwdTJzUpgvMhYqS4zrza1PUp9up0HRfvjfObL6u7A8SdzCNxSkQ2gruPVnozaGu1eu6eacppy5WlZ/F5YnE3CEb7CFQO8xh2tiuJMUnBJ+RlrcwnxF2VCwDzonNbxtU9Su/u78dZeCU8LiJHEL82WuMi7jhnlm/7Wrv9cUnJLV34W1SeIW5stYBcoGnQZ77jVKcqDXaON674tk9XdhadLHLcyXkQox/SG0pedeY+IcKDkAL/aErsNh5N2eey0hPExa3MJcWhtdJZ2HNrxtc7VsD9FdD3fNCGFRkriFuQ7thrJDnjsxeay4RGNIoBAWJolbmCvD2fr1VuLuPAwO7YHSg955PSE8QBK3MFfGKmN4XgcvFX+Kc3bJSKtbWJgkbmGutBVgC4bSPO+8XidjrUPp5xZWJolbmMdeDfvXQkWR96ahR7SGdidBhrS4hXXJcEBhDk+Xcj2RuKGw51fPvoYQHiQtbmGOqeug88mHb3tzGnrcMCjKhKIsz7+WEB4giVuYIzIWqkqN68EeKuVan87OqfXvjZcSr8KSJHEL8xRnGZNibvFQKdf6xA4yLnO3SIlXYUnSxy3MoTVoB/Qd55lSrvUxs29dCDeRFrcwR2GGMWOypvXrLbUlXm3GbSnxKixIErcwR9YG4zLWgxUB61Jb4tVh3K6WEq/CeiRxC3NkOVu4Hft7/7VLcqDPhcb13qOlxKuwHOnjFubIWgdtexitX2+bOMfopnk+AbqOgNPv9X4MQjSDtLiFObI3QMeB5r1+RBuI7nK4y0YIC5HELbyvoggO7vJ+//axYgcd7rIRwkIkcQvvy95kXHp7RMmxYgdB3naoLDU3DiEaSRK38L6sdcZlrIldJWAkbu2AnM3mxiFEI0niFt6Xtd7oY47qbG4cNS3+mi8SISxCErfwvszV4LBDsckzFVt3M8ZwSz+3sBiXE7dSyqaUWq2UmufJgISfs1dD9kaoKDS/TohSxsiWbBlZIqylMeO4pwKbgSgPxSL8nS/WCYkdBKtng8MBQfIPqLAGl35TlVLxwFjgLc+GI/za1HXQZfjh275QJyR2EFSVGKvNC2ERrjYxXgXuBxz1baCUmqKUSlFKpeTm5rolOOFnImOhusy47u0a3PWpOUH54USpzS0so8HErZS6GMjRWqeeaDut9Zta62StdXKHDh3cFqDwM0VZEBbt/Rrc9enQF1BwYJv5fe5CuMiVPu6RwDil1EVAOBCllJqttb7Os6EJvxTSwqgP4s0a3PXxxT53IVzQYItba/2Q1jpea50ATAQWStIWTVJZAof2QIwJFQHrIrW5hUXJaXThPTlbAG1OKde6SG1uYVGNStxa68Va64s9FYzwcznOGiW+0uIGozZ3r9HG9T4Xmd/nLoQLpB638J6cTUZ3RJvuZkdy2MQ5ULgfXu4LPc6CU6eYHZEQDZKuEuE92Rshpq/vTXSJjDVqp8gMSmERPvYXJPxazmaIGWB2FMernfq+0exIhHCJJG7hHSUHjP5kXzkxeayOA4wvFke9c8yE8BmSuIV31LRmY/qZG0d9YvobU9/z95gdiRANksQtvKNmsQJf7CqBw+tfSneJsABJ3MI7cjZCi3bQKsbsSOoW45z6XrOsmhA+TBK38LyiLNjwBbQ9yTgR6ItCW0Lb7jKyRFiCJG7heYufh8piqCwyO5IT6zhAukqEJcgEHOE5xxZxytkE06J9t4hTx4GweZ6x6ntoC7OjEaJe0uIWnlNTxMkWaty2hfl2EaeY/oCGXFn1Xfg2SdzCc2qKONmrjNv2St8u4tTROeLli1tlUQXh0yRxC88qyYG2PaBVR0i+2beLOLXpbpR4zdspiyoIn6a01m7faXJysk5JSXH7foVFvTESouLg2k/NjqR+x/bH1/DV/njhd5RSqVrrZFe2lRa38Cx7lbEsmK/OmKwhiyoIC5HELTzr4C6jb9uXanDXRRZVEBYiiVt4Vu3iCT7e4gajP77Phcb13mN8uz9eBDQZxy08K2cLqCBo39vsSBo2cQ6UHoQXukO3kTDybrMjEqJO0uIWnpWzyRhVEhJhdiSuadEWIjsd/k9BCB8kiVt4Vs5ma3STHCmmv9QsET5NErfwnKpyOLgTOlgscXccALlbwV5tdiRC1EkSt/CcA9uMURpWa3F3HGCMhMnbYXYkQtRJErfwnNrFE3x8KOCxaqa+50ilQOGbJHELz8nZBEEh0K6n2ZE0TvvexkQcKfEqfJQkbuE5OZuNJGgLMTuSxgkOM+KW1XCEj5LELTyjKAt2LTYKN1lRx/7S4hY+SxK38IyFT4O9Aor2mx1J08T0h4I0KC8wOxIhjtPgzEmlVDiwFAhzbv+Z1vpxTwcmLOrYKnuZqb696k19alZ9z9kMXYebG4sQx3ClxV0BnKO1HgIkAhcopeQ3WdStdtUbZ792cLg1q+x1dI6E+fIOWVRB+JwGE7c2FDtvhjh/3F/EW/gHq616U5/oLsaImIO7ZFEF4XNcKjKllLIBqcBJwL+11r/Xsc0UYApA165d3RmjsJqSHGjVCVq0hq6nWa/K3rHdPSlvGz9W6+4RfqtRK+AopVoDc4G7tNb1FnOQFXB8k92hyS2qIK+kgrziSuwOTdd2LejZoZV7X0hreLEn9LkIxk937769oSgLfnwENn4B2m4sqtDvYhj9jPX+cxCW0ZgVcBpV1lVrna+UWgRcAEgVHguoqLYTFmzD4dCMfG4hWYXlRz0+bkgcr/9pKFprbn0vld4dW3H6Se0Z3qMdQUGqaS9anAOleYdnIFqNLKogfJwro0o6AFXOpB0BnA9Ip58P01qzbMcB3v11D7nFFXz1l5EEBSmmjOpBaHAQ7VuF0a5VKLYgRasw41egsLyaA8UVLN6aw4zFO4lvE8HVyV24MrkLsdHhjQugZqq41aa6H6kkB/peAlu+hpPOtV53j/BrrrS4OwGznP3cQcAnWut5ng1LNNXqtEM89tVG1mcU0L5VKNec2o1qhybEprj59Ponw0RHhPDlX0ZSWlnNgs05fPxHGi/9tI3uHVpy8eC4xgVRM+PQqi1uMBZVqCyBf3wD8afAWQ+YHZEQtRpM3FrrdcBQL8QimmnptlxueGclMZFhvDBhMOMT4wgLtjVqHy1Cgxk3JI5xQ+LYd7CUzq2NBRDeXrab4vJqbj+rR8P7zNkELTtAy/ZNPRTfENrSWAQi22JDGYXfk6XL/EB5lZ3wEBsjerbj/8b0YdKIhNoukObo0rZF7fWtWYV8kpLOt+szeXHCEIZ0aV3/E3M2Wbub5EixA2H/OrOjEOIoMuXdwrTWvL9iD+e9vITcogpCbEH8+ayT3JK0j/XChCG8c9PJFJZVc9mMX3nhhy1U2R3Hb+iwG+tMWrmb5EgdB8Gh3VBRZHYkQtSSxG1R1XYHD3y+jke/2kjvjpGE2jz/UZ7dJ4b5fx3FhKR4/rNkJ+vS66jjcWgPVJf5T4u75gtIKgUKHyJdJRZUWe1g6ker+X5DFnefcxL3nNe76UP3GikqPIQXJgxhyqgenBQTCUBOYTkxUc6RJzUV9Tr6SeKOddYsyV4PXU81NxYhnKTFbUGv/byN7zdk8djF/fnr6D5eS9pHqknay3cc4PQXFvHl6gzjgZxNgLLeOpP1ie4C4dGQJdMWhO+QFrcF3X5mTwbERXPRoE5mh8KAztEkdW3DPR+vIauwnNuyN6LadofQFg0/2QqUMioFyqrvwodIi9siHA7N28t2U1pZTWR4iE8kbTDGf79788lcMiSOt79fQdWWH9Fte5gdlnt1HGj0cTvqOBkrhAkkcVvE8z9u4al5m5i31vcWJggLtvHa1Yn8N/5HQhwVFB7wvRibJXYgVJUYo0uE8AHSVWIBs5bv4b9LdnHd8K5cmRxvdjjHezqGoOoKhgEoiM7faM3FE+pTs6hC9gbrLXws/JK0uH3c8p0HeOKbjZzXryNPjBuIUt4/Edmg4xZPCKOg12X8s99n2B1+ULo9ph+g4Me/y6IKwidI4vZhDofm8a820r19S16dmIjNhNEjLjlu8YQqMkuDmf5HEQ/PXU9jSgf7pJAIozpgQbosqiB8gnSV+LCgIMW7N59CeZXdI7Mh3eqYxRP6FWdz59knMX3RDqIjQnjwwr6++d9CQ2RRBeGDfDwbBK5l2w9wWs92tUWefN7Vs+H5BIgfDRe/DMB9WlNYXsV/l+4iKiKEv5x9krkxNsXUdcaiCpu+BEeVsYZmv0uMRRWEMIl0lfigRVtzuO7t3/nwjzSzQ3FdYQaU50PsoNq7lFJMu2QAlybGsXzngbprm/i6mm4gR7Vxu7pCFlUQppMWt485WFLJ/Z+to0/HSK4Y5oMjSOqT5Sx9ekTiBqO758Urh6A1hHihnopHlORA4rWwZjbEJ8uiCsJ0Fv1L8k9aax75cj35pZW8cnUi4SGNq6Vtqpop4XVUBQyxBREaHMTBkkquf/t3NmTUUZzKl02cA5f+G6K7Quuuxm0hTCSJ24d8uSaD79Znce/5vekfF2V2OI2TvR7adDe6FepRZXewM6eYm9/9g8z8Mi8G5yadBsP+tWZHIYQkbl/SKTqCcUPiuG2UBSd5ZK0/rpvkWB2jwnn35lMoq7Rzy6wUSiurvRScm3RKhLwdUF5odiQiwEni9iHDe7Tj9T8N9d3x2vWpKIKDuxtM3AC9O0byr2uGsiWrkHs/XoPDShN0Og0xLqXglDCZJG4fsGBTNs9+v5nyKrvZoTRN9iZAH54a3oCz+sTw8Nj+bMkq4kBJRcNP8BU1iTtzjblxiIAnidtkxRXVPPrVBhZvybVeS7tGdt0jSk7k5pEJfHf3GcREhnsoKA+I7AiRnaSfW5hOErfJ/vnjVrIKy/nH5YOsO1wua72x2EC068MXlVK0DAumotrOY19tYGOmRUaadBoiiVuYzqKZwj9syChg1oo9XHdqN5K6tTE7nKbL2gCxg41FBxqpqLya+Ruzue39VA6VVHogODfrNAQObIXKUrMjEQFMEreJnv52E21bhPK3MX3MDqXpCjIgcxW0adriCe1bhfHf65PIKargzg9XUe3rsys7DQHtOLy2phAmkMRtoqfGD+SfVw0hOiLE7FCa7qfHjER2aFeTdzGkS2ueuXQgv+7I48Uft7oxOA+oOUH5xRQp8SpMI1PeTeBwaIKCFL06RtKrY/0TVnzasVXz9vzSrMUTrkzuwrr0Aub8nsbkM7r77knLqM5gCzO+qJY8X1tQSwhvUp6olZycnKxTUlLcvl9/8cIPW0g7WMqrVycSbNUTkkVZzqp5c40CTMER0O9io2peEwswVVY7yMwvI6F9SzcH6ybHflnVkBKvwg2UUqla62RXtm0wayiluiilFimlNimlNiqlpjY/xMCVllfKW7/sJjQ4yLpJG46pmqfA3vyqeaHBQSS0b4nWms9S0ymu8LGZlTUr/QQ5/1ENjoBBV8LU9ebGJQKOK10l1cB9WutVSqlIIFUp9ZPWepOHY/NLz/+4BVuQ4oEL+podSvOV5EBQKPS9EFq0d1vVvG3Zxdz/2VqWbMvl9YmJvrMAQ+2XlXOiVHW5lHgVpmiwyae13q+1XuW8XgRsBjp7OjB/tDrtEN+u28+tZ3SnY5SP9uE2xvlPgqMSTjrP6Ot1U9W8PrGR3De6D9+szeT93/a6ZZ9uU5IDQyYa17ucKiVehSkadXJSKZUADAV+r+OxKcAUgK5du7ohNP/z+s/bad8qjClnWrCIVF0yVxuXcUPdvus7zuxJ6t5DPDVvE4M6RzO0q4+Mc6/5ctq91JhwNOFtc+MRAcnlTlalVCvgc+AerfVx5dG01m9qrZO11skdOnRwZ4x+459XDuGN64b5/vqRrspcbSzl1cH93T5BQYqXrxpCx6hw7vxgte/VcemcBBlyAl6Yw6UMopQKwUjac7TWX3g2JP9jd2iCFLRrFUa7VmFmh+M+mauN+iQ2z4xDb90ilBnXDiOvuNL3FpWIT4bNX0PJAWjZ3uxoRIBxZVSJAt4GNmutZdBqE3z0RxqXzlhOfqkFpnS7ymE3anZ4oJvkSIPjW3N23xgAsgvLPfpajdLZOWorI9XcOERAcqWrZCRwPXCOUmqN8+ciD8flN8qr7Lz+83aCg5S1Z0geK28HVBZ7PHHXWLApmzNeWMRvu/K88noNiksEZYN06S4R3ufKqJJlWmultR6stU50/nznjeD8wazle8gurOD+MX18Z1ibO3jwxGRdhvdsR3zrCO76cDU5RT7Q8g5tCTH9pZ9bmMLCM0B8X0FZFTMW7+TM3h04tUdJ3rGeAAAW6ElEQVQ7s8Nxr8zVENIC2vf2ysu1CgtmxnXDKCqv4p6P1mD3hZVz4pOMrhKHjxfGEn5HErcHzf5tLwVlVfyflav/1SdztVFwKch7Jw37xkbx5PiBLN+Zx2s/b/fa69arczKUF8DBnWZHIgKMn4xL8023nNGdfp0iGdg52uxQ3Cs/HfathGE3eP2lr0ruwqq9h2gR6gOjTOKdJyjTU6B9L3NjEQFFEreHaK0JC7ZxTl8/nA790yOAhkO7TXn5Zy8f5BvnC9r3hpCWsGAa9DxHpr4Lr5GuEg/ILizn3JeXsGKnj4yAcJenY4zSrRvnGrd3LzFuPx3j1TBqkvbirTlMeS/FvMUXgmwQHgXFWUaJVyG8RBK3B8xYtIO0vFI6t44wOxT3qqmOp5zdFCZXxysoq2L+pmxenG/C4gs1X2JF+43bKW+b8iUmApMkbjfLyC/jw5X7uDI5nq7tWpgdjnvVVMfTdlBBbinl2hzjEztz7ald+e+SXfy0ycvFnmq+xGyhxm1bmJR4FV4jidvNpi/cAcCd5/jpyarCDOMyeTIk3WR6dbxHL+7PwM5R3PfJGvYd9OICvjVfYvYq47bJX2IisEjidqOM/DI+TdnHxFO6+F83SY2TJxuX/ce7tZRrU4WH2JhxTRIa+Cw13bsvXpIDyTcbJylbdTL9S0wEDhlV4kZx0eHMuHYYQ7q0NjsUz0n7zVgBpnOS2ZHU6tquBd/dfQbxbbz8ZVnzpfXDQ5AyEya8493XFwFLWtxupJRi9IBY/1gkoT77fofYwRDqW/33Xdq2QCnFztxi7/d3dzvNWA2npgyAEB4midtNnvhmI9MX+sBsPk+yVxlTvLsONzuSej373Wbu+nAVW7KOKxnvOV1HGJd7l3nvNUVAk8TtBnvzSnhvxV7ySvyobGtd9q8zWpZdTjU7knr94/JBRIWHcPv7qRSWV3nnRVu2NxaT2LvcO68nAp4kbjf418IdBAcp7vCXJcnqs+8349KHE3dMZDj/vnYY6YfKuO+TtTi8VYyq20hI+x3sPrYyvfBLkribac+BEuauzuDaU7sR489922D0b7fuClGdzI7khE5OaMvDY/vx06ZsPk7Z550X7XYaVBZBtozjFp4no0qaafoio7V9+1k9zA7Fswr3w5bvoPcFZkfikhtPSyAixMalQzt75wW7nWZcfnoj3DxfxnMLj5IWdzPdMCKBpy4dSEykn7e2f3oUHFVQao36K0opJp7SlfAQGwVlVWTkl3n2BaPiIDQSDu2RuiXC45TW7u8DTE5O1ikpsjKIX3g6Bqorjr8/OAweyfF+PI2ktebyN5ZTVmnniz+fRotQD/yTafH3SPgGpVSq1jrZlW2lxd1Eew6U8LdP1/rWArae4GOFpRpLKcW95/VmW3YR//fZOjzRUDlct8S5pqjULREeJom7iaYv2sE3azPxhbLQHhUZC2GtfKawVFOM6t2BBy7oy7fr9vPGEg+sVlNbt8Q5osReabn3SFiLnJxsgpqRJDeeluD/fdtg9NsCnH4vlOVbsibHlFE92JBZyIs/bqVvbKT7F7ioqVuyeylUFFryPRLWIYm7CWrGbd92pp+PJKnRdyzsWgSJ10I7a45VV0rx4oTBtAqzMTDOA0vJ1dQt+flJWPYqXDrD/a8hhJN0lTTSrtxi5q5O5/rh3QKjtQ2wa7Exfruttb+owkNsPHv5YGKiwqm2Ozwzs/Kk84xupV2L3b9vIZwkcTdSq/Bgrh/ejdv8fZZkDXs17P4FepyFv3Toa635yweruOXdFCqq7e7defwpEBYNOxa4d79CHEESdyPFRIbzxPiBdIgMMzsU78hcDRUF0ONssyNxG6UUYwfHsXLPQR78fL17R5rYgqHHmbB9AXhiBIsQSOJulDcW7yR170Gzw/CuXYsBBd3PNDsStxo3JI77zu/N3NUZvP7zDvfu/KTzoCgTcja7d79COEnidtH27CJe+HELCzYH0ISKoixY/jrE9IOW7cyOxu3uPOckrhgWzysLtvHl6gz37fik84zLDydCkYwuEe7XYOJWSs1USuUopTZ4IyBf9cqCbbQIsXHrGdY+QdcoC58xhrYF2cyOxCOUUjx7+SAuHtyJnh1auW/H0Z0hvA3k75Xp78IjGpzyrpQaBRQD72mtB7qyU3+b8r4ho4CL/7WMu8/txV/P7212OJ4XwFO4D5ZU0rZlaNN3EMDvnWget05511ovBQKsY/doL83fSnRECLec0d3sULyjZgp3kDWnuTfVf5bs5IJXlzZvtfja6e/Ok9e20IB474R3SR93A7TWnNy9LfeN7k1UeIjZ4XhHzTR3h3WnuTfFOX1jKK+yc8PMlRworqPV7Iqa6e8O5xhxmf4uPMBtiVspNUUplaKUSsnNzXXXbk2nlOLPZ53EpBEJZofiXQec62ee/XdIuikgpnD37hjJ2zeeTGZBGTfMXNn0CTolOcZ7NnACEASFbjzxKQRuTNxa6ze11sla6+QOHTq4a7emStlzkK/WZHhv+StfEp8EQcFw8q1w8cuHp3T7uZMT2vKf65LYll3ELbNSmvbZT5xjvGen3Ao4YOAVbo9TBDapVVIPh0PzxDebyCuuYMyAWML9dGRFnbSGTV8bY7cjWpsdjded1SeGV68eSrXDQVBQM2aLxp8CkXGwcS4Mvsp9AYqA58pwwA+BFUAfpVS6Umqy58My37z1+1mfUcB9o/sQHhJASRsgeyMc2g39x5kdiWnGDu7E+ERj2bNVaYeaNjU+KAgGXGZMfy/Ld3OEIpC5MqrkT1rrTlrrEK11vNb6bW8EZqbKagf/dJb/9Nqahb5k89eAgj5jzY7EdJn5ZUx88zfu/GA1ldWOxu9gwGXGCcq3zpPJOMJtZFRJHT74fS9pB0t54MK+2Jrzr7IVFWXB8n9BfDK08o9zFc0R1zqCR5wrxt/14arGJ+/4ZAhtCXnbZTKOcBvp465DXOsIrkyK56zeAZi4fnwYqkoPL1UmmDQiAYdDM+2bTdw+O5UZ1w5zrfvs2Mk4KW8bPzIZRzSTLBYsDDLjr0Ef/J7Gw1+u58nxA7l+eLeGn1CUBT8+YnQ92SsgKAQGXAqjn5Fx3eI4slhwE+07WMq/F+2gvMrNNZqtYOo657hjZ9dQgMyWbIxrTu3KZ7eP4LpTu7r2hKMm4yjjMjRSkrZoNkncR3j2+81MX7iDQ6WVZofifZGxUFkMaKNlGCCzJRsrqVtblFKk5ZVy7Vu/kZlfduIn1EzGOecR43bOJs8HKfyeJG6n33bl8d36LO44qyedoiPMDscc+9cZtTVu/iFgZks2VXZROev2FXDFG8vZnl1U/4Y1k3FG3Akt2kHL9t4LUvgtSdyA3TnZpnPrCKaMCqCyrUcqOQAluZA82RgJEUCzJZvi5IS2fHTbcKodmiveWM7yHQdO/ISQcBh6PWz5Fv53rgwNFM0iiRv46I80Nu8v5O8X9Qu8yTY11swx+mCTbjQ7EssYEBfNF3ecRseocCbNXMmSbQ3U6Em+CdCQkSJDA0WzyHBAYEh8a248LYGLBsWaHYo5CvfDon9A5ySI6Wt2NJbSpW0LPv/zabz603aSu7Wpf0MZGijcSFrcwMDO0UwbNwDlJ6uYN9q390J1uTFRRDRaVHgIj13Sn5ZhwZRUVPPQF+uPLwtbW6fbuUiDLURG7YgmC+jEvXRbLv/36VqKmlq+0+qejoFp0bD1e+P27qXG7adjzI3Lwtam5/P5qnQu+dcyVqcdOvxA7dDAakCBvcr4opRRO6IJAjZxl1XaefjL9aTuPURocIC+DVPXQc9zD9+WsdvNdlrP9nxxx2nYghRX/XcFb/2y63Bp2JqhgWNfNm7v+8O8QIWlBWjGMhb/3XewjGcuG0RYcICekIyMhdwtxnVbmIzddpOBnaOZd9fpnNUnhqe/3cwLP241HqgZGph8EyScAUX74e0xMsJENFpAJu7fduXxv1928adTujKiZzuzwzFHURb853RjdZbOyXDrzzJ2241atwjlzeuTePrSgVzrnGlZWxpWKTjnUSg7CPt+lxEmotECrlaJ1pqxry+jtLKab+8+g5ZhATqwZt5fnaMawuH+XXJi0sO01tz6Xgotw4J5dfsYlF3qwoijNaZWScBlLaUUb92QTEFZVWAm7WOHpVWXwz/iJGl4mEMbXSjTF+5gS8S/eavzXOKzFhgJXNlg4OVG8SkhXBBQXSW7cotxODRxrSPo1ynK7HDMMXWdcw1EKSblTbYgxT3n9earO0cS2roTS/aWo+2VaGUDbTdWyJFzC8JFAZO4t2cXccm/lvH8D1vMDsVckbGQtxMpJmWOAXHRfPmXkZwRp/k6+AKKr/sBItqidy+FnC3wzoVyslI0KCD6CgrKqrj1vRRahAVz08juZodjnqIsmHMVZK2Fdj3hylmQ8o6ckPQyW5Ci25/nEu/Q2IIU1dd9Cf87m5J3JhBVnoFa8rwx+kSIevh94rY7NFM/Wk1Gfhkf3jqc2Ohws0Myz+JnjaQd0gJu+Rki2kiCMFHNsni2meejsBNdts94QKbDiwb4feJ++aetLN6ayzOXDSQ5oa3Z4Zjj2BOSVaXwfIIkBh+h7lmH/uHvsHEuCgcVOpjlQUmMaFdFeFG2dGOJ4/h94j6nbwx2B1x7qgtLTfmjoizoOAgqSyB3s3FfcAT0u1hGMfiKyFhUeDQojdaKUFXNQFsaYTnZsOR5fu37d7q2bUGXti3MjlT4CL89ObkxswAwVix58MIArni35AXISD0iacsMSZ9UkgNJN6NswSigQ/V+lHZAytuMnN2TDq915cZ3VvLTpmyq7Y1caV74Hb9scb+9bDdPzdvE/yYlc37/AE1O9S3+qx0yQ9IX1Sxaceb98M1U2PaDcdsWSmnCueRn7yM7Yy+3vpdL+1ahPDK2P5cO7WxevMJUfpW4q+wOnvhmI7N/S+OCAbGc1aeD2SGZo6Z7pDwf8nYY9x3ZPSItbd8VGQuRcdSOs7dX0mL/H7QoO8i8Ycv5pdNNxC34CwUhbwKd2ZRZyOer0hkzIJakbm1qT3gK/+Y3iftQSSV/nrOKFbvyuP3MnvzfmD6B9UtclAWf3QQT3oUFTxirrNSQ7hFrKcmB5Jth9ftgr4RSY1k0W+pMzmKmsc2e/8HAl9mYWcD7K/by9rLdRIUHc3qv9ozq1YFLh3YO3NWcAoDf1CpZtCWH295P5bkrBnH5sHivvrZpjkzWS543hpDVxRYCQycZ3SOyjqR1FGXBjw/D5q+M+t11UpTctZEVOw7Qfcld3Fl1F7vLW7Hu8TGEBgfx5eoMCsqqGBwfTb9OUZLMfVjA1CrZfaCEden5jE/szNl9Y5h/7ygS2vthsaQjEzT6iGT9AuxdDi/1rvt50j1ibZGxxn9JDjsEhYKj8vBjKghaxkBxNi1XvMR5AGXr+S5pOfsT7yb0/bEw4V0Wpqzj2n3TmFJ5N/m2NvSKieSM3u156MJ+AOQUldOuZVhg/XfqB1waVaKUukAptVUptUMp9aDHoinKOjzlt57r9oL9HPr3eUx96weufmku8V9OoPRgBhRlkfDNhBM+16vXXTwel64veQHSfjNa1Ucm6zpb2M6PVOpr+4eaxRemLIT2R4yO0g4ozgL04Qk72oFKmUncW4nG78iS53mt03xOCdrK5/2XMvXUKF4ufQh7QVbt79fk6d8y6tEPWfvUaUx96wee/3QxB6efW/s7WP7mGApz09GF+635t2Pm378H2aZNm3bCDZRSNuAHYAzwLPD6E088sXTatGn1Lmn95ptvTpsyZUrjo/npMdgyDypL0HuXweZ5FBbmo9KWE7z1W9Kycvnp+y8YVrKM0uJCbuycQWLxMkLtZbD319rn+sT13mOOOp66tyuGPctgy7dQdgh2LTaWESvKNEYV7FgAK/4FmasBbVxmrq7rU3J+mmGgq6FDX7juc+O+4mxnUSlhSQOvMH6XWnU0fj8SzoBzH4d9K6G8ADhBV2fmalTmahQQfWgDp+z/gPb2HEbFVhn72vYjw+OCODd8G4NLfoXyIlpk/U5S6a+oqlIce5Zh2zqPT1dsZcOK7xlYuJR5qTvQe5bRPv0nqiuKSV0yj077F7BxbxaHNi+kXdp8CgrzCc/4Db1lHrl5BynbvpSInd9RWFhA5c5fCN3+HaqqBMce42+8qrwIx55lBG39FkdFMfQag1rwuPl/w839+2+kJ554Yv+0adPedGXbBvu4lVIjgGla6zHO2w8BaK2fre85je7jrm/omjiBIMBxuGXdoS9c8dbh2iPSl+3fvrkXVr1rLD5cXW7cp4KNL24UJ0zowjsaOTO5MX3criTuCcAFWutbnLevB07VWt95zHZTgJpmdh9gq6sBh9oI6RKl4qPDVGulCNKglbMZacnrGke1pjpYEeyu46nrvgo7ZbsPOXZ3aKk6hASpkO0HHTtdfc99QHvggNlBeJFbj7dX26CeVQ5dlVuic3u0CeoJsOuQY2f3NkHdw2xEnOh3yCf+Rrz4t2PGMRRU6Py0Ar2vykF1Iz7Wblprl8Ywu+3kpNb6TcClZv6JKKVSXP3W8RdyzP4v0I4X5Jg9yZWTkxlAlyNuxzvvE0IIYQJXEvcfQC+lVHelVCgwEfjas2EJIYSoT4NdJVrraqXUncCPgA2YqbXe6MGYmt3dYkFyzP4v0I4X5Jg9xiMzJ4UQQniO35Z1FUIIfyWJWwghLMa0xN3QNHqlVJhS6mPn478rpRK8H6X7uHC8f1VKbVJKrVNK/ayUsvySPa6WSlBKXaGU0kopyw8dc+WYlVJXOT/rjUqpD7wdo7u58LvdVSm1SCm12vn7fZEZcbqLUmqmUipHKbWhnseVUup15/uxTik1zO1BaK29/oNxknMn0AMIBdYC/Y/Z5s/Af5zXJwIfmxGrF4/3bKCF8/odVj5eV4/ZuV0ksBT4DUg2O24vfM69gNVAG+ftGLPj9sIxvwnc4bzeH9hjdtzNPOZRwDBgQz2PXwR8jzEhZzjwu7tjMKvFfQqwQ2u9S2tdCXwEjD9mm/HALOf1z4BzlVJWLWHW4PFqrRdprUudN3/DGC9vZa58xgBPAc8D5d4MzkNcOeZbgX9rrQ8BaK2tvlqzK8esgSjn9Wgg04vxuZ3Weilw8ASbjAfe04bfgNZKqU7ujMGsxN0Z2HfE7XTnfXVuo7WuBgqAdl6Jzv1cOd4jTcb4xrayBo/Z+S9kF631t94MzINc+Zx7A72VUr8qpX5TSl3gteg8w5VjngZcp5RKB74D7vJOaKZp7N97o1m6Hrc/UkpdByQDZ5odiycppYKAl4EbTQ7F24IxukvOwvivaqlSapDWOt/UqDzrT8C7WuuXnEXr3ldKDdRay6rHTWRWi9uVafS12yilgjH+xcrzSnTu51LZAKXUecDDwDittdXLJTZ0zJHAQGCxUmoPRl/g1xY/QenK55wOfK21rtJa7wa2YSRyq3LlmCcDnwBorVcA4RhFt/yVx8uEmJW4XZlG/zVwg/P6BGChdvb8W1CDx6uUGgr8FyNpW73fExo4Zq11gda6vdY6QWudgNGvP05r7d0179zLld/rLzFa2yil2mN0nezyZpBu5soxpwHnAiil+mEk7nrr+fuBr4FJztElw4ECrfV+t76CiWdmL8JobewEHnbe9yTGHy8YH+6nwA5gJdDD7LPJHj7eBUA2sMb587XZMXv6mI/ZdjEWH1Xi4uesMLqINgHrgYlmx+yFY+4P/Iox4mQNMNrsmJt5vB8C+4EqjP+gJgO3A7cf8Rn/2/l+rPfE77VMeRdCCIuRmZNCCGExkriFEMJiJHELIYTFSOIWQgiLkcQthBAWI4lbCCEsRhK3EEJYzP8D3jZ3r4bCUZkAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "y5 = [beta_pdf(i, 4, 8) for i in x]\n", "y6 = [beta_pdf(i, 23, 27) for i in x]\n", "\n", "plt.plot(x, y5, '--', label = 'Beta(4, 8)')\n", "plt.plot(x, y6, '-*', label = 'Beta(23, 27)')\n", "plt.legend(loc=0,numpoints=1,fontsize=10)\n", "plt.ylim(0, 6)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "- And if you started with a Beta(30, 10) \n", " - expressing a belief that the coin was biased to flip 75% heads\n", " - your posterior distribution would be a Beta(33, 17), centered around 0.66. \n", " - In that case you’d still believe in a heads bias, but less strongly than you did initially. " ] }, { "cell_type": "code", "execution_count": 58, "metadata": { "ExecuteTime": { "end_time": "2018-10-27T04:56:07.907170Z", "start_time": "2018-10-27T04:56:07.738499Z" }, "slideshow": { "slide_type": "skip" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsnXlclOX6/9/PsO8oq4gsLrgggkCm2Z5mZVmWlWX7dsoW61fflrO0n/Y903M6p9LK9o6VtqlppmkaKKKiIqAiIPu+MzP374+HQVSQAWbmmRnu9+vFa2ae5X4uhuF67rnu6/pcihACiUQikTgOOq0NkEgkEknvkI5bIpFIHAzpuCUSicTBkI5bIpFIHAzpuCUSicTBkI5bIpFIHAyzHLeiKIGKonylKMpeRVH2KIoyxdqGSSQSiaRrXM087k3gJyHEHEVR3AFvK9okkUgkkpOg9FSAoyhKAJABDBeyWkcikUg0x5wZdyxQBnygKEoikA4sEEI0dD5IUZQ7gDsAfHx8UsaMGWNpWyUSicRpSU9PLxdChJhzrDkz7lTgD2CqEGKLoihvArVCiH90d05qaqpIS0vrjc0SiaQXNOub2Vu5lzGDx+Dp6qm1ORILoChKuhAi1ZxjzZlxFwAFQogt7a+/Ah7tq3ESiaRvCCH4KOsjPt7zMUcajpyw/67Eu5ifNF8DyyS2pkfHLYQoVhTlsKIoo4UQ+4DzgCzrmyaRSEwsyljE4h2Lu9x338T7uH3C7Ta2SKIl5maV3Assa88oyQNutp5JEonkeG6Kv4mvs78m2DuYT2d+ik7RIYTgkQ2P8Pb2t4kPiue0oadpbabERpjluIUQGYBZsZfuaGtro6CggObm5v4MI+kHnp6eREZG4ubmprUpkl7y353/pbSplFfPfhWdopZfKIrCk1OeZH/Vfh7Z8AifX/w5Eb4RGlsqsQXmzrj7TUFBAX5+fsTExKAoiq0uK2lHCEFFRQUFBQXExsZqbY6kFxTUFbB091JmDp9JUmjSMfuW7F5CTnUOADO+ntGxXca7nRubOe7m5mbptDVEURSCgoIoKyvT2hRJL3kt/TVcdC7cn3z/CfvmJ81nftJ8/rbxb6zLX8e6q9fh4eKhgZUSW2JTrRLptLVFvv+OxaKMRSQsTWD1odU06ZuY/tV0EpYmsChj0QnHzhw+k7q2On4r+E0DSyW2xmYzbolE0jvmJ81HQWHxjsWsvWotwV7B3R57avipBHsFsyJ3BdOjp9vQSokWSHVAicSOWXVoFclhySd12gAuOhcuir2IDYUbqG6utpF1Eq0YUI7bxcWFpKQkEhMTSU5OZtOmTSc9vrq6mkWLTvxa2hVNTU2cddZZGAyGjm21tbVERkZyzz339Hh+RkYGkydPJikpidTUVLZu3QrAypUrefzxx82yYcBRVwwfXAh1JVpbYhXyavLIqc4xewZ9yYhL0Bv1rDq0ysqWSbRmQDluLy8vMjIy2LFjB88//zyPPfbYSY/vjeN+//33ufzyy3FxcenY9o9//IMzzzzTrPMffvhhnnjiCTIyMnj66ad5+OGHAZg5cyYrVqygsbHRrHGcns7O+tcXIf8PWP+i1lZZhTWH1gAwLWqaWcePHjSakYEjWZm30ppmSewAzWLcV/978wnbLp4whOunxNDUauCmD7aesH9OSiRXpg6jsqGVuz5OP2bf53/pnUR4bW0tgwYN6nj98ssv88UXX9DS0sLs2bN56qmnePTRR8nNzSUpKYnp06fzxBNPcOmll1JVVUVbWxvPPvssl156KQDLli3jk08+6RgvPT2dkpISLrjgAszRbVEUhdraWgBqamqIiIjo2H722WezcuVKrrrqql79jk7J+pfg0CZ4Ne7otrT31B8UeHAf+IVpZp4lWX1oNUkhSYT5mPf7KIrCzOEzeXPbmxyuO8wwv2FWtlCiFQNqcbKpqYmkpCSam5s5cuQIa9euBWDVqlXs37+frVu3IoRg1qxZ/Pbbb7zwwgvs2rWLjIwMAPR6PcuXL8ff35/y8nImT57MrFmzaGtrIy8vj5iYGACMRiMPPvggH3/8MWvWrDHLtjfeeIMZM2bw0EMPYTQajwnjpKamsmHDhoHtuJ8NBX1LNzsVCIyC6nx19n3xazY1zRocrj3M3sq9PJT6UK/OmxmrOu4fD/zIHRPusJJ1Eq3RzHGfbIbs5e5y0v2Dfdx7PcOGo6ESgM2bN3PDDTewa9cuVq1axapVq5g4cSIA9fX17N+/n6ioqGPOF0Lw17/+ld9++w2dTkdhYSElJSUYjUYCAwM7jlu0aBEXXXQRkZGRZtu2ePFiXn/9da644gq++OILbr311g6nHxoaSlFRUa9/X6diQSZ8di0Utn/TUlxAmNYTBFQfUp+aZt+uHvD3Uk1MtQSr81cD9CpDpLOeydvb3+bt7W8DshjHGRlQM+7OTJkyhfLycsrKyhBC8Nhjj/GXv/zlmGMOHjx4zOtly5ZRVlZGeno6bm5uxMTE0NzcTEBAwDGl/Js3b2bDhg0sWrSI+vp6Wltb8fX15YUXXujWnqVLl/Lmm28CcOWVV3Lbbbd17GtubsbLy8sCv7UD4+oBxbuOPte3QMgYuOAF+PZuqC1s3+cFYy+G8/+pna0WYPXB1YwPGt+rEnZTMc6LW1/kq+yv2HTNJtxcpLyBMzKgFic7s3fvXgwGA0FBQcyYMYP333+f+vp6AAoLCyktLcXPz4+6urqOc2pqaggNDcXNzY1169Zx6JA6yxs0aBAGg6HDeS9btoz8/HwOHjzIK6+8wg033NDhtG+44YaOjJHOREREsH79egDWrl3LqFGjOvZlZ2czfvx467wRjkBdMSyaAoYWGHcZ3PYLpN4KQSNhxDkwagbQXlykbwYPf4eOcxc3FLOrYhfTos1blDyelLAUmg3NZFVKEU9nZUDNuE0xblDDHkuXLsXFxYXzzz+fPXv2MGWKGn7x9fXl448/ZsSIEUydOpXx48dz4YUX8sgjj3DJJZeQkJBAamoqnbv8nH/++WzcuJFp007+z5aZmdmx8NiZ//znPyxYsAC9Xo+npyfvvvtux75169bx/PPPW+ItcEx+ehTqjqgz7KuWqts6x7EbSiH5Bsj6Ftx9oN5x0wM7hzve2PYGb2x7A+hduGNiqBrySy9JJzEk0TqGSjSlxw44faGrDjh79uxh7NixFr+WvbBt2zZef/11Pvroo26Pqa2t5dZbb+XLL780e9ySkhKuvfZafvnlF0uY6Vh/h+4WJLuLX294DX55Cu74FSImWts6q/HkpidZdWgVG+du7FAC7C2zvpnFML9hvHPeOxa2TmItetMBZ8CGSixNcnIy55xzzjEFOMfj7+/fK6cNkJ+fz6uvvtpf8xyTBZkwstPinKsXJFwJC3Z2ffwpt4JHAKx9zqELc9JL0pkYOrHPThvUcMn2ku0YjN1/HiWOi3TcFuSWW245pgDHEpxyyikd4Z0Bh1841Berz1081Bj3yeLXngFw6h2QswoObXbIwpzypnIO1h4kJSylX+MkhyZT11bXIfkqcS4GVIxb4mAIAeU54D8Urv0c0j44efz6mNCKcMjUwO2l2wHV8faH1DD1G3daSRqjB4/ut10S+0LOuCX2S0Ea6Jvg3L9DeIK6IDl3WffHL8iE8VeqOd7Qc2jFDtlWsg1PF0/ig+L7Nc4Q3yFE+ESQXpLe88ESh0M6bon9susrNUQyZqZ5x/uFg4cfCKP62gFTA9NL0pkQMsEi+dfJYcmkl6RjjQQEibZIxy2xT4wG2L0cRk1XY9fm0lAKydeDzhXC4h0qNbC+tZ59Vfv6Hd82kRKWQmVzJQdrD1pkPIn9YN+O28KynbaSdc3IyGDKlCnEx8czYcIEPv/8847jbr31VhITE5kwYQJz5szpKPrpjtWrV5OSkkJCQgIpKSkd+ip1dXUkJSV1/AQHB3P//Wprq4ULF/L++++bZbfdkvWd6nRH9bIIZe4ymPU2jDgXWurg6o+tY58VyCjLwCiMJIf1L75twnQD2FayzSLjSewIIYTFf1JSUsTxZGVlnbCtR1Y8IMSTgeqjBfDx8el4/tNPP4kzzzzzpMcfOHBAxMfHmzX2woULxRtvvCGEEGLfvn0iOztbCCFEYWGhCA8PF1VVVUIIIWpqajrOeeCBB8Tzzz9/0nG3bdsmCgsLhRBC7Ny5U0RERHR5XHJysli/fr0QQoiGhgaRlJTU5XF9+jtowdunCPGEvxDf3tu389M+UM8/stOiZlmTN9LfEElLk0RDa0O/x3pn+zti/JLxJ/y8s/0dC1gqsQZAmjDTx2qTVfLjo1B8kgWj/N/VjAITpuwARYGoqV2fE54AF3avBXI81pR1jYs7KjkaERFBaGgoZWVlBAYG4u/vD6g3zKamph77QJqErwDi4+NpamqipaUFD4+jDWGzs7MpLS3ljDPOAMDb25uYmBi2bt3KpEmTzH5P7ILji262LVV/epsZEnchcD/s/R7CHUMuYFvJNsYFjcPbzbvfY5l0S+5cfScVzRV8eUnv6gck9o19pgNGnAJVB6CpQl1oUnTgHQSDYvs1rK1kXTuzdetWWltbGTFiRMe2m2++mR9++IFx48b1qrjm66+/Jjk5+RinDfDZZ59x9dVXH3MTMEnBOpzjXpAJy++EvHXq676KRvmFQeQpapz8wK8wZ4ldL1K2GFrYWb6TeWPnWXTc0YNH82HWh7QZ2qTglBOhjeM2Z2a84gHYtgRcPcHQCmNn9Vtn2VayriaOHDnC9ddfz9KlS9Hpji4nfPDBBxgMBu69914+//xzbr755h5t3717N4888girVp3Yluqzzz47odQ+NDSUvXv39vym2Bt+4dBcoz43p+jmZIyZCWueABS71unurE+yZPcSluxeAlhGjnXs4LHojXpya3IZM3hMzydIHAL7nHGDmh2QcjOk3txz4UUfsKasK6ihmJkzZ/LPf/6TyZMnn3B9FxcX5s6dy0svvdSj4y4oKGD27Nl8+OGHx8zcAXbs2IFerycl5dhMBIeWgq0+BF6D4cbv+v63d6BinPlJ8/F39+fFP19k7ZVrCfEOsdjYpuKbvZV7peN2Iuw3q2TuMnWGZE7hRR+wpqxra2srs2fP5oYbbmDOnDkd5wshyMnJ6Xj+3XffdSgMLl++vMsemNXV1cycOZMXXniBqVNPjO9/+umnXHPNNSdsd1gp2LZmaKmHxGv697c3FeOY5F7tvBgnqyKLEK8QizptgCi/KLxcvdhb6YDfviTdYtaMW1GUg0AdYAD0wkwFK3vDVrKuX3zxBb/99hsVFRUsWbIEgCVLljBhwgRuvPFGamtrEUKQmJjI4sXqV+Tc3NyOhcvOLFy4kJycHJ5++mmefvppQI3Jh4aGAvDFF1/www8/nHDe77//zpNPPmmx985mFPyphkdiz+jfOKZiHNoXufsTcrEBWRVZjA2yvGqji86FuEFx0nE7G+akngAHgWBzU1Uslg7oQKSnp4vrrruuz+fPmzdPlJaWWsSWbdu2dWuL3f8d1v5TTQFtqu7/WJ9eK8TSS9W0wM+uV1/bIQ2tDWLC0gli4faFVhn/mc3PiMnLJguj0WiV8SWWgV6kA9pvqMTBMEfW9WR8/PHHhIRY5mtyeXk5zzzzjEXGsjkHNsCQxN5VS3bH3GVwWXsB1bBJFg+3WYrsqmyMwsi4weOsMv7owaOpb6unoL7AKuNLbI+5jlsAqxRFSVcURbaO7gZryLr2henTp3eZmmj3tDaqoZKYfoZJOuMfoaaRHjp5layW7K7YDcC4IOs47rGD1RDMvsp9VhlfYnvMddynCyGSgQuBuxVFOfP4AxRFuUNRlDRFUdLKysosaqRkgHB4CxjbIPaEj1f/iJ4K+ZvBaLTsuBYiqyKLIM8gQr1DrTL+yMCR6BSdjHM7EWY5biFEYftjKbAcOKGqQwjxrhAiVQiRaqmv/JIBRvaP6uPg4ZYdN3oKNFVCuX3OOLMqshgXNK7HKtq+4unqSax/rJxxOxE9Om5FUXwURfEzPQfOB3ZZ2zDJAGTXcvVxs4X7JEafpj4e+t2y41qAJn0TeTV5VguTmBgTNIY9lXuseg2J7TBnxh0GbFQUZQewFfheCPGTdc2SDCieDYUnA9SiK1ALZZ4MULdbgkGx4DfELuPcpoVJa6QCdmbMoDGUNJZQ1Vxl1etIbEOPjlsIkSeESGz/iRdC9FI0wn6wlazroUOHSE5OJikpifj4eP71r391HHfBBReQmJhIfHw8d955Z49ZKHv37mXKlCl4eHjwyiuvdGzft2/fMbKu/v7+vPHGGwA89NBDHTosDsGCTIg96+hrSxfLKIo66z606VjxMjsgqyILoN8db3rCVEG5r0qGS5yBAZUOaNIq2bFjB88//3yXlYqd6Y3jfv/997n88stxcXFhyJAhbN68mYyMDLZs2cILL7xAUVERoBbM7Nixg127dlFWVtZj1/fBgwfz1ltv8dBDDx2zffTo0WRkZJCRkUF6ejre3t7Mnj0bgHvvvZcXXjBfKVFz/MKhtUF93l99ku6IPg3qjsB/zrWr7u9ZFVkM9hxMmLd1C4NM5e57K+QCpTOgiVbJi1tftPgK95jBY3hk0iNmH29NWVd3d/eOcVtaWjB2ymYwVUfq9XpaW1t7XJAKDQ0lNDSU77//vttjfvnlF0aMGEF0dDQA0dHRVFRUUFxcTHh4uNnviabUFIC7L9zyk1W0aYhulwso2m5XglOmiklrLUyaGOQ5iDDvMDnjdhLsUmSqsL6QIw1HTtg+xGcIQ32H9nlcW8q6Hj58mJkzZ5KTk8PLL79MREREx74ZM2awdetWLrzwwmO0TPrKZ599doJeSXJyMr///jtXXHFFv8e3CS5uEDfjqD6JJbFTwalmfTO51bmcFXlWzwf3g87qgyvzVrIybyVgGfVBiTZo4rh7MzO2JLaUdR02bBiZmZkUFRVx2WWXMWfOHMLC1K/DP//8M83NzcybN4+1a9cyffr0Pv9Ora2tfPfddzz//PPHbA8NDe0Iz9g9dSVQcxgm32Wd8Rdkws9/h91fq/rufdX4tjD7q/ZjEAarZ5SYmiq8uPVFvsr+ii3ztqBTBlSU1OkYsH+9rmRdTTHjnJwcbr311hPO6SzrmpGRQVhYWId86vGyriYiIiIYP348GzZsOGa7p6cnl156Kd9++22/fo8ff/yR5OTkjpuCCYeSdS1q74kYYZleiyfQ0f3dJDilbff3RRmLSFiawLU/XAvAA78+QMLSBBZlmLee0ldGBI6g2dDc5bdZiWMxYB23NWVdCwoKaGpqAqCqqoqNGzcyevRo6uvrOXJE/afR6/V8//33HQqDCxcuZOHChb3+PZxC1rVwGyguMGSC9a7RUAqjzlefx12kaff3+Unz2XnjTq4dcy3ert7suGEHO2/cafWwxfAAtbAptzrXqteRWB+7jHFbC1vJuu7Zs4cHH3wQRVEQQvDQQw+RkJBASUkJs2bN6liwPOecc7jzzjsB9UbSld52cXExqamp1NbWotPpeOONN8jKysLf35+GhgZWr17Nv//972POaWtrIycnh9RUB1HfLUyH0HHg7mO9a8xdBg3l8PIIiJoMU++z3rXMJLsqm1GDRtksbGFy3AdqDnBmpIVlBSQ2ZUA57pPlTC9YsIAFCxacsN2UKWJi8+bNXZ5/99138/rrrzNt2jSmT59OZmbmCceEhYXx559/dnn+wYMHee21ExflwsPDKSjoWtXNx8eHioqKE7avXLmSOXPm4OrqAH9eIdRQydhLrH8tn2AIiDoamtEQIQT7q/czPbrv6xu9JdAzkCDPIDnjdgIc4D/bMegs69oXhcCVK1dazBa9Xs+DDz5osfGsStUBaKqCoSk9H2sJIpLUlECNKW0spaalhrhBcTa97vDA4eTV5Nn0mhLLY9MYt7CzqjVLYy+yrldeeWWXzYvt8v0vbJ/92spxD02GqoPQWGmb63VDdlU2AKMCR9n0usMDhpNXnWefnwWJ2djMcXt6elJRUSE/MBohhKCiogJPT0+tTTmWwm1qel6IdbU6OohQUz61nnXvr94PwKhBtnfcdW11lDeV2/S6Estis1BJZGQkBQUFSK1u7fD09CQyMlJrM44lfzO4uEJjhW3S84aoi9MUbYOR51n/et2QXZVNmHcYAR4W6PTTC0YEjgAgtybX4o2JJbbDZo7bzc2N2NhYW11O4ggY9HBkBwiD7crQvQIhaCQUajzjrtpv8/g2HJsSOHnIZJtfX2IZ5OKkRBuOKUPHtmXoERPhoHba3G3GNvJq8jh96Ok2v3awVzB+7n4cqDlg82tLLMeALcCRaMyCTBh6ytHXlpZyPRkRyVBXBHXF1r9WFxysOYjeqNdkxq0oCsMDhsuUQAdHOm6JNviFQ1uj+tzVSlKu3TG0vbT+w0s1kXjtyCix8cKkiRGBI2RKoIMjHbdEO+qL1aKY236BlJttV4YenqA+lu1VY+s2JrsqG1edK7H+2qz5DA8YTmVzJdXN1ZpcX9J/ZIxbog1CqEp9Y2ZZR8q1O7SMrbezv2o/sQGxuLm42eR6x2NaoMyrySPZ00rCXhKrImfcEm2oLVQrJk2zX1uxIBPGX6mKWoFtY+vtZFdlaxLfNtE5JVDimEjHLdGG4l3qY7gVFQG7okPitb0rkd62Eq81LTWUNJbYvGKyM+E+4Xi5epFXLePcjop03BJtKG6f4YZZt4lAlzSUwugL1edx59tU4nV/lVoxqdWMe1HGIhI/TKRJ38THez4mYWmCTbTAJZZFsUYJempqqkhLS7P4uBIn4vProWQX3KdRIUxTFbwYA9OehNMfsPrlOrcP64xW7cMe/u1hMssy+emKn2x+bUnXKIqSLoQwS4tZzrgl2lCyC8I0bPTgNQgChh0N2VgZU/OEq0dfja+bL5k3ZNqkeUJ3xPjHUFRfRKuhVZPrS/qHdNwS29NSB5V5to9vH094wtGQjY3Irc5leOBwq3d174lo/2gEgsN1hzW1Q9I3pOOW2J6SLPXR1hklxxOeABX7obXRJpcTQpBTnaPpwqSJmIAYQK3ilDge0nFLbE9xe3egcI17YoYnqNklpXtscrmK5gqqW6o70vG0JNovGoCDtQe1NUTSJ6Tjltie4p1qjNl/qLZ2mGb8xSe2mbMGJn0Qe3Dcvu6+BHsFc6j2kNamSPqAdNwS21O0HYwGqLdNpWK3BEarOdw2inPnVOcAMDJwpE2u1xPR/tHScTsoZjtuRVFcFEXZriiK5ZojSgYeBj2U7IaWWk10Qo5BUdTMlhLbZJbkVufi5+5HiJd9NDCI8Y+RoRIHpTcz7gWAbYKBEufk2VB4JkhtnACqRsiTAep2rQhPUFMCjUarXyq3OpeRgSM1zygxEeMfQ2VzJTUtNVqbIuklZjluRVEigZnAf61rjsSpWZAJwzp1XdFAJ+QEwhOgrUHtNm9FTBkl9hImATVUApBfm6+xJZLeYu6M+w3gYaDbaYmiKHcoipKmKEqa7Csp6RK/cNA3qc9trcHdHaYFyk/nWlWbu6ypjNrWWrtYmDQRHSAzSxyVHh23oigXA6VCiPSTHSeEeFcIkSqESA0JsY8YnsQOqSsGjwDba3B3R8gYQIHybKvG3O1tYRJgmO8wXBQXuUDpgJijxz0VmKUoykWAJ+CvKMrHQojrrGuaxClx84aoKbbV4O4OG2pz21MqoAk3FzeG+g6VM24HpMcZtxDiMSFEpBAiBpgLrJVOW9InWhug6iCEaqAI2BU21ObOrc4l0COQIM8gi4/dH2RKoGMi87gltqN0LyC0kXLtChtqc5sWJu0lo8SEyXFbQyVUYj165biFEL8KIS62ljESJ6e0XaPEXmbcoGpzjzpffT76IqvE3IUQ5Fbn2lWYxERsQCxN+iZKGzUuhpL0CjnjltiO0iw1HDFImya5XTJ3GVzypvp8+NnqawtT0lhCfVu9XS1MmjClBMo4t2MhHbfEdpTshtAxoLOzj51fuKqdYqUKSlNGiT3OuE2OW8a5HQvZ5V1iO0r3HA1L2BMdpe+7LTrs8V1vbvn5FkC7rjddEeodiperl5xxOxjScUtsQ0O5Gk+2l4XJ4wmLh20fqaXvFvpGMD9pPvOT5vP474+zvmA9669eb5FxLYlO0RHtHy11uR0MO/vOKnFaTLPZ0LHa2tEdoePU0vfqgxYf2qRRYm8sylhEwtIE9lbuZUPhBtk42IGQjltiG0zNCkLjtbWjO0z9Ly0cLhFCkFtjnxklpj6Ytyfcjqviyvbrt2vaB1NiPtJxS2xD6W7wDgJfDZUAT0Zoe+m7qa2ahShuKKahrYERAfbnuE1E+UehF3qK6ou0NkViJtJxS6xPXTHs+h8MHqkuBNoj7j4wONbimSW5NfZX6n48MrPE8ZCOW2J9fn0RWuuhtU5rS05OWLzFQyUmjRJ7jHGbkI7b8ZBZJRLrcbyIU2mW2jjBCiJOFiFsPOxZqXZ9d/e2yJA51TkEeQYR6BlokfGswSCPQfi5+UnH7UDIGbfEephEnFzc1dcuHto3TjgZoeMAAWWWa/RkrxklnVEUhSj/KPLrZEMFR0E6bon1MIk4GdrU14ZW7RsnnIyw9oyX/91ukaYK9qxRcjxR/lFyxu1ASMctsS4NpTB4OPiGQeot2jdOOBmDYlWJ14pcizRVKG4oplHf6BCOO9o/miMNR2g1tGptisQMZIxbYl3mLoPFUyFohPaNE06GFZoq2LNGyfFE+0djFEYK6goYHjhca3MkPSBn3BLrYmhT24LZa8WkCSs0VXCEjBIT0X4ys8SRkI5bYl0q89TYtj1pcHeFFZoq5FTnEOwVTIBHgIWMtB5R/lEAcoHSQZCOW2JdOpon2PmMG9R4/OgL1edxM/odj3eUhUmAAI8AAj0CpUqggyAdt8S6lO4FRQfBcVpb0jNzl8Gl76jPo6f2q6mCSaPEEcIkJqL8o8ivlTNuR0A6bol1Kc1Ss0rcvLS2xDy8B4PfkKPfFPrIkYYjNOmbGB7gOAt9Mf4xMsbtIEjHLbEupXscI0zSmdBx/dYsMWWUONSM2y+KksYSmvRNWpsi6QHpuCXWo60ZKnMhxMEcd1g8lO0Dg77Xp5o0ru/+5W4AbvzpRofRuDZplhyuO6wfk8eqAAAgAElEQVSxJZKekI5bYj3Ks9UsDUebcYfFq5kwFTm9PtWkcT1rxCxCvUPZeeNOh9G4NmWWyHCJ/SMdt8R6dDRPsPNUwOMxlb6X9l0pcH/VfkYFjrKQQbZBqgQ6DtJxS6xHaRbo3NSqSUciOE4txOmjxKvBaCCvJs+h4tsAPm4+BHkGycwSB0A6bon1KN2jOkEXN60t6R2uHqrdfeyGU1BfQIuhhZGDHMtxgzrrljNu+0dqlUisQ10x5P0KI6drbUnfCBsHh//s06k5VWps3JFCJYsyFrF4x+KO1wlLEwC4K/Euh4jPDzTkjFtiHdY+C4YWqDuitSV9I3Qc1ORDc02vT91fvR8FhdiAWCsYZh1Mi6oPpDwAwKZrNjnMoupApMcZt6IonsBvgEf78V8JIZ6wtmESB+V4lb2idPvuetMdpq7vpXsganKvTs2pziHSLxJvN8t00bElpgXK/Np84oPjNbZG0h3mzLhbgHOFEIlAEnCBoii9+yRLBg4dXW/a49qunvbd9aY7wtozYb65q9dNFXKqchxuYdJEjH8MgNQssXN6dNxCpb79pVv7j7CqVRLHxdG63nRHwDA1I6Yyr1dNFVoNrRyqPeSwjnuY3zAUFLlAaeeYtTipKIoLkA6MBN4RQmzp4pg7gDsAoqKiLGmjxNFoKAXfIeAdCFGn2XfXm67oR1OFg7UH0Qs9owY5zsJkZ9xd3InwjZAzbjvHrMVJIYRBCJEERAKTFEUZ38Ux7wohUoUQqSEhIZa2U2IBDEZBcU0zu4tq+C27jHV7S8ktq+/5xN5y9cdgbIWhqWrXm36o7GlCP5oqmDJKHHXGDVJsyhHoVTqgEKJaUZR1wAVA/1R4JDahRW/Aw9UFo1Ew9YW1FNc2H7N/VmIEb10zESEEt3+YTlyYL6ePDGby8CB0OqVvF60vhcaKoxWIjkY/mirsr96Pq+LaESt2RKL9o8nIzUAIgaL08TMgsSrmZJWEAG3tTtsLmA70v5OqxGoIIdiYU86S3w9SVt/Ct3dPRadTuOPM4bi76gj29SDI1x0XnYKvh/oRqG3WU17fwq/7Sln0ay6Rg7y4OnUYV6YOIzzAs3cGmErFHa3UvTMNpTDmEtj7HYw8z+xwT05VDjEBMbg5WtFRJ6L9o2loa6CiuYJgr2CtzZF0gTkz7iHA0vY4tw74Qgix0rpmSfrK9vwqHv92NzsLawj2defaU6PRGwVuLgq3nN59XnGAlxvf3D2VxlY9a/aU8vmf+by6OpvYEB8unhDROyNMFYeOOuMGNbzT2gDPrYDISXD2I2adtr96PwnBCVY2zrp0ZJbUHJSO207p0XELITKBiTawRdJPfssu48YPthLq58FLcyZwaVIEHq4uvRrD292VWYkRzEqM4HBlI0MD1QYI7208QH2znjvPHt7zmKVZ4BMCPg7+T+/uozaBKDEvlbGxrZHC+kJmj5xtZcOsS3TAUbGp1PBUja2RdIUseXcCmtsMeLq5MGVEEP83YzQ3TInpCIH0h2GDjxaQ7Cuu5Yu0Ar7fWcTLcxJJHBbY/YmlWY4dJulM+Hg4kmnWoR1d3R1Qo6Qz4d7huOvc5QKlHSNL3h0YIQQfbT7ItNfWU1bXgpuLjvlnj7SI0z6el+Yk8sHNp1DbpGf2ot956ae9tBmMJx5oNKh9Jh05TNKZsASoOgAtdT0e6ohdb7rCRedClH+UTAm0Y6TjdlD0BiOPfJ3JP77dTVyYH+4u1v9TnjM6lFX/70zmpETyr/W5ZBZ0oeNRdRD0Tc4z4zbdgE6iFGjqevP4pscBuHj5xQ7T9aY7pEqgfSNDJQ5Iq97Igs+28+OuYu47dyT3T4vre+peL/H3dOOlOYncceZwRob6AVBa20yof3vmiUnDOsxJHHd4e8lCyU6IOrXLQ+YnzWd+0nxu+fkWWvQtLJvpYHnrXRDtH836gvUYjAZcdL1bJ5FYHznjdkDe/CWbH3cV8/jF4/h/54+2mdPujMlpb8op5/SX1vHN9kJ1R2kWoDhen8nuCBgGngFQfPKyBSEE+yr3ETc4zkaGWZdo/2j0Rj1FDUVamyLpAjnjdkDuPGsE8REBXJQwRGtTiB8aQErUIO7/PIPi2mb+UrIbZXAsuDueMl6XKIqqFNhD1/eSxhJqW2sZPWi0jQyzLp3bmA3zG6axNZLjkTNuB8FoFLy38QCNrXr8PN3swmmDmv+95JZTuCQxgvd+3Ezb3p8Rg4drbZZlCRuvxriNXSzGtrOvch8Aowc7n+OW2B/ScTsIL/68l2dWZrFyh/01JvBwdeHNq5P4d+TPuBlbqC23Pxv7Rfh4aGtQs0u6YV+V6rgdqevNyQjyDMLXzZeDNQe1NkXSBTJU4gAs3XSQf6/P47rJUVyZGqm1OSfybCg6fQvJAAoEVO92zOYJ3WFqqlCyq9vGx/sq9zHUdyi+7r42NMx6KIoiM0vsGDnjtnM25Zbz1IrdTBsbxlOzxtun6M8JzRM8qBk1m1fGfoXB6ATS7aFjAQV+/mu3TRWyq7KdJr5tSm/cXbGbzUc2k7A0weHTG50N6bjtGKNR8MS3u4kN9uGNuUm4aJA9YhYnNE9oo6jRlYV/1vG35TsRwsGdt5uXqg5YU9BlU4UmfRP5dflOE9829Z+8J+keALZcu0X2n7QzZKjEjtHpFJbcMonmNoNVqiEtynHNE8bWl3DPOSNZuC6HAC83Hr1wjH1+W+gJM5oq5FTlYBRGp5lxmxgeqC4yH6g9QHyQk1TCOglyxm2nbNxfjtEoGBroxYgQB4ibXv2xWjEZOamjecKD58dxw5Ro/v1bHot+zdXawr5hCgPpuu+haVqYdJYcbhPDA9odd033i7ISbZCO2w5Zt6+U697bwqd/5mttivnUFkJzNYQflTRVFIUnL4nnsqQINuWWd61tYu+YwkBGvfpa33JCU4V9lfvwcfNhqO9QjYy0DlF+UbgoLuRV52ltiuQ47Pz798CjsqGVh7/KZHSYH1ck22EGSXcUt89Aw4/VotbpFF6+MhEhwM0GeipWoaEUkuZBxscQmXpCU4XsqmziBsWhUxz09+sGNxc3hvkNk2JTdohzfdIcHCEEf/9mJ9WNrbx+dRKebg6kEWEqCe9CFdDNRYe7q47Khlauf28Luwq7EKeyZ+Yug8vegYAoCIw6poemEKLDcTsjMQExcsZth0jHbUd8k1HIDzuLeWB6HOMi/LU2p3eU7IRBsWpYoRvaDEZyS+u5ZcmfFFU32dA4CzFkAhzZccymooYi6tvqndZxDw8YzqG6Q+hNoSKJXSAdtx0xJMCLWYkR/OXMros87JrinSeESY4nzN+TJbdMoqnVwG1L02hsdTBnMCQJKnKgubZjk7OVuh/P8IDh6I16CuoKtDZF0gnpuO2IycODeOuaifabr90dLXVQeaBHxw0QF+bH29dOZG9xLQ98noHRkQp0hiSqjyW7OopUFqxbAMB1P1znlEUqsQFqn9K8GhkusSek47YD1mSV8PyPe2huM2htSt8oyQLE0dLwHjh7dCh/mzmOvcV1lDe09HyCvWBy3EUZHUUqZ0eezfCA4ey8cadTFqmYHLdMCbQvpOPWmPoWPf/4dhe/7i1zvJm2iZKuM0pOxi1TY/jhvjMI9fO0klFWwC8M/IYcE+fOqshiXJCTNI3oAj93P0K9QuWM286QjltjXvl5H8W1zTx3eYLjpssV71SbDQSYn76oKAo+Hq606A08/u0udhc5SKbJkMQOx13WWEZpU6lTO25QZ91SJdC+cFBP4RzsKqxh6eaDXHdqNCnRg7Q2p+8U74LwCWrTgV5S16xn1e4S/vJROlUNrVYwzsIMSYTyfdDaSFaF2ofS2cvBYwNiyavJc3zNGSdCOm4Nefb7LAZ7u/PQDAfOSKgphKJtMKhvzROCfT349/UplNa1cM+n29Dbe3XlkEQQRijZTVZFFgoKYwaP0doqqxIbEEt9Wz3lTeVamyJpRzpuDXnm0vG8clUiAV5uWpvSd1Y/rjqyqr7HQBOHBfLPy8bze04FL/+8z4LGWQHTAuX/7iCrZDuxAbF4uzlJm7ZuMIlNyTi3/SBL3jXAaBTodAqjwvwYFdZ9wYpdc7xq3sEN/WqecGXqMDILali2JZ9bz4i130VL/6Hg4gFVeewu0TE5dobWFlmdzmJTpw7putO9xLbIGbcGvLJqH/d84gBhgZPRoZrXfu939TpBNa+3/OPicay893T7ddrPhsJTgWBoocxFR5loZdyfH6rbnZgQrxB83HzkjNuO6NFxK4oyTFGUdYqiZCmKsltRlAW2MMxZya9o5L8bDuDuqsPVUbNI4DjVPAUMJ6rm9RZ3Vx0xwT4IIfgqvYD6FjurrOx0s8pydwdgXOTUft2sHAFFURgeMFw6bjvCnFCJHnhQCLFNURQ/IF1RlNVCiCwr2+aUvPjzXlx0Co9c4AQLWg2loHOHMReCd/AJqnl9Jbuknoe/2sH67DLemptkPw0YOm5WBrI8vNEJwRivsH7drOydRRmLWLxjccfrhKVqrv5diXc5XbGRI9HjlE8IcUQIsa39eR2wB3Au4WEbsT2/iu8zj3D7GbGE+dtpOKA3TH8ajK0wclpH8wRLMDrcjwfPH82KHUV89IedNattKIXEuez28CAWN7wbKrS2yKqYKkQfSn0IgA1Xb3DKClFHo1ff1RVFiQEmAlu62HeHoihpiqKklZWVWcY6J+OtX/YT7OvBHWc5oIhUVxRtVx8jJlp86LvOGsG5Y0J5ZmUW2/OrLD5+n5m7DGb/iyxPT8a5BVjsZmXvjAocBcD+6v0aWyKBXjhuRVF8ga+B+4UQtcfvF0K8K4RIFUKkhoSEWNJGp+GVKxNZfF2y/fePNJei7WorrxDLh310OoXXrkokzN+Tez7Zblc6LqWNpZTpFOJrK7U2xWaY2rJlV2VrbIkEzEwHVBTFDdVpLxNC/M+6JjkfBqNAp0CQrwdBvh5am2M5irar+iQu1slDD/R2Z9G8ZCrqW+2qqYSpYnJc9RFoKAefYI0tsj5BnkEM8hjE/io547YHzMkqUYD3gD1CiNesb5Lz8dmf+Vy2aBPVjQ5Q0m0uRoOq2WGFMElnJkQGcs4YNd2upLbZqtcyl6yKLHQojG5tg8J0rc2xCYqiMGrQKOm47QRzQiVTgeuBcxVFyWj/ucjKdjkNzW0G3vplP646xbErJI+nIgda663uuE2sySrhjJfW8UeedouBJg3uxTsWY0RwaswwEjY/5HQa3N0RNyiO/dX7MQoHrj9wEszJKtkohFCEEBOEEEntPz/YwjhnYOmmg5TUtvDwjNH2k9ZmCay4MNkVk0cEERnoxb2fbqe0TpuZ9/yk+WRcn4Gvmy9Xj76anY0B7NSNHDAZFqMGjaJJ30RhXaHWpgx4HLgCxP6paWpj0a+5nBUXwqnDg7Q2x7IUbQc3bwi2Ta9FXw9XFl2XTF1zG/d/loFBo845uTW51LfVkxiSCJEpaqjEODBmoKbMkuxquUCpNdJxW5GP/zhETVMb/+fI6n/dUbRdFVzS2W7RcEy4P09fOp5NuRW8+Ys2sdaM0gwAkkKSYGgqNNdAZa4mttiaEYEjUFBkZokd4CR5afbJbWfEMnaIH+OHBmhtimWpLoDDWyH5Rptf+qrUYWw7VIW3uzZZJjvKdjDYczCRfpEQmapuLEiD4FGa2GNLvN28GeY3TC5Q2gHScVsJIQQeri6cO8YJy6FX/x0QUKVNH8LnL0/QbL1gR9kOkkLay/CD48DNB9Y8CSPOderSdxMys8Q+kKESK1BS28x5r61nc66TlUM/G6pKt+5err4+sF59bWN1PJPT/nVfKXd8mGYzlcXK5koO1R4iKTRJ3aBzAU9/qC+G9S/axAatGTVoFPl1+TTr7SM1c6AiHbcVWLQuh/yKRoYGemltimUxqeMp7WEKC0i59oeapjZWZZXw8irbNF/YUar2mkwMSTx6E6s7ou5Me0+Tm5itiRsUh1EYya0ZGHF9e0U6bgtTWN3Ep1sPc2VqJFFBTtYZxaSOJwyg6Cwi5dofLk0ayrxTo/j3+jxWZ1lGmfBk7CjbgavOVW0ObLqJuajyrrh4aHoTsxUdmiUyXKIp0nFbmIVrcwC451wnXayqbc/hTb0VUm62mJRrX/nHxeMYP9SfB7/I4HBlo1WvlVGWwbjB4/B09Tx6EzO0qTs1vonZimF+w/B08ZSZJRojHbcFKaxu4su0w8ydNMz5wiQmTrlVfRx3qUWlXPuKp5sLi65NQQBfpRdY7TptxjZ2le8iMTTx6MaGUki9RV2k9B2i+U3MFrjoXBgROELOuDVGZpVYkIgATxbNSyZxWKDWpliP/D/UdmVDU7S2pIOoIG9+uO8MIgdZ72a5r3IfLYYWNb5twnTT+ukxSHsf5nxgtevbE3GD4lhfsB4hhHNVAzsQcsZtQRRF4fz4cOdoktAdh7dA+ARwt6/4/bDB3iiKQm5ZvVXi3abCm2Mct4no00DffFQGwEkxabUsz1lOZXMlEz6cQMLShAGj1WJPSMdtIZ5asZuFa53866OhXQ0varLWlnTL8z/s4d5Pt7G3+ATJ+D5hclYv/qmm+03/avqJzipqivp4aKNFrmmvmLrhfHLRJwC8dvZrshuORkjHbQEOVTTw4eZDVDQ4kWxrVxzJVGeWw07V2pJuee7yBPw93bjzo3Rqm9v6PZ5JWMrf3Z/ZI2ez88adJzorn2C1mcShTf2+niMwevBoXHWu7CrfpbUpAxbpuC3A22tzcNUp3OUsLcm64/Af6qMdO+5QP0/emZdMQVUTD36xA6MFxKj2Vu2ltrWWU4ec5PeOngr5W8BgZ53prYC7iztxg+LYXbFba1MGLNJx95OD5Q0s317IvFOjCXXm2Dao8e3AKPAforUlJ+WUmMH8beZYVmeV8Hna4X6Pt/XIVgAmhU/q/qDo06C1DkqcO4/bRHxQPFnlWVKbWyOk4+4nC9eps+07zx6utSnWpfYI7P0BwrtYnLNDbjothhcuT2D2xKH9HmvLkS2MCBhBiPdJeqlGn6Y+fnkT1Dl/WuD44PHUtdWRX5uvtSkDEum4+8mNU2J45rLxhPo5+Wx79T/A2AaNjqG/oigKcydF4enmQk1TG4XVTX0ap83QxrbSbScPkwD4R4C7H1QdHBC6JfFB8QDsqpBxbi2Qedz9JCEygIRIJ5Nt7cyzoaBvOfo6f5OqyeHqAX8v1c4uMxFCcNMHW2lqNfC/+afh7d67j/yOsh006ZuYNOQkYZLj36O099QfB3mP+sKIwBF4uniyu3w3Fw+/WGtzBhxyxt1HDpY38NCXO+ymga3VsDNhqd6iKAoPTIsju6SO//sqEyF6t1i5tXgrOkXHKeGndH9Qh25Je0/RAaBb4qpzZczgMXKBUiOk4+4jC9flsGJHEU5fOOYXDh6+diMs1RfOjAvhkQvG8H3mERav752q3ZYjWxg3eBz+7v7dH9ShW9KeUWJodbj3qC+MDx7Pnoo96I3On0ljb0jH3QdMmSTXTY52/tg2qHFbgNMfsAthqb5wx5nDuSQxgpd/3sfavebZ39jWSGZZZs/xbTiqWxI0CnxDHfI96i3xwfE0G5rJrZYSr7ZGxrj7gClv+y9nOXkmiYkxMyFvHSTNgyDHzFVXFIWX50zA18OF8RHmrUlsK92GXuhPHt82YdIt+eVp2PgGXOb8ZeCmBcqsiixGD3bCvqp2jJxx95K8snqWby/g+oEy2wbI+1XN3x7s2DcqTzcXnr98AqH+nugNxh4rK/8o+gM3nRsTQyeaf5GR09SwUt6v/TPWAYj2j8bXzVdWUGqAdNy9xNfTlesnR/MXZ6+SNGHQw4ENMPxsnCWgL4Tg7k+2cduSNFr0hhP2m/RJlmYtpc3YxqRlk8wXU4qcBB4BkLPGCpbbFzpFR3xQvEwJ1ADpuHtJqJ8nT106nhA/D61NsQ1F26GlBoafo7UlFkNRFGZOiGDrwUoe/XrnCZkm85Pm8+1l3wLw11P/2rU+SXe4uMLws2D/GuhlBosjYbq5bSneQlZFFglLE6RSoA2RjrsXLP41l/RDlVqbYVvyfgUUiD1La0ssyqzECB6cHsfy7YW89UvOCfvX5a8D4JxhfbhhjZwGdUVQuqe/ZtotJqXAxdMWA/Du9HelUqANkY7bTPaX1PHSz3tZs8c5Cyq6pK4YNr0FoWPBJ0hrayzOPeeO5IrkSF5fk8032wuP2bfu8DrGDh5LuE947wceOU19/HSu05e/TwydiIviwp/Ff2ptyoCiR8etKMr7iqKUKooyoANZr6/JxtvNhdvPcOwFul6x9p/QUgs6F60tsQqKovD85QlcPGEII0J8O7aXN5WTWZbJOVF9DA8FDAXPQVB9yOnL333cfIgPipeO28aYkw64BFgIfGhdU+yXXYU1/LCzmPvOG8VgH3etzbE+x5dwF+90qDL33uDuqmPhtckdrysbWllfuB6B4Nxh5/Z+wAFY/n5K+Cks3b2UxrZGvN3sqzOSs9LjjFsI8RswwAK7x/Lqqn0EeLlx2xmxWptiG0wl3DrHLHPvK/9an8sFb/zGD3lriPCJIG5QXO8H6Sh/b1+8dnF3+vfulPBT0As9GWUZWpsyYJAx7h4QQnBK7GAePD8Of083rc2xDaYyd6Pjlrn3hXPHhNKkb2Rr8RYmh5/Zt0a4pvJ3Y3uO+AAofzfFudOK07Q2ZcBgscpJRVHuAO4AiIqKstSwmqMoCvPPHqm1GbanvL1/5jl/VbW4B0AJd1yYH3fPNPLO7jY27gijNrWtbzfrhlJVGqC5Bnb9D2oLez7HgfF28yY+WMa5bYnFZtxCiHeFEKlCiNSQkJMIzjsQaQcr+Taj0CLtrxyOyBTQucIpt8PFrx0t6XZyDjVtwcvFl/yiMG5bmta3v/3cZep7Nul2wAjjr7C4nfbGKWGnsKt8F41tjVqbMiCQoZJuMBoFT63I4sUf99JqGGDtmYSArO/U3G2vQK2tsQmmgpKVeStpMtTjGfdX9njewb8yF/d90MhJ4BcBu5dbzlA7ZVL4JBnntiHmpAN+CmwGRiuKUqAoyq3WN0t7Vu48ws7CGh48fzSebs6ZDtctJbuh6gCMm6W1JTZjftJ8npn6DAAfXfjRMdWS2/KruiyN7xGdDuJnq+XvTdUWtti+SApNwlVxlXFuG2FOVsk1QoghQgg3IUSkEOI9WximJa16I6/8vI8x4X5cZoGehQ7Hnu8ABUbP1NoSm/JtzrdE+0eTGHK0r2ZRdRNz3/2Dez7ZTqu+D9+84merC5T/nebUxTgyzm1bZKikCz7Zcoj8ykYeuXAMLjrnEFYym7pi2PQ2RKaCr3OsVZhDQV0BaSVpzBox65hskohAL/7e3jH+3k+39d55R6aCuw9U7HfaYhxTmGlH2Q4yyjKkbokNkHrcXRAR6MWVKZGcHTdwHFcHP/8N2hqPtiobIKzIXYGCwiXDLzlh3w1TYjAaBU+uyOLOj9NZNC/ZvPDZACnGmZ80n/lJ88ksy2TeD/N47vTnuGTEie+jxHLIGXcXnB8fzstXJvYtj9dReTZUrY7c9ZX6+vAf6utnQ7W1ywYYhZFvc79l0pBJDPEd0uUxN02N5bnZCazbV8qX6QXmDXx8MY7OzamLccYHjyfUK5S1+Wu1NsXpkY67E4crG3lnXQ7NbX1YiHJ0FmTC+DlA+81qgFRLAmwr2UZhfSGXjrj0pMdde2oUX905hetONbNO4ZhiHEV9dPdz2mIcnaLj3Khz2Vi4kSZ9k9bmODXScXfi+R/3sHBtDlWNrVqbYnv8wqG1HhDqzHCAVEsCLM9Zjo+bD+dFndfjsSnRg1EUhfyKRub99w+KqntwUKZinHP/rr4uzbKAxfbLedHn0WxoZlPRJq1NcWqk427nj7wKfthZzF1nj2BIgJfW5mjDkUxVW+OWnxy2KXBvMC2qfZf7HQ1tDZz6yalmL6qV1DWTebiGKxZvYn9JXfcHmopxptwD3kHgE2zB38D+SAlLIcAjgF8O/aK1KU6NXJwEDO3FNkMDvbjjzAEk29qZhnJoKINTblMzISJTtbbI6sxPmk+zvpmlWUv5fvb3RPpFmn3uKTGD+ewvk7npgz+5YvEm/nVdCqeNPIlTdvOEidfD72/Cf86DuZ845bcZN50bZ0WexbrD62gztuGmGyD6PjZGzriBz/7MZ8+RWv560diBV2xjImOZGoNNuUlrS2xGbWstX2R/wYzoGb1y2ibiIwL4312nEebvyQ3vb2V9dtnJT0i9GRBQmOa0qYEA06KmUddaJ3O6rYh03EBiZCA3nRbDRQl96HbiDNQegXXPwdAUCB2jtTU244t9X9DQ1sDN42/u8xjDBnvz9fzTuGFKDKnRg7o/8NlQePNoYQ9p7zlt1s6UiCl4uXrJcIkVkY4bGD80gCdnxQ+s9L/OfP8A6JvVQpEBQouhhWV7lnFaxGmMDRrbr7H8Pd14/JJx+Hi40tCi57H/7aS8vuXYgzpSA9sbcbg4b2qgp6snpw89nbWH12IwDsAMLRswoB33b9ll/N+XO6hrbtPaFG0w5W7v+1F9feA3p50FHs+K3BWUN5X3a7bdFTsKqvl6WwGXvL2R7flVR3d0pAbqAQUMbeqN0sni3KYF39WHVlPeVE7SR0myitIKKEJYXrI0NTVVpKXZt9hMU6uB899Yj5tOx4/3n4GH6wCMbdcVwzfzIbf9K62rF4y9GM7/p9M5FBOLMhaxeMeJin93Jd5lsQ7luwpruPPjdEpqm3nkgjHcMjUWnU6Bz+aBbxiEjVe/5YTGw3znTJvTG/Vc8PUFxAbE8p/z/6O1OQ6BoijpQgizsgIGbFbJ62uyOVzZxKe3Tx6YThvUWWDZXvW5i8eAyN02OefFOxbzwYwPSA23fPbM+KEBrLz3dP7vq0ye/X4P5fWtPHrhmKOa5kLA7v+pKozvzYCrPnS699xV58rVo6/mre1vkVedx/DAAZqtZSUGZJ7QA/0AAA4rSURBVKjkj7wK/rMhj2smRTFlRJDW5mhDXTH863S1O8vQVLj9lwGRu11YX8j7u97nwtgLreK0TQR6u/Pu9Sk8e9l45rVXWnZIwyoKnPsPaKqEw1ucNsPkirgrcNO58eneT7U2xekYcKESIQQz39pIY6ue7+87Ax+PAfqlY+X/axc88oSH8wbMwuQD6x7g96Lf+e6y7wj3sV0WkRCC2z9Mw8fDlTf2z0AxtJx4kJOJTwH8bePfWHNoDb9c+Qu+7r5am2PX9CZUMuBm3Iqi8N8bU1l8XcrAdNqmBcm0dll1fTM8FzEgFiQ3F21mTf4a7phwh02dNoBRqCGU7zOPcAHvcHjoTIRJfEpxcdoMk2vHXEujvpFvc7/V2hSnYkA57ryyeoxGQUSgF2OH+GttjjYsyGzvgTiwxKTqWut45o9nGOY3jOvHXW/z67voFO6fFse390zFPXAI6w81IwytCMUFhEHtkONkcW6A+OB4JgRP4NO9n2IUA6wFoBUZMI57f0kdl7y9kRd/2qu1KdriFw4VuQwkMal3tr/DaZ+exuG6wxyuO0zqx6mapajFRwTwzd1TOSNC8J3rBdRf9xN4DUYc+A1K98IHFzpNpxxTamBmeSaHag+R+GGiTA20EAMiVlDT1MbtH6bh7eHKzVNjtTZHO+qKYdlVULwDgkbAlUsh7QOnX5D091C/XT2U+hA3xt+osTXq7Dt6/nIijQIXnYL+um/gP+fQ8MEc/JsLUda/qApTOTimBgtCCG766SYO1BxgxewVBHgEaG2aw+P0M26DUbDgs+0UVjexeF4y4QGeWpukHb8+rzptN2+47RcIT1AdhClNzQnJKM3gtbTXOHfYudww7gatzTkGU1s8l/en44qBgKbDKMLodOXwiqLw2KmPUdNawzsZ72htjlPg9I77tdX7+HVfGU/Oiic1ZrDW5miDaUEyfYn6uq0RXoxxGsfQHbnVudy/7n7CfcJ55vRn7FbSQLk/ExF/BaL937FFuLJOOZXmoHinCZuMGTyGq+Ku4vN9n7Ovcp/W5jg8Tu+4zx0Typ1njWDeqdFam6INdcUQlgAhnfQ4BsCC5DObn+Gyby+jormCgvoCpn461X7jq37hKJ4BKIpAoOCu6Bnvko9HaQasf5Hfc8o5XNmotZX95p6J9xDgHsBzW56TC5X9xGlj3LuLaoiPCCAlejAp0QN0pg2w/iUoTAfa8/Vdnb9CMrsqmzX5awjxCuG9Ge8RG+AA6xoNpZByC8r2D8HQRoj+iLo97T2mpr1Hs3DjppgfmHdqNOeMDsHVxbHmXJ2lBqpKq0j8UFVKtKTUwEDCKQtw3tt4gGdWZvGfG1KZPs45nVOPHN9h3ISLG0y8QV2QdMLY9k8Hf+KpTU/h7ebN+zPeJ9rfwb5p1RXDigWQ/ZP62sWdxpjzqC45zK3N97On3ptgX3f+PnMcl00cqq2tfUAIwT+3/JPP933O30/9O1ePuVprk+yGAatV0mYw8tSK3Xz8Rz4XxIdz9ugQrU3SBlN4pLkaKnLUbU4uINWkb+LGH29kT+UeAOrb6rl4+cWAg83q/MLBL4KOPHtDK95H/sS7qZKVyZvYMORmItbcTY3bu8BQsopq+XpbATPiw0mJHtSx4GmvKIrCo5MepbihmOe2PkeYTxhnDztba7McDqeZcVc1tDJ/2TY251Vw51kj+L8Zo+3+Q2xR6orhq5thzhJY8yTs+OToPlcPVUY05WanSDPrjFEYWX1oNW9vf5v82nxuGX8Ld0+827FbZplUBLd/BIZuGlen3goXv8aXaYf52/JdtBqM+Hu6cvqoYM4cFcJlE4fadTenxrZGbv75ZvKq83jytCeZOXym1iZpTm9m3E7juNftLeUvH6XzwhUJXJ7c+zZUDklnZ73+xaNl7MfjhOERg9HA70W/88SmJyhvKj9hv0PNsrujrhh+/hvs+Va98XaJQsO9u9mcU07s+nu5p+1eDjT7kvnEDNxddXyzvZCapjYmRAYwdoi/XTjz7qR1bxt/GwtSFmhgkX0wYBz3gfIGMguquTRJjfUdLG8gJtgJxZI6O2iEec4anC48IoQgtzqXp/94mu2l20/Y7xTO+nhWPADbloDiCsZOs29FBz6h6s049RZ1W/oHiJSbOJJ0HxGr58OcJdz32TbmHX6Se1rvo9plEKNC/TgjLpjHLlSzjErrmgny8dDs26neqOdfO/7Fu5nvEu0fzZ2Jd3J+zPmO/Y2pj1hcZEpRlAsURdmnKEqOoiiP9s+8k1BXfLTkt5vnhpojVL0zjQX//YmrX11O5DdzaKwshLpiYlbMOem5Nn1u5u9j1vP1L0H+H6qjXv8SHNoEr8Z147Tb/6ROoK/d2NbIzrKdfJX9FY9teIxpX05j9nezySzL5MzIM3n1rFdJvy6dnTfuZOeNO53PaUN7tsnNcMdaCO7UD1QYob4YEOrnIO09EEaUtPeJ+G+S+hlZ/yJvDlnFJN0+vh73GwtO9ee1xscw1BR3fL5uXfg9Z/7jU3Y8cxoL/vsTL375K5ULz+v4DDa/O4PasgJE7RGr/L+8++GZ/Dvz3wgEB2sP8uiGR0n+KJlbv7+OfR+ch7H2iH38P/fl/9+K9DjjVhTFBcgGpgMFwJ/ANUKIrO7O6fOMe+X/g/QPIOVmBALSllAbfx2e7i54ZCwlP/Zq/sirYI5Yzf905zMqzJcJxctRUtvbT7WfaxfPL37tmN+n6+NuUkX1ty2FpHnqP2PGJzD+cvWr8Z7vzHzjFEAcddYhY+CK/x4tZ7dxeEQIgVEYMQojeqHHYDSgN+ppNbbSYmihzdBGo76RxrZGGtoaqGmtoaalhqrmKkoaSyhuKOZIwxEK6wt7vJZTzrK7wxT7HnuJ+tmqOgj0IR96/BWqIuHOL8kbdhk1TQYSy1fwk/sM6loNXMkadKm3YBQC0j5gmeE8dDqFa3Rr+N79AkaF+jKm8Gv0yTex/VAVKeXfkBUxBzcXiDv8FTXx1xHo5Y5I/+D/t3dmsVGVYRh+3plSFq1SSkyURZYUY8WlSBRvRINRQwzESAwaFg1KiuCNd4QbohfqhV6YkCgXBNwQ9cI0UTRRICTE4gayRRNAotUGGkAwJuBM+3lxTklTCnOGnqWnfE8y6f/P+Wf6vvOf+eacf6Vz2tMUBA2/fMDZ2xYjQd2B99HMZ+k2Qz9upNS8FAkKezaxuvFutpYuDnx3nTvPYw13clNhJKN//ZK6pse5VjWM2LuZ2uYlDEPop42D6/tfJbE2lUi6D1hrZo+E+dUAZvbqpV5TdeAOh67NnjiOc4N0dtugR4Vgbe2uEmBBugqMy58HEATk3uV78j1pIwjYUd6rWlruaGFl88rY3ze39DShFGuDpXkhaE6xcE/LBOogTTqKRV5rqGfbNaMiv6ZgwdzTgoEwxIWxORdIOrqM6epia3s4Br/K9dXjDtwLgEfN7Lkwvxi418xW9Sm3HFgeZm8BIs9rHd4wbLzqivm8n3cGTPlsuaN8qvxX1joSZixwcS/qFdI4pjC11G2lzn+tc0p9YSrA0dPdRybXFyYPLzLSwBTGqZ5037+9jw2atNFdNso1okai0F+5zmKREznYbtD+6Tp+/mSpvYqX3GxmkcYwxzaO28zWA+sH+j6Sfoj6qzNUcM9Dn6vNL7jnJInSOfknMKFXfnz4nOM4jpMBUQL390CjpMmSaoGFQNReM8dxHCdmKjaVmFlZ0irgK6AIbDCzgwlqGnBzSw5xz0Ofq80vuOfESGQCjuM4jpMc+Vob0nEcx/HA7TiOkzcyC9yVptFLGi5pS3h8t6RJ6auMjwh+X5J0SNI+Sd9IytlC0hcTdakESU9IMkm5HzoWxbOkJ8O6Pijpw/7K5IkI5/ZESdsl7QnP77lZ6IwLSRsknZB04BLHJemt8PPYJ2lG7CLMLPUHQSfnEWAKUAv8DDT1KfMC8HaYXghsyUJrin4fBEaF6RV59hvVc1iuDtgJtAEzs9adQj03AnuA+jB/Q9a6U/C8HlgRppuAY1nrHqDn+4EZwIFLHJ8LbCWYODQL2B23hqyuuO8BDpvZUTP7D/gImN+nzHxgU5j+FJijwbrba2Uq+jWz7WbWs7FgG8F4+TwTpY4BXgFeB86lKS4honh+HlhnZqcBzCz6nOjBSRTPBlwXpq8Hcj1L1sx2AqcuU2Q+8K4FtAGjJd0Yp4asAvc44I9e+fbwuX7LmFkZOAM0pKIufqL47c0ygl/sPFPRc3gLOcHMPk9TWIJEqedpwDRJuyS1SXo0NXXJEMXzWmCRpHbgC+DFdKRlRrXf96oZUluXDQUkLQJmArOz1pIkkgrAm8AzGUtJmxqC5pIHCO6qdkq63cz+zlRVsjwFbDSzN8JF696TNN3Mt3q/UrK64o4yjf5CGUk1BLdYJ1NRFz+Rlg2Q9BCwBphnZv3s9JsrKnmuA6YDOyQdI2gLbM15B2WUem4HWs2sZGa/ESyZ3JiSviSI4nkZ8DGAmX0LjCBYdGuokvgyIVkF7ijT6FuBpWF6AbDNwpb/HFLRr6Rm4B2CoJ33dk+o4NnMzpjZWDObZGaTCNr155lZunvexUuU8/ozgqttJI0laDo5mqbImIni+XdgDoCkWwkCd2eqKtOlFVgSji6ZBZwxs45Y/0OGPbNzCa42jgBrwudeJvjyQlC5nwCHge+AKVn3Jifs92vgOLA3fLRmrTlpz33K7iDno0oi1rMImogOAfuBhVlrTsFzE7CLYMTJXuDhrDUP0O9moAMoEdxBLQNagJZedbwu/Dz2J3Fe+5R3x3GcnOEzJx3HcXKGB27HcZyc4YHbcRwnZ3jgdhzHyRkeuB3HcXKGB27HcZyc4YHbcRwnZ/wPuf4Q4oI1L0wAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "y5 = [beta_pdf(i, 4, 8) for i in x]\n", "y6 = [beta_pdf(i, 23, 27) for i in x]\n", "y7 = [beta_pdf(i, 33, 17) for i in x]\n", "\n", "plt.plot(x, y5, '--', label = 'Beta(4, 8)')\n", "plt.plot(x, y6, '-*', label = 'Beta(23, 27)')\n", "plt.plot(x, y7, '_-', label = 'Beta(33, 17)')\n", "plt.legend(loc=0,numpoints=1,fontsize=10)\n", "plt.ylim(0, 6)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "\n", "\n", "What’s interesting is that this allows us to make probability statements about hypotheses: \n", "\n", "> “Based on the prior and the observed data, there is only a 5% likelihood the coin’s heads probability is between 49% and 51%.” \n", "\n", "This is philosophically very different from a statement like \n", "> “if the coin were fair we would expect to observe data so extreme only 5% of the time.”" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "If you flipped the coin more and more times, the prior would matter less and less until eventually you’d have (nearly) the same posterior distribution no matter which prior you started with.\n", "\n", "Using Bayesian inference to test hypotheses is considered somewhat controversial— in part because of the subjective nature of choosing a prior. " ] } ], "metadata": { "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Python 3", "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.7.6" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autoclose": false, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": false, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": { "height": "calc(100% - 180px)", "left": "10px", "top": "150px", "width": "165px" }, "toc_section_display": true, "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 2 }