{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Probability Theory Review" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Preliminaries\n", "\n", "- Goal \n", " - Review of Probability Theory as a theory for rational/logical reasoning with uncertainties (i.e., a Bayesian interpretation)\n", "- Materials \n", " - Mandatory\n", " - These lecture notes\n", " - Optional\n", " - Bishop pp. 12-24\n", " - [Ariel Caticha, Entropic Inference and the Foundations of Physics (2012)](https://github.com/bertdv/BMLIP/blob/master/lessons/notebooks/files/Caticha-2012-Entropic-Inference-and-the-Foundations-of-Physics.pdf), pp.7-56 (ch.2: probability)\n", " - Great introduction to probability theory, in particular w.r.t. its correct interpretation as a state-of-knowledge.\n", " - Absolutely worth your time to read the whole chapter, even if you skip section 2.2.4 (pp.15-18) on Cox's proof.\n", " - [Edwin Jaynes, Probability Theory--The Logic of Science (2003)](http://www.med.mcgill.ca/epidemiology/hanley/bios601/GaussianModel/JaynesProbabilityTheory.pdf). \n", " - Brilliant book on Bayesian view on probability theory. Just for fun, scan the annotated bibliography and references.\n", " - [Aubrey Clayton, Bernoulli's Fallacy--Statistical Illogic and the Crisis of Modern Science (2021)](https://aubreyclayton.com/bernoulli)\n", " - A very readable account on the history of statistics and probability theory. Discusses why most popular statistics recipes are very poor scientific analysis tools. Use probability theory instead!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### [Data Analysis: A Bayesian Tutorial](https://www.amazon.com/Data-Analysis-Bayesian-Devinderjit-Sivia/dp/0198568320)\n", "\n", "- The following is an excerpt from the book [Data Analysis: A Bayesian Tutorial](https://www.amazon.com/Data-Analysis-Bayesian-Devinderjit-Sivia/dp/0198568320) (2006), by D.S. Sivia with J.S. Skilling:\n", "\n", "

\n", "\n", "- Does this fragment resonate with your own experience? \n", "\n", "- In this lesson we introduce *Probability Theory* (PT) again. As we will see in the next lessons, PT is all you need to make sense of machine learning, artificial intelligence, statistics, etc. \n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Example Problem: Disease Diagnosis\n", "\n", "- **Problem**: Given a disease with prevalence of 1% and a test procedure with sensitivity ('true positive' rate) of 95% and specificity ('true negative' rate) of 85% , what is the chance that somebody who tests positive actually has the disease?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- **Solution**: Use probabilistic inference, to be discussed in this lecture. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### The Design of Probability Theory\n", "\n", "- Define an **event** (or \"proposition\") $A$ as a statement, whose truth can be contemplated by a person, e.g., \n", "\n", "$$𝐴= \\texttt{'there is life on Mars'}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- If we assume the fact $$I = \\texttt{'All known life forms require water'}$$ and a new piece of information $$x = \\texttt{'There is water on Mars'}$$ becomes available, how _should_ our degree of belief in event $A$ be affected (if we were rational)? " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "- [Richard T. Cox, 1946](https://aapt.scitation.org/doi/10.1119/1.1990764) developed a **calculus for rational reasoning** about how to represent and update the degree of _beliefs_ about the truth value of events when faced with new information. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- In developing this calculus, only some very agreeable assumptions were made, e.g.,\n", " - (Representation). Degrees of rational belief (or, as we shall later call them, probabilities) are represented by real numbers.\n", " - (Transitivity). If the belief in $A$ is greater than the belief in $B$, and the belief in $B$ is greater than the belief in $C$, then the belief in $A$ must be greater than the belief in $C$.\n", " - (Consistency). If the belief in an event can be inferred in two different ways, then the two ways must agree on the resulting belief." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "- This effort resulted in confirming that the **sum and product rules of Probability Theory** [(to be discussed below)](#PT-calculus) are the **only** proper rational way to process belief intensities. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- $\\Rightarrow$ Probability theory (PT) provides _the_ **theory of optimal processing of incomplete information** (see [Cox theorem](https://en.wikipedia.org/wiki/Cox%27s_theorem), and [Caticha](https://github.com/bertdv/BMLIP/blob/master/lessons/notebooks/files/Caticha-2012-Entropic-Inference-and-the-Foundations-of-Physics.pdf), pp.7-24), and as such provides the (only) proper quantitative framework for drawing conclusions from a finite (read: incomplete) data set." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Why Probability Theory for Machine Learning?\n", "\n", "- Machine learning concerns drawing conclusions about model parameter settings from (a finite set of) data and therefore PT provides the _optimal calculus for machine learning_. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- In general, nearly all interesting questions in machine learning can be stated in the following form (a conditional probability):\n", "\n", "$$p(\\texttt{whatever-we-want-to-know}\\, | \\,\\texttt{whatever-we-do-know})$$\n", "\n", "where $p(a|b)$ means the probability that $a$ is true, given that $b$ is true." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Examples\n", " - Predictions\n", " $$p(\\,\\texttt{future-observations}\\,|\\,\\texttt{past-observations}\\,)$$\n", " - Classify a received data point $x$ \n", " $$p(\\,x\\texttt{-belongs-to-class-}k \\,|\\,x\\,)$$\n", " - Update a model based on a new observation\n", " $$p(\\,\\texttt{model-parameters} \\,|\\,\\texttt{new-observation},\\,\\texttt{past-observations}\\,)$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Frequentist vs. Bayesian Interpretation of Probabilities\n", "\n", "- The interpretation of a probability as a **degree-of-belief** about the truth value of an event is also called the **Bayesian** interpretation. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- In the **Bayesian** interpretation, the probability is associated with a **state-of-knowledge** (usually held by a person, but formally by a rational agent). \n", " - For instance, in a coin tossing experiment, $p(\\texttt{tail}) = 0.4$ should be interpreted as the belief that there is a 40% chance that $\\texttt{tail}$ comes up if the coin were tossed.\n", " - Under the Bayesian interpretation, PT calculus (sum and product rules) **extends boolean logic to rational reasoning with uncertainty**. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "- The Bayesian interpretation contrasts with the **frequentist** interpretation of a probability as the relative frequency that an event would occur under repeated execution of an experiment.\n", "\n", " - For instance, if the experiment is tossing a coin, then $p(\\texttt{tail}) = 0.4$ means that in the limit of a large number of coin tosses, 40% of outcomes turn up as $\\texttt{tail}$. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- The Bayesian viewpoint is more generally applicable than the frequentist viewpoint, e.g., it is hard to apply the frequentist viewpoint to events like '$\\texttt{it will rain tomorrow}$'. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- The Bayesian viewpoint is clearly favored in the machine learning community. (In this class, we also strongly favor the Bayesian interpretation). " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Aubrey Clayton, in his wonderful book [Bernoulli's fallacy](https://aubreyclayton.com/bernoulli) (2021), writes about this issue: \n", " > “Compared with Bayesian methods, standard [frequentist] statistical techniques use only a small fraction of the available information about a research hypothesis (how well it predicts some observation), so naturally they will struggle when that limited information proves inadequate. Using standard statistical methods is like driving a car at night on a poorly lit highway: to keep from going in a ditch, we could build an elaborate system of bumpers and guardrails and equip the car with lane departure warnings and sophisticated navigation systems, and even then we could at best only drive to a few destinations. Or we could turn on the headlights.” " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- In this class, we aim to turn on the headlights and illuminate the elegance and power of the Bayesian approach to information processing. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Probability Theory Notation\n", "\n", "##### events\n", "- We define an **event** $A$ as a statement, whose truth can be contemplated by a person, e.g.,\n", "\n", "$$A = \\text{`it will rain tomorrow'}$$\n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- We write the denial of $A$, i.e. the event **not**-A, as $\\bar{A}$. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Given two events $A$ and $B$, we will shortly write the **conjunction** \"$A \\wedge B$\" as \"$A,B$\" or \"$AB$\". The conjunction $AB$ is true only if both $A$ and $B$ are true. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- We will write the **disjunction** \"$A \\lor B$\" as \"$A + B$\", which is true if either $A$ or $B$ is true or both $A$ and $B$ are true. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note that, if $X$ is a variable, then an assignment $X=x$ (with $x$ a value, e.g., $X=5$) can be interpreted as an event. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "##### probabilities\n", "\n", "- For any event $A$, with background knowledge $I$, the **conditional probability of $A$ given $I$**, is written as \n", "$$p(A|I)\\,.$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- All probabilities are in principle conditional probabilities of the type $p(A|I)$, since there is always some background knowledge. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "##### Unfortunately, PT notation is usually rather sloppy :(\n", "\n", "- We often write $p(A)$ rather than $p(A|I)$ if the background knowledge $I$ is assumed to be obviously present. E.g., $p(A)$ rather than $p(\\,A\\,|\\,\\text{the-sun-comes-up-tomorrow}\\,)$." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- (In the context of variable assignments) we often write $p(x)$ rather than $p(X=x)$, assuming that the reader understands the context. " ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- In an apparent effort to further abuse notational conventions, $p(X)$ denotes the full distribution over variable $X$, i.e., the distribution for all assignments for $X$. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- If $X$ is a *discretely* valued variable, then $p(X=x)$ is a probability *mass* function (PMF) with $0\\le p(X=x)\\le 1$ and normalization $\\sum_x p(x) =1$. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- If $X$ is *continuously* valued, then $p(X=x)$ is a probability *density* function (PDF) with $p(X=x)\\ge 0$ and normalization $\\int_x p(x)\\mathrm{d}x=1$. \n", " - Note that if $X$ is continuously valued, then the value of the PDF $p(x)$ is not necessarily $\\le 1$. E.g., a uniform distribution on the continuous domain $[0,.5]$ has value $p(x) = 2$.\n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Often, we do not bother to specify if $p(x)$ refers to a continuous or discrete variable. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Probability Theory Calculus\n", " \n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Let $p(A|I)$ indicate the belief in event $A$, given that $I$ is true. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- The following product and sum rules are also known as the **axioms of probability theory**, but as discussed above, under some mild assumptions, they can be derived as the unique rules for *rational reasoning under uncertainty* ([Cox theorem, 1946](https://en.wikipedia.org/wiki/Cox%27s_theorem), and [Caticha, 2012](https://github.com/bertdv/BMLIP/blob/master/lessons/notebooks/files/Caticha-2012-Entropic-Inference-and-the-Foundations-of-Physics.pdf), pp.7-26)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- **Sum rule**. The disjunction for two events $A$ and $B$ given background $I$ is given by\n", "$$ \\boxed{p(A+B|I) = p(A|I) + p(B|I) - p(A,B|I)}$$\n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- **Product rule**. The conjuction of two events $A$ and $B$ with given background $I$ is given by \n", "$$ \\boxed{p(A,B|I) = p(A|B,I)\\,p(B|I)}$$\n", " \n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- **All legitimate probabilistic relations can be derived from the sum and product rules!**" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Independent and Mutually Exclusive Events\n", "\n", "- It will be helpful to introduce some terms concerning special probabilistic relationships between events. \n", "\n", "- Two events $A$ and $B$ are said to be **independent** if the probability of one is not altered by information about the truth of the other, i.e., $p(A|B) = p(A)$\n", " - $\\Rightarrow$ If $A$ and $B$ are independent, given $I$, then the product rule simplifies to $$p(A,B|I) = p(A|I) p(B|I)$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Two events $A_1$ and $A_2$ are said to be **mutually exclusive** if they cannot be true simultanously, i.e., if $p(A_1,A_2)=0$.\n", " - $\\Rightarrow$ For mutually exclusive events, the sum rule simplifies to\n", " $$p(A_1+A_2) = p(A_1) + p(A_2)$$\n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- A set of events $A_1, A_2, \\ldots, A_N$ is said to be **collectively exhaustive** if one of the statements is necessarily true, i.e., $A_1+A_2+\\cdots +A_N=\\mathrm{TRUE}$, or equivalently \n", "$$p(A_1+A_2+\\cdots +A_N)=1$$\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note that, if $A_1, A_2, \\ldots, A_n$ are both **mutually exclusive** and **collectively exhausitive** (MECE) events, then\n", " $$\\sum_{n=1}^N p(A_n) = p(A_1 + \\ldots + A_N) = 1$$\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### The Sum Rule and Marginalization\n", "\n", "- We mentioned that every inference problem in PT can be evaluated through the sum and product rules. Next, we present two useful corollaries: (1) _Marginalization_ and (2) _Bayes rule_ " ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- If $X \\in \\mathcal{X}$ and $Y \\in \\mathcal{Y}$ are variables over finite domains, then it follows that \n", "$$\n", "\\sum_{Y\\in \\mathcal{Y}} p(X,Y) = p(X) \\,.\n", "$$\n", " - Proof:\n", " $$\\begin{align*}\n", " \\sum_{Y\\in \\mathcal{Y}} p(X,Y) &= \\sum_{Y\\in \\mathcal{Y}} p(Y|X) p(X) \\\\\n", " &= p(X) \\underbrace{\\sum_{Y\\in \\mathcal{Y}} p(Y|X)}_{=1} \\\\\n", " &= p(X)\n", " \\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Summing $Y$ out of a joint distribution $p(X,Y)$ is called **marginalization** and the result $p(X)$ is sometimes referred to as the **marginal probability**. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note that this is just a **generalized sum rule**. In fact, Bishop (p.14) (and some other authors as well) calls this the sum rule.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Of course, in the continuous domain, the (generalized) sum rule becomes\n", "$$p(X)=\\int p(X,Y) \\,\\mathrm{d}Y$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### The Product Rule and Bayes Rule\n", "\n", "- Consider two variables $D$ and $\\theta$; it follows from symmetry arguments that \n", "$p(D,\\theta)=p(\\theta,D)$, and hence that\n", "$$p(D|\\theta)p(\\theta)=p(\\theta|D)p(D)$$ \n", "or, equivalently,\n", "$$ p(\\theta|D) = \\frac{p(D|\\theta) }{p(D)}p(\\theta)\\,.\\qquad \\text{(Bayes rule)}$$ " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- This last formula is called **Bayes rule** (or Bayes theorem). While Bayes rule is always true, a particularly useful application occurs when $D$ refers to an observed data set and $\\theta$ is set of model parameters. In that case,\n", "\n", " - the **prior** probability $p(\\theta)$ represents our **state-of-knowledge** about proper values for $\\theta$, before seeing the data $D$.\n", " - the **posterior** probability $p(\\theta|D)$ represents our state-of-knowledge about $\\theta$ after we have seen the data." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "$\\Rightarrow$ Bayes rule tells us how to update our knowledge about model parameters when facing new data. Hence, \n", "\n", "
\n", "
\n", "Bayes rule is the fundamental rule for learning from data!\n", "
\n", "
" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Bayes Rule Nomenclature\n", "- Some nomenclature associated with Bayes rule:\n", "$$\n", "\\underbrace{p(\\theta | D)}_{\\text{posterior}} = \\frac{\\overbrace{p(D|\\theta)}^{\\text{likelihood}} \\times \\overbrace{p(\\theta)}^{\\text{prior}}}{\\underbrace{p(D)}_{\\text{evidence}}}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note that the evidence (a.k.a. _marginal likelihood_ ) can be computed from the numerator through marginalization since\n", "$$ p(D) = \\int p(D,\\theta) \\,\\mathrm{d}\\theta = \\int p(D|\\theta)\\,p(\\theta) \\,\\mathrm{d}\\theta$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "- Hence, having access to likelihood and prior is in principle sufficient to compute both the evidence and the posterior. To emphasize that point, Bayes rule is sometimes written as a transformation:\n", "\n", "$$ \\underbrace{\\underbrace{p(\\theta|D)}_{\\text{posterior}}\\cdot \\underbrace{p(D)}_{\\text{evidence}}}_{\\text{this is what we want to compute}} = \\underbrace{\\underbrace{p(D|\\theta)}_{\\text{likelihood}}\\cdot \\underbrace{p(\\theta)}_{\\text{prior}}}_{\\text{this is available}}$$ \n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- For a given data set $D$, the posterior probabilities of the parameters scale relatively against each other as\n", "\n", "$$\n", "p(\\theta|D) \\propto p(D|\\theta) p(\\theta)\n", "$$\n", "\n", "- $\\Rightarrow$ All that we can learn from the observed data is contained in the likelihood function $p(D|\\theta)$. This is called the **likelihood principle**." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### The Likelihood Function vs the Sampling Distribution\n", "\n", "- Consider a distribution $p(D|\\theta)$, where $D$ relates to variables that are observed (i.e., a \"data set\") and $\\theta$ are model parameters." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- In general, $p(D|\\theta)$ is just a function of the two variables $D$ and $\\theta$. We distinguish two interpretations of this function, depending on which variable is observed (or given by other means). " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- The **sampling distribution** (a.k.a. the **data-generating** distribution) $$p(D|\\theta=\\theta_0)$$ (which is a function of $D$ only) describes a probability distribution for data $D$, assuming that it is generated by the given model with parameters fixed at $\\theta = \\theta_0$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- In a machine learning context, often the data is observed, and $\\theta$ is the free variable. In that case, for given observations $D=D_0$, the **likelihood function** (which is a function only of the model parameters $\\theta$) is defined as $$\\mathrm{L}(\\theta) \\triangleq p(D=D_0|\\theta)$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note that $\\mathrm{L}(\\theta)$ is not a probability distribution for $\\theta$ since in general $\\sum_\\theta \\mathrm{L}(\\theta) \\neq 1$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Code Example: Sampling Distribution and Likelihood Function for the Coin Toss\n", "\n", "Consider the following simple model for the outcome (head or tail) $y \\in \\{0,1\\}$ of a biased coin toss with parameter $\\theta \\in [0,1]$:\n", "\n", "$$\\begin{align*}\n", "p(y|\\theta) &\\triangleq \\theta^y (1-\\theta)^{1-y}\\\\\n", "\\end{align*}$$\n", "\n", "We can plot both the sampling distribution $p(y|\\theta=0.5)$ and the likelihood function $L(\\theta) \\triangleq p(y=1|\\theta)$." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[32m\u001b[1m Activating\u001b[22m\u001b[39m project at `~/github/bertdv/BMLIP/lessons`\n" ] } ], "source": [ "using Pkg; Pkg.activate(\"../.\"); Pkg.instantiate();\n", "using IJulia; try IJulia.clear_output(); catch _ end" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd0AT5/sA8DeDvfeIgCxFAREVlKGgiAMXCOIeiButs1pHlX617tHaakVtqy3OOttqFUEBAUGmMlRU9pYdIEByud8f194vhmEEQhLyfP5KLm8uz12Se+597733peA4jgAAAABpRRV1AAAAAIAoSXAirK2tTU5Ojo2Nffv2LYvFEnU4HcrJyQkMDDx37hy55OHDh4GBgZGRkaIL6tOCg4MDAwNbWlqE9xGhoaGBgYFZWVnkkm3btgUFBQnvEwkxMTGBgYH3798X9gcBwW3dunXlypXk08LCwpUrV/7000/dXC2TyQwMDNy/fz+5JCUlJTAw8MaNG+SSX3/9deXKle/evevmZ3Xf1atXV65cmZGRIepAekZISEhgYGB+fr6oAxEALoGSkpLc3d0pFAq5FXQ6feTIkefPnxd1aO149uwZQmjevHnkkiNHjiCETp8+3ZthTJkyhdxdMjIympqaAwcOnDlz5nfffffhw4e25YcMGYIQYjKZAq6/oKAgJCQkOjpa8JBWrVqFEHr48CG5xNDQUE5OTvA1dO7Zs2chISE5OTl8y8+fP48Q+vbbb3vqg0Dnrly5IiMj4+fn10kZIyMjOp1OPk1JSUEI+fj4dPOjy8vLEUKjRo0ilxAp8MsvvySXLFq0CCH0WT9dIVm3bh1C6N69e52USUlJkemAqH7Sf/zxR0hISGtrK99yf39/hFBiYqJIovos9F7LuD3l2bNnHh4eLBZrxIgRkyZN0tLSKioqyszMjIiIuHfvXmBgoKgD/DRjY2N3d3cGg9H7H21ra6unp4cQwjCsuLj49u3bt27d2r59+759+zZu3Mh7buHo6KipqUmnC/oLyczMXLly5apVq0aPHi3gWwYOHOju7q6pqfm5WyGgmzdvHj169MaNG6amprzLDQ0N3d3dTUxMhPS5gA+GYWw2m81md1LG2dm5qqqqF4LR1dV1d3e3sLDohc8SBi6Xy2azFRQURo4cyfeSkZGRSELau3fvy5cv586dKyMjw7vc2tra3d1dVVVVJFF9FslLhNu2bWOxWFu3bj106BDv8oqKivT0dFFF9Vn8/f2Jc6Xet3PnztmzZ5NPKyoqTp06deDAgc2bNzOZzD179pAv8bblCsmGDRs2bNgg7E9pa/LkyZMnT+79zwWduHr1au980OjRo588edI7nyU8RkZG4r8Vu3fv3r17t6ijEIiEJUIMw+Li4hBCGzdu5HtJV1fXw8ODb2FGRkZGRkZJSYm8vLydnZ2TkxOV+tFl0dzc3OrqaisrKwUFhaioqJcvX6qoqEyaNMnQ0JAokJOT8+TJk/r6+uHDh48ZM4b3veXl5UVFRUZGRrq6usnJyfHx8RiGOTk5OTg4dL4VxBtNTEy0tbWJJW/evGloaLC1taXRaBEREVlZWQoKCuPGjbO0tGz79sLCwrCwsPr6enNz84kTJ9JotBcvXigpKVlZWX1i97Whq6v7zTff2NnZ+fn57d27d86cOQMHDiReevXqVVNTk729Pe8ey8jIePnyZWlpqaqqKoPBcHZ2VldXJ3YjcYnlw4cPycnJRGFizzQ1Nb169UpdXd3c3Ly8vDwsLKysrGzSpEm2trYFBQUfPnwYMGCAiooKX2CNjY33798vKCjQ19efPHkyX60xIyOjpaVl+PDhfG95/fq1pqYmUf978eJFWVkZQignJ4cMydraWl5evqqqKi8vj8Fg6Ovr866htrb20aNHhYWF8vLyDg4Obb/H5ORkOTk5GxubxsbGe/fuFRQUGBgYTJw4kfweQZelp6ez2exhw4Z1Xoz46g0MDMh/KJPJjIiIyMvLo9Fo9vb2Li4uvA0bbdXW1r5//15fX7/dJpnY2Njk5GQqlerq6jp06NC2BVpaWh4/fpydnU2hUAYNGuTu7s5XDSJUVFSEh4eXlJSoqqo6OTnZ2tq2LcPlch8/fpyRkaGoqDh27Nh2/+xdkJOTU1NTM2jQIEVFRd7PSk1NVVBQGDx4MLHkw4cPBQUFxL/gzZs3kZGRLBbLxsZm3LhxfAdJQn5+fnR0dHl5uZaWlqWlpZOTE41GI/7dRBeNtLQ04hNlZWWJ7SWOrnyREMujoqIqKir09PTc3Nz69+/P+2p9ff3bt291dHSMjY2LiooePnxYW1s7YMCAiRMnysrK9sguaoeo22Y/D4fDkZeXJ76Vzkvm5uYaGxvzbaydnd27d+94i82ZMwch9Pfffzs7O5PF5OXlb9++zeVyt2/fzvubmD9/PpfLJd9LXOo7fvy4n58f76fMmjWrubmZLCbINUKiOfHZs2e8BwIajXb48GG+7Tp+/Djvr8HMzIzodMN7FaRdxDXCq1evtvvqxIkTEUIbN24kl/BdI2xsbJwxYwbf/qTRaK9evcJxfN68eW1/Wt9//z2O46mpqQghHx+f06dPk5EfO3YM7/gaYWxsrK6uLrkedXX1O3fu8EZLHDIwDONdSOzn+fPnE09510DKzMzEO7hGeP78eb42nNGjR5eUlJAFuFwuQsjU1DQ6Opp35aqqquHh4Z3vfCkXGhqKEJoxY0YnZT55jZDL5RItFk5OTuRV7bNnzxKnYiQHBwfeg4Pg1wj//vvvCRMm8K5q7dq1fEFGRUXxHVUsLCyeP3/OV+zAgQPEYYo0ffr02tpa3jIVFRWjRo0iC1AolK1btwpyjTApKQkhNGDAgI4K+Pr6IoRSU1N5FzY0NCCEhgwZQi45ffo0QujgwYN810TGjBnD1zOgvr5+wYIFfNnRyMgI/+/fzadfv37EG9teI2xtbV21ahXvqqhUalBQEJvNJss8ePAAIbRq1arjx4/zXpoZPHhwcXFxJ3umOyQsEeI4TrSMe3p6Zmdnd1Ls5cuXzs7OZ86ciYmJeffuXVRU1MKFCxFCNjY2HA6HLEYkQhMTE2dn5z///DMpKWnv3r00Gk1dXX3v3r1aWlo//fRTUlLSjRs3+vXrhxC6cuUK+V4inxkaGvbv3//u3bsFBQVRUVGOjo4IoeXLl5PFBE+EpqamHh4ef/75Z3Jy8okTJxQUFKhUKu8P+o8//kAI6erqXrp0qaCgICUlxd/fnzg17mYiPHv2LEJo2LBh5BK+RPj1118jhLy8vKKiogoKCjIyMm7fvj1v3rzXr1/jOJ6enk50zJs6deqj/xQUFOD//VX69eunoKCwa9eusLCwx48fE8eOdhMhnU7X09NbunTpixcv8vLyTpw4IScnJysrm5aWRhYTJBFGRUXNmjULIbR7924ypIaGBry9REjsWDU1tTNnzmRnZ8fHxxNZ38bGhsViEWWIRKiurq6pqbl69eqIiIj4+Hiij6u+vj5ZDLTV/UTY3Nw8d+5chNDMmTObmpqIhSEhIQghBoNx/vz5ly9fJiQkBAUFUSiUwYMHk2UET4SmpqZDhw79448/UlNTf/75Zy0tLYQQ7xlYZmamoqIijUb7+uuvMzIy0tPTN23aRKFQ1NTUcnNzyWLff/89QsjAwCA0NJQ48ri5uSGE3N3dyV8shmGurq7EBiYlJeXn558/f15NTY34L/dmIjQ1NdXT0zt79mxSUtKff/5J1Be3bdtGFmOz2cTRadSoUXfv3n337t3z589DQkLGjBmD43h9ff2jR4/MzMwQQn/++SfxLyO7HbVNhCtWrEAIWVtb37t37/3793/99RfRjrV69WqyDJEI+/fvr6KicvTo0YSEhLCwMCKGWbNmdbJnukPyEmFERISCggJ5OrZw4cJz586VlpYK8l7i8lhYWBi5hEiE1tbWvF2eiL8clUpNSUkhF965cwch5OvrSy4h8hmVSk1PTycX1tTUaGhoUCgUMk8Lngjd3Nx4D+7BwcHEcZx4yuVyBwwYgBB69OgRWQbDMCL1djMRxsbGIoS0tLTIJXyJ0MXFBSFUVVXV0fr/+ecf4jyObzl5zhgSEsL3UruJECE0ceJE3mLHjx9HCE2fPp1cIkgixHF8y5YtCKEbN27wfS5fIuRwOETjzO3bt8kyHA6HOGH/8ccfiSVEIkQI7dy5k3dtRIM875cC+HQzEVZXVxO55IsvviC/9OrqamVlZTU1tby8PN71rFmzBiF05swZ4qngidDCwoL3bObixYsIoQULFpBLfHx8EEI7duxo+3FLliwhnjKZTDU1NQqFkpCQQJZhsVjm5ua8PzDiYGJvb897Un79+nXiByZIIpSTkxv+MbJT7mclQjk5ubdv35ILX79+TaVSzc3NySXEKfKoUaN4W7n4EMeK+vp6vuV8iTArK4tCoSgrK5eVlZFliouLFRUVKRQKcUqN/5cIKRRKZGQkWaympkZFRUVGRoa37tiDJO8+wnHjxqWnp8+dO1dNTe3du3e///778uXLDQ0Nvb29KyoqOn8vcZr//PlzvuXr1q3jbegnz+Ds7e35Fubm5vK9d8KECTY2NuRTdXX1wMBA/L/f+mfZvHkzb6OBp6cn7ydmZ2dnZ2fb2NiMHz+eLEOlUnukvwnRKlhfX99RAQ0NDYQQkWy6QFdXd+nSpQIW3rx5M+/TFStWqKio/PPPP0K6WzQ1NTUvL2/QoEG8bb80Gm3r1q0IoVu3bvGVJ/IrifiacnJyhBEbyMvLc3FxiY6OPnTo0Pfff0/+QW7fvt3Q0LBw4UK+3r9EZurCTaLr1q3jbc8kmknJf19LS8v9+/dlZWU3bdrE+66vvvqKQqHcuXOHOE+KiIioq6sbN24ccXpKkJeXJ/6k5G+JODhs2LCBRqORxXx9ffm6N3eCw+Hkfay4uPhzNxkh5OPjw9uBduDAgSYmJgUFBRiGEUsuX76MENqzZ4+cnFwX1s+LOA9YunQp0XGdYGhouGjRorYHTEdHR+KQS1BXV3d0dGSz2UVFRd0Mo10S1lmGYG5ufvnyZTabnZCQkJiY+ODBg/Dw8Lt3775//z4pKYn8wlJTU48cOZKWllZQUNDY2Ei+vbKykm+FRE2LpKOjgxDiu3atrq4uIyND9L/gZWdnx7eEuMaemZn5udtFdlQhED8X8hPfvHmDECKvdZOsra0/94PaIk4Y23ZaIS1duvTvv/+eNm2ai4uLp6enh4fHqFGjeP/GnRswYIDgd2Lw7VIlJaUBAwYkJydnZ2e33dvdR3xT9vb2fJ0siM44fN+jnp4e30Up4msiah6gZ719+9bZ2bmmpub69et8V+KJloa3b99+9dVXvMuJmzS6cBM3379PV1eXQqGQ/77379+3tLRYWloSTaYkIyMjHR2dioqK4uJiIyMjYnSItl1++H5Lr169Qm1+51QqdciQIW1Ptdtlbm5OHBC6iW+rEUJ6enq5ublVVVXEhfAXL14ghHirBF1G7Jy2q2r3j9ZuYAih8vJyvs41PUIiEyFBRkbG1dXV1dV148aNUVFREydOzMjIuHHjxvz58xFCERERXl5eCCEPD4/p06cTFZrXr19fuHCBPNkh8XVqIg6IfAsRQlQqFW8zNCuRNXkRPyAmk/m5W8T3icTJL/mJzc3N6L+aGa+2S7rg7du3CCEDA4OOCvj4+Ny7d+/QoUOxsbExMTF79uzR0dH55ptvVq9eLcj6Be9XSaVS2xYmdnIXdqkgiJOAtp1riD8e34e2+6tACJENp6AHlZWVVVVVMRiMtkfPmpoahFBsbGzbBh4NDQ3Bz7pInf/7OvqRIIT09PQqKiqI1hQBf0tEsY4OHb2JvMxE4ttwJpNJo9F6JDBiq3mrgwRx+KNJXtNou9zc3IgWfLKj/I4dO1pbW+/du3f//v2DBw9u27Zt27ZtfPc/9Ii27bFE5aDHbyMlEl5JSQnf8q41ifAhmpKIC/gdIXrKVFRU3L59e/ny5Uwmc82aNb///rsg6++8RzsvLpf74cMHvoXETiZ3abt/Cd5K/2ch6sFtq3REbUBNTa1rqwXdN3r06LNnz5aUlLi6uvINPEZ8a0ePHq1uD3EhrQcRH9fuxRfe30nnvyXyB6ysrNzu2nqkXYH4r/XUv0NNTQ3DsLYtYV0g4M4RiT6SCBFCRIso2Vj34sULHR0d3stpiCdN9qC2HYiJ6/y8Fw57xNChQ6lUakJCAt+lsu6PWfr8+XOiB0FAQMAnC2tqanp7e589e5ZIgeSYjcRFVg6H081gUJtdymQys7OzZWVlyRZs4v4/vn8U75ilvCG1bQDgQ9zzRHSM4l2emJiIeqjlGXTZsmXLLl26VFlZ6eHhQTTTEYg6ItHJqxeYm5vLy8vn5+fznaXl5eVVVlZqamoSdyUS//q2aZj4LZHHBOJHxfc7xzAsLS2t+6EK+O8QELGfiWNaRwT87xOb3/YgzLdzRELyEuGvv/7adqymkpISotcieZO1trZ2XV1dbW0tWSY/P/+XX37p8XjCw8N5f77V1dW//PILlUolaqg9SEdHZ9KkSRUVFUQvSkJpael3333X5XViGHb16lUvLy8OhxMQENB20CYS0azBi+jhSY7KTRwICgsLuxwMibjLkHwaEhLS0NAwdepUsi8D0aeAt0NEU1NT2/1ARPjJkOzs7CwsLN68ecPbL4bD4Rw+fBghRNyDAURozpw5t27dqqurGzt2bEJCArHQz89PVVX1+vXr7R6j2/5cu0lWVnb69Omtra28/z6E0IEDB3Ac9/X1Jeph48aN09TUjIyMjI+PJ8s0NTUR91SQv6WZM2cihL777jve5HH9+vWCgoLuh0rczED0vSTgOH7w4MGurY246yw4OLiTrmrEH+2T3VhmzpxJpVIvXLjAW78sLCwMDQ2lUqlEZ1dRkbxrhMuXL9++fbufn5+zs7OhoSGTyXzx4sWZM2cqKiqGDBlC/MIQQmPHjg0NDfX29t6zZw+DwUhMTNy5c6eBgUGP9+4zMjKaNm3akSNHhg8f/vbt2507d9bV1QUFBRE/x5514sSJZ8+e7dq1Kykpyc3Nrby8/MKFC3Z2dmFhYQKu4cqVK8R5KJPJLC0tffbsWVlZGYVCWbZs2alTpzp547BhwxwdHadNm2Zubq6kpPTq1SvizkJyrDhTU1MtLa2IiIi1a9daWVnJyso6Ozt34SyPTqdnZWXNnz9//fr1ysrKf/75Z3BwsLy8/N69e8kyc+bMuXDhwpYtW1paWuzt7XNycg4ePNj2shAxNMzhw4fr6uoYDAaFQpk1a1bbS6pUKvX48eMzZswICAgoKioaP358eXn50aNHExMT7e3tlyxZ8rmbANr18uXLtWvXtl1+7NixT/ZInDZt2u3bt319fSdOnHj//n1nZ2dNTc0ff/xx8eLF7u7umzdvdnZ21tfXLywszMjIuHjx4po1a3p8GpN9+/bdv3//8OHDGIbNnj0bw7CLFy+ePXtWU1OTHEhMUVFx//79q1atmjFjxoEDB0aOHFlYWPi///0vNzd34sSJRK8FhJCXl9e4ceMeP348ffr07du36+johIeHb9++3djYuPu50MfHZ/v27adPn1ZTU5swYUJZWdmZM2e6vNoFCxaEhoY+evTIxcXlq6++sra2rq6uTktLCw0NJU9KHBwc/vrrr4CAAH9/f3V1dSUlJaKjBh9LS8ugoKAffvjB3d19//79AwcOfPXq1VdffcVisTZs2CCMA+ZnEMY9GUIVGBjYtk8HjUabPXt2eXk5WayiomLEiBG8ZXx9fYl2vHXr1pHFiPsI4+PjeT+CqBnwDrNCkJOTMzAwIJ8StwOeOHGC71tfvHgx712Jgt9HWFhYyPtxRM4eP34878KXL1+OGTOGOP1UVlZev349McIqX7G2eGefIHfaoEGDVq1a1e7w8Hz3EY4dO5bvOp+8vPz//vc/3qF2Hj58yNv/lm9kmbYf0dHIMs+fPycH0EIIaWtrP3jwgO+9O3bs4L3VZMqUKcTZAO99hDiOHzhwgLebXycjy1y+fJmvQyBR/yYLkCPL8EVC3HAWHBzc4a6XesR9hB0h7j8TZPaJyMhIFRUVJSUl8q7NW7dute1DaGFhQY71053ZJygUipmZGe+SZ8+e8XUmt7GxefHiBd8bT548SVwFJFAolLlz5/IN11JVVTV27FiyDJVK3bNnT4+MLIPj+MWLF3lvBRk6dCjRJ7PdkWX43kuMscV7q19jY+Py5cv5uohbWVmRBRoaGubOnUuezXQysgyHw9m0aRPvOauMjMzWrVt576ckR5bhC2zBggUIobi4uE42vMsouATOUI/jeFZWVk5OTmlpKULIyMhoxIgRbbtgcbnc+Pj4N2/e0Gi0ESNGEONNlJWVqaqqkv0SKyoqGhoaGAwG7zlpY2NjeXm5mpoa35ExNzeXSqWS9y0dPXr0yy+//OGHH9auXZuRkZGcnIxh2MiRI/muKjU3NxcWFqqqqpLdpWpraz98+KCrq0t2xCgpKWlubjY2Nub9iXA4nIKCAgUFhbaJv66urr6+Xl9fX0ZG5p9//vHy8goMDCSO7x0pKytramoinxI3I3dyJl5QUNDS0mJhYUHmP2Ic0bKyMi6Xa2xsPHz48HY7rDY3NxNltLS01NTUWltbi4qKlJSU2vYW+/DhQ21tLYPBIHuI5efnc7lcU1PT5ubmiIiIoqIiHR0dT0/Pdu/rePv2bUxMDJfLtbOzGzFiREtLS3FxsbKyctsebnV1dcTMBv369ZOVla2vry8vL9fW1uaLv6mpKTo6Oi8vT0FBwcHBoe2dKjk5OXQ6nW+QrYaGhoqKCg0NjR7pvtsnMZnMTnpbEL+xvLw8DodD3tNGfJttfzbEH1ZOTo4cKZTD4SQmJr5586a1tdXAwMDS0pJ30F0Mw3JychQUFIiRoRBCjY2NJSUlmpqa5L+bWKehoSHfuGg5OTkyMjJ8UzpwOJxnz569efOGGMLG0dGx3ZuI6urqoqKiSkpKlJWVnZ2d263u4DiekJCQnp6uoKAwZswYY2PjyspK4n/dts8kqaWlpaCgQFZWtvPpU0pLSx8/ftzY2Ghpaenm5kahUHJzc2VlZcn9UF9fT1zd5LsdiDgWmZiY8G1XWVnZ06dPq6ur1dTUrKys2h2IlTjIkP+RsrIyJpNpbGzMd5wpKyuLioqqqanR1NR0c3Pj+4pZLBYxmjFf7/GOvqYeIZGJUEzwJkIRhuHj43Pnzp1Lly61O+AnAACAzkneNUIp5+HhsWTJEjs7OxUVlXfv3oWEhNy5c2fw4MGivdQMAACSCxKhhElMTHz8+DHvEmdn58uXL3d/ACQAAJBO0DTadc+fP3/y5Imnp+cnJ1HrQS0tLfHx8fn5+dXV1aqqqsOGDWu3sR4AAICAIBECAACQapJ3Qz0AAADQgyARAgAAkGrC6ixTVlaWlJT0+vVrNzc3YoCPtojhCahU6uLFi2FERwAAACIhrEQ4Z84cNpudl5dHo9HaTYRpaWljxozZuXMnm812cXFJSEhoOwEVAAAAIGzC7Szj5eXl6em5cePGti8tXrxYT0+PGNd47dq1OI53PtYlAAAAIAwiu0b49OlTDw8P4rGHh8fTp09FFQkAAABpJrIb6svKysjRQXV1ddvON0s6dOjQ06dPyfHoqFTqN998wzc+nvjgcrm8g0GLMwkKFSGE47jgs/uKFpVKFcZwiGLi+++/X7x4Md8fUPy/HfH/tUOE3dflCEWWCGVkZMi5uFpbWzsZGCU3N5dOp5MXGuXk5DQ0NGRlZXsjys/HZDLbHSFaDElQqBwOh81mS0p2aTtfZl/y888/jxs3jjcRcrlcFoulpKQkwqg+qbGxUcx/7RBh93U5QpElQgaDUVRU5OjoiBAqKioix5JvS11d3djYmJiyR/zRaLR2R6MXQxIUKo7jXC5XUqLFMEzUIQAAPkOv1nMrKysjIiKIxzNmzLh69Srx+Pr16zNmzOjNSAAAAACCsBJhcHDwiBEj4uLijh8/PmLECGKuxcTERD8/P6LA+vXrk5OTJ02aNH78+Pfv369evVpIkQAAAACdEFbT6JIlS6ZNm0Y+NTU1RQi5urqSvUP19fXT09MjIyOpVKq7u7ukXP4BQMzV1NS8fPlSR0en7cTChObm5vv37zc0NHh6evLO+fzmzZvY2FgGg+Hp6SnmfSIA6FnCSoT9+/fv378/30IVFRUbGxvyqaKiopeXl5ACAEAKbdiw4aefflJSUpo5c+b58+fbFmCxWC4uLurq6iYmJhs3boyIiCBmL7l9+/by5cv9/PwSEhLOnTt348aNXo8dAJGB8z4A+o6tW7cymczAwMCOCly5coVGoz169OjXX38NCgrau3cvQgjH8Z07d/74449nzpyJjo5++vRpfHx8L0YNgIhBIgSg7zA0NOz8zqJ79+55e3sT/W99fX3v37+P4/j79+/fvn1LdFhTUVGZMGHCvXv3eiliAMQAzFAPgBQpLi4mL94zGIzm5uaqqqqSkhJNTU0FBQVyeScDXNTX158/f54c4MLU1HTmzJlsNlvM756ECLtP/CN8WIRP6M9W+Dit0Wi0T17zhkQIgBThHXqDqBdyOBwMw3jHhaFSqeRgF21xOJy6ujoZGRniqby8PPYfYQbeXRBh94l5hLEVlLXJci8YmOzHYxwJ0vMLEiEAUkRfX7+8vJx4XFZWRqfTdXV1a2trq6urORwOnU5HCJWXlxsaGna0Bk1Nzc2bN9va2pJLuFwujuNi3vFb/Ecmggi7o64VLXvG+cGhVU9VsQtvh2uEAPRxXC63oqKCy+UihMaOHRsWFkYsDwsLc3Nzo1KpFhYWenp6jx8/RghxOJyIiIixY8eKMmIAPlNQHOZlRJlk2MUKK9QIAeg7wsLCbt68GR8f39LSsnLlSi8vrxkzZhQVFZmYmBQXFxsaGi5ZsuTEiRPLli0zNTU9evTozZs3EUJ0On3Hjh2BgYHr16+PiooyMDAYP368qDcFAEHdyOWmVOJJ3nSM1cU1QCIEoO/Q09MbPnz48OHDiadEC6eWltbPP/+soaGBENLQ0EhKSvr999+ZTObjx4/t7e2JkqtXr7a0tIyMjPTy8lq8eDHcUA8kRVEjHhSH/Vi+PaUAACAASURBVD2BrkhHzK6uBBIhAH2HnZ2dnZ0d30IlJaWlS5eST3V1dTdv3tz2vePHj4eKIJAsXBwtisQ22dAcdLo1Cxic9wEAAJBIR15yOTjaMqS7iQxqhAAAACRPahV+LB1L9KbTuj0nNNQIAQAASJhmDC2KxL5zopkodzsNQiIEAAAgcbYkYHZalHnmPZPCoGkUAACAJHlYhP9VgKf59Fj+gkQIAABAYnxoRkujsUtjaRpyPbZOaBoFAAAgGXCEAqI4AQMo7gY9cGmQBDVCAAAAkuGnLG4pC90aRuvZ1UIiBAAAIAFe1+J7UrCYaXTZnm7KhKZRAAAA4o7NRYujsH0jaAPVerJRlACJEAAAgLjbnYxpyaMVVkLJWdA0CgAAQKzFlOEX33JTfWR6vjKIEIIaIQAAAHFW24oWRGLnR9P1FIT1EZAIAQAAiK+gWGyqMcXLSEi1QYSgaRQAAIDYuvSOm1qFJ3kLN1VBIgQAACCO8hvwTQnYg0l0RSFnKmgaBQAAIHa4OFoShW2xpdlrCbFRlACJEAAAgNg5+IKL4WiTbW8kKWgaBQAAIF6SKvHvM7Gknph0VxBQIwQAACBGmjhowRPs+1E0I6VeSYOQCAEAAIiVTfHYSF3KnB6adFcQ0DQKAABAXPxTiIcV42kzezU3QSIEAAAgFj40o2VPsctjaaoyvfq50DQKAABA9HCElkZzlg6kuPXopLuCgEQIAABA9H7K4paz0G77Hp50VxDQNAoAAEDEXtXiwSnY02l0GVHUzqBGCAAAQJSISXf3Owhl0l1BQCIEAAAgSruSsH5KlGUDRZaPoGkUAACAyDwpxS+9x9N8RJmMoEYIAABANGpbUUAUFuJK05YXZRiQCAEAAIjG6lhshgllijAn3RUENI0CAAAQgd/fcV8Kf9JdQYg+AgAAANKmsBH/MgF7MImuIAZZCJpGAQAA9CoMR/OeYFuG0IYKf9JdQUAiBAAA0KsOvuDKUtEmG3FJQGJQKQUAACA1kivxk5lY4gw6VSxqgwhBjRAAAECvaeSgeU+wH51pxspikwYhEQIAAOg1G+MxJ13KLFPxSj3QNAoAAKA33M3nRhTjqb076a4gxC4gAAAAfU8ZC62KwW6Op/fypLuCEK/6KQAAgL4HRyggirPCiuqsJ0aXBkmQCAEAAAjXD5ncyma0SxST7goCmkYBAAAIUWYNvi8Vi5sumkl3BSGucQEAAJB8LRia/wQ74ECzUBXHRlECJEIAAADCsjMJM1WhBIpu0l1BCLFp9Nq1a1FRUQYGBkFBQZqamm0L/PPPP2FhYQghT09PLy8v4UUCAACg9z0uwa/l4KkinXRXEMLK0kePHt21a9eIESPevHnj5ubG4XD4CoSEhCxbtmzQoEGDBg1avnx5SEiIkCIBAADQ+6pb0JIo7OcxIp50VxBCSdRsNvvYsWNXrlxxd3dfsmSJlZXVvXv3ZsyYwVvmzz//XL9+/YoVKxBC9fX1d+/eXblypTCCAQAA0PvWxGK+ppQJDPG9NEgSSo3w/fv31dXVo0ePRghRqdSxY8fGxMTwlRk6dOjz5885HA6Hw0lISLC3txdGJAAAAHrfxbfcrFr8gIOY3i/BRyg1wrKyMg0NDRrt312go6OTl5fHVyY4OHj69Ol6enoIoREjRgQHB3e0tpycnKysrOTkZOKpnJzckSNHNDQ0hBF597FYLHLDxZwEhcrhcNhsNo7jog5EIBiGycrKijoKAEQml4l/mYCFe9HlJeMAI5xEKCsry3tRsLW1VV6ev5F4165dra2tsbGxCKG1a9fu3Lnz8OHD7a5NS0vL1tbW29ubeEqlUrW0tGRkxG+UHoQQQq2trXJycqKOQiASFCqNRqNSqZISLZvNFnUIAIgMF0cB0dj2obQhmhLQKEoQSiJkMBg1NTUNDQ3KysoIoaKiIgsLC74yFy5cuHTpkpWVFUJox44dc+fO7SgRqqmpGRkZzZ49Wxih9jgajSYp1SwJChXHcS6XKynRYhgm6hAAEJl9aVwZKlpvLdb3S/ARSqwmJiZDhgy5du0aQqi6ujosLIyoz1VWVj548IAoo6en9/r1a+Lxq1evdHV1hREJAACAXpNQgZ/Owi6MoYnPpLuCENbtHYcPH547d+7Dhw9TU1NnzJgxbNgwhFBaWpqvr29jYyNCaP/+/YsWLYqMjEQIPX78+LfffhNSJAAAAHpBIwctisJOOtEYShKVBoWXCD09PTMyMhISErZt2zZ8+HBioZOTU2JiIvF46tSp2dnZqampCKEzZ85oa2sLKRIAAAC9YP0zzFWP4m8mSY2iBCHe8K+vr89376CSktLgwYPJp9ra2p6ensILAAAAQO+4nceNLJWAQWTaJZFBAwAAEB/lLLQ2jntjPE1FTLvzf4Lk1WEBAACIDxyhpdGcFVZUJ10JuzRIgkQIAACg605mcGta0M6hEpxNoGkUAABAF2XV4vvSsLhpdLoE50GoEQIAAOiSFgzNe4wdcaRZqklqoygBEiEAAICu2J6IWahRlgyQ+DwCTaMAAAA+W3gx/kcunjazLyQRic/kAAAAelltKwp8ip11pWlJxkj4nwCJEAAAwOdZ9hTz7U+ZbCTZlwZJfaFWCwAAoNf8ms19XYuHuved9NF3tgQAgBB69OjRjz/+yGKx5s+fv3jxYr5Xnz59Ghoayrvkm2++0dfXv3HjxqNHj8iFP/zwA8wtDNqVw8S3PZekSXcFAYkQgL4jIyPD19f3zJkzOjo6S5YsUVJS8vPz4y3AYDDGjx9PPI6Ojr558+apU6cQQvHx8ZWVlXPmzCFekpSpH0Ev43DRgifYTomadFcQkAgB6Dt++umnBQsWzJs3DyH09ddfnzx5ki8RmpmZmZmZEY8vX768aNEiOv3fg8CgQYNmzZrVywEDybIvDVOWQV/Y9LXOJX1tewCQZqmpqU5OTsTjUaNGpaSkdFSyvLz8/v37S5YsIZc8fPhw9uzZX375ZU5OjrDjBJIotYZ69jX3ghutT1UGEUJQIwSgLykvL9fQ0CAea2lpNTY2NjQ0KCsrty3522+/jRw50srKing6YsSIAQMGaGpqPn782M7OLiUlxdLSst2PKC4u9vf3V1BQIJ5aWlqeO3euubkZx3EhbFCPaWxspFDE+gAu5hE2cCiLY+nHhrFVuS0NDaKOpgPt7kN5eXmy2aMjkAgB6DuUlZVZLBbxuLGxkU6nkxmLz6+//rp161byKXl10M/Pr7y8/OzZs0eOHGn3jbq6usHBwRYWFsRTDQ0NZWVlGo2mpKTUY5shBDiOt3tCID7EPMIvojE3PfZcK0VRB9KZLu9DSIQA9B0mJiZkw2ZOTk6/fv3a7fYSFxdXWFjId/mQZGpqWlVV1dFHyMjIDB482NbWllzC5XK7FzUQd7fzuNFl+FNPjqgDERa4RghA3zFnzpzffvutqakJx/GQkBCynnf+/PnXr1+TxX7++ec5c+bwnjtnZWURD/Lz869du+bq6tqbYQNxVtyIr47FLrnTlOli3frdHZAIAeg7/P39hwwZYm5uPmDAgA8fPpCNn/v37yc7zjQ0NPzxxx+BgYG8b/T19TUwMLCxsRk8ePDMmTN5O9EAaYYjFBCNrRlMGymxk+4KAppGAeg76HT6lStXSkpKWCyWubk5uTw7O5tsI1VWVq6vr+d746tXr4qKihoaGvr37y8vL997EQPx9l0Gt5Ej2ZPuCgISIQB9jaGhId+ST/aaQwj169dPOOEASZVejR9Iw+Jn0PvgDRMf6+N5HgAAQBc0Y2h+JHZkJM1Mpa+nQUiEAAAA2voqEbNUpSy2lIocAU2jAAAAPhJWjN/qK5PuCkJathMAAIAgKptRYDR20Y2m2Scm3RWEVFR7AQAACCgoDptjRhln2PcvDZKgRggAAOBfP7/hvqnFf3OTrtQgXVsLAACgI+/q8e2J2JMpdDkpm48SmkYBAAAgDhctjMS+tqdZa0hRoygBEiEAAAC0NxVTl0VrraUxKUDTKAAASLu4cvzsa26Kj4zUVQYRQlAjBAAAKdfARkuisR+caQZiPdugEEEiBAAAqbY2DhtrQPEzld50AE2jAAAgvf7I5T6rwFN8pDoXSPXGAwCANCtuxNfFYXc96UrSnQqkty4MAADSjIujxVHYOus+PumuICARAgCANDqRwW3hoq/sIAtA0ygAAEifzBr84Avs2fS+P+muIOBcAAAApAuLg2Y/xo6PolmoQhpECBIhAABIm22JmI0GZaEFHP//BU2jAAAgRcKK8Tt5UjTpriBgXwAAgLSobEZLo7HfpGnSXUFA1RgAAKTFsqfYfHPpmnRXEFAjBAAAqXD2NbegAb/uAYd9frBHAACg73tfj+9Mwh570WWhHbAN2CUAANDHcbhofiQWPIxmqwmNou2ARAgAAH1ccAqmLovWDIYDfvugaRQAAPqy2HL85zfc1JlSOumuIOAEAQAA+qx6NloYiZ0dTdNXEHUoYgwSIQAA9FlrY7EJDMo0YzjUd0agplEWi6WgAKcTAAAgSW7mcp9V4KnSPemuINrfQWVlZRcvXgwPD3/37h2TycRxHCGkra09bNgwPz+/yZMnKyoq9m6cAAAAPkNhI74mDrs3ka4sI+pQxB5/ImxpaTl06FBycvKUKVO+++47S0tLWVlZ4qXGxsa0tLSYmJipU6cGBAQsWLCAQoGLrwAAIHa4OFoShW2woY3QhqP0p33UcFxZWblhwwYfH5+7d++uWLHC2tqazIIIISUlJRcXl23btkVERKioqGzZsqW1tbXXAwYAAPAJR9O5rVy0dQhcGhTIRzXCqKioEydOyMvLd/4eCoXi7e3t6OgYFRXl6enZbhkul3vhwoX4+HgjI6MvvvhCTU2tbZn6+vqQkJC3b9/q6OgEBARYWFh0eTMAAAAQUqvwY+nY8xkw6a6gPjpf8PX1/WQWJBkaGnaUBRFCu3btOnny5OjRo9PT0ydNmkRcZeRVXV3t4OCQmJjo4OCgqKj49u3bzw0dAAAAn2YMLY7CToyimShDGhTUp3sTNTc3p6en19fXDx06VEtLS5CVNjQ0nDp1Kjo62s7Obu7cucbGxtHR0W5ubrxlvv32W2tr6+vXr3cxcAAAAG1sScAGqlHmmUOj6Gfg31kxMTG8lbPo6GgrKytHR8fx48czGIxjx44JstL09HQZGRk7OzuEEJ1OHz16dExMDF+ZsLAwX1/fc+fO7d27Nz4+vntbAUCfUlNTk5WV9eTJk6SkpJKSEg6HI+qIgGR4WIT/mY+HuNJEHYiE4a8RGhoaTps27eXLlzQaraKiYt++fSdOnDAzM6NSqcXFxT/88IOJiYmfn1/nKy0rK9PW1iaf6ujolJaW8pXJy8sLDg6ePXu2iorKtGnTTp48OXfu3HbX9v79+1evXiUnJ/8bMZ1+7NgxTU3Nz9vQ3sJisWg0yfgVSlCoHA6HzWa3bWAXTxiG8fYyE9yrV69OnTp1//59ZWVlDQ0NdXX1xsbGmpqaysrK0aNH+/v7T5s2Dbpqg44Qk+6GjoVJdz8bfyI0NTXNycmJiIiYMGFCTEzMpUuXdHR0iJdsbW09PT2/+uqrTyZCOTk5NptNPm1tbW3bWUZGRmbatGn79u1DCGlpaR08eLCjRKijo0On0318fIinNBpNR0dHRkZMb41hs9mCX2cVLQkKlcPh0Gg0SYmW98cvoObm5v3795eUlCxYsODw4cN89+lyudzMzMw7d+789ttve/fuHTRoUM8FC/qOwKfYIkvKWAM4Vfps/ImQy+VyOJy8vDyEkJKSEl+NgUaj6evrf3KlDAajvLyczWYT6aqwsHDw4MF8Zfr160d2Ex0wYEBJSUlHa1NVVe3Xr5+/v/+nt0YMUKlUKlUyWuclLlRJifZzK2319fW7d+/+4osvzMzM2i1ApVJtbW1tbW1ZLNb+/fs9PDzc3d17IFDQh5x9zS1qxP+ASXe7hP/IkpiYyOFw1NXVa2pq7O3tz58/z+VyEULx8fEVFRWNjY2CNE/Z2trq6+vfuXMHIVRSUvL06dPp06cTjx8+fEiU8fPzi46OJh5HRkYOGTKkB7cKAAmSn59/8ODBjrIgLwUFhb1798rJQcsX+Mj7enx3MnbJnQaT7nYN/+mDlpZWSEhIbW3tH3/8gRBSV1evqqrS0dHx9vYePnz4qFGjNm3a9MmVUqnUY8eOBQYGXr9+PSEhISgoiPiTx8XFrV+/vri4GCG0bt06d3f3sWPHqqmpJSUl/f3330LYOgAkgK2t7WeVd3JyElIkQBKRk+5aqUOjaBfxJ0JLS0tLS8u25SIjI2VlZQU5aSXMmDHDwcEhJSVlz549NjY2xMKJEydGRUURjzU0NJKSkuLi4jAMc3R0VFZW7uomACDxOBzOpUuXHjx4UFlZqaysPGjQoNGjR3t4eHSt0w2QKntSMF15yspBUBnsOkEblK2srD531YaGhoaGhrxLVFRUVFRUyKcyMjJ8NxcCIJ1OnDghIyPj7e1dUVHx6tWre/fuHThwQF1dfcWKFVu2bCE7rAHAJ6YM/wUm3e02uLIKgOgpKyuvXr2afJqRkREVFaWtrX3z5k07O7uzZ89OnTpVhOEB8VTXihZGYedG02HS3W76jNp0dXX1kCFD7t27J7xoAJBOBgYGu3fv5r1xXkNDY/bs2devX3/+/Pkvv/wC/zvQVlAc5mVEmWoMtcHu+oxEKC8vb2FhoaqqKrxoAJBO3t7eurq6Q4cOPX36dFVVFe9L/fr1u3XrVkZGhqhiA+LpRi43pRI/4igZY2KIuc9oGlVUVLx165bwQgFAmq1du9bFxWXLli1ffPHFwIEDra2tlZSUbG1ttbW1GxoaampqRB0gECNFjXhQHPb3BLoiXN3qCe3sRRzHb926dfPmzdTU1KqqKhaLpaGh0a9fv4kTJ65YscLAwKD3owRAGtjb20dERKSkpNy4cSM8PHzRokX19fWysrIuLi5HjhwRdXRAXHBxtDAS22hDc9CBRtGe0U4iDAgIaGpqsrCwsLOzU1FRoVKpdXV1dXV1JSUlEyZMuHz58ufe9gQAENywYcOGDRu2f/9+hBCLxVJQgI4Q4CNH07kIwaS7PYk/EUZFRfn7+3t5ebVbmslkHjp0CBIhAL0DsiDgQ066S4XaYM/hT4T5+fmLFi3qqLSKioq5ubmQQwJAurx48SIhIUHw8pqamp8c+B70SU0cNP8JdtIJJt3tYfyJ0MDA4P79+x3VCOvr63Nzc4UfFQBShMFgmJubCz7JFO8cZ0CqbEnAhmtTZptBo2gP40+E48ePX7BgwZkzZ2xtbVVVVYnpYBoaGurq6iorKxMSEq5evSqKOAHos7S1tT08PEQdBRB3D4rwe4V4mg/0E+15/PuUQqGEhobevXv35s2bycnJ1dXVra2tWlpahoaGU6dO/fbbb/X09EQSKAAASK0PzSgwGrs0lqYBU48IQTsnFxQKxdvb29vbu/ejAQC0Kzk5efjw4aKOAogGjlBgNBYwgOIOk+4KB7Q1AyABiGnRgHT6KYtb3ITvHgaDyAgLNDcDIHosFuv+/fvEJNhtcTic169f93JIQEy8rsX3pGAx0+gw6a7wQCIEQPSoVOr+/ftra2vbfZXFYgk+FSjoS9hctDgK2zeCNlANGkWFCBIhAKInJyf3xRdfODg4DB48uN0Cq1at6uWQgDjYnYxpyaMVVlAZFC7YvwCIBWtr65cvX3b0qomJSW8GA8RBdBn+21v8whg6VAaF7ROJ8Pnz58SDkpKSoqIi4ccDgJQaNmyYj49PR69u3769N4MBIlfXihZHYedH03RhlD3h+0QizMzMJB58+PChoqJC+PEAIKWoVKqcHNwjBv61JhabZkyZbAS1wd4A1wgBAEC8XH7PTavCk7zh+NxLYEcDAIAYyWPiG+Oxh5PoCnB47i3QWQaAvqakpKSuru5z38XlcgsKCpqamoQREhAQF0cB0diXQ2hDtaBRtPdAIuwxSUlJqzdv91u6Zs+3hyorK0UdTt9x49btxWs2zl2+7qez51tbW0Udjlj78OHDqFGjnJycTE1NN23a1LbAtWvXKDwiIyOJ5VlZWQMHDhw3bhyDwfjhhx96NWjA4+ALLoajjTZwZO5VsLt7xhdbd01as+cM1zXMZvO3Odo2rhMehUeIOiiJx2KxnDynLj//+LKm312zFVvCy6xGuEpD7+Wmpqb3798Tj9PT0wV/4zfffGNqapqfn5+dnX39+vWIiHZ+hJ6envh/3N3diYXr16+fN2/eu3fv4uLiduzYkZ+f3+2NAJ8tpRL/LgMLdafRoDbYuyAR9oDo6OjQmFdVy28j6wmIYY05LSpf+ffitVtaWlpEHZpkCz5wNFXPo3bmcWThjIztmyZsy5t8aP7K9aKOS+hqamrIRJiUlCT4Gy9durR27VqEkLa29pw5cy5fvtxusebmZt6nZWVljx8/DgoKQggNGjTI3d392rVrXQwddBULo8yPxE460Yxh0t1eB4mwB/xy9XaN80pE4fn5qui0mLt+1rTjoK0/7t5rcQnkXYJbur7KLYAG0nbV1NTU1tZaWloSTy0tLfPy8toWi4yM1NXVVVNTW716NYvFQgjl5eWpqanp6up2/kYCm83OyspK/k9hYWGPb4h02p5Kd9ShzDGHY7IIQLekHlBRVYP0dPkWtijqVFVViSSePqOVw0Ey/LcTU5U0mEymlpaWSEISZ0wmEyGkoPDvHlNSUmrbZWb06NElJSXa2tp5eXne3t7BwcGHDh1iMpnku4g3lpeXd/QpFRUVwcHB8vLyxFMLC4vz5883Nzd3NGK4mGhoaBB1CJ35p4QaXkp/NpnFZOKijqVDYr4PUQcRysvLy8jIdP7GTyRCCk8th/cx4GU3yOJhcTqXYc27ULEsfcAAX1GF1DfoaGkW1xQjDcb/L+JiqL5cQ0NDdEGJL6JKV1tbq6KighCqqalpO422oaEh8aB///5bt249dOjQoUOHdHV1ecf7bveNJAaDcenSJVtbW3IJl8ul0+lKSko9uC3CQOwWMVTBQhuTOedHthhqKIs6lk8Q231I6lqEn6iGk9PzWllZDRo0qAsfIA3WLg/QijqOaorJJbT0eyZ0prW1dSfvAp+0Z8s69ZsbEJv173Ocq3Tvm0X+vlQqNB+1Q15efuDAgeSwiAkJCXZ2dp2Ur6mpIbKXubk5lUrNyMgQ8I2gB+EIBT7lLB1IcdUV6yp13/aJGqG6ujrxAAZ/6gSDwbj96+mFq2fXmbnXum/QCg0YwlC7+keoqOOSeN7Tp5VWVH5z1K3BYwtbs7/arU1zpnoeCN4p6rjEV1BQ0Ndff21kZPT+/fu//vorNTUVIVRWVjZhwoSIiAgdHZ3Tp08zGAwjI6OXL18GBwfv378fIaSsrLxkyZINGzYcPXo0IiKisLBw1qxZot4UaXE6i1vOQrvtac2Nog5FisE1wp7h4uz0NiXur8Tsta/k426eNTY2FnVEfcTqZQGL5vp/86Q4gykbuvEfTU1NUUck1tauXdvS0rJx40ZVVdU7d+4QsxjKyMhYW1sTl0kUFRXPnj1bUVHBYDBOnz7t5+dHvPHw4cPBwcGrVq1iMBjh4eGKioqi3Ayp8aoWD07BYqbRZaio+dPFgbBAIuwxNBrNyspKqYBtbAxXsHqSkpKSmZlZVQWmqQnNEp9AoVC2bNmyZcsW3oVaWlpXrlwhHi9ZsmTJkiVt36igoHDo0KFeiBCQiEl3DzjApLuiB9daAABABHYmYQxFyrKBcBAWvXZqhB8+fEhOTmaxWNra2nZ2dqqqqr0fFgCAAL21+6SnZXjoO26qzye69YPewZ8I//e///3111+ysrLFxcUFBQUUCsXBwWHu3LlLly4V/46zAPQNOjo6w4cPJx5PmjRJtMGAHlfbihZGYudG0/Vg0l3x8FGtnMlkmpqaJiYmxsbG5uXlVVVVrVq1atiwYSEhIf379z979iyOi+/NngD0GbKysuSIAfr6+qINBvS4NbHYDBPKFJh0V2x8lAiVlJRSU1MbG//txquhoWFiYrJ///6srKywsLDw8PDly5dDLgQAgC77/R33RRV+0IEm6kDA//soEVKp1CVLlri7u9+5c4dvwKThw4dfv3596tSpHQ3jCwAAoHOFjfiXCdilsTSYdFes8HdYGjJkyOXLl48cOWJmZrZp06YXL17wjt7m7e1N1hcBAAAIDsPRvCfYFph0V/y003PX0tIyJibm1KlTmZmZN27cMDExsbGxmTFjxuLFi6dNm2ZkZNT7UQIgPd6+fRsYGHjhwoVOpoAAkujgCy4VoU0w6a74ab9+TqFQpkyZMmXKlIaGhri4uKKiopqaGkNDw/Hjx+vo6PRyiABIFQMDAzMzs0uXLgUFBeno6Li7u48dO9bNza1///6iDg10XXIlfjITS5xBp0JtUPx8oqFaWVl5woQJvRMKAAAhpKysvHPnzp07d7a2tj5//jwyMvL48eNLly51dnY+duyYo6OjqAMEn62Rg+Y9wX50hkl3xRRU0gEQU7Kysq6urrt27UpMTDx+/PimTZs2b94cHR0t6rjAZ9sYjznrUWaZwvFWTH1UIzxy5EhiYqLgb549e7avL0y5B0BPamxs3L9/v6Kior+/PzHdvKysrKqqqo+Pz5QpUzZu3DhmzBhRxwg+w918bkQxnjoT+omKr4++my+//FJUcQAACDdv3tTT03v+/LmDg4OZmZmrq6u+vn5WVlZAQICsrOzMmTPz8/NNTExEHSYQSBkLrYrBbo6nq8JgamIMquoAiBcbG5uhQ4eGhoaWlZUFBwerqKg0NDQcOnQoKSlp4MCB+fn5MNihpMARWhbNWWFFddaDS4NirbPa+okTJ27cuOHm5ubu7u7i4kJMZg0AEKphw4bFxsYeO3YsICBg+vTp06dPJ5ZraWnNnz/fyMgIJmWUFD9mcstZaJc9DCIj7jpLhIGBgQwG3V+NAQAAIABJREFUIyoqauPGjW/fvnVwcHB3dx8/frybmxuVClVJAITFxcXFxcWFb6G8vPzu3btFEg/ogswafG8qFjedLgMHS7HX2Vekqqrq7+9P3FlfWFi4fv36R48eBQQEDBo0KCUlpddCBKBvi4mJefnypYCFm5qafv75Z6HGA7qvBUPzn2AHHGgWqtAoKgEEPVfR09Pz9/ePiIjYv3//b7/99uWXX757906okQEgJVxdXf/5558ff/yxubm585LPnj1bu3bt1KlTeycw0GU7kzBTFUogTLorITr7nu7fv+/o6Lh8+fKIiAgMwxBCKioqTCZz5MiRd+/evXLlSm8FCUAft23bNltbWy8vr02bNt25cycvL48Y45fNZpeXl8fFxR05cmT8+PHh4eFnzpzR09MTdbygM49L8Gs5+LnRcGlQYnR2jfD9+/dHjx599OjRmjVrqqqq7OzsEEILFy5kMpk4jhsbG/dWkAD0fW5ubuHh4fHx8Tdv3jxy5EhJSUl5ebmampq+vr6ZmZmHh8f169ehm4z4q25BS6Kwn8fQtOVFHQoQWGeJ0MrKisvl7t27d+/evZmZmZmZmcbGxqNGjfr6668fPHiwaNGiXosSAGlApVKdnZ2dnZ1FHQjoujWxmK8pZQIDLg1Kks4SoaenZ3Jy8vnz55ctW2ZtbW1tbU0s37Jli6Oj4/jx43slQgCkUXV1dVNTk76+Pp0OI5JIjItvuVm1+AU3+MokzEfXCNvOPj98+PBly5bxLVRTU5s2bZqCgkIns9W3trZ+9913AQEBBw4c6GQKQxzHT506de3atc+PHIC+icPhLFiwQEtLy8jIyMDAYOfOnSwWS9RBgU/LYeJfJmCX3GnycHFQ0nyUCENCQtLT0wV8599///3PP/909OqqVatu3bo1ceLEuLi4WbNmdVTsl19+2bNnz/nz5wX8UAD6vFu3bi1btqy6ujopKWnfvn3JycmjR4+uqakRdVygM1wcLY3Gtg+l2WpCo6jk+SgRrly58u7du99++211dXUn78nNzd28eXN1dbWXl1e7BUpLSy9fvnz16tU5c+Zcu3bt6dOn7ebXkpKS48ePr1+/vjsbAEAfw+VyHR0dNTQ0hg8fvnLlygcPHnz//ff79+8XdVygM9+mcWWoaANMuiuZPmrLplAou3btev78+YIFC+h0uouLy8CBAzU0NJSVlevq6qqqqtLT0+Pi4jQ1NQ8dOmRqatrRShMTE83MzAwNDRFCioqKjo6OcXFxtra2fMWCgoK+/fbbgoICYWwYABJq+vTpv/zyy5o1a8jxm1xcXARvqgG9L6ECP5WFJXvToTIoodq5qOvo6Hj//v2ioqI7d+5EREQUFhYymUwdHR0DAwNra+sVK1b069ev85WWlZVpaWmRT3V0dEpLS/nKhIaGysjIeHt7nzx5svO1vXv37vXr18nJycRTGo12/Phx8exH3txMwXFqU1OTqAMRCIvFotEk42pGayvicPCmJkzUgQgEwzBZWdkuv72pqen69eshISEzZsxwd3d3dnam0+nEbYUIITabLSMDExmIkUYOWhSFnXSiMZQgD0qqDns39evXb+3atV1bqby8PJvNJp+2tLTIy390T01lZeU333wTGRkpyNp0dXVlZGR8fHyIp3Q6XVdXVzy70sk24xQKxrexYovNZktKqHQ6h07H5eUlIwHw/vi74M6dO7t27aqpqXn69On69evfvXunpKQ0f/788PBwJyenI0eOBAcH91CkoAesf4a56lH8zaBRVIJ9Ip2UlZXFxMTIyck5OTlpa2sLuFIGg1FYWIjjOIVCQQgVFhbyzd8bGRlZWlpKzC9aV1fX1NRkbW2dmZnZ7tpUVVX79evn7+8v4KeLEJWKI4RJyojkVCpVgkKVoB1L/Oy7bPTo0X///fe4ceNmz56NECovL4+MjIyMjFy3bt27d+/U1NQgEYqP23ncyFI8zUccz8uB4Dr7/qKiogIDA3Ecz8vLo9PpK1euPHDggCCTMbm6urLZ7CdPnowbNy4rK+vVq1eTJ09GCOXk5OTm5np4eEyePJkcZfjixYsRERG///57j2wPAJJu4MCB5ubmT548qamp0dDQ0NPTmz17NpEUS0tLg4KCRB0g+FdJE746FrvtSVeWjKYK0KHOTrHT0tKys7Pfv3/f0NDw119/NTY2Tp8+vbW19ZMrlZOTO3bs2OzZs2fOnOnh4bF3717ikuHDhw83b96MEFJSUjL7j5aWloKCQiddbwCQNnQ63dPTU0NDg2+5gYEBzMQkJnCEAqKw1YNoTrpwaVDidVYjVFBQIBp5FBQUJkyYMGHChKioqF9//XXlypWfXO/ChQvd3d0zMzOPHDlibm5OLJw3b96UKVP4Si5evLiTGw0BALyGDh0q6hAAQgidzOAy2WjnUMlorged6+xb9PHxOXPmDO+Vfzc3N8F7qRgZGU2aNInMggghNTW1tkN1q6mpGRgYCBwwAACIWFYtvi8Nu+hGo0Me7BM6+xr/+uuvjRs3MhgMf3//n376KSsr6+HDh/r6+sSrcP8fAEAKtWBo7mPssCPNUg0aRfuIzqp39fX1NTU12dnZkZGRDx8+3LlzJ5vNnjRpUn5+vqura2ho6OHDh3stUAAAEAc7kjBLNUrAAKgM9h2dfZeurq4nT54sLCxcv379nTt3Kisro6OjnZ2dHzx44ObmlpOT02tRAgCAOAgvxq/n4CGukjESBRBQZzXCESNGjBgxory8nHhKpVLt7e3t7e03btzI4XC2bt3aKxECAIBYqGpBS6OxX8fQtOREHQroUZ+u3evp6bVdSKfTd+zYIYR4AABATK2MwfxMKeNh0t0+p+sDIgg+0AwAAEi6X7O5r2vxUHcYRKYPgi8VAAA+IYeJb3uOhXvRYdLdPgk6PgEAQGc4XLTgCbZzKG0ITLrbR0EiBACAzuxLw5Rl0Bcw6W7fBU2jAADQoaRK/OxrbhJMutunwTkOAAC0r4GN5jzGTjnTDBUhD/ZlkAgBAKB9XzzDxhpQfPrDcbKPg6ZRAABox+08bnQZngqT7koB+I4BAIBfcSO+Oha760lXgUl3pQBU+QEA4CM4QstjsKDBtJEw6a50gEQIAAAfOZHOrWtFO2DSXakBTaMAAPD/XlbjB19gCTPoNKgNSg045QEAgH+1YGhhJHZ0JM1UBdKgFIFECAAA/9qWiFmoUhZZwoFRukDTKAAAIIRQWDF+KxdPmwlHRakDJz4AAIBqWtDyp9ivbjRNmHRX+kAiBAAAtDIG8zeleBjCpUFpBI0AAABp9/MbbnYd/jtMuiut4IsHAEi1HCa+PRGL8KLLwaS70gqaRgEA0ovDRfOfYF/b02xh0l0pBokQACC99qZi6rJorTUcCaUaNI0CAKRUXDl+9jU3xUcGKoNSDs6DAADSqIGNlkRjPzjTDBRFHQoQNUiEAABptDYOG2tA8TOFYyCAplEAgPT5I5f7rAJPgUl3AUIIEiEAQNoUN+Lr4rC7nnQlOP4BhBA0jQIApAoXR4ujsLUw6S7gAYkQACBFTmRwW7hoO0y6C3hA0wAAfUpFRcXVq1dZLNbMmTMtLS35XsUwLDY2NikpiUKhjBs3zs7Ojlj+4sWL7OxsstjMmTNptD44zkpmDX7wBfZsOky6Cz4Cp0UA9B0VFRX29vZpaWnV1dUODg6pqal8BX788cegoKCioqKCggI3N7fz588Ty3///fcDBw6E/wfDsF6PXehYHDT7MXZ8FM1CFdIg+AjUCAHoO86ePTtixIhffvkFISQrK3vo0KGrV6/yFli8ePH69euJx7a2tocPH162bBnx1MvLa9++fb0ccG/alojZaFAWWsDZP+AHvwkA+o7w8PApU6YQj728vB49esRXQF1dnXwsKysrIyNDPs3Ozj537lxYWBiHw+mFUHvZ43LanTz8tEsfbO8F3Qc1QgD6jpKSEn19feKxvr5+dXU1i8VSUFBoW5LJZP7vf//bunUr8VRNTQ0hlJSUdOLECXl5+aioKBUVlXY/orq6+tixYzo6OsTTfv36rVixoqWlhU4X34NJZQtldQL94miuEmppaRF1NB1oaWmRlZUVdRSdkdAI6XT6Jy94i+9vFwDwuWg0GpfLJR5zuVwKhUKlttPq09zcPHPmTFdX18DAQGLJ119/TTzgcDhOTk4//PDDjh07OvoINTU1DQ0N4qmGhgb1Pz28MT1n9TM0y6jV3UAWIfG9Oijm+xBJbIQUyqe/dEiEAPQdBgYGpaWlxOPS0lItLS05OTm+Mq2trb6+vtra2ufOnWt7jKDT6ePGjXvz5k1HH6GmprZs2TJbW1tyCZfL5XA4vK2sYuXsa25RE/fXUVyxjZAgIyMDEXZTlyMU6/QOAPgskyZNunv3LvH47t27kydPJh6/fv26trYWIcRms/39/eXl5X///Xfe9iKymyiGYVFRUQMGDOjdwIXlfT2+6//au/O4Jo7+ceBDEu5DCAmnAUEFRRSVAPogl6CAVgSCKFgtXlBt7aFtPav25VOVPh48HrVaj0fr+bS1PhUVUBBpqyiIVg5RURTlvsMRINns74/5Pvvkx61lswE+7782yyT5ZJid2Z2dnckkTnixNaCqA12DK0IABo6lS5cePnxYJBKZmJj8+9//TktLw/vxiNDIyMh//vOfly5dCgoKioyMRAix2eyzZ88ihCZMmDB+/Hgul5uSkoIQ+vDDDxn8FX1FJkfzU4nNE9ljuWoNDUxHA1QYNIQADBxcLjcrK+vixYsSieTLL7+0sLDA+48fP44v8kJDQ8ePH0+lp7pGz5w5k5mZ2djYGBgY6Ovrq8ojX3rvq/uEoQZa4QAXg6AHA6G4AwAoBgYGCxcubLfTy8sLb9ja2tra2nZ8l6Ojo6OjI+3BKdEf5eTRx/L7sOgu6AU4VwIADDRiKVqYShyawjbt5MkRANqDhhAAMNCsvEX4D1WbZQX1G+gV6BoFAAwo55/L71aS94KhcgO9BWUFADBwvGoiP7pNXPHn6EDdBnoNug4AAAOEnERRN4lPHdnOPBgiA94ANIQAgAFiZ7a8TY4+HwfVGngz0H0AABgI7leTOx8Sd2fDorvgjcGpEwCg32sh0Hs3ibjJ7GH60AyCNwYNIQCg3/vsDmE/RC1yOFRo4G1A1ygAoH+7+oqMLyIfhEJtBt4SjUXn/v37hw8fbmtrmz9//tSpU9v9taqq6ueff87MzGSz2dOmTQsNDe3NqlEAAKCosgUt/Y045cM2VOklY4FKo6sn4enTp97e3iNGjPjb3/4mEolSU1PbJTh58mRKSoqLi4uTk9PHH3/89ddf0xQJAGAAW/obsXCkmo85nEaDt0fXFeGBAwfmzJmzevVqhFBVVdWuXbu8vb0VE3zyySfUUsKGhoabN2/euHEjTcEAAAakQ/ny103kj77QKQr+ErquCNPT06mWz9vb+/bt2+2/mPW/r66oqODz+TRFAgAYkArE5KZ7xGlvWHQX/FV0nUmVlZUZGxvjbR6PV1NT09raqqmp2TFlUVHR1q1bT58+3dVHPX36ND8/PzMzE79ksVi7du2iPlylNDeryeXsxsZGpgPplaampv5yX7alBclkao2NUqYD6RW5XK6hATes6CWTo3dTia8mskcZ9o8yDFQZXQ2hlpZWW1sb3m5tbeVwOOrq6h2TlZeXT58+fc2aNf7+/l19lJmZmYaGRkhICPXJ5ubmqrlwqFYbyWIROjo6TAfSKwTRb0LV0JBxOKSOTidFSAVJpf2jwe7XNt0j+FooZjRcDII+QFdzYmlpWVRUhLeLioosLCwU+0KxqqqqadOmRUZGfv755918lL6+vqWlZXh4OE2h9iEWi0SI6PhLVROLxepHofajjO0v19n91+9l5PEn8vuhsOgu6Bt01Swikejs2bMEQSCETp06JRKJ8P5Lly69fv0aIVRXVxcQEODv779p0yaaYgAADDz1bWjBTeKwB9sMFt0FfYSuhjAqKorFYgmFQk9Pz4yMjM8++wzvX7FiBR44s3v37qysrAsXLgwfPnz48OFjxoyhKRIAwEDywS1ihgAW3QV9ia6uUR0dnbS0tIyMDKlU6urqSo0duHXrFh7n8sknn0RFRVHpoTcJANCjnwrl96pg0V3Qx2gsTywWy83Nrd1OgUCAN7hcLpfLpe/bAQADzOsm8oNbxKXpsOgu6GPQvQAA6AfkJFqQSqxyZLvyofcI9DFoCAEA/cA/HsoJEn0Gi+4CGkAXAwBA1d2vJnfnwKK7gC5wegUAUGnNMhR5g9g7mW2tB80goAU0hAAAlbb6DiHkqc21hcoK0AW6RgEAqivhNXnlFfkgBGoqQCMoXgAAFVXZgpakEad92EadTNcPQJ+B3gYAgCoiEVqSRkTZqXnDoruAZnBFCABQRQfz5MXN5E9+UEcB2kEhAwConPw6cksW8dssDiy6C5QAShkAQLW0ydH8VGKrkG0/BDpFgTJAQwgAUC2b7hGm2ih6FNROQEmgaxQAoEJulpI/PCUfhHLgYhAoDZxzAQBURV0beu8mccSDzddiOhQwmEBDCABQFSv+IGZZqQUK4GoQKBV0jQIAVMLpAvmf1WQmLLoLlA7KHACAeS8byVV3iMQAjjbUSUDpoGsUAMAwOYmibhKfjWWPN4ZOUcAAaAgBAAzb8aecINGqsVAdAWZANwQAgElZVWRcDpEZDIvuAsbAKRgAgDHNMjQ/ldg7mW0Fi+4C5kBDCABgzOo7hAtPbd5wqIgAk6BrFADAjPgiMuk1eT8UaiHAMCiCAAAGVLagmN+JMz5sA3WmQwGDHvRIAACUjUQo6qZsib2aFyy6C1QANIQAAGU7kCevakGbJrCZDgQAhKBrFACgZI/qyK+yiD9mcThwHg5UA5REAIDySOXovZvEdhe2HSy6C1QGNIQAAOXZkElY6KgttYeaB6gQ6BoFACjJb2XkqQL5/RAYJwpUC5yXAQCUobYVvZtKHPPkmGozHQoA/z9oCAEAyrD8D2K2tVrAULg1CFQOdI0CAGj3Q4E8uwYW3QUqCsolAIBehQ3k6nTi2gxYdBeoKOgaBQDQiCDRwpvEGie2Exc6RYGKgoYQAECj7Q/kbDX0qSNUNUB1QVcFAIAu96rIfXlExmwOC64GgQqD0zQAAC2aZCjiBnHgb7DoLlB10BACAGjxaTrhbqoWZgOVDFB10DUKAOh7/3kpTy6GRXdB/wDFFADQxyokaPkfxLmpHFh0F/QL0GsBAOhLJEKL02TL7FmeZnBrEPQP0BACAPrS/lx5uQRthEV3Qf8BXaMAgD6TW0tuvU/cCuKowzk26D+gtAIA+kYrgd5NJba7sEcYQKco6E+gIQRg4CAIYuPGjQ4ODpMmTfr11187TfPLL79MmjRpzJgxmzZtIggC72xoaIiJibG3t586deqdO3fe7ts3ZBLD9NSWwKK7oL+BrlEABo49e/ZcuXLl4sWLBQUFERERd+/etbe3V0yQl5cXFRV17tw5W1vbuXPnGhkZffrppwihVatWlZeXJyQkJCcnz5w5s7CwUF9f/42+OqWEPP+cvB8CVQrof+DcDYCB4+DBg5s3b7azs5sxY0ZISMjRo0fbJThy5IhIJAoMDLS3t9+0adN3332HEBKLxadPn/7mm29sbGyWLl1qZ2d3/vz5N/reeqna4jTi0BQ2T6vPfgsASgMNIQADhEQief78ubOzM34pFAqzs7PbpcnJyVFM8PTp09bW1ufPn7NYrFGjRuH9zs7OOTk5b/TVH91lhwxTmyGAW4OgX4J+DAAGiMrKSoTQkCFD8EtDQ0O8p10aQ0NDKgFJkpWVlZWVldS7EEJGRkYFBQVdfUtRUZG7uzub/X9PR0ycOHHd9xcf1akfdG1uaOjDX9PHGhsbmQ6hBxDhX9dphFpaWurqPczsAA0hAAMEbuGamprw7b2GhgYjI6OOaajKoqGhAe8xNDRsamqi0ojFYi6X29W3DB069NChQw4ODvilpqamppbWbzwJ1+DN7ikq35ve9VQ+iPCve7sIoWsUgAHCwMDA2Nj48ePH+OXjx49tbGzapbG1tVVMYGpqqqenZ21t3djYWFpa2s0bKSwWy8DAwOi/dHR01BDShKfnQX8GDSEAA8fChQvj4uIIgiguLj537tyCBQsQQnV1datXrxaLxQihd99999y5cyUlJQRBxMXF4QQmJib+/v67d+9GCGVlZf3+++8RERHM/hAAlKkfNIRVVVU1NTVMR9ErUqm0ubmZ6Sh6Ky0tTbFDTJXV1tb2lzKAELpy5QpTX71p0yaJRMLn88eMGbNixQoPDw+EUFNT08mTJyUSCULIx8dn6dKlDg4OfD4fP3SI37hv374bN27w+fzp06fv37/fwsKi919aWVl59+5dOn5OX5FIJKmpqUxH0R2SJBMTE5mOogcJCQkkSTIdRXdu3LjR0tLyNu8k6SGTyQ4ePBgREfHFF1+UlZV1miY1NXXRokVLliy5detWNx81adIkf39/esLsY5fvFWjEPmY6it5ydXXtPudVx/xDN+z/nsx0FL3S1tamrq7ObAyNjY1tbW3dJGhra2tqauq4XywWy2Sy7j987NixDx8+VNzz448/BgcHv0WcSpOZmTl+/Himo+hOTU0NHrukyoYMGVJbW8t0FN0ZN25cVlbWW7yRrivCLVu2fPfddyKRqLa21sfHRyaTtUuQnp4eFBQ0efLkiRMnBgQE/PnnnzRFAsBgo6ur2/0wOXV1dR0dnY779fX1qeGgvUeq9lUCAD2iZdSoRCI5cOBAYmKii4tLaGionZ3d5cuXZ8+erZhmz549H3300bJlyxBCz54927t3b8eHfwEAAAC60XJF+OTJk7a2NqFQiBBSU1Pz9PRMT09vl+bOnTteXl5429PT862nNwQAAAD+ClquCMvKyrhcrpra/00zwePxqJHZlPLycmNjY7zN5/M7JqC8evXq3r17VGI1NTUHB4ceH5BkRC3LQD56ga/vcqYD6ZXHjx9/9NFHBgYGTAfSs8c6o8R6Q319v2Y6kJ6RJCmXy5mOgkbl5eUxMTG6urqKe8rLy6dNm8ZgVN1raGh49uyZKkcok8mamppUOUKEUHNzc0hICIejuk+fFxYWLl++vN2jhCEhIStWrOj+jbT8JG1t7dbWVuplS0uL4mGDaWlpUWlaWlo6vWOB7d69OzMzk8fj4ZcsFmvcuHEsliqOdyVJ8uXLl8OC1jEdSK+8evXK3NxclYs1pbGxUSKR8Pn9I2OXLl3KdAg0OnTokI6OjuIB2NLSUltba25uzmBU3SMI4vXr19bW1kwH0p3CwsJuHt9UBeHh4Soe4cuXLwUCQbvWoTcx01IJWlpa1tTUNDY26unpIYSKiopcXFzapRk6dGhRUdGkSZNwgqFDh3b1aeHh4eHh4XTECQB4U8HBwUyHAEAfo+W6avjw4Y6OjmfOnEEIlZSUJCcni0QihFBpaelPP/2E04hEoh9++AH3I506dQonAAAAAJRMjaahz6mpqXPmzHF2ds7Ozn733XdjY2MRQlevXo2MjKytrUUI1dTU+Pj4aGpqymQydXX169evq/4sdgAAAAYeuhpChFBdXd2DBw+srKxsbW3xnpaWlurqaktLS/ySIIjMzEwWi+Xs7Kya9/wAAAAMeDQ2hAAAAIDqU+kRg62trYmJifX19X5+fh3HpJEkmZycTL0UCAT29vbKDK+goOD3338XCARTp06lnhVRVF9fj+cPDAgIYPYpBalUmpiYWFNT4+vrS12RK0pOTqZOiSwtLUePHq3cAP9HJpPl5OTgCYm6SpORkZGTkzN27Fj8rCqDKioq8vLybGxsOh2RWFhY+OzZM+qlh4eHpqamEqPrMy9evLh586aZmZmfn1+nU880NjYmJCRIpdKAgADFtZ9ycnIyMjJGjhw5ZcoUWiMUi8UJCQmo62MtPz///v37mpqanp6e1BD03Nxc6sEtNpvdTZHrqwhJkgwICFBc/RFrbGxUfNh61KhR1PjBsrIyfOcoICCA1vJDEMT169fLysq8vLyGDRvW7q/19fUZGRmKe8aOHWtqalpSUpKXl0ftdHV1pa+uk0gkDx8+ZLFYHUdfYiRJpqSkvHr1asqUKSNGjKD295yHfTPFGw0kEomrq6uHh0dUVJSxsfG9e/faJZBKpQghb29vPz8/Pz+//fv3KzO8ixcvGhsbL1u2zMnJKTw8vGOCkpISgUAQEhISEhJibW1dWlqqzPAUtbW1ubu7u7u7L1q0iMvlpqend0zDZrM9PT1xTu7Zs0f5QWJpaWlaWlo8Hk9XV7erNFu3bhUIBDExMQKBYPv27coMr52ZM2dqa2vr6urGxsZ2muCrr74SCAR+/1VdXa3kCPvEtWvXuFzukiVLhEJhYGCgXC5vl6C6unrkyJEzZsyYO3euubl5YWEh3v/999+bmpq+//77I0eOXLlyJX0RlpSUWFlZ4WPNysqq47GGy8y8efOCgoIMDQ3T0tLw/gULFtjb2+P/zqxZs+iLsLS01NraOjg4ODQ01MrKqqSkpF2C7OxsTU1Nqqj85z//wfsfPHhgbGy8cOFCb29voVDY3NxMU4RyuTwwMNDZ2XnJkiXGxsbXrl1rlyA3N5cKD5/W4JmKjx8/zufzqT/l5+fTFOGhQ4c0NDS4XO7kyZO7ShMeHj5u3Lhly5bxeLyLFy/inb3JQ9VtCE+ePDl+/HipVEqS5JYtW2bPnt0uAW4I6+vrmYiOdHBwOH36NEmS9fX1fD7/zp077RKsXbt27ty5eHvu3Lnr169Xdoj/df78+TFjxuBZmHfs2BEQENAxDZvNrqioUHpo7dXV1VVUVGRmZnbVEFZXV+vo6OTl5ZEkmZubq6urW1dXp9wY/+fFixdSqXTWrFndNIQffvihkqPqc5MmTTp48CBJks3NzcOGDUtKSmqXYNu2bYGBgXh76dKlH3zwAUmSra2tpqamKSkpJEmWlpa7qj9LAAANQklEQVTq6OhQDWSfW7duHXWszZs3b926de0SFBYW4pqEJMn169f7+vri7QULFsTFxdEUlaL169fPmTMHb0dERKxZs6ZdguzsbAsLi45vFIlEGzZsIElSJpMJhcJjx47RFOG1a9cEAgGeiv3w4cNubm7dJD5x4oSdnR0+JTp+/HhQUBBNUSkqLy8Xi8XHjh3rqiG8c+cOj8fDLcLZs2cdHBzw/rCwMCoPXVxcOs1D1R2iEh8fHxwcjB/3DgsLu3LlSqcTdty6devGjRt4JKrSPHv27MmTJyEhIQghAwOD6dOnx8fHt0tz6dKlsLAwvC0SiTomUJr4+PigoCA8F09YWFhSUlJbW1vHZOnp6SkpKcyudjRkyBA+n99NguTkZFtbW9xz6+DgYGVllZKSoqzo2rO2tu5xOoLKysqrV69mZ2eT/fNmfGVlZXp6On66SVtbe8aMGR1Lcnx8PFXUw8LCcIKMjAyZTObt7Y0QMjMzmzx58uXLl2kK8tKlS9TzV50ea8OGDaP+U+bm5orTfRQVFV29erWgoICm2DDFLOqqNpDJZNevX7916xa1lBtJkpcvX8ZvZLPZISEh9FUj8fHxM2bMwBObhIWF3blzp6KioqvER48eXbx4MXU/CPf6ZmVlEQRBU3gIIRMTk+6fLIiPj58+fTrumJ09e/aTJ0+ePXtGkiSV+d3koeo2hMXFxdTdLEtLS6lU2vEfY2Zmtn///o0bN9rY2Fy4cEGZsXG5XG1tbSq84uLijmkU4++YQGnaRSKXyztOaGdiYnLo0KHNmzfb2NicO3dO6TH2VnFxseLcC8xmbI/YbHZhYeHBgwf9/f2nTp3aX1Z/VFRSUqKhoUGdnfSmqJeUlMjlcryTqivxfpqC7P2xJhaL4+LiqKl/1NXVMzIyDhw44OLismDBAvrmxutNhHp6env37l2xYsXIkSPx3Ms1NTUtLS3KqUYUIzQyMtLW1u7qu54/f3779m28pDNWV1f37bffhoWFCYXCsrIymiLskWLloK2tzeVyi4uLe5mHTA6WuX379vLlnUzLeenSJYFAQBAE9UwFvj/fbi0nDofz+vVr/KfTp08vXrx45syZyhmMQBCE4ugYNpvdcZ0pxTSdJlAaxZzEGx2DefXqFc7Jn376KSoq6p133sGzAqmadjnP4XAYzNgerV27dsOGDQghiUTi6em5c+fOzZs3Mx3Um1EsPKjroq54qOIes94cI30YZG+OtdbW1vDwcHd394ULF+I9hw8fxsW+srJywoQJ58+fj4iIYCTC0aNHU+Oqvvzyy5iYmAcPHuALLOVUI+3+0d0cWUeOHAkMDKSWbl6wYEFUVBRCSCaThYSEbNy48ciRIzQF2b1OK4de5iGTV4Rjx4491RlTU1OEkLm5OXUJWF5ezmKxzMzM2n0CNYBt7ty5eF5d5URubm5eW1tLZWh5eXnHQa3m5uaVlZVUgjda8rtvKeYk3ugYDJWTIpFIJpM9ffpUmRH2nuJvQUxnbI+oXNXW1g4ODr5//z6z8bwFMzOzlpaWhoYG/LKroq54qJqamrLZ7I7/KfomI+3NsSaVSsPDw/X19Y8ePapYLeINPp/v6+tL3z+oxwgVx+JGRERkZ2fLZDIej8fhcJRTjSj+v5qbmxsaGjr9LplMduLEicWLF3eMnMPhzJkz58GDBzRF2CPFnyCTyaqrqy0sLHqZh0w2hHp6eo6d0dDQQAh5e3snJSXhlElJSVOmTMG9/PX19S0tLe0+Cud+pw8G0GHEiBEmJiY3btxACMlkspSUFDzwWiqVUvfYfHx88LMTOH58s4QR7XLSzc0Nd+qKxeKOOZmTk9PW1tbN1K+MoEL18PDIyckpLy9HCJWWlubl5dE9Lv9NKZYBRVlZWVZWVsqP5y8yNze3t7fH5Ucul1+/fh0XdVzR4DTe3t4di7qzs3NjY2N2djZCqLm5+Y8//qDv4QQfHx/FEk4dazU1NXhIHUEQCxcuZLFYZ86c6fS2LkEQf/75p0AgoCnCTrNIMUJFWVlZeDZ8Fovl6enZ6U+jI8Lr16/jzuHExEQ7OzvcYDQ0NEgkEipZQkICQRAzZszo9EPu3btHXx52hcpDb2/vlJQUfH2SmprK5/NHjBjBYrG8vLx6zkNaxvf0hdraWisrqyVLluzYscPIyAg/gkOSpFAo3Lt3L0mSZ8+ejYyM3LFjx9q1a/l8vpKHZe7fv18gEOzatSsoKMjNzY0gCJIkk5KSqOGOjx49GjJkyLp169atW2doaEjfqOIeNTQ02Nravvfee7GxscbGxr/++ive7+7u/s0335Ak+fPPP8+bN2/79u3r1q0zNTVdvXo1U6HW1NRER0eHhoZyOJzo6Gjqfzp58uSdO3fi7cWLF7u4uMTFxQmFwpiYGKZCJUny5MmT0dHR1tbWbm5u0dHRt2/fJkkyKSlJT08PJwgKClq7dm1sbGxQUBCPx6Nv2CStTp48aWZmtnPnzvDwcEdHx9bWVpIk8U0sXOxfvHjB5XJXr169efNmAwODrKws/MaNGzeOHj16z549Pj4+tD6ckJ+fj4+19evXKx5rOjo6+DGATZs2cTicqKio6Ojo6OjoL774giRJuVzu5ub25Zdfbt++3cPDw97eXiwW0xehoaHh2rVrcYR42DNJknp6engU7vbt25ctW/bNN9+sXLlSX1//X//6F05w7do1Q0PDbdu2RUdHDx06lL4ncNra2hwdHcPCwnbu3Glubn7ixAm8f9asWYqjcIODg3HuURYtWrRq1ap//OMfkZGR+vr6GRkZNEWYnZ0dHR3t5eVlamoaHR397bff4v1UHhIE4ebmNmvWrN27d1tZWVEP1F2/fr3HPGRv2bJFWS33m9HS0oqMjCwuLm5oaPj73//u6emJ95uYmAiFQhMTEz6f39raWl5erqen9/nnny9atEiZ4bm6ujo4ODx69Gj8+PF79uzB9ya1tLRGjRo1fvx4hBCPxxOJRPn5+erq6nv37rWzs1NmeIo0NDTmz59fUlIiFou3bNni6+uL9/P5fGdnZ1NTUx6PJ5VKy8vLdXV1V61aFR0dzVSocrm8qqpqxIgRM2fOtLCwsLKycnR0xKEKhULcZ/7OO+/o6uo+f/48ODh4zZo1nU5loBz19fUaGhoeHh4TJkywsLBwcHAwNjZWLAO4r0YsFguFwkOHDqnyQkXdcHJyEgqFubm59vb2+/btw0uqaWhojBw5Ek9oYGhoOHfu3KdPn5IkuXv37nHjxuE3Tp061cLC4smTJ15eXl9//XWnT+L3CR6PFxYWlp+fz+FwFI81S0vLSZMm6evrq6uru7q6CgQCCwsLCwuLoUOHjhs3Tk1Njc/nV1RUtLa24geRO64W17cRPn78mM1m7927l5r6w9LS0s3NzcDAwNzcvKGhAXcgb9u2bfr06TiBra2tn59fbm6umZnZgQMHuh9T/Vew2ez58+dXV1dXV1d/8cUXs2fPxvuNjY0nTpxIdSdKJJKIiAjFR+YtLS0rKyurq6sdHBwOHjxI36wmra2tTU1NY8eO9fHxsbCwsLW1HT58OFLIQzU1tYiICJyNK1eupG732traTps2LScnp5s8hCnWAAAADGqq+/gEAAAAoATQEAIAABjUoCEEAAAwqEFDCAAAYFCDhhAAAMCgBg0hAACohI4TXADlgIYQAAAY9uLFi08++eTHH3/Ek9MCJVPpFeoBAGDAq6ioCA0NTUxM5PP5/v7+BQUFiqurAyWAK0IAAGDS8uXLY2Ji8Iwnmpqad+/eZTqiQQeuCAEAgDH37t1LTk4+deoUfllSUtIf163s7+CKcHD57rvveDze69evFXcuWrTIycmJvlVJAQBdiYuLe+edd/CCMFKpNC8vj74JRUFXoCEcXObNmyeRSI4ePUrtqa6uPnv2rEgkUlyWEwCgBARBXL16VSwWx8bGxsbGbtiwQSKRODk5MR3XoAN13+BiaGg4Z86cI0eO4IWbEULHjx8nCGLJkiXMBgbAIJSdnV1dXb1161a8PpSmpqatra2NjQ3TcQ060BAOOsuXL3/9+nVCQgJ+eezYsZkzZyptTWMAAKWoqMjU1NTJycnIyMjIyCgtLY1aPAgoEzSEg46bm5uzs/Phw4cRQjdu3Hj06FFMTAzTQQEwGNXV1VFPSrx48SIzM3PFihXMhjQ4QUM4GL3//vuXL18uKir6/vvvraysqFVAAQDKZG5uPmTIELy9Z8+eNWvWUEvgAmWChnAwioyMNDAw2LVr1y+//LJs2TL6lg4HAHTD3d1dIpEghG7dulVUVLR+/XqmIxqkYIX6Qerjjz/et28fm81++fIlnIQCwJSUlJRHjx5pampGRUVxOPBgNzOgIRykcnNzHR0dQ0JCLly4wHQsAADAJOgaHaRKS0sRQtHR0UwHAgAADIMrwsGIIAg/P7+qqqqHDx+qqakxHQ4AADAJuqQHnfnz56ekpFRVVSUkJEArCAAAcEU46Fy4cKGhocHd3R2WegEAAAQNIQAAgEEOBssAAAAY1KAhBAAAMKhBQwgAAGBQ+39eX4ol5mfRtwAAAABJRU5ErkJggg==", "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/html": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "using Plots\n", "using LaTeXStrings\n", "\n", "f(y,θ) = θ.^y .* (1 .- θ).^(1 .- y) # p(y|θ)\n", "\n", "θ = 0.5\n", "p1 = plot([0,1], f([0,1], θ), \n", " line=:stem, marker=:circle, xrange=(-0.5, 1.5), yrange=(0,1), title=\"Sampling Distribution\", xlabel=\"y\", ylabel=L\"p(y|θ=%$θ)\", label=\"\")\n", "\n", "_θ = 0:0.01:1\n", "y=1\n", "p2 = plot(_θ, f(y, _θ), \n", " ylabel=L\"p(y=%$y | θ)\", xlabel=L\"θ\", title=\"Likelihood Function\", label=\"\")\n", "\n", "plot(p1, p2)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "The (discrete) sampling distribution is a valid probability distribution. \n", "However, the likelihood function $L(\\theta)$ clearly isn't, since $\\int_0^1 L(\\theta) \\mathrm{d}\\theta \\neq 1$. \n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Probabilistic Inference\n", "\n", "- **Probabilistic inference** refers to computing\n", "$$\n", "p(\\,\\text{whatever-we-want-to-know}\\, | \\,\\text{whatever-we-already-know}\\,)\n", "$$\n", " - For example: \n", " $$\\begin{align*}\n", " p(\\,\\text{Mr.S.-killed-Mrs.S.} \\;&|\\; \\text{he-has-her-blood-on-his-shirt}\\,) \\\\\n", " p(\\,\\text{transmitted-codeword} \\;&|\\;\\text{received-codeword}\\,) \n", " \\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- This can be accomplished by repeated application of sum and product rules." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- In particular, consider a joint distribution $p(X,Y,Z)$. Assume we are interested in $p(X|Z)$:\n", "$$\\begin{align*}\n", "p(X|Z) \\stackrel{p}{=} \\frac{p(X,Z)}{p(Z)} \\stackrel{s}{=} \\frac{\\sum_Y p(X,Y,Z)}{\\sum_{X,Y} p(X,Y,Z)} \\,,\n", "\\end{align*}$$\n", "where the 's' and 'p' above the equality sign indicate whether the sum or product rule was used. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- In the rest of this course, we'll encounter many long probabilistic derivations. For each manipulation, you should be able to associate an 's' (for sum rule), a 'p' (for product or Bayes rule) or an 'm' (for a simplifying model assumption) above any equality sign." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Working out the example problem: Disease Diagnosis\n", "\n", "- **Problem**: Given a disease $D$ with prevalence of $1\\%$ and a test procedure $T$ with sensitivity ('true positive' rate) of $95\\%$ and specificity ('true negative' rate) of $85\\%$, what is the chance that somebody who tests positive actually has the disease?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- **Solution**: The given data are $p(D=1)=0.01$, $p(T=1|D=1)=0.95$ and $p(T=0|D=0)=0.85$. Then according to Bayes rule,\n", "\n", "$$\\begin{align*}\n", "p( D=1 &| T=1) \\\\\n", "&\\stackrel{p}{=} \\frac{p(T=1|D=1)p(D=1)}{p(T=1)} \\\\\n", "&\\stackrel{s}{=} \\frac{p(T=1|D=1)p(D=1)}{p(T=1|D=1)p(D=1)+p(T=1|D=0)p(D=0)} \\\\\n", "&= \\frac{0.95\\times0.01}{0.95\\times0.01 + 0.15\\times0.99} = 0.0601\n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Note that $p(\\text{sick}|\\text{positive test}) = 0.06$ while $p(\\text{positive test} | \\text{sick}) = 0.95$. This is a huge difference that is sometimes called the \"medical test paradox\" or the [base rate fallacy](https://en.wikipedia.org/wiki/Base_rate_fallacy). \n", "\n", "- Many people have trouble distinguishing $p(A|B)$ from $p(B|A)$ in their heads. This has led to major negative consequences. For instance, unfounded convictions in the legal arena and even lots of unfounded conclusions in the pursuit of scientific results. See [Ioannidis (2005)](https://journals.plos.org/plosmedicine/article?id=10.1371/journal.pmed.0020124) and [Clayton (2021)](https://aubreyclayton.com/bernoulli)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Inference Exercise: Bag Counter\n", "\n", "- **Problem**: A bag contains one ball, known to be either white or black. A white ball is put in, the bag is shaken,\n", " and a ball is drawn out, which proves to be white. What is now the\n", " chance of drawing a white ball?\n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- **Solution**: Again, use Bayes and marginalization to arrive at $p(\\text{white}|\\text{data})=2/3$, see the [Exercises](https://nbviewer.org/github/bertdv/BMLIP/blob/master/lessons/exercises/Exercises-Probability-Theory-Review.ipynb) notebook.\n", "\n", "- $\\Rightarrow$ Note that probabilities describe **a person's state of knowledge** rather than a 'property of nature'." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Inference Exercise: Causality?\n", "\n", "- **Problem**: A dark bag contains five red balls and seven green ones. (a) What is the probability of drawing a red ball on the first draw? Balls are not returned to the bag after each draw. (b) If you know that on the second draw the ball was a green one, what is now the probability of drawing a red ball on the first draw?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- **Solution**: (a) $5/12$. (b) $5/11$, see the [Exercises](https://nbviewer.org/github/bertdv/BMLIP/blob/master/lessons/exercises/Exercises-Probability-Theory-Review.ipynb) notebook.\n", "\n", "- $\\Rightarrow$ Again, we conclude that conditional probabilities reflect **implications for a state of knowledge** rather than temporal causality." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Moments of the PDF\n", "\n", "- Consider a distribution $p(x)$. The **expected value** or **mean** is defined as \n", "$$\\mu_x = \\mathbb{E}[x] \\triangleq \\int x \\,p(x) \\,\\mathrm{d}{x}$$ " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- The **variance** of $x$ is defined as \n", "$$\\Sigma_x \\triangleq \\mathbb{E} \\left[(x-\\mu_x)(x-\\mu_x)^T \\right]$$ " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- The **covariance** matrix between _vectors_ $x$ and $y$ is defined as\n", "$$\\begin{align*}\n", " \\Sigma_{xy} &\\triangleq \\mathbb{E}\\left[ (x-\\mu_x) (y-\\mu_y)^T \\right]\\\\\n", " &= \\mathbb{E}\\left[ (x-\\mu_x) (y^T-\\mu_y^T) \\right]\\\\\n", " &= \\mathbb{E}[x y^T] - \\mu_x \\mu_y^T\n", "\\end{align*}$$\n", " - Clearly, if $x$ and $y$ are independent, then $\\Sigma_{xy} = 0$, since $\\mathbb{E}[x y^T] = \\mathbb{E}[x] \\mathbb{E}[y^T] = \\mu_x \\mu_y^T$.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "### Linear Transformations \n", "\n", "- Consider an arbitrary distribution $p(X)$ with mean $\\mu_x$ and variance $\\Sigma_x$ and the linear transformation $$Z = A X + b \\,.$$ \n", "\n", "- No matter the specification of $p(X)$, we can derive that (see [Exercises](https://nbviewer.org/github/bertdv/BMLIP/blob/master/lessons/exercises/Exercises-Probability-Theory-Review.ipynb) notebook)\n", "$$\\begin{align}\n", "\\mu_z &= A\\mu_x + b \\tag{SRG-3a}\\\\\n", "\\Sigma_z &= A\\,\\Sigma_x\\,A^T \\tag{SRG-3b}\n", "\\end{align}$$\n", " - (The tag (SRG-3a) refers to the corresponding eqn number in Sam Roweis [Gaussian identities](https://github.com/bertdv/BMLIP/blob/master/lessons/notebooks/files/Roweis-1999-gaussian-identities.pdf) notes.)\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### PDF for the Sum of Two Variables\n", "\n", "\n", "- Given eqs SRG-3a and SRG-3b (previous cell), you should now be able to derive the following: for any distribution of variable $X$ and $Y$ and sum $Z = X+Y$ (proof by [Exercise](https://nbviewer.org/github/bertdv/BMLIP/blob/master/lessons/exercises/Exercises-Probability-Theory-Review.ipynb))\n", "\n", "$$\\begin{align*}\n", " \\mu_z &= \\mu_x + \\mu_y \\\\\n", " \\Sigma_z &= \\Sigma_x + \\Sigma_y + 2\\Sigma_{xy} \n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Clearly, it follows that if $X$ and $Y$ are **independent**, then\n", "\n", "$$\\Sigma_z = \\Sigma_x + \\Sigma_y $$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- More generally, given two **independent** variables\n", "$X$ and $Y$, with PDF's $p_x(x)$ and $p_y(y)$. The PDF $p_z(z)$ for $Z=X+Y$ is given by the **convolution**\n", "\n", "$$\n", "p_z (z) = \\int_{ - \\infty }^\\infty {p_x (x)p_y (z - x)\\,\\mathrm{d}{x}}\n", "$$ " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- **Proof**: Let $p_z(z)$ be the probability that $Z$ has value $z$. This occurs if $X$ has some value $x$ and at the same time $Y=z-x$, with joint probability $p_x(x)p_y(z-x)$. Since $x$ can be any value, we sum over all possible values for $x$ to get\n", "$\n", "p_z (z) = \\int_{ - \\infty }^\\infty {p_x (x)p_y (z - x)\\,\\mathrm{d}{x}}\n", "$ \n", " \n", " - Note that $p_z(z) \\neq p_x(x) + p_y(y)\\,$ !!\n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- [https://en.wikipedia.org/wiki/List_of_convolutions_of_probability_distributions](https://en.wikipedia.org/wiki/List_of_convolutions_of_probability_distributions) shows how these convolutions work out for a few common probability distributions. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- In linear stochastic systems theory, the Fourier Transform of a PDF (i.e., the characteristic function) plays an important computational role. Why?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Code Example: Sum of Two Gaussian Distributed Variables\n", "\n", "- Consider the PDF of the sum of two independent Gaussian distributed $X$ and $Y$:\n", "\n", "$$\\begin{align*}\n", "p_X(x) &= \\mathcal{N}(\\,x\\,|\\,\\mu_X,\\sigma_X^2\\,) \\\\ \n", "p_Y(y) &= \\mathcal{N}(\\,y\\,|\\,\\mu_Y,\\sigma_Y^2\\,) \n", "\\end{align*}$$\n", "\n", "- Let $Z = X + Y$. Performing the convolution (nice exercise) yields a Gaussian PDF for $Z$: \n", "\n", "$$\n", "p_Z(z) = \\mathcal{N}(\\,z\\,|\\,\\mu_X+\\mu_Y,\\sigma_X^2+\\sigma_Y^2\\,).\n", "$$" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd3xUVdoH8HNumZn0Xie9ESAk9B4CVkSFtS0iRUFZBVREkK5SBMFFRbDyCmtZ3VUUKyCd0EMLNb0nECC9z9x23j/CBggpk2Rm7szc5/vhDzK5ufeZTDK/nHNPwYQQBAAAACgVJXcBAAAAgJwgCAEAACgaBCEAAABFgyAEAACgaBCEAAAAFA2CEAAAgKJBEAIAAFA0CEIAAACKBkEIAABA0SAIAQAAKJoMQXjhwoWvvvrK8ON5njdZLRZNsU9cEAS5S5CHIAjKXPJQFEVlPnFJkiRJkrsKGUiSJIqi3FXcIkMQXrp0ac+ePYYfr9PpTFeMJVPsE9fr9cp8W+Q4TplvizzPW9TbotnwPK/MP/tEUbSoJw5dowAAABQNghAAAICiQRACAABQNMbwQ8+dO/frr7/a29tPnjzZz8/v7gOysrK2bdtWV1c3cODAMWPGYIyNVycAAABgEoa2CA8fPpyQkEDTdE5OTv/+/cvKypodcObMmT59+lRWVnp4eLz++usLFiwwdqkAAACA8RnaIlyzZs3ixYsb462wsHDz5s3z58+//YBt27Y9/PDDq1evRghFRUVNnz79vffeM3q5AAAAgHEZ2iI8cODA6NGjG/8/evToAwcONDsgKioqOzub4ziEUEpKSlRUlBGrBAAAAEzEoBZhZWVlQ0ODt7d344c+Pj5Xr15tdsyUKVOSk5ODg4Pd3d0RQm3MFLx+/frJkydfeOGFxg9pmp41a1YbwanT6ViWNaROG6PkJ05RlALvMet0OkIITdNyF2JuOp2OYRiG6cCQBdug1+sxxgqcPMrzvCRJ5pkurFKpKKqdJp9BP3mNP6BNr5YgCHe/Qf/++++///775s2b/fz81qxZM3v27K1bt7Z4No1G4+bm1q9fv6ZHXF1d2/jlp2lagW8NSMFPnKqv5jNOUfZO6h4D5a7FrOj/kbsQc1PyE8cYK/CJS5JktiduyJ/UBgWho6Ojk5PT1atXtVotQqi4uPjuUaOfffbZjBkzxowZgxDauHGjj4/PtWvXfH197z6bi4tLZGTkjBkzDLk0QohlWWU2jBT1xCVdnT7trC7znD49WWyoUYf2FIpzhezzLuP+gWmlNBQaX3EFvi2KoqjMFmFjHijn1/x2kiRZzhM39Cfv4Ycf3rZt24ABAyRJ+uWXXyZNmoQQEgTh1KlTffr00Wg0rq6uhYWFjQcXFhbSNO3s7GyqqoFtIfqGkvWvUc4e6vBebuNf5Rw97Oztib6+6pdNJR/N9Zi6hHbzlrtGAKzPsWPH1q9fL3cVLSCEEELa7bFs5o033hgwYIAp6jE0CBcvXjxq1Kji4uKrV6/W1tZOnjwZIVRRUTF06NDMzMyIiIh58+Y9+OCDJSUl3t7eW7duXbRokb29vSkqBraGkPL/fMAGRrk8+nzjA1x9PUKI0ji4Pf1aXdKuGx+86vbMPE33/rJWCYD1SU5Orq6ufv755+UuxAg+//zzCxcuyByEvXr1unz58u7dux0cHB566CE7OzuEkKur68GDBwMCAhBC/fv3z8jISExMrK2tnT59ekxMjCnKBbanZv9PQulVj2lvtfA5jB0Gj2Z9gyv+84Fj/Fin+582e3UAWLeIiIinnnpK7iqMYNeuXaY7eQc65X18fBobgk1Ylk1ISGj60MPD4/HHHzdaaUAB9Jnnag9u83hxBWZavVugCunu8eKqsi+WqKP7qgJhWg4AwMhgrVEgG7GypPybNS6PvUg7e7R9JO3k4njPU5U/fYoUuT0TAMCkIAiBPIjAl21eYT/0YXVErCHH2/cZgQS+4dxhUxcGAFAaCEIgj8pfPqdcPByHjjH0CzB2evCZqt+/JDxnyroAAIoDQQhkIJRebTh32PVv/0AdWT5GFdKd8Q+rPbjNdIUBABRIcTNYgSWoPfizw4B7sUrT0S90euDpsi/etB/0AO3sborCALBtZ0vJfTsFs91pV1HowMNMD9fmf+9eu3bNwcHByckpMzMzJCRE9pn1EITA3KT6mvqziZ6zOrM5CePmbd/vnurtX7tNmGP0wgCweX09cdJYVjRXEjIUirhrYZWampqjR49u2bJl2rRp/fv3nzZt2jfffCPv2sIQhMDcag//oYnuTzu5du7LHRPGlWyYxxVmwFQKADrBVYXM1iRkW7r5lpmZ+dBDD61evfqxxx6jKCo7O7tp/U65wD1CYFZE4OuO/OEw5KFOnwGrNI6jnqj6dZMRqwIAmE3fvn1zc3NjYmIoihJFMTs7283NTd6SIAiBWdWf2c/4BTM+gV05iX2fBLGqjMtLNVZVAABzOnjwYONiLPv373/44YdlX48TghCYESG1B39xHPpwV89DUXZ9EupOtrrnJQDAkh04cKC+vv7ixYtbt2798MMP5S4HghCYkS7tNJIkVWiPrp/KLi6+4dwhwum7fioAgDkRQnJzc1944QWWZb/44gsXF5fGx3U6XUNDQ01NjSiKZi4JghCYT83+nxyGPdKhuYOtoZ3d2ICIhkvHu34qAIA5nTt3LioqSqPRREdHNw0WvXbtWkpKysqVK48dO/bNN9+YuSQIQmAm/NVc4XqhptdgY53QLm54fdJuY50NAGAGHMclJyffc8892dnZtz/u6OjYt2/f+vr6Bx98cOrUqWauCqZPADOp2fej/eDRRtxuXhM9oHr712JlCe3qZaxzAgBMSqVSTZs27e7Hi4qKPD09KYqqr6+vqKgw82wKaBECc5BqK3WpJ+37jTLiOTHLamIG1Z/ca8RzAmDzMEYYme+fgTIyMs6dO3fPPfekpKSYf04htAiBOTRcOKaOiKPsHIx7Wvu4EZW/fO50/9NGue8IgBJ4dnhlQ3MYO3asjFeHFiEwh4ZzhzU9Bhr9tGxgBKIomFAIgHVpGhdq/gGiLYIWITA5qb6WK0h3Hf+qKU5uFzes7tReo0zJAMDmcYWZpZ8tQshsy25jr1c/YH2D7qiB477//vv9+/c/+eSTzs7OO3fuXLt2rbnqaRkEITC5hovHVOG9MKs2xcnt4kaUfrrA9W8vYpVJzg+ALVEFRnrN3YgEwTyXwyo14+bd7MGzZ8+OGzduw4YNmzZt0mg0S5cuvX79uo+Pj3lKahEEITC5hnOHNTGDTHRy2tmNDYhsuHjMuCNxALBVlNoesebqkGxplPjgwYNzc3PDw8M1Gg0hJDs7W6OR+b4l3CMEpkX0DVzuZXVUH9Ndwq43TCgEwJocPHhw5MiRCKEzZ87079+/aXEZuUAQAtNquHSCDe5Oqe1MdwlNt/5cUZZYWWK6SwAAjOjAgQNqtbqmpuajjz769NNP5S4HghCYWMO5w5oeA0x6Ccyymu79G84fNelVAADGcvny5dGjRx89evTzzz8PDOzSXjRGAUEITIhwen3mOU10X1NfSB3ZR5dy0tRXAQB0XVZWVnBwcEBAwOjRox0cbs4tvnTp0vvvv19WVvb7779nZGSYuSQIQmBCutSTbGAkZedo6gupI2K4/DSibzD1hQAAXaHX63fv3j148OCsrKzbH/fy8iooKHB2dj506FBYWJiZq4JRo8CEGs4d1nTvb4YLYZWG0YbpMs/bxRhtUW8AgNGp1eqZM2fe/biPj49Op2NZlud5hjF3MEGLEJgKEXhd+hlNd9PeIGyijojTpZ4yz7UAAEY3YMCA3377zfwLjSIIQmA6urQzjE8w5eBsnstpovroLiUhYrYlMwCwThib7V/TdoPt2rt379ChQwsKCqZMmWLSZ98i6BoFpmKi9UVbw3j5Y4bmr+WzfiFmuygAVoc219+mHdK/f/+8vLzHHnvM19fX/FeHFiEwCSIKupQk89wgbKKO6A1jRwGwRq6urr179w4ICJDl6tAiBCahzzxPe/rTzu7mvKg6Kq7u6A6ne/9uzosCYEUyyrOn75xjtstRGG95eGOoS1CzxzMzMx0dHT08PM6ePRsTE+PoaPKB5W2DIAQmoU8/q46IM/NFVSE9K7d+LDXUmmHCBgDWKMo9/LcnvhWJmdYaZSjGRd28J7aioiInJ+eDDz6YOXPmfffd9+yzz/74448UJWf3JHSNApPQpSdrwmPMfFHMsmxwtD7trJmvCwAwXGFhYUJCQllZ2dixYx0cHIqKiq5evSpvSRCEwPikumqx/BrjH2r+S6sj4+A2IQCWLDY2Ni0trXfv3hhjnudzc3M9PT3lLQmCEBifLiNZFdIdt7QDi6lpuvVpSDkFkygAsGQHDx6Mj49HCO3cuXP8+PGwDROwQfqMc6owc/eLNqJdPGlHJ64wU5arAwAMceDAgdra2gMHDuzfv1/27ekRDJYBpqBPP+v2zFy5rq6O7K27nKQKipKrAABAG0RRvHLlysyZM0tLS0eNurWfdm1tLcdxdnZ2dnYm3LWtRdAiBEYmlF8nnJ7xkmGdpEbqyN661NNyXR0A0LbGZbV5nvfy8rr98Z9++mnr1q0HDhwwf0kQhMDI9GlnVBExyOCllYyODYoSSorEmgq5CgAAtIbjuLKysvHjxxcWFjb71JNPPhkVFTVmzBjzVwVBCIxMl5GsCpXnBmEjTDOqsBh92hkZawAAtEilUj355JNPPPFEeHj47Y/X1tb+9NNPw4YNO3bsmPmrgiAERkWIPuu8OrynvFWoI2J16TCbEIAWUJgy5z8DqyoqKhoxYsTy5ct79Ohh0qffIhgsA4yJv5pDaRxpZw95y1CFdK89+Iu8NQBgmZzVTnKX0ILo6GiE0KpVq2S5OrQIgTHp0pPVZl9Q5m6Mhy8iklh+Q+5CAABWAIIQGJNe7huETVTB0frsC3JXAQCwAhCEwGiIKHB5qarQaLkLQQghNribPguCEADQPghCYDRcbirj6WchOz+oQ7rrsy7KXQUAwApAEAKj0Wcmy7Wy2t0YL62kqxeryuQuBABg6SAIgdHo0pPVFhOECGNVcDd9NjQKAQDtgOkTwDiIvoEvzmMDLWiFT1VwtD77on3fkXIXAoA8VCrVt99+u3PnTrkLMYKSkpKEhAQTnRyCEBiHPus8GxCBWVbuQm5Rh0ZX/vyZ3FUAIJvnnnvu3nvvlbuKFgiCIEmSSqXq0FcFBQWZqB4IQmAcuozz6jCZF5RphvEJFqsrpNpKytFV7loAkAHLsmFhYXJX0QKe5yVJUqvVchdyE9wjBMbBZV9SB3eTu4o7YcwGd9NnX5K7DgCARYMgBEZAOB1/o5Dxs7i/PVVB3WBaPQCgbRCEwAi4/HTGN9iibhA2UofCbEIAQDsgCIER6HNT1EGRclfRAtYvRCy/LtVVy10IAMByQRACI+ByLrGBFnaDsBFFs4GR+twUuesAAFguCELQZYRwBekqi2wRIoRUwdEcTKsHALQOghB0FX8tn7J3ohyc5S6kZeqQaH0mjJcBALQKghB0FZdzWRVkkf2iCCGEGG04X1Ik6erkLgQAYKEgCEFX6XMuWdTKas1gmmG14VxuqtyFAAAsFAQh6CouN0UVZLlBiBCC1bcBAG2AIARdIlaXS7p6xtNP7kLaog7tAZv0AgBaA0EIuoTLvcwGdUMYy11IW1htOH81lwi83IUAACwRBCHoEn1OiiooQu4q2oFZFe3py1/JlrsQAIAlgiAEXcLlXFJZ5lT6O6kDIrm8NLmrAABYIghC0HmE5/jrhay/xa21fTc2IILLg/VlAAAtgCAEncflp7G+QRa41vbd2EBoEQIAWgZBCDpPn3vZwidONGHcfSR9g1hVJnchAACLA0EIOo/LvmzJU+nvgDEbGMkVpMtdBwDA4kAQgs4ihCtIYwMsdK3tu7HacOgdBQDcDYIQdBJ/LZ+yc6SdXOQuxFCqQBgvAwBoAQQh6CQuN4UNtoKJE03YgAiuKBtJotyFAAAsCwQh6CR9ziVVoNX0iyKEKLUd4+rFX82VuxAAgGWBIASdxOWlqaznBmEjNiBcD7cJAQB3Ygw/tKCgYO/evS4uLo888oharW7xmAsXLpw8edLJySkhIcHX19dIRQKLI9XXSLUVjHeA3IV0DBsQyeWloOGPyF0IAMCCGNoiPH78eFxcXFJS0ieffJKQkMBx3N3HLFq0aPTo0UePHv311183btxo1DqBZeEKMhj/cAtfa/tuqsAIGDgKAGjG0Bbh8uXLFy5cuGDBAkEQ+vXr9/PPP0+YMOH2A/bu3bt58+ZLly55e3uboE5gWbj8NFVAuNxVdBjjpZXqqqS6asrBWe5aAACWwqAWIcdxe/bsefzxxxFCDMOMHTt2x44dzY754YcfpkyZUltbu2PHjitXrhi/UmBJuPx01t/6ghBhzGrDuXxoFAIAbjGoRVhcXCxJklarbfxQq9UeO3as2THZ2dkcxx08eDA6Onry5MkbNmyYOHFii2errKxMSUlZvXr1zQoYZvz48f7+/q1dned5nlfiTnKW/MS5ggz7B54xUXmCIAiCYIozI4Ro/zBd9iU6so+Jzt8VPM9TFCVJktyFmBvP84QQQojchZgbz/MYY2xttxi6jud5SZIoyhyjNRmGafc7bFAQNv5mNp2Loqi736f0en1dXd2ZM2domv7jjz+effbZCRMmtPg8OY7T6/Xl5eVNj9TU1Ihiq7O7RFFs47M2zGKfuFh+HSGEHFxM9Jbd+MRN9O5A+4XqTu/TWOY31lJfcVNrfLkVmAeNT1yBL7ooipIkmeeJ0zRtnCD08/NDCN24cSM4OBghdO3atbsbcP7+/p6enjRNI4RGjBhRUVFRXFzc1Ii8nbe3d58+fdatW2fIpRFCPM9rNBoDD7YlFvvE66/nqYIiWxs53HWiKKrVahO9LbJh3et++0KjUiGz/DXaIZIkqdXqxl8ipWEYhmE6MIjdNjTGv0qlkrsQc6NpuvGnXe5CbjLovUCj0QwdOnTnzp2NH+7ateuee+5BCImiWFxc3Nihcd9996Wn31zROC0tTaPRwKgZW8UXZLBaK9iDsEWUnSPl6MpfL5C7EACApTD0T7AlS5ZMmjSptLQ0NTW1pKTkmWeeQQjl5ORERUWVlJR4enpOmjTpww8/nDp1amxs7CeffPLmm2+y1rBNHegELi/VceRjclfReWxQJJeXyvqFyF0IAMAiGNo7NGbMmD179oiiOGjQoBMnTjg4OCCEfH19v/rqKycnJ4SQg4NDUlLSgAEDamtrv/zyy8WLF5uwaiAjSeKu5jLWsCt9a1htGJeXKncVAABL0YFO+X79+vXr1+/2R5ycnJ599tmmD11cXGbOnGm00oBF4otzGRcPSm0ndyGdpw6MrDy5V+4qAACWwuLGCwALxxVksFornEF4G8Y7UKgskXT1chcCALAIEISgY7j8NCbAivtFEUKIolnfEL4wU+46AAAWAYIQdAyXl6a2tk0n7sZqw7iCDLmrAABYBAhC0AGE0wkVN2hr23Tibqw2jMtPl7sKAIBFgCAEHcAVZLK+QZi2+onPrDacK4QgBAAgBEEIOoTLT7PKtbbvwrh5E71Oqq2UuxAAgPwgCEEHcPnprLWPlGmEMesfyhXAeBkAAAQh6AiuIN3a5040YbVhXAH0jgIAIAiBwcTqCqJvYNxsZAlZ1h8GjgIAEIIgBIbj8tPYwEhkK3vlsNpwGDgKAEAQhMBwfEGGKsBG+kURQrSzG6ZpsbJE7kIAADKDIASG0uelsta81vbdYFo9AABBEAJDEcJfybTebQhbxPiFQhACACAIgUGE0quUxpFycJa7EGNSacO4/DS5qwAAyAyCEBiEK8hgbKs5iBBiA8L5oixEiNyFAADkBEEIDMIVZrK+IXJXYWSUnSNl5yCUXpW7EACAnCAIgUH4/HSVbawpcydGGw7T6gFQOAhCYABC+OJcxi9E7jqMj/ULhSAEQOEgCEH7+OsFlKMLpXGQuxDjU2nDuXwYOAqAokEQgvbxBRk2NoOwCasN5YvzkCTKXQgAQDYQhKB9XFEm4x8qdxUmgVUaxsWDL86XuxAAgGwgCEH7uPwMlc3NnWjCBIRzhbAfEwDKBUEI2iNJ/LV81hZHyjRi/WE/JgAUDYIQtIMvzmNcPLFKI3chpqLyD4VtKABQMghC0A6uMJOxxRmETRjfYKHkCuE5uQsBAMgDghC0g7fFNWVuhxmW9fLnr+bIXQgAQB4QhKAdXH66DY+UacT4h8J4GQAUC4IQtIWIAn+jgPENlrsQ02L9Q3nYjwkApYIgBG3hr+YwHn6YVcldiGmp/GGHXgCUC4IQtIUvyGRtdCr97WjvAKHiBuF0chcCAJABBCFoC1eoiCDENMN6afkrMF4GACWCIARt4QoyWFsfKdOIhfEyACgVBCFoFeE5ofQq4x0odyHmwPhBEAKgUBCEoFX8lRzWS4sZVu5CzEGlhYGjACgUBCFoFVeQzmht/wZhI6ZxvIy+Qe5CAADmBkEIWsUVZrJ+SglCRNGsbyAH42UAUB4IQtAqviBDpZgWIUKI9QvjC6F3FADFgSAELSOcXqi4QXsFyF2I+TD+oTCtHgAFgiAELeOKslifQEwzchdiPiqYQQGAIkEQgpbxhRmsvyJmEDZhvLRiVZmkq5e7EACAWUEQgpZxBRmMAtaUuQNFsb7BfFG23HUAAMwKghC0jCvMVCktCBFi/UK4IugdBUBZIAhBC4i+QawqY7z85S7E3GA/JgAUCIIQtIArzGJ8ghFFy12IuTGwHxMAygNBCFrAFSprBmET1stfqq2QGmrlLgQAYD4QhKAFfEGGgtaUuR3GjE8IjJcBQFEgCEELuMJM5awy2gyrhdmEACgLBCFoTtLViTXlrKfiRso0Yv1DYLwMAIoCQQia4wszWb8whLHchciD9QvjYMVRAJQEghA0xxVkqvxD5K5CNoynn1RXI9XXyF0IAMBMIAhBc1yh8taUuR3GrF8wX5Qldx0AADOBIATN8YWZSltltBkWVt8GQEkgCMEdpPpaqa6a8fCVuxA5Mf6hXD7cJgRAKSAIwR1u9osqdaRMI1YbDjv0AqAcEITgDrwi19puhnHzlnT1Um2V3IUAAMwBghDcQYm7L90NYwZuEwKgGBCE4A4wUqaRSgurbwOgFBCE4BaprlpqqGPcvOUuRH6sfyjcJgRAISAIwS1cQQajVe6aMrdjYT8mABQDghDcwhXBSJmbaFdPIopiVZnchQAATA6CENzCF2SwEIT/w2rDYLwMAEoAQQhu4QoyYchoE5U2lIcgBEABIAjBTVJtFeF0jKuX3IVYCsYvjCtIl7sKAIDJQRCCm7iCDDYgHEbKNGG14bAfEwBKAEEIbuIKM1m/ELmrsCC0kwumGLH8htyFAABMC4IQ3MTlp7H+4XJXYVlYLWzSC4DtgyAEN/GFmWwArClzB8YPFloDwPZBEAKEEBIrSxAitLO73IVYFpU2lIdp9QDYOghCgNDNkTIRcldhcVhtOFeUgQiRuxAAgAlBEAKEbm46Af2izVH2TlhtL5QVy10IAMCEIAgBQghx+WkqLQRhC2DRUQBsHgQhQIgQ/koWLK7WItY/hC/KkrsKAIAJQRACJJRcoTROlL2T3IVYIpU2nMtPk7sKAIAJQRACxBVksNAv2gpGG8ZfyYHxMgDYsA4EYWVl5f79+y9evNj2YSUlJfn5+V2rCpgVV5AOa223hlLbUw7O/I1CuQsBAJiKoUF44sSJqKiodevWjR079oUXXmjtsPLy8ri4uAEDBhipPGAOXH46jJRpA6sN4wtgWj0ANsvQIFy4cOH8+fN37Nhx5syZP/74IykpqcXDXnvttbFjxxqvPGB6kshfy4eRMm1g/EK5QtiGAgCbZVAQlpaWJiYmTpkyBSHk7u7+8MMPb9u27e7DduzYUVpaOmHCBCPXCEyJL85nXDywSiN3IZZLpQ3j8mEGBQA2izHkoKKiIo1G4+3t3fhhUFBQZmbznqLq6uq5c+fu2LGjoKCg7bM1NDQUFBT8+OOPTY/cd999rq6urR0vSZIkSYbUaWPM88T1+WmMNsyivsONTxxbzIZQtF8Ify1P5DlMG/T70mmW9sTNRvofuQsxt8aXW5lP3GyvOEW1394z6Bdbp9OpVKqmD9VqdUNDQ7Nj5s2bN2PGjNDQ0HaDsKysLD8//4cffmh6JDg4uFevXq0dr9frWZY1pE4bY54n3pCbSnkH6fV6U1/obiJBhbUoqxZn11BZNahMj3gJ6UTcwLM8kmiMvDXEU4P87IiHGmntSawb8ZGp4Uq5etXlpTNa0+7OodPpCCE0TZv0KhZIp9MxDMMwpv07wwLp9XplBiHP85IkEbMMxtZoNO1moUE/eX5+fjU1NRzHNcZhSUmJn5/f7QcUFhZ+++23L7/88sKFCwsLC+vq6hYuXDhv3jxPT8+7zxYQEBAfH//dd98Z+DREUbS3tzfwYFtinidecyXbqW8Ca2dn6gs14iV0uoTsL5YSi0lWNfG2w6FOKNQRx3pidw1WYaShERH0jnY0Qfh6AyrXkRIdyq4jiTfIomRiz+D+Hri/Fx7ghXu5m6/pxAVEUDfy7SNb/XPNWNRqtQKDkKIoZQYhTdMY49ubGQrRGIRqtVruQm4y6CcvICBAq9UmJibef//9CKHExMQ5c+bcfoCDg8OyZcsa/19ZWYkxdnNzM6RBCuRFeE4ovcr4BJr6QjU82l4o7b1CjlyTgpzwMG/8Riwd7YJVLb3nNzQQOw1GGEc6I4TuCLuCGnK+glwoI//KkPQS+lswfiyE6uFq8kBktWF8Xhoa9oipLwQAMD+DgpCm6Tlz5syaNWv58uXHjh0rLy9/8sknEUInTpwYOXKkTqdzd3dfsGBB48GJiYnbtm1r+hBYMv5KNuulxYwJO2AvV5CvMqU/C6SBXlSCH57bi3Xvwl+BQU44yAk/GoQQQpnVZFeh9FyiaEejccHUxAjsY2eqRFQFhNcn7TbRyQEA8jK0L2LOnDk+Pj47d+709/c/cuSIRqNBCAUEBCxevLjZkWFhYW+++aaRywSmYbo1ZTgR/VYgfZ0pXatHj4dQP9/HelAfd6YAACAASURBVBi7FyTSGUf2pF/uiS6Ukx2F0sjt4phAalZ3KszZ+HHIeAcJlaWSrp7SKLGXHgDbZmgQYownTpw4ceLE2x8MCAh46623mh0ZGBj4yiuvGKc6YGJcfhobGGXcc4oE/ZQrrbsohTigZyOp4T4UZeKey1h3HOtOz+xB/5AtjdsjDPahXu5Bxbkb9aoUxfqG8IUZ6sjexjwtAMACKO7uNLgdV5DhMHi0sc5GENpZKK09L7mo0ep+dKyHWacBOLNoejQ1KYL6NV+adkjs6YqW96VDjdc6VAWEcfnpEIQA2B4IQuWSdHVidTnj5W+Us524QZYni5yIXouhh/rINhPOjkETwqmnQqn/5EiP7hGeDqNei6EdjXEPlNGG6VNPG+FEAAALAwM7lYsvzFT5hSKqq4P1q3k076Q486j4TBj1/ShGxhRswlBocgS19V72Sj2K3y5szZW6Pl+JDYjg8mA/JgBsEAShcnEFGay2q0uMbi+QErYLooR+uo95IMCyfpzc1WhZX/r9gfTmdGncHiGvpktpyLh6IULEqjJjlQcAsBCW9c4FzInLz+jK7kvXG8jUQ+K756W1/emFcbSDpfay93TDXycw9/pRj+wWvs3qUtOQCQjn8mH1bQBsDQShcnF5KWxAROe+dleRdP9OIdQR/XcUE2feQTGd83Q4tTme+SZDmnxAuKHr5ElY/1CuAIIQAFsDQahQYmUpEUXGzbujX8hL6O0z4tIz0geDmBe706z1LAcW4oS/SmAinfH9O/jtBZ1Z3VEFtwkBsEWW2p8FTIzLT1UFRXb0q/JryUtHRU81/s8oxskKF0JnKPRSD3q4H7XktHiihLzdh2Y68qcgqw3jizIRIUh5G0QAYMOgRahQXH4628G9FP4skB7ZLTwYQK0bRFtjCjaJccPfjWQyK8nf9wulHekmpewcKUcX/nqhyUoDAMgAglChuNwUVaChNwgJQmvPiyuTpY1DmAlhtvAz48ii9UOYWDc8epeQXNaBATSsNhxuEwJgY2zhTQ10FBEF7mou42/QKqN6Ec06KiZeI18lMN1Nv8+D2WCMZvSgF8bRzx4Uvssy9JYh6x/Gw8BRAGwLBKES8VdzGDdvSt3+HoTXG8i4PYJE0KZhTFd2jbBY8T74yxHMZ6nS8rOiIVMrVAERXD6MlwHApsBgGSXi8tNZA/pFUyvJlETx0UD8YnfrGRvaccGO+OsEZk6SMP2I+MkQWtPm7wTjF8zfKCI8h1nFbabaRdfqbmRX5F2tvVbWUF6mq6hsqCqtL2sQdAghQgjCNwcgOWkc3dWurnauXnYebhqXQGdtuGuoh52bvMUD2wZBqERcXmq7Mwj3XiVzTghL4uhR/rbfbeDIok+HMG+eFZ85KPxrBOPSesZhhmW9tPyVbFVIdzMWaJXKGipOFyenlmdkV+RlV+ZpaE2oS5Cfo4+rxrmbW4SLn5ObxtWe0SCE9BxH0zRD0wihWq6+Ul9Vpasu11VmV+YlFhzLrSrACIe7hkS5h3f3jOrnG+ekcpT7yQGbAkGoRFxemv3gh9o4YFuetPyMuH4IE+NmOzcF28bS6N3+9AeXxHG7he9G0VqHVp84GxDG5adBELZIJ+gvlFw+VXzu1NWzpbryOO+eEW5hj0c9Euoa1EZ66Rk9TdMMwyCEfBxaOKBcV5lXVZBfVfhn1q61JzYGO2sH+vft79u7p1c0jW25uwKYBwSh4kh11VJtJdv6phNb0qWPU6TP45kwJ6WkYCOM0dxe9HfZ0tg94r9H0q2NDGL9w7kCuE14B5GIJ68m7847cOLK6XCXkFjvnjP6Tg1zDaGMNOHSXePqrnHt6xOLEOIlIaMs61zJ5Y/ObLpRV5oQNOzB0JExXt0xUtaPKzAiCELF4fLTGW14a1PC378k/ZQjbY6n/VtvEtm2ieGUlxo9vV/4diQT29LuvkxAeO2h38xfmGVKKc3YnXfgQP5hf0e/+MBBz/Ya78Satt+SpZieXtE9vaIn9niirKH8cOGJdUmfNAi6B0JHPRA6Msg5wKRXBzYJglBxuPxUdUtryhCEVpwV918lX8bTnhqFpmCjBwIoexZPOihsiWf6ezX/VrCe/lJ9jVhTQTspdwSHIIkHCg7/mPpbHV+fEDRkzci3vO09zV+Gh53736LG/C1qTF514eGC46/tWRLsEji+x98G+feDBiIwHASh4nC5qfYD72v2oETQ3CQxq5psGcEYZRtbazfcB6/oy0w9LGwaxgxptsMixmxgBJeXZtdriEzVyamOr9+RtfeHtF+87DzGRo4e7N+PwvIPpwpxDgyJCZzQ44mTxWe3nP/+o9Obnoh69JGI+zWMRu7SgBWAIFQYQriiDJfHX7rjMYQWnRYzq8nHgxk7SMH/GeKD3xvAvHhU2DiUSfC9IwvZgAguL1VpQVilr/735Z92ZO8Z6N9v0dDXQpwD5a6oOYaih2oHDNUOuFyatj17z7eXf3iq27gnox+FOARtgyBUFv56IWXnRNk7NT0iEfTGSTGzmnw8hLGDH4c79fXE6wYyLx8VPhzC3Od/KwvVQVG1iQq6TVjH1/+Q+uu2jD/jA4dsuG+1q8ZF7ora0dMzuqdn9NWaaz+k/Trh939M6vn3sREPsjT8lQdaBu98ysLlp94+lZ4gNP+UmFFFNkIKtiLOA380hJl9XPh42K12IRsQwV3NIaKAaRv/rukE/S8Z2/+buq2fb+91o5Z5yXEjsNP8nXznDHgpr7rwvynb/pv6y9ReEx4MGwXTLcDd5O/cB+bE5aez2ptBSBBafFpMrSAbIAXb1NMNfzCYefmYcOz6zUXYsErDuPvyV7LlLczUDhYcnfznjAs3UlbGL57Vd5p1pWCTEOfAhYNnvz7gpZ05e6ftmH3uxiW5KwIWB97/lIXLTXGJHYoQIggtPCWmVJBPIAUNEOuO3x3AvHhU+GoE088TI4RUQRFcbooqqJvcpZlEYfWVj05vKq69PqP31DifnnKXYwRR7uHLhs8/fS151dEPwt1C5gyY4ePgJXdRwFJAi1BBCKcTyq8zPkEIodXnxHOlZMMQGB1jqAGeeEVfZuoh4WI5QQixAZH63BS5izI+naDbcuH7l/csiPXq8eF9K20jBZv09+2z4f53g5wDn98xe8v573mRl7siYBEgCBWEy09nfYMxzWxMkXYVkY3DGHtoC3bEUB+8pDc9OVFIqyJsUCSXc1nuiozsxNXTk/6YUVh95cN733k44n6bvJ2motnx3cetHfnWxdLUF/6ak16WJXdFQH7wRqggXH6aKjDi31nSvzOlLSNoF2gLdtxIP0ovoAkHxF/u9bYnklhxg3bzlrsoI6jj6z87+6+TxWdn9X0+1ruH3OWYnK+j99Khc45dObUgccVDofdOi30GxpQqGbQIFYTLS72oCVt3Qfx4iNLXjumKBwOp6d2oCQdESWsjvaMnr559bvvLelH/wT0rlZCCTYZqB3xwz8rcqoIXdryWVpYpdzlANtAiVJDa3LQ3xb9vHMUEKWw1baN7IoQq16Hvz4RPyU6x7ztS7nI6r55v2HDm/5KvXXi57wsxXtFylyMDF7XTG4NmHSk8seDgikcjHnyu1wSGssEOYdA2aBEqxem0ogoOLxrhE+UCKWgE06Mpoo3MvJCiF+UupbPSy7Ke3zlbksQP7l2pzBRsMjxw8If3rkwry3x5z/zi2utylwPMDYJQEbKryaZdF6XA7r1a2k4BdM7E4RH+9YWT99SJRO5SOoggsjX19zcOLHu6+2Mv9nlWw6jlrkh+LmrnRUNm3x888qVd8/bkHZS7HGBW0DVq+0p06KFd4hY21SfUNie9yQUzLOsb5FGaOetYr8+HWU1/WqW+avXR9ZX6yjUj34S5dM2MDB4W5h6y/uTnSVfPzh040w4WKVUGaBHauFoePfSXMC6ICim7hAMgCI2M0kbOd0s/dp2sPi/JXYtBzl67MG37q8Guge8kLIYUbFGQk/bdhKUSEafvnJNbmS93OcAcIAhtmUjQxINitAt+PawG1VchL4vbLsDa4YBItiDlmwTm8xTp2yyLzkKCyPeXf155bN3s/i8+3f1vNjlH0FjUjHpW3+cfj3p49t4l+/IPy10OMDnoGrVZBKHph8V6gXw8hCGpl7C2W2u70oPO849Cf23x0ZBvR9FP7RW09vgef0v8JjcIunePry+qubo6YaksO+hao4SgocGuge8d33j22vk5A2bAaFIbBi1Cm/VOsnSujHw+jGEoRPIvo4AouSuyRU5uWKVBZVeinPHHQ5kJB4S0SosbOVNYfeWlv+YihFeNWAIp2CEhzoHr7lle2lA+e+/isoYKucsBpgJBaJu25kpfpEqb4+mbi6jlXkRwg9BEtFFSQSpCKN4Xv92XHrNLvN4gd0m3OVx4YtbuBQ+HPzC7/3QVLJ7Scfas3byBs3p6Rv/jr9dTSjPkLgeYBAShDTpVQl4+Kn4zkva1xwghxOlISRH2C5O7LhuljUQFqY3/fSyYeiwYP7xLqBPkremmf1/a+tGpL5YOff3ekHi5a7FiFMbju4/7R+/JixJX7Ms7JHc5wPggCG1NTg0Zu0dYP4Tp7vq/m1WFqdgnGDHQGjAJHBCF8m6tvj03lg53xs8myjy3kBf5d4+v35t/aNXIJeFuIbLWYiP6+/ZeNnz+pnPffJH8tUQsrgMcdAUEoU2p5tHY3eKrPemRfreGbJC8yygQ+kVNxiuQVF5HutrGjzBC/xxIX68nb56WbcmZan3N3P1vVetrV49Y4mnnLlcZtifIOWDtqLcu3EhZdmStTtDLXQ4wGghC28FL6PE9Qrwvfi7yjpdVyr2EtZFyVWX7KBppI0j+rdW3WQp9Ec/8J4fIMqEipzJ/+l9zoj2i5g6aoWZU5i/AtjmpHN8a9jqNqVf2LiytL5O7HGAcEIS2Y+YxkaXw0t53DvKWRHwlA4aMmhQO6E5yL97+iJsKbRlBzzshHr1u1j6008Xn5uxb8kyPJ8Z3H4eRJU7ksAEMzb7c74VBfv1m7H4jB2bc2wQIQhvx0SXpcDH5aChN3fnuR65mIxdPpHGUqS5lCI5GOReaPRbljD8ayjy+V8iuNlMW7srZv/Lo+/MGzhoeMMg8V1Syx6LGPBvz9Gt7l5y8elbuWkBXQRDagj1XyOpz4lcJtNPdCyTkXYIbhKaG/SPIjQLENZ82keCLZ/ekH9ktVnEmr+G7lJ82X/hu5YiFPTzh5TaTIdr+C4e8uvrE+j+zdsldC+gSCEKrl1ZJJh0U/i+eCXZsoSuM5F1CWugXNTGaxX6ht98mbDI1ihrijSfsF0w3ilQi0rqTn+7PP7xm5JtaJz9TXQa0pJt7xKr4xd9d/nnLhe/lrgV0HgShdavQo0f3iPNj6QFerdwQKkjB0CI0g6Ae5K7e0UYr+tEcQQtPmWQQqU7QzT+wvLi2eOXwhS5qZ1NcArTN19F75YhFRwpP/DPpE4lY9HqzoDUQhFZMkNCTe4UHtPiZ8FZex9IrhGGRM6yqZXI4MBrlnG/xUzRGHw9htuWRzelGfpes1te8tm+Jq8Z5weBX1bCnoHzcNC4rRyy6UlP81uE1vMjLXQ7oMAhCKzb7hEhTaFFcq2sBkzzYeslctJHkej7idC1+0kWF/jWCXnRKPHHDaD2kZQ0Vs/cu7u4e9VKf52ArCdlpGPXCIa8iQubuf7uOr5e7HNAxEITWalOatO8K+XgYQ7c+SJ7kXcJwg9A8GBb7hJCCFm4TNop0xu8PZp7YJxbVGSELC6qLXto1Lz5wyKSYp2CahIVgKea1gS8FuwS+tndJha5K7nJAB0AQWqVj18lbZ8R/jWhpmOhtSN4lAjcIzSYoGuVcbOPz9/njaVHUo7vF+q6tRJpeljV77+Knuo0dFzm6SycCxoYRfq7X+EH+/WbsmnelpljucoChIAitT0EteXKf+NEQJsSpraYAqa1ADTXYM8BshSkcDupBclu+TdhkZncq0hlPO9z5MaRnrp2ff3DZzL5TYR1ti/VY1JixkQ/O3rs4t6pA7lqAQSAIrUy9gMbtFl/qTo3wba9DLPs8DoyGzXjNRxtJinMR184SlO8NpDMrydrznRk4c7To5Iqj/3xj0Mt9feI6VSIwkwdD75nc86k5+5aml2XJXQtoHwShNSEITTssRrrg6d3af+FI9jkU3NMMVYGbWBXyDiJFqW0fpabR5hHMx5elPws61izcl3/4vaQNCwe/1t0D7vtageGBg2f2mTrvwNtnrrXTTwBkB0FoTVYnS1lV5L2Bhg0RzE5GITEmrgjcAQd1b/s2YSMfO/T5cPr5w0JGlaFZ+EfWrk9Of/n28PlR7rCvpNXo5xv3xqBZK46ug2XYLBwEodX4s4B8mip9GU+rDcnBihuI02MPf5OXBW4XGN3ubcJG/Tzxglh67B6DVl/7Ke2Pby/9uHLEomBnuONrZXp6Ri8ZMufd4+sPFhyVuxbQKghC65BeRZ4/LGwaTvvYGXTPj2Qno9AYuEFoZjggilzJRrxBS4tOCKeGeeO/72tn9bUtF77/NXPHqhGLfR29jVMlMK9wt5Alw17fcGrTXzn75a4FtAyC0ArU8OixPeL8WLqvp6HBRrLOouAeJq0KtEClwd4BpDDdwMOX9aXrBfTWmVZXX9t8/rtDhcdXjljgbudmpBKBDEJdgpbHL/jy3Lfb0rfLXQtoAQShpZMIenq/MNSn9XXU7kYIyjkPI2XkEdjdwN5RhBBDoc+HM//OIltzmw8iJYhsPPN/x4pOLhv+hrMKFhG1ev5Ovu8kLPoh9Zcf036TuxbQHAShpVt6WizTo7f6dGANLXIjH6nU2MXLdFWBVgVFo9z2x8s0cVejLSPoGUfFs6W3ekgJIutPfZFSkr4sfr6TCvaStBFe9p7L4xf8kr79q4v/lbsWcAcIQou2LU/6Lpt8Gc+wHXqhss6hYBgvKg8cEE2KMpDQgZWXe7jitQPox/eJpXqMEJKItPb4hszy7CVDX7dn7UxWKZCBl73HOyMW7c8//EXy13LXAm6BILRclyrxi0fEzfG0Rwf3FSDZyXCDUDYqDfb0J0UZHfqiMYHUIwF4yhFKL4orjq67Xlfy5rB5dqzGRDUCGblqXJYNn3/sysn/O/eN3LWAmyAILVS5Hk08qlrZl+7p1sGRn5JE8i9huEEoo+CeJKvD88YW9aZpJE7Ysa6Wq1s4+BUVzZqiNGAJXNROy+MXHLtyatN5yEKLAEFoiUSCntkvPOAnjgvp8AtErmQgZ09kD8MrZINDeqGOB6FEhO7q96/X1wd7zWAgBW2dk8pxefyC86WXPzqziSCjbc4FOgeC0BLNSxJ5gl6L7tSe5rCymuyCuqHr+aihxvCvECRh/ck1NJYWDHr53Qv4fDm8M9o+B9Z+0aDZ6RXZ75/8TCLwissJgtDifJsl/V5APhna1kaDbZCyk6kQCEJZUQwK7EayDZ1EIYj8ByfeJQRN6vF8pDO9tDc97ZB4o+UtfoFNsWfs3hwyJ7cyf93JTyALZQRBaFmSy8jcJHFzPO2i6tTX8xwqzCAB0UYuC3RUaCzJPGPIgZzIrTmxgqXZab1fojGFEBrpR40Lws8fErhO9QgA66Km1YsGz86tLIAslBEEoQW53oDG7RHfG0B3c+nk0mikMAV5B2I1jLmXGQ6LJekn2z1ML+jXHl/uqHKaEvsChW/9Mr7YnXZh8ZLWV5wBtsSO1bw5dE5+VeF7SRshC2UBQWgpOAk9vlcYH0o9GND5F4VknYPxohbB3Q/TNCotauMQvaB/78QKB9Zxcsy021Ow0cr+dNIN8nVmZ7YtBFZHzaiXDJlztaZ45dF/igT+ADI3CEJL8fIx0YXFc3p17RXJSsaw9ZKFCI4hGa32jjamoIva9dnYf9ydggghewatH0x/eEk8dh2aCIqgZlSLhrxWrqtcefR9yEIzgyC0CB+nSInFZP3Qzo2P+R+ugdwoQP4RxqoKdElorJRxqsXP6AXd2uPLXdVuE2OmUa3vEKJ1wCv6Mi8dFQpqIQsVQUWziwfPruXqVhxdB1loThCE8kssJu+cFb9OoJ2Yrp0o5zzSRiAGpqBZBBzSE+en3L3WWoPQsOroW572XhN7TW0jBRsN8sZTIunpR0SdYKo6gUVhaHbewJlV+pp3oF1oRhCEMsuvJRMOiBuHMcGOXd07kKSfwqG9jFIVMAKNA/EKQvmXb3+sXqhfffQtP0f/p3tOwcigV3xyBBXhjF8+0fauhcB2sDS7cNArVfoaaBeaDQShnBoE9PgecWYPariPEXbQldJP4oi+XT8PMBYcFiPdNomiXqhffeQtXwe/8T0nG5iCjRbG0kW1ZGMKDJxRCpZm5w9+pZ6vX3bkPUGCLDQ5CELZEISmJIpRrvj5KGO8CtdyMaaRh78RTgWMBIfEoozTjf9vTEF/J+2EmGc7lIIIIRWN3h9E/ytd2nMFmoVKwVLM3IGzOIFbdmQtZKGpdewtuLy8PDs7W5Ja/cv0xo0bBQUFoggvW/uWn5Xya8ia/h3YaLANJP0kiuhtlFMBo/ELI9WlpKa8XqhffeRNrVOA4T2izXhq8PrB9OtJQloVZKFSsBQzZ+AMXuSXHHqHFzuwsRfoqA4E4dtvvx0REfHoo4927949Kyur2WcvXrwYFRXVo0ePUaNGBQUF7d6926h12ppf8qT/S5M2xdMq4+QgklKTUBgEoYWhKBzcsy7j+MojS0JcIzqdgo2iXfHrMfSziWK53oglAovWmIWIoLePrOUlGDFlKoYG4YULFzZs2HD+/PmUlJSxY8e+8cYbzQ7QaDRbtmwpLS3Nzs5evHjxpEmTCCyR0IpzZeTFo+JXCbSPnRFuDSKEUEMNup6Hg7sb52zAeGpDur2T/+8I16gnosd3/WxjAqlRvnjGEUGA24WKwVLM6wNnSpL05qF3oV1oIoYG4X/+859HHnkkMDAQITRjxow///yzpuaOxfUjIyOHDx/e+P/Ro0eXlpbW1dUZt1bbcKMB/W2vuKofHdPRjQZbR9JP4ZCeCPbusTB1Qt27DYcjKnSPd/u7sc45O4YmCL19Fu4+KAhD0XMHzVRR7KLEdziRk7scG2RoEObl5YWHhzf+PzQ0FCFUVNTq8lHffPNNfHy8o6Nji58VRbG8vPzMbfR6pfT18BJ6ar/wZAh+JMiYw5RI2gkU3seIJwRdV81Xrzj9bqRbt8frXNGNfGOdlsJo7UDmyDXyVQa0ChWEwtQr/V9wYO0WJ67SQxYam6FTuGtra+3sbi7ljDHWaDTNWoRN/vjjj88+++zIkSOtnSo/Pz8pKWn69OlNj6xbt27gwIGtHV9XV4fbm3dsLWYmMfYYvxgmGNJarq+vN+SJY0lUZSbzA8aS+nojlGgBGnQ6IklW/aJX8zX/vPxRtHPkaN/7eL96lHZKdPJu96sadDqBZWm6nfvGGKFVcXhmEhOg0g/xsoU41Ov1DMO0+8RtD8dxGGOWNbQvZ2r0hE2X/v36nqUrhixQ02qT1mZSPM9LksTz5ujptbe3p6h2Gh6GBqG3t3dlZWXj/zmOq6ur8/HxufuwPXv2vPDCC3/++WdUVFRrpwoLC3vooYe+++47Ay9NCGmtcWld/nlBOlcp/XI/Y88Y9BNMCHFwcGj/sNyLxNVL421DEycwttNokNUGYRVXs+7Chp7uPR4LfQQhRMJj0amdqoSn2v1CTFEqlqUMyINIe7RmIJl7Gv12PxPiZK3fqCaNKcgwXVxayfqwLNuhIEQIzRn04idntyw7+c/3Ri2zYzSmq82kGoNQrbaULDe0g653794nTpxo/H9SUpKXl5dWq212zJEjRyZNmrR169Y2mneKtbOQrL8k/SuBtjf6L3taEoJ59BajQl+5/PSqvh69G1MQIYQDIkl5MaouN+6F+nnhl3vQkxPFahg/oSQUpmb1neZp57Hg4AqdANs3G4ehQThp0qQLFy6sX7/++PHjc+fOnTFjRuOfby+99NJHH32EELp48eLo0aOfeuqp69evb926devWra31nSrQ+XLy7CFh8wja3974f7yT9CQUDhMnLEKZvnzFmXf7e/d9KOj+W49SDA7pSbIM2qe3Q8YGU4O98YswiFRhKEzN6jdN6+g7d/9bdbyN3BCRl6FB6Obmtm/fvmPHji1YsOCxxx5bunRp4+PR0dGNQ0k5jps4cSLP83v/p95Wbll10bUGNHa3uKofHedughQsv0bqqrFvqNHPDDqqVFe28sya4b5DxwQ+0PxzYb1Ress7UXTR6zG0KKE3YRCpwmCEn4+bGOIS9NreJdUcNDm6Cpt/tt/333+/fft2w+8R1tTUODk5mbQk09GJaOR2YZQf9WrPDg8Tra2tbffmKDn2GylMw2Omt32YdWloaLC6e4QlutJ3zqwdqY0f5TeihU8LnLR5ETVrA9K0ddNXp9MZeI/wdvUCmnZIeCaMmt7dWldM1Ov1yrxH2NHBMs0QRL6+8N+UsowP733HWW1Nb5LWeo8QdAJB6NlEMcABv9LxFDT0EmnQLyq/4rpry06tukc7ouUURAgxKhwQhbKTTXH1xi18P02T9lyBHlJlwQg/Fzuhj0/sa/uWVuqr5C7HikEQmtDS02JuDVk3sGvb7baB05GCVBzS00SnB4YorLuy8uzaR4IfSvCLb+MwHB5H0k6aqAZfe/z+IPr1JPFSBSznpDgTejzWx6fX7L1Lyhsq5K7FWkEQmsqWDOn7LLI5nlGbbHIUyUpG2giktjfVBUB78mryV519b1zImCE+7Y2UDu2F8i4j3lRToXu64cVx9JRE8Wo9ZKHiPNPj8fjAwTN3zy+uvS53LVYJgtAkdl8hi0+K/x5Je5pyng9JS8KwoIx8cqrz1pz74OnwJwZ5D2j/aLUD8QlGuRdNV889/tTTYdRziWIdTKhQnscixzwa8eDsvYuv1BTLXYv1gSA0vovlZNJB4YvhTJizKYd7SCJJPY6j+pvwEqB1XdzDZwAAIABJREFUaRXpa5Lffyb8qd4esQZ+CQ6LkzJM1Tva6NlIqpsLfukoTKhQotFh9zzR7eFX9izMrSqQuxYrA0FoZMX16JHd4lt96IHeph30SLLOIjc/5OJp0quAFqVUpH148ZNp0ZNjPWIM/yocHoeykpGJN1ld3JsWCZ5/CiZUKNF9ISOnxPx97r43sypy5a7FmkAQGlMtjx7eLUyMoJ4IMfk3llw4hLsPMvVVwN1O3Tj70cVPX+wxLdq11XUEW+bgil08UWGaaeq6icJodT/6fBlZfxlahUoUHzjkuV4T5h14O60sU+5arAYEodHwEnpyrxDjil/pYfrvqiiQtBO4GyxlZ25Hrx3/Mu1fM3o8H+YU0pmvD4sjpplZfzs7Fm0YwnyXJf2UB1moRMMCBs7qO/WNA8vOXDsvdy3WAYLQOAhC/zgiEoxW9TfHCvok/RTyDkJObma4Fmiyt2j/vzN/mB0zM8QpuHNnwOF9SPpJZPpVLDw1aOMQesUZ8dA1GESqRH194uYPenn5kX8eLjwudy1WAILQOJacEi+Xk8+GMoxZvqPkwkEcPdgcVwL/83ve9t/yt8+NfcXfwa/zZ3H1RioNumaO+zehTnjtQGbWMSG1ErJQiXp4dls67PX3T372V85+uWuxdBCERvDxZemnPPLVSMbOPEtEcXqScRp3M2DIPjCSn3J+TSw+PLfXq54ajy6eCpuld7RRX088vxc98aBYWAtZqEThriHLh8//8ty3P6f/IXctFg2CsKt+z5fWnJf+nUC7qcx0RZJxEmsjkL2zma6nbBKR/i/lX8mlya/HvuKqdjHCGSN6k/QkI5zHMA8EUM9HUeP3iyWwY48iBTj7rxix4MfU375L+UnuWiwXBGGXHLpGph8Wvx5JBzqab4Vocv4ggn5RsxCIsPHSZ1cbil+JmenAtL9JsiGwVzBCCBVnG+VshngylHpQiycfEGphor0i+Th4r05Ysicn8aPT/yeZfZcFqwBB2HlnS8mT+4RPhjM9XM24T4K+AWUl46h+5ruiUulE/XvJH+pEbmaP6RramMvk424D0YVDRjxhu17qQUe74mmJAgfTCxXJVeOyYsSCtLKM1cc/EEw8k9UaQRB2UmYVeWS3uKofPczEE+ebIWnHSVB3pGlneybQRVVc9fLTqz007tOjn2OpTu6S0xrcbaCUehwJZm2gLYyjHVT4xaOiCE0CRXJg7ZcNm1fL1S9JfEcn6OUux7JAEHZGQS25b6e4II56JMjc30DpXCLMoze1Gw0ly06t6uXec3z4E5QptkV0ckeeWpRt1jleGKOVfelKPVl8GqJQoRiafX3gS/as/dz9b9ZwtXKXY0EgCDvsRgO6f6f4Qjfq76Fm/+7palHeRRQB/aImlFOTt+zUqnu0CWOC7tpo3nhw9CByMdF0528RS6MPBjPnSsk7ydA5plA0pmf1mxbpFjZr9/zrdTfkLsdSQBB2TDWPHtoljA2iXugmw7eOXD6GQ3thtZ35L60Ql8pT1px9/8nwv8X7DTXphXBEH1KYhuqrTXqVu9kz6NNhzMFr5P1LsOiMQmGEJ8f8/f6QkTN2z88oN9+gLUsGQdgBdQIas0vo64Ff7yXP941cSMTR0C9qKoeKj268+NmLPab29ext8osxahwag1JkWPXDkUUfD2F+zpE+T4UsVK4x4fe9EDtx3v63T149K3ct8oMgNFSDgB7dJQQ64GV9zbGIWguqSsiVDBTRV56r27rf87dvzdk2J+6VcOcw81wRdxtMLpp17GgTdzXaNJz+V4b0TSZkoXIN8u+3aMjs1SfW/5m1W+5aZAZBaBBOQk/sFbzs8LqBNGXWUaK3kFN/4R5DEWuuefuKIRJxU+q/jhUnzY97zdfO23wXDogidVWopNB8V7yNlx3eNJzecFn6bw5koXJFuYcvGzb/64v//fbyVrlrkRMEYfv0IvrbbsGOwe8Pki0FkSRKp3bi2ASZLm+zGoSGNcnvl+lK58S+7MQ6mfXaGONuA8ilI2a96G187fHHQ+k158TfYJMKBQt09l89cunB/CPvndio2CmGEITt4CX09/0ihfH6wTQtVwoiRNJOIhdv5BMiWwW2qFxfvuz0u24q1xe7P6+mZWhq4+5DyKXDSJIth0Kd8BfDmWXJ0s+QhQrmrnFdGb+wrKF89t7FVXpzD+CyBBCEbeEl9PQBkRfRp8No82wr0RpycjvufY+cFdicrKqcpUkrBvsMeCby7xSW6dV19UaObij/kjxXRwghFOKEPxlKrzwrQbtQydSM+o3BL/f0jH5p17zC6ityl2NuEISt4iT01D6xjkefD6dZeb9PlTdIUTqOhm14jSbpxql/nv9wYtT4Uf7x8laCowfJNWSmSbgz/nQY/TZkobJhhMd3Hzc24sFX9y66WJIqdzlmZZ59g6yPXkRP7Rf1AtokewoiRE7uxD3jYZiMURBEfs759cDVw6/GzNA6+MtdDsJR/aTjv2N9HVIbZ1Hvzgl3xl8Mp186KgoIPREi9088kM8DoaO87b2WHlr1Up+pD4XdK3c5ZgI/8S2oF9AjuwQGoy/jaZVMcyVukUTpzC7Ue6TcddgCvahff+Hjc2UXFsS9ZgkpiBBCagccGkPOHZC7jlt9pFtzoV2oaL19YlaNWPLd5Z/+mfSxQobPQBA2VyegR3cJzir80WCZ7ws2IinHsbsv9gyQuxCrV64vX376XYTw7JiZzioL2s0R97mXnNmFLOAdJ9wZb4qn15yXNsFce2Xzc/R5N2Hp9bqSBQeXK2FVUgt4p7ckVRx6YKfg74A3DLGIFEQIkaTtCIbJdFl6ZebSpBW9vXpN7TbJ6LtJdJVnIHZwR+mn5a4DIYSCHfGXw+ktmdLGFMhCRbNn7RYMfiXISfuPv17Pq5JntqvZWMabvWW43oBG/in0cMVr5Zs13wwpLybXsnHUALkLsW67i/Z9cGHjlG7PPKC11Hseve+RTm6Xu4ib/B3w5nhma4707jn5G6lARhSmJsU89UTUo7P3LjpSlCR3OSYEQXhTXg0Z/qcw3Bev6CfjdMHmyMkduNcIxFhYC8Z68BL/RcrmvUX734ibHe0aJXc5rcJhsbiuChfnyF3ITZ4atGUEc+gaWXBKlGDTJmUbFTxs6dC5G05v+iL5a4nYZj8BBCFCCF0oJ8P+FF+Iohb3ln1szG0kAZ3dg+OgX7STrtffWHpyeT3fMC92tqfGQ+5y2oQxih1Jn90ldx23OLHok2FMSgV55bjI2+a7HzBUmGvwe6PeTi3LeG3f0gpdpdzlGB8EIUq6QR78S3irNzUl0rK+G8yFg8gnFLn7yl2IVUouPf/mqZX9vPpMjZ6kkmPVmI7CPYfiglRUUyZ3Ibc4MOizoUy9gCbsE6p5uasBsnJSOS4ZOifaPfIff72eUpohdzlGZllv/eb3W7706B7ho8HM2GAL+1ZIEnP8VzRkrNx1WB+RiP/N/unL1K9m9Hz+gQBLvSl4N1ZNIvujs3vkruMOLI1W96cDnfBje4TrDdBJqmgUpsZ3Hze114RFiSt+y9wpdznGZGHv/ua1/qI046j0TQIT72s5twVvIuf3E0d3HNhN7kKsTLm+fOWZtVmV2Qv7zA11CpG7nI4RY+LRhUTE6eQu5A4URgvj6Pu11KN7xMxqyEKlG+Tf7534xb9m7Hj7yNo6vl7ucoxDoUEoEvTKcXFTuvTb/XScu8WlICKEJP4oDHhI7jqszKXylKUnV0Y4h8/s+Q8n1lHucjqMOLkj/0hy8bDchbRgahT1UjT1xN7/b+++w6K41j+Av2dmti+7lKUvHaUoqFiwoyhRUzXW2FCjuak3XpOoMSFNb0xMzE0xzdQfxpbExwZqRAUlwShg7yCICkpddlm2zsz5/UHi9SYWVNjZZc/n8Y9dnAe/8+w675wzp7AHa0gtdHdBHgHvpLyqYBQzs587VXdW6DhtwB2XWDPaYVIua7TDxjRG5ZTjMfHxfBBLeG2s0EFcBoe5jeVb86r2PR47zWE767aLbkNw7hqUNByQ092fPRBCaSRodj6bkUSPj3DTe2iiBUOLHu82uatv7Mt7F0+IHT25y1gETveNbT23+zZfMUFKNqsSocwUJ62CgDGftxb1HyN0Dpdx1VT9euGSEn3py91fdO0qCACBkUgig3PFQue4sWQ/9O0g5sMTfEYxmVZBQHJQ0tKUjLyLv72ct0Rn0Qsd5+65VyE8UIP7bGZHaqnlycIvpX0z+PTvABgiuwkdxDXsrszLKFzc26/nU/GzFSK50HHaAOpzP5+3xhlWXLuhEA/0XQpzvAHPzOeayVBSt+ev8F0yeFGQ0n/Wtuf2VzrF6kh3wVmrQTtYXco/lMO+1ZN6Nt6pzxrnrkH9H3XCnjFnY7Ab3j/60baLO+cmPJ0SOFDoOG0nrCtSqOGEMz4pbKEWwWf9GR8xPJTDXjKShqG7Yyh6SpdxLyU/+2HhF0sK/mNhnWu0V2s4dUloK3Yeni3g3jjEbxjGjNQ69SnjkmKwmVHnnkIHcXZH6o8u3P+aRuLzco8XnGUfibaD+j/C7/0R7Dahg9wUTcGCbvTDYdRDOeyv1aQWEhDjHf3+sDftnHX29rln60uFjnNnOv5gmToLTNjN0hRkjXDWh4LXwXvWQP/RpDl4Cya7aVXJ2uMNJ2fETu2sjhY6TvvwC0f+4bj4F9T3IaGj3MpjkVSMCj1XwE2JRvO6OssKvYRQ5Izsn72e+PXygfm5bzzUaWR614ki2umvuQDQ4VuEBdU4aRPbzQd9P9gVqmDZEdzUgGL7Ch3EeR2tOz7/wKscz72atKDDVkEAAED9HsYHssDaLHSQ20jSoFUpdG4VnpnPkdVnCAAYqE3+z/Alpbry2dvmnq53jTVoOmwhxAAfHucfyWEzulPzE13hXpXn8NbPqSETgOqwH8q9MLPmr0599/XZ76d2mjS50wQpLRE6UTvzCkARiXh/ltA5bs9Xhr4ayGjlMGwbe6SedJMSoJaoXkp+Zlzcwy/vXfJx0VcW1ip0otvomNfcOgvcv4NdU8bvGMk8GOoa54j3b8EyFcQkCx3EGR2oKXxh/yIe4Ywe8515E4m2hZIfwEd2g8GJVh+9GZqCf3ah/xlPTdvLZpaQJboJAID+wb0/SF18xVg9Z8e/jtacFDrOrXTAZ4T5V/HkXO6RMLQy0Vk2170tbNTxeWupKa8JHcTp1FrqvjuTWWOqdfmZ8ndB4Yni++PfNqJRs4WO0ir3aakYT/RyIZdbhT/oR3u5wFLnRPtSSzxeTH7696riN399v09Q96d7zFJJPIQOdQMuUihax87Da8XchN3sO73pV7q7TBUEALz9G5SQAj4dbfTjveAwt/1izisH3gySBy/s8YLbVUEAAEA978PniqCuUuggrRWmRJkpTIAMhm0jo0mJP/QN6rkibamMlk3d+tTW0l8wON0Xo+O0CE814ml5nEYKO0aJ/KRCp7kjFSfx+cPUnPeEzuFEzjSe++ZspqdYPb/7XGffSrBdSeSoZxrelYkmLnSVscQMBXMT6H7++LkCbmwEWpDovItXEA4jYcRTu47rp+315eH/21W+97mkOSFKJ7rv7wjfUB7Df47zKVnsY5HU94MZF6uCPM9t+RQNfQzErpW7vdRbG744+/WKk1+MCkl7pssTbl0FAQAAJQ4BcxMcyxM6yJ1J9kNrU5lTOhi5gz2hc7oWACGIKM/wd4Zk9Azs/kJuxqdHvmuyGYVO9AeXL4QXmvDwbez6cn7rfcyUaNc7Hf5gNmIkKK6f0EGEZ+PtWyqyFx143YPxyOixMEnTXehEzoGiYfg0Pm89GOqEjnJnvMTwn770lGjqsT3ssmNkm3sCAIBC1ANRw5cPfcvG2aZsefLnM1t5LPw3w/UqxzUchuXH+F6b2YEB1M/DmDCla3Qc/Q+TAfasRvelu0qvVzvBgPdXH3ihYOG5xtIF3f/1QMhIiSvsKe8wyCsQdRvCZ38J2PWaVg+EUOtSRUfr4b7t7NEG18tPtAelSJHeZcJrA1/Mvfjr7G1zD109JmweV31GeLgez97HKUSwNY0J93DVKsJv+RTiB4BfqNBBhHSy4fSa0h85zE3r/FjLHHmzxfXWKmxvKCkN/3QEju+FxCFCZ7ljGil80JfefpmflstOiKTmJdByV73wEG0pQh361qAFBZWFyw58ovUIerLHjGivCEGSuN730cTC68VcZim/qDs9PsIFJsrfDC7aAVfKUPpioYMI5lJz5c/nN5YaykaFpPX370u5d7P4NigapaXzGz+iwhNA5ZLPTUdpqWRf6qMTXEo2+1YSNSrEhbujiDbUP7h3n8Ck3Iv58/e80dUv7skeM4KUAQ7O4GKFcH0Z/9JBPtkX7b5f5OPKS4vgmgr+l2/R5AwQuWMf4FVT9YayTccbTo4MSZveeTKNyDXx9pB3ICSm8NlfUJMWuWhfurcE3uxJH6rD7xzlfijl/92Ldt3uHKINMRSdFj5kkLbf5pIdT+yYNyJi6JT4cd4yL4cFcJkL0OF6PDiL/fdh/uO+9Mf9aJeugmCz8KuXUEOnIE2w0FEc7aqp+rOTX71etMRb4v1mr1eHBA0iVbD1UM80MBng+D6hg9yTJA1aM5RJ0lAP7mTfP8aZWKEDEc5Bykgmxj3y8fC3baxtevYzK4q/bjDrHPNPu8A1qM6KnviVG7WDfTiU2j6SSfZz+VtIfvMKFBgJCYOEDuJQtZa6r05993rREpVI9UbPRaNC08iImDtGMdTwdD53LVRXCB3lnjAUTI2m1g1lzhpgQBa7rhyR/e6JFmqJakbiY5+kvc3x3NSsp947sKLO3NDe/6gLFMKXDjEUwN4HRVOiKRd+JPgnXLwTV5xEaelCB3GcsqYLHx3/7NWDbypFyjd6vfJg2AgZIxM6lMvyCUIpE/mflkGTg26W24+vDC3pSX+QTP9cge77BedeIcWQ+INKrJradfwnw5dSgGZkPbv84OeVTVfa759D2OEDstesWZOdnb169epWHl/b2MSKnXF5uruAayr4lS9SkxaBf9htDzY1N8sVCgekaj9nG89tvpBVZrgwMKD/MO0QGd2qRQPMFotMInHRx2D3wmK1ikUMRdG3PRIX7YDy42jq6x3jGbPNZvutBn1yBmkVsLAb3cPHXT56m82GEBKJnH6LuLZmt9sxxmJxq769Bpth+/ncHWW74306T0+YGK+JafM8LjZYxrU16/nVi9HQya2pgi6N5e0F1QezK37hMJumHfp47AzyILBtoZ4jsO4q3vYFevi5jnHHMDgADQ5iNl3k5+RzsWp4MZHu7jblkLg1lVg1Me6R0Z1G5lzY93r+Mq0q8O2UV2VMWy7FRQqhg2BrM/7uFRSTjBJThM7SjuqtDTmXdudW5Ycogx8Juz/OOwYBuZy1A4RQ6hR+40fw6wY0aJzQadoGTcHYcOqRUGrHZf7J37hQBbzc3Y1ah8StSRjJg9Fpo6JST9eXiNt643tSCB3CZsXfvwbB0WjQWKGjtAsM+FTDmZ2Xd53UnUn26/VC4nN+Ml+hQ3V0tIh68En803vgEwjxA4RO02YYCh4MpUZqqY0V/Ox8LtoDnu5CDw7oEM1e4p7RiO7p341Gt398cEdIIWx/dhufmQHeAWjYNKGjtD2dRZd3JT+vKl9CSwYG9JscPZGMBXUcqRLd/w9+08eUSAqdegqdpi0xFIyPoMaEUblX+H8f4d7kYFYMNT6CkrTxBZAgAEghbG+YY/HaJSCWo/tmdoxnOS3svP1w3bHcqr2ljeeTfLvPip0epgwROpRb8glCDz+Nt36OTE3QbYjQadoYQ0FaMJUWTO2vxqtKueUn+JmdqMlRlIbs1EK0KVII2xPH8uveQRij0c8A1RFGi/AYn9ad+fVqQVHNIa1HcF+/3jNjpokptxvz5lSQbyg8OpffsoKyNkOfB4SO0y76+aN+/kyJAa8/zw/KsqcEUunRVD//jnNnSQiLFML2gpt0/Oq3kFyFHnkWWjEg3pnxGJ83nP+9pmj/ld9VEo9evkmvJM33lKiFzkX8ydOfGjuP37QCmY0oZaLQadpLJxV6tQf9rwQ66yK/sJDDANOiqTHhpIFI3CtSCNsFvnyW/2Ex6j4U9R/tuj2iHOZO684eqCkqrjkkF8m7+ST+M+HJALmj18MlWkXhRY39F976GTYb0X0zO0YPxA0pGJgYSU2MpA7X4U0V/PIT9t4+aFwENUJLScn1jLgr5IvT9vhDOXj719T9syHaJccvGFnjsfqTh2uPHKk/rpH69NAkPp/wtL/cT+hcxO1IlWj083jH13j1YvTw06Du4AN3e2hQDw1tYek9V/jV5/kFRdwoLfVQKDUoAIk67G0A0S5IIWxLmGPx9q/wmYPUlAzwCRI6zh3gMb5ovHi0/vih2qOXmi93Vkd38Yq7P3SEl8RT6GjEnRBJ0INP48O7+O9fpYZNg64DhQ7U7qQM3B9C3R9C1Vngl8v8B8e5ZwrwsEDqwTA0JIC0EYlWIV+TNoPLj+PNn4CnPzVjMUjkQsdplaum6hMNJ483nDqlO+MhVsaqY0aEpHZSR4vI+BfXhRBKSkOhcTjn/6C0GI2cDVLXXqivlTRSmBJNTYmm6iywp4r/8jT//H5ugD81LAilBqFAuas+oSAcgBTCttDcyG//Bp8/goZNRTG9hU5zKzzGlc2VpxvPnmk8d0Z3DgBiPTvFe8WMiXjIU0wGv3QgGi2aMB8XbOa/XkCNmNnBZhnemkYKEyKpCZGU3g4F1XzeVfz2US5AilKDUUog1csHkWYi8RfkG3FvMOYPZuNdq1DXQdScZSByxm0SjXbjeX35eUNZif78OX2JSqyK9AjvpI4apU3zlWmETke0G1qEBo1D4Qk4dx3ev4VKmQRhcUJncii1CEZpqVFawJg+ocO/VfPvHOHO6HGiNxrgRw0IQEk+SOzaA7qJtkEK4d1i7fjIbj5/A5J5oMcWIV8nmk6utzVdbLp4wVhRbqg4byhrsjWFq8LClCHJfr0nd5rgIVIKHZBwoJAY9NgiOFeIt30BXgFoyCQIiBA6k6MhBAneKMGbfjIOzCwcqceFdfwbh/gSPY7zRL00qJcv6qVBAaT71F2RQnjnTAb+9yz4fQsOiKDSpkNYF2Hj2Hh7ZXPVZWPlJePlS8bLFU0XbbxNqwzWKoI7q6NHhKT6Sv0ol53CQbQBhFBMH+iUBCcL+J/eQ4FRqMcwiEjswFMsbkHGtEzPpwHAzMIpHT6qw2vL8MuFnIhGCV4o0RslekNXL/JY0Y2QQthqPA/lx/hje/GJfahzb3hsEaXROjgCBlxvabhqqr5qunrFXFPVXFXVfKXR2ugvDwiQ+wXJAvr69xkfNdpb4u3gYIQLoBhIGEzF98NnD/L5P8P2r6jEFEhMAU9/oZMJRsZAT1/U07el4NGXm/GZRnxGj786g083YkAQq0KxXihWjWLUqLMnUpExZB3UnRVCu91+2z0kW3OMC8EcC2VH8Ylf8anfkNoXYvqg2e8jZbuPK9HbDJXNlc3NplpzbbWltsZcW2uurTHXKRi5n9zPV6rxk/km+/YOCg/QSH0ostsf0Uq0CMUPQPEDcMMV/lQB/N/ryCcYOiehyG7g8Bs7Z6NVIK0CDQ/+422tGZ9vglI9LqjBmaV8mQHLGYhSoQgPFKVCEUqIUKEQOVJ0nKud+2ptITQajenp6Tk5OTRNL1iwYOHChX8/ZtmyZUuXLuU4btiwYZmZmR4eLrutPGuHynO44iRffgIungSfINS5D0pfjNp0hrKJNeuteoPdoLM26qx6nbVBZ9XXW+sbLI0N1gYpLfUUq3yk3j5SHx+pd7Qqwkeq0Ui8xWRvB6ItIO9ANHAs9H8ELp7CFadwUQ7m7SiyG4roBsHRoPIROqDwfGXIVwZ9/f7bQVptxhVGuNSMLzbhX6/CJSNfacIyBkLkKESJQhQQJEeBcgiQo2A5+MoQTfpWXURrC+Hbb79tNBrr6uqqqqqSk5MHDx7cv3//6w84cODAsmXLioqKQkJCxowZs2TJknfffbcdArcPcxOuroCai7j6AlSdx1fKkE8gaGNQfH8YMRMp7qD9hwE3s83NdlMzazLZTUa70Wg3NrX8sRmNdqPBbmi06pvsRgpRarFKJVZ5ilUqsUotVkerI3qJe3iK1V5SLzElMptMMrlrzEckXBXFQHgiCk8EANRYgy+e5o/moZ3fYQDkH44CIyEgEny1oNYATR6jgL8M+cugj+//1LcGK1SZcFUzvmKGM3q87yrUWPirJtDZsJcEfMTIXw6+UuQrBY0UeUnAW4K8JeAtBm8pkgl1JsT/Qhjj1hwXHBz8/fffp6WlAcDcuXNNJtPKlSuvP+Cpp56iaXrFihUAkJubO2nSpOrq6hv+qjVr1mRnZ69evbqVEWsbm1hxWzQurWZsMkBzIzTpQFeNDXXQWIMba1FDFWbt4KtFGi3yDrT5Btv9Qyw0sDxrYc0cxibWZOftNt5mYs0sbzdzFrPdZOasNt7WbDdZWIuJNVs4s4k1m1mLmTMrGLmMkSlEcjkjl9MyuUihFCkUjFzBKBQiuZJRqiUqpUh5200b3LYQmi0WmUTiumu03jWL1SoWMZQzLNHe3IhrLkFNBdReBt0VbGxECjX28kNegUjtC0pPUKhA6Q1yNShU9/5J2Ww2mqJopkPVWh5DgxV0NlxjhgYrrrdAgw3rrdBow402aLRBoxUb7KASgacYqcWgFoNKgjwYUIpAJQKFCKlEoGBAxoCCQR4ikNAgZ0DOIDEFahfvGLLb7RhjsfhuTkNCS5TiNl4jolXfPJPJVFVVFR8f3/I2Pj5+3bp1fzmmtLR0zJgx1w6oqakxGAwqleqGv9Bms+l0umtvPT090c3/L5VUHa5uqml5zdnMFt4KAJjngbUBgJmz8KwNMAbOxvK8xW4CnsecneWsVs4GHAc8a+FtHGcHhMxihqdoTIvMIgCKsSvyKwLFAAAKOElEQVSQTYkhVGbiEIZqC1fBGXiRUSSuEIlpCYNoKSOlECVnZCKKESGRlJExFCOhJDJaohKpJDQjY+RSSiplpFJaImWkMloqY8hNHuH6FJ4owhMiElreIZ6DpgbQ14K+Dht1uKYCmZuwsRGZDNjUBGIZksqxRIFkCpDIQCQFsRTEUkSLQCIFWgSMGABAqgQAYJj/Trdl/vwrmw1oGhgGJLf770MxIHaNzSYoAA2ARgydbn61t9nsRhaZMd1kh2Y7NLHYxEKzHZotUNuEL7Bg4sDKgYWDJhvYMFhZMLPYDtBkAwqBQgQUgJJBFAVKGhCAlAERDTQCBQMIQESDjAYAEFMgE0HLRVZB/7fPViEC5roBBn8ZDUQhUDI3uDJLGRDfbliCjAbRze/oeI7jeZ4R2W7zWwDkDBL9bwSpiFf6CFEIW4qWUvnH/DOVStXQ0PD3Y64d0PJ0UKfT3bAQlpSUbNmyZdeuXdd+snbt2oEDb7oo4uHiny9YqwEDAFA0LcM0AABCmKIRgBRomhZhhIBiGETJRFJgKKA9RIxULJICI0K0SCpVMxI50LSEljCIAQAZI0UIMYhp2U5dLpIjQFJKQjvDzfifKDtWiN2xrDIckklc42LXthgOScQSinLKprBUCb6hN/wbbG4GqwmsZrAYsdUMrA1sVmyzAGuHpibM2YBlAWOwNANCwNqx/c/LH8ci1gYAEowRQhgwsphunQJzHLJb2vTEhCQB8AC4u88bA2AAwC2XRrjWs9fSx4evO+z6n1//k7+/xbi1WVrTk3iLQ1p/nbUA/OXzNiARzPtYrWzttVEul1O3mynUqkKo0WgQQnq9Xq1WA4BOp/P3/+uQa19fX71e3/K6sbGx5Sc3/G2dOnUaN25c67tGp4/+twuPu7kHTU1N7nnizc3Ncrn8Fp0EHZXJZJJIJDTtRHdjjmGxWBiGYTpW12hrWK1WhNDd9RC6NLvdzvO8ROIsS3G1ati9RCKJjIw8fPhwy9sjR47Exsb+5Zi4uLjrDwgPD5e75fMtgiAIwrW0dv7ZnDlzlixZcunSpX379v3444+zZ88GgPr6+hEjRtTW1gLA448/vmHDhtzc3MuXLy9evHjOnDntmJogCIIg2khr+yJeeOGFurq6QYMGeXp6fv755wkJCQCAMbZarS3jTrt06bJy5cp58+bpdLqxY8e+9NJL7ZiaIAiCINpIa6dPtKE7nT6xfPnyWbNmeXl5tWsqJ7RixYpx48YFBAQIHcTRvvnmm9TU1IgIt1sbeu3atV27dm25y3QrmzdvDgwM7NOnj9BBHC0nJ4dhmKFDhwodxNF+++03g8EwatQooYP8wQWW5srMzKyoqBA6hQDWr19/9uxZoVMIYOPGjceOHRM6hQCysrKKi4uFTiGAnTt37t+/X+gUAsjNzd23b5/QKQSQn59//cQBwblAISQIgiCI9kMKIUEQBOHWSCEkCIIg3JoAg2WWLl26dOnSm023/7vKyko/P7+OtLVTK125csXb29t55pw6THV1tUqlksncblWd2tpauVyuULTx8lHOr76+XiwWu+HyETqdDiHk6ekpdBBH0+v1HMd5ezti59TJkycvXrz41scIUAh5ni8pKWl9YbNarW5YDMCNT9xms4lEIjdcWcZut9M0fdvloDoelmURQm64pA7HcQDgnieOMXbMWkKBgYG3vasWoBASBEEQhPNwu3tPgiAIgrgeKYQEQRCEWyOFkCAIgnBrpBASBEEQbs3ZNwDbsmVLXl5ecHDwnDlzbrbffcdTV1dXXFxcUVExcODA+Ph4oeM4TllZWVZWVkVFRXBw8PTp0zUajdCJHOTgwYO7d++ura318/ObPHlyaOiNt8DtwDZt2mSz2SZMmCB0EAc5fvz49avKTZo0yX2ubzabbdWqVSdOnNBoNBMnToyOjhY6kXO3CD/66KO5c+dGRUUVFBSkpqbyPC90IgcZNWrUa6+9lpGRsXfvXqGzONTo0aNPnDgRGhpaWFjYtWvXqqoqoRM5yKZNm4xGY2RkZFlZWUJCgrutMbtv374ZM2a88cYbQgdxnJycnBUrVpT9yW63C53IQSwWS0pKyurVq0NDQ00mU1FRkdCJAAAAOyu73R4UFLRr1y6MMcuy4eHh27ZtEzqUg7RMshk6dOhnn30mdBaHMpvN11737t37k08+ETCMUIYPH/7+++8LncJxmpubExMT33333bi4OKGzOM7y5ctnzJghdAoBLF26dODAgS2XOOfhvC3CkpKSurq6lJQUAKBpOjU11X2aR244pbqFVCq99tpqtSqVSgHDCKK6urqkpKRr165CB3GcjIyMadOmRUVFCR3E0UpLS5ctW5aZmanX64XO4jjbtm2bOnXqxo0bP/zww0OHDgkd5w/Oe8G9evWqt7f3taUH/P393aejjFi5cqXRaBw/frzQQRzn008/DQ0NDQkJSU9PHzFihNBxHOTAgQN5eXnPP/+80EEczcfHp1OnTgaDITMzMzY29sKFC0IncpDy8vIPPvhg586dOp1uxIgRq1atEjoRgDMPlmEYhmXZa2/tdrt7rjfmhn755ZeMjIzt27e71aqbjz/++Pjx448ePTpz5sxu3bo9+uijQidqd1ar9R//+Me3337rhisJp6enp6ent7yeMGHC0qVLv/zyS2EjOQZFUcnJyS0nGxMT09IfIHQoJ24RBgUF6XQ6k8nU8raysjIwMFDYSIQD5OTkTJ8+fevWrUlJSUJncSipVOrn55eWljZr1qz169cLHccRDh06VFpa+sQTT/Tq1Wv+/Pnl5eW9evVyq37CFgMGDCgrKxM6hYNotdprI+G7dOly6dIl7ATLfDpvizAqKio2Nnbjxo1TpkzR6/U5OTnz5s0TOhTRvvLz86dOnfrTTz/16dNH6CwOZTKZ5HI5AGCMi4qKunTpInQiR0hMTLz24H/Pnj2ff/75l19+6SbdAGazuWUlaI7jsrOz3eep8OjRowsKClpe79+/Py4uzhmW13fqRbezsrJmzpz5wAMPFBYW9ujR44cffhA6kYMsWrRo586d586d8/Ly8vX1fe+994YOHSp0KEfw9/fHGF+bRTd16tS5c+cKG8kx/Pz8kpOTvby8iouLaZrevXt36/cp6xg2bNiQkZFx6tQpoYM4SN++fdVqtZ+fX2FhoUQi2b17t5vMmtXr9YMHDw4ICAgKCsrOzl63bl1qaqrQoZy7EAJARUVFQUGBVqsdOHCgM9w4OEZ5eXlDQ8O1t1FRUW6yY9mRI0daNqZp4e/vr9VqBczjMDU1NYWFhU1NTWFhYcnJyW44bFin01VVVblJUxj+/MQNBkNYWFjfvn3d6hO3Wq179uyxWCwDBgzw8/MTOg6A8xdCgiAIgmhXbnQbQhAEQRB/RwohQRAE4dZIISQIgiDcGimEBEEQhFsjhZAgCIJwa6QQEgRBEG6NFEKCIAjCrZFCSBAEQbg1UggJgiAIt0YKIUEQBOHWSCEkCIIg3Nr/Awc05w/QP25CAAAAAElFTkSuQmCC", "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/html": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "using Plots, Distributions, LaTeXStrings\n", "\n", "μx = 2.\n", "σx = 1.\n", "μy = 2.\n", "σy = 0.5\n", "μz = μx+μy; σz = sqrt(σx^2 + σy^2)\n", "x = Normal(μx, σx)\n", "y = Normal(μy, σy)\n", "z = Normal(μz, σz)\n", "range_min = minimum([μx-2*σx, μy-2*σy, μz-2*σz])\n", "range_max = maximum([μx+2*σx, μy+2*σy, μz+2*σz])\n", "range_grid = range(range_min, stop=range_max, length=100)\n", "plot(range_grid, pdf.(x,range_grid), label=L\"p_x\", fill=(0, 0.1))\n", "plot!(range_grid, pdf.(y,range_grid), label=L\"p_y\", fill=(0, 0.1))\n", "plot!(range_grid, pdf.(z,range_grid), label=L\"p_z\", fill=(0, 0.1))" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### PDF for the Product of Two Variables\n", "\n", "- For two continuous **independent** variables\n", "$X$ and $Y$, with PDF's $p_x(x)$ and $p_y(y)$, the PDF of \n", "$Z = X Y $ is given by \n", "$$\n", "p_z(z) = \\int_{-\\infty}^{\\infty} p_x(x) \\,p_y(z/x)\\, \\frac{1}{|x|}\\,\\mathrm{d}x\n", "$$\n", "\n", "- For proof, see [https://en.wikipedia.org/wiki/Product_distribution](https://en.wikipedia.org/wiki/Product_distribution)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Generally, this integral does not lead to an analytical expression for $p_z(z)$. For example, [**the product of two independent variables that are both normally (Gaussian) distributed does not lead to a normal distribution**](https://nbviewer.jupyter.org/github/bertdv/BMLIP/blob/master/lessons/notebooks/The-Gaussian-Distribution.ipynb#product-of-gaussians).\n", " - Exception: the distribution of the product of two variables that both have [log-normal distributions](https://en.wikipedia.org/wiki/Log-normal_distribution) is again a lognormal distribution.\n", " - (If $X$ has a normal distribution, then $Y=\\exp(X)$ has a log-normal distribution.)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Variable Transformations\n", "\n", "- Suppose $x$ is a **discrete** variable with probability **mass** function $P_x(x)$, and $y = h(x)$ is a one-to-one function with $x = g(y) = h^{-1}(y)$. Then\n", "\n", "$$\n", "P_y(y) = P_x(g(y))\\,.\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- **Proof**: $P_y(\\hat{y}) = P(y=\\hat{y}) = P(h(x)=\\hat{y}) = P(x=g(\\hat{y})) = P_x(g(\\hat{y})). \\,\\square$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- If $x$ is defined on a **continuous** domain, and $p_x(x)$ is a probability **density** function, then probability mass is represented by the area under a (density) curve. Let $a=g(c)$ and $b=g(d)$. Then\n", "$$\\begin{align*}\n", "P(a ≤ x ≤ b) &= \\int_a^b p_x(x)\\mathrm{d}x \\\\\n", " &= \\int_{g(c)}^{g(d)} p_x(x)\\mathrm{d}x \\\\\n", " &= \\int_c^d p_x(g(y))\\mathrm{d}g(y) \\\\\n", " &= \\int_c^d \\underbrace{p_x(g(y)) g^\\prime(y)}_{p_y(y)}\\mathrm{d}y \\\\ \n", " &= P(c ≤ y ≤ d)\n", "\\end{align*}$$\n", "\n", "- Equating the two probability masses leads to identificaiton of the relation \n", "$$p_y(y) = p_x(g(y)) g^\\prime(y)\\,,$$ \n", "which is also known as the [Change-of-Variable theorem](https://en.wikipedia.org/wiki/Probability_density_function#Function_of_random_variables_and_change_of_variables_in_the_probability_density_function). \n", "\n", "- If the tranformation $y = h(x)$ is not invertible, then $x=g(y)$ does not exist. In that case, you can still work out the transformation by equating equivalent probability masses in the two domains. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Example: Transformation of a Gaussian Variable\n", "\n", "- Let $p_x(x) = \\mathcal{N}(x|\\mu,\\sigma^2)$ and $y = \\frac{x-\\mu}{\\sigma}$. \n", "\n", "- **Problem**: What is $p_y(y)$? \n", "\n", "- **Solution**: Note that $h(x)$ is invertible with $x = g(y) = \\sigma y + \\mu$. The change-of-variable formula leads to\n", "$$\\begin{align*}\n", "p_y(y) &= p_x(g(y)) \\cdot g^\\prime(y) \\\\\n", " &= p_x(\\sigma y + \\mu) \\cdot \\sigma \\\\\n", " &= \\frac{1}{\\sigma\\sqrt(2 \\pi)} \\exp\\left( - \\frac{(\\sigma y + \\mu - \\mu)^2}{2\\sigma^2}\\right) \\cdot \\sigma \\\\\n", " &= \\frac{1}{\\sqrt(2 \\pi)} \\exp\\left( - \\frac{y^2 }{2}\\right)\\\\\n", " &= \\mathcal{N}(y|0,1) \n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A Notational Convention\n", "\n", "- Finally, here is a notational convention that you should be precise about (but many authors are not.)\n", "\n", "- If you want to write that a variable $x$ is distributed as a Gaussian with mean $\\mu$ and covariance matrix $\\Sigma$, you can write this properly in either of two ways:\n", "$$\\begin{align*} \n", "p(x) &= \\mathcal{N}(x|\\mu,\\Sigma) \\\\\n", "x &\\sim \\mathcal{N}(\\mu,\\Sigma)\n", "\\end{align*}$$\n", "\n", "- In the second version, the symbol $\\sim$ can be interpreted as \"is distributed as\" (a Gaussian with parameters $\\mu$ and $\\Sigma$).\n", "\n", "- Don't write $p(x) = \\mathcal{N}(\\mu,\\Sigma)$ because $p(x)$ is a function of $x$ but $\\mathcal{N}(\\mu,\\Sigma)$ is not. \n", "\n", "- Also, $x \\sim \\mathcal{N}(x|\\mu,\\Sigma)$ is not proper because you already named the argument at the right-hand-site. On the other hand, $x \\sim \\mathcal{N}(\\cdot|\\mu,\\Sigma)$ is fine, as is the shorter $x \\sim \\mathcal{N}(\\mu,\\Sigma)$.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Summary\n", "\n", "- Probabilities should be interpretated as degrees of belief, i.e., a state-of-knowledge, rather than a property of nature." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- We can do everything with only the **sum rule** and the **product rule**. In practice, **Bayes rule** and **marginalization** are often very useful for inference, i.e., for computing\n", "\n", "$$p(\\,\\text{what-we-want-to-know}\\,|\\,\\text{what-we-already-know}\\,)\\,.$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Bayes rule $$ p(\\theta|D) = \\frac{p(D|\\theta)p(\\theta)} {p(D)} $$ is the fundamental rule for learning from data!" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- For a variable $X$ with distribution $p(X)$ with mean $\\mu_x$ and variance $\\Sigma_x$, the mean and variance of the **Linear Transformation** $Z = AX +b$ is given by \n", "$$\\begin{align}\n", "\\mu_z &= A\\mu_x + b \\tag{SRG-3a}\\\\\n", "\\Sigma_z &= A\\,\\Sigma_x\\,A^T \\tag{SRG-3b}\n", "\\end{align}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- That's really about all you need to know about probability theory, but you need to _really_ know it, so do the [Exercises](https://nbviewer.org/github/bertdv/BMLIP/blob/master/lessons/exercises/Exercises-Probability-Theory-Review.ipynb)." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "open(\"../../styles/aipstyle.html\") do f\n", " display(\"text/html\", read(f,String))\n", "end" ] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "anaconda-cloud": {}, "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Julia 1.9.2", "language": "julia", "name": "julia-1.9" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.9.2" }, "widgets": { "state": { "0c9c7079-a918-4ed4-ad45-8de3d7874019": { "views": [ { "cell_index": 12 } ] }, "261fa34c-df0a-4604-a2ff-35a9cb4d9fb2": { "views": [ { "cell_index": 20 } ] }, "42bb27af-6a2c-4a24-bcd3-1d8b8aec035c": { "views": [ { "cell_index": 20 } ] }, "50e83b63-be9e-4e83-987a-8b3ee9bbabd6": { "views": [ { "cell_index": 20 } ] }, "8053293d-ebc1-46b4-ac1c-afeb3a53d190": { "views": [ { "cell_index": 20 } ] }, "f496660b-92e4-41a4-a67b-b3ca19795f03": { "views": [ { "cell_index": 12 } ] } }, "version": "1.2.0" } }, "nbformat": 4, "nbformat_minor": 4 }