{ "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!\n", " - [Joram Soch et al ., The Book of Statistical Proofs (2023 - )](https://statproofbook.github.io/)\n", " - On-line resource for proofs in probability theory and statistical inference." ] }, { "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": [ "### Challenge: 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 that can be considered for its truth by a person. For instance, \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'}$$ as background information, 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, including:\n", " - (Representation). Degrees of rational belief (or, as we shall later call them, probabilities) about the truth value of propositions 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)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Why Probability Theory for Machine Learning?\n", "\n", "- Machine learning concerns updating our beliefs about appropriate settings for model parameters from new information (namely a data set), 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 (and information processing in general) 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", "\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", "\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": [ "### Events\n", "\n", "- Technically, a probability expresses a degree-of-belief in the truth value of an event. Let's first define an event. \n", " \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": [ "- Events can be logically combined to create new events. Given two events $A$ and $B$, we will shortly write the **conjunction** (logical-and) $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** (logical-or) $A \\lor B$ also as $A + B$, which is true if either $A$ or $B$ is true or both $A$ and $B$ are true. (Note that the plus-sign is not an arithmetic here but rather a logical operator to process truth values). " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Probability\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": {}, "source": [ "- $p(A|I)$ indicates the degree-of-belief in event $A$, given that $I$ is true. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- In principle, all probabilities are conditional probabilities of the type $p(A|I)$, since there is always some background knowledge. However, we often write $p(A)$ rather than $p(A|I)$ if the background knowledge $I$ is assumed to be obviously present. E.g., we usually write $p(A)$ rather than $p(\\,A\\,|\\,\\text{the-sun-comes-up-tomorrow}\\,)$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- The expression $p(A,B)$ is called the **joint probability** of events $A$ and $B$. Note that $p(A,B) = p(B,A)$, since $AB=BA$. Therefore the order of arguments in a joint probability distribution does not matter: $p(A,B,C,D) = p(C,A,D,B)$, etc." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note that, if $X$ is a variable, then an _assignment_ $X=x$ (where $x$ is a value, e.g., $X=5$) can be interpreted as an event. Hence, the expression $p(X=5)$ should be interpreted as the degree-of-belief that variable $X$ takes on the value $5$. " ] }, { "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 $p(x)$ is not necessarily $\\le 1$. E.g., a uniform distribution on the continuous domain $[0,.5]$ has value $p(x) = 2$." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- The notational conventions in PT are unfortunately a bit sloppy:( For instance, in the context of a variable $X$, we often write $p(x)$ rather than $p(X=x)$, assuming that the reader understands the context. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Probability Theory Calculus\n", " \n", " " ] }, { "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 of two events $A$ and $B$ with given background $I$ is \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\n", "\n", "$$ \\boxed{p(A,B|I) = p(A|B,I)\\,p(B|I)}$$\n", " \n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- PT extends (propositional) logic in reasoning about the truth value of events. Logic reasons about the truth value of events on the binary domain {0,1} (FALSE and TRUE), whereas PT extends the range to degrees-of-belief, represented by real numbers in [0,1]." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "- **All legitimate probabilistic relations can be derived from the sum and product rules!**\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Independent, Exclusive and Exhaustive Events\n", "\n", "- It will be helpful to introduce some terms concerning special relationships between events. \n", "- Two events $A$ and $B$ are said to be **independent** if the probability of one event is not altered by information about the truth of the other event, i.e., $$p(A|B) = p(A)\\,.$$\n", " - $\\Rightarrow$ If $A$ and $B$ are independent, then the product rule simplifies to $$p(A,B) = p(A) p(B)\\,.$$\n", " - $A$ and $B$ with given background $I$ are said to be **conditionally independent** if $p(A|B,I) = p(A|I)$. In that case, the product rule simplifies to $p(A,B|I) = p(A|I) p(B|I)$.\n", "- Two events $A_1$ and $A_2$ are said to be **mutually exclusive** ('disjoint') if they cannot be true simultanously, i.e., if $p(A_1,A_2)=0$.\n", " - $\\Rightarrow$ For mutually exclusive events, probability adds (this follows from the sum rule): $$p(A_1+A_2) = p(A_1) + p(A_2)$$\n", "- 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 $$p(A_1+A_2+\\cdots +A_N)=1$$\n", "- If a set of events $A_1, A_2, \\ldots, A_n$ are both **mutually exclusive** and **collectively exhausitive** events, then we say that they **partition the universe**. Technically, this means that $$\\sum_{n=1}^N p(A_n) = p(A_1 + \\ldots + A_N) = 1$$\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- We mentioned before 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_. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Marginalization\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Let $A$ and $B_1,B_2,\\ldots,B_n$ be events, where $B_1,B_2,\\ldots,B_n$ partitions the universe. Then\n", "\n", "$$\\sum_{i=1}^n p(A,B_i) = p(A) \\,.\n", "$$\n", "- This rule is called the [law of total probability](https://en.wikipedia.org/wiki/Law_of_total_probability). Proof:\n", " \n", " $$\\begin{align*}\n", " \\sum_i p(A,B_i) &= p(\\sum_i AB_i) &&\\quad \\text{(since all $AB_i$ are disjoint)}\\\\\n", " &= p(A,\\sum_i B_i) \\\\\n", " &= p(A,\\text{TRUE}) &&\\quad \\text{(since $B_i$ are exhaustive)} \\\\\n", " &= p(A)\n", " \\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- A very practical application of this law is to get rid of a variable that we are not interested in. For instance, if $X$ and $Y \\in \\{y_1,y_2,\\ldots,y_n\\}$ are discrete variables, then\n", "$$p(X) = \\sum_{i=1}^n p(X,Y=y_i)\\,.$$" ] }, { "cell_type": "markdown", "metadata": {}, "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** of $X$. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note that marginalization can be understood as applying a \"generalized\" sum rule. Bishop (p.14) and some other authors also refer to this as the sum rule, but we do not follow that terminology." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Of course, in the continuous domain, marginalization becomes\n", "$$p(X)=\\int_Y p(X,Y) \\,\\mathrm{d}Y$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### 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**, named after its inventor [Thomas Bayes](https://en.wikipedia.org/wiki/Thomas_Bayes) (1701-1761). While Bayes rule is always true, a particularly useful application occurs when $D$ refers to an observed data set and $\\theta$ is set of unobserved 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) = \\theta^y (1-\\theta)^{1-y}\\\\\n", "\\end{align*}$$\n", "\n", "- Next, we use Julia to plot both the sampling distribution $$p(y|\\theta=0.5) = \\begin{cases} 0.5 & \\text{if }y=0 \\\\ 0.5 & \\text{if } y=1 \\end{cases}$$ and the likelihood function $$L(\\theta) = p(y=1|\\theta) = \\theta \\,.$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "using Pkg; Pkg.activate(\"../.\"); Pkg.instantiate();\n", "using IJulia; try IJulia.clear_output(); catch _ end" ] }, { "cell_type": "code", "execution_count": 2, "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+ceBDEu5DCAmnAUUFRRSUAPogl6CgVgSCIHgUL6i29tC2nlX78mmVPh48HrVYj0fr+bS1PhUVDxBpiyiIVg5RURTlvsMRINns74/5Pvvkx61lswE+7782yyT5ZJid2Z2dnckkjnuxNaCqA12DK0IABo5ly5YdOnRIJBKZmJj8+9//Tk1NxfvxiNDIyMh//vOfFy9eDAwMjIyMRAix2ewzZ84ghCZMmODk5MTlcpOTkxFCH3zwAYO/oq/I5Gh+CrFlInscV62hgelogAqDhhCAgYPL5WZlZV24cEEikXzxxRcWFhZ4/7Fjx/BFXkhIiJOTE5We6ho9ffp0ZmZmY2PjjBkzfH19VXnkS+99eZ8w1EAr7eFiEPRgIBR3AADFwMBg0aJF7XZ6eXnhDRsbGxsbm47vcnBwcHBwoD04JfqjnDzyWH4fFt0FvQDnSgCAgUYsRYtSiPgpbNNOnhwBoD1oCAEAA82qNMJ/qNpsK6jfQK9A1ygAYEA591x+t5K8FwSVG+gtKCsAgIHjVRP54W3isj9HB+o20GvQdQAAGCDkJIq6RXziwHbmwRAZ8AagIQQADBA7s+VtcvTZeKjWwJuB7gMAwEBwv5rc+ZC4OwcW3QVvDE6dAAD9XguB3r1FxE1mD9OHZhC8MWgIAQD93qd3CLshapEjoEIDbwO6RgEA/duVV2RCEfkgBGoz8JZoLDr3798/dOhQW1vb/Pnzp06d2u6vVVVVP//8c2ZmJpvNnjZtWkhISG9WjQIAAEWVLWjZb8RJH7ahSi8ZC1QaXT0JT58+9fb2Hjly5N/+9jeRSJSSktIuwYkTJ5KTk11cXBwdHT/66KOvvvqKpkgAAAPYst+IRaPUfMzhNBq8PbquCA8cODB37tw1a9YghKqqqnbt2uXt7a2Y4OOPP6aWEjY0NNyyZcumTZtoCgYAMCDF58tfN5E/+kKnKPhL6LoiTE9Pp1o+b2/v27dvt/9i1v++uqKigs/n0xQJAGBAKhCTm+8Rp7xh0V3wV9F1JlVWVmZsbIy3eTxeTU1Na2urpqZmx5RFRUXbtm07depUVx/19OnT/Pz8zMxM/JLFYu3atYv6cJXS3Kwml7MbGxuZDqRXmpqa+st92ZYWJJOpNTZKmQ6kV+RyuYYG3LCil0yOFqQQX05kjzbsH2UYqDK6GkItLa22tja83drayuFw1NXVOyYrLy+fPn362rVr/f39u/ooMzMzDQ2N4OBg6pPNzc1Vc+FQrTaSxSJ0dHSYDqRXCKLfhKqhIeNwSB2dToqQCpJK+0eD3a9tvkfwtVDMGLgYBH2ArubE0tKyqKgIbxcVFVlYWCj2hWJVVVXTpk2LjIz87LPPuvkofX19S0vLsLAwmkLtQywWiRDR8ZeqJhaL1Y9C7UcZ21+us/uv38vIY0/k90Ng0V3QN+iqWUQi0ZkzZwiCQAidPHlSJBLh/RcvXnz9+jVCqK6uLiAgwN/ff/PmzTTFAAAYeOrb0MJbxCEPthksugv6CF0NYVRUFIvFEgqFnp6eGRkZn376Kd6/cuVKPHBm9+7dWVlZ58+fHzFixIgRI8aOHUtTJACAgeT9NGKmABbdBX2Jrq5RHR2d1NTUjIwMqVTq6upKjR1IS0vD41w+/vjjqKgoKj30JgEAevRTofxeFSy6C/oYjeWJxWK5ubm12ykQCPAGl8vlcrn0fTsAYIB53US+n0ZcnA6L7oI+Bt0LAIB+QE6ihSnEage2Kx96j0Afg4YQANAP/OOhnCDRp7DoLqABdDEAAFTd/Wpydw4sugvoAqdXAACV1ixDkTeJvZPZ1nrQDAJaQEMIAFBpa+4QQp5auA1UVoAu0DUKAFBdia/Jy6/IB8FQUwEaQfECAKioyha0NJU45cM26mS6fgD6DPQ2AABUEYnQ0lQiylbNGxbdBTSDK0IAgCo6mCcvbiZ/8oM6CtAOChkAQOXk15Fbs4jfZnNg0V2gBFDKAACqpU2O5qcQ24RsuyHQKQqUARpCAIBq2XyPMNVG0aOhdgJKAl2jAAAVcquU/OEp+SCEAxeDQGngnAsAoCrq2tC7t4jDHmy+FtOhgMEEGkIAgKpY+Qcx20pthgCuBoFSQdcoAEAlnCqQ/1lNZsKiu0DpoMwBAJj3spFcfYe4GsDRhjoJKB10jQIAGCYnUdQt4tNxbCdj6BQFDICGEADAsB1/ygkSrR4H1RFgBnRDAACYlFVFxuUQmUGw6C5gDJyCAQAY0yxD81OIvZPZVrDoLmAONIQAAMasuUO48NTmjYCKCDAJukYBAMxIKCKvvSbvh0AtBBgGRRAAwIDKFhTzO3Hah22gznQoYNCDHgkAgLKRCEXdki21U/OCRXeBCoCGEACgbAfy5FUtaPMENtOBAIAQdI0CAJTsUR35ZRbxx2wOB87DgWqAkggAUB6pHL17i9juwraFRXeByoCGEACgPBszCQsdtWV2UPMAFQJdowAAJfmtjDxZIL8fDONEgWqB8zIAgDLUtqIFKcRRT46pNtOhAPD/g4YQAKAMK/4g5lirBQyFW4NA5UDXKACAdj8UyLNrYNFdoKKgXAIA6FXYQK5JJ67PhEV3gYqCrlEAAI0IEi26Rax1ZDtyoVMUqChoCAEANNr+QM5WQ584QFUDVBd0VQAA6HKvityXR2TM4bDgahCoMDhNAwDQokmGIm4SB/4Gi+4CVQcNIQCAFp+kE+6maqHDoZIBqg66RgEAfe8/L+VJxbDoLugfoJgCAPpYhQSt+IM4O5UDi+6CfgF6LQAAfYlEaEmqbLkdy9MMbg2C/gEaQgBAX9qfKy+XoE2w6C7oP6BrFADQZ3JryW33ibRAjjqcY4P+A0orAKBvtBJoQQqx3YU90gA6RUF/Ag0hAAMHQRCbNm2yt7efNGnSr7/+2mmaX375ZdKkSWPHjt28eTNBEHhnQ0NDTEyMnZ3d1KlT79y583bfvjGTGKanthQW3QX9DXSNAjBw7Nmz5/LlyxcuXCgoKIiIiLh7966dnZ1igry8vKioqLNnz9rY2ISHhxsZGX3yyScIodWrV5eXlycmJiYlJc2aNauwsFBfX/+Nvjq5hDz3nLwfDFUK6H/g3A2AgePgwYNbtmyxtbWdOXNmcHDwkSNH2iU4fPiwSCSaMWOGnZ3d5s2bv/vuO4SQWCw+derUN998M3z48GXLltna2p47d+6NvrdeqrYklYifwuZp9dlvAUBpoCEEYICQSCTPnz93dnbGL4VCYXZ2drs0OTk5igmePn3a2tr6/PlzFos1evRovN/Z2TknJ+eNvvrDu+zgYWozBXBrEPRL0I8BwABRWVmJEBoyZAh+aWhoiPe0S2NoaEglIEmysrKysrKSehdCyMjIqKCgoKtvKSoqcnd3Z7P/7+mIiRMnrv/+wqM69YOuzQ0Nffhr+lhjYyPTIfQAIvzrOo1QS0tLXb2HmR2gIQRggMAtXFNTE76919DQYGRk1DENVVk0NDTgPYaGhk1NTVQasVjM5XK7+pahQ4fGx8fb29vjl5qamppaWr/xJFyDN7unqHxvetdT+SDCv+7tIoSuUQAGCAMDA2Nj48ePH+OXjx8/Hj58eLs0NjY2iglMTU319PSsra0bGxtLS0u7eSOFxWIZGBgY/ZeOjo4aQprw9Dzoz6AhBGDgWLRoUVxcHEEQxcXFZ8+eXbhwIUKorq5uzZo1YrEYIbRgwYKzZ8+WlJQQBBEXF4cTmJiY+Pv77969GyGUlZX1+++/R0REMPtDAFCmftAQVlVV1dTUMB1Fr0il0ubmZqaj6K3U1FTFDjFVVltb21/KAELo8uXLTH315s2bJRIJn88fO3bsypUrPTw8EEJNTU0nTpyQSCQIIR8fn2XLltnb2/P5fPzQIX7jvn37bt68yefzp0+fvn//fgsLi95/aWVl5d27d+n4OX1FIpGkpKQwHUV3SJK8evUq01H0IDExkSRJpqPozs2bN1taWt7mnSQ9ZDLZwYMHIyIiPv/887Kysk7TpKSkLF68eOnSpWlpad181KRJk/z9/ekJs49dulegEfuY6Sh6y9XVtfucVx3z42/a/T2J6Sh6pa2tTV1dndkYGhsb29rauknQ1tbW1NTUcb9YLJbJZN1/+Lhx4x4+fKi458cffwwKCnqLOJUmMzPTycmJ6Si6U1NTg8cuqbIhQ4bU1tYyHUV3xo8fn5WV9RZvpOuKcOvWrd99951IJKqtrfXx8ZHJZO0SpKenBwYGTp48eeLEiQEBAX/++SdNkQAw2Ojq6nY/TE5dXV1HR6fjfn19fWo4aO+Rqn2VAECPaBk1KpFIDhw4cPXqVRcXl5CQEFtb20uXLs2ZM0cxzZ49ez788MPly5cjhJ49e7Z3796OD/8CAAAAdKPlivDJkydtbW1CoRAhpKam5unpmZ6e3i7NnTt3vLy88Lanp+dbT28IAAAA/BW0XBGWlZVxuVw1tf+bZoLH41Ejsynl5eXGxsZ4m8/nd0xAefXq1b1796jEampq9vb2PT4gyYhalgExeoGv7wqmA+mV/Pz8Dz/80MDAgOlAevZYZ7RYb6iv71dMB9IzkiTlcjnTUdCovLw8JiZGV1dXcU95efm0adMYjKp7DQ0Nz549U+UIZTJZU1OTKkeIEGpubg4ODuZwVPfp88LCwhUrVrR7lDA4OHjlypXdv5GWn6Strd3a2kq9bGlpUTxsMC0tLSpNS0tLp3cssN27d2dmZvJ4PPySxWKNHz+exVLF8a4kSb58+XLYnPVMB9Irr169Mjc3V+ViTWlsbJRIJHx+/8jYZcuWMR0CjeLj43V0dBQPwJaWltraWnNzcwaj6h5BEK9fv7a2tmY6kO4UFhZ28/imKggLC1PxCF++fCkQCNq1Dr2JmZZK0NLSsqamprGxUU9PDyFUVFTk4uLSLs3QoUOLioomTZqEEwwdOrSrTwsLCwsLC6MjTgDAmwoKCmI6BAD6GC3XVSNGjHBwcDh9+jRCqKSkJCkpSSQSIYRKS0t/+uknnEYkEv3www+4H+nkyZM4AQAAAKBkajQNfU5JSZk7d66zs3N2dvaCBQtiY2MRQleuXImMjKytrUUI1dTU+Pj4aGpqymQydXX1GzduqP4sdgAAAAYeuhpChFBdXd2DBw+srKxsbGzwnpaWlurqaktLS/ySIIjMzEwWi+Xs7Kya9/wAAAAMeDQ2hAAAAIDqU+kRg62trVevXq2vr/fz8+s4Jo0kyaSkJOqlQCCws7NTZngFBQW///67QCCYOnUq9ayIovr6ejx/YEBAALNPKUil0qtXr9bU1Pj6+lJX5IqSkpKoUyJLS8sxY8YoN8D/kclkOTk5eEKirtJkZGTk5OSMGzcOP6vKoIqKiry8vOHDh3c6IrGwsPDZs2fUSw8PD01NTSVG12devHhx69YtMzMzPz+/TqeeaWxsTExMlEqlAQEBims/5eTkZGRkjBo1asqUKbRGKBaLExMTUdfHWn5+/v379zU1NT09Pakh6Lm5udSDW2w2u5si11cRkiQZEBCguPoj1tjYqPiw9ejRo6nxg2VlZfjOUUBAAK3lhyCIGzdulJWVeXl5DRs2rN1f6+vrMzIyFPeMGzfO1NS0pKQkLy+P2unq6kpfXSeRSB4+fMhisTqOvsRIkkxOTn716tWUKVNGjhxJ7e85D/tmijcaSCQSV1dXDw+PqKgoY2Pje/futUsglUoRQt7e3n5+fn5+fvv371dmeBcuXDA2Nl6+fLmjo2NYWFjHBCUlJQKBIDg4ODg42NraurS0VJnhKWpra3N3d3d3d1+8eDGXy01PT++Yhs1me3p64pzcs2eP8oPEUlNTtbS0eDyerq5uV2m2bdsmEAhiYmIEAsH27duVGV47s2bN0tbW1tXVjY2N7TTBl19+KRAI/P6rurpayRH2ievXr3O53KVLlwqFwhkzZsjl8nYJqqurR40aNXPmzPDwcHNz88LCQrz/+++/NzU1fe+990aNGrVq1Sr6IiwpKbGyssLHmpWVVcdjDZeZefPmBQYGGhoapqam4v0LFy60s7PD/53Zs2fTF2Fpaam1tXVQUFBISIiVlVVJSUm7BNnZ2ZqamlRR+c9//oP3P3jwwNjYeNGiRd7e3kKhsLm5maYI5XL5jBkznJ2dly5damxsfP369XYJcnNzqfDwaQ2eqfjYsWN8Pp/6U35+Pk0RxsfHa2hocLncyZMnd5UmLCxs/Pjxy5cv5/F4Fy5cwDt7k4eq2xCeOHHCyclJKpWSJLl169Y5c+a0S4Abwvr6eiaiI+3t7U+dOkWSZH19PZ/Pv3PnTrsE69atCw8Px9vh4eEbNmxQdoj/de7cubFjx+JZmHfs2BEQENAxDZvNrqioUHpo7dXV1VVUVGRmZnbVEFZXV+vo6OTl5ZEkmZubq6urW1dXp9wY/+fFixdSqXT27NndNIQffPCBkqPqc5MmTTp48CBJks3NzcOGDbt27Vq7BF9//fWMGTPw9rJly95//32SJFtbW01NTZOTk0mSLC0t1dHLkrTOAAANQklEQVTRoRrIPrd+/XrqWJs3b9769evbJSgsLMQ1CUmSGzZs8PX1xdsLFy6Mi4ujKSpFGzZsmDt3Lt6OiIhYu3ZtuwTZ2dkWFhYd3ygSiTZu3EiSpEwmEwqFR48epSnC69evCwQCPBX7oUOH3Nzcukl8/PhxW1tbfEp07NixwMBAmqJSVF5eLhaLjx492lVDeOfOHR6Ph1uEM2fO2Nvb4/2hoaFUHrq4uHSah6o7RCUhISEoKAg/7h0aGnr58uVOJ+xIS0u7efMmHomqNM+ePXvy5ElwcDBCyMDAYPr06QkJCe3SXLx4MTQ0FG+LRKKOCZQmISEhMDAQz8UTGhp67dq1tra2jsnS09OTk5OZXe1oyJAhfD6/mwRJSUk2Nja459be3t7Kyio5OVlZ0bVnbW3d43QElZWVV65cyc7OJvvnzfjKysr09HT8dJO2tvbMmTM7luSEhASqqIeGhuIEGRkZMpnM29sbIWRmZjZ58uRLly7RFOTFixep5686PdaGDRtG/afMzc0Vp/soKiq6cuVKQUEBTbFhilnUVW0gk8lu3LiRlpZGLeVGkuSlS5fwG9lsdnBwMH3VSEJCwsyZM/HEJqGhoXfu3KmoqOgq8ZEjR5YsWULdD8K9vllZWQRB0BQeQsjExKT7JwsSEhKmT5+OO2bnzJnz5MmTZ8+ekSRJZX43eai6DWFxcTF1N8vS0lIqlXb8x5iZme3fv3/Tpk3Dhw8/f/68MmPjcrna2tpUeMXFxR3TKMbfMYHStItELpd3nNDOxMQkPj5+y5Ytw4cPP3v2rNJj7K3i4mLFuReYzdgesdnswsLCgwcP+vv7T506tb+s/qiopKREQ0ODOjvpTVEvKSmRy+V4J1VX4v00Bdn7Y00sFsfFxVFT/6irq2dkZBw4cMDFxWXhwoX0zY3Xmwj19PT27t27cuXKUaNG4bmXa2pqWlpalFONKEZoZGSkra3d1Xc9f/789u3beElnrK6u7ttvvw0NDRUKhWVlZTRF2CPFykFbW5vL5RYXF/cyD5kcLHP79u0VKzqZlvPixYsCgYAgCOqZCnx/vt1aThwO5/Xr1/hPp06dWrJkyaxZs5QzGIEgCMXRMWw2u+M6U4ppOk2gNIo5iTc6BvPq1Suckz/99FNUVNQ777yDZwVSNe1ynsPhMJixPVq3bt3GjRsRQhKJxNPTc+fOnVu2bGE6qDejWHhQ10Vd8VDFPWa9OUb6MMjeHGutra1hYWHu7u6LFi3Cew4dOoSLfWVl5YQJE86dOxcREcFIhGPGjKHGVX3xxRcxMTEPHjzAF1jKqUba/aO7ObIOHz48Y8YMaunmhQsXRkVFIYRkMllwcPCmTZsOHz5MU5Dd67Ry6GUeMnlFOG7cuJOdMTU1RQiZm5tTl4Dl5eUsFsvMzKzdJ1AD2MLDw/G8usqJ3NzcvLa2lsrQ8vLyjoNazc3NKysrqQRvtOR331LMSbzRMRgqJ0UikUwme/r0qTIj7D3F34KYztgeUbmqra0dFBR0//59ZuN5C2ZmZi0tLQ0NDfhlV0Vd8VA1NTVls9kd/1P0TUbam2NNKpWGhYXp6+sfOXJEsVrEG3w+39fXl75/UI8RKo7FjYiIyM7OlslkPB6Pw+EopxpR/H81Nzc3NDR0+l0ymez48eNLlizpGDmHw5k7d+6DBw9oirBHij9BJpNVV1dbWFj0Mg+ZbAj19PQcOqOhoYEQ8vb2vnbtGk557dq1KVOm4F7++vr6lpaWdh+Fc7/TBwPoMHLkSBMTk5s3byKEZDJZcnIyHngtlUqpe2w+Pj742QkcP75Zwoh2Oenm5oY7dcVicceczMnJaWtr62bqV0ZQoXp4eOTk5JSXlyOESktL8/Ly6B6X/6YUy4CirKwsKysr5cfzF5mbm9vZ2eHyI5fLb9y4gYs6rmhwGm9v745F3dnZubGxMTs7GyHU3Nz8xx9/0Pdwgo+Pj2IJp461mpoaPKSOIIhFixaxWKzTp093eluXIIg///xTIBDQFGGnWaQYoaKsrCw8Gz6LxfL09Oz0p9ER4Y0bN3Dn8NWrV21tbXGD0dDQIJFIqGSJiYkEQcycObPTD7l37x59edgVKg+9vb2Tk5Px9UlKSgqfzx85ciSLxfLy8uo5D2kZ39MXamtrraysli5dumPHDiMjI/wIDkmSQqFw7969JEmeOXMmMjJyx44d69at4/P5Sh6WuX//foFAsGvXrsDAQDc3N4IgSJK8du0aNdzx0aNHQ4YMWb9+/fr16w0NDekbVdyjhoYGGxubd999NzY21tjY+Ndff8X73d3dv/nmG5Ikf/7553nz5m3fvn39+vWmpqZr1qxhKtSampro6OiQkBAOhxMdHU39TydPnrxz5068vWTJEhcXl7i4OKFQGBMTw1SoJEmeOHEiOjra2trazc0tOjr69u3bJEleu3ZNT08PJwgMDFy3bl1sbGxgYCCPx6Nv2CStTpw4YWZmtnPnzrCwMAcHh9bWVpIk8U0sXOxfvHjB5XLXrFmzZcsWAwODrKws/MZNmzaNGTNmz549Pj4+tD6ckJ+fj4+1DRs2KB5rOjo6+DGAzZs3czicqKio6Ojo6Ojozz//nCRJuVzu5ub2xRdfbN++3cPDw87OTiwW0xehoaHhunXrcIR42DNJknp6engU7vbt25cvX/7NN9+sWrVKX1//X//6F05w/fp1Q0PDr7/+Ojo6eujQofQ9gdPW1ubg4BAaGrpz505zc/Pjx4/j/bNnz1YchRsUFIRzj7J48eLVq1f/4x//iIyM1NfXz8jIoCnC7Ozs6OhoLy8vU1PT6Ojob7/9Fu+n8pAgCDc3t9mzZ+/evdvKyop6oO7GjRs95iF769atymq534yWllZkZGRxcXFDQ8Pf//53T09PvN/ExEQoFJqYmPD5/NbW1vLycj09vc8++2zx4sXKDM/V1dXe3v7Ro0dOTk579uzB9ya1tLRGjx7t5OSEEOLxeCKRKD8/X11dfe/evba2tsoMT5GGhsb8+fNLSkrEYvHWrVt9fX3xfj6f7+zsbGpqyuPxpFJpeXm5rq7u6tWro6OjmQpVLpdXVVWNHDly1qxZFhYWVlZWDg4OOFShUIj7zN955x1dXd3nz58HBQWtXbu206kMlKO+vl5DQ8PDw2PChAkWFhb29vbGxsaKZQD31YjFYqFQGB8fr8oLFXXD0dFRKBTm5uba2dnt27cPL6mmoaExatQoPKGBoaFheHj406dPSZLcvXv3+PHj8RunTp1qYWHx5MkTLy+vr776qtMn8fsEj8cLDQ3Nz8/ncDiKx5qlpeWkSZP09fXV1dVdXV0FAoGFhYWFhcXQoUPHjx+vpqbG5/MrKipaW1vxg8gdV4vr2wgfP37MZrP37t1LTf1haWnp5uZmYGBgbm7e0NCAO5C//vrr6dOn4wQ2NjZ+fn65ublmZmYHDhzofkz1X8Fms+fPn19dXV1dXf3555/PmTMH7zc2Np44cSLVnSiRSCIiIhQfmbe0tKysrKyurra3tz948CB9s5q0trY2NTWNGzfOx8fHwsLCxsZmxIgRSCEP1dTUIiIicDauWrWKut1rY2Mzbdq0nJycbvIQplgDAAAwqKnu4xMAAACAEkBDCAAAYFCDhhAAAMCgBg0hAACAQQ0aQgAAAIMaNIQAAKASOk5wAZQDGkIAAGDYixcvPv744x9//BFPTguUTKVXqAcAgAGvoqIiJCTk6tWrfD7f39+/oKBAcXV1oARwRQgAAExasWJFTEwMnvFEU1Pz7t27TEc06MAVIQAAMObevXtJSUknT57EL0tKSvrjupX9HVwRDi7fffcdj8d7/fq14s7Fixc7OjrStyopAKArcXFx77zzDl4QRiqV5uXl0TehKOgKNISDy7x58yQSyZEjR6g91dXVZ86cEYlEistyAgCUgCCIK1euiMXi2NjY2NjYjRs3SiQSR0dHpuMadKDuG1wMDQ3nzp17+PBhvHAzQujYsWMEQSxdupTZwAAYhLKzs6urq7dt24bXh9LU1LSxsRk+fDjTcQ060BAOOitWrHj9+nViYiJ+efTo0VmzZiltTWMAAKWoqMjU1NTR0dHIyMjIyCg1NZVaPAgoEzSEg46bm5uzs/OhQ4cQQjdv3nz06FFMTAzTQQEwGNXV1VFPSrx48SIzM3PlypXMhjQ4QUM4GL333nuXLl0qKir6/vvvraysqFVAAQDKZG5uPmTIELy9Z8+etWvXUkvgAmWChnAwioyMNDAw2LVr1y+//LJ8+XL6lg4HAHTD3d1dIpEghNLS0oqKijZs2MB0RIMUrFA/SH300Uf79u1js9kvX76Ek1AAmJKcnPzo0SNNTc2oqCgOBx7sZgY0hINUbm6ug4NDcHDw+fPnmY4FAACYBF2jg1RpaSlCKDo6mulAAACAYXBFOBgRBOHn51dVVfXw4UM1NTWmwwEAACZBl/SgM3/+/OTk5KqqqsTERGgFAQAArggHnfPnzzc0NLi7u8NSLwAAgKAhBAAAMMjBYBkAAACDGjSEAAAABjVoCAEAAAxq/w8wq4og7bTROAAAAABJRU5ErkJggg==", "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": [ "### Revisiting the Challenge: 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 information is $p(D=1)=0.01$, $p(T=1|D=1)=0.95$ and $p(T=0|D=0)=0.85$. We are asked to derive $p( D=1 | T=1)$. We just follow the sum and product rules to derive the requested probability:\n", "\n", "$$\\begin{align*}\n", "p( D=1 &| T=1) \\\\\n", "&\\stackrel{p}{=} \\frac{p(T=1,D=1)}{p(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", "- Distributions can often usefully be summarized by a set of values known as moments of the distribution. \n", "\n", "- Consider a distribution $p(x)$. The first moment, also known as **expected value** or **mean** of $p(x)$ 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 second central moment, also known as **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 a mixed central moment, 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", "\n", "\n", "\n", "- Clearly, if $x$ and $y$ are independent, then $\\Sigma_{xy} = 0$, since in that case $\\mathbb{E}[x y^T] = \\mathbb{E}[x] \\mathbb{E}[y^T] = \\mu_x \\mu_y^T$.\n", "\n", "- Exercise: Proof that $\\Sigma_{xy} = \\Sigma_{yx}^{T}$ (making use of $(AB)^T = B^TA^T$)." ] }, { "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 &\\qquad \\text{(SRG-3a)}\\\\\n", "\\Sigma_z &= A\\,\\Sigma_x\\,A^T &\\qquad \\text{(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, assume two jointly continuous variables $X$ and $Y$, with joint PDF $p_{xy}(x,y)$. Let $Z=X+Y$, then\n", "$$\\begin{align*}\n", "\\text{Prob}(Z\\leq z) &= \\text{Prob}(X+Y\\leq z)\\\\\n", "&= \\int_{-\\infty}^\\infty \\biggl( \\int_{-\\infty}^{z-x} p_{xy}(x,y) \\mathrm{d}y \\biggr) \\mathrm{d}x \\\\\n", "&= \\int_{-\\infty}^\\infty \\biggl( \\int_{-\\infty}^{z} p_{xy}(x,t-x) \\mathrm{d}t \\biggr) \\mathrm{d}x \\\\\n", "&= \\int_{-\\infty}^z \\biggl( \\underbrace{\\int_{-\\infty}^{\\infty} p_{xy}(x,t-x) \\mathrm{d}x}_{p_z(t)} \\biggr) \\mathrm{d}t\n", "\\end{align*}$$\n", "\n", "- Hence, the PDF for the sum $Z$ is given by $p_z(z) = \\int_{-\\infty}^{\\infty} p_{xy}(x,z-x) \\mathrm{d}x$.\n", "\n", "- In particular, if $X$ and $Y$ are **independent** variables, then\n", "$$\n", "p_z (z) = \\int_{-\\infty}^{\\infty} p_x(x) p_y(z - x)\\,\\mathrm{d}{x} = p_x(z) * p_y(z)\\,,\n", "$$ \n", "which is the **convolution** of the two marginal PDFs. " ] }, { "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": [ "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Code Example: Sum of Two Gaussian Distributed Variables\n", "\n", "- Consider two independent Gaussian-distributed variables $X$ and $Y$ (see [wiki:normal-distribution](https://en.wikipedia.org/wiki/Normal_distribution) for definition of a Gaussian (=Normal) distribution):\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", "$$\n", "\n", "- We illustrate the distributions for $X$, $Y$ and $Z$ using Julia:" ] }, { "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/jlp2lagW8NSMFPnKqv5jNOUQ4u6u795a7FrOj/kbsQc1PyE8cYK/CJS5JktiduyJ/UBgWho6Ojk5PT1atXtVotQqi4uPjuUaOfffbZjBkzxowZgxDauHGjj4/PtWvXfH197z6bi4tLZGTkjBkzDLk0QohlWWU2jBT1xCVdnT7trC7znD49WWyoUYf2FIpzhaxBLuP+gWmlNBQaX3EFvi2KoqjMFmFjHijn1/x2kiRZzhM39Cfv4Ycf3rZt24ABAyRJ+uWXXyZNmoQQEgTh1KlTffr00Wg0rq6uhYWFjQcXFhbSNO3s7GyqqoFtIfqGkvWvUc4e6vBebuNf5Rw97Oztib6+6pdNJR/N9Zi6hHbzlrtGAKzPsWPH1q9fL3cVLSCEEELa7bFs5o033hgwYIAp6jE0CBcvXjxq1Kji4uKrV6/W1tZOnjwZIVRRUTF06NDMzMyIiIh58+Y9+OCDJSUl3t7eW7duXbRokb29vSkqBraGkPL/fMAGRrk8+nzjA1x9PUKI0ji4Pf1aXdKuGx+86vbMPI3CukkB6Lrk5OTq6urnn39e7kKM4PPPP79w4YLMQdirV6/Lly/v3r3bwcHhoYcesrOzQwi5uroePHgwICAAIdS/f/+MjIzExMTa2trp06fHxMSYolxge2r2/ySUXvWY9lYLn8PYYfBo1je44j8fOMaPdbr/abNXB4B1i4iIeOqpp+Suwgh27dplupN3oFPex8ensSHYhGXZhISEpg89PDwef/xxo5UGFECfea724DaPF1dgptW7BaqQ7h4vrir7Yok6uq8qEKblAACMDNYaBbIRK0vKv1nj8tiLtLNH20fSTi6O9zxV+dOnSJHbMwEATAqCEMiDCHzZ5hX2Qx9WR8Qacrx9nxFI4BvOHTZ1YQAApYEgBPKo/OVzysXDcegYQ78AY6cHn6n6/UvCc6asCwCgOBCEQAZC6dWGc4dd//YP1JHlY1Qh3Rn/sNqD20xXGABAgRQ3gxVYgtqDPzsMuBerNB39QqcHni774k37QQ/Qzu6mKAwA23a2lNy3UzDbnXYVhQ48zPRwbf737rVr1xwcHJycnDIzM0NCQmSfWQ9BCMxNqq+pP5voOaszm5Mwbt72/e6p3v6124Q5Ri8MAJvX1xMnjWVFcyUhQ6GIuxZWqampOXr06JYtW6ZNm9a/f/9p06Z988038q4tDEEIzK328B+a6P60k2vnvtwxYVzJhnlcYQZMpQCgE1xVyGxNQralm2+ZmZkPPfTQ6tWrH3vsMYqisrOzm9bvlAvcIwRmRQS+7sgfDkMe6vQZsErjOOqJql83GbEqAIDZ9O3bNzc3NyYmhqIoURSzs7Pd3NzkLQmCEJhV/Zn9jF8w4xPYlZPY90kQq8q4vFRjVQUAMKeDBw82Lsayf//+hx9+WPb1OCEIgRkRUnvwF8ehD3f1PBRl1yeh7mSre14CACzZgQMH6uvrL168uHXr1g8//FDuciAIgRnp0k4jSVKF9uj6qezi4hvOHSKcvuunAgCYEyEkNzf3hRdeYFn2iy++cHFxaXxcp9M1NDTU1NSIomjmkiAIgfnU7P/JYdgjHZo72Bra2Y0NiGi4dLzrpwIAmNO5c+eioqI0Gk10dHTTYNFr166lpKSsXLny2LFj33zzjZlLgiAEZsJfzRWuF2p6DTbWCe3ihtcn7TbW2QAAZsBxXHJy8j333JOdnX37446Ojn379q2vr3/wwQenTp1q5qpg+gQwk5p9P9oPHm3E7eY10QOqt38tVpbQrl7GOicAwKRUKtW0adPufryoqMjT05OiqPr6+oqKCjPPpoAWITAHqbZSl3rSvt8oI54Ts6wmZlD9yb1GPCcANg9jhJH5/hkoIyPj3Llz99xzT0pKivnnFEKLEJhDw4Vj6og4ys7BuKe1jxtR+cvnTvc/bZT7jgAogWeHVzY0h7Fjx8p4dWgRAnNoOHdY02Og0U/LBkYgioIJhQBYl6ZxoeYfINoiaBECk5Pqa7mCdNfxr5ri5HZxw+pO7TXKlAwAbB5XmFn62SKEzLbsNvZ69QPWN+iOGjju+++/379//5NPPuns7Lxz5861a9eaq56WQRACk2u4eEwV3guzalOc3C5uROmnC1z/9iJWmeT8ANgSVWCk19yNSBDMczmsUjNu3s0ePHv27Lhx4zZs2LBp0yaNRrN06dLr16/7+PiYp6QWQRACk2s4d1gTM8hEJ6ed3diAyIaLx4w7EgcAW0Wp7RFrrg7JlkaJDx48ODc3Nzw8XKPREEKys7M1GpnvW8I9QmBaRN/A5V5WR/Ux3SXsesOEQgCsycGDB0eOHIkQOnPmTP/+/ZsWl5ELBCEwrYZLJ9jg7pTaznSX0HTrzxVliZUlprsEAMCIDhw4oFara2pqPvroo08//VTuciAIgYk1nDus6THApJfALKvp3r/h/FGTXgUAYCyXL18ePXr00aNHP//888DALu1FYxQQhMCECKfXZ57TRPc19YXUkX10KSdNfRUAQNdlZWUFBwcHBASMHj3aweHm3OJLly69//77ZWVlv//+e0ZGhplLgiAEJqRLPckGRlJ2jqa+kDoihstPI/oGU18IANAVer1+9+7dgwcPzsrKuv1xLy+vgoICZ2fnQ4cOhYWFmbkqGDUKTKjh3GFN9/5muBBWaRhtmC7zvF2M0Rb1BgAYnVqtnjlz5t2P+/j46HQ6lmV5nmcYcwcTtAiBqRCB16Wf0XQ37Q3CJuqIOF3qKfNcCwBgdAMGDPjtt9/Mv9AogiAEpqNLO8P4BFMOzua5nCaqj+5SEiJmWzIDAOuEsdn+NW032K69e/cOHTq0oKBgypQpJn32LYKuUWAqJlpftDWMlz9maP5aPusXYraLAmB1aHP9bdoh/fv3z8vLe+yxx3x9fc1/dWgRApMgoqBLSTLPDcIm6ojeMHYUAGvk6urau3fvgIAAWa4OLUJgEvrM87SnP+3sbs6LqqPi6o7ucLr37+a8KABWJKM8e/rOOWa7HIXxloc3hroENXs8MzPT0dHRw8Pj7NmzMTExjo4mH1jeNghCYBL69LPqiDgzX1QV0rNy68dSQ60ZJmwAYI2i3MN/e+JbkZhprVGGYlzUzXtiKyoqcnJyPvjgg5kzZ953333PPvvsjz/+SFFydk9C1ygwCV16siY8xswXxSzLBkfr086a+boAAMMVFhYmJCSUlZWNHTvWwcGhqKjo6tWr8pYEQQiMT6qrFsuvMf6h5r+0OjIObhMCYMliY2PT0tJ69+6NMeZ5Pjc319PTU96SIAiB8ekyklUh3XFLO7CYmqZbn4aUUzCJAgBLdvDgwfj4eITQzp07x48fD9swARukzzinCjN3v2gj2sWTdnTiCjNluToAwBAHDhyora09cODA/v37Zd+eHsFgGWAK+vSzbs/Mlevq6sjeustJqqAouQoAALRBFMUrV67MnDmztLR01Khb+2nX1tZyHGdnZ2dnZ8Jd21oELUJgZEL5dcLpGS8Z1klqpI7srUs9LdfVAQBta1xWm+d5Ly+v2x//6aeftm7deuDAAfOXBEEIjEyfdkYVEYMMXlrJ6NigKKGkSKypkKsAAEBrOI4rKysbP358YWFhs089+eSTUVFRY8aMMX9VEITAyHQZyapQeW4QNsI0owqL0aedkbEGAECLVCrVk08++cQTT4SHh9/+eG1t7U8//TRs2LBjx46ZvyoIQmBUhOizzqvDe8pbhToiVpcOswkBaAGFKXP+M7CqoqKiESNGLF++vEePHiZ9+i2CwTLAmPirOZTGkXb2kLcMVUj32oO/yFsDAJbJWe0kdwktiI6ORgitWrVKlqtDixAYky49WW32BWXuxnj4IiKJ5TfkLgQAYAUgCIEx6eW+QdhEFRytz74gdxUAACsAQQiMhogCl5eqCo2WuxCEEGKDu+mzIAgBAO2DIARGw+WmMp5+FrLzgzqkuz7rotxVAACsAAQhMBp9ZrJcK6vdjfHSSrp6sapM7kIAAJYOghAYjS49WW0xQYgwVgV302dDoxAA0A6YPgGMg+gb+OI8NtCCVvhUBUfrsy/a9x0pdyEAyEOlUn377bc7d+6UuxAjKCkpSUhIMNHJIQiBceizzrMBEZhl5S7kFnVodOXPn8ldBQCyee655+699165q2iBIAiSJKlUqg59VVBQkInqgSAExqHLOK8Ok3lBmWYYn2CxukKqraQcXeWuBQAZsCwbFhYmdxUt4HlekiS1Wi13ITfBPUJgHFz2JXVwN7mruBPGbHA3ffYluesAAFg0CEJgBITT8TcKGT+L+9tTFdQNptUDANoGQQiMgMtPZ3yDLeoGYSN1KMwmBAC0A4IQGIE+N0UdFCl3FS1g/ULE8utSXbXchQAALBcEITACLucSG2hhNwgbUTQbGKnPTZG7DgCA5YIgBF1GCFeQrrLIFiFCSBUczcG0egBA6yAIQVfx1/IpeyfKwVnuQlqmDonWZ8J4GQBAqyAIQVdxOZdVQRbZL4oQQojRhvMlRZKuTu5CAAAWCoIQdJU+55JFrazWDKYZVhvO5abKXQgAwEJBEIKu4nJTVEGWG4QIIVh9GwDQBghC0CVidbmkq2c8/eQupC3q0B6wSS8AoDUQhKBLuNzLbFA3hLHchbSF1YbzV3OJwMtdCADAEkEQgi7R56SogiLkrqIdmFXRnr78lWy5CwEAWCIIQtAlXM4llWVOpb+TOiCSy0uTuwoAgCWCIASdR3iOv17I+lvcWtt3YwMiuDxYXwYA0AIIQtB5XH4a6xtkgWtt340NhBYhAKBlEISg8/S5ly184kQTxt1H0jeIVWVyFwIAsDgQhKDzuOzLljyV/g4Ys4GRXEG63HUAACwOBCHoLEK4gjQ2wELX2r4bqw2H3lEAwN0gCEEn8dfyKTtH2slF7kIMpQqE8TIAgBZAEIJO4nJT2GArmDjRhA2I4IqykSTKXQgAwLJAEIJO0udcUgVaTb8oQohS2zGuXvzVXLkLAQBYFghC0ElcXprKem4QNmIDwvVwmxAAcCfG8EMLCgr27t3r4uLyyCOPqNXqFo+5cOHCyZMnnZycEhISfH19jVQksDhSfY1UW8F4B8hdSMewAZFcXgoa/ojchQAALIihLcLjx4/HxcUlJSV98sknCQkJHMfdfcyiRYtGjx599OjRX3/9dePGjUatE1gWriCD8Q+38LW276YKjICBowCAZgxtES5fvnzhwoULFiwQBKFfv34///zzhAkTbj9g7969mzdvvnTpkre3twnqBJaFy09TBYTLXUWHMV5aqa5KqqumHJzlrgUAYCkMahFyHLdnz57HH38cIcQwzNixY3fs2NHsmB9++GHKlCm1tbU7duy4cuWK8SsFloTLT2f9rS8IEcasNpzLh0YhAOAWg1qExcXFkiRptdrGD7Va7bFjx5odk52dzXHcwYMHo6OjJ0+evGHDhokTJ7Z4tsrKypSUlNWrV9+sgGHGjx/v7+/f2tV5nud5Je4kZ8lPnCvIsH/gGROVJwiCIAimODNCiPYP02VfoiP7mOj8XcHzPEVRkiTJXYi58TxPCCGEyF2IufE8jzHG1naLoet4npckiaLMMVqTYZh2v8MGBWHjb2bTuSiKuvt9Sq/X19XVnTlzhqbpP/7449lnn50wYUKLz5PjOL1eX15e3vRITU2NKLY6u0sUxTY+a8Ms9omL5dcRQsjBxURv2Y1P3ETvDrRfqO70Po1lfmMt9RU3tcaXW4F50PjEFfiii6IoSZJ5njhN08YJQj8/P4TQjRs3goODEULXrl27uwHn7+/v6elJ0zRCaMSIERUVFcXFxU2NyNt5e3v36dNn3bp1hlwaIcTzvEajMfBgW2KxT7z+ep4qKLK1kcNdJ4qiWq020dsiG9a97rcvNCoVMstfox0iSZJarW78JVIahmEYpgOD2G1DY/yrVCq5CzE3mqYbf9rlLuQmg94LNBrN0KFDd+7c2fjhrl277rnnHoSQKIrFxcWNHRr33XdfevrNFY3T0tI0Gg2MmrFVfEEGq7WCPQhbRNk5Uo6u/PUCuQsBAFgKQ/8EW7JkyaRJk0pLS1NTU0tKSp555hmEUE5OTlRUVElJiaen56RJkz788MOpU6fGxsZ+8sknb775JmsN29SBTuDyUh1HPiZ3FZ3HBkVyeamsX4jchQAALIKhvUNjxozZs2ePKIqDBg06ceKEg4MDQsjX1/err75ycnJCCDk4OCQlJQ0YMKC2tvbLL79cvHixCasGMpIk7mouYw270reG1YZxealyVwEAsBQd6JTv169fv379bn/Eycnp2WefbfrQxcVl5syZRisNWCS+OJdx8aDUdnIX0nnqwMjKk3vlrgIAYCksbrwAsHBcQQartcIZhLdhvAOFyhJJVy93IQAAiwBBCDqGy09jAqy4XxQhhCia9Q3hCzPlrgMAYBEgCEHHcHlpamvbdOJurDaMK8iQuwoAgEWAIAQdQDidUHGDtrZNJ+7GasO4/HS5qwAAWAQIQtABXEEm6xuEaauf+Mxqw7lCCEIAAEIQhKBDuPw0q1xr+y6MmzfR66TaSrkLAQDID4IQdACXn85a+0iZRhiz/qFcAYyXAQBAEIKO4ArSrX3uRBNWG8YVQO8oAACCEBhMrK4g+gbGzUaWkGX9YeAoAAAhCEJgOC4/jQ2MRLayVw6rDYeBowAABEEIDMcXZKgCbKRfFCFEO7thmhYrS+QuBAAgMwhCYCh9XiprzWtt3w2m1QMAEAQhMBQh/JVM692GsEWMXygEIQAAghAYRCi9SmkcKQdnuQsxJpU2jMtPk7sKAIDMIAiBQbiCDMa2moMIITYgnC/KQoTIXQgAQE4QhMAgXGEm6xsidxVGRtk5UnYOQulVuQsBAMgJghAYhM9PV9nGmjJ3YrThMK0eAIWDIAQGIIQvzmX8QuSuw/hYv1AIQgAUDoIQtI+/XkA5ulAaB7kLMT6VNpzLh4GjACgaBCFoH1+QYWMzCJuw2lC+OA9JotyFAABkA0EI2scVZTL+oXJXYRJYpWFcPPjifLkLAQDIBoIQtI/Lz1DZ3NyJJkxAOFcI+zEBoFwQhKA9ksRfy2dtcaRMI9Yf9mMCQNEgCEE7+OI8xsUTqzRyF2IqKv9Q2IYCACWDIATt4AozGVucQdiE8Q0WSq4QnpO7EACAPCAIQTt4W1xT5naYYVkvf/5qjtyFAADkAUEI2sHlp9vwSJlGjH8ojJcBQLEgCEFbiCjwNwoY32C5CzEt1j+Uh/2YAFAqCELQFv5qDuPhh1mV3IWYlsofdugFQLkgCEFb+IJM1kan0t+O9g4QKm4QTid3IQAAGUAQgrZwhYoIQkwzrJeWvwLjZQBQIghC0BauIIO19ZEyjVgYLwOAUkEQglYRnhNKrzLegXIXYg6MHwQhAAoFQQhaxV/JYb20mGHlLsQcVFoYOAqAQkEQglZxBemM1vZvEDZiGsfL6BvkLgQAYG4QhKBVXGEm66eUIEQUzfoGcjBeBgDlgSAEreILMlSKaREihFi/ML4QekcBUBwIQtAywumFihu0V4DchZgP4x8K0+oBUCAIQtAyriiL9QnENCN3IeajghkUACgSBCFoGV+YwforYgZhE8ZLK1aVSbp6uQsBAJgVBCFoGVeQwShgTZk7UBTrG8wXZctdBwDArCAIQcu4wkyV0oIQIdYvhCuC3lEAlAWCELSA6BvEqjLGy1/uQswN9mMCQIEgCEELuMIsxicYUbTchZgbA/sxAaA8EISgBVyhsmYQNmG9/KXaCqmhVu5CAADmA0EIWsAXZChoTZnbYcz4hMB4GQAUBYIQtIArzFTOKqPNsFqYTQiAskAQguYkXZ1YU856Km6kTCPWPwTGywCgKBCEoDm+MJP1C0MYy12IPFi/MA5WHAVASSAIQXNcQabKP0TuKmTDePpJdTVSfY3chQAAzASCEDTHFSpvTZnbYcz6BfNFWXLXAQAwEwhC0BxfmKm0VUabYWH1bQCUBIIQ3EGqr5XqqhkPX7kLkRPjH8rlw21CAJQCghDc4Wa/qFJHyjRiteGwQy8AygFBCO7AK3Kt7WYYN29JVy/VVsldCADAHCAIwR2UuPvS3TBm4DYhAIoBQQjuACNlGqm0sPo2AEoBQQhukeqqpYY6xs1b7kLkx/qHwm1CABQCghDcwhVkMFrlrilzOxb2YwJAMSAIwS1cEYyUuYl29SSiKFaVyV0IAMDkIAjBLXxBBgtB+D+sNgzGywCgBBCE4BauIBOGjDZRaUN5CEIAFACCENwk1VYRTse4esldiKVg/MK4gnS5qwAAmBwEIbiJK8hgA8JhpEwTVhsO+zEBoAQQhOAmrjCT9QuRuwoLQju5YIoRy2/IXQgAwLQgCMFNXH4a6x8udxWWhdXCJr0A2D4IQnATX5jJBsCaMndg/GChNQBsHwQhQAghsbIEIUI7u8tdiGVRaUN5mFYPgK2DIAQI3RwpEyF3FRaH1YZzRRmIELkLAQCYEAQhQOjmphPQL9ocZe+E1fZCWbHchQAATAiCECCEEJefptJCELYAFh0FwOZBEAKECOGvZMHiai1i/UP4oiy5qwAAmBAEIUBCyRVK40TZO8ldiCVSacO5/DS5qwAAmBAEIUBcQQYL/aKtYLRh/JUcGC8DgA3rQBBWVlbu37//4sWLbR9WUlKSn5/ftaqAWXEF6bDWdmsotT3l4MzfKJS7EACAqRgahCdOnIiKilq3bt3YsWNfeOGF1g4rLy+Pi4sbMGCAkcoD5sDlp8NImTaw2jC+AKbVA2CzDA3ChQsXzp8/f8eOHWfOnPnjjz+SkpJaPOy1114bO3as8coDpieJ/LV8GCnTBsYvlCuEbSgAsFkGBWFpaWliYuKUKVMQQu7u7g8//PC2bdvuPmzHjh2lpaUTJkwwco3AlPjifMbFA6s0chdiuVTaMC4fZlAAYLMYQw4qKirSaDTe3t6NHwYFBWVmNu8pqq6unjt37o4dOwoKCto+W0NDQ0FBwY8//tj0yH333efq6tra8ZIkSZJkSJ02xjxPXJ+fxmjDLOo73PjEscVsCEX7hfDX8kSew7RBvy+dZmlP3Gyk/5G7EHNrfLmV+cTN9opTVPvtPYN+sXU6nUqlavpQrVY3NDQ0O2bevHkzZswIDQ1tNwjLysry8/N/+OGHpkeCg4N79erV2vF6vZ5lWUPqtDHmeeINuamUd5Berzf1he4mElRYi7JqcXYNlVWDyvSIl5BOxA08yyOJxshbQzw1yM+OeKiR1p7EuhEfmRqulKtXXV46ozXt7hw6nY4QQtO0Sa9igXQ6HcMwDGPavzMskF6vV2YQ8jwvSRIxy2BsjUbTbhYa9JPn5+dXU1PDcVxjHJaUlPj5+d1+QGFh4bfffvvyyy8vXLiwsLCwrq5u4cKF8+bN8/T0vPtsAQEB8fHx3333nYFPQxRFe3t7Aw+2JeZ54jVXsp36JrB2dqa+UCNeQqdLyP5iKbGYZFUTbzsc6oRCHXGsJ3bXYBVGGhoRQe9oRxOErzegch0p0aHsOpJ4gyxKJvYM7u+B+3vhAV64l7v5mk5cQAR1I98+stU/14xFrVYrMAgpilJmENI0jTG+vZmhEI1BqFar5S7kJoN+8gICArRabWJi4v33348QSkxMnDNnzu0HODg4LFu2rPH/lZWVGGM3NzdDGqRAXoTnhNKrjE+gqS9Uw6PthdLeK+TINSnICQ/zxm/E0tEuWNXSe35DA7HTYIRxpDNC6I6wK6gh5yvIhTLyrwxJL6G/BePHQqgeriYPRFYbxueloWGPmPpCAADzMygIaZqeM2fOrFmzli9ffuzYsfLy8ieffBIhdOLEiZEjR+p0Ond39wULFjQenJiYuG3btqYPgSXjr2SzXlrMmLAD9nIF+SpT+rNAGuhFJfjhub1Y9y78FRjkhIOc8KNBCCGUWU12FUrPJYp2NBoXTE2MwD52pkpEVUB4fdJuE50cACAvQ/si5syZ4+Pjs3PnTn9//yNHjmg0GoRQQEDA4sWLmx0ZFhb25ptvGrlMYBqmW1OGE9FvBdLXmdK1evR4CPXzfayHsVBiNiMAACAASURBVHtBIp1xZE/65Z7oQjnZUSiN3C6OCaRmdafCnI0fh4x3kFBZKunqKY0Se+kBsG2GBiHGeOLEiRMnTrz9wYCAgLfeeqvZkYGBga+88opxqgMmxuWnsYFRxj2nSNBPudK6i1KIA3o2khruQ1Em7rmMdcex7vTMHvQP2dK4PcJgH+rlHlScu1GvSlGsbwhfmKGO7G3M0wIALIDi7k6D23EFGQ6DRxvrbAShnYXS2vOSixqt7kfHeph1GoAzi6ZHU5MiqF/zpWmHxJ6uaHlfOtR4rUNVQBiXnw5BCIDtgSBULklXJ1aXM17+RjnbiRtkebLIiei1GHqoj2wz4ewYNCGceiqU+k+O9Oge4ekw6rUY2tEY90AZbZg+9bQRTgQAsDAwsFO5+MJMlV8ooro6WL+aR/NOijOPis+EUd+PYmRMwSYMhSZHUFvvZa/Uo/jtwtZcqevzldiACC4P9mMCwAZBECoXV5DBaru6xOj2AilhuyBK6Kf7mAcCLOvHyV2NlvWl3x9Ib06Xxu0R8mq6lIaMqxciRKwqM1Z5AAALYVnvXMCcuPyMruy+dL2BTD0kvnteWtufXhhHO1hqL3tPN/x1AnOvH/XIbuHbrC41DZmAcC4fVt8GwNZAECoXl5fCBkR07mt3FUn37xRCHdF/RzFx5h0U0zlPh1Ob45lvMqTJB4Qbuk6ehPUP5QogCAGwNRCECiVWlhJRZNy8O/qFvITePiMuPSN9MIh5sTvNWs9yYCFO+KsEJtIZ37+D317QmdUdVXCbEABbZKn9WcDEuPxUVVBkR78qv5a8dFT0VOP/jGKcrHAhdIZCL/Wgh/tRS06LJ0rI231opiN/CrLaML4oExGClLdBBAA2DFqECsXlp7Md3EvhzwLpkd3CgwHUukG0NaZgkxg3/N1IJrOS/H2/UNqRblLKzpFydOGvF5qsNACADCAIFYrLTVEFGnqDkCC09ry4MlnaOISZEGYLPzOOLFo/hIl1w6N3CcllHRhAw2rD4TYhADbGFt7UQEcRUeCu5jL+Bq0yqhfRrKNi4jXyVQLT3fT7PJgNxmhGD3phHP3sQeG7LENvGbL+YTwMHAXAtkAQKhF/NYdx86bU7e9BeL2BjNsjSARtGsZ0ZdcIixXvg78cwXyWKi0/KxoytUIVEMHlw3gZAGwKDJZRIi4/nTWgXzS1kkxJFB8NxC92t56xoR0X7Ii/TmDmJAnTj4ifDKE1bf5OMH7B/I0iwnOYVdxmql10re5GdkXe1dprZQ3lZbqKyoaq0vqyBkGHECKEIHxzAJKTxtFd7epq5+pl5+GmcQl01oa7hnrYuclbPLBtEIRKxOWltjuDcO9VMueEsCSOHuVv+90Gjiz6dAjz5lnxmYPCv0YwLq1nHGZY1kvLX8lWhXQ3Y4FWqayh4nRxcmp5RnZFXnZlnobWhLoE+Tn6uGqcu7lFuPg5uWlc7RkNQkjPcTRNMzSNEKrl6iv1VVW66nJdZXZlXmLBsdyqAoxwuGtIlHt4d8+ofr5xTipHuZ8csCkQhErE5aXZD36ojQO25UnLz4jrhzAxbrZzU7BtLI3e7U9/cEkct1v4bhStdWj1ibMBYVx+GgRhi3SC/kLJ5VPF505dPVuqK4/z7hnhFvZ41COhrkFtpJee0dM0zTAMQsjHoYUDynWVeVUF+VWFf2btWntiY7CzdqB/3/6+vXt6RdPYlrsrgHlAECqOVFct1VayrW86sSVd+jhF+jyeCXNSSgo2whjN7UV/ly2N3SP+eyTd2sgg1j+cK4DbhHcQiXjyavLuvAMnrpwOdwmJ9e45o+/UMNcQykgTLt01ru4a174+sQghXhIyyrLOlVz+6MymG3WlCUHDHgwdGePVHSNl/bgCI4IgVBwuP53Rhrc2Jfz9S9JPOdLmeNq/9SaRbZsYTnmp0dP7hW9HMrEt7e7LBITXHvrN/IVZppTSjN15Bw7kH/Z39IsPHPRsr/FOrGn7LVmK6ekV3dMremKPJ8oayg8XnliX9EmDoHsgdNQDoSODnANMenVgkyAIFYfLT1W3tKYMQWjFWXH/VfJlPO2pUWgKNnoggLJn8aSDwpZ4pr9X828F6+kv1deINRW0k3JHcAiSeKDg8I+pv9Xx9QlBQ9aMfMvb3tP8ZXjYuf8taszfosbkVRceLjj+2p4lwS6B43v8bZB/P2ggAsNBECoOl5tqP/C+Zg9KBM1NErOqyZYRjFG2sbV2w33wir7M1MPCpmHMkGY7LGLMBkZweWl2vYbIVJ2c6vj6HVl7f0j7xcvOY2zk6MH+/Sgs/3CqEOfAkJjACT2eOFl8dsv57z86vemJqEcfibhfw2jkLg1YAQhChSGEK8pwefylOx5DaNFpMbOafDyYsYMU/J8hPvi9AcyLR4WNQ5kE3zuykA2I4PJSlRaEVfrqf1/+aUf2noH+/RYNfS3EOVDuippjKHqodsBQ7YDLpWnbs/d8e/mHp7qNezL6UYhD0DYIQmXhrxdSdk6UvVPTIxJBb5wUM6vJx0MYO/hxuFNfT7xuIPPyUeHDIcx9/reyUB0UVZuooNuEdXz9D6m/bsv4Mz5wyIb7VrtqXOSuqB09PaN7ekZfrbn2Q9qvE37/x6Sefx8b8SBLw195oGXwzqcsXH7q7VPpCULzT4kZVWQjpGAr4jzwR0OY2ceFj4fdaheyARHc1RwiCpi28e+aTtD/krH9v6nb+vn2XjdqmZccNwI7zd/Jd86Al/KqC/+bsu2/qb9M7TXhwbBRMN0C3E3+zn1gTlx+Oqu9GYQEocWnxdQKsgFSsE093fAHg5mXjwnHrt9chA2rNIy7L38lW97CTO1gwdHJf864cCNlZfziWX2nWVcKNglxDlw4ePbrA17ambN32o7Z525ckrsiYHHg/U9ZuNwUl9ihCCGC0MJTYkoF+QRS0ACx7vjdAcyLR4WvRjD9PDFCSBUUweWmqIK6yV2aSRRWX/no9Kbi2uszek+N8+kpdzlGEOUevmz4/NPXklcd/SDcLWTOgBk+Dl5yFwUsBbQIFYRwOqH8OuMThBBafU48V0o2DIHRMYYa4IlX9GWmHhIulhOEEBsQqc9Nkbso49MJui0Xvn95z4JYrx4f3rfSNlKwSX/fPhvufzfIOfD5HbO3nP+eF3m5KwIWAYJQQbj8dNY3GNPMxhRpVxHZOIyxh7ZgRwz1wUt605MThbQqwgZFcjmX5a7IyE5cPT3pjxmF1Vc+vPedhyPut8nbaSqaHd993NqRb10sTX3hrznpZVlyVwTkB2+ECsLlp6kCI/6dJf07U9oygnaBtmDHjfSj9AKacED85V5veyKJFTdoN2+5izKCOr7+s7P/Oll8dlbf52O9e8hdjsn5OnovHTrn2JVTCxJXPBR677TYZ2BMqZJBi1BBuLzUi5qwdRfEj4cofe2YrngwkJrejZpwQJS0NtI7evLq2ee2v6wX9R/cs1IJKdhkqHbAB/eszK0qeGHHa2llmXKXA2QDLUIFqc1Ne1P8+8ZRTJDCVtM2uidCqHId+v5M+JTsFPu+I+Uup/Pq+YYNZ/4v+dqFl/u+EOMVLXc5MnBRO70xaNaRwhMLDq54NOLB53pNYCgb7BAGbYMWoVKcTiuq4PCiET5RLpCCRjA9miLayMwLKXpR7lI6K70s6/mdsyVJ/ODelcpMwSbDAwd/eO/KtLLMl/fML669Lnc5wNwgCBUhu5ps2nVRCuzeq6XtFEDnTBwe4V9fOHlPnUjkLqWDCCJbU39/48Cyp7s/9mKfZzWMWu6K5Oeidl40ZPb9wSNf2jVvT95BucsBZgVdo7avRIce2iVuYVN9Qm1z0ptcMMOyvkEepZmzjvX6fJjV9KdV6qtWH11fqa9cM/JNmEvXzMjgYWHuIetPfp509ezcgTPtYJFSZYAWoY2r5dFDfwnjgqiQsks4AILQyCht5Hy39GPXyerzkty1GOTstQvTtr8a7Br4TsJiSMEWBTlp301YKhFx+s45uZX5cpcDzAGC0JaJBE08KEa74NfDalB9FfKyuO0CrB0OiGQLUr5JYD5Pkb7NsugsJIh8f/nnlcfWze7/4tPd/2aTcwSNRc2oZ/V9/vGoh2fvXbIv/7Dc5QCTg65Rm0UQmn5YrBfIx0MYknoJa7u1tis96Dz/KPTXFh8N+XYU/dReQWuP7/G3xG9yg6B79/j6opqrqxOWyrKDrjVKCBoa7Br43vGNZ6+dnzNgBowmtWHQIrRZ7yRL58rI58MYhkIk/zIKiJK7Ilvk5IZVGlR2JcoZfzyUmXBASKu0uJEzhdVXXvprLkJ41YglkIIdEuIcuO6e5aUN5bP3Li5rqJC7HGAqEIS2aWuu9EWqtDmevrmIWu5FBDcITUQbJRWkIoTiffHbfekxu8TrDXKXdJvDhSdm7V7wcPgDs/tPV8HiKR1nz9rNGzirp2f0P/56PaU0Q+5ygElAENqgUyXk5aPiNyNpX3uMEEKcjpQUYb8wueuyUdpIVJDa+N/HgqnHgvHDu4Q6Qd6abvr3pa0fnfpi6dDX7w2Jl7sWK0ZhPL77uH/0nrwoccW+vENylwOMD4LQ1uTUkLF7hPVDmO6u/7tZVZiKfYIRA60Bk8ABUSjv1urbc2PpcGf8bKLMcwt5kX/3+Pq9+YdWjVwS7hYiay02or9v72XD5286980XyV9LxOI6wEFXQBDalGoejd0tvtqTHul3a8gGybuMAqFf1GS8AknldaSrbfwII/TPgfT1evLmadmWnKnW18zd/1a1vnb1iCWedu5ylWF7gpwD1o5668KNlGVH1uoEvdzlAKOBILQdvIQe3yPE++LnIu94WaXcS1gbKVdVto+ikTaC5N9afZul0BfxzH9yiCwTKnIq86f/NSfaI2ruoBlqRmX+Amybk8rxrWGv05h6Ze/C0voyucsBxgFBaDtmHhNZCi/tfecgb0nEVzJgyKhJ4YDuJPfi7Y+4qdCWEfS8E+LR62btQztdfG7OviXP9HhifPdxGFniRA4bwNDsy/1eGOTXb8buN3Jgxr1NgCC0ER9dkg4Xk4+G0tSd737kajZy8UQaR5nqUobgaJRzodljUc74o6HM43uF7GozZeGunP0rj74/b+Cs4QGDzHNFJXssasyzMU+/tnfJyatn5a4FdBUEoS3Yc4WsPid+lUA73b1AQt4luEFoatg/gtwoQFzzaRMJvnh2T/qR3WIVZ/Iavkv5afOF71aOWNjDE15uMxmi7b9wyKurT6z/M2uX3LWALoEgtHpplWTSQeH/4plgxxa6wkjeJaSFflETo1nsF3r7bcImU6OoId54wn7BdKNIJSKtO/np/vzDa0a+qXXyM9VlQEu6uUesil/83eWft1z4Xu5aQOdBEFq3Cj16dI84P5Ye4NXKDaGCFAwtQjMI6kHu6h1ttKIfzRG08JRJBpHqBN38A8uLa4tXDl/oonY2xSVA23wdvVeOWHSk8MQ/kz6RiEWvNwtaA0FoxQQJPblXeECLnwlv5XUsvUIYFjnDqlomhwOjUc75Fj9FY/TxEGZbHtmcbuR3yWp9zWv7lrhqnBcMflUNewrKx03jsnLEois1xW8dXsOLvNzlgA6DILRis0+INIUWxbW6FjDJg62XzEUbSa7nI07X4iddVOhfI+hFp8QTN4zWQ1rWUDF77+Lu7lEv9XkOtpKQnYZRLxzyKiJk7v636/h6ucsBHQNBaK02pUn7rpCPhzF064PkSd4lDDcIzYNhsU8IKWjhNmGjSGf8/mDmiX1iUZ0RsrCguuilXfPiA4dMinkKpklYCJZiXhv4UrBL4Gt7l1ToquQuB3QABKFVOnadvHVG/NeIloaJ3obkXSJwg9BsgqJRzsU2Pn+fP54WRT26W6zv2kqk6WVZs/cufqrb2HGRo7t0ImBsGOHneo0f5N9vxq55V2qK5S4HGAqC0PoU1JIn94kfDWFCnNpqCpDaCtRQgz0DzFaYwuGgHiS35duETWZ2pyKd8bTDnR9Deuba+fkHl83sOxXW0bZYj0WNGRv54Oy9i3OrCuSuBRgEgtDK1Ato3G7xpe7UCN/2OsSyz+PAaNiM13y0kaQ4F3HtLEH53kA6s5KsPd+ZgTNHi06uOPrPNwa93NcnrlMlAjN5MPSeyT2fmrNvaXpZlty1gPZBEFoTgtC0w2KkC57erf0XjmSfQ8E9zVAVuIlVIe8gUpTa9lFqGm0ewXx8WfqzoGPNwn35h99L2rBw8GvdPeC+rxUYHjh4Zp+p8w68feZaO/0EQHYQhNZkdbKUVUXeG2jYEMHsZBQSY+KKwB1wUPe2bxM28rFDnw+nnz8sZFQZmoV/ZO365PSXbw+fH+UO+0pajX6+cW8MmrXi6DpYhs3CQRBajT8LyKep0pfxtNqQHKy4gTg99vA3eVngdoHR7d4mbNTPEy+IpcfuMWj1tZ/S/vj20o8rRywKdoY7vlamp2f0kiFz3j2+/mDBUblrAa2CILQO6VXk+cPCpuG0j51B9/xIdjIKjYEbhGaGA6LIlWzEG7S06IRwapg3/vu+dlZf23Lh+18zd6wasdjX0ds4VQLzCncLWTLs9Q2nNv2Vs1/uWkDLIAitQA2PHtsjzo+l+3oaGmwk6ywK7mHSqkALVBrsHUAK0w08fFlful5Ab51pdfW1zee/O1R4fOWIBe52bkYqEcgg1CVoefyCL899uy19u9y1gBZAEFo6iaCn9wtDfVpfR+1uhKCc8zBSRh6B3Q3sHUUIMRT6fDjz7yyyNbf5IFKCyMYz/3es6OSy4W84q2ARUavn7+T7TsKiH1J/+THtN7lrAc1BEFq6pafFMj16q08H1tAiN/KRSo1dvExXFWhVUDTKbX+8TBN3Ndoygp5xVDxbequHlCCy/tQXKSXpy+LnO6lgL0kb4WXvuTx+wS/p27+6+F+5awF3gCC0aNvypO+yyZfxDNuhFyrrHAqG8aLywAHRpCgDCR1YebmHK147gH58n1iqxwghiUhrj2/ILM9eMvR1e9bOZJUCGXjZe7wzYtH+/MNfJH8tdy3gFghCy3WpEr94RNwcT3t0cF8Bkp0MNwhlo9JgT39SlNGhLxoTSD0SgKccofSiuOLouut1JW8Om2fHakxUI5CRq8Zl2fD5x66c/L9z38hdC7gJgtBClevRxKOqlX3pnm4dHPkpSST/EoYbhDIK7kmyOjxvbFFvmkbihB3rarm6hYNfUdGsKUoDlsBF7bQ8fsGxK6c2nYcstAgQhJZIJOiZ/cIDfuK4kA6/QORKBnL2RPYwvEI2OKQX6ngQSkTorn7/en19sNcMBlLQ1jmpHJfHLzhfevmjM5sIMtrmXKBzIAgt0bwkkSfotehO7WkOK6vJLqgbup6PGmoM/wpBEtafXENjacGgl9+9gM+Xwzuj7XNg7RcNmp1ekf3+yc8kAq+4nCAILc63WdLvBeSToW1tNNgGKTuZCoEglBXFoMBuJNvQSRSCyH9w4l1C0KQez0c600t709MOiTda3uIX2BR7xu7NIXNyK/PXnfwEslBGEISWJbmMzE0SN8fTLqpOfT3PocIMEhBt5LJAR4XGkswzhhzIidyaEytYmp3W+yUaUwihkX7UuCD8/CGB61SPALAualq9aPDs3MoCyEIZQRBakOsNaNwe8b0BdDeXTi6NRgpTkHcgVsOYe5nhsFiSfrLdw/SCfu3x5Y4qpymxL1D41i/ji91pFxYvaX3FGWBL7FjNm0Pn5FcVvpe0EbJQFhCEloKT0ON7hfGh1IMBnX9RSNY5GC9qEdz9ME2j0qI2DtEL+vdOrHBgHSfHTLs9BRut7E8n3SBfZ3Zm20JgddSMesmQOVdrilce/adI4A8gc4MgtBQvHxNdWDynV9dekaxkDFsvWYjgGJLRau9oYwq6qF2fjf3H3SmIELJn0PrB9IeXxGPXoYmgCGpGtWjIa+W6ypVH34csNDMIQovwcYqUWEzWD+3c+Jj/4RrIjQLkH2GsqkCXhMZKGada/Ixe0K09vtxV7TYxZhrV+g4hWge8oi/z0lGhoBayUBFUNLt48Oxarm7F0XWQheYEQSi/xGLyzlnx6wTaienaiXLOI20EYmAKmkXAIT1xfsrda601CA2rjr7lae81sdfUNlKw0SBvPCWSnn5E1AmmqhNYFIZm5w2cWaWveQfahWYEQSiz/Foy4YC4cRgT7NjVvQNJ+ikc2ssoVQEj0DgQryCUf/n2x+qF+tVH3/Jz9H+65xSMDHrFJ0dQEc745RNt71oIbAdLswsHvVKlr4F2odlAEMqpQUCP7xFn9qCG+xhhB10p/SSO6Nv18wBjwWEx0m2TKOqF+tVH3vJ18Bvfc7KBKdhoYSxdVEs2psDAGaVgaXb+4Ffq+fplR94TJMhCk4MglA1BaEqiGOWKn48yxqtwLRdjGnn4G+FUwEhwSCzKON34/8YU9HfSToh5tkMpiBBS0ej9QfS/0qU9V6BZqBQsxcwdOIsTuGVH1kIWmlrH3oLLy8uzs7MlqdW/TG/cuFFQUCCK8LK1b/lZKb+GrOnfgY0G20DST6KI3kY5FTAavzBSXUpqyuuF+tVH3tQ6BRjeI9qMpwavH0y/niSkVUEWKgVLMXMGzuBFfsmhd3ixAxt7gY7qQBC+/fbbERERjz76aPfu3bOyspp99uLFi1FRUT169Bg1alRQUNDu3buNWqet+SVP+r80aVM8rTJODiIpNQmFQRBaGIrCwT3rMo6vPLIkxDWi0ynYKNoVvx5DP5soluuNWCKwaI1ZiAh6+8haXoIRU6ZiaBBeuHBhw4YN58+fT0lJGTt27BtvvNHsAI1Gs2XLltLS0uzs7MWLF0+aNInAEgmtOFdGXjwqfpVA+9gZ4dYgQgg11KDreTi4u3HOBoynNqTbO/n/jnCNeiJ6fNfPNiaQGuWLZxwRBLhdqBgsxbw+cKYkSW8eehfahSZiaBD+5z//eeSRRwIDAxFCM2bM+PPPP2tq7lhcPzIycvjw4Y3/Hz16dGlpaV1dnXFrtQ03GtDf9oqr+tExHd1osHUk/RQO6Ylg7x4LUyfUvdtwOKJC93i3vxvrnLNjaILQ22fh7oOCMBQ9d9BMFcUuSnyHEzm5y7FBhgZhXl5eeHh44/9DQ0MRQkVFrS4f9c0338THxzs6Orb4WVEUy8vLz9xGr1dKXw8voaf2C0+G4EeCjDlMiaSdQOF9jHhC0HXVfPWK0+9GunV7vM4V3cg31mkpjNYOZI5cI19lQKtQQShMvdL/BQfWbnHiKj1kobEZOoW7trbWzu7mUs4YY41G06xF2OSPP/747LPPjhw50tqp8vPzk5KSpk+f3vTIunXrBg4c2NrxdXV1uL15x9ZiZhJjj/GLYYIhreX6+npDnjiWRFVmMj9gLKmvN0KJFqBBpyOSZNUvejVf88/LH0U7R472vY/3q0dpp0Qn73a/qkGnE1iWptu5b4wRWhWHZyYxASr9EC9biEO9Xs8wTLtP3PZwHIcxZllD+3KmRk/YdOnfr+9ZumLIAjWtNmltJsXzvCRJPG+Onl57e3uKaqfhYWgQent7V1ZWNv6f47i6ujofH5+7D9uzZ88LL7zw559/RkVFtXaqsLCwhx566LvvvjPw0oSQ1hqX1uWfF6RzldIv9zP2jEE/wYQQBweH9g/LvUhcvTTeNjRxAmM7jQZZbRBWcTXrLmzo6d7jsdBHEEIkPBad2qlKeKrdL8QUpWJZyoA8iLRHawaSuafRb/czIU7W+o1q0piCDNPFpZWsD8uyHQpChNCcQS9+cnbLspP/fG/UMjtGY7raTKoxCNVqS8lyQzvoevfufeLEicb/JyUleXl5abXaZsccOXJk0qRJW7dubaN5p1g7C8n6S9K/Emh7o/+ypyUhmEdvMSr0lctPr+rr0bsxBRFCOCCSlBej6nLjXqifF365Bz05UayG8RNKQmFqVt9pnnYeCw6u0AmwfbNxGBqEkyZNunDhwvr1648fPz537twZM2Y0/vn20ksvffTRRwihixcvjh49+qmnnrp+/frWrVu3bt3aWt+pAp0vJ88eEjaPoP3tjf/HO0lPQuEwccIilOnLV5x5t79334eC7r/1KMXgkJ4ky6B9ejtkbDA12Bu/CINIFYbC1Kx+07SOvnP3v1XH28gNEXkZGoRubm779u07duzYggULHnvssaVLlzY+Hh0d3TiUlOO4iRMn8jy/93/qbeWWVRdda0Bjd4ur+tFx7iZIwfJrpK4a+4Ya/cygo0p1ZSvPrBnuO3RM4APNPxfWG6W3vBNFF70eQ4sSehMGkSoMRvj5uIkhLkGv7V1SzUGTo6uw+Wf7ff/999u3bzf8HmFNTY2Tk5NJSzIdnYhGbhdG+VGv9uzwMNHa2tp2b46SY7+RwjQ8Znrbh1mXhoYGq7tHWKIrfefM2pHa+FF+I1r4tMBJmxdRszYgTVs3fXU6nYH3CG9XL6Bph4Rnwqjp3a11xUS9Xq/Me4QdHSzTDEHk6wv/TSnL+PDed5zV1vQmaa33CEEnEISeTRQDHPArHU9BQy+RBv2i8iuuu7bs1Kp7tCNaTkGEEKPCAVEoO9kUV2/cwvfTNGnPFeghVRaM8HOxE/r4xL62b2mlvkrucqwYBKEJLT0t5taQdQO7tt1uGzgdKUjFIT1NdHpgiMK6KyvPrn0k+KEEv/g2DsPhcSTtpIlq8LXH7w+iX08SL1XAck6KM6HHY318es3eu6S8oULuWqwVBKGpbMmQvs8im+MZtckmR5GsZKSNQGp7U10AtCevJn/V2ffGhYwZ4tPeSOnQXijvMuJNNRW6pxteHEdPSRSv1kMWKs4zPR6PDxw8c/f84trrctdilSAITWL3FbL4pPjvkbSnKef5kLQkDP2i8smpzltz7oOnw58Y5D2g/aPVDsQnGOVeNF099/hTT4dRzyWKdTChQnkeixzzaMSDs/cuvlJTLHct1geC0PgulpNJB4UvhjNhzqYc7iGJJPU4jjLgLRiYQFpF/sTnWwAAIABJREFU+prk958Jf6q3R6yBX4LD4qQMU/WONno2kurmgl86ChMqlGh02D1PdHv4lT0Lc6sK5K7FykAQGllxPXpkt/hWH3qgt2kHPZKss8jND7l4mvQqoEUpFWkfXvxkWvTkWI8Yw78Kh8ehrGRk4k1WF/emRYLnn4IJFUp0X8jIKTF/n7vvzayKXLlrsSYQhMZUy6OHdwsTI6gnQkz+jSUXDuHug0x9FXC3UzfOfnTx0xd7TIt2bXUdwZY5uGIXT1SYZpq6bqIwWt2PPl9G1l+GVqESxQcOea7XhHkH3k4ry5S7FqsBQWg0vISe3CvEuOJXepj+uyoKJO0E7gZL2Znb0WvHv0z714wez4c5hXTm68PiiGlm1t/OjkUbhjDfZUk/5UEWKtGwgIGz+k5948CyM9fOy12LdYAgNA6C0D+OiASjVf3NsYI+ST+FvIOQk5sZrgWa7C3a/+/MH2bHzAxxCu7cGXB4H5J+Epl+FQtPDdo4hF5xRjx0DQaRKlFfn7j5g15efuSfhwuPy12LFYAgNI4lp8TL5eSzoQxjlu8ouXAQRw82x5XA//yet/23/O1zY1/xd/Dr/FlcvZFKg66Z4/5NqBNeO5CZdUxIrYQsVKIent2WDnv9/ZOf/ZWzX+5aLB0EoRF8fFn6KY98NZKxM88SUZyeZJzG3WC8qPn8lPNrYvHhub1e9dR4dPFU2Cy9o436euL5veiJB8XCWshCJQp3DVk+fP6X5779Of0PuWuxaBCEXfV7vrTmvPTvBNpNZaYrkoyTWBuB7J3NdD1lk4j0fyn/Si5Nfj32FVe1ixHOGN6bpCcZ4TyGeSCAej6KGr9fLIEdexQpwNl/xYgFP6b+9l3KT3LXYrkgCLvk0DUy/bD49Ug60NF8K0ST8wcR9IuahUCEjZc+u9pQ/ErMTAem/U2SDYG9gxFCqDjbKGczxJOh1INaPPmAUAsT7RXJx8F7dcKSPTmJH53+P8nsuyxYBQjCzjtbSp7cJ3wynOnhasZ9EvQNKCsZR/Uz3xWVSifq30v+UCdyM3tM19DGXCYfdxuILhwy4gnb9VIPOtoVT0sUOJheqEiuGpcVIxaklWWsPv6BYOKZrNYIgrCTMqvII7vFVf3oYSaeON8MSTtOgrojTTvbM4EuquKql59e7aFxnx79HEt1cpec1uBuA6XU40gwawNtYRztoMIvHhVFaBIokgNrv2zYvFqufkniOzpBL3c5lgWCsDMKasl9O8UFcdQjQeb+BkrnEmEevandaChZdmpVL/ee48OfoEyxLaKTO/LUomyzzvHCGK3sS1fqyeLTEIUKxdDs6wNfsmft5+5/s4arlbscCwJB2GE3GtD9O8UXulF/DzX7d09Xi/IuogjoFzWhnJq8ZadW3aNNGBN010bzxoOjB5GLiaY7f4tYGn0wmDlXSt5Jhs4xhaIxPavftEi3sFm751+vuyF3OZYCgrBjqnn00C5hbBD1QjcZvnXk8jEc2gur7cx/aYW4VJ6y5uz7T4b/Ld5vqEkvhCP6kMI0VF9t0qvczZ5Bnw5jDl4j71+CRWcUCiM8Oebv94eMnLF7fka5+QZtWTIIwg6oE9CYXUJfD/x6L3m+b+RCIo6GflFTOVR8dOPFz17sMbWvp+k3t2LUODQGpciw6ocjiz4ewvycI32eClmoXGPC73shduK8/W+fvHpW7lrkB0FoqAYBPbpLCHTAy/qaYxG1FlSVkCsZKKKvPFe3db/nb9+as21O3CvhzmHmuSLuNphcNOvY0SbuarRpOP2vDOmbTMhC5Rrk32/RkNmrT6z/M2u33LXIDILQIJyEntgreNnhdQNpyqyjRG8hp/7CPYYi1lzz9hVDJOKm1H8dK06aH/ear523+S4cEEXqqlBJofmueBsvO7xpOL3hsvTfHMhC5YpyD182bP7XF//77eWtctciJwjC9ulF9Lfdgh2D3x8kWwoiSZRO7cSxCTJd3mY1CA1rkt8v05XOiX3ZiXUy67Uxxt0GkEtHzHrR2/ja44+H0mvOib/BJhUKFujsv3rk0oP5R947sVGxUwwhCNvBS+jv+0UK4/WDaVquFESIpJ1ELt7IJ0S2CmxRub582el33VSuL3Z/Xk3L0NTG3YeQS4eRJFsOhTrhL4Yzy5KlnyELFcxd47oyfmFZQ/nsvYur9OYewGUJIAjbwkvo6QMiL6JPh9Hm2VaiNeTkdtz7HjkrsDlZVTlLk1YM9hnwTOTfKSzTq+vqjRzdUP4lea6OEEIoxAl/MpReeVaCdqGSqRn1G4Nf7ukZ/dKueYXVV+Qux9wgCFvFSeipfWIdjz4fTrPyfp8qb5CidBwN2/AaTdKNU/88/+HEqPGj/OPlrQRHD5JryEyTcGf86TD6bchCZcMIj+8+bmzEg6/uXXSxJFXucszKPPsGWR+9iJ7aL+oFtEn2FESInNyJe8bDMBmjIIj8nPPrgauHX42ZoXXwl7schKP6Scd/x/o6pDbOot6dE+6MvxhOv3RUFBB6IkTun3ggnwdCR3nbey09tOqlPlMfCrtX7nLMBH7iW1AvoEd2CQxGX8bTKpnmStwiidKZXaj3SLnrsAV6Ub/+wsfnyi4siHvNElIQIYTUDjg0hpw7IHcdt/pIt+ZCu1DRevvErBqx5LvLP/0z6WOFDJ+BIGyuTkCP7hKcVfijwTLfF2xEUo5jd1/sGSB3IVavXF++/PS7COHZMTOdVRa0myPucy85swtZwDtOuDPeFE+vOS9tgrn2yubn6PNuwtLrdSULDi5XwqqkFvBOb0mqOPTATsHfAW8YYhEpiBAiSdsRDJPpsvTKzKVJK3p79ZrabZLRd5PoKs9A7OCO0k/LXQdCCAU74i+H01sypY0pkIWKZs/aLRj8SpCT9h9/vZ5XJc9sV7OxjDd7y3C9AY38U+jhitfKN2u+GVJeTK5l46gBchdi3XYX7fvgwsYp3Z55QGup9zx63yOd3C53ETf5O+DN8czWHOndc/I3UoGMKExNinnqiahHZ+9ddKQoSe5yTAiC8Ka8GjL8T2G4L17RT8bpgs2RkztwrxGIsbAWjPXgJf6LlM17i/a/ETc72jVK7nJahcNicV0VLs6Ru5CbPDVoywjm0DWy4JQowaZNyjYqeNjSoXM3nN70RfLXErHNfgIIQoQQulBOhv0pvhBFLe4t+9iY20gCOrsHx0G/aCddr7+x9OTyer5hXuxsT42H3OW0CWMUO5I+u0vuOm5xYtEnw5iUCvLKcZG3zXc/YKgw1+D3Rr2dWpbx2r6lFbpKucsxPghClHSDPPiX8FZvakqkZX03mAsHkU8ocveVuxCrlFx6/s1TK/t59ZkaPUklx6oxHYV7DsUFqaimTO5CbnFg0GdDmXoBTdgnVPNyVwNk5aRyXDJ0TrR75D/+ej2lNEPucozMst76ze+3fOnRPcJHg5mxwRb2rZAk5vivaMhYueuwPiIR/5v905epX83o+fwDAZZ6U/BurJpE9kdn98hdxx1YGq3uTwc64cf2CNcboJNU0ShMje8+bmqvCYsSV/yWuVPucozJwt79zWv9RWnGUembBCbe13JuC95Ezu8nju44sJvchViZcn35yjNrsyqzF/aZG+oUInc5HSPGxKMLiYjTyV3IHSiMFsbR92upR/eImdWQhUo3yL/fO/GLf83Y8faRtXV8vdzlGIdCg1Ak6JXj4qZ06bf76Th3i0tBRAhJ/FEY8JDcdViZS+UpS0+ujHAOn9nzH06so9zldBhxckf+keTiYbkLacHUKOqlaOqJvcLJG5CFSvf/7d1nWBTX/gfw35mZ7bBLWfqCNBUQULFgR1FsKWqssaFR03PjNYkaExITvTExMTfFNE27GFsS/zawoaKiGAUUewFBLCh1YXfZOjPn/4LE600sGGFnlz2fJy92cR78Tnad35wzpwS6+7+f9KaCUUzPfOls9QWh4zQDV1xizWCDCdmswQYbUxilQ47HxKdyQCzhNVFCB3EaHOY2lm7dV35gRtQUu+2s2yI69sfZa1DCIEAOd3/2SDCllqCZOWxaAj02zEXvoYlGDC2a0XFirE/U6/sXjYsaObHDaAQO941tOpf7Nt8wQlImqxSh9CQHrYKAMb9vLeo1SugcTuOmseLtvMVFdcWvd3rVuasgAASEI4kMLhYInePOEn3R932ZT07zaQVkWgUBiYEJS5LSsq8cfH3fYq25Xug4f59rFcIjlbj7ZnaohlqWKPxS2neDz/0GgCG8o9BBnENO+aG0vEXx6thnYmYoRHKh4zQD1H04v2+NI6y4dkfB7uiHJOZULZ6ewzWQoaQuz0/h817SGyFKzfRtLx26dlToOH+To1aDFrC6mH8si323C/VijEOfNc5eg3o94YA9Y45GZ9N9dOLTTWWZs+OeHxw0kGo1/8faxCKFCk474pPCRioRfNmL8RbDY1nsVQNpGLo6GtHjo0fMTXzxs/wVi3P/bWYda7RXUzh0SWguNh5ezOUWHuM3DGSGahz6lHFRAVhNqF0XoYM4usKaE/MPv6WWeL/e+RVH2Uei+aBeI/j9P4PNKnSQu6IpmNeRfrwN9VgWe7CC1EIC2ntFfjTwHRtnmbl99oWaYqHjPJjWP1im2gzj9rA0BRlDHPWh4G3w3jXQayRpDt6D0WZcVbT2VO2ZaVGT26kihY7TMnxDkV8oLtiJejwmdJR7eTKcaq9EL+VykyLRnFhHWaGXEIqckf2j69MHrx2Zm73wsbZDU2PHi2iHv+YCQKtvEeZW4IRNbEdv9GM/Z6iCJYVYX4uieggdxHGdqD4198ibHM+9mTCv1VZBAABAPR/HRzLA0iB0kPtIUKNVSXR2OZ6ew5HVZwgA6KNJ/PegxcXa0pnbZp+rcY41aFptIcQAn5ziR2SxaZ2oufHOcK/Kc3jrV1T/cUC12g/lYZhY08qzP3x74cfJbSdMbDtOSkuETtTCPP1RWDw+nCF0jvvzkaGVfRiNHAZuYwtrSDcpASqJ8rXEF8ZEP/76/sWf5a80sxahE91H67zmVpth+A52TQm/YyjzaIhznCM+vAXLlNA+UeggjuhIZd4rhxfwCKd1nuvIm0g0L5T4CC7cAzoHWn30bmgK/tGB/kcMNXkfm15ElugmAAB6BXX7OHnRDUPFrB3/PFF5Rug49+IcReKBHKnE3TazGgX6dSATKHf8liAAADZo+X1rqcGpQgdxOFXm6qWFH/9S/H8zoqaMD39C3OobgrdTeKCYXvjQRqFzNNVgDZWexKy7xE/bz9U6ehuAsAeVxP3VxOcnRI965+BH7//2qc6iFzrRnbWqQmjj4a0CbsQudnEX+t0ujjtT8K/w9u9QXBJ4t7bRjw+Dw9z2K1lvHHknUB40v/MrTj9T/m9BXQbji/lQfV3oIE0VpEAr+zIBMhi0nYwmJX7XI7DL8pQlMlo2eetzW4t3YnC4L0brGTV6tg5P2ceppbBzuMhXKnSaB1J2Bl86Ts36UOgcDuR83cXvLqR7iFVzO8129K0EW5REjrqk4N3paPx8ZxlLzFAwO47u6YdfyuVGh6F58c50S0q0EAkjnhw7pqem6zfH/7O7dP9LCbOC3Rzovr81fEN5DP8+xSdlsE+GUz/2Y5ysCvI8t+ULNOBJEDtX7pZSY6n9+sK3y898PSw45YUOT7t0FQQAABTfH0x6OLlP6CAPJtEXrU1mzmph6A72tNbhWgCEICI8Qt/vn9YloNMr2WlfFP6gtxqETvQ7py+El/V40DZ2fSm/dTAzKdL5Toc/mokYCYruKXQQ4Vl525ayzAVH3nZn3NM6z09QdxI6kWOgaBg0hd+3HnTVQkd5MJ5i+HcPelIk9eRedulJss09AQBAIeqRiEHLBrxr5ayTtjz76/mtPBb+m+F8leMWDsOyk3zXzWwff+rXgUwbN+foOPofRh3sXY0GpzpLr1cLwYAPVxx5JXf+xbrieZ3++UjwUIkz7ClvN8gzAHXsz2d+A9j5mlaPBFPrkkUnamDwdvZErfPlJ1qCm0iR2mHcW31ezb5ycOa22cdunhQ2j7M+Izxeg2ce4BQi2JrChLo7axXht3wBMb3BN0ToIEI6U3tuTfHPHOamtHuycY68yex8axW2NJSQgn8phFP7Ib6/0FkemFoKH/egt1/jp2Sz48KpOXG03FkvPERzClOFvNt3Xu71vKVHPte4Bz7beVqkZ5ggSZzv+2hk4e0CLr2YX9CJHhvmBBPl7wbn74AbJSh1kdBBBHO14fqvlzYW60qGBaf08uvRelbNbgkUjVJS+Y2fUqFxoHTK56bDNFSiD/XpaS4pk303gRoW7MTdUUQz6hXUrXtAQvaVnLl7F8b6Rj/beVqgm7+dMzhZIVxfwr92lE/0QXuGi7ydeUYZrizjd36PJqaByBX7AG8aKzaUbDpVe2ZocMrUdhNpRK6J94e8AiA+ic/8mpqwwEn70r0k8E4X+lg1fv8E91Mx/6+utPN25xDNiKHolND+fTU9NxfteHrHnCFhAybFjPGSedotgNNcgI7X4H4Z7L+O85/1oD/rSTt1FQSrmV+9mBowCamDhI5ibzeNFV+eWfl2/mIvidc7Xd/sH9iXVMGmQ11SwKiDUweEDvJQEtRozQAmQU09uov96CRnZIUORDgGKSMZHz3is0HvWVnr1MwXlhd8W2vS2uevdoJrULUFPX2QG7aDfTyE2j6USfR1+ltIfvNyFBAOcX2FDmJXVebqlWd/eDt/sVKkXNhlwbCQFDIi5oFRDDUolc9eCxVlQkd5KAwFkyOpdQOYCzroncGuK0Vkv3uikUqinBb/5Ocp73E8NznjuQ+PLK821bb0X+oEhfC1YwwFsP9R0aRIyokfCf4BF+zCZWdQigutplaiv/zpqS/fPPqOm8htYdc3Hm0zRMbIhA7ltLwDUdJ4/peloLfTzXLL8ZGhxV3ojxPpX8vQ4J04+wYphsTvlGLl5Nixnw9aQgGalvHisqNfXdffaLm/DmG7D8hes2ZNZmbm6tWrm3h8VZ2eFbu3aCS7wZVl/IpXqQkLwK/NfQ82NjTIFQo7pGo5F+oubr6cUaK73Me/10BNfxndpEUDTGazTCJx0sdgD8NssYhFDEXR9z0S5++A0lNo8tut4xmz1Wo9VIk+P480Cpjfke7s7SofvdVqRQiJRA6/RVxzs9lsGGOxuEnfXp1Vt/1S9o6SPTHe7abGjY9Rt2/2PE42WMa5NdTzqxehARObUgWdGsvbciuOZpbt5DCbohkwI2oaeRDYvFCXIVh7E2/7Gj3+Uuu4Y+jnj/oFMpuu8LNyuCgVvBpPd3KZckjcm1KsHB89YmTboVmXD7yds1SjDHgv6U0Z05xLcZFCaCfY0oB/eAO1T0TxSUJnaUE1ltqsq3uyy3OC3YJGtBke7dUeAbmctQCEUPIkfuOncHAD6jtG6DTNg6ZgdCg1IoTacY1/9hAXooDXO7lQ65C4NwkjeTQyZVhE8rmaInFzb3xPCqFdWC34x7cgKBL1HS10lBaBAZ+tPb/r2u4z2vOJvl1fiX/JV+YjdKjWjhZRjz6Lf/kQvAMgprfQaZoNQ8GjIdRQDbWxjJ+Zw0W6w/Md6H7+raLZSzw0GtFd/DrS6P6PDx4IKYQtz2bl09PAyx8NnCJ0lOanNWv33cjZV54joSV9/HtOjBxPxoLaj9QNDX+G3/QZJZJC2y5Cp2lODAVjw6hRbajsG/y/Crl3OHiqPTU2jJI08wWQIABIIWxpmGPx2sUglqPB01vHs5xGNt52vPpkdvn+4rpLCT6dnoqa2sYtWOhQLsk7ED3+PN76FTLqoWN/odM0M4aClCAqJYg6XIFXFXPLTvPT21ITIyg12amFaFakELYkjuXXvY8wRiNfAKo1jBbhMT6nPX/wZm5+5TGNe1AP327T208RUy435s2hIJ8QeGI2v2U5ZWmA7o8IHadF9PRDPf2YIh1ef4nvm2FLCqBSI6mefq3nzpIQFimELQXrtfzqd5FciUa8CE0YEO/IeIwv6S79Vpl/+MZvSol7V5+ENxLmekhUQuci/uDhR42ew29ajkwGlDRe6DQtpa0SvdmZ/mccnXGFn5/HYYApkdSoUNJAJB4WKYQtAl+7wP+0CHUagHqNdN4eUQ5z57QXjlTmF1Qek4vkHb3j/xH3rL/c3uvhEk2i8KRG/xNv/RKbDGjw9NbRA3FHCgbGh1Pjw6nj1XhTGb/stK2bNxoTRg3RUFJyPSP+FvLFaX78sSy8/Vtq+EyIdMrxCwbWcLLmzPGqwsKaU2qpd2d1/Mtxz/vJfYXORdyP1A2NfBlvX4lXL0KPPw+qVj5wt7MadVbTZpbee4NffYmfl88N01CPhVB9/ZGo1d4GEC2CFMLmhDkWb1+Jzx+lJqWBd6DQcR4Aj/EVw5UTNaeOVZ242nCtnSqyg2f08JAhnhIPoaMRD0IkQY+9gI/v5n98kxo4BWL7CB2oxUkZGB5MDQ+mqs2w8xr/8SnuhVw8MIB6tA3q70/aiESTkK9Js8Glp/Dmz8HDl5q2CCRyoeM0yU1jxenaM6dqz57VnncXu0Wp2g8JTm6rihSR8S/OCyGUkIJConHWf6C4AA2dCVLnXqividRSmBRJTYqkqs2wt5z/5hz/8mGutx81MBAlB6IAubM+oSDsgBTC5tBQx2//Dl8qRAMno/bdhE5zLzzG1xuun6u7cL7u4nntRQCI8mgb49l+VNhjHmIy+KUVUWvQuLk4dzP/7TxqyPRWNsvw3tRSGBdOjQun6m2QW8Hvu4nfO8H5S1FyEEoKoLp6I9JMJP6EfCMeDsb80Uy8exWK7UvNWgoiR9wm0WAzXKovvaQrKaq/dLG+SClWhruHtlVFDNOk+MjUQqcjWgwtQn3HoNA4nL0OH95CJU2ANtFCZ7IrlQiGaahhGsCYPq3Fhyr49wu58/U43gv19qV6+6MEbyR27gHdRPMghfDvYm24cA+fswHJ3NGTC5CPA00nr7fqr+ivXDaUlerKLulK9FZ9qLJNG7fgRN9uE9uOcxe5CR2QsKPg9ujJBXAxD2/7Gjz9Uf8J4B8mdCZ7QwjivFCcF/1sNJhYKKzBedX8wmN8UT2O9kBd1airD+qqRv6k+9RVkUL44Iw6/rcM+G0L9g+jUqZCmw7CxrHytusN5dcM168arl01XCvTX7HyVo1bkEYR1E4VOSQ42UfqSzntFA6iGSCE2neHtglwJpf/5UMUEIE6D4Sw+FY8xeIeZEzj9HwaAEwsnNXiE1q8tgS/nseJaBTnieK9ULwXxHqSx4ouhBTCJuN5KD3Jn9yPTx9A7brBkwsotcbOETDgGnPtTWPFTePNG6bK8oby8oYbdZY6P7m/v9w3UObfw6/72IiRXhIvOwcjnADFQFw/KqYnvnCUz/kVtq+k4pMgPgk8/IROJhgZA118UBefxoJHX2vA5+vw+Xq88jw+V4cBQZQSRXmiKBVqr0LtPJCSjCFrpR6sENpstvvuIdmUY5wI5lgoOYFPH8RnDyGVD7TvjmZ+hNxafFxJvVV3veF6Q4OxylRVYa6qNFVVmaoqTdUKRu4r9/WRqn1lPok+3QJD/dVSb4rs9kc0ES1CMb1RTG9ce4M/mwv/eRt5B0G7BBTeEex+Y+doNAqkUaBBQb+/rTLhS3oorse5lTi9mC/RYTkDEUoU5o4ilCjMDcKUKFiOFK3naue6mloIDQZDampqVlYWTdPz5s2bP3/+X49ZunTpkiVLOI4bOHBgenq6u7vTbivP2uD6RVx2hi89DVfOgHcgatcdpS5CzTpD2cia6i31OptOa6nTWuq1llqtpb7GUlNrrqu11EppqYdY6S318pZ6e0u9IpVh3lK1WuIlJns7EM0BeQWgPqOh1wi4chaXncX5WZi3ofCOKKwjBEWC0lvogMLzkSEfGfTw/W8HaYUJlxngagO+oscHb8JVA3/diGUMBMtRsBsKVkCgHAXIwV+OguTgI0M06Vt1Ek0thO+9957BYKiuri4vL09MTOzXr1+vXr1uP+DIkSNLly7Nz88PDg4eNWrU4sWLP/jggxYI3DJMelxRBpVXcMVlKL+Eb5Qg7wDQtEcxvWDIdKR4gPYfBtzANjTYjA2s0WgzGmwGg82gb/zPajDYDDqbrs5Sr7cZKESpxEqlWOkhVirFSpVYFakK6yru7CFWeUo9xZTIZDTK5M4xH5FwVhQDofEoNB4AUF0lvnKOP7EP7foBAyC/UBQQDv7h4KMBlRpo8hgF/GTITwbdff6nvtVaoNyIyxvwDROcr8cHbkKlmb9pBK0Ve0rAW4z85OAjRT5SUEuRpwS8JMhLAl5i8JIimVBnQvwvhDFuynFBQUE//vhjSkoKAMyePdtoNK5YseL2A5577jmappcvXw4A2dnZEyZMqKiouOOvWrNmTWZm5urVq5sYsapOz4qbo3FpMWGjDhrqQK8FbQXWVUNdJa6rQrXlmLWBjwapNcgrwOoTZPMLNtPA8qyZNXEYG1mjjbdZeauRNbG8zcSZTTajibNYeWuDzWhmzUbWZOZMRtZkYs0mzqRg5DJGphDJ5YxcTsvkIoWbSKFg5ApGoRDJ3Rg3lUTpJnK776YNLlsITWazTCJx3jVa/zazxSIWMZQjLNHeUIcrr0JlGVRdA+0NbKhDChX29EWeAUjlA24eoFCCmxfIVaBQPvwnZbVaaYqimVZVa3kMtRbQWnGlCWotuMYMtVZcb4E6K66zQp0V6ixYZwOlCDzESCUGlRiUEuTOgJsIlCJQiJBSBAoGZAwoGOQuAgkNcgbkDBJToHLyjiGbzYYxFov/zmlIaImbuJnXiGjSN89oNJaXl8fExDS+jYmJWbdu3Z+OKS4uHjVq1K0DKisrdTqdUqm84y+0Wq1arfbWWw8PD3T3f0tF5ccr9JWNrzmrycxbAADzPLBWADBxZp61AsbAWVmeN9uMwPOYs7GcxcJZgeOAZ828leNsgJBJzPAUjWmRSQRAMTYFsrphCJEOOrzHAAAKNElEQVQZOYShwsyVcTpeZBCJy0RiWsIgWspIKUTJGZmIYkRIJGVkDMVIKImMlihFSgnNyBi5lJJKGamUlkgZqYyWyhhyk0c4P4UHCvOAsLjGd4jnQF8L9VVQX40NWlxZhkx6bKhDRh026kEsQ1I5liiQTAESGYikIJaCWIpoEUikQIuAEQMASN0AABjmv9NtmT/+yGoFmgaGAcn9/vlQDIidY7MJCkANoBZD27tf7a1Wm4FFJkzrbdBgAz2LjSw02KDBDFV6fJkFIwcWDswc6K1gxWBhwcRiG4DeChQChQgoADcGURS40YAApAyIaKARKBhAACIaZDQAgJgCmQgaL7IK+r99tgoRMLcNMPjTaCAKgRtzhyuzlAHx/YYlyGgQ3f2Ojuc4nucZkfU+vwVAziDR/0aQing3byEKYWPRcnP7ff6ZUqmsra396zG3Dmh8OqjVau9YCIuKirZs2bJ79+5bP1m7dm2fPnddFPF4wa+XLRWAAQAompZhGgAAIUzRCEAKNE2LMEJAMQyiZCIpMBTQ7iJGKhZJgREhWiSVqhiJHGhaQksYxACAjJEihBjENG6nLhfJESApJaEd4Wb8D5QNK8SuWFYZDskkznGxa14MhyRiCUU5ZFNY6gY+IXf8E2xqAIsRLCYwG7DFBKwVrBZsNQNrA70ec1ZgWcAYzA2AELA2bPvj8sexiLUCgARjhBAGjMzGe6fAHIds5mY9MSFJANwB/t7njQEwAODGSyPc6tlr7OPDtx12+89v/8lf32Lc1CxN6Um8xyFNv86aAf70eeuQCOZ8pnJr6rVRLpdT95sp1KRCqFarEUL19fUqlQoAtFqtn9+fh1z7+PjU19c3vq6rq2v8yR1/W9u2bceMGdP0rtGpI//lxONuHoJer3fNE29oaJDL5ffoJGitjEajRCKhaQe6G7MPs9nMMAzTurpGm8JisSCE/l4PoVOz2Ww8z0skjrIUV5OG3UskkvDw8OPHjze+LSwsjIqK+tMx0dHRtx8QGhoqd8nnWwRBEIRzaer8s1mzZi1evPjq1asHDhz4+eefZ86cCQA1NTVDhgypqqoCgBkzZmzYsCE7O/vatWuLFi2aNWtWC6YmCIIgiGbS1L6IV155pbq6um/fvh4eHl999VVcXBwAYIwtFkvjuNMOHTqsWLFizpw5Wq129OjRr732WgumJgiCIIhm0tTpE83oQadPLFu27KmnnvL09GzRVA5o+fLlY8aM8ff3FzqIvX333XfJyclhYS63NvTatWtjY2Mb7zJdyubNmwMCArp37y50EHvLyspiGGbAgAFCB7G3Q4cO6XS6YcOGCR3kd06wNFd6enpZWZnQKQSwfv36CxcuCJ1CABs3bjx58qTQKQSQkZFRUFAgdAoB7Nq16/Dhw0KnEEB2dvaBAweETiGAnJyc2ycOCM4JCiFBEARBtBxSCAmCIAiXRgohQRAE4dIEGCyzZMmSJUuW3G26/V9dv37d19e3NW3t1EQ3btzw8vJynDmndlNRUaFUKmUyl1tVp6qqSi6XKxTNvHyU46upqRGLxS64fIRWq0UIeXh4CB3E3urr6zmO8/Kyx86pEydOXLRo0b2PEaAQ8jxfVFTU9MJmsVhcsBiAC5+41WoViUQuuLKMzWajafq+y0G1PizLIoRccEkdjuMAwDVPHGNsn7WEAgIC7ntXLUAhJAiCIAjH4XL3ngRBEARxO1IICYIgCJdGCiFBEATh0kghJAiCIFyao28AtmXLln379gUFBc2aNetu+923PtXV1QUFBWVlZX369ImJiRE6jv2UlJRkZGSUlZUFBQVNnTpVrVYLnchOjh49umfPnqqqKl9f34kTJ4aE3HkL3FZs06ZNVqt13LhxQgexk1OnTt2+qtyECRNc5/pmtVpXrVp1+vRptVo9fvz4yMhIoRM5dovw008/nT17dkRERG5ubnJyMs/zQieyk2HDhr311ltpaWn79+8XOotdjRw58vTp0yEhIXl5ebGxseXl5UInspNNmzYZDIbw8PCSkpK4uDhXW2P2wIED06ZNW7hwodBB7CcrK2v58uUlf7DZbEInshOz2ZyUlLR69eqQkBCj0Zifny90IgAAwI7KZrMFBgbu3r0bY8yybGho6LZt24QOZSeNk2wGDBjw5ZdfCp3Frkwm063X3bp1+/zzzwUMI5RBgwZ99NFHQqewn4aGhvj4+A8++CA6OlroLPazbNmyadOmCZ1CAEuWLOnTp0/jJc5xOG6LsKioqLq6OikpCQBomk5OTnad5pELTqluJJVKb722WCxubm4ChhFERUVFUVFRbGys0EHsJy0tbcqUKREREUIHsbfi4uKlS5emp6fX19cLncV+tm3bNnny5I0bN37yySfHjh0TOs7vHPeCe/PmTS8vr1tLD/j5+blORxmxYsUKg8EwduxYoYPYzxdffBESEhIcHJyamjpkyBCh49jJkSNH9u3b9/LLLwsdxN68vb3btm2r0+nS09OjoqIuX74sdCI7KS0t/fjjj3ft2qXVaocMGbJq1SqhEwE48mAZhmFYlr311mazueZ6Yy5o586daWlp27dvd6lVN2fMmDF27NgTJ05Mnz69Y8eOTzzxhNCJWpzFYnnmmWe+//57F1xJODU1NTU1tfH1uHHjlixZ8s033wgbyT4oikpMTGw82fbt2zf2BwgdyoFbhIGBgVqt1mg0Nr69fv16QECAsJEIO8jKypo6derWrVsTEhKEzmJXUqnU19c3JSXlqaeeWr9+vdBx7OHYsWPFxcVPP/10165d586dW1pa2rVrV5fqJ2zUu3fvkpISoVPYiUajuTUSvkOHDlevXsUOsMyn47YIIyIioqKiNm7cOGnSpPr6+qysrDlz5ggdimhZOTk5kydP/uWXX7p37y50FrsyGo1yuRwAMMb5+fkdOnQQOpE9xMfH33rwv3fv3q+++uqbb75xkW4Ak8nUuBI0x3GZmZmu81R45MiRubm5ja8PHz4cHR3tCMvrO/Si2xkZGdOnT3/kkUfy8vI6d+78008/CZ3IThYsWLBr166LFy96enr6+Ph8+OGHAwYMEDqUPfj5+WGMb82imzx58uzZs4WNZB++vr6JiYmenp4FBQU0Te/Zs6fp+5S1Dhs2bEhLSzt79qzQQeykR48eKpXK19c3Ly9PIpHs2bPHRWbN1tfX9+vXz9/fPzAwMDMzc926dcnJyUKHcuxCCABlZWW5ubkajaZPnz6OcONgH6WlpbW1tbfeRkREuMiOZYWFhY0b0zTy8/PTaDQC5rGbysrKvLw8vV7fpk2bxMREFxw2rNVqy8vLXaQpDH984jqdrk2bNj169HCpT9xisezdu9dsNvfu3dvX11foOACOXwgJgiAIokW50G0IQRAEQfwVKYQEQRCESyOFkCAIgnBppBASBEEQLo0UQoIgCMKlkUJIEARBuDRSCAmCIAiXRgohQRAE4dJIISQIgiBcGimEBEEQhEsjhZAgCIJwaf8PgCbn+taLk5kAAAAASUVORK5CYII=", "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)$. \n", "\n", "- For example, [the product of two independent variables that are both Gaussian-distributed does not lead to a Gaussian 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. (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 identification 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.10.2", "language": "julia", "name": "julia-1.10" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.10.5" }, "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 }