{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Logistic regression"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"remove-cell"
]
},
"source": [
"Think Bayes, Second Edition\n",
"\n",
"Copyright 2020 Allen B. Downey\n",
"\n",
"License: [Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0)](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
]
},
{
"cell_type": "code",
"execution_count": 91,
"metadata": {
"tags": [
"remove-cell"
]
},
"outputs": [],
"source": [
"# If we're running on Colab, install empiricaldist\n",
"# https://pypi.org/project/empiricaldist/\n",
"\n",
"import sys\n",
"IN_COLAB = 'google.colab' in sys.modules\n",
"\n",
"if IN_COLAB:\n",
" !pip install empiricaldist"
]
},
{
"cell_type": "code",
"execution_count": 92,
"metadata": {
"tags": [
"remove-cell"
]
},
"outputs": [],
"source": [
"# Get utils.py\n",
"\n",
"import os\n",
"\n",
"if not os.path.exists('utils.py'):\n",
" !wget https://github.com/AllenDowney/ThinkBayes2/raw/master/code/soln/utils.py"
]
},
{
"cell_type": "code",
"execution_count": 93,
"metadata": {
"tags": [
"remove-cell"
]
},
"outputs": [],
"source": [
"from utils import set_pyplot_params\n",
"set_pyplot_params()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Generational Changes\n",
"\n",
"As a second example of logistic regression, we'll use data from the [General Social Survey](https://gss.norc.org/) (GSS) to describe generational changes in support for legalization of marijuana.\n",
"\n",
"Since 1972 the GSS has surveyed a representative sample of adults in the U.S., asking about issues like \"national spending priorities, crime and punishment, intergroup relations, and confidence in institutions\".\n",
"\n",
"I have selected a subset of the GSS data, resampled it to correct for stratified sampling, and made the results available in an HDF file."
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"remove-cell"
]
},
"source": [
"The following cell downloads the data."
]
},
{
"cell_type": "code",
"execution_count": 61,
"metadata": {
"tags": [
"remove-cell"
]
},
"outputs": [],
"source": [
"# Load the data file\n",
"\n",
"import os\n",
"\n",
"datafile = 'gss_eda.hdf5'\n",
"if not os.path.exists(datafile):\n",
" !wget https://github.com/AllenDowney/ThinkBayes2/raw/master/data/gss_eda.hdf5"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use Pandas to load the data."
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(64814, 169)"
]
},
"execution_count": 62,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gss = pd.read_hdf(datafile, 'gss')\n",
"gss.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The result is a `DataFrame` with one row for each respondent and one column for each variable.\n",
"\n",
"The primary variable we'll explore is `grass`, which encodes each respondent's answer to this question ([details here](https://gssdataexplorer.norc.org/variables/285/vshow)):\n",
"\n",
"> \"Do you think the use of marijuana should be made legal or not?\"\n",
"\n",
"This question was asked during most years of the survey starting in 1973, so it provides a useful view of changes in attitudes over almost 50 years.\n",
"\n",
"Here are is the distributions of responses:"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"NaN 27268\n",
"2.0 25662\n",
"1.0 11884\n",
"Name: grass, dtype: int64"
]
},
"execution_count": 63,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gss['grass'].value_counts(dropna=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The value 1.0 represents \"yes\"; 2.0 represents \"no\"; `NaN` represents peope who were not asked the question and a small number of respondents who did not respond or said \"I don't know\".\n",
"\n",
"To explore generational changes in the responses, we will look at the level of support for legalization as a function of birth year, which is encoded in a variable called `cohort`. Here's a summary of this variable."
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"count 64586.000000\n",
"mean 1948.846069\n",
"std 21.262659\n",
"min 1883.000000\n",
"25% 1934.000000\n",
"50% 1951.000000\n",
"75% 1964.000000\n",
"max 2000.000000\n",
"Name: cohort, dtype: float64"
]
},
"execution_count": 64,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gss['cohort'].describe()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The oldest GSS respondent was born in 1883; the youngest was born in 2000.\n",
"\n",
"Before we analyze this data, I will select the subset of respondents with valid data for `grass` and `cohort`:"
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(37427, 169)"
]
},
"execution_count": 65,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"valid = gss.dropna(subset=['grass', 'cohort']).copy()\n",
"valid.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are about 37,000 respondents with the data we need.\n",
"\n",
"I'll recode the values of `grass` so `1` means yes and `0` means no."
]
},
{
"cell_type": "code",
"execution_count": 66,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0 25572\n",
"1.0 11855\n",
"Name: y, dtype: int64"
]
},
"execution_count": 66,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"valid['y'] = valid['grass'].replace(2, 0)\n",
"valid['y'].value_counts()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, for this problem, I'm going to represent the data in a different format. Rather than one row for each respondent, I am going to group the respondents by birth year and record the number of respondents in each group, `count`, and the number who support legalization, `sum`."
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plot_data(data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There is a strong relationship between birth year and support for legalization. People born before 1920 are the least likely to say \"yes\"; people born after 1990 are the most likely. \n",
"\n",
"There are substantial departures from the long-term trend for people born in the 1950s and late 1960s. If you want to conjecture about the causes, it might help to think about what was happening when each group turned 18. People born in 1950 turned 18 during [the counterculture of the 1960s](https://en.wikipedia.org/wiki/Counterculture_of_the_1960s). People born in the late 1960s turned 18 during the \"[Just Say No](https://en.wikipedia.org/wiki/Just_Say_No)\" era of the War on Drugs and the peak in the AIDS epidemic in the U.S."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Point estimates\n",
"\n",
"I'll use StatsModels again to generate point estimates for the slope and intercept of a logistic model.\n",
"\n",
"As we did with the previous problem, I'll center the values of the explanatory variable so the mean is 0."
]
},
{
"cell_type": "code",
"execution_count": 70,
"metadata": {},
"outputs": [],
"source": [
"offset = valid['cohort'].mean()\n",
"valid['x'] = valid['cohort'] - offset"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here are the results from StatsModels."
]
},
{
"cell_type": "code",
"execution_count": 71,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"Intercept -0.849094\n",
"x 0.031691\n",
"dtype: float64"
]
},
"execution_count": 71,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import statsmodels.formula.api as smf\n",
"\n",
"formula = 'y ~ x'\n",
"results = smf.logit(formula, data=valid).fit(disp=0)\n",
"results.params"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To visualize the results, I'll use these parameters to estimate the probability of support in each cohort."
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {},
"outputs": [],
"source": [
"inter = results.params['Intercept']\n",
"slope = results.params['x']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"I'll shift the birth years in `data` by `offset`."
]
},
{
"cell_type": "code",
"execution_count": 73,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"probs.plot(label='Logistic model', color='C1')\n",
"\n",
"plot_data(data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With these parameters, the model captures the long term trend in the data."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Computing likelihoods\n",
"\n",
"Before we do the Bayesian update, let's compute the probability of the data with the estimated parameters.\n",
"\n",
"From the data, we know how many people there are in each group and how many of them support legalization. From the model, we have an estimate for the probability of support in each group.\n",
"\n",
"So we can use the binomial distribution to compute the probability of the data given the estimated probabilities."
]
},
{
"cell_type": "code",
"execution_count": 76,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(116,)"
]
},
"execution_count": 76,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from scipy.stats import binom\n",
"\n",
"ks = data['sum']\n",
"ns = data['count']\n",
"likes = binom.pmf(ks, ns, probs)\n",
"likes.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For each group `likes` contains the probability of the outcome, `k`, given the group size, `n`, and the estimated probability, `p`.\n",
"\n",
"The likelihood of the data is the product of these likelihoods:"
]
},
{
"cell_type": "code",
"execution_count": 77,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.5665006815705933e-209"
]
},
"execution_count": 77,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"likes.prod()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This likelihood is very small, for two reasons:\n",
"\n",
"* The dataset is large, which means that there are many possible outcomes, so the probability of *any* particular outcome is small.\n",
"\n",
"* The data deviate substantially from the model, so the probability of *this* particular outcome is small.\n",
"\n",
"In theory, it's not a problem if the likelihood of the data is small. We might not get a model that fits the data perfectly, but we'll get the parameters that come as close as possible.\n",
"\n",
"However, in practice small likelihoods can be problematic. With floating-point numbers, the smallest positive number we can represent is about `1e-1021`."
]
},
{
"cell_type": "code",
"execution_count": 78,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-1021"
]
},
"execution_count": 78,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import sys\n",
"sys.float_info.min_exp"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Any number smaller than that \"underflows\"; that is, it gets rounded down to 0. When that happens, we lose the ability to distinguish between parameters that make the model fit the data or not. In the worst case, if all likelihoods underflow, all probabilities in the posterior distribution would be 0.\n",
"\n",
"In this example, the likelihoods are big enough that we can still do a Bayesian update, so we'll do that next.\n",
"Then I will demonstrate a trick we can use to avoid underflow: computing likelihoods under a log transformation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The update\n",
"\n",
"I'll use uniform priors for the parameters, with locations centered around the point estimates."
]
},
{
"cell_type": "code",
"execution_count": 79,
"metadata": {},
"outputs": [],
"source": [
"qs = np.linspace(-0.95, -0.75, num=51)\n",
"prior_inter = make_uniform(qs, 'Intercept')"
]
},
{
"cell_type": "code",
"execution_count": 80,
"metadata": {},
"outputs": [],
"source": [
"qs = np.linspace(0.025, 0.035, num=51)\n",
"prior_slope = make_uniform(qs, 'Slope')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"I'll make a joint prior."
]
},
{
"cell_type": "code",
"execution_count": 81,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
Intercept
\n",
"
-0.950
\n",
"
-0.946
\n",
"
-0.942
\n",
"
-0.938
\n",
"
-0.934
\n",
"
-0.930
\n",
"
-0.926
\n",
"
-0.922
\n",
"
-0.918
\n",
"
-0.914
\n",
"
...
\n",
"
-0.786
\n",
"
-0.782
\n",
"
-0.778
\n",
"
-0.774
\n",
"
-0.770
\n",
"
-0.766
\n",
"
-0.762
\n",
"
-0.758
\n",
"
-0.754
\n",
"
-0.750
\n",
"
\n",
"
\n",
"
Slope
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0.0250
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
...
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
\n",
"
\n",
"
0.0252
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
...
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
\n",
"
\n",
"
0.0254
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
...
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
\n",
"
\n",
"
0.0256
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
...
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
\n",
"
\n",
"
0.0258
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
...
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
0.000384
\n",
"
\n",
" \n",
"
\n",
"
5 rows × 51 columns
\n",
"
"
],
"text/plain": [
"Intercept -0.950 -0.946 -0.942 -0.938 -0.934 -0.930 \\\n",
"Slope \n",
"0.0250 0.000384 0.000384 0.000384 0.000384 0.000384 0.000384 \n",
"0.0252 0.000384 0.000384 0.000384 0.000384 0.000384 0.000384 \n",
"0.0254 0.000384 0.000384 0.000384 0.000384 0.000384 0.000384 \n",
"0.0256 0.000384 0.000384 0.000384 0.000384 0.000384 0.000384 \n",
"0.0258 0.000384 0.000384 0.000384 0.000384 0.000384 0.000384 \n",
"\n",
"Intercept -0.926 -0.922 -0.918 -0.914 ... -0.786 -0.782 \\\n",
"Slope ... \n",
"0.0250 0.000384 0.000384 0.000384 0.000384 ... 0.000384 0.000384 \n",
"0.0252 0.000384 0.000384 0.000384 0.000384 ... 0.000384 0.000384 \n",
"0.0254 0.000384 0.000384 0.000384 0.000384 ... 0.000384 0.000384 \n",
"0.0256 0.000384 0.000384 0.000384 0.000384 ... 0.000384 0.000384 \n",
"0.0258 0.000384 0.000384 0.000384 0.000384 ... 0.000384 0.000384 \n",
"\n",
"Intercept -0.778 -0.774 -0.770 -0.766 -0.762 -0.758 \\\n",
"Slope \n",
"0.0250 0.000384 0.000384 0.000384 0.000384 0.000384 0.000384 \n",
"0.0252 0.000384 0.000384 0.000384 0.000384 0.000384 0.000384 \n",
"0.0254 0.000384 0.000384 0.000384 0.000384 0.000384 0.000384 \n",
"0.0256 0.000384 0.000384 0.000384 0.000384 0.000384 0.000384 \n",
"0.0258 0.000384 0.000384 0.000384 0.000384 0.000384 0.000384 \n",
"\n",
"Intercept -0.754 -0.750 \n",
"Slope \n",
"0.0250 0.000384 0.000384 \n",
"0.0252 0.000384 0.000384 \n",
"0.0254 0.000384 0.000384 \n",
"0.0256 0.000384 0.000384 \n",
"0.0258 0.000384 0.000384 \n",
"\n",
"[5 rows x 51 columns]"
]
},
"execution_count": 81,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"joint = make_joint(prior_inter, prior_slope)\n",
"joint.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And stack it into a `Pmf` with a two-column index."
]
},
{
"cell_type": "code",
"execution_count": 82,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Slope Intercept\n",
"0.025 -0.950 0.000384\n",
" -0.946 0.000384\n",
" -0.942 0.000384\n",
" -0.938 0.000384\n",
" -0.934 0.000384\n",
"dtype: float64"
]
},
"execution_count": 82,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"joint_pmf = Pmf(joint.stack())\n",
"joint_pmf.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's the update, using the binomial distribution to compute the likelihood of the data in each group."
]
},
{
"cell_type": "code",
"execution_count": 83,
"metadata": {},
"outputs": [],
"source": [
"likelihood = joint_pmf.copy()\n",
"\n",
"xs = data['x']\n",
"ks = data['sum']\n",
"ns = data['count']\n",
"\n",
"for slope, inter in joint_pmf.index:\n",
" ps = expit(inter + slope * xs)\n",
" likes = binom.pmf(ks, ns, ps)\n",
" likelihood[slope, inter] = likes.prod()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Again, the likelihoods are small."
]
},
{
"cell_type": "code",
"execution_count": 84,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"8.656375286337876e-208"
]
},
"execution_count": 84,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"likelihood.sum()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But we can do the update in the usual way."
]
},
{
"cell_type": "code",
"execution_count": 85,
"metadata": {
"tags": [
"hide-output"
]
},
"outputs": [
{
"data": {
"text/plain": [
"3.328095073563197e-211"
]
},
"execution_count": 85,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior_pmf = joint_pmf * likelihood\n",
"posterior_pmf.normalize()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And there are enough non-zero elements to get a useful posterior distribution.\n",
"\n",
"Here's what it looks like."
]
},
{
"cell_type": "code",
"execution_count": 86,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxi0lEQVR4nO3deXwV9fX/8dchYQ1LWAKEfQcVFBABC+6o4AbuWFuXutS2Wq2tVatttda617pVv1qtS1XUCoqIRVRcQNlBdiSyQ4AQ1iQQCDm/P+7E3yXcbHAvmcD7+Xjcx713Pp+ZOTOZ5J2ZO3fG3B0REZGwqVbZBYiIiMSigBIRkVBSQImISCgpoEREJJQUUCIiEkoKKBERCSUFlFQ6M8sxsw6VXUeiJXI5zexkM1sd9X6+mZ0cp2lfbmYfR713M+sUj2kH0zssfv5ScQooSRgz+9zMri2rn7vXdfel5ZxmXP84lnOe7YL5Jh/IdCqynAfK3Y9y989L61Pe5XL31939jHjUFWubOJjrRaoWBZRIgh1osB3o+FV13iIKKDkozOw6M8sws01mNtrMWkS1/bBXZGYvm9kzZvahmW03sylm1jFo+zIY5dvgsNClMeZzlZlNMrOnzGyrmS0ys9Oi2lsE898U1HNdVFtfM5tuZtvMbL2Z/T1oKprvlmC+xwf9f2ZmC81ss5mNM7O2xZbpV2a2BFgSYzkbmNmrZpZlZivM7G4zq1ZsGR43s03APTGWs3awrjab2QLguGLty81sUEWXK9a8g2ETi5VwlpktNbONZvZIVO33mNl/our4YS/NzO4HTgCeDub39H6sl4lm9miw3MvMbEjxdSOHDgWUJJyZnQo8AFwCpAMrgBGljHIZcC/QEMgA7gdw9xOD9mOCw0JvlTB+P2Ap0AT4MzDSzBoFbW8Cq4EWwEXA36IC7AngCXevD3QE3g6GF803NZjvN2Y2DPgDcAGQBnwVTDvasKCWI2PU+BTQAOgAnARcAVwdYxmaFi1/MX8OauwInAlcGXtVVGy5yjlvgPOBPkBvYCjws1LmD4C730VkPd0YzO/GGN3Ks14WE/nZPgy8aGZW1rylalJAycFwOfCSu89093zgTuB4M2tXQv+R7j7V3QuA14GeFZzfBuAf7r47CLHFwNlm1hoYCNzu7jvdfTbwL+CnwXi7gU5m1sTdc9x9cinz+DnwgLsvDOr8G9Azei8qaN/k7juiRzSzJOBS4E533+7uy4HHouoAWOvuT7l7QfHxA5cA9wfTXwU8WUqtFVmu8swb4KFg3iuBfxD5p+KAlHO9rHD3F9x9D/AKkX94mh3ovCWcFFByMLQgstcEgLvnANlAyxL6r4t6nQfUreD81vjeV0FeEdTQAtjk7tuLtRXVcQ3QBVhkZtPM7JxS5tEWeMLMtpjZFmATYOy9TKtKGLcJUIOodVKsjtLGLdKiWJ8VJXWkYstVnnkX71O0fg9UedbLD9uGu+cFLyu6fUgVoYCSg2EtkT/oAJhZCtAYWJOg+bUsdtinTVDDWqCRmdUr1rYGwN2XuPtlRA5tPQT8N6g11iX/VwE/d/fUqEdtd/86qk9JtwrYSGSvJnpv64c6yhi3SCbQutj4MVVwucozb2LMe23wOheoE9XWvALTLs96kcOIAkoOhjeAq82sp5nVJHI4bEpwCKei1hP5fKI0TYFfm1l1M7sYOAIYGxwK+xp4wMxqmdnRRPYuXgcws5+YWZq7FwJbgmntAbKAwmLzfQ6408yOCsZtEMyrTMHhqbeB+82sXnBY8FbgP6WPuZe3g/k3NLNWwE0ldazgcpXXbcG8WwM3A0WfB84GTjSzNmbWgMjh3Ggl/vzitF7kEKKAkkRzd/8U+CPwLpH//DsCw/dzevcArwSH1i4poc8UoDOR/8jvBy5y9+yg7TKgHZH/+EcBf3b38UHbYGC+meUQObFgePBZVV4wnUnBfPu7+ygieyMjzGwbMA+oyBllNxHZ21gKTCQS4i9VYPx7iRz+WgZ8DLxWSt9yL1cF5v8+MINIIH0IvAgQrMu3gDlB+5hi4z0BXBSchRfrc7MDXS9yCDHdsFASxcxmAn9x9/cO4jyvAq5194EHa54ikhjag5KECA59HQHMquxaRKRqUkBJ3JnZQ0QOO93u7qWdXSYiUiId4hMRkVDSHpSIiITSYXEhyCZNmni7du0quwwRkUPGjBkzNrp7WiLncVgEVLt27Zg+fXpllyEicsgws4R/vqxDfCIiEkoKKBERCSUFlIiIhJICSkREQkkBJSIioaSAEhGRUFJAiYhIKCmgREQklBRQIiISSgooEREJJQWUiIiEkgJKRERCSQElIiKhpIASEZFQUkCJiEgoKaBERCSUFFAiIhJKCQ0oMxtsZovNLMPM7ojRbmb2ZNA+x8x6B8NrmdlUM/vWzOab2b0xxv2dmbmZNUnkMoiISOVIWECZWRLwDDAEOBK4zMyOLNZtCNA5eFwPPBsMzwdOdfdjgJ7AYDPrHzXt1sDpwMpE1S8iIpUrkXtQfYEMd1/q7ruAEcDQYn2GAq96xGQg1czSg/c5QZ/qwcOjxnsc+H2xYSIicghJZEC1BFZFvV8dDCtXHzNLMrPZwAZgvLtPCYafB6xx929Lm7mZXW9m081selZW1gEtiIiIHHyJDCiLMaz4Hk+Jfdx9j7v3BFoBfc2su5nVAe4C/lTWzN39eXfv4+590tLSKla5iIhUukQG1GqgddT7VsDaivZx9y3A58BgoCPQHvjWzJYH/WeaWfM41i0iIiGQyICaBnQ2s/ZmVgMYDowu1mc0cEVwNl9/YKu7Z5pZmpmlAphZbWAQsMjd57p7U3dv5+7tiARcb3dfl8DlEBGRSpCcqAm7e4GZ3QiMA5KAl9x9vpndELQ/B4wFzgIygDzg6mD0dOCV4EzAasDb7j4mUbWKiEj4mPuhfyJcnz59fPr06ZVdhojIIcPMZrh7n0TOQ1eSEBGRUFJAiYhIKCmgREQklBRQIiISSgooEREJJQWUiIiEkgJKRERCSQElIiKhpIASEZFQUkCJiEgoKaBERCSUFFAiIhJKCigREQklBZSIiISSAkpEREJJASUiIqGkgBIRkVBSQImISCgpoEREJJQUUCIiEkrJlV2AyKGksLCQNUsyWbVoLXnbd7AjZyc7c/PZmbOTHTk72JmbT2GhUze1DnUb1qVuagr1GqZQt2EKdVNTaNo2jYZNG1T2YoiEggJKZD8VFhayNmMd381YynfTv2fJzKVkzFxG/cZ1aXtUa1Ia1KF2Si1q1a1FrZSaNEhrQPP2tQDI3ZrH9uztrM1YR87WXHI255KzOYd1yzaQXCOZ9j3a0L5HWzocHXm0OaIlNWrVqOQlFjm4FFAiFbR+RRajnhzLx698Tkr92nQ+tgOde3fkx3+4gM69O1C/cb39nra7s3HNJpbNXcnSOSuYMf5b3nlsNGsz1nHk8V04+dIBnHBh/wOah0hVYe5e2TUkXJ8+fXz69OmVXYZUcQsmf8fIf4xh5idzOfOqUxh642Cat2t6UOadvyOfaf+bzYQRk5g+bjbdB3bjlOED+dHQ46hTr/ZBqUEkmpnNcPc+CZ2HAkqkZO7O5DEzGPHgKDZlbub8X5/NmT87hZT6dSqtprztO/hm9HQmjJjI3K8WcuwZx3DuDWfQ85TumFml1SWHFwVUnCigZH9s3rCVh654io2rs7ninksYMKwvSclJlV3WXrZt2s4Xb3/D+09/RLWkalxw89mc+uOB+rxKEk4BFScKKKmojFnLuOeCRzj1xwO58t5L9yuYdu/eQ2bWVtas28KadVtYnbmFNes2sy1nJ7VqVqd2rep7PdeqWZ3mafU5olNzOrRpQnIF5unuzBg/h5FPfMiSGUs5+/pBnPuLM2mc3rDCdYuUhwIqThRQUhGfvzWJp296kZuevpaTLvlRhcb9fkUWYz6dy6Tp35OVnUOTRnVplZ5Ky+aRR6vmDalfrxb5+QXsyN/Nzp27f3jemb+bNeu2sDAjk/Ubt9OxbRpHdGoePNJplZ5arkN4Kxet4b0nxzJhxCSOP68PV957Kc3apu3v6hCJSQEVJwooKY/CwkJe/uMIPntjIveMuo1OPduXa7yc3Hw+mbiQMZ/OJXtzLkNO6c6ggV1pnd6I6tX375Bgbl4+i5euZ2HGOhZmrGPBkkxqVE/mrFOO4syTjqRZk/plTmPbpu2MemIs7z/zPy685Rwu/t25OvQncaOAihMFlJTHa/e+w/Tx33LvqNtITSv7y7I7du7i7TEzeXvMDHp1b805p/XguKPbkpQU/wu0uDvzl2Tyvwnz+ezrxXTt2Iwhp3TnxL6dqFWzeqnjrlu+gedufZmVi9Zyx2s30eXYjnGvTw4/Cqg4UUBJWdYt38CvjruDZ2c8RNM2ZR8OW7Akk7sfGU2Pbi247rKBtDqIn/Xk5+9m4vTvGTthHguWrOPk/l24fNhxZdYwYcQk/nnzSwz79VkMv2MYSUnhOuFDqhYFVJwooKQs9170KJ16tufyuy8ss+9XUzN46Nlx3PHLwQw8ruy9kfyCAqauWc0XK5bz5YplZGzetFd79KdKzVLqckzzdHo1T6dn83R6NG1G7eol7yFlZW/ng0/n8u7YWQw55SiuvKg/9VJqldh/w6qNPHL1M+zauZs7Xr2J9A7NyqxfJBYFVJwooKQ0Mz+Zw+M//z9enP94mZ/RvDt2Jq+NmsqDtw+jW6fmJfZbunnTD4E0Y+1aujZpwklt23Ni23YcldaUasHJDtG/fe7Omu3bmLUuk2/XZTJrXSbfZW+kQ8NG9GyezpkdOzOgdZuYJ0ps2pLLC29OYuK0DK66+HiGnnEMySUcaiwsLGTUE2MZ8eAobnrmWk686PiyV5JIMQqoOFFASUncnet63MrVf72MAcP6ltr31XcnM+6LBTxy1wW0aJYas8+O3bu59eOPmLVuLae068CJbdsxoHUb6tcsea+mNPkFBczP2sDMzLW8u3A+u/bs4Zpex3LhEUdRM3nfK5VlLM/iqZcnkL05h9t/cSY9urUscdqLp3/PXy95jJOHD+Rn91+mL/lKhSig4kQBJSVZsWAVd539AK8tfabUP9Cbt+by41//m9cev4omjerG7JNfUMB1H7xHkzp1eOC0M2IGyIFwd6atXcOz06fwXfZGrut9HMO796BWcvV9+n05ZQmPvfAJQ08/hqsuPr7EEze2ZW/n96f/hX5n9ebqv14W13rl0HYwAiqh94Mys8FmttjMMszsjhjtZmZPBu1zzKx3MLyWmU01s2/NbL6Z3Rs1ziNmtijoP8rMUhO5DHJom/XZPHqdWvYlgt4eM5PTftS1xHDatWcPvxz7AQ1q1eTh0weXGU6FXsj23dtZnbeG+VsX8PXGyXyUOY4vsr4ic8c6Yv3jaGb0bdmKfw+9kH+ePZRvVq/kpJdf5IWZ08gvKNir30n9u/DSo1cwe8Fq7nzoPXLz8mPWUb9xPR4cdzcTR03hjb+NLLVmkYMtYXtQZpYEfAecDqwGpgGXufuCqD5nATcBZwH9gCfcvZ9F/lqkuHuOmVUHJgI3u/tkMzsD+MzdC8zsIQB3v720WrQHJSW558JHOOGC/px2+Qkl9snJzefSX/2L5x+8nJbNU/dpLygs5KaPxlDohTw95Fyql3B23JZdWxi99kNmbp7FtoJt1KpWkwbVG/zwqF+9PtsLtvPd9u/YVbibLvU6Rx51O9OmTmuSq+0beguzNvDYN5NYs30bfz/zLI5osvcZiAUFe/j7vz5l3uK1PHTn+aSXcK+p7MzN3HrSnxj6y8FccMvZpawxkYiDsQeVyNtt9AUy3H0pgJmNAIYCC6L6DAVe9UhKTjazVDNLd/dMICfoUz14OIC7fxw1/mTgogQugxzCCgsLmfPFAm586ppS+40aN5t+vdrFDKdCd3778UfsLNjNc2cPjRlOO/fsZPTaMXy+4QtOSBvInUf8nkY1GlK9Wsln52XnZ/NdTgZLti/hq6xJbMzfSL/GxzG0xXk0rtnoh35HpDXlhXOH8e7C+fx05Dtcf+xxXNe7zw97hMnJSdz289P579hZ3HDnG9x323kcHeNzqcbpDXl4/J/47cl/plZKTc66blCp60TkYEhkQLUEVkW9X01kL6msPi2BzGAPbAbQCXjG3afEmMfPgLdizdzMrgeuB2jTps3+1C+HuJwtuezZvYcmLRqV2m/KrGVceVH/mG1z169jzvp1fHT5FSUe1ntj5Vts272Nv/a4l0Y1Sp9XkcY1G3N8zcYc3zjyK7Nt93bGrfuYe+bfx1Xtf8qxDXv/0NfMuOjI7vRv1Zpfjv2A7B153DHgxB9Cysy4+OzetEpP5Q8PvccLD/0k5p5Us7ZpPDT+j9wy4G6OOL4L7bvr90YqVyI/g4p1UL/48cQS+7j7HnfvCbQC+ppZ971GNLsLKABejzVzd3/e3fu4e5+0NF2HTPZVr2Fddu8qIH9H7M9nihTsKaRmCVdrWJy9kd7pLfY5UaHIspzlzN7yLdd1uKbc4RRL/er1uLj1hdzc5UbeWDGC/6x4g92Fu/fq06p+A14ddiFfLF/GP6dP3Wcax/fuwI+H9eWex8dQULAn5nxadkrninsv5bFr/snuXbtj9hE5WBIZUKuB1lHvWwFrK9rH3bcAnwODi4aZ2ZXAOcDlfjichigJYWY0Tk8le+3mUvvt2VNIUrXYvyrfZWfTpXHjmG2FXshrK17nolYXkJIcn/tHdarbkb90/zObdm3mvgUPsH7n+r3aU2vV5pVhF/HO/Ln8Z87sfcYffm4f6qXU4oU3J5U4j3N+fjoNm6Xy4p1vxKVmkf2VyICaBnQ2s/ZmVgMYDowu1mc0cEVwNl9/YKu7Z5pZWtHZeWZWGxgELAreDwZuB85z97wE1i+HgcYtG7FxzaZS++wpdJKSYp/l9132Rro0bhKzbeLGrwEY2KRiV0QvS0pyCjd1+iUnpg3kvgV/45vsvY9+N6tbl1fPv4h/TpvC+4sX7tVWrZpx101D+PjLBUyZtSzm9M2M2/79K7787zd884FOLpLKk7CAcvcC4EZgHLAQeNvd55vZDWZ2Q9BtLLAUyABeAH4ZDE8HJpjZHCJBN97dxwRtTwP1gPFmNtvMnkvUMsihL611E1Z/l1lqn2oWubdTLKu2baVFvdhXFh+z9kMubn0h1Sz+v2ZmxqBmp3Jb11t5d/VIpm6atld7mwap/HvYhfz1ywksyc7eq61hgzr88eaz+NvT/2PHzl0xp1+/cT3+8MYt/P2659iRsyPu9YuUR0K/B+XuY929i7t3dPf7g2HPuftzwWt3918F7T3cfXowfI6793L3o929u7v/JWqandy9tbv3DB43xJ67SNlOGT6AD58fH/N7R0V692jD1zOWxmzr06IlE1euiNnWJqUNmTvXxaXOkrRNacslrS/i0/UT9mnr2rgJP+5xDK/Pnb1PW+/ubTiic3M+/nLhPm1Fug/oRrd+nfji7W/iWbJIuSU0oETCrt/ZvcnZnMPCKUtK7DNoQDc+mbgoZoid17UboxfH/iM/oPHxfL0x8X/ce6X2ZPWONWTlZ+3TdulRPXh/8SLydu97wsOFQ3rz7kezSg3ns687nbH/+iSu9YqUlwJKDmtJSUmcdvmJfPlOyUHSuX1TkpOrsWDJvntDP2rVhnU5OSzdvO/nWD0adGfdznVs2LlvcMRT9WrV6d+oLxOzvt6nrUW9+vRp0YIx3y3ap63P0W0oLCxk1vxV+7QVOW5wT7JWZbNs3sq41ixSHgooOewNOL8vk0ZNobCwMGa7mTFoQDc+/nLBPm1J1apxdpeujFy4b1tytWT6NurLl1lfxb3m4k5IG8jEjbHPzLus+zG8OW/OPsPNjAuG9GLkR7NKnG5SchJnXn0KY1/QXpQcfAooOex1OLotjdIbMu7f+36OU+Tc04/ms68XM3fRmn3arjy6F2/Nnxvzs6jBzc9gUvbXjF/3aVxrLs7dsRJOxujapAnrc3JitvXp0ZaMFRtLnfZxg3uxaGrGAdcoUlEKKDnsmRk3PnUN/777TbZt2h6zT9PG9bj9F2dyz+Nj2Lp977Pa2qam8vSQc/jNuA/5LnvvP/ZNa6XxhyPuYPz6T3hvzehSP+85EF9kfcmJaQNjtq3Zto2W9WOfaZiUVK3EPccijdJT2ZRZ+nfFRBJBASUCdO7dgYHn9+PlP8a8chYAA4/ryKk/6sr9T31EYeHeQdOvVWvuOuFkrhk9iqzc3L3a0mo24Q9H3MH0TTN4Y+UICr30QKio/D35TN00jYFNBsRsX7t9G+n16sVsq1bN9lmW4hqnN2Tzui0JC1eRkiigRAJX/XU4E0dOJqOEL7AC/PzyE9iWs5MRo6ft0zas25FcfGR3rv1g1D5nzaXWaMCdR9zOstzlPJXxT7LySz+sVhFTN02jU91ONKrRMGb7mu3baVnCd7WSkqqxZ0/pgVmjVg1q1qnJ9k2xDxOKJIoCSiRQv1E9rn3wJ/x1+ONsXJMds09ychL33noOb42Zwejx+554cFPf/nRrksbw/45g7oa9L0OUklyH33f7LW3qtObP8+7lscX/YEr2VHYVVvyad4VeyMJti3hx2cu8sXIEpzc7LWa/gsJCPlv2PV0bx74e5YaN26lTu/Tb3APUrFODHTk7K1ynyIHQHXVFinnr4fcZ+69PePSze0hrFfs6e6vWbuZ397/LSf068/PLT9jrjrXuzn8XzueRr79icMfO/Pb4gTSotfct3/P35DNj8ywmbpzE8twVHNeoDyc0+REd63Ys9eaJq/NW83X2ZL7Jnkzd5Loc37g//Rv3K3Hv6fHJk5iVmcnLwy6kWozp3v3I+xx9RCsuOefYEue5JWsrV3e9mZHZ/9Zt4eUHuuV7nCigpKLeeXQ0Hzz3MY9+9meatom997F1+w7ufmQ0tWtV58+3nE1KnZp7tW/ZuYNHvp7IJ0u/5/YBJ3B+tyNj/oHPzt/E19nfMHXTNFbmraJB9fqkJKVQJzmFlKQ6pCTXoXZSHZbkLCGnIIfjG/fn+Mb9aV2nVanL8PLsmbwwcxqjLr2cpin73gl45dpN/PKuN3nn2euoXavkvaipH83iv3//gIfH/6nU+cnhRQEVJwoo2R/vPj6G957+iMcm3FNiSBUU7OEfL37G7AWreejO82Pe1PDb9ev444RPqJ2czF9OGUTXEi4uC7CrcDd5BbnkFOSStyeP3IJccvfkkVeQS+s6relar0uZ1/Zzdx75eiLjvl/CK8MupFX92HfRfejZj2nSMIVrhsc+uaLIa395h/y8fK598Cel9pPDS1W/o65IlXbhb84B4LZBf+HeUb+n3VGt9+mTnJzE735+OiP/N4tf3PUGv7n2NE7u32WvPaVjmjVn1CU/5s15c7j83bc5qV17zu3SjeNbtd7nJoc1qlWnRo1UUmuk7lfNizZm8Y8pX7M+J4d3Lh5Oo9qxb/OxMCOTCd8s5s0y7iYMMPerhZx7wxn7VY/IgdAelEgZxr08gRd+/xpX3XcZZ18/qMTPYeYsWsPDz35Mi2YNuPW6QTRP2/fMuU078hi5cAHjvl/Ckk3ZnNi2HWd26MxJ7dpTt0bZJyvE4u5MWbOa/5sxjYUbN3DVMb254phe1Kke+yaKH3+1kCdf+ozbf3EmJ/TtVOq0vx49jedve43/m/0INWvXLLWvHF50iC9OFFByoFYuWsP9wx+nZZd0bn3+BuqmpsTst3v3Hl5/fypvfTCDk/t3Zvi5fWhbwokWWbm5jF+awcffZzAzcy3HtWxF7/QWNE1JoVlKXdJSUmhaJ4WGtWtTzQx3Z+OOPJZv2cyKLVtYsXULK7ZsYXH2RgoKC7m+dx+GdTuyxFvPFxY6L7w5kU8mLuSB28+nU7vS7zSduzWX63r8lttfu4ljTjqqYitMDnkKqDhRQEk87Nq5i+dve40pH87gjv/czFE/6lpi3y3b8hj5v9mM+t9sjuqSzmVDj+Pobi1L3Pvalp/P58uXsmBjFlm5uWzIzSErN5f1ubnk7d5F49p12L4rnxpJSbRt0JB2qam0TU2lbYOGtE9NpUez5jHP0iuSt2MX9/7jQ3Lz8rnvd+fRsEHZd/h98pcvsKdgD795Xne0kX0poOJEASXx9PX703j85//HsBuHcMnvz6N6jdiH0gB25u/mownzGfHBdFLr1+ay845jYN9OJCeV/yuI+QUFbMzLo26NGvucrl4eK1Znc/ejo+netSW3Xnsa1asnlTnOnC8X8MDlT/DC3L+XuLcohzcFVJwooCTeNqzayD9ueJ7Vi9dyzd9+zIkXH1/qd4T27Cnkq6kZvPXBdFZlbqZ/7/b079WBvse0pX692nGvL2/HLmbNX8WEbxYzeeYyrrl0AMPOPKZc32Oa+9VC7rvkMX730q/oO6RX3GuTQ4MCKk4UUJIoMz+dy79uf41qSdW47uGfluuzmrXrtzB51jImz1zG7AWr6dimCf16t6d/r/Z0ad+MatUq/mVYd2fpyo1MmbWMKbOXs2BJJkd0ak7/3h04d1AP6qWUvefl7rz31Ee88beR/P6VGznuzJ4VrkMOHwqoOFFASSIVFhby+YhJ/PvuN2nXvQ3XPHB5zFPSY8nfVcCchav5ZmYksDI3bKVxwxSaNqlHWqN6NGtSj7TG9WjauC6YkbdjF3l5+eTu2EVu3i5yd+SzPWcn3y5cQ43qSfTr1Z5+PdvTu3vrcl3CqMjOvHwev/45ls9fxT3v3kZ6h2b7uzrkMKGAihMFlBwMu/J388E/xzHiwVH0GdyTs64dRPeB3Sp0eaD8XQVs3JTDhuztZGVvZ/3GyHNWduRCrXXq1CCldg3q1K5BSp2a1Kldg7p1anJUl3Rapce+3FFZ1n6/jnsvfJQOx7Tl5mevp1YdnU4uZVNAxYkCSg6mnC25fPj8J4x/9XN25+9myLWDOPOqk2nYLLWyS9tL/o583nvqf7zz6Pv85E8XM/RXg3WtPSk3BVScKKCkMrg7i6ZmMPb58UwcNZVeg3pw9nWD6HVaD6pVq5wbCbg7y+et5It3vuHjVz6ny7EduOaBy2ndtWWl1CNVlwIqThRQUtlyt+by2RsT+fCFT8heu5keJ3Sj+4AjOGpgNzr1bEdSctmnfu8vd2f5/FV88fbXfPnfyezM3cmJFx3PqT8eSJdjOyZsvnJoU0DFiQJKwiRz2XrmT1rM/EmLmDdpEeuXZ9G1bye6D+hGh6Pb0rhlI9JaNaZR89QKB9eePXvYlLmF9cs3sG55FisXrmbSe1PZkRMJpRMvPp5ufTtV2h6cHDoUUHGigJIw2745hwXffMe8iYtYtXgNG1dns3HNJrZmbaN+k/qktWpEk5aNqFV0qnjwMZGZYWYUFhb+EEpZq7Kp16guzdql0axdU1p0aEa/c45VKEnc6WrmIoeBeg3r0u+s3vQ7q/dew/cU7GHTui1krc4me80m8nfs4od/KB2i/7lslN6Q5u3SaNqmCTVKubeTSFWigBIJqaTkJNJaNS7xrr4ihzrt84uISCiVK6DMrIuZfWpm84L3R5vZ3YktTUREDmfl3YN6AbgT2A3g7nOA4YkqSkREpLwBVcfdpxYbVhDvYkRERIqUN6A2mllHwAHM7CIgM2FViYjIYa+8Z/H9Cnge6GZma4BlwOUJq0pERA575Qood18KDDKzFKCau29PbFkiInK4K+9ZfI3N7EngK+BzM3vCzPTlDBERSZjyfgY1AsgCLgQuCl6/VdZIZjbYzBabWYaZ3RGj3czsyaB9jpn1DobXMrOpZvatmc03s3ujxmlkZuPNbEnwvH83wRERkVArb0A1cvf73H1Z8PgrkFraCGaWBDwDDAGOBC4zsyOLdRsCdA4e1wPPBsPzgVPd/RigJzDYzPoHbXcAn7p7Z+DT4L2IiBxiyhtQE8xsuJlVCx6XAB+WMU5fIMPdl7r7LiJ7YUOL9RkKvOoRk4FUM0sP3ucEfaoHD48a55Xg9SvAsHIug4iIVCHlDaifA28Au4LHCOBWM9tuZttKGKclsCrq/epgWLn6mFmSmc0GNgDj3X1K0KeZu2cCBM9NY83czK43s+lmNj0rK6t8SykiIqFRroBy93ruXs3dk4NHtWBYPXevX8Jose4dXfzeHiX2cfc97t4TaAX0NbPu5ak1qubn3b2Pu/dJS0uryKgiIhIC5b6auZmdB5wYvP3c3ceUMcpqoHXU+1bA2or2cfctZvY5MBiYB6wPDgNmmlk6kT0sERE5xJT3NPMHgZuBBcHj5mBYaaYBnc2svZnVIHLtvtHF+owGrgjO5usPbA2CJ83MUoN51wYGAYuixrkyeH0l8H55lkFERKqW8u5BnQX0dPdCADN7BZhFKWfQuXuBmd0IjAOSgJfcfb6Z3RC0PweMDaadAeQBVwejpwOvBGcCVgPejtpjexB428yuAVYCF5d3YUVEpOqoyA0LU4FNwesG5RnB3ccSCaHoYc9FvXYil1EqPt4coFcJ08wGTitXxSIiUmWVN6AeAGaZ2QQiJzacSOT2GyIiIglR3mvxvRmcqHAckYC63d3XJbIwERE5vJUaUEWXHoqyOnhuYWYt3H1mYsoSEZHDXVl7UI8R+V5S0feVin+P6dS4VyQiIkLZAXU7sKroyg1mdiWRC8YuB+5JaGUiInJYK+t7UM8RuXArZnYikZMlXgG2ErmBoYiISEKUtQeV5O5Fp5ZfCjzv7u8C7wbXyRMREUmIsvagksysKMROAz6LaqvId6hEREQqpKyQeRP4wsw2AjuI3FEXM+tE5DCfiIhIQpQaUO5+v5l9SuTSQx8HV36AyJ7XTYkuTkREDl9lHqYLbiRYfNh3iSlHREQkorw3LBQRETmoFFAiIhJKCigREQklBZSIiISSAkpEREJJASUiIqGkgBIRkVBSQImISCgpoEREJJQUUCIiEkoKKBERCSUFlIiIhJICSkREQkkBJSIioaSAEhGRUFJAiYhIKCmgREQklBRQIiISSgooEREJJQWUiIiEkgJKRERCSQElIiKhpIASEZFQUkCJiEgoJTSgzGywmS02swwzuyNGu5nZk0H7HDPrHQxvbWYTzGyhmc03s5ujxulpZpPNbLaZTTezvolcBhERqRwJCygzSwKeAYYARwKXmdmRxboNAToHj+uBZ4PhBcBv3f0IoD/wq6hxHwbudfeewJ+C9yIicohJ5B5UXyDD3Ze6+y5gBDC0WJ+hwKseMRlINbN0d89095kA7r4dWAi0DMZxoH7wugGwNoHLICIilSQ5gdNuCayKer8a6FeOPi2BzKIBZtYO6AVMCQbdAowzs0eJBOyPYs3czK4nsldGmzZt9nMRRESksiRyD8piDPOK9DGzusC7wC3uvi0Y/AvgN+7eGvgN8GKsmbv78+7ex937pKWlVbh4ERGpXIkMqNVA66j3rdj3cFyJfcysOpFwet3dR0b1uRIoev8OkUOJIiJyiElkQE0DOptZezOrAQwHRhfrMxq4Ijibrz+w1d0zzcyI7BktdPe/FxtnLXBS8PpUYEniFkFERCpLwj6DcvcCM7sRGAckAS+5+3wzuyFofw4YC5wFZAB5wNXB6AOAnwJzzWx2MOwP7j4WuA54wsySgZ0EnzOJiMihxdyLfyx06OnTp49Pnz69sssQETlkmNkMd++TyHnoShIiIhJKCigREQklBZSIiISSAkpEREJJASUiIqGkgBIRkVBSQImISCgpoEREJJQUUCIiEkoKKBERCSUFlIiIhJICSkREQkkBJSIioaSAEhGRUFJAiYhIKCmgREQklBRQIiISSgooEREJJQWUiIiEkgJKRERCSQElIiKhpIASEZFQUkCJiEgoKaBERCSUFFAiIhJKCigREQklBZSIiISSAkpEREJJASUiIqGkgBIRkVBSQImISCgpoEREJJQUUCIiEkoKKBERCSUFlIiIhFJCA8rMBpvZYjPLMLM7YrSbmT0ZtM8xs97B8NZmNsHMFprZfDO7udh4NwXTnW9mDydyGUREpHIkJ2rCZpYEPAOcDqwGppnZaHdfENVtCNA5ePQDng2eC4DfuvtMM6sHzDCz8e6+wMxOAYYCR7t7vpk1TdQyiIhI5UnkHlRfIMPdl7r7LmAEkWCJNhR41SMmA6lmlu7ume4+E8DdtwMLgZbBOL8AHnT3/KB9QwKXQUREKkkiA6olsCrq/Wr+f8iUu4+ZtQN6AVOCQV2AE8xsipl9YWbHxZq5mV1vZtPNbHpWVtb+L4WIiFSKRAaUxRjmFeljZnWBd4Fb3H1bMDgZaAj0B24D3jazfabj7s+7ex9375OWlrY/9YuISCVKZECtBlpHvW8FrC1vHzOrTiScXnf3kcXGGRkcFpwKFAJN4ly7iIhUskQG1DSgs5m1N7MawHBgdLE+o4ErgrP5+gNb3T0z2CN6EVjo7n8vNs57wKkAZtYFqAFsTOByiIhIJUjYWXzuXmBmNwLjgCTgJXefb2Y3BO3PAWOBs4AMIA+4Ohh9APBTYK6ZzQ6G/cHdxwIvAS+Z2TxgF3Cluxc/dCgiIlWcHQ5/2/v06ePTp0+v7DJERA4ZZjbD3fskch66koSIiISSAkpEREJJASUiIqGkgBIRkVBSQImISCgpoEREJJQUUCIiEkoKKBERCSUFlIiIhJICSkREQkkBJSIioaSAEhGRUFJAiYhIKCmgREQklBRQIiISSgooEREJpcPihoVmth1YXNl1lFMTqsYt7KtKnaBaE0W1xl9VqROgq7vXS+QMEnbL95BZnOg7P8aLmU2vCrVWlTpBtSaKao2/qlInRGpN9Dx0iE9EREJJASUiIqF0uATU85VdQAVUlVqrSp2gWhNFtcZfVakTDkKth8VJEiIiUvUcLntQIiJSxSigREQklKpsQJlZIzMbb2ZLgueGJfS72czmmdl8M7slavg9ZrbGzGYHj7Oi2u40swwzW2xmZ4ag1kfMbJGZzTGzUWaWGgxvZ2Y7opbhuRDUWuL4lbhefxPUOc/M3jSzWsHwt6LW3XIzmx0Mr8z1WlKtB2V7jUOdYdxWS6o1VNuqmXWNWj+zzWxb0e9W2LbVMmqN37bq7lXyATwM3BG8vgN4KEaf7sA8oA6R73x9AnQO2u4BfhdjnCOBb4GaQHvgeyCpkms9A0gOXj9UND7QDpgXsvUac/xKXK8tgWVA7eD928BVMfo9BvypktdribUerO01DnWGbVstrdZQbavF+icB64C2YdxWS6s1nttqld2DAoYCrwSvXwGGxehzBDDZ3fPcvQD4Aji/HNMd4e757r4MyAD6Vmat7v5xMAxgMtDqAOtJWK2ljF9Z6xUiIVrbzJKJhOra6EYzM+AS4M0DrKc0cam1hOnGc70eUJ0h3Fah5HUaxm21yGnA9+6+InpgyLbVIjFrLWG6FVqvVTmgmrl7JkDw3DRGn3nAiWbW2MzqAGcBraPabwwORbwUtRvbElgV1Wd1MKyyay3yM+CjqPftzWyWmX1hZiccYJ3xqLWk8Stlvbr7GuBRYCWQCWx194+LdTsBWO/uS6KGHfT1Wo5aD8b2Gq91CiHYVsuoNVTbajHDiR1CodhWy1FrXLbVUF/qyMw+AZrHaLqrPOO7+0IzewgYD+QQ2b0s+u/uWeA+wIPnx4j8QlmsSVVyrUXzuCsY9nowKBNo4+7ZZnYs8J6ZHeXu2yq71lizjTWpsuZ1oLUGvxxDiRxS2AK8Y2Y/cff/RHW7jL1/wSplvZZRa9y214OxTsOyrZbz57/PaDGGJXxbjZpODeA84M4YzaHYVsuoNW7baqgDyt0HldRmZuvNLN3dM80sHdhQwjReBF4MxvkbkdTG3ddHTesFYEzwdjV777m0ouzDLAmtNXh/JXAOcJoHB3TdPR/ID17PMLPvgS5AqdfISnCtJY1fWet1ELDM3bOCcUYCPwL+E7xPBi4Ajo2aZ2Wt1xJrjef2ehDWaZi21dJqDdu2WmQIMDP6Zx5MI0zbaom1xnNbrcqH+EYDVwavrwTej9XJzJoGz22I/HDfDN6nR3U7n8hhq6LpDjezmmbWHugMTK3kWgcDtwPnuXteVP80M0sKXncIal1ambWWMn5lrdeVQH8zq2NmRuR4+cKo9kHAIneP/megstZribUexO31QOsM27Za2s8/bNtqkeJ7SUXCtK2WWGtct9XSzqAI8wNoDHwKLAmeGwXDWwBjo/p9BSwgchjqtKjhrwFzgTnBikuParuLyBkmi4EhIag1g8ix29nB47lg+IXA/KD/TODcENQac/xKXq/3AouCX5TXgJpRbS8DNxSbbmWu15i1HqztNQ51hnFbLanWMG6rdYBsoEGMaYRtW41Zazy3VV3qSEREQqkqH+ITEZFDmAJKRERCSQElIiKhpIASEZFQUkCJiEgoKaBE9pOZ5ZSjzy0WuRzUQWdmqWb2y8qYt0g8KKBEEusWIt8XKbeiL17GQSqggJIqSwElcoDM7GQz+9zM/muReyG9bhG/JvLlxglmNiHoe4aZfWNmM83sHTOrGwxfbmZ/MrOJwMVmNjjo862ZfRr0SbHIxTenWeTioEOD4VeZ2ftm9j+L3Gfnz0FpDwIdLXJPnkcO+ooROUChvhafSBXSCziKyLXFJgED3P1JM7sVOMXdN5pZE+BuYJC755rZ7cCtwF+Caex094FmlkbkqgAnuvsyM2sUtN8FfObuP7PIjQCnWuSinxC5bUF3IA+YZmYfErmXT3d375ngZRdJCAWUSHxM9eAaaRa522k7YGKxPv2J3LRtUuSycNQAvolqfyuq35ceuWcO7r4pGH4GcJ6Z/S54XwtoE7we7+7ZwfxHAgOB9+KwXCKVRgElEh/5Ua/3EPt3y4gEyWUlTCM3ql+sa5AZcKG7L95roFm/GP11DTOp8vQZlEhibQfqBa8nAwPMrBNAcIXtLjHG+QY4KbjiM1GH+MYBNwVX5cbMekWNc7qZNTKz2kTugDqp2LxFqhwFlEhiPQ98ZGYTPHJPoquAN81sDpHA6lZ8hKDf9cBIM/uW/3/o7z6gOjDHzOYF74tMJHIV6dnAu+4+PTjkN8nM5ukkCamKdDVzkSrOzK4C+rj7jZVdi0g8aQ9KRERCSXtQIiISStqDEhGRUFJAiYhIKCmgREQklBRQIiISSgooEREJpf8HAbSJeXhPhQgAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"joint_posterior = posterior_pmf.unstack()\n",
"\n",
"plot_contour(joint_posterior)\n",
"decorate(title='Joint posterior distribution')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can confirm that the parameters with maximum posterior probability are consistent with the point estimates."
]
},
{
"cell_type": "code",
"execution_count": 87,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(0.0318, -0.85)\n",
"[ 0.0316906 -0.84909407]\n"
]
}
],
"source": [
"print(posterior_pmf.max_prob())\n",
"print(results.params.values[::-1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here are the means of the marginal distributions."
]
},
{
"cell_type": "code",
"execution_count": 88,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-0.8491556076374563, 0.03169425569628929)"
]
},
"execution_count": 88,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"marginal_inter = marginal(joint_posterior, 0)\n",
"marginal_slope = marginal(joint_posterior, 1)\n",
"\n",
"marginal_inter.mean(), marginal_slope.mean()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Recall that the intercept indicates the log odds of the hypothesis at `x=0`.\n",
"To make the distribution of intercepts easier to interpret, I'll use `expit` to transform the values to probabilities."
]
},
{
"cell_type": "code",
"execution_count": 89,
"metadata": {},
"outputs": [],
"source": [
"marginal_probs = transform(marginal_inter, expit)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And here's what it looks like."
]
},
{
"cell_type": "code",
"execution_count": 90,
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1C0lEQVR4nO3deXwdd3nv8c+jzbLlRV7kTd6XLM6emCQkrFm4SYCE9tLbhB1uSXPLFkpLA5cWSmkLF9qylhD2hCXQEMAFl4RQEsiK7ax27DjyLm+SF8myZFnbc/+Y0dHoWMuRdObMOUff9+ull2afZ+Ysz/n95je/MXdHREQk35QkHYCIiMhAlKBERCQvKUGJiEheUoISEZG8pAQlIiJ5SQlKRETykhJUnjOzN5vZ/TFu/ztm9qlw+OVm9kIWt/1fZvb2cPgdZvZwFrcd63kZYr+Xm9mLZnbczN4wwPxNZvaqXMeVa+PlOKPM7HQze8rMWszs/THva4mZuZmVjXJ9N7MVg8zr99mJLmtmt5vZ344u6uwz3QfVx8x2AnOAbqAVWAu8z92Pj3J7nwBWuPtbshVjtpnZd4B6d//YCNb5BCM8LjN7B/Bn7v6yUcS4BNgBlLt710jXzyYz+w2wxt2/kIVt7SQ4Jw+MObAsivt9a2YOrHT3uji2Hxcz+yZwzN0/mIN9LWEM7/mRnOPBlg1/gHzP3ReMdP/ZohLUqV7v7pOBC4GXABl/cWfbaH89heuamSXy+ia57xxYDGxKOogiP8eDGstnIguy9tonfByFw931F/4BO4GrIuOfBX4RDl9P8OZsAh4Ezows9zfAXqAFeAG4ErgG6AA6gePAM+Gy04BvAvvDdT4FlIbz3gE8AvwbcCSc9w7g4ci+LgPWAc3h/8si8x4E/jHcxgmCX8Hpx3gB8GQY64+Au4FPhfNeRVCaGs1xnbLvcNqfpR3bl8LYtwBXDnHuP0Hw6w1gN+Dh/o4DLx3FefmHcP8twP3ArCHeB+8G6sLXYA0wP5y+DegJj+84MGGo91B4DD8G7gz3uwlYHc67K21bHw6nXwo8SvA+ewZ41VCvL3AW8Osw1oPAR8NlS4DbwpgPh3HMCOctCc/nzcA+gvfih8J5A76+oz3OAdb7Xbjv1nD7fxpOfx3wdHjcjwLnpu3rb4BngZPhcTvwTmAPcBS4heAH5bPhNr4cWX8F8FD43jgE/GiI137Azznw3wQ1K+1h3KcNsO6DwD8Dfwj39fMBzvn/Jng//y58jT4G7AIawvM3bbjXKJx/MfBYGOd+4MtARWS+A+8HtofH/FmgJPJZfDht2RXh8HcIvneqCN5jPfR97uYDbcDMyLoXAY0EJb3sfyfHsdFC/aP/h25h+Eb9B+C08AN1NVAOfJjgC6wCOD38kPR+iS0Blkc+uN9L28fPgK+Fb4DZ4Zv5zyNvnC7gfUAZMDH6ZgJmEHwY3xrOvykcnxn5gOwm+NIqS3/ThPHuAj4YHscbCb6ITklQoziuU/bNqQmqK7LvPyX4EM9IP/fp+6Dvw1oWmT/S87ItfB0nhuOfHuQ9cAXBB/pCYAJBQv3dQO+RDN5DnyD4QrsOKCX48np8sG0BtQTJ5DqCL6+rw/GaQc7xFMIvLqAyHL8kXPZW4HFgQXgcXwN+mHY+f0jwPjyH4EvmqvRzn43jHGDd1BdiOH4hwRf0JeH6bw+3PyGyr6cJPpMTI/HfHh73a8L9/4zgM1Ubbu+V4fo/BP5veE4rgZcNEtegn/PI+f+zIY7rQYIfdGeH5/UnnPoevjOcNxF4V7j9ZcBk4F7grgxfo4sIfsyUhctuBm5NO8e/JfhsLAK20v+zOGSCSv8+iCy7Fvg/kfF/A74U13fyuKsiyMDPzKwJeJjgV9c/EXyZ/tLdf+3uncDnCN5glxH8qpoArDKzcnff6e7bBtqwmc0BriV4I7W6ewPBC3xjZLF97v4ld+9y9xNpm3gt8KK73xXO/yFBSeT1kWW+4+6bwvmdaetfSvDB+7y7d7r7PQSljYFkfFwZ7huCL43eff+IoFT22mG2mYlMzsu33X1reE5/DJw/yLbeDHzL3Z9095PAR4CXhtcERuNhd1/r7t0Epabzhlj2LcDacPked/81sJ7gi79X6hwTlDoOuPu/uHu7u7e4+xPhcn8O/F93rw+P4xPAG9Oqlv4+fB8+B3ybILGP1kiOM927ga+5+xPu3u3u3yUoKV0aWeaL7r4n7TPxD+Fx30+QWH7o7g3uvhf4PUFtAQQ/whYT/Nhqd/fBGusM9TnP1F3uvtHdW4G/Bf6XmZVG5n8iPOcnCN5r/+ru2z24zv0R4MZMXiN33+Duj4fv950EP0BemRbLZ9z9iLvvBj7P2F7fXt8leJ8SHtdNBK93LJSgTvUGd69298Xu/hfhG2k+QckDAHfvIShd1HpwYfFWgi+ABjO728zmD7LtxQQJYr+ZNYWJ8GsEv/p67Rkitn5xhHYR/GLMdP29Hv70iax/ihEeVyb7ZpB9D7fNTGRyXg5EhtsIfrEOu63wi+Nw2rZGIn2/lUNcf1gM/EnveyN8f7wMmBdZJnqOFxKUDAfb1k8j29lM8KNjziDbGutrMZLjHCjWD6Ud98K0eAZ6bx2MDJ8YYLz3Nf4wYMAfwtaH7xokjkE/5xkeR3qcuwg+77MGmZ/+vt1FUCIa9jUys9PM7BdmdsDMjhH8kI7uZ9B1x+jnBD9alxGUNJvd/Q9Z2O6AlKAys4/gQwQEF6gJPkB7Adz9Bx60TltMUFz+TLiop21nD8Evw1lhEqx296nuflZkmfR1Bo0jtKg3jgzW3w/UhvFH1x/QCI4rk30zyL73hcOtwKTIvLkj2G4m5yVT6a91FTBzlNsazkDvj7si741qd69y908Pss4eYPkg294DXJu2rcqwdNFrYWQ4+loMd76zbQ/wj2mxTgpLwr1GHZO7H3D3d7v7fIKS5b8P0gR7yM95htLPaSdBlXEqnMH2Fy7fRf9EO9hr9FWCWoKV7j4V+ChBEh4qln2MzCnn3N3bCWog3kxQpR5b6QmUoDL1Y+C1ZnalmZUT1PmfBB4N7424wswmENSDnyD4pQrBG21Jb2srd99PcIH+X8xsqpmVmNlyM0svmg9mLXCamb3JzMrM7E+BVcAvMlz/MYIPwPvD9f+Y4GLrKUZyXCMwO9x3uZn9CXBmeEwQXGO4MZy3muD6WK9Ggou1ywbZ7ljPS9QPgHea2fnhsf8T8ERYjZJtB+l/TN8DXm9m/8PMSs2s0sxeZWaDNfP9BTDXzG41swlmNsXMLgnn3Q78o5ktBjCzGjO7IW39vzWzSWZ2FkGDgx9F4hrN65up9OP+OnCLmV0Stk6sMrPXmtmUbOzMzP4kcg6PEnzxdg+w6KCf8xHs7i1mtsrMJgGfBO4Jqz0H8kPgg2a21MwmE7zXfuT9m5UP9hpNAY4Bx83sDOD/DLD9vzaz6Wa2EPhAZN1MHQRmmtm0tOl3ElzHup7gPRsbJagMuPsLBPWuXyL4NfR6guboHQTXaT4dTj9A8CX80XDV/wj/HzazJ8PhtxE0Vnie4MNyD/2rcIaK4zDBdYcPEVQ7fRh4nbsfGnLFvvU7gD8meHMdJahzv3eQxUd6XJl4AlgZbvMfgTeGxwRBff3yMK6/J0gUvXG3hcs/ElYBRa9NjPm8pG3rN2EsPyEocS6n/zXCbPpn4GPhMf2Vu+8BbiA4z40EJYu/ZpDPqbu3EFSzvJ7gNXoReHU4+wsELRDvN7MWggYTl6Rt4iGCi/S/AT4XXsuB0b++mfoE8N3wuP+Xu68nuA71ZYLXv47gPZotLwGeMLPjBOfkA+6+I32hYT7nmbqLoKHBAYIGGUPd0PutcPnfEdzz1E7QQCpqsNfor4A3EbSa/DoDJ5+fAxsIfvz9kqD1cMbcfQtBEt0evlbzw+mPEPxgfDKmH24pulFXZJyxPLrxuZiY2YMErfa+kXQscTOz/wZ+EPex6mYxERHJmJm9hODWgPQq46xTFZ+IiGTEzL4LPEBwq0xL7PtTFZ+IiOQjlaBERCQvFdU1qFmzZvmSJUuSDkNEREZgw4YNh9y9Jn16USWoJUuWsH79+qTDEBGRETCzAXu0URWfiIjkJSUoERHJS0pQIiKSl5SgREQkLylBiYhIXlKCEhGRvKQEJZKHTrZ3ol5eZLwrqvugRIrB+ge2ce8XH6dmwVTe+rFXMmv+1KRDEkmESlAieaSzo4tfffsp3J2GPc3c/uH72bf9SNJhiSRCCUokjzzzu120HmtPjbc2t/P1jzzA9o0Hh1hLpDgpQYnkCXfn0TVbTpne3tbBdz7+W55/Yk8CUYkkRwlKJE/s2NTA/h1HASifUMa7//lqps6YBARVf9//p9+z4TfbkgxRJKeUoETyxGP/+UJq+IJXL2XZ2XO4+TNXM3PeFAB6enq45/OP8fufbU4qRJGcUoISyQNHG46z6bG+KrzLrz8DgJlzp/Dnn3kN85ZOT81b+80N/P6nz+c8RpFcU4ISyQOP/WJr6r6nFefPY/bCaal5U6ZP5N3/fDVLzpqdmvbAD56jq6s753GK5JISlEjCTrZ3su7+utR4b+kpamJVBe/65BVMm1UFQEd7J3tfVPNzKW5KUCIJe/q3O2hv7QBg5rwpnL56/oDLlVeUsfKCeanxumcO5CQ+kaQoQYkkyN15NNI44qWvOx0zG3T5ZefOSQ1vf1YJSoqbEpRIguqeOUDDnmYAJkws56Krlg25/PJIgtq1+RCdHV2xxieSJCUokQRFb8y96KrlVE6qGHL5qTMmpRpQdHd1s/P5xljjE0mSEpRIQg7tO8aWdXsBMDNe+rrTMlpv2blzU8Pbn1UXSFK8lKBEEvLYL7amhk9fPT/jXstXnNdXzbdN16GkiClBiSSgva2DDQ/0dVt02etPbVo+mKVnz0k1pKjfepgTYQtAkWKjBCWSgI2P7ObkiU4AZi+cxorz5w6zRp9JUyYwb1nQs4S7s3NTQywxiiRNCUokAbu3HEoNX3DFsiGblg9kxXl9CW2brkNJkYo1QZnZNWb2gpnVmdltA8w/w8weM7OTZvZXkekLzey3ZrbZzDaZ2QfijFMk1+pfPJwaXnjazBGvH20osU037EqRii1BmVkp8BXgWmAVcJOZrUpb7AjwfuBzadO7gA+5+5nApcB7BlhXpCB1nOzi4K7g3iczo3bFjBFvY8lZNZSUBh/fAzuPcry5fZg1RApPnCWoi4E6d9/u7h3A3cAN0QXcvcHd1wGdadP3u/uT4XALsBmojTFWkZzZv/0oPT09AMyqnTrsvU8DmVBZzqLTZ6XGtz+naj4pPnEmqFog+gjQekaRZMxsCXAB8MQg8282s/Vmtr6xUTctSv6LVu8tWDny0lOvaLdHquaTYhRnghroqq+PaANmk4GfALe6+7GBlnH3O9x9tbuvrqmpGUWYIrnVP0GN/PpTr+Xn6TqUFLc4E1Q9sDAyvgDYl+nKZlZOkJy+7+73Zjk2kcT0S1CnzRpiyaEtOn0W5RVlABze30JTY+uYYxPJJ3EmqHXASjNbamYVwI3AmkxWtKDN7TeBze7+rzHGKJJTJ1o7OLQ3qAwoKS1h3tLqUW+rrLyUxav6ag3Uq4QUm9gSlLt3Ae8F7iNo5PBjd99kZreY2S0AZjbXzOqBvwQ+Zmb1ZjYVuBx4K3CFmT0d/l0XV6wiubK3ru8hg3MXV6dKQKOlaj4pZmP7dAzD3dcCa9Om3R4ZPkBQ9ZfuYQa+hiVS0PpX743++lOv6OM3tj17EHcf8U2/IvlKPUmI5NDeLDWQ6DV/+Qwqq4Jm6scOt3FoX8uYtymSL5SgRHJoz9bsJqjS0hKWnj07Na5qPikmSlAiOdLSdILmQ0FLu/KKMmYvmpaV7S6PdnukhhJSRJSgRHJk74t9DSTmL59OaWl2Pn7RhhLbw+tQIsVACUokR/Zs7evBPBvVe73mLJpG1bRKANpaTrJ/x9GsbVskSUpQIjkSbcFXm8UEZWantOYTKQZKUCI54O79qviyWYKC/o/f2L1FfVJKcVCCEsmBpsZWWo8Fj8SorKpg1vwpWd1+7fK+TmcP7GzK6rZFkqIEJZID0ebltStmZP1m2tmLpqW2eXhfC50dXVndvkgSlKBEciDbN+imq5hQxsx5QanM3Tm4uznr+xDJNSUokRzI1iM2hjJ3SXVq+OCuplj2IZJLSlAiMXP3fp3EZqMPvoFEE9T+HU2x7EMkl5SgRGLWWH+Mkyc6AZhcPZFpMyfFsh+VoKTYKEGJxCz9Ee9x9TY+d8n01LBu1pVioAQlErNsPUF3ODPmTqZ8QvAEndbmdo43t8e2L5FcUIISiVkuGkhA0KPE3MXVqXHdDyWFTglKJEZdXd3s395X3bZg5Ywhlh67Of0SlKr5pLApQYnE6OCuZro6uwGYPnsyVVMrY93fvKXVqWGVoKTQKUGJxCjbj3gfTrQEpZZ8UuiUoERiVJ/lJ+gOp39T82Z6enpi36dIXJSgRGKUqwYSvaqmVjJ1RnCfVWdHF4f3H499nyJxiTVBmdk1ZvaCmdWZ2W0DzD/DzB4zs5Nm9lcjWVck33Wc7KIh7BPPzKhdEW8DiV5zFvc9Sl4NJaSQxZagzKwU+ApwLbAKuMnMVqUtdgR4P/C5Uawrktca9vRVsc2cP4UJE8tzst/oDbtqKCGFLM4S1MVAnbtvd/cO4G7ghugC7t7g7uuAzpGuK5LvGiI9is9ZNG2IJbNLXR5JsYgzQdUCeyLj9eG0rK5rZjeb2XozW9/YqCeJSv5o2NOXoGYvqs7ZfqMJSiUoKWRxJqiBOhzzbK/r7ne4+2p3X11TU5NxcCJxi5ZeclmCqlkwlZKS4KN95MBxTranV1CIFIY4E1Q9sDAyvgDYl4N1RfJC9KGBs3OYoMoryphV2/fwwgY9vFAKVJwJah2w0syWmlkFcCOwJgfriiTuZHsnRw8GTbxLSkqoqZ2a0/2rmk+KQWwJyt27gPcC9wGbgR+7+yYzu8XMbgEws7lmVg/8JfAxM6s3s6mDrRtXrCLZ1rjnWGp4xrzJlJWX5nT/asknxaAszo27+1pgbdq02yPDBwiq7zJaV6RQHNzdlBqO9jCeK7oXSoqBepIQiUFDpASVy+tPveYt7V+Ccs+0fZJI/lCCEolBQ6QElcsWfL2qa6pSNwa3tZyk5eiJnMcgMlZKUCIxOLgrmXugepmZGkpIwVOCEsmyk+2dHG3oa8E3a/6UROLQozek0ClBiWRZtAXfzPlTct6Cr1e0BLV/hxpKSOFRghLJsoMJX3/qFW1qHq1yFCkUSlAiWdaQUA8S6aJNzQ/ubqa7Ww8vlMKiBCWSZdEujuYkcA9Ur0mTJzBtVhUA3V3dHN7XklgsIqOhBCWSZf1KUAuTK0FBepdHug4lhUUJSiSL+rXgKy1JddqaFDU1l0KmBCWSRdEWfLPmT6GsLJkWfL2UoKSQKUGJZFG0BV+SDSR6Ra+BKUFJoVGCEsmiaHPuOQn0IJGuZsFUSsNS3NGG47S3dSQckUjmlKBEsijfSlBlZaXULOh7FpXuh5JCogQlkkXRFnxJ3qQb1e/RG+rySAqIEpRIlpw80UlTYysQtOCbmVAffOmiVY16/LsUEiUokSxp2NP35Z8PLfh69e9Roim5QERGSAlKJEvypQeJdCpBSaFSghLJkny8/gQwfU4V5RVlALQcPUFby8mEIxLJjBKUSJZEq89qEu7iKKqkpKR/Sz6VoqRAKEGJZEm0CffcPKrig/7XoRp0HUoKRKwJysyuMbMXzKzOzG4bYL6Z2RfD+c+a2YWReR80s01mttHMfmhmlXHGKjIW7W0dNB8KWvCVlpUyY97khCPqL/rYeZWgpFDElqDMrBT4CnAtsAq4ycxWpS12LbAy/LsZ+Gq4bi3wfmC1u58NlAI3xhWryFg15FkffOmi18TUUEIKRZwlqIuBOnff7u4dwN3ADWnL3ADc6YHHgWozmxfOKwMmmlkZMAnYF2OsImPSkGc9SKSLxqSbdaVQxJmgaoE9kfH6cNqwy7j7XuBzwG5gP9Ds7vcPtBMzu9nM1pvZ+sbGxqwFLzISB/O0BV+vGXMmUz4haMnX2txO67H2hCMSGV6cCcoGmOaZLGNm0wlKV0uB+UCVmb1loJ24+x3uvtrdV9fU1IwpYJHRypfHvA/GzPo9PFHXoaQQxJmg6oGFkfEFnFpNN9gyVwE73L3R3TuBe4HLYoxVZEz6l6CqkwtkCLoOJYUmzgS1DlhpZkvNrIKgkcOatGXWAG8LW/NdSlCVt5+gau9SM5tkZgZcCWyOMVaRUcv3Fny9oiW7g7oOJQWgLK4Nu3uXmb0XuI+gFd633H2Tmd0Szr8dWAtcB9QBbcA7w3lPmNk9wJNAF/AUcEdcsYqMRb8WfLX514Kv1xw1NZcCE1uCAnD3tQRJKDrt9siwA+8ZZN2PAx+PMz6RbOjXgi+PepBI1/9mXSUoyX/qSUJkjA7sys9OYtNV11RRUVkOQOuxdo43qyWf5DclKJExipag8rGJea+gJV/06bpNyQUjkgElKJExOlggJShIayihaj7Jc0pQImPQ1nKSY0faACgrL2Vmnrbg6xXtxFbXoSTfKUGJjEG026DZC6dRUpLfH6n+ncY2JRaHSCby+9Mkkuei13HmLqlOLI5Mpd+sGzSkFclPSlAiY3BgZ1NqON+vPwFMmzWJCRODlnxtLSc53qSWfJK/lKBExiBagiqEBGVmaighBWPIBGVm90eGPxJ/OCKFw937fcEXQhUf9O9RQk/XlXw2XAkq2j34n8QZiEihaT7URntrBwCVVRVMnTEx4YgyoxKUFIrhEpSuoIoM4kBaA4mgX+P8N0edxkqBGK4vvmVmtobguU29wynufn1skYnkuWgDibkFcP2pV3oJyt0LJrnK+DJcgoo+ov1zcQYiUmgKrYFEr2kzJ1E5qYL2tg7aWztoOXqCqTMmJR2WyCmGTFDu/lDvsJnVhNP0XHURCu8eqF69Lfl2bwk+ygd3NytBSV4arhWfmdnHzewQsAXYamaNZvZ3uQlPJD91dXX3ew5UIZWgQNehpDAM10jiVuBlwEvcfaa7TwcuAS43sw/GHZxIvjqy/zjdXd0ATJtVxcSqioQjGpnZevy7FIDhEtTbgJvcfUfvBHffDrwlnCcyLh3od/0pfx+xMZhoiU9NzSVfDZegyt39UPrE8DpUeTwhieS/fi34lkxPLpBRUp98UgiGS1Ado5wnUtQO7DyaGp5bgCWoKdMnUhlWS7a3ddB8uC3hiERONVyCOs/MjplZS/h3rHccOCcXAYrko+h1m0JrIAFBS770UpRIvhkyQbl7qbtPdfcp4d/UyLiq+GRcOtneyZEDxwEoKSmhZsHUYdbIT/2uQ6kln+Sh4ZqZV5rZrWb2ZTO72cyGu7E3ff1rzOwFM6szs9sGmG9m9sVw/rNmdmFkXrWZ3WNmW8xss5m9dCT7FolL9JrNzPlTKK8Y0ccib8xRn3yS54ar4vsusBp4DrgO+JdMN2xmpcBXgGuBVcBNZrYqbbFrgZXh383AVyPzvgD8yt3PAM4DNme6b5E4FeoNuunm6PHvkueG++m3yt3PATCzbwJ/GMG2LwbqwmbpmNndBF0nPR9Z5gbgTg9+jj4elprmAa3AK4B3ALh7B2qUIXniwK7Cvv7US33ySb4brgTV2Tvg7l0j3HYtsCcyXh9Oy2SZZUAj8G0ze8rMvmFmVQPtJKx6XG9m6xsb1QuTxK/QW/D1mjytkqqplQB0RK6rieSLTFvx9bbcOzfSqu/YMOsO9FMs/WaLwZYpAy4EvuruFxCUqE65hgXg7ne4+2p3X11TUzPQIiJZdXBX9CGFhXcPVC8zY/7yvvj3bjuSYDQip8q0FV9vy72yyPBwTZfqgYWR8QXAvgyXqQfq3f2JcPo9BAlLJFHHm9s53nQCgPKKMmbMnZxwRGMzf/mM1PA+JSjJM8OVoMZiHbDSzJaaWQVwI7AmbZk1wNvC1nyXAs3uvt/dDwB7zOz0cLkr6X/tSiQRB9O6OCr0aza1K/oS1N46JSjJL7G1j3X3LjN7L3AfUAp8y903mdkt4fzbgbUErQPrgDbgnZFNvA/4fpjctqfNE0lEoT4DajDpJSg1lJB8EusNHO6+liAJRafdHhl24D2DrPs0QRN3kbwR7YOvGBLUjDmTqayqoL21g7aWkzQfaqO6ZsD2SCI5F2cVn0jROVAk90D1MjNqI6UoNZSQfKIEJZIhd+/fgq8ISlDQv5pP16EknyhBiWToaEMrHe3BrYGTpkxgcnVlwhFlR7ShhFryST5RghLJUL8bdJdUF01jAjU1l3ylBCWSoWK5QTfdrPlTqKgMHk7QcvQEx47o2VCSH5SgRDIUbWIe7ceu0JkZ85f1Jdx9244OsbRI7ihBiWSo2FrwRfW7YVfVfJInlKBEMtDV1U1jfV/3k3OKqAQF6S35DicYiUgfJSiRDDTWH6OnuweA6bMnUzmpIuGIsqt/QwlV8Ul+UIISyUD/HiSKq/QEULNgaurJwM2HWjne3J5wRCJKUCIZqd/aV+0VLW0Ui9LSEuYurU6Nq7m55AMlKJEM7HnhUGp40RnF+dyxWt0PJXlGCUpkGJ0dXezb3nddZuFpMxOMJj5qySf5RglKZBj7th+lu6sbgFm1U5k0ZULCEcVDffJJvlGCEhlGv+q902clGEm8Zi+aRmlZKQBHDx6n7fjJhCOS8U4JSmQYu7f0JaiFRZygyspK+92AvH+7mptLspSgRIbRv4FE8SYoUMexkl+UoESGcOxIG02NrQCUV5QVxVN0h1K7vK9PPl2HkqQpQYkMIVq9t+C0mZSWFvdHZr6erit5pLg/bSJjFK3eK+brT73mLqmmpCT4Wji09xjtbR0JRyTjmRKUyBB292vBV5z3P0WVV5T1e5TI/h1NyQUj416sCcrMrjGzF8yszsxuG2C+mdkXw/nPmtmFafNLzewpM/tFnHGKDKSrq5u9L/ZVc42HEhTA/H7XodSzuSQntgRlZqXAV4BrgVXATWa2Km2xa4GV4d/NwFfT5n8A2BxXjCJDObirmc6OLgCqa6qYOmNSwhHlhro8knwRZwnqYqDO3be7ewdwN3BD2jI3AHd64HGg2szmAZjZAuC1wDdijFFkUOPt+lOv2hV9VZl69IYkKc4EVQvsiYzXh9MyXebzwIeBnqF2YmY3m9l6M1vf2Ng4poBFonZv6Xs/Ffv9T1Fzl1ZjZgA07Gmm42RXwhHJeBVngrIBpnkmy5jZ64AGd98w3E7c/Q53X+3uq2tqirOXaUnGnhf6rr+MpxLUhMpyahZMBcDd2b9DpShJRpwJqh5YGBlfAOzLcJnLgevNbCdB1eAVZva9+EIV6a+t5SSH9gWPeC8tK2X+sunDrFFc1KOE5IM4E9Q6YKWZLTWzCuBGYE3aMmuAt4Wt+S4Fmt19v7t/xN0XuPuScL3/dve3xBirSD/R5uXzl01PPW12vIg+ekPXoSQpsX3q3L3LzN4L3AeUAt9y901mdks4/3ZgLXAdUAe0Ae+MKx6RkRivDSR6RUtQO59vSDASGc9i/Vno7msJklB02u2RYQfeM8w2HgQejCE8kUFFuzgaTw0kei08bSblFWV0dnRxaO8xmhpbqa6pSjosGWfUk4RIGnenfuv4bCDRq7yijCVn9TU6qnt6f4LRyHilBCWSprG+rw+6qmmVTJ89PksOKy+Ylxre+qQSlOSeEpRImvT7n3rvCRpvVkQS1LZnDtDTM+QtiSJZpwQlkmZ35P6nRWeM33vr5i6uZnL1RCBodq/WfJJrSlAiafo9QXccXn/qZWasvGBualzXoSTXlKBEIk6e6OTgriYg+IKuXTlj6BWK3Irz+6r56p4+kGAkMh4pQYlE7Nl6mODuh+DhfRMqyxOOKFkrzu8rQe18vpGT7Z0JRiPjjRKUSMR4v0E33dQZk5izuBqA7q5udm5Sh8ySO0pQIhHjtQfzoUSr+V58StehJHeUoERC7j5uezAfSvR+qDolKMkhJSiR0J6th2k91g5A1dRKamqnJhxRflh69mxKy0oBOLi7iWNH2hKOSMYLJSiR0MZHdqeGz7xkwbi9QTddxYQylqzqux9M1XySK0pQIgTVe5se63u481mXLRxi6fEn2quEmptLrihBiQAHdjZx5EALABMmlrPivLnDrDG+RJub1z19INUUXyROSlAiwMZH+6r3znhJLWXlpQlGk39ql89g0pQJABxvOsGBnU3JBiTjghKUCPSr3jv78kUJRpKfzCytFKXrUBI/JSgZ9xr3Hkt1b1ReUcbKC+cNvcI4pcdvSK4pQcm4Fy09rbxw3rjv3mgw0Rt2dz3fSGdHV4LRyHigBCXj3qbI9aez1XpvUNU1VcwK7w3r7Ohi5/Pq9kjipQQl41pTYyv1Lwa9R5SUlnD6S2oTjii/qVcJySUlKBnXnn+8r3pvxXlzmTR5QoLR5L+Vuh9KcijWBGVm15jZC2ZWZ2a3DTDfzOyL4fxnzezCcPpCM/utmW02s01m9oE445Txa+OjfQlq1UtVvTecZefMoaQ0+NrYt/0ILU0nEo5IillsCcrMSoGvANcCq4CbzGxV2mLXAivDv5uBr4bTu4APufuZwKXAewZYV2RMWppOsHNTAxA0o1516YKEI8p/EyaW93vKsKr5JE5xlqAuBurcfbu7dwB3AzekLXMDcKcHHgeqzWyeu+939ycB3L0F2Azo4oBk1eYn6lM9IixeVcOU6okJR1QYTrtofmp43f3bEoxEil2cCaoW2BMZr+fUJDPsMma2BLgAeGKgnZjZzWa23szWNzaqVZFkrl/fe6rey9iFVy5LVfPt2HiQfduPJByRFKs4E9RAXUGnd+A15DJmNhn4CXCrux8baCfufoe7r3b31TU1NQMtInKKE60dbHvmYGr87MvUe0Smps2cxDmR3jYeWbMlwWikmMWZoOqB6M/SBcC+TJcxs3KC5PR9d783xjhlHNqybi/dXd0ALFg5k+qaqoQjKiyXXX9GaviZh3ZxvLk9wWikWMWZoNYBK81sqZlVADcCa9KWWQO8LWzNdynQ7O77LXgQzzeBze7+rzHGKONU9OZcVe+N3KLTZ6WeONzd1c0T//ViwhFJMYotQbl7F/Be4D6CRg4/dvdNZnaLmd0SLrYW2A7UAV8H/iKcfjnwVuAKM3s6/LsurlhlfDnZ3snWDX2tz85S9d6oXB4pRT2xditdYYlUJFvK4ty4u68lSELRabdHhh14zwDrPczA16dExuzFJ/en+pGbvXCaHu0+SmddtpCpMyZx7EgbLUdP8Nzvd3HBq5clHZYUEfUkIeOKu/P42q2pcT1aY/TKykq59LWnpcYf/c8X9CBDySolKBlXXnxqP9ueCbroKSkp4YJXL004osJ28TUrKK8IKmLqXzzM7i2HEo5IiokSlIwb7s6vvvNUanz1a5Yza76q98aiamol571ySWpcTc4lm5SgZNx4+sEd7N9xFIDyCWVcedM5CUdUHC67/vTU8KZH99DU2JpgNFJMlKBkXOjs6OL+u55Jjb/sDWcydcakBCMqHvOWTGf5ucHj4Ht6enjsl1uHWUMkM0pQMi48vvbF1C/7qqmVvOKPz0w4ouISLUWtu6+Ok+2dCUYjxUIJSoreidYOHvzxxtT4FTeeTeWkigQjKj5nvKSWGXOnAHDi+EmefnBnsgFJUVCCkqL30D2baGs5CcCMuVO4+NqVCUdUfEpKSnjp6yJNztdsUZNzGTMlKClqTYdaeXTNC6nx17z1PMrKShOMqHitvno5EyaWA9Cwp5mH7tmUcERS6JSgpKj95gfPpXqNqF0xk3NfvjjhiIpX5aQKLnt937Wo++96hhc2pPcPLZI5JSgpWgd3N7Hhgb4H6l37zgsI+iGWuFxx0zksOWs2ENx3dvdnH+bw/paEo5JCpQQlRetX33kqdR3ktIvmp5pCS3zKykp5020vZ9qs4PEl7a0d3PWph9SqT0ZFCUqKjrtz351Ps2XdXgDMjGvefkHCUY0fU6on8paPvoKy8uBa38HdTdzz+cfUaEJGTAlKioq784uvb+DB/+hrVn7RVcuZt3R6glGNPwtWzuQN77k4Nb7xkd089JPnE4xICpESlBSNnp4efvrlJ3j0P/v6gzvz4gVcf8vqBKMavy66cjmXvjbSaOLOp9VoQkZECUqKQnd3Dz/+10dZd39dato5L1vMmz7y8lRv25J7r/2zC9VoQkZNCUoKXldXN3f/v4d55qGdqWkXXrmMG//6ct3zlLDeRhNTZwb9Hra3dvDVv76PDb/ZpmtSMiwlKCloHSe7uOtTD7Hx0d2paZdedxpv/MBLKSnR2zsfpDeaaG1u557PP8bX/uZ+9m0/knB0ks+smH7FrF692tevX590GJIDx4608fjaF/nDf71I67H21PSX/9Eq3e+Up17YsI97v/Q4xw63paaZGZdcdxpXv+VcJk2ekGB0kiQz2+Dup1wsVoKSgrJ32xEe+flmnv39brq7uvvNu+LGc7jqTecqOeWxk+2d/PZHG3n4Z1v6vX5VUyv5H28/n3NfsZgJleUJRihJUIKSgtTV2U1j/TH2bjvChge2sXNTwynLTJ89mStuPJvVV69IIEIZjca9x1hz+zrqnt7fb3pJaQkLT5vJsnPmsPScOSw6Y5YS1jiQSIIys2uALwClwDfc/dNp8y2cfx3QBrzD3Z/MZN2BKEEVnq6ubtpbOzlxvIMTxztobW7n4K4mDuwM/hr3HqOnu2fAdZesms3lN5zBmZcsoLRU15sKjbuz6bE9/PIbGwZ9Cm9vwlqwciZTZkxkyvSJTK6uDP5Pr6Rq6gRdaywCOU9QZlYKbAWuBuqBdcBN7v58ZJnrgPcRJKhLgC+4+yWZrDuQ0Saohj3N/Pp7zwy/4DgWfZ9E3zKp6R4Me/gfh54ep7urh57unvC/09XZTU+309HexYnjHamOXDNVUlrCuS9fzOXXn8GClTOzcWiSsI6TXfz+p5vZ+MhuDuw8OqJ1zYwJE8spLS+hrLyU8opSyipKKQvHrcQoKTGsxDCjb9wMrP92gv/9xyUzq1+zgtMvmj/q9QdLUHHeIHIxUOfu28MA7gZuAKJJ5gbgTg++5R43s2ozmwcsyWDdrDne3N6vFZjklxlzpzBn8TQWnjaLi65apke1F5mKCWVceeM5XHnjORxvbmfnpga2P3eQ7c8d5OCupiHXdXfa2zpyE6gMatm5c2LZbpwJqhbYExmvJyglDbdMbYbrAmBmNwM3AyxatGhsEUvOlZSUMHFyBZVV5eH/CmoWTGXO4mrmLZ3OnEXTUs8YkuI3eVolZ1+2iLMvCz7LvQnr8P4WWo62c7zpBMePttNy9ATHm9pTD6KU4hRnghqojJxenzjYMpmsG0x0vwO4A4IqvpEE2Gv2wmm86baXj2bVopNp1Ua/xSysPrG+qhOzoCqlpNQoLSuhpLSEsvLgf2lZCeUTSpk4uYKKCWWqTpFB9SaswXR1ddN5spuuzm66Orrp6uxJDXd2dOM9fdXOPd0eVEP3eKpq2p1UnXUmVzuKqVFZNs1fNiOW7caZoOqBhZHxBUB6R1yDLVORwbpZM3laJedcrgfZiRSasrJS9RZSxOJs/rIOWGlmS82sArgRWJO2zBrgbRa4FGh29/0ZrisiIkUsthKUu3eZ2XuB+wiain/L3TeZ2S3h/NuBtQQt+OoImpm/c6h144pVRETyj27UFRGRRA3WzFx3uImISF5SghIRkbykBCUiInlJCUpERPJSUTWSMLNGYFfMu5kFHIp5H8VM52/0dO7GRudv9OI+d4vdvSZ9YlElqFwws/UDtTaRzOj8jZ7O3djo/I1eUudOVXwiIpKXlKBERCQvKUGN3B1JB1DgdP5GT+dubHT+Ri+Rc6drUCIikpdUghIRkbykBCUiInlJCSrCzK4xsxfMrM7Mbhtg/pvN7Nnw71EzOy8y74NmtsnMNprZD82sMrfRJ2uM5+4D4XnbZGa35jTwPJHB+bshPHdPm9l6M3tZpusWuzGeu2+ZWYOZbcxt1PljtOfPzBaa2W/NbHP42f1A1oNzd/0F1+FKgW3AMoIHJj4DrEpb5jJgejh8LfBEOFwL7AAmhuM/Bt6R9DEVyLk7G9gITCJ4/MsDwMqkjykPz99k+q4ZnwtsyXTdYv4by7kLx18BXAhsTPpYCu38AfOAC8PhKcDWbL/3VILqczFQ5+7b3b0DuBu4IbqAuz/q7kfD0ccJnvTbqwyYaGZlBF+2sT0BOA+N5dydCTzu7m3u3gU8BPxRjuLOF5mcv+MefhMAVYBnum6RG8u5w91/BxzJVbB5aNTnz933u/uT4XALsJngx3rWKEH1qQX2RMbrGfpk/2/gvwDcfS/wOWA3sJ/gycD3xxRnPhr1uSMoPb3CzGaa2SSCB1gujCXK/JXR+TOzPzKzLcAvgXeNZN0iNpZzJ1k6f2a2BLgAeCKbwSlB9bEBpg3YBt/MXk3wJfs34fh0gl8dS4H5QJWZvSWmOPPRqM+du28GPgP8GvgVQRVDVzxh5q2Mzp+7/9TdzwDeAPzDSNYtYmM5d5KF82dmk4GfALe6+7FsBqcE1aee/r/cFzBANZ2ZnQt8A7jB3Q+Hk68Cdrh7o7t3AvcSXHMZL8Zy7nD3b7r7he7+CoLqlhdjjjffZHT+eoXVUsvNbNZI1y1CYzl3MsbzZ2blBMnp++5+b9ajS/oiXb78EVxD2k5QCuq9WHhW2jKLgDrgsrTplwCbCK49GfBd4H1JH1MhnLtw3uzIMlsIG1OMl78Mz98K+i5UXwjsDd9rw65bzH9jOXeR+UsYv40kxvLeM+BO4PNxxVeWaSIrdu7eZWbvBe4jaNnyLXffZGa3hPNvB/4OmAn8u5kBdLn7and/wszuAZ4kqJ56inHUrcpYzl24iZ+Y2UygE3iP9zWmGBcyPH//E3ibmXUCJ4A/9eAbY8B1EzmQBIzx3GFmPwReBcwys3rg4+7+zQQOJRFjOX9hc/O3As+Z2dPhJj/q7muzFZ+6OhIRkbyka1AiIpKXlKBERCQvKUGJiEheUoISEZG8pAQlIiJ5SQlKxiUz6w57Z95oZv8RdrOU6brvMLMvj3B/xweZ/kkzuyocftDMVofDa82sOvz7i5Hsa4gYsrattO0uNbMnzOxFM/uRmVVkex8yPilByXh1wt3Pd/ezgQ7gluhMMyvNRRDu/nfu/sAA069z9yagGshWUsnmtqI+A/ybu68EjhJ0ZSUyZkpQIvB7YIWZvSp8vs0PCG4+rDSzb5vZc2b2VNiPYK+FZvar8Dk6H++daGY/M7MN4fNxbo7uxMz+xcyeNLPfmFlNOO07ZvbG9IDMbGfYncynCbqWedrMPmtmd5nZDZHlvm9m16etOzncx5Nh7L3L99tW2jovCZ/5U2lmVWH8Zw934iy46/oK4J5w0ncJ+msTGTP1JCHjWvh4lGsJOqqF4PEDZ7v7DjP7EIC7n2NmZwD3m9lp0eWANmCdmf3S3dcD73L3I2Y2MZz+Ew/6HawCnnT3D5nZ3wEfB96bQYi3hfGcH8b7SuCDwM/NbBpBn49vT1unHfgjdz8WJrnHzWxN+rai3H1duMyngInA99x9o5lNIUjgA3kT0AA0efCoFBh/valLjJSgZLyaGOme5ffANwm+7P/g7jvC6S8DvgTg7lvMbBfQm6B+HSYezOzecNn1wPvNrPd5VguBlcBhoAf4UTj9ewQdCo+Yuz9kZl8xs9nAHwM/iSSHXgb8k5m9ItxvLTAng81/ElhHkODeH+6vBTh/sBV6S4LpYWawL5FhKUHJeHUivSQR9hHYGp00xPrpX8JuZq8i6Nn+pe7eZmYPApUZrj8SdwFvBm5k4GcbvRmoAS5y904z2zlEHFEzCJ6eWh4u35pBCWozUG1mZWGiHG+9qUuMdA1KZHC/I/iyJ6zaWwS8EM672sxmhFV5bwAeAaYBR8PkdAZwaWRbJUDvtaY3AQ9nGEMLweO0o74D3AowSMew04CGMDm9Glg8xLai7gD+Fvg+QcMH3L0lbEwy0N/zYaerv40c29uBn2d4bCJDUoISGdy/A6Vm9hxB9dw73P1kOO9hgpLM0wTVbOsJrmOVmdmzBA91ezyyrVbgLDPbQNCo4JOZBBBWIz4SNof/bDjtIEHJ5duDrPZ9YLWZrSdIsFsG21YvM3sbQQ/zPyBoTPESM7sikxgJHj75l2ZWR9Bj/bjpDVzipd7MRQpMeM/Wc8CF7t6cdDwicVEJSqSAhDf1bgG+pOQkxU4lKBERyUsqQYmISF5SghIRkbykBCUiInlJCUpERPKSEpSIiOSl/w/MIpz/lv+qPwAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"marginal_probs.plot(color='C4')\n",
"decorate(xlabel='Probability at x=0',\n",
" ylabel='PDF',\n",
" title='Posterior distribution of intercept in terms of probability')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The mean of this distribution is about 24%, which is the predicted probability of supporting legalization for someone born around 1949."
]
},
{
"cell_type": "code",
"execution_count": 91,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.29961602953948796, 1948.7242)"
]
},
"execution_count": 91,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"marginal_probs.mean(), offset"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The estimated slope is the log of the likelihood ratio for each additional year of birth. To interpret slopes as likelihood ratios, we can use `np.exp` to transform the values in the posterior distribution."
]
},
{
"cell_type": "code",
"execution_count": 92,
"metadata": {},
"outputs": [],
"source": [
"marginal_lr = transform(marginal_inter, np.exp)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And here's what it looks like."
]
},
{
"cell_type": "code",
"execution_count": 93,
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA23ElEQVR4nO3deZxcZZn3/8+3t3R3esnW2XcIhIAIGNlVBtQBXHBcHsEVHAeZAZVRx1FnfiO/GR11xnHceEREXB5RXECfjKKgsoggmBDWEAKdhaSTTtLZekkn6e16/jinKqcq1d3V3XW6lr7er1e9+tRZr3PqdF117nOf+5aZ4ZxzzhWasnwH4JxzzmXiCco551xB8gTlnHOuIHmCcs45V5A8QTnnnCtInqCcc84VJE9QBUzSOyXdE+P6vyvpM+HwKyRtyOG6fy3pveHwlZL+mMN1x3pchtjueZJekNQl6U0jWO4GST+IMbTBttslael4bzefhvuMJG2R9Opw+FOSbgmHF0sySRUj3F7KcnGe98PEYZKOj3H9yWM1nkb0YZQySVuAWUA/cBC4C/igmXWNcn03AMeb2btGG5OZ3QbcNtrlR7itB4ETh5sv2/0ys0tyEZekxcBmoNLM+sJ1j9txSfOvwNfN7Ct52PaImVndaJbLdMyLSNafkZn9e643nqvzPp8kXQD8wMzmJ8bFcayy4VdQqd4Q/lOfAbwc+Od8BTLSX3Jpy0pSXj7bfG57HCwC1uU7iEI3lnM3B/wzGkKx/X8WTaDjycy2A78GTgGQ9EZJ6yQdkHS/pJMS80r6R0nbJXVK2iDpIkkXA58C3h4WNTwZztso6duSWsNlPiOpPJx2paSHJP23pH3ADelFBJLOlbRaUnv499zItPslfVbSQ0A3cEzRjqTTJa0NY/0xUB2ZdoGkllHu1zHbDse9P3Xz+loY+3OSLopMSBa7hO+jRWJ/CP8eCLd5ziiOy7+Fx7ZT0j2SZgzy0SPpbyQ1S9onaZWkueH4jeEx/Z8wjkkZlj3mmA2yjaHOpy2SPinpWUn7JX1HUvRzer2kJ8JlH5Z06hD7kiz2UVCce6OkX4XxPSrpuEEWPeaYh+t4n6T1YVx3S1qUtq1rJb0AvJA4nyR9XNLu8Jx/k6RLJT0fHt9PRZY/U9IaSR2Sdkn60hD7NerPKG09gxa9SnpL+FmcIqlM0ickbZS0V9JPJE0bZLn08x5JXwyP2WZJl0TGzw3j3xfuz99Epk2S9GVJO8LXl6P7I+kfwmO6Q9L7htnPTP+fV4WfZaekTZI+EM47meC7b254DLvCOFOO1VDncE6Zmb+C5p62AK8OhxcQ/Ar7N+AEgiK/1wCVwMeBZqCKoEhsGzA3XG4xcFw4fAPBZXJ0G78AvglMBmYCfwY+EE67EugDPkhQ9FoTjvtjOH0asB94dzj9ivD99HD6/cBW4ORwemXatquAF4G/D/fjrUAv8Jlw+gVASzg80v06ZtvhuPen7Vti228H2oFp6cc+fRvhtg2oiEwf6XHZGH6ONeH7zw9yDlwI7CG4gp4EfA34Q6ZzJMOyWR0zhjifItt4huAcnAY8FPmMzgB2A2cB5cB7w/knDRKTERTHAnwX2AecGR6n24DbB1ku0zF/UxjnSeHy/ww8nLat34Yx1xCcT33Av4T7+TdAG/BDoJ7gXDkMLA2X/xPw7nC4Djg7159Rhv/z6OeS3GfgqnBfE8fueuARYH64zW8CP8p0rDj2vO8N970c+FtgB6Bw+gPA/yb4oXhaeHwuCqf9a7jNmUAT8DDwb+G0i4FdBD+gJ4fHNPlZZ9jn+zn2//N1wHGAgFcRJK4z0r8LBvmfHPIczuXLr6BS/ULSAeCPBCfPvxN8mf7KzH5rZr3AFwn+Ac8luF81CVghqdLMtpjZxkwrljQLuAS43swOmtlu4L+ByyOz7TCzr5lZn5kdSlvF64AXzOz/hNN/BDwHvCEyz3fNbF04vTdt+bMJTqYvm1mvmf0MWD3Icch6v7LcNgRfrIlt/xjYEO7TWGVzXL5jZs+Hx/QnBF8GmbwTuNXM1prZEeCTwDkK7skMJ9tjNtT5lPB1M9tmZvuAzxIkXQi+6L5pZo+aWb+ZfQ84QvDZZuNOM/uzBfeVbmPw45DJB4DPmdn6cPl/B06LXkWF0/dFzt1e4LPhft4OzAC+YmadZraO4EfgqZF5j5c0w8y6zOyRQeIYy2eUjeuBfwAuMLPmcNwHgH8ys5ZwmzcAb1V2RZkvmtm3zKwf+B4wB5glaQFwPvCPZnbYzJ4AbiH4oQXBfv6rme02szbg/49M+18E5/QzZnYwjGc4Kf+fZvYrM9togQeAe4BXZLEeyO4czglPUKneZGZTzGyRmf1d+I82l+DKAwAzGyD4pTwvPIGvJzhBdku6PVHckMEiggTRGl4WHyD4JTYzMs+2IWJLiSP0IjBvBMtvt/AnUGT5Y4xwv7LZNoNse7h1ZiOb47IzMtxN8At92HVZUEFmb9q6MhrBMRv0fIrMEz2W0eO0CPho4vwJz6EFZH8csz0OmSwCvhLZ7j6CX99DnX97wy9mgETS2hWZfigSw18T/DJ/TkEx7esHiWPUn1GW/gG40cxaIuMWAT+P7Pt6gh8ks7JYX/KYm1l3OFhHsB/7zKwzMm/0vE0/r6PnwVyOPUeGk/LZSLpE0iNh8eIB4FKCHxDZyOYczglPUMPbQXCCAsGNFIIvhe0AZvZDMzs/nMeAL4SzpjcTv43g1+6MMAlOMbMGMzs5Ms9QTcunxBFamIgji+VbgXlh/NHlMxrBfmWzbQbZ9o5w+CBQG5k2ewTrzea4ZCv9s54MTM92XUMcs6G2kXI+hRZEhqPHaRvBFcmUyKs2vGrMpUzHfBtBcXR02zVm9vAwy2W3QbMXzOwKgh9sXwB+Fh7/dGP6jLLwWuCfJb0lMm4bcEnavldbcK96tHYA0yTVR8ZFz9v08zp6HrRy7DkynORnE97LuoPgymeWmU0hqLWs9HmHiH24czgnPEEN7yfA6xRUEqgEPkqQaB6WdKKkC8MP/DDBL8LEL8ZdwGKFNWbMrJXgMvq/JDWEN16Pk/SqLOO4CzhB0jskVUh6O7AC+GWWy/+J4J7Ah8Ll30xwP+IYI9mvEZgZbrtS0tsI7mXcFU57Arg8nLaS4P5YQhswQIZKH6GxHpeoHwJXSTot3Pd/Bx41sy3DLTjMMYsa9HyKzHOtpPnhjfhPAT8Ox38LuEbSWQpMlvS6tC+5XMh0zG8CPinpZEhW+HlbrjYo6V2SmsJf4wfC0ZmO36g/oyytI7jHc6OkN4bjbgI+myjOlNQk6bKxbMTMthF85p+TVK2gsstfc/TxiR8RJMomBZV6/gVIVFL4CXClpBWSaoFPj3DzVQTF0W1An4KKG6+NTN8FTJfUOMjy2ZzDOeEJahhmtgF4F8HN2D0E9zbeYGY9BB/y58PxOwm+hBM1k34a/t0raW04/B6Ck+NZghv5PyMok84mjr3A6wlOhr0ENyZfb2Z7sly+B3gzwY3b/QTlyHcOMvtI9ysbjwLLwnV+FnhruE8A/x/BDdv9BGXtP4zE3R3O/1BYxJJyv2WsxyVtXb8PY7mD4FfqcaTeIxzKUMcsuo2hzqeEHxL8mNkUvj4TLruG4D7U1wmOVTPB55lTmY65mf2c4MrmdkkdBBU5cvnMz8XAOkldwFeAy83scIbYxvIZZcXMniQ4p74Vfnl/BVgF3COpk6Dywlk52NQVBJUsdgA/Bz5tZr8Np30GWAM8BTwNrOXoefBr4MvAvQTnwL0j2WhYrPghgkSzH3gHwf4lpj9HkCA3hZ//3LTlszmHcyJRm8Q5VwAUPDD+fjP7Xb5jcS7f/ArKOedcQfIE5ZxzriB5EZ9zzrmC5FdQzjnnClJJtWY+Y8YMW7x4cb7DcM45NwKPPfbYHjNrSh9fUglq8eLFrFmzJt9hOOecGwFJGVvD8CI+55xzBckTlHPOuYLkCco551xB8gTlnHOuIHmCcs45V5A8QTnnnCtInqCcy7PDfT14iy7OHauknoNyrtj8cctabl19B3Pqm/jgee9kdv0xzyo6N2H5FZRzedLb38dPn/oNhrGjczefu+9bvLh/x/ALOjdBeIJyLk8e3fYUHUe6ku87jnTx+fu/xXO7N+UxKucKhyco5/LAzPjtC8f2kH247whfevC7PL59fR6icq6weIJyLg+a925l64GgOK+qvJJ/fNX7mVLdAEDvQB9fe/gHPLj5sXyG6FzeeYJyLg+iV09nL3wpy2cu5VMXXs3MuukAGMata+7gNxsezFeIzuWdJyjnxtm+7nYe274u+f7Vx58LQNPkaXzqLz7AgsY5yWk/furXnqTchOUJyrlxdu/GRxmwAQCWNy1lwZTZyWmN1XV84oL3c8KMxclxv3j29/QN9I93mM7lnSco58ZRT38vD2z6c/L9a5ade8w8tVU1fPSVVzGtphGAI309bN7XMm4xOlcoPEE5N44e2fokXT3dAMyoncppc5dnnK+qvJKTZy1Lvl+/e+O4xOdcIfEE5dw4MTN+F6kccdHxZ1Omwf8FV8w6Ljm83p+NchOQJyjnxsnze7awrX0nAFXlVbxiycoh51/etCQ53Lz3RY709cQan3OFxhOUc+Pkty/8KTl83qLTmVxVM+T8U2oamFs/E4C+gX6a926NNT7nCo0nKOfGwZ6D+1kbrVq+7JysljspUsz37C6/D+UmFk9Qzo2Dezc+ghF0qXHyzOOZ2zAzq+VOalqaHH6uze9DuYnFE5RzMTvc18MDm9Yk32eqWj6Y5TOXIgTA5n0tdPccynl8zhUqT1DOxezJHevp7g0Sy8y66Zw658Ssl51cVcOiqXOBoPmjDW1b4gjRuYLkCcq5mDXv3ZYcPmfhaUga0fIrZkbuQ/nzUG4CiTVBSbpY0gZJzZI+kWH6ckl/knRE0sci4xdIuk/SeknrJH04zjidi9PmfUcT1NJpC0a8/PKZkftQnqDcBBJbgpJUDtwIXAKsAK6QtCJttn3Ah4Avpo3vAz5qZicBZwPXZljWuYLXN9CffPYJYHFYXDcSy2YspqKsHICWjl20H+4aZgnnSkOcV1BnAs1mtsnMeoDbgcuiM5jZbjNbDfSmjW81s7XhcCewHpgXY6zOxWJ7+y56+oPTe0btVBqq60a8juqKqpQrL7+KchNFnAlqHrAt8r6FUSQZSYuB04FHB5l+taQ1kta0tbWNJk7nYrNl//bk8OJpo/+NddJMb/bITTxxJqhMd4JtRCuQ6oA7gOvNrCPTPGZ2s5mtNLOVTU1NowjTufhEWyFfPHUsCerofSivKOEmijgTVAsQvSM8H9iR7cKSKgmS021mdmeOY3NuXGyOXEEtnTZ/1Os5bvpCqsqrAGg7uI89B/ePOTbnCl2cCWo1sEzSEklVwOXAqmwWVFAP99vAejP7UowxOhebnv5etrfvSr5fNGXkFSQSKsrKUzox9KsoNxHElqDMrA+4DriboJLDT8xsnaRrJF0DIGm2pBbgI8A/S2qR1ACcB7wbuFDSE+Hr0rhidS4O2w600m9BT7iz65uoHaZx2OGclFLd3O9DudJXEefKzewu4K60cTdFhncSFP2l+yOZ72E5VzQ2Re4/LRnD/aeE9Ad2zWzED/06V0y8JQnnYvJijmrwJSycOofayuAqrP1wJ62dXmvVlTZPUM7FZPO+owlqydTRV5BIKFNZSieGXt3clTpPUM7F4HDvkeQVTpnKWDhlTk7Wm/o8lFeUcKXNE5RzMXjxwI5k/0/zGmYyqaIqJ+tdEenA8LndmxiwgZys17lC5AnKuRikPqA79uK9hDn1TTRW1wNwsPcQW/e35mzdzhUaT1DOxSClBl8OKkgkSEqpbr6+zYv5XOnyBOVcDKJt8OWigkRU9D7Uxr3bhpjTueLmCcq5HOs60k3bwX1A0ALE/Cmzc7r+aIsU0a48nCs1nqCcy7HN+48W7y1onJPsyylX5jQ0UabgX7etax+H+3pyun7nCoUnKOdyLOX5pxzef0qoKq9kdt0MAAxju19FuRLlCcq5HNuyP1pBYuRdvGcjWmwYbZDWuVLiCcq5HIteQY2mi/dszG88mqC2eYJyJcoTlHM5dOBQBwcOB31rTqqoYm7DzFi2s6BxVnK4xYv4XInyBOVcDkWvnhZNmZuszJBr0SuolvadmI2os2rnioInKOdyaHPK/afcPv8UNb12CtUVkwDo6umm/XBnbNtyLl88QTmXQ6ktmOe+Bl+CpJSrqK0HvJjPlR5PUM7liJml1OBbHOMVFMCCtGI+50qNJyjncmRP9366eroBqK2sYebkabFub36kosT2Dq/J50qPJyjnciS9eC/u7tijz0JtO+CtmrvS4wnKuRxJ6WIj5uI9gPkNR6+gWjvb6Bvoj32bzo0nT1DO5UhqC+bxVZBIqK2qYXrtFAD6BvrZGfbg61ypiDVBSbpY0gZJzZI+kWH6ckl/knRE0sdGsqxzhSSoIBFtgy/+KyhIfR7KmzxypSa2BCWpHLgRuARYAVwhaUXabPuADwFfHMWyzhWMvd0HONx3BID6qslMrWkYl+0u8KrmroTFeQV1JtBsZpvMrAe4HbgsOoOZ7Taz1UDvSJd1rpBEr17mNs6MvYJEgtfkc6UszgQ1D4h299kSjsvpspKulrRG0pq2Ni+Dd/mxo3N3cnheTO3vZbJgypzksNfkc6UmzgSV6Sdktg2GZb2smd1sZivNbGVTU1PWwTmXS9vbjyaouZHadXGbWTc92SHivkPtHOw5NG7bdi5ucSaoFiDaGc58YMc4LOvcuNvRkZ8rqIqy8pSE6C1KuFISZ4JaDSyTtERSFXA5sGoclnVuXJlZShFfXF1sDCa16w2/D+VKR0VcKzazPknXAXcD5cCtZrZO0jXh9JskzQbWAA3AgKTrgRVm1pFp2bhidW4s9nYf4EhfDxDU4GuorhvX7ad3veFcqYgtQQGY2V3AXWnjbooM7yQovstqWecK0fZI8d7cxvG9egJPUK50eUsSzo1Ra57uPyXMTyvi884LXanwBOXcGEXv+8ypH/8E1VhdT33VZAAO9x1hT/f+cY/BuTh4gnJujFojbeDNy0MRn6SUls1bDnhFCVcaPEE5NwZmltKCw3g+AxUVLebb5vehXInwBOXcGOw71J6swVdXVUvDpMl5icMrSrhS5AnKuTFIbUFi/NrgS+fdv7tS5AnKuTHYkVK8N/73n6LbVthC2M7OPfT0p7e/7Fzx8QTl3BikPAOVxwQ1qaKKmXXTATAspekl54qVJyjnxiD1Gaj8VJBISL0P5TX5XPHzBOXcKAU1+ArjCgpg4RS/D+VKiyco50Zp/6GOZC+6kytraBznNvjSRa/gWrx3XVcCPEE5N0otKb3ozspbDb6E6MO6/iyUKwWeoJwbpdY89aI7mJmTpzGpogqAjiNdtB/uynNEzo2NJyjnRin9Gah8k8R877zQlRBPUM6NUqE8AxU1L9Lk0XavyeeKnCco50Yh6EU30khsoSSoyBVUtI1A54qRJyjnRmH/oQ4O9R4GoLayhsbq+jxHFPArKFdKPEE5Nwo7OgqjDb500Su5HZ1t3nmhK2qeoJwbhR157kV3MI3V9UyurAHgUO9h9h/qyHNEzo2eJyjnRiGlD6g8dFI4GEkpxXze5JErZp6gnBuFHQXUBl+6aDw7vKKEK2KxJihJF0vaIKlZ0icyTJekr4bTn5J0RmTa30taJ+kZST+SVB1nrM5ly8zY0VF4NfgSolXe/QrKFbPYEpSkcuBG4BJgBXCFpBVps10CLAtfVwPfCJedB3wIWGlmpwDlwOVxxercSBw41EF37yGgsGrwJcyLFDl6txuumMV5BXUm0Gxmm8ysB7gduCxtnsuA71vgEWCKpDnhtAqgRlIFUAvsiDFW57IWff6pkGrwJaQ+C7Xba/K5ohVngpoHbIu8bwnHDTuPmW0HvghsBVqBdjO7J9NGJF0taY2kNW1tbZlmcS6nos8XzW1oymMkmTVU19EwKWhZvae/hz3d+/MckXOjE2eCyvSzMv2nXMZ5JE0luLpaAswFJkt6V6aNmNnNZrbSzFY2NRXel4UrPYVcQSIheh9qR7sX87niFGeCagEWRN7P59hiusHmeTWw2czazKwXuBM4N8ZYnctaIXVSOJiUquZek88VqTgT1GpgmaQlkqoIKjmsSptnFfCesDbf2QRFea0ERXtnS6pVUMB/EbA+xlidy0pQg68IElSDV5Rwxa8irhWbWZ+k64C7CWrh3Wpm6yRdE06/CbgLuBRoBrqBq8Jpj0r6GbAW6AMeB26OK1bnstV+uDNZg6+mspqpNQ15jigzr2ruSkFsCQrAzO4iSELRcTdFhg24dpBlPw18Os74nBuplOK9+qaCq8GXEL03trOzjQEboEz+XL4rLn7GOjcCKU0cFWgFCYC6SbVMqQ6u7nr6e9ndtS/PETk3cp6gnBuB1mgLEgXUBl8mfh/KFTtPUM6NQPR+TqFWMU+INmK73ROUK0KeoJzLkpmlFPHNbyzsBDU/2qJE+848RuLc6HiCci5L+w61J3vRnVyAbfClm9sYbdXcr6Bc8fEE5VyWUor3GmcVbA2+hLn1R1tWae1so2+gP4/RODdynqCcy1K0Db5CL94DqK2qYVpNIwB9A/3s7tqb54icGxlPUM5lKXr/qdArSCREmzzyihKu2AyZoCTdExn+ZPzhOFe40ov4ikFK77reooQrMsNdQUWbB39bnIE4V8gGbIDWImjFPF20ySO/gnLFZrgE5T2dOQfs7tpL70AfAFOqG6ibVJvniLITvVe23a+gXJEZri2+pZJWEfTblBhOMrM3xhaZcwWkJaWCxOw8RjIycyJXULu69tI30E9FWXkeI3Iue8MlqGgX7V+MMxDnClmx1eBLqK6oYkbtVPZ076ff+tnZ2VZUCdZNbEMmKDN7IDEsqSkc5/2quwkn2ulfsVSQSJjXOCvZ7fv29l2eoFzRGK4WnyR9WtIe4DngeUltkv5lfMJzrjBsT2mDr7AbiU0XrdDhFSVcMRmuksT1wPnAy81suplNBc4CzpP093EH51wh6OnvZVf4kKtQyn2dYhBtdd2bPHLFZLgE9R7gCjPbnBhhZpuAd4XTnCt5rR1Bh38ATXXTqK6oynNEI+NVzV2xGi5BVZrZnvSR4X2oynhCcq6wpLRgXiTPP0XNbZiJCNoN3NW5h57+3jxH5Fx2hktQPaOc5lzJ2N4eeUC3yCpIAFSVV9JUNw0Aw1I6XXSukA2XoF4qqUNSZ/jqSLwHXjIeATqXb9uLuAZfgveu64rRkAnKzMrNrMHM6sNXQ+S9F/G5CaEl0tlfMRbxQWpNvhZvUcIVieGqmVdLul7S1yVdLWm4B3vTl79Y0gZJzZI+kWG6JH01nP6UpDMi06ZI+pmk5yStl3TOSLbtXC4c6j3M3u4DAFSUlTOrfkZ+Axql+d55oStCwxXxfQ9YCTwNXAr8V7YrllQO3AhcAqwArpC0Im22S4Bl4etq4BuRaV8BfmNmy4GXAuuz3bZzuRKt9Tanvqlomwmam/IslF9BueIw3BXRCjN7CYCkbwN/HsG6zwSaw2rpSLqdoOmkZyPzXAZ838wMeCS8apoDHAReCVwJYGY9eKUMlwfRB3TnFtnzT1Gz62dQpjIGbIA9B/dzuK+n6KrLu4lnuCuoZH1UM+sb4brnAdsi71vCcdnMsxRoA74j6XFJt0ianGkjYdHjGklr2tq8dpLLrdQ+oIq3iaDK8gpm1wXFk4al3FdzrlBlW4svUXPv1Eitvo5hllWGcenddww2TwVwBvANMzud4IrqmHtYAGZ2s5mtNLOVTU1NmWZxbtSK/RmoqEVT5yaHX9y/I4+ROJedbGvxJWruVUSGG4ZZdwuwIPJ+PpD+XzHYPC1Ai5k9Go7/GUHCcm5cbS/CXnQHs3jq0QKMLfu35zES57Iz3BXUWKwGlklaIqkKuBxYlTbPKuA9YW2+s4F2M2s1s53ANkknhvNdROq9K+di13G4i44jXQBMqqiiafLUPEc0NqlXUJ6gXOEbUbXxkTCzPknXAXcD5cCtZrZO0jXh9JuAuwhqBzYD3cBVkVV8ELgtTG6b0qY5F7uUB3QbZiFlKpEuHgunzEkOb+/YTU9/L1Xl/jijK1yxJSgAM7uLIAlFx90UGTbg2kGWfYKgirtzeZFSQaLI7z8B1FRWM7u+iZ2dQeO3LQd2snT6guEXdC5P4izic66oFWsvukNZFLmK8vtQrtB5gnJuENGq2MVeQSIhWlHixQNek88VNk9QzmVgZqndvBfxQ7pR0YoSW7yquStwnqCcy2Bv9wGO9AWNl9RV1dJYXZ/niHJj0ZSjCWpHxy56+0f6/L1z48cTlHMZpD//VOw1+BJqq2qYWTcdgL6Bfm9RwhU0T1DOZdDSUVo1+KKiV1HeooQrZJ6gnMugFGvwJXhFCVcsPEE5l0FLCTVxlC61ooRXNXeFyxOUc2n6B/pp7TzaD1SxNxKbLlrE19K+k76B/jxG49zgPEE5l2ZX197kl/a0mkZqq2ryHFFu1U2qpWnyNCCoKLHdu4B3BcoTlHNpNu1rSQ4viLS8UEoWTvGGY13h8wTlXJpNe4/2oXnc9IV5jCQ+i6Mtm3tFCVegPEE5l2bjvqMJaum0+XmMJD6pfUN5gnKFyROUcxFH+nqSD68KsSTyRV5KojX5th1o9YoSriB5gnIuYsv+7QzYAABz6ptKroJEQv2kyUyvnQJA70AfrR1t+Q3IuQw8QTkXsTlSQaLU+0ryLuBdofME5VxEc7SCxLTSTlDRmnyeoFwh8gTlXMTmaAWJCXQFtdVr8rkC5AnKudD+7nb2HWoHoKq8qmT6gBpMtKLE1gOt9HtFCVdgPEE5F4pWL18ybR7lZeV5jCZ+jdV1TK1pBKCnv5fWTq8o4QqLJyjnQtEKEqV+/ylhsfew6wpYrAlK0sWSNkhqlvSJDNMl6avh9KcknZE2vVzS45J+GWeczkH6FVRpPqCbLqXrDU9QrsDElqAklQM3ApcAK4ArJK1Im+0SYFn4uhr4Rtr0DwPr44rRuYT+gX427ztak22iXEEt8pp8roDFeQV1JtBsZpvMrAe4HbgsbZ7LgO9b4BFgiqQ5AJLmA68DbokxRucA2N6xm57+HiBowXxqbWOeIxof6RUlEg8pO1cI4kxQ84Btkfct4bhs5/ky8HFgyP8YSVdLWiNpTVub3+R1oxNtIHbJBLl6AphS08CU6gYAevp72Nm5J88ROXdUnAlKGcZZNvNIej2w28weG24jZnazma00s5VNTU2jidO5lPtPx02fGPefEryHXVeo4kxQLUD0p+h8IP0u7GDznAe8UdIWgqLBCyX9IL5Q3UQXvYJaOoGuoCD1PpRXlHCFJM4EtRpYJmmJpCrgcmBV2jyrgPeEtfnOBtrNrNXMPmlm881scbjcvWb2rhhjdRNYd8+h5DNAZSpLqdk2EaTU5PMWJVwBqYhrxWbWJ+k64G6gHLjVzNZJuiacfhNwF3Ap0Ax0A1fFFY9zg9m8fzsWlj7Pb5zNpIqqPEc0vqLPQm3e10JPfy9V5ZV5jMi5QGwJCsDM7iJIQtFxN0WGDbh2mHXcD9wfQ3jOAaldvE+U6uVRU2sbmV3fxM7ONnr6e2nes5UVs47Ld1jOeUsSzm3cuzU5XOoNxA7m5FnHJ4fX7Xohj5E4d5QnKDehmdmEbOIo3SmRBPX0Tk9QrjB4gnIT2p7u/XQc6QKgtrKG2fUz8hxRfixvWkq5gsZxt7W3cuBQZ54jcs4TlJvgUh7QnToPKdOjeaWvunISx89YmHz/7O7mPEbjXMATlJvQUh/QXTjEnKXv5FnLksPrdnmCcvnnCcpNaKkP6E6sFiTSRe9DPbPzBYJKts7ljycoN2H19velPJg6UbrYGMyiqXOpq6oFoONIFy3tu/IckZvoPEG5CWvbgVb6wm7OZ9ZNp6G6Ls8R5VeZylgRvYra9Xweo3HOE5SbwKL3nyZ68V5CyvNQO/0+lMsvT1Buwkp5QHeCPv+U7pRIRYnn92zhSF9PHqNxE50nKDch9Q/0p9RUWzZ9UR6jKRzTahuZWz8TgN6BPp7fsyW/AbkJzROUm5Ce37OFrp5uAKbWNKb0iTTRnTLbq5u7wuAJyk1Ij21/Njl8+tyTJuwDupmkPA/lzR65PPIE5SYcM+PxSIJ62byT8xhN4TmhaTEVZUGzRy0du9jf3Z7niNxE5QnKTThb9m9n36HgS3dyZQ0nNC3Ob0AFprqiimUzFiffr9u9MX/BuAnNE5SbcKLFe6fNXZ68WnBHnZzSqoQ/D+XywxOUm3DWbl+XHD7Di/cyilY3f3bXRm/2yOWFJyg3obR2tNHa2QZAVXllypWCO2rhlDk0TApa1ujsOZjSJJRz48UTlJtQHotcPb1k9glMqqjKYzSFS1JaL7te3dyNP09QbkJZG7n/dMa8FXmMpPCdnNa6uXPjzROUmzD2dbezeX/QvXu5ynnpnOV5jqiwRRNU894XOezNHrlxFmuCknSxpA2SmiV9IsN0SfpqOP0pSWeE4xdIuk/SeknrJH04zjjdxBC9elretITJVTV5jKbwTalpYEHjbAD6BvrZsHtTniNyE01sCUpSOXAjcAmwArhCUnqZyiXAsvB1NfCNcHwf8FEzOwk4G7g2w7LOjcjaHV68N1LRViWe8urmbpzFeQV1JtBsZpvMrAe4HbgsbZ7LgO9b4BFgiqQ5ZtZqZmsBzKwTWA/MizFWV+K6jnSzoW1z8v3pc0/KYzTF4yWzT0gOP7L1SS/mc+MqzgQ1D9gWed/CsUlm2HkkLQZOBx7NtBFJV0taI2lNW1vbWGN2JeqJ1ucYsAEAjpu2gKm1jXmOqDgsn7mEmXXTAejuPcTDW9bmOSI3kcSZoDK1vpn+tN+Q80iqA+4ArjezjkwbMbObzWylma1samoadbCutEUfzj3di/eyVqYyXn38Ocn397zwsD+068ZNnAmqBYj2AjcfSH/ab9B5JFUSJKfbzOzOGON0Je5wX0/Kczwv8wQ1Iq9YspLayqBCya6uPTzVuiHPEbmJIs4EtRpYJmmJpCrgcmBV2jyrgPeEtfnOBtrNrFVB3wffBtab2ZdijNFNAM/sfJ6e/l4A5tbPZHa9X2mPRHVFFa9csjL5/u4XHspjNG4iiS1BmVkfcB1wN0Elh5+Y2TpJ10i6JpztLmAT0Ax8C/i7cPx5wLuBCyU9Eb4ujStWV9oe37E+Oey190bnouPPRmGJ/PrdG9l2YGeeI3ITQUWcKzezuwiSUHTcTZFhA67NsNwfyXx/yrkR6Rvo54kdzyXfe4IanRmTp7Jy/imsbnkagHteeIi/fvlb8hyVK3XekoQraau3PU137yEAptU0sniqP60wWq894bzk8CNbn6D9cFceo3ETgScoV7L6Bvr5+brfJd+fv/hl3rX7GBw3bQFLpwV1mvoG+rlvY8YnP5zLGU9QrmTdv/HPtB3cB0BdVS1/GbkCcCMnidcuOzf5/r6NjyYrnzgXB09QriQd7uvhf9bfl3z/uuWvotbb3huzl80/hWk1wUPOHUe6eHTrU3mOyJUyT1CuJN3z/B/pOBLcI5lW08iFx5+d54hKQ0VZORdFHtz97QsP+YO7LjaeoFzJ6TxykF9veDD5/rKTL6KqvDKPEZWWVy19OVXlQUeP29p3st5bOXcx8QTlSs6vnnuAw31HAJhT38R5i07Pc0SlZXJVDecvPnpM7/EHd11MPEG5krK3+wC/b/5T8v1bTnkt5WXleYyoNL0mUlniydbn2NnpDTW73PME5UrKL9b9nr6BfgCWTlvgD+bGZHZ9U0qPxN9efSe9/X15jMiVIk9QrmTs6NjNQ5HuIN76ktf6c08xev1JFySbP2re+yI/eHyVV5hwOeUJypWMO56+Bwt7azll1gmcNPO4PEdU2o6fvpD/derFyfd/2LzGH951OeUJypWEjXu3pnTp/paXvCaP0Uwcf3nC+Zyz8LTk+x8+8auUnoudGwtPUK7o7eray//+04+S789ccKq3uTdOJHHlyr9KHu9+6+fGh3/InoP78xyZKwWeoFxR29nZxufv+xb7DrUDUFVeyZtP8aun8VRVXskHz30XDZPqAOjsOcjXHr6Nw309eY7MFTtPUK5o7ejYzefu+xYHDncAwRflh857N7Pqpuc5solnWm0j1577DirCKv1bD+zg1tV3eKUJNyaeoFxR2nZgJ1+4/5Zkc0aTKqr4+/Pfy8mzjs9zZBPXCTMW867T35h8v7rlaX713AN5jMgVO09Qruhs2b+d/3ggNTl95BVXsnzm0jxH5l619OVceNxZyfd3PHMPt/z5Zxw41JnHqFyx8gTlisrGvVv5zwdupaunG4Caymr+4ZXv44QZi/MbmEu64rTXc2LTkuT7h15cyyd/8yV+s+HB5EPUzmVDpVRGvHLlSluzZk2+w3AxaN67ld9s+CNrt69LPutUW1nDx155FUumzc9zdC5d15FuvrPmzpSq/xC0QHHFS1/HqXNOyFNkrhBJeszMVh4z3hOUK1QDNsBj25/l7g0PsnHftpRpdVW1fOyV72PR1Ll5is5lY92uZn74+C/Z0bk7ZfxL5yznzae8hgWNs721D+cJyhWHw71HaOnYRfOeF7l346PJHnGjTpm1jCtOex1zG2bmIUI3Un0D/dzb/Ai/ePb3HOo9nDKtYVIdJzQt5sQZSzixaQnzG2d5wpqA8pKgJF0MfAUoB24xs8+nTVc4/VKgG7jSzNZms2wmnqAKX99APwd7DnGo9xAHew5x4FAnLe072da+i23trezu2ptxuYqycs5a8FL+8oTzWTBl9jhH7XKh43AXdzxzDw9ufixZTJuurqqWZTMWMad+Jo3VdTRW19FQXUdjdT0Nk+qYXFXjCawEjXuCklQOPA+8BmgBVgNXmNmzkXkuBT5IkKDOAr5iZmdls2wmo01QOzp2c+czvxvxcqXv6LkRPU8Mwyzx15Lvwei3AXr7+xiwAfoG+ukf6KffBjjS18PBnkP09I/s4c26qlouWHomFx1/NlNqGnK0Xy6fNu9r4VfP/YHndm/kYO+hES1brnImVVRRUVZOZXkFFWUVVIV/y8vKKS8ro0xllCGk4FWmMoRIz2sK50m8c6P3yiUrx3RfcbAEVTGmqIZ2JtBsZpvCAG4HLgOiSeYy4PsWfPs9ImmKpDnA4iyWzZnOIwd5bPszcazajVCZyphT38T8xlksb1rK2YtOo7qiKt9huRxaMm0+1537DsyMlvZdbGjbzIY9m3m+bUvy0YHB9Fs/3SNMai5+J8X0iEecCWoeEL2z3UJwlTTcPPOyXBYASVcDVwMsXLhwbBG72JWpjNrKamqraphcWcPkqlrmNjSxYMocFjTOZk5Dk3fPPkFIYsGU2SyYMptXLzsHM2Nn5x5e2PsiBw510H64i44jXcHfcDj9HpYrbXEmqEzXzOnliYPNk82ywUizm4GbISjiG0mACXPqm/i7c94xmkVLQraFG4rMWaYypHCclJxWWV5OmcqoKCsPilxUTmV5OZXllUyurGFSRZXfQ3AZSWJOQxNzGpoGnae3v4+e/l56+/voGwhePf199PX30TsQFC0nip8HbIABMzAL/kLy3leiaHo4pVOFLF4Lp8yJZb1xJqgWYEHk/XxgR5bzVGWxbM40VNfx8vmnxLV651yOVJZXUFke59eWKyRxtiSxGlgmaYmkKuByYFXaPKuA9yhwNtBuZq1ZLuucc66ExfZTxMz6JF0H3E1QVfxWM1sn6Zpw+k3AXQQ1+JoJqplfNdSyccXqnHOu8PiDus455/JqsGrm3lisc865guQJyjnnXEHyBOWcc64geYJyzjlXkEqqkoSkNuDFfMdRQGYAe/IdRBHy4zY6ftxGx48bLDKzY57QLqkE5VJJWpOpZowbmh+30fHjNjp+3AbnRXzOOecKkico55xzBckTVGm7Od8BFCk/bqPjx210/LgNwu9BOeecK0h+BeWcc64geYJyzjlXkDxBFSFJF0vaIKlZ0ieGmO/lkvolvXWky5ai0R43SQsk3SdpvaR1kj48flHn31jOt3B8uaTHJf0y/mgLxxj/T6dI+pmk58Lz7pzxibqweIIqMpLKgRuBS4AVwBWSVgwy3xcIuiwZ0bKlaCzHDegDPmpmJwFnA9f6ccs4X/pxS/gwsD7OOAtNDo7bV4DfmNly4KVMsOOX4Amq+JwJNJvZJjPrAW4HLssw3weBO4Ddo1i2FI36uJlZq5mtDYc7Cb4s5sUfckEYy/mGpPnA64Bb4g60wIz6uElqAF4JfBvAzHrM7EDsERcgT1DFZx6wLfK+hbQvS0nzgL8CbhrpsiVsLMctOs9i4HTg0dyHWJDGety+DHwcGIgpvkI1luO2FGgDvhMWjd4iaXKcwRYqT1DFRxnGpT8r8GXgH82sfxTLlqqxHLdgBVIdwa/d682sI7fhFaxRHzdJrwd2m9ljMcVWyMZyvlUAZwDfMLPTgYPAhLpfnBBbl+8uNi3Agsj7+cCOtHlWArdLgqAhyksl9WW5bKka9XEzs19IqiRITreZ2Z3jEXCBGMv5dhbwRkmXAtVAg6QfmNm74g8778Zy3B4BWswscZX+MyZogsLM/FVEL4IfFZuAJUAV8CRw8hDzfxd462iWLaXXGI+bgO8DX873fhTTcUsbfwHwy3zvT7EcN+BB4MRw+AbgP/O9T/l4+RVUkTGzPknXEdT6KQduNbN1kq4Jpw96/2SwZccj7nwby3EDzgPeDTwt6Ylw3KfM7K44Yy4EYzxuE1YOjtsHgdskVREkuqtiDbhAeVNHzjnnCpJXknDOOVeQPEE555wrSJ6gnHPOFSRPUM455wqSJyjnnHMFyROUy4qkrgzjrpH0nnD4fkkrs1zXYknPhMMrJX01HL5B0sdyGXcu1ivptPBh08T7N+aiJXhJkyT9TtITkt4+1vUNso1jPrcxrCvjZyzpSklfD4ej58SVkuZG5rsl143sRrftSo8/B+VGLRfPwJjZGmBNDsIZE0kVZtY3yOTTCJ76vwvAzFYBq3Kw2dOBSjM7LQfrKghp58SVwDOELSiY2fvzEVOuDHOOuBj4FZQbtUxXJpLKJH1P0mfCfoD+U9JqSU9J+kCGdVyQ1k/QivCX+iZJH4rM9xFJz4Sv67MY/09hXzy/A04cJP7vSvqSpPuAL0g6U9LDYQOdD0s6MXxQ8l+BtyeudNKuGBZJ+n24f7+XtDDDdqZJ+kU4zyOSTpU0E/gBcFq43uPSljlO0m8kPSbpQUnLw/FvkPRoGOPvJM0Kx9dJ+o6kp8PtvCWyrs9KejLc9qwM8R2z3+H4Gkm3h+v7MVATWeYqSc9LeoDgQebE+BskfUxB30YrCR42fSJcV/IKTNIVYazPSPpCZPmuTPEOtt+DfK5lkl6Q1BR53yxphqQmSXeE5+RqSecNcwyulPRTSf8D3DPYNl1M8t2Uhb+K4wV0ZRh3A/CxcPh+gr6SfgT8UzjuauCfw+FJBFdKS4DFwDPh+AsIm8AJ1/dwOO8MYC9QCbwMeBqYDNQB6wiuPoYbXws0AM2JONPi/y7wS6A8fN8AVITDrwbuCIevBL4eWS75Hvgf4L3h8PuAX2TYzteAT4fDFwJPpO97hmV+DywLh88C7g2Hp3L0Afv3A/8VDn+BSFNMwNTwrwFvCIf/I/F5pG1rsP3+CEELCACnEvSLtRKYA2wFmgia8XkocjzSz4mVke3cHy4/N7J8BXAv8Kah4h1iv1M+m8i2Pk3QqC/AayP79EPg/HB4IbA+i8++BZiW7//BifjyIj6XS98EfmJmnw3fvxY4VUd7Cm0ElgHPD7GOX5nZEeCIpN3ALOB84OdmdhBA0p3AKwjayMs0viwc3x2OH6o47qd2tDXpRuB7kpYRfFFWZrHP5wBvDof/D8GXarrzgbcAmNm9kqZLahxshQpaTT8X+KmUbBR7Uvh3PvBjSXMIksPmcPyrgcsTM5vZ/nCwhyAJAzwGvCbDJgfb71cCXw3X95Skp8LxZwH3m1lbGO+PgRMG258MXp62/G3htn4xRLyD7fdgbgX+L0GL4e8DvhOOfzXBVXpivgZJ9Qz92f/WzPaNYP9cjngRn8ulh4G/kFQdvhfwQTM7LXwtMbPhikmORIb7CX5hZ+q6gCHGQ/bdiByMDP8bcJ+ZnQK8gaAF7pHKtN2RdnNSBhyIHLfTLOjNF4Krsa+b2UuAD0Ri1CDr7LXwUoCjxzPdUPs9WJxjaSNtqM9tsHgH2+/MwZltA3ZJupAgof46nFQGnBM5rvMs6IRyqGMQPUfcOPIE5XLp2wQVCX4qqYKgocy/VdBVBZJO0Og6XvsD8CZJteHyf0XQ2vNQ4/8qvO9RT/CFk41GYHs4fGVkfCdQP8gyD3P0yuWdwB8Hif+dENxzA/bYEP1JhdM2S3pbuIwkvTRDjO+NLHYPcF3ijaSpg60/g8H2Oxr3KQTFfBB01nhBeCVYCbxtkPUOdtweBV4V3hMqB64AHhhBjO8dasaIWwju8/0kcpWcfpxOy7D+K7Ncv4uZJyiXrVpJLZHXRzLNZGZfAtYSFHfdAjwLrFVQrfybjKLmqAXdrX8X+DPBl9stZvb4MON/DDxB0IfTg1lu6j+Az0l6iKAF6oT7CIqFMlUH/xBwVVj89W7gwxnWewOwMpzn82T3BftO4K8lPUlwby3RXfgNBD8AHgT2ROb/DDA1rHTwJPAXWWwjYbD9/gZQF8b9cYLjjJm1hnH8CfgdweedyXeBmxKVJBIjw+U/SXBcnwTWmtn/HSbGG8i830NZRXBv8juRcR8i/CwkPQtcE44f7Bi4PPLWzJ1zJSmsMfjfZvaKfMfiRscrSTjnSo6CB6n/lrCI0hUnv4JyzjlXkPwelHPOuYLkCco551xB8gTlnHOuIHmCcs45V5A8QTnnnCtI/w8Jj0OKUa/JtAAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"marginal_lr.plot(color='C2')\n",
"\n",
"decorate(xlabel='Likelihood ratio of each additional year',\n",
" ylabel='PDF',\n",
" title='Posterior distribution of slope in terms of likelihood ratio')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The mean of this distribution is about 0.43, which indicates that each additional year is evidence that the respondent will say \"yes\", with a a likelihood ratio (or Bayes factor) of 0.43."
]
},
{
"cell_type": "code",
"execution_count": 94,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.42780656870625144"
]
},
"execution_count": 94,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"marginal_lr.mean()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Later we will use the joint posterior distribution to generate predictions, but first I'll show how to compute likelihoods under a log transform."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Log Likelihood\n",
"\n",
"Because of the problem of underflow, many likelihood computations are done under a log transform. That's why the distributions in SciPy, including `binom`, provide functions to compute logarithms of PMFs and PDFs.\n",
"\n",
"Here's a loop that uses `binom.logpmf` to compute the log likelihood of the data for each pair of parameters in `joint_pmf`: "
]
},
{
"cell_type": "code",
"execution_count": 95,
"metadata": {},
"outputs": [],
"source": [
"log_likelihood = joint_pmf.copy()\n",
"\n",
"for slope, inter in joint_pmf.index:\n",
" ps = expit(inter + slope * xs)\n",
" log_likes = binom.logpmf(ks, ns, ps)\n",
" log_likelihood[slope, inter] = log_likes.sum()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`log_likes` is an array that contains the logarithms of the binomial PMFs for each group.\n",
"The sum of these logarithms is the log of their product, which is the log-likelihood of the data.\n",
"\n",
"Since the likelihoods are small, their logarithms are negative. The smallest (most negative) is about -610; the largest (least negative) is about -480."
]
},
{
"cell_type": "code",
"execution_count": 96,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-607.8199099184933, -480.8081879180253)"
]
},
"execution_count": 96,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"log_likelihood.min(), log_likelihood.max()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So the log likelihoods are comfortably with the range we can represent with floating-point numbers.\n",
"\n",
"However, before we can do the update, we have to convert the logarithms back to a linear scale. To do that while minimizing underflow, I am going to shift the logs up toward zero.\n",
"\n",
"Adding a constant to the `log_likelihood` is the same as multiplying a constant by `likelihood`. \n",
"We can do that without affecting the results because we have to normalize the posterior probabilities, so the multiplicative constant gets normalized away."
]
},
{
"cell_type": "code",
"execution_count": 97,
"metadata": {},
"outputs": [],
"source": [
"shifted = log_likelihood - log_likelihood.max()\n",
"likelihood2 = np.exp(shifted)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"After subtracting away the largest element in `log_likelihood`, the range of values in the result is from -127 to 0."
]
},
{
"cell_type": "code",
"execution_count": 98,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-127.01172200046801, 0.0)"
]
},
"execution_count": 98,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"shifted.min(), shifted.max()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So the range of likelihoods is from near 0 to 1."
]
},
{
"cell_type": "code",
"execution_count": 99,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(6.910508383385572e-56, 1.0)"
]
},
"execution_count": 99,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"likelihood2.min(), likelihood2.max()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can use them as likelihoods in a Bayesian update."
]
},
{
"cell_type": "code",
"execution_count": 100,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.021604219991658798"
]
},
"execution_count": 100,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior_pmf2 = joint_pmf * likelihood2\n",
"posterior_pmf2.normalize()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"To confirm that we get the same results using likelihoods or log-likelihoods, I'll compute the mean of the marginal posterior distributions:"
]
},
{
"cell_type": "code",
"execution_count": 101,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-0.8491556076374563 0.031694255696289274\n"
]
}
],
"source": [
"joint_posterior2 = posterior_pmf2.unstack()\n",
"\n",
"marginal2_inter = marginal(joint_posterior2, 0)\n",
"marginal2_slope = marginal(joint_posterior2, 1)\n",
"\n",
"print(marginal2_inter.mean(), marginal2_slope.mean())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And compare them to what we got using (non-log) likelihoods."
]
},
{
"cell_type": "code",
"execution_count": 102,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-0.8491556076374563 0.03169425569628929\n"
]
}
],
"source": [
"print(marginal_inter.mean(), marginal_slope.mean())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"They are the same except for small differences due to floating-point approximation.\n",
"\n",
"In this example, we can compute the posterior distribution either way, using likelihoods or log likelihoods.\n",
"But if there were more data, the likelihoods would underflow and it would be necessary to use log likelihoods."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Making predictions\n",
"\n",
"As we did with the previous example, we can use the posterior distribution of the parameters to generate predictions, which we can use to see whether the model fits the data and to extrapolate beyond the data.\n",
"\n",
"I'll start with a sample from the posterior distribution."
]
},
{
"cell_type": "code",
"execution_count": 103,
"metadata": {
"tags": [
"remove-cell"
]
},
"outputs": [],
"source": [
"np.random.seed(42)"
]
},
{
"cell_type": "code",
"execution_count": 104,
"metadata": {},
"outputs": [],
"source": [
"sample = posterior_pmf.sample(101)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And a range of `xs` that extends 20 years past the observed data."
]
},
{
"cell_type": "code",
"execution_count": 105,
"metadata": {},
"outputs": [],
"source": [
"xs = np.arange(1880, 2021) - offset"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use the sampled parameters to predict probabilities for each group."
]
},
{
"cell_type": "code",
"execution_count": 106,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(101, 141)"
]
},
"execution_count": 106,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ps = np.empty((len(sample), len(xs)))\n",
"\n",
"for i, (slope, inter) in enumerate(sample):\n",
" ps[i] = expit(inter + slope * xs)\n",
" \n",
"ps.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But that only accounts for uncertainty about the parameters.\n",
"\n",
"We also have to account for variability in the size of the groups. Here's the distribution of group size, dropping the groups smaller than 20."
]
},
{
"cell_type": "code",
"execution_count": 107,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"count 104.000000\n",
"mean 358.567308\n",
"std 232.895275\n",
"min 29.000000\n",
"25% 162.500000\n",
"50% 336.500000\n",
"75% 492.750000\n",
"max 846.000000\n",
"Name: count, dtype: float64"
]
},
"execution_count": 107,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"not_small = (data['count'] >= 20)\n",
"counts = data.loc[not_small, 'count']\n",
"counts.describe()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To simulate variation in group size, I'll use `np.random.choice` to resample the group sizes; that is, I'll draw from `counts` a sample with the same length as `xs`, sampling with replacement."
]
},
{
"cell_type": "code",
"execution_count": 108,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([718, 330, 842, 31, 31, 197, 782, 260, 78, 155])"
]
},
"execution_count": 108,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ns = np.random.choice(counts, len(xs), replace=True)\n",
"ns[:10]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Even if we know how many people are in each group and their probability of saying \"yes\", there is still uncertainty in the outcome. We can use the binomial distribution to simulate this (final) source of uncertainty.\n",
"\n",
"Putting it all together, the following loop combines these sources of uncertainty to generate predictive distributions for each group."
]
},
{
"cell_type": "code",
"execution_count": 109,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(101, 141)"
]
},
"execution_count": 109,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pred = np.empty((len(sample), len(xs)))\n",
"\n",
"for i, (slope, inter) in enumerate(sample):\n",
" ps = expit(inter + slope * xs)\n",
" ns = np.random.choice(counts, len(xs), replace=True)\n",
" ks = binom(ns, ps).rvs(len(xs))\n",
" pred[i] = ks / ns\n",
" \n",
"pred.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The result is an array with one row for each pair of parameters in the sample and one column for each value in `xs`.\n",
"\n",
"Now we can use `np.percentile` to compute percentiles in each column."
]
},
{
"cell_type": "code",
"execution_count": 110,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(141,)"
]
},
"execution_count": 110,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"low, median, high = np.percentile(pred, [5, 50, 95], axis=0)\n",
"median.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And use them to plot a 90% credible interval for the predictions."
]
},
{
"cell_type": "code",
"execution_count": 111,
"metadata": {
"scrolled": true,
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAByKElEQVR4nO29eXxcZ3no/33OrNLMSBrtmxdZsp14i504G4EkdgiEJWwtNxDIAi0tbYFC7+Ve2tL99l5K+2tpC7cQKM0CZSmUlkIaCraTEEhInDjeY1teZMuSrF2a0exznt8fZzQe7ZJj2ZL8fj+f+eic97znPc85mplnnud93ucRVcVgMBgMhoWGdbkFMBgMBoNhMoyCMhgMBsOCxCgog8FgMCxIjIIyGAwGw4LEKCiDwWAwLEiMgjIYDAbDgsQoKMOiQESKROQ/RGRIRP7lIox3u4i0XwzZZrjOgyLyzHxfZ5rrHxSR23PbvyciX7lcsixEROSUiLz+csthmByjoJYwIvJaEfl57ku9X0R+JiLXX265xjNLZfHLQA1QoarvvgRiLQlUdb2qPpnb/j+q+quXWaQljYioiLRcbjmWCu7LLYBhfhCREuAHwG8A3wa8wOuA5OWUazwiMtv34ArgqKpmLuQaF3LeYuZKvOfLiXne84OxoJYuawBU9RuqmlXVuKr+l6ruAxCRPxaRr412FpGVuV9/7tz+kyLyf0Xk+ZwF9u8iUj6u76+JSIeIdIrIfy8Yyycin8sd68ht+3LHbheRdhH5XyLSBXwD+E+gXkSiuVd94Y2IyJ8Afwjckzv+KyJiicinRaRNRLpF5FERKR0n36+IyGlg50wPS0TqReS7ItIjIidF5GMFx4pE5BERGRCRwyLyPwstPhH5lIgcF5GIiBwSkXfO5h9UIOcHRORMbvwPi8j1IrJPRAZF5PMF/ZtFZKeI9IlIr4h8XUTKCo6fyj3XfcCIiLgLXViF//PJrNZxfW8QkWdzMnSKyOdFxFvQV3OyHsvJ/QURkdnIOe6aXxSRvxrX9u8i8ju57f8lImdzz/aIiNwxxThFIvL/5d4PQyLyjIgU5Y69TRxX52DufX31uNM35573kIh8S0T8BeN+SERaxfFAfL/wvZl7Br8lIseAYyLydO7Q3tz79J7JZDXMAVU1ryX4AkqAPuAR4E1AeNzxPwa+VrC/ElDAndt/EjgLbAACwHdH+xf0/Ubu2EagB3h97vifAs8B1UAV8HPgz3LHbgcywF8APqAo19Y+w/2Ml/eDQCuwCggC/wo8Nk6+R3PyFU0yXv6aOD/UXsRRgt7cmCeAN+aOfwZ4CggDjcC+QnmBdwP1uXHuAUaAutyxB4FnprinUTm/CPiBNwAJ4N9yz64B6AZuy/VvAe7MPbcq4GngcwXjnQJeBpaN3nOu7fXjn+Fkz3xc3+uAm3C8LCuBw8DHC/oqjoVeBizP/f/vmo2c4655K3AGkNx+GIjnnufa3LH6gufVPMU4X8B5zzYALuA1ueuvyf0/7gQ8wP/Eed94C+75+dz1ynP3+eHcse1AL3Btbqy/B54e9wx+nDuvqKCt5XJ//pfK67ILYF7z+M+Fq4GHgXYcpfB9oCZ3LP9lldtfyUQF9ZmC4+uAVO7DP9r3qoLjnwX+Mbd9HHhzwbE3Aqdy27fnxvEXHJ/wZTnJvYyXdwfwmwX7a4E0579QFVg1zXj5awI3AqfHHf9d4J9y23llldv/1enkxVESb89tP8jMCqqhoK0PuKdg/7sUKIZx578D2FOwfwr44Lg+p7gABTXJtT4OfK9gX4HXFux/G/jUbOQcd0yA08Ctuf0PATtz2y04Cvr1gGea523hKLVrJjn2B8C3x/U9C9xecM/vH/c+/mJu+x+BzxYcC+beYysLnsH2cdczCuoivoyLbwmjqodV9UFVbcSxhOqBz81hiDMF2204v0Arpzk+6v6oz+1PdgygR1UTc5BjMia7hhsnkGIy+aZjBY6LcXD0BfxewVj148YaM66I3C8iLxecu4Gxz2kmzhVsxyfZD+auUy0i38y5vIaBr01yndne87SIyBoR+YGIdOWu9X8muVZXwXZsjnICoM63+jeB9+aa7gW+njvWiqMY/xjozo1ZP8kwlTgW6PFJjo15n6iqjfOMGma6j0nOjeL8gCg896I8b8PkGAV1haCqr+BYUxtyTSNAcUGX2klOW1awvRzn12PvNMc7ctsdOF/6kx0D51cm0+zPhsmukWHsl/tsxz0DnFTVsoJXSFXfnDveiePaGyV/3yKyAvgy8BGcCMMy4ACOZXCx+b8497RJVUuA909yndne85j/v4i4cNxxo/wD8AqwOnet35vkWq9GzkK+Afxy7lneiGM1AqCq/6yqr8X5XyuOa3g8vTiu0eZJjo15n+TmyZbhWFEzMf7cAFAx7lxTDmIeMQpqiSIiV4nIfxeRxtz+Mpxfqc/lurwM3Coiy8UJLvjdSYZ5v4isE5FinHml76hqtuD4H4hIsYisBz4AfCvX/g3g0yJSJSKVOHM7X2NqzgEVOTlmyzeAT4hIk4gEcX7hf0svLJLqeWA4NyFfJCIuEdkg50Pyvw38roiERaQBRxmNEsD5kuoBEJEPcP5HwMUmBESBwZwcn3wVYx0F/CLyFhHxAJ/GmWcpvNYwEBWRq3CiQedFTlXdg/P8vgL8SFUHAURkrYhsFyfAJoFjTWYnOd8Gvgr8tTjBLi4RuTl33reBt4jIHbn7/O84kaw/n8V9/DPwARHZnBvr/wC/UNVT05xzDmcO03ARMApq6RLB+TX6CxEZwVFMB3A+oKjqj3EUyj6cAIEfTDLGYzhWVxeOC+Vj444/hTPhvAP4K1X9r1z7/wZ258beD7yUa5uUnHX3DeBEzk02mRtnPF/Nyfc0cBLnC+yjszhvsutngbuBzbmxenG+LEcV5p/izOOdBH4CfIdcuL6qHgL+P+BZnC+njcDPLkSOWfAnOBP2Q8APcQJDLghVHQJ+E+c+z+JYVIVRff8Dx90WwbEQvzV+jIss5zdw5pr+uaDNhxOg0ovzHqzGseQm43/gvNdeAPpxLC1LVY/gWHB/nxvnbuBuVU3NJJCq7sCZw/oujhXdDLxnhtP+GHgk9z7+bzNdwzA9o5EzBsMYRORJnAn1CZkHRGQlzpe15wItlkWNiPwG8B5Vve1yyzIXRORPgUZV/eDllsVgmA3GgjIYZkBE6kTkFnHWXq3FsUK/d7nlmgu5uZd1OD8sDIZFgckkYTDMjBf4EtAEDOJEnf2/yynQBfASjlvyIzN1NBgWCsbFZzAYDIYFiXHxGQwGg2FBsuhcfJWVlbpy5crLLYbBYDAYLhIvvvhir6pWjW9fdApq5cqV7N69+3KLYTAYDIaLhIi0TdZuXHwGg8FgWJAYBWUwGAyGBYlRUAaDwWBYkCy6OajJSKfTtLe3k0i82gTZhvnG7/fT2NiIx+O53KIYDIYFzpJQUO3t7YRCIVauXImzYN6wEFFV+vr6aG9vp6mp6XKLYzAYFjhLQkElEgmjnBYBIkJFRQU9PT2XWxSDYV5obe9nx542Onsj1FWGuGPLCloayy+3WIuWJTMHZZTT4sD8nwxLldb2fh5+Yh/RWJKGqhDRWJKHn9hHa3v/5RZt0bJkFJTBYDBcTnbsaSMc8lMW9GOJUBb0Ew752bFn0iU+hllgFNRFIhgMztxpBnbv3s3HPja+5NJ5Tp06xT//8z/Puv988uCDD/Kd73znVfcxGJYKnb0RSgK+MW0lAR9dfdHLJNHiZ0nMQS0Vtm7dytatW6c8Pqqg7r333ln1NxgMl466yhDDI0nKgv582/BIktqKV//jdSGg2QwkhsHjB0/RJXHXGwtqHnn55Ze56aab2LRpE+985zsZGBgA4IUXXmDTpk3cfPPNfPKTn2TDBqdC+JNPPslb3/pWAJ566ik2b97M5s2b2bJlC5FIhE996lP89Kc/ZfPmzfzN3/zNmP7RaJQPfOADbNy4kU2bNvHd7353gjwrV67k937v97j55pvZunUrL730Em984xtpbm7mi1/8IuBE2o3KtHHjRr71rW/l2z/ykY+wbt063vKWt9Dd3Z0f98UXX+S2227juuuu441vfCOdnZ3z91ANhgXKHVtWMBBJMBhNYKsyGE0wEElwx5YVl1u0V42qYncfJ3t6L9lTu7FPPIc93D3zia+SJWdBJR/99Xkb23f/l+bU//777+fv//7vue222/jDP/xD/uRP/oTPfe5zfOADH+Chhx7iNa95DZ/61KcmPfev/uqv+MIXvsAtt9xCNBrF7/fzmc98hr/6q7/iBz9wqrM/+eST+f5/9md/RmlpKfv37wfIK8PxLFu2jGeffZZPfOITPPjgg/zsZz8jkUiwfv16PvzhD/Ov//qvvPzyy+zdu5fe3l6uv/56br31Vp599lmOHDnC/v37OXfuHOvWreODH/wg6XSaj370o/z7v/87VVVVfOtb3+L3f//3+epXvzqnZ2UwLHZaGst58K5N7NjTRkdvlNqKIG+/Zc2SiOLTgbPoUAcSqkRE0Gwau+MQOjKAVb0Kcc3PusYlp6AWCkNDQwwODnLbbU5V8AceeIB3v/vdDA4OEolEeM1rXgPAvffem1c4hdxyyy38zu/8Du973/t417veRWNj47TX+8lPfsI3v/nN/H44HJ6039ve9jYANm7cSDQaJRQKEQqF8Pv9DA4O8swzz/De974Xl8tFTU0Nt912Gy+88AJPP/10vr2+vp7t27cDcOTIEQ4cOMCdd94JQDabpa6ubo5Py2BYGrQ0li8JhVSIxgawe1qRQHnerScuDwQrYLgL9fqRivmxEo2CusTMtkDkpz71Kd7ylrfw+OOPc9NNN/GTn/xkxnFn4xP2+ZxJXMuy8tuj+5lMZlr5JhtfVVm/fj3PPvvsjNc2GAwLG03FQG3E58ybaSJKtv0A4i9BLNeYviKCur2obc+bPEtOQc3VDTdflJaWEg6H+elPf8rrXvc6HnvsMW677TbC4TChUIjnnnuOm266aYzVU8jx48fZuHEjGzdu5Nlnn+WVV15h2bJlRCKRSfu/4Q1v4POf/zyf+9znAMfFN5UVNR233norX/rSl3jggQfo7+/n6aef5i//8i/JZDJ86Utf4v7776e7u5tdu3Zx7733snbtWnp6enj22We5+eabSafTHD16lPXr18/52gaD4fKh6QT26b1oNolVvRoJhLHb9yEeP+L2XhaZlpyCulzEYrExbrjf+Z3f4ZFHHuHDH/4wsViMVatW8U//9E8A/OM//iMf+tCHCAQC3H777ZSWlk4Y73Of+xy7du3C5XKxbt063vSmN2FZFm63m2uuuYYHH3yQLVu25Pt/+tOf5rd+67fYsGEDLpeLP/qjP+Jd73rXnO/jne98J88++yzXXHMNIsJnP/tZamtreec738nOnTvZuHEja9asybsuvV4v3/nOd/jYxz7G0NAQmUyGj3/840ZBGQwLFE1GnSi8AotIs2nss4dAQIrD6LlW1BJweRFv0WWTVWbrcloobN26VccXLDx8+DBXX331ZZJo7kSj0fy6qc985jN0dnbyt3/7t5dZqkvHYvt/GQxLBU0nsU+9gJQvw8rNG6ltY3e+ArF+pLjsfN9sesbgB01EoKQWV9Wry60pIi+q6oQ1M/MaZi4id4nIERFpFZEJ4WoiUioi/yEie0XkoIh8YD7lWSj88Ic/ZPPmzWzYsIGf/vSnfPrTn77cIhkMhisAu68NVLF7T6HJEQB04AxEe8YoJ2Bm5ZQcIbP7O+hQ13yJO38uPhFxAV8A7gTagRdE5Puqeqig228Bh1T1bhGpAo6IyNdVNTVfci0E7rnnHu65557LLYbBYLiC0PgwOtiBBCuQ5Ah2z3GscCN2z0kkMLf5as2kyL7wbbT/DNlnvooVqsSqWX3RZZ7POagbgFZVPQEgIt8E3g4UKigFQuKEhwWBfiAzjzIZDAbDFYdmUtjnjiG+gBON6w+ikV7s2CBSFJoQoTftWHaW7IvfRYc6c2MnHVffPDCfCqoBOFOw3w7cOK7P54HvAx1ACLhHVSfELIrIrwG/BrB8+fJ5EdZgMBiWGpocwR7ogOEuQJBA2fmDgTBk04jbN9Xp58cZ6cfuOASJCDrcjQ6ezR9zXfNWXCuuvfjCM78KarJFOeMjMt4IvAxsB5qBH4vIT1V1eMxJqg8BD4ETJHHxRTUYDIalhcaHyJ7Z51hHxaWIjA05EMsF01hOmk6gPSex2/eiPScm7WM134y18vqLKnch86mg2oFlBfuNOJZSIR8APqNOKGGriJwErgKen0e5DAaDYUljR/uwzx5AfEHEM7OFNObcwQ7sI0+ivW1MtCnOYzVdjzSNd4pdXOYziu8FYLWINImIF3gPjjuvkNPAHQAiUgOsBSZX1Qucc+fOce+997Jq1Squu+46br75Zr73ve8Bzhqp973vfWzcuJENGzbw2te+lmjUScH/53/+56xfv55NmzaxefNmfvGLX0x7nVOnTuWTy07Xp7Ash8FgWPpoNoM93E32zF7s9v2IPzQn5aTpJNmD/0X2Zw+jvacYr5ykugVr/RtxbXkH7tf9Kq51d857RvN5s6BUNSMiHwF+BLiAr6rqQRH5cO74F4E/Ax4Wkf04LsH/paq98yXTKBe7LLOq8o53vIMHHnggrxja2tr4/vcdffy3f/u31NTU5BO5HjlyBI/Hw7PPPssPfvADXnrpJXw+H729vaRSrz6AcXxZDoPBsLjRdMLZsFyThn9rOuEstE0OgzeAFaqc09j2qd3YJ1+AdHzMMSmtQ6qbsRo2IIFLn2NwXjNJqOrjwOPj2r5YsN0BvGE+ZRjPaFnmcMhPQ5VTv+XhJ/bx4F2bLlhJ7dy5E6/Xy4c//OF824oVK/joRz8KQGdnJytWnE+muHbt2nx7ZWVlPideZeXkb6oXX3yRD37wgxQXF/Pa1742337q1Cnuu+8+Rkac9Qyf//zn8xnSDx8+zObNm3nggQd45zvfOWk/g8Gw8NFEhOzpPbkdsCpXIuXL8taLJiLY7ftBBAnOXjEB2P3tZHf/y0TFVLkS14a7LotSKuSKS3VUWJYZyP/dsaftghXUwYMHufbaqaNYPvjBD/KGN7yB73znO9xxxx088MADrF69mje84Q386Z/+KWvWrOH1r38999xzTz6FUCEf+MAH8mU7PvnJT+bbq6ur+fGPf4zf7+fYsWO8973vZffu3RPKcsRisUn7GQyGS8uFeG/s4XOI5UaKSlA7i91zAsmmscIN2ANn0YEziDcw55REdl8b2Re+Ddn0+cbiMlxrbkPq112SgoQzccUVLLwUZZl/67d+i2uuuYbrr3eiWzZv3syJEyf45Cc/SX9/P9dffz2HDx8mGAzy4osv8tBDD1FVVcU999zDww8/PGas8WU77rvvvvyxdDrNhz70ITZu3Mi73/1uDh06xGTMtp/BYJg/Rr030ViShqoQ0ZjjvWlt75/yHM1mYLATfAHAibyTYDkMtJM98QsY7HTKYMxBOakqdschss9/67xy8hbjuuZu3Lf9OlbD+gWhnOAKtKDmoyzz+vXrx1Sw/cIXvkBvb++YcuzBYJB3vetdvOtd78KyLB5//HGuvvpqXC4Xt99+O7fffjsbN27kkUce4cEHH8yfN10Zjb/5m7+hpqaGvXv3Yts2fr//VfUzGAzzxwV5b+JDYGfHLKQVsZxaTLMssVOIxofIHvgR2t16vtEXxH3TvXN2D14KrjgLaj7KMm/fvp1EIsE//MM/5NtisVh++2c/+1m+wm0qleLQoUOsWLGCI0eOcOzYsXy/l19+ecxcFUBZWRmlpaU888wzAHz961/PHxsaGqKurg7LsnjsscfIZrMAhEKhMWU5pupnMBguHRfivbEHzuatp/HMRTlpJkX2yFNknvzSWOVUVIr75vdfuHKamFfhonLFWVDzUZZZRPi3f/s3PvGJT/DZz36WqqoqAoEAf/EXfwE4tZ1+4zd+wzGtbZu3vOUt/NIv/RIvvfQSH/3oRxkcHMTtdtPS0sJDDz00Yfx/+qd/ygdJvPGNb8y3/+Zv/ia/9Eu/xL/8y7+wbds2AgHnjbxp06YxZTmm6mcwGC4dc/XeaCqOxvqxXoVlo5kk9umXsU/8ApJjFaG1fAvWVdsQz4V5VDSbRm0bV0nVBcs3E6bchuGSY/5fhiuRwgjikoCP4ZEkA5HEpBHEo8EQDHXNKZGrxgbJHnkSjfQgYqGxQcgkx/SRkhqs9W/EKm+cfJDZXEdtNNqH1bBxTiHtUzFVuY0rzoIyGAyGy8FsvTeaTjj1mRLDUDwH5TTUReaFb8FoGY3xHXxBXGtvRxo3TEh7NBc0FUcTUazqVRdFOU2HUVAGg8FwEVBV7HOtWNXNiDW5AmhpLJ9yOkGTUezhbhjscBbkzmENkt19nOye70FmkoX+xWFcq25EGjfOWOMpL0s2jUb7kWB5/hxNJ9FkBPGXYDVuuCRBFUtGQU0X7WZYOCw2l7LBMMqMa5jSCYh0Q7gOfDNHBWsmhUZ60NgAxIZQO+NE6/lnLn+hqo7rLhEhe+ynaOcr5w+6fbiueWtuHDeEKuduMcWHkYoV6EA7FJdCNgOZJK5lmycUNpxPloSC8vv99PX1UVFRYZTUAkZV6evrM2HuhkVH4fyRx22xY/cJvvHjA9y5tYl7tq9zFFUqjsYH0XgEmUZBqW2jgx1o3yknCs5TBP4g1ixrMtln9pE9+F+QncRa8gVx3/AepKT6Au/UUZy4/VhVTWhxqZN01uXFWrYZ8V/4cpwLYUkoqMbGRtrb2+np6bncohhmwO/309h44ZOzBsPlYHQNUyqd5ReHOij2u6kpL2bviW7iqQwP3rWJJt8Q4vahI31QVjflWDp8Dj13FALliGtuX8Ea7SV74D/BnrhUROrX47r6jjkrEc2mIT4M3mJnwW8iitRf7SwKDlVB4ybE40emCHefT5aEgvJ4PDQ1NV1uMQwGwxKlszdCQ1WIZ/adodjvpsjnQRWGognCIT879rTxq5sUikphZAC17UnnoTSbQXtPQiA8a+V03i2uZPf+4LxystzgCyCBMFbLLVgVc1/LqckRNJPCqmpBh7uwIz1O1d1ARb6PFayYZoT5ZUkoKIPBYJhPRtcwDY0kKQs6i20TqTQlQSdk/GzPMCQFisNoOg6pGExiyejwObAzsw5WsM8dI/vS98CykGAVOpgrqScW7lsefHWuvOQIWC5cK65FfAG0rA6NdIOnaMogj0vNwpDCYDAYFjCjGWi8bot4MkM8mSaWyHDVsnKGR5LUlXpRnEX7gqCJyIQxNJt26iz5S2Z1TU3Fye77IdgZyKTGlFm3Vr/2VSknZ/wEVs3qvOtOLAurtBaruPRVjXsxMQrKYDAYZmB0DdPG5hq6B0awFW5aV89wLMmul9o4eOIcX3m6i9auKHhy81DkFs62/pzs6Zexzx3LWU+zc1zZR592LLFxSEkNVvPNc5JfE1EnYjDnLtRMCvEVzVpZXi6WRCYJg8FguFSMhpsfOtnN8Y5BNq2qYkUwTWRwiIGkxf2vq6elJIu1fDPZ0y8j3mLnxGzKCUSYJlpPsxlIRNDEMNnnzufdtK5+vZN5PJPAWnXTnAIWNB4Blwt8IRjpR4pL0WgfUrMGa5pgjkuJySRhMBgMF4HRxbZf+o89LK8ppSzox+45SWmoCPHYPHl4gOYbgthn9jnRb6Nl12cov253HCZ74IlJigc2YTVdP/fM5Zmkk1XCG8RqWO9cY6TPcTWKdVmDH2aLUVAGg8FwAYxG9qmdRe0MlsdHqMiioz+BuCtBmLJOU2tXlF2H+ugYSFAf9nNraSerunZO7CgWrg1vnFvmcjuLjvQjRSVIVQsSqkTcXme4qmbs9r1IZXO+bSFjFJTBYDBcAKORfaVeZVR9ROIZ6sJ+xD+1C661K8ojT7dTHvRQH/YxfPooj3ae4b4mL82hFLh9Tgi5ZeFa87q5l12PDWJVr540GayU1CDhZQvGtTcTRkEZDAbDBXDHlhU8/PjLqN1PyOciGkvTH01z97U1wEQradu6Clpqg+w61Ed50EOJT7DPvExJpAv1WuzqDtLS6MN1wz3n563miKYT4PEjUyggsSxcOXffYsBE8RkMBsMF0Fxfyn3X+gl6lM6IEvS7eeDWRlpqg3krKZrI0FDuJ5rI8MjT7bR2RenoHSYwcga79Wcw3AVAyGNzjipcN73vgpUTONF6Vu2aGXP5LRaMBWUwGAxzRO0sdvdxWkJpVt81sbbZqJVUWuwsyB39u+Pn+6nu2s9wGkq956vRRoMrqF9+1aznhVRtSCfB5XbcgekEJKNIqBqZQ4mOhY6xoAwGg2EOaCaFffaAY/1MMT/UMZAgVDT293/QnaGj9QjbqocZSLkZSlnYuBiuWM9gsIntG+ZQmTY2DGI55TVig+D1I3XrsGrXvIo7W3gYC8pgMFzxaDaD3XEIq3bt+bBwQFMxJ/WPiLPINT6I3XUU7Oy0wQv1YT+ReCZvOQEMnT5KrT9BcyjFfVeleCqznnOUU1cR5O25+alZyaqKahZX/TrEW7SkSw0ZBWUwGK54dKQfHezA9hThql0NgB3p5eie3Tx5PE5n3Etd0Ob2lW5WL6+acaHstnUVPPJ0OwChIjdD/f0M9PTw1qYoAGte+0auqrrABNfJESRUlQ9hX6rKCYyLz2AwXOGoKtp/GglVokNnnbRA6STH9u7hsRcGiSZt6n0jRKMJHntxhON9menHiw/TFNnL+2tPUtR9gPYDL1J8bi/3NfXTHEohNaux5qCc1M46ltzofiaBFW644PtdTBgLymAwXNkkhtHECFaoAuwsdu8JEDe7jgxSXlKUd9OV+UFiaXYd6pvSHadDXWSe/RpkU6wCVlUBhVNLYuG6+o5Zi6bJETSdQPwhNNIL/qAT5bfAc+hdLIyCMhgMVzT2QEd+3kl8AexILwCdIxYN5WO/IkNFbjr6E5OOo7FBMi98a0Kl2+MRL7u6g3TFvdSvaOKOiJeWGVLpaSpXsqOoBFfDBvAUoQNnsLtbseo3LGm3XiFGQRkMhisWTSfQSPeYgAcJOGHa9eGRCYEOo5kiRrF7TmCfeA4yaTQ2CKkR54Dbh7X6Fo4Pe/j6niRl9X4aQ8VEsx4eebo9v14qL0cuSSyqTtkOfxCrYSMUl+VrM0nFCkfOV7FOarFhFJTBYLhi0eFuRKwxFsnoItfxgQ6ReGZMpgi7r43s7n+ZWH5dLFxbfxmrYgVP7WyjvOa8khuttDTBTRgfgoqVWP4QeLzgDUxqJYk/dJHufHFgFJTBYFiSjJbF6OyNUFcZ4o4tK2hpPG8pqW2jg2dpHRSe/EXbmJRE4CiR4Xia070xwgEP65eVcPe1NbTUBtFID9nd35lcOW15e778esdAgoZy/5gu492Emk2Dy4NV3rhkMkBcLIyCMhgMS47W9n4efmIf4ZCfhionqevDT+zjwbs2nVdSiWGOnR3isRejlAc9NJQ7a5c+9/gJEGiuCbBxeUnechrNpafxITLPfxMySWccXwDXNW8Dt9fJIF5g5Uy2Hmq8m5BEFKlsMsppEoyCMhgMS44de9oIh/yUBf1obIiSzAiqaXb89CWa33MHIsLRV47x5/95lu4Rm9pSP5VBN73RDHtODeJzu2gs91NafD5d0a5DfTSHhcxz/+zMFwG4PLivv4fj8SC79vfRMdAxJjHsTG5CVRtVGys0hywSVxBmHZTBYFhydPZGKAn40FQMe7gLMilCHpvOs53Yfac41naOR3/yCj0jNrWlPnqHk/zw5R76IkncLotEOsOzxwbpHnJccaEiN509Q2R+8c8QG3AuIhau636J4/HglIlhW2qDPHBrI0G/49YrTCgLOItuS2vHZK8wnMdYUAaDYclRVxlieDhGSfKcU9XW5SEaS1NbXYb2trHjhRHCxS5qSv0k0zZD8QwBr4vBWAaPy8JtCQGfiyOdI1R5UwyeaaUqcQ6CfbkrCK4t78CqWsWunW2TJoYdDYQYfY1H7SyaSeEqq79Uj2XRYSwog8Gw5Ni+eTn95zoZiqVRy81QrlbT9vVVSCBMR1cfoVCAtXUBYskMkXgGv9diOJ4m4HUR8rtQlMGhKAOHf8FAXz/bqiP58V3XvBWr7ipg8sSwoSI3nQOTr5eC81VvrapVV1xk3lwwCspgMCw5Wqp83Lc1SKgkOMG1JpaLhtoKoimoLvVz0+owQb+LwZE0IS9sr4+xfXUxtip2tIeAO30+TVH5clw33IPVuDF/rdFAiEImBEIUoNmMo5xqr8IqXzavz2GxY1x8BoNhyaEjfayuDbFmVemkxwuDFypLfNx6dQU/O9DBzf4zVESTRAZO0uKp5r71J5wy7GLhuun+ScuozxQIMUYuO4vGB7Hq1mGVTjxuGMu8WlAicpeIHBGRVhH51BR9bheRl0XkoIg8NZ/yGAyGpY+qooNd4Js648L44IXlrgF+b8VBlgeTdMS9BDxZ7qs+5ignwGq+eVLlNNlYQb+b+1/XQHO5Cx0ZQKN9TgLaUbde9WqjnGbJvFlQIuICvgDcCbQDL4jI91X1UEGfMuD/AXep6mkRqZ4veQwGwxVCcgTNJrD8ExPetXZF2XWob8yi3FXpVrL7noIi2MbIxPGCFVgtt0x7yfGBEHakF9wBJNyIuH3Yw+cg0o1UrrpiMpFfDObTxXcD0KqqJwBE5JvA24FDBX3uBf5VVU8DqGr3PMpjMBiWIE65jHaktAZxe9GRPpzfx2Np7YryyNPtYxblPvKfe3lf0bM05+IUpKQG13W/RPbYM2j7PrDcuDa9FXHN/qtSY0NIqBqr/up8uiJXsBy115jFuHNkPhVUA3CmYL8duHFcnzWAR0SeBELA36rqo+MHEpFfA34NYPny5fMirMFgWKQkItjdR5FID1bjBnSwEyYpKLjrUN+YcPCQPUy2/xC7PEGaQ/2OcrrpfYjHj/uat6JNN4DHjxTNvrSFZpIgglXTMiGXnlFOc2c+56Amywev4/bdwHXAW4A3An8gImsmnKT6kKpuVdWtVVVmxbXBYDiPPdCOeAOQjmO370czScTlmdCvMBxcU3H09B5CnixdcS8Uh3Hd8B7Ecz7yTkqq56acsmk0PoxVdxXi9r76GzPMqwXVDhTGUDYCHZP06VXVEWBERJ4GrgGOzqNcBoNhEaLZzARXm6ZiuXIZFYgIGh8eo2QKGQ0HL/G7sE/vgWyKSNqiNmThvvG9M5Zxn1IuO+tkI7fcWA0bxpTuMLw6plVQuUCHR1T1/Rcw9gvAahFpAs4C78GZcyrk34HPi4gb8OK4AP/mAq5lMBiWEOMzkW9fV8EqTx9W4zX5+kgA9lAXYnny7rTpLJ7RcHD73DGCsUEiaYuBlJu3vfUWpLjsguRUtdGRAaRiOVb5skktN8OFM62LT1WzQJWIzNleVdUM8BHgR8Bh4NuqelBEPiwiH871OQw8AewDnge+oqoH5notg8GwdBjNRB6NJWmoChGNJXn4By9yrLUNjfTk+2kmBQNnwT95+fXxNNcEeP+Kboojp3Kh5Db3b29hzdrVFySnqkK0H6lcgatqlVFO88BsXHyngJ+JyPfhfAymqv71TCeq6uPA4+Pavjhu/y+Bv5yNsAaDYelTmIkcoLTYi7oSPNlWzOoVJ9FQJYiF3XMckFkFH6idxX5lF6t6n2dVs9MmVc24rn/dhQs60g9l9VgVKy98DMO0zEZBdeReFk6kncFgMMwbnb0RGqoKvmrScUJ+i85IBjIpdLgbTcdhuBsJVkw4P7/WqT9OrTvK7aUdrEofh3Q830fKl+Pa8o5Jq9bOBo1HoLgMq7r5gscwzMyMCkpV/wRARELOrkbnXSqDwXDFMpqJvKzUCVrQkQEiaZeT264ohJ47iqoiwYnBCM5ap9OUpc5RGzlNJJbk0ZSb+5qy59c61a7FtfltF+yS03QSUKzaq0zo+Dwzo4ISkQ3AY0B5br8XuF9VD86zbAaD4Qpk+8Y6Hv6XHWgiTEl5OUNDwwymPbzthgrE5UG9RYjbh8jEKfRdh/ooi5ykJOYswSz1AmT49ukyqoJuznkaqPe3sL0nSUvt3BWU2lk0GcG1fIup4XQJmI2L7yHgd1R1Fzi584AvA6+ZP7EMBsOVSnNYuP/GcnYdjdBxYpDaEg9vv3lZPpWQeIvzbrwDZ4YZHEkTDnhYv6yEA6+cZF3mzPlVmC4PiVAdT3YU89b1jTQWe5wMEk+3jy0cOAmaTaPRPsRTBEUlkE6gqRGs2rVI0eRJaA0Xl9koqMCocgJQ1SdF5MIWDBgMBsMM6PA5WpZVsrp5mZOZwXKPcaWNpiyy1eZkTwyXCO39cY6dHeBkez+twRJuqYpRVVOJ1XgNBw70UV0OZQEnGHl8QcEpSYwg1S3OvNdgB+ItwrViKzLLqEHDq2c2CuqEiPwBjpsP4P3AyfkTyWAwXCk4odq9EChHLBeaSaGxgfxiV3Gfd6ONWk0/fKkLv9dFJmsT9LlJZ216h5P40wM0BZKcGfHylOXltuY1FCVsuoeS3Lmxcsx1Q0VO5vG8HNkMGulBQpX5uSnVLK6SGsRbhIYbwO01oeSXmNmkOvogUAX8K/A9oBL4wHwKZTAYrgx0sIPsmZfRoS6nIT6MwITIuFGrKZrIICIIcLgjSjqTpWsoQVE2SjarVBVlqSnKUlZdxwsnowT9bu7YWInfO/a3+ISCgvEhpKwO4sOOXOkE4i9BvEUAiC9glNNlYDYW1EpV/di8S2IwGK4oNDaI3d2KBCqwe44jgXKnLMUkqYoKE72WBTwk0zYlRR7O9CfIJGK4MwmK3JDICg311bx2SyMd/Qk+tH1FXrnB5AUFNRUHXwCrZg12/EU0m4ZUDKlZe0mfh2Eis7Gg/lpEXhGRPxOR9fMukcFgWPJoIop99iDiDyJuL2K5sbtbOXainS8/08OffPcIX97ZRmuXs6qlMNHr2roAIyMxynWQ4f4+NDlMLGNR5s4S84S5qrlhjIU0WUHB0QAJVXUCH2pWO3JUroTYEApIIHyZno5hlNmsg9omIrXAfwMeEpES4Fuq+r/nXTqDwbCk0Gwae6Ad7W1z3Ga5OSYpKuHo8TM89lwv5eGS8/WactF2+USvRW4q013cYB1lb9ZPvd9N0JPlXMJDkd/D5k3N+DyuCSXXxxcUzBMfRkpq81F5UlKN9p1y3HsmI/llZ1bZzFW1C/g7EdkF/E/gDwGjoAwGw6zRRBS74yBkUkiwfMI6pidPZgiXBvNRdoXRdrc1Znl01wmymX5C2QF8LovVoRT3beqiudTmhLuZp+wNdEUh6Hdz97U100foUVC7qWpVvk0sF1K71sw3LRBms1D3auAe4JeBPuCbwH+fZ7kMBsMSwh7sxD53BPEGpnSddQ6laCgfO/8U9LvYv+8AmQNnGBosoi1tEfb6ubo0wd1rLVbf/N+QcANrLReFM0aqNhobnDJLuao6tZsaN02wlCxTLmPBMBsL6p+AbwBvUNXx9ZwMBoNhWjQ+hN31ClIcnrZ0+qgbb9RyAjjxyhFOdA6zvNpiQ1k8XyJj+7VNrLn5jqnHS8XRbAbSyckzPsSGkPAyrEly+RkWDrOZg7rpUghiMBiWJvZAB+LxT6uc4Hy9JoCgz2K4/QT7T/WxsSxBqdeGYAVlZY1YVilPpwKsnW68dAIprYVIN4xTUJrNgIBVueJV35thfpkxik9EVovId0TkkIicGH1dCuEMBsPiRtNJp4aTb+bsC6PRdsWpftr3PUfx0DFWBdOsDKYgWIW14nqscAMlJQE6BxJTjqOqqAhW+TIUQdUe2yExjFSsNPNMi4DZuvj+CKfS7TacRbomv7zBYJgRjfQ4C2tnWZJiVfoEKzw/hjXO/leOlxNxVxBecW2+ku6ERbbjSSeQojInSrCk2qnb5HdSmaudBbGcdsOCZzbroIpUdQcgqtqmqn8MbJ9fsQwGw2JHbRsdaJ91xVtNxcge/sn5BpeX7TdczVDFOoYTNrYqQ7E0/dE029ZNM3eUijnuPcAqrXWi9UaJR5DyZSaEfJEwGwsqIU486DER+QhwFjA/PwwGw/TEBiCTnHVyVfuVXZDOue6Ky3C/5gHW+gI8sHK0AGGCurB/xhByFcEazTbuL0HcRWh8GERQFKu07tXemeESMRsF9XGgGPgY8Gc4br4H5lEmg8GwSNCRfvAG8pFyqgrxIez+0+jIAOKbuvCBqqI9J9BoDwD2mb35Y651d+bPnXKR7WRjppOIL5SXR0SQ2jWOnIBlFuAuKqZUUCLymKreB7xGVV8AopgksQaDIYem4mTb94NYWBUrwONH+06hqRjiKZo2hFsTUbL7H0e7Wycck5rVWDWr5y5PJoUmIli1Y3PoWYEwmLRFi5LpLKjrRGQF8EEReZRxgRGq2j+vkhkMhgWNPdiBWG7wB9G+U6CAP4AVrJz+vN5TZF/6HqTjEw9ablzr7pzV9dXOQioGdtZ5ubxYDRsQs7ZpyTCdgvoi8ASwCniRsQpKc+0Gg+EKRDMpdPAsUlTqFBOcZfYFe7CD7O5vQzaTb5O6dWC5wE5jLds8ZfaHCcQGoaQG8TkJZwmETej4EmNKBaWqf4eTf+8fVPU3LqFMBoNhgaNDnQjWmEq3M54z0k/2hQLl5AviuuZurKqmuV/fzoLlwqpaZeaUljCzySRhlJPBYABG3WpxtO80FJXM3F8VPXsAu/sY2nMSRkO+PUWcWvVOntxv0zFwhPqwn23rKmYdDEEiCqV1RjktcWaVzdxgMFzZaCqG3dsG0R6nVpLLMyvryd7/+JjoPAAsNydX3M1ju0coD3omlNaYjZJSO4Mrt9bJsHQxCspgMEyJppPYQ51oX5tjrRSHsWaZFcLuODRROflDuDa+macOuCkPyqSlNSZTUJpJovFILnxckEB42hB2w9LAKCiDwTABTY5g951GI92IuJDisrnNN8UGye7/z/y+1KzBteZWCFUhInQMHJlQWiNU5FS8nZREFKludqL2hs9hlV91QfdlWFzMph7Uu4C/wMkeIbmXqurMDmiDwbDo0EwKu30/oEigfNZ59DQ+RPbwTnSwAxIRGE3SWlSK65q7x5S9mKy0xlQ59jSbAZcbq6wesVxoZZOZe7pCmE0uvs8Cb1PVUlUtUdWQUU4Gw9JEVbHPtYKdRfyh2SunoS4yP3sE7TwM8aHzygnBteUdE2oybVtXQX80zVAsPXOOvUQUCS/LW3BGOV05zEZBnVPVw/MuicFguCyobaPJEec10A6RbqS4dNbn2z0nyDz7GCSjYw94inBteCNWuGHCOaOlNYJ+x60X9LvHBEhoNu38VRvFRkpqLvwGDYuW2cxB7RaRbwH/BuTTAqvqv86XUAaD4dLgWExH0OFupwyF2nNKC2R3HSX70r+et5jcPsdimkXG8Kly7GkmhcaHURHEziDh5ZNXxTUseWajoEqAGPCGgjYFjIIyGBY5dt8pGO6+oNLndsdhsi//+3nl5A/hvuE9SKhq2vNau3LZyQcSk69/SsacBbihKifhrMmjd8Uym4W6JkGswbAEsYe60N5TF5S7zm7fT3bvf5xvKA7jvulepGh612BrV5RHnm6nPOihPuwlEhnhkR293L+thdX1zrmKjRUoR7xFiLdozrIZlg7TZTP/n6r6WRH5exyLaQyq+rF5lcxgMMwbmoxidx1BisM45d5mj336ZbL7Hz/fEKhwlFOuau107DrUR3nQQ4k7g6ZsSstKwFvMt545RU1lGR39MeqCwh1vSLJ62SyzShiWLNO9M0cDI3bjJIsd/zIYDIsQtbPYnUcRjx9xzX4ppGYzZI88NVY5hapw3/z+WSkngI6BBEErCW4fVlUTVmktCXeAnYcHiMTT1IeEKEU88qP9tLabgglXOtMli/2P3N9HLp04BoNhvtGBdkhG5uTas/vPkN33Qxg5rzSkpAbXjffO6IZTVSfCL5OmtihDRIsoK6tHLOf38b6T/VSFQ5R6bcgqZeVhJGWxY08bLY2zy5JuWJrMzbY3GAyLFrVt7MFO7N6TMNuSFoB97hjZ574+VjmVL8d10/tmN0eUjIK3GKldy+u338SghhiKpbBVGYwm6B6IseWqZU5ouQBuPyUBH1190RmHNixtTKojg2GJo6qQGObo/n3s2ttBZ8yivnxkVtnDneKCBWHkLi/WVbdjrbh21nNXmk7hql+PFJWwuhQe9IXYsaeNjt4otRVB7tzahM/nRlIBsFyIZTEUTVBbYeagrnRmk+roFlX92UxtU5x7F/C3gAv4iqp+Zop+1wPPAfeo6ndmJbnBYJgWtbPowFl0sJNjZwd47IV+yksCNFS4p80erpEe7DMvo9E+p6yGnXUOFJfhvvHeKQsKThY+3lzpQXxFUDBH1dJYPsZ119rez8NP7IOiEkIBL4PRBAORBG+/Zc1FfyaGxYWoTgjQG9tB5CVVvXamtknOcwFHgTuBduAF4L2qemiSfj8GEsBXZ1JQW7du1d27d08rs8FwJdDa3s+OPW0cPNHNQDRBuEhY39zAHdc1oWqz48nn6TzXR11VCT3RDMVe95jcd0OxNEG/mw9tX5Fvs3tPkn3hX8DOjL2YP4T75vumVU6j4eOhIkcB9kfT3L81yOqN12CV1c3qXrr6HKvqji0rzPzTFYSIvKiqW8e3TxdmfjPwGqBKRH6n4FAJjkU0EzcArap6IjfeN4G3A4fG9fso8F3g+lmMaTAYOG91ZG2bk12DWHaagc5hirMj/PXRdkgnWFWSpaG2nEg8w479vdy5sXKMghqfPdzua+Pojn9n17kSuuJeaotSbKuO0lwTwH39u6ctxT4aPl5YPkPVZueRCGtursgroM7eCHWVoQkKaLxVZTDA9C4+LxDM9SmMIR0GfnkWYzcAZwr224EbCzuISAPwTmA70ygoEfk14NcAli9fPotLGwxLk9Ev+sd/fgyfz002kyXgdeHPxIgX++kYTBBPDIMKW5Y3IOLUXKou9bH3dIQ3lJ0PaijMHq5DXRzd8T0eO1FC2JuhPmQRDW/i6xkf969tYXVw+gW4HQOJMeUz1M4ScqXpjPs43hXl4Sf2EQ75aagKMTyS5OEn9vHgXZuMUjJMy3Rh5k8BT4nIw6radgFjT5YGebw/8XPA/1LV7HRZk1X1IeAhcFx8FyCLwbDoGbWawiE/CAjK4bY+rqorwu9RivweBkfS2Op26uIUfKY2LQ/xk/29DMXSY1xwd19bg9pZMnv/g11dfsJepwSGtepGwr4AVizNk68M5rM8TMVo+YwSv4Wm4ojLQ8T2U9dQw449bYRDfsqCjgIb/WvCyA0zMZsoPp+IPASsLOyvqttnOK8dWFaw3wh0jOuzFfhm7oNUCbxZRDKq+m+zkMtguKIo/KIvC/pJpjOE/Bbt5wZZt6KcRMqmtNhDPJXNn9M9lOBI5wjdQ0nKit3EUhki8TQ1dg9vXpaiuaQW+/hzEOmhK15LfXEWq+nGfLXaaYsIFrBtXQWPPN2OnUxRWtPIsO1lMJPgHdev5ms/PkBD1diFvCUBHx29JozcMD2zUVD/AnwR+AqQnaFvIS8Aq0WkCTgLvAe4t7CDqjaNbovIw8APjHIyGCbS2t7PD589hgBloSKqSnwcOdFHuTfNqagyEE0zHE9T5HXRF0kS8LvZ3dpP+2ASlwguS7i6IUQikeb95YdYlW6Ffsg8+UL+GrVFKaJlawj7z0f1TVVEcDwttUHuv76UXSfTdEShrtLH229ZQ0tjOXWVjltv1HICGB5JmjByw4zMRkFlVPUf5jqwqmZE5CPAj3CCKr6qqgdF5MO541+c65gGw1JjpuCB0T4PP7EPv8eFiJBMpDhy7hyrK92ciRTToGkEJZbMUlXi45Y1YaLJLI/v6aY04KGpqpi1dQGq/FkGju1lZ3eEVc25wfOFBWF7czFfG6nFmsQNOBmqikZ7wHJcii21AdbctHVMmY3W9n66B0bYsfsk1eUBNjdX4/O6TRi5YVbMJsz8j4Fu4HuMrQd1WRJlmTBzw1KhcE6pJOBjeCTJQCTBts3Lae0cyiut3oER/D43qXSWZ/e3UWTHnMW3Lg+rc4X/dh3qI5oYW0L9u7/ooKrEy61XV6LxIeyTL2BnUhwYLGJdWZyuTCm17iEnUq8kg/u1H+R4rJhdh/roHEjgsgCBTJZJy2JofBgC5VjhRjQZRXxBpMD6Kry/RCrD3tZzdA/GuPO6Ju7Zvs7MPxnyzDnMvIAHcn8/WdCmwKqLIZjBcKUyWfBA71CMv/nOC2zbsiIf8fbj3Sd5w/VNVJd4ubFeONLtYiiukM3mF9p+/Wdnx0TRAVSWeOkeTmJH+9BTL4JmOBnxcHLEy8qW1TTW1jPc3cXX+vu5f3Mza0qqaSlx3HWF65riqQw7DvTwjZ+1c8fGSu65qYHmmgCaTeMqb5xUMRVGGm5pqaG2PEjtDUEGowmCxT6jnAyzYjb1oJpm6mMwGOZOZ29kQvDA2Z4I2aw9RmlVh4t5+Xg3d7b4qC71UVNZml9kO2rRjEbRFVpQjeV++nv7GTxygJAnSyRtcSASZNPVKwjX1wMQrq3HKqniqQ43a646L8fouqZkOsvzrYMU+9xUl3j4+ZF+njs2yKoKN+tXVvL6QIqWxvPnTRZp+OzBs9y8voHqcMAERxjmxIzJtESkWEQ+nYvkQ0RWi8hb5180g2Hp0trez/GOAb771Cs8ve8M3QMjAPQMxqgKOxF03QMjPL33NP3Dcfa3dnG8vRd1+xiKpemPptm27nw28m3rKuiPphmKpZ0krCNJpP80H607SMCTpSPuJVDkpbm5iVUrxmZ1CBW56RwYG6nXMZAgVOTmSOcIRR7wS5pMOsOp7hHcZDnTG2Pn0RF+9bM/5M8ffSZfGmN8pKGIUOx388oZ57gJjjDMhdlke/wnIIWTVQKc8PH/PW8SGQxLnFEro7EqhMdtMRiJ8/MD7Rxv78fKpqgvznCuZ4BnD54lmc5QUeJnZdhif2eKA6eHCfrdE3LoteTmooJ+Nx19cYp7DvL+8AG21Y3wq839fPrmLB9+8B1saK4hEj+fxqh7KMFP9vfw0qkhvryzjdYux7oZtcgGown8HgsJhDkbcxMqCeJyuzjak8Zyu6kpL2bviW4efmIfre39dPZGKAn4ALhqeQWxRAZVGIrE8zn27tiyAoNhNsxmDqpZVe8RkfcCqGpcpltVazAYpqXQyigp9vHK6T66+0c4097Jx15XwtPHY+zZ34Pf5UYTbmJJm9c2h/AHghNy5xXSUhukucpH9qXvod2t+XapWY1r89sRt5dt69w88nQ7APFUhqdf6QeFDY3BMfNMNzWHefJwP17LJmEVI1khEktx9cpK2ruHCZYEKfJ5HOUTTRAO+dmxp21MSHl1OMDN6xvYc6wLBYLF50PPDYbZMBsFlRKRInJZIESkmYJoPoPBMDcK556qwwGqyoroPHOG3Yc7+PlpH36/n4FElIDXptRjs2llEdXhICoy7aJZTY44yqn/dL7NWrkVa92d+awSLQVRf08f6qWs2MOych/HumIU+9zUlPrYfzpCIqXc1lKMzwO7jsWpCnu4ankFmaydV1QAiVSakqA/P7f0vtevdzKT4yzG9XpcrFlWYdIaGS6I2SioPwKeAJaJyNeBW4AH51Mog2EpU2hlaCbNuTOn+OmhHkqDRTRUFBGJZ/B53KxrDNFc4MYbjqXHLJrVSA9224vgDSDFZWSPPAWJ4fxxq/lmrLW3M97h0VIbpKU2mM+f97Mj/RT73BR5XSjK4Eia8qCH410R/uBX7uLeviw79rRx6GQ3xzsGaaovxe2yiCfTxBIZGitD/PiFEyTSWWorgvkw+dF6T8ZqMlwos4ni+7GIvATchJNf77dVtXfeJTMYlih3bFnhWBnpOMFULy+fGgLLw5aVpVi55K6blofYfzpCZYlv0kWzGh8m8+zXIB2f9BrW2ttxtbxm0mOjjM4zDcXSlOWi/0bTJQXdWTqHXUhRKS2N5BVMa3s/39h5iJ0vnqQqHGB1Qxn7TvQAcPvm5URjSXa9fNpYTIaLwmyi+N6Jk03ih6r6AyAjIu+Yd8kMhiVKS2M5D2xrojjeQ2fEJmlb3HZ1OdWl562jpuoAq6qLnaCH/sSYwAi1s06V28mUk9uHa+u7Z1ROcD7yz+u2iKezxFNZRpJZ1tQ4i4brGhsmlf0P7n8tX/7kW7jjuiZaO4YoC/nZtmUFNeVByoL+/HyUwfBqmZWLT1W/N7qjqoMi8kfAv82bVAbDEkbjQzRl2vi1NzQjbh9f3tlGNDG2QGAknmH9spIJARGqin3oJ+jgaN5lwVp2DTrSD24vrnWvRwKzs1xG56O+9exZdhzopbrUx41NxfhEGXRX8M4b1059bq5+0+h8mlXgRjRrnQwXi9koqMmsrNmcZzAYxqGZFPbZg4gvgLidcOzRTODAtDnwNBUne+A/0c5X8m3W1dtwrbpp7nKoQnyI5uoQv//OtdxzcwM7D3TT1Rch1NjMO26Y3byRSQRrmE9mo2h2i8hfA1/AieT7KPDivEplMCwBNJPC7juNFJdhhZyoN7v3FNg24vHl+xVG1nX0J6gL+7n72pox65zs3lNk9/5gTBCE1K7FahpTA3R2ctlZdKQf8YcgEYXiUidEPZRGKjZjVcy+KGh+Pg3G5BM0iWANF4PZKKiPAn8AfCu3/1/Ap+dNIoNhCaDxIeyOw5BNYw+2Q+3ViMeLDnYgwYoJ/Ucj6yaMk0liH96JfXrPmHZr+Rasda+fEKE3pTzZNGRSYGfRTBKrqhkJVZE9+TzkEkar2lgl1XO6z5bGch68axM79rSZqD3DRWdaBSUiLuDfVfX1l0geg2HRo4kI2dMvO248fxCyGezOQ4jLi/hDs1Iqdu9J7DP70HNHIZs+f8Djx7XxzVh1V019cqEs2TTEh8Hth0AY8fiw/KH8PJWU1sJIPyBIqArxzFz7aTyj81EGw8VmWgWVK8UeE5FSVR26VEIZDIsVzWawOw9zfEB58tg5OgYS1If93H5VGS0V7jGuvanIntqNffC/JrRLzWpcG+5yXHOzkSWTRJMxrJo1SEk1Yrkm9LHKGsgOdgKCq35i1J7BcDmZjYsvAewXkR8DI6ONqvqxeZPKYFjgaHwYTScmuMTsvjZazw7w6AsRyoMeGsqdtUaPPtOZCxOfYdzYIPbhnWMbA+W41rwOqVs3a5eeI2MEq3Ej1iQuxVHEH0SC5ZBOgL9k1mMbDJeC2SioH+ZeBoMBsKN92GcPALkveG8xABobQPtPs+tEivKgJ1/6YvTvrkN9k84zjaKqZA/+F9i5kPNQFe5r7oaSmjkpJnDSHkmwfFYh51bVKjSTnvM1DIb5ZjaZJB7J5eJbrqpHLoFMBsOCxY70Yp89gBSXQjqJ3d2K1bARUiNO+HhRCZ2D/ROKB4aK3BPy6KmdRc/uxz71IpqKIYEKtO9U/rhr45udOaI5oqpoOoGrYf2slI74gsjMnkeD4ZIzo4ISkbuBvwK8QJOIbAb+VFXfNs+yGQwLCk1GsTsOIsWliMsDLg8a6cXuPwMDZ8HtQ9y+SYsHRuIZXBZ8eWcbHf1xaq0BbpO9NHv7z4+fiOS3reVbsMIXOCcUG0DCDYjPrEUyLG5m4+L7Y+AG4EkAVX1ZREyVXcMVhRP88AriLXKU0yjFpei5o44V4i2itStK91Ain5nhmuUh/F43x7tGQKDIylIbOUKkv4fHUl7ua/LSHEqNvZi3GGvt7XOXUdWJyAtWYFWYmkuGxc9sFFRGVYfGuQp0nuQxGBYkdt8pSMUmzOmIy5OfI2rtivLI0+2UBz28fmMl+05H+PH+Xu7YWElduQ9/cojQ2QOQTVHqBciwq6+c1VuvcuaBhjrR+DBWwwbEWzStPBobAl8AcTkfYVVFo31IeSNW5SrEmk0tUoNhYTMbBXVARO4FXCKyGvgY8PP5FctgWDhobBDtP4MEKyc9PvrjbdehvnxwRGmxh5oSL4O9PQQZorN7gMrESaceQI6Smga6fCtxtax3xpnlfJNmcuXY4kMwGqEXG0LK6nFVt1zYTRoMC5DZZpL4fZwihf8M/AhT8t2wRFE7iw6dQ0prEMuF2lnsc8ccF940AQetXVF+uOccoIQDXtZUWFT27iGYTtAR91JblCKStij12s5cVeM1RF2l1PkvIK1lIorUrkUHO9BUDFweEMGqXHnB920wLESm/HSIiB/4MNAC7AduVtXMVP0NhsWIpmL5MHEAHTjrzDUlo1g1q9Hhbse1N81aolHXnt9jIQKJeIJnXzrDTeUZvJZFbVGKbdVRHjtZDsUllDZvIpJxTZoQdkZ5sxlwuZFgJeILkm17EXQEq3494vZe8HMwGBYi0/18ewRIAz8F3gRcDXz8EshkMMw7moph97ahw11I5SonQWoqht17EimtgcEObMuCoS4oKp12rFHX3uYVJfz8UAfFsXMUSZqXB/ysLs1w96YAzeVV3L+6gacGa+gcTFEXdk9ICDsrEhGkssmZe3IFkXAjJKNIqOpVPA2DYWEynYJap6obAUTkH4HnL41IBsP8ovFhJ1eey40EytHek9ioM4/j9jmuvWA52t+OuDz5QISp6BhIUF+UItj3Cje6hjgqPgZtFyA88LbrWXO1U1dpbe41rWxqAyAyMchBsxkUxmSvsCqbQLNmka1hSTLdJy+foVJVM+YDYFgq2ANnEbfXSeQKECxHe08hCBJyXHkiFhKaPChiPLWuQYZeOUCpJ011EVQXZRhKewituoY1V6+ltSvqlNLI5eXbtq5iastpZNBRUr4A4i1C7SwkoqhmEJfXyUJe4MpzovVMxJ5haTLdO/saERnOvSLAptFtERme5jyDYcGi6SQa6QZfIN8mYjn1moJzy8itqmRPPMdtiZ8ykBSGUha2CkNFyxiqvY47broqPz8VTWRoKPcTTWR45Ol2WrsmVpxVVVTAtWwTiGBHep16TeEGXCuuw1p1E1a4/lU/A4NhsTClBaWqE1MfGwyLHI32OtbRJB6B0bbZWDx232nswzvQoU6aQ3BfUz+7Bqo5V7aRuppy3pE758s722afly+X7kgC5VhFpZCIgD80aRZyg+FKwJRuNyxpNBXHHj7npA0SF9rfDv6pAxMKF9uOZiJ/5On2XCbyIJqKkT34Y7Tj4JjzWpZXs/ZdvzwmIhCc+anZ5OUDIJ3Ayq1jEssFxWUXdtMGwxLBKCjDksFJ9TMARSWIy42mk9jt+5zovEg3El6GZhJY/sCUYxQutoXzFs/Og72sstvIHvgvSMfPn2C5sJpuwFr92rEpkHJMlZevLux35FXbCcpQG0WgyJS8MBhGMQrKsHSID5I98zLi8SOVTehAu6MAQlVOWHnnYcQ3tXKCyS2eoDvD2UP7yCYOj2mXuqtxXbUNmcbS2baugkeebgccyykSz5xf/xQbcJSULwB21qloO4mSMxiuVIyCMiwJ1Laxu48jRSVguZxS6S4vUuRUnxVvMXiKZgzHrg/7GY6lKbUSaKQHjQ0w1D9Ajbug7Lo/hGvjmzhh17LruT46Bs5NOVfVUhvkgVsbnTmt/gR1YT93X1tDc4UHbMFVuwb77EGnPEb16ov+XAyGxYxRUIYlgUZ7ITlyPuPDJJkfZlJOmk1za/EpHn3pLLZECXlsImmLgZSbt9Y7UXfWims5UXYD3/5ZDzsO7M1nLB+NzhudqyqkpTZIc00g784Dp66Ua9kmpDiMtXwLdu8p494zGMZhFJRh0aPZDNp9/FV9waudJfvS91jV18p9jV52dQfzOfTe2tBPc10I14a7OJGu4NGn2znWFaWmzIcgPNc6xM2ryygPeth1qI/m6iInkasqo9lhFSfwQbMZEBx3Xi4zuvgCuBrWX4QnYTAsLYyCMiw61M6CWKCKjvSjfW1gZy54/kZVye79AdrdCkBzKEVzWRSpbsGqWumkEwpUICLsyoWNpzI2ZcWevFV2pHOEW9aWO9F5IwNQsQKrqATUBrcXvAFH5tggduQcVsXKi/U4DIYli1FQhkWB2jbEB7H72yE+6FgnIo5p4itGAuG5j5mKoZ2vYHccQvtP59utphux1t46qcIbDaIoLfaQSNsUeV34vRaDI2ki8Qy1xTaEqrEqV07uUgyW45rjgmCD4UrFKCjDgkazGSdYoe80ZJLg9U8oGjgVky24BSeU/GzHOWoG97GtcnBMRVtr+Rasq7dPOV81Gja+ti7Ac8cGHBlV8bot+odjvPU1tVg1LSY3nsFwEZjXJF4icpeIHBGRVhH51CTH3yci+3Kvn4vINfMpj2FxoakY9qkX0O5j4PEiwfIJC2GnYrIUQ597/ASf+88TRAYHqBvcx0jS5rGT5RyPOLntrOVbsDa8cVrlsm1dBf3RND6PixtayrBVOTeUZOOyEPddX8aaLVtN2QuD4SIxbxaUiLiALwB3Au3ACyLyfVU9VNDtJHCbqg6IyJuAh4Ab50smw8JD40NoJoU1rlyEqmKfOwbIrC2mQgoX3GomRSjRTU/HIADXJM6BZpyy6y6L78Supbqsns6jQn3PmUnDxQutMb9HiKUyZLNwx4Yqtq2roLkkC6FqxB+6wCdhMBjGM58uvhuAVlU9ASAi3wTeDuQVlKoWlo5/DmicR3kMCwx7sBO76wigaHgZVlVTPgxbh7pgZGDWGcXH09Eboc41SPbcOYj2AkoiWgIqUJormW65SVZvZOfRBHfXemgod09IbdTaFeVbz55lx4HefEi53+umP5o+n/4om4Fk1KkpZTAYLhrzqaAagDMF++1Mbx39CvCfkx0QkV8Dfg1g+XLzJbDYUVXsnhNo/xknuEEsGDyLHR9AQjWItxi7+9i0GRqmHNvOYh/6MdVdbQynxSmxnsNnnd9GLGTFdew/kaW61DdpMleAR55u52jHENUBQdCJIeU1AYgPItWrEY/vwh+KwWCYwHwqqMkc+TppR5FtOArqtZMdV9WHcNx/bN26ddIxDAsPTScnfGmrKnbvSeg/7ZQtH53vCVY4pTD6T6O2jXh8iMs961pKrV1Rdh7s4ewr+6nNdtBcnOSpnhCQcRbcusJUVZeCuBguq6ekoooIxXQPdXHnRsdK6x5KcKRzhMGRNKjSPZwgHHCTSqUpKw0i6QR4XRzpHOE1q0vp7BlCR/ocpVpaO9+P02C44phPBdUOLCvYbwQ6xncSkU3AV4A3qWrfPMpjuIRofBj7zMsQqMAqXwYeJ7+dDnWifaeRYMWEYATx+KBAoc2UWXxMv6faKI2coD7bQSRt8VRPiNubvJywmjmnYeqqSvlEQRRf50CCurCbOzZW4ve66R5K8NyxAYp9bvweCwV27O/l9VcFKC0LkbSKKPL48cUGGRxOEYn6qa2twtV0w6wDNwwGw9yYTwX1ArBaRJqAs8B7gHsLO4jIcuBfgftU9eg8ymK4xNh9p8Dtg0SE7Ok9jhsPELWdaLxJSpqPZ6rM4oW1lDSbYcdP91DadZJSK0Z3ws3RqI9uLae9v5Y/eNfa/Dlf/9nZCeHm3UMpTnYPYqsSKnI+DrGUzc2ry9hzcoC9ZxNs2biW5w53gt+NesrwlcKgt5J33LrJKCeDYR6ZNwWVKxP/EeBHgAv4qqoeFJEP545/EfhDoAL4f7lf0xlV3TpfMhkuDRobRKP9TpVaQJi8/tJM7rvRRbFq22hsAOLDFMcjnE24sa+KotF+7LaX6DgRor4oRXfczS/6AgSCJdRU1tI9nOJzj58AgeaaQN4KK2zbuLyEoN/Fv73QRUO5n4byIjYtD1FV4mNjnZedJ7P4fB5uXFfP3tZzdA/FufO6Ju7Zvo6WRrPg1mCYT0R1cU3pbN26VXfv3n25xTBMgao6rr1selrrotB9V1iGotB99+WdbQz391HSux/ScbrjbvYO+knYFm+qH2ZbdZTmUIqvHC9nJG2xfzhEqqiSonAliZSNz2MRT2UBeMOm6vy1/2tf94xtmk4wlHYR91VQHQ7Q1ReltiLIHVtWGMVkMFxkROTFyYwTk0nCcNFQ24ZoDxobyltP4xm1mn74Uhd+r4vNK0ooLfZMcN+pKreGzvLo80dQT5pE1s1PuwOAcGt1lJG0xWMny7mtKkJPuphdA9X0p72sLgsgKZuRZJZNy0P8onWQ8bE5ybSdb1Pb2d7YWMSOg4MMxdIE/S6GowmGPFU8aCwlg+GyMa+ZJAxXDna0D7ttN3bnK0hx6aR9CrM7iAgCPHtskO4hp/x5qMhN50ACe6Cd7M8fYVX3k9y3so+Ax+aF/iClpSFu21BDzfIVlFZUki0q4++7NhBcezN3bm3C73Vz4EyUVDbLzavLqC714/NY+DyuMXL4PBZel6CJqJM+SW2KPG62rykm4M7S2RMhVFrGg2+5zigng+EyYiyoeaK1vZ8de9ro7I1QVxlaNK6hucrd2t7PjmcPcrbtJPUVpWzfVENLcPKs4oVBD2UBD4lYjOJUhFdOxqm6ZiXDkTjV0aNkf36A4xEvu7rL6Yp7qQsXsXJ1CxtXVWGJOOHggyPs6x4mm7VJZZXasiLeubWWJw/34bYsKkt8DMXSVAa9IDAUSxMqcjM8kqTSb4PLw7CnglC4jEgszWAqwQPvuJFVVjeaGMbVtBXxTT53ZjAYLg3GgpoHWtv7efiJfURjSRqqQkRjSR5+Yh+t7f2XW7Rpma3cmhxBsxmn/+N7iHR30FBVykjGWdja2hWddPyOgQShIjeqymr/ICPnTmNHexno7GDg0LP0vbKb231HOR7x8tjJckbSLhqWNxKrv44T/VlOnhvJh4Mn0zZZ28bvtfJWWHWpn1uvKieZytLRnyDod/Pbd9bz27dVEJAkHV0DBFxZPv6e1/CJD7yJYLiczr4YwWIfD961idWrGrGWX4Nr+WajnAyGBYCxoOaBHXvaCIf8lAWdtT+jf3fsaVuQVtSo1fT4z4/h87nZ0lKDJTKp3JqKYbe9BJbFjhdTlOkwpQEP4vZQmns3FYaBF1JX6mboXAcl8Q6qoj3cWOEEPYxkhL1dacIem13dQXoTbsKV5YRXXIV4iykDNi0Psf90hGKfiyKf47JziUV1iY+Az1k8W13qp8jr5s3X1vKh7SvQRAQUrLpNrLkWyKahqBTJrclavWziPJlYLriADBYGg+HiYyyoeaCzN0JJYGwGhZKAj66+yS2Ly0mh1YSAoDx78CzdAyPAWLnVzjq589xe8BbTceYsIU86/4UP5+eRxmO37+fW/h/Qf/IwQ/192Ao+l01lsZs1JWmuKYuzvizOiBaxK7qSdNXVY6IAm6oDrKouJpG2SaSy+DwW29eXY4mgKIMjKYZiafqjabatq0Djw+DyYa3YggTKnVdJzRhZDQbDwsZYUPNAXWWI4ZFk3gIBGB5JUltxYW6j+ZzPKrT2yoJ+kukMxT43r7T1UB0OjJHb7j8NiUg+u3h9TZhoIkNhSEQknqEuPFYJZE+9iH3wRzQH4L6mdL6cen1DHQ3FKylyKSXxswCEK5uoPjTI3tMR3lBWNGbc9ctKWL+sxLlmLuqvIuTj5bZhQAj63dx9bQ3NVT5IxbEaN5jSFwbDIsYoqIvIqCI5eKKbE52DbFpVRVN9mOGRJAORBG+/Zc0FjfnwE/sIh/w0VDmK7+En9vHgXZsuipLq7B2mvsSFjgywpkJ47pUBiqwsA2mbwd4gA2k316yq4kvf3MnZ9nbqq8NsX++lpTbItnUVPPJ0O8CYtUxvXR8ke+RJNB6BbArtOpK/XnN1Eau3XI1Vvw4JVfEn3z1CSYkfq3Rtvs+m5SF+sr83H9gwOu7d19YA5K8Z9Ak+j4vVtcEx66fsSC9Ww3qjnAyGRY5RUBeJQkWysbmaQJGHfcd7GEmkWddUzdtvWXNBCuVizGepKqQTiLdobLttU+tNMNTZR1nAS7XP4qZVQV5uj4GdpTjVy8amRnY98yJhT4aG6jDRpD0mH94DtzY62SD6E9SWunnzyl5WHv4P7Gx6ghxSVo/rhveMcbONVqgdtYgAirxu7thQSdDvpqM/QV3Yz93X1uQV0P03VbDzQDedvVlqQy7e+rrl51MfJSJOnr/ghZXpMBgMCwejoC4S4xVJS0M5laXFBIt9/PrdW2Y8fyo3XmdvhIaqsUXwSgI+OnpnP5+lQ53YXUeR8mVOzSLLDZkkdu9Jti2zebTbg5VxMjr4/S7W1Lp44NZGmmuK+fITrxD2QllZGQClxU6Aws4DPaxy97Iq2U9T/RCtqS52Horyz3EXtUUhtlU78u3qDjqh4hUB7rjmTlaPmwOaygobnxAWnDkwRgZoWVbJms3Xgi+I3XMSHWhHU25IJ0AEq7rZlFw3GJYARkFdJF6NIpnOjTfVfJZlCV/6jz0zzktpcgT7XKuzeHawA3uoy2m3M4hYrG5exgOBaN4KGm+tdMZcNJSPVSpBYpzdf5Bs8gRAPiw87FXqi1JE0hZ/d7wRK1RNU6WXhqCfqKuUR3/ezQO3esconvFW2Pjr5+8jk0QTEazqFiTckFdAVnUz6vGhIwNIuBEpLptgKRoMhsWJUVAXiVcTGDGdG++OLSt4+Il9gKPwhkeSHD87gCVQ7HNPUGjNlR7srmNIWR0SrMA+dzRXW8kDgbBT/VUEyzqfXaGlNjhpWDiMdcFpKo72n2Ho7ElqPNl8n13dQcLejFMc0FtMKtTIi8ddJCI2I94AawMBqgM+RNKThqCPXn80DVJh1vGW2qBTJyoVw7Vs84QihiKClC+D8mUYDIalhVFQc2QqV9xkimQgkmBzc/WMls501ldLYzkP3rWJHXva6Oh1EpY2VIbw+9wTFdrzR2jamAWXBz13FO1xO0lbC/LiiWtu//Jt6yp45MeHyba1kYgMsm/QT0+ihO21EY5Hi2lZu4pzHUXULwtiBSvoicHzxwZIZmJ4LCGZtnn22CBraovojqQ52xfPjzuhptPT7YSLLeqLs0SGhnnkyQj3v66B5lLFtfwapGjyFEoGg2FpYtZBzYHpMi2MKpJgsaNYgsU+tm1ezq6XT8+YmWHU+ipkOBKnNuysA2ppLOfX797CHz34On797i2kszYlAR+aSmD3tWH3txPUKJ2nToLbh/gCSLACfMUQnGUgRSKCffplskefJrv/P8ke3oHdvp+mvp/zPs9TxCKD7OhylOjrayMUBYJ8PXMrp2q309C8hhFfNeLxc6RzhGKfmyKvG7fLosjrwrZtdh7qYziWobGiiGgiwyM7jnHsTG/++rsO9RH22ZR4srhKayirayRc5GLn3k5cyzYa5WQwXIEYC2oOzBRRN/oa5Uv/sWdWEXgTrK9ogv6uDt56Sw2ajE5Iu1NXXsRwTzcheygXSp0l0tdHbbh4zPyLuCbPiVeIpuLYrT/HbnsR7MykfZpDUOXP8JaGIUrLw0hZMz2Uc+xMhP/x9UNcv6qUzsEkzTUBBkfS+D0WAa8LESWeyjIYS5NMZ7EVrqoPUuJVtKSYXYf7aanwgMfP2XP9NFSXYZXW5MPDS5YF6OgeRIrDM96HwWBYehgFNQfmGggx2/4T3HgBeMt1AVoqvWRPvYiEGyAXKk5yhNsqhnjsSDdaGqTE5yESzzCQcvG2m2pnfS+aHME++Tz2qRchm5pw3EnW6kTg1RalOJyoZuOmq3D5g3QPJfjFsQGKfC4EpdjnBk0SS2VAFQVev9FxKx7pHKE3kqYi6M1nGLcTUUqqG+kYSEFxGJIRGlasJJp1U1awdikSS1FXbZSTwXClYhTUHJhrIMRc+o9aX5rNYJ/aDe5qxO0Djx8iPU7JdLHA42dNU5AHikonRL6BU+SvYyBBXamXbc0+WhrKQCy08xXs9n1opNsphZRNM75OkpTWIlWraI34+NqZCGW+KPVBJeqv5eRZN6FBaK4l78YDKAt4KS320FwbIOh389l7mnn0mQ68bouSYg8+j4veSIpNy0JUl/rRdBLxFTOctqirLsXVsA5V5Q7vgGNFioyZw7uQxc0Gg2FpYBTUHJgqEGKqL9HZ9lfbzikMOPrKUXY910Zn3F0QyRaaMPb4yLvWriiPPHmCslQPtel+htuGefgF4b6mfppDEy2kMQQrcK25Daldi4jw9M42KpafXzwbBja5ouw/HaGyxJd348VSNpuWO7KFitx09MdpDge4//Ymdu45TUd3nLqqEj5+VxNPHu5ncCRF0Eox4ikf8xxEZNJgkAtd3GwwGJYGRkHNgdl8iaqdRQc7QG2aK8t54A3r2Ln7OGfPnqOuOszbcymKVBWSUexIDwx1onaW1u44j/38HOHSEA3lXiLxzJisDZOhmRQ63M2OHQcoPXeWEo+j6Eo9gFrs6g7SHJq8zIeU1mE135xXTKN0DCQmrH1qqg4wksgS9LvzbrxRlx04ufJqA4qUNbJ2bQtrNmx0EstmU4g/xIqqYnbuOeMs2q0p4e23TYxmHD+HZzAYrmyMgpoj032JaiLifCknR8ByYfeeognhVzaAaBFKHHH3YA8m0YEzaCruhH37gliWi6dOtlEeLs1bLuPLoIOzYNVuP4AOnkWHOiHaB0BHey31RWPTC4X8HjpiOddgoBxr2TVY9evA5QWRKXPVTZZ+aDRZ64e2r8hnf/B5XNiquewPSd5yfSlWeSMA4g9iNazHbtuDpuI0l2RofuMarPr1cw51NxgMVybmm2KOOJbPCJqIIsHy/Je8xofJnt6DeIudEG9gQrIdVYj2ocPnwBfAyvUbZTLLxXGdOeUrNDZI5vlvwshEi6g2l8GhNBRAKpuQYAXRjJtGvxv3tuVzSv0zVfqhu6+tQe0szaEM928NseuVYTq6ItSWF/HWrSWsvvqqMXn2xOPHWrYRu20P+ENOglijnAwGwywx3xazRFMxNNKLDnai2QRi22hxGKthPa1n+tix4xk6h7PUVwYnLEIdRUSgqGTy8ZMj1KY7GDqZcCwXl5vuuLC3yyaRUb70r/3cpnto9g9OPDlYyfaNdXytrRKrqoqSYg/DBUpltspJ1QmamCr9UHNY0NggUrGSNY0lrNnqgmwGe6QfUlGktG7iPfuCWCu2OOuzZhH2bjAYDKPI6JfSYmHr1q26e/fueRtfbRtSMXSkD430OmuD1EYzKURc4A/kv2g1NkTrADz6zFnCPigpDUyb7HSU0ZQ+Hf1x6gI2t5V2sGr4JY4PWrmcdhkSWYufdgcA4dbqKD6XzUDK7QQ9lGSw1tyKlC9HSqrzVtzouJ0DjlLZts6x0HYd6qNjIDEmfdAYOQYS1JW6ub3Jw+qaQN4CPP9MnCStFIexappNOXSDwXBREZEXVXXrhHajoECTUezek5AYgUwSFRxl5C3OhXeLUwp8Eh564hWi8TRl4fOZDoZiaYJ+Nx/avuL8NVTR3lMcO3SIR/emCXuzhOwokXjyvOIJpfLrj57oKMFn2WwOJ6jyOwtoh1IWAZ+LX3/vnU5W8hkYTR9UHnQylQ/H0gwMRrj/thWI25s/FvQqw5EYg1LG/deHaKn0Ir6AI3c6gSZHsCqbnGSslkk+YjAYLi5TKagr0sWndhayOcso0oPde9KZO/EWIf7gxLmjaegcERrKz7vtNJ0kkOjh7OkhMsXPOI2+AER60Ug3O4+XE05blIoNQKkXIONE29W7WL12PavtDOd2jlBXnMGyA4yuVypx+ziUqOYf9yodA0cmWETj2XWoj/Kg53zQhSsF4VK+9bM2jvXZdEdS1ITcrK32U7NiFVYSnjyjNJfGwVvkBHuo4lq+BZnCNWkwGAzzxRWpoOzOV9Bor2MZgVOiYQoLaTwaG0R7TqBqY9WsyUe8lbjTaK420XBSqfHYaN/EYIauuJf6ooJ1SeImVF5Ol1WD+3Wvyc8XNbS3jSltDnD8XJQTnRFWJDI0lPunDEPXTApcnjFBF5pKgLeYhJaw8+UefCSpKfGRTNn8ogNurshSWVZMR28UqViJnjuGBCucwIZxNZwMBoPhUnBFKihnbU7QydQwDZocwe44hA6chdQIGhuE+FD+uH3wx9zqWcWjh22yRAh5bCJpi4GUm7c2TLL2yOWmvqGOqDsXSu72If4Sooks9X73mGCGySLp9p2OsGlZaPow9JFBWntT7DzUy57Wfg65LTavCFBdWoyU1rLvxTaqq8L4NUkyHqc4XI6kbV45088mj4vaiiBWuBG13E5mCRN1ZzAYLhNX/LePxgaxu45CJomUVIM/hA6cRXtPoj0nQO3pzmZV+jj3NTjzRh25vHV3r7VoWXMjUlIDLrfjKgOkuoU7+rM88nQ7ljs3L1QQbVfIZJF0zdXFNNUE0FTCKTjoLRoXhj5E65Dw2AEhHGzghmvKeHp/J0+dstnYXMyZk228eKSLa5qrqa6r5OjZPiRt4/O46e6PMlAV4u23rEFcbiS3nslgMBguF1dUkITGBtH+02RPPI/GBpw1Sf2n5zaI5UYqVkA2PfZcsZCKlVirbnDWIU0T2j1ZtN1U80iFfHlnG5GROKU+QUJVaKSboViWoN/iQ6+thqISvrInQzSRzef/6x4Y4ad7T3Oya5CNzTVk0lmyqqjCmsYwPUNxugdGqCwr5o8eeJ3J5GAwGC45JkgCsNteIvPCt9D4ECizikiTsgakYT0SKHci2wLl58PMIz1o/xkoKkXKl02ZmQFAs2lHiVmuMXn0VG3IRcrh8U+YC9NMygl19/i5/eowj+7sQ+qXUeIPMZyyGBju5pqmOr6yP0vn0BAvHe/mpqvr8+dXhwNUlhWjwBuvX0X3wAjPHjyLCHQPxbmmuZqBSIIHcymYDAaDYaFwRSmowsqyY9qrViHBCnS4G40PO2uLKlZgVa1CAlN/aUuoCglVzXhdTSectEaWhZ3NIAiIE5t3vDvOrpMZOgdT1BVnuX1NiJaagppO3mIn3H1kgOZghvvftJknT2XzuQA3r9nMrpdPEw75aagJcOh0H0++fJptW1ZQHXZCxXsGY1TltqvDAW5e38Dhtl7O9gxzy8ZlJimrwWBYkFxZCqqkBqvuajQdR/wlSKjSWezqn5gt/NWg2QygiMuDZpJoKoZr+RbwBiAZgXQS3B5au0Z47PAxwiWlNDY5hQq/djjGA40ttNSXgstzPpWSnYVUnDXeYtZuPG/5jS+KuGV1Lbv2tLGn9Rx3bm1ieCSJy7JorDzvQqwOB/B6XNyyaTm/fveWi3rvBoPBcLG4olZdSkkNnjs/jmvDXVjNNzmJSy9AOanaaDrJ+Pk7VUVjQ05QRCaFPdKPJkdwNW5C/CHEspCiUsdCKw6z82AP4ZIiyoJ+LBHKQkWESwPsPHDOKdte4DIUy+VEHo5zS3b2RigJnI9GrA4HuG3zcpKpbL70/CfefT2WZTEYTWCrMhhNMBBJcMeWFRgMBsNC5YqyoC4EtbMQH0ZVEQFFESzwBdGRAbAscHkgk0LVdkpYVK5E3N7c/FF2TBn2QuZaoXcyJiuK6Pe6efPNLWOsoxW1ZabWksFgWFQYBTUFqraTtTybwapahVVUihNZIeApRizLsaKGzznzS8VlWP7AmDx10wVNwNwr9E7GbIsimlpLBoNhsXFlKiiXB432oTJCzixySmGIE7yArU40XaAcV+XKKS0g8fiQWeTEm4q5VuidDFOJ1mAwLFWuqHVQo2g2nc/FB5pLdaeQi67D8oDbO6caShdKa3s/O/a00dXnKJc7tkysNGswGAxLGbMOqgBxeZx5owWAcb0ZDAbD5MyrghKRu4C/BVzAV1T1M+OOS+74m4EY8KCqvjSfMs2GUaumszdCXWUoH+02m7bplM1k4xrlZDAYDJMzby4+EXEBR4E7gXbgBeC9qnqooM+bgY/iKKgbgb9V1RunG3e+Cxa2tvfz8BP7CIf8+Xmh42cHsASa6sPTtk2XkWGycU0GB4PBYJjaxTef66BuAFpV9YSqpoBvAm8f1+ftwKPq8BxQJiIT64ZfQnbsacsvfLVEKAv66R2K0T0Un7EtHPKzY0/brMedrr/BYDBc6cyngmoAzhTst+fa5toHEfk1EdktIrt7enouuqCFjF/4CpBMZ0mlMjO2lQR8dPVNvoZpsnGn628wGAxXOvOpoCYLgRvvT5xNH1T1IVXdqqpbq6pmzn33ahhdm1SIz+PC63XP2DbdGqbJxp3rmieDwWC4kphPBdUOLCvYbwQ6LqDPJeWOLSsYiCTGpAWqLC2murRoxrbp0gdNNq5JN2QwGAxTM59BEm6cIIk7gLM4QRL3qurBgj5vAT7C+SCJv1PVG6Ybd76DJGDytUnArNpmE8Vn1jwZDAbDeS75OihVzYjIR4Af4YSZf1VVD4rIh3PHvwg8jqOcWnHCzD8wX/LMhanWJs22ba7jGgwGg2Ei87oOSlUfx1FChW1fLNhW4LfmUwaDwWAwLE6uqHIbBoPBYFg8GAVlMBgMhgWJUVAGg8FgWJAYBWUwGAyGBcmiK7chIj3Aq80PVAn0XgRxLieL/R4Wu/yw+O9hscsPi/8eFrv8cHHuYYWqTsjCsOgU1MVARHZPFnO/mFjs97DY5YfFfw+LXX5Y/Pew2OWH+b0H4+IzGAwGw4LEKCiDwWAwLEiuVAX10OUW4CKw2O9hscsPi/8eFrv8sPjvYbHLD/N4D1fkHJTBYDAYFj5XqgVlMBgMhgWOUVAGg8FgWJAsGQUlIl8VkW4ROVDQtllEnhORl3MVeW/ItXtE5BER2S8ih0XkdwvOuS7X3ioifycikxVVvFTyXyMiz+bk+Q8RKSk49rs5GY+IyBsvt/xzvQcRuVNEXsy1vygi2y/3Pcz1f5A7vlxEoiLyPy63/BdyDyKyKXfsYO64/3LewxzfQwvxc7xMRHbl5DkoIr+day8XkR+LyLHc33DBOQvqszzXe5jXz7KqLokXcCtwLXCgoO2/gDfltt8MPJnbvhf4Zm67GDgFrMztPw/cjFPt9z9Hz79M8r8A3Jbb/iDwZ7ntdcBewAc0AccB1+WU/wLuYQtQn9veAJwtOGfB/w8Kjn8X+Bfgf1xu+S/gf+AG9gHX5PYrLvf7aI7yL8TPcR1wbW47hFMTbx3wWeBTufZPAX+R215wn+ULuId5+ywvGQtKVZ8G+sc3A6O/Fks5X61XgYA4RRWLgBQwLCJ1QImqPqvO030UeMd8yw5Tyr8WeDq3/WPgl3Lbb8f5YCZV9SROPa0bLqf8MLd7UNU9qjr6/zgI+EXEt4j+B4jIO4ATOPKPti2a/wHwBmCfqu7NndunqtlF9D9YiJ/jTlV9KbcdAQ4DDTif2Udy3R4pkGfBfZbneg/z+VleMgpqCj4O/KWInAH+Chh1AXwHGAE6gdPAX6lqP84/ob3g/PZc2+XiAPC23Pa7gWW57QbgTEG/UTkXmvww9T0U8kvAHlVNsvDuYVL5RSQA/C/gT8b1X2jyw9T/gzWAisiPROQlEfmfufaFdg9Tyb+gP8cishLHuvgFUKOqneAoAKA6121Bf5ZneQ+FXNTP8lJXUL8BfEJVlwGfAP4x134DkAXqcczq/y4iq3DM0PFczjj8DwK/JSIv4pjaqVz7VHIuNPlh6nsAQETWA38B/Ppo0yRjLMT/wZ8Af6Oq0XH9F5r8MPU9uIHXAu/L/X2niNzBwruHqeRfsJ9jEQniuH8/rqrD03WdpG1BfJbncA+j/S/6Z3leK+ouAB4Afju3/S/AV3Lb9wJPqGoa6BaRnwFbgZ8CjQXnN3LeLXjJUdVXcNwwiMga4C25Q+2MtURG5WxnAckP094DItIIfA+4X1WP55oX1D1MI/+NwC+LyGeBMsAWkQTOB3rByA8zvo+eUtXe3LHHceZ/vsYCuodp5F+Qn2MR8eC8D76uqv+aaz4nInWq2plzfXXn2hfkZ3mO9zBvn+WlbkF1ALfltrcDx3Lbp4Ht4hAAbgJeyZmtERG5KRdtcj/w75da6FFEpDr31wI+DXwxd+j7wHtyft4mYDXw/EKTH6a+BxEpA34I/K6q/my0/0K7h6nkV9XXqepKVV0JfA74P6r6+YUmP0z7PvoRsElEinPzOLcBhxbaPUwj/4L7HOeu94/AYVX964JD38f5wUzu778XtC+oz/Jc72FeP8sXI+pjIbyAb+D4otM4mvtXcNwWL+JEyfwCuC7XN4hjUR0EDgGfLBhnK47P+zjweXLZNi6T/L+NE0FzFPhMoSzA7+dkPEJBZMzlkn+u94DzRTMCvFzwql5M/4OC8/6YsVF8i+J/kOv//tzn4ADw2ct9D3N8Dy3Ez/FrcdxY+wre12/GiZDcgfMjeQdQXnDOgvosz/Ue5vOzbFIdGQwGg2FBstRdfAaDwWBYpBgFZTAYDIYFiVFQBoPBYFiQGAVlMBgMhgWJUVAGg8FgWJAYBWUwzEBunc0zIvKmgrb/JiJPXIJrv1ucrNK7xrXfLiI/mOKcx3NrUyY79nsF2yulIGu4wbDQMArKYJgBddZifBj4axHx5xaF/jnwWxcynoi45tD9V4DfVNVtsz1BVd+sqoPjrim5ha6/N/lZBsPCwygog2EWqOoB4D9wEsT+EU46oN8XkRdEZI+IvB3yVslPc8lXXxKR1+Tabxenxs4/A/vHjy8i7xWnbs4BEfmLXNsf4iya/KKI/OUkYpWIyPdE5JCIfDGngBCRUyJSmZPlsIj8P+AlnOwAReLUR/t6bgyXiHxZnLo//yUiRRfxsRkMrwqzUNdgmCU5y+klnGSlPwAOqurXcu6053GyPitgq2pCRFYD31DVrSJyO046mA3qlFUoHLceeA64DhjAqWP2d6r6byLyJE6Wit3jzrkdeAKnTk9bbvtLqvodETmFs4I/iFMO5DWq+lzuvKiqBnPbK3HKO2xV1ZdF5NvA91X1axflgRkMrxJjQRkMs0RVR4BvAY8BdwKfEpGXgScBP7Ac8ABfFpH9OGl41hUM8fx45ZTjepximj2qmgG+jlO4byaeV9UTqprFSRH02kn6tI0qpyk4qaov57ZfBFbO4roGwyVhqWczNxguNnbuJcAvqeqRwoMi8sfAOeAanB+AiYLDI1OMeaGlvMe7PyZzh0x1zVGSBdtZnMJ/BsOCwFhQBsOF8SPgo7kszYjIllx7KdCpqjZwHzCbgIhfALfl5o1cwHuBp2Zx3g0i0pSbe7oHeGYW56TFKaVgMCx4jIIyGC6MP8Nx5+3LhWr/Wa79/wEPiMhzOBVrZ7JgUKcswe8Cu3Ay77+kqrMpS/AsTnbvA8BJnHo8M/FQTuavz9jTYLjMmCAJg8FgMCxIjAVlMBgMhgWJUVAGg8FgWJAYBWUwGAyGBYlRUAaDwWBYkBgFZTAYDIYFiVFQBoPBYFiQGAVlMBgMhgXJ/w/7y/iRWqrV/gAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.fill_between(xs+offset, low, high, \n",
" color='C1', alpha=0.2)\n",
"\n",
"plt.plot(xs+offset, median, label='Logistic model', \n",
" color='C1')\n",
"\n",
"plot_data(data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The model fits the long-term trend of the data, with most data points falling within the 90% credible interval despite the apparent deviations from the trend.\n",
"\n",
"The model predicts that people born between 2000 and 2020 will be more likely to support legalizing marijuana (when they are old enough to be respondents in the General Social Survey)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise:** The GSS dataset we used to explore generational changes contains more than 100 other variables related to the attitudes and beliefs of the respondents. You can [read about these variables here](https://gssdataexplorer.norc.org/projects/52787). Choose one and run an analysis similar to what we did with `grass`.\n",
"\n",
"How well does the logistic model describe the long-term trend for the variable you chose?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"celltoolbar": "Tags",
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.6"
}
},
"nbformat": 4,
"nbformat_minor": 1
}