{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"tags": [
"remove_cell"
]
},
"outputs": [],
"source": [
"# HIDDEN\n",
"import warnings\n",
"warnings.filterwarnings('ignore')\n",
"from datascience import *\n",
"from prob140 import *\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"plt.style.use('fivethirtyeight')\n",
"%matplotlib inline\n",
"from scipy import stats"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Chi-Squared Distributions ##"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The gamma family has two important branches. The first consists of gamma $(r, \\lambda)$ distributions with integer shape parameter $r$, as you saw in the previous section.\n",
"\n",
"The other important branch consists of gamma $(r, \\lambda)$ distributions that have *half-integer* shape parameter $r$, that is, when $r = n/2$ for a positive integer $n$. Notice that this branch contains the one above: every integer $r$ is also half of the integer $n = 2r$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Chi-Squared $(1)$ ###\n",
"\n",
"We have already seen the fundamental member of the branch. Let $Z$ be a standard normal random variable and let $V = Z^2$. By the change of variable formula for densities, we found the density of $V$ to be\n",
"\n",
"$$\n",
"f_V(v) ~ = ~ \\frac{1}{\\sqrt{2\\pi}} v^{-\\frac{1}{2}} e^{-\\frac{1}{2} v}, ~~~~ v > 0\n",
"$$\n",
"\n",
"That's the gamma $(1/2, 1/2)$ density. It is also called the *chi-squared density with 1 degree of freedom,* which we will abbreviate to chi-squared (1)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"tags": [
"remove-input",
"hide-output"
]
},
"outputs": [
{
"data": {
"text/html": [
"\n",
" VIDEO"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# VIDEO: Chi-Squared Distributions\n",
"from IPython.display import YouTubeVideo\n",
"\n",
"YouTubeVideo('TXM-yzqYmwg')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### From Chi-Squared $(1)$ to Chi-Squared $(n)$ ###\n",
"\n",
"When we were establishing the properties of the standard normal density, we discovered that if $Z_1$ and $Z_2$ are independent standard normal then $Z_1^2 + Z_2^2$ has the exponential $(1/2)$ distribution. We saw this by comparing two different settings in which the Rayleigh distribution arises. But that wasn't a particularly illuminating reason for why $Z_1^2 + Z_2^2$ should be exponential. \n",
"\n",
"But now we know that the sum of independent gamma variables with the same rate is also gamma; the shape parameter adds up and the rate remains the same. Therefore $Z_1^2 + Z_2^2$ is a gamma $(1, 1/2)$ variable. That's the same distribution as exponential $(1/2)$, as you showed in exercises. This explains why the sum of squares of two i.i.d. standard normal variables has the exponential $(1/2)$ distribution.\n",
"\n",
"If $Z_1, Z_2, Z_3$ are i.i.d. standard normal variables, then:\n",
"\n",
"- $Z_1^2$ has the gamma $(1/2, 1/2)$ distribution\n",
"- $Z_1^2 + Z_2^2$ has the gamma $(1/2 + 1/2, 1/2)$ distribution\n",
"- $Z_1^2 + Z_2^2 + Z_3^2$ has the gamma $(1/2 + 1/2 + 1/2, 1/2)$ distribution\n",
"\n",
"Now let $Z_1, Z_2, \\ldots, Z_n$ be i.i.d. standard normal variables. Then $Z_1^2, Z_2^2, \\ldots, Z_n^2$ are i.i.d. chi-squared $(1)$ variables. That is, each of them has the gamma $(1/2, 1/2)$ distribution. \n",
"\n",
"By induction, $Z_1^2 + Z_2^2 + \\cdots + Z_n^2$ has the gamma $(n/2, 1/2)$ distribution. This is called the *chi-squared distribution with $n$ degrees of freedom,* which we will abbreviate to chi-squared $(n)$.\n",
"\n",
"In data science, these distributions often arise when we work with the *sum of squares of normal errors*. This is usually part of a *mean squared error* calculation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Chi-Squared Distribution with $n$ Degrees of Freedom ###\n",
"For a positive integer $n$, the random variable $X$ has the *chi-squared distribution with $n$ degrees of freedom* if the distribution of $X$ is gamma $(n/2, 1/2)$. That is, $X$ has density\n",
"\n",
"$$\n",
"f_X(x) ~ = ~ \\frac{\\frac{1}{2}^{\\frac{n}{2}}}{\\Gamma(\\frac{n}{2})} x^{\\frac{n}{2} - 1} e^{-\\frac{1}{2}x}, ~~~~ x > 0\n",
"$$\n",
"\n",
"Here are the graphs of the chi-squared densities for degrees of freedom 2 through 5."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"tags": [
"remove_input"
]
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEZCAYAAABFFVgWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeXhU1fnA8e+5s2WyEiAQSNgCYQcRZJNVZFNRinvd6tqqRbu7VaxVi7a1/dVWq3WlaitWcRcElB0EEWTfwhIgCYTsZJ/t/P64A4aQkJlkQjLk/TzPPGRmzr3nvczyzlnuuUprjRBCCNEYjKYOQAghxLlLkowQQohGI0lGCCFEo5EkI4QQotFIkhFCCNFoJMkIIYRoNJJkhBBCNBpJMkIIIRqNJJl6UEotU0q9WkeZOUqpL89WTM2RUupxpdTeAMpFK6UylVJDg9j3v5RSzzYswqYR6HvjbL2HlFKzlVLZSimtlLq1sesTLYskmWqUUm2UUn9SSu1WSlUopY4ppVYopW5RSlmD2NXPgGvqqMuplHpSKZWmlCpXSuUppdYrpe5v2FGEnQeBb7XW64PY5gngHqVUypkK+b+otf/mUUrlK6W+Vkr9TinVukFR198p7w2l1JdKqTl1lWsMSqnhwMPAj4EOwLuNWV9jUkr9xv/aFiilCpVSq5RSU4PY/mal1Ab/9uVKqZ1KqV8ppVQ945mglPIG8kOryjaPV3m/Vr31qE8MzUEwX5rnPKVUMrAa8ACPAd8BbuBC4NfAFmBTIPvSWhcFUOxF4CLML5PNQCxwPtA52Ngbg1LKrrV2NXIdEcA9wC3BbKe1zlRKfQXci/nanMlK4FrMH1XxwHDgAeBupdQ4rfWeoANvgADfGwGXa6BUwKe1/ri+Ozgb75MATQBeB9YD5cBdwGf+13h1ANsfA54EdgOVwBjgn5jfB88FE4hSqj3wb2AxEGyCSAdGVnssJ8h9NB9aa7n5b8CnwFEgrobnbECU/+9lwKvALH/5fGDOief9ZeYAX9ZRXyEws44yDsxkVAQU+P9+Gthbpcwy4NVq2z0KpFe5P8lfLt+/r+XAsGrbLANew/ygHQFyqjx3H7ALqADSgN8C1mDirOX4fgCUVttXPKCBG4H3gONANvCTatveBhytY/81vg6YCX0fsKSG52o91kBee3+50Zg/WIr9t83AlOox+f/W1W7ja4s9gNeh1npr+b85pe4q7/VngEzABewAbgjkfVJDHQG9lo34md4K/KUB238IfBjkNgbwJfAQ8Hhdn4Fq2wZVPhxu0l3m5+86uRR4XtfwC1Jr7dZal1Z56GqgNTAeuAHzy/KBIKs9Akyto9vmGeAqzF/6IzG/kH8aZD0A0cALwAjMllka8IVSqk21ctcCCcDFmL8MUUo9jtlaeBjog9ny+gnwuxDEOQ74TmvtqfLY+f5/fwm8DQwC3gSeV0pFVSm3DmivlOoTQD2n0Fofx0yE45VSCSceD/BYz/jaK6UswCf++Ab7b48DZTWE8jPMltb/MLurOgBraoq5rtiCrPdE3T8HvFXqBpiN2Qr4OdAf8zV4Wyl1cZVtT3uf1CLQ1/LEMT6ilCqp4/bIGeqrui8DiAFyAylfbVullBoGjAKWBrn5LMzE+qdg6/VLVkpl+G8LlFIX1nM/zUNTZ7nmcgOGYb4xrgyg7DJgS7XHXgK+rnJ/DnW3ZEYBBzE/5FuAl4HpgPI/H4X5i/Wuatt9S5AtmRrqNjBbHDdW288ewKjyWCTml9TUatvfAhQGE2ctcXwEvFvtsV9hdlOmVnlsgP/1Sa7yWKz/scvOsP9aXwdgqn/7YUEcayCv/Ylf7+MDiQnzV++cM5ULMLYz1ltLLLcCnmqvdyVwb7VyH+Jv9dX0PjnD/gN6Las81xqze+lMt9YBHtujmL0Fp9Vzhm3igBLMFpwXeCzQbf3bX4T54zHRf//xuj4D1ba/BDOBD8TsrvuvP45JwcTRnG4yJvO9E4N7gV77oPrYTCYwucYdKzUGWFDlodla69la69VKqe6YCW4kMBaYByxQSl0BdMfshqr+y3YVMC3AOE/E0A1zsHwk0A4zyUQCXaoV3aC19lW53w9wAvOUUlX/byxAhL8V0KEBcToxu9iqOh9YqrVOq/JYKuaXbFaVxyqq7KM+qr/mgRwr1PHaa60L/LMPFyqllmB2TX6otd5dzzgDik1rnROCensAdmBFtceXY7agTqj+PqlNoK8lAFrrfMwuyAZRSt0LPAJcobXOCGLTYszWViRmi/9ppVSW1vqMs0n9dbbFbK3drrU+Wo+w0VovqPbQSqVUEvAbzPGdsCNJ5ntpgA/zw/xhAOWrD3Rqap+t9y3mG/eEkx8ibXYTrfHf/qKUugl4CzPhFFbZ95n4+P4L8wRbtfufYXYb/BQ47I9/FeYXSlWl1e6fOKZrMH+9VpcPdAwwzprkYP56rep8zP77qgZjtiCqfrGd2K6+g6L9MWPe778fyLFCAK+91voupdRzmMlnEvCkUmqm1vpf9Yw1oNhCWG/111JVe6z6+6Q2gb6WZiVmV1hd3WGztdaza3tSKfVr4PeYCSaoKeD+mE7MBtuilIoHnsIch6tLf8zPwqdVJqQZZkjKA9yitf5vMPH4fQ1cWY/tmgVJMn5a63yl1AJgplLqH7rauIxSygbY9anjMoHuu5zv37h12en/tx1mcnJhdqvtqFKmeh/tMb7/oj9h8Ik//OMufYFLtdYL/Y8l++uoy3bMFkOK1np+TQX8UzQDibMmG4GZVfblBHphHntVg/1lqxqA2ZXwXQD1nEIpFYs5q+0rrXWe/+FAjjXgOrTW24BtwF+VUi9hThOu6cvehdkiOZM6Y6tHvTXZi9ldNs5f5wljq92vU5Cv5QkvYY5PnUmtLR2l1BPALzDf68sDDPVMDMxWeiDWY74nq7oXszV/KeaPu/o4vwHbNjlJMqe6F3Nmzgal1GOY3SIuzMHy3wA/IsApzIFQSi0H3sH8EOZgdlXMxmzBLNVal/q/JJ5SSmVjTq28A+iNmVhO+BJ4USl1LeaH92rM/twTLaEC//7vUkrtA9pgDkqW1xWj1rpEKTUbmO3/gl2M+b4ZAJyvtX4wiDhrsgCzBddJa30Ysy/aAmyoVm4wZldiVeOBVdocxD8Tu1IqEfPXeDzm6/kA5pfHPcEcax31AOA/p+EuzNmKhzF/AIyh9i/WA8BF/q7TIqBIa+2uWiCQ2OpR72m01mVKqb9jtoByMN/v12COFU4KdD9+wbyWJ+qvd3eZUupvmBMhfgjs9r/mAOXVfzTWsv3vMSdh7MfsCRiLeQ7XG4HU7/8Buq3aPo8BLn/iD+QY/orZ65COOeZ4F+b/+/RAtm+WmnpQqLndMGfM/AWzS6IC80tyOXAT1aaxVtuu+pThOdQ98P8Q5pv6mL+uQ5h9un2rlHFi/got8t9e5vQpzDbgb/79FGLOInuiWjzjMKezVmAmgaswf7U+XqXMacdV5bk7ML9wKjCT1jrgnmDiPMP/w1LgEf/fdwOZ1Z5PwuyqOb/KYwrzy/mHdex7Dt9P0fX4Y1+LeR5UfLDHGuBr3wH4AMjAbBVkAa/gnxpf/b0BpGCOgZRQ9xTmM8V2xnprOdZbqTLwX+X9VNcU5hrfJ9X2E9BrGcLPrq7lNqfa8Wqgaw3b/x/mZ6Lc/3+7AbN72RLI9rXE9Hj1z0AdMbxT5fU7hvkDckKo/6/O5u3ELCYRRvxTWW/SWoftWcBV+SdGzMWcgVTbdNvq21yLOVV0kNba25jxiXOHvzvtKuA8feq0+bOyfaj2EU7kPBnR5LTWKzEHarsFsZkDuE0SjAjSNMwToOv75d7Q7UO1j7AhLZkwdK61ZIQQ5y5JMkIIIRqNdJcJIYRoNGdtCnNRUZE0mYQQ4hwXFxd3yslk0pIRQgjRaCTJCCGEaDRhl2TS0tLqLtTMhfsxhHv8IMfQXMgxNA+NeQxhl2SEEEKED0kyQgghGo0skCmEECGmtaakpASfL5BL7jS9iIgIiorqXEMUAMMwiI6ODnhF8oCSjFJqKvAc5oqqr2qtn6n2/HjgY8wFCwE+0Fo/EVAEQghxjikpKcHhcGC3V79cU/PkcDiIiIgIqKzL5aKkpISYmJiAyteZZPzXDX8Bc7npDGC9UuoTrfWOakVXaq2DulqjEEKci3w+X9gkmGDZ7XbKy+u8SshJgYzJDMNcqnq/1tqFuVpu+F7bQAghxFkTSJJJ4tSrsmX4H6tupFJqs1JqgVKqX0iiq8KnNVvyXCzKqesCgkII0bJlZGQwbdo0hg0bxogRI3jxxRcD2i4p6fuv9lmzZjFixAhmzZrVoFjqXCBTKXUNMEVrfaf//s3AMK31fVXKxAI+bV6971LgOa11atX9VF1Wpj5zso97YOJaJ1YFS0aUEyG5RgjRTEVERJCQkNBk9WdnZ5Odnc3AgQMpKSlh8uTJvPHGG/Tq1euM26WkpLB//34AevTowfbt23E4Tr/6dE5ODhUVFSfvp6Z+/3VffVmZQAb+M4BOVe4nY15x7yRd5fK3Wuv5Sql/KqXaaq1za9ph1YCCMTDtGJvz3OTFdmJ8x8AGqZqjtLS0ev8fNAfhHj/IMTQX5+oxFBUVBTyQ3hi6dOlCly5dADPh9erVi/z8/NNiSk9P56677sLlcjFp0iSUUkRERHD99ddTVlbGZZddxi9/+UuuvPLKU7aLjY2lU6dOBCKQJLMeSFVKdcO8HOv1wA1VC/ivpZ2ttdZKqWGY3XB5AUUQhDGJDjbnuVl5pDKsk4wQomVp9UZmSPdXeFtNIxY1O3jwIFu3bmXIkCGnPffQQw9x++23M2PGDN56662Tj8+dO5ekpCRWrVrV4FjrHJPxX71tJrAQ2An8T2u9XSl1t1Lqbn+xq4FtSqnNwN+B63UjXKhmbAez2bbiSGWody2EEOeckpISbrnlFmbPnk1sbOxpz69bt46rr74agOuuu65RYgjoPBmt9XxgfrXHXqry9/PA86EN7XQjE+1Y0GzMdVPs9hFjkwULhBDNXzAtj1Bxu93ccsstXHPNNVxxxRW1lgv0pMr6Cqtv6RibQd8YH14NXx91NXU4QgjRLGmtmTlzJj179mTmzJm1lhs+fDjz5s0D4L333muUWMIqyQAMiTOXaVh5VLrMhBCiJmvXruXdd99lxYoVjB49mtGjR7No0aLTyj3zzDO8+uqrTJkyhePHj9ewp4YLu7XLLmjlZU6GTcZlhBCiFiNHjqSwsLDOcl27dmXx4sVUVFQQERHBL37xi5PPZWaGZrJC2LVkzovxYTNgS56bwsrwWHxOCCFaqrBLMhEWuCDBjgZWS5eZEEI0a2GXZECmMgshRLgIyyQzxp9kZPBfCCGat7BMMkMT7ERYYEeBh5xyb1OHI4QQohZhmWQcFsWI9tJlJoQQzV1YJhmACR3NJLMkS5KMEEJUVVFRwYQJExg1ahQjRoxg9uzZAW3XGEv9h915MidMSIrgsW+PsySzAq11oy+NIIQQ4cLhcPDJJ58QHR2N2+1m6tSpTJo0iaFDhwa8jzlz5rB3794al/oPRti2ZPrFW2nnNDhS5mNXoaepwxFCiGZDKUV0dDRgrmHmdrtr/CGenp7OpEmTmDJlCk899dTJx6+//npKS0u5+OKL+eCDDxoUS9i2ZJRSXNTRwbv7ylmSVUmfeFtThySEEDWK/tH4kO6v5N/L6izj9XoZN24cBw4c4M477+SCCy44rUyzWOq/OZuQZF5TZmlmRR0lhRCiZbFYLKxatYrt27ezYcMGduzYcVqZZrPUf3N1kX/wf/VRFxUeTYRVxmWEEM1PIC2PxtKqVStGjx7NV199Rd++fU97Xpb6P4N2TgsDWtso92rWHpNZZkIIAZCbm3tygczy8nKWL19e42WuZan/AJycypwpSUYIIQCOHj3K5ZdfzoUXXsiECRMYP348U6dOPa2cLPUfgAlJETy3rYQlWZU80dTBCCFEM9C/f39WrlxZZzlZ6j8AI9rbcVoU2/LdZJfJEjNCCNGchH2ScVgUoxPtACyVs/+FEKJZCfskA3CRfyrzVzKVWQghmpVzIslMTjYH/7/MrMDr000cjRBCiBPOiSTTI85GSoyFgkrN+hxXU4cjhBDC75xIMgBTOpldZosypMtMCCGai3MuyXxxWJKMEEKAuX7ZmDFjAl4yRpb6P4ML2zuItip2FHg4XOKhU/Q5c2hCCFEvL774Ir169aK4uDjobVv8Uv/V2S2Ki5LM/4zFGTKVWQjRsmVmZrJo0SJuvvnmWsvIUv9BmpwcwacHK1h4uJzbe0c1dThCCAFA6ZLTl3RpiKgJX9RZ5uGHH+aJJ544YytGlvoP0uRkc1xmxREX5R6ZyiyEaJm++OILEhISGDRo0BnLyVL/QWofaeH8tja+y3Wz8kglk/2TAYQQoikF0vIIpXXr1rFgwQIWLVpEZWUlxcXF/PjHP+bll18+raws9R+kE62ZhTKVWQjRQv3ud79jx44dbN26lddee42xY8fWmGCazVL/SqmpSqndSqm9SqmHzlBuqFLKq5S6OnQhBmeqv/Wy8HAFWkuXmRBC1KZZLPWvlLIALwCTgAxgvVLqE631jhrK/RFY2BiBBuq8NjbaOw0ySr1szXczsI29KcMRQogmNWbMGMaMGVPjc81lqf9hwF6t9X6ttQuYC0yvodx9wDzgWEgiqydDKS7tbLZmPjskXWZCCNGUVF1dSv6ur6la6zv9928GhmutZ1YpkwT8F5gAvAZ8prV+v+p+ioqKTlaUlpYWsgOoydcFBvdvj6BHpI93BkuiEUKcXRERESQkJDR1GI0mJyeHiorvv1urXto5Li7ulJkEgcwuq2nqQfXM9DfgQa21N5CZCjVdazpQaWlpdW7fxat5dM8R9pYZWNt3o1ts85pEF8gxNGfhHj/IMTQX5+oxFBUVERERPrNbT3SXBSo2NpZOnToFVDaQ7rIMoOrekoGsamUuAOYqpdKBq4F/KqV+EFAEjcBuUSenL392qLypwhBCiBYvkCSzHkhVSnVTStmB64FPqhbQWnfTWnfVWncF3gfu1Vp/FPJogzCtsxOAzw9Kd5kQQjSVOpOM1toDzMScNbYT+J/WertS6m6l1N2NHWB9TUx24LDAumMujpV7mzocIYRokQIarNBazwfmV3vspVrK3trwsBou2mYwvmMECw9XMP9QBbf2krXMhBAtx4ABA4iJicEwDKxWK8uWLatzm6SkpJNTl2fNmsXixYuZNGkSTz75ZL3jaF4j4iE2rbOZZD47WC5JRgjR4nz66ae0adOmXtvKUv8BuKRzBIaC5UcqOe7yNXU4QgjRrMhS/w3UNsLCiHZ21mS7WJRRwdUpkU0dkhCiBZo150ch3d+Tt/67zjJKKWbMmIFSittuu41bb731tDKy1H8IXNHVnGX20QGZyiyEaDkWLlzIihUreP/993nllVdYvXr1aWVkqf8QuKKLk4fXFbE4s4Jit48Y2zmfV4UQzUwgLY9Q69ChAwAJCQlMmzaNjRs3MmrUqNPKyVL/DdQxysKI9nYqvbBA1jITQrQApaWlJ6+IWVpaytKlS+nTp89p5ZrNUv/h7spuZpfZB9JlJoRoAXJycpg6dSqjRo3i4osvZvLkyUycOPG0cs1iqf9zwRVdnDy4roglmRUUVvpo5WgRuVUI0UJ17dq1xjGYmso1h6X+w177SAuj2ttx+WC+rGUmhBBnTYtIMgBXdjOnL38oXWZCCHHWtJgkc3nXCCwKlmZVUlApJ2YKIcTZ0GKSTNsIC2M7OPBo+PSgtGaEEOJsaDFJBmCGf5aZdJkJIRqTYRi4XK6mDqNRuFwuDCPw1NEiZpedcHkXJ7/6upDlRyo5WuYlMdLS1CEJIc5B0dHRlJSUUF4eHj9ojx8/TmxsbEBlDcMgOjo64H23qCQT7zCYnBzB54cqeH9/GTP7xzR1SEKIc5BSipiY8Pl+OXbsWMCXUw5Wi+ouA7iuuznL7N194fELQwghwlmLSzJTOkUQZ1dszXezo8Dd1OEIIcQ5rcUlGYdFMcO/MvO7e8uaOBohhDi3tbgkA3BdD7PL7L39ZXh9uomjEUKIc1eLTDIj2tnpEm0hq8zHqqOVTR2OEEKcs1pkklFKca1/AsBcmQAghBCNpkUmGYDr/Unm0/RyyjyyzIwQQjSGFptkusdZuSDBRolH8+lBuZiZEEI0hhabZABu7BEFwFt7Sps4EiGEODe16CRzVYqTSKti1VEX+497mjocIYQ457ToJBNrN/iB/5yZt9OkNSOEEKHWopMMwM09zQkA/0krwyPnzAghREi1+CQzop2d1Dgr2eU+FmXIBAAhhAilFp9klFLcnGq2Zt7aI8vMCCFEKLX4JANwfY9IrAoWZVRwpMzb1OEIIcQ5I6Ako5SaqpTarZTaq5R6qIbnpyultiilNimlvlVKjQ59qI2nndPCJZ0j8Gp4RxbNFEKIkKkzySilLMALwCVAX+CHSqm+1Yp9BZyntR4E3A68GupAG9stPc1zZubsLpVFM4UQIkQCackMA/ZqrfdrrV3AXGB61QJa6xKt9Ylv5igg7L6lJ3R00DnawqESL4szZQKAEEKEQiBJJgk4XOV+hv+xUyilZiildgGfY7ZmworFUNzR22zNvLpTzpkRQohQUN83QGopoNQ1wBSt9Z3++zcDw7TW99VSfizwmNZ6YtXHi4qKTlaUlpbW0LgbRaEbLvvGiUsrPhxSTrIz7BpkQghx1qWmpp78Oy4uTlV9zhrA9hlApyr3k4Gs2gprrVcopborpdpqrXPrCihYaWlpDdq+LlflFfDO3jK+qkjgqYFxjVJHYx9DYwv3+EGOobmQY2geGvMYAukuWw+kKqW6KaXswPXAJ1ULKKV6KKWU/+/BgB3IC3WwZ8Nd/i6zt9NK5RIAQgjRQHUmGa21B5gJLAR2Av/TWm9XSt2tlLrbX+wqYJtSahPmTLTrdF39cM3U4AQ7g9vaKHRp5u2XC5oJIURDBNJdhtZ6PjC/2mMvVfn7j8AfQxta07mzdxT3rirk1V2l3JQaib+RJoQQIkhyxn8NruwWSWuHweY8N+uOuZo6HCGECFuSZGoQYVXc3sscm3l+W0kTRyOEEOFLkkwt7uwThd2Azw9VyAXNhBCiniTJ1CIx0sLVKZFo4MUd0poRQoj6kCRzBj/tFw2YFzQrqJTpzEIIESxJMmfQr7WNizo6KPNo3tgtS80IIUSwJMnUYWZ/szXz8o4SXN6wPPVHCCGajCSZOkzo6KBPKytHy328v1+uNSOEEMGQJFMHpRQ/9bdmnttagi88FzIQQogmIUkmANemRJIcZWF3kYfPDsq1ZoQQIlCSZAJgtyju97dm/rKlmDBdlk0IIc46STIBurlnFO2c5lIzX2VWNnU4QggRFiTJBMhpVSfPm/nLluImjkYIIcKDJJkg3N47ilZ2xdfZLlYfldaMEELURZJMEGJsBj/p62/NbJbWjBBC1EWSTJDu7htNtFWxJKuStdnSmhFCiDORJAOgNSr3KMaB3Rhp2zAyDkB5zSdexjsM7vaPzfxh4/GzGaUQQoSdgK6MeU7yebFsXodt1RdYdn6HKj21+0srhW6fjLfnADznj8Lb/wKwOwCY2S+aV3aWsPKoi+VZlYzr6GiKIxBCiGavRSYZy87vsP/3eSyH9p18zBfTCh3fFmw2VFkJ6tgRjKOHMY4exrZiPjo6FvfYS3FPmE6rhA7c1z+GpzYe5w8bjzO2Q1u5RLMQQtSgZSUZnxf7B29g++w/KK3xtW6He+IMPMPGo9smQtVE4XFjHN6PZcs6rN+uwHJoL/b5c7F98T88Yy7h3ktv5EWHwTc5Lr7MrGRSckTTHZcQQjRTLSfJeDw4Xp6Nbd0StDJw/eAWXJfdcLIL7DRWG75uvfB164V7+i0Y+3Zi+/IDrGu/wrb8cxJWLeSjIdOZ4riUpzYeZ2KSQ1ozQghRTcsY+NcaxytPmwkmIpKKB57FNeO22hNMDXzd+1D5k99S9vSbuEdOBJ+Xkd/MY8f6B+m6aw2fpJc34gEIIUR4ahFJxvbZf7Ct/QodEUn5A8/i7Tu43vvSiclU3v0o5bNewNslleSKXOZt/xvxLz+Bu6gwhFELIUT4O+eTjGXnd9jnvYZWiop7HsXXvW9I9uvr3pfyx1+i7Mb7KbVEcGnWWmyP3I5ly7qQ7F8IIc4F53aSqazA8dqfUVrjvuJmvIMuDO3+DQu+yVey7P5/sjq2J9El+Tj/8iD2//wDPO7Q1iWEEGHonE4y9o/mYORk4U1OwXXFzY1Wz9jzuvHbKU/wcMr1eAwr9kXzcM7+GSovu9HqFEKIcHDOJhmVl41t0TwAKm//DVhtjVeXUjwxvDV/7nw5Fw+ehSs+Acu+HUTOugvLZuk+E0K0XOdskrF//CbK48Y9fAK+7n0avb7BCXauSXGyOroHP73kWTwDh6NKj+P864PYPvo3yIXOhBAt0DmZZFR2JtaVC9CGgevK285avY8OjsVhgTeyrCy74XEqr74TrRSOD9/A8cLvoVKmOQshWpZzMsnYvvwQ5fPhGTkJndjprNXbJcbKff1iAPjNN8VUXHYjFb94Gu2MwrZ+Gc6n7pNxGiFEi3LuJZmKMmwrFwDgnnzVWa/+l+dFkxxlYWu+mzd2l+I9bwRlj/0TX/skLIf24nz8bqIO7z3rcQkhRFMIKMkopaYqpXYrpfYqpR6q4fkblVJb/Lc1SqnzQh9qYKyrF6PKS/H26I+va8+zXn+k1WD2sDgAntx4nNwKL7pjF8oeexFPvyEYxwvo8dazWFcvOuuxCSHE2VZnklFKWYAXgEuAvsAPlVLVz2g8AIzTWg8EngReDnWggbIt/xwA98QZTRUCl3eJ4KKODopcmic2+K85Ex1Lxa/+iGvSlRg+LxEvz5YJAUKIc14gLZlhwF6t9X6ttQuYC0yvWkBrvUZrXeC/uxZIDm2YgVFZB7Ec3IN2RuEZMropQjDjUIo/jYjDZsBbe8rYkOMyn7BYcd10P4enXI9Whjkh4NVn5MRNIcQ5S+k6fkkrpa4Gpmqt7/TfvxkYrrWeWUv5XwO9T5Q/oaio6GRFaWlpDY27RonLPqbDqs/IO28Uhy6/tVHqCMY/0pLZajIAACAASURBVG28mWGjV5SPOYMqsFZZpDl2zya6fvgKFreL4q69OXD1PXgjIpsuWCGEqKfU1NSTf8fFxZ2yHH0gS/3XtH59jZlJKXURcAdwxmZE1YCClZaWVvP2WhP5yiYAIqdc2aA6QmV2Vx9LPzrG7hJYVJnIzwaYM8/S0tJof9k1VPYdSMT/PUxM+i76/ff/qPjl0+iEDk0cdd1qfQ3CiBxD8yDH0Dw05jEE0l2WAVSdB5wMZFUvpJQaCLwKTNda54UmvMCpI4cwsjPQ0bF4+ww629XXKNpm8NyFrQB4+rvj7CvynPK8r1svyh/7J96krliy0nE+eS/GgV1NEaoQQjSKQJLMeiBVKdVNKWUHrgc+qVpAKdUZ+AC4WWu9J/Rh1s266WsAPANHgGFpihBqNCEpguu7O6nwws/WFFC9e1K3TaT8t//A03cwRlEBztk/x7JxdRNFK4QQoVVnktFae4CZwEJgJ/A/rfV2pdTdSqm7/cUeA9oA/1RKbVJKfdtoEdfCutlMMt5BI8921XWaPSyOthEGq466eCut7PQCUTFU/OqPuEdPRbkqiPj7o9gWf3D2AxVCiBAL6PLLWuv5wPxqj71U5e87gTurb3fWlBZj7NmKtljw9L+gycKoifa5aOXL4aVBuTy3KZvFW7wMSvHgLahA2WJQjjZgjUZZbVTe+SC+dh1xfPA6jrf/jjqWheuH9zSrlpkQQgQjoCTT3Fl2b0H5fHh7DoComCaLQ2sfvuI0fIVb8RZuw1dyAF1xDNCMBEa29xcshIrvqmxo2DEik1BRXXENSsHd+iYi35yLfdH7GLlHqbj7UXBEnP0DEkKIBjpHksxmALy9m2bA31eyH/eRL/EeW4GuzD31SWWgHO1RthhcKpI1OT6Uz0PvOE07aym6Mg+8ZfhKDkDJAbzZS3EDpT90YMuxYs9ai/X5e/Hc8Wdo1aZJjk8IIerr3Egyu8ypy97eZ281G601voLvcB18D1/B980S5UjA0nowRqsBWGJ7oZyJKMO8lo0TcKWX86Ol+UTmK1ZckUCPOBvaU4qv9BC+knR8JfvxFe0wE1cCuBOsQBbG6puwdBiLpeslGPEDUercW3ZOCHHuCf8kU1aCcXAv2mLF26PfWanSV7KfyrRX8RVsNB+wRGBNnIg1cQJGbO8zJoDpXZ1ckuBhQY6Vn6woYOFlCVitUVji+mCJ+/66N9pdbHa5HV2DN2MJPqcXX+Fy3JuWoxwJWBMvxtphIkZkkyyuIIQQAQn7JGNJ24bSPrzd+oLD2ah1aZ8L94G3cR98H/CBNQpb52uwJV2GsgU+FvSb7i62ljnYkOvmL1uKeXBQ7GlllC0Ga8JISBiJTr0X25uP4CnbSnmKBR85uA/OxX1wLkZcX2xJl2NpN/pki0kIIZqLsO9zsezbCYA3tXFbMb6SdMrX34/74P8AjTX5CiJHvoG96/VBJRiAGCu8MDoegD9tKmZtduUZy6sIJ547n8We8APazqsk/otK7O4UsDjxFe2gcscfKV9zC679b+GrPOvnwQohRK3CPskY6bsB8HXr3Wh1eHJWU77hF+jSdJSzIxFD/oKj570o2+ktkECN6+jg/v7ReDXcsayA/ArvmTcwLLhunInrxvuwHYP4/+4gfu8F2Hv8FBXVBe0qwJ3+H8rX/IjKHX/GV5Je79iEECJUwru7TGuMA2aS8Xbr1Qi717gPvot7/xwALO3H4+j9c5QlNNOJZw2JZW22i29yXNyzsoB3JrbBUDUtFfc99+Sr8LVtT8SLT2JfsRgjL5/ynz6Lz7Ufd8YneHO+xnP0KzxHv8LSdgS2LtedMtYjhBBnU1i3ZFRBDsbxAnRUTMgXltRa4973mj/BKOw97sTR98GQJRgAm6F4bXw8reyKhRmVvLCtJKDtvINHU/7wc/ji4rFu30DkH+7D6kskYsAsnCNfx5p0ORh2vLlrqdjwC8o3PoAnb8NpS9oIIURjC+skc7IV07UX1NECCIbWGteef+I+9D4oC45+D2HrfDUqhHWc0CnayotjzPGZxzccZ10d4zMn+FJ6Uz7rn/g6dsGSmY7z9/dgpO/BcCbi6PVTIi/8N7Yu14ElEl/hFio3/5aKjb/Cm78p5McghBC1CeskYzlwYjwmtF1l7gNv4cn8FAwbjgGPYW0/LqT7r+6Szk5m9jPHZ25blk92WR3jM346oQNljz6Pp/cgjKJ8nH+4H8v65QAoezz27rcROeotbCm3gS0WX9EOKjY9RPnGB/AWbmvMQxJCCCDMx2SMzAMA+Dp3D9k+3Rmf4E7/LygDR79HsLYdHrJ9n8nvLohlQ66Lr7Nd3LI0n0+ntsVuqbnl5PF6KCzJpbisgNLKYsqmX4q7nQUO7cX3+V9x7f0S3fM87HYndqsdu7UjUZ1/jbP4O5w5i3EWbKFi468x4gdjT7lZxmyEEI0mzJPMQQB8SV1Dsj9P3re49rwIgL3Xz83zVM4Sm6H490WtueiTHNYdc/HA2kL+NiqesooSMnP3czh3H5m5B8gpyqKoJA+f9p26AwvQzbx2DZWHYOuhWmqKwlDRxFs9tM7ZR5sDs0ho1YWknleT2HEoFiOs3xJCiGYmfL9RXJWoY1loZeBr3/Cz3n1lWVRufwbQ2LrdhK3j5IbHGKR2TgtvTWjFrfO/Y+ee3TyVtY/K0tOThULRKqotcdGtiYqIJdIRjdMRhcWwYj2WhX3dMnxeNxWtWlN+/nAqDEVpxXGOlxVSXF5AeWUpeW4LeW4LaeXA8Vw49BJW9S86tO5MUrueJLXtRtf2vWkVLeulCSHqL2yTjJGdidI+fImdwGZv0L60t4KKrU+Ap8Sc9tv1hhBFGbicwiw27V/D5n1rmK7zQEFlKRiGleS23UhO6E5y2xTaxycTH52AzVr7MasBl+P82yMY6Qfw7Sug4v4n8fUccPJ5l7uSvONHySk6Qk7+frKPrOdIUQ6FHguH8w5yOO/gybLxMQmkJPahW2IfunXoQ2xkfKP+Pwghzi3hm2Sy0gHwdezS4H259r5inmgZmYyj72/O2uKTPp+PXYc3smbHQg5mf39B0bio1hTa+/FpfgolRgoLLkyiZ6vAl4zRHbtQ9ruXiHjhcazbN+B85hdU/ugXeMZdBoDd5qBDmy50aNMFUkYAN+Ary6Joz+tkZn7DEZeVTLeTjMoICopz2FCcw4a0FQAktu5MQmRnIuIUSW1TMIywnjsihGhk4ZtkTozHNDDJeHLX4cn8HJQNR79HUNaoUIR3Rj6fl2/3LGPl1s/JLz4GgMMWQb+uwxiUciFdEnuhteLA0nw+P1TB1Yvz+HJaAu2cQVy8zH+1TfvcF7EvmkfE63/GfWAXlTfeV2PLz4jsSPygR4lNSaP7vjfwFWzEp+GYbkOGcwgHy3ykZ+/maP4hjuYfYmvGKqIiYkhNGkjvzoPpmTTwjK0rIUTLFLZJRmU1PMloVyGuXX8DwJZyC5aYlJDEVhuf9rHtwDd88d27FFfkAxAfncDIvpMZnDoGh63KAp8KXhkXz+ULctmQ6+b6L/P4dGpbomxBtBwsVlw33oevU3ccb/4ftqWfYqSnUXHf79Ft2te8SWwqzvNn483fiGvf6yQW7yWxchFDYzphDLiZQ55WfLN9OdnF6RSW5LJp32o27VuN3eqgV6dB9Os6TBKOEOKksE0yxhFzQNzXsXO991GZ9i+0qwCj1QBsna8MVWg1yso7yKdf/5uM3H0AtI5pz8Xnz6B/1+G1djlFWg3mTmzDxM9y2Jjr5s7lBbw9oTUWI7iTQj1jL8XXqTsRzz+G5cAuIh+7i4p7HsN7hktVW1oPJiJ+EN5jK3Htm4MuO4x3x2w6xfXFkTSZzv1mklOUxa7Dm9hxcD2ZuQfYemAdWw+sO5lw+ncdRs/k87BaZHVoIVqq8EwyWmPkHgHA1y6pXrvw5n+HN3spGHYcfX6JUkF0RQWh0l3OVxs/YO2uxWitiXbG0b/jKKaOuiqg6cIJTgvvT27DpM9yWHC4gp+tKeTvo1rVucZZdb5uvSj7/ctEvPQHrFu/IeLZ3+C68nbc026EWpKcUgbW9uOwJFyIJ2sBrgP/wVe0g7bsoNK3lrbdb2PsgMsYO+AyCopz2H5wPdsOfENm3vcJx2mPon+3YQzqPopOCT0aZdUEIUTzFZ5JpqQIVVGOjoyCqOCW2QfQXheVu58HwNb1BgxnaNc9OyE9ezfzVr5MYUkuSilG9p3MhEFXcvhgRlDno6TG2Zg7sQ0zFubxdloZMTbF7GFxwX9hR8dR8cunsX/8JraP38Qx7zUs+3ZS8eOHz/j/qAwbtuQrsCZOxH1oHpUH34PctZTnfoO1wyRs3W4iPiaB0f0vZXT/SykozmFb+jds2b+WowWHWL97Ket3L6V1TDsGdR/Fed0vpHVMu+BiF0KEpbBMMsYxfyumbf2Sg/vQe+jyTFRkZ2ydrwplaIB5Rv6STR+waut8NJoOrbswY9Qd5myuehrR3sHbF7fm+i/zeHFHKbF2g4fPr8elBgwLrhm34U3pQ8S//oB10xoiZ91JxT2z8KX2P+OmyhqJPeVmDrv6kqy+xpO1AM+RhXiyl2JL/gG2LteibNHExyQwZsBljBlwGUfzD7Fp3xq27P+a/OJjLNn0IUs2fUiXdj05v8do+ncbjsMWukVHhRDNi+Xxxx8/KxVVVlaGpKL8/Hza5WViXb8cX49+eIZPCGp7X2Ueldtmg/YS0f+3GJEdQxHWScfLCnjry7+y9cA6lFKMHXg5V4/5CXFRrU85hjZtgj/JMSXWSu9WNj4+WM6qoy6ibIrh7Rz1ilMnJuMZdhGWtO1Ysg5iXfUFWCz4UvtBHVO48wpKadfzEqztx6NdheiS/fiKtuPOmo9SFozoHijD7H6MdsbRI6k/I/tOoUu7VDSa/OJj5BcfY9fh71i780sKinOIcsYSGxl/1rrT6vsaNCdyDM2DHMOpIiIifl/1fli2ZFSOvyVTj+X93QfeBl8llrYXYokfGNK40rN38+7SFyipKCI2sjXXjruHLu17hrSOK7o6eX5UK+5dVcis9cexKsU9/aLrtS+d0IHyR/+Bfd5r2OfPxfH+q1h2bKTyJ79Ft6r7DWdEJhHR/xG8x6/Ctfc1fIVbcO19Bffhj7Gl3II18aKTY12GYdAjaQA9kgZQ6S5ne/q3bExbwcFje9iQtpwNactp1yqZC3qO47zuFxLpqN8xCSGal7BMMkbOUQB028SgtvOVHsSTtRCUgb37bSGNacOe5Xzy9b/xaS/dEntz7bifEu2s/5Uzz+SG1CjKvZpffV3Ew98U4fFp7hsQ/NgUAFYbruvuxtt3MI6Xn8a6YyPGo3dQedfDeM8LbHFQS2wvIs7/I978b3HtfQ1dmo5r57O4D72PvfvtWNoMPaWF4rA5GZw6hsGpY8gpzGJD2go27VvFscIM5n/zHxZ9+z/6dBnCBT3H0TWxN8ZZOjlWCBF6YZlk6tuSce17A/Bh7TgNI6pTSGLRWrNk04cs2/wxABf2m8rkIddiMRpnttoJd/SOxqoUP19TyKxvj+Pywa/Oq2eiAbwDhlH+5Ks4Xp5trhLw1wdxTb4K1zU/BnvdXXJKKaxthmJpPRjP0aW497+JLk2ncstjGK0GYu9xB5bY0y/JkNCqI1OHXs/EwVez+/B3bEhbzt7MbWw9sJatB9YSH5PAkNRxDO4xhpjIVvU+PiFE0wjLJGMU5ACgWwc+Q8lbnIY3dy0YDuzdbgxJHF6fh0/WzGHj3pUopZg2/BaG9Q5ujKghftQrCqsBM1cV8uTG47h8mocGxdR7XEO3akPFr/+Mbf472D94HfuieVi3rqfix4/gS+kd0D6UsmDrMBFru7F4Mj/FlT4XX+EWKr79GZaE0di734YRefq0c6vFSr+uQ+nXdSiFJXls3LuCjWkrKSjO4cuN77Pkuw/omTyIob3G06PjAFnORogwEZZJRhXmAeCLD3ygyn3gHQCsSdNQ9oYv8ujxunl32QvsOvwdNouda8ffS+9O5zd4v8G6MTUKq6G4Z2UBf9xUTH6lj2eGxQV9wuZJhoF72o14+11AxMuzMbIO4nzyXtxX3Izr8pvBGthbRlns2DpfhbXDFNyH/of78Ed4c1ZRnrsGa8dLzanjjtY1btsqug0TBs1g/MDp7Mvaxrdpy9l16Dt2Hd7IrsMbiYtqzeAeYxmcOlZWiRaimQu7JKPcLlRZCdpihei4gLbxlezHm7sGDHtIpiy7PS7mLnuePRmbcdqjuHnSr+iUELoLpwXruu6RRFgUdy3P55WdpeSU+/jX2HgctVz0LBAnTt60z3sN28L3sH/0byybvqbix48EtR9li8be/XasyVfg3v82niOL8GR+hufIYmzJ07B1vgZlr7kbzDAMUpMHkpo8kJLyIjbuXcWGPcvJL85m6eaPWLb5Y1KTBzIkdRy9Op0n18IRohkKu0+lraQIAB3XGgLsFnKl+1sxHS+t9ddzoNweF+8s/TtpmVuJdERz6+QHGnT+S6hM7+qktaMtN36Vx0fp5eRWePnPxW2IszegW8nuwPXDe/GcfyERrzyNJX0Pkb+7i3ZjLoeUbmAJ/O1jONri6PNzbJ2vxLVvDt7cNbgPzcOd+Tm25OnYOl+FstU+USLaGcfYAZcxuv8lpB/dxbd7lrHj4Ab2ZGxmT8Zmop1xDE4dy5DUsXKipxDNSEDfQEqpqUqp3UqpvUqph2p4vrdS6mulVKVS6tehD/N7J5NMAFNsAXxlmXiPrQJlw9bl6gbV7fF6qiSYGG6b8mCzSDAnjOng4PNLE0h0Gqw66uLS+TlklHgavF9f70GUPfU67rGXotxukpZ8gPP392Ck76l742qMqM5EDHyMiAv+gaXNMPBW4D74LmVrbsW1/y20u+TM2yuDlA59uXbcvfzm2r8xdegPaRvXgZLyIlZs+ZT/m/cb5iz8E1sPrMPjbfixCyEaps6foso80eEFYBKQAaxXSn2itd5RpVg+cD/wg0aJsgrrySQTWIvEffgjQGNNnIDhaFvven3axwerXj6ZYG6f+hDt4xt+Rc5QG9DaxsLLErhqUR7bCzxM+CyH/0xow9B2DVwV2RlF5R0P4Bl2EcYrz+A4mIbz8btxT70G14zbwBHcWfuW2FQs5z2Bt2gX7gNv4s3fiDv9P7gzPsbW6UpsnabXedmFqIgYRvWbyoV9p3DoWBrr9yxle/p69h3Zzr4j24mKiGFQ99Fc0HM8beOCm+4uhAiNQFoyw4C9Wuv9WmsXMBeYXrWA1vqY1no94G6EGE/xfXdZ3S0Z7S7Gc2SRuV2nGfWuU2vN/HVv+1cYjuCWSb9qlgnmhC4xVhZPS2BsBwfHyn1M+yKHuXvLQrJv74Ch7PrJ73FNvRYA+4J3ifzt7Vi2f1uv/VniehMxaDYRg5/FiB8EnhLcB940Wzbpc9Ge0jr3oZSiS/ueXD3mJzxw7XNcNvwm2scnU1pRzOrtC3juwwd5bcHTbN63BrfHVa84hRD1E0inehJwuMr9DCCws/QawYkk4wugu8yd9QX4KjHiB2NEd613nUs3f8S6XV9hMazcePHPSWrbrd77OlviHQbzJrfh4XVFvLqrlLtXFrCr0M2swbH1n3nm5zsxVjN8Ao7X/4zl8D6cf/o17hEX47rubnTrhKD3aWnVH+f5z+At2Ixr/5vmMjX75+A+9B62pMuxdZqBstc90cPpiGJEn0kM7z2RjNz9bNizjC0H1pKevYv07F04173NoO6jSIjoBqTW4+iFEMFQWuszF1DqGmCK1vpO//2bgWFa6/tqKPs4UKK1frb6c0VFRScrSktLq3fAnT+dQ5vNqzl06c3kDR5be0Htpf2Rx7F4C8lrezeVzn71qm/fsS2sTvsEhWJc76vo3Caw80Wak/eOWPnLPhteFMNbeXmiZyWtQ3VNMa+H9l8vJHHV5xgeN16bg6NjppEzfKI5A7A+tMZeuZuY4wtxVO4FwKfslEWNoiT2YnyWwGYVnuDyVJKeu420o5vIKz1y8vGEmGRS2w+iS9u+2CxykTUh6is19fsfbHFxcaf8ig0kyYwEHtdaT/HffxhAa/10DWUfJ4Ak0xC+J2YSu28b5b+YjXfQhbWW8xxbReW2p1CRyTiHv4yqx9IkB7P38MbCP+L1eZg2/GaG95nYkNBPSktLO+VFORuWZ1Vyx/J8cit8dIg0eH18a0a2r9/imjXFr3KO4Hjnn1g3rATAl9iJypvuxztgaIPi9hZux31wLt689f6KbOblBbpcU69LNGTlpfPtnmVs2rsat9fsOnPYnJyXMpIhPcfTsRlN5KhLU7yPQk2OoXkI5TFUTzKB/NRcD6QqpboBmcD1wA0hiaYeLGXFAOiYMy8x4slaAIAtaVq9EkxBcQ7/XfJ3vD4Pw3tPDFmCaSrjOjpYcUU77liez9fZLqYtyOXxIbHM7B8dkpWPdUIHKu5/EsvWb3C8/Q+Mo4dxPvsbPEPGUHndT9Dt6zeGZWnVD0urJ/EW78WdPhdvzmo8WfPxHPkCS7vx2DpfHdRlszu26coVI2+lR/xQyi15bNizjMM5+/hm9xK+2b2EpDbdGNJzHANTRpx6OWwhRL3UmWS01h6l1ExgIWABXtdab1dK3e1//iWlVCLwLRAL+JRSPwf6aq2PhzzgcnMAW0fXfk6Fr/wo3vyNYNiwJl4cdB0VrnLe+uqvlFUW06PjAC4Z1mQ5NaQ6Rln4ZGpbntxwnL9vK2HWt8dZfqSS50fHkxgZmrXWvAOGUfaH180TOD9+E+uGlVg2rcE9YTqu6bdAHT8OamOJ6YFlwKP4Sg/jPvgunuwleP03I34Qts5XYWl9QcAJ02ax09d/Xs3RgsNs2LOcTftWk5l3gMyvD/DF+ncY0G04F/QcT1LbFLmipxD1FFCnudZ6PjC/2mMvVfn7KHBWpltZKszZRjqq9iTjyfoC0FgSxqBswS0aqbXmg1WvkFOYRUKrjlw3/t5GX+zybLIZiieGxjGivZ2frirgy8xKRn6Uzd8ujGd61xD9crfacF92A56Rk7DPew3r6oXYF3+AbdVCXNNuwD356oAW3ayJEdUJR99fY+t2E+6Mj/FkfYGvYBOVBZtQUZ2xdZqBtf3FqCDGWBLjO3HZ8JuYPORath9cz7d7lnEwew8b0lawIW0F7VolmVf0TBlJbFTDTuYVoqUJr4uW+bw45r2OAlzX3FnjxbW0z4tr57PgLcfR814MZ/ugqli9fQFrdy4mwhbJ7VMfIjay4eucVdccLnKUGmfjuu6R7Cxws7PQw0fp5aQXexjTwUFEHcvRBBy/MwrvkNF4B49B5R7BkpmOdcdGrKsXoaNj8CV3q/MCabVRtmisbS7AljQNbDHossPo8iN4c9fhzlpgziqM6oyy1Hz+Tk3HYDEsJLbuzODUsQzoOhyrxUbe8WwKS3PZd2Q7X+9YRHr2bgBax7TDWt+JDSHSHN5HDSXH0Dw05kXLwivJlJbgmP8OOjIa9+U31VjEm7sOz5EvUJHJ2HvcGVQ3x4EjO5m36mUArh1/L53bNc5gXnN5U8bYDK7t7qRthMHKI5VsyvPw3r5yusdZ6BFnq3W7YOPXca3xXDgZb2p/jIx9WI5mYN24Cus3S9FRcfiSutQ/2VjsWFr1w5p8OUZkMro8G11xBF/hFtwZH6PLslCONqediFvXMURFxNAjaQAj+04mqW03vD6v/4qe2ew8tJGvdy4ip/AIdquDVtEJTdKd1lzeRw0hx9A8yJUx/VSpOcRzxq6yo18BYO0wJagP/vGyAv63/EW01owdMI0+nQc3LNgwoZTirj7RjOvg4O6VBWzMdXP9l/nM6OrkmeFxtA/RWA2At/8FlPd9GeuaL7F/NAfjyGEiXnoS76dv4ZpxK94hY6GeS/gr//ibpf0EM8Ecmoc3bz2eo1/iOfolRkwq1uTLsbYbh7IE3lVntVjp03kwfToPpryylG3p37Bp32oOHUtj8/41bN6/hhhnKwamjOC87qNIjO8k4zdCVBFeSabEP7MsquZL82p3Cd7cdYDCmnhRwPv1+rz8b9k/KakoIqVDXyacf2Uowg0rPVvZWHRZAv/aWcofNh7nw/RylmRV8OTQOG5KjcQI1RenYcEzegqeERdjXfUF9k/ewpKZjvP5x/F26o5r+i14h4ypf7JRCkv8eVjiz8NXloUn63PcWYvwFafh2vlXXGmvYOs4BYunH8GejOl0RDG010UM7XUR+cXH2LxvDZv2rSG/OJvV279g9fYvSIjrSP9uwxjQdTgJrTrW6xiEOJeEVXeZkXkA25rF+JK64Rk1+bTnPdlL8OauwYgfhD35ioD3u2LLZ3y3bxUxka24dfIDRNgbd+pqc21eG0oxrJ2dq1Oc7C3ysLPQw4LDFXyVWUG/eBsdo8xWTUjiNwx8XXvivvgH+OLbYhzaiyU7E9s3y7B+swRtj8DXsQtY6t+SUrYYLK2HYEuejnJ2RFfmml1pRTuIKlmBr2gHKBsqsgPmEn2Bczqi6JbYmxF9JtIjaQBWw0pBSQ5FpXmkH93Ful1fsePgBspdpcQ6W+F01PzDqCGa6/soGHIMzYOMyfhZ9u/C+u0KvCm98A4df9rzrrRX0BXZ2Lr+EEtMj4D2mZl7gHkrX0ajueGi+0lsHZrLMp9Jc39TtnIYXJPipEeslW+OuUgr8vLmnjIOlni5IMGOq7ggdPEbFnzdeuOeMB0d1xojKx3jWBbW71ZjXbEAtA9fcgrYah8jqosyrFhiumNLuhRLm6Hg8+ArOYQuz8SbsxJ3xmfoyjwMR5ugL2inlCIuqjW9Og3iwn5T6NyuJ4ZhUFicS1FpHvuP7GDtzsXsPryJSnc5sZHxRNgj630sVTX391Eg5BiaB0kyfpbdm7FuWYev13l4B4085TlfRQ7utJfAsOHo80uUUfcUVrfHxZuL/0JpxXFzvauzdMJlOLwplVL0a23j1l5RaA0bc11sznMzaeD6vwAAIABJREFUZ3cpWmtGdorD1sA10E5hseLr3gf3xTPwJSajsjOwHMvCuv1bbEs+RpWX4uvQCZxnXpm5LoajLdaECzlU2Zf/b+/Mo+M47jv/qT7mHgzuGwRIALxJUSd1y5ZlSba8lhPLsbKxE+e0s95NvLvJJt68t2/fy0vWm93Em/eStZO1s3ZsxV7FUmLJR3RYJ2WRosSbAC+AJADiHAAzwGDO7q79oxsXTwAECI5Zn/fqVXV3dfevZ3rm21Vd9ftV1K1H5saQ2UGcieNY537odbdKtFDDgu6heccWGhUlNWxacyt3b3mExip3fs34ZNwdodbvjlA7de4w6VyKSLCE0FW0cIrhProS6hquD5TIeOiH92J07sfechv2llvnbSuc+xHO+D70yrsx6xYmFi/s/S7H+w5QGavjk+/7/DWLrFhMN6VfF7yvPsDH14XoSdl0JCz2JnW+fTJNUBdsLTev2uHmPDQNZ00r1oOPY6/diDY2jDbYh37iMOZLz6D1nUbGypEV1QsOWncxRscnqVp7N2bDY+iVO0HoOJlzyOwg9ug7FHr/GSfV7XanBWsX3Z2maTqVsTq2NN/G3Zsfoa6iBXA9SYyn3CHRuztf5siZd5hMJ/CbQaLB0kUNGiim++hSqGu4PlCjyzxELgOADFzY3WAPvQaAUfvggo7V1X+UtztfRBM6T9z3WXzG0iYH3iisKzH4zkMVvNaf5YtvDdOZgt/bneSvjqb44s0lPLE2uLxiIwT2jrvI7LgL7eQRzBefwXj3dcy9r2HufQ17TSuFD34c684PLHli5zR6tB19Qzu+tt/EHnmLQv8LOImD2MNvYg+/CUYEo/o+jJoH0Uq3LNpNkWn42NJ8G1uabyNfyHHy3CE6et7jRO9BRhL9vJ7o5/VDzxMLl7Npza1sWnMrzTXrf6YmAStuXIpKZMhl3fy8AFlOuh8n1QV6CL3itiseJpOb4tldXwPgfTseLwrX/dcL76sP8M2bcnT6GvmTfROcSFp89o1x/uLgJF/YHuWJdcHl7UYDnPat5Nq3kh8bwXz1OYxXn0fv6UL/+p8hv/sVCnc/jHX/h3HWtF7VeYTux6h9EKP2QZzsCNbQa9hDr+CkTmP1/xir/8cIfzVGzfvQax5Aiyze3YzP9LOl5Xa2tNyO7VicHjxGZ88+OnveIzk1xu7Ol9jd+RJBf5i2+m1saLyJtoZthAOL81yhUFwvFJXICE9k5HkiY8d/CoBeuXNB/eg/2PMtJtJjNFa2cv+2jyy/oT/jCAGPtwR5bE2A/9eV5ksHJjmetPjtN8f50/0T/M7WCJ9qDxM0lldsZHkV+Y//Ovl/9SmMd17DfPlZ9NPH8b30DL6XnsFeu5HCAx92WzdX++4mUIWv+RPQ/AlXZAZfxRp6FZkbptDzNIWepxGBOozqe9Cr7kUr2bBowdE1g7b6rbTVb+WxnZ/iXPw0nT3v0Xl2H/GJAQ6f3s3h07sRCBqrWlnfuJ31jTdRV96s5uIoioaiEplLtWSs4bcAMKruueIhDp/ew6HutzENHx+/77dUl8RVYGiCX2oP8wutIf6xK82XD6c4mbT4/d1J/uzAJJ/bHOEzG0JUBJb5M/b53bk29z6CdvYkxus/xHz7JfTTx9BPH0P+w//GuuMBrHsewd54E1zld6xF1uJrW4vZ+hmcxFGsoVexRn6KzA5Q6PkehZ7vIfyV6FX3YFTfhxbbtPh3OEKjqaqVpqpWHr71F4gnBznRd5AT5w5yZvA4vSOn6B05xU/2P0s0WEp743aiWhWNa+oJ+q9OUBWKlaSoRGbmnYx/dh6LkxvFmegEzXfFrrKJ9DjPv/1NAB697UkV932ZMDXBv24P88nWED/oyfIXhyY5OFrgj/dN8GcHJ3hiXYjf2hTmporlDwzmNLeT/+UvkH/ytzH2vo7xxo8wjh3A3PUC5q4XcEorse58EOvOD+C0rL+qwQJCaOhl29DLtuHb8HmcRAfWyC7skbeQuThW3/ex+r6P8JWhV9yBXrkTvexmhLH4eVeVsVoqY7XcveURcoUs3QMdruj0HWIiPca+k28A8MaJZ2moWMe6+s201m1hTXUbhr704d4KxXJTlCIztyVjj3hdZeW3XtIZIrjelf9p19fJ5Kdob9jG7RsWNkBAsXB0TfB4S5CPNgd4tT/HVztSvNiX46mTaZ46meauGh+/uTHMY81B/FdwwrlofH6sex7GuudhxFAf5q4XMN7+CdpIP75/eRrfvzyNU9tE4c4PYN21+PAP5yOEPiM4sv2zOBPHsUd2YQ2/hcwOYg28gDXwAggTvWz7jOhowcU/2PjNwIxrGyklQ+O9nOg7xKFTe4inztEX76Iv3sUbh57H1H0016yntX4L6+q2UFvehLZEv3AKxXJQVCLDzDuZ2SdDa8TtKtOv0FX2zvFXONV/mKA/zMfu+XXVp72CCCF4sCHAgw0BupIW/+dYin84mebtoTxvD+Up9yf5hdYgn2oPs7V8+Z+6ZU2j++7m538NrbsT4+2XMfa8ijbYi/+fv4H/n7/Bxqp69Ls+gHXrfTjN7VffwoltQo9twmz9DZzUaezRd7Dje3AmjmGPvYc99h6c/AoivAaj4g708lvQYlsW5UfNPZegtnwNteVrqAtuYE1LE2eHjtM10EFX/1GGxns51X+EU/1HAAj5I7TUbqSlZgPNNRuoLWtCW6LLHoViKRSVyIjz3snIwgRO4hAIHaNy5yX3iycHeWHvdwH46F2fWRH3/YqL0xoz+NLOUv7olhK+eyrNN0+kOTJW4KsdU3y1Y4odFSafXh/i51qClC/3uxshcFo3k2/dTP4X/w16x35XcPbvIjjSD899C99z38KprMG65V5XcNZvu6p3OEII9Og6N1pny5PIfAJr9F1XdEbfRU71UJjqodDzPdBMtNhW9PKb0ctvRou0Lnp4tN8MsL7xJtY33gRAKjNB90AH3QNHOdV/lOTUKB1n36Xj7Lte/SBrqttpqdlAS+0G6ivWrnrIAsXPNsV1d503usyKvwPSQSu7+ZLByWzH5pk3/4aCneemdXezteWOa2auYpaoqfGbmyL8xsYwB0cLPHUyzdPdaQ6MFjjwdpI/2J3k/fV+fm5tkMeag8R8y/y0rRvY227H3nY7OavA4Ms/pHmwG33fLrT4EL4Xn8H34jPIcAnW1tuwt+/E3nY7MnZ1QcqErxSz7iHMuoeQjoWTPIo9uhd7bD9OqgtnfD/O+H4KXYBZ4jn3vBm9dBsi1LjoFnckWML2dXeyfd2dSCkZmxzmzOAxzg6d4MzQccZTI5w8d4iT5w4BYOgmTVWtrKleT1NVK41VrWq4tGJZKSqREdn5L/7t0b0Al23FvHHoefri3ZSEynls58Vj0CiuHUIIdlT62FHp449vj/GDngzfOZXmtf4cL51zk++nCT7YGODn1wb5YGOAkuUWHMNksnULuUc/Br/8BbdL7b03Md57E23oHOaeVzD3vAKA3bwee/sdWNt34rRugqt46heaMeMhGkDmE9jjB7HH9mOP70Nmh2cngAKYpeilW9FLt6KVbvHm5Sy8lSWEoKKkhoqSGm5d/wAAyakxzg4d58zQcc4OnWA4cY7Tg8c4PXhsZr/yaDWN3ki3xqpWasvWqNaOYskU152Tn+0uk47t9nMDesXFWyd98W5eO/h9AH7+3t9QQz2vMwKG4Il1IZ5YFyKetXn+TJZnTqd5azDPD3uy/LAni6nBvbV+Hm0K8GhTgOboMt+ymobTtoV82xbyn/wcYrAP49Ae9MPvoHfuRz97Av3sCXzPfxsZDGOv34a96WbsTTe7kz+vpmvNV4pR8wBGzQNIKd3InuP73VZO8igyP449sgt7ZJe7gx5Cj21GK92KHtuEcBZ/7li4nO3r7mL7Otf331R2krNDJ+gdOUXfSBfn4qe94GzDHOp+GwBDM6mrWENjVSsNFWupq2ihsqRWvdtRLIjiEhnLcnPDwJk4BlYKEWxAC10YtyNv5Xjmjb/FkQ53bX6Y1vot19hYxWKoDOj86sYwv7oxzEDa5vtnMnz/TIY9w3le7c/xan+OP9iTZHOZwYea3EEFt1f58C3zKDVZ20ihtpHCwx+HfA792AH0w+9gHHoHbbAX4+BujIO73bqhCPaGm7A37cDeuAOnad2SRUcIgQjVo4XqMRse80SnHztxGCdxBDtxxPWrNvYu9ti7FIBaBOlEM3rJBrSSjeixDYhw86JaO+FAlM3Nt7K52fUFaDs2Q+N99I100TfSRW+8i3hygN6RLnpHumb2Mw0ftWVrqK9opq6ihbryZqpLG1SLR3EBxXNHOA5COm5Z02e6yi41N+bFd58mPjFAVayeD97yiWtlpWIZqAvpfG5zhM9tjjCWtXnpXI4f97hxbTrGLTrGU/z5oRRhQ3BPrY8H6gO8r87P5jJjeUcN+vzuu5ntO8n/EoixYfTOA67wdO5HGxlwQxLsd0c4ykAQe90mnNbN2G1bsNs2QyS2pFO7otOAFmqA+kcB19O4kziCnTzqjlqb7EZOncGaOgMDL7g76gG0aLsrPNE2dzBBqH7BwqNrOvUVzdRXNHPHRneYfyY3xbl4N73xbgZGz9A/epbk1OjMBNG5+1aXNlJf0UxNWRPVpQ3UlDUSDpSo0Zw3MEUlMgBSaK7zxBmRuf2Cqqf6j7Dn2Muu88v7P4tpLP8kQMW1oTyg88nWEJ9sDZG3JW8N5nihL8vr/Tk6ExYv9uV4sS8HQHVQ495aP3dW+7izxseWsuX1EC3Lq2fm4gCI+OCM4OjHD7mi07EPOvbN7OPUNmG3bcZetxmnZb3b2lmiQ08tUIVW+/6ZqK8njx9lXa3ASR7DnjiOM3HcDVuQOIyTODxnRz9aZC1atBUtsg4t0ooWabnsvLK5BP1h2hq20dawbWZdOptiYOws/Z7oDIydZXRikIExtzyXkD9KTVkD1aWN1JQ1Ul3aQHVpg+q+vkEoIpGxAZCa5s7yT3WB5kcv3T6vWjqX4p8855fv3/Ex6j0X64rix6cL3t8Q4P0N7p/jYNrmtf4cr/Vnea0/x2DG4dnTGZ497Q4QiZqC26tcwdlZ7eeWKpOouXzvEWRlLda9j2Ld67Y0RGIUrasD/VQH+qmjaKePoQ32og32Yu5yWxpS03DqW3Ba2nGa12O3tOOsaYOLeBa/IpoPPdaOHtvM9GwjmU94gnPCHb022Y3MDeNMHHO7mGd3dltK4eY5aY3b6lmA/79QIEJr/ZZ53dDZfIbB8R4Gx3oYGu9jONHH0Pg50rnJCwYXAJSEygmbMY6PrqOipJbKWB2VJbXEwhXqfc/PEEUoMjr2qDvmXy+7CaHP/iCklDz3028wkR6nqaqN+7Y9tiqmKq4NtSGdJ9tCPNkWQkrJiaTlTfjMsWc4z5lJm1f6c7zSnwMmEUB7zGBHhUmDNHi4JMe2cpPIMgmPLK3AvvU+7Fvvc1dYFlrvKfRTHWjdx9DOnkDr70Hv60bv64Zp4RECWduI07AWp2EtduNanIYWZE0jGIv7iQpfqTvacs6IS1mYwJnsxkl1YU924aS6kekeZLoXO907O7AAQGjue87wGi81I0JNaMH6K7rHCfiC7vybmg2z55aSifSYJzrnZsRnONHPRHqMCcYYSJ6edxxDMykvqaEyVktFSQ2VJXVUxmopj1arrrcipHhExp5tyVyqq2z/qV0cPbsXnxHgifs/q5xf3kAIIdhQarKh1I3mCTCQttkzlGf3sCs6R8YKnEhanEhagI8vn44jgA2lBtvLTTaVmWwqM9hcZtIU1q/+z8wwcNZuxFm7cXZdLovW24V29hT62RNoZ06i9XWjDfSiDfTCu2/MVJW6gVPbiNPoio9T34ysacSpabjASexlPxuzBL18B3r5jtkWj53HmTqLM3UWme7Bmepxy5nBOeLz1vzj+MpmBtqIYL0bPTTYgBasu6QAueGpK4iFK2YmjAI4jsN4aphDnfvwhQXxiUHiyUFGJwaZzCQ8Ieq74Him4aM0Ukl5pJqyaBVlkSo398p+c+Gfi+LaUDwiM92S0TXs8QOA669smrGJIX6459sAfOTOT1Merb72NiquK+pCOh9bG+Rja90/wJwt6RwvsD9e4PXuON2FIB3jBY4lLI4lLCAzs2/UFGwsNVzhKTXZUGqwrsSgKaxf3XsefwCnbQtO2xas6XVWAe3cGTf1nZ4pi/gA+rkz6OfOAK/OO4xTVklbSQX+lnZXiGoaXAGqrl/QOx+h+9BL2tFL2uetl3YOJ92H9ETHmepxlzMDyPw4Mj+Okzxy4fF85YhgHSJQgxaodvNgDSJQjfBXz+txANA0jYqSWhrL22lvn29DrpCZEZx4cpD4hFsenxwhk59iJNHPSKL/otcVDkQpi1R5wlZOLFxOSbicWMjNo8FS1RV3jSkakRFeS8YqE+7Q5UDNzNBl27H4xzf+hryVZWvLHexovbLLf8WNh1+fnQh6r95Pe3szWUvSMV7gyHiBjvECneMWnYkCwxmHvSMF9o4U5h3Dp8HaqCs4bTGD1hI3rSsxqAtpaEtp/RgmTnO760NtLrkMWv/ZWQHqP4s21IcYHkAbjxMdj8PZ4/N2kUIgYxXIylqcyhpkRY2Xe8uVNeC/dLeX0P3o0VaIzg8AJ6WNzMaRmXM4mQGc9Dlkph8n3e8J0BgyPwbJo9gXO66vzBWcGRGqRvgrMHNpnGypu90Lf+43gzRUrr1oMMFsPs345AhjqRESXu6GtB4hMRlnKjvJVHaSvnj3Ra9PExrRUCklofJZEQqVEw2VEg3GiARLiQRj+M2A6pZbJopGZKZbMoUad3F61jTAy/ueoS/eRUmonI/e9Rl1cygWTMAQ3FLl45aq+U/a8axNx7hFpyc+pyYsuicsBtIOx5MWx5MW9M4/lqlBQ1inKazTFDFoiug0RXTWRAzWRHTqQ/ri5vX4gxd2twHYFmJ0mIF9e2jSJWLoHNpQnytAIwNoiTgk4uinLmxxAMhoDKeiBllejVNWiSytQJZWIsvc3CmrgHDJPKehQuiIYA0Eazi/E3pWgAZwssPI7BAyO+SVh5G5kZlWEBPH54lQFZAZBhAIX6nbIvJXzCZfBcJf7iazFL8Zo66imbqK5guuy5EOqUyS8ckRklOjTKTHSU6NueUpt5zKJr11Y/SOXPqjN3UfkWBsJkWn85ArQpFACeFACaFABCnl5b7FG57iERmvJVOodr9QvexmAI6e2cuuIz9CExqfuP+zalikYlmoDOjcX6dzf938rqdUwaF7wqJ7wubUhEWXJz7dExYjWYczkzZnJm0gf8ExBVAV1KgJ6tSFNGpDulfWqQlq1IV0akM61UEN43JdcrqBrK5nsnUrhfO6mrAtxHgcER9Eiw8hRofQ4oPu8ugQYnQYMZlEn0zCmROXPIU0TE983OTEypHRUmRJKTJaCtEYMlqKU1KGiEQvKUDgiVBuzBOf4VnxyY+SSfbjFylkPjErRKmuixxl7vWHEL4Ywix1hcmMeQJVSsiMEfbFaKquAbMdYURA8888eFp2gcl0guTUKElPhCamRpnMJEl5aTKToGDlGU+5LaQroQmd8IEoIX+UcCBKyB8h5OXzl91y0B/CZ9w4LaXiERnHRmpgVbhio5VtZzjRz7PecOVHbvskLbUbL3cEheKqiZga2yt8bK+4cFvGkvRNWfSmbHpSNr0pt9w7ZdObsulP2wxnHIYzDofHLn0OAZT7NSoCXvJrVE6XA/rMunRK4E9ZVPg1QoZw/7R0A1lZ63aXXezgjoOYGEfEBxHjI2jjo4hEHDGdJ8bQEnFEOoWID0J88IqfiRQCIiWuCEVLkdGYm0JRZNhLoQiEI8hQAzK8EVkZgVCY3q5u2tvbkY6NLIwjc6OuIOXibjnvlfMJNxWSYKeRmTQyM7CwL00zXbExowgjQtCMEjKi1JlRRCSCKCtFGE1gRhBGBGGEyUmdqUKBVC5LKjvhic+sEKUyCaayKdK5SQpWnsl0gsl0YmH24Hbb+X1Bgr4wAV9oJgXnledvc1MQnxHAZwaKZmBT8YiMbVOoFGCACK9hytb51st/PvMe5q7Nj6y2hYobnKAhaI+ZtMcuHiPHciTDGYfBtM1gxmYw7Xj5dHKXRzIOozk3kbzsGeHAEOB21ZWYGjGfIObXiPm8sk+bk6aXw5RE2wiXtRMxNcKGIGwKQoaYfaeUyyISo+7cn0QckRxHTCa8lERMJGaXUxMwmURMJoGzlzb3PKQQbPMF0KIxZCiCDEchGEIGQhBwcxmohOAarxxCRgNIv0AaNo5RwBE5JBmkk0IWksh80hOjBBQmkdYkOAW3hZQfZzEdW0EgiKBKDyKMkNuC8ocgFEIYIYTRCHqI4cQUgVgVGUcjYwvSlk3GskhbBdL5POlCjkw+SzqfJp1Lk8lPUbDyZHJTZHJTi7BoPoZu4jcDM6IzXfZPl73kNwL4TD9+M+iuM/z4DD+m4cc0fESDS/NKsWA7F1JJCPEo8JeADnxNSvml87YLb/uHgTTwGSnlvgsOdBUIxyZf544KyUe38tTLf04iFaehcq0KQqYoCgxNUB/WqQ9f/gnUciRjOYd41mE063hlm9Gsu256W38yQ0oajOYccjaeMAGTF3v1vjCmBcfNTSJGHWGznrAhiFRrhBsEIV0QMARBXeDXBSHNIZabJJadoCSXJJKZIJyZIJCbIpBL4cu4ychMYmSm0NKTiPQkIj2FkctALnNlw66A1A3Xca7PDz4/0hcAXwzpq0YGfMiAhh3QkQGBY4L0gWM6SMNG6hZSs3BEHkkOyCNlDilzIPNuy8lOu+e5yLnLAVJwyem0AvB7KQoIHVv4yYsAWXzk8JGTPrJSJyd1so5GzhFkbUnWluRsh6ztClfett1kWVh2AcsuMMXkVX12D93yBPXBlesFuqLICNfp0V8DHwT6gL1CiOeklB1zqn0IaPfSTuArXr582Db5Wp2kpfFM52mGJ8eoKKnh0w/9BzU2XvEzhaEJqoM61cHLi9HJkydnhv9mLclEwSGZd5jIS5J5x0sXKeccJgqSKUsyVXC83Fv2kstixarES+cx/QdbOrsqoENAk5QVUtSQpcJKU2FPUWJnKXGylNgZonaWiJUlYmcIW1nClpuHCmmChQyBQpZAIUOgkEG3LUinEOnUIm2+PFKANECaIH2eQJnCXTYFjl/D9mkQ0HF8mrdNIA1AF0hDInWJ1CToDlJzABtdpgmS5qLj/DQvXebfWUqwJOSlRt4R5KW4SK5dZL1OQQoKUvNyQTB1ClZTZIA7gFNSym4AIcR3gceBuSLzOPD30h1msVsIUSqEqJNSLrDT9MpIK01HJMAPBmLk5BiVsTo+8/DvEw5c5KZWKG4wAoYgYFxZmC6HIyVpT3BSBUnKcmbFZ85y2nKfsLOWJGNLcrabZ63psjsnKTNdz56tn7UleQeyNmRtQYIop4m6/0RL7byXEr9TIOTkCTp5QnaekJMjaOcJOXlCdo6AUyDk5AjZ03VyM/WD59X3OxY+aeF3CvgcN/c7Fr60V57eJm0WL8QgNZC6J1yGcMuGVzaYFSjdrYcOUjtv2duObs8c74J9DIHUpusDlxhMMj44TKZqiZ/9AljI19rA/MGafVzYSrlYnQbgoiJz8uTJRZjoYqQyhPObyckBmsrXc1fbRxjuH2OYy7xBvY5ZymdwPVHs9oO6hisR8tK8/x/TS1eBLSHvQEFCwQFLCgoSrJl1AktyhXXesrfdLfuxpB9buudwvHPZUpCTkJEwPLPOq+Mdxzlv/fQ2W7oyMrvOPZctBVJKDMfC5xQwHQvTLmBK280dC59dwJQWpm2hSRtDOhhebkprZtmUtps7c7ZZDmbBnrddlw6mY3l15m8zve0zCWfesoZElza6cDCEja456MJNRzevZR1Xdy+dP6F2LgsRmYvJ3/ldkwupsyCDLk07csdOHjvwJjt33FfU72DmdnMUI8VuP6hruF5wr6Fttc24Ktxr2LygulJKHE/UHDmd5MyylO4fpzNdb15dd53Eredum1Nvept3HuntZzN7XDlTljPrtkUMpvq7V+xeWojI9AFNc5YbgfN9OiykzlUjhKAiUlfUAqNQKG5chBDogvPmE63+/9lKtukX4sRnL9AuhFgrhPABTwLPnVfnOeCXhcudQHI538coFAqFoji5YktGSmkJIf4t8AKuAP+dlPKoEOJz3vavAj/CHb58CncI86+unMkKhUKhKBYWNJ5DSvkjXCGZu+6rc8oS+PzymqZQKBSKYkf5vFYoFArFiqFERqFQKBQrhhIZhUKhUKwY4lrFQkgmkyrogkKhUPyME4vF5o3JVi0ZhUKhUKwYSmQUCoVCsWJcs+4yhUKhUNx4qJaMQqFQKFaMohEZIcSjQojjQohTQog/XG17FosQokkI8aoQolMIcVQI8burbdNSEULoQoj9QogfrLYtS8ELRfE9IcQx7/u4a7VtWixCiH/v3UdHhBDfEUJc90GVhBB/J4QYFkIcmbOuXAjxkhDipJeXraaNV+IS1/A/vHvpkBDin4QQpZc7xmpyMfvnbPs9IYQUQlQu5zmLQmTmBE77ELAZ+EUhxMLcnl4/WMB/lFJuAu4EPl+E1zDN7wKdq23EVfCXwL9IKTcCN1Fk1yKEaAB+B7hNSrkV193Tk6tr1YL4BvDoeev+EPiJlLId+Im3fD3zDS68hpeArVLK7cAJ4IvX2qhF8A0utB8hRBNuYMqe5T5hUYgMcwKnSSnzwHTgtKJBSjkwHZJaSjmJ+8fWsLpWLR4hRCPwGPC11bZlKQghSoD7ga8DSCnzUsrE6lq1JAwgKIQwcMO+LLvX8+VGSvkGXBAA6nHgm175m8DHrqlRi+Ri1yClfFFKaXmLu3G90F+XXOI7APgy8J+4TIiWpVIsInOpoGhFiRCiBbgZ2LO6liyJ/4V7MzqrbcgSWQeMAP/X6/L7mhAivNpGLQYp5Tngf+I0hMULAAADJ0lEQVQ+dQ7gej1/cXWtWjI10x7bvbx6le25Wn4N+PFqG7EYhBAfBc5JKQ+uxPGLRWQWFRTtekYIEQGeAb4gpZxYbXsWgxDiI8CwlPK91bblKjCAW4CvSClvBqa4/rto5uG9t3gcWAvUA2EhxKdW1yqFEOKPcLvFn1ptWxaKECIE/BHwX1bqHMUiMtckKNpKI4QwcQXmKSnls6ttzxK4B/ioEOIMbpflg0KIb6+uSYumD+iTUk63Ir+HKzrFxEPAaSnliJSyADwL3L3KNi2VISFEHYCXD6+yPUtCCPErwEeAX5LFNS+kFfdh5aD3u24E9gkhapfrBMUiMgsJnHZdI9xwnl8HOqWUf7Ha9iwFKeUXpZSNUsoW3O/gFSllUT1BSykHgV4hxAZv1QeAjlU0aSn0AHcKIULeffUBimzwwhyeA37FK/8K8P1VtGVJCCEeBf4A+KiUMr3a9iwGKeVhKWW1lLLF+133Abd4v5NloShExnupNh04rRN4Wkp5dHWtWjT3AJ/Gffo/4KUPr7ZRNyj/DnhKCHEI2AH86Srbsyi8Vtj3gH3AYdzf8d+uqlELQAjxHeBtYIMQok8I8evAl4APCiFO4o5u+tJq2nglLnENfwVEgZe83/VXL3uQVeQS9q/sOYurZadQKBSKYqIoWjIKhUKhKE6UyCgUCoVixVAio1AoFIoVQ4mMQqFQKFYMJTIKhUKhWDGUyCgUCoVixVAio1AoFIoVQ4mMQqFQKFYMJTIKxTIhhIidH/RJCLFZCDHkhRhQKG44lMgoFMuElDKJ6/tpbjC6PwH+e7F53FYolgslMgrF8nIE2AQghLgD18PzX6+qRQrFKqJERqFYXo4w25L5b8B/lVLmVtEehWJVUSKjUCwvR4DNQoiHgDrg71fZHoViVVEio1AsL9MtmT8F/rOU0l5lexSKVUW5+lcolhEhRBBIAe9IKe9abXsUitVGiYxCoVAoVgzVXaZQKBSKFUOJjEKhUChWDCUyCoVCoVgxlMgoFAqFYsVQIqNQKBSKFUOJjEKhUChWDCUyCoVCoVgxlMgoFAqFYsVQIqNQKBSKFeP/Ayn8rQeSFwKtAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# NO CODE\n",
"x = np.arange(0, 14, 0.01)\n",
"y2 = stats.chi2.pdf(x, 2)\n",
"y3 = stats.chi2.pdf(x, 3)\n",
"y4 = stats.chi2.pdf(x, 4)\n",
"y5 = stats.chi2.pdf(x, 5)\n",
"plt.plot(x, y2, lw=2, label='2 df')\n",
"plt.plot(x, y3, lw=2, label='3 df')\n",
"plt.plot(x, y4, lw=2, label='4 df')\n",
"plt.plot(x, y5, lw=2, label='5 df')\n",
"plt.legend()\n",
"plt.xlabel('$v$')\n",
"plt.title('Chi-Squared $(n)$ Densities for $n = 2, 3, 4, 5$');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The chi-squared (2) distribution is exponential because it is the gamma $(1, 1/2)$ distribution. This distribution has three names:\n",
"\n",
"- chi-squared (2)\n",
"- gamma (1, 1/2)\n",
"- exponential (1/2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Mean and Variance ###\n",
"You know that if $T$ has the gamma $(r, \\lambda)$ density then \n",
"\n",
"$$\n",
"E(T) ~ = ~ \\frac{r}{\\lambda} ~~~~~~~~~~~~ SD(T) = \\frac{\\sqrt{r}}{\\lambda}\n",
"$$\n",
"\n",
"If $X$ has the chi-squared $(n)$ distribution then $X$ is gamma $(n/2, 1/2)$. So\n",
"\n",
"$$\n",
"E(X) ~ = ~ \\frac{n/2}{1/2} ~ = ~ n\n",
"$$\n",
"\n",
"Thus **the expectation of a chi-squared random variable is its degrees of freedom**.\n",
"\n",
"The SD is\n",
"$$\n",
"SD(X) ~ = ~ \\frac{\\sqrt{n/2}}{1/2} ~ = ~ \\sqrt{2n}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"### Estimating the Normal Variance ###\n",
"Suppose $X_1, X_2, \\ldots, X_n$ are i.i.d. normal $(\\mu, \\sigma^2)$ variables, and that you are in a setting in which you know $\\mu$ and are trying to estimate $\\sigma^2$. \n",
"\n",
"Let $Z_i$ be $X_i$ in standard units, so that $Z_i = (X_i - \\mu)/\\sigma$. Define the random variable $T$ as follows:\n",
"\n",
"$$\n",
"T ~ = ~ \\sum_{i=1}^n Z_i^2 ~ = ~ \\frac{1}{\\sigma^2}\\sum_{i=1}^n (X_i - \\mu)^2\n",
"$$\n",
"\n",
"Then $T$ has the chi-squared $(n)$ distribution and $E(T) = n$. Now define $W$ by\n",
"\n",
"$$\n",
"W ~ = ~ \\frac{\\sigma^2}{n} T ~ = ~ \\frac{1}{n} \\sum_{i=1}^n (X_i - \\mu)^2\n",
"$$\n",
"\n",
"Then $W$ can be computed based on the sample since $\\mu$ is known. And since $W$ is a linear tranformation of $T$ it is easy to see that $E(W) = \\sigma^2$. \n",
"\n",
"So we have constructed an unbiased estimate of $\\sigma^2$. It is the mean squared deviation from the known population mean.\n",
"\n",
"But typically, $\\mu$ is not known. In that case you need a different estimate of $\\sigma^2$ since you can't compute $W$ as defined above. You showed in exercises that\n",
"\n",
"$$\n",
"S^2 ~ = ~ \\frac{1}{n-1}\\sum_{i=1}^n (X_i - \\bar{X})^2\n",
"$$\n",
"\n",
"is an unbiased estimate of $\\sigma^2$ regardless of the distribution of the $X_i$'s. When the $X_i$'s are normal, as is the case here, it turns out that $S^2$ is a linear transformation of a chi-squared $(n-1)$ random variable. We will show that later in the course."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### \"Degrees of Freedom\" ###\n",
"The example above helps explain the strange term \"degrees of freedom\" for the parameter of the chi-squared distribution. \n",
"- When $\\mu$ is known, you have $n$ independent centered normals $(X_i - \\mu)$ that you can use to estimate $\\sigma^2$. That is, you have $n$ degrees of freedom in constructing your estimate.\n",
"- When $\\mu$ is not known, you are using all $n$ of $X_1 - \\bar{X}, X_2 - \\bar{X}, \\ldots, X_n - \\bar{X}$ in your estimate, but they are not independent. They are the deviations of the list $X_1, X_2, \\ldots , X_n$ from their average $\\bar{X}$, and hence their sum is 0. If you know $n-1$ of them, the final one is determined. So you only have $n-1$ degrees of freedom."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"celltoolbar": "Tags",
"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.8.4"
}
},
"nbformat": 4,
"nbformat_minor": 1
}