{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### Programme inspecteur modèles BNP Paribas - module Monte-Carlo methods - 12 décembre 2019" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# TP 1\n", "\n", "## Warm-up: Simulation of random variables with Python and fundamental limit theorems\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 1. Simulation of a probability density with the rejection method\n", "\n", "We consider the probability distribution with support $[-2, 2]$ and density function\n", "$$\n", "f(x) = \\frac 1{2\\pi} \\sqrt{4 - x^2} 1_{x \\in [-2,2]} ,\n", "$$\n", "usually called the Wigner law.\n", "\n", "$\\blacktriangleright$ Implement a simulation method allowing to draw i.i.d. samples from the Wigner law. Simulate a large number of samples, plot the histogram and compare with the density $f$.\n", "\n", "_Reminder on the rejection method_: in order to simulate a random variable with density $f$ of support $[x_0, x_1]$ and such that $0\\leq f\\leq M$, one can use a sequence of independent couples of random variables $(Z_k, U_k)_{k\\geq 1}$ such that, for every $k$, $Z_k$ and $U_k$ are independent and uniformly distributed respectively on $[x_0, x_1]$ and $[0, 1]$.\n", "\n", "If we set\n", "\n", "$$ \\nu = \\min\\{k\\geq 1 : \\, f(Z_k) > M U_k \\}, $$\n", "\n", "the random variable $Z_\\nu$ is distributed according to the density $f$." ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEICAYAAACzliQjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4FNX6wPHvmxASeglB6QFphmKA0KSqSBEEC6iIFBURkYtiuRfwioh6xY5eFMWGqAiIBUT8KRfpRQioIE2KCBGU0AIBAiQ5vz/OJC5LygaSzG7yfp5nn+z0N7Oz7545c+aMGGNQSilVOAS5HYBSSqn8o0lfKaUKEU36SilViGjSV0qpQkSTvlJKFSKa9JVSqhApUElfRBJFpJYfxBEpIkZEirgdi/IPzvFQ24XtVne+F8FZzJNnsfmyfX9ab17L7dxwIZ+d3yZ9ERktIvO9xm3PZNxtAMaYksaYXfkZp78QkcUikiQi1TzGdRKR3XmwrSkisk1EUkVkUG6v39lGtIisE5GTzt/oDOYpKiJbRSQuL2IoCIwxe5zvRQqkHyeD3dr+hRKR3SLSKbfXm9e84/YHfpv0gaVAm7RfchG5FAgBmnqNq+3M6wo/K82fAB7Ph+38DAwD1l/MSpwE1DGD8UWBOcBHQDngA2COM97To8CBi4lBXRw/O/6VL4wxfvkCigIngWbO8C3A+8ASr3E7PJYxQG3nfTjwFXAMWAs8DSz3mncosB04ArwOiMf0u4AtzrRvgRpey97vLPtbBrFHOvMUcYbvdNZ1HNgF3Osx7xLgZud9W2e565zhTsBPPu6vxcATzjZqeyy/Ow8/o+XAIK9xQcAoYCdwCJgFlM8i5o4ZjO8M/OH1eewBunoM13T2aTcgLosYKwDzgKPAYWAZEORMS4vzOLAZuNFjuUHACuAVZ9ldwJXO+L3YH5uBHvNPBd4EFjjrW5LBMZP2uYQCLzr/01/OcsWyi9fr/3oS+K/zPgT7g/+8M1wMSML+YKYfi8AzQIozLRGY5Mt3wWu744DZ2B/kY8DgrD5zzv8ulAHeBfY7n/HTQLDH+u/h7+/KZqAp8CGQCpxy4v5nBuutDMx19tkO4B6vmGcB05z1bgJisjhmDLZQs92Z/yngMmCV8z/PAop6zN8D+Mn5zFYCjZ3xWcU90Pn8DwKPeawrFJgI7HNeE4FQj+mPOvtuHzZHpR9XPn9v8yoh5FJSWQSMdN5Pcv7JZ7zGvZfJF2uG8yoORGG/qN5Jfx5QFqgOxOMkFeAG58C5HPtl+Tew0mvZBUB5nC+rV9zeB2R356ARoAP2x6ypM208f395x2C/OM95THvVx321GPsFfBn4yBmXZdIHNjgHakavN3zYZkZJ/0FgNVDVOYDfAj7JIuaMkv5I4BuvcfOAh72GbwQ6knXSfxabVEOcVzuchAb0wSaLIOBWbOKs5EwbBCRjf7CDsclpDzYhhmJ/mI4DJZ35pzrD7Z3pr2ZwvKUdmxOxCao8UApbOHk2u3i9/q+rgY3O+yud4+YHj2k/Z3IsLgYGe60r0+9CBtsdB5zFfkeCsD8wmX7mGWz/S2d6CaAisAanEOR8Hn8AzbHfldo4P5zAbqBTFt+xJcAbQBgQ7fwP13jEnARc53yWzwKrszhmjPP5lAYaAKeBhUAt7I/WZpwffOyP0gGgpbPugU6sodnE/baz765w1n+5x3d+tbNvIrA/Ik8507piCwkNnf03nQKY9McBXzjvfwbqOP+457iB3l8sZ+efBep5TMuopN/WY3gWMMp5/w1wt8e0IGyiruGx7NVZxH3OAZnB9C+BB5z31wAbnPf/h03cqz0O5Jt83FeLnWUjgATnYHWjpL8l7cvmDFdyPovz9gWZJ/3HgRle4z4GxjnvbwT+z3nfkayT/nhsVVG2Xwxsaa2X834QsN1jWiPnM73EY9whINp5P9UzZqAktlRdzevYFOyPy2Ue87bGOWP0NV7+Ls2HY0vZY4A4Z7tPAq9ldCySedLP8LuQyXdyqa+fOeeeaVyCTXDFPObtCyxy3n+L873IYLu7ySTpA9WcfV3KY/qzwFSPmP/nMS0KOJXFvjVAG4/hdcC/PIZfAiY67yfjJGWP6duADtnEXdVj3BrgNuf9TpwzfWe4C853GHgPmOAxrS4XkPT9uU4fbF19WxEpB0QYY7Zjf/mudMY1JOP6/AjswbDXY9zeDOb70+P9SewXBqAG8KqIHBWRtNNsAapks74MiUg3EVktIoed9V2HPY0He8pYV0QuwZZQpgHVRKQC0CKT/y9Txph47BnQ+Jwsl4tqAF947Lst2C/kJQBp451pbYF5HuNGOetIxJayPJUGjotICeB54B8+xvMC9qztOxHZ5bENRGSAiPzkEU9D/v5cwJaq0pwCMMZ4jyvpMZx+TBhjErHHTWWveCKwZ5/rPLb7f874LOP1ZIw5BcRizxzbYwsIK4E2zrglmeyPzGT2XciI97Gf5WfuNV8IsN9j3rewpVqwyXtnDuMGu48PG2OOe4z7nXO/r97/X1g21yO8P+fMPvcawMNex3U1zv/cvWW2vys7sXv+H5U9pu31mpZj/p70V2FPp4Zg61cxxhzD1mcNAfYZY37LYLl47Kl5VY9x1TKYLzN7saecZT1exYwxKz3mMb6sSERCgc+wdbiXGGPKAvOxPyIYY05iSxIPAL8YY85gv7wPATuNMQdzEHeaF4CrgGbZxLbJafaW0evNC9gu2H3XzWvfhRlj/gDwHI89U+jhMW6Cs45NQGMREY/1NnbG18GWlpaJyJ/A50AlEflTRCK9gzHGHDfGPGyMqQVcDzwkIteISA3sKfZwINyJ5xecz+UCebacKomtvtnnNc9BbNJo4PF/lzHGlMwq3ky2twRbldMEe91qCbZkmFVhwafjNhve68jyM/ea7zRQwWO+0saYBh7TL7uAuPcB5UWklMe46tiqory2F3jG638vboz5xJme0/29D/tDkqY6fx9D+zk3j1W/kID9Oul7lGYewl7QSrPcGZfhgW1sM67PgXEiUlxE6gMDcrDpN4HRItIAQETKiEifC/gXwF6QDsX5IRKRbtj6YE9LsMknrXS22GvYs31vZHYbNMYcxZ6C/jOb+RoY2+wto9fQzJZzmkqGYRNkiIiEiUjasfQm8IyTVBGRCBHplV3MXhZjS4ojRCRURIY747/HJuZq2LOiaGyV1l/O+/POvkSkh4jUdn5AjjnrTcHWiRrs54KI3Ikt6V+M60SkrdPK6ClsHfs5MRljUrE/Nq+ISEVn21VEpEs28WZkCfa43uwUFhZj98dvzhlfRv7C1k3nJp8+c2PMfuA74CURKS0iQSJymYh0cGZ5B3hERJqJVTttnVnF7ezjlcCzzrHYGLgbWyWY194GhopISyfmEiLS3eMHKKf7+xPg384+rACMxV40B1vtNkhEokSkOLbhRo75ddJ3LMGe/i33GLfMGZdV1cdw7FnCn9ir6J9gSxnZMsZ8ATwHzBCRY9hE0y3Hkdt1HQdGYD+wI8Dt2ItEnpZgL+gtzWQYbKL7Hd9LL6+SebK4WN9hS6tXAlOc9+09tjsXWz1xHHtRqmVOVu4ksBuwCe0o9gL+DcaYM8aYZGPMn2kvbBVKqjOc0f9bB/gftspoFfYC9WJjzGbsD+Mq7BezEc7Z5EWYjv0iHsaeZfXLZL5/YatwVjvH1/+AelnFm8l6VmLr9tOOk83Yev6svhevAr1F5IiIvObD/+SLnHzmA7AFoc3Y78Ns7DUAjDGfYhtqTMdeFP8Se7YEto7+304VyiMZrLcv9gxwH/AF8IQxZsFF/2fZMMbEYlscTcL+Pzuw14PSZBe3t6exBd0NwEZss+innW19g20E8L2zne8vJOa0VgwFnog8B1xqjBnodiwXQkT+DcQbY95yOxZ1PhGZir2g/G+3Y/EnYu+Q3469kFw4ko2fK7A3VjhVOkWxv5bNsad7+XYnYm4zxjztdgxKXYCG2NYnmvD9RIFN+tjqkU+wV7wPYE/l57gakVKFiIg8hL2u5GtLK5UPCk31jlJKqcC4kKuUUiqX+F31ToUKFUxkZKTbYSilVEBZt27dQWNMRHbz+V3Sj4yMJDY21u0wlFIqoIiIT3foavWOUkoVIpr0lVKqENGkr5RShYjf1eln5OzZs8TFxZGUlOR2KEoBEBYWRtWqVQkJCXE7FKVyJCCSflxcHKVKlSIyMpJzO15UKv8ZYzh06BBxcXHUrFnT7XCUypGAqN5JSkoiPDxcE77yCyJCeHi4nnmqgBQQSR/QhK/8ih6PKlAFTNJXSil18TTp+yg4OJjo6GgaNGjAFVdcwcsvv0xqamqubuPNN99k2rRpAEydOpV9+7wfunRhli1bRoMGDYiOjubUqVPp448ePcobb7yRPrx48WJ69OhxQds4evQo4eHhac/uZNWqVYgIcXFxACQkJFC+fHlSU1MZO3Ys//vf/y7iP8q5QYMGMXv27HzdplL+yKcLuSLSFfughGDgHY/H2qVNHwrcj31oRyIwxBiz2XnK0xbsg4LBPvA70ycy+bNixYrx008/AXDgwAFuv/12EhISePLJJ3NtG0OH/r1rpk6dSsOGDalcObtHbWbv448/5pFHHuHOO+88Z3xa0h82bNhFb6Ns2bJceumlbNmyhaioKFauXEmTJk1YuXIlt9xyC6tXr6Zly5YEBQUxfnzeP743JSWF4ODgPN+O2yJHfZ3ptN0TuudjJCpQZFvSF5Fg4HXsk6OigL4iEuU123RjTCNjTDT2odUve0zbaYyJdl4BmfC9VaxYkSlTpjBp0iSMMaSkpPDoo4/SvHlzGjduzFtv2eecLF68mI4dO9K7d2/q169Pv3790kvCo0aNIioqisaNG/PII/aBOuPGjePFF19k9uzZxMbG0q9fP6Kjo/n666+58cYb07e/YMECbrrppvPiWrhwIU2aNKFRo0bcddddnD59mnfeeYdZs2Yxfvx4+vU790FOo0aNYufOnURHR/Poo48CkJiYmGG869ato0OHDjRr1owuXbqwf//+87bfpk0bVq60jxFeuXIlI0eOPGf4yiuvBM4tdc+fP5/69evTtm1bRowYkX6mMW7cOO666y46duxIrVq1eO21vx/y9NFHH9GiRQuio6O59957SUmxD8wqWbIkY8eOpWXLlqxatSrTz2/8+PE0b96chg0bMmTIEIwxHDhwgGbN7COFf/75Z0SEPXv2AHDZZZdx8uTJTNenVCDxpaTfAthhjNkFICIzgF7Yx50B6Q8rT5P27NG88eCD4JS4c010NEycmKNFatWqRWpqKgcOHGDOnDmUKVOGtWvXcvr0adq0aUPnzvYxuD/++CObNm2icuXKtGnThhUrVhAVFcUXX3zB1q1bERGOHj16zrp79+7NpEmTePHFF4mJicEYw8MPP0x8fDwRERG8//7755Xak5KSGDRoEAsXLqRu3boMGDCAyZMn8+CDD7J8+XJ69OhB7969z1lmwoQJ/PLLL+lnMIsXL84w3pYtW/KPf/yDOXPmEBERwcyZM3nsscd47733zlnflVdeydKlSxk8eDC7du2iT58+6T+AK1euZPTo0efFfO+997J06VJq1qxJ3759z5m+detWFi1axPHjx6lXrx733XcfO3bsYObMmaxYsYKQkBCGDRvGxx9/zIABAzhx4gQNGzbM9kxi+PDhjB07FoD+/fszb948rr/+epKSkjh27BjLli0jJiaGZcuW0bZtWypWrEjx4sWzXKdSgcKXOv0qnPvA6Thn3DlE5H4R2Ykt6Y/wmFRTRH4UkSUi0u6iovUzaaXg7777jmnTphEdHU3Lli05dOgQ27dvB6BFixZUrVqVoKAgoqOj2b17N6VLlyYsLIzBgwfz+eefZ5tQRIT+/fvz0UcfcfToUVatWkW3buc+snfbtm3UrFmTunXrAjBw4ECWLs3qUakZyyjebdu28csvv3DttdcSHR3N008/nV5X7ymtpP/bb78RGRlJWFgYxhgSExNZt24dLVq0OGf+rVu3UqtWrfS27t5Jv3v37oSGhlKhQgUqVqzIX3/9xcKFC1m3bh3NmzcnOjqahQsXsmvXLsBed7n55puz/R8XLVpEy5YtadSoEd9//z2bNm0C7I/WihUrWLp0KWPGjGHp0qUsW7aMdu0K1GGrCjlfSvoZtU07ryRvjHkdeF1Ebgf+DQwE9gPVjTGHRKQZ8KWINPA6M0BEhgBDAKpXr551NDkskeeVXbt2ERwcTMWKFTHG8N///pcuXbqcM8/ixYsJDQ1NHw4ODiY5OZkiRYqwZs0aFi5cyIwZM5g0aRLff5/1M47vvPNOrr/+esLCwujTpw9Fipz70eXWw3AyitcYQ4MGDbKsMgGoU6cOR44c4auvvqJ169YANGvWjPfff5+aNWtSsmTJHMWcWSwDBw7k2WefPW/+sLCwbOvxk5KSGDZsGLGxsVSrVo1x48alt7dv164dy5Yt4/fff6dXr14899xziMgFX9xWyh/5UtKPA6p5DFfFPnE+MzOAGwCMMaeNMYec9+uAnUBd7wWMMVOMMTHGmJiIiGy7g3ZdfHw8Q4cOZfjw4YgIXbp0YfLkyZw9exaAX3/9lRMnTmS6fGJiIgkJCVx33XVMnDgxvXrFU6lSpTh+/Hj6cOXKlalcuTJPP/00gwYNOm/++vXrs3v3bnbs2AHAhx9+SIcOHbL8P7y3kZl69eoRHx+fnvTPnj2bXjr21rp1a1599dX0pN+6dWsmTpyYXp/vHfOuXbvYvXs3ADNnzsw2lmuuuYbZs2dz4MABAA4fPszvv/vUoyxAeoKvUKECiYmJ57Toad++PR999BF16tQhKCiI8uXLM3/+fNq0aePz+pXyd76U9NcCdUSkJvAHcBtwu+cMIlLHGLPdGewObHfGRwCHjTEpIlILqAPsyq3g89OpU6eIjo7m7NmzFClShP79+/PQQw8BMHjwYHbv3k3Tpk0xxhAREcGXX36Z6bqOHz9Or169SEpKwhjDK6+8ct48gwYNYujQoRQrVoxVq1ZRrFgx+vXrR3x8PFFR3tfRbSn3/fffp0+fPiQnJ9O8efNzWgNlJDw8nDZt2tCwYUO6detG9+4Zt/YoWrQos2fPZsSIESQkJJCcnMyDDz5IgwYNzpu3TZs2zJ8/n5iYGMAm/V27dmWY9IsVK8Ybb7xB165dqVChwnnVPxmJiori6aefpnPnzqSmphISEsLrr79OjRo1sl0WbCuje+65h0aNGhEZGUnz5s3Tp6U9vKd9+/YAtG3blri4OMqVK+fTupUKBD49I1dErgMmYptsvmeMeUZExgOxxpi5IvIq0Ak4CxwBhhtjNonIzcB4IBnbnPMJY8xXWW0rJibGeD9EZcuWLVx++eU5/+8KmOHDh9OkSRPuvvtut0PJNYmJiZQsWRJjDPfffz916tRh5MiRboflE384LrXJpkojIuuMMTHZzedTO31jzHxgvte4sR7vH8hkuc+Az3zZhspas2bNKFGiBC+99JLboeSqt99+mw8++IAzZ87QpEkT7r33XrdDUqpAC4heNpVtJ18QjRw5MmBK9koVBNoNg1JKFSKa9JVSqhDRpK+UUoWI1ukrV22IO5rl9MZVy+ZTJEoVDgGZ9LNqpnYhfGnaFhwcTKNGjdKHb7vtNkaNGnXR277uuuuYPn06Zcuen9yympaVqVOnEhsby6RJk3war5QqPAIy6bvBs2vl3DR//vzzxhljMMZkOE0ppS6G1ulfpMjISMaMGUPr1q2JiYlh/fr1dOnShcsuu4w333wTsH3wtG/fnhtvvJGoqCiGDh2a/gCWyMhIDh48yO7du7n88ssZNmwYTZs2Ze/evenTAKZNm0bjxo254oor6N+/PwBfffUVLVu2pEmTJnTq1Im//vrL57gzW7ZRo0YcPXoUYwzh4eHpD3Xp379/vj/4RCmV+7Sk76O0bhjSjB49mltvvRWAatWqsWrVKkaOHMmgQYNYsWIFSUlJNGjQIL0rhDVr1rB582Zq1KhB165d+fzzz8/r6njbtm28//775zzNCmDTpk0888wzrFixggoVKnD48GHAdhOwevVqRIR33nmH559/3uebtzJbNq075Ro1alCrVi2WLVvGgAEDWL16NZMnT77g/VcY6N2xKhBo0vdRVtU7PXv2BGwpOTExkVKlSlGqVCnCwsLS+8pv0aIFtWrVAmwXwsuXLz8v6deoUYNWrVqdt/7vv/+e3r17U6FCBQDKly8PQFxcHLfeeiv79+/nzJkz6V0U+yKzZdu1a8fSpUupUaMG9913H1OmTOGPP/6gfPny5/WSqZQKPFq9kwvSugAOCgo6pzvgoKAgkpOTAdsnvifvYYASJUpkuH5jTIbz/+Mf/2D48OFs3LiRt956K70HSV9ktmz79u1ZtmwZy5Yto2PHjkRERDB79mztU16pAkJL+vlkzZo1/Pbbb9SoUYOZM2cyZMgQn5e95ppruPHGGxk5ciTh4eEcPnyY8uXLk5CQQJUq9nk2H3zwQY7iyWzZatWqcfDgQc6cOUOtWrVo27YtL774ol+2+MmquWdWTT0vdDmlCoKATPpu1I961+l37dqVCRMmZLHEuVq3bs2oUaPYuHFj+kVdXzVo0IDHHnuMDh06EBwcTJMmTZg6dSrjxo2jT58+VKlShVatWvHbb7/5vM6slm3ZsmX6c2fbtWvH6NGjadu2rc/rDgRiDEVSUwhJTaZIagrBqakEp6bAHycgNRVSUuzf1FQwxr7OWYFAfDw88wyULAmlSjFi9X4SwkpypFgpjoaVIr5kOQ6UKM/h4qWzjCW7Jsh6PUDlJp+6Vs5PBbFr5cWLF/Piiy8yb948t0PxOxdzc1a2JfazZ+H0aUhKsq/Tp+HMGc6eSiIkNSXzjQYFQXCw/RsUZBO8d/WaMWzZv5/L77kHTpyAxETI5OHpZ4OCCalRHWrUgMhIqF0b6tSBunWhXj0in8z6qWlZJX29eKzS5GrXykr5u9DkMxQ7e5piyacJO3sGDv1uk34aEShaFEJDOR5agrPBRTgbVITk4GCSg+wrJSiYBlXLnZ/gMyMCu/5+JtBl/5xLmaREyp06RrlTx4hIPELFE0e4JPEQlY/FU2X7AarFbuDSxMPpy6QifF+uEr9G1GBLRE02VKrDL5fUJr6kPrhF5Q1N+vmgY8eOdOzY0e0wCo6UFEhM5JLEwxQ/k0Txs0kEG3vfQ6oISUWKQunSUKyYfYWG2peTzOOyOrvwNeFnFFZQMIeLl+Fw8TJZzhd2NonII/updfgP6h78nbrxv1Pv4O90/nU1Qc7jp/eVqsD6KpezvnJ9iL0EoqOhiH5d1cULmKMosxYsqhBITbXVJ8eOwfHjtjoFqAgkFQnlaLFSnAwJ41RIKKeLFMWQ9xdkL6ZaNCkkjK0Va7K1Yk3m8/e1khKnTxJ1YBeN/tzJFft/pdkfW+ixdRk0f9v+iLVrB1ddBV26QIMGF/UDpQqvgEj6YWFhHDp0iPDwcE38hURISjKlT5+A7fE20aem2iRXogRUqgSlSrEpIZlUyd1Wx75cYzDGcOjQIcLCwnJ12ydCi7O2WkPWVmuYPu6S4wf5oX0YLF4MixbB11/DI49A5crQpQudj1ZmWWQTThXN3VhUwRUQSb9q1arExcURHx/vdigql/115FT6+yKpyen18qQkcww4VqSIraIJC7MvY2yJ/9gx9nss623L8WI+bTOn0tYbFhZG1apVL3g9vvqrVAW4tTs4d3+zdy989x18+y18/jlTEhI4HRzCihpX8E29NnxbtzXHwvQmOpW5gEj6ISEhObrbVOWNvGgpcvew9+m1eQk9Ny/h8vjdAKyvXI+5dVqxoHZL/vf20EyrMbpdYDxZLZcd11vEVKsGd99tX2fP0nfAi3Ta8QNdfl3F1bti+c+3k1geGc2cqI6Q2ME2J1XKg09JX0S6Aq8CwcA7xpgJXtOHAvcDKUAiMMQYs9mZNhq425k2whjzbe6FrwJSYiJ8+il88AErlywBYF3l+oy/+h6+qXcl+0tH/D2vVudl3Y6/RmNW1WjMU1cPptGfO+i+bTk9tixj4ryX4JLJcMMNcOedcPXVtvmpKvSyTfoiEgy8DlwLxAFrRWRuWlJ3TDfGvOnM3xN4GegqIlHAbUADoDLwPxGpa4zJopG0ulh+23Z77Vp46y2YOdMm/rp1ebHdHcyJ6sjespe6F1dBIMLGSnXYWKkOz3UYSLM/tjC75C67r6dPt/cH3HWXPUOoXNntaJWLfPnpbwHsMMbsMsacAWYAvTxnMMYc8xgsAaQ1begFzDDGnDbG/AbscNanCotTp+D996F5c2jRAmbMgFtugeXLYetWJl15myb8XGYkiNiqDeDNN2H/fpv0L7sMxo6F6tWhTx97UdjPbsxU+cOX6p0qwF6P4TigpfdMInI/8BBQFLjaY9nVXstWuaBIVWDZvx/eeAMmT4ZDhyAqCl5/He64wzY/DFC5/dS2PBcWBn372tfOnfaH4L33YPZsaNQIRo6E22+39zGoQsGXpJ9Rpep5RQRjzOvA6yJyO/BvYKCvy4rIEGAIQPXq1X0ISfmtrVvh+efho48gORl69oQHH4QOHbR+Pp9l+AMV3JHQQa3Z1jABXnnFVvmMHg0jRsCwYZDDR3OqwONL0o8DqnkMVwX2ZTH/DCDtaRs+LWuMmQJMAdv3jg8xqQuUZ/X969fDf/4Dn39uS5dDhthkX7v2ha9T5YnTIaH24u6gQbBwIbz0Ejz2GDz3nE38I0dCxYpuh6nyiC9Jfy1QR0RqAn9gL8ze7jmDiNQxxmx3BrsDae/nAtNF5GXshdw6wJrcCFz5h6i/dvHgiunw3GooUwbGjIEHHoCIv1vguFElEnDVMG4QgU6d7OvHH2HCBJv4X3sNhg+HRx8F58E9quDI9kKuMSYZGA58C2wBZhljNonIeKelDsBwEdkkIj9h6/UHOstuAmYBm4H/A+7XljsFw2WH9vLGF/9h/tQRtNqzEZ58En7/HZ5++pyErwJEkya2pc+WLbaZ5wsv2BY/jz9ub4ZTBUZAdK2scuZCS7nZVe9EjvqaS44f5IEVn3DLhgUkhYTybswNvNu8Fxsm3prr8Sj37B5YC8aNg1mzIDxua5wTAAAgAElEQVTcVv8MG6YXfP2Yr10r690ayjcnTvDg8o9ZPOVeem9cyLSmPWh/7zu80q6f3vZfAEV+sIvImgPoMXAiS0tVh4ceYs8lNbjvhtHa1DPAadJXWTPGtsSpV48HV3zCwtotuPqeNxnfaUi2XQirwPfLpbUZcOtT3HHLU5wMCWPynAm2Jdb69W6Hpi6QJn2VuY0boX176N8fKlWid7/nGN7rX8TpzVSFzvKaTeh+52uM6XK/rfePiYH774cjR9wOTeWQJn11vuPH4aGH7MW9LVvgnXfghx/sXZ6q0EoJCmZ6dDfYvt227nnzTahXD6ZN0yqfAKJJX53r66/tAzomToTBg+HXX21/LdpZl0pTtqxt1rlunb0PY+BA6NzZ3vGr/J5+kxUA5U8m2Fv1e/Sw3SSsWGFLcuXLux2a8lfR0bYPpcmTYc0a263DCy/Yx1kqv6VJX9Fl20q+e3cYfPYZjB9vL9K1bu12WCoQBAXB0KGwebN9jOM//wlt28K2bW5HpjKhSb8QK52UyKtzX+CtL//D/lIVbLJ//HEoWtTt0FSgqVLFdsExfbpN+NHRtm+f1FS3I1NeAuLJWSr3tdqzgZfnvUzEiSO81LYfk1v1YUfDhtkvqAq9rPtv6gsdO8K999rGAN98A1Onah/+fkRL+oVMkZRk/rlkKtM/eYykIkW56Y4X+W+bviQH6++/yiWVKsGcOTBlir021LgxfPml21Ephyb9QqRKwgFmTf8Xw1bPZsYVnek+6DU2VqrjdliqIBKBe+6xVYaRkXDjjbb75tOn3Y6s0NPiXSHRafsPvDj/FYJTUxjWaxTz67c9bx7tI0flunr1YOVKGDXK1vGvXGn786lVy+3ICi1N+gVccGoKjyz9kPt+mM0vl1zG/b3+xe/lcr9+VX8wVKaKFoWXX7Z3d995JzRtCh9+CNdf73ZkhZJW7xRg5U4mMHXWE9z3w2w+ju7KzXe8kCcJXymf3HCD7bf/ssvsE9WeeEJb97hAS/oFVIO/dvLW508TceIoj3YbwaeNO7sdklK2fn/5crjvPntPSGysbeZZRjvvyy+a9ANUVtUpXbat5JWvX+JoWCn69HuODZXq5mNkSmWjWDF4/31o0cI+Za11a/jqK3sGoPKcJn0/dUF15MYwfNVMHln2ET9WqseQm/5NfMlyuR+cUhdLxD6UJSoKbr7Z/gB89plt46/ylCb9AiIk5SwTvnmNmzct4ouojozqNoLTRfTOWuU/MivI1Og9gXdnj6fGNZ0Y3fUfzG7U6bx5snuqm/KdXsgtAEonJfLBrCe4edMiXmrbj5E9HtaErwLG7+Uqc+OAl1hdrREvzp/IA8una1fNeUiTfoCrdCyeTz/+J83jNjGy+0P8t01fe+qsVAA5HlqCu/o8wWcNr2bkiuk8/82rFElJdjusAkmrdwLYZYf28uHMxyl5+iQD+zzJyshot0NS6oKdDQ7h4etGElf6Eh5Y+QnlTyZwf69RnA4Jzaa/H636yQmfSvoi0lVEtonIDhEZlcH0h0Rks4hsEJGFIlLDY1qKiPzkvObmZvCFWaP925n18b8ISU3mttsnaMJXBYMIr7Trx2Odh3H1zlimzRpL6aREt6MqULJN+iISDLwOdAOigL4iEuU1249AjDGmMTAbeN5j2iljTLTz6plLcRdqrfZs4JMZYzhZtBi9+z3P5kv0lnZVsHzc5DpG9HyUJvu2MeOT0YSfOOp2SAWGLyX9FsAOY8wuY8wZYAbQy3MGY8wiY8xJZ3A1UDV3w1Rp2v22nqmfjmNfqQh693tO77BVBda8y9sz+ObHiTyyjxmfjCYi8bDbIRUIviT9KsBej+E4Z1xm7ga+8RgOE5FYEVktIjdktICIDHHmiY2Pj/chpMLpqp1reeez8ewsX5Vbb3+Wv0pVcDskpfLU0lrNGNTnSSofi2fm9FFceuyg2yEFPF+SfkZNQTJsTyUidwAxwAseo6sbY2KA24GJInLebXfGmCnGmBhjTExERIQPIRU+nX9dxVufP8O2iEhuv+0ZjhTX29ZV4bCmWkP63/IUESeOMPOTUVRJOOB2SAHNl6QfB1TzGK4K7POeSUQ6AY8BPY0x6Z1mG2P2OX93AYuBJhcRb6F09Y41TJrzHJsvqcUdtz5NQrFSboekVL5aX/Vy7rj1acqdOs70GWO45LiW+C+UL0l/LVBHRGqKSFHgNuCcVjgi0gR4C5vwD3iMLycioc77CkAbYHNuBV8YtPttPZO//A9bK0Yy4JbxHAsr6XZISrni58r1GHDLeMqfTGD6jMeISDzidkgBKdt2+saYZBEZDnwLBAPvGWM2ich4INYYMxdbnVMS+FTsjUF7nJY6lwNviUgq9gdmgjFGk76PWv++gbc/f5qd4dXof8tTmvCV38vr5yr8VLked/YZxweznuDjGY9x2+3P5un2CiIxfna7c0xMjImNjXU7DNf1HPgKn3wyhj9KV+S225/lsNbhK5Wu1Z4NvP/pk+wIr0qjX9dD6dJuh+Q6EVnnXD/NknbD4I+2bmXqp+M4VLwMd9z6lCZ8pbysrt6Y+24YRf343fbhLElJbocUMDTp+5u9e+Haa0kJCqL/rU9xoFS42xEp5ZcWX9acR657EBYtgr59IVn76vGFJn1/cuQIdOkCx44xsM94vfFKqWzMaXAVvPoqfPml7Z/fz6qr/ZF2uOYiz4teRZPPMm3W4zTZt50BtzylXSso5asRI+DPP+HZZ6FmTRg92u2I/JomfT8gJpXnv5lIq72/MOL6R/mheiO3Q1IqYESO+hrMlUyM6sANY8YwYtVh5kZ1BLQHzoxo0vcDDy/7iBs2L+H59gOYG9XB7XCUCjwi/LPbg1x6/BAvzJ/IXyXDtfCUCa3Td9mNv3zP8FWzmH5FF95o1cftcJQKWGeKhDDkpn+zt8ylvPnFf6h29E+3Q/JLmvRdFL1vGxP+77+sqt6Isdfep0+8UuoiHQsryd29xyIY3vlsPBw75nZIfkeTvlv++IO3vniGv0qWZ1ivUSQHa02bUrnh93KVGdZrFJcdioM77oCUFLdD8iua9N2QlAQ33ECJM6cYfPPj2mOmUrlsZWQ046+5B776Ch5/3O1w/IomfTcMHw6xsYzs8TC/RkS6HY1SBdK0pj3gnntsU84vv3Q7HL+hST+/vfuufY0Zw4I6rdyORqmCSwReew1iYmDgQNi+3e2I/IIm/fy0bh3cfz906gTjx7sdjVIFX1gYzJ4NRYrAzTfDiRNuR+Q6Tfr55cgR6N0bKlaE6dMhONjtiJQq8CJHfU3k5F/o3+lBUjf+wmctryfyX/PyvAtof6ZJPz8YA4MHQ1wcfPop6CMhlcpXy2o25dU2fbl50yJu/uV7t8NxlbYTzGORo76m34/zeea7z/lPxzuZ8sVB+KLwljKUcst/r7yV1ns2MH7BZNZXqe92OK7Rkn4eqxu/m8e/f4clNZvydosb3Q5HqUIrNSiYB3s8wukiRfnv3Ofh9OnsFyqANOnnpZMnmTTneY6HFufh7iMxortbKTf9WboCj173AA3/2gn/+pfb4bhCs1BeGj2auof2MLL7wxwsUc7taJRSwMLaLXm/2fW2H/4FC9wOJ99p0s8rCxfCa6/xfrPrWV6zidvRKKU8TOgwCOrXhzvvtC3rChGfkr6IdBWRbSKyQ0RGZTD9IRHZLCIbRGShiNTwmDZQRLY7r4G5GbzfSkiwB1O9ejzXoXD8y0oFktMhofDhh/bhKyNGuB1Ovso26YtIMPA60A2IAvqKSJTXbD8CMcaYxsBs4Hln2fLAE0BLoAXwhIgU/HqOBx6Afftg2jSSQsLcjkYplZGYGNsvz0cf2Ru4CglfSvotgB3GmF3GmDPADKCX5wzGmEXGmJPO4GqgqvO+C7DAGHPYGHMEWAB0zZ3Q/dTcufDBB/DYY9CihdvRKKWyMmaMTf5Dh8KBA25Hky98aadfBdjrMRyHLbln5m7gmyyWreK9gIgMAYYAVK9e3YeQ/Eva3X2lTp9gwTv3cSQikutPNiW5EN/1p1RACAmxhbToaBg5Ej7+2O2I8pwvJf2MnuyR4SPnReQOIAZ4ISfLGmOmGGNijDExEQF8t+qoxe8TceIo/+o2QvvHVypQREXZM/Pp0+Hrgl9Q8yXpxwHVPIarAvu8ZxKRTsBjQE9jzOmcLFsQtNyzkX4//R/vxfRkQ6W6boejlMqJ0aOhQQO47z44ftztaPKUL0l/LVBHRGqKSFHgNmCu5wwi0gR4C5vwPSvGvgU6i0g55wJuZ2dcgRKafIZn/++/7ClzCS+3vcPtcJRSOVW0KLz9tu0fa8wYt6PJU9kmfWNMMjAcm6y3ALOMMZtEZLyI9HRmewEoCXwqIj+JyFxn2cPAU9gfjrXAeGdcgTJs1afUOrKPMV2Gc6qottZRKiC1bg3/+Ae8/jqsWeN2NHlGjMmwet41MTExJjY21u0wfLdzJ6frXc63dVszouc/3Y5GKeWj3RO6nz/y2DGoVw+qVoXVqwOqC3QRWWeMicluPr0j92IYAyNGcDa4CM9cdZfb0SilLlbp0vDSSxAba59wVwBp0r8YX30F8+czsU1f/ipVwe1olFK5oW9f6NDBXtw9eNDtaHKdJv0LdeqUvfM2KoqpzXpmP79SKjCIwKRJtjuVAnhRV5P+hXrhBdi9GyZN0jb5ShU0DRvaQt0778D69W5Hk6s06V+I/fvh+efhppvgqqvcjkYplRfGjoXwcHj4YXv9roDQpH8hxo61T9157jm3I1FK5ZUyZWDcOFi8GObNczuaXKNJP6c2boT33oP774fatd2ORimVl4YMgbp14dFH4exZt6PJFZr0c+rRR22zrrFj3Y5EKZXXQkJsVe62bfaO3QJAk35OfPcdfPut7YO7fHm3o1FK5YeePW0TzieesDdvBThtduIrY2zzrRo1bNWOUqpAi/ToGr1xjV7MXbKEl7sP47U2fTO+mzdAaNL31Zw5sG6drc8PDXU7GqXURYrMwfMuNlSqy7d1WjF4zRd80LRHHkaV97R6xxepqbYOv04d6N/f7WiUUi54pW0/Sp85yT1rv3A7lIuiJX1ffPopbNzIiOsfYe6/C1zP0EopH2ytWJOv6rfjzti5NB0xncPFy2Q4n79X/WhJPzspKTBuHL+GV2de/XZuR6OUctHEtrcTlnyGoasD90HqmvSzM306bN3Ky+36kRoUON2sKqVy387wanzZoCMDfvyaiMTAfDSIVu84MrqoIyaVBe+M4WxEJN/Wbe1CVEopf/PqlX25YdNi7l77JRMCsEt1Leln4drtP1D7cByTW/XBiO4qpRTsKVeJr+u35Y6fvqF0UqLb4eSYZrLMGMOw1Z/ye9lL+bp+W7ejUUr5kcmtelPyzCnu+HG+26HkmCb9TLTau5Ho/b8ypcVNpGhdvlLKw5aKtVhUqxl3xc4h9Oxpt8PJEU36mbhv9WziS5RldqNOboeilPJDk1v1ocLJBPps/J/boeSIT0lfRLqKyDYR2SEiozKY3l5E1otIsoj09pqWIiI/Oa+5uRV4Xmrw5w46/Lae92J6cbpIUbfDUUr5oTVVG7Cucn3uXfM5wakpbofjs2yTvogEA68D3YAooK+IRHnNtgcYBEzPYBWnjDHRzisgnis4ZM0XHC9ajI+aXOd2KEopfyXCG637UC3hL7pvXe52ND7zpaTfAthhjNlljDkDzAB6ec5gjNltjNkApOZBjPkqIvEw121bzqzGnTkeWsLtcJRSfuz7y5qzq1xlBqwPnIes+JL0qwB7PYbjnHG+ChORWBFZLSI3ZDSDiAxx5omNj4/PwapzX7+fviEkNYVpTf37VmqllPuMBPFh0+7E/LGFBn/ucDscn/iS9CWDcTl5YGR1Y0wMcDswUUQuO29lxkwxxsQYY2IiIiJysOrcFZJyltt/+j8W1WrG7+UquxaHUipwzG7UiRMhYQwMkNK+L0k/DqjmMVwV2OfrBowx+5y/u4DFQJMcxJevum1bScUTR/ig6fVuh6KUChDHQ0vwecOr6bV5CeVOJrgdTrZ8SfprgToiUlNEigK3AT61whGRciIS6ryvALQBNl9osHltwPp5/FauEktqNXU7FKVUAJnWpDuhKWe5dcMCt0PJVrZJ3xiTDAwHvgW2ALOMMZtEZLyI9AQQkeYiEgf0Ad4SkU3O4pcDsSLyM7AImGCM8cuk3+DPHcT8sYUPm/TQLheUUjmyPaIGK2o05o4fv4bkZLfDyZJPHa4ZY+YD873GjfV4vxZb7eO93Eqg0UXGmC8Grp/HyZBQZje6xu1QlFIBaFrTHrz1xX9g3jy4IcM2K35Bi7QAiYl037qcuZd34FhYSbejUUoFoP/VbsmfJcvbR6r6MU36AJ9+SomzSXyqXS4opS5QSlAwXzS4GubPhz//dDucTGnSB5g6lV3lKrOuyuVuR6KUCmCzG11jn7b38cduh5IpTfo7d8LSpbZjNcnolgSllPLNzvBq0KoVTJ0KJie3M+UfTfoffABBQXze4Gq3I1FKFQR33gm//ALr1rkdSYYKd9JPTbVJ/9pr+bN0BbejUUoVBLfcAmFhtrTvhwp30l+0CPbsgUGD3I5EKVVQlC0LN94I06dDUpLb0ZyncCf9qVOhTBm/blOrlApAd94JR47AV1+5Hcl5Cm/ST0qCL7/8+1RMKaVyy9VXQ6VKMHOm25Gcp/Am/e++g8RE6N07+3mVUiongoPhpptsm/0TJ9yO5hyFN+nPng3lysFVV7kdiVKqIOrdG06dgm++cTuScxTOpH/mDMydC716QUiI29EopQqidu0gIgI++8ztSM5ROJP+woWQkKBVO0qpvBMcbFvxzJvnV614CmfS/+wzKF0aOmlfO0qpPHTzzfba4XffuR1JusKX9JOTbaud66+H0FC3o1FKFWRXXWWvHc6e7XYk6Qpf0l+yBA4dsr/ASimVl0JC7LXDuXPttUQ/UPiS/uzZUKIEdO3qdiRKqcKgd297DXHhQrcjAQpb0jcG5syBbt2gWDG3o1FKFQadOkGpUrZa2Q8UrqS/aRPs32+TvlJK5YfQUHuH7nff+UV3y4Ur6addQb/2WnfjUEoVLtdeC7t32+d3uMynpC8iXUVkm4jsEJFRGUxvLyLrRSRZRHp7TRsoItud18DcCvyCLFgA9etDtWquhqGUKmQ6d7Z/FyxwNw6gSHYziEgw8DpwLRAHrBWRucaYzR6z7QEGAY94LVseeAKIAQywzln2SO6EnwOnT9uWO4MH5/umlVKFR+Sor88faQzLS1ek6oIFcN99+R+Uh2yTPtAC2GGM2QUgIjOAXkB60jfG7HampXot2wVYYIw57ExfAHQFPrnoyHNqxQo4dYq795VlYUYfilJK5RURlkVG03fhQnuvUBFfUm/e8KV6pwqw12M4zhnnC5+WFZEhIhIrIrHx8fE+rjqHFizgbFAwq6s1ypv1K6VUFpZHNoFjx2DtWlfj8CXpZ/S0cF8vQfu0rDFmijEmxhgTExER4eOqc2jBAtZXrs+J0OJ5s36llMrCisgrQMT1en1fkn4c4Hnlsyqwz8f1X8yyuefgQVi/nuWR0fm+aaWUAjharDQ0a+Z6Pzy+JP21QB0RqSkiRYHbgLk+rv9boLOIlBORckBnZ1z+WrjQXkiJbJLvm1ZKqXTXXgurV9tqHpdkm/SNMcnAcGyy3gLMMsZsEpHxItITQESai0gc0Ad4S0Q2OcseBp7C/nCsBcanXdTNVwsWQNmybKhUJ983rZRS6Tp3hpQUWLzYtRB8uoRsjJkPzPcaN9bj/Vps1U1Gy74HvHcRMV68BQvg6qtJCQp2NQylVCHXujUUL25zUs+eroRQ8O/I3b8f9uyxT7FRSik3hYZCy5bwww+uhVDwk/66dfZvTIy7cSilFNhc9PPPrnW1XPCT/tq1EBQETfQirlLKD8TE2IS/caMrmy/4ST82FqKibB/6SinltubN7d/YWFc2X7CTvjF2x2rVjlLKX0RGQvnymvTzxN69cOCAJn2llP8QsTnJpe4YCnbST/slTTudUkopf9C8OfzyC5w6le+bLvhJv0gRaNzY7UiUUupvMTH2Jq2ff873TRfspL92LTRqBGFhbkeilFJ/S6tydqGKp+Am/bSLuFq1o5TyN1WqwKWXunIxt+Am/V274OhRvYirlPI/aRdzNennIr2Iq5TyZ82bw5YtkJiYr5stuEl/7Vrbz0WDBm5HopRS54uJsdXQ69fn62bde1BjXouNhehoCAlxOxKllALOfWh6+ImjrAOeGjeNd1scB2D3hO55HkPBLOkbYzta0/p8pZSfOlSiLHGlI2j854583W7BTPpHj9p6stq13Y5EKaUytadsJSodj8/XbRbMpB/v7MS8esi6UkrlgsPFyxB+Mn8fnVgwk/6BA/avJn2llB87VLw04SeP5us2C2bSTyvpV6zobhxKKZWFQ8XLUjYpkSIpyfm2zYKd9LWkr5TyY4eLlwGg3Knj+bZNn5K+iHQVkW0iskNERmUwPVREZjrTfxCRSGd8pIicEpGfnNebuRt+JtKqdypUyJfNKaXUhThUrDRAvlbxZNtOX0SCgdeBa4E4YK2IzDXGbPaY7W7giDGmtojcBjwH3OpM22mMic7luLMWHw9lytibs5RSyk8dKlEWgPInE/Jtm76U9FsAO4wxu4wxZ4AZQC+veXoBHzjvZwPXiIjkXpg5FB+vVTtKKb93qJit3gnPx6Tvyx25VYC9HsNxQMvM5jHGJItIAhDuTKspIj8Cx4B/G2OWeW9ARIYAQwCqV6+eo3/AW+Sor/lo5WaKnS3CzR53vymllL85XDytese/SvoZldiNj/PsB6obY5oADwHTRaT0eTMaM8UYE2OMiYnIhRJ6+MmE9AskSinlr44WK0WKBPld9U4cUM1juCqwL7N5RKQIUAY4bIw5bYw5BGCMWQfsBOpebNDZKX/qGAc16Sul/JyRIA4XK034Kf9K+muBOiJSU0SKArcBc73mmQsMdN73Br43xhgRiXAuBCMitYA6wK7cCT0TxlBeS/pKqQBxuHhpyufjXbnZ1uk7dfTDgW+BYOA9Y8wmERkPxBpj5gLvAh+KyA7gMPaHAaA9MF5EkoEUYKgx5nBe/CNpSp8+QUhqiiZ9pVRAsF0x+FGTTQBjzHxgvte4sR7vk4A+GSz3GfDZRcaYI2kXRLR6RykVCA4VL8vlB37Lt+0VuDty0y6IaElfKRUIDhUv7XcXcgNKBec0SZO+UioQHC5WhnJJx/Ot/50Cl/TTLoho9Y5SKhCk3ZWbX/3vFMCkb0+TjhTTpK+U8n9p/e+Uz6dmmwUu6Vc4eZRjRYtzpog+G1cp5f/SqqLDT+RPC54Cl/TLnzym9flKqYBxqHj+9r9TAJN+QvpOVEopf5eWr8qfyp8btApc0q9w8qiW9JVSASOt/x2t3rlA5U8d05K+UipgGAniSLFS+db/TsFK+sZQ7qQmfaVUYDlUvEy+9b9TsJJ+QgJFU5O1ekcpFVAOFy+Tb3flFqyk7zwbV0v6SqlAcqhYGSpo0r8A8fGAJn2lVGA5VEJL+hfGSfpavaOUCiRp/e+QnPf97xSspJ9WvaNdMCilAkh67cTBg3m+rYKV9LWkr5QKQOk5y8lheanAJX3td0cpFWgOadK/QAcOaClfKRVw0pO+U0WdlwpW0o+P53Dx0m5HoZRSOaLVOxcqPp5Dxcu6HYVSSuXI0bCSpEiQ/yR9EekqIttEZIeIjMpgeqiIzHSm/yAikR7TRjvjt4lIl9wLPQMHDmgbfaVUwEkNCuZIsVL+Ub0jIsHA60A3IAroKyJRXrPdDRwxxtQGXgGec5aNAm4DGgBdgTec9eU+Y+DgQa3eUUoFpMPFyvhNSb8FsMMYs8sYcwaYAfTymqcX8IHzfjZwjYiIM36GMea0MeY3YIezvtyXkABnz2obfaVUQDpUIn+SfhEf5qkC7PUYjgNaZjaPMSZZRBKAcGf8aq9lq3hvQESGAEOcwUQR2eZT9BlZ9G4FFr2b93c45FwFQOPyncaVMxpXzvhdXFcC7KECIhcaVw1fZvIl6UsG44yP8/iyLMaYKcAUH2LJlojEGmNicmNduUnjyhmNK2c0rpwpzHH5Ur0TB1TzGK4K7MtsHhEpApQBDvu4rFJKqXziS9JfC9QRkZoiUhR7YXau1zxzgYHO+97A98YY44y/zWndUxOoA6zJndCVUkrlVLbVO04d/XDgWyAYeM8Ys0lExgOxxpi5wLvAhyKyA1vCv81ZdpOIzAI2A8nA/caYlDz6X9LkSjVRHtC4ckbjyhmNK2cKbVxiC+RKKaUKg4J1R65SSqksadJXSqlCJOCTvoi8ICJbRWSDiHwhIhl2vpNdVxJ5EFcfEdkkIqkikmkTLBHZLSIbReQnEYn1o7jye3+VF5EFIrLd+Vsuk/lSnH31k4h4NyjIzXguuOuRvORDXINEJN5jHw3Oh5jeE5EDIvJLJtNFRF5zYt4gIk3zOiYf4+ooIgke+2psPsVVTUQWicgW57v4QAbz5N0+M8YE9AvoDBRx3j8HPJfBPMHATqAWUBT4GYjK47guB+oBi4GYLObbDVTIx/2VbVwu7a/ngVHO+1EZfY7OtMR82EfZ/v/AMOBN5/1twEw/iWsQMCm/jidnm+2BpsAvmUy/DvgGe99OK+AHP4mrIzAvP/eVs91KQFPnfSng1ww+xzzbZwFf0jfGfGeMSXuw5GrsvQDefOlKIrfj2mKMufA7i/OIj3Hl+/7i3K48PgBuyOPtZeViuh5xO658Z4xZim21l5lewDRjrQbKikglP4jLFcaY/caY9c7748AWzu+pIM/2WcAnfS93YX8dvWXUlcR53UG4xADficg6pzsKf+DG/rrEGLMf7JcCqJjJfGEiEisiq0Ukr34YfPn/z+l6BEjreiQv+fq53OxUCcwWkWoZTM9v/vz9ay0iP4vINyLSIL837lQLNie+xxQAAAJySURBVAF+8JqUZ/vMl24YXCci/wMuzWDSY8aYOc48j2HvBfg4o1VkMO6i26r6EpcP2hhj9olIRWCBiGx1SihuxpXv+ysHq6nu7K9awPcistEYs/NiY/NyMV2P5CVftvkV8Ikx5rSIDMWejVydx3Flx4195Yv1QA1jTKKIXAd8ib2BNF+ISEngM+BBY8wx78kZLJIr+ywgkr4xplNW00VkINADuMY4FWJe8qQ7iOzi8nEd+5y/B0TkC+wp/EUl/VyIK9/3l4j8JSKVjDH7ndPYDDsW99hfu0RkMbaUlNtJPyddj8R5dT2Sl7KNyxhzyGPwbZxuzl3ml92xeCZaY8x8EXlDRCoYY/K8IzYRCcEm/I+NMZ9nMEue7bOAr94Rka7Av4CexpiTmczmS1cS+U5ESohIqbT32IvSGbY0yGdu7C/PrjwGAuedkYhIOREJdd5XANpg7/bObRfT9UheyjYur3rfntj6YrfNBQY4LVJaAQlpVXluEpFL067DiEgLbD48lPVSubJdwfZisMUY83Ims+XdPsvvK9e5/cL20b8X+Ml5pbWoqAzM95jvOuxV8p3Yao68jutG7K/1aeAv4FvvuLCtMH52Xpv8JS6X9lc4sBDY7vwt74yPAd5x3l8JbHT210bg7jyM57z/HxiPLVwAhAGfOsffGqBWXu8jH+N61jmWfgYWAfXzIaZPgP3AWefYuhsYCgx1pgv2QUw7nc8t09Zs+RzXcI99tRq4Mp/iaoutqtngkbeuy699pt0wKKVUIRLw1TtKKaV8p0lfKaUKEU36SilViGjSV0qpQkSTvlJKFSKa9JVSqhDRpK+UUoXI/wM2hlRAezcw4gAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def wigner_density(x):\n", " return np.sqrt(4 - x*x) / (2*np.pi)\n", "\n", "N = int(1.e4)\n", "\n", "########################################\n", "# To Do: complete with an upper bound M\n", "# for the density f\n", "########################################\n", "upper_bound = 1 / np.pi\n", "\n", "X = np.array([])\n", "\n", "for _ in range(N):\n", " # There is no \"while ... do\" in Python\n", " # A possible solution is given below\n", " \n", " ########################################\n", " # To Do: sample the variables\n", " # z according to uniform distribution over [-2,2]\n", " # u according to uniform distribution over [0,1]\n", " ########################################\n", " z = np.random.uniform(high=2, low=-2, size=1)\n", " u = np.random.uniform(high=1, low=0, size=1)\n", " \n", " v = upper_bound * u # THIS IS M * U\n", " \n", " while v >= wigner_density(z): # WHILE THE CONDITION FAILS\n", " #################################\n", " # To Do: complete the code below\n", " #################################\n", " z = np.random.uniform(high=2, low=-2, size=1)\n", " u = np.random.uniform(high=1, low=0, size=1)\n", " \n", " v = u * upper_bound\n", "\n", " X = np.append(X, z)\n", "\n", "##############################\n", "# Plot the histogram\n", "##############################\n", "n_columns = 2 * int(N**(1/3))\n", "\n", "plt.hist( X , density=True, bins=n_columns, label=\"Empirical law\")\n", " \n", "########################################\n", "# Plot the true density for comparison\n", "########################################\n", "x = np.linspace(-2., 2., 100)\n", "\n", "f_x = wigner_density(x)\n", "\n", "plt.plot( x, f_x, \"r\", label=\"Density of the Wigner law\", linewidth=1.5)\n", "\n", "plt.legend(loc='best')\n", "plt.title(\"Wigner law, N = %1.0e samples with rejection method\" %N)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 2. Law of large numbers.\n", "\n", "Let $(X_i)_{i\\ge1}$ be a sequence of i.i.d. random variables such that $\\mathbb{E}\\bigl[ \\bigl|X_1\\bigr| \\bigr]<\\infty$.\n", "We remind that the law of large numbers states that the sequence of _empirical means_\n", "$$\n", "\\hat m_n = \\frac1n \\sum_{i=1}^n X_i\n", "$$\n", "converges almost surely to $\\mathbb{E}[X_1]$ as $n\\to \\infty$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 1.1 Convergence of the sequence $\\hat m_n$ as $n$ becomes large\n", "\n", "Draw $n$ independent random samples uniformly distributed $[0,1]$ and plot the sequence $n \\mapsto \\hat m_n$ for $n$ ranging from $1$ to $N$.\n", "\n", "The following functions from `numpy` can be used: `numpy.random.rand` to draw the $n$ samples from the uniform distribution on $[0,1]$, and the function `numpy.cumsum` to compute the cumulated sum up to $n$." ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt8VNW5//HPk0CAKlWupyh3BRVUQGMEr1S5ROoBUVFs64GjlmpFW2mrqLRysB4vtbU9LUotgtZfPYhKS7Ra73iKihC8tIAiFxGjKCFcBCGQkOf3x5okk8mQTEKuM9/36zWZvddee89aszPPXrP2nrXN3RERkdSQ1tgFEBGRhqOgLyKSQhT0RURSiIK+iEgKUdAXEUkhCvoiIilEQV9EJIUo6IuIpBAFfRGRFNKisQsQq2PHjt6zZ8/GLoaISLOyfPnyLe7eqbp8TS7o9+zZk9zc3MYuhohIs2JmHyeST907IiIpREFfRCSFKOiLiKSQJtenLyIHr6ioiLy8PAoLCxu7KFLHWrduTdeuXWnZsmWt1k8o6JtZNvBbIB2Y7e53xclzCTAdcOA9d/92JH0/8K9Ito3uPrpWJRWRhOXl5dG2bVt69uyJmTV2caSOuDsFBQXk5eXRq1evWm2j2qBvZunATGA4kAcsM7Mcd18VlacPcDNwurtvM7POUZvY4+4Da1U6EamVwsJCBfwkZGZ06NCB/Pz8Wm8jkT79LGCtu693933APGBMTJ7vATPdfRuAu2+udYlEpE4o4Ceng92viQT9I4FPoubzImnR+gJ9zex1M1sS6Q4q1drMciPpFxxUaavx1lvwzjv1+QoiIs1bIkE/3mEl9sa6LYA+wFDgMmC2mR0eWdbd3TOBbwO/MbOjKr2A2aTIgSH3YL62DB4MJ51U69VFpA6lp6czcODAssddd1U6FVgro0aNYvv27TVeVpWHH36YyZMnH2zRmoVETuTmAd2i5rsCn8XJs8Tdi4CPzGw14SCwzN0/A3D39Wa2CBgErIte2d0fBB4EyMzM1J3aRZJAmzZtePfdd+t8u88++2ylNHfH3eMuk4oSaekvA/qYWS8zywDGAzkxef4KfBPAzDoSunvWm1k7M2sVlX46sAoRSVk9e/bklltuYciQIWRmZvL2228zcuRIjjrqKGbNmgXAokWLOOussxg7diz9+vXj6quvpqSkpGz9LVu2sGHDBo477jh+8IMfcNJJJ/HJJ5+ULQP405/+xIknnsiAAQO4/PLLAXj66ac59dRTGTRoEMOGDeOLL76osqzTp09nwoQJjBgxgp49e7JgwQJuvPFGTjjhBLKzsykqKgJg+fLlnH322Zx88smMHDmSTZs2AfDHP/6RU045hQEDBnDRRRexe/duACZOnMj111/PaaedRu/evXnyySfr/o0+gGpb+u5ebGaTgecJl2zOcfeVZjYDyHX3nMiyEWa2CtgP/NTdC8zsNOAPZlZCOMDcFX3Vj4jUvx/9COq6wT1wIPzmN1Xn2bNnDwMHll+4d/PNN3PppZcC0K1bN958801uuOEGJk6cyOuvv05hYSH9+/fn6quvBmDp0qWsWrWKHj16kJ2dzYIFC7j44osrvMbq1auZO3cu999/f4X0lStXcscdd/D666/TsWNHtm7dCsAZZ5zBkiVLMDNmz57NPffcw69+9asq67Fu3TpeffVVVq1axZAhQ3jqqae45557GDt2LH/729/41re+xXXXXcfChQvp1KkTjz/+OLfeeitz5szhwgsv5Hvf+x4A06ZN46GHHuK6664DYNOmTSxevJgPPviA0aNHV6pbfUnoOn13fxZ4Nibt51HTDkyJPKLzvAGccPDFFJHmpqrundGjw891TjjhBHbt2kXbtm1p27YtrVu3LuuTz8rKonfv3gBcdtllLF68uFJg7NGjB4MHD660/VdeeYWLL76Yjh07AtC+fXsg/H7h0ksvZdOmTezbty+ha93PO+88WrZsyQknnMD+/fvJzs4uK/uGDRtYvXo1K1asYPjw4QDs37+fLl26ALBixQqmTZvG9u3b2bVrFyNHjizb7gUXXEBaWhr9+vWr9htHXdIvckWSXHUt8sbQqlUrANLS0sqmS+eLi4uBypcmxrtU8ZBDDom7fXePm/+6665jypQpjB49mkWLFjF9+vQalbVly5Zl2y0tq7vTv39/3nzzzUrrTpw4kb/+9a8MGDCAhx9+mEWLFlXabml5G4rG3hGRJmnp0qV89NFHlJSU8Pjjj3PGGWckvO65557L/PnzKSgoACjr3tmxYwdHHhmuOH/kkUfqpJzHHHMM+fn5ZUG/qKiIlStXArBz5066dOlCUVERf/7zn+vk9Q6Wgr6I1IvSPv3Sx9SpU2u0/pAhQ5g6dSrHH388vXr1YuzYsQmv279/f2699VbOPvtsBgwYwJQpoed5+vTpjBs3jjPPPLOs6+dgZWRk8OSTT3LTTTcxYMAABg4cyBtvvAHA7bffzqmnnsrw4cM59thj6+T1DpY15NeKRGRmZnptb6JS+m2uiVVJpMG9//77HHfccY1djFpbtGgR9957L88880xjF6VJird/zWx55DdRVVJLX0QkhehErog0OUOHDmXo0KGNXYykpJa+iEgKUdAXEUkhCvoiIilEQV9EJIXoRK6I1KmCggLOPfdcAD7//HPS09Pp1KkTEH5wlZGR0ZjFO2ivvPIKX/va1+IO/1BVvpkzZ3L44Yfzne98pyGKeUAK+iJSpzp06FA25s706dM59NBD+clPflIhT+lQyGlpza+z4ZVXXqFjx44JBf3ofNdee21DFK9aze8dF5Fmae3atRx//PFcffXVZUMhH3744WXL582bx1VXXQXAF198wYUXXkhmZiZZWVksWbKk0vaKi4uZMmUKWVlZnHjiicyePRuAJ554omxgs08//ZS+ffuyefNmZs+ezdixYxk5ciTHHHMMv/jFL8q29cgjj5CVlcXAgQP5wQ9+UDaM89/+9jdOOukkBgwYwIgRI1i3bh2zZ8/ml7/8ZdkvbxcuXFg2XPOIESPYvHlz3HzTpk3jN5GBkN5++21OPfVUTjzxRC666CJ27NgBhFFAp06dSlZWFsccc0zZL3vrklr6IsmuscZWjmPVqlXMnTuXWbNmlQ2sFs/111/PjTfeyODBg9mwYQPnn38+K1asqJDnwQcfpHPnzixdupS9e/cyePBgRowYwbhx43jqqaeYNWsWCxcu5I477qBz585A6F5asWIFGRkZnHLKKZx//vm0aNGCv/zlL7zxxhu0aNGCSZMmMW/ePM455xyuueYa/vGPf9CjRw+2bt1K+/btueqqq+jYsSM/+tGPANi2bRujR4/GzJg1axa/+tWvuPvuuyvli77By3e/+10efPBBzjjjDG655RZuv/127r33XiB8C1q6dCk5OTnMmDGDv//97zV+n6uioC8iDeaoo47ilFNOqTbfSy+9xOrVq8vmt23bxp49e2jTpk1Z2gsvvMD777/PvHnzgDCY2po1a+jevTszZ87k+OOP56yzzmLcuHFl64wcOZJ27doBYWjjxYsXU1xczLJly8jMDCMY7Nmzh27dutGmTRu++c1v0qNHD6B8eOZYGzdu5JJLLuHzzz9n79699O3bt8q6FRQUUFhYWDaA3IQJE8pu8gJw4YUXAnDyySezYcOGat+rmlLQF0l2TWhs5eihkNPS0ioMKVxYWFg2Xdrareqkr7tz//33l500jpaXl0d6ejqff/55hWGW4w3X7O5cccUV3H777RWWLViwIO7wzLGuvfZabrnlFkaNGsVLL71U7b2AqxvvrHTI5fT09Cq/DdWW+vRFpFGkpaXRrl071qxZQ0lJCX/5y1/Klg0bNoyZM2eWzce7GcvIkSO5//77ywLj6tWr2bNnD0VFRVxxxRXMnz+f3r1789vf/rZsnRdeeIHt27eze/duFi5cyOmnn86wYcOYP39+2W0WCwoK2LhxI6effjqvvPIKH3/8MVA+PHPbtm3ZuXNn2TZLh2t29wrDNcfmK9WxY0fatGlT1l//6KOPcvbZZ9f8DawltfRFpNHcfffdZGdn0717d/r168fevXuBcHnjNddcw9y5cykuLuab3/xmhYMAwPe//302btxYdkvGzp07s3DhQu68807OPfdcTjvtNPr3709WVhajRo0CwonSb3/726xbt47LL7+8bN3bbruNYcOGUVJSQsuWLZk1axannHIKDzzwAGPGjMHdOeKII3juuecYM2YM48aNY8GCBcycOZPp06czduxYunbtSlZWVtn9cWPzRXv00Ue55ppr2LNnD0cffTRz586t1/c5moZWFklCzX1o5fowe/ZsVqxYUXYFTXOmoZVFRCQhCQV9M8s2s9VmttbM4t7+xswuMbNVZrbSzB6LSp9gZmsijwl1VXARkZq46qqrkqKVf7Cq7dM3s3RgJjAcyAOWmVmOu6+KytMHuBk43d23mVnnSHp74DYgE3BgeWTdbXVfFRGJdqCbg0vzdrBd8om09LOAte6+3t33AfOAMTF5vgfMLA3m7r45kj4SeNHdt0aWvQhkH1SJRaRarVu3pqCg4KADhDQt7k5BQQGtW7eu9TYSuXrnSOCTqPk84NSYPH0BzOx1IB2Y7u5/P8C6R8a+gJlNAiYBdO/ePdGyi8gBdO3alby8PPLz8xu7KFLHWrduTdeuXWu9fiJBP973w9jmQwugDzAU6Ar8w8yOT3Bd3P1B4EEIV+8kUCYRqULLli3p1atXYxdDmqBEunfygG5R812Bz+LkWejuRe7+EbCacBBIZF0REWkgiQT9ZUAfM+tlZhnAeCAnJs9fgW8CmFlHQnfPeuB5YISZtTOzdsCISJqIiDSCart33L3YzCYTgnU6MMfdV5rZDCDX3XMoD+6rgP3AT929AMDMbiccOABmuPvW+qiIiIhUT7/IFRFJAvpFroiIVKKgLyKSQhT0RURSiIK+iEgKUdAXEUkhCvoiIilEQV9EJIUo6IuIpBAFfRGRFKKgLyKSQhT0RURSiIK+iEgKUdAXEUkhCvoiIilEQV9EJIUo6IuIpBAFfRGRFKKgLyKSQhT0RURSSEJB38yyzWy1ma01s6lxlk80s3wzezfyuCpq2f6o9Jy6LLyIiNRMi+oymFk6MBMYDuQBy8wsx91XxWR93N0nx9nEHncfePBFFRGRg5VISz8LWOvu6919HzAPGFO/xRIRkfqQSNA/Evgkaj4vkhbrIjP7p5k9aWbdotJbm1mumS0xswsOprAiInJwEgn6FifNY+afBnq6+4nAS8AjUcu6u3sm8G3gN2Z2VKUXMJsUOTDk5ufnJ1h0ERGpqUSCfh4Q3XLvCnwWncHdC9x9b2T2j8DJUcs+izyvBxYBg2JfwN0fdPdMd8/s1KlTjSogIiKJSyToLwP6mFkvM8sAxgMVrsIxsy5Rs6OB9yPp7cysVWS6I3A6EHsCWEREGki1V++4e7GZTQaeB9KBOe6+0sxmALnungNcb2ajgWJgKzAxsvpxwB/MrIRwgLkrzlU/IiLSQMw9tnu+cWVmZnpubm6t1rXI2YcmViURkXpnZssj50+rpF/kioikEAV9EZEUoqAvIpJCFPRFRFKIgr6ISApR0BcRSSEK+iIiKURBX0QkhSjoi4ikEAV9EZEUoqAvIpJCFPRFRFKIgr6ISApR0BcRSSEK+iIiKURBX0QkhSjox/j4Y7jtNti/v7FLIiJS9xT0Y0yZAjNmwD/+0dglERGpewr6MUpvtfj++5WXPf88PPVUw5ZHRKQuVXtj9FSze3d43rq18rLs7PCse/CKSHOVUEvfzLLNbLWZrTWzqXGWTzSzfDN7N/K4KmrZBDNbE3lMqMvC17WiInjttTAdG/TXrm348ohI01dQAG+/DSUlVefbswc2bYKpU+Gaa+DRR2HRIvjRj+DUU+Eb34B///f6L2+1LX0zSwdmAsOBPGCZmeW4+6qYrI+7++SYddsDtwGZgAPLI+tuq5PS17GCAigsDNOxQb9Pn/LpffsgIyP+NnJzYdAgSE+vnzKKNGcrVsCHH4bPz3nnVfycuMOXX0KbNtCyZUgzC8/FxaFRtn07dOkCL70Eq1bB0UfDkUfCgAH1U95t22DlSvj0U/jiC1i+PJSxQ4fQECxtJJaaMAF++lM47jhIS4ONG2Hu3LDeq6/Crl3leWfNCs/p6SFmZGdD3771U49oiXTvZAFr3X09gJnNA8YAsUE/npHAi+6+NbLui0A28L+1K2792rKl8vSWLXDRRRXzXXZZ/L79m2+Gu+6C66+H3/62/sop0lw89hj88Y9w2mlw1FFw5ZXly848EyZNCoHzzTdDcI3Wti1ceGEIlhs3Vv9ao0fDf/5nOJgUF8PXvlZ+0Ii1dy88+STk50O/fqGVvXIlPP54CPTf/ja8/DLk5IS8B5KRAQMHwtixcO+9ofX+yCPQokVYVtpdDHDSSTBkCAwfDoccEg5yBQVwxhnQtWv19asr5tV0UJvZxUC2u18Vmb8cODW6VW9mE4E7gXzgQ+AGd//EzH4CtHb3X0Ty/QzY4+73Huj1Mtu29dyTT65VZRZFjrpDz67V6mzfDu++F1WWzNByL/X1tvDlzjDdtw+0aw9tWld+/Tatw9c1kVTgwNYC+Pzz8Blq3RrS0kN3Rum34n37Qt5DvgaHHQZFxWGd/ZEukZYtQlrnzrBvbwjYO3fBzp2h1d++Xci3vwRaZYClhddxD9/OP/20crnatoUjuoQ8RcWwYQN06hQ+x59+CnsKK69jkfpAaKl3aA//9m8hiKelhW1iId++fdAyI0yX2rsPPt8UnktKQtk7tA/fXjJaVcxb1+y115a7e2Z1+RJp6ccrZ+yR4mngf919r5ldDTwCnJPgupjZJGASwImtWiVQpPpRVFRx/v2Y7zJHdoUvI1f1fLgm/DOfNiTMR1dq777wj3agVoZIMti9Gz76CPK3VEwv2gVpFgJip05w7LGwbSt8tRu6HlnepVO4Fwq2hEBf2p0TzSN/Evkc9Tk6tMg3b4YvNodulF27YPWHFfNt3hweh3wttPD37w/B24C2Xw8HpL17Ycf20KhrdYBuXIjfxdsqA3r0qL68jcrdq3wAQ4Dno+ZvBm6uIn86sCMyfRnwh6hlfwAuq+r1Tj75ZK+tEGprvbo/8EBY/89/Ds+tWpVvE9xff919+vSKaaXuuivMDx4cnl96qXzZ9u3uJSW1L5c0P/v2lU8/+qj7Qw+F/4NoO3e633ef+wcfhP+tRYuaz//JP/7h3rp1+efg0kvdX3jBvbjYfdWqxi1bSUl4bN/uPnOm++9+5/7ii+7797tv2eL+xhthOtkAuV5NPPdI2Kou6LcA1gO9gAzgPaB/TJ4uUdNjgSWR6fbAR0C7yOMjoH1Vr9dQQX/cOPcbbqiYdvvtYf3CQvcjjqgY3ME9Pz98eKPTNm6s+Np//rN7y5buAwa4T5nifvnl5cs++6zWVWsSmktAagwffOC+cKH7pEnuRx8d9veFF7qfdFLF/5dTT3V/6in33//evV27yv9j3bu79+oV/o9q6kD7Z/1693vvdf/DH9y3bStP37nT/eab3Q8/3P2009z//nf3goKwbPNm97/+1f35591zctzPO8/9ggvcTz7ZfeRI97Q0986d3Z94oubllPpRZ0E/bItRhL76dcCtkbQZwOjI9J3AysgB4VXg2Kh1rwDWRh7/Wd1rNVTQj5f3hz90b9s2TJe22IcPDx/QHj3Ch+rxxyt/UD/+uHz6vfdCwI/NA+6zZ4eW0Lhx7hdfHFodzcGMGSFYgfv777vv2eO+dWvi6//rX+633BLq++WXB1eWkpLGa6Xt3RuCe6zFi91btCjfz7HB/Kqr3K+7LgTJ6PTMTPc77nA//XT3s84KaR07ln/D7N49LJs8OQTe0qC+c6f7//yP+9Ch7llZ7t/7nnvPnu5f/7p7797uP/6x+9y57rfdFoK0WcWDyoIFIcB36xb//zT2G270o7R1f/nl7jt2NOS7L9Wp06DfkI/GCvrbt4f5ww4L8xdeGOanTau43pYt7iee6H799eXb+NOfwvOVV4YP5oIF8T8wvXq5v/lm+fwll9S6qg3m//2/yvVo0cL9G99I7EMffUCMDhz//Gf8/Lt2hQNMmzbut94a9seQIe5ffeX+9NNh/fHjK5ZtzpzK2/ngA/eiojD905+G4LhrV+V8X37p/s477h9+GPK/9577qFEhSF9xhft//Edo2X7jG+Xlv/RS9y5dwjfF7OzyfTtxYjiwu4funblz3TdsqPh677wTGgXXXFOxCyjanj3hNWPft0GDwnsR73+rc2f3ESPiLxs/3n35cvfnngv/u6XpnTq533hjOCivWRO+pUR/I7nttlCnm25yX7nS/e23wwFX3/iaJgX9Gub95S8rpn3nO2H697+Pv/7One4ZGSHPD34QnnfvLl8+dGj59r7/fff58yt/GPv1q109G8qaNe7p6aEuDz0UWpCxdViypPJ6u3e7L1sWguaBWozgPmZM1curehx7bMX5hx4Kr71qlXv79uUH1XvuqZjvyivd164NQfDnPw/f7GpbhtJHRkY4uNWl/fvDQaGgIHT13H9/aKWX/s9Nmxa6FhcsCAG5NBDv3h3qPH16SI89MO/ZE/5fb731wAft3bvDN1JpXhT0D2DrVvc776yc9777KqZdeWWYvv/+A29r2bKKH/5oX30VWpjgfvfd4UOUlhbmb7ghfO0+5JDy1mhTUVTkvm5dmB49OrTKN20K8/v3h+C6fn0of3Tw7to1dAOtXVu5G2PUqBCUXnwxtCqjW83xHldc4f7RR+Eg89pr4dxIhw6h5fn00+7nnhvyDR4cTh6ec07V2+vfP3R/xFvWr5/7hAkhD7gPHOj+zDPu114b+rvvuy8E15KSEDDdQ1B86qlwsv53v2v4czWl5RCJpqB/ANF9r1DeQpo5s+L6v/51mK7qRFVh4YGDvnsI9HPmhL5g93DOAEJ/6mOPhel33ql5PatSUlLxZN2BrFoVrsDYurW8j/yZZyoHxVtuib/+7t0hmFcVbIcPD91hhYUV1y0sDN0bs2eHg8CHH4aWclFROKBU131QWBi6yUrzrVtX/potWrjPmxcC4+DB4fzJ3r3hPXntNffc3NB1N3Cg+89+Fg7OIslAQb+aPKWP0g/9vfeG+QsuCPPFxRVPniWyveoUF7u/+mrY5vr1Xu03idoYMyacuHv77XBAiQ64S5aElunatRXLPWqU+3//d/zAXVUrNj8/9FtDuLqjdJ0//SkE2thgX58+/TS8v9HfnNT3LKlEQT+OkpIDB7Wf/SzM17Qv87nnwnrDhtVsvZKS0CVy/PE1W68q0S3e0kfpSc/bbitPmzIlfoCH0G+8Zk3oukjkG4N7aM3v3Rv6kH/zGwVbkcaQaNBPqfH0442hMWdOeN6+Pfwar6YDpWVnh/Eznn66ZuuZwRVXhAGoFiyo/GvgmnrttTC2CcDvfleePm8eXHop/Nd/laf9+tdhLJSf/xxOOaU8/Re/CHcOO/poOPdcOPzwxF67Q4fw68R+/eCHP9QvkUWasqQM+u7x06MHPyo1bVp4/uyz8JPx2mjfPowDUlOlQwxddFEYxK0mCgvh8svh+98PI/oNHRrSr7oKJk8ON3y5556QNn9+5fXHjQsHgqVLw0/Rn3wSfvzjmtdBRJqXagdca2iZmZmeGz3KWQ2UtjBLSuK3Nv/2Nzj//Mrpn3wC3bqF6YZ8O3btigzgFLFjB3z961Wv88tfwo03Hnj5V1+FwapKXXZZaO0ffXQ4ELzzThi/+6234IgjDq78ItJ0mFlCA66lVEs/OuAfd1z5IE+NdYOUQw8N3Ur/939h/oknqs6/d2/FgH/DDfD738Ptt4f5xYsrBnwIw7yOHRvy9e4dvlV88okCvkiqSsrbJSbSWk9LC+Nf//CHYdQ9CNMN7bDD4PTTQ//50qUVxxuP9fDD5dNTp8Kdd4Zp93B+IF4gz8gI5wxERCCFgn5sWnTXSmlLf8KE+i3XgaSlhROqS5dWne+558LNFl57DXr1Kk83U8tdRBKTlN078Xz1VcX50nG0AW69NTxX159en7Ky4F//CjeeiGfbtnAXn3HjQjeNrpARkdpImZb+9u3h+YEHQsv+u98Nd8OJFn1StaGdcko4EL3zTricsrg4lK+kBMaPD5eEuodbwYmI1FbKBP0dO8Jzu3ahL79UixYhwEK49LKxZGWF5//7v3C/0J/8BO6+G957r+IJ3pNOapzyiUhySLmgf9hh8fM+/njlln9D6tIFTjgh3Fy91E03VcxzzDHhh1AiIrWVlH36NQn6JZEbMyf669P6NHFi/PR//3fYsgVef71BiyMiSSjlg35p3qYQ9MePD8+jR8PKleXpTzwRWvhq5YvIwUqZ7p0VK8JzbHDv1w9WrYJDDqn/clXniCNg587wA6u0tDD+TXo6tGrV2CUTkWSRMkH/rrvCc2zQf/758GvVY4+t/3Il4tBDy6dfeEGXZopI3UrKoF/aTx89v38/nHlm5WEKunYtPyA0NWlJ2fkmIo0pobBiZtlmttrM1prZ1CryXWxmbmaZkfmeZrbHzN6NPGbVVcGrEtvSz88Pz5dc0hCvLiLSdFXb0jezdGAmMBzIA5aZWY67r4rJ1xa4HngrZhPr3H1gHZU3IbFBP/oafRGRVJZISz8LWOvu6919HzAPGBMn3+3APUBhHZavVmK7dw505Y6ISKpJJOgfCXwSNZ8XSStjZoOAbu7+TJz1e5nZO2b2mpmdGe8FzGySmeWaWW5+aV/MQThQS78xx9YREWkKEgn68a4fKQurZpYG3AfEu+/SJqC7uw8CpgCPmVml0OvuD7p7prtndqrt7asqbK/ivFr6IiJBIkE/D+gWNd8V+Cxqvi1wPLDIzDYAg4EcM8t0973uXgDg7suBdUDfuih4VWK7d778MjyrpS8iqS6RoL8M6GNmvcwsAxgP5JQudPcd7t7R3Xu6e09gCTDa3XPNrFPkRDBm1hvoA6yv81rEiG3pl94QvTb3sRURSSbVXr3j7sVmNhl4HkgH5rj7SjObAeS6e04Vq58FzDCzYmA/cLW7b62Lgldd5orzpePm65etIpLqEvpxlrs/CzzrdN+RAAALv0lEQVQbk/bzA+QdGjX9FPDUQZSvVmK7d0pb+hkZDV0SEZGmJSl/86mWvohIfCkR9Etb+o05Xr6ISFOQlEE/tntn377QtaPBy0Qk1SVl0I/X0lfXjohIigT9+fPDOPUiIqkuKYN+dPeOO3z22YHzioikkqQM+tEt/eLixiuHiEhTk/RBv6io8cohItLUJGXQj+7eUdAXESmXlEE/Xkt/8uTGKYuISFOSMkG/f//GKYuISFOSlEE/XvdOy5aNUxYRkaYkKYN+dEu/dNwdDbYmIpICQV8tfRGRcgr6IiIpJCmDfmmffn4+bNkSphX0RUQSvIlKc1Pa0u/cuTxNQV9EJIla+tFdOrEDroGCvogIJFHQjxY7nj7o6h0REUgw6JtZtpmtNrO1Zja1inwXm5mbWWZU2s2R9Vab2ci6KHQ8aumLiFSv2j59M0sHZgLDgTxgmZnluPuqmHxtgeuBt6LS+gHjgf7AEcBLZtbX3ffXXRUqixf0W7euz1cUEWkeEmnpZwFr3X29u+8D5gFj4uS7HbgHKIxKGwPMc/e97v4RsDayvXr1wQeweHHFtPbt6/tVRUSavkSC/pHAJ1HzeZG0MmY2COjm7s/UdN26Et26v/xyOPPMissV9EVEEgv68W4nXhZizSwNuA/4cU3XjdrGJDPLNbPc/Pz8BIpUMy1bwqGH1vlmRUSanUSCfh7QLWq+KxB9A8K2wPHAIjPbAAwGciInc6tbFwB3f9DdM909s1OnTjWrQdk2DrysfXuweIcfEZEUk0jQXwb0MbNeZpZBODGbU7rQ3Xe4e0d37+nuPYElwGh3z43kG29mrcysF9AHWFrntahGhw4N/YoiIk1TtVfvuHuxmU0GngfSgTnuvtLMZgC57p5TxborzWw+sAooBq6t7yt34lHQFxEJEhqGwd2fBZ6NSfv5AfIOjZm/A7ijluVLWHXdOyIikqS/yI2loC8iEiRN0K+qpa8rd0REgqQJ+lXRr3FFRIKkCfpVtfQ12JqISJA0Qb8q8UbdFBFJRSkR9Pc3+EWiIiJNU9IE/aq6d4qLG64cIiJNWdIE/aqoe0dEJEiaoF9VS1/dOyIiQdIE/aq0SMrbv4uI1FxSh8OhQ2HgQPh53AEjRERST9IE/XjdO+npcN99DV8WEZGmKqm7d776qrFLICLStCRN0I/X0t+1q+HLISLSlCVN0I9n587GLoGISNOS1EG/b9/GLoGISNOStCdyX34ZBg1qnLKIiDRVSRP0Y51zTmOXQESk6Uma7p2qfpErIiJB0gR9ERGpXkJB38yyzWy1ma01s6lxll9tZv8ys3fNbLGZ9Yuk9zSzPZH0d81sVl1XoJRa+iIi1au2T9/M0oGZwHAgD1hmZjnuvioq22PuPiuSfzTwayA7smyduw+s22KLiEhtJNLSzwLWuvt6d98HzAPGRGdw9y+jZg8B1O4WEWmCEgn6RwKfRM3nRdIqMLNrzWwdcA9wfdSiXmb2jpm9ZmZnHlRpq6DuHRGR6iUS9C1OWqUQ6+4z3f0o4CZgWiR5E9Dd3QcBU4DHzOzrlV7AbJKZ5ZpZbn5+fuKlFxGRGkkk6OcB3aLmuwKfVZF/HnABgLvvdfeCyPRyYB1Q6Xey7v6gu2e6e2anTp0SLXvMNmq1mohISkkk6C8D+phZLzPLAMYDOdEZzKxP1Oy3gDWR9E6RE8GYWW+gD7C+LgouIiI1V+3VO+5ebGaTgeeBdGCOu680sxlArrvnAJPNbBhQBGwDJkRWPwuYYWbFwH7ganffWh8VERGR6pk3sX6RzMxMz83NrfF627ZB+/bl802sWiIi9crMlrt7ZnX59ItcEZEUkjRBXy17EZHqJU3QFxGR6iVN0FdLX0SkekkT9EVEpHoK+iIiKSRpgr66d0REqpc0QV9ERKqXNEFfLX0RkeolTdAXEZHqKeiLiKSQpAn67dvDZ1UN+CwiIskT9NPToUuXxi6FiEjTljRBX0REqqegLyKSQhT0RURSiIK+iEgKUdAXEUkhCvoiIilEQV9EJIUkFPTNLNvMVpvZWjObGmf51Wb2LzN718wWm1m/qGU3R9ZbbWYj67LwIiJSMy2qy2Bm6cBMYDiQBywzsxx3XxWV7TF3nxXJPxr4NZAdCf7jgf7AEcBLZtbX3ffXcT3KPPww9OpVX1sXEWneqg36QBaw1t3XA5jZPGAMUBb03f3LqPyHAKVjXo4B5rn7XuAjM1sb2d6bdVD2uCZMqK8ti4g0f4kE/SOBT6Lm84BTYzOZ2bXAFCADOCdq3SUx6x4ZZ91JwCSA7t27J1JuERGphUT69C1OWqXR6919prsfBdwETKvhug+6e6a7Z3bq1CmBIomISG0kEvTzgG5R812BqsaznAdcUMt1RUSkHiUS9JcBfcysl5llEE7M5kRnMLM+UbPfAtZEpnOA8WbWysx6AX2ApQdfbBERqY1q+/TdvdjMJgPPA+nAHHdfaWYzgFx3zwEmm9kwoAjYBkyIrLvSzOYTTvoWA9fW55U7IiJSNfMmdnPZzMxMz83NbexiiIg0K2a23N0zq8unX+SKiKQQBX0RkRTS5Lp3zCwf+LiWq3cEttRhcZoD1Tk1qM6p4WDq3MPdq73mvckF/YNhZrmJ9GklE9U5NajOqaEh6qzuHRGRFKKgLyKSQpIt6D/Y2AVoBKpzalCdU0O91zmp+vRFRKRqydbSFxGRKiRN0K/u7l7NlZl1M7NXzex9M1tpZj+MpLc3sxfNbE3kuV0k3czsfyLvwz/N7KTGrUHtmFm6mb1jZs9E5nuZ2VuR+j4eGQeKyLhOj0fq+5aZ9WzMcteWmR1uZk+a2QeRfT0kBfbxDZH/6RVm9r9m1joZ97OZzTGzzWa2IiqtxvvWzCZE8q8xs1rfOSQpgn7U3b3OA/oBl0XfsrGZKwZ+7O7HAYOBayN1mwq87O59gJcj8xDegz6RxyTggYYvcp34IfB+1PzdwH2R+m4DroykXwlsc/ejgfsi+Zqj3wJ/d/djgQGEuiftPjazI4HrgUx3P54wrtd4knM/Pwxkx6TVaN+aWXvgNsK9TLKA20oPFDXm7s3+AQwBno+avxm4ubHLVU91XUi4deVqoEskrQuwOjL9B+CyqPxl+ZrLgzAE98uEm/E8Q7gvwxagRez+JgwEOCQy3SKSzxq7DjWs79eBj2LLneT7uPTmTO0j++0ZYGSy7megJ7CitvsWuAz4Q1R6hXw1eSRFS5/4d/eqdIeu5i7ylXYQ8Bbwb+6+CSDy3DmSLRnei98ANwIlkfkOwHZ3L47MR9eprL6R5Tsi+ZuT3kA+MDfSpTXbzA4hifexu38K3AtsBDYR9ttykns/R6vpvq2zfZ4sQT+hO3Q1Z2Z2KPAU8COveE/iSlnjpDWb98LMzgc2u/vy6OQ4WT2BZc1FC+Ak4AF3HwR8RfnX/XiafZ0jXRNjgF7AEYR7a58XJ2sy7edEHKiedVb/ZAn6SX2HLjNrSQj4f3b3BZHkL8ysS2R5F2BzJL25vxenA6PNbAPhLmznEFr+h5tZ6f0foutUVt/I8sOArQ1Z4DqQB+S5+1uR+ScJB4Fk3ccAw4CP3D3f3YuABcBpJPd+jlbTfVtn+zxZgn61d/dqrszMgIeA993911GLcojcrCbyvDAq/T8iVwEMBnaUfo1sDtz9Znfv6u49CfvxFXf/DvAqcHEkW2x9S9+HiyP5m1UL0N0/Bz4xs2MiSecSbjyUlPs4YiMw2My+FvkfL61z0u7nGDXdt88DI8ysXeRb0ohIWs019gmOOjxRMgr4EFgH3NrY5anDep1B+Br3T+DdyGMUoT/zZcKtKV8G2kfyG+FKpnXAvwhXRzR6PWpZ96HAM5Hp3oRbba4FngBaRdJbR+bXRpb3buxy17KuA4HcyH7+K9Au2fcx8F/AB8AK4FGgVTLuZ+B/Cectiggt9itrs2+BKyL1Xwv8Z23Lo1/kioikkGTp3hERkQQo6IuIpBAFfRGRFKKgLyKSQhT0RURSiIK+iEgKUdAXEUkhCvoiIink/wPy4wPSIes1PgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "N = 1000 #sample size\n", "\n", "########\n", "# Complete with N draws from the uniform distribution on [0,1]\n", "#\n", "X = np.random.uniform(0, 1, N)\n", "#\n", "########\n", "\n", "integers1toN = np.arange(1,N+1) \n", "# A numpy array containing the integers from 1 to N\n", "\n", "############\n", "# Stock in the variable 'empMean' the sequence of empirical means for n ranging\n", "# from 1 to N\n", "#\n", "empMean = np.cumsum(X) / integers1toN ## WANT: AN ARRAY OF SIZE N\n", "############\n", "\n", "############\n", "# Plotting\n", "############\n", "plt.plot(integers1toN, empMean, color=\"b\", label=\"Empirical mean\")\n", "\n", "# A horizontal line for the reference theoretical value\n", "plt.axhline(0.5, color=\"r\", label=\"True expectation\")\n", "\n", "plt.legend(loc=\"best\")\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 1.2 Another example\n", "\n", "The Cauchy distribution with parameter $a > 0$ is the probability distribution on $\\mathbb R$ with density function\n", "$$\n", "f(x) = \\frac a{\\pi(a^2+x^2)}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $\\blacktriangleright$ Question (a):\n", "\n", "$F(x)=\\int_{-\\infty}^x f(y) dy $ is the cumulative distribution function (cdf) of the Cauchy distribution.\n", "Let $U$ be a r.v. uniformly distributed on $[0,1]$. \n", "\n", "What is the law of $F^{-1}(U)$?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $\\blacktriangleright$ Question (b):\n", "\n", "Noticing that $F(x) = \\frac1{\\pi}\\arctan(\\frac xa)+\\frac12$, generate $N$ i.i.d. samples distributed according to the Cauchy distribution wit parameter $a=1$, still using the function `numpy.random.rand`.\n", "\n", "Plot the sequence of the empirical means $\\hat m_n$ for $n$ ranging from $1$ to $N$.\n", "\n", "What do you observe?" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsXWd4VEXbfs629N5DEjoJIRBKAOlNpEkRRLHQFZCXIiKi2BBFRcDewNeCLwqCKE0UEQRE6U16Cx0CqRAgdTPfj5v55uxmy9nNbhLCua/rXGfLKXPmzNzzzNNGYoyRChUqVKio+tBUdAFUqFChQkX5QCV8FSpUqLhLoBK+ChUqVNwlUAlfhQoVKu4SqISvQoUKFXcJVMJXoUKFirsEZSZ8SZJiJUn6U5KkI5IkHZIkaeLt34MlSVonSdKJ2/ugshdXhQoVKlQ4C6msfviSJEURURRjbI8kSX5EtJuI+hHRMCLKYoy9LUnS80QUxBibWtYCq1ChQoUK51BmCZ8xdpkxtuf251wiOkJE1YioLxEtuH3YAsIgoEKFChUqKghllvBNLiZJNYhoMxElEdE5xlig7L9sxlgptY4kSaOIaBQRkY+PT7OEhASH73vqFFFOjviu0xElJ5sec+4cUXo6UVISUVoaUUYGUf36RN7epscdOUKUl0cUGkoUF4fvGg1RfDzRgQNEhYU4LiSEqEYNh4tabjhxgqi4GM94J6GggOjgQdRtSEhFl0bg2jWikyeJEhKIfHwsH3PlCtGFC0RNmqDNVEYwhnJmZxPdukVkMBBFRaGuJamiS6fCWezevTuDMRZm90DGmEs2IvIlqHP63/6eY/Z/tr1rNGvWjDmDfv0YQ1PGFhJS+pjOnfHf2rWMDRyIz999V/q4pk0Z0+sZGzsW38eMYczfnzGjkbFBg0zvU5nRsydjKSkVXQrHceIE6vZ//6vokphi9WqUa8cO68fMno1jcnPLr1zOoqSEsd9+Y6xFC5S5Vi3GvvmGsaKiii6ZCmdARLuYAp52iRwiSZKeiJYR0XeMsZ9u/3zltn6f6/mvuuJeSsClcDnCw7G/elVI5qdP279WUhLR9euQitq2Fb/rdJCQVKi4EyFJRN26EW3bRrR6NVFgINGwYUSJiUQbNlR06VS4C67w0pGI6EsiOsIYe1f210oiGnr781AiWlHWe1mD+fTZEuFHRGB/5QpUNUTKCL9ePeyPHCFq0wafExOhLtmyxbnyqlBRWSBJRL16Ee3aRbR8OVF+PlHXrkRPP010+TLR998T7dhR0aVU4Sq4QsJvQ0SDiaizJEn7bm89iehtIuoqSdIJIup6+3u5wB7hGwz4rITwmzXDfscOothY8bteT7R+fdnKqUJFZYEkEfXtC4l/6FCijz8mio4meuwxopYtMQhs2gRlpoo7F67w0tnCGJMYY40YY41vb2sYY5mMsS6Msbq391muKLAlmBubLDVKf3/sr1wRv505Y//awcEw2G7ZIq57+jRR69ZEa9c6VVwVKiotoqOJvvqK6NgxokcewW8dO8JhoWNHCEtvvEHUrx/Rhx8S5eZWZGlVOIpK6kvgevBBQU74SiR8IjT0zZuhxiGCF0/jxkT79xOdPevSYqpQUSlQuzbUOUVF0OmfPg2CNxqJXn6ZaMUKookTMet9+WVV8r9TUCUI3xEXODnhMwY3QHvo2BGSzL59pe+5apXye6tQcadBp4Ow5OVFNH48UUkJ0fbtRM8/T/T770RNm0LiT0+v6JKqUIIqQfiO+A/LCZ9ImYTeujX227djHx0N6T4hgWjlSuX3VqGiKqBFC6K33oJef8gQ/HbzZsWWSYUy3PWEr0StExtL5OdH9M47+F6/PtFff8GtbeNGBOWoUHE3ggcuqrr8OwN3HeGXlJh+V2K4lSQYbrnffcOG0G2GhGD/22/K769CRVVCVBT2nTsTTZuGSGMVlRdVgvDLAqWG27feEp9r1QLZHztGFBamqnVU3L1o1w6z3PbtiWbNQlDjoEFEW7eqhtzKiCpB+I7mAOENMTQUunglkEfZarVEPXsS/for9mvWQNJXoeJuRIcORD/9hFxDTz+NGW/r1kTNmxN9+aWq369MqBKE7ywaN4Z/sRJ4euJ4jt69ibKyEImbkwOdvgoVdzNq1iSaMwdqnU8+QdTuE08QVasGD5+DByu6hCqqBOFrtfaPsTS9TEwkunhRudG1a1fsL16EwVavR24eLy9IOCpUqCDy9SUaOxbC1F9/Ed1/P9H8+bB9tWlDtGCBmoeqolAlCN8SlOgPExOxP3JE2TXvuw/7U6cQuduhA9Q6vXoR/fgjglJUqFABSBJUoQsXQkiaO5coMxNJ2qpVI5owQZX6yxtVgvAt6fAt5dMxR4MG2B8+rOw+nTsTvf02Gi4R0YABRMePI9/OlStEu3cru44KFXcbQkOJnnkGwtXGjbB9zZsHqb91a6Kvv0aOfhXuRZUgfEtQMmWsUQO6+UOHxG+2ZgYaDdHUqZBOiIgGDoRaJzUV31U9vgoVtiFJmBl/952Q+rOyiEaMQILDfv2gHlUSAa/CcVQJwrdE0nl59s/TahEtKyd8RxASQtSjB9Evv0DKX7zYueuoUHE3Qi71b95MNG4cotkHDEAK84yMii5h1UOVIHxLKh0lhE+EBU4sEb5SV8/HHiO6dIkoJQU5xY8eVXaeChUqAEmCP/+77xKdPw+16dWrmAGocC2qBOFbglIvgAYN4EbmbHqE3r2RduHaNcwYvv7aueuoqNxQg4jKBzqdsK2psS2uR5UlfKUSvqOGW3N4eWEKumYN1DsLFqgNtSpDXejb/dDpsFf7ketRJQjfUnpkRwnfWT0+EdQ616/DzfPKFTW3jgoVZYFejz1ff0KF61AlCN8SlKp0atRAxr+y+AN36oRF0s+ehafBV185fy1rmDULC06oqgUVVR2c8FUJ3/VwCeFLkvSVJElXJUk6KPttuiRJF83WuXULHJXw5aSp0SDdcVkkfK0WQVlr1xI9+ijR6tWl0zCXFc8/jxWHvvnGddc8epTojz9cdz0VKlwBTvhKZ+kqlMNVEv43RNTdwu/vyde5ddG9FEGJhM/1sQ0alI3wiYhGjkROnbAwTEUXLizb9czRrh3248Ypjwy2h6efRrqIH35wzfVUqHAFAgOxf+ABonvvJfr0U3VFLVfBJYTPGNtMRG5bpNwZOJKhLymJ6PLlskX6deiAnPkrVxK1aoUsga5Uv2i1RHXqEPn4ED38sH3pR8m9PT2xHzIEwS8qVFQG1K9PtG0b0aRJ8KD7z3+Qd79nT/yuwnm4W4c/TpKkf2+rfILcfC8T2CJ8czJ0heFWkojGjEGD7N4dUvi6dc5fzxKio4m+/RZJqZ55puzX8/LCvrCQ6H//K/v1VKhwFVq2hN3q6FGif/8lmjIFkezt22MN3VOnKrqEdybcSfifEVFtImpMRJeJaK6lgyRJGiVJ0i5JknalOzlvs+QqJyf8c+cwLbQGVxA+ESRlT08Yb6OjIaG4Gt27o/F//jkSttmCEhfCOnWQwXDmTOQFUqGisqFhQyxAtGMHBoKXX0a75YPCsWMVXcI7B24jfMbYFcaYkTFWQkRfEFELK8fNZ4ylMMZSwsLCnLqXPcJ/801TP3tzCT8uDildy5q5LzgYLpqLFiEnyOHDrjWycrzxBhaSfuKJsks6Gg0GDyKopVasgIrr99/LXk4VKlwJvpb0uXMg+pISODMkJOC/F15AagbzZUxVCLiN8CVJipJ9fYCIyjURqpzwY2JsHytJWNxk586y37drV+jXR47E97kW5zVlg8EAQ6tGg6CvsnozJCUhgyFjGKgOHRILtKtQUdkQG0v03HPor+fOEX38Mfr4nDlE99yDz2PGIB5GTcJmCle5ZS4ioq1EFC9J0gVJkkYS0TuSJB2QJOlfIupERG5QcPD7l/5NTvixsfav0aYN0Z49ZZcOuCFUkpAT5OBB9+gba9SA3n3/fnjulBWJiUSbNkEVxdGpE9Hs2arvv4rKi9hYGHXXrUP+nYUL0ZcXLkTke1gY1thdtMj59ClVCa7y0nmEMRbFGNMzxmIYY18yxgYzxhoyxhoxxvowxi674l6WYInwb9wQn3187F+jTRsEepR1JR5O+Pn50Ol7eIA0ywpLpNurF9GLLyLQ68svy36P+HhkLXzuOXSemBh8njJFnSarqPwICoJKdelSZNr85ReQ/Z9/Ij4mLAwz108+gb3qbhRkqmykrSNeOkRYhIHIdKBwBpzwb92CK9mwYUiodtkFw52lge211+CrPGYM0b59Zb9H7drQj4aFwfg8bhzUUg89pKarVXHnwNMTbpzz5yOb7d9/I+7k9Gm06fh4zA6GDIGd7ezZii5x+aBKEL49lY4ShITA+OPoeeaoWRN77i88eTLcHpcsKdt1rUGrhT4/KIjokUcgmbsKGg2ie2fPRnxB06Zl92RSoaK8odVCoHvnHXj0HD+O1bbatoWef/hwqEhr1yZ68kmof9LSKrrU7kGVIHxLcIa4W7fGeWWZ6sXFYdr48su4Vt26RLVqQcp315q3wcFw0Tx7FikeXBlEJUlEzz5L9M8/iCBu04ZowwbXXV+FivKEJKFPjhqFBYuuXEFcywcfwP1z6VKof6Ki4K49fjzRzz9XncDEKkv4ttQP1gi9TRvXkPJLL+Ee77yD72+/DePqvHllv7Y1tG8Pl8qDBzFbyc21f44jA1tKCmYtMTGIBVi0yPmyqlBRWSBJ8FKbMIFo+XIssr5zJ9SasbGwj/Xvj9W5mjaF8LNmjbL+VRlRJQjfkkonJ8f68ZzozM9r0wb7kpKy5T1v2xbG2h9/xL0efBALoL/4ontzgnTtKlRHW7Yoa5SOPGdcHK6bnEw0dCg6iAoVVQlaLYSb556Duic7G23+tdeIAgKIPvoIzhJBQdAITJiAuJitW903g3clqgThW4Iz/rf16mHvCo+Ujz9G4NWqVSDVjz+GQfill8p+bVvo318kVzt61PVZOwMDiZYtI2rUCMmtJkyAR5IKFVURBgMEwZdfhrdPTg4yzE6dCmHuk0/wX+vWSJHerx+Waty+Hba7yoYqS/jOVLYk4QUTld1la+hQomrViPr2BdHXrw994BdfwN/fnUhIwPSTyPmVvGwhLg46/UmTIPF4ecEFToWKqg4vL6IuXZCKZOtW2LWyskR0/aFDcNS45x7MCNq2xfeffhKeejt2wKtu/nzkCSrPmUGVIHxLErnR6FxFhodjf/162cqk12PkJyKaNg37V16BLrA8FjLhsQeLF7vHh95ggCSzejW+338/0eDB0IGqUHG3QJKg3hk0CLEwJ05g8fWlS0HqfBYwYACCGqOjkQNo3jyi0aOhHg0MxCDy2WfuL2+VIHxrrog8iMqcXG2RLSdKa25ZmzYRnTmjrFyjR0PtMX8+zgkMhGSwZQtcHt2pCvH1xepb8+djtuGu1YN69cJzvPIKBpfQ0KqZefNuDNJR4Ryio2G3e+89+P9fu4b9e++B2GNi4Ml34gT6ytChUBW5ap0LW6gShG9pxSsi666ZtoiWB05t3mz5/44d4Ws/ZYqysn34IQxB/PgRI8R/r72m7BpEzhFOXBwGmIUL4WO8aZPj11ACDw88C3fXHDIEks/QoVVv4Qp1EXMVjsLDAzr+p58GwZ8/D4NwnTpEjz8O+97u3XANdTeqBOFbk/CtEb6t4CTeoQsKbJPVnDnKVCUxMTDw/PgjkpFptRjxW7cm+u9/7Z9vqWyOYNo0SPnnz2Ow4pkx3YF27bBgxbBh+P7tt8jRs3ChKiGrUGEP5SFMVAnCt6Zvt0b49jxX+Ixh+HDbx/36q+3/OaZMgbQ9fjyMyf7+kIIzMogWLFB2jbLgySeRXKp2baKnnoINobjYPfeqVg1BZtu3w22zVi3o9lu1Ut04VaioaFQJwrcmaVvzQ7dH+Fot9r/8UtrwGx5O1KwZPt9/vzIp38sL07aDB4mmT8dvI0ci//xTT5U9D78S3Hsvwsqffhpqpj593Os21qIFPJT++QeLz+zeDXtGly7IzaOmrVWhovxRJQjfmg7fWji0knwzSUnY//VX6f+aN4fUTIQcM0rQuzdmDG+/jWvqdHDl8vNDSHd5LNmm1cJwNG8eJP71693vK6zVYlC7cgX7DRsQrZicDPWPChUqyg9VgvCteaBYcxFUEozUpg2Rt7f1pGeffgoVyRtvKNdPf/ghDL6PPYbBKCpKpCioU6fsiduUYtQoorVrYbw+dw5eQ+5GcDDq7MwZZCs8fRqxAl99paZeVqGivHBXEr6STHh6PVQ2y5ZZ1nfrdFhSbfdu5NZQAl9fuC6mpUHaZwwpF2bOxP/Dh5cf+XXuDPWKRgNja//+7nPdlKN6dQRr/fMPBrmRI+GXvGaNathVocLdqBKEb826bS9/uz2r+EMPQf1jzUVzyBCkVZ0+XTlZNW+OpGorVyJwiQieNDNmIFhj/Xpl13EF/PxQ/oQEZARs27b80sI2awbf5IULUce9ehE1aaJ88FShQoXjqNKEby/q017gU48eCMSSq3XkxK7XIzfOrl2OEdXEiYi8e+456NKJkJPG3x//lWcqVq0W6WFnz0aYd8eOyOxZHpAkqLdOnIBnT34+np8Ig2xlzEWiQsWdjCpB+NZgj/DtuQl6e8PYaq7WkQ8wXMp/4QXl5ZIkrLKTmIj89adOIe/GypX43Lt36YXJ3anu0OlgSF20CNJ248bICV5ehGswwHf/0CHYOYiQc6hOHah/ysu2oUJFVYerFjH/SpKkq5IkHZT9FixJ0jpJkk7c3ge54l6WYE3vbc84q0SSHTgQqiG+gpU59Hrkhz9wAEupKYWvL4KxiCDt5+XBTfO775CUqV+/0jMQVwRmDB4Mo6klt8h+/YhSUzHbWLSIKDISGQLLC1otZlVEiF2oUQNliY7GXl1iUYWKssFVEv43RNTd7LfniWg9Y6wuEa2//d0tsBdRay4d84AqS+vMmh+bmIi9rTUvx46FlDxpkv2yyhEfD1///fshzeblIQfH/PlEv/+Oz85K2dZmBAsXIplTp06WjwkMRIj3smUoT+fOqC9HBjNXoFEjqHX+/hszno8+wszjm2/Kx7isQkVVhEsInzG2mYjMNc99iYjHkS4gon6uuJcjsKYLb9AA+9RU+9eoXl2oW6yhYUNkxlyyBETtCHr2REDSpUuQvo1GoieeQOa8X37Bb6703ImIwH7rVgw01gaG/v0xQ5owASRbrVrFJEVr3RqD1Lp1SMw2fDhUQBMmKE9ip0KFCsCdOvwIxthlIqLb+3BLB0mSNEqSpF2SJO1KdzLTli2jrS3dt5IUyF5e0NMvWQKjpjVMnQq/+m7dlA0kcvz0Ezx2li0D2ZeUILXqrFm479NPO3Y9Dmv1Mnw43CGPH8fM5euvLQ8q/v6Q9jdsQFK5qVOdK4crcO+9RHv3Ih1zfDxmKfXqwaV08ODyiVZWoeJOR4UbbRlj8xljKYyxlLCwMKeuwVMhWIItkk5Ls27YlZMlz3Q5Y4b1a3l4iHzWAwZYP84aJk0ievVVSNN6PdRUU6aIRUaOHkXSNVcYMA0GJG5r1Qp2ghEj4At/8aLl4zt1gkdRWhoCzdy5YIOtAVqS4L559CgGqrFjETS2cCFmWfXrYzArjzSzKlTciXAn4V+RJCmKiOj2XkFCA+fAUxpbgrm3ixyMKTNKxsbCC2fZMmTQtJYHpm9f+OTv2weViaN49VWkFC4pgQH31i3knXnzTdx3zx4Ye209kyOIjUW0cEgIVsZq2NB6ZssJE6BLf/llGJbdDXsG6pgYovffR138+y+yl2q1wvupSxfrQXMqVNytcCfhrySiobc/DyWiFe66kS1yOH269G+c0Hx8lAc6PfOM+PzTT9aPmzwZXiUjRzpOzJIE9cqsWVBfdO0Kz5QXXiCqW1cc9/jjrnHTZAyziYwMDFLx8VCP9OxJdPKk6bEhIXBj1ekwCJWn944teHpioJo8GWqdtDRELp84AaO3Xo8cPvv2VXRJVaioeLjKLXMREW0lonhJki5IkjSSiN4moq6SJJ0goq63v7sFtgif56qxhObNlRN+aCh85omwOo01Pb2vL0j7yBHndN6SBPXJ0qWYJYSHg6wiI6GvfuUVDDgDBljPBmoL5gMFr7u6daEeef99eMYkJaH8166ZHvvVV5htdO5M9PDDlY9IIyIQuZyaigGqcWOsAdCkCdRWX30lVkJToeJug6u8dB5hjEUxxvSMsRjG2JeMsUzGWBfGWN3be7fFj9oifEupAjjpbdwISdBcmrWGX3/F+pVaLXTZ1nDffVCBfPSR84TYv79I6dCuHWwNOh1URu++C6+hFi2QT9/RTJvW6kurRaTrkSMg89mz4S768cfCFXLwYBh7p07Fqj1NmsAGcO6cc8/pLuh0ULHt3Yu6e/99DJAjRyKVxKxZla/MKlS4GxVutHU3lKQp4IuN24NGA1Js0ACrOR07Zv3Y55+HOmHKFOf9xtu1g3GyWjXo2DMzQdaTJhH98QfiDIYNAym7Mr1ytWoYSHbtgrpk/HhI/CtWYLD09UWa57Nn8Xx8CcXRox33UCoPBAdjIDt0CMFu0dF4P7VqwVaSkgLJ31ZktprYTUVVQJUnfEuul+b5cBwly2bN4KmSkGCdCKKikHf+jz+QjtjacbdugbBffNHyMXFx8Mf38YFxctYs/M5z3gQG4nvjxiAtVxJT06ZQea1ejYGuXz8Meh9/DMN1YCASwZ08iWf8+mu4Sg4dats7qqIgSVCFbduGd/7ssyjn7t2Q/KOiMCtYv966J5K6pq2KOxlVgvBtkZw9fa2XF9HOnY5Fknp5wT2QCC6C1jB8OFQw33wjVroyR3o6yOfNN0GmlhYFqV0bC4aEhEAy5V5CMTFE2dkwTKekCNKyFRXsKLgr5L//wu1UkiDx16+PQCyjEYPSJ58gEGrCBNgfkpNhf/joI6IbN1xXHlehVi3MUrKyMFPavRuzgI0b4fMfEYF8QkuXmtoxVKi4k1ElCN/aildElj1l5DlZ+Axg1SrH7jlhAvb2VEavvAId94wZRF9+Wfp/Plj5+2PfoIHlICKNBqoIIuTpP39e/FejBqTS119HdGz//q53R9TpEAxmNGLxlMBABKQ1agQjckkJyvfuuyD+Bx7A7GnCBLh/TptWfqmXHYEkEYWFYTYzezbKvmgRBrl165AiOzQUWVGJMMCqUHGnosoSPg/GKikpTfqzZ5ueGxgIslQSecvh54e9veUSJQleIt26QcdtnkaZE/4HH0DK9PFBXvrffit9rbAwSNl//AGp+ssvxfkaDUhp1iz462/ejFmLPRVPYSFmIps22X9m/jz33Qf9/pIlGAAGDIDEPH06Zhvh4RgELl6E+qRLF0jT1asjkvjECWX3qggEBRENGgQbRloaPJeefVZI+ffeiwRvs2ejDtwZhKZChatRJQjfkl5V3hGtZWzQauFzXlAAcgoIsH8vTqBJSZBoX3jBfvSrXg/VQKNG0IOvXl36epIEKXPrVpB5z57Qj5sT9pgxwvPniScwe5C7Zz73HEhWq4V6qEMH0xmD+fVu3oTKqWNHuFpu326/DogwwAwciGsvWgTd/YwZUD/16SM8n1q2hKH0+HGonL77Dsd26IABqzKrS7RaLHX51ltQTRHhmc+cQT03bw41W58+mMmpEb4qKjuqBOFbQ/Xq2NvSIY8YYToDUOJRI0mQ8P/3PxDZmDH2JWk/P+SkqVMHEatcrSQnfF7m7duhSpg6FcFQ5mmek5Nhm7j3XpC1vz90zt9+C1VOy5bw8KlWDYTcuDGuxdMtWxogX3gBevp77iF65BFIr0qg00Ei/v13EOHzz2NAq18fA9K2bXjGOnWwpm1qKlxa09Lwf2Qkzl+z5s6Iin3+eRD7pUsYvAYOhAfV66/DaE4kUlGsXo2YDRV3Bu6K2RpjrNJszZo1Y84gKYkx0IrlbdIky79rtYwVFDAWHCx+CwxkTK9nbPx4y/cKCWHsP/8R32fMwHkvvaSsrCtX4vju3RkzGhk7cQLfv/3W9DijkbFXXsF/AQGMJScz1rlz6ett3cqYn5/pc9WsyViHDoy1aMFYRgZjTzyB3+vXx/M99RRjr77KmJeXODc9nbHLlxkbMEBc5/77Gdu5U9lzyXHgAGMjRzLm44PrJCQwNm0ars9RUsLYtm2oS17/4eGMDRuGz//7n+P3dSdWrEC5du+2/P/Nm4yNGoVj7rkHbYiIMUlirHFjxqZMYezw4fItc1XEkSOMvfYaY998g7afleWa637wAd5XUhJj/fvjfc2bx9gffzB25gxjxcWuuY+7QES7mAKOrXCSl2/OEn6jRrYJ//HHTb+HhQnCZ4yxZ55hTKdjrF49/K7TKSf8khKQGxFjDz6orLwzZ4pyHTpkmfA5du0CEer1jDVpYvmYoiKQSefO4hk1GsaqVxfH/PYbY5GR+C8mhrGmTU3rJCNDHHviBGOTJzMWFCSI3xrR2cL164x98QVjbdoI4jdHcTGO+/lnxh54QBBlTAxjb77J2PHjjt/XHbBH+IwxNns2jsnNZezWLcb+/BPk1Lkz2hQRY3XrYsBdv56xa9fKrfhVBhMnlu7fYWGMtW2LfvjOO4y98QZjixYxtmOH8gFhwgRcq1cvxuLjGTMYTO9hMAh+qFULx02cyNjUqYzNn4+BITUVfbEioBK+bOvRw/R7x47s/6UvxiCREjH29NOCLJUSPmOMFRaKa6em2i9vSYkg/erVmV2JdudOxjw9xSCRlmb92OvXGfv3XyE137xp+p+PDyR7/pyBgZgBlZSUvta1a+g8nPj79mVs40bMPhzF5Mm4xt69pr/36iUGy+++Q8chYqxOHfGOevfGM1UkHCV8c1y5wtj776MtarWivcTHM/bYY4y99x5jf/3F2I0b7nuGqoBx4zArPXoU7+Sdd0D0bdsKQc58CwpirHlzxgYNwkz8668xu8zIEO1+wgT0BY7iYsbOnmVswwYQ+tSpaKMNGjAWFcVYw4aiT5prDWrWxCA/ciT6+fffu7/9KiV8XQVrlFwCe8Ew5iH0iYnwt2YMPusNGsBgy710Skoc0yfr9fA8qVsXOt4ffrBf3mnTYPTlq2/ZeoaUFOjh09Jw7VVcGEU6AAAgAElEQVSr4P44YkTpY/38EB0bFQWX0cuXYUjl/3l7w3C8bRsiT3NyYFSNjYVRlbuHEuHziy9iScQPPsA9V6yA0XXKFCRxs5WpVI6BA+Fh1KQJFjXp0wfePTwyd8sWlIN7V7VtC6PuihVY53bVKjxv9+6wXQS5bcFM9yA8HH7+EycionfXLrFt3CgykGo0aJ8pKdjXqoXPcXFq0BeHXg/bVnx86f+ys9EXvbzQtk6eRJzLyZOwjS1ZYrr2g68v+kujRqbX0WpR53FxsMlYQkkJ7pebC++01FTT/apVpl58ISEI1qxVS2w1a2Lv44Oy6NzNyEpGhfLanJXwk5NLj7RyvX5EhOl/XLomgrqBMUivRELtUa0aJGJzWJLwOV57Dee+9ZbysnOd9aJFto9r04axLl0g2SQm4pzJkxnLy7N8/D334JiHH2YsO1v8HhYGlcL998M2EBfHWPv2ONbPD9PUv/6yLPHfvImZSOPG7P+n0i+9ZKqbt4XMTMZmzRLvS6uFlNS/P2YN//wj7A18a9mSscGDsQ8MFDOT1q1R39u3Q31SWKisDM7CEQnfUruxh0uXGFu1CraVXr2gxpPXg78/2uh770HyvFsxbhxmr86ioABqwuXLMePq3x/126EDZgKuxo0bkO7feYexESOgXYiLw8zVnLP69XP+PnQ3qXQsEb5843phS4RPhEZw5Ag+e3iYvgxz4rNF+MXFjD36KM774gtlZS8qYmzxYstqADk44TMGkn/qKVHGqVNhdJWjWzeUlZPFSy/hmNBQnNurFwg0MRHH79gB1QK/ZqNGjJ0+bbksJSWY6nbtKgi4e3cMBnJbgC0cPcpYu3Y4f9Qo8fvx4/ht1iyok5o1M31XPXpgap2SUrrTJCTADvDjj/br01EoIfw5c5wnfHOUlIDY161j7PPPUUdc/detW9mvf6eirIRvjoMHTdtQSgra0MSJjM2dy9jSpVD/XLgAgceSIOQMCgpgK1u7lrHPPmPs5ZcZ++UX5693VxE+l3iVbuaEP3AgrjN5MkhEoxH/rVljeq/gYOuEzxgIvEMH9v8SuKskT074v/8OIn/rLcbefdf0OR59lLGTJ3F89+7w0tm7FwQpSdDfe3pCYu7Z05TwOc6fhzRNxNiCBfbLdeQIPHBiY9n/S+3dukFPKp9ZWMPFi6Z2Bk74CxeK3/btg51h6FAMXkSMeXvDo2jaNMYeeYSxVq0EIRKhjm7dsn9/pShvwreG++9n/z/z6d0bz//f/8Jj5W4wArua8BlDG2zTBsJet27wZuMeZuablxeMto0aob9Nn462unUr2sbp0+6fbVqCUsKvEjp8R4N3zMPj161DZGhODl4rEdIGbNwI/+oePUyPt6VL1emQ7GziRCwUsncvgq6Cgx0rozUcOgQd8AsvwLf9m29EXpjvv4eOf+hQ4XPfuDHuf/gwonC//VYsRu7rW/r6MTHwHw8LQ0Rvjx74bA0JCVhw5PXXYRdYuRJlGD4c25AhWJzEXEfKwdNF2EJyMuwHRERffAF9/6JF0HsvWwZ9fteusEEkJ8Mf/rPPkGLi0UdhD6hRo2rowCdNgn76zBm80zVrTP3HY2Kg+69fH1t8PFJDxMdD962iNKKjYVs6ckREuDMGPjh3DtvFi4idyclBDMa1a1g3YtEiwRkckoT4krg42MZiY3GP6GjY+RISRKR+uUPJqFBem7MSPldduGLz9sZ+zBhY2X19oVflCA6GlKEECxYId7xnnrGub1cCLuG//z6u9913kES4vWHtWuiBJ0zAPfV6eLqYIzhYeDVFRJSW8BmDPp17z2i1UKMsXKhcci0pYWzzZtP4hkGDMD22p2qxJOFbw/Xr8IAYPtxU520w4L3J32tUFGYE774LF9WTJx3zNqosEr458vJQZytWYNb3+ONwueWeWHwbNqz8yuRO8LiNynLd/HzMcletggryo4/AFyNGMHbvvfDCMn8XRJgRd+8OXvjyS6iNytJu6G5S6YSGKid0TsB802hMyYH/37QpXDP5761a4V6OED5jMIDK9eLOkr454WdlYer40Ufi+vPn49idO0UjW7asdF2NHQvXzj59LBM+x4EDsA9wdQ0nzkmT8J89lJSgLPJ6HT3a9jmOEL4cRiM63vffM/bss4x16mT6Xj08xGDON09PGPefeAKd7tAh64NAZSV8azAaYQNYu1a056qA//wHAt6dcl3G0A9yctC+fvoJA8Ojj8L5wcNDtMe+fZ2/h1LCrxKpFRiz/p98LVgiTHM5JAlJzby8xG/cHXPfPlPXzK1bRUZNW/czR9u2mP51747UBRs3Kj9XDkv31OvhMpmfD/fIFbdXDU5JQYoEjQZqHEsZHiMiiAwG2/dMSoKq6MwZor/+wtqwly8TvfceXNk6dkRiuKNHLZdPklAWxpCzp0YNrBHQsSOSjx054lhd2oJGg6nyI4/g2hs2YPr977+456BBSDXBodNBzZabC1XYyJFwz42Lgxrq889x7p0abq/R4Fnuuw9pOvbsgeqvTx/Ux5Ytpllj73a4qh1agiTB7TsxEVlkX3wR6si9e9Evjh/HcpxPP+2+MnC4XYcvSdIZIsolIiMRFTPGUlx/D+v/paXhf/5CmzYlOnBA/N+xI5Yu1GpNO3dJSelMmIMG4bhPPoH++623lOmF/f1x7G+/Cd26M5A/hxweHsif89dfSPLVsiUILTycaMcOYT9Yv965hq3RYOBq2xb5cNLToSOfNw+DABH0lEYj6tKSvt5ggM7zs8+g53/uOWzh4bARjBuH5G2uhFaLgalhQyzQQoS8RH//DcL7+28sIsPXFyBCe1m6VNg5/PywlCTXuR45goHF29u1ZXUnZs+GoHPgAOxV8lTgISFC11+7tthq1brzYh3Kioqw8Wi1EErNBVN3obyMtp0YY26TJ2y9qNxckB+X1uWBQowJcpKTvSTBwLJsmem15IupzJoFIpw5U1lD4QEVzi53aA7ze86Zg/TLPE+/VgupgmdxXLoUaYq9vUVgGZ9MOoqwMFz35ZcR1PLHHzDUbtyIxGhLllg+Lzoaxt3XX4chbM0akM+aNWLgCA/HfssWDCI7diBjZfPmrglKiYiAMbd/f3wvKcHaAsePIzjn7bdNA/UYw3e+qMzjj6Pu69TBDIgPKElJ7pUSy4K4OLEAT0kJZmzHjqFdHD2KvXmQEBEIPyYGMzM+CNSujWCh2FjLRn8VlRxK9D5l2YjoDBGFKjnWWR2+eZCKJVcq/nn0aNM8GVzHbL716SM+V6vG2KZNpv9zt62mTZEHxh6OHsXx33/v1COy1q1hBHrvPVzHkstjSQn8ij/4AAFk3L00PBw+xVwXb+6/zo2Z3KXTGRiNuNbgwdBXOnrugQNIVtW9u+X3YTDAGHb1qvNlVILiYtgRfvoJel1+f15ntWrBDhIZWTppXUAA9nPmII9OWcvqKp9vpcjNZWz/frTnOXNEvEZSUmn7BxHcehs2hIvvk0/CRXT+fMQOnDhh6m7rKrhL1/7UU3ivdypIoQ5fYm4WSyRJOk1E2UTEiGgeY2y+2f+jiGgUEVFcXFyzs06szxcRYX8hEo7Ro7G4BVetjB6NlMTTpkFq4StJ+flBr9ajBxYJefFFrGXLJUMiuADu34/P990Hd0xrUuipU5AKv/2WaPBgy8esXYtjeCoEOdq0gXTeqxdc87KzxXq2ltC9O3TYr7+O1Mh798JWUVSEtAuenmIhkuhoscTjyJGYJTRs6PgUNykJbqMGA1H79ihD9+7QXdq61unT0JsnJiKt8+DBWPLxhx9Qv56epqowLy/oQps3h9tp48a266IsOHAAS2D+8gtcd5OS0EaUugJz3W2TJrAR8C001PZ5BQVoB1lZOPf++/GMjz1mmv6ivMAY+tipU3hf589jvYXz58XnrKzSNg9/f7QvgwFpyKOj0beMRtRBRAS2yEj0v7p10T6ttZdx44gWL3a9/WHsWKT2UMojlQ2SJO1mStTlSkaFsmxEFH17H05E+4movbVj3SXhy7fRo0snPVq+HN435pLvuHHI2siPr13b1Kru7Y2kTfz7sGHWs+WdOYNjvvrK8v83b4rrvPMOIvHk4BI+D7ayF9TUrRuCcxiD1LpxI6Qj/owBAZBW69fHMampCO3mib3q10eK2PXrlUuahYWIwH32WdPUFjExkACXL4frqDl++qn0e4qPR/kCApA+YdSo0hK1XPru1AkeSzt3uifwZfly3GvPHvFbcbFIwHX2LALDiNAmeICYtS0gAJJqbCzSYDz5pEj5m5aGYCB+LE/3wbcaNZCc64knEES4aBHq6OrV8p8VyFFUhHa+cSMC7958E55uDz6I9hsfj+jq2rWRYKx6dcsJyPR61EuLFvBcGTMGaTQ++wxtwtMT/fjBB3EvV+Cpp5Aq5E4FVRYJXw5JkqYT0Q3G2BxL/6ekpLBdSlfekCEysvQiIX5+pitBcYwaBSnb3Hiq10PqkCdWIoJRq2tXSJuMQW8vPyYyEgEwU6Yg2VfbtgjGiIkxvc6lS/ASmTdPGBDlyM01ldwiIyFtDxwIO0PbtpDwe/YkeuYZSO+2VujiEv62baWfp00bzFrWroVEdfy4+D89HbaLRYuwTCIRJK/u3THb6dhRuTHvwgXc47ffoKfnNpCYGBhC27SB18iePUQPP4xlFg8eJPrPf6A3PnNGXMvHB0bd3r1R9sOHkQzr8OHSS1j6+WFG1qQJZmEREcrKawsrViDp3J49uK4lzJ2L5RCvX0cZMjMx4zl0CM+1cqXlReptQafDMwQGov1ptWi72dlo8+arufn6Irldw4Yi2KdaNWwxMagLnqCuImA0IhHf5ct4p2FhqCutFrOaq1fxbGlpOObyZfQda2tH+/nB0y4yEs92+jRmX0FB+C0qynQfGWnqlcfx1FNo91Vdwner0VaSJB8i0jDGcm9/vo+IZrjznhzWGrU5oXNYMqbqdOi06ek4LzhYuDi2bAnCSUtDR//vf5FVb8wYqBgWLID6hcPHB2VavBjZIpOSTO/Fx905c0DwH3wAz56ZM3Gv3FwMNpYWZSeCwbFaNcuN2RwXLsBd0lJ9hIXhGcaMQXThN99grd3FixHlSgSyTkpCB23XDt/r1MGgFxYGwo6MBMGMHIktLw+urjt2YNu+HSqSyZPFvd96C2RFhCn2unVQIbz+OgaD5ctNPUx0OhxfuzbqNzsbqqurV5G11PyZ+vbFAFdeKpGQEKi22rfH948+gvojLQ0Elp8v3gNfh/jgQQyMZ87gPfXogXrmhmNL7/+rr9B+liyBMf7UKbHalrmKRasVg0BMDLbAQAg8BgP+q1ED/2u1+M/Hx3UeLIcOmb5zS9DrYbxPSiJq1Qplql4dy2K++CJUL8XFIPv8fDgL3LhhKsSFh6PfWpJnAwJKDwaff66s79zpcKuEL0lSLSL6+fZXHRF9zxibae14ZyX8wMDSOlVzN0uOzp3hAVJYaPla5hI8ERr8jRtCSicqrVdesQLS6vHj8Hvevx8dcOBAccx//0s0fjzOGz5cLIlXXAy/3GHDMHg88wx+P3cOPuJz55bWWd5/PwapnTtxPS49jx4N99G33kKdmEv45nWl0UC6HjsWswhrKCpCLMKmTSDdnTutH0uETtSiBTpsq1b4bJ5K+fRpPN+PP2IwiIgwnanp9ZBYp0+H62BiIqTn3btFiPuBA6hrbnshQocODsb1LcHfHx4nycm4blISvoeECC8hczgj4bsajGGgOHsWbWPjRggGRHiH1tJy9+iB93vhQunt/HlT7zNLkCS0m/Bw65ufH95XdDTqPjjYsj1r924Mct9/j76RmYl3fvWq6f7MGfSlS5fwG2MQFL75hmj+fKTI7tMHAztPa26OgAC8ax8fCAXVq6Md5+WhP2dkYPC9eBG/xcaWTqV+p0CphF+uKh17cCXhW0N4OBoZHwysqX6IQDCHD+Pz448Tff21kCo45ANEWBiIu107+DVnZYm8L1euwB9aktDQjx8X0/vff8e0lEiokCZNAkkS4X6NG6Pxx8dbX282ORmG2Fu3IK0lJaGDWaqrZ55Bx5Gv9xsZiQ7QowcMom3aiLUCHnoI0vHkyXDzy8lBvRUVCYn93DkY1bZtg6vftm3CMNy1K57TEpYuxfUPHEBdJidD7fXzz6UNgR07gnh79MCsgpf/4YcxiAcFCXfayEgQRe3aUIdlZ0OClvvdmyMqSrgcpqUhbiMtDc+6fTsG5ZYt4ZpYsyZIxMMD57qb8C3hvvswE+I4ehTvISwMZD5qFJ6FvwdL4G66+fkgvzNnhOtmtWpoL9eugXjlW2amdVdUjUbEWISFoV5jYnCv2bMxA5HPgG3hzz8hqG3YAOFg/nzhJEGENpiVhfJkZKCfHDuG37KzMZO8edP0mno93l+NGhik+Ay+Xj0xYAUFoT+GhOB9VuZcTHcV4QcFObZYtJykq1cXPtbm0GrRoPmx8fHoEJxkYmMxfW7UCB2No3dvSIEzrCivOnfGQHL5Mr6vWgVpxRJ69cLCH7Nno9F17w7S7d4dki1jIKSQEESZ9ukDffn77+P85GRIVM8+i4AhTviLF4NojxwhevBBlLV5cyG5WwvyIkIn6dBBbDVris7AmGnHyMhA1K9OZ1pHcnDCP3gQA1W9ekQLF2KGdPYs7ABHjuC3X37B89qDhweIJjfXVBioXx+DT2wsYikyMpBsjnuSHDgAguDQ603VffKYDl5PXA2Sl4cZwMSJGGRiYnCfxET3BmqdPYsgstBQEKEcgwfjP77QjCtRXAySvXoV6pOiInzPyoKAk5aG39PT0dYvXBB1uWULBAol2LYNs8Q1azBQzJ8PAWn1amXnG42mA9WVK3jPqakYCLOysJkPCnJotRg4AwLQh/gWEIB36++Pz/KN/883f3/32U8qhQ6/vOBo+LtcZWMrEMr8useOCTJbvBidS68HGe3eDV3trVsg8PXrTc+VJKhM/voLkootgysRSCQwEATH3T3r1xdl0ushhezbB8n522+hY//4Y6hGwsMhyd68CWPyl1+CEPj0XaMR5erbF4SfnAxpnQhS6o4dGFR27cKM5KWXIDVu2oQyLVhgueydOmHQmzQJ92zbFlKaNfCBxZIEJUnCda9jRxzLXSWPHsWzt2uH4955B+U8cACRsr//DhfTV18FEf/zD2ZgH36Ia0dEQG8rfw4eANa6NY7n7SMwEELF8OEgcIMBhHf+POr2wAExc+BqFg4e9Sx3QeSfo6MxMFSrhndpMGD24og0Wb06NkvgA9SyZRj8PDzwTvgWFqZ81TJL1+bPoQQlJZC+c3Iwc1QKPoO6dEm0f0fqR6vFDCMqyvZxBQWYEfABgM8aMjPx+40bEB5yckQmzZwc9LHcXGVR9J6eaEsREaj7oCCxpaRA+HInqoSE7+tre3S2BXOdvfy7l5d1I+nYsfAV/+cfEJGXFwhn6lTT43r2hHRz7RokCm9vXJMx+J6PGoWGZUnCr1cPEa1aLaTQwkIQx9WrINTTp5HvhQgd4Pp1EPHy5ZBS5ZJoeDiIZc8efF+7FpLS0aOQrKtXx302boTKxh5KSjDQbdqEbckSSPJym0FUFH67eBH3+PprqIXq1YNelQgdZtUq1MWhQxjIuIT/ww8g1L177ZfHEpKS8Hz16wtPpyZN0F6WL8cAfuIEBgnzfEMJCagHroc+edLUm8kSzKX/oCAMylot/mMMpHLjBu5nTdiwZSdwFE88AUnY3ItNDh8f00EgMBAkGxgo7BpcxVG9OmYv5enpc+YMZpFy9OkjckdVFhQUCPWXte3GDbR57mGVnY3v2dlQ9S5c6Ny97yqVjqenbb2sI/D3x2jtSLV4esJImp+PPPUdO5omSdNqsX7r+vWiwxMJtcmTTwoPmIQEEN/ChfBuMRhwfOvWKNf58yB2f39Iij4+GADGjcMMIyUF53Trho4yZAim0jwZWIcO+O+779DAjh0DGT/0EFQrRCDHVq1A1m3aQFfuiER18SJsBFotJPGTJ03/lySQh7+/qevlyJEg5eHDBeFfuCAGKUexaxcGvgMHMGhcvCj+CwtDXdepA6lx3TqUiTGQtocH6vnKldJtiyfD8vFBXWu1kMq5LWPoUFzz5k10Zq7OkNtLzOHnh0GUJ5kbMULovbl7ZUgIjtFolL+P0aMh3Wdmwvjdrx+eKyPDdEtPF58PHbJvyNXrIeR4eeE9hoQIfX10NAbtRo0wcwkIKHsu/n/+QVvl7XjlSqjmqgoYQ5txNn3IXUX4Op3zWQ3NpTJRFmEc5cSckgLVjXmVNWokJO2yIioKknNAAGYQL74I0ujRA5Jply7C5dDcU4gjKQnShJcX9Ldbt0Ii2rkTZP/oo3ARHDBAEH5REchq3Tocv3278H6IjITnT8eOIJy4OEzJlaoCsrORVG35ctgurlwBoeTnW3eT9fUVJDF1KiTV4GBBdAUFeA5PT3gsxcejXCEheJaCAlPDKWOYYR0+jGc+ehT7U6eELcXXF+UxJztPT5QlN1cYabltp6AAM7bcXOueX3o93mtEBMiRq/s8PND+/P0hTRcWojz2IEl4t0FBkL5jY/Hc1arhc40aQoU0YwYGzuxs0V4CA0HMXKXD95GRGAB798Z9Zs/GQHXunFj8g89MuBuk0n6n0aAePD0xy/XxEc8dFCTKUKcOyp+fL4ynwcHoDxoNZtZLl5aOP7jbcVcRflms59bcN5OT0QHWrhW/eXgIaW/ePOESZqvR5+aKoCxHEqcNGoSG//bbkMrT0iDt9egBA6xc9eTjI1RawcFwMdy7V5SLHzt4MCTefv0gQQ8YIHTe5igpweD299/wltm61bT8Wi3qZ+5ceMi89x5IICkJ0l14OCTLbt1KT8flYAyqqfR0nLd0KaTSTp0wiGZliQE2IAB1wtNZWILcyBocjEGgenXL+/Bw1I3cp5tDkoRLn06HQSAjQ5zDXfssvfugINMZAH9PhYVoL5cv49rdu+N8nsjMVYn1zOujpESUMzgYKUKKi8UzZGUJoyvHgAHwiLEHoxFtnOu1r14VLpUHD2Jmcf067sMH+YICx2bQHHygLS6GHSUsTBhG+cDBvWqCgzHI8mO4cGLuVGANRUU4zhUJ+8oDKuE7AXN9/oEDILD580FCcnC3skOHkN/EGgkNH45pqNEIgrTmM2wLljxmgoPRwUpK0JgLC0tLy7Vqgdx/+MFUnUEE6XjzZlzj0UfRgdq0EYFARBhgLl6E3r1+fdyzuBjHpKcjKIwI5GXruUJDRSBWZCQ6aEIC7Bvm9oJ//kE5OnWCkS4vD2qe48fxfIWFIIzISBDXpEkgrFOnII1eugTCuX4dAw33WT97trT7rYcHpOK4OKGC8/bGM1qT1q3B0xPt59YtSKk6Hcqbl4fNlgtjYCAEDB6Mxxh+q18fEnV6Okhs2DCxvN716yDUkyeh5rt8WbjgJiVBz56Xh4E/I8P6vSUJs7W6dQVZ1q2Lgdo8WtxVYAxly84WW3o62trhw3iec+eE9w9jGDx5fTpDWRqNGNy1WgyEfPP2xoAfEIBBw9MTQg6R8MDhXji1a6N+g4MxmHP1algY6i40FHsPD5SZO0e4GyrhuwjLl8OL5exZqDS4zjkoCJIp97b58EO4423eDPdIc4Lt0gWDQkQEBgkidOgTJ0qrlMxdAZVAp4Ov/v79pud6e9vXx8rRqBHcM3v0EB4DAwaAXI4dE8dxG0JAAAyM8+eDIG7dsh4Gbwm8wwUGQiLLzRX1Yw08yrd5c9QzDzwrLEQHCw0Fib/7rgik+vVXEKJeD4K5dUvo2HmOfGvl8/W1vIiMo+jaFW0jOhrlOnsWA9XatRjozB0P/PxQJ1lZIOWEBLSfuDgM0K1amR7/8MMwnlsKHuTw9cX97RmgiXAdb29BiPn5eN88SjcgAKQXFIR98+YoW2Sk8KwpK8x12zzJ2enTwuWWDwwZGcKj5vp1Mdu4fh11m5eHvlFUJNKDK4VWKwYMR87z9sZg7uODz76+qMugIPzu4SHSZVSvjhgPZ6ASvhPw9bVsWGvQAAt/vPQS3Co5tFpkbRw/HpLU00+jwX3/PX4bPx5EKDf6+fjASLt5szJjJJcUrXVgOazp9JX+bwnyKMvatYWHTU4OSDQ1FQR88CBcID09od5JTYURNj0dEhuPFDYY0Mi5pG4LPj44xmgE0Wk0IqunOaKi0BGTk0WEcd26IKHYWHgIWcuTIrfjDByIZw4NFWqPa9cwMB85AmP3wYOIc7DmwVWvHp4zPx/nWnLZ8/REnXIf7aAgvOMOHXDdCxdAUteuiQCjsDBTaX3QIMxyuFpj+3YIHq1a4Z72vJskCeXw8MCWkwPbAB8wiotR93yvpA1yaLWoAy8v7L29EUjIffS9vUG83t7CVubnB0HCzw+bry/KFRAgiH/XLpB4//74HhODmTYfmOSbjw9m4tOm4bunJ67Dr+/nB1VpRgbaAK8H3r4jIzHgFxZiEOF2G952uQ2juBgCBO+rubn47cYNxwS3li1LR8YrhUr45YwRI5AqIT0dZDBzJqaFEydaDtfWaEBEzZtDwrMW7WsLvGPqdJhG+voKo58lKc9gKK2qsBVgZQ/+/sJwyRs2fxcBAZaD4Xx9xZS3oAAdTomHldxOERWFerZkbOdh9PaM6Ho9OrynJ87hkai2pGMikGv79nh33JjL8zyuXg2VRHKy0FdnZFgfGIiU17/cs8lgwIBq6fnliIwE+aemCpVjTAzKfemSIHJnHR7CwiA9W3p/8kA8JeBGXe6BxNU+1uBou42Px6ymIulOo0Ff1WpFP+TPwRgS3u3b59y1VcJ3Aby9IemVZ36N4GBTlYgz6h0iENPNmzjXmmG6IiBJ2MxJlTd8Lm0HB6P+b96EmqJaNfjrf/opiCw8XHgicV9/vmqTRoPn9fZG/V2/jmvHxuK306cd19E7A09PbFyS5p5HJ04I46G8LozG8ickSRL1VVb4+qJetTAq4ZAAACAASURBVFoRe0CEAaGgwLEZAge/Ft/y89E+AgJM7UacA5TcIyAAZTUYoFaLiIAzwM2bYsvLE4IM94riBnB5BD43iPON/8bfIxduuLHaVryQNSFJCVTCdyF69oREKU8YFhUF6c0dnhVKEBQEwnOmEzkCexJvRUE+EMbFIeeNJKFDcXtA375Cf1pYiE7OCTY9HbOL06ehJnjzTfzXuTPeKzeM5uU5ZgO5W8BdZuXEVxHgA5Z84zAY4LjABx+dTnhN5eZC4ndmZi2HRgNVb5MmWOozMhLtcsYMeMU98ghIPzAQs7OmTdG+Zs2CYMdzEmVlQcI3j9JWXg8q4ZcbuFogM7Ni7p+UBAnCPMCpqkKS0Im4m6MS6HSQ0nx9RY6bqChchzGkYNZqEaim0QipLzAQQT5//IEBIiwMNp1hw3DetWv4/cYNuK7u3g2Dd3Exzr1xQ+TRycmByufaNZTdnkqGw9MTMx5fXww+XB1lNOI6+fmoh1u3hGTJpUxr3dvcXdMWHFXPqFAGuVpKoxFBl85dSyX8CkF5ScS21DQ8EI2/2gkTYCeQJ36rrJJ7VcDDD0MVyCPAly4VnkRlmRFyaZa/V/P9nYSy2I6qMpytk7sqeVplgi0SlXsBBQWh88oDixyBLclMLjlqtUgkZh4gZq2clgy7cqgd1T6s5aUvK7inSlWA2oYqBuUUFqCCSJC9vz8yPPIl/Zo2dd89jUa4lCklCnvGzMrWUbmaxscHKpq6daGOMRgqumQqygpHcgbdCeBGbK3W1Cef20OULh1aFqgqHRUqVFiEOpsrX1jyXlN+rjKVjirhq1ChwiJUsncdNBr7KaV5ynC3lsPdN5AkqbskScckSTopSdLz9s9QoUKFiqoFJR5RjkbBOwO3Er4kSVoi+oSIehBRIhE9IklSojvvqUKFChV3IpS66ZYF7pbwWxDRScZYKmOskIgWE1FfN99ThQoVKlRYgLsJvxoRnZd9v3D7t/+HJEmjJEnaJUnSrnR1VQMVKlSocBvcTfiW/GdMTEGMsfmMsRTGWEpYWJibi6NChQoVdy/cTfgXiChW9j2GiKwkuFWhQoWKuxflsTC8uwl/JxHVlSSppiRJBiIaREQr3XxPFSpUqLjjUB5R1G4lfMZYMRGNI6K1RHSEiJYwxpxMD6RCxZ0BvR7LK/KFty1BDRYUcFay1euRephHqt7p4OvuuhNuz6XDGFtDRGvcfR8VVRt8ZaZr1yq6JPZRVIS0y7ZQFYKa/PyQmlqnw/MUFIh1ltPTie6/HyvCXbmCLJBr1iDlRdOmYhWpPXuw1oFOh/Nv3cLvWq1YGc3aOgEVlZrcXSgPCV9NraDi/yFf6u9OgLy8ahoAU9hLaWztf56/Rr6Ih4ryg5otU4VNaDREbdpA2rlxQ6Tg5WuyHjig/FqVhez5QhYeHsgpn5dnugoYJ3d5ed1NTtYGFI0GyfCio7GwBV9xq6JTT9urD2v/V3S5VbgXKuHfYfDxgWQbFob89nl5pgurE4Eos7JsL6dWUVBChnwFpeJi556BqwNu3UIGTV9fou7dsQIRY1ikJDUVapdz56A2uHoVhF1cLMjQy0usq2qLIHNyTJemU0lThTMoD02FqtKppOjWDettEsEodfEiiP7q1dLHVrQ6w98futzUVKKOHZET5OhRLJTN4edXejk5f3/lK1YphV6Pja8opdWifipivVgVlsFTVzuyrrC/v1g5LDOTqHp1vF8PDwzMBQVEp07ht5ISpBrmRtCCAtEueIpioxG/Z2SIdWq1WpRNq3Ve2CgLPD1tL9xuC+qKV3cgYmOhwlCqhomKIrp8WXzX69HQq1fH1rw50dSpRHXqYPlDvg4sl7ItDRTmul0PD8wm6tXDup0nTxKtWIH/6tSBGqm4GB2Hn1+7No5t0gTqpgYNiH76iWjUKKL167E8YI0aRPXrw5h38ybRv/8KyZiTtFxl4+kJNdWFCxhciotNB5TKCO6hUxESP19429dXDMjh4Zj5GQyYKW7fjrqfMkWUlS/ZWFwsFvBetQrf+/bFe87PB3FKEtpG3bpCnXj1KtGGDc4vxu3I81Ui6nIJfH2dX2NXJfw7HM2aoQFs3249i15cHFGLFkSvvIJOZzCgLg4eJPrnH6gsZs0CaRcUiPPkRFS3LtHzz2Nd3MRE3JMIHTsoCOSbmUn01ltEc+daLoeHB47T6yHlHzqE87OzQcqpqUR79xLt2weiP3vWdmflAddNmqBcx46BTLZtc14CchZ+fiA7T0/Ub1ERvvPNESlVKTQaSLNeXrivpyfKERwMoSAmBoP9/PmolyVLcPz8+URz5hB9+CGIOzsbwsPx4ygnX/9WieRaUYOVJAkJ2xbM3V1dUU6t1nQJSV4G/t41GrwTXpeuAL+nVos+v2WLc9dRjbaVFFzyWrAAq9T37QtCIyJ66SVIXi+8AD2zOYKC0Im/+AILZRcUENWqhQZ64gSWMvzxR6hT5JCTvSQRRUYSDRiARZOHDBEkz3HhAgaM48eJFi82XVj5lVcgAa5YIabYnEAyM7EouDnCw8XnRo0gZdapA2n/xg2isWOJ/vwTn7VaPGNxMdHvv2NzJTw8QKBGIzpySIhIXXvzJsomd4/jEpcrBhofHzxXQQHRxIkY/Fq0wKDm64v3GBMDd8Yvvih9vtGIesrNRV3n5qLMn31GtG4d3hsR1jAmwrs2GNDetFrc12jEwGw02iZJa//xdXXtuRBqNBh4tFoILQUFcMssLsbs02gEiXLy5ANpUZEyd0vz8mm1qF8uLJw6hf3UqWjv/N3m5BD9/DMGRPNrGI2Wn4uTu9EoVq1zFfg9i4qIzp+3f3xZUSUkfL5ot7shSeiYU6dCcnZk+uXriw5erx7RkSNEhw/j90uX4IdMBLWGhwfRoEFEXbsSDR1K1Ls3GqncS6VWLUjNHFxVYwvBwURnzkC90qUL0ddf4/eMDNx3xw7sDQYMRhwhITB4Zmai4WdlQVKXQ6+HBJqXR3TffXiGvDwsNShJ8MPetw8G0jlziF59FaSr1wtDqVIEBYlAm1u3QJZGI9HOnfht0yZ09hkz8LzuBCc/cxdGSUI5CwrQNnmAUKNG+D0tDUFZXO2RmYl2wNViPXqgHg8fRr1duQJbx7VrZVNjSJKod70e7bekBKoy7kefm2u7L3EdN1epMIZrFBU5VzZehzqdmCVqNBj0Kxp8sJTHflhzOuDPwY/h753bqIKDRRwJV4cRmdZhSgrRxo3OlfWuUunwKZe78dBDRB98AIkhORnTZVvVN3480Ucf4XNAAFHr1iD+AwdA+kTo/FFR+LxpE6Tzb74h2roV0t/OnaXvkZKCgBXe8AICoI65cQPX7tiR6LffcJ0xY6C2+fBDEHdYGIjl3nuh+klLM722vz/sCEFBYuaQl4eGHx+PhuvrC/VKURFIVz6D4NBqnRuENRqixx7D1Pb8eQwG5vpaX18MPvz5ucRZVt9xrRbPFxaGThkejsG1fn0MLPn5kCLr1oUbJu/gRiPUZ++/T/TJJ0Tffw9d9vHjIPKLF0Ha164J8pAk55a0UxIrodGgT+h0KLNS9Yi9snA1k5z05aotrjLS6UqrRtyhGuLqH64W4TEEfAZz/To+y++t06H8vMzOtNGwMEH8jzwCTrCEtDSiJ58kWr0a3zUaDK68D/n4YOOfW7XC8c7griJ8b2/lU+7wcMueLvZQsyak6vx8omeewTRajvHjIf3euCGMmtu2Ed1zDz5rtURvvw21y/HjRO++S/TUUzje3x+dw5wko6NBqA0aEP3yC9QbDz0EPe7585hlREaCaPbtEx1sxAiiL7/EZ40GqqIZM/B92DCitWtx7YYNxfb770JH7+WFRlmvHrZTp0Boly+jPEqh1YoISlsID4e0W1iI6X5BATqLvz9mPwUFpjMaRxEYCMKuWRMqk4gIzFxCQ033fn6m9qDCQgysO3cS7d8P1da5c5A+Cwtd6/3DicpgEITq7Y364wNZfr5Q4yi9psGAsrqym/MyuYPEOWlzqdcaOEnq9WLwNBohzPDZhtJnliQISQ8+SDRzJgbp++7DwN6wIfrKxIngjaIi9LmwMKhdU1LE7NIS8vOFKogxcER6Or4nJYG38vLQrwIC8JszuKsI39fXPS5UY8dCvdK+PYhv0yZ4nZiTT2goXiKRqcS+ezcMMUREw4dDd8in8ffeSzR7NtHnn8PYxl/DDz8QDRwIksvKguTO3Re9vHD98HCikSOJPv0UEmOtWlD/tG+P89q2FetjajREL76Ia37xBYgrN5do9GgMJPn5UOF89x0Ir3lzkJwtLwvuBudMh4+IwLPExGCW9OOPjl9DDi6Vp6ej49atC9VJly4YUOfOxQBtC4xhJvbzz3i3V66gYyt5Pls6bb0endjXF/XFJWROaNeugVysSZmenjhXp0N58vNBDLytBATYTzXBBxIiEdtgC/7+mN35+hI1box24esLgnLUM8ZetK9GA/IkEnXDBz0fH2G09vISNpZduyB4mL8bjQZtKzoaQkJUFJ4lNhYCUWoq+qlej2tzW0FhoVClcTzwALzK5EhKgs1L3l7btcNMdOZMvIecHDGjun4dzgnc405J3T38MGxmzuCuMtr6+LiW8A0GqDvq1hVW88JCkG92NhqVpyc+MwaitIQPPxSfv/0WjbZePTS+3bvhhWLpWc6fR4ObOhX35Z2a769ehdeMry908Y89ZppA6soVoo8/RuNhDJ/feMP0Pk89VfrehYUYUDp0wL127hT1yjs9kZhN2WvEPj4YjGJjYbAjEvr9Eyew2YLBgGmu0QgCDw8XRufq1TFoBASgPk6exCwkMxODGy8XJ9ZjxzDYnTpFtHw5OieX1K2BqwzkagDzGQ5j+D8uDiR09iwIRG685pIgR2Ag6qRBA0iK//5rmYjz80FksbHCO4dvyckY6FatIpo+3bIxkXuV8PKnpOC4ggJcOzUVn/398Z7z80FUXO/MjfX82rbetU6He3Ef+LAwGObDwrCFhkL1ERAAQn3pJaHL5t5ctqDVopwBAZhN16uHfsgJPiJC+NibQ68HUfv741m4776Xl/Buuuceol698B569rRdFg5ue3vxRdRvYCDK4O0t3GHr18fsgdcB3x86hGO8vcVgFxtr83YuQZWQ8KtVU+6TnZgoDKaW4OEBCZgT6N9/Q2LmSEqClNG5MwirqAgvmifLOnsWxllzJCQQLVwIrxj5/R95BGqC7dutS19RUehscn376NFEzz0HQi0qwjXnzYNtwNzwExKCMnLPBTl5e3tjwBg3DrOLVatwLaW6Te45ZA333y90mErQpAmed98+fF62DGR94gSI/cgRSE1btoDQ7LnHOZLmwHwAMxjElJ5Lo5xwOdlER6MOJAmziWefBWH6+mLgPXUK6ri33kJd6PWQUM+dw54xvIO5cyFgnDuH569WDYRRUoLnvnUL90pLw3l8v2+fbU8mSwFvziAsDLPbxEQQeXQ0yhkVhee3RrZytGqF8nzxBWbPBgPOlW8hIRjIOLnzzdvbfe7XBgPR5Ml4R9ZgScIvLoaAERBgmumSq994jEpxsbD/REe75xnuKgnfw0P5seZk7+1tKrXVq2c93WpEBAxyI0ZgBsDvyxjRr79Cp75pkzi+ZUsQOREItFMndD6dDhLagw/CzdGSwblBAxD4o49icGnaFOS3aRNUOTNn4v+EBEi38sGidWvcl5N2ZiaeMzkZaiLu/tWoEQaqr7+27AbIpZbcXCHVt2mD+uGDCif7gABh3PX2Fl49lsjeYIABeMYMEOrMmRhoMjIwK2rbFqqwX34BCcjtM56e+M1gAClyD6HCQvj6m8Me2et0ILAGDfDuo6IgadWsibq11baKi2GvOXMGLpEbNuD3pCTMIMyjiNeuhTAQEwOJ0stL5OE5cwbCxcWL2C5csG8vkQ9Q3HHBXH67eRPP5O2N+goLE8du2YIyVq+OAUQ+cNeogTIuXgwJ+IknXEe41avj3ZY3Skpg20pNxSCano52WlwMlV5xMdojt5Nw77ghQ/D9/Hn0vcxMseXk4PiMDGw5ObiONYGpVy+ixx9HXV+5gn1WFvatWxNNm+beOqgShM9DtZ3BnDmQNjjMc3Nz32YiqCXGjIGPOhEkEd4wLE0DOdkTQULt0AHEtHUrpvnc116vR+e/fp3ovfeInn5anFdSAuly/Hh837TJtLxHj4ppPydyXj45zp9HAzUYhBvnv/+WPk5uOC4uBiF6euIeN2+ClEJDIa1oNCDFc+fQeInwfE2agDBDQuBplJwMPepnn0Fa79ULpNivH87hRjd+/tat4j0YDCB8HjyWny+ii7lKyFwyN/+u1YJk69eHFF27Nkg+IQHEpiQfe24uCHLNGhDy+fOY6XH4+AgpLzIS5CpJqLP9+/F7ixaYiW7aVHog4l5B3t74z9tbGFx5mghz8GfU6TBTvHIF3lmffIJ3lJ+PgXnZMtQZn+FxGAx4h7VrQ+0QHIwBfsAAvL8DB0D4wcGuIfuKDpDcuxfR3ryuuXqFCO9pzhzL582bh/2ZM3BfJoKQERIiZiDx8SBsf3+0VT8/zPJ4uga9Hna2X34xHez8/VG/5oKNu1AlCD8oyP4x1qb2XLf988/QLcoj+A4cgNHW0xONgxtgmzaFYfPaNZCjTkf03/+i02VliWmuXGJ66CGodIYPB6ERQaf63XdE77wDvTIRSEiOkhJMDR98EB3XHFotjGuWArXMj7t1y1Rq5MYrrhapWxdSXUICBo1t24SqrE8fkN7x4yA8InTgoiK4e27fDonp3DkRkCUH12F7euK8/HyQ/+nTmC1xFdOUKahbfv3CQgyG+fnoWDdvliY/c6k2OhpE/vffmEW89JLjZJOfj9ng/v2Ypfz8s7hP/fqY7fXujffv7w9VHh/gOSnIwQfMhg1Rvzx4igdQ5edbVkt6eAhDanw8JO3iYqH6CA3FAKPTYZb2229E//lP6etMnw5i57aA6Gio8k6exOzKFlyp9S0PDXJBASR4vl25gj3vI5s2gZw5PDyg0hw8GETt54c27OODd19UBJVPnTrwsgsOdk7IbNsWZQkKwsAaElL+i7dUCcJXoqOV52mRT7cuXUJHkUubRNAVN2qEz97eQnpev55o82aQktEIAouIwOgtR3Y2OiLXu5eUgChOncKAsXq18ObhrmhEpsbGZctw/pIl+M1ggKQ/dSo66dixaNzymUJ0NCSVzp1NDYdGI+7Trh3qYMMGPNdDD8Ezp0YNZN38/XcMTIyh/N26gfhWrgSpdOoEF7XvvwcRcRXWkCE43xLZm6OoCAT688+l/+Nkz+tALvXcuAHJrE4deCQ1aQJjaVEROm2NGpCqPDxwrJ+f8Iqxh+vXobrYvh0kn5oq2onBAFuJTocB/cQJEUchB79P+/Yi8Co7W6gP9u+3HLQWEACJunZtDAgpKZgNREY6thpU375oz/7+aJM1a6JOIiOdk67tedk4ez05uFE9KwvvsWZNDH7p6dgyMvAf96bJycE5XN3CVWfyzVo0rJeX0MWbQ6+H4GSO3r2xf/VV1Cv3KnIGtWphq0hUCcJ3BE2bmvrNajTwTefSmZcXRuHERHwfOhRui15e6LweHsJzhfsB8w6RmUnUvz8+JyZiOswHjR9/xAxh+XI0IvlMQt6h/vgDhDJ9upghaLXwqmnZEoOTPFUBEQhm4kRI4dx4xj16AgJAPBER8Azw8YGaat8+SHrLl2OGwW0bNWuCSC9dwmCTng6Sf+01DIp8NrV6tWkHdoRQ3n8f5FezJraPPsIAxsndwwP12KgR3ldyMgYSa95QjqKgAOqYTZswwBw9ik0uCHA1Fp8BcXWIXo+BNCoKuujQUBEgdvAgyH3zZtP7aTQoe0KC8M7p0wezKVcuaxceDmJyFVxN+EYj2vdDD0HFdOkSNmu5oqxBr4eUHRyM9s0TBvr7YwsJQV1wtU1EBDbzFCIc5aVq4oZco1E8/+XLYmvSBDM4d8JthC9J0nQiepKIbnuo07Tbyx26HEos31yvm5RkSvhr14Ik//4b3+vWRWckAmkvWCA8ErgBr3dvGC+PHxcGV6MRkhkP5x87FtM2jgEDiJYuNW1cXGe/dq3we3/3XdNye3qCKDZvhnslh0YDop40yfK6qXyqyP30588H2ZeUoByTJ6M+JAnl7tIFOurTp0FaXbqAdHv3RsdyJfr2xcbx1VfopFu2INIwJwczCI6cHMcM83LcvAkX1507oaLbsUPETNgC976qUQNbdLTwsMnKgoRvLq1zVKuGOh85EoNlbGzF66+dgavL3LgxBtr9+1Gf99wjPH0CAiCQ6PUgZ7k7J5+phYUJlWBlq89du9DG+GxOrlLKzEQ/tBWAyNdeuGMJ/zbeY4xZMYW4DkpePpdSuO+sObh0d/o0SOGDDyA1E0HCtBWdm5OD6TjPX09UmoSbNjUt5x9/wNuHq4rk/4WEQFK7dAmeP/v2CR9fHuHYsSNI2xq8vSGFh4VBZVNcjFnG229DBRQSgsHt5k00VEkCyb/wAtQ+9nSL1qbnSnDrFsh9/XoYsA4dQsfmUbaZmfBOyspCR9m2DW6on39u+Xr5+ZihbdgA4ti0SRhKp0+3Xg6NBvfkWTkTE4Xb4ahRSEgmN67zBF08dJ+rB3180GZWrULm0SNHXDcbqQyw9V6LivCeIiLsX+fzz62/w/JGURGImQdxye1CN25gRnfoEOxMqanYsrPhtBEYKFJmXLqEPiyPA+JpLYhEYj6etI7HDOTlQUWVmYk2n5dXPvn3q4RKx5Epp7XgBk74+/aBvBcswIs9dAid39Z9ue7w3nvR2blR0xp++AEGs1q1QBDvvy+mtf37wz3w+edNDaw5OdBXN2kCV8C4OPvP2qsXnmv1arh7HTwIabV6dQxOmZkYDN57D9PssvgIWxt0i4tRdh5x2Ls3ylFSgg7AYxYyMqDq4Kql1FTMLBo2xPdz58Syhjt3Qv/P/fGVdBQ/Pwy6bdtC0uTELje+FRfjvqtXw44RGoo6SU8HQfD8LBERmBU1aoStUyd05F9/darqKi240DJoEN4hT7h28iS2EycwozUaYfvp2rVCi1sKXNd/4QJI+dw5zHTPnCkdOzJ3Lma+Fy5YtwneuiU8dmzlH5KnTzYYhE2ppARlKC5Ge2nfHnzE1U+WbAuuhrsJf5wkSUOIaBcRTWaMlQrRkSRpFBGNIiKKU8JiFuBIiL+1UHQ+wufkgBD37IGx1BIJ7toFqZO7zxFhkFi3ThAxV5eY49gxdKCEBKgvpk8H2XNXyb/+EmHdXA3Vqxe8ex54AF4nK1ZAAraFrCx4Dn3yCRpZeDhmBXv2gLS4h1C9evZqzDHk5KDj/PorPEbMXc08PTGYtW8vUkAEBIBI4uJA3uHhwqYiSZhh/forjpVH+cpTGoSEYADj/yUlCS+rN97A4CrXSV++jFnWrl0YpLk0Jw/kMhpxby75x8ejjDExlU+l4A7UqQPj8Y4dkGw5/PwwO2zWDCT/+eew9xiNECSOH0c9JSaiTRcUWNeflxU3bsD+cvgwBqHz57FduIDNXBjQatGfmjXDgO7rCzudp6dw9w0PRxsxGkWuG3Po9Tg3MBDHR0eDN+rWFcn25E4Tp08jAjwkBCrcikKZCF+SpD+IyJLd+kUi+oyIXicidns/l4hGmB/IGJtPRPOJEGlblvIoAXeJNAcn/PBwoQ548MHSx23fLhKiEeEFFhaiEdnDn39Cd80NgpMng/giI/Ef90zg6NkTnSkmRvz23HNooDVrlr5+SQmmop9+CuPuzZsw9Pr7C4MiEQy1ch26s5DPrHhagbAw1GV4OCT2tDR0gIYNYZ9YsgQdIy8P72LDBkhBjGGACA9HJ5ITat26mDrLvS/Cw0Xq2uxskL0k4R48SMjXF1J97dog9fXr4W20b59pPUdEYOAOCRGdXqeDl1SnTmWvpzsVOh3e0ZYtGGxDQlCnISHi/Zw5gzY6ZMj/tXfm4VFUWRt/LyFEgbANSgBBQBaBkV0UQUXIiOACCIygg4g6igs6DooL4AIz6rggzriAOEEcBATEAVF2+URklS1CAAmyhVUICGQjy/3+eOtOVXdXp6vTadJdfX/P0091V1VXV3VVnTr3rPRZWJMI1UM4Lo7C2C4D3Qlnz/KcKWG+Zw9HYunpvmbU2rV5DSYl8SGv2hpWqsT9yMriOZ461dNhbHUcZ2dz5KZ8d3Xq8J5q1YrXr9PILysNG1LTL+t+xyEJfCllspP1hBCTAQSRYB8cwfyJ/gS+eoqr7Mjdu31P6tGjdL4CZlLMwYMUVj16FP+7mzez5kaNGryA9+6lc3j2bAp2pZ0CvMCmTqWJyJuKFX2FvZRMarLGXw8eTIfurFm0248dS83s5pvNCKRQ8E5uuukmJnL17EmzTadOnv+fqr8P8KH1r3/xJlOJKYrjx/mAOHeOPoePP6ZDPT6eo6KCAmrix45Ri0xO5kNNJXhVqcLtKM1u5EgKJFXw7rLLGLJ68iTP9YYN3FZcHB3xgwaxBMI11zgLMXU75cpRKfFHvXoMFVYjsyZNeK2vWsXzkZfHaLVjx+wFflER74f9++1f+/b5+s+qVTOzhuvVoyauWjMq27o3qtZPUpLpMFZ1eM6fp5C/8krG59s52T/7jAI7lGvCX+2pwkIqIIWF/n2MpUU4o3RqSylVx9W+ALaF67ecCPwaNWjm8K7/rlAXSW4uzQeNG/v+Ru/e1CY3b+Z7ZVdu2DBwwaUVKzhVDq7cXGqbd91FzVNdYN26MbkjkD39xAneEFddRZPFuHGc37gxnYdr1jCFOy2NGu6YMcVvL1SGDuUrECNH8iF07bVMiLr+ej6Yliyho/rAAQqLatXMrlSqKuRPP/G/Hj6cv9W6tf/fsd6w1aub4Ypr1tCHAvAGHzSI2mnHjs4S+DSexMV5Fgn0ZvFiCvyZM+lMP3iQ99ovv1DpsWtiHxfHUXBSEq+DypUpkFWzo5dTkQAAIABJREFUnNOnzai2ChX4IGnUyBwtq8J6SkOvWjVyInuWL6eikpBAmaRCMwsLeS1ao9PCQTht+G8IIdqAJp19AB4O1w85cdq2asU0c3+NOVSIXfny9tr6gQPUCufOpdNPaaXx8WZIpzfWC0xpnP36UdArp4/KsJWSETlPP20fZqk4coTlflUZ1UsuMc0T/fuz6cnYsQzZbNeODtlBg/xv70Kh4s1nzeKNmZLCm3TbNrOrk2r+UKmSeY7i4/kQ/MMfODq54orAv7VvH51zip9/Zh4BwN8cNow+kOuvL7XDi1rCnfmqRlwTJnCqGqsnJJhtHVWteIUqhbxnDwV2/foU6g0b0ud04gTPY8eO1IiLu18iiQEDeGxr13KUrRoXqZGGytkJJ2ET+FLKweHatjcXXxx4HSXwmzXzLaC2cKFZCvbuu301gYwMCvsxY+g4BTh8Vdm0KmMW8K9F5Obyif7FF4y0UX4CgNsoV47arz9ycym8X33V05ZduTJHBDffbDYwT0mhSWLNmvBqNcEIiz59+HCrW5eO0muv9a1DrqhTh+au3r1pGnJSiXH/fjrNp0/3rVXz+98zKqprV95gkaDpRQLh+h/y8iisd+2iQ7VnT2r1qammE7R2bTPxTmWgqvdJSVw/Kal0E9PKmiefNEO9ywpXhGVaE5z80awZp02aeAr8zZs9Ba13klFqKl+JiZ5ZjNYelnZISc1SodLya9Qwqzrefz8TtP71L9PkY8e8eSyotm8fheBbb/HmUHXWy5XjEPmll8xGJs8+GxmCraiIjtlFizi832Zj2KtenUP1KVMYuXPppWZDdLtj+O03bisjg//JkiVmI3iADvRnnqGAB2gyGjAgHEenyc6mv2vXLtZfWriQzlTrA7d2bd5/Dz9M/063boGVtJI6eDXF4wqBH6h2u+p3Cfja5q3C3k64vP8+pw0bBlfXZPVqFmRSFBRQW8nMpKnFWnfGXwgnQOHWvz+F4rJlTI5SdOxIbWn8eJos4uIo6J59tvSzY71x+jBJSaGpJj6eGpsqR129Ok1Nd91FE9WWLWZUlOrwtXgxRy5SUlNUlQZXrTKjqipWpGAfNoxmH1XeurCQoXdCeEZVaYKjqIhmRBV7v3cvR2bp6VRoVOKgonVr5nw0a2a+lFlHU/a4QuAHMi1Uq0aBAfjagGvVMkv7WuviABQ6n3xCr723WcGfwFMPHzvnS24utfO//MX5w0NVhxw71lPYS8kIn2eeoX/hzjvpPAu3l98J69bxAfTDD6ZTLj+fwqFHD/oxDh2iL+ODD7hcCWmANvwzZzjq+eorCnnVYOaqq+jnqFePzrnkZPuyC3FxnuWLNcVz9ixzNH78kaOwnTvpGN2/3zcOvUoVCvIbb+RU9T5u2lRHNkU6rhD4gTR81UYO8E15V8K+QgVfe2GvXgzZSkoK3AsUoJCyhpBNmOBZ2757d99yCFJSYwr0ALGWOti1iyGYy5fTgfzJJ2UfL372LB9Kkyf7Jrc1bmw2LfnlFzNzuXVrvv/hB0+TjIqFf/11Dv2TkzkS69XLWYZxWRNBTeQ8ULkSP/1Eob55M6+9xo15XtR+16rFKJemTfmfN25svurV48M0EsyFmuBxhcAPFJZ5ySXF28gB2nit3ZlWrTKjb2rW9B/OqcjPp0lCJZ7ceCM1UYUQ9hl2c+awqcYTT9hvVz1oypenM+y11+i4vfhipoM/+WRwpqbSREpq39OnM/qmoICmsypVaEapX5/aviqrAFCQ3HADNfytW03ndYUKHKG0bMkEnn37zCbkTpzykUCkCcGiIj6AN2yggN++3dPhLwT/9/btmbdx9dX0CXlXY9W4h5gQ+P6iQayMHWsK/OPHzZC9Y8cogIpDStqply41swu/+44OYtWVSdnYraxbx3jyK69kGKUdSuCfO0dNd9UqRhKNH++sYFU4OH/edJjedpvnsqIijnRUn9X27fkgTEqi32LrVn73uutoe7cm6qiuXkuXcn6nTtEj7COF/HyOkH76Cfj0U7OXwk038Vq76iq+WrTg559/NvMSNO7HFQI/0BDaidBQjQlyckxBOmdOYG1HNUE5f55C/ZNPKPC7dWO8bZs2dEh6R/ScP09zT0ICo1j8lf9VJp3nn+d3PvnEt9lKuMnNZbbqRx9xtGSt937RRWaD5nvuYSRGejpNNPfcwxHJxx9z3ZYtWYW0Tx9P04x3SdhI05Qjhfx8Xktr1/KVl0cfzubN1ODXr+cDVdUDEoLhyCtWhN+Jr4kOXCHwA9nwrfZ3u+JpKlQwN9czYUc1MykO5QyuWpVx+lOnMlZ//XqagubP97U7nznDba9dSzu/NY7fSmGh6dSsWJFa8IVIzlCcO0eTyqRJpl3du9NTbi4fAkePmg+1du04ldIU9mvXMqpIC3PnFBRQU1+7lqPOKVM8a8cAZtvLhASaZIYPpzO7TRu+3FSmWRM6rhD4gWzY1lKo1lBJgAk+LVvy5rI2KDh2LLBwslZW3LXLrJXx/ffc1qpVvi3Rfv2VPWBTU4vX1k+cYKXH5cv5+a23wiPsCwpoWmrf3nRanz3LcNS33jIbnwNmk4aEBK7722980G3dap+ToApSvfoqE8GCJVKdn+Hi7Fmas5QGv3GjWSJbCJpi3nyToazly5s+jy5dKOT1w1QTCFcI/EBaTEaG/2WPPMKptUH48OG+phw74aPKJWRmMgqlqIiabk4O489btfKsHqhs/du2UfPv2dN3mytXUtBu3cqHzqBBwIwZFKz5+cCLLzIcc/360IfpubkUxKmp9CHcey+TwD78kMekBL16sMXHs7Txo49S6Jw5w4eFt7DPzKTpZuxYfq5aNbj9iiXBtX8/nflLljBpKTub/3vbtrxWmjfnCLFzZ9949tdeK5t9jlXcoIC4QuA7KZ7WtKln5qsiLo4n8u23+bl5c1NQKewEkIqPL1fOLLr16qsU9vXqsU+uFSlp8pk3j4XO7IT96tWM7gGoQf/1r6yv06oV4+s7djRr/hw9GprALyhgtIwKVx050kxCUyGgylTWpQv35dZbzYeAENT0rcL++HE6k99/34wG6d6dIxUNz9mKFYxsqlSJpT7UNXnZZXy49+9PB2tJWzpqNMXhCoHv5Mn7wAPMQLVjwwa+LrqI2ZqBSjUUFprlk5s04bwVK1jaoHJlRvh4PySWLKGpp3p1asneLFtmll6uU4fa25AhfDA89BAjWipUoPBUhdNCYc8eHrPqfgXwQXXwIIXR6dN0ZH/+ObVNf2RmcjsvvMDEHSGo/Y8axRo2oRDtGtXp04zWWr6cjnlVrwngqOf66zla6tGDCUyxNLKJNtxyblwh8J0kRakoHDsmT6aQc1KkC6DJ49w5DrHj4ujQ/OMfedPm5tpfHN9/z+muXb6/s2YNBXvz5ozJ37mT2nTr1hT0/fpRy58/n+vOnBmaMNyyhfsLeDoBMzK474WFNCsNH+5Zq97K8eN0KFq7+vTuzWSpC9GqLdLIzTWTx4YPp3N740azK9r119Nk1rmz2U+4rPInNLGLKwS+nfBT9e8VVsFkpaCAJY9796awDcS5czT5VK1qCsP77qOJ5/vvGZdu3R9VurhaNT4YvLN5N21iNmPdutQCd+9mBM/ll5uJSX37Mqa6cmU68/zx7rt8WHz4of91/vtfhksmJvI/ycujjXjHDu73kCE0IwWK8VctBGvUYCRRnz6lZ4aIFm3qwAE+hBcsoCavnNSffkoz2JgxNGldc43/B6dGcyFxhcC3s+Fbe1kmJpodj7zZsoUPhj/+0ZnAnzCBQrxNG2aL5uXRXPPaa76abUEBzRsATTHewn73bkbsVKli2nNvuYXhnCpZ7JZbmA8QqOb38eNmGQd/Ar+oiI7AwkLud1YWHyxpaRxRjBtXvPnGyuTJDNcMZ9hfpJl0cnP5UF+6lL4YZX9v2pRmt+xshqEeOhRaQ3iNJly4QuDbCQZriGVOjm9yj2LlSgqtQC0KAT4Y3nyTo4HDhxlGl5PDIbty0lq109dfZ2gm4NmXFqADT1WCXLqU2+vRg5p1QgI/v/EGyzPYabzex+yvObvi+HGWplUPkpYtKaB27qRGX1wtfjsqVQpu/WglN5dRNNOnUyE4d45O7a5dKeRvu80svT1+PKfhatit0YSKawW+FauN39ohCmC9nDvucNZoYf58hiKOGcNMWmucuYq3P32aZpNvvgFefpkavnfqelYWhe+vv1KzP3GCwr52bdp6U1L4W7ff7rsP/swd3hUNrSxdSrOQGvWoZiTx8TRnqaYukUJZm3ROn6aQ/+orvs6c4XVz990U8DfcEHyoqUYTCURJc7Di8Weft8NbmJw9SwHohK1bzXBE6wjir3/ldOdOCu/sbJpICgsZ126lsJA29E2b6HxdsoSOvNq1afZJSeFowbtGTXGcP2/W7ff+rRdf5EgiK4vJVQCFfY0aDMmMNGFfVhQVMZrmnnvMc7FoEcMkFy2imWbSJD6EtbDXRCshafhCiAEAXgbQHEBHKeWPlmXPA3gAQCGAJ6SUNrUiS4dgoh2s5YsB2saTkwN/7+hR2u+rVGHGo0I1JFZJVVaWLPEN8Rw5kvbfd9/lqGLUKM5/+GGabwYMoB04kJZrHdU89ZRZEkLF0O/ezY5ayqTUsSOTtQA6WYcNK3tNOhAXwoa/fz8znlXZgmrV+L/dfTcrfupIGo2bCNWksw3AnQAmWWcKIVoAGAigJYA6AJYJIZpKKQNUvSkZgWrpFEf9+oHj7nNzzVoyqqEHQE1P2Wu/+ILCtUIFatwffsiYfpVpW1TEaJyFCxm299tvXN6iBW3ovXtTE582rXghYyek09I4SujShX1vFy1iKGd2tlm7fP16PlzsqnZGGuF+EOXkcJSTkmKWrkhOps+lTx939VHVRD4nT3K0ffHF4e/OFpLAl1LuAADhe4f2BjBTSpkHYK8QIh1ARwBrQvm9cLBvn+dnO63Sex2FCrU7d45adps2fCA0aECNHTCF18SJdMRWrsziYkOHcv7s2TTfNGhAe3pJwveyskyz1vnzfLBYj6NSJfa6vfXW4LftJtLSaPqaPp12+gYN6GcZMoTRShpNOMnPZ66G6pOdmkoz8eHDXN63L2VAOAmX07YuAGvEeIYxzwchxEMAHgKA+iVsZxTK0P+FF6z74rs8J4cCX/ViHTGCZRisMed16tAX8PnntL/XquW7rcOHaT765hs6fLt1A0aPZrQMwPnBRL5IyQuob19mug4aZI4mrP9H8+ZM5Y/0TlF257C0TDpS0oT10Uc8b/36MfO6a9fA4a4aTUlJS+OIXgn2HTs861K1aME8jVat+GrdOvz7FFDgCyGWAUiyWTRKSjnP39ds5tnevlLKjwB8BAAdOnQo0S0eimBo2rT45dOnU5AmJlLgt2jB+UqgFxZS2ANs6uHNgQOcVqtGJ+0NN9Ckcu+9psbdr599bR071O9mZTG6aNEifh4yxBw1VKtGDRZgZm6khwl6PxydmnTy85nNOncuna3eN8zJkyxXPWkSY+b796f/4pJLSme/NRpFTg4TETdvNl/KZ1anDgV6jx6mYG/WzLNt6YUioMCXUjpwafqQAaCe5fNlAA6XYDuOCEXge5cv9mb1appZzp1jyKRKqClXjr+r0umtbfys+6VMOwMGMIkK4EPivvvoAN64kVp4sKiHyxNP0CyhatYPHUoHZKdODPl0Y4ZnVhYdrS+8YPpUhDAF/ubNjI6aMYP+l+uu42hq0CDn5TM0muKQkvfZ//0fr7cdO0xfYtWqTGAcMYLVeK+4okx31YNwXf7zAUwXQowHnbZNAKwP02+FRHHapJTUkC+6iIJFlVJW31O1ygGzUFheHh2CK1fSFLRsGecvXkyz0Lp1ZpTPnDnBC/vVq833zz/PHADAzAmYMoXa/qRJ0S/svR/kZ87QBj9+PMNfa9fmA/Wf/6TJ7J13gP/8hzegSoZ75JEL2zQm2oi0bOZIpKCA99fevVSoVJ2kggJeZ1270tnfti1fDRpEbgRcqGGZfQH8C8AlAL4WQmyRUvaQUm4XQswCkAagAMBj4YrQAYK7aK1lCwKxcaPZ4Skujvbyr77i53LlzBLAynwCMJkqJ4fRH6+/boZDHjhAx6kyCd13H6N0guHbbxkaCvBBpDz6+fl8eADUKt58M3IvOCd473tmJoX6u+/yv+7Vi9p9585c/uabjG6aNo1NuN99lyazQNFXsU40XyPhpKCACtrKlaxd9eOPpnJ36BCVu4ceYo2kwYOj638MNUrnSwBf+ln2dwB/D2X7TrGrpaO6T3nTqBEFvl0GrDcvvmi+r16dmr5Ko09MpO2+QQMzEefkSTPjVdXUv/NOCvxrrmFst9qmU5u9lWnTGLp17pwZOnjkiGlmio+PfmFvZf9+JrVNnsxj7tOHphmVQKYYO5b/++DBJTOPaWKbjAwqUOq1aROdq3FxVCAefBD48kteW4sWRff95VqLZny8ZwtCgBrhW2/Rvh3IcbdyJUsSKFTYY6tW3G5iIu3B1gbpn37quY1mzcyerj17sozDuHFM7BkwILjjKSqiv+B3vzNHFkeOmFE+N97IQnDRfDF6M24cb7q77mIPAWvCm5UxYy7sfmmiFymZlPjDD7zHly0zO+IlJFDAP/44+wPfcQdNNgD7XVSqFP33lysEvp0mb5dctGqVGVFTHDk5ZucphVWwz5tHW713yznvBtP9+wN/N8Y4QvBCqleP5olgL5w//5n2+/vvZ8JQXh4FYW4uMGsWk6rcYsJISqKDtVEj2ujr1Qv8HY3GDimpCH37Le//H34wa2nVqMHw6OuvpxLYunX0+70C4QqBb0dJQ/2k9NTU//Y3XihHj5rLBwygNmCNxS8qou0YoCYwciRLJt95J8MGN23ihff5587j7adOZfz8gQOeNXZSUmhDTE1lzRfVzGThQmfbjXTi4xkOq9E4RUqG3i5bRv/Zl1/S1GqtItu4Mf0/XbpwtN+sWezlYbhC4NvZ8L3npaUF3o56KCgzDMD6Nust8UUrVnBap45pWgHMevsVKlAob9jAUcbbb1Pgb9sGNGzo3JSzbZtZclmFEr73nhn1o3rRqnaHc+aYYZ8ajRvJz+d9kZZGx2nNmhTu6nXokOf69etToHftymAGPVJ0icC3M+lYSyK/8ILztnuZmfTKA9QAlBavfmPCBF5odeuaMfiA6aytWZMPiL176UBVdfD37KFN2slIIy+PDl5Fixa0N1aq5Hus3bpRs3f7UFQTm2Rl8dpfsIAjb6uSpfjd75ixql6NGoXH1u6GENaYEPhW+3sgrKYEFcWjLp7du3nhjR7NlnZWVFmDkyfN0snDh5v7UbUqtQwnPPOMZyLXwoVmJJC1k9eSJcGHdmo00cC+fQyaGDWKNvfy5ZmlPngwkw3PnuV9vXo1AylizTRTUlwr8K3s2VOy7Xon7Lz3Hi+8Rx7xL/CtdfJTU5loBdDc4uTBM2cOs0SfeoqvihXNCKGzZzlqqFiRzqc2bUp2XBpNpCElgyG++YZmy717OT8piYrXrbeavq+DB+nf+tOfLtw9EO3ROYqYEPgqdDEYRo/2PcmzZzNUq3Zt3/WVwLeSmcmyB4CZcFUcR48yGufqq5m0ZTXTFBSw5s62bRxlaGGviXZOnaJPbOlSBjSsXcv8kh49qMVffTXzV7zLYYwZo0NxS4orBH4gnnjC2XqZmeZ7FfmiyMpi3Lu/7lh2ETJ/+xtNPEBgDSEtzXwwPfOMr01+xAjeGB9/rJ2zmuhESnaFU60jV682gytq1WIY82efaX9UOHGl5UuVH1DYXUDFjQqqVTNr4yhOnaJJRgl8bwGuatpY66qvWsW6NoEoKvLsltW/v+fyHj0Yu//UUyzrq9FECwUFNH+OGMHKtC1aAM8+SwXqhRfokN29m8rU7Nla2IcbVwh87xDM4upKO7HFtWvnuV5REWu49O5tX2pYtREEGLnzzTfm53HjAv/eBx9Q2xk1itE+1t9eupTOWQB4443A29JoyprcXGrwQ4dSc+/alf6vJk14rR88yLyUceOY9NS4sXts5JGOK006dqFbweBdI//kSZY+HTjQfv2nn+Z0xQqGcap44O7dzTo3/sjOZtXLm2/2Ddv8/nvOr1uXSSW6tK8mUsnJoaIzcyYb7uTkMLLs9ts5Ku7RI/L7MsQCrhAh3hp+qALf206vttehg++6UpqRODfcYH7/66/Z1jAQx49z+wMHegr7334zt7dwoVnTw624IcY51jh/nqPPmTMZYXPuHHDppTRj9ukD3HSTNtFEGq4Q+N44qZejsBM03qUPMjOZNWunrSun7P33m7HANWsytRsI3GBdxdV7/+YTT3B7q1b5LxrmFvRwPnooKGDTj5kzmUF+6hQryQ4aRKXlxhvt61hpIgNXCHxvoV0SDd/6kLBesFKynHJior1gyszk/JLa1zdt4rRGDXPenDnMKhwzhkWdNJqyJiuLPqYZMzgqTUykFj9wIJCcrDX5aMEVAt8bJfCvuorOIidY6+VYBX5aGp1QNWv6fuf0aWrw11xjJkcFw6FDbNTRti2HvwA7Nz38MGOQdayx5kIiJe+DqVN5zU+axBo1s2dTo8/OZp2ooUNZ7juYDHZNZOAKge9tw1d9Tvv0Me3ggVizxnxvTdNW8fWJib7fUclWfy9hmxdVauHf/zYfMk8/TYfXtGll0+RYE3tkZLA15KefMk5e8eWXVGoSE9nt7d57GUSgiV5cIfC9US0HvevVF4e1/IJVw//mG27HW/iqlmflyjEaJ1hWrWLK+Jgx1PABlkuYMYMalHekkEZTmmRns+z2vn2sKiklywZPnsyosGHDaI8fMID1mlSHNU10E2pP2wEAXgbQHEBHKeWPxvwGAHYAUPUk10oph4XyW8VRv77nZ5UxG4zAt5ZWVQL/zBkK5gYNfNffsIFTVdQsGIqKgL/8heVan3uO8woL2dAEMNsjajSliZS8bidNopnm7FkqMmPGUHu/4gpzXe9mPhp3EKqGvw3AnQAm2SzbI6W8IBVfvNsVqqYHdmYYO1RXHIUS+MuX02xTq5ZpJvLG6mx1yrRpbJA+bZoZbjllCh86s2aZJZU1mtIgO5s2+A8+4HVXqRI1919+YZTZK6+U9R5qLhQhZdpKKXdIKXcFXjO8eLf2U8LZqYZ/6JDZ9gwwE0SUOcdOqKt1nNrZVSRRdjZTyjt2ZCgbwAiIF19kRI53WQWNpqScOsXOa5ddxpIcubnA++8zMGDKlMB9nTXuI5w2/IZCiM0AzgAYLaX8Plw/5O20VSGWTgW+Co1U1K9PE8vHHzNaxq7WdqdOvt+zwzuU8+23+YD5/HNzu++8w1oic+bomHRNaJw5A3z4Ifs6ZGRQ6A8YADz2GMsY6Osrtgko8IUQywAk2SwaJaWc5+drRwDUl1KeFEK0B/BfIURLKaWPYUQI8RCAhwCgvrcxvoQcOcJpMAJfCA5vy5Xje2Xi6d7dM3JBkZcXfNOFw4dZ9rh/f/bUBOhgfucdpqBfd11w29NoFPv3U9BPnGiaNNu2ZSltncuhUQQU+FLK5GA3KqXMA5BnvN8ohNgDoCmAH23W/QjARwDQoUOHEiXYeydeqQveTuDbaTibNrEFYvXq5jzVkPzRR33LK69cyVetWsHt56uvMlPx9dfNeSNG0MmsHbWakrBhAxvmfPYZ74N+/dgVqnx5xsoH0uh1SYvYIiwmHSHEJQAypZSFQohGAJoA+CUcv1UcTp22a9eyuJOV9HRO7Rof33gjpzVqOL9hCgposrn9djMaYv58ICWFDxrd0ETjlKIi+pf+8Q9GkVWsyKgvFfnlFG3eiT1CctoKIfoKITIAdALwtRBisbHoBgCpQoitAOYAGCalzPS3nVDxJ3SdFBzbsYMOW+8440OH7B2o1s5WTZo438etW4FjxxjfrHj8cU779XO+HU3sUlhIpaFtWyoOBw5wJHrkCH1DwQh7TWwSapTOl1LKy6SUCVLKWlLKHsb8L6SULaWUraWU7aSUX5XO7gaHk9TvGTM4vf12c96hQ7yZunQx56mHyveG63nu3OCSUVav5o2qkrRWrGBd8EGDGC6n0fgjK4vRNc2a0QGbm8us2PR0mhuDyTfRlBw3mL9c2QBFCIacORmyqpBOq8BXZRaUs8u6nfnzKeiDTTE/cYLFp4Tg/nbrxtr5KSm6hILGP+PHM2rs8cdZz2n2bNZ3GjxYXzea4HFlaQUpnVfvO3CAHXmsQn3NGgp1b7u6lKz7nZzsW87YH9YqnH37crpgAaejR+uUdY09UlKLHzGCORvjxzOKS9vdywa3/O+u0PDthloJCc6+u3cv0KqV57w1a4D27T0fGkVFFPT79rHVoVOmTeP03nsZxikli601bGiWVdBorKxZw2Y7991HB/+CBQzjdYvQ0ZQdrhD4djgV+Hl5bKxs/bxxo2/s8vbtwLff8r3V/BOIWbM4VZE5333HErQjR+qWhVbcYB8NlaNHgQcfpO/o8GGa+3bs0BmxmtLDFSLH24YPOBf4gGdlyq1b2brt2mv9r2+Nvy9OUKWkMF7fysSJrJ0/ZIjz/XM7sa65nj/PWPpXXqFDdvhw+nu0oNeUNq7T8F98kdNgbOPW8EqVYduunf26L71kvg8kqKzrAgz/nDePXYJ08wgNACxaRJPi008zv2P7dmDCBC3sNeHBdQL/8ss5dSrwExI8e9Vu2cIwN7uSyL16AS+/7Gy7ubmMu3/6aXPexIk0GT32mLNtaNxLYSETp3r25Puvvwa++iq43A6NJlhcIfCtZpXz5zl1knQFUNhba+Js2cLoHKv2rippOu2eBbA4Wn4+m0eoffz5Zz6Qmjd3vh2N+1i8mNfYc8+xZeCWLVQmNJpw4wqBb0UJ/EAmE/WQsGr3hYVAaqrZgUqxahWnwQj8iRMp2JMtlYgKC7WjNlYpLGRJYiGAW25hG8s5c5jX4TTEV6MJFVcI/MaNzfe5uZwkKj8rAAAIyklEQVT60/CV5q4cvXXrmsvS05nV6K+uTfv2zvYnI4P1ef70J8+RQlGRZ/tEjftJT2eSXfnywP33c97bbzN5ql8/7bDWXFhcoW9ab5qsLE4DaU0nTnBqjbhRDltvgT9zJhO0nCZzzZ3LqXeNnMJCLfBjiZdeAl57zay/NGMG21hqIa8pK1wh8K29OJ0KfFUz31vgx8d7xuUDZq9Zp3zxBdCyJWufWP0LWuDHDp99xpLXyclspFO3rjbnacoeV1yCV17JOiMnTpgCX7Ug9Mfhw5wmWVq7bNlCYe9Uk7fj2DEWWBszxneZFvixwaWXMhqrbVtq9TVrlvUeaTTEFTZ8wKxTs2cPp+vWFb++0vAvvdSct2iRpxO3JMybZzai8EYLfPeSlcVzD1DYjx5NP44W9ppIwhUavhVlH+3Qofj1lHNXVRxUNn3vhujB8sUXdCJfdZXvMi3w3cm33wJ//jPwi9Hi5+BBNg7XaCIN12j4ChVTH2wHqe3bOQ2l5MGpU7z5/UVfFBYG3wdXE7mcOsXaN92787w++ijn6/r0mkjFdeJHae5O2xsqlMD//e9L/tvz57OVob8OVlrDdw9Ll3IU98knLISXmuoZPKDRRCKuE/h5eZwGctp6k5ZGzSwUG/7cuWwz58+clJ2tIzWinexsNiO5+WZeL+vWsUSCro2kiQZcJ/BVpm2w2YvbtzOUsqQx0tnZTJm/8077bRw6xLLInTuXbPuasmfdOkbevP8+G4Zv3Og8GU+jiQRCbWL+phBipxAiVQjxpRCimmXZ80KIdCHELiFEj9B31RlKww+2k1Ramm/8fTB89x1/219NlKlTmWn7wAMl/w1N2ZCfzyqsnTuzJMLy5cA772itXhN9hKrhLwXweyllKwA/A3geAIQQLQAMBNASwC0APhBCXBDrtcpqDCaW/sQJ4PhxavglZfFiCgB/9XbOn+eyRo1K/huaC8+OHWyGM24ccM89tNV361bWe6XRlIyQBL6UcomUssD4uBaACkbrDWCmlDJPSrkXQDqAjqH8llOUScdpAxQpqd0DJdPwVSbtokWsZ17cyGLw4OC3HytEYser//yH/pj9+1nobOpUZ2G7kXgs/oimfS1r3PBfCVlKRyGE+ArA51LKaUKI9wCslVJOM5b9G8BCKeUcm+89BOAh42MzALtC2I2aAE6E8P1oI9aOF9DHHCvoYw6Oy6WUAdvmBIwZEUIsA5Bks2iUlHKesc4oAAUAPlNfs1nf9skipfwIwEeB9sMJQogfpZQBUq7cQ6wdL6CPOVbQxxweAgp8KWVyccuFEEMA3AaguzSHCxkA6llWuwzA4ZLupEaj0WhCJ9QonVsAPAvgDilltmXRfAADhRAJQoiGAJoAWB/Kb2k0Go0mNEJNA3oPQAKApYLB52ullMOklNuFELMApIGmnseklIUh/pYTSsU0FEXE2vEC+phjBX3MYaDUnLYajUajiWxcl2mr0Wg0Gnu0wNdoNJoYwRUCXwhxi1HCIV0I8VxZ708oCCHqCSFWCCF2CCG2CyGeNObXEEIsFULsNqbVjflCCPFP49hThRDtLNsaYqy/24imiliEEHFCiM1CiAXG54ZCiHXGvn8uhKhgzE8wPqcbyxtYtlEm5TxKghCimhBijlGaZIcQolMMnOOnjGt6mxBihhDiIredZyFEihDiuBBim2VeqZ1XIUR7IcRPxnf+KUSQ1b+klFH9AhAHYA+ARgAqANgKoEVZ71cIx1MbQDvjfSJYsqIFgDcAPGfMfw7AP4z3vQAsBHMfrgWwzphfA8AvxrS68b56WR9fMcf9VwDTASwwPs8CMNB4PxHAI8b7RwFMNN4PBJP9YPxHW8EggobGNRFX1sdVzPFOBfCg8b4CgGpuPscA6gLYC+Biy/m9z23nGcANANoB2GaZV2rnFYx27GR8ZyGAnkHtX1n/QaXwB3cCsNjy+XkAz5f1fpXi8c0D8AcwA7m2Ma82gF3G+0kABlnW32UsHwRgkmW+x3qR9ALzNJYD6AZggXExnwBQ3vscA1gMoJPxvryxnvA+79b1Iu0FoIoh/ITXfDef47oADhpCrLxxnnu48TwDaOAl8EvlvBrLdlrme6zn5OUGk466kBQZxryoxxjGtgWwDkAtKeURADCmqhuvv+OPpv9lAoCRAIqMz78DcFqadZqs+/6/4zKW/2asH03H2wjArwCmGGasj4UQleDicyylPATgLQAHABwBz9tGuPs8K0rrvNY13nvPd4wbBL7jMg7RhBCiMoAvAPxFSnmmuFVt5sli5kcUQojbAByXUm60zrZZVQZYFhXHa1AeHPZ/KKVsCyALHOr7I+qP2bBb9wbNMHUAVALQ02ZVN53nQAR7jCEfuxsEvuvKOAgh4kFh/5mUcq4x+5gQoraxvDaA48Z8f8cfLf9LZwB3CCH2AZgJmnUmAKgmhFCJgdZ9/99xGcurAshE9BwvwH3NkFKuMz7PAR8Abj3HAJAMYK+U8lcpZT6AuQCug7vPs6K0zmsGzIrE1vmOcYPA3wCgieHtrwA6eOaX8T6VGMPr/m8AO6SU4y2L5gNQ3vohoG1fzb/X8PhfC+A3Y9i4GMDNQojqhnZ1szEvopBSPi+lvExK2QA8d99KKe8BsAJAf2M17+NV/0N/Y32JKCrnIaU8CuCgEKKZMas7mJXuynNscADAtUKIisY1ro7ZtefZQqmcV2PZWSHEtcZ/eK9lW84oawdHKTlJeoHRLHvAKp5lvk8hHEsXcJiWCmCL8eoF2i+XA9htTGsY6wsA7xvH/hOADpZt3Q/2IkgHMLSsj83BsXeFGaXTCLyR0wHMBpBgzL/I+JxuLG9k+f4o43/YhSCjF8rgWNsA+NE4z/8FozFcfY4BvAJgJ4BtAP4DRtq46jwDmAH6KPJBjfyB0jyvADoY/98esLSNCGb/dGkFjUajiRHcYNLRaDQajQO0wNdoNJoYQQt8jUajiRG0wNdoNJoYQQt8jUajiRG0wNdoNJoYQQt8jUajiRH+HxOCVydxGSK0AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "N = 10000\n", "M = 100\n", "\n", "U = np.random.rand(M, N)\n", "\n", "a = 1.\n", "\n", "########\n", "# Stock in the variable X the N draws from the Cauchy distribution with parameter a,\n", "# exploiting the inverse of the cdf \n", "#\n", "X = a * np.tan( np.pi * (U - 0.5) )\n", "#\n", "########\n", "\n", "integers1toN = np.arange(1,N+1)\n", "\n", "########\n", "# Stock in the variable 'empMean' the sequence of empirical means for n ranging\n", "# from 1 to N\n", "# \n", "empMean = np.cumsum(X, axis=1) / integers1toN\n", "########\n", "\n", "# Plot the sequence of the empirical means\n", "#\n", "plt.plot(integers1toN, empMean.T, color=\"b\")\n", "\n", "#plt.legend(loc=\"best\")\n", "plt.ylim(-20, 20)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 3. Central limit theorem and confidence intervals.\n", "\n", "Let $(X_i)_{i\\ge1}$ be a sequence of i.i.d. random variables with $\\mathbb E[X_1^2]<\\infty$, and denote $\\hat m_n = \\frac1n \\sum_{i=1}^n X_i$ the sequence of their empirical means.\n", "\n", "The central limit theorem states that the sequence of _renormalized errors_\n", "\n", "$$\n", "e_n = \\frac{\\sqrt n}{\\sqrt{Var(X_1)}} \\bigl(\\hat m_n - \\mathbb E[X_1] \\bigr)\n", "$$\n", "\n", "converges in law towards a Gaussian distribution $\\mathcal{N}(0,1)$ as $n \\to \\infty$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 2.1 Convergence in law of the error" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The goal of this first part is to display the closeness of the distribution of the renormalized error $e_n$ with the Gaussian distribution, when $n$ is sufficiently large.\n", " \n", "We consider a sequene of i.i.d. random variables with exponential distribution pf parameter $\\lambda = 2$.\n", "\n", "For fixed $n$, draw a sample of $M$ values $ \\bigl(\\hat m_n^j \\bigr)_{j=1,\\dots,M}$ of the empirical mean.\n", "\n", "Plot the histogram of the corresponding sample of values of the renormalized error\n", "\n", "$$\n", "e^j_n = \\sqrt{\\frac n{Var(X_1)}} \\bigl(\\hat m^j_n - \\mathbb E[X_1] \\bigr),\n", "\\qquad j = 1, \\dots, M\n", "$$\n", "\n", "and compare with the standard Gaussian distribution." ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8FPX5wPHPkxAIch9BJQE2yBkCBAi3cghyCAasoniCF7WK2FoPrK1S1NaKVati1SqtVv0holQuBeQWEEi4IXIEA4QzEE4h5Pr+/pjdJYQcG7LJ7PG8X6997c7Md2aezfHsd78z84wYY1BKKRUcQuwOQCmlVMXRpK+UUkFEk75SSgURTfpKKRVENOkrpVQQ0aSvlFJBRJO+UkoFEU36SikVRDTpK6VUEKlkdwAF1a9f3zgcDrvDUEopv5KUlHTUGBNRUjufS/oOh4PExES7w1BKKb8iIns8aafDO0opFUQ06SulVBDRpK+UUkHE58b0lVIXy87OJi0tjczMTLtDUT4gPDycqKgowsLCLmt9TfpK+bi0tDRq1KiBw+FAROwOR9nIGMOxY8dIS0sjOjr6srahwztK+bjMzEzq1aunCV8hItSrV69M3/o06SvlBzThK5ey/i1o0ldKqSCiSV8pVSoTJkzgtdde8+o2b7zxRk6cOMGJEyd49913vbrt8vLee+/xySef2B1GqemBXKXycYyfc9F06itDbIokuMydOxeA1NRU3n33XR555BGbIyrZww8/bHcIl0V7+kqpEr388su0bNmS/v37s337dvf8lJQUBg0aRKdOnbjuuuv46aefABg9ejTjxo2jR48eNG3alOnTpwNw8OBBevXqRVxcHLGxsSxfvhywyq8cPXqU8ePHk5KSQlxcHE899RQAkyZNonPnzrRr144XXnih0Pg++ugjWrRoQZ8+fXjooYcYO3YsALNmzaJr16506NCB/v37c/jwYeDSbyuxsbGkpqbyyy+/MGTIENq3b09sbCxffPEFAOPHjycmJoZ27drx5JNPXrKN0v4clixZQp8+fbj11ltp1aoVd911F8YYAJKSkujduzedOnVi4MCBHDx4sEy/u4K0p6+UPymvA7rOhFOYpKQkpk6dyvr168nJyaFjx4506tQJgDFjxvDee+/RvHlzVq9ezSOPPMKiRYsAK8H/8MMP/PTTTyQkJHDrrbfy+eefM3DgQJ577jlyc3M5e/bsRft65ZVX2LJlCxs2bABg/vz57Ny5kzVr1mCMISEhgWXLltGrVy/3OgcOHODFF19k3bp11KhRg+uvv5727dsDcO211/Ljjz8iInz44Ye8+uqr/P3vfy/yvX733Xc0bNiQOXOsb3wnT54kIyODGTNm8NNPPyEinDhx4pL1SvtzAFi/fj1bt26lYcOG9OzZkxUrVtC1a1cee+wxvvnmGyIiIvjiiy947rnnmDJlSvG/v1LQpK+UKtby5cu5+eabueKKKwBISEgA4MyZM6xcuZIRI0a4254/f979evjw4YSEhBATE+PuYXfu3Jn777+f7Oxshg8fTlxcXLH7nj9/PvPnz6dDhw7ufe7cufOipL9mzRp69+5N3bp1ARgxYgQ7duwArGscbr/9dg4ePEhWVlaJ57a3bduWJ598kmeeeYahQ4dy3XXXkZOTQ3h4OA8++CBDhgxh6NChF61zOT8HgC5duhAVFQVAXFwcqamp1K5dmy1btnDDDTcAkJuby9VXX11szKXlUdIXkUHAP4BQ4ENjzCtFtLsV+BLobIxJdM57FngAyAXGGWPmeSNwpbyh4Bh+adtX+Jh/MT3y8lTYaYJ5eXnUrl3b3SsvqEqVKu7XrqGLXr16sWzZMubMmcM999zDU089xb333lvkfo0xPPvss/z6178utk1RHnvsMZ544gkSEhJYsmQJEyZMAKBSpUrk5eW527nOe2/RogVJSUnMnTuXZ599lgEDBvD888+zZs0aFi5cyNSpU3nnnXfcvfjL/TkUnB8aGkpOTg7GGNq0acOqVauKfE9lVeKYvoiEApOBwUAMcIeIxBTSrgYwDlidb14MMBJoAwwC3nVuTynlJ3r16sWMGTM4d+4cp0+fZtasWQDUrFmT6OhovvzyS8BKaBs3bix2W3v27KFBgwY89NBDPPDAA6xbt+6i5TVq1OD06dPu6YEDBzJlyhTOnDkDwP79+zly5MhF63Tp0oWlS5dy/PhxcnJy+Oqrr9zLTp48SWRkJAAff/yxe77D4XDve926dfz888+ANVR0xRVXcPfdd/Pkk0+ybt06zpw5w8mTJ7nxxht58803L0nul/NzKErLli1JT093J/3s7Gy2bt16Wdsqiic9/S7ALmPMbgARmQoMA7YVaPci8CrwZL55w4CpxpjzwM8issu5vfL7GFNKeVXHjh25/fbbiYuLo0mTJlx33XXuZZ999hm/+c1veOmll8jOzmbkyJHu8fTCLFmyhEmTJhEWFkb16tUvOeWxXr169OzZk9jYWAYPHsykSZNITk6me/fuAFSvXp1PP/2UBg0auNeJjIzkD3/4A127dqVhw4bExMRQq1YtwDrYOmLECCIjI+nWrZs7ud9yyy188sknxMXF0blzZ1q0aAHA5s2beeqppwgJCSEsLIx//vOfnD59mmHDhpGZmYkxhjfeeOOS91Xan0NRKleuzPTp0xk3bhwnT54kJyeH3/72t7Rp06bU2yqKFPfVCNxDNoOMMQ86p+8BuhpjxuZr0wH4ozHmFhFZAjxpjEkUkXeAH40xnzrbfQR8a4yZXtT+4uPjjd5ERVWU4oZ3OqVt46ucJDh3zppRpQp3hrRnpePCOHRFDO8kJyfTunXrct+PPztz5gzVq1cnJyeHm2++mfvvv5+bb77Z7rDKTWF/EyKSZIyJL2ldT3r6hZ0u4P6kEJEQ4A1gdGnXzbeNMcAYgMaNG3sQklLlp2pWJk8t+4TRSbMo+Of6OV/yWdwg/trnfs5UucKeANUlJkyYwPfff09mZiYDBgxg+PDhdofkszxJ+mlAo3zTUcCBfNM1gFhgifNgz1XATBFJ8GBdAIwxHwAfgNXTL0X8SnlVi/RUPvj6ZRwnDpIjIbzX9RbWN2wFQJvDKTyy6kvu2vAdfVKSePjmP9gcrXLx9hXCgcyTpL8WaC4i0cB+rAOzd7oWGmNOAvVd0wWGd84Bn4vI60BDoDmwxnvhK+U9NTPP8OFXL9L45GGSIxw8eeNv2XpVM/fyBc27MbdlT16b+ybtDu3iX1+/CJPuBOf4su1n9ijlgRLP3jHG5ABjgXlAMjDNGLNVRCY6e/PFrbsVmIZ10Pc74FFjTG7Zw1bKy4zh73PeoPHJw2y+8hqG3/P3ixK+y44IB7+6+zXWRMVw1ZkMuPNOyNU/aeU/PCrDYIyZa4xpYYy5xhjzsnPe88aYmYW07eM6R985/bJzvZbGmG+9F7pS3jNmzdfcsGs1J6tU4zfDn+V8WJUi2+aEVuKxhKc5ekUtWLgQ/vznCoxUqbLR2jsq6HVK28bTS61zuH8/5AnSal9V4jqHa9Rn3E1PWWURXnwR5uk1h8o/aBkGFdTE5PHn79+nksnjva638H3zrh6vu9IRZ/Xyn38exo0jdNgkckPK/9rD0l5FXBI99uC50aNHM3ToUG699VYefPBBnnjiCWJiLrlW1adpT18FtcHbVxJ7OIVD1evyRs87S16hoPHjoWlT2LGDX21ZVHL7AJZb4NhGwemybq+scnJyvLq9Dz/80O8SPmjSV8EsJ4ffL/8UgLd63lHsOH6RwsJg4kQAHl/xOZVzsr0Zoc/49NNP6dKlC3Fxcfz61792J+Tq1avz/PPP07VrV1atWoXD4WDixIlce+21fPnll8WWHHaVGXZtB6wrdvv27cudd95J27ZtL4mjevXqPPfcc7Rv355u3bq5C5jt2bOHfv360a5dO/r168fevXvd+3niiSfo27cvzzzzDBMmTGDUqFEMGDAAh8PB119/zdNPP03btm0ZNGgQ2dnW72/ixIl07tyZ2NhYxowZU2h9nz59+pCYmEhubi6jR48mNjaWtm3buq/YLeq9p6enc8stt9C5c2c6d+7MihUrvPI78pQmfRW8/vtfrslIY0/tq5jW9obL387IkRAbS9SpdO7Y+J334vMRycnJfPHFF6xYsYINGzYQGhrKZ599BsAvv/xCbGwsq1ev5tprrwUgPDycH374gZEjRzJmzBjefvttkpKSeO211zy6OcqaNWt4+eWX2batYKUXa3/dunVj48aN9OrVi3/9618AjB07lnvvvZdNmzZx1113MW7cOPc6O3bs4Pvvv3eXVE5JSWHOnDl888033H333fTt25fNmzdTtWpVd0nlsWPHsnbtWrZs2cK5c+eYPXt2kfFu2LCB/fv3s2XLFjZv3sx9990HUOR7f/zxx/nd737H2rVr+eqrr3jwwQdL/Jl4k47pq+B0/jw4Ky6+ce1d5ISW4V8hNBReegmGD2fsyi+Y1vYGzlUO906cPmDhwoUkJSXRuXNnAM6dO+eufRMaGsott9xyUfvbb78dKLnkcFG6dOlSZAnkypUru0sbd+rUiQULFgCwatUqvv76awDuuecenn76afc6I0aMIDT0wrGWwYMHExYWRtu2bcnNzWXQoEGAVVY5NTUVgMWLF/Pqq69y9uxZMjIyaNOmDTfddFOhMTVt2pTdu3fz2GOPMWTIEAYMGFDse//+++8v+kA7deoUp0+fpkaNGiX+bLxBk74KTh9+CHv3sr1+Y2a27lVy+5IkJLDh6hbEHdzB6HWz+Ge3ESWv4yeMMYwaNYq//vWvlywLDw+/KKECVKtWDSi+5HD+0sbGGLKysi5ZvzBhYWHuMs+ucsSFyV8KuuD2XCWNXUXVXG1DQkLIyckhMzOTRx55hMTERBo1asSECRPcpZcLU6dOHTZu3Mi8efOYPHky06ZN48033yzyvefl5bFq1SqqVq1a5DbLkw7vqOBjDPzjHwC82fNO8rxxxo0Ir197FwD3Js0mNC9wLtjq168f06dPd5c0zsjIYM+ePSWuV1zJYYfDQVJSEgDffPONeyz9cvXo0YOpU6cCVsVL11DT5XAl+Pr163PmzJmLjj0U5ujRo+Tl5XHLLbe47+BV3HsfMGAA77zzjnv9ourwlxft6avgs3gx7NwJkZHMb9Hda5tdFt2RlLqRXJOxn+tT1gLFXrB+2Sr6FMuYmBheeuklBgwYQF5eHmFhYUyePJkmTZqUuG5RJYcfeughhg0bRpcuXejXr1+xvXtPvPXWW9x///1MmjSJiIgI/v3vf1/2tmrXrs1DDz1E27ZtcTgc7mGtouzfv5/77rvP/c3F9Y2oqPf+1ltv8eijj9KuXTtycnLo1asX77333mXHW1olllauaFpaWZW7226DL7+ECRNwnCuxEm2x8idgx/g5PLBmBn9a/BFLojvRZ7d3/o61tLIqqCyllXV4RwWXQ4dgxgzr4KsXzppwjJ/jfgB81bYf50PD6PXzOti9u8zbV8rbNOmr4DJlCuTkwE03gfM2et50ompNZre6lhAM7971zCUfCkrZTZO+Ch65ufDBB9brhx8ut918FncjALdtWkBYrncu1vK1YVhln7L+LWjSV8Hju+9gzx6IjoYbynAxVgnWRbYiOcJB/bMnGbR9ZZm3Fx4ezrFjxzTxK4wxHDt2jPDwy78ORM/eUcHjo48A+FvjXvzzD+VY5VuEz+IG89KCf3LbpgXMiuldps1FRUWRlpZGenq6lwJU/iw8PJyoqKjLXl+TvgoOp07B3LnkIXzdpm+5725W6168sPADeuzdRL1fTnCsWu3L3lZYWFiRV6gqVVoeDe+IyCAR2S4iu0RkfCHLHxaRzSKyQUR+EJEY53yHiJxzzt8gIhV3MqpS+X3zDZw/z9pGbThco37J7cvoZNUaLIvuSKjJY/D2ii2opVRxSkz6IhIKTAYGAzHAHa6kns/nxpi2xpg44FXg9XzLUowxcc5H+R09U6o4X3wBwKxW11XYLmc793XTT8srbJ9KlcSTnn4XYJcxZrcxJguYCgzL38AYcyrfZDVAjzgp33H8OMyfDyEhfNeyR4XtdkHzbpwPDaPzvq00OH2swvarVHE8SfqRwL5802nOeRcRkUdFJAWrpz8u36JoEVkvIktFpNBuloiMEZFEEUnUg1XK62bMgOxs6NuXo9XqVNhuz1S5giVNOxGCYcj2Hypsv0oVx5OkL4XMu6Qnb4yZbIy5BngG+KNz9kGgsTGmA/AE8LmI1Cxk3Q+MMfHGmPiIiAjPo1fKE85CXIwcWeG7nuWs4HlT8rIK37dShfEk6acBjfJNRwEHimk/FRgOYIw5b4w55nydBKQALS4vVKUuQ3o6LFoElSrBr35V4btfeE0XzoZVoeOB7dY1AkrZzJOkvxZoLiLRIlIZGAnMzN9ARJrnmxwC7HTOj3AeCEZEmgLNAS1IoirOV19ZV+LecAPUrVvhuz9XOZxF13SxJqZNq/D9K1VQiUnfGJMDjAXmAcnANGPMVhGZKCKu2rFjRWSriGzAGsYZ5ZzfC9gkIhuB6cDDxpgMr78LpYriqoXuvJuTHdxnDJVQl12piqCllVXgOnkS6teHvDw4cgTq1bOl8FnVrEw2vHUHVXKz4eBBuOqqCo9BBT4trazUvHlWRc2ePaFePdvCOFc5nBVN2lsTc7TaprKXJn0VuGbNsp6LuKF1RVrYzDmuP3u2vYGooKdJXwWmnByYO9d67QtJ33Uwd/58KOYm20qVN036KjCtWgUZGdCsGbRsaXc0HKpZHzp0gLNnrXv0KmUTTfoqMOUf2pHCri+0gesbhys2pWygSV8FJtfYuQ8M7bgNHWo9z54NPnbWnAoemvRV4ElJgeRkqFULrr3W7mgu6NTJOl1z3z7YtMnuaFSQ0qSvAo9r+GTwYAgLszeW/EJCLvT2dYhH2USTvgo8rnPhfWlox8UVk566qWyit0tUgeXsWVi2zDp4O2CA3dFcxDF+DldkZbEhpBKha9YSmpFhSz0gFdy0p68Cy9KlkJVljZ/XL//bIpbW2cpVSYpqTajJg4UL7Q5HBSFN+iqwzJ9vPQ8caG8cxVgW3dF64YpVqQqkwzsqsLgSqY8N7eS3zNGBZ5Z+TNq0b7i27jD3dQSprwyxOTIVDLSnrwLHvn2wbRtUrw7dutkdTZG2XdmUo1fUIupUOtdkpNkdjgoymvRV4FiwwHq+/nqoXNneWIphJIQfHHEAXPfzepujUcFGh3eUX8tfH/+db/7DULhoaMeO+vmeWBbdkeHbltLr53X8Jz6h5BWU8hKPevoiMkhEtovILhEZX8jyh0Vks4hsEJEfRCQm37JnnettFxHfPbqm/FpIXi49UzdYEz48nu+y3NEBgG77NlM5J9vmaFQwKTHpO+9xOxkYDMQAd+RP6k6fG2PaGmPigFeB153rxmDdU7cNMAh413XPXKW8qe2hXdTJPA3R0VZlTR+XXr0uyREOrsg+T6f92+wORwURT4Z3ugC7jDG7AURkKjAMcP+lGmNO5WtfDXBVkxoGTDXGnAd+FpFdzu2t8kLsSrldl2qNjX9WqyXPPTvX5mg8syy6I63TU+n183pWue6spVQ582R4JxLYl286zTnvIiLyqIikYPX0x5Vy3TEikigiienp6Z7GrpSb64DoMkdHmyPxnGuIx/WBpVRF8CTpF1aM/JK6sMaYycaYa4BngD+Wct0PjDHxxpj4iIgID0JS6oKqWZl0OLCdXAlhVZN2dofjsbVRMZwPDSPm8G5qnztV8gpKeYEnST8NaJRvOgo4UEz7qcDwy1xXqVLrnLaVynk5bLnyGk6FV7c7HI+dD6tCUmRrQjB027vZ7nBUkPAk6a8FmotItIhUxjowOzN/AxFpnm9yCLDT+XomMFJEqohINNAcWFP2sJW6oMdeqzb9Sj8cF1/hjLnnno02R6KCRYkHco0xOSIyFpgHhAJTjDFbRWQikGiMmQmMFZH+QDZwHBjlXHeriEzDOuibAzxqjMktp/eiglQPZ8Jc6UdDOy6rmrSD5dBjj95URVUMjy7OMsbMBeYWmPd8vtePF7Puy8DLlxugUsWpmXmG2EMpZIVUYm1UwTOJfd/Gq1twunJVqxzD/v0Qecl5Dkp5lZZhUH6t+95NhGBYH9mKzLBwu8MptdyQUNY0irUmFi2yNxgVFDTpK7/W3TksssIPx/Nd3McitL6+qgCa9JVf6+nH4/ku7tgXLgRzyRnNSnmVJn3lvw4epPmxffwSFs7Gq1vYHc1l+ynCwbGqNSEtDXbtsjscFeA06Sv/5RwDXxvVhuzQMJuDuXxGQljVOF9vX6lypElf+S9n0vfn8XyXlQ7ne9CDuaqcadJX/suZIP15PN9lpaunv3gx5OXZG4wKaJr0lX9KTYXUVE5WqUZyg2i7oymz1DoNrXP0jx61bvmoVDnRpK/805IlAKxpFEteSADcokEE+vSxXi9ebGsoKrBp0ld+xTF+Do7xc5j+908BLhwADQR9+1rPzg80pcqDJn3lf4yhm7PI2o+N29ocjPf0WmOdo5/x7fdEPzPL5mhUoNKkr/xO1MnDRJ1K50R4dZIbOOwOx2v21r6K/TUiqHvuFC3T99gdjgpQmvSV3+nu7OWvbhSLkQD6Exbhx8ZWHR6tr6/KSwD9x6hg4UqIATWe7/Sj8z25PtiU8jZN+sq/GEN3Z9IPpPF8l1XO99R13xY9X1+VC036yq80PnGIhqePklG1Jtsjmtgdjtel1bqStJoR1M48A5t1iEd5n0dJX0QGich2EdklIuMLWf6EiGwTkU0islBEmuRblisiG5yPmQXXVao0AnY830XEPcSj5+ur8lDif42IhAKTgcFADHCHiBS8RdF6IN4Y0w6YDryab9k5Y0yc85HgpbhVkOoWwEM7Lu73pufrq3LgSVepC7DLGLPbGJMFTAWG5W9gjFlsjDnrnPwRiPJumErhPD8/8JO++wD10qWQq7eUVt7lSdKPBPblm05zzivKA8C3+abDRSRRRH4UkeGFrSAiY5xtEtPT0z0ISQWllBSuPnOMY1VrsqN+Y7ujKTf7azVgX60r4cQJHddXXudJ0pdC5hV6ex8RuRuIByblm93YGBMP3Am8KSLXXLIxYz4wxsQbY+IjIiI8CEkFJedwR8CO5+fzYyMd4lHlw5P/nDSgUb7pKOBAwUYi0h94Dkgwxpx3zTfGHHA+7waWAB3KEK8KZs4EGMhDOy46rq/KiydJfy3QXESiRaQyMBK46CwcEekAvI+V8I/km19HRKo4X9cHegJaN1aVnjFBlfRXO6/MZdkyPV9feVWJSd8YkwOMBeYBycA0Y8xWEZkoIq6zcSYB1YEvC5ya2RpIFJGNwGLgFWOMJn1VeikpsH8/x6rWZGcAj+e7pNW6Epo0gePHYZNenau8p5InjYwxc4G5BeY9n+91/yLWWwkEfrdMlb8gGs9369MHPv7Yeu9xcXZHowJEkPz3KL+3dCkQHEM7bq6bqui4vvIiTfrK9wXZeL6bK+nruL7yIk36yvft3g1paVC/flCM57s43ttKWs0GcPw4Nz4w2e5wVIDQpK98n2t4o3fv4BnPd3J9s9H6+spbgus/SPknV9J3DXcEEXfS36dJX3mHJn3l2/KN59O7t62h2MGV9Lvu3ax1eJRXaNJXvi0lxT2eT5s2dkdT4az6+g2odf4XPV9feYUmfeXbXDXle/eGkOD8c9WSDMqbgvO/SPkPV6Lr29fWMOy0Sm+qorxIk77yXcZcSHRBnPTdPf1ly3RcX5WZJn3lu3bsgIMHoUEDaN3a7mhss79WA/bWuhJOnoQNG+wOR/k5TfrKd7l6+X36gBR2W4fgoUM8yls06SvfpeP5bquaaNJX3qFJX/mm/OfnB+FFWQW576S1fDnk5NgbjPJrmvSVz3GMn0P/h96Dw4c5Uq0OtGxpd0i2O1SzPjRrBqdPw7p1doej/JgmfeWTuu+1LkRa1bhd0I/nu7mGuXSIR5WBR0lfRAaJyHYR2SUi4wtZ/oSIbBORTSKyUESa5Fs2SkR2Oh+jvBm8ClyuAmM/Nm6LY/wc9yOoaX195QUlJn0RCQUmA4OBGOAOEYkp0Gw9EG+MaQdMB151rlsXeAHoCnQBXhCROt4LXwUiMXnupL8qmOrnl8TV01++HLKz7Y1F+S1PevpdgF3GmN3GmCxgKjAsfwNjzGJjzFnn5I9AlPP1QGCBMSbDGHMcWAAM8k7oKlC1Sk+l3rlTHKhRn9Q6De0Ox3dcfTW0agW//AJr1tgdjfJTniT9SGBfvuk057yiPAB8W5p1RWSMiCSKSGJ6eroHIalA1iN1IwArm7TX8fyC+vWznhcutDcO5bc8SfqF/deZQhuK3A3EA5NKs64x5gNjTLwxJj4iIsKDkFQg6+E8iLvSdW66usCV9BctsjcO5bc8SfppQKN801HAgYKNRKQ/8ByQYIw5X5p1lXLLzqbrvi0ArGjS3uZgfFDv3ta3n1Wr4OzZktsrVYAnSX8t0FxEokWkMjASmJm/gYh0AN7HSvhH8i2aBwwQkTrOA7gDnPOUKtzatVTPOkdK3SgO16hvdzS+p25d6NgRsrLghx/sjkb5oRKTvjEmBxiLlayTgWnGmK0iMlFEEpzNJgHVgS9FZIOIzHSumwG8iPXBsRaY6JynVOGcwxY6tFMMHeJRZVDJk0bGmLnA3ALzns/3un8x604BplxugCrIOA9Q6tBOMfr1g1df1YO56rJ4lPSVqhBnz8LKleQhF2rIKzfXxWlVszLZGFKJyklJcPw41NFLX5TntAyD8h0rV0JWFtuubMqJqjXtjsZnnasczvrIVlZRuqVL7Q5H+RlN+sp36NCOx9w/Ix3iUaWkSV/5DmcCW6lJv0TuA92a9FUpadJXvuH4cUhKgrAw1kYVLO2kCtpwdUuoVg2Sk2H/frvDUX5Ek77yDYsWQV4e9OjB2cpV7Y7G5+WEVrpQgG3BAnuDUX5Fk77yDfPnW88DBtgbhz9x/aw06atS0KSv7GcMzHNeqK1J33P5k35enr2xKL+hSV/Zb9cu2LMH6tWDDh3sjsZ/tGgBjRtDejps3Gh3NMpPaNJX9nMN7fTvD6Gh9sbiT0Qu9PZdP0OlSqCFhQfOAAAUyklEQVRJX9lPx/MvnyZ9VUpahkHZKzv7wo2+b7jB3lj8jGP8HGqfy2EdQvbSZVT55RfrNE6liqE9fWWv1avh9Gl21muEY/Imvfl5KZ2oWpNNVzejSm4OLFtmdzjKD2jSV/ZyDkssd+gB3Mu13NHReqFDPMoDmvSVvZyJalm0Jv3Ltdz1s9OkrzzgUdIXkUEisl1EdonI+EKW9xKRdSKSIyK3FliW67yxivvmKkoBcOwYrF3L+dBKrG6kpZQv1/qGLTlTuSps2wZ799odjvJxJSZ9EQkFJgODgRjgDhEpWBxlLzAa+LyQTZwzxsQ5HwmFLFfB6rvvIC+P1Y3acq5yuN3R+K3s0DB+cMRZE99+a28wyud50tPvAuwyxuw2xmQBU4Fh+RsYY1KNMZsAvSxQeW6udTO2JU3jbQ7E/y12/Qzn6IFwVTxPTtmMBPblm04DupZiH+EikgjkAK8YY/5XinVVAHKMn0NIXi5JX8+iDrD4Gk36ZbWkaSfrxcKFkJkJ4frNSRXOk56+FDLPlGIfjY0x8cCdwJsics0lOxAZIyKJIpKYnp5eik0rfxV3YAd1Mk/zc52r+blupN3h+L3DNepDXJx1y0m9m5YqhidJPw1olG86Cjjg6Q6MMQecz7uBJcAlp2kYYz4wxsQbY+IjIiI83bTyY9fvXgvo0I5X3Xij9ewcNlOqMJ4k/bVAcxGJFpHKwEjAo7NwRKSOiFRxvq4P9AS2XW6wKnD0TUkEYNE1nW2OJIAMGWI9z5ljVS5VqhAlJn1jTA4wFpgHJAPTjDFbRWSiiCQAiEhnEUkDRgDvi8hW5+qtgUQR2QgsxhrT16Qf5K48fZQ2R3ZzNqwKaxrF2h1O4OjaFerWhZQU2LnT7miUj/Ko9o4xZi4wt8C85/O9Xos17FNwvZWAnoCtLuLq5a9oEsf5SpVtjiaAhIbCoEHw+edWb79FC7sjUj5Ir8hVFa7vbivp61k75UDH9VUJNOmripWZSc891g0/FutBXO8bNMiqs790KZw6ZXc0ygdp0lcVa9EiqmedY1uDaA7W1DO1vK5ePejZ0ypZ/d13dkejfJAmfVWxZswAYF7z7jYHEsCGD7eenT9rpfLTpK8qTm4uzLTO9p3XQpN+ubn5Zut5zhw4f97eWJTP0aSvKs6qVXDkCHtrXclPEQ67owlcTZtCu3Zw+vSFu5Ip5aRJX1Wc/1lll+a16G4dbFTlR4d4VBE06auKYYw7Ac1v3s3mYIKAa4jnm28gT4vfqgv0xuiqYmzZArt3Q0QESZGt7Y4mIF10f2FjSHU4IDUVfvwRevSwKyzlY7SnryqGa5ghIYG8kFB7YwkGIjrEowqlSV+VO8f4OWx59xMA7jt5SbUOVV5cQzwzZmgBNuWmSV+Vu0YnDhF7OIUzlauyskl7u8MJHj17QkSEVYBt0ya7o1E+QpO+Knc3JS8DYEGzrlpgrSKFhsItt1ivp061NxblMzTpq3KXsM26k9Os1r1sjiQIjRxpPU+dqkM8CtCkr8rbli20OrqHE+HVWR59yU3TVHm77jpo2NA6i2f1arujUT5Ak74qX85hhW9b9CA7NMzmYIJQSAjcfrv1+v/+z95YlE/wKOmLyCAR2S4iu0RkfCHLe4nIOhHJEZFbCywbJSI7nY9R3gpc+QFj3El/Zkxvm4MJYnfcYT1Pm2bVP1JBrcSkLyKhwGRgMBAD3CEiMQWa7QVGA58XWLcu8ALQFegCvCAidcoetvILiYmQksKRanVYrbdFtE98PFxzDRw6ZNXZV0HNk55+F2CXMWa3MSYLmAoMy9/AGJNqjNkEFLzeeyCwwBiTYYw5DiwABnkhbuXDHOPn4Bg/h389+hcA5rS6Vi/IspPIxQd0VVDzJOlHAvvyTac553miLOsqPyYmj6HJywE9a8cnuJL+9OmQlWVvLMpWntTeKawcoqfnfnm0roiMAcYANG7c2MNNK1/WY88mrj5zjH21rmRdw1Z2hxOULqrFA3wb4aB1eiq/Hvln3v/6ZZuiUnbzpKefBjTKNx0FHPBw+x6ta4z5wBgTb4yJj4jQW+gFgts3zQdgemw/LaPsI6bH9gPgtk0LbI5E2cmTpL8WaC4i0SJSGRgJzPRw+/OAASJSx3kAd4Bzngpgtc6dZuCOVeQhTG/b3+5wlNOM2OvJDgmlz+4kOOBpv00FmhKTvjEmBxiLlayTgWnGmK0iMlFEEgBEpLOIpAEjgPdFZKtz3QzgRawPjrXAROc8FcCGbVtCldxsfnDEsb9WA7vDUU4ZV9Ti+2ZdCTV58PHHdoejbCLGxy7Njo+PN4mJiXaHocpgy1XNiD2cwtiEp5mtB3F9Sp+Utfxn+p+hWTPYsUOH3gKIiCQZY+JLaqc3UVFllv+AYZvDKcw5nMKJ8Oos0Dtk+Zzl0R05VL0uV+3aBcuXQy/9UA42WoZBedUI50HCGW36akVNH5QbEnrhOMuUKfYGo2yhSV95TZWcLIZvWwLAtHY32BuMKtKXrqT/5Zdw8qS9wagKp0lfeU3CtiXUzjzDpquakdygqd3hqCLsqdMQ+vaFs2fhP/+xOxxVwTTpK+8whvsTrTN5/9PpJpuDUSV67DHr+e23tQhbkNGkr7yi277NtE5PJb1abWa30oODPi8hARwO61aKc+faHY2qQJr0lVe4evmfxt1IViWtm+/zQkNh7Fjr9Ztv2huLqlCa9FWZNTpxiP47V3M+tBKfdRhsdzjKUw88ANWqwaJFsHmz3dGoCqJJX5XZqKRZhGCY3boXR6vp7RL8Ru3aMHq09fqtt2wNRVUcTfqqbE6fdhfwmtIpweZgVKm5Duh++imkp9sbi6oQmvRV2bzzDjWzzrK6USxbr2pmdzSqtFq2hCFDIDMT3njD7mhUBdCkry7f6dPw2msAvNVjpM3BqMv2xz9az2+/DceO2RuLKndae0ddvsmTISODtZExrGjS3u5oVCkUvMFK6sCBMG8evP46vKw3WAlk2tNXl+fMGXcv/x8979Bqjf7uhRes57ffhgytfh7INOmry/Puu9ZQQPfu/OCIszsaVVbdu8OAAdaQnY7tBzRN+qr0zpyBSZOs1xMmaC8/ULh6+//4h/b2A5hHSV9EBonIdhHZJSLjC1leRUS+cC5fLSIO53yHiJwTkQ3Ox3veDV/Z4e0BD8LRoyQ1bIVjYZbd4SgvcIyfg2PmcZY5Oli9/T//2e6QVDkpMemLSCgwGRgMxAB3iEhMgWYPAMeNMc2AN4C/5VuWYoyJcz4e9lLcyi6pqYxZ8zUAL1//gPbyA8xf+t5ProRYB+m3bbM7HFUOPOnpdwF2GWN2G2OygKnAsAJthgGum25OB/qJaDYISE89RZXcbGbE9GFdZGu7o1Fe9lODaP6v/UDIzWXp0HtwPDP7kjN9lH/zJOlHAvvyTac55xXaxnkj9ZNAPeeyaBFZLyJLReS6Msar7LR0KUyfztmwKvyt92i7o1Hl5PXr7uZUlWr0/nkdfXfr/aoDjSdJv7Aee8G7qRfV5iDQ2BjTAXgC+FxEal6yA5ExIpIoIonpeim4b8rNhd/+FoD3ut7KoZr1bQ5IlZeMK2pZp+ECf1z0IWG52TZHpLzJk6SfBjTKNx0FHCiqjYhUAmoBGcaY88aYYwDGmCQgBWhRcAfGmA+MMfHGmPiIiIjSvwtV/t58EzZsgMaN+aDLzXZHo8rZJx2HkFI3imsy9vPoqml2h6O8yJOkvxZoLiLRIlIZGAnMLNBmJjDK+fpWYJExxohIhPNAMCLSFGgO7PZO6KrCJCfDc89ZrydPJjMs3N54VLnLDg3jDwMfBbCS/rp1NkekvKXEpO8cox8LzAOSgWnGmK0iMlFEXGUVPwLqicgurGEc12mdvYBNIrIR6wDvw8YYPQHYn+TkwKhRcP483HcfDB1qd0Sqgqxu3JZ/d7qJsLxcuPde629A+T0xpuDwvL3i4+NNYqIePPIZf/mL1ctv1Mi60UatWno2RxAJz85k7r/H0fT4ARg/Hv76V7tDUkUQkSRjTHxJ7fSKXFW0pCTriluAKVOgVi1bw1EVLzMsnCdv/B2EhMCrr8KyZXaHpMpIk74q3OHDMHw4ZGfzccchOL4/b121qb38oLMuqjU8/TTk5cGIEbBvX8krKZ+lSV9dotlT37C6c39ISyMxsjUv933Q7pCU3V58Efr1gyNH4Fe/sm66ovySJn11iT8u+pCu+7ZwqHpdfjP8WbIqhdkdkrJbpUrwxRfgcEBiIjz8MPjY8UDlGU366mKvv87odbM5H1qJ3wz/A+nV69odkfIV9erB//4HV1wBH39s9f6V39Gkry549134/e8BeHbQY6yPbGVzQMrntG8Pn3xiHdh94QXr4K7yK3q7RGWZMgUetS7G+eMNv+Hr2H42B6R8ycUH8MNJ/fe/YfRoeOYZCA+HcePsCk2Vkvb0FX8a8BvyHrAO1r54/YN82nGIzREpn3fvvfCe8/YYjz8Of/ubjvH7CU36wSwnBx57jBcXvEcIhr/1HsVHnYfbHZXyA47xc3DsjuRPNzhvkTF+PF+0H0jzp/6np/X6OB3eCVbHj8Ndd8G333I+tBLPDH6c/7Xpa3dUys/8t+NQjlSry5uz/87tmxfQ6ORhxg57xu6wVDG0px+MFiyAtm3h22+hXj3uGvmyJnx12ea17MFtd77CkWp16LF3E/OmPAqzZtkdliqCJv1gcvo0jB0LAwbA/v2sa9iS6279G4lRbeyOTPm5zVc3J+HeN1jVuC0Rv5yAhAR44AHrG6XyKVpwLRjk5MC//mXV0TlyBCpV4tUed/J+11vIDQm1OzoVQMTkcX/iTP604r9WVc66deFPf6L5fgfZoRcu8kt9RU8W8DYtuKasf7qPP7aGch55xEr43bvDmjW82/02TfjK64yEWCcDrFsHffpARgb87ncs+PARRm74jirZWp7ZbtrTD0R79vD63c9xz/o51ldtILX21bzSZzTftegBes96VRGM4fqUtfxh8RSaZaQBkFG1Jp/FDebLtv3ZW+dqd1Pt+Zedpz19TfqBYvdumD3bqo+ycqV7dnKEgynxCfyvTd+Lvl4rVVFC83K5KXkZ9yd+Q7tDu9zzN1zdnFmterGoWRcWv/+QdkbKyKtJX0QGAf8AQoEPjTGvFFheBfgE6AQcA243xqQ6lz0LPADkAuOMMfOK25cmfQ/k5DDw1+/T7uAOOu7/iZ57NtL45GH34nOVqrCwWRc+ixvMqsZt9Z9J+QZj6LQ/mbvXz+WGXaupnnXuwrJGjaB/f+jRA+LjafbZHnJCrTPK9VuAZ7yW9J33uN0B3IB1A/S1wB3GmG352jwCtDPGPCwiI4GbjTG3i0gM8H9AF6Ah8D3QwhiTW9T+NOk7/fILHDxo1S7fuxdSU+Gnn2D7duv53LmLmp8Ir86KJu35rkUPFjbrwtnKVe2JWykPhGdn0jclkRu3r6DHno3UO3fqouWZlSqzu24ku+tGMXREH2ja1PpgaNwYrroKqlfXzkwB3kz63YEJxpiBzulnAYwxf83XZp6zzSoRqQQcAiJw3ivX1TZ/u6L2d9lJPyPDGuIoSVHvN/9812tjSn7k5VmP3NyLHzk51iM7G7KyrMf581Yd8sxMOHvWSuxnzlinUp44YT0yMqwDrgWSekF7al/F5quas+mqZqxq3I6tVzYlTw/MKj8kJo/WR1Lpvmcj7Q/tpO2hnUQfP1j8SuHhEBFhVf6sXdt61KwJ1apZHwhVq1qP8HCoUgUqV7YeYWFWmehKlSA01HqEhFx4dj1Ein/AxR86RX0AefLB5HBA/foe/ayK42nS9+SK3Egg/61y0oCuRbUxxuSIyEmgnnP+jwXWjfRgn6W3cCHcdlu5bNoW4eFw5ZVWz6ZRIybvzubnupGk1I0ipV4Up8Kr2x2hUl5hJIRtVzZl25VN3fNqnP+FpsfSuCYjjeiMA0SeOkLkqXSuPpVOxC8nqJqZaX0LDoC7eD09aBzT2g8AKmYoy5OkX9hHVcHuclFtPFkXERkDjHFOnhGR7R7E5WvqA0e9trXMTNizx3r4Ju++X9+n77eCbanY3dn3fr97y3oA8rcybamJJ408SfppQKN801HAgSLapDmHd2oBGR6uizHmA+ADTwL2VSKS6MlXq0Ch7zew6fsNXJ5cnLUWaC4i0SJSGRgJzCzQZiYwyvn6VmCRsQ4WzARGikgVEYkGmgNrvBO6Ukqp0iqxp+8cox8LzMM6ZXOKMWariEwEEo0xM4GPgP+KyC6sHv5I57pbRWQasA3IAR4t7swdpZRS5cvnLs7yVyIyxjlMFRT0/QY2fb+BS5O+UkoFES24ppRSQUSTfjkQkSdFxIhI2a+48GEiMklEfhKRTSIyQ0Rq2x1TeRCRQSKyXUR2ich4u+MpTyLSSEQWi0iyiGwVkcftjqkiiEioiKwXkdl2x1LeNOl7mYg0wipZsdfuWCrAAiDWGNMOq1THszbH43XOMiSTgcFADHCHs7xIoMoBfm+MaQ10Ax4N8Pfr8jiQbHcQFUGTvve9ATxNIRehBRpjzHxjTI5z8kes6zACTRdglzFmtzEmC5gKDLM5pnJjjDlojFnnfH0aKxGWz1X0PkJEooAhwId2x1IRNOl7kYgkAPuNMRvtjsUG9wPf2h1EOSisDElAJ0EXEXEAHYDV9kZS7t7E6qjl2R1IRfDkilyVj4h8D1xVyKLngD8AAyo2ovJV3Ps1xnzjbPMc1rDAZxUZWwXxqJRIoBGR6sBXwG+NMadKau+vRGQocMQYkyQifeyOpyJo0i8lY0z/wuaLSFsgGtgoVmW9KGCdiHQxxhyqwBC9qqj36yIio4ChQD8TmOf/elRKJJCISBhWwv/MGPO13fGUs55AgojcCIQDNUXkU2PM3TbHVW70PP1yIiKpQLwxJmCLdDlvrvM60NsYk253POXBWUtqB9AP2I9VluROY8xWWwMrJ2L1WD4GMowxv7U7nork7Ok/aYwZancs5UnH9FVZvAPUABaIyAYRec/ugLzNeaDaVYYkGZgWqAnfqSdwD3C983e6wdkLVgFCe/pKKRVEtKevlFJBRJO+UkoFEU36SikVRDTpK6VUENGkr5RSQUSTvlJKBRFN+kopFUQ06SulVBD5f0LiUznug9D3AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "n = 3000 # n draws of the exponential distribution to compute the empirical mean\n", "M = 5000 # number of samples of the empirical mean\n", "\n", "####################################################\n", "# n*M independent draws of the exponential distribution\n", "# with parameter lambda=2\n", "#\n", "lambd = 2\n", "\n", "# Pay attention to the arguments of the function random.exponential \n", "X = np.random.exponential(1/lambd, (M, n))\n", "####################################################\n", "\n", "# True known values of expectation and variance (as a benchmark)\n", "esp = 1/lambd\n", "var = (1/lambd)**2\n", "\n", "####################################################\n", "# Stock in the variable 'empMean_n' the sample of M independent \n", "# values of the empirical mean\n", "\n", "# WANT: M independent values of the empirical mean for fixed n \n", "empMean = np.sum(X, axis=1) / n # SAMPLE OF SIZE M\n", "\n", "# Stock in the variable 'renormError_n' the sample of M independent\n", "# values of the renormalized error \n", "renormError_n = np.sqrt(n / var) * (empMean - esp)\n", "\n", "####################################################\n", "# Plotting\n", "plt.hist( renormError_n , density=\"True\", bins=int(np.sqrt(M)), label=\"erreur normalisee\")\n", "\n", "####################################################\n", "# Standard gaussian density for comparison\n", "x = np.linspace(-5, 5, 100)\n", "\n", "gaussianDensity = np.exp( - x*x / 2) / np.sqrt(2*np.pi)\n", "\n", "####################################################\n", "# Plotting\n", "plt.plot(x, gaussianDensity, color=\"red\", label=\"densite gaussienne\", linewidth=2.0)\n", "\n", "plt.legend(loc=\"best\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 2.2 Confidence intervals\n", "\n", "Contrary to the previous example, in practice one does not know the value $Var(X_1)$ (which cannot be used, then, to build the confidence intervals).\n", "\n", "Yet, it is possible to estimate $Var(X_1)$ with the same random sample used to estimate $\\mathbb E[X]$, using the empirical variance as an estimator:\n", "\n", "$$\n", "\\sigma_n^2 = \\frac1n \\sum_{i=1}^n X_i^2 - (\\hat m_n)^2\n", "$$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $\\blacktriangleright$ Question (a):\n", "\n", "Which result allows to state that the random sequence\n", "$$\n", "\\frac{\\sqrt n}{\\sigma_n} \\left(\\hat m_n - \\mathbb E[X_1]\\right)\n", "$$\n", "converges in law towards $\\mathcal{N}(0,1)$, when $n \\to \\infty$?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $\\blacktriangleright$ Question (b):\n", "\n", "Defining for $\\delta > 0$ the (random) interval\n", "$$\n", "I_n^{\\delta}=\\Bigl[\\hat m_n -\\delta\\frac{\\sigma_n}{\\sqrt n},\\ \\hat m_n +\\delta\\frac{\\sigma_n}{\\sqrt n}\\Bigr],\n", "$$\n", "the previous question entails that\n", "\n", "$$\n", "\\mathbb P\\Bigl(\\mathbb E[X_1] \\in I_n^{\\delta}\\Bigr)\n", "=\n", "\\mathbb P\\biggl(\\frac{\\sqrt n}{\\sigma_n} \\bigl|\\ \\hat m_n - \\mathbb E[X_1]\\bigr| \\le \\delta \\biggr)\n", "\\longrightarrow \\mathbb P \\bigl(|\\mathcal{N}(0,1)|\\le \\delta \\bigr) = 2\\int_0^{\\delta} \\frac{e^{-x^2/2}}{\\sqrt{2\\pi}} dx\n", "$$\n", "as $n \\to \\infty$.\n", "\n", "When $\\delta=1.96$, the last term on the right hand side is approximately equal to $0.95$, which allows to identify $I_n^{1.96}$ as the $95\\%$ _asymptotic_ confidence interval for the true value of the expectation $\\mathbb E[X_1]$.\n", "\n", "In the cell below: give an estimate of the expectation of a random variable uniformly distributed on $[0,1]$ with the $95\\%$ _asymptotic_ confidence interval obtained above." ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "True expectation: 0.5000 \n", "\n", "Empirical mean: 0.4931 \n", "\n", "0.95 Confidence interval: [0.4841, 0.5020] \n", "\n", "Relative error*100: 1.8\n" ] } ], "source": [ "N = 4000\n", "X = np.random.rand(N)\n", "\n", "empMean = np.sum(X) / N\n", "\n", "################################################################\n", "# Stock in the variable 'mean_of_squares' the empirical mean\n", "# of the X_i^2\n", "#\n", "mean_of_squares = np.mean(X*X)\n", "\n", "# Stock in the variable 'emp_variance' the empirical variance,\n", "# and in radius_CI the half lenght of the 95% confidence interval\n", "# for the expectation of X\n", "#\n", "emp_variance = mean_of_squares - empMean*empMean\n", "# \n", "radius_CI = 1.96 * np.sqrt(emp_variance / N )\n", "################################################################\n", "\n", "# Theoretical true value for comparison\n", "esp = 0.5\n", "\n", "print(\"True expectation: %1.4f \\n\" %esp)\n", "print(\"Empirical mean: %1.4f \\n\" %empMean)\n", "\n", "print(\"0.95 Confidence interval: [%1.4f, %1.4f] \\n\" %(empMean - radius_CI, empMean + radius_CI) )\n", "\n", "print(\"Relative error*100: %1.1f\" %(radius_CI/empMean * 100) )" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }