{
"cells": [
{
"cell_type": "code",
"execution_count": 2,
"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": "code",
"execution_count": 3,
"metadata": {
"tags": [
"remove_cell"
]
},
"outputs": [],
"source": [
"# NO CODE\n",
"\n",
"def simulate_T1_T2(N, n):\n",
" tanks = np.arange(1, N+1)\n",
" sample = np.random.choice(tanks, size=n, replace=False)\n",
" t1 = 2*np.mean(sample) - 1\n",
" t2 = max(sample)*(n+1)/n - 1\n",
" return [t1, t2]\n",
"\n",
"def compare_T1_T2(N, n, repetitions):\n",
" tbl = Table(['T_1 = 2*Mean-1', 'T_2 = Augmented Max'])\n",
" for i in np.arange(repetitions):\n",
" tbl.append(simulate_T1_T2(N, n))\n",
" \n",
" tbl.hist(bins=np.arange(N/2, 3*N/2))\n",
" plt.title('$N =$'+str(N)+', $n =$'+str(n)+' ('+str(repetitions)+' repetitions)');\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Additivity ##"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calculating expectation by plugging into the definition works in simple cases, but often it can be cumbersome or lack insight. The most powerful result for calculating expectation turns out not to be the definition. It looks rather innocuous:\n",
"\n",
"### Additivity of Expectation ###\n",
"Let $X$ and $Y$ be two random variables defined on the same probability space. Then\n",
"\n",
"$$\n",
"E(X+Y) = E(X) + E(Y)\n",
"$$\n",
"\n",
"Before we look more closely at this result, note that we are assuming that all the expectations exist; we will do this throughout in this course. \n",
"\n",
"And now note that **there are no assumptions about the relation between $X$ and $Y$**. They could be dependent or independent. Regardless, the expectation of the sum is the sum of the expectations. This makes the result powerful."
]
},
{
"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: Additivity of Expectation\n",
"from IPython.display import YouTubeVideo\n",
"\n",
"YouTubeVideo('HzZEhM4NHUQ')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Additivity follows easily from the definition of $X+Y$ and the definition of expectation on the domain space. First note that the random variable $X+Y$ is the function defined by\n",
"\n",
"$$\n",
"(X+Y)(\\omega) = X(\\omega) + Y(\\omega) ~~~~ \\text{for all }\n",
"\\omega \\in \\Omega\n",
"$$\n",
"\n",
"Thus a \"value of $X+Y$ weighted by the probability\" can be written as\n",
"\n",
"$$\n",
"(X+Y)(\\omega) \\cdot P(\\omega) = X(\\omega)P(\\omega) + \n",
"Y(\\omega)P(\\omega )\n",
"$$\n",
"\n",
"Sum the two sides over all $\\omega \\in \\Omega$ to prove additivty of expecation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{admonition} Quick Check\n",
"Let $X$ and $Y$ be random variables on the same space, with $E(X) = 5$ and $E(Y) = 3$.\n",
"\n",
"(a) Find $E(X-Y)$.\n",
"\n",
"(b) Find $E(2X-8Y+7)$.\n",
"\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{admonition} Answer\n",
":class: dropdown\n",
"(a) $2$ because $X-Y = X+(-Y)$\n",
"\n",
"(b) $-7$\n",
"\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"By induction, additivity extends to any finite number of random variables. If $X_1, X_2, \\ldots , X_n$ are random variables defined on the same probability space, then\n",
"\n",
"$$\n",
"E(X_1 + X_2 + \\cdots + X_n) = E(X_1) + E(X_2) + \\cdots + E(X_n)\n",
"$$\n",
"\n",
"regardless of the dependence structure of $X_1, X_2, \\ldots, X_n$.\n",
"\n",
"If you are trying to find an expectation, then the way to use additivity is to write your random variable as a sum of simpler variables whose expectations you know or can calculate easily. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### $E(X^2)$ for a Poisson Variable $X$ ###\n",
"\n",
"Let $X$ have the Poisson $\\mu$ distribution. In earlier sections we showed that $E(X) = \\mu$ and $E(X(X-1)) = \\mu^2$.\n",
"\n",
"Now $X^2 = X(X-1) + X$. The random variables $X(X-1)$ and $X$ are both functions of $X$, so they are not independent of each other. But additivity of expectation doesn't require independence, so we can use it to see that\n",
"\n",
"$$\n",
"E(X^2) ~ = ~ E(X(X-1)) + E(X) ~ = ~ \\mu^2 + \\mu\n",
"$$\n",
"\n",
"We will use this fact later when we study the variability of $X$. \n",
"\n",
"It is worth noting that it is not easy to calculate $E(X^2)$ directly, since\n",
"\n",
"$$\n",
"E(X^2) ~ = ~ \\sum_{k=0}^\\infty k^2 e^{-\\mu}\\frac{\\mu^k}{k!}\n",
"$$\n",
"\n",
"is not an easy sum to simplify."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Sample Sum ###\n",
"Let $X_1, X_2, \\ldots , X_n$ be a sample drawn at random from a numerical population that has mean $\\mu$, and let the sample sum be \n",
"\n",
"$$\n",
"S_n = X_1 + X_2 + \\cdots + X_n\n",
"$$\n",
"\n",
"Then, regardless of whether the sample was drawn with or without replacement, each $X_i$ has the same distribution as the population. This is clearly true if the sampling is with replacement, and it is true by symmetry if the sampling is without replacement as we saw in an earlier chapter.\n",
"\n",
"So, regardless of whether the sample is drawn with or without replacement, $E(X_i) = \\mu$ for each $i$, and hence\n",
"\n",
"$$\n",
"E(S_n) = E(X_1) + E(X_2) + \\cdots + E(X_n) = n\\mu\n",
"$$\n",
"\n",
"We can use this to estimate a population mean based on a sample mean."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Unbiased Estimator ###\n",
"\n",
"Suppose a random variable $X$ is being used to estimate a fixed numerical parameter $\\theta$. Then $X$ is called an *estimator* of $\\theta$. \n",
"\n",
"The *bias* of $X$ is the difference $E(X) - \\theta$. The bias measures the amount by which the estimator exceeds the parameter, on average. The bias can be negative if the estimator tends to underestimate the parameter.\n",
"\n",
"If the bias of an estimator is $0$ then the estimator is called *unbiased*. So $X$ is an unbiased estimator of $\\theta$ if $E(X) = \\theta$.\n",
"\n",
"If an estimator is unbiased, and you use it to generate estimates repeatedly and independently, then in the long run the average of all the estimates is equal to the parameter being estimated. On average, the unbiased estimator is neither higher nor lower than the parameter. That's usually considered a good quality in an estimator.\n",
"\n",
"In practical terms, if a data scientist wants to estimate an unknown parameter based on a random sample $X_1, X_2, \\ldots, X_n$, the data scientist has to come up with a *statistic* to use as the estimator. \n",
"\n",
"Recall from Data 8 that a statistic is a number computed from the sample. In other words, a statistic is a numerical function of $X_1, X_2, \\ldots, X_n$.\n",
"\n",
"Constructing an unbiased estimator of a parameter $\\theta$ therefore amounts to finding a statistic $T = g(X_1, X_2, \\ldots, X_n)$ for a function $g$ such that $E(T) = \\theta$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Unbiased Estimators of a Population Mean ###\n",
"\n",
"As in the sample sum example above, let $S_n$ be the sum of a sample $X_1, X_2, \\ldots , X_n$ drawn at random from a population that has mean $\\mu$. The standard statistical notation for the average of $X_1, X_2, \\ldots , X_n$ is $\\bar{X}_n$. So\n",
"\n",
"$$\n",
"\\bar{X}_n = \\frac{S_n}{n}\n",
"$$\n",
"\n",
"Then, regardless of whether the draws were made with replacement or without,\n",
"\n",
"$$\n",
"\\begin{align*}\n",
"E(\\bar{X}_n) &= \\frac{E(S_n)}{n} ~~~~ \\text{(linear function rule)} \\\\\n",
"&= \\frac{n \\mu}{n} ~~~~~~~~~ \\text{(} E(S_n) = n\\mu \\text{)} \\\\\n",
"&= \\mu\n",
"\\end{align*}\n",
"$$\n",
"\n",
"Thus the sample mean is an unbiased estimator of the population mean.\n",
"\n",
"It is worth noting that $X_1$ is also an unbiased estimator of $\\mu$, since $E(X_1) = \\mu$. So is $X_j$ for any $j$, also $(X_1 + X_9)/2$, or any linear combination of the sample if the coefficients add up to 1.\n",
"\n",
"But it seems clear that using the sample mean as the estimator is better than using just one sampled element, even though both are unbiased. This is true, and is related to how variable the estimators are. We will address this later in the course."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{admonition} Quick Check\n",
"Let $X_1, X_2, X_3$ be i.i.d. Poisson $(\\mu)$ random variables, and suppose the value of $\\mu$ is unknown. Is $0.4X_1 + 0.2X_2 + 0.4X_3$ an unbiased estimator of $\\mu$?\n",
"\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{admonition} Answer\n",
":class: dropdown\n",
"Yes\n",
"\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"tags": [
"remove-input",
"hide-output"
]
},
"outputs": [
{
"data": {
"text/html": [
"\n",
" VIDEO"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# VIDEO: Example of an Unbiased Estimator\n",
"\n",
"YouTubeVideo('ruEpGZJwHmw')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### First Unbiased Estimator of a Maximum Possible Value ###\n",
"\n",
"Suppose we have a sample $X_1, X_2, \\ldots , X_n$ drawn at random from $1, 2, \\ldots , N$ for some fixed $N$, and we are trying to estimate $N$. \n",
"\n",
"How can we use the sample to construct an unbiased estimator of $N$? By definition, such an estimator must be a function of the sample and its expectation must be $N$.\n",
"\n",
"In other words, we have to construct a statistic that has expectation $N$.\n",
"\n",
"Each $X_i$ has the uniform distribution on $1, 2, \\ldots , N$. This is true for sampling with replacement as well as for simple random sampling, by symmetry. \n",
"\n",
"The expectation of each of the uniform variables is $(N+1)/2$, as we have seen earlier. So if $\\bar{X}_n$ is the sample mean, then\n",
"\n",
"$$\n",
"E(\\bar{X}_n) = \\frac{N+1}{2}\n",
"$$\n",
"\n",
"Clearly, $\\bar{X}_n$ is not an unbiased estimator of $N$. That's not surprising because $N$ is the maximum possible value of each observation and $\\bar{X}_n$ should be somewhere in the middle of all the possible values.\n",
"\n",
"But because $E(\\bar{X}_n)$ is a linear function of $N$, we can figure out how to create an unbiased estimator of $N$. \n",
"\n",
"Remember that our job is to create a function of the sample $X_1, X_2, \\ldots, X_n$ in such a way that the expectation of that function is $N$.\n",
"\n",
"Start by inverting the linear function, that is, by isolating $N$ in the equation above.\n",
"\n",
"$$\n",
"2E(\\bar{X}_n) - 1 = N\n",
"$$\n",
"\n",
"This tells us what we have to do to the sample $X_1, X_2, \\ldots, X_n$ to get an unbiased estimator of $N$.\n",
"\n",
"We should just use the statistic $T_1 = 2\\bar{X}_n - 1$ as the estimator. It is unbiased because $E(T_1) = N$ by the calculation above."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{admonition} Quick Check\n",
"In the setting above, what is the bias of $2\\bar{X}_n$ as an estimator of $N$? Does it tend to overestimate on average, or underestimate?\n",
"\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{admonition} Answer\n",
":class: dropdown\n",
"$1$; overestimate\n",
"\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Second Unbiased Estimator of the Maximum Possible Value ###\n",
"\n",
"The calculation above stems from a problem the Allied forces faced in World War II. Germany had a seemingly never-ending fleet of Panzer tanks, and the Allies needed to estimate how many they had. They decided to base their estimates on the serial numbers of the tanks that they saw.\n",
"\n",
"Here is a picture of one from [Wikipedia](https://en.wikipedia.org/wiki/Panzer_IV).\n",
"\n",
"![Panzer Tank](panzer.png)\n",
"\n",
"Notice the serial number on the top left. When tanks were disabled or destroyed, it was discovered that their parts had serial numbers too. The ones from the gear boxes proved very useful.\n",
"\n",
"The idea was to model the observed serial numbers as random draws from $1, 2, \\ldots, N$ and then estimate $N$. This is of course a very simplified model of reality. But estimates based on even such simple probabilistic models proved to be quite a bit [more accurate](https://en.wikipedia.org/wiki/German_tank_problem#Specific_data) than those based on the intelligence gathered by the Allies. For example, in August 1942, intelligence estimates were that Germany was producing 1,550 tanks per month. The prediction based on the probability model was 327 per month. After the war, German records showed that the actual production rate was 342 per month.\n",
"\n",
"The model was that the draws were made at random without replacement from the integers 1 through $N$. \n",
"\n",
"In the example above, we constructed the random variable $T$ to be an unbiased estimator of $N$ under this model.\n",
"\n",
"The Allied statisticians instead started with $M$, the sample maximum:\n",
"\n",
"$$\n",
"M ~ = ~ \\max\\{X_1, X_2, \\ldots, X_n\\}\n",
"$$\n",
"\n",
"The sample maximum $M$ is a biased estimator of $N$, because we know that its value is always less than or equal to $N$. Its average value therefore will be somewhat less than $N$.\n",
"\n",
"To correct for this, the Allied statisticians imagined a row of $N$ spots for the serial numbers $1$ through $N$, with marks at the spots corresponding to the observed serial numbers. The visualization below shows an outcome in the case $N= 20$ and $n = 3$.\n",
"\n",
"![gaps](all_gaps.png)\n",
"\n",
"- There are $N = 20$ spots in all. \n",
"- From these, we take a simple random sample of size $n = 3$. Those are the gold spots.\n",
"- The remaining $N - n = 17$ spots are colored blue.\n",
"\n",
"The $n = 3$ sampled spots create $n+1 = 4$ blue \"gaps\" between sampled values: one before the leftmost gold spot, two between successive gold spots, and one after the rightmost gold spot that is at position $M$.\n",
"\n",
"A key observation is that because of the symmetry of simple random sampling, **the lengths of all four gaps have the same distribution.** \n",
"\n",
"But of course we don't get to see all the gaps. In the sample, we can see all but the last gap, as in the figure below. The red question mark reminds you that the gap to the right of $M$ is invisible to us.\n",
"\n",
"![mystery gap](mystery_gap.png)\n",
"\n",
"If we could see the gap to the right of $M$, we would see $N$. But we can't. So we can try to do the next best thing, which is to augment $M$ by the estimated size of that gap.\n",
"\n",
"Since we can see all of the spots and their colors up to and including $M$, we can see $n$ out of the $n+1$ gaps. The lengths of the gaps all have the same distribution by symmetry, so we can estimate the length of a single gap by the average length of all the gaps that we can see.\n",
"\n",
"We can see $M$ spots, of which $n$ are the sampled values. So the total length of all $n$ visible gaps is $M-n$. Therefore\n",
"\n",
"$$\n",
"\\text{estimated length of one gap} ~ = ~ \\frac{M-n}{n}\n",
"$$\n",
"\n",
"So the Allied statisticians decided to improve upon $M$ by using the *augmented maximum* as their estimator:\n",
"\n",
"$$\n",
"T_2 ~ = ~ M + \\frac{M-n}{n}\n",
"$$\n",
"\n",
"By algebra, this estimator can be rewritten as\n",
"\n",
"$$\n",
"T_2 ~ = ~ M\\cdot\\frac{n+1}{n} ~ - ~ 1\n",
"$$\n",
"\n",
"Is $T_2$ an unbiased estimator of $N$? To answer this, we have to find its expectation. Since $T_2$ is a linear function of $M$, we'll find the expectation of $M$ first.\n",
"\n",
"Here once again is the visualization of what's going on.\n",
"\n",
"![gaps](all_gaps.png)\n",
"\n",
"Let $G$ be the length of the last gap. Then $M = N - G$.\n",
"\n",
"There are $n+1$ gaps, made up of the $N-n$ unsampled values. Since they all have the same expected length,\n",
"\n",
"$$\n",
"E(G) ~ = ~ \\frac{N-n}{n+1}\n",
"$$\n",
"\n",
"So\n",
"\n",
"$$\n",
"E(M) ~ = ~ N - \\frac{N-n}{n+1} ~ = ~ (N+1)\\frac{n}{n+1}\n",
"$$\n",
"\n",
"Recall that the Allied statisticians' estimate of $N$ is\n",
"\n",
"$$\n",
"T_2 ~ = ~ M\\cdot\\frac{n+1}{n} - 1\n",
"$$\n",
"\n",
"Now\n",
"\n",
"$$\n",
"E(T_2) ~ = ~ E(M)\\cdot\\frac{n+1}{n} - 1 ~ = ~ (N+1)\\frac{n}{n+1}\\cdot\\frac{n+1}{n} - 1 ~ = ~ N\n",
"$$\n",
"\n",
"Thus the augmented maximum $T_2$ is an unbiased estimator of $N$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{admonition} Quick Check\n",
"A gardener in Berkeley has 23 blue flower pots in a row. She picks a simple random sample of 5 of them and colors the selected pots gold. What is the expected number of blue flower pots at the end of the row?\n",
"\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{admonition} Answer\n",
":class: dropdown\n",
"$3$\n",
"\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Which Estimator to Use? ###\n",
"The Allied statisticians thus had two unbiased estimators of $N$ from which to choose. They went with $T_2$ instead of $T_1$ because $T_2$ has less variability.\n",
"\n",
"We will quantify this later in the course. For now, here is a simulation of distributions of the two estimators in the case $N = 300$ and $n=30$. The simulation is based on $5000$ repetitions of drawing a simple random sample of size $30$ from the integers $1$ through $300$."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"tags": [
"remove_input"
]
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlUAAAEICAYAAAB2/gEGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU1f3/8ddnJishQlhlFUQUEK2C4IK71qpV8Fer/eoXv34t9SdVf7QutYK1de23tdVuX1uqQt3qrq0tirXWpW4gCi7IFkSBhADZyJ6ZZOb8/riTOIQkTMgkk5m8n4/HPDJz7rn3fk5uQj6cc+655pxDRERERDrHl+gARERERFKBkioRERGROFBSJSIiIhIHSqpERERE4kBJlYiIiEgcKKkSERERiQMlVSIiIiJxoKRKREREJA6UVIkkiJn1NbNCM5uW6FhSgZk9aGavxKtePJnZH83sl915ThHpfkqqkpCZPW1mzsxua1E+PFL+1Tie6xIz+8DMys2szszWmtl1ZmYt6p1tZh+aWcDMvjCza1s51l7r9CTxbHsbfgi875xbEXWsWyLXsOXroI6eMxWvSRMze8XMHmxR/D3ggr3U2aNeN7kN+K6ZHdjN5xWRbpSW6ABknxwFbAO+Bfw4qrypx+P9OJ5rJ3A7sB4IACcAvwcagd8AmNlRwPPA3cBFwNHAQjOrdc4tjLVODxSXtrfGzLKA7wL/1crmL4BjW5QVR+0bl+93vK6JmWU454Kx1u8qzrmKeNaLJ+dcoZn9C7gSuL67zy8i3cQ5p1cSvYCBgAP+LxACvhK17Q5gYzfE8BfgL1GfHwPeaVHnF8DnHanTgfPnRb4H/wk8DVQCO4AremLb2zjOeUANkNai/Ja9XcN4fb87EfvrwCK8hLMIKI7a9v+AdUA9kA/cFN3GyL6LgZ8BJZFr9wCQ3eIcbR4HeDBy/aNfJ0fKX2mvTtS2V6LOlR6JpxAIAmuAi1tp8wPAzcB2oCxynJyoOscDbwNVkddHwNeitl8GbO/qn1G99NIrcS8N/yWfpt6ov+H9A35hi2179FKZ2QIzq97La8HeTmye6cAM4LWoTTOAl1pUfwkYY2YjO1AnVkdGvl4LPAocATwM/K+Z5bSIuSe0vTUnAaucc42tbBtpZgWR11IzO67F9nh9vztzTS4EBgOnAaeCN3SJ1wszH5iIN8x2BfCTFvt+E+8/ByfgJcYzgZ83bYzhON8D3gSeAoZFXu+0OEcsdZr8FLgc+D4wGe9n6lEzO62VuAfgJXAX4yXGN0Ri9uP9Ti4HpkRetwC1UfsvB4aa2cQ24hCRJKfhv+RzFFDonNtuZk/h/fG4KbJtKvA/reyzEO+PS3vK2tpgZv3w/hefAfiBW51zv42qMgzvf+/RtkdtK4ixTqyOxBuC+w/nXH4kxofx/hDn4fUANekJbW/N2MhxW1qONyS4DuiHN0T4ppmd6Zz7ZwfO2dXXpAi40jkXBjCzPngJxjecc02J2udm9iPgt3g9PE3KgLnOuRCwNlLnd2Y2H69Hqd3jOOcqzCwI1DnnmuOPnurWVp2WInHPA65xzj0dKf6peTcP3AT8K6r6FufcNZH368zsCeAMvGRvP7yfvb81/Uzi9bBFa/p+HgisbSsmEUleSqqSzzSgaWLzM8BvzGwKsAvvf/979FQ558poJ3GIQRVeb1Af4Djgf8xsm3PugRj2dXGqE+1I4LWoP14A4/F6BbbtduCe2/ZsYI+5Pc65pS2K3jSzEcAPgH+2rN/Bc8azzgdNCVXEoXhtetbMovf1A1lmNtg51zQv7L1IQtXkbbykdRyQ2YHjxMNBkXP/u0X5G3g9ZdE+bPG5EC+pwjlXbmYPAP8ws1cj+//FObc+qn595Gt2PAIXkZ5Hw3/J5ygiiVPkf+Bv4Q3FTAPCwMqWO3R2CMw5F3bObXTOfey8Ccx34c3falIE7N9it6GRr9s7UCdWR7LnUM4U4OMWf+h7SttbU4w3lBSLd4ExHTxnV1+Tmhafm/4tuQAvCW16HYaX8LaX2EbfTdmZ43RGyyTSWilrORnfEfVvqHPucrze4n/iDe+uNrMrouo3Xe94JoUi0oOopyqJmNkwYDhf9lSBN7R1Hd61XO+cq2pl104NgbXCh9ej0ORt4Gt4t403ORPY7Jwr6ECdvTKzbOAQ9uyRm0IrCSU9o+2tWQlcHeM5jwS2dvCc3XZNIj7F64k50Dn34l7qTjMzf1Rv1bF4CctneMlMLMcJ4vVetSeWOhvx7uw8Ca8NTU5s8TkmzrnVwGrgHjNbiHdDyR8jmw/Du7lkVUePKyLJQUlVcmltyYRn8OaazAZebm2nzgyBmdmteBN+N+HdJXUi3vpKf4qq9ivgHTO7E3gEmI5399Y1HawTi8Px/lB+0KJ8CvBsy8o9pO2tWQrcbWajnHPNCZOZ3QMswVtWYT+8CdRfBWZ18JzdeU1wzlWb2U/x5iOB11uThpdIHOmc+2FU9YHAvWb2G7z5RbcD9zvnaiLfg1iO8zlwipmNwxtGbW2ZhD3qOOcaWsRda2a/BW43s2K8Ib4L8L7fMa/3Zt46YpcDf8dLgIfjTcSPTvRPBt5yzlXGelwRSTKJvv1Qr9hfeD0Km1opfw1vKGJeF5zzV3j/m68DyvGSmasAf4t6X8e7hTwAbAaubeVY7dYB/jvSjjHtxDMXb6J+dNmIyH5H9tS2t3H814AFLcoex5vQHMBbJ+sV4NSOfi/jWaeVfV4HHmhj2xy8xKQ+8j1bDny3xb6L8ZZuKMWbs7YY6NPB4xyINw+qmlaWVGirTqS8Zb2Yl1RoUfYj4IvI+2HAc1HXbhtwP9Avst3wkryL4v07qpdeevWclznX0TnCIl3DvBXiz8dbe6u1pQZSipmdADwBjHfO1e6tfiows9fx1uH6TqJj6U5mdiHeHZBHuN0n6YtICtFEdelJzgGu7g0JFYBz7k3gVrzlFSS1ZQKXKaESSW3qqRKRbtNbe6pEpHdQUiUiIiISBxr+ExEREYmDbllSwcwW482X2emcmxwpGwA8ibeo4RfAhc658uj9Kioq1I0mIpLi+vXrZ3uvJdLzdVdP1YN4iwpGuxH4l3NuPN7ztW7splhERERE4q5bkirn3L/ZcwHGWcBDkfcP4T3xXURERCQpJXJO1VDnXBFA5OuQrjxZfn7LB8anHrUxNaiNqUFtFOl9uu3uPzMbAyyJmlO1yznXP2p7uXMuL3qf6DlV+uUVEUkd48ePb36vOVWSKhL57L8dZjbMOVcUeVDwzvYqR/8C7ov8/PxOH6OnUxtTg9qYGtRGkd4nkUnV34BL8Z65dSnwfAJjERGRFPPBBx/4srKybkhPT5+IlhCSzgs3NDSsra+vv2vq1Knh1ip015IKj+M98HSQmRUAP8FLpp4ysznAFrwnw4uIiMRFVlbWDYMHD74wMzOz1T+AIh0VCAQOKy4uBi+H2UO3JFXOuYva2HRad5xfRER6n/T09IlKqCSeMjMzw5Gez1apO1RERFKV/sZJV2jz50o/cCIiIiJxkMiJ6iLSCda4GSwL5x+a6FBEkkJRUc2o4uK6jHgdb/Dg7OCwYTlb29q+c+dO/9e//vVDAEpLS9N9Pp/Ly8trBHjzzTfXZmVl7bam0aWXXjrmjTfe6DdgwIDGlStXftrZ+DZt2pR++eWXjy0pKUn3+XxcdNFFxTfccEPznfYLFy4ceNppp1WNGzcu6PP5uOSSS8a8/PLLeevXr/+of//+YYArr7xy1GOPPTZk3bp1H+2///6NnY2pLXfffffgxYsXDy0sLMzs6nN1JSVVIknK37iBUJpuZxeJVXFxXcZNNy2L25pYd955TMawYTltbh8yZEhoxYoVawAWLFgwPCcnJ3TTTTftaKv+7NmzS+bOnbvzqquuGhuP+NLT07njjjsKjj322Npdu3b5TjjhhElnnHFGZV5eXugnP/nJ8JEjRwbfeOONvj/72c9yFy1atBlg5MiRgWeeeab/d77znbJQKMTy5ctzBw8e3BCPeNpz/PHHV8+cObPi3HPPPaSrz9WVNPwnIiLSA3z1q1+tHjRoUNx6aEaNGtVw7LHH1gL0798/PHbs2LqCgoKMAw44oOHOO+8sfOaZZwb99a9/HXDfffdtbtrn3HPPLXv++ecHALz88su5U6ZMqfb7/c09ag888MCA4447buK0adMmzZkz54DGRi/cyy+/fPQxxxwz8cgjjzx0/vz5w5vqT5gw4bD58+cPnz59+sQpU6ZM+vjjj7Nai/Xoo4+uGz9+fDBebU8UJVUiIiJJYtGiRQOmTZs2qeXr/PPPP7C9/fLz8zPWrVvXZ8aMGdVbtmxJv/nmm0d885vfLJk1a1bZ3LlzRzfVO+iggwJlZWVpJSUl/qeffnrABRdc0Pzc3o8++ijr+eefH/D666+vW7FixRq/3+8WL148EOCnP/1p4bJly9auWLHi0+XLl+e+//772U37DRw4sPG9995be8kllxT/6le/Sun5Chr+E0kmLoA/8BahjGMTHYmIJMCcOXPK5syZU7b3ml+qrKz0zZ49e9wtt9yyNS8vL5yXlxdevHjx5oULFw48+eSTq1se78wzzyx/9NFHB3z00Uc5999/f3Mv1iuvvJK7Zs2aPjNmzJgIEAgEfE09a4899tiAP//5z4NCoZCVlJSkr169Ouuoo46qA7jgggvKAaZOnVr74osv7vY4ulSjpEokmTiHv3GNkiqRXmrRokUDFi5cuH/L8tGjR9c/++yzm1qWB4NBu/DCC8edd955ZRdffPGu6G1z584tbe0cF198cfmpp546cdasWaV+v7+53Dln5513Xundd99dGF1/w4YNGffdd9/QN954Y+2gQYNCl1xyyZj6+vrmkbDMzEwH4Pf7XSgUMoAzzzxzfGlpafrkyZNr/vSnP20mRSipEhERSRId6akKh8NcdtllB4wbN65+/vz5bU6Qb2ncuHHB6667rvCMM86oii4//fTTK2fPnn3Qddddt2P48OGNxcXF/oqKCn9FRYU/Ozs7nJeXFyosLEx76623+s2YMaOqreMDvPTSS/mxxpNMlFSJiEivMHhwdvDOO4+J65IK8ToWwLe+9a2x77//fu6uXbvSxo8ff/j3v//9bVdddVXJvh7vtdde6/vCCy8MHDduXN20adMmAcyfP7/wG9/4RsXe9p03b94e5z3iiCPqf/CDHxTOmjXr4HA4TFpamrvrrru2nHTSSTUTJ06snTJlyqEjR44MHHHEEdUdjfUXv/jFkPvuu2//0tLS9BkzZkw68cQTK5KxB8ucc3uvlSAVFRVxC643PE1dbUwN7bYxXE9G7f0E+8whLfg2obTxuLQx3RpfPPT665gi4tXGfv36xW2Zg2gbNmx4ZOjQoW0+UkRkX+zYsWPtwQcffElr23T3n4iIiEgcaPhPRESkm0Svsh7thRdeWD9kyJBQImKS+FFSJSIi0k2iV1mX1KPhPxEREZE4UFIlIiIiEgdKqkRERETiQEmViIiISBxoorqIiPQKFioaZaHiuC3+6fyDg84/bGtb26Pv9CstLU33+XwuLy+vEeDNN99cm5WV1bwW46ZNm9Ivv/zysSUlJek+n4+LLrqo+IYbbtgZjzgff/zx/t/97nfH/fvf//708MMPr4/HMePpvffeyy4oKMiIZVHSaKeeeuoht99++9YZM2bUtiwvKCjIWLdu3Sc+n9d3NGvWrHHvvffefkVFRaviGPoelFSJiEivYKHijOyqm+K20Ghd7p0Zzj+sze3Rd/otWLBgeE5OTuimm25q9XEx6enp3HHHHQXHHnts7a5du3wnnHDCpDPOOKPyiCOO6HQS9Oyzzw74yle+Uv3EE08MOPzww7d19njxtnLlyj4ffvhhTkeTqvbk5uaGXnvttb6nnXZadWlpqb+4uDg9Xsduj4b/REREEmzUqFENxx57bC1A//79w2PHjq0rKCjodK9aZWWl78MPP+x77733frFkyZK8pvJ//OMfueecc85BTZ/nzp07euHChQMBnnvuuX6HH374oSeddNIhV1555aimegsWLBg+e/bsMV/72tfGT5gw4bDHHnus/zXXXDNyypQpk84666zxwWDQAN55550+p5xyyiFHH330xDPPPHP81q1b08HrQbr22mtHHHfccRMPO+ywya+88krf+vp6u+eee4YvXbo0b9q0aZMeeuihvKqqKt+ll1465thjj504ffr0SU8++WR/gJqaGrvwwgsPnDp16qQLLrjgwPr6+jYT5HPPPbfs6aefHgDw5JNP9j/zzDObHyZdWVnpO/300w+ePn36xClTpjQf/6233uozderUSbW1tVZVVeU78sgjD125cmVWR77fSqpERER6kPz8/Ix169b1mTFjxh7P0Fu0aNGAadOmTWr5Ov/88w9s7VhPPfVU/+OPP75i8uTJgX79+oXefffdPu2du7a21m688cYDnnnmmfw33nhjfWlp6W4jWlu3bs38+9//vvGRRx7ZeM0114w98cQTK1euXLkmMzMz/Nxzz/ULBoP2wx/+cPRjjz322fLly9dedNFFJT/60Y9GNO3f2Nho77zzztpbb711689//vPhWVlZ7tprr9121llnla9YsWLNpZdeWn7rrbcOO+GEEyrffffdtUuXLl1/xx13jKyqqvL99re/HZKdnR3+4IMP1vzwhz8sWrduXU5b7Tj11FOrVqxY0bexsZHnn39+wEUXXdT8EOrs7Ozw008/vfG9995b++KLL2647bbbRobDYY4//vja0047bdf8+fNHXH/99SNnzZpVOmXKlA71FGr4T0REpIeorKz0zZ49e9wtt9yyNS8vL9xy+5w5c8rmzJlT1tq+rfnLX/4y4Lvf/e5OgJkzZ5Y9+eSTA5p6xFqzevXqrBEjRgQOPvjgIMD5559f9vDDDw9u2n7yySdXZGRkuClTptSFw2GbNWtWJcCECRPqNm/enLF69erMzz77LHvmzJkHA4RCIQYPHtzQtP95551XDnD00UfX3Hzzza32xL311lv7vfrqq/3/+Mc/7g8QCARs06ZNGcuWLes7d+7cnQBHHXVU3UEHHdRmO/x+v5s6dWr1ww8/PKC+vt43fvz45odfO+ds/vz5I1esWNHX5/NRXFycsW3btrSRI0c23nbbbUXHH3/8xIyMjPC99967JZbvcTQlVSIiIj1AMBi0Cy+8cNx5551XdvHFF+9qrc6iRYsGLFy4cP+W5aNHj65/9tlnN0WX7dy5079ixYr9Nm7cmH399dcTDocNcPfcc09BWlqaC4e/zNkCgYABOOdoT0ZGhgPw+/2kpaW5pongPp+PxsZGc87ZgQceWPfWW2+ta23/zMzM5v1DoVCrw3fOOR566KGNhx12WKDlNrPYp8RdcMEFZXPmzDlo3rx5u80jW7x48YDS0tK0d999d21GRoabMGHCYXV1dT6AkpISf21tra+xsdHq6up8ubm5eyS27dHwn4iISIKFw2Euu+yyA8aNG1c/f/78Viezg9dTtWLFijUtXy0TKoAnnngi79xzzy1dv379J+vWrftkw4YNH48YMSL46quv9h07dmxg06ZN2XV1dVZWVuZftmzZfgCTJ0+uLywszMzPz88Ar6erI+049NBD68vLy9Nef/31HPASxVWrVrU7Lyk3NzdUXV3dnI+ccMIJlffee+/QpqRv+fLl2QDHHHNM9ZNPPjkAYOXKlVkbN25sdyjztNNOq77iiiuKZs+evVvPXmVlpX/QoEENGRkZ7qWXXsrdvn17c4/ZlVdeOeb666/fNmvWrNIf/OAHIzvSdlBPlYiI9BLOPzhYl3tnXJdUiNexXnvttb4vvPDCwHHjxtVNmzZtEsD8+fMLO3NH3F//+teB8+bNK4ouO/vss8uffPLJAffff/+Ws846q/zoo48+dPTo0fUTJkyoBcjJyXF33nnn5vPPP398Xl5e4+GHH17TkXNmZWW5P/3pT5/dcMMNo6urq/2hUMjmzJmz48gjj2xzbtIZZ5xR9b//+7/Dpk2bNunqq68uuvXWW7fNmzdv9FFHHTUJsOHDhweWLFmycd68eTsvu+yysVOnTp00YcKE2kMPPbTd2Hw+H63dbflf//VfZd/85jcPOuaYYyZOnDix9oADDqgHuO+++wb6/X737W9/u6yxsZETTzxxwtKlS3PPOuusqljbb3vr6kukioqKuAWXn5/P+PHj43W4HkltTA3ttjFcT0bt/QT7zCEt+DahtPG4tDHdGl889PrrmCLi1cZ+/frFbZmDaBs2bHhk6NChE7vi2KmssrLSt99++4XD4TBz584dfeCBB9bfeOONcVkzKxXs2LFj7cEHH3xJa9vUUyWSLMLV+ELbMbfHDUEiInHz+9//ftCzzz47qKGhwSZOnFh79dVXlyQ6pmShpEokSZirJ7PqNpx/8N4ri0iPFL3KerQXXnhh/ZAhQ0KJiKmlG2+8cad6pvaNkioREZFuEr3KuqQe3f0nIiKpqkO3w4vEqM2fKyVVIiKSkhoaGtYGAgH9nZO4CQQCvoaGhrVtbdfwn4iIpKT6+vq7iouLSU9Pn4g6EaTzwg0NDWvr6+vvaquCkioREUlJU6dODQM/S3Qc0nskPHM3s2vM7FMzW21mj5tZh54ILSIiItITJDSpMrMRwDzgKOfcZMAP/EciYxIRERHZFwnvqcIbgsw2szSgD7BtL/VFREREepyEJlXOuULgl8AWoAiocM69nMiYRERERPZFQp/9Z2Z5wLPAt4BdwNPAM865R2H3Z//l5+cnJEaRnmJQHvSvnw/+IRTbDeRlfkyAsRSVaBqiJJ/oZwZ21bP/RLpbou/+Ox343DlXDGBmzwHHAY+2rNjZh3bq4aapoTe30UIlZFX0wfmzGZI7hLRgNhlpufTNG9P9QXZSb76OqaQ3tFGkIxI9p2oLcIyZ9TEzA04D2lxUS0RERKSnSvScquXAM8BK4JNIPPclMiYRERGRfZHo4T+ccz8BfpLoOEREREQ6I9HDfyIiIiIpQUmViIiISBwoqRIRERGJAyVVIiIiInGgpEpEREQkDpRUiYiIiMSBkioRERGROFBSJSIiIhIHSqpERERE4kBJlYiIiEgcKKkSERERiYOYkiozK2ujfGd8wxERERFJTrH2VKW3LDCzdMAf33BEREREklNaexvN7E3AAVlm9u8Wm0cC73RVYCIiIiLJpN2kCngAMGAasCiq3AE7gFe7KC4RERGRpNJuUuWcewjAzJY559Z1T0giIiIiyafNpMrMLnHOPRL5eJyZHddaPefc4i6JTERERCSJtNdTdRHQlFRd0kYdByipEhERkV6vzaTKOXd21PtTuiccEWmNhYrxhT4HGhMdioiItGFvE9V3Y2ZDgL7RZc65TXGNSET25OrIrL7be5vgUEREpHUxJVVmdibe3X/DWmxyaK0qERERkZgX/7wXuB3Icc75ol5KqERERESIffgvD/ijc04jDyIiIiKtiLWnahFwWVcGIiIiIpLMYu2pOgaYZ2Y3AtujNzjnTox7VCIiIiJJJtak6oHIS0RERERaEVNS1fS4GhERERFpXaxLKny7rW16TI2IiIhI7MN/LR9Tsz8wDngbPaZGREREJObhvz0eUxPpvZoY94hEREREklCsSyq05kFgTpziEBEREUlqsc6papl89QFmA7viHpGIiIhIEop1TlUjez7HtRC4PL7hiIiIiCSnWJOqsS0+1zjnSuIdjIiIiEiyinWi+uauCsDM+uMtLDoZrzfs2865d7vqfCIiIiJdIdaeqq70G+Al59w3zSwDb76WiIiISFJJaFJlZvsBJwL/DeCcCwLBRMYkIiIisi9iWlKhlbv/4uVAoBj4k5mtMrMHzCyni84lIiIi0mXMuZY39bWoYOYHqoH+zrlAXE9udhSwDJjhnFtuZr8BKp1zNwNUVFQ0B5efnx/PU4skjcF5YbIzysiouNMrSBtNsd1AXubHBBhLUUlWYgMU2Qfjx49vft+vXz9LYCgicbPX4T/nXMjMNgADgW1xPn8BUOCcWx75/AxwY2sVo38B90V+fn6nj9HTqY2poWUbfQ1ryKr6NeR4nbhhfzZDcoeQFswmIy2XvnljEhTpvuuN1zEV9YY2inRErHOq/gwsifQkFRC1ZpVz7tV9PblzbruZbTWzQ5xz64HTgDX7ejwRERGRRIk1qfpu5OstLcod3ryozvh/wJ8jd/5tAi7r5PFEREREul2s61S1XPwzbpxzHwJHddXxRURERLpDzHf1mVm6mZ1gZt+KfM7RnXoiIiIinliXVDgM2ADcDyyKFJ8ELO6iuER6PV/DOixUmOgwREQkRrH2VP0B+LFzbgLQECl7Azi+S6ISEXyhreDCdKBDWUREEijWieqHAo9G3jsA51yNmWV3SVQiAkBW5U1E3WwrIiI9WKz/Bf4CmBpdYGbTgY3xDkhEvmQEMD25SUQkKcTaU3Uz8IKZLQQyzGw+MBe4vMsiExEREUkiMfVUOeeWAGcBg/HmUh0AfMM593IXxiYiIiKSNGLtqcI5txK4sgtjEREREUlasS6pkGFmt5lZvpnVRL7ebmZ6kquIiIgIsfdU/QE4BJgHbMYb/psPjAC+3TWhiYiIiCSPWJOq84Bxzrldkc9rzGw53t1/SqpERESk14t1SYXtQJ8WZdlAUXzDEREREUlOsfZUPQK8ZGa/AwqAUcBVwMNmdmpTJefcq/EPUURERKTnizWpuiLydUGL8rmRF3jLPh8Yj6BEREREkk1MSZVzbmxXByIiIiKSzPSkVhEREZE4UFIlIiIiEgdKqkRERETiQEmViIiISBzE+piaVW2Uvx/fcERERESSU6w9VQe1LDAzQ0soiIiIiAB7WVLBzB6OvM2Iet9kDPBpVwQlIiIikmz2tk7VZ228d8DbwNNxj0ikN3MNjBwSgHB9oiMREZEOajepcs7dCmBmy5xz/+iekER6MRciq+4BMq0vvtD2REcjIiIdEOuK6v8ws0OArwB9W2xb3BWBifReYfyN+YkOQkREOiimpMrMFgA/Bj4CaqM2OUBJlYiIiPR6sT5Q+fvAdOfcx10ZjIiIiEiyinVJhTpgXVcGIiIiIpLMYk2qbgZ+Z2bDzMwX/erK4ERERPm6kqcAABZ3SURBVESSRazDfw9Gvn4nqszw5lT54xmQiIiISDKKNaka26VRiIiIiCS5WJdU2AwQGe4b6pwr6tKoRERERJJMrA9U7m9mjwH1wMZI2Uwzu6MrgxMRERFJFrFONF8IVAAHAMFI2bvAt7oiKBEREZFkE+ucqtOA4c65BjNzAM65YjMbEo8gzMwPvA8UOufOiccxRURERLpTrD1VFcCg6AIzGw3Ea27V94C1cTqWiIiISLeLNal6AHjWzE4BfGZ2LPAQ3rBgp5jZSODrkXOI9F7hXfgb3ofQzkRHIiIi+yDW4b+f401SvxdIx3ve3x+B38Qhhl8DNwC5cTiWSNIy10hGzX00hquBnESHIyIiHWTOucSd3Owc4Gzn3JVmdjJwffScqoqKiubg8vPzExChSPcZlAf96+eDq2+/Ytpoiu0G8jI/JsBYikqyuidAkTgaP3588/t+/fpZAkMRiZuYeqrM7EbgX865FVFl04GTnXN3deL8M4CZZnY2kAXsZ2aPOudmt6wY/Qu4L/Lz8zt9jJ5ObUxuFiohq6IPtTX15OS03VMV9mczJHcIacFsMtJy6Zs3pvuCjJNUvo5N1EaR3ifWOVXfA9a0KFsDfL8zJ3fOzXfOjXTOjQH+A3i1tYRKREREpKeLNanKABpalAXxepdEREREer1Yk6oPgCtblM0FVsYrEOfc61qjSkRERJJVrHf/XQP808wuAT4DDgKGAl/tqsBEREREkslekyozM6AOOBg4BxgFPAcscc5Vd214IiIiIslhr0mVc86Z2SdArnPuiW6ISURERCTpxDqnahVeT5WIdIVwGbhdQCjRkYiIyD6KdU7V68BLZvYgsBVoXpTTObc4/mGJ9C4Wria78qZEhyEiIp0Qa1I1A/gcOKlFucN7ZI2IiIhIrxZTUuWcO6WrAxERERFJZrHOqcLMBprZJWb2g8jn4WY2sutCExEREUkeMSVVZnYSsB74T+DHkeLxwB+6KC4RERGRpBJrT9WvgW85584EGiNly4HpXRKViIiISJKJNaka45z7V+R9051/QWKf6C4iIiKS0mJNqtaY2ddalJ0OfBLneERERESSUqw9TdcBS8zsBSDbzP4InAvM6rLIRERERJJITD1VzrllwOHAp3jrUn0OTHfOrejC2ERERESSRrs9VWbWB/gRMBlYCfyPcy7QHYGJiIiIJJO99VT9L94w3zrgm8AvuzwiERERkSS0t6TqLOAM59wNkffndH1IIr1MuAbMn+goRESkk/Y2UT3HOVcE4Jzbamb9uiEmkV7FF9pKZvVdiQ5DREQ6aW9JVZqZnQJYG59xzr3aVcGJpCp/cAVh/4E4GrBwCeZqEh2SiIh00t6Sqp14d/s1KW3x2QEHxjsokVTna/iYsH8MvtAWMmt+l+hwREQkDtpNqpxzY7opDhEREZGkFuuK6iIiIiLSDiVVIiIiInGgpEpEREQkDpRUiYiIiMSBkioRERGROFBSJZIQDghhrj7RgYiISJzsbZ0qEekKLkR2xfcSHYWIiMSReqpERERE4kBJlYiIiEgcKKkSERERiQMlVSIiIiJxoInqIt3JNWLhcszVJjoSERGJMyVVIt3JNZJZ/Ut8oS8SHYmIiMRZQof/zGyUmb1mZmvN7FMz0z3mIiIikpQS3VPVCFznnFtpZrnAB2b2T+fcmgTHJSIiItIhCe2pcs4VOedWRt5XAWuBEYmMSURERGRfmHMu0TEAYGZjgH8Dk51zlQAVFRXNweXn5ycmMJE46puTzv6+e6Bh074fJG00xXYDeZkfE2AsRSVZ8QtQpJuMHz+++X2/fv0sgaGIxE2ih/8AMLO+wLPA95sSqpaifwH3RX5+fqeP0dOpjT1cuA5z1WRW+fFl5LRZraamhpyctreH/dkMyR1CWjCbjLRc+uaN6YJgu1ZSX8cYqY0ivU/CkyozS8dLqP7snHsu0fGIdAnngCBZFddj6CHKIiKpKKFJlZkZsAhY65y7J5GxiHQlCxWQUXsfdCChamx0+Hzg82lkREQkGSR6RfUZwCXAqWb2YeR1doJjEukCDn/jBjqSHlVWBgkEQl0WkYiIxFdCe6qcc29Bh/7OiCQdCxVj4YoO7xcKOXrIfSQiIhKDRPdUiaQ8C5eQVX1HosMQEZEupqRKREREJA6UVIl0k4aGMLt2BRMdhoiIdBElVSLdxDmoqFBSJSKSqpRUiXShgoJq6uoaO7SPt9KIiIgkGyVVIl1o27YaGhq+XBahujrIzp11hEKt39YXDkNdnbW5XUREei4lVSLdqKEhTHFxXfPnUMi1SKCc5l2JiCSphD+mRiTVlZcHSCOAc45wePceqOrqBvx+o2/f9H06dm1tI2VV1YwY0TceoYqISCcoqRKJo0AgxKpVxXzlKwPJzvYSpaLttVTXffmc8HDYEQiEyMjw45w3gT0QCFFaWo9zUF8f+xys2tpGyssDSqpERHoADf+JxJFzjuXLdxAMhgkGW3/ETDAYZt268t16rcJhKCqqZfv22j16s3Y/vvdMwI8+LiHYEI57/CIisu+UVInEWVFRDTfdtKxDyye0d8NfQ0OYhkgC5ZyjsLCG99/fuUe9urpGduyo3a2soKCawsLqmOMQEZF9p6RKZB9s3VpFdfWXSVN5eT3bttUAXm9SWVkg5mPt3Fm3WzIUCITYvr2WkpI6GhvDVFU1UF3dQENDmB076qitbSDcSidVTU0Dzz//+W5l27bVUF4eeywiIrLvlFSJ7IN//GMrdXUhamoaqK1toKqqgaee2sjGjRUdXuBz164AJSX1zZ/DYcf27bXU1jaSn1/B559787GcIzI86PVYBQMhqqpiO1d9feNuSaCIiMSfJqqLxKiiIsDWrdUMH57TvKDnhg27SE/3kZubzrvvbufdd7cD4Pd743kj9g9QV5lGQ12bh43JF19U4fdb83wr52BrQTXby6rJG7L3/YuL68jPr+DUU0d2LhAREWmTkiqRNmzcuIs+fdIJhcI0Njr22y+dX/xiFfX1Iczg7LMPYNu2Gh5+eP0e+5pBv5xSQm4lX6z/7R7bnYOdO2tjXm09HN5zOYYmGRl773BueVfhli1V5OVlkpubEdP5RURk7zT8J9KGL76oIhgMUV3dsMcEcOdgwYJluyVUp568H7f9aACDB2cDEK79hPodD7R5/B076qisbOh0nH/4w6esXl26W1llZWCPuw/z8yvYsGEXNTUNLF26hWAwvNflG1avLmXDhl2djlFEpDdQT5VIDMzA52v/mXxDh6aRFXqNIUNO4Rszcyks+JBgcW27+8TDtm01lJbWM3hwNs55vVkvvLCZ008fBcCqVcX075/JO+9sZ9myHeTlZbJrV4CVK3eyYMFU0tP9hMOOkSP3XOuqtrZRzyIUEYmRkirplbZuraJ//92Hv2pqGigtrWf06NzmsowMP42NYZYu3cL69bsIBHbv/bn8soFU14R4/CmvN8c1lvONczPIDj5KsPiD7mkM8Pvfrwbg+OP3B7y2NFmyZHPz+3DYUVrqTYrftSuIc95E+aqqYKtJlYiIxE5JlfRKL7+8lZkzx5Kb682dqq8PMWxYH5Yu3cwVV0ymqipIQ0OY3/72IxoawhQU1LB2bfkex9l/YDGBqMQsUPkx/tByEnWf3bZtNaxZU8aWLdVs3VrNxo0dG7qrr2/kkUfWU1JSz3e+M6mLohQRSU1KqqTXa2gI8+qrBRx22MDmRTb//vcvePHFzc2fk8WmTVXceusKAH7+85X7dIz8/AotGCoisg+UVIkAb7+9nbff3k7//hk8+WQ+b765LekSKhERSSwlVdJrlJbWs2jRGgKBEJs3VzFz5tg96uzaFeS55zZ16Lh9swq54ZrD2Lot9gch9xRbtlRTXR1k4MCsDu1XWRlk7doyhg3L2W0OmohIb6YlFaRX2Lmzlo8/LuHDD0tYvbqM2lovASorqycUavsBxrEIlz1IdnoRWVnp8Qi1W9177yc89JC3LMSuXQE++KC4eSJ7k5UriykoqKaurpHt22spLq7jiSfy+d3vPsF17lsnIpJS1FMlvUIgEGLhwk/3KP/sswp++csPO3w8M/j+1UPI8H9KCEgPvMT4/ccTLE3O5+z5/YYZ/PGPnxIIhEhL+3IZhSVLvqB//wwmTRrAj3+8nMZGZVIiIq1RUiUpq6CgmowMP5mZPvz++HbK+nzGwL4baSxeDEBjzWqoWR3Xc3SnBx9cx8CBWXssGdHk0Uc3tLt/INBIQ0OYvn33XKG9rq6RcNiRk5N8PXkiIh2hpEpSRkVFgB07ajn44DwA1q8v5/HH88nI8O8x6TwcdnzySSnr1++5TEJvVFxcT3Hxl8N+oZBj9epS1qzZ+/dn27YaVq0qJicnjZNP3vPZglu2VFFT08CUKTE8pFBEJIkpqZKktmNHLT6fMXhwNo2NjtdeK2xOqgCqqhqAPR8F45w31CWtc45Wh0tbs2tXgIcfXs9///eEVreHQk53UopIr6CkSpJCINBIMLjnH+YtW6rIzk5rft7e9u21rF5dSmOjY/36+D6z7tST9+OkGT4+XpPBG2/qeXjtqaoKkpnpJyPD31wWDIYIBEJ6iLOIpCwlVZIUiovrWb++nDFjvpwbVVpaRyAQJjsbqquDhEJhPvuskttvf79LYhg6NA0rvZ2jxx9M+a5vASVdcp5k8tFHJVRVeevHv//+TsaOzWXChAEsWfIFp58+qjnZBS/RevnlrVx00cHNZWVl9fh8Rv/+md0eu4hIvCmpkh5t585a/H4jPd1HXd2X60Bt3VrFG29s4+9//4Kbbz6K8vIA11//TtcH5EI0Bkq7/jxJ4s9//nIC++rVZeTlZTJhwgCqqhqorAw230VYXFzH4MHZVFR4CVjTQ5o3bqwgJyddSZWIpAQlVdJjlJbWs2zZdr7+9TEArFpVzJ/+tJa0NB9ZWWkcd9zQ5rpLl27hX/8q6PKYrp03mNraEAsfKGsuCwUrOOmYOlxgfZefP1mFQmEWLFjG3XfPwO83HnnES76ys/2YwRln9N+tfiDQSHV1427XX0Qk2Sipkh5l69YvnzlXVlbPjh11zZ+zsnxkZg5k+/YiNmz48q40n8/w+YzOmnBINifOyOGxpyrw+43/M7MfA/tuI43s3eqFGnZB8S20vvhA77ZpUyXvvLOd/PwKAG655T3C4S/XtaqrC/Hpp+Wcddaw3Xoed+yo46c//YDJkwe0etyqqiDBYIiBA7Nb3S4i0hMoqZKE2rGjlrQ0Y+DAbNLTv5wvtXVrNUVFtbvV/fTTct57r4CcnJzdypcu3cx++3Vu8vMZp+dy4nFZpAeWcvjh53D+18upLfo9DTs30XfQ97liznAqq7Xo5d4UFtbwm9981PzZu/tydzt21PK9773Jgw+ehd//ZTJcXh6gvDzAli1VDBvWh/R0P9u315KebgSDYT78sISzzjqgW9ohIrIvlFRJl6iqClJeHtjtzryysnrefdcb3mm6O2zLlioefzyfAQMyCYUc/fp5ydFLL23mlVdiG957772dMcc14ZC+/OeF6Tz910Y+/qSqufygsT4at10PfQ7iqydWUb/zGUL13jMAwyW/Zv+B11FZrXWW4qWhIcyvfvUhpaUBrrrqsObkavXqMn7964/4yU+mkZ3t3d35xBP5ZGX5OeGE4c37v/jiF0yfPpQdO2rJyPAzfnz/tk4lItJtEp5UmdmZwG8AP/CAc+5nCQ5JOqCtu7eWLPmCF1/czO23H83gwdmUlNSxenUpn33mDQu98MJmTjllBMFgmMLCGgoLawDYf/9s3nlnOxs2xG/Jgq+eth8D8tJ48pky+vdPI8f9kwEDTgfg/PPyGDjAx8hRmQS3QGPtRqz2x62sbCXxtnatd43vumvVbj1WhYU13HTTMs4/fxxr1pQ3/2xEJ1VNZdXVjfTt6+1bV9fYPBF+0KAs0tL0aFMR6V4JTarMzA/cC3wVKABWmNnfnHNrEhmXtK2xMUxhYTXDh/eltLSOVatKGDmyL4FAiE2bKhgxoi+jR+dSXd1AMBgmPd1HUVENL7ywmX/+cyvHHONNNq+qCjJv3pt7HH/79rrdho864qgpOZxyYiaLH6ll0KAsiopqqawMcvjEAHk5Gxh+9WFsLcomHKrjpBkZ9O07hMmHlGCl/0NwS6e+LdIJNTV7prDFxfV7LD766qsFDByYxfTp3s9Qerqv+c7QiooApaX13Hzzcvx+HzNnjuWMM0ZRXh4gPd3H8OE5e5wjWnl5PfX1IYYNa7tedbXX+zpqVO4+tFJEegNzCXzMvJkdC9zinPta5PN8AOfc/wBUVFRoEouISIrr169f5+80EekBEt0/PgLYGvW5IFImIiIiklQSnVS19r8T9U6JiIhI0kn0RPUCYFTU55HAtqYP6hIWERGRZJHonqoVwHgzG2tmGcB/AH9LcEwiIiIiHZbQpMo51whcDfwDWAs85Zz7tP299mRmi81sp5mtjiq7xcwKzezDyOvsqG3zzWyjma03s6/Foy1dzcxGmdlrZrbWzD41s+9FygeY2T/NLD/yNS9qn6RqZzttTJlraWZZZvaemX0UaeOtkfJUuo5ttTFlrmMTM/Ob2SozWxL5nDLXsUkrbUy56ygSLwm9+y9ezOxEoBp42Dk3OVJ2C1DtnPtli7qTgMeB6cBw4BXgYOdcj37qiJkNA4Y551aaWS7wAXAe8N9AmXPuZ2Z2I5DnnPthMraznTZeSIpcS/OeJJzjnKs2s3TgLeB7wDdInevYVhvPJEWuYxMzuxY4CtjPOXeOmd1FilzHJq208RZS7DqKxEuih//iwjn3b6BsrxU9s4AnnHMB59znwEa8fwR6NOdckXNuZeR9FV7P3gi89jwUqfYQXhICSdjOdtrYlmRso3PONT3gMD3ycqTWdWyrjW1JujYCmNlI4OvAA1HFKXMdoc02tiUp2ygSTymRVLXjajP7ODI82NQNn/TLOJjZGOBIYDkw1DlXBF5SAjQ9SyWp29mijZBC1zIynPIhsBP4p3Mu5a5jG22EFLqOwK+BG4BwVFlKXUdabyOk1nUUiZtUTqr+AIwDjgCKgLsj5Um9jIOZ9QWeBb7vnKtsr2orZUnRzlbamFLX0jkXcs4dgXe363Qzm9xO9VRqY8pcRzM7B9jpnPsg1l1aKUvWNqbMdRSJt5RNqpxzOyL/sIeB+/myG7rdZRx6ssj8lGeBPzvnnosU74jMRWqak9T0dOGkbGdrbUzFawngnNsFvI431yilrmOT6Dam2HWcAcw0sy+AJ4BTzexRUus6ttrGFLuOInGVsklV0z9sEf8HaLoz8G/Af5hZppmNBcYD73V3fB0Vmfy7CFjrnLsnatPfgEsj7y8Fno8qT6p2ttXGVLqWZjbYzPpH3mcDpwPrSK3r2GobU+k6OufmO+dGOufG4C0F86pzbjYpdB3bamMqXUeReEv04p9xYWaPAycDg8ysAPgJcLKZHYHX/fwFcAWAc+5TM3sKWAM0Alclyd0pM4BLgE8ic1UAFgA/A54ysznAFuACSNp2ttXGi1LoWg4DHjLvYeI+vGVElpjZu6TOdWyrjY+k0HVsSyr9Prblrl5wHUX2SUosqSAiIiKSaCk7/CciIiLSnZRUiYiIiMSBkioRERGROFBSJSIiIhIHSqpERERE4kBJlYiIiEgcKKkSERERiQMlVSIiIiJx8P8BAf70ZzXQzaIAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"compare_T1_T2(300, 30, 5000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can see why $T_2$ is a better estimator than $T_1$. \n",
"\n",
"- Both are unbiased. So both the empirical histograms are balanced at around $300$, the true value of $N$.\n",
"- The emipirical distribution of $T_2$ is clustered much closer to the true value $300$ than the empirical distribution of $T_1$.\n",
"\n",
"For a recap, take another look at the [accuracy table](https://en.wikipedia.org/wiki/German_tank_problem#Specific_data) of the Allied statisticians' estimator $T_2$. Not bad for an estimator based on a model that assumes nothing more complicated than simple random sampling!"
]
}
],
"metadata": {
"anaconda-cloud": {},
"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.7.11"
}
},
"nbformat": 4,
"nbformat_minor": 4
}