nhanes_confidence_intervals.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Confidence intervals case study using NHANES data\n",
    "\n",
    "\n",
    "This notebook demonstrates how to use Python and its statistical libraries to construct confidence intervals for proportions and means.  We will also cover some important points relating to the properties of confidence intervals, and discuss how to use and interpret confidence intervals in practice.  We will use the 2015-2016 wave of the [NHANES](https://www.cdc.gov/nchs/nhanes/index.htm) data for all the analyses below.\n",
    "\n",
    "It is important to note that the NHANES data are a \"complex survey\".  The data are not an independent and representative sample from the target population.  Proper analysis of complex survey data should make use of additional information about the manner in which the data were collected.  Since complex survey analysis is a somewhat specialized topic, we ignore this aspect of the data here, and analyze the NHANES data as if it were an independent and identically distributed sample from a population.\n",
    "\n",
    "The following module import statements and data reading statement are identical to what we have used previously in course 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now read the data into Python using the [read_csv](https://pandas.pydata.org/pandas-docs/stable/generated/pandas.read_csv.html) function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "da = pd.read_csv(\"nhanes_2015_2016.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confidence intervals for one proportion\n",
    "\n",
    "In this section, we demonstrate the construction of confidence intervals for the proportion of people who smoke.  The specific definition of \"smoker\" used here ([SMQ020](https://wwwn.cdc.gov/Nchs/Nhanes/2015-2016/SMQ_I.htm#SMQ020)) identifies a person as being a smoker if they self-report as having smoked 100 or more cigarettes in their lifetime.  It is more accurate to refer to this as a measure of \"lifetime smoking\" rather than \"current smoking\".  Recall that the definitions of these and other NHANES variables can be found using the NHANES code books, or by searching using the link below.\n",
    "\n",
    "https://wwwn.cdc.gov/nchs/nhanes/search/default.aspx\n",
    "\n",
    "We will calculate the proportions of smokers separately for females and for males.  Initially we can compare these two proportions and their corresponding confidence intervals informally, but later we will discuss methods to compare two proportions formally using confidence intervals.\n",
    "\n",
    "First we replace the numeric codes in the variables of interest with text labels, and set the rare answers other than \"yes\" and \"no\" to be missing (so they will automatically be omitted from all the analyses below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SEQN</th>\n",
       "      <th>ALQ101</th>\n",
       "      <th>ALQ110</th>\n",
       "      <th>ALQ130</th>\n",
       "      <th>SMQ020</th>\n",
       "      <th>RIAGENDR</th>\n",
       "      <th>RIDAGEYR</th>\n",
       "      <th>RIDRETH1</th>\n",
       "      <th>DMDCITZN</th>\n",
       "      <th>DMDEDUC2</th>\n",
       "      <th>...</th>\n",
       "      <th>BPXSY2</th>\n",
       "      <th>BPXDI2</th>\n",
       "      <th>BMXWT</th>\n",
       "      <th>BMXHT</th>\n",
       "      <th>BMXBMI</th>\n",
       "      <th>BMXLEG</th>\n",
       "      <th>BMXARML</th>\n",
       "      <th>BMXARMC</th>\n",
       "      <th>BMXWAIST</th>\n",
       "      <th>HIQ210</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>83732</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>62</td>\n",
       "      <td>3</td>\n",
       "      <td>1.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>...</td>\n",
       "      <td>124.0</td>\n",
       "      <td>64.0</td>\n",
       "      <td>94.8</td>\n",
       "      <td>184.5</td>\n",
       "      <td>27.8</td>\n",
       "      <td>43.3</td>\n",
       "      <td>43.6</td>\n",
       "      <td>35.9</td>\n",
       "      <td>101.1</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>83733</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>6.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>53</td>\n",
       "      <td>3</td>\n",
       "      <td>2.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>...</td>\n",
       "      <td>140.0</td>\n",
       "      <td>88.0</td>\n",
       "      <td>90.4</td>\n",
       "      <td>171.4</td>\n",
       "      <td>30.8</td>\n",
       "      <td>38.0</td>\n",
       "      <td>40.0</td>\n",
       "      <td>33.2</td>\n",
       "      <td>107.9</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>83734</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>78</td>\n",
       "      <td>3</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>...</td>\n",
       "      <td>132.0</td>\n",
       "      <td>44.0</td>\n",
       "      <td>83.4</td>\n",
       "      <td>170.1</td>\n",
       "      <td>28.8</td>\n",
       "      <td>35.6</td>\n",
       "      <td>37.0</td>\n",
       "      <td>31.0</td>\n",
       "      <td>116.5</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3 rows × 28 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "    SEQN  ALQ101  ALQ110  ALQ130  SMQ020  RIAGENDR  RIDAGEYR  RIDRETH1  \\\n",
       "0  83732     1.0     NaN     1.0       1         1        62         3   \n",
       "1  83733     1.0     NaN     6.0       1         1        53         3   \n",
       "2  83734     1.0     NaN     NaN       1         1        78         3   \n",
       "\n",
       "   DMDCITZN  DMDEDUC2  ...  BPXSY2  BPXDI2  BMXWT  BMXHT  BMXBMI  BMXLEG  \\\n",
       "0       1.0       5.0  ...   124.0    64.0   94.8  184.5    27.8    43.3   \n",
       "1       2.0       3.0  ...   140.0    88.0   90.4  171.4    30.8    38.0   \n",
       "2       1.0       3.0  ...   132.0    44.0   83.4  170.1    28.8    35.6   \n",
       "\n",
       "   BMXARML  BMXARMC  BMXWAIST  HIQ210  \n",
       "0     43.6     35.9     101.1     2.0  \n",
       "1     40.0     33.2     107.9     NaN  \n",
       "2     37.0     31.0     116.5     2.0  \n",
       "\n",
       "[3 rows x 28 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "da.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "da['SMQ020x'] = da.SMQ020.replace({1: 'Yes', 2: 'No', 7: np.nan, 9: np.nan})\n",
    "da['RIAGENDRx'] = da.RIAGENDR.replace({1: 'Male', 2: 'Female'})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now [tabulate](https://pandas.pydata.org/pandas-docs/stable/generated/pandas.crosstab.html) the numbers of female and male smokers and non-smokers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>RIAGENDRx</th>\n",
       "      <th>Female</th>\n",
       "      <th>Male</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SMQ020x</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>No</th>\n",
       "      <td>2066</td>\n",
       "      <td>1340</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Yes</th>\n",
       "      <td>906</td>\n",
       "      <td>1413</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "RIAGENDRx  Female  Male\n",
       "SMQ020x                \n",
       "No           2066  1340\n",
       "Yes           906  1413"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dx = da[['SMQ020x', 'RIAGENDRx']].dropna()\n",
    "pd.crosstab(dx.SMQ020x, dx.RIAGENDRx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The confidence interval (CI) is constructed using two inputs: the sample proportion of smokers, and the total sample size for smokers and non-smokers combined.  We calculate these values next."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Proportion</th>\n",
       "      <th>Total_n</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>RIAGENDRx</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Female</th>\n",
       "      <td>0.304845</td>\n",
       "      <td>2972</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Male</th>\n",
       "      <td>0.513258</td>\n",
       "      <td>2753</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           Proportion  Total_n\n",
       "RIAGENDRx                     \n",
       "Female       0.304845     2972\n",
       "Male         0.513258     2753"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dz = dx.groupby(dx.RIAGENDRx).agg({\"SMQ020x\": [lambda x: np.mean(x==\"Yes\"), np.size]})\n",
    "dz.columns = ['Proportion', 'Total_n']\n",
    "dz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Confidence intervals are closely connected to standard errors.  Recall that the standard error essentially tells you how\n",
    "far you should expect an estimate to fall from the truth.  A confidence interval is an interval that under repeated sampling covers the truth a defined proportion of the time.  In most settings, this \"coverage probability\" is set to 95%.  \n",
    "\n",
    "It turns out that in many settings, a 95% confidence interval can be constructed as the interval consisting of all points that are within two (or 1.96) standard errors of the point estimate.  More concisely, the confidence interval approximately spans from e - 2•SE to e + 2•SE, where e is the point estimate and SE is the standard error.\n",
    "\n",
    "Since the standard error plays such an important role here, we calculate it separately first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.008444152146214435\n",
      "0.009526078653689868\n"
     ]
    }
   ],
   "source": [
    "p = dz.Proportion.Female\n",
    "n = dz.Total_n.Female\n",
    "se_female = np.sqrt(p * (1-p) / n)\n",
    "print(se_female)\n",
    "\n",
    "p = dz.Proportion.Male\n",
    "n = dz.Total_n.Male\n",
    "se_male = np.sqrt(p * (1 - p) / n)\n",
    "print(se_male)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the standard errors for the estimated proportions of females and males who smoke are similar, and are each around 1% (since we are studying a proportion here, 0.01 corresponds to a 1 percentage point change in the smoking rate).  \n",
    "\n",
    "The standard error for a proportion is maximized when the true proportion is around 1/2, and gets smaller as the true proportion approaches either 0 or 1.  The estimated male smoking proportion is closer to 1/2 than the estimated female smoking proportion, and the male sample size is smaller than the female sample size.  Both of these factors lead to the male standard error being larger than the female standard error, although the difference is very small in this case.\n",
    "\n",
    "Next we calculate the 95% confidence intervals for the proportions of female and male smokers using the formula for the one-sample confidence interval for a proportion:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.288294683866098 0.32139576027925865\n"
     ]
    }
   ],
   "source": [
    "p = dz.Proportion.Female\n",
    "n = dz.Total_n.Female\n",
    "lcb = p - 1.96 * np.sqrt(p * (1 - p) / n)\n",
    "ucb = p + 1.96 * np.sqrt(p * (1 - p) / n)\n",
    "print(lcb, ucb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results above indicate that any population proportion (for female lifetime smokers) between 0.288 and 0.321 would be compatible with the data that we observed in NHANES."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.49458714955108174 0.531929377873546\n"
     ]
    }
   ],
   "source": [
    "p = dz.Proportion.Male\n",
    "n = dz.Total_n.Male\n",
    "lcb = p - 1.96 * np.sqrt(p * (1 - p) / n)\n",
    "ucb = p + 1.96 * np.sqrt(p * (1 - p) / n)\n",
    "print(lcb, ucb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These results indicate that any population proportion (for male lifetime smokers) between 0.493 and 0.531 would be compatible with the NHANES data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a routine data analysis, we do not need to calculate these intervals manually.  We can use the Statsmodels library to calculate the CI for us in one line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.2882949879861214, 0.32139545615923526)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 95% CI for the proportion of females who smoke (compare to value above)\n",
    "sm.stats.proportion_confint(906, 906+2066)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.49458749263718593, 0.5319290347874418)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 95% CI for the proportion of males who smoke (compare to value above)\n",
    "sm.stats.proportion_confint(1413, 1413+1340)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confidence intervals comparing two independent proportions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The confidence intervals for the proportions of female and male smokers shown above are quite narrow and do not overlap.  This suggests that there is a substantial difference between the lifetime smoking rates for women and men. However there is no explicit information here about how different the two population proportions might be.  To address this question, we can form a confidence interval for the difference between the proportion of females who smoke and the proportion of males who smoke.  \n",
    "\n",
    "The point estimate of the difference between female and male smoking rates is -0.208 (0.305 - 0.513).  That is, the smoking rate is about 20 percentage points higher in men than in women.  This difference of around 20 percentage points is only a point estimate of the underlying true value -- it is not exactly equal to the difference between the unknown proportions of females and males who smoke in the population.  A confidence interval helps us assess how far the estimated difference may be from the true difference.\n",
    "\n",
    "As above, we start with the standard error.  The difference between two sample proportions based on independent data has a standard error that reflects the combined uncertainty in the two proportions being differenced.  This standard error can be calculated very easily.  If SE1 and SE2 are the standard errors for two proportions, then sqrt(SE1^2 + SE2^2) is the standard error for the difference of these proportions (`sqrt` is the square root function).  Note that this formula is only accurate if the two sample proportions being differenced are independent.\n",
    "\n",
    "In the next cell we calculate the standard error for the difference between the proprotion of females who smoke and the proportion of males who smoke."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.012729881381407434"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "se_diff = np.sqrt(se_female**2 + se_male**2)\n",
    "se_diff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The standard error of around 0.013 indicates that the estimated difference statistic -0.208 is expected to fall around 0.013 units from the true value.  We do not know in which direction the error lies, and we do not know that the error is exactly 0.013, only that it is around this large on average.  For most purposes, a standard error of 0.013 relative to an observed difference of -0.21 would be considered very small.  That is, we have a very accurate estimate of the difference between smoking rates in women and in men.\n",
    "\n",
    "Now that we have the standard error, we can construct a 95% confidence interval for the difference in proportions by taking the estimate and subtracting and adding two (or 1.96) standard errors from it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.2338728044024504 -0.18295327887682067\n",
      "-0.20841304163963553\n"
     ]
    }
   ],
   "source": [
    "d = dz.Proportion.Female - dz.Proportion.Male\n",
    "lcb = d - 2 * se_diff\n",
    "ucb = d + 2 * se_diff\n",
    "print(lcb, ucb)\n",
    "print(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 95% confidence interval above shows us that any value for the difference of population proportions (between females and males) lying between -0.233 and -0.183 is consistent with the observed data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Confidence intervals for subpopulations\n",
    "\n",
    "Since smoking rates vary strongly with age, it might be more informative to stratify the data into homogeneous age bands and compare the proportions of female and male smokers within each age band.  We can also calculate the 95% confidence interval for this difference within each age band.  These data can be displayed as a plot, with the difference in proportions plotted as a curve.  The confidence intervals can then be used to construct a \"confidence band\" around the estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEKCAYAAAA1qaOTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsnXdYVEfbh++hWVFUmiJ2RWNX7L13jUaNGnuMqSZ5ExPTi1/eaExMNOWNJjH2XpJoLMSCvQtiV2woNkCKgHTm++OsG5GlKHt2F5j7uvaCPWfOPg9l93dm5ilCSolCoVAoFHnFztoOKBQKhaJgoARFoVAoFGZBCYpCoVAozIISFIVCoVCYBSUoCoVCoTALSlAUCoVCYRaUoCgUCoXCLFhFUIQQZYUQW4UQwYavZbIYN8YwJlgIMeah4zuFEOeFEMcND3fLea9QKBQKUwhrJDYKIWYAkVLK6UKI94AyUsopj4wpCxwFfAEJHAOaSimjhBA7gclSyqMWdl2hUCgUWeBgJbsDgI6G7xcCO4Epj4zpAWyVUkYCCCG2Aj2B5U9q1NXVVVapUuVJL1coFIpCybFjxyKklG45jbOWoHhIKW8BSClvZbFk5QVcf+h5qOHYA+YLIdKAtcAXMhdTrSpVqnD0qJrUKBQKxeMghAjJzTjdBEUIsQ3wNHHqw9y+hIljD0TjOSnlDSGEM5qgjAIWZeHHRGAiQKVKlXJpWqFQKBSPi26CIqXsmtU5IcQdIUR5w+ykPBBmYlgo/y6LAVREWxpDSnnD8DVWCLEMaE4WgiKl/AX4BcDX11dVwlQoFAqdsFbY8HrgQdTWGOAvE2P8gO5CiDKGKLDugJ8QwkEI4QoghHAE+gKnLOCzQqFQKLLBWoIyHegmhAgGuhmeI4TwFUL8BmDYjP8/4IjhMdVwrAiasJwAjgM3gF8t/yMoFAqF4mGsEjZsLXx9faXalFcoFIrHQwhxTErpm9M4lSmvUCgUCrOgBEWhUCgUZsFaeSgFlqDr0Sw6EMKZW/coVdSBfg0rMLhpRYo62lvbNYVCodAVJShmZOmhED764xQP70oduhLJ6qPXWTKhBc5FHa3mm0KhUOiNWvIyE1cj4vn4z4xi8oCg0Bi+8TtvcZ8UCoXCkihBMROrjl4nPZuAubUBN0hKTbOcQwqFQmFhlKCYiZDI+9mej0tKJTI+2ULeKBQKheVRgmIm3J2LZHvewQ5cijlZyBuFQqGwPEpQzMQzTSpme764kwNR99UMRaFQFFyUoJiJel6lealD9SzP30tMZdD/9nPu9j0LeqVQKBSWQwmKGZnS04cfhjfGt3IZSjjZ41mqKM+1qEQdT2cAbt9LZMicAxy4dNfKnioUCoX5UbW8LEB8UiqvLgtg5/lwAJzs7fhmaEP6N6xgcV8UCoXiccltLS+V2GgBShRx4NfRvnyw7iSrj4WSnJbO68sDCbuXyIR21aztnqKQcC8xhRWHr7HtbBjJqen4Vi7D6FZVqFSuuLVdUxQQ1AzFgkgp+W5bMN9vDzYeG9+mKh/1qYOdnakGlQqFebgZncCwXw5y7ZHw9mKO9vw62pe2NV2t5JkiP6CqDdsgQgje6laLaYPq80A/ft93hUkrAklMUUmPCv14f93JTGICkJCSxmvLAkhIVv9/iryjBMUKDG9eiV9H+1LMUDBy44lbjP79MDH3U6zsmaIgEhp1n10XwrM8H52QwuZTtyzokaKgogTFSnSp48HyiS0pW0JLdjx8JZLBc/ZzMzrByp4pChKxiSksP3wtx3FXI+It4I2ioKMExYo08nZh3cutqWzYFA0Oi2Pg//Zx9pbKVVE8OZHxyaw6cp3xC47Q9P+28ZP/pRyvccuh0oNCkRtyjPISQtQC3gEqPzxeStlZR78KDVVcS7D25dY8v+AIQaEx3LmXxNA5B5g7qimta6iNUkXuuB2TiN/p22w5dZtDV+5mW6j0URzsBX0aqBB2Rd7JTdjwamAO8Cugdu50wLVkEZZPbMmrSwPwPx9ObFIqY+Yf5pshDRnQyMva7ilslJC78Ww5dZstp28TeC3a5Jg65UvRs64n1dyKM2XtCe4np2caU7qoI04OarFCkXdyDBs2hIs1tZA/umLtsOGcSE1L56M/T7HiyHXjsQ961+aFdtUQQoUVF3aklJy/E6uJyKnbnLsda3Jck0ou9KznSY+6nlQuV8J4/GJYHP/beZGtZ+6QlJJGiSIORBkCQQY18eLboY0s8nMo8h+5DRvOjaB8BoQBfwBJD45LKSPz4FxZYCVQBbgKDJVSRpkYtwVoCeyVUvZ96HhVYAVQFggARkkpc6y8aOuCAtqHxuztwcza9m+uytjWVfi471PYq1yVQoeUkqDQGLacuo3f6dtcMbF5bm8naFG1LD3redL9KU88SxfN1WtHxCXRa/YewmO1t/WsZxvxdGM1I1ZkxpyCcsXEYSmlfOIUbyHEDCBSSjldCPEeUEZKOcXEuC5AceDFRwRlFbBOSrlCCDEHCJJS/pyT3fwgKA9YeeQaH/xxijTDYnivep5892wj1Zu+EJCWLjlyNdIoIrdiEjONcbK3o11NV3rU86RrHQ9jtODjsic4nFHzDgNQsogDG19vm2FWo1CAGQVFD4QQ54GOUspbQojywE4ppU8WYzsCkx8IitDWfsIBTyllqhCiFfCZlLJHTnbzk6AA+J8L45WlASQYkh6bVynLL6Ob4lJc9VUpaCSlprH/0l38Tt1m65k73DXRjK24kz2dfNzpUc+TTj5uOBd1NIvtaZvPMnfXZQAaVizN6pdaqz0VRQbMVstLCOEIvAy0NxzaCcyVUuYlC89DSnkLwCAq7o9xbTkgWkqZangeChTIeXqn2u6smNiS8QuOcDc+mcNXIxk85wALxjWjYhlVfym/cz85lV3nw9ly+jY7zoYRm5SaaUzpYo50reNBz3qetKvpqssM9e1uPhy8dJeg0BiCQmP4dusF3utV2+x2FAWf3ER5/Qw4Av8zPB9lODYhu4uEENsATxOnPnwcB029tIljWU6zhBATgYkAlSpVyqNpy9PQ24V1r7RmzO+HuXr3PhfD4hj0v/0sGNecpyqUsrZ7isckJiGFHefusOXUbXZdCCcxJXPUlWvJIvSoq4lIy2rlcLTXd7bg5GDH98Mb03v2HuKT05i7+xJta7iq+l6KxyY3eyhBUsqGOR17LKNqyeuxuRuXxPiFRwm6roWHliziwNxRTWmjclVsnoi4JLaeucPmU7fZfzGCVBNJIl4uxehVz5Oe9TxpXKmMVQIw/ggM5T8rgwAt0XHLG+0oV1IlPCrMW74+TQhRXUp5yfDC1ch7Psp6YAww3fD1r9xeKKWUQgh/YDBapNdjXZ9fKVeyCMtfaMGkZYFsPxdGXFIqY+cf5uvBDVVkjg1yIzoBP0OOyNGrkSYTDWu4l6RnXU1E6lYoZfXQ8IGNK7LnQgTrAm8QHpvEO2tOMG+Mr9X9UuQfcjND6QLMBy6jLTdVBsZJKf2f2KgQ5YBVQCXgGjBEShkphPAFXpJSTjCM2wPUBkoCd4HnpZR+BlF7EDYcCIyUUiaZMJWB/DxDeUBqWjof/3U6Q32mKT1r81IHlatibS6Hx7HZEJl1IjTG5Jh6XqXoVa88Pep6UMPd2cIe5kxcUip9vt9DyF2tMvGn/Z5iXJuqVvZKYW3MGuUlhCgC+KAJyrncfHjbIgVBUEDLTfhxx0Vmbr1gPDamVWU+6VdX5apYECklZ27dM85ELtyJyzRGCPCtXIYedbVEQ++yth9MEXQ9mmd+3k9qusTJ3o4/Xm1N3Qqlre2WworkWVCEEJ2llDuEEINMnZdSrsujjxanoAjKA1Ydvc77604ac1V61PVg9rDGKlflCYlPSmXJwRD+PH6TmPvJ1PRwZnSrynSp42Eck54uCbwebaybZarHiIOdoFX1cvSs50m3pzxwd85doqEtMXfXJaZtPgdANbcS/D2pLcWdVIPXwoo5BOVzKeWnQoj5Jk5LKeX4vDppaQqaoADsPK/lqtw3NEjyrVyGX0f7UuYJE90KK/cSUxj+y0FO38xc6fnVjtVpXcPVmGgYFpt5gl7EwY72tdzoWdeTLnXc832uUHq6ZMz8w+wJjgBgWDNvpj/TwMpeKayFOTPlq0opr+R0LD9QEAUF4GRoDOMWHCYiTkuGq+ZWgoXjmueL5RVbYeqGM/y+7/H+pUs42dO5jge96nnSoZYbJYoUrDv4sHuJ9Jq9x5hk+eOIxvRVVYkLJeYUlAApZRMTL57vCkZaSlDi4+O5fPkyLi4ueHt7624P4Nrd+4yZf9hY68nNuQgLxjVTa9+5IC1d0njqP9xLzJxY+ChlijvS7SktR6R1dX0SDW0J/3NhjFtwBADnog5ser2dulEphOS5p7wQorYQ4hmgtBBi0EOPsUD+WxS2AAkJCbz55pt4enrSoEEDKlWqRJs2bThw4IDutiuVK87al1vTuJILAOGxWl+VPcFZt35VaFyJiMtRTLxcirLshRYc+bArMwY3pHNtjwIvJqBVaxhviPKKTUzlzZXHSU3LnIypUED2HRt9gL6AC9DvoUcT4AX9XctfpKenM2jQIGbPnk1c3L/RPvv376dz584cPnxYdx/KlnBi2YSWdDVsIscnpzFu/hHWBYTqbju/kZKWzj+nbzNh4RF6zNqT4/gOPu60ru6Kg85Z67bIlF4+PFVeq8pwLCSK77cH53CForCS7ZKXEMIemCKl/NJyLumHnkteW7ZsoVevXlme79q1K1u3btXF9qOkpqXz6frTLD30b67KOz18eKVj9UKfq3I5PI5VR0NZGxBqLNueG/56tQ0NvV109My2uRQeR9/v95KQkoYQsPyFlrSsVs7abiksRJ6XvACklGlAN7N5VYBZs2ZNtue3bdtGdLTprnrmxsHeji+ersc7Pf6tZvO133k+/uvfcviFifvJqaw5FsrQOQfoPHMXc3ZdyiAmtT2dmdLTBx/Pkiavn9S5RqEWE4DqbiX5fEBdAKSE/6w8TpSJisiKwk1uwlL2CyF+RGuIZezuI6UM0M2rfEh8fObGR49y//59XFws88EkhODVTjXwKFWU99aeIDVdsuTgNcLuJfH98IKfq/KgMdXKI9fZEHSTuEcq+ToXcaB/owo828yb+l6lEUIwulUVlhwM4a/jN4k25KGMaV2ZzrU9srBSuBjStCK7L4Tz94lb3IpJZMraE8wd1bTQz3oV/5KbKC9TJVaklLKzPi7ph55LXjNnzmTy5MlZnvf09OTGjRvY2Vl+DX73hXBeXnKMeEOuSpNKLswb06xA5qpExSfzR+ANVh29brJFbvOqZXnW15ve9ctTzKlgi6oexCSk0Of7PYRGJQDwf0/XY1TLylb2SqE3Nt1gy1roKSh3796lRo0aWS5rlStXjsDAQIuFET/KqRsxjJ1/hIg4bamnIOWqpKdL9l6MYOXR62w9fYfkR6KQ3JyLMLhpRYb6elPVVXUjzCvHQqIYOvcAaemSIg52rH+tLT6etleXTGE+zJmHUhr4lH8bbO0CpkopTVe/s2H0zkPZu3cvAwYMIDIy0uT5atWq4e/vb7W+LNcjtVyVy+Ha8pxrSS1XpZ5X/sxVCY26z+qjoaw5FsqN6IQM5+ztBJ183BnWzJuOPm6FMjpLT37cEcw3/2i15Gp5lGT9a20L/DJqYcacgrIWOAUsNBwaBTSUUpqs8WXLWCKxMTY2lqVLl3LixAlKly5Nr169mDx5MkeOaMlhVapUwd/fnypVqujqR1ZExSfz/MIjBFzTZlIlnOz538imdKjlZhV/Hpek1DS2nrnDyiPX2Xsxgkf/fau6lmCorzfPNPHCvZRKl9KLtHTJc78d5OBl7eZpZMtKfPF0fSt7pdALcwrKcSllo5yO5QesVXolOjqanj17cujQIUDrHLlz506qVrVOWfDElDReXx7IP2fuAFoxw+nPNGBw04pW8Sc3nLt9j5VHrvNn4A2i7mfsPl3U0Y7e9cvzrK83zauWVZvEFuJ2TCI9Z+8m2vD3mDuqKT3qmmrSqsjvmFNQDgDvSCn3Gp63Ab6RUrYyi6cWxJq1vGJiYujZsycHDx4EwNvbG39/f6pXr24Vf9LSJZ+tP83igyHGY5O71+LVTjVs5gM5NjGFDUG3WHn0urFT5cM0rFiaoc286dewAqWKOlrBQ8U/p28zcfExAEoXc2TzG+2o4FLMyl4pzI05BaUR2nJXabR+KJHAGCnlCXM4akmsXRzy3r179OrVi/379wNQsWJF/P39qVGjhlX8kVLy865LzNhy3nhsRItKTO1f12p7DlJKjlyNYuWR62w6eYuElIzNQV2KOzKwsRfPNvOmtmcpq/ioyMgnf51i0QHtxqRF1bIse6Gl6stTwDB7lJcQohSAlDJzfe98grUFBbQ9lt69e7N3714AKlSogL+/P7Vq1bKaT+sCQnl3zQljr/OudTz4YXhji4bVhsUmsi7gBquOXOdyRMacHiGgbQ1Xhvp6072uB0Uc1OavLZGYksbTP+0zhmm/3a0Wk7rUtLJXCnNizhlKObQor7aABPaiRXndNYejlsQWBAUgLi6OPn36sHv3bgDKly/Pjh07qF27ttV82hMczstLAowJgI0NuSpldcxVSU1LZ+f5cFYevc6Oc2GZsvi9XIoxuGlFhvhWpGKZ/B/eXJAJvhNL3x/2kpSajr2dYOXElvhWKWtttxRmwpyCshXYDSwxHHoO6Cil7JpnLy2MrQgKaJn1/fr1w99fyxv18PDA39+fOnXqWM2n0ze1XJUHZUmqumq5KpXKmffD/EpEPKuOXmftsdBMzaoc7QXdn/Lk2WbetKnhqpZO8hFLD4Xw4R+nAO1mYNMb7ShdTO1tFQTMKSiZep8IIY7m5sVtDVsSFNBKsfTv35/t27cD4O7uzo4dO6hbt67VfLoeeZ+x8w9zyZir4sT8sc2pXzFvuSoJyWlsPnWLlUeuc+hK5jwdHw9nhjbzZmBjL11nRQr9kFLy8pIAtpy+DUCf+uX5cURjmwnyUDw55hSUb4CjwCrDocFAXSnlp3n20sLYmqCA1kNlwIABxkrEbm5ubN++nfr1rRfTH30/mQkLj3I0JAqA4k72/PRcEzr5uD/W60gpOXlDq6e1/vhNYh+pp1WyiAP9Gmr1tBpWLK0+eCxEYmIiqamplCxpuhhmXoi+n0yv2Xu4FZMIwFfP1OfZZtZJ5FWYD3MKSixQAnhQz8KOf4tESinlY4faCCHKohWbrAJcBYZKKaNMjNsCtAT2Sin7PnR8AdABeJCtP1ZKeTwnu7YoKKCJysCBA/Hz8wPA1dWVbdu20bBhQ6v5lJiSxpsrjhvvNu3tBNMG1WeorzeXw+O4EZ1A+dJFqeGeueRG9P1k/gy8wYojputpNatShqG+3vRpUJ7iTgWrba4tc+jQIT777DP8/PyQUtKgQQOmTJnCiBEjzGvn8l2G/3qQdAnFHO3ZMKktNdzNL14Ky2HTtbyEEDOASCnldCHEe0AZKeUUE+O6AMWBF00Iyt9Syuxrxj+CrQoKaHeNzzzzDJs2bQKgbNmybN++nUaNrJc/mpYumbrhNAsP/Jur4uVSLEOZk8aVXJg2qD613J3Zf+kuK49ex+/0bZJTM9bTci3pxDOGelrV3dSHi6Xx9/enZ8+eJCdnLjn/5Zdf8v7775vV3ndbLzDb0IirTvlS/PFKa1WaJR9jVkERQvTn31peO6WUf+fRufNoG/u3hBDlDa/pk8XYjsDkgi4oAElJSQwePJi//9Z+vWXKlGHbtm00adLEaj5JKZm7+zLTN5/LckxRBzvKlHAyLnM8wE5AJx93hjbzpnNtdxxVPS2rIKWkXr16nDlzxuR5BwcHQkJCqFChgtlspqalM+yXg8Zl03FtqvBpP+vtDSryhlkabBleaDrwBnDG8HjDcCwveEgpbwEYvj7e4rzGf4UQJ4QQ3wkhimQ1SAgxUQhxVAhxNDzctvurFylShDVr1tC/f38AoqKi6NKlC9YUQSEEL3WoTrensv4TJaamZxCTyuWK804PHw6834V5Y5vRo66nEhMrEhQUlKWYAKSmprJ69Wqz2nSwt2PWsEaUKqotac7fd5Ud5+6Y1YbC9sjNu7w30E1K+buU8negp+FYtgghtgkhTpl4DMir08D7QG2gGVAWyLRc9gAp5S9SSl8ppa+bm+0XQCxSpAirV69m4MCBgFYHrGvXrhbpSZ8d1+4m5DhmYGMvVkxsyc7JHY3NvRTWIzw8nAULFjBp0qQcx969a/60soplijP9mQbG55NXnyDsXmI2VyjyO7m9bXy4zWCu4kellF2llPVMPP4C7hiWujB8DXscp6WUt6RGEjAfaP4419s6Tk5OrFy5kmeeeQbQ6oB169bNWAfMGqQ80mPkUVxLOvHds41oWa2citayElJKzp07x4wZM2jbti0eHh6MGzfOWJUhO/QKVe9dvzzDm2s9gCLjk/nPquOkF8I21IWF3AjKNCBQCLFACLEQOAZ8mUe764Exhu/HAH89zsUPiZEAnkYrr1+gcHR0ZPny5QwZMgTQ6oB1797dWAfM0vhWKZPt+RZVy1nIE8XDpKamsnv3biZPnoyPjw916tRhypQp7Nu3j4f3R52ds26A5eLiwtNPP62bj5/0rWuM8tp38S5zd1/WzZbCumQrKIYP7L1oobvrDI9WUsoVebQ7HegmhAgGuhmeI4TwFUL89pD9PcBqoIsQIlQI0cNwaqkQ4iRwEnAFvsijPzaJo6Mjy5YtY9iwYYBWB6xHjx65uuM0N+PbVsUhi6x1eyGY0M46pfgLI7GxsaxZs4bRo0fj4eFBhw4dmDlzJsHBwcYxdnZ2tG/fnm+++Ybz589z6dKlLCMGk5KSCAkJMXnOHBRzsuf7YY1xctA+bmb+c57jJqpHK/I/T5Qpn1+x9SivrEhNTWXs2LEsXboUgBIlSrBp0ybat2+fw5XmZeuZO7y16jixif8mKBZ3smf6Mw3o39B8EUKKzFy/fp0NGzawfv16/P39TYb/lixZkp49e9K/f3969+5NuXIZZ41JSUmsXLmS9evXk5SURHx8vLH0T+PGjTlw4ABFimQZ35JnFu6/yqfrTwNQqWxxNr7eFmfVdiBfYM7Exp+ABVLKI+ZyzlrkV0EBSEtLY/z48SxatAiA4sWLs3HjRjp27GhRP+KTUtl08hY3ohPwLFWU3g3Kq14kOiClJDAwkPXr17N+/XoCAwNNjqtYsSL9+/enf//+dOzY8bEEITU1lfbt23PgwAEA/vOf//Dtt9+axX9TSCl5YdFRtp3VtkyfblSBWcMa62ZPYT7MKShnAB+0jPZ4tJ4oUkrZILvrbJH8LCigicqECRNYsGABAMWKFePvv/+mc+fO1nVMYRaSkpLw9/dn/fr1bNiwgdDQUJPjmjRpYhSRRo0a5SkI4urVqzRq1IiYGK3oxMaNG+ndO8cgzicmMj6ZnrN2G4uCfju0IYOa2G6nUIWGOQWlsqnjUkr9Fl11Ir8LCkB6ejoTJ05k3rx5ABQtWpQNGzbQtWu+K/6sACIiIti0aRPr16/Hz8+PuLi4TGOcnJzo3Lkz/fv3p2/fvnh7e5vVh9WrVzN06FBAqyUXFBRE+fLlzWrjYfZfjOC5eYeQEko42fP36+2o6lpCN3uKvGPuTPkm/NsPZZ+UMiDvLlqegiAooInKyy+/zC+//AJoovLXX3/RvXt3K3umyA0XLlwwLmXt27eP9PTMIdnlypWjT58+9O/fn+7du2cbpWUOXnjhBX77TYuH6dq1K35+ftjZ6ZeMOmPLOf638xIA9b1Ks/bl1sZNe4XtYc4ZyifAELQIL9DCdFdLKfNdZFVBERTQROW1117j559/BrSEyD///JOePXta2TPFo6SlpXHgwAGjiJw/f97kuJo1azJgwAD69+9Pq1atcHCwXOHM+Ph4fH19OXdOK7Hz1Vdf8e677+pmLyUtnSFzDhijvV5sX433e1uvF5Aie8wpKGeBxlLKRMPzYkCAlDLf/fUtISiJKWn4nb7NmZv3KFXMkT71y1NFp+m8lJJJkybx008/AdrSyLp16+jTp48u9hS5Jy4ujn/++Yf169ezceNGIiIiMo2xs7OjdevWxv0QHx+T5ewsRlBQEM2bNyc5ORkHBwf27dtH8+b65Qxfj7xP79l7jG0NFo1vTvtatl/NojBiTkHZDAyXUkYbnrsASx4u1phf0FtQTobGMGHREe7cy9iF8IV2Vfmgdx1dMsillLz55pt8//33gJa7snbtWvr162d2W4rsuXHjhjG0d/v27SZDe0uUKEGPHj2Mob22Vg7ohx9+4PXXXwegWrVqBAYGUqrUY3eoyDV/Hb/BGyu0zhOuJYuw5c12uJbUL3RZ8WSYU1D+RKuZtRVtD6UbWrJjGICU8vU8e2sh9BSUe4kpdPx6J5HxmT9EAD7vX5cxravoYltKydtvv813330HaKKyevVqBgwwR9m0wkdKSgr379+nVKlS2d4ESCkJCgoyLmUdO3bM5LgKFSoYZyGdOnWiaFHbrXEmpaR///7GitcjRoxgyZIlupbTmbw6iDXHtIi2DrXcmD+2GXaq9bNNYU5BGZPdeSnlwsf0zWroKSjz913h8w1ZV3StWKYYu9/ppNsbRUrJu+++yzfffANoJclXrlzJoEGDdLFXEAkJCeHjjz9m1apVJCUlUbFiRV599VUmT55s3M9ISkpi165dRhG5fv26yddq1KiRUUSaNGmSr+qbRURE0KBBA27dugXAwoULGT16tG724pNS6fvDXq5EaH37PupThwntqulmT/H42HSDLWuhp6C8vjyQ9UE3sx1z6IMuulbglVLy/vvv89VXXwFgb2+foR6YImtCQkJo2bIlt2/fznSuX79+DBkyhA0bNrBlyxZiYzN3oXR0dKRTp07079+ffv36UalS/m57u2PHDrp27YqUkhIlShAQEECtWrV0s3fqRgwD/7ePlDSJo73gj1faUM8rV3VoFRbAbP1QFLmjSC5CHnMzJi8IIZg2bRoffPABoEUXDR8+nJUrV+pqtyDw8ccfmxQTgA0bNjB69GhWr16dQUzKlCnDyJEjWbVqFREREfj5+fHqq6/mezEB6Ny5s7GLY3x8PMOZMGcxAAAgAElEQVSHDze5J2Qu6nmVZkrP2gCkpEleXx5IfFJqDlcpbA0lKGaie13PbM838CqNS3En3f0QQvDFF1/w8ccfA5qojBgxgmXLluluO7+SkpLCqlWrcjW2evXqvPXWW+zcuZOwsDAWL17MkCFDdN24thafffYZrVq1AiAgIMDsbYIfZXybqnQwRHldjojnM0PdL0X+QQmKmehc253mVctmef5GdAI3o3NuUmUOhBBMnTqVzz77DNByVkaNGsWSJUssYj8/kZKSwtq1a0lKSsp2nI+PD2fOnCE4OJiZM2fSoUMHi+aJWIMH1a4fiOW3337L5s2bdbNnZyf4ZkhDY5TX6mOhOS4jK2yL3LQA3iCEWP/IY7EQ4g0hhO2Gq1gYezvB72ObMdS3Io72/27AFneyB+BufDKj5h3KMgpMDz799FOmTp0KaKIyevRoFi7MNzEUunLq1CkmT56Mt7c3w4cPz3F8//79qVNHn9BvW6ZKlSrGigwAY8aMyXJp0By4ORdh5tCGxucfrjvJ9cj7utlTmJfcRHnNBtyA5YZDzwK3gWJAKSnlKF09NCOWypSPik/mckQ8pYo64Fm6KCN/O0RQqFZ8r6G3C8smtKBEEcvd3X755Zd8+OGHgDZ7+e233xg/frzF7NsKkZGRLF++nPnz52cZ4msKBwcHzp49S40aNXT0zraZMGGCsX5ct27d2LJli66lWb7cdJZfDI24GldyYdWLrXC0Vwsq1iK3m/JIKbN9ALuzOgaczul6W3o0bdpUWoO7cUmy0zf+svKUv2XlKX/Lkb8dlIkpqRb1Yfr06RItj0gC8tdff7WofWuRkpIiN27cKIcMGSKdnJwy/A4A6eLiIl9++WW5b98+OWjQoEznHRwc5JIlS6z9Y1iduLg46ePjY/y9zJgxQ1d7SSlpsu/3e4zvmRlbzupqT5E9wFGZi8/Y3JZe6SGlvGZ4XgnYIqV8SggRKKXMNw0NrFnL60Z0AoN/3s+tmEQA+tQvz/fDG2NvwQSub775hnfeecf4fM6cObz44osWs29Jzp49y4IFC1i8eLExn+IBdnZ2dO/enbFjxzJgwABjoqGUEj8/P5YuXUpkZCR169Zl4sSJhXpm8jDHjx+nRYsWxtIs+/fvp1mzZrrZuxIRT5/v93A/OQ0hYOmEFrSu7qqbPUXWmHOG0hu4BvgDO4EQoA9QAngzN6plKw9rzVAeEHznnmz0uZ/xruuDdSdkenq6RX349ttvM9yB//TTTxa1rydRUVFyzpw5skWLFplmGoCsVauWnDZtmgwNDbW2q/mW2bNnG3+f1atXlzExMbraW3P0uvH90vy/W+XduCRd7SlMQy5nKLn6IAaKAA2BRkDR3Fxjiw9rC4qUUgZei5J1Pt5sfJPM9DtncR9mzZqV4YP2hx9+sLgP5iI1NVX6+fnJ4cOHy6JFi2YSEWdnZ/nCCy/I/fv3W1y8CyLp6emyT58+xt/vyJEjdbf3xvIA4/vl+QVH1N/RCuRWUHLbD6U1UAUw7iRLKRfleKGNYSvl6/cEhzN+wRFS0rTf/af9nmJcm6oW9eHHH39k0qRJxuezZs3ijTfesKgPeSE4OJgFCxawaNGiTJ0NhRB06dKFsWPHMnDgQIoXL24lLwsm4eHhNGzY0LiUuGjRIkaN0i82JzYxhT7f7+WaIdpr6oC6jG5VRTd7isyYs5bXYqA6cBxIMxyWMh8VhXyArQgKwN8nbjJpeSAPfv2zhzViQCMvi/rw888/88orrxifz5w5k7feesuiPjwO9+7dY9WqVSxYsIB9+/ZlOl+9enXGjh3L6NGjC0S2ui2zfft2unXrhpSSkiVLEhAQQM2aNXWzF3gtiiFzDpCaLnFysOOvV9tQp3zBSya1Vcy5h3IWg/CY6wGURateHGz4WsbEmEbAAeA0cAJ49qFzVYFDhutXAk65sWsLS14Ps+jAVeNUvvr7G+WOc3cs7sPcuXMzLBHpHb3zuKSlpcnt27fLkSNHymLFimVa0ipZsqQcP3683L17t1oKsTDvv/++8e/QtGlTmZSk7/7G//wvGt8vXWfulPeTLBspWZjBXHsowGqgfG5eLLcPYAbwnuH794CvTIypBdQ0fF8BuAW4GJ6vAoYZvp8DvJwbu7YmKFJKOXvbBeObxOejTfLo1bsW9+HXX3+VQgjjh8O0adMs7sOjXLp0SX7yySeycuXKJjfYO3bsKBcuXChjY2Ot7WqhJTk5WbZs2dL4N3n77bd1tZeWli5H/HrA+H55f90JXe0p/sWcguIPRAF+wPoHj9y8eDavef6BSAHlgfO5uCYIqAkIIAJwMBxvBfjlxq4tCkp6err89K9TxjdJ/U+3yHO37lncj99//z2DqHzxxRcW9yE2NlbOnz9fdujQwaSIVK5cWX766afy8uXLFvdNYZrLly/LUqVKGf9Gmzdv1tXe7ZgE2XjqP8b3y+aTN3W1p9DIraDkZg+lg6njUspd2V6Y/WtGSyldHnoeJaUsk8345sBCoC7actlBKWUNwzlvYLOUsl4W104EJgJUqlSpaUhIyJO6rRvp6ZL/rDrOX8e1ukUepYqw5qXWeJe17GbywoULGTdu3AMB5/PPP+eTTz7R1aaUkj179jB//nxWr15NfHx8hvPFixdn8ODBjB07lg4dOuiana14MlauXMmwYcMAcHd3JygoCE/P7Iul5oUd5+4wfoG2F1qqqAOb32yPl0sx3ewpzLiH8qQPYBtwysRjABD9yNiobF6nPNqMpqXhuRtw8aHz3sDJ3PhkizOUBySlpMnR8w4Z77w6zNghw2MTLe7H4sWLpZ2dnfGO89NPP9Vlb+Lq1aty6tSpslq1aiZnI23btpXz5s2T9+5ZframeHzGjx9v/Nt1795dpqWl6Wrv4Vn9kJ/3y5RUfe0Vdsjrkhew1/A1Frj30CMWuJebF8/mtXO15AWUAgKAIQ8dK1BLXg8Tn5QiB/601/hG6T17t7yXkGxxP5YuXZpBVD766COziEp8fLxcvHix7Ny5s0kR8fb2lh999JEMDg42w0+hsCSPlmb5+uuvdbWXkJwqe87abXyvfLf1vK72Cjt5FhQ9H8DXZNyUn2FijBOwHRPZ+GiBAg9vyr+SG7u2LihSShkVnyS7fbvT+EZ5du5+mZBs+WiWFStWSHt7e+MHxPvvv/9EopKeni737dsnJ0yYIJ2dnTOJSNGiReWIESPk1q1bZWqqitrJzwQEBBjrpTk4OMjDhw/rai/4Tqys/ZGWJFz1vb/l4SuWD2gpLJhNUIDnTRybnpsXz+Y1yxnEItjwtazhuC/wm+H7kUAKWv7Lg0cjw7lqwGHgokFciuTGbn4QFCmlvBWdIFtP224UlRcWHrHKlH7VqlUZROXdd9+V6enpMj09XSYkJGQrMNevX5dffvmlrFWrlsnZSKtWreTcuXNldHS0BX8ihd48XIWhevXqui9ZLj8UYnyftPpym4yOt/yMvjCQW0HJzab8ZmCJlHKp4fn/0Mqv5Lv657aU2JgTVyLiGfzzfu4a+qc86+vN9GfqW7wfx7p163j22WdJTdXasbZo0YIrV64QFhZGuXLlGD9+PB9++CGlS5cmISGBv/76iwULFrB161bS09MzvFaFChUYPXo0Y8aMoXbt2hb9ORQa95NTOXwlkuTUdBp5u+BeyrwtjaSU9OvXj40bNwIwatQoFi3Sr6iGlJLXlgWy8aSWte9dphgO9gJHezs6+bgzrk1VPEurtk15xZyZ8sXQQoV/B3oBkVLKN83ipYXJT4ICcOpGDMN+OUicobf2Sx2q814vy38Q//nnnwwdOpSUlBST52vVqkW7du1Ys2YNMTExGc45OTnx9NNPM27cOLp164a9vb0lXFY8gpSSeXuvMHt7MLGJ2v+TvZ1gUGMvpg6oRzEn8/1dwsPDadCggbERl96lWWISUugycycRcZmb15Ut7sTyiS3x8XTWzX5hIM+CIoR4uJ+tM/AnsA/4BEBKGWkGPy1KfhMUgP2XIhj7+xGS07S7/Q971+GF9tUs7seLL76YoXNfTjRr1oyxY8cybNgwypbNujWywjLM33eFzzecMXmuR10P5o7KOSL0cXi0NEtgYKBubQASU9Jo/t9t3DMI5aM0qFiav15tU+i6bZqT3ApKdkH9x4Cjhq/+gAta2foHxxUWoHV1V74f3pgHbVP+u+ksq49et7gfQUFBOY7x8PBg8uTJnDx5ksOHD/PKK68oMbEBklLT+GHHxSzP+52+w6kbMVmefxK6dOnClClTAIiLi2P48OEkJ+vT/nrb2TtZignAidAYzt6K1cW2IiNZCoqUsqqUstojXx88LH+LXIjpWc+TaYPqG5+/t+4kW8/csagPERER2Z53cXEhNDSUr7/+mnr1TOaYKqzEkSuRRMZn/2G+83yY2e1OnTqV5s2bA3D06FFjG2pzExqVkOOYa3fjcxyjyDs5ph0LIRyFEK8LIdYYHq8JIRwt4ZziX55tVokpPbX9k7R0yavLAjh0+a7F7NepUyfb8w0bNsTBwSHbMQrLcSUint/3XmHUvEOMnX8kx/Fp6TkOeWwcHR1Zvnw5zs7a/sU333yDn5+f2e2Uz8Wm++d/n2HJwRASU9JyHKt4cnKzKf8b4IhW+gRgFJAmpZygs29mJz/uoTyMlJIvN53l1z1XAHAu4sCKF1tSt0Jp3W1v3ryZ3r17Z3l+1apVDBkyRHc/FKZJTEnj4OW77Dwfjv/5MELu3n+s61dObEmLauV08W358uWMGDEC0EqznDhxAg8PD7O9fkJyGi2nbScmwXTQyMOUK+HEmNZVGNWyMmVKOJnNh4KOOaO8gqSUDXM6lh/I74ICmqhMXn2CtQFaUynXkkVY+3IrKpcrobvtjz76iP/+97+Zjr/++uvMmjVLbXpamOuR9/E/H8bO8+HsvxRBYorpaUZ9r9IUdbDjSEiUyfPuzkU49EEXXf9+48ePZ/78+QD06NGDTZs2mbUum9/p27y6NIDU9IyfZ8Wd7PGtXJZ9lyJIe+hcMUd7hvpWZEK7ahavmZcfMaegPCh9csnwvBqwRkrZxCyeWpCCICgAqWnpvLTkGNvOauve3mWLsfal1mbPKTDFgQMH+O233wgJCcHLy4vnn3+e9u3b625XoW2uH70ahf+5MPzPh3Ep3PS+gHNRB9rXdKOjjxsdfNxwdy5KWrpk6obTLD4YQrqJt/xvo33p+pT5Zg2PEhcXR9OmTblw4QKgLX+9/fbbZrVx6kYMv+y+zKErd3Gws6NzbXcmtKtK5XIluB55n3l7r7DyyHUSHlr2shPQu355XmxfnfoV9Z/p51fMKShdgPnAZbQ6WpWBcVJKf3M4akkKiqCAtsQxet5hDl/VordrezqzcmIrShdX21sFiZvRCcZlrH0XI7ifbHoPoLanM51qu9PJx53GlVxwtDd9938jOoEd58JITk0nJS2d6ZvPAVCmuCOb32ivaxJgYGAgLVu2JDk5GUdHR/bv34+vr3nDlXMiKj6ZJQdDWHjgaqa8ldbVyzGxfTU61HJTs+1HMJugGF6sCOCDJijnpJRJeXfR8hQkQQG4l5jCs3MPcvbWPQB8K5dh8fMtzJqkprAsKWnpHAuJwv98GLvOh3Putulw1xJO9rSt6UpHH3c6+rhRvvSTlW//6M+TLDl4DYCW1cqydEJL7O30+zCdNWsW//nPfwCoUaMGAQEBxk17S5KYksa6gBv8uucyVyIyzvRqezozsX01+jWskKUwFzbMOUOxR8s/qQIYw3iklN/m0UeLU9AEBSAsNpEhcw4YN2E713Zn7qim6o2Qjwi7l8jOC+HsPB/GngsRxCaZzqmo4V6STj5udPJxx7dKWZwc8v43TkxJY8CP+zh/RxOuyd1r8Vpn/XrDSynp27cvmzZtAmD06NEsXLgwh6v0Iy1dsvXMHebuvkTgtegM58qXLsr4NlUZ1twb56KFe+ZvTkHZBCQCJwHjrp+U8vO8OmlpCqKgAFy7e59n5uwnPFabOA5q7MU3Qxpip+OdpuLJSUuXHL8ehf85bSnr9M17JscVdbSjTXVXOtZ2p2MtN902jy/ciaXfD3tJSk3H3k6w6sWWNK2sX0JqWFgYDRs2NJZmWbx4MSNHjtTNXm6QUnI0JIq5uy6z7WzGHC/nog4816Iy49pUwcMC+5S2iDkF5YSUsoHZPLMiBVVQAM7eusfQuQeMdZqeb1uVj/rUUWvBNsLduCR2XQhn5/lwdgeHE33fdIhrVdcSdKjlRqfa7rSoWpaijpZZvlx6KIQP/zgFgJdLMTa90Y7SxfS7K9+2bRvdu3c3lmY5fvw41atX183e43AxLJZfd1/hj8AbxpJHAI72goGNvZjYvho13AtXbTBzCspXwHYp5T/mcs5aFGRBATh8JZJR8w6RlKq9Cd7p4cOrnfSpn6TInvR0yckbMfifD8P/fDgnQqMx9VZzcrCjZbVydPJxo6OPO1Vd9Q//NoWUkleWBrD5lDZr6F3fk59GNNH1huS9997jq6++ArTab3v37sXJyXZyQ8LuJTJ//1WWHAwx3qg9oEttd17sUJ1mVcoUips2cwrKQGAJWlZ9CtrGvJRSljKHo5akoAsKwPazd5i4+Jgx5n7aoPoMb17Jyl4VDqLvJ7M7OIKd58LYdSHc2HrgUbxcitG5tjudarvRslo5ijvZRoWBmPsp9P5+DzeitVImev/vpKSk0LZtWw4fPgzAO++8w4wZM3Sz96TEJaWy4vA1ft97hZsxiRnONfJ24aUO1ej2lKeuwQzWxpyCchl4Gq1ve84hYTZMYRAUgHUBoby1SivmaCfgpxFN6FW/vJW9yh9IKbkUHkf0/RSquJbAtWSRbMeevnmPnYbkwoBrUSZzPBztBc2rlqVjLU1EqruVtNm72qNXIxk69wDpUtvD2fBaW2p66Le8c+nSJRo3bkxsrBYU4OfnR/fu3XWzlxdS0tL5+8RN5u66nCn6rkq54kxoV43BTStabJnSkphTUPyAXlJKHar9WJbCIigAv+25zBcbzwLgZG/HgnHNaF3D1cpe2TZHrkby2frTxk1yBztBnwblmdq/njG/515iCvuCI4wZ6mGxpiPoPUsVpVNtbRmrTQ1XShaxjVlIbvh+ezDfbtUSEGt7OvPnq210/ZBctmwZzz33HKBVrA4KCjJraRZzI6VkT3AEc3dfYt/FjPX0CmppF3MKygK0lrubAeO7R4UN2z4ztpzjfzsvAVrewvKJLWlQ0cXKXtkmJ0KjGTLngHH/6WFqupfk6cZe7L4QzrGQqEzlPUBrVtW0chk6+WizEB8PZ5udheREWrpkxK8HOXRFS5od1bIy//e0vhWkx44dawwf7tmzJxs3bjRraRa9eJCdv/HkrQJd2sWcgvKpqeMqbNj2kVLywR8nWX5Y659StoQTq19qRXW3klb2zPYYN/8w/ufDH+sa15JFjJvpbWu66hoVZWluxSTQa/YeYzTa3FFN6VHXUzd7cXFxNGnShODgYABmzpzJW2+9pZs9c1PQS7uYNVP+oRe1A0pKKU0Hzts4hU1QwFDqfmkAW05r0TteLsVY83KrJ86sLogkpqRR5+Mt5PROsBPaJqw2C3HnqfKlCnSuz9Yzd3hhkfZ+KV3Mkc1vtKOCi37/NwEBAbRs2ZKUlBQcHR05cOAATZs21c2eHhTU0i7m6Nj44IWWCSFKCSFKAGeA80KId/LoXFkhxFYhRLDhaxkTYxoJIQ4IIU4LIU4IIZ596NwCIcQVIcRxw6NRXvwpyNjbCWYNa0Tr6lpp8hvRCYyed5ioHBouFRYuhsXxzT/ncxSTuhVKceyjbqx7pQ2TutSknlfpAi0mAN2e8mBMq8qA1rf9zRXHMyzrmJsmTZoYw4hTUlIYNmyYcbM+v1CmhBOTutRk75TOfDmwfoYw8P2X7jJ2/hF6zd7DuoBQUvRoQmNlcrPkdVxK2UgI8RzQFJgCHMtLsqMQYgYQKaWcLoR4DygjpZzyyJhaaOHJwUKICmith+tIKaMN+zp/SynXPI7dwjhDeUBcUirDfznISUOr10beLiyd0IIS+Wiz2FxExiezIegm6wJCCQrNXevbwprTk5iSxtM/7TNGNb3ZtSZvdq2lm7309HT69u3L5s2bARgzZgwLFizQzZ7ePCjt8svuSwTk49Iu5txDOQ00ApYBP0opd+W1H4oQ4jzQUUp5SwhRHtgppfTJ4ZogYLBBYBagBOWxuRuXxJA5B7hsKIbXrqYr88Y0M0tNKFsnKTWNHWfDWBtwg53nwzJtrBdxsDO5IQ9a6Q3/yR2zDSEuyFwMi6XfD/tISEnDTsCKia1oXlXf0iwNGjTgzh2tBMqSJUuMUWD5maNXI5mTT0u7mFNQXkeblQShFYmsBCyRUrbLg3PRUkqXh55HSSkzLXs9dL45WsfIulLKdIOgtEKLOtsOvJebCsiFXVAAQqPuM/jnA9y+pyVo9WtYgdnPNiqQyzdSSgKuRbE24AYbT9zK1NHPTkDbmm4808SLrnXc+WX3FX7YEZwhl6RcSSd+GdVU19pW+YGVR64xZe1JACqULsqmN9rhUly/sNitW7ca81GcnZ0JDAy0mdIseSU/lnbRZVPe8MICsJdSmi6J+u+4bYCpsJAPgYW5FZQHMxhgjJTy4EPHbgNOwC/AJSnl1CyunwhMBKhUqVLTkJCQ7H/AQsCFO7EMmXPA+AE7ulVlPu9fN99tFGbFtbv3+SPwBusCQ022wq3t6cygJl4MaOSV6Y4wNOo+f5+4RfT9FGp5lKR3/fIFMlHtcZFS8tryQDaeuAVAj7oezBnZVNf/mSlTphgz522xNEteCbuXyIL9V1n8GKVdbsUkcP52LKWKOdKwoovFsvN1ExRzkNslLyFEKTQxmSalXJ3Fa3UEJksp++ZkV81Q/iXgWhTP/XrIGOKo99q43sQkpLDp5C3WBYRy5GrmVreuJYvwdKMKDGpSkacq5LuqQTZBTEIKvWf/W5rl/56ux6iWlXWzl5ycTNu2bTly5AgA7777rnHTviCRm9IuzaqU5eO/TrHl1G3jDNq7bDE+6VuXbjp22nyArQvK18Ddhzbly0op331kjBNaMuUGKeWsR86VN4iRAL4DEqWU7+VkVwlKRnZdCGfCwiOkpGn/A1MH1GV0qyrWdeoxSElLZ/eFcNYF3mDrmTskP7IHUsTBjh51PRnYxIt2NVxxUD1i8syxkCiGzj1AWrrEycGO9a+1obanfgKdn0qz5JXsSrs4Odhl+v8Gbdl20fgWtK2pbxUMWxeUcsAqtP2Ya2g96yOFEL7AS1LKCUKIkWith08/dOlYKeVxIcQOwA2tUOVxwzVxOdlVgpKZ9UE3eWNFIFKCEDDr2UYMaORlbbey5EH9rLUBoaw/ftNkAcYWVcvyTJOK9KzvSSkbj57Jj/zkf5Gv/c4DWhWB9a+11bVL6NKlS439Ujw8PDhx4gTu7u662bM22ZV2MUWTSi6se6WNrj6Zc1O+OPA2UElK+YIQoibgI6X82zyuWg4lKKZZdOAqn/yl6baDnWDe2GZ0qOVmXace4VZMAn8G3uSPwFAu3Ml871DNtYRxX6QglLqwZdLSJaPmHWL/Je3DbkSLSnw5sL6uNseMGcOiRYuA/FWaJa+cuhHDy0uOcT0qIdtxgR9307V2mNkSG9FmCUloUVUAocAXefBNYWOMblWFN7pobV9T0yUvLT5GwLXM+xCWJj4plXUBoYz87RCtp+/gqy3nMoiJS3FHRrWszB+vtGb72x14rXNNJSYWwN5O8N2zjShr+ABbdugam0/e0tXmjz/+SI0aWh7Qli1bmD17tq72bIV6XqVpUjnLAFgjtpIkmRtBqS6lnIHWCwUpZQLaUpOiAPFm15qMNmRFJ6SkMX7BEYLvWD5LOS1dsjc4grdWHqfZf7fx1qog9l6MMDancrQX9KjrwdxRTTn8QVf+7+l6NK5UOJoc2RIepYryzZB/c5unrD1BaFTmiDpz4ezszIoVK3B01JYwp0yZQkBAgG72bAnfKtmHrFcpVxw3Z9vIkcqNoCQLIYqBVp1CCFGdh6oOKwoGQgg+61eXvg20vinR91MYNe+wrh8SD3PhTizTN5+jzfQdjJx3iHWBN7if/G+RvcaVXPi/p+tx+IOuzB3lS4+6noUiIdOW6Vzbg3FtqgBwLzGVN1ccJ1XHO+WmTZsyffp04N/SLHFxOW6d5nsGNvbCM5uEx5c6VLeZG6rc7KF0Az4CngL+AdqgbY7v1N07M6P2UHImOTWd5xceYU9wBKDtTax+qRXldMgSj4hLYv3xm6wLDOXUjcz1Rr1cijGoiRcDG3tRTVVItkmSUtMY+NN+ztzS/n6vd67BW92zLXqRJ9LT0+nTpw9btmwBtLL38+fP182erXDhTiwvLDqaIa/KwU4wqXNNXu9SQ3dBMWuUlyEqqyXaUtdBKWVE3l20PEpQcsf95FRG/HqI49e12kP1vUqzfGJLszSJSkxJY9vZO6wLuMGuC+GZig2WLOJAn/rlGdjEi+ZVyhbIDP6CxqXwOPr9sJf7yWkIAcsmtKSVoRipHty5c4eGDRsaS7MsXbqUESNG6GbPVkhNS2fXhXDO3rqHc1FHetXzxN1CpVryLChCiCbZXSilzHcLmEpQck9UfDJD5h7gYpi2pNC6ejl+H9vsibLGpZQcuRrFH4Gh/H3iVqasYHs7QbuargxqUpFudTx0DUFV6MPqo9d5Z80JADxKFWHzG+2Nm/Z68M8//9CjRw9A2185fvw41apV081eYcccguKfzXVSStn5SZ2zFkpQHo9bMQkM/vmAMTO6Z11PfnquSa7LPVyNiGdd4A3+CAzlemTmsMenypdiUBMv+jeqgLuzbRbFU+QOKSVvrDjO+qCbAHSt486vo311XYp59913+frrrwFo3rw5e/fuNW7aK8yLTSc2WgslKI/PpfA4hsw5QKQhgbBDLTeKO9lxIzoRz1JFGeRt+z4AACAASURBVOrrTZc67sYPjuj7yfx9QiuB8mi5bgB35yIMbOzFwCZeumZYKyxPbGIKvb/fY7x5+Lx/Xca0rqKbveTkZNq0acOD9/SUKVOMm/YK82LuPZR6aJvyxttIKeWiPHloBZSgPBknQqMZ/stB4h+KunqYZxp70aOeJ38E3mD72bAMFVQBijra0bOuJ4OaVKRNDVeLFbRTWJ7j16MZ/PN+Ug2lWf58pY2utdMuXrxI48aNjdFea9asoWTJkjg4ONCqVSuKF1d5SebA3D3lO6IJyiagF7BXSjnYDH5aFCUoT86sbReYtS041+OFgFbVyjGwsRe96pc3y4a+In8wZ9clpm8+B0B1txJsmNSW4k76/f2XLFnCqFGjMh13cXHhgw8+YPLkyTYTVptfMWem/GCgC3BbSjkOaAjYRhaNwmKYWr4yRXW3ErzTw4e9Uzqz7IWWDPH1VmJSyJjYrhrtDMUKL4XH8/n6M7raGzlypDGL/mGio6Mz7LMo9Cc3gpIgpUwHUg3l5MMAFU5RyAi5G5/t+aKO9qx/rQ3b3urAq51q4OVSzEKeKWwNOzvBzKENKWeI8lp59DobDJv1enDx4kUuXryY5fkvv/yS+Pjs/38V5iE3gnJUCOEC/IrW1z0AOKyrVwqbwy2HxMYq5YrToKKLWlpQAODuXJSZQ//tEv7BupNcj9Sn6sKD/vNZERMTw759+3SxrchIjoIipXxFShktpZwDdEPrnDhOf9cUtsSgJhWzPT+4afbnFYWPjj7uvNCuKgCxSalMWh6oSxHDlJSUHMcEBQWZ3a4iM7kqhiSEaCCE6A80AWoIIQbp65bC1hjctKJxXfxRfCuXYaSOnfsU+Zd3etSmvldpQIsA+27rBbPbaNeuXY5j3n33XTp06ICfnx+FKVXC0uQoKEKI34HfgWeAfoZHju12FQULJwc7fhvjy/u9alPVtQRO9nZ4ly3G5O61WPx8C9V3XWESJwc7vh/emBKG6gc/77rEvovmrdzUrFkzunTpkuO43bt307NnT5o2bcqaNWtISzMdBq94cnITNnxGSvmUhfzRFRU2rFBYh3UBoby1Slt2cnMuwpY32pm14GhkZCRDhgxhx44dGY4/99xzDB06lK+//pq9e/dmOOfj48OUKVN47rnncHLSr0xMQcCceSjzgJlSSn1j/yyAEhSFwnq8tfI46wJvANDJx415Y5qZtfinlJLDhw+za9cuHB0d6dWrF7Vr1zae37NnD9OmTcu0ie/t7c0777zD888/rxIhs8CcgtIe2ADcRuuDItBqeTXI9kIbRAmKQmE94pJS6fv9Hq4aSrB/3Pcpnm9b1eJ+BAYGMn36dFavXp1hP8XNzY0333yTV155BRcXF4v7ZcuYU1AuAm8BJwFjiIaUMiSvTloaJSgKhXU5GRrDoJ/3kZImcbQX/PFKG+oZNu0tzfnz55kxYwaLFi0iNfXfCtilSpXi1Vdf5c0338Td3d0qvtka5syUvyalXC+lvCKlDHnwMIOPCoWikFG/Ymne7aEtQ6WkSSYtDyQ+KTWHq/TBx8eHefPmcfnyZd544w2KFdOSce/du8e0adOoXLkykyZNIiREfdzlltwIyjkhxDIhxHAhxKAHj7waFkKUFUJsFUIEG76WMTGmshDimBDiuBDitBDipYfONRVCnBRCXBRCfC9URp1CkS94vm1VOtRyA+BKRDyf/HXaqv54e3sza9YsQkJC+PDDDyldWpsxJSYm8uOPP1KjRg3GjRvHuXPnrOpnfiA3glIMbe+kO+YNG34P2C6lrAlsNzx/lFtAayllI6AF8J4QooLh3M/ARKCm4dHTDD4pFAqdsbMTfDOkIa6GKK+1AaH8adistyZubm588cUXhISEMH36dONyV2pqKgsWLOCpp55i8ODBHDt2zMqe2i5W64cihDgPdPz/9u48zorqzP/459s0+75vstMGZRUIi0JYWiUxMq5JRIPOJHFDYEaCiWNmoq/4m4wJMVFxDzHgSFyJwohRBAF1MCwKjYABGhDZlB1lkfX5/VHVze3mdvft7rt187xfr/vqe6vqVp3TdW89t06deo6Z7ZDUElhgZkUORh0OQ7ycYChiA+abWZdw3qhwXbcWt02/huJc+nhv/S5G/ynI4lSneiazxw+iXePaKS7VaUeOHOGZZ55h0qRJZzR7jRgxgnvuuYfBgwefFemG4nkNJVGam9kOgPBv1KtfktpIWglsAX5jZtuB1sDWiMW2htOccxXE4Kym3DokyDN78OgJxj+/nGMn4p+apaxq1qzJHXfcwfr165k2bVqBLshvvfUWQ4YMYdCgQcyePdvvvg8lNKBImitpVZTHFbGuw8y2hF2UOwM3SWpO0HX5jEWLKMMtkpZJWrZr166yVcQ5lxATL/0GPdsEXXRzth7gwTlrU1yiM1WtWpUbb7yR1atXM2PGDPr06ZM/b9GiRVx++eX06tWLF1544ay/+z6hAcXMLjazblEeM4EvwqYuwr87S1jXdmA1MJjgjCQyG+E5QNT82Gb2tJn1NbO+TZs2jUe1nHNxUrVKBpOvuyB/zJyn3t3Iu+vS84dfRkYGV199NUuXLmXOnDkMHTo0f97KlSsZNWoUXbp0YcqUKRw9ejR1BU2hUgUUSa/HcduzgJvC5zcBM6Ns7xxJNcPnDYGLgLVhE9lXkgaEvbtujPZ+51z6a9u4Fv91Vbf81xNeymHXV+l7QJbEJZdcwvz581m0aBEjR47Mn5ebm8vNN99Mx44d+cMf/nDWjcNS2jOUeF6neAC4RNJ6grT4DwBI6itpSrjMecBiSTnAQuB3ZvZxOO92YAqQC2wAih8UwTmXtq7o1Tp/CITdB4/y05dzOHUq/a9LDBw4kFmzZpGTk8OoUaPIyAgOqdu3b2fChAm0a9eO+++/n3379qW4pMlRql5ekp4xsx8lsDwJ5b28nEtfh46eYOTk99m4O/hVf89lXbjlW51SXKrSyc3NZdKkSUydOpVjx47lT69Tpw633347d955Jy1btkxhCcsmIb28KnIwcc6lt9rVM3lk1AVUqxIcln775lpytuxPcalKp3Pnzjz11FNs3LiRCRMmULt20A364MGDTJo0iQ4dOjBmzBg2bdqU4pImRiq7DTvnXAHdWtfn7u8E3XNPnDLGv7Ccr74ueUTGdNO6dWsefPBBNm/ezL333kvDhkEikKNHj/LEE0+QlZXF6NGjWb06tVkC4s0DinMurfzLRe0Z3iW4LW3znsMpT81SHo0bN+a+++5j8+bNTJo0iRYtWgBw8uRJnnvuObp168ZVV13FkiVLUlzS+Ig5oEhKn1tYnXOVliQmXduDZnWD1CyvLt/GjA+3lvCu9Fa3bl0mTpzIpk2bePLJJ+nQ4XTa/tdee43+/ftz8cUX884771TomyRjGQL4QklrgE/C1z0lPZ7wkjnnzlqN61TnoR/0Ii+ryX/OXMXGXQdTW6g4qFGjBrfeeivr1q1j+vTpdO3aNX/evHnzyM7OZuDAgcycOZNTp9Ina0CsYjlD+QMwAtgDYGY5wLcSWSjnnLuwcxPGDA16eR0+dpLxLyzn6InKcSd6ZmYm119/PStXrmTmzJn0798/f97ixYu58sor6dGjB9OnT88fq2XNmjWMGzeOwYMHM3LkSKZPn87x4+l1fSmWAbYWm1l/ScvN7IJwWo6Z9UxKCePIuw07V7EcP3mK7z/1Acs/C3p7/WRQB/7j8vNTXKr4MzMWLFjAr3/9a+bOnVtgXocOHRgyZAjPPvvsGWctQ4YMYfbs2fm9yRIlnt2Gt0i6EDBJ1SRNJGz+cs65RKpaJYNHrruAumFqlinvb2L+2mKzNFVIkhg2bBhvv/12/hlKnk2bNjF16tSoTWALFy7k3nvvTWZRixVLQLkNuIPTGX57ha+dcy7h2jSqxX9f0z3/9cSXctj55dcpLFFi9evXj1dffZVVq1YxevToEtPjT5kypcBNlKlUYkAxs91mdoOZNTezZmb2QzPbk4zCOeccwOU9WnHdN9sAsOfQMSa8VDFSs5RH165defbZZ7nsssuKXe7AgQPs3JkeZ22ZRc2QNJkiUsIDmNn4hJTIOeei+OXI81n66V427DrE+7m7eerdjdw+tGKlZimLyC7G0WRmZubfOJlqxZ2hLAM+LObhnHNJU6taJpNH9aZaZnDYenDOWpZ/VvmTLo4ePbrY+VdffXXCL8rHKmVDAKeC9/JyruKbtuhT7p0V3D3fplFNZo8fTL0aVVNcqsQaO3Ysjz322BnTW7RowQcffED79u0Tuv249fKS1FTS7yS9IemdvEd8iumcc6Vz48B2XHxecwC27D3CL15dVaHvLo/F5MmT+eMf/0j37kHnhHr16nHLLbewZMmShAeT0ojlPpQ5wIvARIIeXzcBu8zs54kvXnz5GYpzlcO+Q8f4zsPv8XnY2+u31/bg+33bpLhUyXHq1Kn8cVeSJZ73oTQ2sz8Bx81sYZjCfkC5S+icc2XUsHY1HrrudGqWe2euJndnxU/NEotkB5PSiKVkeff275D0XUkXUHA8d+ecS7oBHRszblhnAI4cP8m455fz9fHKkZqloiqy23CE/yepPvBTYDJQD7gzoaVyzrkYjM/OYtGGPSzbvI9PdnzJDVP+Tt0aVcnMEEO/0Yyre7emVrVYDnMuHryXl3OuQtu67zAjHnqXQ0fPPDvp2KQ2f7l5AC3q10hBySqPePby6iDp95L+KmlW3iM+xXTOufJpUa8GNapWiTpv4+5D3PVKTpJLdPaK5VzwNeBPwP8CFS9Bv3OuUlu4bhd7Dhady+q99bvZuOsgHZvWSWKpzk6xBJSvzeyReG5UUiOCrsjtgU+B75vZvkLLtAP+ClQBqgKTzezJcN4CoCVwJFz8UjNLj2Q2zrmkiqV3V+5ODyjJEEsvr4cl3StpoKTeeY9ybvduYJ6ZZQHzwteF7QAuNLNeQH/gbkmtIubfYGa9wocHE+fOUg1rVStxmZc/3MqWvYeTUJqzWyxnKN2B0cBwTjd5Wfi6rK4AhobPpwELgAI3SppZ5DlsdWILfs65s8yIri345axVfH286Bb5t9d8wfx/7OSa3ucwZlgn2jVOj9xXlU0sB+mrgI5mNsTMhoWP8gQTgOZmtgMg/Nss2kKS2khaCWwBfmNm2yNm/1nSCkn/qZIGDHDOVVr1a1XlF9+NPopjhqBBrSDP14lTxovLtjD8wYX89KWcSjFGfbqJJfXKi8C40jYrSZoLtIgy6xfANDNrELHsPjMrMv9y2NT1GjDSzL6Q1NrMtkmqC8wAnjOzZ4t47y3ALQBt27bts3nz5tJUwzlXQcxZ/TmPL9jAii37keBbWU0Zn92Zrq3q8/KyLTy+YAM7DpwemCtD8E89WzF2eGc6N6ubwpKnv1i7DccSUBYAPYClwNG86Wb2T+Uo3FpgqJntkNQSWGBm3yjhPX8GZpvZK4Wm/zPQ18zGlrRdvw/Fucrv6ImTZEhUrZJxxvQZH27jsfm5bNt/JH+6FAzgNW54Z85t7oElmngGlCHRppvZwjKWDUmTgD1m9oCku4FGZvazQsucEy5zRFJDYDFwDcF49g3MbLekqsDzwNy8HmDF8YDinDt+8hSvfrSNR+fn8lmhC/WXdW/B2GFZnN+qXopKl57iFlDClbUDssxsrqRaQBUz+6ochWsMvAS0BT4DvmdmeyX1BW4zs59IugR4kKADgIBHzexpSbWBdwm6ElcB5gITzKzEJD4eUJxzeY6fPMXMFdt5bH4um3YfKjDv0vObMz47i26t66eodOklnmcoNxNcg2hkZp0kZQFPmll2fIqaPB5QnHOFnTh5itdX7mDyO+vZsKtgYMnu0oxx2Vn0atOgiHefHeIZUFYA/YDFZnZBOO1jM+sel5ImkQcU51xRTp4y3vg4CCzrvijYA2zIuU0Zn51Fn3bpMXZ7ssUaUGK5D+WomR3L65krKZOgGco55yqNKhliZM9WfLd7S95a/TkPz1vPPz4PWvYXrtvFwnW7GNS5CeOzs+jXoVGKS5ueYgkoCyXdA9QMr2uMIcjr5ZxzlU5GhvhO95aM6NqCuZ98wSPvrGfVti8BeD93N+/n7mZAx0b8a/a5DOjYCL8N7rRYmrwygB8DlxJcHH8LmGIVMO+9N3k550rLzJi/dicPz11PztYDBeb1a9+I8dlZXNS5caUOLOW+hiKprZl9FveSpZAHFOdcWZkZC9ft4uF561n+2f4C83q3bcD47CyGnNu0UgaWeASUj8ysd/h8hpldE+cyJp0HFOdceZkZ/5e7h4fnrWPppwWSpNPznPqMz85ieJdmlSqwxOOifOR/o2P5i+SccxWfJAZlNeGizo35+8a9PDJvPR9s3ANAztYD/HjaMrq2qsf47CwuPb95pQosJSkuoFgRz51z7qwniYGdGjOwU2OWbAoCy/u5uwFYvf1Lbv2fDzmvZT3GD+/MiK4tyMio/IGluCavk8AhgjOVmkBejgIBZmYVLjeBN3k55xLpw817eWReLgvX7Sow/dzmdRg3PIvLurekSgUMLHFNvVJZeEBxziXDii37mTxvPfP+UTBJe6emtRk3PIvLe7Qks0rFGeLJA0oUHlCcc8m0atsBHpm3njlrvigwvUOT2twxrDNX9mpVIQKLB5QoPKA451JhzfYveXT+et74+PMC09s2qsXYYZ25qnfrM9LtpxMPKFF4QHHOpdLaz7/i0fm5vL5yO5GH3tYNanLHsM5c06c11TOrpK6ARfCAEoUHFOdcOsjd+RWPzd/AzBXbOBVxCG5Vvwa3D+3E9/q2oUbV9AksHlCi8IDinEsnm3Yf4tF3cnltxTZORkSW5vWqc9uQTozq17ZAYNlz8Cgbdh2ibo1MurSom7R7XDygROEBxTmXjjbvOcTj8zcw46OtnIgILE3qVOe2IR0Z2bMlv31zHbNytnH8ZDA/q1kdfjnyfAZnNU14+TygROEBxTmXzrbsPcwTCzfw8rIt+YEDIDNDBQJN5PTnftKfAR0bJ7RcsQaU9O1W4JxzZ5k2jWrx66u6s/CuYdw4sB3Vwp5f0YJJ3vTfz1mXzCIWywOKc86lmVYNavKrK7rx7s+G0blZnWKXXfLpXvYfPpakkhXPA4pzzqWpFvVrcH7LkrNcHTtxKgmlKZkHFOecS2O92zYodn6bRjVpUqd6kkpTvJQFFEmNJL0taX34t2Exy9aTtE3SoxHT+kj6WFKupEd0NuWIds6dNa7ucw6Na1crcv7NgzumTSbjVJ6h3A3MM7MsYF74uij3AwsLTXsCuAXICh/fTkQhnXMulerVqMq0H/WjZf0aBaYLuHVIR0YPaJeagkVR3HgoiXYFMDR8Pg1YAPy88EKS+gDNgTeBvuG0lkA9M/sgfP0scCXwt0QX2jnnkq1b6/osuGsoc1Z/wSc7vqRujapc3qMlbRrVSnXRCkhlQGluZjsAzGyHpGaFF5CUATwIjAayI2a1BrZGvN4aTjuDpFsIzmRo27ZtfErunHNJVj2zCiN7tmJkz1apLkqREhpQJM0FWkSZ9YsYVzEGeMPMthS6RBKtwTBqR20zexp4GoIbG2PcrnPOuVJKaEAxs4uLmifpC0ktw7OTlsDOKIsNBAZLGgPUAapJOgg8DJwTsdw5wPY4Ft0551wppfKi/CzgpvD5TcDMwguY2Q1m1tbM2gMTgWfN7O6wqewrSQPC3l03Rnu/c8655EllQHkAuETSeuCS8DWS+kqaEsP7bwemALnABvyCvHPOpZQnh3TOOVcsTw7pnHMuqTygOOeci4uzqslL0i5gcxI32QTYncTtJVtlrl9lrht4/Sq6ZNevnZmVOJLXWRVQkk3SsljaHSuqyly/ylw38PpVdOlaP2/ycs45FxceUJxzzsWFB5TEejrVBUiwyly/ylw38PpVdGlZP7+G4pxzLi78DMU551xceEBxzjkXFx5QCpFUU9JCSVXC129K2i/p9ULLZUv6SNIKSe9L6lzCevuFy66QlCPpqoh535a0NhzO+O6I6dMl7ZV0bSLqJ6mdpA/DMq2WdFvEcmUaYlnSNyWdjCyzpJvCoZ7XS7opYvp8SQclxa37Y+H9F04r9xDSkoZKOhCxD38ZMS8p+y/KZ/NkRHlmRSzXQdLi8P/9oqSix489/Z4ekj4IPwcfS6oRTo/6f5I0SdLnkibGo25F1K+tpDmSPpG0RlL7stRP0g0R/6cVkk5J6pXK+kkaVqhMX0u6soz1qyppWliPTyT9e8S8pB1bADAzf0Q8gDuAf414nQ2MBF4vtNw64Lzw+RhgagnrrQVkhs/z0vVnAlUIklt2BKoBOcD5Ee+bClybiPqF26sePq8DfAq0Cl8vIRg+QASJN78Tw7qrAO8Ab+SVGWgEbAz/NgyfN4x4zwKgb6L2XzjtYeAvwKMR00pVP4LRRV+PMj1p+y/KZ/NgEcu9BFwXPn8SuL2E9WYCK4Ge4evGQJWS/k/AfcDERO278LNxScTns1ZZ6ldoG92BjbF8DhJdv4jpjYC9Za0fcD3wQvi8Vvg9bp/Mz2bew89QznQDEanwzWwe8FWU5QyoFz6vTwnjsZjZYTM7Eb6swekBwfoBuWa20cyOAS8QDI+cKPn1M7NjZnY0nF6d8IxVEUMsW/DJyxtiuSTjgBkUHNtmBPC2me01s33A28C341KT6ArsP50eQnpOxLSy1i+aZO6/AnWLJvyFPRx4JZw0jZLrdimw0sxyAMxsj5mdjPP/KRb59ZN0PsEPsLfDMh00s8NlrF+kUcDz4TZSVr9CrgX+Vo76GVBbUiZQEzgGfEnyjy0eUCKFp5YdzezTGBb/CfCGpK0EQxQ/EMP6+0taDXwM3BYGmNbAlojFihzOuLyi1U9SG0krwzL8xsy2U4ohliPW0xq4iuAXVaSU1U+nh5C+K0qZSlW/0EAFzZV/k9Q1Yl0Jr18Rn80akpZJ+ntecwnB2cX+iB8vsZTnXMAkvaWgGfdn4fSy/p9KLUr9zgX2S/qrpOVhE1QVyla/SD8gDCiktn6RrosoU1nq9wpwCNgBfAb8zsz2ksTvXp5UjimfjpoA+2Nc9k7gMjNbLOku4PcEQaZIZrYY6CrpPGCapL9RiuGM4+CM+pnZFqCHpFbAa5JeKWOZHgJ+Hv6yjZyeyvqVewjpCB8R5DM6KOky4DUgq4zrKoton822ZrZdUkfgHUkfE/wyLW15MoFBwDeBw8A8SR+WcV1lVbh+mcBg4AKCg+SLwD8TDMxXpjJJ6g8cNrNVeZPKuq4yiHpsCc+SugNvlaNM/YCTQCuCZuX3FAy/nsz6AX6GUtgRguaoYklqStDevDic9CJwYawbMbNPCH5RdCP41dAmYnYihzMusn7hmclqgi/xVko/xHJf4AVJnxKcwj8e/mpOZf0GAmPDMv0OuFHSA5Shfmb2pZkdDJ+/AVSV1ITk1e+MfRfuM8xsI8H1hgsIEgY2CJs/Yi3PVmChme02s8ME18B6U7bPQVkVrt9WYHnYXHOCIID3pmz1yxN5JpC3jVTVL8/3gVfN7Hj4uiz1ux5408yOm9lO4P8Ivo/J/O4BHlAKCNv4qyjs4VKMfUB9SeeGry8BPgGQdJWk/y78hrDnRmb4vB3wDYKLZ0uBrHB+NYIPfbRfYeVWuH6SzpFUM3zeELgIWGvFDLEsaayksVHW3cHM2lswXPMrwBgze43gl9elkhqG27iU07/GElo/K8MQ0kXVT1KLiB5A/Qi+O3tI0v6Lsu8aSqoePm9CsO/WhNcC5hMEdYgYXruozybB/ughqVb4GR0SritpQ21H+e4tBRqGP94guK5Q1vrlNX9+j+A6Qt42U1m/PPnXdMLlylK/z4DhCtQGBgD/IInHljweUM40h+D0HwBJ7wEvA9mStkoaEf5iuhmYISmH4BpKXjt9J6I3FQwCciStAF4lOODuDtc1luBL/QnwkpmtTlDdoGD9zgMWh3VYSND2+nE4r6ghlrsQHEhjErbl3k/w4V4K/CqcligF9l8xSlu/a4FV4f/qEYJeOJbk/Vd43y0LyzMfeMDM1oTzfg5MkJRL0Cb/p3B61M9meLD7PcH+WQF8ZGazw9nJHGo7v35mdpLgR8C8sClPwB/D5UpVv9C3gK3h2VyklNQPQEE36DYE371Ipa3fYwS94FYR7MM/m9nKFBxbvNtw4QdBs8H/lOP9zwFN41ieqcS323B56/c6UC2O5VlAfLsNp1v94rb/0vCzeR/x7Vbr9Utu/eJ6bDHzbsNnMLPlwHxF3BhXyvf/0Mx2xaMskqYTND98HY/1QVzqd7kFXRDLTdJ8gj7yx0taNlZpVr+47r80+2xOAn5IcC0wLrx+Jb4/rY8t4MkhnXPOxYmfoTjnnIsLDyjOOefiwgOKczEIu2yapC6pLotz6coDinOxGQW8T9CXP2nKegHXuVTwgOJcCSTVIbhx8MdEBBRJGZIeV5Dy/XVJbyhMB64gLfpCBcMDvBWm2Ci83k4K8nAtlfQrSQfD6UMVpPb/C0HeNyRNkLQqfPxbOK29pFUR65so6b7w+QJJD0laFL6nX8L+Qc6FPKA4V7IrCVJbrAP2SuodTr+aIE14d4I8bgMhGJ8CmEzQx78P8AzwX1HW+zDwsJl9kzNTYvQDfmFm5yvImPwvQH+Cu6BvlnRBDOWubWYXEuQ0eybWyjpXVh5QnCvZKE6n7HghfA3BXc8vm9kpM/uc4I51CNLqdAPeDjMj/AcFc0blGUiQhQGC8VoiLTGzTRHbedXMDlmQT+yvBDnXSvI8gJm9C9ST1CCG9zhXZp5t2LliSGpMkEeqmyQjGLTIFKR4L2qURwGrzWxgOTYdeUNdUds5QcEfhYXzRBW+ycxvOnMJ5WcozhXvWoKkku0sSH7ZBthEcNbwPnBNeC2lOcGojgBrgaaS8pvAdHr8lEh/B64Jnxd3sf9d4MoweWNtgnFn3gO+AJpJahwmiry80Pt+EG5/EHDAzA6UKmezDAAAAONJREFUqubOlZKfoThXvFGcOXjaDIKU4XcQDBG9imBI6MUEB+5j4cX5RyTVJ/iePUQwPECkfwOek/RTYDYQ9YBvZh9JmkowXC3AlDCNB5J+FW53E0GG2Uj7JC0iGFn0R6WptHNl4alXnCsHSXUsGHSrMcEB/6Lwekos760FHDEzk3QdMMrM4jJEq6QFBIkNl8Vjfc7Fws9QnCuf18OL3dWA+2MNJqE+wKPhWBz78bMIV8H5GYpzzrm48Ivyzjnn4sIDinPOubjwgOKccy4uPKA455yLCw8ozjnn4sIDinPOubj4/+MefD2pHTaCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Calculate the smoking rates within age/gender groups\n",
    "da[\"agegrp\"] = pd.cut(da.RIDAGEYR, [18, 30, 40, 50, 60, 70, 80])\n",
    "pr = da.groupby([\"agegrp\", \"RIAGENDRx\"]).agg({\"SMQ020x\": lambda x: np.mean(x==\"Yes\")}).unstack()\n",
    "pr.columns = [\"Female\", \"Male\"]\n",
    "\n",
    "# The number of people for each calculated proportion\n",
    "dn = da.groupby([\"agegrp\", \"RIAGENDRx\"]).agg({\"SMQ020x\": np.size}).unstack()\n",
    "dn.columns = [\"Female\", \"Male\"]\n",
    "\n",
    "# Standard errors for each proportion\n",
    "se = np.sqrt(pr * (1 - pr) / dn)\n",
    "\n",
    "# Standard error for the difference in female/male smoking rates in every age band\n",
    "se_diff = np.sqrt(se.Female**2 + se.Male**2)\n",
    "\n",
    "# Standard errors for the difference in smoking rates between genders, within age bands\n",
    "\n",
    "# The difference in smoking rates between genders\n",
    "pq = pr.Female - pr.Male\n",
    "\n",
    "x = np.arange(pq.size)\n",
    "pp = sns.pointplot(x, pq.values, color='black')\n",
    "sns.pointplot(x, pq - 2*se_diff)\n",
    "sns.pointplot(x, pq + 2*se_diff)\n",
    "pp.set_xticklabels(pq.index)\n",
    "pp.set_xlabel(\"Age group\")\n",
    "pp.set_ylabel(\"Female - male smoking proportion\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot above shows for each age band, the point estimate of the difference in smoking rates between genders (black dot), and the lower and upper end points of the 95% confidence interval (blue points).  Based on this plot, we see that in the United States, smoking is more common in men than in women, not just overall, but also in every one of the age bands.  The difference is largest for older people -- for people older than 60, the smoking rate for males is around 30 percentage points greater than the smoking rate for females, while for people younger than 30, the smoking rate for males is only around 15 percentage points greater than the smoking rate for females.\n",
    "\n",
    "Also note that the 95% confidence bands shown above are much wider than the 95% confidence intervals for the data that were not stratified by age.  Stratifying by age leads to smaller sample sizes, which in turn results in wider confidence intervals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confidence intervals for the mean\n",
    "\n",
    "In this section, we discuss how to construct confidence intervals for the mean.  First note that the proportion discussed above is also a mean -- for example, if the data are 0, 1, 0, then the mean is 1/3, which is also the proportion of 1's in the data.  However the proportion has the special property that the variance is completely determined by the mean. That is why we constructed the standard errors for the sample proportion above using p•(1 - p) as the variance.  In general, the variance of quantitative data will not be a function of the mean, as this is a very special property of binary data.  Therefore, in general we must estimate the variance as a separate step after estimating the mean. \n",
    "\n",
    "To illustrate the construction of confidence intervals for the population mean of a quantitative variable, we will use the body mass inde (BMI) data from NHANES.  To begin, we calculate the mean BMI for all women and for all men in the NHANES sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>BMXBMI</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>RIAGENDRx</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Female</th>\n",
       "      <td>29.939946</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Male</th>\n",
       "      <td>28.778072</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              BMXBMI\n",
       "RIAGENDRx           \n",
       "Female     29.939946\n",
       "Male       28.778072"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "da.groupby('RIAGENDRx').agg({'BMXBMI': np.mean})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The numbers in the first column of the table above are estimates of the population mean BMI for all women and for all men in the United States (the population that the NHANES study represents). As with the sample proportions, these numbers are not exactly equal to the mean BMI for all women and men, they are only estimates.  To establish the uncertainty for these estimates, we can use the standard errors for these two estimated means.\n",
    "\n",
    "The standard error for the mean based on an independent and identically distributed sample is equal to the standard deviation of the variable divided by the square root of the sample size.  We next calculate all the relevant values needed to compute the standard error. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr:last-of-type th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th colspan=\"3\" halign=\"left\">BMXBMI</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>size</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>RIAGENDRx</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Female</th>\n",
       "      <td>29.939946</td>\n",
       "      <td>7.753319</td>\n",
       "      <td>2976.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Male</th>\n",
       "      <td>28.778072</td>\n",
       "      <td>6.252568</td>\n",
       "      <td>2759.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              BMXBMI                  \n",
       "                mean       std    size\n",
       "RIAGENDRx                             \n",
       "Female     29.939946  7.753319  2976.0\n",
       "Male       28.778072  6.252568  2759.0"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "da.groupby('RIAGENDRx').agg({'BMXBMI': [np.mean, np.std, np.size]})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now calculate the standard error of the mean BMI for women and for men:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.14211938534506902 0.119045388988243\n"
     ]
    }
   ],
   "source": [
    "sem_female = 7.753 / np.sqrt(2976)\n",
    "sem_male = 6.253 / np.sqrt(2759)\n",
    "print(sem_female, sem_male)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the sample mean BMI for women is expected to be off by around 0.14 relative to the population mean BMI for women, and the sample mean BMI for men is expected to be off by around 0.12 relative to the population mean BMI for men.\n",
    "\n",
    "The standard error of the mean for women is slightly larger for women than for men.  The reason for this is that even though the NHANES sample size for women is slightly larger than that for men, the data for women appears to be more spread out.  The greater standard deviation for the female BMI values leads in turn to less precision when estimating the population mean BMI for females. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As was the case for proportions, the 95% confidence interval for the mean can be calculated by taking the estimate plus and minus 2 (or 1.96) times the standard error.  The 95% confidence interval for female BMI is thus calculated as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "29.661446004723665 30.218553995276338\n"
     ]
    }
   ],
   "source": [
    "lcb_female = 29.94 - 1.96 * 7.753 / np.sqrt(2976)\n",
    "ucb_female = 29.94 + 1.96 * 7.753 / np.sqrt(2976)\n",
    "print(lcb_female, ucb_female)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we show how the one-sample confidence interval can be calculated using Statsmodels.  The numbers differ slightly\n",
    "due to rounding in the calculation above.  The result below is more exact."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(29.659875498090155, 30.22001580625768)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "female_bmi = da.loc[da.RIAGENDRx==\"Female\", \"BMXBMI\"].dropna()\n",
    "sm.stats.DescrStatsW(female_bmi).zconfint_mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Confidence intervals for the difference between two means\n",
    "\n",
    "Now we turn to studying the difference between two means, taking the difference between mean female and male BMI for illustration.  As discussed above, the standard error for the difference of two means taken from independent samples is sqrt(SE1^2 + SE2^2), where SE1 and SE2 are the standard errors for the two means being compared.  Below we see that this gives us a value around 0.19 when comparing the female BMI to the male BMI.  This is substantially larger than either the SEM for estimating the female mean (0.14) or the SEM for estimating the male mean (0.12).  It is expected that the standard error for the difference between two means is greater than the standard errors for estimating a single mean, since the uncertainty of both gender-specific proportions impacts the statistic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.18539073420811059"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sem_diff = np.sqrt(sem_female**2 + sem_male**2)\n",
    "sem_diff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can can now construct a 95% confidence interval for the difference between the female and male mean BMI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.789218531583779, 1.5307814684162213)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bmi_diff = 29.94 - 28.78\n",
    "lcb = bmi_diff - 2*sem_diff\n",
    "ucb = bmi_diff + 2*sem_diff\n",
    "(lcb, ucb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This finding indicates that while the point estimate shows that the women in our sample have around 1.1 unit greater BMI than the men in our sample, the true difference between the mean for all women in the population and for all men in the population could fall between 0.79 and 1.53, and still be consistent with the observed data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Age-stratified confidence intervals** As a final example, we refine the analysis above by considering the difference of mean BMI values between females and males within age bands.  We see below that the overall average difference of 1.1 units results from differences that are very different based on age.  Specifically, the difference between female and male BMI is much smaller than 1.1 for younger people, and much larger than 1.1 for older people.  \n",
    "\n",
    "Since the confidence bands for people under 40 contain 0, the data are consistent with there being no difference between female and male BMI in this age range.  For people older than 40, a hypothetical zero difference between the mean BMI values for females and males is not very consistent with the data.  Informally, we can say that the data strongly suggest that the female mean BMI is greater than the male mean BMI in this age band, with the difference being anywhere from 0.5 to 2 units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd4VGX2wPHvoYQqRUB6EaWoYIGAIqgo6wIqsiio4KprX5qy2HHX/lNsLCjYFlDRFbuISBFFRVERECUgsBQLKL0jJCHk/P5470wmZSaTZFom5/M882RumZlzk8k9975VVBVjjDEGoFy8AzDGGJM4LCkYY4zxs6RgjDHGz5KCMcYYP0sKxhhj/CwpGGOM8bOkYIwxxs+SgjHGGD9LCsYYY/wqxDuAoqpbt662aNEi3mEYY0ypsmTJku2qWq+w/UpdUmjRogWLFy+OdxjGGFOqiMgv4exnxUfGGGP8LCkYY4zxs6RgjDHGz5KCMcYYP0sKxhhj/CwpGGOM8bOkYIwxxs+SgjEmYWRkwAsvwNlnw8knw5VXwsKF8Y6qbCl1ndeMMclp/37485/h669z1v3wA7zyCjz9NAwbFr/YyhK7UzDGJIR7782dEALddBOsWhXbeMoqSwrGmLjLyoJJk4JvV4WJE2MXT1lmxUfGmLg4fBiWLYP58+Hjj2HPntD7r1sXm7jKOksKxpiYyMiARYtcEvjiC/jqK9i7N/zXH3lk9GIzOSwpGGOiYt8+d+L/4guXCL791iWGYKpXd5XNwcya5ZJKp06Rj9XkiFpSEJGmwBSgAZANvKCq4/Ls0x14H/jJW/Wuqj4QrZiMMdGzbZtLAL7H0qWQnR18/3bt4Iwz4Mwz3c/y5aFrV1i/vuD9N22Cbt1gzBgYMgREonMcZV007xSygFtU9TsROQJYIiJzVfXHPPt9oaoXRDEOY0wU/PJLzl3AF1+Ebh1UoQJ07OhO/mec4U7+derk3++bb+Dxx2HqVNi1C048Ea6+GmbMgOnTITPTNU2dPx/+8x+oUSN6x1dWiarG5oNE3gfGq+rcgHXdgVuLkhRSU1PVJtkxJrZUYeXKnLuA+fNhw4bg+1epAqedlnMXcNppUK1ayT7/ySfhzjtdBTVAq1bw9tsucZjCicgSVU0tdL9YJAURaQHMB9qp6t6A9d2Bd4CNwO+4BLGigNffANwA0KxZs46//BLWBELGmGLKynLFP74k8OWXsH178P1r1XJFO74k0KEDpKREPq4FC+DSS+G339xy5cowYYK7m7DipNASJimISHXgc+D/VPXdPNtqANmqul9EzgPGqWqrUO9ndwrGRN7Bg64i2FcU9PXXoSt9GzXKKQo680w44QQoF6NeT9u2wV//Ch99lLPuqqtccijJ3UiyS4ikICIVgRnAHFUdE8b+PwOpqhr0msSSgjElt2ePu+r2FQUtWgSHDgXfv1WrnCRwxhnQsmV8r8wPH4aHH3a9oH2nsBNOcMVJbdvGL65EFm5SiGbrIwEmASuDJQQRaQBsUVUVkc64HtY7ohWTKVsyM6FixeQrVkhPhw8+gF9/hWbNoE8fV4wSyubNuVsG/fBDzsk0LxFXTu8rCurWDRo2jPxxlET58vCvf8Hpp8OgQbB1K6xYAampbkC9QYPiHWHpFc3WR12BK4A0EfneWzcKaAagqs8B/YHBIpIFHAQu01jVfJuklJ0NzzwD48fD6tWuOOGyy9wVZdOm8Y6u5GbNciOHBpbv160LU6ZA795uWRV++il3y6A1a4K/Z8WKru2/ryjo9NNdHUFp0KOHq/sYONAd6x9/wOWXu2P+978LT5Ymv5i1PooUKz4ywajCtdfCiy/m39awoSsnb9489nFFSlqauxLOzMy/rUIFuP12NxTEF1/A778Hf59q1dyJ31cUdOqprrVQaZaVBffcA488krOuQwd46y1X1GUSpE4hGiwpmGDmz4ezzgq+feBAlzCys8N7qIa/b3EfRfmMiRNdD+GiqlMnd33AKae4JJKMPvwQrrjC9XEAqFnT/c379YtvXInAkoIpc264wXVoKuuaNs3dU7ht29i1DEoEv/zimq0GTs7zj3/A6NHRaSZbWsS9otmYWNu2Ld4RxNdRR7lmpaW5iCwSmjd3d4133AFjx7p1//636y39xhvJUbcUTZYUTNI45pjQ21NS4IILXMuVcuUKfogE31bcR6Tec/JkVzwSzMCBlhB8UlJcIujWDa65xo3G+vXXrujslVdyKuVNfpYUTFLYvBnmzg29z+23w4MPxiaeaOjQwdUp+MrLA9WuDSNGxD6mRHfxxXDSSTBgAHz/PezYAeedB6NGwf33J2/dSkmUoZJGk6zS0lwLmmXLgu/TsyfcfXfsYoqGFi3gk0/cSS7QSSe59S1axCOqxHfssS6Z3nBDzrqHH4Zzz3Ujr5rcrKLZlGqzZ8Mll7ix+wHat3flyO+/75JF7dquWKVfP1dslAxU4bvvcjqvdeiQfB30ouXVV+HGG+HAAbdcv74bkfXss+MbVyxY6yOT9J591g2j7Buzv3dveP11G07ZhPbjj6446UdvEP9y5VxR0qhRyd1KK9ykkMS/ApOsDh+GkSPdRCu+hDB0qBtv3xKCKczxx7tWWldc4Zazs92QGeedF3ok2LLCkoIpVfbvh4suci1LwBWbjB0LTz9tlYYmfNWqwcsvu34tlSq5dXPmuNZJxekgmEwsKZhS47ffXIes6dPdcrVqru7g5putTN0UnQhcd53rv3DssW7dxo2uV/yYMcEHDEx2lhRMqfD9966F0dKlbrlRIzfGT58+8Y3LlH4nnwyLF0P//m45KwtuucU1Tti9O76xxUNYSUFEqohIm2gHY0xBZsxwnZB8s22dcoorEz7llPjGZZJHzZrw5pvw1FNu1Fhwd6EdOsCSJfGNLdYKTQoi0gf4HpjtLZ8sItOjHZgx4P5J+/Z1QyKDuzOYPx8aN45vXCb5iMDw4e4OtFkzt+6nn9yIss8+W3aKk8K5U7gP6AzsBlDV74EW0QvJGHcLP3y4qy/wtTAaMQLeew+qV49vbCa5+Yopzz/fLWdmupZugwbl9IdJZuEkhSxV3RP1SIzx7NsHF17oJsoB13Z8wgTX4ihZOqCZxHbkka5Bw6OP5nznXn/dTUaUlhbf2KItnKSwXEQGAeVFpJWIPA2U8UZbJlo2bHD1B7NmueXq1V2dwpAh8Y3LlD3lyrnxsubNy5mOdPVqdyfx0ktxDS2qwkkKw4ETgAzgNWAPYENvmYhbvBg6d84Zw6hpUze5vI1oaeLpzDNd67cePdzywYNw9dVu9FXfcBnJpNCkoKoHVPVuVe3kPf6pqumxCM6UHdOmuX++zZvdcmqqmyTlxBPjG5cx4OaqmDPHzfXt6xPz4otw2mnu7iGZhNP6aK6I1ApYri0ic6IblikrVOHJJ10v5YMH3bp+/eDzz3Nu2Y1JBOXLw333ueRQr55b55s3+4034hpaRIVTfFRXVf1dOFR1F3BU9EIyZcWhQzB4MNx6a05zv9tug7ffhqpV4xubMcGce65rndStm1vevx8uu8yNv5WREd/YIiGcpJAtIs18CyLSHCgjLXZNtOzZ45r8Pf+8Wy5fHl54AR57LLlHqjTJoXFjVwF9++056555Brp2dX0bSrNw/v3uBr4UkVdE5BVgPnBXdMMyyeznn90/j2+mtBo13LwI118f17CMKZKKFV2T1enT3bwd4Ho/d+jgekOXVuFUNM8GOgBvAG8CHVXV6hRMsSxc6Jr0rVjhllu0cHPn/ulPcQ3LmGLr08dNetSpk1vevRv+8hdXLHroUHxjK45wb9QrATtxzVGPF5EzoxeSSVZvvw3du8PWrW751FPdCJXHHx/XsIwpsRYt3PAYw4fnrHvySfd937gxXlEVTzitjx4FFuCKkW7zHreG8bqmIvKpiKwUkRUicnMB+4iIPCUia0VkmYh0KMYxmASnCqNHu9mu0r3GzAMGwKefuukQjUkGlSq5sbrefBOOOMKt++orN3DjnFJUthLOtCR/AdqoalHr1bOAW1T1OxE5AlgiInNV9ceAfXoDrbzHqcCz3k+TJDIzXQujyZNz1o0aBQ8+aBXKJjkNGAAnneR+LlvmZnPr3Rvuvts1af3pJ5g50/1vnHGG67CZSPOBhJMU1gMVcT2aw6aqm4BN3vN9IrISaAwEJoW+wBR1E0V/IyK1RKSh91pTyu3a5caonzfPLVeo4FoYXX11fOMyJtpat3ZFozfdBBMnurvlhx6CSZNgU56zW/fu8NZbULduXELNJ5xrtQPA9yLyvFfU85SIPFWUDxGRFsApwMI8mxoDGwKWN3rr8r7+BhFZLCKLt23bVpSPNnGybh106ZKTEGrVgo8+soRgyo4qVdx0ny+/nNPvJm9CAPjsM3fxlChDc4dzpzDdexSLiFQH3gFGqOrevJsLeEm+X42qvgC8AJCampogvzoTzIIFrvWFbxL0Y46BDz+ENjZNkymDrrzSTffZrVvwE//nn7s7iy5dYhtbQQpNCqr6sohUAZqpapFG+RCRiriE8F9VfbeAXTYCTQOWmwC/F+UzTGKZOtXdDfh6dnbt6sY1SpRbY2Pi4Y8/Cr8T+PzzxEgKUZt5TUQEmASsVNUxQXabDlzptUI6Ddhj9Qmlk6qrPB40KCchDBoEH39sCcGYCmGUyYSzTyyEE8Z9uJnXPgM385qIHB3G67oCVwBpIvK9t24U0Mx7n+eAmcB5wFpc3YWVOJdCGRmuN/Irr+Ssu/fe3CNKGlOWnXaam7hn587g+/hmeou3cJJClqrukdz/3YWW66vqlxRcZxC4jwJDw4jBJKgdO9wIp/Pnu+WUFNfC4q9/jW9cxiSSKlXgn/+EkSML3j5oEBx3XGxjCiacpJBr5jXgJmzmNQOsWeOubtascct16rg5lM84I75xGZOIRoxw840/9JAbCgPcRdR118GYYAXscSBaSO2HiFTF9Wb+s7dqDvBQvCbaSU1N1cWLF8fjo02A+fPdvAe+2+FWrVyHnGOPjW9cxiS6AwdcC71Dh9x4Sb65GaJNRJaoamph+4W8UxCR8sD9qnobLjEYw5Qp7urGN9jXWWfBu++6MlNjTGhVq7o5GRJVyNZHqnoY6BijWEyCU4V77oGrrspJCFdd5TqlWUIwJjmEU6ew1GuC+hbwh29lkH4HJkmlp7uJyqdOzVn30ENuHCNrYWRM8ggnKRwJ7ADOCVingCWFMmLbNtdD+SuveUGlSvDSS24KQmNMcgmnR7P1HSjDVq1yLYzWr3fLdeu6WaVOPz2+cRljoiOcHs2tReQTEVnuLZ8oIv+Mfmgm3ubNc93ufQmhbVs3c5olBGOSVzijpP4HNyfzIQBVXQZYwUGSmzwZevbMaU/do4ebNrNly/jGZYyJrnCSQlVV/TbPuqxoBGPiLzsb7roLrr0Wsry/8rXXwqxZbvhrY0xyC6eiebuIHIM3tIWI9MebPMckl4MH3TC/b7+ds+7RR+G226yFkTFlRThJYShuLoO2IvIb8BNweVSjMlG3caMrDqpY0c38lJEBF14I33r3hJUrw6uvwsUXxzVMU8aoKnPnzmXq1Kns3r2b9u3bc91119GsWbN4h1ZmBB3mQkRuVtVxItJVVReISDWgnKrui22IudkwFyVz4ICbM/nVV11REbgEULlyTv1B/fowfbqbO9aYWDl06BADBw7knXfeybW+UqVKTJ06lX79+sUpsuQQ7jAXoeoUfE1RnwZQ1T/inRBMyV1xhRumwpcQwHVM8yWEE05wLYwsIZhYe/TRR/MlBICMjAwGDhzIr7/+Goeoyp5QSWGliPwMtBGRZQGPNBFZFqP4TAQtW+bGKAqmQgU3ZEXz5rGLyRiA7OxsJkyYEHR7RkYGEydOjGFEZVfQOgVVHSgiDXCjol4Yu5BMJKSnu/4F69bB2rXu8fHHoV+TleU6qzVqFJsYTdm2Y8cOFi1axKJFi1iwYAGbN28Ouf/y5ctjFFnZFjQpiMgnqtpDROao6i+xDMqEZ+/enJN+3p+//Vb4nLAFCSxWMiZS9u3bx3fffedPAosWLeKnn34q0nvs2bMHVUWsKVxUhWp91FBEzgL6iMhU8syipqrfRTUygyps317wSX/dOjcmUVGUKxf6pF+tmtUlmJJLT0/nhx9+yJUAVq1aRWFzt1SrVo0//vgj6PZ58+Zx3nnn8fTTT3OsTdwRNaGSwj3AnUATIO+8QEruAfJMMWVnu6v6YFf8+4pYtZ+SAkcf7Sa7OeaY3D+bN4c+fWDu3IJfe9NNUKNGyY/JlB1ZWVmsWLGCRYsWsXjxYhYtWsSyZcvIygrdv7VRo0Z06tTJ/+jYsSM7duzg9NNPZ8eOHUFfN3v2bNq1a8ddd93FHXfcQeXKlSN9SGVeODOv/UtVH4xRPIWKRZPU7GxX/j5jBmRmuukl+/d3o4MWx6FD8MsvBZ/01693fQSKolq1gk/6xxwDTZpA+fLBX7t7t5s/+cMPc9aVLw9Dh7opAUO91pRt2dnZrF27NtcdwNKlSzl48GDI1x155JF06tSJ1NRUfxJoFKTias2aNYwaNYpp06aRlZVFo0aNGDJkCF27dmX48OG56hWOOeYYJkyYQM+ePSN6nMkq3CapofoptFXVVSLSoaDt8So+inZS2LvXdeL6/PPc61u2hDlzgk83eeBA7ordwJ+//AKHDxctjjp1Cj7pH3ssHHVUyXsYL18OX37pOq/16gWNG5fs/UxyUVU2bNjgP/kvXryYxYsXs2fPnpCvq1atGh07dsx1F3D00UcXuR4gPT2d/fv3c+SRR1KunGskeejQIZ566inuvffeXMVM/fv3Z+zYsTS2L3FIkUgK/1HV60Xk0wI2q6rGpfgo2knhr3+F//634G1t2sArrxR88v/996J/VqNGwa/4bZwhE0vbtm3LdQewaNEitm7dGvI1KSkpnHzyybnuAtq2bUv5KN9ubty4kX/84x+8HTAeS/Xq1bn//vsZPnw4FStWjOrnl1YlTgqJKppJYdMmaNq06Ff1wZQrBy1aFHzSb9nSzdVqTHHs3r2bzZs306BBA2oV8Qpiz549LFmyJFcCKKxjWLly5TjhhBNy3QG0b9+elJSUkhxGicyePZthw4axbt06/7r27dvz7LPP0rVr17jFlagicadwUagXxms6zmgmhTlzXFFKUVSq5E7weU/8vopdu2gxkfT7778zcuRI3nnnHbKysqhQoQIXX3wxY8aMKbCc/uDBg3z//fe5EsDq1asL/ZxWrVrlSgAnn3wy1apVi8YhlUh6ejqjR49m9OjRZARUzl199dU8+uij1KtXL47RJZZIJIUXvadHAacD87zls4HPVDVk0hCRycAFwFZVbVfA9u7A+7gB9gDeVdUHCgs4mklhwQLo1i30Pj16uGkofSf/xo3dHYEx0bZjxw46d+7Met+sRwFatmzJggUL2LRpU64EsHz5cg4XcuvbtGnTXJXAHTt2pHbt2tE6jKhYs2YNw4cPZ86cOf51tWvXZvTo0Vx33XX+eomyLNykgKqGfAAzgIYByw1xJ/DCXncm0AFYHmR7d2BGYe+T99GxY0eNlkOHVJs0UXU9BPI/ypVT/fXXqH28MSHde++9imsOXuCjQoUKIbcDWrduXe3du7fec889+sEHH+imTZvifVgRk52drW+99ZY2btw41zGfeuqp+t1338U7vLgDFmsY59hwhs5uoaqB8ydsAVqHkWzmi0iLMN4/YVSoAE88EXxC+ltvdXUOxkRDRkYGW7duZfPmzf7Hli1b/M9nzpwZ8vV5+wYcccQRpKam5roLaN68edL2CBYR+vfvT8+ePbn//vsZO3Yshw8fZuHChaSmpjJ06FAefPBBatasGe9QE1o4/RTGA62AqbjMexmwVlWHF/rmLinM0ODFR+8AG4HfgVtVdUVh7xmLfgrvvQf/+hes8KJp1MglhBEjbLIZUzRZWVls27Yt18k98Hng8q5du0r0WSkpKdx4443+BNC6desyXWySlpbG4MGDWbBggX9dgwYNePLJJxk4cGDSJsdgItr6SET64YqDAOar6nthBtGC4EmhBpCtqvtF5DxgnKq2CvI+NwA3ADRr1qzjL79EfygmVfj5Z9fxrGVLdxdhEp+qMnnyZMaPH8+yZcuoXbs2AwcOZNSoUTRs2DAin5Gdnc3OnTvDOtFv27at0OEdwlWhQoWQPYUvuuiiAoeeLsuys7OZMmUKt912G9u3b/evP+ecc5gwYQJt27aNY3SxlRBNUkMlhQL2/RlIVdXtofazSXZMKEOGDOHZZ5/Nt75Zs2YsWLCAJk2aFPg6VWXv3r1BT+6Bz7ds2VLoMA7hqlSpEvXr16dBgwb+R+Cy73n9+vVZsmQJ3bt3D/pen332GWeddVZE4ko2O3fuZNSoUbzwwgv+JF2xYkVuu+027r77bqqWgfbhCZ8UvGG5t6iqikhn4G2guRYSkCUFE8yCBQvoFqL5WPfu3bnkkkuCXuGnp6dHJI7y5ctTv379Ak/ueZdr1qxZpGKMCRMmcPPNN+dqUVS+fHnGjRvH0KFDIxJ/Mlu4cCGDBw9m6dKl/nXNmzfn6aefpk+fPnGMLPrinhS8kVW7A3VxldP3AhUBVPU5ERkGDAaygIPASFX9qrD3taRggvn73//O888/H7X3r1evXlhX9XXq1IlqWf7PP//MlClT2LBhA02bNuWqq66iuc2MFLbDhw/z7LPPcvfdd7N3717/+gsvvJBx48bRokWL+AUXRXFPCtFiScEE069fP6ZNm1ak19SqVSvkid63XK9ePRs+Icls3ryZW2+9lf8GjGtTpUoV/vWvf3HLLbfEtbd2NISbFEJNspOGa22UbxNu7KMTSxCfMRFX2Bj7VatW5dVXX6Vx48b+4h0bernsatCgAa+++irXXnstQ4YMYdWqVRw8eJBRo0YxZcoUJkyYwDnnlL0ZAkK1qbkgZlEYU0L79+9n0aJFIfcZPnw4/fr1i1FEprQ4++yz+eGHHxgzZgwPPPAABw8eZNWqVfTo0YNBgwbxxBNPRKzlWmlgxUem1Nu8eTPnn38+330XfDT37t278+GHH5aJViam+H7++WdGjBjB+++/719Xo0YNHnroIQYPHkyFUtw2Pdzio6C1YSKyT0T2Bjz2Bf6MbLjGFM/KlSs57bTT/AmhcePGvPHGG9xwww2cdtpp9O7dmylTpvDRRx9ZQjCFatGiBdOmTWP69On+Cue9e/dy00030blzZxYuXBjfAGMg1IB404AGwLvA66oaemzdGLE7BeMzf/58+vbty+7duwE3bPLMmTOD9kUwpigOHDjAww8/zGOPPcahQ4cAN5TG9ddfzyOPPMKRRx4Z5wiLpsR3Cqr6F6AnsA34j4h8LiJDRKR0/SZMUnr99dc599xz/QmhR48efPHFF5YQTMRUrVqVhx56iGXLlvkrnFWVF154gTZt2vDiiy+SnZ0d5ygjL2RjalXdo6ovAr2B54AHgL/FIC5jCqSqPP744wwcOJDMzEwArrzySmbOnGkDnZmoaNu2LR9//DGvvfYaDRo0AGD79u1cc801nHXWWaSlpcU5wsgKmRRE5HQReRr4DugK9FPVMTGJzJg8Dh8+zLBhw7j99tv96+655x5eeumlpGtTbhKLiDBw4EBWrVrFTTfd5O+c+OWXX3LKKadw6623sm/fvjhHGRmh6hR+BnYDr+Mm2Mk12IuqBm/qEUVWp1A2/fHHHwwcOJAPPvgAcEM7PP/881x77bVxjsyURUuXLmXw4MG5Kp4bN27M2LFjufjiixNyBNZIzLz2GTmd1xTXac1HVTUuvTosKZQ9W7ZsoU+fPv5+CNWrV+ftt9+mZ8+ecY7MlGXZ2dlMmjSJO+64I9ew5z179mT8+PGFdqaMNRvmwiSF//3vf/Tq1YuffnKztjZs2JAPP/yQU045Jc6RGeNs27aNO++8k8mTJ/vXVapUiTvvvJM777wzYXrNl7j1kfcmzUWkrvf8NBG5VUT+EqkgjQllwYIFdOnSxZ8QTjjhBL755htLCCah1KtXj0mTJvHll1/Svn17wM2id//999OuXTtmz54d5wiLJlTntXtwdQnfiMhDwFjciKc3i8jYGMVnyqh33nmHHj16sHPnTsD1SP7yyy9p1qxZnCMzpmBdu3ZlyZIlPPnkk1SvXh2AdevW0bt3b/r378+GDRs4dOgQb7zxBldccQWXXXYZEyZMyDVSa0IINnkz8COQAtQC9gJVvfUVgOXhTAAdjUfHjh1DTE1tksGYMWNURPwTrw8aNEjT09PjHZYxYduwYYMOGDDA/x0GtFq1atq4ceNc6wBt2LChpqWlRT0mYLGGcY4NVXyUrqqZqrobWKeqB7wkkgVkRj49mbLu8OHDjBgxgpEjR/pnx7rrrrt45ZVXqFSpUpyjMyZ8TZo04c0332T27Nn+Cuc//viD3377Ld++mzZtom/fvhGbza+kQiWFWiJykYhcDNTwnvuWrZeQiaiDBw8yYMAAxo0bB0C5cuV47rnnePjhh8v05POmdOvZsydpaWmMHDky5H7r169n1qxZMYoqtFD/bZ8DfXBDaM/3ngcuGxMR27dvp0ePHrz33nuAG15g+vTp3HjjjXGOzJiSq1y5Mueee26h+y1btiwG0RQu6Diwqnp1LAMxZdPatWvp3bs3a9euBaB+/frMmDGD1NRCW84ZU2rUqlUrIvvEgt2Xm7hZuHAhXbp08SeENm3a8PXXX1tCMEmnU6dOHH300UG3V6hQgYsvvjiGEQVnScHExbRp0zj77LPZvn07AN26deOrr74K+Y9jTGlVvnx5nnrqKcqXL1/g9vvuu88/2F68WVIwMTd+/HguuugiDh48CMAll1zC3LlzS9349MYUxQUXXMBHH31Et27d/OvatGnD5MmTGTVqVBwjy63QueVEpCpwC9BMVa8XkVZAG1WdEfXoTFLJzs7mjjvu4IknnvCvu+222xg9erS1MDJlwjnnnMM555zD7t27yczMpF69egk3eF44E46+CCwBunjLG4G3AEsKJmzp6elcddVVvPnmm4BrcvrUU08xdOjQOEdmTOwlSqVyQcJJCseo6qUiMhBAVQ9KoqU2k9B27txJ3759+fLLLwHfNz4YAAAgAElEQVSoUqUKU6dOpW/fvnGOzBiTVzhJIVNEquANoy0ixwAZUY3KJI2ffvqJ3r17s3r1asANHvbBBx9w6qmnxjkyY0xBwinIvReYDTQVkf8CnwC3h34JiMhkEdkqIsuDbBcReUpE1orIMhHpUKTITcJbvHgxXbp08SeEVq1a8fXXX1tCMCaBFZoUVHUucBFubuapQKqqfhbGe78E9AqxvTfQynvcADwbxnuaUmLGjBmcddZZbNmyBYAuXbrw1Vdfccwxx8Q5MmNMKKGGzu7gewDNgU3A70CzcK7qVXU+sDPELn2BKd4Aft/gxlpqWLTwTSJ6/vnn6du3LwcOHADgoosu4pNPPqFu3bpxjswYU5hQdQpPhtimQEmn42wMbAhY3uit21TC9zVxkp2dzd13383o0aP960aMGMETTzwRtNOOMSaxhBr76Owof3ZBLZgKnBtURG7AFTHZJCsJKiMjg2uuuYbXXnsNABFhzJgxjBgxIs6RGWOKIpzWR4hIO+B4wD/ZqKpOKeFnbwSaBiw3wRVP5aOqLwAvgJujuYSfayJs9+7d9OvXj88++wxwo0K++uqrCTOWizEmfOH0aL4X6I5LCjNxFcRfAiVNCtOBYSLyOnAqsEdVreiolPn111/p3bs3P/74IwB16tRh+vTpnH766XGOzBhTHOHcKfQHTgKWqurVIlIfmFjYi0RkKi6Z1BWRjbimrRUBVPU5XII5D1gLHABsqO5SZunSpZx//vls2uRyecuWLZk1axatW7eOc2TGmOIKJykcVNVsEckSkRrAVqBlYS9S1YGFbFfAxjgopWbPns2AAQPYv38/AJ07d+aDDz7gqKOOinNkxpiSCKfz2mIRqQX8BzcG0nfAt1GNyiS0SZMmccEFF/gTQt++ffn0008tIRiTBAq9U1DVId7T50RkNlBDVRNj3jgTU6rKvffey4MPPuhfN3ToUMaNG2dNTo1JEuG2PjoRaOHbX0SOVdV3oxiXSTCZmZlcf/31TJmS077giSeeYOTIkQk39K8xpvjCaX00GTgRWAFke6sVsKRQRuzZs4f+/fvz8ccfA5CSksKUKVO49NJL4xyZMSbSwrlTOE1Vj496JCYhbdy4kfPOO4+0tDQAateuzfvvv88ZZ5wR58iMMdEQTlL4WkSOV9Ufox6NSSjLli3jvPPO47fffgOgRYsWzJw5k+OOOy7OkZlktno1vPEG7NoFJ54Il1wC1arFO6qyI5yk8DIuMWzGzaMguBalJ0Y1MhNXH3/8MRdffDF79+4FoGPHjsyYMSNhJhc3yUcVbr8dAmZrBeCOO2DaNLD+kLERTlKYDFwBpJFTp2CS2Msvv8x1111HVlYWAOeffz6vv/461atXj3NkJpk991z+hACwbRucfz6sWQM20G70hdNP4VdVna6qP6nqL75H1COLs4yMDL744gvmzZvH7t274x1OTKgqDz74IH/729/8CeHGG29k2rRplhBMVKnCkyHGZd69G158MXbxlGXh3CmsEpHXgA8ImIYzWZukqirPPPMM999/P9u2bQPcnMKDBw9m9OjRVKxYMc4RRk5GRgblypWjYsWKHDp0iMGDBzNp0iT/9kceeYQ77rjDmpyaiEtPh5UrIS0Nli+HJUtg3brQr/nmm9jEVtaFkxSq4JLBnwPWJW2T1GeeeYZhw4blWnfw4EHGjBnDrl27mDx5cpwii5x3332XRx99lG+//ZZy5crxpz/9ib179/KN919XsWJFXnzxRS6//PI4R2pKu8OH3cned/L3/VyzBrKLWBj91Vcwdy786U9g1ynRI24IotIjNTVVFy9eHJX3Tk9Pp2nTpmzfvj3oPh999BGtW7cmJSWFlJQUKlas6P9ZGnr1jhs3LuQcBzVr1uS9997j7LOjPZ2GSSaq8NtvuU/8aWnubiA9vfDXi0CVKuBN1hfSySe7CukBA6BCWN1vDYCILFHV1EL3s6SQ47PPPivRybBcuXIFJou860ryvCTvsXfvXlJTUzl06FCB8VeoUIGlS5fSrl27Yv8OTPLbtSvnxB+YBMKtemvUCNq1g/btc34ed5x7j7POKjiJVK8O3lBbfs2bw8iRcO211mQ1HOEmBcuzAYKdLMOVnZ1Neno66eFcGiWgrKwsDh48GO8wTII4eBB+/DH/1f/vBU6FlV/NmvlP/iecAHXqFLx/587w8ccwYgT4rvsqVYIrrnCV0N9/D489Bh9+6Lb98gvcfDPcfz8MHQrDhoGNyVhydqcQYOfOnTRq1IiMjIyg+wwYMIB69eqRmZnJoUOHyMzMDOt5sO2+Vj6JYtasWfTq1SveYZgYysqCtWvzn/zXrQuv3L9SJXelH3jyb9cOmjQpftn/2rXujuTYY6F27dzbli93TVf/+18Xu0/lynD11e7u4dhji/e5ySwqxUciMkNVLyhRZCUUzaQAMHz4cMaPH1/gti5durBgwYKItsZR1UITR3ETTt7nq1evZu7cuSHjWb9+PUcffXTEjs9E1uHD7mQ4eTL8+is0awbXXAOXXw6FVWmpwsaN+St9V66EENdBfuXKuZNt4Im/XTu3Lh5l+xs3wtix8MILsG9f7jgvvhhuuw06dYp9XIkq3KSAqob9wM2+VqTXRPrRsWNHjab09HS9/PLLFdfCyv/o0qWLbt68OaqfHW0ZGRnapEmTfMfme5x//vnxDtGEkJWlOmCAqju9535cconb7rNjh+pnn6mOH696442qXbuq1qxZ8GsLejRurNqzp+qtt6q+9JLqkiWqBw7E79hD2bVLdfRo1QYN8h9H9+6qM2eqZmfHO8r4AxZrGOfYot4pTFbVa4qYoCIq2ncKPitWrODDDz8kMzOTM888kzPOOCMp2usvWrSIXr16sXPnzlzr27Zty7x582jYsGGcIjOFeeklVzwSTM+errhn+XLYFOZs57VqFVzuf+SREQk5pjIy4NVX4fHH3fhJgdq3h1tvhcsug5SU+MQXb9b6yAS1detWJk6cyPz580lJSaFPnz5cfvnlVK1aNd6hmSAOH4ZTT3WdvIqjcmU4/vj8RT+NGydfm//sbJgxw1VKL1iQe1uTJvCPf8D118MRR8QnvnixpGBMKbVtmyvvX7Ys5+eKFa41UGHKlYNWrfJf/R9zTOF1DslowQJ35/D++7nX16wJQ4bATTdBWRnj0ZKCMQnON9RD4Mk/LQ02by7e+514ohsKokqVyMaZDFatci2WXnkFMjNz1qekwFVXwS23QJs28YsvFiKeFESkmqr+UeLISsiSgiltVF2b+rwn///9zxULFaZ6dXfFn5IC8+cH32/iRNeRywS3aRM89RQ8+yzs2ZOzXgT69nU9pbt0iV980RSxpCAipwMTgeqq2kxETgJuVNUhkQm1aCwpmES2e7c74Qee/NPScjeZDMZX9HPiia7Ix/ezRQu3LSvLNbWcPj3/a/v2hbfftmEfwrV3L/znP/Dvf7vhOQJ16+aSw/nnu997sohkUlgI9Aemq+op3rrlqhqXsRAsKZhEcOiQa+GS9+p/w4bwXl+/fu4Tf/v2riK4sKKfrCw3hPSkSTn9FK691rVKsoRQdJmZMHWqq3dYsSL3tuOOcy2WLr/cddAr7SKaFFT1VBFZGpAUflDVkyIUa5FYUjCxpOqGdch78l+50iWGwlSu7Jp45r36t+EYEkt2Nsya5ZLD55/n3tawoRt648YbXQV1aRXJpPA2MAYYD5wG3ASkquplYQTRCxgHlAcmquroPNv/BjwO+G7gxqvqxFDvaUnBhKLqrvyeftqdvGvXhoEDXe/WevVCv3b//pyevoFJYNeu8D67Zcv8J/9jjy2brX5Ks4ULXXJ49133ffI54gj4+9/deEuNG8cvvuKKZFKoizux/wk3P/NHwM2quqOQ15UH/gecC2wEFgEDVfXHgH3+hkswwwp8kwJYUjChjBzpyonzatnSNU9s0MBV7q5dm//kv359eJ9Ru3b+k3+7dq5C2CSPNWvcQHwvvZR7GJCKFeGvf3VFS8cfH7fwiizuTVJFpAtwn6r29JbvAlDVRwL2+RuWFEyELFwIp50WfHurVlCjhis7Dmcg24oVXbly3gTQqFHydfgywW3ZAuPHw4QJ+e8aL7jAVUp365b434kSD50tIk/jxsQpkKreVMh7NwYCq902AqcWsN/FInIm7q7iH6oaZlWdMbm9/HLo7WvWBN/WtGn+k3+bNi4xmLKtfn148EG44w5XwT9mjKvkB9dzesYMdzFy++1w4YWlv7gwVHuFkl6OF5Q38yaZD4CpqpohIn8HXgbOyfdGIjcANwA0a9ashGGZZBXOeD/Vq7uTfmCrn/bt3RhAxoRSvbqrTxgyBN580w2jsWyZ2/bNN3DRRdC6tStWuuIK18igNIpr8VGe/csDO1U1ZP2+FR+ZYG65xV3FBXPEEbBjh139m8hQdXNGP/YYfPJJ7m3167shNAYPzj8fRLyEW3xUaNcMEaknIk+IyEwRmed7hBHDIqCViBwtIinAZUCubjciEjgk54XAyjDe15h8srJcr+FQbrjBEoKJHBH485/dbHGLF8Oll+Z0dtuyBe6+2xVLjhyZU9y0caNrCXfccW48qr/9DX74IW6HULDCxtbGtTa6FnfCPguYDDwazrjcwHm4uoJ1wN3eugeAC73njwArgB+AT4G2hb1ntOdTMKXPH3+oXnBB6PkBTj9ddd++eEdqkt26darDhqlWqZL7+1ehgvuO1q6d/7tZsaLq++9HPzYiNZ+Cd8vRUUSWqeqJ3rrPVfWsCOensFjxkQm0Ywf06QNff+2WjzgCRo92Zby+fgqXXeYGPUuGXqmmdNi+3bVWevpp9x0tTK1a7i6iWrXoxRSx4iPA129zk4icLyKnAE1KFJ0xEfDrr64poC8hNGjgBowbMgSmTIGlS2HePFdsZAnBxFLdunDvve47On68m8chlN27XWe5RBDOaCkPiUhN4BbgaaAG8I+oRmVMIdLSoFcvNwQFuFYfc+a4weOMSRRVq8LQodC8ubujDaWwOrFYKTQpqOoM7+ke4OzohmNM4ebPd+3BfUMfd+7s2ooXNoyFMfESzsVK06ZRDyMshSYFETkaGA60CNxfVS+MXljGFOydd9yolb5hB3r3hrfeim5ZrDEl1a6du3j59tuCt9es6YZFTwThFB9NAybhOpplRzccY4J79ll3K+5rG3HVVW5MfGtmakqDiRPh7LPzVzxXqOCGQ0+UsbPCSQrpqvpU1CMxJghVuOceeOihnHV33QX/93+JP96MMT7t28N337kBGz/4wA29fsYZrh9Dhw7xji5HOE1SBwGtcP0V/GMFqup30Q2tYNYktWzJynLDFU+a5JZFYNw4GD48vnEZU9qUeEC8AO2BK3BjEvmKj5QCxigyJpIOHHB9DD74wC2npLiJ1y+5JL5xGZPMwkkK/YCWqpoZ7WCM8SmoU9r777syWWNM9ISTFH4AagFboxyLMYDr8NOzJ6xa5ZYbNHBTJZ58cnzjMqYsCCcp1AdWicgictcpWJNUE3F5O6W1auU6pR19dHzjMqasCCcp3Bv1KIzBOqUZkwjC6dH8uYg0B1qp6sciUhUo5XMLmUTz7rswaJB1SjMm3sKZT+F64G3geW9VY1yHNmMi4tlnoX//nIRw1VWuUtkSgjGxF84oqUOBrsBeAFVdAxwVzaBM2aAK//qXG9XU113mzjtd707rpWxMfIRTp5ChqpnidR0VkQrkn2vZmCIpqFPa2LFuCkNjTPyEkxQ+F5FRQBURORcYghsHyZhisU5pxiSucIqP7gS2AWnAjcBM4J/RDMokrx074E9/ykkIRxzh+iBYQjAmMQS9UxCRZqr6q6pmA//xHsYUm3VKMybxhbpT8LcwEpF3YhCLSWJpadClS05CaNUKvvrKEoIxiSZUUggclLhltAMxyWv+fDdEsK+XcqdOsGCB9VI2JhGFSgoa5LkxYXv3Xfjzn3N6KffuDZ9+ar2UjUlUoVofnSQie3F3DFW853jLqqo1oh6dKdVspjRjSp+gSUFVbSgLUywFzZR2553w8MM2U5oxiS6cfgrGhM06pRlTuoXTT6HYRKSXiKwWkbUicmcB2yuJyBve9oUi0iKa8ZjoOnAALrooJyGkpMDUqZYQjClNopYURKQ8MAHoDRwPDBSR4/Psdi2wS1WPBf4NPBqteEx0BeuUduml8Y3LGFM00bxT6AysVdX13lSerwN98+zTF3jZe/420EPESp1Lm19/dU1OfVNnNmjgmqGeY7N4G1PqRDMpNAY2BCxv9NYVuI+qZgF7gDpRjMlEmK9T2sqVbtk6pRlTukUzKRR0xZ+3v0M4+yAiN4jIYhFZvG3btogEZ0rOOqUZk3yimRQ2Ak0DlpsAvwfbxxuSuyawM+8bqeoLqpqqqqn1rNdTQiioU9q8edYpzZjSLppJYRHQSkSOFpEU4DJgep59pgNXec/7A/NUNe69p9esgeuvh4YNoW5d6NfPXQEbJ9hMadWrxzcuY0zJRS0peHUEw4A5wErgTVVdISIPiMiF3m6TgDoishYYiRumO66WLIHUVJg4ETZvdq1qpk2DM8+E116Ld3TxZTOlGZP8JAEuzIskNTVVFy9eHJX3VoWOHWHp0oK3V68Ov/0GNcrgAB9ZWTB4sEuWYJ3SjCltRGSJqqYWtl9UO6+VNsuWBU8IAPv3w3PP5VwllxW+Tmm+hFCxonVKMyZZ2TAXAX7PWw1egDvugEcegRNPdI+TTnKPdu2gSpXoxxhrO3ZAnz45fRCOOMIVp1kfBGOSkyWFAOE2pdy92zXHnD8/Z125ctC6de5EcdJJ0Lhx6R0E7tdfoVevnD4INlOaMcnPkkKAtm1du/svvih4e9Wqrhhl+XL48UfIzMzZlp3tZhVbtQrefDNn/ZFHuuQQmCyOPx4qV47usZRUWppLCL67p1atYM4c64NgTLKzpJDHpEnQvXv+oqTKleG991zbfIBDh2D1avjhh9yPLVtyv27nTjepzKef5qwrX94loLzJokGDxLirmD8fLrwwpw9Cp07w4YfWB8GYssBaHxVg61bXFv+DD9zdwJlnukrV1q0Lf+2WLa7COjBRrFzpWu8Upl69/IniuOPcaKOx8u67MGhQTh+EXr3grbesD4IxpV24rY8sKcRARoYrVsp7V7F9e+GvrVjRJYa8dRVHHRX5OPPOlHblla7FkfVBMKb0s6SQ4FRh0yaXHALvLFavhsOHC399gwY5CcKXMNq0Ce8EvnWrO9l/8YXb//zz4eefYfTonH18rawSoTjLGFNylhRKqfR0WLEif7LYtavw16akwAkn5C+CqhMw7uyiRa5IaGe+EaYcEfj3v+HmmyNzPMaYxGBJIYmowsaNuYueli2D//0vvI50jRu75HDCCa4iPVRCmDrVJsYxJhlZUigDDhxwzWPzJou9e4v/nuvXW7NTY5JRuEnBmqSWYlWrQufO7uGjCr/8kr9Se9268N5z9WpLCsaUZZYUkowItGjhHn0DJj/dtw8efBAefzz06+vYvHfGlGk2IF4ZccQRcNttofs8tGnjhg03xpRdlhTKkHr14LHHCt6WkgITJlgTVGPKOksKZczNN8M77+TUQ5Qr5/opzJ8PPXrENzZjTPxZnUIZdNFF7pGR4ZKC9Vg2xvhYUijDKlWKdwTGmERjxUfGGGP8LCkYY4zxs6RgjDHGz5KCMcYYP0sKxhhj/CwpGGOM8St1o6SKyDbglxh+ZF0gjDnSSi07vtItmY8vmY8NYn98zVW10JnWS11SiDURWRzOcLOllR1f6ZbMx5fMxwaJe3xWfGSMMcbPkoIxxhg/SwqFeyHeAUSZHV/plszHl8zHBgl6fFanYIwxxs/uFIwxxvhZUjDGGOOXtElBRKqIyOciUt5bni0iu0VkRp79eojIdyLyvYh8KSLHFvK+nb19vxeRH0SkX8C2XiKyWkTWisidAev/KyI7RaR/pI9NRJqLyBIvnhUi8veA/TqKSJoXz1Mi4c2rJiKdRORwYLwicpWIrPEeVwWs/1RE9otIxJrW5f3beetqiMhvIjK+uMcnIt1FZE/A3++egG0x+dsVdHze79oX0/SA/Y4WkYXe7/wNEQkxmar/NSeKyNfedyFNRCp76wv8XYnI4yKyWURujeLxNRORj0RkpYj8KCItinN8InJ5wO/pexHJFpGT43l8InJ2npjSReQvxTy+iiLysnccK0XkroBtMft+oqpJ+QCGAjcHLPcA+gAz8uz3P+A47/kQ4KVC3rcqUMF73hDYipuXojywDmgJpAA/AMcHvO4loH+kj837rEre8+rAz0Ajb/lboAsgwCygdxjvXR6YB8z0xQscCaz3ftb2ntcOeM1nQGq0/nbeunHAa8D4gHVFOj6ge96/f8Axx+RvF+S7uT/Ifm8Cl3nPnwMGF/K+FYBlwEnech2gfGG/K+A+4NYoHt9nwLkB39GqxTm+PJ/RHlgfznch2scXsP5IYGdxjw8YBLzuPa/q/S+3iPX3M2nvFIDLgfd9C6r6CbCvgP0UqOE9rwn8HupNVfWAqmZ5i5W91wN0Btaq6npVzQReB/oWP/yQ/MemqpmqmuGtr4R39yciDYEaqvq1um/OFOAvYbz3cOAdXLLz6QnMVdWdqroLmAv0isiRFCzX305EOgL1gY8C1hX3+AoSy78d5Dm+gnhXuucAb3urXqbw4/szsExVfwBQ1R2qejjCv6tw+I9PRI7HXUTN9WLar6oHinl8gQYCU73PiNvx5dEfmFWC41OgmohUAKoAmcBeYvz9TMqk4N2mtVTVn8PY/TpgpohsBK4ARofx/qeKyAogDfi7lyQaAxsCdtvorYuogo5NRJqKyDLv8x9V1d+9z95YlHhEpDHQD3dVEygmx+bFkOv4RKQc8CRwWwExFen4PF3EFfvNEpETAt4rLsfnqSwii0XkG1/RA+4qf3fABUg4MbUGVETmiCsSvd1bX9zfVZEVcHytgd0i8q6ILPWKc8pTvOMLdCleUiC+xxfosoCYinN8bwN/AJuAX4EnVHUnMfx+QvJOx1kX2B3mvv8AzlPVhSJyGzAGlyiCUtWFwAkichzwsojMwt225tu1CDGHK9+xqeoG4EQRaQRME5G3ixnPWOAO7+oycH2sjg3yH98QYKaqbohATN/hxn/ZLyLnAdOAVsV8r+Iq6LvZTFV/F5GWwDwRScNdIRY1pgpAN6ATcAD4RESWFPO9iivv8VUAzgBOwZ3o3gD+BkzP98owYxKRU4EDqrrct6q471UMBZ5bvLuV9sCcEsTUGTgMNMIV034hIh8X872KLSnvFICDuKKdkESkHq78daG36g3g9HA/RFVX4jJ7O1z2bhqwuQmFFEUVU9Bj8+4QVuD+CTd6MRQlnlTgdRH5GXcr/Ix35RqrY4P8x9cFGObF9ARwpYiMphjHp6p7VXW/93wmUFFE6hLf4/P93VDV9bjy91NwA6XV8ooSwo1pI/C5qm5X1QO4eqEOFO+7UFx5j28jsNQr+sjCJeIOFO/4fAKvyH2fEa/j87kEeE9VD3nLxTm+QcBsVT2kqluBBbj/yVh+P5MzKXjl3uXFa3kRwi6gpoi09pbPBVYCiEg/EXkk7wu8FgUVvOfNgTa4CqFFQCtvewrui1vQ1VCJ5D02EWkiIlW857WBrsBqVd0E7BOR07zyzSvJKecdJiLDCnjvo1W1haq2wN3KDlHVabirnz+LSG3vM/5MzhVRVI9PVS9X1WZeTLcCU1T1zuIcn4g0CGiV0hn3/d9BjP52BR2f9zut5D2vi/v7/eiVjX+KS84AVwUcX4HfTdzf5EQRqep9R8/y3ivo7yrax4f73db2LsDAlbMX9/h8xYkDcOXqvs+M5/H5+Os4vP2Kc3y/AueIUw04DVhFDL+fkKRJwfMR7lYaABH5AngL6CEiG0Wkp3flcj3wjoj8gKtT8JVdH0PBt93dgB9E5HvgPdyJc7v3XsNw/5grgTdVdUUMju04YKEX/+e4csg0b9tgYCKwFtd6YZa3vi3uZBgWr1zzQdyXcxHwgLcuWnL97UIo6vH1B5Z7v6uncC1DNMZ/O8j/91vsxfQpMFpVf/S23QGMFJG1uDLqSd76Ar+b3glrDO5v9D3wnap+6G0O9ruKBv/xqephXDL/xCsWE+A/3n5FOj7PmcBG764qUFyOD0BcE9umuP+/QEU9vgm41lnLcX/DF1V1Wcy/n8VttpToD9wt+CsleP2rQL0IxvMSkWuSWtJjmwGkRPDYPiOyTVIT7fgi9reL0PFF+rt5H5FtsmnHF9vji+j3M2nvFFR1KfCpBHSAKuLr/6qq2yIRi4j8F3crnx6J94vAsV2grmlbiYnIp7j204cK2zdcCXZ8Ef3bQcJ9Nx8H/oqrG4sIO75CX5+w5xawAfGMMcYESNo7BWOMMUVnScEYY4yfJQVTpnjNAVVE2sY7FmMSkSUFU9YMBL7EtfWOmeJWShoTa5YUTJkhItVxncOuJSApiEg5EXlG3HDTM0RkpnhDEYsbkvlzccOTz/GGM8j7vseIG7dokYg8ICL7vfXdxQ0t/hpunCxEZKSILPceI7x1LURkecD73Soi93nPPxORsSLylfeazlH7BRmDJQVTtvwFN4zA/4CdItLBW38Rboji9rhxr7qAG98eeBrXBrwjMBn4vwLedxwwTlU7kX/4gc7A3ap6vLjRXq8GTsX1Vr1eRE4JI+5qqno6bhyoyeEerDHFYUnBlCUDyRke4XVvGVzv1LdUNVtVN+N6FoMbwqQdMNfrwf5Pco+x49MF11se3JwPgb5V1Z8CPuc9Vf1D3RhM7+LGqSrMVABVnQ/UEJFaYbzGmGJJ1lFSjclFROrgxt1pJyKKm7hExQ0vHWzGNgFWqGqXEnx0YKepYJ+TRe4LtLzj6uTtTGSdi0zU2J2CKSv64wbTa65u0L+mwE+4q/cvgYu9uoX6uBnaAFYD9UTEX5wkOXMwBPoGuO4xrrEAAAEISURBVNh7HqoCez7wF2/Aumq4uSu+ALYAR4lIHW9wvAvyvO5S7/O7AXtUdU+RjtyYIrA7BVNWDCT/BErv4IYrHoqbrnU5bnrWhbiTb6ZX4fyUiNTE/b+MxQ1PHmgE8KqI3AJ8CBR40lbV70TkJdzUkQATvSETEJEHvM/9CTcyZqBdIvIVbobAa4py0MYUlQ1zYQyuZZK6yXfq4E7aXb36hXBeWxU4qKoqIpcBA1U1ItMlishnuMHcFkfi/YwpjN0pGOPM8CpwU4AHw00Ino7AeG8s/93Y1bwpxexOwRhjjJ9VNBtjjPGzpGCMMcbPkoIxxhg/SwrGGGP8LCkYY4zxs6RgjDHG7/8BFuHH+bleDHAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Calculate the mean, SD, and sample size for BMI within age/gender groups\n",
    "da[\"agegrp\"] = pd.cut(da.RIDAGEYR, [18, 30, 40, 50, 60, 70, 80])\n",
    "pr = da.groupby([\"agegrp\", \"RIAGENDRx\"]).agg({\"BMXBMI\": [np.mean, np.std, np.size]}).unstack()\n",
    "\n",
    "# Calculate the SEM for females and for males within each age band\n",
    "pr[\"BMXBMI\", \"sem\", \"Female\"] = pr[\"BMXBMI\", \"std\", \"Female\"] / np.sqrt(pr[\"BMXBMI\", \"size\", \"Female\"]) \n",
    "pr[\"BMXBMI\", \"sem\", \"Male\"] = pr[\"BMXBMI\", \"std\", \"Male\"] / np.sqrt(pr[\"BMXBMI\", \"size\", \"Male\"]) \n",
    "\n",
    "# Calculate the mean difference of BMI between females and males within each age band, also  calculate\n",
    "# its SE and the lower and upper limits of its 95% CI.\n",
    "pr[\"BMXBMI\", \"mean_diff\", \"\"] = pr[\"BMXBMI\", \"mean\", \"Female\"] - pr[\"BMXBMI\", \"mean\", \"Male\"]\n",
    "pr[\"BMXBMI\", \"sem_diff\", \"\"] = np.sqrt(pr[\"BMXBMI\", \"sem\", \"Female\"]**2 + pr[\"BMXBMI\", \"sem\", \"Male\"]**2) \n",
    "pr[\"BMXBMI\", \"lcb_diff\", \"\"] = pr[\"BMXBMI\", \"mean_diff\", \"\"] - 1.96 * pr[\"BMXBMI\", \"sem_diff\", \"\"] \n",
    "pr[\"BMXBMI\", \"ucb_diff\", \"\"] = pr[\"BMXBMI\", \"mean_diff\", \"\"] + 1.96 * pr[\"BMXBMI\", \"sem_diff\", \"\"] \n",
    "\n",
    "# Plot the mean difference in black and the confidence limits in blue\n",
    "x = np.arange(pr.shape[0])\n",
    "pp = sns.pointplot(x, pr[\"BMXBMI\", \"mean_diff\", \"\"], color='black')\n",
    "sns.pointplot(x, pr[\"BMXBMI\", \"lcb_diff\", \"\"], color='blue')\n",
    "sns.pointplot(x, pr[\"BMXBMI\", \"ucb_diff\", \"\"], color='blue')\n",
    "pp.set_xticklabels(pr.index)\n",
    "pp.set_xlabel(\"Age group\")\n",
    "pp.set_ylabel(\"Female - male BMI difference\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Inter-group and intra-group differences:** As the sample size grows, estimates become increasingly precise, but it is important to remember that a highly precise estimate for the mean does not imply that individuals within a population do not vary from each other.  To put the differences shown above in context, below we show the underlying summaries on which the plot above was based.  Note that the standard deviation of BMI within both females and males ranges from around 5 to around 8 depending on the age band.  This means, for example, that two randomly-selected males will tend to have BMI values that differ by around 6 units.  This is a far greater difference than the mean difference of up to around 1.5 BMI units between females and males.  Thus, while there is a tendency for females to have slightly higher BMI than males, the heterogeneity within genders is substantially greater than the difference of means between genders."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              BMXBMI                                                         \\\n",
      "                mean                  std             size              sem   \n",
      "RIAGENDRx     Female       Male    Female      Male Female   Male    Female   \n",
      "agegrp                                                                        \n",
      "(18, 30]   28.123881  27.391822  7.745893  6.649440  609.0  544.0  0.313879   \n",
      "(30, 40]   30.325586  29.611726  8.315608  6.622412  474.0  458.0  0.381949   \n",
      "(40, 50]   31.160643  29.724623  8.076195  6.407076  502.0  401.0  0.360458   \n",
      "(50, 60]   30.743777  29.231486  7.575848  5.914373  470.0  454.0  0.349448   \n",
      "(60, 70]   31.074828  29.392488  7.604514  5.933307  441.0  437.0  0.362120   \n",
      "(70, 80]   29.138213  27.957692  6.284968  4.974855  410.0  402.0  0.310392   \n",
      "\n",
      "                                                             \n",
      "                    mean_diff  sem_diff  lcb_diff  ucb_diff  \n",
      "RIAGENDRx      Male                                          \n",
      "agegrp                                                       \n",
      "(18, 30]   0.285092  0.732059  0.424026 -0.099032  1.563150  \n",
      "(30, 40]   0.309445  0.713861  0.491570 -0.249616  1.677338  \n",
      "(40, 50]   0.319954  1.436019  0.481976  0.491347  2.380692  \n",
      "(50, 60]   0.277575  1.512291  0.446275  0.637591  2.386991  \n",
      "(60, 70]   0.283829  1.682340  0.460097  0.780550  2.584130  \n",
      "(70, 80]   0.248123  1.180521  0.397377  0.401662  1.959380  \n"
     ]
    }
   ],
   "source": [
    "print(pr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Confidence intervals and sample size\n",
    "\n",
    "Confidence intervals reflect the precision of an estimate, which is largely driven by the amount of data used to construct the estimate.  We can explore the relationship between precision and sample size by subampling data from NHANES and calculating confidence intervals for the subsamples.  Below we calculate confidence intervals based on subsamples of size 100, 200, 400, and 800.\n",
    "\n",
    "A wider confidence interval implies that we have less precision in our estimate.  In the simulation below, we calculate the average width of the confidence intervals constructed for each sample size.  We see that the confidence interval steadily becomes shorter as the sample size grows.  For most settings, the confidence interval will become around half as wide when the sample size is increased by a factor of 4.  Below we see this scaling when the sample size increases from 100 to 400, and when it increases from 200 to 800, both of which are increases by a factor of 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 3.020821172332006\n",
      "200 2.1396241183284594\n",
      "400 1.5173780418171638\n",
      "800 1.073859723411907\n"
     ]
    }
   ],
   "source": [
    "dx = da.loc[da.RIAGENDRx==\"Female\", [\"RIAGENDRx\", \"BMXBMI\"]].dropna()\n",
    "\n",
    "all_cis = []\n",
    "for n in 100, 200, 400, 800:\n",
    "    cis = []\n",
    "    for i in range(500):\n",
    "        dz = dx.sample(n)\n",
    "        ci = sm.stats.DescrStatsW(dz.BMXBMI).zconfint_mean()\n",
    "        cis.append(ci)\n",
    "    cis = np.asarray(cis)\n",
    "    mean_width = cis[:, 1].mean() - cis[:, 0].mean()\n",
    "    print(n, mean_width)\n",
    "    all_cis.append(cis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also informative to plot the individual confidence intervals, computed for 500 subsamples of size 100, to see how they vary.  The vertical grey bars below each correspond to a confidence interval.  The red horizontal line is the mean BMI calculated using the entire data set, which can be taken as a proxy for othe population mean.  While the individual intervals are quite different from each other, it appears that the vast majority of them cover the population value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fbffe409cf8>]"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD8CAYAAAB6paOMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztnXuQJVd52H/f3dlZaSPwMmjISEibxSXKYYuYFUwROUoc2MjOhHJcqRiMcUJQAlmSMmWwKRzIi1Iq5QTHEa8Qmy0rQKoE5iGpTJRosQqkwkoZCS2sxGpXgHChikoaJDMSYtlEo5358sftXnp7+93ndJ/u/n5VUzNzb997Hn3O9zrfOS2qimEYhjFdZn1XwDAMw+gXUwSGYRgTxxSBYRjGxDFFYBiGMXFMERiGYUwcUwSGYRgTxxSBYRjGxDFFYBiGMXFMERiGYUychb4rUIWLL75Y9+3b13c1DMMwBsXRo0f/XFWXy64bhCLYt28f9957b9/VMAzDGBQi8nCV6yw0ZBiGMXFMERiGYUwcUwSGYRgTxxSBYRjGxDFFYBiGMXFMERiGYUwcUwSGYRgTx5siEJELROQeEblPRB4QketS739YRE75Kn8sHDlyhCNHjvRdDcMwRozPDWXPAAdV9ZSI7ATuEpHbVPUrIrIK7PFY9mhYX1/vuwqGYYwcbx6Bzokt/p3Rj4rIDuA/Ab/lq2zDMAyjOl7XCERkh4gcAx4HblfVu4G3AZ9X1cdKPntIRO4VkXufeOIJn9U0DMOYNF4VgapuqeoB4DLglSLys8DrgA9X+OxhVV1V1dXl5dIzkwzDMIyGdJI1pKpPAXcCrwauAB4Ske8Cu0XkoS7qYBiGYWTjM2toWUT2RH9fCFwDHFXVFVXdp6r7gNOqeoWvOhjjwDKnDMMvPrOGLgE+ES0Oz4DPqOqtHsszRoplThmGX7wpAlW9H7iy5JqLfJVvGKETezlra2s918SYOoN4MI1hjBHzdIxQsCMmDMMwJo4pAsMwjIljisAwDGPimCIwnGKpnoYxPGyx2BNTzQixBVDDGB6mCDxhAtEwjKFgoSHDMIyJY4qgARYHN4xq2FwZBhYaaoCFfQyfjGl9yebKMDCPYABsbGyYVTUh1tfXTYAanWKKYABsbm4OQjAcOXKEjY2NvqthGMEwlNCYKQLDGevr62xubnr5bvOKwiNLyA1F8HXFULw7WyMwBsFQvKIpkXU/7B4NE/MIamLhD8OYDlPxcMwjqInP8IdhGGExFQ/HPAKjlKlYRYbRlqHOFW8egYhcAHwZ2BWV8zlVfa+I3ACsAgJ8C7hWVU/5qofRnqlYRaEzpv0FY2Woc8VnaOgZ4KCqnhKRncBdInIb8Buq+jSAiFwPvA34jx7rYRjO6UMoD1XIGPnkjaOux5fPZxYrEFv6O6MfTSgBAS4E1FcdjPo0HYBTW0Q3oTxOuhbAeeOo6/HldY1ARHaIyDHgceB2Vb07ev1jwDrwl4EP+6zDkAghvtg079nXIvrUFIzxY/qYD13m/Yc0tr1mDanqFnBARPYAt4jIS1X1uKr+YxHZwVwJvB74WPqzInIIOASwd+/exnUYUlzVrMzzsSyt6dLnfKiqgNrIl5DGdidZQ6r6FHAnsJZ4bQv4NPBLOZ85rKqrqrq6vLzcuOyh7OwzDCMcqsqNscgXb4pARJYjTwARuRC4BvimiFwRvSbA3wUe9FUHwzDcEkL40jcbGxvBhGy6wmdo6BLgE1EIaAZ8BvifwJ+IyHOZp4/eB/xzj3UwDMMhda3fPgVq07BNKOGaLvGZNXQ/cGXGW1f7KnNshLSY5JIhrdsMgZD7s0uhmu6HPKXlq78++tGP8v3vf59LL7208Xf05W3ZERMBE9Jikku6jqmGLChdMIYYtQvK+iEeB77668knn+TMmTOtvqOve2mKYMKMXUDGZE2uqbTd+DFthOzYx4spggkzBkuy6QQdQ9vHRsjC1uV4CXGx3RTBwAl58nSBCfTxMJV7GWI7TREMnBAH1RiYuoKdCnaf55giMIwM6irYMoESosAJsU5dY4bUnMk/j2AKG2SmSNf3tWyHaWg7UI8cOcKJEyeCqpML7NnWzZi8RzC2iWDMsftazFhTk0N5tnVVZXTmzJkg9gpNXhGEQrx5bGlpqe+q5GKhhO6ILdu1tTWzcAdIVWW0vb0dhEI2RRAIQ7DQ+ra0ksJx7CQt27773TifsSlnUwQdEqpFHWq90oTi9kO/Z+gY/RPKOHSFKYIcfAjHUAdPqPUKma69t1jxpEOHeSHFoSh3IwxMEeQwJeFo1m0xIQjVPMWTF1Kc0vg12jP59FGfDCWVbXNz09tjJtPt7/pEVRflnThxghMnTjiq0XiY4rn9Y8U8ghKaWIOx8Akppp3ElXIq65v19fWzgiJ5NHCXYRUX5YW+iN8X1i/jwTyCHGJrp8lGIFfCbmNjo/Wxtlkk29SmjCp9U6YMbUNfWNj9KMe1V3vmzJnMPu8yomAeQQ4hWDubm5tsb2+f/d9HrDpdRtfUUbIuNt+YkCum7H6EGAqK76kvwylNVUMvb4E/zfb2dma/dxlRMEUQEcKCYBGxhxLyhjPfuNh801eoLkQB2oQ+sqXK9o7E97QPo6bovoZgTFbFmyIQkQuALwO7onI+p6rvFZEbgVXgWeAe4K2q+qyvelQlxFh+kr4HVVXrxsim7/tXhAsjyJenFeo6W0zI97UOPj2CZ4CDqnpKRHYCd4nIbcCNwD+Mrvkk8Bbg9zzWIxOXA7evcEOXx1KMZcAPHR+7q10I2qrfkaV0LFzXPz4fXq/AqejfndGPqur/iq8RkXuAy3zVoYjkwG1j7canOPZhKTddlK7S3q5DZeny2ggHFwqyj1BOlTaHZiFnLXTmKau8uVLUnlAOZcsiXpNYWBh+hN1rC0RkB3AUuAL4iKrenXhvJ/BG4O0+61CFNtbuEM4ISlOlvl0Im6TATpfXpvyuU0ZdeWZt+7yt8m7SjqyFzjxl1eS+hHIoWxZ9J1q4xKsiUNUt4ICI7AFuEZGXqurx6O3/CnxZVf8k67Micgg4BLB3716f1TR6YohKNIt0O/pIPHCRTNDl/ajbR6Ekc9RNHQ3Zo0nSiU+jqk+JyJ3AGnBcRN4LLANvLfjMYeAwwOrqqnZRzyEwpRM4iwg5rlxm2fuo+9AUapO9Ob6oM6fqKsuQPZokPrOGloFnIyVwIXAN8D4ReQvwt4G/papB+lVHjhwJNvYXWozYNUeOHOHhhx/mRz/6UaF1m9UHISuHJGO+f0PE9ZwaggeQxqekuwT4RLROMAM+o6q3isgZ4GHgT0UE4GZV/Xce61Gb9fX10cT+hsb6+jpPPvlk4892Sd6O0CK6yPQqKsOXsgw9vdgy7IrxmTV0P3Blxuvhmdkt6Go3Y5quD28zzidvR2gRvuPw8VpB3ppFnfrWicv3JfyqCvixrEf5YlRCuQquLKI4rlg1cyA5YF0sfGUN7KIc7ZDWFOIFtFCtx6ZUORvGt+HQ9kjq5OJm3ufyPKH0wmgXIdY8Ad/EA/BtXBWNjb5PKp6cIiiaFPHN2NjYYL4NIp+6ccXkgPUVwshL2QuNrAW0vNBCl7naeYq0avlVxkTacOhi8b+OcKuyuJn0hFT17PenP9tniDWtIKoo4PX1dU6fPg3gZbzFR5lnKae+1/4mpwiKiG9GGxdyCuGapOXiypLJ6/Myj8tlf+cp0qbCrIo31oUA8B0SGULIparnXnZNG8Mk7qcQ11NMETgmhEnhMxyUzlfPE2RdrWGE0N95VBXwfYcFpkQdIZx1rYtNZCGOWVMEjondzz5TT+vGhMvCYEmqDuJQFueKQgKhCOCkMvWlPMsynELwZEMzHLoav6raS8JJElMEjklaC3V2FfaVfVQWdglBULahyILrOy6bhS/hU5bhFILS7qIObeZYPEdjwe3S2Os7XX0STyjLE2hlz1xtuz08b1E06ztDOLckHc7Z3NzkxIkTQViLxnw81hVk8WK3MWd7e7vxPEvO0b7nqmsm4RHEAi3r9SJ8bA8PwfLKIyucE3J9fRCSF5SuSxPhM4bNkaGmGxeFVYemfCehCKCbOOyUaZNNUeezTXbzlpUNlC5+d0G6H/qsS1c7cauUk2WQ+VTWVUM/RQp2aMp3EqGhNJubm+cNrCZud1vG9KDw2G1uEk6rExZrspu3iNOnT5/NHe+bEMKDMUnv0Oc4TXuhVRX9+vp6rYfh1J3bodyHrpiMR1BG3o33eYxsaAuVVRasy96vE04L4cEeU5vwTUiP0+QmsjyazpksRd82q2YM4THfmCIoIRZsXVnuPnealrUhzyJNpsS2mVBZi9FdTtAuDx7zSQhHkZcp+83NzdqpyUnSYyUeJ1VPp21SxpQxRVCRptZ7FVc36bq6jAunBV+V742tveQkcyWs+95b0Hf5MVl9XIcQ016zaDNu8u5Vm9Npq5bRllgBRqcrDwJTBJ5Ju7pZLq4L1zUt9LNOoaxKCMKyiL636LtIyQy9j33i+/51beWnjb0hhqFMEXSMr0GStm7GLGjK2la21pEnKKpu6nMZcx5amqELqozNrHtRVcBX/X5XuE5g6IPJZA21WfSNP1u2Ac0lrjcCTSkeWrb2kJU1lvW55JjJC/FVGRN593JjY4Otra3z6uozg83VGPadZZd1D/PuW9Pvr/Ndfe3874rJeARtNodtbW1x+vTpTrNb8qzOpptruo6Pd5kR5HpvQUxyzMRWX9aCdxl59zLvmGSfoYW847/rMsTwRxtcJzaE9lD7ySiCtvQ18NOpc0N5GHbViZNsX9PJkeea+7DifCnUdH/VSZlsc9BhKGNp7BZ3mtDmsbfQkIhcICL3iMh9IvKAiFwXvf42EXlIRFRELvZV/piIz0cZqxWWPL+l6eTICnmEtEGrCXU22Q25neDnXoVmdScJ4cTRJD49gmeAg6p6SkR2AneJyG3A/wZuBe70WHYjqmyUSVpffVsxybL7rkvf9GldhXD0eB7x0/ZEJIgNfF0SmtWdJiTl7fPh9Qqciv7dGf2oqn4dCDbHtspBdPDjxbLt7e1WeeFtTzeNqfIkr83NzWAtpCGT1e9ZKZJFirpu7nnVw/HSmWQhCR+X+DKEqnoVQzfEvJoGIrIDOApcAXxEVe/2WZ5vkpM1K6Mhj6KBlNyB6UI55pW1ubmJqjq1kIY2+NtaxHWUaFY/NzmkrOh+dp2yGMKR1mkvOHlYoA8lt729XeksqqErWa+KQFW3gAMisge4RUReqqrHq3xWRA4BhwD27t3rrY51Bnabc8x9fG+TsupSpKTKjqQIiY2NDU6fPl2pr/MUXB9hhip7JrraxRrCmT1pLxj8j7e+29wFnewjUNWnmK8JVD4cRVUPq+qqqq4uLy97q9sYFtrKaGO5N+kbV33q0uOoY7E1yXjqi5Bj4Gl8eZC+5nCbs5KGhs+soeXIE0BELgSuAR70VZ6Rjw+3NTmpfU3wIbjbIdWvzSYvVxvEir4n636GHF7s6952uXE1xqdHcAlwh4jcD3wVuF1VbxWRXxeRR4DLgPtF5A881mGQ9DEQ6nL69Gm2traAagvVVTd8JQWJK4ssabmXWfEhWPlNaWMZu7Kq63xPnXDdlHC5g7oqPrOG7geuzHj9Q8CHfJU7FMoWkIsoy2ToIp2xzuRNLmyWCffk97oUEOnvLRL2Jpjck7XQPASPr0v6NP6mkVAcIG00fll+dIiTKx7kvurm4wlUTbyDvj25qo9Z7BpfC81dLpa3IVaERfelz/WesEaLZ6rkBLsKC/R9JnmXAiGvrck6+B7kvhRM3e8NYfE2REPAF229566oogj73Jg4KUVQZaehq0nUx2TMOpcoD5cTZAwP8Q550bIPxtIfIewurjrX+pwrkzmGGoa9EFiVOufT9D1BQiLUs26KsnDqCus6i+8Wv3fHEObapDwCGI6FanSHLwOhSABULS9vvCaPOKlTn1AYu0GWR6jtLlQEIvKbRe+r6vVuq9Mdod6QrkgLv5D6I0sw+15zqSokXfVTG6F85swZtra2CvsilNh4Hk3aX2UMhLpYHhOSMk5S1lvP6aQWPRDCDWkqVFxZsL5SNdM0WQRL1yeE+wX91aPonKssqoQjQj41NYs6Yc8qpE8SDomuDbPCEaCq13VVkSnSRqiEIhir4Gv7f2iT1ye+DlSrQtnD5kPyJuuQdW5RKHQ9v8tCQ4Ubv1T1191Wxy2hPKc3JMtrqJM2TWgLcKGHJNrQ5aGJfRDS/OyLspb/M+A48BngUSDsXRspun5Obx5tJ8psNhtMWutsNqtUzlgUUpKhC0RfhO692X0rVwSXAK8DXg+cAT4N3KSqT/qumDFHREo3ZIWU8x1bVWO3Io3qhOa9GedTuI9AVb+vqr+vqq8GrgX2AA+IyBu7qJxRjqrawV0ZTGHPSJe4tOq7Pt55SMdJ9zVuKwXFROTlwBuAnwNuY/7UMSMQpqoEyibMVPvFBy6t+i7uS90sq7YC2GV6cx/jtmyx+DrgF4CTwB8C71FVM7OMIDBBb+TR9IFKXZYXEmUewb8B/gx4WfTz25HGE+YPov9pv9UzDKMpRRZuyKd2Vl3zstCfO8oUwYs6qYVhOCKkhfO+KbJSQ168rXrO0dCt8JAo21D2cPo1EbkY+L4OZfWlQ4qsrKkvXnbVfjssrTtms9nZp9S1pe9j26dOYdaQiFwlIneKyM0icqWIHGe+r+B7IlL5QfRTwVImi5li+0PPoW/DwsKCM8E9xbEREmXHUP8X4LeBTwFfAt6iqivAzwL/oeiDInKBiNwjIveJyAPRwjMi8iIRuVtEvi0inxaRRQftCI54Y9WUmM1mLC6O8nY2Zmtri9OnT/ddjd4JLYXTUljPpUxaLajqH6vqZ4F1Vf0KgKo+WOG7nwEOqurLgAPAmohcBbwPeL+qvhh4Enhz8+qHSxtrqegM+pBZWFjIPY9mypi1G14fdF2f0NqfpkwRJGv/f1PvFao3nXMq+ndn9KPAQeBz0eufAP5etapOh+3t7eAHztSZosdnjJeyrKGXicjTzNNFL4z+Jvr/grIvF5EdzDefXQF8BPgO8FRiL8IjwAtzPnsIOASwd+/esqI6xeXZP8YwWVhYyBwDdQ4wc7nY6pupJzs0perZW31TljW0o82Xq+oWcEBE9gC3AC/Juizns4eBwwCrq6uNgmuv/OQnueihhwDYVgVVEGEmMv///EIhDufk/Q3nfz6R7aB5ZVT4jry6FX5Xyd+z6P9k+8s+F8f6Nzc3i/ut4ufq9G3VNmXev2RftSwv2iiTe88q39+C8vK+I6vsKv1StX3nCKeM8pJjZseOHWxvbZ1bn5r3QUTmk7xiPSuPd8f94uNzTeZfsrwLdu2CY8fgAx/AJ534t6r6FHAncBWwR0RiBXQZ81NNB49SEitrgcxmwzr21QNdt19EzhUAqtUX++pcW6VszxTVVLe3241rEUsJHQDeDuAWkWXgWVV9SkQuBK5hvlB8B/Ba5kdWvAn4I191uOdXf5X19XVgntoZ5ynHVmuaZB5z3t/AeZ/PurbomqLviMMFyWv37t17zpHaVeqZrEe6/WWf27VrFysrK2fLrNKmos/V6dus93fs2MHu3bs5depUZt8lid9rU16b+1envLzvaDM+q7Rv165dADzzzDOZn0uOmaL+rHIfkn+37Zeq7Wtz311/rsn8S8/9a6+9NrOPXeLTI7gEuENE7ge+CtyuqrcC/wL4TRF5CHg+cIPHOgyO0GOJfWDZSP1g6wL1mc1mg0wk8OYRqOr9wJUZr/8Z8Epf5VbBdjHWZzabWb9NEDNM6rF79242NzcHkwQQMzzV5YC2g3txcXGQWr8p8cNxTCgYxvnEskBEBuu5TkeaOWRpacnL801DcitDPrwt9F2aQ6fvDY2h3d+iORkbSV2V54swpI4BzGPhrgdVk0l9NgU0UA8g1HqNhb43NPZZdiyEY6PMh6DPKi/Gd3m59ei8RKNT6k7qIbu3trgZBm3uQ98e8cLCwtkstS4EssuD+9pgisAYFeYthEGT+9CVNVwWgu0ySy1vvbHrMHH3PogRPCGvDwyVkI4amPq9jZVNCJk9S0tLbG5unrcfo+vwkHkEI8C15RDa+sAYjrb2sf7TlL7XAHzTd3jJFYuLi6ysrHRS1jh6bML0tbjUFW3WLMYiEJL03aa+yy9iNpuxY8eO3PkQUlZeFZaWllhb6+b5X8PplR7wOWhCSJEb2sTIIq8NZQpyiG0valNXbelycbPuPSrzutp4ZS7HS4jjLrwaBYJvS7uKa54efK6Fl+twRR8hnKZtCClU05axeoUu7pGr+eJyvISSKZTEFEHApAdfyMJryGmnoTNE7yUExqogfWCjyxgMIQnELo8ZCdkASBPK/THqYXetBBvY3VBFyIckELOOGbGxEmbYI1Rms1lnWUFl2Mgtoc+Bvbi4GFTqpEtBl/4u30Let5B2HYao6v3kXVf2+bEprSG2Z2FhobOsoDKG13sTYmlpyUvcPU9wpP/PEtYulGJTodlUMQ4xVlxVMeZdVyWDZiyWeyj3NyQLvy6mCCZIWghkTSRXFnrduH7RtU0UY0jrCsa4CcnCr4vNEMMrdRWKa0s1pHWFsTG153LUZUh9462mInK5iNwhIidF5AEReXv0+stE5E9F5Bsi8j9E5Lm+6gD9n63uE1fW7mw2C2otYsp0KTzajp8qz+UoGltDEpR1CSVcVRWfd+IM8E5VfQlwFfBrIrIf+APg3ar6V4BbgHd5rEMpQw4duLJ27ZnAYdC18OjCW8obW6EIyhDmfwielbfSVfUxVf1a9PcPgZPAC4GfAr4cXXY78Eu+6gB+t53nMSQLO7TMJKN76gghH0KrT2FcNP/L2uqq3r6eeFiHTnpfRPYxf5D93cBx4Bejt14HXN5FHWJ8Drr4e4dkYfvKTDLOFbAhK9w66zJ5QqvNvAp1HadIQM9ms2Dr3QTvikBELgJuAt6hqk8D/4R5mOgo8BxgM+dzh0TkXhG594knnnBWn7yb50I5jCklrwv6sgKrlOvC8k2OB98Kt+/wwpiEYhmhhLVc4nXkiMhO5krgRlW9GUBVH1TVn1fVVwCfAr6T9VlVPayqq6q6ury87LOawLmT1pfXEEI80jV125Psgz4UZ9VJHIK7XodkfZuOsS7GZ90yQvWiXHt4fStyn1lDAtwAnFTV6xOvvyD6PQP+NfD7vupQRl7H17FuXB6TOzSqCNX0ZBliHwxJebexVrtaPK5TxtLSEiISnBHl2sPr2/Dw2bNXA28EDorIsejnNcAbRORbwIPAo8DHPNYhF1funUurNlTrpyljOZHUh+dSV6j1bTHm0VW9qiiQ0NZhQqtPEd5UkKreBeTNng/6KtcV8Q1MP0vUF7HQXF9f76S8unTdH0OiriBsYoTEz7b1zeLiImfOnKn8KMuu6lWF2OhoOofibD9X7cmrT9zHIRGeiREIY8mmyXKpm1hxyf5w/bSmoacj1g0l9m3Zlx3j4eJhMH23sQk+sv1WVlbO8Qpms1nvYaAshne3PJB24bp8aHRd6rqaWUKqzUCczWbs3r3b6dOa2nxXlvsdwjqEi/UnH7jOeEkKuvhelOXm90Ve2T7rtLa2dla5hJxtZIqA863/Kg+NbhMbzbKYyk4ujCdZcvEsawD7tsRC2iMRW1eh1CcmxAnva1wkBV3Zvehzzaio7HhONSHLEAl1PaeIYdXWMW0Wc2Kruu5NX1xczLSoy04uTE+yLIEcogDyxVjbWjYmmwiZIffVyspK0EI1S/ktLS2xe/fuYKMKWYTbwx1Qx5rMi3vWDbOEaMEW0VW8Ny0A88ptIxiG8LCWovGxuLjI/v37zxtv6VBmsi/7jNe7OGplbW3tbHuHkoED2VGFkI+e6X/k90CeVVU0aarGdkO+2XXJ816KSC+OVSXL48kqNykY6tL0YS19ufrpcvNClunXk33pY02ian+4DiO2CeGEQEhh1TSTVAR5VryLSRPyzS4jrQireC9poZCMGbuqUx0Xu423kEdVr6+twki3ta636etAuDQhZr34JvayfHhYKysrvYeRJqkIXBPa4lCXh3/5Fgp1nvoUx8Lrrv3UiaEnBUGWxV5301N8eFn8u80Trlzci6RHl9cvTb2+PqlrUKSJjSIfHtba2lrvTzYLR3p5pkw4thncPk5kbEMIh98Vtd33jkuf6zBJQdBE8GaFwC699NLeLew4pFnFo3Pt9bUhb5ylX2+rZMfOZBRBmXD0Mbh3797N7t27z/4fL+r17Qa2oeoeiyLLqUxQd7nOMtTNT64pCmn21T9VxkDeOMt7fUjHPnSJzQCPJAVenMccghtYRNlEqbLHoi2u11ny2pS1GG6K4Xz6OiXWh9exf/9+9u/fn/t+WWRgrIrERvyAcRGrTXoocXpiegJ2tdPaR+w5r00wV2r79+8/Z2esy13TXdDXLnhfnq3v9bYyQ6wsMjC09O+qmCIYMEWDtqplm5wYVdMTm1IW8vERniure52dsT7IEuRVFWK8s7rNvcmzcMvGjy/PdooZSSFgiqAheZM1PbGrbICqEhOva4H1faZNFkNOrfVFliCvohBd7RbOU351x08Tz6SrbLu6a05151rou5+rMOza16Dt2UDJPGLIn6zpiV1lA1QVAenKAvM1YJuk5zWNt/aRrpsss204po/0S999VuaZZHkYXVn/dQ2QunMtPceHmBAyGUVQZ1NQfBPjCRsPJB8DN29yNtlIlXXcdFLg+Dxzpkl6XtNQTB/hg2SZSaHXRJn1kX4Z17+rjKx0vywsLJw9f6doZ3/WmA/lyIw0eQI/9ISQLMLo0YBITvK1tbVzFhNdUySY6wrWvOOmfQqctpt0YD6Znve857X+HpcCoo7FX6WPXfSTK7oKz+UdxhYrwbyUz6wxn8z0cR3ybDO3fQn8PjKTJqEI2rjyaeut6yyNkN3MMmVVpa/W1tZ461vf2npC1REQdVJkXQjxNpuZXAiFkJ+vUYUsgVtlw2JZ31VJUXXR/3XXEPpIWvD58PrLReQOETkpIg+IyNuj1w+IyFeiZxjfKyKv9FWHmKL4ZV1B20UefZK8SZB1Bnpo+c3pvgqlfnUmWt87Ul0IBZdZX3UXUcuuL8paKgoTpRXwuwi/AAALBklEQVR/UtnFfeaq7+LviNuT1a6i+Ze1ThjafPUZaD0DvFNVvyYizwGOisjtwO8A16nqbdHD7H8HeJXHehTSZoLEg6Htc4brWmxZ7n38f1ZdFhcX2bFjB88++2yrerah6gahlZUVNjY2nJcff2/e82jje5Dsvy4narosV2PLJXWVYnztxz/+8dxr8sZtXllLS0usrKxw4sSJc+5lrOyKykqyuFj/ucFlzwupQ2jZc948AlV9TFW/Fv39Q+Ak8EJAgedGl/0E8KivOvgmttabZIEkF8zSi4+uWVpa4l3vehcvf/nLg7JCsshbSG0bIquyUSg90V276Ok2ZFmxyfrmCZ46GUBDDwtl4WKx3fYrnEsnPSEi+4ArgbuBdwBfEJHfZa6I/loXdQB/Ft7a2hrr6+u1rLelpaVM63RpaYlTp07lfq7MYo7bmPXdcT1d06Rf6wqnKo8Ohex2J8uM+86XcIwtzaz6ZCmaJh5pHSHYdSjTKCZETw86UAQichFwE/AOVX1aRP498BuqepOI/DJwA3BNxucOAYcA9u7d66QufbljRW5oXYutTJgXhYh80aRfXQunKu1O9p0v4Rgr+Lg+RUq9bTlFZBkMIXuDIdbNh9CuEi7rA69ZQyKyk7kSuFFVb45efhMQ//1ZIHOxWFUPq+qqqq4uLy/7rKZ3itzQviy20BarQqRqOCovBbbPHadZ4ZMQz8mJ+7hu3brYlDfE/QBN8eYRyPy4whuAk6p6feKtR4G/CdwJHAS+7asORj5VJ11VQeiDqhM9z3Kr6m3lXVMmBOLPJa9LWnpra2scO3astHxXtFlo72stoamFnAzHpjeBwo/Hgq/kg7HhMzR0NfBG4BsiEs+Gfwn8U+CDIrIA/D+i8E8I5E2EUON6ebgc/EULlkXXuBAsVZVVnjCp6m01tfq6tBar9GWbNaCuPdMm4yPv+vQmUPjxWKjaJyE8a3xxcZHt7e3zste6wJsiUNW7gLxDzF/hq9wkdQVi3kQINa4H2W30tSicpEhIx+7+2Nxql4KiqhDM8jq6pMoifJKqc66J4nHdB8mU4RAORIzLr5MG64pR508VCUSXFkmf5LWxTy8mPWGb5GxDeP3tUlCk+6hpeMo3dRe9uzBCkrTxPGNldOTIkULlFdo49MGoFUERIVgkPsnzYvpQEPv37+fEiRO1P+ejv8smdV+Tvouxlb73WZvohoaLkFaZ8nJxb6p4k30qnMkqgi4JacGqjzBXXSvR54Qom9RVF4iHSPre192NazSnijfZp6FpiqADunaXh07InlcIdRta8kJd+l60nSKmCHpkjNv/y3DZ3qL+G3O/du3Vrays8Mgjj3gvI6bvRdspMilFEJollRffbJrK1vXx2E1waVEXxYdDsNzHQhf7Ifq+X10bDqGtVU1KEYScBpqkbiqbixRD12f/GP4Zs9fjkqp7MLqk7VqVayalCMZCemC7GDQm2P3gM94dyj2rImizdv92RSj9FDKmCDzRZ+aLEQ5DjXfXUWBVxmPW7t8xMXTvzBSBJ0If7EMfuHlMcQE+SZO2Z30mrcAsk6eY0Od7GaYIHDIkITT0gZtHaOfvhxoGSY7Vqha9MV5METgkNCFkuKOOQK8rZPug7lgdQkaa0RxTBEYmoXs3XdevjtAco0HQZXtclBXy2A0RUwRGJqELs1DrZwIoDEIbG6GPC1MEARK6NW7kE5oACgEby+GPC1MEjnA52EO1do1h0XRMut6B72ssm4JxhykCR5jgNkKj7ZPXQt+B3+VO+rFjiqAjxjDwxtAGwwAz3NL4fHj95cB/B1aAbeCwqn5QRD4N/FR02R7gKVU94KseodDXwHMpvG3yjANT6EYanx7BGeCdqvo1EXkOcFREblfV18cXiMh/Bn7gsQ7B0tVknJLwNgFXjapjwvpzOvh8eP1jwGPR3z8UkZPAC4ETACIiwC8DB33VIWSmJKC7wvrULdaf06GTNQIR2QdcCdydePlvAN9T1W/nfOYQcAhg7969nmtoDAmzVA3DLd4VgYhcBNwEvENVn0689QbgU3mfU9XDwGGA1dVVdVknEyTFhN4/Zqkahlu8KgIR2clcCdyoqjcnXl8A/j7wCp/l52GCpBjrH8OYFjNfXxytAdwAnFTV61NvXwM8qKp+H4RqGIZhlOJNEQBXA28EDorIsejnNdF7v0JBWMgwDMPoDp9ZQ3cBkvPetb7KHTqhx+ebMtZ2GcYYsJ3FgTHW+PxY22UYY8BnaMgwDMMYAOYRGIbhFAsDDo/RKwIblIbRLRYGHB6jVwQ2KA3DMIqxNQLDMIyJY4rAMAxj4pgiMAzDmDimCAzDMCaOKQLDMIyJY4rAMAxj4pgiMAzDmDimCAzDMCaOqDp9+JcXROQJ4OGGH78Y+HOH1RkC1uZpYG2eBm3a/JdUdbnsokEogjaIyL2qutp3PbrE2jwNrM3ToIs2W2jIMAxj4pgiMAzDmDhTUASH+65AD1ibp4G1eRp4b/Po1wgMwzCMYqbgERiGYRgFjFoRiMiaiHxTRB4SkXf3XR9XiMh/E5HHReR44rUlEbldRL4d/X5e9LqIyIeiPrhfRF7eX82bISKXi8gdInJSRB4QkbdHr4+5zReIyD0icl/U5uui118kIndHbf60iCxGr++K/n8oen9fn/Vvg4jsEJGvi8it0f+jbrOIfFdEviEix0Tk3ui1Tsf2aBWBiOwAPgL8HWA/8AYR2d9vrZzxcSD9xJ13A19U1RcDX4z+h3n7Xxz9HAJ+r6M6uuQM8E5VfQlwFfBr0b0cc5ufAQ6q6suAA8CaiFwFvA94f9TmJ4E3R9e/GXhSVa8A3h9dN1TeDpxM/D+FNr9aVQ8k0kS7HduqOsof4GeALyT+fw/wnr7r5bB9+4Djif+/CVwS/X0J8M3o748Cb8i6bqg/wB8BPzeVNgO7ga8Bf5X5xqKF6PWzYxz4AvAz0d8L0XXSd90btPUy5oLvIHArIBNo83eBi1OvdTq2R+sRAC8E/k/i/0ei18bKX1TVxwCi3y+IXh9VP0Tu/5XA3Yy8zVGI5BjwOHA78B3gKVU9E12SbNfZNkfv/wB4frc1dsIHgN8CtqP/n8/426zAH4vIURE5FL3W6dge8zOLJeO1KaZIjaYfROQi4CbgHar6tEhW0+aXZrw2uDar6hZwQET2ALcAL8m6LPo9+DaLyC8Aj6vqURF5VfxyxqWjaXPE1ar6qIi8ALhdRB4suNZLm8fsETwCXJ74/zLg0Z7q0gXfE5FLAKLfj0evj6IfRGQncyVwo6reHL086jbHqOpTwJ3M10f2iEhswCXbdbbN0fs/AWx0W9PWXA38ooh8F/hD5uGhDzDuNqOqj0a/H2eu8F9Jx2N7zIrgq8CLo4yDReBXgM/3XCeffB54U/T3m5jH0ePX/1GUbXAV8IPY5RwKMjf9bwBOqur1ibfG3OblyBNARC4ErmG+gHoH8NrosnSb4754LfAljYLIQ0FV36Oql6nqPubz9Uuq+g8YcZtF5C+IyHPiv4GfB47T9djue6HE8yLMa4BvMY+t/qu+6+OwXZ8CHgOeZW4hvJl5bPSLwLej30vRtcI8e+o7wDeA1b7r36C9f525+3s/cCz6ec3I2/zTwNejNh8H/m30+k8C9wAPAZ8FdkWvXxD9/1D0/k/23YaW7X8VcOvY2xy17b7o54FYTnU9tm1nsWEYxsQZc2jIMAzDqIApAsMwjIljisAwDGPimCIwDMOYOKYIDMMwJo4pAsMwjIljisAwDGPimCIwDMOYOP8fc54O1xRDDZsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ci = all_cis[0]\n",
    "for j, x in enumerate(ci):\n",
    "    plt.plot([j, j], x, color='grey')\n",
    "    plt.gca().set_ylabel(\"BMI\")\n",
    "mn = dx.BMXBMI.mean()\n",
    "plt.plot([0, 500], [mn, mn], color='red')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can calculate the fraction of the 500 simulated confidence intervals that did not cover the target value.  This is called the \"non-coverage probability\".  There are two ways for an interval to fail to cover the target -- either the upper limit of the interval can fall below the target, or the lower limit of the interval can fall above the target.  We calculate each of these below.  The sum of these two probabilities should be around 0.05, which is the allowed proportion of the time that a 95% confidence interval does not cover its target."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.036\n",
      "0.024\n"
     ]
    }
   ],
   "source": [
    "print(np.mean(ci[:, 1] < mn)) # Upper limit falls below the target\n",
    "print(np.mean(ci[:, 0] > mn)) # Lower limit falls above the target"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}