{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "$\\newcommand{\\si}{\\sigma}\n", "\\newcommand{\\al}{\\alpha}\n", "\\newcommand{\\tta}{\\theta}\n", "\\newcommand{\\Tta}{\\Theta}\n", "\\newcommand{\\bet}{\\beta}\n", "\\newcommand{\\Ga}{\\Gamma} \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{\\E}{\\op{\\mathbb{E}}}\n", "\\newcommand{\\Var}{\\operatorname{Var}}\n", "\\newcommand{\\ds}{\\displaystyle}\n", "\\newcommand{\\bE}{\\mathbf{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", "$\n", "\n", "\"Open\n", "\n", "\n", "# Simulation par inversion de la fonction de répartition\n", "\n", "On rappelle que l'algorithme de simulation par inversion de la fonction de répartition consiste à simuler une v.a. $X$ de fonction de répartition $F_X$ dont on sait calculer la pseudo-inverse $F_X^{-1}$ en \n", " 1. Simulant $U\\sim \\mathcal{U}([0,1])$.\n", " 2. Puis en posant $X=F^-(U)$.\n", " \n", "\n", "\n", "## Loi exponentielle\n", "Soit $X\\sim \\mathcal{E} (1)$, sa densité s'écrit $f_X(x) = e^{-x}\\mathbf{1}_{x\\geq 0}$ et sa fonction de répartition s'écrit $$F_X(x)=1 - e^{-x}.$$ \n", "Son inverse est $$F_X^-(u)= -\\log (1-u).$$ \n", "On en déduit que si $U$ suit une loi uniforme sur $[0,1]$, $-\\log U$ suit une loi $\\mathcal{E} (1)$, ce qui nous permet de simuler facilement selon la loi exponentielle." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.stats as sps\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from TP_simulation import *" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHMtJREFUeJzt3Xl8VPW9//HXJwlhERSVKFxCEqwsBgWUVHApitUKtUj7w32pii0Pb61Lq7Zardfqr+21elu7WKtt1VqtltJeQYrSXsVabwUJVrCAQIAEIqIBFUEjS/jcP2YCQwxkOLOcMzPv5+Mxj9lOZt6O8Oabz5w5Y+6OiIjkl6KwA4iISPqp3EVE8pDKXUQkD6ncRUTykMpdRCQPqdxFRPKQyl1EJA+p3EVE8pDKXUQkD5WE9cS9evXyqqqqsJ5eRCQnzZ8/f727l3W0XWjlXlVVRW1tbVhPLyKSk8ysIZntNJYREclDKncRkTykchcRyUMdlruZPWhmb5vZv/Zwv5nZT8yszswWmtkx6Y8pIiL7IpmV+8PA2L3cPw4YED9NBu5LPZaIiKSiw3J39xeAd/ayyQTgEY+ZA/Q0sz7pCigiIvsuHTP3vsCahOuN8dtERCQk6Sh3a+e2dr+7z8wmm1mtmdU2NTWl4alFRKQ96Sj3RqBfwvVyYG17G7r7A+5e4+41ZWUdfsAq2ubPhz59YN26sJOIiHxMOsp9OvDF+F4zo4CN7v5mGh432lasiBX73LlhJxER+ZhkdoV8HHgJGGRmjWZ2uZldYWZXxDeZCawE6oBfAl/JWNoo+le7e4iKiISqw2PLuPv5HdzvwJVpS5QrPP62wqJF4eYQEWmHPqGaKq3cRSSCVO6pev112LYt7BQiIrtRuQfVOpbZtg3q6sLNIiLShso9HTSaEZGIUbkH5Qmf09KbqiISMSr3VJWUaOUuIpGjck/VwIFauYtI5Kjcg2odywwZAsuXw5Yt4eYREUmgck/VkUdCSwssXRp2EhGRnVTuQbWu3I88MnauubuIRIjKPVWDBsXeVNXcXUQiROWeqs6dYcAArdxFJFJU7kEl7ud+5JEqdxGJFJV7Ohx1FKxcCZs2hZ1ERARQuafODIYNi11+7bVws4iIxKncg0ocy7SW+4IF4WQREWlD5Z4OFRXQsye8+mrYSUREAJV7cK0rd7Ndoxmt3EUkIlTu6TJsWGzm3tISdhIREZV7ysxi58OGwYcfwooV4eYREUHlHlziG6rAGTcPB+DsQQuoqgohj4hIApV7mjy7rhqKi/nDzQtoaAg7jYgUOpV7quJjmS10gcGDtceMiESCyj2o+Fimqn+s3ysr0R4zIhIZKvcU1a+K9Xx9PTB8ODQ2ciDvhB1LRAqcyj2oxP3cW8U/qToMrd5FJFwq93RSuYtIRKjcU5W4cj/0UDj0UIajN1VFJFwq96Da7Oe+0zHHcAyvZDeLiEgbKvd0GzGCahbHPq0qIhISlXuqEscyACNGUEKLdokUkVCp3IPa01hmxIjY+fz52csiItKGyj3dyst5i0NU7iISqqTK3czGmtlSM6szsxvbub/CzGab2T/NbKGZfTb9USOmvf3c49fnM0LlLiKh6rDczawYuBcYB1QD55tZdZvNbgGmuPvRwHnAz9MdNJfUHTCC7a8tpqs16wiRIhKKZFbuxwJ17r7S3bcCTwAT2mzjwP7xywcAa9MXMeLartyBqx+Ovana/JKOECki4Uim3PsCaxKuN8ZvS3QbcJGZNQIzgavSki7K9vSGKux6U7W2NjtZRETaSKbcP740ja3UE50PPOzu5cBngd+a2cce28wmm1mtmdU2NTXte9pcUV4OZWWau4tIaJIp90agX8L1cj4+drkcmALg7i8BXYBebR/I3R9w9xp3rykrKwuWOGraGctgFlu9q9xFJCTJlPs8YICZ9TezUmJvmE5vs81q4NMAZnYEsXLP46U5ex/LANTUwOLFdKE5O3lERBJ0WO7uvh34KjALWEJsr5hFZna7mZ0Z3+w64MtmtgB4HLjUvaP2y3MjRkBLi44QKSKhKElmI3efSeyN0sTbbk24vBg4Ib3RIm5P+7m3qqkB4JPMA0ZlJ5OISJw+oZopfftCnz6MZG7YSUSkAKncM8UMRo1SuYtIKFTuQXU0lgEYOZIB1MH69dnJJCISp3LPpFHxWftcrd5FJLtU7kEls3IfMYIWilTuIpJ1KvdM6t6d1zgK5swJO4mIFBiVe4bNYRS8/DLs2BF2FBEpICr3oJIZywBzGQkbN8LSpVkIJSISo3LPsDnoTVURyT6Ve4Z9VDGI9ziAX1w2R1/cISJZo3IPKsmxzKqGInqedixXDJ+rL+4QkaxRuWfDyJGwcCH7sTnsJCJSIFTuQSW5cgfguONgx474QcRERDJP5Z4Nxx8PZpzIi2EnEZECoXLPhp494aij+BR/DzuJiBQIlXtQ+zKWATjxRI7jJdi+PXOZRETiVO7ZcuKJ9GAzLFwYdhIRKQAq92w58cTY+Yuau4tI5qncg9rXsUy/fjRQoXIXkaxQuWfRK/t9ijf/8HfMXJ9WFZGMUrkHta8rd+ALd59IH9bhdSv1aVURySiVezZp7i4iWaJyz6bq6tg+7yp3EckwlXtQAcYyFBXBCSfACy9kJpOISJzKPdtOOgmWLaMPa8NOIiJ5TOWebaecAsDJPB9uDhHJayr3oIKMZQCGD4eePRnD7PRnEhGJU7lnW3ExjB6tcheRjFK5BxV05Q5wyikczgpYvTq9mURE4lTuYRgzJnY+W6t3EckMlXsYjjySd4oO5uFLZ2OGDkUgImmncg8qlbFMUREH/b8xXNrvOXyH61AEIpJ2KvewjBkDa9bAypVhJxGRPJRUuZvZWDNbamZ1ZnbjHrY5x8wWm9kiM/tdemPmIc3dRSSDOix3MysG7gXGAdXA+WZW3WabAcBNwAnuPgS4NgNZoyWVsQzA4MHQuzc8+2z6MomIxCWzcj8WqHP3le6+FXgCmNBmmy8D97r7uwDu/nZ6Y+YhMzjtNPjrXzF2hJ1GRPJMMuXeF1iTcL0xfluigcBAM/tfM5tjZmPbeyAzm2xmtWZW29TUFCxxVKS6cgc4/XTYsIFjeCU9mURE4pIp9/bay9tcLwEGACcD5wO/MrOeH/sh9wfcvcbda8rKyvY1a/457TQATmdWyEFEJN8kU+6NQL+E6+XwsUMaNgLT3H2bu68ClhIre9mbQw6BY45RuYtI2iVT7vOAAWbW38xKgfOA6W22eRIYA2BmvYiNafJ7H790jGUAPvMZjuMleP/91DOJiMR1WO7uvh34KjALWAJMcfdFZna7mZ0Z32wWsMHMFgOzgRvcfUOmQueV00+nE9u1S6SIpFVJMhu5+0xgZpvbbk247MDX4yfZF8cfzya602PWLJjQdickEZFg9AnVoNI1liktZTZjYJbm7iKSPir3CJjF6bHDENTVhR1FRPKEyj2odK3ciZc7wNNPp/xYIiKgco+EFRwOAwfCjBlhRxGRPKFyj4DKSrh72Xi2/OV5hlRsCjuOiOQBlXtQaRzL1NfD9c+PpzNbGbTmryk/noiIyj0qTjgBevZkPE+FnURE8oDKPShve3idFJWUwLhxnMGfoaUlvY8tIgVH5Z6qNIxldho/nkNogpdfTt9jikhBUrlHydixbKdYe82ISMpU7kGl8Q3VnQ48kL/zKXhKc3cRSY3KPWJePHA8vPYa/W0VVVVhpxGRXKVyj5hvz/8CAKvu/hMNDSGHEZGcpXIPKhNjGYD+/eHoo+GPf0zv44pIQVG5R9HEifDSS/SlMewkIpKjVO5BpXs/90QTJwLwBf47c88hInlN5Z6qdI9lAAYPhiFDmIhGMyISjMo9qiZOZDQvwFtvhZ1ERHKQyj2oTL2h2mriRIpwePLJzDy+iOQ1lXtUHXUUyzkcpk4NO4mI5CCVe1SZMYVz4LnnNJoRkX2mcg8q02MZ4HdcADt2wJQpGXsOEclPKvcI+6ByCAsYyktX/06HIhCRfaJyDyqT+7nH1dfDsDsv5DjmUNSwMuPPJyL5Q+WeqgyOZQA47zwALuB3mX0eEckrKveoq6iA0aO5kMey8tuCiOQHlXtQ2SzaCy7gCF6HBQuy95wiktNU7rngrLPYRgk8+mjYSUQkR6jcc8HBB/MU42Plvm1b2GlEJAeo3INyZwcZfjM1wYyySfDWW5xZ+rR2ixSRDqncc8SDa8dC795Mn/CgvqFJRDqkcg8q23uulJTAF78IM2ZwKOuy+9wiknNU7inwLI5lALjsMmhp4SL0xqqI7F1S5W5mY81sqZnVmdmNe9nuLDNzM6tJX0TZafBgOP54JvGg9nkXkb3qsNzNrBi4FxgHVAPnm1l1O9v1AK4G5qY7ZCSFVa6TJlHNEpgzJ5znF5GckMzK/Vigzt1XuvtW4AlgQjvb3QH8APgojfkiLetjGYBzzuF9esAvfpH95xaRnJFMufcF1iRcb4zftpOZHQ30c/cZe3sgM5tsZrVmVtvU1LTPYQXo0YNH+CL8/vewfn3YaUQkopIp9/aWpztnEmZWBPwIuK6jB3L3B9y9xt1rysrKkk8ZRe7hrNyBaX3+HbZs4Yayh7TPu4i0K5lybwT6JVwvB9YmXO8BHAk8b2b1wChgut5UzZy/rh0CJ53EXf3vY01DS9hxRCSCkin3ecAAM+tvZqXAecD01jvdfaO793L3KnevAuYAZ7p7bUYSR0XYe6t85SuwahWnMyvcHCISSR2Wu7tvB74KzAKWAFPcfZGZ3W5mZ2Y6YJSFNZYB4POfh969+Qo/Dy+DiERWSTIbuftMYGab227dw7Ynpx5LOlRaCpMn89nb74C6Ojj88LATiUiE6BOqQYU9lgG44gq20Ql+9KOwk4hIxKjcUxDqWAagTx8e5SJ46CHYsCHcLCISKSr3oKKwcgee6PN1aG7mll73abdIEdlJ5Z6C0FfuxHeLHDeO/3/IT1nXUDAfDhaRDqjc88F118Hbb8e+RFtEBJV7cBEZywBwyikwfDjXczfs2BF2GhGJAJV7CqIwlgHADK6/niN4HaZNCzuNiESAyj1fnHsuyzkcbr89Wr9ViEgoVO5BRa1AS0r4LjfDq6/CU0+FnUZEQqZyT0FkxjJxj3IRHHaYVu8ionIPLILl2UIJ3HwzzJ8PM2d2/AMikrdU7imI2sodgIsvhqoqrd5FCpzKPY9UVoKVduLL9d+Cl1+GGXv9YiwRyWMq96AiuCqur4/F+uXWS1nGALjxRmjRl3mIFCKVewoiOZYB6NSJb/E9WLwYfvObsNOISAhU7nmqtmIicxhJ4+W3Mrjiw7DjiEiWqdyDiuBYJlF9gzHq+Tsp5w0mrPlp2HFEJMtU7imI7Fim1UknwRlncBPf1/HeRQqMyj2oiK/cd7rzTrqzGW65JewkIpJFKvcURH7lDjBkCPdyJdx/P7zySthpRCRLVO4F4D/4DvTqBVddlTu/cYhISlTuQeVQSW6kJ9x5J/zjH/Doo2HHEZEsULmnICfGMsQ+uVo06RLmMJKmS2+AjRvDjiQiGaZyLwD19bDDixg172ccvONtuOmmsCOJSIap3IPKobHMTjU13MO1cN998MILYacRkQxSuacgV8Yyib7NHdC/P3zpS9DcHHYcEckQlXtQubhyBz5kP/jlL2H58thhgUUkL6ncU5CLK/fKSrBTP82vuJzt/3kX1NaGHUlEMkDlXmBaDwv8pXfv5k36wIUXwgcfhB1LRNJM5R5Ujo5ldurZk28c+gg7li3nF92vwyz2BU4ikh9U7inIxbFMosfXjaHoGzdwBffjT06joSHsRCKSLir3QnfHHXD00XD55fRhbdhpRCRNkip3MxtrZkvNrM7Mbmzn/q+b2WIzW2hmz5pZZfqjRkyuj2ValZbCY49BczNPcB5s2xZ2IhFJgw7L3cyKgXuBcUA1cL6ZVbfZ7J9AjbsPBaYCP0h30CjK9bHMTkccAQ88wGj+rk+viuSJZFbuxwJ17r7S3bcCTwATEjdw99nu3vpdbnOA8vTGjKB8Wbm3uvBCftPjSviv/+Ism6o3V0VyXDLl3hdYk3C9MX7bnlwOPJ1KqFyRNyv3uEvW/xBGjWJq98vo2rAk7DgikoJkyr29Bmt32WpmFwE1wF17uH+ymdWaWW1TU1PyKSU7SkvhD3+Abt2Ywedg/fqwE4lIQMmUeyPQL+F6OXx8twozOxW4GTjT3be090Du/oC717h7TVlZWZC80ZFvY5lW5eUwbRp9eQO+8AXY0u7/ShGJuGTKfR4wwMz6m1kpcB4wPXEDMzsauJ9Ysb+d/pjRlG9jmZ1GjeL6Xr+BF1/kkS5fpqoyT/8hE8ljHZa7u28HvgrMApYAU9x9kZndbmZnxje7C+gO/MHMXjWz6Xt4OMkRP2s6F26/nS/yWyat/o+w44jIPipJZiN3nwnMbHPbrQmXT01zrujL17FMoltugYYGbv31HXDPQXDttWEnEpEk6ROqKcjbsUwrM7j/fmZ2mwhf+xqX2UPaRVIkR6jcgyqElTtAcTGffecxOO00Hir6EiMa/hh2IhFJgso9BXm/cm/VuTP86U8wcmTsEAVTp4adSEQ6oHKX5HTvDs88w6udR9Jy9rlcaI9pRCMSYSr3oAplLJNo//355PpnKD55NI/ZxYxpeCjsRCKyByr3FBTMWCZR9+7w5z/HZvBMgrvuKsx/6EQiTuUeVCEXWrduMG0av+cc+MY34KqroKUl7FQikkDlLsF06cKNFY9zN9fBvffylx4T4cMPO/45EckKlXsKCnIsk2BVQxHX+93wk59wavN0OPHE2Ddwi0joVO5BFfJYpq2rrmI8T8HKlVBTA88+G3YikYKncpe0WFR5BgM2zmPRhkNpOfUz8IMfwI4dYccSKVgq9xQU+lgmUX09LPcBDNk0l6e7TYRvfpO/FI/lk+Vvhh1NpCCp3IPSWKZ93bvzuc2/h/vv5zNdX2TmG0NhxoywU4kUHJV7CrRy3wMzmDwZ5s+nqVNfGD+eB2wyQyveCzuZSMFQuQellXvHjjiC6k1z4YYbmFz0a55ZUw1PPhl2KpGCoHKXzOrcOfbm6ty5NFEW++q+s8+GNWs6/lkRCUzlngKNZfZBTQ0TK2q5ie/RPHUGzZWD4Dvf0QefRDJE5R6UxjL7rK6hE9/3m+ha/zrPdh0Pt93G6v0Gc3XZ49ptUiTNVO6SfZWVfO6D38Pf/kbF0b34yfoL4JhjYNo0/aMpkiYq9xRoLJOi0aNh3jyuPfi3LF/wAXz+8yzocizMnKmSF0mRyj0olU96FBdzz/qLGLBtCfz61xywtQnOOAOGDYNHHoGtW8NOKJKTVO4p0Mo9jUpKYNIkBrIMHnwwNoO/5BLo3z+2t82GDWEnFMkpKvegtHLPiH+rLMUmXYYteo2xPM3/rD0CvvlNtvTqCxddBC+8oNdeJAkqd4mU+vpYd7sbz/hYTvX/gQULeLzHl3nvsRlw0kmsKD0i9g1Q2ldeZI9U7inQWCZLhg7l0vd/Ss8P1sLDD/Ne8UGxb4CqqGBul9Hw859DU1PYKUUiReUelEYD2detG1xyCSM++gcsWwZ33EH3LRvgyiuhTx849VT48Y9hxYqwk4qETuUuuWnAALjlFs6o+BdDWcD3W25g8bNr4dpr4fDDWV5aHVvdP/88bNkSdlqRrFO5p0BjmfDVNxgLfSg3+fep9sVQVwf33MO64r5sveseGDOGj7r2jK3qv/td+Mc/YNu2sGOLZFxJ2AFylsYy0fSJT8A11/Cpa66B99+H2bP53cWzGfHsbIY9e0tsm/32g5Ejdz/17h1ubpE0U7mnQCv3iNt/f5gwgUnvT4hdb2ri36v/xpD1zzPquTkMe+4uOrE9dl9lZazkhw+HoUNjp/Ly2LHpRXKQyj0ordxzT1kZ9zWdBZwFwKCKZnqteYWRzGVkw1yObXiZ/lOm7Nr+wAN3FX11NQwcGJv19+0LRZpoSrSp3KVgLV3dFTghforbuJGJg/7FIW8tZNi7Cxj6t4UMfeEhuvvmnZt8SFeWM4DGbgM549qBcNhhUFGx69S1a9b/W0TaUrmnQGOZPHTAAfxx3e6F379yB9tWr2UgyzjuoGV895JlDFu2jO6zFrLte0/uGu20KiuLlXxlZey8b1849NDYXL9379jlXr20+peMSqrczWws8GOgGPiVu/9nm/s7A48AI4ANwLnuXp/eqBGjsUzBWNVQBJTHT6fsvP0TENvz5o03OGfUajq/1UAFq6loWk1lUwMV85dQZc/Qzdv5QpLiYjjkkF2l36sXHHTQ3k89e8Z+TiQJHZa7mRUD9wKnAY3APDOb7u6LEza7HHjX3Q83s/OAO4FzMxFYJFI6dYKqKqasq2r37qpK553VmziUt+jNOnqzLna5ZR2934yd+pW+xbC+S+Gdd2Djxr0/X/fu0KNH+6f27uvWLTYm6toVunTZdbm9k36TyCvJrNyPBercfSWAmT0BTAASy30CcFv88lTgZ2Zm7vm9vNVYRjpS32DA/vHTgHa3qaqChlWxy4dVbGfF/PdiRZ9wuu3qd+Ddd9h/8/v02LyJHm9uogetpzUJlzfRjeZgYTt1gq5dadrclQ93dGYrpdCpEwOqY+eU7uV8T/eVlMR+22jvVFS05/uSPRUVxfZoMtv9ctvre7sv3T/beoL2z8125c+gZMq9L5B4hKZGYOSetnH37Wa2ETgYWJ+OkLu55x749rfT/rD7rLmZHfQPO4Xkgfr6xGslQK/YmCbBbRftwwNu3w6bN8OmTdDczLiTm3n3zWa60kwXPqIrscvlBzVz53c+guZmfvjdZrZsbKbrtmYO6d7MBWdtha1bmfHkNhYt2EontlFK6/nmhMsfP0+8XERer++Cu+8+uOKKjD6FdbS4NrOzgdPd/Uvx6xcDx7r7VQnbLIpv0xi/viK+zYY2jzUZmBy/OghYmq7/kJD0IhP/gOUuvR676LXYnV6P3aXyelS6e1lHGyWzcm8E+iVcLwfW7mGbRjMrAQ4A3mn7QO7+APBAEs+ZE8ys1t1rws4RFXo9dtFrsTu9HrvLxuuRzDso84ABZtbfzEqB84DpbbaZDlwSv3wW8Fy+z9tFRKKsw5V7fIb+VWAWsV0hH3T3RWZ2O1Dr7tOBXwO/NbM6Yiv28zIZWkRE9i6p/dzdfSYws81ttyZc/gg4O73RckLejJjSRK/HLnotdqfXY3cZfz06fENVRERyjz61ICKSh1TuAZnZWDNbamZ1ZnZj2HnCYmb9zGy2mS0xs0Vmdk3YmaLAzIrN7J9mNiPsLGEzs55mNtXMXo//OTku7ExhMbOvxf+e/MvMHjezLpl6LpV7AAmHZBgHVAPnm1l1uKlCsx24zt2PAEYBVxbwa5HoGmBJ2CEi4sfAM+4+GBhGgb4uZtYXuBqocfcjie2gkrGdT1Tuwew8JIO7bwVaD8lQcNz9TXd/JX55E7G/uH3DTRUuMysHzgB+FXaWsJnZ/sBoYnvU4e5b3f29cFOFqgToGv88UDc+/pmhtFG5B9PeIRkKutAAzKwKOBqYG26S0N0DfAPYEXaQCDgMaAIeio+pfmVm+4UdKgzu/gZwN7AaeBPY6O5/ydTzqdyDae+IYQW925GZdQf+CFzr7u+HnScsZvY54G13nx92logoAY4B7nP3o4EPgIJ8j8rMDiT2G35/4N+A/cxsX44atE9U7sEkc0iGgmFmnYgV+2Pu/qew84TsBOBMM6snNq47xcweDTdSqBqBRndv/W1uKrGyL0SnAqvcvcndtwF/Ao7P1JOp3INJ5pAMBcHMjNg8dYm7/zDsPGFz95vcvdzdq4j9uXjO3TO2Oos6d18HrDGzQfGbPs3uhwsvJKuBUWbWLf735tNk8M1lfc1eAHs6JEPIscJyAnAx8JqZvRq/7VvxTzWLAFwFPBZfCK0ELgs5Tyjcfa6ZTQVeIbaX2T/J4CdV9QlVEZE8pLGMiEgeUrmLiOQhlbuISB5SuYuI5CGVu4hIHlK5i4jkIZW7iEgeUrmLiOSh/wN1Yg6x+9km0wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "simu_loi_expo(int(1e5))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loi de Cauchy\n", "La *loi de Cauchy de paramètre $a$* est la loi de densité \n", "$$f_a(x) = \\frac 1 {\\pi}\\frac{a}{x^2+a^2}$$ sur $\\mathbb{R}$. \n", "\n", "1. Calculer la fonction de répartition $F_a$ de cette loi, ainsi que son inverse $F_a^{-1}$.\n", "2. Simuler un grand nombre de variables aléatoires de loi de Cauchy de paramètre $a$, représenter l'histogramme associé et le comparer à la densité. Que se passe-t-il lorsque $a$ s'approche de $0$ ?" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt8VPWd//HXhwQSSLgTRS4y4S4IyGVBqtJSxAtasatssfZXb1131+VXW9u6Vne1Vbu/6nZ1t9W2Wi9ru1qldndFpXVR8VIWWCJyrQQSTCCCEiFcwy3J9/fHOQPTIZBJMjNnZs77+XjMY86cc2bO5zwI7znzPd/zPeacQ0REwqFD0AWIiEj6KPRFREJEoS8iEiIKfRGREFHoi4iEiEJfRCREFPoiIiGi0BcRCRGFvohIiOQHXUC8Pn36uEgkEnQZIiJZ5b333vvUOVfS0noZF/qRSISysrKgyxARySpmVp3IemreEREJEYW+iEiIKPRFREJEoS8iEiIKfRGREFHoi4iEiEJfRCREFPoiIiGi0BcRCRGFvkgLIhEwO/Gh0UIkG2XcMAwimaa6Gpw7cb5Z+msRaS8d6YuIhIhCX0QkRBT6IiIhotAXEQkRhb6ISIgo9EVEQkShLyISIgp9EZEQUeiLiISIQl9EJEQU+iKJWLoUPvMZ+M53mh+TQSRLaOwdkRZ0Yw9ceSXs3euFf2kp3HJL0GWJtImO9EVa8DWegB074O23Ydo0+OEPoaEh6LJE2kShL9KC6/k3r2ln8mS49VbYuhXefDPoskTaRKEv0ozoGPpDrJIxrIM5c7wFl14KXbrASy8FWp9IWyn0RZoRHUO/8meLvBmXXeY9d+4MM2bAokXBFSfSDgp9kVP5wx+gb18YOvT4vAsugE2bOI1PgqtLpI0U+iKnsmQJnHfen94m6/zzAfgM/xNQUSJtl1Dom9klZlZuZhVmdkczy6eZ2UozazCzq+OWXWdmm/zHdckqXCTltm2Dqiov9GNNmAAFBVzAu4GUJdIeLYa+meUBjwKXAqOAa8xsVNxqW4Drgefi3tsLuAeYAkwG7jGznu0vWyQNVqzwns8990/nFxTAxImcV1Cmm6RL1knkSH8yUOGc2+ycOwI8D8yOXcE5V+WcWwM0xb33YmCRc26Xc64OWARckoS6RVJv9Wov0ceOPXHZuHFMKVyDa3I45534FckGiYR+f2BrzOsaf14i2vNekWCtXu2dwC0qOnHZ2LGwZ4/SXrJOIqFvzcxLdPCRhN5rZjebWZmZldXW1ib40SIptno1jBvX/LLo/DVr0lePSBIkEvo1wMCY1wOAbQl+fkLvdc497pyb5JybVFJSkuBHi6ROMfugsvLkoT9mjNf0s3p1egsTaadEQn8FMMzMSs2sEzAXWJDg578GXGRmPf0TuBf580Qy2kg2eBNnn938CsXF3tnb9evTVpNIMrQY+s65BmAeXlh/AMx3zq03s3vN7AoAM/szM6sB5gCPmdl6/727gPvwvjhWAPf680Qy2gjK/YkRp1hpBJSXp6cgkSRJaGhl59xCYGHcvLtjplfgNd00996ngKfaUaNI2g1nI3ToAIMHn3ylESPgnXegqQld5yjZQn+pIs0Yzkav+aag4OQrjRgB9fXw0Udpq0ukvRT6Is0YzsZTN+3A8eVq4pEsotAXieecF/rDh596PYW+ZCGFvki8bdso5kDLod+vn9eLR6EvWUShLxJv40bvuaXQN/PWUehLFlHoi8QrT6C7ZtTw4bBpU2rrEUkihb5IvI0bqacz9E9gmKjSUtiyhTx0o3TJDgp9kXgVFVQyxOun35LBg6GxkQHUpL4ukSRQ6IvEq6riQ0oTW9e/eGswm1NYkEjyKPRFYjkHVVVUEUls/VLvy0GhL9lCoS8Sq64O9u1LPPQHDoS8PEr5MKVliSSLQl8kVlWV95Ro6Ofnw5ln6khfsoZCXyRWa0MfYPBgRnb6UPfLlayg0BeJ1cbQH999M86h++VKxlPoi8SqqoKuXamjZ+LvKS2F2lrYvz9lZYkki0JfxBeJwEs/qWbNvgiDBjV3e+eTiI65/6FO5krmU+iL+KqrYfbYKsZ+IRJt5UmM321ToS/ZQKEvcozXR7/VZ2KjR/qb1YNHMp9CX8TXg92wd2/rQ793b+jaVaEvWUGhL+KLUOVPRFr3RjOviadVbUIiwVDoi/jaHPrR9yj0JQso9EV8x0J/0KA2vDninch1LpkliSSdQl/EF6HKu/1hr15teHPE66e/a1eyyxJJKoW+iC9ClRfe1oo++lHRbptq4pEMp9AX8R0L/Ta92X+fQl8ynEJfxDeIaoW+5DyFvgjA7t30YE/bQ79HD+jeXaEvGU+hLwLHw7o94yJHe/CIZLCEQt/MLjGzcjOrMLM7mlleYGYv+MuXm1nEn9/RzJ4xs7Vm9oGZfTe55YskSbJCX0f6kuFaDH0zywMeBS4FRgHXmNmouNVuAuqcc0OBh4EH/PlzgALn3BhgIvBX0S8EkYySjNA/dlWu+upL5krkSH8yUOGc2+ycOwI8D8yOW2c28Iw//SIww8wM76+/yMzygc7AEWBvUioXSaaqKvbRxj76UZEIHDhAb3YmrSyRZEsk9PsDW2Ne1/jzml3HOdcA7AF6430BHAC2A1uAHznndPWKZJ6qKu9uWW3pox/l/0o4dmWvSAZKJPSb+18Q//v1ZOtMBhqBfkAp8C0zG3zCBsxuNrMyMyurra1NoCSRJIuGfnso9CULJBL6NcDAmNcDgG0nW8dvyukO7AK+DPzeOXfUObcDWAJMit+Ac+5x59wk59ykkpKS1u+FSHtVVVFNG8bcieWP2VOKevBI5kok9FcAw8ys1Mw6AXOBBXHrLACu86evBt50zjm8Jp3Pm6cIOBfYkJzSRZJk927Ys6f9R/o9ekCPHjrSl4zWYuj7bfTzgNeAD4D5zrn1ZnavmV3hr/Yk0NvMKoDbgGi3zkeBYmAd3pfH0865NUneB5H2qa4GaH/oA5SWKvQlo+UnspJzbiGwMG7e3THTh/C6Z8a/b39z80Uyit9dMymhH4kQeb+8/Z8jkiK6Ilck2aFPlcbVl4yl0JdQi0Tg4W9UsZ8iis/snZQPLKKe0zrUYta+a71EUkGhL6FWXQ3fvLKK4lGDqKpuRx/9KD/ldyyvwrljpwtEMoZCX6Sq6vhNUNpLQyxLhlPoi1RVJa8dRqEvGU6hL6HWnd1eP/1khX63bt74PQp9yVAKfQm1Y33qk3nGVUMsSwZT6EuoHQv9ZLXpg26mIhlNoS+hltIjffXVlwyk0JdQi1AFxe0cR/+ED43AoUOwY0fyPlMkSRT6EmoRqryQbs84+vGiTUVq15cMpNCXUDsW+kn9UP/zFPqSgRT6EmopCX1/XH2dzJVMpNCX8Nq9mx7sSX7od+0KvXvrSF8ykkJfwit6JJ6KUdHUV18ylEJfwisaygp9CRGFvoRXNJSTeWFWVGmpP8Sm+upLZlHoS3hVVbGXrtCzZ/I/2++rfzqfJP+zRdpBoS/hVVXl3S0rmX30o/wmo1LUg0cyi0Jfwisa+qngh75uki6ZRqEv4eRcakPf76uv0JdMo9CXcKqrg717Uxf6xcXQp49CXzKOQl/CqbLSe2JI6rZRWqrQl4yj0JfQiUTgmskVABw8I4WhH4ko9CXjKPQldKqr4df3eUf6/10xOHUbioZ+U1PqtiHSSgp9CafKSujXD7p0Sd02IhEKOAIff5y6bYi0kkJfwqmiAoaksGkHNMSyZKSEQt/MLjGzcjOrMLM7mlleYGYv+MuXm1kkZtlYM1tqZuvNbK2ZFSavfJE2qqyEoUNTuw1/eIcvn1eFWWqG+BFprRZD38zygEeBS4FRwDVmNiputZuAOufcUOBh4AH/vfnAvwN/7ZwbDXwOOJq06kXaoAsHYPv21B/p+331n/tBFc75Q/GIBCyRI/3JQIVzbrNz7gjwPDA7bp3ZwDP+9IvADDMz4CJgjXNuNYBzbqdzrjE5pYu0zWA2exOpDv0uXeC003QzFckoiYR+f2BrzOsaf16z6zjnGoA9QG9gOODM7DUzW2lmt7e/ZJH2GYLXcyflzTvgteko9CWDJBL6zY1GFT9e7MnWyQfOB671n79oZjNO2IDZzWZWZmZltbW1CZQk0nZD8frop/xIH2DYMNi0KfXbEUlQIqFfAwyMeT0A2Haydfx2/O7ALn/+2865T51z9cBCYEL8BpxzjzvnJjnnJpWUlLR+L0RaYQiV0KtXaoZUjjd8OGzZAgcPpn5bIglIJPRXAMPMrNTMOgFzgQVx6ywArvOnrwbedM454DVgrJl18b8MPgv8MTmli7TNECrTc5QPXuiD10VUJAO0GPp+G/08vAD/AJjvnFtvZvea2RX+ak8Cvc2sArgNuMN/bx3wEN4XxypgpXPu1eTvhkjihpKGPvpR0dDfuDE92xNpQX4iKznnFuI1zcTOuztm+hAw5yTv/Xe8bpsiwTt6lDPZAkOvTc/2hg3znhX6kiF0Ra6ES3U1+TSm70i/a1c44wyFvmQMhb6ES3m59xw9Ak+H4cMV+pIxFPoSLhs2eM8jR6Zvmwp9ySAKfQmX8nJ2UAK9e6dvm8OHw6ef0oO69G1T5CQU+hIuGzZQzoj0btPvwTMMXaQlwVPoS7hs2MAG0ti0A8dCfzhq4pHgKfQlPHbuhNra9If+4MHQoYNCXzKCQl/Cw++5k/bQ79QJSksV+pIRFPoSHn7PnbSHPsDw4YygPP3bFYmj0Jfw2LABOnWiikj6tz1ypBf6ukm6BEyhL+GxYQMMH04Teenf9ujRdOGgxtaXwCU09o5INotEvFsVbqCcNYyN3sUwvc4+23tety59Q0CINENH+pLzqqvBHTrMiLxK5tw1gqqqAIoY5d9Wev36ADYucpxCX8JhwwZobDx+xJ1uXbtSk3cmz961HjPv14dIEBT6Eg5r13rPY8YEVsKAi0Zz7bj1OOf9+hAJgkJfwmHtWujY8fhNTYIwerT3i6OhIbgaJPQU+hIOa9fCWWd5wR+U0aPh8GGorAyuBgk9hb6Ew9q1gTbtAMfPJ+hkrgRIoS85rwd1UFMTfOifdZb3rNCXACn0JeedzTpvIujQLyqC0lKvr75IQBT6kvPGEHzPnWPGjYNVq4KuQkJMoS85byxroHt3GDAg6FJg/HjYtIli9gVdiYSUQl9y3hj8k7hmQZcCEyaAc4xjddCVSEgp9CW3NTZ6R/rjxgVdiWf8eO+J9wMuRMJKoS+5beNGurIfJk0KuhJPv35QUqLQl8Ao9CW3lZV5z5kS+mYwYYJCXwKj0Jfc9t571NMZRgZwt6yTGT/e60Z6+HDQlUgIKfQlt5WV8T7jIT+Dbh0xfjwdaWBCoUbclPRLKPTN7BIzKzezCjO7o5nlBWb2gr98uZlF4pafaWb7zezbySlb5NQiEcizRg4seZ+NXTOkaSdqwgQAVj7xvkbclLRrMfTNLA94FLgUGAVcY2aj4la7Cahzzg0FHgYeiFv+MPC79pcrkpjqamhct4Ei6rnhkYlBl/OnBg/2rhtYsSLoSiSEEjnSnwxUOOc2O+eOAM8Ds+PWmQ0840+/CMww8zpFm9mVwGZAA45IemXaSdyoDh1gyhRYtizoSiSEEgn9/sDWmNc1/rxm13HONQB7gN5mVgT8HfD99pcq0krvveeNdzNiRNCVnOjcc72RP/fpylxJr0RCv7nLGF2C63wfeNg5t/+UGzC72czKzKystrY2gZJEErBsGUycCHl5QVdyoqlToalJTTySdomEfg0wMOb1AGDbydYxs3ygO7ALmAI8aGZVwDeAO81sXvwGnHOPO+cmOecmlZSUtHonROJ14QCsXAnnnRd0Kc2bMsV7VhOPpFki/dhWAMPMrBT4CJgLfDlunQXAdcBS4GrgTeecAy6IrmBm3wP2O+ceSULdIqc0mf/1boR+/vlBl9K8nj29aweWLg26EgmZFo/0/Tb6ecBrwAfAfOfcejO718yu8Fd7Eq8NvwK4DTihW6dIOp3HEm9i6tRgCzmVqVP9I/341lKR1EnoihXn3EJgYdy8u2OmDwFzWviM77WhPpE2OY8l3u0Je/YMupSTO/dcePpphlAJDA26GgkJXZEruaexkc/wP5nbnh/l/wo59qtEJA0U+pJ71q+nO3szP/RHj4bevZnO4qArkRBR6EvuWeIfOWfqSdyoDh3gc5/jwrzFmDmNwyNpodCXnBGJeCMXz79lMdvyBmRHgk6fzoDGLbjKDzUOj6RFBg09KNI+1dXgGhrhtDfgiisy4/aILZk+3XtevNgbk0ckxXSkL7ll1SrYtQtmzgy6ksScdRacfroX+iJpoNCX3LJokfc8Y0awdSTKzDvaf/NNcOqvL6mn0Jfc8vrrMGaMd/ScLaZPh+3bobw86EokBBT6kjMKOQh/+ANceGHQpbTOxRd7zwsXnno9kSRQ6EvOuIB3vfvOZkt7ftSgQd7Vw6++GnQlEgIKfckZl/MKFBbCtGlBl9J6l10G77xDV/YGXYnkOIW+5AbnmM1L3lF+UVHQ1bTeZZdBQwMzWRR0JZLjFPqSG1avZhBbYHb8nTyzxNSp0LOn92tFJIUU+pIbFiygCYPLLw+6krbJz4eLL2YWC737AIikiEJfslp06IX37nmJlQVTs6urZrwrr+R0djAtf4nG4ZGUUehLVquuBle9hYmsZNL3r2j5DZnsssugc2feueUFjcMjKaPQl+z3/PPe89VXB1tHexUXe8H/4otq4pGUUehL9nvuOe9G40OGBF1J+33pS7BjB7z9dtCVSI5S6EtWG8V6WL0arr026FKSY9Ysr8vp/PlBVyI5SqEvWe1anoW8PPiLvwi6lOTo0sUbFnr+fAo4FHQ1koMU+pK9mpr4Ms95Y+1kc6+deDfcAHV1XMl/BV2J5CCFvmSv118nQjV89atBV5JcM2bAoEHcxJNBVyI5SKEv2euxx6ilD1x1VdCVJFeHDnDjjczkdaiqCroayTEKfclO27bBSy/xNDdAQUHQ1STf9dd7Vxg/9VTQlUiOUehL1olE4B/6PwmNjbzS76+CLic1zjyTxZ1n8cl9j1Foh3R1riSNQl+yzrbqI9zX7+cwcybvfJQDffNPYsYrt3E6Ozj01K91da4kjUJfss61POs173zrW0GXklrTp8PYsfDQQ4DunyvJodCX7NLUxO08COecAxddFHQ1qWUGt90G69ZpnH1JmoRC38wuMbNyM6swszuaWV5gZi/4y5ebWcSfP9PM3jOztf7z55NbvoTOggWcxQa4/XYvFHPd3LnQrx9/z/3gdLQv7ddi6JtZHvAocCkwCrjGzEbFrXYTUOecGwo8DDzgz/8U+IJzbgxwHfCrZBUuIdTUBPfdx2ZKYc6coKtJj4ICuPNOpvEuF3Z4Q0MuS7slcqQ/Gahwzm12zh0Bngfib080G3jGn34RmGFm5px73zm3zZ+/Hig0sxzsXydp8dvfwsqVfI/veTcdCYuvfQ0GDOD1qXfjmpxO6kq7JBL6/YGtMa9r/HnNruOcawD2AL3j1rkKeN85dzh+A2Z2s5mVmVlZbW1torVLmBw9CnfdBaNH8yw5MrhaogoK4O//HpYuhVdfDboayXKJhH5zDafxjYunXMfMRuM1+TTbqdo597hzbpJzblJJSUkCJUnY3NX3Sdi0iSvW/yMDB+UFXU763XADDB8O3/oWHTkSdDWSxRIJ/RpgYMzrAcC2k61jZvlAd2CX/3oA8J/AV51zle0tWELo00+5bdddMG0aC5q+EM6RCTp1gocfho0bmccjQVcjWSyR0F8BDDOzUjPrBMwFFsStswDvRC3A1cCbzjlnZj2AV4HvOueWJKtoCZk776Q7e+DRR8PRY+dkZs2CSy/lHr4P27cHXY1kqRZD32+jnwe8BnwAzHfOrTeze80selPSJ4HeZlYB3AZEu3XOA4YC/2Bmq/zHaUnfC8ldS5fCE0/wr9wKZ58ddDXB+5d/oRNH4G/+Rl04pU3MZdgfzqRJk1xZWVnQZUgGOOvMA7y0dTydOMJlA9eyfkvXoEvKCD/o+SPu2v0druE5lg66JpzNXXICM3vPOTeppfV0Ra5krFu23sFwNhF582kFfoy7Pv0mTJnCr3vN43C1mnmkdRT6klEiEa/Z/gv2Mv+XR+DWW70xaOS4vDx4+mk4eJDn+DI0NARdkWQRhb5klOpqcJsqeLn7/4GJE+GHPwy6pMx01lnw058ynbfgnnuCrkayiEJfMkpn6r07YeXlwYsvQmFh0CVlruuv54WiG+Ef/5HZ9pKGZ5CEKPQlcNEmnXxr4KXOc2HdOnj2WQ0yk4Av1T4CkybxUudrKKleEXQ5kgUU+hK46mpwTY6Gv7yFmQdfhkcegUsuCbqs7NC5M7z8Mpx+Oq9wOWzeHHRFkuEU+pIBnDdU8i9+AXfe6fVBl8T17Qu/+x35NMDnPw8ffhh0RZLBFPoSLOe8C69+9CP427+F++8PuqLsNHIkX+27iF3Ve6ke/Dku6K8jfmmeQl8CEYlARzvKLzrczNf5iXeHqJ/8JNzDLLTTq9sn0GvlGwzqtZ/5n0xjnK3W+PtyAoW+BGJPdR1HL5zFX/KEN2Tyj36kwE+G8ePhrbc4oy+sLj4f9+pCjb8vf0KhL+n3xz+ylKnw9tveRUb336/AT6YxY2D5chg2DL7wBW7nAe+uYyIo9CWNIoMcf2WPcXD0RPp02AWvvw7XXx90Wbmpf3945x246ioe4A4W5l1OidWqqUcU+pIm27fz0JareIy/pvPMC+jz0RqYNi3oqnJbcTG88AL89KfMKniT2tPHMKn6RY3OGXIKfUmtpib42c9g5EhmsRAefBB+/3uvm6GknpnXBXb5cujfnxeZA3/+5/DRR0FXJgFR6EvSRa+wnWFvUJY3GW65hTf2TuLSfmvgO9+BDvqzS7tx42D5cv5fjwc4+F+/p37AMP61x92wb1/QlUma6X+fJF1J9QrczIt4gwuZNHAH/OpXzGh6ncUfDQ+6tHDLz+e7dbfTefMf6TJ3NrfuuY+Puw3j6/Zjuli92vtDQqEvSREZ5JhlC1ls01nBZFi5Eh56CDZuhK98Rb1zMklpKfz617B8OX0/O5Ifcyv1fQYxb/f99LQ69e3PcQp9aZ+dO+HHP+aVLWNYyGVM778J/umfvDFgvvlNjZKZySZPhrfegnffhSlT+Paef6CuywDcDTfSt3qZTvjmKIW+tN6RI7BwIa8UfYnDffrBrbfS1KkQnnnGC/tvfxu6dQu6SknU+efDK6/A6tXer7Lf/IZlTIWxY+Gf/5nzBlRjhn4B5AjdI1cSc/AgvPYa//GV3/L5Ay/Tgz3UdehFz3lfgRtv9E4USm7Yt487Is/z57t+wWT84ZqnTIE5cxjy7SupdEOCrU+aleg9chX60rymJu/Ib9Ei3r1nEX926F0KOewF/XWzvW5/M2dCQUHQlUoqVVbCb34D8+fD++8DUMEQXuNi3i+5mCcqp0NX3b84Eyj0pXXq66GsDJYtg6VL+XTBEvo01QKwoePZjJw3E2bNgs9+Fjp2DLhYCURlJfzud951FosXQ309jXRgFeewput53PCL8+C882DAgKArDSWFvpxcXR2sXQtr1niPlSu9o3r/BtubGMqaos9w1U9nwIUXQr9+ARcsGefwYfjDH7wTwUuWUP/Wcrq4egBq6M/7jGcV51BTMp7Hlp3j9RhSD66UUuiHXVMTbN0KFRXeY9Mm+OADL+y3bj22Wl2HXqxsOodlnEtlyVSe+uO50KdPgIVLVjp6FFatgiVLvF+Mq1bBhg3Q2AjAHrqxgZGUM4JyRrCrzwh+9uYIb1A49fBKCoV+rjtyBLZt8wK8pgZqanj6vhp67atiGJsYzGYKOXxs9YMUsolhrGHsnzw6nXkGVdU6ApMUOHjQu9/xqlXeL8kNG6C83Pt7jTJjW4f+bG4cRDWDqCJCNd50Q79BvLHpTOjSJbh9yCIK/Wx09CjU1nqPHTu8R/x0NOg/+eSEt++1bnQbfaZ39DR06PHnoUO9URc1/IFkgv37uWzYRrp97B33jyn6kKsmVnk3S66pOfbr4JiuXeGMM7zxmqKP2Nd9+kDv3tCrl9dVOKTNSImGfn46igkF57yj7/374cAB2LsXdu9O/LFrl9fW3owG8qilhFpK2EY/ahjHVgZSwwCOnDaQXy0eAAMG0E194yUbFBfz6vYJwIQTlzU0eAc21dV844vVdN65ldP3fULffR9zxsbt9GUVffmY7uxt9qMbyGMXvdib35uhk3sd/zLo3Rt69vS+QLp1O/6If11cnPMHRwkd6ZvZJcC/AnnAE865H8YtLwB+CUwEdgJfcs5V+cu+C9wENAJfd869dqptJfVIPxrEhw/DoUPNP0617ODB4yEefY6djpnXsGc/+TS2WNIBurCbHic86ujJDk6jsddp/ODxEjjtNO9RUgI9euT8H6JIq9TXw8cfe4+dO73Hrl3Hnl/55U4K63fRm530wnsu5kBCH73fiik+w/9CKCrympdiH507nzivuWWdO3tdmgsKvPMW8dOdOiX1V0nSjvTNLA94FJgJ1AArzGyBc+6PMavdBNQ554aa2VzgAeBLZjYKmAuMBvoBr5vZcOdcy+nYWqtWwRe/eGJwt9NBCjlAEfsp5gBFMdP9j83bTzH53YqY93dFUFzM7d8vompXV+roSee+PVjwTg8vuLt3p6hTJ4qA/u3fY5Hw6tIFBg/2Hs24/GfNzDx61BtVdO/e48/RR8zr5x/ai9u2j27spQv1/mMvXfg45nU9nTn4J+fN2iT6RRD9Mpg927tXdAol0rwzGahwzm0GMLPngdlAbOjPBr7nT78IPGJm5s9/3jl3GPjQzCr8z1uanPJjdO8O06bx3H8U8un+Ag5R2OzjMAV07VPIz/+tEAoLueraQqo+8Zb16VfA28u9+dFH5/x8OgOt6c/y4NeTvnci0l4dO3pNPb16nXK1r32vFZ/Z2AgHDzLxrHp21nhfBkP61vPyC/VeS8Hhw8dbE2KmH7z/MAd3H6bw8CEKDh+mgMMUcojNz43kvtRmfkKh3x/YGvO6BphysnWccw1mtgfo7c9fFvfe1BzklpbCM89w7S9bN07Ubz9dTbaBAAAE/UlEQVROSTUiEgZ5eVBczHtbi1v1ttu/1fx8M7gvCWWdSiKh31yjU3ysnmydRN6Lmd0M3Oy/3G9m5QnUdTJ9zPi0He/PFH0gJ/YDtC+ZKFf2A3JsX9qRX4MSWSmR0K8BBsa8HgBsO8k6NWaWD3QHdiX4XpxzjwOPJ1JwS8ysLJGTGZkuV/YDtC+ZKFf2A7QvrZVIl5AVwDAzKzWzTngnZhfErbMAuM6fvhp403ndghYAc82swMxKgWHA/yandBERaa0Wj/T9Nvp5wGt4XTafcs6tN7N7gTLn3ALgSeBX/onaXXhfDPjrzcc76dsA/G1Keu6IiEhCEro4yzm3EFgYN+/umOlDwJyTvPcHwA/aUWNrJaWZKAPkyn6A9iUT5cp+gPalVTJuGAYREUkdXeYpIhIiORH6Znafma0xs1Vm9t9m1s+fb2b2YzOr8Jc3M9hHZjGzfzKzDX69/2lmPWKWfdffl3IzuzjIOhNhZnPMbL2ZNZnZpLhl2bYvl/i1VpjZHUHX0xpm9pSZ7TCzdTHzepnZIjPb5D/3DLLGRJnZQDNbbGYf+H9bt/rzs2p/zKzQzP7XzFb7+/F9f36pmS339+MFv/NMcjnnsv4BdIuZ/jrwc396FvA7vOsFzgWWB11rAvtyEZDvTz8APOBPjwJWAwVAKVAJ5AVdbwv7chYwAngLmBQzP6v2Ba8DQyUwGOjk1z4q6LpaUf80vNHN1sXMexC4w5++I/p3lukP4Axggj/dFdjo/z1l1f74mVTsT3cElvsZNR+Y68//OfA3yd52ThzpO+dih9wr4vgFYLOBXzrPMqCHmZ2R9gJbwTn33865Bv/lMrxrGyBmSAvn3IdAdEiLjOWc+8A519yFdtm2L8eGInHOHQGiQ5FkBefcO3i96mLNBp7xp58BrkxrUW3knNvunFvpT+8DPsC7yj+r9sfPpP3+y47+wwGfxxvKBlK0HzkR+gBm9gMz2wpcC0R7FjU3hEQ2jXV2I94vFcj+fYmVbfuSbfUm4nTn3HbwghQ4LeB6Ws3MIsB4vKPkrNsfM8szs1XADmAR3q/J3TEHfSn5O8ua0Dez181sXTOP2QDOubuccwOBZ4F50bc181GBd1dqaV/8de7Cu7bh2eisZj4qK/alubc1My/wfTmFbKs355lZMfBb4Btxv/SzhnOu0Tl3Dt6v+cl4zaEnrJbs7WbNTVSccxcmuOpzwKvAPSQ4DES6tbQvZnYdcDkww/mNe2TpvpxERu7LKWRbvYn4xMzOcM5t95s8dwRdUKLMrCNe4D/rnPsPf3bW7o9zbreZvYXXpt/DzPL9o/2U/J1lzZH+qZjZsJiXVwAb/OkFwFf9XjznAnuiPwEzlXk3rPk74ArnXH3Molwa0iLb9iWRoUiyTezQKdcBLwVYS8LMzPBGAPjAOfdQzKKs2h8zK4n2zDOzzsCFeOcnFuMNZQOp2o+gz2In6Uz4b4F1wBrgZaB/zBnyR/HaytYS04MkUx94JzW3Aqv8x89jlt3l70s5cGnQtSawL1/EO0o+DHwCvJbF+zILr6dIJXBX0PW0svZfA9uBo/6/x014Q5+/AWzyn3sFXWeC+3I+XpPHmpj/I7OybX+AscD7/n6sA+725w/GOwCqAH4DFCR727oiV0QkRHKieUdERBKj0BcRCRGFvohIiCj0RURCRKEvIhIiCn0RkRBR6IuIhIhCX0QkRP4/G/kX08/B4LcAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "simu_loi_cauchy(int(1e5),3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulation par la méthode du rejet\n", "\n", "On rappelle le principe de la simulation par méthode du rejet. On note $f$ la densité sous laquelle on cherche à simuler. On fait l'hypothèse qu'on connaît une densité $g$ sous laquelle on sait facilement simuler (par exemple une loi uniforme), et telle qu'il existe une constante $c$ telle que $$f(x) \\leq c g(x) \\;\\;\\forall x.$$\n", "L'algorithme du rejet poru simuler suivant $f$ consiste alors en les étapes suivantes :\n", "1. Simuler $X$ suivant $g$\n", "2. Simuler $U$ avec une loi uniforme sur $[0,c]$\n", "3. Accepter la valeur de $X$ uniquement si $U\\leq \\frac{f(X)}{g(X)}$ et sinon recommencer à l'étape 1.\n", "\n", "En pratique, pour que cet algorithme soit efficace, on a intérêt à choisir $c$ le plus petit possible tel que $f\\leq c g$ soit vérifiée. L'algorithme est d'autant plus efficace que le taux d'acceptation est grand, il perd vite de son intérêt quand la dimension de l'espace augmente." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Méthode de rejet pour la loi de Wigner\n", "La *loi de Wigner* est la loi de support $[-2, 2]$ et de densité $$\\frac 1 {2\\pi} \\sqrt{4 - x^2}.$$ \n", "1. Simuler un grand nombre de variables aléatoires de loi de Wigner grâce à la méthode du rejet\n", "2. Représenter l'histogramme associé et le comparer à la densité. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Correction \n", "D'abord, on remarque que la densité est maximale quand $x=0$ et vaut $\\frac 1 {\\pi}$. Si on choisit \n", "$$g(x) = \\frac 1 4 \\mathbf{1}_{[-2,2]}(x),$$ la densité de la loi de Wigner est donc majorée par $\\frac{4}{\\pi}\\times g$, donc on peut appliquer l'algorithme de rejet avec $c = \\frac{4}{\\pi}$." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VOXZ//HPRSDgLhJQ9iDgglpBItaiYFUUV0RF0dKiRakKamv1p9ZW++Daamut4oL7jqi1UlRQUfSRPoCBgoKALLJTQdlEkUC4fn/cEwwhITMkM2eW7/v1mtdsZ+ZcOUm+c+Y+97lvc3dERCQ31Im6ABERSR2FvohIDlHoi4jkEIW+iEgOUeiLiOQQhb6ISA5R6IuI5BCFvohIDlHoi4jkkLpRF1BRQUGBFxYWRl2GiEhGmTx58lfu3ri65dIu9AsLCykuLo66DBGRjGJmC+NZTs07IiI5RKEvIpJDFPoiIjlEoS8ikkMU+iIiOUShLyKSQxT6IiI5RKEvIpJDFPoiIjkk7c7IFUm6LVtgxQpYvBiWLAmXFStg9eptLxs2QElJuGzaFK7r1IH8fKhXL1zy82GvvWDvvaFhw63X19/XjGmrWrCEFiylOWvYG7DtSmndGhYsSPkWkBym0JfstWYNfPopzJ4Ns2b9cL1gAWzevO2yZj8Ed9mlSZMQ6mWXevUY/uIWNq7fRD4l1GMT9dnInqyjIfPZmzXszRr25Bv+VLGW3XeH9u3hwAPhgAPCdYcOtO/VAbP8an8UfThIbTF3j7qGbRQVFbnG3pGEffMNTJwIkyfDlCnhet68H55v0CCE7kEHQdu20LIltGgRLi1bQqNGYS++GmZQ7b9MSQksXx6+QSxdGq4XLoTPPw8fPAsW/PAm+flw2GFwxBHhcuSRcPjhUHfb/bG41is5zcwmu3tRtcsp9CUjrVoFH30EH3wAH34Ygn7LlvBcmzYhQDt35qK/dWTcioNZRCt8B4ew4t2TrpXw/f778IE0fXqou+yyalV4fo894Cc/gW7d4NhjoUsXrEF9hb7skEJfssuWLSEY33wzXCZNCulbvz78+MchII85BoqKYJ99tr4s3pCu7eUS5h6+DUyYED7EPvwQZswIz+22GyO/PZ4zHzoVTjklfEKJVBBv6KtNX9JXSQm88w688koI+hUrQuoedRTccgscf3xoDmnQoMarat06vHU8yyWFGRQWhkvfvuGxr7+G//1fePttOg17Cy7/FwCfcTCv04sRnMdUOlL+ALHa/qU62tOX9FJSAmPHwogR8NprsHZtOMB66qnhcvLJUFAQ99tlTVu4ezge8NZb8MYbMG4clJaG4xPnnQd9+kDHjlgdy46fVxKm5h3JLNOmweOPw/PPh7btvfaCs86C88+HE04IBzx3QtaEfkVffRU+FF9+Gd57L3wAdOjAbct/yf2r+7GCfat9C30ryC4KfUl/q1fDiy+GsJ8yJQT7WWfBz38OPXqE9vpKFBaG5u945ESwffUVvPoqPPVUOCZQty6cdhr88pfh21Hdyltxs/YDMUfFG/pxnZFrZj3NbLaZzTWzGyp5/jIz+9TMpprZR2bWodxzN8ZeN9vMTk7sx5BsUFgYAqbs0sE+42G7jO/2aQ6DBjF1SilX8ncOb7IcXnoJTj+9ysCHEPju8V2yPvAhNHf96lfwf/8XDv7+5jch/Hv1Cs0/f/7zDz2DRNx9hxcgD5gH7A/kA9OADhWW2bPc7TOB0bHbHWLL1wfaxN4nb0fr69y5s0t2AXcvLXUfNcq9R4/wQIMG7gMGuE+evHW51q3ji/LWraP6STJISYn7P/7h/tOfho22yy7uAwe6f/rp1kUgwvqk1gHFXk2eu++o4/IPugBz3X2+u5cAw4FeFT441pW7uxtQ9qWxFzDc3Te6+xfA3Nj7Sa7YvJmf8wx06BD24GfMgNtvD0MgPPZY6E8fU3bOkvbea0G9etC7d2jv/+QT+NnP4JlnwolgZ5wRvglIToon9JsDi8vdXxJ7bBtmNsjM5gF/Bq5K5LWShTZuhEcfhQMO4Bn6h+aa558Pif273yXUA0dq6LDDwu9iyRIYMgT+/W84+mje4UR4/3017OeYeEK/st7L2/2VuPtQd28LXA/8PpHXmtlAMys2s+KVK1fGUZKkg4pt9WZQzzZxuT3E4gbtYOBAJn1RwCWNX4f//AcuvDDsgUo0GjWCP/whHBS55x4OYUY41+GYY8KZzZIT4gn9JUDLcvdbAMt2sPxw4KxEXuvuw9y9yN2LGjduHEdJkg62OaBaugV/aQSb2nXgIa6gZdfWMGYMXbZM5LEVZ8Y1ro2kyO67w29/Sxu+gKFDwy/yuONCT59p06KuTpIsnv/Ej4H2ZtbGwnCAfYGR5Rcws/bl7p4GzIndHgn0NbP6ZtYGaA9MqnnZklbGjg1nyZ5/PuyyC4waFc4kPemk+E5zlUhspAFccQXMmRN6+EyYAJ06Qb9+8MUXUZcnSVJt6Lv7ZmAwMAaYCYxw9xlmNsTMzowtNtjMZpjZVOAaoH/stTOAEcBnwGhgkLuXJuHnkFpUWbNNZZduzeaGg7MnnhiGSHj66dCMc9ppCvtMsssucN11YRC4668Pff4POghuugnWr4+6OqllOjlLtlPtSTvffgt33AH33BMO0N58MwweXCtj4EjqVPl7XroUbrwRnn0WmjeHu+8O4wHpgzyt1erJWSJASIiXXw57gXfcEZpzZs+Ga69V4GeT5s1D987x42HffcMB+O7dw1DQkvEU+hKfRYvCgb7zzgvdLT/6KARD06ZRVyY7qWxk0SovXX9C3pRJ3LDPMJg5M7T333xz6I4rGUuhLzu2ZUvo4XHIIeHg7H33QXExdO0adWVSQ/GcDFfqeQzf41IKvprJ05svhFtvZWaDjhxjH233IVFYGPVPJPFQ6EvVZs8OX+sHDw4zOU2fDlddBXl5UVcmKbRgAXzlBfT3p2H0aA5uvYGPOBa/YhD+zfqtHxDxDoIn0VLoSyUcHngAOnYMwyY89RSMHq1dOQnzGUyfDr/+NTz0UGjy0ZAOGUWhL9v67395k1PhyivD2ZozZkD//uq5IT/YfXe4994wkcumTeGM3ltuoS6boq5M4qDQlx+MHAmHHcZxjAvt+KNG6UCtVK1bt3AG74UXwpAhjKcrfP551FVJNRT6EnpjDBoUxl9v2ZIjmBLO1NTevVRnr71CL64RI2jLPOjcOcyJIGlLoZ/rFiwIX88ffBB++1uYMIFZHBx1VZJp+vShI1PhRz8KJ3INGqSunWlKoZ/L3ngjjGf/+edhvtV77tnpuWhFltAytPP/9rdhJ6JrV5g/P+qypAKFfi4qLQ3jqpx+ejhDZ8qUMDetSE3Vqxd2Hv75zzCWzxFHwFtvRV2VlKPQzzXr1oW2+zvugAEDwoQabdtGXZVkgW3O8D2rF23WTOE/a9uw5dTTuNbuwcx1ElcaUOjnkGObz2fGXkez+Y3RXMFQ7PHHsF132e7Mytato65UMlHFM3y/8DZ0Wv8Rdc49h3u4Dv95f3zD9zqJK2IK/VwxbhyvLevCIQ2XU3fs2zzoV2gOWkm+3XaDESPgf/4njNrZvTtNdzgHkySbQj8XPP449OjBShrDpEnhpCuRVDELA7W9+ipMn85EjtKInRFS6Gcz9zAR9iWXwPHH82MmQLt2UVcluerss2H8ePIoDd2Ex42LuqKcpNDPVqWlcPnlcMstYRiFUaNYx15RVyW5rmPHsPPRrFkYx2f48KgryjkK/Wy0YQOccw488kiYAenJJ0NXOpE0sJhWYYKWo46CCy6Av/wl6pJyikI/C5Sf07ahrWb8riey5fWRDOZ+7M47sDqmXjmSXho2hLffhj59wsxrN91UzRydUlvqRl2A1NzChbH/l6++gh494LPP4LmXeKBPHx6IujiRqjRoAC++GMbvueOOMPfyvfdqzKckU+hniy+/hBNOCGdBvv469OwZdUUi1cvLg2HDQtfO++6D9etDs6Qm6kkahX4WaMZS6H4CLF4cxtNRl0zJJGZhD3+PPeC228Ie/zPP6DhUkqhNP9MtXMiHdINly2DMGAW+pL1KJ2SvY9htt3I9d8Hw4byU34+2rTdHXWpWimtP38x6AvcBecBj7n5XheevAS4BNgMrgV+6+8LYc6XAp7FFF7n7mbVUuyxdCscfzz6sgnffhS5doq5IpFo7PuP7ergnj/Ovu45Ni+pB6dNq6qll1e7pm1keMBQ4BegAXGBmHSos9h+gyN1/BLwC/LnccxvcvWPsosCvLWVt+CtXcjJjFPiSPa69Fm6/nX48D5deClu2RF1RVomneacLMNfd57t7CTAc6FV+AXd/392/i92dALSo3TJlG6tWhV46ixbBG2/wMQp8yTK/+x1/5JZwjskVV6g7Zy2KJ/SbA4vL3V8Se6wqA4DyA2g3MLNiM5tgZhq0vabWrg1nMn7+eZjT9thjo65IJCn+h1vCyYWPPALXXKPgryXxtOlX1mm20q1vZv2AIqB7uYdbufsyM9sfeM/MPnX3eRVeNxAYCNCqVau4Cs92hYVsNwRtAzbwNqdzFNPozWu82eNEQCddSbYyuP12+O47+NvfYL/94Prroy4q48Wzp78EaFnufgvYfmxUMzsRuAk40923To7p7sti1/OBcUCniq9192HuXuTuRY0bN07oB8hWZSdcbb1sLmVD759xrI0n/6XneMNP01DIktVatw69eurc91de4AK44QYutie36/mjSVkSE0/ofwy0N7M2ZpYP9AVGll/AzDoBjxACf0W5xxuaWf3Y7QKgK/BZbRWfM9zhqqvCPLb33gvnnRd1RSJJVzYpyxavw4Ubn4IePXgy71J85L+22SHSpCyJqTb03X0zMBgYA8wERrj7DDMbYmZlvXHuBnYHXjazqWZW9qFwMFBsZtOA94G73F2hn6i77goTTV93HVx9ddTViKRefn4Yj79Tp7DT8+9/R11RxjJPs4MjRUVFXlxcHHUZkTOLHbd6+mm46CK48MIw81AdnU8nOWzFijAW/1dfwcSJ0L79D/8rOc7MJrt7UXXLKUHS2fvvhwlQTjghdF1T4Euua9IE3nor/C+cfjqsXh11RRlHKZKm2jIXzj0X2rcPX2vz86MuSSQ9tG0L//gHfPEFnHsuddkUdUUZRaGfjtauZSSxwyX/+lcYelZEftCtGzz6KLz3Hg8wWO07CVDop5vSUujbl/bMCXv4bdtGXZFIeurfH268kV8xLPTjl7go9NPNddfB6NEMYigcd1zU1Yikt9tu41XODuP1vPde1NVkBIV+OnnhhdAP/8oreTScoCwiO1KnDhfxFBx4IJx/fphTQnZIoZ8uZswIIwoee6wmihZJwHr2CCcubtwI55wD338fdUlpTaGfDtatg7PPDjMHvfSSZgwSSdSBB4ZzWj7+OJy9LlVS6EfNHQYMCHPbjhgBTZtGXZFIZurdO4zK+eij8PjjUVeTthT6KVZYuO1gUb+p8zd45RWuLb0L695t6+MaOVMkPuWnX8y781bepgcbLhnMoTZdA7NVQsMwpNg2p4xPmgRdu8IZZ4TumVbZKNYikpAvv4TDD4dGjUJzz667AmT9cA0ahiHdffNNGE+nWTN44gkFvkht2XdfeOYZ+Owz+M1voq4m7Sj0o3LlleE08ueeg733jroakexy0klhwpVhw+Dll6OuJq0o9KPw4ouhp8Hvf6/pDkWS5dZb4aijQldozTS0lUI/xVqzAC67DI4+Gv7wh6jLEcle9eqFHSx36NePOpRGXVFaUOinUmkpz9Ev3H7+eagbzxTFIrLT2rSBoUNh/Hh+w71RV5MWFPqpdO+9HMP48EfYpk3U1Yjkhp/9DHr35nZuCme+5ziFfqrMmgW//z3/pFf4IxSR1DCDhx9mHXvCL34Bm3J7/H2FfiqUlsLFF8Nuu3EZD6t7pkiqNWnCr3gEpkyBO+6IuppIKfRT4a9/hQkT4P77+ZL9oq5GJCe9xtnhW/Ztt4Xwz1EK/WSbOTP00undGy64IOpqRHLb/fdDQUHoxrl5c9TVREKhn0xbtoQ/rt12g4ceUrOOSNQaNgzBP2UK/P3vUVcTCYV+Mj3xBIwfH8bH33ffqKsREQhj7p9xRvgGnoMnbSn0k2XFCvh//w+6dw9zeYpIejAL3abr1IHLL8/uUdgqEVfom1lPM5ttZnPN7IZKnr/GzD4zs0/MbKyZtS73XH8zmxO75E76XXstrF+vZh2RdNSyJdx+O4weHSYuyiHVhr6Z5QFDgVOADsAFZtahwmL/AYrc/UfAK8CfY6/dB7gFOAroAtxiZg1rr/w09f778OyzYU//4IOjrkZEKjNoEBx5JFx9NaxZE3U1KRPPnn4XYK67z3f3EmA40Kv8Au7+vrt/F7s7AWgRu30y8I67r3L31cA7QM/aKT1NbdwYvjK2bQs33RR1NSJSlbw8eOQRWLkS/vjHqKtJmXhCvzlQfor5JbHHqjIAeCuR15rZQDMrNrPilStXxlFSGrv3Xpg9O7QZ7rJL1NWIyI506gQDB8IDD+TMEA3xhH5lDdKVHvkws35AEXB3Iq9192HuXuTuRY0bN46jpPRTWAhNbTnf3Hg7/6QX1vPkbaZq0zSIItEqP61i+UvBI7exqnRP3j30asw866dWjCf0lwAty91vASyruJCZnQjcBJzp7hsTeW02WLgQll98E3vU28hZc+7BnUovOdhDTCQtLFhQ+f/kV17APg/cyomMxV99Dffw/5yt4gn9j4H2ZtbGzPKBvsDI8guYWSfgEULgryj31BjgJDNrGDuAe1LssazTmWJ48kn49a+hXbuoyxGRRPzqV3DYYXDNNfDdd9Uvn8GqDX133wwMJoT1TGCEu88wsyFmdmZssbuB3YGXzWyqmY2MvXYVcCvhg+NjYEjssezizt/4NTRpEmbDEpHMUrduOFN34cIwVlYWM0+zExOKioq8uLg46jISM3x4GFfnscdgwICoqxGRndW7N4wdS5Nv5rLCm0RdTULMbLK7F1W3nM7IramSErjxRv5DR7jooqirEZGauPNO+O47/sCtUVeSNAr9mnrkEViwgBu4K/T7FZHMddBBcMklYd6LuXOjriYpFPo18c03cOut8NOf8jYnRV2NiNSGP/6RjdSH3/0u6kqSQqFfE/feG87mu/NOKj8lQUQyzn778Rd+Cy+/DBMnRl1NrVPo76yVK+Gee+Dss+Goo6KuRkRq0T1cG3rj3bDd+JIZT6G/s+64A779Nky9JiJZZT17hOadcePCJYso9HfG0qXw4IOht45G0RTJTgMHQtOmcMstWTXmvkK/GoWF24/VcV+LP7OpZAuFT/xBY+qIZKtddoEbb4QPPwzDpWcJhX41Fi6sMFbHsuVc3WAY9Qb0Z4EXakwdkWx26aXQvHlW7e0r9BN1992waVPYAxCR7NagQfhf/+gjGDs26mpqhUI/EV9+CQ8/DP36hUlSRCT7XXIJtGiRNROtKPQT8Ze/hJmxsvSkDRGpRP36YerT8ePDHn+GU+jHa9Wq0GOnb1844ICoqxGRVPrlL6FRI/jTn6KupMYU+vF66KHQL19t+SJZb7tZtnbfjZu/vgpGjeJQm7718UycYUuhH4/vvw9jbffsCYceGnU1IpJklc2yNeSrQbDrrkz/xd1bH8vEGbYU+vF4/vlwEPfaa6OuRESi0qhR6ML5wguwaFHU1ew0hX41jC3hAG7HjnD88VGXIyJRuuaacJ3Bs2sp9KtxCm/BzJlw3XWhEU9EclerVmGWvMcfh7Vro65mpyj0q3Edd0PLltCnT9SliEg6uPpqWL8ennoq6kp2ikJ/R6ZO5Tg+CL/kevWirkZE0kHnzvCTn8D994fm3wyj0N+RoUP5jl1CH10RkTJXXQXz5oXm3wyj0K/KmjXw/PO8wIXQsGHU1YhIOjn7bGjWjCu5P+pKEqbQr8pTT8GGDQxlUNSViEi6qVcPLr+cnoyBWbOiriYhcYW+mfU0s9lmNtfMtps/zMy6mdkUM9tsZudWeK7UzKbGLiNrq/Ck2rIlDLlw9NFMpVPU1YhIOho4kI3kw9ChUVeSkGpD38zygKHAKUAH4AIz61BhsUXARcALlbzFBnfvGLucWcN6U+Pdd2HOHBikvXwRqUKTJrzKOfDcc7BhQ9TVxC2ePf0uwFx3n+/uJcBwoFf5Bdx9gbt/Ahl4KLsyDz4IjRvDuedWv6yI5KxHuTQc/3v11ahLiVs8od8cWFzu/pLYY/FqYGbFZjbBzM5KqLoo/Pe/MGoUXHxxGFJVRKQK4zgO2rWDxx6LupS4xRP6lZ2Gmsi8Ya3cvQi4EPibmW03+4iZDYx9MBSvXLkygbdOgmefhdJSddMUkTgYDBgAH3wAn38edTFxiSf0lwAty91vASyLdwXuvix2PR8YB9sfGXX3Ye5e5O5FjRs3jveta587PPlkOPHiwAOjq0NEMsdFF0FeXhiaIQPEE/ofA+3NrI2Z5QN9gbh64ZhZQzOrH7tdAHQFPtvZYpNu0qQwzs7FF0ddiYhkiv32gzPOCN28S0qirqZa1Ya+u28GBgNjgJnACHefYWZDzOxMADM70syWAH2AR8xsRuzlBwPFZjYNeB+4y93TN/SfeAJ22QXOOy/qSkQkk1x6KaxYAW+8EXUl1TL3RJrnk6+oqMiLi4tTv+LvvoOmTeGss+Dpp7c+bBZafUREKtqaD5s3h8nTu3aNrCePmU2OHT/dIZ2RW+Yf/4B169S0IyKJq1s3DLk8ahSsXh11NTuk0C/z/PNhwstu3aKuREQyUb9+oU0/zfvsK/QBVq6Ed96Bvn2hjjaJiMRnmwnUi45gFgcy7tLntp1UPc0mUM/ZhCss/OEXclmTV6G0lMPv6rvdL6t166grFZF0te0E6sZBt/bjOD7AFy7aZlL1dJpAPWdDf+HCH34hD3cfDgcfzLQtP9rmF+UefqkiInG58MJw/eKL0daxAzkb+lstXQoffhiadjQHrojUxP77h5M7X6hs7Mn0oNAfMSLs0vftG3UlIpIN+vSBTz6BuXOjrqRSCv3hw+GII+CAA6KuRESywdlnh+s07cWT26G/aFEYekFn4IpIbWnVCrp0gVdeibqSSuV26L/+ergu+2QWEakN55wDxcXp1W0nJrdD/5//hA4doH37qCsRkWxyzjnhOg2beHI29BuyKoyB3atX9QuLiCSibVvo2FGhn05O440wWcpZ6T+Zl4hkoN694f/+L5zxn0ZyNvR78To0awZF1Q5KJyKSuNNOC93B33or6kq2kZuh//339GR0aNrRWDsikgydOoUJVtJsjP3cTLwPPmB3vg2z3YiIJEOdOmFvf/Ro6rIp6mq2ys3QHzOG76kP3btHXYmIZLPTToN16+jK+Kgr2So3Q3/0aD6gO+y6a9SViEg2O/FEqFcvdBxJE7kX+gsXwsyZjKZn1JWISLbbYw/o3l2hH6kxY8IVJ0dciIjkhFNPpQMzw7AvaSA3Q79lS2ZycNSViEgu6NEjXI8dG20dMbkV+ps2wbvvQs+egMbOF5EUOOQQ/su+IXvSQG6F/sSJsG4dnKymHRFJETPe4/iwp+8edTXxhb6Z9TSz2WY218xuqOT5bmY2xcw2m9m5FZ7rb2ZzYpf+tVX4Thk3Llwfd1yUVYhIjnmXE+HLL2HGjKhLqT70zSwPGAqcAnQALjCzDhUWWwRcBLxQ4bX7ALcARwFdgFvMrGHNy95J48bBj34EjRpFVoKI5J6xnBC7EX27fjx7+l2Aue4+391LgOHANkNTuvsCd/8E2FLhtScD77j7KndfDbwDEfWVLCmBf/9bJ2SJSMotojW0a5cxod8cWFzu/pLYY/GoyWtrV3ExbNigph0Ricbxx4fh3EtLIy0jntCvrJtLvEcj4nqtmQ00s2IzK16ZrGFIy9rzu3VLzvuLiOzIMceEjiQRt+vHE/pLgJbl7rcAlsX5/nG91t2HuXuRuxc1btw4zrdO0AcfwKGHQkFBct5fRGRHunYN1+OjHYcnntD/GGhvZm3MLB/oC4yM8/3HACeZWcPYAdyTYo+l1qZNYUOrPV9EotKmTRhqOd1D3903A4MJYT0TGOHuM8xsiJmdCWBmR5rZEqAP8IiZzYi9dhVwK+GD42NgSOyx1PrkE/j2Wzj22JSvWkQEALPQxBNx6NeNZyF3fxN4s8JjN5e7/TGh6aay1z4BPFGDGmtuwoRw/eMfR1qGiOS4rl3hlVdg6VJoHk2fltw4I3fixPC1qlWrqCsRkVx2zDHhOsK9/dwI/QkTwl6+abwdEYnQ4YeHeTw++iiyErI/9L/+GubMUdOOiESvXj048kiYNCmyErI/9CdODNcKfRFJB507w9SpoVdhBLI/9CdMCBMUFxVFXYmISMiijRvhs88iWX32h/6UKdChA+y2W9SViIiEPX2AyZMjWX32h/7UqdCpU9RViIgE7dqFuXMV+kmwcmXoD9uxY9SViIgEderAEUeEQSCjWH0ka02VqVPDtUJfRNJJ584wbVokB3OzLvQLC0N3fDO47qQQ+o1OOHzrY2WX1q2jrVNEcljnzpEdzI1rGIZMsnBhuWkofzYV/rclXy/STFkikkYOPzxcT5/+w+0Uybo9/W1MnaqmHRGJXOvW27Y05B/anhLqcWe/6ds8XliY/FqyN/Q3bIBZs9RzR0Qit2BBaIEou5R4PvmHHsiNp0/f5vGFC5NfS/aG/qxZsGULHHZY1JWIiGzv0END806KZXfoAxx8cLR1iIhU5tBDw1eAb75J6WqzO/Tr1AknQoiIpJtDDw3XKe7Bk72hP3Mm7L8/1K8fdSUiItsrC/0UN/Fkb+jPmgUHHRR1FSIilSssDEMtf/55SlebnaFfWho2pNrzRSRd5eWF1oi5c1O62uwM/YULw9lu2tMXkXTWvr1Cv1aU9dw58MBo6xAR2ZF27ULobx1GIPmyM/TLPjnbt4+2DhGRHWnXDr77Dv7735StMjtD/4svwuTDjRtHXYmISNXKupTPmZOyVWZv6LdpEwazEBFJV2WtESls148r9M2sp5nNNrO5ZnZDJc/XN7OXYs9PNLPC2OOFZrbBzKbGLg/XbvlVKAt9EZF01qoV1K2b0tCvdmhlM8sDhgI9gCXAx2alKn+tAAAJUElEQVQ20t3Ln0Y2AFjt7u3MrC/wJ+D82HPz3D2FQ116CP3u3VO3ShGRnVG3buivP29eylYZz55+F2Cuu8939xJgONCrwjK9gKdjt18BTjCLpm1lH1aFsSy0py8imaBVK1iyJGWriyf0mwOLy91fEnus0mXcfTOwFiibuaSNmf3HzD4ws2MrW4GZDTSzYjMrXrlyZUI/QEVt+CJ2Q6EvIhmgefO0C/3K9tgrdiqtapnlQCt37wRcA7xgZntut6D7MHcvcveixjXscdOC2MZr2bJG7yMikhItWsCyZWEkgRSIJ/SXAOUTtAWwrKplzKwusBewyt03uvvXAO4+GZgHHFDTonekKcvDjWbNkrkaEZHa0aIFbN4MK1akZHXxhP7HQHsza2Nm+UBfYGSFZUYC/WO3zwXec3c3s8axA8GY2f5Ae2B+7ZReuWYsC0MqN2mSzNWIiNSOFi3CdYqaeKrtvePum81sMDAGyAOecPcZZjYEKHb3kcDjwLNmNhdYRfhgAOgGDDGzzUApcJm7r0rGD1KmKctD4OflJXM1IiK1Y5vQPzLpq6s29AHc/U3gzQqP3Vzu9vdAn0pe9yrwag1rTEhTlkPTpqlcpYjIzisL/aVLU7K6rDsjtynL1Z4vIpmjoADy82Hx4uqXrQXZGfra0xeRTFF2DDKNDuRmjtJSmrBCoS8imaVxY6jhOUrxyq7QX7GCPLYo9EUkszRuDF99lZJVZVfof/lluN5332jrEBFJREGB9vR3ypo14XqffaKtQ0QkEdrT30mrV4frvfeOtg4RkUQUFMC6deSzMemryq7QL9vTV+iLSCaJjTlWQPL39rMz9Bs2jLYOEZFEFBQA0Iivk76q7Ar91avZgsEee0RdiYhI/PbaC4A9WZf0VWVX6K9Zw1r2Cic7iIhkij3DiPN7sTbpq8qudFy9mtWoaUdEMkws9LWnn6g1a1iDDuKKSIaJNe9oTz9RZc07IiKZRM07O2nDBr5lt6irEBFJTIMGANRXP/0ElZRQQn7UVYiIJCYvD/LyyKck6atS6IuIpIP8fIV+whT6IpKpFPo7QaEvIplKob8TFPoikqkU+jtBoS8imUqhvxMU+iKSqRT6O0GhLyKZKp1C38x6mtlsM5trZjdU8nx9M3sp9vxEMyss99yNscdnm9nJtVd6BaWlUFqq0BeRzJQuoW9mecBQ4BSgA3CBmXWosNgAYLW7twPuBf4Ue20HoC9wCNATeDD2frVv0yYAhb6IZKZ0CX2gCzDX3ee7ewkwHOhVYZlewNOx268AJ5iZxR4f7u4b3f0LYG7s/WpfSdhYCn0RyUhpFPrNgcXl7i+JPVbpMu6+GVgLNIrztbUjFvqbqJeUtxcRSaoUhX7dOJaxSh7zOJeJ57WY2UBgYOzuejObHUddVbiq4AG7KjXTyiemAFIwAWbiVFdiVFdiVFdiCjDb2bpax7NQPKG/BGhZ7n4LYFkVyywxs7rAXsCqOF+Luw8DhsVTcHXMrNjdi2rjvWqT6kqM6kqM6kpMLtcVT/POx0B7M2tjZvmEA7MjKywzEugfu30u8J67e+zxvrHePW2A9sCk2ildREQSVe2evrtvNrPBwBggD3jC3WeY2RCg2N1HAo8Dz5rZXMIeft/Ya2eY2QjgM2AzMMjdS5P0s4iISDXiad7B3d8E3qzw2M3lbn8P9KnitbcDt9egxkTVSjNREqiuxKiuxKiuxORsXRZaYUREJBdk1zAMIiKyQxkf+mZ2t5nNMrNPzOw1M9u7iuV2OJREEurqY2YzzGyLmVV5NN7MFpjZp2Y21cyK06iuVG+vfczsHTObE7tuWMVypbFtNdXMKnYoqM16dnrokWSKo66LzGxluW10SQpqesLMVpjZ9CqeNzP7e6zmT8zsiGTXFGddx5nZ2nLb6ubKlktCXS3N7H0zmxn7X7y6kmWSt83cPaMvwElA3djtPwF/qmSZPGAesD+QD0wDOiS5roOBA4FxQNEOllsAFKRwe1VbV0Tb68/ADbHbN1T2e4w9tz4F26janx+4Ang4drsv8FKa1HUR8ECq/p5i6+wGHAFMr+L5U4G3COft/BiYmCZ1HQeMSuW2iq23KXBE7PYewOeV/B6Tts0yfk/f3d/2cBYwwATCuQAVxTOURG3XNdPda3CSWXLEWVfKtxfbDuXxNHBWkte3IzUZeiTqulLO3T8k9NqrSi/gGQ8mAHubWdM0qCsS7r7c3afEbn8DzGT7kQqSts0yPvQr+CXh07Gi1A0HkTgH3jazybEzk9NBFNtrX3dfDuGfAmhSxXINzKzYzCaYWbI+GGoy9Egyxft7OSfWJPCKmbWs5PlUS+f/v6PNbJqZvWVmh6R65bFmwU7AxApPJW2bxdVlM2pm9i6wXyVP3eTur8eWuYlwLsDzlb1FJY/VuNtSPHXFoau7LzOzJsA7ZjYrtocSZV0p314JvE2r2PbaH3jPzD5193k1ra2Cmgw9kkzxrPNfwIvuvtHMLiN8Gzk+yXVVJ4ptFY8pQGt3X29mpwL/JJxAmhJmtjvwKvBrd19X8elKXlIr2ywjQt/dT9zR82bWHzgdOMFjDWIVxDUcRG3XFed7LItdrzCz1whf4WsU+rVQV8q3l5l9aWZN3X157Gvsiireo2x7zTezcYS9pNoO/ZoMPZJM1dbl7l+Xu/sosWHOI5aUv6eaKh+07v6mmT1oZgXunvQxecysHiHwn3f3f1SySNK2WcY375hZT+B64Ex3/66KxeIZSiLlzGw3M9uj7DbhoHSlPQ1SLIrtVX4oj/7Adt9IzKyhmdWP3S4AuhLO9q5tNRl6JJmqratCu++ZhPbiqI0EfhHrkfJjYG1ZU16UzGy/suMwZtaFkIdf7/hVtbJeI4xiMNPd/1rFYsnbZqk+cl3bF8IY/YuBqbFLWY+KZsCb5ZY7lXCUfB6hmSPZdfUmfFpvBL4ExlSsi9ALY1rsMiNd6opoezUCxgJzYtf7xB4vAh6L3f4J8Glse30KDEhiPdv9/MAQws4FQAPg5djf3yRg/2RvozjrujP2tzQNeB84KAU1vQgsBzbF/rYGAJcBl8WeN8JETPNiv7cqe7OluK7B5bbVBOAnKarrGEJTzSflcuvUVG0znZErIpJDMr55R0RE4qfQFxHJIQp9EZEcotAXEckhCn0RkRyi0BcRySEKfRGRHKLQFxHJIf8fc7QYrtl7VOsAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "simu_loi_wigner(int(1e5))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Méthode de rejet pour les lois beta\n", "\n", "Le but de cet exercice est d'utiliser la méthode du rejet pour simuler selon la loi Beta de paramètres $\\alpha>1$ and $\\beta>1$ définie sur l'intervalle $[0, 1]$ par sa densité \n", "$$\n", "f(x; \\alpha, \\beta) =\\frac{ x^{\\alpha-1}(1 - x)^{\\beta-1}}{B(\\al,\\bet)},\\qquad\\qquad \\text{ où } \\qquad \\;\\;B(\\al,\\bet):=\\f{\\Ga(\\al)\\Ga(\\bet)}{\\Ga(\\al+\\bet)}\n", "$$\n", "en utilisant une méthode de rejet. Rappelons que l'espérance, le mode est la variance de cette loi sont respectivement :\n", "\\begin{eqnarray}\n", "\t\\E[X] &= &\\frac{\\alpha}{\\alpha + \\beta}\\\\\n", "\t\\argmax_{x\\in[0,1]} \\ f(x; \\alpha, \\beta)& =& \\frac{\\alpha-1}{\\alpha + \\beta-2}\\\\\n", "\t\\Var(X) &=& \\frac{\\alpha\\beta}{(\\alpha +\\beta)^2(\\alpha+\\beta+1)}\n", "\\end{eqnarray}\n", "\n", "1. Comme la distribution de la loi est définie sur l'intervalle [0,1], on se propose tout d'abord d'utiliser comme loi $g$ la loi uniforme sur cet intervalle. \n", " + Prouver que cela est possible et coder l'algorithme de rejet pour la loi Beta. Simuler un grand nombre de variables aléatoires et comparer avec la densité cible.\n", "\t + Donner une condition sur les paramètres de la loi cible Beta pour laquelle la méthode est la plus efficace.\n", "2. Lorsque les paramètres satisfont $$\\alpha= \\beta,$$ la densité de la loi Beta ressemble à  une cloche. Pour cette raison, on peut dans ce cas utiliser la loi normale comme loi de proposition. \n", " + Quelles valeurs de moyenne et variance $\\mu$, $\\sigma^2$ semblent raisonnables pour cette loi normale ?\n", " + Pour ces valeurs, comment choisir la constante de domination $c$ de la méthode de rejet en fonction de $\\alpha = \\beta$ ?\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Correction avec la loi uniforme\n", "On pose $g$ la densité de la loi uniforme sur [0,1]. Soit \n", "$$c = f\\left(\\frac{\\alpha-1}{\\alpha + \\beta-2},\\alpha,\\beta\\right),$$ la densité de la loi Beta est dominée par $c\\times g$." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XucjeXawPHfNQwj55wqg0EoxuQwieRQJCS0o7Q70IGwS7vDlI677Np7l3a1e18dd1QipVey0S5ElENGDoUIoUE51EzEMGOu9497mQZzWGbWWs9aa67v5/N81qy1nlnP9ayZudY993Pf1y2qijHGmOgS43UAxhhjAs+SuzHGRCFL7sYYE4UsuRtjTBSy5G6MMVHIkrsxxkQhS+7GGBOFLLkbY0wUsuRujDFRqKxXB65Zs6YmJCR4dXhjjIlIK1as2KuqtYraz7PknpCQQGpqqleHN8aYiCQi2/zZz7pljDEmCllyN8aYKGTJ3RhjopBnfe75ycrKIi0tjczMTK9DiSpxcXHEx8cTGxvrdSjGmBAJq+SelpZG5cqVSUhIQES8DicqqCr79u0jLS2Nhg0beh2OMSZEwqpbJjMzkxo1alhiDyARoUaNGvbfkDGlTFgld8ASexDYe2pM6RN2yd0YY0zJWXIvxGOPPcYzzzwT0Nfs3bs36enppKen8+KLL57y96ekpNCiRQtSUlICGpcxJrpYcg+x2bNnU61atWIn91deeYWvvvqKsWPHBiE647WEBBAperPKHaYoltxP8OSTT9KsWTO6d+/Ohg0bch/fvHkzPXv2pG3btnTq1Ilvv/0WgCFDhjBq1CguvPBCGjVqxPvvvw/Arl276Ny5M61atSIxMZFFixYBruzC3r17GT16NJs3b6ZVq1akpKRwww038OGHH+Ye77rrrmPGjBnHxda3b19+++03LrjgAt59991gvxWlmldJdts2UC162+bXBHRTmoXVUMjj/PnPsGpVYF+zVSt4/vkCn16xYgVTpkxh5cqVZGdn06ZNG9q2bQvAsGHDePnll2nSpAnLli1j5MiRfPrpp4BL5J9//jnffvstffv2ZcCAAUyePJnLLruMhx56iKNHj3Lw4MHjjvWPf/yDb775hlW+c/zss8947rnn6NevHxkZGSxevJg333zzuO+ZMWMGlSpVyv0eEzzHkmxRin2tOjsbfv4Z9u07bruLdHjyEBw6BJmZv9+quoPFxEBMDC8hcEcsVKoElSu722Nf16wJZ5wBdepA9eolCNJEsvBN7h5YtGgRV155JaeddhrgWsoABw4cYPHixQwcODB338OHD+d+3b9/f2JiYmjevDk//fQTAOeffz4333wzWVlZ9O/fn1atWhV67C5duvCnP/2J3bt3M23aNK666irKlrUfT8Q6fBg2bYKNG90nxfbt8MMPv9/++GO+nx7PAjyMS+IVKkBcnNtiYiAnx31PTg5XkgNvH4EDB9wHRQGOEMtP1GE79dlCIzbTOPd2A83YR00AGjSArVuD8k4Yj4Rv9iikhR1M+Q0bzMnJoVq1agW2mMuXL5/7tfr+YDt37szChQuZNWsWN9xwAykpKdx4442FHvuGG25g0qRJTJkyhfHjx5fgLEyoxHIEVq6FlSth3Tr49lvYsAG2bHHJ+JgKFaB+fahXD3r1gvh418KuWRNq1MjdqjWsRvrh0yA2ttAW9wUJv3fNxHKEyuynEgeowq/UZC9n8CN1+ImmVX5i5JW7qLd9Ox03L4QfJh3/oXLmmZCUxNMfJ8HbSdC2LTRr5j5MTEQL3+Tugc6dOzNkyBBGjx5NdnY2//nPf7jtttuoUqUKDRs2ZOrUqQwcOBBVZc2aNZx33nkFvta2bduoW7cuQ4cO5bfffuOrr746LrlXrlyZ/fv3H/c9Q4YMoV27dpxxxhm0aNEiaOdpiunIEZfEly93t199xQHWQpss93xcHDRtCq1bw7XXuiTZrBk0bAinn+5X90gGQLmiQzm+lV0OqOHbinD4sPtU2LwZ1q+HNWtgzRru5F9wwxG3T/Xq0KEDXHih29q3dx9OJqIUmdxFZDzQB9itqon5PH8dcL/v7gFghKquDmiUIdKmTRuuueYaWrVqRYMGDejUqVPuc5MmTWLEiBE88cQTZGVlMWjQoEKT+4IFCxg7diyxsbFUqlSJt95667jna9SoQceOHUlMTKRXr16MHTuWOnXqcO6559K/f/+gnaM5BenpsGQJfP45fPEFfPml6wMH1+Ju04Z/0pMH3m3tEnrjxgW2eBMS/LsI2qBB4MLPV/ny7gOoaVP3H4RPJcki6+sN7hwXL3bnPXv279/TuTP06OG2li2tHz8SqGqhG9AZaAN8U8DzFwLVfV/3ApYV9ZqqStu2bfVE69atO+mx0uS3337TRo0aaXp6esBfO5rf2wYN/Blf4vYr1KFDqnPnqt5/v6bSRlXEfWOZMqrJyap33qk6darq9u2qOTmq6p72h7/7eSW/97AaP2tvZuo/uUu/pkXuEzs4U8cxQrsxR8uQVfz32xQLkKp+5NgiW+6qulBEEgp5fnGeu0uB+GJ+zpRqc+fO5eabb+buu++matWqXocTUYo9skXV9ZN/9BHMmQMLF7qRKbGxHKADPPYYdOoE7dpBxYrBCD1s5H8xtTpwuW8DduyAOXM4a9YsRs5+k5EHX3LdTf36waBB0K0blCmT+93WuPdWoPvcbwE+CvBrlgrdu3dn+/btXocR/Y4edd0OH37otk2b3OPNm8Ntt8Gll0KXLnStXAl9tOiXa9DAvyQW9O6WUKhbF4YMcdvBg/DJJ/B//wfTpsGECe75G290zzdtekrvjY3UCbyAJXcRuRiX3C8qZJ9hwDCA+vXrB+rQxhQuKws++YR3K/4fu8vOpDZ7OEIsn3IJ07mHmfRhx7p4WAf8y32Lv8m41Cal006D/v3ddvgw/Oc/8MYb8NRT8Pe/w4UXsvXJkTBwIJQr/AqxtfCDIyDjnUQkCfg30E9V9xW0n6q+qqrJqppcq1aRi3cbU3w5OfDZZzB8uJvQ06cPV5edRu0/Xgrvvku5jL301P/ysg4nTeNP6jEutUm7OMqXhwEDYOZMSEuDp5+GvXvh+uvdp+SYMeCb/2FCp8TJXUTqA9OAG1R1Y8lDMqYEVq+GlBSXVLp2hYkT4bLLXMty926YNAmuvhqqVPE60uh05pnu/V+/Hv77XzeK6C9/cWP8hw514/9NSBSZ3EXkHWAJ0ExE0kTkFhEZLiLDfbs8ihtg+6KIrBKR1CDGa8xJqpABL78MycmuxMS//uWSyuTJLqFPngx9+hTZPWACKCbGfajOnu0mdd16q/ugbdrU9clvtHZg0PkzpCYYmz9DIf0d4ubv5s/QrO+//15btGhx0uOPPPKIzpkzp8Dv++CDD3Tt2rVFH8AjUTcUMidHdeFC1Rtv1N+o4H7ASUmqL7ygunev19GZ/OzcqXrXXaoVKqjGxKhef73q1q1hP0w03ODnUMiwnmPsb4U8f7eSVNIbM2YM3bt3L/D56dOns27duuIfII/sQmqFlHoHDsCLL7rRLZ07wwcf8BY3ulmjq1bBHXe4qfwm/Jx5Jjz7rLugce+98P770KwZT3Ef/PKL19FFnbBO7l45evQoQ4cOpUWLFvTo0YNDhw4xZMiQ3HK+o0ePpnnz5iQlJXHvvfeyePFiZsyYQUpKCq1atWLz5s2sWrWK9u3bk5SUxJVXXskvvl/e5cuXk5SURIcOHUhJSSEx0U36feONNxg4cCBXXHEFPXr04MCBA3Tr1o02bdrQsmXL3HLAW7du5ZxzzuHWW28lMTGR6667jrlz59KxY0eaNGnCl19+6c2bFmxbtsA997iaLH/6k6uAOGEC7NrFCHxdMjbsIjLUru1G1WzcCIMGcS/PwNlnw//8T6FF0Mwp8qd5H4zNn26ZQP+75s/rff/991qmTBlduXKlqqoOHDhQJ06cqIMHD9apU6fqvn37tGnTpprjm6H4yy+/qKrmPn9My5YtdcGCBarqunTuvPNOVVVt0aKFfvHFF6qqev/99+d2AU2YMEHr1q2r+/btU1XVrKwszcjIUFXVPXv2aOPGjTUnJyc3vjVr1ujRo0e1TZs2etNNN2lOTo5Onz5d+/Xrl+95RWS3TE6O6qefqvbt62aLli2rOmiQ6pIluTNEVQP/e2JCK4lVqt27ux9k69aqy5Z5HVJYIxq6ZbzSsGHD3BK9bdu2ZWuecXFVqlQhLi6OW2+9lWnTpuWWB84rIyOD9PR0unTpAsDgwYNZuHAh6enp7N+/nwsvvBCAP/7xj8d936WXXsrpp58OuA/dBx98kKSkJLp3786OHTtyywk3bNiQli1bEhMTQ4sWLejWrRsiQsuWLY+LNWLl5MAHH8AFF8All7g6Jw895P6df+cdV8jKWulRYw3nuQlR773nhky2bw8jRlhXTQlZcs9H3hK+ZcqUOa4PvGzZsnz55ZdcddVVTJ8+nZ49e/r9uu5Dt2AV80xxnzRpEnv27GHFihWsWrWKOnXqkJmZeVJ8MTExufdjYmIiu78+K8tNhGnRAv7wB7at3MdtvEzcnu3IE39F4uvmuxpSVMz+LMUaNACJEeTqgVTZuZ7ndRRHX36VHacn0ks+sqUFi8mS+yk6cOAAGRkZ9O7dm+effz63xnveEr5Vq1alevXquUvrTZw4kS5dulC9enUqV67M0qVLAZgyZUqBx8nIyKB27drExsYyf/58tkXzumqZma6/tXFjuOkmN2TxnXdonL2BV/Q2MjWu0Avl0fDPSmm2devvP8tftQp/1ucpk/oldVtU5yN6o8NuQ3/db0sLnqKwrufub22KU3m9ktq/fz/9+vUjMzMTVeW5554DYNCgQQwdOpQXXniB999/nzfffJPhw4dz8OBBGjVqxIQJEwB4/fXXGTp0KBUrVqRr164FFgm77rrruOKKK0hOTqZVq1acc845JQ8+3Bw5AuPHwxNPuKJUF13kxqv36gUiHL3W6wCNZ9q2hdRUePRReOYZmDOHC5gMtPc6ssjhT8d8MLbSWvJ3//79uV///e9/11GjRoXkuOH03jaqn6VDGK9bSFAF/ZwL9WLmWclYk79Fi1QTEvQIZVWfffa4i+mlEXZBNTzNmjWLVq1akZiYyKJFi3j44Ye9Dil0cnJgyhQ+2t6cCdxMw7Y14KOP6JjzOZ/qJdbdYvJ30UWwciUz6QN33w1XXeUWUjGFCutumWh0zTXXcM0113gdRugtWuTGqS9fTiYtYfp06NvXRr0Y/1Srxh+Yhj77PNx3n+u2mTHDXXw3+Qq7lrsWMaLEnDpP39ONG+EPf3CzSXfuhDfeoBWr3AIPltjNKRG46y63qMrBg26d149s+YiChFVyj4uLY9++fZbgA0hV2bdvH3FxcaE98L59MGqUa1nNmeMumm7cCIMHY72BpkQ6dHBrvTZu7ArCPf+8f0txlTJh1S0THx9PWloae/bs8TqUqBIXF0d8fIhWPzx6FP79b3jwQdcvOnSoW67ujDNCc3xTOtSr5xYuv+EG15rfuNENp82zzF9pF1bJPTY2loYNG3odhimuZctc3ZcVK6BLF/jf/wVf7RxjAq5iRVd8bPRoGDvW/bc4caKVdvax/49Nye3d61ro7duzc8VOrmUy8tl8pGWizSg1wRUT41Z+GjvWlS+44gr47TevowoLltxN8anC66+7BRjeeIOx3MtZv27gHb0WVbEZpSagjk1qzHdLuZebGM/RT+ayosalkJHhdbies+RuimfTJujWza2w07IlrFrFfYyFypW9jsxEqbxlCvLbJuhNlJn2PkmHl7tZzr5yIKWVJXdTpISE31tIZSWbe+UZDjZJImP+CobxCjEL5yOJLay7xXjvyiu5mvfcaJrevd3iLqWUJXdTpNwVsVatJrtte54hhdP6XkrVtHW8qsPI0RjrbjFhYzpXutLQS5bA5ZeX2j54S+6mSGXIhr/+1a129MMP8O67boZp3bpeh2ZM/gYOhLffdsMlr7mmVK7wZMndFO677/ici1x1voEDYd06uPpqm11qwt+gQW693VmzYPjwUjfRKazGuZswogqvvAL33EMzXH11Bg3yOipjTs1tt0FampshXbcuPP641xGFjLXczcl27XJ9lSNGQMeOJPKNJXYTucaMcYvAjBnjFlUvJSy5m+PNmuWGNi5Y4KZz//e/7MT61k0EE3H/hXbv7rpnlizxOqKQKDK5i8h4EdktIt8U8LyIyAsisklE1ohIm8CHaYIh7xDHcnKEsZICffqwel9dzjn0FXLH7UiZGBviaCJfbKwbCBAf76qU7tjhdURB50/L/Q2gsFWgewFNfNsw4KWSh2VCIXeI4/dbOXJBZ1J4BkaO5LxDy/hWz7EZpSYiFTiTtcbptNgyg/0/HmBZ/B9o1iDT61CDqsjkrqoLgZ8L2aUf8JZvBailQDUROTNQAZogmzYNWreG9eth6lQYNw5CXR7YmAAqbCbrWm1B5Q8mcgFfcvf2O70ONagC0edeF/ghz/0032MmnGVl8U98S5Y1aQIrV8KAAV5HZUzw9e8Po0dzG6+6rpooFYjknt+A53wHlIrIMBFJFZFUq9nuod274dJLuZvn4I473ESPRo28jsqY0BkzhsV0cNVMN2/2OpqgCERyTwPq5bkfD+zMb0dVfVVVk1U1uVatWgE4tDlly5e79SeXLeN6JsILL1j9a1P6xMZyLe+4xT2uuQaOHPE6ooALRHKfAdzoGzXTHshQ1V0BeF0TaOPHQ6dO7hd68WImcb3XERnjme00cH8TK1a4SU5Rxp+hkO8AS4BmIpImIreIyHARGe7bZTawBdgEvAaMDFq0pniysmDkSLjlFpfcU1PdRVRjSrsrr4TBg+Fvf3NJPoqIV4tRJycna2pqqifHLlV++cXVhJk3D1JS3C9xWVd1QqTUldswJlfu7396ulvIvXp1l+DLl/c6tEKJyApVTS5qP5uhGs02b3YrxS9c6KZdP/10bmI3xvhUq+YWdV+71i3mHiUsuUerRYvgggtgzx6YOxeGDPE6ImPCV69ecPPNbi3Wr7/2OpqAsOQejd56yy2BV6MGLFsGnTt7HZEx4e/pp13XzIgRkJPjdTQlZsk9mqi6q/6DB7sLp0uXwtlnex2VMZGhRg146in44gvXQIpwltyjxdGjbkLSI4/A9dfDRx+5VogxJl/51aCJuWUIi+nA7pvu43T5GRFXYC8SWXKPBpmZrt76uHFw773w5ps2McmYIuRXgyZHY7hw5YvUjtnHzyMfQdUV2ItEltwjRN7yvHm3qpLB/Aq94P33uZt/Is+MRcrE5F8V74TNSvkak49WrVzd91degQ0bvI6m2Cy5R4jc8rx5t10/knFeFy4u+zm8/TbP6t0FVsPLb7NSvsYU4C9/gdNOg9GjvY6k2Cy5R6q0NOjSBb77DmbOhOuu8zoiY6JH7dpw//0wfTod+dzraIrFknsk2rrVDW/ctQs++QQuu8zriIyJPnfdBWedxVhSInIqtyX3SPPdd26YY3q6KynQsaPXERkTnU47DR5/nA4sdWsLRxhL7pFk7VrXYj98GObPh/PP9zoiY6Lb4MF8TwKMGRNxrXdL7hEika+ha1c3zOWzz+C887wOyZjoFxvL33jQrYPw8cdeR3NKLLlHgnXrmEc3V61u4UI491yvIzKm1HiTwVC/Pjz+eES13i25h7uNG6FbN45SxnXFWDkBY0Iqi3LwwAOunMfcuV6H4zdL7uFs82a45BI4epRuzHMLWRtjQu+mm+Css1xxsQhhyT1cbdvmEntmJsybx3qaex2RMaVX+fKudtPcubBmjdfR+MWSezjatcsl9l9/hTlzoGVLryMyxgwb5oZHPv+815H4xZJ7uElPd5OSfvrJXZ23tU6NCQ+nn+7KaU+a5P4+w5wl93By6BBccQV8+y1Mnw7t2nkdkTGlXt7SwM1euhOOHOHxM17MtxhfOJUHtuQeLrKz4Zpr3EIBkyZB9+5eR2SM4fjSwBu0GfTqxV/O+jealX1SMb5wKg9syT0cqMKtt8J//uNqsg8c6HVExpiCDBsGO3fC7NleR1IoS+7h4MEH3QIbjz/u1m80xoSvyy+HM86A117zOpJCWXL32muvwT/+4ZL6I494HY0xpiixsW7c++zZrvR2mPIruYtITxHZICKbROSk6vUiUl9E5ovIShFZIyK9Ax9qFJozh+xhI5hNL8q+9AISI7ZqkjGR4JZbICcHJkzwOpICFZncRaQMMA7oBTQHrhWRE2fUPAy8p6qtgUHAi4EONOp88w0MGMBaWtD713fJ1rK2apIxkaJxY+jWDcaPD9t6M/603NsBm1R1i6oeAaYA/U7YR4Eqvq+rAjsDF2IU+vFH129XsSKXMwsqV/Y6ImPMqbrxRtfqWrrU60jy5U9yrwv8kOd+mu+xvB4DrheRNGA2cEdAootGBw9C376wdy/MnMkO4r2OyBhTHP37Q1wcTJ7sdST58ie5Sz6Pnfh/yLXAG6oaD/QGJorISa8tIsNEJFVEUvfs2XPq0UY6VRg6FFJT4Z13oE0bryMyxhRXlSpu0uG777p5KmHGn+SeBtTLcz+ek7tdbgHeA1DVJUAcUPPEF1LVV1U1WVWTa9WqVbyII0RCwskXRe+JeRYmT+ZBfQLp19culBoT6f74R9izxy15GWb8Se7LgSYi0lBEyuEumM44YZ/tQDcAETkXl9xLYdP8d9u2nXBR9JM5/DPmPhgwgL/lPGAXSo2JBr16QdWqYdk1U2RyV9Vs4HbgY2A9blTMWhEZIyJ9fbvdAwwVkdXAO8AQ1TC9hOyFLVtcaYHmzd3QKcmvp8sYE3HKl3fX0GbODLuumbL+7KSqs3EXSvM+9mier9cBHQMbWpQ4cMBdeAFXDKxSJW/jMcYEVr9+MHEifP450NXraHLZDNVgUoXbboO1a2HKFDc21hgTXS67zLXgp0/3OpLjWHIPptdec31xY8ZAjx5eR2OMCYZKlVwV1w8/5OSBhN6x5B4k57EKRo1ySf2BB7wOxxgTTP37w9atJBE+S/BZcg+GX39lKgOhRg14+22IsbfZmKh2xRUgQt+TBhJ6x7JOoPkmKjXke9fPHuXj+Y0xQJ060KYNlzLH60hyWXIPtJdfhvfe42GegE6dvI7GGBMqPXrQgSVuYfswYMk9kNatg7vvhp49eZr7vI7GGBNKPXoQSzYsWOB1JIAl98A5cgSuu85dOZ8wAbW31pjSpUMHDlARPvnE60gAS+6B8+ijsGoV/PvfbgkuY0zpUr48C+hqyT2qfPYZPP20q/jY78RS98aY0mIe3eC772DHDq9DseReYunprmh/48bw7LNeR2OM8dAifIMoFi3yNhAsuZfcqFHuU3rSJKsbY0wpt4pWLg9Yco9ws2a5gkEPPQTt2nkdjTHGY0cpCx06WHKPaOnpMGwYJCa65G6MMeDmt3zzDfzyi6dhWHIvrpQUt9D1hAlQrpzX0RhjwkWnTm6m+hdfeBqGJffimDPHDXlMSYHkZK+jMcaEkwsugNhYX31371hyP0XN6x9ga4+hfEszKjz1l5PWST222dqoxpRSFSpAUhIsX+5pGH6txGR+d+sPj5Ig22HRIg51rOB1OMaYcHT++fDOO5CT41lVWGu5n4rVqxnFC+5CakdbVdAYU4DkZMjIgE2bPAvBkru/cnJgxAh+5nT429+8jsYYE87OP9/dpqZ6FoIld3+NHw9LlpDCWDj9dK+jMcaEs+bNXd+7h/3ultz9sXcv3H8/dOrEW9zodTTGmHBXtiy0bm0t97B3//2uAP9LLwHidTTGmEiQlOQmM6k3i2Zbci/Kl1+6Lpm77oIWLbyOxhgTKVq2dDPZPaoQ6VdyF5GeIrJBRDaJyOgC9rlaRNaJyFoRmRzYMD2i6lZWqlMHHn7Y62iMMZGkZUt3+/XXnhy+yOQuImWAcUAvoDlwrYg0P2GfJsADQEdVbQH8OQixht7UqW4K8RNPQJUqXkdjjIkkiYnuNlyTO9AO2KSqW1T1CDAFOHFFiqHAOFX9BUBVdwc2TA9kZsJ998F558FNN3kdjTEm0lSvDvHxYZ3c6wI/5Lmf5nssr6ZAUxH5QkSWikjP/F5IRIaJSKqIpO7Zs6d4EQdJQsLx5QMeqPAcbNvGxaufQ8qWsbICxphT17KlZ8ndn/ID+Q0POfHyb1mgCdAViAcWiUiiqqYf902qrwKvAiQnJ3tzCbkA27bluaj944/Q5G/QrR/zp1/saVzGmAiWmAjz5sHRo1CmTEgP7U/LPQ2ol+d+PLAzn30+VNUsVf0e2IBL9pFpzBjXLTN2rNeRGGMiWdOmcOQI/PBD0fsGmD/JfTnQREQaikg5YBAw44R9pgMXA4hITVw3zZZABhoymzfDa6+5xa6bRO7nkzEmDJx9trv97ruQH7rI5K6q2cDtwMfAeuA9VV0rImNEpK9vt4+BfSKyDpgPpKjqvmAFHVSPPeZqMT/yiNeRGGMiTIMGx1+7q3uxayCO7PHdcY8nJAQ/FlGPZk8lJydrqodTc08kArrmazc6JiUFnnrK65CMMZFO1S2YPWwYPPdc7sMixZ+4KiIrVLXIVYJshmpeDz/sxrPff7/XkRhjooGI65oJx26Z0uIClsKMGa7VblUfjTGB0qSJJ3XdLbn7PMoYqFkT7rzT61CMMdHk7LNhyxbIzg7pYS25A6Sm0puPXB2ZSpW8jsYYE02aNIGsLNi+PaSHteQO8OST/EI1+NOfvI7EGBNtGjZ0t1u3hvSwlty//hqmT+df3GnFwYwxgVe/vrsN8UQmS+5PPgmVKvECo7yOxBgTjeLj3a0l9xDasAHeew9uv51fsBEyxpggiIuDWrUsuYfUM89A+fJulSVjjAmWevXsgmrI7N4NEyfCkCFQu7bX0Rhjoln9+tZyD5mXXoLDh+HP0bFolDEmjNWrZ8k9JA4dgnHjoE8faNbM62iMMdGuXj349Ve3hUjpTO6TJsGePXDPPV5HYowpDc46y93++GPIDln6krsqPPsstG4NXbp4HY0xpjQ4dl3vp59Cdkh/ltmLLnPnwvr18NZbrmKbMcYE27Hkvnt3yA5Z+lruL73kCoRdfbXXkRhjSos6ddytJfcg2bHDlfWsFQo6AAALxklEQVS9+WY3vt0YY0KhZk13a8k9SF5/3a1CPmyY15EYY0qTsmWhRo2Q9rmXnuSenQ2vvgqXXQaNG3sdjTGmtKld21ruQTFrluuWGT7c60iMMaVRnTqW3IPitdfcWNM+fbyOxBhTGlnLPbASEuAM+ZHsWf/l7ztvRGLLIsJJW4MGXkdqjIlqtWtbn3sgbdsGP/5zMmU5ygPrbkSVfLcQL5JijCltatWC9HS35F4IRH1yB+DNN+H88+Hcc72OxBhTWlWv7m4zMkJyOL+Su4j0FJENIrJJREYXst8AEVERSQ5ciCWTxGpYswYGD/Y6FGNMaVa1qrsNl+QuImWAcUAvoDlwrYg0z2e/ysAoYFmggyyJwbwJsbEwaJDXoRhjSrNq1dxtenpIDudPy70dsElVt6jqEWAK0C+f/f4KPA1kBjC+ksnO5jomwRVXuAkExhjjlWMt9zBK7nWBvFXm03yP5RKR1kA9VZ1Z2AuJyDARSRWR1D179pxysKdswQLqsBuuvz74xzLGmMIca7mHS7cMkF/pRM19UiQGeA4osji6qr6qqsmqmlyrVi3/oyyuqVPZTyXo2TP4xzLGmMKEYbdMGlAvz/14YGee+5WBRGCBiGwF2gMzPL+omp0N06Yxkz5QoYKnoRhjTDh2yywHmohIQxEpBwwCZhx7UlUzVLWmqiaoagKwFOirqqlBidhfn30Ge/cylYGehmGMMQBUqeJmTIZLt4yqZgO3Ax8D64H3VHWtiIwRkb7BDrDY3nsPKlbkI3p5HYkxxkBMDFSuHLKWu18rManqbGD2CY89WsC+XUseVgn5umTo04fMd61LxhgTJqpWDZ+We0T64gvYuxcGDPA6EmOM+V2lSvDbbyE5VHQm95kzoVw5V7vdGGPCRcWKcOBASA4Vvcm9a1fXv2WMMeGiYkVruRfbpk3w7bdWt90YE36sW6YEZs1yt5df7m0cxhhzImu5l8DMma60b6NGXkdijDHHsz73YjpwwE1esla7MSYcWcu9mD7/3K1y0qOH15EYY8zJ4uLg8OGQHCq6kvu8eW4IZMeOXkdijDEny03uWuSuJRVdyf3TT6FDBzjtNK8jMcaYk8XFQU4OZckO+qGiJ7n//DOsXAmXXOJ1JMYYk7/y5QGIC8GaRtGT3BcsAFXo1s3rSIwxJn9xce7GkvspmDfPXYk+/3yvIzHGmPz5knt5gn9RNXqS+/z50KmTu6BqjDHhyFrup+jnn2H9erjoIq8jMcaYglmfe+ESEtyCJse2njW+BOCShzsc97gINGjgbazGGJMrhC13vxbrCDfbtrlrp7n+sgSeiOHTjHZQybOwjDGmcL5u43IcCfqhIrLlfpIlS6BlS1dxzRhjwlWZMgDEkBP0Q0V+cs/JgWXLoH17ryMxxpjCiQCW3P2zfj38+qubmWqMMeEsxqVcsfIDfli2zN1ay90YE+58yd1a7v5Ytcr1tTdp4nUkxhhTOEvup2DVKkhKyn3TjDEmbIVbn7uI9BSRDSKySURG5/P83SKyTkTWiMg8EQnN6HJVWL0azjsvJIczxpgSCaeWu4iUAcYBvYDmwLUi0vyE3VYCyaqaBLwPPB3oQPO1dau7mNqqVUgOZ4wxJRJmF1TbAZtUdYuqHgGmAP3y7qCq81X1oO/uUiA+sGEWYPVqd2vJ3RgTCcKp5Q7UBX7Icz/N91hBbgE+KklQflu1yr1ZiYkhOZwxxpRICJO7P+UHJJ/H8v2fQkSuB5KBLgU8PwwYBlC/fn0/QyzE2rXQqJGtvGSMiQxhdkE1DaiX5348sPPEnUSkO/AQ0FdV8y1WrKqvqmqyqibXqlWrOPEeb8MGOOeckr+OMcaEQph1yywHmohIQxEpBwwCZuTdQURaA6/gEvvuwIeZj5wc+O47aNYsJIczxpgSC6fkrqrZwO3Ax8B64D1VXSsiY0Skr2+3sbh6jFNFZJWIzCjg5QJn+3bIzLTkboyJHCEcLeNXyV9VnQ3MPuGxR/N83T3AcRVtwwZ3a8ndGBMpwqnlHraOJXfrczfGRIowu6AanrZscTVlAnFh1hhjQsFa7n7Yvt2toSf5jdQ0xpgwZMndD9u3QyDGyhtjTKiEWfmB8GTJ3RgTaazPvXBxHII9eyy5G2Mii3XLFK7esVI39eoVvqMxxoQTS+6Fq8923xfWcjfGRBBL7oXLTe7WcjfGRBK7oFq42vjK15x5preBGGPMqbALqoWryV5X5rdCBa9DMcYY/1m3TOFqshdq1PA6DGOMOTWW3AtXg31Qs6bXYRhjzKmx5F44a7kbYyKSXVAtXA32WXI3xkQeu6BauJrstW4ZY0zksW6ZQhw9SjXSreVujIk8ltwLcfAgMair5W6MMZHEknshDh92t+XLexuHMcacKl+fu11Qzc+x5B4X520cxhhzqqzlXghruRtjIpUl90JYcjfGRCpL7oWw5G6MiVSW3Athyd0YE6nC7YKqiPQUkQ0isklERufzfHkRedf3/DIRSQh0oLksuRtjIlU4zVAVkTLAOKAX0By4VkSan7DbLcAvqno28BzwVKADzWXJ3RgTyWJiwiO5A+2ATaq6RVWPAFOAfifs0w940/f1+0A3Ed9HVKBZcjfGRLIwSu514diK1ACk+R7Ldx9VzQYygODUB7DkboyJZCFK7mX92Ce/FviJVwP82QcRGQYM8909ICIb/Dh+fmpy3nl7i/m9kaomYOcc/eycS4Wnaj4gTxX3nBv4s5M/yT0NyLsSdTyws4B90kSkLFAV+PnEF1LVV4FX/QmsMCKSqqrJJX2dSGLnXDrYOZcOoThnf7pllgNNRKShiJQDBgEzTthnBjDY9/UA4FNVDf5YH2OMMfkqsuWuqtkicjvwMVAGGK+qa0VkDJCqqjOA14GJIrIJ12IfFMygjTHGFM6fbhlUdTYw+4THHs3zdSYwMLChFarEXTsRyM65dLBzLh2Cfs5ivSfGGBN9Iq/8gDHGmCKFdXIPq7IHIeLHOd8tIutEZI2IzBMRv4ZFhbOizjnPfgNEREUk4kdW+HPOInK172e9VkQmhzrGQPPjd7u+iMwXkZW+3+/eXsQZKCIyXkR2i8g3BTwvIvKC7/1YIyJtAhqAqoblhrt4uxloBJQDVgPNT9hnJPCy7+tBwLtexx2Cc74YOM339YjScM6+/SoDC4GlQLLXcYfg59wEWAlU992v7XXcITjnV4ERvq+bA1u9jruE59wZaAN8U8DzvYGPcPOE2gPLAnn8cG65h1fZg9Ao8pxVdb6qHvTdXYqbdxDJ/Pk5A/wVeBrIDGVwQeLPOQ8FxqnqLwCqujvEMQaaP+esQBXf11U5eT5NRFHVheQz3yePfsBb6iwFqonImYE6fjgn9/AqexAa/pxzXrfgPvkjWZHnLCKtgXqqOjOUgQWRPz/npkBTEflCRJaKSM+QRRcc/pzzY8D1IpKGG513R2hC88yp/r2fEr+GQnokYGUPIojf5yMi1wPJQJegRhR8hZ6ziMTgKo0OCVVAIeDPz7ksrmumK+6/s0Uikqiq6UGOLVj8OedrgTdU9Z8i0gE3dyZRVYNfiMUbQc1f4dxyP5WyBxRW9iCC+HPOiEh34CGgr6oeDlFswVLUOVcGEoEFIrIV1zc5I8Ivqvr7u/2hqmap6vfABlyyj1T+nPMtwHsAqroEiMPVnYlWfv29F1c4J/fSWPagyHP2dVG8gkvskd4PC0Wcs6pmqGpNVU1Q1QTcdYa+qprqTbgB4c/v9nTcxXNEpCaum2ZLSKMMLH/OeTvQDUBEzsUl9z0hjTK0ZgA3+kbNtAcyVHVXwF7d6yvKRVxt7g1sxF1lf8j32BjcHze4H/5UYBPwJdDI65hDcM5zgZ+AVb5thtcxB/ucT9h3ARE+WsbPn7MAzwLrgK+BQV7HHIJzbg58gRtJswro4XXMJTzfd4BdQBaulX4LMBwYnudnPM73fnwd6N9rm6FqjDFRKJy7ZYwxxhSTJXdjjIlCltyNMSYKWXI3xpgoZMndGGOikCV3Y4yJQpbcjTEmCllyN8aYKPT/LlN9E2wiwUoAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "simu_loi_beta(int(1e5),1.5,1.3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "L'algorithme est le plus efficace lorsque la probabilité de rejet est faible, donc lorsque l'aire entre les courbes de $f$ et $c\\times g$ est faible. C'est le cas lorsque $\\alpha$ et $\\beta$ tendent tous les 2 vers $1$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Correction avec la loi normale lorsque $\\alpha = \\beta$\n", "On pose $g$ la densité de la loi normale centrée en $\\mu = \\frac{\\alpha}{\\alpha + \\beta}$ et de variance $\\sigma^2 = \\frac{\\alpha\\beta}{(\\alpha +\\beta)^2(\\alpha+\\beta+1)}$. Soit \n", "$$c = f\\left(\\frac{\\alpha-1}{\\alpha + \\beta-2},\\alpha,\\beta\\right),$$ la densité de la loi Beta est dominée par $c\\times g$." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd0FFUbx/HvTQiErtJrAgooNRSlgwqIIEURpElHQBRBOgQQBaUpGAsgvQqhN/UFQaoUpUrHBEIIICVCpISQ8rx/bIyUQBayu5PyfM7Zc5Kd2ZnfZJMnd+/cuWNEBKWUUimLm9UBlFJKOZ4Wd6WUSoG0uCulVAqkxV0ppVIgLe5KKZUCaXFXSqkUSIu7UkqlQFrclVIqBdLirpRSKVAaq3acPXt28fb2tmr3SimVLO3Zs+eyiORIaD3Liru3tze7d++2avdKKZUsGWNO27OedssopVQKpMVdKaVSIC3uSimVAlnW566USpoiIyMJCQnh1q1bVkdJ1Tw9PcmfPz8eHh6P9Xot7kqpu4SEhJA5c2a8vb0xxlgdJ1USEUJDQwkJCaFQoUKPtQ3tllFK3eXWrVtky5ZNC7uFjDFky5YtUZ+etLgrpe6jhd16iX0PtLgrdQcRISoqiujoaKujKJUoWtxVqhUTE8OmTZsYOnQotWrVomDBgqRJkwYPDw88PDzIlSsXVapUoVevXqxatUpPMFpk+PDhfP755w7dZv369bl69SpXr15l4sSJj/z6fv36UaJECfr16+fQXI6kJ1RVqnPhwgW++uor5syZQ0hICO7u7vj4+PDiiy9SoEABPD09iYyM5Pz58xw7dowpU6bg5+dH5syZadmyJb1796ZYsWJWH4ZKhB9//BGAoKAgJk6cSPfu3R/p9d999x2XLl0iXbp0zojnGCJiyaN8+fKilCtduXJFevfuLenTpxc3NzepX7++LFiwQP7555+Hvu727duydu1a6dChg6RLl04AeeuttyQoKMhFyV3ryJEjVkeQkSNHStGiRaVWrVrSokULGTdunIiIBAQESN26daVcuXJSrVo1OXr0qIiItGvXTnr06CGVK1eWQoUKyeLFi0VE5Ny5c1K9enUpU6aMlChRQrZs2SIiIl5eXnLp0iVp3ry5eHp6SpkyZaRv377y9ttvy4oVK+JytGrVSlauXHlXtoYNG4qbm5uUKVNGFi5c6NSfQ3zvBbBb7KixWtxVihcTEyPz58+XXLlyiZubm7Rr106OHz/+WNu6cOGC+Pr6Svr06cXT01NGjBghkZGRDk5srbsKSs+eIjVrOvbRs+dD9797924pWbKk3LhxQ8LCwuTpp5+OK+4vv/yynDhxQkREdu7cKS+99JKI2Ip706ZNJTo6Wg4fPixPP/20iIh8/vnnMnLkSBERiYqKivtH/m9xP3XqlJQoUSJu35s2bZLGjRuLiMjVq1fF29s73vc3Y8aMCf0YHSIxxV373FWKFhYWRosWLWjdujUFCxZk9+7dzJo1i6JFiz7W9nLmzMnIkSM5duwYDRs2ZOjQoVSvXp3AwEAHJ0+9tm7dyhtvvEGGDBnIkiULjRo1AuD69ets376dZs2a4ePjQ9euXTl//nzc615//XXc3NwoXrw4Fy5cAOD5559n5syZDB8+nIMHD5I5c+aH7rtmzZoEBARw8eJFFixYwJtvvkmaNMmz9zp5plbKDgcPHuSNN94gKCiIUaNG0a9fP9zd3R2y7YIFC7Jo0SL8/f3p2rUrZcuWZeHChdSvX98h208yvvzSkt3GNwwwJiaGJ554gv3798f7mjv7v20NXKhRowZbtmzhhx9+oE2bNvTr14+2bds+dN9t2rRh/vz5LFy4kBkzZiTiKKylLXeVIq1fv55q1aoRHh7O5s2bGThwoMMK+52aN2/OgQMHePrpp2nYsCETJkyIKyzq8dSoUYPly5cTHh7OtWvXWL16NQBZsmShUKFCLF68GLAV8AMHDjx0W6dPnyZnzpy88847dOrUib179961PHPmzFy7du2u59q3b8+Xsf/USpQo4ajDcjkt7irFmTt3LvXq1aNgwYLs3LmTqlWrOnV/Xl5ebNu2jcaNG9O7d28GDx6sBT4RypUrR/PmzfHx8eHNN9+kevXqccvmz5/P9OnTKVOmDCVKlGDlypUP3damTZvw8fGhbNmyLF26lJ49e961PFu2bFStWpWSJUvGDWvMlSsXzz33HB06dHD8wbmSPR3zznjoCVXlDDNmzBBjjLz88sty9epVl+47OjpaunbtKoD07NlTYmJiXLp/R0kKo2WsdOPGDSlcuLDLf3/ik5gTqtrnrlKMuXPn0qlTJ2rXrs2qVavw9PR06f7d3NyYNGkSnp6e+Pn5ERMTg5+fn17Kn4ysX7+ejh070rt3b7JmzWp1nETR4q5ShOXLl9O+fXteeuklVq5c6fLC/i9jDBMmTMDd3Z3x48eTM2dOhgwZYkkW9ehq165NcHCw1TEcQou7SvZ27dpFq1ateOGFF1i9ejXp06e3NI8xhnHjxnH58mWGDh1Kzpw56dKli6WZVOqjxV0la4GBgTRs2JC8efOyatUqMmTIYHUkwNZFM23aNEJDQ3n33XcpVKgQderUsTqWSkV0tIxKtsLCwnjttdeIjo7mp59+IkeOHFZHuouHhwcLFy6kRIkSNG/enICAAKsjqVREi7tKlkSEDh06EBAQwLJlyx77ilNny5QpEytXrsQYQ+PGje8bU62UsyRY3I0xBYwxG40xR40xh40xPeNZ50VjTJgxZn/sY5hz4ipl8/nnn7N8+XLGjRtHzZo1rY7zUIUKFWLRokUcP36c9u3b6xj4R5QUp/y116xZswgKCrLkPben5R4F9BGR54BKwHvGmOLxrLdVRHxiH584NKVSd9i0aRMDBw6kadOm9OrVy+o4dqlVqxZjxoxh2bJlTJo0yeo4qd6PP/7IE0884bTifvbsWTp16kRwcDDbtm2jW7duDt9HQhIs7iJyXkT2xn59DTgK5HN2MKXiExoaSqtWrXjmmWeYPn16shpD/uGHH1KvXj169+6d4GXzqd2nn35KsWLFqF27NsePH497PjAwkFdffZXy5ctTvXp1jh07BtimDPjggw+oUqUKhQsXZsmSJQCcP3+eGjVq4OPjQ8mSJdm6dSsA3t7eXL58mYEDBxIYGIiPjw/9+vWjTZs2d1312rp1a1atWnVXtpiYGLp3706JEiVo0KAB9evXj9vfv/Lly8dnn33GjBkzWLhwYbz/0C9cuMAbb7xBmTJlKFOmDNu3b3fMDy/WI42WMcZ4A2WBXfEsrmyMOQCcA/qKyOFEp1PqDiJCt27duHz5Mj/88ANZsmSxOtIjcXNzY9asWfj4+NCiRQt2795NxowZrY71UL169XrgRF2Py8fHJ27ulvjs2bOHhQsXsm/fPqKioihXrhzly5cHoEuXLkyePJkiRYqwa9cuunfvzi+//ALYCvm2bds4duwYjRo1omnTpnz//ffUrVsXX19foqOjuXnz5l37Gj16NIcOHYo7xs2bNzNhwgQaN25MWFgY27dvZ/bs2Xe9ZtmyZQQFBXHw4EEuXrzIc889R8eOHe9a59y5c3z00Ud07NiRQoUK8d57791X4D/44ANq1qzJ8uXLiY6O5vr164/3A30Au4u7MSYTsBToJSL/3LN4L+AlIteNMfWBFUCReLbRBegCtln1lHoUc+fOZcmSJYwaNYqyZctaHeex5MyZk3nz5lG7dm369u2rXTTxuHPKXyDeKX//FREREff1g6b87dixI5GRkbz++uv4+Pg8dN81a9bkvffe4+LFiyxbtizeKX+3bdtGs2bNcHNzI3fu3Lz00kv3bSdv3rxMnTqVWbNmUb16dd5+++371vnll1+YM2cOAO7u7g6/Itau4m6M8cBW2OeLyLJ7l99Z7EXkR2PMRGNMdhG5fM96U4ApABUqVNCzSspuQUFBvP/++1SvXj1J37fSHi+//DK9e/fmiy++oGnTptSqVcvqSA/0sBa2MyXlKX8f5eRo+/bt7V7X0ewZLWOA6cBRERn/gHVyx66HMeaF2O2GOjKoSr1iYmLi/kjmzJnjlKl7XW3EiBEULVqUTp066fDIeyT1KX+rVavG0qVLiYmJ4cKFC2zatOmxjrNWrVpxn9yio6P55597O0QSx57RMlWBNsDLdwx1rG+M6WaM+fcUcFPgUGyf+1dAC9HxXspBpk2bxubNmxk/fjze3t5Wx3GI9OnTM3PmTIKDg5P9JxFHS+pT/r755pvkz5+fkiVL0rVrVypWrPhYXSp+fn5s3LiRUqVKUb58eQ4fdvBpSnumjnTGQ6f8VfY4e/asZM2aVV566aVkO4Xuw/Tp00cA+fnnn62OEken/E14yt9r166JiMjly5elcOHCcv78eadk0XuoqhSrR48eRERE8N133yWrYY/2GjFiBEWKFOHdd9/l1q1bVsdJ9davX8+zzz5Ljx49Htoab9CgAT4+PlSvXp2hQ4eSO3duF6a0j04cppKs5cuXs2zZMkaNGkWRIvcNvkoR0qdPz8SJE6lTpw6jR49m+PDhVkdK1eyd8vdx+9ldSVvuKkm6fv06PXr0oEyZMvTp08fqOE5Vu3ZtWrZsyahRo/jzzz+tjqNSCC3uKkn69NNPOXv2LJMmTcLDw8PqOE43fvx40qdPT/fu3XXuGeUQWtxVknPixAm++OIL2rVrR+XKla2O4xK5c+fm008/Zf369fj7+1sdR6UAWtxVkiIi9OzZk/Tp0zN69Gir47hUt27dKF++PP369bvvMnmlHpUWd5WkrF69mv/9738MHz48SY5AcCZ3d3e+/PJLQkJCGDdunNVx4nh7gzGOe9hzqUJQUBAlS5a87/lhw4axfv36B75uxYoVHDly5PEPNgXR4q6SjFu3bvHhhx9SvHhx3n//favjWKJatWq89dZbjBkzhjNnzlgdB4DTp0HEcY/Tpx8/yyeffELt2rUfuNyRxT0qKsoh27GKFneVZHz11VecPHkSPz+/VHES9UHGjh1LTEwMgwYNsjqKpaKjo3nnnXcoUaIEr7zyCuHh4bRv3z5uet2BAwdSvHhxSpcuTd++fdm+fTurVq2iX79++Pj4EBgYyP79+6lUqRKlS5fmjTfe4MqVKwD8/vvvlC5dmsqVK9OvX7+4TwmzZs2iWbNmNGzYkFdeeYXr169Tq1YtypUrR6lSpeKuiA0KCuLZZ5+lc+fOlCxZktatW7N+/XqqVq1KkSJF+O2336z5od3JniudnPHQK1TVnS5fvixZs2aV+vXrWx0lSfD19RVAduzY4fJ933tVJDh2+/Zs79SpU+Lu7i779u0TEZFmzZrJ3LlzpV27drJ48WIJDQ2VokWLxl21fOXKFRGRuOX/KlWqlGzatElERIYOHSo9e/YUEZESJUrIr7/+KiIiAwYMkBIlSoiIyMyZMyVfvnwSGhoqIiKRkZESFhYmIiKXLl2Sp59+WmJiYuLy/fHHHxIdHS3lypWTDh06SExMjKxYsUIaN26c2B+TiOgVqioFGDlyJNeuXWPs2LFWR0kSBg4cSJ48eejVq1eqHRpZqFChuCl6y5cvT1BQUNyyLFmy4OnpSefOnVm2bFnc9MB3CgsL4+rVq3G3YWzXrh1btmzh6tWrXLt2jSpVqgDQqlWru15Xp04dnnrqKcDW+B08eDClS5emdu3anD17Nm464UKFClGqVCnc3NwoUaIEtWrVwhhDqVKl7spqFS3uynInT57k22+/pWPHjvHOwpcaZcqUiU8//ZRdu3axbNl9s2ynCndO4evu7n5XH3iaNGn47bffePPNN1mxYgWvvvqq3dtN6J/lnTdQmT9/PpcuXWLPnj3s37+fXLlyxU0TcWc+Nze3uO/d3NySRH+9FndlucGDB+Ph4cHHH39sdZQkpW3bthQvXhxfX98kUSySkuvXrxMWFkb9+vX58ssv4+Z4v3MK36xZs/Lkk0/G3Vpv7ty51KxZkyeffJLMmTOzc+dOABYuXPjA/YSFhZEzZ048PDzYuHEjpxNzNtjFtLgrS/3222/4+/vTp08f8ubNa3WcJMXd3Z3PPvuM48ePM3PmTMtyeHk5diikl1fiM127do0GDRpQunRpatasyYQJEwBo0aIF48aNo2zZsgQGBjJ79mz69etH6dKl2b9/P8OGDQNg+vTpdOnShcqVKyMiD5wkrHXr1uzevZsKFSowf/58nn322cSHdxFjVX9ehQoVZPfu3ZbsWyUNIsKLL77IsWPHCAgIIHPmzFZHSnJEhGrVqhEUFMSff/4Zb9+yox09epTnnnvO6fux0vXr18mUKRNgu4/q+fPn8fPzszjV/eJ7L4wxe0SkQkKv1Za7ssxPP/3Eli1b+Oijj7SwP4AxhtGjR3Pu3Dm+/vprq+OkGD/88AM+Pj6ULFmSrVu3MmTIEKsjOZy23JUlRIQKFSpw9epVjh07lqrHtdujQYMG/Prrr5w8eZInn3zSqftKDS335EJb7irZWbFiBXv37mXYsGFa2O0watQowsLCGDNmjEv2l1qHXyYliX0PtLgrl4uOjmbo0KEUK1aM1q1bWx0nWShVqhRvv/02X331FX/99ZdT9+Xp6UloaKgWeAuJCKGhoXh6ej72NvROTMrlFi1axOHDh1m4cCFp0uivoL2GDRvG999/z5gxY+JGhzhD/vz5CQkJ4dKlS07bh0qYp6cn+fPnf+zXa5+7cqmoqCiKFy+Op6cn+/fvx81NPzw+io4dO7JgwQJOnjxJnjx5rI6jLKB97ipJmjdvHn/++ScjRozQwv4YhgwZQlRUVKqb6149Om25K5e5ffs2xYoVI3v27Pz2228YY6yOlCx17tyZefPmERgYSL58+ayOo1xMW+4qyZk5cyZBQUGMGDFCC3siDBkyhOjoaEaNGmV1FJWEaXFXLhEZGcmoUaOoWLEidevWtTpOsubt7U2HDh2YOnVqkrmhh0p6tLgrl5g3bx6nT59m6NCh2mp3AF9fX0REW+/qgbS4K6eLioris88+o1y5ctSvX9/qOCmCl5cXnTp1Ytq0adp6V/HS4q6czt/fn4CAAIYMGaKtdgcaNGgQIpKkbqatkg4dLaOcKiYmhpIlS+Lu7s6BAwd0+KODderUie+//56goCBy5cpldRzlAg4bLWOMKWCM2WiMOWqMOWyM6RnPOsYY85UxJsAY84cxptzjBlcpy9KlSzl69Ci+vr5a2J1g4MCB3L59m/Hjx1sdRSUxCbbcjTF5gDwistcYkxnYA7wuIkfuWKc+0AOoD1QE/ESk4sO2qy33lC8mJoayZcsSERHB4cOHcXd3tzpSitSyZUvWrFnD6dOn4+79qVIuh7XcReS8iOyN/foacBS498qJxsCc2Jtz7wSeiP2noFKxNWvW8McffzB48GAt7E40ePBgrl+/rvO9q7s80udkY4w3UBbYdc+ifMCdp+xDuP8fgEpFRIQRI0ZQqFCh++4urxyrVKlSNGrUCD8/v7j7hypld3E3xmQClgK9ROSfexfH85L7+nuMMV2MMbuNMbt1xrmUbe3atezevZtBgwbpzI8u4Ovry5UrV5g8ebLVUVQSYddoGWOMB7AGWCsi9525McZ8B2wSkQWx3x8HXhSR8w/apva5p2w1atTg1KlTBAYGkjZtWqvjpAp16tTh4MGDnDp1ivTp01sdRzmJI0fLGGA6cDS+wh5rFdA2dtRMJSDsYYVdpWzbt29n69at9O3bVwu7C/n6+nLhwgVmzJhhdRSVBNgzWqYasBU4CMTEPj0YKAggIpNj/wF8A7wK3AQ6iMhDm+Xack+5GjduzLZt2wgODiZjxoxWx0k1RIRq1aoREhJCQECA3r4whbK35Z5gZ6iIbCP+PvU71xHgPfvjqZTqyJEjrFq1imHDhmlhdzFjDL6+vrz22mvMnz+f9u3bWx1JWUivKlEONW7cONKnT0+PHj2sjpIq1atXjzJlyjB27FhiYmISfoFKsbS4K4c5c+YM8+bNo3PnzmTPnt3qOKmSMYb+/ftz9OhR1qxZY3UcZSEt7sphJkyYgIjQu3dvq6Okam+99RZeXl6MHTvW6ijKQjoAWTnE33//zZQpU2jRogXe3t5Wx7HP7dsQHAx//WV7/P237bnISBCBzJkha1Z44gnw9gYvL0iXzurUCUqTJg19+vThgw8+4Ndff6Vq1apWR1IW0OKuHGLixIncuHGD/v37Wx0lfpGRsHcvbNkCv/8Ohw/DiRMQFWX3JmIwhJCffZRlNxX4necJLlCNI8GZnBj88XTs2JGPP/6YMWPGsGrVKqvjKAtocVeJdvPmTfz8/Khfvz6lS5e2Os5/Ll2C1ath+XL45Re4eROAkxTiECU5TCOOU4xz5OUvchNKNiJIRxRpKFgA/th+HcLCIDQUgoJwO3WKgidOUHDvXhofXw0i3D7jAS9Xg7p1oWlTePppiw/aJmPGjLz//vt8/PHHHDlyhOLFi1sdSbmaiFjyKF++vKiU4ZtvvhFANm/ebMn+vbxEbP0oIukIl+YskLXUkSjcREBO4SVf8568m32RyPnzjtlpWJjIzz/LpCz9ZD+l4wLsoKL0wE+e4nJcJrBldLVLly5J+vTppX379q7fuXIaYLfYUWP1Zh0qUaKioihSpAi5c+dm+/btltxpyRiQMyHg5wczZtj6zgsWhLZtoUkT8PGxreRMwcHg7w/ffw/794OnJ7RqBe+/D2XL2jJa8Kf2wQcfMHnyZE6ePEn+/PldH0A5nL0XMWnLXSXK999/L4CsWLHCmgCHD8sM2ot4eIi4u4s0ayaydq1IVJQ1eUREDhwQ6dpVJEMGW7O9Th2pwjZLopw6dUrc3d2ld+/eluxfOR52tty1uKvHFhMTI2XKlJHnnntOoqOjXbvzoCCRtm1FjJEbpBfp0UPk5EnXZkjIlSsiY8eK5Mxp+1OrVUtk+3aXx2jVqpVkypRJ/v77b5fvWzmevcVdx7mrx7Z27VoOHDhAv379XHcLvStXoE8fKFrU1g3Sty8FCYavvoJChVyTwV5PPAH9+sGpU/TmCzh0CKpUsXXXBAe7LEb//v25fv06kyZNctk+VRJgz38AZzy05Z78vfjii5IvXz6JiIhw/s5iYkTmzrW1gt3cRDp2FAkOFpHYz59JnJeXSAauy3CGyU085Qbp5SM+knSEu+TEa926dSVnzpxy8+ZN5+xAuQzaclfOtGvXLjZt2kTv3r2dNq2vt7ftPGgxc5xf3GpBmzbsvFiIsjG7MTOmYwoWwBjbtUVJXVAQ3JCMfCQfkz7oGBmaN2I4H3OrmA+y7de48n76tHP2P2DAAC5evMicOXOcswOV5OhoGfVYmjRpwsaNGwkODiZz5sxO2YebiSHG7xsYMMA2+mTMGOjcGVzVBeRs69ZBly62Lpru3WH0aEzmTE4ZVSMiVKxYkStXrnDs2DG9p20y5rCbdSh1rxMnTrBixQq6d+/utMLOmTOs4xXo2RNq1YKjR22FMKUUdoBXXrH1w/foARMnQrlylGOPU3b174RiAQEBLFu2zCn7UElLCvpLUa7yxRdfkDZtWj744APn7GDpUihVikrshClTbFeZ5s7tnH1ZLVMm2/j8TZsgPJwdVIYvvgAnTNf7xhtvUKRIEcaMGYNVn9iV62hxV4/kr7/+Yvbs2bRv355cuXI5duORkdC7t+0y/mefpQwH4J13nH8BUlJQowYcOMAaGkDfvvDaa7ZpDxzI3d2dvn37smfPHjZu3OjQbaukR4u7eiRff/01t2/fpk+fPo7d8Nmz8NJLMGGCrZtiyxZOkjTmaXGZp57iTZbC5MmwcSNUqAAHDjh0F23btiVXrlw6HXAqoMVd2e369etMnDgx7uO9w+zYAeXK2S7bX7DANmY91d5Y20DXrrB1q+2TTOXKsHChw7bu6elJz549465RUCmXFndlt2nTpnH16lXHTuu7cKGtxZ4pE/z2G7Ro4bhtJ0NeXrZeKPPC8+Q6u4ct4RWgZUtGm4G4mRjbMmMbJvq4unXrRqZMmbT1nsJpcVd2iYyMZPz48VSvXp2KFSsmfoMi8Mkn0LIlPP887NoFOi0tQUH/XdJ0QXJRI2I9vPsuAxlDTPNWSPitRI+Hf/LJJ+nSpQv+/v4EBQU5KrpKYrS4K7ssWrSIM2fOOKTVXsTrNnPd2sJHHzGHNqTbth6TI3tcq/TfR3K4OMnp0qaFb7+FceNs0y288opt1stE+vDDDzHGMH78eAeEVEmRXsSkEiQi+Pj4EBUVxcGDBxM3j8yNG/yUqSn1+B+MGAG+vqljNIwjLFwI7dpBoUIUOv4TpyRxc+m0b9+eRYsWERwcrDc0T0b0IiblMOvWreOPP/5I/ARhV67AK6/wCutg6lQYMkQL+6No0cJ2VeuFC2yjmu3CrkTo168f4eHhfPvttw4KqJISbbmrBNWuXZujR49y6tSpx59H5vx5263ojh+n6e35LJGmjg2Zmhw6xF+lapM7R4yt2Pv4PPamGjZsyI4dOwgODiZDhgwODKmcRVvuyiH27NnDhg0b6NWr1+MX9jNnoHp1OHkSfviBpWhhT5SSJanBFtt8Oy+9ZDsZ/ZgGDBhAaGgoM2fOdGBAlRRocVcPNW7cOLJkyUKXLl0ebwNnztgK0KVLsH491K7t2ICp1J8UtY2Fz5bN9jPdsuWxtlO1alUqV67M559/TlRUlINTKitpcVcPdPLkSRYvXky3bt3ImjXro28gJOS/wr5uHVSq5PiQqZmXl62oFygA9evDr78+8ib+nVAsKCiIJUuWOCGksooWd/VAEyZMwN3dnZ49ez76i0NC4MUX/yvsjhgbr+6XNy/88gvkywf16tkuBHtEjRo1olixYjqhWAqTYHE3xswwxlw0xhx6wPIXjTFhxpj9sY9hjo+pXO3y5ctMnz6dt99+m7x58z7ai8+d+6+wr12rhd3Zcue2FfgcOWwnrffufaSXu7m50a9fP/bv38/69eudFFK5mj0t91nAqwmss1VEfGIfnyQ+lrLat99+S3h4OH379rX7Nd7e8JT5m0P5XuFa4AUq/bMWU7mSXpzkBHHTFPz7yJ8Pr5O/cPpqFkLL16GUOfhI0xS8/fbb5MmThzFjxjg1t3KdBIu7iGwBEn9JnEqcq/+aAAAgAElEQVQ2bt68yTfffEODBg0o/ghTAlw+fZ2/K71GybR/knnDSnZKpTvuDvrfQ694T7w7pymIu0WfeOEV8AvZ8qXnYM7aSECg3dMUpEuXjl69erFhwwb27HHODUOUazmqz72yMeaAMeYnY0wJB21TWWTWrFlcvnz50aYaiIhgGU1sfb7+/vDyy84LqB7s6adto5Kio+GVV8jJBbtf2rVrV7JkycK4ceOcGFC5iiOK+17AS0TKAF8DKx60ojGmizFmtzFm96VLlxywa+Vo0dHRfPHFF1SqVIlq1arZ+yJo04ZX+BmmTYPXX3duSPVwzz4LP/wAf/3FT9SDf/6x62VZs2alW7duLF68mMDAQCeHVM6W6OIuIv+IyPXYr38EPIwx8U5UISJTRKSCiFTIkSNHYnetnGDZsmWcPHmS/v37Y+yZGkAE3nsPFi+mN19Ahw7OD6kSVrEiLFlCKQ7CG29ARIRdL+vZsydp0qTRCcVSgEQXd2NMbhNbBYwxL8Ru07H3B1MuISKMGTOGIkWK0KhRI/teNHYsfPcdDBzIBHo7N6B6NPXq0YnptpE0bdrYPmElIG/evLRp04YZM2Zw8eJFF4RUzmLPUMgFwA6gmDEmxBjTyRjTzRjTLXaVpsAhY8wB4Cughehg2WRp48aN7Nmzhz59+uDu7p7wC/z9YeBA25zsn37q/IDqkc2lrW264MWLbfdmtUPfvn2JiIjgm2++cXI65Uw6cZiKU6dOHQ4dOsSpU6fw9PR8+Mq//gq1asELL8DPP0O6dBhj66VRSUfce9Kzp+32hRMnwrvvJvi6N954g82bNxMcHEymTJmcH1TZTScOU4/k999/Z/369fTu3Tvhwv7nn9C4MRQsCMuXQ7p0rgmpHt/48fDaa7abj69dm+Dq/fv358qVK0yfPt0F4ZQzaMtdAdCkSRM2btxIcHAwmTNnfvCKly/bbtp89artxtbPPBO3SFvuSc9d78m1a1Ctmm2Q/PbtUOLho5Zr1KjB6dOnCQgIwMPDw+lZlX205a7sdvToUZYvX87777//8MJ++zY0aWKb6XHlyrsKu0oGMmeGNWsgQwZo0AASOGHav39/goOD8ff3d1FA5Uha3BVjxowhffr0fPDBBw9f8YMPbNPMzpwJVaq4JpxyrAIFYPVquHDB1rUWHv7AVevXr0/x4sUZO3asTiiWDGlxT+WCg4OZP38+77zzDg+69sDbG941k+C77xjFQEyrlvfNF6NzxiQjFSrAvHmwcyd07frAvjQ3Nzf69+/PwYMH+d///ufikCqxtLincp9//jkAffr0eeA6Xqc3MynNB/DaawyKGhnvfDE6Z0zSdN8EY/8+3mzCUD6BuXPp6fbVAycYa9myJfnz52fs2LEuza0ST4t7Knbp0iWmTZtGmzZtKFiwYPwrBQWxhKa2/vX588Ge8e8qyYhvgrF/HyOifeH11/Fz70Oh0xvjfX3atGn58MMP2bRpEzt27HBteJUoWtxTMT8/P27dusWAAQPiX+H6dWjcGA8ibSdQH+duTCrpcnODOXOgaFEW0+yBH726dOlCtmzZ+FQvVEtWtLinUv/88w/ffPMNTZo0oVixYvevIGKbJ+bQIZrjD0WLuj6kcr7MmWHFCtIQZZuD5ubN+1bJlCkTH374IT/88AP79u2zIKR6HFrcU6lJkyYRFhbGoEGD4l9hwgRYsgRGj2YddV0bTrlW0aK0Zj4cOADvvBPvCdb333+frFmzaus9GdHingqFh4czYcIE6tSpQ/ny5e9fYetW6N/fNqb9Ee7EpJKvH3kNRoyA77+HL7+8b3nWrFnp0aMHS5cu5fDhwxYkVI9Ki3sqNGPGDC5cuBB/q/2vv+Ctt6BwYZgxwza0QqUOgwfb5uLv3982d9A9evbsScaMGRk1apQF4dSj0uKeykRERDB69GiqVavGiy++ePfCqCho3hzCwmDpUj2BmtoYY7tArWBB2+/BPTfUyZ49O++++y4LFiwgICDAopDKXlrcU5lZs2YREhLCsGHD7r8Zx+DBsGULTJkCpUpZE1BZ64knbOdaLl+G1q3vmwO+T58+eHh4MHr0aIsCKntpcU9Fbt++zWeffUalSpWoXbs23t7/XdTSxCyDceOYyLuYNm/rlaepzF0XO5UryzsRX8PPP/NRmpF3/S7kyZObtGnfYfbs2QQHB1sdWz2EFvdUZM6cOQQHB8e12k+fjr2g5fgJlmVuD88/T/dbE/TK01To3oudpsZ0hrZt+dh8jKxdd9eya9f6YYzRq1aTOJ3yN5WIjIykWLFiZM+enV27dmGMsU0HG34LKlWyzfS4d68209V/btyw3Yv1wgXYtw/y5wdsLfjOnd9h7ty5nDp1ijx58lgcNHXRKX/VXebPn8+pU6fu72vv29c2vnnOHC3s6m4ZM9pOrN+6ZTvBGhkZt2jgwIFERkbGzU2kkh4t7qlAVFQUI0eOpGzZsrz22mtxz7/BMvj2W+jd23aXHqXuVawYTJtmu7nHkCFxTz/99NO0bt2ayZMnc+HCBQsDqgfR4p4KLFiwgMDAwLtb7adPM51OtulfddyyepjmzW1TA48dC+vWxT09dOhQIiIiGDNmjIXh1INon3sKFx0dTfHixfH09GTfvn24ubnZPl7XrMk/Ow6RJWAfPP201TFVUhceDs8/D5cvk+vCAS5ILgA6dOjAwoULCQwMJG/evBaHTB20z10B4O/vz4kTJxg2bJitsAMMGwY7dvAOU7WwK/ukTw/+/hAWxhzaQkwMYGu9R0VF6bj3JEiLewoWFRXF8OHDKVWqFG+88YbtyXXrYPRo6NyZRTS3NqBKXkqUgC+/pC7rIPZEauHChenQoQPfffcdZ86csTigupMW9xRszpw5/Pnnn4wYMcLWav/rL2jTBooXBz8/q+Op5KhLF5bwJvj6wq5dAPj6+iIifPbZZxaHU3fSPvcUKiIigqJFi5IrVy7buHYRqFsXtm2D33+HkiVt49z1vsfqET1prnCloI/trlz79kHWrHTv3p1p06Zx4sQJvB90zz7lENrnnspNnz6d4OBgRo4caRshM2YMrF8PX30FJUtaHU8lY1d5EhYsgOBg6NYNRBg8eDBubm6MHDnS6ngqlhb3FOjmzZuMHDmS6tWrU6dOHdi923YStVkz6NzZ6ngqJahSBT7+GBYuhJkzyZ8/P127dmXWrFkEBgZanU6h3TIpjrc3nD79BdAX2EwGyrOXcmTgJqX5w9bqiuXlpfPGqEcX150XHQ116tj63vfu5XyWLBQuXJjmzZsza9Ysq2OmWNotk0qdPn2N7NlH88orryBSgxvd+lLM/EmBDbO5Ik/qhGAq0eJmkEzjTt6N8wi96cnuZ1vjlTcbt269x+zZczHmCNr1bq0Ei7sxZoYx5qIx5tADlhtjzFfGmABjzB/GmHKOj6ns58fly5dtfZ+rV8PkydCnD7z8stXBVApx5wyS5yQv2ZZOpQJ7uD1oOJcuDSRz5oy8/rovp09bnTR1s6flPgt49SHL6wFFYh9dgEmJj6UeR2hoKPA5jRs35vmCBaFTJyhTBvQkl3KmJk2gY0cYPZrsR44wYMAAVqxYAWy3OlmqlmBxF5EtwN8PWaUxMEdsdgJPGGN0DlAL2MYZX2PkiBG2P7Z//oH58yFdOqujqZTOz8923902bejVvj25c+cGBmDVOT3lmD73fMCdl6aFxD53H2NMF2PMbmPM7kv33J9RJU5QUBDffPMN0J6S27bBjz/aJnoqUcLqaCo1yJTJ1pA4e5aM/fvz0UcfAdtYs2aN1clSLUcUdxPPc/H+uxaRKSJSQUQq5MiRwwG7Vv8aMmQI7u7uFKaNrY+9bl3o0cPqWCo1qVjRNuT2++/plCEDUIRBgwYRfc99WJVrOKK4hwAF7vg+P3DOAdtVdtq3bx/z58+nV48eLKIPZMhgu4v9vTfAVsrZBg+GKlXw6NGD7HzI4cOHmTt3rtWpUiVHFPdVQNvYUTOVgDAROe+A7So79e/fn2zZsjHg9m3Ks9d2cwW99ZmyQpo0MHcuiLCE73m+QgWGDh1KeHi41clSHXuGQi4AdgDFjDEhxphOxphuxphusav8CJwEAoCpQHenpVX3WbduHevXr2doy5Zk9fNjKp3h9detjqVSs8KF4ZtvqMk2xpQpQ0hISOz5IOVKeoVqMhYdHU358uW5FhbG0eho0qZLR6aAfVyXTFZHU6mdCP5uLWieZhmvVazItoMH+fPPP8mZM6fVyZI9vUI1FZg7dy4HDhzg07x5SXvuHMybxw20sKskwBi6MRly5+aLs2e5ceMGw4YNszpVqqLFPZn6559/GDhwIJWLFKH59u22UQoVK1odS6k4V3kS5szh2dOnee+555g6dSp//PGH1bFSDS3uydRnn33GhQsX8Dt3DlOlim2UglJJiJcXmJdfYqz05aNDh3CPyUiZMh9ijNjmpol96Bw0zqHFPRkKCAhgwoQJtMuVi+eNsY1OSJPG6lhK3eXfOWj63xrBUz4+jM8YA/zCihWr7prATuegcQ4t7smEt/d/LZ0iRfoit2HUhQu0v/415unCccu8vKxOqtQ90qWD+fPpGhXFc5ky0adPHyIiIqxOleJpcU8mTp+2tXJ+/nk9sJJP3KLI07Qps2La6TS+KukrXhyPzz9n/PXrBAYG6tBIF9ChkMmEMRAZGYVP6dKEBwRwOHt2PA8ehGzZrI6mlH1EoH596q9bx7b06Tl24gR58+bVe/k+Ih0KmQJ98803HD56lM8jI/GcM0cLu0pejIGZM/kqa1Zu37xJ3969rU6UomlxTzbOMnTwYOoBr/fqBbVrWx1IqUeXOzfPzJzJQBEW+PuzYcMGqxOlWFrck4m0vEtUeDjfFCuGGTXK6jhKPb7GjRnQoQOFgfc6dABuW50oRdLingz876efuM1qfN3dKbx4MXh6Wh1JqURJ//XXfJM3L8fPnMETvVOYM2hxT+LCw8N5r21bigH9Ro+GUqWsjqRU4mXMSL0VK2gCxPAZQadOWZ0oxdHinsR91qcPJy9fpiPlSacnoFRK8vzzfNm3Lx5E0+PNN/WWfA6mxT0JO3bwIGMnT+bttGn5klXgpm+XSlkKjB5Ne7xZs28fi7791uo4KYpWiyQqOjqajvXqkUmEzydN4jx5rY6klOO5u7OGdTzv5kaP3r25fOGC1YlSDC3uSdTXvXqx4+xZvqpZk1wdO1odRymnOU0RZowYwdXISHrWrWt1nBRDi3sSFLh/P4O//ZYGGTLQavVqq+Mo5XQlBw3Ct3hxvj9wgDUTJlgdJ0XQ4m6xOycEsz1iaFC2Lh4inL25ALcsmXVCMJXyGcOgDRsomSYN3fr3J+zcOasTJXta3C3274Rg/z4mdujMMS4yvlEj9kojnRBMpRppc+dmhp8f56Oi6KdXYCeaFvckJGDLFvrNnEntrFnpuGSJ1XGUcrnnu3en7/PPM/XoUdbobfkSRYt7EhEVEcHbDRrgAcxYswbj4WF1JKVcwsvr7q7JCb+vpyiedBgxkhzmkN6x6TFpcU8iPm3QgF3XrjG5WzcKVKtmdRylXObfOzb9+7gtWVi6fAHXEKrkrEVMdIzesekxaHFPAnbOns2I9et5u2BBmk+caHUcpSxX8vXXGfX666y6eJHpbdpYHSdZ0pt1WCyTOU9uDy+iYmI48OefZC1UyOpISiUJMVFRvJIrFzv//pv9P/1EkXqv6k090Jt1JBs+VOdkZCRzJ0zQwq7UHdzSpGHWunWkNYbmTZqQhn+sjpSsaHG30MwuXfiVQIbWqEH1Hj2sjqNUkpO/fHlmDxzI3vBwnqem1XGSFS3uFvlj9Wq6T51KeZ5g2Lp1VsdRKslq+Nln9ClVih3sx3/QIKvjJBta3C1wLTSUZs2a8aQxnGMD7unSWR1JqSRt1NatlCE974wZw5/btlkdJ1mwq7gbY141xhw3xgQYYwbGs7y9MeaSMWZ/7KOz46OmDCJCl6pVCYiIYMGIEZynnNWRlEryPLJm5Qar8RCh2auvcvPaNasjJXkJFndjjDvwLVAPKA60NMYUj2dVfxHxiX1Mc3DOFOPLrl1ZePw4IypWpKavr9VxlEo2AqjF3F69+OPGDTpXqaI390iAPS33F4AAETkpIreBhUBj58ZKmdYuWEDfqVN5M2tWBv7yi9VxlEp26o8fz6c+Piw4dIix3bpZHSdJs6e45wPO3PF9SOxz93rTGPOHMWaJMaaAQ9KlIMePHKF527aUMobZGzbgliGD1ZGUSn6MYeCWLTTPnJlBU6bw44IFVidKsuwp7iae5+79PLQa8BaR0sB6YHa8GzKmizFmtzFm96VLlx4taTJz91S+V6laojppo6J4RsaSqUL5uGU6la9Sj8ZkzsyMtWvxMYaWbdty7MgRqyMlSfYU9xDgzpZ4fuCuyZZFJFREImK/nQqUj29DIjJFRCqISIUcOXI8Tt5k49+pfG/fjuSVcrUI42+WvvYaS6TvXfNo6FS+StnnzgnGMlapTHEZgWdUFJVL1MCYizrB2D3sKe6/A0WMMYWMMWmBFsCqO1cwxuS549tGwFHHRUy+RIR3WrZk3d69TM6fn+qLF1sdSalk694JxubFDGblSy8RQSjPP1uT69dv6ARjd0iwuItIFPA+sBZb0V4kIoeNMZ8YYxrFrvaBMeawMeYA8AHQ3lmBk5OhgwYxe+lShqdLR6eNGyF9eqsjKZVyGEOl5ctZmCcPe44d463GjYmKirI6VZKhE4c5iTGTgXfpDExZuhTTpInVkZRKmQ4dYnK5crwbGUnnjh2ZNmMaIvGdKkwZdOIwCy1atAhDd14DJn34oRZ2pZypZEm6zZnDYGDajBmAr46BR4u7w61cuZJWrVpRCYN/pUqkGTPG6khKpXwtWjCyZ0/eAWAUI0eOtDiQ9bS4O9BPP/1Es2bNqJAmDTPJQcalS0Fvl6eUS5hx45hctSqtcWfYsGGMHTvW6kiW0uLuIOvXr6dJkyaUSJeOn2JiaMcKyJvX6lhKpR4eHrgtXsxostMiUyYGDBiAn5+f1akso8XdAVavXk2DBg14JmNGfr5+nScnT2YXlayOpVTqkycPLVjCnPBwmuTMSa9evVJtC16LeyL5+/vTpEkTSuXNy6bQULK//z507Gh1LKVSrV+phsfkySy8eJEWRYsyYMAABg8enOpOsmpxT4Tp06fTsmVLqpQqxYZz58hWsyaMH291LKVU58549OrFvBMn6FqjBqNGjeK9994jJibG6mQuo8X9MYgIw4cPp3PnztStUYOfzp0jS+7csHixnkBVKqkYNw73V19l0q+/MqBlSyZNmkSrVq24deuW1clcQov7I7p9+zbt2rXj448/pkPr1qy6fJkMt27BDz9ACp8vR6lkJU0aWLgQU6QIo9euZeyAAfj7+/Pyyy9z8eJFq9M5nRb3RxAaGkqWLHWZO3cubnxM8/kX4PBxaoUtxZQscccskDrbo1JWuXOCMfNEVp45tprQv6HhmBVkYiY7duwnV66K5M172OqoTqXF3U579+6lfPnyRERsZ97cuUR3PkNd1uMxYwobpNZdExrpbI9KWefeCcYC5BmybV7Os2lPca3KVH7fuo7cuW9x/nwVVq1aleD2kist7naYOXMmVapUiT0Zs43WZ87AtGng6wsdOlgdTymVkBo1YP582LGDCl98wW87dgDP0LhxYwYMGJAiJxzT4v4QN2/epEuXLnTs2JFq1aqxZ88eOrMfBg+Gli3hk0+sjqiUslfTpuDnBytWUGDsWGAbXbt2ZezYsdSqVYvz589bndCxRMSSR/ny5SUp2717txQrVkwAGThwoERGRoosWiTRGJFXXxWJiLA6olLqcQwYIALiywgREZk3b55kyJBBcubMKStWrLA4XMKA3WJHjdXifo+oqCj59NNPJU2aNJIvXz7ZsGGDbcHatSIeHrKVqiI3blgbUin1+KKjRdq2tZW/CRNEROTw4cPi4+MjgLRv316uXr1qccgH0+L+GPbv3y8vvPCCAPLWW2/J33//bVvw668iGTKIlC4tWblibUilVOJFRsoimtpK4KRJIiISEREhvr6+4ubmJgUKFJC1a9daHDJ+WtwfwY0bN6R///7i7u4uOXLkkO+//15iYmJsC7dtE8mUSaRIEZHz520/MaVUsudBhEjDhrYyOGNG3PM7d+6UokWLCiDNmzeXs2fPWpjyflrc7RATEyNLliyRQoUKCSCdOnWS0NDQ/1bYuvW/wh4SIiKixV2pFAJEJDxcpE4dEWNEZs+OWxYeHi4ff/yxpEuXTjJnzixffvml3L5927qwd9DinoCdO3dK1apVBZASJUrI5s2b717h38JetKjIHf+5tbgrlTLE/S3fuCHy8su2JyZOvGudgIAAefXVVwWQokWLytKlS//7VG8RLe4PcOjQIXnrrbcEkFy5csmUKVMkMjJSvLwk7rKH+qyRG6SXoxSTPJy96/IkLy9LYiulHOyuhlp4+H9dNGPH3rVeTEyMrFy5Up577jkBpFKlSvc3Bl1Ii/s99u/fL02bNhVAMmbMKEOGDJF//vknbnncGz1rloi7u0j58iIXLrg0o1LKde5s0IFIGm7LApqLgIxksEDMXQ26yMhImTp1quTNm1cAqV69uvz0008ub8lrcRfbf9yff/5ZGjZsKIBkyZJFhgwZIpcvX75vXYgRGTPG9iOpXVvkjsKvlEoloqJE3nnHVgdathQJD7+vK/bGjRvi5+cn+fPnF0DKli0r/v7+LuuTT9XF/dq1azJx4sS4j1E5cuSQ4cOH/ze08V63bsl0Oth+HC1a6AVKSqVmMTEio0bZ6kHVqpKdi/GuFhERIdOnT5ciRYoIIHny5JGPPvpIQmIHXzhLqivu0dHRsmHDBmnbtq1kzJhRAClfvrzMnj1bwsPDH/zC8+dFKle2/SiGDrVd4KCUUosWiXh6SgCFRfbvf+BqUVFRsmbNGqlfv74YY8Td3V1ef/11WbJkycNrz2NKFcU9OjpafvvtNxkwYIAUKFAgruulc+fOsn379oT7wjZvFsmXTyRDBmnKokTnUUqlMDt3Sgh5RdKlE5kyxdaqf4jAwEDp37+/5M6dWwDJmjWrdOzYUdavX2+bwsQBUmxxj4iIkHXr1kn37t0lX758Aoi7u7vUq1dPFi5cKDdv3kx4I5GRtla6m5vIM8+I7NunQxyVUvHKwQXbeTgQefttkbCwBF8TGRkp69atk7Zt20qmTJkEkKeeekratGkjixcvvmswx6NKscV9xowZAkiGDBmkSZMmMnv27HhPkD7QwYMilSrZDr1du7gTp1rclVLxAbGdaB0+3NYgLFBA5H//s/v1N27ckCVLlkjbtm3lqaeeEkDee++9RORJocX98uXLsnLlSvta6PLfcKd0hMsIfOU2aeQS2aQ5C3T8ulIqQXc1/HbuFHnuuf8ah+fPP9K2IiMjZdOmTXLkyJFE5LGvuBvbuq5XoUIF2b17t9P3426iiZ41D4YNg+BgaNsWvvgCsmd3+r6VUsmftzecPv3f9+m4xUd8TB++IIJ0jGIQE/iQW6THy8v5d2EzxuwRkQoJrWfXzTqMMa8aY44bYwKMMQPjWZ7OGOMfu3yXMcb70SM7WEQEzJrFAcpA+/a2m1f/8gvMnq2FXSllt3tv23dLPBkko0h74jCZG9fiM3wJz1UIGTWaK6fDrI4bJ8HiboxxB74F6gHFgZbGmOL3rNYJuCIizwATgDGODmoXEdi3DwYMsN0lt0MHBAP+/vDbb/DSS5bEUkqlQEWKwIoVsHkzlCkDgwYRQn7brTc3bYKYGEvjJdgtY4ypDAwXkbqx3w8CEJFRd6yzNnadHcaYNMBfQA55yMYT3S0THQ1XrkBgIBw9Ctu2wcaNcPIkpEkDr74KPXti6tRCxDz+fpRSyh579zK9/Ld0yrwYrl2DbNng5ZehZk0oWRKKFrU9lzZtonZjb7dMGju2lQ84c8f3IUDFB60jIlHGmDAgG3DZvriPYNkybjZtQwa5edfTf/MkW6nOj/RnadSbhK7JDmtsDXillHK6cuUY4TWdHqe/phGrqBu6ljqLfyb/4sV3rRaOJ1Oz9uWDqyOcGseelnszoK6IdI79vg3wgoj0uGOdw7HrhMR+Hxi7Tug92+oCdIn9thhw/DFzZ8cZ/ziSNj3m1EGPOXVIzDF7iUiOhFayp+UeAhS44/v8wLkHrBMS2y2TFfj73g2JyBRgih37fChjzG57PpakJHrMqYMec+rgimO2Z7TM70ARY0whY0xaoAWw6p51VgHtYr9uCvzysP52pZRSzpVgyz22D/19YC3gDswQkcPGmE+wDaZfBUwH5hpjArC12Fs4M7RSSqmHs6dbBhH5EfjxnueG3fH1LaCZY6M9VKK7dpIhPebUQY85dXD6MVt2hapSSinnsesKVaWUUslLki7uyXLag0Sy45h7G2OOGGP+MMZsMMYk+5H8CR3zHes1NcaIMSbZj6yw55iNMW/FvteHjTHfuzqjo9nxu13QGLPRGLMv9ve7vhU5HcUYM8MYc9EYc+gBy40x5qvYn8cfxphyDg1gz+xiVjywnbwNBAoDaYEDQPF71ukOTI79ugXgb3VuFxzzS0CG2K/fTQ3HHLteZmALsBOoYHVuF7zPRYB9wJOx3+e0OrcLjnkK8G7s18WBIKtzJ/KYawDlgEMPWF4f+AkwQCVglyP3n5Rb7i8AASJyUkRuAwuBxves0xiYHfv1EqCWMSY5zzWQ4DGLyEaRuMtzd2K77iA5s+d9BhgBjAVuuTKck9hzzO8A34rIFQARuejijI5mzzELkCX266zcfz1NsiIiW4jnep87NAbmiM1O4AljTB5H7T8pF/f4pj3I96B1RCQK+Hfag+TKnmO+Uyds//mTswSP2RhTFiggImtcGcyJ7HmfiwJFjTG/GmN2GmNedVk657DnmIcDbxtjQrCNzutByvaof++PxIXfP9cAAAHBSURBVK6hkBaJrwV+79Aee9ZJTuw+HmPM20AFoKZTEznfQ4/ZGOOGbabR9q4K5AL2vM9psHXNvIjt09lWY0xJEbnq5GzOYs8xtwRmicgXsRMWzo09ZmunV3Qep9avpNxyf5RpD3jYtAfJiD3HjDGmNuALNBKRCBdlc5aEjjkzUBLYZIwJwtY3uSqZn1S193d7pYhEisgpbPMwFXFRPmew55g7AYsARGQH4IltDpaUyq6/98eVlIt7apz2IMFjju2i+A5bYU/u/bCQwDGLSJiIZBcRbxHxxnaeoZGIOP82Xs5jz+/2CmwnzzHGZMfWTXPSpSkdy55jDgZqARhjnsNW3C+5NKVrrQLaxo6aqQSEich5h23d6jPKCZxtrg+cwHaW3Tf2uU+w/XGD7c1fDAQAvwGFrc7sgmNeD1wA9sc+Vlmd2dnHfM+6m0jmo2XsfJ8NMB44AhwEWlid2QXHXBz4FdtImv3AK1ZnTuTxLgDOA5HYWumdgG5Atzve429jfx4HHf17rVeoKqVUCpSUu2WUUko9Ji3uSimVAmlxV0qpFEiLu1JKpUBa3JVSKgXS4q6UUimQFnellEqBtLgrpVQK9H/0AMu92nvP0gAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "simu_loi_beta2(int(1e5),3)" ] } ], "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.8.5" } }, "nbformat": 4, "nbformat_minor": 2 }