MonteCarloExample.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example Script for Monte Carlo and Importance \n",
    "The purpose of this notebook is to suplement my tutorial on Monte Carlo and importance sampling. Feel free to use and adapt any of these functions for your own numeric simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import random\n",
    "from matplotlib import pyplot as plt\n",
    "from IPython.display import clear_output\n",
    "\n",
    "PI = 3.1415926\n",
    "e = 2.71828"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Baseline Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's define some helper functions for later\n",
    "\n",
    "def get_rand_number(min_value, max_value):\n",
    "    \"\"\"\n",
    "    This functions gets a random number from a uniform distribution between\n",
    "    the two input values [min_value, max_value] inclusively\n",
    "    Args:\n",
    "    - min_value (float)\n",
    "    - max_value (float)\n",
    "    Return:\n",
    "    - Random number between this range (float)\n",
    "    \"\"\"\n",
    "    range = max_value - min_value\n",
    "    choice = random.uniform(0,1)\n",
    "    return min_value + range*choice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets also define the function we want to integrate over f(x)\n",
    "\n",
    "def f_of_x(x):\n",
    "    \"\"\"\n",
    "    This is the main function we want to integrate over.\n",
    "    Args:\n",
    "    - x (float) : input to function; must be in radians\n",
    "    Return:\n",
    "    - output of function f(x) (float)\n",
    "    \"\"\"\n",
    "    return (e**(-1*x))/(1+(x-1)**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Crude Monte Carlo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def crude_monte_carlo(num_samples=5000):\n",
    "    \"\"\"\n",
    "    This function performs the Crude Monte Carlo for our\n",
    "    specific function f(x) on the range x=0 to x=5.\n",
    "    Notice that this bound is sufficient because f(x)\n",
    "    approaches 0 at around PI.\n",
    "    Args:\n",
    "    - num_samples (float) : number of samples\n",
    "    Return:\n",
    "    - Crude Monte Carlo estimation (float)\n",
    "    \n",
    "    \"\"\"\n",
    "    lower_bound = 0\n",
    "    upper_bound = 5\n",
    "    \n",
    "    sum_of_samples = 0\n",
    "    for i in range(num_samples):\n",
    "        x = get_rand_number(lower_bound, upper_bound)\n",
    "        sum_of_samples += f_of_x(x)\n",
    "    \n",
    "    return (upper_bound - lower_bound) * float(sum_of_samples/num_samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_crude_MC_variance(num_samples):\n",
    "    \"\"\"\n",
    "    This function returns the variance fo the Crude Monte Carlo.\n",
    "    Note that the inputed number of samples does not neccissarily\n",
    "    need to correspond to number of samples used in the Monte\n",
    "    Carlo Simulation.\n",
    "    Args:\n",
    "    - num_samples (int)\n",
    "    Return:\n",
    "    - Variance for Crude Monte Carlo approximation of f(x) (float)\n",
    "    \"\"\"\n",
    "    int_max = 5 # this is the max of our integration range\n",
    "    \n",
    "    # get the average of squares\n",
    "    running_total = 0\n",
    "    for i in range(num_samples):\n",
    "        x = get_rand_number(0, int_max)\n",
    "        running_total += f_of_x(x)**2\n",
    "    sum_of_sqs = running_total*int_max / num_samples\n",
    "    \n",
    "    # get square of average\n",
    "    running_total = 0\n",
    "    for i in range(num_samples):\n",
    "        x = get_rand_number(0, int_max)\n",
    "        running_total = f_of_x(x)\n",
    "    sq_ave = (int_max*running_total/num_samples)**2\n",
    "    \n",
    "    return sum_of_sqs - sq_ave"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perform the Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Monte Carlo Approximation of f(x): 0.6994168954395196\n",
      "Variance of Approximation: 0.2663694734689711\n",
      "Error in Approximation: 0.005161099431990931\n"
     ]
    }
   ],
   "source": [
    "# Now we will run a Crude Monte Carlo simulation with 10000 samples\n",
    "# We will also calculate the variance with 10000 samples and the error\n",
    "\n",
    "MC_samples = 10000\n",
    "var_samples = 10000 # number of samples we will use to calculate the variance\n",
    "crude_estimation = crude_monte_carlo(MC_samples)\n",
    "variance = get_crude_MC_variance(var_samples)\n",
    "error = math.sqrt(variance/MC_samples)\n",
    "\n",
    "# display results\n",
    "print(f\"Monte Carlo Approximation of f(x): {crude_estimation}\")\n",
    "print(f\"Variance of Approximation: {variance}\")\n",
    "print(f\"Error in Approximation: {error}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Importance Sampling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determine Optimal Weight Function Template"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAH35JREFUeJzt3Xt4lPWd9/H3dyYnEggBEpJAwklQCCggAU8VT6B4KOqjVl3rah+r7aNW2+7VXbvdbbt2n8N2n7Z2t7qrW+1qt0qttRYVT9RzUSEggoBAQIVwkHAOh5DDfPePDDjGQCbJJHdm5vO65Jq57/nN5DN48eHmN/fcP3N3REQktYSCDiAiIomnchcRSUEqdxGRFKRyFxFJQSp3EZEUpHIXEUlBKndJK2Z2gpm9a2Z1ZnaHmZ1vZk/F+dyFZja+uzOKJILpPHdJJ2b2ILDX3b8V3a4Cbnf3t+N47peAq939im6OKdJlOnKXdDMcWAFgZlOB/vEUe9Rc4BwzK+2ucCKJonKXtGFmLwPnAL8ws33AhcBrMY+fbmbbzaw8uj3RzHab2VgAd68HFgPn93x6kY5RuUvacPdzgTdomYbpC5wIrI55fAFwP/CwmfUBfg38nbt/EPMyq4CJPZdapHNU7pLOCoC6Vvt+CPQHFgKbgXtbPV4XfZ5Ir6Zyl3S2C+gXu8PdG4H/BCYAP/HPn3HQD9jdI+lEukDlLulsGXB87A4zGwr8APgV8BMzy271nHHAez0TT6TzVO6SzuYBZx3eMDOj5aj9QeAmYAvwo5jHs4EpwEs9mlKkE1TukrbcfQmwx8xOie66AygG/j46HfMV4Ctmdmb08dnAq+6+uefTinSMvsQkac3MzgdudffL4hj7DnCTu7/f/clEukblLiKSgjQtIyKSglTuIiIpSOUuIpKCMoL6wYWFhT5ixIigfryISFJavHjxdncvam9cYOU+YsQIqqqqgvrxIiJJycw+jmecpmVERFKQyl1EJAWp3EVEUpDKXUQkBancRURSUFzlbmazzGy1mVWb2V1tPH6jmdWa2dLor68mPqqIiMSr3VMhzSxMy2o0M4EaYJGZzXX3la2G/tbdb++GjCIi0kHxnOc+Dah29/UAZjYHuBRoXe49YtFHO3lj7XYyQkY4+isjZITMyAgbOZlhBuZmMbBvFgNzsyjpn0NOZjiIqCIigYmn3IcCG2O2a4BT2hh3hZlNB9YA33L3ja0HmNktwC0Aw4YN63haYMnHu/iXP62Ne7wZDBuYy5jBfTmprIAzRhcysaw/GWF93CAiqavdS/6a2VXABe7+1ej29cA0d/9GzJhBwD53P2RmXwe+FF1p/qgqKyu9s99QdXciDk2RCJFIy21zxGmKOAcbmtm5v4Gd+xvYsb+BjTsPUF27jzVb66iu3Yc79M3O4Nyxg7liShlfGF1IOGSdyiEi0tPMbLG7V7Y3Lp4j9xqgPGa7jJZV4Y9w9x0xm/8B/FM8ITvLzAgbhEOHp1s+O+1SPjC3zeft2t/AW+t38PqaWp57fytz39vMqMI8vnbWKC6fXEZWho7mRSQ1xHPknkHLVMt5wCZgEfAX7r4iZkypu2+J3r8c+Bt3P/VYr9uVI/dEONTUzIsrPuH+19fx/qa9jCrK4wdfHM9Zx7d7PR4RkcDEe+Te7qGquzcBtwMvAKuAx919hZndbWazo8PuMLMVZvYeLetQ3tj56D0jOyPMFycO4enbv8CDN1TiDjc8tJBv/3YpdfWNQccTEemSwJbZC/rIvbVDTc3c98o6/vXltZQNyOXfvzyFiiH5QccSEfmMhB25p4vsjDDfmnk8v/3aaTQ0Rbjq3xfwyuptQccSEekUlXsrU0cM5KnbzmBEYR5ffbiKPy7dFHQkEZEOU7m3oaR/Do9/7TSmjhjAtx9/j2eXbQk6kohIh6jcjyIvO4MHb5jKycMKuHPOu7y+pjboSCIicVO5H0NedgYP3TiV0YP7ctujS6jeVhd0JBGRuKjc29EvJ5MHb5xKdkaImx6uYs8BnSYpIr2fyj0OQwv6cP/1lWzadZC7nlxGUKePiojES+UepynDB/CdC07gufe3MmfR566JJiLSq6jcO+DmM0dx5phC/uHpFXy4fX/QcUREjkrl3gGhkPH/r5pIZjjE3z65XNMzItJrqdw7qDg/h+9eOI631u/gd1U1QccREWmTyr0TrplazrQRA/nHZ1eyY9+hoOOIiHyOyr0TQiHjf18+gf0NzdwzP/5VoUREeorKvZPGFPfjulOG8ejCDaz9RF9uEpHeReXeBd+ccTy5WWH+8dlVQUcREfkMlXsXDMzL4o5zx/DamlreWrej/SeIiPQQlXsXXX/acIrzs/nZ/DU6NVJEeg2VexflZIa57ZzRLPxwJwt09C4ivYTKPQGunlpOaf8cfvLiah29i0ivoHJPgOyMMLeeM5olG3bzzoc7g44jIqJyT5SrppRR2DeL+19bF3QUERGVe6LkZIa58fQRvLK6lg+27g06joikOZV7An351OHkZoW5/7X1QUcRkTSnck+ggtwsrpk6jKff28y2uvqg44hIGlO5J9j1pw2nKeLMWagFPUQkOCr3BBtZmMf044t49J0NNDVHgo4jImlK5d4Nrj91OFv31jN/1SdBRxGRNKVy7wbnjh3M0II+PPLWx0FHEZE0pXLvBuGQ8RenDGPBuh1Ub9PlgEWk56ncu8k1U8vJCof4tY7eRSQAKvduMqhvNhedWMKT726ivrE56DgikmbiKnczm2Vmq82s2szuOsa4K83MzawycRGT11WV5dTVN/HSSn2wKiI9q91yN7MwcC9wIVABXGtmFW2M6wfcAbyT6JDJ6rRRgxjSP4cnFtcEHUVE0kw8R+7TgGp3X+/uDcAc4NI2xv0I+DGgr2ZGhULGFVPKeGNtLVv36LdFRHpOPOU+FIj9umVNdN8RZjYZKHf3Z471QmZ2i5lVmVlVbW1th8MmoytOLiPi8OS7OnoXkZ4TT7lbG/uOrEhhZiHgZ8BftfdC7v6Au1e6e2VRUVH8KZPYiMI8po4YwBOLa7SQh4j0mHjKvQYoj9kuAzbHbPcDJgCvmtlHwKnAXH2o+qkrp5SxvnY/727cHXQUEUkT8ZT7ImCMmY00syzgGmDu4QfdfY+7F7r7CHcfAbwNzHb3qm5JnIQuOrGUnMyQPlgVkR7Tbrm7exNwO/ACsAp43N1XmNndZja7uwOmgn45mVwwvoRnl22hoUkXExOR7pcRzyB3nwfMa7Xv+0cZe3bXY6We2ROH8Melm3mzupZzxxYHHUdEUpy+odpDzhxTRP8+mTz93pago4hIGlC595CsjBAXTijhxRVbOdigyxGISPdSufegL04cwv6GZl5ZvS3oKCKS4lTuPejUUYMo7JvN3KWb2x8sItIFKvceFA4Zl5xUysurt1FX3xh0HBFJYSr3HvbFiUNoaIrw4gpdKVJEuo/KvYedPKyAoQV9eGaZpmZEpPuo3HuYmXHhhBLerN7OXk3NiEg3UbkH4MITS2hsdl75QGfNiEj3ULkHYHL5AIrzs3lu+dago4hIilK5ByAUMi4YX8Kra7ZxoKEp6DgikoJU7gGZNb6E+sYIr61Oj0VLRKRnqdwDMm3kQAbkZvL8Ck3NiEjiqdwDkhEOcX5FCS+v2sahJl1rRkQSS+UeoFknllB3qIk/V28POoqIpBiVe4BOP24Q/bIzeP59Tc2ISGKp3AOUnRHmvHGDeWnlJzQ1a4UmEUkclXvAZk0oZdeBRt75cGfQUUQkhajcA3bW8UX0yQxrakZEEkrlHrA+WWGmH1/I/FWf4O5BxxGRFKFy7wVmjCtmy5563t+0N+goIpIiVO69wHnjigkZvLRSUzMikhgq915gYF4WlcMH8uJKLeAhIomhcu8lZlYU88HWOjbuPBB0FBFJASr3XmJmRTEA81fp6F1Euk7l3kuMKMxjzOC+vKSpGRFJAJV7LzKzoph3PtzJngNafk9Eukbl3ovMrCimOeK8slrL74lI16jce5GJZQUM7petqRkR6TKVey8SChnnjSvm1dW6xruIdE1c5W5ms8xstZlVm9ldbTz+dTNbbmZLzexNM6tIfNT0cH5FMfsbmnlr3Y6go4hIEmu33M0sDNwLXAhUANe2Ud6PuvuJ7j4J+DHw04QnTROnHTeI3KywpmZEpEviOXKfBlS7+3p3bwDmAJfGDnD32Iui5AG6AlYn5WSGmT6miPmrPiES0W+jiHROPOU+FNgYs10T3fcZZnabma2j5cj9jrZeyMxuMbMqM6uqra3tTN60MLOimE/2HmL5pj1BRxGRJBVPuVsb+z53SOnu97r7ccDfAH/X1gu5+wPuXunulUVFRR1LmkbOHTuYcMg0NSMinRZPudcA5THbZcDmY4yfA1zWlVDpbkBeFpXDB+hSBCLSafGU+yJgjJmNNLMs4BpgbuwAMxsTs3kxsDZxEdOTLiQmIl3Rbrm7exNwO/ACsAp43N1XmNndZjY7Oux2M1thZkuBbwM3dFviNHF+RQmALgMsIp2SEc8gd58HzGu17/sx9+9McK60N2xQLicU9+OllVu56Qsjg44jIklG31DtxWZWFLPoo13sPtAQdBQRSTIq917s8IXEXv5AFxITkY5RufdiJw7tT3G+LiQmIh2ncu/FQiFjxrhiXltTS32jLiQmIvFTufdyMyqKOaALiYlIB6nce7nTjxtE3+wMXlixNegoIpJEVO69XHZGmLNPKOKllZ/QrAuJiUicVO5J4ILxJezY38CSDbuCjiIiSULlngTOPqGIrHCIF97X1IyIxEflngT65WRyxuhBvLByK+6amhGR9qnck8QF40vYuPMgq7bUBR1FRJKAyj1JzKgoJmTorBkRiYvKPUkU9s2mcvhAlbuIxEXlnkTOH99yjfcNO3SNdxE5NpV7ErlgfMs13nX0LiLtUbknkfKBuYwrzVe5i0i7VO5J5oLxxSzesIvaukNBRxGRXkzlnmQuGF+CO1o8W0SOSeWeZMaW9GPYwFxNzYjIManck4yZccH4YhZU76CuvjHoOCLSS6nck9CsCSU0NEf40yotvycibVO5J6HJ5QMo7Z/DM8u2BB1FRHoplXsSCoWMi04s5fU1tezV1IyItEHlnqQuPqmUhuYI87V4toi0QeWepCaXFzC0oA/PampGRNqgck9SZsZFJ5bw+tpa9hzQ1IyIfJbKPYldfNIQGpudF1fqnHcR+SyVexKbWNa/ZWpmuaZmROSzVO5JzMy45KRS3ly7nd0HGoKOIyK9iMo9yV1y0hCaIs6LK3TWjIh8SuWe5CYMzWfYwFye0dSMiMSIq9zNbJaZrTazajO7q43Hv21mK81smZn9ycyGJz6qtMXMuPikUv5cvZ0d+3QZYBFp0W65m1kYuBe4EKgArjWzilbD3gUq3f0k4Angx4kOKkc3e+IQmiOuD1ZF5Ih4jtynAdXuvt7dG4A5wKWxA9z9FXc/vLDn20BZYmPKsYwrzWdsST/+8O6moKOISC8RT7kPBTbGbNdE9x3NTcBzbT1gZreYWZWZVdXW1safUtp1+eShvLthNx9u3x90FBHpBeIpd2tjn7c50OzLQCXwz2097u4PuHulu1cWFRXFn1LademkoZjBUzp6FxHiK/caoDxmuwzY3HqQmc0AvgfMdnd9stfDSvrncPpxg3hq6Sbc2/y7V0TSSDzlvggYY2YjzSwLuAaYGzvAzCYD99NS7FpBIiCXTRrKxzsOsGTD7qCjiEjA2i13d28CbgdeAFYBj7v7CjO728xmR4f9M9AX+J2ZLTWzuUd5OelGsyaUkJMZ0tSMiJARzyB3nwfMa7Xv+zH3ZyQ4l3RCv5xMZlaU8PSyzfz9JRVkZeg7aiLpSn/6U8yVU8rYfaCR+at0OQKRdKZyTzFfGF3IkP45zFm0sf3BIpKyVO4pJhwyrqws5421tWzafTDoOCISEJV7CrpqSssXhJ+oqgk4iYgEReWegsoH5nLGcYU8XrWRSETnvIukI5V7ivrS1HI27T7In9dtDzqKiARA5Z6izq8opiA3k9/qg1WRtKRyT1E5mWEunzyUF1ZspbZOV4MQSTcq9xR23SnDaWx2frtoQ9BRRKSHqdxT2OjBfTlzTCH/9fYGmpojQccRkR6kck9x1586nK1763lppb6xKpJOVO4p7rxxxQwt6MMjb30cdBQR6UEq9xQXDhlfPnU4b63fwZpP6oKOIyI9ROWeBq6eWk5WRoiHF3wUdBQR6SEq9zQwMC+L2ROH8OSSTeza3xB0HBHpASr3NHHzmaM42NjMr9/W3LtIOlC5p4kTSvpx7tjB/OeCj6hvbA46joh0M5V7Gvna9FHs3N/A7xbrapEiqU7lnkamjRzIpPIC/uP19TTrapEiKU3lnkbMjK+fNYoNOw/w/Ptbg44jIt1I5Z5mZlaUMLIwj3tfqcZdR+8iqUrlnmbCIeO2c0azcsteXlihSxKIpCqVexq6bNIQRhbmcc/8NVqpSSRFqdzTUEY4xJ3njeGDrXU8p7l3kZSkck9TX5w4hNGD+3LP/DU6c0YkBanc01Q4ZHxzxhjWbtvHM8s2Bx1HRBJM5Z7GLppQytiSfvzkxTUcatK3VkVSico9jYVCxncvGseGnQd4ZIGuOSOSSlTuae6s44s4+4Qi/uXltezUFSNFUobKXfjeReM40NDMPfPXBB1FRBIkrnI3s1lmttrMqs3srjYen25mS8ysycyuTHxM6U5jivtx7bRyfvPOBqq3abUmkVTQbrmbWRi4F7gQqACuNbOKVsM2ADcCjyY6oPSMb804ntzMMD+cu1KXJRBJAfEcuU8Dqt19vbs3AHOAS2MHuPtH7r4MiHRDRukBg/pm851ZJ/Bm9XaeWrop6Dgi0kXxlPtQYGPMdk10X4eZ2S1mVmVmVbW1tZ15CelG150ynEnlBfzomVVajk8kycVT7tbGvk79u93dH3D3SnevLCoq6sxLSDcKh4z/+z9OZO/BRv7PvFVBxxGRLoin3GuA8pjtMkBfaUxR40rzuXn6KH63uIYF67YHHUdEOimecl8EjDGzkWaWBVwDzO3eWBKkO88bw4hBuXznd8vYW98YdBwR6YR2y93dm4DbgReAVcDj7r7CzO42s9kAZjbVzGqAq4D7zWxFd4aW7pWTGeZnV09i6956fvBH/a8USUYZ8Qxy93nAvFb7vh9zfxEt0zWSIiYPG8Ad547hZ/PXcM7YwcyeOCToSCLSAfqGqhzVbeccx+RhBXzvD8vZtPtg0HFEpANU7nJUGeEQ91w9iUjEufU3S3TlSJEkonKXYxo+KI+ffGki723czQ/nav5dJFmo3KVdsyaUcuvZx/HYwo08tnBD0HFEJA4qd4nLX51/AmeOKeQHf1zB4o93Bh1HRNqhcpe4hEPGv147mSEFOdz0cBXravcFHUlEjkHlLnEryM3i4f85jYyQccNDC9m2tz7oSCJyFCp36ZDhg/J46Map7NzfwI2/WqRvsIr0Uip36bCTygq477qTWfNJHX/54EIVvEgvpHKXTjn7hMHce93JvL9pjwpepBdSuUunXTC+hPuuO5kVm/dw/YML2X1A14AX6S1U7tIl548v4b7rprBq816u+LcFbNx5IOhIIoLKXRJgZkUxv75pGrV1h7j8vgUsr9kTdCSRtKdyl4Q4ZdQgnrz1dLIzQlz9wFs8/Z7WcxEJkspdEmb04H784dbTGVeazzcee5e7n15JY7PWTBcJgspdEmpwfg6P3XwqN54+gof+/CHXPvC25uFFAqByl4TLygjxw9nj+fk1k/hgax0X/vwNHq/aiHun1lUXkU5QuUu3uXTSUJ6780zGD8nnr59Yxs2PVFGzS0fxIj1B5S7dqnxgLo/dfCp/d/E4/ly9gxk/fY37Xq2moUlz8SLdSeUu3S4UMr565ihe+vZ0po8p4sfPr2bWPa/z3PItmqoR6SYqd+kxZQNyeeAvK3noxkpCIeN//WYJl923gDfXblfJiySYBfWHqrKy0quqqgL52RK8puYITy7ZxE9fWsPWvfWcVNafr591HBeMLyEcsqDjifRaZrbY3SvbHadylyDVNzbz+yU1PPD6ej7ecYCRhXl85YwRXDppKP37ZAYdT6TXUblLUmmOOM+/v5X7X1/Hspo95GSGuPjEIVw7rZwpwwdgpqN5EVC5SxJbXrOHxxZtYO7Szew71MTwQblcfGIpF59USkVpvope0prKXZLe/kNNPLtsC08v28yCdTtojjgjC/OYMW4w048vYuqIgeRkhoOOKdKjVO6SUnbub+CFFVuZt3wL76zfSUNzhJzMEKeOGsT0MUWcMmogY0vy9WGspDyVu6SsAw1NvLN+J6+tqeW1NbV8uH0/AHlZYSYPG8DJwwcwZfgAxg/Jp7BvdsBpRRIr3nLP6IkwIomUm5XBOWMHc87YwQDU7DrA4o93UfXRLhZ/vItfvLyWSPSYpahfNhWl+YwrzWdcaT/GleYzbGCupnMk5ancJemVDcilbEAul04aCkBdfSPLa/awcsteVm2pY9WWvSxYt57G5pbGN4Mh/fswqiiPEYPyGFGYx/CBuZQW5DCkfx8KcjP1oa0kPZW7pJx+OZmcPrqQ00cXHtnX0BRhXe0+1nxSx4fb9/PR9v18uH0/Ty3dRF1902ee3yczfKToS/vnUJyfw6C+WQzqm01hXsvtoL5ZDMjN0hy/9FpxlbuZzQJ+DoSBX7r7/2v1eDbwCDAF2AFc7e4fJTaqSOdlZYSiUzP5n9nv7uzc38DGXQfZsvsgm/fUs3n3QbbsOcjm3fW8vraW7fsaaI58/rMpMxiQm0VBn0z69ckkPyeD/JxM8vu03PbLySC/T+aR+32ywvTJDJOblUFuVvjIdp/MMCH9JSEJ1m65m1kYuBeYCdQAi8xsrruvjBl2E7DL3Ueb2TXAPwFXd0dgkUQys+iReDaTygvaHBOJOHsONrJj/yG272tgx76GmPuH2HOwkb31TdTVN7J590H21jex92Ajhzpw5cuczNCR4u+TFSYnM0RWOERmOERWRojsjJbbzHDL/sP3D+/PCofIPHwbNkIhIyNkhEMhwiEIh0LRbSNsRjgcfdxa9mWEo2M/s/3p42Ytv1chg5AZxqfbsbcWfTxkYHy6/Zn9mvLqEfEcuU8Dqt19PYCZzQEuBWLL/VLgh9H7TwC/MDNzXQ1KUkAoZAzIy2JAXhajB8f/vENNzdRFi35vfRMHG5o52NjEgYbm6P3mVvdbHqtvbKa+MUJDU4SG5gh19U3siN5vbI7uj24fvk22P2kW85fE4fJvvY+W/6Lj7cjzoI39Ma8b81OO8pzD28d+zdZ5O/LcNmJ8Zswd541h9sQhn/s5iRRPuQ8FNsZs1wCnHG2MuzeZ2R5gELA9dpCZ3QLcAjBs2LBORhZJDtkZYbL7hrv9dEx3pyniR4q/KeI0x/xq2Y7QHIGmSKSNx9rejh3r0Z/jDhGHiH+6LxJ9POKxY7xlm+j2kTEt+4neHn6dSPR5HrvfP/8+AfzIdvQ2uid2fOsxtBpztOfG/sjWY1rdHCPP0TMfvlPQA9dNiqfc2/o3VOvjhHjG4O4PAA9Ay3nucfxsEWmHmZEZNjLDIXKzgk4jvUU813OvAcpjtsuAzUcbY2YZQH9gZyICiohIx8VT7ouAMWY20syygGuAua3GzAVuiN6/EnhZ8+0iIsFpd1omOod+O/ACLadCPuTuK8zsbqDK3ecCDwK/NrNqWo7Yr+nO0CIicmxxnefu7vOAea32fT/mfj1wVWKjiYhIZ2kNVRGRFKRyFxFJQSp3EZEUpHIXEUlBgS3WYWa1wMedfHohrb79moT0HoKX7Pkh+d9DsueHnn8Pw929qL1BgZV7V5hZVTwrkfRmeg/BS/b8kPzvIdnzQ+99D5qWERFJQSp3EZEUlKzl/kDQARJA7yF4yZ4fkv89JHt+6KXvISnn3EVE5NiS9chdRESOQeUuIpKCkq7czWyWma02s2ozuyvoPB1lZg+Z2TYzez/oLJ1hZuVm9oqZrTKzFWZ2Z9CZOsrMcsxsoZm9F30P/xB0ps4ws7CZvWtmzwSdpTPM7CMzW25mS82sKug8nWFmBWb2hJl9EP0zcVrQmQ5Lqjn36GLda4hZrBu4ttVi3b2amU0H9gGPuPuEoPN0lJmVAqXuvsTM+gGLgcuS7P+BAXnuvs/MMoE3gTvd/e2Ao3WImX0bqATy3f2SoPN0lJl9BFS6e9J+icnMHgbecPdfRte7yHX33UHnguQ7cj+yWLe7NwCHF+tOGu7+Okm8SpW7b3H3JdH7dcAqWtbQTRreYl90MzP6K3mOcgAzKwMuBn4ZdJZ0ZWb5wHRa1rPA3Rt6S7FD8pV7W4t1J1WxpBIzGwFMBt4JNknHRac0lgLbgJfcPdnewz3AXwORoIN0gQMvmtliM7sl6DCdMAqoBX4VnR77pZnlBR3qsGQr97gW4pbuZ2Z9gd8D33T3vUHn6Sh3b3b3SbSsCTzNzJJmiszMLgG2ufvioLN00RnufjJwIXBbdMoymWQAJwP/5u6Tgf1Ar/kcMNnKPZ7FuqWbReepfw/8xt2fDDpPV0T/Gf0qMCvgKB1xBjA7Omc9BzjXzP4r2Egd5+6bo7fbgD/QMu2aTGqAmph/9T1BS9n3CslW7vEs1i3dKPph5IPAKnf/adB5OsPMisysIHq/DzAD+CDYVPFz9++6e5m7j6Dlz8DL7v7lgGN1iJnlRT+QJzqVcT6QVGeQuftWYKOZnRDddR7Qa04siGsN1d7iaIt1BxyrQ8zsMeBsoNDMaoAfuPuDwabqkDOA64Hl0TlrgL+NrrObLEqBh6NnX4WAx909KU8nTGLFwB9ajhXIAB519+eDjdQp3wB+Ez3YXA98JeA8RyTVqZAiIhKfZJuWERGROKjcRURSkMpdRCQFqdxFRFKQyl1EJAWp3EVEUpDKXUQkBf034D97aW5uFTIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x210e46c85f8>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the function\n",
    "xs = [float(i/50) for i in range(int(50*PI*2))]\n",
    "ys = [f_of_x(x) for x in xs]\n",
    "plt.plot(xs,ys)\n",
    "plt.title(\"f(x)\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is the template of our weight function g(x)\n",
    "def g_of_x(x, A, lamda):\n",
    "    e = 2.71828\n",
    "    return A*math.pow(e, -1*lamda*x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl8VPW9//HXJ/vOloQl7AIKIghExFqtiChat6q1arV1qf5sa+1yb2vrba21t9f2dq/aW7UudV9bRStad+uGBBCQPSxqCJAEsu+Z+f7+OAPGkGWASc7M5P18POYxc875zszny+h7Tr5zzvmacw4REYkvCX4XICIikadwFxGJQwp3EZE4pHAXEYlDCncRkTikcBcRiUMKd/GdmR1qZsvNrNbMrg2tO9nMngrz+e+Z2eG9W+Wn3s+Z2YSDeP5bZjYjjHZnmtkjB/o+0r8p3CUa/AB4zTmX7Zz7U2jd/wC/DPP5vwFu6pXKIszMzgBqnXPLe2rrnFsITDWzab1fmcQbhbtEgzHA6j0LZnYUMMA5926Yz18IzDWz4b1RXIRdDdy/H+0fBq7qpVokjincxVdm9gowF7jVzOrMbBJwKvB6uzafMbMKMxsVWp5uZlVmdhiAc64JWAqc3MV7HGJmr5jZrtDrPGhmA9tt32pm/2lmK82s2sweNbO0dtu/b2bbzazUzC7voT/jzOyN0BDTS2Z2m5k9ENqWApzYoW/Pmdlv2y0/amZ3t3vJ14DP9/DPKLIPhbv4yjl3IvBv4BrnXJZzbgNwBLC+XZu3gduBv5lZOt6e74+dc+vavdRaYHoXb2PAzcAIYDIwCrixQ5vzgQXAOGAacCmAmS0A/hOYD0wETuqhSw8B7wFDQu9xSbttE4Ggc66k3brLgUvM7EQz+zJwFPDtDv0aa2Y5PbyvyKck+V2ASCcGArUd1t0IvIsXnNuA2zpsrwU6HZZxzhUDxaHFcjP7HfDTDs3+5JwrBTCzZ4AjQ+vPB+5xzn0Q2nYjcGFn72Nmo/HCeZ5zrgV408wWdtcv59wOM7sa+BuQDpztnGvfprbdc2s6e1+RzmjPXaJRJZDdfoVzrhW4F5gK/Nbte8W7bKCqsxczs3wze8TMtplZDfAAkNuh2Y52jxuArNDjEcDH7bZ92E3dI4DdzrmGduvaP3effoU8CyQC651zb3bYtqd9p30T6YrCXaLRSmBS+xVmVoC3t30P8FszS+3wnMnAii5e72bAAdOccznAxXhDNeHYjjeMs8foHtoONrOMduvaP3cjYKG+tPcLvOGX4WbW8a+CycBW55z22mW/KNwlGj0HfG7PgpkZ3l77XcAVeCH683bbU4FZwItdvF42UAdUhYL1+/tRy2PApWY2JRTaHYdz9nLOfQgUATeaWYqZHQOc0W57K/BSh74dD1wGfCV0u6VD+H8OWLQf9YoACneJQs65ZUC1mR0dWnUtMBT4SWg45jLgMjM7LrT9TLzj5Eu7eMmfATOBauCfwN/3o5ZFwB+AV/DG7V/p4SlfBo4BdgH/DTwKNLfbfjuhH1lDP5Leh/dj8rbQkMxdwD2hLzTwxvdvD7dekT1Mk3VINDKzk4FvOOfODqPtYuCKPT96RhMzexRY55z7abt1bwLf6ulEptAJT5c4587v5TIlDincRSIodALWbmAL3nH3TwHHhHNGqkgk6VBIkcgahjfsMwQoAb6uYBc/aM9dRCQO6QdVEZE45NuwTG5urhs7dqxfby8iEpOWLl1a4ZzL66mdb+E+duxYioqK/Hp7EZGYZGbdnSW9l4ZlRETikMJdRCQOKdxFROKQwl1EJA4p3EVE4lCP4W5md5tZmZl1e90OMzvKzAJmdl7kyhMRkQMRzp77vXjTj3XJzBKBXwEvRKAmERE5SD2Gu3PuDbwLIXXnW8CTQFkkiupW+Xp4/kfQ1tLrbyUiEqsOesw9NLHAF4C/hNH2KjMrMrOi8vLyA3vDyg/h3T9DcVfzMoiISCR+UP0DcJ1zLtBTQ+fcHc65QudcYV5ej2fPdu6QuZCRCysfPbDni4j0A5G4/EAh8Eho4phc4DQza3POPRWB195XYjJMPReW3guNVZA+sFfeRkQklh30nrtzbpxzbqxzbizwBN7sOb0T7HtM+xIEmmHtwl59GxGRWBXOoZAPA+8Ah5pZiZldYWZXm9nVvV9eFwpmwpAJsEJDMyIinelxWMY5d2G4L+acu/SgqgmXmbf3/uovoOpjGDiqT95WRCRWxO4Zqkd80btf9bi/dYiIRKHYDffB42DUHO+oGU0VKCLyKbEb7gDTzofydbBjpd+ViIhEldgO98O/AAnJsPIxvysREYkqsR3uGYNh0ineuHugze9qRESiRmyHO3hHzdTthM2v+l2JiEjUiP1wn7QAMobA8vv9rkREJGrEfrgnpXh77+ueg/pdflcjIhIVYj/cAY78MgRbdcy7iEhIfIT7sKkw/EhY/oDflYiIRIX4CHeAGRfDzlWwfYXflYiI+C5+wn3quZCYqr13ERHiKdwzBsPk070Tmlqb/K5GRMRX8RPu4P2w2lQF65/zuxIREV/FV7iPPwFyRsL7D/pdiYiIr+Ir3BMS4ciLoPhlqPrI72pERHwTX+EOMPMS737Zff7WISLio/gL94GjYeJ8WHY/BFr9rkZExBfxF+4Asy6Duh2w4Xm/KxER8UU4E2TfbWZlZvZBF9u/bGYrQ7e3zWx65MvcTxNPhuwRUHSP35WIiPginD33e4EF3WzfAnzOOTcN+DlwRwTqOjiJSTDzK7DpFdi9xe9qRET6XI/h7px7A9jdzfa3nXOVocV3gZERqu3gzPwKmMGyv/ldiYhIn4v0mPsVwKIIv+aBGVDgXet9+QPQ1uJ3NSIifSpi4W5mc/HC/bpu2lxlZkVmVlReXh6pt+7arMugvhzW/7P330tEJIpEJNzNbBrwV+As51yXM2Y45+5wzhU65wrz8vIi8dbdmzAPBoyCort7/71ERKLIQYe7mY0G/g5c4pzbcPAlRVBCIsy6FLa8AeXr/a5GRKTPhHMo5MPAO8ChZlZiZleY2dVmdnWoyQ3AEODPZva+mRX1Yr37b9al3qWA3/P/IB4Rkb6S1FMD59yFPWz/GvC1iFUUaZm53rXe338Y5t0AaQP8rkhEpNfF5xmqHR19FbTWw3JdLVJE+of+Ee4jZsDI2bDkTggG/a5GRKTX9Y9wBzj6/8HuzVD8kt+ViIj0uv4T7pPPhKxh8N7tflciItLr+k+4J6VA4eXennvFRr+rERHpVf0n3ME7LDIhGd670+9KRER6Vf8K9+yhMPUcb47Vxiq/qxER6TX9K9wBjrkGWupg6b1+VyIi0mv6X7gPnwbjPgeL/6KrRYpI3Op/4Q7wmW9B7XZY/Xe/KxER6RX9M9wnnAR5h8Hbt4BzflcjIhJx/TPczbyx950fwObX/K5GRCTi+me4A0w7HzLz4Z1b/a5ERCTi+m+4J6V6FxQrfgl2rvG7GhGRiOq/4Q5QeAUkpWvvXUTiTv8O94zBMOursPJRqPrY72pERCKmf4c7eIdFgnfkjIhInFC4DxgJ0y+AZX+DunK/qxERiQiFO8Cx34W2Znj3z35XIiISEQp3gNwJMOUsWPJXXVBMROJCj+FuZnebWZmZfdDFdjOzP5lZsZmtNLOZkS+zDxz3PWiu8QJeRCTGhbPnfi+woJvtpwITQ7ergP87+LJ8MHw6TJjvDc20NPhdjYjIQekx3J1zbwC7u2lyFnCf87wLDDSz4ZEqsE8d9x/QsMv7cVVEJIZFYsy9AGh/kHhJaN0+zOwqMysys6Ly8ig8MmXMMTDms/Dm76G10e9qREQOWCTC3TpZ1+mlFp1zdzjnCp1zhXl5eRF4615wwg+hbicU3eN3JSIiBywS4V4CjGq3PBIojcDr+mPccTD2OG/vXWPvIhKjIhHuC4GvhI6amQNUO+e2R+B1/TP3eqgvg6K7/a5EROSAhHMo5MPAO8ChZlZiZleY2dVmdnWoyXPAZqAYuBP4Rq9V21fGfMabiu+tP0BLvd/ViIjst6SeGjjnLuxhuwO+GbGKosXc6+HuU2DJXXDstX5XIyKyX3SGaldGz4Hxc+GtP0Jznd/ViIjsF4V7d+ZeDw0V8N4dflciIrJfFO7dGTUbJp7ijb03VvpdjYhI2BTuPZl3AzTVeIdGiojECIV7T4ZN9SbTXnw71MTu4fsi0r8o3MMx93oIBuC1X/pdiYhIWBTu4Rg0Fgovh+UPQMVGv6sREemRwj1cx38fktLglZ/7XYmISI8U7uHKyoPPXANrnoZtS/2uRkSkWwr3/XHMNZCRCy/8GFynF74UEYkKCvf9kZYDJ/4XfPQ2rF3odzUiIl1SuO+vGV+BvMnw4g3Q1ux3NSIinVK476/EJDjlF1C51Tv2XUQkCincD8SEed5k2m/8Guor/K5GRGQfCvcDdfJ/e9d6f+1mvysREdmHwv1A5R8GhZd5c62WrfW7GhGRT1G4H4wTfgSpWbDoBzo0UkSiisL9YGTmwok/gS1vwOp/+F2NiMheCveDVXg5DDsC/vVjzdgkIlEjrHA3swVmtt7Mis3sh51sH21mr5rZcjNbaWanRb7UKJWQCKf9Bmq2wb9/43c1IiJAGOFuZonAbcCpwBTgQjOb0qHZj4HHnHMzgAuAP0e60Kg2eg5MvxDevhUqiv2uRkQkrD332UCxc26zc64FeAQ4q0MbB+SEHg8A+t+sFvNvguR0/bgqIlEhnHAvAD5ut1wSWtfejcDFZlYCPAd8q7MXMrOrzKzIzIrKy8sPoNwolpXvTeqx6WVY85Tf1YhIPxdOuFsn6zruml4I3OucGwmcBtxvZvu8tnPuDudcoXOuMC8vb/+rjXZHXQnDp8Oi66Cxyu9qRKQfCyfcS4BR7ZZHsu+wyxXAYwDOuXeANCA3EgXGlMQkOOOPUF8OL//M72pEpB8LJ9yXABPNbJyZpeD9YNrxercfAfMAzGwyXrjH2bhLmEbMgKO/DkV3w0fv+l2NiPRTPYa7c64NuAZ4AViLd1TMajO7yczODDX7D+BKM1sBPAxc6lw//lVx7vUwYBQ8821oa/G7GhHph8yvDC4sLHRFRUW+vHef2PACPHQ+nPhjb/5VEZEIMLOlzrnCntrpDNXeMukUOPwL8PqvoXyD39WISD+jcO9NC34FKRnw9DcgGPC7GhHpRxTuvSl7qHdpgpIl8M6tflcjIv2Iwr23TT0XDjsdXvkFlK/3uxoR6ScU7r3NDD7/O2945ikNz4hI31C494U9wzPbiuDtW/yuRkT6AYV7X5l6Lkw+A179Bexc7Xc1IhLnFO59xQxO/wOkDYQnvwatTX5XJCJxTOHelzJz4azboGwNvHyT39WISBxTuPe1SSd7V4989zbY9Irf1YhInFK4+2H+TZB7qHf0TMNuv6sRkTikcPdDSgaceyfUV8Az12rmJhGJOIW7X4ZPh3k/gbXPwJK/+l2NiMQZhbufjvkWTJgPL1wP21f4XY2IxBGFu58SEuALt0NGLjx+KTTV+F2RiMQJhbvfMofAeXdB5Yfw7Hc0/i4iEaFwjwZjPuPN3vTBk7D0Xr+rEZE4oHCPFp/9HhxyIiy6DrYt87saEYlxCvdokZAA5/wVsobCo5d4h0mKiBwghXs0yRwCX7ofGirgicsg0OZ3RSISo8IKdzNbYGbrzazYzH7YRZvzzWyNma02s4ciW2Y/MuJIOP33sOUNePlnflcjIjEqqacGZpYI3AbMB0qAJWa20Dm3pl2bicCPgGOdc5Vmlt9bBfcLR14E25bC23+CETNg6jl+VyQiMSacPffZQLFzbrNzrgV4BDirQ5srgducc5UAzrmyyJbZD51yM4w62rv+TOn7flcjIjEmnHAvAD5ut1wSWtfeJGCSmb1lZu+a2YLOXsjMrjKzIjMrKi8vP7CK+4ukFPjSA95lgh++EGq2+12RiMSQcMLdOlnX8UybJGAicAJwIfBXMxu4z5Ocu8M5V+icK8zLy9vfWvufrHy48GFoqoZHLoLWRr8rEpEYEU64lwCj2i2PBEo7afO0c67VObcFWI8X9nKwhh3hXUGydLk3RKMzWEUkDOGE+xJgopmNM7MU4AJgYYc2TwFzAcwsF2+YZnMkC+3XDvs8nPRTWP13eO1mv6sRkRjQ49Eyzrk2M7sGeAFIBO52zq02s5uAIufcwtC2k81sDRAAvu+c29Wbhfc7x34HKjbC67+CAaNg5iV+VyQiUcycT3/mFxYWuqKiIl/eO2YFWuGh82Hz63DRYzDxJL8rEpE+ZmZLnXOFPbXTGaqxJDEZzr8Phk6Bx76iQyRFpEsK91iTmg1ffgIyBsODX4TKrX5XJCJRSOEei7KHwcVPQqAF7v8C1O70uyIRiTIK91iVdyh8+XGo3QEPnAONlX5XJCJRROEey0bNhgsehIoN8OD50FLvd0UiEiUU7rHukBPh3LtgWxE8ejG0NftdkYhEAYV7PJhyJpx5C2x6BR6/DNpa/K5IRHzW40lM0WbJ1t3832ubSDAAI8HADBLMMAMzw/hkOSG0bHuXwTASEj55fnJiAilJCaQmJZCSmEBq8p77xL3LqUmJZKYmkpOWTHZaEtmh++TEKPl+nHGxd+2Z5/7Tm+jji/d6h06KSL8Uc+He1BqgrLYJ5yDowDmHc+BwHZYhGHq8594512G9t64t6GhuC9DcFtzvS7ekJiWQnZZMTnoSuZmp5GankJuVuveWn53KiIHpjBycTk5aL4ft7CshGIDnr4Mnr/CGaxTwIv1SzIX7cRPzOG5i71xRck/Qt7QFaW4Lhu4DtLQFaWoNUtfcRm1TK7VN7e5D66obW6moa2Hdjloqaiuoadp3irzstCRGDsqgYGA6owdnMCE/a+9tcGZKZDox52pwAXjherDQvKyJMfcxi8hB0v/17ZgZyYlGcmICmakH91rNbQF217ews6aZbZWNbKtqYFtlIyWVjZRUNvBWcQWNrYG97YdkpnBIfhaHDctm2siBTB85gPF5WSQmdHbF5R4c801vD/7Fn3j3597lXR9eRPoNXVvGJ8Ggo7S6kY1ldWwqq6O4rI6NZXWs215DfYsX+pkpiRxeMIDpIwcwa8wgjh43hEH7s4f/zp/hhR/BxJO9yxYkp/dSb0Skr4R7bRmFe5QJBB2by+tYUVLNqpIqVpRUs2Z7DS1tQcxg8rAcjjlkCMeMH8Ls8YN7Hscvugee/S6MOw4ueBhSs/qmIyLSKxTucaSlLcjKkire2bSLtzftYulHlbS0BUkwmDVmEPOnDOWkyUMZn9dFcK94FJ66GkYeBRc9CumD+rYDIhIxCvc41tQaYPlHVby9qYKX15axZnsNAIfkZTJ/yjDmT8lnxqhBJLQfr1/zNDxxBQyZ4F2XZkDHaXBFJBYo3PuRksoGXl5bxotrdvLu5l20BR0FA9M568gRnDOzgAn52V7Dza/DI1+GtBzvypJDp/hbuIjsN4V7P1XT1Mora8t46v1t/HtjBYGgY2pBDmcfWcCZ00eQ37ARHjgP2hq9Mfixx/pdsojsB4W7UF7bzLMrS/nH8m2sLKkmMcE4aXI+X5uaTOFbX8Mqt8I5d8LhZ/tdqoiESeEun1JcVscTS0t4rOhjdte3MH1IgDuSfkN+9UpswS+9k59EJOpFdJo9M1tgZuvNrNjMfthNu/PMzJlZj28sfWtCfhY/PPUw3vnRifzhS0eSlJXL8Tu/y8tuFjx/HVWPX+vN0SoicaHHM1TNLBG4DZgPlABLzGyhc25Nh3bZwLXA4t4oVCIjNSmRs2cUcPaMAtaU1vDgu+PY+v5v+Nrqv7F20woaz76bGYeOx+wAzowVkagRzp77bKDYObfZOdcCPAKc1Um7nwP/CzRFsD7pRVNG5PCLc47knOvu5vlJP2N802pyHzqF797yEC+s3kEw6M+QnYgcvHDCvQD4uN1ySWjdXmY2AxjlnHu2uxcys6vMrMjMisrLy/e7WOkdgzNTWHDRd3Bf/Se56fA/u7/H4w/ewfzfv87CFaUEFPIiMSeccO/s7/O9/7ebWQLwe+A/enoh59wdzrlC51xhXl7vXNlRDlzauKPJ+OYbpA+fzJ0pv+Pilsf59sNLWfCHN3h2Zan25EViSDjhXgKMarc8Eihtt5wNTAVeM7OtwBxgoX5UjVE5I7DLF2FHfJHLmh/g3TF3kB2s4ZqHlnPqH//NolXbFfIiMSCccF8CTDSzcWaWAlwALNyz0TlX7ZzLdc6Ndc6NBd4FznTO6TjHWJWcDufcAaf9hqEV7/JkwnXcN9/RGgzy9QeXccatb/Lmxgq/qxSRbvQY7s65NuAa4AVgLfCYc261md1kZmf2doHiEzNvZqcr/oUlJnP8W1/lpTkr+d0Xp1Hd2MrFdy3mkrsWs7q02u9KRaQTOolJetZYBU9/E9Y9C4eeRvPpt3L/+9Xc+mox1Y2tfOHIAr538iRGDsrwu1KRuKczVCWynIPFf4F//QSyhsLZf6Z6+Gf4v9c2cc9bW3DApZ8ZyzdPmMCADM3bKtJbFO7SO7Ythb9fBbuKYc43Yd4NlNY7fvfiBp5cVkJ2ahLXzpvIV44ZS0pSWCdAi8h+ULhL72lpgBdvgCV3Qt5kOOd2GD6dtdtruHnROt7YUM743Ex+fPpk5h6ar7NdRSIooteWEfmUlAz4/G+8ST8aK+HOefDv3zJ5aCb3XT6bey49Cgwuv7eIr96zhOKyWr8rFul3tOcuB6dhtzdH65qnoKAQzrwFhk6hpS3Ife9s5Y8vb6ShJcAlc8bw3ZMmaTxe5CBpWEb6jnPwwZOw6AfQVA3HfgeO/z4kp7GrrpnfvriBR977iAHpyXxv/iQunD2apET90ShyIBTu0vfqd8G//gtWPOzN1XrGn/bO9LSmtIabnl3Nu5t3c+jQbH5y+hQ+OzHX54JFYo/G3KXvZQ6BL/wFLvmHd234e0+DhddCYyVTRuTw8JVz+MvFM2lobePiuxZz5X1FbK2o97tqkbikPXfpHS318NrN8M5tkD4I5t0AMy6BhESaWgPc9eYWbnu1mLaA47LPjuWauRPITtN4vEhPNCwj0WHHKnjuB/DR2zD8SDj1f2H00QDsrGnif59fz5PLSsjNSuUHpxzKebNGkpCgQydFuqJwl+ix5wfXf/0Eakth2gVw0o2QMxyA9z+u4qZnVrPsoyqOKBjAT8+YQuHYwb6WLBKtFO4SfZrr4M3fwdu3QGIKHPttmPMNSM3COcfCFaXc/Nw6dtQ0ccb0Efzw1MMoGJjud9UiUUXhLtFr92ZvL37ds5CZDydcBzO/ConJNLS08ZfXN3P765swg/93/CFc/blDSE9J9LtqkaigcJfo9/F78OJPvfH4wYfAvJ/AlLPBjJLKBm5etI5/rtzOiAFpXHfqYZw5fYQuZSD9nsJdYoNzsOEFeOlGKF8LI2bC3OthwklgxntbdvOzZ1azurSGwjGDuOGMKUwbOdDvqkV8o3CX2BIMwIpHvMMnqz+Gglnwuetg4skEHDyx9GN+/cJ6dtW3cN7MkXx/waHkZ6f5XbVIn1O4S2xqa4EVD8G/fwtVH8GIGV7IT1pAbXMbt75SzN1vbSElMYErjx/P144bT1Zqkt9Vi/QZhbvEtkCrtyf/799A5VYYNs07umbK2WypbObXL6zjuVU7yM1K4dvzJnLB7NEk63o10g8o3CU+BFph5WPw5u9h10YYONqbJGTGxSzf2covF61j8ZbdjB2SwX+eciifP2K4fnSVuBbRa8uY2QIzW29mxWb2w062f8/M1pjZSjN72czGHEjRIvtITIYZX4ZvvgcXPAw5BfD8dfD7w5mx8RYeuWgc91x6FKlJiVzz0HLOvu0t3t5U4XfVIr7rcc/dzBKBDcB8oARYAlzonFvTrs1cYLFzrsHMvg6c4Jz7Unevqz13OWAfL4G3/whrn4WEJJhyFoGjruQf5QX87sUNlFY3cfykPL570kRmjB7kd7UiERWxYRkzOwa40Tl3Smj5RwDOuZu7aD8DuNU5d2x3r6twl4O2axO8dye8/yA018CwI2ideQX31R/Frf8upbKhlRMPy+e7J03iiJED/K5WJCIiOSxTAHzcbrkktK4rVwCLwnhdkYMz5BA49ZfwvbVw+u8hGCT5ue9wxeLTWDzzJW4+NoGlH1Zyxq1vcuV9Rawurfa7YpE+E84xZJ39OtXp7r6ZXQwUAp/rYvtVwFUAo0ePDrNEkR6kZkHh5TDrMvjwbVhyJynL7ubC4O2cP3wGr2cu4L+KW/j8n3Zy8pShfGPuBI4cpROhJL5FbFjGzE4CbgE+55wr6+mNNSwjvaq+wjvKZvn9ULYGl5TO2kEn8vvyWbzcNIk5h+TxjRMmcOyEITq6RmJKJMfck/B+UJ0HbMP7QfUi59zqdm1mAE8AC5xzG8MpUOEufcI5KF0Gy+73LjvcXEN9Si4L2+bwSONsGDGTr8+dyMlThuo68hITInqcu5mdBvwBSATuds79wsxuAoqccwvN7CXgCGB76CkfOefO7O41Fe7S51obYcPzsOoJ3MZ/YYEWSmwYf2+dQ1H2POYddzznzRpJps54lSimk5hEutNYBeueJbjycWzLGxhB1gTH8HLCHFIPP4NTTzyRUUMy/a5SZB8Kd5Fw1e6E1f+gbtmjZJUtA2CrG8rGgccz9OjzmDr7JBKStDcv0UHhLnIgandQtfxpdi37B6OqlpBCG5UMYMfwuYyYcx4DJp8IKdqjF/8o3EUOUlNdJStfe5KWD55hWuNicqyRVkumfuhscqYuIGHiSZA/GXS0jfQhhbtIBG0s3cXiV58hsPEljg6+z2EJ3nl9bZnDSJp0EhwyD8afABma2Ft6l8JdpBc0tQZ4YfUOXlvyPqlbX+WzCSs5IekDslw9DsOGHQFjj4Oxn4Uxx0C6rm0jkaVwF+llZbVNLHy/lKeXfUTyjuUcl/gB8zOLOax1LUnBZsCgfdiPnqM9ezloCneRPrRhZy1Pv7+NRat2UFJRxZEJmzh38BaOT1nPsOqVWKDJa5g7CUbOhlFHefd5h0GCJhmR8CncRXzgnGPDzjqeW7WdRR9sZ8POOlJo5Qv52zl94EdMYwM5Fcuxhl3eE1JzvPliR82GgkIYPh2yh/rbCYkPkGsSAAALCElEQVRqCneRKFBcVsuLa8p4Zd1Oln5YSdBBbmYy541r4ZQBHzG5bT1pO4qgbA24oPekrGEw4kgv6IeH7nNG6KgcARTuIlGnsr6FNzaW8/LaMl5bX0ZNUxsAU4bncOL4dOYN3MnhtoWUslWw/X2o2PBJ4GfmhcJ+Ogw9HPKnwJAJ3kxV0q8o3EWiWFsgyMpt1bxdXMGbxRUs+7CKlkCQ5ERj2siBFI4ZxOyCVArTSxlQuQa2r/Bu5Wsh6H0pkJAMuRO9cfv8Kd4x9/mTYdBYSEj0tX/SexTuIjGksSXAkq27eWtTBUu27GbVtmpaA97/m+NzMykcO4jCMYOZNSqDca6UhPJ13lBO2VrvvurDT14sKR3yJsGQid7e/ZAJ3sQmQyZAWo5PPZRIUbiLxLCm1gCrtlWzZOtulm6tZOlHlVQ1tAKQnZrElBE5HFEwgKmh27gcR2LFBm/Pvix021UMVR/xqbl1MvM/CfvciTB4PAwcAwNHQ7omMIkFCneROBIMOjZX1LH0w0pWbatm1bYa1m6voaXNG5PPSEnk8BE5HD5iAJOGZjNpaBYTh2YzICkAlVu9oN972+Td13eYUydtwCdBP2hsu8ehe11TJyoo3EXiXGsgyKbyOlaVVLO6tIZV26pZU1pDY2tgb5uhOalMGprNxHwv8CfkZzE2N5MhmSlYcw3s3uzt3Vd+6N1XffjJ47bGT79hRi4MGAk5BZAz3DuCJ6cAsod/sk5fAL1O4S7SDwWDjm1VjWzYWcuGnXVs3FnLhrJaisvqaGoN7m2XnZrEmNwMxg7J9G65mYwdksHoIRnkZqaSYEB9eSjoP/wk9GtKQ7dt0FS1bwFpA9oF/nDIGuoNBWXlhe7zvSN/0gfp0M4DpHAXkb0CQUdJZQObyuvYWtHA1l31bN3VwNaKekoqGwi2i4GUxASGD0yjYGA6IwamUxC6jRiYTsGgdIYPSCMtORFaGqB2uxf0NXvuS9utK/XmsnWBfQtKSPZCvmPoZ+V7yxmDvC+A9MHeJRtSc/RlEBJuuGsGApF+IDHBGDMkkzGdzC7V0hZkW1UjWyvq+biygW1VjWyrbKS0qpE3N1aws7aJjvuAgzKSyctOJT87jfzsHPKy88jLPor8kWnkZaWSn5NKXnYq2SkJWGMl1JV5Y/x15aH7Mu8vgz3rd672loOtnXcgIenTYZ8+OPQF0G45faD3l0Nqjne/53FSSi/8i0Y/hbtIP5eSlMC43EzG5XY+Xt7SFmRnTRMlocDfVtVIWW0TZTXNlNU2s6WinvLaZloCwX2em5xoDMpIYXBmCgMzUhmcOZ5BGYd5y7kpDB6T/Mn2tGRyqCMrUElScxU07IbG3R3uK71b5VZv4vOG3RBo7r6DSWntAj+nk8cDvccpmZCSFbpltruFlpMzYuo6QAp3EelWSlICowZnMGpwRpdtnHNUN7ZSXusFflltE+W1zeyub6WqoYXd9S1UNrSwfkctVQ2tVDa0fGooqKP05ESy01LISR9Fdto4stOSyU5LIictiezsZLJTk8hOSyI7NYmcpFZyXC3Z1JERrCc9WE96oI6UQB0pbXUkNFdDUw0010BT6HF1SWi5Zt8fjrtkHUK/XfDv/VLIgOR074sgOd075yA5/dPrkjO8H6YHFOzfB7Gfwgp3M1sA/BFIBP7qnPtlh+2pwH3ALGAX8CXn3NbIlioi0crMGJiRwsCMFCYOze6xfTDoqGlqpbKh1Qv+UPjXNrWFbq3efbN3X93YSkllw95t7X8c3lcKMDh0g9SkBDJSEslISSI9JZGMlETvPstbl5EUYIA1kmUtZCU0kkkzGTSRQRPpNJLumkhzjaQGG0kJNpESaCA52EhyoIGktgYSa3eR0Poxia11WGsDtDVibU3d/wMc+x2Y/7Ow/30PRI/hbmaJwG3AfKAEWGJmC51za9o1uwKodM5NMLMLgF8BX+qNgkUk9iUkfPJl0NVwUHda2oLUNXtBX9fcRmNLgIaWAA0tbaH7QOfrWkOPmwOU1zXT0NJAc2uQlkCQ5tYAzW1JNLclAOkH1T8jSFZCK9mJbWQntpCd2EpWQiuZCa1kJrRwVGAaXzyod+hZOHvus4Fi59xmADN7BDgLaB/uZwE3hh4/AdxqZub8OhRHROJaSlICg5O8sfpIc855Yd8WpKXNu/eC/5PH3pdBaLkt8Kn1rQFHayAYujla2oKfWq4PBEke2vuXdQ4n3AuAj9stlwBHd9XGOddmZtXAEKCifSMzuwq4CmD06NEHWLKISO8xM1KTEklNiu2Lr4Xz029nB5d23CMPpw3OuTucc4XOucK8vLxw6hMRkQMQTriXAKPaLY8ESrtqY2ZJwABgdyQKFBGR/RdOuC8BJprZODNLAS4AFnZosxD4aujxecArGm8XEfFPj2PuoTH0a4AX8A6FvNs5t9rMbgKKnHMLgbuA+82sGG+P/YLeLFpERLoX1nHuzrnngOc6rLuh3eMm6PUje0REJEyxcy6tiIiETeEuIhKHFO4iInHIt+u5m1k58GGPDTuXS4cTpGKQ+hAd1Af/xXr90Ld9GOOc6/FEId/C/WCYWVE4F6uPZupDdFAf/Bfr9UN09kHDMiIicUjhLiISh2I13O/wu4AIUB+ig/rgv1ivH6KwDzE55i4iIt2L1T13ERHphsJdRCQORXW4m9kCM1tvZsVm9sNOtqea2aOh7YvNbGzfV9m9MPpwqZmVm9n7odvX/KizK2Z2t5mVmdkHXWw3M/tTqH8rzWxmX9fYkzD6cIKZVbf7DG7orJ1fzGyUmb1qZmvNbLWZfbuTNlH9OYTZh2j/HNLM7D0zWxHqwz6ToEZVJjnnovKGdwXKTcB4vBlvVwBTOrT5BvCX0OMLgEf9rvsA+nApcKvftXbTh+OBmcAHXWw/DViEN2HLHGCx3zUfQB9OAJ71u85u6h8OzAw9zgY2dPLfUVR/DmH2Ido/BwOyQo+TgcXAnA5toiaTonnPfe/crc65FmDP3K3tnQX8LfT4CWCemXU2K5RfwulDVHPOvUH3E6+cBdznPO8CA81seN9UF54w+hDVnHPbnXPLQo9rgbV4U1u2F9WfQ5h9iGqhf9u60GJy6NbxiJSoyaRoDvfO5m7t+B/Dp+ZuBfbM3RotwukDwLmhP6WfMLNRnWyPZuH2MdodE/pze5GZHe53MV0J/Zk/A2+vsb2Y+Ry66QNE+edgZolm9j5QBrzonOvyc/A7k6I53CM2d6uPwqnvGWCsc24a8BKffOvHimj/DMKxDO96HdOBW4CnfK6nU2aWBTwJfMc5V9NxcydPibrPoYc+RP3n4JwLOOeOxJtudLaZTe3QJGo+h2gO93iYu7XHPjjndjnnmkOLdwKz+qi2SAnnc4pqzrmaPX9uO29immQzy/W5rE8xs2S8UHzQOff3TppE/efQUx9i4XPYwzlXBbwGLOiwKWoyKZrDPR7mbu2xDx3GRc/EG4uMJQuBr4SO1pgDVDvntvtd1P4ws2F7xkXNbDbe/xe7/K3qE6Ha7gLWOud+10WzqP4cwulDDHwOeWY2MPQ4HTgJWNehWdRkUljT7PnBxcHcrWH24VozOxNow+vDpb4V3AkzexjvKIZcMysBfor3QxLOub/gTb94GlAMNACX+VNp18Low3nA182sDWgELoiynYRjgUuAVaHxXoDrgdEQM59DOH2I9s9hOPA3M0vE++J5zDn3bLRmki4/ICISh6J5WEZERA6Qwl1EJA4p3EVE4pDCXUQkDincRUTikMJdRCQOKdxFROLQ/weyjWmo4E8bkQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x210e49ee550>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xs = [float(i/50) for i in range(int(50*PI))]\n",
    "fs = [f_of_x(x) for x in xs]\n",
    "gs = [g_of_x(x, A=1.4, lamda=1.4) for x in xs]\n",
    "plt.plot(xs, fs)\n",
    "plt.plot(xs, gs)\n",
    "plt.title(\"f(x) and g(x)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determine the Optimal Parameters for Our Weight Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def inverse_G_of_r(r, lamda):\n",
    "    return (-1 * math.log(float(r)))/lamda"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_IS_variance(lamda, num_samples):\n",
    "    \"\"\"\n",
    "    This function calculates the variance if a Monte Carlo\n",
    "    using importance sampling.\n",
    "    Args:\n",
    "    - lamda (float) : lamdba value of g(x) being tested\n",
    "    Return: \n",
    "    - Variance\n",
    "    \"\"\"\n",
    "    A = lamda\n",
    "    int_max = 5\n",
    "    \n",
    "    # get sum of squares\n",
    "    running_total = 0\n",
    "    for i in range(num_samples):\n",
    "        x = get_rand_number(0, int_max)\n",
    "        running_total += (f_of_x(x)/g_of_x(x, A, lamda))**2\n",
    "    \n",
    "    sum_of_sqs = running_total / num_samples\n",
    "    \n",
    "    # get squared average\n",
    "    running_total = 0\n",
    "    for i in range(num_samples):\n",
    "        x = get_rand_number(0, int_max)\n",
    "        running_total += f_of_x(x)/g_of_x(x, A, lamda)\n",
    "    sq_ave = (running_total/num_samples)**2\n",
    "    \n",
    "    \n",
    "    return sum_of_sqs - sq_ave"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal Lambda: 1.6500000000000001\n",
      "Optimal Variance: 0.0465612105426072\n",
      "0.002157804683992673\n"
     ]
    }
   ],
   "source": [
    "# get variance as a function of lambda by testing many\n",
    "# different lambdas\n",
    "\n",
    "test_lamdas = [i*0.05 for i in range(1, 61)]\n",
    "variances = []\n",
    "\n",
    "for i, lamda in enumerate(test_lamdas):\n",
    "    print(f\"lambda {i+1}/{len(test_lamdas)}: {lamda}\")\n",
    "    A = lamda\n",
    "    variances.append(get_IS_variance(lamda, 10000))\n",
    "    clear_output(wait=True)\n",
    "    \n",
    "optimal_lamda = test_lamdas[np.argmin(np.asarray(variances))]\n",
    "IS_variance = variances[np.argmin(np.asarray(variances))]\n",
    "\n",
    "print(f\"Optimal Lambda: {optimal_lamda}\")\n",
    "print(f\"Optimal Variance: {IS_variance}\")\n",
    "print((IS_variance/10000)**0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmcHHWd//HXp3vuK5PJTM7JSRJCAsRADMiNogRcCF4sKJdyqLt4rbqLuj9/yq676m9XcRc8AA8OBcGDRRcIIiCQmJAAJiSZECYh5JpJZiaZzH329/dHVYdOO5PpSXq6unvez8djHtNdVd396eqqd1d/q+pb5pxDRESySyjoAkREJPkU7iIiWUjhLiKShRTuIiJZSOEuIpKFFO4iIllI4T7CzGyambWZWTjoWhJhZoVm9jszO2hmDwddT7KZ2Q/N7P/E3P+kme31P6NxZnammb3u378syFqDZGbbzeyCYUzvzGz2SNY0yOt+zczuT/XrZgKFewwzW25mtw4wfJmZ1ZtZznCf0zm3wzlX4pzrT06VI+6DwARgnHPuQ/Ej/ZXJmdmn44Z/1h/+tZhhZWZ2m5nt8MOy1r9fOdyiElmJ/UDqNLNWM2s2s5Vm9gkzO7ScO+c+4Zz7F3/6XOA7wHv8z6gJuBW43b//yHDrPBZmdp2ZvTDENM+a2Q2pqmkkmdkUM+szs+MGGPdbM/uPIOrKFgr3w/0MuNrMLG741cDPnXN9w3myo/kySAPTgS1DvNctwLVxw67xhwNgZnnAH4EFwFKgDDgDaAKWJLPgOJc450rx3sc3gX8CfjzItBOAAmBjzLDpcfcTlqGfd2Ccc7vxlpGrY4ebWQVwMXBPEHVlDeec/vw/oBA4CJwTM2ws0AUs9O+/F3gFaAF2Al+LmXYG4IDrgR3AczHDcvxpPgrUAK3ANuDjMY8/D9gFfB7YB9QBH42r7z+BN/06XwAK/XGnAyuBZmAdcN4R3ucJwLP+tBuBS/3hXwd6gF6gDbh+gMd+Dbjffw8L/GEL/Pv3R+cHcAOwFygZxvz/nj9PW4CXgLP94Uvj6lo3yOO3AxfEDVsCRIAT/fs/A/4VmAu0+59NG/A0sNWfttMflg+MwftyqAN2+48N+891HbAC+C6wH/hXf/jH/PlxAFgOTI+pxwGfAF73x98BmP+ZdAH9/ms3D/IenwVuGGTcw0C9v2w8F/18Yt7394HH/edfAUwEbvPr2AwsipuXXwI2+eN/ChTEjP+iP0/2+O/XAbOHWkcGqPnDwNa4YX8HvDzUchG7PMauP4MtE3gbs7f4n3MT8BBQ4Y8rwFt+m/DWizXAhKAz6Vj+Ai8g3f6Au4C7Y+5/HPhLzP3zgJP8BeVkvAC7zB83w1/I7wWK8cI4Oiwa7u8FjvNX6HOBDuCUmOfuw2sayMXbeukAxvrj7/BX7ilAGG9LON+/3+RPHwLe7d+vGuD95QK1wJeBPOCdeF80x/vjD60sg8yfr/krwZeBb/nDvo0XBLHh/iBwzzDn/VXAOCAH7wuuHj9QhqrLn+bQihw3fAfwSf/2z3grhA/7bAZ6DuAR4Ef+5zkeeBH/Cxkv3PuAT/k1FwKX+fP3BH/YPwMrY57PAb8HyoFpQAOwNOb5XhjiPT7L4OH+MaDUXyZu4/Dl9mdAI3AqXpA9DbyB94srjPel9UzcfNgATAUq8L4MovNtKd5yf6I/X37B4eF+HoOsIwPUHN2gOitm2J+Bzw53uWDocP8ssAqo9ufRj4AHYtbz3wFF/vw4FSgLOo+O5S/wAtLtDzjLX9iiW8QrgM8dYfrbgO/6t2f4C/msmPHRYTmDPP4R4DP+7fPwthpjw2Yf3lZ5yB+3cIDn+Cfgvrhhy4FrB5j2bH/lCMUMe4C3QvnQyjJIvV/DC/FpeKGZ6/+fyuHh/gfgm8f4WRzgrV9MR6zLn+bQihw3fBXwFf/2z0gw3PGabbqjy4I/7Er8EMQL4x1xr/U4Mb94/M+tA3/r3X+92CB7CLgl5vmOOtzjpiv3X2tMzPu+K2b8p4CamPsnEfNrwZ8Pn4i5fzH+Fjbwk9jPFu9X0KFwP9I6Msj4u4E7/dtz8H6ljR/ucsHQ4V4DvCtm3CS8X4M5eF+MK4GTj2WZTac/tbnHcc69gLc1tczMZgFvx9syAcDMTjOzZ8yswcwO4v3Ejt9BuHOw5zezi8xslZntN7NmvJUm9vFN7vD27g6gxJ+mAO8nZbzpwIf8nYjN/vOehbfwxpsM7HTORWKGvYm39Z8w59wOvC3UfwNed87Fv+emQV5/UGb2eTOr8Y/UacZrEhn2ztcBTMFrNhmu6XhfXnUx8/VHeFvwUfHvezrwvZjp9+P9Soudv/Uxt6Of7zExs7CZfdPMtppZC16oweHzb2/M7c4B7sfXEfve3sRbdvD/x4+LrSWRdSTWPcDlZlaA1/7+hHNuX8zzJWu5mA78NuazqcFrBpsA3Ie3QfSgme0xs2/7O9wzlsJ9YPfi/Vy9GnjSORe7EvwCeBSY6pwbA/wQb+WN5QZ6UjPLB34N/Adee1458NgAjx9II16b7F8dWYC3ot3nnCuP+St2zn1zgGn3AFNjjyDB2wrfnUAN8e7F+5l87wDjngIuNLPiRJ7IzM7G+wVyOV4zVDneL6jovBlwnibwvG/HC9YjHoUyiJ14W+6VMfO1zDm3IGaa+Lp24jXbxH4Whc65lQm83lG9R9+HgWXABXjhN8MfnsiyNZipMben4S074LW1x4+Llcg6cohz7nm8jYFleE0wh5anBJaLWO14zSrRx4aBqpjxO4GL4j6bAufcbudcr3Pu6865+XjNnX+DlwEZS+E+sHvxVpIb+es99qXAfudcl5ktwVupEpWH19bXAPSZ2UXAexJ5oL+l/RPgO2Y22d9Se4f/hXE/cImZXegPLzCz88yseoCnWo23EvyjmeWa2XnAJXht5MP1S7/+hwYYdx/eyvRrM5tnZiH/OPIvm9nFA0xfitd+3QDkmNlX8Y6widoLzIj7UhqUfxjm3+C9r/udc68m/rY8zrk64EngP/3nC5nZcWZ27hEe9kPgS2a2wK9jjJn91SGlg9gLVPtHGh1Jjv8ZR/9y8eZfN15IFuH9ojpWf29m1f7RK1/G+7zB+7yvM7P5ZlYE/N+4xx3NOnIv8C285qTfxT3XkZaLWFuAAjN7rz9P/hlvfYv6IfANM5sOYGZVZrbMv32+mZ3kfyG04DXXZMrhywNSuA/AObcdr/2tGG8LJNbfAbeaWSvwVQYOtsGetxX4tP+YA3gLffzzH8kXgFfx9uTvx1sZQn6TyDK8FbABL1S/yACfr3OuB7gUuAjv18D3gWucc5uHUUf0uTqdc0855zoHGNeN9wW5Ga/9vQVvZ2Ql3hdMvOV47dVb8H7md3H4T//oCVVNZvbyEcr6nf/Z7AS+gncc+0eH877iXIP3pRw9auRXHKG5yTn3W7zP5UG/eWQD3rxOxNN4Ry/Vm1njEab7AV4zSvTvp3jh+CbeL7BNePsZjtUv8L7ctvl//wrgnHscrx39abymuafjHnc068i9eL8AfukvO1FDLReHOOcO+q99N958aMc7+izqe3jr25N+bauA0/xxE/E+2xa85po/4W00ZSzzdyyIiEgW0Za7iEgWUriLiGQhhbuISBZSuIuIZKHAOjqqrKx0M2bMCOrlRUQy0ksvvdTonKsaarrAwn3GjBmsXbs2qJcXEclIZvbm0FOpWUZEJCsp3EVEspDCXUQkCyncRUSykMJdRCQLKdxFRLKQwl1EJAsp3EVEUqS7r59vPbGZdTubR/y1FO4iIinS0NrND57dymv1rSP+Wgp3EZEU2dfqXYekqjR/iCmPncJdRCRFGhTuIiLZR+EuIpKFGlq7MYOK4qGugX7sFO4iIinS0NZNRVEeueGRj16Fu4hIijS0dqekSQYU7iIiKaNwFxHJQgp3EZEs45yjoU3hLiKSVVq6+ujpi1BVonAXEckaqTzGHRTuIiIpsa+1C1C4i4hkleiW+3iFu4hI9jjULFNSkJLXU7iLiKRAQ1s3eeEQZYU5KXk9hbuISApEj3E3s5S8nsJdRCQFGlq7qUxRezso3EVEUqKhtTtlx7iDwl1EJCUa27oZX6ZwFxHJGn39EZrae7TlLiKSTfa39+Bc6k5gAoW7iMiIS+WFsaMU7iIiIyzV/cqAwl1EZMS9dXaqwl1EJGs0tGnLXUQk6zS0dlNakENBbjhlr6lwFxEZYam8vF6Uwl1EZISl+uxUULiLiIy4VF47NSqhcDezpWb2mpnVmtktA4yfZmbPmNkrZrbezC5OfqkiIpmpobWb8aWp6cc9ashwN7MwcAdwETAfuNLM5sdN9s/AQ865RcAVwPeTXaiISCbq6OmjrbsvLbfclwC1zrltzrke4EFgWdw0Dijzb48B9iSvRBGRzNXY2gOk9jBISCzcpwA7Y+7v8ofF+hpwlZntAh4DPjXQE5nZTWa21szWNjQ0HEW5IiKZJdUXxo5KJNwHumyIi7t/JfAz51w1cDFwn5n91XM75+50zi12zi2uqqoafrUiIhkmiLNTIbFw3wVMjblfzV83u1wPPATgnPszUABUJqNAEZFMFsTZqZBYuK8B5pjZTDPLw9th+mjcNDuAdwGY2Ql44a52FxEZ9RpauwkZVBTnpfR1hwx351wfcDOwHKjBOypmo5ndamaX+pN9HrjRzNYBDwDXOefim25EREadhtZuxpXkEw6l5sLYUTmJTOScewxvR2nssK/G3N4EnJnc0kREMl8QZ6eCzlAVERlRQZydCgp3EZERFUSnYaBwFxEZMZGIo1Fb7iIi2eVgZy+9/Y7xCncRkewR1DHuoHAXERkx+1qCOTsVFO4iIiOmoS2YfmVA4S4iMmIO9SujcBcRyR4Nrd0U5IYoyU/ofNGkUriLiIyQ6DHuZqntegAU7iIiI6ahLZiuB0DhLiIyYoI6OxUU7iIiI0bhLiKSZXr6Ihzo6KWqpCCQ11e4i4iMgKb24A6DBIW7iMiIiJ6dGkS/MqBwFxEZEUGewAQKdxGRERFkp2GgcBcRGRHRLfdxJam9MHaUwl1EZAQ0tHZTXpRLfk44kNdXuIuIjICgLowdpXAXERkBQV0YO0rhLiIyAoI8OxUU7iIiSeecU7OMiEi2ae/pp7O3X1vuIiLZJOgTmEDhLiKSdPtagrt2apTCXUQkyaJnp44vDaZHSFC4i4gknZplRESyUENrNzkho7wwN7AaFO4iIknW0NpNZUk+oVDqL4wdpXAXEUmyoM9OBYW7iEjSBX12KijcRUSSLuizUyHBcDezpWb2mpnVmtktg0xzuZltMrONZvaL5JYpIpIZ+iOOpvaewLfcc4aawMzCwB3Au4FdwBoze9Q5tylmmjnAl4AznXMHzGz8SBUsIpLODnT00B9xgYd7IlvuS4Ba59w251wP8CCwLG6aG4E7nHMHAJxz+5JbpohIZkiHY9whsXCfAuyMub/LHxZrLjDXzFaY2SozW5qsAkVEMsm+NAn3IZtlgIEO1HQDPM8c4DygGnjezE50zjUf9kRmNwE3AUybNm3YxYqIpLvolvv4DNhy3wVMjblfDewZYJr/cc71OufeAF7DC/vDOOfudM4tds4trqqqOtqaRUTSVjTcKzPgaJk1wBwzm2lmecAVwKNx0zwCnA9gZpV4zTTbklmoiEgmaGjtpjgvTHF+Ig0jI2fIcHfO9QE3A8uBGuAh59xGM7vVzC71J1sONJnZJuAZ4IvOuaaRKlpEJF2lw9mpkFibO865x4DH4oZ9Nea2A/7B/xMRGbUaWrvSItx1hqqISBKlQ9cDoHAXEUmqdOh6ABTuIiJJ09XbT0tXn7bcRUSySWNbepzABAp3EZGkSZezU0HhLiKSNIf6lSkJ7sLYUQp3EZEkSZdOw0DhLiKSNA2t3ZjBuJK8oEtRuIuIJEtDWzcVRXnkhoOP1uArEBHJEulyAhMo3EVEkkbhLiKShdLl7FTIwHBfUdvIFx5eh9dXmYhIenDOpU2PkJCB4b67uZNfvbSL1W/sD7oUEZFDGtt66OmLMKEs+GPcIQPD/ZKTJzOmMJf7V70ZdCkiIofU1LUAMG9SacCVeDIu3Avzwnzw1Gqe2FDPvtauoMsREQHeCvf5k8oCrsSTceEO8JHTptEXcTy0ZmfQpYiIALCproXJYwooLwr+BCbI0HCfVVXCWbMr+cXqHfRHtGNVRIK3aU8L8yenx1Y7ZGi4A1x1+nT2HOzi6c37gi5FREa5rt5+tjW2c0KaNMlABof7BSeMZ2JZAfdpx6qIBGzL3lb6Iy5t2tshg8M9JxziyiXTeG5LA282tQddjoiMYtGdqdpyT5IrlkwlHDJ+vnpH0KWIyCi2aU8LxXlhplUUBV3KIRkd7hPKCnjP/Ak8tHYnXb39QZcjIqNUTV0r8yaVEQpZ0KUcktHhDnD16dNp7ujlf9fXBV2KiIxCkYhjU11LWrW3QxaE+zuOG8esqmLuX60dqyKSersOdNLW3ZdWh0FCFoS7mXHVadN5ZUczG3YfDLocERllNqXhzlTIgnAH+MCp1RTkhvi5tt5FJMU21bUQMjh+Qnr0KROVFeE+pjCXZQun8Mgre2jp6g26HBEZRWrqWphZWUxhXjjoUg6TFeEO3hmrnb39/OalXUGXIiKjiNftwJigy/grWRPuJ1WPYeHUcu5fvUMX8hCRlDjY0cvu5k5OSJNufmNlTbgDXHXaNGr3tbFqmy7kISIjr6Y+vbr5jZVV4X7JQl3IQ0RSZ9MeP9zT7DBIyLJwL8gN86FTq1m+sZ59LbqQh4iMrJq6FipL8hhfmh6X1ouVVeEO8JHTp9MXcTyoC3mIyAjbVNeSdse3R2VduM+sLObsOZU88OIO+vojQZcjIlmqtz/C63vb0rK9HRIMdzNbamavmVmtmd1yhOk+aGbOzBYnr8Th+8hp06g72MWftjQEWYaIZLGtDW309EfSsr0dEgh3MwsDdwAXAfOBK81s/gDTlQKfBlYnu8jhetcJE6gsyeeBF9UVsIiMjOjO1ExullkC1DrntjnneoAHgWUDTPcvwLeBwPdk5oZDXL64mqc376P+YODliEgWqqlrIS8nxKzK4qBLGVAi4T4FiN07ucsfdoiZLQKmOud+f6QnMrObzGytma1taBjZJpO/fftUIg4eWqsdqyKSfJvqWpg3sZSccHruukykqoF6nz90CqiZhYDvAp8f6omcc3c65xY75xZXVVUlXuVRmD6umLNmV/LLNTvpj+iMVRFJHuccNXWtnDAxPZtkILFw3wVMjblfDeyJuV8KnAg8a2bbgdOBR4PeqQpw5ZJp7G7u5PnXtWNVRJJnb0s3+9t70nZnKiQW7muAOWY208zygCuAR6MjnXMHnXOVzrkZzrkZwCrgUufc2hGpeBjePX8C44rztGNVRJJqU5137Yh03ZkKCYS7c64PuBlYDtQADznnNprZrWZ26UgXeCzyckJ88NRqnqrZpzNWRSRpaupaAZiXhh2GRSW0J8A595hzbq5z7jjn3Df8YV91zj06wLTnpcNWe9Tfvn0q/RHHw+oKWESSZNOeFqZWFFJWkBt0KYNKz928STSrqoTTZ1Xw4JodRLRjVUSSoCYNL4gdL+vDHbwdqzv3d7Jia2PQpYhIhmvv7uONpnbmT0q/C3TEGhXhfuGCiYwtyuXBF3XMu4gcm831rThHWl6gI9aoCPeC3DDvP8XrCrihtTvockQkg9XUpW8f7rFGRbgDXLlkKn0Rx69f1o5VETl6m+paKCvIYUp5YdClHNGoCffZ40tZMqOCB1/UNVZF5OjV+H24mw108n76GDXhDnDFkqlsb+rgz9uagi5FRDJQf8Sxua41rU9eihpV4X7xSZMoK8jhAe1YFZGj8GZTO529/Wnf3g6jLNwP7VjdUM/+9p6gyxGRDLMpujNVW+7p58ol0+jpj/Ab7VgVkWHatKeFnJAxZ0JJ0KUMadSF+/ETSzllWjm/0I5VERmmmroWZo8vIT8nHHQpQxp14Q7e1vu2hnbWbD8QdCkikkE2+UfKZIJRGe7vPXkSpfk56gpYRBLW1NbN3pbujGhvh1Ea7kV5Obz/lCn8bt0eVumwSBFJQLSbX225p7nPX3g808YV8fc/f5ndzZ1BlyMiae6tC3Skd58yUaM23MsKcrnrmsX09EX4+H1r6ertD7okEUljNXWtTCwrYFxJftClJGTUhjvAcVUlfO/Kt7FxTwtf+s2rOnpGRAa1aU9Lxmy1wygPd4B3zpvA5989l9++spsfv/BG0OWISBpq6eqltqGNBZPTuw/3WKM+3AH+/vzZXHTiRP7tsRpeeF0X9BCRwz2zeR/9Ecf586qCLiVhCnfAzPiPDy1kzvhSbn7gZXY0dQRdkoikkeUb66kqzWfR1LFBl5IwhbuvOD+HO685FefgpvvW0tHTF3RJIpIGunr7eWZzA++ZP4FQKL27+Y2lcI8xfVwx/33lIrbsbeWLD6/XDlYR4bktDXT29rP0xIlBlzIsCvc458yt4paL5vG/r9bx/We3Bl2OiATsiY31lBXkcPqscUGXMiwK9wHcePYsLl04mf948jWe3rw36HJEJCC9/RGe2rSXC+ZPIDecWXGZWdWmiJnxrQ+czPxJZdxwz1o+9cArbK5vCbosEUmxVduaaOnqY+mCzGqSAYX7oArzwtx//WnceM4snq7Zy9LbnueGe9bwyg71JCkyWjyxoZ7C3DDnzM2cQyCjFO5HMLY4jy9ddAIrbnknn7tgLmvfPMD7vr+Sj9y9ipVbG7XDVSSL9Uccyzfu5fx5VRTkpn//7fEU7gkoL8rjMxfM4YV/eidfvngeW/a28eG7VvOBH6zk6c17FfIiWeiVHQdobOvmwgxskgGF+7CU5Odw0znH8fw/ns+/XHYie1u6+djP1vKRu1fT3aeOx0SyyfKN9eSFQ7xz3vigSzkqCvejUJAb5urTp/PsF8/j65cuYOXWJr78mw3aghfJEs45nthYz5mzx1FakBt0OUdF4X4McsMhrj1jBp+9YA6/fnkXdz+vjsdEssGmuhZ27u/MuBOXYinck+DT75zDRSdO5N8fr+GZ1/YFXY6IHKPlG+oJGVxwwoSgSzlqCvckCIWM/7x8IcdPLOPTv3iF2n1tQZckIsfgiY31LJlZkTEX5hiIwj1JivJyuOuaU8nLCXHjvWs52NEbdEkichS2NrSxZW9bxh4lE6VwT6LqsUX86OpT2XWgg5sfeJm+/kjQJYnIMC3fWA8wOsLdzJaa2WtmVmtmtwww/h/MbJOZrTezP5rZ9OSXmhkWz6jgG5edxPOvN/KNx2qCLkdEhmn5hnoWVo9hcnlh0KUckyHD3czCwB3ARcB84Eozmx832SvAYufcycCvgG8nu9BMcvnbp/LRM2fw0xXb+eWaHUGXIyIJ2tPcybpdB7kwg4+SiUpky30JUOuc2+ac6wEeBJbFTuCce8Y5F7180SqgOrllZp6vXHwCZ8+p5J8f2cCa7fuDLkdEEvCk3ySTiR2FxUsk3KcAO2Pu7/KHDeZ64PGBRpjZTWa21szWNjQ0JF5lBsoJh7j9ylOoHlvEJ+57iV0HdOk+kXT3xMZ65k4oYVZVSdClHLNEwn2g60oNeCqmmV0FLAb+30DjnXN3OucWO+cWV1VlXi9rwzWmKJe7rllMT3+E9/7XC3zh4XU8tWkvXb3qqkAk3TS1dfPiG/uzYqsdICeBaXYBU2PuVwN74icyswuArwDnOue6k1Ne5ps9voRf3HA6P13xBss31vOrl3ZRnBfm/HnjWXriRM47fjwl+Yl8DCIykp6q2UvEkRXt7ZBYuK8B5pjZTGA3cAXw4dgJzGwR8CNgqXNOp2jGOal6DN/527fR0xfhz9uaeGJDPX/YVM/v19eRlxPinDmVXLhgIu+eP4HyorygyxUZlZ7YUM/UikLmTyoLupSkGDLcnXN9ZnYzsBwIAz9xzm00s1uBtc65R/GaYUqAh80MYIdz7tIRrDsj5eWEOHduFefOreJfLzuRtdv388TGepZvqOepmn3khUOcP6+K9y2awvnzxpOfk3l9SItkopauXlbUNnHtGdPxMyzjWVA9GS5evNitXbs2kNdON8451u86yKPr9vA/f9lDY1s3ZQU5vPfkybz/lCksnj42axY4kXT0P3/ZzWce/Au//uQ7OHV6RdDlHJGZveScWzzUdGrsTQNmxsKp5SycWs6XLprHiq1N/PblXTzyym4eeHEHUysKuextU3jfoilZsRdfJN0s31hPVWk+i6aODbqUpFG4p5mc8FtNN+3dfSzfWM9vX9nNHc/U8t9P17JkRgU3njOLd80bTyikrXmRY/VGYzt/rNnH5YunZtU6pXBPY8X5Obz/lGref0o1e1u6eOSV3dz75ze58d61HFdVzI1nz+KyRVMy8vqOIumgP+L4wsPrKMgNc/M7ZwddTlKp47AMMaGsgI+fexx/+uJ5fO+Kt1GQG+aW37zKWd96htuffp3mjp6gSxTJOD954Q1eevMAX790ARPKCoIuJ6m0QzVDOedYubWJO5/bxp+2NFCUF+byxVO5/qyZTK0oCro8kbRXu6+Vi//rBc6bW8WPrj41Yw5a0A7VLGdmnDm7kjNnV7K5voU7n9vG/ave5L5Vb3Le3CouXDCRC+ZPoKJYx82LxOvrj/D5h9dTnBfmG+87KWOCfTgU7llg3sQyvnP52/jihcdzz8o3+d26Pfxx8z5Cv4G3z6jgwgUTufDEiUzJ8C5MRZLlR89tY93OZm7/8CKqSjP3aktHomaZLOScY+OeFp7cWM8TG+vZste77N9JU8Zw4YIJXLhgItPHFRNxjv6Io985IhH/dvS+g4llBYSz6OgBEYDN9S1c8t8v8J4FE7njw6cEXc6wJdoso3AfBd5obGf5xnqWb6znlR3NCT9u+rgiPnnucbz/lGrycrTvXTJfb3+Ey+5Ywd6WLp783LkZ2WypcJcB7W3p4o81+9jf3k0oZITNCIfe+gv59/v6Izz80i7W7zrIpDEFfPycWVyxZJoOu5SMdttTW7jtqdf54VWnsjRDOwhTuMsxc87x/OuN3P50LS9u309lST43nD2Tq06frp4sJeNs2H2Qy+5Ywd+cPInbrlioUh/uAAAM8ElEQVQUdDlHTeEuSbV6WxO3P1PL8683MqYwl4+eOYOPnjGTMUW59PRFqD/Yxa7mDvY0d7H7QCd7mjvZc7CTuoNdVBTlMauqmJmVxcyqKmFWVTHTKorIDaupR1Kju6+fZbevYH97D09+7pyM7n1Vh0JKUp02axynzRrnHWHwTC23PfU6dz23jZKCHPa1dhO/jVBVms+U8kKOqypmf3sPf9i0l6b2t060CoeMaRVFzKosZkZlMaUFOeSGQ+SGjZxQiNycEHkxt/NzQsweX8LMccVZdYq4pMZ//7GWzfWt/PjaxRkd7MOhcJdhWTi1nLuuWUxNXQv3rNxOf8QxZWwhk8sLqS73/k8qLxiwu+KDHb1sa2xjW0M72xrbeKOxnW0N7azY2khXbySh1y/Nz+HEKWM4uXoMJ1eXc3L1GKrHFmblccqSHOt2NvODP23lg6dW864TJgRdTsqoWUbSQiTi6I1E6O139PVH6Ol/63Zvf4SOnn4217eyflcz63cdpKauhd5+b9kdW5TLSdXlnDxlDHMmlDCrsoSZVcXaLyA0tnXzoR/+mc6efpZ/7hzGFOYGXdIxU7OMZJRQyMgPhTlSHp9cXc7li70rPnb39fNafSvrdx08FPg/qG2kP/LWxsqEsnxmVXpt/NG2/uMqS6geW6imnVFgf3sPH7lrNXUHO7nv+tOyItiHQ+EuGSk/J+w3y5QD0wHo6u3nzaYOtjW0sa2xna0NXtPP79fXcbCz99Bji/LCzJ1QyryJpRw/sZR5E8uYN7GUsRl4zLMMrLmjh6vuXs32pnZ+et3befuM9L4Ax0hQuEvWKMgNc7wf2LGcc+xv7/ECf18bm+tbea2+leUb63lwzc5D040vzWfepDLmji+hsjSf8sJcyotyKS/Ko7wol7FFeYwpzP2rY/0jEUd3X4Su3n66+vrp7Ok/tA9hZmUxhXk6NyCVDnb2cvWPX6R2Xxt3X7uYM2ZXBl1SIBTukvXMjHEl+YwryT9sC845R0NrN5vrW9lc3+L9r2tl9bYmuvsG38FbmBumpCCHHj/QjzRtyGDGuGLmTfJ+IRw/sZQTJpYN2TTUH3G0dvXS0tlHXyTC5PJCnUCWgNauXq75yYteZ3pXL+acuVVBlxQYhbuMWmbG+LICxpcVHBYCzjk6evpp7uyluaOH5o5emjt6OdDRw8HOXg6099DW3Ud+ToiC3DD5uWEKc8MU5Hr3C3JDFOaG6Ys4Xt/bxub6FjbuaeHxDfWHDhktzgszd2Ip0yqKaO/up6Wrl5bOXlq7+rz/3X1/VW9lST5TKwqpHlvE1LH+f//+lPLCUd9FRFt3H9f9dA0bdx/kB1edyvnzxgddUqAU7iJxzIzi/ByK83OS2pNme3cfW/a2HmoWqqlr4eUdByjJz6WsIIepFUWUFeRSVpjDmMJc/3YuIYM9zZ3s3N/JruYO1u1s5vFX6+iL2Xlcmp/Dh0+fxsfOnJl1F51IREdPHx/76Rr+srOZ269cxLvnj55DHgejQyFFMlB/xFHf0sWu/R3sPNDJs6/t47FX6wiHjPctmsJN58xi9vjSoZ8oC3T29POxn61h9RtNfO+KRVyycHLQJY0odT8gMsrsaOrgxy9s45drd9LVG+GCE8bz8XOPY/H0sVl7kldXbz833ruWF2ob+c7lC3nfouqgSxpxCneRUWp/ew/3/nk796zczoGOXk6ZVs7Hzz2Od58wIWuO72/p6mX5hnruX72D9bua+fYHTuZD/jkQ2U7hLjLKdfb08/BLO7nr+W3s3N9JSX4OIQPnIOJfkCXiXMx9R044REVRHhXFeYwr8f5XFOcxrjiPiuJ8KopzmT6umOMnlKb8i6Krt59nNu/jUf9KYz19EaaPK+JzF8zlskVTUlpLkHSGqsgoV5gX5pp3zODDS6bx+IZ61mzfT8gMMw712x+9HfL/d/dF2N/ew/72Hprae9je1M7+th7ae/oPe+7KkjzOOK6Ss+ZUctbsSiaP0CUc+/oj/HlbE//zlz0s31BPa3cflSX5fOS0aSx72xQWVo/J2ianY6VwF8lyOeEQlyycfEw7Grt6+znQ0UNTWw81dS2sqG3khdomHl23B4BZVcWcNdsL+tOPG0dZQWKn+jvnaO7opaGtm4bWt/72tXaxt6WblVubaGzrpiQ/h6UnTmTZ2ybzjlnjyFF30UNSs4yIHBXnHK/tbeWF1xt5obaR1dv209nbTzhkTBpTgBkY5v/3DjH1/nu327v7aGzrPtQBXKz8nBBVpfmcOHkMy942mfPnjddJXD41y4jIiDIzv1+eMm44exY9fRFe3nGAFbWN7D7QicP7AvD+E3ffUZibw/iyfKpK8qkqPfyvND9HzS3HSOEuIkmRlxPi9FnjOH3WuKBLEUANVyIiWUjhLiKShRTuIiJZSOEuIpKFFO4iIlkooXA3s6Vm9pqZ1ZrZLQOMzzezX/rjV5vZjGQXKiIiiRsy3M0sDNwBXATMB640s/lxk10PHHDOzQa+C3wr2YWKiEjiEtlyXwLUOue2Oed6gAeBZXHTLAPu8W//CniX6QwEEZHAJHIS0xRgZ8z9XcBpg03jnOszs4PAOKAxdiIzuwm4yb/bZmavDfB6lfGPywCqOTUyreZMqxdUc6ocS83TE5kokXAfaAs8vjOIRKbBOXcncOcRX8xsbSL9JqQT1ZwamVZzptULqjlVUlFzIs0yu4DYXvCrgT2DTWNmOcAYYH8yChQRkeFLJNzXAHPMbKaZ5QFXAI/GTfMocK1/+4PA0y6o7iZFRGToZhm/Df1mYDkQBn7inNtoZrcCa51zjwI/Bu4zs1q8LfYrjqGmIzbbpCnVnBqZVnOm1QuqOVVGvObA+nMXEZGRozNURUSykMJdRCQLpTTcE+jG4B/MbJOZrTezP5rZ9Jhx/Wb2F/8vfodukDVfZ2YNMbXdEDPuWjN73f+7Nv6xAdX73Zhat5hZc8y4oObxT8xsn5ltGGS8mdl/+e9pvZmdEjMuiHk8VL0f8etcb2YrzWxhzLjtZvaqP49Tdp3JBGo+z8wOxnz+X40Zd8RlKsCavxhT7wZ/+a3wx6V8PpvZVDN7xsxqzGyjmX1mgGlStyw751Lyh7czdiswC8gD1gHz46Y5Hyjyb38S+GXMuLZU1TrMmq8Dbh/gsRXANv//WP/22KDrjZv+U3g7yAObx/7rngOcAmwYZPzFwON451OcDqwOah4nWO8Z0Trwuu1YHTNuO1CZhvP4POD3x7pMpbLmuGkvwTtKL7D5DEwCTvFvlwJbBsiLlC3LqdxyH7IbA+fcM865Dv/uKrxj6oOUSNcLg7kQ+INzbr9z7gDwB2DpCNUZNdx6rwQeGOGahuSce44jnxexDLjXeVYB5WY2iWDm8ZD1OudW+vVAeizHiczjwRzLOnBMhllz4Muyc67OOfeyf7sVqME7ez9WypblVIb7QN0YxL/xWNfjfcNFFZjZWjNbZWaXjUSBA0i05g/4P7F+ZWbRE76G+36TIeHX9Ju8ZgJPxwwOYh4nYrD3FcQ8Hq745dgBT5rZS+Z1x5FO3mFm68zscTNb4A9L+3lsZkV4QfjrmMGBzmfzesZdBKyOG5WyZTmVF8hOqIsCADO7ClgMnBszeJpzbo+ZzQKeNrNXnXNbR6DOw0oZYFh8zb8DHnDOdZvZJ/A6UHtngo9NtuG85hXAr5xz/THDgpjHiRjsfQUxjxNmZufjhftZMYPP9OfxeOAPZrbZ30IN2svAdOdcm5ldDDwCzCHN57HvEmCFcy52Kz+w+WxmJXhfNJ91zrXEjx7gISOyLKdyyz2RbgwwswuArwCXOue6o8Odc3v8/9uAZ/G+FUfakDU755pi6rwLODXRx46A4bzmFcT9jA1oHidisPcVxDxOiJmdDNwNLHPONUWHx8zjfcBv8Zo9Aueca3HOtfm3HwNyzaySNJ7HMY60LKd0PptZLl6w/9w595sBJkndspzCnQ05eDsJZvLWjpkFcdMswtt5Mydu+Fgg379dCbxOCnbqJFjzpJjb7wNWubd2kLzh1z7Wv10RdL3+dMfj7XCyoOdxzOvPYPCdfe/l8J1QLwY1jxOsdxpQC5wRN7wYKI25vRJYmibzeGJ0ecALwh3+/E5omQqiZn98tB+r4qDnsz+/7gVuO8I0KVuWU/IBxbyxi/H2IG8FvuIPuxVvKx3gKWAv8Bf/71F/+BnAq/6C9SpwfRrV/O/ARr+2Z4B5MY/9mL+S1wIfTYd6/ftfA74Z97gg5/EDQB3Qi7cFcz3wCeAT/njDu2DMVr+2xQHP46HqvRs4ELMcr/WHz/Ln7zp/mflKGs3jm2OW41XEfDENtEylQ83+NNcBD8Y9LpD5jNf85oD1MZ/9xUEty+p+QEQkC+kMVRGRLKRwFxHJQgp3EZEspHAXEclCCncRkSykcBcRyUIKdxGRLPT/AW4qI8Qa+LSqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x210e4a8c5f8>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(test_lamdas[5:40], variances[5:40])\n",
    "plt.title(\"Variance of MC at Different Lambda Values\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perform Importance Sampling Monte Carlo with Optimized Weight Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def importance_sampling_MC(lamda, num_samples):\n",
    "    A = lamda\n",
    "    \n",
    "    running_total = 0\n",
    "    for i in range(num_samples):\n",
    "        r = get_rand_number(0,1)\n",
    "        running_total += f_of_x(inverse_G_of_r(r, lamda=lamda))/g_of_x(inverse_G_of_r(r, lamda=lamda), A, lamda)\n",
    "    approximation = float(running_total/num_samples)\n",
    "    return approximation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Importance Sampling Approximation: 0.6983629517447016\n",
      "Variance: 0.0441790662875613\n",
      "Error: 0.002101881687620911\n"
     ]
    }
   ],
   "source": [
    "# run simulation\n",
    "num_samples = 10000\n",
    "approx = importance_sampling_MC(optimal_lamda, num_samples)\n",
    "variance = get_IS_variance(optimal_lamda, num_samples)\n",
    "error = (variance/num_samples)**0.5\n",
    "\n",
    "# display results\n",
    "print(f\"Importance Sampling Approximation: {approx}\")\n",
    "print(f\"Variance: {variance}\")\n",
    "print(f\"Error: {error}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}