Skip to content

Instantly share code, notes, and snippets.

@jonathan-taylor
Created May 29, 2019 05:43
Show Gist options
  • Save jonathan-taylor/1af6485a97f6b526eae3a5a35236e4cd to your computer and use it in GitHub Desktop.
Save jonathan-taylor/1af6485a97f6b526eae3a5a35236e4cd to your computer and use it in GitHub Desktop.
ROSI fails under model misspecification
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The rpy2.ipython extension is already loaded. To reload it, use:\n",
" %reload_ext rpy2.ipython\n"
]
}
],
"source": [
"import numpy as np\n",
"from selection.algorithms.lasso import ROSI\n",
"from scipy.stats import norm as normal_dbn\n",
"import regreg.api as rr\n",
"%load_ext rpy2.ipython\n",
"import seaborn as sns\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>variable</th>\n",
" <th>pval</th>\n",
" <th>lasso</th>\n",
" <th>onestep</th>\n",
" <th>sd</th>\n",
" <th>lower_confidence</th>\n",
" <th>upper_confidence</th>\n",
" <th>lower_truncation</th>\n",
" <th>upper_truncation</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0</td>\n",
" <td>0.207549</td>\n",
" <td>1.112094</td>\n",
" <td>1.271107</td>\n",
" <td>0.092361</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>-0.067633</td>\n",
" <td>0.141366</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1</td>\n",
" <td>0.014068</td>\n",
" <td>1.210254</td>\n",
" <td>1.379464</td>\n",
" <td>0.091534</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>-0.050427</td>\n",
" <td>0.154843</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>2</td>\n",
" <td>0.202229</td>\n",
" <td>1.081264</td>\n",
" <td>1.267556</td>\n",
" <td>0.088498</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>-0.020597</td>\n",
" <td>0.171282</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td>8</td>\n",
" <td>0.568667</td>\n",
" <td>0.021850</td>\n",
" <td>0.097024</td>\n",
" <td>0.075877</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>-0.066078</td>\n",
" <td>0.074976</td>\n",
" </tr>\n",
" <tr>\n",
" <th>10</th>\n",
" <td>10</td>\n",
" <td>0.702063</td>\n",
" <td>0.002883</td>\n",
" <td>0.080366</td>\n",
" <td>0.072552</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>-0.051511</td>\n",
" <td>0.077452</td>\n",
" </tr>\n",
" <tr>\n",
" <th>13</th>\n",
" <td>13</td>\n",
" <td>0.379092</td>\n",
" <td>0.029393</td>\n",
" <td>0.114922</td>\n",
" <td>0.076726</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>-0.059319</td>\n",
" <td>0.084909</td>\n",
" </tr>\n",
" <tr>\n",
" <th>14</th>\n",
" <td>14</td>\n",
" <td>0.078252</td>\n",
" <td>-0.083292</td>\n",
" <td>-0.161437</td>\n",
" <td>0.073898</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>-0.076718</td>\n",
" <td>0.057073</td>\n",
" </tr>\n",
" <tr>\n",
" <th>18</th>\n",
" <td>18</td>\n",
" <td>0.322493</td>\n",
" <td>0.051861</td>\n",
" <td>0.122945</td>\n",
" <td>0.077180</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>-0.075555</td>\n",
" <td>0.070382</td>\n",
" </tr>\n",
" <tr>\n",
" <th>19</th>\n",
" <td>19</td>\n",
" <td>0.021591</td>\n",
" <td>-0.104132</td>\n",
" <td>-0.196766</td>\n",
" <td>0.074534</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>-0.091214</td>\n",
" <td>0.044888</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" variable pval lasso onestep sd lower_confidence \\\n",
"0 0 0.207549 1.112094 1.271107 0.092361 NaN \n",
"1 1 0.014068 1.210254 1.379464 0.091534 NaN \n",
"2 2 0.202229 1.081264 1.267556 0.088498 NaN \n",
"8 8 0.568667 0.021850 0.097024 0.075877 NaN \n",
"10 10 0.702063 0.002883 0.080366 0.072552 NaN \n",
"13 13 0.379092 0.029393 0.114922 0.076726 NaN \n",
"14 14 0.078252 -0.083292 -0.161437 0.073898 NaN \n",
"18 18 0.322493 0.051861 0.122945 0.077180 NaN \n",
"19 19 0.021591 -0.104132 -0.196766 0.074534 NaN \n",
"\n",
" upper_confidence lower_truncation upper_truncation \n",
"0 NaN -0.067633 0.141366 \n",
"1 NaN -0.050427 0.154843 \n",
"2 NaN -0.020597 0.171282 \n",
"8 NaN -0.066078 0.074976 \n",
"10 NaN -0.051511 0.077452 \n",
"13 NaN -0.059319 0.084909 \n",
"14 NaN -0.076718 0.057073 \n",
"18 NaN -0.075555 0.070382 \n",
"19 NaN -0.091214 0.044888 "
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def simulate_logit(n=1000, p=20, s=3, signal=2, lam=10):\n",
" \n",
" X = np.random.standard_normal((n, p))\n",
" beta = np.zeros(p)\n",
" beta[:s] = signal / np.sqrt(s)\n",
" \n",
" linpred = X.dot(beta)\n",
" pi = np.exp(linpred) / (1 + np.exp(linpred))\n",
" Y = np.random.binomial(1, pi)\n",
" \n",
" L = ROSI.logistic(X, Y, lam, approximate_inverse=None)\n",
" L.fit(solve_args={'min_its':200, 'tol':1.e-12})\n",
" return L.summary(truth=beta)\n",
" \n",
"simulate_logit()"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"P = []\n",
"for _ in range(100):\n",
" P.extend(simulate_logit()['pval'])\n"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAYAAAB91L6VAAAEGWlDQ1BrQ0dDb2xvclNwYWNlR2Vu\nZXJpY1JHQgAAOI2NVV1oHFUUPrtzZyMkzlNsNIV0qD8NJQ2TVjShtLp/3d02bpZJNtoi6GT27s6Y\nyc44M7v9oU9FUHwx6psUxL+3gCAo9Q/bPrQvlQol2tQgKD60+INQ6Ium65k7M5lpurHeZe58853v\nnnvuuWfvBei5qliWkRQBFpquLRcy4nOHj4g9K5CEh6AXBqFXUR0rXalMAjZPC3e1W99Dwntf2dXd\n/p+tt0YdFSBxH2Kz5qgLiI8B8KdVy3YBevqRHz/qWh72Yui3MUDEL3q44WPXw3M+fo1pZuQs4tOI\nBVVTaoiXEI/MxfhGDPsxsNZfoE1q66ro5aJim3XdoLFw72H+n23BaIXzbcOnz5mfPoTvYVz7KzUl\n5+FRxEuqkp9G/Ajia219thzg25abkRE/BpDc3pqvphHvRFys2weqvp+krbWKIX7nhDbzLOItiM83\n58pTwdirqpPFnMF2xLc1WvLyOwTAibpbmvHHcvttU57y5+XqNZrLe3lE/Pq8eUj2fXKfOe3pfOjz\nhJYtB/yll5SDFcSDiH+hRkH25+L+sdxKEAMZahrlSX8ukqMOWy/jXW2m6M9LDBc31B9LFuv6gVKg\n/0Szi3KAr1kGq1GMjU/aLbnq6/lRxc4XfJ98hTargX++DbMJBSiYMIe9Ck1YAxFkKEAG3xbYaKmD\nDgYyFK0UGYpfoWYXG+fAPPI6tJnNwb7ClP7IyF+D+bjOtCpkhz6CFrIa/I6sFtNl8auFXGMTP34s\nNwI/JhkgEtmDz14ySfaRcTIBInmKPE32kxyyE2Tv+thKbEVePDfW/byMM1Kmm0XdObS7oGD/MypM\nXFPXrCwOtoYjyyn7BV29/MZfsVzpLDdRtuIZnbpXzvlf+ev8MvYr/Gqk4H/kV/G3csdazLuyTMPs\nbFhzd1UabQbjFvDRmcWJxR3zcfHkVw9GfpbJmeev9F08WW8uDkaslwX6avlWGU6NRKz0g/SHtCy9\nJ30o/ca9zX3Kfc19zn3BXQKRO8ud477hLnAfc1/G9mrzGlrfexZ5GLdn6ZZrrEohI2wVHhZywjbh\nUWEy8icMCGNCUdiBlq3r+xafL549HQ5jH+an+1y+LlYBifuxAvRN/lVVVOlwlCkdVm9NOL5BE4wk\nQ2SMlDZU97hX86EilU/lUmkQUztTE6mx1EEPh7OmdqBtAvv8HdWpbrJS6tJj3n0CWdM6busNzRV3\nS9KTYhqvNiqWmuroiKgYhshMjmhTh9ptWhsF7970j/SbMrsPE1suR5z7DMC+P/Hs+y7ijrQAlhyA\ngccjbhjPygfeBTjzhNqy28EdkUh8C+DU9+z2v/oyeH791OncxHOs5y2AtTc7nb/f73TWPkD/qwBn\njX8BoJ98VQNcC+8AAEAASURBVHgB7d0HmBRF2sDxl5xzTnoERSUICod4iKDAoQKKiChJMZAEJYgo\noigGclKCBEXhFJFwgGRBFBQUQUCJgkRBQCRIUOJ8+9Z3u26Y3Z2ZnVDd8+/nWZjpWPWr2Xm3qyuk\n88QswoIAAggggAACYRVIH9arcTEEEEAAAQQQMAIEYD4ICCCAAAIIRECAABwBdC6JAAIIIIAAAZjP\nAAIIIIAAAhEQIABHAJ1LIoAAAgggQADmM4AAAggggEAEBAjAEUDnkggggAACCBCA+QwggAACCCAQ\nAQECcATQuSQCCCCAAAIEYD4DCCCAAAIIRECAABwBdC6JAAIIIIAAAZjPAAIIIIAAAhEQIABHAJ1L\nIoAAAgggQADmM4AAAggggEAEBAjAEUDnkggggAACCBCA+QwggAACCCAQAQECcATQuSQCCCCAAAIE\nYD4DCCCAAAIIRECAABwBdC6JAAIIIIAAAZjPAAIIIIAAAhEQIABHAJ1LIoAAAgggQADmM4AAAggg\ngEAEBAjAEUDnkggggAACCGSEAAEEokPg/PnzsmTJEsmWLZvUr1/fZPrs2bOyePFiOXjwoFSsWFHu\nuOMOs37SpElSqlQp+fe//y0bN26U/fv3xyFlypRJcuXKJVWqVJGcOXOa9WvXrpVVq1ZJjx49JF26\ndHH78gIBBFIQ8LAggEBUCBw+fNgT81Xgueqqq+Ly26BBA7NO19etW9esjwmmZt28efPM+3bt2sXt\no/vF/sQEaM/nn39u9tm8ebNZP3v27Lhz8wIBBFIWSKebY36hWBBAwOUCR44ckaJFi0pMAJZ9+/bJ\npUuXzN2wZnv79u2SP39+yZcvn7k71rveQ4cOid7tPvbYYzJ58mRp3LixxARp0bvmuXPnyrp166R0\n6dKye/duI3fTTTfJxYsXZdOmTZI+PU+3XP5xIntBEKAKOgiInAKBUAj8+eef8u6778rWrVtNoGvV\nqpUUL148waW+/fZbmT9/vui+Wn1cp04dyZ49e9w+s2bNki+++MJUJ997771x63///Xd5++23TRDW\nauSZM2fK3XffbYLpsmXLpGvXrib4xh0Q86JWrVrSvXt3s6pbt25SuHBh2bNnj2zbtk2uv/56adu2\nrdmu12zevHn8Q3mNAAJeBAjAXlBYhUCkBU6fPi233HKLCb76zFYD7KBBgySmylcqV65skqcBVAPh\nlStXJGPGjDJs2DAThJcvX262P/vss2advsmaNauMGjXKrNd/9G741VdfNe/PnDkjzz//vBQqVEj0\nObEuMVXT5v/k/jl58qS529XtmTNnNrvFPlfW6xOAk5NjPQJ/C1BP9LcFrxCwRmDw4MEm+LZp00aO\nHz8u06ZNE71rffHFF00aNYD27t3bNHjSO9xTp07FNZjS91rFPHLkSFMVrO+PHj1q7mBjM1i+fHn5\n6aefzFsNvHq+li1bmupjXakNsBIvH3/8sbRu3VqaNm0qVatWNXfP2nCrbNmyZlet2tYl5nmw+Z9/\nEEAgZQHugFP2YSsCERFYuXKlua4Gx+nTp5u73CxZssiKFSvM+u+//97cFVevXl1uv/12sy6m0VTc\n3ag+o718+bK5i47d/uijj8onn3xi9s2QIYMULFjQvNbntVqdrMsPP/xg/o8NpubN//7ZsGGD6I8u\neoze8U6YMOF/W8W0jM6TJ49s2bIlbh0vEEAgeQECcPI2bEEgYgJaLazLnDlzZNGiReZ1mTJlzP/a\nCErvinXRhlOxS2xVsL7/448/zGptdBW7XH311bEvk/1f75S1G5EG0sRLnz59pFOnTqbrkW731t1I\n03PgwAET/DXIsyCAQPICVEEnb8MWBCImUK9ePXPtvn37mqpoveMdMmSILFy4UHLkyCHlypUz29es\nWWPujvXNmDFjTDX0ggULJDZY63F6J6xL7N2zeZPMP1o1rR0jtLo78aJBt2TJkpI3b16vwVf312fX\n2jKa4JtYj/cIJBUgACc1YQ0CERfQFsm69OvXT0aMGCExfXGlUaNGMnDgQLNeq56rVatm7nRr1Kgh\n2uDqpZdektWrV4t2B9IGXBUqVDADaOiAGb169Yp7fmxOkMw/N9xwg9miXZD8XfSu+9ixY3Lttdf6\neyj7IxCVAgTgqCx2Mm27gD63HTdunGhVtI4upXe+Dz74oAnGmnZ9BjtjxgzTL3f9+vVmvd69arel\nYsWKmTtQ3a5BWBtF6bm0lXTsyFXJ5V+DtS7ffPNNcrsku177DutSs2bNZPdhAwII/C3AQBx/W/AK\nASsFdJjIIkWKmK5G3hKo1b4XLlyQAgUKeNtsBtTQRlbaVSm1Rbs7aRVynTp1RFs9+7Nodbm2vNYW\n2MmlxZ/zsS8CbhcgALu9hMkfAn4KDB06VPr372/Gh9Yxn31Z9Lmx3oFr1bkGYRYEEEhdgCro1I3Y\nA4GoEujYsaPpzjRlyhSf860TOuzdu1d69uzp8zHsiEC0C3AHHO2fAPKPgBcBHdNZuxn5Um2th2tL\na/2J3xXKy2lZhQAC8QQIwPEweIkAAggggEC4BKiCDpc010EAAQQQQCCeAAE4HgYvEUAAAQQQCJcA\nAThc0lwHAQQQQACBeAIE4HgYvEQAAQQQQCBcAgTgcElzHQQQQAABBOIJEIDjYfASAQQQQACBcAkQ\ngMMlzXUQQAABBBCIJ0AAjofBSwQQQAABBMIlQAAOlzTXQQABBBBAIJ4AATgeBi8RQAABBBAIlwAB\nOFzSXAcBBBBAAIF4AgTgeBi8RAABBBBAIFwCBOBwSXMdBBBAAAEE4gkQgONh8BIBBBBAAIFwCRCA\nwyXNdRBAAAEEEIgnQACOh8FLBBBAAAEEwiVAAA6XNNdBAAEEEEAgngABOB4GLxFAAAEEEAiXQMZw\nXciG68ycOVMuXbpkQ1JIAwIIIICABQKFCxeWO+64IyIpSeeJWSJy5TBfdNasWTJs2DB55JFHwnxl\nLocAAggg4KvAL7/8IosWLZLGjRtL0aJFfT0s4P3eeust+fDDD6VKlSoBnyPQA6PmDljvfNu2bSsd\nOnQI1IrjEEAAAQRCKLB582bp27evvPfee1K5cuUQXunvU//0009y5cqVv1eE8VXUBOAwmnIpBBBA\nAIEABCpWrChz5swJ4EhnHkIjLGeWG6lGAAEEEHC4AAHY4QVI8hFAAAEnC1y+fFl+++03uXDhgpOz\nEVDarQvA+qz2xIkTAWWGgxBAAAEEnCOwf/9+adGihWzZskUyZ87snIQHKaVWBGD9y6dPnz5SqlQp\nUwj58+eXHDlyiD4PmDx5cpCyymkQQAABBGwROHTokPTs2dM0jK1Tp44tyQprOqxohNW1a1c5fPiw\nLFiwQMqUKWOC7x9//CFbt26Vbt26yV9//SWdOnUKKwwXQwABBBAIjYDWck6aNEl69eol//znP0Nz\nEQec1Yo74KVLl8r48eNNs/OcOXNKunTpJE+ePFKzZk0ZNWpUVLWKc8BnhiQigAACaRLIly+fvPzy\ny1EdfBXQigCsVc0rVqzwWqDz58+XQoUKed3GSgQQQAABBJwqYEUVdP/+/aVly5YyYsQIKVu2rOTO\nnVtOnTol27ZtM0NHLly40Km+pBsBBBBAIEbg4sWL8u2330qJEiWkdOnSmMQIWBGAq1atKhs2bJA1\na9bI3r17zfNgvevV5761a9c2VdKUFgIIIICAMwXOnDkjzzzzjNx2221Sq1YtZ2YiBKm2IgBrvrJm\nzSp169ZNksUdO3bIuXPnRIN0asvKlSvlm2++8brbqlWrzFifelcdu5QvX17Sp/+7Fl6vFX9IMrbj\nw+eD34/Y7wu+H/z/ftQeLtqY9oknnhD9Pj179qyMGTPGkOrr+Iv2fIm/hGv79u3b4182rK+tCcDJ\n5XrGjBmyb98+mThxYnK7xK3XgbsrVaoU9z7+i9WrV8uRI0ckV65c8VcneK0NwFKam4Lt+PD5SH7u\nFn4/+P2I//uhwVcb0movFp3o4KGHHkrwfZspU6YEtZuJB+II1/aNGzcmSFc430TNbEjdu3c3VdvT\npk0Lpy/XQgABBKJKQGss33nnHRk3bpwULFhQvvjiC8mSJYu1BtoXuVWrVnLTTTeFPY3W3wGHXYQL\nIoAAAgj4LaBVxvoIUGst9ZHf0KFD5d577/X7PNF0AAE4mkqbvCKAAAIhEtCeLDqXb4MGDUyvFh3P\ngSVlASsC8LBhw0wT9eSSet1118l9992X3GbWI4AAAghESEC7F919992mB4s2hNVGVvo8t1GjRhFK\nkXMua0UA1q5Ho0ePlkceecQMQ5mYj4E4EovwHgEEEIiswMmTJ001s47Xrw1gv/rqK/PsVxti6ShX\n8XsQRDal9l7digD89ttvm+4/2gUotom6vWSkDAEEEECgYcOGBmHu3LlSrVo1OXjwoLzwwgtROatR\noJ+Gvzv5BXqGIB03aNAg0QkYtMM2CwIIIICAvQI6cNLOnTtFu3dq8NVFR7iKxikF01JKVtwBawa0\nD+GHH36YlrxwLAIIIIBAGAQ6dOggTZo0oZo5jdbW3AGnMR8cjgACCCAQBgGdoe7nn3+W4cOHywcf\nfGDmctfBNlj8F7DmDtj/pHMEAggggEC4BLSLUceOHUVHjvrkk09k1qxZsnbtWhOIdShhFv8FuAP2\n34wjEEAAgagSOHDggJm7V8dr3rJlixlVUGc20i6k+viQJTAB7oADc+MoBBBAIGoEmjZtKrfccotM\nnz7d5FmHbtQflrQJEIDT5sfRCCCAgKsFFi9eLHv27DEtnl2d0QhkjiroCKBzSQQQQMAJAtq3t02b\nNtK1a1czWuHhw4fl8uXLTki6I9JIAHZEMZFIBBBAIPwCM2fOlGLFiplpBdu1a2ema82QIUP4E+LS\nKxKAXVqwZAsBBBBIi8BPP/0kI0aMkBtvvFEmTJhgXmswZgmeAM+Ag2fJmRBAAAFHC2h/Xq1uzps3\nr3z00Udy2223Sfny5aV9+/ZSuHBhR+fNxsQTgG0sFdKEAAIIhFng6NGjMnDgQNPY6r333pMnnnjC\nBN8wJyOqLkcAjqriJrMIIICAd4Hbb79dypUrZwbWqFGjhvedWBtUAZ4BB5WTkyGAAALOEtAx+HVC\nhXTp0sns2bPN6yVLlsiJEyeclREHppYA7MBCI8kIIIBAMAR0hKsuXbqY+XvXr18v2vDq8ccfN/P7\n5suXLxiX4BwpCBCAU8BhEwIIIOBGAY/HI/PmzZO77rrLNLTSmY327dsnffv2lZdeesm0fHZjvm3L\nE8+AbSsR0oMAAgiEWKBRo0aydetW6d27t+jUgidPnpTt27fL2LFjTb/fEF+e0/9PgADMRwEBBBCI\nEoFNmzbJ008/bYaW1ICbPXt2k3PtdnTfffdFiYI92aQK2p6yICUIIIBAyAQuXrwoDz74oFx//fWy\nc+fOuOAbsgty4lQFCMCpErEDAggg4HyBl19+WXSgDa1mzpIli2nlrCNd/fjjj87PnENzQAB2aMGR\nbAQQQMBXgWXLlokOrqFdjtKnTy9HjhyRTp06ybXXXiuVKlXy9TTsF2QBngEHGZTTIYAAArYIaLWz\nBlidwejNN9+UWrVqybFjx6R79+7y0EMPyT333GNLUqMyHQTgqCx2Mo0AAm4X0KkEO3fuLNmyZZO1\na9dKpkyZTCA+ffq0jBkzRujnG/lPAFXQkS8DUoAAAggEVeDbb78V7Wp0/vx50SkFNfjqolMJli5d\nmuAbVO3AT0YADtyOIxFAAAGrBC5duiRTpkyRe++9V2rXri1z586VsmXLWpVGEvO3AAH4bwteIYAA\nAo4WGDx4sBnJavTo0TJq1CjT2vnPP/+Unj17ytSpUx2dNzcmnmfAbixV8oQAAlElsGLFCnO3O2PG\nDNPNSO+AddHg++yzz5q+v23atIkqEydkljtgJ5QSaUQAAQRSEGjXrp1pbKXPe2OD74ULF6Rfv35S\noEABM+FCCoezKUIC3AFHCJ7LIoAAAsEQ0KEltbvRgAEDEpwuc+bMolXSLPYKEIDtLRtShgACCCQr\noNXL06ZNM4NrfPfdd8nuxwZ7BQjA9pYNKUMAAQSSCFy5ckXatm0rq1evlgoVKsj48eOlTJkycfsd\nPXrUvC5cuHDcOl7YKUAAtrNcSBUCCCCQREBHtLrttttEg+yaNWukSJEicftoYH799dfNeM+vvfZa\n3Hpe2CtAALa3bEgZAgggkECgf//+cvjwYTN3rz7jjb8MHDhQzpw5I2+88YYZcCP+Nl7bKUAAtrNc\nSBUCCCAQJ6B3tzqb0bvvvmsmVUgcfOfMmSP16tWTm266STJm5Gs9Ds7yF5SU5QVE8hBAAIF58+bJ\n5MmTZenSpV5nL7rvvvtAcqAA/YAdWGgkGQEEoktgyJAhZmhJpg50V7lzB+yu8iQ3CCDgIgFtdKXP\nfXfv3i16Fxx/2bp1q3kerI2yYidbiL+d1/YLcAdsfxmRQgQQiEIBnVRBJ1L4/PPPRYeY1BGtYhcd\n11nviqtVq0bwjUVx4P/cATuw0EgyAgi4V2D48OFxEye8+eab0rJlywSZnT59unz22WdmsoXcuXMn\n2MYbZwkQgJ1VXqQWAQRcLPDHH3+YwNq7d2958sknk9zd/vzzz5InTx6ZMGGCZM2a1cUS0ZE1AnB0\nlDO5RAABywU2bNggDRo0kFtvvVU6deok6dKlS5JirZJmft8kLI5dwTNgxxYdCUcAAbcIXLp0STp3\n7iy1atUy0wp6C75uySv5+FuAAPy3Ba8QQACBiAjoABvHjh2TSZMmJbn+5s2bzaxGv//+e5JtrHC2\nAAHY2eVH6hFAwAUC77//vjRs2DBBS2fN1qpVq8wIWC1atEiyzQXZjvos8Aw46j8CACCAQKQEtmzZ\nIuvWrZOdO3earkbx0/HNN9+ItogeOnSoXH311fE38dolAtwBu6QgyQYCCDhLQFs833XXXfLDDz/I\nmDFjpGTJknEZ0Ll+ixUrJtrliEZXcSyue8EdsOuKlAwhgIDNAjqb0aeffiovvfSSVK9e3dzhJm50\nlS1bNu56bS7EIKWNABwkSE6DAAIIpCSgd7X/+c9/ZNy4cZI/f37Tl7dJkyYpHcI2lwtQBe3yAiZ7\nCCAQeYHVq1ebYSN1CMlnnnlGlixZIomDrz4H1lGvNm7cGPkEk4KwCFh7B6zzX547d05y5swZFggu\nggACCIRKoEuXLvLPf/5TJk6c6HW+3l27domOfvXcc89JlSpVQpUMzmuZgBV3wNoYQVv6NW7cWFas\nWGFm/ShSpIgUL15c2rdvL2fOnLGMjeQggAACqQvonW/Tpk3l+PHjMmzYMK/B98CBA/LGG2+YAHzL\nLbekflL2cI2AFQF44MCBsnbtWrn77rtN9czzzz9vgvDevXvlwoULSZrnu0afjCCAgKsF9AYiX758\n8t1335nnvt4yW6pUKZk8ebLUqFHD22bWuVjAiirouXPnmgCcI0cOOXLkiBkRpmbNmob9hRdekB49\neki7du1cXAxkDQEE3CYwaNAg0dGrtNFVlixZ3JY98hMEASvugK+//nozvdapU6dk5cqVsn79+ris\naR+5m266Ke49LxBAAAHbBebPn2/m6503b57X4KtjP+/Zs0f0O48legWsCMA9e/YUvdMtXLiw3Hzz\nzXLjjTeKBuWHHnpIunfvLo888kj0lhA5RwABRwl8/vnn8uijj0q/fv1MP9/Eidcxn/U7TVs769SC\nLNErYEUVtFY3b9261TRUKFCggJw/f9400z958qR5NqKd0lkQQAAB2wXWrFljHpc1a9ZMunbtmiS5\nJ06cMI/UtAuSNs5iiW4BKwKwFoGOBKPBVxd9XhLbR27Hjh2mO1LVqlXNtpT+0dbUp0+f9rqLtqTW\nah8WBBBAIBQC2rdXG5DecccdMnbs2CSX0IE4pk2bZnp26LSDLAhYE4CTK4oZM2bIvn37TP+55PaJ\nXb9o0SIzxFvs+/j/64DnOrYqCwIIIBBsAa1Wbt26tRleUqufM2TIkOQSWpOnc/6yIBArkM4Ts8S+\ncfP/+ixZx2DVv0BZEEAAgWAK6CAbOm7BnDlzgnlazhUGAW2D1KpVq4g09rWiEVZ8Y60m1uckLAgg\ngIATBB5//HE5evSovPPOO0mSe/nyZVm+fLls2LAhyTZWIGBFANbBNvr06SPaIT1z5symw7r2Ca5Y\nsaJphEUxIYAAArYJaHCtXbu2LF68WPTxV9GiRRMkUb/XdAyDTZs2iS9tWBIczJuoELAiAGtrQZ2Y\nesGCBaINqXQc6EOHDpnnvvpXpXZkZ0EAAQRsERg+fLhce+215lmv9uDQbpPxl4sXL5quldruRIMw\nCwLeBKxohLV06VLR5vvx/4LU/nHaPWnUqFGmP12nTp28pZ91CCCAQNgF3nrrLRNgO3TokOTa2qzm\n+++/N2M769gGLAgkJ2DFHbBWNeskDN4WHVGmUKFC3jaxDgEEEAi7QKNGjUSD7GOPPeb12tqlUsd1\nJvh65WFlPAEr7oD79+9v5sEcMWKElC1bVnLnzm2GaNu2bZvpu7tw4cJ4SeYlAgggEH4BndFIv6u0\nS6N2jcyUKVP4E8EVXSVgRQDWBgraSlCroffGzICk3YX0rlernbWRg/5FyYIAAghEQkC/j3788Udz\nx6vTBWqrZm+TK0yfPt1MnZpcP+BIpJ1r2i1gRQBWoqxZs0rdunXt1iJ1CCAQdQJ6E1C5cmV5+umn\npVevXl7zr/OZHzx4UHQGJG+DcHg9iJVRL2BNAI76kgAAAQSsEtA73dmzZ4uOSf/xxx9Lxozevy61\noahWSQ8ZMsR0o7QqEyTGagHvnyirk0ziEEAAgdAL6B1vnTp1zBSpyQVfnU5Qu1HqYzIelYW+TNx2\nBQKw20qU/CCAQJoEdCyC++67T3TyhGHDhpnHY8mdkOkEk5NhvS8CBGBflNgHAQRcL6Bdi3SKwF9+\n+UXy5s1r+vJq2xQWBEIlYEU/4FBljvMigAACvghs3rzZdIXU/7Xbo/5oEPa2zJw507SI1mfDLAik\nRYAAnBY9jkUAAccL6NC3DRo0MF0f582bZwbQ0DHpvS3//e9/RX8GDx6cbID2dhzrEPAmQBW0NxXW\nIYCA6wU08NavX1/OnTsnBQoUEB1eMqVFW0XrsLm6n+7PgkBaBQjAaRXkeAQQcKSATvSya9cu2blz\np0/dh+68807RHxYEgiVAAA6WJOdBAAHHCGhVsw4rOXDgQJ+Cr2MyRkIdJcAzYEcVF4lFAIG0Crz6\n6qvy5JNPmlGrdNjIlJbTp0/Ljh07ROf2ZUEg2AIE4GCLcj4EELBWQO96J02aJB988IE88sgjKabz\nm2++Mfvo0JLJNcpK8QRsRCAVAaqgUwFiMwIIuENAh4vUGde0IVX16tVTzNT69etlwIABprVzuXLl\nUtyXjQgEKkAADlSO4xBAwFECOpGCzryWWvDdv3+/rFq1SsaOHSslSpRwVB5JrLMECMDOKi9SiwAC\nAQiMHDlSVqxYIYsXL0716Kuuukq6deuW6n7sgEBaBQjAaRXkeAQQsFpAG12NGzdOtNvRzTffbHVa\nSVx0CdAIK7rKm9wiEFUC2oJZpwucNm2aNGvWLNm867CSOsSkzunLgkC4BAjA4ZLmOgggEHYBHTij\ndevWUrdu3WSvrY2zHnvsMSlatCjPfJNVYkMoBAjAoVDlnAggEFEBHS5SWy8XLlw4xSEmdeajHj16\nmDl9a9WqFdE0c/HoE+AZcPSVOTlGwLUCOpevBtRFixbJ+PHjpVGjRsnmVef71WEo3377bSlevHiy\n+7EBgVAJEIBDJct5EUAg7AI6tOTq1atNa+dKlSqleP1s2bKlWDWd4sFsRCAIAgTgICByCgQQiLzA\n2rVrTWvn6dOnS2rBN/KpJQUIiPAMmE8BAgg4XuDnn3+W5s2bm3l9dW7f5JYzZ86YCRi++OKL5HZh\nPQJhEyAAh42aCyGAQCgE9u7dK/fee6+ZKvDjjz9O9hKnTp0yja10dKs6deokux8bEAiXAAE4XNJc\nBwEEgirg8Xikb9++5jlumTJlZOjQoZIuXTqv19A73549e8ptt90mbdq08boPKxEItwDPgMMtzvUQ\nQCDNAtddd50UKVJEjh8/bp77NmzYMMVznj9/XiZMmCDp03PPkSIUG8MqQAAOKzcXQwCBtArs3r1b\ntLvR9u3bfT5VgQIFfN6XHREIlwB/DoZLmusggECaBfROVoeUrFGjRprPxQkQiLQAATjSJcD1EUDA\nJwEdr1lnKdLnudrVKKXl0qVLotMPar9gFgRsFSAA21oypAsBBOIEVq5cKbVr15b169fL1KlTJXPm\nzHHbEr+4fPmyvPTSS5I7d2557rnnEm/mPQLWCPAM2JqiICEIIOBNYOvWrfL0009L9erVZcyYMZI1\na1Zvu5l1V65cMbMf6TjQ7dq1o9FVslJssEGAAGxDKZAGBBDwKqBjOf/444/SuHFjefPNN1MMvnoC\nbeWsY0GzIOAEAQKwE0qJNCIQZQLa0lm7FukdrbZ21nGbWRBwmwAB2G0lSn4QcLCADq6hg2WcPn1a\nKleuLDNnzvQpNwcOHBBtIV22bNlkB+Pw6UTshEAYBQjAYcTmUgggkLLA4MGD5dChQ6LPfVN61hv/\nLDr3r04rqMcmNxJW/P15jYAtAgRgW0qCdCAQxQJHjx6VJ554Qr777jszP6+vwXfcuHGybds2GTZs\nGNXUUfz5cWrWCcBOLTnSjYCLBDp16mSqnTdt2iSFCxf2KWdfffWVZM+eXUaOHClZsmTx6Rh2QsAm\nAQKwTaVBWhCIMgF95tu6dWtZtWqVLFiwwOfgq0y1atUyP1FGRnZdJEAAdlFhkhUEnCbw6aefyrJl\ny2TLli1SqFAhpyWf9CKQJgFGwkoTHwcjgEAgAjpa1ezZs+WVV14x3Y18Db7aJUlbRv/555+BXJZj\nELBKgDtgq4qDxCDgfoF9+/bJvffeKxcvXpSuXbtKq1atfMq0VlFPmTLFNNKiX7BPZOxkuQAB2PIC\nInkIuE1Ah5PUALpx40afs/bZZ5/J+++/b4aZ9LWRls8nZ0cEIiRAFXSE4LksAtEooJMpzJgxQ+65\n5x6fs3/kyBHTJ/iDDz6Q4sWL+3wcOyJguwB3wLaXEOlDwEUCjz32mDRo0ED69Onjc66KFCki+sOC\ngNsEuAN2W4mSHwQsFRgxYoTo3eyoUaOYpcjSMiJZ4RUgAIfXm6shELUCEydONNMK+jLK1Y4dO8xd\n8v79+6PWi4y7X4AA7P4yJocIRFzg+PHjcvjwYROAU0uMjobVu3dv0zr6qquuSm13tiPgWAFrA/Bf\nf/1luik4VpaEI4CAEdBWzzVq1JAbbrhBcubMmaKKTsLQr18/0z+4QoUKKe7LRgScLmBFANZqprZt\n28q6devkt99+k8cff1yKFi0qefPmFW20ceHCBac7k34Eok5An/dqP9833njD9N/VsZtTWnRwjmLF\niplW0lWqVElpV7Yh4AoBKwLwyy+/LFrVpH/xvv3223Lp0iXZvHmz/PDDD2aA9tdee80V2GQCgWgQ\n0LvY5s2byz//+U/ReXp15KqaNWummvUMGTJIvnz5JFOmTKnuyw4IuEHAim5IK1euFB1iLnPmzPLf\n//5X5syZIyVLljS+Gnw7duzoBmvygICrBc6cOSM6YEbnzp2lTZs2ZpaiEiVK+JRnvfvVAMyCQDQJ\nWHEHfO2115oqKoWvU6eOLFy4MK4M5s+fL9dcc03ce14ggIB9AufPn5e77rpLtKvRSy+9JIMHDxZf\ngu8vv/wiDzzwgHz55Zf2ZYoUIRBiASvugLWRRqNGjeTdd9+VcuXKybPPPivvvfee6Sv4xx9/iN4h\nsyCAgJ0C586dk3r16ok2nNQ/mH0drerXX3+VHj16yJNPPil33HGHnZkjVQiEUMCKAFy2bFnR50Za\nfaX9//R5sD4L0jtfHbIuY0YrkhnCYuDUCDhPQIPtsGHDRO9+td2GNrLKnj27TxnRbklDhgyRp556\nSm6//XafjmEnBNwmYE1kS5cunRmiToepY0EAAXsFdu7cKQ8//LDoHeygQYNM1bN2L8qSJYvPic6f\nP78MHz7c5/3ZEQE3ClgTgJPD1TtireKqWrVqcrvErZ8wYYJ89NFHce/jv9i1a5eULl06/ipeI4CA\nnwJ6t6stnLXWSu94fRnVys9LsDsCUSNgfQDWmVN0/lAdxi61pX379qI/3pbu3bubkXi8bWMdAgik\nLqA9FRo2bGh6KOjUgP4GX23pvG3bNtG7X1+fE6eeKvZAwLkCVrSCjs+nz5JOnDgRt6pv374+Bd+4\nA3iBAAJBF/jpp59Me4zatWubO99cuXL5dQ3totShQwf59ttvCb5+ybGzmwWsCMA60pVOT1aqVCnT\nF1j/Qs6RI4dUrFhRJk+e7GZ/8oaA9QJnz56Ve++91wysMWXKFL/T++eff0rPnj2lWrVqZpQ7v0/A\nAQi4VMCKKmgdrk4Hal+wYIGUKVPGBF/tfqQto7t162a6N3Tq1MmlRUC2ELBXQH8vBwwYYHoiTJo0\nye+Eao3WvHnzzKOhm2++2e/jOQABNwtYcQe8dOlSGT9+vFSuXNkM1q4tovPkyWOGr9O5Q3VkLBYE\nEAivgNZMaa8ErX4eOXKk+cPY3xRoF8IWLVoIwddfOfaPBgErArBWNa9YscKrt/Y1LFSokNdtrEQA\ngdAJxPbx1eFh77zzztBdiDMjEKUCVlRB9+/fX1q2bGmGsdPuDblz55ZTp06ZFpNahRV/aMooLSey\njUBYBa5cuWL+KNZGV/62dvZ4PLJkyRLTnoMRrsJabFzMYQJWBGDt47thwwZZs2aN7N271zwP1rte\nfe6rXwBaJc2CAALhE/j8889l06ZN8vrrr/t1Ue1q9Oqrr5pjdF5fFgQQSF7AigCsydO/suvWrZt8\nStmCAAJhEdBnvzrCVa1atUzLZ18vqne+OvfvxYsXRWcxY3YjX+XYL1oFrAnA0VoA5BsBWwR0xLnl\ny5dLr169zLCSOrKcP4s21tLJFbQLIbVW/sixb7QKEICjteTJNwKJBB599FHRcZ61W6BOkuDvUr58\neX8PYX8EolqAABzVxU/mEfh/gQMHDpi737Vr15pxnnFBAIHQC1jRDSn02eQKCCCQkkCbNm1Mv3vt\nheDPMnv2bBk4cKAZLMef49gXAQRECMB8ChCIUgHtarRs2TIzBKw2vBo3bpxfEjpMrHYR7NKli99d\nlfy6EDsj4FIBqqBdWrBkC4GUBHRSBL3rTZ8+vRmH3d+hXnXazy+++ELeeustM3pdStdiGwIIeBcg\nAHt3YS0CrhZo1aqVmWBh8ODBfncX0rvlJk2amMFzXI1E5hAIsQABOMTAnB4BmwR0ZDmtatYZioYO\nHRpQd6HMmTObUa5syhdpQcCJAgRgJ5YaaUYgAAENuo8//riZk3fIkCF+B18d5UoH2fB3aMoAksoh\nCESFQLIBWKuZdHjILVu2yO7du6VcuXJmRpPrr7/eTE0WFTpkEgEXCbz00kuyefNmWbx4sVxzzTV+\n5eyzzz4THZhDZycrXry4X8eyMwIIeBdI0gpah5ObNm2aVKhQQbp37y46Jqz+5aytHfW5UcmSJUVn\nSdEAzYIAAs4RmDVrlhlkw9/gq42tdLpQ/b0n+DqnvEmp/QIJ7oD/+usvadasmdSvX99MjFCwYMEk\nOThx4oSMHTtWdJaTTz75hF/IJEKsQMAuge3bt8vTTz9tEvXEE0/4lbiNGzfKvHnzTGtngq9fdOyM\nQKoCCQKwTp798ccfS65cuZI9MF++fPLiiy+aMV+1HyELAgjYK/Dbb79Jw4YN5e6775aZM2f6/dy3\nSpUqoj8sCCAQfIEEVdAagGOD76effirHjh1LcMWff/7ZVEPpymzZsplB1xPswBsEELBK4Mknn5Qi\nRYrI6NGjzTzbViWOxCAQ5QIJAnB8iyNHjojO06vPf3R59913zdRksQHarOQfBBCwUkBrp9q2bSvr\n1q2TKVOmmAE3fE3o6dOn5bvvvpMzZ874egj7IYBAAAIJqqDjH6/Piq6++mrzS1ysWDHRLggajCtV\nqhR/N14jgIBlAnv37pVXXnnFTK6ggdSfZ7dbt241I2Pp8Tlz5rQsZyQHAXcJJHsHHJvNTJkymYHW\ndX5PHbaOBQEE7BVYvXq11KpVyzw+WrlypV/Bd8eOHaZ9R9++fXnua28RkzIXCSQbUceMGWPGitV+\nf5s2bZKOHTvK7bffLpMmTXJR9skKAu4R0EkRHnjgAenVq5fMnz/fr2kFf//9d/n6669l5MiRUq1a\nNfegkBMELBZItgr6qquukh9++EEKFy5skq8j6NSpU0fmzJljcXZIGgLRKbB+/XqZMWOGLF26VG68\n8Ua/EQoUKCCPPfaY38dxAAIIBC6Q4A5YB9xYsWKFOVvjxo3jgm/s6XWu0J49e5q3q1atkl9//TV2\nE/8jgECEBHROXp0coXPnzgEF3wglm8siEPUCCQKwDrKuDa2aNm1qRsM6depUAqBffvlFdDSd2rVr\nm+20iE7AwxsEwi7w7LPPmtGtBgwYIDrUpD+Ltnb+4IMPRAfqYEEAgfALJKiCzpAhg7z66qty6NAh\n6devn/nF1u4M2o9Qg2/evHmlZs2aos+HaQ0d/sLiigjEFzh69Kj5Q1iHjNUuR/4sOkBH165dpXXr\n1nLdddf5cyj7IoBAkAQSBODYc2q3hYkTJ5qfw4cPx03GEPs8OHY//kcAgcgInD9/3kyMoM9uNQD7\ns2iDq27dusmDDz4ojRo18udQ9kUAgSAKeA3A8c9ftGhR0R8WBBCIvID269XZjL766itTU9W7d2/R\nmitfF50PeP/+/aLTEfrTP9jX87MfAgj4LpBsAD558qRp1PHjjz8mmPnorrvuMl0VfL8EeyKAQDAE\ntNpY71h1kBydXEEfB+XPn9+vU+twszrCHQsCCEReINkAPHjwYNFGWG+99VaCEXH8/YWPfBZJAQLO\nF1izZo20aNFCbr31VnnjjTecnyFygAACkmwAPnjwoLkDrlu3LkwIIBAhgT179phntTouc8uWLWXg\nwIF+p0SfF2uVc/ny5aV58+Z+H88BCCAQGoEE3ZDiX+L++++XqVOnira0ZEEAgfAL6EQK9erVM4Fz\nw4YNAQVfneNbuyppl0GCb/jLkCsikJJAsgFYuyItXLhQdCKGa665xnRV0O4K2nqSBQEEQiugvRD0\nee8jjzwi7733nmTNmtXvC164cEFeeOEFM6nKM8884/fxHIAAAqEVSLYKWn/5vY0JyzPg0BYIZ0dA\npwLVfvivvfaa6Hy+aVlGjBiRlsM5FgEEQiiQJAA3bNhQPvroIylVqpQULFjQzKqir1kQQCD0Au+/\n/75oA8ibb77ZtHZOyxV1ZDsWBBCwVyBJFbQ+d7p48aJJ8dq1a03DD3uTT8oQcIeAzrf95ptvmrHW\nX375ZZk3b57oFKCBLDqmu/6wIICA3QJJ7oDtTi6pQ8B9AufOnZMKFSpInjx5ZO7cuWY+30By6fF4\n5PXXXxc9n1ZfsyCAgN0CBGC7y4fURYHA8OHDTSvljRs3pim3gwYNMn339U5aB9xgQQABuwW8/pbq\nxAvafUHHgdY+hPv27YvLRY4cOcyz4bgVvEAAgYAF9HHPlClTpFWrVgGfQw/88MMPpUqVKnLnnXdK\npkyZ0nQuDkYAgfAIeA3AiVs//+Mf/4hLjfYl/OSTT+Le8wIBBAIXmDBhghnqVfvqpmVJawBPy7U5\nFgEEAhNIEoC1C0RKS6ANQ1I6J9sQiEYB7Wevc/jqnLxas8SCAALRJZAkAPszs0p0UZFbBIInoCNb\nde7c2Qy00bhx44BOrI+Gjh07Zqqe+b0NiJCDEIioQJJuSBFNDRdHIAoEtm3bJnfffbc0bdrUjNEc\nSJa1r752VypXrpxf0xEGci2OQQCB0AgkuQMOzWU4KwIIqMDOnTtN8G3SpIkEOkrVjBkzZNGiRTJq\n1CjTdQlZBBBwpgB3wM4sN1LtQIGePXvK7bffbsZ4fueddwLKwdatW0UH7Rg/frzfcwEHdEEOQgCB\nkAlwBxwyWk6MwN8C69evN7OLLV++XCpVqvT3Bj9f3XDDDaI/LAgg4HwB7oCdX4bkwAEC+rz21ltv\nTVPwdUA2SSICCPghQAD2A4tdEfBX4I8//jDz8G7evNmM9ezv8bq/PjfW/sInTpwI5HCOQQABSwUI\nwJYWDMlyh8D9998vu3btkjVr1gRUdfz1119L7969pX79+pIvXz53oJALBBAwAjwD5oOAQAgELl26\nJN9++63o+M67d++W3Llz+30VnY1s6NChpqtS6dKl/T6eAxBAwG4BArDd5UPqHCigDa6eeuopMyVg\nly5dAgq+p0+fFp3PV8d4zp49uwMVSDICCKQmYG0A1skgLly4ENCXV2qZZjsCoRL47bffpF27dmZ0\nKu1qFGjwzJUrlzlHqNLJeRFAIPIC1j4DnjVrlvTo0SPyQqQAAR8FdBz1tm3bmme1I0eODDj4+ng5\ndkMAAYcLWBGAr7nmGvOlpY1MYn/at29v+k3qe72jYEHAZgF9XluzZk3RyUomTpwY0CAZe/fula5d\nu4q2mGZBAAH3C1gRgCdPniyFChWS7t27m0Yr2nBFJxXXsXL19ZAhQ9xfEuTQsQI6IYIOLdmyZUv5\n9NNP5dprr/U7L9pQS2t89BwVK1b0+3gOQAAB5wlYEYBr1aol69atM9019EtIp2YrWLCg5MyZU66+\n+mrz2nm0pDgaBPSOtU6dOlKjRg15/fXXA5oYYf/+/aLzAXfr1s3cRUeDG3lEAAERaxphaTeNKVOm\nyCeffCK1a9c2X2hMscZH1GYBHRJSuxvNnz/fzEoUaFpLlCgh2uaBubYDFeQ4BJwpYMUdcHy6Bx98\nUJYuXWrmOS1atGj8TbxGwBoBnZHo+PHjsn37dlPlnD594L9K+ocmwdeaoiUhCIRNwJo74Pg5Llmy\npHmWput27Ngh586dk6pVq8bfxevr7777zjwz9rbxxx9/lCxZsnjbxDoE/BLQ57UdOnSQcePGSaCB\n1+PxyKlTpyRbtmx8Lv3SZ2cE3CMQ+J/tYTLQO42xY8f6dLVMmTKZ58f6DDnxj27jLsMnRnZKQUAD\npz73rVChgrRo0SKFPZPfpHfOrVu3Fp0ZiT8Kk3diCwJuF7DyDjg+et++feO/TfF1lSpVkh28QO+O\nDx8+nOLxbEQgJQFt7azz8P7www+iYzQHspw8edI0tmrcuLE0a9YskFNwDAIIuETAugCsjVp0GD4G\nnnfJJ8wl2Th69KjUrVtXtF2Cdpu76qqr/M6ZPkrR0bG0q1HDhg39Pp4DEEDAXQJWVEHrkJN9+vSR\nUqVKmfFv8+fPb6qQtT+kftmxIBBpAZ3LV6ud582bJ/Xq1QsoOTos5fPPP0/wDUiPgxBwn4AVd8D6\nTE2rhxcsWCBlypQxwVfnUd26dauprtNxoTt16uQ+fXLkCAGdSlAbTGkXORYEEEAgWAJW3AFrtyN9\ntla5cmUz+IY2lsqTJ48ZlGDUqFEyZ86cYOWX8yDgl8DOnTvNKFe9evXy67jYna9cuWIGmdHW/CwI\nIIBAfAErArBWNa9YsSJ+uuJe6yAHOkwlCwKRENDhUO+77z557rnn/L68PlrRkd2++OILKV++vN/H\ncwACCLhbwIoq6P79+5uGKSNGjJCyZcuaKQi1ym/btm1mpKGFCxe6uxTInZUC2tL50KFDpr+vvwm8\nePGiaAv+IkWKSM+ePf09nP0RQCAKBKwIwDrIxoYNG0SftemMMPo8WO969bmvDktJ/90o+CRalkX9\no3DMmDGmcWDGjP79mmhf4ZUrV5rPb+nSpS3LGclBAAFbBPz7ZglhqrNmzWq6eYTwEpwaAZ8FJk2a\nJG+99VZAg23oH4x33nmnz9diRwQQiE4BawJwdPKTa5sEtOW9DjOpQ0zq8sADD9iUPNKCAAIuE7Ci\nEZbLTMmOAwV0dKu2bdtK586dRauQteGUv7NxLVq0SHTo1MuXLztQgCQjgEC4BQjA4RbnetYJaOM/\nHZlKB4J5//33ZcKECaY/uj8JHT58uCxZssR0WfI3cPtzHfZFAAH3CFAF7Z6yJCcBCGj/3IEDB4oG\n0FatWgVwBpHRo0fLzz//LEOHDmVyhYAEOQiB6BQgAEdnuZPrGAHtZtS8eXNp06aN6QYXCMqvv/4q\n7du3F20pHejUhIFcl2MQQMD5AlRBO78MyYGfAtpHV+9WdTYi7eqmrwPt6lasWDEzfjnB189CYHcE\nEBDugPkQRI3Arl275IknnjCzbZ05c8YMcXrLLbdETf7JKAII2CXAHbBd5UFqQiSwceNGadCggWlc\npXND67PfQIOvtnbWISZ1tDYWBBBAIFABAnCgchznGIE333xTdExnHRxDB9hIS3WxTkc4ZcoUeeGF\nF8yEIY5BIKEIIGCdAFXQ1hUJCQqWgD7jPXDggBw7dsxMJVi9evU0nXrx4sXy0Ucfic7QxQQhaaLk\nYAQQiBEgAPMxcJ2AVg2//fbb8tVXX5kxxosXLx6UPGpfYf1hQQABBIIhQAAOhiLnsEZA+/NOnTrV\nzKild6vBCr7WZJCEIICAawQIwK4pSjKisxfp814dyeqee+4JyqAY58+fN42tChQo4PfQlJQIAggg\nkJIAjbBS0mGbIwTOnj1rpg585ZVXZNq0aXL//fcHJfiuXbtWHnroIROAGV7SER8FEomAowQIwI4q\nLhKbWECD7nXXXSfTp083d77169dPvEtA73V+6jfeeEMGDBgg11xzTUDn4CAEEEAgJQGqoFPSYZv1\nAu+++66MHTtWGjduHLS07tmzxwzSMXLkSCldunTQzsuJEEAAgfgCBOD4Grx2lICO5XzhwoWgBl8F\n0KD76quvOsqCxCKAgPMEqIJ2XpmR4v8JtG7dWlq0aIEHAggg4EgBArAji41Ez507V7SF8ltvvRUU\nDB0bevny5XL06NGgnI+TIIAAAqkJEIBTE2K7dQLHjx83sxj16dMnKGnbv3+/tG3bVjJlyiSFCxcO\nyjk5CQIIIJCaAAE4NSG2WyWg8+9qH99y5cpJ586d05y2gwcPSs+ePeWpp56S2rVrp/l8nAABBBDw\nVYBGWL5KsZ8VAn379pUrV66IjsuclkkVNDPaf1j7+g4aNMjMkmRFBkkEAghEjQABOGqK2vkZ/euv\nv2ThwoVmNqLs2bOnOUM5cuQwsySl+UScAAEEEAhAgCroANA4JDICWuWsz2jr1asXmQRwVQQQQCCI\nAgTgIGJyqtAJaHWx3v2+9tprki5duoAvdO7cORk/frx88803AZ+DAxFAAIFgCBCAg6HIOUIqsH37\ndqlYsaIZErJJkyYBX+v06dPSpUsXyZo1q9xyyy0Bn4cDEUAAgWAI8Aw4GIqcI2QCOrHCTz/9JJUr\nVxbt+xvoonfQPXr0kFtvvVUeeeSRQE/DcQgggEDQBAjAQaPkRMEWWLZsmaxatUp27Ngh+fPnD/j0\nHo/HDLAxevTooMySFHBCOBABBBCIJ0AAjofBS3sEHnzwQdNFqF+/fmkKvpojfWbMpAr2lC0pQQCB\n/xcgAPNJsEpA+/jqQBubNm2S1atXyz/+8Q+r0kdiEEAAgWAJEICDJcl5giKwbt062bhxo+zcuVNy\n5swZ8DkvX75s5vLNkyePdO3aNeDzcCACCCAQKgFaQYdKlvP6LfDHH3+ItnJu2bJlmoPvyy+/LBqE\ndYhJFgQQQMBGAe6AbSyVKE2TDjNZtGhRGTZsWMAC2uDqzTfflFy5cslzzz2X5uEqA04IByKAAAKp\nCBCAUwFic3gE1q9fL9OmTUtT8NWUaoOrl156KTyJ5ioIIIBAGgSogk4DHocGR2Dz5s3SoEEDefbZ\nZ820gME5K2dBAAEE7BYgANtdPlGRulGjRkmFChWkd+/eAef3xIkTpq9vwCfgQAQQQCDMAlRBhxmc\nyyUU6NChg8yZM0dmzZqVcIMf73SAjW3btsnQoUP9OIpdEUAAgcgKEIAj6x+1V9dJEdq3by8rVqwQ\n7Xp09dVXB2QxYcIE+fHHH2X48OGSLVu2gM7BQQgggEAkBAjAkVCP8mseOnRIqlataoKuPv/Nly9f\nQCI6VKUG8pEjRxJ8AxLkIAQQiKQAz4AjqR+F19aWztWqVTP9fdeuXRtw8FU6nRe4W7duBN8o/ByR\nZQTcIMAdsBtK0SF52L9/vzz99NMyZMgQefTRRx2SapKJAAIIhEaAO+DQuHLWRAI6HaAG3UqVKqUp\n+O7bt0+06vmvv/5KdAXeIoAAAs4SIAA7q7wcm1pt7Xzs2DGZPn16wHlYtGiRGd2qYsWKkjVr1oDP\nw4EIIICADQJUQdtQCi5PQ82aNUXHeZ49e7YUKlQooNzqXe97771nGlzpcJUsCCCAgNMFCMBOL0EH\npP+nn36SI0eOSMaMgX3cfvnlF9PaWQOwjvHMggACCLhBILBvRDfknDyEVEBnIvrzzz+lRYsWUqJE\niYCDryayZMmS5iekCebkCCCAQJgFrH8GrF/k58+fDzMLl0urgM5sVL16ddHZiT799NO0no7jEUAA\nAdcJWBGADxw4YAbh1wnY69evL7t27YqDnjFjhrRp0ybuPS/sF3jnnXfM89ouXbrIwoULAxrlas+e\nPaa70sGDB+3PMClEAAEEAhCwIgCPGDFCihUrZoYk1AY7tWvXFn1uyOI8AW3p3L9/fxM8n3rqqYAy\noENL6gAbDRs2NNXXAZ2EgxBAAAHLBax4Bqx3SRs2bDAjGumX9w033CD//ve/5auvvrKcj+TFF1i5\ncqUJnJUrVw54WkGdVEHn833llVdMn+H45+c1Aggg4CYBK+6ANeDqgPyxy0MPPSRdu3aVu+66S37/\n/ffY1fxvsYBOB6h3vDVq1BB9bBDIcuHCBcmSJYt89NFHZqzoQM7BMQgggIBTBKwIwB07dpTmzZvL\noEGD4tx69OghzZo1k+7du8et44W9AvqcXrsIDR48OOCuQpkzZ5YyZcpI9uzZ7c0oKUMAAQSCJGBF\nFXSDBg3k559/lt27dyfIVr9+/eT222832xJs4I1VAvr4YP369abhFf10rSoaEoMAAhYLWHEHrD45\ncuTw+sxPG2fddNNNPhFevHjRDNigU9Ql/tFtV65c8ek87OS7gFYbv/jii+a5vTaa8nc5fPiwtGvX\nTr788kt/D2V/BBBAwNECVtwBpySozxN1AP6JEyemtJvZ9uGHHyb7/HHr1q1y1VVXpXoOdvBd4NSp\nU3L33XeL/nHz8ccfS7p06Xw/OGZPDb7PPPOMmZxBazpYEEAAgWgSsD4A64AOvi46247+eFv0WbJ+\n4bMET2DcuHFmgoWvv/5aChYs6NeJtXFdnz59THlpYzsWBBBAINoErAvAly5dktOnT6dpovZoK8RI\n5FcD6MiRI+WFF17wO/hqegsUKGCeGUci7VwTAQQQsEHAimfA+hxR74ZKlSol2hI2f/785pmwTjs3\nefJkG5xIQzwBbTCnVca33XabqUKOt4mXCCCAAAI+ClhxB6x9frV6eMGCBaYbijbI0unr9Lmtjoik\nk6936tTJxyyxWygFtK+vju2so5Xp7ET+LNoI7tdffzV/XOXNm9efQ9kXAQQQcJ2AFXfAS5culfHj\nx4uOoKTjQWtjnjx58ogOSzlq1CiZM2eO6+CdmCF9NKBlpf21//Of/5jaCl/zcfbsWenQoYMsWrRI\nCL6+qrEfAgi4WcCKAKxVzStWrPDqPH/+/IAncfd6QlYGLKB/EGlL8iZNmvh1Dp2WsFevXmZ0qyee\neMKvY9kZAQQQcKuAFVXQOv5zy5YtRSdlKFu2rOTOnVu0i4uOC6yNsnSsaJbICuhjAK0+3rhxo19z\n+2r5afewpk2bmpmuIpsLro4AAgjYI2BFAK5ataqZjGHNmjWyd+9e8zy4UKFC5rmvPmv0t3+pPbzu\nSYl27ypdurRfwVdznzFjRuGu1z2fA3KCAALBE7AiAGt2smbNKnXr1g1ezjjwJ23qAAAWzklEQVRT\n0ARWrVolS5YsMc9/g3ZSToQAAghEuYAVz4CjvAyszv6uXbukVatWZqaj6tWr+5xWrc349ttvfd6f\nHRFAAIFoEyAAR1uJ+5lfnRBDG169/vrrPh2pXY30GO1SVq1aNZ+OYScEEEAgGgWsqYKORnzb8zx6\n9Gj57LPPZPny5T4l1ePxyIABA0y/7TfeeEMyZMjg03HshAACCESjAAE4GkvdhzxPnz5dXnnlFZk6\ndarXWaq8nUKnJWzfvr0ZZjJ9eipXvBmxDgEEEIgVIADHSvB/nMD27dtNC/QhQ4aIPxMl+DptZNyF\neIEAAghEsQC3KVFc+MllvUWLFvLQQw/J448/ntwurEcAAQQQSKMAATiNgG47/Pz583Lo0CEzcpUv\nedOhKceOHSt6HAsCCCCAgO8CBGDfrVy/5/79+6VcuXJSvnx5M+hGahn+4IMP5OOPP5Y2bdpIlixZ\nUtud7QgggAAC8QQIwPEwov3lc889Z4LvV199lSqFBl5tHT18+HDJlStXqvuzAwIIIIBAQgEaYSX0\niNp3Ol6zVievXLkyVQOdKrJWrVrSvHlzuhqlqsUOCCCAgHcBArB3l6hau3v3bnn66adl4sSJojNT\npbboZBn6w4IAAgggELgAVdCB27niSJ2H+dZbbzUtnps1a+aKPJEJBBBAwAkCBGAnlFII0zh06FDp\n27evDB48OMWrfPnll9K2bVvTQjrFHdmIAAIIIOCTAFXQPjG5dyedd7lDhw4pZlCfC48aNUpGjhwp\nxYsXT3FfNiKAAAII+CbAHbBvTq7cS2cs0rmWM2XKlGz+1q1bJ++9954JvjopAwsCCCCAQHAEuAMO\njqMjz9K7d29p0KBBimnXGY3ef//9FPdhIwIIIICA/wIEYP/NXHHEkiVLZOvWrQRXV5QmmUAAAScK\nUAXtxFILQppnzpwpN998s5QpUybJ2f7880/ZtWuX6P8sCCCAAAKhESAAh8bV6rMuWrRIli1bJjrp\nQuJFZ0J6+OGHRQfbyJYtW+LNvEcAAQQQCJIAAThIkE46jQ668cgjj0jr1q0TJHvnzp3ywgsvSJ8+\nfYSpBRPQ8AYBBBAIugDPgINOavcJ3333XTl79qw8//zzkjlz5rjE/vbbb/Lpp5/KoEGD5Nprr41b\nzwsEEEAAgdAIEIBD42rtWV988UVzl5s1a9YEaSxUqJD06NEjwTreIIAAAgiEToAq6NDZWndmvfvN\nmDGjdO3a1bq0kSAEEEAg2gQIwFFS4nv27BHt96t9etOn//9iP3funMybN090MgYWBBBAAIHwChCA\nw+sdkavpOM//+te/5IEHHpB69eqZNOgz38cff9y89tYVKSIJ5aIIIIBAFAnwDDgKCvudd96RMWPG\nSNOmTU1ujx8/Lt27dxed/ahJkyZRIEAWEUAAAfsECMD2lUlQU6QjXp08eVIaNWpkznvhwgXZtGmT\nvPLKK1KuXLmgXouTIYAAAgj4LkAA9t3KcXvq7EX688QTT8RNuKBdj+rWreu4vJBgBBBAwG0CPAN2\nW4nGy8+sWbPknnvuSXWu33iH8BIBBBBAIEwCBOAwQYfrMleuXJGDBw+Kjna1Y8cO879WO48ePVrm\nz58frmRwHQQQQACBVAQIwKkAOW3z3XffLXfddZdoQ6vPP/9c/vGPf8izzz4rFy9ejHsO7LQ8kV4E\nEEDAjQI8A3ZRqV66dMk0sNKJFipUqCB656tjO5csWdK0enZRVskKAggg4HgBArDji/D/M6CNrcaO\nHSv58+c3wVfX6oxG2gc4Q4YMLskl2UAAAQTcI0AAdnhZ6sQKOsLVJ598YiZTqFGjRlyOChYsGPea\nFwgggAACdgkQgO0qD79T06ZNG9FhJvV5b8WKFf0+ngMQQAABBCIjQACOjHtQrqrjOn/99deyYcMG\nKV68uDnnwIED5fTp0/Lqq6+aiReCciFOggACCCAQdAECcNBJw3dCnVpwwIABccF3yJAhomM86zqd\n9YgFAQQQQMBeAb6l7S2bFFO2ZcsWuXz5srRr187sN2nSJNPa+ZlnnhEd7YoFAQQQQMBuAQKw3eXj\nNXXff/+93HvvvXLnnXdKunTpzD463CQLAggggIBzBBiIwzllFZfSli1bmkE1pk6dGreOFwgggAAC\nzhLgDthB5bV3716ZMmWKGeVKpxf8/fff5dSpU1K6dGn6+jqoHEkqAgggoALcATvoczBs2DBZuHCh\nzJw5U2bMmCE9e/YU7evLQBsOKkSSigACCPxPgDtgh3wUPB6PCbwTJkwwLZ3nzZsnb7/9tuTNm9ch\nOSCZCCCAAALxBbgDjq9h8evnn39ecuTIIaVKlZIjR47IO++8Y4adtDjJJA0BBBBAIAUB7oBTwLFl\nk04v+NFHH0n//v2lSpUq5seWtJEOBBBAAIHABLgDDswtbEetXbtWqlWrZqYYbN26ddiuy4UQQAAB\nBEIrQAAOrW+azq5Vzffff780bdpUGjRoIOfOnUvT+TgYAQQQQMAeAesCsM5pe+LECXuEIpSSw4cP\ny0033WS6GMW+zpMnT4RSw2URQAABBIItYEUA1onj+/TpYxoY6TCKOqetNjjS2X0mT54c7Dxbf76j\nR4+aKucSJUqYbkY6wUKZMmWsTzcJRAABBBDwXcCKRlhdu3YVvctbsGCBCTQafHUy+a1bt0q3bt3k\nr7/+kk6dOvmeK4fvqXnVyRRef/110fl9ufN1eIGSfAQQQMCLgBV3wEuXLpXx48dL5cqVJWfOnGZ8\nYw06NWvWlFGjRsmcOXO8JN19q7T6/dFHHzVTDP7nP/8xz30Jvu4rZ3KEAAIIqIAVAVirmlesWOG1\nRObPny+FChXyus1tKwcNGiSrV6+WVatWSfny5d2WPfKDAAIIIBBPwIoqaO3fqhMMjBgxQsqWLSu5\nc+c2Yxxv27ZN9K5Qh190+7Jy5UoZO3as3HjjjWaaQbfnl/whgAAC0S5gRQCuWrWqbNiwQdasWSM6\n4YA+D9a7Xn0WWrt27bgp99xcWA8++KDJ83PPPSfXXXedm7NK3hBAAAEEYgSsCMBaElmzZpW6desm\nKZQdO3aY/q8apFNbdHq+2bNne93thx9+MPPnasOm2EXn1dUxlmMXvUb69H/Xyodru/bv1VmN+vXr\nZxy++eYbk6Tq1asnmGhBB+W4cuVKbHKF7fjEn4iDzwe/H3w/BP79GPfFGsYX6WIC0N8RKIwX9vVS\nGjD37dsnEydOTPUQbS2tP94WnT1I+xe3atUqbrO2to6/aCCMzxGO7ZcvX5aTJ0+aavds2bIluH6u\nXLniJ0/OnDnD9ngfV3z4fMT/BeH3g++H+N/fvn4/6KxyGhd03IVwL9YH4GCBTJ8+3QTgjh07BuuU\nnAcBBBBAwOECkQzAf9e3WoLISFiWFATJQAABBBAIqYAVATjaRsI6fvy4tGnTRrRanAUBBBBAIDoF\nrGiEFU0jYWljq+7du5uhJps3bx6dnzpyjQACCCBgx0Ac0TIS1tmzZ0XHdW7SpInp98znDwEEEEAg\negWsqIKOlpGwtFX1gAEDpFmzZtH7iSPnCCCAAAJGwIoqaEbC4tOIAAIIIBBtAlYEYDePhKX90nbt\n2mUG2ChVqlS0fb7ILwIIIIBAMgJWBGBNW3IjYSWTbkesvnjxovTu3VuKFi0qOsQkCwIIIIAAArEC\nVjwDjk2Mm/7X/sx9+/aVAgUKSK9evdyUNfKCAAIIIBAEAWvugIOQF6tOodMoPvzww1KlShWr0kVi\nEEAAAQTsEIiaoSg3btwo99xzj/gyqUMkimb58uWSPXv2SFw6otfUsbszZcqUYNKJiCYoTBfXQfPP\nnz8vOv53tC3aHS/xOOvRYKCPpLRNSObMmaMhuwnyqJ/1OnXqJFhny5vdu3fLZ599JiVKlAh7kqIm\nAIdd1s8L6ofziy++8PMo5+/+wgsvyH333Sc1atRwfmb8yMGePXvktddek/fee8+Po9yxa/369WXR\nokWSMWN0VcDpyHe//fabdO7c2R0F6UcuovX7LTUingGnJsR2BBBAAAEEQiBAAA4BKqdEAAEEEEAg\nNQECcGpCbEcAAQQQQCAEAgTgEKBySgQQQAABBFITIACnJsR2BBBAAAEEQiBAAA4BKqdEAAEEEEAg\nNQG6IaUmFKbtv/76qxQrVixMV7PnMsePHzd9QrNkyWJPosKQEh0p7cSJE1KoUKEwXM2uSxw+fNgM\nz2pXqkKfGu3/fPnyZcmdO3foL2bZFaL1+y21YiAApybEdgQQQAABBEIgQBV0CFA5JQIIIIAAAqkJ\nEIBTE2I7AggggAACIRAgAIcAlVMigAACCCCQmgABODUhtiOAAAIIIBACAQJwCFA5JQIIIIAAAqkJ\nEIBTE2I7AggggAACIRAgAIcAlVMigAACCCCQmgABODUhtiOAAAII+Czg8XjMgCM+HxDFOxKAw1j4\nX3zxhdSqVUtKly4tTZs2NSMhebu8r/t5O9bGdTri04MPPijXXHONVKpUSVavXu01mVu3bpWHH35Y\nbrzxRrnzzjtl+vTpXvdz0soBAwZI5cqVTZnr69SW9u3bS4cOHVLbzfrtvn6Gv/32W6lWrZpUqFBB\nGjVqJNu2bbM+b6kl0Ncyf/XVV6V69epyyy23yNChQ1M7rSO2X7lyxfyuDxkyJNn0+vrZSPYEbtoQ\n89cKSxgEfvvtN0/MUJOeTZs2eS5cuODp3r27p127dkmu7Ot+SQ60eEXz5s09r732mifml9OzYsUK\nT5EiRTznzp1LkuL69et7PvjgA7P+4MGDnsKFC3tihi1Msp9TVnzyySeef/3rX56TJ096Yobi88T8\nYeFZuHBhssmfP3++J3/+/J6YIJzsPk7Y4Otn+K+//vKUKVPGs2bNGpOtadOmeZo1a+aELCabRl/L\nfMmSJZ6Y4Gu+C9Thuuuui3NI9uSWb1i3bp0n5gbDky9fPk/MHyFeU+vrZ8PrwS5cyR1wmP6aivlw\nyvXXX2/uhjJlyiRdu3aV2bNnJ7m6r/slOdDiFYsXL5bOnTtLunTppE6dOlKyZEn56quvEqRY/3LW\nffQOWJfixYtLrly55Pvvv0+wn5PeaL5bt24tefLkMWMfa97++9//es3C77//Lq+//rr5XHjdwUEr\nff0Mx/wxIuXKlTN3gKdOnZKHHnpIZs6c6aCcJk2qr2V+7NgxSZ8+veh3gY6DnjlzZjl06FDSEzpo\nTcwfz/L000/H/Q57S7qvnw1vx7pxHQE4TKW6f//+BJMtxNwFin7pnD9/PkEKfN0vwUEWv9HqZ81j\nzJ1dXCqLFi0qR48ejXuvL/TL6L777jNfSPp++fLlpoq+Zs2a+taRS+Ky1HwfOXLEa146deokr7zy\niuTMmdPrdietTJzv5D7r+/btM5+L2rVrm0kpypYtK1u2bHFSVpOkNXHekytz/ayXKFFCbrvtNrn1\n1lsl5g5Y7rnnniTnc9KKt956S2Jqu1JMcmKf5D4bKZ7ERRsJwGEqTL3DyZEjR9zVsmXLZl7HVMXG\nrdMXvu6X4CCL3yTOjyZV837mzJlkU/3TTz9JmzZtZPTo0ZI3b95k97N9Q+K8Z8+eXXRGnMTLRx99\nZEz+/e9/J97kyPeJ853cZz2mOlJiqmylY8eO5nPfsGFDGTRokCPzHJvoxHlPrsx37dol+jm/4YYb\nTNsIfX3gwIHY07j2/8Q+yX02XAuQKGME4EQgoXpbsGBB+eOPP+JOf/r0acmaNavEPC+JW6cvfN0v\nwUEWv0mcH02qOmgVs7dl+/btppr65ZdfTrEqy9uxtq1LnHdv+dYvpG7duskdd9whMc+ATSMkvTOM\neS5qW3Z8Tk/ifCf3Wdc/rrTxVcuWLc3jht69e8u8efMkpo2Ez9eybcfEefdW5prm4cOHS+PGjWX8\n+PGiVbd6J6zr3L4k9knus+F2h9j8ZYx9wf+hFdDnnnv37o27iL4uVapU3PvYF77uF7u/7f/rl6z+\nlfvLL7+YZ7+aXs37VVddlSTpu3fvlnr16smLL75o7oqS7OCwFVqWGkxjF29lHtNAyzwH1S9iXfQ5\nYEyjHJk6dao4tfrd18+w7qfP+WMXfR76559/irYHcOriS5lr3rQq9v7774/LpraGjmmEFvferS98\n/Wy4Nf9J8uXChmVWZklbOmqr3mXLlnn0ddu2bT3PP/+8SWvMpPSezz77zLxOaT8rM+ZDoh577DFP\nTKMzz8WLFz0xjWxMi09tCa6LtorWlpG6xDwL8zz33HOemLvCuJ+Y58dmmxP/WbRokSemC5JHW3Tv\n2bPHE9PgyPPdd9+ZrMR0ufL8+OOPSbIV033D8a2gU/oMx/+sx9wdegoUKOCJ6YpkHGIaoXlingcn\nMXHSCl/LfMKECZ4HHnjA9AzQHgExjx/iegA4Kb/e0hrTmDJBK+j4ZZ7SZ8Pbudy+TtyeQZvyp10U\nYhrZeGIaX3jq1q3rial+Mcn78ssvPTHPiuKSmtx+cTs47IUGn4oVK3piqp09MQ1tTNCNzYJ2SVqw\nYIFn7dq1npi/DpP8vP/++7G7Ou5/7XalXc1iagE8MY1xPP369YvLQ5cuXTyPPvpo3PvYF24IwJqX\n5D7DiT/rc+fONTb6uShfvrwnphYklsKR//ta5vrHR8yzb/MHWkz/eI8GrZj2AY7Mc+JEJw7Aics8\nuc9G4vNEw/t0mskkt8WsCJnApUuXRJ97JH72m/iCvu6X+Dib32ujm0KFCtmcxJCkTZ8DalcT/Ymm\nxdfPsH4FabccN302fC1z7SGg3fO0G1I0Lb5+NtxuQgB2ewmTPwQQQAABKwVoBW1lsZAoBBBAAAG3\nCxCA3V7C5A8BBBBAwEoBArCVxUKiEEAAAQTcLkAAdnsJkz8EEEAAASsFCMBWFguJQgABBBBwuwAB\n2O0lTP4QQAABBKwUIABbWSwkCgEEEEDA7QIEYLeXMPlDAAEEELBSgABsZbGQKAQQQAABtwsQgN1e\nwuQPAQQQQMBKAQKwlcVCohBAAAEE3C5AAHZ7CZM/BBBAAAErBQjAVhYLiUIAAQQQcLsAAdjtJUz+\nEEAAAQSsFCAAW1ksJAoBBBBAwO0CBGC3lzD5QwABBBCwUoAAbGWxkCgEEEAAAbcLEIDdXsLkDwEE\nEEDASgECsJXFQqIQQAABBNwuQAB2ewmTPwQQQAABKwUIwFYWC4lCIPgCe/bskUqVKsmuXbvMySdP\nnizNmzcXj8cT/ItxRgQQSFUgXcwvH799qTKxAwLuEOjRo4fs3LlTxo8fL5UrV5bFixdLtWrV3JE5\ncoGAwwQIwA4rMJKLQFoEzp49KxUqVJDcuXPLPffcIwMGDEjL6TgWAQTSIEAVdBrwOBQBpwnkyJFD\nOnXqJJs3b5YuXbo4LfmkFwFXCXAH7KriJDMIpCxw8uRJueGGG8xPsWLFZOrUqSkfwFYEEAiZAHfA\nIaPlxAjYJ9CzZ09p2LChzJo1S5YtW2aeAduXSlKEQHQIZIyObJJLBBD4/PPPZe7cubJ9+3bJkyeP\nDBs2TDp27Giqo3PmzAkQAgiEWYAq6DCDczkEEEAAAQRUgCpoPgcIIIAAAghEQIAAHAF0LokAAggg\ngAABmM8AAggggAACERAgAEcAnUsigAACCCBAAOYzgAACCCCAQAQECMARQOeSCCCAAAIIEID5DCCA\nAAIIIBABAQJwBNC5JAIIIIAAAgRgPgMIIIAAAghEQIAAHAF0LokAAggggAABmM8AAggggAACERAg\nAEcAnUsigAACCCBAAOYzgAACCCCAQAQECMARQOeSCCCAAAIIEID5DCCAAAIIIBABAQJwBNC5JAII\nIIAAAgRgPgMIIIAAAghEQOD/AJ8t4ZLNaVpdAAAAAElFTkSuQmCC\n"
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%%R -i P\n",
"plot(ecdf(P))\n",
"abline(c(0,1), lty=2)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def simulate_misspec(n=1000, p=20, s=3, signal=2, lam=10, truth=None):\n",
" \n",
" X = np.random.standard_normal((n, p))\n",
" beta = np.zeros(p)\n",
" beta[:s] = 2 / np.sqrt(s)\n",
" \n",
" linpred = X.dot(beta) \n",
" pi = normal_dbn.cdf(np.log(1 + np.exp(linpred)))\n",
" Y = np.random.binomial(1, pi)\n",
" \n",
" L = ROSI.logistic(X, Y, lam)\n",
" L.fit(solve_args={'min_its':200})\n",
" return L.summary(truth=beta)\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"import regreg.api as rr\n",
"n, p, s = 100000, 20, 3\n",
"beta = np.zeros(p)\n",
"beta[:3] = 2 / np.sqrt(s)\n",
"\n",
"beta_star = []\n",
"nsim = 100\n",
"for i in range(nsim):\n",
" X = np.random.standard_normal((n, p))\n",
" linpred = X.dot(beta)\n",
" pi = normal_dbn.cdf(np.log(1 + np.exp(linpred)))\n",
" Y = np.random.binomial(1, pi)\n",
" loss = rr.glm.logistic(X, Y)\n",
" beta_hat = loss.solve()\n",
" beta_star.append(beta_hat)\n",
"beta_bar = np.mean(beta_star, 0)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"P = []\n",
"for _ in range(100):\n",
" P.extend(simulate_misspec(truth=beta_bar)['pval'])"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAYAAAB91L6VAAAEGWlDQ1BrQ0dDb2xvclNwYWNlR2Vu\nZXJpY1JHQgAAOI2NVV1oHFUUPrtzZyMkzlNsNIV0qD8NJQ2TVjShtLp/3d02bpZJNtoi6GT27s6Y\nyc44M7v9oU9FUHwx6psUxL+3gCAo9Q/bPrQvlQol2tQgKD60+INQ6Ium65k7M5lpurHeZe58853v\nnnvuuWfvBei5qliWkRQBFpquLRcy4nOHj4g9K5CEh6AXBqFXUR0rXalMAjZPC3e1W99Dwntf2dXd\n/p+tt0YdFSBxH2Kz5qgLiI8B8KdVy3YBevqRHz/qWh72Yui3MUDEL3q44WPXw3M+fo1pZuQs4tOI\nBVVTaoiXEI/MxfhGDPsxsNZfoE1q66ro5aJim3XdoLFw72H+n23BaIXzbcOnz5mfPoTvYVz7KzUl\n5+FRxEuqkp9G/Ajia219thzg25abkRE/BpDc3pqvphHvRFys2weqvp+krbWKIX7nhDbzLOItiM83\n58pTwdirqpPFnMF2xLc1WvLyOwTAibpbmvHHcvttU57y5+XqNZrLe3lE/Pq8eUj2fXKfOe3pfOjz\nhJYtB/yll5SDFcSDiH+hRkH25+L+sdxKEAMZahrlSX8ukqMOWy/jXW2m6M9LDBc31B9LFuv6gVKg\n/0Szi3KAr1kGq1GMjU/aLbnq6/lRxc4XfJ98hTargX++DbMJBSiYMIe9Ck1YAxFkKEAG3xbYaKmD\nDgYyFK0UGYpfoWYXG+fAPPI6tJnNwb7ClP7IyF+D+bjOtCpkhz6CFrIa/I6sFtNl8auFXGMTP34s\nNwI/JhkgEtmDz14ySfaRcTIBInmKPE32kxyyE2Tv+thKbEVePDfW/byMM1Kmm0XdObS7oGD/MypM\nXFPXrCwOtoYjyyn7BV29/MZfsVzpLDdRtuIZnbpXzvlf+ev8MvYr/Gqk4H/kV/G3csdazLuyTMPs\nbFhzd1UabQbjFvDRmcWJxR3zcfHkVw9GfpbJmeev9F08WW8uDkaslwX6avlWGU6NRKz0g/SHtCy9\nJ30o/ca9zX3Kfc19zn3BXQKRO8ud477hLnAfc1/G9mrzGlrfexZ5GLdn6ZZrrEohI2wVHhZywjbh\nUWEy8icMCGNCUdiBlq3r+xafL549HQ5jH+an+1y+LlYBifuxAvRN/lVVVOlwlCkdVm9NOL5BE4wk\nQ2SMlDZU97hX86EilU/lUmkQUztTE6mx1EEPh7OmdqBtAvv8HdWpbrJS6tJj3n0CWdM6busNzRV3\nS9KTYhqvNiqWmuroiKgYhshMjmhTh9ptWhsF7970j/SbMrsPE1suR5z7DMC+P/Hs+y7ijrQAlhyA\ngccjbhjPygfeBTjzhNqy28EdkUh8C+DU9+z2v/oyeH791OncxHOs5y2AtTc7nb/f73TWPkD/qwBn\njX8BoJ98VQNcC+8AAEAASURBVHgB7Z0FeBTX18YPTiAES3AprsGKFilQ3IpDoVCkLV68uHtxKRS3\nQilOi1uAQqFAcW+A4C7BJWG/Pff/7bJJdpPdZHfmzsx7nydkdubO3HN+d8K7186NZTInQgIBEAAB\nEAABEFCUQGxFS0NhIAACIAACIAACggAEGC8CCIAACIAACKhAAAKsAnQUCQIgAAIgAAIQYLwDIAAC\nIAACIKACAQiwCtBRJAiAAAiAAAhAgPEOgAAIgAAIgIAKBCDAKkBHkSAAAiAAAiAAAcY7AAIgAAIg\nAAIqEIAAqwAdRYIACIAACIAABBjvAAiAAAiAAAioQAACrAJ0FAkCIAACIAACEGC8AyAAAiAAAiCg\nAgEIsArQUSQIgAAIgAAIQIDxDoAACIAACICACgQgwCpAR5EgAAIgAAIgAAHGOwACIAACIAACKhCA\nAKsAHUWCAAiAAAiAAAQY7wAIgAAIgAAIqEAAAqwCdBQJAiAAAiAAAhBgvAMgAAIgAAIgoAIBCLAK\n0FEkCIAACIAACECA8Q6AAAiAAAiAgAoEIMAqQEeRIAACIAACIBAXCEAABOQi8PbtW9q2bRt5eXlR\n5cqVrcb99ddfdOzYMUqTJg3VrFmTvL296cCBA3TkyBHq1q0bnThxgq5fv27NHy9ePEqSJAkVKlRI\n5LVc2LhxI3348EF8zJEjB6VMmZIOHTpkuUxx4sShhAkTUs6cOSljxozW8+46eP78ubCZ7StatKjw\nk5998OBBevDggSgmRYoUVKZMGVq8eDH5+flRjRo13FU8ngMC8hAwIYEACEhF4O7duybz/xCmTJky\nWe2aO3euOMfn+efhw4em9+/fm8wCavr6669FvtatW4fJY8lrFlHT7t27rc+KHz++NV/fvn1Nmzdv\ntn623MO/zQJpGjZsmMks1tZ7Y3pw/PhxU7p06azlsf2XL18Wj61UqZL1fOnSpcW5IUOGmMxibAoO\nDo5p0bgfBKQjgC5o8/80SCAgO4Ht27cLEydPnkz37t0TrdYFCxbQf//9Ry1btgxjfu3atWnSpEk0\nYsQI0cK8ceMGtW3bNkwebmHu3buXvv/+e+t5bu1OmTKFxo4dS2ZRp9DQUDILIP3777/WPLYHn332\nGa1YsYJCQkJsT0d63LNnT7p9+zb179+f2rVrJ+wfPHiwuIfLXrduXZj72Y7Hjx8Lf8JcwAcQ0AMB\n6b4SwCAQ0CABc7euydwNbDILi8ncVRzBA3MXr2ngwIEmswCZNm3aZHr58mWYPKtXrzZ17tzZNG7c\nONOFCxdES9DSAp43b54pe/bs4pxZtEwrV64U92bJkkW0Js1CKT5bWsD8DEsyd/eazF3Z4t5z586J\n09wCNndjW7JYW8Dm7mDrOT6oX7++uG/o0KFhzls+mP//E9fTp09vGjVqlMncfWy5ZPf31atXRX62\nm9ObN29M5i5y0dJmOzndvHlT5LG0gPlcqVKlTD4+PqanT5/yRyQQ0A0B0o0ncAQEVCLQo0cPIRpx\n48Y18U+sWLFM5hao1Zpp06aZYseObc3DwlWxYkXrdRZli5iZx16tXbQWAc6XL5/1OuerVq2aEDs+\nbtGihfU59gTY3PoVNnHewMBAkdcZAWZRL1++vCh39OjR1jJsD06dOmXq2rWr6CLm57Ptbdq0MZnH\nom2zWY+5G5zzValSxXquQIEC4tz58+fFOXsCzN3QfN/OnTut9+EABPRAAF3Q5r9sJBCILoGTJ08S\ndwtz9+2VK1fI3Mqj1KlT06BBg8g8Tiu6i/v06UNmUaY9e/aQeSyTqlatKiZM8edr166Jbl+zQIvr\n9+/fF5OPbO3Zt28flStXTpzasGGD6PY1i5/4bBZp26zimLuFueu2Xr16VLhwYdFFnD9/fsqWLVuE\nvLYn2H6zoFOTJk2oYMGCwh6eKOVoApS/v7+wnbuUucySJUsSd4vzpC/uNg6f7ty5I04lSpTIeokn\nmnGyXLNesDmw+Hj27FmbszgEAe0TwCxo7dchPFCRAM9MNn8Tp08++YTMLTxhSdasWenvv/+mw4cP\nC+F9/fo1FStWjD7//HNx/Y8//iBzK1Qcs6DyWCuLl+V6q1atyNzNbPWKx2st+ZMnT05JkyYlFn5O\n9mYpmyc6Ef9wYmHnmdRz5swRnyP7h0Xz119/tWZhwR4/frwQY+tJOwc8xszlmVux4irPvDb3BETI\nybOrOTEvS7Ic28tvyZMhQwZxeObMGcsp/AYBXRCI+FeiC7fgBAgoQ+DFixeiIBYf89irtdA8efKQ\neYyTWHw5sYhakkVM+fOzZ8/EaV5aZEmZM2e2HDr8zS1lTizI4RNPcOrQoYNYesRiza1vZ5K5O5jM\n49NC7Pm53PqNLC1fvpzMs7PFZC4WUv7iwb0B5m5oMo/ZRrjV4qOFGWfgJUmczDOjxW97//AyKU7c\ne4AEAnoiAAHWU23CF8UJmJfOUL9+/UQL17ycR5TPs4uTJUtGefPmFet2+SSvceW1t9wi/fnnn4lb\nwT/88IMQLb7O63u5JcytxICAAD4VacqVK5e4blk3a5uZRdfSarQ9H9UxfzFw5b7mzZuLR1aoUIHM\nY8HEs6/ZP0fJPJFMXOeW7Lt378QXFO6y5/XM5olcjm6zfknhNctIIKAnAhBgPdUmfFGcwKeffirG\nfHfs2CGCYXCXsHltrega5m5Z7nrmYBNHjx6lEiVKiG5mHic1r+GlRYsWUapUqcg8yYp4fJPHTs0T\nrESrMipH+B5OkY2dRvWMmF7nli4LL7ecnUksshxA5M8//6SvvvpKtH65l4CfwYE/HCUem+YEAXZE\nCOe1SgACrNWag91SEODuXV6jy63BqVOnihZs8eLFxeQkbolyWrVqleiW3WOedMUtXb7evXt3Sps2\nrfV6o0aNiFuG3CLk53Bkq8gSd3FzF7FtBKvI8nvi2vz5811+7OzZs8m8nIjWrl0r7uUvHJZ1wI4e\nxl9kmDOPkyOBgJ4IxDKP3XycEaEnz+ALCChMgIWFE3c/20s83sldr5YxzfB5eDYxt4gjm5Bke893\n331Hy5YtoydPnlCCBAlsL0V6zHl5TFrN1jPPEOfu9vBj2Ldu3RLd4OZ1wLR//37hB3e38w932yOB\ngJ4IOB6w0ZOX8AUEFCDAwutIfLl4nh3sSHz5Ok9EclZ8OT93dbOgR0eY+MtCs2bN6Pfff+dHKZ58\nfX0jiO+YMWOoU6dOYWzh1u+lS5eIl3IhgYDeCECA9Vaj8McwBHiZEI+lzpgxwyWfeTYyf1HgyV7m\nOMwu3evJzLyZxD///CM2m7B8UZk+fbpYF80tYiQQ0BsBdEHrrUbhj6EI8MxpjsXsShe0lgBxC59n\nVrvSM6Al/2CrsQlAgI1d//AeBEAABEBAJQLoglYJPIoFARAAARAwNgEIsLHrH96DAAiAAAioRAAC\nrBJ4FAsCIAACIGBsAhBgY9c/vAcBEAABEFCJAARYJfAoFgRAAARAwNgEIMDGrn94DwIgAAIgoBIB\nCLBK4FEsCIAACICAsQlAgI1d//AeBEAABEBAJQIQYJXAo1gQAAEQAAFjE4AAG7v+4T0IgAAIgIBK\nBCDAKoFHsSAAAiAAAsYmAAE2dv3DexAAARAAAZUIQIBVAo9iQQAEQAAEjE0AAmzs+of3IAACIAAC\nKhGAAKsEHsWCAAiAAAgYmwAE2Nj1D+9BAARAAARUIgABVgk8igUBEAABEDA2AQiwsesf3oMACIAA\nCKhEIK5K5apS7OrVqykkJESVslEoCIAACICAfARSpUpFFStWVMWwWCZzUqVkhQtds2YNTZw4kb75\n5huFS0ZxIAACIAACrhDghtK6devI19eXKlSo4MqtLuedNm0aLVu2jAoVKuTyvTG9wTAtYK7Qli1b\nUrt27WLKDPeDAAiAAAh4iMC7d++oT58+VLVqVerdu7eHSvn42EuXLtGHDx8+nlDwyDACrCBTFAUC\nIAACIBANAqGhoaLl+9VXX1Hx4sWj8QRt3QIB1lZ9wVoQAAEQ0C2BOHHiUJMmTXTrX3jHMAs6PBF8\nBgEQAAEQAAEFCEgnwDxW++TJEwVcRxEgAAIgAAIyENi7dy/t2rWLDDIn2IpcCgHmQff+/ftTxowZ\nKX78+JQiRQpKnDgx5c+fnxYuXGg1FgcgAAIgAAL6IcCCO2rUKNq0aROVK1eOYsWKpR/nnPBEijHg\nLl260N27d0UlZM2aVYjvs2fP6Ny5c9StWzd68+YNdejQwQl3kAUEQAAEQEArBMaPH0+PHz+mMWPG\nULx48bRittvslEKAt2/fTgcPHqQ0adJYHUuaNCmVKlWKpk6dSkOGDIEAW8ngAARAAAS0T+DChQvU\nvn178vb2ptixpeiMVRyqFF5zV3NAQIBd5zdu3Eh+fn52r+EkCIAACICANgnkzp2bfHx8DCu+XGtS\ntICHDx9OzZo1o8mTJ1O2bNlEpQQHB9P58+dF6MjNmzdr8w2D1SAAAiAAAiDggIAUAly4cGE6fvy4\n6IYOCgoS48Hc6uVxXyMOzDuoK5wGARAAAU0T4B7NM2fOUNeuXcnLy0vTvrjDeCkEmB1JmDCh3Zif\nFy9epFevXhGLdFRp3759dOjQIbvZ/vrrLxHrk1vVlpQrV64w3R9clm1IMlwHH9uxKbwf+PvA/w8f\nQzY6+/8j92byJKs7d+4Q/w01bNiQFixYQC9fvrT8Vyx+88oX26TUdR6LVitJI8COAKxatYquXbtG\nc+fOdZTFep4ncfn7+1s/2x78/fffdO/ePUqSJInt6TDHPBkgsnVouA4+eD8c792Cvw/8fdj7+2ja\ntKkYTvzkk0+IV7zwe8KJZz3bLjvi5ai2SanrJ06csC1W0WPD7IbUvXt30bX922+/KQoYhYEACICA\nUQlMnz6dBgwYIHom8+TJE0ZwZWHSs2dPat68ORUpUkRxk6RvAStOBAWCAAiAAAjEmEDp0qXp6tWr\ntH79esqbN2+Mn6fHB0CA9Vir8AkEQAAEVCLAgTWWLl0qJlvxVn+pU6dWyRL5i5VCgCdOnEjv3793\nSIvXi9WtW9fhdVwAARAAARBQn8CpU6eoUqVKxOO9+fLlE8tI1bdKXgukEGBeejRjxgz65ptvRBjK\n8LgQiCM8EXwGARAAAfkItGnThjJlyiTi+k+ZMoXSp08vn5ESWSSFAPNAPU/v55+ff/5ZIjwwBQRA\nAARAICoCjx49ovr169OtW7eoXr16NHDgQEqXLl1Utxn+uhShKLkWxo0bR7wBw4sXLwxfKQAAAiAA\nAlogwEtEeRtBHibkYUQe8505cybE18nKk0aAeW3YsmXLrGvEnLQf2UAABEAABBQmwGLLQ4aff/45\nzZ8/n3r16kUcayGyOAsKm6iJ4qTogtYEKRgJAiAAAiBAT58+FdvE7tmzh7Zu3SpmOQcGBoreS0uQ\nDWByjoA0LWDnzEUuEAABEAABNQjwlrG8RSzPbuawklu2bBGznFu3bi3iOkN8Xa8VCLDrzHAHCIAA\nCBiKAC8v4qWgHFyD4zlv27ZN+D9o0CAaNmyYwxDAhoIUDWfRBR0NaLgFBEAABIxA4PXr19SkSRP6\n559/qEWLFjRhwgTh9oMHD+jw4cNi+WiGDBmMgMIjPkKAPYIVDwUBEAABbRN4+/YtlS9fXuxGx9vF\n2i4r4tgMrVq10raDElgPAZagEmACCIAACMhGgLcN5K0EuaXr4+Mjm3m6sAdjwLqoRjgBAiAAAu4j\nwHvx8rIiDoxkEd8nT56IJUc84xnJPQQgwO7hiKeAAAiAgC4I8Ljv0KFDxZreL774QvjEEa6+//57\nEWYye/bsuvBTBifQBS1DLcAGEAABEJCEQJ06dcTe6SNGjBAW3b17V6z7bdu2LVWuXFkSK/VhBgRY\nH/UIL0AABEAgxgS+++47OnbsGP37779iR6N3797RvXv3RHhJbIoTY7wRHgABjoAEJ0AABEDAOARC\nQ0PpypUrNGvWLFqzZg3t379fiC8TiB8/PhUsWNA4MBT2FAKsMHAUBwIgAAKyEOAWbpEiRYQ5xYoV\no927d1PevHllMU/3dkCAdV/FcBAEQAAEIhK4fv069enTh3jS1ZkzZ0Q4Sc7FO9KNGjVKhJ3k8WAk\nzxHALGjPscWTQQAEQEBKAkFBQcQtXt7VaMOGDVbxff78uZhwlTNnToL4er7q0AL2PGOUAAIgAAJS\nEGDB5ahWHTp0oMKFC9Pq1autdvHa3549e1Lx4sWJN1hA8jwBCLDnGaMEEAABEFCVAHczjx49mlat\nWkUJEiQgXt/Le/jaJi8vL5ozZ47tKRx7mAAE2MOA8XgQAAEQUJMA72Q0efJkMcFq6dKlVK5cObvm\nxI6NEUm7YDx4EsQ9CBePBgEQAAE1CfAaXh7L5Y0VuLvZnvg+ffqUeOwXSXkCEGDlmaNEEAABEPAo\ngQ8fPoitA3k8N1euXLRkyRIx6cq2UB4P5m7oX375hby9vW0v4VghAuiCVgg0igEBEAABpQgMGzZM\njOfOnTuXqlWrRnHjhv2vPiQkhAYPHkxJkyalH3/8kWLFiqWUaSjHhkDYWrG5gEMQAAEQAAFtETCZ\nTLR3716aNGkSbdmyhcqUKRPBAc6zePFi0R3N4gzxjYBIsRMQYMVQoyAQAAEQ8ByBtWvXWluz3Lq1\nJ75cOgsub6yApD4BCLD6dQALQAAEQCDaBJ49eyaWFd28eVME0eDoVkjaIAAB1kY9wUoQAAEQsEtg\n4MCBInzkhQsXxJiu3Uzmk6dPnxZhJ4sWLUpYcuSIkrLnIcDK8kZpIAACIOAWAhxY47fffiNuAXPs\nZp5Q5ShNnDiRbty4QePGjYP4OoKkwnkIsArQUSQIgAAIxIQAr+9lAd62bZsIHRkvXjyHj5s2bRpd\nvXqVxo8fL6JgOcyIC4oTgAArjhwFggAIgED0Cfzzzz9Us2ZNatiwIZUuXTrSB3Hc5yZNmlDKlCkj\nLEWK9EZcVIQABFgRzCgEBEAABGJOgMNKsvDWq1ePeI1vVIk3XECSlwAiYclbN7AMBEAABKwEHj16\nRPXr16cSJUrQzJkzredxoF0CEGDt1h0sBwEQMAiBI0eOUJUqVcjPz0+ElYxszPfgwYM0a9YsevXq\nlUHoaNdNCLB26w6WgwAIGIDAlClThPhWrFhRRLdKlCiRQ695w4UZM2aIcd/I8jl8AC4oSgBjwIri\nRmEgAAIg4ByB4OBgKlmyJPGmCStWrKCqVatGeuOGDRtozZo1xLOeU6RIEWleXJSDAARYjnqAFSAA\nAiAgCBw+fJhatmxJL1++pCJFihCHmIwTJ06kdFis/f39xezo8BsvRHojLqpKAAKsKn4UDgIgAAIf\nCYwZM4bmzZtHn332GS1duvTjhSiOOAhHZIE4orgdl1UiAAFWCTyKBQEQAAELAd6/l3coGjt2rNif\nt3HjxpZL+K1jApiEpePKhWsgAAJyE3j9+jX9+uuvlD17dpo6dSrNnz+fvvrqqyi7nNmrQ4cOUatW\nrYg3YUDSJgG0gLVZb7AaBEBABwQKFiwoolR16NCBevfu7bRHvNTop59+Io7xnCFDBqfvQ0a5CECA\n5aoPWAMCIKBzArdu3RLLia5fv04PHz6kEydOkCtLho4dOyaWGnFs56xZs+qclr7dQxe0vusX3oEA\nCEhCgLubeUZzoUKFaM+ePZQuXTratGmTS+LLrvDM6GXLlolua0lcgxnRJIAWcDTB4TYQAAEQcIUA\nj+2eP3+ehgwZQp07d3blVuTVKQEIsE4rFm6BAAjIQ4CXFAUEBNCZM2coY8aMLhv29u1bun37NqVJ\nk4a8vLxcvh83yEkAXdBy1gusAgEQ0DABXlbUtm1bMUZbvnx5McN5woQJ0RLfixcvUtOmTYk3Y4D4\navilsGM6WsB2oOAUCIAACESHAE+qWrBggRjjPXnyJO3bt4+yZcsWnUeJewIDA6lPnz7Ur18/Klq0\naLSfgxvlJAABlrNeYBUIgIDGCLx48YIqV65MCRIkoAEDBohoVilTpoy2F/fv36f169fTuHHjKFeu\nXNF+Dm6UlwAEWN66gWUgAAIaIbB8+XKxjjdLliz0xx9/uGUzhFSpUlGvXr00QgBmRocABDg61HAP\nCIAACJgJ8FjvunXriANpcAznRo0agQsIOE0AAuw0KmQEARAAgf8R4M3ueS0v79V74cIF0fp1h/jy\nDki7du0SOxtxaxpJ3wSkFWD+Zskvube3t75rAN6BAAhojgCL7bVr18S2gUuWLBHLg2LqBI/5/vDD\nD9S8eXOC+MaUpjbul0KAnz17RnPmzKG9e/dSjx496Pnz52IKP6994+n3kyZNghBr432ClSCgewKD\nBg0i3rOXf9wllLzEqFu3btSkSROqXbu27hnCwf8RkEKAeQsunm5fo0YN6tq1K4WEhIiJDDzzjwV5\n1apV1Lp1a9QZCIAACKhGgBsGvE8v/+ZWr7vE9927d8RLlkaOHInYzqrVrjoFSyHAGzZsEN8mEydO\nTPfu3RMBykuVKiWI8Po3FmEIsDovCEoFAaMT4BjOHEaSo1hlzpxZbAPI/1e5K8WPH58qVqzorsfh\nORoiIIUA58mTh3bs2EEVKlQQC9f5hbekU6dOieDjls/4DQIgAAJKEeCeuVq1alHSpElp69at2ABB\nKfAGKUcKAe7Zsye1adOGrly5IiYhcBcPizLvlbl//34x29Ag9QE3QQAEJCAQHBxM3AvHw2H8e/Hi\nxW61iue3zJo1SwTYqF69ulufjYdph4AUAswv+Llz5+jx48dic2p+Obdt20ZPnz6lhQsXIv6pdt4n\nWAoCmifAEa1y585NxYoVozVr1lC8ePHc6hP38PGwmr+/P0F83YpWcw+TQoCZWqxYsYT48jGHcqtT\npw4fEgci5+VIhQsXFp8j+4dnU3Pr2V7iPyr+NosEAiAAAuEJhIaGiv9n+vbtK/boZQHmiFbuTm/e\nvBFrhvn5HTt2dPfj8TyNEZBGgB1x4xnQvN5u7ty5jrJYz2/ZsoX+/PNP62fbg6NHj1LatGltT+EY\nBEAABEQ0K55kdfr0abHccefOnR4b6+XGxLRp0yh2bGxEh1ePSHoBHjhwoNP1xGvo+Mde6t69O929\ne9feJZwDARAwMAEOIfnXX3/R2rVrxXivJ1GkSJHCk4/HszVGQLqvYdxN/OTJE41hhLkgAAJaJHDo\n0CGaOXMmlSxZ0uPiq0U+sNmzBKQQYF6I3r9/f7FZNa+J42+JvM4uf/78YhKWZxHg6SAAAkYjcPv2\nbRo6dKiIOlWzZk2nhriiw4hD6nIv3rBhw0RXd3SegXv0S0CKLuguXbqI7uFNmzaJSDAsvjyhimdG\nc3g2nrjAu40ggQAIgEBMCPzyyy80ffp0MplM4gv/b7/9RpUqVYrJIx3ey+LLwsuJRRjjvg5RGfaC\nFAK8fft2OnjwYJiA5rzwnZcnTZ06lYYMGQIBNuwrCsdBwD0EOJIVR9YbMWIEde7c2T0PjeQpvM43\ne/bs1KxZM4oTJ04kOXHJqASkEGDuag4ICBDh3sJXxMaNG8nPzy/8aXwGARAAAacJnDhxgsqXL0/t\n2rVTRHzZsE6dOjltHzIak4AUAjx8+HDxLXHy5MmULVs28vHxIY5Ec/78ebF2d/PmzcasHXgNAiAQ\nYwLc3VyvXj3iLQTHjRsX4+fhASDgLgJSCDAH2Th+/Ljohg4KChLjwdzq5XHfcuXKiSAd7nIYzwEB\nEDAOgQULFtCUKVNENKuJEyd63PEbN26IgB45cuTAmK/HaWu/ACkEmDEmTJhQbMagfaTwAARAQAYC\nHPBi9OjRxDEA+vTp43GTeMyXN49hoceEK4/j1kUB0giwLmjCCRAAASkIcEjaCRMmiHCPvXv39rhN\nHKmPe/EmTZpEiRIl8nh5KEAfBCDA+qhHeAECIGBDoFq1aiKWwIABAzzeGuVgHhy/nlvc3JOHBALO\nEoAAO0sK+UAABKQncOTIEbHM6Pr166I7WInlPxxFi3+QQMBVAlJEwnLVaOQHARAAAVsC//77r9hB\nrXbt2mIrQe4OTp48uW0WHIOAdATQApauSmAQCICAKwR4mVGrVq2E8O7fv99jOxnZ2sRR+k6ePElf\nfvklxnxtweDYJQJoAbuEC5lBAARkI8ABNniT+59++kkR8eWQuRy7oGLFihBf2V4GjdmDFrDGKgzm\nggAIfCSwYcMG4ihXN2/epCRJkny84KGjbdu20aJFi8SEq9SpU3uoFDzWKAQgwEapafgJAjojcO/e\nPfrmm29EvHglxPfBgwfEMeqXLFlCXl5eOqMJd9QgAAFWgzrKBAEQiBYB3mHo119/pTFjxoj7K1eu\nLMZ/o/UwF2/i6HyIS+8iNGSPlAAEOFI8uAgCICATgZYtW9KBAwdo8ODBxMdKLDOSyX/Yoi8CEGB9\n1Se8AQHdEXjy5Anxhiw8zrtlyxZavnw5Va1aVRE/z549SwsXLqQffviBMmXKpEiZKMQ4BDAL2jh1\nDU9BQHMEQkJCxPpe3sWIo0ytXLlSMfE9duwY9e/fn7799luIr+beHG0YjBawNuoJVoKA4Qjw+l7e\nTOHy5csiqpWvr69iDHiNLy81GjVqlFhfrFjBKMhQBCDAhqpuOAsC2iAwfvx4Wrt2Ld26dYuWLl1K\nSoovT/TKmTOnKB+7GmnjfdGqlRBgrdYc7AYBnRLgMV5ueXJ3Mwe7iBtX2f+mWHSxzEinL5dkbin7\nZkvmPMwBARCQh8CjR4+E8C5evFiMvVapUkVR47jl++rVKyG+mF2tKHrDFoZJWIatejgOAnIQCA4O\nJt42MF++fMSzjjme848//qiocVevXqWGDRuK+M4QX0XRG7owtIANXf1wHgTUJzBy5Egx3rps2TL6\n4osvFDeIty7s1asXde3alUqXLq14+SjQuAQgwMate3gOAqoTuHPnjljX2717d1XEl8NLzpw5k/r2\n7UvFihVTnQcMMBYBCLCx6hvegoA0BA4fPizW9HKrlwVYjcShJceOHatG0SgTBAhjwHgJQAAEFCVw\n9+5dGjp0KFWrVo369etHq1evRkhJRWsAhclCAC1gWWoCdoCAzgn07NmTdu7cSS9evKAcOXLQH3/8\nQWXKlFHc63fv3hFHucqYMSOlT59e8fJRIAhYCECALSTwGwRAwGMEvvvuO1q/fj3NmDGDmjRp4rFy\nonrw48ePRVznL7/8kkqWLBlVdlwHAY8SgAB7FC8eDgLGJhAYGEhNmzYVEa2OHj1KmTNnVg0IL3fi\nseYaNWpQo0aNVLMDBYOAhQDGgC0k8BsEQMCtBKZPny4mWSVLloyCgoJUFV/udt61a5dYatSsWTO3\n+omHgUB0CUCAo0sO94EACERKYPLkyWLP3u3bt1OCBAkizevpi/Hjx6f69etTkSJFPF0Ung8CThOA\nADuNChlBAASiIsDrankThQoVKhB3+fJ2ftjQICpquG5UAhBgo9Y8/AYBNxNg4S1VqpQIbFGnTh3i\nIBvx4sVzcynOP473EuaNHXbv3u38TcgJAgoSwCQsBWGjKBDQK4FNmzZR27ZtRUhHDuuodpczj/n2\n7t1bjDtjzFevb532/YIAa78O4QEIqEKAlxX9+uuvxJOseILTt99+KzZVUMUYm0JZfDnAR9q0aalH\njx42V3AIAnIRgADLVR+wBgQ0QeDmzZtiglWfPn3EMiMWvCxZskhhO9vG+wknTJhQCntgBAg4IgAB\ndkQG50EABCIQeP36NbVp04YOHjxI9erVk6LFG97IrFmzhj+FzyAgJQEIsJTVAqNAQD4C586dowYN\nGlDixIlpw4YNVLBgQfmMhEUgoCECEGANVRZMBQG1CJQvX55u375NZcuWpV9++UXV2c32GHDQj1ev\nXolJYHHixLGXBedAQDoCEGDpqgQGgYAcBM6fP09z5swRUaxOnTpFJ0+eFBsYyGHdRyvGjBlDHOOZ\nx30hvh+54Eh+AlgHLH8dwUIQUJwAR6+qW7cunTlzhrp160ZXrlyRUnwnTZpEiRIlEuLL0a6QQEBL\nBNAC1lJtwVYQ8DCBR48e0YIFC2j8+PHUokULGjlyJHl5eXm41Og/HsuMos8Od6pPAAKsfh3AAhBQ\nnQDv0Xv48GFq3rw5FStWjGbNmiUmXKluGAwAAR0TgADruHLhGgg4S6BEiRLk7e1NnTt3lnJpka0f\nz549o7dv35Kvry/FihXL9hKOQUBTBCDAmqouGAsC7iVgMpmoXLly9PLlSzp27JjqISSj8m7RokW0\nb98+mjp1KsQ3Kli4Lj0BCLD0VQQDQcBzBObNm0cXLlyg//77T3rx5Y0VAgICaNq0aZQkSRLPQcGT\nQUAhAhBghUCjGBCQhcD79++pUqVKxFsHchduz549RTxnWeyzZ8f+/fvp1q1b9PPPP4uucnt5cA4E\ntEYAAqy1GoO9IBADAkFBQWLThIcPH9KhQ4dES1IL46hlypQh/kECAT0RgADrqTbhCwg4IMBRoniz\nBN4hyN/fnxYuXEg+Pj4OcuM0CICAEgQgwEpQRhkgoCIBnmDVuHFjSpcuHR0/flxFS1wr+saNG2J8\nunTp0iLYhmt3IzcIyE8AkbDkryNYCAIxIvDNN9+ISVZLliyJ0XOUvJn3F+7evTvlzZsX4qskeJSl\nKAFpW8Bv3rwRcV3jxYunKBAUBgJ6InDgwAHasWMHcWhJ7nrWQtq7dy/NnDmTJk+eTOnTp9eCybAR\nBKJFQIoW8PXr18Xm3kePHhUzM9u2bUtp0qQRMzN579F3795FyzncBAJGJcCbE7Dg8vaBgwcPJg60\noYXEOy4FBweLcJiZMmXSgsmwEQSiTUAKAeb/IPiPLV++fMTbioWEhIgg8LwDy/Pnz2nEiBHRdhA3\ngoDRCAwZMoQKFSpE2bNnF1sI8jIjrSQep65Tpw4lTZpUKybDThCINgEpuqA5sg0HA+DdTNatW0fr\n16+nDBkyCKdYfNu3bx9tB3EjCBiNAH+JHTt2LH3//fdGcx3+goCmCEjRAs6ZMydZJojwxt+bN2+2\nQty4cSPlyJHD+hkHIAAC9glwWMn+/fuLL7JaEl9em8xfGHjWMxIIGImAFC1gjm5Tq1Ytmj9/vug2\n69WrlxgDih07NnHgdW4hI4EACDgmwOLbsWNHWrVqFf3555+OM0p2hfcbHjhwIA0bNkzK/YYlwwVz\ndEZACgHOli0bnTt3TszWvHjxohgPTp48uWj51qxZk+LGlcJMnVU93NEDAd4VaOvWrWKiFR9zj1HJ\nkiU14RoPOw0YMIB4zLpgwYKasBlGgoA7CUijbBwOr0qVKuLHnQ7iWSCgVwK8I9Ds2bNFl3OrVq3E\nVoJaWbbHKxu8vLxoxYoV4rde6wh+gUBkBKQRYEdGcouYw+gVLlzYURbr+Tlz5hDvmGIvBQYGilB8\n9q7hHAhoiQB3N/NwzdChQ8X7XqFCBUqYMKGWXBBfGjJnzqwpm2EsCLibgPQCzGNa165do7lz50bp\nO088cTT5hKPq3L17N8pnIAMIyEyAl+Z16tSJbt68SaNHj6bq1avLbC5sAwEQiISAFLOgbe3jNcBP\nnjyxnuIJGs6Ir/UGHICATgnwZMUvvviCihcvLtbJd+jQQVOe8naCLVq0II50hQQCIEAkhQDzeBAv\nn8iYMaPomkqRIgUlTpyY8ufPL3ZtQUWBgNEJXL582bq2d+LEieLvQ0tM7ty5Q926daOvv/6aPv/8\ncy2ZDltBwGMEpBDgLl260NmzZ2nTpk1i2dGHDx9EBB9u+f7yyy80a9YsjwHAg0FAdgKjRo0Se+GW\nK1eOBg0aJLu5Eex78OCBGK9u164dVa1aNcJ1nAABoxKQYgyYA8UfPHhQxH+2VASHoitVqhTxTE9e\npqC17jaLH/gNAtElwF9KeVOC3377jZYtW6bZ8V4/Pz8xWzu6HHAfCOiVgBQtYO5qDggIsMuY1zXy\nHzASCBiBAG9EcPjwYRETnVu8POOZA9FgspURah8+Go2AFC3g4cOHU7NmzcT2YxyUw8fHR+yIcv78\nebExg21oSqNVEPw1BoGdO3dS165dheByHHTeTIGHXho3bqxJAKGhocQhJpMlS0YpU6bUpA8wGgQ8\nTUAKAeY1vsePHxfd0PxHy8uFuNXL3c7cCuAgHUggoFcCvH6dw6/ysqLOnTtr3k0OH8tfJnh9csuW\nLTXvDxwAAU8RkEKA2TkOJMB/sEggYCQCU6ZMEdttckSoGjVqaN71Fy9eUI8ePcSkMYiv5qsTDniY\ngDQC7GE/8XgQkIYAr3UfP348/f7779bN5/UgvryccOXKlWKtL5YaSfO6wRCJCUCAJa4cmKY/Avfv\n3xcTqp4+fSoCzJQtW5a0Er85qtrg/bzbtGkTVTZcBwEQ+H8CUsyCRm2AgN4JHDhwgNKmTUu8uxfv\nf33y5EmqWLGibsRX7/UH/0DAEwTQAvYEVTwTBGwIcAjGunXrUsOGDXUXVIaD5uzatYu8vb3Fun0b\nt3EIAiAQBQG0gKMAhMsgEBMCz58/Fzt5cQxnvUV047Hsvn37inXLWtmDOCZ1iXtBwN0EIMDuJorn\ngcD/E+AgGl9++aWIcc6znPWUeJ3v4MGDRUzqfv36YamgnioXvihGAF3QiqFGQUYhwGLL8Zt5QxHe\nSpPHe/WWTp8+LdYu88YpSCAAAtEjAAGOHjfcBQIRCHCrkIV33rx5IgDFyJEjI+TRywmO1IUEAiAQ\nMwIQ4Jjxw90gYCXAovTq1SuxhWb58uWt53EAAiAAAvYIOBRgXlTP4SF5R5YrV65Q9uzZ6dNPP6U8\nefJQ3LgOb7NXBs6BgC4J8Bgvh12cNGkSLV68mBIkSECBgYG6HQ/lrvWHDx9S+/bt8X+ALt9oOKU0\ngQiTsPg/Fd7+LF++fNS9e3favXs3vX79mnhDhObNmxMHiucNwVmgkUDAqAQeP35M9erVo4IFC4qd\nvHi7wIsXL+pWfDlkJq9l/vbbbyG+Rn3p4bfbCYRpyr5584YaNGhAlStXFhsj+Pr6RijwyZMnYo9S\nDiLAYefSpUsXIQ9OgICeCfB2gbxpAm8SwsepUqXSs7s0Y8YMunTpkvjizTHbkUAABNxDIIwAc9cy\ndzMlSZLE4dOTJ09OAwYMEAHXeRE+EggYicCpU6dEUA1u/fLfgd7Fl8e09bBDk5HeUfiqHQJhuqBZ\ngC3i++eff4rxHltXLl++LL4F8zkvLy+xzML2Oo5BQK8ENmzYQFmzZqUmTZqIXqKpU6caovcnUaJE\neq1S+AUCqhMII8C21ty7d09E8NmzZ484PX/+fCpevLhVoG3z4hgE9Exg+PDh1KpVKxo2bBidOXOG\npk+frvtxUJ4LwsuqkEAABDxHIEwXtG0xPNkic+bMYj0jB5HnP0YWY39/f9tsOAYB3RLgyYe5cuUS\nvT1btmwho4Rb5Lkdf/zxh/iiwUNOSCAAAp4h4FCALcXxVmk8OYu3Gosd22GD2ZIdv0FAFwR4uGXI\nkCFiadGFCxd0O7s5fGWtW7eO1q9fD/ENDwafQcADBBwq6s8//yw21uaxLg6lx2v/eJNtjvKDBAJ6\nJ1CuXDniSYYsSDzb2QiJe7h4qRHPek6ZMqURXIaPIKAqAYct4EyZMhHP+LTM8mzbti1xdB/+dowE\nAnolwHMfypYtK3p8li5dSnHixNGrqxH84r9v/kECARBQhkCYFjCPeQUEBIiSa9eubRVfiynZsmWj\nnj17io9//fUX3blzx3IJv0FA8wQ4ylOlSpWIv3xevXrVUOKr+cqDAyCgQQJhBJjHebkbitc4cjSs\n4ODgMC7dvHmT1qxZQ9w9x9ctS5bCZMIHENAggZkzZ1KBAgXE0iKegGSUxBG9jh07JqLdGcVn+AkC\nshAI0wXN3W281OL27dtiAkqXLl3EOFjq1KmJxTdZsmRUqlQp4vFhzIaWpQphR0wJ8JfJPn360KJF\ni8Qa35g+Tyv379+/X8Sx5tCyvK4fCQRAQFkCYQTYUjSHl5w7d674uXv3rnUzBst4sCUffoOAlgk0\na9aMeD9bFuBp06YZSnz/+ecfmjBhgvjJkiWLlqsRtoOAZgnYFWBbb9KkSUP8gwQCeiLA3c3c48MT\nrXr16kWffPKJntyL1JcbN26IgCJz5syJMM8j0htxEQRAwK0EHArw06dPqWPHjnT69OkwOx9Vr16d\neGcUJBDQIgGeaMgTDLlnh4dVeN6D0VLGjBmJVzUggQAIqEvAoQD/9NNPYhIWd815e3tbreQuOyQQ\n0CqBbt26EQfZ4H2ujSi+Wq032A0CeiTgUIBv3bolWsAVKlTQo9/wyWAEOJRqzpw5RUjVtWvXkp+f\nn6EI3L9/n3iDlWrVqlH69OkN5TucBQFZCYRZhmRrZP369cX4GP/hIoGA1gnw0rnEiRNTUFAQFSlS\nROvuuGQ/7+Xbrl07sXIB4usSOmQGAY8ScCjAvBRp8+bNxBsx5MiRg3Lnzi1+uAsPCQS0ROCLL74Q\nS+uMtL7XUj9XrlwRS6x+/PFHsZuZ5Tx+gwAIqE/AYRd0rVq1qGjRohEsxBhwBCQ4ISkBjmHet29f\nEcucg01whCsjJZ5wxoE2eCtR/N0aqebhq1YIRBBgHiNavnw58UxJX19f4vB8fIwEAlohcP78eerU\nqZNYatO6dWvxPhtxWz0OrmHvS7RW6hF2goDeCUQQ4KNHj9L79++F34cPH6aBAwcSx31GAgHZCfA4\n55EjR4SZPN7L2wii5Sd7rcE+EDAugQgCbFwU8FzLBJYtWyYiWnHr18gTjbjLedSoUVS3bl2xq5OW\n6xS2g4DeCTichKV3x+Gffgjs3LmTevToIXprjCy+Dx48EF3vPOObt1REAgEQkJuA3RYwRwh68+aN\niBb09u1bunbtmtUL7trjsWEkEJCBwI4dO6hx48ZCfI08Q59bvux/gwYN6Msvv5ShamADCIBAFATs\nCnD4iRu2cXIbNWpEK1eujOKxuAwCnifAW2c2bdpUtH4t+1R7vlQ5S/Dx8RHr9mPHRqeWnDUEq0Ag\nIoEIAnzv3r2IuWzOxIoVy+YTDkFAeQLPnj2jBQsW0NSpU4ljkw8aNEh5IyQrMW7cCH/KklkIc0AA\nBMITiPBXyzvEIIGArATOnTtHHFgje/bsNGvWLKpYsaKspipiV3BwsIhpjf18FcGNQkDArQTQX+VW\nnHiYpwhwSFSe2fv5558Tr+3lpXG8Zt2oGyq8evVKhJdcsWIFQXw99dbhuSDgWQIRWsCeLQ5PBwHX\nCXAIyQEDBlCyZMnELkapUqVy/SE6uoMjXPEexv7+/kKEdeQaXAEBQxGAABuqurXl7PPnz+nTTz8l\nnnfQsmVLMcuXZ+EbOXGQnDlz5oiu94YNGxoZBXwHAc0TgABrvgr16UDhwoWJJxZxa3ffvn2E2b3/\nq+d48eJR165d9Vnp8AoEDEYAAmywCpfd3W3bttHp06fp8uXLxLOdkUAABEBArwQwCUuvNatBvzj4\nC68z5y7ndevWadADz5l8/PhxOnXqlOcKwJNBAAQUJ4AWsOLIUaA9AkFBQVSjRg1KkyYNGT2ohi2f\n0NBQGjZsGJlMJho6dKjtJRyDAAhonABawBqvQD2Yz5sH8P7TmTNnJt5MAel/BD58+CA2VggJCaEh\nQ4YQ1ujjzQABfRFAC1hf9ak5b1hYpk2bJvbs5aAaEJmPVfj333+LNc+8wQQmoX3kgiMQ0AsBaQWY\nxwPfvXtHHOMWSX8EeIJV27ZtieM5b9myhUqWLKk/J2PoUZkyZWL4BNwOAiAgMwFpu6DXrFkjguzL\nDA+2RY8Aiy8vM+Jdtw4dOgTxjR5G3AUCIKBxAlK0gHPkyEEPHz4Mg5Jbvzz2xULMIQgXLlwY5jo+\naJMATyriABLJkyenAwcOoGs1XDXy9op37twRuzwZNcxmOCT4CAK6JSCFALO4tmnThr7++mv65ptv\nBOz169fTwYMHady4cWT06Ed6eftYfPPnzy/2mmahwbhm2JqdPXs2HTt2jCZNmmTYGNdhieATCOib\ngBRd0DzWdfToUQoMDBTdziy4vr6+5O3tLWbG8jGS9gksXbqUnjx5QlevXhW7GWnfI/d5wNsrHjly\nhCZMmIAvnO7DiieBgNQEpGgBMyGebLVkyRJauXIllStXjkqUKIEZsVK/Os4bt2rVKho8eLAQX/6N\nFJYAD7/wO8+9P5gFHpYNPoGAnglI0QK2Bdy4cWPavn27GBPmoAxI2iXAM9l5EwXeuadevXp09+5d\n6tixo3Yd8pDl3MPD+xtDfD0EGI8FAUkJSNMCtuWTIUMG+vPPP8WpixcvEu99yrNmo0rchXfixAm7\n2Ti+cIIECexew0n3E+BWXYECBShTpkzE47080Q4JBEAABEDgIwEpBfijeUTcfXnt2jWaO3eu7Wm7\nx7xTjKMJW3yNYwwjeYYAb5PH+/byEqN///1XxHIuWrSoOOeZErX9VO7lWbt2rYh0lTJlSm07A+tB\nAASiRUB6AR44cKDTjhUqVIj4x17i1jF3gSK5nwCHTOQeCl421qNHD9Hy7dSpE+XJk8f9hengiVu3\nbqX58+fT9OnTCeKrgwqFCyAQTQLSCTD/J84bsfM6USS5CXBd3bt3j8aMGUOPHj2i27dvo5chiirb\nvXs3/frrryL8JuY4RAELl0FA5wSkmITFQTf69+9PGTNmFOsfU6RIIbqSec0oAnDI+wa2a9eOypYt\nSxcuXBAT59DFH3VdcbxrFuC0adNGnRk5QAAEdE1AihZwly5dRPfwpk2bKGvWrEJ8eSzx3Llz1K1b\nNxG4oUOHDrquCK05N3HiRDG+GxAQIIJraM1+2AsCIAACahOQogXME1I4ChDPmuXgG9ySSpo0KZUq\nVYqmTp1KHBULSR4C3I3KWwjOmDED4utEtbx+/ZquX78uNhdxIjuygAAIGISAFALMXc3ckrKXNm7c\nSH5+fvYu4ZwKBE6ePEnNmzen7777jpo0aaKCBdoqkiO8ffXVV8QijNjO2qo7WAsCniYgRRf08OHD\nqVmzZjR58mTKli2biIoVHBwsNmfniT6bN2/2NAc8PwoCHJebeyl4Bm+DBg3ExKsobjH8ZV6TPnLk\nSBo9ejTlypXL8DwAAARAICwBKQSYl7AcP35cbL4QFBQkxoO51cvjvhyiD5N7wlaa0p84fnOdOnWo\nevXqoqcCy4uirgHucuZ9jqdNmyaCkUR9B3KAAAgYjYAUAszQEyZMSBUqVDAaf2n9PX/+PM2cOZN4\nByMe8y1WrBgtWrQIOxg5WWMcAaxfv35O5kY2EAABIxKQYgzYiOBl9pn3oy1durSYNPT999+LPZl5\nLB7bB8pca7ANBEBAawQgwFqrMQ/bO2/ePOIQktz1z2O+HFksX758EF8nuPO8Bf6igohrTsBCFhAA\nAYIA4yWwEmjVqhX9+OOPNGvWLCz9slJx7oDnLrRu3ZqSJUtGiHDlHDPkAgGjE5BmDNjoFaGm/yaT\niZYvX04cCIWXGXFEMiTnCdy4cYN69uxJHFCmTJkyzt+InCAAAoYmAAE2dPWT2L2I1/Tu3buXevfu\nDfF18X14+fIlnT17lqZMmQJ2LrJDdhAwOgEIsIHfAN5zmXcv4v2XDxw4INZgGxhHtFzn7S+rVasW\nrXtxEwiAgLEJYAzYoPXP61R5zLdFixa0Y8cOiK9B3wO4DQIgoB4BCLB67FUr+c2bN9SwYUMRe3vw\n4MEUNy46QlypDN4uc/z48fT333+7chvyggAIgEAYAhDgMDj0/2HJkiUiqAYH2OBjJNcIPH36lDp2\n7CiiW3322Weu3YzcIAACIGBDAE0fGxh6Pfzw4YPY2nHs2LEilCTHJ+YNAjj6GJLzBHidb/fu3alK\nlSrYiMJ5bMgJAiDggAAE2AEYvZxm8f3222/FOG/NmjXFxhYFCxbUi3uK+cFLtbjrngOVxIkTR7Fy\nURAIgIB+CUCA9Vu3wjPes5d3k9q5cyf27o1BXfOGIKlTp47BE3ArCIAACIQlAAEOy0M3n06fPk3d\nunWjy5cvi5nOvOcyEgiAAAiAgDwEIMDy1IVbLeElRilTpqRTp06J/ZXd+nCDPOzdu3ciNGeWLFmo\na9euBvEaboIACChFAAKsFGkFy5kzZw4FBgbSvn37iANFILlOgMW3T58+lD59eoiv6/hwBwiAgBME\nIMBOQNJSFl5axBsqLFy4EOIbzYoLCQkR63yzZctGnTt3juZTcBsIgAAIRE4AAhw5H01dvX37NvXq\n1YtGjBhB9evX15TtMhnLgUkGDBggk0mwBQRAQIcEEIhDJ5V669YtsROPv7+/2JVHJ27BDRAAARDQ\nLQEIsA6q9syZM8RRmXLlykXbtm3TgUfquMDxsa9evapO4SgVBEDAcAQgwBqu8hcvXtDUqVOpfPny\nVKdOHdqyZQviOkejPjnIxqhRo4jXTPPOUEggAAIgoAQBjAErQdlDZTRo0IDu3LlD06dPF6ElPVSM\n7h/LGys8fvyYxowZQ/HixdO9v3AQBEBADgIQYDnqwSUrOCZxhQoV6P79+/TXX38Rr1NFih4B3oqx\nVq1alDNnTvQeRA8h7gIBEIgmAQhwNMGpedunn35KadKkoaCgIIhGDCuicuXKMXwCbgcBEACB6BGA\nAEePm2p3rVy5kp48eUKXLl2i2LExhK9aRaBgEAABEIghAfwPHkOASt7O8Z15nLJu3boQ3xiA//ff\nf2nDhg3E0a6QQAAEQEAtAhBgtci7WC5vhVexYkUqVaoUDRw40MW7kd1CYOnSpWK2M88cjx8/vuU0\nfoMACICA4gTQBa04ctcLfPToEVWvXp1SpEhBM2fOdP0BuEMQWLFihdgXedq0aZQ0aVJQAQEQAAFV\nCaAFrCr+qAvfv38/5c2bl5InT06HDh2K+gbksEvg5s2blDlzZpo3bx4lS5bMbh6cBAEQAAElCaAF\nrCRtF8viscomTZpQmzZtxNivi7cjuw0BDrCBIBs2QHAIAiCgOgG0gFWvgogG8G48U6ZMoZo1a4o1\nqhylCQkEQAAEQEBfBCDAktXnyZMnKUeOHGI7Qe4unT17NmY8R7OODh8+TP3796cHDx5E8wm4DQRA\nAAQ8RwAC7Dm2Lj95woQJVKlSJapXrx4dO3ZMtH5dfghuEAQCAgJo3Lhx1KFDB/Lz8wMVEAABEJCO\nAMaAJakSHuvdtWsXLV68WHQ9S2KWJs3giWscH3vSpEmUMWNGTfoAo0EABPRPAAIsSR1v3bqV/vvv\nP0qVKpUkFmnTDB4/L1iwIK1evRpd99qsQlgNAoYhAAGWpKrjxIkD8XVDXcSNG5eSJEnihifhESAA\nAiDgWQIYA/YsX6efHitWLKfzImNEAh8+fKDXr18T/0YCARAAAS0QgABLUksQjuhXBMfI5olrgYGB\n6HaOPkbcCQIgoDABCLDCwB0Vx13QSK4TOHfuHA0aNIiGDRtG/v7+rj8Ad4AACICASgQwBqwS+PDF\nhoaGhj+Fz1EQuHbtmggtOXr0aBGuM4rsuAwCIAACUhGAAEtSHWgBu14RHNuZlxohgQAIgIAWCaAL\nWpJaQwtYkoqAGSAAAiCgEAEIsEKgoyomdmxURVSM+DrPdOYQkw8fPnQmO/KAAAiAgLQE8L++JFWD\nWdBRV8StW7eoRYsW9OzZM/L19Y36BuQAARAAAYkJQIAlqRy0gCOviLt371L37t3p22+/FfGyI8+N\nqyAAAiAgPwFMwpKkjjAG7LgiXr58Sbt376ahQ4ditrNjTLgCAiCgMQIQYEkqDLOgHVdE4sSJqVmz\nZo4z4AoIgAAIaJCA9F3Q3DJ8+/atBtG6ZrLJZHLtBuQGARAAARDQNAEpBPjGjRvUsmVL8vb2psqV\nK4uQghaqq1atEhNvLJ/1+htjwGFr9s2bN7RkyRIx4znsFXwCARAAAX0QkEKAJ0+eTGnTpqWjR49S\nqVKlqFy5cnTp0iV9EHbSC8yC/gjq+fPn1LlzZ+Lej+LFi3+8gCMQAAEQ0BEBKcaAN2/eTMePHycv\nLy8aPny4mGhTtWpV4o3VjZKwG9L/aponXPXs2VMIb+vWrY1S/fATBEDAgAThic0dAAAXX0lEQVSk\nEOC8efOK1m/ZsmVFFTRt2pRu375N1atXp3bt2hmiWtACJrGVIA9HTJkyhRIlSmSIeoeTIAACxiUg\nRRd0+/btqVGjRjRu3DhrTfTo0YMaNGgg1n5aT+r4AC1gElsJ5s6dG+Kr4/ccroEACHwkIEULuEqV\nKnT58mW6cuXKR8vMR0OGDKHPP/9cXAtzAR9AAARAAARAQOMEpBBgZshrPe3t58qTs5ImTeoU5vfv\n3xP/2Et8XuZuXqO2gENCQmjChAmUPHlywww32Hs/cQ4EQMB4BKQRYEfoeRkS7/s6d+5cR1ms55ct\nW0ac317ijdszZcpk75IU54y4DpjFd9CgQaLL+bvvvpOiHmAECIAACChFQHoBHjhwoNMsWrVqRfxj\nL3EcYY4njCQHAV5ixKEl48ePT/379xfjv3JYBitAAARAQBkC0gkwt4p4HSh3SRopGa0LmkNvjhw5\n0khVDF9BAARAIAwBKWZBv3v3TrSCMmbMKFpEKVKkEGPC+fPnp4ULF4YxWK8fZB6f1itz+AUCIAAC\nahKQogXcpUsX0T28adMmypo1qxBf3vOVx227detGHJawQ4cOanLyeNlG2YwhODhYRLjiL1lIIAAC\nIGBkAlK0gLdv306zZ8+mAgUKiHjQ3B3LM585LOXUqVNp/fr1uq8jI2xHOHHiRLG0jGe8I4EACICA\n0QlIIcDc1RwQEGC3LjZu3Eh+fn52r+nppN5bwNOmTaOrV6/SmDFjKEGCBHqqOvgCAiAAAtEiIEUX\nNMd/5v1eeVOGbNmykY+PD3FX5fnz54knZXGsaL0nPbeA+UtUiRIlxDpfiK/e32T4BwIg4CwBKQS4\ncOHCYjOGgwcPUlBQkBgP5lYvj/vyzkhGmCGs5+0Ia9Wq5ez7iHwgAAIgYBgCUggw006YMCFVqFDB\nMODDO4pZ0OGJ4DMIgAAI6JuAFGPA+kbsnHd6awFzXO+9e/cSLzFDAgEQAAEQiEgAAhyRiSpn9NQC\nXr16tQgxybG9OdIVEgiAAAiAQEQC0nRBRzTNWGf00gLesGEDrVmzhnjWM9b6GusdhrcgAAKuEYAA\nu8bLY7n1sBkDdzvzhLn58+djT1+PvSl4MAiAgF4IQIAlqUk9zPTmKGb8gwQCIAACIBA1AYwBR81I\nkRx6aAErAgqFgAAIgIBOCECAJalIrbaAL1y4QBxi8uHDh5KQhBkgAAIgoA0CEGBJ6kmLLWAOnNKv\nXz+qV68e+fr6SkISZoAACICANghAgCWpJ621gA8fPkw//fQTjRs3DuO+krxDMAMEQEBbBDAJS5L6\n0lIL+PXr16LFu2LFCmysIMn7AzNAAAS0RwACLEmdaakF7OXlhVavJO8NzAABENAuAXRBS1J3WmoB\nS4IMZoAACICApglAgDVdfcoZHxgYSC1atKATJ04oVyhKAgEQAAEdE4AAS1K5MndBs/j26tWLOnfu\nTIUKFZKEGMwAARAAAW0TgABLUn+ybsZw/fp1GjFihFhuVKJECUlowQwQAAEQ0D4BTMKSpA5l3Ywh\nU6ZMtHjxYkkowQwQAAEQ0A8BtIAlqUtZW8CS4IEZIAACIKA7AhBgSapUphZwSEgI8c5Gz549k4QO\nzAABEAAB/RGAAEtSp6GhoVJYcv/+ffr666/p7Nmz5OPjI4VNMAIEQAAE9EgAAixJrcaJE0d1Sx49\nekRdu3alJk2aUO3atVW3BwaAAAiAgJ4JYBKWJLWr9hjwy5cviUNLsgCXLFlSEiowAwRAAAT0SwAC\nLEndqj0GnDhxYurUqZMkNGAGCIAACOifALqgJaljtVvAkmCAGSAAAiBgGAIQYEmqWo0W8Pv372nr\n1q108uRJSSjADBAAARAwDgEIsCR1rfRmDLyl4A8//CCWGxUsWFASCjADBEAABIxDAAIsSV0rGQv6\nzZs31Lt3b8qdOzd17NhREgIwAwRAAASMRQCTsCSpb6VawDzWzDsaDRkyhPz8/CTxHmaAAAiAgPEI\nQIAlqXOlWsA81oxlRpJUOswAARAwNAF0QUtS/Uq1gCVxF2aAAAiAgOEJQIAleQU82QJmcV+wYAEt\nWbKEIPSSVDjMAAEQMDwBCLAkr4CnhJHHfIcOHSpmOzdv3pw8KfSSoIQZIAACIKAJAhgDlqSaPCGM\nLOqjR48mnvU8atQokiHetCS4YQYIgAAIqE4AAqx6FfzPAE+0gHk7wQEDBqDVK0kdwwwQAAEQsCUA\nAbalobPjpEmT6swjuAMCIAAC+iGAMWBJ6tITXdCSuAYzQAAEQAAE7BCAANuBosYpd3VBz5w5U0S3\nevXqlRpuoEwQAAEQAAEnCaAL2klQns7mjhbw3LlzRZSrSZMmUaJEiTxtMp4PAiAAAiAQAwJoAccA\nnjtvjel2hOvWrRNLjVh8vb293WkangUCIAACIOABAmgBewBqdB4Z0+0I69WrR/yDBAIgAAIgoA0C\naAFLUk8xbQFL4gbMAAEQAAEQcJIABNhJUJ7OFp0W8IMHD+j06dP0/v17T5uH54MACIAACLiZAATY\nzUCj+7jQ0FCXbt20aRN16dKFUqVKRfHixXPpXmQGARAAARBQnwDGgNWvA2GBK2Eit23bRosWLaJp\n06ZR6tSpJfEAZoAACIAACLhCAALsCi0P5nV2DPjixYt048YN4iVHyZIl86BFeDQIgAAIgIAnCUCA\nPUnXhWc72wLOlSsX8Q8SCIAACICAtglgDFiS+nO2BSyJuTADBEAABEAghgQgwDEE6K7bI4uExV3O\nixcvpocPH7qrODwHBEAABEBAZQLSCXBISAg9efJEZSzKF+8oFvSxY8eoc+fOVKJECfL19VXeMJQI\nAiAAAiDgEQJSCPC7d++of//+lDFjRoofPz6lSJGCEidOTPnz56eFCxd6xHEtPPTkyZM0fPhwGjVq\nFOXOnVsLJsNGEAABEAABJwlIMQmL17PevXuXeG1r1qxZhfjyZvLnzp2jbt260Zs3b6hDhw5OuqSP\nbM+fPxd+c9cz9vXVR53CCxAAARCwJSBFC3j79u00e/ZsKlCggNhIgMdDWXRKlSpFU6dOpfXr19va\nrMvj8GPASZIkEd3OEF9dVjecAgEQAAGSQoC5qzkgIMBudWzcuJH8/PzsXtPTSUdjwHryEb6AAAiA\nAAh8JCBFFzSPczZr1owmT55M2bJlIx8fHwoODqbz588TT8ravHnzR4t1fHTr1i3i7QRbtGhBhQoV\n0rGncA0EQAAEQEAKAS5cuDAdP36cDh48SEFBQWI8mFu9PO5brlw5Ct89q8dq43XAPBbetWtXiK8e\nKxg+gQAIgEA4AlIIMNuUMGFCqlChQjjziDj04qtXr4hFOqq0dOlSWrt2rd1sp06doi+++IJGjhxp\nvc5LfGy7frkM212JlLr++vVrMeGK7UuQIAEdOnRI2FisWDGyjZB1+PBhsg3Ygevgg/cjjvXvGX8f\n+P8hJv8/Wl8kBQ9imQXIpGB5LhfFgnnt2jUR+ziqm3m2NP/YS6tWrRLri5s3b269zEudbBMLvS0O\npa6zzfwFJG7cuGHK54lYtunFixe4bvO6gg/eD/x9fCSA/x+i9/9jz549iXWhSJEiH2EqdCS9ALuL\nw++//y4EuH379u56JJ4DAiAAAiCgcQJqCrAUs6Bt688okbC4pc1rfTkICRIIgAAIgIDxCEghwEaL\nhPX48WMx03nr1q0i8pfxXjt4DAIgAAIgIMUkLCNFwuLlVd27d6caNWpQo0aN8AaCAAiAAAgYlIAU\nLWCjRMLiSV4zZsygJk2aiHXPBn3n4DYIgAAIgICZgBQCbJRIWIkSJaIBAwaI1i/ePhAAARAAAWMT\nkKILGpGwjP0SwnsQAAEQMCIBKQRYz5GweGE4R/niNas5c+Y04jsGn0EABEAABOwQkEKA2S5HkbDs\n2KyZUzy7u3fv3pQ5c2bq0aOHZuyGoSAAAiAAAp4nIMUYsOfdVL4EFt9+/fpR2rRpIb7K40eJIAAC\nICA9AWlawNKTctHAPXv2iM0ksmfP7uKdyA4CIAACIGAEAoYJRXnixAmqWbOmU5s6qFHxu3fvJi8v\nLzWKVrVMjnzGPzwEYbTEm3DEjx8/zIYbRmHw8uVLCh9r3Qi+h4aG0vv37w35vnPMe39/f/L19ZWq\nqq9cuUI7duyg9OnTK26XYQRYcbIuFli+fHniVrPR0r59+2jXrl00bNgwo7kutp9s164d8TI8oyWj\nvu/79++nbdu20YgRI4xW5dStWzdq3bo1FSxY0HC+O3IYY8COyOA8CIAACIAACHiQAATYg3DxaBAA\nARAAARBwRAAC7IgMzoMACIAACICABwlAgD0IF48GARAAARAAAUcEIMCOyOA8CIAACIAACHiQAATY\ng3DxaBAAARAAARBwRADLkByRUfj8nTt3RNQshYtVvTheG8jrYZMnT666LUob8OjRI/Lx8aF48eIp\nXbTq5eF9N+b7zjHxee070v8IQIDxJoAACIAACICACgTQBa0CdBQJAiAAAiAAAhBgvAMgAAIgAAIg\noAIBCLAK0FEkCIAACIAACECA8Q6AAAiAAAiAgAoEIMAqQEeRIAACIAACIAABxjsAAiAAAiAAAioQ\ngACrAB1FggAIgAAIgAAEGO8ACIAACICA2wm8f//e7c/U2wMhwArW6J49e6hMmTKUJUsWqlevHj15\n8sRu6c7ms3uzpCfHjBlDBQoUEL7zsb3Ef7C9e/emokWLip9+/frRu3fv7GXVzDmu48aNG1OOHDnI\n39+f/v7770htf/bsGWXOnJl27twZaT4tXHT2Pf7nn39EfefLl49q1apF58+f14J7kdrozPvODxg2\nbBgVK1aMSpYsSRMmTIj0mVq6+Ntvv1GpUqUcmuzq34XDB2n9gglJEQIPHjwwpU2b1nTy5EmTWVRM\n3bt3N7Vu3TpC2c7mi3CjxCdWrlxpKl26tOnp06cmcwhCU8GCBU2bN2+OYPHcuXNN5i8mgg8zqlOn\njonPaTk1atTINGLECNOHDx9MAQEBptSpU5tevXrl0CV+J5IlS2basWOHwzxauODse2wORWrKmjWr\n6eDBg8It83/cpgYNGmjBRYc2Ovu+b9u2zWQWX/G+M4fcuXNbOTh8uOQXHj9+bOrUqZPJz8/PVKRI\nEYfWuvp34fBBGr+AFrBC36COHj1KefLkEa1Ajv3bpUsXWrt2bYTSnc0X4UaJT2zdupW+/vprSpo0\nKaVJk4a++uorWrduXQSLzcJM48ePF7GRmVHevHnpwIEDEfJp6QT73rFjR4oVKxaVL1+eMmTIQPv3\n77frwvr16ylhwoSUM2dOu9e1dNLZ99j8RYyyZ88uWoDBwcHUtGlTWr16tZZcjWCrs+/7w4cPKXbs\n2OJ9T5AggYiRfPv27QjP09KJXbt2UaJEiWjx4sWRmu3K30WkD9L4RQiwQhV4/fr1MJstmFtCxP/h\nvH37NowFzuYLc5PkH8L7xCJ87969CFZzV1y2bNnE+ZcvX9Ly5ctFl2SEjBo5wd1sXL8pUqSwWsy+\n379/3/rZcsDnzC1lGjt2rOWUpn+Hr3NH7/u1a9cEn3LlypG51STq/+zZs7ry3dH7XrduXUqfPj2V\nLVuWPvvsMzK3gKlmzZqa9r1hw4b0008/kZeXl0M/XPm7cPgQnVyAACtUkbzzTeLEia2lWV5Qc3ek\n9RwfOJsvzE2SfwjvE39DZoF1lHjcl1tCLMjm7khH2aQ/H95vNpjr/cWLFxFsb9eunRgP5N2R9JDC\n++7ofTd3VZO5y5bat28v3v1q1arRuHHjNI0gvO+O3vfAwEC6dOmS6OnhOQJ8fOPGDU377ozx4fnw\nPY7+Lpx5npbzQIAVqj1fX1/iCTaW9Pz5c9HdGH4bPmfzWZ6jhd/hfWIO6dKls2s6i2/9+vUpNDRU\ntIDtZtLIyfB+s9n2fF+xYgVdvnxZeLVx40Yyj5XToUOHKCgoSJzT4j/hfXf0vpvHu4knXzVr1ox4\nq7o+ffrQH3/8oenJd+F9t1fnXKeTJk2i2rVr0+zZs0WXLbeE+ZzeU3g+7K8jRnpnAQFWqIZ57M/2\nP1Q+zpgxY4TSnc0X4UaJT7BP3NVoSY58DwkJES1fFl8eH9f6vqEsLvzN/ubNmxbXxTuQKVMm62c+\nYHHy9vam0aNHix8eB+RZpDyOqtXk7HvM+Vh4LYnH/nl/aPOkNcspzf129n3nbnqe/WxJ3OPDfxt6\nT87+Xeidg/BP45PINGM+z3JMlSqVyby8xMTHLVu2NPXt21fYzzMHLbNeI8unGWfDGbplyxaTeQmS\n6datW6arV6+azJNuTEeOHBG5zp07Zzp9+rQ4njx5sshnFiCTuZtK/JjFKdzTtPWxTZs2JvOEO5N5\niZXJPLlIzHTlGd6ceFY0zxYOn4oXL259H8Jf08rnyN5j2/fd3PIxpUyZ0mReiiRcGzlypMk8HqwV\nN+3a6ez7PmfOHJN5zFTMkOeZ8VWrVjWZJy/ZfabWTvK7HX4WtO37HtnfhdZ8jYm9FJObca9rBHh5\ngrmlYzJPvDBVqFDBZBGXvXv3mszjRNaHOcpnzaCxA16CY1leY56QYhoyZIjVg86dO5tatWolPmfO\nnNlk/lYY5qdGjRrWvFo84C8c+fPnN5m73E3mCWZCdC1+8JKkTZs2WT5af+tBgNkZR+9x+Pd9w4YN\nJn4vmE+uXLlMV65csbLQ4oGz7zt/+TCPfYsvneYxYJN5trzJPDdCiy5HsNmeANu+75H9XUR4mI5P\nxGLfDNHUl8RJ7mblLsfwY7/hzXM2X/j7ZP7M4zy83IJ/jJZ4shHP8jVacvY95v+GeFmOnhg5+77z\nTHlepqb1IZfovNtG/buwsIIAW0jgNwiAAAiAAAgoSACTsBSEjaJAAARAAARAwEIAAmwhgd8gAAIg\nAAIgoCABCLCCsFEUCIAACIAACFgIQIAtJPAbBEAABEAABBQkAAFWEDaKAgEQAAEQAAELAQiwhQR+\ngwAIgAAIgICCBCDACsJGUSAAAiAAAiBgIQABtpDAbxAAARAAARBQkAAEWEHYKAoEQAAEQAAELAQg\nwBYS+A0CIAACIAACChKAACsIG0WBAAiAAAiAgIUABNhCAr9BAARAAARAQEECEGAFYaMoEAABEAAB\nELAQgABbSOA3CIAACIAACChIAAKsIGwUBQIgAAIgAAIWAhBgCwn8BgEQAAEQAAEFCUCAFYSNokAA\nBEAABEDAQgACbCGB3yAAAiAAAiCgIAEIsIKwURQIgAAIgAAIWAhAgC0k8BsEQAAEQAAEFCQAAVYQ\nNooCATUJXL16lfz9/SkwMFCYsXDhQmrUqBGZTCY1zULZIGBYArHMf3z46zNs9cNxoxHo0aMH/fff\nfzR79mwqUKAAbd26lYoWLWo0DPAXBKQgAAGWohpgBAgoQ+Dly5eUL18+8vHxoZo1a9KYMWOUKRil\ngAAIRCCALugISHACBPRLIHHixNShQwc6c+YMde7cWb+OwjMQ0AABtIA1UEkwEQTcReDp06eUN29e\n8ZM2bVpaunSpux6N54AACLhIAC1gF4EhOwhomUDPnj2pWrVqtGbNGtq5c6cYA9ayP7AdBLRMIK6W\njYftIAACzhPYvXs3bdiwgS5cuEBJkyaliRMnUvv27UV3tLe3t/MPQk4QAAG3EEAXtFsw4iEgAAIg\nAAIg4BoBdEG7xgu5QQAEQAAEQMAtBCDAbsGIh4AACIAACICAawQgwK7xQm4QAAEQAAEQcAsBCLBb\nMOIhIAACIAACIOAaAQiwa7yQGwRAAARAAATcQgAC7BaMeAgIgAAIgAAIuEYAAuwaL+QGARAAARAA\nAbcQgAC7BSMeAgIgAAIgAAKuEYAAu8YLuUEABEAABEDALQQgwG7BiIeAAAiAAAiAgGsEIMCu8UJu\nEAABEAABEHALAQiwWzDiISAAAiAAAiDgGgEIsGu8kBsEQAAEQAAE3EIAAuwWjHgICIAACIAACLhG\nAALsGi/kBgEQAAEQAAG3EIAAuwUjHgICIAACIAACrhGAALvGC7lBAARAAARAwC0E/g8EenPNd+NC\nFAAAAABJRU5ErkJggg==\n"
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%%R -i P\n",
"plot(ecdf(P))\n",
"abline(c(0,1), lty=2)\n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0.42954652, 0.43039195, 0.42918563, -0.00043435, 0.00073129])"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.mean(beta_star, 0)[:5]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## CLT seems to hold but ROSI doesn't work"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.axes._subplots.AxesSubplot at 0x1228ef1d0>"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0XNWdJ/Dvr1YtVVqrtFi2LFuSd8sLwg62McEQAgZC\naEjaGEKHZALJJBkyp9MTcjJJJ5mZTmemOxPICkNCQgdCEmwDNsaAAQNmMZYXWZZkbFm2ZO2y9tJS\n650/JIMxkvUkV9WrV/X9nKNjLc+q7ztV9T1P9913nyilQERExmHSOwAREU0Ni5uIyGBY3EREBsPi\nJiIyGBY3EZHBsLiJiAyGxU1EZDAsbiIig9FU3CKSISJPi8gxEakVkSsiHYyIiMZn0bjdgwB2KaVu\nFxEbgJSLbexyuVRRUdGlZiMiShgHDhw4q5Rya9l20uIWkXQA6wF8EQCUUj4Avov9n6KiIlRUVGh5\nfCIiAiAiDVq31TJUMgdAJ4DHROSQiDwqIqnTTkdERJdES3FbAKwE8Bul1AoAgwAeuHAjEblXRCpE\npKKzszPMMYmI6Bwtxd0EoEkptW/s66cxWuQfoZR6RClVrpQqd7s1DdMQEdE0TFrcSqk2AGdEZP7Y\nt64BUBPRVERENCGts0q+CeCJsRkl9QDuiVwkIiK6GE3FrZQ6DKA8wlmIiEgDXjlJRGQwLG4iIoNh\ncRMRGYzWk5NECePJfY2attu8ujDCSYjGxyNuIiKDYXETERkMi5uIyGBY3EREBsPiJiIyGBY3EZHB\nsLiJiAyGxU1EZDAsbiIig2FxExEZDIubiMhgWNxERAbD4iYiMhgWNxGRwbC4iYgMhsVNRGQwLG4i\nIoNhcRMRGQyLm4jIYHjPSYoLvE8kJRIWNyW0YV8Qrx/vQG3rADoGRuDxBuEZ8cPttKMkxwmHnW8R\nij18VVJCquvw4JevnsBLNe0Y8gUBAC6HDSk2C5p6hhBSgNkkKCtIx7pSF/LTk3VOTPQhFjcllJ5B\nH/7v7uN4Yl8jUqxm3LK8ADeV5ePyoizYLKOnfB5/5zQ6+r2oaOjGwcZeHD7TiytL3bhmYQ6sZp4W\nIv1pKm4ROQ1gAEAQQEApVR7JUESR0NI7jJt+sRdt/SO4c3Uh7r+mFNkO+8e2s5hMmJGRjM9kFOBT\nC/Owq7oVb5zoxLG2ftyxqhC5aUk6pCf60FSOuK9WSp2NWBKiCKpq7sPTB87A5bBj69fWYNmsDE3/\nL9lmxq0rZmLxjHRsOdCE375+EptXF6I0xxnhxEQT4999FPeONvfhqfcakZ+ejGe/sVZzaZ9vXq4T\nX/tkMTJTbPjj26dRcbo7AkmJtNFa3ArAbhE5ICL3RjIQUTg1dA3irxVnMDMzGV9aOwc5zukPc2Sk\n2HDf+rkoyXFg66Fm/Me7DWFMSqSd1uJep5RaDuAGAF8XkfUXbiAi94pIhYhUdHZ2hjUk0XR0ebx4\n/J0GZKRYcfcVRR+cfLwUdqsZd62ejQV5Tnz/maN4/J3Tl/w7iaZK0ytZKdU89m8HgG0AVo2zzSNK\nqXKlVLnb7Q5vSqIpCimFvx1oAgB8cc0cpIZxPrbFbMLm1YX41KJc/ODZajz21qmw/W4iLSYtbhFJ\nFRHnuc8BXAfgaKSDEV2Kd052obF7CDeV5SMr1Rb2328xmfCrzSvx6cW5+NH2Gjz6Zn3YH4NoIlqO\nuHMB7BWRSgDvAXheKbUrsrGIpq/L48VLNW2Yn+vE8mmciNTKZjHhl5tX4oYlefifz9fi93t55E3R\nMenfj0qpegDLopCFKCyeq2yBSQSfXVEAEYnoY1nNJjx0xwp888lD+PGOGjjsFnz+8lkRfUwiTgek\nuHLq7CBOdHiwYUEO0pOtUXlMq9mEB+9YjitLXXhg6xE8f6Q1Ko9LiYvFTXFld207nHYLVs/Jjurj\n2i1mPPyFy7CiMBP/9S+H8W59V1QfnxIL1yqhuFHf6cGps4O4cWn+hFP/tC7/Oh0pNgt+9w/luO03\nb+Pexyuw5WtrUJrLKywp/HjETXFBKYXdtR1IS7Jg1Zws3XJkpNjwh3tWwW4144uP7UfHwIhuWSh+\nsbgpLjR2D+F01yDWz3PrvoLfrKwUPPbFy9E16MU3njwEfzCkax6KPyxuigv7TnXDbjGhfLZ+R9vn\nW1KQjp/83VK8d6ob/3vXMb3jUJzhGDfFvMnGpT3eAKqa+7DqvDW1o0HLePkn5mbh/715Ch5vEEsL\n0ifcjrdUo6ngETcZ3oGGHgRDStex7YlsXJqPmZnJeOZQMwZG/HrHoTjB4iZDCymF9051YY4rNSZv\ncGAxmXD7ZTPhD4bw7OEWKKX0jkRxgMVNhnai3YOeIT9Wx+DR9jk5ziRcuzAXNa39qGru0zsOxQEW\nNxlaRUM3Um1mLJqRpneUi1pX6sLMzGQ8V9mCIV9A7zhkcCxuMqwRfxDvtw2gbGYGLKbYfimbRHDr\nigIM+4LY8z7Xq6dLE9uvdqKLqG7pQyCkIroCYDjlpydj5exMvFPfhe5Bn95xyMBY3GRYlWf6kJVq\nw8zMZL2jaHbtwlyYBHippk3vKGRgLG4ypP4RP052erBsZkbEl24Np/RkK9aVuHCkqQ9nuof0jkMG\nxeImQ6pq6oMCsGzWxBe1xKr1pW6k2Mx47f0OvaOQQbG4yZAqm3oxIyPpku7arhe71Ywr5mbjWNsA\n2vu5CBVNHYubDKdn0IemnmGUFRjjpOR4PjE3G1azYO+Js3pHIQNicZPh1LT2AwAWx/jc7YtJtVtw\n2ewsHD7Ti75hXgpPU8PiJsOpbulDXloSsh12vaNcknUlLoSUwtsnedRNU8PiJkPxeANo6BqK+Ssl\ntchKtWFJQTr2n+7GsC+odxwyEBY3GUptaz8UgEX5xi9uAFg9Jwsj/hBeOMobDJN2LG4ylJqWfmSm\nWJGfbrzZJOOZ40pFdqoNT+0/o3cUMhAWNxnGiD+Iuk4PFuWnGeqim4sREZQXZeG9U9042enROw4Z\nBIubDON4+wCCIYXFM4x30c3FrCzMgNkk+CuPukkjFjcZRnVLP1LtFhRmp+gdJaycSVZcsyAHWw42\nwRfgjYVpcixuMgR/MIT32wewKN8JU5wMk5xv06pZOOvxYQ8vgycNNBe3iJhF5JCI7IhkIKLxnOz0\nwBcIYVF+fA2TnHNlqRuZKVbsOMLZJTS5qRxx3w+gNlJBiC6mpqUfdosJxe5UvaNEhNVswvVL8rC7\ntp1zumlSmopbRGYCuBHAo5GNQ/RxIaVQ29qP+XlOWMzxO7p3c9kMDPmCXDWQJqX1XfBzAP8NAM+c\nUNQ1dA1h0BeMm4tuJrJ6bjZcDjt2HGnROwrFuEmLW0RuAtChlDowyXb3ikiFiFR0dvKeehQ+NS19\nsJgE83OdekeJKLNJsHFpHl491gGPlzcUpolpOeJeC+AzInIawFMANojIny7cSCn1iFKqXClV7na7\nwxyTEpVSCtWt/Sh2O2C3mvWOE3E3lc3AiD+EV2rb9Y5CMcwy2QZKqe8C+C4AiMgnAXxbKXVXhHMR\nARidu9075MeG+Tl6R4moJ/c1Ahgdz09LsuDh1+sx6P34ScrNqwujHY1iUPye6aG48GJ1GwTAgjgf\n3z7HJIKF+Wk40TEAf5CnlGh8UypupdQepdRNkQpDdKEXq9tQ5EqFwz7pH4dxY2F+GvxBhXquXUIT\n4BE3xaxTZwdxvN0T97NJLjTXlQqbxYTa1gG9o1CMYnFTzHqxug0A4uKmCVNhMZtQmuPAsbZ+KKX0\njkMxiMVNMevF6jYsKUhDZopN7yhRtzAvDf0jAbT08i7w9HEsbopJbX0jONTYi08vytM7ii7m5Tkh\nAGrb+vWOQjGIxU0x6fmq0cWWbliar3MSfTjGlq+tbWVx08exuCkmba9swaL8NJTkOPSOopuFeWlo\n7RtB37Bf7ygUY1jcFHPOdA/h8Jle3Lxsht5RdDVv7BL/E+2cXUIfxeKmmLN9bJGlm8oSc5jknNw0\nO9KSLDjRwfnc9FEsboo52ytbsaIwA7Oy4usWZVMlIijJcaKuw4MQpwXSeVjcFFPqOgZQ29qPm8sS\ne5jknHm5Dgz7g2jqGdY7CsUQFjfFlO2VrRDhMMk5JW4HBKN3uCc6h8VNMUMphe1HWvCJOdnISUvS\nO05MSLFbMDMzmSco6SNY3BQzalr7Ud85mPCzSS5UmutEU88whny8uQKNYnFTzNhe2QqLSXD9ksS8\nWnIi83IcUADqOLuExrC4KSYopbC9sgXrSl3ISk28tUkupiAzBUlWE4ubPsDipphw6EwvmnuHOZtk\nHGaTYK7LgbpOD1cLJAAsbooRzx1ugc1iwqcW5+odJSaV5DjQO+RHY/eQ3lEoBrC4SXfBkMLzVa24\ner4baUlWvePEpBL36Jote+vO6pyEYgGLm3S371QXOge8nE1yEdkOG9KTrdh7gsVNLG6KAdsrW5Fi\nM2PDgvi+k/ulGL383YG3T3YhGOI4d6JjcZOu/MEQXjjaimsX5iLFljg3BJ6OErcDfcN+VLf06R2F\ndMbiJl3trTuL3iE/h0k0mOtOBcBxbmJxk862V7bAmWTB+nkuvaPEPGeSFQvynBznJhY36WfEH8RL\n1e24fnEe7Baz3nEM4cpSFypO92DYF9Q7CumIxU262fN+JzzeAIdJpmBtiQu+YAgVDd16RyEd8WwQ\n6eZXr9Uh1WZGQ9cQntzXqHccQ1g1JwtWs2Bv3VlcWerWOw7phEfcpItBbwDH2vqxpCAdZpPoHccw\nUmwWrCzMxFs8QZnQJi1uEUkSkfdEpFJEqkXkR9EIRvFtd207/EGFspkZekcxnCtLXahu6Uf3oE/v\nKKQTLUfcXgAblFLLACwHcL2IfCKysSjeba9sRVqSBbOzE/u+ktOxtsQFpYC3T/KoO1FNWtxq1Ln1\nJK1jH7x0i6atb8iP1493oGxmBkzCYZKpWlqQDmeShcMlCUzTGLeImEXkMIAOAC8rpfaNs829IlIh\nIhWdnZ3hzklx5MWatrFhknS9oxiSxWzCFXOzeSFOAtM0q0QpFQSwXEQyAGwTkSVKqaMXbPMIgEcA\noLy8nEfkCe5is0Qee+sUslJtKMhIjmKi+LKu1IWXatrR0DWI2dmpesehKJvSrBKlVC+A1wBcH5k4\nFO883gBOdnpQVpAO4TDJtK0rGb3SlEfdiUnLrBL32JE2RCQZwKcAHIt0MIpPR5v7EFLgbJJLNMeV\nihnpSRznTlBahkryAfxRRMwYLfq/KqV2RDYWxasjTb3IcdqRm2bXO4qhiQjWlrjwcm07giHFufAJ\nRsuskiNKqRVKqTKl1BKl1I+jEYziT9+wH6e7hlA2k8Mk4bCu1IXeIS7zmoh45SRFTVVTLwCgrIDD\nJOGwppjj3ImKxU1Rc6S5DzMykuBycpgkHNxOOxbkOTnOnYBY3BQV3YM+NPUM82g7zNaVuLD/dA9G\n/FzmNZGwuCkqjjaPjsMuLeBFN+G0ttQFXyCEitM9ekehKGJxU1RUNfdhZmYyMlNtekeJK6vHlnl9\ns45XKycSFjdFXM+gD829w1gyg0fb4cZlXhMTi5si7ujYdLUlHCaJiHUlXOY10bC4KeKqxmaTZHGY\nJCLWlo4u8/rOyS69o1CU8NZlFFE9Q6OzST69KFfvKHFhvMW7giGFJKsJf3j7FPqG/di8ulCHZBRN\nPOKmiKpu5jBJpJlNgrkuB+o6PJNvTHGBxU0RdbSlH/npSch28KKbSCrOcaBnyM9x7gTB4qaI6R3y\nobF7iEfbUVDidgAAj7oTBIubIqa6pR8AsJTTACPO5bAhPdmKuo4BvaNQFLC4KWKONvchL41rk0SD\niKDE7cDJzkEEQ7wBVbxjcVNE9A370dA9hCUFaXpHSRgluQ4M+4M4fKZX7ygUYSxuiohqXnQTdfNy\nnDAJ8Eptu95RKMJY3BQRR5v7kOO0I8eZpHeUhJFsM2N2dipeqe3QOwpFGIubwq6jfwQNXUNcCVAH\nC/OceL99AGe6h/SOQhHE4qaw21XdBgUOk+hhQf7oOYVXj/GoO56xuCnsdla1wu20IzeNwyTR5nLY\nMdedit0c545rLG4Kq84BL9471c0lXHV0zYIc7Kvvhscb0DsKRQiLm8Lqxeo2hBTvdKOnaxbmwhcM\n4c3jvLlCvGJxU1jtrGrFXFcqctN40Y1eymdnIjPFiher2/SOQhHC4qawOevx4t36Lmxcmg8R0TtO\nwrKYTbhuUR5213bwJsJxisVNYXNumOSGpXl6R0l4G8vy4fEGsPcEb2kWj1jcFDY7KkeHSRbl8zJ3\nva0pzkZ6shU7q1r1jkIRwOKmsOgYGMG+U124sYzDJLHAajbhukW5eLm2Hd4Ah0vizaTFLSKzROQ1\nEakRkWoRuT8awchYdh0dHSa5qWyG3lFozMal+RgYCfAO8HFIyxF3AMA/KqUWAfgEgK+LyKLIxiKj\n2VHZitIcB+bnOfWOQmPWlrjgTLJgZxVnl8SbSYtbKdWqlDo49vkAgFoABZEORsbR1jeC/Q3duLEs\nX+8odB6bZXR2yYtH2zi7JM5MaYxbRIoArACwb5yf3SsiFSJS0dnJif+JZGdVKxSHSWLSrSsKMOAN\ncMXAOKO5uEXEAWALgG8ppfov/LlS6hGlVLlSqtztdoczI8W4HUdasCDPiZIch95R6AJXFGcjN82O\nbYea9I5CYaSpuEXEitHSfkIptTWykchImnuHcbCxFzdxmCQmmU2Czy4vwJ73O9Hl8eodh8JEy6wS\nAfA7ALVKqZ9FPhIZyc4jo/OEOUwSu25dWYBASGF7ZYveUShMtBxxrwXwBQAbROTw2MfGCOcig9hR\n1YolBWkocqXqHYUmsCAvDQvz07DtULPeUShMtMwq2auUEqVUmVJq+djHzmiEo9h2pnsIlWd6ceNS\nHm3HuttWFqCyqQ91HQN6R6Ew4JWTNG07Phgm4fh2rLtleQEsJsFf9p/ROwqFAYubpu35qhYsm5mO\nWVkpekehSbiddly3OBdPH2jinO44wOKmaanv9OBocz9PShrIHasK0TPk5zrdccCidwAypm2HmmES\n4DPLWdyx5sl9jeN+P6QUMlOs+PnuExj0BrF5dWGUk1G48IibpiwUUth6sBnrSt28IbCBmERweVEW\nTp0dxNkBzuk2MhY3Tdl7p7vR3DuM21ZyyRqjuWx2Jkwy+hyScbG4acq2HmyCw27BdYt4pxujcSZZ\nsXhGOioaujHIu8AbFoubpmTYF8TOqjZsXJqHZJtZ7zg0DWuLszHiD2HLQa5fYlQ8OUlT8s/PVcPj\nDSAt2TrhSTCKbbOyUjAzMxmPvXUad62eDZOJdywyGh5x05RUNHQjM8WKomxe4m5UIoK1JS6cOjuI\nPce53KsRsbhJs1NnB1HfOYjLi7Jg4n0lDW3JjHTkpSXhd3tP6R2FpoHFTZo9tb8RJgFWzs7UOwpd\nIrNJcPea2XirrgvH2j62vD7FOBY3aeILhPB0RRMW5KUhLcmqdxwKg82rCpFkNeGxvaf1jkJTxOIm\nTV6uaUfXoA+r5mTpHYXCJCPFhttWzsS2w828yYLBcFYJafLkew0oyEjm7cniyJP7GuF22OELhPCd\nLVXYsCDnY9vwsvjYxCNumtSJ9gG8VdeFTZfP4knJOJOTloR5uQ7sq+9CIBjSOw5pxOKmST329mnY\nLCYefcWpNcUuDHgDONLUp3cU0ojFTRfVO+TD1oNNuHV5AbIddr3jUASU5jiQl5aE1090IqSU3nFI\nAxY3XdSf3zuDEX8I96wr0jsKRYiI4Kp5bnQOeHGslVMDjYDFTRPyB0N4/J3TWFOcjQV5aXrHoQha\nUpCOrFQb9hzvhOJRd8xjcdOEdla1orVvBPesnaN3FIows0mwvtSNpp5hnOwc1DsOTYLFTeMKhRR+\n/dpJlOY4cM0408Qo/qwszIAzycL1SwyAxU3j2l3bjvfbB/Cfry7m6nEJwmI2YV2JC/WdgzjTPaR3\nHLoIFjd9jFIKv3qtDoVZKbiZNwNOKKuKspBsNWPP8U69o9BFsLjpY/bWnUVlUx++elUxLGa+RBKJ\n3WrGFcXZqG3tR1v/iN5xaAJ8V9JHKKXw4O4TyE2z47bLeE/JRLRmbjZsZhPe4FF3zJq0uEXk9yLS\nISJHoxGI9PVKbQcqGnrwzQ2lsFt4a7JElGK34PKiTBxp6kVjF8e6Y5GWRab+AOCXAB6PbBTS25/e\nbcBDr5xAdqoNSoG3JktgV5a6se9UN37x6gn8n88t0zsOXWDSI26l1BsAuqOQhXR2uLEXHQNeXLc4\nD2bOJEloaclWrJ6Tha2HmnHqLOd1xxqOcRMAYMQfxO7adhRkJGPJDF4lScD6eW5YzYKHXjmhdxS6\nQNiKW0TuFZEKEano7ORJDaN5+PV69A77cf2SPAiXbiUAziQr7r6iCM8ebkZdx4Deceg8YStupdQj\nSqlypVS52+0O16+lKGjsGsKv99RhaUE6it28UQJ96L71c5FkNeNnLx/XOwqdh0MlhB/vqIbZJNi4\nNF/vKBRjsh12fOXKudhZ1YaDjT16x6ExWqYD/hnAOwDmi0iTiHw58rEoWl6uacfu2g5869pSpCfz\nJsD0cV9ZPxcuhx0/2VnLlQNjhJZZJXcopfKVUlal1Eyl1O+iEYwir3fIh+9tq8KCPCdXAKQJOewW\nfOvaUuw/3YOXatr1jkPgUElC+/6z1ege9OHfP78MVl7aThex6fJZmOtOxU9fOAZfgPem1BvfrQlq\nx5EWbK9swf3XlGLxjHS941CMs5hN+O83LkT92UH8/q1TesdJeCzuBNTYNYTvbTuKZbMy8LVPFusd\nhwxiw4JcXLswFw/uPoGW3mG94yQ0FneCGfYF8dU/HYBSCg9tWs7V/2hK/vnmRVBQ+PH2Gr2jJDS+\naxOIUgrfe6YKNa39+Pmm5Zidnap3JDKYWVkp+OaGUuyqbsOrx3iiUi8s7gTyyBv12HqwGfdfU4oN\nC3L1jkMG9Z+unIN5uQ48sKUKvUM+veMkJBZ3gthyoAk/eeEYbizLx/3XlOodhwzMbjHjZ59fju5B\nH77/bLXecRKSlmVdyeB+8OxR/OndBhS7U7G6KAtP7T+jdyQyiIst7Xv1ghxsr2xBis2Mn95WFsVU\nxCPuOPdSdRue2NeIvPQk3Ll6Nk9GUtisL3VjVmYynjnUjPpOj95xEgrfxXHsucoWfO2Jg5iRnoQv\nrx1dLIgoXMwmwaZVhTCbBPf9xwF4vAG9IyUMFnccUkrhkTdO4v6nDuGy2Zn40to5SLaxtCn8MlNs\nuGNVIU52evBPf6vkWiZRwuKOM75ACA9sqcK/7DyGjUvy8cd7VsHOI22KoGK3Aw/csAAvHG3Dv+46\npnechMCTkwZ3/smjvmE//vxeIxq7h3D1fDeuKM7GtkPNOqajRPGVK+fiTPcwHn69HtmpNty7nlfk\nRhKLO07UdXjwl/2N8IcUNl0+C2UzM/SORAlERPDDzyxGz5AP/7LzGJxJVtyxqlDvWHGLxW1wgVAI\nL9e0Y++Js3A77di8uhA5ziS9Y1ECMpsEP/v8cni8AXx3axVG/EEuFxwhLG4Dq+vw4Ld7TqKlbwSr\nirKwcWk+bBaetiD92CwmPPyFy/Bf/nwIP9peg0FvAF+/uoT3MQ0zvssNSCmFJ/Y14KZfvIneYT/u\nWj0bn11RwNKmmGC3mPGrzSvxdysK8G8vHccDW6q4hneY8YjbYLo8XnxnSxV217bjylIX1pa4kJbE\nW45RbLGYTfi3zy3DzMxkPPRqHU53DeI3d12GrFSb3tHigkRi3mV5ebmqqKgI++9NdG8c78Q//q0S\nfUN+fOeGBbhnTREvX6eYd/hMD7YebEaO045f3rkSKwsz9Y4Uk0TkgFKqXMu2/NvaAAZG/Pju1irc\n/fv3kJFsxTNfX4svr5sDk4njhhT7ls/KxH3ri2E2C/7+4Xfw6Jv1CIV4oc6l4BF3jDo3P/t4+wC2\nHWpG/7Af60pcuHZRLu8PSYY07Avi6YNNqG3tR0mOA7etnIn05I8P821enZjTCHnEHQeGfUFsOdCE\nP7x9GjaLCV+9qhg3LM1naZNhJdvMuGt1IW5ZPgMNXYN46JUTONDQw8vkp4EnJ2NMKKTw9MEm/Hz3\ncQz6ArhqnhsbFuSwsCkuiAhWz8lGscuBLQebsOVgEw6f6cEtywvgctj1jmcYLO4YoZTC68c78e8v\nHUdVcx8Ks1Jw9xVFKMhM1jsaUdi5nHZ8Zf1c7D/djV1H2/Dg7hO4ojgbV8/P0TuaIbC4deYLjF75\n+PAbJ3GkqQ8FGcl4cNNyeEYCvGiB4ppp7Oh7UX4aXq5px1t1Z3GgoQeDvgD+YU3RuOPfNIonJ3UQ\nDCkcauzBi9Vt2HqwGV2DPhRmpeDrVxfj1hUzYbOYLnrnEaJ41NI7jN217TjWNgCn3YLby2fiztWF\nKMlx6h0tKqZycpJH3FEwMOLH8XYPDp/pxYGGbuyr70bXoA8Wk+CahTnYtKoQ60vdMHN6HyWwGRnJ\nuPuKIiyblY7fvl6PP73bgMfeOo3lszLw6cV5uG5xLua6UvmXKDQWt4hcD+BBAGYAjyql/jWiqQwm\nEAyh0+PF4283oHvIh55BH7oHfege8qHb40PvsP+DbTNSrLh6fg42LMjBVfPdvOqR6AKLZ6TjF3es\nwFnPImw50IQdR1rx013H8NNdx+B22nF5USb8QQWXww63w45sh23Ck/fxOrVw0uIWETOAXwH4FIAm\nAPtF5DmlVE2kw0WTUgreQAhDviAGvQEM+gKj/3rPfT36r8cbQHv/CNr6RtDeP4LWvhGc9Xhx4fUE\nqXYLslKsmJ2dglVpSchNS8KMjGSkJ1vj9sVEFE4uhx33XVWM+64qRnPvMPa834GK0z3Yf7obTT3D\nH2wnANKTrUhLtsKZZIHDboEzyQKn3Yocpx3usQ+Xwx436/lMOsYtIlcA+KFS6tNjX38XAJRSP5no\n/4RrjFsphZACQkohGFJQCgiOfe71BzHiD8EbGP13JBDEyNj3hnznF27gg+Id8gbg8QbHfj76vc4B\nL7yBIHyB0MfKdyJpSRbkpY+WcV5a0gefv982gMxUGzJTrLBbeNcZounQcmDz2Fun0OXx4azHi06P\nF10eHwbzXjicAAAFOUlEQVRG/BgYGT24GvIFx/1/6clWZKfaxt6nNmSlWpGZakN2qg0ZyTak2M1I\nsZmRbLUg9dznNgtSrGaYzQKTCMwiMJkAswjMJgnb0E24x7gLAJy/IEYTgNXTCTaZ5T9+CUPeIIJK\nIaRGizockq3msSfBglS7Bak2MzJSbCjINCM9yQqbxQS7xXTev+YPPj/3/c9dNgupdjNS7ZYJb7rL\nE4pE0WG3mDEjIxkzMsafLhsIhTDoDWJtSTY6B7wffni86Bnyo3vQi+beYRxt7kP3oA++4PRXLxTB\nB4Xudtrx1gMbpv27tArbyUkRuRfAvWNfekTk/Yts7gJwNlyPHQ3fHv/bhtuPcXAfYgP3YcydYQhy\nCS5pH04AGB2TmJbZWjfUUtzNAGad9/XMse99hFLqEQCPaHlQEanQ+idBLIuH/eA+xAbuQ2wwyj5o\nGanfD6BUROaIiA3AJgDPRTYWERFNZNIjbqVUQES+AeBFjE4H/L1SqjriyYiIaFyaxriVUjsB7Azj\n42oaUjGAeNgP7kNs4D7EBkPsQ0QueSciosiJj9noREQJJKLFLSLXi8j7IlInIg+M8/MFIvKOiHhF\nZIIZd/rSsA93isgREakSkbdFZJkeOS9Gwz7cMrYPh0WkQkTW6ZFzMpPtx3nbXS4iARG5PZr5tNDw\nXHxSRPrGnovDIvIDPXJejJbnYWw/DotItYi8Hu2Mk9HwPPzTec/BUREJikiWHlnHpZSKyAdGT2Se\nBDAXgA1AJYBFF2yTA+ByAP8LwLcjlSXC+7AGQObY5zcA2Kd37mnsgwMfDpuVATimd+7p7Md5272K\n0XMyt+udexrPxScB7NA76yXuQwaAGgCFY1/n6J17Oq+l87a/GcCreuc+/yOSR9yrANQppeqVUj4A\nTwG45fwNlFIdSqn9APzj/YIYoGUf3lZK9Yx9+S5G57nHEi374FFjr1AAqQBi8cTHpPsx5psAtgDo\niGY4jbTuQyzTsg+bAWxVSjUCo+/zKGeczFSfhzsA/DkqyTSKZHGPd6l8QQQfLxKmug9fBvBCRBNN\nnaZ9EJFbReQYgOcBfClK2aZi0v0QkQIAtwL4TRRzTYXW19OasaGrF0RkcXSiaaZlH+YByBSRPSJy\nQETujlo6bTS/r0UkBcD1GD0YiBlcjztMRORqjBZ3TI4PT0YptQ3ANhFZD+B/ALhW50jT8XMA31FK\nhQy8ZvNBjA4xeERkI4BnAJTqnGmqLAAuA3ANgGQA74jIu0qp4/rGmpabAbyllOrWO8j5Ilncmi6V\nj3Ga9kFEygA8CuAGpVRXlLJpNaXnQSn1hojMFRGXUiqW1s7Qsh/lAJ4aK20XgI0iElBKPROdiJOa\ndB+UUv3nfb5TRH4dY8+FluehCUCXUmoQwKCIvAFgGYBYKe6pvCc2IcaGSQBE9OSkBUA9gDn48ATA\n4gm2/SFi8+TkpPsAoBBAHYA1eue9hH0owYcnJ1di9EUsemef7utpbPs/IPZOTmp5LvLOey5WAWiM\npedC4z4sBPDK2LYpAI4CWKJ39qm+lgCkA+gGkKp35gs/InbErSa4VF5Evjr289+KSB6ACgBpAEIi\n8i2Mnt3tn/AXR5GWfQDwAwDZAH49dqQXUDG0SI3GfbgNwN0i4gcwDODv1dgrN1Zo3I+YpnEfbgfw\nNREJYPS52BRLz4WWfVBK1YrILgBHAIQweteso/ql/qgpvJZuBfCSGv3LIabwykkiIoPhlZNERAbD\n4iYiMhgWNxGRwbC4iYgMhsVNRGQwLG4iIoNhcRMRGQyLm4jIYP4/hSq6WgMTmIkAAAAASUVORK5C\nYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x1228ef8d0>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import regreg.api as rr\n",
"n, p, s = 1000, 20, 3\n",
"beta = np.zeros(p)\n",
"beta[:3] = 2 / np.sqrt(s)\n",
"\n",
"beta_star = []\n",
"nsim = 1000\n",
"for i in range(nsim):\n",
" X = np.random.standard_normal((n, p))\n",
" linpred = X.dot(beta)\n",
" pi = normal_dbn.cdf(np.log(1 + np.exp(linpred)))\n",
" Y = np.random.binomial(1, pi)\n",
" loss = rr.glm.logistic(X, Y)\n",
" beta_hat = loss.solve()\n",
" beta_star.append(beta_hat)\n",
"beta_bar = np.mean(beta_star, 0)\n",
"beta_star = np.asarray(beta_star)\n",
"sns.distplot(beta_star[:,0])\n"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.axes._subplots.AxesSubplot at 0x1224e6b00>"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd0m/d5L/Dvg0USACcAbnGJFLWHTQ1bViwpHrE8Gid2\n4zhxTnJz4yRNmyatb9s0TZq29570nnOdum3ipmqSOtPOsJ3hWI7tWLYkyxqUrEEtkqJIiksEJwgO\nEON3/wAlS44kvgQBvBjfzzk8kiiQ/L4C+OjH5/0NUUqBiIiSh0HvAERENDcs3ERESYaFm4goybBw\nExElGRZuIqIkw8JNRJRkWLiJiJIMCzcRUZJh4SYiSjKmWHxSp9OpqqqqYvGpiYhS0qFDhwaUUi4t\nj41J4a6qqkJjY2MsPjURUUoSkQ6tj2WrhIgoybBwExElGRZuIqIkw8JNRJRkWLiJiJIMCzcRUZJh\n4SYiSjIs3ERESUbTAhwRaQcwBiAIIKCUaohlKCIiura5rJzcopQaiFkSoij7yf5OTY97eH1FjJMQ\nRRdbJURESUZr4VYAXhWRQyLy6NUeICKPikijiDS63e7oJSQioitoLdy3KKVWA7gLwOdE5D3vfoBS\nartSqkEp1eByadrgiijufP4gJnwBKKX0jkIUMU09bqVU98yv/SLyPIB1AHbFMhjRfIWUwrmBcTT3\njaGl34uh8WlMB0MAALNRkG+1oL44G6sW5GJpSQ5EROfERNrMWrhFxAbAoJQam/n9HQD+MebJiCI0\nPD6NZw6ex3/uOouRCT+MIqh0WrHQlY/sTDOMBsHIxDT6x3x4s3UAd//bHtxQkYev3LMUayry9Y5P\nNCstI+4iAM/PjEZMAH6ilHoppqmIIuCZ8uM7u9rw3T3nMD4dRI3ThruWl2BRkR0ZJuNVP2bcF4DF\nZMA3d7bi/if34gNryvD39y1DbpY5zumJtJu1cCul2gCsikMWoogEQwo/OdCJx18+g5EJP7atKMbn\n31uHwx0js36sLcOEh9dX4IM3luPJna3YvqsNhzqH8e2P3oglJTlxSE80dxKLmzQNDQ2KJ+BQrFw+\nP7treAK/PNKNnpEp1Dht2LaiBKV5WRF/7o7BcTx9oBMT00E82LAAK8pyL/0d53tTLInIIa2LG2Ny\ndBlRrAVDCjvP9OP1M/2wZ5jw0NpwkZ3vDcZKhw2f21KLn+zvxDMHOjF9QzlurGTfmxILCzclHc+k\nHz/e34Hzw5NYvSAP964sRZbl6j3sSGRnmvGJjdX40f4OPHu4C/5gCBtqHJpWYnJUTvHAwk1J5XjX\nKJ58vRVT/hAeWrsAK8vzYvJ1LCYDHtlQiacPdOLXR3tgyzBd0TYh0hOXvFPS2N3ixoP/uRcGEXz6\n1pqYFe2LzEYDPryuApUFVvy88Tw6Bsdj+vWItGLhpqRwqGMIj/7gEKocNnx280KU5EZ+A3IuzEYD\nPrqhErlZZvxwXweGxqfj8nWJroeFmxLeyR4PPvHfB1Gcm4kffnI9sjPjO8falmHCx2+uQkgp/PRg\nJ4IhLpcnfbFwU0IbnfTjUz9ohC3DhB9+ch1c2Rm65HDYM/D+1WU4PzyJ35++oEsGoot4c5LiIpIZ\nGUopfPn547jgmcIvPnszyvOtsYqnycryPLRc8OKNM27Uuuyocdl1zUPpiyNuSljPHu7GC8d68cXb\nF2H1gtjeiNTqnlUlKLBZ8Nzb3QjMbFhFFG8s3JSQekcn8fe/asK66gJ85taFese5JMNkxH2rSzE0\nPo09rTwQivTBwk0J6fGXm+EPKjz+4CoYDYm13WpdYTaWluRg55l+jE769Y5DaYg9bkoYF/vgvaOT\nePZQFzbWOrG7JTFHtdtWlOCJV8fwUlMvPrSWqyUpvjjipoTzuxN9yDQbsaW+UO8o11Rgs2BTnQtH\nu0bRPTKpdxxKMyzclFBa+71ovuDF5npXVPcfiYVNdU5kmAx4o5lnrFJ8sXBTQnm9uR+5WWZsqHHo\nHWVWmWYjbqpx4ET3KAbGfHrHoTTCwk0Jo88zhTb3ODbUOGA2JsdL8+ZaJ4wGwRstHHVT/CTHdwel\nhbfODsJkEKxNov2v7RkmrK0qwNudwxiZ4D4mFB+cVUIJYWI6gCPnh7F6QR6sGcn1stxU58T+c4N4\nq20QeVaLpo/hvt00HxxxU0JobB+GP6hw80Kn3lHmLM9qweLiHBzuHOEGVBQXLNyku5BS2HduEDVO\nG4pzM/WOE5GGqnyM+wI43efROwqlARZu0l374DhGJvxoqCrQO0rE6gqzkZNpQmP7sN5RKA2wcJPu\njp4fhcVowNKSHL2jRMxoENxQmY/mC2NcBk8xx8JNugqEQmjqHsXS0hxYTMn9cmyoLIACcKiDo26K\nreT+TqGk13LBi0l/EKvKk/8g3gKbBTUuGw53DkMp3qSk2GHhJl0dOT8Cq8WI2sJsvaNExeryPAyN\nT6NndErvKJTCWLhJNz5/EKf7PFhRlptwW7dGaklJDgwCnOge1TsKpTAWbtLN6b4x+IMKq8oT43Sb\naLBlmFDltKGpx8N2CcUMCzfp5lSfB7YMEyoc+p4lGW3LS3Mx4PWhnxtPUYywcJMugiGF5gtjWFyc\nDYOkRpvkoqWlORAATT1sl1BssHCTLtoHxzHlD2FJcWrclLxcTqYZFQVWnOjmKkqKDRZu0sXpXg9M\nBkmZ2STvtqwsF32eKQx62S6h6NNcuEXEKCJvi8gLsQxEqU8phVN9Y6hx2ZJ+0c21XFwFerpvTOck\nlIrm8l3z5wBOxSoIpQ/3mA9D49NYksRL3GdTYLPAac9ASz8LN0WfpsItIuUA7gbwndjGoXRwcRS6\nuDh1CzcALCqyo809Dn8wpHcUSjFaR9xPAPgrANd8BYrIoyLSKCKNbjePcaJrO93nQWluJnKzzHpH\nialFRdkIhBTODYzrHYVSzKyFW0TuAdCvlDp0vccppbYrpRqUUg0ulytqASm1+PxBdA5NoK4oNW9K\nXq7aaYPJIGi5wHYJRZeWEfdGAPeJSDuAZwBsFZEfxTQVpaxzA+MIKaC20K53lJgzGw2odtrQfMGr\ndxRKMbMWbqXUl5RS5UqpKgAPAXhNKfXRmCejlNTS74XZKKgsSK3VkteyqCgbbq8Pw+M8SJiiJzXn\nYlHCau33osphg8mYHi+9uqLwTxbNnF1CUTSn7x6l1OtKqXtiFYZS2+ikH26vLy3aJBe57BnIs5rR\nwnYJRVF6DHsoIbTOjDrTqXCLCBY67TO9fe4WSNHBwk1x09LvhT3DhOKc5DzJPVI1Lhsm/UH08XAF\nihIWboqLkFI42+9FbaEdkmK7Ac6mxhX+CaON87kpSli4KS4ueKYwPh1ErSt92iQX5WaZ4bBZ0OZm\nn5uig4Wb4uLi6sFql03nJPqocdnQPsg+N0UHCzfFRZt7HPlWM/KtFr2j6KLaaceUP4TeEfa5af5Y\nuCnmQjP7ddQ4069NclHNzE8abQNsl9D8sXBTzDX3j2HSH0S1Mz3bJED4VBynPQNtbt6gpPlj4aaY\n23d2EADSunAD7/S5gyH2uWl+WLgp5va1DYX727b07G9fVO20wRcIoXd0Uu8olORYuCmmQiGF/ecG\nUZ3G/e2Lqhzhnzg6Bid0TkLJjoWbYqq5fwzDE37UpHmbBAjP5863mtE+yD43zQ8LN8XU/rYhAOxv\nX1TlsKF9YByK87lpHli4Kab2tQ2iLC8r7fvbF1U5bBifDvI4M5oXFm6KmXB/ewgbahx6R0kYlc7w\nARIH24d0TkLJjIWbYqal34uh8WlsqCnQO0rCcNkzYLUYceDcsN5RKImxcFPM7D8Xnr/NEfc7RARV\nDhsaOzjipsixcFPMXOxvL0iT8yW1qnJY0TE4gX4P9y2hyLBwU0wopbCvbQjr2Sb5A5Uz87kPsM9N\nEWLhpph4p7/NNsm7leZlIctsRGM7+9wUGRZuiol9beH+9k0s3H/AaBDcUJmHA+c44qbIsHBTTOxv\nG0JZXhbK87P0jpKQGioLcLrPA8+UX+8olIRYuCnqwv3tQayvLki78yW1WlddgJACDnewXUJzx8JN\nUdfS78Xg+DQ2LGSb5FrWVOTBaBAuxKGIsHBT1L11lv3t2VgtJiwvzcFB3qCkCLBwU9S9dZbzt7VY\nW1WAI+dH4AsE9Y5CSYaFm6Lq4v7bN7FNMqu11QWYDoRwvGtU7yiUZFi4KarOXAjvv83527NrqMwH\nALZLaM5YuCmqLvW3OeKelcOegYUuG29Q0pyxcFNU7WsbREWBFWV5nL+txdqqAjS2DyHEA4RpDli4\nKWre2X+b+5NotbaqAJ6pAJr7x/SOQkmEhZui5mSvB6OTfrZJ5mBddfg/uYNc/k5zMGvhFpFMETkg\nIkdF5ISI/EM8glHyubg/CW9Maleen4WinAzeoKQ5MWl4jA/AVqWUV0TMAPaIyA6l1L4YZ6Mks69t\nEFUOK0py2d/WSkSwtqoAB9uHoJTiFgGkyawjbhXmnfmjeeaNd1LoCsGZ/jbbJHO3rroAvaNT6Bqe\n1DsKJQlNPW4RMYrIEQD9AF5RSu2/ymMeFZFGEWl0u93RzkkJ7kTPKMamAmyTRKChMtzn5nFmpJWm\nwq2UCiqlVgMoB7BORJZf5THblVINSqkGl8sV7ZyU4Lj/duTqi7ORnWniAcKk2ZxmlSilRgDsBPC+\n2MShZPXW2UHUuGwozMnUO0rSMRoEDZX5XIhDmmmZVeISkbyZ32cBuB3A6VgHo+QRCIZwsH2Yo+15\naKgqQOvMcW9Es9Ey4i4BsFNEjgE4iHCP+4XYxqJk0tTjgdfH/vZ8XJzP3chRN2kw63RApdQxAGvi\nkIWS1MX9SVi4I7eyPBcWkwEH24dwx7JiveNQgtMyj5vouvaeHUBtoR2u7Ay9oySNn+zv/IP3leRm\n4qWmPlQ77Zfe9/D6injGoiTBJe80L1P+IA6cG8IttU69oyS9KocN3SOTmA6E9I5CCY6Fm+blUMcw\nfIEQNtWxcM9XlcOGkALOD0/oHYUSHAs3zcvulgGYjcL+dhRUFFghANoHx/WOQgmOhZvmZXeLG2sq\n8mHL4O2S+cqyGFGcm4mOAY646fpYuClig14fTvR48B62SaKm0mFF59AEgjxYga6DhZsi9ubMNMBb\n6rjFQbRUOWyYDobQO8oNp+jaWLgpYrub3cjNMmNFWa7eUVJGpcMGAGgfZLuEro2FmyKilMKe1gFs\nrHXAaOAe0tGSm2VGgc2C9gHeoKRrY+GmiJx1e9E7OoVbatkmibYapw1tA16EFPvcdHWcCkAR2Xk6\nvOf6yMT0VVcBUuRqC+1o7BhGNw9WoGvgiJsi8npzP+qLspFntegdJeXUuMJL3lvd3lkeSemKhZvm\nzOsL4MC5IWxezDZJLNgzTCjNzURrPws3XR0LN83Zm60D8AcVNi8q1DtKylpYaEfn0AQmpgN6R6EE\nxMJNc/b6GTfsGSY0VOXrHSVl1brsCIYUDrbzODP6QyzcNCdKKbx+ph+31DphNvLlEyuVDhuMBsGb\nrQN6R6EExFkldFXXminS55lC7+gUbqoxcjZJDFlMBlQWWLG7hYWb/hCHTDQnZ/rGAAB1Rdk6J0l9\ntYV2nOr1oH9sSu8olGBYuGlOTvd6UJqbidwss95RUt6imf8cXz/j1jkJJRoWbtLM6wugc2gCi0ty\n9I6SFkpyM1Gck4mdp/v1jkIJhoWbNGvuG4MCsISFOy5EBFsWF2J3ywCPM6MrsHCTZqf6PMjJDC8O\nofjYurgQXl8AB9uH9I5CCYSFmzTxB0NoueDF4pIciHA3wHjZWOuAxWTAa2yX0GVYuEmTNvc4poMh\nLClmmySerBYTbqpxsHDTFVi4SZPTfR5YjAbUuGx6R0k7WxcX4tzAONq46RTNYOGmWYWUwqleD2oL\n7VwtqYOti8N7wvz+FEfdFMbvQppV9/AkPFMBLCtlm0QPCwqsWFqSgx1NvXpHoQTBwk2zOtHjgUGA\nxexv6+bulSU43DmCnhEerkAs3DQLpRRO9IyixmVHlsWod5y0ddfyYgDAjqY+nZNQImDhpuvqH/Nh\ncHyabRKd1bjsWFycjR3H2S4hFm6axYmeUQi4WjIR3L2iBI0dw+gb5aZT6W7Wwi0iC0Rkp4icFJET\nIvLn8QhGieFkjwcLCqzIyeSmUnrbtrIEAHiTkjSNuAMA/lIptRTABgCfE5GlsY1FiWBofBo9o1Ns\nkySIhTPtkt8eY+FOd7MWbqVUr1Lq8MzvxwCcAlAW62Ckv6buUQDAstJcnZPQRfeuKkVjxzA6Byf0\njkI6mlOPW0SqAKwBsD8WYSixHOseQXl+FgpsFr2j0Iz715RBBHj2cJfeUUhHmgu3iNgBPAvgC0op\nz1X+/lERaRSRRrebG78nuwGvDz0jU1hZnqd3FLpMaV4WNi504rm3uxAKKb3jkE40FW4RMSNctH+s\nlHruao9RSm1XSjUopRpcLlc0M5IOjnWF2yQrytgmSTQfvLEM54cmudVrGtMyq0QAfBfAKaXUN2If\niRLB8e4RVDmsPKIsAd25rBg2ixG/OMR2SbrSMuLeCOARAFtF5MjM27YY5yIdnekbwwWPDyvYJklI\nVosJd68swYvHezExHdA7DulAy6ySPUopUUqtVEqtnnl7MR7hSB8vHOuBAFjOaYAJ64EbF2B8OogX\nODUwLXHlJF1BKYUXjvWixmVDNhfdJKy1VfmoK7Tjx/s69I5COmDhpiuc6PHg3MA4Z5MkOBHBIzdV\n4mjXKI6eH9E7DsUZCzdd4TdHe2AyCFdLJoH715TBajHiRxx1px0WbrrkYptkU50TVotJ7zg0i+xM\nM96/pgy/PtqDkYlpveNQHLFw0yWHO0fQPTKJe1aW6h2FNPro+kr4AiFODUwzLNx0yQvHemAxGXD7\nsiK9o5BGS0tzsLYqH0/tbUcgGNI7DsUJCzcBAIIhhd8e68WWehe3cE0yn9pUg67hSbzI03HSBhuZ\nBAB46+wg+sd8uHcV2yTJ5rYlRXDaM/DPO05hbNKP8GLnq3t4fUUck1GscMRNAMK7zWVnmnDbErZJ\nko3BINhU60TPyBTaBsb1jkNxwMJN8PoCeKmpD/esLEWmmQcCJ6PVFXmwZ5iwq5k7c6YDFm7Ci8d7\nMekP4oEbeT5GsjIbDbh5oQMt/V70jk7qHYdijIWb8NzhLlQ7bbihIl/vKDQP66sdsBgN2NMyoHcU\nijEW7jR3fmgC+9qG8IE1Zde9qUWJL8tiRENVPo52jXBBTopj4U5zzx3uBgC8fw3bJKlgY60TALD3\n7KDOSSiWWLjTWDCk8LPG87il1okFBVa941AU5FstWFGWiwPtQ5icDuodh2KE87jT2K4WN7pHJvG3\n25boHYWu4Sf7O+f8MZvqXDjaNYr95waxub4wBqlIbxxxp7FnDnTCYbPg9qWcu51KSvOysKjIjj2t\nA/AFOOpORSzcaarfM4VXT/XjgRvLYTHxZZBqttYXYmI6iAPneKBwKuJ3bJr6+aEuBEMKH1q7QO8o\nFAMVDhtqXXbsbhnAdICbT6UaFu40FAwpPHOwExtqClDjsusdh2Jky+JCeH0BHGznqDvVsHCnoZ2n\n+3F+aBKPbKjSOwrFULXThmqnDbtb3PBzy9eUwsKdhp7a247inEzcwX23U96W+kJ4pgI41DGsdxSK\nIhbuNNNyYQx7WgfwyE2VMBv59Ke6hS4bKgqseKPZjUCIo+5UwXncaebvftkEk0FgNhoimiNMyUVE\nsHVxIZ7a2463O0bwsZv0TkTRwCFXGhmd9ONw5zBWloe3AKX0UFdoR3l+Fl5v7mevO0WwcKeRpw90\nwh9UuHmhQ+8oFEcigi31hRie8ONXR3r0jkNRwMKdJqb8QXx3zznUuuwozcvSOw7F2eLibJTkZuJb\nO1sRDCm949A8sXCniecOd8M95sN7Frn0jkI6uDjqPjcwjheOcdSd7Fi400AwpLB911msLM/FQpdN\n7zikk6WlOagvysa/v9aKEEfdSY2FOw281NSH9sEJfObWhTwsIY0ZRPCnW2vR2u/Fi029eseheWDh\nTnGhkMI3d7ai2mnDncuK9Y5DOtu2ogS1hXb866st7HUnMRbuFPe7E3041evB599bC6OBo+10ZzQI\nvnjbIrT0e/GrI916x6EIzVq4ReR7ItIvIk3xCETREwop/Murzahx2XDfKh5NRmF3LS/GstIcPPFq\nC3cOTFJaRtxPAXhfjHNQDLxwvBfNF7z4wm2LONqmSwwGwWN31KNzaAI/azyvdxyKwKzL55RSu0Sk\nKvZRaL4uX8IeDCn82+9bUJidAc+kn8vb6Qqb611oqMzHv7/WggduLEem2ah3JJoD9rhT1OGOYbi9\nPty2pAgGziShdxERPHZnPS54fPjhWx16x6E5ilrhFpFHRaRRRBrdbne0Pi1FwBcI4pVTF1BZYMWy\n0hy941CC2lDjwKY6J558vRVjU36949AcRK1wK6W2K6UalFINLhdX5+lpV/MAvL4Atq0o4bxtuq7H\n7qjH8IQf39vTrncUmgO2SlLM6KQfe1rdWFmeiwUFVr3jUIJbtSAPdy4rwnd2t2F4fFrvOKSRlumA\nTwN4C0C9iHSJyCdjH4si9bsTfQgp4I6lXGxD2vzlHfXwTgfwzZ2tekchjbTMKvlwPILQ/LW5vThy\nfgRb6gtRYLPoHYcS0LVmF91YkY+n3mxHbpYZn39vXZxT0VyxVZIipgMh/OpoD/KtZmyu5z0Gmpvb\nlxbBaBTsaOrTOwppwMKdIr735jm4x3y4d2Upz5KkOcvONGPzIhdO9Xqw9+yA3nFoFvwOTwEdg+N4\n4tVmLCnOxuISTv+jyGysdSLPasY//uYkAjziLKGxcCe5UEjhr35xDGaDAfet5n4kFDmz0YBty0tw\num8MT+1t1zsOXQcLd5L78f4O7D83hK/csxS5WWa941CSW1aag62LC/GNV5rRMzKpdxy6BhbuJNY5\nOIGv7ziN9yxy4cGGcr3jUAoQEfzDfcsQUgr/+JuTeseha2DhTlL+YAiff+ZtGA2Cr39gBVdIUtQs\nKLDiz7bW4aUTfXiJs0wSEgt3knr85WYcOT+C//vBlSjjqe0UZZ/aVINlpTn48vPHMeD16R2H3oWF\nOwntbnHj22+cxYfXVWDbihK941AKspgM+MYfr8bYVABffv44lOIxZ4mEhTvJdA5O4PNPv41FRXZ8\n9Z6lesehFFZfnI2/vGMRfnfiAp49zGPOEgkLdxLx+gL41A8aEVLA9kcakGXh5vcUW/9zUw3WVRfg\nK79swpm+Mb3j0IxZ9yqhxBAMKXzxp0fQ6vbi+59YhyqnTe9IlKLevZ/J1sWFONXjwUe+sw9/srn2\n0mk5D6+v0CMegSPupKCUwt/98jheOXkBX7l7CW6pc+odidJITqYZD62rwND4NH5xqAsh9rt1x8Kd\n4JRS+PqO03j6wHl8bstCfHxjtd6RKA1VO224a3kJTvZ6OEUwAbBwJzClFB5/uRnbd7XhkQ2VeOyO\ner0jURq7eaEDN9U4sKd1ALtbeDyhntjjTlDBkMJXf9WEH+/vRENlPuqLs/H0gfN6x6I0JiK4e2UJ\nxnwB7Gjqw08PduJDa9nn1gMLdwKamA7gsZ8fxYvH+/DpW2tQkW/lykhKCAYRPHhjOXz+IP762eOY\nnA6yfacDtkoSzPmhCXzgyb3Y0dSHL29bgi/dtYRFmxKK2WjAIxsqccfSInztNyfx779v4QKdOGPh\nTiAvn+jDvd/cg56RSfz3x9fiU++p0TsS0VWZjAZ86yM34P41ZXj8lWZ84adHMOUP6h0rbbBVkgAm\npgP4pxdO4ekDnVhWmoNvPXwD52lTwjMbDfjGH69CbaEd/+/lM2hzj+NbD9+ACodV72gpjyNune1q\nduPOJ3bhmYOd+MytC/H8n2xk0aakISL43JZa/NcjDWgfHMe2f9uN5w53sXUSYxxxx8C1TtK+nNcX\nwKleD55/uxs1Thue+dQGrK9xxCEdUfTdtrQIO/58E/7ip0fxFz87it+d6MPX7luGklzuXBkLLNxx\nppTC2+dH8OLxXkwHQvizrbX43JZ3lhETJavyfCuefnQDtu9qwxOvNuO2x9/AF29fhI/fXAUTD7CO\nKhbuOOr3TOHXx3rQ5h5HRYEV//WxBtQXZ+sdiyhqjAbBZzcvxN0rSvDVXzfhf//2FL675xzev7oM\nCwqu3/vm3ifasXDHgS8QxM7T/djTOgCLyYD7VpViXXUBizalrAqHFf/98bV4qakPf/3sMXz7jbNY\nU5GH25YUIc9q0Tte0mPhjiGlFE70ePDb470YnfTjxop83Lm8GPYM/rNT6hMR3LWiBH2jU3jtTD/e\nOjuIY12j2FjrxK2LXGwPzgMrSIwMeH34zdEetPR7UZKbiYfWLkCl48rZIlpuYhIluwyzEXctL8GG\nGgdePXkBu5rdONg+hC31hVhXXQAz+99zxsIdZaMTfrzU1Ic3zw7AZBDcs7IE66sdMBq4+pHSW77V\nggcbFmBjrRMvNfXht8d7savZjU11Tqyr5oyquZBYzLdsaGhQjY2NUf+8iWzKH8RTe9vxH6+fhWfS\nj9UL8nDn8mLkZJr1jkaUkNrcXrx2ph9t7nFkmg34yPpKPLKhMm3XMYjIIaVUg5bHcsQ9T4FgCM8e\n7sK/vNKCPs8UNte7sKIsl/NXiWZR47KjxmVH59AE9p4dwPf3tuO7e85h1YI83LeqFFsXF6LKwQ3W\nroYj7giNTvjxzMFOfH9vO3pGp7B6QR7+5q7F2FDjYO+aKAK3LSnE829341dHenCy1wMAKMvLwvrq\nAiwtzcHS0hxUOmzYebofhlmKeTJOLeSIO0am/EHsaRnA80e68crJC5gOhHBTjQP/9P7l2Lq4kCMD\nonkozMnEp29diE/fuhDtA+PY3TqAPS1u7GkdwHNvv3PKvEGA3Cwz8q0W5FstyLOaZ94syMsyp8V0\nQ02FW0TeB+BfARgBfEcp9c8xTZUgJqeDaOoZxdudw9h7dhBvnR2ELxBCgc2Ch9dV4MGGciwrzdU7\nJlHKqXLaUOW04ZENlQAA95gPZ/rGcH54Ai819WF4YhrD49NovjCGMV/gio81CPDU3nOoLbRfeqsr\nzEZtoT3O5fufAAAGLUlEQVRlpiDOWrhFxAjgWwBuB9AF4KCI/FopdTLW4WJJKQXPZABurw/uMR8G\nvD50j0yiY3Ac7QMT6BgcR8/o1KXH1zht+PC6CmxZXIibFzo4hYkojlzZGXBlZwAA3t3dDQRDGJ30\nY2TSj5EJPwa9PmSYDWjt9+LVU/0IhsIfYJDwfwiLi7OxqCgb9UXZqC/ORkWBNemW5M/a4xaRmwB8\nTSl158yfvwQASqmvX+tjotXjDoUUQkohpDDzq0IwpBAKAVOBICang5gKBDHlD2HKH7z05pkMYHTS\nD8+UP/zrZPjX8PvCfzc64cd0MPQHX9NmMcJhz4DDZkGB3YKy3CyUF1i5aIYoiVzscU8HQugYHMeZ\nC2No7hvDmQtjONM3ho6hiUv/ARgNguKcTJTlZ6E8Lwtl+Vlw2CzItZqRm/XOW4bJCLPRAJNRYDYa\nYDEaYDYKjAaJSps02j3uMgCXH3bYBWB9JMFms+ofXsbEdCBcnKNwz1QEyMkM/6PnZJmQm2VGcW4m\ncjLDfTCn3QJXdgac9gzsPzeE3Ewzsiyp8aMUEQEWkwF1RdmoK8oGVr7z/snpIFr6x3C6bwydgxPo\nGp5A98gk9rUNos8zNaf6IxI+0s0ggMuegb1fem/0L+RdojaMFJFHATwKwAnAKyJnovW5deIEMKB3\niHniNSQGXkOcfeTq7475NbQCkL+N+MMrtT5QS+HuBrDgsj+Xz7zvCkqp7QC2i0ijUqpKa4BENXMd\nmn5sSVS8hsTAa0gMqXANF2npyB8EUCci1SJiAfAQgF/HNhYREV3LrCNupVRARP4UwO8Qng74PaXU\niZgnIyKiq9LU41ZKvQjgRY2fc3vkcRJKKlwHryEx8BoSQypcA4AYLXknIqLYSa5Z50REFHnhFpH3\nicgZEWkVkb+5yt8vFpG3RMQnIo/NL2ZsaLiGj4jIMRE5LiJ7RWSVHjmvR8M1/NHMNRwRkUYRuUWP\nnNcz2zVc9ri1IhIQkQfimU8rDc/FZhEZnXkujojIV/XIeT1anouZ6zgiIidE5I14Z5yNhufhf132\nHDSJSFBECvTIGjGl1JzfEL5JeRZADQALgKMAlr7rMYUA1gL4PwAei+TrxPJN4zXcDCB/5vd3Adiv\nd+4IrsGOd1piKwGc1jv3XK/hsse9hvC9lgf0zh3hc7EZwAt6Z53nNeQBOAmgYubPhXrnjuT1dNnj\n7wXwmt655/oW6Yh7HYBWpVSbUmoawDMA/ujyByil+pVSBwH4I/wasablGvYqpYZn/rgP4TnsiUTL\nNXjVzCsUgA1Aot3UmPUaZvwZgGcB9Mcz3BxovY5EpuUaHgbwnFKqEwh/n8c542zm+jx8GMDTcUkW\nRZEW7qstgy+bf5y4mus1fBLAjpgmmjtN1yAi94vIaQC/BfA/4pRNq1mvQUTKANwP4D/imGuutL6e\nbp5pXe0QkWXxiaaZlmtYBCBfRF4XkUMi8rG4pdNG8/e1iFgBvA/hAUFS4c5JGojIFoQLd8L1h7VQ\nSj0P4HkReQ+AfwJwm86R5uoJAH+tlAol+Z7nhxFuMXhFZBuAXwKo0znTXJkA3AjgvQCyALwlIvuU\nUs36xorIvQDeVEoN6R1kriIt3JqWwSc4TdcgIisBfAfAXUqpwThl02pOz4NSapeI1IiIUymVKPtO\naLmGBgDPzBRtJ4BtIhJQSv0yPhE1mfU6lFKey37/oog8mYTPRReAQaXUOIBxEdkFYBWARCncc/me\neAhJ2CYBEPHNSROANgDVeOcGwLJrPPZrSMybk7NeA4AKhPeNuVnvvPO4hlq8c3PyBoRfxKJ39khe\nSzOPfwqJeXNSy3NRfNlzsQ5AZ7I9FwCWAPj9zGOtAJoALNc7+1xfTwByAQwBsOmdOZK3iEbc6hrL\n4EXkMzN//20RKQbQCCAHQEhEvoDw3V3PNT9xHGm5BgBfBeAA8OTMaC+gEmiTGo3X8EEAHxMRP4BJ\nAB9SM6/cRKDxGhKexut4AMBnRSSA8HPxULI9F0qpUyLyEoBjAEIIn4jVpF/qK83h9XQ/gJdV+CeH\npMOVk0RESYYrJ4mIkgwLNxFRkmHhJiJKMizcRERJhoWbiCjJsHATESUZFm4ioiTDwk1ElGT+P+Yx\nO1h4pnN9AAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x121be0ba8>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"sns.distplot(beta_star[:,1])"
]
}
],
"metadata": {
"jupytext": {
"cell_metadata_filter": "all,-slideshow"
},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment