From aa9a60fd95755136e988566dbff90ffc740c3e93 Mon Sep 17 00:00:00 2001 From: George Ho <19851673+eigenfoo@users.noreply.github.com> Date: Wed, 26 Jun 2019 00:47:21 +0000 Subject: [PATCH] BLD: begin porting ANOVA chapter (#11) * MAINT: update toc function * BLD: building * BUG: add d as a good name to pylint * MAINT: remove patsy as req --- .pylintrc | 2 +- Makefile | 2 +- index.html | 346 ++++++++++++++++++++++++++++++++++-------- plots.py | 51 +++++++ requirements-dev.txt | 5 + requirements.txt | 5 - tests-as-linear.ipynb | 297 ++++++++++++++++++++++++++++++------ utils.py | 7 +- 8 files changed, 602 insertions(+), 113 deletions(-) create mode 100644 requirements-dev.txt diff --git a/.pylintrc b/.pylintrc index 399caca..7c70d31 100644 --- a/.pylintrc +++ b/.pylintrc @@ -94,7 +94,7 @@ evaluation=10.0 - ((float(5 * error + warning + refactor + convention) / stateme bad-functions=map,filter,input # Good variable names which should always be accepted, separated by a comma -good-names=a,b,c,f,i,j,k,df,x,y,y2,_,fig,ax +good-names=a,b,c,d,f,i,j,k,df,x,y,y2,_,fig,ax # Bad variable names which should always be refused, separated by a comma bad-names=foo,bar,baz,toto,tutu,tata diff --git a/Makefile b/Makefile index 3f18912..7f342eb 100644 --- a/Makefile +++ b/Makefile @@ -16,7 +16,7 @@ venv: # Set up Python virtual environment. python -m venv ${VENV_PATH}; \ source ${VENV_PATH}/bin/activate; \ pip install -U pip; \ - pip install -r requirements.txt; \ + pip install -r requirements-dev.txt; \ deactivate; \ ) @printf "\n\nVirtual environment created! \033[1;34mRun \`source ${VENV_PATH}/bin/activate\` to activate it.\033[0m\n\n\n" diff --git a/index.html b/index.html index baae204..6e4e239 100644 --- a/index.html +++ b/index.html @@ -13187,7 +13187,7 @@

Com
-

Last updated: June 25, 2019

+

Last updated: June 26, 2019

@@ -13203,74 +13203,33 @@

Contents

Common statistical tests are linear models (or: how to teach stats))
  • 1 The simplicity underlying common tests
  • 2 Settings and toy data
  • -
  • 3 Pearson and Spearman correlation -
  • +
  • 3 Pearson and Spearman correlation
  • 4 One mean
  • 5 Two means
  • 6 Three or more means
  • 7 Proportions: Chi-square is a log-linear model
  • 8 Sources and further equivalences
  • 9 Teaching materials and a course outline
  • -
  • 10 Limitations
  • +
  • 10 Limitations +
  • @@ -14755,7 +14714,7 @@

    6 Three or more means6.1 One-way ANOVA and Kruskal-Wallis

    6.1.1 Theory: As linear models

    Model: One mean for each group predicts $y$.

    $y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 +... \qquad \mathcal{H}_0: y = \beta_0$

    where $x_i$ are indicators ($x=0$ or $x=1$) where at most one $x_i=1$ while all others are $x_i=0$.

    -

    Notice how this is just "more of the same" of what we already did in other models above. When there are only two groups, this model is $y = \beta_0 + \beta_1*x$, i.e. the independent t-test. If there is only one group, it is $y = \beta_0$, i.e. the one-sample t-test. This is easy to see in the visualization below - just cover up a few groups and see that it matches the other visualizations above.

    +

    Notice how this is just "more of the same" of what we already did in other models above. When there are only two groups, this model is $y = \beta_0 + \beta_1*x$, i.e. the independent t-test. If there is only one group, it is $y = \beta_0$, i.e. the one-sample t-test. This is easy to see in the visualization below - just cover up a few groups and see that it matches the other visualizations above.

    @@ -14769,19 +14728,41 @@

    6.1 One-way ANOVA and Kruskal-Wall
    -
     
    +
    plots.one_way_anova_plot()
    +plt.show()
     
    +
    +
    + +
    +
    + + +
    + + + + + +
    + +
    + +
    +
    -

    A one-way ANOVA has a log-linear counterpart called goodness-of-fit test which we'll return to. By the way, since we now regress on more than one $x$, the one-way ANOVA is a multiple regression model.

    +

    A one-way ANOVA has a log-linear counterpart called goodness-of-fit test which we'll return to. By the way, since we now regress on more than one $x$, the one-way ANOVA is a multiple regression model.

    The Kruskal-Wallis test is simply a one-way ANOVA on the rank-transformed $y$ (value):

    -

    $rank(y) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 +...$

    +

    $\text{rank}(y) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 +...$

    This approximation is good enough for 12 or more data points. Again, if you do this for just one or two groups, we're already acquainted with those equations, i.e. the Wilcoxon signed-rank test or the Mann-Whitney U test respectively.

    @@ -14803,10 +14784,116 @@

    6.1.2 Example data
    -
     
    +
    num_points = 20
    +df = pd.DataFrame()
    +df["y"] = np.concatenate([
    +    np.random.normal(0.0, 1, num_points),
    +    np.random.normal(1.0, 1, num_points),
    +    np.random.normal(0.5, 1, num_points),
    +])
    +
    +df["group"] = list("".join([num_points * char for char in "abc"]))
    +df = df.join(pd.get_dummies(df.group, prefix="group", drop_first=True).astype(np.float64))
     
    +
    +

    + +
    + + + +
    +
    + +
    +
    +
    df.head()
    +
    + +
    +
    +
    + +
    +
    + + +
    + + + + +
    +
    + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
    ygroupgroup_bgroup_c
    01.080830a0.00.0
    1-2.316629a0.00.0
    2-0.138365a0.00.0
    32.003775a0.00.0
    42.051200a0.00.0
    +
    +
    + +
    +
    @@ -14834,10 +14921,97 @@

    6.1.3 Python code: one-way ANOVA
    -
     
    +
    F, p = scipy.stats.f_oneway(df[df["group"] == "a"].y,
    +                            df[df["group"] == "b"].y,
    +                            df[df["group"] == "c"].y)
    +
    +res = smf.ols("y ~ 1 + group_b + group_c", df).fit()
     
    +
    +

    + + + + + +
    +
    + +
    +
    +
    # FIXME what to tabulate here?
    +utils.tabulate_results([None, p, None, None, None],
    +                       res,
    +                       ["scipy.stats.f_oneway", "smf.ols (y ~ 1 + group_b + group_c)"],
    +                       coeff="group_b")
    +
    + +
    +
    +
    + +
    +
    + + +
    + + + + +
    +
    + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
    valuep-valuest-values0.025 CI0.975 CI
    scipy.stats.f_onewayNaN0.015156NaNNaNNaN
    smf.ols (y ~ 1 + group_b + group_c)0.8018630.0121262.5911780.1821821.421545
    +
    +
    + +
    +
    @@ -14866,7 +15040,11 @@

    6.1.4 Python code: Kruskal-Wallis + + + +
    +
    + +
    +
    +
    !cat requirements.txt
    +
    + +
    +
    +
    + +
    +
    + + +
    + + + +
    +
    jupyter==1.0.0
    +matplotlib==3.1.0
    +numpy==1.16.4
    +pandas==0.24.2
    +patsy==0.5.1
    +scipy==1.3.0
    +statsmodels==0.10.0
    +
    +
    +
    + +
    +
    + +
    diff --git a/plots.py b/plots.py index 6b45fba..a9f0aba 100644 --- a/plots.py +++ b/plots.py @@ -182,3 +182,54 @@ def dummy_coding_plot(): ax.legend(fontsize="large") return fig, ax + + +def one_way_anova_plot(): + a = np.random.normal(0, 1, 20) + b = np.random.normal(-2, 1, 20) + c = np.random.normal(3, 1, 20) + d = np.random.normal(1.5, 1, 20) + + df = pd.DataFrame() + df["y"] = np.concatenate([a, b, c, d]) + df["group_2"] = np.concatenate( + [np.zeros_like(b)] + [np.ones_like(b)] + 2 * [np.zeros_like(b)] + ) + df["group_3"] = np.concatenate( + 2 * [np.zeros_like(c)] + [np.ones_like(c)] + [np.zeros_like(c)] + ) + df["group_4"] = np.concatenate(3 * [np.zeros_like(d)] + [np.ones_like(d)]) + + res = smf.ols("y ~ 1 + group_2 + group_3 + group_4", df).fit() + beta0, beta1, beta2, beta3 = res.params + + fig, ax = plt.subplots(figsize=[10, 8]) + ax.scatter(0 * np.ones_like(a), a, color="k") + ax.scatter(1 * np.ones_like(b), b, color="k") + ax.scatter(2 * np.ones_like(c), c, color="k") + ax.scatter(3 * np.ones_like(d), d, color="k") + + ax.axhline(beta0, color="b", label=r"$\beta_0$ (group 1 mean)") + + ax.plot([0.7, 1.3], 2 * [beta0 + beta1], color="navy") + ax.plot( + [0, 1], + [beta0, beta0 + beta1], + color="r", + label=r"$\beta_1, \beta_2, ...$ (slopes/differences to $\beta_0$)", + ) + + ax.plot( + [1.7, 2.3], + 2 * [beta0 + beta2], + color="navy", + label=r"$\beta_0+\beta_1, \beta_0+\beta_2 ...$ (group 2, 3 ... means)", + ) + ax.plot([1, 2], [beta0, beta0 + beta2], color="r") + + ax.plot([2.7, 3.3], 2 * [beta0 + beta3], color="navy") + ax.plot([2, 3], [beta0, beta0 + beta3], color="r") + + ax.legend(fontsize="large") + + return fig, ax diff --git a/requirements-dev.txt b/requirements-dev.txt new file mode 100644 index 0000000..ccac16d --- /dev/null +++ b/requirements-dev.txt @@ -0,0 +1,5 @@ +-r requirements.txt +black==19.3b0 +nbdime==1.0.6 +nbinteract==0.2.4 +pylint==2.3.1 diff --git a/requirements.txt b/requirements.txt index 487e94f..48e8d1b 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,11 +1,6 @@ -black==19.3b0 jupyter==1.0.0 matplotlib==3.1.0 -nbdime==1.0.6 -nbinteract==0.2.4 numpy==1.16.4 pandas==0.24.2 -patsy==0.5.1 -pylint==2.3.1 scipy==1.3.0 statsmodels==0.10.0 diff --git a/tests-as-linear.ipynb b/tests-as-linear.ipynb index 6a50d7d..2a184c6 100644 --- a/tests-as-linear.ipynb +++ b/tests-as-linear.ipynb @@ -17,7 +17,7 @@ { "data": { "text/markdown": [ - "Last updated: June 25, 2019" + "Last updated: June 26, 2019" ], "text/plain": [ "" @@ -34,49 +34,22 @@ "- [1 The simplicity underlying common tests](#1-The-simplicity-underlying-common-tests)\n", "- [2 Settings and toy data](#2-Settings-and-toy-data)\n", "- [3 Pearson and Spearman correlation](#3-Pearson-and-Spearman-correlation)\n", - " - [3.0.1 Theory: As linear models](#3.0.1-Theory:-As-linear-models)\n", - " - [3.0.2 Theory: rank-transformation](#3.0.2-Theory:-rank-transformation)\n", - " - [3.0.3 Python code: Pearson correlation](#3.0.3-Python-code:-Pearson-correlation)\n", - " - [3.0.4 Python code: Spearman correlation](#3.0.4-Python-code:-Spearman-correlation)\n", "- [4 One mean](#4-One-mean)\n", " - [4.1 One sample t-test and Wilcoxon signed-rank](#4.1-One-sample-t-test-and-Wilcoxon-signed-rank)\n", - " - [4.1.1 Theory: As linear models](#4.1.1-Theory:-As-linear-models)\n", - " - [4.1.2 Python code: One-sample t-test](#4.1.2-Python-code:-One-sample-t-test)\n", - " - [4.1.3 Python code: Wilcoxon signed-rank test](#4.1.3-Python-code:-Wilcoxon-signed-rank-test)\n", " - [4.2 Paired samples t-test and Wilcoxon matched pairs](#4.2-Paired-samples-t-test-and-Wilcoxon-matched-pairs)\n", - " - [4.2.1 Theory: As linear models](#4.2.1-Theory:-As-linear-models)\n", - " - [4.2.2 Python code: Paired sample t-test](#4.2.2-Python-code:-Paired-sample-t-test)\n", - " - [4.2.3 Python code: Wilcoxon matched pairs](#4.2.3-Python-code:-Wilcoxon-matched-pairs)\n", "- [5 Two means](#5-Two-means)\n", " - [5.1 Independent t-test and Mann-Whitney U](#5.1-Independent-t-test-and-Mann-Whitney-U)\n", - " - [5.1.1 Theory: As linear models](#5.1.1-Theory:-As-linear-models)\n", - " - [5.1.2 Theory: Dummy coding](#5.1.2-Theory:-Dummy-coding)\n", - " - [5.1.3 Theory: Dummy coding (continued)](#5.1.3-Theory:-Dummy-coding-(continued))\n", - " - [5.1.4 Python code: independent t-test](#5.1.4-Python-code:-independent-t-test)\n", - " - [5.1.5 Python code: Mann-Whitney U](#5.1.5-Python-code:-Mann-Whitney-U)\n", " - [5.2 Welch’s t-test](#5.2-Welch’s-t-test)\n", "- [6 Three or more means](#6-Three-or-more-means)\n", " - [6.1 One-way ANOVA and Kruskal-Wallis](#6.1-One-way-ANOVA-and-Kruskal-Wallis)\n", - " - [6.1.1 Theory: As linear models](#6.1.1-Theory:-As-linear-models)\n", - " - [6.1.2 Example data](#6.1.2-Example-data)\n", - " - [6.1.3 Python code: one-way ANOVA](#6.1.3-Python-code:-one-way-ANOVA)\n", - " - [6.1.4 Python code: Kruskal-Wallis](#6.1.4-Python-code:-Kruskal-Wallis)\n", " - [6.2 Two-way ANOVA](#6.2-Two-way-ANOVA)\n", - " - [6.2.1 Theory: As linear models](#6.2.1-Theory:-As-linear-models)\n", - " - [6.2.2 Python code: Two-way ANOVA](#6.2.2-Python-code:-Two-way-ANOVA)\n", - " - [6.3 ANCOVA](#6.3-ANCOVA)\n", "- [7 Proportions: Chi-square is a log-linear model](#7-Proportions:-Chi-square-is-a-log-linear-model)\n", " - [7.1 Goodness of fit](#7.1-Goodness-of-fit)\n", - " - [7.1.1 Theory: As log-linear model](#7.1.1-Theory:-As-log-linear-model)\n", - " - [7.1.2 Example data](#7.1.2-Example-data)\n", - " - [7.1.3 Python code: Goodness of fit](#7.1.3-Python-code:-Goodness-of-fit)\n", " - [7.2 Contingency tables](#7.2-Contingency-tables)\n", - " - [7.2.1 Theory: As log-linear model](#7.2.1-Theory:-As-log-linear-model)\n", - " - [7.2.2 Example data](#7.2.2-Example-data)\n", - " - [7.2.3 Python code: Chi-square test](#7.2.3-Python-code:-Chi-square-test)\n", "- [8 Sources and further equivalences](#8-Sources-and-further-equivalences)\n", "- [9 Teaching materials and a course outline](#9-Teaching-materials-and-a-course-outline)\n", - "- [10 Limitations](#10-Limitations)" + "- [10 Limitations](#10-Limitations)\n", + " - [11 Computing Environment](#11-Computing-Environment)" ], "text/plain": [ "" @@ -626,7 +599,9 @@ { "cell_type": "code", "execution_count": 13, - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [ { "data": { @@ -1412,25 +1387,41 @@ "\n", "where $x_i$ are indicators ($x=0$ or $x=1$) where at most one $x_i=1$ while all others are $x_i=0$. \n", "\n", - "Notice how this is just \"more of the same\" of what we already did in other models above. When there are only two groups, this model is $y = \\beta_0 + \\beta_1*x$, i.e. the [independent t-test](#5.1.4-Python-code:-independent-t-test). If there is only one group, it is $y = \\beta_0$, i.e. the [one-sample t-test](#4.1-One-sample-t-test-and-Wilcoxon-signed-rank). This is easy to see in the visualization below - just cover up a few groups and see that it matches the other visualizations above." + "Notice how this is just \"more of the same\" of what we already did in other models above. When there are only two groups, this model is $y = \\beta_0 + \\beta_1*x$, i.e. the [independent t-test](#5.1-Independent-t-test-and-Mann-Whitney-U). If there is only one group, it is $y = \\beta_0$, i.e. the [one-sample t-test](#4.1-One-sample-t-test-and-Wilcoxon-signed-rank). This is easy to see in the visualization below - just cover up a few groups and see that it matches the other visualizations above." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 28, "metadata": {}, - "outputs": [], - "source": [] + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAHNCAYAAAAOvD9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdf3zO9f7H8cc200yS30l2zZEfo04ZnRT6zRH6hVIujYZVOnJOh4OGTljkR5KOZTlquET5FY4KOSGxMhXZj5BdfnUoyq/FbPt8/3i3fY1tdu3XZ9f2vN9uu2XvXdfn87o+Xex5vd/vz/vtY1kWIiIiIpI7X7sLEBERESnLFJZERERE8qGwJCIiIpIPhSURERGRfCgsiYiIiOSjUkkdOD4+XrfZiYiIiNdo3bq1T27tJRaWfj9pSR4egMTEREJCQkr8POWJrpnndM08p2vmGV0vz+maeU7XLG/x8fF5/kzDcCIiIiL5UFgSERERyYfCkoiIiEg+FJZERERE8qGwJCIiIpIPhSURERGRfCgsiYiIiORDYUlEREQkHwpLIiIiIvlQWBIRERHJh8KSiIiISD4UlkRERETyobCUhwMHDjBgwABuueUWOnTowJIlS+wuSURERGygsJSHIUOG0K5dO7Zu3cr48eOJjo726PlTpkxh06ZNJVRd4VmWxfTp0/n3v/9tdyk5nDlzhgEDBnD27Fm7SxEREclBYSkXSUlJ/Prrrzz11FP4+fkBULNmzQI//5tvvmHPnj106NChpEoslL1799K3b182b95sdymXqFq1Kt26dWP69Ol2lyIiIpJDJbsLmDsX5swp/PNTU4MIDMz/MeHhEBZW8GNu376d0NBQMjMzSUhIYMKECTz99NMFfv6MGTPo06dP9vcxMTEsXryYqlWr0qZNGz799FPWr19PXFwcUVFRBAYGkpqayuLFi1m2bBnz5s3D19eX2rVrM3r0aBo1akRcXBzjxo1j1apVANnfjx49mkmTJlGvXj0OHDhAQEAAEydOpHHjxpfU5XK56N69O4H5XLC4uDhee+016taty+7du6lSpQqDBw9m3rx57Nu3j06dOvHiiy8CsH79eqKjozl//jwBAQEMHz6cVq1akZmZySuvvMK3337LmTNnsCyL8ePH07p1a+Li4pg2bRoNGzZk9+7dpKWlMWbMGNq2bcv999/PlClT6N+/P7Vr1y7w9RYRESlJ6lnKRVJSEjfccANhYWH06NGDKlWq0LFjRwAmT55M7969GTZsGOfPn7/kuSdPniQ+Pp527doBsGnTJpYuXcrixYtZunQpZ86cyfH43bt3M3XqVFasWEF8fDyzZ89m7ty5rFixgm7duvHcc89hWVa+9SYkJBAeHs7KlSvp3r07w4YNy/VxY8aM4eGHH77s69+5cyfPPvssH3/8MbVq1SImJoZZs2axdOlSFixYwJEjR0hJSWHatGnExMSwfPlyxo0bx+DBg0lNTeXbb7/l6NGjLFq0iNWrV/PII4/w9ttvZx9/x44dhIeHs3z5cnr27Mmbb74JwBVXXEFoaCgbNmy4bI0iIiKlxfaepbAwz3p9LpaYuJ+QkJDiKwhITEzkoYceIiwsjIMHDzJmzBgmT57ME088wZEjR1iwYAHR0dF88skndOvWLcdz3W43derUoXLlygBs2LCBzp07c9VVVwHgdDrZunVr9uPr169PgwYNABOsunTpkj3k1717d6Kiojh48GC+9TZv3pw2bdoA0KNHD8aOHcsvv/xCjRo1CvX6r7vuOlq0aAFAUFAQ1apVo3LlytSsWZOqVaty4sQJvvrqK44ePUq/fv2yn+fj48P+/ftp1aoV1atXZ+HChRw4cIC4uDiqVq2a/bhrr702+/9ZixYtWLZsWfbPgoKC2LdvX6HqFvFGLpeLyMhI9u/fT1BQEFFRUTidTrvLEpELqGfpIhkZGezdu5cWLVrg6+tLUFAQoaGhgBmea9++PQAdOnRg+/btlzzf19eXjIyM7O8rVaqUo2coaw5UlguHxHLrQbIsi/T0dHx8fHL8/MJerYuPaVnWJW2eyAp6WSpVujRTZ2Zmctttt/Hhhx9mf73//vs0adKEzz77LHvY8t577+WJJ57I8dyAgIDsP1/8ujIyMopUu4g3cblcRERE4Ha7sSwLt9tNREQELpfL7tJE5AIKSxfZt28fZ8+eZePGjWRkZJCYmMjixYt55JFHOHnyJFdeeSUA1apV48SJE5c8v2HDhhw/fpxz584BcOedd7JmzRpOnToFwOLFi/M8d/v27Vm9ejXHjx8HYMmSJVx99dU4HA5q1qzJ4cOHOXbsGJZlsW7duuznJSUlkZSUBMCiRYsIDQ3N7skqKW3btmXz5s3s3bsXMD1oDz74IOfOnWPz5s3cfffd9O7dmxtvvJF169blCJD5OXjwII0aNSrJ0kXKjMjISFJTU3O0paamEhkZaVNFIpIb24fhypqEhAQaN27Mq6++yogRIwgKCmLUqFHcfPPN7Nq1i9OnTwNw6tQpqlevfsnzr7rqKlq3bs3WrVu58847ue2223jsscfo1asXAQEBNGnShCpVquR67nbt2tGvXz/69u1LZmYmNWvWZNasWfj6+nL99dfz+OOP06NHD+rUqcNdd92V/bzatWvz+uuvc+jQIWrWrMmkSZNK5NpcqEmTJowdO5YXXngBy7KoVKkS0dHRBAYG8vjjjzN06FAeeOAB/Pz8aNOmDWvWrCEzMzPfY6alpfH1118TFRVV4vWLlAX79+/3qF1EbGJZVol8bdu2zSoNCQkJxXq8iRMnWm+99Vae5xo2bJhlWZYVHR1trVy5MtfHxcfHWwMHDrQsy7J27NhhxcbGZv9szpw51pAhQ4qt3q1bt1pdu3b16DnFfc2Ky5IlS6yJEyfaXUauyuo1K8t0zS7P4XBYwCVfDofD7tK8gt5jntM1y9vvuSXXTKNhuIskJibmets9QEhICLVq1aJ3797s3r2bTp065fq40NBQGjVqxMaNG2nUqBHbtm2jW7duPPDAA2zZsoWRI0eW5EvwSqdPn2bVqlUMHjzY7lJESk3W0iEXCgwMVO+qSBmjYbiLJCUl5TtnZvjw4QU6zoWB6I033ihyXXm59dZbs9de8mZXXnklc4qy4JaIF8q66013w4mUbQpLF7nwtn4RkZLmdDpxOp0kJiYW+zIoIlI8NAwnIiIikg+FJREREZF8KCyJiIiI5ENhSURERCQfCksiIiIi+VBYEhEREcmHwpKIiIhIPhSWRERERPKhsCQiIiKSD4WlPBw4cIABAwZwyy230KFDB5YsWVImjiXiCb33RESKTmEpD0OGDKFdu3Zs3bqV8ePHEx0dXSaOVVKmTJnCpk2b8vx5XFwc3bp1K8WK/t/zzz/P999/n6Pt448/5sknnwRg586dPP/889k/e+mll7jnnnuYNm1ajj97q/DwcI4fP16o53ry3jtz5gwDBgzg7NmzhS1VRKRc0t5wuUhKSuLXX3/lqaeeym6rWbOm7ccqKd988w179uxh6NChdpdyibS0NNxuN02bNs3zMTfeeGOOzYoXLVrEZ599xjXXXEPz5s2z/+ytNm/eXKjnefreq1q1Kt26dWP69OkF3jBaRKQisD8szZ0LRdhtPig1FQID839QeDiEhRX4mNu3byc0NJTMzEwSEhKYMGECTz/9dKHqK85jlZQZM2bQp08fwPQujBw5Erfbja+vLy1btmTs2LGXPGfRokXMmzcPX19fateuzejRo2nUqBFxcXFMmjSJevXqceDAAQICApg4cSKNGzdm/fr1REdHc/78eQICAhg+fDitWrXK85y+vr588cUX3HbbbQBMnz6dlStXcvXVV+NwOLJriYuLY9y4caxatYrevXtjWRYDBw7M7o0aOHAgL730EidPnsz1/HFxcURFRREYGEhqaiqLFy/m888/z/HYXr16ERISQlxcHNOmTaNhw4bs3r2btLQ0xowZQ9u2bVm8eDHvvPMOvr6+1KhRg1dffZX69esX6nVnGTlyJAB9+/YlJiaG+vXr53ntL1aY997999/PlClT6N+/P7Vr1/bwnSQiUj5pGC4XSUlJ3HDDDYSFhdGjRw+qVKlCx44dOXXqFD179qRVq1aXDAt5eiyAHTt20KtXL5xOJy+88ALnz58vyZeVq5MnTxIfH0+7du0AWLt2LWfOnOHDDz9k8eLFgJn3cqEtW7Ywe/Zs5s6dy4oVK+jWrRvPPfcclmUBkJCQQHh4OCtXrqR79+4MGzaMlJQUpk2bRkxMDMuXL2fcuHEMHjyY1NTUfM/56aefct9997Fu3TrWrFnD8uXLWbhwIadPn8719SxYsACA2NhYkpOTs/9cu3btPM8PsHv3bqZOncqKFSs4fPjwJY999dVXsx+7Y8cOwsPDWb58OT179uTNN98kKSmJKVOmMHv2bFauXMk999xDdHR0oV93lgkTJmS/hvr161/22l8ov/fe5MmT6d27N8OGDcvxvrviiisIDQ1lw4YNl3nniIhUHEXqWWrWrFldIB7omJycnFSog4SFedTrc7H9iYmEhIQU+vm5SUxM5KGHHiIsLIyDBw8yZswYJk+ezKhRo4iJiWHSpElFPtbLL7/MNddcQ2xsLAEBAUydOpVPP/2Uzp07F+truRy3202dOnWoXLkyAK1bt2batGk8+eST3H777fTt2xeHw8H//ve/7Ods2rSJLl26ZA/pdO/enaioKA4ePAhA8+bNadOmDQA9evRg7NixfPTRRxw9epR+/fplH8fHx4f9+/fnec7MzEy++eYb/vnPf/LKK6/QsWNHrrzyyuzjzps3r8Cvc/PmzXmeH6B+/fo0aNCgQI+99tprs99zLVq0YNmyZWzZsoX27dtTv359gOznulwuj193fvK79g0bNszx2Lzee0888QRHjhxhwYIFREdH88knn+SYjxYUFMS+ffsKeGVFRMq/QoelZs2a+QOzgN+Krxz7ZWRksHfvXlq0aIGvry9BQUGEhoZy7Ngx/P39PZpvlN+xAOrWrZv9WH9//xzDL6XF19eXjIyM7O8bNmzI2rVriYuLY+vWrTz11FOMGjWKGjVqZD8mt14My7JIT08HwM/P75KfZWZmctttt/H6669nt//444/UrVsXPz+/XM9Zr149brjhBvz8/PDx8clx3ovPcTn5nX/btm0EXjCUm9tjN27cSJMmTdi2bRsBAQHZ7Vl1ZdWY5ezZsxw6dKhQrzu/wHy5a58lv/fe9u3bad++PQAdOnRg6dKlOcJSRkZGdngWEZGiDcNNAd4CDhdTLWXCvn37OHv2LBs3biQjI4PExEQWL17MI488ku/zRowYwYgRIwp1rEOHDrF582buvvvuPI9VUho2bMjx48c5d+4cYIaxRo4cSfv27Rk2bBjt27dn9+7dOZ7Tvn17Vq9enX2H1pIlS3LMI0pKSiIpyXQ0Llq0iNDQUDp16sTmzZvZu3cvABs2bODBBx/k3LlzeZ5z3bp13HvvvYD5pf7xxx9z8uRJMjMz+fDDDz16nW3bts3z/AV57F//+tdcH5vl1ltvZcuWLRw9ehSAhQsXMnny5HzPW5BrDSYYZoWhy137LPm9906ePJndQ1etWjVOnDiR47kHDx7MdQ6UiEhFVaiepWbNmvUDfkpOTv6kWbNmI/N6XGJiYmHrKrCzZ88W63k2bNjAddddx7hx4/jHP/7BNddcQ79+/bjiiiuyz/Prr7/yww8/5OiR2bNnDx06dMhRS0GOlZqayvjx43nuuefYs2dPnsfKMnbsWEJCQnj00UcL3X7xNWvevDkffPABrVu3JiQkhHXr1nHfffdxxRVXUKdOHXr06MG+ffs4d+4ciYmJ1KxZk86dO9OrVy8sy+Kqq65i+PDhJCcn43a7qV69OuPGjePo0aNUr16dv/zlL6Snp/P0008zaNAgwPRoDR8+HLfbnec5x4wZQ8eOHUlMTKRu3brccccdPPDAA1x55ZUEBwdz5swZEhMTcbvd2bVl2b17N1dddVWOP+d1/tyef/Fjhw4dmutjs77PzMzE6XRmL2dQo0YNBg8eXKjXffH/91tvvZWePXvy4osv4nA48rz2F8rvvXfmzBm+//57GjRowJ49e8jMzMw+5/nz5/nqq68ICwsr8t+r4v67Wd7penlO18xzumaFZFmWx19Nmzbd2LRp0w1Nmzb9rGnTpr82bdr0y6ZNm15z4WO2bdtmlYaEhIRiPd7EiROtt956K9/HDB8+3EpOTs7+/ty5c1bnzp2ttLQ0j451/vx5a8CAAdYXX3xx2WMVp4uvWXx8vDVw4MBiOfbWrVutrl27FsuxypLifp+VtPzeewkJCdawYcMsy7Ks6Ohoa+XKldk/W7JkiTVx4sRiqcHbrpnddL08p2vmOV2zvP2eW3LNPYUahktOTr4jOTn5zuTk5LuAb4Cw5OTk/13maV4hMTGRxo0b5/nzgQMH8vnnnzN69GiWLl0KQOXKlfnoo4/w9/f36FirVq1ix44dzJw5kyeffJLVq1fneaySFBoaSqNGjdi4cWOpnVNKVn7vvZCQEGrVqkXv3r3ZvXs3nTp1AuD06dOsWrWKwYMHl2apIiJlnv3rLJUxSUlJ+c7XePvtt4vtWA8//DAPP/ywR/WVlKz1fIrq1ltvZdWqVcVyLCm8y733clt08sorr2ROEdY8ExEpr4ocln7vXSo3tm7dWiaPJeIJvfdERIqPFqUUERERyYfCkoiIiEg+FJZERGzkcrkIDg6mZcuWBAcH43K57C5JRC6iCd4iIjZxuVxERERk7zvodruJiIgAwOl02lmaiFxAPUsiIjaJjIzMDkpZUlNTiYyMtKkiEcmNwpKIiE2yNmcuaLuI2ENhSUTEJkFBQR61i4g9FJZERGwSFRVFYGBgjrbAwECioqJsqkhEcqOwJCJiE6fTSUxMDA6HAx8fHxwOBzExMZrcLVLGKCzl4cCBAwwYMIBbbrmFDh06sGTJEq85nzfXLiIiUtYoLOVhyJAhtGvXjq1btzJ+/Hiio6M9PsaMGTOYMWNGqZ2vuI9V0PqLs/aimjJlCps2bbLt/Ln58MMPefDBB3nooYd4/PHH2blzZ4GeN3/+fLp27Uq3bt149tlnOXbsWAlXao8zZ84wYMAAzp49a3cppS5r6QC3241lWdlLB2itJZGyRWEpF0lJSfz666889dRT+Pn5AVCzZk2vOJ83115U33zzDXv27KFDhw62nD83P/zwA5MnT2b27Nl8+OGHPPvsswwePPiyz/vuu++YM2cOCxcuZNWqVQQHBzN9+vRSqLj0Va1alW7dupXb15cfLR0g4h1sX5Ry7txvmTPn60I/PzU1lcDAuHwfEx7eirCwmwp8zO3btxMaGkpmZiYJCQlMmDCBp59+utA1lub5vLn2opoxYwZ9+vTJ/j4mJobFixdTtWpV2rRpw6effsr69euJi4vLnlibmprK4sWLWbZsGfPmzcPX15fatWszevRoGjVqRFxcHOPGjWPy5MkA2d+PHj2aSZMmUa9ePQ4cOEBAQAATJ06kcePGOWqqXLky48ePp27dugDccMMN/Pzzz6SlpVG5cuU8X8sNN9zAJ598gr+/P+fOnePIkSNcd911xXKd4uLieO2116hbty67d++mSpUqDB48mHnz5rFv3z46derEiy++CMD69euJjo7m/PnzBAQEMHz4cFq1akVmZiavvPIK3377LWfOnMGyLMaPH096ejrTpk2jYcOG7Ny5E19fX8aMGUPbtm05c+YMI0eOxO124+vrS8uWLRk7diy+vr7cf//9TJkyhf79+1O7du1ieZ3eQEsHiHgH9SzlIikpiRtuuIGwsDB69OhBlSpV6NixIwCTJ0+md+/eDBs2jPPnz5fo+U6dOkXPnj1p1aoV33//fZFrz6p/5MiRxVZ/fufbsWMHvXr1wul08sILLxTb9crNyZMniY+Pp127dgBs2rSJpUuXsnjxYpYuXcqZM2dyPH737t1MnTqVFStWEB8fz+zZs5k7dy4rVqygW7duPPfcc1iWle85ExISCA8PZ+XKlXTv3p1hw4Zd8pjrrruOu+66CwDLspgwYQL33HNPvkEpi7+/P+vWreOOO+7gq6++onv37gW8Gpe3c+dOnn32WT7++GNq1apFTEwMs2bNYunSpSxYsIAjR46QkpLCtGnTiImJYfny5YwbN47BgweTmprKt99+y9GjR1m0aBGrV6/mkUce4e233wbM//fw8HCmTZtGz549efPNNwFYu3YtZ86c4cMPP2Tx4sWAme8GcMUVVxAaGsqGDRuK7TV6Ay0dIOIdbO9ZCgu7yaNen4slJiYSEhJSjBWZYz700EOEhYVx8OBBxowZw+TJk3niiSc4cuQICxYsIDo6mk8++YRu3brleO7TTz9NfHw8AOfOnQMgNjYWgNatWzNr1qwCn2/UqFHExMQwadKkItf+8ssvk5SUxJEjR5gwYQKfffZZsdSf3/muueYaYmNjCQgIYOrUqXz66ad07ty5wK/FE263mzp16mSHkA0bNtC5c2euuuoqwNx1tHXr1uzH169fnwYNGgAmWHXp0iV7+LB79+5ERUVx8ODBfM/ZvHlz2rRpA0CPHj0YO3Ysv/zyCzVq1LjksampqYwYMYL//e9/zJ49u8Cv67777uO+++7j/fffp3///qxduxZf36J/xrnuuuto0aIFYH4xV6tWjcqVK1OzZk2qVq3KiRMn+Oqrrzh69Cj9+vXLfp6Pjw/79++nVatWVK9enYULF3LgwAHi4uKoWrUqANdeey0hISEkJibSokULli1bBpj3z7Rp03jyySe5/fbb6du3Lw6HI/vYQUFB7Nu3r8ivzZtERUXx1FNP5fgg4e/vr6UDRMoY28NSWZORkcHevXtp0aIFvr6+BAUFERoayrFjx9i+fTvt27cHoEOHDixduvSSsHFhmMiaHJ3fHJX8zufv7+/R/J/8jgUUe/2XO1/W0BOYXwDF8Us+L76+vmRkZGR/X6lSpRw9Q1nzqbJcuLZNbj1IlmWRnp6Oj49Pjp9f+Evt4mNalnVJG8Dhw4d55plnaNy4MXPnziUgIOCyr8ftdvPTTz/lCGMvvfQSJ06cyDWMeerinq1KlS79pyAzM5PbbruN119/Pbvtxx9/pG7dunz22WfZv+jvvfde/vCHP7BixQqAHK/vwuvXsGFD1q5dS1xcHFu3buWpp55i1KhR2QE6IyOjQD1u5Y2Pj0++34uI/TQMd5F9+/Zx9uxZNm7cSEZGBomJiSxevJhHHnmEkydPcuWVVwJQrVo1Tpw4UaLny8+IESMYMWKER8cq7voLWvuhQ4fYvHkzd999d5HOl5+GDRty/Pjx7N6wO++8kzVr1nDq1CmA7GGf3LRv357Vq1dz/PhxAJYsWcLVV1+Nw+GgZs2aHD58mF9//RXLsli3bl3285KSkkhKSgJg0aJFhIaGZvdkZfn111/p06cPnTp1Ytq0aQUKSgA//fQTL7zwQnZNK1eupEmTJsUSlAqqbdu2bN68mb179wKmt+7BBx/k3Llz2f8/e/fuzY033si6detyhNXcLFiwgJEjR9K+fXuGDRtG+/bt2b17d/bPDx48SKNGjUr0NZU1kZGRpKWl5WhLS0vTBG+RMkY9SxdJSEigcePGvPrqq4wYMYKgoCBGjRrFzTffzK5duzh9+jQAp06donr16iV6vvz8+OOPdO3a1aNjVatWrVjrL0jtp0+f5h//+AcTJkzA398/u33AgAG0adOGZ555JscxPW3PctVVV9G6dWu2bt3KnXfeyW233cZjjz1Gr169CAgIoEmTJlSpUiXX57Zr145+/frRt29fMjMzqVmzJrNmzcLX15frr7+exx9/nKFDh3Lttddmzz8CqF27Nq+//jqHDh2iZs2auQ6Xvvfee/z444+sXbuWtWvXZre/++671KhRg4EDB/L4449z77335nhe1msNCwvDz8+PunXr8q9//SvX+i88Rl5/LowmTZowduxYXnjhBSzLolKlSkRHRxMYGJh9TR544AH8/Pxo06YNa9asITMzM8/jPfzww3z55Zd06dKFKlWqcO211xIWFgaYgPD1119XuOEnTfAW8RKWZZXI17Zt26zSkJCQUKzHmzhxovXWW2/lea5hw4ZZlmVZ0dHR1sqVK0v0fFmGDx9uJScnZ39/7tw5q3PnzlZaWppHx8qqPyEhoVjqv9z5zp8/bw0YMMD64osvinSegoqPj7cGDhxoWZZl7dixw4qNjc3+2Zw5c6whQ4YU+tgXv8+2bt1qde3atdDHy7Jo0SJrzZo1RT5OWeTJ380lS5ZYEydOLMFqyiaHw2EBl3w5HA67S/MKxf3vf0Wga5a333NLrplGw3AXSUxMvOT27ywhISHUqlWL3r17s3v3bjp16lSi5wPTa/D5558zevRoli5dCpj5Jh999FGOnpqCHCur/pEjRxZL/Zc736pVq9ixYwczZ87kySefZPXq1UU63+WEhobSqFEjNm7cSKNGjdi2bRvdunXjgQceYMuWLYwcObJEz18Yfn5+OXqrKqLTp0+zatWqAq0/Vd5obzgRL5FXiirql7f2LN16663Wnj17ivWYpXW+gh6ruK5ZaV8rO+nTmOd0zQpm/vz5lsPhsHx8fCyHw2HNnz/f7pK8ht5jntM1y1t+PUuas3SRC28v97bzeXPtIhWV0+nE6XSWyDIoIlI8NAwnIiIikg+FJREREZF8KCyJiIiI5ENhSURERCQfCksiIjZyuVwEBwfTsmVLgoODcblcdpckIhfR3XAiIjZxuVxERESQmpoKmD0BIyIiAHOXnIiUDepZEhGxSWRkZHZQypKamqq94UTKGIUlERGbaG84Ee+gsCQiYpOgoCCP2kXEHgpLIiI20d5wUloGDRpEpUqVaNGiBZUqVWLQoEF2l+RVFJZERGzidDqJiYnB4XDg4+ODw+EgJiZGk7ulWA0aNIjo6GgyMjIAyMjIIDo6WoHJAwpLIiI2cjqdpKSksGvXLlJSUhSUpNjFxMR41C6XUlgSEbGR1lmSkpbVo1TQdrmU1lkSEbGJ1lmS0uDn55drMPLz87OhGu+kniUREZtonSUpDVkBvKDtcin1LImI2ETrLL6ojWQAACAASURBVElpmDlzJmDmKGVkZODn50dERER2u1yeepZERGyidZaktMycOZP09HQSEhJIT09XUPKQwpKIiE2ioqLw9c35z7Cvr6/WWRIpYxSWRERssnnzZjIzM3O0ZWZmsnnzZpsqkvJKd10WjcKSiIhNZs2a5VG7SGFk3XXpdruxLCv7rksFpoJTWBIRscnFvUqXaxcpDN11WXQKSyIiIuWY7rosOoUlERGbVK1a1aN2kcKoWbOmR+1yKYUlERGbzJo165JVlP38/DRnSaSMUVgSEbGJ0+kkNjYWh8OBj48PDoeD2NhYbXUixer48eMetculFJZERGzkdDpJSUlh165dpKSkKChJsdPip0WnsCQiIlKORUVFERgYmKMtMDBQi596QGFJRESkHHM6ncTExOQY7o2JiVEvpge0ka6IlBtz537LnDlf212Gx6qdP03dX1LYW/cGu0vxWHh4K8LCbrK7DLkMp9OJ0+kkMTGRkJAQu8vxOgpLIiI2anFyLy8lzOSKjHN0r/MGmT7q8BcpaxSWRKTcCAu7yXt6OSwL3niDzL9P4KCPD13T0/nZPZaoqCgNj4iUMQpLIiKl7cQJ6N8flizhP35+hKWn8yvA73t2AQpMImWI+ntFRErTN99AmzawfDlRV1/NgxkZJij9Tnt2iZQ9CksiIqXBsmD2bGjbFn77DTZsYPSJE7k+VHt2iZQtCksiIiXtzBno1w8GDoQ77oCvv4Z27bRnl4iXUFgSESlJSUlw660wbx7885/w0UdQp47dVYmIBzTBW0SkpLz3nulNCgyETz6Bjh1z/Fh7dol4B/UsiYgUt3PnYNAg6N0bbr7ZDLtdFJRAe3aJeAuFJRGR4rRvH7RrB9HRMHQo/Pe/0KBBrg/t0qWLR+0iYg+FJRGR4rJiBYSGwt69sHw5TJ4M/v55Pvz999/3qF1E7KGwJCJSVOfPw7Bh8NBD0LgxbN9u/nwZx44d86hdROyhCd4iIkVx6BA8/jh8/rmZpzR1KgQE2F2ViBQjhSURkcJau9ZM4v7tN1iwAJ54wu6KRKQEaBhORMRTGRnw8svw5z9DvXqwbVuhglKtWrU8ahcReygsiYh44uhRuP9+s8Dkk09CXBw0b16oQ02fPh3/iyaA+/v7M3369GIoVOT/uVwugoODadmyJcHBwbhcLrtL8ioKSyIiBfX559CqFWzcCG+/De++C1WrFvpwTqeTd955B4fDgY+PDw6Hg3feeQen01l8NUuF53K5iIiIwO12Y1kWbrebiIgIBSYPKCyJiFyOZcGUKXDXXWY17q1bYcAA8PEp8qGdTicpKSns2rWLlJQUBSUpdpGRkaSmpuZoS01NJTIy0qaKvI/CkohIfn75BR5+2CwN8PDDZn7SzTcX2+E1PCIlze12e9Qul9LdcCIieYmPh0cfhQMH4PXX4fnni6U3KUvW8EjWp/6s4RFAPUxSbPz8/MjIyMi1XQrGa3uW9GlMREqMZZntSm6/HdLTYdMmGDKkWIMSaHhESkduQSm/drmUV4YlTVYTkRJz+jQ4nWaByXvvNZvgtm1bIqfav3+/R+0iheFwODxql0t5ZVjSpzERKRHffQe33AKLFkFUFKxaBSW45lFQUJBH7SKFERUVRWBgYI62wMBAoqKibKrI+3hlWNKnMREpdnPnwp/+ZCZ0r1sHL74IviX7T6R+iUlpcDqdxMTE5FiiIiYmRvPiPOCVYUmfxkSk2Pz2GwwcCH37mrD09ddw992lcmr9EhPxDl4ZlvRpTESKxZ49ZhL37NmmJ2ndOqhfv1RL0DpLUtI0z7fovDIs6dOYiBTZkiUQGgr798N//mPmKFXSaipS/mieb9F5ZVgCfRoTkUJKS4O//hV69oSQENi+Hbp0sbsqkRKjeb5F57VhSUTEY/v3w513wvTpZoHJTZtAt09LOVezZk2P2uVS6nMWkYrho4+gTx84fx7ef9+szC0iUgCFCkvNmjXzB+YAwcAVwPjk5OQVxViXiEjxSE+Hf/7TzEn64x/hgw+gaVO7qxIpNcePH/eoXS5V2GG4PsCx5OTkDkBn4M3iK0lEpJj873/QsaMJSv37w9atZS4oaesmKWlabqfoCjsM9wGw+Pc/+wDpuT0oMTGxkIcvuLNnz5bKecoTXTPP6Zp5zu5rFvjllzQYOhTf06f53yuvcOLhhyElxbZ6crNq1SrGjBnD2bNnAbOR7oABAzh8+DDdunWzubqyz+73mLd47rnncrzPAAICAnjuued0/QrKsqxCfzVt2rRa06ZN/9u0adPeF/9s27ZtVkmaP3++5XA4LB8fH8vhcFjz588v0fOVJwkJCXaX4HV0zTxn2zXLyLCsV16xLF9fy2rWzLJ27LCnjgJwOBwWcMmXw+GwuzSvoL+XBaffmZf3e27JNe8UeoJ3s2bNGgLLgJnJyckLiiW5FVDWAltZ60ZkLbAFaAkBkYrs2DEIC4PVq6FXL3j7bahWze6q8uR2uz1qFyksp9OJ0+kkMTGRkJAQu8vxOoWas9SsWbN6wBpgeHJy8pziLenytMCWiFwiLs4sMrluHfzrX/Dee2U6KAH45rH3XF7tImKPwv6NfBGoAYxu1qzZZ79/VSnGuvKlBbZEJJtlwRtvQIcOZuPbzZth0CDw8bG7ssvKzMz0qF1E7FGoYbjk5OQhwJBirqXAgoKCcu2m1sx+kQrm5EkYMMAsB/DAAxAbCzVq2F2ViJQzXtnXq410RYRvv4XWrWHpUpg0CZYv97qg5JNH71de7SJiD68MS9pIV6SCmzMH2raFM2fgv/+FYcPMEJyXsSzLo3YRsYf3/evyO22kK1IBpabCU0+ZBSbbtYNvvjFzlbyUI4996fJqFxF7eG1YEpEKJjkZbr3VzEsaMwY++QTq1rW7qiLRlAIR76CwJCJl36JF0KaN2b7k44/h5ZfBz8/uqopMUwpEvEOhF6UUESlx587BCy/AzJlw++0mNF13nd1VFSstFihS9qlnSUTKppQUaN/eBKW//x0++6zcBSUR8Q7qWRKRsmflSrNtiWXBsmXw8MN2VyQiFZjCkoiUHenpEBlp1k1q1QoWL4Y//MHuqkRymDv3W+bM+druMjxS5+wxOh3ZwrorQzhSq7Hd5XgsPLwVYWE32XZ+hSURKRsOHYInnoBNm+Dpp+H11yEgwO6qRLxa49P76XXgY+4+GocPFqu4m61X/ESjRo2oV6+e3eV5DYUlEbHfunXQu7dZR2n+fNDdYFKGhYXdZGsvx2VZlvk7NXkyxK/lfEAAb1XyZXJ6OgdYD+fWs39/IKNG6c7LgtIEbxGxT2YmjB0LnTpBnTrw1VcVLii5XC6Cg4Np2bIlwcHBuFwuu0sSb3X+PLhcEBpq/k7t3AkTJtC6dm2eT0/nwAUPTU1NJTIy0rZSvY3CkojY46ef4P774aWXTED68kuoYLfOu1wuIiIicLvdWJaF2+0mIiJCgUk8c+oUvPYaNG4MffpAWprZEiglBUaM4LtDh3J92v79+0u3Ti+msCQipW/zZjOBe8MGmDUL5s6FqlXtrqrURUZGkpqamqNNn/ilwA4fhhEjoGFDs7zGH/4Aq1aZHqWnnoIrrgCgZs2auT49r3a5lOYsiUjpsSzzCXjECHA4YMsWE5oqqLw+2esTv+Rr1y6YMsUMuWVkQM+eMHQo3HKL3ZWVW+pZEpHS8euv0L27+Uf9gQcgPr5CByWAoKAgj9qlArMsszBr165www3w/vvmrtHdu83K9vkEpePHj3vULpdSWBKRkrd9O7RubYYIXnsNliyB6tXtrsp22khXLis93YShP/0J7r4btm2DceNg/36YMaNA65AplBedwpKIlBzLMnOSbr/dTDrdsAH+9jfw8bG7sjJBG+lKns6cgTfegCZN4PHH4eRJ83cpJQVGjYJatQp8KIXyolNYEpES4XPmjLkz55ln4K674OuvTWiSHJxOJykpKezatYuUlBQFpYruyBEThho2hCFDoEEDs+VPYiJERECVKh4f0ul00rdvX/z8/ADw8/Ojb9++eq95QGFJRIpfQgKNevWChQvNkMHq1VC7tt1ViZRdyckmDDkc8Mor5gPGF1/A55+bvRF9C//r2uVyERsbS0ZGBgAZGRnExsZqiQoPKCyJSPGaPx9uuQW/Eydg7VrzKbkI/9CLlFuWZcLQQw9B8+Ywb5655T8pCZYuhdtuK5bTaImKotO/YCJSPM6eNXfnPPkktGnDviVL4J577K5KpOzJyDBh6PbboUMHs+7YmDHgdkN0NDRtWqync7vdHrXLpRSWRKTo9u41n4JjYmD4cPj0U9Lr1rW7KpGy5bffTBhq3hx69ICjR+HNN82dbS+/DCX0d8Ynjxsq8mqXS2lRShEpmqVLzdCBnx+sXAndutldkUjZ8vPP8K9/mWD0889mGYAPPoBHHjF/b0qYZVketculFJZEpHDS0sxK3NOmmQXx3n8fgoPtrkqk7Nizx6wr9u67plfpgQfMoqwdOmj5DC+jsCQinjtwAHr1MtuVDB4Mkydn70MlUuHFxZm/E0uXgr+/mcf397/btlF0rVq1OHbsWK7tUjCasyQinvn4Y7NNyc6dZmXhN95QUBLJzIQVK+COO6BtW/j0U9PzmpICs2fbFpQApk+fTuXKlXO0Va5cmenTp9tUkffx2rDkcrkIDg6mZcuWBAcHa70IkZKWkQGjR0OXLnDttWZvt8ces7sqEXudPWvCUMuWZgmA/fvh9ddN7+srr0D9+nZXiNPppH///jkWpezfv78WpfSAVw7DuVwuIiIisteNcLvdREREAOh/vkhJOHIEeveG9evNZO4334SLtk8QqVCOHzd3ts2YYf5+tGoFCxbAo49CpbL1qzWvRSnbtWun35kF5JU9S1pgS6QUbdxofhF88QXMmWO+FJSkokpJMduQBAWZBVdbtYJ160xP6xNPlLmgBPqdWRzK3v/VAtACWyKlIDMTJk2CyEho3NjMVfrjH+2uSsQe8fEwZYq55d/Hx/S0Dh0KN95od2WXtX//fo/a5VJe2bPkl8e6FHm1i4iHjh838y9GjjSL523bpqAkFY9lwUcfmZXo27Qxexy+8ALs2wexsV4RlACCgoI8apdLeWVYyhp3LWi7iHjgyy8hNBQ++cTMx1i0CK66yu6qREpPWppZG+mPfzQ3NHz/vVkKYP9+09t63XV2V+iRqKgoAi8aOg8MDCQqKsqmiryPV4alvNaG0JoRIkVgWWbidvv25vvPP4e//EWL50nFceKECUONGpkbGXx8YO5c+OEHM+RWvbrdFRaK0+kkJiYGh8OBj48PDoeDmJgYTe72gFfOWRKRYnbyJAwcaFbh7trV/IKoWdPuqkRKx4ED5nb/t9+GU6fg3nvNjQydOpWbDwtOpxOn00liYiIhNq755K28MiwdP37co3YRyceOHdCzp9kMd8IE+Mc/wNcrO51FPPPtt2bS9sKFpmf1scdMD1JoqN2VSRnjlf8iarKaSDF55x249VY4fdqsoTRihIKSlG+WZW71//Of4eabYdkyM9y8d69ZJ6mcBiUt5Fw0XvmvoiariRRRaiqEh5uv22+Hr7+GO++0uyqRknP+PLhcJgx17Gh6VCdMMENw06aBw2F3hSXG5XIRHh6O2+3Gsizcbjfh4eEKTB7wyrCkyWoiRfD992bvqnffNduXrFkD9erZXZVIyTh1Cl57zawV1qcPnDsH//63WVxyxAioUcPuCkvckCFDSEtLy9GWlpbGkCFDbKrI+3jlnCXQZDWRQnn/fejf32x8u3o1dO5sd0UiJePwYbPJ81tvmbvc7rwTZs40SwFUsKHmY8eOedQul/LasCQiHjh3zkxcffNNuO02s3ZSw4Z2VyVS/HbtgqlTYf58s/lzjx7mvf+nP9ldmXgxhSWR8s7tNpt7fvUV/O1vMHEiVK5sd1UixceyYMMGs3Dk6tVQpQpERJj3e+PGdldnu1q1auXai6S1CQuuYvVFilQ0q1aZjT6Tk2HJEjN3Q0FJyov0dDO0/Kc/wd13mw8EY8eaSdtvvqmg9Lvp06fj7++fo83f35/p06fbVJH3UVgSKY/S083k1QcegOBg2L4dune3uyqR4nHmDMyYQeMuXaBXLzMn6a23TC/q6NGgHpMcnE4nAwYMyN4/1c/PjwEDBuimKA8oLImUN4cPmxWIX33VDEV88YU+YUv5cOSICUNBQfD886TXqWPWSUpKgqefNsNvcgmXy0VsbGz2/qkZGRnExsZq6QAPKCyJlCfr15tht23bYN48mDULAgLsrkqkaJKTTfB3OCAqCu64AzZvxu1ywcMPV7i72zwVGRlJampqjrbU1FQiIyNtqsj76B0mUh5kZsL48WaxvVq14MsvzZoyIt5s82YThkJCzH6F/fqZXqRly8xiqlIg+/fv96hdLqW74US83c8/m2D0ySfQu7fpTbrySrurEimcjAz48EOzZ9uWLWZD51GjzJYkdevaXZ1XCgoKwu1259ouBaOeJRFvtmWLGXb773/NBNf58xWUxDv99pt5DzdvbtZGOnLE3NG2f7+5w01BqdC0RVjRKSyJeCPLMvtZ3XEH+Pub0PT00+DjY3dlIp75+Wd4+WUzafvZZ832I++/b7blee45qFrV7gq9nrYIKzoNw4l4m19/NRvgLltm5nO88w5cfbXdVYl4Zu9es+7XO++YXqVu3WDYMOjQQaG/BGiLsKJRWBLxJl9/DT17mqGJqVPNCsX6xSLeJC7OrLS9dKnpFe3TB/7+d2jRwu7KRPKksCTiDSwL3n4bnn8eateGzz6Ddu3srkqkYDIz4T//MSFp0yaoXh2GDzfv5/r17a5O5LIUlkTKujNn4JlnzOTtTp3Mf+vUsbsqkcs7e9a8X6dONbf8BwWZuXb9+0O1anZXJ1JgCksiZVlCgtkENzHRTIKNjITftywQKbN++QWio+GNN8xdbTffDC6XeS9ftEeZiDfQ3XAiZZXLBbfcAj/9BGvWwJgxCkpStqWkwJAh0LChCfY33wzr1pm9CXv3VlCy0aBBg6hUqRItWrSgUqVKDBo0yO6SvIp6lkTKmrNn4a9/NYtLtm8PCxdCgwZ2VyWSt+3bzXykDz4wNxz07m0mbf/xj3ZXJpigFB0dnf19RkZG9vczZ860qyyv4rU9Sy6Xi+DgYFq2bElwcLA2BJTy4YcfzMTtWbPgH/8wi00qKElZZFnw0Udm0+bWrc0E7r/9Dfbtg9hYBaUyJCYmxqN2uZRX9iy5XC4iIiKyNwZ0u91EREQAaJEt8V7Ll5u9r3x8zHYPDz5od0Uil0pLg/feM9uRfPedCfOTJpmNbqtXt7s6yUVGRoZH7XIpr+xZ0g7KUq6cP2+GLB55BJo0MUMaCkpS1pw4YUJRo0Ym1IPpQfrhB7OYpIJSmeWXx1zHvNrlUl4ZlrSDspQbBw/CXXeZlYwHDYLPPze/jETKigMHYOhQM2l7+HCzd9tHH8GOHRAWBpUr212hXEbWyEtB2+VSXjkMpx2UpVxYswacTjOh+7334PHH7a5I5P/t2GGG2t57z8xPeuwxE5pCQ+2uTDyUNYk7JiaGjIwM/Pz8iIiI0ORuD3hlz1JUVBSVL/o0U7lyZe2gLN4hIwNeegk6d4ZrroFt2xSUpGywLHOr/5//DDfdZLYkee452LMHFixQUPJiM2fOJD09nYSEBNLT0xWUPOSVPUsAlmXl+71ImXT0qLmt+tNPoW9fmDkTAgPtrkoquvPn4f33TU/SN9+YEP/KK2bl+Bo17K5OxHZe2bMUGRnJ+fPnc7SdP39eE7ylbNu0ySzSt3kz/PvfZrd1BSWx06lTZvuR6683G9qePQuzZ5vFJUeOVFAS+Z1X9ixpgrd4lcxM84n9xRfN5O2PPjJDHCJ2OXzYbEXy1lvmLrc77oA334SuXcHXKz9Di5QorwxLmuAtXuOXX8xw28qV0LOn6VG66iq7q5KKKiHBBPf5883cue7dzaTtW2+1uzKRMs0rP0JERUUReNHwRWBgoCZ4S9ny1VdmQuzHH8P06WZOiIKSlDbLgg0boFs3aNnSbJ8TEQHff2+2J1FQErksrwxLTqeTvn37Zi+o5efnR9++fbV6t5QNlgX/+pfZ1y0jw8xVev55szK3SGlJTzcB/dZbzVpeX34JL78M+/ebIbfGje2uUMRreGVYcrlcxMbGZi/VnpGRQWxsrPaHE/udOmXudvvLX+C+++Drr/XJXUrXmTMwYwY0bQq9esGvv5q5SW43jBkDtWvbXaGI1/HKsKTtTqRM2rkT2rQxn+ZfecXMU6pVy+6qpKI4cgRGj4agINOTec01Zp2kxER4+mmoUsXuCkW8lldO8M5tcnd+7SIlLjYWnn3WzEn69FMz7CFSGpKTzXY5sbFmk9sHHzR7tbVrZ3dlIuWGV4YlPz+/XHdL1qaAUup++w0GDzZ3ud11l9ka4ppr7K5KKoLNm2HyZFixwuzP1rcvvPACNGtmd2Ui5Y5XDsPlFpTyaxcpEbt3Q9u2Jii9+CKsXaugJCUrIwOWLYPbbzc3EGzaBJGRZj7SrFkKSpInl8tFcHAwLVu2JDg4WHN8PeSVPUu+vr5kZmbm2i5SKj74APr3B39/WL0a7r/f7oqkPPvtNzPM9tprJqQ3amQmcT/1FFStand1Usa5XC4iIiKy5/q63W4iIiIAdBd5AXllusgtKOXXLlJs0tLM5NnHHoMWLczdbgpKUlJ+/hnGjgWHw8yJq14dFi0yayT95S8KSlIguimq6LyyZ0nEFm63CUlffglDhsCkSWauiEhx27vX9CK9847pVera1UzavuMOrdclHtMWYUXnlT1LtfK4HTuvdpGiunLDBmjVytyG/cEH8PrrCkpS/L78Eh591KyR9Pbb8Pjj8N13sGoV3HmngpIUSl5bgWmLsIIr0Z6lkrp7uk6dXRw/nkR16yTPsZz5dGS/T33q1GmuO7YLIDU1SJvde+DRA6/x3A9/Z0/Vm3ip5WIOvXk9vGl3VWWf3mcFVzPtf4z6biCc2sxpv+p82OAfLG0wmGM/XAvP2V1d2aX3WMFUqRKHr29yjqkqvr6+VKnSTL8zLzB1at4/88phuHr16gFQee9/iTw/n1HMZ36NgSyvNZ7TNtcm5c91v+3m/boRzG76Oml+WthPil+9s24qZ57jzcav8Z9rBvBbpWp2lyTlSNbvzH379nHu3FmuuCKARo0aZbfL5flYllWoJzZr1swXmAncBJwDBiQnJ+/J+nl8fLzVunXrYikyP7vXr6fJ3Lkwdy7UrGmW83/mGQ2R5CMxMZGQkBC7y/Aqumae0zXzjK6X53TNPKdrlrf4+Hhat26d61h3UeYsPQwEJCcn3waMAPLpwCo56fXrw7vvQnw83HSTmXjbsqVZ5r+QQVBEREQkS1F6ll4DvkxOTl74+/eHkpOTG2T9PD4+3goshcHks2fPEhAQYL6xLKpu3Ei9KVO4Yu9eUkNDOTJsGGdvuqnE6/AmOa6ZFIiumed0zTyj6+U5XTPP6ZrlLTU1Nc+epaLMWboKOHHB9xnNmjWrlJycnJ7VUBpdfZd0KbZoAQMGwJw5BI4ZQ6MnnjA7b7/yCvzhDyVejzdQN6zndM08p2vmGV0vz+maeU7XLG/x8fF5/qwow3AngQtnIfpeGJRsVakSRESYlW7HjDF7JzVvDn//Oxw/bnd1IiIi4kWKEpY2A10AmjVr1hbYWSwVFadq1eDll01oevJJmDYNrr/e/PfcOburExERES9QlLC0DDjbrFmzL4BpwN+Kp6QS0KCB2ez0m2/gllvMztwtWpjFBTUJXERERPJR6DlLycnJmcAzxVhLyfvjH+GTT8zXsGFm64q2bc1KVLffbnd1IiIiUgZ55XYnYHZRDg4OpmXLlgQHB+NyuQr+5D//2WyA+u9/m/2+2rWDnj1hz57LP1dEREQqFK8MSy6Xi/DwcNxuN5Zl4Xa7CQ8P9yww+flBeLiZz/Tyy/Dxx2Zo7q9/hWPHSq54ERER8SpeGZaGDBlCWlpajra0tDSGDBni+cGqVjV3zO3eDf36wYwZ0LgxTJkCZ88WT8EiIiLitbwyLB3Lo+cnr/YCqV8fYmJgxw4zLDdsmFlu4L334ILNB6ViKdJwr4iIlAteGZZKVMuW8J//wNq1cPXV0Lu3mQS+caPdlUkpc7lcRERE5BjujYiIUGASEalgvDIs1apVy6P2QrnvPrPf3LvvwuHDcOed8Mgj8P33xXcOKdMiIyNJTU3N0ZaamkpkZKRNFYmIiB28MixNnz4df3//HG3+/v5Mnz69eE/k5wd9+5qANH48rFtnep4GD4affirec0mZs3//fo/aRUSkfPLKsOR0OnnnnXdwOBz4+PjgcDh45513cDqdJXPCwECIjDRLCwwcCNHRZiXwiRPht99K5pxiu6CgII/aRUSkfPLKsAQmMKWkpLBr1y5SUlJKLihdqF49mDkTdu40w3IjR0KzZjB/viaBl0NRUVEEBgbmaAsMDCQqKsqmikRExA5eG5ZsFRJiNuddvx7q1DH7zt1yC/z3v3ZXJsXI6XQSExOTowczJiamdIK5iIiUGQpLRXH33fDVVzBvnpnDdM898OCDkJhod2VSTGzpwRQRkTJFYamofH2hTx9ITjZzmDZsgBtvhGefhSNH7K5OREREishrw1KZWyywShUYPtxMAn/2WZg920wCj4qCi24/FxEREe/hlWGpTC8WWKeO2TLlu+/MWk2jRkHTphAbCxkZdlcnIiIiHvLKsOQViwU2awbLlpmVv6+91uw716aNWatJREREvIZXhiW32+1Ru606dICtW80ec7/8Ah07QpcupudJyrwyN9wrIiKlzivDkp+fn0fttvP1hccfh6QkmDwZvvgCbroJIiLgxx/trk7yUKaHe0VEpNR4ZVjKyGPuT17tZUZAAAwdCnv3mi1T3n0XmjSBsWPhzBm7q5OLeMVwr4iIlDivDEsOh8OjUEYsWwAAFOhJREFU9jKnVi14/XVISID774eXXjKh6d//1iTwMkR7w4mICHhpWIqKisp1I12v24bi+uvhgw9g82ZwOGDAAGjVCj75xO7KBO0NJyIihleGJQAfH598v/cqt99u5jG9/74ZjuvcGf78Z9ixw+7KKjTtDSciIuClYSkyMpK0tLQcbWlpad49l8THBx591AzNvfaa2Ubl5puhf384dMju6iok7Q0nIiLgpWGpXM8lueIK+NvfzCTwF16A+fPNfKYxY+DUKburq3C0N5yIiHhlWKoQc0lq1IApU8ymvA8+COPGmdAUEwPp6XZXJyIiUmF4ZViqUHNJ/vAHWLgQtmwxE8Kfftqs0bR6NViW3dWJiIiUe14ZlpxOJ3379s1ehNLPz4++ffuW7yGStm1h0yZYsgTS0qBrV7Ma+Ndf212ZiIhIueaVYcnlchEbG5u9CGVGRgaxsbHlf2VlHx/o3h127YI33oBvvoHWraFvXzhwwO7qyiVtdyIiIl4Zlir8ysqVK5sVwPfsgWHDYNEiaNoUIiPh5Em7qys3tN2JiIiAl4alcn03nCeuvhpefdXsOde9O7zyipnXFB0N58/bXZ3Xq/ChXEREAC8NSxXibjhPBAeDy2XWZgoJgUGD4I9/hBUrNAm8CBTKRUQEvDQsVai74TzRpg189hl8+KEJSQ89BHffDdu22V2ZV1IoFxER8NKwpJWV8+HjY9Zl2rkT/vUvMxn8llugTx9wu+2uzqsolIuICHhpWAKtrHxZ/v5mOG7PHhg50iw50KwZjBiBr1YCL5AKuUSFiIhcwmvDkhRQ9epm4vf330OvXvDqqzT+859hxgxNAr+MCrtEhYiI5KCwVFE0bAixsRAfz7lmzeD556FlS1i2TJPA86C74UREBBSWKp7QUPbPmQOrVkGlSmbJgTvugLg4uysrc9x5zPHKq11ERMonhaWKyMfHbJeyYwe89ZYZomvbFp54Avbts7u6MiNrrlJB20VEpHxSWKrIKlUyG/Pu2QOjR5slB5o3h6FD4Zdf7K7OdllzlQraLiIi5ZPCkkC1ajB2LOzeDU4nvPYaNG4Mr79uNu2toBwOh0ftIiJSPiksyf9r0ADmzIGvvzYb9P7tb9CiBSxeXCEngXfp0sWjdhERKZ8UluRSN90Ea9bARx9BlSrw6KPQrh1s2WJ3ZaVq9erVHrWLiEj5pLAkufPxgc6d4ZtvYPZsM/H79ttNcNq71+7qSoX2hhMREVBYksvx84P+/c18pn/+0/Q2hYSYIbpjx+yurkRpbzgREQGFJSmoK6+El14yoalfP3jjDbj+epg6Fc6ds7u6EqG94UREBBSWxFP160NMDHz7Ldx2m1lmoHlzWLiw3E0C14bNIiICCktSWDfcAKtXm4ng1aubBS3btoVNm+yurFhpw2YREVFYkqLp2BHi4+Hdd+HQIbN1SvfuZlVwERGRckBhSYrOzw/69jUBafx4WLvWbNL7/PPw8892V1ckLpeL4OBgWrZsSXBwMC6Xy+6SRESklCksSfEJDITISLN9yoABMHOmWQn81Vfh7Fm7q/OYy+UiPDwct9uNZVm43W7Cw8MVmEREKhiFJSl+9epBdLTZqPeOO2DECGjWDFwuyMy0u7oCGzJkCGkXbfeSlpbGkCFDbKpIRETs4LVhScMjXqBFC1i5Etavh9q1oU8f+NOf4LPP7K6sQI7lsY5UXu0iIlI+eWVYcrlc/F979x8jeV3fcfw5ByfnUMByNWJSdibo3YftVhBWWxoIVbBIFHKVRP5wCuhBB9QiHCGkMIiJOmC1iNfSxgzJXYR+VfxRKKHQ0lzRCoVC1tLIuXxSbHa3KNQK1h5M7k5g+8d3j+4tO5+7WfbmOzP7fCSE3c93+M47n3zYfe3n+5nPp16v7/V4pF6vG5j61bvfDY8+CrfdBj/9af79hg3wxBNFVyZJ0j4NZFhqNBq02+292trtNo1Go6CKtE+rVuUzSzHCDTfA/ffn2w987GN5gOpDa9eu7apdkjScBjIseWbXAHv96/M1TD/6EVxySb7B5VvfCtdfDwsCcNHOPffcrtolScNpIMPSkUce2VW7+tAb3wg33wzbt8Ppp+efogsBbr21bxaB33PPPV21S5KG00CGJQ2REOCOO+C7382PUrngAhgfh23biq6M6enprtolScNpIMPSc88911W7BsCpp8LDD8NXvwo//zm85z3w/vfnM08FOeigg7pqlyQNp4EMSyMjI121a0CsWpWfMffEE/CFL8CDD8Jxx8HFF8Mzz/S8nJdeeqmrdknScBrIsNRsNimXy3u1lctlms1mQRVpWa1ZA1demS8Cv/RS2LIlXwT+mc/ACy/0rAw/DSdJggENS7VajVarRaVSoVQqUalUaLVangg/bNauhS99CX74QzjzTLjuOli/Pg9Pzu5IknpkIMMS5IFpamqK7du3MzU1ZVAaZuvWwbe+BQ88AEcfDRdeCCecAPfdd0Df1rVxkiQY4LCkFejkk+Ghh+Ab34Dnn4f3vjefcfrBDw7I27k2TpIEhiUNmlIJPvhBmJyEL34RHnkE3v52uOgi+MlPlvWtXBsnSQLDkgbVIYfApk3w5JNw+eX5Zpbr1sGnPpXPOi0D18ZJksCwpEF35JFw4435dgNnnQWf/nQemm65BV588TXf3rVxkiTDkobDMcfA7bfna5re8hao1/PHc/feC7OzS75tlmVUq1XGxsaoVqtkWbaMRUuSBoFhScPlpJPge9+Db38bdu2C970PzjgDHnus61tlWUa9Xmd6eprZ2Vmmp6ep1+sGJklaYQxLGj6lEpxzTn5UyubN8P3vw4knwoc/DE89td+3aTQatNvtvdra7TaNRmOZC5Yk9TPDkobX614Hn/hEvhP4lVfC176Wb2p57bWwY8c+//OZmZmu2iVJw8mwpOH3hjfA5z8PMcIHPgDNZn58ype/nFwE7j5LkiQwLGklqVYhy/K9mY49Fj76UXjb2+DuuxddBO4+S5IkMCxpJXrnO+E734E774SXX4azz4bTToOJib1e5j5LkiQwLGmlKpVgwwZ4/HG4+eb83+94B5x3Hsxbk+Q+S5Ikw5JWttWr4eMfz3cCv/rq/MDe9evzr3/xi6KrkyT1AcOSBHDEEXD99fki8HPPhc99Ll8EfvPN8MtfFl2dJKlABxddgFaOW2/9N7Zs+deiy9gPp7HuxGP46H/czomXXsrMVU0+e8xH+PGvrS+6sK5t3HgC559/fNFlSNJAc2ZJWsS/H1bliuOu4o9/83J2HHwo4YX/LLokSVJBnFlSz5x//vEDOstxE6smJ7ludLToQiRJBXBmSZIkKcGwJEmSlLCkx3AhhCOAvwIOB14HXBFjfGg5C5MkSeoHS51ZugLYFmP8XeDDwF8sW0VSH8myjGq1ytjYGNVqlSzLii5JktRjS13gfROwa949di5POVL/yLKMer1Ou90GYHp6mnq9DuBO3pK0gpRmFzlAdL4QwoXApgXNH4kxPhpCOAq4F7g8xvjd+S+YmJiYXXgI6YGwc+dO1qxZc8DfZ5jYZ/vn9NNP5+mnn35V+5vf/Ga2bdtWQEWDxXHWHfure/ZZ9+yzztrtNuPj46XFru0zLHUSQngb8HXgyhjjvQuvT0xMzI6Pjy/p3t2YnJxk1I90d8U+2z+rVq1isf8/SqUSL7/8cgEVDRbHWXfsr+7ZZ92zzzqbmJjoGJaWtGYphPAbwDeBDy0WlKRhMDIy0lW7JGk4LXWB9w3AGmBzCOE7IYS/WcaapL7QbDZZ+Ci5XC7TbDYLqkiSVIQlLfCOMW5Y7kKkfrNnEXej0WBmZoaRkRGazaaLuyVphfG4EymhVqtRq9V8zi9JK5g7eEuSJCUYliRJkhIGNiy5s7J6wXEmSRrINUvurKxeyLKMjRs3snv3biAfZxs3bgQcZ5K0kgzkzFKj0XglKO3RbrdpNBoFVaRhdNlll70SlPbYvXs3l112WUEVSZKKMJBhaWZmpqt2aSmeffbZrtolScNpIMOSOytLkqReGciw5M7K6oW1a9d21S5JGk4DGZZqtRqtVotKpUKpVKJSqdBqtVx0q2W1efNmVq9evVfb6tWr2bx5c0EVSZKKMJBhCfLANDU1xfbt25mamjIoadnVajW2bt26VyjfunWrY02SVpiB3DpA6hWPO5EkDezMkiRJUi8YliRJkhIMS5IkSQmGJUmSpATDkpTgQbqSJD8NJ3Xggc2SJHBmSerIA5slSWBYkjrywGZJEhiWpI48sFmSBIYlqSMPbJYkgWFJ6sgDmyVJ4KfhpCTPhpMkObMkSZKUYFiSJElKMCxJkiQlGJYkSZISDEuSJEkJhiVJkqQEw5IkSVKCYUmSJCnBsCRJkpRgWJIkSUowLEmSJCUYliRJkhIMS5IkSQmGJUmSpATDkiRJUoJhSZIkKcGwJEmSlGBYkiRJSjAsSZIkJRiWJEmSEgxLUkKWZVSrVcbGxqhWq2RZVnRJkqQeO7joAqR+lWUZ9XqddrsNwPT0NPV6HYBarVZkaZKkHnJmSeqg0Wi8EpT2aLfbNBqNgiqSJBXBsCR1MDMz01W7JGk4GZakDkZGRrpqlyQNJ8OS1EGz2aRcLu/VVi6XaTabBVUkSSqCYUnqoFar0Wq1qFQqlEolKpUKrVbLxd2StML4aTgpoVarUavVmJycZHR0tOhyJEkFcGZJkiQpwbAkJbgppSTJx3BSB25KKUkCZ5akjtyUUpIEhiWpIzellCSBYUnqyE0pJUlgWJI6clNKSRIYlqSO3JRSkgR+Gk5KclNKSZIzS5IkSQmGJUmSpATDkiRJUoJhSZIkKcGwJEmSlGBYkiRJSjAsSZIkJRiWJEmSEgxLkiRJCYYlSZKkBMOSJElSgmFJkiQpwbAkSZKUYFiSJElKMCxJkiQlGJYkSZISDEtSQpZlVKtVxsbGqFarZFlWdEmSpB47uOgCpH6VZRn1ep12uw3A9PQ09XodgFqtVmRpkqQecmZJ6qDRaLwSlPZot9s0Go2CKpIkFeE1zSyFEI4F/gV4U4xx5/KUJPWHmZmZrtolScNpyTNLIYTDgRuBXctXjtQ/RkZGumqXJA2nJYWlEEIJaAHXAO19vFwaSM1mk3K5vFdbuVym2WwWVJEkqQil2dnZ5AtCCBcCmxY0TwNfjzHeFkKYAo5d+BhuYmJiduEvmgNh586drFmz5oC/zzCxz/bf3XffzU033cQzzzzDUUcdxaZNmzjrrLOKLmsgOM66Y391zz7rnn3WWbvdZnx8vLTYtX2GpcWEEJ4Enpr79iTgkRjjqfNfMzExMTs+Pt71vbs1OTnJ6OjoAX+fYWKfdc8+65591h37q3v2Wffss84mJiY6hqUlLfCOMb51z9dzM0tnLKkySZKkPufWAZIkSQmveVPKGGN1GeqQJEnqS84sSZIkJRiWJEmSEgxLkiRJCYYlSZKkBMOSJElSgmFJkiQpwbAkSZKUYFiSJElKMCxJkiQlGJakhCzLqFarjI2NUa1WybKs6JIkST32mo87kYZVlmXU63Xa7TYA09PT1Ot1AGq1WpGlSZJ6yJklqYNGo/FKUNqj3W7TaDQKqkiSVATDktTBzMxMV+2SpOFkWJI6GBkZ6apdkjScDEtSB81mk3K5vFdbuVym2WwWVJEkqQiGJamDWq1Gq9WiUqlQKpWoVCq0Wi0Xd0vSCuOn4aSEWq1GrVZjcnKS0dHRosuRJBXAmSVJkqQEw5IkSVKCYUmSJCnBsCRJkpRgWJIkSUowLEmSJCUYliRJkhIMS5IkSQmGJUmSpATDkiRJUoJhSZIkKcGwJEmSlGBYkiRJSjAsSZIkJRiWJEmSEgxLkiRJCaXZ2dkDcuOJiYkDc2NJkqQDYHx8vLRY+wELS5IkScPAx3CSJEkJhiVJkqQEw5IkSVLCwUUXsL9CCKuAvwSOB3YBF8UYn5x3/Q+Bi4EXgc/GGO8upNA+sh99thk4Bdgx17QhxviLnhfaZ0IIvw38SYzxXQvazwauIx9jW2KMtxRQXl9K9Nkm4CLgv+eaLo4xxh6X11dCCKuBLUAVOIT859Vd8647zubZj/5yjC0QQjgIuAUIwCxwSYzx8XnXHWNdGpiwBPw+sCbG+DshhJOAG4ENACGEo4BPAO8A1gAPhBD+Ica4q7Bq+0PHPpszDrw3xvizQqrrQyGEq4DzgBcWtK8GbgLeOXftwRDCXTHG/+p9lf2lU5/NGQfOjzFO9LaqvvYHwLMxxvNCCEcCjwF3geOsg479Nccx9mpnA8QYTw4hvAto8v+/Lx1jSzBIj+FOAf4OIMb4MHkw2uO3gAdjjLvmZkaeBI7rfYl9p2Ofzc06rQNaIYQHQwgbiymx7/wIOGeR9lHgyRjjz2OMu4EHgFN7Wln/6tRnkP8iuzqE8EAI4eoe1tTPvgl8cu7rEvlf93s4zl4t1V/gGHuVGOOdQH3u2wrwP/MuO8aWYJDC0uHA/EdEL4UQDu5wbQdwRK8K62OpPjsU+HPyv9rOBD4WQljxATPG+G3gl4tccox1kOgzgK8DlwCnAaeEEM7qWWF9Ksb4fIxxRwjhMOBbwLXzLjvOFthHf4FjbFExxhdDCF8h/zmfzbvkGFuCQQpL/wscNu/7VTHGFztcO4y9k/RKleqzNrA5xtiOMe4A/pF8bZMW5xjrUgihBHwpxvizub9g/xY4oeCy+kII4WjgfuC2GONX511ynC2iU385xtJijBcA64FbQgiHzjU7xpZgkNYsPUj+HPYbc+tvfjDv2iNAM4SwhnwB4Cjw+KtvseKk+mw9cHsI4QTy0HwK8JXelzgwJoF1c2smnieftv7TYkvqe4cDj4cQRsnXRpxGvlB3RQshvAm4D/ijGOO2BZcdZwvso78cY4sIIZwH/HqM8QbyP4xfnvsHHGNLMkhh6Q7g90II/0z+3PojIYQryJ+93hVC+DPge+S/+Bsxxp0F1tov9tVntwEPkz9CuTXGuL3AWvtSCOFDwK/EGFtzfff35GNsS4zxx8VW158W9Nk15DMCu4BtMcZ7iq2uL1wD/CrwyRDCnrU4twCHOs4Wta/+coy92l8DW0MI/wSsBi4HPhBC8GfZEnnciSRJUsIgrVmSJEnqOcOSJElSgmFJkiQpwbAkSZKUYFiSJElKMCxJkiQlGJYkSZIS/g8XUAxT1QeBNQAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
    " + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plots.one_way_anova_plot()\n", + "plt.show()" + ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "A one-way ANOVA has a log-linear counterpart called [goodness-of-fit](#goodness) test which we'll return to. By the way, since we now regress on more than one $x$, the one-way ANOVA is a **multiple regression** model.\n", + "A one-way ANOVA has a log-linear counterpart called [goodness-of-fit](#7.1-Goodness-of-fit) test which we'll return to. By the way, since we now regress on more than one $x$, the one-way ANOVA is a **multiple regression** model.\n", "\n", "The **Kruskal-Wallis** test is simply a **one-way ANOVA** on the rank-transformed $y$ (`value`):\n", "\n", - "$rank(y) = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2 + \\beta_3 x_3 +...$\n", + "$\\text{rank}(y) = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2 + \\beta_3 x_3 +...$\n", "\n", "This approximation is [good enough for 12 or more data points](https://lindeloev.github.io/tests-as-linear/simulations/simulate_kruskall.html). Again, if you do this for just one or two groups, we're already acquainted with those equations, i.e. the [Wilcoxon signed-rank test](#4.1-One-sample-t-test-and-Wilcoxon-signed-rank) or the [Mann-Whitney U test](#5.1-Independent-t-test-and-Mann-Whitney-U) respectively. " ] @@ -1446,10 +1437,111 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 29, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "num_points = 20\n", + "df = pd.DataFrame()\n", + "df[\"y\"] = np.concatenate([\n", + " np.random.normal(0.0, 1, num_points),\n", + " np.random.normal(1.0, 1, num_points),\n", + " np.random.normal(0.5, 1, num_points),\n", + "])\n", + "\n", + "df[\"group\"] = list(\"\".join([num_points * char for char in \"abc\"]))\n", + "df = df.join(pd.get_dummies(df.group, prefix=\"group\", drop_first=True).astype(np.float64))" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
    \n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
    ygroupgroup_bgroup_c
    01.080830a0.00.0
    1-2.316629a0.00.0
    2-0.138365a0.00.0
    32.003775a0.00.0
    42.051200a0.00.0
    \n", + "
    " + ], + "text/plain": [ + " y group group_b group_c\n", + "0 1.080830 a 0.0 0.0\n", + "1 -2.316629 a 0.0 0.0\n", + "2 -0.138365 a 0.0 0.0\n", + "3 2.003775 a 0.0 0.0\n", + "4 2.051200 a 0.0 0.0" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "df.head()" + ] }, { "cell_type": "markdown", @@ -1469,10 +1561,93 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 31, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "F, p = scipy.stats.f_oneway(df[df[\"group\"] == \"a\"].y,\n", + " df[df[\"group\"] == \"b\"].y,\n", + " df[df[\"group\"] == \"c\"].y)\n", + "\n", + "res = smf.ols(\"y ~ 1 + group_b + group_c\", df).fit()" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
    \n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
    valuep-valuest-values0.025 CI0.975 CI
    scipy.stats.f_onewayNaN0.015156NaNNaNNaN
    smf.ols (y ~ 1 + group_b + group_c)0.8018630.0121262.5911780.1821821.421545
    \n", + "
    " + ], + "text/plain": [ + " value p-values t-values 0.025 CI \\\n", + "scipy.stats.f_oneway NaN 0.015156 NaN NaN \n", + "smf.ols (y ~ 1 + group_b + group_c) 0.801863 0.012126 2.591178 0.182182 \n", + "\n", + " 0.975 CI \n", + "scipy.stats.f_oneway NaN \n", + "smf.ols (y ~ 1 + group_b + group_c) 1.421545 " + ] + }, + "execution_count": 32, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# FIXME what to tabulate here?\n", + "utils.tabulate_results([None, p, None, None, None],\n", + " res,\n", + " [\"scipy.stats.f_oneway\", \"smf.ols (y ~ 1 + group_b + group_c)\"],\n", + " coeff=\"group_b\")" + ] }, { "cell_type": "markdown", @@ -1494,10 +1669,16 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 33, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "_, p = scipy.stats.kruskal(df[df[\"group\"] == \"a\"].y,\n", + " df[df[\"group\"] == \"b\"].y,\n", + " df[df[\"group\"] == \"c\"].y)\n", + "\n", + "res = smf.ols(\"y ~ 1 + group_b + group_c\", df).fit() # TODO rank" + ] }, { "cell_type": "markdown", @@ -1850,6 +2031,36 @@ "\n", "4. Several named tests are still missing from the list and may be added at a later time. This includes the Sign test (require large N to be reasonably approximated by a linear model), Friedman as RM-ANOVA on `rank(y)`, McNemar, and Binomial/Multinomial. See stuff on these in [the section on links to further equivalences](#links). If you think that they should be included here, feel free to submit \"solutions\" to [the github repo](https://github.com/lindeloev/tests-as-linear/) of this doc!\n" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 11 Computing Environment" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "jupyter==1.0.0\r\n", + "matplotlib==3.1.0\r\n", + "numpy==1.16.4\r\n", + "pandas==0.24.2\r\n", + "patsy==0.5.1\r\n", + "scipy==1.3.0\r\n", + "statsmodels==0.10.0\r\n" + ] + } + ], + "source": [ + "!cat requirements.txt" + ] } ], "metadata": { diff --git a/utils.py b/utils.py index 0f4e69f..624cace 100644 --- a/utils.py +++ b/utils.py @@ -58,7 +58,7 @@ def tabulate_results(test_values, ols_results, names, coeff="x"): return table -def generate_toc(notebook="tests-as-linear.ipynb"): +def generate_toc(notebook="tests-as-linear.ipynb", max_header_levels=2): """ Generates a table of contents in Markdown. @@ -70,6 +70,9 @@ def generate_toc(notebook="tests-as-linear.ipynb"): ---------- notebook : str Path to notebook for which to generate a table of contents. + max_header_levels : int + Maximum number of header levels to show in table of contents (i.e. the + depth of headers to display). Returns ------- @@ -83,7 +86,7 @@ def generate_toc(notebook="tests-as-linear.ipynb"): for cell in cells: if cell["cell_type"] == "markdown": for line in cell["source"]: - match = re.search(r"^#+ ", line) + match = re.search(r"^[#]{{1,{0}}} ".format(max_header_levels), line) if match: level = len(line) - len(line.lstrip("#")) link = line.strip(" #\n").replace(" ", "-")