Skip to content

Instantly share code, notes, and snippets.

@exaucae
Last active January 10, 2021 20:23
Show Gist options
  • Save exaucae/858701a8eb91b9ae78e9ff33c2b7c580 to your computer and use it in GitHub Desktop.
Save exaucae/858701a8eb91b9ae78e9ff33c2b7c580 to your computer and use it in GitHub Desktop.
monte-carlo-simulation
"""
we tackle the acceptance rejection alogorithm with the beta law
It is well established that the rejection method allows to simulate data
of the normal law N (0; 1) according to a certain algorithm.
"""
# We illudtrate how to calculate a double integral with Monte Carlo's method
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Random variables Simulation\n",
"\n",
"It is well established that the rejection method allows to simulate data\n",
"of the normal law based on another well known random variable.\n",
" Here, we use the uniform law to simulate 1000 observations of X$\\sim$ N( $\\mu$ = 0,$\\sigma^2$ = 1).\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"\n",
"import matplotlib.pyplot as ploter\n",
"\n",
"from numpy import random, zeros, sqrt, pi, exp, log, linspace, arange, mean\n",
"from seaborn import histplot"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"\"\"\"\n",
"Generate random samples for the gaussian law N(0,1)\n",
"The algorithm is based on rejection method and the uniform law\n",
"\n",
"@param time: the number of samples generated\n",
"@return samples: an array of randomly generated samples \n",
"\n",
"\"\"\"\n",
"def get_gaussian_samples(time):\n",
" samples = zeros(time)\n",
" counter = 0\n",
" while (counter!=time-1):\n",
" U = random.uniform(size=3)\n",
" Y = -log(U[0])\n",
" if(U[1] <= exp(-(pow((Y-1),2)/2))):\n",
" Z=Y if U[2]<=1/2 else -Y\n",
" samples[counter] = Z \n",
" counter = counter+1\n",
" return samples"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmoAAAFNCAYAAACwk0NsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABoaUlEQVR4nO3dd3zV1f348dc7A8IIYYWRsJcsBxCWIC7EWaHWXUdbq7Vqa7G2tbXVarVfbb+K9Vdbd92KdXylTkDrQFkBkSnKEAhhz4SR+f79cT4XLjEJAXLvueP9fDzyuJ95P+988rmf+84553OOqCrGGGOMMSb2pPgOwBhjjDHGVM8SNWOMMcaYGGWJmjHGGGNMjLJEzRhjjDEmRlmiZowxxhgToyxRM8YYY4yJUQmTqInIIhE5yXccPonId0VkjYgUi8iAOmx/kogURCO2RCIifxSR53zHURci8o6IXOk7jkgQkS4ioiKS5juWWBbt81Qf9+JD/YwFv1+PIznmQd7/BBFZGqn3r3Kso0RknogUicjP67hPRH9/X0TkYRH5QwTeN66+++IiURORb0RkdJVlPxCRaaF5Ve2nqh8e5H0S/cb+v8ANqtpUVT+vujJRP8yHqrrrKVGp6pmq+rTvOEz0+L6+63Ivjjeq+omqHhWaj/A5/jXwX1XNVNUHq64UkQ9F5McROnZMUdVrVfVPvuPwLS4StXgRAwlgZ2CR5xi8ioG/gTFJyT579Sbp7+PmQAmTqIX/hyMiQ0QkX0R2isgGEbk/2Ozj4HV7UD04XERSROT3IrJKRDaKyDMikhX2vlcE67aIyB+qHOePIvKKiDwnIjuBHwTHni4i20VknYj8XUQahL2fish1IvJ1ULT9JxHpLiKfBfG+HL59ld+x2lhFpKGIFAOpwBcisryafUO/+xfB735R2LpfBu+3TkR+GLa8oYj8r4isDs7jwyLSqIbYUkXkPhHZLCIrReSG8NLLIM4ngmOsFZG7RCQ1WPcDEZkWHGtbsP+ZYe99sH0/FZEJIrIF+GNwPj8I/mabReR5EWkebP8s0An4T3Aefh0sHxb8DbaLyBcSVnUjIl1F5KPg7zUFaF3dOQjb/tdBrIUi8mMJK8kUkbNF5PPgb71GRP4Ytt+3iuOlDte1iGQE1+CWIP7ZItI2WLfvv+/azkvYsW4WkfkiskNEJopIRg2/Y4/gnOwI3mti2Lq/Bb/bThGZIyInhK37o4j8O4i3SEQWiEgvEfltcA2uEZExYdt/KCL/IyKzgvd7Q0Ra1hBTbddJjfEmkpqu78D3xX2WN4vIrWH7pIjILSKyPLg2Xg4/xyJyrrjqzO3B36NP2LpvROQ3IjIf2CUiaVWu2VQR+V3w3kXB9dAxWFfjdVKH3/NXYZ+xH1VZV+N9K/QZk5rveWeJyOIg1rUicnP4fjWdYxF5S0R+ViWO+SLy3Rrir/acisgHwMnA34P37lVlv7uBE8LW/z1s9Whx3yvbReQhEZGw/X4kIkvE3V/fE5HOtZzb2r7zavx+k2pqrOTA+0+1n0FxJgR/j53i7gn9g3VPichdwXQLEXlTRDYFv8ebItKhyrH+JO77oEhEJotIrffqsH1D139R8Pf/bti6VSIyKJj+fvA79gvmrxKR/6vLMY6Iqsb8D/ANMLrKsh8A06rbBpgOXB5MNwWGBdNdAAXSwvb7EbAM6BZs+xrwbLCuL1AMjAQa4KoWy8KO88dgfhwu6W0EDAKGAWnB8ZYAvwg7ngJvAM2AfkAJ8H5w/CxgMXBlDeehxljD3rtHLefxgPXASUA5cCeQDpwF7AZaBOsnAJOAlkAm8B/gf2p472uD2DsALYCp4ecaeB14BGgCtAFmAT8J+1uWAVfjks2fAoWA1HHfcuBnwTlvBPQATgMaAtm4BP2Bmq4nIBfYEvz+KcG+W4DssOvp/uD9RgFFwHM1nIczgPXB37Yx8Fz4eQ/O+dHBcY4BNgDjwtYV1HTtU/N1/ZPgb9M4OH+DgGbBug+BHwfTdTkvs4Cc4G++BLi2ht/zReDW4PfIAEaGrbsMaBX8PX4ZnI+MsM/MXuD0YP0zwMrgvdKDa2Bl2Ht9CKwF+gd//1dD554qn2dqv05qjDfRfqq5vkPn6THc5+NY3H2nT7D+RmAG7rPbMDiHLwbregG7gusmHVcttwxoEHaseUBHoFE11+yvgAXAUYAEx25Vx+ukts/YhrBr4gUO/IzVeN/i4Pe8dcAJwXQLYGB1n81qzvGFwMyw+WNx95AG1cR/sHP6IcFntobf/1vrg9//TaA5LoncBJwRrBsbvH+f4Fz/Hvishvc+2Hdejd9vVP/9ui9WavgM4u4Fc4LYJYizfbDuKeCuYLoV8D3cfS4T+Dfwf1WOtTw4v42C+Xtq+D2r/j0vwN33UoCLgr9PKIZngF8G048Gx/hp2LrxEf9M+76p1ClI96EoBraH/eym5kTtY+AOoHWV96nuQnofuC5s/qjgwkwDbiO4YQXrGgOlHJiofXyQ2H8BvF7lAzUibH4O8Juw+fsI+/Ks8l41xhr23oeaqO2pcj424j6IElys3cPWDSfsS7TKe39A8KUYzI8OnWugLe6LoVHY+ktw7TDAJVvLqpxnBdrVcd/VB/kbjAM+r+5aCeZ/Q1jCGyx7D7gSd9MrB5qErXuBmr9EniQsmcUlRzX+XYAHgAlhf4/aErWarusfAZ8Bx1Tz/h9Sw02/hvNyWdj8X4CHa9j3GdxNq0Nt5z7YdhtwbNhnZkrYuu/gPtupwXxmcL6ah8V/T9j2fXGfwVTCPs91uE7qHG+8/1RzfYfOU4ewZbOAi4PpJcCpYevas/8e+Afg5bB1KbjE+aSwY/2olmt2KTC2jnFXvU5q+4yFXxO9Qp8xDnLfopZ7XjC9GvePT7MqxzyJ2hO1jCD+nsH8/wL/qCH+g53TDzm8RC38n6WXgVuC6XeAq6ocbzfQuZr3rvU7r5rtf0Hw/cbBE7VqP4PAKcBXuO+dlCrrniJI1Ko59nHAtirH+n3Y/HXAuzXse8Dfs5r180LXLXAVMCnss/Jj4KVgfhVBMh/Jn3iq+hynqs1DP7g/Qk2uwn14vxRXDXROLdvm4E52yCr23/hzgDWhFaq6G/dfUrg14TPiqnHeFJH14qpD/8y3q8o2hE3vqWa+6WHEeri2qGp52Pzu4PjZuA/pnKCYezvwbrC8ptjCz0X4dGfcf47rwt7rEVypR8j60ERwngniqMu+Vf8GbUXkpaDqYieuVKu2IvDOwAWh9w+OMRL3hZWDuxnsCtt+VTXvEVLbeUBEhorIf4Pi+x24ksg6Fc9T83X9LC6xfElcVdBfRCS96s51PC/rw6ZD10J1fo37UpwVVOHsq34SV326JKji2I4rKQ4/TtXrfbOqVoTNU+W44edwFe56qBr3wa6TGuNNIjX9bTsDr4edtyVABfvvgfuud1WtxP09csPe64BrvIqOuBKIb6nDdVKTqp+x8M9jXe5bNd3zwJXYnAWsCqrphtchHlR1LzARuExEUnD/JDxbS/wHO6eHo7a/79/CzsdW3GehuuPV+p1Xx++3mlT7GVTVD4C/Aw8BG0XkURFpVnVnEWksIo8EVZE7cf+4NpegecNBzkGtxFX3zgs7R/3Dfq+PgBNEpD3uH8SXgREi0gV3zc6r269/+OIpUaszVf1aVS/B3aTvBV4RkSa4bL+qQtyFHBIqQdmAKwYPrwNvhCt+PeBwVeb/CXyJ+8+qGfA73MVZH2qLtb5txn1p9gtLkLNUtaYL/4BzhbtBh6zBlXa0DnuvZqrarw5x1GXfqn+DPwfLjg7+Bpdx4N+g6vZrcCVqzcN+mqjqPcHv1SK4fkI61RJvbecBXGncJKCjqmYBD4fFtgv3JQO49j2EfcHUdF2rapmq3qGqfYHjgXOAK6qJ7WDnpc5Udb2qXq2qObgSiH+Ia4NyAu6GfCGuOqk5sONwjxMIP4edcKU9m6tsU+t1UlO8RxBTLKvuPlebNcCZVa7/DFVdS5V7jogI7u+xto7HWwN0r7rwCK+TdXz7mgg51PvWAVR1tqqOxX3G/g/3pVztptUsexr4PnAqsFtVp9ewb13Oaa1h1nG7kDW42o7wv28jVf2smm0P9p1X2/db6J/ZxmHbt9sXdC2fQVV9UFUH4UrMe+GqzKv6Ja4WaWhw7FGhMA9+Cmomrr3eY8ANuGr55sDC0Puq6jJc0vczXA3aTlxCeA2uVq/ySI5fFwmZqInIZSKSHZzA7cHiSly9fSWujVfIi8B4cQ3Gm+K+zCYG/3G9AnxHRI4X12Dyjxz8osgEdgLFItIb196qvtQWa11s4MDfvUbBuXsMmCAibQBEJFdETq9hl5eBG4NtmuOqE0PvtQ6YDNwnIs3ENV7uLiIn1iGOw9k3E1edtkNEcvn2h77qeXgO93c+XVzj5wxxjYc7qOoqIB+4Q0QaiMhIXHVdTV4GfigifUSkMa6ao2psW1V1r4gMAS4NW/cVkCHugYN0XFuShqGVNV3XInKyiBwdJHY7cYlMdTePg52XOhORC2R/Q95tuC+PyuAY5bjPWpqI3IZrj3kkLhORvsH5vBN4JawEDjj4dVJLvImozp/zwMPA3cEXFiKSLSJjg3UvA2eLyKnBNflLXEJc3Zd8dR4H/iQiPcU5RkRacWTXycu4B7dC18TtoRWHcd/aJ/h8f19EslS1DPdZquka+dY5DhKzSlzzlZpK00LxH8k5PZy/729lfwP4LBG5oIZtD/adV+P3m6puwiWblwX30R8RlqTX9BkUkcHiahrSccneXmq+f+3BPQzYkrC/+xEKFeJsCuL8Ia5ELdxHuETuo2D+wyrzEZWQiRqusekicU9C/g3XFmNPUIx7N/BpUMQ5DNfe4VlcMepK3EXyMwBVXRRMv4T7T6MY156hpJZj34z78i3C3TDq8+myGmOtoz8CTwe/+4V12P43uEaoM8QVNU/F/UdTncdwX5Tzgc+Bt3E34tAX6hW4xqmLcR/SV3BVi3VxqPveAQzE/Yf+Fu6hi3D/A/w+OA83q+oaXIPb3+E+rGtwSUzo83EpMBRXZXA7rq1FtVT1HeBB4L8E5y5YFbpmrgPuFJEiXHuQl8P23RGsfxx3w9sFhD8FWu11jfuv9RXcDXQJ7uZR3RfFwc7LoRgMzAximQTcqKorcFWw7+KSzlW4a7S2qrG6eBbXVmU9ri1QTZ2A1nad1BRvIjrg+q7D9n/DnZPJwXU5A3e9o6pLcSWv/w9XWvUd4DuqWlrHWO7HXeOTcdfnE7iG3od9nQSfsQdw7WKXBa/hDuW+VdXlwDfBftfiSsiqU9M5fgb3sFCNnfXWwzn9G3C+uCcfv9XPWjXHex1XAv9S8HstBM6sYduDfecd7Pvtaty9cwvugarw5LOmz2Cz4L224a6FLcBfqwnvAdy1sxl3jb57sN+9LlR1MS65no5Lgo8GPq2y2Ue4RPHjGuYjKvRUnamDoBRrO67Yd6XncGKauO41HlbVzgfdOIGJe+x+IdDwEEo+TUBEPsQ1Kn/cdyzGHIyIXAFco6ojfcdSH+w7LzYkaolavRGR74hrxNgE9yTPAtwTPyaMiDQS1wdRWlCtdjuuu4SkI24or4Yi0gL3n+x/LEkzJrEF1bDX4Z5sjFv2nRd7LFE7uLG4xp+FQE9cdZMVQ36b4KrWtuGqPpfgqvaS0U9w1QXLcVW/9dlO0RgTY4I2cJtwVWcveA7nSNl3Xoyxqk9jjDHGmBhlJWrGGGOMMTEqYomaiDwpbuyuhWHL/ioiX4obA+11OXCcwd+KyDIRWRr+KLWInBEsWyYit0QqXmOMMcaYWBOxqk8RGYV7tPcZVQ0NsDoG+EBVy0XkXgBV/Y2I9MX1ETYE1zPyVFynd+Ae3z4N103BbOCS4HHaGrVu3Vq7dOlS/7+UMSZmzZkzZ7Oq1jRyRlyxe5gxyaW2+1dadQvrg6p+LG6IhfBlk8NmZwDnB9NjcWNnlQArRWQZLmkDNwbkCgAReSnYttZErUuXLuTn5x/5L2GMiRsiUtvQXnHF7mHGJJfa7l8+26j9CDdYLLgxx8I7OywIltW03BhjjDEm4XlJ1ETkVlyv9c/X43teIyL5IpK/adOm+npbY4wxxhhvop6oicgPcINGfz+sb5a1HDjIbodgWU3Lv0VVH1XVPFXNy85OiGYqxhhjjElyUU3UROQM4NfAucG4myGTgIuD3ty74jrZm4V7eKCnuEHIGwAXB9saY4wxxiS8iD1MICIvAicBrUWkADek0G+BhsAUEQGYoarXquoiEXkZ95BAOXC9qlYE73MDbgDfVODJYNBYY4wxxpiEF8mnPi+pZvETtWx/N3B3NcvfBt6ux9CMMcYYY+KCjUxgjDEHUdeOt0XkeyKiIpIXtqzazryNMaYuIlaiZowxiUBEUoGHCOt4W0QmVe14W0QygRuBmWHL+uLa1vYj6MxbRHqFmnYYY8zBWImaMcbUbghBx9uqWgqEOt6u6k/AvcDesGX7OvNW1ZVAeGfexhhzUFaiZowxtauu4+2h4RuIyECgo6q+JSK/qrLvjCr7Wqfd5kCVlTB1KkyeDFu2QMeOMG4cDBzoOzITAyxRM8aYIyAiKcD9wA+O8H2uAa4B6NSp05EHZuLDokXwgx9A1SHD/vQnGDsWHn0U2rTxEpqJDZaombjWp38fCgsLa90mJyeHJQuXRCkik4AO1vF2JtAf+DDodqgdMElEzq3Dvvuo6qPAowB5eXla3TYmwUyd6krOdu2CnBz44Q+ha1f4/HN46il44w2YOxc++AB69PAdrfHEEjUT1woLCxn/2vhat5lw3oQoRWMS1L6Ot3FJ1sXApaGVqroDaB2aF5EPgZtVNV9E9gAviMj9uIcJQp15m2Q3Ywaccw6UlMCll8I//wnNmu1ff8stcMEFbruTToKZMyHXas2TkT1MYIwxtVDVciDU8fYS4OWgk+47g1Kz2vZdBIQ6836XsM68TRIrLITzznNJ2tVXw7PPHpikAXTo4NqsjRgBa9fC+ee77U3SsRI1Y4w5iOo63lbV22rY9qQq89V25m2SlCpcdRWsW+dKyh56CFJqKDPJzITXX4dBg1zJ2l13ubZrJqlYiZoxxhgTLc8+C+++Cy1awIsvQnp67dtnZ7vtROCee2D+/OjEaWKGJWom4agqlVrpOwxjjDnQzp3wy1+66fvvh3bt6rbfiBFw3XVQXg7XX+9K5UzSsKpPkzBUlekF05leMJ09ZXvo3LwzY4+qrl9SY4zx4H//FzZvhuOPhyuvPLR9//xnmDgRpk2D//wHzq21eaRJIFaiZhKCqvLW128xZcUU2jRpw5DcIRTsLODROY9SmWmla8YYzzZscKVoAPfe66oyD0WzZnBb0Czyllugwp5JSRaWqJmEMLtwNnPWzeH4jsdz2dGXMab7GK4acBVllWXsOWMPZRVlvkM0xiSz++5z/aV95zswcuThvcdPfgJdusCSJe4hA5MULFEzcW/Trk1MXj6Zni17MrrraIJOR2nTpA3n9DyHipwKHpv7mOcojTFJa+dOeOQRN31btQ8L102DBnDzzW763nutrVqSsETNxDVFeWfZO6SnpjP2qLH7krSQ/m36k1qQyh0f3UFxabGnKI0xSe2JJ1yyNmoU5OUd2Xv98IfQurUbcuq//62f+ExMs0TNxLXy7uWs3L6SU7qcQpMGTb61XkRoOK0hG3dt5KFZD3mI0BiT1MrL4YEH3HSoNOxING4MP/+5mw69r0lolqiZuFVRWUHJ8BJaN27NoJxBNW6XtiGN0d1G8+CsBymtKI1ihMaYpPfmm7B6NfTqBWefXT/vec01kJYGb70FBQX1854mZlmiZuLWiwtfpLJVJSd1OYkUqf1SvmnYTRQWFfLyopejFJ0xxgCPP+5er7225hEIDlXbtvDd70JlpatWNQnNEjUTl8oqyrj9w9tJ2ZhC39Z9D7r9GT3O4KhWR/GP2f+IQnTGGIMbo/Odd9zoA5dfXr/v/ZOfuNfHH3fVqyZhWaJm4tJT855ixbYVNJze8FsPEFRHRPjJoJ8wvWA6CzYsiEKExpik99RTrtRr3Dj3AEB9Ovlk6NHDVX2+9179vreJKZaombhTWlHKXZ/cxbAOw0j7pu6Da1x53JU0TG3II3MeiWB0xhiD6zrjySfd9FVX1f/7p6TAj37kpp9/vv7f38QMS9RM3Hlq3lOs3rGaP574R4S69+7dslFLLuh3Ac/Of5ZdpbsiGKExJunNmgUrVkBODoweHZljXHyxe33jDdeZrklIlqiZuFJaUcqfP/kzQ3OHMqb7mEPe/yeDfsLOkp1MXDQxAtEZY0zgpZfc60UXQWpqZI7RtasbN3T3bpg0KTLHMN5ZombiytPznmbVjlXcfuLtdWqbBlBcXExWyyyyWmZx1nFnkbIlhasfvXrfsqyWWfTp3yfCkRtjkkZFhRtAHVyiFkmXXupeX3ghsscx3tS9gY8xnpVVlPHnaX9mcM5gzuhxRp33q6ysZPxr4/fNzyiYwXvL3+OSZy6hXdN2AEw4b0K9x2uMSVLTpsG6dW5cziFDInusCy6AG2+Ed9+FLVugVavIHs9EnZWombjxzBfP8M32bw6pNK06x7Y9lrSUNOasm1OP0RljTCBUmnbxxXAE96o6adPGtYErL7eB2hOUJWomLpRVlHH3J3eTl5PHWT3POqL3apTeiH7Z/Zi/Yb6NVGCMqV8VFfDqq2460tWeIeed517feCM6xzNRZYmaiQvPzX+OldtXHnFpWsig9oMorShl4caF9RCdMcYEZs6EjRuhWzc49tjoHPPcc13J3ZQpUFwcnWOaqLFEzcS8sooy7vrkLga1H8TZPetnrLwOzTrQpkkbq/40xtSv0NOXoeQpGtq1g2HDoKQEJk+OzjFN1FiiZmLey4teZsW2Fdx24m31UpoGbqSCQe0HUVhUyLqidfXynsYYc0CiFk3jxrnX//u/6B7XRJwlaiamqSoPzHyAo1odxTm9zqnX9z6m7TGkpaSRvy6/Xt/XJB4ROUNElorIMhG5pZr114rIAhGZJyLTRKRvsLyLiOwJls8TkYejH72Jmq+/hiVLoHlzGDkyusceO9a9vvkmlJVF99gmoixRMzFtesF08gvz+fnQn5Mi9Xu5ZqRl0L9NfxZuXIima72+t0kcIpIKPAScCfQFLgklYmFeUNWjVfU44C/A/WHrlqvqccHPtVEJ2vgRKk07+2w3EHs0HXUU9O4N27bBxx9H99gmoixRMzGpT/8+ZLXM4sTfnAglcMs5txzQQW3op7joyBrOhh4qKOtj/4GaGg0BlqnqClUtBV4CxoZvoKo7w2abAJb5JyNf1Z4hoVK1t9/2c3wTEdbhrYlJhYWFXPXiVTww4wGGdxjOmInVDxd1x+g7jug4HZp1oGOzjhQMLKC8spy0FPtImG/JBdaEzRcAQ6tuJCLXAzcBDYBTwlZ1FZHPgZ3A71X1kwjGanzZts11dJueDqef7ieGM8+Ee+91nd/ed5+fGEy9sxI1E7NmF84GYEhuZHv2Pr7j8WiW8sriVyJ6HJPYVPUhVe0O/Ab4fbB4HdBJVQfgkrgXRKRZdfuLyDUiki8i+Zs2bYpO0Kb+vP8+VFbCiBGQleUnhuHDITMTFi+G1av9xGDqXcQSNRF5UkQ2isjCsGUtRWSKiHwdvLYIlouIPBg01J0vIgPD9rky2P5rEbkyUvGa2KKifLH+C3q27EnzjOYRPdZRrY4iZWsK9356L6pWY2W+ZS3QMWy+Q7CsJi8B4wBUtURVtwTTc4DlQK/qdlLVR1U1T1XzsrOz6yNuE02hbjHGVF/6HxUNGsCpp7rp997zF4epV5EsUXsKqDog4y3A+6raE3g/mAfXSLdn8HMN8E9wiR1wO66aYQhweyi5M4mtomMFRaVFHNsu8h1GiggN5jRg3vp5TF0xNeLHM3FnNtBTRLqKSAPgYmBS+AYi0jNs9mzg62B5dvAwAiLSDXePWxGVqE30qO5PjHxVe4acEXztvvuu3zhMvYlYgxxV/VhEulRZPBY4KZh+GvgQV00wFnhGXXHGDBFpLiLtg22nqOpWABGZgkv+XoxU3CY2lPUpIyMtg16tqi18qHclc0qQ4cKZfz6TJq83qXabnJwclixcEpV4TOxQ1XIRuQF4D0gFnlTVRSJyJ5CvqpOAG0RkNFAGbANCpf+jgDtFpAyoBK4N3c9MAvnqK1fVmJ0Nxx3nN5ZQojh1quumI9pPn5p6F+2W021VNdS76HqgbTBdXWPd3FqWmwS2t3wvZd3KODr76Kg17tcyZfQxo5m6YioXPX0ROZk539pmwnkTohKLiT2q+jbwdpVlt4VN31jDfq8Cr0Y2OuNdqDTttNMgxXPT7y5dXFcdS5fCjBlwwgl+4zFHzNsVFZSe1VuDIGuImzimrpgKDaBP6z5RPe6g9oNomNqQz9Z8FtXjGmPiXCy0Twtn1Z8JJdqJ2oagSpPgdWOwvKbGunVuxGsNcRPHa0tegxLo2rxrVI+bkZbBoJxBLN60mK17rHbKGFMHJSXw3/+66VhL1Gzcz4QQ7URtEvvbblwJvBG2/Irg6c9hwI6givQ9YIyItAgeIhgTLDMJqryynElLJ5G+Mp3UlNSoH39Y7jBSJIXpBdOjfmxjTBz67DPYvRuOPhrat/cdjXPCCa5t2ty5sH2772jMEYpk9xwvAtOBo0SkQESuAu4BThORr4HRwTy4th8rgGXAY8B1AEGj2z/hnrqaDdxpDXET2yerPmHLni2kLfPT8Wxmw0yOaXsM89bPY1fpLi8xGGPiSKxVewI0aQJDh7p+3Ww4qbgXsURNVS9R1faqmq6qHVT1CVXdoqqnqmpPVR0dSrrUuV5Vuwfj5eWHvc+Tqtoj+PlXpOI1seH1L18nIy2DtFX+Rgg4vuPxlFeWM3PtTG8xGGPiRKjac/Rov3FUdUowOMYHH/iNwxwxG5nAxAxV5fUvX+eMHmcg5eItjtaNW9O7VW9mF86mtKLUWxzGmBhXVAT5+ZCa6kYkiCWWqCUMS9RMzJi3fh4FOwsYd9Q436EwtMNQ9pbvZclm6zfNGFODTz+FigrIy3NDN8WSYcMgIwMWLADrCSGuWaJmYsZ7y91zIqf38NyzN9A5qzPNM5ozf8N836EYY2JVqNrz5JP9xlGdhg33l/J9+KHXUMyRsUTNxIzJyydzbNtjade0ne9QEBGOaXMMK7etpKikyHc4xphYFEqATjrJZxQ1s+rPhGCJmokJxaXFTFs9jTHdY+fJqaPbHo2iLN682HcoxphYs3MnzJkTm+3TQixRSwiWqJmY8NE3H1FWWcbp3f1Xe4a0btya7MbZfLnpS9+hGGNiTah92uDB0LSp72iqF2o799VXsLbavuJNHLBEzcSEycsn0yitESM6xdZ/pr1b92bVjlXWp5ox5kCh9mmxWu0JkJa2f6xPa6cWtyxRMzHhveXvcVKXk8hIy/AdygH6tO6Dony15SvfoRhjYkko8YnFBwnCjRrlXj/5xG8c5rBZoma8W7V9FUu3LI2p9mkh7Zq2I7NBJsu2LvMdijEmVoTap6WlwfHH+46mdiNHutdp0/zGYQ6bJWrGu8nL3RAssZioiQg9WvZg+bblqKjvcIwxsWDaNDc8Uyy3TwvJy3NddSxaBFu2+I7GHAZL1Ix3k1dMpkOzDvRp3cd3KNXq0bIHJRUlVLSr8B2KMSYWfPSRez3xRL9x1EXDhm7cT3APQJi4Y4ma8aq8spypK6YyptsYRPwNG1Wbbi26IQjlXcp9h2KMiQWhhCfUUD/WheK0dmpxyRI141V+YT7b926PidEIapKRlkFuZi4VHaxEzZikt3cvzJ7tpocP9xtLXVmiFtcsUTNeTV4+GUE4teupvkOpVefmnaloW8Hust2+QzHG+DRnDpSWQr9+0KKF72jqZvhwSElxse+yrobijSVqxqv3lr/H4NzBtGrcyncoteqc1RlSYfqa6b5DMcb4FKr2DD1NGQ+aNYNjj4Xycpg503c05hBZoma82b53OzMLZjKmW+w97VlVp6xOUAkffvOh71CMMT6FErVYHTaqJqHqT+umI+5Yomairk//PmS1zCLnhBwqtIL7f34/WS2zDvgpLir2HeYBGqY1JGVTCh+t+sh3KMYYX1Tjs0QNrJ1aHEvzHYBJPoWFhYx/bTxvfvUmCzYu4KYHbyI1JfWAbe4YfYen6GqWVpDGzJyZ7CnbQ6P0Rr7DMcZE29Klri+y9u2hSxff0RyaUKI2fbqrAk2zr/94YSVqxgtVZfm25XRt3vVbSVqsSl2bSmlFKTPXWhuPZCMiZ4jIUhFZJiK3VLP+WhFZICLzRGSaiPQNW/fbYL+lIhK7jzebgwuv9ozR7oRq1LYt9OzpHib4/HPf0ZhDYIma8WLrnq1s37ud7i27+w6lztLWpiEIH31j1Z/JRERSgYeAM4G+wCXhiVjgBVU9WlWPA/4C3B/s2xe4GOgHnAH8I3g/E49C7bvirdozxNqpxSVL1IwXy7ctB6B7i/hJ1KRUOK7dcdZOLfkMAZap6gpVLQVeAsaGb6CqO8NmmwCh8cbGAi+paomqrgSWBe9n4lG8PkgQEur3bbo9vR5PLFEzXizftpwWGS1o2ail71AOyajOo5hRMIOyijLfoZjoyQXWhM0XBMsOICLXi8hyXInazw9l32D/a0QkX0TyN23aVC+Bm3q0cSN8/TU0buy6uohHoQHkLVGLK5aomajTFOWb7d/QrUU336EcsuEdhrOnfA/zN8z3HYqJMar6kKp2B34D/P4w9n9UVfNUNS87O7v+AzRHJlSaNmwYpKf7jeVw9e4NzZtDQYH7MXHBEjUTdRVtKyitKI2ras+QYR2GATCjYIbnSEwUrQU6hs13CJbV5CVg3GHua2JVvFd7ghudIDRAu5WqxQ1L1EzUVXRyY2Z2ad7FbyCHoVNWJ9o3bc+MtZaoJZHZQE8R6SoiDXAPB0wK30BEeobNng18HUxPAi4WkYYi0hXoCcyKQsymvsVr/2lVWTu1uGMdqZioK+9YTk5mTlz2RSYiDOswzErUkoiqlovIDcB7QCrwpKouEpE7gXxVnQTcICKjgTJgG3BlsO8iEXkZWAyUA9eraoWXX8Qcvr173TiZIq7qM56FErXPPvMbh6kzS9RMVBWXFlPRroJuzeOvfVrIsA7DeP3L19m8ezOtG7f2HY6JAlV9G3i7yrLbwqZvrGXfu4G7IxedibjPP4eyMujf342bGc+GDnUJ59y5LgHNyPAdkTkIq/o0UfXxqo8hFbq26Oo7lMMWaqc2s8A6vjUmKcwIStDjvTQNICsL+vZ1iefcub6jMXVgiZqJqqkrpkJ5MMh5nBrUfhCpkmrVn8Yki0RK1MDaqcUZS9RMVE1dMZXUwlTSUuK31r1JgyYc0/YYe6DAmGSRaIma9acWVyxRM1GzoXgDCzYuIG11/CZpIcM6DGNmwUwqKq1duDEJbd06WL0aMjNdP2SJILxETbX2bY13lqiZqPlg5QcApK1JjEStqLSILzd/6TsUY0wkzQzaog4ZAqkJMkxrr17QogUUFsKaNQff3nhliZqJmqkrptIiowUpm+L/srOOb41JEolW7Qmu49vQ72PVnzEv/r8xTVxQVaasmMIpXU9BVHyHc8R6tuxJi4wWlqgZk+gSMVED608tjliiZqJi2dZlrNm5htHdRvsOpV7s6/jWHigwJnGVl8Ps2W46NPRSorAnP+NG/DcWMnHh/ZXvA3Bq11M9R3L4iouLyWqZtW++ZEgJJcNKaNauGVK6v5QwJyeHJQuX+AjRGFOfFi2C3buhWzfIzvYdTf0aMsR1fPv557BnDzSKv5FikoWXRE1ExgM/BhRYAPwQaI8bzLgVMAe4XFVLRaQh8AwwCNgCXKSq3/iI2xy+j1d9TE5mDj1a9vAdymGrrKxk/Gvj980v37qc5xY8x7h/jqNbi/0jLUw4b4KP8Iwx9S1Rqz3BjbBw9NEwf74bHivexzBNYFGv+hSRXODnQJ6q9seNnXcxcC8wQVV74MbKuyrY5SpgW7B8QrCdiSOqyserPmZU51GIxH/7tJDcZrkAFOws8ByJMSYiEjlRg/2/10wbZSWW+WqjlgY0EpE0oDGwDjgFeCVY/zQwLpgeG8wTrD9VEunbPgms2LaCtUVrObHzib5DqVcZaRm0bNSSdUXrfIdijImEUKKWaO3TQkK/lyVqMS3qiZqqrgX+F1iNS9B24Ko6t6tqebBZAZAbTOcCa4J9y4PtW0UzZnNkPl71MQCjOo/yHEn9y83MZW3RWt9hGGPq2/bt8OWX0LAhHHec72giwxK1uOCj6rMFrpSsK5ADNAHOqIf3vUZE8kUkf9OmTUf6dqYefbz6Y1o3bk2f1n18h1LvcjJzKCotoqikyHcoxpj6NGuWex04EBo08BtLpPTu7UZcWL0a1q/3HY2pgY+qz9HASlXdpKplwGvACKB5UBUK0AEIFVOsBToCBOuzcA8VHEBVH1XVPFXNy060p3Pi3EfffMQJnU5IqPZpITmZOQAUFhV6jsQYU68SvX0auJEWBg9201aqFrN8JGqrgWEi0jhoa3YqsBj4L3B+sM2VwBvB9KRgnmD9B6o2OFms6tO/D1kts/b9ZHbKZOX2lbz9yNv7lhUXFfsOs960a9oOQSxRMybRJHr7tBCr/ox5Ue+eQ1VnisgrwFygHPgceBR4C3hJRO4Klj0R7PIE8KyILAO24p4QNTGqsLDwgC4sFmxYwGtfvsaV46+kfWZ7AO4YfYev8Opdg9QGZDfJtkTNmESiuj9xSeQSNbBELQ546UdNVW8Hbq+yeAUwpJpt9wIXRCMuU/9W7VhFw9SGtG3a1ncoEZOTmcPSzUtR1YSs3jUm6SxbBlu3Qrt20KmT72giK5SozZ4NFRWJM/B8ArEhpExErd6xmo7NOpIiiXup5Wbmsqd8D9v3bvcdijGmPoS3T0v0f75CyWhRkXvK1cScxP32NN7tLd/Lpt2b6NCsg+9QIsoeKDAmwSRL+7SQIUFlllV/xiRL1EzErN3pHtxN9EStbZO2pEqqJWrGJIpkeOIznLVTi2mWqJmICQ2tlOiJWmpKKm2btrVEzZhEsHu3G/8yJQXy8nxHEx2WqMU0S9RMxBTsLKBNkzY0TGvoO5SIy8nMobC4EOs5JjGJyBkislRElonILdWsv0lEFovIfBF5X0Q6h62rEJF5wc+k6EZuDtncuVBeDv37Q9OmvqOJjkGD3EMECxbArl2+ozFVWKJmIkJVKSgqSPjStJCczBxKK0rZvHuz71BMPRORVOAh4EygL3CJiPStstnnQJ6qHoMbk/gvYev2qOpxwc+5UQnaHL5kq/YEaNwYjj4aKithzhzf0ZgqLFEzEbF592b2lu+lQ2ZyJGq5mW5oWqv+TEhDgGWqukJVS4GXcMPg7aOq/1XV3cHsDNzoKiYeJWOiBlb9GcMsUTMRkSzt00JaN25Nekq6DdCemHKBNWHzBcGymlwFvBM2nxGMQzxDRMZFID5Tn0KJSrI88RliiVrM8tLhrUl8a3auISMtg9aNW/sOJSpSJIX2me1ZV7TOdyjGIxG5DMgDTgxb3FlV14pIN+ADEVmgqsur2fca4BqAToneyWqsWrsWCgqgWTM3YHkysUQtZlmJmomItUVr6ZDZIal66s/JzGH9rvVoij1QkGDWAh3D5jsEyw4gIqOBW4FzVbUktFxV1wavK4APgQHVHURVH1XVPFXNy87Orr/oTd2FkpTBg91Tn8mkd2+XoBYUQKE14YglSXYlmmjYW76Xjbs2Jk21Z0hOZg7lleVUtqr0HYqpX7OBniLSVUQa4MYbPuDpTREZADyCS9I2hi1vISINg+nWwAhgcdQiN4cmWas9wSWmgwe7aStViymWqJl6lywd3VYVeqCgom2F50hMfVLVcuAG4D1gCfCyqi4SkTtFJPQU51+BpsC/q3TD0QfIF5EvgP8C96iqJWqxKpkTNbDqzxhlbdRMvQs9SJDbrLb21omnRUYLMtIyLFFLQKr6NvB2lWW3hU2PrmG/z4CjIxudqRcVFZCf76YtUfMbhzmAlaiZelews4DsxtlkpGX4DiWqRISczBwq2liiZkzcWbTIdfbapQu0bes7Gj9CiVp+vktcTUywRM3UKyW5OrqtKiczh8rWlewp2+M7FGPMoUj2ak9wCWrnzlBcDIuthj5WHDRRE5E5InK9iLSIRkAmvlW2qGRv+V46Nut48I0TUE5mDqTAFxu+8B2KMeZQWKLmhH7/WbP8xmH2qUuJ2kVADjBbRF4SkdMlmfpcMIck1D4r2dqnhYQeKJi9drbnSIwxh8QSNcfaqcWcgyZqqrpMVW8FegEvAE8Cq0TkDhFpGekATXypbFtJekp60nR0W1Vmg0xkl5C/Lt93KMaYuioqcm3U0tJgQLXd3CUPS9RiTp3aqInIMcB9uEfQXwUuAHYCH0QuNBOPKtpW0D6zPSmSnM0fRYTUjankF1qiZkzcyM8HVTj2WGjUyHc0fg0c6BLWhQtdWzXjXZ3aqAETcJ0+HqOqP1fVmap6H7Ai0gGa+FFWUUZFdgU5TXN8h+JV6oZUlmxaQnGp3eSMiQtW7blfo0ZwzDFQWQlz5viOxlC3ErULVPVUVX0hfFgUAFU9L0JxmTi0aNMiSIOcZsmdqKVsSEFR5q6b6zsUY0xdWKJ2IKv+jCl1SdR+LCLNQzPBkCh3RS4kE69CDehDDeqTVerGVACr/jQmHqhaolaVJWoxpS6J2pmquj00o6rbgLMiFpGJW/mF+bDX9dCfzFJ2p9CxWUdL1IyJBwUFsG4dNG8OPXv6jiY2WKIWU+qSqKWGBhUGEJFGQMNatjdJanbhbFI3pmK9t0BeTh6zC62LDmNiXigZGTLEDUxuoFcvyMqCtWvdj/GqLlfl88D7InKViFwFTAGejmxYJt7sLd/Lgo0LSN2Q6juUmJCXk8eyrcvYtmeb71CMMbWxas9vS0lxiStYqVoMqEs/avcCdwN9gp8/qepfIh2YiS9frP+C8spyS9QCg3MGA9gDBcbEOkvUqmfVnzEjrS4bqeo7wDsRjsXEsVA1nyVqzqCcQYBrt3dqt1M9R2OMqVZ5+f4uKEIlSMaxRC1m1KUftfNE5GsR2SEiO0WkSER2RiM4Ez9mF86mbZO2SLG1TwNo2agl3Vp0s3ZqxsSyhQth927o1g2ys31HE1tCiVp+PlRU+I0lydWljdpfgHNVNUtVm6lqpqo2i3RgJr7kF+YzOHcwgiVqIYNzBtuTn8bEMqv2rFl2NnTtCrt2ueG1jDd1SdQ2qOqSiEdi4lZRSRFLNi0hr32e71BiSl5OHqt2rGLTrk2+QzHGVMcStdpZ9WdMqEuili8iE0XkkqAa9DwRsREJzD5z181FUQbnDvYdSkzJy3GJq5WqGROjLFGrnSVqMaEuiVozYDcwBvhO8HNOJIMy8SWUiIQSE+MMbD8QQSxRMyYW7dwJS5ZAejocd5zvaGLTsGHudcYMv3EkuYM+9amqP4xGICZ+zS6cTaesTrRp0sZ3KDGlWcNmHNX6KPLXWaIWK0SkDTACyAH2AAuBfFWt9BqYib7Zs93wUccdBxkZvqOJTQMGQIMGsHgx7NjhOsE1UVeXpz57icj7IrIwmD9GRH4f+dBMvJhdOHtfv2HmQHk5eVaiFgNE5GQReQ94CzgTaA/0BX4PLBCRO0TEHpJKJlbteXANG8LAgQeOh2qiri79qD0G/Ap4BEBV54vIC4ANzG7YumcrK7at4JqB1/gOJWYUFxeT1dL951lybAklJ5WQ2TGTlF37/y/KyclhyUJ7RieKzgKuVtXVVVeISBquOcdpwKvRDsx4Yola3Qwf7qo+p0+HMWN8R5OU6pKoNVbVWVXGbyyPUDwmzlj7tG+rrKxk/GvjAVizYw1PznuScyacw1Gtj9q3zYTzJvgKLymp6q9qWVcO/F9t+4vIGcDfgFTgcVW9p8r6m4Af4+6Nm4AfqeqqYN2VuJI7gLtU1Ybg8y28hMgStdoNHw4TJrhEzXhRl4cJNotId0ABROR8YN2RHFREmovIKyLypYgsEZHhItJSRKYEnetOEZEWwbYiIg+KyDIRmS8iA4/k2KZ+zV7rOnQN9cRvDtSuaTsEobCo0HcoBhCRZ0UkK2y+i4i8f5B9UoGHcFWmfYFLRKRvlc0+B/JU9RjgFVz/k4hIS+B2YCgwBLg9dG8zHq1eDRs2QMuW0KOH72hiW+iBgpkzodKacvpQl0Ttely1Z28RWQv8AvjpER73b8C7qtobOBZYAtwCvK+qPYH3g3lwN8eewc81wD+P8NimHuWvy6dXq140z2juO5SYlJ6aTpsmbSxRix3TgJkicpaIXA1MBh44yD5DgGWqukJVS4GXgLHhG6jqf1V1dzA7A+gQTJ8OTFHVraq6DZgCnFE/v4o5bKHSoaFDQayT7lp17Ai5ubB9Oyxd6juapFSXQdlXqOpoIBvoraojVfWbwz1g8N/sKOCJ4P1LVXU77sYXqhJ4GhgXTI8FnlFnBtBcRNof7vFN/covzLdqz4Non9mewuJCVNV3KElPVR/BVVG+AdwJjFLV/xxkt1xgTdh8QbCsJlexf2zkQ93XREMoUTv+eL9xxIvhw92rVX96UZenPm8TkduAXwLjw+YPV1dcG45/icjnIvK4iDQB2qpqqEp1PdA2mK7TjU5ErhGRfBHJ37TJeoKPhvXF6ynYWWBPfB5ETmYOu8t2s6Nkh+9Qkp6IXA48CVwBPAW8LSLH1uP7XwbkAX89jH3tHhYtoYQjlICY2lmi5lVdqj53hf1U4KoiuxzBMdOAgcA/VXVA8L63hG+grujhkIofVPVRVc1T1bxsG1w3KuxBgrrJzXT/V1j1Z0z4HjBSVV9U1d8C1+ISttqsBTqGzXcIlh1AREYDt+LGRi45lH3B7mFRs2cPfP45pKTAkCG+o4kP1vGtV3Wp+rwv7Odu4CSg2xEcswAoUNVQpyyv4BK3DaEqzeB1Y7C+zjc6E135hfmkSArHtTvOdygxrU2TNqRIiiVqMUBVx6nqxrD5WbiG/rWZDfQUka4i0gC4GJgUvoGIDMC15T03/P2B94AxItIieIhgTLDM+JKfD+XlcPTRkJnpO5r4MHCgG8Fh0SLX8a2JqrqUqFXVmP0NZQ+Zqq4H1ohIqK+CU4HFuBvflcGyK3FtSAiWXxE8/TkM2BFWRWo8yi/Mp0/rPjRt0NR3KDEtLSWNtk3aWqLmkYj8PngC81tUtVREThGRaofGC7rvuAGXYC0BXlbVRSJyp4icG2z2V6Ap8G8RmScik4J9twJ/wiV7s4E7g2XGl88+c69W7Vl3GRn7O76dNct3NEnnoP2oicgC9ldDpuIeKrjzCI/7M+D54L/TFcAPcUnjyyJyFbAKuDDY9m1cZ5XLcGOO2pBWMUBVyS/M54we9gBbXeQ2y2XBhgWoKmJPmfmwAPiPiOwF5uLayWbgniY/DpgK/LmmnVX1bdy9KHzZbWHTo2vZ90lcuzgTC+xBgsMzfLjromP6dDjtNN/RJJW6dHgb/l9mObAh+A/zsKnqPFyD26pOrWZbxXURYjzr078PhYWuVKiyaSXFVxUz8YGJvP6j1/dtU1xU7Cu8mJabmUt+YT6bd28mu4m1P/LgfFUdISK/xjWraA/sBJ4DrlHVPV6jM9GhaiVqh2v4cHjgAXugwIO6JGpFVeabhZcIWDF+8igsLNzX4/6Xm79k4qKJfP+m79Oh2f6a8DtG3+ErvJgWOkcFOwssUfNjkIjkAN8HTq6yrhFugHaT6FasgE2bIDsbunf3HU18qdrxbcrhtJwyh6MuidpcXGP+bYAAzYHQeHnKkT1YYOJUYVEhgtC2SduDb2xo1agVGWkZFBQVMKD9AN/hJKOHcR1pdwPyw5YLdh9LHuHdclgThEPTsSPk5EBhIXz1FfTu7TuipFGXlHgK8B1Vba2qrXBVoZNVtauq2s0tSRUWFdKmSRvSU9N9hxIXRITczFzW7rQHln1Q1QdVtQ/wpKp2C/ux+1gysWrPwydi/al5UpdEbVjQkBYAVX0HsFaYSUxVKSwqJCczx3cocaVDsw5s3LWRkvKSg29sIkJVj3T4OxPP7EGCI2OJmhd1qfosFJHf4xrdgmvjYf0MJLEdJTvYU77HErVD1KFZBxS1bjqM8aGoCObPh7Q0yLNOug+LdXzrRV1K1C7BdcnxOvBaMH1JJIMysS2UaFiidmhCIxQUFBV4jsSYJDR7tmsEf9xx0Lix72ji06BBruPbhQth507f0SSNuoxMsFVVb8QNuzJQVX9hT3omt8KiQlIkhTZN2vgOJa40Sm9Eq0atrJ2aMT7Y+J5HLiMDBgywjm+jrC6Dsh8vIotxPXIjIseKyD8iHpmJWYVFhbRt0pa0lLrUnJtwHZp1oGBnAXpoQ9kaY45U6EECa592ZEKJbuh8moirS9XnBOB0YAuAqn4BjIpkUCZ22YMERya3WS67ynahmZaoGRM1lZX721VZidqRGTHCvX76qd84kkideqxT1TVVFlVEIBYTB7bt3UZJRYklaoepQ6br+LaivX2EjImar76CrVuhfXvo1Ml3NPEtlKh99pkb3N5EXF0StTUicjygIpIuIjcTVIOa5GMPEhyZNk3akJaSRkU7S9SMiZpQ6c/xx1tHt0cqJwe6dYPiYliwwHc0SaEuidq1uLE2c4G1uAGMbezNJFVYVEiqpJLd2IZBOhypKankZOZYomZMNH3yiXs94QS/cSSK0HkMnVcTUbUmaiKSCvxNVb+vqm1VtY2qXqaqW6IUn4kxhUWFtGvajtSUVN+hxK3czFwqsius41tjoiWUUIwc6TeORBE6j9Om+Y0jSdSaqKlqBdBZRBpEKR4TwxRlXfE6q/Y8Qh2adYA0+Hz9575DMSbxFRa6wdibNoVjj/UdTWIIT9TUHoyKtLr0r7AC+FREJgG7QgtV9f6IRWViUmWLSkorSi1RO0Idm3UEYPqa6QzrMMxzNMYkuFCpz/HHu1EJzJE76iho3RrWrXNJcPfuviNKaDWWqInIs8HkucCbwbaZYT8myVS0ce2qLFE7MpkNM5EdwrQ1Vm1gTMRZ+7T6J2LVn1FU278Xg0QkB1gN/L8oxWNiWEW7CtJT0mnduLXvUOJeWmEa01ZPQ1URewrNmMixRC0yRo6E//s/l6hdeaXvaBJabYnaw8D7QFcgP2y5AAp0i2BcJgZVtKugU2YnUqRO3e+ZWqQWprJx10aWbV1Gz1Y9fYdjTGLavt0NxJ6eDkOG+I4msViJWtTU+I2rqg+qah/gX6raLeynq6pakpZk9pTtoTK7ktxmub5DSQiphe6p2Wmr7SZnTMR89plr7J6XB40a+Y4msQwc6M7pl1/Cpk2+o0lodRmU/afRCMTEts/Xfw6pwROL5oilbE2hZaOWlqgZE0lW7Rk56ekwLHgYyoaTiiirwzJ1MqPAjZNniVr9EIQRHUfYAwXGRFKoWs4Stciw6s+osETN1MmMghnIDqFpg6a+Q0kYIzuN5KstX7Fx10bfoZiDEJEzRGSpiCwTkVuqWT9KROaKSLmInF9lXYWIzAt+JkUv6iS3dy/MmuWmQ+NTmvoVStRshIKIskTN1MmMghmkrrfRCOrTyE7uJvfpaqs2iGXBCC0PAWcCfYFLRKRvlc1WAz8AXqjmLfao6nHBz7kRDdbsN3s2lJZC//7QooXvaBLT8OGQkgJz58KuXQff3hwWS9TMQa3duZY1O9dYolbPBrUfRMPUhtZOLfYNAZap6gpVLQVeAsaGb6Cq36jqfKDSR4CmGtY+LfIyM+G446C8HGbM8B1NwrJEzRzUzLUzAUhdZ4lafWqY1pAhuUOsnVrsywXWhM0XBMvqKkNE8kVkhoiMq9fITM0sUYuOk05yr//9r9cwEpklauagZhTMoEFqA1I3W6JW30Z2GsncdXPZVWrVBgmss6rmAZcCD4hItePtiMg1QUKXv8m6OzgyZWX7G7iPGuU3lkR38snu9cMPvYaRyCxRMwc1vWA6A9oNQCqsB/36NrLTSMory5m1dpbvUEzN1gIdw+Y7BMvqRFXXBq8rgA+BATVs96iq5qlqXnZ29uFHa2DOHCguhp49Idf6foyoE05w7dRmzbJ2ahFiiZqp1d7yvcxaO4sTOln1QSQM7zAcQaydWmybDfQUka4i0gC4GKjT05si0kJEGgbTrYERwOKIRWqcUDXcKaf4jSMZZGW5zm/Lyqw/tQixRM3UatbaWZRWlDKqs1UfREKLRi3o36Y/n6y2x9tjlaqWAzcA7wFLgJdVdZGI3Cki5wKIyGARKQAuAB4RkUXB7n2AfBH5AvgvcI+qWqIWaR984F5D1XImskLn2dqpRURtY30aw8erPkaQfV1JmPo3qvMonpr3FGUVZaSnpvsOx1RDVd8G3q6y7Law6dm4KtGq+30GHB3xAM1+JSX7S3ZCDd1NZJ18Mvz1r5aoRYiVqJlafbL6E/q36U+LRtYPUaSc3OVkdpXtIr8w33coxsS/WbNgzx7o2xfatvUdTXIYORJSUyE/H4qKfEeTcCxRMzUqryzn09WfWrVnBBQXF5PVMouslln84LQfAHDyVSfvW5bVMos+/fv4DdKYeGTt06IvM9MNfF9RYcNJRYBVfZoafb7uc3aV7bJELQIqKysZ/9r4ffP/zP8nTc9pyuW/u3zfsgnnTfARmjHxLZSoWfu06Dr5ZJg5053/M8/0HU1CsRI1U6OPV30MYE98RkGX5l1Ys2MNFZUVvkMxJn7t2QPTp4MInHii72iSiz1QEDGWqJkafbz6Y3q27En7zPa+Q0l4XbK6UFZZxtqiOnfPZYypavp09zDBMcdAq1a+o0kuxx8PaWlu3M8dO3xHk1C8JWoikioin4vIm8F8VxGZKSLLRGRi0F8RItIwmF8WrO/iK+ZkUqmVfLLqEytNi5LOzTsD8M32b/wGYkw8s2pPf5o2hSFDoLJy//Bdpl74LFG7EdcnUci9wARV7QFsA64Kll8FbAuWTwi2MxE2f8N8tu3dxoldrPogGhqnN6Ztk7as2r7KdyjGxC97kMCvUIL8/vt+40gwXhI1EekAnA08HswLcArwSrDJ08C4YHpsME+w/tRgexNBU5ZPAWB0t9GeI0kenZt3ZvXO1ZRXlvsOxZj4s3Ona8yemmoDsfty2mnudcoUv3EkGF8lag8AvwYqg/lWwPagB3CAAiA0QFsusAb29RC+I9jeRNDkFZPpl92PnMwc36Ekja7Nu1JeWU5hUaHvUIyJPx9+COXlMHQoNG/uO5rkNHw4NGkCixbBWmtvW1+inqiJyDnARlWdU8/ve42I5ItI/qZNm+rzrZPOnrI9fLLqE8Z0H+M7lKTSOcu1U1u5faXnSIyJQ5Mnu9cxdt/ypkGD/dWfVqpWb3yUqI0AzhWRb4CXcFWefwOai0ioX7cOQCgdXwt0BAjWZwFbqr6pqj6qqnmqmpednR3Z3yDBfbL6E0oqSjit22m+Q0kqjdIbWTs1Yw7Xe++519NP9xtHsgslyqHE2RyxqCdqqvpbVe2gql2Ai4EPVPX7uAGLzw82uxJ4I5ieFMwTrP9AVTWKISedKcun0CC1gXV060GX5l1Ys3ONtVMz5lCsWAHLlrkqz7w839Ekt1CiNmWKewLUHLFY6kftN8BNIrIM1wbtiWD5E0CrYPlNwC2e4ksak1dMZkTHETRp0MR3KEmnS/MulFeWU7CzwHcoxsSPUOnNqae6vryMP716QadOsHkzzJvnO5qE4DVRU9UPVfWcYHqFqg5R1R6qeoGqlgTL9wbzPYL1K3zGnOjWF69n/ob51j7Nky7NuyAIK7bZZW5MnYUSNav29E/Eqj/rWSyVqBmP+vTvQ1bLLLqe1hWAu3581wEDhGe1zKK4qNhzlIkvIy2D3Ga5lqgZU1fl5fv77TrN2tXGBEvU6pWVERsACgsLGf/aeF5Z/ArfbP+GX/7zl1Ttru6O0Xd4ii65dGvRjU9WfULThk19h2JM7Js50/Wh1qsXdOniOxoDrgpaBKZNg127XJcd5rBZiZrZp6KygmVbl9GzVc9vJWkmerq36I6ilHe0BwqMOSir9ow9LVvC4MFQVub6tzNHxBI1s8/qHaspqSjhqFZH+Q4lqeVm5tIgtQHlnSxRM+agQt1yWP9psSWUOL/zjt84EoAlamafpVuWkiqpdGvRzXcoSS01JdWNUmCJmjG127gRZs2Chg3hpJN8R2PCnXOOe33zTbAetY6IJWoGAEX5astXdG3RlQapDXyHk/S6teiGZinLty73HYoxsevtt10ScPLJ0NTadMaUvDxo0wZWrXJDSpnDZomaAaCyRSXb9m6jV6tevkMxuHZqAJOX21NTxtTozTfd63e+4zcO820pKXDWWW76rbf8xhLnLFEzAJR3ddVsvVpaohYLWjZqiewQ3llm7TuMqVZJyf72aWef7TcWU73w6k9z2CxRMwCUdy+nXdN2ZGVk+Q7FACJC2so03l/5PnvK9vgOx5jY8/HHUFwMRx8NnTv7jsZU57TTID0dPvsMtnxriG5TR5aoGdYVraOifQW9W/X2HYoJk74ynd1lu/nwmw99h5L0ROQMEVkqIstE5FvD2InIKBGZKyLlInJ+lXVXisjXwc+VVfc1h+k//3GvVu0Zu5o1gxNPdGN+vvuu72jiliVqhte/fB0E+mT38R2KCZO6NpUm6U148yurNvBJRFKBh4Azgb7AJSLSt8pmq4EfAC9U2bclcDswFBgC3C4iLSIdc8JT3V+dFqpeM7HJqj+PmCVqhleXvErK1hSyG2f7DsWEkQrhtO6n8Z+v/oPa4+0+DQGWBeMRlwIvAWPDN1DVb1R1PlBZZd/TgSmqulVVtwFTgDOiEXRCW7IEVq6E1q1hyBDf0ZjahBK1d95xHeCaQ2aJWpLbvHszH33zEWnL0mw0ghh0Ts9zWLNzDQs2LvAdSjLLBdaEzRcEyyK9r6lJqNrz7LMhNdVvLKZ23btD796wY4drV2gOmSVqSe6NL9+gQitIX5buOxRTjbN6usfbrfoz8YnINSKSLyL5mzZt8h1ObHvtNfd67rl+4zB1c9557vXVV/3GEacsUUtyry55la7Nu5KyyS6FWNQ+sz2Dcwbzn6/+4zuUZLYW6Bg23yFYVq/7quqjqpqnqnnZ2dYMoUarV7vRCBo3hjOsFjkufO977vX1192DBeaQ2LdzEtu+dztTV0zle32+h2DVnrFq7FFjmVEwg4KdBb5DSVazgZ4i0lVEGgAXA5PquO97wBgRaRE8RDAmWGYOV6g07cwzXbJmYt+AAdClC6xf77rqMIfEErUk9uZXb1JWWcb3+n7PdyimFhf2uxCAfy/6t+dIkpOqlgM34BKsJcDLqrpIRO4UkXMBRGSwiBQAFwCPiMiiYN+twJ9wyd5s4M5gmTlcoeqz88+vfTsTO0T2l6pZ9echs0Qtib265FVyM3MZkmtPTcWynq16MqDdACYumug7lKSlqm+rai9V7a6qdwfLblPVScH0bFXtoKpNVLWVqvYL2/dJVe0R/PzL1++QENavh08/dYOw22gE8SWUqL32mg3SfogsUUtSxaXFvLvsXc7rcx4pYpdBrLuw34XMXDuTb7Z/4zsUY/x5/XX3JT9mDGRm+o7GHIqhQyEnx7UxzM/3HU1csW/oJPXmV2+yt3wv3+tj1Z7xwKo/jcGqPeNZSoo9/XmYLFFLUhMXTSQnM4eRnUb6DsXUQbcW3RicM9iqP03y2rQJPvwQ0tJs2Kh4Far+/Pe/rfrzEFiiloR2luzkna/f4YK+F5CaYp1FxosL+13InHVz+HrL175DMSb6Jk6Eigo4/XRoYaNwxaUTToD27WHFCpg503c0ccMStST0xpdvUFJRwkX9LvIdiqlFcXExWS2z9v3ccdEdUAn9r+i/b1mf/jY+q0kSzz/vXr//fb9xmMOXmgoXX+ymQ39Pc1BpvgMw0Tdx0UQ6ZXViWIdhvkMxtaisrGT8a+MPWPbCghdYN3IdN958IymSwoTzJniKzpgoWr4cZsyAJk1sNIJ49/3vw4QJroT0/vsh3UbFORgrUUsy2/ZsY/LyyVzY90Ib2zMODWw/kOLSYqv+NMklVPpy3nkuWTPxa+BAN/bnpk0wdarvaOKCJWpJ5vUvX6essoyL+lu1Zzzq2bInTdKb8Pn6z32HYkx0qFq1ZyIR2f93tOrPOrFELclMXDSRbi26Maj9IN+hmMOQmpLKce2O46stX1FUUuQ7HGMib84c+OoraNMGTj3VdzSmPlx6qXt9/XUoLvYbSxywRC2JbNq1ifdXvM9F/S6yas84NqDdABRl3vp5vkMxJvKefda9Xnyx65rDxL9u3WD4cNi92yVrplaWqCWRV5e8SoVWcHH/i32HYo5Aq8at6Nq8K7MLZ6Mp1heRSWB798Jzz7npK67wG4upX1de6V6feMJvHHHAErUE16d/n31dOdzw8A2kbE1hZJ+RB3T7kNUyi+IiK36OJ0Nzh1JUWkR593LfoRgTOa+/Dlu3woABMMiaaySUSy6Bxo3ho49c1bapkZUjJ7jCwkLGvzaeopIi7p9xPyd2PpGTvnvSt7a7Y/QdUY/NHL5erXrRIqMFO4/b6TsUYyLn8cfd69VX+43D1L9mzVx19pNPur/zX/7iO6KYZSVqSWLx5sUA9Mvu5zkSUx9EhKG5Q6nIqWDW2lm+wzGm/i1fDh98AI0audIXk3hCCfhTT0FpqddQYpklakli0cZFtGnShuwm2b5DMfXkuHbHQQn8bebffIdiTP0LtV264AJo3txrKCZChg6Ffv1cn2qTJvmOJmZZopYEtu3Zxpqdazi6zdG+QzH1qGFaQxosbsDLi16msKjQdzjG1J/SUvjXv9y0VXsmLpH9f9/HHvMbSwyzRC0JzN8wH8AStQTU4IsGVFRW8M/Z//QdijH159//hvXrXWnLiBG+ozGRdPnlkJEBkyfbQwU1sEQtwSnK/I3z6dK8C1kZWb7DMfUsZUcK3znqOzw852H2lO3xHY4xR04V/hZU5994oyt1MYmrZUu47DI3/eCDfmOJUVFP1ESko4j8V0QWi8giEbkxWN5SRKaIyNfBa4tguYjIgyKyTETmi8jAaMcczyraVbB1z1aOaXuM71BMhIwfNp7Nuzfz7PxnfYdizJGbORNmz3Zf4DZkVHK48Ub3+q9/wbZtfmOJQT5K1MqBX6pqX2AYcL2I9AVuAd5X1Z7A+8E8wJlAz+DnGsDqeA5BWe8y0lLS6Nu6r+9QTISc2PlEBrYfyP3T76dSK32HY8yRCZWmXX2162fLJL7+/WH0aDdSgXWA+y1RT9RUdZ2qzg2mi4AlQC4wFng62OxpYFwwPRZ4Rp0ZQHMRaR/dqONTaUUp5b3K6d2qNw3TGvoOx0SIiHDz8JtZumUpb331lu9wjDl8BQXwyiuQmgrXX+87GhNNv/iFe/3736HcOvIO57WNmoh0AQYAM4G2qrouWLUeaBtM5wJrwnYrCJaZg3j767fRRmrVnkng/L7n07FZR+6bfp/vUIw5fPfd576kv/c96NjRdzQmms48E3r2hFWr3MMkZh9viZqINAVeBX6hqgd0r66qChzSIIYico2I5ItI/qZNm+ox0vj17Pxnkd1C95bdfYdiIiw9NZ1fDPsFH636iPzCfN/hJBwROUNElgZtZW+pZn1DEZkYrJ8Z/BOKiHQRkT0iMi/4eTjqwceLTZvgkUfc9C3fOsUm0aWkwK9+5ab//GeotGYcIV4SNRFJxyVpz6vqa8HiDaEqzeB1Y7B8LRD+r1WHYNkBVPVRVc1T1bzsbOvUdduebbz51ZukL00nRezh3mTw44E/plnDZlaqVs9EJBV4CNdeti9wSdCuNtxVwDZV7QFMAO4NW7dcVY8Lfq6NStDx6IEHYM8eOOssN7anST5XXAEdOsDChdYBbhgfT30K8ASwRFXvD1s1CbgymL4SeCNs+RXB05/DgB1hVaSmBi8ufJHSilLSl6T7DsVEUHFxMVkts8hqmUXH9h3Z+9leXpr/EpmdM/ct79O/j+8w490QYJmqrlDVUuAlXNvZcOFtbF8BTg3udaYutm93bZMAbr3VayjGo4YN95eq3XWX66rFeBmUfQRwObBAROYFy34H3AO8LCJXAauAC4N1bwNnAcuA3cAPoxptHFJVHpv7GAPaDWDZpmW+wzERVFlZyfjXxu+b37F3Bw/OepD+f+jP6T1OB2DCeRN8hZcoqmsnO7SmbVS1XER2AK2CdV1F5HNgJ/B7Vf0kwvHGn7//HXbuhJNPhuOP9x2N8enHP4a774Y5c+Ddd13btSTn46nPaaoqqnpMWHXA26q6RVVPVdWeqjpaVbcG26uqXq+q3VX1aFW1BjgHMWfdHOatn8fVA69GsH/qk0lWRhb9svsxd/1c9pbv9R2OgXVAJ1UdANwEvCAizarbMGnb2W7ZAn/9q5v+wx/8xmL8a9wYfvlLN/2HP1hbNWxkgoT02JzHaJzemEuPvtR3KMaD4R2GU1pRypx1c3yHkijq0k523zYikgZkAVtUtURVtwCo6hxgOdCruoMkbTvb//kfV5o2ZowrUTPm+uuhfXtXqvbyy76j8c4StQRTXFrMCwtf4MJ+F9qQUUmqfWZ7ujTvwqy1s6iorPAdTiKYDfQUka4i0gC4GNd2Nlx4G9vzgQ9UVUUkO3gYARHphuu4e0WU4o59q1fvb5v2P//jNxYTO5o0gTvucNO33gqlpX7j8cwStQTz0sKXKC4t5uqBV/sOxXh0fIfj2Vmyk0WbFvkOJe6pajlwA/AeroPul1V1kYjcKSLnBps9AbQSkWW4Ks5Q/xKjgPlBe9xXgGtDzToMcNttUFICF10EA210QBPmhz+E3r1hxQp4OLl7tfHxMIGJEFXlH7P/Qb/sfgzvMNx3OMajHi170Lpxa6YXTEcPrUtCUw1VfRv3YFP4stvCpvcCF1Sz36u4rohMVdOnw9NPQ4MG7gk/Y8KlpcE998C4ca507dJLoXVr31F5YSVqCWTa6ml8vv5zfjbkZ1jPAMlNRBjeYTjri9dT0cGqP02MqaiAG25w0zffDD16+I3HxKZzz4VTToGtW+G3v/UdjTeWqCWQv838Gy0yWnD5sZf7DsXEgGPaHkOT9CaUDkru9h0mBj32GMyd6zo3/d3vfEdjYpUIPPQQpKfD44+7UtgkZIlaHOvTv8++Tk0zO2fy6qJX2TVtF+3btt+3vLio2HeYxpO0lDSG5A6hvEs5izct9h2OMU5h4f7Skfvvdw3HjalJ7977O8H96U+TcsB2S9TiWGFhIeNfG8/418ZzzO3HICnCdddft2/Z+NfGU6nWB00yy8vJgzK4f/r9B9/YmEhThWuucSMRnHUWnH++74hMPLj1VujSBb74wrVbSzKWqCWAkvIS5q6bS+/Wva1LDnOAxumNSV+czrPzn2V98Xrf4Zhk9+yz8NZbkJUFjz7qqraMOZjGjV3VJ7gHC+YkVx+RlqglgNmFs9lbvpcRHUf4DsXEoIafN6SsooyHZj3kOxSTzFauhJ//3E0/8ADk5noNx8SZU0911095OVx2GezZ4zuiqLFELc6VVpQyvWA6PVr2ILeZ3fjMt6XsSGFs77H8I/8f7Crd5Tsck4xKS11faTt2wNixcOWVB9/HmKruuce1WfvyS7jpJt/RRI0lanEuvzCf3WW7GdVplO9QTAy7efjNbN2zlcfmPuY7FJOMfv1rmD0bOneGJ5+0Kk9zeBo1guefh4YNXSe4Tz3lO6KosEQtjmma8tmaz+jWohsdszoefAeTtEZ0GsEpXU/h3k/vZU9Z8lQZmBjwzDPwt7+5LhYmToSWLX1HZOLZwIGuyw6Aa6913bwkOEvU4ljpgFJ2le3ixM4n+g7FxIHbT7yd9cXreWTOI75DMcnik0/gxz920w88AEOHeg3HJIirrnJPD5eUuKr0ggLfEUWUJWpxas2ONZQMLqFv6750yurkOxwTB0Z1HsUpXU/hnmn3sLtst+9wTKL78kv47nehrMw1Ar/uOt8RmUTy4IMwYoRL0k4/3Y1ekKAsUYtTv5ryKxA4rftpvkMxceSPJ/6RDbs28M/Z//Qdiklky5e7p/S2bHH9pd1v/fiZetawIUyaBH37wuLF8J3vwK7EfFjKErU49NE3HzFx0UQa5jekeUZz3+GYGFdcXLxvpIpzBpxD6qpUbv7PzTRr32zf8qyWWfTp38d3qCYRrFzpxmcsLIQTT4R//xtSU31HZRJRy5bw7rtuKLLPPnMlazt2+I6q3qX5DsAcmp0lO/nhGz+kS/MubJmzxXc4Jg5UVlYy/rXx++Y3FG/gkTmP0O9P/Ti9x+n7lk84b4KP8EwimTcPzjwT1q+H4cPhP/9xnZUaEykdO8L777sS3E8/hdGjXfLWqpXvyOqNlajFmeveuo5VO1bx3HefQ8rtEXdz6No2bcuA9gOYVTiLLbst2Tf1ZOpUGDXKJWmnnOK+LDMzfUdlkkGvXu7BlW7dID8fhg2DJUt8R1VvLFGLI8/Nf47nFzzP7SfezohONgqBOXwndzmZtJQ0pqyY4jsUE+8qK11HpKefDkVFcOGF8Pbb0KyZ78hMMunSBT7+GAYMgGXL3BPGb77pO6p6YYlanJhTOIefvPkTTuh0AreecKvvcEyca9qgKSd0OoGlW5ayeNNi3+GYeLVhg+se4be/dQnbrbfCiy+6ht7GRFturitZu+AC90/Dd77jnjiO8+GmLFGLA6u2r+KcF88hu3E2L1/wMqkp1jDXHLnjOx5P+6bteevrt2xoKXNoVF0P8X37ulKL5s1de7S77oIU+1oxHjVp4jpWvvdeSEuD//f/IC/PJXBxyj5RMW7bnm2c/cLZ7Cnbw1uXvkW7pu18h2QSRIqkMK73OErKS3j767d9h2PiRX4+nHyyGxh761YYMwa++ALOOcd3ZMY4Im7Yshkz4KijXPcdo0bB5Ze7p5HjjCVqMapP/z40a9eM1r9szaJ1iyh7oYzjex9/QHcKxUXFvsM0ca5Nkzac2OVEFm9eTOnRpb7DMbHs88/h4oth8GD46CP3VN0TT7iHBjpZp9smBg0a5K7b22931fHPPeceOPj5z2HtWt/R1Zl1zxGj1m5eS7NfN2NX8S4u6XcJvU7t9a1t7hh9h4fITKIZ0XEEa3as4euTvmbqiqmM7jbad0gmVpSWwjvvuLEVpwQPnjRoADfeCL/7navyNCaWNWoEf/wjXHEF/OY38Morrjr0kUfg/PPhpz91IxxI7PaiYCVqMWjTrk3sOm8X64rXcUHfC+jV6ttJmjH1JUVS+F6f75GyNYUL/n0BX27+0ndIxqfycvf03I03usbZ48a5JK1pUxg/Hr7+Gv7yF0vSTHzp1s11vvzFFy5BKyuDF16AE05wbS1vuw0WLHDtL2OMJWox5pvt3zDyXyOpbFnJRf0uonfr3r5DMkmgYVpDGk9qTIPUBpz89Mks2LDAd0gxRUTOEJGlIrJMRG6pZn1DEZkYrJ8pIl3C1v02WL5URE6vum/MmDjRVW1mZ7sRBR58EDZvhn794K9/hdWr3VBQVs1p4tkxx7iEbeVKVyrcpo0bl/ZPf3Lrund3pWwffeQ70n0sUYshn6z6hGGPD2Pjro00fr2xlaSZqEopSuHDKz8kRVI48akTmbV2lu+QYoKIpAIPAWcCfYFLRKRvlc2uArapag9gAnBvsG9f4GKgH3AG8I/g/WLPk0+6ZG37dtcA+6ab3IMDCxbAzTdDixa+IzSm/nTuDHff7QZ1nzwZrrkGWrd2CdzDD7sOnGOEJWoxQFV5YMYDnPz0yTRr2IxPf/Qpaeus+aCJvj7ZfZj2w2m0aNSCk546iefmP+c7pFgwBFimqitUtRR4CRhbZZuxwNPB9CvAqSIiwfKXVLVEVVcCy4L3iz3XXw8PPOCqNr/8Eu67zzXGjuG2O8YcsfR0OO0012Zt/XqYOdOVrn3ve74j28eyAc+Wb13ONW9ewwcrP2DsUWN5etzTZGVk+Q7LJKHQ4O0AlY0rKT2zlMtfv5wf3f4jMqZlIBVCTk4OSxYmztAsdZQLrAmbLwCG1rSNqpaLyA6gVbB8RpV9cyMX6hE491zfERjjV2oqDBnifmKIJWqeFJUUcf/0+7n303tJT03nn2f/k2sGXUOKWCGn8aPq4O0VlRVMXTmVGcyg+fHNGdd7HBOvnOgxwsQmItcA1wB0snZgxpiAJWpRtrNkJz0u6cHmnpvRxkra12mkfpzKb/7yG37Db/ZtZ32kGd9SU1I5vfvpdG/RnUlLJ/H43MdJH5ZOaUUpDVIb+A4vmtYCHcPmOwTLqtumQETSgCxgSx33BUBVHwUeBcjLy4u9R8+MMV5YohYla3as4cGZD/Lo3EfZeexOujbvyqldTyX3xFz48be3tz7STKzo0bIH1w2+jveWvce8ofM49uFjuX/M/ZzZ80zfoUXLbKCniHTFJVkXA5dW2WYScCUwHTgf+EBVVUQmAS+IyP1ADtATsKc0jDF1ZvVsEaSqTF8zncteu4xuD3ZjwowJnNXzLJq82IQrjr2C3Gax2VTFmKoy0jIY23ssjd5oREVlBWe9cBanP3c6n67+1HdoEaeq5cANwHvAEuBlVV0kIneKSKhh1xNAKxFZBtwE3BLsuwh4GVgMvAtcr6oV0f4djDHxy0rUImDrnq08N/85Hpv7GAs3LqRpg6b8bMjPuHHojXRu3pmsq+1hAROf0r9JZ+F1C/l/M/8f93x6DyP/NZITOp3AT/N+yrje42iU3sh3iBGhqm8Db1dZdlvY9F7gghr2vRu4O6IBGmMSliVq9WTjro0ce/GxbGq9iYqOFZAKKRtSyFiQgXwlPFH2BE/wBGDtz0x8a5DagF8e/0uuzbuWx+c+zoQZE7j0tUvJapjF2N5jObvn2YzpPobmGc19h2qMMXEvbhI1ETkD+BuQCjyuqvf4ikVVWbVjFTMKZjBt9TSmrZ7G/A3z0TylRUYL+rTuQ/82/Wmf2R4u/Pb+1v7MxKvwLjxCFKVxx8bs7rObZ3Y8wzNfPEOqpDKi0whGdRrFCZ1PYHiH4WQ2zPQUtTHGxK+4SNTCegY/DdcP0WwRmaSqiyN53PLKclbvWM3yrctZumUpCzYsYMHGBSzcuJCi0iIAmjZoyvAOw7njpDu498f38rNHf4ZYB5EmQVXtwuNb67WS+264jxsfvJHJyyfz52l/pvKTSlIllV6tetGvTT/6tu5Lr1a9aNe0HW2atCG7STatGrUiPTU9ir+JMcbEh7hI1AjrGRxAREI9gx9xovbhNx+yYtsKNu7auO9nffF6Vm5fyartq6gIb/e7F1I3p5KyJYWMzRmkbkxFNgkzdSYzmcmeoj2WpJmkliIp7P16Lw9d8BAATdKbUNG+gvKccr5u/TVLWy7llaxXqn+MqQRkrxzwk7sml28++Saqv4MxxsSSeEnU6tIz+GG59YNb+WzNZ26mFGSPkLI7BdkppO5IJX1HOik7UthTsIc/vPaHWhMxq9I05uClbmUVZfz5+3/myoevZHfZbnaV7mJ3+W72lO1hT/ke91q2h93lu9m0cFMUIzfGmNgjqrHfr6KInA+coao/DuYvB4aq6g1h2+zr1Rs4Clga9UCd1sBmT8eOhePHQgzJfvxYiMHH8TuranaUjxkRIrIJWOXh0L6vm1iIIdmPHwsxJOPxa7x/xUuJ2kF79w7v1dsnEclX1bxkPX4sxJDsx4+FGHwfP975Sjhj4e/mO4ZkP34sxJDsx68qXjq83dczuIg0wPUMPslzTMYYY4wxERUXJWqqWi4ioZ7BU4Engx6/jTHGGGMSVlwkalB9z+Axynf1q+/jg/8Ykv344D8G38c3hycW/m6+Y0j244P/GJL9+AeIi4cJjDHGGGOSUby0UTPGGGOMSTqWqEWAiPxJROaLyDwRmSwiOVE+/l9F5MsghtdFpHmUj3+BiCwSkUoRieqTMyJyhogsFZFlInJLlI/9pIhsFJGF0Txu2PE7ish/RWRxcP5v9BBDhojMEpEvghisc8E4k+z3ryAGL/cwn/ev4PhJfQ+L1fuXVX1GgIg0U9WdwfTPgb6qem0Ujz8G+CB4CONeAFX9TRSP3weoBB4BblbV/CgdNxX4irChxoBLIj3UWNjxRwHFwDOq2j8ax6xy/PZAe1WdKyKZwBxgXLR+/yAGAZqoarGIpAPTgBtVdUa0YjBHJtnvX0EMUb+H+b5/BTEk9T0sVu9fVqIWAaGbXKAJENVsWFUnq2p5MDsD1+9cNI+/RFV9dDi8b6gxVS0FQkONRYWqfgxsjdbxqjn+OlWdG0wXAUtwo3pEMwZV1eJgNj34sf8G40iy37+CGHzcw7zev8DuYbF6/7JELUJE5G4RWQN8H7jNYyg/At7xePxoqm6osagmKrFCRLoAA4CZHo6dKiLzgI3AFFWNegzmyNj9ywu7f4XxdQ+LxfuXJWqHSUSmisjCan7GAqjqraraEXgeuKH2d6v/4wfb3AqUBzFE/fjGDxFpCrwK/KJK6UhUqGqFqh6HKwkZIiJRr0IxtUv2+1ddYzB++LyHxeL9K276UYs1qjq6jps+j+v/7fZoHl9EfgCcA5yqEWiIeAi/fzQddKixRBe0q3gVeF5VX/MZi6puF5H/AmcAXhonm+ol+/2rLjF4kPT3L4ide1gs3b+sRC0CRKRn2OxY4MsoH/8M4NfAuaq6O5rH9iyphxoLGsI+ASxR1fs9xZAdekpPRBrhGkZH9fo3R8buX94k9f0L/N/DYvX+ZU99RoCIvAochXtqaBVwrapG7T8jEVkGNAS2BItmRPmpre8C/w/IBrYD81T19Cgd+yzgAfYPNXZ3NI4bHPtF4CSgNbABuF1Vn4ji8UcCnwALcNcewO+CUT2iFcMxwNO4858CvKyqd0br+ObIJfv9K4jByz3M5/0rOH5S38Ni9f5liZoxxhhjTIyyqk9jjDHGmBhliZoxxhhjTIyyRM0YY4wxJkZZomaMMcYYE6MsUTPGGGOMiVGWqJmkISIfikie7ziMMcaYurJEzRhjjDEmRlmiZrwSkSYi8paIfBGMs3eRiNwmIrOD+UeD3qpDJWITRCRfRJaIyGAReU1EvhaRu4JtuojIlyLyfLDNKyLSuJrjjhGR6SIyV0T+HYwth4jcIyKLRWS+iPxvdM+GMSYZBfey+SKSEdwTF8XCGJMmNliiZnw7AyhU1WNVtT/wLvB3VR0czDfCjfkXUqqqecDDwBvA9UB/4Aci0irY5ijgH6raB9gJXBd+QBFpDfweGK2qA4F84KZg/+8C/VT1GOCuyPzKxhizn6rOxg0XdRfwF+A5VbXxcQ1giZrxbwFwmojcKyInqOoO4GQRmSkiC4BTgH5h208K22+Rqq5T1RJgBfsHNF6jqp8G088BI6sccxjQF/hUROYBVwKdgR3AXuAJETkPSKZxBo0xft2JG1syD5esGQNAmu8ATHJT1a9EZCBwFnCXiLyPKyXLU9U1IvJHICNsl5LgtTJsOjQfup6rjotWdV6AKap6SdV4RGQIcCpwPnADLlE0xphIawU0BdJx97xdfsMxscJK1IxXIpID7FbV54C/AgODVZuDdmPnH8bbdhKR4cH0pcC0KutnACNEpEcQQxMR6RUcLysYAHg8cOxhHNsYYw7HI8AfgOeBez3HYmKIlagZ344G/ioilUAZ8FNgHLAQWA/MPoz3XApcLyJPAouBf4avVNVNIvID4EURaRgs/j1QBLwhIhm4UrebDuPYxhhzSETkCqBMVV8QkVTgMxE5RVU/8B2b8U9Uq9YKGRO/RKQL8GbwIIIxxhgT16zq0xhjjDEmRlmJmjHGGGNMjLISNWOMMcaYGGWJmjHGGGNMjLJEzRhjjDEmRlmiZowxxhgToyxRM8YYY4yJUZaoGWOMMcbEqP8PfTCxFk/zkk4AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 720x360 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"samples = get_gaussian_samples(10000)\n",
"\n",
"fig, ax = ploter.subplots(nrows=1, ncols=2, sharex=True, figsize=(10,5))\n",
"\n",
"histplot(ax= ax[0], data=samples, bins=25, color='g', kde=True)\n",
"ax[0].set(xlabel='samples', ylabel='frequency',title='Histogram of the generated gaussian samples')\n",
"\n",
"x = linspace(-3,3,1000)\n",
"mu, sigma = 0, 1\n",
"ax[1].plot(x,(1/(sigma * sqrt(2 * pi))* exp( - (x - mu)**2 / (2 * sigma**2))), linewidth=2, color='r')\n",
"ploted = ax[1].set(xlabel='x', ylabel='f(x)',title='theorical density of the gaussian law ')\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The two graphs are pretty close in shape.\n",
"We have well generated randm variables. That's a graphical way to verify our result.\n",
"\n",
"\n",
"We will now use the previous function to estimate moments of a guassian variable X$\\sim$ N( $\\mu$ = 125,$\\sigma^2$ = 100)\n",
"\n",
"The first moments have their exact value given: E[X] =125, E[X$^2$] =15725\n"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [],
"source": [
"\"\"\"\n",
"Compute an approximation of gaussian law's N(sigma, mu) moments \n",
"A linear transformation help gets sample from the normal law N(0,1).\n",
"\n",
"@param order: the order of the moment\n",
"@param observation: the number of samples generated\n",
"@param repetition: number of oparations to generate an average moment\n",
"\n",
"@return moment: the moment of the gaussion variable for the given order\n",
"\n",
"\"\"\"\n",
"def calc_gaussian_moment(order, observation = 1, repetition = 1):\n",
" mu, sigma = 125, 10\n",
" moments = [0 for x in range(repetition)]\n",
" for i in range(repetition):\n",
" normal_samp = get_gaussian_samples(observation)\n",
" actual_samp = sigma * normal_samp + mu\n",
" for j in range(observation):\n",
" moments[i] += (1/observation) * (pow(actual_samp[j],order)) \n",
" \n",
" return mean(moments) \n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can compute the first two moments and compare them with the real values."
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [],
"source": [
"real_moment_1 = 125\n",
"real_moment_2 = 15725\n",
"\n",
"computed_moment_1 = calc_gaussian_moment(order=1,observation=10000, repetition=100)\n",
"computed_moment_2 = calc_gaussian_moment(order=2,observation=10000, repetition=100)\n",
"\n",
"\n",
"relative_err_1=100* abs(real_moment_1 - computed_moment_1)/real_moment_1\n",
"relative_err_2=100* abs(real_moment_2 - computed_moment_2)/real_moment_2"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(125, 124.99187088760463, 0.006503289916292942)"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"real_moment_1, computed_moment_1, relative_err_1"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(15725, 15725.074837779624, 0.0004759159276593404)"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"real_moment_2, computed_moment_2, relative_err_2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As swe can see the _relative error is less than 0.005%!_ Sounds great!\n",
"\n",
"We can rely on the rejection method to approch certain values, with a small risk of error.\n",
"\n",
"\n",
"Estimating the fifth moment gives the result below:"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"32508844761.31357"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"computed_moment_5 = calc_gaussian_moment(order=5,observation=100, repetition=100)\n",
"computed_moment_5"
]
},
{
"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.8.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
# Let's estimate gauss integral via Monte Carlo' method.

Using python capabilities, we are simulating a set of 6 mathematical problems resolvable with Monte carlo's methods. They take advantage of the strong law of big numbers. The work on these examples is on progress. Fell free to leave a comment if you want.

# Monte carlo's method appliyied to Markov's chains
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment