{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "$\\newcommand{\\si}{\\sigma}\n", "\\newcommand{\\al}{\\alpha}\n", "\\newcommand{\\tta}{\\theta}\n", "\\newcommand{\\Tta}{\\Theta}\n", "\\newcommand{\\Si}{\\Sigma}\n", "\\newcommand{\\ld}{\\ldots}\n", "\\newcommand{\\cd}{\\cdots}\n", "\\newcommand{\\cN}{\\mathcal{N}}\n", "\\newcommand{\\R}{\\mathbb{R}}\n", "\\newcommand{\\p}{\\mathbb{P}}\n", "\\newcommand{\\f}{\\frac}\n", "\\newcommand{\\ff}{\\frac{1}}\n", "\\newcommand{\\ds}{\\displaystyle}\n", "\\newcommand{\\bE}{\\mathbf{E}}\n", "\\newcommand{\\E}{\\mathbb{E}}\n", "\\newcommand{\\bF}{\\mathbf{F}}\n", "\\newcommand{\\ii}{\\mathrm{i}}\n", "\\newcommand{\\me}{\\mathrm{e}}\n", "\\newcommand{\\hsi}{\\hat{\\sigma}}\n", "\\newcommand{\\hmu}{\\hat{\\mu}}\n", "\\newcommand{\\ste}{\\, ;\\, }\n", "\\newcommand{\\op}{\\operatorname} \n", "\\newcommand{\\argmax}{\\op{argmax}}\n", "\\newcommand{\\lfl}{\\lfloor}\n", "\\newcommand{\\ri}{\\right}\n", "\\newcommand{\\supp}{\\operatorname{supp}}\n", "$\n", "\n", "\"Open\n", "# TP Introduction MCMC\n", "\n", "On commence par importer les librairies qui nous serons utiles dans le TP." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import scipy.stats as sps\n", "import matplotlib.pyplot as plt\n", "from time import time" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from TP_MCMC import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Metropolis-Hasting versus rejet\n", "\n", "### Rappel sur l'algorithme de rejet\n", "Soit $f$ la densité sous laquelle on cherche à  simuler, appelée *densité cible*. On considère une autre densité $g$ (par rapport à $\\mu$), appelée *densité de proposition*, telle que:\n", "- il est aisé de simuler suivant $g$,\n", "- il existe une constante $M$ telle que $f(x)\\leq M g(x)$ pour tout $x$ (ce qui implique que $M\\ge 1$ et $\\supp (f) \\subset \\supp(g)$). \n", "\n", "On peut alors générer un échantillon suivant l'algorithme suivant:\n", "\n", "**Algorithme de rejet de simulation de $X$ de densité $f$**\n", "1. Générer $X$ suivant la densité $g$.\n", "2. Générer $U$ suivant une loi $\\mathcal{U}([0,M])$\n", "3. Accepter la valeur $X$ si $U\\le \\frac{f(X)}{ g(X)}$, sinon, recommencer.\n", "\n", "### Algorithme de Metropolis Hasting\n", "\n", "On a par ailleurs vu en cours l'algorithme de Metropolis-Hasting.\n", "Soit $f$ une densité de probabilité sur $\\R^d$. Posons $\\supp_{>}(f):= {\\{x\\ste f(x)> 0\\}}$.\n", "\n", "**Algorithme de Metropolis-Hastings de simulation de $X$ de densité $\\approx f$**\n", "\n", "1. Choisir un noyau de transition $Q$ irréductible apériodique sur $\\supp_{>}(f)$.\n", "2. Choisir $x\\in \\supp_{>}(f)$.\n", "3. Répéter un grand nombre de fois :\n", " + tirer, indépendamment, $Y\\sim Q(x,\\cdot)$ et $U\\sim\\mathcal{U}([0,1])$,\n", " + remplacer $x$ par $Y$ si $U <\\rho(x,Y)$ et le laisser invariant sinon, où $$\\rho(x,y):= \\f{f(y)Q(x,y)}{f(x)Q(y,x)}.$$\n", "4. Retourner $X:=x$.\n", "\n", "Remarque : on a vu en cours que $\\rho$ devait plutôt s'écrire $\\min(1,\\f{f(y)Q(x,y)}{f(x)Q(y,x)})$, mais comme dans l'algorithme on va utiliser $\\rho$ uniquement dans l'étape 2 de l'algorithme et que $U$ est toujours inférieur à 1, on peut se contenter de choisir $\\rho$ comme plus haut.\n", "\n", "### Exercice\n", "Le but ici est de générer des échantillons selon la loi sur $\\R$ de densité :\n", "$$\\pi_\\al(x) =\\ff{\\sqrt{2\\pi}} \\me^{-x^2/2}(1 + \\sin(\\al x))\\,$$\n", "où $\\al>0$ (on pourra prendre $\\al=2$)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Question 1** : Tracer la densité $\\pi_\\al$ pour plusieurs valeurs de $\\al$. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGxCAYAAACwbLZkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAACWNklEQVR4nOzdeXxU5fX48c/sM9k3sgAhYV9lESwC4g4uKF+Xqq0L1WK/8qWuVNtSW7X82mLrUmgrilWKtVSpdalWqmIVARFFFpU9LCEBErLvy2z398fNnSQkM5mZzGQm4bxfr+lMZu7c+yTSzMl5znMenaIoCkIIIYQQEaKP9ACEEEIIcWaTYEQIIYQQESXBiBBCCCEiSoIRIYQQQkSUBCNCCCGEiCgJRoQQQggRURKMCCGEECKiJBgRQgghRERJMCKEEEKIiJJgRAjhsXr1anQ6nedmtVrJzMzkoosuYunSpZSUlER6iADodDoee+wxz9d79+7lscceIz8/P2JjEkIET4IRIUQHf/nLX/jss89Yv349zzzzDBMnTuS3v/0to0eP5sMPP4z08Pjss8+48847PV/v3buXX/7ylxKMCNFLGSM9ACFE9Bk3bhxTpkzxfH399dfzwAMPcN5553HdddeRl5dHRkZGxMZ37rnnRuzaQojQk8yIEMIvgwYN4qmnnqK2tpaVK1d6nv/yyy+ZO3cuKSkpWK1WJk2axD/+8Y9279Wmfz7++GP+7//+j7S0NFJTU7nuuus4efJku2M/+ugjLrzwQlJTU7HZbAwaNIjrr7+ehoYGzzFtp2lWr17NDTfcAMBFF13kmWJavXo1/+///T+MRiOFhYUdvp/vf//7pKam0tTUFKofkRAiSBKMCCH8duWVV2IwGNi4cSMAH3/8MTNmzKCqqornnnuOf/3rX0ycOJGbbrqJ1atXd3j/nXfeiclk4u9//zu/+93v2LBhA7feeqvn9fz8fObMmYPZbGbVqlW89957PP7448TGxmK32zsd05w5c/jNb34DwDPPPMNnn33GZ599xpw5c7jrrrswGo3tgieAiooKXn31VebPn4/Vag3RT0cIESyZphFC+C02Npa0tDRPNmPhwoWMHTuWjz76CKNR/XVy2WWXUVZWxs9+9jPmzZuHXt/6N8/ll1/OH/7wB8/XFRUV/PjHP6a4uJjMzEy2b99OU1MTTzzxBBMmTPAcd/PNN3sdU79+/Rg+fDgAY8aM6TCF853vfIc///nPPPLII5jNZgBeeOEFmpubWbhwYTd/IkKIUJDMiBAiIIqiAHDo0CH279/PLbfcAoDT6fTcrrzySoqKijhw4EC7986dO7fd1+PHjwfg2LFjAEycOBGz2cz//u//8tJLL3HkyJFuj/e+++6jpKSE1157DQC3282zzz7LnDlzyM3N7fb5hRDdJ8GIEMJv9fX1lJeX079/f06dOgXAgw8+iMlkanfTMg5lZWXt3p+amtrua4vFAkBjYyMAQ4cO5cMPPyQ9PZ0f/vCHDB06lKFDh7J8+fKgxzxp0iRmzpzJM888A8C///1v8vPzufvuu4M+pxAitGSaRgjht3fffReXy8WFF15IWloaAIsXL+a6667r9PiRI0cGfI2ZM2cyc+ZMXC4XX375JX/84x+5//77ycjI4Dvf+U5Q47733nu54YYb2LFjB3/6058YMWIEs2bNCupcQojQk2BECOGXgoICHnzwQRITE7nrrrs8tRpfffWVp4A0lAwGA1OnTmXUqFGsWbOGHTt2eA1GTs+wnO7aa69l0KBB/OhHP+KTTz7h97//PTqdLuRjFkIER4IRIUQHu3fv9tR/lJSUsGnTJv7yl79gMBh488036devHwArV67kiiuu4LLLLuP2229nwIABVFRUsG/fPnbs2OGp0/DXc889x0cffcScOXMYNGgQTU1NrFq1CoBLL73U6/vGjRsHwPPPP098fDxWq5XBgwd7poUMBgM//OEP+clPfkJsbCy33357ED8VIUS4SDAihOjgjjvuAMBsNpOUlMTo0aP5yU9+wp133ukJREDt6/HFF1/w61//mvvvv5/KykpSU1MZM2YMN954Y8DXnThxIh988AGPPvooxcXFxMXFMW7cON5++21mz57t9X2DBw9m2bJlLF++nAsvvBCXy8Vf/vKXdkHHTTfdxE9+8hNuu+02EhMTAx6bECJ8dIpWGi+EEH3YH//4R+699152797N2LFjIz0cIUQbEowIIfq0nTt3cvToUe666y5mzJjBW2+9FekhCSFOI8GIEKJPy83Npbi4mJkzZ/Lyyy+TmZkZ6SEJIU4jwYgQQgghIiqopmcrVqxg8ODBWK1WJk+ezKZNm3wev2bNGiZMmEBMTAxZWVnccccdlJeXBzVgIYQQQvQtAQcja9eu5f777+fhhx9m586dzJw5kyuuuIKCgoJOj9+8eTPz5s1j/vz57Nmzh9dee41t27Zx5513dnvwQgghhOj9Ap6mmTp1KmeffTbPPvus57nRo0dzzTXXsHTp0g7HP/nkkzz77LMcPnzY89wf//hHfve733W6rbcQQgghziwB9Rmx2+1s376dn/70p+2enz17Nlu2bOn0PdOnT+fhhx9m3bp1XHHFFZSUlPDPf/6TOXPmeL1Oc3Mzzc3Nnq/dbjcVFRWkpqZK10QhhBCil1AUhdraWvr3799uB+/ODvTbiRMnFED59NNP2z3/61//WhkxYoTX97322mtKXFycYjQaFUCZO3euYrfbvR7/6KOPKoDc5CY3uclNbnLrA7fCwkKf8UVQHVhPz04oiuI1Y7F3717uvfdeHnnkES677DKKiop46KGHWLBgAS+++GKn71m8eDGLFi3yfF1dXc2gQYMoLCwkISEhmCELIYQQoofV1NSQnZ1NfHy8z+MCCkbS0tIwGAwUFxe3e76kpISMjIxO37N06VJmzJjBQw89BMD48eOJjY1l5syZ/OpXvyIrK6vDeywWi2fjq7YSEhIkGBFCCCF6ma5KLAJaTWM2m5k8eTLr169v9/z69euZPn16p+9paGjoME9kMBgANaMihBBCiDNbwEt7Fy1axAsvvMCqVavYt28fDzzwAAUFBSxYsABQp1jmzZvnOf7qq6/mjTfe4Nlnn+XIkSN8+umn3HvvvXzrW9+if//+oftOhBBCCNErBVwzctNNN1FeXs6SJUsoKipi3LhxrFu3jpycHACKiora9Ry5/fbbqa2t5U9/+hM/+tGPSEpK4uKLL+a3v/1t6L4LIYQQQvRavaIdfE1NDYmJiVRXV0vNiBBC9HKKouB0OnG5XJEeiugmg8GA0Wj0WhPi7+d3UKtphBBCiGDY7XaKiopoaGiI9FBEiGhbvZjN5qDPIcGIEEKIHuF2uzl69CgGg4H+/ftjNpulkWUvpigKdrud0tJSjh49yvDhw303NvNBghEhhBA9wm6343a7yc7OJiYmJtLDESFgs9kwmUwcO3YMu92O1WoN6jzBhTBCCCFEkIL961lEp1D895R/EUIIIYSIKAlGhBBCCBFREowIIYQQ3ZCfn49Op2PXrl1+v2f16tUkJSWFbUy9jQQjQgghhADUxqU333wzI0eORK/Xc//99/fIdSUYEUL0Kg6Xg1f2v8I/DvwDl1uaZgkRSs3NzfTr14+HH36YCRMm9Nh1JRgRQvQqS79Yym8+/w3/b+v/4487/xjp4YhuUBSFBrszIrdAmo+/9957nHfeeSQlJZGamspVV13F4cOHvR6/YcMGdDod7777LhMmTMBqtTJ16lS++eabDse+//77jB49mri4OC6//HKKioo8r23bto1Zs2aRlpZGYmIiF1xwATt27Ajshxyg3Nxcli9fzrx580hMTAzrtdqSPiNCiF7jVP0p3sh7w/P1y3tf5o5xd5Bo6blfmiJ0Gh0uxjzyfkSuvXfJZcSY/fsIrK+vZ9GiRZx11lnU19fzyCOPcO2117Jr1y6fy1ofeughli9fTmZmJj/72c+YO3cuBw8exGQyAequ9k8++SQvv/wyer2eW2+9lQcffJA1a9YAUFtby/e+9z3+8Ic/APDUU09x5ZVXkpeXR3x8fKfXXLNmDXfddZfP72flypXccsstfn3vPUWCESFEr/Hu0XdxKS7OTj+bOkcdBysP8uGxD7l+xPWRHprow66/vv2/rxdffJH09HT27t3LuHHjvL7v0UcfZdasWQC89NJLDBw4kDfffJMbb7wRAIfDwXPPPcfQoUMBuPvuu1myZInn/RdffHG7861cuZLk5GQ++eQTrrrqqk6vOXfuXKZOnerz+8nIyPD5eiRIMCKE6DU+O/kZALNzZ1NnV4ORT09+KsFIL2UzGdi75LKIXdtfhw8f5he/+AVbt26lrKwMt9sNQEFBgc9gZNq0aZ7HKSkpjBw5kn379nmei4mJ8QQiAFlZWZSUlHi+Likp4ZFHHuGjjz7i1KlTuFwuGhoaKCgo8HrN+Ph4r1mTaCbBiBCiV7C77Owq2QXA1Myp1DpqYRdsP7UdRVFkj5NeSKfT+T1VEklXX3012dnZ/PnPf6Z///643W7GjRuH3W4P+Fxt/51q0zVtX2tby3L77bdTWlrKsmXLyMnJwWKxMG3aNJ/XlWkaIYQIowMVB2hyNZFsSWZo0lAcbgcWg4WKpgqO1RwjNzE30kMUfVB5eTn79u1j5cqVzJw5E4DNmzf79d6tW7cyaNAgACorKzl48CCjRo3y+9qbNm1ixYoVXHnllQAUFhZSVlbm8z0yTSOEEGF0sPIgAKNSRqHT6TAbzIxMGcnXpV+zv3K/BCMiLJKTk0lNTeX5558nKyuLgoICfvrTn/r13iVLlpCamkpGRgYPP/wwaWlpXHPNNX5fe9iwYbz88stMmTKFmpoaHnroIWw2m8/3hGKaRmveVldXR2lpKbt27cJsNjNmzJhundcXWdorhOgVDlQeAGBkykjPcyOSRwBwsOJgRMYk+j69Xs+rr77K9u3bGTduHA888ABPPPGEX+99/PHHue+++5g8eTJFRUW8/fbbmM1mv6+9atUqKisrmTRpErfddhv33nsv6enpwX4rfps0aRKTJk1i+/bt/P3vf2fSpEme7Ey4SGZECNEraJkRLQBp+zivMi8iYxJnhksvvZS9e/e2e65tbUdubm6nfUvOO+88du/e3ek5b7/9dm6//fZ2z11zzTXtzjNp0iS2bdvW7phvf/vbgQ4/YIH0YAkVyYwIIaKeoig+gxHtNSFE7yTBiBAi6lU1V1FrrwUgJyHH8/zgxMEAFNUX0exqjsjYhBDdJ8GIECLqFdYWApAek47VaPU8n2xJJs4Uh4LCidoTkRqeEO1ceOGFKIoiu/IGQIIRIUTUO157HICBcQPbPa/T6ciOzwagoNZ7IyghRHSTYEQIEfW0zIgWeLSlPacdI4TofSQYEUJEveN1LZmR+IEdXvNkRmokMyJEbyXBiBAi6vnKjAxKGNTuGCFE7yPBiBAi6nlqRjrJjAyIGwCoK2qEEL2TBCNCiKjmdDspbSwFoH9s/w6vZ8So+2wU1xdHpFmTEKL7JBgRQkS18sZy3Iobg85AijWlw+sZsWow0uBsoM5R19PDE4L8/Hx0Op1nTxd/rF69Wpb+tiHBiBAiqpU0lACQZkvDoDd0eN1mtJFoSQTU7IgQInhvvPEGs2bNol+/fiQkJDBt2jTef//9sF9XghEhRFTTghFtOqYz2munGk71yJiE6Ks2btzIrFmzWLduHdu3b+eiiy7i6quvZufOnWG9rgQjQoiopgUY/WL6eT3GE4zUSzDSqygK2Osjcwugvui9997jvPPOIykpidTUVK666ioOHz7s9fgNGzag0+l49913mTBhAlarlalTp/LNN990OPb9999n9OjRxMXFcfnll1NU1FqIvW3bNmbNmkVaWhqJiYlccMEF7NixI7CfcYCWLVvGj3/8Y8455xyGDx/Ob37zG4YPH84777wT1uvKrr1CiKimZUbSY7xvnZ4ZmwlAcYNM0/Qqjgb4Tcei5B7xs5NgjvXr0Pr6ehYtWsRZZ51FfX09jzzyCNdeey27du1Cr/f+N/1DDz3E8uXLyczM5Gc/+xlz587l4MGDmEwmABoaGnjyySd5+eWX0ev13HrrrTz44IOsWbMGgNraWr73ve/xhz/8AYCnnnqKK6+8kry8POLj4zu95po1a7jrrrt8fj8rV67klltu8et7d7vd1NbWkpLSsV4rlCQYEUJENX+CEcmMiHC6/vrr23394osvkp6ezt69exk3bpzX9z366KPMmjULgJdeeomBAwfy5ptvcuONNwLgcDh47rnnGDp0KAB33303S5Ys8bz/4osvbne+lStXkpyczCeffMJVV13V6TXnzp3L1KlTfX4/GRnepzxP99RTT1FfX+8Zc7gEFYysWLGCJ554gqKiIsaOHcuyZcuYOXNmp8fefvvtvPTSSx2eHzNmDHv27Anm8kKIM0hJox81I7Gty3tFL2KKUTMUkbq2nw4fPswvfvELtm7dSllZGW63G4CCggKfwci0adM8j1NSUhg5ciT79u3zPBcTE+MJRACysrIoKSnxfF1SUsIjjzzCRx99xKlTp3C5XDQ0NFBQ4L3bcHx8vNesSaBeeeUVHnvsMf71r3+Rnu79j4FQCDgYWbt2Lffffz8rVqxgxowZrFy5kiuuuIK9e/cyaNCgDscvX76cxx9/3PO10+lkwoQJ3HDDDd0buRDijOBPZiTdpr6m9SMRvYRO5/dUSSRdffXVZGdn8+c//5n+/fvjdrsZN24cdrs94HPpdDrPY226pu1rbXvl3H777ZSWlrJs2TJycnKwWCxMmzbN53VDNU2zdu1a5s+fz2uvvcall17q89hQCDgYefrpp5k/fz533nknoBa7vP/++zz77LMsXbq0w/GJiYkkJiZ6vn7rrbeorKzkjjvu6MawhRBnCn+CkVRbKgAVTRU9MiZx5igvL2ffvn2sXLnSMwOwefNmv967detWzx/plZWVHDx4kFGjRvl97U2bNrFixQquvPJKAAoLCykrK/P5nlBM07zyyit8//vf55VXXmHOnDl+j7c7AgpG7HY727dv56c//Wm752fPns2WLVv8OseLL77IpZdeSk5OjtdjmpubaW5u9nxdU1MTyDCFEH1Eg6OBekc9AP1s3lfTaMFIZVMlTrcTo17K4URoJCcnk5qayvPPP09WVhYFBQUdPgO9WbJkCampqWRkZPDwww+TlpbGNddc4/e1hw0bxssvv8yUKVOoqanhoYcewmaz+XxPd6dpXnnlFebNm8fy5cs599xzKS5Wpz5tNlu7xEKoBbS0t6ysDJfL1SGqysjI8AzYl6KiIv7zn/94sireLF261JNRSUxMJDu74+ZYQoi+r7K5EgCLwUKsyXs6P9mSjF6nR0Ghsqmyp4YnzgB6vZ5XX32V7du3M27cOB544AGeeOIJv977+OOPc9999zF58mSKiop4++23MZvNfl971apVVFZWMmnSJG677TbuvffesNdurFy5EqfTyQ9/+EOysrI8t/vuuy+s1w3qz4e2c14AiqJ0eK4zWvvbriLDxYsXs2jRIs/XNTU1EpAIcQaqaFSnXVKsKT5/xxj0BpItyZQ3lVPeVO6zJ4kQgbr00kvZu3dvu+fa1nbk5uZ2ui/Seeedx+7duzs95+23387tt9/e7rlrrrmm3XkmTZrEtm3b2h3z7W9/O9DhB2TDhg1hPb83AQUjaWlpGAyGDlmQkpKSLuegFEVh1apV3HbbbV1GhhaLBYvFEsjQhBB9kFYD0tmeNKdLtaVS3lROWaPvOXUhRPQJaJrGbDYzefJk1q9f3+759evXM336dJ/v/eSTTzh06BDz588PfJRCiDOSFowkW5O7PDbVqtaNlDeWh3VMQojQC3iaZtGiRdx2221MmTKFadOm8fzzz1NQUMCCBQsAdYrlxIkT/PWvf233vhdffJGpU6f6XJMthBBtlTepgYU/mZE0W1q79wgRKRdeeGGn0zbCu4CDkZtuuony8nKWLFlCUVER48aNY926dZ7VMUVFRR0aslRXV/P666+zfPny0IxaCHFG0DIjWtbDF21FjUzTCNH7BFXAunDhQhYuXNjpa6tXr+7wXGJiIg0NDcFcSghxBgukZsSTGZFpGiF6Hdm1VwgRtbRluim2roMRLWCRYESI3keCESFE1Ap0NQ1IzYgQvZEEI0KIqKX1GZHVNEL0bRKMCCGikqIoARWwagFLtb0at+IO69iEEKElwYgQIirV2GtwKk7Av8xIkiUJALfiptZeG86hCdFOfn4+Op2OXbt2+f0erSO5UEkwIoSISlpWJM4Uh8XQdUdms8Hs2b+mqrkqnEMTos/asGEDOp2uw23//v1hva5sbSmEiEqelTR+FK9qkixJ1DvqqWyqJCfB+87gQgjfDhw4QEJCgufrfv3Cu9+TZEaEEFFJ27FXm37xh3asZEZ6B0VRaHA0ROQWSIfU9957j/POO4+kpCRSU1O56qqrOHz4sNfjtezCu+++y4QJE7BarUydOpVvvvmmw7Hvv/8+o0ePJi4ujssvv5yioiLPa9u2bWPWrFmkpaWRmJjIBRdcwI4dOwL7IQcpPT2dzMxMz81gMIT1epIZEUJEpZrmGgASLYl+vyfJmgS0ZlVEdGt0NjL171Mjcu3Pb/6cGFOMX8fW19ezaNEizjrrLOrr63nkkUe49tpr2bVrF3q997/pH3roIZYvX05mZiY/+9nPmDt3LgcPHsRkMgHQ0NDAk08+ycsvv4xer+fWW2/lwQcfZM2aNQDU1tbyve99jz/84Q8APPXUU1x55ZXk5eURHx/f6TXXrFnDXXfd5fP7WblyJbfccovPYyZNmkRTUxNjxozh5z//ORdddJHP47tLghEhRFTSshuBBCPJlpYVNc3V4RiSOENdf/317b5+8cUXSU9PZ+/evT73W3v00UeZNWsWAC+99BIDBw7kzTff5MYbbwTA4XDw3HPPMXToUADuvvtulixZ4nn/xRdf3O58K1euJDk5mU8++YSrrrqq02vOnTuXqVN9B3gZGRleX8vKyuL5559n8uTJNDc38/LLL3PJJZewYcMGzj//fJ/n7Q4JRoQQUUkLKIKZptGmeER0sxltfH7z5xG7tr8OHz7ML37xC7Zu3UpZWRlut7p0vKCgwGcwMm3aNM/jlJQURo4cyb59+zzPxcTEeAIRUAOBkpISz9clJSU88sgjfPTRR5w6dQqXy0VDQ0OH/d/aio+P95o18cfIkSMZOXJku++hsLCQJ598UoIRIcSZp9quBiMJloQujmylLQGWmpHeQafT+T1VEklXX3012dnZ/PnPf6Z///643W7GjRuH3W4P+Fw6nc7zWJuuafta21qW22+/ndLSUpYtW0ZOTg4Wi4Vp06b5vG6opmnaOvfcc/nb3/7m9/HBkGBECBGVtMxIojmAmhEtMyI1IyJEysvL2bdvHytXrmTmzJkAbN682a/3bt26lUGDBgFQWVnJwYMHGTVqlN/X3rRpEytWrODKK68EoLCwkLIy37tSd3eapjM7d+4kKysroPcESoIRIURUCqqAVVbTiBBLTk4mNTWV559/nqysLAoKCvjpT3/q13uXLFlCamoqGRkZPPzww6SlpXHNNdf4fe1hw4bx8ssvM2XKFGpqanjooYew2XxPL3V3mmbZsmXk5uYyduxY7HY7f/vb33j99dd5/fXXgz6nP2RprxAiKmnTNAEVsLadpmmqgc3L4J/fhw9+DmWHwjBK0dfp9XpeffVVtm/fzrhx43jggQd44okn/Hrv448/zn333cfkyZMpKiri7bffxmw2+33tVatWUVlZyaRJk7jtttu49957SU9PD/Zb8YvdbufBBx9k/PjxzJw5k82bN/Puu+9y3XXXhfW6OiWQxdYRUlNTQ2JiItXV1e2asAgh+q5Z/5xFcX0xr8x5hXFp3osE28qrzOO6t68j2RTPxuIqqDnR+qLBAnP/ABO+E54Biy41NTVx9OhRBg8ejNVqjfRwwmbDhg1cdNFFVFZWnhEt3339d/X381syI0KIqBRMzUjrZnk1uGpOQPJguOQRGHIhuJrhrf+DA++FY7hCiG6QYEQIEXXsLjuNzkYgsNU0iSb1WLdOR22/UfCDj2Dmj+C2t+DseaC41YCk9lQ4hi2ECJIEI0KIqFNjV4tXdeiIN/tfjGc6/F/iWnpAVF25FGJa9rXR6eDKJyHzLGisgI+W+DiLEN1z4YUXoijKGTFFEyoSjAghoo42RZNgSUCv8/PXlNsFHz5GoqslGIk9bYM9owXmPK0+3rkGTu0J1XCFEN0kwYgQIuoE032VvA+gdD8JqE2ltOxKO9nfgtFzAQU+Xd79gYqg9IJ1EyIAofjvKcGIECLqePalCaB4lc+fAyAhNhPwsT/NeQ+o97tfh5qTwQ5RBKHtBnGi79D+e57eUTYQ0vRMCBF12k7T+KX0IBzZADo9iWkj4WRp55kRgAFnQ84MOPYpfLkKLv55aAYtumQwGEhKSvLsvxITE9OuPbroXRRFoaGhgZKSEpKSkjAYDEGfS4IRIUTU0QIJvxuefb1WvR8+m4S4zHbn6NSU76vByFevwoU/Ax/bwIvQysxU//u03RBO9G5JSUme/67BkmBECBF1Auoxoiiw+5/q47NuIMFeCLS2k+/UqDlgSYTqQji2GQaHbzdS0Z5OpyMrK4v09HQcDkekhyO6yWQydSsjopFgRAgRdQIqYD2xAyrzwRQDI68g4cCrQBeZEZMNxl4DO15SsyMSjPQ4g8EQkg8x0TdIblIIEXW0fWn8qhnZ84Z6P/JKMMd63uMzGAEYf6N6v/9dcDmDHaoQIgQkGBFCRB3Pahp/akYOvq/ej74agARzSzDia5oGYNA0iEmFpiq1fkQIETESjAghok6tvRZoDSy8qjgC5XmgN8LQi9q9p8vMiN4AI69QH+9/t1vjFUJ0jwQjQoioowUjXbaCz/tQvR80DaxqFkXLpnSZGQEYdZV6v/9dtRBWCBEREowIIaKOJxgxdRWMfKDeD5/lecrvzAiou/maYqDmOJTsDWaoQogQkGBECBFVFEWhzl4HdJEZcTRB/ib18fDZnqe1AtYmVxN2l933xUw2yJmuPj78cdBjFkJ0jwQjQoio0uhsxKmoq1t8BiMnvgRnE8RlQL9RnqfjTHHofO1Pc7qhF6v3hz8KesxCiO4JKhhZsWIFgwcPxmq1MnnyZDZt2uTz+ObmZh5++GFycnKwWCwMHTqUVatWBTVgIUTfVudQsyIGnQGb0eb9wPyWFTA5M6BNS3G9Tu8JYvyqGxmiFr5ybIuabRFC9LiAm56tXbuW+++/nxUrVjBjxgxWrlzJFVdcwd69exk0aFCn77nxxhs5deoUL774IsOGDaOkpASnU9b1CyE6alu86nPfkmOb1fvcGR1eSjAnUGOv8S8zkj4a4jKhrhgKP4chFwQzbCFENwScGXn66aeZP38+d955J6NHj2bZsmVkZ2fz7LPPdnr8e++9xyeffMK6deu49NJLyc3N5Vvf+hbTp0/v9uCFEH2PFozEmeK8H+S0Q+E29XHOeR1e1upGvO7c25ZO51kWzBGpGxEiEgIKRux2O9u3b2f27Nntnp89ezZbtmzp9D1vv/02U6ZM4Xe/+x0DBgxgxIgRPPjggzQ2Nnq9TnNzMzU1Ne1uQogzg1/Lek/uAGcjxKRBv5EdXtb2tPErMwKQ2xLQHPssoLEKIUIjoGmasrIyXC4XGRkZ7Z7PyMiguLi40/ccOXKEzZs3Y7VaefPNNykrK2PhwoVUVFR4rRtZunQpv/zlLwMZmhCij/ArGMlvmaLJmd6uXkTjd0t4zaBp6v3JHWrdiMnq93iFEN0XVAHr6fO4iqJ4ndt1u93odDrWrFnDt771La688kqefvppVq9e7TU7snjxYqqrqz23wsLCYIYphOiF/ApGjmtTNJ1P9/rdEl6TMkRdleOyw4ntfo9VCBEaAQUjaWlpGAyGDlmQkpKSDtkSTVZWFgMGDCAxsXWPidGjR6MoCsePH+/0PRaLhYSEhHY3IcSZodbRRTCiKHD8S/XxgCmdHhJQ4zNQsytadqSg8ylnIUT4BBSMmM1mJk+ezPr169s9v379eq8FqTNmzODkyZPU1dV5njt48CB6vZ6BAwcGMWQhRF/WZQFrVQE0lIHeBJlndXpIwNM00JplKdjq/3uEECER8DTNokWLeOGFF1i1ahX79u3jgQceoKCggAULFgDqFMu8efM8x998882kpqZyxx13sHfvXjZu3MhDDz3E97//fWw2Hz0EhBBnJK37qtdN8rRplMxxXms7Ap6mgdbMSOEX4Hb5/z4hRLcF3Gfkpptuory8nCVLllBUVMS4ceNYt24dOTk5ABQVFVFQUOA5Pi4ujvXr13PPPfcwZcoUUlNTufHGG/nVr34Vuu9CCNFndFkzogUjAyZ7PUfA0zQAGWPBkgDNNXBqN2RN8P+9QohuCTgYAVi4cCELFy7s9LXVq1d3eG7UqFEdpnaEEKIzNQ41gIgze5mm8QQjndeLtH2vVn/iF71BDXCOfKxeQ4IRIXqM7E0jhIgqPjfJcznh5C71sY/MiLbbr3Yuv2nnlBU1QvQoCUaEEFHFM01j6iQYKdmrNjuzJELqMK/n0AIZ7Vx+8wQjOwJ7nxCiWyQYEUJEFZ+ZkZM71fv+E0Hv/deXNk1T76jHrbj9v/iAs9X7kn3QHGAgI4QImgQjQoio4rPPSPE36n3WeJ/n0N6roFDvqPf/4vGZkDAQUKDoK//fJ4ToFglGhBBRw+F20OhUOzP7DEYyfQcjFoMFs94MBFM30pIdkboRIXqMBCNCiKjRNnCINcW2f9HthlN71McZ47o8lzZVE9DyXpAiViEiQIIRIUTU0ApOY02xGPWndR6oygd7LRgskDa8y3NpmZU6R7AraqSIVYieIsGIECJq+GwFr03RpI8Gg6nLc2nnCHiapv9EQAfVhVBXGth7hRBBkWBECBE1fBev7lbvM7ueoml7joCnaSzxkDq05ZpfB/ZeIURQJBgRQkQNLTPS6b40fhavaoKepoHWDfgkGBGiR0gwIoSIGtqUSqet4LVgxI/iVejGNA20BjzaNYUQYSXBiBAiamhTKh2maRoqoOa4+tjPaRrP/jSBdmGF1mCkSDIjQvQECUaEEFHDawHrqZZ6kaQcsCb6dS5PS/hANsvTaE3Vyg+BPYCmaUKIoEgwIoSIGlp9R4fMSMl+9T59jN/nCnqzPIC4dIjLAJTW3iZCiLCRYEQIETU8NSOnZ0ZK96n36aP8PpdnmiaYzAi0qRuRqRohwk2CESFE1GhwNgCddF8tPaDe9/M/GAl6516NtqJG6kaECDsJRoQQUUPLjHQMRlqmaQIJRrozTQOtdSOyokaIsJNgRAgRNeqdarFou2CkrhQaygEdpI3w+1zaNE3QwYg2TVOyF1zO4M4hhPCLBCNCiKhR37JypV3NiJYVSc4Bc4zf59IyI0HXjCQPBnMcOJugPC+4cwgh/CLBiBAianSaGQliigZaa0YanY043I7AB6PXt67eKdkb+PuFEH6TYEQIETW0zEjnwcjIgM4Va249R32wvULSR6v3pyQYESKcJBgRQkQFRVG8ZEa0lTSjAzqfSW/CZrQB3VhR48mM7Avu/UIIv0gwIoSICo3ORtyKGwhNZgRCUDeiZUZkmkaIsJJgRIg+bkPhBh7b8hhbTmyJ9FB80nqM6HV6T0aD+nKoL1UfB7CSRtPtFTUZY9X7ynxpCy9EGEkwIkQf9knhJ9zz0T28nvc6//ff/2Nb8bZID8krT48RYyw6nU59UsuKJA0CSyc7+Xah211YY9Mgth+gtI5FCBFyEowI0Ue5FTdPfvkkAEa9Ebfi5ukvn47wqLzT6kViTG2W7wa5kkbjmaYJtmYE2kzVSN2IEOEiwYgQfdTnRZ+TX5NPrCmWt695G5PexO7y3RysPBjpoXWq0x4jZS1jDWKKxu1WaGo2A/Dcxj18f/U2lq7bx67CqsBOJEWsQoSdBCNC9FHrj60HYM7gOWTHZ3P+wPMBeO/oe5Ecllf1jk5W0pQfUu/Thgd0rm+OV3PNik/57JB6ziMVZXy0v4SVG49wzTOfctuLn1NY0eDfyaTXiBBhJ8GIEH2QoihsPrEZgAuyLwDg4kEXA/DZyc8iNi5f6hyd7EujBSOpw/w+z2tfFnL9s1v4+ng1Jp065XP+qDh+c+1Z/M/E/pgMOjbllXHF8k1syivt+oSSGREi7CQYEaIPOlp9lKL6Isx6M+dkngPA1MypAOyt2EuNvSaSw+tUg+O0HXudzVBVoD5O9S8z8soXBTz0z6+xu9zMGpPBnTPUeo/+yTpunjqI5d+ZxIeLLmByTjJ1zU6+v3obH+wp9n3S9JZ6ldoiaKgI+PsSQnRNghEh+qBdpbsAmJA+wbNMNiM2g9yEXNyKm10luyI3OC86ZEYq80Fxgzke4tK7fP+He0/xszfVHXbnnzeYlbdOJiMuud25AXJSY/n7D6Zy1fgsHC6Fe17ZyfZjld5PbIlXV/OAZEeECBMJRoTog74pUz+Uz0o7q93zY9PUvhl7y6Ov/qFDzUhZy+Z0qUNBW+rrRX5ZPQ/8YxeKAjdPHcTP54xGr9d5imFP7zNiMRpYdtNELhmVTrPTzV0vb6e0ttn7BaRuRIiwkmBEiD5od9luAManjW/3/JgU9UO1VwQjftaLuNwK963dRW2Tk7MHJfHY1WM9fUq0c2nnbsto0PPHmycxMiOesrpmHnztK9xupfOLSCdWIcIqqGBkxYoVDB48GKvVyuTJk9m0aZPXYzds2IBOp+tw279fGggJEQ6NzkbyKtWswri0ce1eG5Pa94KRl7bk81VhFfEWI8/ccjZmY+uvNW3n3rbTNG3FmI384buTsBj1fHKwlFe2FXR+Ea3PSWl0LosWorcLOBhZu3Yt999/Pw8//DA7d+5k5syZXHHFFRQUePk/cYsDBw5QVFTkuQ0fHthSPSGEfw5WHsSluEizpZERm9HutdGp6l/4pxpOUdVUFYHReacFI54+I+WH1Xsfy3pPVDXy5AfqRnqLrxxNVqKt3etaYOMtGAEYmRnPjy9Xg43f/mc/ZXWdTNdofU7KDnT5fQghAhdwMPL0008zf/587rzzTkaPHs2yZcvIzs7m2Wef9fm+9PR0MjMzPTeDwRD0oIUQ3h2pOgLA8KSOH+KxpliyYrPU46qP9Oi4uqIFI54OrJ7MyFCv73nqgwM02F2ck5vMd87J7vC6t5qR031vWg5j+ydQ0+TkN+92UqSqBSP1pbKiRogwCCgYsdvtbN++ndmzZ7d7fvbs2WzZ4nsTrkmTJpGVlcUll1zCxx9/7PPY5uZmampq2t2EEP45XKVmFIYmdf4hPiRxCKAu/40m7TIjTdVQX6K+kNL597H3ZA1v7jwBwM/njEGv71jkqmVGGpwNuNwur9c2GvT8+lq12PeNnSfYfaK6/QGWOEgYqD4uk6kaIUItoGCkrKwMl8tFRkb71G9GRgbFxZ2v1c/KyuL555/n9ddf54033mDkyJFccsklbNy40et1li5dSmJioueWnd3xLx4hROcOV6vByJCkIZ2+PjhxMBC9mZFYU2xrViQuA6wJnR7/u/f3oygwZ3wWE7KTOj1G2ygPWncF9mZidhJzJ/QH8Ez9tNOvJTtSKlM1QoRaUAWsutOW2SmK0uE5zciRI/nBD37A2WefzbRp01ixYgVz5szhySef9Hr+xYsXU11d7bkVFhYGM0whzkjaNI2WATmdFoxEW2bE02fEHNtaL+Kl2dk3x6vZcKAUg17HQ7NHej2nWW/GqDcCna+oOd2iWSMw6nVsOFDK50fK27+Y1nIdyYwIEXIBBSNpaWkYDIYOWZCSkpIO2RJfzj33XPLy8ry+brFYSEhIaHcTQnStwdHAyfqTAAxN7Hx6I1ozI54OrMbYLutFnvtEDVbmTuhPblpsp8eA+oeTv3UjALlpsdzUUnuy7MPTfkdJZkSIsAkoGDGbzUyePJn169e3e379+vVMnz7d7/Ps3LmTrKysQC4thPCDlu1IsaaQZE3q9BgtGDlZdxKHy9FTQ/NJURTqnS01I+Y4n8t6j5bV85/dRQDcdUHn2Z+2/FlR09bCi4Zh1Ov47Eg5Xx+van3BkxmRYESIUAt4mmbRokW88MILrFq1in379vHAAw9QUFDAggULAHWKZd68eZ7jly1bxltvvUVeXh579uxh8eLFvP7669x9992h+y6EEAAcqzkGQG5CrtdjUq2pWA1WFBSK6ot6aGS+NTobcStuAGKMMW26r3YMRv686QhuBS4elc6ozK6zplqvEX+maQAGJNk8tSNaBgaAfi3BSFUh2P3c8VcI4RdjoG+46aabKC8vZ8mSJRQVFTFu3DjWrVtHTk4OAEVFRe16jtjtdh588EFOnDiBzWZj7NixvPvuu1x55ZWh+y6EEAAcrzsOwMD4gV6P0el09I/rz5HqIxyvO86ghEE9NTyvtEBBr9NjM1jb1Iy0D0Zqmhy8uUNdQfO/53edFYHAMyMAd10wlDd2nuA/u4vJL6tXp4Ji08CWAo0VUJ4HWRP8Pp8QwregClgXLlxIfn4+zc3NbN++nfPPP9/z2urVq9mwYYPn6x//+MccOnSIxsZGKioq2LRpkwQiQoTJ8dqugxGAAXEDAHWqJhp4VtIYY9HVl4CjHnR6SM5td9ybO07Q6HAxIiOOqYNT/Dq3VjPib2YE1EZoF49KR1Fg9Zb81he07Ih0YhUipGRvGiH6kBN1atZgYJx/wYh2fKR5ghFzLFSqU00kDASj2XOMoiis+Vx97ZapOV5X8J1Oy4zU2msDGtPt03MBeH37cRrsTvVJ6cQqRFhIMCJEH6JlRrLjfffm8QQjtVEWjBhjoTJffTI5p90xXx6r5OCpOmwmA9eePcDvcweTGQE4b1gauakx1DY7+deulgySJzMiwYgQoSTBiBB9hMPloLhBXXbf5TRNfHRlRtr1GKlqyYwktQ9G1mxVn587oT8JVpPf5441B14zAqDX67hlqjqGlz87hqIo0mtEiDCRYESIPqKovgi34sZqsJJqTfV5rJYZ0QpeI619ZqQlGGlTL1LX7OS9PWqg9d2pgRXcBpsZAfj25IFYjHr2FtWwo6CqtddI+WFwOQM+nxCicxKMCNFHaFM0A+IGdFlPoQUjFU0VNDmbwj62rnj2pTHHtWZG2kzTvL+7mCaHmyFpsUwYmBjQuT2rafxoena65FgzV41Xl/mu3Vag1rGYYsDtgMro6mArRG8mwYgQfYQ/y3o1CeYErAYrAKUNpWEdlz+0KZQYY0xrZqTNNM1bu9TppGsmdR1ona47mRGAG6aoP8913xTT6FQgraVFvdSNCBEyEowI0UdowYiW9fBFp9OREatu4aDVmUSS1go+zhgDNS1TRy2ZkVM1TXx6qAyAayb6X7iq0TbLC7RmRPOt3BQGJtuoa3bywd7i1rqR0v1BnU8I0ZEEI0L0EcX1alDRP66/X8enx6QDUNJQErYx+cuTGXE7QXGD0aru2Au889VJ3ApMzklmUGpMwOfubmZEr9dx3dlqduSf24+31o2Ued9fSwgRGAlGhOgjTtWfAiAjxr9NK7XjTjWcCtuY/OWpGXE0q08k5UDLdMybO1unaIKhBSOB9hlp6/qWpcSfHiqj0parPqntnyOE6DYJRoToI7SgQpt+6YqWGdGCmEjyrKZpbsletEzRHCuvZ8/JGgx6HXPOCm5zTa2ANdjMCEBOaizn5CbjVuCDU+peN5TngaIEfU4hRCsJRoToA9yK2xOMZMZk+vUeLTMSDdM0nmCkqUZ9oqV49b3d6tTTuUNSSIk1d/rermg1I/WOes9mfMG4vmWqZk2eAdBBUzXUlwV9PiFEKwlGhOgDKpoqcLqd6NCRFpPm13u0DEo0TdPE1leqT7RkRv7TEoxcPi64rAi0ZkYUFBqdjUGf5/JxmRj1Or4+ZcfR0jROpmqECA0JRoToA7SpljRbGia9f91Jo7JmpL5lmXFyLkXVjewqrEKng8vG+Df11BmrwYpBZwCC6zWiSYoxM3O4GuidNLQsny6XIlYhQkGCESH6AG15bmasf1M00BqMlDWW4XRHtpuoZzVNTcsy46Qc3m/JikwelEx6gjXoc+t0upDUjQDMaWmAtrOhJfskK2qECAkJRoToAwJdSQOQYk3BoDPgVtyUN5aHa2h+8fQZaWidptHav18+zv8AyxttRU2wvUY0s8ZkYDbo2VHf0m5fpmmECAkJRoToA7TMiL8raQAMegNpNvUv/EgWsSqK0lozorjBmkS508oXRysAuGxsCIKRbjY+0yTaTJw/Io0jSksvFwlGhAgJCUaE6AO0zIi/K2k0qTb1L/zypshlRhqdjSioS2Rj3Qok5/Lf/SW4FRg3IIHslMAbnZ2uu43P2pozPosjbrWgVqk4KhvmCRECEowI0QcE2mNEo2VGIjlNo2Ur9OiwKgok5/DxfjVTc8mo4AtX2+rOZnmnu3R0BuXGNBoVMzq3o3VjPyFE0CQYEaIP0FrBB1IzAq3BSFlj5PpleKZodAZ0gCtxEJvy1PFcPCo9JNcIVc0IQLzVxHnD0slXWrJQMlUjRLdJMCJEL+dW3J6aj0BW0wCkWiM/TdNaL6K2f8939aOu2UlanJmzBiSG5Bqx5pbMSAiCEVALWQ8rLb1PZEWNEN0mwYgQvVxlUyUOtwOAfrZ+Ab1XqxmJhsxInNsFwBdVahbjghHp6PW6kFzDUzNi737NCMAlozPIbwlGGopk914hukuCESF6OS2QSLYkYzL41/BM4ylgjYKakZiWTfLeO6H2FAnVFA20qRkJUWakX7wFd8owAGpOSDAiRHdJMCJEL6cFI/62gW8rzdpSwBrBaRpPjxGXEwUdn5XHYtDrOG944N+PN6FcTaMZOHw8AJaqwyE7pxBnKglGhOjlPMGINfAP72jKjMQqCg2WftgxMSUnmURbYFkeX0LVZ6StSZOmAJDsrqCmuiJk5xXiTCTBiBC9nBaM9IsJrF4EWlfT1DnqaHI2hXRc/vIUsLrdHEedmgnlFA2EJzMyeOAAKnRqge2undtDdl4hzkQSjAjRy2nBiJblCEScKQ6z3gxEbqqmNRhR2NeYDMBFIQ5GQl0zoqmPGwzAkQO7QnpeIc40EowI0ct5MiMBrqQBdRO5SPcaaZsZyXf3IzPByvD0uJBew9NnJARNz9qyZY4EoPHkfhwud0jPLcSZRIIRIXq50sZSoHXKJVCRrhtpuy9NoTudGcPS0OlCs6RXE+o+I5qUnLEADHQfZ/uxypCeW4gziQQjQvRyWhDR3WAk0pmROLdCodKPmSFcRaNpWzOiKErIzqtPGw7AYF0xHx+I3GaDQvR2EowI0ct5VtMEG4xEuAurNnUS43ZTqKQzfVjgtS9d0YIRt+Km0dkYuhOnqr1GBuuK2LBPghEhgiXBiBC9WKOz0TP1EEzNCER+mqahSZ3esLh1JKVnkx5vDfk1bEYbOtSpn1CuqCE5F0VnIFbXTHXJMU5WhTDQEeIMIsGIEL2YlhWxGqyeFSOBSrGmAGpb+Uioa64GoMGVwPQRodml93Q6nS6km+V5GM3oknMBGKwvZsOB0tCdW4gziAQjQvRibZf1Blv0mWxRl9NWNVeFalgB0TIVte5kzhsW+noRjdb4LKSZEfBM1QzVnZS6ESGCFFQwsmLFCgYPHozVamXy5Mls2rTJr/d9+umnGI1GJk6cGMxlhRCn6c6yXk2yVQ1GKpoi00W0zqU2W6t2p/KtwSlhu064eo3QUsQ6RFfEp4fKaHa6Qnt+Ic4AAQcja9eu5f777+fhhx9m586dzJw5kyuuuIKCggKf76uurmbevHlccsklQQ9WCNFeaUP3lvVCazASicyIW3HT6HYCYIobRKzFGLZrhavXiJYZGWU6RYPdxZf5ssRXiEAFHIw8/fTTzJ8/nzvvvJPRo0ezbNkysrOzefbZZ32+76677uLmm29m2rRpQQ9WCNFed1fSQJtpmqaqkC579UejsxGlZXYpPWtUWK8Vrl4jWjAy0ngKgI/3y1SNEIEKKBix2+1s376d2bNnt3t+9uzZbNmyxev7/vKXv3D48GEeffRRv67T3NxMTU1Nu5sQoiNtOW4oMiNOxUmtozYk4/JXTbMaGBgUhZEjxof1WuHYnwbwTNOkOIox45C6ESGCEFAwUlZWhsvlIiOjfcV7RkYGxcXFnb4nLy+Pn/70p6xZswaj0b8U7NKlS0lMTPTcsrOzAxmmEGcMbZommE3yNGaD2VNP0dMrar5pmd6NcSsMHzEurNcK2zRNXAaY49HhZrC+hMOl9RRWNIT2GkL0cUEVsJ5eta8oSqeV/C6Xi5tvvplf/vKXjBgxwu/zL168mOrqas+tsLAwmGEK0eeFYpoGIMmSBPR8MLLv8DcAxCoKxrjQNztrSwu4Qp4Z0ekgdSgAl6arWdyNebLEV4hABBSMpKWlYTAYOmRBSkpKOmRLAGpra/nyyy+5++67MRqNGI1GlixZwldffYXRaOSjjz7q9DoWi4WEhIR2NyFER6EKRiLVa+REUR4AVozqh3oYhaXPiKZlqmZ6kvrz23QwMq31heitAgpGzGYzkydPZv369e2eX79+PdOnT+9wfEJCAt988w27du3y3BYsWMDIkSPZtWsXU6dO7d7ohTiDuRW3Zzmu1tI9WFpmpCdX1LjdCnVVatYzzhD6rqunC9vSXmi3ogZgy+EynLKLrxB+C3gd3aJFi7jtttuYMmUK06ZN4/nnn6egoIAFCxYA6hTLiRMn+Otf/4per2fcuPbzwOnp6Vit1g7PCyECU9Ncg0tRe1pomY1gRaLXyL7iGqyKmkGIt8SF/Xrx5nggDNM04AlGUpoKSbAaqWly8vWJas4elBz6awnRBwUcjNx0002Ul5ezZMkSioqKGDduHOvWrSMnJweAoqKiLnuOCCG6r6JZDRzizfGYDKZunSsSXVg/P1JBrEGtsYizhH8q1pMZCXUBK3iCEV15HjOGpfGf3cVsziuTYEQIPwVVwLpw4ULy8/Npbm5m+/btnH/++Z7XVq9ezYYNG7y+97HHHmPXrl3BXFYI0UZFY2imaCAymZGtR8qx6tWlxLHdzOz4I6w1Iy3BCA3lXJJjBmCTFLEK4TfZm0aIXkoLHLRAoju0aZ6eyoy43Qpf5Fdg0qtLYOO6sTTZX1rTs7BM01jiIL4/ADNTqgDYUVBFbZMj9NcSog+SYESIXkpb+dLdehHo+aW9B07V4m6oxGFQa15i4sKzW29bYc2MgGd5b4ajkNzUGFxuha1HIrPfjxC9jQQjQvRSocyM9PQ0zedHyhmkK6GhZTlvXAi+h660bXoWlrb3Lct7KT/EzOFqpkemaoTwjwQjQvRSWuAQisxIT2+Wt/VIBdm6Uur06q8grbg0nOLMajDiUlw0tewUHFJa3UhZHucNV/u+bM6TfiNC+EOCESF6qXAEI/WOeuwue7fP54tWL5KtK6G+B4ORGGMMOtRMTHhW1GiZkcNMG5qKQa/jSFk9xyulNbwQXZFgRIheKpTBSLwpHqNOXekf7rqRvJI6Kurt5BrKqNerwUFPBCM6nS7MK2rUmhEqDpNgNjApOwmQ7IgQ/pBgRIheKpQFrDqdjiRrknre5vAGI1uPqDsNj7NVUq/rucwIhHlFTVIO6E3gbIKa456pmk0SjAjRJQlGhOilQlnACj23oubzo2owkq0v6dHMCLQWsdbaa0N/coMRUoaoj8vyPEWsmw+V4XKHoWBWiD5EghEheiGX2+UpNg1FZgQg0ZIIQLW9OiTn64yiKGzLr0SHm8TmYk/NiBYkhJt2nbBkRqC1iLX8EBMGJhJvNVLd6OCbE+H7mQrRF0gwIkQvVNVchYL617aW0egu7Tw1zTUhOV9nCioaKK1tZqChBtx2TzASY4oJ2zXb0lbUhCUzApDWGowYDXqmDVG74356SKZqhPBFghEheiFtiibJkoRRH/AWU53SMiPhXN67LV+dAprZr57Glh4j0AczI2V5AMwYptaNbDkswYgQvkgwIkQvFMriVY1nmqY5fFMK24+pQdS0lFpPjxGDzoDFYAnbNdvSMiPh68LaurwXYMYwNTPyZX4lTQ5XeK4pRB8gwYgQvVCoi1cBEs09lxkZY6tsV7yqa5MlCae2XVjDQsuMVBeCo5Gh/eLISLDQ7HSz41jPtNoXojeSYESIXqi8SV2REsrMSLhrRirq7RwqUYOAgZT0+LLettcKW2YkNg2siYACFUfQ6XTMGKpO1XwqUzVCeCXBiBC9UDinacKVGdnekhkY2i8WS93xHl/WCxBvjgfCGIzodK1TNS11I9Nb6kY+PVQenmsK0QdIMCJELxTK7quacC/t/bKlXuSc3BSozO/RVvCasGdGoN3yXmitG/n6eBU1TY7wXVeIXkyCESF6od5YwPplS73IOQNjoeZkj/cYAbXtPUC9PUyraaDd8l6ArEQbQ9JicSvw+ZGe2RVZiN5GghEheqFwFLBqNSPVzdUoSmg7hjY5XHxzXA1ypqY2AAp1RnUFTU/1GIHWdvA9khlpmaYBmD5M+o0I4YsEI0L0QuGcpnEprpB/WH9zohq7y01anIUBlABQH6MGUj2ZGQnrRnkaz/LePGgJ6jxFrBKMCNEpCUaE6IXCEYxYDBZsRhsQ+qmabflavUgyuqpjANTbEoCerRkJ+9JeaN2fpqkaGtSi1WlDU9Hp1B2LS2qawndtIXopCUaE6GUcbgc1dnX5bSiDEYAEsxoghDoY0epFJuckQ2U+AHUtUyY9GoyYWzuwuhV3eC5ijoHEbPVxS91IUoyZsf3Vn+2Ww7KqRojTSTAiRC9T1VQFgF6n90ythEo4iljdbsWzrFddSaNmRhpMViAymREFhUZnY/gu1EndiEzVCOGdBCNC9DJt96XR60L7f2GtiDWUvUYOldZR3ejAZjIwpn8CtEzT1BnUPXV6MhixGCwYdep1w7ZZHnRY3gut/Ua2HC4PeYGwEL2dBCNC9DLh6L6qCUevEa1eZGJ2EiaDvjUz0tL0rCcLWHU6nWdFTdg2ywNI04pYW4ORc3KTMRl0nKhq5Fh5Q/iuLUQvJMGIEL1MOHqMaMLRhXW71l8kNxmaaqBRDU7qUDeO68nMCLQGP+HNjAxV79tM08SYjUwapK4gktbwQrQnwYgQvUw4VtJotM3yQrk/zbaWzqtTclM8UzTYUqh3qqtKIhWMhDUzoi3vrTgC7tbderW6kS3SGl6IdiQYEaKXCWdmJNQ1I8XVTRRWNKLXwaRBSZ4pGpJzPMGAtsKlp2jXC2uvkcRsMFjA7WgNwGhtDb/lcBlut9SNCKGRYESIXiYc3Vc1oV5No+1HMyozgXirqfWDOTnXEwzEGiOTGQlrrxG9vnWqpvyw5+kJ2UnEmg1UNjjYVxye3ZGF6I0kGBGil+mRAtZQBSNt60XAkxlxJQ7yLK3VCkp7So9kRqDT5b0mg55vDVb/u8lUjRCtJBgRopfpiQLWUK2m0TIjk3NbxtrS8KwhMctzTKRqRnosGGmzogZgRssSXyliFaKVBCNC9DLhLGANZc1IXbOTvSfVqQhPZkRrBR/XDwCj3ohZb+72tQKhBT9hnaaBNst789o9Pb2liPWLoxXYnWHqAitELyPBiBC9TE/UjNQ013S7XfqugircCgxIspGVaFM3jasqAKA+Vv1AjjPFodPpujfoAMWb49UxhHM1DbSZpmmfGRmVGU9KrJkGu4uvjleFdwxC9BISjAjRizS7mj0fouFc2qugdLsPR9vN8QCoLwVHA6CjzqoGBD09RdP2mj02TVN7Eppbr6XX65g2VF1VI63hhVAFFYysWLGCwYMHY7VamTx5Mps2bfJ67ObNm5kxYwapqanYbDZGjRrF73//+6AHLMSZTKsXMeqMnk3tQslkMBFjjAG6X8TqrV6EhAHUu5qByAQjPbKaBiAmBWLUoIOKw+1ekn4jQrQXcDCydu1a7r//fh5++GF27tzJzJkzueKKKygoKOj0+NjYWO6++242btzIvn37+PnPf87Pf/5znn/++W4PXogzTdspmnBNb4SibsTpcrOzQH3/lJz2K2lIzqHe2dJjpAdbwWt6rIAVfBSxqkHKzsJKGuzO8I9DiCgXcDDy9NNPM3/+fO68805Gjx7NsmXLyM7O5tlnn+30+EmTJvHd736XsWPHkpuby6233spll13mM5sihOhcOItXNaFY3ru/uJYGu4t4q5GRGeqUDFX56n1SjicrEWOK6c5Qg9JjS3uhtRPraXUjg1JiGJBkw+FS+OJoRfjHIUSUCygYsdvtbN++ndmzZ7d7fvbs2WzZssWvc+zcuZMtW7ZwwQUXeD2mubmZmpqadjchRHiLVzWh2J/my5Z6kbMHJaNv2RCvNTOSS4NT3SgukpmRenuYC1ihTeOz9sGITqdr041VpmqECCgYKSsrw+VykZGR0e75jIwMiouLfb534MCBWCwWpkyZwg9/+EPuvPNOr8cuXbqUxMREzy07OzuQYQrRZ4Wzx4gmFJmRL4+p4/RM0UCb7qutmZFI1ozUOsK4UZ7Gy/JeaO03sjlPiliFCKqA9fS5akVRupy/3rRpE19++SXPPfccy5Yt45VXXvF67OLFi6murvbcCgsLgxmmEH1OOLuvajyb5dmDz0hubwlGJrcNRrQC1qTWfWkiEoy0TNM0OhtxtdnELiw8NSOH1aXNbWj9RvYW1VBRbw/vOISIcsZADk5LS8NgMHTIgpSUlHTIlpxu8ODBAJx11lmcOnWKxx57jO9+97udHmuxWLBYLIEMTYgzQk9mRoKdpjlR1UhRdRMGvY6Jg5LUJ11OqD6hPk7Oof7kh0BkpmnaBkB1jjrP9xsWKUNAp4fmGqgrgfjW35P94i2Myoxnf3EtWw6XcdX4/uEbhxBRLqDMiNlsZvLkyaxfv77d8+vXr2f69Ol+n0dRFJqbmwO5tBCC3lHAqtWLjMlKIMbc8vdOzXFQXOpOtnGZnuLRSBSwmg1mT9fXsDc+M1ogaZD6uJOpmvNkqkYIIIhpmkWLFvHCCy+watUq9u3bxwMPPEBBQQELFiwA1CmWefPmeY5/5plneOedd8jLyyMvL4+//OUvPPnkk9x6662h+y6EOENUNPZcAWuw+9N0PkXTUi+SNAj0ehockStghZ5eUdP58l6AGcPVYGRTXhnKadM4QpxJApqmAbjpppsoLy9nyZIlFBUVMW7cONatW0dOTg4ARUVF7XqOuN1uFi9ezNGjRzEajQwdOpTHH3+cu+66K3TfhRBniMrmHpimMbe2hA+GtlPvlNxO6kWS1d8TWhDQ0zv2auJMcVQ0VYS/8Rmoy3sPfdhu917Nt3JTMBl0nKhqpKCigZzUyPw8hIi0gIMRgIULF7Jw4cJOX1u9enW7r++55x7uueeeYC4jhDhNtE/T1DU72V+sBjFTctqMUVtJk6QGI54CVmOEgpEezYxoy3sPd3gp1mJk0qBkvjhawaa8MglGxBlL9qYRopdodDbS6GwEwhuMdKcDa9vN8TITra0vtOm+Cq3BiBYU9LQeawkPPpf3AsxsqRuRfWrEmUyCESF6CW0ljUlvCuuS2ASLuudNrb024KWv2n407aZooE2PkVygzTRNBJb2Qg+3hE8bod5XHAVnx8J9rW5ky+FyXG6pGxFnJglGhOgl2k7RhGtfGmi/c2+gH9bbO2t2Bm0KWE+bpolUMNKSkQn7ahqA+CywJKiriTqZqhk/IJF4i5HqRgd7TnZvc0IheisJRoToJXqiXgSC37nX5VY8m+NNblsvYq+H+hL1cXIODreD5pZdeyO1mkYLgmrtPdCFVaeDfiPVx2UHOrxsNOg5d6jaGn6TLPEVZygJRoToJXpiXxqNVjcSSDCyv7iGumYn8RYjIzPjW1/QsiLWRLAle5b1QmT6jECb/Wl6IjMCkNYSjJR2DEYAZg6XuhFxZpNgRIheoie6r2qC6cKqTdFMHJSEQd9mGum0ehEtALAYLJj0pm6PNRg9upoGWjMjpfs7fVnbp+bL/Eoa7WFuUS9EFJJgRIheQgtGeiIzohWxBtL4zNNfJOe0YMnTYyQXiHzxKvTwahpoE4wc7PTlIWmxZCVasbvcniJgIc4kEowI0Uv0xCZ5Gq2INZBpmk47r0KHYCTSxasQgWkaLRgpz1P36TmNTqeTXXzFGU2CESF6iUhM0/jbhbWoupETVY3odbRujqdps1svtOkxEqHiVWidpql19EABK0DiIDDawGVvnbY6jWefGqkbEWcgCUaE6CU80zSWHixg9XOaRpuiGZ2VQJzltMbOUTxN02OZEb2+tflZF3Uje07WUFFv75lxCRElJBgRopfwLO21RV8Bq9f+IorSpvtqLoBnNU0kg5EeXdqr6TdKvfcSjPSLtzCqZRXSlsOSHRFnFglGhOglPJvkWcIfjCSYWwpY/awZ8dSL5J42troScDaCTg+J2epT9shnRuLN6od+j2VGAPq1dGL1UsQKSN2IOGMFtVGeEKJnNTgaWvel6cHMiD81I/XNTvYWaZvjeSleTRgIRrN6fBTUjGiBULOrGYfLgcnQA0uMvWRGFEXh30f+zRfFX2CMywYy2JRXhqIoYe20K0Q0kWBEiF5Ay4qY9WZPd9Rw8uzc60fNyFeFVbjcCv0TrfRPsrV/0VMvkuN5KhpW07S9dp2jjmRD+OtwPI3PyvLA7VbrSIBff/5r1h5Y6zksZtAIThR8j4KKBtnFV5wxZJpGiF6gbY+RnvhrOZAOrF96m6KBToORaChgNeqN2Iy2duMJu5TBoDeBox5qjgPw4bEPWXtgLXqdnuuHX4/NaMMQexBz2gZZVSPOKBKMCNEL9NS+NBrPNI29Brfi9nnsl96KV6HDShqIjswIRKDxmcEEqcPUx6UHcbldLN+xHID54+bz2PTH+MW5vwDAnLqBjw523FRPiL5KghEheoGeDka0Ala34va54sTlVtjprdkZtGkFP9jzVNQEIz3dEh7aFLHu5+PCj8mvySfBnMD8s+YDcNWQqxgSPxqd3sEXFW/hdPkOBIXoKyQYEaIX6MmGZwBmg9kzjeGriPXgqVpqm53Emg2eZant+MiMaMFApMSb1PFGannvW4feAuDGkTd6AjOdTscPJ/0AAHfcNrYXyFSNODNIMCJEL9CTO/Zq/Cli1aZoJg1Kxmg47deJowlqTqqPkzopYDVGNjMSb4lAMJKmZkaqyvbz6YlPAbh6yNXtDrko50LMJKE31rHmm//03NiEiCAJRoToBSIRjGhFrL4an32Zr46r0yma6kJAAVMsxKZ5nvYUsJojG4wkmNSpqEhkRv5bdxSn4mRUyiiGJA1pd4hJb2Jy6qUAfFHycc+NTYgIkmBEiF5AC0ZSrak9ds2uNstTFIUvjqrjmjrY10qaXGizAiga+oxAa+OzHg1GUoeBTs+nLW1NLh50caeH3TLuKgDqDN9QVNOD4xMiQiQYEaIXaLu0t6ckWHx3YT1e2UhRdRNGvY5Jg/xbSaMoSlQs7YXWYKTG7t9mgCFhsuJMzmGr1QrAjP4zOj3s/JyzMbiT0OntvLTz/Z4bnxARIsGIEL1ATxewQtc1I9tapmjGDUjEZjZ0PKCTYKTZ1YzT7QTO0MwIsCc1l1qDnni9hbGpYzs9RqfTMSR2KgCfFG7uyeEJERESjAjRC0SkgNXsuyW8Fox8q7MpGug0GNGyIjp0xJjC30nWl4hkRoCtMeoqpXMN8Rj0nQRxLWblzgTgRNPXuN1Kj4xNiEiRYESIKNfgaKDJ1QT0bGakqwLWz1vqRb7VWedVaLNbb+tKGi0LEWeKQ6+L7K8frZdKT2dGdtEMwOTGJp/H3TjuQhRFj2Iq5ZMj3jfXE6IvkGBEiCin7UtjMVh6ZF8ajWeappOakbK6Zo6UqoWoU3I7ydYoSueZkZZup5HuMQKRCUYUReHrxiIAJpYfV/eo8SI1JpF4nbrS5vV9sqpG9G0SjAgR5SoaW6doenIXV08Bayc1I9qS3pEZ8STFmDu+uaECtA/5pEGep2sdLZmRKAhGIlEzcqzmGDWOOixuhRENNa0dar0Yn3oOADtLtvXE8ISIGAlGhIhyWmakJ6dowHfNiGeKxlu9SFW+eh+fBabWnXy1zIjW/TSSIhGMfFX6FQBjFCMmgJK9Po+/ari62qZKOUhNkyPMoxMiciQYESLKlTeWAz1bvAq+p2m04tVzuipebdN5FVoLWKMpM1LnqOtyM8BQ+br0awAmWDPUJ075DkYuGXIOKDr0pirW7d0X7uEJETESjAgR5TyZEUvPZka0AtZqe3W7D+vaJgd7T6rZEq/FqxVH1PuU9t1F2xawRpoWjCgoPbZZ3jdl3wBwVupo9YmSPT6PjzHFkGRUA7r/5H0e1rEJEUkSjAgR5SLRYwRaa0bcirvdh/X2Y5W4FRiUEkNmorXzN1ccVe9PC0a082iBQCSZDWasBnX8vjYDDBWn28mhqkMAjO5/rvpkF5kRgLPSJgDwTdkuFEWW+Iq+SYIRIaJcJHqMgLp6R9u5t+1UjWeKxltWBKD8sHqfMrjd057VNFGQGYGerRvJr87H4XYQY4xhQHZL59XyQ+Bs9vm+2UPV5mdNhiMcPNUzGRwhepoEI0JEOS0Y6enMCLQuf22bOdh2VM3UfGuwj+Coq2maKKgZgZ5d3nug8gAAI5JHoE8cCNZEUFxQesDn+87JOhsAvfUk/9lTEPZxChEJQQUjK1asYPDgwVitViZPnsymTZu8HvvGG28wa9Ys+vXrR0JCAtOmTeP992WvBSH8Fclg5PQi1iaHi12FVYCPzEhzLdSXqI+9TdNEwWoa6NnMiBaMjEwZqW4cmN7SCr7Ed2Fq/9j+xBqS0OncvJ+3PdzDFCIiAg5G1q5dy/3338/DDz/Mzp07mTlzJldccQUFBZ1H7Bs3bmTWrFmsW7eO7du3c9FFF3H11Vezc+fObg9eiDNBJDbJ05zehfXr49XYXW7S4iwMTvOy0Z1WLxKTCrakdi9FU9Mz6NmW8Acr1S6qI5JHqE9kjFHvuyhi1el0jEtTA5fDNfspq/M9rSNEbxRwMPL0008zf/587rzzTkaPHs2yZcvIzs7m2Wef7fT4ZcuW8eMf/5hzzjmH4cOH85vf/Ibhw4fzzjvvdHvwQvR1iqJQ1lgGQJotrcevf/pmea370fhowOZligZam55FQwEr9Gxm5GDFacFIeksw4kcR66TMswDQW07w8f6SsIxPiEgKKBix2+1s376d2bNnt3t+9uzZbNmyxa9zuN1uamtrSUnxnnJubm6mpqam3U2IM1GtoxaHW212lWpL7fHrazUV2jTN1iNqzxOvS3rBZzASrQWs4c6MVDRVUNpYCrTNjGjTNF0HI9ruvnrrCf67T4IR0fcEFIyUlZXhcrnIyMho93xGRgbFxcV+neOpp56ivr6eG2+80esxS5cuJTEx0XPLzs4OZJhC9Blaw7M4UxwWg6XHr9+2ZsTudPNlvjplNG2ojyxNhbaSppNgJIqankHPFbAeqlSX9A6MG9i6W3F6S6+RmhPQWOnz/WNS1SyK3lLCpkMnaHa6wjZWISIhqALW09OziqL4tWfGK6+8wmOPPcbatWtJT0/3etzixYuprq723AoLC4MZphC9nhaMRCIrAq3BSI29hq+OV9HocJESa2Z4uo9gwtNjZGiHl7QP/WgpYO2pYORotfozGZrU5mdiTYTElj+0upiqSY9JJ82Whk6n0Kg/zudHKsI1VCEiIqBgJC0tDYPB0CELUlJS0iFbcrq1a9cyf/58/vGPf3DppZf6PNZisZCQkNDuJsSZqLypJRixRiYYaVvA+tlhdSznDklBr/fxx4eXaZpmV7NnyilaMiM9VTNytEYNRgYntu+7QkstCMVfd3kObarGYD3Of/edCun4hIi0gIIRs9nM5MmTWb9+fbvn169fz/Tp072+75VXXuH222/n73//O3PmzAlupEKcgbTi1YhlRsyt0zRavci0IT7GYq+H2iL18WkNz7QPfB06Yk1eVuL0sJ6qGdEyIx2DkfHqfVHXwYg2VWOwnuDDfSXSjVX0KcZA37Bo0SJuu+02pkyZwrRp03j++ecpKChgwYIFgDrFcuLECf76178CaiAyb948li9fzrnnnuvJqthsNhITE0P4rQjR93imaSKUGdFawlc1V5N3TKsX8TEWbYM8axLEtC9y1YpXY02x6HXR0W/Rkxlx9Mw0TYdgJKslGPEjM+IJRmwnOVHUyIFTtYzKlKyx6BsC/o1w0003sWzZMpYsWcLEiRPZuHEj69atIydH3cypqKioXc+RlStX4nQ6+eEPf0hWVpbndt9994XuuxCij9IankW6ZqSisYpmp5t+8RaG9vMxxVLee4pXoWdqRhocDRTVq9miwQleMiOl+7tsCz8qZRQABnMp6Jyyqkb0KQFnRgAWLlzIwoULO31t9erV7b7esGFDMJcQQhAFBawt0zR1jhpA4dwhqb6L1bV6kVTvxavRsqwX2kzThHGjvGM1xwBItiSTZE1q/2LiQLAlq6tpSvZC/0lez5MRk0GiJZHq5mr05hLW703jhxcNC9u4hehJ0ZErFUJ0KtIFrFpmRMEN+mbf9SLgu8dIFO3Yq9HG0uBswOl2huUaXqdoQG0L72fdiE6n8/Qo0VuL2FVYRVF1Y0jHKkSkSDAiRBSLdAGr1WjFYrACoDM0cO6QLvbH6UUNz6B9YKSNL9S8rqTRBFA3MjJ5pPqWfur03Qd7ZFWN6BskGBEiSimK4pmmiUQreI3NoAYPaQlO7/vRaDw9RjppBR9lO/YCGPVGYoxqE7Jw1Y34zIwAZE5Q7/1YUaNlRuLj1W6u/9ld1P0BChEFJBgRIkrVOeqwu+1A5KZpAHRu9cN61ACj73oRRyPUHFcf+5qmiZKGZxpP3YgjPHUjXQYjWmbk1G5w++6sOiJFDUaq3QWAwhdHKyiXjfNEHyDBiBBRSsuKxJpisRqtERtHU7N67cHpXXRZ1qZorInqjr2nicbMCLQGI9r+O6HkVtyeAtbchNzOD0odBqYYcDS0rkbyYmjiUPQ6PTX2KkYNALcC6/fKVI3o/SQYESJKRbp4FaC+2UldgxmAzGS374PL1F1pSR2uFmaeJhoLWKG1y2w4VtSUNJTQ7GrGqDPSP65/5wfpDa2b5nVRN2I1Wj1Bzfgh9QC8t8e/fcGEiGYSjAgRpSJdvArw+dFyXE4bAAZjFys3ytTN4Egb3unL0VjACu03Awy1wlp1X62suCyMeh+dFDwrar7q8pxaEWu/VLWI9dNDZVQ3Oro3UCEiTIIRIaJUpLuvAmzKK0NxqTUj1fYuPqzL89R7L8FItE7TtN1/J9S0YCQ7voudxwNYUaPVjZTZ8xmWHofDpfDxfmmAJno3CUaEiFKeaZoIZkY25ZWBW82MdJk5aDtN0wlt/xetkVq0aG15XxXyc/sfjGgrar6CLvac0VbUHKw8yOVjMwFZVSN6PwlGhIhSke6+WlTdyKGSOtAyI76CEUXpcppGC0a0D/9o4akZCcNmeX4HI+ljwWBWO7FWHvV5qBaMHK0+yiVj1L4vnxwspcEenqZtQvQECUaEiFKRLmDdlKfWrAxK6gd0EYzUnQJ7Lej0nS7rhdYCUW0/mGihZWrCmRkZGD/Q94FGc2vdyIkdPg/V2sK7FBcWWyk5qTE0Odyyqkb0ahKMCBGlIp0Z2dwSjEwckAV0UTOiTdEk5YDR0uFll9vl2Rk32oIRLTMSzgLWLjMjAAMmq/cntvs8rG1b+INVB7l6vLpK552vTgY/UCEiTIIRIaJUJAtY3W6FzYfUYOTcXPWD1OeHdZnv4lVtWS9E3zSNNp5QByPVzdWeot2BcV1kRgAGnK3edxGMQOuKmrzKPOZOVIORTw6WUtVgD26wQkSYBCNCRCFFUTzTNJFoBb+3qIaKejuxZgNTc9RgpKa5BsVbcWW5Vi8yotOXtSkam9GGSW8K+Xi7I1yZES0rkmZLI8YU0/UbtMxI0Vfg8r1Ut20R64iMeEZlxuNwKby3W3qOiN5JghEholCNvYZml9rmu19Mvx6/vicrMiSVtJgkAJyKk3pHfedv8Kyk6XxLe0/xapRN0UCbPiP2atxKF43dAhDQFA1AylCwJIKzCUr2+Ty0bTCiKApXT1CzI2/LVI3opSQYESIKlTaoG6ElWhKxGDrWYITbpjz1+jOHp2Ez2jxj8Fo30sU0jfa+aJuigdZgxK24200ndVfAwYheDwMmqY+7mKoZmqS2ha9oqqC8qZy5LcHIZ0fKKalpCnrMQkSKBCNCRKGSRrWJVT9bz2dFGu0utuVXAnDecPX62oqTTqcyHE1QVaA+7qLHSDRmRiwGCzZjSy+VptBN1QQcjIDfRaxWo5WchBwADlYcJDslhrMHJaEo8O+vpeeI6H0kGBEiCmmZkfSY9B6/9mdHyrA73QxIsjG0XyzQRZFnxWFAUacY4jofb7Qu69W0naoJlaCCkf5aEavv5b3QfqoG8GRHZKpG9EYSjAgRhUob1WAkEpmRj1pai180qh+6lg3vfH5Ye6ZohnW6QR5Ed2YEwtNrpFuZkdJ9YPdSn9Pi9GBkzvj+6HWwq7CKY+W+3ytEtJFgRIgoVNLQMk3Tw8WriqLw8X41ELp4VGuWwzNN09k0hla86mUlDURv91VNqFfUNDmbPP8NAwpGErIgvj8obji5y+ehWjByoPIAAP3iLcwYpq68en378cAHLUQESTAiRBTSpml6OjOSV1LHiapGLEY904a0LilOsiYBXjIj2sqPfqO8njfap2lCvT/N8Vo1GIgzxXkCHb8NbMmOFH7u8zAtGDlSfQRHy1LgG6aogc/rO07gdvve40aIaCLBiBBRSCtg7emaEW2KZtrQVGxmg+d5nwWspfvV+/TRXs8b7dM0nv1pmkOzP03bKRqdl6krrwZNazmJ72AkKzaLeFM8TreTozXqfjazx2QQbzVyoqqRz46UBzxuISJFghEhopAnM9LD0zRaMNJ2igZ8FLC6HK01I74yI71kmiZUmRG/96TpTPa56n3BVnB773ui0+kYnqyuXtLqRqwmg6eQ9bUvCwO/thARIsGIEFHGrbg9Bazptp7LjFQ3ONh+TF3Se9HI9tf1WsBafhjcDjDHQaL32ohon6bRvr9QByMB1YtossaD0QZNVa31OF6cXsQKrVM17+0ppqbJdydXIaKFBCNCRJmq5iqcbnU7+J5sBb8xrxSXW2F4ehzZKe3bl3udpinV6kVGqk27vIj2aZpQL+0trOtGMGIwwcAp6uOCz3weOiKlJRipaA1GJgxMZFh6HE0ON+9KzxHRS0gwIkSU0aZoUqwpmAw9t4/Lxwe0Jb0dszHaNEZlU2X7F0pa6kX6ea8XgTaZkSidpvG5WigIWgFrUMEIwKCWqRo/i1jbZkZ0Oh03TFanh2SqRvQWEowIEWUi0WPE5VbYcEC97ulTNKAGRgCVzacFI1pmJN17vYjL7aLWoe5eG62ZEZ+rhQLkcrs4UXcCCEEw0kVmZHiSWjNS2lhKRVOF5/lrJw3AoNexo6CKvFO1wY1BiB4kwYgQUSYSxavb8iuoqLeTFGPinNzkDq8nW9XnqpurPVNIgF+Zkbb7vWgZiGjjqRlpqur2uU41nMLpdmLUG8mIyQjuJAPPAXRQmQ+13nfijTHFeAKevMo8z/PpCVZPEfKazwuCG4MQPUiCESFCQVHUrd83PQ3/+B68dDW8NBfeWgjbX/L5gXI6T8OzHsyMvL9HHd8lozIwGjr+WkiyJKFDXaLqKfJ02ltaweMzM6JN0diMth6ddgpEqjUVgFpHradnR7C0KZoBcQMw6A1dHO2FNREyxqmPC7b6PHRk8kig/VQNwG3nqnvXvL79OA12Z4f3CRFNjJEegBC9mtsNe95Qg5CSPZ0fs2sN6I0w7nq46GFIzvF5Ss80TQ9lRhRF4YM9pwC4fFxmp8cY9AaSrcnqLrGN5WphbfkhcDvBkgAJA7ye37Njb5RO0QDEm+Mx6Ay4FBcVTRVkxAaZ0QCO16nBSFDLetsaNBVOfaMGI2Ov8XrYiOQRfFjwIQcqDrR7/rxhaeSkxnCsvIG3d53kO98a1L3xCBFGkhkRIlhFX8Hz58Pr89VAxGCBUVfBrCVw/Ytw7Uo4/8fqfiNuJ3y9Fp6ZClufVTMpXpyqVwODnlrWu/tEDSeqGokxG5g53PvqnWSLOlXjqRspbdN51UdjL63oVas7iUZ6nd4zFdW29iIYnh4jcd0MRnJmqPf5m3we1lkRK4Ber+OWqWoA8vLWYyg+/s0JEWmSGREiUG43bH4aNixtyQwkwox74JwfgC2p4/EXP6zuwvrBL+DYZnjvp3B8G/zPM2CydTi8uEGdMsmKywrzN6LSpmguGNEPq8n7tEKKLYXD1YepaGz5sC7pungVWqd1tLqMaJViTaGssazbwUi3V9Jocmeq96d2Q30ZxHYeKGrByOGqw55aFc0Nk7N58oOD7DlZw67CKiYN6lgPJEQ0kMyIEIGw18Nr34OP/p8aiIyeC/dsh/Mf6jwQ0Qw4G27/N1zxO3XKZvfr8Mp3wN7Q4dDiejU4CLr4MUBaMHLZ2M6naDRaZsPzYV28W73Xahu80IIRLbMSrTp8f0HSgpFuZ0bi+kH6WPWxj+zIgPgB2Iw27G47BTXti1WTY81cNV4Nav+2VQpZRfSSYEQIf9WXwV+ugH1vg94Ec/8IN/5V/dDwh04HU++C294EUywc2QCvflctBG3R6Gz0fHj3RGbkcGkdeSV1mAy6TvuLtKUFE63ByDfqfeZZPt+nTdNoy2ejVciCkVDVjAAMPl+9P/KJ10P0On2HtvBt3dpSyPrO1ycprW3u/piECIOggpEVK1YwePBgrFYrkydPZtMm71F7UVERN998MyNHjkSv13P//fcHO1YhIqf2FKy+Sq0TiUlTsxxnz/NZK+HV4PPhtjfUFupHNsC7D3hqSLSsSIwxhnhTfAi/gc69t1u93rShaSTafK90SbG1+bBuqICalm3qM8b6fN+ZlBmptdd6vt+QBCNDLlDvj270eZg2VXOg8kCH1yZlJzExOwm7083Ln+V3f0xChEHAwcjatWu5//77efjhh9m5cyczZ87kiiuuoKCg8xRgc3Mz/fr14+GHH2bChAndHrAQPa6mCFZfqRZsxmfB999rbUoVrEHnwrf/Ajo97PwbfPYM0BqMZMVmBb7baxDe+eokAFd6WUXTVoqlzYe1lhVJzlWXofrQm2pGoHvBiDZFk2JNIdYU2/1B5UxX/41UHIbq414P87a8F9SOrP97/hBALWRttLu6Py4hQizgYOTpp59m/vz53HnnnYwePZply5aRnZ3Ns88+2+nxubm5LF++nHnz5pGY6N8vo+bmZmpqatrdhIiIxir42/XqMtbEbLhjHaQND825R8yGyx9XH3/4KJzY7glGMmO7Dg66K+9ULfuLazEZdF6X9LalZUYqmyr9nqKBNpkRa9/PjHimaLpbL6KxJkL/SepjH9kRbytqNJeNzSQ7xUZlg4N/7vAe1AgRKQEFI3a7ne3btzN79ux2z8+ePZstW7aEbFBLly4lMTHRc8vO7mZVuhDBcDTCK99Vl+3GZcD33oGUIaG9xrf+F8ZcoxbD/vP7FFcfA3omGNGyIucP70dSjLnL49t9WHuCkfFdvs9TM9Kyv0208rS8P33/nQB4ildDMUWj8aNuRKsZKa4v7riZIWDQ65g/YzAAqzYfxeWWZb4iugQUjJSVleFyucjIaF/ln5GRQXGx/x0mu7J48WKqq6s9t8JC2exJ9DC3C16/Ewq2qE29bn0dUgaH/jo6HVy9HJIGQWU+RXnrgPAHI4qi8HZLMDJ3Yn+/3tOuD0dfzIzYQjdNE9JgZMiF6v2Rj732p4k3x9M/Vv3v2LYtfFs3TMkmwWrkaFk9H+47FbrxCRECQRWwnj6XrShKSOe3LRYLCQkJ7W5C9Kj//hL2/xsMZvjO3/360A2aLQmuUac5i6uPAGrNSDh9c6Ka/PIGrCY9l472bwmx1jK9zlGHvaylULKLn4uiKJ79XnpLZqQ7wUjIGp61NWiauvqq7hQUf+31MF9FrACxFqNnZc2KDYelCZqIKgEFI2lpaRgMhg5ZkJKSkg7ZEiF6rd1vwKfL1cfXPAuDZ4b/mrnnweTbKTaoDasyLeHtVqpN0VwyOoNYi3+9D+PN8Rh16rEVOgVsyT7bwIMauDgVdV+U3lLA2uhspMHRsf+LP0K6rFdjtLRmR/I+8HqYr+W9mjtmDMZq0vNVYRWfHCwN3RiF6KaAghGz2czkyZNZv359u+fXr1/P9OnTQzowISLi1B741w/Vx9PvhbO+3WOXVi55jGJTSzBywPuHTne53Ar//roIgLkT/JuiAbWfhdYrpMKgV7MiXWREtSkam9GGzdix22w0iTHGYDFYgOCyI063k6I69efa7e6rpxs+S73PW+/1kNGp6s7J+8r3eT2mX7yFW6eq2ZHl/82T7IiIGgFP0yxatIgXXniBVatWsW/fPh544AEKCgpYsGABoNZ7zJs3r917du3axa5du6irq6O0tJRdu3axd+/e0HwHQoRKQwW8ejM4GtS/RC95tEcvX2PQ09jy4Z7x+QvqkuIw+PRQGUXVTSTaTFwwIrDN+DxFngaDX8Wr2hRNtGdFQJ1+7k4R66mGUzgVJya9ifSYEO8rpAUjx7ep/047MTZV7feSV5lHs8t7c7P/vWAIFqOenQVVbD5UFtpxChGkgIORm266iWXLlrFkyRImTpzIxo0bWbduHTk5arRdVFTUoefIpEmTmDRpEtu3b+fvf/87kyZN4sorrwzNdyBECLhdDt77500sNNczZ9Agvptk5rndL1Bnr+uxMZysU6dOUhQ9VkcD/HdJWK7z2nZ1KuF/Jvb3uRdNZzxFrAa9ugFgF7RN9aK94ZmmO3UjWr3IgLgB6HUhbm6dOFBtu6+44dB/Oz0kKzaLFGsKTsXZYQffttLjrdyiZUc+lOyIiA5B/T9m4cKF5Ofn09zczPbt2zn//PM9r61evZoNGza0O15RlA63/Pz87oxbiJApayzjB/+4jIeUYjbF2CgwwO7K/Tyz6xn+563/YU/Znh4Zh6feILFl1c5Xf4cT20N6jeoGh2cvmhunBD6VkGJOAqBCb/ArGNGmaaK9eFWTalOLdMsaA88YhGyDPG88UzWdT+HpdDrGpI4BYE+573+zd10wBLNRz5fHKtmUJ9kREXmyN404o52sO8n3/nU9X9hLsbnd3JV5PqsuW8X/m/H/GBQ/iJLGEu54/w72lod/WtGzLDR1JIz/jvrk+tBOFb391QnsTjejMuMZ2z/wVWppbvW+zBqrLkfuQm/Zl0bTz6ZOW5U0lgT83rAs621reEt/p7wPwOXo9JBxaeqmhbvLdvs8VUaCldtaVtYs/c9+3NJ3RESYBCPijFVrr+X/3vs+Bc0VDHA4+Uf6pdx92TOck3kO1wy7hrVXreWczHNodDZyz0f3dHsDta60+zC7+OfqsuL8Ter+NSHyjy/Va9wwJTuo5fj9mtRpq9K4NL/25SlvLAdalwVHO63Wo7Qh8JUmIe++errsqRDbD5qqvO7iq9WN+BM8333RMOKtRvYV1fDWrhOhHKkQAZNgRJyR3Iqbn2xYxJH6E6Q7naw2DSb38qfaHRNnjmP5RcsZnDiYkoYSfvP5b8I6pnYfZknZMPl29YWPfuW12VUg9hXV8M2JakwGHdf42ejsdP1q1YxBqdm/lTHlTWowkmZLC+p6Pa1fjJoZCSYY8fQYCVdmRG+AUVepj/f+q9NDtGDkSPWRLpcnJ8eaWXjhMACefP8ATQ7Zs0ZEjgQj4oz0yt41bCraitXt5g/1RjK//TIYOvbbiDfH8/jMxzHoDLyf/z4bCjeEbUwd0vwzHwSjTV1BcfD9bp//b1vVVvOXjs4gNc4S1DnSy/MBKNH7FxxptRdaLUa0S7epmZFTDYF1KFUUhYIatXA/JyEn5OPyGDNXvd/3b7VL8Gn6xfQjPSYdt+JmX4X3Jb6aO2bkkpVo5WR1E6u35Id4sEL4T4IRccY5Wn2U33+pZkEerKpn7A1/g1jvH5ZjUscwb6y6XP3p7U/jdDtDPiaX2+VZTeMpgIzPgKn/qz7+6Ffgdgd9/upGB2/sUFPxt00L8sOyqZq0SjWgKXPU+/UWbZqmt2RGPNM0jYFlRsqbyqlz1KFDF77MCEDuTLXZXEMZHOt8P7BxqWrdiD+F11aTgR/NVnf8feajQ5TUNIVurEIEQIIRcUZRFIXHPryHZlxMa2zkxot/C1kTunzfD876AcmWZI5WH+XNQ2+GfFxte1RoRZQAzLhf3Rvn1Dew/52gz//69uM0OlwMT49j2pAgsxQndpDuVP8ar3PW+9WlVMuM9JZgRJumKW8sDyjo1LIi/eP6exqnhYXBBCPnqI+9TdWkqVM1u8t9F7Fqrps0gAkDE6ltdvKbdV1nU4QIBwlGxBnlw69Xs6PuGDa3m19mX4Vuwk1+vS/eHM//jlezFC9+82LIsyPaFM2AuAEY9G16f8SkwNS71McbnwyqdsTtVjxTNPOm5wa/j1TBZ8QqCraWXxtdZQ9cbpenz0hvKWBNsaZg0BlQUDxZHX8cq1F/voPiu15h1G3aVM3ef4Gr479DbUXN16Xe97FpS6/X8f+uGYdOB2/tOsmWw7LUV/Q8CUbEGcNRV8Lvtz8NwPd0yWRd/mRA779+xPUkW5I5UXeC9ce8t+UOhla8OiC+k71epv6fulFa8dc+24F7s/lQGUfK6om3GLluku+9ZHw6tgUdkG5WlwR3VeRZ1VyFW3GjQxf1O/Zq9Dq9J4tT0uD/8l4tGAlrvYhmyEVgS4H6kk5XWo1PG49ep+dE3Qm/v4fxA5O4ZaoaSD3yrz3YncFPCQoRDAlGxJnB7eLVN26i0KD2yrjj2lc6LVj1xWa08d3R3wXgL7v/EtLOlVqav9NlobGpcM731ccbfxdwduTFzUcBuH7yQL83xevAaYfjXwKQ1rKjcFeZEW2KJtmajFEf5HUjQKsbCaTXiBaM5CbmhmNI7RnNcNYN6uOv/t7h5ThznGcH350lO/0+7UOzR5Eaa+ZQSR3PfXI4JEMVwl8SjIgzQvX6h3nOqa6QuHvcncQkBtcl87sjv4vNaGNfxT52lOwI2fiOVqsBw2Ct++rppt0DRqu6suboRr/Pu+dkNZ8cLEWvg+/P8HJufxR9Bc5GsKWQnqj+9d/VX93ast7espJGE0yvkfyafKCHpmkAJrQ0xdv/LjRVd3h5UvokILBgJDHGxCNXqx1c//DfPPac7HheIcJFghHR933zT/584BVqDAaGWdO5ZvLdQZ8qyZrEFYOvAOC1g6+FaoQcrekiGInPgLO/pz7e+ITf533ukyMAXDW+P4NSY4IfYEHLyo1B0+jX8mHdVcv03tbwTOPpwurnFIdbcXt6jOQm5IZrWO31nwT9RoGzCfa81eHls9PPBmDHqcAC5rkT+nPZ2AycboUf/eMrma4RPUaCEdG3ndxJ4b/v5e8J8QD86Lxfti8QDcINI9QU+fr89Z5dabvD4XZQWKN+mA1JHOL9wBn3gt6kdt8s2NrleY+V1/Pu1+py4QUXDO3eII99pt7nTGudxugqM9LLlvVqtO/P314jp+pP0exqxqgzkhWXFc6htdLpYII6ZciuNR1enpg+EYADlQeo93MZtnpaHb++9ixSYs3sL67lD//NC8VoheiSBCOi76oqgL/fxPIEKw6djmlZ5zKj/4xun3Zs6lhGp4zG7rbzr8OdL68MxInaEzgVJzajzffW84kDYeLN6uONXRffrtx4BLcCF4zox5gg9qHxcLugoCUYGTTdE1x0VTOi1Vz0tmBECyiK6ov8Ov5YrVovMjB+YM/Wxkz4DuiNUPg5FLVfOZMZm0n/2P64Fbffq2o0aXEWfnWNuiLnmQ2H2HJIVteI8JNgRPRNTdWw5kZ2OSp5Py4WHTp+NOXB4Je1tqHT6bhhpJod+efBf3a7kFWrF8lNyO166/nzHgCdAQ6thxPeU/AF5Q38Y5uabVl4YTezIkW71P1QLAmQNcHvzMipejWzkBmb2b3r9zCtiPhErX/7tRyr7sGVNG3FZ8LolmW+2/7c4eVJGWrdSDC1TVeelcWNUwaiKHDvq7soqZVmaCK8JBgRfY/TDmtvQyndx1P91A/Oa4Zdw8iUkSG7xJWDr8RmtJFfk883Zd9061xd1ou0lTK4dSWFj+zIsv8exOlWmDk8janBNjnTHP5IvR98PhiMZLWspimqK8KteK8pKG4oBiAjJqN71+9h/ePUfXtONZzyq5/MkWq1LqfH6kXa+lZLh96vX4OG9hs5Ts6YDMDnRZ8Hdepfzh3HqMx4yuqaue+VXbhkZ18RRhKMiL7F7YZ37oWjn/BhQjK7THpsRhs/nPjDkF4m1hTLJYMuAeDtw29361xdrqQ53fkPAjo48G6H9DzAwVO1vLlT/av+octCEIAd/li9H3oxABmxGeh1euxuu8/GYMX1ajDS2zIjabY0THoTLsXlV93I4Sp1GezQpG5moIIx6FzIOEtd6bTzb+1empY1DVCbn9XZ6wI+tc1s4JlbzibWbOCzI+XSnVWElQQjou9QFHjvJ/DVKzh0Bn7fPxeA7439Hhmxof/r/OqhVwPwXv572F32oM9zqPIQEEAwkjYcxl2nPu5kZc2T7x9AUeDysZmMH5gU9LgAaK5VaxIAhl4EgElv8mQ7Ttaf7PRtTrfTs9qmt2VG9Dq9Jzui7Rfky6Eq9b/fsKRhYR1Xp3Q6+NYP1MefP6dmBVsMjB9Idnw2LsXFtuJtQZ1+aL84fvdtdbuEFzcf9XTyFSLUJBgRfcdH/w++eB7Q8er02ylsLifNlsYdY+8Iy+WmZk4l3ZZOdXM1m45vCuocTreTvCp1xcKolFH+v/H8h9T7fW/Dqb2epzfnlfHB3lMY9Dp+NHtEUGNqJ/9TcDshORdSWlf6aFM13j6syxrLcCtujDpjr+szAmpbfoATdb7rRqqaqjz9VCKSGQEYfxPEZULNCfjqlXYvTe8/HYAtJzvfVM8fc8Zn8aNZ6r+lR9/ew8aDgW0iKIQ/JBgRfcOmp9QbUHXZr3i2TF36+sOJPyTG1I3+Gj4Y9AbmDFU3LQt2qqagpoBmVzM2o611t15/pI+GMf+jPt6k1o44XG4ee0fdqfW2c3MYnhEf1JjaOfieet8yRaPp6sNam6JJj0nvuig3CmmZka6CES0r0j+2f9j+nXXJZFWXfQNsfrrdfjXaVM3Woq6Xgvty98XDuG7SAFxuhQV/286OgspunU+I0/W+3xJCtKUo8PFS+O8S9etLf8mz+lpq7bWMSB7BtcOuDevlrx6iTtVsPLGRyqbAf0Hvr9gPwIjkEYF/aGvZkd1vQOlBXtqSz6GSOlJizTxwaQiyIm43HFinPtZ2im2hfVgX1XW+/FUrXu1t9SIaLdjqapomovUibU2+HWJSoTIfvmltxndO1jnodXrya/K7DKx80el0LL3+LGYMS6XB7uJ7q75g9wnp0CpCR4IR0XspCqz/BXzyuPr1JY9w5Ky5rD2wFoCHznmo2w3OujI8eTijU0bjdDt5L/+9gN9/oPIAEOAUjSbzrJYgQaH+v7/l9+sPAvDjy0aSGGMK/HynO/El1J1Sl/QOPr/dS57MSH3nH3Dast7eVi+i0b4/rbOqNxGtF2nLHAvT71Eff/xrcDQCkGBO8LSG/7jg425dwmI08Od5UzgnN5naJie3vfi5tIwXISPBiOidXE545z7Y8kf168sfh5k/4qkvn8KluLhw4IWcm3Vujwzlf4ap0yVvHwp8quZAhRqMBL3s+AI1O2Lb/wZpjhOck5vMjVOC23eng33vqPfDZ6ubs7WhNQbzljnQPsQHxney8V8voC3T1VY6eXO4OkoyIwDfugsSBkB1IWx91vP0pYMuBeDDgg+7fYkYs5FVt5/DhIGJVDY4+M7KrWw94n1FlRD+kmBE9D5NNfD3G2HHS4AOrl4O5/4fHxd8zMbjGzHqjCyasqjHhnPF4Csw6ozsLt/tWRnjD0VR2FOu1niMSg4iMwLQfxLH+81Ej5v7TW/xxLcnoNd3v7EbigL7/60+HjWnw8sDYlunMTrrNaIFIwHVwUSR3MRcdOioaq6ioqmi02MURWF/uTrNFlRmK9TMMXDJI+rjTU9DndqU7uJBar3PzpKdPpdi+yveauLlO6cydXAKtc1O5q36gvd2+9etVghvJBgRvUvlMVh1GRz+L5hi4DtrYPLt1Dvq+fXnvwZg3th5/i+TDYEUawoXZF8AwFuH3vL7fcdqjlHVXIXFYAn6w2x/cQ2LimcD8D/6zeS6QrT08viXUHFE/RkPn9Xh5ay4LEx6E82uZk+xalsFNQVA7w1GbEabpy7mSNWRTo8prC2k1lGLWW9mSJKPPYV60lk3QtZEsNfCez8F1PqeMaljcCtuNhRuCMllEqwmXvr+t5g1JgO7082Cv+3g9+sP4pbGaCJIEoyI3mP/Olg5E0r2QlwG3LHO81f7n3b+iVMNpxgQN4AFExb0+NC0Qtl3jryDw+Xw6z27SncB6l43JkPgNR61TQ7+7287+MIxlC9sanaEDx8L+Dyd0paIjr4aLB1X5Rj1Rk/7c60Dqcbhdnj2demtwQi09n05/fvT7C1Xl1SPTBmJSR+CGp1Q0Ovhqt+DTg+7X1f/P0PrVM26o+tCdimrycCzt5zN96ap/w6W/zeP/335S6ob/Pv3L0RbEoyI6Oe0wwc/h1e/q+45M2Ay/OAjdRt11PTz3/f/HYBfnPsLbEZbjw9xxoAZpNnSqGiqYOPxjX69Z1fJLqB1h9VAuNwKD772FUfL6slKtDLi5ifUTdPy3oej/l3fK2ez+kEGrTvDdkL7sD69rqKorgiX4sJqsNIvpl/3xhJB2g7K3upGtGBkTOqYHhuTXwac3VrM+u8HoL6cq4ZchQ4dXxR/0WVRbiCMBj2//J9xPPHt8ZiNej7cV8JlyzayKU96kYjASDAiotuJHfD8Ba2FqucuhDveU3ewBWrttSzetBi34ubqIVczY0D3d+UNhlFvZO5QddOyNw+96dd7dpbsBGBiv4kBX+/X7+7j/T2nMBv0/Onms0nKHq0u7wT44Bfqstxg7XtH3Rgvvn+HVTRtaUWep2cOCmrVKZqB8QN7ZY8RjRaMaCtmThe1wQjAhYshdTjUFcMbd5IVk+4p6P7Xoe7vNH26G6Zk888F0xicFktxTRO3vbiVh97YxP6yYxTXF/udLRRnrh7c71qIADTXqq3Ot/wJFBfEpMHVy9RpgxaKovCrrb/iRN0JBsQN4GdTfxa58aJuxrdq9yo2ndhESUOJZ3fbzhTXF3Ok+gh6nd6z9NJfKz85zKpP1b/Wn7xxApNzktUXLvgJfPWqusvurr/B2fMC/yYUBT57Rn08+XbwsTRaq5M4PXOg9d7oybqdcBiVqtbx7C3fi6Io7XZ8drqdng0Sx6aOjcj4fDLZ4MaX4M+XqBsdbnica4dfy2dFn/HWobe4a8JdIZ1aUhQFxXKM2TO/4N2Dm6lyFvBerYv33lVf1+v0ZMZkMr7feM7OOJsLB17oWZElBEhmREQbtwt2/BX+OBk+Xa4GIuOuhx9+3i4QAfjr3r+y7ug6DDoDj898nDhzXIQGrRqcOJhJ6ZNwK27eyHvD57Fae+5xaeNIsib5fY0VGw6x9D/qCo6fXjGKuRP6t74Ylw4XqkWLfPALqAsiVV74BZzcAQYLTPm+z0O1YONw1WEUpbVw8WCl2u9kRHIIGq9F0IikEZj0JmrsNR2mNg5UHKDB2UC8OT7yPUa8yRir1o8AbPwdF5edJMWawqmGU7x3NPCeOJ1xuV38+8i/ue7t67h13a38ff/LVLuPotO7QNGhuE0oih634uZk/Uney3+P33z+G2a/Pptb1t3Cmn1rqLHXhGQsoneTYEREB5cTvloLz06Ht+9Rm22lDIHvvgrfXgWxae0O33h8I09vfxqAB6c8GFTdRTjcNPImAF7Z/wrNrmavx20+sRmA8wac59d5XW6Fpf/Zx+/eU/uSPHDpCO46v5MVHFP/T22G1lQF7weRKdI23ht/A8T5rvcYljQMo95IVXNVuw3z+kowYjKYGJ0yGoDdZbvbvfblqS8BODv97LA31uuWid+FGfcDYFn3ILelqFm4VbtXtQsgg/FF0Rd8+51vs3jTYg5VHcJisDBnyBx+O/O3/Oe6//D5zV9y39B/wtHfUnfwYRqO/YD4xqvIiRmHDh1fl37N4188ziX/uIRHPn2EPWV7uvvdil5MghERWQ0VaoOmP02GN/8XSveDJRFm/xoWfg4jr+jwlq1FW1m0YRFuxc31w6/nltG3RGDgnZudO5uMmAwqmip498i7nR7T4GjwBCPnD/Bek6GpaXLwv3/9kpWfqLUZP758JPddOrzdtIGHwaj2XUEH3/wD9vhXvwJA/mY4tF4thD2v6z4tFoOFkclqszZtysLhcnimaXp7MAIwNk2dgtG+P82XxWowMiVjSo+PKWCXPgaTbgPFzY2fvUSs3sShqkNBdQwGqHfU86utv2L+B/M5VHWIeFM89066l49u/IjHZz7OlUOuZGD8QGLMZn5w/hA2/+Ri7r5gEjbXSE7mn8fu7beiP/4Ik+NuZ1DcEJpcTbx56E2+8+53+M6/v8ObeW/S6GwM7c9ARD0JRkTPczZD3np4/U54apTaD6EyX91b4+JfwP1fw/S7O3T9BNhQuIF7/nsPza5mLhx4IQ9PfbjzD+UIMelN3DbmNgD+svsvON3ODsf8t+C/NDobGRQ/qMvix815ZVz++438d38JFqOe5d+ZyMILu5gWGDAZzntAffz2vVDhu4soAC4HvLdYfXz2PEj1r6PouLRxAOwuVTMHe8r34HA7SLIkefp09GZnp58NtN9ortHZ6Pl6atbUiIwrILqWxoBTvk+C283t5er03VNfPkmDoyGgU3128jOu+9d1ni0XbhhxA/+5/j/8YPwPSDAndPqepBgzP5o9ks0/uYgHZ4+gf6KVqlobG7aNYs+2H5Bas4hhMRdg1JnYU76HR7Y8wiX/uITHv3jca48X0fdIAavoGdUn4Nin6i6weeuhuc08ceZZcPb3YOItahfJFm63Qr3dSU2Tk6rGJv5x8K+8mf8iCgpDYqcwXLeQJz84RG2Tk5pGB7VNTpqdrk4vbzLoiTEbsJkM2MxGbCYDcRYDKbFmUuMspGr3cWZSYszd6mJ6/fDreeGbF8ivyeeNvDe4ceSN7V7/58F/AjBnyByvgVRhRQNPfnCAf+1Spz8GpcTwx+9OYkJ2kn+DuOhn6s+78HN49Wa1J4st2fvxm5dB8dfqMRcu9u8awPh+41l7YC3bT20HYEfJDgAmpU/q1StpNNP6T0Ov03Oo6hDF9cVkxmay9eRWmlxN9I/tHx2dV/2hN8CcpyGhP7d//BveiovjREMJT3y0iEcve67Lt9fZ63jyyyd5PU9d8j0gbgC/nP7LgIKxpBgzd188nAUXDOWj/SWs3VbIxrxS8k+kk3/iCnSGmSSm78Kc9AW1jhLW7FvDmn1rOCfzHK4ddi0XZl9IvDkEO1GLqBRUMLJixQqeeOIJioqKGDt2LMuWLWPmzJlej//kk09YtGgRe/bsoX///vz4xz9mwYKeb0wlekhDBcqpPThOfoPr5FcYCz/DVNO+M2i9OY285AvYmnQFB3TDqN3npGbHV9Q0qUFFTZODumYnigJ6ayHWjHcwxKhLRu2V5/LVvqv5ivD81WQ26hmQZGNgsq31PtlGdnIM2Skx9Iuz+AxW4sxxLJiwgMe/eJxndj3DpTmXkmJNAWDHqR3sKNmBUW/k+uHXd3jv7hPVrN6Sz9u7TmJ3qctzbz13EIuvGE2sJYD/uxpMcP2L8MKlapO4v38Hbl4LtqSOxx78ADb8Rn18+eNqIayftC3qd5fvpqyxjM9OfgbA5IzJ/o81iiVaEhmfNp5dpbv48NiH3DrmVt45ou7Zc/Ggi6MqK9clnQ7OfwjrwHN45J27WJCo8M/iTxnz8ixuOO8XkDtTPeY0m09s5pef/dLTafe7o77L/WffT4wppsOx/jAa9Mwem8nssZlUNzj4YG8x/9ldzGeHDVQVnQdF0zHEHsKUvBVj3D62FW9jW/E29BgZkzSZS3MvYdbgGWTHZ/eun7/wSacEWMW0du1abrvtNlasWMGMGTNYuXIlL7zwAnv37mXQoEEdjj969Cjjxo3jBz/4AXfddReffvopCxcu5JVXXuH66zv+Mu5MTU0NiYmJVFdXk5DQeSpQBEZRFBqdjVQ3V1PeWE15YxWVjdVUNNZQ1VxNdXMNNfYaau211DvqaHI20uRqwu6sw+5qxOFuQlFc6FAwKG4Mihuz20WS20mK206Cy02S202m00WG00k/p5saRxZ5jjH813kOO5VhKL5mCfVNGOP2Y0rcjjEuDwCd20pq8w30N55PgtVEvNVIvNXU5rERq8nQ4fepooDD5abB7qLR7qLR4aLB7qKu2UFFvZ2yOjsV9XbK65qp9KN7pNmoZ6AnOLExMDnG8zg7OYakGBNOxckNb9/A4erDzOg/gz9e/EecipNb193KwcqDXD/8eh6b/hjVDQ6+PlHFlsPlfLCnmMOl9Z7rzBiWyuIrRjNuQGJQ/40BKN4Nf7kSmqshbQRcu1JtiqX9YL7+h1ow7GpW6wrm/rHTDyRfvvvv77K7fDc/OOsHrNq9CpfiYt2168hO6L3dV9t6df+r/PrzX5ObkMufLvkT//PW/+BSXLw+9/XeWxfTWMkz//4+zzWoPVTurqxivj4V4+i5MOQiGHA2e+qO88yuZ9h0YhMAA+MGsmTGEs7JPKfzcyqKuiS/oRwaK6GxAhq0+4o2z1WomVFHEzhbb4qjEcXlwOUGtwJuRaHYYOCdeBsfxlk5am4fjKc5YaTdwGCXmcHYGKSLpZ8pDqs5FqPZhskag9kag8kSg9ESg8FsA6NVXfZstILJqt4bbW0et3ndaFU72opu8ffzO+BgZOrUqZx99tk8+2zrrpCjR4/mmmuuYenSpR2O/8lPfsLbb7/Nvn37PM8tWLCAr776is8++6zTazQ3N9Pc3LoSoaamhuzs7JAHI8v+cQ/51ftQAIX2P4bOvlbUB9r/oOhO+9rzqP1z7c94+vPtr6WeTgEd7Y5r+59JvW7ba7Y5p9I6VicKdp2CXefGoVMfO3QKdr362B2hPyqMxGLRJRJjSCLOmEKMyYbJYMCoV2h0V1FpL6aoocDzfel1eq4achX3TLqHzNjMsI7N4XJTXN3E8cpGjlc2cKKq0fP4eGUjRdVNuLrYfyPOYmRgso34hFIOGn6NGweJhhwUBWrcxzASR3b9I5RWmymqbmr3XpNBxxXjsvje9FzOHpQUmr/8ir+BNTdC7UlAB7nnQcpgKPpKvQGMugq+/ZdO63S68sr+V/jN57/xfD2+33jWXLmm++OOEvWOema9NotaRy0mvQmH28H0/tNZOWtlpIfWLYqi8LtNP+NvR9UNETOcTqY2qv8e91gsHDarfUiMwHfN/bk7biQxelNL8NAMjkZorGoJMMrVIMMdvuZmh01G1sfGsMVm5RuLBaeX/2+kuFxkOF3Eud3EKAoxbjc2RUGngB4FHXhunt+v7W46FJ362IUel06PGwNunR6XzoAbHYpOp963HOtG53leafnlraBv+T2uXk3nuVe1/arjfcf/bX+vPT796+65Ytz3uWzazd08S3thCUbsdjsxMTG89tprXHvttZ7n77vvPnbt2sUnn3zS4T3nn38+kyZNYvny5Z7n3nzzTW688UYaGhowmTo23nnsscf45S9/2eH5UAcjdz4/jc8tdSE7X29kVBQS3G4SXG4S3G7i3ep928exbgVby/+xbYqCHit6fSwuYzzNlkTstn7YYzKwx6XSEJNAU2wCdoNCg7uKyqZSTjWc4lTDKYrri2lyNXU9qBY5CTlclnsZ1wy9Jmr+ytaClcKKBgorGyisaGy5V4OVktr2y3kNsXnYBryMzmAHQHFZaCy8A1djrueYnNQYJmYncfGodC4cmU6iLQz7nNSVwPsPqyts2g3QAjN/BOc/6LPBmS+NzkZufOdG8mvy0ev0PD/r+d5R2BmAtw+/zcObHwYg3hzPmivX9PqmbqAGJP8+8m+e3PYEFc2V7V4zKQqz6xv4v8pqcpwdC7G9MtogJkWtP7IltzxOaf/Ymth5ZkJvVDNzStswAfVrlwOcTTgdTRwvL+WL0r3sqcknv6mI485SyqnBpetG52HBgvjL+eF1T4T0nP4GIwHVjJSVleFyucjIyGj3fEZGBsXFHXfuBCguLu70eKfTSVlZGVlZHbvwLV68mEWLWpcWapmRUDs7bSapVWrG5vQIE9pkqxXt606iU11rnKtrOabNGdr8b5t36TqLZnUtpzotDta1fZ/ecz71GT16nXpNfctxep0enQ70Oh0mnQGb3ohVb8KqN2HTm7DpjdgMZmx6I4lGK7FGM2aDHqPBgNGgR2eytaQptV8Qbe6tieovFENwdc+KolBjr6GssazdrdHZiIKCQWcg1ZpKRmwGo1NGk2pLDeo64WQy6MlOUWtHOtPkcHG8Ug1QKurs1NvHUtpwHntrNqLXwci4maSNTCMjwUpWko3c1BiSYgLPRgQsLh2u/zNc/HN1+W59mdpSf8TlHXq4BMpmtLH68tW8e+RdxqSOYUpmL1juGqC5Q+eSFZvFN2XfMCtnVq/eALAtnU7H1UOvZnbubLac2MKhqkMoKAy2pPItUwqJDZVqIOtsVKdVFLf6+8BgAaOlJcBIVYOMmFQ10DAHV0viLyOQOwhymdPueUVRqGquoqShhJKGEuod9dQ56qhprqfO3kCz04XT7cbhcuNwuXC4XDjd6rekoMY7CqC4FHRuF7id4Hahczs9N4PiBkVR8yCKgk5R0KPedIobnYJ6j6KeGHebey2PrbTJbHty7i33Wga+7XPqow7HdMi2n57PD9zogZH7IyKoT5XTU8ent0r25/jOntdYLBYsFkswQwvIwut+F/ZriFY6nY5ESyKJlkSGJvm3dLS3sZoMDEuPY1h6226wuUCUfEAn58A5d4b8tKm2VOaNDaL9fC9yTuY53uslejmLwcJFgy7iokEXRXooQdPpdCRbk0m2JjMyZWSkhyMCFFB1TlpaGgaDoUMWpKSkpEP2Q5OZmdnp8UajkdTU6PvLVwghhBA9K6BgxGw2M3nyZNavX9/u+fXr1zN9+vRO3zNt2rQOx3/wwQdMmTKl03oRIYQQQpxZAl63tGjRIl544QVWrVrFvn37eOCBBygoKPD0DVm8eDHz5rWmaxcsWMCxY8dYtGgR+/btY9WqVbz44os8+OCDofsuhBBCCNFrBVwzctNNN1FeXs6SJUsoKipi3LhxrFu3jpycHACKioooKCjwHD948GDWrVvHAw88wDPPPEP//v35wx/+4HePESGEEEL0bQH3GYkEaXomhBBC9D7+fn5LezkhhBBCRJQEI0IIIYSIKAlGhBBCCBFREowIIYQQIqIkGBFCCCFEREkwIoQQQoiIkmBECCGEEBElwYgQQgghIiq4veB7mNaXraamJsIjEUIIIYS/tM/trvqr9opgpLa2FoDs7OwIj0QIIYQQgaqtrSUxMdHr672iHbzb7ebkyZPEx8ej0+kiPZyIqqmpITs7m8LCQmmNH2bys+4Z8nPuGfJz7hnyc25PURRqa2vp378/er33ypBekRnR6/UMHDgw0sOIKgkJCfIPvYfIz7pnyM+5Z8jPuWfIz7mVr4yIRgpYhRBCCBFREowIIYQQIqIkGOllLBYLjz76KBaLJdJD6fPkZ90z5OfcM+Tn3DPk5xycXlHAKoQQQoi+SzIjQgghhIgoCUaEEEIIEVESjAghhBAioiQYEUIIIURESTAihBBCiIiSYKSPaG5uZuLEieh0Onbt2hXp4fQp+fn5zJ8/n8GDB2Oz2Rg6dCiPPvoodrs90kPr9VasWMHgwYOxWq1MnjyZTZs2RXpIfc7SpUs555xziI+PJz09nWuuuYYD/7+9ewmF9Y/jOP4ZNm6Toowkl6RcJompSSFKgyyYhSgpilJDNAvCRlIWiFJuG0QuIaQoNsZCymWklFPIGRk0Y+G2MBnzX53pPzmdv07O/3eex+dVFs93LN5Nevr2zA/fvonOkrXOzk4oFAo0NDSITpEMLiMy0djYiLCwMNEZsnRycoK3tzcMDw/j+PgYvb29GBoaQktLi+g0SZudnUVDQwNaW1thNpuRkZGB/Px8WCwW0WmyYjKZYDAYsLOzg42NDby+vkKn0+H5+Vl0mizt7u5iZGQESUlJolMkhX9nRAbW1tZgNBqxsLCAxMREmM1mJCcni86Sta6uLgwODuL8/Fx0imRptVqkpKRgcHDQPYuPj0dRURE6OzsFlsmbzWZDSEgITCYTMjMzRefIytPTE1JSUjAwMICOjg4kJyejr69PdJYk8MmIxN3e3qK6uhoTExPw8/MTnfNl3N/fIygoSHSGZDkcDuzv70On03nMdTodtre3BVV9Dff39wDAn98/wGAwoKCgADk5OaJTJEcS/7WXfs7lcqGiogI1NTXQaDS4uLgQnfQlnJ2dob+/Hz09PaJTJMtut8PpdEKlUnnMVSoVbm5uBFXJn8vlgtFoRHp6OtRqtegcWZmZmcHBwQF2d3dFp0gSn4z8hdra2qBQKH75tbe3h/7+fjw8PKC5uVl0siR99H3+N6vViry8PBQXF6OqqkpQuXwoFAqPa5fL9W5Gn6e2thZHR0eYnp4WnSIrl5eXqK+vx+TkJHx8fETnSBLPjPyF7HY77Hb7L78nKioKpaWlWFlZ8bh5O51OeHt7o6ysDOPj4386VdI++j7/uLlYrVZkZ2dDq9VibGwMXl7c5X+Xw+GAn58f5ubmoNfr3fP6+nocHh7CZDIJrJOnuro6LC0tYWtrC9HR0aJzZGVpaQl6vR7e3t7umdPphEKhgJeXF15eXjxeo/e4jEiYxWLBw8OD+9pqtSI3Nxfz8/PQarUIDw8XWCcvV1dXyM7ORmpqKiYnJ3lj+QRarRapqakYGBhwzxISElBYWMgDrJ/I5XKhrq4Oi4uL2NzcRGxsrOgk2Xl8fMT37989ZpWVlYiLi0NTUxM/EvsAnhmRsIiICI/rgIAAAEBMTAwXkU9ktVqRlZWFiIgIdHd3w2azuV8LDQ0VWCZtRqMR5eXl0Gg0SEtLw8jICCwWC2pqakSnyYrBYMDU1BSWl5ehVCrdZ3ICAwPh6+sruE4elErlu4XD398fwcHBXEQ+iMsI0X9YX1/H6ekpTk9P3y15fLD4+0pKSnB3d4f29nZcX19DrVZjdXUVkZGRotNk5cevTmdlZXnMR0dHUVFR8f8HEf0EP6YhIiIioXgCj4iIiITiMkJERERCcRkhIiIiobiMEBERkVBcRoiIiEgoLiNEREQkFJcRIiIiEorLCBEREQnFZYSIiIiE4jJCREREQnEZISIiIqH+AVICkTZVx+VLAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def f(x,alpha):\n", " return sps.norm.pdf(x)*(1+np.sin(alpha*x))\n", "\n", "plt.figure()\n", "x=np.linspace(-5,5,int(1e4))\n", "for alpha in [1,2,5]:\n", " y=f(x,alpha)\n", " plt.plot(x,y,label='alpha = %s'%alpha)\n", "plt.title(\"Density\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Question 2** : Simuler un échantillon de taille $n=1000$ de vaiid de loi de densité $\\pi_\\al$ via la méthode de Metropolis Hasting avec une loi de proposition $Q(x,.)=\\cN(x,\\si^2)$ (pour $\\si=1$ par exemple). On mesurera le temps T que prend la simulation avec la commande `time()` de la librairie `time`." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Estimation de la moyenne : 0.17483162236420008\n", "\n", "Erreur sur l'estimation de la moyenne : 0.06460838228941505\n", "\n", "Temps nécessaire pour tirer 1000 echantillons: 1.865004062652588\n" ] } ], "source": [ "Nsteps = 20 # nombre d'itérations de l'algorithme\n", "# en pratique il faudrait plus d'échantillons pour bien converger, mais cela prendrait trop de temps\n", "# on prend donc Nsteps=20 ici et on initialise HM avec une loi initiale gaussienne ce qui permet de vite converger car elle ressemble à f\n", "sigma = 1\n", "alpha = 2\n", "\n", "# On echantillonne 1000 fois avec la fonction HM (à écrire)\n", "n=int(1e3)\n", "Ech=[] \n", "t = time() \n", "for i in range(n):\n", " Ech.append(HM(Nsteps,sigma,alpha))\n", "tempsHM=time()-t\n", "\n", "# Estimation de l'espérance de la loi, \n", "mHM=np.mean(Ech)\n", "eHM=1.96*np.std(Ech)/np.sqrt(n)\n", "print(\"\\nEstimation de la moyenne : \"+str(mHM))\n", "print(\"\\nErreur sur l'estimation de la moyenne : \"+str(eHM))\n", "print(\"\\nTemps nécessaire pour tirer 1000 echantillons: \"+str(tempsHM))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGxCAYAAACwbLZkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABM6ElEQVR4nO3de1xUdf4/8NcADowgF0FQRwQUL+TkDdRQyVKj8LJdtqLaNAvajNKQzDK2Tf22S61puLZQ/tTMzJbS7tEFd1NRa1cJ1sxLoSI4oAgiCCIozO+PDzMwMiAzzMyZy+v5eMzjHA7nzHlDhC8+53ORaTQaDYiIiIgk4iJ1AUREROTcGEaIiIhIUgwjREREJCmGESIiIpIUwwgRERFJimGEiIiIJMUwQkRERJJiGCEiIiJJMYwQERGRpBhGiBzIpk2bIJPJcODAAYOfnzVrFkJDQ/WOhYaGQiaT4ZZbbjF4zebNmyGTySCTybBz507zFmwjtN+3oqIiqUshckoMI0SEXr16Yffu3Th+/Hi7z23cuBHe3t4SVEVEzoJhhIgwefJkKJVKbNy4Ue/48ePHsXv3bsTHx0tUGRE5A4YRIoKLiwvmzp2Ld999F83NzbrjGzduRHBwMKZPn96l97l06RIWL16MsLAweHh4oHfv3oiKisIHH3ygO+fAgQN44IEHEBoaCoVCgdDQUDz44IM4deqU3ntpH538+9//xuOPPw5/f394e3tj7ty5qKurw5kzZ3D//ffD19cX/fr1w+LFi3HlyhXd9UVFRZDJZPjb3/6Gv/zlLxg4cCA8PDwQFRWFf/3rX136enbs2IFp06bB29sbPXv2xKRJk7p8LRF1HcMIkQNqamrC1atX2706W6T7scceQ2lpKb799lvde7z77ruYN28eXFy69qsiJSUFmZmZWLhwIb755hu89957uO+++1BZWak7p6ioCMOGDUN6ejq+/fZbvPbaaygrK8O4ceNQUVHR7j0TExPh4+ODf/7zn/jTn/6ErVu34vHHH8fMmTMxatQobNu2DY888ghWrVqFtWvXtrv+zTffxDfffIP09HRs2bIFLi4uiIuLww8//NDp17JlyxbExsbC29sb7777Lj788EP07t0bt99+OwMJkblpiMhhvPPOOxoAnb5CQkL0rgkJCdHMnDlTo9FoNFOmTNHce++9Go1Go/nqq680MplMc/LkSc1HH32kAaD5/vvvO72/SqXS3HXXXUbVfPXqVU1tba3G09NTs2bNmnZfy4IFC/TOv+uuuzQANKtXr9Y7Pnr0aM3YsWN1H588eVIDQNO/f39NfX297nhNTY2md+/emunTp7e718mTJzUajUZTV1en6d27t2b27Nl692hqatKMGjVKM378eKO+RiLqHFtGiBzQ5s2bsX///navyZMnd3rdY489hs8//xyVlZXYsGEDbr311najbzozfvx4fP3113jhhRewc+dO1NfXtzuntrYWzz//PMLDw+Hm5gY3Nzd4eXmhrq4OR44caXf+rFmz9D6OiIgAAMycObPd8Wsf9QDAPffcAw8PD93HvXr1wuzZs7F79240NTUZ/Dr27duH8+fP45FHHtFrWWpubsYdd9yB/fv3o66u7vrfECLqEjepCyAi84uIiEBUVFS74z4+PigpKenwunvvvRcLFizAG2+8gS+++AKbNm0y6r5///vfMWDAAGRlZeG1116Dh4cHbr/9dqxcuRJDhgwBADz00EP417/+hZdeegnjxo2Dt7c3ZDIZZsyYYTC89O7dW+9juVze4fHLly+3u75v374GjzU2NqK2thY+Pj7tPn/27FkA4vvRkfPnz8PT07PDzxNR1zGMEJFOz5498cADDyAtLQ3e3t645557jLre09MTy5cvx/Lly3H27FldK8ns2bNx9OhRVFdX48svv8TLL7+MF154QXddQ0MDzp8/b+4vBwBw5swZg8fkcjm8vLwMXhMQEAAAWLt2LW666SaD5wQFBZmvSCInxzBCRHqefPJJnD17FlOmTNF7vGGsoKAgzJs3D//73/+Qnp6OS5cuQSaTQaPRwN3dXe/c9evXd/jIpLs+/vhjrFy5Uve1XLx4EV988QViYmLg6upq8JpJkybB19cXhw8fxtNPP22RuoioFcMIEekZPXo0Pv30U5OunTBhAmbNmoWRI0fCz88PR44cwXvvvYfo6Gj07NkTAHDzzTdj5cqVCAgIQGhoKHbt2oUNGzbA19fXfF9EG66urrjtttuQkpKC5uZmvPbaa6ipqcHy5cs7vMbLywtr167FI488gvPnz+Pee+9FYGAgzp07h//97384d+4cMjMzLVIvkTNiGCEis5k6dSo+//xzvPHGG7h06RKUSiXmzp2L1NRU3Tlbt27FM888gyVLluDq1auYNGkScnJy2nVINZenn34aly9fxsKFC1FeXo4RI0bgq6++wqRJkzq97uGHH8bAgQPxt7/9DU888QQuXryIwMBAjB49GvPmzbNIrUTOSqbRdDLxABGRnSoqKkJYWBhWrlyJxYsXS10OEXWCQ3uJiIhIUgwjREREJCk+piEiIiJJsWWEiIiIJMUwQkRERJJiGCEiIiJJ2cU8I83NzSgtLUWvXr0gk8mkLoeIiIi6QKPR4OLFi+jfvz9cXDpu/7CLMFJaWorg4GCpyyAiIiITlJSUYMCAAR1+3i7CSK9evQCIL8bb21viaoiIiKgrampqEBwcrPt3vCN2EUa0j2a8vb0ZRoiIiOzM9bpYsAMrERERSYphhIiIiCTFMEJERESSYhghIiIiSTGMEBERkaQYRoiIiEhSDCNEREQkKYYRIiIikhTDCBEREUmKYYSIiIgkxTBCREREkmIYISIiIkkxjBAREZGkTAojGRkZCAsLg4eHByIjI5Gbm9vp+e+//z5GjRqFnj17ol+/fnj00UdRWVlpUsFERF2hvlCPQ+pqHFJXQ32hXupyiKgTRoeRrKwsJCcnIzU1Ffn5+YiJiUFcXByKi4sNnr9nzx7MnTsXCQkJ+OWXX/DRRx9h//79SExM7HbxRESGqC/UY/qqXZi1dg9mrd2D6at2MZAQ2TCjw8jq1auRkJCAxMREREREID09HcHBwcjMzDR4/o8//ojQ0FAsXLgQYWFhmDx5Mp544gkcOHCg28UTERlSVdeI+itNSI8fjfT40ai/0oSqukapyyKiDhgVRhobG5GXl4fY2Fi947Gxsdi3b5/BayZOnIjTp08jOzsbGo0GZ8+exbZt2zBz5swO79PQ0ICamhq9FxGRscIDvRAe6CV1GUR0HUaFkYqKCjQ1NSEoKEjveFBQEM6cOWPwmokTJ+L9999HfHw85HI5+vbtC19fX6xdu7bD+6SlpcHHx0f3Cg4ONqZMIiIisiMmdWCVyWR6H2s0mnbHtA4fPoyFCxfiz3/+M/Ly8vDNN9/g5MmTmD9/fofvv3TpUlRXV+teJSUlppRJREREdsDNmJMDAgLg6urarhWkvLy8XWuJVlpaGiZNmoTnnnsOADBy5Eh4enoiJiYGr7zyCvr169fuGnd3d7i7uxtTGhEREdkpo1pG5HI5IiMjkZOTo3c8JycHEydONHjNpUuX4OKifxtXV1cAokWFiIiInJvRj2lSUlKwfv16bNy4EUeOHMGiRYtQXFyse+yydOlSzJ07V3f+7Nmz8fHHHyMzMxMnTpzA3r17sXDhQowfPx79+/c331dCREREdsmoxzQAEB8fj8rKSqxYsQJlZWVQqVTIzs5GSEgIAKCsrExvzpF58+bh4sWLePPNN/Hss8/C19cXU6dOxWuvvWa+r4KIiIjslkxjB89Kampq4OPjg+rqanh7e0tdDhHZuEPqasxauwdfLpgMALp9ldJH4sqInEtX//3m2jREREQkKYYRIiIikpTRfUaIiGyV+kI9quoaUVheK3UpRGQEhhEicgjaxfHqrzQBABQ9XOHnKeeaNER2gGGEiBxC28XxwgO94Ocph9JXwTBCZAcYRojIoYQHenHUDJGdYQdWIiIikhTDCBEREUmKYYSIiIgkxTBCREREkmIYISIiIkkxjBAREZGkGEaIiIhIUgwjREREJCmGESIiIpIUwwgRERFJimGEiIiIJMUwQkRERJJiGCEiIiJJMYwQERGRpBhGiIiISFIMI0RERCQphhEiIiKSFMMIERERSYphhIiIiCTFMEJERESSYhghIiIiSTGMEBERkaQYRoiIiEhSDCNEREQkKYYRIiIikhTDCBEREUmKYYSIiIgkZVIYycjIQFhYGDw8PBAZGYnc3NwOz503bx5kMlm714gRI0wumoiIiByH0WEkKysLycnJSE1NRX5+PmJiYhAXF4fi4mKD569ZswZlZWW6V0lJCXr37o377ruv28UTERGR/TM6jKxevRoJCQlITExEREQE0tPTERwcjMzMTIPn+/j4oG/fvrrXgQMHUFVVhUcffbTbxRMREZH9MyqMNDY2Ii8vD7GxsXrHY2NjsW/fvi69x4YNGzB9+nSEhIR0eE5DQwNqamr0XkREROSYjAojFRUVaGpqQlBQkN7xoKAgnDlz5rrXl5WV4euvv0ZiYmKn56WlpcHHx0f3Cg4ONqZMIiIisiMmdWCVyWR6H2s0mnbHDNm0aRN8fX1x1113dXre0qVLUV1drXuVlJSYUiYRERHZATdjTg4ICICrq2u7VpDy8vJ2rSXX0mg02LhxI+bMmQO5XN7pue7u7nB3dzemNCIiIrJTRrWMyOVyREZGIicnR+94Tk4OJk6c2Om1u3btQmFhIRISEoyvkoiIiByWUS0jAJCSkoI5c+YgKioK0dHRWLduHYqLizF//nwA4hGLWq3G5s2b9a7bsGEDJkyYAJVKZZ7KiYiIyCEYHUbi4+NRWVmJFStWoKysDCqVCtnZ2brRMWVlZe3mHKmursb27duxZs0a81RNREREDsPoMAIASUlJSEpKMvi5TZs2tTvm4+ODS5cumXIrIiIicnBcm4aIiIgkxTBCREREkmIYISIiIkkxjBAREZGkGEaIiIhIUgwjREREJCmGESIiIpIUwwgRERFJimGEiIiIJMUwQkRERJJiGCEiIiJJMYwQERGRpBhGiIiISFIMI0TkuOrrgd9+k7oKIroOhhEicjzNzcCrrwJ9+gAPPCCOpaYC589LWxcRGcQwQkSORaMBnnwSWLoUqKsDPD3F8W++AW69lYGEyAYxjBCRY/n0U2DdOsDFBXj7bWDXLnHc3x84eBB49FERWIjIZjCMEJFjef11sX31VeCPfwRkMvHx2rWAXA58/jmwZYt09RFROwwjRORYLl8GJk4Enn1W//iwYcCyZWJf+wiHiGwCwwgROYbTp1v3V64Uj2mutWgREBoKqNXA//t/ViuNiDrHMEJEjuHDD8U2Olq0jBji4QG8+KLYf+MN4MoV69RGRJ1iGCEi+1dXB3z2mdjXDuXtyJw5QGAgUFwMbN9u+dqI6LoYRojI/mVlAbW1Yr+jVhEtDw9g/nyxv2GDZesioi5hGCEi+7d1a+u+ob4i15o3T2z/9S/RQkJEkmIYISL7duYM8P33xl0TFiYmQNNogHfftUxdRNRlDCNEZN+2bRPTv6tUxl03d67Yaju+EpFk3KQugIioW7ZtE9vYWKCTmd4Ly2t1+36ecih/9zvAzQ04dEgspjdkiIULJaKOsGWEiOxXdTWwZ4/Yv2WKwVP8POVQ9HBFclYBZq3dg1lr92D6ql1QuyjEoxoA+OQTKxVMRIYwjBCR/dqxA2hqErOrKgcYPEXpq8COZ6fgywWT8eWCyUiPH436K02oqmsE7rlHnPTxx1YsmoiuxTBCRPYrO1ts4+I6PU3pq4BK6QOV0gfhgV6tn7jzTrH9z3+As2ctVCQRXQ/DCBHZJ40G+PprsT9jhmnv0a8fMGaM2M/JMU9dRGQ0hhEisk+HDgFlZUDPnsDNN5v+PrffLrbffWeeuojIaCaFkYyMDISFhcHDwwORkZHIzc3t9PyGhgakpqYiJCQE7u7uGDx4MDZu3GhSwUREAICdO8V20iTA3d3094mNFdvvvhOtLURkdUYP7c3KykJycjIyMjIwadIkvP3224iLi8Phw4cxcOBAg9fcf//9OHv2LDZs2IDw8HCUl5fj6tWr3S6eiJzY7t1iO8XwKJoumzhRtK6cPQscPAiMGtX92ojIKEa3jKxevRoJCQlITExEREQE0tPTERwcjMzMTIPnf/PNN9i1axeys7Mxffp0hIaGYvz48Zh4vfUjiIg6otGYL4y4uwO33CL2+aiGSBJGhZHGxkbk5eUhVtus2SI2Nhb79u0zeM3nn3+OqKgo/O1vf4NSqcTQoUOxePFi1NfXd3ifhoYG1NTU6L2IiHSOHQPKy8Wid+PGdf/9pk0TW23AISKrMuoxTUVFBZqamhAUFKR3PCgoCGfOnDF4zYkTJ7Bnzx54eHjgk08+QUVFBZKSknD+/PkO+42kpaVh+fLlxpRGRM5k1y6xvemm7vUX0YqJEds9e8TU8l1ZbI+IzMak/+NkMpnexxqNpt0xrebmZshkMrz//vsYP348ZsyYgdWrV2PTpk0dto4sXboU1dXVuldJSYkpZRKRo9KGke4+otEaMwbw9AQuXBCjdIjIqowKIwEBAXB1dW3XClJeXt6utUSrX79+UCqV8PHx0R2LiIiARqPB6dOnDV7j7u4Ob29vvRcRkc4PP4ittkWju9zcREdWgI9qiCRgVBiRy+WIjIxEzjWTA+Xk5HTYIXXSpEkoLS1FbW3rIlW//vorXFxcMGCA4embiYg6VF4OFBUBMhkQFWW+99XOVXKdqQqIyPyMfkyTkpKC9evXY+PGjThy5AgWLVqE4uJizJ8/H4B4xDJXuzQ3gIceegj+/v549NFHcfjwYezevRvPPfccHnvsMSgUCvN9JUTkHPbvF9vhw4E2La7dpm1lyc3lfCNEVmb0PCPx8fGorKzEihUrUFZWBpVKhezsbISEhAAAysrKUFxcrDvfy8sLOTk5WLBgAaKiouDv74/7778fr7zyivm+CiJyHv/9r9iOH2/e9x0/HpDLxayux48D4eHmfX8i6pDRYQQAkpKSkJSUZPBzmzZtands+PDh7R7tEBGZxFJhRKEAxo4FfvxRLJzHMEJkNRy/RkT2Q6OxXBgBgAkTxFZ7DyKyCpNaRoiIJHHiBHD+vHicMnIkAEB9oR5VdY0oLK+9zsVdoA04//lP99+LiLqMYYSI7Ie2xWLMGEAuh/pCPaav2oX6K00AAEUPV/h5yk1/f23LSH4+0NBgngnViOi6GEaIyH5oR9K0TAFfVdeI+itNSI8fjfBAL/h5yqH07cYovUGDAH9/oLJSLJpnjqnmiei62GeEiOxHQYHYjh2rdzg80AsqpU/3gggg5i7hoxoiq2MYISL7oNEA//uf2B81ynL30T6qYRghshqGESKyD6dPi86rbm7ADTdY7j7alhGOqCGyGoYRIrIP2laR4cMBDw/L3UcbRn79VSycR0QWxzBCRPZB219k9GjL3sffHxg4UP+eRGRRDCNEZB+s0V9ES9tBNj/f8vciIoYRIrIT1moZAcQ8JgDDCJGVMIwQke27eBEoLBT71mgZYRghsiqGESKyfT//LLb9+wN9+lj+ftowcuQIUF9v+fsROTmGESKyfdpHNGZsFSksr8UhdTXUFwyEDaUSCAgAmpqAQ4fMdk8iMoxhhIhsn7ZlpGVxvO7w85RD0cMVyVkFmLV2D6av2tU+kMhkfFRDZEUMI0Rk+w4fFtsRI7r9VkpfBXY8OwVfLpiM9PjRqL/ShKq6xvYnasPITz91+55E1DkulEdEts+MYQQQgeS669iwZYTIatgyQkS27dw5oKJCPDoZPtx699XONXLwIHD1qvXuS+SEGEaIyLb98ovYhoYCPXta777h4YCnJ3D5spganogshmGEiGybmR/RdJmLC6BSiX1tB1oisgiGESKybdow0malXvWFehxSV6OwvNay977xRrHl8F4ii2IHViKybdeEEfWFekxftQv1V5oAAIoervDzlFvm3towwpYRIotiGCEi26btM9ISRqrqGlF/pQnp8aMRHugFP0/59UfGmIqPaYisgmGEiGxXRQVQXi72IyL0PhUe6AWV0sey99e2jJw4AdTWAl5elr0fkZNinxEisl1HjohtSIg0QaBPHyAoSOxrHxcRkdkxjBCR7TLQedXq2G+EyOIYRojIZmhHyegWsNP2F7H2sN62GEaILI59RojIJhgaJbOjsARKoF1/EatiJ1Yii2PLCBHZhLajZHQL2KlbOq9acxr4a3GuESKLY8sIEdmU8MA2HVXPnhXboUOlKQYQj4hkMjGqp7wcCAyUrhYiB8WWESKybb17AwEB0t2/Z09g8GCxz0c1RBbBMEJEtk3KVhEtbb8RPqohsgiGESKybbYQRjiihsiiTAojGRkZCAsLg4eHByIjI5Gbm9vhuTt37oRMJmv3Onr0qMlFE5ETsYUwop3nRDsJGxGZldFhJCsrC8nJyUhNTUV+fj5iYmIQFxeH4uLiTq87duwYysrKdK8hQ4aYXDQRORFbCCPaocVHjgAajbS1EDkgo8PI6tWrkZCQgMTERERERCA9PR3BwcHIzMzs9LrAwED07dtX93J1dTW5aCJyIsOGSV2BCEQuLkBVVetaOURkNkaFkcbGRuTl5SE2NlbveGxsLPbt29fptWPGjEG/fv0wbdo0fP/9952e29DQgJqaGr0XETmR6urW/fBw6erQUiiAsDCxzzVqiMzOqDBSUVGBpqYmBGkXjmoRFBSEM2fOGLymX79+WLduHbZv346PP/4Yw4YNw7Rp07B79+4O75OWlgYfHx/dKzg42JgyicjenToltkFBYmitLWj7qIaIzMqkSc9kMpnexxqNpt0xrWHDhmFYm2bW6OholJSU4PXXX8fNN99s8JqlS5ciJSVF93FNTQ0DCZEzOXUKgJ9YrRdiqviqukYAQGF5rVVKaHtPP085lDfcAHz5JVtGiCzAqDASEBAAV1fXdq0g5eXl7VpLOnPTTTdhy5YtHX7e3d0d7u7uxpRGRI5EG0YGDmy3Zg0g1q3x85Rb7PYG18kJHyHWyWHLCJHZGRVG5HI5IiMjkZOTg7vvvlt3PCcnB3feeWeX3yc/Px/9+vUz5tZE5EyKi4Hg0UBoqN6aNdqp4v085VD6Kix2+7b3BIDkrAJUDRzMMEJkIUY/pklJScGcOXMQFRWF6OhorFu3DsXFxZg/fz4A8YhFrVZj8+bNAID09HSEhoZixIgRaGxsxJYtW7B9+3Zs377dvF8JETmOU6eAYAADB+oOhQd6QaX0sWoZeuvkaDuwlpUBFy4Avr5WrYXIkRkdRuLj41FZWYkVK1agrKwMKpUK2dnZCGl5tltWVqY350hjYyMWL14MtVoNhUKBESNG4KuvvsKMGTPM91UQkeNobhYtI4Cuz4hN8PQElEpArRatI9HRUldE5DBM6sCalJSEpKQkg5/btGmT3sdLlizBkiVLTLkNETmjs2eBRtFxFLb2OPeGG0QYOXyYYYTIjLg2DRHZFu2wXgCw0uSIheW1OKSuhvpCfecncngvkUWY1DJCRGQxJSVWu5WfpxyKHq5IzioA0DJq5tkpHV/AMEJkEQwjRGRbSkutdiulrwI7np2CqrpGFJbXilEzLXOLGKRdMI9zjRCZFcMIEdkWtdqqt1P6Kro+TFjbMnLqFFBXJzq1ElG3sc8IEdkWK4cRo/TpAwQEiJV7jx2Tuhoih8EwQkS2xZbDCMB+I0QWwDBCRLbl4kWpK+icNoyw3wiR2TCMEJHt6e0ndQUdGz5cbPmYhshsGEaIyPYoB0hdQce0q5AzjBCZDcMIEdme/v2lrqBj2paR334Dmpo6P5eIuoRhhIhszwAbbhkJCQHc3YGGhtY1dIioWxhGiMj22HLLiKsrEB4u9vmohsgsGEaIyPbYcssIwH4jRGbGMEJEtqG5uXXflltGAIYRIjNjGCEi23DunNi6uABBQdLWcj0MI0RmxTBCRLZBu0Bev36iX4Yt04aRo0elrYPIQTCMEJFtOH1abJVKaevoCm0YKS21/RljiewAwwgR2QZty4g9hBE/P7FoHgD8+qu0tRA5AIYRIrIN2gXy7CGMAOw3QmRGDCNEZBsYRoicFsMIEdkGbRix9WG9Wlwwj8hsGEaISHqXL7cO7bX1Cc+02DJCZDYMI0QkvaKi1n0fH8nKMIo2jPz6q/6EbURkNIYRIpLeyZOt+zKZdHUYIywMcHMDLl1qfcRERCZhGCEi6bUNI/aiRw9g8GCxz0c1RN3CMEJE0jtxQuoKrquwvBaH1NU4pK6G+kK9OMh+I0Rm4SZ1AUREttwy4ucph6KHK5KzCnTHFD1csePZKVAyjBCZBcMIEUnPhltGlL4K7Hh2CqrqGgGIFpLkrAJU1TUyjBCZCcMIEUlLoxFhRNFH6ko6pPRVQOmraP8JLphHZBbsM0JE0qqqAmpqpK7CNNowUlwsRtUQkUkYRohIWtr+Iv7+0tZhioAAoHdvsf/bb9LWQmTH+JiGiKSl7S9iL2vStCgsrxU7kVOAgweB/N/gFzLU8OMcIuoUwwgRSUvbMmInYaTd6JqxCcBYAEcBxapdYpQNAwmRUUx6TJORkYGwsDB4eHggMjISubm5Xbpu7969cHNzw+jRo025LRE5Il3LiH0skKcdXfPlgsni5fUbvtz0DNIr96H+SpNu1A0RdZ3RYSQrKwvJyclITU1Ffn4+YmJiEBcXh+Li4k6vq66uxty5czFt2jSTiyUiB6RrGbGTBfIgAolK6SNeI0KgOnsc4cd/lrosIrtldBhZvXo1EhISkJiYiIiICKSnpyM4OBiZmZmdXvfEE0/goYceQnR0tMnFEpED0raM9LePlpF2tCNq2i72R0RGMSqMNDY2Ii8vD7GxsXrHY2NjsW/fvg6ve+edd3D8+HG8/PLLXbpPQ0MDampq9F5E5ICamoBTp8T+APtpGdEzeDDg4sKhvUTdYFQYqaioQFNTE4KCgvSOBwUF4cyZMwav+e233/DCCy/g/fffh5tb1/rLpqWlwcfHR/cKDg42pkwishdqNXDlilh0rk/7Sc8Ky2tbR63YKnd3sYIvEZnMpA6ssmuW+NZoNO2OAUBTUxMeeughLF++HEOHDu3y+y9duhTV1dW6V0lJiSllEpGt0/YXCQkBXF11h9uOWEnOKoCihyv8POUSFdkF2kc1RGQSo4b2BgQEwNXVtV0rSHl5ebvWEgC4ePEiDhw4gPz8fDz99NMAgObmZmg0Gri5ueG7777D1KlT213n7u4Od3d3Y0ojInuk7S9yTcvCtevB+HnKbXu47NChQB7XpyEylVFhRC6XIzIyEjk5Obj77rt1x3NycnDnnXe2O9/b2xs//6zfwzwjIwP//ve/sW3bNoSxaZPIuWlbRgYNavepDteDsUXDhgH4QuoqiOyW0ZOepaSkYM6cOYiKikJ0dDTWrVuH4uJizJ8/H4B4xKJWq7F582a4uLhApVLpXR8YGAgPD492x4nICXXQMmJ3+JiGqFuMDiPx8fGorKzEihUrUFZWBpVKhezsbISEhAAAysrKrjvnCBERgE5bRuxK2zBy5Yp0dRDZKZOmg09KSkJSUpLBz23atKnTa5ctW4Zly5aZclsicjSO0jLSrx+gaHmkdPo0EBogbT1Edoar9hKRNOrrAW1neHtvGZHJgNBQsX+qSMpKiOwSwwgRSUM7Y6m3N+DnJ2kpZtHyqBpFp6Stg8gOMYwQkTTaPqIxME+R3QltCSOnGEaIjMUwQkTScJTOq1ohoWLLNWqIjMYwQkTScJTOq1ohbBkhMhXDCBFJw9FaRgYOFNvqaqCyUtpaiOwMwwgRWZ36Qj0OVVzGoaDBUPcLlboc81C0mS32GKeGJzKGSfOMEBGZSn2hHtNX7UL9pKeBSYDiJxl2TK2XuizzOnYMmDhR6iqI7AbDCBFZVVVdI+qvNCH9i9cBAMmzF+sWxHMYbBkhMgrDCBFJIryyBOjTR+oyLINhhMgo7DNCRNLp31/qCiyDYYTIKAwjRCSdAQOkrsAyCguBq1elroLIbjCMEJF0WlpGCstrUVheK3ExZiKXi5V7Od8IUZexzwgRScYvpD8Up1yRnFUAAFD0cIWfp1zaorpr4ECg5Ih4VDN4sNTVENkFhhEikoxyaAh23DtON5rGz1MOpa/iOlfZuJAQYC9EGJkxQ+pqiOwCwwgRWVdTU+v+oEFQ+irsP4C0FRoqtuzEStRl7DNCRNZ19qzY9ujhmKNptGvUMIwQdRnDCBFZV2mp2PbrB7g44K8ghhEio/ExDRFZ1+nTAAIApVLqSgzq9qge7WOasjKgpgbw9u52TUSOjmGEiKyrtBRAgM09ovHzlEPRwwwje7y8gKAg8Tjq11+BqCjzFkrkgBhGiMi6Tp8G+o20uQnPlL4K7Hh2inlG9gwbJsLIsWMMI0RdwDBCRNZVWgr0g821jAAw38ieoUOB3bvZb4Soixyw9xgR2TS1WmxtrGXErIYNE1uGEaIuYRghIuupqwPOnxf7NtgyYjYMI0RGYRghIuspKmrdd+RRJtow8uuvQHOztLUQ2QGGESKynhMnpK7AOsLCADc3oL6+ZSgzEXWGYYSIrOfkSakrsI4ePVoXyeOjGqLrYhghIutxlpYRQP9RDRF1imGEiKzHWVpGAHZiJTICwwgRWY8ztowwjBBdF8MIEVmHRsOWESIyiGGEiKzj3Dkxz4hMJnUl1qENI8XFYlQNEXXIpDCSkZGBsLAweHh4IDIyErm5uR2eu2fPHkyaNAn+/v5QKBQYPnw43njjDZMLJiI7pW0V6dNH2jqsJSAA8PMTLUK//SZ1NUQ2zegwkpWVheTkZKSmpiI/Px8xMTGIi4tDcXGxwfM9PT3x9NNPY/fu3Thy5Aj+9Kc/4U9/+hPWrVvX7eKJyI5o+4soldLWYS0yGR/VEHWR0WFk9erVSEhIQGJiIiIiIpCeno7g4GBkZmYaPH/MmDF48MEHMWLECISGhuLhhx/G7bff3mlrChE5IG3LiLOEEUAsmAcwjBBdh1FhpLGxEXl5eYiNjdU7Hhsbi3379nXpPfLz87Fv3z5MmTKlw3MaGhpQU1Oj9yIiO+dsLSMAW0aIusioMFJRUYGmpiYEBQXpHQ8KCsKZM2c6vXbAgAFwd3dHVFQUnnrqKSQmJnZ4blpaGnx8fHSv4OBgY8okIlukbRlx5AXyrsUwQtQlbqZcJLumN7xGo2l37Fq5ubmora3Fjz/+iBdeeAHh4eF48MEHDZ67dOlSpKSk6D6uqalhICGyM+oL9aiqa9R97HemCkoAGDAAOFkrWV3Wor5Qj6qAUCBoMPzUFVBqNM4zkojISEaFkYCAALi6urZrBSkvL2/XWnKtsLAwAMCNN96Is2fPYtmyZR2GEXd3d7i7uxtTGhHZEPWFekxftQv1V5p0xxR3pGJHyZMtLSOOPUW63tc/bw0UjZex47cSKIcOlLo0Iptk1GMauVyOyMhI5OTk6B3PycnBxIkTu/w+Go0GDQ0NxtyaiOxIVV0j6q80IT1+NL5cMBnpt/ZHfQ8PVPn4iyGvDq7t15/+47uol3ug6mih1GUR2SyjH9OkpKRgzpw5iIqKQnR0NNatW4fi4mLMnz8fgHjEolarsXnzZgDAP/7xDwwcOBDDhw8HIOYdef3117FgwQIzfhlEZIvCA72gUvoA/60QB/r3B1ycZ67F8EAvoLeH+ODUKWmLIbJhRoeR+Ph4VFZWYsWKFSgrK4NKpUJ2djZCQkIAAGVlZXpzjjQ3N2Pp0qU4efIk3NzcMHjwYLz66qt44oknzPdVEJFtU6sB9AEGONFIGq2W340MI0QdM6kDa1JSEpKSkgx+btOmTXofL1iwgK0gRM5OG0b6O2kYOQOGEaJOOE97KRFJp7QUAFDYdxAKyx1/JI2e0FCxLSqSsgoim8YwQkQW53fyVygaLyP5Yn8kZxVA0cMVfp5yqcuyjpaRhDh9Grh8WdpaiGyUSY9piIiMoTxSgB0FT6Lqs6+AoUPh5ymH0lchdVnW4e8vthoN8OuvwMiR0tZDZIMYRojIsmpqgIoKKAEoI4cB3t5SV2RdbSc6O3yYYYTIAD6mISLL0k4DHxDgfEHkWocPS10BkU1iGCEiy9IukDdokLR12AKGESKD+JiGiCyLYaTV4cPt1+xxpv4zRB1gGCEiy2IY0VGfqcL0VTtRf6VZd0zRwxU7np3CQEJOjWGEiCyLYURQKFAFoP5KM9LjRyM80AuF5bVIzipAVV0jwwg5NfYZISLLYhgR2nz92jV7wgO9JCyIyHYwjBCR5TQ1tc48yjAidQVENothhIgsp+Ic0NgIuLkBAwZIXY20GEaIOsQwQkSWc1ottqGhgKurpKVIjmGEqEMMI0RkOeqWMMJ/iFvXqAHE4ysi0mEYISLLUZ8WW4YRoF8/wN1d7J8+LW0tRDaGYYSILOc0W0Z0XFxaW0dOnpC2FiIbwzBCRJbDxzT6tN+HEyelrYPIxjCMEJHl8DGNPl0YYcsIUVsMI0RkOeerxJZhRNA+pmEYIdLDMEJEltW7N+DjI3UVtkEbyoqKgObmTk8lciYMI0RkWWwVaaVUim1DA3DqlLS1ENkQhhEisiyGkVZtJ347dEi6OohsDMMIEVkWw4hhDCNEOgwjRGRZDCOGHTwodQVENsNN6gKI7J36Qj2q6hoBAH6ecih9FRJXZGMYRgAAheW1+gd+/rnd5/jzQ86KYYSoG9QX6jF91S7UXxFrjSh6uGLHs1P4D0rbkSJOHkb8POVQ9HBFclYBAEDhJoPfpRrg2Cn4uUH/c/z5ISfFMELUDVV1jai/0oT0+NEAgOSsAlTVNfIfk4oKsXVxAYKDpa1FYkpfBXY8O6W19axnDyjfbASuXoXyTJHuc4Xltfz5IafFMEJkBuGBXlKXYFu0C8H16we48deM0lehHzBGjgRyc4Gff4Zy1CiGD3J67MBKROZXXCy2AwdKW4etuvFGsW3Tb4TImTGMEJH5lZSIrZM/oumQNoxwRA0RAIYRIrIEbRhhy4hhI0eKLVtGiAAwjBCRGakv1OOQuhqF5y+LA2wZMUylElu1GqiqkrYWIhvAnmVEZBZ6w5xH3Q9F42X4hQ+Tuizb5O0NhISI9Wl+/hm4+WapKyKSlEktIxkZGQgLC4OHhwciIyORm5vb4bkff/wxbrvtNvTp0wfe3t6Ijo7Gt99+a3LBRGSbdMOcbw/Fl5uewY6NSVDeOFTqsmyX9lEN+40QGR9GsrKykJycjNTUVOTn5yMmJgZxcXEo1vaev8bu3btx2223ITs7G3l5ebj11lsxe/Zs5Ofnd7t4IrI94XXnoDp7HMrenoBcLnU5tosjaoh0jA4jq1evRkJCAhITExEREYH09HQEBwcjMzPT4Pnp6elYsmQJxo0bhyFDhuCvf/0rhgwZgi+++KLbxRORDdJ2Xh0yRNo6bB3DCJGOUWGksbEReXl5iI2N1TseGxuLffv2dek9mpubcfHiRfTu3bvDcxoaGlBTU6P3IiI7oW0lDQ+Xtg5b13ZETdvp84mckFFhpKKiAk1NTQgKCtI7HhQUhDNnznTpPVatWoW6ujrcf//9HZ6TlpYGHx8f3SuYPfLJxuhGjVy7+Bm1tow4aRgpLK/t2s/FkCHiMVZtrejISuTETBpNI5PJ9D7WaDTtjhnywQcfYNmyZfjss88QGBjY4XlLly5FSkqK7uOamhoGErIZhhbH8/OU69YecXralhEne0zTbkG8lp+LDvXoAdxwA1BQIF7jp1qjTCKbZFQYCQgIgKura7tWkPLy8natJdfKyspCQkICPvroI0yfPr3Tc93d3eHu7m5MaURW03ZxvPBAL92y7wwjLZy0ZaTdgngtPxedGjtWBJH8fIYRcmpGPaaRy+WIjIxETk6O3vGcnBxMnDixw+s++OADzJs3D1u3bsXMmTNNq5TIxoQHekGl9OEiZ9eqrxer9YaFSV2J1Sl9FVApfbr+czFmjNj+9JNlCyOycUY/pklJScGcOXMQFRWF6OhorFu3DsXFxZg/fz4A8YhFrVZj8+bNAEQQmTt3LtasWYObbrpJ16qiUCjg4+Njxi+FiGzGwIEAWzevb+xYseVUB+TkjA4j8fHxqKysxIoVK1BWVgaVSoXs7GyEhIQAAMrKyvTmHHn77bdx9epVPPXUU3jqqad0xx955BFs2rSp+18BEdkeJ+svYrKRIwGZDCgtBSorpa6GSDImdWBNSkpCUlKSwc9dGzB27txpyi2IyE6oL9Sjqq5RfwSJk/UXMZmXFzBsGHD0KHDsmNTVEEmGa9MQkcnajSxqvgK/SzUMI8YYM0aEkaNHAXD6fHJOXLWXiEzWdmTRlwsmY8fev0N58Rwf0xhD22/k6FFp6yCSEMMIEXVbeKAXVP29ofz5QMsBtox0mXZEDcMIOTGGESIyj7Iy4OJFwNUVGDxY6mrshzaMqNXS1kEkIYYRIjIPbQfMQYO4Wq8xevcGWkYjEjkrdmAl6iKDo0YsdA+gizN42hLtY4Zhw6Stwx6NHQvsOyh1FUSSYRgh6oKO1qOx9D12PDvFfgKJtmVk+HBp67BHY8YwjJBTYxgh6oKO1qOx1D0AIDmrAFV1jfYTRtgyYroxYwC8K3UVRJJhGCEygnY9Gkvfwy6xZcR0UVGt+3V1ALhUBjkXdmAlou67fBk4dUrss2XEeH37AtqVz48ckbYWIgkwjBBR9xUXAxoN4OcHBARIXY19UqnE9pdfpK2DSAJ8TENkJXY9UuZ6tK0iw4eLhd/IeCoVUA3g0CGpKyGyOoYRIiuw+5Ey11N0Umz5iMZ0I0YA++oZRsgp8TENkRW0HSmTHj8a9VeadK0kDqGoTcsImSYiQmzLy4HSUmlrIbIyhhEiKwoP9LLf0TKdKSoSW7aMmK5nz9b9//5XujqIJMAwQkTdd4otI2bFMEJOhmGEiLqvvl4skDdokNSVOAaGEXIyDCNEZB6DB3OBPHPZvx9obpa6CiKr4WgackqWHGbbdiE9hxvC2xk+ojEPd3fgbI2Y0VbbqZXIwTGMkNOx1DBbP085FD1ckZxVoDumfW+noJ20i7rnhhuA4sPADz8wjJDTYBghp2OpBemUvgrseHaKrsWlsLxW995OYcQIqStwDKNGAd9sA/buBR57TOpqiKyCYYScliWG2Cp9Fc7zWAYQU8BrsWXEPMaMEdvcXGnrILIidmAlItOdOye2Li6cY8RcRo4U299+A86elbYWIithGCEi0504IbbBwaLjJXWft3drK9PevdLWQmQlfExD1EbbUTaAZUfDtB11Y+hju1BYCCAUCA+XuhKHUVheC0yZCZyrB/b+D5gwzblGZZFTYhghanHtKBvAMqNhDI26aXs/P0+5/XR6PXEC8AvlZGdmoPdz4RUDzIsRn1i7x/EWViS6BsMIUYu2o2zCA70sNhrm2lE3bWn/ArabMHL8OBA1FQgfLHUldk/v5+LMGWDmTMDFBYUff4PkT4+YZcQXka1iGCG6RnigF1RKH+MvbGoSy7/n5wMVFUCTAkCo6IQYFKR3qkOMutFoWsIIgEEMI+ag+7lQ+gDyRqCkBDhzQuqyiCyOYYSou06fBv7+d+C998RftFpBg4F5a4AZM8REVuOfEIHFURQXizVpAGDgQGlrcUSTJwMffCDCLW6Quhoii+JoGiJT1dUBzz0nOm+uXCmCSK9ewNSpwMMPA7ffLs6TyYDDh8X+3LnATz9JV7M5/fJL674b/64xu5iWPiMHDkhbB5EV8DcIkSn27gUeeUQ8pgDEPxzPPgvExbUuFqeuBtbuAXJygE8+AeoAHD0KREcDa9cCf/yjZOWb4tqRRsj/DYX+wdIV5OimThXbgweBG6UthcjSGEaIjKHRiCCxaJF45BIcDGRmikcxMpnha/z8xLTea/cAU24GPnwHeOIJ8Xhn+fKOr7MhhkYaAeHA7MVQoBl+nlyt1+yGDgWUSuDKFakrIbI4hhEiY6SlARl/E/sPPQRkZAA+RnR2fX0VoAoD/vxn4P/+T/xDk5ZmmVrN6NqRRgCA3/8eKCqC31tr7b8zri2SyYBp04BvOfEZOT6T+oxkZGQgLCwMHh4eiIyMRG4nayiUlZXhoYcewrBhw+Di4oLk5GRTayWSztWrYrt9u/hHYuVKYMsW44IIIKZNf+klYM0a8fGrr4pAYye0I41Uvm5Q/fffUJ09DuVNY6Quy3FNmyZ1BURWYXQYycrKQnJyMlJTU5Gfn4+YmBjExcWhuLjY4PkNDQ3o06cPUlNTMWrUqG4XTGR1jY3A88+LfTc3EUgWL+7e45WFC0XLCAAsWAB8803367SmQ4fEI6ugIKBvX6mrcVxtw8jFi9LVQWRhRoeR1atXIyEhAYmJiYiIiEB6ejqCg4ORmZlp8PzQ0FCsWbMGc+fOhY+xf0US2YIXXwR27hT7q1cDd99tnvdNTRV9SZqbgTlzALXaPO9rDQUFYjt6tJRVOD6lEggJEft5HFVDjsuoMNLY2Ii8vDzExsbqHY+NjcW+ffvMVlRDQwNqamr0XkRWp9GI7fffAz16iP1Jk8z3/jIZ8I9/iCXjKyqABx9sfRzUorC8FofU1Tikrob6Qr357t1d+fliO4aPaCxu/Hix/e9+AKIzsfZnwuZ+LohMZFQH1oqKCjQ1NSHomtkkg4KCcKbtZE/dlJaWhuXLl5vt/YhMsm4dgBEiNPz1r8Av173CeB4ewIcfAmPHArm5ouVlyRKD69fY1PokbBmxnvHjxc/eDz90un6STfxcEJnIpA6ssmuelWs0mnbHumPp0qWorq7WvUpKSsz23kRdsn59SxgB8MILrXM+WEJ4uJjBFQBefhn49VfdOiVfLpiMLxdMRnr8aNRfabKNNWuamsTcFwDDiDVoW0aKi1F17IRuVJPN/VwQdYNRLSMBAQFwdXVt1wpSXl7errWkO9zd3eHu7m629yMySkEBkJQE9G6Z4vzeey1/z0ceEVN/f/cdkJAA7Nplu+vX/PabmAa+Z08RpMiyvLxa93NzAQSbvn4SkY0yqmVELpcjMjISOTk5esdzcnIwceJEsxZGJJnnnhPzf9x2m/XuKZMBb78NeHoCe/YA77xjvXsbS9tfZNQowNVV2lqczZ49UldAZBFGP6ZJSUnB+vXrsXHjRhw5cgSLFi1CcXEx5s+fD0A8Ypk7d67eNQUFBSgoKEBtbS3OnTuHgoICHNau1UFkKxoaxPb8eWDkSPHIxJpCQ8WMrIAYwVNdbd37d9V+0ZESY8dKW4czysuTugIiizB6Btb4+HhUVlZixYoVKCsrg0qlQnZ2NkJahp+VlZW1m3NkTJse93l5edi6dStCQkJQVFTUveqJzEWjETOh+k8HvL2BT7cDcvGIpLC8tltvbdT1CxaIviq//gq88oqYXM3WaMPIhAnS1uFkCoePBc6dk7oMIoswaTr4pKQkJCUlGfzcpk2b2h3TaIdIEtmqDz8EvvgCmDddzIoaFga/C/V6I1oUPVyNWoPl2hExXbpeLgfeeAOYOVPM0vr442KNElvR1NT617m2YyVZlO7n6KZHAAAKzVWuBUQOh2vTEAHA66+37rf8xa8d0aIdqeDnKTeqQ6nJ18+YIV7Z2WJitI8+6vrXYWnHj4vOq97ewJAhUlfjFHQ/Rzk7gQUL4NfLA8q02VKXRWRWDCPk3LQjw5qbgbi4dp/u7ogWk69/7TXg66+BbdtES0RfGxm18kvLZCvjxol1dsgqlL4KKH83FUg8B5ytAX74wbwT8BFJjL9NyHldvgykpIj94cPFAna2QqUC/vAHsZ+aKm0tbR06JLZ8RGN97u7AnXeK/W3bpK2FyMwYRsh5vfIKcOyY2F+1SvyytyXLl4uF+b791nZGUWhbRhhGpKGd82bbNtGaR+QgGEbILrVdn8PktTm+/lr3qKHQxavbo2bMbtAg0YEVAN580+q3b/s91n1vCgvFlmFEGrGxYhK006eB//5Xd1i7hhHXqSF7xT4jZHeuXZ/D6LU5/vMf3a7fc8lQnDd9xIzFvfQSsGmTmH7ditN6GFwDxRXwq6sWK8n272+9YqiVhwcwe7aYrfejj+D30ivtRmxxnRqyRwwjZHeq6hp163MAQHJWAarqGrv2C/j4ceD554HfvwLMng1l8pPYUX3Z5BEzFtevHzB/PrD1c/GxlYbJt/0ehweK6cj91mdCefEcMCPeKjVQB+69VxdGlCtX6kZsFZbXGvf/ApENYRghu6X9R7LLLl4UHQAvXhYfL10KyGS2uwaM1nPPAR99LfYPHAAGTLfarfXWQNn7b7GNibHa/cmAuDjAxwcoKQF27oRy6lTb/vkl6gL2GSHn0NwsFqP75RfA318cs7UOqx3p1w+45x6x///+nzQ1XL0K7Nsn9hlGpKVQAPEtrVMGJpkkskcMI+QcXnkF+OQTMcNp2wnO7MUjLes95eW1rNxqZfn5QF0d4Osrhh2TtObNE9vt20WLH5GdYxghx/fpp62L3mVmikXw7E1Q39b9//s/699fG4AmT+ZkZ7bgppvEMgGXLnHOEXII/K1CDk29Zz8OLXoJh4IGQ73wOeCxx6QuqXtcXYGcHODHH3WHzDLM+Xq0YYSPaGyDTNbaOvLOO5KWQmQO7MBKDkv9SyGmf1KC+vhXAQCKHi7YYe/zMMyaBaxbLVpHvvqq+8Ocu6KpCdi9W+wzjNiOuXPF0O/cXDH02z9E6oqITMaWEXJM1dWo+uNTqO/hjvT8D5A+eyjqrzTrhvDarUcfFY9JsrOBAwf0huCmx49G/ZUm83+NP/0EnD8vFscbN868702mUypbOzZLMCkekTkxjJDjuXJFzMVw/DgAIPyvLyE8NFDioswkOLh1zZoVK3SHwwO9jB/q3FXffSe206aJ6enJdixYILZbtgDV1dLWQtQNDCPkWJqbxbP0HTvEEEgA6Nu300vsTmqqaB354gvg6FHL308bRmJjLX8vMs7kycCoUUB9PfD551JXQ2QyhhFyHBoNkJQEbN0q/oJ/9VWpK7KMYcOABx4Q+5aed6SurnV+EYYR2yOTtbaObH1f2lqIuoFhhByDRgMsWQL11o9xqG84Dq3bisKho9udVlhea3sL4plAnfwCDvUNR+HPxy17owMHxIRngweLhfvI9jz8sOg/Un5O6kqITMYHwOQY3noL6rffxfTETNTLPYBjAI4V6C18d+2CYja1IJ4R1BfqMf1zNeofSQcAKJqvwM9TbpnOubt2ie0dd5j/vck83N3Fekt/WSM+vnJF2nqITMAwQvZNu3Dc+vWo6umNermH/uJubRa+0y4odu1xe6MbQTMpAOF/fBh+ddVQzvkXqnoPNP/NtEN677rL/O9N5pOYCGS+J/a//BJ4ep6k5RAZi49pyH41NwNpaa0fP/ccgNbF3VRKH73AofRVGDxur8LHDocqZoxYSddSs7JWVQF+fsCUKZZ5fzIPhaJ1ErS33gJq7f9RJDkXhhGyXy++KNbmAIA//7m1U6czeeklsd22TTeU2exmzQJ69LDMe5P53Hef2FZUACtXSlsLkZEYRsj+VFaKbU5O67wXd94pXT1SuvFGMfGVRmPekTXNza37fERjH+Rt+kCtXAmUlEhXC5GRGEbIosy+bkpBATDvEbHfqxfwj38AcJxRMiZ5+WUxxDMnx6TLDf43+uknsfXyAuLizFQoWcXYsWLekSef1PWpssr6RUTdwA6sZDFmXTdF+5f/woWA7wBx7N134XfDECgO7HKIUTImGzkSmDMH+Hav+FjbqbcLOvpvhK++AnpPA267rXXyOLIPS5cC//pM/Dd87z2of3ef5dcvIuomtoyQxZht3ZQLF8RcCk88ATQ0AJMmieMhIVD6KrDj2Sn4csFkfLlgsvP+kn3lldZmeu0ImC4w+N+oorq1lWXWLAsUSxY1aBCwbJnYf+YZVBWesvz6RUTdxDBCFtetdVO++AK44QYxq6qrK/Daa0B6ut4pjjZKxiRt16xZs8bouSb0/ht9/bVo5gfEVONkf557TixqeOECsGQJAAuvX0TUTQwjZJuKi8XomN/9DigrA4YOFX/xL1ki1mWh9h59VGxPnQLeeMP099m6tXVfJuteTSQNNzfgww+B3r2Bw4fFMSMe3xFZG3+rk225eFEsBDdsGJCVJYLHkiWi4+rEiVJXZ9s8PVv3ly0DTpww7X1OngR69jRLSSSh0FDg/fdbA+Xbb0taDlFn2IGVzE59oR5VdY2djm7RnqNTXQ188AHwz3/C76waysuXgVtuAVavBsaM6fL7OrqufO2FMbcDud8Ci14So41a/jG67qyzLX85F/oHc/p3O6b3M3JjNApT/gRUQnQAd3cHMMTqNbX9/72rsx+bco2xtbR17T26cn9L1eiMGEbIrAyNzrh23ZRrz2mlAn7/ChRXG7BjvCuU9/1O9w9pR+/rLPw85dddW0d3zuAZwOAZ4uCbe3Wfv94oCr8930PR2Izk2Ys7vAfZrmt/RloFQIEm+F2qQdWbbwLz1gBN1/6/ZzmmjKoz60i8Tt63rbb36Mr9LVWjs2IYIbNqOzojPNBL99dC2zBSdeyEOOfgRwj/aU/rxcOGofChRCSfdEfVpMlQtumv0NH7OgvtqKHO/grTO+fDD0VnXzc3YMMGFAaGIjmrAFV1jYa/b3V1UD6fjB1V9ahasAhISnK677G9u/ZnpC0/TzmUQUtR9ZqYlwfz5wPrVgMDBli8rrb/7wLo/OewG9cYW0vbzryF5bV69+jK/S1Vo7Myqc9IRkYGwsLC4OHhgcjISOTm5nZ6/q5duxAZGQkPDw8MGjQIb731lknFkv3Qrg+j9FUANTXA3pa/0OPjdbOlhv+0B6qqEqhuHQfV1nVQ7fwS4bOmdv19nUxXRg3pzklOhCpmDFTqY1A98nuEX7nQ+ZuvWgUUF0MZ4AXViwud9nts79r+jLR9KX0VwKJFwGuvihN/+knM3mvCyCtTmTKax1IjgNquX6VS+nR4j67cn6OUzMPoMJKVlYXk5GSkpqYiPz8fMTExiIuLQ3FxscHzT548iRkzZiAmJgb5+fl48cUXsXDhQmzXrilCjkOjAc6fF/sffAAkJIjZIP38xGRlAFBY2Doa5vnngdJSce7UqRy5YU4yGbB+vRiae/YsMP/Jzs//7LPWa9p2hCXHMv02sb3hBjHsNzlZdBZPTxeLIlqTRgNcvizqKCsTHaePHROvU6fEOcXFrdPanzkDlJeLP24aGjg6yMEY/Zhm9erVSEhIQGJiIgAgPT0d3377LTIzM5HWdgXVFm+99RYGDhyI9Ja5ISIiInDgwAG8/vrr+P3vf9+96slytL8o6utbXxcuiHVh2r7Ky4GiIvHLo6gI6NVPPJN+/XXgbJuF25RKsV2xApg8Gdj8M3D//YC/jwRfnJPw8QG++UZ8v0tLxbFPPgGe+INY+K62VnQQ7jFefO7VV0UoJMe3aRPw7XYxcu3kSdFq8txzwK23ill3x40DVCrA3//6fyQ0NorfBRUV+tvz58W2plnM5puQIDqqT0kGoqOBkiMdv2fQYPF75O67xcfz1gAzZ+r/TgEADw/RKdfDo/WlUIiXof2egYBioljZ2FPT+jl4AfAFcnOBAA+goWVhyJMnWycTrK0F6uXi0aerq9HfcuqcUWGksbEReXl5eOGFF/SOx8bGYt++fQav+eGHHxAbG6t37Pbbb8eGDRtw5coV9DCwGmhDQwMaGhp0H1dXVwMAampqjCn3+tLSgO+/b03Y126v3Td0jqFzr7c197ldqa8r79/c3BpALl9u/95dUOvZjOaGS6idGI2aIb8HRowAxo1DrcwLzW//iIMRY4HKS+KcizWoqdH/RVd7sQbNDZdw8EQZai+2/vc+ca6uw2uoEz17At9+i9rHk8X39c0NqH3jDaBvEHD6NE4o/NEcp0LtnIdR88QT4q9Ocli6/79OlaM25g7gu1uB774DPv5YDAUvKBQvZIoL3OVAQADg7gHIe4h/hBsaRABpaBD/QNdd6vSeJ3or0RwXjdqW9Y6ab7qEgz39UdsnrPUkmUzcq2XhyxO+fcX/7woPQKMRNQcEo7btAo5tNQO4BOBSI4BGANWd1DIatRvfRs25k63flz5haP7Dazi47O+oPa9uOW8hau+5R7z9H17DwTsfQO15det7+SvRfMdC1I4fB7i5ofme5ai9+WbUVKtbA4urq9h3cWltEZbJWl9tP+7scx19rD127fuY6sUXRSA1I+2/25rrtWRpjKBWqzUANHv37tU7/pe//EUzdOhQg9cMGTJE85e//EXv2N69ezUANKWlpQavefnllzUA+OKLL7744osvB3iVlJR0mi9MGk0juyZ9aTSadseud76h41pLly5FSkqK7uPm5macP38e/v7+nd7HGdTU1CA4OBglJSXw9vaWuhyHxu+1dfD7bB38PlsHv8/6NBoNLl68iP79+3d6nlFhJCAgAK6urjhz5oze8fLycgQFBRm8pm/fvgbPd3Nzg7+/v8Fr3N3d4e7urnfM19fXmFIdnre3N3/QrYTfa+vg99k6+H22Dn6fW/n4+Fz3HKNG08jlckRGRiJHu6Jni5ycHEzsYKru6Ojodud/9913iIqKMthfhIiIiJyL0UN7U1JSsH79emzcuBFHjhzBokWLUFxcjPnz5wMQj1jmzp2rO3/+/Pk4deoUUlJScOTIEWzcuBEbNmzA4sWLzfdVEBERkd0yus9IfHw8KisrsWLFCpSVlUGlUiE7OxshISEAgLKyMr05R8LCwpCdnY1FixbhH//4B/r374+///3vHNZrInd3d7z88svtHmOR+fF7bR38PlsHv8/Wwe+zaWQazfXG2xARERFZjknTwRMRERGZC8MIERERSYphhIiIiCTFMEJERESSYhghIiIiSTGMOIiGhgaMHj0aMpkMBQUFUpfjUIqKipCQkICwsDAoFAoMHjwYL7/8MhobG6Uuze5lZGQgLCwMHh4eiIyMRG5urtQlOZy0tDSMGzcOvXr1QmBgIO666y4cO3ZM6rIcWlpaGmQyGZKTk6UuxW4wjDiIJUuWXHfufzLN0aNH0dzcjLfffhu//PIL3njjDbz11lt48cUXpS7NrmVlZSE5ORmpqanIz89HTEwM4uLi9OYpou7btWsXnnrqKfz444/IycnB1atXERsbi7q6OqlLc0j79+/HunXrMHLkSKlLsSucZ8QBfP3110hJScH27dsxYsQI5OfnY/To0VKX5dBWrlyJzMxMnDhxQupS7NaECRMwduxYZGZm6o5FRETgrrvuQlpamoSVObZz584hMDAQu3btws033yx1OQ6ltrYWY8eORUZGBl555RWMHj0a6enpUpdlF9gyYufOnj2Lxx9/HO+99x569uwpdTlOo7q6Gr1795a6DLvV2NiIvLw8xMbG6h2PjY3Fvn37JKrKOVRXVwMAf34t4KmnnsLMmTMxffp0qUuxO0ZPB0+2Q6PRYN68eZg/fz6ioqJQVFQkdUlO4fjx41i7di1WrVoldSl2q6KiAk1NTe1W+w4KCmq3yjeZj0ajQUpKCiZPngyVSiV1OQ7ln//8J3766Sfs379f6lLsEltGbNCyZcsgk8k6fR04cABr165FTU0Nli5dKnXJdqmr3+e2SktLcccdd+C+++5DYmKiRJU7DplMpvexRqNpd4zM5+mnn8bBgwfxwQcfSF2KQykpKcEzzzyDLVu2wMPDQ+py7BL7jNigiooKVFRUdHpOaGgoHnjgAXzxxRd6v7ybmprg6uqKP/zhD3j33XctXapd6+r3WfvLpbS0FLfeeismTJiATZs2wcWFWd5UjY2N6NmzJz766CPcfffduuPPPPMMCgoKsGvXLgmrc0wLFizAp59+it27dyMsLEzqchzKp59+irvvvhuurq66Y01NTZDJZHBxcUFDQ4Pe56g9hhE7VlxcjJqaGt3HpaWluP3227Ft2zZMmDABAwYMkLA6x6JWq3HrrbciMjISW7Zs4S8WM5gwYQIiIyORkZGhO3bDDTfgzjvvZAdWM9JoNFiwYAE++eQT7Ny5E0OGDJG6JIdz8eJFnDp1Su/Yo48+iuHDh+P555/nI7EuYJ8ROzZw4EC9j728vAAAgwcPZhAxo9LSUtxyyy0YOHAgXn/9dZw7d073ub59+0pYmX1LSUnBnDlzEBUVhejoaKxbtw7FxcWYP3++1KU5lKeeegpbt27FZ599hl69eun65Pj4+EChUEhcnWPo1atXu8Dh6ekJf39/BpEuYhghuo7vvvsOhYWFKCwsbBfy2LBouvj4eFRWVmLFihUoKyuDSqVCdnY2QkJCpC7NoWiHTt9yyy16x9955x3MmzfP+gURGcDHNERERCQp9sAjIiIiSTGMEBERkaQYRoiIiEhSDCNEREQkKYYRIiIikhTDCBEREUmKYYSIiIgkxTBCREREkmIYISIiIkkxjBAREZGkGEaIiIhIUv8f0oXomMhGe0EAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# affichage de l'histogramme des echantillons\n", "\n", "alpha = 2\n", "plt.figure()\n", "x=np.linspace(-5,5,int(1e4))\n", "y=f(x,alpha)\n", "plt.plot(x,y,\"r\")\n", "plt.hist(Ech, bins=100,density=1, histtype=\"step\")\n", "plt.title(\"HM sample\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Question 3** : Simuler ensuite, en dédiant à la tâche le même temps T que pour la simulation précédente, un échantillon par méthode de rejet (on notera que $\\pi_\\al(x)\\le 2g(x)$ pour $g$ la densité de $\\cN(0,1)$). \n", "\n", "**Question4** : Comparer les histogrammes et afficher le rapport entre les tailles des deux échantillons ainsi simulés. Quelle méthode permet de simuler le plus grand échantillon pendant le temps T ? \n", "\n", "**Question 5** : Calculer ensuite la valeur théorique de l'espérance ainsi que sa valeur empirique calculée grâce aux deux méthodes." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Rapport nombre d'échantillons Rejet / HM pour le même temps de calcul : 1184.865\n" ] } ], "source": [ "# Echantillonnage pendant le même temps que pour HM\n", "# Remarque : pour le rejet, on pourra remarquer que le quotient à calculer se simplifie, ce qui permet de faire un calcul très rapide \n", "Ech=[] \n", "t = time() \n", "while time()-t" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "Estimation la moyenne :\n", "Valeur reelle : 2.707e-01\n", "Methode rejet : 2.717e-01 +/-1.733e-03\n", "Methode HM : 1.748e-01 +/-6.461e-02\n" ] } ], "source": [ "# Affichage de l'histogramme\n", "plt.figure()\n", "x=np.linspace(-5,5,int(1e4))\n", "y=f(x,alpha)\n", "plt.plot(x,y,\"r\")\n", "plt.hist(Ech, bins=100,density=1, histtype=\"step\")\n", "plt.title(\"Reject sample\")\n", "plt.show()\n", "\n", "# estimation de la moyenne\n", "mR=np.mean(Ech)\n", "eR=1.96*np.std(Ech)/np.sqrt(m)\n", "mT = alpha*np.exp(-.5*alpha**2)\n", "print(\"\\nEstimation la moyenne :\")\n", "print(\"Valeur reelle : %.3e\" % mT)\n", "print(\"Methode rejet : %.3e +/-%.3e\" % (mR,eR))\n", "print(\"Methode HM : %.3e +/-%.3e\" % (mHM,eHM))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.5" } }, "nbformat": 4, "nbformat_minor": 2 }