10_01_montecarlo_integration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Monte Carlo integration\n",
    "=======================\n",
    "\n",
    "Imagine that we want to measure the area of a pond with arbitrary shape.\n",
    "Suppose that this pond is in the middle of a field with known area $A$.\n",
    "If we throw $N$ stones randomly, such that they land within the\n",
    "boundaries of the field, and we count the number of stones that fall in\n",
    "the pond $N_{in}$, the area of the pond will be approximately\n",
    "proportional to the fraction of stones that make a splash, multiplied by\n",
    "$A$: $$A_{pond}=\\frac{N_{in}}{N}A.$$ This simple procedure is an example\n",
    "of the “Monte Carlo” method.\n",
    "\n",
    "Simple Monte Carlo integration\n",
    "------------------------------\n",
    "\n",
    "More generaly, imagine a rectangle of height $H$ in the integration\n",
    "interval $[a,b]$, such that the function $f(x)$ is within its\n",
    "boundaries. Compute $n$ pairs of random numbers $(x_i,y_i)$ such that\n",
    "they are uniformly distributed inside this rectangle. The fraction of\n",
    "points that fall within the area contained below $f(x)$, <span>*i.\n",
    "e.*</span>, that satisfy $y_i \\leq f(x_i)$ is an estimate of the ratio o\n",
    "fthe integral of $f(x)$ and the area of the rectangle. Hence, the\n",
    "estimate of the integral will be given by:\n",
    "$$\\int _a^b{f(x)dx} \\simeq I(N) = \\frac{N_{in}}{N}H(b-a).\n",
    "$$\n",
    "\n",
    "Another Monte Carlo procedure is based on the definition:\n",
    "$$\\langle g \\rangle=\\frac{1}{(b-a)} \\int _a^b{f(x)dx}.\n",
    "$$ In order to determine this average, we sample the\n",
    "value of $f(x)$:\n",
    "$$\\langle f \\rangle \\simeq \\frac{1}{N}\\sum_{i=1}^{N}f(x_i),$$ where the\n",
    "$N$ values $x_i$ are distributed unformly in the interval $[a,b]$. The\n",
    "integral will be given by $$I(N)=(b-a) \\langle f \\rangle .$$\n",
    "Monte Carlo error analysis\n",
    "--------------------------\n",
    "\n",
    "The Monte Carlo method clearly yields approximate results. The accuracy\n",
    "deppends on the number of values $N$ that we use for the average. A\n",
    "possible measure of the error is the “variance” $\\sigma^2$ defined by:\n",
    "$$\\sigma ^2=\\langle f^2 \\rangle - \\langle f \\rangle ^2,\n",
    "$$ where\n",
    "$$\\langle f \\rangle = \\frac{1}{N} \\sum_{i=1}^N f(x_i)$$ and\n",
    "$$\\langle f^2 \\rangle = \\frac{1}{N} \\sum_{i=1}^{N} f(x_i)^2.$$ The\n",
    "“standard deviation” is $\\sigma$. However, we should expect that the\n",
    "error decreases with the number of points $N$, and the quantity $\\sigma$\n",
    "defines by (\\[mc\\_sigma\\]) does not. Hence, this cannot be a good\n",
    "measure of the error.\n",
    "\n",
    "Imagine that we perform several measurements of the integral, each of\n",
    "them yielding a result $I_n$. These values have been obtained with\n",
    "different sequences of $N$ random numbers. According to the central\n",
    "limit theorem, these values whould be normally dstributed around a mean\n",
    "$\\langle I\n",
    "\\rangle$. Suppouse that we have a set of $M$ of such measurements\n",
    "${I_n}$. A convenient measure of the differences of these measurements\n",
    "is the “standard deviation of the means” $\\sigma_M$:\n",
    "$$\\sigma_M ^2=\\langle I^2 \\rangle - \\langle I \\rangle ^2,\n",
    "$$ where\n",
    "$$\\langle I \\rangle = \\frac{1}{M} \\sum_{n=1}^M I_n$$ and\n",
    "$$\\langle I^2 \\rangle = \\frac{1}{M} \\sum_{n=1}^{M} I_n^2.$$ \n",
    "It can be proven that\n",
    "$$\\sigma_M \\approx \\sigma/\\sqrt{N}.\n",
    "$$ This relation becomes exact in the limit of a very\n",
    "large number of measurements. Note that this expression implies that the\n",
    "error decreases with the squere root of the number of trials, meaning\n",
    "that if we want to reduce the error by a factor 10, we need 100 times\n",
    "more points for the average.\n",
    "\n",
    "### Exercise 10.1: One dimensional integration \n",
    "\n",
    "1.  Write a program that implements the “hit and miss” Monte Carlo\n",
    "    integration algorithm. Find the estimate $I(N)$ for the integral of\n",
    "    $$f(x)=4\\sqrt{1-x^2}$$ as a function of $N$, in the interval\n",
    "    $(0,1)$. Choose $H=1$, and sample only the $x$-dependent part\n",
    "    $\\sqrt{1-x^2}$, and multiply the result by 4. Calculate the\n",
    "    difference between $I(N)$ and the exact result $\\pi$. This\n",
    "    difference is a measure of the error associated with the Monte\n",
    "    Carlo estimate. Make a log-log plot of the error as a function of\n",
    "    $N$. What is the approximate functional deppendece of the error on\n",
    "    $N$ for large $N$?\n",
    "\n",
    "2.  Estimate the integral of $f(x)$ using the simple Monte Carlo\n",
    "    integration by averaging over $N$ points, using (\\[mc\\_integral2\\]),\n",
    "    and compute the error as a function of $N$, for $N$ upt to 10,000.\n",
    "    Determine the approximate functional deppendence of the error on $N$\n",
    "    for large $N$. How many trials are necessary to determine $I_N$ to\n",
    "    two decimal places?\n",
    "\n",
    "3.  Perform 10 measurements $I_n(N)$, with $N=10,000$ using different\n",
    "    random sequences. Show in a table the values of $I_n$ and $\\sigma$\n",
    "    according to (\\[mc\\_integral2\\]) and (\\[mc\\_sigma\\]).\n",
    "    Use (\\[mc\\_sigmam\\]) to estimate the standard deviation of the\n",
    "    means, and compare to the values obtained from (\\[mc\\_sigma2\\])\n",
    "    using the 100,000 values.\n",
    "\n",
    "4.  To verify that your result for the error is independent of the\n",
    "    number of sets you used to divide your data, repeat the previous\n",
    "    item grouping your results in 20 groups of 5,000 points each.\n",
    "\n",
    "### Exercise 10.2: Importance of randomness \n",
    "\n",
    "To examine the effects of a poor random number generator, modify your\n",
    "program to use the linear congruential random number generator using the\n",
    "perameters $a=5$, $c=0$ and the seed $x_1=1$. Repeat the integral of the\n",
    "previous exercise and compare your results.\n",
    "\n",
    "#### Challenge 10.1: \n",
    "\n",
    "Exercise 10.2\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 75 3.0 0.14159265359\n",
      "200 159 3.18 0.0384073464102\n",
      "400 318 3.18 0.0384073464102\n",
      "800 636 3.18 0.0384073464102\n",
      "1600 1269 3.1725 0.0309073464102\n",
      "3200 2500 3.125 0.0165926535898\n",
      "6400 5045 3.153125 0.0115323464102\n",
      "12800 10017 3.1303125 0.0112801535898\n",
      "25600 20188 3.154375 0.0127823464102\n",
      "51200 40219 3.142109375 0.000516721410207\n",
      "102400 80318 3.137421875 0.00417077858979\n",
      "204800 161064 3.14578125 0.00418859641021\n",
      "409600 321654 3.14115234375 0.000440309839793\n",
      "819200 643676 3.14294921875 0.00135656516021\n",
      "1638400 1286191 3.14011474609 0.00147790749604\n",
      "3276800 2575084 3.14341308594 0.00182043234771\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEHCAYAAABGNUbLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYFNW9//H3gWEQWQRNVEB4iBFlEYkmV1RQBzBKNIJK\njCIoboAM4pKfuTHROCDmGp/oNahsohjUKyoBBfXmGkGHgEvUKIgBUYjIMGxG2QRZhjm/P8401TM9\nSy9VvdXn9TzzUFtXny7bb1Wf5XuMtRYREQmXRpkugIiIpJ+Cv4hICCn4i4iEkIK/iEgIKfiLiISQ\ngr+ISAgp+IuIhJCCv4hICBUEeXJjTHNgMrAXKLXWPhPk+4mISHyCfvK/BHjeWjsSGBjwe4mISJwS\nDv7GmBnGmM3GmOU1tg8wxnxijPnMGPOrqs3tgbKq5QMpllVERHySzJP/E8CA6A3GmMbAI1XbuwFD\njDFdgfVAhxTeS0REApBwQLbWLga21th8KrDaWrvWWrsfeBYYBMwFBhtjJgPzUy2siIj4w68G3+jq\nHXBP/L2stbuBa+t7oTFGaUVFRJJgrTXJvtavqpiUAri11te/kpISX4+t65h4t9e3XteyrkV+XYt4\ntulapP9aJHq+bLoWqfIr+Jfj1e1Ttbw+3hePGzeO0tJSn4oCRUVFvh5b1zHxbq9vPZGyJkPXou73\nTvXYRK5FPNt0LWpfD/JaJHrubLgWpaWljBs3rsFyNCiZuyXQCVgetV4ArKnaXggsBbrGeS4rTklJ\nSaaLkDV0LTy6Fh5dC09V7Ez6V08yXT1nAW8Bxxtjyowx11hrK4AbgVeBFcBz1tqVKd6XQifop71c\nomvh0bXw6Fr4x1gf6o5SKoAxtqSkhKKiIv2HFRFpQGlpKaWlpYwfPx6bQoNvVgT/TJdBRCTXGGNS\nCv4aeCUiEkJZEfz97u0jIpKv/Orto2ofEZEcpGofERFJmIK/iEgIZUXwV52/iEh8VOcvIhJiqvMX\nEZGEKfiLiISQgr+ISAhlRfBXg6+ISHzU4CsiEmJq8BURkYQp+IuIhJCCv4hICGVF8FeDr4hIfNTg\nKyISYmrwFRGRhCn4i4iEkIK/iEgIKfiLiISQgr+ISAhlRfBXV08Rkfioq6eISIipq6eIiCRMwV9E\nJIQU/EVEQkjBX0QkhBT8RURCSMFfRCSEFPxFREJIwV9EJISyIvhrhK+ISHw0wldEJMQ0wldERBKm\n4C8iEkIK/iIiIaTgLyISQtkR/N9+O9MlEBEJlezo7dO9O3zwARQWZrQsIiK5Ij96+/zzn3DffZku\nhYhIaGTHkz+4p/5ly6BLl4yWR0QkF+THkz/Avn0wciRUVma6JCIieS87gn9Bgft38WJ47LHMlkVE\nJAQCDf7GmO8ZYx4zxsyu98Bf/tJb/s//hA0bgiyWiEjopaXO3xgz21p7aR37rN29G046CVavdhsv\nuQTmzAm8XCIiuSotdf7GmBnGmM3GmOU1tg8wxnxijPnMGPOrZAtBs2bw6KPe+ty58OKLSZ9ORETq\nF2+1zxPAgOgNxpjGwCNV27sBQ4wxXY0xVxpjHjTGtEuoJH37wrXXeutjxsD27QmdQkRE4hNX8LfW\nLga21th8KrDaWrvWWrsfeBYYZK19ylp7q7V2gzHmcGPMVOAHcf0y+MMf4Kij3PKGDfDrXyfwUURE\nJF4FKby2PVAWtb4e6BV9gLX2a+CGhk4UPTFB0ahRFN19t1uZMgWGDoXevVMopohI7istLfV10qu4\nG3yNMZ2Al6y1ParWBwMDrLUjqtaHAb2stWMTKkDNyVyshYED4eWX3XqXLrB0KTRtmshpRUTyWiYH\neZUDHaLWO+Ce/hNWbRpHY2DyZGjRwq1/8gnce28KxRQRyR9pn8axlif/AmAV0B/YALwLDLHWrkyo\nAMbYRYssZ57p4v5BDz8MN93klps0cU//3bolcmoRkbyVrq6es4C3gOONMWXGmGustRXAjcCrwArg\nuUQDf8TZZ0PPnjB1KnzzTdXG4mI47TS3vH8/jBih1A8iIj7JisRuUAIUAUW0bAnDh7vY3/XAx3DK\nKS74g6sOGj06g6UVEcmsSMPv+PHjU3ryz4rgf+ihlt27Y/f16wfFbWYxaM6VFHAAWraEFSvgmGPS\nX1ARkSySarVPVgT/bdssM2e6B/tVq2KPaV+wiVEVk7iex2g7qBe88EKNBgIRkXDJi5TODz44jpNO\nKmXlSliwAC6+GBpFlay84mjuYgIdWcfl8y5n8T2LyPA9S0QkI9Le2ycoMf38q5SVuXQ/06fD5s2x\nr+vR7QDFYxszbJjXK1REJCzyotqnvjLs2+cSfE5+qIIl78QOSG7VyjUQjx4NXbsGWVIRkeyR98E/\n2rIHFjD5tjU8zTB20zxmf79+Lh/cwIHe/DAiIvkoL+r8q43wrUfPX/Rn2qC/UE57JnITJzT5V7X9\nr78OgwdDp04wYQJs2hRMeUVEMiXv6/zrVF7u6nd27sQCC6+YwaTd1zB/fuwYsIICdzMYMwb69FEH\nIRHJH6Gq9jloyhQ3CgxchP/wQ8oOO5Fp01wD8ZYtsS/p0cO9RA3EIpIPwhn8KyvhrLPgzTfdeq9e\nbrlx44MNxJMmebujqYFYRPJBqOr8D2rUyPUDbdLErf/9726EGFBYCEOGwJIlLhfcyJFw6KHeS3fs\ncDnjunWD/v3djJEVFf59FhGRIIW3zj/auHEwfrxbbtHCpX7o0CHmsG3b4Mkn6xlB3B5GjXK5444+\nOrmiiIikUzirfSL27oWTT4aVVclEf/pTmD+/zpZda2HhQncTmDdPDcQikrvCHfzBVez36eOtP/cc\n/PznDb6srAw1EItIzlLwBxelp0xxy0ce6X4JHH54XC9VA7GI5CIFf4Dt210L7oYNbr1zZzjqqIRP\ns+yb7zN508U8veVcdlc2i9nf77B/MKbtXAYe8SYF5kDy5W3XDvr2dS3Oxx2n+iURSVheBP+SkhKK\nioooKipK/kQvvujSgfpgG4fxJFcxmWJW0SVmf3vWM4ppjGA6R1NL1rlEdOjgbgL9+7v8FO3apXY+\nEclreTWZi29lGD7cdevxiQUW0p/JFDOPQVTSuNr+AvYzmDmMYRJ9WIIvz+9du3o3gqIiaNPGj7OK\nSJ7Jiyd/38pw4AB89FHURMD+KdtcyLT5bZn+8tFs2VoYs7/HsbsovngDw368hRaH1jPX8IEDsGyZ\n63ZUWgo7d9Z9bKNGbhrLyC+D3r2rD1oQkdBS8E8zXxuIKyrg/ffdjWDhQnfCffvqPr6wEM44w/tl\n8B//4Q10E5FQUfDPoGXL3JiBp5+m7jmIi2HQoDhTTH/7rbsBRG4G//hH7GCEaC1bujQXF1wA114L\nTZsm/VlEJLco+GeBeEcQX389tG2b4IkXLfJuBitW1H3sj37kxjgce2zC5ReR3KPgn0UCH0G8caOb\ntCByM1i3rvr+Vq1gxgz3JiKS1/Ii+PvS1TPLBD6C2FpYs8Y97Y8fD/v3e/tuvBHuv1/VQCJ5SF09\nc0RDDcQtW7oG4uLiFEYQv/eeS2mxdq237ZRT4Pnn4fvfT/KkIpLN8uLJP9NlSBffG4ijbdvmGn1f\neMHb1qoVPPYYXHppSuUWkeyj4J+DAmsgttZNVnDbbdWrgYqL4YEH4JBDUi67iGQHBf8cFlgD8fvv\nu2qgzz/3tv3gB64aqHNn38ovIpmj4J8nfG8g3rbN/XSYM8fb1rKle4PLLvOt3CKSGQr+ecbXBmJr\n3c+KX/yi+sjhUaPgwQehWWzmUhHJDQr+ecy3BuIPPnDVQGvWeNt69nTVQMcf73u5RSR4Cv4h4EsD\n8fbtbpLi2bO9bS1auLqmK64IpNwiEhwF/xBJuYHYWpg6FW691c1/HDFiBEycqGogkRySavBv5Gdh\nkjVu3DhKS0szXYysZwyccw7MnevGc91xh5u1MqKiwg34PessV6szdWqN7NbGuHSjb7/tZhCLmD4d\nevWCTz5J10cRkSSVlpYybty4lM+jJ/8cl3QD8Y4dMHKku1tENG/u7hjDhgVebhFJjap95KCEG4it\nhUcfhZtvrl4NdO21blBY69ZpK7uIJEbBX2Ik3EC8bJnrDfTpp95BzZvD1VfDTTepR5BIFlLwlzol\n1EDccyfmhlEwa1bsiS64AG65xc0glnAeahEJgoK/xKWszNXwTJ8OmzfH7u/RA4pHW4Y1nU2LByfA\nxx/HHtS9u7sJDB2qnkEiGabgLwmJNBBPngxLlsTub9kShl9lKT7lHbq+8F/w8suxBx1xBNxwg2tA\naNcu+EKLSAwFf0laXA3EF21g0Cf3UTDzcdi1q/oBBQUuT9Att7hpJEUkbRT8JWXbt8PMmQ00EA//\nlusLZtJ25u/hiy9iD+rd290ELrooickIMsha+PZbOPTQTJdEJCEK/uKbuBqIL6lkzIl/o8+rv8W8\nWUu9UceOMHYsXHcdtGmTnoLHY/du+Owz16Np1arq/27bBkVF8Ne/QpMmmS6pSFwU/CUQcTUQX/AF\nwz6fQIu5T1afPAYy01X0wAE3qX3N4L5qlftADZk/Hy68MPhyivhAwV8CFVcD8c92Udz0cbr+eQL8\n+9+xB513nptLuGnT6n+FhbHbGtpXWOi6m37+eWyQX726+mC1RN1yi0t1LZIDsj74G2MGARcArYDH\nrbWv1div4J8jGmwgLjpAcddSBi2+jYKPl6a/gPVp3BiOPdb9CjnhBO/fsjK46ip3zEknuQ8pkgOy\nPvgffCNjWgP3W2uvr7FdwT/HNNxAbBn148+5vnw8bRc85RoT0uWoo6oH98i/3/ue+9VQ0zffuLaJ\nigq3vmULfPe76SuvSJLSFvyNMTNwT/BbrLU9orYPAP4INAYes9beV8fr7weettYurbFdwT9HWQuv\nv+6SytXZQHzeTsZ0X0Sf9p9j9u9z1TLRf/uS2FZRAR06xAb4zp3hsMMS/yB9+nhZ8Z57zqW6EMly\n6Qz+ZwLfAE9Ggr8xpjGwCjgHKAfeA4YAPwJOAf4AbAR+D/zVWruwlvMq+OeBuBqIE5mDOJ3uugsm\nTHDLo0a5zKYiWS6t1T7GmE7AS1HB/3SgxFo7oGr9dgBr7e+jXnMTcBXuxrDUWjutxjltSUnJwfWi\noiKKioqS+zSScXE1EMc7B3G6lJZC375uuXPn6gnuRLJEaWlptXlPxo8fn9Hg/zPgPGvtiKr1YUAv\na+3YBM6pJ/885dscxEHbu9elr96zx62vW+eqlUSyWKZn8lLUljr17OmmCN6wwc0SecIJ1fe//jr8\n7GfQqZOrddm4MSPFdF1I+/SpXjCRPJdq8C8Hoh+ROgDrEz2JpnHMb4cd5sZ6rVwJCxbAxRdDo6hv\nXnm5q3bv2BEuvxwWL05vByHA/QyJWBjTNCWSNTIyjWMt1T4FuAbf/sAG4F1giLV2ZQLnVLVPCGVd\nA/G777p5jMElMyor09wFktXS2dtnFnA2cASwBbjLWvuEMeYneF09H7fW3ptQART8Qy1rGogrKlyq\n6h073PqqVZrBTLJa2ur8rbVDrLXtrLVNrbUdrLVPVG3/i7X2BGvtcYkG/ghV+4RXYSEMGeKqepYu\ndXPKRyfY3LkTHnkEunVzE4nNmeONx/JVQQGcfba3rqofyVIZqfYJgp78paa4UkxHz0Hsl4kTXX4f\ncC3Rs2f7eHIRf+VMeoc6C6DgL3WIawRxZA7iPj5U0S9f7vL7gKsC2rKlesu0SBbJdFdPX6jaR2pj\njKvqmTsX1q6FO+90qXsiKipcNoazznLdSqdOdal6knbiiV5en6++go8+SqX4IoFQtY+EUuANxJdf\n7u4oAA88AL/4RUrlFQlKXjz5i8Qr8Abi6P7+GuwleUxP/pLzfG0gXrMGjjvOLbdoAV9/rakdJSvl\nxZO/6vwlFb6OID72WHcguAaE998PvPwiiVCdv0g9UhpBfM018Kc/ueV77oE77gi6uCIJy4snfxG/\ndejgksWtWwfPPFM9bxu4Xp2jR0O7djB2rPvVcFD//t6yBntJntKTv4RG3Cmmf1ROQadj3MamTWHr\nVmjWLL2FFWlAXjz5q85f0iHuFNO92zPhiD+ykaNdrv+3385MgUVqoTp/kRQ1OIKY/QxmDmOG7aDP\nkyOV5FOyitI7iPgg61JMJ8pad/eyNr7l5s3VhTXHKfiL+OjgCOKH9rPkndjgmPY5iLdsgYEDXaqJ\n2oJ4zZ8r8WrVCp59Fn7yE3/LK2mj4C8SkGVdLmPyqn48zTB20zxmf9++LqlcoHMQ33wzPPRQMOc+\n6STXCi45ScFfJCi33QYPPMB2WjGz6E9M3nhxnSOIR46EESN8TjH97beuL+q2bfEdb4z316hR3eu7\ndnmj3JYudS3hknPU20ckKFX9/Q9jBzd9Pa7eEcQlJd4I4r/9zac5iP/8Zy/wH3us65+6Z4/rgbR/\nv0tadOBA9SqgAwfc9n373HF79ribyK5dbsTyzp2ukBEzZ/pQUEkn9fYRCdrOnXD44V5muC1bDqZ8\nbqiB+MQTvQbili2TfP8zz/RSl957L9x+e5InquH//s+r6z/ySHf3CqzeSoKiah+RIPXuDW+95Zaf\nfx4uvbTa7sBSTK9YAd27u+WCAne3Ofro5D5DTRUV7mfKxo1u/eWX4YIL/Dm3pE1eVPuIZK0GUj3U\nTDE9alTdKab79UsgxfRjj3nLgwb5F/jB3UyGDvXWn3zSv3NLztCTv0h9Sktdtx6Azp3h008bfEk8\nKabrbSDes8cd9PXXbv3VV+Hcc5P+CLX6+GM3eAFcCouNG6FNG3/fQwKlah+RIO3Z44Linj1ufd06\nlzUuDvHOQVxc7Kr3D44gfuYZ78m8Uyc3x0AQcwmfcgp8+KFbnjbN3ZEkZ6jaRyRIhxxSPSVoArN7\nxTsH8dlne3MQ79yJa0mOGDEiuEnkhw/3ltXrJ3SyIvirq6dkNR+mdow3xXT7tgcYu2gwK+kCjRu7\nuQWCMmSI18vnrbdg9erg3kt8o66eIuny7rvQq5dbPuYYF8F9yPLWUIrpvt9ZzpipPYIdQTxwILz0\nklv+7W/h7rsDeiPxm+r8RYJWUQFHHAE7drj1Vavg+ON9O/3BBuJJlaz6NPbHeGAjiMENJIt0Xw2y\nfUF8pzp/kaAVFLiK+Ygkq37qcnAO4nHPs4D+XMxcGnHg4P7ARhADXHih18tn7VrXZ1VCQcFfJB5p\nmNrRTH+U/rzOXAaz9v89ElcD8TffpPimTZtWT/egPv+hoWofkXgsX+6yYIKrAtqyxd/qkU8/9aYW\na9TItSu0bx/cCOJo77wDp5/unXDTpuoj1SQrqdpHJB26dz+Y14evvnI3Az9Fj+j96U9dRT+xI4hH\njqx7BHH//gmMII7Wq5fXhrFzJ7z4YmqfRXKCgr9IPBo1qt7l08+qn7174YknvPU6BlvFPQdxJ9et\ndNOmON/fGLjqKm9dff5DQcFfJF4+9Pev1bx58O9/u+VjjoEBA+o9/GADcT0ppu+6y40tuPxy96uh\nwZrVK6/0lhcscCeRvJYVwV+DvCQnRAf/RYtcTn0/RI/ovf56N7grDvGOID7rLO9XQ50NxB07ejmM\nKivhf/4nqY8SKpEbdpppkJdIulnr6lTWrXPrb73lNZQma/VqlzAO3OP72rVx5w6qzb597mYwaVLt\nDcStWnkNxF261Ng5cyZcfbVb7tbNJX/zYTBb3ti82d30S0vd38qVbtuRR2akOGrwFUkXY/yv+olu\n6D3//JQCP7gG4khVT20ppnfsgIcfdj2DIr8aDjYQX3KJd/CKFfDBBymVJedt2QKzZ7uJmrt3d2m1\nL7sMpkxxgR/czSBHKfiLJCK6v3+qwX/fvrgaepMVGQtQXwPx4MFRDcS7WroNEWHr818z2B91FPz8\n566f7YoVsccXFnq/AnOQqn1EElFe7hplwQ2Q2roVmjVL7lzRqRXat3dVPgFOp9hQiukmTWBw700U\nl15KH5ZgvvMd93kLCwMrU0Z9+WX1apx//rP+45s0gdNOg6Ii93faaRkdD6HcPiLp1qWLN0vLwoXV\nq4ISce658Nprbvmuu2D8eH/KF4eG5iDuwUeMYRJDnx1Ii8vyZIrHL790uTHeeCP+YN+rlxfsTz89\nqwa/KfiLpNuYMa4qAOCOO+CeexI/x7/+Bd//vls2xj31d+zoWxHj1WADccEuho9uXnsDcTaprHTB\nfcMG91deXv3fzz/36unr0qQJnHqq6/WUhcG+JgV/kXSbM8eNpgL30//ttxM/x29+A/fe65bPPx9e\necW/8iVp2TLXlvnUk5Xs/ja2ObBfP3ffGzgw0Nqp6qx1rdTRgby25Y0bEx/aXFAQ+2TfvHkQnyIQ\nCv4i6fbVVy7Vg7WuT/7XX7s+lPHav9/16onUt7z4opukPUts3w4zT7qfyet+yipiH/fbt3e9iEaM\n8Hde+YPmz4eHHnKNqRs2wK5d/py3oMA92UeC/Rln5FSwr0nBXyQToue/fekll48nXnPner1q2rZ1\nQS5tj9JxmjQJe+ONvE4/JrW5k3nb+9beQFw1B3GfPj4NCfjHP9wTeDID6Nq0cXemdu28f6OXu3aF\nFi18KGR2UPAXyYTbboMHHnDLt94K//3f8b92wAB49VW3fOedrp9ltvnqK3djqgrCZa9/xrSFxzF9\nuusRWVOPHu4mMGxYCvF1xw53U12zpvr2Zs1cAK8toEeW27ZNvtdVjko1+GOtzeifK4JIjnnlFWtd\nxY+1PXvG/7rPP7fWGPc6Y9x6trr4Yu8z/vrX1lpr9+619plnrO3d29sV/deqlbVjx1q7cmWC71VZ\nae2QId6JWra0dskSa7dudfskRlXsTDr2BjrIyxjTxRgzxRjzvDHmuiDfSyStzjzTq6pZtsz1NInH\n4497WdbOPdeNsMpWw4d7y089BZWVB1NML1lSe4rpekcQ12fGDJg1y1ufNg1694bWrZViIiip3Dni\n/cONJH6+jn1B3BRFgnfGGd6T6vPPN3z8/v3Wtm3rvWbOnODLmIq9e6094givvAsW1HrY1q3WTpxo\n7Qkn1P5roH17a+++29qNG+t4n48/trZZM+8F110X3GfKI6Tjyd8YM8MYs9kYs7zG9gHGmE+MMZ8Z\nY35Vx2svBF4Bnk3lJiWSdRJN9fDKK65LIrjUARdeGEy5/FJYCFdc4a3Xke6hdWsvxfRrr8FFF9We\nYrpjR29imoPNfLt3uxQK337r1rt1cz19JHjx3CGAM4GTgeVR2xoDq4FOQBNgKdAVuBJ4EGhX4xzz\n6jh3wPdHkYC88Yb3tNq5c8PHn39+TB161nvvPa/MzZtbu3NnXC9bt87aO+6w9sgja/810KOHtVOn\nWrtz+BhvY7Nm7leAxIUUn/zj7u1jjOkEvGSt7VG1fjpQYq0dULV+e1Uk/33Ua84GLgEOAVZaa/9Y\ny3ltvGUQySp79rjuhXv2uPV16+rOyrlunavfj3zX16yBY49NSzFTYi2ceKKX2GzmzOqzfjUgMgfx\npEnw5pux+1uxneHMpJjJdJl+m5vPQOKSam+fVDoXtwfKotbXA72iD7DWLgIazHkaPTFBUVERRUVF\nKRRLJE0OOcR1cF+wwK2/8UbdgTG6offHP86NwA/eFI+33+7WEwz+kQbiIUNcu/jkyfD00662B2AH\nh/EwN/EwN9FvlmXM4WkeQZxDSktLfZ30KpUn/8HAAGvtiKr1YUAva+3YhAqgJ3/JZffe61I1gAuK\ntc1/W1HhnvojUyPOnu2lh8gF5eXuF421vuQh2rZ5LzN/9DCT11/Ip5wQsz/wEcR5IpOTuZQD0b9x\nO+Ce/kXCo+bkLrU9yPzlL17gP/JI92ibS9q3h3POccvWukf3FLS+91fcvP6XfEIXXms8gIuKtsbf\nQCy+SSX4vw90NsZ0MsYUApcB85M5kebwlZz1wx96eX3Wr3fTMtYUPUfvNdfkZn786D7/Tz6ZfDSe\nP9/NLAMY4JwHfsILb7Rh7VqXIDV6RsT9++HZZ+OcgzhE/JrDN97ePrOADcBeXD3/NVXbfwKswvX6\n+XUyLc6ot4/kugsv9HqsTJlSfd+6ddY2auTt/+yzzJQxVd98Y22LFt7neOedxM/xxRfWtmnjnWPg\nwJjRu4GNIM5DpKOfv7V2iLW2nbW2qbW2g7X2iartf7HWnmCtPc5ae2/qtyKRHFTfvL4zZnhTZvXv\nD8cdl75y+al5c2/WMUh8iseKCjdmYOtWt96hg5vCssbo3cBGEEuMrJjDV9U+ktNqDvaKBPsDB6pP\n0O7zHL1pF93LZ9Ys2Ls3/teWlHh9PRs3dq8//PB6XxKp6ikvhz/+EY4/vvr+mDmIN8VfnFzmV7WP\nsnqKpKqy0nVLieT3WbrURa5XXvFSPX/3u65NIBfr+yMqK10X1S++cOtz5sAllzT8utdeg/PO89oJ\nfvc7r4dUAqx1s2ZOmuSaDmqmmC4ocDeDMWN8TDGdxTLZ20dEwOUyqK3qJ7qh9+qrczvwg/ucV17p\nrdfWrbWmTZtcnudI4D/nHG/MQIKMcS9/4QVqbSCuqIDnnvMaiKdOVQNxfbIi+KvaR3JedPBfuNDV\nVbz8srctX0auRlf9/O//1p/NtLLS3SwiEwAcdZTLDtoo9bDToYObOrmsDJ55xiUAjbZ8OYwe7Xqp\nRvIO5Yu09vYJ8g/19pF88NlnXpeUFi2s/e1vvfW+fTNdOn+dfrr32SZOrPu43/3OO84Ya197LdBi\nLV1q7ciR1h56aO09hfr1s/bPf3bJVfMB6crtExTV+UtesNa1PK5b59YPOcTL+TNrFlx+ecaK5rtp\n0+CGG9zyD38I778fe8ySJW6e3AMH3PpvfuPq+tNg2zZXIzV5Mnz6aez+fBlBrGkcRbLFNdfAn/5U\nfdsRR7gqoKZNM1KkQGzd6qZNjPT2+fhj6N7d2//VV/CDH7gGbnB1MqWlaU/YE2kgnjwZ5s3Lvwbi\nvGjwVZ2/5IXoev+Iq6/Or8APLpNpdIqK6D7/1sK113qB//DD3S+fDGRqizQQz52bXw3E6uopkm3K\ny+GYY6pvW7kSunTJTHmCFN2NtV07V93VuLFL3XDLLd5x8+ZlVS6jBlNMt3KZLEaPdgPJspmqfUSy\nSZcusGr66TNdAAAF8klEQVSVWz7rLFjUYEbz3LR/v7vRRXryvPqqe8o/4wy3D+Dmm93orCxVW4rp\naP36uSqhbE0xnRfVPiJ5I3pqxuLizJUjaE2awNCh3vqkSa5ROxL4TzkF7rsvM2WLU/QI4okT4YQa\n2aXzfQRxVjz5l5SUaBIXyQ87drig166dC/651oqYiKVL4eSTY7e3bAkffJBzeYysdQF/0qTaG4ib\nNHE3g+LizDYQRyZ1GT9+vKp9RCRDevaEjz6qvi0PuraWlbkB2tOnw+bNsft79HBVQkOHQosW6S8f\nqNpHRDKp5pSO11+f84Ef3AjiCRNcO/asWe5JP9ry5W6ow623ZqZ8ftCTv4gkb/NmV1m+fbvr6//u\nu9VzMOeRZctgyhSXoSLSQPz3v8Opp2amPOrtIyKZ9eGHbkTvlVdC69aZLk3gtm93I4jfftv9KsiU\nvAj+avAVEYmPGnxFREJMDb4iIpIwBX8RkRBS8BcRCSEFfxGREFLwFxEJoawI/srnLyISH+XzFxEJ\nMXX1FBGRhCn4i4iEkIK/iEgIKfiLiISQgr+ISAgp+IuIhJCCv4hICCn4i4iEUFYEf43wFRGJj0b4\nioiEmEb4iohIwhT8RURCSMFfRCSEFPxFREJIwV9EJIQU/EVEQkjBX0QkhBT8RURCSMFfRCSEFPxF\nREIo8OBvjGlujHnPGHNB0O8lIiLxSceT/38Cz6XhfXKektt5dC08uhYeXQv/xBX8jTEzjDGbjTHL\na2wfYIz5xBjzmTHmV7W87sfACuBLf4qb3/TF9uhaeHQtPLoW/on3yf8JYED0BmNMY+CRqu3dgCHG\nmK7GmCuNMQ8aY9oBZwOnAVcAI4wxSWegS0QiX5B4jq3rmHi317ce9JdZ16Lu90712ESuRTzbdC1q\nXw/yWiR67ny6FnEFf2vtYmBrjc2nAquttWuttfuBZ4FB1tqnrLW3Wms3WGvvtNbeCjwDPJqu3M0K\neHW/d6rH6lo0fEy2/U9eG12L5M6dT9ci7nz+xphOwEvW2h5V6z8DzrPWjqhaHwb0staOTagAxiiZ\nv4hIElLJ51+Qyvum8FrvJCkUXkREkpNKb59yoEPUegdgfWrFERGRdEgl+L8PdDbGdDLGFAKXAfP9\nKZaIiAQp3q6es4C3gOONMWXGmGustRXAjcCruO6cz1lrVwZXVBER8UvGJ3AXEZH0U24fEZEQyrrg\nb4wZZIx51BjzbNUI4dAyxnQxxkwxxjxvjLku0+XJNOWJcowxRcaYxVXfjbMzXZ5MMs7vjDEPGWOu\nynR5MskY06fqOzHdGPNmQ8en0tUzENbaecA8Y0xr4H7gtQwXKWOstZ8Ao40xjXCD6B7PcJEyTXmi\nnEpgJ9AU9bC7CGgP/JuQXwtr7RJgiTFmEPBuQ8en5ck/ydxAd+LSR+SVRK+FMeZC4BVc8M8riVyL\nfM8TleD3YrG19nzgdmB82gsbsASvxfHAm9ba24DRaS9swJKMnVfgsirUz1ob+B9wJnAysDxqW2Ng\nNdAJaAIsBboCBrgP6J+OsqX7L5FrUeN18zJd9gx/L+4BHsT1LnuRqs4K+fKXzPcCKARmZ7rsGf5e\nDAUurTrmuUyXPdPfC6AjLpVOg+dOS7WPtXZxVXqIaAdzAwEYY54FBgHnAP2BVsaY46y109JRxnRJ\n5FoYY44ELgEOAd5IYzHTIpFrYa29s2p9OPClrfqm54sEvxddgPOA1sDDaSxmWiQYLyYCDxtjzgRK\n01fK9EjwWqwErgVmxHPuTNb5twfKotbX4+UGyrsvdAPquhaLgEWZKVLG1HotIivW2plpL1Hm1PW9\n+D3wQmaKlDF1XYtvgeszU6SMqfP/EWvtuHhPksnePnn15JYiXQuProVH18Kja+Hx5VpkMvgrN5BH\n18Kja+HRtfDoWnh8uRaZDP7KDeTRtfDoWnh0LTy6Fh5frkW6unoqN1AVXQuProVH18Kja+EJ8loo\nt4+ISAhlXXoHEREJnoK/iEgIKfiLiISQgr+ISAgp+IuIhJCCv4hICCn4i4iEkIK/iEgI/X/A4SAI\njO26MgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108623310>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot\n",
    "\n",
    "# Hit and miss Monte Carlo integration\n",
    "ngroups = 16\n",
    "\n",
    "I = np.zeros(ngroups)\n",
    "N = np.zeros(ngroups)\n",
    "E = np.zeros(ngroups)\n",
    "\n",
    "n0 = 100\n",
    "for i in range(ngroups):\n",
    "\n",
    "    N[i] = n0\n",
    "    x = np.random.random(n0)\n",
    "    y = np.random.random(n0)\n",
    "    I[i] = 0.\n",
    "    Nin = 0\n",
    "    for j in range(n0):\n",
    "        if(y[j] < np.sqrt(1-x[j]**2)):\n",
    "            Nin += 1\n",
    "            \n",
    "    I[i] = 4.*float(Nin)/float(n0)\n",
    "    E[i] = abs(I[i]-np.pi)\n",
    "    print n0,Nin,I[i],E[i]\n",
    "    n0 *= 2\n",
    "    \n",
    "            \n",
    "pyplot.plot(N,E,ls='-',c='red',lw=3);\n",
    "pyplot.plot(N,0.8/np.sqrt(N),ls='-',c='blue',lw=3);\n",
    "\n",
    "pyplot.xscale('log')\n",
    "pyplot.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 3.20920994888 0.0676172952856\n",
      "200 3.19023195241 0.0486392988246\n",
      "400 3.13539345091 0.00619920268006\n",
      "800 3.18467753241 0.0430848788188\n",
      "1600 3.1148121176 0.0267805359933\n",
      "3200 3.14792033551 0.00632768192034\n",
      "6400 3.14749188649 0.00589923289621\n",
      "12800 3.13735813292 0.00423452066773\n",
      "25600 3.13359006584 0.00800258774604\n",
      "51200 3.1398924647 0.00170018888973\n",
      "102400 3.1435290202 0.00193636660599\n",
      "204800 3.14206545238 0.00047279879124\n",
      "409600 3.13999805065 0.00159460293962\n",
      "819200 3.14154966288 4.29907095283e-05\n",
      "1638400 3.14063731607 0.000955337516072\n",
      "3276800 3.14223148552 0.000638831930158\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEHCAYAAABGNUbLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XecFdX9//HXoYMiIoJGBEGkGruAJj8jsUQQ80CjQTEU\nJREV7AUTU3aN5iuoUWMoKqGpUdEkRlGxs2pQBEWUWFmKihgp0pe2u5/fH2eXudtvmVv23vfz8dgH\nM3fmnjk7rp8793OaMzNERCS3NEh3BUREJPUU/EVEcpCCv4hIDlLwFxHJQQr+IiI5SMFfRCQHKfiL\niOQgBX8RkRyU1ODvnOvsnPubc+7JZF5HRERik9Tgb2YrzOxXybyGiIjELqrg75yb5pz71jm3pNLr\n/Z1znzrnljrnbkpOFUVEJGzRPvlPB/pHvuCcawhMKHu9FzDEOdfTOTfMOXePc+6gcKsqIiJhiSr4\nm9mbwIZKL/cBCs1spZntBh4HBpnZw2Z2rZmtds7t55y7Hzha3wxERDJHowTe2x74KmJ/FdA38gQz\n+w64rLZCnHOaVlREJA5m5uJ9byINvqEFbTML9ScvLy/Uc2s6J9rXa9uvaVv3IrvuRTSv6V6k/l7E\nWl4m3YtEJRL8vwY6ROx3wD/9xyw/P5+CgoIEqlJRv379Qj23pnOifb22/VjqGg/di5qvnei5sdyL\naF7Tvah+P5n3ItayM+FeFBQUkJ+fX2c96hTtJx7QCVgSsd8IWFb2ehNgMdAz1k9ewNatMzGzvLy8\ndFchY+heBHQvAroXAR++4//WE21Xz8eAt4BuzrmvnHMXm1kxcAXwIvAxMMvMPonnA+jgg+Hii+Hd\nd+N5d/ZI9tNefaJ7EdC9COhehMdZCLmjhCrgnEEe0A/oR+/eMGYMnH8+NGuW1qqJiGScgoICCgoK\nuOWWW7AEGnwzJPhXrUObNjByJFx+OXTunIaKiYhkMOdcQsE/IyZ2m1+wg2HDoEmT4LX16+HOO6FL\nFxg4EJ5/HkpL01dHEZFskhHBf84Z+zJyyzmseuG/jBsHhxwSHDPzgX/gQOja1X8grF+fvrqKiKRT\nWL19MiLtU6EGffpQ8stRzGl9IROnNeeFF6q+p1kzuOACGD0aevdOVU1FRDJHommfzAv+5fbeGy68\nkML+V3D/f77PtOmODZUnmIA9DcSDB0Pz5smurYhIZsiO4P/aa/Dgg/Cvf8GuXVVPOuYYikZczqxG\nv2Di9Ba8917VU/bbD375SzUQi0huyIoG3/zXX6fg0kvh66/hz3+GHj0qnvD++7S4ZhQXj23LwiNG\nMn/KEoYPN5o2DU757js1EItI9suunH/lOpjBvHn+28CTT8KOHVXfePjhrLvwKqbuHMrkmS344ouq\np3Tu7L8JjBzpu46KiGSL7Ej71FaHDRvg73/3HwRLllQ93rQpJecOZs4RY5n4+uG88ELVe6EGYhHJ\nNtkf/MuZwYIF/kPg8cehqKjqOd27s+ycG5i8dSjT/t5MDcQikrVyJ/hH2rwZHn0UpkyBRYuqHu/S\nhaIF/2XW082YOBE1EItI1smOBt9Yp3TeZx+47DIf1d99Fy69FFq2DI4vW0aLOf/k4oth4UKYPx+G\nD0cNxCJS72V3g288tm6FsWNh8mS/f9JJ8MYbFU5Ztw6mTvWnVNdAfOih/jMlqQ3ECxb4Rogjj0zS\nBUQkF+Rm2qcm33wDHTtCcbHf/+9/4fDDq5xWUgJz5sDEiaR2BPHf/w5Dh/rtc84JvnqIiMQoK9I+\nofne92DQoGD/wQerPa1hQzjrLP8BUFgI118PrVsHx3fsgBkzoE8f/zNzJmzfnmDdSkvh1luD/aee\ngl69/LeVzZsTLFxEJDbZ9eQP8PLL8JOf+O199/UDx1q0qPNtRUUwaxbJayB+6SU444zqj7VrB7fd\n5vNNDRvGUbiI5JqsePIPdQ3fU08NUikbN8ITT0T1thYtSG4D8X33Bdtnnw0nnhjsr1kDo0bBccdB\niGsZi0j2UYNvbe64A266yW+fcAK8/XZcxdTVQBz1COKlS6FbN7/tHHz+uf8Uefxxn/ZZVWnde7UH\niEgd1OBbnTVr/MLAu3f7/cWL4aij4i4u4Qbia66Bv/zFbw8cCM8+GxwrKoK77oLx4ysOXGvSxL/v\nt7/1XVtFRCIkGvzjXvk9rB9fhSS44AIzPy7Y7PLLQyu2sNDs+uvNWrcOio/86d3bbPp0s6Kisjds\n3mzWsmVwwosvVl/wV1+ZDR1atcB27cwefNCsuDi030FE6r+y2Bl37M3OJ3/wufMf/9hvt2wJq1f7\nNQJCEm0D8WUtHuLQW0b4F7t3h48/hga1NLW8845/4p8/v+LrRx0F994L/fqF9juISP2ltE9NzKBn\nT/jsM7//4INwySVJucyCBTBpkv8w2Lmz4nFHKQOYwxgm0v++gTS4ckx0hao9QERqoeBfm3vugeuu\n89vHHeengkiiOhuIO5Vy+egG0Y8gVnuAiNQgK3L+eXl5Nnfu3FDyYBWsX2/WtGmQP1+4MPxrVKO4\n2Gz2bLP+bRdW2y7QtKnZiBFmCxZEWaDaA0SkzNy5cy0vL085/zoNHw4PP+y3f/UrPxNoKpR171zG\noUxmNNNaXcOGTVUHcPXu7XsJnX9+FFNM19QecMwx8MwzvoeTiOSErBjklVSXXhpsP/ZY6qZSmDAB\ngC4s566zCli1uiHTpvnsU6SFC/3gsoMP9in+5ctrKbNvX3jrLT9HUGSgf/9937qc5g9yEak/sv/J\n3wyOOAI++sjvT5rkR2Yl05Yt0L69/xfgxRf3TDlRZwOxgwED/IIz/fvX0jGoqMg3/P7xj8FQ49mz\n/aRFIpL11OAbjQkT4Mor/faRR/pBXy7+dpKYrtejh+/eWc31QhlBfPnlcP/9fvuww/yHXJMm4fwe\nIpKxFPyjsXEjHHRQMDXn22/7aR+SobTUdzH9/HO/P3GiT+rXoq4RxE2b+hHEY8ZUM4J47Vo/dcTG\njX7/zjvhhhsS/z1EJKMp5x+NffeFIUOC/QceSN61XnopCPz77OMbnOtQ1xTTO3f6aaXLp5ieMSNi\nium2bSEvLzj51lvh229D+3VEJDvlxpM/+ER7375+u1kzP+I3MsKG5cwzfRQHuPZauPvuuIqJegTx\nZXBoh90+nfXpp/5gKns1iUhaKO0TLTM49lif7wc/0dpVV4V7jc8/91M4gM/xL12a8EjcqBuI+75L\n/7w+NMD8i++957uAikhWUtonWs75x+RyDzwQftfIsu6dgJ+9M4QpGJzzX1hmzvQzPYwbB4ccEhw3\n82sLDMw7nsOar+ZObmC9tYarr1bXTxGpUUY8+efl5dGvXz/6JXvSsi1bfMPv1q1+/403/ELvYdi8\n2XfvLC/7pZfg9NPDKbuSOhuI2cEFPM6YP7Wn983JqYOIpEdBQQEFBQXccsstSvvE5NJLg7V9f/EL\neOSRcMr961+DNFLPnr7LZTK7k5YpLPQ9PadNgw0bqh7vfVwpo69oEN0IYhGpN5Tzj9X77/vcP/j+\n8F9/Dfvvn1iZpaW+P//SpX4/FQPJKtnTQHxfMe8tblTleIUG4kNTWjURSQLl/GN1zDFBZ/ldu3wy\nPVEvvhgE/latYNiwxMuM0Z41iBc14p2b/sVwZtKUHXuOl69BfNhhca5BLCJZJfeCP4Tf8Bu5OPvI\nkaEuGhMr56DPnwYx88i7WcXBjGcsnfZas+f4ngbigf6D4M47Yf36tFVXRNIk99I+ANu2+Ybf8kne\nXn0VTjklvrI++8ynfCC07p2hiFjJrIQGzLnzv0x6reeeIQiRah1BLCIZSWmfeOy1V8XUTCIjfiO7\nd551VmYEfvDLPZ57LgANKeWsJ0bw/LOlsY8gFpGslJtP/gBLlvhRsQCNG8NXX8EBB8RWRuXunS+/\nDKedFm49E7Fihe95VD4ybMYMGOHXE96+3a8UGdUIYjUQi2QcPfnH64gj4Ac/8Nu7d8P06bGXMWNG\nEPh79oRTTw2teqHo3Nk/5pf79a/3TDPdvHlZA/FCv0bM8OE+/VNODcQi2S13gz9UXOhlypTYoltp\nqe/bX+7KK1PSrz9mv/kNfO97fvt//4Pbb69w2Dmf6ikfQTx+PHTqFBxXA7FIdsrdtA/43Ef79sHo\nqIhFV+pUHhHBd+9ctSqtvXxq9fDDweyiTZrAJ5/UmsspH0E8aZIfQVz5P48aiEXST2mfRDRvvicH\nDgSLokQjsnvnL3+ZuYEf/EjmPn389q5dcOONtZ5ePsX088/7zktqIBbJPkl/8nfODQIGAvsAU83s\n5UrH0/fkD34a5J49/XbDhvDll74baG0qd+8sLMz8VtH58+HEE4P9GLu3qoFYJLNk/JO/mT1tZqOA\ny4Dzk329mPXoASef7LdLSvwkOXWJ7N7505/Wj2h3wgkwdGiwf801UFwc9dvVQCySXaJ+8nfOTcM/\nwa8xsyMiXu8P3As0BP5mZuNreP9dwCNmtrjS6+l98gd47DG48EK/3aGD7yLZsGH1527aBAcfHPTy\neeWVzOvlU5Ovv/ZLPhYV+f3JkyuOdo7RunX+s3LyZFi5surxqNYgFpG4JPrkj5lF9QOcBBwDLIl4\nrSFQCHQCGgOLgZ7AMOAe4CDAAeOBU2so19Juxw6z/fc3822bZs8+W/O5994bnNerl1lpaerqGYbb\nbgvq36aN2XffJVxkcbHZ7NlmAwaYORcUX/7TtKnZiBFmCxYkXn0R8cpiZ9QxvPJP1GkfM3sTqDxp\ncB+g0MxWmtlu4HFgkJk9bGbXmtlq4ErgVOA859ylZKKmTX1Oo1xNI37rS/fO2lx3XbAazPr18Mc/\nJlykGohF6p+YGnydc52A2VaW9nHOnQecYWaXlO0PBfqa2ZUxlGl5EQuQp2RRl+oUFkLXrn67QQOf\nx+jQoeI5zz3noxz47p1ff+2niqhv/vEP+PnP/XajRvDhh0Gjd0jUQCwSrvJFXMoluphLTF8T8Omd\nyLTPucCUiP2hwF9jLDPUr0IJOe20IFfxhz9UPf6TnwTHr78+9fULS2mp2Y9+FPwu/fsn9VLvvGM2\nfLhP/1ROCTlnduaZZs89Z1ZSkrRqiGQdUpX2qcHXQOTjcQdgVayF5OfnV/hES5vIEb9/+1vF3jCf\nfuqXZgSf6hkzJrV1C5NzfgH78pTVCy/4nE2SLqURxCLhKSgoID8/P/GCYvmkoOqTfyNgWdnrTShr\n8I2xzOR9NMZq1y6zAw8MHkufeio4NmZM8PqgQemrY5hGjQp+p27dzHbuTMllM6qB+JtvzMaPN3vt\ntRRcTCQ8JPjkH0uQfgxYDewEvgIuLnt9APAZvtfPb2KuQCYFfzOzm2+umg7ZuNFsr72C1195Jb11\nDMu335rts0/we919d8qrUFjoM2itW1f9EACz3r3Npk83KypKwsVffdWsXTt/oYYNzRYtSsJFRJIj\nZcE/WT+A5eXl2dy5c0O/OXFZsSJ4HHXObPlys3vuCaJRfezeWZs//zn43Vq18h8IaVBUZDZtmtlx\nx1X/IbDffmY33mi2bFkIFysp8V1eGzSoeJEzzgihcJHkmjt3ruXl5WVH8M84AwYEAeGmm8wOPTTY\nv//+dNcuXDt3+pRP+e83alRaq5P0BuJ16yr+96388+qrof9OIsmQaPDP7Vk9a/L003D22X67UaOg\n4XfffX2rZX3s3lmbyC6szsGiRXD00emtE0kYQbxgge/i+uWXwWsnnQTt2sE//+n3jz/en1ffxm9I\nzknZCN9k/ZBpaR8zs927zdq3r/pUeMMN6a5ZcpSW+vaN8t/z5JMzKrWVcANxaanZhAlmjRtXfOPY\nsf6/9ZdfVvya8cQTqfz1RGISVtpHT/41yc+HW24J9p2DZcv842Y2+uQTv7pZSYnf79vXz27arh20\nbev/jdxu29Y/bjdqlNJqLlvmvwlMmxYswxCpd28YPRrOP99PRseWLTBqlB9xVq5VK9/3dNCg4LWx\nY30/U/CD/T76yC/vKZKhEn3yV/CvyapVfhqE8qkpBw2Cf/87vXVKtmuu8f3/o+WcH6pb3QdD5L8H\nHeQ78YeYSolqBPGgdVz2+hAOXf5KcOCYY/wI58rDir/7Drp0gY0b/X6Ck96JJJuCfzL9/Oc+UAAU\nFARTP2erDRv8nP+ffRZ+2d27+yfw4cNh//1DK9bMTzM9cSLMmhWsVV/OUcoA5jCGifS/pCMN7rsX\nmjWrvrDx4/06xwAHHuin/Mi29h3JGlkR/PPy8tI3p09ttmyBu+6CI4+Ec89Nd21So7jYB/+1a2HN\nGv9Tvl353+++i738Jk38vRw1yn+YhvhtYN06mPbAbiaP38TKLVU/YOpsIC4q8imf1av9/p/+BDff\nHFr9RMJQPsdPonP7ZETwT3cdJE67d/u5GKr7YIjc/vDDYP2DSN26+Q+BESPC+TawYgWcdx4lixYz\nhwFMYjQv0J/Ks5jUugbx3/4Gl1zit/fZB5Yv12IEkpGyorePZLktW8ymTPHDdavrW9+kidkFF/gp\nFuLtZfT002b77lux3MGDrXDxlthGEO/ebdajR3DCddeFdRdEQoV6+0i98v77MGUKPPKIT6tV1rWr\nf/IeMcI3FteluBh+9zufry/XuDHcfbd/tC9LK8U0xfQHT8HPfuYPNGkCn38erIEgkiGy4sk/4/r5\nS/Jt3Wo2dapZ377VP5I3bmw2eLCfR6mmobyrV1ecmhrMOnQwmz+/xstGN4K41J7rfq2VUDaoYMSI\n5NwDkTion79kjw8+8N8GHn4YNm+uerxLF/9t4KKL4IAD/GsFBT5x/+23wXn9+/tvFFHm6OscQcxy\nLmcyI5lOmw/n+nEQIhkiK3r7pLsOkiG2bYMnn/TLaM6fX/V4o0Z+2o0uXfyArPIxGM755Shvvtmv\nwhajkhKYMwcmTfJLG1T+c2zGdi44eB6j/3Va1QZikTRR8Jfs9OGHwbeBTZtqPq9tW3j0UTjttFAu\nG80I4jFjYPDgshHEImmi4C/ZrajIfxt48EF4662Kx374Qz+yq3370C+7p4H4hhW8913VKT3KG4gv\nvzx7Z/yQzKbgL7njv//13wYKCvx0G7//fdLn37EVK1nY7RdMLB7FLM5nJxVHBzsHAwb4bwP9+8eV\ndRKJS1YE/4wd4SsCcO21cO+9rKMN09r9msktrmflyqr/zx16qP8mcPHFGhcmyaMRviKpsm6dj+xl\n4xJKHpzKnO+NrLmBuJnviDR6dDUjiEVCkhVP/umug0idbrvNp5nAtzEsXQrNm0fdQHz++TXPJycS\nDwV/kVTYts13MS0fVzB+vF8DoExdI4jbtPETyqmBWMKi4C+SKpMn+1wO+CU9ly+H1q0rnGLmV4Gc\nNKmGKabVQCwhUfAXSZXdu6FXLz/PP8BNN8G4cTWeXtcIYjUQSyIU/EVS6YknfAIffBJ/6VI4+OBa\n3xI5gnjOnKrHyxuIx4zx68eLRCPR4J8RXzrz8/MpKChIdzVE6nbeeXDccX57x46K6zzXoGFDOOss\neP55/6Xh+usrZot27IAZM3zjcJ8+fnnhHTuSU32p/woKCsjPz0+4HD35i8Tq1VeD6SQaNPCDz3r2\njKmIoiLfJqAGYomX0j4i6fCTn8DLL/vtc86Bf/0rrmLKG4jL1yDetavicTUQS00U/EXSYdGiIP0D\nft6hE09MqMi1a4MG4i++qHpcDcQSScFfJF2GDPGd+wFOOglefz2UBelLSnz7QPkI4srUQCyg4C+S\nPoWFPtdfXOz3n30WBg4M/RKTJ8P06RpBLBUp+Iuk0xVX+IQ9+JW+3n/fd+8JWVFRMIJ40aKqx9VA\nnHsU/EXS6dtv/bQP27b5/ZkzYfjwpF1ODcRSTsFfJN3y8vwykgAdO8Jnn4WXgzHzHzB77+1/IqiB\nOLcp+Iuk2+bN/ul/3Tq/f/fdfg2AWOza5deQ/PTTqj+bN/s1I+fMgZNPrvLW8hHEEyeqgTiXZEXw\n12IuUu/ddx9cfbXfbtPGB/JWraqet2FD9QF+2TIfxWvTtSssWQJNm9Z4SmEh3H+/ppjOZlrMRSST\n7NwJPXoEM7hddZVPulcO8mvWJHadceP8hHJ1UANx9suKJ/9010EkFH//OwwdGv/7DznEf4BU/pk1\nC665xp+z116+TSHKRevVQJy9FPxFMkVpKRx7LHzwQc3nNGsG3btXDfDdukGLFtW/Z/duOOYY+Ogj\nv3/hhf6DJkZJayD+8kuYNw/OPLP6VJckhYK/SCZZtAgGD4atW6t/iu/YMb7H67lz4ZRTgv033vCj\niuMQagNxYaGf5mLzZj/i+dFH46qTxE7BXyRXDB4MTz7pt486yk8HmuCAsoQaiEtK/AfQ22/7/aZN\nYdOmWhukJTxZMZ+/iEThrrt8l0/wqaUHHki4yMMO88WuWuU/ACLnqgNYuBAuusivVzN2LKxYEXHw\nzjuDwA++0bu6+aklI+nJX6Q+ue02+P3v/Xbr1vD557D//qEVH3UD8Zkr6H9NDxoUVzrhjjvgxhtD\nq4/UTGkfkVyyYwccfrhfPB7g0kt93iYJ6mwgZhmXM5mLmz5Gm52r/Ytnnw1PPZWU+khFCv4iueaZ\nZ2DQIL/tHLz7ru9llCR1NhA3KeGCXQ8xmkn0bvuFn44ihKmtpXYK/iK5xsx3qyyPxD/4AfznPykJ\nuIVPvs/9g19jGhezgf2qHO/NAsaM68j5Vx+oEcRJpuAvkos+/xy+/30/BgDgoYdg2LDkXrOoCI4+\nGpYupYjmzOqRz8S9buS996rGH40gTr6M7u3jnOvhnJvsnHvCOffLZF5LJKd061Zx8rixY31f+2S6\n6SZYuhSAFi0bcfEL57NwoWP+fBh29BKasHPPqevX+85AXbr49W2ef96PgZPMkZInf+dcA+BxMxtc\nzTE9+YvEY8sWP3BsdVlj6w03+IibDK++CqedFuxPneof7cu9/jpr+53HNEYyufFVfLG76vQTmmI6\nXClJ+zjnpgEDgTVmdkTE6/2Be4GGwN/MbHw17/0pMBqYYmb/qua4gr9IvB59FH7xC7/dqJGf9bNH\nj3CvsWmTX6Xsq6/8/lln+UbnyDaGoiI/tUNxMSU0YM6jG5n4UEtNMZ1EqUr7TAf6V7pwQ2BC2eu9\ngCHOuZ7OuWHOuXuccwcBmNlsMxsAjIi3kiJSgyFDgmkeiov9tNJhP0xdfXUQ+Nu0gSlTqjYut2ix\np8dRQ0o5q+XrzJnjs0TXX++HJJTbsQNmzPCjh/v08Yuf7dgRbpWr0ANmFVGnfZxznYDZ5U/+zrkT\ngTwz61+2/2sAMxsX8Z6TgZ8BzYBPzOzeasq1vLy8Pfua118kRh984ANveVL9qad8f/sw/PvfcM45\nwf6TT8J551V/7nXXwT33+O1f/xpuv33PoaIiP2hs4sTqBwEntYF42jS48kro1cvX65xz6uX0peXz\n+JdLdD5/zCyqH6ATsCRi/zx8Kqd8fyjw12jLi3ifiUiCRo8288+3Zp06mRUVJV7mt9+atW0blHvh\nhbWf/49/BOeedFK1p5SWms2fbzZsmFmTJsHp5T/OmZ15ptlzz5mVlCT+K1hpqdn3vlfxIr16mT3y\niNnu3SFcIH3KYmdM8TbyJ5GPP32PEskUt94atKKuXJl4w68ZXHaZH+YLcNBBMGFC7e/54Q+D7YUL\nq84Ngc8W9e3re6auWuXXpjnkkIqXff5530Ooa1c/79D69Qn8Hh99BN98U/G1jz/26y706OEbrqup\nZy5IJPh/DXSI2O8ArIqnoPz8/ApfZ0QkRvvtB3/6U7B/++3Vz8kQrUceqThNw9SpFRP31TnwQN+l\nB3wSv7olxCK0bet7jy5bBrNn+8VkIi1f7qcJOvhg30Po3Xfj+D1efjnY7twZ9t472F+2DH71Kz+7\n3YQJsH17HBdIvYKCAvLz8xMvKNqvCFRN+zQClpW93gRYDPSM9asHSvuIhKO42OzYY4P0xrnnxlfO\nl1+atWoVlHPZZdG/d9iw4H133RXzpZcuNbv+erPWraumhMCsd2+zGTPMtm+PssAzzwzefP/9ZuvX\nm+XnV3+BAw4wu/NOsy1bYq53OpBg2ifaAP0YsBrYCXwFXFz2+gDgM6AQ+E1cFVDwFwnPvHkVA9or\nr8T2/tJSs9NPD95/6KGxBcP77w/ee845sV07wrZtZtOmmR13XPUfAm3amN14o9ny5bUUsnOn2V57\nBW9atiw4tmmT2bhxZu3aVS18v/3Mbr3VbMOGuOufCikJ/sn8ASwvL8/mzp0b+s0RyUmRT989e5rt\n2hX9eydODN7rnNmbb8Z27Q8/rPgkXVoa2/srSaiBuKAgOLFz5+ovsG2b2V/+Yta+fdXC99nH7Le/\nNVu7NqHfIWxz5861vLy87Aj+IhKi1avNWrYMgtjdd0f3vs8/N2vePHjfjTfGfu2Skoopo8LC2Muo\nwZo1/mH9kEOqxuny+H7HHWbr1pW94Xe/Cw6OGlV74Tt2mD3wgC+kcsEtWphdd52/rxlEwV9Eqrrz\nzopPsP/7X+3nFxebnXhi8J7DD48hsV5J//5BOTNnxldGLYqLzWbPrniZyJ9mzcwuushsweEXBS8+\n+WR0he/ebfbQQ2bdu1ctuGlT36V25crQf6d4ZEXwV9pHJGQ7d1YMYBddVPv5t98enNuokdmiRfFf\n+9Zbo3/iTlCdDcS8Y9O5yIq+Wld3YZGKi82eeMLsyCOrFtqokdnIkf6bUhoo7SMitXvxxYpB6+23\nqz/vgw/MGjcOzrv11sSu+9prQVmHH55YWVGqq4F4v/2iaCCuTmmp2TPPmPXpU7XQBg18ailNFPxF\npGZnnx0Eq+OPr9oqunOn2VFHBef06ZP4yNetW80aNgzK/O67xMqLwZ4G4m7zrQk7Ymsgrqvgl14y\n+9GPKhb46KNJ+13qouAvIjVbvtwnwcuD1ZQpFY/ffHNwrFkzs08+Cee6xx8flPvcc+GUGYvDDrM1\n7G/jGGuHHFAUXQNxtN54w+yMM8y6dvXpoTTJiuCvnL9IEv3hD0HE23//4En8rbd86qL82L33hnfN\nq68Oyr355vDKjcaKFcG1mze34m07omsgXhDjddI0DkA5fxGJzrZtZh07BtHuyit9aqZr1+C1H/84\npJnUyjy8VMSaAAAM70lEQVTxRFD2ySeHV240pkwJrn3GGRUOhT6COI0SDf71b15TEYlNixZw993B\n/sSJfh2AsiUZadkSpk8Pd5rjyEneFiwI1hpOhcj5fCJXH8NP43PXXX5SualT9yxBsMfChXDRRX4+\nobFjYcWK5Fc3XbSAu0guMIPTT/fLMVY2bZqfOS1snTv7GUYB3nnHr9ySbKWl0K5dMBXo4sVw1FE1\nnm7mP5smTvTrDVSe4NM5GDDArzrWv39mLQOQ0Qu4R0uzeookmXNw331+qcdIP/2pf9RNhsin/3nz\nknONyhYvDgJ/u3Z++clapHyK6RCkfFbPZP2gnL9I6lx7bZDgbtPG7JtvknetSZOCa8U7w2isxo0L\nrjlkSFxFRDuCeOHCkOseI5TzF5Go5ef79E/Hjj7PceCBybtW5Sf/VKR3X3kl2D799LiKaNjQr1Ef\nzRrE116bWHXTSTl/EUmOkhK/yMzmzX5/+fIkLNAbYft2H6V37vT7X34JHTrU/p4oFRXB44/7toHI\nNWqefdanhNIhK3L+IpKFGjaEE04I9pOd9//Pf4LA3717aIEffIepkSP9amLz58OwYX4VyMqrj9Un\nCv4ikjypbPQNIeVTl8gG4iVL/OdbfZURwV+9fUSyVCqDf2T//iQF/0iVO06lSli9fZTzF5Hk2boV\n9t3X5/+dg+++8/thW7vWd+0E/zi+fj20ahX+dTKIcv4ikrn23jsYZGXmE+bJ8NprwXbfvlkf+MOg\n4C8iyZWK1E+KUz7ZQMFfRJIr2cHfrNb5fKR6yvmLSHKtWhV0u2zRAjZuhMaNwyt/6VLo1s1vt2zp\n8/1hlp+hlPMXkcx28MF+RDH40VIffBBu+ZFP/f365UTgD0NGBH919RTJcslM/eRYykddPUWk/pg4\nEa64wm///OfwxBPhlFtcDPvvD5s2+f2PP4aePcMpO8Mp7SMimS9Zk7y9+24Q+Nu393MuSFQU/EUk\n+Y44wjfGAqxe7SddC0PllI+L+0E45yj4i0jyJWuStxTM55OtFPxFJDXCbvTduhXefjvYz4HG3jAp\n+ItIavzgB8F2GMH/9deDheGPOAIOOCDxMnOIgr+IpMYJJwQroC9ZEizyEi+lfBKi4C8iqdGyJRx5\npN8uLU18kjfN55OQjAj+GuQlkiPCyvuvXg0ffeS3mzSBk05KrF71iAZ5iUj989hjcOGFfvvUUyum\nbmLx8MMwfLjf7tcP5s4NpXr1iQZ5iUj9EfnkP3++H6EbD6V8EqbgLyKp07Gjn+gNYNs2+PDD2Msw\nq/iNQV0846LgLyKplWje/+OP4Ztv/Hbr1nDcceHUK8co+ItIaiUa/CNTPqec4kcPS8wU/EUktcIM\n/kr5xE29fUQktYqLYd99fc4f/CRv5St91WXXLthvv+C9hYXQpUty6pnh1NtHROqXRo3in+Rt/vwg\n8HfunLOBPwwK/iKSevGmfpTyCY2Cv4ikXrzBX/P5hEY5fxFJvc2bfd7fzE/2tnFjsNhLTTZuhDZt\n/LxAzsHatX4/R2V8zt85t5dzbqFzbmCyryUi9cQ++/hpmMEH83feqfs9c+f6cwGOPTanA38YUpH2\nGQvMSsF1RKQ+iTX1o5RPqKIK/s65ac65b51zSyq93t8596lzbqlz7qZq3nc68DGwNpzqikjWiDX4\naz6fUEWV83fOnQRsBR4ysyPKXmsIfAacBnwNLASGAMcDxwJ3AqOBvYBewHbgnMoJfuX8RXLUypW+\nuyb4fP+GDTWP1v3iC+jUyW83a+bPbdYsFbXMWInm/BtFc5KZvemc61Tp5T5AoZmtLKvI48AgMxsH\nPFx2zu/Kjo0A1irKi8gehxwCBx3k5+bfssWv7nX00dWfG5ny+dGPcj7whyGq4F+D9sBXEfurgL7V\nnWhmM2srKHJhgn79+tGvX78EqiUi9YJzPvXz5JN+f968moO/Uj4UFBSEuuhV1F09y578Z0ekfc4F\n+pvZJWX7Q4G+ZnZlTBVQ2kckd/3lL3DNNX57yBB49NGq55SW+sXZ163z+++/X/OHRA5JZ1fPr4HI\nCTk64J/+RUSiE02j7+LFQeBv2zZYB1gSkkjwfxfo6pzr5JxrApwPPBNPQVrDVyRHHXUUtGjht7/8\nElZV8/xYeeGWBrk9MUFK1/B1zj0GnAy0AdYAfzCz6c65AcC9QENgqpndHnMFlPYRyW2nnBKswTtr\nFgweXPH46acHHwBTp8LIkamtX4ZKVW+fITW8PgeYE+/FRUT44Q+D4D9vXsXgv307vPlmsJ+jjb3J\nkBHfn5T2EclhteX9582DnTv9dvfu0c/7n8VSmvZJJqV9RHLcpk1+LV4zP8hr40bYe29/7Kab4I47\n/PaYMTBhQvrqmWEyfmI3EZFatWoF3/++3y4pgQULgmPq3580GRH8lfYRyXHVpX7WrfN9+sF/I9Dg\nT0BpHxHJJg8/DMOH++0zzoAXXvA9fy64wL924onw1lvpq18GUtpHROq/yCf/t9/26R+lfJJKwV9E\n0q9zZzjwQL+9eTN89JGCf5JlRPBXzl8kx5VP8lZu5kw/4hd8z5++1c4ZmZOU8xeR7HLPPXDddX67\neXM/wAvgrLNg9uz01StDKecvItkh8sm/PPCDUj5Joid/EckMu3f7Pv+RgR98/r9Xr/TUKYNlxZO/\ncv4iQuPG0KdPxdcOOgh69kxPfTKUcv4ikn1++1v4v/8L9keMgBkz0ladTJYVT/4iIkDFvD/4+fsl\nKfTkLyKZY8MGaNPGT/IG8M03Qf9/qUBP/iKSPVq3hquu8qt1jRmjwJ9EGfHkn5eXR79+/einiZtE\nBHyPn+bN012LjFRQUEBBQQG33HJLQk/+GRH8010HEZH6RmkfERGJmYK/iEgOUvAXEclBCv4iIjlI\nwV9EJAdlRPDX3D4iItHR3D4iIjlMXT1FRCRmCv4iIjlIwV9EJAcp+IuI5CAFfxGRHKTgLyKSgxT8\nRURykIK/iEgOyojgrxG+IiLR0QhfEZEcphG+IiISMwV/EZEcpOAvIpKDFPxFRHKQgr+ISA5S8BcR\nyUEK/iIiOUjBX0QkByn4i4jkIAV/EZEclNTg75zr55x70zk32Tl3cjKvJSIi0Uv2k38psAVoCqxK\n8rXqPU1uF9C9COheBHQvwhNV8HfOTXPOfeucW1Lp9f7OuU+dc0udczdV89Y3zexM4NfALSHUN6vp\nDzugexHQvQjoXoQn2if/6UD/yBeccw2BCWWv9wKGOOd6OueGOefucc4dFDFd50b8039KxPIHEs25\nNZ0T7eu17Sf7j1n3ouZrJ3puLPcimtd0L6rfT+a9iLXsbLoXUQV/M3sT2FDp5T5AoZmtNLPdwOPA\nIDN72MyuNbPVzrlznHP3Aw8Bfw2t1nVQwKv52omeq3tR9zmZ9j95dXQv4is7m+5F1PP5O+c6AbPN\n7Iiy/fOAM8zskrL9oUBfM7sypgo4p8n8RUTikMh8/o0SuW4C7w0KSaDyIiISn0R6+3wNdIjY74B6\n9IiI1AuJBP93ga7OuU7OuSbA+cAz4VRLRESSKdquno8BbwHdnHNfOecuNrNi4ArgReBjYJaZfZK8\nqoqISFjSvoC7iIiknub2ERHJQRkX/J1zg5xzDzrnHnfOnZ7u+qSTc65H2bxITzjnfpnu+qSbc24v\n59xC59zAdNclnTRnVsB5f3LO3eecG57u+qSTc+7/lf1NTHHOzavr/ES6eiaFmT0NPO2c2xe4C3g5\nzVVKGzP7FLjcOdcAP4huapqrlG5jgVnprkQG0JxZgbOB9sA6cvxemNl/gP845wYBC+o6PyVP/nHO\nDfQ7/PQRWSXWe+Gc+ynwHD74Z5VY7kXZt8CPgbXpqGuyxfh3kdVzZsV4L7oB88zsBuDylFc2yeKM\nnRcCj9ZZuJkl/Qc4CTgGWBLxWkOgEOgENAYWAz0BB4wHTk1F3VL9E8u9qPS+p9Nd9zT/XdwG3IPv\nXfZvyjorZMtPPH8XQBPgyXTXPc1/F78Afl52zqx01z3dfxdAR+DBaMpOSdrHzN4smx4i0p65gQCc\nc48Dg4DTgFOBfZxzh5nZA6moY6rEci+cc+2AnwHNgLkprGZKxHIvzOx3ZfsjgLVW9peeLWL8u+gB\nnAHsSwrnzEqVGOPFX4C/OudOAgpSV8vUiPFefAKMBKZFU3Y6c/7tga8i9lcRzA2UdX/QdajpXrwO\nvJ6eKqVNtfeifMfMZqa8RulT09/FOOCp9FQpbWq6F9uBX6WnSmlT4/8jZpYfbSHp7O2TVU9uCdK9\nCOheBHQvAroXgVDuRTqDv+YGCuheBHQvAroXAd2LQCj3Ip3BX3MDBXQvAroXAd2LgO5FIJR7kaqu\nnpobqIzuRUD3IqB7EdC9CCTzXmhuHxGRHJRx0zuIiEjyKfiLiOQgBX8RkRyk4C8ikoMU/EVEcpCC\nv4hIDlLwFxHJQQr+IiI56P8DlFWaM7iih3kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10969b650>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Simple Monte Carlo Integration\n",
    "ngroups = 16\n",
    "\n",
    "I = np.zeros(ngroups)\n",
    "N = np.zeros(ngroups)\n",
    "E = np.zeros(ngroups)\n",
    "\n",
    "n0 = 100\n",
    "for i in range(ngroups):\n",
    "\n",
    "    N[i] = n0\n",
    "    r = np.random.random(n0)\n",
    "    I[i] = 0.\n",
    "    for j in range(n0):\n",
    "        x = r[j]\n",
    "        I[i] += np.sqrt(1-x**2)\n",
    "        \n",
    "    I[i] *= 4./float(n0)\n",
    "    E[i] = abs(I[i]-np.pi)\n",
    "    print n0,I[i],E[i]\n",
    "    n0 *= 2\n",
    "    \n",
    "            \n",
    "pyplot.plot(N,E,ls='-',c='red',lw=3);\n",
    "pyplot.plot(N,0.8/np.sqrt(N),ls='-',c='blue',lw=3);\n",
    "pyplot.xscale('log')\n",
    "pyplot.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.14093947721 10.6643035615 0.798802762018\n",
      "3.1499382371 10.7055506255 0.783439727913\n",
      "3.13319009661 10.6235086554 0.806628473869\n",
      "3.13795800203 10.6515853995 0.804804977031\n",
      "3.13531651984 10.6418893471 0.811679667517\n",
      "3.13173378372 10.6092906027 0.801534110608\n",
      "3.13380941075 10.6278067095 0.807045286608\n",
      "3.16020428513 10.7543778246 0.767486700818\n",
      "3.14413190095 10.6786028396 0.793037429062\n",
      "3.13419696292 10.6266539235 0.803463321146\n",
      "=============================\n",
      "3.15153961803 10.7326143302 0.800412366189\n",
      "3.13033933639 10.5959927928 0.796968431876\n",
      "3.13355454219 10.6259457425 0.806781673606\n",
      "3.16632193201 10.7851555084 0.759560931301\n",
      "3.12325523145 10.5737464366 0.819023195819\n",
      "3.14312496178 10.6732708742 0.794036348827\n",
      "3.12971418414 10.6123016964 0.817190822021\n",
      "3.14620181992 10.6908691026 0.792283210975\n",
      "3.13347136259 10.6306027368 0.811959956686\n",
      "3.13716167709 10.6531759573 0.811392569138\n",
      "3.13254956736 10.6148269262 0.80196013421\n",
      "3.13091800008 10.6037542792 0.801106756001\n",
      "3.13850922008 10.6499766584 0.799736533867\n",
      "3.12910960142 10.6056367606 0.814309862934\n",
      "3.16758589926 10.7971696007 0.763569171544\n",
      "3.15282267101 10.7115860485 0.771295253638\n",
      "3.1307753286 10.626421502 0.824667343858\n",
      "3.1574884733 10.7307841772 0.761050718215\n",
      "3.11465580988 10.5383340428 0.837253228776\n",
      "3.15373811595 10.7149738042 0.768909700192\n",
      "=============================\n",
      "3.14014186763 10.6583569489 0.797866000074\n"
     ]
    }
   ],
   "source": [
    "n0 = 100000\n",
    "I = np.zeros(n0)\n",
    "r = np.random.random(n0)\n",
    "for j in range(n0):\n",
    "    x = r[j]\n",
    "    I[j] = 4*np.sqrt(1-x**2)\n",
    "\n",
    "def group_measurements(ngroups):\n",
    "    global I,n0\n",
    "    \n",
    "    nmeasurements = n0/ngroups\n",
    "    for n in range(ngroups):\n",
    "        Ig = 0.\n",
    "        Ig2 = 0.\n",
    "        for i in range(n*nmeasurements,(n+1)*nmeasurements):\n",
    "            Ig += I[i]\n",
    "            Ig2 += I[i]**2\n",
    "        Ig /= nmeasurements\n",
    "        Ig2 /= nmeasurements\n",
    "        sigma = Ig2-Ig**2\n",
    "        print Ig,Ig2,sigma\n",
    "        \n",
    "group_measurements(10)\n",
    "print \"=============================\"\n",
    "group_measurements(20)\n",
    "print \"=============================\"\n",
    "group_measurements(1)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Variance reduction\n",
    "------------------\n",
    "\n",
    "If the function being integrated does not fluctuate too much in the\n",
    "interval of integration, and does not differ much from the average\n",
    "value, then the standard Monte Carlo mean-value method should work well\n",
    "with a reasonable number of points. Otherwise, we will find that the\n",
    "variance is very large, meaning that some points will make small\n",
    "contributions, while others will make large contributions to the\n",
    "integral. If this is the case, the algorithm will be very inefficient.\n",
    "The method can be improved by splitting the function $f(x)$ in two\n",
    "$f(x)=f_1(x)+f_2(x)$, such that the integral of $f_1(x)$ is known, and\n",
    "$f_2(x)$ as a small variance. The “variance reduction” technique,\n",
    "consists then in evaluating the integral of $f_2(x)$ to obtain:\n",
    "$$\\int _a^b{f(x)dx}=\\int _a^b {f_1(x)dx} + \\int _a^b{f_2(x)dx} = \\int\n",
    "_a^b{f_1(x)dx}+J.$$\n",
    "\n",
    "Importance Sampling\n",
    "-------------------\n",
    "\n",
    "Imagine that we want to sample the function $f(x)=e^{-x^2}$ in the\n",
    "interval $[0,1]$. It is evident that most of our points will fall in the\n",
    "region where the value of $f(x)$ is very small, and therefore we will\n",
    "need a large number of values to achieve a decent accuracy. A way to\n",
    "improve the measurement by reducing the variance is obtained by\n",
    "“importance sampling”. As the name says, the idea is to sample the\n",
    "regions with larger contributions to the integral. For this goal, we\n",
    "introduce a probability distribution $P(x)$ normalized in the interval\n",
    "of integration $$\\int _a^b{P(x)dx} = 1.$$ Then, we can rewrite the\n",
    "integral of $f(x)$ as $$I=\\int _a^b{\\frac{f(x)}{P(x)}P(x)dx}\n",
    "$$ We can evaluate this integral, by sampling\n",
    "according to the probability distribution $P(x)$ and evaluating the sum\n",
    "$$I(N)=\\frac{1}{N} \\sum_{i=1}^N \\frac{f(x_i)}{P(x_i)}.\n",
    "$$ Note that for the uniform case $P(x)=1/(b-a)$, the\n",
    "expression reduces to the simple Monte Carlo integral.\n",
    "\n",
    "We are free to choose $P(x)$ now. We wish to do it in a way to reduce\n",
    "and minimize the variance of the integrand $f(x)/P(x)$. The way to to\n",
    "this is picking a $P(x)$ that mimics $f(x)$ where $f(x)$ is large. if we\n",
    "are able to determine an apropiate $P(x)$, the integrand will be slowly\n",
    "varying, and hence the variance will be reduced. Another consideration\n",
    "is that the generation of points according to the distribution $P(x)$\n",
    "should be a simple task. As an example, let us consider again the\n",
    "integral $$I=\\int _0^1 {e^{-x^2}dx}.$$ A reasonable choice for a weigh\n",
    "function is $P(x)=Ae^{-x}$, where $A$ is a normalization constant.\n",
    "\n",
    "Notice that for $P(x)=f(x)$ the variance is zero! This is known as the\n",
    "zero variance property. There is a catch, though: The probability function\n",
    "$P(x)$ needs to be normalized, implying that in reality, $P(x)=f(x)/\\int f(x)dx$, which\n",
    "assumes that we know in advance precisely the integral that we are trying to calculate!\n",
    "\n",
    "### Exercise 10.3: Importance sampling \n",
    "\n",
    "1.  Choose the weight function $P(x)=e^{-x}$ and evaluate the integral:\n",
    "    $$\\int _0^{\\infty} {x^{3/2}e^{-x}dx}.$$\n",
    "\n",
    "2.  Choose $P(x)=e^{-ax}$ and estimate the integral\n",
    "    $$\\int _0^{\\pi} \\frac{dx}{x^2+\\cos ^2{x}}.$$ Determine the value of\n",
    "    $a$ that minimizes the variance of the integral.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trapezoids:  1.32934018965\n",
      "Simple Monte Carlo:  1.32519811567\n"
     ]
    }
   ],
   "source": [
    "# Trapezoidal integration\n",
    "def trapezoids(func, xmin, xmax, nmax):\n",
    "    Isim = func(xmin)+func(xmax)\n",
    "    h = (xmax-xmin)/nmax\n",
    "    for i in range(1,nmax):\n",
    "        x = xmin+i*h\n",
    "        Isim += 2*func(x)\n",
    "\n",
    "    Isim *= h/2\n",
    "    return Isim\n",
    "\n",
    "def f(x):\n",
    "    return x**1.5*np.exp(-x)\n",
    "    \n",
    "    \n",
    "print \"Trapezoids: \", trapezoids(f, 0., 20., 100000)\n",
    "\n",
    "# Simple Monte Carlo integration\n",
    "n0 = 1000000\n",
    "r = np.random.random(n0)\n",
    "\n",
    "x = -np.log(r)\n",
    "Itot = np.sum(x**1.5)\n",
    "print \"Simple Monte Carlo: \", Itot/n0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trapezoids:  1.58118797085\n",
      "Simple Monte Carlo:  1.58232612038\n",
      "Importance Sampling:\n",
      "0.1 1.51733652865\n",
      "0.2 1.50271113189\n",
      "0.3 1.49497765112\n",
      "0.4 1.49760515522\n",
      "0.5 1.51440492709\n",
      "0.6 1.53705127814\n",
      "0.7 1.55998782484\n",
      "0.8 1.57856622469\n",
      "0.9 1.5850784562\n",
      "1.0 1.57501692127\n",
      "1.1 1.54617823139\n",
      "1.2 1.4988445746\n",
      "1.3 1.43637280277\n",
      "1.4 1.36362800283\n",
      "1.5 1.28410162522\n",
      "1.6 1.20147734863\n",
      "1.7 1.12083134139\n",
      "1.8 1.04338661046\n",
      "1.9 0.970940233688\n",
      "2.0 0.904242852582\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VdW5//HPQxLCJKBSBlHBxgksCKio9BJiKwmatlpr\n1WptHVqrfWnw1tZaAiZeTf3ZWhVCtWjFWey1zkYZemuII3VARq0SioLWoEKRMCdZvz9W4iHJSXKS\nnJx9hu/79Tov99p7Z+fJZvtk59lrr2XOOUREJLl0CzoAERGJPiV3EZEkpOQuIpKElNxFRJKQkruI\nSBJSchcRSUJtJnczm2tmVWa2opV9csxsqZmtNLPyqEYoIiLtZm31czeziUA1cL9zblSY7f2Bl4E8\n59wGMxvgnPusS6IVEZGItHnn7px7Edjcyi7nAo855zbU76/ELiISsGjU3A8D9jOzF8zsDTM7PwrH\nFBGRTkiPwjEygHHAN4FewKtm9ppz7v0oHFtERDogGsl9PfCZc24HsMPMKoCjgUbJ3cw0iI2ISAc4\n56y9XxONssxTwH+ZWZqZ9QKOB1aH29E5p0+UPkVFRYHHEOTn2WcXk5tbyKRJReTmFvLss4t1LuPk\no/MZ3U9HtXnnbmbzgEnAADNbDxThSzE45+Y45941s/nAcqAOuMs5Fza5i0RDWVkFU6cuoLKy5Mt1\nlZWFAOTnZwcVlkhcaTO5O+d+EME+NwM3RyUiSQllZRXMmrWQXbvSycysoaAgt1li3rYNPvgA1q3z\nn4bl+fMX8sUXJY32raws4dJLZ1BYmM24cTBqFPTsGbMfRyTuRKPmLgHIyckJOoQOC3fn/dZbhUyc\nCGlp2V8m889a7FQb/rLdsCGNyy7zy2lpcOSRMHZs6DNmDOy7b/OvS+RzGY90PuNDmy8xRe0bmblY\nfS+JX1u3woQJ01m58oYwW2cA10dwlOlAx75++PBQsq+traC8fCHOpdOjR/i/HkSCZma4DjxQ1Z27\ndLndu2H+fHjoIXj6adi5s6XLLq1RKyMDDj7YJ+SGz7Bh8O9/53L77YV88EHozn/IkGnk5k5h505Y\nuhTefx/C3Us0/FXwxBMVwAIgdIw1a1S3l+Sh5C5doq4OXnrJJ/RHH4XNjd5xrgn7NYcdVktRUSiR\nDxkC3cL258rmqKOgtHQGO3em0aNHLVdcMaVRUq6uhmXLfKJv+KxcCXv2NOyxkL0TO8DatSVMmzaD\nU07JbuH7iiQOlWUkapyD5cvh4Ydh3jxYvz78fsOGVbB16wI2bQol16ysacycOaVL75p374bVq32i\nnz69mI8/Lg6zVzHHHlvM738PKh1LPFBZRmKmaU+Xs8/O5ZNPsnn4YVi1KvzXDBsG557rP1/7WjZl\nZa3feXeF7t39Q9UxY+CRR2r4+ONwe9Xyxhtw0knwrW/BTTfByJFdGpZIl9Cdu7RLuJ4uUAjkAY2T\n8/77w1ln+YQ+YUJLJZZghPs5+vefxrZtU9izJ/RzdOsGF18M113ny0QisdbRO3cld2mXCROm8+qr\nLfdU6dULTjsNzjsPJk/2d8vxqqysgtLSRXv99TCZ0aOzmTED7r+/8QPZXr3gl7/0n332CS5mST1K\n7tJlamt9L5eZM2Hx4mKguNk+++5bTGlpMaedBn36xDrC6Hv7bbj6ali0qPH6QYOguBh+8hNIV1FT\nYqCjyT2O/lCWeLNlC9x6Kxx2GJxxBixeDC31dBk/vpbzzkuOxA6+Lr9wISxYAKNHh9ZXVcFll8HX\nvgZPPRW+u6VIPNCduzSzZg2UlsLcub5L4d66daugV68FVFfHtqdLkGpr4cEHYfp02LCh8baRIyvY\nZ5+F9OjR8jAKIp2hsox0inNQXg633QbPPNP8jnS//eBnP4Of/xyWLWteq06FhLZjhz8/N97o37SF\n5i9CZWUVMnNmXkqcD4kNJXfpkJ07fb/0mTN9H/WmRo6EK6/0D0h79Yp9fPHo00/h+uuhtDT8MAh5\neTOYPz+SYRRE2qZ+7hKRhj7qX3yRzsaNNXz2WS5ffNH8LvPUU31SP/lksHZfVsntK1+BWbNgyZJ0\n/vGP5ts3bUprvlIkxpTcU0hZWQWXXrqADRua9lEHyKZXL7jgAigogCOOCCDABNO/f/iHy0uX1rJk\nCRx/fIwDEtmLesukgD174C9/gfPOW9gksQOU0KPHIn73O/+w8I9/VGKPVEFBLllZhU3WTqOmZjLf\n+AY8/3wgYYkAunNPap9+CnfdBbffDh99BC39c48fn8avfhXT0JJCw0PThmEUdu2qZdWqKWzdms32\n7fDtb8Pdd8OPfxxwoJKSlNyT0Ntv+5rwww/Drl17bwlfRujZszYmcSWj/PzsRj1j3nsP8vL8sMK1\ntb7M9e9/w69/rWcXElsqyySJmhp47DGYNMlPRHHPPY0T+6BBcO65uQwf3riMkJU1jSuumBzjaJPX\n4YfDK6/A0UeH1v3mN/7hdF1dcHFJ6kmIrpCRzLeZqjZtgj//2dfKP/yw+fZjj4WpU/0AXt27hx9P\nRecy+rZsgdNP9+8ONDjrLD9mTWZmYGFJAkrafu5lZRWcf/4CNm/WiyIQ+kW3aVM6n3xSw8aNueze\n3fg8pKfDmWf6Xi8nnKByQFB27YLzz/eTlTT4xjfgiSegb9/g4pLE0tHkjnMuJh//rdovN7fQ+fcl\nG3+ysqa7l192bs+eDh02IT3yyGI3cOC0JudimoPFDpwbMMC5wkLnNmwIOlJpUFPj3OWXN752x4xx\n7t//DjoySRT1ubPdOTfua+6bN4d/5ltZmcbXv+5fKDnzTF+aaGnmn0RWVQV33glTpsA55yxk48bm\nXRn32WcRc+f6n/+GG2Do0EBClTDS0vzD7ZK9/tneftuPb//++8HFJcmvzd4yZjYXyAc2OudGhdme\nAzwFrK1f9ZhzLtyA3x3S0osi4Ht4/Oc//kHiY4/5tSNG+ESYlwfZ2dCzZ7QiiZ0PPvB/uj/+uJ+H\nNFTNCv/PNXZsGhdeGLPwpJ3MYNo0GDwYLrnE96L51798gn/uOTjuuKAjlGQUSVfIe4BS4P5W9lns\nnPtOdEJqbOrUXNauLWw0Y86gQdMYPXoKq1c39N8Oeecd/7n1VujRw/ceOfDACt55ZyHp6en06BGf\nD2Tffdcn88cfhzffbGkvdWVMZBddBAMH+gerO3bAZ5/56fwee8zfjIhEU0QPVM1sOPBMK3fuVznn\nvt3GMVwk3yuclnp4OOcnPF6wAObPh4qKpv26IdzIfX37FnLOOXmccUY2hx8OBx/s/3zuak17/Zx6\nai6ffprN44/7X0jhmMHEiX489b59KygpaTw1XLIPt5uMXn3Vz8+6aZNvd+tWwciRC9l/f/UGk+a6\ntLdMG8l9EvA4sAH4CPilc251mP06nNwjtX27T/Dz5/uE/+67AOFH7muYFg58F8GsLN9H+bDD/H8b\nlocM8Qm2I90xnfN3aFu3wtNPV/A//xNuXJfmc49mZMA3v+kT+ne+4/uoN1BXxuTwzjv+bn39eg0b\nLK0LMrnvA9Q657ab2SnATOfc4WH2c0VFRV+2c3JyyMnJaW+87fLBBzB5cjHvv18cZmsx4aaLa6p3\nbxg0qIKqqgVs27b3ZMqFnHhiHvvtl83WrT6BV1fz5XJDO/TiSuu/ZHr1glNO8Qk9Px/69WvfzyqJ\n56OP4Mgjp1NdrWGDJaS8vJzyvV6QuO6667quKyQwHFgR4b7/AvYLsz56fYPaoaWulEOHTnfZ2c4N\nHtx8W/NP+GPA9Ai+tuFTFHb9wIFF7oknnNu2LZDTIwGbMCH8dTFpUlHQoUmcoINdITs9toyZDcL3\npHFmNh7/18Cmzh43WgoKcqmsLGyhTu3bX3zhp5Z77z3/ef/90PJ//gMtP3duu1Dfowfssw9UV9ew\nY0fz7WPH1nL66e3+sSRJ9OkT/iF5XZ0ekkvnRNIVch4wCRhgZuuBIiADwDk3BzgTuMzMaoDtwDld\nF277NR25z9epGz+A7NsXxo3zn705B59/Dvn5NWEnZRgxopZrrvHJu+mnTx//3/T6M1xWlsvUqc1/\nyVxxxZSo/8ySOMLdfMA01q2bwpYtKs9Jx8X98APxoKysgqlTO99LRQ9DJZyG66KqKo1ly2pxbjKQ\nzZQpfj7bdI3dmtKSdmyZeKHELLHw0EPwwx+G2gUFfn5bSV1K7iJJYsYMP4xEgzvugEsvDS4eCZaS\nu0iSqKuDs8+Gv/7Vt9PS/LsbJ58cbFwSDCV3kSSyfbsfOuONN3y7Xz947TU48shg45LYU3IXSTIf\nfwzjx4fGT8rKgiVLYP/9g41LYqujyT3uh/wVSVUHHABPPw29evl2ZSV873uwe3ewcUliUHIXiWPj\nxsGDD4baixfDz3++9zDQIuEpuYvEue9+F268MdS++2645Zbg4pHEoJq7SAJwDi64wE+wDX6k0qee\ngm+3OtC2JAM9UBVJcrt2+e6QL73k2717w8svw9FHBxuXdC0ld5EU8OmncPzxfpo+gIMOgn/8w0/h\nJ8lJvWVEUsBXvuLHm+nb17fXr4fTTyfsiKOS2pTcRRLMUUfBX/4C3er/712yBC6+WD1opDGVZUQS\nVGmpH1jMzxO8kGHD0jniCM3BmmxUcxdJMc7Bt75VwXPPaQ7WZKaau0iKMYM9exayd2IHqKwsobR0\nUTBBSdxQchdJYLt3h5/JY/v2tqeAlOSm5C6SwDIzw8/Bun695mBNdUruIgmsoCCXrKzCJmunsW7d\nZBapMpPS9EBVJME1TAG5Y0caq1bV8vnnfg7WwYNh+XLfN14Sl3rLiAhVVX44gqoq387P9y89WbtT\ng8QL9ZYREQYNgnvvDbXLymD27MDCkQDpzl0kCf3iF3DrrX45M9OPPzN6dLAxScd02Z27mc01syoz\nW9HGfseZWY2ZndHeIEQkum68EcaO9cu7dsE55/h5WSV1RFKWuQeY0toOZpYG3ATMB1TdEwlYZibM\nmxeaou+dd+Cqq4KNSWKrzeTunHsR2NzGblcAfwU+jUZQItJ5RxwBs2aF2n/6EzzxRHDxSGx1+oGq\nmQ0FTgPuqF+lwrpInLjoIjjzzFD7Jz+BDRuCi0diJ/y7y+1zG3CNc86ZmdFKWaa4uPjL5ZycHHJy\ncqLw7UWkJWZw553+geqHH8KmTXD++fC3v0GaRiiIS+Xl5ZSXl3f6OBH1ljGz4cAzzrlRYbatJZTQ\nBwDbgZ86555usp96y4gE5KWXYNIkqKvz7ZISmDYt2JgkMoH1c3fOfdU5d4hz7hB83f2ypoldRIL1\nX/8FM2aE2tdeC6+9Flw80vUi6Qo5D3gFOMLM1pvZRWb2MzP7WdeHJyLRMn06fP3rfrm2Fs49F7Zs\nCTYm6Tp6iUkkhXzwgR+eoCGpn3cePPhgsDFJ6zT8gIi0adgw/4C1wUMPwQMPBBePdB3duYukoIsv\nhrlz/XKfPrB0KRx6aLAxSXgaFVJEIlZdDcccA++959vHHed71HTvHmxc0pzKMiISsT59/PAEGRm+\n/frrUFQUbEwSXbpzF0lhf/gD/PKXABXAQkaPTmfw4BoKCnLJz88OODoBlWVEpAPq6uDYYytYunQB\nUPLl+qysQmbOzFOCjwMqy4hIu3XrBn37LmTvxA5QWVlCaakmYU1kSu4iKS/8EFM7d2rwmUSm5C6S\n4jIza8Ku79atNsaRSDQpuYukuIKCXLKyCpusnUZt7eRA4pHo0ANVEaGsrILS0kV89FEaK1fWApOB\nbJ54Ak4/PejoUpt6y4hIVFxwAdx3n18eNAhWr4b99gs0pJSm3jIiEhW33gpDhvjlqiq48spg45GO\nUXIXkUb23dfPt9rggQegrCy4eKRjVJYRkbDOOw8eftgvH3AArFoF/fsHG1MqUllGRKJq1iwYONAv\nf/wxXHVVsPFI+yi5i0hY++8Pt98eas+dCwsWBBePtI/KMiLSqrPOgkcf9csHHQQrV0LfvsHGlEpU\nlhGRLjF7NgwY4JfXr4errw42HomMkruItGrgQCgtDbXnzIG//z24eCQyKsuISJucgzPOgCef9O3h\nw2HFCj/ph3QtlWVEpMuYwR13+D7wAOvWwW9+E2hI0gYldxGJyODBMHNmqD17NlRUBBePtK7N5G5m\nc82sysxWtLD9NDNbZmZLzexNM/tG9MMUkXjwwx9Cfn6ofdFFsH17cPFIy9qsuZvZRKAauN85NyrM\n9t7OuW31y6OAJ5xzh4bZTzV3kSTw0Udw1FGwZYtv//d/wy23BBtTMuuymrtz7kVgcyvbt+3V7AN8\n1t4gRCRxDB3aOJnfdhu88kpw8Uh4Uam5m9npZvYO8DxQEI1jikj8uvBCyMvzy8758syOHcHGJI2F\nnzyxnZxzTwJP1pdwHgCOCLdfcXHxl8s5OTnk5ORE49uLSIyZwZ13wte+Blu3wj//CcXFcNNNQUeW\n+MrLyykvL+/0cSLq525mw4FnwtXcw+xbCYx3zn3eZL1q7iJJZs4cuPRSv9ytG7z6KowfH2xMyaZL\nZ2JqLbmbWRaw1jnnzGwc8KhzLivMfkruIknGOTj55IY3Vivo3Xsh48al07NnDQUFueTnZwcdYsLr\naHJvsyxjZvOAScAAM1sPFAEZAM65OcD3gB+Z2R58r5pz2huEiCQmM/jzn2HEiAp27VrAtm0lvPii\n31ZZ6SfdVoIPhoYfEJFOGzFiOu++e0Oz9Xl5M5g///oAIkoeGn5ARAIzcGD4IsDOnWkxjkQaKLmL\nSKf16FHTwvraGEciDZTcRaTTCgpyycoqbLJ2Gt///uRA4hHV3EUkSsrKKpg1axEvvZTG9u21wGTy\n8rJ5/nn/4FU6pku7QkaDkrtIaliyBE480XeTBHjoITj33GBjSmR6oCoiceH44+Hyy0PtK6+Ezz9v\neX/pGrpzF5Go27rVjxy5fr1v//jHcO+9gYaUsHTnLiJxY5994I9/DLXvuw/+7/+CiycV6c5dRLrM\nWWfBo4/65awsP+9qz57BxpRodOcuInFn1izo398vV1bCddcFG08qUXIXkS4zeDD8/veh9s03w7Jl\nwcWTSlSWEZEuVVcHJ50Umkz7uOP80MBpGpkgIirLiEhc6tbNT+yRmenbr78OpaXBxpQKlNxFpMsd\ncQRMnx5qT58OH3wQXDypQGUZEYmJ3bth3DhYtcq3Tz0Vnn1WQxO0RWUZEYlr3bvDXXeFkvlzz8Ff\n/hJsTMlMyV1EYubEE+Gyy0LtqVNh06bg4klmSu4iElM33ghDh/rljRvhV78KNp5kpeQuIjHVt2/j\noQnmzoUXXggunmSlB6oiEojvfQ8ef9wvH3ooLF+uoQnC0QNVEUkopaXQr59fXrMGbmg+v7Z0gpK7\niATigAPgpptC7d/9zt+9S3SoLCMigamrg0mT4KWXfPv44+HllzU0wd66rCxjZnPNrMrMVrSw/Twz\nW2Zmy83sZTMb3d4gRCQ1NQxN0L27by9ZArffHmxMyaLNO3czmwhUA/c750aF2X4isNo5t8XMpgDF\nzrkTwuynO3cRCeu666C4GKCCbt0Wctxx6fTrV0NBQS75+dkBRxesLp0g28yGA8+ES+5N9tsXWOGc\nOzDMNiV3EQlr1y449NAKNmxYAJR8uT4rq5CZM/NSOsHHS2+Zi4HnonxMEUlymZkwZMhC9k7sAJWV\nJZSWLgomqASXHq0DmdlJwEXA11vap9j/3QVATk4OOTk50fr2IpLgevUKn4527kytp6vl5eWUl5d3\n+jhRSe71D1HvAqY45za3tN/eyV1EZG+ZmTVh1/foURvjSILV9Mb3ug7OTdjpsoyZHQw8DvzQObem\ns8cTkdRUUJBLVlZhk7XTOP74yYHEk+gi6S0zD5gEDACqgCIgA8A5N8fM/gx8F/iw/kv2OOfGhzmO\nHqiKSKvKyiooLV3E0qVpbNxYC0xm6NBsVq0Kvc2aarq0t0w0KLmLSKQ+/xxGjIBPP/XtSy6BOXOC\njSko8dJbRkSk0/bfH2bPDrXvvFMjR7aX7txFJC45B2ecAU8+6dtf/SqsWAG9egUbV6zpzl1EkoqZ\nH4qgoda+di3MmBFsTIlEyV1E4taQIXDLLaH2bbf58WekbSrLiEhccw7y8mBR/YuqI0fCW2/5t1pT\ngcoyIpKUzPwD1d69fXv1avjtb4ONKREouYtI3Bs+vHFC/+1vNbFHW1SWEZGEUFcHEyfCK6/49jHH\nwGuvQXrURsiKTyrLiEhS69YN7r47NLHHm2/CrbcGG1M8U3IXkYRx5JFQVBRqX3stvPdecPHEM5Vl\nRCSh7NkD48fD22/79sSJUF7u7+yTkcoyIpISMjJg7tzQJNovvpi64860RsldRBLO2LFw9dWh9tVX\nw4cftrx/KlJZRkQS0s6dPsm/+65vT5kCzz3n+8UnE5VlRCSl9Ojhe880JPP58+HBB4ONKZ4ouYtI\nwpowAa64ItS+8kqoqgounniisoyIJLTqahg1Ctat8+0zz4RHHw00pKjSTEwikrIWLYLcXIAKYCEj\nR6Zz4IE1FBTkkp+fHXB0ndPR5J7kL+6KSCqYPBlOPrmCv/1tAVDC6tV+gLHKSj/hdqIn+I5QzV1E\nkkJt7UKgpNG6ysoSSksXBRNQwJTcRSQp1NWFL0Ts3JkW40jig5K7iCSFzMyasOvT0mpjHEl8UHIX\nkaRQUJBLVlZhk7XTqKmZHEg8QWuzt4yZzQXygY3OuVFhth8J3AOMBQqdc39o4TjqLSMiXaqsrILS\n0kVs2JDGqlW1wGQgm//9X/j+94OOrmO6rCukmU0EqoH7W0juXwGGAacDm5XcRSQeXHgh3HuvX95/\nf1i1CgYNCjSkDumy4Qeccy8Cm1vZ/qlz7g1gT3u/uYhIV7ntNjjoIL/8+edw6aV+su1UoZq7iCSl\nfv382DMNnnwytcaeielLTMXFxV8u5+TkkJOTE8tvLyIpZvJkuOwyuOMO377iCjjpJDjwwGDjak15\neTnl5eWdPk5Eww+Y2XDgmXA19732KQKqVXMXkXhSXQ1HHw1r1/p2Xh48/3ziDA0cD0P+JsipEpFU\n0qePf7DakMwXLIC77go0pJiIpLfMPGASMACoAoqADADn3BwzGwy8DvQF6oCtwEjnXHWT4+jOXUQC\nc9VVcMstfrlPH1i+HA45JNiYIqFRIUVEWrFjB4wbF5q5adIk+Pvf439i7Xgoy4iIxK2ePeG++0LJ\nfPFimD072Ji6kpK7iKSM8ePhmmtC7WuugffeCy6erqSyjIiklF27fJJfvty3TzgBXnoJ0uJ08EiV\nZUREIpCZ6csz6fVv+bz2Gtx8c7AxdQUldxFJOWPGQFFRqH3ttbByZXDxdAWVZUQkJdXUwIknwhtv\n+PbYsbBkCWRkBBtXUyrLiIi0Q3q6L89kZvr20qXw298GG1M0KbmLSMoaORJuuCHUvuEGeOut4OKJ\nJpVlRCSl1db6F5peftm3jzoK3nwzdEcfNJVlREQ6IC3Njz3Tq5dvr1rV+GFrotKdu4gI8Mc/wuWX\nA1QACxkzJp2BA2soKMglPz87sLg0toyISCfU1cHYsRUsX74AKPlyfVZWITNn5gWW4FWWERHphG7d\noF+/heyd2AEqK0soLV0UTFCdoOQuIlKvW7fwk9Pt3BmnYxO0QsldRKReZmZNC1tqYxpHNCi5i4jU\nKyjIJSursMnaaWzaNJnaBMvveqAqIrKXsrIKSksX8cknaSxbVgtMBrIpKYFp02Ifj3rLiIhE2bXX\nwvXX++W0ND808AknxDYGJXcRkSirqfFvr77yim8PHw5vvw39+sUuBnWFFBGJsvR0ePjhUDJftw4u\nvRQS4T5VyV1EpBXDhsFdd4XajzziR5OMdyrLiIhE4Cc/gbvv9su9e/vRIw8/vOu/r2ruIiJdaNs2\nOOYY+Oc/fXvcOF+L7+rRI7us5m5mc82sysxWtLLPLDN738yWmdnY9gYhIhLvevf2JZnu3X37rbeg\nsGmX+DgSSc39HmBKSxvN7FTgUOfcYcAlwB1Rik1EJK6MGQO/+12o/Yc/wPz5wcXTmjaTu3PuRWBz\nK7t8B7ivft8lQH8zGxSd8ERE4ktBAZx6aqj94x9DVVVw8bQkGr1lhgLr92pvAA6MwnFFROKOGdxz\nDwwe7NsbN/oEX1cXbFxNhR8Crf2aFvvDPjktLi7+cjknJ4ecnJwofXsRkdgZOBDuvx9yc317wQK4\n7Tb4xS86f+zy8nLKy8s7fZyIesuY2XDgGefcqDDb/gSUO+ceqW+/C0xyzlU12U+9ZUQkqfz616Ea\nfEYGvPqq71ETTUG+ofo08KP6IE4A/tM0sYuIJKPrr4djj/XLe/bAD34A1dXBxtSgzTt3M5sHTAIG\nAFVAEZAB4JybU7/PbHyPmm3Ahc65t8IcR3fuIpJ01qyBsWNDSf2CC3xNPlr0EpOISEAeeAB+9KNQ\n++GH/V18NCi5i4gE6Pzz4cEH/XLfvn70yEMO6fxxldxFRAL0xRe+PLN2LUAF/fotZPTodHr2rKGg\nIJf8/OwOHbejyT1aXSFFRFJa374wbx6ceGIFdXUL2LKlhBdf9NsqK/04BR1N8B2hIX9FRKJk/HjI\nyloIlDRaX1lZQmnpopjGouQuIhJFQ4aEL4js3JkW0ziU3EVEoqhHj5qw6zMyamMah5K7iEgUFRTk\nkpXVdCzgaTg3OabT86m3jIhIlJWVVVBauojKyjTWrKkFJgPZ3HGHn4O1PdQVUkQkzjjn31i9/37f\nzsiAigo44YTIj6HkLiISh3bsgAkT/EtNAAcc4GdxGhThrBdBDhwmIiIt6NkTHn8c9t3Xtz/+GM4+\n2w801pWU3EVEutghh/gXnKz+/nvxYj9ccFdSchcRiYG8PD9EcINbb/UTbncV1dxFRGKkrg7OOAOe\nesq3e/WC116DUc2mQQrRA1URkQSwZYsfpuC993z70EPh9dehf//w++uBqohIAujXzz9g7d3bt9es\n8cMFR3uCbSV3EZEYO+qoxrM1PfsslJS0vH9HqCwjIhKQX/0Kbr7ZL5tBWRmcckrjfVRzFxFJMDU1\nkJsLL7zg2/37wxtvQFZWaB8ldxGRBLRxIxxzDGzY4NtHHw2vvOJ70oAeqIqIJKSBA+Gxx6B7d99e\ntgwuuYRSNfPjAAAD3klEQVROjyCp5C4iErDx42H27FD7oYcatzuizeRuZlPM7F0ze9/Mmr0wa2b7\nmtkTZrbMzJaY2VGdC0lEJPX89Kdw8cWh9pVXVjB+/PQOH6/V5G5macBsYAowEviBmY1osts04C3n\n3NHAj4CZHY5GIlZeXh50CElD5zK6dD47bvZsOPZYAD/J9uuv39DhY7V15z4eWOOcW+ec2wM8ApzW\nZJ8RwAsAzrl/AsPN7Csdjkgiov+BokfnMrp0PjuuRw9ff8/IaD7Jdnu1ldyHAuv3am+oX7e3ZcAZ\nAGY2HhgGHNipqEREUtTBB8OIEeEn2W6PtpJ7JM9r/x/Q38yWApcDS4HYzgQrIpJEBg8OP8l2e7Ta\nz93MTgCKnXNT6tu/Aeqccze18jX/AkY556qbrFcndxGRDuhIP/e27v3fAA4zs+HAx8DZwA/23sHM\n+gE7nHO7zeynwOKmib2jwYmISMe0mtydczVmdjmwAEgD7nbOvWNmP6vfPgffi+be+jvzlcDFLR5Q\nRERiImbDD4iISOxE/Q3Vtl56qt9nVv32ZWY2NtoxJIsIXiDLMbMtZra0/tPxNx6SnJnNNbMqM1vR\nyj66LiPU1vnUtRk5MzvIzF4ws1VmttLMClrYr33Xp3Muah986WYNMBzIAN4GRjTZ51Tgufrl44HX\nohlDsnwiPJc5wNNBx5oIH2AiMBZY0cJ2XZfRPZ+6NiM/l4OBMfXLfYB/RiNvRvvOPZKXnr4D3Afg\nnFuC70Y5KMpxJINIziWAHlRHwDn3IrC5lV10XbZDBOcTdG1GxDn3iXPu7frlauAd4IAmu7X7+ox2\nco/kpadw++ilp+YiOZcOmFD/Z9pzZjYyZtElH12X0aVrswPqeyaOBZY02dTu67Pzr0E1FunT2aa/\n0fVUt7lIzslbwEHOue1mdgrwJHB414aV1HRdRo+uzXYysz7AX4GpLkx3ctp5fUb7zv0j4KC92gfh\nf8O0ts+B9euksTbPpXNuq3Nue/3y80CGme0XuxCTiq7LKNK12T5mlgE8BjzonHsyzC7tvj6jndy/\nfOnJzLrjX3p6usk+T+NHj2x4A/Y/zrmqKMeRDNo8l2Y2yMysfnk8vmvrptiHmhR0XUaRrs3I1Z+n\nu4HVzrnbWtit3ddnVMsyLoKXnpxzz5nZqWa2BtgGXBjNGJJFJOcSOBO4zMxqgO3AOYEFHOfMbB4w\nCRhgZuuBInwvJF2XHdDW+UTXZnt8HfghsLx+jC7wQ6kfDB2/PvUSk4hIEtI0eyIiSUjJXUQkCSm5\ni4gkISV3EZEkpOQuIpKElNxFRJKQkruISBJSchcRSUL/H2HtqHzXaAG7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10741c250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Trapezoidal integration\n",
    "\n",
    "def g(x):\n",
    "    return 1./(x**2+np.cos(x)**2)\n",
    "\n",
    "print \"Trapezoids: \", trapezoids(g, 0., np.pi, 1000000)\n",
    "\n",
    "# Simple Monte Carlo integration\n",
    "n0 = 100000\n",
    "a = np.arange(0.1,2.1,0.1)\n",
    "I = np.arange(0.1,2.1,0.1)\n",
    "\n",
    "r = np.random.random(n0)\n",
    "\n",
    "I0 = np.sum(1./((r*np.pi)**2+np.cos(r*np.pi)**2))\n",
    "print \"Simple Monte Carlo: \", I0/n0*np.pi\n",
    "\n",
    "# Importance Sampling\n",
    "print \"Importance Sampling:\"\n",
    "x = -np.log(r)\n",
    "i = 0\n",
    "for ai in a:\n",
    "    norm = (1.-np.exp(-ai*np.pi))/ai\n",
    "    x1 = norm*x/ai\n",
    "    Itot = 0.\n",
    "    Nin = 0\n",
    "    for xi in x1:\n",
    "        if(xi <= np.pi):\n",
    "            Nin += 1\n",
    "            Itot += g(xi)*np.exp(xi*ai)\n",
    "    Itot *= norm\n",
    "    I[i] = Itot/Nin\n",
    "    i += 1\n",
    "    print ai,Itot/Nin\n",
    "    \n",
    "pyplot.plot(a,I,ls='-',marker='o',c='red',lw=3);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 10.4: The Metropolis algorithm \n",
    "\n",
    "Use the Metropolis algorithm to sample points according to a ditribution\n",
    "and estimate the integral $$\\int _0^4 {x^2e^{-x}dx},$$ with\n",
    "$P(x)=e^{-x}$ for $0 \\leq x \\leq 4$. Plot the number of times the\n",
    "walker is at points $x_0$, $x_1$, $x_2$, ... Is the integrand sampled\n",
    "uniformly? If not, what is the approximate region of $x$ where the\n",
    "integrand is sampled more often?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trapezoids:  1.52379338694\n",
      "Metropolis:  1.58122841473\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEACAYAAACuzv3DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADvFJREFUeJzt3X+sZGddx/H3hy0FUbCamqrbNRip2poUiro2YmQIINfG\n0iDRuoqkaKR/uEr8Q2v5w97EqOlfEtKkViwNxoQlATStVFYEJhDELU36A+wu7oY02W2xgvyIIia7\n9usfM9lcZufOzN175s59zn2/kknnnHnuOc/cvf3c537P85xJVSFJatdzVt0BSdL2GOSS1DiDXJIa\nZ5BLUuMMcklqnEEuSY2bG+RJ3p3kmSSfndHmnUlOJnksyXXddlGSNMsiI/L7gLXNXkxyA/CSqroK\neCtwd0d9kyQtYG6QV9Unga/OaPJ64D3jtseAy5Jc0U33JEnzdFEj3w+c3rB9Briyg+NKkhbQ1cXO\nTGy77l+SdsglHRzjKeDAhu0rx/u+RRLDXZIuQlVNDpa/RRcj8vuBNwMkuR74WlU9s0lnevu44447\nVt4H35/vzffXv8ci5o7Ik7wXeCVweZLTwB3Ac8fBfE9VPZjkhiSngG8Ab1nozJKkTswN8qo6tECb\nw910R5K0Va7s7MhgMFh1F5aqz++vz+8NfH97QRatwWz7REnt1LkkqS+SUDtwsVOStEJdTD8Uo9+a\n0/hXiKRlM8g7NRnaM/8akqROWFqRpMYZ5JLUOINckhpnkEtS4wxySWqcQS5JjTPIJalxBrkkNc4g\nl6TGGeSS1DiDXJIaZ5BLUuMMcklqnEEuSY0zyCWpcQa5JDXOIJekxhnkktQ4g1ySGmeQS1LjDHJJ\natwlq+5A3yWZur+qdrgnkvrKIF+6aYE9Pdwl6WJYWpGkxhnkktQ4g1ySGmeQS1LjDHJJapxBLkmN\nM8glqXEGuSQ1bm6QJ1lLciLJySS3TXn98iQfTvJoks8luWUpPZUkTZVZS8WT7AM+D7wGeAr4DHCo\nqo5vaLMOPK+qbk9y+bj9FVV1buJY1edl6aOl+JPvb9q+0f4+fy8kdScJVTVzOfi8EflB4FRVPVlV\nZ4EjwE0Tbb4IvGj8/EXAf06GuCRpeebda2U/cHrD9hngpybavAv4WJKngRcCv9xd9/pr2s20HKVL\nuhjzgnyRZHk78GhVDZL8EPCRJC+tqv+abLi+vn7++WAwYDAYbKGrfTOtDCNprxsOhwyHwy19zbwa\n+fXAelWtjbdvB56tqjs3tHkQ+JOq+tR4+6PAbVX18MSxrJHP3G/dXNKFuqiRPwxcleTFSS4Fbgbu\nn2hzgtHFUJJcAfwI8IWL67Ikaatmllaq6lySw8BRYB9wb1UdT3Lr+PV7gD8F7kvyGKNfDH9QVV9Z\ncr8lSWMzSyudnsjSypz9llYkXaiL0ookaZczyCWpcQa5JDXOIJekxhnkktQ4g1ySGmeQS1LjDHJJ\napxBLkmNM8glqXEGuSQ1ziCXpMYZ5JLUOINckhpnkEtS4+Z9Zqd20LQPZAY/lFnSbAb5rrLZh1BI\n0uYsrUhS4wxySWqcQS5JjTPIJalxBrkkNc4gl6TGGeSS1DiDXJIaZ5BLUuMMcklqnEEuSY0zyCWp\ncd40qwHeFVHSLAZ5E7wroqTNWVqRpMYZ5JLUOINckhpnkEtS4+YGeZK1JCeSnExy2yZtBkkeSfK5\nJMPOeylJ2lRmTWFLsg/4PPAa4CngM8Chqjq+oc1lwKeA11XVmSSXV9WXpxyr+jxdbjRFcPL9Tdu3\n2f6ttB3t7/P3U9JIEqpq5jS1eSPyg8Cpqnqyqs4CR4CbJtr8KvCBqjoDMC3EJUnLMy/I9wOnN2yf\nGe/b6Crgu5N8PMnDSX69yw5KkmabtyBokb/dnwu8HHg18ALg00n+papObrdzkqT55gX5U8CBDdsH\nGI3KNzoNfLmqvgl8M8kngJcCFwT5+vr6+eeDwYDBYLD1Huu8aUv3rZtLbRsOhwyHwy19zbyLnZcw\nutj5auBp4CEuvNj5o8BdwOuA5wHHgJur6omJY3mxc+b+rV/snHaMPn+Ppb1okYudM0fkVXUuyWHg\nKLAPuLeqjie5dfz6PVV1IsmHgceBZ4F3TYa4JGl5Zo7IOz2RI/I5+x2RS7pQF9MPJUm7nEEuSY0z\nyCWpcQa5JDXOIJekxhnkktQ4g1ySGmeQS1LjDHJJapxBLkmNM8glqXEGuSQ1ziCXpMYZ5JLUuHmf\nEKTGTPvUIPCTg6Q+M8h7Z7N7l0vqK0srktQ4g1ySGmeQS1LjDHJJapxBLkmNM8glqXFOP9wjps0v\nd2651A8G+Z4xGdrOLZf6wtKKJDXOIJekxhnkktQ4g1ySGmeQS1LjDHJJapxBLkmNM8glqXEuCNrD\n/DQhqR8M8j3NTxOS+sDSiiQ1ziCXpMbNDfIka0lOJDmZ5LYZ7X4yybkkv9htFyVJs8wM8iT7gLuA\nNeAa4FCSqzdpdyfwYSyyStKOmjciPwicqqonq+oscAS4aUq73wHeD3yp4/5JkuaYF+T7gdMbts+M\n952XZD+jcL97vMu5a5K0g+YF+SKh/A7gD2s0+ThYWpGkHTVvHvlTwIEN2wcYjco3+nHgyHhxyeXA\nzyc5W1X3Tx5sfX39/PPBYMBgMNh6jyWpx4bDIcPhcEtfk1mr+JJcAnweeDXwNPAQcKiqjm/S/j7g\ngar64JTXqs8rBke/yKZ9nNpmi26207aLY2zets//TlJrklBVMysdM0fkVXUuyWHgKLAPuLeqjie5\ndfz6PZ31VruGS/eltswckXd6Ikfkc/bvnhG5I3Vp91hkRO7KTklqnEEuSY0zyCWpcd7GVgubdhHU\nurm0ega5tmDaBVNJq2ZpRZIaZ5BLUuMMcklqnEEuSY0zyCWpcQa5JDXOIJekxhnkktQ4FwRpW7zl\nrbR6Brm2abNb4UraKZZWJKlxBrkkNc4gl6TGGeSS1DgvdmopvHe5tHMMci2J9y6XdoqlFUlqnEEu\nSY0zyCWpcQa5JDXOIJekxjlrRTvGG2xJy2GQawd5gy1pGSytSFLjDHJJapylFa2cy/ml7THItQu4\nnF/aDksrktQ4g1ySGmeQS1LjrJFrV3LxkLS4hUbkSdaSnEhyMsltU17/tSSPJXk8yaeSXNt9V7W3\n1JSHpGnmBnmSfcBdwBpwDXAoydUTzb4A/GxVXQv8MfCXXXdUkjTdIiPyg8Cpqnqyqs4CR4CbNjao\nqk9X1dfHm8eAK7vtpiRpM4sE+X7g9IbtM+N9m/lN4MHtdEqStLhFLnYuXJxM8irgN4BXTHt9fX39\n/PPBYMBgMFj00JK0JwyHQ4bD4Za+JvNmASS5HlivqrXx9u3As1V150S7a4EPAmtVdWrKcarPMw5G\nsyymrVDc7I5/22nbxTF28/k2b9vnnyFpmiRU1czlzouUVh4Grkry4iSXAjcD90+c6AcYhfibpoW4\n1JUkUx/SXja3tFJV55IcBo4C+4B7q+p4klvHr98D/BHwXcDd4/+pzlbVweV1W3uX9zSXJs0trXR2\nIksrc/bv5lLHTp9v633r88+W9rZFSiuu7FQveCtc7WUGuXrCW+Fq7/KmWZLUOINckhpnaUW95R0U\ntVcY5Ooxpypqb7C0IkmNM8glqXEGuSQ1zhq59hwXD6lvDHLtQS4eUr8Y5BJOVVTbDHIJcKqiWubF\nTklqnCNyaQZLLmqBQS7NZMlFu59BLl0EpzBqNzHIpYty4RRGyzBaFYNc6oxlGK2Gs1YkqXEGuSQ1\nztKKtGReGNWyGeTS0nlhVMtlkEsrMf3CqKN3XQyDXNpVvDOjts4gl3Y5yzCaxyCXdr3FyzBgwO9F\nBrnULBcgacQgl3pms5H6NI7e+8Egl3pn2gVTR+99ZpBLe5jTHfvBIJf2NBcr9YFBLmnC1mbJTD2C\nob+jDHJJC7L2vlsZ5JI65+h9Z829jW2StSQnkpxMctsmbd45fv2xJNd1301Jbakpj2n7R6G/3cde\nNzPIk+wD7gLWgGuAQ0munmhzA/CSqroKeCtw95L6KqmXFg/9zdvubfNG5AeBU1X1ZFWdBY4AN020\neT3wHoCqOgZcluSKznsqSZpqXpDvB05v2D4z3jevzZXb75okaRHzLnYu+nfLZJFq6tfdeOON37L9\nxje+kVtuuWXBU0iSpsmsK8ZJrgfWq2ptvH078GxV3bmhzV8Aw6o6Mt4+Abyyqp6ZOJbFLEm6CFU1\n84ruvBH5w8BVSV4MPA3cDByaaHM/cBg4Mg7+r02G+CIdkSRdnJlBXlXnkhwGjgL7gHur6niSW8ev\n31NVDya5Ickp4BvAW5bea0nSeTNLK5Kk3W/ugqAuJfmlJP+a5P+SvHwnz70siyyYalWSdyd5Jsln\nV92XZUhyIMnHxz+Tn0vyu6vuU5eSPD/JsSSPJnkiyZ+tuk9dS7IvySNJHlh1X7qW5Mkkj4/f30Oz\n2u5okAOfBd4AfGKHz7sUiyyYatx9jN5bX50Ffq+qfgy4HvjtPv37VdX/Aq+qqpcB1wKvSvIzK+5W\n194GPEE/VwYVMKiq66rq4KyGOxrkVXWiqv5tJ8+5ZIssmGpWVX0S+Oqq+7EsVfXvVfXo+Pl/A8eB\n719tr7pVVf8zfnopo+tcX1lhdzqV5ErgBuCv6O+duhZ6Xzs9Iu+bRRZMqQHjmVnXAcdW25NuJXlO\nkkeBZ4CPV9UTq+5Th/4c+H3g2VV3ZEkK+KckDyf5rVkNO7/7YZKPAN875aW3V1Xf6lh9/HNuz0ny\nHcD7gbeNR+a9UVXPAi9L8p3A0SSDqhquuFvbluQXgP+oqkeSDFbdnyV5RVV9Mcn3AB9JcmL8V/IF\nOg/yqnpt18fcxZ4CDmzYPsBoVK5GJHku8AHgb6rq71bdn2Wpqq8n+RDwE8Bwxd3pwk8Drx/ftO/5\nwIuS/HVVvXnF/epMVX1x/N8vJflbRqXcqUG+ytJKH2pa5xdMJbmU0YKp+1fcJy0oo/uf3gs8UVXv\nWHV/upbk8iSXjZ9/G/Ba4JHV9qobVfX2qjpQVT8I/ArwsT6FeJIXJHnh+Pm3Az/HaLLIVDs9/fAN\nSU4zmiHwoST/sJPn71pVnWO0qvUooyvn76uq46vtVXeSvBf4Z+CHk5xO0rfFXq8A3sRoNscj40ef\nZul8H/CxcY38GPBAVX10xX1alr6VOa8APrnh3+7vq+ofN2vsgiBJapyzViSpcQa5JDXOIJekxhnk\nktQ4g1ySGmeQS1LjDHJJapxBLkmN+3/RNEP9/eirFgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1072f74d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "delta = 2\n",
    "xmin = 0.\n",
    "xmax = 4.\n",
    "\n",
    "def f(x):\n",
    "    return x**2*np.exp(-x)\n",
    "\n",
    "def P(x):\n",
    "    global xmin, xmax\n",
    "    if(x < xmin or x > xmax):\n",
    "        return 0.\n",
    "    return np.exp(-x)\n",
    "\n",
    "\n",
    "def metropolis(xold):\n",
    "    global delta\n",
    "    xtrial = np.random.random()    \n",
    "    xtrial = xold+(2*xtrial-1)*delta\n",
    "    weight = P(xtrial)/P(xold)\n",
    "    xnew = xold\n",
    "    if(weight >= 1): #Accept\n",
    "        xnew = xtrial\n",
    "    elif(weight != 0):\n",
    "        r = np.random.random()\n",
    "        if(r <= weight): #Accept\n",
    "            xnew = xtrial\n",
    "    return xnew\n",
    "\n",
    "\n",
    "xwalker = (xmax+xmin)/2.\n",
    "for i in range(100000):\n",
    "    xwalker = metropolis(xwalker)\n",
    "\n",
    "I0 = 0.\n",
    "N = 300000\n",
    "\n",
    "x = np.zeros(N)\n",
    "x[0] = xwalker\n",
    "for i in range(1,N):\n",
    "    x[i] = metropolis(x[i-1])\n",
    "    I0 += x[i]**2\n",
    "\n",
    "binwidth=0.1\n",
    "pyplot.hist(x,bins=np.arange(xmin-1, xmax+1, 0.1),normed=True);\n",
    "\n",
    "print \"Trapezoids: \", trapezoids(f,xmin,xmax,10000)\n",
    "print \"Metropolis: \", I0/(1.-np.exp(-4.))/N"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Challenge 10.1\n",
    "\n",
    "- Calculate the integral $\\int_0^1 x^2 dx=1/3$ using simple MC integration and importance sampling with $P(x)=x$.\n",
    "\n",
    "- Calculate the integral $\\int_0^1 \\sqrt{x}dx=2/3$ using simple MC integration and $P(x)=1-e^{-ax}$. Find the values of $a$ that minimizes the variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}