{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Bootstrap Tutorial\n", "\n", "This notebook contains a tutorial on how to use the bootstrap functionality provided by estimagic. We start with the simplest possible example of calculating standard errors and confidence intervals for an OLS estimator without as well as with clustering. Then we progress to more advanced examples.\n", "\n", "In the example here, we will work with the \"exercise\" example dataset taken from the seaborn library.\n", "\n", "The working example will be a linear regression to investigate the effects of exercise time on pulse." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import seaborn as sns\n", "import statsmodels.api as sm\n", "\n", "import estimagic as em" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Prepare the dataset" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
iddietpulsetimekindconstant
01low fat851rest1
11low fat8515rest1
21low fat8830rest1
32low fat901rest1
42low fat9215rest1
\n", "
" ], "text/plain": [ " id diet pulse time kind constant\n", "0 1 low fat 85 1 rest 1\n", "1 1 low fat 85 15 rest 1\n", "2 1 low fat 88 30 rest 1\n", "3 2 low fat 90 1 rest 1\n", "4 2 low fat 92 15 rest 1" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df = sns.load_dataset(\"exercise\", index_col=0)\n", "replacements = {\"1 min\": 1, \"15 min\": 15, \"30 min\": 30}\n", "df = df.replace({\"time\": replacements})\n", "df[\"constant\"] = 1\n", "\n", "df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Doing a very simple bootstrap" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first thing we need is a function that calculates the bootstrap outcome, given an empirical or re-sampled dataset. The bootstrap outcome is the quantity for which you want to calculate standard errors and confidence intervals. In most applications those are just parameter estimates.\n", "\n", "In our case, we want to regress \"pulse\" on \"time\" and a constant. Our outcome function looks as follows:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def ols_fit(data):\n", " y = data[\"pulse\"]\n", " x = data[[\"constant\", \"time\"]]\n", " params = sm.OLS(y, x).fit().params\n", "\n", " return params" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In general, the user-specified outcome function may return any pytree (e.g. numpy.ndarray, pandas.DataFrame, dict etc.). In the example here, it returns a pandas.Series." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we are ready to calculate confidence intervals and standard errors." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(constant 90.795969\n", " time 0.169532\n", " dtype: float64,\n", " constant 96.648003\n", " time 0.628462\n", " dtype: float64)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results_without_cluster = em.bootstrap(data=df, outcome=ols_fit)\n", "results_without_cluster.ci()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "constant 1.486144\n", "time 0.120032\n", "dtype: float64" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results_without_cluster.se()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The above function call represents the minimum that a user has to specify, making full use of the default options, such as drawing a 1_000 bootstrap draws, using the \"percentile\" bootstrap confidence interval, not making use of parallelization, etc.\n", "\n", "If, for example, we wanted to take 10_000 draws, while parallelizing on two cores, and using a \"bc\" type confidence interval, we would simply call the following:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(constant 91.357344\n", " time 0.197786\n", " dtype: float64,\n", " constant 96.346417\n", " time 0.609604\n", " dtype: float64)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results_without_cluster2 = em.bootstrap(\n", " data=df, outcome=ols_fit, n_draws=10_000, n_cores=2\n", ")\n", "\n", "results_without_cluster2.ci(ci_method=\"bc\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Doing a clustered bootstrap" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the cluster robust variant of the bootstrap, the original dataset is divided into clusters according to the values of some user-specified variable, and then clusters are drawn uniformly with replacement in order to create the different bootstrap samples. \n", "\n", "In order to use the cluster robust boostrap, we simply specify which variable to cluster by. In the example we are working with, it seems sensible to cluster on individuals, i.e. on the column \"id\" of our dataset." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "constant 1.164278\n", "time 0.100821\n", "dtype: float64" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results_with_cluster = em.bootstrap(data=df, outcome=ols_fit, cluster_by=\"id\")\n", "\n", "results_with_cluster.se()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the estimated standard errors are indeed of smaller magnitude when we use the cluster robust bootstrap. \n", "\n", "Finally, we can compare our bootstrap results to a regression on the full sample using statsmodels' OLS function.\n", "We see that the cluster robust bootstrap yields standard error estimates very close to the ones of the cluster robust regression, while the regular bootstrap seems to overestimate the standard errors of both coefficients.\n", "\n", "**Note**: We would not expect the asymptotic statsmodels standard errors to be exactly the same as the bootstrapped standard errors.\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: pulse R-squared: 0.096\n", "Model: OLS Adj. R-squared: 0.086\n", "Method: Least Squares F-statistic: 13.75\n", "Date: Thu, 30 Jun 2022 Prob (F-statistic): 0.000879\n", "Time: 21:04:30 Log-Likelihood: -365.51\n", "No. Observations: 90 AIC: 735.0\n", "Df Residuals: 88 BIC: 740.0\n", "Df Model: 1 \n", "Covariance Type: cluster \n", "==============================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "constant 93.7611 1.205 77.837 0.000 91.400 96.122\n", "time 0.3873 0.104 3.708 0.000 0.183 0.592\n", "==============================================================================\n", "Omnibus: 20.828 Durbin-Watson: 0.827\n", "Prob(Omnibus): 0.000 Jarque-Bera (JB): 26.313\n", "Skew: 1.173 Prob(JB): 1.93e-06\n", "Kurtosis: 4.231 Cond. No. 31.7\n", "==============================================================================\n", "\n", "Notes:\n", "[1] Standard Errors are robust to cluster correlation (cluster)\n" ] } ], "source": [ "y = df[\"pulse\"]\n", "x = df[[\"constant\", \"time\"]]\n", "\n", "\n", "cluster_robust_ols = sm.OLS(y, x).fit(cov_type=\"cluster\", cov_kwds={\"groups\": df[\"id\"]})\n", "\n", "print(cluster_robust_ols.summary())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Splitting up the process" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In many situations, the above procedure is enough. However, sometimes it may be important to split the bootstrapping process up into smaller steps. Examples for such situations are:\n", "\n", "1. You want to look at the bootstrap estimates\n", "2. You want to do a bootstrap with a low number of draws first and add more draws later without duplicated calculations\n", "3. You have more bootstrap outcomes than just the parameters\n", "\n", "### 1. Accessing bootstrap outcomes\n", "\n", "The bootstrap outcomes are stored in the results object you get back when calling the bootstrap function. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[constant 93.732040\n", " time 0.580057\n", " dtype: float64,\n", " constant 92.909468\n", " time 0.309198\n", " dtype: float64,\n", " constant 94.257886\n", " time 0.428624\n", " dtype: float64,\n", " constant 93.872576\n", " time 0.410508\n", " dtype: float64,\n", " constant 92.076689\n", " time 0.542170\n", " dtype: float64]" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "result = em.bootstrap(data=df, outcome=ols_fit, seed=1234)\n", "my_outcomes = result.outcomes\n", "\n", "my_outcomes[:5]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To further compare the cluster bootstrap to the uniform bootstrap, let's plot the sampling distribution of the parameters on time. We can again see that the standard error is smaller when we cluster on the subject id. " ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "result_clustered = em.bootstrap(data=df, outcome=ols_fit, seed=1234, cluster_by=\"id\")\n", "my_outcomes_clustered = result_clustered.outcomes" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABA3ElEQVR4nO3dd3yV9dn48c91Vk72DiEkAcIGQZYIuHCjrVvrqtbRWqptrT62T2ufp7WtT8evw1ZttVpba4e2jjpxoIgLFJlhQwgjO4Hsfcb398d9QEbGCeSMJNf79Tovzrnv732fKyHJdb5bjDEopZQaumyRDkAppVRkaSJQSqkhThOBUkoNcZoIlFJqiNNEoJRSQ5wj0gH0VUZGhhk1alSkw1BKqQFl9erV+4wxmV2dG3CJYNSoUaxatSrSYSil1IAiInu6O6dNQ0opNcRpIlBKqSFOE4FSSg1xmgiUUmqI00SglFJDnCYCpZQa4jQRKKXUEKeJQCmlhjhNBEopNcQNuJnFSnXngSXbgy5717njQxiJUgOL1giUUmqI00SglFJDXMgTgYjYRWStiLzaxTkRkQdFpEhECkVkZqjjUUopdbhw1AjuBLZ0c+4CYFzgcRvwSBjiUUopdYiQJgIRyQU+B/ypmyKXAE8Zy8dAiogMD2VMSimlDhfqGsFvge8A/m7OjwBKDnldGjh2GBG5TURWiciqmpqafg9SKaWGspAlAhH5PFBtjFndU7EujpmjDhjzmDFmtjFmdmZmlxvsKKWUOkahrBGcAlwsIruBZ4CzROTvR5QpBfIOeZ0LlIcwJqWUUkcIWSIwxnzPGJNrjBkFXAMsNcZ88YhiLwM3BkYPzQUajDEVoYpJKaXU0cI+s1hEFgEYYx4FFgMXAkVAK3BzuONR6iiF9x3bddOO8TqlIiwsicAYswxYFnj+6CHHDXBHOGJQSinVNZ1ZrJRSQ5wmAqWUGuJ09VEVtfqymqhS6thpjUAppYY4TQRKKTXEadOQGpQa2jwUVTezv7mDDq+fOJed7GQ3o9LjiY/RH3ulDqW/EWpQaenw8lHRPrZUNgGQEOMgxmGjvL6NjeWN2AQmDEvkylm55KXFRThapaKDJgI1aFQ0tPFqYQUdHj+zRqYybUQySbFOAIwx7GvuZFN5A5vKGznnN+/xzbPHseiMMdhtXS15pdTQoYlADQqlda28uK6chBgHl80YQUZCzGHnRYTMxBgWTMhi1shUKhvb+eWb23hvWw2/v34mmYkx3dxZqcFPE4Ea8GqaOnhlfQXJbidXzBpBnKvnH+tEt5P//fxkXlxXxvde2MClv/+IP990EhOyE7u9ZkXx/l7j+LjKGu5617nj+/YFKBVhOmpIDWgdHh+vFpbjcti4bEbvSeAAEeGyGbk8+9X5eHx+rnhkOcu2VYc4WqWikyYCNWAZY3h7SzXNHV4unJpNgrvvFdypucm89PVTGJkex1eeWsVbmypDEKlS0U0TgRqwiqqbKappZl5BOsOTY4/5PsOTY3n6trlMyUnmjn+u4a3SlP4LUqkBQBOBGpDaPT6Wba8hKzGGmfmpx32/JLeTp26dwwkjkrn9gzEs0WSghhBNBGpAWrm7ltZOH2dPzMLWT8M/k9xOnrplDlPSWvn6R2NYsy++X+6rVLTTRKAGnPrWTtaX1DMlJ4msJHe/3jvR7eTPZ+wgO7aTL783jj1NOqxUDX6h3LzeLSIrRWS9iGwSkR91UWaBiDSIyLrA4wehikcNHst37sduE+YVpIfk/uluL39ZsB0D3LRsPLXtOspaDW6hrBF0AGcZY04EpgMLA/sSH+kDY8z0wOPHIYxHDQL7mjvYUd3M9LyUkK4ZVJDUwZ9OL6KsxcW3VhTgNyF7K6UiLpSb1xtjTHPgpTPw0F8ndVw+2VWLy27rlw7i3szKbOYHs/byfkUyL5cXhPz9lIqUkPYRiIhdRNYB1cASY8wnXRSbF2g+el1EpnRzn9tEZJWIrKqpqQllyCqK1bV0UlTdzIl5ybid9rC85/Vja7h45H7+XTKOzQ1pYXlPpcItpInAGOMzxkwHcoE5InLCEUXWACMDzUcPAS92c5/HjDGzjTGzMzMzQxmyimJrS+qx24QTc1PC9p4i8NM5uxnubuGhomk0epxhe2+lwiUso4aMMfXAMmDhEccbDzQfGWMWA04RyQhHTGpgaev0saWikYnZiWHfTyDB6eeb49bR5HXx1O5JYX1vpcIhlKOGMkUkJfA8FjgH2HpEmWwRkcDzOYF4el/dSw05G8oa8PoNM/JSIvL+I+ObuXTETj7an8OaOq2VqsEllDWC4cC7IlIIfIrVR/CqiCwSkUWBMlcCG0VkPfAgcI0xRjuU1WG8fj/rS+sZmRZHekLkxvVfmlNMbmwTT+yaTKs3PH0USoVDyOrYxphCYEYXxx895PnDwMOhikENDturmmnt9DEjPyWicThshq8WbOQHm+by9N4J3FqwOaLxKNVfdGaximrGGNbtrSc93kV+FGwtOTaxgYXZe3inOo9dLd3vX6DUQKKJQEW1qqYOapo7mJabTKA7KeKuyC0i3uHhH3smog2ZajDQRKCi2qayBhw26XH3sHCLd3i5IreITY3prK3XjmM18GkiUFGr0+tnW1UT44clEuOIrs7Zc7JKGO5u4R97JuAz0VFTUepYaSJQUWtbVRMen+GEEUmhfzNjYP8qKPwh7P037PkXlL8GjdvB7zuquMNmuC5/G+XtCSytyg19fEqFkC6rqKLWxrIG0hNcZPfzUtNHqf4Q1twNtZ+C2MCZCmKHlmKoXQXOVFJlHnWOw9cbmpVazcTEWl4oG8MZWWWhjVGpENIagYpKG8saqG7q4IScEHYS+9rh09vh7dOgvRJmPwxX7IfxX4dxX4OJ34b8q8HmYGLnYkZ2fgjGf/ByEbgqt4h6j5ul1VorUAOX1ghUVHp65V7sNmFiqDqJ2yrhvYusT/wT74ZpPwbHETuS2RyQNBESx1Ox7SVyvOtwmWZ2uM6zag7ApKRaJibW8nJZAdeN8uOw62crNfDoT62KOq2dXl5aV874rISQrDKa7NsDb82Hxi1w+ksw89dHJ4FDiY3drtPZ45xPhq+I0Z4PODBuVMQaTlrncbO5orHfY1UqHLRGoKLOq+sraO7wMmVEcr/fO9m3l6vqbwCnF85aChlzgr623DkTp2kjx7uWFlsG1Q5r1fQpSbWMT6jj090OJueEoWNbqX6miUBFnac/3cvYrARykvu3kzjOX8MV9TfhMB38Le6v7FubAmw/qtzclu7XPdzjnE+cfx+jO9+n2TaMVlsGInB57k5+vjWVLRVN/RqzUuGgTUMqqmypaGTt3nqunZPfr53EdtPJRQ3fIM6/nxeSn2CfY+Kx3UiEHTHn4iWGsR1LEGMNLZ2WvI+sxBjW7K3Dr/taqgFGE4GKKs+s3IvLbuPyGSP676bGcGbzj8nxruWtxJ9R7Txyf6S+8Uocxa4FxJv9DPeuA6y+gpn5qdS3eli6tbofglYqfDQRqKjR1unjP2vLuGBqNqnxrn6777T2fzK1/Vk+iVvEdveF/XLPOkcB++0F5HlWEuNvAGBsVgIJMQ4e/6C4X95DqXDRRKCixuINFTS2e7nmpPx+u+eIzpUsaP4pxa4FLI+7s9/uC7DbeToGId+zAgC7TZiRl8Inu2rZUNrQr++lVChpIlBR4+mVexmdEc/cgv7ZJD7RV8bnG79Jgz2P1xN/dXDsf3/ptCVQ7pxJhq+IRF8FAFNGJGmtQA04odyq0i0iK0VkvYhsEpEfdVFGRORBESkSkUIRmRmqeFR021HVxKo9dVw7J69fOokdppWLG+/AhpeXk/5Apy00E9PKHTPokHjyPR+BMcQ47FxzUh6vbaigvL4tJO+pVH8LZY2gAzjLGHMiMB1YKCJzjyhzATAu8LgNeCSE8ago9vTKEpx24YqZ/bBUgzGc1/R9Mr1beT3x10etEdSf/OKkzDGbJH/lwSaim08dDcCTy3eH7H2V6k8hSwTG0hx46Qw8jhxXdwnwVKDsx0CKiAwPVUwqOrV7fDy/ppTzpmT3y57EJ7U9zoSOxXwYfze7Y87ohwh7Vu2YTIckMLflYTCGESmxXDh1OE9/spemdk/I31+p4xXSPgIRsYvIOqAaa/P6T44oMgIoOeR1aeDYkfe5TURWiciqmpqakMWrIuOV9eU0tHn44skjj/teozuWcUrLb9ga8zlWxX6lH6LrnRE7ZY6ZjPCuJs/zMQBfOW00TR1e/vVpSS9XKxV5IU0ExhifMWY6kAvMEZEjB3B31Rh81GwcY8xjxpjZxpjZmZm6I9Rg8/eP9zA2K+G4O4lTvTu5oOm/qHZMYkni/1mD+8Ok2jGZZlsWc1utWsG03BTmjE7jLx/txuvz934DpSIoLEtMGGPqRWQZsBDYeMipUiDvkNe5QHk4YlLRYUNpA+tLG7jvosnH1Ukc42/k4sbb8eHilaTf45VY5rY81I+R9syIg0/jbuPM5vuhehkMO5NbTx3NV/+2mrc2V3HhVG3xVNErlKOGMkUkJfA8FjgH2HpEsZeBGwOjh+YCDcaYilDFpKLP3z/eQ6zTzuWzjr2TWIyPhU33kOwr5dWkB2my5/RjhMHb4P4CzbZM2PgTAM6ZNIyR6XH8SYeSqigXyqah4cC7IlIIfIrVR/CqiCwSkUWBMouBYqAIeBy4PYTxqCjT0OrhpfVlXDojhyS385jvM7/ldxR0vseyhO9T5jqpHyPsG5/EsCb2Jqh6F2rXYrcJN88fxZq99azZWxexuJTqTShHDRUaY2YYY6YZY04wxvw4cPxRY8yjgefGGHOHMWaMMWaqMWZVqOJR0ef5NaW0e/xcfxydxOPbFzOn7Y8Uuq+m0H1tP0Z3bDa6vwCOBNj6AABXzc4j0e3giQ93RTgypbqnM4tVRPj8hieX72bWyFROOMZ9B7I8mziv6XuUOWbybsL/hLVzuDsdtiQYcyvseRpay4iPcXDdnHze2FhJaV1rpMNTqku6H4GKiCWbK9lb28r3Lji25aDjfVVc0riINlsaryY/hF/6b5G64/VEzUXcbB5i1Zs/5qOE/8Lj8+M3hm8+vZbTxh096u2uc8dHIEqlPqM1AhURj3+wi7y0WM6bkt3na52mhUsaF+EyzbyU/AittowQRHjsGu15FLnOZVr7MzhNC4luJ+OyEthY1kinV4eSquijiUCF3Zq9dazeU8ctp4zGbutbc44YHwsbv02mdyuvJf722DeYCbHVcbfgNo1Mbv8PADPyUun0+XVfYxWVtGlIhd0TH+4i0e3gqtl5vReGw+YDjOz8iBzvWnY5TyPbW0i2tzBUYR6XSud0yh0zmNn6JIXua8lOdjM82c3avXVMy03GFgX9GUodoDUCFVYlta28vqGC607OJyGmb59DsrybyPGupdIxlUrHtBBF2H/WxN1Eir+Egs6lgLWDWWO7l+KalghHptThNBGosPrDsiLsNuGm+aP6dF2qt5iCzmXU2fLZ5TwtKkYI9abIdS4Ntlxmtf4ZgILMeJLcDtbqnAIVZTQRqLDZva+Ff68q5fqTRzI8OTbo65J8JYzvfINmWxbbYxb2+wYzoWLEztrYGxnhXUO2Zz02EabnpVDe0E5lQ3ukw1PqoIHxG6UGhQfe3o7TLtx+5pigr8n2rGdix2LaJYWtMRdF1TDRYGx0X0G7JDKz7S8ATMlJJsZhY9We2ghHptRnNBGosNhW2cTL68u5af5oshLdQV2T7t3OZQ1fwSNxbHZfgleCuy6aeGwJbHBfzbiON0nyleBy2Jiel8LOmhb2NXdEOjylAE0EKkx+/dY2ElwOFp0R5G5hTTu5vOEWvOJic8wleCQ+tAGG0LrYGzDYmNH2FADT81Jw2oVVu7WvQEUHTQQq5NaX1PPW5iq+cnoBKXFBNO00boe3z8BuPLyQ/Bdr2YYBrNmezfaYCzmh/Tli/I24nXamjUhhe1UT9a2dkQ5PKU0EKvR+9dY20uJd3BLYy7dHDVvhnQXg7+TZlKfY7xgX8vjCYXXcLbhMK1Pb/wXAjPwUbDZh1R6tFajI00SgQurj4v18sGMfXztjTO/zBuo3wTtngPHDOcvY75gQniDDoMYxib3OeUxv+xs200l8jIMpOUlsqWikpFYXo1ORpYlAhYwxhl+9uY1hSTHcMK+XpabrCq2agNjh7GWQPDkMEYbX6tibSfRXMaFjMQCzR6YiIvz27R0RjkwNdbrEhAqZZdtrWLWnjvsvPQG30w7AA0u2H1Uu07OZKxpuxitunkt+kvqVNuDoctGu160xjaFV0jit+Zcke/eCCOtyr+E/a0tZdEYB44YlhidQpY4Qyq0q80TkXRHZIiKbROTOLsosEJEGEVkXePwgVPGo8DLG8Ou3tpGXFssXelhTKMuzkSsbbsIjcTyb8nfqHaPCF2S4iVDhnE682UeSvxSA2SPTiHM5+PVbAy/xqcEjlDUCL/Bfxpg1IpIIrBaRJcaYzUeU+8AY8/kQxqFCqKtP+AA7qpvYWNbIeZOH8ft3i7osk+1Zz2UNt9IhyTyX8lca7ce+b/FAUWMfTx4ryPGspdGeR6zLzldOK+CBt7ezvqSeE/NSIh2iGoJClggCm9BXBJ43icgWYARwZCJQg4zfGD7eWUtanIsJ2Yc3dxxoPknwVTKp42W8EkuR6+yDyzUPdkYcVDpPJN/zMfH+agBuPW00f12xm1+8sZV/fPlkZACso6QGl7B0FovIKGAG8EkXp+eJyHoReV1EpnRz/W0iskpEVtXU1IQyVNUPtlU2UdvaydwxaV0ut3xoEtgUcxmdtqHVNl7pmIqXGHI91hbdCTEOvnnWWJbv3M+SzVURjk4NRUElAhF5XkQ+J9L31b5EJAF4HviWMebIXTnWACONMScCDwEvdnUPY8xjxpjZxpjZmZlHb/WnoofPb/i4eD9ZiTGMzUw46nyCr+qIJHB0mcHOJzFUOKaR5ismw7sVgOvnjmRcVgL3v7aFdo8vwhGqoSbYP+yPANcBO0Tk5yIS1LZQIuLESgL/MMa8cOR5Y0yjMaY58Hwx4BSR6Np3UPXJpvIGGtu9zCtIP6qJY5inkEkdL+EVN5tiLh2SSeCACueJeHFycusjADjtNn5w0WT21rbyxIe7IhydGmqC6iMwxrwNvC0iycC1wBIRKQEeB/5ujPEceY1YfwWeALYYY37T1X1FJBuoMsYYEZmDlZj2H9uXosKq8D4A5rZ89t/V6bfxVPHpTEhs42r3G8gh+6/E+6uZ3P5iIAkMveagI/nETaVjGuM63oSGLZA8idPGZXJuoHP9ylm5DEsaeIvsqYEp6KYeEUkHbgK+DKwFfgfMBJZ0c8kpwA3AWYcMD71QRBaJyKJAmSuBjSKyHngQuMYYY47tS1GR9lZlPnUeN1fnbT9s3xi3v45J7S/jlZhATWBoJ4EDKpzT8RALG+8/eOx/PjcJr8/w08VbIhiZGmqCqhGIyAvAROBvwEWBEUEA/xKRVV1dY4z5EOhx+IMx5mHg4eDDVdGq3Wfn5fICpiXvY1LSZ+vnOP3NTOp4GYAtMRfTOcAXkOtPXollfex1nLTnCZjyXUiZysj0eBYtGMOD7+zg8pm5nDFe+8RU6AVbI/iTMWayMeZnB5KAiMQAGGNmhyw6NWC8VZVPk9fFVXmfLZdgNx1M6ngFp2lnS8xFtNtSIxhhdPo07ivgTIZ13z147PYFYyjIjOf7/9lAa6c3gtGpoSLYRHB/F8dW9GcgauBq99l5rXwU05JrGJvQYB00fsZ3vE6sqWNbzAW02IdFNsgo1WFLgSn3QvliqFoGgNtp52eXTaW0ro3f6TpEKgx6TAQiki0is4BYEZkhIjMDjwVAXDgCVNFvaXUujd4YrsjdefDYKM+HpPhLKXYtoMGeH8HoBoDxX4e4XFj7HQh0kZ1ckM61c/L404e72FjWEOEA1WDXW43gfOBXQC7wG+DXgcfdwL2hDU0NBJ1+G6+Uj2ZK0n7GJ9YDkOXdxHBvIeWOE6lxDL5VRPudIxam/QRqP4W9zx48/N2Fk0iLd/Gd5wrx+PwRDFANdj0mAmPMX40xZwI3GWPOPORxcVfzAtTQs7Q6l3qPm8tzrfWEEn0VjO58j3pbHnucp0Q4ugFk1A2QMhXW3wu+dgCS45zcf+kJbK5o5NFlO3u5gVLHrsdRQyLyRWPM34FRInL3kee7mx+ghgaPX3ilfDQTE2uZnFSHw7QxvvMNOiSR7THnQ98nog9JBxbuyzd3c0XzzSx/5b/5JP6Og+fHD0vggbe3U9XYTnpCDHedOz5SoapBqrff1AM7hicAiV081BD2+t5UajtjuShnFxjD2I63cZh2tscsxCc6Gaqv9rrmsy3mAua0Pkqyb+/B42eMzyTGYWfJlir8fp1mo/pfjzUCY8wfA//+KDzhqIHCGHhi2zCGu1uYnlLDcO86Uv17KHaeTqtNx74fq/fiv8fozvdY0Hw/LyX9EUSIczlYMCGT1zdWsqZE9zhW/S/YRef+n4gkiYhTRN4RkX0i8sVQB6ei15p98azfn8DC7N0k+qvI96xgv72AKsfUSIc2oLXYh7E87k4KOt9jTOfbB4+Py0pgTGY8HxfXsrOmOYIRqsEo2Ebc8wIrh34eKAXGA98OWVQq6v15WzZJTi8LMvYwrvMtPBLPTtdZoGvpH7d1sV+kxj6BM5vvx+W3/uiLCGdOyMJpE77zXCE+bSJS/SjYROAM/Hsh8LQxpjZE8agBoKy+jTdKUrl2bA3j/B8Raxoocp2t/QL9xIiDtxN/TLy/mtNbfnbweHyMgzPGZ7J6Tx1/Xb47cgGqQSfYHcpeEZGtQBtwu4hkAu2hC0tFs6cCf4S+nL+SjIr1VDqmDoltJkOhpw3vKxzTmdr+HA7TTr19FAAnJ4InZy7/782tnD0pi5Hp8d1er1Swgl2G+rsi8gug0RjjE5EW4JLQhqaiUYfXx7OrS/lcbiWZ+5+nXRLZ45wX6bAGpRLnyaT6djOmcynr3dfhFTcicEX2Wj6sOpvrHv+EK2aO6HFrSx1qqoLRl4Hek4CrReRGrOWjzwtNSCqaLdlcRW1LJ98e/nforGWn62z84op0WIOSETtFMefgMO2M6nz/4PE0Vwenjc+grL6NDbr8hOoHwS5D/TdgDLAOOLCPngGeCk1YKpIOTHDqyn/WlnFS8m5y29+n0jFFm4RCrMWWRZlzNnmeldR7R7LPMQGAKcOT2FHVzIdF+xiVHk9SrLOXOynVvWD7CGYDk3XTmKGtoc1DaW0Tf57yBzy42eucH+mQhoRSx2ySfaUUdC6j2ZZFuy0VEeHsiVn8/ZM9vLO1mkun5/TYRKRUT4JtGtoIZPflxiKSJyLvisgWEdkkInd2UUZE5EERKRKRQhGZ2Zf3UOG1ubyRa9PeZKx9M3tcp+KztqRQoSY2drjOw4+d8R1vIsbaoyAp1smpYzPYW9vKporGCAepBrJgE0EGsFlE3hSRlw88ernGC/yXMWYSMBe4Q0SOXIryAmBc4HEb8EgfYldh5PcbKqt3872cv7LXOZd9du2EDKdOWwJFMWcTb/Yx0vPRweNTRySTmxLLB9v30dyhm9ioYxNsIrgPuBT4KZ8tRf3rni4wxlQYY9YEnjcBW4ARRxS7BHjKWD4GUkRkeNDRq7DZXdvCnemPE2PrZGnCD3XiWATU20dT7jiR4d4NjG9/DbAmmp09KQufMXy4Y1+EI1QDVVCJwBjzHrAbcAaefwqsCfZNRGQUMAP45IhTI4CSQ16XcnSyQERuE5FVIrKqpqYm2LdV/cjsW8kVqe+yOvZm6hwFkQ5nyNrrnE+jbTjnNX2fdO82AFLiXMzKT2VbVRNldW0RjlANRMGuNfQV4Dngj4FDI4AXg7w2AXge+FZgmYrDTndxyVEd0saYx4wxs40xszMzdUGzcOv0+PiS+2Hq/amsivtqpMMZ0ozY2R6zkA5bAhc3fJ0Yv/UrNXtUKoluB+9ur9YVSlWfBds0dAdwCtAIYIzZAWT1dpGIOLGSwD+62cimFMg75HUuUB5kTCpMkve/wknxm1jivAOPLSHS4Qx5Honn1aQHSfRXsLDp22D8OO02Th+Xyf7mTgp1boHqo2ATQYcxpvPACxFx0MUn90OJNZbtCWBLDxvYvAzcGBg9NBdoMMZUBBmTCgO76eQi8zu2d4yiNOWaSIejAiqcM1mWcC8FncuY2/p7AMZkxpOfFseK4v20dfp6uYNSnwk2EbwnIvdibWJ/LvAs8Eov15wC3ACcJSLrAo8LRWSRiCwKlFkMFANFwOPA7X3/ElQoTW7+GyMc5Tzt+SbYdNJSNCl0X8ummMuZ1/owBR1LERFOH5eBx+tn5S5dF1IFL9gJZd8FbgU2AF/F+gP+p54uMMZ8SNd9AIeWMVjNTioKuf21zG/7A8uaZtGeeTZJkQ5IHU6EdxJ/SIZvGwubvs3T9mchoYApOUkUltUzLS850hGqASLYUUN+rM7h240xVxpjHtdZxoPf3NY/EEMrD9XfRmaiTh6LRj5x80rSw/hwcXHj7bj8TcwtSMduE5YX7Y90eGqA6DERBNru7xORfcBWYJuI1IjID8ITnoqUVG8x09qe5una87GnTtXlC6JYkz2HV5MeJNlXwsKmbxPvsjFrZCpFNc2s2q1NRKp3vdUIvoXV1n+SMSbdGJMGnAycIiJ3hTo4FTmntfySTuPigcrrmTAsMdLhqF6UuU7ivYR7GdP5LvNaH2JmfirxMXZ+ungLWnlXvektEdwIXGuM2XXggDGmGPhi4JwahPI6VzCmcyl/rr8WW+ww0uJ1memBYL37Oja6r2Bu6x+Y5H2bk0ens2ZvPcu26SRM1bPeOoudxpij5q0bY2oCcwTUQFd43+GvjZ8LGh+hlSQeLP0cl+cVMrflpUhEpvpKhKUJPyTdW8T5jf/N/sxn2FEdz6+XbGPBhExt3lPd6q1G0HmM59RAVV9IvNnHi80X02FczEvXaR0DiU9ieCXpITy2eC5t+jr3LMhiY1kjb26qinRoKor1lghOFJHGLh5NwNRwBKjCyN8JVe/QZBvGI5UXMj6hjswY3Zp6oGmxD+OVwMzji5rvZmymmweWbNelJ1S3emwaMsbYwxWIigI1y8HbzEoupaQ1iZtGbY50RKobPW16f8Bu16mMqVzCI/l+zl19F6+++hAXX/zNMESnBpq+7FmsBjNPI+xbDkmTeWX/SQiGk9MqIx2VOg7VjimQNptxnne4bfgSfrsxB6/PH+mwVBTSRKAsVe8Cfsywc1ixL5spyftJcWk30ICXvRDi8vnvzD/g7qzgxXW6pqM6miYCBW0VUL8O0udQ2JxLVUc887WTeHCw2SH/KmwON38p+Al/fXclHq0VqCNoIhjqjIHKt8AeC5mn8/LuNOziZ06ajjIZNBwJSP7VZDjq+F7Sj3hl7Z5IR6SiTLCLzqnBqmk7tOyG4RfgEzev7k1jekoN8Q7d/3YwWFF8YL0hNxnOBcxPeJt/rLyb39T+vMt5BXedq3tRD0VaIxjK/B6oXAKudEibxac1iVS1ubRZaJDa55zIUnMt16e8SHr1M5EOR0URTQRD2Y5HoXM/ZJ8HYuflPWnE2n3MTNUlCQardRnf55PW6VzD/5HVWRjpcFSU0EQwVHXWwYb7IH40JI7D4xdeL0nl3Nx63Hbd3WqwEpuTZ5y/pMabyuca7iDOf9QKMmoIClkfgYj8Gfg8UG2MOaGL8wuAl4ADC9q9YIz5cajiUUfYeL+VDPKuAhE+LE+irsPJRSNrdfGQQWxuy0PMThS+s+Ue/jzyf7im7gtsjrkEJPCZsDC96wun3Re2GFX4hbJG8CSwsJcyHxhjpgcemgTCpWknbH8ICm6G2GwAXtmTRpLTy+nDdePzwc5hM0xMN/xP2e0k+8vI86yMdEgqwkKWCIwx7wO6K0Y0WvcdECdM+wkA7V7hzdJULsivI8au69EMBWdklfF286m82XQ6ud5VpPh2RzokFUGR7iOYJyLrReR1EZnSXSERuU1EVonIqpoa7cg8LlXvQckLMPm7EJcDwNLyFFq8di4eqVsbDhUum5/PDd/NN3ffSa3JYmzHElz+xkiHpSIkkolgDTDSGHMi8BDWnshdMsY8ZoyZbYyZnZmZGa74Bh+/D9bcBXG5MOm/Dh5+eU8aGW4Pc7OaIhicCrdzhpXgsAk/qPwWgmF855vg1/kjQ1HEEoExptEY0xx4vhhwikhGpOIZEnY9BXVrYfovwBEHQJPHxtKyFD6fX4s90vVDFVZuu4/zs/fyas00VnIhif4qa5a5GnIi9qsvItkSmNooInMCsWjbRKh4mmH9vZB+Moy89uDhN0pS6fTbuEibhYak87P3EGPz8kj5BZQ7ToTaT6FhY6TDUmEWyuGjTwMLgAwRKQV+CDgBjDGPAlcCXxMRL9AGXGN0l+3Q2fwLaK+E016AQ5YW+M+uDEYmtDMzoyWCwalISXJ6ODOrlCVV+azOPYucuFooewXc2RCjFfShIpSjhq41xgw3xjiNMbnGmCeMMY8GkgDGmIeNMVOMMScaY+YaY5aHKpYhr2UvbP0VjLwGMucdPFzR0MaKqkQuHbUf3c526Pr88N0AvFo5FvKuBHHA3n9bO9apIUEXnRssCu/r/lzJC1YnYMyww8q9tDkbQx6XjdZmoaEsPaadUzPKebc6l/2+WtLzLofdf4fyxZB7aaTDU2Gg3YODXWspNGyAjLngSjns1Iu705mR3syoxI7IxKaixkU5u/D4bTy5bRgkjIGsM6B+PdSti3RoKgw0EQxmxkDFm+BIgIxTDzu1uS6WrfVxXD5a15pRMCK2hZPSqvjr9iyaPDbIPB3iR1m1gg79GRnsNBEMZg2boK0Uss4Ee8xhp17cnY5D/Hx+pE7+VpaLc4pp9Dh4uijLWnso93KwOaHkOfC1Rzo8FUKaCAYrvweq3rZGf6ROP+yUzw8v7U5nQU4DqTG60qiyjElo5JRhDfxp6zA6fALORKuPoL0K1twT6fBUCGkiGKz2rQBPA2Sf/9nKkgHvVSRT1ebiygLtJFaH+9qUSqrbXLywK7AKaeI4SJ8HO34PJf+JbHAqZDQRDEaeJtj3ISROhIRRR51+emcmGW4PZ4+oD3toKrqdMqyRaWkt/HHzcA7ucT/sbEg7CT6+BVp0v+PBSBPBYFS1FIwPss896lR1m5OlZSlcWbAPp03n76nDicDXJlewu9nN4pI066DNDqc+A/jho2utZkc1qGgiGGzaKqB+HaTNgZi0o04/W5yBzwjXjNFVXFXXzs+royCpjUc2Z3Nwrn9CAcx5zGpyLPxhRONT/U8TwWBijLVomD3WGgd+BL+BZ4oymDesUecOqG7ZBBZNqmRzXTzvVSR9dmLk1TDmK7D551CxJHIBqn6niWAwadoGLbshawHY3UedXl6VREmLW2sDqleXjtrP8LhOHtk8/PATs34LyZNgxQ3QVhWR2FT/00QwWPh9ULnEWigsbXaXRZ4uyiTF5eX8vLowB6cGGpfd8OWJlXxSncTqmvjPTjji4JR/g6cRVtwIxt/9TdSAoYlgsKhdCZ21XQ4XBahodfJmidVJ7NbtKFUQrhlTQ4rLe3StIGUKzHzAaobc+pvIBKf6lSaCwaB9H1S/BwljIXFsl0We2j4MP8KXxmt1XgUn3unnpglVvF2WyrbKI3avG3sb5F1u7XGxf1VkAlT9RhPBYLDhh9aSwdnndXm61Wvjn0WZnJ9bR16CLi2sgvel8VXEOXw8+t7Ow0+IwJzHwT3MGlLq0W1OBzJNBANdXSEUPWr1C7i73s/5+V3pNHQ6uHWi1gZU36TG+Lh2TA0vry+npLb18JMxaTD/H9BSDKu+EZkAVb8I5Q5lfwY+D1QbY07o4rwAvwMuBFqBm4wxa0IVz6BkDKy+E1yp1sJyXfAb+Mu2YZyY1sysjOYwB6gGmhXFRy87MiOumb+YbO745xrOnJB18Phd546HrNNhyv/Cxh9ZExhHXx/OcFU/CWWN4ElgYQ/nLwDGBR63AY+EMJbBqeQ5qF4G034Cjtgui7xXnkxxYyy3TKzSXcjUMUmP6WBidhKbyhtp7fQeXeCE/4HMU+DTr0FzcfgDVMctZDUCY8z7IjKqhyKXAE8F9in+WERSRGS4MaYiVDENKt5Wa0XIlGkw5jbY+JMuiz2yJZs0Vztp3iJWFOtoIXVsbs16jXsqTqVm59tcnb/DOliY/lmB1NmwfzUsOR0Kbgaxw7T7IhKr6rtI9hGMAEoOeV0aOKaCseWX0LoXZj1orQXThRVViaysTuKinGIcuq6QOg45gY1r3qrKp9Xbxc+bKwVGXARtZVD1btjjU8cnkomgq4aKLv9aichtIrJKRFbV1OisWFr2wuZfQP4XYNjRS0kc8LsNOWTFdnJWVmkYg1OD1cU5xbT6nLxdnd91geTJkDoT9n2kTUQDTCQTQSmQd8jrXKC8q4LGmMeMMbONMbMzM7seGTOkrP229e+MX3Zb5JPqBD6uTmLRpEpcNp39qY7fmIRGTkjex+sVI+n0d/OnY/j51uz20v9Au35oGygimQheBm4Uy1ygQfsHglC1DPb+Gyb/N8R388kM+N2GEWS4PVw3tjp8salB75KcYuo9bt6v6aYV1+aC3CvA1wYf38xny5eqaBayRCAiTwMrgAkiUioit4rIIhFZFCiyGCgGioDHgdtDFcug4fdaw0XjR8Kk73RbbGV1Asurklg0qQK3Q38RVf+ZklTLmPh6Xikfjbe7imZstjWUtPw12P5QWONTxyaUo4au7eW8Ae4I1fsPSjsfh/pCOPXZboeLGgM/W5tHVmwn14/TqrnqXyJw6Yhifr19Ji/syuALY/Z1XTBtDojDasbMOgNSTwxvoKpPdGbxQNFWBevuhWFnQt4V3RZ7bW8qa/cncM+0MmId2jeg+t+s1GrGxNfz2w051ib3XRGBuX+GmHT46BrwtoQ3SNUnIasRqONQeN/Rx0peAG8jJE2GDT/q8rIOn/CLdblMTGnlitHdfFJT6jiJwNX5O/jplpP4Z1EmN0/oph/KnQnz/gZLz4XVd8HJj4U3UBU0rREMBM3F0LABMk61RmR046ntWZS0uPn+jBLs+j+rQmhq8n7mD2vk4Y05tHh6+GHLPtsa2LDzcdj7bPgCVH2ify6ind9rdbq5UiHztG6L1bY7eGhjDguG13Pa8MYwBqiGqntOLGV/h5O/bBvWc8FpP4b0OfDJV6BlT3iCU32iiSDa1XxobTiT8zmwdd+S9/N1ubR6bdw7o6TbMkr1p5kZLZwzoo4/bslmX3sPrcw2J8z/p7Wb2fLrrQ83KqpoH0E0a6uEmg8geSokjDl4+MgVIrc2pvLv4kwuzilmf20ZK2rDHagaqr47vZTzF5/AbwpH8NM5PXzaTxwDcx61EsG678LMX4UvSNUrrRFEK78Pyl60hokO734RV69feGLXZDJcbVw2Yme35ZQKhbHJ7dwwrppndmaytb7rIc0HjboOxt0BW38Nu58JT4AqKJoIolXNB9BeBTmftzYM78biilGUtiVy0+jNuO2+MAaolOXOqWUkOn3cvyav94nEM39jLVn9yS3WpkoqKmgiiEZtFYEmoWmQNLHbYhVtcTxfNpbZqVXMStXJYyoyUmN83HlCOR9WJrO0PLnnwnYXnPqctVrpB5dBh7ZjRgNNBNHG22It2OWI67FJyG/gkZ1TcYqPm0dvDmOASh3thvHVFCS18ePV+bR7e9kBKTYbTn0eWksCncdak400TQTRZtXXoaMGci/tdhkJgFcrRrOjOZWbRm8hzdURvviU6oLTZrh/9h72NLt5cGNO7xdkzoNZD0HFG1D4v6EPUPVIRw1Fk+InrUfm6YeNEjpSSWsCz5aMY05aJaek64KtKjrMz27iitH7eGxLNpeMqmVC4X09X2CMtX/B5p9BwyZInWEd153Nwk5rBNGifhN8eru1llBW95vNePzCIzunEmv3csvozboPsYoq359RQqLLx70rR+LvreNYBHIutD70lL0KTTrqLVI0EUSDjv3wweXgTLIm3kj3/y1/2DScXS3JfLlgE8nOzjAGqVTv0txevj+jhNX7Evn7jqzeLxA75F1lrUtU8m9r7owKO00EkeZrh/cvsaben/qc1ZHWjfX743lo43BOSS9nTlpVGINUKnhXjN7P6cMb+OnaXHY2unu/wB4DI68DWwzs+YfWDCJAE0EkGT+suBFqPoL5f4OsU7st2uSx8Y2PChgW6+EmHSWkopgI/PLkXbjthruWF+DxB9F+6UyCUV8E44OlZ1v7cquw0c7iSDEG1txjrcg441eQf1WPRe9dOYqylhj+dc5WPE26VouKvCOXOjnSzaMMD2yfwXc+SOELeUXMK0jv+YbuLBh1A+z9Fyw9B855v8casuo/Ia0RiMhCEdkmIkUi8t0uzi8QkQYRWRd4/CCU8UQNY2D1t2DbAzD+mzDx7h6LP1ucwSt70rl7WhmzM5vDE6NSx2lOWhVnZJbyYtkYNjemBndR7HBY8Dq0lcM7Z0FreWiDVEBo9yy2A78HLgAmA9eKyOQuin5gjJkeePw4VPFEDeOHlV+F7Q/ChG/BrN/S09CfogY3P1yVz/xhjSyapENF1cDypVFbGO5u4Xc7plPZ6gzuosx5cMZr1oSzt0+H5t0hjVGFtkYwBygyxhQbYzqBZ4BLQvh+0c/XDstvsDbpmPJ9a92VHpJAu0/4+kdjiHP4eWBesW42owacWLuPuyespdNn52sfjO1+a8sjDTsDznrbGlH39mnQuC20gQ5xofzTMgI4dHH80sCxI80TkfUi8rqITOnqRiJym4isEpFVNTUDdE2dlhJYchrs+Sec+DM48f4ek4Ax8INPR7K1Po5fzdvFsDhPGINVqv+MiG1h0ZgNrN2fwI9W5wd/YcbJcM4y8HVYvzv7Pg5ZjENdKDuLu/ord+QUkzXASGNMs4hcCLwIjDvqImMeAx4DmD17dm/TVKJP9fvw4VXgbYPTX4TcoytGDyzZfvD53Jb9vFOVy793ZXLZiCLc7cWsKA5jvEr1s5PTq1jkrODRzcMpSGzny5N6GP5ceN/hr/Ovht3/sJJB7mWQ3FULMzoj+TiEskZQCuQd8joXOKznxxjTaIxpDjxfDDhFpPtNeQcaXyes/z68cyY4U+D8T7pMAkfa0ZTMX3ZP5sTkGq7MLQp9nEqFwbenlXJhXi33r83n5d1pwV8Ykw5jbgV3NpQ8C/uW0/t616ovQpkIPgXGichoEXEB1wAvH1pARLJFrPYREZkTiKfnMWkDRf0GeHMObPopjP4SLPwUkif1ellLh5ff7phBuqudr48rxKZLSKhBwm6D38wvZk5mE/d8PJrlVYnBX+yIh9E3QtJkqFxirdDr1+bS/hKyRGCM8QJfB94EtgD/NsZsEpFFIrIoUOxKYKOIrAceBK4xZoCnem+LtRXf6zOhvQJOfwnm/tmaMNMLn9/w+sZKmr1O7h6/lgSH/qCrwcVtNzx+xg5GJnbwlffGsaomIfiLbU7IuxKyzoKGDVD8BHTWhS7YIUQG2t/d2bNnm1WrVoX3TQvvC65c4zaoeB08DZAyHbLPsT7JBGFF8X7+smsSb1WN5PYx6zktU4eKqsHl0AllVa1Orn1nAlVtLp5csJ2Tsvo4P6apCEqft3od8y6HxHHaR9ALEVltjJnd1TkdkNgfOuthzzOw9xmwuWD0TVZfQJBJAOCNynzeqhrJhcN3aRJQg96wOA/PnLONYXGdfGnZeD6p7kPNACBxLIy5DVzJ1ki8qnfBrzPuj5UuMXE8jA/2rbBGBQEMOwcy5lorKvbB0rJknto9jtmpVVyfr+Ol1eDU1ZIU94yt4/7NJ/HFd8Zzx9hCTk63RhP1uhwFgCsVCm6F8sVQ8z68fQbM/wckjOrnyAc/rREcq5Y9UPRHqHrHWk993B3Wptx9TAIba+P4+kdjGBXfyB1jtXNYDS2prg7um/IJo+Mb+d2O6SyuGNm3G9icVu0793Jo2AivT4c9/w5JrIOZJoK+8rZA6Uuw60nwd0L+tTDyaquK2kdFDW5ufHc8qS4v90xYg9uue7eqoSfR6eH7kz/lpLQq/rZnEn8qnhL8DOQDUqbCBesgaRJ8dDV8fCt4GkMS72CkiSBYxkDtGtjxe6gvhIxTrFpA0vhjul1pi4sb3p2ATeDvZ23TfYfVkOay+blz3Douzinmneo8vrBkImUtrr7dJGE0nPs+TLnX+qD22lSoWBKSeAcbTQTBaK+CXX+B8lcgJgvGftUaEWQLchGtI1S2Orlh6QSaPTaeOnMbo5M0CShlE7g2fzt3j19DcZObz78xmTdKgly19OBNnHDi/8E5H4IjFt49z1rk0dMUmqAHCU0EPfE0W3sGFP3RWvxqxCXW5DB3EFvwdaOk2cVVSyZS0+bkyQU7mJza1o8BKzXwnZRWzUvnb2ZEXCeLPhjLXctH09DZt743MufBwrUw6R4oehxeOwHK3wxNwIOAziPoijFQ+iKs/ia0lkLqTBh2Njjijuu2RQ1uvrh0Am0+G0+duZ0T01sOnuttkw+lhpJ5Bel4/MLDG4fz8KYc0t0evj+jhItH1va0VmPXWkusfr3O/VYfQvb5h/fpDZH5BzqPoC+aiuC9i6zN5F1pcO5HMOKi404CK6oSuertiXiN8MzZWw9LAkqpozlthrumlfPi+ZsZFuvhzuVjuPqdCWyui+3bjeLyYOwiyDoTmnZY/Xw1H4JfB2ccoIngAG8LrP8feG0KVL9nbR+5cDVkzj/uW/99RyY3LB1PWoyX587dwiRtDlIqaFPTWnnxvM38bM4udjTEcuHrJ3DHh2MoanAHfxObA7JOtwZ4JBRYw753PmqtBjDAWkVCQSeUGQMlz8Oau60q5Kgvwoz/Z22Zd5xaPDZ+siafZ3ZmcmZOPb+bX0ySSz+FKNVXdhtcO3YfF+bX8fiWbP6ybRiL96ZyQV4dt06sYmZGc3BNRq4UGHkNNG2Hijet1QDaSmD6L6x5QEPU0E4EVctg/b3W7OCUadasxKzT+uXWq2sSuHvFaPY2x/C1yRXcM61UdxhT6jglu3zcc2IZt0yo4vGt2fxjRyaLS9KYltbCNWNr+Fx+LcnBfNhKHG9NBK1bC7WrYcmpVj/glHth2Jk9bho1GA29zmJjrPbBTf9nfSKIHQFTfwAFt1jVx64U3tent3hqexY/XJVPRkwbXxuzgUlJukKiUqHQ7rPzwb4cPqgtYEdDLC6bn7NH1LMwr44FOQ3BJYXJ34btf4Ctv4H2Skg/GSZ8E/KuAHtM6L+IMOmps3jo1Aj8PmsN8y2/gv2fWJtdzPil1Wbo6GPnUy9mZDSzIKuUL+ZvJc6hTUFKhYrb7uPcYSX879xWNtbF8XxxBq/uTeP1kjTsYpiV0cycrCZmZzYxM6Ol66ZZR7yVDCZ8A4qfhC2/huXXQ8ydUHCz1VycMnVQ1xKGTo1g5xPwyZet6uDEu6HgpuBHAhXe1+e30+GgSoXPoYvU+Q2s3x/PktIUPqpMYmNdPD4jCIYJKW1MSG5jfEob45LbGJfURs7c7xDjOGSegvFD5TtQ9Kg17NT4IGkC5F0FORdC+kndtx5EsZ5qBEMnEXhbrKagEZeArY+TUwrv6/PbaSJQKjq0++wUNSeztSmVoqYUytoS2Nd5eCtAnMtOQoyDRLeDxBgn8TF2Yl12Mp0NfG3MBjLrXsJd+wGCHxyJkHWGtdJw6gxIm2FtoxnlNYaIJQIRWQj8DrADfzLG/PyI8xI4fyHQCtxkjFnT0z2Pq4+g8L5ju+4YaCJQKnq1+eyUtSVQ1hbPOv8Cmjq8NLd7D/7b6fMfdU2KvZHTkjZyZkohc2LXk+soOXiug3iaHHm0ufLwx2Rhi0nF7k7DGZtOTGwKLqcLh8OBfe9ziNgCScN2SPIQ63HoazlwzA6TvwOJY47ra45IH4GI2IHfA+dibWT/qYi8bIzZfEixC4BxgcfJwCOBf5VSKmRi7T7GJjQwNqGBM/jnUefbfXaaPC4avU4aPS4avS7rtSebl5rzeajpvxBPIyPtOxjrKCLHUUaeq5I81zbS7J+SZG8mxtZ/G+Ws2biCRZW/4fYFY7jplNH9dt8DQtnQNQcoMsYUA4jIM8AlwKGJ4BLgqcA+xR+LSIqIDDfG6BZdSqmIcdt9uO1tZNL15M+P488EcoCJANQZQ43fsMLr5+qT8tjT2klzWyPtzfvpaKun0+vB6/HgqVmJzxi8PmuP8poWL2AwxmDVBQyCwXCgjmBw4KVCxhAf4yA+JjR/skOZCEYAJYe8LuXoT/tdlRkBHJYIROQ24LbAy2YRifQ2XhnAvgjH0JuBECMMjDg1xv4zEOIMIsaHuz3zg/6N5TDLDn/Z1+9lt7v+hDIRdNVzcmSHRDBlMMY8BjzWH0H1BxFZ1V1bW7QYCDHCwIhTY+w/AyHOgRAj9G+coZzrWgrkHfI6Fyg/hjJKKaVCKJSJ4FNgnIiMFhEXcA3w8hFlXgZuFMtcoEH7B5RSKrxC1jRkjPGKyNeBN7GGj/7ZGLNJRBYFzj8KLMYaOlqENXz05lDF08+ippmqBwMhRhgYcWqM/WcgxDkQYoR+jHPATShTSinVv3Q9TKWUGuI0ESil1BCniaAbIrJQRLaJSJGIfLeL8xNFZIWIdIjIPZGIMRBHb3FeLyKFgcdyETkxCmO8JBDfOhFZJSKnhjvGYOI8pNxJIuITkSvDGV/gvXv7Xi4QkYbA93KdiIRyWPsxxXhInOtEZJOIvBfuGAMx9Pa9/PYh38eNgf/ztCiLMVlEXhGR9YHv5bH1sxpj9HHEA6tzeydQALiA9cDkI8pkAScB/wfcE8VxzgdSA88vAD6JwhgT+Ky/ahqwNRq/l4eUW4o10OHKaIsRWAC8Gomfxz7EmIK1wkB+4HVWNMZ5RPmLgKXRFiNwL/CLwPNMoBZw9fW9tEbQtYPLYxhjOoEDy2McZIypNsZ8CngiEWBAMHEuN8Yc2BnnY6y5GtEWY7MJ/CQD8XQxqTAMeo0z4BvA80B1OIMLCDbGSAomxuuAF4wxe8H6XQpzjND37+W1wNNhiewzwcRogMTAAp4JWImgz4scaSLoWndLX0SbvsZ5K/B6SCM6WlAxishlIrIVeA24JUyxHarXOEVkBHAZ8GgY4zpUsP/f8wJNBa+LyJTwhHZQMDGOB1JFZJmIrBaRG8MW3WeC/t0RkThgIdYHgHAKJsaHgUlYE3E3AHcaY45eOrUXA293hfAIaumLKBB0nCJyJlYiCHf7e7DLiPwH+I+InA78BDgn1IEdIZg4fwv8tzHGJ5FZez6YGNcAI40xzSJyIfAi1uq+4RJMjA5gFnA2EAusEJGPjTHbQx3cIfryO34R8JExpjaE8XQlmBjPB9YBZwFjgCUi8oExprEvb6Q1gq4NlKUvgopTRKYBfwIuMcaEe6OEPn0vjTHvA2NEJCPUgR0hmDhnA8+IyG7gSuAPInJpWKKz9BqjMabRGNMceL4YcIb5exns0jJvGGNajDH7gPeBcA9i6MvP5TWEv1kIgovxZqxmNmOMKQJ2cWBJ1L4IdyfNQHhgfWIpBkbzWSfNlG7K3kfkOot7jRPIx5q5PT+KYxzLZ53FM4GyA6+jKc4jyj9J+DuLg/leZh/yvZwD7A3n9zLIGCcB7wTKxgEbgROi7XsZKJeM1e4eH874+vC9fAS4L/B8WOB3J6Ov76VNQ10wQSyPISLZwCogCfCLyLewevT7VCULdZxYq+KmY316BfCaMK6sGGSMV2CtOeUB2oCrTeAnO8rijKggY7wS+JqIeLG+l9eE83sZTIzGmC0i8gZQCPixdi/cGK4Yg40zUPQy4C1jTEs44+tDjD8BnhSRDVhNSf9trFpWn+gSE0opNcRpH4FSSg1xmgiUUmqI00SglFJDnCYCpZQa4jQRKKXUEKeJQKleiEiKiNweeJ4jIs9FOial+pMOH1WqFyIyCmtFzxMiHYtSoaATypTq3c+xlr1YB+wAJhljThCRm4BLsSb7nAD8GmsG6A1AB3ChMaZWRMYAv8daJrgV+IoxZmu4vwiluqNNQ0r17rvATmPMdODbR5w7AWtZ5TlYe1O0GmNmACuAA6tqPgZ8wxgzC7gH+EM4glYqWFojUOr4vGuMaQKaRKQBeCVwfAMwTUQSsDYHevaQFUtjwh+mUt3TRKDU8ek45Ln/kNd+rN8vG1AfqE0oFZW0aUip3jUBicdyYWARwl0ichWAWMK+b7RSPdFEoFQvjLWHw0cishH45THc4nrgVhFZD2wi+raXVEOcDh9VSqkhTmsESik1xGkiUEqpIU4TgVJKDXGaCJRSaojTRKCUUkOcJgKllBriNBEopdQQ9/8BM+4QwlGWvn8AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# clustered distribution in blue\n", "sns.histplot(\n", " pd.DataFrame(my_outcomes_clustered)[\"time\"], kde=True, stat=\"density\", linewidth=0\n", ")\n", "\n", "# non-clustered distribution in orange\n", "sns.histplot(\n", " pd.DataFrame(my_outcomes)[\"time\"],\n", " kde=True,\n", " stat=\"density\",\n", " linewidth=0,\n", " color=\"orange\",\n", ");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Calculating standard errors and confidence intervals from existing bootstrap result" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you've already run ``bootstrap`` once, you can simply pass the existing result object to a new call of ``bootstrap``. Estimagic reuses the existing bootstrap outcomes and now only draws ``n_draws`` - ``n_existing`` outcomes instead of drawing entirely new ``n_draws``. Depending on the ``n_draws`` you specified (this is set to 1_000 by default), this may save considerable computation time. \n", "\n", "We can go on and compute confidence intervals and standard errors, just the same way as before, with several methods (e.g. \"percentile\" and \"bc\"), yet without duplicated evaluations of the bootstrap outcome function. " ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(constant 90.709236\n", " time 0.151193\n", " dtype: float64,\n", " constant 96.827145\n", " time 0.627507\n", " dtype: float64)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "my_results = em.bootstrap(\n", " data=df,\n", " outcome=ols_fit,\n", " existing_result=result,\n", ")\n", "my_results.ci(ci_method=\"t\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can use this to calculate confidence intervals with several methods (e.g. \"percentile\" and \"bc\") without duplicated evaluations of the bootstrap outcome function." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Extending bootstrap results with more draws\n", "\n", "It is often the case that, for speed reasons, you set the number of bootstrap draws quite low, so you can look at the results earlier and later decide that you need more draws. \n", "\n", "In the long run, we will offer a Dashboard integration for this. For now, you can do it manually.\n", "\n", "As an example, we will take an initial sample of 500 draws. We then extend it with another 1500 draws. \n", "\n", "*Note*: It is very important to use a different random seed when you calculate the additional outcomes!!!" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(constant 90.768859\n", " time 0.137692\n", " dtype: float64,\n", " constant 96.601067\n", " time 0.607616\n", " dtype: float64)" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "initial_result = em.bootstrap(data=df, outcome=ols_fit, seed=5471, n_draws=500)\n", "initial_result.ci(ci_method=\"t\")" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(constant 90.689112\n", " time 0.128597\n", " dtype: float64,\n", " constant 96.696522\n", " time 0.622954\n", " dtype: float64)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "combined_result = em.bootstrap(\n", " data=df, outcome=ols_fit, existing_result=initial_result, seed=2365, n_draws=2000\n", ")\n", "combined_result.ci(ci_method=\"t\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Using less draws than totally available bootstrap outcomes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You have a large sample of bootstrap outcomes but want to compute summary statistics only on a subset? No problem! Estimagic got you covered. You can simply pass any number of ``n_draws`` to your next call of ``bootstrap``, regardless of the size of the existing sample you want to use. We already covered the case where ``n_draws`` > ``n_existing`` above, in which case estimagic draws the remaining bootstrap outcomes for you.\n", "\n", "If ``n_draws`` <= ``n_existing``, estimagic takes a random subset of the existing outcomes - and voilà! " ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(constant 90.619182\n", " time 0.130242\n", " dtype: float64,\n", " constant 96.557777\n", " time 0.625645\n", " dtype: float64)" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "subset_result = em.bootstrap(\n", " data=df, outcome=ols_fit, existing_result=combined_result, seed=4632, n_draws=500\n", ")\n", "subset_result.ci(ci_method=\"t\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Accessing the bootstrap samples" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is also possible to just access the bootstrap samples. You may do so, for example, if you want to calculate your bootstrap outcomes in parallel in a way that is not yet supported by estimagic (e.g. on a large cluster or super-computer)." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "ename": "TypeError", "evalue": "get_bootstrap_samples() got an unexpected keyword argument 'seed'", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", "Input \u001b[0;32mIn [16]\u001b[0m, in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mestimagic\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01minference\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m get_bootstrap_samples\n\u001b[0;32m----> 3\u001b[0m my_samples \u001b[38;5;241m=\u001b[39m \u001b[43mget_bootstrap_samples\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdata\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdf\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mseed\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1234\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m 4\u001b[0m my_samples[\u001b[38;5;241m0\u001b[39m]\n", "\u001b[0;31mTypeError\u001b[0m: get_bootstrap_samples() got an unexpected keyword argument 'seed'" ] } ], "source": [ "from estimagic.inference import get_bootstrap_samples\n", "\n", "my_samples = get_bootstrap_samples(data=df, seed=1234)\n", "my_samples[0]" ] } ], "metadata": { "interpreter": { "hash": "5cdb9867252288f10687117449de6ad870b49795ca695c868016dc0022895cce" }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.13" } }, "nbformat": 4, "nbformat_minor": 4 }