Monte Carlo and Importance Sampling.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import seaborn as sns\n",
    "from tqdm import trange\n",
    "from scipy.stats import norm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# plotting params\n",
    "%matplotlib inline\n",
    "plt.rcParams['font.size'] = 10\n",
    "plt.rcParams['axes.labelsize'] = 10\n",
    "plt.rcParams['axes.titlesize'] = 10\n",
    "plt.rcParams['xtick.labelsize'] = 8\n",
    "plt.rcParams['ytick.labelsize'] = 8\n",
    "plt.rcParams['legend.fontsize'] = 10\n",
    "plt.rcParams['figure.titlesize'] = 12\n",
    "plt.rcParams['figure.figsize'] = (15.0, 8.0)\n",
    "sns.set_style(\"white\")\n",
    "\n",
    "# path params\n",
    "plot_dir = './plots/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question Prompt\n",
    "\n",
    "Given the following current equation\n",
    "\n",
    "$$I(\\Delta L, \\Delta V_{TH}) = \\frac{50}{0.1 + \\Delta L} (0.6 - \\Delta V_{TH})^2$$\n",
    "\n",
    "* $\\Delta L \\sim \\ N(0, 0.01^2)$\n",
    "* $\\Delta V_{TH} \\sim \\ N(0, 0.03^2)$\n",
    "\n",
    "We would like to calculate $P(I > 275)$ using direct Monte-Carlo and Importance Sampling."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Direct Monte-Carlo Estimation\n",
    "\n",
    "In MC estimation, we approximate an integral by the sample mean of a function of simulated random variables. In more mathematical terms,\n",
    "\n",
    "$$\\int p(x)\\ f(x)\\ dx = \\mathbb{E}_{p(x)} \\big[\\ f(x) \\big] \\approx \\frac{1}{N} \\sum_{n=1}^{N}f(x_n)$$\n",
    "\n",
    "where $x_n \\sim \\ p(x)$.\n",
    "\n",
    "A useful application of MC is probability estimation. In fact, we can cast a probability as an expectation using the indicator function. In our case, given that $A = \\{I \\ | \\ I > 275\\}$, we define $f(x)$ as\n",
    "\n",
    "$$f(x) = I_{A}(x)= \\begin{cases} \n",
    "      1 & I \\geq 275 \\\\\n",
    "      0 & I < 275 \n",
    "   \\end{cases}$$\n",
    "   \n",
    "Replacing in our equation above, we get\n",
    "\n",
    "$$\\int p(x) \\ f(x) \\ dx = \\int I(x)\\ p(x) \\ d(x) = \\int_{x \\in A} p(x)\\ d(x) \\approx \\frac{1}{N} \\sum_{n=1}^{N}I_{A}(x_n)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def monte_carlo_proba(num_simulations, num_samples, verbose=True, plot=False):\n",
    "    \n",
    "    if verbose:\n",
    "        print(\"===========================================\")\n",
    "        print(\"{} Monte Carlo Simulations of size {}\".format(num_simulations, num_samples))\n",
    "        print(\"===========================================\\n\")\n",
    "    \n",
    "    num_samples = int(num_samples)\n",
    "    num_simulations = int(num_simulations)\n",
    "    \n",
    "    probas = []\n",
    "    for i in range(num_simulations):\n",
    "        mu_1, sigma_1 = 0, 0.01\n",
    "        mu_2, sigma_2 = 0, 0.03\n",
    "\n",
    "        length = np.random.normal(mu_1, sigma_1, num_samples)\n",
    "        voltage = np.random.normal(mu_2, sigma_2, num_samples)\n",
    "\n",
    "        num = 50 * np.square((0.6 - voltage))\n",
    "        denum = 0.1 + length\n",
    "        I = num / denum\n",
    "        \n",
    "        true_condition = np.where(I >= 275)\n",
    "        false_condition = np.where(I < 275)\n",
    "        num_true = true_condition[0].shape[0]\n",
    "        proba = num_true / num_samples\n",
    "        probas.append(proba)\n",
    "        \n",
    "        if plot:\n",
    "            if i == (num_simulations - 1):\n",
    "                plt.scatter(length[true_condition], voltage[true_condition], color='r')\n",
    "                plt.scatter(length[false_condition], voltage[false_condition], color='b')\n",
    "                plt.xlabel(r'$\\Delta L$ [$\\mu$m]')\n",
    "                plt.ylabel(r'$\\Delta V_{TH}$ [V]')\n",
    "                plt.title(\"Monte Carlo Estimation of P(I > 275)\")\n",
    "                plt.grid(True)\n",
    "                plt.savefig(plot_dir + 'monte_carlo_{}.pdf'.format(num_samples), format='pdf', dpi=300)\n",
    "                plt.show()\n",
    "    \n",
    "    \n",
    "    mean_proba = np.mean(probas)\n",
    "    std_proba = np.std(probas)\n",
    "    \n",
    "    if verbose:\n",
    "        print(\"Probability Mean: {:0.5f}\".format(mean_proba))\n",
    "        print(\"Probability Std: {:0.5f}\".format(std_proba))\n",
    "    \n",
    "    return probas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===========================================\n",
      "10 Monte Carlo Simulations of size 10000\n",
      "===========================================\n",
      "\n",
      "Probability Mean: 0.00196\n",
      "Probability Std: 0.00056\n"
     ]
    }
   ],
   "source": [
    "probas = monte_carlo_proba(10, 10000, plot=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MC_histogram(num_samples, plot=True):\n",
    "    \n",
    "    num_samples = int(num_samples)\n",
    "    \n",
    "    mu_1, sigma_1 = 0, 0.01\n",
    "    mu_2, sigma_2 = 0, 0.03\n",
    "\n",
    "    length = np.random.normal(mu_1, sigma_1, num_samples)\n",
    "    voltage = np.random.normal(mu_2, sigma_2, num_samples)\n",
    "\n",
    "    num = 50 * np.square((0.6 - voltage))\n",
    "    denum = 0.1 + length\n",
    "    I = num / denum\n",
    "    \n",
    "    if plot:\n",
    "        n, bins, patches = plt.hist(I, 50, normed=1, facecolor='green', alpha=0.75)\n",
    "        plt.ylabel('Number of Samples')\n",
    "        plt.xlabel(r'$I_{DS}$ [$\\mu$A]')\n",
    "        plt.title(\"Monte Carlo Estimation of P(I > 275)\")\n",
    "        plt.grid(True)\n",
    "        plt.savefig(plot_dir + 'mc_histogram_{}.pdf'.format(num_samples), format='pdf', dpi=300)\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4AAAAHtCAYAAABIwRD6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3X2YV2W9L/73KMxQyIAkOKOgpZJa\nYyIaZU9H2obUttplaGwkNNOSLaVWkog2Cm5HO5hmSZGGHDXxISVLz2l33D6046i7MatxqyVHLfvJ\nwITS+AAzMOv3R5dzJMQBne+MuF6v6+pqZs19r/VZc99+5/vmXt+1qoqiKAIAAMDr3nb9XQAAAAB9\nQwAEAAAoCQEQAACgJARAAACAkhAAAQAASkIABAAAKAkBEKDE7rnnnuy999655ZZbNtr+0Y9+NF/7\n2te2en8PP/xw/vM//3OL2//v//2/M23atEybNi2TJ0/O//pf/2uL+15yySW55pprtqjttGnT8qlP\nfar7WNOmTcsdd9yx2fbXXnttOjs78+CDD+bb3/72Ftf0Un7+85+ntbU1q1atSmNj46va1+aceuqp\nOeKII7J8+fLubZdcckkOO+yw7vP99Kc/nXvuuaf75/PmzcuKFSty44035r//9/++1cd88MEH88//\n/M+ZNm1ajjvuuLS1teXBBx/c6He833775a677srTTz+dd73rXd3bFy9enKIoMmvWrKxdu7ZXfgcA\nbJkB/V0AAP1rjz32yC233JJ//Md/TPK3EPf888+/on3927/9W3baaae8853v7LHtfffdlyuuuCLf\n+973Mnjw4Dz11FM56qijstdee2WvvfZ6Rcd/Oeeff3723HPPLWr7ve99L//0T/+UfffdN/vuu++r\nOu7/+B//I42Njdlzzz0rFgCXLVuWu+++e5PtxxxzTKZMmZIkWb58eb7yla/kpptuyv33358BAwak\nrq7uJff3u9/9LmvWrMl73/veVFVVvWSbc889N2eeeWb23XffLFmyJN///vdz+umn58orr0yS/M//\n+T8zcuTIfOADH8iyZcty+OGH58wzz9xoH4cffnguu+yynHTSSa/m9AHYCgIgQMnts88+efTRR9Pe\n3p4hQ4bk5ptvzkc/+tE8+eSTSZKbb745ixcvTnV1dd785jfnnHPOyU9+8pPceeedWbt2bf74xz/m\n+OOPz3vf+97cdNNNGThwYN7+9rdn7dq1+eY3v5ntt98+o0ePzjnnnJOBAwd2H/f666/P9OnTM3jw\n4CTJjjvumOuvvz61tbVZsWJFGhsbs27duqxatSonn3xyDj300Bx++OF585vfnIEDB2aPPfbo3ldT\nU1Oam5uT/C1UTJ8+fYvOffXq1Tn55JNTFEXWrVuXs88+Oy0tLVm1alVOOeWUTJ8+PUuWLMk3v/nN\nfOhDH8oBBxyQxx57LAcffHDa29vz29/+Nm95y1vyjW98I7///e/T1NSUDRs25KmnnkpjY2P++te/\n5sEHH8ysWbPyjW98I7Nmzcp1112XX/7yl7noootSU1OTYcOG5V//9V/z4IMP5vvf/34GDhyYJ554\nIh/5yEdy4oknblTvS/W78MIL88wzz+TEE0/MggULNnuuTz/9dN74xjcmSa688soce+yxm207cuTI\n3HrrrbnkkksyYcKEfPKTn8zIkSM3anPhhRd2b9uwYUNqamq6f/bcc8/lkksuyVVXXZUkaWlpyQMP\nPJCjjz46w4cPz5w5czJy5Mi85z3vSVNTU2bMmJHttnNREkCfKAAorbvvvrs4+eSTi+985zvFDTfc\nUHR1dRVHH310ceeddxazZs0qVq9eXRx66KFFe3t7URRFce655xZXXnll8aMf/aj47Gc/WxRFUTz6\n6KPFYYcdVhRFUXzrW98qfvjDHxZdXV3FxIkTi7a2tqIoiuKb3/xmce2112507M9+9rPFgw8++JJ1\n/fKXvyzuvvvuoiiKorm5uTjmmGOKoiiKCRMmFA888MBGx/r3f//34l/+5V+Krq6uoqOjo/jUpz5V\nPPTQQxvt7+ijjy6OOOKI4uijj+7+31/+8pfi9ttvL2bOnFk8//zzxe9+97viV7/6Vfdx1q5d2/37\nKYqi2HfffYs///nPRUdHRzF27NjiD3/4Q9HV1VVMmDChWLNmTXHLLbd0H/fmm28uzjjjjO5jP/LI\nI8Wf/vSnYvLkyd19VqxYURRFUVxxxRVFU1NTcffddxcf/vCHi87OzuLZZ58txo0bt9E5bK5fURTF\ne97znk1+h9/61reKiRMnFkcffXTxmc98pjjppJOK5cuXF0VRFO973/uKdevWFUVRFD/60Y+Kb3zj\nGy85Dp2dncXSpUuL8ePHF7fddttLtmlubi4mTZpU/OUvf+netnjx4uLiiy/u/v7nP/958ctf/rIo\niqL48Y9/XMycOXOjsdncPACg91kBBCAf/ehH09jYmNGjR+eggw7q3v6nP/0pe+21V3bYYYckyTvf\n+c78x3/8R/bff//ss88+SZL6+vp0dHRstL/Vq1dn5cqVOfnkk5Mka9euzXve856N2uyyyy558skn\nu/eTJM3Nzdlpp50yYsSILFiwIDfccEOqqqqyfv367jZvectbNtrP8uXLc9BBB6WqqioDBw7M/vvv\nn+XLl2fvvffeqN1LXQL6gQ98II899lhmzJiRAQMGbLLi9mLDhg3LLrvskiR54xvf2H2Z6pAhQ7Ju\n3bqMHDkyl156aQYNGpRnn322+3f295566qnssMMO2Xnnnbt/pxdeeGEOOeSQvPWtb82AAQMyYMCA\nDBo0aIv6vZwXXwL6Yl1dXamurn7Zvi0tLbnuuuvy6KOP5tRTT8273/3uTdrceuutWbBgQRYuXJjh\nw4d3b//JT36Sb33rW93fv/vd784b3vCGJMmHPvShjX42cuTIPP300y9bCwC9x/UWAGT06NF57rnn\ncuWVV+ZjH/tY9/ZRo0Zl+fLlee6555Ik9957b3cAe6nPhlVVVaWrqys77rhj6urqcumll+bKK6/M\nF77whU0CxCc/+clcfvnl3fv+y1/+ktmzZ+f555/PxRdfnI9//OP5xje+kXe9610piqK7399fKrjn\nnnt2X/7Z2dmZX//619l999236LzvueeejBw5Mj/4wQ9y4okndgeqF87j78/t5Zx77rn54he/mPPP\nPz9vfetbu2uuqqraqP4dd9wxzzzzTFauXJnkb7/TN7/5zT0e4+X6ba2ampps2LBhsz//93//91x7\n7bU58sgjc+WVV+aoo47qvnz0BT/+8Y9z1VVX5corr8zo0aO7t7e3t6ejoyP19fXd2+bMmZOf/exn\nSZL/83/+T97+9rd3/2zNmjV505ve9IrOA4CtZwUQgCTJRz7ykfz4xz/OW97ylvzpT39KkgwfPjwz\nZ87MZz7zmWy33XbZbbfd8pWvfGWTu4a+oKGhIRdccEH23HPPnHHGGTnhhBNSFEUGDx6cCy64YKO2\nBxxwQI488sh89rOfzYABA7J27dqceuqp2WeffTJp0qRccMEFWbhwYerq6vLUU09ttu4JEybk3nvv\nzVFHHZXOzs5MmjRpo4DxglmzZnWvQiXJhz/84Xz4wx/OqaeemmuuuSbr16/Pv/zLvyRJDjrooJxw\nwgnd32+Jj33sY/nSl76U2trajWo+4IADctppp2Xu3LlJ/hby5s2bl5kzZ6aqqipDhw7Neeedlz/8\n4Q8vu//N9Xslxo0blwceeCDveMc7XvLnH/zgB/PBD35ws/03bNiQc889N/X19Zk5c2aSv61IfvGL\nX8yjjz6aXXfddaP2X/7ylzN79uxcc801ecMb3pB58+Yl+dtKZGtra0Vu+gPAS6sqXvzPkgDA696v\nf/3r3HLLLZkzZ06/1nHnnXfmgQceyIwZM/q1DoAycQkoAJTMAQcckA0bNmTFihX9VkNRFPnJT36S\nY445pt9qACgjK4AAAAAlYQUQAACgJARAAACAkthm7wK6du3atLS0ZMSIEdl+++37uxwAAIA+tWHD\nhqxatSoNDQ2bPD92c7bZANjS0pKpU6f2dxkAAAD96uqrr85BBx20RW232QA4YsSIJH872bq6un6u\npn+0tLSkoaGhv8ugnxh/zAHMgXIz/pgD5dbS0pKddtopU6dO7c5GW2KbDYAvXPZZV1eXUaNG9XM1\n/aO1tbW0547xxxzAHCg74485UG6tra3Zeeedk2SrPhLnJjAAAAAlIQACAACUhAAIAABQEgIgAABA\nSQiAAAAAJSEAAgAAlIQACAAAUBICIAAAQEkIgAAAACUhAAIAAJSEAAgAAFASAiAAAEBJDKjETltb\nW9PU1JShQ4dmzJgxmTp1apJk2bJlWbp0aYqiyJQpUzJu3Li0t7dn9uzZOe644zJ27Njcddddue22\n27Jhw4YccsghOfTQQytRIgAAQOlUZAVwyZIlmTZtWhobG3PHHXeks7MzSbJo0aLMmzcvc+fOzcKF\nC5MkCxYsyKBBg7r7Xn/99Rk2bFjWr1+ft73tbZUoDwAAoJQqsgLY1taW+vr6JEltbW3a29szfPjw\nFEWR6urqJOkOhaeddlouueSS7r4PPPBALrjggrS1teWiiy7KBRdc8LLHamlpSWtrayVOY5vQ3Nzc\n3yXQj4w/5gDmQLkZf8yBcmtpadnqPhUJgPX19VmxYkXq6+uzZs2a1NbWJklqamrS0dGRrq6u7iD4\n93bdddfU1NRk2LBhW3SshoaGjBo1qtdq35Y0NzfnwAMP7O8y6CfGH3MAc6DcjD/mQLk1NzenoaFh\nq/tVJABOnjw5TU1NuemmmzJx4sScd955mTVrVqZPn545c+Zk/fr1mTFjxkv2nTp1ambNmpX169fn\nxBNPrER5AAAApVSRADhixIjMnz9/k+3jx4/P+PHjN9k+c+bM7q8nTZqUSZMmVaIsAACAUqtIAATo\nbRMWT3jFfW+ffnsvVgIAsO3yHEAAAICSEAABAABKQgAEAAAoCQEQAACgJARAAACAkhAAAQAASkIA\nBAAAKAkBEAAAoCQEQAAAgJIY0N8FAFTahMUTXnHf26ff3ouVAAD0LyuAAAAAJSEAAgAAlIRLQIE+\n8WouwwQAoHdYAQQAACgJARAAAKAkBEAAAICSEAABAABKQgAEAAAoCQEQAACgJARAAACAkhAAAQAA\nSkIABAAAKAkBEAAAoCQEQAAAgJIQAAEAAEpCAAQAACgJARAAAKAkBEAAAICSEAABAABKQgAEAAAo\nCQEQAACgJARAAACAkhAAAQAASkIABAAAKAkBEAAAoCQEQAAAgJIQAAEAAEpCAAQAACgJARAAAKAk\nBEAAAICSEAABAABKQgAEAAAoiQH9XQDAa9mExRNecd/bp9/ei5UAALx6VgABAABKQgAEAAAoiYoE\nwNbW1pxyyilpbGzM1Vdf3b192bJlOe200/LVr3419913X5Kkvb09M2fOzP3339/drr29PR/5yEey\natWqSpQHAABQShUJgEuWLMm0adPS2NiYO+64I52dnUmSRYsWZd68eZk7d24WLlyYJFmwYEEGDRrU\n3berqysXXnhhdtttt0qUBgAAUFoVuQlMW1tb6uvrkyS1tbVpb2/P8OHDUxRFqqurk6Q7FJ522mm5\n5JJLuvt++9vfzlFHHZUrrrhii47V0tKS1tbW3j2BbUhzc3N/l0A/2pbGv729vb9L6HN9MT7b0hyg\nMsyBcjP+mAPl1tLSstV9KhIA6+vrs2LFitTX12fNmjWpra1NktTU1KSjoyNdXV3dQfDFVq9end/8\n5jdpa2vLr3/961x22WU5/fTTX/ZYDQ0NGTVqVCVO4zWvubk5Bx54YH+XQT/Z1sZ/SMuQ/i6hz1V6\nfLa1OUDvMwfKzfhjDpRbc3NzGhoatrpfRQLg5MmT09TUlJtuuikTJ07Meeedl1mzZmX69OmZM2dO\n1q9fnxkzZmzSb/jw4bn88suTJF/72tfyuc99rhLlAQAAlFJFAuCIESMyf/78TbaPHz8+48eP32T7\nzJkzN9nW1NRUidIAAABKy4PggS32ah6KDgBA//McQAAAgJIQAAEAAEpCAAQAACgJARAAAKAkBEAA\nAICSEAABAABKQgAEAAAoCQEQAACgJARAAACAkhAAAQAASkIABAAAKAkBEAAAoCQEQAAAgJIQAAEA\nAEpCAAQAACgJARAAAKAkBEAAAICSEAABAABKQgAEAAAoCQEQAACgJARAAACAkhAAAQAASkIABAAA\nKAkBEAAAoCQEQAAAgJIQAAEAAEpCAAQAACgJARAAAKAkBEAAAICSEAABAABKQgAEAAAoiQH9XQDA\n69WExRNecd/bp9/ei5UAAPyNFUAAAICSEAABAABKQgAEAAAoCQEQAACgJARAAACAkhAAAQAASkIA\nBAAAKAkBEAAAoCQEQAAAgJIQAAEAAEpCAAQAACgJARAAAKAkBEAAAICSGFCJnba2tqapqSlDhw7N\nmDFjMnXq1CTJsmXLsnTp0hRFkSlTpmTcuHFpb2/P7Nmzc9xxx2Xs2LH54Q9/mN///vd57rnnMmnS\npHzwgx+sRIkAAAClU5EVwCVLlmTatGlpbGzMHXfckc7OziTJokWLMm/evMydOzcLFy5MkixYsCCD\nBg3q7ltbW5vGxsacfvrp+elPf1qJ8gAAAEqpIiuAbW1tqa+vT/K3QNfe3p7hw4enKIpUV1cnSXco\nPO2003LJJZd09z388MPz7LPP5vzzz88JJ5zQ47FaWlrS2tpagbPYNjQ3N/d3CfSjvh7/9vb2Pj1e\nmW3p2HoNwBwoN+OPOVBuLS0tW92nIgGwvr4+K1asSH19fdasWZPa2tokSU1NTTo6OtLV1dUdBP/e\n8uXLs2DBgnzpS1/K6NGjezxWQ0NDRo0a1av1byuam5tz4IEH9ncZ9JP+GP8hLUP69HhltiVj6zUA\nc6DcjD/mQLk1NzenoaFhq/tV5BLQyZMn56qrrspZZ52ViRMn5rzzzktHR0emT5+eOXPmZPbs2Zkx\nY8Ym/YqiyBe+8IWsW7cuF198cfdlogAAALx6FVkBHDFiRObPn7/J9vHjx2f8+PGbbJ85c2b31z//\n+c8rURIAAEDpeQwEAABASQiAAAAAJSEAAgAAlIQACAAAUBICIAAAQElU5C6gwGvXhMUT+rsEAAD6\niRVAAACAkhAAAQAASkIABAAAKAkBEAAAoCQEQAAAgJIQAAEAAEpCAAQAACgJARAAAKAkBEAAAICS\nEAABAABKQgAEAAAoCQEQAACgJARAAACAkhAAAQAASkIABAAAKAkBEAAAoCQEQAAAgJIQAAEAAEpC\nAAQAACgJARAAAKAkBEAAAICSEAABAABKYkB/FwDApiYsntBjm/b29gxpGfKSP7t9+u29XRIA8Dpg\nBRAAAKAkBEAAAICSEAABAABKQgAEAAAoCQEQAACgJARAAACAkhAAAQAASkIABAAAKAkBEAAAoCQE\nQAAAgJIQAAEAAEpCAAQAACgJARAAAKAkBEAAAICSEAABAABKQgAEAAAoCQEQAACgJAZUYqetra1p\namrK0KFDM2bMmEydOjVJsmzZsixdujRFUWTKlCkZN25c2tvbM3v27Bx33HEZO3bsS7YBAADg1atI\nAFyyZEmmTZuWcePG5fjjj8+RRx6ZgQMHZtGiRfnOd76Trq6unHzyyfnud7+bBQsWZNCgQd19X6oN\nAAAAr15FLgFta2tLfX19kqS2tjbt7e1JkqIoUl1dnUGDBqWzszNJctppp2W33Xbr7vtSbQAAAHj1\ntmoF8Mknn+wOdi+nvr4+K1asSH19fdasWZPa2tokSU1NTTo6OtLV1ZXq6uqX7LslbV6spaUlra2t\nW3MaryvNzc39XQL96JWM/wv/IMPrw+bG02tDeRjrcjP+mAPl1tLSstV9egyAl112WWpra/PXv/41\nN954Y97//vfn9NNPf9k+kydPTlNTU2666aZMnDgx5513XmbNmpXp06dnzpw5Wb9+fWbMmPGSfbek\nzYs1NDRk1KhRPbZ7PWpubs6BBx7Y32XQT17p+A9pGVKBaugP7e3tGTLkpcfTa0M5+DtQbsYfc6Dc\nmpub09DQsNX9egyA//Zv/5arrroqn/vc53LrrbfmM5/5TI87HTFiRObPn7/J9vHjx2f8+PGbbJ85\nc2aPbQAAAHh1evwM4HbbbZe2trbstNNOSZK1a9dWvCgAAAB6X48B8F3velemTZuWo48+Ov/6r/+a\n//bf/ltf1AUAAEAv6/ES0FNOOSWnnHJKnn766XzlK1/ZohuzAAAA8NrTYwD8z//8z5x99tnZsGFD\nJk2alF122SWTJ0/ui9oAAADoRT1eAnrRRRflqquuyk477ZQvfOELueaaa/qiLgAAAHrZFt0EZtiw\nYamqqkpNTU0GDx7cF3UBAADQy3oMgLvttlvmz5+fp59+OgsXLswuu+zSF3UBAADQy3oMgGeffXZ2\n2WWXHHjggXnjG9+YuXPn9kVdAAAA9LLN3gTmP/7jP7q/Hj16dEaPHp0kuffee/O+972v8pUBAADQ\nqzYbAG+55ZbNdhIAAQAAtj2bDYDnnXde99f/9V//lUcffTR77bVX9t577z4pDAAAgN7V43MAL7ro\notx99915xzvekSuvvDKHHnpoPve5z/VFbQAAAPSiHgPgXXfdlRtuuCHbbbddNmzYkKOOOkoABAAA\n2Ab1eBfQurq6PPvss0mS9evXZ6eddqp4UQAAAPS+HlcAV65cmcMOOyz77LNPHnnkkQwcODCf/vSn\nkyRLliypeIEAAAD0jh4D4MUXX9wXdQAAAFBhPQbAv/zlL7nllluybt267m2NjY2VrAkAAIAK6DEA\nzpo1K8cff3xqa2v7oh4AAAAqpMcAuPvuu+eTn/xkX9QCAABABfUYAA877LCccsop2XPPPbu3nXTS\nSRUtCnh5ExZPSHt7e4a0DOnvUgAA2Ib0GACvvvrqTJw40SWgAAAA27geA+CwYcNywgkn9EUtAAAA\nVFCPAXDHHXfMWWedlbe97W2pqqpKkhx11FEVLwwAAIDetUU3gUmStra2ihcDAABA5fQYAE866aSs\nXLky69evT1EUWblyZV/UBQAAQC/rMQDOnj07999/f55//vmsXbs2o0ePznXXXdcXtQEAANCLtuup\nwUMPPZRbbrkl73vf+3LLLbekpqamL+oCAACgl/UYAHfcccdUVVXlueeey/Dhw/uiJgAAACqgx0tA\n3/72t+fyyy/PyJEjc8opp2Tt2rV9URcAr8KExRNecd/bp9/ei5UAAK8lPQbAU089Nc8880wGDRqU\nu+66K+94xzv6oi4AAAB62WYvAV23bl0WL16coijyzDPP5NRTT82tt96aoij6sj4AAAB6yWYD4Lx5\n8/L//X//X7q6unL22Wdnn332yWGHHZbGxsY+LA8AAIDestlLQP/whz9kyZIlWbduXZqbm/Otb30r\nAwcOzA9+8IO+rA8AAIBestkVwMGDBydJ7rvvvuy3334ZOHBgkr9dGgoAAMC2Z7MrgIMHD861116b\nn/3sZzn88MPT1dWVm2++OfX19X1ZHwAAAL1ksyuAjY2N+eMf/5j3v//9+cQnPpF77rknP/vZz3wG\nEAAAYBu12RXA4cOH56tf/Wr39wcffHAOPvjgPikKAACA3rfZFUAAAABeXzYbANvb2/uyDgAAACps\nswHwhBNOSJJ8/etf77NiAAAAqJzNfgZwwIABOeKII/L444/n4YcfTpIURZGqqqosWbKkzwoEAACg\nd2w2AF5xxRVpbW1NY2NjGhsbUxRFX9YFAABAL9tsANx+++2zyy675NJLL821116bRx55JG9+85sz\nZcqUvqwPAACAXtLjXUDPOuus/PGPf8x73/ve/PnPf86cOXP6oi4AAAB62WZXAF/w+OOP5+qrr06S\nHHroofn0pz9d8aIAAADofT2uAK5bty7PP/98kmTt2rXZsGFDxYsCAACg9/W4AviZz3wmH//4xzNm\nzJg88sgj+eIXv9gXdQEAANDLegyAH/vYx/KBD3wgf/rTnzJq1KjsuOOOPe60tbU1TU1NGTp0aMaM\nGZOpU6cmSZYtW5alS5emKIpMmTIlY8eOzZlnnpnBgweno6MjjY2N+cUvfpE777wzHR0dOeCAA/KJ\nT3zi1Z8lAAAAPV8CmiTDhg3Lfvvtt0XhL0mWLFmSadOmpbGxMXfccUc6OzuTJIsWLcq8efMyd+7c\nLFy4MPfcc09Gjx6d2bNnZ/jw4bnvvvvyq1/9Kv/1X/+Vxx57LHV1da/8zAAAANhIjyuAr0RbW1vq\n6+uTJLW1tWlvb8/w4cNTFEWqq6uTJJ2dnWlra+sOeXV1dVm5cmXe85735MQTT8zatWvzta99LQcf\nfPDLHqulpSWtra2VOI1tQnNzc3+XQD9ob2/f6P8pr0rMAa8r2xbjVW7GH3Og3FpaWra6T48B8PLL\nL89xxx23VTutr6/PihUrUl9fnzVr1qS2tjZJUlNTk46OjnR1daW6ujr19fXdk3bFihXZa6+98s1v\nfjOLFi3KkCFDtuiGMw0NDRk1atRW1fd60dzcnAMPPLC/y6AfDGkZkvb29gwZMqS/S6EfVWoOeF3Z\ndvg7UG7GH3Og3Jqbm9PQ0LDV/XoMgHfeeWeOOeaYbL/99lu808mTJ6epqSk33XRTJk6cmPPOOy+z\nZs3K9OnTM2fOnKxfvz4zZsxIQ0NDbr755sybNy9JMm7cuHzyk5/Mqaeemje+8Y055phjtvqEAAAA\neGk9BsCnnnoq73//+zNq1KhUVVWlqqoqS5Ysedk+I0aMyPz58zfZPn78+IwfP36jbeecc85G33/i\nE59w4xcAAIAK6DEAfve73+2LOgAAAKiwHgPggAED8o1vfCOrV6/OpEmTsvfee2fXXXfti9oAAADo\nRT0+BuLMM8/MEUcckc7Ozhx00EE599xz+6IuAAAAelmPAXDt2rU5+OCDU1VVlT322CM1NTV9URcA\nAAC9rMcAWFNTk1/84hfp6urK/fff3/0cPwAAALYtPQbAuXPn5sYbb8xTTz2VH/zgB2lsbOyDsgAA\nAOhtPd4Epq6uLp///Ofz2GOPZcyYMRk9enRf1AUAAEAv6zEAXnrppfnFL36R/fbbL1dccUUmTZrk\nAe0AAADboB4D4J133plrrrmQVGCBAAAWyUlEQVQm2223XdavX59//ud/FgABAAC2QT1+BvBNb3pT\nnn/++SRJZ2dnhg8fXvGiAAAA6H2bXQE86qijUlVVlb/85S857LDDsvfee2f58uUZNmxYX9YHAABA\nL9lsALzwwgv7sg4AAAAqbLMBcNddd02S/Pa3v80tt9ySdevWdf/MoyAAAAC2PT3eBGbWrFk5/vjj\nU1tb2xf1AAAAUCE9BsDdd989n/zkJ/uiFgAAACqoxwB42GGH5ZRTTsmee+7Zve2kk06qaFEAAAD0\nvh4D4NVXX52JEye6BBQAAGAb12MAHDZsWE444YS+qAUAAIAK6jEA7rjjjjnrrLPytre9LVVVVUn+\n9oxAAAAAti1bdBOYJGlra6t4MQAAAFROjwHQHUABAABeH3oMgKecckqqqqrS1dWVJ554Irvvvnuu\nueaavqgNAACAXtRjALz22mu7v/7rX/+aM888s6IFAQAAUBnbbU3jIUOG5E9/+lOlagEAAKCCelwB\nPOqoo1JVVZWiKLJ69eocfPDBfVEXAP1kwuIJr7jv7dNv78VKAIDe1mMAvPDCC7u/rqmpyU477VTR\nggAAAKiMzQbApUuXbrbTP/3TP1WkGAAAACpnswFw+fLlG31fFEVuvPHGDBo0SAAEAADYBm02AH75\ny1/u/vqPf/xjZs2alUMOOSSzZ8/uk8IAAADoXT1+BvDqq6/O4sWLc/rpp2fChFd+YwAAAAD612YD\nYGtra04//fQMHTo0119/fYYOHdqXdQEAANDLNhsA//Ef/zHV1dV597vfnXPOOWejn82fP7/ihQEA\nANC7NhsAL7300r6sAwAAgArbbAAcP358X9YBAABAhW3X3wUAAADQN3q8CyjQ+yYsdkddAAD6nhVA\nAACAkhAAAQAASkIABAAAKAkBEAAAoCQEQAAAgJIQAAEAAEpCAAQAACgJARAAAKAkBEAAAICSEAAB\nAABKYkAldtra2pqmpqYMHTo0Y8aMydSpU5Mky5Yty9KlS1MURaZMmZKxY8fmzDPPzODBg9PR0ZHG\nxsbcddddue2227Jhw4YccsghOfTQQytRIgAAQOlUZAVwyZIlmTZtWhobG3PHHXeks7MzSbJo0aLM\nmzcvc+fOzcKFC3PPPfdk9OjRmT17doYPH5777rsv119/fYYNG5b169fnbW97WyXKAwAAKKWKrAC2\ntbWlvr4+SVJbW5v29vYMHz48RVGkuro6SdLZ2Zm2trbU1dUlSerq6rJy5co88MADueCCC9LW1paL\nLrooF1xwwcseq6WlJa2trZU4jW1Cc3Nzf5fAK9De3v6a2g/brtfaHPCa1Pf8zsvN+GMOlFtLS8tW\n96lIAKyvr8+KFStSX1+fNWvWpLa2NklSU1OTjo6OdHV1pbq6OvX19d2TdsWKFdlrr72y6667pqam\nJsOGDduiYzU0NGTUqFGVOI3XvObm5hx44IH9XQavwJCWIa96H+3t7Rky5NXvh23Xa3EOeE3qW/4O\nlJvxxxwot+bm5jQ0NGx1v4oEwMmTJ6epqSk33XRTJk6cmPPOOy+zZs3K9OnTM2fOnKxfvz4zZsxI\nQ0NDbr755sybNy9JMm7cuEydOjWzZs3K+vXrc+KJJ1aiPAAAgFKqSAAcMWJE5s+fv8n28ePHZ/z4\n8RttO+ecczb6ftKkSZk0aVIlygIAACg1j4EAAAAoCQEQAACgJARAAACAkhAAAQAASkIABAAAKAkB\nEAAAoCQEQAAAgJIQAAEAAEpCAAQAACgJARAAAKAkBEAAAICSEAABAABKQgAEAAAoiQH9XQAArx8T\nFk94xX1vn357L1YCALwUK4AAAAAlIQACAACUhAAIAABQEgIgAABASQiAAAAAJSEAAgAAlIQACAAA\nUBICIAAAQEkIgAAAACUhAAIAAJSEAAgAAFASAiAAAEBJCIAAAAAlIQACAACUhAAIAABQEgIgAABA\nSQiAAAAAJSEAAgAAlIQACAAAUBICIAAAQEkIgAAAACUhAAIAAJSEAAgAAFASAiAAAEBJCIAAAAAl\nIQACAACUhAAIAABQEgIgAABASQiAAAAAJSEAAgAAlMSA/i4AtlUTFk/o7xIAAGCrWAEEAAAoiYqs\nALa2tqapqSlDhw7NmDFjMnXq1CTJsmXLsnTp0hRFkSlTpmTs2LE588wzM3jw4HR0dKSxsTFJ0t7e\nnqOOOiqLFy/OiBEjKlEiAABA6VRkBXDJkiWZNm1aGhsbc8cdd6SzszNJsmjRosybNy9z587NwoUL\nc88992T06NGZPXt2hg8fnvvuuy9dXV258MILs9tuu1WiNAAAgNKqyApgW1tb6uvrkyS1tbVpb2/P\n8OHDUxRFqqurkySdnZ1pa2tLXV1dkqSuri4rV67Mt7/97Rx11FG54oortuhYLS0taW1trcRpbBOa\nm5v7u4TSam9v7+8SXhM10L9eT3PgoG8f9Kr6f+/g7/VSJdsWfwfKzfhjDpRbS0vLVvepSACsr6/P\nihUrUl9fnzVr1qS2tjZJUlNTk46OjnR1daW6ujr19fXdk3bFihV585vfnN/85jdpa2vLr3/961x2\n2WU5/fTTX/ZYDQ0NGTVqVCVO4zWvubk5Bx54YH+XUVpDWob06/Hb29szZEj/1kD/Mgc2VsbXQ38H\nys34Yw6UW3NzcxoaGra6X0UC4OTJk9PU1JSbbropEydOzHnnnZdZs2Zl+vTpmTNnTtavX58ZM2ak\noaEhN998c+bNm5ckGT9+fMaPH58k+drXvpbPfe5zlSgPAACglCoSAEeMGJH58+dvsv3FAe8F55xz\nzkvuo6mpqRKlAQAAlJbHQAAAAJSEAAgAAFASAiAAAEBJCIAAAAAlIQACAACUhAAIAABQEgIgAABA\nSQiAAAAAJSEAAgAAlIQACAAAUBICIAAAQEkIgAAAACUhAAIAAJSEAAgAAFASAiAAAEBJCIAAAAAl\nIQACAACUhAAIAABQEgIgAABASQiAAAAAJSEAAgAAlIQACAAAUBICIAAAQEkIgAAAACUhAAIAAJTE\ngP4uAAB6w4TFE15x39un396LlQDAa5cVQAAAgJIQAAEAAEpCAAQAACgJARAAAKAkBEAAAICSEAAB\nAABKQgAEAAAoCQEQAACgJARAAACAkhAAAQAASkIABAAAKIkB/V0A9KcJiyf0dwkAANBnrAACAACU\nhAAIAABQEgIgAABASQiAAAAAJSEAAgAAlIQACAAAUBICIAAAQEkIgAAAACVRkQfBt7a2pqmpKUOH\nDs2YMWMyderUJMmyZcuydOnSFEWRKVOmZOzYsTnzzDMzePDgdHR0pLGxMT/84Q/z+9//Ps8991wm\nTZqUD37wg5UoEQC6TVg84RX3vX367b1YCQBUVkVWAJcsWZJp06alsbExd9xxRzo7O5MkixYtyrx5\n8zJ37twsXLgw99xzT0aPHp3Zs2dn+PDhue+++1JbW5vGxsacfvrp+elPf1qJ8gAAAEqpIgGwra0t\n9fX1SZLa2tq0t7cnSYqiSHV1dQYNGpTOzs60tbWlrq4uSVJXV5eVK1fm8MMPz7PPPpvzzz8/J5xw\nQiXKAwAAKKWKXAJaX1+fFStWpL6+PmvWrEltbW2SpKamJh0dHenq6kp1dXXq6+vT3NycJFmxYkX2\n2muvLF++PAsWLMiXvvSljB49usdjtbS0pLW1tRKnsU144ffHK/PCP05sq7b1+nn1zIH+19+vw/19\nfPqX8cccKLeWlpat7lNVFEXR24WsWrUqTU1NGTx4cBoaGvLwww9n1qxZuf/++3PDDTdk/fr1OfbY\nY9PQ0JCvf/3rqa6uTpKcccYZmThxYvbZZ5/U1NTkrW9962ZXAZ944on8wz/8Q2677baMGjWqt09h\nm9Dc3JwDDzywv8vYpr2az/30t/b29gwZMqS/y6AfmQOvDf35GUB/B8rN+GMOlFtzc3N23nnnrc5E\nFVkBHDFiRObPn7/J9vHjx2f8+PEbbTvnnHM2+v7nP/95JUoCAAAoPY+BAAAAKAkBEAAAoCQEQAAA\ngJIQAAEAAEpCAAQAACgJARAAAKAkBEAAAICSEAABAABKQgAEAAAoiQH9XQAAbMsmLJ7wivvePv32\nXqwEAHpmBRAAAKAkBEAAAICSEAABAABKQgAEAAAoCQEQAACgJARAAACAkvAYCLZ5r+YW7AAAUCZW\nAAEAAEpCAAQAACgJARAAAKAkBEAAAICSEAABAABKwl1AAaCfvJq7GN8+/fZerASAsrACCAAAUBIC\nIAAAQEkIgAAAACUhAAIAAJSEAAgAAFASAiAAAEBJCIAAAAAl4TmAALANmrB4Qtrb2zOkZcgr6u85\nggDlZAUQAACgJARAAACAkhAAAQAASsJnAHlNmLB4Qn+XAAAAr3tWAAEAAErCCiAAlNCrufLCHUQB\ntl1WAAEAAEpCAAQAACgJARAAAKAkfAYQANgqPj8IsO2yAggAAFASAiAAAEBJuASUXuFB7gBsCZeP\nAvQvK4AAAAAlYQUQANgmWD0EePUEQADgde/VflRBgAReLyoSAFtbW9PU1JShQ4dmzJgxmTp1apJk\n2bJlWbp0aYqiyJQpUzJ27NiceeaZGTx4cDo6OtLY2LhJm3HjxlWiRF6Cz/EBwEuz+gi8XlQkAC5Z\nsiTTpk3LuHHjcvzxx+fII4/MwIEDs2jRonznO99JV1dXTj755EyfPj2jR4/OF77whXzrW9/Kfffd\nt0mb7373uy95jA0bNiRJVqxYUYlT2CasWrUqTzzxxEbbpvxoSj9VQ1/rerYr6zes7+8y6EfmAObA\ntuH9F7+/Ivt99tlnM/g/Bldk30lyzRHXVGzf9I6Xei9IeaxatSpFUST5f9loS1QkALa1taW+vj5J\nUltbm/b29gwfPjxFUaS6ujpJ0tnZmba2ttTV1SVJ6urqsnLlyk3abM6qVauSpHt1Ecro2Tzb3yXQ\nz8wBzIFyq+T4/8Ol/1CxfQO9a9WqVdl99923qG1FAmB9fX1WrFiR+vr6rFmzJrW1tUmSmpqadHR0\npKurK9XV1amvr09zc3OSv63k7bXXXpu02ZyGhoZcffXVGTFiRLbffvtKnAYAAMBr1oYNG7Jq1ao0\nNDRscZ+q4oV1w160atWqNDU1ZfDgwWloaMjDDz+cWbNm5f77788NN9yQ9evX59hjj01DQ0O+/vWv\ndwe9OXPm5N57792ozX777dfb5QEAAJRSRQIgAAAArz0eBA8AAFASAiAAAEBJCIAAAAAlUZG7gNL7\nHn/88XzpS1/K0qVLc9lll+XPf/5z2tvbM3v27HR2dqapqSlDhw7NmDFjPBrjdejF4z979uxs2LAh\nVVVVmTJlSmpqarJw4cIMGjQoEyZMyIc+9KH+LpdedN9992XJkiUZPHhw3vSmN+W5555LR0dHnnnm\nmZxzzjn5v//3/xr/17m/nwMtLS0ZNmxYkmTGjBl54oknsnTp0hRFkSlTpmTcuHH9XDG97bHHHsvF\nF1+cHXfcMQ0NDVm9erX3ASXy9+P/q1/9yvuAkvryl7+cD37wg3nyySdf1WuAALgNWLVqVa6//vq8\n4Q1vyLp16/KrX/0q3/3ud3P33Xfnuuuuy7p16zJt2rSMGzcuxx9/fI488sgMHDiwv8uml7x4/JPk\n4Ycfzn777Zftttsue+21V84+++zMmjUrO++8cz772c964X+d+etf/5qzzjorO+ywQz772c9m1KhR\nOeecc/KjH/0oP/vZz/KLX/zC+L/OvXgOfOYzn8mqVasycuTI7LDDDtl1110zd+7cfOc730lXV1dO\nPvnkfPe73+3vkull7e3tOfXUU1NXV5fPf/7zqa6u9j6gRF48/jNnzsyqVau8DyihRYsWZfDgwUny\nqrOAALgNGDFiRL7yla/kuOOOy9NPP53hw4cnSerq6rJy5cp0dnamvr4+SVJbW5v29vbuNmz7Xjz+\nRVHktNNOy7ve9a7cfvvtufLKK7N69ersvPPOSZKqqqp+rpbedsghh6QoiixYsCDjxo3rHuO6uro8\n9NBDxr8EXjwHPvGJT2SPPfbI/vvvn6uuuiq33npriqLofpxSZ2dnP1dLJey3335pbW3N5z//+Ywf\nPz5//OMfk3gfUBYvHv93vvOdGTdunPcBJXPbbbdlyJAhGTt2bLq6ul51FvAZwG3Mm970pjz99NNJ\nkhUrVmTkyJGpr6/PihUrkiRr1qxJbW1tf5ZIBT377LN55JFHkiTDhg1LZ2dndt5556xcuTJJ4qku\nrz/PPPNMzjjjjIwdOzZHHHFEWltbk/y///6N/+vfi+fAO9/5zjzxxBNJ/t9rQE1NTTo6OrJ27dru\nIMjry4MPPpjq6ur84Ac/yAMPPJCnnnoqifcBZfHi8W9pacl9992XxPuAMvnJT36S3/72t7npppty\n3XXXZfXq1Ule+WuA5wBuQ4477rhcfvnlWbx4cR599NH89a9/zdlnn521a9emqakpgwcPTkNDQ448\n8sj+LpUKeGH8582bl+Rvl4R89atfzdNPP50FCxZk4MCBOeywwzJhwoR+rpTedPrpp+fxxx/PLrvs\nku233z4777xznnvuuTz77LM555xz8vjjjxv/17kXz4GiKDJgwIAMGzYszzzzTObMmZPf/e53ueGG\nG7J+/foce+yx2W+//fq7ZHrZb3/721x22WUZOXJkBg4cmLq6Ou8DSuTvx/+FlX7vA8rnxhtvTE1N\nTdra2l7Va4AACAAAUBIuAQUAACgJARAAAKAkBEAAAICSEAABAABKQgAEAAAoCQEQAACgJARAAEpt\nyZIlOeuss/q7DADoEwIgAKX2+9//PnvvvfdW97vxxhtzyCGHZNGiRVvU/vvf/37e9773Zd26dUmS\nM844IwcddFCWL1++1ccGgFdKAASg1B5++OG89a1vfUV9Dz/88Bx77LFJkmnTpnWHuaeeeiqHH374\nRm1vvvnmfOQjH8ktt9ySJDn33HOz7777vorKAWDrCYAAlNof/vCHV7QC+Pcef/zxvOUtb0myaai8\n5557sttuu+XTn/50rr766ld9LAB4pQb0dwEA0F+efPLJvOENb0htbW2OOOKI7L///nnmmWcyfvz4\nfOpTn0pHR0fOP//8bLfddmlvb88pp5ySnXfeeZP9/PnPf87OO++c7bb727+rPvzwwxuFyuuvvz6T\nJ0/OHnvskerq6vzmN7/J/vvv32fnCQAvEAABKK0XPv/35JNPZv/99+++Gcy0adPyiU98Ivfee292\n3HHHnHTSSS+7n4ceemijwPfAAw/kIx/5SJJkzZo1ueuuu7J69epceeWVeeaZZ3LVVVcJgAD0C5eA\nAlBaL1yq+cADD+Ttb3979/Y3vvGN2W677XLAAQeks7MzX/3qV/PjH/94s/t58MEHu2/u8thjj+W2\n227rvgT05ptvzhFHHJEf/OAHufzyy3Pdddfll7/8ZVavXl3ZkwOAlyAAAlBaL1yq2dLS0h0AH3ro\nodTX16eqqiqDBw/OKaeckvPPPz8//OEPX3Y/XV1d+djHPpbvfOc72WuvvXLTTTcl+dvlnx//+Me7\n277hDW/IxIkTc91111X25ADgJbgEFIDSmj9/fpLk+OOPz9NPP53q6upsv/32+cpXvpIkaWxszIAB\nA/Lss8/m85///Gb38/DDD+fGG2/MDjvssMnPbr755k22NTY2Jkl++ctf9sJZAMCWEwABKL3vf//7\nL7n9haC2OT/96U8zaNCgVFVVvWT4ezlnnHFGHnzwwa3qAwCvVlVRFEV/FwEAAEDl+QwgAABASQiA\nAAAAJSEAAgAAlIQACAAAUBICIAAAQEkIgAAAACUhAAIAAJSEAAgAAFAS/z9abtsadoxVGwAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1075abc88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "MC_histogram(1e6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 1/4\n",
      "Iter 2/4\n",
      "Iter 3/4\n",
      "Iter 4/4\n"
     ]
    }
   ],
   "source": [
    "num_samples = [1e3, 1e4, 1e5, 1e6]\n",
    "num_repetitions = 25\n",
    "\n",
    "total_probas = []\n",
    "for i, num_sample in enumerate(num_samples):\n",
    "    print(\"Iter {}/{}\".format(i+1, len(num_samples)))\n",
    "    probas = monte_carlo_proba(num_repetitions, num_sample, verbose=False)\n",
    "    total_probas.append(probas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4AAAAHwCAYAAAD3iQG/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3XucVXW9//H3gAwqMgiGzCSYF0xN\nNAJCzfSUqZlZmUVKxDGzQ+WPfqWnRPFGQIEe8fLQY0oqaKljpJKKVkdNJe+iplhmB/BCOsMgioDK\nDDP794fH+UUIox1mD7Cez8fDR7PW3sz67K3f6uVae+2KUqlUCgAAAJu8Th09AAAAAOUhAAEAAApC\nAAIAABSEAAQAACgIAQgAAFAQAhAAAKAgBCAA/2sPPvhg9t1334wcOTJf+9rXcvTRR+fWW29Nkvz5\nz3/ORRdd9L/6/dddd12amppW23fDDTdk1113zeOPP966r6mpKXvvvXcuvPDC93yMhx9+OE8//fS7\nfv7SpUszduzY1td7wgknZNmyZe/6z++3337v6nl//96+/df//b//d63Pf/HFF3PnnXcmSX784x/n\nxRdffNcz/aNXX301N998c5Jk6tSpeeKJJ/7p3wXAhmGzjh4AgE3DPvvsk/POOy9JsmLFiowcOTI7\n7rhjdt999+y+++7/q9996aWX5ogjjlhj/0477ZRZs2Zl4MCBSZLZs2ene/fu/9Qxrr/++hx22GHZ\nbbfd3tXzTzzxxBx99NE5+OCDkyTTp0/PGWec0foerE9//9625YEHHsj8+fNz4IEH5tRTT/1fHfcv\nf/lL7rzzznzuc5/LqFGj/le/C4ANgwAEYL3r1q1bjjrqqPzmN7/Ja6+9ltra2px33nn55Cc/mZ12\n2ik777xzjj322Jx++ulZuXJlunbtmgkTJqSmpiYXX3xxbr/99jQ3N2f48OHp3LlzGhoacsIJJ+Ti\niy9e7TgHHHBA/vCHP6SlpSWdOnXKrFmz8tnPfrb18SuuuCKzZs3KZpttliFDhuSHP/xhLrzwwixc\nuDAvv/xyXnzxxZxyyinp2bNnZs+enaeeeir9+/fPH//4x0yfPj2dOnXK4MGD84Mf/GC14/7tb3/L\n4sWLW+MvSUaOHJkvfelLSZJf/OIX+d3vfpc33ngjPXv2zEUXXZRbbrkl119/fVpaWlY7g/enP/0p\nEyZMSOfOnVvfh/e///3v6n2++uqrM3PmzHTq1Cl77rlnTjnllEydOjVvvvlmPvKRj2T69OkZN25c\nbr311jz33HN55ZVX8uqrr2bEiBH53e9+lwULFuSss87KwIEDM2XKlMydOzevvvpqdtttt0yaNCmX\nXHJJnn766Vx33XV57LHHcthhh2XffffNKaeckoULF6a5uTnHHntsDjvssIwcOTK77bZb/vrXv2b5\n8uW54IILst12273nf3YAaF8uAQWgXWyzzTZ55ZVXVtv30ksv5ZxzzsnYsWNz1llnZeTIkfn5z3+e\n4447Luecc07+9Kc/5Z577smMGTMyY8aMPPvss/nyl7+c3r17v+MZsC5dumTgwIF56KGHsnz58ixf\nvjzV1dVJ3jp7ddttt6W2tja1tbV57rnn8vvf/z5JUllZmcsuuyynnnpqpk+fngEDBmT//ffPD3/4\nw2y55Za58MILM3369Fx77bWpr6/Pvffeu9pxFy1alL59+662r3PnzunevXtaWlry6quvZvr06Zkx\nY0aam5vz5JNPJkmqqqpy7bXXZt999239c6eddlrOOOOM/OIXv8jw4cMzefLkNV7nAw88sNoloJdd\ndlmSty6DPf3003Pddddlp512SqlUyqhRo3L44YfnU5/61Gq/Y/PNN8/ll1+eT3/607n77rtzySWX\nZNSoUZk1a1aWL1+eqqqqTJs2Lddff30ef/zx1NfX59vf/nb22WefHHXUUa2/57rrrkuvXr1SW1ub\nadOm5fzzz8+SJUuSJHvttVemT5+e/fbbL7NmzVrHPx0AdBRnAAFoFy+++GJrjL2tZ8+e6dmzZ5Lk\nmWeeyaWXXprLLrsspVIpm222WRYsWJC99tornTt3TufOnXPyySe3eZzDDz88s2bNyksvvZSDDz64\n9bOC8+fPz4c//OF06dIlSTJkyJD89a9/TZLWS1Krq6vT2Ni42u97/vnns2TJktZLHlesWJHnn38+\nt956a55//vn07Nkzp556aurq6lb7c01NTbntttvy+c9/Pl26dMmJJ56YLbfcMnV1dVm1alWSZMcd\nd1xj/kWLFrXO89GPfjRTpkxZ4zlruwR00qRJueKKK3L22Wdn4MCBKZVKa32fPvShDyVJunfvnv79\n+ydJevTo0XoGdsmSJa0zv/7662t85vJt8+bNy8c+9rEkyVZbbZWdd945L7zwwmrHqK6uzuLFi9c6\nCwAdRwACsN4tX748M2bMyAUXXJCGhobW/Z06/f8LT3baaad84xvfyKBBgzJv3rw8/PDD2WmnnXLt\ntdempaUlzc3NGTVqVC699NJUVFSkpaXlHY+199575yc/+UkWLVqUKVOmtN60ZKeddsq0adOyatWq\ndO7cOQ8//HCOOOKIPP3006moqFjj91RUVKRUKqVv376pqanJFVdckS5duuSGG27I7rvvnuHDh6/2\n/J49e+b222/PQQcdlCS56qqr8sQTT+SDH/xgbr/99syYMSNvvPFGjjzyyNYw+/vX/7Ztt902Tz/9\ndHbbbbc8/PDD2WGHHd71+/zLX/4yP/rRj9K1a9ccd9xxeeyxx9KpU6d3fK/e6TW/7Z577slLL73U\nejbvv/7rv1Iqld7xd+2888555JFHcvDBB2f58uV55pln1jgbCsCGSwACsF68fZlip06d0tzcnO9+\n97vZaaedVgvAvzdmzJiMGzcuK1euzJtvvplTTz01u+++e/bff/8MHz48LS0tGT58eCorKzNkyJCM\nGjUqV1111Roh06lTp+y333556aWXstVWW7Xu33XXXfOZz3ym9XcNHjw4Bx100Frv9PnhD38455xz\nTs4///x8/etfz8iRI9Pc3Jztttsun/nMZ9Z4/tlnn53x48fniiuuSFNTU7bffvtMnDgxm222WbbY\nYoscffTRSZLevXtn0aJFa33fJk6cmAkTJqRUKqVz5875yU9+stb39u/97Gc/y6677pqvfvWr6dat\nW/r06ZMPf/jD2WqrrfLTn/40e+yxx1qP+Y/22muvXHzxxRkxYkQqKirSr1+/LFq0KNtvv32eeeaZ\nTJ8+vfW5X/nKV3L66adn+PDhWblyZUaPHp1tttnmXR8LgI5VUVrX9SIAAABsMtwEBgAAoCAEIAAA\nQEEIQAAAgIIQgAAAAAWxSd0F9M0338zcuXPTu3fvdO7cuaPHAQAAKKvm5uY0NDRkwIAB2Xzzzdd4\nfJMKwLlz52bEiBEdPQYAAECHuvrqqzNkyJA19m9SAdi7d+8kb73Y6urqDp5mdXPnzs2AAQM6egyg\njKx7KBZrHopnQ1z3dXV1GTFiRGsb/aNNKgDfvuyzuro6ffv27eBpVldfX7/BzQS0L+seisWah+LZ\nkNf92j4S5yYwAAAABSEAAQAACkIAAgAAFIQABAAAKAgBCAAAUBACEAAAoCAEIAAAQEEIQAAAgIIQ\ngAAAAAUhAAEAAApCAAIAABSEAAQAACgIAQgAAFAQAhAAAKAgNivHQerr6zN58uT06NEju+yyS0aM\nGJEkue+++zJz5syUSqUMHz48gwYNyrJlyzJ27Ngcd9xxGThw4Ds+BwAAgPeuLAFYW1ubkSNHZtCg\nQfm3f/u3fOUrX0mXLl0ybdq0/Od//mdaWlry/e9/P5dcckl++tOfZvPNN2/9s+/0HAAAAN67sgTg\n4sWLU1NTkySpqqrKsmXL0qtXr5RKpVRWViZJmpqakiQnnXRSLrzwwtY/+07PacvcuXNTX1+/Pl/C\nejFnzpyOHgEoM+seisWah+LZ0NZ9Q0PDOh8vSwDW1NSkrq4uNTU1Wbp0aaqqqpIkXbt2TWNjY1pa\nWloj7x+9m+f8owEDBqRv377rbf71Yc6cORk8eHBHjwGUkXUPxWLNQ/FsiOt+4cKF63y8LAE4bNiw\nTJ48OTfeeGMOOeSQTJo0KWPGjMkxxxyT0047LatWrcrxxx//jn/23TwHAACAtpUlAHv37p0pU6as\nsX/o0KEZOnToGvu/+93vtvkcAAAA3htfAwEAAFAQAhAAAKAgBCAAAEBBCEAAAICCEIAAAAAFIQAB\nAAAKQgACAAAUhAAEAAAoCAEIAABQEAIQAACgIAQgAABAQQhAAACAghCAAAAABSEAAQAACkIAAgAA\nFIQABAAAKAgBCAAAUBACEAAAoCAEIAAAQEEIQAAAgIIQgAAAAAUhAAEAAApCAAIAABSEAAQAACgI\nAQgAAFAQAhAAAKAgBCAAAEBBCEAAAICCEIAAAAAFIQABAAAKQgACAAAUhAAEAAAoCAEIAABQEAIQ\nAACgIAQgAABAQQhAAACAghCAAAAABSEAAQAACkIAAgAAFIQABAAAKAgBCAAAUBACEAAAoCAEIAAA\nQEEIQAAAgIIQgAAAAAUhAAEAAApCAAIAABSEAAQAACgIAQgAAFAQAhAAAKAgBCAAAEBBbFaOg9TX\n12fy5Mnp0aNHdtlll4wYMSJJct9992XmzJkplUoZPnx4Bg4cmNNPPz3dunVLY2Njxo0bl9mzZ+fu\nu+9OY2NjPvKRj+SLX/xiOUYGAADY5JTlDGBtbW1GjhyZcePG5a677kpTU1OSZNq0aZk4cWImTJiQ\nqVOn5sEHH0y/fv0yduzY9OrVK48++mgeeeSR/OlPf8qzzz6b6urqcowLAACwSSrLGcDFixenpqYm\nSVJVVZVly5alV69eKZVKqaysTJI0NTVl8eLFrZFXXV2dRYsW5WMf+1i+853v5M0338zJJ5+cfffd\nt83jzZ07N/X19e33gv5Jc+bM6egRgDKz7qFYrHkong1t3Tc0NKzz8bIEYE1NTerq6lJTU5OlS5em\nqqoqSdK1a9c0NjampaUllZWVqampaX0D6+rq0r9//5x33nmZNm1aunfvnubm5nd1vAEDBqRv377t\n9nr+GXPmzMngwYM7egygjKx7KBZrHopnQ1z3CxcuXOfjZQnAYcOGZfLkybnxxhtzyCGHZNKkSRkz\nZkyOOeaYnHbaaVm1alWOP/74DBgwIDfddFMmTpyYJBk0aFCOPPLInHjiidlyyy3z9a9/vRzjAgAA\nbJLKEoC9e/fOlClT1tg/dOjQDB06dLV948ePX237i1/8ohu/AAAArAe+BgIAAKAgBCAAAEBBCEAA\nAICCEIAAAAAFIQABAAAKQgACAAAUhAAEAAAoCAEIAABQEAIQAACgIAQgAABAQQhAAACAghCAAAAA\nBSEAAQAACkIAAgAAFIQABAAAKAgBCAAAUBACEAAAoCAEIAAAQEEIQAAAgIIQgAAAAAUhAAEAAApC\nAAIAABSEAAQAACgIAQgAAFAQAhAAAKAgBCAAAEBBCEAAAICCEIAAAAAFIQABAAAKQgACAAAUhAAE\nAAAoCAEIAABQEAIQAACgIAQgAABAQQhAAACAghCAAAAABSEAAQAACkIAAgAAFIQABAAAKAgBCAAA\nUBACEAAAoCAEIAAAQEEIQAAAgIIQgAAAAAUhAAEAAApCAAIAABSEAAQAACgIAQgAAFAQAhAAAKAg\nBCAAAEBBbFaOg9TX12fy5Mnp0aNHdtlll4wYMSJJct9992XmzJkplUoZPnx4Bg4cmNNPPz3dunVL\nY2Njxo0bl3vuuSd33HFHmpub84lPfCIHHXRQOUZer56vfy33P70svfu+lu37VHX0OAAAQEGV5Qxg\nbW1tRo4cmXHjxuWuu+5KU1NTkmTatGmZOHFiJkyYkKlTp+bBBx9Mv379Mnbs2PTq1SuPPvpoZsyY\nka233jqrVq3Khz70oXKMu149X/9afvSzB/LbR5fmRz97IM/Xv9bRIwEAAAVVljOAixcvTk1NTZKk\nqqoqy5YtS69evVIqlVJZWZkkaWpqyuLFi1NdXZ0kqa6uzqJFi/LUU0/l7LPPzuLFi3P++efn7LPP\nbvN4c+fOTX19ffu9oPfg/qeXZdErbyRJFr3yRm66/dHsu1v3Dp4KKJc5c+Z09AhAGVnzUDwb2rpv\naGhY5+NlCcCamprU1dWlpqYmS5cuTVXVW5dBdu3aNY2NjWlpaUllZWVqampa38C6urr0798/2223\nXbp27Zqtt976XR9vwIAB6du3b7u8lveqd9/X8tiCB7LolTeybc8t8vmDBrkMFApizpw5GTx4cEeP\nAZSJNQ/FsyGu+4ULF67z8bIE4LBhwzJ58uTceOONOeSQQzJp0qSMGTMmxxxzTE477bSsWrUqxx9/\nfAYMGJCbbropEydOTJIMGjQoI0aMyJgxY7Jq1ap85zvfKce469X2fapy5r/tk5tuf1T8AQAAHaos\nAdi7d+9MmTJljf1Dhw7N0KFDV9s3fvz41bYPPfTQHHrooe06X3vbvk9V9t2tu/gDAAA6lK+BAAAA\nKAgBCAAAUBACEAAAoCAEIAAAQEEIQAAAgIIQgAAAAAUhAAEAAApCAAIAABSEAAQAACgIAQgAAFAQ\nAhAAAKAgBCAAAEBBCEAAAICCEIAAAAAFIQABAAAKQgACAAAUhAAEAAAoCAEIAABQEAIQAACgIAQg\nAABAQQhAAACAghCAAAAABSEAAQAACkIAAgAAFIQABAAAKAgBCAAAUBACEAAAoCAEIAAAQEEIQAAA\ngIIQgAAAAAUhAAEAAApCAAIAABSEAAQAACgIAQgAAFAQAhAAAKAgBCAAAEBBCEAAAICCEIAAAAAF\nsdm7edL999+f559/Ph/+8Iez4447pmvXru09FwAAAOtZmwF47rnnpq6uLvPmzUtlZWWmTp2ac889\ntxyzAQAAsB61eQnonDlzcvbZZ2fLLbfMF7/4xSxcuLAccwEAALCetRmAzc3NWblyZSoqKtLc3JxO\nnXxsEAAAYGPU5iWgxxxzTI488sgsWbIkw4YNy7HHHluOuQAAAFjP2gzAj3zkI7nmmmvy3HPPpW/f\nvnn11VfLMRcAAADr2Vqv53zmmWcye/bsfOtb38qTTz6Z1157LXPnzs0JJ5xQzvkAAABYT9Z6BvC1\n117LrbfempdffjmzZs1KklRUVOSrX/1q2YYDAABg/VlrAA4ZMiRDhgzJU089lT322KOcMwEAANAO\n2vwMYF1dXc4999w0NTWlVCrl1Vdfzc0331yO2QAAAFiP2vxOh/PPPz+jR49OTU1NvvjFL2bXXXct\nx1wAAACsZ20G4LbbbpuPfOQjSZIjjzwy9fX17T4UAAAA61+bAdilS5c8/PDDWbVqVWbPnp1XXnml\nHHMBAACwnrX5GcAf/ehHmT9/fr7zne/kggsuyPHHH/+eD1JfX5/JkyenR48e2WWXXTJixIgkyX33\n3ZeZM2emVCpl+PDhGThwYE4//fR069YtjY2NGTduXJJk2bJlOeqoo3LllVemd+/e7/n4HW3ylQ/k\n/ifrs+/cppx8zD4dPQ4AAFBQbZ4BfN/73peePXvm9ddfzze/+c307dv3PR+ktrY2I0eOzLhx43LX\nXXelqakpSTJt2rRMnDgxEyZMyNSpU/Pggw+mX79+GTt2bHr16pVHH300LS0tOffcc7P99tu/91e3\nAZh85QO594n6tJSSe5+oz+QrH+jokQAAgIJq8wzgqFGj0tjYmKqqqiRvfRfgRRdd9J4Osnjx4tTU\n1CRJqqqqsmzZsvTq1SulUimVlZVJkqampixevDjV1dVJkurq6ixatCgXXXRRjjrqqEyfPv1dH2/u\n3LkbzGcV73+yfo3tOXPmdNA0QLlZ71As1jwUz4a27hsaGtb5eJsBuHLlyvziF7/4Xw1RU1OTurq6\n1NTUZOnSpa0x2bVr1zQ2NqalpSWVlZWpqalpfQPr6uqyww475I9//GMWL16cxx57LJdddllOOeWU\nNo83YMCAf+pMZXvYd25T7n3i/0fgvnv2yeDBgztwIqBc5syZY71DgVjzUDwb4rpfuHDhOh9vMwCH\nDBmS2bNnZ+edd27d9/73v/89DTFs2LBMnjw5N954Yw455JBMmjQpY8aMyTHHHJPTTjstq1atyvHH\nH58BAwbkpptuysSJE5MkQ4cOzdChQ5MkJ598cr75zW++p+NuCE4+Zp///xnAPfv4DCAAANBh2gzA\nl19+OT/5yU9WuwS0trb2PR2kd+/emTJlyhr7/z7w3jZ+/Ph3/B2TJ09+T8fckJx8zD4b5L8dAAAA\niqXNAJw/f35uu+22cswCAABAO2rzLqC77rprHn/88TQ2Nrb+BQAAwManzTOADz/8cO66667W7YqK\nitxxxx3tORMAAADtoM0AvPnmm8sxBwAAAO1srQE4fvz4nHHGGTnqqKNSUVGx2mPv9SYwAAAAdLy1\nBuDxxx+fJDnrrLPSpUuX1v1Lly5t/6kAAABY79Z6E5hSqZQFCxbkpJNOSlNTUxobG/Pmm2/mjDPO\nKOd8AAAArCdrPQP4xz/+MVdeeWUWLFiQM844I6VSKZ06dcrHP/7xcs4HAADAerLWADzooINy0EEH\n5e67786//Mu/lHMmAAAA2kGb3wPYpUuX3HPPPbn77rtz0EEHuSsoAADARqrNADzvvPOyww475Kqr\nrsq1117rDqAAAAAbqTYDcPPNN88222yTzTbbLL17917jKyEAAADYOLQZgFtttVW++c1v5jOf+Uyu\nvvrq9OrVqxxzAQAAsJ6t9SYwb7vgggvy/PPPp3///nnmmWcybNiwcswFAADAerbWM4DXX399kqSy\nsjKlUilJ8sEPfjBTp04tz2QAAACsV2sNwF//+tetP0+YMKH154ceeqh9JwIAAKBdrDUA3z7rt66f\nAQAA2HisNQD//m6fa/sZAACAjcdabwLz6quv5t57701LS0uWLl2aP/zhDymVSlm6dGk55wMAAGA9\nWWsA7rHHHrnllluSJB/60Icya9as1p8BAADY+Kw1ACdNmlTOOQAAAGhnbX4RPAAAAJsGAVgGP7/t\nqfzH9S/m57c91dGjAAAABdZmAI4fPz5//vOfyzHLJunntz2VX97+31mxsiW/vP2/RSAAANBh2gzA\nT3ziE7nkkkty9NFH55prrsny5cvLMdcm43cPPL/ObQAAgHJpMwAPOOCAXHDBBbn44oszZ86c7L//\n/jn55JPz/PNC5t04ZJ/t17kNAABQLmu9C+jb5s2blxtuuCG///3vs/fee+fqq6/OqlWr8v3vfz83\n3HBDOWbcqI38zB5Jklmz5+ez++/Uug0AAFBubQbgaaedlmHDhmX06NHZYostWvd/6UtfatfBNiUj\nP7NHPrTtmxk8WPwBAAAdp81LQPfff/8ceeSRrfE3ZcqUJMmIESPadzIAAADWq7WeAZwxY0Z+9atf\nZd68ebnnnnuSJM3NzVm1alX+/d//vWwDAgAAsH6sNQC/8IUvZN99982ll16ab3/720mSTp06ZZtt\ntinbcAAAAKw/aw3Av/zlL9lzzz1zyCGHZMGCBa37582bl49//ONlGQ4AAID1Z60BeP/992fPPffM\nrbfeusZjAhAAAGDjs9YA/PrXv57Gxsb86Ec/Kuc8AAAAtJO1BuChhx6aioqK1faVSqVUVFTkjjvu\naPfBADZmtU/8Or+Zf1cO7bIwR+/1hY4eBwAgyToC8M477yznHACbjNonfp0b/vybJGn9TxEIAGwI\n1hqA48ePzxlnnJGjjjpqjTOBtbW17T4YwMbqjnn3rrEtAAGADcFaA/D4449Pkpx77rllGwZgU/CB\nHn3zRMOfV9sGANgQrDUA3/e+9yVJWlpacvbZZ+fZZ5/NLrvskh/+8IdlGw5gY7TkzVfXuQ1sela8\nsDCr7n8wK7btk279/EsfYMPVqa0njB07Nl/+8pdzzTXX5PDDD8/YsWPLMRfARuuQ/gescxvYtKx4\nYWGeHHtaVv3XHXly7GlZ8cLCjh4JYK3aDMDOnTvnX/7lX9K9e/cceOCBaWlpKcdcAButAX12Tbcu\nWyZJunXZMgP67NrBEwHtaeHMX6f5tWVJkubXlmXhr2/q4IkA1m6tl4D+4Q9/SJJsscUW+dnPfpaP\nfvSjeeKJJ1ovDQXgnT1e91RWNL2eJFnR9Hoer/tT+vao6eCpgPbyxnMv/MP28x00CUDb1hqAs2bN\nSpJsvfXWmT9/fubPn58kqaysLM9kABupgdV7ZGbX3+W1lctS1bV7BlZ/qKNHAtpRl216JX/9h22A\nDdRaA3DSpEnvuH/RokXtNgzApqBueUNeW/nW5WCvrVyWuuUNzgDCJqxTxT9sp+KdnwhsUm6e8Wj+\n+PCLeXF+RT43bFBHj/OutfkZwAsuuCD77LNPBg8enD322CPHHntsOeYC2GjdOf/edW4Dm5amZctX\n316+rIMmAcrl5hmP5rEH/paW5uSxB/6Wm2c82tEjvWttBuCdd96Ze+65J5/73Ody6623pk+fPuWY\nC2CjdeBO+6Xif84AVKQiB+60XwdPBLSn5lWrVt9uWrWWZwKbisce+Ns6tzdkbQZg7969U1lZmRUr\nVuQDH/hAmpqayjEXwEareqve2ep/7gK6VZctU71V7w6eCGhPrz/9l3VuA2xI2gzA6urq/OpXv8oW\nW2yRKVOm5LXXXivHXAAbrcfrnsqyphVJkmVNK/J43Z86eCKgXXXqtO5tgA1Im/8NNX78+HzsYx/L\nSSedlG233TZTpkwpx1wAG62B1Xuke5duSZLuXbq5Cyhs4nof/Kl1bgObns26rHt7Q7bWu4C+benS\npbnqqqvy7LPPZpdddvEZQIA21C1vyPL/+R7A5U2vuwsobOKa/+HqqH/cBjY9fXfYJs/+9eXVtjcW\nbQbgmDFj8olPfCJHHHFEHnnkkYwZMyYXX3zxezpIfX19Jk+enB49emSXXXbJiBEjkiT33XdfZs6c\nmVKplOHDh2fgwIE5/fTT061btzQ2NmbcuHG55ppr8swzz+T111/PoYcemgMPPPCfe6UAZXLn/HtT\nSilJUkopd86/N0O226uDpwLay/IFz66+/eyz7/g8YNPx/r49VgvA9/ft0YHTvDdtXgK6cuXKfPWr\nX81uu+2Wr33ta1m27L3f2rhj00kyAAAY8ElEQVS2tjYjR47MuHHjctddd7XeSGbatGmZOHFiJkyY\nkKlTp+bBBx9Mv379Mnbs2PTq1SuPPvpoqqqqMm7cuJxyyim55ZZb3vsr3ABcNOPR/Pi6hbloI7o9\nLPDPcxdQKJbtjvj86ttf+PxanglsKj780X7pXtU1SdK9qms+/NF+HTzRu7fWAFywYEEWLFiQnj17\n5rbbbktDQ0PuuOOO9O3b9z0fZPHixampeevyp6qqqtaILJVKqayszOabb56mpqYsXrw41dXVSd66\n+cyiRYty+OGHZ8WKFTnrrLMyatSof+Y1dqiLZjya3z7wQpqak98+8IIIhIL4+zOAwKbt/Z85ND33\n2zfp1Ck999s37//MoR09EtDOevepyte+vU92H9Q9X/v2Pundp6qjR3rX1noJ6BlnnNH68zXXXJNr\nr702pVIpFRUV7/kgNTU1qaurS01NTZYuXZqqqrfeoK5du6axsTEtLS2prKxMTU1N5syZkySpq6tL\n//79M2/evPz0pz/N9773vfTr9+7Keu7cuamvr3/Pc7aHOx9e+A/bL2TfnfwfQtiU3fDi71bbvvHR\nWamo8xU6sKlqfOiRtNx7f5LklXvvz/0/vTSVQ4d08FRAOey0W/c8v/CveX5h288tl4aGhnU+XlEq\nldqskVdeeSUvvPBC+vbtm169ev1TQ0yePDndunXLgAED8pe//CVjxozJ448/nl/96ldZtWpVjj32\n2AwYMCBnnnlmKisrkySnnnpqDjnkkOy2227p2rVrPvjBD67zLODChQvzqU996p8+U9ke3j4D+LZP\n79Mvo4cN6sCJgPb2yN+eyDl/uDQtaUmndMoPPv4tnwGETdij3/1+3nj+//9v/Rbb98ugC8/vwImA\ncpkzZ04GDx7c0WOspq0mavMmMLfddlvOP//87LzzzvnrX/+a0aNH5wtf+MJ7GqJ3797v+PURQ4cO\nzdChQ1fbN378+NW2/+u//us9HWtD83bs3fnwCznwo+IPimDIdnvlBx//Vm58dFa+OOiz4g82cb32\n2Tt/+7sA7LXP3h04DcC6tRmA06dPzw033JBu3bpl+fLlOeaYY95zABbd6GGDsu9OpQweLP6gKP77\n5QX525uL8t8vLxCAsInrstVWq293795BkwC0rc27gFZUVKRbt7e+0HirrbZK165d230ogI1Z7RO/\nzg1//k1eb3kzN/z5N6l94tcdPRLQjrYe9JF03bZ3kqTrtr2z9UcGdvBEAGvX5hnAfv36ZfLkyRky\nZEgeeeSRbL/99uWYC2Cjdce8e9fYPnovV07Apqpbv77Z/YzTMvemm7P75z+Xbv02jPsQALyTNs8A\n/vjHP06/fv1y3333pV+/fpkwYUI55gLYaH1q5/3WuQ1serr165vN9hkq/oANXptnAL/97W/niiuu\nKMcsAJuEj39gaG575vd5o3lltujcNR//wNC2/xAAQBm0GYBVVVW54447ssMOO6RTp7dOGO64447t\nPhjAxurxuqfyRvPKJMkbzSvzeN2f0rdHTQdPBQDwLgLw5ZdfzvTp01u3KyoqctVVV7XnTAAbtYHV\ne+S2LX+fhteXpPeWvTKw+kMdPRIAQJI2AnD58uWZOnVqtthii3LNA7DR69ujJqccMDq3PPzbHP7R\nTzv7BwBsMNZ6E5hf/OIX+fznP58vfOELmT17djlnAtjo9e1Rk4/23FP8AQAblLUG4C233JLf/OY3\nqa2tzZVXXlnOmQA2eguXvpSHX3kyC5e+1NGjAAC0WmsAVlZWprKyMr169UpTU1M5ZwLYqC1c+lIm\n3XNR7nz5wUy65yIRCABsMNr8HsAkKZVK7T0HwCbj8bqn0vD6kiRJw+tL8njdnzp4IgCAt6z1JjD/\n/d//nX//939PqVRq/fltU6ZMKctwABujgdV75MYuv8myphXp3qWbu4ACABuMtQbg+eef3/rz0Ucf\nXZZhADYFdcsbsrzp9STJ8qbXU7e8wc1gAIANwloDcOjQoeWcA2CTcef8e1PKW5fOl1LKnfPvzZDt\n9urgqQAA3uVnAAF497bv8f51bgMAdBQBCLCePb/0xXVuAwB0FAEIsJ4duNN+6fQ///XaKZ1y4E77\ndfBEAABvEYAA69mQ7fbKDz7+reyy5fb5wce/5fN/AMAGY603gQHgnzdku71SUdeUweIPANiAOAMI\nAABQEAIQAACgIAQgAABAQfgMYBmcfukf8vgzL2fgI29kwrc+3tHjAAAABeUMYDt7O/6S5PFnXs7p\nl/6hgycCAACKSgC2syf+J/7Wtg0AAFAuArCd7fXBbda5DQAAUC4CsJ1N+NbHM/B/om/gB7fxGUAA\nAKDDuAlMGUz41sczZ86cDB48uKNHAQAACswZQAAAgIIQgAAAAAUhAAEAAApCAAIAABSEAAQAACgI\nAQgAAFAQAhAAAKAgBCAAAEBBCEAAAICCEIAAAAAFIQABAAAKQgACAAAUhAAEAAAoCAEIAABQEAIQ\nAACgIAQgAABAQQhAAACAghCAAAAABSEAAQAACkIAAgAAFIQABAAAKAgBCAAAUBCbleMg9fX1mTx5\ncnr06JFddtklI0aMSJLcd999mTlzZkqlUoYPH56BAwfm9NNPT7du3dLY2Jhx48at8ZxBgwaVY+T1\n6qGnXsov716c5s1fytA9ajp6HAAAoKDKcgawtrY2I0eOzLhx43LXXXelqakpSTJt2rRMnDgxEyZM\nyNSpU/Pggw+mX79+GTt2bHr16pVHH310jedsbB566qX8eNpD+cvf3syPpz2Uh556qaNHAgAACqos\nZwAXL16cmpq3znxVVVVl2bJl6dWrV0qlUiorK5MkTU1NWbx4caqrq5Mk1dXVWbRo0RrPeTfmzp2b\n+vr6dngl790v716cltJbP7eUkl/+9ol0fvPFjh0KKJs5c+Z09AhAGVnzUDwb2rpvaGhY5+NlCcCa\nmprU1dWlpqYmS5cuTVVVVZKka9euaWxsTEtLSyorK1NTU9P6BtbV1aV///5rPOfdGDBgQPr27dtu\nr+e9aN78rTOALaWkU0XylU/vlcEuA4VCmDNnTgYPHtzRYwBlYs1D8WyI637hwoXrfLwsAThs2LBM\nnjw5N954Yw455JBMmjQpY8aMyTHHHJPTTjstq1atyvHHH58BAwbkpptuysSJE5MkgwYNyqpVq1Z7\nzsZm6B41OfXYofnlb5/IVz69l88AAgAAHaYsAdi7d+9MmTJljf1Dhw7N0KFDV9s3fvz4Np+zsRm6\nR006v/miM38AAECH8jUQAAAABSEAAQAACkIAAgAAFIQABAAAKAgBCAAAUBACEAAAoCAEIAAAQEEI\nQAAAgIIQgAAAAAUhAAEAAApCAAIAABSEAAQAACgIAQgAAFAQAhAAAKAgBCAAAEBBCEAAAICCEIAA\nAAAFIQABAAAKQgACAAAUhAAEAAAoCAEIAABQEAIQAACgIAQgAABAQQhAAACAghCAAAAABSEAAQAA\nCkIAAgAAFIQABAAAKAgBCAAAUBACEAAAoCAEIAAAQEEIQAAAgIIQgAAAAAUhAAEAAApCAAIAABSE\nAAQAACgIAQgAAFAQAhAAAKAgBCAAAEBBCEAAAICCEIAAAAAFIQABAAAKQgACAAAUhAAEAAAoCAEI\nAABQEAIQAACgIAQgAABAQQhAAACAghCAAAAABSEAAQAACkIAAgAAFMRm5TjIzJkz88gjj+SNN97I\n6NGjs+OOO7Y+dvbZZ6exsTHLly/P+PHjM3/+/EydOjWbb755PvnJT+aAAw7IuHHjsvXWW+ell15q\n/RkAAID3pmwBOH369LzwwguZOnVqJkyYkCR54YUXWsPv+uuvz29/+9vMnj07Y8aMSZ8+ffKNb3wj\ne+65Z770pS9lyJAhueKKK/Lkk09m//33L8fYAAAAm5R2CcDrrrsut9xyS+t2165dkyTV1dVZtGhR\n6/6Ghob06dOn9bGnn346S5Ysad1XUVGR6urqVFdX54knnsjTTz+df/3Xf23z+HPnzk19ff36fEnr\nxZw5czp6BKDMrHsoFmseimdDW/cNDQ3rfLxdAvCoo47KUUcd1bo9atSoJEldXV223Xbb1v01NTWt\nofb2Y3369MmiRYuy7bbbplQqJUlqa2tTV1eXn/zkJ9lss7ZHHjBgQPr27bs+X9L/2pw5czJ48OCO\nHgMoI+seisWah+LZENf9woUL1/l4WS4BPeKII3Lqqadm+fLlOemkk1JfX58ZM2Zk9OjR2XrrrTNx\n4sSsWLEi48ePz+67756zzjorXbp0yciRI/PQQw/lP//zP7P33nvn5JNPzle/+tUMGjSoHGMDAABs\nUsoSgIcddlgOO+yw1faNHj06SXLiiSeutr9///6ZMmXKavtmz57dvgMCAAAUgK+BAAAAKAgBCAAA\nUBACEAAAoCAEIAAAQEEIQAAAgIIQgAAAAAUhAAEAAApCAAIAABSEAAQAACgIAQgAAFAQAhAAAKAg\nBCAAAEBBCEAAAICCEIAAAAAFIQABAAAKQgACAAAUhAAEAAAoCAEIAABQEAIQAACgIAQgAABAQQhA\nAACAghCAAAAABSEAAQAACkIAAgAAFIQABAAAKAgBCAAAUBACEAAAoCAEIAAAQEEIQAAAgIIQgAAA\nAAUhAAEAAApCAAIAABSEAAQAACgIAQgAAFAQAhAAAKAgBCAAAEBBCEAAAICCEIAAAAAFIQABAAAK\nQgACAAAUhAAEAAAoCAEIAABQEAIQAACgIAQgAABAQQhAAACAghCAAAAABSEAAQAACkIAAgAAFIQA\nBAAAKAgBCAAAUBCbleMgM2fOzCOPPJI33ngjo0ePzo477tj62Nlnn53GxsYsX74848ePz/z58zN1\n6tRsvvnm+eQnP5mDDz44SXLPPffk6quvzqWXXlqOkQEAADY5ZQvA6dOn54UXXsjUqVMzYcKEJMkL\nL7zQGn7XX399fvvb32b27NkZM2ZM+vTpk2984xs5+OCDM2/evDz22GNZtWpVOcYFAADYJLVLAF53\n3XW55ZZbWre7du2aJKmurs6iRYta9zc0NKRPnz6tjz399NNZsmRJ676KioosXbo006dPz5lnnplv\nfetb7+r4c+fOTX19/fp6OevNnDlzOnoEoMyseygWax6KZ0Nb9w0NDet8vF0C8KijjspRRx3Vuj1q\n1KgkSV1dXbbddtvW/TU1Na2h9vZjffr0yaJFi7LtttumVCrl9ttvz+uvv57x48dn3rx5+f3vf59P\nfvKT6zz+gAED0rdv33Z4Zf+8OXPmZPDgwR09BlBG1j0UizUPxbMhrvuFCxeu8/GyXAJ6xBFH5NRT\nT83y5ctz0kknpb6+PjNmzMjo0aOz9dZbZ+LEiVmxYkXGjx+f3XffPWeddVa6dOmSkSNH5pOf/GS+\n9KUvJUmOO+64NuMPAACAd1aWADzssMNy2GGHrbZv9OjRSZITTzxxtf39+/fPlClT3vH3XH755e0z\nIAAAQAH4GggAAICCEIAAAAAFIQABAAAKQgACAAAUhAAEAAAoCAEIAABQEAIQAACgIAQgAABAQQhA\nAACAghCAAAAABSEAAQAACkIAAgAAFIQABAAAKAgBCAAAUBACEAAAoCA26+gB1qfm5uYkSV1dXQdP\nsqaGhoYsXLiwo8cAysi6h2Kx5qF4NsR1/3YLvd1G/2iTCsCGhoYkyYgRIzp4EgAAgI7T0NCQD3zg\nA2vsryiVSqUOmKddvPnmm5k7d2569+6dzp07d/Q4AAAAZdXc3JyGhoYMGDAgm2+++RqPb1IBCAAA\nwNq5CQwAAEBBCEAAAICCEIAAAAAFIQABAAAKQgACAAAUxCb1PYAbi9tvvz2LFi1Kknz1q1/t4GmA\ncvnlL3+ZioqKDBs2rKNHAdrZjBkzUl9fn/79++fQQw/t6HGAMpg9e3aefvrpVFZW5phjjunocdbK\nGcAOcNBBB6VHjx5ZuXJlR48ClMkjjzzyjt/FA2ya5s2bl6222iqbbebftUNRPPTQQ+natWv69evX\n0aOskwDsAPPmzctnP/vZvPLKKx09ClAms2fPzoIFC/Loo4929ChAGRxzzDH5+te/nkceeaSjRwHK\nZMWKFfnXf/3XPPTQQx09yjr511Id4Mknn8wdd9yR6urqjh4FKJMTTjghCxcuzP3339/RowBlcPvt\nt6e5uTkf+MAHOnoUoEz222+/XHXVVampqenoUdapolQqlTp6iE3Nc889l+9973uZOXNm6uvrM3ny\n5PTo0SO77LJLRowY0dHjAe3AuodiseaheDaVde8S0PWsoaEhM2bMyBZbbJEkqa2tzciRIzNu3Ljc\nddddaWpq6uAJgfXNuodiseaheDaldS8A17PevXvnBz/4QbbccsskyeLFi1tPA1dVVWXZsmUdOR7Q\nDqx7KBZrHopnU1r3ArCd1dTUpK6uLkmydOnSVFVVdfBEQHuz7qFYrHkono153bsJTDsbNmxYJk+e\nnBtvvDGHHHKI20FDAVj3UCzWPBTPxrzu3QQGAACgIFwCCgAAUBACEAAAoCAEIAAAQEEIQAAAgIIQ\ngAAAAAUhAAEAAApCAAIAABSEAARgg/bggw9m8ODBeemll1r3nXPOObnhhhv+6d+5cOHCfOUrX1kf\n461h1apVGTlyZI4++ugsXbq0df+SJUvy3e9+N9/4xjdy9NFH59RTT82bb7653o9/8skn55577lnv\nvxeATYMABGCDV1lZmVNOOSWlUqmjR2nTokWLsmLFitTW1qZHjx6t+y+77LJ87GMfyxVXXJHa2tps\nueWWqa2t7cBJASiizTp6AABoyz777JOWlpZcffXV+drXvta6f+HChTnxxBPzy1/+Mknyla98Jeee\ne25uvPHGPPfcc3nllVfy6quvZsSIEfnd736XBQsW5Kyzzsr73ve+LFmyJN/+9rfz8ssv5xOf+ET+\nz//5P3nppZdy+umnZ+XKlenatWsmTJiQ5ubmfOc738nWW2+dAw44IP/2b//WevybbropV155ZSor\nK7PDDjtk/PjxOfPMM/Pss8/mjDPOyPjx41uf+773vS+//e1v84EPfCCDBg3KmDFjUlFRkSSZMmVK\n5s6dm1dffTW77bZbJk2alAsvvLDN1/C9730vvXv3Tn19fQ444ICccMIJrcdramrKmWeemeeeey4t\nLS35/ve/n7333jvnnXdeHnzwwaxatSqHHHJIRo0a1d5/+wDYgDgDCMBGYdy4cZk+fXqee+65d/X8\nzTffPJdffnk+/elP5+67784ll1ySUaNGZdasWUmS119/Pf/xH/+R2trazJ49O08//XTOOuusjBw5\nMj//+c9z3HHH5ZxzzkmSNDQ05PLLL18t/l555ZVceOGFufLKK3Pttdeme/fuue6663LmmWemf//+\nq8Vfknz961/P4Ycfnssvvzz7779/Ro8enUWLFmX58uWpqqrKtGnTcv311+fxxx9PfX39u3oNf/vb\n3zJ58uT86le/ygMPPJCnnnqq9XgzZsxIz549c/XVV+fiiy9unefmm2/OOeeck2uuuSZVVVX/5N8N\nADZWzgACsFHo2bNnxo4dmzFjxmTQoEHv+Jy/v0T0Qx/6UJKke/fu6d+/f5KkR48eWblyZZJkt912\nS/fu3ZMke+65ZxYsWJBnnnkml156aS677LKUSqVsttlb/zPZt2/fVFZWrnasF154If37989WW22V\nJPnoRz+aP/zhD/nEJz7xjrM98MADOeKII/LlL385jY2N+dnPfpaf/OQnOffcc7NkyZKceOKJ2XLL\nLfP666+nqanpXb+GrbfeOkmy1157ZcGCBa3He+aZZzJnzpw88cQTSd76bOKSJUvyH//xH5kyZUoW\nL16c/ffff91vOgCbHGcAAdhoHHjggdlxxx1z4403Jkm6du2al19+Oc3NzXnttdf+Xzv3z5JaHMdx\n/O2BxCGQApdQI0noOfgEDrqGSInPIKIQGmzwPAHHlpbQB+DS0NwUje5uRlDRlP+guHe4XOH+8eLq\nPe/X+vvH7yyHD9/f78doNFr0/Xm8cpnhcMh4PObz85PBYECxWKRQKNBsNun1ekRRRBiGAATBn7/L\nbDbLcDhkMpkA8Pj4yN7e3tL1ut0ut7e3wI87jcVikWQyyf39Pc/Pz3Q6Hc7Pz5nNZosgu8oeptMp\nX19fDAaDRUgEKBQKVCoVer0e19fXhGHI5uYmd3d3dDodut0u/X6fp6enf64hSfq/WAGUJK2VVqvF\nw8MDAJlMhlKpxOHhIblcjt3d3ZXnSafTnJ2d8f7+TrlcZn9/n4uLC9rtNvP5nNlsRqvVWjp+e3ub\nk5MTGo0GQRCQz+dpNpu8vr7+tX8URURRxM3NDalUiq2tLdrtNkEQcHV1xfHxMYlEglwux8vLy0p7\n2NjY4PT0lLe3N8Iw5ODgYNFWq9W4vLykXq/z8fHB0dERyWSSdDpNtVollUpRKpXY2dlZ+ZtJktZf\n4ts6PKkmSZJ+8fsDOJIkrcIjoJIkSZIUE1YAJUmSJCkmrABKkiRJUkwYACVJkiQpJgyAkiRJkhQT\nBkBJkiRJigkDoCRJkiTFhAFQkiRJkmLiO9dhzShQ/M3LAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x107a06c88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "y_axis = np.asarray(total_probas)\n",
    "x_axis = np.asarray(num_samples)\n",
    "\n",
    "for x, y in zip(x_axis, y_axis):\n",
    "    plt.scatter([x] * len(y), y, s=12)\n",
    "\n",
    "plt.xscale('log')\n",
    "plt.title(\"Direct Monte-Carlo Estimation\")\n",
    "plt.ylabel(\"Probability Estimate\")\n",
    "plt.xlabel('Number of Samples')\n",
    "plt.grid(True)\n",
    "plt.savefig(plot_dir + 'monte_carlo_convergence_speed.pdf', format='pdf', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importance Sampling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "With importance sampling, we try to reduce the variance of our Monte-Carlo integral estimation by choosing a better distribution from which to simulate our random variables. It involves multiplying the integrand by 1 (usually dressed up in a “tricky fashion”) to yield an expectation of a quantity that varies less than the original integrand over the region of integration. Concretely,\n",
    "\n",
    "$$\\mathbb{E}_{p(x)} \\big[\\ f(x) \\big] = \\int f(x)\\ p(x)\\ dx = \\int f(x)\\ p(x)\\ \\frac{q(x)}{q(x)}\\ dx = \\int \\frac{p(x)}{q(x)}\\cdot f(x)\\ q(x)\\ dx = \\mathbb{E}_{q(x)}  \\big[\\ f(x)\\cdot \\frac{p(x)}{q(x)} \\big]$$\n",
    "\n",
    "Thus, the MC estimation of the expectation becomes:\n",
    "\n",
    "$$\\mathbb{E}_{q(x)}  \\big[\\ f(x)\\cdot \\frac{p(x)}{q(x)} \\big] \\approx \\frac{1}{N} \\sum_{n=1}^{N} w_n \\cdot f(x_n)$$\n",
    "\n",
    "where $w_n = \\dfrac{p(x_n)}{q(x_n)}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our current example above, we can alter the mean and/or standard deviation of $\\Delta L$ and $\\Delta V_{TH}$ in the hopes that more of our sampling points will fall in the failure region (red area). For example, let us define 2 new distributions with altered $\\sigma^2$.\n",
    "\n",
    "* $\\Delta \\hat{L} \\sim \\ N(0, 0.02^2)$\n",
    "* $\\Delta \\hat{V}_{TH} \\sim \\ N(0, 0.06^2)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def importance_sampling(num_simulations, num_samples, verbose=True, plot=False):\n",
    "    \n",
    "    if verbose:\n",
    "        print(\"===================================================\")\n",
    "        print(\"{} Importance Sampling Simulations of size {}\".format(num_simulations, num_samples))\n",
    "        print(\"===================================================\\n\")\n",
    "    \n",
    "    num_simulations = int(num_simulations)\n",
    "    num_samples = int(num_samples)\n",
    "    \n",
    "    probas = []\n",
    "    for i in range(num_simulations):\n",
    "        mu_1, sigma_1 = 0, 0.01\n",
    "        mu_2, sigma_2 = 0, 0.03\n",
    "        mu_1_n, sigma_1_n = 0, 0.02\n",
    "        mu_2_n, sigma_2_n = 0, 0.06\n",
    "        \n",
    "        # setup pdfs\n",
    "        old_pdf_1 = norm(mu_1, sigma_1)\n",
    "        new_pdf_1 = norm(mu_1_n, sigma_1_n)\n",
    "        old_pdf_2 = norm(mu_2, sigma_2)\n",
    "        new_pdf_2 = norm(mu_2_n, sigma_2_n)\n",
    "\n",
    "        length = np.random.normal(mu_1_n, sigma_1_n, num_samples)\n",
    "        voltage = np.random.normal(mu_2_n, sigma_2_n, num_samples)\n",
    "        \n",
    "        # calculate current\n",
    "        num = 50 * np.square((0.6 - voltage))\n",
    "        denum = 0.1 + length\n",
    "        I = num / denum\n",
    "        \n",
    "        # calculate f\n",
    "        true_condition = np.where(I >= 275)\n",
    "\n",
    "        # calculate weight\n",
    "        num = old_pdf_1.pdf(length) * old_pdf_2.pdf(voltage)\n",
    "        denum = new_pdf_1.pdf(length) * new_pdf_2.pdf(voltage)\n",
    "        weights = num / denum\n",
    "\n",
    "        # select weights for nonzero f\n",
    "        weights = weights[true_condition]\n",
    "\n",
    "        # compute unbiased proba\n",
    "        proba = np.sum(weights) / num_samples\n",
    "        probas.append(proba)\n",
    "        \n",
    "        false_condition = np.where(I < 275)\n",
    "        if plot:\n",
    "            if i == num_simulations -1:\n",
    "                plt.scatter(length[true_condition], voltage[true_condition], color='r')\n",
    "                plt.scatter(length[false_condition], voltage[false_condition], color='b')\n",
    "                plt.xlabel(r'$\\Delta L$ [$\\mu$m]')\n",
    "                plt.ylabel(r'$\\Delta V_{TH}$ [V]')\n",
    "                plt.title(\"Monte Carlo Estimation of P(I > 275)\")\n",
    "                plt.grid(True)\n",
    "                plt.savefig(plot_dir + 'imp_sampling_{}.pdf'.format(num_samples), format='pdf', dpi=300)\n",
    "                plt.show()\n",
    "    \n",
    "    \n",
    "    mean_proba = np.mean(probas)\n",
    "    std_proba = np.std(probas)\n",
    "    \n",
    "    if verbose:\n",
    "        print(\"Probability Mean: {}\".format(mean_proba))\n",
    "        print(\"Probability Std: {}\".format(std_proba))\n",
    "    \n",
    "    return probas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===================================================\n",
      "10 Importance Sampling Simulations of size 10000\n",
      "===================================================\n",
      "\n",
      "Probability Mean: 0.0021612338069065583\n",
      "Probability Std: 6.967937733495277e-05\n"
     ]
    }
   ],
   "source": [
    "probas = importance_sampling(10, 10000, plot=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def IS_histogram(num_samples, plot=True):\n",
    "    \n",
    "    num_samples = int(num_samples)\n",
    "\n",
    "    mu_1_n, sigma_1_n = 0, 0.02\n",
    "    mu_2_n, sigma_2_n = 0, 0.06\n",
    "\n",
    "    length = np.random.normal(mu_1_n, sigma_1_n, num_samples)\n",
    "    voltage = np.random.normal(mu_2_n, sigma_2_n, num_samples)\n",
    "\n",
    "    # calculate biased current\n",
    "    num = 50 * np.square((0.6 - voltage))\n",
    "    denum = 0.1 + length\n",
    "    I = num / denum\n",
    "\n",
    "    if plot:\n",
    "        n, bins, patches = plt.hist(I, 50, normed=1, facecolor='green', alpha=0.75)\n",
    "        plt.ylabel('Number of Samples')\n",
    "        plt.xlabel(r'$I_{DS}$ [$\\mu$A]')\n",
    "        plt.title(\"Importance Sampling of P(I > 275)\")\n",
    "        plt.grid(True)\n",
    "        plt.savefig(plot_dir + 'is_histogram_{}.pdf'.format(num_samples), format='pdf', dpi=300)\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3sAAAHtCAYAAACgbWeUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3X+clnWd7/H3IAwqMMAoMLdCpcHq\n2pgIhtr2Q8xFH2a16kElDotu2prHbdPsQRLoCCRY4TE17eAP5CEk/gjJ1D1t+QuTwhqyDVdNLS1c\nZ2BWwYGEGZg5f/RwToQwODAgF8/nX8x139d9fa6Zrze+uO6577LW1tbWAAAAUChddvUAAAAA7Hhi\nDwAAoIDEHgAAQAGJPQAAgAISewAAAAUk9gAAAApI7AGwieXLl+eMM87o9OP84he/yLPPPtvpx0mS\nlpaWTJ8+Peecc07Gjh2bc889N3/84x93+HEWLFiQb33rW1m5cmVqamp2+OMnycUXX5zTTz89L774\nYtu26667LieeeGLGjRuXcePG5ayzzsqSJUvabp82bVrq6ura5nunnnnmmXz2s5/NuHHj8rnPfS4N\nDQ155pln2o43bty4HH744Vm0aFFWrVqVo48+um37nDlz0tramgkTJmTdunU75HsAwLbpuqsHAGDP\n9P3vfz8nn3xyDj300E4/1uOPP54VK1Zk9uzZSZKf/OQnufLKK3PjjTd2yvH69evXabG3ePHi/Pzn\nP99s+9lnn50xY8YkSV588cVccskluffee/PUU0+la9euqaqqetvH+81vfpPVq1fn7/7u71JWVva2\n9/n617+eyZMn52//9m8zf/783HTTTbn00ktz++23J0n+7d/+Lf3798/HPvaxLF68OKecckomT568\nyWOccsopufnmm3PhhRduz+kD8A6IPQC2aNy4cTnkkEPy/PPPZ999981RRx2Vn/70p3njjTdy6623\n5qGHHspPfvKTrF27Nq+//nr+1//6XznxxBPzxBNP5Jprrkn37t3Tp0+fXHnllXnmmWfyrW99K926\ndcuHP/zhPP7443n66aczePDgPPzww/n3f//3vPnmm+nbt2+uv/763H///Xnssceybt26/OEPf8h5\n552X0047Lb/+9a9z5ZVXpqWlJQMGDMi3vvWtvPzyy5k2bVqStB2vV69ebefRt2/fLFu2LA8++GCO\nOeaYfOITn8jHPvaxJMn//b//N/PmzcuGDRtSVlaW66+/Ps8//3xmzZqVbt26pa6uLmeddVZ+/vOf\n59lnn80//uM/5rOf/WxOPvnkHHXUUXn++efTu3fvXH311W3HW758eS6++OLcdddd+dSnPpURI0bk\nueeeS1lZWW644Yb07NkzV1xxRZYtW5b9998/r7zySm688cYMHDiw7THe7nt49dVXZ82aNfnCF76w\n1VBdtWpV9t133yTJ7bffnnPOOWeL9+3fv38efPDBXHfddRk5cmROO+209O/ff5P7XH311W3bNm7c\nmO7du7fd9qc//SnXXXdd5s6dmyRZtmxZnn766fzP//k/U1lZmUmTJqV///758Ic/nBkzZuSCCy5I\nly5eWASwM3i2BWCrPvjBD2bOnDlpamrK3nvvndmzZ2fw4MH5xS9+kSR58803M3v27Nx6662ZMWNG\nmpubM3ny5Fx//fWZO3duPvShD7WFyfr16/O9730vF154YT760Y/mK1/5SqqqqrJq1arcdtttufvu\nu7Nx48b85je/SZKsWbMm/+f//J/ceOONmTVrVpLksssuy5VXXpm77747H//4x/Piiy9m8uTJufzy\ny3P77bfnYx/7WG6++ebNzmHq1Kn5yU9+klNOOSWnn356nnrqqSTJSy+9lFmzZuWOO+7I4MGD89Of\n/jRJUldXl+uuuy41NTW58cYb841vfCM33XRT7rzzziTJunXr8qlPfSp33HFHDj744Lbtf23t2rX5\n5Cc/mblz56Z///5ZtGhRHnrooaxatSr33HNPrrzyyrz66qub7NPa2vq238Oampr07t37bUPvtttu\ny7hx4zJ+/PjcdtttmTp1apLkySefzN/8zd9s8ec7YMCATJgwIfPmzUupVMqnPvWpPPzww5vc563Q\nW7p0aebOnZuzzz677bZ77rknJ510UiorK5MkBx98cL74xS9m7ty5OeGEE9oifK+99kplZWV++9vf\nbnEWAHYsV/YA2KoPfOADSZKKiooMHjy47c/r169PknzoQx9Kly5dsv/++6eioiINDQ3p2bNnBgwY\n0Hb71VdfneOOOy4HHXTQZo/fpUuXdOvWLRdffHH23Xff1NXVZcOGDUnS9hLPUqmUpqamJElDQ0Pe\n//73J0lGjx6d5M8vW7ziiiuSJM3NzXnf+963yTGeffbZHHTQQbn66qvT2tqaJ554Il/60pfyxBNP\nZL/99suECRPSo0eP/O53v8vQoUOTJEOGDEm3bt3Sq1evvOc970l5eXl69+7ddt5du3bNhz70oSTJ\nsGHDsmjRorZ9/9phhx3Wdh7r16/PK6+80nbfysrKHHzwwZvc//XXX3/b7+HW/OXLOP9SS0tLysvL\nt7rvsmXLctddd+X3v/99Lr744hxzzDGb3efBBx9si+63wi5JfvjDH+baa69t+/qYY47JPvvskyT5\n+7//+01u69+/f1atWrXVWQDYccQeANvl6aefTvLnCFuzZk369++fNWvWZMWKFenfv3+efPLJtvj6\ny5fvlZWVpbW1Nc8++2x+8pOf5O67786bb76Z0047La2trW33+Wv9+/fPSy+9lPe9732ZNWtWDjro\noBx00EG56qqrcsABB6S2tjYrV67cZJ+f/exneeGFFzJ16tR06dIlQ4YMyT777JM1a9bk2muvzaOP\nPpokOeecc7Z67L+0YcOGPPvsszn00ENTW1vbFsJv568fa8iQIfnBD36QJFm9enVeeumlTW7v27fv\nFr+H71T37t2zcePG7LXXXm97+8MPP5xHHnkkZ555Zqqrq9/2Pj/4wQ9y55135vbbb0+fPn3atjc2\nNqapqSmlUqlt26RJkzJq1KicfPLJ+dnPftb2jwVvnet+++3XofMA4J0TewBsl4aGhowfPz6NjY25\n/PLLs9dee2XatGn5l3/5l5SVlaV3796ZPn16nn/++U32O+KII/Ktb30rV199dfbZZ5+cddZZSf78\n5iYrVqzY4vGuuOKKTJw4MV26dEm/fv1y9tlnp1QqZcKECW2/d/f1r399k33GjRuXq666Kp/5zGfS\ns2fPdOnSJd/4xjfSs2fPDBs2LGeeeWa6du2aioqKrFixYpPfnduam266Kf/1X/+VAw44IBdddFHu\nv//+bdrvuOOOy6JFi3LWWWdl//33z957751u3bq13V5WVva238OOGDZsWJ5++ul88IMffNvbjz/+\n+Bx//PFb3H/jxo35+te/nlKplH/5l39J8ucrjV/84hfz+9//PgceeOAm9//yl7+ciRMn5o477sg+\n++zT9jLOlpaW1NfXbzWKAdixylrf+idMAHiHFixYkN/97ne55JJLdvUoO93xxx+ff/u3f9vkzUq2\n1Ysvvphnn302n/zkJ/P666/nlFNOySOPPNLuyy074le/+lUeeOCBTJo0aYc/9jvx2GOP5emnn84F\nF1ywS+cA2JN4gxYA2MlKpVLuv//+nHHGGTn33HNzySWXdEroJcmRRx6ZjRs3pq6urlMef1u0trbm\nhz/84SZv7AJA53NlDwAAoIBc2QMAACggsQcAAFBAu+27ca5bty7Lli1Lv379tvh20gAAAEW1cePG\nrFy5MtXV1dl77703u323jb1ly5Zl7Nixu3oMAACAXWrevHk56qijNtu+28Zev379kvz5xKqqqrJs\n2bItfhgsbI21Q0dZO3SUtUNHWTt0hHVTXHV1dRk7dmxbG/213Tb23nrpZlVVVQYOHJj6+vpt/hBc\n+EvWDh1l7dBR1g4dZe3QEdZN8W3p19q8QQsAAEABiT0AAIACEnsAAAAFJPYAAAAKSOwBAAAUkNgD\nAAAoILEHAABQQGIPAACggMQeAABAAYk9AACAAhJ7AAAABST2AAAACkjsAQAAFJDYAwAAKCCxBwAA\nUEBiDwAAoIDEHgAAQAGJPQAAgALq2hkPWl9fnxkzZqR3794ZMmRIxo4dmyRZvHhxFi5cmNbW1owZ\nMyZDhw7N5MmT06NHjzQ1NaWmpiaPP/54HnvssTQ1NeXII4/Mqaee2hkj8i4zcs7IDu/7yPhHduAk\nAABQDJ1yZW/+/PkZN25campq8uijj6a5uTlJMnv27EybNi1Tp07NrFmzsmTJkgwaNCgTJ05MZWVl\nli5dml/+8pf5z//8z7z00kupqqrqjPEAAAAKr1Ou7DU0NKRUKiVJKioq0tjYmMrKyrS2tqa8vDxJ\n0tzcnIaGhragq6qqyooVK/LhD384X/jCF7Ju3bp89atfzbHHHrvVYy1btiz19fVJktra2s44HXaC\nxsbGDu+7I37u1g4dZe3QUdYOHWXt0BHWTTGtXLlyq7d3SuyVSqXU1dWlVCpl9erVqaioSJJ07949\nTU1NaWlpSXl5eUqlUtvCq6ury+DBg/O///f/zuzZs9OrV69s3Lix3WNVV1dn4MCBqa2tzfDhwzvj\ndNgJei3r1eF9t/fnbu3QUdYOHWXt0FHWDh1h3RTX8uXLt3p7p8Te6NGjM2PGjNx7770ZNWpUpk+f\nngkTJmT8+PGZNGlSNmzYkAsuuCDV1dW57777Mm3atCTJsGHDctppp+Xiiy/Ovvvum7PPPrszxgMA\nACi8Tom9fv36ZebMmZttHzFiREaMGLHJtilTpmzy9amnnupNWQAAALaTj14AAAAoILEHAABQQGIP\nAACggMQeAABAAYk9AACAAhJ7AAAABST2AAAACkjsAQAAFJDYAwAAKCCxBwAAUEBiDwAAoIDEHgAA\nQAGJPQAAgAISewAAAAUk9gAAAApI7AEAABSQ2AMAACggsQcAAFBAYg8AAKCAxB4AAEABiT0AAIAC\nEnsAAAAFJPYAAAAKSOwBAAAUkNgDAAAoILEHAABQQGIPAACggMQeAABAAYk9AACAAuq6qweA7TVy\nzsgO7/vI+Ed24CQAAPDu4coeAABAAYk9AACAAhJ7AAAABST2AAAACkjsAQAAFJDYAwAAKCCxBwAA\nUEBiDwAAoIDEHgAAQAF13dUDUAwj54zc1SMAAAB/wZU9AACAAhJ7AAAABST2AAAACkjsAQAAFJDY\nAwAAKCCxBwAAUECd8tEL9fX1mTFjRnr37p0hQ4Zk7NixSZLFixdn4cKFaW1tzZgxYzJ06NBMnjw5\nPXr0SFNTU2pqajJ16tSsXbs29fX12X///fPNb36zM0YEAAAotE65sjd//vyMGzcuNTU1efTRR9Pc\n3JwkmT17dqZNm5apU6dm1qxZWbJkSQYNGpSJEyemsrIyS5cuzeTJkzNlypT07ds3l19+eWeMBwAA\nUHidcmWvoaEhpVIpSVJRUZHGxsZUVlamtbU15eXlSZLm5uY0NDSkqqoqSVJVVZUVK1YkSRYsWJBP\nfvKT6dmzZ7vHWrZsWerr65MktbW1nXE6bIPGxsZdPUKHvLVmrB06ytqho6wdOsraoSOsm2JauXLl\nVm/vlNgrlUqpq6tLqVTK6tWrU1FRkSTp3r17mpqa0tLSkvLy8pRKpbaFV1dXl8GDBydJHn300dx4\n443bdKzq6uoMHDgwtbW1GT58eGecDtug17Jeu3qEDhk+fLi1Q4dZO3SUtUNHWTt0hHVTXMuXL9/q\n7Z3yMs7Ro0dn7ty5ueyyyzJq1KhMnz49TU1NGT9+fCZNmpSJEyfmggsuyPDhw/PKK69k2rRpeeON\nNzJs2LCsXr065eXlKSsr64zRAAAA9gidcmWvX79+mTlz5mbbR4wYkREjRmyybcqUKZt83bt371x7\n7bWdMRYAAMAew0cvAAAAFJDYAwAAKCCxBwAAUEBiDwAAoIDEHgAAQAGJPQAAgAISewAAAAUk9gAA\nAApI7AEAABSQ2AMAACggsQcAAFBAYg8AAKCAxB4AAEABiT0AAIACEnsAAAAFJPYAAAAKSOwBAAAU\nkNgDAAAoILEHAABQQGIPAACggMQeAABAAYk9AACAAhJ7AAAABST2AAAACkjsAQAAFJDYAwAAKCCx\nBwAAUEBiDwAAoIDEHgAAQAGJPQAAgAISewAAAAUk9gAAAApI7AEAABSQ2AMAACggsQcAAFBAYg8A\nAKCAxB4AAEABiT0AAIACEnsAAAAFJPYAAAAKSOwBAAAUkNgDAAAoILEHAABQQGIPAACggMQeAABA\nAXXtjAetr6/PjBkz0rt37wwZMiRjx45NkixevDgLFy5Ma2trxowZk6FDh2by5Mnp0aNHmpqaUlNT\nk0WLFuWhhx7Kxo0bc9xxx+WEE07ojBEBAAAKrVOu7M2fPz/jxo1LTU1NHn300TQ3NydJZs+enWnT\npmXq1KmZNWtWlixZkkGDBmXixImprKzM0qVLc/fdd6dPnz7ZsGFDDjvssM4YDwAAoPA65cpeQ0ND\nSqVSkqSioiKNjY2prKxMa2trysvLkyTNzc1paGhIVVVVkqSqqiorVqzI008/nW984xtpaGjINddc\nk2984xtbPdayZctSX1+fJKmtre2M02EbNDY27uoROuStNWPt0FHWDh1l7dBR1g4dYd0U08qVK7d6\ne6fEXqlUSl1dXUqlUlavXp2KiookSffu3dPU1JSWlpaUl5enVCq1Lby6uroMHjw4Bx54YLp3754+\nffps07Gqq6szcODA1NbWZvjw4Z1xOmyDXst67eoROmT48OHWDh1m7dBR1g4dZe3QEdZNcS1fvnyr\nt3dK7I0ePTozZszIvffem1GjRmX69OmZMGFCxo8fn0mTJmXDhg254IILUl1dnfvuuy/Tpk1Lkgwb\nNixjx47NhAkTsmHDhnzhC1/ojPEAAAAKr1Nir1+/fpk5c+Zm20eMGJERI0Zssm3KlCmbfH3SSSfl\npJNO6oyxAAAA9hg+egEAAKCAxB4AAEABiT0AAIACEnsAAAAFJPYAAAAKSOwBAAAUkNgDAAAooE75\nnD3YXYycMzKNjY3ptazXO973kfGPdMJEAACwY7iyBwAAUEBiDwAAoIDEHgAAQAGJPQAAgAISewAA\nAAUk9gAAAApI7AEAABSQ2AMAACggsQcAAFBAYg8AAKCAxB4AAEABiT0AAIACEnsAAAAFJPYAAAAK\nSOwBAAAUkNgDAAAoILEHAABQQGIPAACggMQeAABAAYk9AACAAhJ7AAAABST2AAAACkjsAQAAFJDY\nAwAAKCCxBwAAUEBiDwAAoIDEHgAAQAGJPQAAgAISewAAAAUk9gAAAApI7AEAABSQ2AMAACggsQcA\nAFBAYg8AAKCAxB4AAEABiT0AAIACekex9+qrr3bWHAAAAOxAXdu7w80335yKioq88cYbWbBgQT76\n0Y/m0ksv3eo+9fX1mTFjRnr37p0hQ4Zk7NixSZLFixdn4cKFaW1tzZgxYzJ06NBMnjw5PXr0SFNT\nU2pqanLvvffm/vvvT79+/XL00Ufn1FNP3TFnCgAAsAdp98rev//7v+cf/uEfsmjRojz44IN55pln\n2n3Q+fPnZ9y4campqcmjjz6a5ubmJMns2bMzbdq0TJ06NbNmzcqSJUsyaNCgTJw4MZWVlVm6dGl+\n8YtfZMCAAWltbc0RRxyx/WcIAACwB2r3yl6XLl3S0NCQ/fffP0mybt26dh+0oaEhpVIpSVJRUZHG\nxsZUVlamtbU15eXlSZLm5uY0NDSkqqoqSVJVVZUVK1bk9NNPz+GHH541a9bka1/7Wm688cYOnxwA\nAMCeqt3YO/roozNu3Lh885vfzJVXXpmPf/zj7T5oqVRKXV1dSqVSVq9enYqKiiRJ9+7d09TUlJaW\nlpSXl6dUKqW2tjZJUldXl8GDB6e2tjZHHnlkevbsmdbW1naPtWzZstTX1ydJ22Ox8zU2Nu7qEbZL\nR+a33kisAzrO2qGjrB06wropppUrV2719rLWbSmqJKtWrcq+++7bdmWuvYPOmDEjPXr0SHV1dZ57\n7rlMmDAhTz31VO65555s2LAh55xzTqqrq3P55Ze3PeakSZPy/e9/P7/85S/T0tKSs846K0ceeeTb\nHmP58uX5xCc+kYceeigDBw5MbW1thg8fvi2nQicYOWfkrh6hwxobG9OrV693vN8j4x/phGnYnXje\noaOsHTrK2qEjrJvi+usm+mvtXtn7xS9+kSuuuCIbN27MSSedlAMOOCCjR4/e6j79+vXLzJkzN9s+\nYsSIjBgxYpNtU6ZM2eTr008/Paeffnp7YwEAALAV7b5ByzXXXJO5c+dm//33z/nnn5877rhjZ8wF\nAADAdmg39rp06ZI+ffqkrKws3bt3T48ePXbGXAAAAGyHdmPvPe95T2bOnJlVq1Zl1qxZOeCAA3bG\nXAAAAGyHdmPviiuuyAEHHJDhw4dn3333zdSpU3fGXAAAAGyHLb5By09/+tO2Pw8aNCiDBg1Kkjz5\n5JP5yEc+0vmTAQAA0GFbjL0HHnhgizuJPQAAgHe3Lcbe9OnT2/78n//5n/n973+fwYMH55BDDtkp\ngwEAANBx7X7O3jXXXJOf//zn+eAHP5jbb789J5xwQs4999ydMRsAAAAd1G7sLVq0KPfcc0+6dOmS\njRs35swzzxR7AAAA73LtvhtnVVVV1q5dmyTZsGFD9t9//04fCgAAgO3T7pW9FStW5MQTT8yhhx6a\nF154Id26dctZZ52VJJk/f36nDwgAAMA7127sffvb394ZcwAAALADtRt7//3f/50HHngg69evb9tW\nU1PTmTMBAACwndqNvQkTJuS8885LRUXFzpgHAACAHaDd2Hvve9+b0047bWfMAgAAwA7SbuydeOKJ\nueiii/L+97+/bduFF17YqUMBAACwfdqNvXnz5mXUqFFexgkAALAbaTf2+vTpk89//vM7YxZ2sZFz\nRu7qEQAAgB2k3djr27dvLrvsshx22GEpKytLkpx55pmdPhgAAAAdt01v0JIkDQ0NnT4MAAAAO0a7\nsXfhhRdmxYoV2bBhQ1pbW7NixYqdMRcAAADbod3YmzhxYp566qm8+eabWbduXQYNGpS77rprZ8wG\nAABAB3Vp7w7PPvtsHnjggXzkIx/JAw88kO7du++MuQAAANgO7cZe3759U1ZWlj/96U+prKzcGTMB\nAACwndqNvQ984AO55ZZb0r9//1x00UVZt27dzpgLAACA7dDu7+xdfPHFWbNmTfbee+8sWrQoH/zg\nB3fGXAAAAGyHLV7ZW79+febMmZPW1tasWbMmF198cR588MG0trbuzPkAAADogC3G3rRp0/Jf//Vf\naWlpyRVXXJFDDz00J554YmpqanbieAAAAHTEFl/G+fzzz2f+/PlZv359amtrc+2116Zbt2659dZb\nd+Z8AAAAdMAWr+z16NEjSbJ06dIcfvjh6datW5I/v7wTAACAd7ctXtnr0aNH7rzzzvzoRz/KKaec\nkpaWltx3330plUo7cz4AAAA6YItX9mpqavKHP/whH/3oR3PqqadmyZIl+dGPfuR39gAAAHYDW7yy\nV1lZma985SttXx977LE59thjd8pQAAAAbJ92P1QdAACA3c8WY6+xsXFnzgEAAMAOtMXY+/znP58k\nufzyy3faMAAAAOwYW/ydva5du+b000/Pyy+/nOeeey5J0tramrKyssyfP3+nDQgAAMA7t8XYu+22\n21JfX5+amprU1NSktbV1Z84FAADAdthi7O2111454IADcsMNN+TOO+/MCy+8kPe9730ZM2bMzpwP\nAACADmj33Tgvu+yy/OEPf8jf/d3f5ZVXXsmkSZN2xlwAAABshy1e2XvLyy+/nHnz5iVJTjjhhJx1\n1lmdPhQAAADbp90re+vXr8+bb76ZJFm3bl02btzY6UMBAACwfdq9sveP//iP+cxnPpMhQ4bkhRde\nyBe/+MWdMRcAAADbod3Y+/SnP52Pfexj+eMf/5iBAwemb9++O2MuAAAAtkO7sZckffr0SZ8+fTp7\nFtitjJwzssP7PjL+kR04CQAAbK7d39kDAABg99Pulb1bbrkln/vc597Rg9bX12fGjBnp3bt3hgwZ\nkrFjxyZJFi9enIULF6a1tTVjxozJ0KFDM3ny5PTo0SNNTU2pqalJkjQ2NubMM8/MnDlz0q9fv3d+\nVgAAAHu4dq/sPfbYY+/4HTjnz5+fcePGpaamJo8++miam5uTJLNnz860adMyderUzJo1K0uWLMmg\nQYMyceLEVFZWZunSpWlpacnVV1+d97znPR07IwAAANq/svf666/nox/9aAYOHJiysrKUlZVl/vz5\nW92noaEhpVIpSVJRUZHGxsZUVlamtbU15eXlSZLm5uY0NDSkqqoqSVJVVZUVK1bk+uuvz5lnnpnb\nbrttm05g2bJlqa+vT5LU1tZu0z68vcbGxl09wi6zs8/dWi0OP0s6ytqho6wdOsK6KaaVK1du9fZ2\nY++73/3uOz5oqVRKXV1dSqVSVq9enYqKiiRJ9+7d09TUlJaWlpSXl6dUKrUtvLq6urzvfe/Lr3/9\n6zQ0NORXv/pVbr755lx66aVbPVZ1dXUGDhyY2traDB8+/B3Pyv/Xa1mvXT3CLtHY2JhevXbuuVur\nxeB5h46ydugoa4eOsG6Ka/ny5Vu9vd3Y69q1a775zW/mtddey0knnZRDDjkkBx544Fb3GT16dGbM\nmJF77703o0aNyvTp0zNhwoSMHz8+kyZNyoYNG3LBBRekuro69913X6ZNm5YkGTFiREaMGJEk+epX\nv5pzzz13W88TAACAv9Bu7E2ePDnnnHNObrjhhhx11FH56le/mrvuumur+/Tr1y8zZ87cbPtfxtxb\npkyZ8raPMWPGjPZGAwAAYAvafYOWdevW5dhjj01ZWVkOPvjgdO/efWfMBQAAwHZoN/a6d++exx9/\nPC0tLXnqqafa3mAFAACAd692Y2/q1KlZsGBBXn/99dx6661tn4UHAADAu1e7v7NXVVWVf/7nf85L\nL72UIUOGZNCgQTtjLgAAALZDu7F3ww035PHHH8/hhx+e2267LSeddFLOPvvsnTAaAAAAHdVu7D32\n2GO544470qVLl2zYsCGf/eyvPPATAAAVH0lEQVRnxR4AAMC7XLu/s7fffvvlzTffTJI0NzensrKy\n04cCAABg+2zxyt6ZZ56ZsrKy/Pd//3dOPPHEHHLIIXnxxRfTp0+fnTkfAAAAHbDF2Lv66qt35hwA\nAADsQFuMvQMPPDBJ8h//8R954IEHsn79+rbbfPwCAADAu1u7b9AyYcKEnHfeeamoqNgZ8wAAALAD\ntBt7733ve3PaaaftjFkAAADYQdqNvRNPPDEXXXRR3v/+97dtu/DCCzt1KAAAALZPu7E3b968jBo1\nyss4AQAAdiPtxl6fPn3y+c9/fmfMAgAAwA7Sbuz17ds3l112WQ477LCUlZUl+fNn8AEAAPDutU1v\n0JIkDQ0NnT4MAAAAO0a7seedOAEAAHY/7cbeRRddlLKysrS0tGT58uV573vfmzvuuGNnzAYAAEAH\ntRt7d955Z9uf33jjjUyePLlTBwIAAGD7dXknd+7Vq1f++Mc/dtYsAAAA7CDtXtk788wzU1ZWltbW\n1rz22ms59thjd8ZcAAAAbId2Y+/qq69u+3P37t2z//77d+pAAAAAbL8txt7ChQu3uNM//MM/dMow\nAAAA7BhbjL0XX3xxk69bW1uzYMGC7L333mIPAADgXW6LsfflL3+57c9/+MMfMmHChBx33HGZOHHi\nThkMAACAjmv3d/bmzZuXOXPm5NJLL83IkSN3xkwAAABspy3GXn19fS699NL07t07d999d3r37r0z\n5wIAAGA7bDH2PvnJT6a8vDzHHHNMpkyZssltM2fO7PTBAAAA6Lgtxt4NN9ywM+cAAABgB9pi7I0Y\nMWJnzgEAAMAO1GVXDwAAAMCOJ/YAAAAKSOwBAAAUkNgDAAAoILEHAABQQGIPAACggMQeAABAAYk9\nAACAAhJ7AAAABST2AAAACkjsAQAAFJDYAwAAKCCxBwAAUEBiDwAAoIC6dsaD1tfXZ8aMGendu3eG\nDBmSsWPHJkkWL16chQsXprW1NWPGjMnQoUMzefLk9OjRI01NTampqcmPf/zjPPzww2lpacn48eNz\n2GGHdcaIAAAAhdYpV/bmz5+fcePGpaamJo8++miam5uTJLNnz860adMyderUzJo1K0uWLMmgQYMy\nceLEVFZWZunSpSkrK8uUKVMyduzY/PjHP+6M8QAAAAqvU2KvoaEhpVIpSVJRUZHGxsYkSWtra8rL\ny7P33nunubk5DQ0NqaqqSpJUVVVlxYoVOeGEE/KrX/0qkydPzjHHHNMZ4wEAABRep7yMs1Qqpa6u\nLqVSKatXr05FRUWSpHv37mlqakpLS0vKy8tTKpVSW1ubJKmrq8vgwYPzs5/9LMccc0wWLFiQc889\nN0cfffRWj7Vs2bLU19cnSdtj0TFvRfmeaGefu7VaHH6WdJS1Q0dZO3SEdVNMK1eu3OrtnRJ7o0eP\nzowZM3Lvvfdm1KhRmT59eiZMmJDx48dn0qRJ2bBhQy644IJUV1fnvvvuy7Rp05Ikw4YNy913350J\nEyZk7733zqc//el2j1VdXZ2BAwemtrY2w4cP74zT2WP0WtZrV4+wSzQ2NqZXr5177tZqMXjeoaOs\nHTrK2qEjrJviWr58+VZv75TY69evX2bOnLnZ9hEjRmTEiBGbbJsyZcomX48ePTqjR4/ujLEAAAD2\nGD56AQAAoIDEHgAAQAGJPQAAgAISewAAAAUk9gAAAApI7AEAABSQ2AMAACggsQcAAFBAYg8AAKCA\nuu7qAWBPNHLOyO3a/5Hxj+ygSQAAKCpX9gAAAApI7AEAABSQ2AMAACggsQcAAFBAYg8AAKCAxB4A\nAEABiT0AAIACEnsAAAAFJPYAAAAKSOwBAAAUkNgDAAAoILEHAABQQGIPAACggMQeAABAAYk9AACA\nAhJ7AAAABST2AAAACkjsAQAAFJDYAwAAKCCxBwAAUEBiDwAAoIDEHgAAQAGJPQAAgAISewAAAAUk\n9gAAAApI7AEAABSQ2AMAACggsQcAAFBAYg8AAKCAxB4AAEABiT0AAIACEnsAAAAFJPYAAAAKSOwB\nAAAUkNgDAAAoILEHAABQQF0740Hr6+szY8aM9O7dO0OGDMnYsWOTJIsXL87ChQvT2tqaMWPGZOjQ\noZk8eXJ69OiRpqam1NTU5Hvf+15++9vf5k9/+lNOOumkHH/88Z0xIgAAQKF1ypW9+fPnZ9y4camp\nqcmjjz6a5ubmJMns2bMzbdq0TJ06NbNmzcqSJUsyaNCgTJw4MZWVlVm6dGkqKipSU1OTSy+9NPff\nf39njAcAAFB4nXJlr6GhIaVSKUlSUVGRxsbGVFZWprW1NeXl5UmS5ubmNDQ0pKqqKklSVVWVFStW\n5JRTTsnatWtz1VVX5fOf/3y7x1q2bFnq6+uTJLW1tZ1xOnuMxsbGXT3CLrO7nbu1/u7hZ0FHWTt0\nlLVDR1g3xbRy5cqt3t4psVcqlVJXV5dSqZTVq1enoqIiSdK9e/c0NTWlpaUl5eXlKZVKbQuvrq4u\ngwcPzosvvpgbb7wx//qv/5pBgwa1e6zq6uoMHDgwtbW1GT58eGeczh6j17Jeu3qEXaKxsTG9eu1e\n526tvzt43qGjrB06ytqhI6yb4lq+fPlWb++Ul3GOHj06c+fOzWWXXZZRo0Zl+vTpaWpqyvjx4zNp\n0qRMnDgxF1xwQYYPH55XXnkl06ZNyxtvvJEjjzwy559/ftavX59vf/vbmTVrVmeMBwAAUHidcmWv\nX79+mTlz5mbbR4wYkREjRmyybcqUKZt8/eMf/7gzRgIAANij+OgFAACAAhJ7AAAABdQpL+Nk1xk5\nZ+SuHgEAAHgXcGUPAACggMQeAABAAYk9AACAAhJ7AAAABST2AAAACkjsAQAAFJDYAwAAKCCxBwAA\nUEA+VB12QyPnjOzwvo+Mf2QHTgIAwLuVK3sAAAAFJPYAAAAKSOwBAAAUkNgDAAAoILEHAABQQGIP\nAACggMQeAABAAYk9AACAAhJ7AAAABST2AAAACkjsAQAAFJDYAwAAKCCxBwAAUEBiDwAAoIDEHgAA\nQAGJPQAAgAISewAAAAUk9gAAAApI7AEAABSQ2AMAACggsQcAAFBAYg8AAKCAxB4AAEABiT0AAIAC\nEnsAAAAFJPYAAAAKSOwBAAAUkNgDAAAoILEHAABQQGIPAACggLru6gGAnWvknJEd3veR8Y/swEkA\nAOhMruwBAAAUkNgDAAAooE55GWd9fX1mzJiR3r17Z8iQIRk7dmySZPHixVm4cGFaW1szZsyYDB06\nNJMnT06PHj3S1NSUmpqaJMmiRYtyzz335Nprr+2M8QAAAAqvU2Jv/vz5GTduXIYNG5bzzjsvZ5xx\nRrp165bZs2fnO9/5TlpaWvKlL30p48ePz6BBg3L++efn2muvzdKlS9PU1JQ//vGPWbt2bWeMBgAA\nsEfolNhraGhIqVRKklRUVKSxsTGVlZVpbW1NeXl5kqS5uTkNDQ2pqqpKklRVVWXFihU56aSTcswx\nx+Thhx/epmMtW7Ys9fX1SZLa2tpOOJvdS2Nj464eYbfk+7Zt/De2Od8TOsraoaOsHTrCuimmlStX\nbvX2Tom9UqmUurq6lEqlrF69OhUVFUmS7t27p6mpKS0tLSkvL0+pVGpbeHV1dRk8ePA7PlZ1dXUG\nDhyY2traDB8+fIeex+6o17Jeu3qE3U5jY2N69fJ92xb+G9uU5x06ytqho6wdOsK6Ka7ly5dv9fZO\neYOW0aNHZ+7cubnssssyatSoTJ8+PU1NTRk/fnwmTZqUiRMn5oILLsjw4cPzyiuvZNq0aXnjjTcy\nbNiwzhgHAABgj9MpV/b69euXmTNnbrZ9xIgRGTFixCbbpkyZ8raPccstt3TGaAAAAHsEH70AAABQ\nQGIPAACggMQeAABAAYk9AACAAhJ7AAAABST2AAAACkjsAQAAFJDYAwAAKCCxBwAAUEBiDwAAoIDE\nHgAAQAF13dUDALuPkXNGdnjfR8Y/sgMnAQCgPa7sAQAAFJDYAwAAKCCxBwAAUEBiDwAAoIDEHgAA\nQAGJPQAAgAISewAAAAUk9gAAAApI7AEAABSQ2AMAACggsQcAAFBAYg8AAKCAuu7qAYA9w8g5I7dr\n/0fGP7KDJgEA2DO4sgcAAFBAYg8AAKCAvIzzXWh7X+4GAADgyh4AAEABiT0AAIACEnsAAAAFJPYA\nAAAKSOwBAAAUkHfjBHYL2/MutT6QHQDYE7myBwAAUEBiDwAAoIDEHgAAQAGJPQAAgAISewAAAAXk\n3TiBwvNOngDAnsiVPQAAgAISewAAAAXkZZwAW7G1l4A2Njam17JeW7zdS0ABgF3JlT0AAIACcmUP\noJN4YxgAYFcSewDvQtsTiolYBAA6Kfbq6+szY8aM9O7dO0OGDMnYsWOTJIsXL87ChQvT2tqaMWPG\nZOjQoZk8eXJ69OiRpqam1NTUbHafYcOGdcaInWp7/ycNYHvtquchkQkA7x6dEnvz58/PuHHjMmzY\nsJx33nk544wz0q1bt8yePTvf+c530tLSki996UsZP358Bg0alPPPPz/XXnttli5dutl9vvvd777t\nMTZu3JgkqaurS5KsXLkyy5cv74zTecc2rNqwq0fgHWhZ25ING/3MeOesnc199Nsf3dUj7HR3nH7H\nO97n3fR3FrsXa4eOsG6K660WequN/lqnxF5DQ0NKpVKSpKKiIo2NjamsrExra2vKy8uTJM3NzWlo\naEhVVVWSpKqqKitWrNjsPluycuXKJGm7agjbY23W7uoR2E1ZO3zihk/s6hEA2MOtXLky733vezfb\n3imxVyqVUldXl1KplNWrV6eioiJJ0r179zQ1NaWlpSXl5eUplUqpra1N8ucqHTx48Gb32ZLq6urM\nmzcv/fr1y1577dUZpwEAAPCutXHjxqxcuTLV1dVve3tZa2tr644+6MqVKzNjxoz06NEj1dXVee65\n5zJhwoQ89dRTueeee7Jhw4acc845qa6uzuWXX94WdZMmTcqTTz65yX0OP/zwHT0eAABA4XVK7AEA\nALBr+VB1AACAAhJ7AAAABST2AAAACqhT3o1zZ9rSB7jDX1q6dGnmz5+fHj16ZL/99suf/vSnNDU1\nZc2aNZkyZUp+97vfZdasWdl7770zcuTI/P3f//2uHpl3mS9/+cs5/vjj8+qrr+aVV15JY2NjJk6c\nmObmZs9BvK3ly5fnhhtuSM+ePdO7d++sXbvW8w7bpK6uLtdff3169uyZJNl///0979Cul19+Of/6\nr/+ahQsX5uabb253zfz1fSorK3f1KdAJdvvY29IHuMNfeuONN3LZZZelZ8+e+ad/+qcMHDgwU6ZM\nyfe///386Ec/yuOPP54JEyZkwIAB+ad/+if/08UmZs+enR49eiRJfvnLX+a73/1ufv7zn+euu+7K\n+vXrPQfxtmbPnp1Bgwbl5ZdfzpAhQ/Kb3/zG8w7b5He/+11+9rOfpbq6OoceeqjnHdq1cuXK3H33\n3dlnn32yfv36dtfMZz7zmc3uc/755+/q06AT7PYv43y7D3CHv3bcccelR48eufHGGzNs2LAMGDAg\nSVJVVZUVK1bktddea9tWVla2K0flXeahhx5Kr169MnTo0LS0tLT9y+dba8dzEFvy8ssv5+Mf/3im\nTp2a++67z/MO26yqqiq33XZbrrnmmjz55JPp27dv23bPO7ydfv365ZJLLsm+++6bVatWtft31Rtv\nvLHZfSim3T723voA9ySbfIA7/KU1a9bka1/7WoYOHZrTTz899fX1Sf78Upn+/ftnwIABbU90Po2E\nv/TDH/4w//Ef/5F77703d911V1577bUk/3/teA5iS/r165eePXu2XXHxvMO2mjdvXhobG1NWVpZe\nvXrllVdeSeJ5h22z3377ZdWqVUm2vGb69++/2X0opt3+c/b++gPczzjjjF09Eu9Cl156aV5++eUc\ncMAB2WuvvTJgwID86U9/ytq1azNlypS8/PLLufHGG9OtW7eceOKJGTly5K4emXeZBQsWpHv37mlo\naMjvf//7vPHGG7niiiuybt06z0G8rRdffDHXXnttevbsmaOPPjovvPCC5x22ydNPP53rrrsupVIp\nBx54YLp16+Z5h23yuc99LrfcckvmzJnT7pr56/v06tVrV49PJ9jtYw8AAIDN7fYv4wQAAGBzYg8A\nAKCAxB4AAEABiT0AAIACEnsAAAAFJPYAAAAKSOwBsEebP39+Lrvssl09BgDscGIPgD3ab3/72xxy\nyCHveL8FCxbkuOOOy+zZs7fp/jfddFM+8pGPZP369UmSr33taznqqKPy4osvvuNjA8C2EHsA7NGe\ne+65/M3f/E2H9j3llFNyzjnnJEnGjRvXFm6vv/56TjnllE3ue9999+Xkk0/OAw88kCT5+te/nr/9\n27/djskBYOvEHgB7tOeff75DV/b+2ssvv5yDDjooyeYBuWTJkrznPe/JWWedlXnz5m33sQBgW3Td\n1QMAwK7y6quvZp999klFRUVOP/30HHHEEVmzZk1GjBiR//E//keamppy1VVXpUuXLmlsbMxFF12U\nAQMGbPY4r7zySgYMGJAuXf78b6jPPffcJgF59913Z/To0Tn44INTXl6eX//61zniiCN22nkCsGcS\newDssd76fb1XX301RxxxRNsbtYwbNy6nnnpqnnzyyfTt2zcXXnjhVh/n2Wef3STunn766Zx88slJ\nktWrV2fRokV57bXXcvvtt2fNmjWZO3eu2AOg03kZJwB7rLdebvn000/nAx/4QNv2fffdN126dMmR\nRx6Z5ubmfOUrX8kPfvCDLT7OM8880/bGKy+99FIeeuihtpdx3nfffTn99NNz66235pZbbsldd92V\nJ554Iq+99lrnnhwAezyxB8Ae662XWy5btqwt9p599tmUSqWUlZWlR48eueiii3LVVVfle9/73lYf\np6WlJZ/+9Kfzne98J4MHD869996b5M8v4fzMZz7Tdt999tkno0aNyl133dW5JwfAHs/LOAHYY82c\nOTNJct5552XVqlUpLy/PXnvtlUsuuSRJUlNTk65du2bt2rX553/+5y0+znPPPZcFCxakZ8+em912\n3333bbatpqYmSfLEE0/sgLMAgLcn9gDY4910001vu/2tKNuS+++/P3vvvXfKysreNvS25mtf+1qe\neeaZd7QPALwTZa2tra27eggAAAB2LL+zBwAAUEBiDwAAoIDEHgAAQAGJPQAAgAISewAAAAUk9gAA\nAApI7AEAABSQ2AMAACig/wfr/7APpe/P/wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10b039588>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "IS_histogram(1e5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 1/4\n",
      "Iter 2/4\n",
      "Iter 3/4\n",
      "Iter 4/4\n"
     ]
    }
   ],
   "source": [
    "num_samples = [1e3, 1e4, 1e5, 1e6]\n",
    "num_repetitions = 25\n",
    "\n",
    "total_probas = []\n",
    "for i, num_sample in enumerate(num_samples):\n",
    "    print(\"Iter {}/{}\".format(i+1, len(num_samples)))\n",
    "    probas = importance_sampling(num_repetitions, num_sample, verbose=False)\n",
    "    total_probas.append(probas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4AAAAHwCAYAAAD3iQG/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xuc1nWd///nIAygOBxygHEBxaD4\nJpYhkmWaa/7M1E4WShKhfdM1s1vqutF64riKJqxtbhaWx1QKCzPLbMVV8KywpXgLNUQRcw4IclKY\nAeb3h+t8lwwGW+a6YD73++3Wrfl8rmvm87qG21t98Lmuz6eiubm5OQAAALR7Hco9AAAAAKUhAAEA\nAApCAAIAABSEAAQAACgIAQgAAFAQAhAAAKAgBCAAO51ly5blxBNPbPPjPP7441m0aFGbHydJNm/e\nnEsvvTSnnnpqRo8ena9+9at56aWXdvhxfvGLX+SKK65IQ0NDJkyYsMN/PgC7NgEIQGH9/Oc/T319\nfUmONW/evNTX1+e6667LzTffnFGjRuWSSy5ps+NVV1cLQADepmO5BwCAbRkzZkze+9735rnnnsvu\nu++e4cOH54EHHsjq1atz7bXXZs6cObnnnnuybt26rFy5Ml//+tfziU98Ig8++GCuvPLKdO7cOT16\n9Mgll1ySP/7xj7niiivSqVOnfOQjH8m8efPy9NNPZ9CgQbn33nvzu9/9Lm+88UZ69uyZq666Knfe\neWfuv//+rF+/PkuXLs1pp52WE044IX/4wx9yySWXZPPmzenTp0+uuOKKvPjii5kyZUqStBxvzz33\nbHkdPXv2zMKFC/Ob3/wmhxxySD7+8Y/n8MMPT5L89re/zc0335yNGzemoqIiV111VZ577rnMmDEj\nnTp1Sm1tbUaNGpVHHnkkixYtype//OWcfPLJOfbYYzN8+PA899xz6d69e6ZPn95yvGXLluXcc8/N\nz372s3zqU5/KiBEj8swzz6SioiLf//73061bt0ycODELFy7MXnvtlZdffjlXX311+vXrV9o/YABK\nyhlAAHZ673//+3PDDTeksbExXbp0yXXXXZdBgwbl8ccfT5K88cYbue6663Lttddm6tSpaWpqykUX\nXZSrrroqP/nJT3LwwQfn6quvTpJs2LAht9xyS84666wcdthh+ad/+qf07ds3r732Wq6//vrMmjUr\nmzZtylNPPZUkWbt2bX74wx/m6quvzowZM5IkF198cS655JLMmjUrH/vYx7J48eJcdNFFGT9+fG66\n6aYcfvjh+dGPfvS21zB58uTcc889Of744/P5z38+v//975MkL7zwQmbMmJFbb701gwYNygMPPJAk\nqa2tzfe+971MmDAhV199dS6//PJcc801+elPf5okWb9+fT71qU/l1ltvzX777dey/y+tW7cuxx13\nXH7yk5+kd+/emTt3bubMmZPXXnstt912Wy655JK88sorO/hPDYCdkTOAAOz09t9//yRJVVVVBg0a\n1PL1hg0bkiQHH3xwOnTokL322itVVVVZvnx5unXrlj59+rQ8Pn369BxxxBEZOHDg235+hw4d0qlT\np5x77rnZfffdU1tbm40bNyZJhgwZkiSpqalJY2NjkmT58uV597vfnSQZOXJkkmTx4sWZOHFikqSp\nqSn77rvvFsdYtGhRBg4cmOnTp6e5uTkPPvhgzj777Dz44IN517velXHjxmWPPfbI888/nwMPPDBJ\nMnjw4HTq1Cl77rlnBgwYkMrKynTv3r3ldXfs2DEHH3xwkmTYsGGZO3duy/f+pfe9730tr2PDhg15\n+eWXW57bq1ev7Lffftv5pwHArkwAArDLe/rpp5O8GWZr165N7969s3bt2tTX16d379557LHHWoKs\nQ4f/9+aXioqKNDc3Z9GiRbnnnnsya9asvPHGGznhhBPS3Nzc8py/1Lt377zwwgvZd999M2PGjAwc\nODADBw7MZZddlr333jvz589PQ0PDFt/z8MMP509/+lMmT56cDh06ZPDgwenatWvWrl2bf/u3f8t9\n992XJDn11FO3eez/aePGjVm0aFGGDBmS+fPnt8TxX/OXP2vw4MH55S9/mSRZtWpVXnjhhW0eC4D2\nQQACsMtbvnx5xo4dmzVr1mT8+PHZbbfdMmXKlHzjG99IRUVFunfvnksvvTTPPffcFt/3gQ98IFdc\ncUWmT5+erl27ZtSoUUnevIDKti4OM3HixJx//vnp0KFDqqurc8opp6Smpibjxo1r+Rzfv/zLv2zx\nPWPGjMlll12Wz3zmM+nWrVs6dOiQyy+/PN26dcuwYcNy0kknpWPHjqmqqkp9ff12fxbvmmuuyZ//\n/OfsvffeOeecc3LnnXdu1/cdccQRmTt3bkaNGpW99torXbp0SadOnbbrewHYdVU0v/XXjACwC/rF\nL36R559/Puedd165Rym5I488MnfddVc6d+78jr938eLFWbRoUY477risXLkyxx9/fP7zP/8zlZWV\nbTApADsLZwABoIBqampyxRVX5IYbbsimTZty3nnniT+AAnAGEAAAoCDcBgIAAKAgBCAAAEBBtKvP\nAK5fvz4LFy5MdXV1dtttt3KPAwAAUFKbNm1KQ0NDhg4dmi5durzt8XYVgAsXLszo0aPLPQYAAEBZ\n3XzzzRk+fPjb9rerAKyurk7y5ovt27dvmafZ0sKFCzN06NByjwGUkHUPxWLNQ/HsjOu+trY2o0eP\nbmmjv9SuAvCtt3327dt3u2+gWyp1dXU73UxA27LuoViseSienXndb+0jcS4CAwAAUBACEAAAoCAE\nIAAAQEEIQAAAgIIQgAAAAAUhAAEAAApCAAIAABSEAAQAACgIAQgAAFAQAhAAAKAgBCAAAEBBCEAA\nAICCEIAAAAAFIQABAAAKQgACAAAUhAAEAAAoCAEIAABQEAIQAACgIAQgAABAQQhAAACAghCAAAAA\nBSEAAQAACkIAAgAAFIQABAAAKAgBCAAAUBACEAAAoCAEIAAAQEEIQAAAgIIQgAAAAAUhAAEAAApC\nAAIAABSEAAQAACgIAQgAAFAQAhAAAKAgBCAAAEBBCEAAAICCEIAAAAAFIQABAAAKQgACAAAUhAAE\nAAAoCAEIAABQEAIQAACgIAQgAABAQQhAAACAghCAAAAABSEAAQAACkIAAgAAFIQABAAAKAgBCAAA\nUBACEAAAoCAEIAAAQEEIQAAAgIIQgAAAAAUhAAEAAApCAAIAABSEAAQAACgIAQgAAFAQHUtxkLq6\nukydOjXdu3fP4MGDM3r06CTJQw89lNtvvz3Nzc354he/mAMPPDAXXXRR9thjjzQ2NmbChAmZO3du\n5syZk02bNuWII47IUUcdVYqRAQAA2p2SnAGcOXNmxowZkwkTJuS+++5LU1NTkuS6667LlClTMnny\n5MyYMSOPPvpo+vfvn/PPPz+9evXKggULMmvWrPTo0SMbN27M+973vlKMCwAA0C6V5Azg8uXLU1NT\nkySpqqrKmjVr0qtXrzQ3N6eysjJJ0tTUlOXLl6dv375Jkr59+6a+vj5PP/10Lr/88ixfvjxXXnll\nLr/88laPt3DhwtTV1bXdC/obzZ8/v9wjACVm3UOxWPNQPDvbum9oaNjm4yUJwJqamtTW1qampiar\nVq1KVVVVkqRz585pbGzM5s2bU1lZmZqampZfYG1tbQYNGpS/+7u/S+fOndOjR4/tPt7QoUPTr1+/\nNnktf6v58+fnoIMOKvcYQAlZ91As1jwUz8647pctW7bNx0sSgCNHjszUqVMze/bsHH300bn00ksz\nbty4jB07NhdeeGE2btyYM888M0OHDs0dd9yRKVOmJEmGDRuW0aNHZ9y4cdm4cWO+9rWvlWJcAACA\ndqkkAVhdXZ1p06a9bf+IESMyYsSILfZNmjRpi+1jjjkmxxxzTJvOBwAAUARuAwEAAFAQAhAAAKAg\nBCAAAEBBCEAAAICCEIAAAAAFIQABAAAKQgACAAAUhAAEAAAoCAEIAABQEAIQAACgIAQgAABAQQhA\nAACAghCAAAAABSEAAQAACkIAAgAAFIQABAAAKAgBCAAAUBACEAAAoCAEIAAAQEEIQAAAgIIQgAAA\nAAUhAAEAAApCAAIAABSEAAQAACgIAQgAAFAQAhAAAKAgBCAAAEBBCEAAAICCEIAAAAAFIQABAAAK\nQgACAAAUhAAEAAAoCAEIAABQEAIQAACgIAQgAABAQQhAAACAghCAAAAABSEAAQAACkIAAgAAFIQA\nBAAAKAgBCAAAUBACEAAAoCAEIAAAQEEIQAAAgIIQgAAAAAUhAAEAAApCAAIAABSEAAQAACgIAQgA\nAFAQAhAAAKAgBCAAAEBBCEAAAICCEIAAAAAFIQABAAAKQgACAAAUhAAEAAAoCAEIAABQEAIQAACg\nIAQgAABAQQhAAACAghCAAAAABSEAAQAACkIAAgAAFIQABAAAKAgBCAAAUBACEAAAoCAEIAAAQEEI\nQAAAgIIQgAAAAAUhAAEAAApCAAIAABSEAAQAACgIAQgAAFAQAhAAAKAgBCAAAEBBCEAAAICC6FiK\ng9TV1WXq1Knp3r17Bg8enNGjRydJHnroodx+++1pbm7OF7/4xRx44IG56KKLsscee6SxsTETJkxI\nkqxZsyYnnXRSbrjhhlRXV5diZAAAgHanJGcAZ86cmTFjxmTChAm577770tTUlCS57rrrMmXKlEye\nPDkzZszIo48+mv79++f8889Pr169smDBgmzevDnTp0/PgAEDSjEqAABAu1WSM4DLly9PTU1NkqSq\nqipr1qxJr1690tzcnMrKyiRJU1NTli9fnr59+yZJ+vbtm/r6+lx11VU56aSTcv3112/38RYuXJi6\nurod/jr+t+bPn1/uEYASs+6hWKx5KJ6dbd03NDRs8/GSBGBNTU1qa2tTU1OTVatWpaqqKknSuXPn\nNDY2ZvPmzamsrExNTU3LL7C2tjb77rtv/vCHP2T58uX5r//6r/zoRz/KP//zP7d6vKFDh6Zfv35t\n+preqfnz5+eggw4q9xhACVn3UCzWPBTPzrjuly1bts3HSxKAI0eOzNSpUzN79uwcffTRufTSSzNu\n3LiMHTs2F154YTZu3JgzzzwzQ4cOzR133JEpU6YkSUaMGJERI0YkSb797W/nq1/9ainGBQAAaJdK\nEoDV1dWZNm3a2/b/z8B7y6RJk/7qz5g6dWqbzAYAAFAUbgMBAABQEAIQAACgIAQgAABAQQhAAACA\nghCAAAAABSEAAQAACkIAAgAAFIQABAAAKAgBCAAAUBACEAAAoCAEIAAAQEEIQAAAgIIQgAAAAAUh\nAAEAAApCAAIAABSEAAQAACgIAQgAAFAQAhAAAKAgBCAAAEBBCEAAAICCEIAAAAAFIQABAAAKQgAC\nAAAUhAAEAAAoCAEIAABQEAIQAACgIDpuz5MefvjhLF26NB/4wAcycODAdO7cua3nAgAAYAdrNQCn\nT5+e2traLF68OJWVlZkxY0amT59eitkAAADYgVp9C+j8+fNz+eWXZ/fdd8/nPve5LFu2rBRzAQAA\nsIO1GoCbNm3Khg0bUlFRkU2bNqVDBx8bBAAA2BW1+hbQsWPH5oQTTsiKFSsycuTInHrqqaWYCwAA\ngB2s1QD84Ac/mFtuuSUvvvhi+vXrl9dee60UcwEAALCDbfX9nM8++2zmzZuXf/iHf8hTTz2V1atX\nZ+HChTnnnHNKOR8AAAA7yFbPAK5evTq/+c1v8uqrr+bXv/51kqSioiInn3xyyYZrL5bWrc7Di9ak\nut/qDOhTVe5xAACAgtpqAA4fPjzDhw/P008/nf3337+UM7UrS+tWZ+I1j6R+5Rv5ryWPZPxph4hA\nAACgLFr9DGBtbW2mT5+epqamNDc357XXXsuvfvWrUszWLixYVJ/6lW8kSepXvpEFi+oFIAAAUBat\n3tPhyiuvzFlnnZWampp87nOfy3vf+95SzNVuDBvSO717dk2S9O7ZNcOG9C7zRAAAQFG1GoC9e/fO\nBz/4wSTJCSeckLq6ujYfqj0Z0Kcq4087JJ8Y1t3bPwEAgLJqNQA7deqUxx9/PBs3bsy8efOycuXK\nUswFAADADtZqAE6cODEbN27M1772tfzsZz/LmWeeWYq52o23LgJz94JVmXjNI1lat7rcIwEAAAXV\n6kVg9tprr6xcuTKvv/56vvrVr6aioqIUc7UbLgIDAADsLFoNwNNPPz2NjY2pqnozWioqKnLVVVe1\n+WDtxbAhvfOrec+nfuUbLgIDAACUVasBuGHDhvzkJz8pxSzt0lsXgbnjngX59FHDnP0DAADKptUA\nHD58eObNm5d3v/vdLfv23nvvNh2qvRnQpyofHrKn+AMAAMqq1QB89dVXc8kll2zxFtCZM2e2+WAA\nAADsWK0G4PPPP5+77rqrFLMAAADQhlq9DcR73/ve/P73v09jY2PL/3hnltatzsOL1rgFBAAAUFat\nngF8/PHHc99997VsV1RUZM6cOW05U7vy1n0A61e+kf9a8kjGn3aIzwICAABl0WoA/upXvyrFHO2W\n+wACAAA7i60G4KRJk3LxxRfnpJNOetvN310EZvsNG9I7s//zuaxY05hee1a6DyAAAFA2Ww3AM888\nM0ly2WWXpVOnTi37V61a1fZTtTdvBfRfhDQAAEApbfUiMM3NzVmyZEm+9a1vpampKY2NjVm/fn0u\nvvjiUs63y1uwqD4rVm9IkqxYvSELFtWXeSIAAKCotnoG8A9/+ENuuOGGLFmyJBdffHGam5vToUOH\nfPSjHy3lfLu8YUN651fznk/9yjfSu2dXbwEFAADKZqsBeNRRR+Woo47K/fffn4997GOlnKldGdCn\nKv/wuQPys7ufzImfOMAFYAAAgLJp9T6AnTp1yty5c3P//ffnqKOOclXQd2hp3er8cPZTeebl9fnh\n7KfcCxAAACibVgPwX//1X7PvvvvmxhtvzK233uoKoO/QX7sNBAAAQDm0GoBdunTJu971rnTs2DHV\n1dVvuyUE27b3Xnukw3//yjpUvLkNAABQDq0GYLdu3fLVr341n/zkJ3PzzTenV69epZir3fjz8nXZ\n3Pzm15ub39wGAAAoh61eBOYt3/3ud7N06dIMGjQozz77bEaOHFmKudoNVwGFYlq26pU8vvKp9Fm1\nd/p1ryn3OAAASbZxBvDnP/95kqSysjLNzW+ewnrPe96TGTNmlGaydmJAn6qMP+2QfGJY94w/7RBX\nAYUCWLbqlYyfMy33vvpoxs+ZlmWrXin3SAAASbYRgL/85S9bvp48eXLL14899ljbTtQODehTlQ8P\n2VP8QUHc8cx/ZE3Tm2/3XtO0Lnc88x9lnggA4E1bDcC3zvpt62sA3u4vz/g5AwgA7Cy2GoD/82qf\nW/ua7bO0bnUeXrTGPQChIHp26b7Fdq8uPco0CQDAlrZ6EZjXXnstDz74YDZv3pxVq1blgQceSHNz\nc1atWlXK+XZ5S+tW56IfPJQVqzfk8T89lMlnfMRbQaGdq6nqnfz5/233raou3zAAAP/DVgNw//33\nz5133pkked/73pdf//rXLV+z/e59/KWsWL0hSbJi9Ybc+/hLOeX4/cs8FdCWOlbsts1tAIBy2WoA\nXnrppaWco91avXb9NreB9mfO4gfftj3q/Z8p0zRAKfz5rt9m/W2z8+cvfC57f/KYco8DsFWt3geQ\n/52qbl22uQ20Pxs2bviL7cYyTQKUwp/v+m2W/OCaJGn5fxEI7Ky2ehEYdozddtv2NtD+dOzQ8S+2\nLXxoz166bfY2t4H26d67ns7vfv7n3HvX0+Ue5R1pNQAnTZqUP/7xj6WYpV165KnabW4D7c8+Pfpv\nsb1vjwFlmgQohY1r1mxzG2h/7r3r6Txwz/Np2pA8cM/zu1QEthqARxxxRH7wgx9k1KhRueWWW7J2\n7dpSzNVuHHvovtvcBtqft90ux91zoH3bsGHb20C788A9z29ze2fWagAefvjh+e53v5vvf//7mT9/\nfg477LB8+9vfztKlS0sx3y6vusfuLf/tV/Hf20D71vD68m1uAwCUS6sBuHjx4nznO9/Jl770pVRV\nVeXmm2/OySefnLPPPrsU8+3yfvfoi2n+76+b/3sbaN/67F69zW0AgHJp9SqgF154YUaOHJmzzjor\nXbt2bdn/+c9/vk0Hay8++N7qPPp03RbbQPv24qpl29wG2pkuXZL167fcBtq1jx613xZv+/zoUfuV\ncZp3ptUzgIcddlhOOOGElvibNm1akmT06NFtO1k70bSxeZvbQPvz8Xcfus1toH2p/thh29wG2p8j\nP7l/PnrUfunU+c34O/KT+5d7pO221TOAs2bNym233ZbFixdn7ty5SZJNmzZl48aN+cd//MeSDbir\nGzakd3417/nUr3wjvXt2zbAhvcs9EtDG3rrp+93P3JdPvPcIN4GHdm7jyte23H7tta08E2hPjvzk\n/unee30OOmjXib9kGwH4mc98Jh/+8Ifzwx/+MGeccUaSpEOHDnnXu95VsuHagwF9qjL+tENyxz0L\n8umjhmVAn6pyjwSUwEf3GZHX6lfmo/uMKPcoQBvrc/RRWfnE/GTz5qRDh/T5/44q90gAW7XVAHzm\nmWdywAEH5Oijj86SJUta9i9evDgf/ehH39FB6urqMnXq1HTv3j2DBw9uefvoQw89lNtvvz3Nzc35\n4he/mAMPPDAXXXRR9thjjzQ2NmbChAm55ZZb8uyzz+b111/PMccckyOPPPJvfKnlM6BPVT48ZE/x\nBwWxbNUrmfyfV2blhtX5r//8Yy76+7PTr3tNuccC2si7Dh6eIeePy7Ozbst7Rn4h7zp4eLlHAtiq\nrQbgww8/nAMOOCC/+c1v3vbYOw3AmTNnZsyYMRk2bFhOO+20nHjiienUqVOuu+66/Pu//3s2b96c\ns88+O2PHjk3//v1zxhln5N/+7d+yYMGCVFVVZcKECVm5cmUmT568SwYgUCz3vfBwVm5YnSRZuWF1\n7nvh4XzpAyeUeSqgLb3r4OGp7FCRdx10ULlHAdimrQbgKaecksbGxkycOPF/fZDly5enpubNv/2u\nqqrKmjVr0qtXrzQ3N6eysjJJ0tTUlOXLl6dv375Jkr59+6a+vj7HH3981q1bl8suuyynn376/3qW\nclhatzoPL1qT6n6rnQWEIviLaz1VNLsTPACwc9hqAB5zzDGpqNjyP1qam5tTUVGROXPmvKOD1NTU\npLa2NjU1NVm1alWqqt6MoM6dO6exsTGbN29OZWVlampqMn/+/CRJbW1tBg0alMWLF+fqq6/ON7/5\nzfTv33+7jrdw4cLU1dW1/sQSqF/VmGt/15D1Tc25/6n78pWjq9O7e2W5xwLaUPWGqnTbbfes3fR6\nuu22e/Z6fc+Wf7YB7Zu1DsWzs637hoaGbT5e0dzc3Ob3JWhoaMjUqVOzxx57ZOjQoXnmmWcybty4\n/P73v89tt92WjRs35tRTT83QoUMzfvz4lrOCF1xwQY4++ugMGTIknTt3znve855tngVctmxZPv7x\nj2fOnDnp169fW7+s7XLJ9Y/m4adqW7Y/fEDfnH/Kh8o4EVAKy1a9kjsfvzvHH/wJn/+Dgpg/f34O\n8hZQKJSdcd231kRbPQM4adKkXHzxxTnppJPediZw5syZ72iI6urqlvsH/k8jRozIiBFbXiFv0qRJ\nW2z/x3/8xzs61s5myZ9XbXMbaJ/6da/JwT0PEH8AwE5lqwF45plnJkmmT59esmHao89+7N35wS8W\nbrENAABQDlsNwL322itJsnnz5lx++eV54YUXMnjw4PzTP/1TyYZrD4479N1ZsXp9fj3v+Rx32H45\n7lABCAAAlEeH1p5w/vnn5wtf+EJuueWWHH/88Tn//PNLMVe7sbRude55dGnWbdicex5dmqV1q8s9\nEgAAUFCtBuBuu+2Wj33sY9lzzz1z5JFHZvPmzaWYq92497GlWbGmMUmyYk1j7n1saZknAgAAimqr\nbwF94IEHkiRdu3bNNddck4MPPjhPPvlky1tD2T6rX2/a5jYAAECpbDUAf/3rXydJevTokeeffz7P\nP/98krTcooHt1Lx529sAAAAlstUAvPTSS//q/vr6+jYbpj2q6tZlm9sAAACl0upnAL/73e/mkEMO\nyUEHHZT9998/p556ainmajeOPLh/elV1TpL0quqcIw/uX+aJAACAomo1AO+9997MnTs3n/rUp/Kb\n3/wmffr0KcVc7caAPlX5+hc+kPf+XZd8/QsfyIA+VeUeCQAAKKhWA7C6ujqVlZVZt25d9tlnnzQ1\nuYjJO7G0bnV+OPupPPPy+vxw9lNuAwEAAJRNqwHYt2/f3HbbbenatWumTZuW1asFzDuxYFF96le+\nkSSpX/lGFizyGUoAAKA8Wg3ASZMm5SMf+Ui+9a1vpXfv3pk2bVop5mo3hg3pnd49uyZJevfsmmFD\nepd5IgAAoKi2ehXQt6xatSo33nhjXnjhhQwePNhnAN+hAX2qMv60Q3LHPQvy6aOG+QwgAABQNq2e\nARw3blwGDBiQs88+O3369Mm4ceNKMVe7MqBPVT48ZE/xBwAAlFWrZwA3bNiQk08+OUkyZMiQ3H33\n3W0+FAAAADveVgNwyZIlSZKePXvmrrvuyvDhw/Pkk0+mX79+JRsOAACAHWerAXjxxRe3fH3LLbfk\n1ltvTXNzcyoqKkoyGAAAADvWVgPwpptuavl65cqVeemll9KvX7/06tWrJIMBAACwY7V6EZi77ror\no0aNyg9+8IOcdNJJ+eUvf1mKuQAAANjBWr0IzPXXX59f/OIX2WOPPbJ27dqMHTs2n/nMZ0oxW7ux\ntG51Hl60JtX9VrsSKAAAUDatngGsqKjIHnvskSTp1q1bOnfu3OZDtSdL61Zn4jWP5O4FqzLxmkey\ntG51uUcCAAAKqtUzgP3798/UqVMzfPjwPPHEExkwYEAp5mo3FiyqT/3KN5Ik9SvfyIJF9c4CAgAA\nZdHqGcB/+Zd/Sf/+/fPQQw+lf//+mTx5cinmajeGDemd3j27Jkl69+yaYUN6l3kiAACgqFo9A3jG\nGWfk2muvLcUs7dKAPlUZf9ohueOeBfn0UcOc/QMAAMqm1QCsqqrKnDlzsu+++6ZDhzdPGA4cOLDN\nB2tPBvSpyoeH7Cn+AACAsmoL+ec4AAAQ5klEQVQ1AF999dVcf/31LdsVFRW58cYb23ImAAAA2sA2\nA3Dt2rWZMWNGunbtWqp5AAAAaCNbvQjMT37yk3z605/OZz7zmcybN6+UMwEAANAGthqAd955Z377\n299m5syZueGGG0o5EwAAAG1gqwFYWVmZysrK9OrVK01NTaWcCQAAgDbQ6n0Ak6S5ubmt5wAAAKCN\nbfUiMH/605/yj//4j2lubm75+i3Tpk0ryXAAAADsOFsNwCuvvLLl61GjRpVkGAAAANrOVgNwxIgR\npZwDAACANrZdnwEEAABg1ycAAQAACkIAAgAAFIQABAAAKAgBCAAAUBACEAAAoCAEIAAAQEEIQAAA\ngIIQgAAAAAUhAAEAAApCAAIAABSEAAQAACgIAQgAAFAQAhAAAKAgBCAAAEBBCEAAAICCEIAAAAAF\nIQABAAAKQgACAAAUhAAEAAAoCAEIAABQEAIQAACgIAQgAABAQQhAAACAghCAAAAABSEAAQAACkIA\nAgAAFIQABAAAKAgBCAAAUBACEAAAoCAEIAAAQEEIQAAAgIIQgAAAAAUhAAEAAApCAAIAABSEAAQA\nACgIAQgAAFAQAhAAAKAgBCAAAEBBCEAAAICCEIAAAAAFIQABAAAKQgACAAAURMdSHKSuri5Tp05N\n9+7dM3jw4IwePTpJ8tBDD+X2229Pc3NzvvjFL+bAAw/MRRddlD322CONjY2ZMGHC254zbNiwUowM\nAADQ7pTkDODMmTMzZsyYTJgwIffdd1+ampqSJNddd12mTJmSyZMnZ8aMGXn00UfTv3//nH/++enV\nq1cWLFjwtucAAADwtynJGcDly5enpqYmSVJVVZU1a9akV69eaW5uTmVlZZKkqakpy5cvT9++fZMk\nffv2TX19/duesz0WLlyYurq6Nngl/zvz588v9whAiVn3UCzWPBTPzrbuGxoatvl4SQKwpqYmtbW1\nqampyapVq1JVVZUk6dy5cxobG7N58+ZUVlampqam5RdYW1ubQYMGve0522Po0KHp169fm72ev8X8\n+fNz0EEHlXsMoISseygWax6KZ2dc98uWLdvm4yUJwJEjR2bq1KmZPXt2jj766Fx66aUZN25cxo4d\nmwsvvDAbN27MmWeemaFDh+aOO+7IlClTkiTDhg3Lxo0bt3gOAAAAf5uSBGB1dXWmTZv2tv0jRozI\niBEjttg3adKkVp8DAADAO+c2EAAAAAUhAAEAAApCAAIAABSEAAQAACgIAQgAAFAQAhAAAKAgBCAA\nAEBBCEAAAICCEIAAAAAFIQABAAAKQgACAAAUhAAEAAAoCAEIAABQEAIQAACgIAQgAABAQQhAAACA\nghCAAAAABSEAAQAACkIAAgAAFIQABAAAKAgBCAAAUBACEAAAoCAEIAAAQEEIQAAAgIIQgAAAAAUh\nAAEAAApCAAIAABSEAAQAACgIAQgAAFAQAhAAAKAgBCAAAEBBCEAAAICCEIAAAAAFIQABAAAKQgAC\nAAAUhAAEAAAoCAEIAABQEAIQAACgIAQgAABAQQhAAACAghCAAAAABSEAAQAACkIAAgAAFIQABAAA\nKAgBCAAAUBACEAAAoCAEIAAAQEEIQAAAgIIQgAAAAAUhAAEAAApCAAIAABSEAAQAACgIAQgAAFAQ\nAhAAAKAgBCAAAEBBCEAAAICCEIAAAAAFIQABAAAKQgACAAAUhAAEAAAoCAEIAABQEAIQAACgIAQg\nAABAQQhAAACAghCAAAAABSEAAQAACkIAAgAAFIQABAAAKAgBCAAAUBACEAAAoCAEIAAAQEEIQAAA\ngIIQgAAAAAUhAAEAAApCAAIAABSEAAQAACiIjqU4yO23354nnngib7zxRs4666wMHDiw5bHLL788\njY2NWbt2bSZNmpTnn38+M2bMSJcuXfL3f//3OfzwwzNhwoT06NEjr7zySsvXAAAAvDMlC8Drr78+\nL730UmbMmJHJkycnSV566aWW8Pv5z3+eu+++O/Pmzcu4cePSp0+ffOUrX8kBBxyQz3/+8xk+fHiu\nvfbaPPXUUznssMNKMTYAAEC70iYB+NOf/jR33nlny3bnzp2TJH379k19fX3L/oaGhvTp06flsUWL\nFmXFihUt+yoqKtK3b9/07ds3Tz75ZBYtWpQvf/nLrR5/4cKFqaur25EvaYeYP39+uUcASsy6h2Kx\n5qF4drZ139DQsM3H2yQATzrppJx00kkt26effnqSpLa2Nr17927ZX1NT0xJqbz3Wp0+f1NfXp3fv\n3mlubk6SzJw5M7W1tbnkkkvSsWPrIw8dOjT9+vXbkS/pf23+/Pk56KCDyj0GUELWPRSLNQ/FszOu\n+2XLlm3z8ZK8BfSzn/1sLrjggqxduzbf+ta3UldXl1mzZuWss85Kjx49MmXKlKxbty6TJk3K//k/\n/yeXXXZZOnXqlDFjxuSxxx7Lv//7v+dDH/pQvv3tb+fkk0/OsGHDSjE2AABAu1KSADz22GNz7LHH\nbrHvrLPOSpKce+65W+wfNGhQpk2btsW+efPmte2AAAAABeA2EAAAAAUhAAEAAApCAAIAABSEAAQA\nACgIAQgAAFAQAhAAAKAgBCAAAEBBCEAAAICCEIAAAAAFIQABAAAKQgACAAAUhAAEAAAoCAEIAABQ\nEAIQAACgIAQgAABAQQhAAACAghCAAAAABSEAAQAACkIAAgAAFIQABAAAKAgBCAAAUBACEAAAoCAE\nIAAAQEEIQAAAgIIQgAAAAAUhAAEAAApCAAIAABSEAAQAACgIAQgAAFAQAhAAAKAgBCAAAEBBCEAA\nAICCEIAAAAAFIQABAAAKQgACAAAUhAAEAAAoCAEIAABQEAIQAACgIAQgAABAQXQs9wA70qZNm5Ik\ntbW1ZZ7k7RoaGrJs2bJyjwGUkHUPxWLNQ/HsjOv+rRZ6q43+UrsKwIaGhiTJ6NGjyzwJAABA+TQ0\nNGSfffZ52/6K5ubm5jLM0ybWr1+fhQsXprq6Orvttlu5xwEAACipTZs2paGhIUOHDk2XLl3e9ni7\nCkAAAAC2zkVgAAAACkIAAgAAFIQABAAAKAgBCAAAUBACEAAAoCDa1X0AdxX33HNP6uvrkyQnn3xy\nmacBSuVnP/tZKioqMnLkyHKPArSxWbNmpa6uLoMGDcoxxxxT7nGAEpg3b14WLVqUysrKjB07ttzj\nbJUzgGVw1FFHpXv37tmwYUO5RwFK5Iknnvir9+IB2qfFixenW7du6djR37VDUTz22GPp3Llz+vfv\nX+5RtkkAlsHixYtz3HHHZeXKleUeBSiRefPmZcmSJVmwYEG5RwFKYOzYsTnllFPyxBNPlHsUoETW\nrVuXL3/5y3nsscfKPco2+WupMnjqqacyZ86c9O3bt9yjACVyzjnnZNmyZXn44YfLPQpQAvfcc082\nbdqUffbZp9yjACVy6KGH5sYbb0xNTU25R9mmiubm5uZyD9HevPjii/nmN7+Z22+/PXV1dZk6dWq6\nd++ewYMHZ/To0eUeD2gD1j0UizUPxdNe1r23gO5gDQ0NmTVrVrp27ZokmTlzZsaMGZMJEybkvvvu\nS1NTU5knBHY06x6KxZqH4mlP614A7mDV1dU577zzsvvuuydJli9f3nIauKqqKmvWrCnneEAbsO6h\nWKx5KJ72tO4FYBurqalJbW1tkmTVqlWpqqoq80RAW7PuoViseSieXXnduwhMGxs5cmSmTp2a2bNn\n5+ijj3Y5aCgA6x6KxZqH4tmV172LwAAAABSEt4ACAAAUhAAEAAAoCAEIAABQEAIQAACgIAQgAABA\nQQhAAACAghCAAAAABSEAAdipPfrooznooIPyyiuvtOy74oor8otf/OJv/pnLli3LiSeeuCPGe5uN\nGzdmzJgxGTVqVFatWtWyf8WKFfnGN76Rr3zlKxk1alQuuOCCrF+/focf/9vf/nbmzp27w38uAO2D\nAARgp1dZWZl//ud/TnNzc7lHaVV9fX3WrVuXmTNnpnv37i37f/SjH+UjH/lIrr322sycOTO77757\nZs6cWcZJASiijuUeAABac8ghh2Tz5s25+eab86Uvfall/7Jly3LuuefmZz/7WZLkxBNPzPTp0zN7\n9uy8+OKLWblyZV577bWMHj06v/vd77JkyZJcdtll2WuvvbJixYqcccYZefXVV3PEEUfk61//el55\n5ZVcdNFF2bBhQzp37pzJkydn06ZN+drXvpYePXrk8MMPz2mnndZy/DvuuCM33HBDKisrs++++2bS\npEkZP358XnjhhVx88cWZNGlSy3P32muv3H333dlnn30ybNiwjBs3LhUVFUmSadOmZeHChXnttdcy\nZMiQXHrppfne977X6mv45je/merq6tTV1eXwww/POeec03K8pqamjB8/Pi+++GI2b96cs88+Ox/6\n0Ifyr//6r3n00UezcePGHH300Tn99NPb+o8PgJ2IM4AA7BImTJiQ66+/Pi+++OJ2Pb9Lly758Y9/\nnE984hO5//7784Mf/CCnn356fv3rXydJXn/99XznO9/JzJkzM2/evCxatCiXXXZZxowZk5tuuin/\n9//+31xxxRVJkoaGhvz4xz/eIv5WrlyZ733ve7nhhhty6623Zs8998xPf/rTjB8/PoMGDdoi/pLk\nlFNOyfHHH58f//jHOeyww3LWWWelvr4+a9euTVVVVa677rr8/Oc/z+9///vU1dVt12t4+eWXM3Xq\n1Nx222155JFH8vTTT7ccb9asWenZs2duvvnmfP/732+Z51e/+lWuuOKK3HLLLamqqvob/zQA2FU5\nAwjALqFnz545//zzM27cuAwbNuyvPud/vkX0fe97X5Jkzz33zKBBg5Ik3bt3z4YNG5IkQ4YMyZ57\n7pkkOeCAA7JkyZI8++yz+eEPf5gf/ehHaW5uTseOb/5rsl+/fqmsrNziWC+99FIGDRqUbt26JUkO\nPvjgPPDAAzniiCP+6myPPPJIPvvZz+YLX/hCGhsbc8011+SSSy7J9OnTs2LFipx77rnZfffd8/rr\nr6epqWm7X0OPHj2SJO9///uzZMmSluM9++yzmT9/fp588skkb342ccWKFfnOd76TadOmZfny5Tns\nsMO2/UsHoN1xBhCAXcaRRx6ZgQMHZvbs2UmSzp0759VXX82mTZuyevXqLFu2rOW5b729cmsWL16c\ndevWZePGjXnyySczePDg7LfffjnvvPNy0003ZeLEiTnmmGOSJB06vP1fl/369cvixYvz+uuvJ0ke\ne+yxDBw4cKvHu/HGG3PnnXcmefMzjYMHD05lZWXmzp2bV155JdOnT8+5556b9evXt4Ts9ryGN954\nI5s2bcqTTz7ZEolJst9+++W4447LTTfdlGuuuSbHHHNMunXrlt/+9reZPn16brzxxsyePTsvv/zy\nNo8BQPviDCAAu5QLLrggjzzySJKkuro6hx56aL7whS+kf//+2Weffbb753Tv3j3nnHNOVqxYkWOP\nPTaDBg3KuHHjMmHChGzYsCHr16/PBRdcsNXv79WrV77xjW/ky1/+cjp06JABAwbkvPPOS0NDw199\n/sSJEzNx4sRcf/316dKlS3r27JkJEyakQ4cO+f73v5/Ro0enoqIi/fv3T319/Xa9hk6dOuWb3/xm\nli9fnmOOOSZDhgxpeWzUqFG58MIL86UvfSlr167NySefnMrKynTv3j0nnnhiunTpkkMPPTR77733\ndv/OANj1VTTvCpdUAwC28JcXwAGA7eEtoAAAAAXhDCAAAEBBOAMIAABQEAIQAACgIAQgAABAQQhA\nAACAghCAAAAABSEAAQAACuL/By0lv9jrPplmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108197588>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "y_axis = np.asarray(total_probas)\n",
    "x_axis = np.asarray(num_samples)\n",
    "\n",
    "for x, y in zip(x_axis, y_axis):\n",
    "    plt.scatter([x] * len(y), y, s=12)\n",
    "    \n",
    "plt.xscale('log')\n",
    "plt.title(\"Importance Sampling\")\n",
    "plt.ylabel(\"Probability Estimate\")\n",
    "plt.xlabel('Number of Samples')\n",
    "plt.grid(True)\n",
    "plt.savefig(plot_dir + 'imp_sampling_convergence_speed.pdf', format='pdf', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "* http://ib.berkeley.edu/labs/slatkin/eriq/classes/guest_lect/mc_lecture_notes.pdf"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}