Lec_14.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Physics 300 \n",
    "## Computational Physics I (Fall 2017)\n",
    "## BPB-248, Tues/Thurs 10:00 - 11:15 am \n",
    "\n",
    "|Instructor| Prof. Qiang Zhu|\n",
    "|--|-------------------------------|\n",
    "|Email | qiang.zhu@unlv.edu|\n",
    "|Website|http://www.physics.unlv.edu/~qzhu/|\n",
    "|Office| BPB 232|\n",
    "|Office hours | Tues/Thurs 8:30 - 10:00 |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 14 Monte Carlo\n",
    "\n",
    "We have introduced the Monte Carlo methods in the last lecture. It is not a very good one for doing integrals.\n",
    "There are significantly better ways to perform the calculation, the most common of which is the _mean value method_.\n",
    "\n",
    "## 14.1 The Mean value method\n",
    "\n",
    "Let's take a general integration problem: \n",
    "\n",
    "$$ I = \\int_a^b f(x)dx $$\n",
    "\n",
    "The average value of $f(x)$ in the range from $a$ to $b$ is by definition\n",
    "\n",
    "$$ <f> = \\frac{1}{b-a}\\int_a^b f(x)dx = \\frac{I}{b-a} $$\n",
    "\n",
    "Thus\n",
    "\n",
    "$$ I = (b-a)<f> $$\n",
    "\n",
    "If we can estimate <f> then we can estimate $I$. A simple way to estimate <f> is just to measure $f(x)$\n",
    "at $N$ points $x_1, x_2 .... x_N$ chosen uniformly at random between $a$ and $b$, and then calculate\n",
    "$<f> \\approx N^{-1} \\sum^N_{i=1} f(x_i)$,  \n",
    "\n",
    "$$ I \\approx \\frac{b-a}{N} \\sum^N_{i=1} f(x_i)$$\n",
    "\n",
    "This is the fundamental formula for the mean value method. As with the previous method, it gives only an approximate\n",
    "estimate of the integral. How accurate is it? We can estimate the error using standard results for the behavior\n",
    "of random variable which tell us that the variance on the sum of $N$ independent random numbers - no matter what\n",
    "their distribution - is equal to $N$ times the variance on a single one.\n",
    "\n",
    "$$ f = <f^2> - <f>^2 $$\n",
    "in which\n",
    "$$ <f> = \\frac{1}{N} \\sum_{i=1}^N f(x_i)$$\n",
    "$$ <f^2> = \\frac{1}{N} \\sum_{i=1}^N [f(x_i)]^2$$\n",
    "\n",
    "Then the variance on the sum is $N$ times the variance on a single term. Finally, we get the standard deviation on $I$,\n",
    "\n",
    "$$ \\sigma = \\frac{b-a}{N} \\frac{\\sqrt{\\text{var}f}} {\\sqrt{N}} $$\n",
    "\n",
    "Thus, the errors goes as 1/$\\sqrt{N}$. However, the leading constant is smaller, which means the method is always\n",
    "more accurate than the previous method.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MC intergral\n",
    "import matplotlib.pyplot as plt\n",
    "from math import sin\n",
    "from random import random\n",
    "\n",
    "def f(x):\n",
    "    return sin(1/x/(2-x))**2\n",
    "\n",
    "def MC1(N):\n",
    "    count = 0\n",
    "    for i in range(N):\n",
    "        x = 2*random()\n",
    "        y = random()\n",
    "        if y<f(x):\n",
    "            count += 1\n",
    "    I = 2*count/N\n",
    "    return I\n",
    "\n",
    "\n",
    "def MC2(N):\n",
    "    I = 0\n",
    "    for i in range(N):\n",
    "        x = 2*random()\n",
    "        I += f(x)\n",
    "\n",
    "    return I/N*2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.514\n",
      "1.4552373349084156\n"
     ]
    }
   ],
   "source": [
    "print(MC1(1000))\n",
    "print(MC2(1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 14.2 Integrals in high dimensions\n",
    "\n",
    "In addition to the intergration of pathological functions, Monte Carlo integration is used for performing\n",
    "high dimensional integrals. As we saw in the previous section, performing an integral over two variables\n",
    "by standard grid of points such as trapezoidal rule requires us to take samples on a two-dimensional grid.\n",
    "If we have an integral over four dimensions or even higher, then the number of points on the grid can \n",
    "become very slow and impractical. Instead, Monte Carlo in high dimension is an obvious extension of the \n",
    "technique in one dimension. The integral of a function $f(r)$ over a volume $V$ in a high-dimensional space\n",
    "is given by the generalization:\n",
    "\n",
    "$$ I = \\frac{V}{N} \\sum_{i=1}^{N}f(r_i)$$\n",
    "\n",
    "where the points $r$ are picked uniformly at random from the volume $V$.\n",
    "\n",
    "An important application of this type of integral is in financial mathematics. This is widely used for \n",
    "predicting the values of stocks, which requires the evaluation of integrals over many variables. Such\n",
    "information could be used to instruct the quick trading decisions. This is why stastical physicists\n",
    "were popular in Wall Street!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 14.3 Importance sampling\n",
    "\n",
    "Monte Carlo is useful for integrating pathological functions. However, it is still not good for some cases. \n",
    "In particular, it gives problems if the function contains a divergence. For example,\n",
    "\n",
    "$$ I = \\int_0 ^1 \\frac{x^{-1/2}}{e^x + 1} dx $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/qiangzhu/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:9: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  if __name__ == '__main__':\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGuhJREFUeJzt3Xt4XPV95/H3V3PVzOhmSbYly3dsjLkbl0so4ZItoaS5\ntLRZkidtmk3CQvNs0m03abt9+myzfXaf3T+abpPSzZImTXMHQkNoArtpAsRAjB0bbMAYbPAF27It\nWffRSDMa6bd/zMhXyRoZzZxzNJ/X8+g5Z86cGX1/Hunjn37nd84x5xwiIhIcNV4XICIis6PgFhEJ\nGAW3iEjAKLhFRAJGwS0iEjAKbhGRgFFwi4gEjIJbRCRgFNwiIgETLsebtrS0uBUrVpTjrUVE5qXt\n27efcM61lrJvWYJ7xYoVbNu2rRxvLSIyL5nZwVL31VCJiEjAKLhFRAKmpKESMzsADAHjQN45t7Gc\nRYmIyPRmM8Z9q3PuRNkqERGRkmioREQkYEoNbgf8xMy2m9k95SxIRETOr9Shkl91zh0xs4XAv5rZ\na865TafvUAz0ewCWLVs2x2WKiMikknrczrkjxWUX8APg2in2ecA5t9E5t7G1taQ55Of44s/28vM9\n3Rf0WhGRajFjcJtZ0szqJteB24FXylHMA5v2sUnBLSJyXqUMlSwCfmBmk/t/xzn3f8tRTCIaYjib\nL8dbi4jMGzMGt3NuH3BlBWohFQuTVnCLiJyXr6YDJmNh9bhFRGbgs+AOMZwd97oMERFf81Vwa6hE\nRGRmvgruRDRMJqfgFhE5H18FdzIWJq2hEhGR8/JVcKdimg4oIjITXwV3MhZmZGyc8QnndSkiIr7l\nq+BOxQrTyoc1zi0iMi1fBXdyMrg1XCIiMi1fBXciGgLQXG4RkfPwVXCn1OMWEZmRr4JbQyUiIjPz\nVXBP9rh19qSIyPR8FdxJzSoREZmRz4K7cHBSZ0+KiEzPX8EdLfS4MxoqERGZlq+COxENYaaDkyIi\n5+Or4DYzklFdaEpE5Hx8FdwweTMF9bhFRKbjw+AOk9asEhGRafkuuFO676SIyHn5LrgT0RAZjXGL\niEzLd8Gt+06KiJyf74I7GQvrzEkRkfPwZ3Crxy0iMi3fBbeGSkREzs93wZ2MhhkdmyA/PuF1KSIi\nvuS/4C5eaCozppklIiJT8WFw62YKIiLno+AWEQkY3wV3StfkFhE5L98F9+Q1udXjFhGZmv+CW/ed\nFBE5L98Gt3rcIiJTKzm4zSxkZi+a2Y/KWdDkdMDhnMa4RUSmMpse92eA3eUqZFJKPW4RkfMqKbjN\nrAN4D/AP5S0HaiMhanTfSRGRaZXa4/5fwOeAac9DN7N7zGybmW3r7u6+4IJO3XdSwS0iMpUZg9vM\nfgPocs5tP99+zrkHnHMbnXMbW1tb31ZRukKgiMj0Sulx3wi8z8wOAN8DbjOzb5WzqEQsxLBOwBER\nmdKMwe2c+zPnXIdzbgVwN/Ckc+4j5SwqpZspiIhMy3fzuKFw9qSGSkREpjar4HbOPe2c+41yFTMp\nGQvrWiUiItPwZY87FQupxy0iMg1fBndCs0pERKbly+DWfSdFRKbny+BORsNk87rvpIjIVPwZ3LrQ\nlIjItHwZ3LrQlIjI9HwZ3Lomt4jI9Hwa3JP3nVRwi4iczZ/BffK+kxrjFhE5mz+De3KoRNcrERE5\nhy+DWwcnRUSm58vg1sFJEZHp+TK4J3vcutCUiMi5fBnc8UiN7jspIjINXwa37jspIjI9XwY3FMa5\nM5pVIiJyDt8Gd0NthL7MmNdliIj4jm+De3FDnKMDI16XISLiO74N7vbGOMcGRr0uQ0TEd3wb3G0N\ntZxI58jmNSVQROR0Pg7uOIB63SIiZ/FxcNcC0Nmv4BYROZ1/g7ux0OPWAUoRkTP5Nrjbiz3uoxoq\nERE5g2+DuzYaojERUY9bROQsvg1uKIxzH9UYt4jIGXwe3HE6NVQiInIG3wf3MQ2ViIicwdfB3d5Y\nS19mjJGcTsIREZnk6+CePAlHByhFRE7xeXBrSqCIyNl8HtyFHndnv3rcIiKTfB3ci3W9EhGRc8wY\n3GYWN7OtZrbTzHaZ2ecrURhAPBKiORnVlEARkdOES9gnC9zmnEubWQR41syecM49X+bagMI1S3Rw\nUkTklBl73K4gXXwYKX65slZ1Gp09KSJyppLGuM0sZGY7gC7gX51zW8pb1imFsyfV4xYRmVRScDvn\nxp1zVwEdwLVmdtnZ+5jZPWa2zcy2dXd3z1mBbQ21DI3mSWd1x3cREZjlrBLnXD/wFHDHFM894Jzb\n6Jzb2NraOlf10d44ObNEvW4REShtVkmrmTUW12uBXwNeK3dhk3QnHBGRM5Uyq6QN+CczC1EI+oec\ncz8qb1mnfXOd9i4icoYZg9s59xJwdQVqmdKi+jhmOu1dRGSSr8+cBIiGa2hJxTQlUESkyPfBDdCu\nKYEiIicFIrhXt6bYc3zI6zJERHwhEMG9vr2e44NZTqSzXpciIuK5YAR3Wz0Au48OelyJiIj3AhHc\nlxSD+9VOBbeISCCCuykZpb0hzqvqcYuIBCO4Ada3N7BLPW4RkSAFdz37utO647uIVL3gBHdbPRMO\nXte0QBGpcoEJ7kvbdYBSRAQCFNwdTbXUxcK8enTA61JERDwVmOA2My5pr1ePW0SqXmCCGwrj3K8d\nG2J8omK3vBQR8Z1gBXd7PZncOAd7hr0uRUTEM8EK7uIZlJrPLSLVLFDBvWZRinCN6QxKEalqgQru\nWDjERQtTOkApIlUtUMENcPmSBnYe7tcBShGpWoEL7hsvaqE/M8auTs3nFpHqFMjgBnhm7wmPKxER\n8Ubggru1Lsb6tno27en2uhQREU8ELrgBblrbwgtv9ZHO5r0uRUSk4gIZ3DevaWVs3LFlX4/XpYiI\nVFwgg/uaFU3EIzUa5xaRqhTI4I6FQ1y/qplNezXOLSLVJ5DBDXDTmlb2dQ9zuC/jdSkiIhUV2OB+\n55rCtMBnNVwiIlUmsMF90cIUi+vjGucWkaoT2OA2M25a08Ize7vJ5Se8LkdEpGICG9wAd17exuBo\nnqde7/K6FBGRigl0cN+0poWWVIxHth/2uhQRkYoJdHCHQzV84Kp2nnq9i97hnNfliIhUxIzBbWZL\nzewpM3vVzHaZ2WcqUVip7rqmg7Fxx7/s7PS6FBGRiiilx50H/tg5tx64HviUma0vb1mlu6Stnkva\n6vnnFzRcIiLVYcbgds4ddc69UFwfAnYDS8pd2GzctWEJOw8P8EbXkNeliIiU3azGuM1sBXA1sKUc\nxVyo91+1hFCN8f3tR7wuRUSk7EoObjNLAY8Af+icO+emj2Z2j5ltM7Nt3d2VvYZIa12Mm9e28uiL\nR3RLMxGZ90oKbjOLUAjtbzvn/nmqfZxzDzjnNjrnNra2ts5ljSX54MalHBsc5fGXj1b8e4uIVFIp\ns0oM+Cqw2zn3hfKXdGFuX7+IVa1J/v7pN3FOvW4Rmb9K6XHfCPwucJuZ7Sh+3Vnmumatpsa49+bV\n7D46yNO6rZmIzGOlzCp51jlnzrkrnHNXFb8er0Rxs/WBq5bQ1hDnfz/1pteliIiUTaDPnDxbNFzD\nJ29axdYDvWw70Ot1OSIiZTGvghvg7muX0pSI8PdPq9ctIvPTvAvuRDTMx25cyZOvdbHzUL/X5YiI\nzLl5F9wAH7txBS2pGJ//l11MaF63iMwz8zK46+IRPnfHxbzwVj+P7tDZlCIyv8zL4Ab47Q0dXNnR\nwP944jXS2bzX5YiIzJl5G9w1NcZfvu9Suoay3P/UG16XIyIyZ+ZtcANcvayJuzZ08NVn9vNmd9rr\nckRE5sS8Dm6AP/n1i0nEQvzRQzsZG9dNhUUk+OZ9cC+si/Pff/Nydh7q50tPashERIJv3gc3FO4G\n/1sblvB3T+5l+8E+r8sREXlbqiK4AT7/vktpb6zlPz64Q7NMRCTQqia46+IRvvDBqzjcl+GzD+/U\niTkiElhVE9wA165cwJ/9+iU88cox/vZne70uR0TkgoS9LqDSPnHTSl4/PsTf/mwvaxfV8Z4r2rwu\nSURkVqqqxw1gZvy337yMDcsa+eOHd/Dy4QGvSxIRmZWqC26AWDjEl3/3GpqTMT76j1t5o2vI65JE\nREpWlcENhfnd3/rEdYRqjA9/ZQsHe4a9LklEpCRVG9wAK1uSfOvj1zE2PsGHv7KFzv4Rr0sSEZlR\nVQc3wMWL6/jGv7uOwZExfufLm9mna5qIiM9VfXADXN7RwHc+eT2jY+P8zpc364CliPiagrvo8o4G\nHr73BuKREB/6yvM898YJr0sSEZmSgvs0q1pTPHLfO2hvjPPRr23lm5sP4JzOsBQRf1Fwn2VxQ5xH\n7nsHN69t5S9+uIv//IOXyeV1OVgR8Q8F9xTq4hEe+L2NfOrW1Xx36yHufmAzh3ozXpclIgIouKcV\nqjE+++513P/hDew9nubOLz7D4y8f9bosEREF90zec0UbP/70TaxuTfEH336BP33kJYZGx7wuS0Sq\nmIK7BMuaEzx87w3cd8tqHtp2iHf/zSZ+vqfb67JEpEopuEsUCdXwJ3es45H73kEiFuajX9vKf3p4\nJz3prNeliUiVUXDP0tXLmvjxp3+VT926mkdfPMJtf/1zvrn5AOO6MYOIVIiC+wLEwiE+++51PPGZ\nm7i0vZ6/+OEu3vulZ/mFTtoRkQpQcL8NaxbV8e1PXMf9H97AwMgYH/6HLfz+P27ltWODXpcmIvOY\nlePMwI0bN7pt27bN+fv62ejYON/cfJAvPbmXoWye917RzqfftYaLFqa8Lk1EAsDMtjvnNpay74w9\nbjP7mpl1mdkrb7+0+SseCfHJd65i0+du5d+/czU/3X2c2//m5/zh915kz3HdqEFE5s6MPW4zeyeQ\nBr7hnLuslDetxh732XrSWb7yzH6+sfkAmdw471q3kHtvWc3G5U2YmdfliYjPzKbHXdJQiZmtAH6k\n4J69vuEc39h8kK//Yj99mTGu7GjgYzeu5M7L24iGdYhBRAoU3D40khvn4e2H+PpzB9h3YpjWuhgf\n+pWl/Ntrl7Gksdbr8kTEY54Et5ndA9wDsGzZsmsOHjxYUrHVZmLCsWlvN1//xQF+vqcbA265eCF3\n/8pSbl23kEhIvXCRaqQed0Ac6s3w4C8P8eC2Q3QPZWlORvnA1Uv4rQ1LWN9Wr7FwkSqi4A6YsfEJ\nNu3p5vvbD/PT3ccZG3dctDDFB65q571XtrO8Oel1iSJSZnMa3Gb2XeAWoAU4DvwX59xXz/caBfeF\n6xvO8fgrR/nhi51sPdALwKXt9dx5eRt3Xt7GyhaFuMh8NOc97tlScM+NI/0jPP7SUX788lF2HOoH\nYO2iFLevX8ztly7isvYGamo0nCIyHyi456Ej/SP8v1eO8ZNXj7F1fy8TDlrrYtx28UJuXbeQGy9q\npi4e8bpMEblACu55rm84x5OvdfHk611s2tPN0GiecI2xYXkTN69t5aY1LVza3kBIvXGRwFBwV5Gx\n8Qm2Hehj095uNu3pZldn4QJXjYkIN6xq5h2rm7lhdTOrW1OapSLiYwruKtY9lOUXb57g2b0neO6N\nE3QOjAKFYZVrVy7gupULuHblAtYurNP4uIiPKLgFAOccb/Vm2PxmD5v39bBlXy/HBgtBXh8Ps2F5\nExuXN7FheRNXdjSSjIU9rlikes0muPWbOo+ZGcubkyxvTnL3tctwznG4b4Qt+3vZfrCXbQf6ePr1\nwr0zawwuXlzP1csauaqjkauWNbK6NaVxchEfUo+7yvVncrx4qJ8X3+rnxbf62HGon6HRPACJaIjL\n2hu4bEkDV3Q0cNmSela2KMxFykFDJXLBJiYc+3uG2Xmon52H+nn5yAC7OgfJ5icAqI2EuKStjkvb\nG7ikrZ717fVcvKiO2mjI48pFgk3BLXMqPz7BG91pdh0Z5JXOAXYdGWT30UGGsoWeuRmsbE6yrq2O\ntYvqWLe4sFzenFTvXKREGuOWORUO1bBucT3rFtdz1zUdACfHy3d1DvDasSF2Hx1kV+cgT7xyjMm+\nQCxcw+rWFGsWpVi7qO7k+vIFCcK6CqLIBVNwywUxM5YuSLB0QYI7Lms7uT2Ty7P3eJrXjw2xt2uI\nPcfT/HJ/Lz/c0Xlyn0iocNB0dWuSVa0pVrWcWjYlo140RyRQFNwypxLRMFcubeTKpY1nbE9n87zZ\nlWZvV5p93Wne7E7zRleaJ1/rYmz81HBdQ22ElS1JVrYkWdGcZEVLojAzZkGCxkREJxGJoOCWCknF\npg70/PgEh/tG2Hcizb7uYQ70DLP/xDBb9vXw6I4jnH4Ipi4eZnlzguULkixdkGDZggRLF9SybEGC\n9sZa3YRCqoaCWzwVDtWwoiXJipYkt60787nRsXEO9WY40JPhYM8wb/VmONiTYffRQX7y6rEzeuo1\nBovr43Q0JehoqqWjqZYlTbV0NCVY0ljL4oY48Yhmvsj8oOAW34pHQqxZVMeaRXXnPDc+4Tg2OMqh\n3kzhq2+Ew30ZDveOsHlfD8cHR5k4a8JUSyrGksY47Y21tDfW0tYQP2PZkoppFowEgoJbAilUYyxp\nrGVJYy3Xr2o+5/mx8QmODYxyuG+EI/0jdPaPcKRvhM6BEfYcH+Lp17sZGRs/5z0X1cVY3BBncUOc\nRfVx2orLRfVxFheXmrMuXlNwy7wUCdWcnPUyFeccAyNjdPaPcmxwpLAcGKVzYIRjA6O8dqwQ7pnc\n+DmvrY+HWVQfZ2F9jIV1py3rYrTWxVhYF2NhfZxkNKSDqVIWCm6pSmZGYyJKYyLK+vb6KfdxzjGU\nzdM1OMqxgSzHBkc5PjhK1+AoxwezdA2NsnV/L91DWXLjE+e8vjYSorUY5i2pKC2pyfVYcb2wrTkV\nU8jLrCi4RaZhZtTHI9THI1y08Nxx9knOOfozY3QNFcK8eyhL91CWrqEsJ9KF9X3dw2zd30tfZmzK\n94hHamhOFgK+ORVjQTJKcypKczLKgmSsuIye3J6I6le3munTF3mbzIymZJSmZJSLF08f8FAYe+9J\n5ziRLoT65HrP8OS2HF1Do+w+OkhPOjdlTx4KQb8gEWVBKkpTohDoTYnJ9UihnkSUxkTk5HOaVTN/\nKLhFKigSqjl58HMmzjnS2Tw96Rw9wzl6h3P0DhdCvm84R+/wGL3DWfoyY7zVm6E3nTt5/ZipxCM1\nNNYWwnwy1Bsnl7URGhMRGmqjNJxcLyxrIxrG8RsFt4hPmRl18Qh18QgrWpIlvSaXn6B/JEff8Bh9\nmRz9mULA92VyDIyM0Tecoy8zxsBIjr1dafqL20+fE3+2aKiG+tow9bWFMG+oLQR9Q23k5Lb6eGG9\nvjZMffzUtlQ8rCmWZaDgFplHouGa4gyXmXv0k5xzZHLjJ8N9IDNG/8hYYX1kjP5MYTlYfNyTzrGv\ne7iwbXSMmS4wWhcLUxcvBH9dvBDspz+ui5+1jJ2+LUwyGtZt9s6i4BapcmZGMhYmGQvT0TS7105M\nONK5PAOZMYZG8wyOngr5wdH8ybAfGs0zVHzu6MAoe7omt+UZP/tMqXPqg1Q0TCoeJhU7tSyEfOTU\n9uJzyVgh/AttClEXi5CMhUjGwsTCNfNi2EfBLSIXrKbm1MybCzHZ258M9sHiMp3Nky4G+1C2sG1o\nNM9wNk86m2dwNE9n/wjD2XGGRscYnmK+/VTCNVYI92gh6BOxEKlY4fFk0CdjYZLRyWVheyIWIhkN\nk4ieej4RC1MbCXkyFKTgFhHPnN7bL+WA7XQmJhzDufzJwE9nz1wfzuYZzo2fXE9n82Sy4ydfc3xw\nlOHsqefzM/wVcLp4pKYQ6rEQbfW1PHTvDRfcjlIpuEUk8GpqTh3IpeHtv18uP8FwNk9mbPxk0I/k\nxov/AeQZzo4XHucKQZ/JjZPJjROPVOYKlQpuEZGzRMM1RMNRZjnkXzG6gLGISMAouEVEAkbBLSIS\nMApuEZGAUXCLiASMgltEJGAU3CIiAaPgFhEJGHMzXdrrQt7UrBs4OIuXtAAn5rwQf6vGNkN1trsa\n2wzV2e630+blzrnWUnYsS3DPlpltc85t9LqOSqrGNkN1trsa2wzV2e5KtVlDJSIiAaPgFhEJGL8E\n9wNeF+CBamwzVGe7q7HNUJ3trkibfTHGLSIipfNLj1tEREpUseA2szvM7HUze8PM/nSK52Nm9mDx\n+S1mtqJStZVTCe3+IzN71cxeMrOfmdlyL+qcSzO1+bT97jIzZ2bzYuZBKe02sw8WP+9dZvadStc4\n10r4+V5mZk+Z2YvFn/E7vahzLpnZ18ysy8xemeZ5M7MvFv9NXjKzDXNehHOu7F9ACHgTWAVEgZ3A\n+rP2+QPgy8X1u4EHK1GbD9p9K5Aort8X9HaX0ubifnXAJuB5YKPXdVfos14DvAg0FR8v9LruCrT5\nAeC+4vp64IDXdc9Bu98JbABemeb5O4EnAAOuB7bMdQ2V6nFfC7zhnNvnnMsB3wPef9Y+7wf+qbj+\nfeBdFvzbMc/YbufcU865TPHh80BHhWuca6V81gB/BfxPYLSSxZVRKe3+JHC/c64PwDnXVeEa51op\nbXZAfXG9AeisYH1l4ZzbBPSeZ5f3A99wBc8DjWbWNpc1VCq4lwCHTnt8uLhtyn2cc3lgAGiuSHXl\nU0q7T/dxCv9TB9mMbS7+6bjUOffjShZWZqV81muBtWb2nJk9b2Z3VKy68iilzX8JfMTMDgOPA/+h\nMqV5ara/97Ome076hJl9BNgI3Ox1LeVkZjXAF4Df97gUL4QpDJfcQuEvq01mdrlzrt/TqsrrQ8DX\nnXN/bWY3AN80s8uccxNeFxZklepxHwGWnva4o7htyn3MLEzhz6qeilRXPqW0GzP7N8CfA+9zzmUr\nVFu5zNTmOuAy4GkzO0BhDPCxeXCAspTP+jDwmHNuzDm3H9hDIciDqpQ2fxx4CMA5txmIU7iex3xW\n0u/921Gp4P4lsMbMVppZlMLBx8fO2ucx4KPF9d8GnnTFkf4Am7HdZnY18H8ohHbQxzxhhjY75wac\ncy3OuRXOuRUUxvXf55zb5k25c6aUn/FHKfS2MbMWCkMn+ypZ5Bwrpc1vAe8CMLNLKAR3d0WrrLzH\ngN8rzi65Hhhwzh2d0+9QwSOxd1LoYbwJ/Hlx23+l8EsLhQ/0YeANYCuwyuujxxVq90+B48CO4tdj\nXtdc7jafte/TzINZJSV+1kZhmOhV4GXgbq9rrkCb1wPPUZhxsgO43eua56DN3wWOAmMU/or6OHAv\ncO9pn/P9xX+Tl8vx860zJ0VEAkZnToqIBIyCW0QkYBTcIiIBo+AWEQkYBbeISMAouEVEAkbBLSIS\nMApuEZGA+f/yhepc9NFKKgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1085976d8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Divergence problem\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from math import pi\n",
    "\n",
    "\n",
    "def f(x):\n",
    "    return 1/np.sqrt(x)/(np.exp(x)+1)\n",
    "\n",
    "x = np.linspace(0,1,100)\n",
    "plt.plot(x, f(x))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you try to do it using the mean value method, you will run into problems becuase the value of f(x) diverges when x is approaching to 0. This means that the estimated integral could go up and down without convergence.\n",
    "\n",
    "We can get around these problems by drawing our points x _non-uniformly_. This technique is called _importance sampling_. For any general functions $g(x)$ we can define a weighted average over the interval from $a$ to $b$ thus:\n",
    "\n",
    "$$ \\langle g \\rangle_w = \\frac{\\int_a^b w(x)g(x)dx}{\\int_a^b w(x)dx} $$\n",
    "\n",
    "where $w(x)$ is any function we choose. Now consider the general 1d integral:\n",
    "\n",
    "$$ I = \\int_a^b f(x)dx $$\n",
    "\n",
    "Setting g(x)=f(x)/w(x), we have\n",
    "$$ \n",
    "\\begin{aligned}\n",
    "\\bigg \\langle \\frac{f(x)}{w(x)} \\bigg \\rangle _w= & \\frac{\\int_a^b w(x)f(x)/w(x)dx}{\\int_a^b w(x)dx} \\\\\n",
    "                      = & \\frac{\\int_a^b f(x)dx}{\\int_a^b w(x)dx}\\\\\n",
    "                      = & \\frac{I}{\\int_a^b w(x)dx}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "So\n",
    "$$ I = \\bigg \\langle \\frac{f(x)}{w(x)} \\bigg \\rangle _w \\int_a^b w(x)dx$$\n",
    "\n",
    "This formula is similar to what we have for the mean value method. But it allows us to calculate the value\n",
    "of the integral from a weighted average, rather than a standard uniform average.\n",
    "\n",
    "To calculate it with the weighted average, let's define a probability function\n",
    "\n",
    "$$ p(x) = \\frac{w(x)}{\\int_a^b w(x)dx} $$\n",
    "\n",
    "The average number of samples that fall in this interval is $Np(x)dx$, so for any function g(x)\n",
    "\n",
    "$$\\sum_{i=1}^Ng(x_i) \\approx \\int_a^b Np(x)g(x)dx$$\n",
    "\n",
    "Using this result, we can now write the general weighted average of the function g(x) as\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\langle g \\rangle_w = &\\frac{\\int_a^b w(x)g(x)dx}{\\int_a^b w(x)dx} \\\\ \n",
    "                    = &\\int_a^b p(x)g(x)dx \\\\\n",
    "              \\approx & \\frac{1}{N} \\sum_{i=1}^N g(x_i)    \n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "We then get \n",
    "\n",
    "$$ I \\approx  \\frac{1}{N} \\sum_{i=1}^{N} \\frac{f(x_i)}{w(x_i)} \\int_a^b w(x)dx$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 14.4 Homework\n",
    "\n",
    "Calculate a value for the integral\n",
    "\n",
    "$$ I = \\int_0 ^1 \\frac{x^{-1/2}}{e^x + 1} dx $$\n",
    "\n",
    "using the importance sampling method with $w(x) = x^{-1/2}$, with 1,000,000 random points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}