Lec_12_monte_carlo2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Physics 300 \n",
    "## Computational Physics I (Fall 2018)\n",
    "## BPB-248, Mon/Wed 02:30 - 03:45 pm \n",
    "\n",
    "|Instructor| Prof. Qiang Zhu|\n",
    "|--|-------------------------------|\n",
    "|Email | qiang.zhu@unlv.edu|\n",
    "|Website|http://www.physics.unlv.edu/~qzhu/|\n",
    "|Office| BPB 232|\n",
    "|Office hours | Mon/Wed 03:45 - 05:00 pm |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 12 Monte Carlo\n",
    "\n",
    "We have introduced the Monte Carlo methods in the last lecture. It is not a very good one for doing integrals.\n",
    "There are significantly better ways to perform the calculation, the most common of which is the _mean value method_.\n",
    "\n",
    "## 12.1 The Mean value method\n",
    "\n",
    "Let's take a general integration problem: \n",
    "\n",
    "$$ I = \\int_a^b f(x)dx $$\n",
    "\n",
    "The average value of $f(x)$ in the range from $a$ to $b$ is by definition\n",
    "\n",
    "$$ <f> = \\frac{1}{b-a}\\int_a^b f(x)dx = \\frac{I}{b-a} $$\n",
    "\n",
    "Thus\n",
    "\n",
    "$$ I = (b-a)<f> $$\n",
    "\n",
    "If we can estimate <f> then we can estimate $I$. A simple way to estimate <f> is just to measure $f(x)$\n",
    "at $N$ points $x_1, x_2 .... x_N$ chosen uniformly at random between $a$ and $b$, and then calculate\n",
    "$<f> \\approx N^{-1} \\sum^N_{i=1} f(x_i)$,  \n",
    "\n",
    "$$ I \\approx \\frac{b-a}{N} \\sum^N_{i=1} f(x_i)$$\n",
    "\n",
    "This is the fundamental formula for the mean value method. As with the previous method, it gives only an approximate\n",
    "estimate of the integral. How accurate is it? We can estimate the error using standard results for the behavior\n",
    "of random variable which tell us that the variance on the sum of $N$ independent random numbers - no matter what\n",
    "their distribution - is equal to $N$ times the variance on a single one.\n",
    "\n",
    "$$ f = <f^2> - <f>^2 $$\n",
    "in which\n",
    "$$ <f> = \\frac{1}{N} \\sum_{i=1}^N f(x_i)$$\n",
    "$$ <f^2> = \\frac{1}{N} \\sum_{i=1}^N [f(x_i)]^2$$\n",
    "\n",
    "Then the variance on the sum is $N$ times the variance on a single term. Finally, we get the standard deviation on $I$,\n",
    "\n",
    "$$ \\sigma = \\frac{b-a}{N} \\frac{\\sqrt{\\text{var}f}} {\\sqrt{N}} $$\n",
    "\n",
    "Thus, the errors goes as 1/$\\sqrt{N}$. However, the leading constant is smaller, which means the method is always\n",
    "more accurate than the previous method.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MC intergral\n",
    "import matplotlib.pyplot as plt\n",
    "from math import sin\n",
    "from random import random\n",
    "\n",
    "def f(x):\n",
    "    return sin(1/x/(2-x))**2\n",
    "\n",
    "def MC1(N):\n",
    "    count = 0\n",
    "    for i in range(N):\n",
    "        x = 2*random()\n",
    "        y = random()\n",
    "        if y<f(x):\n",
    "            count += 1\n",
    "    I = 2*count/N\n",
    "    return I\n",
    "\n",
    "\n",
    "def MC2(N):\n",
    "    I = 0\n",
    "    for i in range(N):\n",
    "        x = 2*random()\n",
    "        I += f(x)\n",
    "\n",
    "    return I/N*2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.472\n",
      "1.4606935724606758\n"
     ]
    }
   ],
   "source": [
    "print(MC1(1000))\n",
    "print(MC2(1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12.2 Integrals in high dimensions\n",
    "\n",
    "In addition to the intergration of pathological functions, Monte Carlo integration is used for performing\n",
    "high dimensional integrals. As we saw in the previous section, performing an integral over two variables\n",
    "by standard grid of points such as trapezoidal rule requires us to take samples on a two-dimensional grid.\n",
    "If we have an integral over four dimensions or even higher, then the number of points on the grid can \n",
    "become very slow and impractical. Instead, Monte Carlo in high dimension is an obvious extension of the \n",
    "technique in one dimension. The integral of a function $f(r)$ over a volume $V$ in a high-dimensional space\n",
    "is given by the generalization:\n",
    "\n",
    "$$ I = \\frac{V}{N} \\sum_{i=1}^{N}f(r_i)$$\n",
    "\n",
    "where the points $r$ are picked uniformly at random from the volume $V$.\n",
    "\n",
    "An important application of this type of integral is in financial mathematics. This is widely used for \n",
    "predicting the values of stocks, which requires the evaluation of integrals over many variables. Such\n",
    "information could be used to instruct the quick trading decisions. This is why stastical physicists\n",
    "were popular in Wall Street!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12.3 Importance sampling\n",
    "\n",
    "Monte Carlo is useful for integrating pathological functions. However, it is still not good for some cases. \n",
    "In particular, it gives problems if the function contains a divergence. For example,\n",
    "\n",
    "$$ I = \\int_0 ^1 \\frac{x^{-1/2}}{e^x + 1} dx $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/qiangzhu/.pyenv/versions/3.7.0/Python.framework/Versions/3.7/lib/python3.7/site-packages/ipykernel_launcher.py:10: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  # Remove the CWD from sys.path while we load stuff.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAGuhJREFUeJzt3Xt4XPV95/H3V3PVzOhmSbYly3dsjLkbl0so4ZItoaS5tLRZkidtmk3CQvNs0m03abt9+myzfXaf3T+abpPSzZImTXMHQkNoArtpAsRAjB0bbMAYbPAF27ItWffRSDMa6bd/zMhXyRoZzZxzNJ/X8+g5Z86cGX1/Hunjn37nd84x5xwiIhIcNV4XICIis6PgFhEJGAW3iEjAKLhFRAJGwS0iEjAKbhGRgFFwi4gEjIJbRCRgFNwiIgETLsebtrS0uBUrVpTjrUVE5qXt27efcM61lrJvWYJ7xYoVbNu2rRxvLSIyL5nZwVL31VCJiEjAKLhFRAKmpKESMzsADAHjQN45t7GcRYmIyPRmM8Z9q3PuRNkqERGRkmioREQkYEoNbgf8xMy2m9k95SxIRETOr9Shkl91zh0xs4XAv5rZa865TafvUAz0ewCWLVs2x2WKiMikknrczrkjxWUX8APg2in2ecA5t9E5t7G1taQ55Of44s/28vM93Rf0WhGRajFjcJtZ0szqJteB24FXylHMA5v2sUnBLSJyXqUMlSwCfmBmk/t/xzn3f8tRTCIaYjibL8dbi4jMGzMGt3NuH3BlBWohFQuTVnCLiJyXr6YDJmNh9bhFRGbgs+AOMZwd97oMERFf81Vwa6hERGRmvgruRDRMJqfgFhE5H18FdzIWJq2hEhGR8/JVcKdimg4oIjITXwV3MhZmZGyc8QnndSkiIr7lq+BOxQrTyoc1zi0iMi1fBXdyMrg1XCIiMi1fBXciGgLQXG4RkfPwVXCn1OMWEZmRr4JbQyUiIjPzVXBP9rh19qSIyPR8FdxJzSoREZmRz4K7cHBSZ0+KiEzPX8EdLfS4MxoqERGZlq+COxENYaaDkyIi5+Or4DYzklFdaEpE5Hx8FdwweTMF9bhFRKbjw+AOk9asEhGRafkuuFO676SIyHn5LrgT0RAZjXGLiEzLd8Gt+06KiJyf74I7GQvrzEkRkfPwZ3Crxy0iMi3fBbeGSkREzs93wZ2MhhkdmyA/PuF1KSIivuS/4C5eaCozppklIiJT8WFw62YKIiLno+AWEQkY3wV3StfkFhE5L98F9+Q1udXjFhGZmv+CW/edFBE5L98Gt3rcIiJTKzm4zSxkZi+a2Y/KWdDkdMDhnMa4RUSmMpse92eA3eUqZFJKPW4RkfMqKbjNrAN4D/AP5S0HaiMhanTfSRGRaZXa4/5fwOeAac9DN7N7zGybmW3r7u6+4IJO3XdSwS0iMpUZg9vMfgPocs5tP99+zrkHnHMbnXMbW1tb31ZRukKgiMj0Sulx3wi8z8wOAN8DbjOzb5WzqEQsxLBOwBERmdKMwe2c+zPnXIdzbgVwN/Ckc+4j5SwqpZspiIhMy3fzuKFw9qSGSkREpjar4HbOPe2c+41yFTMpGQvrWiUiItPwZY87FQupxy0iMg1fBndCs0pERKbly+DWfSdFRKbny+BORsNk87rvpIjIVPwZ3LrQlIjItHwZ3LrQlIjI9HwZ3Lomt4jI9Hwa3JP3nVRwi4iczZ/BffK+kxrjFhE5mz+De3KoRNcrERE5hy+DWwcnRUSm58vg1sFJEZHp+TK4J3vcutCUiMi5fBnc8UiN7jspIjINXwa37jspIjI9XwY3FMa5M5pVIiJyDt8Gd0NthL7MmNdliIj4jm+De3FDnKMDI16XISLiO74N7vbGOMcGRr0uQ0TEd3wb3G0NtZxI58jmNSVQROR0Pg7uOIB63SIiZ/FxcNcC0Nmv4BYROZ1/g7ux0OPWAUoRkTP5Nrjbiz3uoxoqERE5g2+DuzYaojERUY9bROQsvg1uKIxzH9UYt4jIGXwe3HE6NVQiInIG3wf3MQ2ViIicwdfB3d5YS19mjJGcTsIREZnk6+CePAlHByhFRE7xeXBrSqCIyNl8HtyFHndnv3rcIiKTfB3ci3W9EhGRc8wY3GYWN7OtZrbTzHaZ2ecrURhAPBKiORnVlEARkdOES9gnC9zmnEubWQR41syecM49X+bagMI1S3RwUkTklBl73K4gXXwYKX65slZ1Gp09KSJyppLGuM0sZGY7gC7gX51zW8pb1imFsyfV4xYRmVRScDvnxp1zVwEdwLVmdtnZ+5jZPWa2zcy2dXd3z1mBbQ21DI3mSWd1x3cREZjlrBLnXD/wFHDHFM894Jzb6Jzb2NraOlf10d44ObNEvW4REShtVkmrmTUW12uBXwNeK3dhk3QnHBGRM5Uyq6QN+CczC1EI+oeccz8qb1mnfXOd9i4icoYZg9s59xJwdQVqmdKi+jhmOu1dRGSSr8+cBIiGa2hJxTQlUESkyPfBDdCuKYEiIicFIrhXt6bYc3zI6zJERHwhEMG9vr2e44NZTqSzXpciIuK5YAR3Wz0Au48OelyJiIj3AhHclxSD+9VOBbeISCCCuykZpb0hzqvqcYuIBCO4Ada3N7BLPW4RkSAFdz37utO647uIVL3gBHdbPRMOXte0QBGpcoEJ7kvbdYBSRAQCFNwdTbXUxcK8enTA61JERDwVmOA2My5pr1ePW0SqXmCCGwrj3K8dG2J8omK3vBQR8Z1gBXd7PZncOAd7hr0uRUTEM8EK7uIZlJrPLSLVLFDBvWZRinCN6QxKEalqgQruWDjERQtTOkApIlUtUMENcPmSBnYe7tcBShGpWoEL7hsvaqE/M8auTs3nFpHqFMjgBnhm7wmPKxER8Ubggru1Lsb6tno27en2uhQREU8ELrgBblrbwgtv9ZHO5r0uRUSk4gIZ3DevaWVs3LFlX4/XpYiIVFwgg/uaFU3EIzUa5xaRqhTI4I6FQ1y/qplNezXOLSLVJ5DBDXDTmlb2dQ9zuC/jdSkiIhUV2OB+55rCtMBnNVwiIlUmsMF90cIUi+vjGucWkaoT2OA2M25a08Ize7vJ5Se8LkdEpGICG9wAd17exuBonqde7/K6FBGRigl0cN+0poWWVIxHth/2uhQRkYoJdHCHQzV84Kp2nnq9i97hnNfliIhUxIzBbWZLzewpM3vVzHaZ2WcqUVip7rqmg7Fxx7/s7PS6FBGRiiilx50H/tg5tx64HviUma0vb1mlu6Stnkva6vnnFzRcIiLVYcbgds4ddc69UFwfAnYDS8pd2GzctWEJOw8P8EbXkNeliIiU3azGuM1sBXA1sKUcxVyo91+1hFCN8f3tR7wuRUSk7EoObjNLAY8Af+icO+emj2Z2j5ltM7Nt3d2VvYZIa12Mm9e28uiLR3RLMxGZ90oKbjOLUAjtbzvn/nmqfZxzDzjnNjrnNra2ts5ljSX54MalHBsc5fGXj1b8e4uIVFIps0oM+Cqw2zn3hfKXdGFuX7+IVa1J/v7pN3FOvW4Rmb9K6XHfCPwucJuZ7Sh+3Vnmumatpsa49+bV7D46yNO6rZmIzGOlzCp51jlnzrkrnHNXFb8er0Rxs/WBq5bQ1hDnfz/1pteliIiUTaDPnDxbNFzDJ29axdYDvWw70Ot1OSIiZTGvghvg7muX0pSI8PdPq9ctIvPTvAvuRDTMx25cyZOvdbHzUL/X5YiIzLl5F9wAH7txBS2pGJ//l11MaF63iMwz8zK46+IRPnfHxbzwVj+P7tDZlCIyv8zL4Ab47Q0dXNnRwP944jXS2bzX5YiIzJl5G9w1NcZfvu9Suoay3P/UG16XIyIyZ+ZtcANcvayJuzZ08NVn9vNmd9rrckRE5sS8Dm6AP/n1i0nEQvzRQzsZG9dNhUUk+OZ9cC+si/Pff/Nydh7q50tPashERIJv3gc3FO4G/1sblvB3T+5l+8E+r8sREXlbqiK4AT7/vktpb6zlPz64Q7NMRCTQqia46+IRvvDBqzjcl+GzD+/UiTkiElhVE9wA165cwJ/9+iU88cox/vZne70uR0TkgoS9LqDSPnHTSl4/PsTf/mwvaxfV8Z4r2rwuSURkVqqqxw1gZvy337yMDcsa+eOHd/Dy4QGvSxIRmZWqC26AWDjEl3/3GpqTMT76j1t5o2vI65JEREpWlcENhfnd3/rEdYRqjA9/ZQsHe4a9LklEpCRVG9wAK1uSfOvj1zE2PsGHv7KFzv4Rr0sSEZlRVQc3wMWL6/jGv7uOwZExfufLm9mna5qIiM9VfXADXN7RwHc+eT2jY+P8zpc364CliPiagrvo8o4GHr73BuKREB/6yvM898YJr0sSEZmSgvs0q1pTPHLfO2hvjPPRr23lm5sP4JzOsBQRf1Fwn2VxQ5xH7nsHN69t5S9+uIv//IOXyeV1OVgR8Q8F9xTq4hEe+L2NfOrW1Xx36yHufmAzh3ozXpclIgIouKcVqjE+++513P/hDew9nubOLz7D4y8f9bosEREF90zec0UbP/70TaxuTfEH336BP33kJYZGx7wuS0SqmIK7BMuaEzx87w3cd8tqHtp2iHf/zSZ+vqfb67JEpEopuEsUCdXwJ3es45H73kEiFuajX9vKf3p4Jz3prNeliUiVUXDP0tXLmvjxp3+VT926mkdfPMJtf/1zvrn5AOO6MYOIVIiC+wLEwiE+++51PPGZm7i0vZ6/+OEu3vulZ/mFTtoRkQpQcL8NaxbV8e1PXMf9H97AwMgYH/6HLfz+P27ltWODXpcmIvOYlePMwI0bN7pt27bN+fv62ejYON/cfJAvPbmXoWye917RzqfftYaLFqa8Lk1EAsDMtjvnNpay74w9bjP7mpl1mdkrb7+0+SseCfHJd65i0+du5d+/czU/3X2c2//m5/zh915kz3HdqEFE5s6MPW4zeyeQBr7hnLuslDetxh732XrSWb7yzH6+sfkAmdw471q3kHtvWc3G5U2YmdfliYjPzKbHXdJQiZmtAH6k4J69vuEc39h8kK//Yj99mTGu7GjgYzeu5M7L24iGdYhBRAoU3D40khvn4e2H+PpzB9h3YpjWuhgf+pWl/Ntrl7Gksdbr8kTEY54Et5ndA9wDsGzZsmsOHjxYUrHVZmLCsWlvN1//xQF+vqcbA265eCF3/8pSbl23kEhIvXCRaqQed0Ac6s3w4C8P8eC2Q3QPZWlORvnA1Uv4rQ1LWN9Wr7FwkSqi4A6YsfEJNu3p5vvbD/PT3ccZG3dctDDFB65q571XtrO8Oel1iSJSZnMa3Gb2XeAWoAU4DvwX59xXz/caBfeF6xvO8fgrR/nhi51sPdALwKXt9dx5eRt3Xt7GyhaFuMh8NOc97tlScM+NI/0jPP7SUX788lF2HOoHYO2iFLevX8ztly7isvYGamo0nCIyHyi456Ej/SP8v1eO8ZNXj7F1fy8TDlrrYtx28UJuXbeQGy9qpi4e8bpMEblACu55rm84x5OvdfHk611s2tPN0GiecI2xYXkTN69t5aY1LVza3kBIvXGRwFBwV5Gx8Qm2Hehj095uNu3pZldn4QJXjYkIN6xq5h2rm7lhdTOrW1OapSLiYwruKtY9lOUXb57g2b0neO6NE3QOjAKFYZVrVy7gupULuHblAtYurNP4uIiPKLgFAOccb/Vm2PxmD5v39bBlXy/HBgtBXh8Ps2F5ExuXN7FheRNXdjSSjIU9rlikes0muPWbOo+ZGcubkyxvTnL3tctwznG4b4Qt+3vZfrCXbQf6ePr1wr0zawwuXlzP1csauaqjkauWNbK6NaVxchEfUo+7yvVncrx4qJ8X3+rnxbf62HGon6HRPACJaIjL2hu4bEkDV3Q0cNmSela2KMxFykFDJXLBJiYc+3uG2Xmon52H+nn5yAC7OgfJ5icAqI2EuKStjkvbG7ikrZ717fVcvKiO2mjI48pFgk3BLXMqPz7BG91pdh0Z5JXOAXYdGWT30UGGsoWeuRmsbE6yrq2OtYvqWLe4sFzenFTvXKREGuOWORUO1bBucT3rFtdz1zUdACfHy3d1DvDasSF2Hx1kV+cgT7xyjMm+QCxcw+rWFGsWpVi7qO7k+vIFCcK6CqLIBVNwywUxM5YuSLB0QYI7Lms7uT2Ty7P3eJrXjw2xt2uIPcfT/HJ/Lz/c0Xlyn0iocNB0dWuSVa0pVrWcWjYlo140RyRQFNwypxLRMFcubeTKpY1nbE9n87zZlWZvV5p93Wne7E7zRleaJ1/rYmz81HBdQ22ElS1JVrYkWdGcZEVLojAzZkGCxkREJxGJoOCWCknFpg70/PgEh/tG2Hcizb7uYQ70DLP/xDBb9vXw6I4jnH4Ipi4eZnlzguULkixdkGDZggRLF9SybEGC9sZa3YRCqoaCWzwVDtWwoiXJipYkt60787nRsXEO9WY40JPhYM8wb/VmONiTYffRQX7y6rEzeuo1Bovr43Q0JehoqqWjqZYlTbV0NCVY0ljL4oY48Yhmvsj8oOAW34pHQqxZVMeaRXXnPDc+4Tg2OMqh3kzhq2+Ew30ZDveOsHlfD8cHR5k4a8JUSyrGksY47Y21tDfW0tYQP2PZkoppFowEgoJbAilUYyxprGVJYy3Xr2o+5/mx8QmODYxyuG+EI/0jdPaPcKRvhM6BEfYcH+Lp17sZGRs/5z0X1cVY3BBncUOcRfVx2orLRfVxFheXmrMuXlNwy7wUCdWcnPUyFeccAyNjdPaPcmxwpLAcGKVzYIRjA6O8dqwQ7pnc+DmvrY+HWVQfZ2F9jIV1py3rYrTWxVhYF2NhfZxkNKSDqVIWCm6pSmZGYyJKYyLK+vb6KfdxzjGUzdM1OMqxgSzHBkc5PjhK1+AoxwezdA2NsnV/L91DWXLjE+e8vjYSorUY5i2pKC2pyfVYcb2wrTkVU8jLrCi4RaZhZtTHI9THI1y08Nxx9knOOfozY3QNFcK8eyhL91CWrqEsJ9KF9X3dw2zd30tfZmzK94hHamhOFgK+ORVjQTJKcypKczLKgmSsuIye3J6I6le3munTF3mbzIymZJSmZJSLF08f8FAYe+9J5ziRLoT65HrP8OS2HF1Do+w+OkhPOjdlTx4KQb8gEWVBKkpTohDoTYnJ9UihnkSUxkTk5HOaVTN/KLhFKigSqjl58HMmzjnS2Tw96Rw9wzl6h3P0DhdCvm84R+/wGL3DWfoyY7zVm6E3nTt5/ZipxCM1NNYWwnwy1Bsnl7URGhMRGmqjNJxcLyxrIxrG8RsFt4hPmRl18Qh18QgrWpIlvSaXn6B/JEff8Bh9mRz9mULA92VyDIyM0Tecoy8zxsBIjr1dafqL20+fE3+2aKiG+tow9bWFMG+oLQR9Q23k5Lb6eGG9vjZMffzUtlQ8rCmWZaDgFplHouGa4gyXmXv0k5xzZHLjJ8N9IDNG/8hYYX1kjP5MYTlYfNyTzrGve7iwbXSMmS4wWhcLUxcvBH9dvBDspz+ui5+1jJ2+LUwyGtZt9s6i4BapcmZGMhYmGQvT0TS7105MONK5PAOZMYZG8wyOngr5wdH8ybAfGs0zVHzu6MAoe7omt+UZP/tMqXPqg1Q0TCoeJhU7tSyEfOTU9uJzyVgh/AttClEXi5CMhUjGwsTCNfNi2EfBLSIXrKbm1MybCzHZ258M9sHiMp3Nky4G+1C2sG1oNM9wNk86m2dwNE9n/wjD2XGGRscYnmK+/VTCNVYI92gh6BOxEKlY4fFk0CdjYZLRyWVheyIWIhkNk4ieej4RC1MbCXkyFKTgFhHPnN7bL+WA7XQmJhzDufzJwE9nz1wfzuYZzo2fXE9n82Sy4ydfc3xwlOHsqefzM/wVcLp4pKYQ6rEQbfW1PHTvDRfcjlIpuEUk8GpqTh3IpeHtv18uP8FwNk9mbPxk0I/kxov/AeQZzo4XHucKQZ/JjZPJjROPVOYKlQpuEZGzRMM1RMNRZjnkXzG6gLGISMAouEVEAkbBLSISMApuEZGAUXCLiASMgltEJGAU3CIiAaPgFhEJGHMzXdrrQt7UrBs4OIuXtAAn5rwQf6vGNkN1trsa2wzV2e630+blzrnWUnYsS3DPlpltc85t9LqOSqrGNkN1trsa2wzV2e5KtVlDJSIiAaPgFhEJGL8E9wNeF+CBamwzVGe7q7HNUJ3trkibfTHGLSIipfNLj1tEREpUseA2szvM7HUze8PM/nSK52Nm9mDx+S1mtqJStZVTCe3+IzN71cxeMrOfmdlyL+qcSzO1+bT97jIzZ2bzYuZBKe02sw8WP+9dZvadStc410r4+V5mZk+Z2YvFn/E7vahzLpnZ18ysy8xemeZ5M7MvFv9NXjKzDXNehHOu7F9ACHgTWAVEgZ3A+rP2+QPgy8X1u4EHK1GbD9p9K5Aort8X9HaX0ubifnXAJuB5YKPXdVfos14DvAg0FR8v9LruCrT5AeC+4vp64IDXdc9Bu98JbABemeb5O4EnAAOuB7bMdQ2V6nFfC7zhnNvnnMsB3wPef9Y+7wf+qbj+feBdFvzbMc/YbufcU865TPHh80BHhWuca6V81gB/BfxPYLSSxZVRKe3+JHC/c64PwDnXVeEa51opbXZAfXG9AeisYH1l4ZzbBPSeZ5f3A99wBc8DjWbWNpc1VCq4lwCHTnt8uLhtyn2cc3lgAGiuSHXlU0q7T/dxCv9TB9mMbS7+6bjUOffjShZWZqV81muBtWb2nJk9b2Z3VKy68iilzX8JfMTMDgOPA/+hMqV5ara/97Ome076hJl9BNgI3Ox1LeVkZjXAF4Df97gUL4QpDJfcQuEvq01mdrlzrt/TqsrrQ8DXnXN/bWY3AN80s8uccxNeFxZklepxHwGWnva4o7htyn3MLEzhz6qeilRXPqW0GzP7N8CfA+9zzmUrVFu5zNTmOuAy4GkzO0BhDPCxeXCAspTP+jDwmHNuzDm3H9hDIciDqpQ2fxx4CMA5txmIU7iex3xW0u/921Gp4P4lsMbMVppZlMLBx8fO2ucx4KPF9d8GnnTFkf4Am7HdZnY18H8ohHbQxzxhhjY75waccy3OuRXOuRUUxvXf55zb5k25c6aUn/FHKfS2MbMWCkMn+ypZ5Bwrpc1vAe8CMLNLKAR3d0WrrLzHgN8rzi65Hhhwzh2d0+9QwSOxd1LoYbwJ/Hlx23+l8EsLhQ/0YeANYCuwyuujxxVq90+B48CO4tdjXtdc7jafte/TzINZJSV+1kZhmOhV4GXgbq9rrkCb1wPPUZhxsgO43eua56DN3wWOAmMU/or6OHAvcO9pn/P9xX+Tl8vx860zJ0VEAkZnToqIBIyCW0QkYBTcIiIBo+AWEQkYBbeISMAouEVEAkbBLSISMApuEZGA+f/yhepc9NFKKgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline  \n",
    "\n",
    "# Divergence problem\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from math import pi\n",
    "\n",
    "\n",
    "def f(x):\n",
    "    return 1/np.sqrt(x)/(np.exp(x)+1)\n",
    "\n",
    "x = np.linspace(0,1,100)\n",
    "plt.plot(x, f(x))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you try to do it using the mean value method, you will run into problems becuase the value of f(x) diverges when x is approaching to 0. This means that the estimated integral could go up and down without convergence.\n",
    "\n",
    "We can get around these problems by drawing our points x `non-uniformly`. This technique is called `importance sampling`. For any general functions $g(x)$ we can define a weighted average over the interval from $a$ to $b$ thus:\n",
    "\n",
    "$$ \\langle g \\rangle_w = \\frac{\\int_a^b w(x)g(x)dx}{\\int_a^b w(x)dx} $$\n",
    "\n",
    "where $w(x)$ is any function we choose. Now consider the general 1d integral:\n",
    "\n",
    "$$ I = \\int_a^b f(x)dx $$\n",
    "\n",
    "Setting g(x)=f(x)/w(x), we have\n",
    "$$ \n",
    "\\begin{aligned}\n",
    "\\bigg \\langle \\frac{f(x)}{w(x)} \\bigg \\rangle _w= & \\frac{\\int_a^b w(x)f(x)/w(x)dx}{\\int_a^b w(x)dx} \\\\\n",
    "                      = & \\frac{\\int_a^b f(x)dx}{\\int_a^b w(x)dx}\\\\\n",
    "                      = & \\frac{I}{\\int_a^b w(x)dx}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "So\n",
    "$$ I = \\bigg \\langle \\frac{f(x)}{w(x)} \\bigg \\rangle _w \\int_a^b w(x)dx$$\n",
    "\n",
    "This formula is similar to what we have for the mean value method. But it allows us to calculate the value\n",
    "of the integral from a weighted average, rather than a standard uniform average. \n",
    "\n",
    "Now we need to solve the $\\langle \\frac{f(x)}{w(x)} \\rangle _w$. Let's define a probability function\n",
    "\n",
    "$$ p(x) = \\frac{w(x)}{\\int_a^b w(x)dx} $$\n",
    "\n",
    "The average number of samples that fall in this interval is $Np(x)dx$, so for any function g(x)\n",
    "\n",
    "$$\\sum_{i=1}^Ng(x_i) \\approx \\int_a^b Np(x)g(x)dx$$\n",
    "\n",
    "Using this result, we can now write the general weighted average of the function g(x) as\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\langle g \\rangle_w = &\\frac{\\int_a^b w(x)g(x)dx}{\\int_a^b w(x)dx} \\\\ \n",
    "                    = &\\int_a^b p(x)g(x)dx \\\\\n",
    "              \\approx & \\frac{1}{N} \\sum_{i=1}^N g(x_i)    \n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "We then get \n",
    "\n",
    "$$ I \\approx  \\frac{1}{N} \\sum_{i=1}^{N} \\frac{f(x_i)}{w(x_i)} \\int_a^b w(x)dx$$\n",
    "\n",
    "In order to calculate it, we need to get the `tranformation` from uniform $x$ to weighted $x_i$.\n",
    "\n",
    "- Search for the weight founction $w(x)$\n",
    "- find the probability function $p(x)$ and then the cumulative distribution function $F(x)$\n",
    "- Inverse transform sampling to $F^{-1}(x)$\n",
    "- averge $f(x)/w(x)$ in the points generated from $F^{-1}(x)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12.4 Homework\n",
    "\n",
    "- 1. Calculate a value for the integral\n",
    "\n",
    "$$ I = \\int_0 ^1 \\frac{x^{-1/2}}{e^x + 1} dx $$\n",
    "\n",
    "using the importance sampling method with $w(x) = x^{-1/2}$, with 1,000,000 random points.\n",
    "\n",
    "- 2. Calculate a value for the integral\n",
    "\n",
    "$$ I = \\int_0 ^1 \\frac{x^{-1}}{e^x + 1} dx $$\n",
    "\n",
    "using the importance sampling method with $w(x) = 1/x$, with 1,000,000 random points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/qiangzhu/.pyenv/versions/3.7.0/Python.framework/Versions/3.7/lib/python3.7/site-packages/ipykernel_launcher.py:16: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  app.launch_new_instance()\n",
      "/Users/qiangzhu/.pyenv/versions/3.7.0/Python.framework/Versions/3.7/lib/python3.7/site-packages/ipykernel_launcher.py:17: RuntimeWarning: divide by zero encountered in true_divide\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8HPV9//HXdw/t6pYsybJs+QIMGNuAjUyAxAQwTQIkBgqhJPCDNMEk6S+kJUl/JZQ2KfzaX9NC0tIc4CYEaGghN4ZwJGATwInBMjaY0xhfyJYtWbZu7f39/TGrlXzIkqW9Rno/H495zOzs7M5nMH7v+Dvf+Y6x1iIiIu7nyXUBIiKSHgp0EZFxQoEuIjJOKNBFRMYJBbqIyDihQBcRGScU6CIi44QCXURknFCgi4iME75s7qy6utrOmjUrm7sUEXG99evX77PW1gy3XVYDfdasWTQ2NmZzlyIirmeM2TGS7dTkIiIyTijQRUTGCQW6iMg4kdU2dBGR0YpGozQ1NREKhXJdSsYEg0Hq6+vx+/2j+vywgW6MuQ/4ONBirZ2fXDcJeASYBWwHrrLWHhhVBSIiI9DU1ERpaSmzZs3CGJPrctLOWktbWxtNTU3Mnj17VN8xkiaX+4GPHbLuFuBZa+0c4NnkaxGRjAmFQlRVVY3LMAcwxlBVVTWmf4EMG+jW2ueB/YesvhR4ILn8AHDZqCsQERmh8Rrm/cZ6fKO9KFprrW1OLu8Baofa0BhzozGm0RjT2NraOrq9vfoINN43us+KiEwQY+7lYp2Hkg75YFJr7QprbYO1tqGmZtgbnY7sjV8q0EUk5+6++27mzp3LNddcw69//Wtuv/32o27/ta99jVWrVmWputH3ctlrjKmz1jYbY+qAlnQWdZiCEgh3Z3QXIiLD+f73v88zzzxDfX0955xzDitXrjzq9jfddBPLly/nggsuyEp9oz1DXwlcn1y+Hng0PeUMIVAKEQW6iOTOF77wBbZu3cpFF13Et771LQKBANXV1QBceumlPPjggwDce++9XHPNNQDMnDmTtrY29uzZk5UaR9Jt8X+A84BqY0wT8A3gn4GfGmM+B+wArspkkQRKINyV0V2IiHv8w2Nv8ObuzrR+5ylTy/jGJ+YN+f4999zDU089xerVq3nsscdYtGhR6r0VK1bwwQ9+kNmzZ3PXXXexdu3a1HuLFi1izZo1XHHFFWmt90iGDXRr7aeGeGtpmmsZWqAMYiGIR8E7ug73IiLp0tzczOBrgrW1tdx+++2cf/75/OpXv2LSpEmp9yZPnszu3buzUpc77hQtKHHm4S4omnT0bUVk3DvamXQ2FBYW0tHRcdC6TZs2UVVVdVh4h0IhCgsLs1KXO8ZyCSQDXe3oIpIH5s6dy5YtW1KvX375ZZ588kk2bNjAnXfeybZt21Lvbd68mfnz52elLpcEeqkzV08XEckD5557Lhs2bMBaSzgcZvny5dx3331MnTqVu+66i89+9rNYa4lGo2zZsoWGhoas1OWSJpf+QNeFURHJne3bt6eWL7zwQp599lkuvPBCXn311dT6ZcuWsWzZMgAef/xxrrzySny+7EStS87Q+5tcFOgikh9uvfVWent7j7pNLBbjq1/9apYqcssZuppcRCTP1NbWps7Eh/LJT34yS9U43HGGPriXi4iIHJE7Ar3/DF29XEREhuSuQNcZuojIkNwR6F4/eAMKdBGRo3BHoIMG6BKRvGat5YILLqCzc+gxZlpbW/nYxw59AFz6uCjQNUCXiOSvJ554gtNOO42ysrIht6mpqaGuro41a9ZkpAb3BHpBqbotikjO/Ou//it33303ADfffHNqjPNVq1ZxzTXX8NBDD3HppZcCsG7dOk499VRCoRA9PT3MmzeP119/HYDLLruMhx56KCM1uqMfOqjJRUQGPHkL7NmU3u+csgAu+uch316yZAl33XUXX/7yl2lsbCQcDhONRnnhhRc499xz+ad/+ifuvfdeABYvXsyyZcu47bbb6Ovr49prr02N59LQ0MBtt92W3tqTXBToJdC9N9dViMgEdcYZZ7B+/Xo6OzsJBAIsWrSIxsZGXnjhBe6++272799PaWlpavu///u/Z/HixQSDwdSZPWR2OF0XBXoptL2X6ypEJB8c5Uw6U/x+P7Nnz+b+++/nnHPO4dRTT2X16tVs2bKFuXPn4vP5SCQSeDxOS3ZbWxvd3d1Eo1FCoRDFxcVAZofTdVEbeomaXEQkp5YsWcKdd97Jueeey5IlS7jnnntYuHAhxhhOOukktm7dmtr285//PHfccQfXXHMNf/M3f5Nan8nhdN0T6IFS9XIRkZxasmQJzc3NnH322dTW1hIMBlmyZAkAl1xyCc899xwADz74IH6/n09/+tPccsstrFu3jlWrVgGwevVqLrnkkozU564ml2gvJOLg8ea6GhGZgJYuXUo0Gk293rx5c2r5hhtu4LrrrkvNr7vuOgC8Xi8vvfRSaruVK1fy6KOPZqQ+95yha4AuEcljdXV1LF++fNgbi77yla9QWVmZkRpcdIY+6DF0hRW5rUVE5Aiuuuqqo75fU1PDZZddlrH9u+cMXWOii0x41tpcl5BRYz0+9wS6HkMnMqEFg0Ha2trGbahba2lrayMYDI76O1zY5KJAF5mI6uvraWpqorW1NdelZEwwGKS+vn7Un3dRoKvJRWQi67+xR4bmoiYX9XIRETka9wR6IDkkpe4WFRE5IhcFus7QRUSOxj2B7guAx69AFxEZgnsCHTQmuojIUbgs0PUYOhGRoYwp0I0xNxtj3jDGvG6M+R9jzOh7xI+EHkMnIjKkUQe6MWYa8GWgwVo7H/ACV6ersCMKlOrGIhGRIYy1ycUHFBpjfEARkJnnKvVTk4uIyJBGHejW2l3AncBOoBnosNb+Nl2FHVFATS4iIkMZS5NLJXApMBuYChQbY649wnY3GmMajTGNYx6DQY+hExEZ0liaXC4EtllrW621UeCXwDmHbmStXWGtbbDWNtTU1Ixhd+gxdCIiRzGWQN8JnGWMKTLGGGAp8FZ6yhpCfz/0RCKjuxERcaOxtKG/BPwceAXYlPyuFWmq68j6B+iK9mR0NyIibjSm4XOttd8AvpGmWoY3eDyX/uF0RUQEcN2doskRF9XTRUTkMO4KdI2JLiIyJHcFen8zi+4WFRE5jMsCvf8MXU0uIiKHclegq8lFRGRI7gp0PYZORGRILgt0naGLiAzFXYHuC4LxKtBFRI7AXYFujB5DJyIyBHcFOmgIXRGRIbg00DtzXYWISN5xX6BrTHQRkSNyX6DrMXQiIkfkwkBXG7qIyJG4NNB1hi4icij3BXpJLXTvhUQ815WIiOQV9wV6eT3YOHTtyXUlIiJ5xYWBPt2ZdzTltg4RkTzj4kB/P7d1iIjkGRcG+jRnrkAXETmI+wI9UArBCjW5iIgcwn2BDk6ziwJdROQgLg30egW6iMghXBzoakMXERnMvYEe6oCQRl0UEennzkCvSHZd7NyV2zpERPKIOwO9vy96u5pdRET6uTTQ65252tFFRFLcGeglteDxqaeLiMgg7gx0jxfKpirQRUQGcWegg24uEhE5hIsDXTcXiYgMNqZAN8ZUGGN+box52xjzljHm7HQVNqzyeqfboh50ISICjP0M/d+Bp6y1JwOnAW+NvaQRKp+efNBFc9Z2KSKSz0Yd6MaYcuBc4EcA1tqItbY9XYUNSw+6EBE5yFjO0GcDrcCPjTEbjDE/NMYUp6mu4aX6oivQRURgbIHuAxYBP7DWLgR6gFsO3cgYc6MxptEY09ja2jqG3R1CD7oQETnIWAK9CWiy1r6UfP1znIA/iLV2hbW2wVrbUFNTM4bdHUIPuhAROcioA91auwd43xhzUnLVUuDNtFQ1UuqLLiKS4hvj528CHjLGFABbgT8fe0nHQOOii4ikjCnQrbUbgYY01XLsKqbDjj/kbPciIvnEvXeKAkw6DsId0LUn15WIiOScuwO97jRn3vxqbusQEckD7g70KQsAo0AXEcHtgR4ohaoTYPfGXFciIpJz7g50gKmnQ7MCXUTE/YFed5oz6mJ3Gu9CFRFxoXEQ6Kc7c7Wji8gENw4C/VRnrmYXEZng3B/owXKnP7oCXUQmOPcHOjjNLrvV5CIiE9s4CfTToGMn9O7PdSUiIjkzPgJ9av+FUTW7iMjENT4CXUMAiIiMk0AvrISKmbpjVEQmtPER6KA7RkVkwnNFoK9+p4XfvNZ89I3qTocD23VhVEQmLFcE+k/+uIPvrt5y9I1mftCZb30u4/WIiOQjVwR6bXmQls7Q0Teqb3Da0t/9XXaKEhHJM+4I9NIgbT0RwrH40Bt5vHD8UtjyO0gksleciEiecEWgTykPANDSGT76hid+FHpadXFURCYkVwT65LIgAC1dwzS7HL8UMGp2EZEJyRWBPiUZ6Hs6hjlDL65y2tLf/W0WqhIRyS+uCPTaZKDvHe7CKMCcj8Cu9dCzL8NViYjkF1cEemWRnwKfZ4SB/ieAhS3PZrwuEZF84opAN8ZQWxYYWaBPOQ2KJ8O7T2e+MBGRPOKKQAen6+KekQS6x+OcpW95FuKxzBcmIpIn3BPo5cHhuy32m/MRCLXDjjWZLUpEJI+4J9CTZ+jW2uE3PvGjzqPpNvwk84WJiOQJ1wT6lPIAvZE43eERNKP4C+HUP4M3H9VgXSIyYbgm0I+p6yLAousgHoZNP8tgVSIi+cN1gT7szUX9piyAqQth/QMwkmYaERGXc12gj/gMHWDR9dDyBux6JUNViYjkjzEHujHGa4zZYIx5PB0FDaW2zBmga0RdF/vNvwL8RfDK/ZkpSkQkj6TjDP0vgbfS8D1HVVTgozToG35c9MGCZTD/T2HTLyDclbniRETywJgC3RhTD1wC/DA95RzdlLIR3lw02KLPQLQHNjyUkZpERPLFWM/Q/w34P0BWnihRWxZk70hvLuo3fTHM/BC8+G2I9mWmMBGRPDDqQDfGfBxosdauH2a7G40xjcaYxtbW1tHuDugP9GM8Qwc4/+vQvRcafzym/YuI5LOxnKF/EFhmjNkOPAxcYIw57NZMa+0Ka22DtbahpqZmDLtzbi5q6QqTSBxjN8RZH4JZS+DF70Ckd0w1iIjkq1EHurX269baemvtLOBqYJW19tq0VXYEtWVB4gnLvp5jbHYBOP9W6GmBxh+lvzARkTzgmn7oMNAXfcSDdA028xw47jx48d8g0pPWukRE8kFaAt1a+5y19uPp+K6jGbhbdBTt6ADn3Qq9++AP301jVSIi+cFVZ+ipZ4uO5sIowIwPwLw/hRfugn1b0liZiEjuuSrQq0sK8BiO7eaiQ33s/4EvCL+5WWO8iMi44qpA93k9VJcERn+GDlA6BS78Bmx7Hl59OH3FiYjkmKsCHUZ5c9GhzvhzqD8Tnr4VetrSU5iISI65NNDHcIYOznNHP/FvEO5U04uIjBuuC/T6ykJ27u8lFh/jaAO18+CC25ynGq3LylA0IiIZ5bpAXzijgt5InLf3pGH0xHP+EuZ8FJ76usZMFxHXc12gL541CYB129PwrFCPBy6/B0pq4WfXQ9+BsX+niEiOuC7Qp1YUMq2ikMYdaQrfoknwyfuhczf84gaIR9PzvSIiWea6QAc4Y2Yljdv3Y9N1MXP6YrjkLtjyDDz2l7pIKiKu5MpAXzyrkr2dYZoOpHF88zM+Ax++BTY+BKv+b/q+V0QkS3y5LmA0zpjptKM37tjP9ElF6fvi826Brt3wwp1Ou/oHbkzfd4uIZJgrz9BPmlJKacDHuu1pvohpDFzyHTjpYnjyr2HtD9L7/SIiGeTKQPd6DItmVrI+3YEO4PXBJx+AuZ+Ap26BF76d/n2IiGSAKwMdoGFmJe/s7aKjNwO9UnwFcOX9MP9KePYf4Nk7dKFURPKeewM92R/9lZ0Z6jvu9cGfroCF/8tpU//FDXrItIjkNdcG+unTK/B5THpuMBqKxwvL/gOWfgNe/znc/3Ho2pu5/YmIjIFrA72wwMu8aeU0ZqIdfTBjYMlX4M9+Ai1vwn+eDzv+mNl9ioiMgmsDHWDxzEo2NrXTHY5lfmdzPwGffRq8BXD/JfD8nZAY4wBhIiJp5OpAv2hBHZFYgsdf3Z2dHdadCp9/HuZdBqvugJ9cDh1N2dm3iMgwXB3oi2ZUMGdyCY80vp+9nQbL4IofwSfuhvdfhu+dBY0/Vi8YEck5Vwe6MYY/WzydDTvb2bw3DcPpjnzHcMb18MU/wNTT4fG/ggcv1YOnRSSnXB3oAJcvnIbfa3hkXRbP0vtNmg3XrYSPfwd2b4DvnwXPfBPC3dmvRUQmPNcHelVJgD85pZZfbdhFOBbPfgEeDzR8Fr7UCAs+CS9+B767GDb8BBI5qEdEJizXBzrAVQ3T2d8T4Zk3W3JXRGktXP4D+NzvoHQKPPq/4QfnwNu/Ufu6iGTFuAj0JXNqmFoezO7F0aFMPxOWr3LGg0nE4OFPO33XFewikmHjItC9HsOVDdN54d1WtrRk8eLoUIxxujb+xVrnTtO+A06w3/MheO1neiqSiGTEuAh0gOvPnklJwMc//uatXJcywOuHRdfBl9bD5fdCPAK/vAH+/XRYczf0tee6QhEZR8ZNoFeVBLjpghNY/U4rz29uzXU5B/P64LSr4S9egk894vSO+d3fwV0nw6NfcnrIiIiMkUnbczlHoKGhwTY2Nmbs+8OxOB/5zvMEfB6e+PISfN48/r1qfg3W/RA2/QyivVB3Oiy8FhZcCYWVua5ORPKIMWa9tbZhuO3yOPGOXcDn5esXnczmvd08nIt+6cei7lRYdjd89W246F+cLo5PfA3uPAl+er1zETUWznWVIuIi4+oMHcBay9Ur1vJuSzerv3oe5UX+jO4vrZpfhQ0POUP19rZBsMIZFGzeZTD7w06bvIhMOCM9Qx91oBtjpgMPArWABVZYa//9aJ/JRqADvLG7g8u+t4YPnziZ/7zuDIwxGd9nWsWjsPU5eO2n8M4TEOl2wv2ki+Hki+H4C6CgONdVikiWZCPQ64A6a+0rxphSYD1wmbX2zaE+k61AB/jxmm38w2Nv8vWLTubzHz4+K/vMiGgI3lsFbz4K7zwJ4Q7wBuC4D8OcjzhT5cxcVykiGTTSQPeNdgfW2magObncZYx5C5gGDBno2fSZc2bRuP0A//L0O5w+vYIPHFeV65JGxx90zspPvtg5c9/xByfY33kC3v2ts031ic5Z+3Hnw6wPQaAktzWLSE6kpQ3dGDMLeB6Yb63tHGq7bJ6hA3SFoiz77hp6wjEev+lDTC4LZm3fGWcttG2Bd38HW37nBH0sBB4fTDsDZi2B2Uug/kwoKMp1tSIyBhlvchm0oxLg98A/Wmt/eYT3bwRuBJgxY8YZO3bsGNP+jtVbzZ1c8YM/MLWikIdvPIvqkkBW95810RC8v9Zpe9/2gtO33cbB43eG+J15Dsw42wn4Ypf+a0VkgspKoBtj/MDjwNPW2m8Pt322z9D7rd3axmd+/DKzqor57+VnMam4IOs1ZF2oE3auhR1rYOcfYdcrkEgOOVB1AtQvds7k6xugdr560IjksWxcFDXAA8B+a+1fjeQzuQp0gDVb9vHZ+9dxfE0JD93wASonQqgPFu1zztrff8l50lLTOuhJ3lHrDcCUBTB1oXM2X3ca1JyskBfJE9kI9A8BLwCbgP6nJd9qrX1iqM/kMtABfr+5leUPNDK1IsgPr2/ghMmlOasl56yF9p2wq9E5e9+9EZo3Ol0kwXkY9uRTnKCfssA5i689RXexiuRA1trQj0WuAx2gcft+vvCT9YSjCf7j0ws576TJOa0nryTisH+rc4NT80ZneIK9rzs3OfUrnQq182DyyVAz15lXn6SeNSIZpEA/iqYDvSx/cD3v7Onk5gtP5IvnHZ/f477kkrXQtccJ9r1vOFPLm7BvszN6ZL+yeqg5yelCWT3HmarmOA/7cNuNXSJ5RoE+jJ5wjFt+uYnHXt3NohkVfPuq05lVrbsvRywegwPboOUt2PcOtG6G1redrpTR3oHt/MVQdbwzTToOJvXPZ0NJrcJeZAQU6CNgrWXlq7v5u1+/TjRu+euPnsR1Z8/U2fpYJBLQtds5g297Lzm96zTlHNjhdKXs5y+CyllQMdOZV86EihnO64oZECzL1VGI5BUF+jFo7ujjll9s4vebWzmxtoRvfmIe55xQneuyxp941LkQu3+bc3a/fxsc2D4wRXsO3j5YARXToXxGcl7vTGXJeclk8HhzcCAi2aVAP0bWWn775l7uePxNmg708Sen1HLzhSdyylSdJWaFtc7F1/Ydzpl8+07oeN+Zt++EjqaBHjj9PD7nIm3ZoKm0btDyFOe1b5zeTCYThgJ9lELROD96cRv3/P49ukIxLl4whS8vncPJUxTsOWUthDqckO/YBZ1Nyflu6NyVnJoh1nf4ZwsrnWAvqXVCfvA8NU2GQKna9CUvKdDHqKM3yo9e3Mp9a7bTHY6xZE41Nyw5jnPnVLtvON6JwloItTsh39XsBHxXs9NLp3tvcnmvs5w4woO6fUEongwlNQfPi2uguHpgXlQNRVXOowVFskCBnibtvREeemknD/xhOy1dYU6YXMLVi6dzxaL6iXe36XiRSEDfAejeA90tyWkv9LQMvO5pTU77Dr6QO1hh5UC4FyfnqWmSMy+clFyeBIFy8OiCuxw7BXqaRWIJHnt1N/+1dgcb32+nwOvhI/NquXzhNJbMqaHAp7+o41J/+Pe0OoHfs89p6+/ZB737Bl4PnhKxI3+X8Tg/AoWTkvNDpqJJzoXgwkoorEguV0CwXMMwTHAK9Ax6e08nD7/8Pr/euIv23igVRX4uml/HxQumcNZxVfjV7XHishbCnclw35+c2qAvudy33/mB6N3vNA/1HYDeAxDpOvr3FpQ4wT445IP983Kni2f/cqBs4HUg+Z5+EFxNgZ4FkViCF7e08ujG3fz2jb30ReOUBX0snVvLBSdP5tw5Ne56pqnkTjzqXPRNBX0y7EPtzvq+9oH1oXZnNM1Qh7McHvIRBAN8hU6wB0qTU3I5WH7IupKB9wKlzg9J6v1S9RjKEQV6loWicV54dx9Pvb6HZ9/eS3tvFK/HsGhGBUvm1PDBE6o4tb5CZ++Sfok4hLucgA/3B32HE/rhzuT80NddyWnQMiPIAo/fCf2C0uS8JDkvHrSu2Fl/0Hv982Ln7uGCQZPuJRiWAj2H4gnLxvfbee6dFp57p5XXd3dgLZQEfCyeVcmZs6s4c/YkFkwrV9u75IdEwunnH+keFPaDpv71kW4Idx/5daRn4PWRehENxRccFPRFyeXB86KB91LzwsOX/YWDtk+u9wXHRVdUBXoeOdAT4Y9b23hxyz5e2trGe63OHZEBn4dT68tZNKOShTMqOG16BVPKguoWKe4XCx8c8NHew5cjPRDpde4QTq3vcab+5WjvwDaR3mP7oQDADIR7aj54KnJCP/V+8JB1hU5zlT94hHlw4EfDX+Q0R2Xo764CPY/t6w6zbtt+1u84wPqdB3hjVyeRuDOkfE1pgFOnlTNvWjnzp5Yxb1o5U8sV8iKAc62hP+ijfYOWk69T6/oOX3+05VhoYN3gUUSPiXHC3RcYCHpfcCD8L7/XGa9oNN88wkDXnRE5UF0S4KIFdVy0oA5w2t/fau7ktaYOXm1q5/VdHax+p4VE8re2LOjj5Loy5k4pZU5tKSdNKeXEyaW64CoTj9fv9PIprMjcPhLxZMAnAz8V9iHnTuQjzkMDn4mFj/x+Fq4VKNDzQNDvZeGMShbOGHgaUF8kzlt7OnljdydvN3fy9p4ufvHKLrrDA32cq0sCnDC5mBMml3BcdQnH1RRzfE0JUysK8Xp0Ri8yKh7vwAVbl1Gg56nCAi+LZlSyaFDIW2vZ3RFi894u3t3bxZaWbt5t6ebRjbvpCg0Evd9rmD6piNlVxcyoKmLmpCJmVBUxY1IR9ZVFBP3qVSAyHinQXcQYw7SKQqZVFHL+oEfnWWtp64mwtbWHra3dbG/rZUdbD9v29fDHrW30Rg6+db2mNEB9pfM99ZVFTKssZFpFkKkVhdSVF1IW9KnNXsSFFOjjgDGG6pIA1SUBzpw96aD3+sN+R1sv7+9PTgd6aTrQx6ZdHTz9xh6i8YMvjBcVeKkrD1JXXkhtWZC68iC15UFqSwPUlgWpLQtSXVKgB4GI5BkF+jg3OOzPmFl52PuJhGVfd5hd7X3sbg+xu72P5o4QzR3O/L339tHSFSaesId8L1QVF1BTGmRyaYCa5OTsq4CakgDVpQGqiguoLCrAozZ9kYxToE9wHo9hclmQyWVBFs448jbxZOi3dIbZ2xliT2eIlq4wrV0hWjrDtHaH2by3i9auMLHE4d1gvR5DZVEBVcUFTCouYFLJwHJVcQGVxQVMKiqgoshZV1HkVzu/yCgo0GVYXo9JNbUsoHzI7ay1dPRFae1yQr6tO0Jbd5h93RHaepzltp4Ib+3upK0nQkff0DeJFPq9VBb5qShyAr6yqIDyIj8VhX7KC/1UFPkpLyygvNBPWaGP8uT6koDa/2XiUqBL2hhjkgFcwJza0mG3j8UTHOiNcqA3wv4eZzrQG6G9N8qBnggHeqN09Dnzt/d00tEXpb03esR/BfTzGChLhntp0EdZ0O9Mhc5yadBZ70x+yoI+SpLLJQFnfcDn0Y+CuJICXXLG5/Wk2t5HylpLTySeDHfnLL+zL0pHcursiznz0MD61q4wXaEYnaHoYT1+jsTvNZQEnKAvCfgpDfgoDngpCfopCXgpLuh/z0dRQfK9gI/igI/i5OvigI+iAmdbXT+QbFGgi6sYkwzbgI9pFYXH/PlYPEF3OJYK+K5QLDlF6QnH6Eq+1x2KJbeL0h2Osa87wva2XrrDMXrCsRH9MPQr9HspDngpKnBCvqjACXxnvY/CAi9Ffmd94aBt+rcPpt7zUjhoOejz6sdCDqJAlwnF5/WkmoXGIp6w9EZi9ITjqZDvCTs/Ar2ROD2RgeDvjcRT7/e/7g7HaO0KJ18PrD9WAZ8nFfT986B/YB70e1LvBQ9Z178c9DnLAb/HWedLrh+0TcDn1d3HLqBAFxkFr8ck2+PTN55OImEJxZxf3y54AAAGH0lEQVRg74vE6YvGU4E/+HVfJE4o6rzuX98XidMbjRNOru+JxGjriaRe90Wdz4SiiVHX5/caAj4vAZ8n9QPQ/zq1zuch4B9YF/D1bzdoW7+HAq+znTP3EPB6KEhuU5D8bMGgKeBzPqNrG0enQBfJEx6PSTazZO6vpbWWcCyRCveBoHdCPxxLEE6+178+FEsQjiYIx5z14Vg89R2peTRBe1+UcDROJLkuEnc+F4rFD7t5bbQKvIOCfrjl5Nzv9Rz0I+H3mtS6g7bxevD7PBR4TWpd/+S8ZwaWvcnv8Q28zod/wSjQRSYQY0yqKSWbEglLJJ4M+ljC+eFI/jCEYwkiyal//eDXkVjC+XEYtF0kHh+0fPC2PZEY7X2Dto0liMQtkeQPSySeOOxGuXQwhkGBb/ANWvZ7Pfzo+sXMqCpK+34HU6CLSMZ5PIagJ/s/JEOJJyzR+MCPQSxuUz8U0dRygmg8kXodG7x9/+eT28XiA9+XWo45n40lnO8J+DM/VIYCXUQmHK/H4M2jH5h0GdNPhjHmY8aYd4wxW4wxt6SrKBEROXajDnRjjBf4HnARcArwKWPMKekqTEREjs1YztDPBLZYa7daayPAw8Cl6SlLRESO1VgCfRrw/qDXTcl1BzHG3GiMaTTGNLa2to5hdyIicjQZv+xqrV1hrW2w1jbU1NRkenciIhPWWAJ9FzB90Ov65DoREcmBsQT6OmCOMWa2MaYAuBpYmZ6yRETkWI26H7q1NmaM+RLwNOAF7rPWvpG2ykRE5JgYa9N/C+yQOzOmFdhxDB+pBvZlqJx8NRGPGSbmcU/EY4aJedxjPeaZ1tphL0JmNdCPlTGm0VrbkOs6smkiHjNMzOOeiMcME/O4s3XMmR9cQEREskKBLiIyTuR7oK/IdQE5MBGPGSbmcU/EY4aJedxZOea8bkMXEZGRy/czdBERGaG8CPThhuE1xgSMMY8k33/JGDMr+1Wm1wiO+SvGmDeNMa8ZY541xszMRZ3pNtIhl40xVxhjrDHG9b0hRnLMxpirkn/ebxhj/jvbNWbCCP4fn2GMWW2M2ZD8//ziXNSZTsaY+4wxLcaY14d43xhj7k7+N3nNGLMorQVYa3M64dyU9B5wHFAAvAqccsg2fwHck1y+Gngk13Vn4ZjPB4qSy190+zGP9LiT25UCzwNrgYZc152FP+s5wAagMvl6cq7rztJxrwC+mFw+Bdie67rTcNznAouA14d4/2LgScAAZwEvpXP/+XCGPpJheC8FHkgu/xxYatz9+O9hj9lau9pa25t8uRZnrBy3G+mQy3cA3wJC2SwuQ0ZyzMuB71lrDwBYa1uyXGMmjOS4LVCWXC4Hdmexvoyw1j4P7D/KJpcCD1rHWqDCGFOXrv3nQ6CPZBje1DbW2hjQAVRlpbrMGNHQw4N8DudX3e2GPe7kP0GnW2t/k83CMmgkf9YnAicaY9YYY9YaYz6WteoyZyTH/U3gWmNME/AEcFN2SsupY/27f0z0TNE8Z4y5FmgAPpzrWjLNGOMBvg18JselZJsPp9nlPJx/iT1vjFlgrW3PaVWZ9yngfmvtXcaYs4H/MsbMt9Ymcl2YW+XDGfpIhuFNbWOM8eH886wtK9VlxoiGHjbGXAj8LbDMWhvOUm2ZNNxxlwLzgeeMMdtx2hhXuvzC6Ej+rJuAldbaqLV2G7AZJ+DdbCTH/TngpwDW2j8CQZwxT8azjA47ng+BPpJheFcC1yeXrwRW2eQVBpca9piNMQuBe3HCfDy0qcIwx22t7bDWVltrZ1lrZ+FcO1hmrW3MTblpMZL/v3+Nc3aOMaYapwlmazaLzICRHPdOYCmAMWYuTqCP98earQSuS/Z2OQvosNY2p+3bc31VeNCV3804V8X/Nrnudpy/zOD8Qf8M2AK8DByX65qzcMzPAHuBjclpZa5rzsZxH7Ltc7i8l8sI/6wNTlPTm8Am4Opc15yl4z4FWIPTA2Yj8JFc15yGY/4foBmI4vzL63PAF4AvDPqz/l7yv8mmdP//rTtFRUTGiXxochERkTRQoIuIjBMKdBGRcUKBLiIyTijQRUTGCQW6iMg4oUAXERknFOgiIuPE/wcv7grvY/ftoQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline  \n",
    "\n",
    "# Divergence problem\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from math import pi\n",
    "\n",
    "f = lambda x: 1/np.sqrt(x)/(np.exp(x)+1)\n",
    "w = lambda x: 1/np.sqrt(x)\n",
    "p = lambda x: x*x\n",
    "\n",
    "x = np.linspace(0,1,100)\n",
    "plt.plot(x, f(x), label='f(x)')\n",
    "plt.plot(x, w(x), label='w(x)')\n",
    "\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "from scipy:  0.8389329600133838\n"
     ]
    }
   ],
   "source": [
    "from scipy.integrate import quad\n",
    "print('from scipy: ', quad(f, 0, 1)[0])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8416504978875713\n"
     ]
    }
   ],
   "source": [
    "from random import random \n",
    "f = lambda x: 1/np.sqrt(x)/(np.exp(x)+1)\n",
    "w = lambda x: 1/np.sqrt(x)\n",
    "p = lambda x: x**2\n",
    "\n",
    "def IMC(N):\n",
    "    const = quad(w, 0, 1)[0]\n",
    "    I = 0\n",
    "    for i in range(N):\n",
    "        x = random()\n",
    "        y = p(x)\n",
    "        I += f(y)/w(y)\n",
    "    return I/N*const\n",
    "\n",
    "print(IMC(1000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0867888907301615\n"
     ]
    }
   ],
   "source": [
    "def MC(N):\n",
    "    I = 0\n",
    "    for i in range(N):\n",
    "        x = random()\n",
    "        I += f(x)\n",
    "    return I/N\n",
    "print(MC(100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x110ac9780>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsnXeYZEXZ9n81PaEn552Ns3nZvGwClhwWWIIgggRBCQqGF/VVUVFQFF9Q+VARFBEVkSBBJCdZgWVhYSNsjrM57+Q80+l8f9Sp7tOnT4eZ6ekZZuu+rrl6+nT36eoT7rrrfp56ShiGgYaGhobG0YG0vm6AhoaGhkbqoElfQ0ND4yiCJn0NDQ2Nowia9DU0NDSOImjS19DQ0DiKoElfQ0ND4yiCJn0NDQ2NowhxSV8I8YgQ4ogQYn2U14UQ4n4hRJUQYq0QYpbltWuFENvMv2uT2XANDQ0Nja4jEaX/KLAgxuvnAePNv5uAPwEIIUqAO4DjgeOAO4QQxT1prIaGhoZGz5Ae7w2GYSwWQoyK8ZaLgccMObV3qRCiSAgxBDgdWGgYRh2AEGIhsvN4Ktb3lZWVGaNGxfo6DQ0NDQ07Vq1aVWMYRnm898Ul/QQwDNhreb7P3BZte0yMGjWKlStXJqFZGhoaGkcPhBC7E3lfvwjkCiFuEkKsFEKsrK6u7uvmaGhoaAxYJIP09wMjLM+Hm9uibY+AYRgPG4YxxzCMOeXlcUcnGhoaGhrdRDJI/2XgS2YWzwlAo2EYB4H/AOcIIYrNAO455jYNDQ0NjT5CXE9fCPEUMihbJoTYh8zIyQAwDOMh4HXgfKAKaAOuN1+rE0L8Alhh7upOFdTV0NDQSAW8Xi/79u2jo6Ojr5uSNLjdboYPH05GRka3Pi/6Wz39OXPmGDqQq6GhkQzs3LmT/Px8SktLEUL0dXN6DMMwqK2tpbm5mdGjR4e9JoRYZRjGnHj76BeBXA0NDY3eQEdHx4AhfAAhBKWlpT0auWjS19DQGNAYKISv0NPfo0lfQ0PDGVvfgoY9fd0KjSRDk76GhoYznrselv25r1vxqUdeXh4Au3btQgjB7bffHnytpqaGjIwMbr755uC2xx57jKlTpzJt2jRmzpzJvffem9T2aNLX0NCIRCAAnhbwDZysl/6A0aNH89prrwWf/+tf/2LKlCnB52+88Qb33Xcfb731FuvWrWPp0qUUFhYmtQ2a9DU0NCLhazcfO/u2HQMMOTk5TJo0KVhq5plnnuHyyy8Pvv7LX/6Se++9l6FDhwKQlZXFjTfemNQ2JKP2joaGxkCDp00++j19244k4uevbGDjgaak7nPy0ALu+MyU+G+04Morr+Tpp5+moqICl8vF0KFDOXDgAADr169n9uzZSW2jHZr0NTQ0IuEdeKTfX7BgwQJ+8pOfUFFRwRVXXJHy79ekr6GhEQmvsncGDul3VZH3FjIzM5k9eza/+c1v2LhxIy+//HLwtSlTprBq1SrOPPPMXvt+7elraGhEIqj0taffG/je977Hr3/9a0pKSsK2/+hHP+L73/8+hw4dAsDj8fDXv/41qd+tlb6GhkYktL3Tq5gyZUpY1o7C+eefz+HDh5k/fz6GYSCE4IYbbkjqd2vS19DQiMQAtHf6Ci0tLYBcIGr9+silxq+77jquu+664PPrr7+e66+/vtfao+0dDQ2NSGh7Z8BCk76GhkYkgimb3r5th0bSoUlfQ0MjEkrp68lZAw6a9DU0NCKhPH1t7ww4aNLX0NCIhFfbOwMVmvQ1NDQioe2dAQtN+hoaGpEI2js6ZbOnEEJwzTXXBJ/7fD7Ky8u58MILg9veeOMN5syZw+TJk5k5cybf+973eq09mvQ1NDQi4WmVj5r0e4zc3FzWr19Pe7vsSBcuXMiwYcOCr69fv56bb76ZJ554go0bN7Jy5UrGjRvXa+3RpK+hoREJq9I3jL5tywDA+eefH6yj/9RTT3HVVVcFX7vnnnu47bbbmDhxIgAul4uvf/3rvdYWPSNXQ0MjEor0QRJ/elbftSVZeONWOLQuufscPA3O+1Xct1155ZXceeedXHjhhaxdu5YbbriB999/H5BKvzftHDu00tfQ0IiEtzX0v7Z4eozp06eza9cunnrqKc4///w+bYtW+hoaGpGwKn2fBwaA0E9EkfcmLrroIm655RYWLVpEbW1tcLsqpzxjxoyUtEMrfQ0NjUiolE3QE7SShBtuuIE77riDadOmhW3//ve/z913383WrVsBCAQCPPTQQ73WDq30NTQ0ImH39DV6jOHDh/Otb30rYvv06dO57777uOqqq2hra0MIEZbOmWxo0tfQ0IiEpw3Ss+UC6bq8co+gSitbcfrpp3P66acHn1944YW9SvRWaHtHQ0MjEt42yC6W/2t7Z0BBk76GhkYkvO2QXST/1/bOgIImfQ0NjXAE/FLduwvl80+5vWMMsMllPf09mvQ1NDTCoTJ33Erpf3rtHbfbTW1t7YAhfsMwqK2txe12d3sfOpCroaERDpW5E7R3Pr3llYcPH86+ffuorq7u66YkDW63m+HDh3f78wmRvhBiAfB7wAX81TCMX9leHwk8ApQDdcA1hmHsM1/zA2ru8x7DMC7qdms1NDR6H6rYmlL6n+LyyhkZGYwePbqvm9GvEJf0hRAu4I/A2cA+YIUQ4mXDMDZa3nYv8JhhGP8QQpwJ/BL4ovlau2EYxya53RoaGr2FCKX/6SV9jUgk4ukfB1QZhrHDMAwP8DRwse09k4F3zP/fdXhdQ0Pj0wJF+u5Pv72jEYlESH8YsNfyfJ+5zYo1wOfM/y8B8oUQpeZztxBipRBiqRDisz1qrYaGRu9DFVvL/vTbOxqRSFb2zi3AaUKIT4DTgP2A33xtpGEYc4AvAPcJIcbaPyyEuMnsGFYOpICLhsanEhFK/9OdsqkRjkRIfz8wwvJ8uLktCMMwDhiG8TnDMGYCt5nbGszH/ebjDmARMNP+BYZhPGwYxhzDMOaUl5d353doaGgkCyplU0/OGpBIhPRXAOOFEKOFEJnAlcDL1jcIIcqEEGpfP0Jm8iCEKBZCZKn3ACcB1gCwhoZGf4Nd6Wt7Z0AhLukbhuEDbgb+A2wCnjUMY4MQ4k4hhEq/PB3YIoTYClQAd5nbJwErhRBrkAHeX9myfjQ0NPobPDZPXyv9AYWE8vQNw3gdeN227aeW/58DnnP43IfANPt2DQ2Nfgyl9DNzIS1dk/4Agy7DoKGhEQ5F+hk54MrU9s4AgyZ9DQ2NcHhbwZUFaS5J+lrpDyho0tfQ0AiHtx0ysuX/6Vma9HsTb/0E3v9tSr9Sk76GhkY4vG3S2gHT3tGk32uoeht2LErpV2rS1+jfaKvr6xYcffC0QaaF9HXtnd6Drx18HSn9Sk36Gv0XhzfAPWPgyKa+bsnRBW3vpA7ejtBkuBRBk75G/0XjfsCApgN93ZKjC9reSR187aFsqRRBk75G/4Ua9mqlmVrYSV/bO70Hb4cmfQ2NIBTp6zzx1MLbHiL99CxdWrm3YBim0tf2joaGhFJAmvRTC29byNPXk7N6D+q4aqWvoWFC3RTaXkgtdPZOauBToqYDAoGUfa0mfY3+C59W+n2CMHsnU9s7vQWvJVUzhWmbmvQ1+i+CSl8HclOKMHsnS3e6vQWfxdZJocWjSV+j/8JrGf5qpAZ+LwS8kJErn+vaO70Hq9JPYTBXk75G/4VSmDpPPHVQ5BOcnKVJv9eglb6Ghg3qptCBxNQhWFbZau9o0u8VaKWvoWGDVvqpR1Dpq+ydDN3p9ha00tfQsEF7+qmHxyT9TMvkLF+nnEikkVxopa+hYYPO0089rKtmgbR3MCDg77MmDVhYxYxW+hoaWPL0tb2TMtgDua4M+ag73uTDp/P0NTTCoYa/mnBSB7unn54lH3WufvJhVffa3tHQQBdc6wtEBHIz5aOelZt8aHtHQ8MGTfqpR0TKpiJ9fQ6SDh3I1dCwQdfTTz2C2TvmjNygvaPPQdLhawdhUnAKlX56yr5JQ6Or8Gqln3JEBHKV0tekn3R4OyA9GzBSqvQ16Wv0X2h7J/VQijNd2zu9Dl87ZLjl/1rpa2hgsXc04aQM3lZJ+Gmm7ZBukr62d5IPpfSF0KSvoSGXktNKP+XwtoesHTAnZ6Htnd6AUvoiTWfvaGjg94JhriakCSd1sC6gAtreiYaGvT3fh1L66W5N+hoafTVb8aiH17JUImh7xwmH1sF9U+HAJz3bj1L6GTk6ZVPDAftWwqH1fd2K1EERfXq2JpxUwtOm7Z14aNwf/thdeDukys/I1kpfw4KAH965C/46H56/sa9bkzoo0ncXamshlfC2RbF3NOkH0dksHz0tPduPz4yfZOT0P9IXQiwQQmwRQlQJIW51eH2kEOJtIcRaIcQiIcRwy2vXCiG2mX/XJrPxAx4tR+Dxz8Lie6CoEo5shNbavm5VaqBy9N0FEPDpKo+pgj2QG7R3dMcbRGej+djcs/2EKf1+ZO8IIVzAH4HzgMnAVUKIyba33Qs8ZhjGdOBO4JfmZ0uAO4DjgeOAO4QQxclr/gDH8zfB3hVw8R/hc3+R2/Z81LdtShVUhU13oflck05KEKH0tb0TgaQr/f5n7xwHVBmGscMwDA/wNHCx7T2TgXfM/9+1vH4usNAwjDrDMOqBhcCCnjf7KEHDbph4Acy8BobOlKpg95K+blVqoEhekX5fWjwH18CL3zg6RhsRpK9KK2vSD0KRfmcPST+o9PtfIHcYYM1P2mdus2IN8Dnz/0uAfCFEaYKf1YgGjyWTIj0Ths89ekjfa1f6fUg629+B1U9Ca3XftSFViLB3dGnlCCRN6XeElP6nsJ7+LcBpQohPgNOA/UDCskgIcZMQYqUQYmV19VFwYyUKe870qJNlulhHY9+1KVVQJJNVIB/7UumrDuhoOO6etlCxNbDYO7q0chDJUvo+i6fv64BAoOdtSwCJkP5+YITl+XBzWxCGYRwwDONzhmHMBG4ztzUk8lnzvQ8bhjHHMIw55eXlXfwJAxj2ofbIE+WEpT3L+q5NqUJ/8vTV0Hugk75hFv6yKv00FyB0BpUVHU3y0dODQG7ALy0zpfQhfKH0XkQipL8CGC+EGC2EyASuBF62vkEIUSaEqhHKj4BHzP//A5wjhCg2A7jnmNs04sHvhYA3nPSHzYG0DNj9Qd+1K1Wwe/p9SvpK6Tf1XRtSAb8XDH846QsRWhxdQ6LTvA56ovSD81DcoXs8RcHcuKRvGIYPuBlJ1puAZw3D2CCEuFMIcZH5ttOBLUKIrUAFcJf52TrgF8iOYwVwp7lNIx7sJW5B+vvDZsPuD/umTamE3dPvS6Wpasx3NPRdG1IBb6t8zMgN3+7K0vaOFcnw9FVKslXppyiYm1DBNcMwXgdet237qeX/54Dnonz2EULKX0NhxyIYcixkFzm/HlzMIid8+8gT4cP7wdMa7r0ONEQo/T4M5B4t9o7HQWiAzOBJpNNt2CuP0eCpyW9bf0KQ9Fu7vw9l5ShPH/qP0tfoBbTVwWOfhU+eiP4e+1qlCiNPkpOV9i7vvfb1B0R4+n1Yf+doCeSq32kXE+lZiXW6b995dMwaDwZye+Dpe53sndQofU36fYG6nYAR8gadEFyr1Eb6lcfLUqwD3eKJyNPvJaX/yZOyzEUsqJsx1vkaCAjaO7ZrzpWZ2PFvq+n5LNVPA9R10BN7R4maDK30jw7U75SPsXr2aEo/Kx+GzBj4pO9tl0HrYGZDL3n6G1+Etc/EacvRbu9kJmbvdDQN/ICvzxMadfYkkOu1FBTsb4FcjV5AnSL9GCfZKZCrUDYBGvckv139Cb5OOfTt7TIA7Q3xvdmjxt6xLYqukJ6ZmL3T2TTwUzuVus8pk7+1uwFuq9JPT+2SiZr0+wL1CZB+tEAugLsI2gc4Aala48EZob3k6bfXxx+mq07haCH9CHsnK7FOt7N54JfBVtZOwRDzeTftLK30jzLU9cDeAZnx09k4sGvBqLokvV0GoKNBdih+X4y2HCVK3xON9BP09DuOAqWvSD5/qHzsrq/v6OnrQG7/wPrnZbXLZKJ+l3yMae9ECeSCVPowsElITVHvzXruhiHtHQgFMZ1wtEzO8kYZXaZnxu90A355DI1A7A70044Ou9LvJumHKX0dyO1f2PSyDPR5ktQLe9uh+YD5fw+UPgzsyUK+FCh9T6uc9az+d4IqTQADu5OFntk71sym/lqRs6MRWmt6to9eUfo6ZbN/QVkxSp33FPW7Q/8nEsiN5ulDSKUORPg65A3h6kXSt3aa0RSb3yNLE8DAJ/2o9k5GfCK3joL6q8Xzxg/hmWt6tg9F+sn09NOzAKGVfr+BCrrWbU/u/vIqEgvkqsi+FUeD0leevisdhKt3iMTaaUZTbKrzzSmV6mwgByq9bZCWHlotSyGR2jtW8uuvx6jpQM/XtVUjmmQqfSFSWlNfk34stNWF1F3djujv83XCxpelFRAPauQwaHJ8eycjR14QdhwtSl91eL1V8Ku9PvR/NHtHdcz5StkNYF/f2xZZdwcSq73T+SlQ+p7W0FKH3UUwe8ck/WR4+pDSmvqa9GNBqXKITfpb3oBnvyjXsI27z12QmQ9FI+LbO05+PhwdSl/ZO5B49khXYT1+0UhfjbjyB5ufGcAWj6fV2U5MpPbOp0Hpe1plOxMRZ9HQ2SxHnnmDzH32QOm7MiHNpOAULpmoST8WlCp3F8YmfXXBN+6Lv8/6nVAySiqqmEq/PTrpH3VK3907Kqgr9k6Q9AfwMbevmqWQyEjr0+Dpe8zsop6UT+hslrPiM/PMffZA6adbjnUKF0fXpB8LivRHnxb63wmqh25KwC+s2wnFo+VJjpUR5G1zVl0gP+vKHOAEZCX9BGeEdhVh9k400rfZOwNZ6Ue1dzITsHcsx6W/lmJQ57gnqbedzXI1t/QsGf/orr2jJh8qaKXfT1C/E/IGQ8VUqeK9UdSmCso0HYi9v4BfLnZePEqq+IA3+s3kaXNWXSB9fnfR0aP0XVm9ox4TsXdU/v5Rbe8kUHvHau/015RNdY57EpdRSl8Iqfa7be90hidp6EBuP0HdTigZDSVjACN62qY3QdJvPihviJLR8SdkxLJ3QPr6A1npWz39REv7dhXtDaEqnokGcgfyBK1Y9k7AF3sN145+nqfv94Y6rp6cw45GSfogH7sdyLUda630+wAfPwZL7g/fVm9aMSVj5PNovn6i9o6yiIoTIf3W2KQ/kJW+YcjjEpa90xuefr0snJXu1vYOxLZ3IDaZ9/dArrVTT4bSB1PpdzNP3zqSBVPpa9JPLT5+DBbfG6pn422XyrxktPyD6KSvCKnpYOzvUNlAJaPjz8KLproUsov6FwFtegW2LUzOvvxewLDZO72UvZNdJKtKRlNs6vzklst1DPrTMU82Ytk7ENvi6e8pm1bS78k57GwGd4H8PytPK/1PNZoOymDUoXXyubJySsZATglkF8dQ+iYxxLN36nbK4E/B8ASUflvs5RDd/czeee8e+OC+5OzLupQcJFb7pTtor5fnNTM3fspmZq4M4A1k0o9l70BsBd/RJNc/gP4ZyO01pd9dT9+u9DXppxYBv1T1ALs+kI+K4ItNlV8yJgbpm0rf0xzbL6zfCYUj5CzTeOVUYwVyQSrU/mTvtDckjxAVaWT0ciC3vUF2nrFuXms9GnfhUTo5K0F7J7c8/vv6Ctbz2+PsHUX6MUaI8eDtCL+/03XKZmrRciRUXyVI+hYrBmKTvs9C3LHUfv2u0P7ilVONF8h1m/ZOrOBasrBtITQfjv2ejiSSvteu9HtpRm7Q3smLE8gVsg3uwoGr9A2jh/ZOI+SWyv8HqtL3e+W9nqXsnfyeTc5SIyjQSj/lUFUv84fIZQgDfqnKswrl8B8k6TfudR7iKmKw7ssJKkcfYit9VdkxXvZOvHV2k4HOZnjy87Dq0ejv8ftkO5Kt9HuT9AMBqfSD9k4Mpa/KYQxk0vd7pfBxtHdM0o9l7/R7pW8l/W4GX9XnrPZOTwquhU3OypGdagrWyNCkDyF1PvXSkK9fZ86cVbVvSsbI2XwNuyM/7+2AwuHh+7KjvV4qy0SUvq8DMGLbO8Ga+r1s8dRslW2JpWgUESZrYZdgMSrz9/dGILezCTBMeyeGp++12GwDmfSDi6JHqb0Dsc9BR5PMhIJ+qvSTYO8ogZVlCeR6WrpX1sFpchakRO1r0ocQUU+/XD7u+iCUrqkQK23T2yYnXFn3ZYeyR1TqX6yTbA0eRoPKL+9tX79mm3yMdSOHlShOwsgjqPRNsumNQK5qc3Y8T789ZHkMZNKPtTxnPE/fMEecuSbp9+fsneyS7l+jTko/4OvetelUhgE06acMTQfkhV0xDUrGwo5F0LAnpMpBbgdn0vd1SPLIKYueq99eJx9zTN8zVspmrEXRFVJVdK16i3yMlSdv7XiSQYpBT9/8/enu5JO+KsGQXWwqtlhK30r6AzSQG2ultqC9E+Uc+Dok+alru7uLhfcm1PktGNoDpW8jffXYVV/fMByUfuoWUtGkD5L08wfLinejTobt78iL2Kr0c0qkx++o9NslQRUMja7022rN/SjSj9Gzx7oBFVJVdK1mq3yMRbrWGjbJIH27p59IGYCuQh23ePaOx2bveJr7bjnADS/Awjt6Z99Be8dJ6cexdxSJZhfJlOR+ae+ochpDul9eucNm76iia1319f1eaRU7Kf0UlFfWpA+SqAuGyf9HnRLK5LEqfSHk82ikn+GW+4g2QatLpB/jBlToT0q/I8lK37rABEibx+/pWUlcO4L2TrG8eaMtju5tD/ncylLrq7TNjS/Dsod6J2OrJ/ZOUAEX9N5Eunjw+2LHkzwtsm05JUlU+qrSZoz1lZ1gv75BK/2Uo/lAaFGEUSeFtisf3/q81mEFLZ+ZXlkwNLq9EyT9EvnoypCTWRztHVsg0wmpUPo+T6iTS1TpJ6M9Ttk7kFwyCdo7RaHYidMw3RrIVQqvr3z9jkbZOcXKEOsuYq3JHM/eUco5q6D3JtLFw19Oh/d+Hf11T2togl23PX3zc26b0u+qvRNcQMVK+ub/2tNPAQzDtHfMAGvBUOnfu7JCS6IplIyWXr9dUagywAVDpXfvdOLa6uQNlWFL0+puIDczVw6le1Pp1+8MjXpSqfTtefrBdXKTOPS12zvgrNjs2TvQd0o/kVXcuotYpJ+ovZOV33sT6eKhdjvsXBz9dU+rJGl3gWxvd0aN0Tz9rk7QsmengVb6KUV7vSQTZe8ATPs8jDk9tKqNQnaxJEHriQn45UWulD44+/ptdSFrRyEjO2TlWJFIIDcV5ZWVtZNVEEfpNxCcp5AUe8emhBIpA9BVdDRIgsrItii2KOci02bv9KXSB+fRZk8RS2i4zPIK8ewdt1L6KbZ3AuY9eXhDdDL3tISUvv0eThSdzbL+kiLo4HXTRU/fUenr7J3UQRF0gUXVn/EjuPrZyPc6KUJFUBnuOKRfG7J2gvuLovQTCeRC75dXrjFJv2Jq/Oyd/MGASC7pW5dLhK4ryP/+XAY/ndBeL4+fqosOzjevtR5NfyH9XlX6sWrvRLN3LAHOvlD6yl7pbHKeRwMhe0dZM93x9TubQrX0IeTpJ1Xpa9LvfTiRfjRkOJC+Nb1QjRaaHYK5bbUyRzhsf9FIP4FALqRA6W+VxeFySuLn6WeXmEPnZNg7dqVvPnbVK/74Mdj4kvNrajYuxLZ3PLaUTRjgpB/L3omSimm1PVx9oPStpHtovfN7rJ4+dM+iU6tmKSTV049TliWJSIj0hRALhBBbhBBVQohbHV6vFEK8K4T4RAixVghxvrl9lBCiXQix2vx7KNk/oMdQgddESN8p4GcNuqq4gFMwtz2avdPNQC6kRumXT4i/Rq1Szcmq/OnrkEHuNJd8Hi+QGA3eNmmrOaGjIRQMj0b6wXIYSun3YSDX2xFS0L1B+gnZO1GOvzWVMb0X0mvjwXo/Ho5F+nmWjjsJpJ9UT78f2TtCCBfwR+A8YDJwlRBisu1ttwPPGoYxE7gSeNDy2nbDMI41/76WpHbHx/Z3oLUm/vuaD0qfLq8i/nuD5GAh6qAVkS2He+7CGPaOE+l3M5ALvav0AwE5G7fsmPiTo5RqTtaMVZ+tAmFQaXaBTAKB2KSvOiqwTLKxkb6vE1kOw1S/WT2wBnoKdVwz8yTpJztt09smEwMUwVuRiL2TkSOrx/ZFymaY0l/n/B6rpw/dy9VX9o6CK0P+3mie/oq/Oa+tbc9Og35n7xwHVBmGscMwDA/wNHCx7T0GoLrAQqAXcsq6AF8nPHEZrHwk/nub9kPuIOeL3Q6ngJ9S6uoEFgyLJH2/V960EaQfZV3MWDegFb2p9Jv2y3aUjY+/cpVSzckkfWsFwkQKftmhjmt7NNJvdLB3bIrNbnmkufqupr76ziEzzLTNOAv2dBXRyiqDJaYSzd5pCpFpXwRyFenmDoqj9Hvq6TeHkz5EL6/c2QKvfRfWPB35mtNI3pUpxWc/If1hwF7L833mNit+BlwjhNgHvA580/LaaNP2eU8IcYrTFwghbhJCrBRCrKyurk689RZUN3fykxfX8/GeeqniDH8oNz4Wmg4kZu1AFHvHovRBWjx2e0flhNsDudGUfqwb0IreLK+sgrjliSr9JJK+vS6J6lC7ovQVYXfF3rHfvE7Bzb6qv6O+c+hM+ViX5AyeaGWVQQYu0zJi2zuKTPsikKvOW+UJsny5E6Enw9PvaIok/awodZtaTR5zes2enQbyGKeopn6yArlXAY8ahjEcOB94XAiRBhwEKk3b57vAP4UQBfYPG4bxsGEYcwzDmFNeXt6tBuRkunhm5V5eWXMgRLKJ9OZNB7tB+tbsHVuvXTA0claufWKWQtRAbpwFVBSyi+R07jgpY4Zh0NLZxdIB1Wb5hbJjYit9n0cGnt1FoU6op7ArfVcce8EJ6hz52sPtOAiVglb2jlOAHkLnxmqz9RvST7KvH295zliL01sVcG8tYh8LHgvpAxzZGP56ICCvUZWnD0lU+lEWR1fWsmMacJSYXYpq6idC+vuBEZbnw81tVnwZeBbAMIyPADdQZhhGp2EYteb2VcB2YEJPG+2E3Kx0ThlXxlsbDmMokk2kN++W0m9fS49TAAAgAElEQVRld20rTy3fQ6DTwd5pORw+FFZqsyuB3ERIP4FZuct21PK5P33I7F8s5GBjFy6omi3S/sg1Fw4PeJ2nuVurVcZbzWvjy/DW7fG/29cRPkW9O4Fc681mt3gUgWYX89Lq/fz2nR3Oi6M7Kf2ezOjsCdRxrpgirYBk5+rHG126MmPk6Vvsnd6okxQPKntIkb7d1/daYmSZedJGSeAcev0BVu6qw1C5/06knxVlcfSg0ncgfSelDylbHD0R0l8BjBdCjBZCZCIDtS/b3rMHOAtACDEJSfrVQohyMxCMEGIMMB7ohdQDiXOmVLC/oZ29B0xPPd6J7WyRAR0b6Vcdaeaj7Q7WkEn6uw9Vc9EflvCj59exeuch+ZryfQuGAgY0Hwp9zl53RyFWIDdeEBdi1t851NjBDY+u4IqHl7KrppVOX4APtiUQ2Fao2QZlE8xhZwylrUheBXK9rc7e745F8Nz18OEDkcrbDvv6od0J5Fo7U7vFo0aC7iJ+//Y27n97Gx5XTuQN6olm7/TBMpWWjoriUclX+rHsHaDJK1i7+0iIAMPa1tQ/lH75RHl+7L6+Oq+ZufJ6zspPSOk/8E4Vlz30ETc8uoIj9c1y1KiyfxQyoyyO3tZdpd8P7B3DMHzAzcB/gE3ILJ0NQog7hRAXmW/7HnCjEGIN8BRwnSGvjlOBtUKI1cBzwNcMw4hisvYcZ02qIE3A5h3mBI14J1YFwyzlFrYebuZzD37Ilx5ZRtUR28lMdxMgjVdWbmNQfhZjynNZuNaMzmdYlD6E+/qK9B3z9NsiZxEmau9EUfqGYXDLv9bw0fZafrhgIktuPZOS3EyW7ujCoa/eIkkfLHnyDhaPSaC+zAK8GeaNbz/uhzfCM18MBQRrt8X+bq+N9LswI/eJpbvZeKApjtKXx+uwN5sd1a240gTVngx87bZ2O+WuW+ydAw3t/PKNTXzjyVV89o9LmP/b93h2xV5nYuwpFOm7C2Mv3RkFR5o6uP7vy7n8zx/x+Yc+5At/WcqOalvqcZRrbm9dG3s9eTQdrOIXr26K/H2dzRZPvy+UfktopmzFtMhcfbNTWLi9lZ+9vEFWy40jCNs8Ph77aBfjBuXx4fZaPn//W/IFR6XvQOxxPX0Ruh9M1Hpc7D7cBWHWTSTk6RuG8bphGBMMwxhrGMZd5rafGobxsvn/RsMwTjIMY4aZmvmWuf3fhmFMMbfNMgzjld77KVCWl8WckSXs2W8Sbjylb8vRP9DQzrWPLCcrw0V2hovbX1wXdoE/8E4VLUYW44sE//7GifxwwUSam82hnSKGokr5WG+ZGRjV049STjXeUokKUZT+K2sP8kFVDbeeN5Gvnz6WnMx0jh9dwtIdCQS2AaO9HtpqaCsw1xDIiEH65nd/7d87eHx1Q2R7mg7K5RYzcuDyx+S2mm08/tEurn1kOV6/QxDarvSDpB+79s6WQ83c/uJ6vvvsavxW9WUP6Jud5PJD8rvv/fx0mvyZVO23rQNsKrKaznT++v4Oqo40g7sQo6ORZ1fs5dzfLeaRD3ay+WAzuVkucjJd/ODfa7np8VXUtCR70ZdGSRLpblkbqm4nbZ0eVu2u47GPdvHgoir21kVXiY8s2cV7W6tJE5Celsaq3fU89J7FIoph7/xnwyGWBiZzXPo2nlyyhV++sTmc+DubJJGCOTmrD/L0M/Okih88VXr61uQGk5Rf2dTIE0t348+Mr/SfXbGXhjYvv750Oq996xQmFMnfu71JhL8xM8o6ufE8/XR3cGavYRj8buFWqhoCtLW04A/0gmiwYMDNyD1nSgXeZvOAx1P6KuBaMJTGNi/XPrKc5g4fj14/l1vPm8TSHXW88InsGJ5duZffLNxKICOXs8flUeDO4JzJFYwukhOIOoWZXhkk/V2h72mvlzdUxHAuSm6utw2fK5ur/7qUX76xKXr7HZR+U4eXX7y6kWnDCrnmhJHB7fPGlkrrKwYxKKzZXAXAg6uaae30xVT6hw5LG6uqOZ2l+81gsTXQueiXknSvfhZGnwoijTWrl/OTlzbw3tZqZ8vJ7uknslwf8NhHuwDYfKiZldusIy1ne2fxXi/jBuVxyczh5OQVUlNXx7p9lrabSv+3i/byf69tYv5vF/PkmgaMjiZ++O/VTB5awNvfPZ13bjmdJ79yAi9+4yRuv2AS722pZsF9i9lyKHqAvbnDy6/f3MyH2xNUdh2NUuWrEt++ds762TNc+qeP+OlLG7jnzS2c+v/e5fq/L2fx1vAMuHaPn6dX7OHcKYN5+qZ5PHXTCVw2ezgvrj5AXat5TGPYO//ZcIg9hXPJNDzcOrWJhxfv4Pdvm6O1gF+SntXeSWKefiBg8OraAzS2RVqGrZ0+2fl0toTSqSumyvNWb8mPN4m31puBL2DQEMiOKQh9/gB/eX8nc0cVM3tkMeMG5fGHS8cD8OpmG8FnRbF3FOlHWw7VvL59/gA/fmEdv397G4UFBRxTloErTUR+JokYcKR/9uQKCrHU4ogFU+kb+UP4zrOr2VXbysNfnM2UoYVcOXcEMyuLuOu1Tby85gA/fn4dp4wvo7CgkDSzTIIQgvnj5MX+1MfmjZbhlnaRlfSdJmZB1KnXhqedlQc6WFJVy2Mf7pbE6wQHpf/bt7ZS09LJXZdMDbt4Thgjvz8Rtf/KUhkIW1+fwc3//Bh/mnMgdfOhJh5/dzUA37xgLnWBnIj2UFslM06GzID0LJqzh7N36xrmT6qgMDuDl9c4TOkwF6Xp8PpZuPEwHtLDvn/9/kZOveddmallorHdy/Mf7+fSWcOZNqyQRessN7219LOlfYv3+Thr4iAAhg4qp8Dl4Qf/XkuH1wxYm+flne3N/M8ZY/n5RVPwZeaThsH/nTeKp248gcrSEFGmpQm+csoYXv7mSfgCBr9+c7Pj8V2+s47zfv8+f1q0nZ++tCExO0iRPrDTGAzARSM6+MuX5vDhrWey5NYz+eYZ49hwoIkvPbKcd7ccCX70pdX7aWjzct2Jo4LbrjtxFB5fgKeW7zF/q7O9c6S5g5W76xk87UwQLq4bvIdLZw3nvv9uk99hLbYGsQO+IOMk+1fR4fXz4KIqzv/9+3zt8VXc//Y2PthWE3EsXlt3kJv/+QlffWIlPsuocP3+Ro6/+21+9cZmGUhVdXAGT5WP1mCuSfqFhcUU5WRwqDMjZgbWa+sOsr+hnZtOHRvcluWX+1hxyMfqvZbrO9MM5NrPoWnvBJw6BDMlubq5ky//YyVPLd/LzWeM45irf0Pahb+J2q5kYcCR/sjSXCqzTUXq94Ty6J3QdACyi3lnezPvbD7CD86dyInj5DqfaWmCuz47jYZ2L9966hPGV+Tz4NWzEJm5YYHIkQVp+HDxwKLdHGo0v6t4lAPpFwef+gMGP3t5A7uazAvFovQNw6CxqZHdTXDVcZW0e/28ud4SFAaWVNXwyAc76RDZIFxBpf/h9hoe+2gX1xw/kunDi8I+M35QHiW5mXwUh/TX7Wtkz145LeOK047l3S3VPLHKtD0sSt/jC/DlR1dSJOTN8NkTJuPKUZ2Q5YZq2ANFMvnr2RV7Wd5cygz3Yf549UzOmzqYtzYcCpGsgq+Tmk7BBfe/z42PreSmp0yP1t/JzppWrn1kOXvq2rjthXUcaZJt+veqfbR7/Vx/0ih+dN5EPO0mGYm0qPZOrT+HM0zSz8jOY0wBbDrYxI9fkLZewDzPxcVFfPPM8Vx74iiuPWMGAFfPKCQtiiKbOLiAG08Zwzubj/DJnlCHYxgGv3lrC1c8/BFpQvDlk0dTdaSFJVUJ2G4m6Xt8AX72gbxevj0zjbMnVzC0KJthRdl895xjWPyDMxg/KI/bnl9Hc4cXwzB49MNdTBpSwHGjQ/bi+Ip8Th5XxuMf7ZYWWxR7560NhzEMOPPYcTBsFmLXe9x1yVQmDs7nu8+s5nC17Fw6XLn8+IV1PLnqEAR8/O397azZ2xBB4sbqf2L8ZT5X/PZF7nlzC1kZaWw+1MRvF27lmr8t41+r9gXfGwgY3P/2NgqzM1i6o04SPNKGveHRFbR0+vj7kl10tDaFlH75JAzS2LVxWXA/h2uk6j5tyihOn1DO7tZ0jCiC0DAM/vzeDsaW5wYFARBKzMjM5w/vVIW2Z+XJtGnbaL2hRgqS1ubGyE7d106bkcGC+xazdEctd18yjVvOPQYxZLoUR72MAUf6ACOzLUojltpvPkggfyh3vrqRcYPyuO6kUWEvTx5awDfPHMf4QXk8ev1c8t0ZkYtoe9sRGVKVfuYPH8ib3JH0Q0r/g6oaHv1wF39aYipVi+9333+3gbeNCSMGcfclU6ksyQlaTCCH6t9+ejV3vrqRs367mM6MAhrra/j6E6v4wl+WMaQwm1vOPSbipwohOGFMCct21MVUln96r4ohmZLszjtuCv9zxlgWbjVJ3KL0X1q9n/0N7SwY64asAlzpGcwYJ+0kX6uphPw+2bEWjuDxj3bxg3+vpbNoHMMDB8hKg8/MGEqrx887m0Oq1OsP0N7ewmub6mn3+PnG6WN5f0cDAQRH6pv44t+WYQB/u3YOnb4At724nkDA4PGlu5lVWcTUYYWcOK6MyWVydBDIG4LRVofHFwj97vZ6OtOyyXG7mT1SzcrNI48O/nf+eJ7/eD9/X7KL9bvkjf6tc2bgzjDrAOWYi3+3xp5EeO2JoyjOyeCBhZvg0Qth1xIeWbKLB96p4tJZw3nj26fwgwXHUJaXyaMfOkzVt8Mk/T+8W8X71dkE0jLIaYmsKOnOcHHPZdM51NTBr97YzNIddWw+1Mz1J41CiPBO6roTR3GoqYP/rD8Y1d75z4ZDjC7LZUJFHow+DfZ/jNvfyoNXz8LjC/CrF1cAcNd/9/P08j3B+Ms9r63l4j8u4fz7P+CJpbtZt6+Re/+zhScXfoggwGSxh8e/fBwvfOMkFn3/DDb8/FzmjCzml69vot60nF5ff5BtR1r4xWen8qV5I/nrBzv557I93PDoCto9fv5+/VwADlVXB+2l93c1syNQwdZ1K1i1W9p6K7fJjuSsGaM5c1IFNV43/jZnpb+kqpaNB5v46qljQ516W51cprJwBCfMO5n/bjrM5kMmrzgUXfvvxsN0NMprOsPfFtaR+fwBNu05wu6mABUFbl795sl84fhKx7b0FgYk6Zenh0hUKZEIBPxwcC27/GXsrm3jjs9MJsMVeTj+d/4E3vrOqVQUmB6zfT1VbzuuzGxe+J+TyM5wccXDS9nQUSxXN1KjDFst/WdX7CUvK52DrSK4D4DHl+7m929vIy/Nw7FjhiKE4LMzh7Fke00wx/7xpbuoaenk9gsmUZSTwf6OLN5bu433t9Xwv/PH8+b/nkJhtnP5hhPGSF9/X72pSnwe+OB3we+vOtLCG+sPcValIrhSbjnnGI4ZLomu0QxaBwIGD723nUlDChie3RmMLcybLFca23PQjJU0HwDDz/vV2fzkpQ3MnzSIs089GeHvhIY9nDCmlLK8LF75ZG8w7/zXb2zG8HYyqqKEN79zKj9YMJG/fGkOHiOd51fsoL7Vw6PXz+WsSRV875wJLNx4mB89v06OACz2xWmjcuk00lnf5Gbxmi1MuP0NrvrLUupbPRjt9dQHcjl1QnnonGfmgqeFb505nnMmV3DX65tYumUvftJYMMMyTSXfrNHUEuW6MpGXlc5XTxvL+m07YNf77Fn9Nne/volzJldwz6XTyc1KJyvdxReOq+TtzUfYXRtnyb2ORhqNHB58t4qLZ44gLUba5szKYr588mieXLaHn760npLcTC6aETkX5YyJg6gsyeGJJVVyBrvN3mlo8/DR9lrOnTJYdhhjTpPv2/0hY8rz+NWl09l7UHaMzUY2z3x1HlefKL3vD753IndfMg0B3P7iej7zhw94cFEVw0xB8X8npnHK+NBEzNysdP7vkqk0d/j41Rubgyp/bHkuF0wbwu0XTGbOyGJ+/MI6qo608OA1szjjmEF84fhKOloaaRXZ1Ld6uOVfa2jIGMTw9Ea++vgq9tS2sWGXFFdlJaWcNr6cFpGD8DgvpPKn96oYlJ/FxTPN4xUIwAtflZl+l/+Da06ZRG6miz++u512j59dLfL6WbJhF1sPN7N4azXf+OcqykQTBgK38HLXK+vYW9dGh9fP15/8mOr6Bgry83nxf05ifEV+RBt6GwOS9LO8TbS4JBF99S/v8NXHV/K7hVv59ZubufOVjfztg53sXvo8NO3j/upZnDO5IuwCtCNMIdlJ3ywONqEin5f+5yRmVxbzl3XmxdRg+qUW0q9t6eStjYe4Yu4I5h87CoA1Ow/w6toD/PSl9Zx9TCnphhdhBnk/N3MYhgEvrT5Aa6ePh97bwSnjy6R3fPPJFJWUM6MMFv/gDP53/gQ5GokC5esHLZ6d78F/fyaL0wF/fm87ma405g4y5JTwzFyEEFx3ykQAXlghSWbhpsNsr27l66ePRbQ3QLb0mudNqsRnpLH3wMGw3//ntT4umDaEB6+eTUbFJPlazTZcaYILpw9h0LanMP4wl+Uff8xfP9iBW3g4bXIlBeZvOXNiBelZ2YzId/Hwl+YErasvnzyGGSOKeGblXsrysjhv6pDgby3L9EJmDjmF5UzI9/K108by8Z4GLv3Thxw+cpj6QA5nTbIM383zmpYm+O0VxzKmLJc0cxlMYV1MRxXmaw633JzwpXkjGZojfej31mxlVGkOv7l8RpgtdPUJI3EJwT8+jFIH3oSntZ6FOzooyc3kjs9MhtKxUBs9bfO7Zx/DqNIcth1p4arjRoRGKha40gRfmjeSjXvkb9lWH+BIc8jCe3vTEXwBgwVTZQyB4cfJoP7O9wA5UvviLGkZ3XXVicwdVRJMQyzPFnzh+Epe+9bJPP+NE/l/l01n2Y/nc8YI2Q5XdWSCwsTBBXz55NE8s3Ivd72+ia2HW/jWWeNxpQky09N48JpZzBlZzD2XTQ/er984Yyx5ooMNNTIgWtfqYcyY8YzPbqbTF+CzDy7B6Azl6RfmZFBUVIrL8EdYMuv2NbKkqpYvnzyarHTzeH3wW9j2Fiz4JQybTVFOJtfMG8kraw4w5Y43ufu/0gq9+8UVnPO7xXzpkeVMLjZIx48wswKz6eR7z67hur8vZ+HGwxxTmsGwshIy0/uGftP75Ft7GaK9nrzBU2D/Sq6YXsQ9W+v4z4bDZLrSyHAJWj1+xmXcT1ZaMW8FZvPmBfaioTEQofTbgnViinMzeezLx3HrfWuhGVoPV5FbMlpOADNJ/4VP9uP1G1wxdwQjPR7YCI8t3szLndnMGVnMA5+fCPcSHGqPKstlVmURz3+8D3/AoK7Vw3fPlvnzrjRBSekgSjoaIDc859cJytdfuqOWy+eMkEFWYO/eXTy0aR0vfLKfa04YSY4vvDjciEHSAlm27QBzDzTy4KLtVJbkcP7UwbAyVJfenZlOsyuf2urDGIbB2vXrmAGMnzCJ2648lnRXmizgBnLG74Rz+MyMoTSvWIUw/Cx+9XEmlV9AWrMRXoYBSM9wc8HEEjBjLur3/7/LpnPxH5Zw/Umjwm8ibxtZ2QWMG1kJe5dx63kTOWvSIG58bCUHmnbTTj6nTbCSfr7ML/d7ycvK4Mkbj8f/0lO4Dtosj1zzM3GUPkBOZjpfml0GKyAv0MyfvzgnolOuKHBz/rQh/GvlXr57zgTyssJvScMweGTJLq5pb8SfVcCzN82jKCdTpm3uWCQtNFfkbZyd6eI3lx/Lb97awrXzRkVt4+VzR/Dhx2uhHv62/AhPf/Q2w4qymTuqmKrqFoYUupk+zEzHzHDLWa873gt+/uKJ+bABcvPNeIFtIp8QglmVxcyqNG00ldUSpTDat+eP59W1B/nbBzsZW57LhdNDI5RB+W6e+/qJYe8flO+mLcPDO3UGbxw5xA8XTKTYVwlVL/DAlTO44R+rKHN7MXAhzLZVDh0MW+Bg9RGGDAsdm4fe206+Oz1kt+xfBe/eBVMvg7lfCb7va6eOpa7Fw5BCN6eke2Ax/O6SsWzKnEZti4eLR7TB35E2b9N+bps/gm++dpj0NMHvrzyWiuVGeHZaijHwlL63Xc6cK5b+8lXTC1l1+9nsuPt8tt51HhvuXMDyr43iNNda1g2+hJ9/9tiwDIy4yMwNX+LQG14GOMOVxlc+cwYAb3+4PJQumF2MYRg8vWIvsyqLmFCRT1a29AP9njbGlufx12vn4jZM39yyz8/NGs7Wwy088M42zpw4iJmVoaBw1NIHhhFRNkH5+ku31/LJnnpWfyL92H8t+ph/f7yPi2YM5dtnjZc3pnVOgZmyWZoZ4KbHVrFmbwM3nTpGkri1cBmAuxCXt5l/Lt/DomWrAPjhFfPle0HuN6cMamRtn1lD3Zzgkqpvnm8F9106MeL3yzZkOeZ/T6jIZ9ltZ/GN08eGv6A86pySYPbO3FEl/PvrJzIyrRZv/ghKrB2lrZjeoHw3Q3ICDu3IlB1iS3ylD3DBREmYJ4/IYNygPMf3XHfSKJo7fTz+Ubjab/f4+d6/1nDPq6vJEl4+O28yo8rMdg49Vo4yHRSzwuyRxfzzxhMYVBCdYArcGTxytcx4+fIZk7n9gknMGFHIku21rN/fxGdmDA0PWI8+FY5sgBYzpqGC9iplM1iRM0oGj5qpWr1Vdlg25GSm8/OLpuBKE3zn7AkJpS9mG+340nM4bnQJN506Rs67MfycPkzw12vnsGBCPkLl8QPjK2VHsnzTruA+dtW08sb6g1xzwshQx7xriQzSnndPaLUspLj7f5+fwXfPOYa5E6T1N6FIcPGxw7jh5NGUYvr9xaMAuHBiAbedP4l/3HAcFx87LHLyYYox8JS+Ss8zDzgdTRFZFoM2Pwlp6Zx9zQ/NZf66gAh7JzLVbdL4cXjTsjiyZzNrtu1gBkBOKR/vqafqSAv3XDpdvtH83K1nVZIzb560M+oiZ4FeOH0Id76ykQ5vgO/Mt5Uuchc5lw/e8Dy8/n34zoaw9s0bU8rr6w5xyYMf8mTmZkiDC8em8ZUvzA/aKREppqZC+szUUp5Y0U5ZXhaXzR4uX1MVNk1k55dQ2NzKt19Yzx9y6wm4B+HOtmWFlE2QZR4Asecj3HjYGBjJCa5NuPJUXZJwpR9rpmeBk6XlaZXHMKdUBvP9XnBlMLbIBdRz0pxZ4e+3FtNTJZejTZLLq4Dmw5HbHZAVkBZCuSt6TZVZlcWcPK6MX7+5mZ01Ldx+4WQaWr189YlVbD7UxI9OrYDlkJVn6eyHzZaP+1bC4GkJtSUqTBEzflgF4yfJuIxhGBxs7KAsz3YeRp8O3CktnmmXOadsQvQJWm21oVnNddtlFVcb5k+u4OPbz6YwJ4Fy5z4Pwu/h8pMmcfUZx8lOQt3TzQc5c+KxsMUfVtZkUJm0hpZt3MnZp55MTmY6D7+/g3RXGtdbkzma9ssRYK5DurWCulasc0FUkL9ICk/hbeXGU2da2pxgba1ewsBT+urgmwc8InvH0warn4BJF3Wd8EFePH5PqCSAml1nhRC4SkZzTFYtf18o1W57RhFPL99LbqaLC6ab3rN54gdnGyHicpj6X5STyRfnjeTq4yuZNtxW+yO3TBKvXTUd3iBvMFvd9YtmDOMrJ4/m91cey7xC2UFOyG0PJ84I0pe/b+7wHC6YPoQfLjgm5BG314cp/fTcIoa6PeS70zlzcAdpRQ6ZCWXjg0qfqrcxXJnUzrsNl+GDLa+ZO7Ir7Djlne3wtsnMCvtN2SgzKTJKRoW/36mCqrfdecJSXoUsqpdQO8z9xanX87fr5vA/Z4zluVX7WPC7xVz4wPvsr2/jkWvnctNc09KyjqhKxsiyHvtXJtaOWAgu2hP6rUIIhhZlR/rOQ4+Vs283vSxHk51NMm1YXa/pMeok+X3yehllVli3Wjy+Trm8pTk6TYjwITgyy8svDvnwqqyKuvZVWWX128z5DnsPHeKkX73Db97awnOr9nHprOEMyrfcy0374xdjLBwh04KtQXU1mjHdhohaU50tkeUcUoiBR/pK6Suysc/KXfcvqTKOu7F7+1cpWupm9nY4qsG0ktHMym/E0yR7/Use3cJzH+/jomOHkqt82+CMXMtFEWVR9J9cOJm7LnFQdDllgBGp9pXasCnSwpwMbr9wMhdPLiKt2UwZbbGlH9qyjdSNnObv5I9fmMXn54wItdXfGSJWAHcRY/L9vPWdU8lpPxDM0Q9D+TGyY2mthe3vICrnccq5l0kSUwuZ25V+epxJP3ZY7R0IHR+1cLa9M3Ja79TbHl3pJ0r6qhOJs8JZVrqL7587kX9//URystKpLM3h1W+eIucRWOvuKAgh1f6+VYm1w4rWWvj9DDi4Rj6PtT6uHWkuOO4rcu3hRb8KVZ5U9kes9XTVOaicJzuKw5YSyKv/CS9/E7Yt7NpvCa7Pa7HOLEofiCB9NSr5+bkjmFVZzAPvVOHzB/jqqWPC951IBd70THktWdc3UHELdY1ZhYRhyPNpXXYxxRiA9o55YeWWmxXwbKT/yeMwaLK88LoDdWMoG8Db5hyUKR5F7s7F3H3uUHgHrjh1Opta8/jaaRbvOT0LEOFZBMGKgAnGGdTQs7UG8iyByVYzQycaOanSvOnZ4e/xe8MCz/I9UcowBCtshnv66Z5GhuRnSVU96TOR360Kue1cJOuknH2nJJPx58Bac6Uh+/DXFWf1LjuUvaOK3CmlrzKq7KSvSMM6g9LTGn5MFfJN0jeMMK83ajsg4cqcMyuLWfidUwFL1pgT6YMk/ar/Opf8jYXaKjmPZOdiORmoK6QPcMbt8ve/9yspOtwWAotVBltNkisYCqXjwuvebzIL9+5dCscsSPy3qE4600L6eRWACJVZUevjKpiEOybPz9+um8vGA03UtnaG4iUKTQdg7KT4bSgZG0yKAKTgcheGxJBVSPg6ZJly+7lMIQau0s8uNpe2s5F+zTYYeWL8mzUa7Ovk+joireIIg/wAACAASURBVAiQMQVvK0WtcuLN9fNnc89lMxhZarmwhIhcMjHRRdEVcs1UUzWkVFDPo5K+WTtlxNzwiUbB2v+WQG5auhzC2m9kS4niIJRf23JYKvNCB6WvMniW/Vk+jj1LPlpvdrtl1tVl+JS9ozovRTgNe+QqUHZrL5q946j0B8vfZi/v4AS1P09L9OUGbRBChKcJqw7DThTD5wAGHPgkof0Godpdba6OluiazAppafCZ+2HK5+R1lmVpV6wy2EoB55ZBxWRpQar27Fws/9+zNPHfAaFO2qr0Xemysw4q/RZHpa8E4eShBZEp236vTMtNZK0NlT6r8v5ba0zR6XBNBTvwvlP6A4/0LdkyuAvCF0D2eeQNlOug3hKF3QaIVpJWBZL3r5KfiZaiZa+p31XVFZwhaiN9ReTxlH7lPPlb1IXpVPtfCNNTtynt4AIqVnunUL5PKR8nT79whNzf3mWSQCumyO1jz5KEDA6k7+5ayd6o9s4eKBwuRxZWOK5/HI30VdpmAhaPtUPv7iL2sZQ+yGBul/ZntsMMpjsuFhMPaS743MMynXHUSaHtsdYzVkIkp1Se84bdcpSy5U0I+GDkSfJ+iVU6xQ61gEmmbaSTPyS6vZOZD4jYBRlbDgMGFNpXhnVAyVjZDnXPtVaHnAb1/QrqO61CKcUYeKTfXi8JIjMnUumrk5LXvSUZgcje29se1d4B4NDayJLKVthXy+kq6eeapG+vL6PsnWhZJjXboGB4qJ0q7zzagi9OKZNR7B0gFKRzUvppLig11f7YM0OjLndBiEDsx9TVRaUfy95RATYrHNc/bnUmwqBnnEDaZtgNn2TSzymRAd39XfT11XlTwfSuXnMKrgy47G9w3q8t22Io/eC1VQaDzI7+yCbY9Ipch+KEb8gR1MHVibfBSemDSfrm+bGTflqatMNilWhpMuNdBQmQfuk4+aiEVGuNvH+crqlo5zKFGICkXxdSnu6CUKAHoNUkNjWrsjuwkr5hyPQrJ3tHKVxfR+TiKVbYV8uJEsiNCrVvq9L3eUIjnKhKvwrKxkVONlI3pupMFNLdkSWgg7aDhfTVsVdVDp0CuRCyeMadFb79mPPlY6btJo61Tq8dfq/0TTPzZOef7g5X+k6jj6hK3+Hc5pmkn8AErfDFXHpA+qqWvh3D5nSD9E17p71OXjddtXdiIdaCN0qI5JRIewdg73LY/raM/ajlDvd8lPj3OXn6AAV2pW973cn6tcK21kZMlJoBYBXMVUrflSmtUSd7pw8DuQOQ9BtCRGhfz1RlqSTL3vFFTqQKvS8nRA5OZZUV7PZOVwO5rnRJtFZP36r6nSYRGYYk/dLxIatCdYjWIbgVjkpfxU8clP6hdWZcJUqAcfBUaeWMOSN8+6xr4bJHQsHe4O/sQp12+zHMKZVK39suO0FH0repMsMw4wIORBi0d1Ko9FUtfTuGz5Hk1mhftjrW/iztqNkqf2daulTuPYXah6PSr5G/w5UBhZXyXlr6oOzMJ10khUbpeNizLPKz0RBU+g72TlutvGbtnj6YgjAG6Td2gfQLK+Xxq62SKaftdfK3CBE5r6dTK/3ko70+pDbta2EGlX5P7B1LmmU8L1RZJzFJP0og12n0EA05ZeFKX9lY2cXOarS1Wl7wpeMi/elgTMQ2OknPjpK9I8IDeepirt7sbO0oHP91+Nr7kRNfMtww9dJIgosyI9cR6iZTo6XsEvm7zBz94BwOK9Jc8jdaYzXgfG6z8uW+E1H63rbQhKXuKv3OpugkoXz9ruTrtzeE2qRIP9ai6F2BK7wMQxjaakMxqLQ0mUXXtF+qYqXyK0+QGTzWla9iIejpO9g7IFewM/yRpJ9VELOmPk0H5DlOxHt3pct7vXa75B8jEEqwyMxzVvo6kJtEtNWFatfbe3N1k/ZI6VvsHUWCPSJ9u73TJofxaV04Nbll4epeqfWKqZLgbeUYggG8snHmTShCo6C2WnlDpNtq+TiRbkeDJCNrWxU5+T3OilohMwcGJZAOZ/3+RAO56ngqIsgpluorWo5+sE254bEacLbZhDCzQxJR+i0hXziRbB8nWBZQicDgaZLAu2LxdDRA2TGyk6veGndR9C4hPUYZhtaacNtQWTwTLwgF1ivnyeOk4g3x0NliWl+267XAJH2VpWbvFOIpfTUxK9Esv9JxcoKWElzqd5rVW4MIBnK10k8ewpS+mUmi/MXW6pDP211kWGyAeKo8SPrxPH1bILerAbWc0vC0S+WdVkyRqsOe2aMya0rHSZWSU2qxd2qd2+uUvdNeH27tQPjFHIv0u4qurL1qt3eyS+TvipajrxBG+nGCm3mDE8ve8bSFSL+n9o4T0rMk8XdlklZ7vewIy8aZSj9KllJ3EGtpS/tM7wpzlSvrXI6u+vpqfVw7lNJX13p3lH4i1o5Cydhw0lcjGvtIvqMxfAZzH2Bgkb5hzkxV1oQtH5eWI6FhV3eRnikJyNNqUYNR0jETIn0He6erF0RuubO9o1Ih7eRUu03enMp+yRsUHsh1GplEy96xD3+tz2PZO12F+v5Elha02zvK06/fLeMIeVHKb2TlhzziePMl8hOcletplR1jRm7PArmxAn/D5shcffuILhrUeSubYLF3kqX0YwRy7dfW9Mvh/HvD4zolY+RIPNF8/c6WyMwdCJG+GtXaSb9iipygdsi52qck/eGJtQFkMNfbFkpgiGbvKKuuu/OEkoCBRfreNqkwgkrfvFFUj956xHmGZVehFGHQ3olyw5SMlo85Zc6vg3Mgt6sjkdwy2dkpH7StRqoJFQyNIP3t8uZSQ+qESD9Knr41Rx9kB6jUXrTMne4gPQswZD53PETYOyWyrfU7ZZuiWWfWobjX1nHYkWjRNa9KHS3qHaUPUD5Bfk+c1bxC+zPPW9kEOfppq02evZPmktee3YozjEh7x10oy6FY50wIAZXHS18/EXhaInP0Qf4+V1YojdI+Gph9neyIP3wg8rMBvwyOd1XpQ6izCpK+3d6Jcy5TgIFF+taJWeCg9Kt7rvRBXiweSyA3WpnU4cfJmYsTYkwrj8jT70YFvpwyaeMoz1jdXCqf3E76Ndvk0F4hd1B4ILcrSt9u70Dook6qvRMjQGiHk71jBKQKixlncPD0o5FhXoXMxLCnsTq1JTPXrIbaS6SvVK2tuJ4jDCNky5WNBwypdpNpNzhdK51NMo02VnxLoXKeVOFNCfyezmZnpS+E9PVroyj9nBKYfS2sfw4a9oa/1nJEBn+7QvqlJunvXQaI0Ojenr3T0dSnQVwYaKSvSC/HkrIJoeBJUpV+S2jmYDSSTkuTF1asBRMy7fZONzIpghO0TIuntUZ2BE6rPPl9UvGqyVEgj0lrtSSErij9ziZn20ERVLLtHega6QftHfN6qNvhnLmj0BVPP1qHGtEWM+2zu0rf2yGPe0zSN9uSCEl62+Vo2F0kg7kQGo0kCy6H4njWiVnxMOJ4+ZhIcDqapw+yM1SjH6fU2xO+Ia/5pX8K3x7M0U9gYpZCwXApTFoOy+tNjV4iSF8r/eSiPYbS93uliu1J5o5C0N7pYp0cJ2RkS8tC1WXxtnVD6VuKroEk/9xSuZ+swvDUwobd8vtKLUo/b5D83tYa+egYyHVSb1EKfbkLw8saJwPBxTkSIP0Ie8fSicVU+vmWWjlx0nGDHWoM0g8EJKFm5spj0Z3snc4Esj26ovStpTNKxwKmt5wsewecg+4qucA+6c8JypZMJIMnmqcPoeMCzqRfNEKuCbDq0fBz05WJWQppadIyhXA3wSllsw8nZsGAI301WchB6bfWAEbPcvQVFOkHs3d6sAqOvbxyd+wdJ6WvLrz8ivBJROpGKrMqfZPAVNVDJzVmV/p+n2yzExnllkPx6OQGq4KVPrtp7yjEUvolo6FxjyzvG29mtDpmsSZo+Sz76K69E8zrjpEvnjtIFsRLJIXUWjojIzvUCSbb3olQ+mrSX4ykBgV3gQy2WytXRkM0Tx9spB+lYzjxW7JjXvHX0LaulGCwQlk81vtHuQIqAaGzqU/r7sBAI/0IT98kpM6mUEpiUpR+XnjKZk+VPoT2Za8TkgjsRdfaakLb8irClf7BNZIgBlnWBVYdxBFz6b1EPH01KcZJ6Z97F3zuz137DfEQK//bDkX6KpU2xzLiiKX0T/ymTCN8/qaQFxzX3okxQSvY+fTA3kmkVosrXV7XiSh9e2VUpaqTbe/YO+eu2DtgLrSzLf77Yin9gjhKH+TM8HHzZcVXlXHUtF9aNYl0UFYElb6N9I1ASDB1NGpPP6mwllWGECF1NIUmHyXF0zd9+KSQvkXpd7bIG7ertYGs9o7PIy8sdeHZJxEd+ETe6NYbRR2T6likb1P6Kk7iRPqlY0PposlCMJCbQP0dFRdRWTphSj8G6Wdkw+f/ITuW938b2uaEnNL46tpK+u4i2a6uFI2D6GWV7cgf3EV7xyR9tVxhMu0dJ6VvLaucCErHRSp9bzs8MBs2m6urGYYUH7E8fYVYndrcG6X3X2Uu4KJy9Ls6UlVK327vgLwW/D4pFrWnn0S018uTqwKnrgz5PEzpJ9HeUcP3rpRMsMOq9HcskjfL+LO7to/0THkhtdVEFkzLGxyuRg98AkNnhn8+aO/EIX3DH1qWsTOG0u8NBO2dBJW+lcTchTKN0JUZv0MtGwcX3S9/K0QnizSXmfWUIOkrko2m9tvqIpfVg8SrMhYM7aK9YwojZfP1eiC3Rp7DRL+nbLyM0VnXnj20XnYEOxbJ5942qaLjefpWAeCEcWfJc7n6n/J50wFZfrurUHGyMKVvWXQpkfhMCjDwSN9eM0bNvFPElxSlbwZnvB2AiFzarysIKv122PqmbG93VvVS9XeC3qlF6XtbJUk3HZTZBXbSV6r1yObQczvSbUo75aTfhUCux5aNIswUusIYOfpWTL0Ujv+aXGvVPr3fCuv8BicEM4ByQyQbzdd//LPw6PkOpS4SJP1ElX4q7B2noH9bnbwmu1LWAMItHlVyWcWl1ES6qErftODi2aWuDDlRbOub8h5KZG1cJ5QdIzs8NSnT+t1W0teB3CTCWoJBQdXYaK2WijzaBdIVBFM2zUybngQsldL3tMC2t6Tq6E61w9wySfj22h9W71mtsDTk2PDPprnkDdnZCAjn3Ht7IDVI+ilSLV21d+znObfcuY5+NJz3a/j2mtjvyR/cNXsHnDN4DEOuYnXgE1j40/DXEiZ9S1XJWOhQRfJM4qmYIkc/yuZJBpyUfmtNZHG9WAjWqLeSvnk+qk3SV5OeogkPpfQTiZEd+wWZ1bb2GSmOukP6eeVw80q5sIyC1d7pB7X0YaCtkdtWF0lYqm6264g8KcnIKFHBmfaGnmXuQIj09yyTKjzWRK5YyCmTE1pabQEzNbJpPiRJRaTJWi125A2SFlh2ceSqUuCg9GN4+r2BWCsy2eE0q/nC33U9QB5L5YMky4NrY7cDZFtUR+5k77RWy+NaMByWPSRXkJp8kfn+Rlk6Il7cyLqwS6zOTU2oUyMedyHckmBxs0SRngVtreHb2moSm5ilUDRS/m6r0j9kHuvmA1J0KOERTchl5oTSh+OhYooUQx89KCeRdTVzR8F+7K0lu33mffVpCOQKIRYIIbYIIaqEELc6vF4phHhXCPGJEGKtEOJ8y2s/Mj+3RQhxbjIbH4H2+siIe1DpH0lO5g6EJk+11fR8WKw+v+EFQMhMgu4gt9T09G0BszzLJKKDq6F8knPQTnUO0W5M++LoqSZ95cU3JVA33m7vgCzk5dTZ9bRNrUei17yxzhcIKn0H0lczQhfcDUNnwUs3Q51cWzlmLX0r8k1lGs/X73Col5RsOK19YC2rnNA+0mU2jArm+jxweKPF9tlqUfoxSD1/aOKd/bFXQ5NZfrs7St8JVnunH1TYhARIXwjhAv4InAdMBq4SQky2ve124FnDMGYCVwIPmp+dbD6fAiwAHjT31zuwrpqlkJUvFUFLdXL8fAidyNaa2LNtE4FScNWbYPjcxLMb7Mgxyyu3HJFBS3VjB/PJD5tB3GOdP58bj/T72NPPHyJ/k1pMOxac7J1eadNg5yqmCoqUVO0dcFb6quRzyVj4/N/l/w/Mgj+eINePTYQkgko/jq/vVBk12XBlOE/O6uq1bU3brN4sFfi0y83nWy2efoxrcNplMPH86K/b36vWaE4a6TvYO58CT/84oMowjB2GYXiAp4GLbe8xAPVLCgFzdgMXA08bhtFpGMZOoMrcX/IRrCniFMhtSl4JBgiRfltNzzJ3IFyRTujBQCi3THqSdTskcQfTFYvlhbz/Y2kj2IO4CvGUvuqcrKQv0pKzxF4iEEIOwdUEslhIZn34WFBEW7MlSjssyxAq4nZS+o2m0i8aIYOAN7wBp37ftAoMGJHALZPorNx2hyJ5yYZ97QNfp0yt7Greu6pR7/eF/PzJF8tVqhJV+qfeAid9O7HvyymBY86T/3fX3rHDau/0k+ydRDz9YYC1ItE+4Hjbe34GvCWE+CaQCyiPYhhgLZe3z9wWBiHETcBNAJWV3SzS1dksSS8ikFso1ZXfkzx7R/XerTU936fVq+2unw+hoXP15vC01LQ0SegqBzku6Ue5Me21bzqapMpPZYnYQZNhzdOyg4/1vU72Tm9g9KlyJPWfH8ON70YG4K2B3DSXVKSOSn+vvE4VGVRM6fo8h5wSGUCNR/odDckthOcE+yL2XZ2YpVA2Xqr7ht3Sz8/Mk9lGJWMl6avfkcxR3Rm3ycV9kpHaDeH2jrL7PgVKPxFcBTxqGMZw4HzgcSFEwvs2DONhwzDmGIYxp7y8mwc74JOLapdPDN+eVSDVqRFIvtL3tPTc3lEjhYJhPZvQpIbOtVWRWRJ5/7+9cw+uq7ru8Lck2Ve25YeMbdnYYBwQYAMhOC6B0KEtmAZIG5iGzNjDNGRKStOGQBmSDH0lDdM0ybRNmsykmTBAQzMdk8bNtC5loCmP6SuADU0C2MEYp2B7MDYJMn6BZHv1j332vece3SsdXZ2XzlnfjEa6515J+2hfrfM7v732WgPuLki62/8Of/Fqdws+ytM/mP2bd+Acpxh9M5R2ZGXv9M519eD3Ptu6RO/I4aALWuBozpjXOnvnwC7XZ3UyiIyfTQTZ2DtRpT/RjVkeXxTwZzuc0l98nhMxCwZdtlMcpT9RFp0Nv/KHyYkZLz68pz+9z61X5EicwLwHCJdLXBYcC3Mj8A8AqvoDoBdYEPN7k2HmfFi/Ac781ebj4ZXyxK7eIRU5WXunu8cFj7Ountwbzf9DnTg2WlF5X3/RqjEKiMX19H3K5pvZ+fkef8Eay+JRzc7eAZdls/LX4fEvwuuRHaTRO4529XeGXkmm98DsJWMrfdXWjW+SprvWKCAIob0jE8jegcbGsf0vuAvrkvPd44VnuUqx/gKaxQW+U7qCLlkjhwtRYRPiBf3NwKCIrBCR6biF2U2R17wCXA4gIitxQX9/8Lp1IlITkRXAIPBUUoOPRViNJqb0Q2+yyZRg8PzWw7D2s5P7GeFAH724zQ6CfrtFXGgsXLWzq4qg9H1P3bEWc0eOApptO7qr/sL9ff7lluaG3sORO45W9XdUnb2TRBnq2YvHLq88fMjtNE7d04/U3vG7aidq78yc79botj/k7t4Wv9MdX3CWEzd7n3Xz3CrFuEhMmxko/aHcrR2IEfRV9RhwM/AwsA2XpfO8iNwpIkEyMbcDvy0iPwI2AB9Rx/O4O4CtwEPAx1U1Zk+3hGhS+gnbO5BM0F+0cvKqOXzrHL2N7osR9BecCdfd28gPj1IEpV+b7fK3xwz6kbLKWTBnCVzxp/Dyf8PL/9U4PnxodDmIqNJ/a8hZVkn47LPHKcUQrrCZJt0158X7C2Cn9g44tf/y/7ivvdL3dwB7ni62yvf4si2+VWLOxDKXVPVB4MHIsc+Evt4KXNLmez8PfH4SY5wcTUo/4cUZSCboJ0FPLagHf3D0bbTPMlnSZhEXnLV07gfH+PktlH7/is7H2ynjZfBEyypnhe/zGu7CNHKk+b0yo3+00h8KZe5MltmL3fy363MQLcGQFuGKqF29QfvOrs5+70mDrhtVd62xa9iXjjjys0ZlyyLjy7a8dWDixRRToFxlGFrhlX53Lblbq2kJevpJ4hdwo4pq1bVw5ZfaZ+7EoVUZhqyVPrh1iddfbF9uINo1KyvC+yHCYwmPo9VCrl+UTsTe8WmbbdR+uIFKmkQb3hza5yzHOHWPovi2ngOrGtlRtb5G0/Ipo/QPBa0S81f65Q/6PtD3LUpuRd4vzsDks3eSxHumUU9/5ny46GOd/dN5opuz8ur1OXCO86X3t8mNz8PeAXdnUZvTXIBt+HDzOHrnub+fb7MJoRz9JOydcTZoZWnvQCNt89Ak9sj4DB5v7Xi8xVMAj3xcvL0zhRZypzb+j5xU5o6nHvSLpPSDoD/RBbM4hBuTHx9xZaXz+IcbL4MnL3sHgqqbIZUdzSJqtSt3aJd7L000s6UVc8YpxVCvzZ+VveOV/mud2xp+8f7k1c3HvdWTZLpmWkyf5XYPt+spnTHlD/phpZ8k3qstkr1TV/ppBP0etxPy2FvZl2AIM/90dwFqt5ibl70Do7uURT39VvV3DrzirJ0k7kLHVfq+yVBWSj9s73Sq9E+HjzzoqmCG8Up/Stg7fW43/IljpvQzoWe686OTVvr+zVYke2fuMlcMLi0l19Pr/pHzDPrdPbDwzPZBPy97BwKlH/X0wwu5rZR+Qjn64OZjel/7tM2jQ+7CnfbfJryQqzr5EiinXTJ6t/OCKab0/V6FnCtsQhWCPsBlfwwX/GayP9MruDwUZTsu/j346Pcn592PRU8tUPo5N4MYOLeg9k5I6dc3iUWyd6BZ6SeVo+8Za4OWr7CZdukMr/SPDzdKoCSdtVJv8zhFgr6nAEq/XPX02/HeTyT/M+v2ToGUfrh+Sxr4Prl5Kn0IavBsCLoxRWoF5WrvLHIXxOEjLkVRj0fy9CONVIYPu8qwSdbCGasUQxYlGCC06D+cbMe6MLMWwnkfaqTKFplw0M+q6dAYVEPpp0ERlX7a+DZ49aCfl9IPKnu3UvsjXunnYe8EavbwvtY2k1f63t4ZSjBzxzOW0s+iwiY0p2ymFfRF4IN3w2CH/SeypGBK34J+p9SDfoGUftoURekvPt8p6e0Pj35u+IjzrcfrepUG9Vz9fc219D29c924fY34JHP0PV7pq45+LosGKtC8e9uvcRRgU1JuWNAvCUXM3kmbnprLMa83g8gp6PcthHN+Azbf06jr4okunmY6rkDNHnqtuZa+p6sbzl8Pz9znmtAfCIJ+Ugu54NI2j78NXzgF/mwAvnxO4yLjWyWmjV90PR6yd5JOpJhKNO3VsIXcqUtd6Vcp6EeUfp5v4Es/6aycJ/6m+fhIZPE0S8K7csO19MNccacLAg9+0in9rmmNlpZJsOpauPhmWP1heM/vuAv0v/2Je+7oGxnZOyGlf3ifO8csfm9RKZjSr8ZCbhrUUzYrFvSHDze6ZuW5nrFoJaz8ADz5TRfkvIIdPpJP5g4E+yTEqduRNkF/1gJXUfWB22DfNpdmm2S21Zwl8L5QqaveefDI52Dn48GO0AztHa/0+waybbZTNPx7oGtaIRI/TOl3ig94BZjEzAgr/ay7ZrXi0k+5bJmn7mocy6prViu6e1xQP7h37Cyi1TfA0ne73O0krZ1WXPS7bs3gX28HNCN7J1hP8Z5+0ou4Uw0vEOM0uM8AC/qdsmiVK2WbxPb5qUI9e6cY28lZ8k448ypn8XjLKauuWe3wufrDY2QRdXXD+//K3S3NW57ueKbNgMs/47pPQcZK/+3J1d0pC17pF8DPBwv6nTO4Fm7flp+VkAdNSr8Yb2B+6VPOq37sz93jLLtmtcLvyh1vk9jJF8D133VrE2lz7nWNCquZpmyOWNCHUNDP388HC/rGRGhS+jll7kRZ+m648Can9p/5dr72DrhF2UPhPP0xFpXPWAv9p6U/pq4uV1p7xvzGTtY08Up/5KirOVPldE1oZJMVJOjbQq4Rn7DST6OSZ6e87wvw+na3ONo9zV0I8sIr/bd9nn5OmURRTn0PfHpnNp6yz945uNftSk6qY91UxV/4C3J3bErfiI9X+m8VSOmDW0D90Legf3ng6eep9Adcq8A39wTZGjlsEmtHVouIXd2ANHoFVN3e6amBdBdG6VvQN+JTV/o5NVAZixn9sP47bmE97cXRsfAB7uc7q7XeE0bEBTpfZqLq9o6I6wPhe0HkjNk7Rnx6aoC6hdMiKX3PgjPgtq0NTzkPfIB746dTowJkWnTXTOmH+dh/5j2COqb0jfj4PQknjhXGnxzFtN58c6F90D+wu1rF+KL0TG8UlrOgXygs6BvxCSvoIir9IuADnJ7IrxxEEfCLudNmVvuOp4BY0DfiE959XFSlnze9cxsBr8pB3y9g9y0qxC5Uo4EFfSM+TUHflH5LRBoWT5WDvr/wVT1ds4BY0DfiY/ZOPLzFU3VPH8zPLyAW9I34mL0Tj7rSr7CX7ZV+1dM1C4gFfSM+YaVftDz9IuHVbVXz9KFRf8eUfuGwoG/Exzz9eJinb/ZOgbGgb8THPP141D39Cgd9s3cKiwV9Iz6+S1jeXbOKzuyg/aEpfcveKSAW9I34eKVfhK5ZRaZu71T4wlhX+hb0i0asoC8iV4rICyKyQ0TuaPH8V0Tkh8HHdhEZCj13PPTcpiQHb2SM9/RrxagWWFhOOt3VyR84L++R5Id5+oVl3IJrItINfB24AtgNbBaRTaq61b9GVW8Lvf4TwAWhH3FUVd+V3JCN3AgrfaM9M/rh1h/lPYp8qc1xTVu8JWgUhjhK/0Jgh6ruVNVh4H7gmjFevx7YkMTgjIJRV/oW9I1xeO8tcP3GvEdhtCBO0F8K7Ao93h0cG4WILAdWAI+GDveKyBYReUJEru14pEb+dJvSN2IyHUZEzAAABSBJREFUZwksy7GDmdGWpOvprwM2qurx0LHlqrpHRN4BPCoiz6rqS+FvEpGbgJsATj311ISHZCRGV5fbdGMbswxjyhJH6e8BTgk9XhYca8U6ItaOqu4JPu8EHqfZ7/evuUtV16jqmoULF8YYkpEbPb3VLi9gGFOcOEF/MzAoIitEZDousI/KwhGRs4F+4AehY/0iUgu+XgBcAmyNfq8xhVj7WVh9Q96jMAyjQ8a1d1T1mIjcDDwMdAP3qurzInInsEVV/QVgHXC/qmro21cC3xSRE7gLzBfDWT/GFOQXPpr3CAzDmATSHKPzZ82aNbply5a8h2EYhjGlEJGnVXXNeK+zHbmGYRgVwoK+YRhGhbCgbxiGUSEs6BuGYVQIC/qGYRgVwoK+YRhGhbCgbxiGUSEKl6cvIvuBlyfxIxYAryc0nKlCFc8ZqnneVTxnqOZ5T/Scl6vquHVsChf0J4uIbImzQaFMVPGcoZrnXcVzhmqed1rnbPaOYRhGhbCgbxiGUSHKGPTvynsAOVDFc4ZqnncVzxmqed6pnHPpPH3DMAyjPWVU+oZhGEYbShP0ReRKEXlBRHaIyB15jyctROQUEXlMRLaKyPMicmtwfL6IfF9EXgw+9+c91qQRkW4R+V8ReSB4vEJEngzm/DtBk59SISLzRGSjiPxERLaJyMVln2sRuS14bz8nIhtEpLeMcy0i94rIPhF5LnSs5dyK42vB+f9YRFZ3+ntLEfRFpBv4OnAVsApYLyKr8h1VahwDblfVVcBFwMeDc70DeERVB4FHgsdl41ZgW+jxl4CvqOoZwBvAjbmMKl2+CjykqmcD5+POv7RzLSJLgVuANap6Lq5x0zrKOdffAq6MHGs3t1cBg8HHTcA3Ov2lpQj6wIXADlXdqarDwP3ANTmPKRVU9VVVfSb4+iAuCCzFne99wcvuA67NZ4TpICLLgPcDdwePBbgM2Bi8pIznPBe4FLgHQFWHVXWIks81rqPfDBHpAWYCr1LCuVbV/wB+Hjncbm6vAf5OHU8A80RkSSe/tyxBfymwK/R4d3Cs1IjIabhG808CA6r6avDUXmAgp2GlxV8DnwZOBI9PAoZU9VjwuIxzvgLYD/xtYGvdLSKzKPFcq+oe4C+BV3DB/gDwNOWfa0+7uU0sxpUl6FcOEekD/hH4fVV9M/xc0Ke4NGlZIvJrwD5VfTrvsWRMD7Aa+IaqXgAcJmLllHCu+3GqdgVwMjCL0RZIJUhrbssS9PcAp4QeLwuOlRIRmYYL+H+vqt8LDr/mb/eCz/vyGl8KXAJ8QET+D2fdXYbzuucFFgCUc853A7tV9cng8UbcRaDMc70W+Kmq7lfVEeB7uPkv+1x72s1tYjGuLEF/MzAYrPBPxy38bMp5TKkQeNn3ANtU9cuhpzYBNwRf3wD8c9ZjSwtV/QNVXaaqp+Hm9lFVvR54DLgueFmpzhlAVfcCu0TkrODQ5cBWSjzXOFvnIhGZGbzX/TmXeq5DtJvbTcCHgyyei4ADIRtoYqhqKT6Aq4HtwEvAH+U9nhTP8xdxt3w/Bn4YfFyN87gfAV4E/h2Yn/dYUzr/XwYeCL5+B/AUsAP4LlDLe3wpnO+7gC3BfP8T0F/2uQY+B/wEeA74NlAr41wDG3DrFiO4u7ob280tILgMxZeAZ3HZTR39XtuRaxiGUSHKYu8YhmEYMbCgbxiGUSEs6BuGYVQIC/qGYRgVwoK+YRhGhbCgbxiGUSEs6BuGYVQIC/qGYRgV4v8BbU3zlMr45qwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 1000\n",
    "res_IMC = []\n",
    "res_MC = []\n",
    "for i in range(100):\n",
    "    res_IMC.append(IMC(N))\n",
    "    res_MC.append(MC(N))\n",
    "\n",
    "plt.plot(res_IMC, label='IMC')\n",
    "plt.plot(res_MC, label='MC')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([93., 38., 37., 26., 21., 14., 10., 24., 23., 15., 23., 14., 19.,\n",
       "        20., 13., 14., 13., 14., 11.,  9.,  3., 12.,  9., 11.,  9.,  6.,\n",
       "        16.,  9., 10., 12., 11.,  4., 13.,  6.,  7., 12.,  7., 13., 12.,\n",
       "         6.,  6., 10.,  7.,  9.,  6.,  9.,  7.,  6.,  6.,  7., 12.,  9.,\n",
       "         5.,  5.,  7.,  8.,  7.,  5.,  6., 10.,  4.,  6.,  3., 10.,  7.,\n",
       "         4.,  5.,  2.,  7.,  5.,  7.,  3.,  5., 14.,  2.,  4., 10.,  7.,\n",
       "         5.,  6.,  1.,  6.,  9.,  2.,  7.,  4.,  5.,  5.,  5.,  3.,  4.,\n",
       "         4.,  2.,  2.,  9.,  7.,  6.,  7.,  5.,  5.]),\n",
       " array([6.52063160e-08, 9.99888324e-03, 1.99977013e-02, 2.99965193e-02,\n",
       "        3.99953373e-02, 4.99941554e-02, 5.99929734e-02, 6.99917914e-02,\n",
       "        7.99906094e-02, 8.99894275e-02, 9.99882455e-02, 1.09987064e-01,\n",
       "        1.19985882e-01, 1.29984700e-01, 1.39983518e-01, 1.49982336e-01,\n",
       "        1.59981154e-01, 1.69979972e-01, 1.79978790e-01, 1.89977608e-01,\n",
       "        1.99976426e-01, 2.09975244e-01, 2.19974062e-01, 2.29972880e-01,\n",
       "        2.39971698e-01, 2.49970516e-01, 2.59969334e-01, 2.69968152e-01,\n",
       "        2.79966970e-01, 2.89965788e-01, 2.99964606e-01, 3.09963424e-01,\n",
       "        3.19962242e-01, 3.29961060e-01, 3.39959878e-01, 3.49958696e-01,\n",
       "        3.59957514e-01, 3.69956332e-01, 3.79955150e-01, 3.89953968e-01,\n",
       "        3.99952786e-01, 4.09951604e-01, 4.19950422e-01, 4.29949240e-01,\n",
       "        4.39948059e-01, 4.49946877e-01, 4.59945695e-01, 4.69944513e-01,\n",
       "        4.79943331e-01, 4.89942149e-01, 4.99940967e-01, 5.09939785e-01,\n",
       "        5.19938603e-01, 5.29937421e-01, 5.39936239e-01, 5.49935057e-01,\n",
       "        5.59933875e-01, 5.69932693e-01, 5.79931511e-01, 5.89930329e-01,\n",
       "        5.99929147e-01, 6.09927965e-01, 6.19926783e-01, 6.29925601e-01,\n",
       "        6.39924419e-01, 6.49923237e-01, 6.59922055e-01, 6.69920873e-01,\n",
       "        6.79919691e-01, 6.89918509e-01, 6.99917327e-01, 7.09916145e-01,\n",
       "        7.19914963e-01, 7.29913781e-01, 7.39912599e-01, 7.49911417e-01,\n",
       "        7.59910235e-01, 7.69909053e-01, 7.79907872e-01, 7.89906690e-01,\n",
       "        7.99905508e-01, 8.09904326e-01, 8.19903144e-01, 8.29901962e-01,\n",
       "        8.39900780e-01, 8.49899598e-01, 8.59898416e-01, 8.69897234e-01,\n",
       "        8.79896052e-01, 8.89894870e-01, 8.99893688e-01, 9.09892506e-01,\n",
       "        9.19891324e-01, 9.29890142e-01, 9.39888960e-01, 9.49887778e-01,\n",
       "        9.59886596e-01, 9.69885414e-01, 9.79884232e-01, 9.89883050e-01,\n",
       "        9.99881868e-01]),\n",
       " <a list of 100 Patch objects>)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAADc1JREFUeJzt3X+MZfVZx/H3p4xYqVgoOxLkh0NTqhKMgUwIDUnVbmOwGJZEQmisbs3GTWusVUxktX/U6D+QaGtNiLop1dXUFsRGNlI1SCHERlZ3AQsstl3pQhcXmGrBX9FC+vjHPa7rdmfvmZn78zvvV7KZc+49997nO/fuZ5/zPeeeTVUhSZp/r5l2AZKk0TDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY1YmOSLbdmypZaWlib5kpI09w4cOPCVqloctt1EA31paYn9+/dP8iUlae4leabPdk65SFIjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIyb6TdGNWNp177Hlw7deO8VKJGk22aFLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIa0SvQk/x8kieTPJHkE0lem+TiJPuSHEpyZ5LTx12sJGl1QwM9yfnAzwLLVXUZcBpwE3Ab8OGqehPwVWDHOAuVJJ1a3ymXBeBbkiwAZwBHgbcBd3f37wGuH315kqS+hgZ6VT0H/DrwLIMgfxk4ALxUVa92mx0Bzh9XkZKk4fpMuZwNbAMuBr4DeB1wTd8XSLIzyf4k+1dWVtZdqCTp1PpMubwd+FJVrVTVK8CngKuBs7opGIALgOdO9uCq2l1Vy1W1vLi4OJKiJUnfqE+gPwtcleSMJAG2AgeBB4Abum22A/eMp0RJUh995tD3MTj4+QjwePeY3cAtwM1JDgHnAHeMsU5J0hALwzeBqvog8METbn4auHLkFUmS1sVvikpSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEb0CPclZSe5O8g9JnkryliRvSHJfki92P88ed7GSpNX17dA/AvxFVX038H3AU8Au4P6qugS4v1uXJE3J0EBP8nrgrcAdAFX1tap6CdgG7Ok22wNcP64iJUnD9enQLwZWgN9L8miSjyZ5HXBuVR3ttnkeOPdkD06yM8n+JPtXVlZGU7Uk6Rv0CfQF4Argt6vqcuA/OGF6paoKqJM9uKp2V9VyVS0vLi5utF5J0ir6BPoR4EhV7evW72YQ8C8kOQ+g+/nieEqUJPUxNNCr6nngy0m+q7tpK3AQ2Ats727bDtwzlgolSb0s9NzufcDHk5wOPA38JIN/DO5KsgN4BrhxPCVKkvroFehV9RiwfJK7to62HEnSevlNUUlqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqRG9Az3JaUkeTfJn3frFSfYlOZTkziSnj69MSdIwa+nQ3w88ddz6bcCHq+pNwFeBHaMsTJK0Nr0CPckFwLXAR7v1AG8D7u422QNcP44CJUn99O3QfxP4ReDr3fo5wEtV9Wq3fgQ4f8S1SZLWYGigJ/kR4MWqOrCeF0iyM8n+JPtXVlbW8xSSpB76dOhXA9clOQx8ksFUy0eAs5IsdNtcADx3sgdX1e6qWq6q5cXFxRGULEk6maGBXlW/VFUXVNUScBPwmar6MeAB4IZus+3APWOrUpI01EbOQ78FuDnJIQZz6neMpiRJ0nosDN/k/1TVg8CD3fLTwJWjL2m4pV33Hls+fOu10yhBkmaO3xSVpEasqUOfRXbrkjRghy5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWrEwrQLGKWlXfceWz5867VTrESSJs8OXZIaYaBLUiMMdElqRFNz6GvlnLukltihS1IjDHRJasSmnnI53mrTL8fffuJ9kjRL7NAlqREGuiQ1YmigJ7kwyQNJDiZ5Msn7u9vfkOS+JF/sfp49/nIlSavpM4f+KvALVfVIkjOBA0nuA94N3F9VtybZBewCbhlfqWvjKYmSNpuhHXpVHa2qR7rlfwOeAs4HtgF7us32ANePq0hJ0nBrOsslyRJwObAPOLeqjnZ3PQ+cu8pjdgI7AS666KL11jnz3COQNG29D4om+VbgT4Cfq6p/Pf6+qiqgTva4qtpdVctVtby4uLihYiVJq+vVoSf5JgZh/vGq+lR38wtJzquqo0nOA14cV5Gz6sRz1CVpmvqc5RLgDuCpqvrQcXftBbZ3y9uBe0ZfniSprz4d+tXAjwOPJ3msu+2XgVuBu5LsAJ4BbhxPiZKkPoYGelX9NZBV7t462nIkSevlN0UlqREGuiQ1wkCXpEZ4+dxNwC89SZuDHbokNcIO/ST8wpCkeWSHLkmNMNAlqRGbYspllAcF1zod0+f/KvVApWaVn9P5YocuSY3YFB36pK3Wxc/7wdaNdmt2e9J42aFLUiPs0OfEuLvbcR1nWO25Zq1bn7V6pPWwQ5ekRtihz4CNdIfj6CztVqX5ZIcuSY0w0CWpEZtuymXeTx0cpVn/XbQ89bPWsc3y72KWa9ts7NAlqRGbrkOfdWs95W9Uj50n476cgh3n9Mza737W6hnGDl2SGmGHLm3QvHVxapcduiQ1wg5dIzHuC5KN6rLFG3nOPo898fZZuEzDJC/hPK69lVnYC5qHS2HboUtSI+zQGzXus1lG+fzj6JRH9fyjtNZap93tjdKp9lxm5f1Zi1mt2Q5dkhphoEtSI5xykU5ho7vWk9w1H9WB4GkdyB7Xa0xy6mq1eiZVgx26JDXCDl2nNKsHf1oy7oPC4zDNz8U4fl+tHIC2Q5ekRtihzzC7481tWvPvrWl5bCeyQ5ekRtihSw0aR1c66U7XPZS121CHnuSaJJ9PcijJrlEVJUlau3UHepLTgNuBHwYuBd6Z5NJRFSZJWpuNTLlcCRyqqqcBknwS2AYcHEVhktrVyhTHrNnIlMv5wJePWz/S3SZJmoKxHxRNshPY2a3+e5LPr/OptgBfGU1Vc8Mxbw6OuXG5bcPj/c4+G20k0J8DLjxu/YLutv+nqnYDuzfwOgAk2V9Vyxt9nnnimDcHx9y+SY13I1MufwdckuTiJKcDNwF7R1OWJGmt1t2hV9WrSX4G+EvgNOBjVfXkyCqTJK3JhubQq+rTwKdHVMswG562mUOOeXNwzO2byHhTVZN4HUnSmHktF0lqxMwF+rDLCST55iR3dvfvS7I0+SpHq8eYb05yMMnnktyfpNcpTLOs72Ujkvxokkoy12dE9Blvkhu79/nJJH806RpHrcfn+qIkDyR5tPtsv2MadY5Sko8leTHJE6vcnyS/1f1OPpfkipEWUFUz84fBwdV/BN4InA78PXDpCdv8NPA73fJNwJ3TrnsCY/5B4Ixu+b2bYczddmcCDwEPA8vTrnvM7/ElwKPA2d36t0+77gmMeTfw3m75UuDwtOsewbjfClwBPLHK/e8A/hwIcBWwb5SvP2sd+rHLCVTV14D/vZzA8bYBe7rlu4GtSTLBGkdt6Jir6oGq+s9u9WEG5/zPsz7vM8CvAbcB/zXJ4sagz3h/Cri9qr4KUFUvTrjGUesz5gK+rVt+PfBPE6xvLKrqIeBfTrHJNuAPauBh4Kwk543q9Wct0PtcTuDYNlX1KvAycM5EqhuPtV5CYQeDf+Hn2dAxd7uiF1ZVCxf96PMevxl4c5LPJnk4yTUTq248+oz5V4B3JTnC4Gy5902mtKka6yVTvB76HEnyLmAZ+P5p1zJOSV4DfAh495RLmaQFBtMuP8BgD+yhJN9bVS9Ntarxeifw+1X1G0neAvxhksuq6uvTLmxezVqH3udyAse2SbLAYFftnydS3Xj0uoRCkrcDHwCuq6r/nlBt4zJszGcClwEPJjnMYK5x7xwfGO3zHh8B9lbVK1X1JeALDAJ+XvUZ8w7gLoCq+hvgtQyu8dKyXn/f12vWAr3P5QT2Atu75RuAz1R3tGFODR1zksuB32UQ5vM+twpDxlxVL1fVlqpaqqolBscNrquq/dMpd8P6fK7/lEF3TpItDKZgnp5kkSPWZ8zPAlsBknwPg0BfmWiVk7cX+InubJergJer6ujInn3aR4VXOQr8BQZHyD/Q3farDP5Cw+BN/2PgEPC3wBunXfMExvxXwAvAY92fvdOuedxjPmHbB5njs1x6vsdhMM10EHgcuGnaNU9gzJcCn2VwBsxjwA9Nu+YRjPkTwFHgFQZ7XTuA9wDvOe59vr37nTw+6s+13xSVpEbM2pSLJGmdDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhrxP+sx5/t5utJxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "res = []\n",
    "for i in range(1000):\n",
    "    x = random()\n",
    "    res.append(x*x)\n",
    "plt.hist(res, bins=100)"
   ]
  },
  {
   "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}