Chapter 13 - Statistical Decision Theory.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 13. Statistical Decision Theory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 13.1 Preliminaries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Decision theory** is a formal theory for comparing between statistical procedures.\n",
    "\n",
    "In the language of decision theory, a estimator is sometimes called a **decision rule** and the possible values of the decision rule are called **actions**.\n",
    "\n",
    "We shall measure the discrepancy between $\\theta$ and $\\hat{\\theta}$ using a **loss function** $L(\\theta, \\hat{\\theta})$.  Formally, $L$ maps $\\Theta \\times \\Theta$ into $\\mathbb{R}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **risk** of an estimator $\\hat{\\theta}$ is\n",
    "\n",
    "$$ R(\\theta, \\hat{\\theta}) = \\mathbb{E}_\\theta \\left( L(\\theta, \\hat{\\theta}) \\right)\n",
    "= \\int L(\\theta, \\hat{\\theta}(x)) f(x; \\theta) dx$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the loss function is squared error, then the risk is just the mean squared error:\n",
    "\n",
    "$$R(\\theta, \\hat{\\theta}) = \\mathbb{E}_\\theta(\\hat{\\theta} - \\theta)^2 = \\text{MSE} = \\mathbb{V}_\\theta(\\hat{\\theta}) + \\text{bias}_\\theta^2(\\hat{\\theta})$$\n",
    "\n",
    "In the rest of chapter, if the risk function is not specified, assume the loss function is the squared error."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 13.2 Comparing Risk Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **maximum risk** is\n",
    "\n",
    "$$ \\overline{R}(\\hat{\\theta}) = \\sup_\\theta R(\\theta, \\hat{\\theta})$$\n",
    "\n",
    "and the **Bayes risk** is\n",
    "\n",
    "$$ r(\\pi, \\hat{\\theta}) = \\int R(\\theta, \\hat{\\theta}) \\pi(\\theta) d\\theta$$\n",
    "\n",
    "where $\\pi(\\theta)$ is a prior for $\\theta$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An estimator that minimizes the maximum risk is called a **minimax rule**. Formally, $\\hat{\\theta}$ is minimax if\n",
    "\n",
    "$$R(\\theta, \\hat{\\theta}) = \\inf_{\\overline{\\theta}} \\sup_\\theta R(\\theta, \\hat{\\theta})$$\n",
    "\n",
    "where the infimum is over all estimators $\\overline{\\theta}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A decision rule that minimizes the Bayes risk is called a **Bayes rule**. Formally, $\\hat{\\theta}$ is a Bayes rule for prior $\\pi$ if\n",
    "\n",
    "$$R(\\theta, \\hat{\\theta}) = \\inf_{\\overline{\\theta}} r(\\pi, \\overline{\\theta})$$\n",
    "\n",
    "where the infimum is over all estimators $\\overline{\\theta}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 13.3 Bayes Estimators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $\\pi$ be a prior.  From Bayes' theorem, the posterior density is\n",
    "\n",
    "$$f(\\theta | x) = \\frac{f(x | \\theta) \\pi(\\theta)}{m(x)} = \\frac{f(x | \\theta) \\pi(\\theta)}{\\int f(x | \\theta) \\pi(\\theta) d\\theta} $$\n",
    "\n",
    "where $m(x) = \\int f(x, \\theta) d\\theta = \\int f(x | \\theta) \\pi(\\theta) d\\theta$ is the **marginal distribution** of $X$.  Define the **posterior risk** of an estimator $\\hat{\\theta}(x)$ by\n",
    "\n",
    "$$r(\\hat{\\theta} | x) = \\int L(\\theta, \\hat{\\theta}(x)) f(\\theta | x) d\\theta$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 13.8**.  The Bayes risk $r(\\pi, \\hat{\\theta})$ satisfies\n",
    "\n",
    "$$r(\\pi, \\hat{\\theta}) = \\int r(\\hat{\\theta} | x) m(x) dx$$\n",
    "\n",
    "Let $\\hat{\\theta}(x)$ be the value of $\\theta$ that minimizes $r(\\hat{\\theta} | x)$.  Then $\\hat{\\theta}$ is the Bayes estimator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proof**.  We can rewrite the Bayes risk as:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "r(\\pi, \\hat{\\theta} &= \\int R(\\theta, \\hat{\\theta}) \\pi(\\theta) d\\theta = \\int \\left( \\int L(\\theta, \\hat{\\theta}(x)) f(x | \\theta) dx \\right) \\pi(\\theta) d\\theta \\\\\n",
    "&= \\int \\int L(\\theta, \\hat{\\theta}(x)) f(x, \\theta) dx d\\theta = \\int \\int L(\\theta, \\hat{\\theta}(x)) f(\\theta | x) m(x) dx d\\theta \\\\\n",
    "&= \\int \\left(\\int L(\\theta, \\hat{\\theta}(x) d\\theta \\right) m(x) dx = \\int r(\\hat{\\theta} | x) m(x) dx\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "If $\\hat{\\theta} = \\text{argmin}_\\theta r(\\hat{\\theta} | x)$ then we will minimize the integrand at every $x$ and thus minimize the integral $\\int r(\\hat{\\theta} | x)m(x) dx$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 13.9**.  If $L(\\theta, \\hat{\\theta}) = (\\theta - \\hat{\\theta})^2$ then the Bayes estimator is\n",
    "\n",
    "$$\\hat{\\theta}(x) = \\int \\theta f(\\theta | x) d\\theta = \\mathbb{E}(\\theta | X = x)$$\n",
    "\n",
    "If $L(\\theta, \\hat{\\theta}) = |\\theta - \\hat{\\theta}|$ then the Bayes estimator is the median of the posterior $f(\\theta | x)$.  If $L(\\theta, \\hat{\\theta})$ is zero-one loss, then the Bayes estimator is the mode of the posterior $f(\\theta | x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proof**.  We will prove the theorem for the squared error loss.  The Bayes rule $\\hat{\\theta}$ minimizes $r(\\theta | x) = \\int (\\theta - \\hat{\\theta}(x))^2 f(\\theta | x) d\\theta$. Taking the derivative of $r(\\hat{\\theta} | x)$ with respect to $\\hat{\\theta}(x)$ and setting it to 0 yields the equation $2 \\int (\\theta - \\hat{\\theta}(x)) f(\\theta | x) d\\theta = 0$.  Solving for $\\hat{\\theta}(x)$ we get the given estimator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 13.4 Minimax Rules"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem of Minimax Rules is complicated and a complete coverage of that theory will not be attempted here, but a few key results will be mentioned.  Main takeaway message from this section:  Bayes estimators with a constant risk function are minimax."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 13.11**.  Let $\\hat{\\theta}^\\pi$ be the Bayes rule for some prior $\\pi$:\n",
    "\n",
    "$$r(\\pi, \\hat{\\theta}^\\pi) = \\inf_{\\hat{\\theta}} r(\\pi, \\hat{\\theta})$$\n",
    "\n",
    "Suppose that\n",
    "\n",
    "$$R(\\theta, \\hat{\\theta}^\\pi) \\leq r(\\pi, \\hat{\\theta}^\\pi) \\;\\text{for all } \\theta$$\n",
    "\n",
    "Then $\\hat{\\theta}^\\pi$ is minimax and $\\pi$ is called a **least favorable prior**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proof**.  Suppose that $\\hat{\\theta}^\\pi$ is not minimax. Then there is another rule $\\hat{\\theta}_0$ such that $\\sup_\\theta R(\\theta, \\hat{\\theta}_0) \\leq \\sup_\\theta R(\\theta, \\hat{\\theta}^\\pi)$.  Since the average of a function is always less than or equal to its maximum, we have that $r(\\theta, \\hat{\\theta}_0) \\leq \\sup_\\theta R(\\theta, \\hat{\\theta}_0)$.  Hence,\n",
    "\n",
    "$$r(\\theta, \\hat{\\theta}_0) \\leq \\sup_\\theta R(\\theta, \\hat{\\theta}_0) \\leq \\sup_\\theta R(\\theta, \\hat{\\theta}^\\pi) \\leq r(\\pi, \\hat{\\theta}^\\pi)$$\n",
    "\n",
    "which contradicts $r(\\pi, \\hat{\\theta}^pi) = \\inf_{\\hat{\\theta}} r(\\pi, \\hat{\\theta})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 13.12**.  Suppose that $\\hat{\\theta}$ is the Bayes rule estimator with respect to some prior $\\pi$.  Suppose further that $\\hat{\\theta}$ has constant risk: $R(\\theta, \\hat{\\theta}) = c$ for some $c$.  Then $\\hat{\\theta}$ is minimax."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proof**.  The Bayes risk is $r(\\pi, \\hat{\\theta}) = \\int R(\\theta, \\hat{\\theta}) \\pi(\\theta) d\\theta = c$ and hence $R(\\theta, \\hat{\\theta}) \\leq r(\\pi, \\hat{\\theta})$ for all $\\theta$.  Now apply Theorem 13.11."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 13.15**.  Let $X_1, \\dots, X_n \\sim N(\\theta, 1)$ and let $\\hat{\\theta} = \\overline{X}$. Then $\\hat{\\theta}$ is minimax with respect to any well-behaved loss function.  It is the only estimator with this property.\n",
    "\n",
    "*Well-behaved means that the level sets must be convex and symmetric about the origin.  The result holds up to sets of measure 0.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 13.5  Maximum Likelihood, Minimax and Bayes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For parametric models that satisfy weak regularity conditions, the MLE is approximately minimax.  Consider squared error loss which is squared bias plus variance.  In parametric models with large samples, it can be shown that the variance term dominates the bias so the risk of the MLE $\\hat{\\theta}$ roughly equals the variance:\n",
    "\n",
    "$$R(\\theta, \\hat{\\theta}) = \\mathbb{V}_\\theta(\\hat{\\theta}) + \\text{bias}^2 \\approx \\mathbb{V}_\\theta(\\hat{\\theta})$$\n",
    "\n",
    "*Typically, the squared bias is of order $O(n^{-2})$ while the variance is of order $O(n^{-1})$.*\n",
    "\n",
    "As seen on the chapter on parametric models, the variance is approximately:\n",
    "\n",
    "$$\\mathbb{V}(\\hat{\\theta}) = \\frac{1}{nI(\\theta)}$$\n",
    "\n",
    "where $I(\\theta)$ is the Fisher information.  Hence,\n",
    "\n",
    "$$ n R(\\theta, \\hat{\\theta}) \\approx \\frac{1}{I(\\theta)}$$\n",
    "\n",
    "For any other estimator $\\theta'$, it can be shown that, for large $n$,  $R(\\theta, \\theta') \\geq R(\\theta, \\hat{\\theta})$.  More precisely,\n",
    "\n",
    "$$ \\lim_{\\epsilon \\rightarrow 0} \\limsup_{n \\rightarrow \\infty} \\sup_{|\\theta - \\theta'| < \\epsilon} n R(\\theta', \\hat{\\theta}) \\geq \\frac{1}{I(\\theta)} $$\n",
    "\n",
    "This says that, in a local, large sample sense, the MLE is minimax.  It can also be shown that the MLE is approximately the Bayes rule.\n",
    "\n",
    "In summary, in parametric models with large samples, the MLE is approximately minimax and Bayes.  There is a caveat:  these results break down when the number of parameters is large."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 13.6  Admissibility"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An estimator $\\hat{\\theta}$ is **inadmissible** if there exists another rule $\\hat{\\theta}'$ such that\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "R(\\theta, \\hat{\\theta}') \\leq R(\\theta, \\hat{\\theta}) & \\quad \\text{for all } \\theta \\text{ and} \\\\\n",
    "R(\\theta, \\hat{\\theta}') < R(\\theta, \\hat{\\theta}) & \\quad \\text{for at least one } \\theta\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A prior has **full support** if for every $\\theta$ and every $\\epsilon > 0$, $\\int_{\\theta - \\epsilon}^{\\theta + \\epsilon} \\pi(\\theta) d\\theta > 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 13.20 (Bayes' rules are admissible)**. Suppose that $\\theta \\subset \\mathbb{R}$ and that $R(\\theta, \\hat{\\theta})$ is a continuous function of $\\theta$ for every $\\hat{\\theta}$.  Let $\\pi$ be a prior density with full support and let $\\hat{\\theta}^\\pi$ be the Bayes' rule.  If the Bayes risk is finite then $\\hat{\\theta}^\\pi$ is admissible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proof**.  Suppose $\\hat{\\theta}^\\pi$ is inadmissible.  Then there exists a better rule $\\hat{\\theta}$ such that $R(\\theta, \\hat{\\theta}) \\leq R(\\theta, \\hat{\\theta}^\\pi)$ for all $\\theta$ and $(\\theta_0, \\hat{\\theta}) < R(\\theta_0, \\hat{\\theta}^\\pi)$ for some $\\theta_0$.  Let $v = R(\\theta_0, \\hat{\\theta}^\\pi) - R(\\theta_0, \\hat{\\theta}) > 0$.  Since $R$ is continuous, there is an $\\epsilon > 0$ such that $R(\\theta, \\hat{\\theta}^\\pi) - R(\\theta_0, \\hat{\\theta}) > v/2$ for all $\\theta \\in (\\theta_0 - \\epsilon, \\theta_0 + \\epsilon)$.  Now,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "r(\\pi, \\hat{\\theta}^pi) - r(\\pi, \\hat{\\theta}) &= \\int R(\\theta, \\hat{\\theta}^\\pi) \\pi(\\theta) d\\theta - \\int R(\\theta, \\hat{\\theta}) \\pi(\\theta) d\\theta \\\\\n",
    "&= \\int \\left[ R(\\theta, \\hat{\\theta}^\\pi) - R(\\theta, \\hat{\\theta})\\right] \\pi(\\theta) d\\theta \\\\\n",
    "&\\geq \\int_{\\theta_0 - \\epsilon}^{\\theta_0 + \\epsilon} \\left[ R(\\theta, \\hat{\\theta}^\\pi) - R(\\theta, \\hat{\\theta})\\right] \\pi(\\theta) d\\theta \\\\\n",
    "& \\geq \\frac{v}{2} \\int_{\\theta_0 - \\epsilon}^{\\theta_0 + \\epsilon} \\pi(\\theta) d\\theta \\\\\n",
    "& > 0\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Hence $r(\\pi, \\hat{\\theta}^\\pi) > r(\\pi, \\hat{\\theta})$.  This implies that $\\hat{\\theta}^\\pi$ does not minimize $r(\\pi, \\hat{\\theta})$, which contradicts the fact that $\\hat{\\theta}^\\pi$ is the Bayes rule."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 13.21**.  Let $X_1, \\dots, X_n \\sim N(\\mu, \\sigma^2)$.  Under squared error loss, $\\overline{X}$ is admissible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The proof of this theorem is very technical and ommitted.  Outline:  the posterior mean is admissible for any strictly positive prior.  Take the prior to be $N(a, b^2)$.  When $b^2$ is very large, the posterior mean is approximately equal to $\\overline{X}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, a rule may be minimax, admissible, both, or neither.  Here are some facts linking admissibility and minimaxity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 13.22**.  Suppose that $\\hat{\\theta}$ has constant risk and is admissible.  Then it is minimax."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proof**.  The risk is $R(\\theta, \\hat{\\theta}) = c$ for some constant $c$.  If $\\hat{\\theta}$ were not minimax then there exists a rule $\\hat{\\theta}'$ that further reduces the risk,\n",
    "\n",
    "$$R(\\theta, \\hat{\\theta}') \\leq \\sup_{\\theta} R(\\theta, \\hat{\\theta}') < \\sup_{\\theta} R(\\theta, \\hat{\\theta}) = c$$\n",
    "\n",
    "But that would imply that $\\hat{\\theta}$ is inadmissible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 13.23**.  Let $X_1, \\dots, X_n \\sim N(\\theta, 1)$.  Then, under squared error loss, $\\hat{\\theta} = \\overline{X}$ is minimax."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Proof**.  According to Theorem 13.21, $\\hat{\\theta}$ is admissible.  The risk of $\\hat{\\theta}$ is $1/n$ which is constant.  According to Theorem 13.22, it is also minimax."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although miminax rules are not guaranteed to be admissible they are \"close to admissible\".  Say $\\hat{\\theta}$ is **strongly inadmissible** if there is a rule $\\hat{\\theta}'$ and an $\\epsilon > 0$ such that $R(\\theta, \\hat{\\theta}') < R(\\theta, \\hat{\\theta}) - \\epsilon$ for all $\\theta$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Theorem 13.24**.  If $\\hat{\\theta}$ is minimax then it is not strongly inadmissible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 13.7 Stein's Paradox"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose that $X \\sim N(0, 1)$ and consider estimating $\\theta$ with squared error loss.  From the previous section we know that $\\hat{\\theta}(X) = X$ is admissible.  \n",
    "\n",
    "Now consider estimating two, unrelated quantities $\\theta = (\\theta_1, \\theta_2)$ and suppose that $X_1 \\sim N(\\theta_1, 1)$ and $X_2 \\sim N(\\theta_2, 1)$ independently, with loss $L(\\theta, \\hat{\\theta}) = \\sum_{j=1}^2 (\\theta_j - \\hat{\\theta}_j)^2$.  Not surprisingly, $\\hat{\\theta}(X) = X$ is again admissible, where $X = (X_1, X_2)$.\n",
    "\n",
    "Now consider the generalization to $k$ normal means.  Let $\\theta = (\\theta_1, \\dots, \\theta_k)$, $X = (X_1, \\dots, X_k)$ with $X_i \\sim N(\\theta_i, 1)$ (independently) and loss $L(\\theta, \\hat{\\theta}) = \\sum_{j=1}^k (\\theta_j - \\hat{\\theta}_j)^2$.\n",
    "\n",
    "Stein astounded everyone when he proved that, if $k \\geq 3$, then $\\hat{\\theta}(X) = X$ is inadmissible.  It can be shown that the following estimator, known as the James-Stein estimator, has smaller risk:\n",
    "\n",
    "$$\\hat{\\theta}_S(X) = \\left(1 - \\frac{k-2}{\\sum_i X_i^2} \\right)^+ X_i$$\n",
    "\n",
    "where $(z)^+ = \\max \\{z, 0\\}$.  This estimator shrinks the $X_i$'s towards 0.  The message is that, when estimating many parameters, there is great value in \"shrinking\" the estimates.  This observation plays an important role in modern nonparametric function estimation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 13.9 Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**13.9.1**.  In each of the following models, find the Bayes risk and the Bayes estimator, using squared error loss.\n",
    "\n",
    "**(a)** $X \\sim \\text{Binomial}(n, p)$, $p \\sim \\text{Beta}(\\alpha, \\beta)$.\n",
    "\n",
    "**(b)** $X \\sim \\text{Poisson}(\\lambda)$, $\\lambda \\sim \\text{Gamma}(\\alpha, \\beta)$.\n",
    "\n",
    "**(c)** $X \\sim N(\\theta, \\sigma^2)$ where $\\sigma^2$ is known and $\\theta \\sim N(a, b^2)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**\n",
    "\n",
    "We can determine the posterior distribution, since its proportional to the likelihood times the prior, $f(\\theta | X) \\propto \\mathcal{L}(\\theta | X) f(\\theta)$.  But when using the square loss, the Bayes estimator is the mean of the posterior (from Theorem 13.9).  So we can then get the Bayes estimator as $\\hat{\\theta} = \\mathbb{E}[\\theta | X]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(a)** \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "f(\\theta | X) &\\propto \\mathcal{L}(\\theta | X) f(\\theta) \\\\\n",
    "&= \\prod_{i = 1}^N \\left( \\binom{N}{X_i} p^{X_i}(1 - p)^{n - X_i} \\right) \\frac{p^{\\alpha - 1}(1-p)^{\\beta - 1}}{B(\\alpha, \\beta)} \\\\\n",
    "& \\propto p^{(\\alpha - 1) + \\sum_i X_i} (1 - p)^{(\\beta - 1) + \\sum_i (n - X_i)} \\\\\n",
    "&= p^{\\left(\\alpha + N \\overline{X}_N \\right) - 1} (1 - p)^{\\left(\\beta + N(n - \\overline{X}_N) \\right)- 1}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "So the posterior is proportional to, and drawn from, a Beta distribution:\n",
    "\n",
    "$$ \\theta | X \\sim \\text{Beta}\\left(\\alpha + N \\overline{X}, \\beta + N(n - \\overline{X}) \\right)$$\n",
    "\n",
    "The mean of the Beta distribution with parameters $\\alpha_p, \\beta_p$ is $\\alpha_p / (\\alpha_p + \\beta_p)$, so the posterior mean (and the Bayes estimator) is:\n",
    "\n",
    "$$ \\hat{\\theta}(X) = \\frac{\\alpha + N \\overline{X}}{\\alpha + \\beta + Nn} $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Bayes risk for an arbitrary estimator $\\tilde{\\theta}$, given a prior $\\pi(\\theta)$ that $\\theta \\sim \\text{Beta}(\\alpha, \\beta)$ is:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "r(\\pi, \\tilde{\\theta}) &= \\int R(\\theta, \\tilde{\\theta}) \\pi(\\theta) d\\theta \\\\\n",
    "&= \\int (\\tilde{\\theta} - \\theta)^2 \\frac{\\theta^{\\alpha - 1}(1 - \\theta)^{\\beta - 1}}{B(\\alpha, \\beta)} d\\theta \\\\\n",
    "&= \\tilde{\\theta}^2 \\int \\frac{\\theta^{\\alpha - 1}(1 - \\theta)^{\\beta - 1}}{B(\\alpha, \\beta)} d\\theta\n",
    "- 2 \\tilde{\\theta} \\int \\theta \\frac{\\theta^{\\alpha - 1}(1 - \\theta)^{\\beta - 1}}{B(\\alpha, \\beta)} d\\theta\n",
    "+ \\int \\theta^2 \\frac{\\theta^{\\alpha - 1}(1 - \\theta)^{\\beta - 1}}{B(\\alpha, \\beta)} d\\theta \\\\\n",
    "&= \\tilde{\\theta}^2 \\mathbb{E}_\\pi[1] - 2 \\tilde{\\theta} \\mathbb{E}_\\pi[\\theta] + \\mathbb{E}_\\pi[\\theta^2]  \\\\\n",
    "&= \\tilde{\\theta}^2 \\cdot 1 - 2 \\tilde{\\theta} \\mathbb{E}_\\pi[\\theta] + \\mathbb{V}_\\pi[\\theta] + \\mathbb{E}_\\pi[\\theta]^2 \\\\ \n",
    "&= \\tilde{\\theta}^2 - 2 \\tilde{\\theta} \\frac{\\alpha}{\\alpha + \\beta}\n",
    "+ \\frac{\\alpha\\beta + \\alpha^2(\\alpha + \\beta + 1)}{(\\alpha + \\beta)^2(\\alpha + \\beta + 1)}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Including the observations from $X$, the prior is modified to a different Beta distribution, replacing $\\alpha$ by $\\alpha + N \\overline{X}$ and $\\beta$ by $\\beta + N(n - \\overline{X})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(b)**\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "f(\\theta | X) &\\propto \\mathcal{L}(\\theta | X) f(\\theta) \\\\\n",
    "&= \\prod_{i = 1}^N \\left( \\frac{\\lambda^X_i e^{-\\lambda}}{X_i!} \\right) \\frac{\\beta^\\alpha}{\\Gamma(\\alpha)}\\lambda^{\\alpha - 1}e^{-\\beta \\lambda} \\\\\n",
    "& \\propto \\lambda^{\\sum_i X_i} e^{-N \\lambda} \\lambda^{\\alpha - 1} e^{-\\beta \\lambda} \\\\\n",
    "&= \\lambda^{\\left(\\alpha + N \\overline{X}\\right) - 1} e^{-(\\beta + N) \\lambda}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "So the posterior is proportional to, and drawn from, a Gamma distribution:\n",
    "\n",
    "$$ \\theta | X \\sim \\text{Gamma}\\left(\\alpha + N \\overline{X}, \\beta + N \\right)$$\n",
    "\n",
    "The mean of the Gamma distribution with parameters $\\alpha_p, \\beta_p$ is $\\alpha_p / \\beta_p$, so the posterior mean (and the Bayes estimator) is:\n",
    "\n",
    "$$ \\hat{\\theta}(X) = \\frac{\\alpha + N \\overline{X}}{\\beta + N} $$\n",
    "\n",
    "The Bayes risk for an arbitrary estimator $\\tilde{\\theta}$, given a prior $\\pi(\\theta)$ that $\\theta \\sim \\text{Gamma}(\\alpha, \\beta)$ is:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "r(\\pi, \\tilde{\\theta}) &= \\int R(\\theta, \\tilde{\\theta})\\pi(\\theta) d\\theta \\\\\n",
    "&= \\int (\\tilde{\\theta} - \\theta)^2 \\frac{\\beta^\\alpha}{\\Gamma(\\alpha)} \\theta^{\\alpha - 1}e^{-\\beta \\theta} d\\theta \\\\\n",
    "&= \\tilde{\\theta}^2 \\int \\frac{\\beta^\\alpha}{\\Gamma(\\alpha)} \\theta^{\\alpha - 1}e^{-\\beta \\theta} d\\theta\n",
    "- 2 \\tilde{\\theta} \\int \\theta \\frac{\\beta^\\alpha}{\\Gamma(\\alpha)} \\theta^{\\alpha - 1}e^{-\\beta \\theta} d\\theta\n",
    "+ \\int \\theta^2\\frac{\\beta^\\alpha}{\\Gamma(\\alpha)} \\theta^{\\alpha - 1}e^{-\\beta \\theta} d\\theta \\\\\n",
    "&= \\tilde{\\theta}^2 \\mathbb{E}_\\pi[1] - 2 \\tilde{\\theta} \\mathbb{E}_\\pi[\\theta] + \\mathbb{E}_\\pi[\\theta^2] \\\\\n",
    "&= \\tilde{\\theta}^2 \\cdot 1 - 2 \\tilde{\\theta} \\mathbb{E}_\\pi[\\theta] + \\mathbb{V}_\\pi[\\theta] + \\mathbb{E}_\\pi[\\theta]^2 \\\\\n",
    "&= \\tilde{\\theta}^2 - 2 \\tilde{\\theta} \\frac{\\alpha}{\\beta} + \\frac{\\alpha(\\alpha + 1)}{\\beta^2}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Including the observations from $X$, the prior is modified to a different Gamma distribution, replacing $\\alpha$ by $\\alpha + N \\overline{X}$ and $\\beta$ by $\\beta + N$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(c)**\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "f(\\theta | X) &\\propto \\mathcal{L}(\\theta | X) f(\\theta) \\\\\n",
    "&= \\prod_{i = 1}^N \\left( \\frac{1}{\\sigma \\sqrt{2 \\pi}} e^{-\\frac{1}{2} \\left(\\frac{X_i - \\theta}{\\sigma} \\right)^2} \\right) \\frac{1}{a \\sqrt{2 \\pi}} e^{-\\frac{1}{2} \\left( \\frac{\\theta - a}{b}\\right)^2} \\\\\n",
    "& \\propto \\exp \\left\\{ -\\frac{1}{2} \\left( \\sum_{i=1}^N \\left( \\frac{X_i - \\theta}{\\sigma}\\right)^2 + \\left( \\frac{\\theta - a}{b}\\right)^2 \\right) \\right\\} \\\\\n",
    "& = \\exp \\left\\{-\\frac{1}{2} \\left( \\theta^2 \\left(\\frac{N}{\\sigma^2} + \\frac{1}{b^2} \\right) - 2 \\theta \\left( \\frac{N \\overline{X}}{\\sigma^2} + \\frac{a}{b^2} \\right) + C\\right) \\right\\} \\\\\n",
    "& \\propto \\exp \\left\\{-\\frac{1}{2} \\left( \\frac{Nb^2 + \\sigma^2}{\\sigma^2 b^2} \\right)\\left(\\theta^2 - 2\\theta \\left(\\frac{N b^2 \\overline{X} + a \\sigma^2}{N b^2 + \\sigma^2} \\right) \\right) \\right\\} \\\\\n",
    "& \\propto \\exp \\left\\{-\\frac{1}{2} \\left( \\frac{Nb^2 + \\sigma^2}{\\sigma^2 b^2} \\right) \\left(\\theta - \\left(\\frac{N b^2 \\overline{X} + a \\sigma^2}{N b^2 + \\sigma^2} \\right) \\right)^2 \\right\\} \\\\\n",
    "& = \\exp \\left\\{-\\frac{1}{2} \\left(\\frac{\\theta - \\left(\\frac{N b^2 \\overline{X} + a \\sigma^2}{N b^2 + \\sigma^2} \\right)}{\\sqrt{\\frac{\\sigma^2 b^2}{Nb^2 + \\sigma^2}}} \\right)^2 \\right\\} \\\\\n",
    "& = \\exp \\left\\{-\\frac{1}{2} \\left(\\frac{\\theta - \\overline{\\theta}}{w} \\right)^2 \\right\\}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where \n",
    "\n",
    "$$\n",
    "\\overline{\\theta} = \\frac{\\sigma^2}{\\sigma^2 + N b^2}a + \\frac{Nb^2}{\\sigma^2 + Nb^2}\\overline{X}\n",
    "\\quad \\text{and} \\quad\n",
    "w^2 = \\frac{1}{N} \\frac{1}{\\frac{1}{\\sigma^2} + \\frac{1}{Nb^2}}\n",
    "$$\n",
    "\n",
    "So the posterior is proportional to, and drawn from, a Normal distribution:\n",
    "\n",
    "$$ \\theta | X \\sim N(\\overline{\\theta}, w^2) $$\n",
    "\n",
    "The mean of the Normal distribution is $\\overline{\\theta}$, so the mean of the posterior (and the Bayes estimator) is:\n",
    "\n",
    "$$ \\hat{\\theta}(X) = \\frac{\\sigma^2}{\\sigma^2 + N b^2}a + \\frac{Nb^2}{\\sigma^2 + Nb^2}\\overline{X}$$\n",
    "\n",
    "The Bayes risk for an arbitrary estimator $\\tilde{\\theta}$, given a prior $\\pi(\\theta)$ that $\\theta \\sim N(a, b^2)$ is:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "r(\\pi, \\tilde{\\theta}) &= \\int R(\\theta, \\tilde{\\theta}) \\pi(\\theta) d\\theta \\\\\n",
    "&= \\int (\\tilde{\\theta} - \\theta)^2 \\frac{1}{b \\sqrt{2 \\pi}} \\exp \\left\\{ -\\frac{1}{2} \\left(\\frac{\\theta - a}{b}\\right)^2 \\right\\} d\\theta \\\\\n",
    "&= \\tilde{\\theta}^2 \\int \\frac{1}{b \\sqrt{2 \\pi}} \\exp \\left\\{ -\\frac{1}{2} \\left(\\frac{\\theta - a}{b}\\right)^2 \\right\\} d\\theta\n",
    "- 2 \\tilde{\\theta} \\int \\theta \\frac{1}{b \\sqrt{2 \\pi}} \\exp \\left\\{ -\\frac{1}{2} \\left(\\frac{\\theta - a}{b}\\right)^2 \\right\\} d\\theta\n",
    "+ \\int \\theta^2 \\frac{1}{b \\sqrt{2 \\pi}} \\exp \\left\\{ -\\frac{1}{2} \\left(\\frac{\\theta - a}{b}\\right)^2 \\right\\} d\\theta \\\\\n",
    "&= \\tilde{\\theta}^2 \\mathbb{E}_\\pi[1] - 2 \\tilde{\\theta} \\mathbb{E}_\\pi[\\theta] + \\mathbb{E}_\\pi[\\theta^2] \\\\\n",
    "&= \\tilde{\\theta}^2 - 2 \\tilde{\\theta} \\mathbb{E}_\\pi[\\theta] + \\mathbb{V}_pi[\\theta] + \\mathbb{E}_\\pi[\\theta]^2 \\\\\n",
    "&= \\tilde{\\theta}^2 - 2 \\tilde{\\theta} a + a^2 + b^2\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Including the observations from $X$, the prior is modified to a different Normal distribution, replacing $a$ with $\\overline{\\theta}$ and $b^2$ with $w^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 13.9.2**.  Let $X_1, \\dots, X_n \\sim N(\\theta, \\sigma^2)$ and suppose we estimate $\\theta$ with loss function $L(\\theta, \\hat{\\theta}) = (\\theta - \\hat{\\theta})^2 / \\sigma^2$.  Show that $\\overline{X}$ is admissible and minimax."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "The risk for an estimator $\\hat{\\theta}$ is the risk for the same estimator using the mean squared error, but scaled by $1 / \\sigma^2$:\n",
    "\n",
    "$$\n",
    "R_L(\\theta, \\hat{\\theta}) = \\mathbb{E}_\\theta[L(\\theta, \\hat{\\theta})] \n",
    "= \\mathbb{E}_\\theta\\left[\\left( \\frac{\\theta - \\hat{\\theta}}{\\sigma} \\right)^2 \\right] \n",
    "= \\frac{1}{\\sigma^2} \\mathbb{E}_\\theta[ L_\\text{MSE}(\\theta, \\hat{\\theta})]  \n",
    "= \\frac{1}{\\sigma^2} R_\\text{MSE}(\\theta, \\hat{\\theta})\n",
    "$$\n",
    "\n",
    "Since $\\overline{X}$ is admissible and minimax for the MSE loss function, it follows it is also admissible and minimax for this rescaled loss function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 13.9.3**.  Let $\\Theta = \\{ \\theta_1, \\dots, \\theta_k \\}$ be a finite parameter space.  Prove that the posterior mode is the Bayes estimator under zero-one loss."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.  The Bayes rule minimizes \n",
    "\n",
    "$$r(\\hat{\\theta} | x) = \\int L(\\theta, \\hat{\\theta}) f(\\theta | x) d\\theta = \\sum_i I(\\hat{\\theta} \\neq \\theta_i) f(\\theta_i | x) = 1 - \\sum_i I(\\hat{\\theta} = \\theta_i) f(\\theta_i | x) = 1 - f(\\hat{\\theta} | X)$$\n",
    "\n",
    "The posterior probability $f(\\theta | X)$ is maximized on the posterior mode, $\\hat{\\theta} = \\text{argmax}_\\theta f(\\theta | X)$, so the posterior mode is also the Bayes estimator under zero-one loss."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 13.9.4 (Casella and Berger)**.  Let $X_1, \\dots, X_n$ be a sample from a distribution with variance $\\sigma^2$.  Consider estimators of the form $bS^2$ where $S^2$ is the sample variance.  Let the loss function for estimating $\\sigma^2$ be\n",
    "\n",
    "$$L(\\sigma^2, \\hat{\\sigma}^2) = \\frac{\\hat{\\sigma}^2}{\\sigma^2} - 1 - \\log \\frac{\\hat{\\sigma^2}}{\\sigma^2}$$\n",
    "\n",
    "Find the optimal value of $b$ that minimizes the risk for all $\\sigma^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.  For an estimator of the form $\\hat{\\sigma}^2 = bS^2$, the risk is\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "R(\\sigma^2, bS^2) &= \\mathbb{E}_{\\sigma^2}[ L(\\sigma^2, bS^2) ]\\\\\n",
    "&=\\mathbb{E}_{\\sigma^2}\\left[\\frac{bS^2}{\\sigma^2} - 1 - \\log b - \\log S^2 + \\log \\sigma^2 \\right] \\\\\n",
    "&= \\frac{b}{\\sigma^2} \\mathbb{E}_{\\sigma^2}[S^2] - 1 -\\log b - \\mathbb{E}_{\\sigma^2}[\\log S^2] + \\log \\sigma^2 \\\\\n",
    "&= b - \\log b + C\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where $C$ does not depend on $b$.  The risk is minimized when the derivative is zero, $1 - 1/b = 0$, or $b = 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 13.9.5 (Berliner, 1983)**.  Let $X \\sim \\text{Binomial}(n, p)$ and suppose the loss function is\n",
    "\n",
    "$$ L(p, \\hat{p}) = \\left(1 - \\frac{\\hat{p}}{p} \\right)^2$$\n",
    "\n",
    "where $0 < p < 1$.  Consider the estimator $\\hat{p}(X) = 0$.  This estimator falls outside of the parameter space $(0, 1)$ but we will allow this.  Show that $\\hat{p}(X) = 0$ is the unique, minimax rule."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "For the estimator $\\hat{p}(X) = 0$, the loss function is always 1, and so the risk is also always 1.\n",
    "\n",
    "For any estimator $\\tilde{p}$ that falls within the parameter space $(0, 1)$, we are interested in calculating its maximum risk, $\\overline{R}(p, \\tilde{p}) = \\sup_p R(p, \\tilde{p})$.\n",
    "\n",
    "$$ \n",
    "\\begin{align}\n",
    "R(p, \\tilde{p}) &= \\mathbb{E}_p \\left[\\left(1 - \\frac{\\tilde{p}(X)}{p} \\right)^2\\right]\\\\\n",
    "&=\\mathbb{E}_p \\left[1 + \\frac{1}{p} (-2 \\tilde{p}(X)) + \\frac{1}{p^2} \\tilde{p}(X)^2 \\right] \\\\\n",
    "&= 1 + \\frac{1}{p} (-2 \\mathbb{E}_p[\\tilde{p}(X)]) + \\frac{1}{p^2} \\mathbb{E}_p[\\tilde{p}(X)^2]\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "For a specific estimator, we have the expectations as constants, $a = \\mathbb{E}_p[\\tilde{p}(X)]$ and $b = \\mathbb{E}_p[\\tilde{p}(X)^2]$.  We are interested in picking a $p \\in (0, 1)$ such that this risk is greater than 1; we can always do $p = a / b < 2a / b$, since $a > 0$ due to the open interval constraint:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "1 + \\frac{1}{p} (-2a) + \\frac{1}{p^2} b > 1 \\\\\n",
    "\\frac{1}{p} \\left(-2a + \\frac{1}{p} b \\right) > 0 \\\\\n",
    "p < \\frac{2a}{b}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Therefore, any estimator within the parameter space can produce a risk greater than 1, and it is not minimax as the estimator $\\hat{p}(X)$ produces a smaller maximum risk."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 13.9.6 (Computer Experiment)**.  Compare the risk of the MLE and the James-Stein estimator by simulation.  Try various values of $n$ and various vectors $\\theta$.  Summarize your results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**.\n",
    "\n",
    "Let's plot the mean square error from running the simulation a large number of times for $k = 3, 10, 100, 1000$, assuming that $\\theta$ is a vector of all ones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from tqdm import tqdm_notebook\n",
    "\n",
    "\n",
    "# MLE estimator: just return the original values\n",
    "def mle(X):\n",
    "    return X\n",
    "\n",
    "# James-Stein estimator: shrink towards zero\n",
    "def james_stein(X):\n",
    "    assert len(X) > 2\n",
    "    return np.maximum(1 - ((len(X) - 2) / sum(X**2)), 0) * X\n",
    "\n",
    "# Calculates the mean square error between two sequences\n",
    "def mse(X_pred, X):\n",
    "    return ((X - X_pred)**2).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_range(k, B=100000):\n",
    "    error_mle = np.empty(B)\n",
    "    error_js = np.empty(B)\n",
    "    for i in tqdm_notebook(range(B)):\n",
    "        theta = np.ones(k)\n",
    "        X = norm.rvs(loc=1, scale=1, size=k)\n",
    "        theta_hat_mle = mle(X)\n",
    "        theta_hat_js = james_stein(X)\n",
    "        error_mle[i] = mse(theta_hat_mle, theta)\n",
    "        error_js[i] = mse(theta_hat_js, theta)\n",
    "\n",
    "    plt.figure(figsize=(12, 8))\n",
    "    plt.hist(error_mle, density=True, bins=100, histtype='step', color='blue', label='MLE, k = ' + str(k))\n",
    "    plt.hist(error_js, density=True, bins=100, histtype='step', color='red', label='James-Stein, k = ' + str(k))\n",
    "    plt.legend(loc='upper right')\n",
    "    plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1a55ac7a4c88417b9210e826735c0ad7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=100000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAHSCAYAAADmLK3fAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dfZSeZX0v+u8PEkFFYRfirhCYoauoYEgCZ3hp9ZyNRQpRhF2PbZEctL6U0gVK3+X0RbtWl13uVc7Bo6gUPWz2tog9KiJVfKlbPNTVVTVsMYVSlOJEIxwNdEt9gULwOn/MkIbJk2SSPDPPzFyfz1r3Yu6X575/yb0Sv175PddVrbUAAEBv9ht1AQAAMAqCMAAAXRKEAQDokiAMAECXBGEAALokCAMA0KVlo3rwYYcd1sbHx0f1eAAAOnHbbbc90FpbMfP4yILw+Ph4NmzYMKrHAwDQiaraNOi41ggAALokCAMA0CVBGACALo2sRxgAYC489thj2bx5cx555JFRl8I8O/DAA7Ny5cosX758VtcLwgDAkrJ58+Y84xnPyPj4eKpq1OUwT1prefDBB7N58+YcffTRs/qM1ggAYEl55JFHcuihhwrBnamqHHrooXv0LwGCMACw5AjBfdrT9y4IAwAMWVXlggsu2La/devWrFixImeffXaS5Nprr80ll1yyw+fGx8dz/PHHZ+3atVm7dm3e+MY37vI5f/zHf5zLL798aHU/8sgjOfnkk7NmzZo8//nPz1ve8pah3Xsh2m2PcFVdk+TsJN9tra0acH59kjdN7/4gya+31r461CoBAPbS+HiyaeByCntnbCyZnNz1NU9/+tNzxx135OGHH85Tn/rU/PVf/3WOOOKIWd3/lltuyWGHHbbvhe6FAw44IJ/73Ody0EEH5bHHHssLX/jCrFu3LqeeeupI6plrsxkRvjbJWbs4/40k/6G1tjrJnyS5egh1AQAMxaZNSWvD22YbqtetW5dPfOITSZLrr78+r3zlK+fwV5m8973vzbp16/Lwww/v9T2qKgcddFCSqdk3HnvssSXdZrLbINxauzXJP+/i/N+21v7H9O7fJVk5pNoAABat8847Lx/84AfzyCOPZOPGjTnllFNm9bkXvehF21ojrrjiill95sorr8xf/dVf5cYbb8xTn/rUJ5277rrrtt1v++0Vr3jFwHs9/vjjWbt2bZ71rGfljDPOmHXdi9Gwp097XZJPDvmeAACLzurVqzM5OZnrr78+L3nJS2b9uT1tjXj/+9+flStX5sYbbxw4f+769euzfv36Wd9v//33z+23357vfe97+YVf+IXccccdWbVqh+7YJWFoQbiqXpSpIPzCXVxzYZILk+Soo44a1qMBABakc845J7/zO7+Tz3/+83nwwQfn5BmrVq3K7bffvtP5c6+77rr82Z/92Q7Hf/qnfzof/vCHd3rfQw45JKeddlo+9alPCcK7UlWrk7wvybrW2k7fcmvt6kz3EE9MTLRhPBsAYKF67Wtfm4MPPjjHH398Pv/5z+/1fa688sokGTjTxAknnJBf//VfzznnnJNPf/rTOfzww590fk9GhLds2ZLly5fnkEMOycMPP5zPfvazedOb3rT7Dy5S+zx9WlUdleSGJBe01r627yUBACwNK1euzKWXXjrw3LXXXpuVK1du2zZv3pzkyT3Cr3rVq5Ik//iP/5hDDz10p8954QtfmMsvvzwvfelL88ADD+x1vffff39e9KIXZfXq1TnppJNyxhlnbJvybSmq1nY9MFtV1yc5LclhSb6T5C1JlidJa+2qqnpfkv81yRPfodzaWpvY3YMnJibahg0b9r5yAIAB7rrrrhx77LHb9kcxfdqwnX322bnhhhvylKc8ZX4fvAjNfP9JUlW3Dcqnu22NaK3tcq6P1trrk7x+T4sEAJgP8x1a58LHP/7xUZewJFlZDgCALgnCAAB0SRBeasbHk6odt/HxUVcGALCgDHtBDUbtiXUkZ1rCyyMCAOwNI8IAAHRJEAYAGLKDDjpopM//zne+k7PPPjtr1qzJcccdt22J58nJyXzgAx/Y7efvu+++vOIVr9jr55922mkZ5jS5X/rSl7bNrbxmzZp89KMfHcp9tUYAAEvbUphIeA+9+c1vzhlnnLFtMY+NGzcm+bcgfP755+/y84cffvgul1+eb6tWrcqGDRuybNmy3H///VmzZk1e9rKXZdmyfYuyRoQBgKXtie/PDGubZaj+wQ9+kNNPPz0nnnhijj/++HzsYx9LMhVGn/e85+X1r399Vq1alfXr1+ezn/1sXvCCF+SYY47Jl770pSTJD3/4w7z2ta/NSSedlBNOOGHb5++8886cfPLJWbt2bVavXp2vf/3rOzz7/vvvz8qVK7ftr169Okly2WWX5W/+5m+ydu3aXHHFFXn88cfzu7/7uznppJOyevXq/Pmf//m2GletWpVkagW8l7/85TnrrLNyzDHH5Pd+7/dm/Vv/4x//OK9+9avzh3/4h7P+zCBPe9rTtoXeRx55JDWk7z4ZEQYAmAMHHnhgPvrRj+aZz3xmHnjggZx66qk555xzkiT33HNPPvShD+Xqq6/OSSedlA984AP5whe+kJtuuil/+qd/mhtvvDFvfetb83M/93O55ppr8r3vfS8nn3xyXvziF+eqq67KpZdemvXr1+fRRx/N448/vsOzL7744vzyL/9yrrzyyrz4xS/Oa17zmhx++OF529velssvv3zbAh1XX311Dj744Hz5y1/Ov/7rv+YFL3hBfv7nf36HoHn77bfnK1/5Sg444IA897nPzRve8IYceeSRu/z1b926NevXr8+qVavyB3/wBzuc/83f/M3ccsstOxw/77zzctlll+1w/Itf/GJe+9rXZtOmTXn/+9+/z6PBiSAMADAnWmv5/d///dx6663Zb7/98u1vfzvf+c53kiRHH310jj/++CTJ85///Jx++umpqhx//PGZnG67+MxnPpObbropl19+eZKpkdBvfvOb+Zmf+Zm89a1vzebNm/Pyl788xxxzzA7PPvPMM3PvvffmU5/6VD75yU/mhBNOyB133LHDdZ/5zGeycePGbW0QDz30UL7+9a/nOc95zpOuO/3003PwwQcnSY477rhs2rRpt0H4137t1/JLv/RLA0NwklxxxRW7/PxMp5xySu68887cddddefWrX51169blwAMP3KN7zCQIAwDMgeuuuy5btmzJbbfdluXLl2d8fDyPPPJIkuSAAw7Ydt1+++23bX+//fbL1q1bk0wF6Y985CN57nOf+6T7HnvssTnllFPyiU98ImeeeWbe97735a677sp73/veJMnNN9+cww8/PD/xEz+R888/P+eff37OPvvs3HrrrTn00EOfdK/WWt75znfmzDPPfNLxyRk90NvXu//++2+rcVd+9md/Nrfcckt++7d/e2Bg3dMR4e1//U9/+tNzxx13ZGJiYrd17IoeYQCAOfDQQw/lWc96VpYvX55bbrklm/bwC3tnnnlm3vnOd6ZNrw/wla98JUly77335qd+6qfyxje+Meecc042btyYiy++OLfffntuv/32HH744fnc5z6XH/3oR0mS73//+/mnf/qnHHXUUXnGM56R73//+096xnve85489thjSZKvfe1r+eEPfzjrGl/1qldt62me6XWve11e8pKX5Bd/8RcHBucrrrhiW83bb4NC8De+8Y1t99i0aVPuvvvujA9hsTAjwgAAQ7R169YccMABWb9+fV72spdlYmIia9euzfOe97w9us8f/dEf5Td+4zeyevXqtNYyPj6ej3/84/nLv/zL/MVf/EWWL1+en/zJn8yb3/zmHT5722235ZJLLsmyZcvy4x//OK9//etz0kkn5bHHHsuyZcuyZs2a/Mqv/EouvfTSTE5O5sQTT0xrLStWrMiNN9446xo3btyYZz/72Ts9/1u/9Vt56KGHcsEFF+S6667Lfvvt3RjsF77whbztbW/L8uXLs99+++Xd7353DjvssL261/aqDVqFbB5MTEy0Yc4vx7Sqna8sN6J3DQDz6a677sqxxx77bwfmefq0r371q/nVX/3VnY6ULhX/8i//kte97nX50Ic+NOpSnmSH95+kqm5rre3QR2FEGABY2uZxzt+rrroq73jHO/L2t7993p45Ks985jMXXAjeU4IwAMCQXHTRRbnoootGXQaz5MtyAAB0SRAGAJacUX0HitHa0/cuCAMAS8qBBx6YBx98UBjuTGstDz744B4tsqFHGABYUlauXJnNmzdny5Ytoy6FeXbggQdm5cqVs75eEAYAlpTly5fn6KOPHnUZLAJaIwAA6JIgDABAlwRhAAC6JAgDANAlQRgAgC4JwgAAdEkQBgCgS4IwAABdEoQBAOiSIAwAQJcEYQAAuiQIAwDQJUEYAIAuCcIAAHRJEAYAoEuCMAAAXRKEAQDokiAMAECXBGEAALokCAMA0CVBGACALgnCAAB0SRAGAKBLgjAAAF0ShAEA6JIgDABAlwRhAAC6JAgDANAlQRgAgC4JwgAAdEkQBgCgS4IwAABdEoQBAOiSIAwAQJcEYQAAuiQIAwDQJUEYAIAuCcIAAHRJEAYAoEuCMAAAXdptEK6qa6rqu1V1x07OV1W9o6ruqaqNVXXi8MsEAIDhms2I8LVJztrF+XVJjpneLkzynn0vCwAA5tZug3Br7dYk/7yLS85N8l/blL9LckhVPXtYBQIAwFwYRo/wEUm+td3+5uljAACwYA0jCNeAY23ghVUXVtWGqtqwZcuWITwaAAD2zjCC8OYkR263vzLJfYMubK1d3VqbaK1NrFixYgiPBgCAvTOMIHxTkldNzx5xapKHWmv3D+G+AAAwZ5bt7oKquj7JaUkOq6rNSd6SZHmStNauSnJzkpckuSfJj5K8Zq6KBQCAYdltEG6tvXI351uSi4dWEQAAzAMrywEA0CVBGACALgnCAAB0SRAGAKBLgjAAAF0ShAEA6JIgDABAlwRhAAC6JAgDANAlQRgAgC4JwgAAdEkQBgCgS4IwAABdEoQBAOiSIAwAQJcEYQAAuiQIAwDQJUEYAIAuCcIAAHRJEAYAoEuCMAAAXRKEAQDokiAMAECXBGEAALokCAMA0CVBGACALgnCAAB0SRAGAKBLgjAAAF0ShAEA6JIgDABAlwRhAAC6JAgDANAlQRgAgC4JwgAAdEkQBgCgS4IwAABdEoQBAOiSIAwAQJcEYQAAuiQIAwDQJUEYAIAuCcIAAHRJEAYAoEuCMAAAXRKEAQDokiAMAECXBGEAALokCO+N8fGkasdtfHzUlQEAMEvLRl3AorRpU9Lajser5r8WAAD2ihFhAAC6JAgDANAlQRgAgC4JwgAAdEkQBgCgS4IwAABdEoQBAOiSIAwAQJcEYQAAuiQI786g5ZTHxkZdFQAA+8gSy7uzs+WUAQBY1IwIAwDQpVkF4ao6q6rurqp7quqyAecPrqq/qqqvVtWdVfWa4ZcKAADDs9sgXFX7J3lXknVJjkvyyqo6bsZlFyf5h9bamiSnJfk/quopQ64VAACGZjYjwicnuae1dm9r7dEkH0xy7oxrWpJnVFUlOSjJPyfZOtRKAQBgiGYThI9I8q3t9jdPH9velUmOTXJfkr9Pcmlr7cczb1RVF1bVhqrasGXLlr0sGQAA9t1sgnANODZzGoUzk9ye5PAka5NcWVXP3OFDrV3dWptorU2sWLFij4td8MbGdpxqrWpqCjYAABaU2QThzUmO3G5/ZaZGfrf3miQ3tCn3JPlGkucNp8RFZHJyaqq1mdumTaOuDACAGWYThL+c5JiqOnr6C3DnJblpxjXfTHJ6klTVv0/y3CT3DrNQAAAYpt0uqNFa21pVlyT5dJL9k1zTWruzqi6aPn9Vkj9Jcm1V/X2mWine1Fp7YA7rBgCAfTKrleVaazcnuXnGsau2+/m+JD8/3NIAAGDuWFkOAIAuCcIAAHRJEAYAoEuCMAAAXRKEAQDokiA8HwatOGe1OQCAkZrV9Gnso8nJHY/VoJWrAQCYL0aEAQDokiAMAECXBGEAALokCAMA0CVBGACALgnCAAB0SRAGAKBLgjAAAF0ShAEA6JIgDABAlwRhAAC6JAgDANAlQRgAgC4JwvtgfDyp2nEbH5/Fh8fG9uHDAADsq2WjLmAx27QpaW3H41Wz+PDk5ODjs/owAAD7yogwAABdEoRnYWctEGNjo64MAIC9pTXiCePjU70OM42N7bQFAgCAxUsQfsKu0q62XQCAJUdrxBwwIQQAwMJnRHgOmBACAGDhMyIMAECXBGEAALokCAMA0CVBGACALgnCAAB0SRAGAKBLgvA8Mr8wAMDCYR7heWR+YQCAhcOIMAAAXRKEF5pB/RN6JwAAhk5rxEIzqH9C7wQAwNAZEQYAoEuCMAAAXRKEAQDokiAMAECXBGEAALokCAMA0CVBGACALgnCAAB0SRAGAKBLgjAAAF0ShAEA6JIgvJ3x8aRqx21sbNSVAQAwbMtGXcBCsmlT0tqoqwAAYD4YEQYAoEuC8GIwNja4Z2N8fNSVAQAsWlojFoPJycHHq+a1DACApcSIMAAAXRKEAQDokiAMAECXBGEAALokCAMA0CVBeAEwOxoAwPybVRCuqrOq6u6quqeqLtvJNadV1e1VdWdV/b/DLXNpm5ycWtFu5rZp06grAwBYunY7j3BV7Z/kXUnOSLI5yZer6qbW2j9sd80hSd6d5KzW2jer6llzVTAAAAzDbEaET05yT2vt3tbao0k+mOTcGdecn+SG1to3k6S19t3hlgkAAMM1myB8RJJvbbe/efrY9p6T5N9V1eer6raqetWgG1XVhVW1oao2bNmyZe8qBgCAIZhNEB60jm+bsb8syf+U5KVJzkzyR1X1nB0+1NrVrbWJ1trEihUr9rhYAAAYlt32CGdqBPjI7fZXJrlvwDUPtNZ+mOSHVXVrkjVJvjaUKgEAYMhmMyL85STHVNXRVfWUJOcluWnGNR9L8j9X1bKqelqSU5LcNdxSAQBgeHY7Itxa21pVlyT5dJL9k1zTWruzqi6aPn9Va+2uqvpUko1Jfpzkfa21O+aycAAA2BfV2sx23/kxMTHRNmzYMJJnD1SVSsuIfjsGqsqu6xl0wc4+tNubAQAsTVV1W2ttYuZxK8sBANAlQRgAgC4JwgAAdEkQBgCgS4LwAjY2NvUdt5nb+PguLhgbG2XJAACLxmwW1GBEJicHH6/azQUAAOyWEWEAALokCAMA0CVBGACALgnCAAB0SRAGAKBLgjAAAF0ShAEA6JIgDABAlwRhAAC6JAgDANAlQRgAgC4JwgAAdEkQBgCgS4LwIjQ2llTtuI2Pj7oyAIDFY9moC2DPTU4OPl41r2UAACxqRoQBAOiSIAwAQJcEYQAAuiQIAwDQJUEYAIAuCcIAAHRJEAYAoEuCMAAAXRKEAQDokiAMAECXBGEAALokCAMA0KVloy6A4RkbS6oGn2vzWwoAwIInCC8hk5O7OLmTgAwA0CutEQAAdEkQBgCgS4IwAABdEoQBAOiSIAwAQJcEYQAAuiQIAwDQpT6D8Pj41MoT229jY6OuCgCAedTnghqbNiVtwFprFp0AAOhGnyPCAAB0TxAGAKBLgjAAAF0ShAEA6JIgDABAlwRhAAC6JAh3ZObUyVVTUyoDAPSoz3mEOzVw6mRzJwMAnTIiDABAlwRhAAC6JAgDANAlQRgAgC4JwgAAdEkQBgCgS4IwAABdEoQBAOiSIAwAQJdmFYSr6qyquruq7qmqy3Zx3UlV9XhVvWJ4JQIAwPDtNghX1f5J3pVkXZLjkryyqo7byXX/Kcmnh10kAAAM22xGhE9Ock9r7d7W2qNJPpjk3AHXvSHJR5J8d4j1AQDAnJhNED4iybe22988fWybqjoiyS8kuWp4pc2tqh23sbFRVwUAwHxZNotrasCxNmP/7Une1Fp7vGrQ5dM3qrowyYVJctRRR822xjnRZv4KAADoymyC8OYkR263vzLJfTOumUjywekQfFiSl1TV1tbajdtf1Fq7OsnVSTIxMSGKAgAwMrMJwl9OckxVHZ3k20nOS3L+9he01o5+4uequjbJx2eGYBamsbGptpBBxycn570cAIB5s9sg3FrbWlWXZGo2iP2TXNNau7OqLpo+v2j6gtnRzsLuLjpcAACWhNmMCKe1dnOSm2ccGxiAW2u/su9lAQDA3LKyHAAAXRKEAQDokiDciye+FTdzGx8fdWUAACMxqx5hlgDfigMAeBIjwgAAdEkQBgCgS4IwAABdEoQBAOiSIMxAJpkAAJY6s0YwkEkmAIClzogwAABdEoQBAOiSIAwAQJcEYQAAuiQIAwDQJUEYAIAuCcIAAHRJEAYAoEuCMAAAXRKEAQDokiAMAECXBGH2yNhYUrXjNj4+6soAAPbMslEXwOIyOTn4eNW8lgEAsM+MCAMA0CVBGACALgnCAAB0SRAGAKBLgjAAAF0ShAEA6JIgDABAlwRhAAC6JAgDANAlQbh3g9ZMtl4yANABSyz3btCaydZLBgA6YEQYAIAuCcIAAHRJEAYAoEuCMAAAXRKEGYpBk0+YgAIAWMjMGsFQDJp8IjEBBQCwcBkRBgCgS4IwAABdEoQBAOiSIAwAQJcEYXZkCggAoANmjWBHpoAAADpgRBgAgC4JwgAAdEkQZk5pNwYAFio9wswp7cYAwEJlRBgAgC4JwgAAdEkQBgCgS4IwAABdEoQBAOiSIAwAQJcEYQAAuiQIAwDQJUEYAIAuCcIAAHRJEGb2xsam1kbefhsfH9qt9uF2AAB7bFZBuKrOqqq7q+qeqrpswPn1VbVxevvbqloz/FIZucnJpLUnb5s2De1W+3A7AIA9ttsgXFX7J3lXknVJjkvyyqo6bsZl30jyH1prq5P8SZKrh10oAAAM02xGhE9Ock9r7d7W2qNJPpjk3O0vaK39bWvtf0zv/l2SlcMtEwAAhms2QfiIJN/abn/z9LGdeV2ST+5LUQAAMNeWzeKaGnCsDbyw6kWZCsIv3Mn5C5NcmCRHHXXULEsEAIDhm82I8OYkR263vzLJfTMvqqrVSd6X5NzW2oODbtRau7q1NtFam1ixYsXe1AsAAEMxmyD85STHVNXRVfWUJOcluWn7C6rqqCQ3JLmgtfa14ZfJgjXkedBMqwYAzJfdtka01rZW1SVJPp1k/yTXtNburKqLps9fleTNSQ5N8u6qSpKtrbWJuSubBWNycvDxGtRRM++3AwDYqWptYLvvnJuYmGgbNmwYybNTNTVpLXNnyL/HXhkAsLeq6rZBg7RWlgMAoEuCMAAAXRKEAQDokiAMAECXBGEAALokCLMomF8YABi22SyxDCNnfmEAYNiMCAMA0CVBmLkxqJdBHwMAsIBojWBuDOpl0McAACwgRoQBAOiSIMyiZjYJAGBvaY1gUTObBACwt4wIM38M3wIAC4gRYeaP4VsAYAExIgwAQJcEYZYkXRgAwO5ojWBJ0oUBAOyOEWEAALokCDN6+hgAgBHQGsHo6WMAAEbAiDAAAF0ShAEA6JIgTFd21o6sJRkA+qNHmK7srB050ZIMAL0xIgwAQJcEYRauQX0Mc9i/YBY3AOiL1ggWrkF9DHPYv2AWNwDoixFhAAC6JAizuOhfAACGRGsEi4v+BQBgSIwIw24YhAaApcmIMOyGQWgAWJqMCAMA0CVBGACALgnCLA3zvPjGzh6pdxgAFg89wiwN87z4xs4eOQ+PBQCGxIgwAABdEoQBAOiSIMzSNaImXr3DALA46BFm6RpRE6/eYQBYHIwI0x9DtgBAjAjTI0O2AECMCAMA0ClBGObJzjoydrbp1ACAuaU1Ap7wRFKdeWxnrRR7aE9vo1MDAOaWEWF4wuRk0tqTt8RwLQAsUUaEYVdG+MW6QQPUTxwf0iA1AHRNEIYFyuQWADC3tEbA3hj0zTftEgCwqBgRhr0xaLjWUC0ALCpGhGFY5mnFOgvjAcBwCMIwLINmndjZzBP7kFr35DECMgDsnCAMc22epmUTkAFgzwjCMAo7S62bNs3boxKr2gHQN0EYFpJ5bADeWUA2igxALwRhWEjmqc94b8uYgwFrABgZ06fBYrBApmuz2h0AS4kRYVisRjCPmi/kAbCUCMKwWO3Nt+DmKKkOoxTBGYD5JgjDUrOrb8HNc+/xnpSiBxmA+TarIFxVZ1XV3VV1T1VdNuB8VdU7ps9vrKoTh18qMHR7MsfxPAzj7qzbwwgyAHNht0G4qvZP8q4k65Icl+SVVXXcjMvWJTlmerswyXuGXCcwX/Z0GDfZs/S6izQ7348WtAH6NpsR4ZOT3NNau7e19miSDyY5d8Y15yb5r23K3yU5pKqePeRagYVoT9PrENPs5KZKS6WNjQ++/dj41Pntt51dO8dBW6gGWHhmE4SPSPKt7fY3Tx/b02sAdrSvQXpXiTXZp3S7LWgPYZvcNIdD17P5tdT4vD9W+AcWutnMI1wDjrW9uCZVdWGmWieS5AdVdfcsnj8XDkvVAyN6NvPjsCTe8dK36/e8adO/BeLubcrgv6rn8InD++3353np8477MMr3PDbo4GyC8OYkR263vzLJfXtxTVprVye5ehbPnFNVtaG1NjHqOpg73nEfvOc+eM9Ln3fch4X4nmfTGvHlJMdU1dFV9ZQk5yW5acY1NyV51fTsEacmeai1dv+QawUAgKHZ7Yhwa21rVV2S5NNJ9k9yTWvtzqq6aPr8VUluTvKSJPck+VGS18xdyQAAsO9m0xqR1trNmQq72x+7arufW5KLh1vanBp5ewZzzjvug/fcB+956fOO+7Dg3nO1tsN32gAAYMmzxDIAAF3qKgjvbqloFr+qOrKqbqmqu6rqzqq6dNQ1MTeqav+q+kpVfXzUtTA3quqQqvpwVf3j9J/pnxl1TQxfVf3m9N/Xd1TV9VV14KhrYt9V1TVV9d2qumO7Yz9RVX9dVV+f/u+/G2WNSUdBeJZLRbP4bU3y2621Y5OcmuRi73nJujTJXaMugjn1fyX5VGvteUnWxPtecqrqiCRvTDLRWluVqS/lnzfaqhiSa5OcNePYZUn+W2vtmCT/bXp/pLoJwpndUtEscq21+1tr/3365+9n6n84rXK4xOJYBkgAAAIuSURBVFTVyiQvTfK+UdfC3KiqZyb5X5L830nSWnu0tfa90VbFHFmW5KlVtSzJ0zJgHQIWn9barUn+ecbhc5P8l+mf/0uS/zivRQ3QUxC2DHRnqmo8yQlJvjjaSpgDb0/ye0l+POpCmDM/lWRLkv883QLzvqp6+qiLYrhaa99OcnmSbya5P1PrEHxmtFUxh/79E+tMTP/3WSOup6sgPKtloFkaquqgJB9J8huttX8ZdT0MT1WdneS7rbXbRl0Lc2pZkhOTvKe1dkKSH2YB/DMqwzXdI3pukqOTHJ7k6VX1v422KnrSUxCe1TLQLH5VtTxTIfi61toNo66HoXtBknOqajJTLU4/V1V/MdqSmAObk2xurT3xLzofzlQwZml5cZJvtNa2tNYeS3JDkp8dcU3Mne9U1bOTZPq/3x1xPV0F4dksFc0iV1WVqZ7Cu1pr/+eo62H4Wmv/e2ttZWttPFN/jj/XWjOCtMS01v6/JN+qqudOHzo9yT+MsCTmxjeTnFpVT5v++/v0+FLkUnZTkldP//zqJB8bYS1JZrmy3FKws6WiR1wWw/eCJBck+fuqun362O9Pr44ILC5vSHLd9ODFvUleM+J6GLLW2her6sNJ/numZv35Shbg6mPsuaq6PslpSQ6rqs1J3pLkbUn+n6p6Xab+T9Avjq7CKVaWAwCgSz21RgAAwDaCMAAAXRKEAQDokiAMAECXBGEAALokCAMA0CVBGACALgnCAAB06f8HJVFLsmpMZ/EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9b7f26e77df64f16a5041d285c6dfb0b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=100000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsQAAAHSCAYAAADrMt2YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de7BlVX0v+u+Ph+JRxCM0CdCwFynbB+HpaUCiJ0ERAUPgHsskIFdEJURLlLzDzcPklmWKU3IvlmIkraHQBCHXaFqCiJIKXpLyeAPElvCIQkhv3UrJQwO+CDaM+0c/3L167b3X3nvt5/x8qmbVXnPOteboyejeX8Yac/yqtRYAAOiq3Za6AQAAsJQEYgAAOk0gBgCg0wRiAAA6TSAGAKDTBGIAADptj6W68H777dd6vd5SXR4AgI64/fbbH26trZnq+JIF4l6vl9tuu22pLg8AQEdU1fh0x02ZAACg0wRiAAA6TSAGAKDTlmwOMQDAXP3oRz/KxMREHn/88aVuCsvIXnvtlbVr12bPPfec1fsEYgBgxZmYmMjee++dXq+Xqlrq5rAMtNbyyCOPZGJiIoceeuis3mvKBACw4jz++OPZd999hWF2qKrsu+++c/rWQCAGAFYkYZh+c+0TAjEAwBxUVV7/+tfveL1ly5asWbMmp59+epLkqquuyoUXXrjL+3q9Xo444ogcffTROfroo/OOd7xj2uv88R//cS699NKRtv1Nb3pT9t9//xx++OE77f/2t7+dk08+OevWrcvJJ5+c73znOyO97nIlEAMAK16vl1SNbhummO4zn/nM3HnnnfnhD3+YJLnpppty0EEHDdXem2++OZs2bcqmTZvyvve9b+5/8Dk677zzcuONN+6y/5JLLslJJ52Ue++9NyeddFIuueSSRW/bUhCIAYAVb3w8aW102/i0dc1+7LTTTsunP/3pJMk111yTs88+ewH/lMmHPvShnHbaaTtC+Fz97M/+bJ773Ofusv9Tn/pU3vCGNyRJ3vCGN2Tjxo3zus5KIRADAMzRWWedlWuvvTaPP/547rjjjhx//PFDve/lL3/5jikTl1122VDvufzyy/O3f/u32bhxY57xjGfsdOzqq6/e8XmTt9e+9rWz+vN861vfygEHHJAkOeCAA/Lggw/O6v0rlWXXAADm6Mgjj8zmzZtzzTXX5NWvfvXQ77v55puz3377DX3+X/zFX2Tt2rXZuHHjwDV2zznnnJxzzjlDfx47E4gBAObhjDPOyG/91m/l85//fB555JEFucbhhx+eTZs2TbnG7tVXX533vOc9u+x/3vOel7/+678e+jo/8RM/kQceeCAHHHBAHnjggey///7zavdKIRADAMzDm970puyzzz454ogj8vnPf37On3P55ZcnycCVKY455pi89a1vzRlnnJHPfvazOfDAA3c6PqoR4jPOOCMf+chHcvHFF+cjH/lIzjzzzHl/5kpgDjEAwDysXbs2F1100cBjV111VdauXbtjm5iYSLLzHOJzzz03SfKv//qv2Xfffae8zste9rJceuml+fmf//k8/PDD82rz2WefnRNOOCFf+cpXsnbt2vz5n/95kuTiiy/OTTfdlHXr1uWmm27KxRdfPK/rrBTVWluSC69fv77ddtttS3JtAGBlu+eee/KiF71ox+teb/iVIYYxNpZs3jy6zxvG6aefnk9+8pN52tOetrgXXmX6+0aSVNXtrbX1U73HlAkAYMVb7PC6EK6//vqlbkJnmTIBAECnCcQAAHSaQMzw+utiDlPXEgBgmTOHmOFtr4u5XdXStQUAYESMEAMA0GkCMQDAHD3rWc9a0ut/61vfyumnn56jjjoqhx122I7y0Zs3b87HPvaxGd//zW9+M6997WvnfP0TTzwxo1xG95FHHsnLX/7yPOtZz9qlQMntt9+eI444Is973vPyjne8I6NcOlggBgBWvv7nXOa7rZDnZN75znfm5JNPzpe//OXcfffdueSSS5IMH4gPPPDAWZV2Xmh77bVX3vWud+XSSy/d5dhb3/rWbNiwIffee2/uvffe3HjjjSO7rkAMAKx8259zGdU2iyof3/ve93LSSSflxS9+cY444oh86lOfSrI1lL7whS/M+eefn8MPPzznnHNO/u7v/i4vfelLs27duvzTP/1TkuT73/9+3vSmN+XYY4/NMcccs+P9d911V4477rgcffTROfLII3Pvvffucu0HHngga9eu3fH6yCOPTLK14tw//MM/5Oijj85ll12WJ598Mr/927+dY489NkceeWT+7M/+bEcbDz/88CRbq+q95jWvyamnnpp169bld37nd4a+B0899VTe8IY35A/+4A+Gfs8gz3zmM/Oyl70se+211y5/zsceeywnnHBCqirnnntuNm7cOK9rTeahOgCAedhrr73yN3/zN3n2s5+dhx9+OC95yUtyxhlnJEnuu+++fPzjH8+GDRty7LHH5mMf+1j+8R//Mdddd13+5E/+JBs3bsy73/3uvOIVr8iVV16Z//iP/8hxxx2XV77ylbniiity0UUX5ZxzzskTTzyRJ598cpdrv+1tb8sv//Iv5/LLL88rX/nKvPGNb8yBBx6YSy65JJdeeumOYh8bNmzIPvvsk1tvvTX/+Z//mZe+9KV51atelep7QH7Tpk350pe+lKc//el5wQtekLe//e05+OCDp/3zb9myJeecc04OP/zw/P7v//4ux3/91389N9988y77zzrrrKFLQ3/jG9/YKfivXbs23/jGN4Z67zAEYgCAeWit5fd+7/dyyy23ZLfddss3vvGNfOtb30qSHHrooTniiCOSJD/90z+dk046KVWVI444Ipu3ldf73Oc+l+uuu27HNIHHH388X/va13LCCSfk3e9+dyYmJvKa17wm69at2+Xap5xySu6///7ceOON+cxnPpNjjjkmd9555y7nfe5zn8sdd9yxY3rEo48+mnvvvTfPf/7zdzrvpJNOyj777JMkOeywwzI+Pj5jIP7VX/3V/NIv/dLAMJwkl1122bTvH8ag+cL9YX4+BGIAgHm4+uqr89BDD+X222/PnnvumV6vl8cffzxJ8vSnP33HebvtttuO17vttlu2bNmSZGvY+8QnPpEXvOAFO33ui170ohx//PH59Kc/nVNOOSUf/vCHc8899+RDH/pQkuSGG27IgQcemOc+97l53etel9e97nU5/fTTc8stt2Tffffd6bNaa3n/+9+fU045Zaf9m/tqXk9u7+67776jjdP5mZ/5mdx88835zd/8zV2mOiSjGSFeu3ZtJiYmdryemJjIgQceONR7h2EOcVcpsgEAI/Hoo49m//33z5577pmbb74547OYf5xsHeV9//vfv2MU9Etf+lKS5P77789P/dRP5R3veEfOOOOM3HHHHXnb296WTZs2ZdOmTTnwwAPz93//9/nBD36QJPnud7+bf/u3f8shhxySvffeO9/97nd3usYHP/jB/OhHP0qSfPWrX833v//9odt47rnn7pjz3O/Nb35zXv3qV+cXf/EXBwboyy67bEebJ2/DhuEkOeCAA7L33nvni1/8Ylpr+ehHP5ozzzxz6PfPZMYR4qq6MsnpSR5srR0+xTknJnlvkj2TPNxa+7mRtZCFocgGAMzLli1b8vSnPz3nnHNOfuEXfiHr16/P0UcfnRe+8IWz+pw//MM/zK/92q/lyCOPTGstvV4v119/ff7qr/4qf/mXf5k999wzP/mTP5l3vvOdu7z39ttvz4UXXpg99tgjTz31VM4///wce+yx+dGPfpQ99tgjRx11VM4777xcdNFF2bx5c1784hentZY1a9bM6qG0O+64IwcccMCUx3/jN34jjz76aF7/+tfn6quvzm67zX3Mtdfr5bHHHssTTzyRjRs35nOf+1wOO+ywfPCDH8x5552XH/7whznttNNy2mmnzfka/WqmNdyq6meTfC/JRwcF4qp6TpIvJDm1tfa1qtq/tfbgTBdev359G+W6dcxS1a6BuL8v9Ho7P2U7NpZM/mpl0HsAYBHcc889edGLXvTjHf2/s+ar/3feAF/+8pfzK7/yK1OOnK4Wjz32WN785jfn4x//+FI3ZSi79I0kVXV7a239VO+ZcYS4tXZLVfWmOeV1ST7ZWvvatvNnDMOsEP2jyACwXM0QXkftiiuuyPve9768973vXdTrLoVnP/vZKyYMz9Uo5hA/P8l/rarPV9XtVXXuCD4TAGDZestb3pK77747r3rVq5a6KYzAKFaZ2CPJf0tyUpJnJPlfVfXF1tpX+0+sqguSXJAkhxxyyAguDQAA8zOKEeKJJDe21r7fWns4yS1Jjhp0YmttQ2ttfWtt/Zo1a0ZwaQCgq2Z6DorumWufGEUg/lSS/15Ve1TVf0lyfJJ7RvC5LKaxsV3ruI+NLXWrAGCgvfbaK4888ohQzA6ttTzyyCMD10KeyTDLrl2T5MQk+1XVRJI/ytbl1dJau6K1dk9V3ZjkjiRPJflwa23XEiksb4v8MAIAzMf2Qg0PPfTQUjeFZWSvvfbaqcTzsGZcdm2hWHZtiY1iyTTLrgEAK8BMy66pVMfcDZpmoeIdALDCCMTM3ebNW0eIJ2+JgAwArCijWHYNfqx/LrKS0ADAMmeEGACAThOIAQDoNIEYAIBOE4gBAOg0gRgAgE4TiAEA6DSBGACAThOIAQDoNIEYRq3XU60PAFYQlepg1MbHf1zGOlGtDwCWOSPEAAB0mkDMwhobM30AAFjWTJlgYW3evPNr0wcAgGXGCDEAAJ0mEHdF/8oHY2NL3SIAgGXBlImu6F/5AACAJEaIAQDoOIEYAIBOE4gBAOg0gRgAgE4TiAEA6DSBGACAThOIAQDoNIEYAIBOE4gBAOg0gRgAgE4TiAEA6DSBGACAThOIWVxjY0nVzluvt9StAgA6bI+lbgAds3nzrvuqFr0ZAADbGSFejXq9XUdhx8aWulUAAMuSEeLVaHw8aW2pWwEAsCIYIQYAoNMEYgAAOk0gBgCg0wRiAAA6TSAGAKDTBGIAADptxkBcVVdW1YNVdecM5x1bVU9W1WtH1zwAAFhYw4wQX5Xk1OlOqKrdk/zPJJ8dQZsAAGDRzBiIW2u3JPn2DKe9Pcknkjw4ikYBAMBimfcc4qo6KMn/SHLF/JsDAACLaxQP1b03ye+21p6c6cSquqCqbquq2x566KERXBoAAOZnjxF8xvok11ZVkuyX5NVVtaW1trH/xNbahiQbkmT9+vVtBNcGAIB5mXcgbq0duv3nqroqyfWDwjAAACxHMwbiqromyYlJ9quqiSR/lGTPJGmtmTcMAMCKNmMgbq2dPeyHtdbOm1drAABgkalUBwBApwnEAAB0mkAMAECnCcQAAHSaQAwAQKcJxAAAdJpAzNIbG0uqfrz1ekvdIgCgQ0ZRuhnmZ/PmnV9vLQMOALAojBADANBpAjEAAJ0mEAMA0GkCMQAAnSYQAwDQaQLxatDr7bxs2djYUrcIAGDFsOzaajA+nrS21K0AAFiRjBADANBpAjEAAJ0mEAMA0GkCMcvP2NjODwn2ekvdIgBgFfNQHcvP5s07v65akmYAAN1ghBgAgE4TiAEA6DSBGACAThOIAQDoNIEYAIBOE4gBAOg0gRgAgE4TiAEA6DSBGACAThOIAQDoNIEYAIBOE4gBAOg0gRgAgE4TiAEA6DSBGACAThOIAQDoNIEYAIBOE4gBAOg0gRgAgE4TiFn+xsaSqp23Xm+pWwUArBJ7LHUDYEabN++6r2rRmwEArE4zjhBX1ZVV9WBV3TnF8XOq6o5t2xeq6qjRNxMAABbGMFMmrkpy6jTH/z3Jz7XWjkzyriQbRtAuAABYFDNOmWit3VJVvWmOf2HSyy8mWTv/ZgEAwOIY9UN1b07ymRF/JgAALJiRPVRXVS/P1kD8smnOuSDJBUlyyCGHjOrSAAAwZyMZIa6qI5N8OMmZrbVHpjqvtbahtba+tbZ+zZo1o7g0AADMy7wDcVUdkuSTSV7fWvvq/JsEAACLZ8YpE1V1TZITk+xXVRNJ/ijJnknSWrsiyTuT7JvkT2vr2rBbWmvrF6rBAAAwSsOsMnH2DMfPT3L+yFoEAACLSOlmAAA6TSAGAKDTBGIAADpNIAYAoNMEYgAAOk0gBgCg0wRiAAA6TSAGAKDTBGIAADpNIAYAoNMEYgAAOk0gBgCg0wTilajXS6p+vI2NLXWLAABWrD2WugHMwfh40tpStwIAYFUwQgwAQKcJxAAAdJpAzMo0NrbzPOpeb6lbBACsUOYQszJt3rzz66olaQYAsPIZIQYAoNMEYgAAOk0gBgCg0wRiAAA6TSAGAKDTBGIAADpNIAYAoNMEYgAAOk0gBgCg0wRiAAA6TSAGAKDTBGIAADpNIAYAoNMEYgAAOk0gBgCg0wRiAAA6TSAGAKDTBGIAADpNIAYAoNMEYlaHsbGkauet11vqVgEAK8AeS90AGInNm3fdV7XozQAAVh4jxAAAdJpADABAp80YiKvqyqp6sKrunOJ4VdX7quq+qrqjql48+mYCAMDCGGaE+Kokp05z/LQk67ZtFyT54PybxQ693q4Pi42NLXWrAABWjRkfqmut3VJVvWlOOTPJR1trLckXq+o5VXVAa+2BEbWx28bHk9aWuhUAAKvWKOYQH5Tk65NeT2zbBwAAy94oAvGgta0GDmlW1QVVdVtV3fbQQw+N4NIAADA/owjEE0kOnvR6bZJvDjqxtbahtba+tbZ+zZo1I7g0AADMzygC8XVJzt222sRLkjxq/jDL0aDnExW0AwBmfKiuqq5JcmKS/apqIskfJdkzSVprVyS5Icmrk9yX5AdJ3rhQjYX5mOr5RAXtAKDbhlll4uwZjrckbxtZiwAAYBGpVAcAQKcJxAAAdJpADABApwnErDrbV5NIVLwGAGYmELPszXa5tMmrSbT2423z5kVqMACwogjELHvbA27/lgwOykaCAYDZmHHZNViujPgCAKNghBgAgE4TiAEA6DSBmNVrbGy4p/CmMNuH+QCAlckcYlav/knG29diG9Lk1Srm8TEAwDInENN52weSB+0HAFY/gZjOs1oFAHSbOcQAAHSaQAwAQKcJxAAAdJpADABApwnEMEv9yxtbnxgAVjaBmEW30gtebN68dX3i/m18fKlbBgDMhWXXWHQKXgAAy4lAzLKhQAYAsBQEYpYNBTIAgKVgDjEAAJ0mEAMA0GkCMQAAnWYOMYzI9ocCW3Z+OHDAghoAwDJihHi56V+k1xILK8b29YmTndcnBgCWNyPEy81Ui/QCALAgjBADANBpAjEsgpVcqhoAVjtTJlgwvd7WGSD9ujgtWqlqAFi+BGIWzLKbDt1fG3psTHk8AEAgpkP6w68hWgAg5hADANBxAjEAAJ0mEAMA0GkCMQAAnSYQAwDQaQIxAACdJhADANBpAjEAAJ02VCCuqlOr6itVdV9VXTzg+D5V9bdV9eWququq3jj6pgIAwOjNGIiravckH0hyWpLDkpxdVYf1nfa2JHe31o5KcmKS/6uqnjbitsKqsr2SdP/W6y11ywCgW4Yp3Xxckvtaa/cnSVVdm+TMJHdPOqcl2buqKsmzknw7yZYRtxVWlf5K0tupKA0Ai2uYKRMHJfn6pNcT2/ZNdnmSFyX5ZpJ/SXJRa+2pkbQQAAAW0DCBeNB4Vet7fUqSTUkOTHJ0ksur6tm7fFDVBVV1W1Xd9tBDD826sQAAMGrDBOKJJAdPer02W0eCJ3tjkk+2re5L8u9JXtj/Qa21Da219a219WvWrJlrmwEAYGSGCcS3JllXVYdue1DurCTX9Z3ztSQnJUlV/USSFyS5f5QNBQCAhTDjQ3WttS1VdWGSzybZPcmVrbW7quot245fkeRdSa6qqn/J1ikWv9tae3gB2w0AACMxzCoTaa3dkOSGvn1XTPr5m0leNdqmQTdtX45t0P6pVqYAAOZuqEAM0+n1kvHxXfePjS16U1YFy7EBwOISiJm38fGk9a87AgCwQgxVuhkAAFYrgRgAgE4TiAEA6DSBGACAThOIAQDoNIEYAIBOE4gBAOg0gZju2l4SbvLW6y11qwCARaYwB901qCSccnAA0DlGiBlar7frgGqVEs0AwMpmhJihKdG8tLbP8Bi0f9BgNwAwHIEYVoipQq9ZHgAwP6ZMAADQaQIxAACdJhADANBpAjEAAJ0mEAMA0GkCMQAAnSYQw2T95ZyVcgaAVc86xEup19ta7WIyZd+WVv9ivxb5BYBVTyBeSkq/AQAsOVMmAADoNIEYAIBOE4gBAOg0gRhWuP6FMSyQAQCz46E6WOH6F8bYzgIZADAcI8QAAHSaQAwAQKcJxAAAdJpADABApwnE7KLXG7xqgarSAMBqZJUJdqGiNADQJUaIYZWyPjEADMcIMaxS1icGgOEYIQYAoNMEYgAAOk0gBgCg0wRi6BgP2wHAzjxUBx3jYTsA2NlQI8RVdWpVfaWq7quqi6c458Sq2lRVd1XV/zvaZgIAwMKYcYS4qnZP8oEkJyeZSHJrVV3XWrt70jnPSfKnSU5trX2tqvZfqAYDAMAoDTNCfFyS+1pr97fWnkhybZIz+855XZJPtta+liSttQdH20wAAFgYwwTig5J8fdLriW37Jnt+kv9aVZ+vqtur6txRNRAAABbSMA/VDXrUpg34nP+W5KQkz0jyv6rqi621r+70QVUXJLkgSQ455JDZtxYAAEZsmBHiiSQHT3q9Nsk3B5xzY2vt+621h5PckuSo/g9qrW1ora1vra1fs2bNXNsMAAAjM0wgvjXJuqo6tKqeluSsJNf1nfOpJP+9qvaoqv+S5Pgk94y2qbAE+hfttVgvAKw6M06ZaK1tqaoLk3w2ye5Jrmyt3VVVb9l2/IrW2j1VdWOSO5I8leTDrbU7F7LhsCj6F+21WC8ArDpDFeZord2Q5Ia+fVf0vX5PkveMrmkstF4vGR/fdf/Y2KI3hWVg+2D4oP1TFfMAgNVApboOGx9PWv/jkXSWCnYAdNVQleoAAGC1EogBAOg0gRgAgE4TiAEA6DSBGACAThOIF1Ovt3ORB+ubAQAsOcuuLSbrnAEALDtGiIFp9VevVsUagNXGCDEwLQU7AFjtjBADANBpAjEAAJ0mEMNsDJpQazItAKxo5hDDbAyaUGsyLQCsaEaIAQDoNIEYAIBOE4gBAOg0gRgAgE4TiAEA6DSBGACAThOIAQDoNIEYmJNBNUrUKQFgJVKYA5iTQTVKEnVKAFh5jBADANBpAnEH9HqDv9oeG1vqlgEALD1TJjpgfDxpbalbAQCwPBkhBgCg0wRiAAA6TSAGAKDTBGJgpKxPDMBK46E6mK/tCXDy6w6zPjEAK41ADPM1VQIEAFYEUyYAAOg0gRgAgE4TiAEA6DSBGACAThOIAQDoNIEYWBTWJwZgubLsGrAorE8MwHJlhHgh9Xo7D4V1vGADAMByZIR4IY2PJ60tdSsAAJjGUCPEVXVqVX2lqu6rqounOe/Yqnqyql47uiYCAMDCmTEQV9XuST6Q5LQkhyU5u6oOm+K8/5nks6NuJAAALJRhRoiPS3Jfa+3+1toTSa5NcuaA896e5BNJHhxh+wAAYEENE4gPSvL1Sa8ntu3boaoOSvI/klwxuqYBAMDCGyYQD1oUqf9Jsfcm+d3W2pPTflDVBVV1W1Xd9tBDDw3bRgAAWDDDBOKJJAdPer02yTf7zlmf5Nqq2pzktUn+tKr+t/4Paq1taK2tb62tX7NmzRybzFT6V3mz2hsAwMyGWXbt1iTrqurQJN9IclaS100+obV26Pafq+qqJNe31jaOsJ0MwSpvrETbK9hNdWyqgh4AMCozBuLW2paqujBbV4/YPcmVrbW7quot246bNwzM2XSBVxU7ABbDUIU5Wms3JLmhb9/AINxaO2/+zQIAgMWhdDMAAJ0mEAMA0GkCMbBsbX/grn/r9Za6ZQCsJkPNIQZYClM9cOdhOwBGyQgxAACdJhADANBpAjEAAJ0mEAMA0GkCMQAAnSYQAwDQaQIxAACdJhADANBpAjGw4qhgB8AoqVQHrDgq2AEwSkaIAQDoNIF4Ber1Bn9dPDa21C0DAFh5TJlYgcbHk9aWuhUAAKuDEWIAADpNIAYAoNMEYgAAOk0gBgCg0wTiURm09INlHwAAlj2rTIyKpR8AAFYkI8TAqqGkMwBzYYQYWDWUdAZgLowQA6uekWMApmOEGFj1jBwDMB0jxAAAdJpADABApwnEAAB0mkAMAECnCcQAAHSaQAwAQKcJxMtYrzd47dSxsaVuGbMyaBFcC+ACwLJhHeJlbHw8aW2pW8G8DVoE1wK4ALBsGCEGOksFOwASI8RAh6lgB0BihBgAgI4TiAEA6DSBGACAThOIAQDotKECcVWdWlVfqar7quriAcfPqao7tm1fqKqjRt9UgMVh9QmAbplxlYmq2j3JB5KcnGQiya1VdV1r7e5Jp/17kp9rrX2nqk5LsiHJ8QvRYICFZvUJgG4ZZoT4uCT3tdbub609keTaJGdOPqG19oXW2ne2vfxikrWjbSasMv1DkIYeAWDJDLMO8UFJvj7p9USmH/19c5LPzKdRsOr1D0EaegSAJTNMIB70m3pgQeGqenm2BuKXTXH8giQXJMkhhxwyZBMBAGDhDDNlYiLJwZNer03yzf6TqurIJB9OcmZr7ZFBH9Ra29BaW99aW79mzZq5tBcAAEZqmEB8a5J1VXVoVT0tyVlJrpt8QlUdkuSTSV7fWvvq6JsJAAALY8YpE621LVV1YZLPJtk9yZWttbuq6i3bjl+R5J1J9k3yp7V1LuSW1tr6hWs2AACMxjBziNNauyHJDX37rpj08/lJzh9t0wCWl+2Lg8zm/KmWcANg+RgqEAMw+3Br8RCAlUHpZgAAOk0gBlggSkADrAymTMBy0D851eTTVUEJaICVQSCG5UDlOgBYMqZMzFWvt/N3oGNjS90iAADmQCCeq/HxpLUfb/P4ers/W8vYDJx8auIpACwIUyaWge3ZGnYY9D9YplGsGlOtZ2zqOMDSEIgBFpmH7QCWF1MmAADoNIEYAIBOE4hhpeh/0M5DdgAwEuYQw0phrWIAWBBGiAEA6DSBGGCZGLT8tNkxAAvPlAmAZcJybJTa79QAAAdsSURBVABLwwjxIlKRDgBg+RGIF1F/tecRVH2my6w60RmmUgAsLFMmYKWy6kRnmEoBsLCMEMNqMWgY0RDiqmbkGGA0jBDDajFoGNEQ4qpm5BhgNIwQAwDQaQIxwCpjKgXA7JgyAbDKmEoBMDtGiGE1szQbAMxIIB5Wf1UN1TRYCTZv3nnR68T36ADQx5SJYW2vqgErmZUoOm37FwaD9isQBHSZQAxd15+SpKNVy9xigMEEYug6Fe86z8gx0HXmEI9Y/1Rj045ZcTyI1zn9U82nm3KuSwCrkRHiETPVmBXPiDHbmGIBdIURYmB6g6o8GCLsNIU/gNVGIAamN+j79EQS6rCppliMjy91ywDmRiAGZm+m9Y0F5E4ycgysVAIxMH8KgJDZP5w31aarAItNIAZGb1Ay8n16Z00VlKfaEkEZWFwC8SCD1k6zZhrMj+XcGJI5ysBiE4gH2b522uStb/2hqdYblpthCv0pZ5h00/8XTYjuNHOUgYUiEM/RoMw8IDcDUxlmxLj/L5ohwk5TQARYKAIxsDSGeRCv/yuXqYYIpaBO8zAfMF8q1QHLwzBfrwxzzkxl1Hq9XUeax8Z8vbMKzfY/qQp80F1DjRBX1alV9ZWquq+qLh5wvKrqfduO31FVLx59U5eGucKwwsw0ipx4WouBZjtHearfD0agYeWZMRBX1e5JPpDktCSHJTm7qg7rO+20JOu2bRck+eCI27mw+v9Vm5R2zRWGFWamNb4G/eUdZiqGdLPqzXbqRTKa5eR0MVh6w4wQH5fkvtba/a21J5Jcm+TMvnPOTPLRttUXkzynqg4YcVsXTn/qlXahW2a7UO5c0o2Es2JN1T1m+6tiVOsxz2XT9WB6wwTig5J8fdLriW37ZnvO8jFgRNjUCGBWFitEL9S2UAmp/x/TYa4zl/d0wFy62Eroesux20+VAVZKV1zp7V8Ohnmorgbsa3M4J1V1QbZOqUiS71XVV4a4/ijsl+ThKY+Oj2fwH2HroRp8qEumv39Mx72bH/dv7mb+d28x/nGby3WWxz+8+t7cLdt7N9uutURdcWT3b3n8VVpU0927aYc4hwnEE0kOnvR6bZJvzuGctNY2JNkwxDVHqqpua62tX+zrrhbu39y5d/Pj/s2dezc/7t/cuXfz4/7N3Xzu3TBTJm5Nsq6qDq2qpyU5K8l1fedcl+TcbatNvCTJo621B+bSIAAAWEwzjhC31rZU1YVJPptk9yRXttbuqqq3bDt+RZIbkrw6yX1JfpDkjQvXZAAAGJ2hCnO01m7I1tA7ed8Vk35uSd422qaN1KJP01hl3L+5c+/mx/2bO/dufty/uXPv5sf9m7s537tqbZdn3wAAoDOGqlQHAACr1aoKxF0uMT1fQ9y7E6vq0aratG1751K0czmqqiur6sGqunOK4/rdNIa4f/reFKrq4Kq6uaruqaq7quqiAefof1MY8v7pfwNU1V5V9U9V9eVt9+7/HHCOvjeFIe+fvjeNqtq9qr5UVdcPODbrvjfUHOKVYFKJ6ZOzdRm4W6vqutba3ZNOm1xi+vhsLTF9/GK3dbkZ8t4lyT+01k5f9AYuf1cluTzJR6c4rt9N76pMf/8SfW8qW5L8Zmvtn6tq7yS3V9VN/t0b2jD3L9H/BvnPJK9orX2vqvZM8o9V9Zlt1Wq30/emNsz9S/S96VyU5J4kzx5wbNZ9bzWNEK/+EtMLZ5h7xxRaa7ck+fY0p+h30xji/jGF1toDrbV/3vbzd7P1l0N/lVD9bwpD3j8G2Nafvrft5Z7btv6HkvS9KQx5/5hCVa1N8vNJPjzFKbPue6spEK++EtOLZ9j7csK2r3c+U1U/vThNWxX0u/nT92ZQVb0kxyT5//oO6X9DmOb+JfrfQNu+st6U5MEkN7XW9L1ZGOL+JfreVN6b5HeSPDXF8Vn3vdUUiEdWYrqDhrkv/5xkrLV2VJL3J9m44K1aPfS7+dH3ZlBVz0ryiSS/1lp7rP/wgLfof5PMcP/0vym01p5srR2drdVpj6uqw/tO0femMcT90/cGqKrTkzzYWrt9utMG7Ju2762mQDyyEtMdNON9aa09tv3rnW3rUu9ZVfstXhNXNP1uHvS96W2bf/iJJFe31j454BT9bxoz3T/9b2attf9I8vkkp/Yd0veGMNX90/em9NIkZ1TV5myd4vmKqvrLvnNm3fdWUyBWYnruZrx3VfWTVVXbfj4uW/vOI4ve0pVJv5sHfW9q2+7Lnye5p7X2f09xmv43hWHun/43WFWtqarnbPv5GUlemeRf+07T96YwzP3T9wZrrf0frbW1rbVetuaVv2+t/e99p826762aVSaUmJ67Ie/da5O8taq2JPlhkrOaqi5Jkqq6JsmJSfarqokkf5StD0jod0MY4v7pe1N7aZLXJ/mXbXMRk+T3khyS6H9DGOb+6X+DHZDkI9tWKdotyf/TWrve79yhDXP/9L1ZmG/fU6kOAIBOW01TJgAAYNYEYgAAOk0gBgCg0wRiAAA6TSAGAKDTBGIAADpNIAYAoNMEYgAAOu3/ByRxdZbMMM17AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e2b8a55a559b4cffb87b0ae60bdfec33",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=100000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArkAAAHSCAYAAADohdOwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5TddX03+vcXEomPIi0QLDF2Nj2hCBIuPgFK8fQoSLmU4jouqmgeLqKltiLU3vT0wulZLlvOKufgUfqo8XKgLUqPt2C9UGiLD3U9tTZ5jDxctKE4owEeCLEFpXIJfM8fkwmTZJLZk9l79uzvvF5r7ZWZvffMfPJjM793Pvv7+3xLrTUAANCSfQZdAAAA9JqQCwBAc4RcAACaI+QCANAcIRcAgOYIuQAANGdRP77pwQcfXDudTj++NQAAJEnWr1//SK116VSP9SXkdjqdrFu3rh/fGgAAkiSllLHdPWa5AgAAzRFyAQBojpALAEBz+rImFwBgT55++uls2rQpTzzxxKBLYQgsWbIky5cvz+LFi7v+GiEXAJhzmzZtyv77759Op5NSyqDLYR6rtWbLli3ZtGlTDjvssK6/znIFAGDOPfHEEznooIMEXKZVSslBBx00466/kAsADISAS7f25rUi5AIAC1IpJRdccMH2z7du3ZqlS5fmnHPOSZJcd911ueyyy3b5uk6nk5UrV+a4447Lcccdl8svv3yPP+cP//APc/XVV/e09ksuuSSHHHJIjj766B3u//73v5/TTz89hx9+eE4//fT867/+6/bH/viP/zgrVqzIEUcckb/+67/uaT3zkZALAAxcp5OU0rtbNxuvvuAFL8idd96ZH/3oR0mSW2+9NS95yUu6qve2227Lhg0bsmHDhrz//e/f+7/4Xrr44otz880373L/VVddldNOOy0bN27MaaedlquuuipJcvfdd+fGG2/MXXfdlZtvvjm/9mu/lmeeeWauy55TQi4AMHBjY0mtvbuN7XYfrB2dddZZ+eIXv5gk+eQnP5k3vvGNffxbJh/5yEdy1llnbQ/We+vnfu7ncuCBB+5y/0033ZSLLrooSXLRRRdl7dq12+8///zzs99+++Wwww7LihUr8vWvf31WNcx3Qi4AsGCdf/75ufHGG/PEE0/kjjvuyEknndTV17361a/evlzhmmuu6eprrr322vzVX/1V1q5dm+c///k7PHbDDTds/36Tb+edd96M/j4PPfRQDj300CTJoYcemocffjhJcv/99+elL33p9uctX748999//4y+97AxQgwAWLCOOeaYjI6O5pOf/GTOPvvsrr/utttuy8EHH9z18//8z/88y5cvz9q1a6ec9bp69eqsXr266+83U7XWXe5r/cI/IRcAWNDOPffc/NZv/Va+8pWvZMuWLX35GUcffXQ2bNiw21mvN9xwQ/7kT/5kl/tXrFiRT3/6013/nBe/+MV58MEHc+ihh+bBBx/MIYcckmS8c/u9731v+/M2bdqUZcuW7cXfZHhYrgAALGiXXHJJrrzyyqxcuXJW3+faa6/NtddeO+Vjxx9/fD784Q/n3HPPzQMPPLDL46tXr95+Idvk20wCbjIe2K+//vokyfXXX5/Xvva12++/8cYb8+STT+Y73/lONm7cmBNPPHGGf8PhIuQCAAva8uXLc8UVV0z52HXXXZfly5dvv23atCnJjmtyL7zwwiTJt771rRx00EG7/TmvfOUrc/XVV+cXfuEX8sgjj8yq5je+8Y05+eST8+1vfzvLly/Pxz72sSTJu9/97tx66605/PDDc+utt+bd7353kuTlL395Xv/61+eoo47KmWeemT/90z/NvvvuO6sa5rsy1RqN2Vq1alVdt25dz78vANCGe+65J0ceeeT2zzud7icidGNkJBkd7d3368Y555yTz372s3ne8543tz94gdj5NZMkpZT1tdZVUz3fmlwAYODmOpD2wxe+8IVBl8AklisAANAcIRcAgOYIuczexF6M3eyhCAAwB6zJZfYm9mJsfKg0ADA8dHIBAGiOkAsALFgvfOELB/rzH3rooZxzzjk59thjc9RRR23fWnh0dDSf+MQnpv36Bx54IOedd95e//xXvepV6eXY1y1btuTVr351XvjCF+ayyy7b4bH169dn5cqVWbFiRS6//PLtWw0/+eSTecMb3pAVK1bkpJNOymiPRm0IuQDA4E1c39Gr25BcJ3LllVfm9NNPzze/+c3cfffdueqqq5J0H3KXLVs2413R+mnJkiV5z3vek6uvvnqXx371V381a9asycaNG7Nx48bcfPPNSZKPfexj+fEf//Hce++9eec735l3vetdPalFyGVmXGQGQD9MXN/Rq9sMdpb44Q9/mNNOOy2veMUrsnLlytx0001JxoPmy172srz1rW/N0UcfndWrV+dv/uZvcsopp+Twww/P17/+9STJ448/nksuuSQnnHBCjj/++O1ff9ddd+XEE0/Mcccdl2OOOSYbN27c5Wc/+OCDWb58+fbPjznmmCTjO5f9/d//fY477rhcc801eeaZZ/Lbv/3bOeGEE3LMMcfkwx/+8PYajz766CTju7O97nWvy5lnnpnDDz88v/M7v9P1MXj22Wdz0UUX5fd///e7/pqpvOAFL8grX/nKLFmyZJe/52OPPZaTTz45pZRceOGFWbt2bZLkpptuykUXXZQkOe+88/K3f/u36cVmZS48Y2ZcZAZAY5YsWZLPfe5zedGLXpRHHnkkP/MzP5Nzzz03SXLvvffmU5/6VNasWZMTTjghn/jEJ/LVr341n//85/NHf/RHWbt2bd773vfm1FNPzcc//vH827/9W0488cS85jWvyYc+9KFcccUVWb16dZ566qk888wzu/zst7/97XnDG96Qa6+9Nq95zWvy5je/OcuWLctVV12Vq6++evsGE2vWrMkBBxyQf/qnf8qTTz6ZU045JT//8z+fstP5eMOGDfnGN76R/fbbL0cccUTe8Y535KUvfeke//5bt27N6tWrc/TRR+f3fu/3dnn8ne98Z2677bZd7j///PO3bxs8nfvvv3+HML98+fLcf//92x+bqHHRokU54IADsmXLlhx88MFdfe/dEXIBgAWt1prf/d3fze2335599tkn999/fx566KEkyWGHHZaVK1cmSV7+8pfntNNOSyklK1eu3L529JZbbsnnP//57W/RP/HEE/nud7+bk08+Oe9973uzadOmvO51r8vhhx++y88+44wzct999+Xmm2/Ol7/85Rx//PG58847d3neLbfckjvuuGP70oRHH300GzduzE//9E/v8LzTTjstBxxwQJLkqKOOytjY2LQh91d+5Vfy+te/fsqAmyTXXHPNHr++G1N1ZicC+p4emw0hFwBY0G644YZs3rw569evz+LFi9PpdPLEE08kSfbbb7/tz9tnn322f77PPvtk69atScZD2mc+85kcccQRO3zfI488MieddFK++MUv5owzzshHP/rR3HPPPfnIRz6SJPnSl76UZcuW5cADD8yb3vSmvOlNb8o555yT22+/PQcddNAO36vWmg984AM544wzdrh/54u0Jte77777bq9xT372Z382t912W37zN39zl2UGSW86ucuXL8+mTZu2f75p06YsW7Zs+2Pf+973snz58mzdujWPPvpoDjzwwK6+755Yk0v/WccLwDz26KOP5pBDDsnixYtz2223ZWwG63mT8W7sBz7wge0dyW984xtJkvvuuy8/9VM/lcsvvzznnntu7rjjjrz97W/Phg0bsmHDhixbtix/93d/l3//939PkvzgBz/Iv/zLv+Qnf/Ins//+++cHP/jBDj/jgx/8YJ5++ukkyT//8z/n8ccf77rGCy+8cPsa4p295S1vydlnn51f+qVfmjIUX3PNNdtrnnzrNuAmyaGHHpr9998/X/va11JrzZ/92Z/lta99bZLk3HPPzfXXX58k+fSnP51TTz1VJ5d5rNN5btH/yIh1vADMO1u3bs1+++2X1atX5xd/8RezatWqHHfccXnZy142o+/zB3/wB/n1X//1HHPMMam1ptPp5Atf+EL+8i//Mn/xF3+RxYsX5yd+4idy5ZVX7vK169evz2WXXZZFixbl2WefzVvf+taccMIJefrpp7No0aIce+yxufjii3PFFVdkdHQ0r3jFK1JrzdKlS7dfuNWNO+64I4ceeuhuH/+N3/iNPProo7ngggtyww03ZJ999r4P2ul08thjj+Wpp57K2rVrc8stt+Soo47KBz/4wVx88cX50Y9+lLPOOitnnXVWkvGQfcEFF2TFihU58MADc+ONN+71z56s9OLqtZ2tWrWq9nLmGvNIKc8F1onXzp7um+prAVjw7rnnnhx55JHP3TG5OdILIyPJNPNWv/nNb+aXf/mXd9vhbMVjjz2Wt7zlLfnUpz416FJmZZfXTJJSyvpa66qpnq+TCwAMXo82AOjWhz70obz//e/P+973vjn9uYPwohe9aOgD7t4QcgGABedtb3tb3va2tw26DPrIhWcAADRHyAUABqIf1wXRpr15rQi5AMCcW7JkSbZs2SLoMq1aa7Zs2TLlDN89sSYXAJhzE5sDbN68edClMASWLFmyw7bA3egq5JZSfizJR5McnaQmuaTW+g8zrhAAIMnixYtz2GGHDboMGtZtJ/f/SXJzrfW8UsrzkvyHPtYEAACzMm3ILaW8KMnPJbk4SWqtTyV5qr9lMe+NjDy3g9nIyGBrAQDYSTed3J9KsjnJ/1tKOTbJ+iRX1Fq73zCZ9szx0G4AgJnoZrrCoiSvSPLBWuvxSR5P8u6dn1RKubSUsq6Uss4icgAABqmbkLspyaZa6z9u+/zTGQ+9O6i1rqm1rqq1rlq6dGkva2RYTCxhKMUSBgBgoKZdrlBr/R+llO+VUo6otX47yWlJ7u5/aQwdSxgAgHmi2+kK70hyw7bJCvcleXP/SgIAgNnpKuTWWjckWdXnWgAAoCds6wsAQHOEXAAAmiPkMncmT1/odAZdDQDQsG4vPIPZmzx9YWK3NACAPtDJZXqdjvm3AMBQ0cllemNjSa2DrgIAoGs6uQAANEfIBQCgOUIuAADNEXIBAGiOkAsAQHOEXAAAmiPkAgDQHCEXAIDmCLkAADRHyAUAoDlCLgAAzRFyAQBojpDLYIyMJKUknc6gKwEAGrRo0AWwQI2Ojv9ZykDLAADapJPL1Dqd8QBaynjXFQBgiOjkMrWxsaTWQVcBALBXdHIBAGiOkAsAQHOEXAAAmiPkAgDQHCEXAIDmCLkAADRHyAUAoDlCLgAAzRFyAQBojpALAEBzhFwGa2QkKWX81ukMuhoAoBGLBl0AC9zo6HMflzKwMgCAtujkAgDQHCEXAIDmCLkAADRHyGVHnc742tiRkUFXAgCw11x4xo7GxpJaB10FAMCs6OQCANAcIRcAgOYIuQAANEfIBQCgOUIuAADNEXIBAGiOkAsAQHOEXAAAmiPkAgDQHCEXAIDmCLkAADRHyAUAoDlCLgAAzRFyAQBojpALAEBzhFwAAJoj5AIA0Bwhl/ljZCQpZfzW6Qy6GgBgiC3q5kmllNEkP0jyTJKttdZV/SyKBWp09LmPSxlYGQDA8Osq5G7z6lrrI32rBAAAesRyBQAAmtNtyK1JbimlrC+lXDrVE0opl5ZS1pVS1m3evLl3FdJ/nc5za2FHRgZdDQDArHW7XOGUWusDpZRDktxaSvlWrfX2yU+ota5JsiZJVq1aVXtcJ/00NpZU/8kAgHZ01cmttT6w7c+Hk3wuyYn9LAoAAGZj2pBbSnlBKWX/iY+T/HySO/tdGAAA7K1uliu8OMnnyvhIp0VJPlFrvbmvVQEAwCxMG3JrrfclOXYOagEAgJ4wQgwAgOYIuQAANEfIBQCgOUIuAADNEXIBAGiOkMv8NDIyvs1wpzPoSgCAIdTttr4wt0ZHx/8cn88MADAjOrkAADRHyAUAoDlC7kLW6YwvBxgZGXQlAAA9ZU3uQjY2ltQ66CoAAHpOJxcAgOYIuQAANEfIBQCgOUIuAADNEXIBAGiOkAsAQHOEXOa3kZHxWb6ljM/1BQDogjm5zG+jo899XMrAygAAhotOLgAAzRFyAQBojpALAEBzhFwAAJoj5AIA0BwhFwCA5gi5AAA0R8gFAKA5Qi4AAM0RcgEAaI6QCwBAc4RcAACaI+QCANAcIRcAgOYIuQAANEfIBQCgOUIuAADNEXIBAGiOkAsAQHOEXAAAmiPkAgDQHCF3oel0klLGbyMjg64GAKAvFg26AObY2FhS66CrAADoK51cAACaI+QyPEZGxpdZdDqDrgQAmOcsV2B4jI6O/1nKQMsAAOY/nVwAAJoj5AIA0BwhFwCA5gi5AAA0R8gFAKA5Qi4AAM0RcgEAaI6QCwBAc4RcAACaI+QCANAcIRcAgOZ0HXJLKfuWUr5RSvlCPwsCAIDZmkkn94ok9/SrEAAA6JWuQm4pZXmSX0jy0f6WAwAAs9dtJ/d9SX4nybN9rAUAAHpi2pBbSjknycO11vXTPO/SUsq6Usq6zZs396xAeqDTSUoZv42MDLoaAIC+66aTe0qSc0spo0luTHJqKeUvdn5SrXVNrXVVrXXV0qVLe1wmszI2ltQ6fhsdHXQ1AAB9N23IrbX+b7XW5bXWTpLzk/xdrfU/9b0yAADYS+bkAgDQnEUzeXKt9StJvtKXSgAAoEd0cgEAaI6Qy/AZGXluWkSnM+hqAIB5aEbLFWBemDwhopSBlQEAzF86uQAANEfIBQCgOUIuAADNEXIBAGiOkAsAQHOEXAAAmiPkAgDQHCEXAIDmCLkAADRHyAUAoDlCLgAAzRFyAQBojpALAEBzhFwAAJoj5AIA0BwhFwCA5gi5AAA0R8gFAKA5Qm7LOp2klGRkZNCVAADMqUWDLoA+GhtLah10FQAAc04nFwCA5gi5AAA0R8gFAKA5Qi4AAM0RcgEAaI6QCwBAc4RchtvIyPgs4FLG5wIDAMScXIbd6OhzH5cysDIAgPlFJxcAgOYIuQAANEfIBQCgOUIuAADNEXIBAGiOkAsAQHOEXAAAmiPkAgDQHCEXAIDmCLkAADRHyAUAoDlCLu0YGUlKSTqdQVcCAAzYokEXAD0zOjr+ZykDLQMAGDydXAAAmiPkAgDQHCEXAIDmCLkAADRHyG1NpzN+4VUp49MGAAAWINMVWjM2ltQ66CoAAAZKJxcAgOYIuQAANEfIBQCgOUIuAADNEXIBAGiOkAsAQHOmDbmllCWllK+XUr5ZSrmrlPJ/zEVhAACwt7qZk/tkklNrrT8spSxO8tVSypdrrV/rc20AALBXpg25tdaa5IfbPl287Wa3AQAA5q2u1uSWUvYtpWxI8nCSW2ut/9jfsgAAYO91FXJrrc/UWo9LsjzJiaWUo3d+Tinl0lLKulLKus2bN/e6TgAA6NqMpivUWv8tyVeSnDnFY2tqratqrauWLl3ao/IA2tHpJKXseut0Bl0ZQHu6ma6wtJTyY9s+fn6S1yT5Vr8LA2jN2FhS6663sbFBVwbQnm6mKxya5PpSyr4ZD8X/X631C/0tCwAA9l430xXuSHL8HNQCvTEyMv4e8MTHo6MDLQcAmHvddHJhuEwOtRNhFwBYUGzrCwBAc4RcgD6YapLCyMigqwJYOCxXAOiDiUkKAAyGTi4AAM0RcgFmydIEgPlHyAWYpak2eZjt5LqpgrOd0QC6Z00uwDw01ZpeE/EAuifkAgzY5P1LJt8HwN4TclvR6Yy3fpwZYejYlA+g94TcVphXBACwnQvPAABojpALAEBzhFwAAJoj5AIMiYkpDGbnAkzPhWcAQ2KqKQxm5wJMTScXWPDsLgbQHiEXWPCm2pY32TX47u5mPDXA/GO5AsAUhmWDhql2S5u4f1j+DgD9IOQCDLHdBVlrdYGFznIFAACaI+QCANAcIRcAgOYIuQANsnEEsNAJucNs8nBPM4ymNvlM7wzPAjI6uutYtLGxQVcFMHdMVxhmE8M92b3Jl5673BwAFgydXAAAmiPkAgDQHCEXAIDmCLkAADRHyAUAoDlCLrCgTJ68ZwIfQLuMEAMWFJP3ABYGnVyABcIuaMBCopMLsEBM3htlgj1SgFbp5AJNmmrtrfW3AAuHTi7QJGtvARY2nVwAAJoj5AIA0BwhFwCA5gi5AAvYVGPFjBYDWuDCM4AFbKqxYonRYsDw08kFAKA5Qi4Lx8T7st6HBYDmWa7AwjHxvqz3YQGgeTq5AAA0R8gFAKA5Qi4w9DqdXUdgjYwMuioABsmaXGDojY0ltQ66CgDmE51cAHYx1SYRBpMAw0Qnd9h0OuNtq8T7sUDfTLVJhMEkwDARcoeN92UBAKZluQIAXbGEARgmOrkAdMUSBmCY6OSy8ExuR2lDAUCTdHJZeCa3o7ShAKBJOrkAADRn2pBbSnlpKeW2Uso9pZS7SilXzEVhAFOxuxkA3ehmucLWJL9Za/1vpZT9k6wvpdxaa727z7UB7MIUPQC6MW0nt9b6YK31v237+AdJ7knykn4XBgAAe2tGa3JLKZ0kxyf5xykeu7SUsq6Usm7z5s29qQ4AAPZC1yG3lPLCJJ9J8uu11sd2frzWuqbWuqrWumrp0qW9rBEAAGakq5BbSlmc8YB7Q631s/0tCYBhYRc0YL6a9sKzUkpJ8rEk99Ra/+/+lwTAsLALGjBfddPJPSXJBUlOLaVs2HY7u891AQDAXpu2k1tr/WoS/y4H5lynMz4ybDIzcQHohm19gXnLTFwA9pZtfQEAaI6QCwBAc4RcAACaI+QCANAcIRcAgOYIuQAANEfIBQCgOUIuAD01MjK+te/Ot05n0JUBC4nNIADoqdHRqe8v9s4E5pBOLjBwnc7UnT9b+AKwt3RygYGzfS8AvaaTCwBAc4TcYTHxfq73b3tr4goZV8QAQFMsVxgW3s/tj4krZFwRAwBN0ckFAKA5Qi4AAM0RcgEAaI6QC8CcmGonNNd8Av3iwjMA5sRUO6G55hPoF51cAACaI+QCANAcIRcAgOYIuQAANEfIBebUxA7Vk292qwag14Tc+WxyGpAC+mvybCMzjfpqYofqybeprroHgNkwQmw+m0gD9N/klGWmEQAMPZ1cAACaI+QCANAcIReAgbHVL9Av1uQCMDC2+gX6RScXAIDmCLlA35iJy96whAHoBcsVgL4xBY+9YQkD0As6uQAANEfIBQCgOUIuAADNEXIBAGiOkAsAQHOEXNjZ5PlF5hYBwFAyQmw+6nTGZy8ZKDoYk+cXmVsEAENJyJ2PDBcFAJgVyxUAmPem2gXNiiJgT3RyAZj3ptoFLbGiCNg9nVwAAJoj5AIA0BwhF5i1Tmfq9ZIGhAAwKNbkArNmIAgA841OLgAAzRFyAQBojpALwNCaan6u2blAYk0u7NnEGXRkZPeDOoGBmep/S7NzgUTIhT2bOIM6awLAULFcAQCA5gi5wIxMNRPXPFwA5hvLFYAZMRMXgGEwbSe3lPLxUsrDpZQ756KgBWtye0xbDABgVrpZrnBdkjP7XAcT7bFaXcUPADBL04bcWuvtSb4/B7UAAEBPuPAMAIDm9CzkllIuLaWsK6Ws27x5c6++LQAAzFjPQm6tdU2tdVWtddXSpUt79W0BYEZs9QskRogB0Bhb/QJJdyPEPpnkH5IcUUrZVEp5S//LAuYDGz8AMKym7eTWWt84F4UA84+NH2jFxBKGqe43tRHaZLkCdGPyGdJZEYbO7v6XtYwB2iXkQjcmnyGdFQFg3jMnFwCA5gi5AAA0R8gdpMmXrrtkfXhMrM81eBMA5i1rcgfJpevDaWJ9rrW5ADBv6eQCANAcIRcAgOYIuUASu5sB0BZrcoEkloizME21E5r9XqANQi4AC9ZUYdY1pdAGyxUAAGiOkAsAQHOEXACYZGKd7uSbvV9g+FiTC3tr8hUrrlSBZlinC20Qcgeh0xm/lN18puE2+UzoDAgA84qQOwhmNQEA9JU1uQAANEfIhQVmqp3N7G4GQGssV4BeGKKL0KyWAWAhEHKhF1yEBk2z/S8MHyEXAKZhrBgMH2tyAQBojpALAEBzhNy5MvmSdpexAwy9qbb/tQUwzB/W5M4Vl7QzABOb603m31jQG7u76MxaXZgfhFzotYn2zjy49Nq/rQBYqIRc6LWJYKudAwADY00uAADNEXKhEVNt12v9LQALleUK0Ajrb2F+sDsazA9CLgD0kN3RYH6wXAEAgOYIudAvkyfFmw4PC9pUG0f4tQD9ZblCP02exO8KoIVn8nuW3quEBc0SBph7Qm4/uRIIAGAgLFeAITPVqDDjwgBgRzq5MGS8QQAA09PJBYABmOpiNBekQe8Iuf0w8X6y94+ZMHE2c/YCthkdHX9XZufbxPXKwOxYrtAP3k9mZxOXVrucGgDmhE5ur0y+GkgHl92Z4ezcqS4y8/KCtpmpC72hk9srurd0Y4azc72sYOExUxd6QycX5gldWwDoHSEXBmWnpQsTXdvJt6k6OsDCYwkDzJzlCjAotv0FumQJA8ycTi4ADCHdXdgzIRcAhtBUc3YTG0zABMsVYAA6nR0Hvte4yAyYvd2t47e0gYVIJxcGYOeLzDIyktEx7RagPyxtYCHSyZ2tiZacNhyzYUc0oI9cuMZCpJO7NyYPNE3MemK3ppp9u8f5t9otANATOrl7wzZUdGnGLxXtFmCOTPybeuf79GxohU4u9EjfdizbadMIgF6YyXQGExsYRjq53Zp8Obz1t0yhbw1+m0YAc2QmXVy/jpjvdHKnM9GeS+y1ynZ969pOx5pdYJ6Y6teRji/zSVcht5RyZinl26WUe0sp7+53UQPnwjK22d2FY8mub/PNyUtkqvcXJw/cBZgjU/06muq2u19RU/1+FYjppWlDbill3yR/muSsJEcleWMp5ah+FzYQurbsZOd5tvP+ZTHxGnamAOaJ3XV8k+7WBPt1xt7qppN7YpJ7a6331VqfSnJjktf2t6w5pGu7IO2uQzuQJQizNfkMktjbE5hXdtfxnepUO9uL4fzaY7JuQu5Lknxv0uebtt03HKZLM8kQtOcWjm7D52xvSXdvsw3FS2LyWWGi4N2dVZK9P0NM9x9nqjPJ5K9xpgH2QrfLInZ3S+bmvIQ8r94AAASGSURBVOJX3PzTzXSFMsV9u1xDXkq5NMml2z79YSnl27MpbM6MjT2Xeuafg5M8MugiGnNwkkfm93/2eWb6gzX+Op3ueQ76TPh/v/cc095zTCfp0a84x3Tmdvueazchd1OSl076fHmSB3Z+Uq11TZI1My6N3SqlrKu1rhp0HS1xTHvPMe09x7T3HNPec0x7zzHtrW6WK/xTksNLKYeVUp6X5Pwkn+9vWQAAsPem7eTWWreWUi5L8tdJ9k3y8VrrXX2vDAAA9lJXO57VWr+U5Et9roVdWf7Re45p7zmmveeY9p5j2nuOae85pj1Ual/2IQUAgMGxrS8AAM0Rcges2y2TSyknlFKeKaWcN5f1DaNujmkp5VWllA2llLtKKf9lrmscNtMd01LKAaWUvyqlfHPbMX3zIOocJqWUj5dSHi6l3Lmbx0sp5f3bjvkdpZRXzHWNw6aLY7p627G8o5TyX0spx851jcNmumM66XnOUV3q5pg6R/WGkDtA3W6ZvO15/2fGL/5jD7o5pqWUH0vyn5OcW2t9eZJfmvNCh0iXr9O3J7m71npsklcl+b+2TWNh965LcuYeHj8ryeHbbpcm+eAc1DTsrsuej+l3kvwvtdZjkrwn1j9247rs+Zg6R83cddnDMXWO6h0hd7C63TL5HUk+k+ThuSxuSHVzTN+U5LO11u8mSa3Vcd2zbo5pTbJ/KaUkeWGS7yfZOrdlDpda6+0ZP06789okf1bHfS3Jj5VSDp2b6obTdMe01vpfa63/uu3Tr2V87jt70MXrNHGOmpEujqlzVI8IuYM17ZbJpZSXJPlfk3xoDusaZt1sQ/3TSX68lPKVUsr6UsqFc1bdcOrmmF6b5MiMbxTz35NcUWt9dm7Ka9Zwb6k+/70lyZcHXcSwc47qC+eoHulqhBh9082Wye9L8q5a6zPFlqjd6OaYLkryH5OcluT5Sf6hlPK1Wus/97u4IdXNMT0jyYYkpyb5n5LcWkr5+1rrY/0urmFdbanOzJVSXp3xkPvKQdfSAOeo3nOO6hEhd7C62TJ5VZIbt/3yODjJ2aWUrbXWtXNT4tDp5phuSvJIrfXxJI+XUm5PcmwSv0Cm1s0xfXOSq+r4TMJ7SynfSfKyJF+fmxKb1NWW6sxMKeWYJB9Nclatdcug62mAc1TvOUf1iOUKgzXtlsm11sNqrZ1aayfJp5P8ml8ee9TNNtQ3JfmfSymLSin/IclJSe6Z4zqHSTfH9LsZ7zqklPLiJEckuW9Oq2zP55NcuG3Kws8kebTW+uCgixpmpZSfTPLZJBfoivWGc1RfOEf1iE7uAO1uy+RSytu2PW6N0wx1c0xrrfeUUm5OckeSZ5N8tNa6x/E4C1mXr9P3JLmulPLfM/42+7tqrY8MrOghUEr5ZMYnURxcStmU5H9PsjjZfky/lOTsJPcm+feMd8vZgy6O6ZVJDkryn7d1HrfWWlcNptrh0MUxZYamO6bOUb1jxzMAAJpjuQIAAM0RcgEAaI6QCwBAc4RcAACaI+QCANAcIRcAgOYIuQAANEfIBQCgOf8/lfQsmlX8DUMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4dec19e9f9214d40a665d087a6af6dbb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=100000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAHSCAYAAAAezFYoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5RlVX0v+u+Ph+BRVB6NNjbUJqGDIGDrbVCiJ0dFhCYEbhxEQQbgK0QjEc1JoieJJI4ME+4I55CjeEQSHT6CkOMLiUGEKLnoPTEI2hAQTROs1gYGIBLwhdo47x9V1RSrq+jd9dq7qj6fMfaovddae69Zq6t2fXvu35yzWmsBAAAescOgGwAAAMNGSAYAgA4hGQAAOoRkAADoEJIBAKBDSAYAgI6dBt2Aqey1116t1+sNuhkAACxhN9xww3dbayum2jeUIbnX6+X6668fdDMAAFjCqmrjdPuUWwAAQIeQDAAAHUIyAAB0DGVNMgBAkvzsZz/Lpk2b8tBDDw26KSxiu+66a1atWpWdd9657+cIyQDA0Nq0aVN222239Hq9VNWgm8Mi1FrLfffdl02bNmX//ffv+3nKLQCAofXQQw9lzz33FJCZsarKnnvuud2fRgjJAMBQE5CZrZn8DAnJAACPoapy2mmnbXm8efPmrFixIscff3yS5IMf/GDOOuusrZ7X6/Vy6KGHZs2aNVmzZk3e9KY3PeZ5/vRP/zTnnXfenLb9Na95Tfbee+8ccsghj9r+ve99L0cffXRWr16do48+Ovfff/+WfX/xF3+RAw44IAceeGA+97nPbdl+ww035NBDD80BBxyQN73pTWmtzWlbh42QDAAsGr1eUjV3t34W+H3CE56Qm2++OT/+8Y+TJFdffXWe/vSn99Xea665JuvXr8/69evzrne9a+bf+Ay96lWvypVXXrnV9nPPPTdHHXVUNmzYkKOOOirnnntukuTrX/96Lr300txyyy258sor89u//dt5+OGHkyRveMMbctFFF2XDhg3ZsGHDlK+7lAjJAMCisXFj0trc3TZOu97ao61bty7/8A//kCS55JJLcsopp8zjd5n89V//ddatW7clmM/Ur/zKr2SPPfbYavunP/3pnHHGGUmSM844I5dddtmW7SeffHJ22WWX7L///jnggANy3XXX5a677sqDDz6YI488MlWV008/fctzliohGQBgG04++eRceumleeihh3LTTTfluc99bl/Pe9GLXrSl3OL888/v6zkXXHBB/v7v/z6XXXZZHv/4xz9q38UXX7zl9SbfTjrppO36fu6+++6sXLkySbJy5crcc889SZI77rgj++6775bjVq1alTvuuCN33HFHVq1atdX2pcwUcAAA23DYYYdldHQ0l1xySY477ri+n3fNNddkr7326vv4j3zkI1m1alUuu+yyKef0PfXUU3Pqqaf2/Xrba6o646qadvtSJiQDAPThhBNOyO/93u/ln/7pn3LffffNyzkOOeSQrF+/fto5fS+++OL85V/+5VbbDzjggHz84x/v+zxPfepTc9ddd2XlypW56667svfeeycZ6yH+zne+s+W4TZs2ZZ999smqVauyadOmrbYvZcotAAD68JrXvCbnnHNODj300Fm9zgUXXJALLrhgyn3Pfvaz8773vS8nnHBC7rzzzq32n3rqqVsGAk6+bU9ATsYC/4c+9KEkyYc+9KGceOKJW7Zfeuml+clPfpJvfetb2bBhQ4444oisXLkyu+22W7785S+ntZYPf/jDW56zVAnJAAB9WLVqVc4+++wp933wgx/MqlWrttwmel0n1ySffvrpSZJvfOMb2XPPPac9zwte8IKcd955+dVf/dV897vfnVWbTznllBx55JH55je/mVWrVuX9739/kuRtb3tbrr766qxevTpXX3113va2tyVJnvnMZ+blL395Dj744Bx77LF5z3vekx133DFJ8t73vjeve93rcsABB+QXf/EXs27dulm1bdjVMM5xt3bt2nb99dcPuhkAwIDdeuutOeigg7Y87vX6n5GiHyMjyejo3L1eP44//vh88pOfzOMe97iFPfEy1/1ZSpKquqG1tnaq49UkAwCLxkIH2vnwmc98ZtBNoA/KLQAAoENIBgCADiGZMRPrfPazPicAwBKnJpkxE+t8LvGJwQEA+qEnGQAAOoRkpqb8AgCSJE984hMHev677747xx9/fJ71rGfl4IMP3rIs9ujoaD760Y9u8/l33nlnTjrppBmf/4UvfGHmcmre++67Ly960YvyxCc+MWedddaj9t1www059NBDc8ABB+RNb3rTluWwf/KTn+QVr3hFDjjggDz3uc/N6KRpTj70oQ9l9erVWb169ZYFUuaCkMzUJsov5nIySgCYrYlOnLm6LYLOoHPOOSdHH310brzxxnz961/Pueeem6T/kLzPPvts94p882nXXXfNn/3Zn+W8887bat8b3vCGXHTRRdmwYUM2bNiQK6+8Mkny/ve/P7vvvntuu+22vOUtb8lb3/rWJMn3vve9vOMd78i//Mu/5Lrrrss73vGO3H///XPSTiEZAFg8Jjpx5urWZ2fQD37wgxx11FF5znOek0MPPTSf/vSnk4wF1Wc84xl53etel0MOOSSnnnpq/vEf/zHPf/7zs3r16lx33XVJkh/+8Id5zWtek8MPPzzPfvaztzz/lltuyRFHHJE1a9bksMMOy4YNG7Y691133ZVVq1ZteXzYYYclGVs174tf/GLWrFmT888/Pw8//HB+//d/P4cffngOO+ywvO9979vSxkMOOSTJ2MqAL3vZy3Lsscdm9erV+YM/+IO+L/3Pf/7znHHGGfnjP/7jvp8zlSc84Ql5wQtekF133XWr7/PBBx/MkUcemarK6aefnssuuyxJ8ulPfzpnnHFGkuSkk07K5z//+bTW8rnPfS5HH3109thjj+y+++45+uijtwTr2TJwDwBgG3bdddd86lOfypOe9KR897vfzfOe97yccMIJSZLbbrstH/vYx3LRRRfl8MMPz0c/+tF86UtfyuWXX54///M/z2WXXZZ3vvOdefGLX5wPfOAD+Y//+I8cccQReclLXpILL7wwZ599dk499dT89Kc/zcMPP7zVud/4xjfmFa94RS644IK85CUvyatf/erss88+Offcc3PeeedtWZzkoosuypOf/OR85StfyU9+8pM8//nPz0tf+tJUZ1D++vXr87WvfS277LJLDjzwwPzO7/xO9t1338f8/jdv3pxTTz01hxxySP7oj/5oq/1vectbcs0112y1/eSTT96y5PW23HHHHY/6z8CqVatyxx13bNk30caddtopT37yk3Pfffc9anv3ObMlJPNoIyOL4qMnAFhIrbX84R/+Ya699trssMMOueOOO3L33XcnSfbff/8ceuihSZJnPvOZOeqoo1JVOfTQQ7fUzl511VW5/PLLt5QYPPTQQ/n2t7+dI488Mu985zuzadOmvOxlL8vq1au3OvcxxxyT22+/PVdeeWU++9nP5tnPfnZuvvnmrY676qqrctNNN20prXjggQeyYcOG/NIv/dKjjjvqqKPy5Cc/OUly8MEHZ+PGjdsMyb/1W7+Vl7/85VMG5CQ5//zzH/P5/ZioP55sIuBPt++xnjNb2yy3qKp9q+qaqrq1qm6pqrPHt+9RVVdX1Ybxr7tP8/xjq+qbVXVbVfX3XwkGZ3RUHTIAdFx88cW59957c8MNN2T9+vV56lOfmoceeihJsssuu2w5bocddtjyeIcddsjmzZuTjIW8T3ziE1m/fn3Wr1+fb3/72znooIPyyle+Mpdffnke//jH55hjjskXvvCFvOc978maNWuyZs2a3HnnnUmSPfbYI6985SvzkY98JIcffniuvfbardrYWsu73/3uLef41re+lZe+9KVbHTe5vTvuuOOWNj6WX/7lX84111yz5Xvuestb3rKlzZNvE/XT/Vi1alU2bdq05fGmTZuyzz77bNn3ne98J8lYr/YDDzyQPfbY41Hbu8+ZrX5qkjcn+a+ttYOSPC/JG6vq4CRvS/L51trqJJ8ff/woVbVjkvckWZfk4CSnjD8XAGDReOCBB7L33ntn5513zjXXXJON29mhdMwxx+Td7373lp7Pr33ta0mS22+/Pb/wC7+QN73pTTnhhBNy00035Y1vfOOWoLvPPvvkC1/4Qn70ox8lSb7//e/n3//937Pffvtlt912y/e///1HneO9731vfvaznyVJ/u3f/i0//OEP+27j6aefvqWGuuu1r31tjjvuuPzGb/zGlKH6/PPP39Lmybd+Sy2SZOXKldltt93y5S9/Oa21fPjDH86JJ56YJDnhhBO2zFzx8Y9/PC9+8YtTVTnmmGNy1VVX5f7778/999+fq666Ksccc0zf53ws2yy3aK3dleSu8fvfr6pbkzw9yYlJXjh+2IeS/FOSt3aefkSS21prtydJVV06/ryvz0HbAQDm1ebNm7PLLrvk1FNPza/92q9l7dq1WbNmTZ7xjGds1+u8/e1vz5vf/OYcdthhaa2l1+vlM5/5TP7u7/4uf/u3f5udd945T3va03LOOeds9dwbbrghZ511Vnbaaaf8/Oc/z+te97ocfvjh+dnPfpaddtopz3rWs/KqV70qZ599dkZHR/Oc5zwnrbWsWLFiy8C3ftx0001ZuXLltPt/93d/Nw888EBOO+20XHzxxdlhh5nP/9Dr9fLggw/mpz/9aS677LJcddVVOfjgg/Pe9743r3rVq/LjH/8469aty7p165KMhfTTTjstBxxwQPbYY49ceumlScZ62N/+9rfn8MMPTzI2E8gee+wx43ZNVlPVckx7cFUvybVJDkny7dbaUybtu7+1tnvn+JOSHNtae93449OSPLe19uhJ8TrWrl3b5nI+Ph5DrzdWXjEyMlZqMbEtGXtcNTb6FwAG4NZbb81BBx30yIaJv1tzZfLfvynceOON+c3f/M1pe1iXigcffDCvfe1r87GPfWzQTZk3W/0sJamqG1pra6c6vu+Be1X1xCSfSPLm1tqDfRZFT3XQlImrqs5McmaS7Lfffv02i9mamEpnssd4swCAgVrAv1EXXnhh3vWud+Wv/uqvFuycg/KkJz1pSQfkmeirn7yqds5YQL64tfbJ8c13V9XK8f0rk9wzxVM3JZk8XHJVkjunOkdr7aLW2trW2toVK1b0234AgHnx+te/Pl//+tenHPzG0tfP7BaV5P1Jbm2t/Y9Juy5Pcsb4/TOSfHqKp38lyeqq2r+qHpfk5PHnAQDA0OqnJ/n5SU5L8uKqWj9+Oy7JuUmOrqoNSY4ef5yq2qeqrkiS1trmJGcl+VySW5P879baLfPwfQAAS9T2jJ+CqczkZ6if2S2+lKlri5PkqCmOvzPJcZMeX5Hkiu1uGQCw7O2666657777sueee87ZIhEsL6213HfffVstg70tVtwDAIbWxAIT995776CbwiK26667PmrJ634IyQDA0Np5552z//77D7oZLEMznwUaAACWKCEZAAA6hGQAAOgQkgEAoENIBgCADiEZAAA6hGQe28hI0usNuhUAAAvKPMk8ttHRxApHAMAyoycZAAA6hGQAAOgQkpezXm+s5hgAgEdRk7ycbdyYtDboVgAADB09yQAA0CEkAwBAh5AMAAAdQjL96fXG5ku2sAgAsAwYuEd/Jgb5WVgEAFgG9CQDAECHkAwAAB3KLdg2C44AAMuMkMy2jY4OugUAAAtKuQUAAHQIyQAA0CEkL0cTcx6rNQYAmJKa5OVoYs5jAACmpCcZAAA6hGQAAOgQkgEAoENIBgCADiEZAAA6hGQAAOgQkgEAoENIXm56vdktIjIyMvYaAABLmMVElpvZLiQyOjq2Wh8AwBKmJxkAADqEZAAA6BCSAQCgQ0gGAIAOIRkAADqEZAAA6NjmFHBV9YEkxye5p7V2yPi2v0ty4PghT0nyH621NVM8dzTJ95M8nGRza23tHLUbAADmTT/zJH8wyQVJPjyxobX2ion7VfXfkzzwGM9/UWvtuzNtIAAALLRthuTW2rVV1ZtqX1VVkpcnefHcNgsAAAZntjXJ/znJ3a21DdPsb0muqqobqurMx3qhqjqzqq6vquvvvffeWTYLAABmbrYh+ZQklzzG/ue31p6TZF2SN1bVr0x3YGvtotba2tba2hUrVsyyWQAAMHMzDslVtVOSlyX5u+mOaa3dOf71niSfSnLETM8HAAALZTY9yS9J8o3W2qapdlbVE6pqt4n7SV6a5OZZnA8AABbENkNyVV2S5J+THFhVm6rqteO7Tk6n1KKq9qmqK8YfPjXJl6rqxiTXJfmH1tqVc9d0AACYH/3MbnHKNNtfNcW2O5McN37/9iTPmmX7AABgwVlxDwAAOoRkAADoEJIBAKBDSF5Oer1kZGT2rzMyMvZaAABL1DYH7rGEbNyYtDb71xkdTapm/zoAAENKTzIAAHQIyQAA0CEkAwBAh5AMAAAdQjIAAHQIyQAA0CEkAwBAh5AMAAAdQjIAAHQIyQAA0CEkAwBAh5AMAAAdQjIAAHQIyQAA0CEkLwe9XlKVjIwMuiUAAIvCToNuAAtg48aktUG3AgBg0dCTDAAAHUIyAAB0CMkAANAhJAMAQIeQDAAAHUIyAAB0CMkAANAhJAMAQIeQDAAAHUIyAAB0CMnMzMhIUpX0eoNuCQDAnNtp0A1gkRodHftaNdBmAADMBz3JAADQISQDAECHkAwAAB1C8lLX640NsgMAoG8G7i11GzcmrQ26FQAAi4qeZAAA6BCSAQCgQ0gGAICObYbkqvpAVd1TVTdP2vanVXVHVa0fvx03zXOPrapvVtVtVfW2uWw4AADMl356kj+Y5Ngptp/fWlszfruiu7OqdkzyniTrkhyc5JSqOng2jQUAgIWwzZDcWrs2yfdm8NpHJLmttXZ7a+2nSS5NcuIMXgcAABbUbGqSz6qqm8bLMXafYv/Tk3xn0uNN49umVFVnVtX1VXX9vffeO4tmAQDA7Mw0JL83yS8mWZPkriT/fYpjaopt007Y21q7qLW2trW2dsWKFTNsFgAAzN6MQnJr7e7W2sOttZ8n+euMlVZ0bUqy76THq5LcOZPzAQDAQppRSK6qlZMe/nqSm6c47CtJVlfV/lX1uCQnJ7l8JucDAICFtM1lqavqkiQvTLJXVW1K8idJXlhVazJWPjGa5LfGj90nyd+01o5rrW2uqrOSfC7Jjkk+0Fq7ZV6+CwAAmEPV2rRlwgOzdu3adv311w+6GUtDVTKf/8bz/foAAPOkqm5ora2dap8V9wAAoENIBgCADiEZAAA6hGQAAOgQkgEAoENIZnZGRpJeb9CtAACYU9ucJxke0+jo2DRwAABLiJ5kAADoEJIBAKBDSAYAgA4hGQAAOoRkAADoEJIBAKBDSAYAgA4hmdmzoAgAsMRYTGQp6/XGAux8s6AIALDECMlL2caNSWuDbgUAwKKj3AIAADqEZAAA6BCSAQCgQ0gGAIAOIRkAADqEZAAA6BCSAQCgQ0gGAIAOIRkAADqEZAAA6BCSAQCgQ0heinq9pCoZGRl0SwAAFqWdBt0A5sHGjUlrg24FAMCipScZAAA6hGQAAOgQkgEAoENIBgCADiEZAAA6hGQAAOgQkgEAoENIBgCADiEZAAA6hGQAAOgQkgEAoGObIbmqPlBV91TVzZO2/WVVfaOqbqqqT1XVU6Z57mhV/WtVra+q6+ey4QyZkZGkKun1Bt0SAIBZ66cn+YNJju1suzrJIa21w5L8W5L/9hjPf1FrbU1rbe3MmsiiMDqatJZs3DjolgAAzNo2Q3Jr7dok3+tsu6q1tnn84ZeTrJqHtgEAwEDMRU3ya5J8dpp9LclVVXVDVZ05B+cCAIB5t9NsnlxVf5Rkc5KLpznk+a21O6tq7yRXV9U3xnump3qtM5OcmST77bffbJoFAACzMuOe5Ko6I8nxSU5trbWpjmmt3Tn+9Z4kn0pyxHSv11q7qLW2trW2dsWKFTNtFgAAzNqMQnJVHZvkrUlOaK39aJpjnlBVu03cT/LSJDdPdSwAAAyTfqaAuyTJPyc5sKo2VdVrk1yQZLeMlVCsr6oLx4/dp6quGH/qU5N8qapuTHJdkn9orV05L98FAADMoW3WJLfWTpli8/unOfbOJMeN3789ybNm1ToAABgAK+4BAECHkAwAAB1CMgAAdAjJzK2RkaTXG3QrAABmZVaLicBWRkeTqkG3AgBgVvQkLzW93lhvLgAAM6YneanZuDGZegFEAAD6pCcZAAA6hGQAAOgQkgEAoENIBgCADiEZAAA6hGQAAOgQkgEAoENIBgCADiEZAAA6hGQAAOgQkgEAoENIBgCADiEZAAA6hGQAAOgQkgEAoENIBgCADiEZAAA6hGQAAOgQkgEAoENIBgCADiEZAAA6hOSlpNdLRkYG3QoAgEVvp0E3gDm0cWPS2qBbAQCw6OlJBgCADiGZuTcyMlb6AQCwSCm3YO6NjiZVg24FAMCM6UkGAIAOIRkAADqEZAAA6BCSAQCgQ0gGAIAOIRkAADqEZAAA6BCSAQCgY5shuao+UFX3VNXNk7btUVVXV9WG8a+7T/PcY6vqm1V1W1W9bS4bDgAA86WfnuQPJjm2s+1tST7fWlud5PPjjx+lqnZM8p4k65IcnOSUqjp4Vq0FAIAFsM2Q3Fq7Nsn3OptPTPKh8fsfSvJ/T/HUI5Lc1lq7vbX20ySXjj8PAACG2kxrkp/aWrsrSca/7j3FMU9P8p1JjzeNbwMAgKE2nwP3aoptbdqDq86squur6vp77713HpsFAACPbaYh+e6qWpkk41/vmeKYTUn2nfR4VZI7p3vB1tpFrbW1rbW1K1asmGGzAABg9mYaki9Pcsb4/TOSfHqKY76SZHVV7V9Vj0ty8vjzAABgqPUzBdwlSf45yYFVtamqXpvk3CRHV9WGJEePP05V7VNVVyRJa21zkrOSfC7JrUn+d2vtlvn5NgAAYO7stK0DWmunTLPrqCmOvTPJcZMeX5Hkihm3jsVrZCSpGvs6Ojro1gAAbJdthmSYkYlgXFON3wQAGG6WpQYAgA4hGQAAOoRkAADoEJIBAKBDSAYAgA4hGQAAOoTkpaDXe2ROYgAAZs08yUvBxo1Ja4NuBQDAkqEnGQAAOoRkAADoEJIBAKBDSAYAgA4hGQDYysTESVVj92G5MbsFALCVyRMnVQ22LTAIepIBAKBDSGZ+jYz4nA5gkZhcYmF9KpY75RbMr9FRn9MBLBLTrU01MvLIW/nIyNhbOyx1QjIALGO93lg4TqbvPZ4civV7sFwIyQCwjE3XewzLnZrkxa7XUzgGADDH9CQvdroAAADmnJ5kAADoEJIBYJkx1Rtsm3ILAFhmVOrBtulJBgD6NjFncpW1olja9CQDAH0zZzLLhZ5kAGBGJvcqT77pYWYp0JMMAMtAPyvrba/plqfWw8xSICQDwDJgsB5sH+UWAMCcMriPpUBPMgAsUfNRYtEPg/tYCoRkAFiilFjAzCm3AACADiEZAAA6hGQAAOgQkhezXm9hR2LM1MiI4c0AwKJi4N5itlhGZIyOGt4MsEAGNaMFLDVCMgAsIYul/wSGnXILAADoEJIBAKBjxiG5qg6sqvWTbg9W1Zs7x7ywqh6YdMw5s28yAADMrxnXJLfWvplkTZJU1Y5J7kjyqSkO/WJr7fiZngcAABbaXJVbHJXk31trG+fo9QAAYGDmKiSfnOSSafYdWVU3VtVnq+qZc3Q+AGARGBkZmwW0ypT5LC6zngKuqh6X5IQk/22K3V9NMtJa+0FVHZfksiSrp3mdM5OcmST77bffbJsFAAyB0dFH7psyn8VkLnqS1yX5amvt7u6O1tqDrbUfjN+/IsnOVbXXVC/SWruotba2tbZ2xYoVc9AsAFgeer1HemstIAJzYy4WEzkl05RaVNXTktzdWmtVdUTGQvl9c3BOAGCcBURg7s0qJFfVf0pydJLfmrTt9UnSWrswyUlJ3lBVm5P8OMnJrfk1BgBguM0qJLfWfpRkz862CyfdvyDJBbM5BwCwNEwM4pu4P7leGYbNXJRbAABsk0F8LCaWpQYAgA4heTGaGMZsCDMAwLxQbrEYGcYMADCv9CSzMCZGa1huCQBYBPQkszAmRmsYqQEwJ3q9sQ8WE9V3MB+EZABYhFTewfxSbgEAAB1CMgCw4CaGqhiuwrBSbgEALDgLizDs9CQDwCIxMU2+qfJh/ulJBoBFwmA9WDh6kgEAoENIBgCADiEZAAA6hGQAAOgQkgEAoENIBgCADiEZABioyavvWYGPYWGeZABgoCavvpdYgY/hoCcZAAA6hGQAGFKTl6G2FDUsLOUWADCkLEMNg6MnGQAAOoRkAGCoTJ7twkwXDIpyi8Wm11OUBsCSNnm2CzNdMChC8mKjQA0AYN4pt2BhjYz47AwAGHp6kllYo6M+OwMAhp6eZAAA6BCSAQCgQ0gGAIAOIRkAhsjkpajN+AmDY+AeAAwRM33CcNCTDAAMLavvMSh6kgGAoWX1PQZFTzIAAHQIyQAA0CEkAwBAh5DMwhsZMfoCABhqBu6x8EZHjb4AAIaanuTFpNczszzAEmQBERg+s+pJrqrRJN9P8nCSza21tZ39leR/JjkuyY+SvKq19tXZnHNZM8M8wJLk7R2Gz1yUW7yotfbdafatS7J6/PbcJO8d/woAsF0mFhaZuD95DmWYa/Ndk3xikg+31lqSL1fVU6pqZWvtrnk+LwCwxFhYhIU025rkluSqqrqhqs6cYv/Tk3xn0uNN49u2UlVnVtX1VXX9vffeO8tmAQDAzM02JD+/tfacjJVVvLGqfqWzf6r/501ZddVau6i1tra1tnbFihWzbBYAAMzcrEJya+3O8a/3JPlUkiM6h2xKsu+kx6uS3DmbcwIAwHybcUiuqidU1W4T95O8NMnNncMuT3J6jXlekgfUIwMAMOxmM3DvqUk+NTbLW3ZK8tHW2pVV9fokaa1dmOSKjE3/dlvGpoB79eyay5IxMUTZ8GQAYAjNOCS31m5P8qwptl846X5L8saZnoMlbCIYG54MAAwhK+4BAECHkAwAAB1CMgAAdAjJADAAvd7YsIyJMcxsn4nx36oW0/8AAAy3SURBVFVj1xLm2nwvSw0ATGHjxqRNubwW/bBENfNNT/JiMNHdoKsBAGBB6EleDHQ3AAAsKD3JAADQISQDAECHkAwALGpmumA+qEkGABY1M10wH/QkD7tez6wWAAALTE/ysDOzBQDAgtOTDAALxCp7sHjoSQaABeLDQVg89CQDAECHkAwA80iJBSxOyi0AYB4psYDFSU8yAAB0CMkAANAhJAMAQIeQDAAsGSMjjwyU7PUG3RoWMwP3AIAlY3T0kftVA2sGS4CeZAAA6BCSGayREZ+HAQBDR7kFgzU66vMwAGDo6EkGAIAOIRkAWJLMdMFsKLcAgDnW640tR52MBTUGw0wXzIaQPMx6veXx7joxeG/yuxnAIrZxY9LaoFsBzIaQPMyWy7uswXsAwJBRkwwAAB1CMgAAdAjJAMCSZ6YLtpeaZACYA2a0GG5mumB7CckAMAeWy1hrWC6UWwAAQIeQPIx6vbHPgnxeBwAwEMothpHP7AAABkpPMgAAdAjJDIeJuXnMywMADIEZh+Sq2reqrqmqW6vqlqo6e4pjXlhVD1TV+vHbObNrLkvW6OhYicnE/EkAAAM0m5rkzUn+a2vtq1W1W5Ibqurq1trXO8d9sbV2/CzOAwAAC2rGPcmttbtaa18dv//9JLcmefpcNQwAAAZlTmqSq6qX5NlJ/mWK3UdW1Y1V9dmqeuZcnA8AYKYsUU0/Zj0FXFU9Mcknkry5tfZgZ/dXk4y01n5QVccluSzJ6mle58wkZybJfvvtN9tmAcC8sxT14mSJavoxq57kqto5YwH54tbaJ7v7W2sPttZ+MH7/iiQ7V9VeU71Wa+2i1tra1traFStWzKZZALAgJqa1b+3RwQtY/GYzu0UleX+SW1tr/2OaY542flyq6ojx890303MCAMBCmE25xfOTnJbkX6tq/fi2P0yyX5K01i5MclKSN1TV5iQ/TnJya5aS4zGMjIx9fqlLBgAYoBmH5Nbal5I8ZiVPa+2CJBfM9BwsQ6OjCsSAoaYOGZaHWQ/cY471et51AYbYRB0yS8PETBcT932QyQQhedh49wWABWOmC6YzJ/MkAwDAUiIkAwBAh5AMAAAdQjLDZ2IaOABYQJarZjID9xg+poEDhoxp35YHg/iYTEgeJqZ/AxhKJh6C5UdIHibehQEAhoKaZAAA6BCSGU4ToyeMnABgAAziQ7kFw2li9ISREwAMgEF86EkGgCn0eo/0JBpTDcuPnmQAmIKx1LC86UkGAIAOIRkA4DEYxLc8KbcAAHgMBvEtT3qSh8HE6BAjQwAAhoKQPGgTn9u09uj/qjJmZMRnWwDAglNuMWiGTz+20VGfbQELptcbe1tOfLjH1Cbqkyfu699auoRkABin34JtUZ+8fCi3YPgpuQAAFpieZIafkgsAYIHpSQZgWbP8NDNl/uSlTU8yi8NEyYUREsAcU4fMTKlPXtqE5EHq9XRb9EvJBTCHzGIBbIuQPEi6LwAGwtsvc83UcEuPkAwAMEtKL5YeA/cAAKBDSAZgWTCLBbA9hORBmHin9i69fSYKvsyzA8zARB1ya+pFmV+mhlsa1CQPghEjMzPxV02xFwBDbPJ/wib6xSYY1Ld46EleaKZ9mz3LVAN9UmLBoI2OPvIJRmuPTD3I8NOTvND0Is+eOZOBPnnLBWZKTzIAAHQIyQtJqcXcUXIBTEOJBTAXhOSFMLlqX7X+3BgdVdgFTMksFgwzM18sHkLyQph4x/ZuDTAv9B6zWEweyJcIzMNMSGbxMm8yME7vMYuRwDzchOT51Ov5SZ9PE+8uiesMy8zknmO9xywFAvPwMQXcfNq4ceyd27v3/DIlHCwLvd4jQxFGRkztxtI1+dMQf94GZ1Y9yVV1bFV9s6puq6q3TbG/qupd4/tvqqrnzOZ8i8rETBajoz77WwiTZ7vQgw9LxuQe40RJBcvP5IF+k2/+zM2/GYfkqtoxyXuSrEtycJJTqurgzmHrkqwev52Z5L0zPd+iYSaLweguWW3mC1i0BGN4RHfFvqlW7pv8OyM8z53Z9CQfkeS21trtrbWfJrk0yYmdY05M8uE25stJnlJVK2dxzuHSLYqb/I7u3XzhTbyTjI6aRxkWganeQgVj6M/kHuZk6npmPc+zM5uQ/PQk35n0eNP4tu09ZnGY+Ah/ui4O7+jDZXLP8rbeLSb+Tb2DwHabrgdrugA8XRj2NgrbZ3IP8+Tfmel6npNt/0728ydzOt3f+aXwJ3U2A/emKiXvDqPo55ixA6vOzFhJRpL8oKq+OYu27ZXku7N4fn82blxqFfULc92GwVT/drP791w+125uuW4zM5TXbXt/hQbwFjqU120RcN1mZslct5n8rs7i93uhr9u0syvMJiRvSrLvpMerktw5g2OSJK21i5JcNIv2bFFV17fW1s7Fay0nrtvMuXYz47rNjOs2M67bzLhuM+O6zcwwXbfZlFt8Jcnqqtq/qh6X5OQkl3eOuTzJ6eOzXDwvyQOttbtmcU4AAJh3M+5Jbq1trqqzknwuyY5JPtBau6WqXj++/8IkVyQ5LsltSX6U5NWzbzIAAMyvWS0m0lq7ImNBePK2Cyfdb0neOJtzzNCclG0sQ67bzLl2M+O6zYzrNjOu28y4bjPjus3M0Fy3apYsAgCAR5nVinsAALAULeqQvK1lsScdd3hVPVxVJy1k+4ZVH8uJv7CqHqiq9eO3cwbRzmHTz8/b+LVbX1W3VNX/u9BtHEZ9/Lz9/qSftZvHf1f3GERbh0kf1+3JVfX3VXXj+M+bMR/p67rtXlWfqqqbquq6qjpkEO0cNlX1gaq6p6punmZ/VdW7xq/rTVX1nIVu4zDq47o9o6r+uap+UlW/t9DtG1Z9XLdTx3/Obqqq/1NVz1roNiZJWmuL8paxwYL/nuQXkjwuyY1JDp7muC9krHb6pEG3e9C3fq5bkhcm+cyg2zpMtz6v21OSfD3JfuOP9x50uwd96/f3dNLxv5bkC4Nu96Bvff68/WGS/2f8/ook30vyuEG3fRFct79M8ifj95+R5PODbvcw3JL8SpLnJLl5mv3HJflsxtY/eF6Sfxl0m4fh1sd12zvJ4UnemeT3Bt3eYbn1cd1+Ocnu4/fXDernbTH3JPezLHaS/E6STyS5ZyEbN8T6vW48Wj/X7ZVJPtla+3aStNb8zG3/z9spSS5ZkJYNt36uW0uyW1VVkidmLCRvXthmDp1+rtvBST6fJK21byTpVdVTF7aZw6e1dm3Gfoamc2KSD7cxX07ylKpauTCtG17bum6ttXtaa19J8rOFa9Xw6+O6/Z/W2v3jD7+csXU2FtxiDsnbXPK6qp6e5NeTXBgm9LtU+JHjH+N+tqqeuTBNG2r9XLdfSrJ7Vf1TVd1QVacvWOuGV99L01fVf0pybMb+U7vc9XPdLkhyUMYWaPrXJGe31n6+MM0bWv1ctxuTvCxJquqIjK22NZA/wItM37/LMMdem7FPMRbcrKaAG7B+lrz+qyRvba09XEtr+ejZ6Oe6fTXJSGvtB1V1XJLLkqye95YNt36u205J/q8kRyV5fJJ/rqovt9b+bb4bN8T6Xpo+Y6UW/19r7bF6s5aLfq7bMUnWJ3lxkl9McnVVfbG19uB8N26I9XPdzk3yP6tqfcb+c/G16IHvx/b8LsOcqKoXZSwkv2AQ51/MIbmfJa/XJrl0PCDvleS4qtrcWrtsYZo4lLZ53Sb/kW2tXVFV/6uq9mqtLYk16Geo32XYv9ta+2GSH1bVtUmelWQ5h+S+l6bP2KqdSi3G9HPdXp3k3DZWtHdbVX0rYzW21y1ME4dSv+9vr07GBqMl+db4jce2Pb/LMGtVdViSv0myrrV23yDasJjLLba5LHZrbf/WWq+11kvy8SS/vcwDctLHdauqp43/8Zj4OHKHJAP5AR0i/SzD/ukk/7mqdhovHXhuklsXuJ3Dpp/rlqp6cpL/krFrSH/X7dsZ+9Qi4zW1Bya5fUFbOXz6eX97yvi+JHldkmuXee97vy5Pcvr4LBfPS/JAa+2uQTeKpamq9kvyySSnDfLT2EXbk9z6Wxabjj6v20lJ3lBVm5P8OMnJ471Vy1Y/1621dmtVXZnkpiQ/T/I3rbUpp7dZLrbj9/TXk1w13gu/7PV53f4syQer6l8z9lH4W5f5pz39XreDkny4qh7O2Gw0rx1Yg4dIVV2SsZmN9qqqTUn+JMnOyZbrdkXGZri4LcmPMt4bv9xt67pV1dOSXJ/kSUl+XlVvztiMK8v6P2Z9/Lydk2TPJP9rvM9uc2tt7YK3c5lnHwAA2MpiLrcAAIB5ISQDAECHkAwAAB1CMgAAdAjJAADQISQDAECHkAwAAB1CMgAAdPz/OJxL0NzNcOAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_range(3)\n",
    "plot_range(10)\n",
    "plot_range(100)\n",
    "plot_range(1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the James-Stein estimator seems to have its mean square error go towards 0.5, while the MLE estimator has its error go towards 1."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}