phys512/montecarlo_solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6d46a72f",
   "metadata": {},
   "source": [
    "# Monte Carlo Integration"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f71e7eb4",
   "metadata": {},
   "source": [
    "Integrate $\\exp(-|x^3|)$ from $0$ to $\\infty$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b53886bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.integrate\n",
    "\n",
    "seed = 239\n",
    "rng = np.random.default_rng(seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "5582478b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrate_MC(f, N, sampling = 'uniform'):\n",
    "\n",
    "    if sampling == 'gaussian':\n",
    "        # We need to normalize the Gaussian from 0 to infinity\n",
    "        p = lambda x: np.exp(-x**2/2) * np.sqrt(2/np.pi)\n",
    "        x = rng.normal(size = N)    \n",
    "    else:\n",
    "        # Generate x values between 0 and xmax\n",
    "        xmax = 10\n",
    "        p = lambda x: np.ones_like(x) / xmax\n",
    "        x = xmax*rng.uniform(size = N)\n",
    "    \n",
    "    # use np.mean to calculate the integral as an alternative to np.sum()/N\n",
    "    # also return the error in the mean, which is an error estimate for the integral\n",
    "    return np.mean(f(x)/p(x)), np.std(f(x)/p(x))/np.sqrt(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "5c39e032",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integral from scipy.quad = 0.89298 with error 2.74557e-09 (135 function evaluations)\n",
      "Uniform: I = 0.892806 +- 0.02417; frac err = 0.0270719; error estimate = 0.0248929\n",
      "Gaussian: I = 0.893114 +- 0.00462334; frac err = 0.00517665; error estimate = 0.00463627\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0lElEQVR4nO3de3wU9b3/8ffesrsBkkCUhECCEVHwCqJgkJaKUURQKdGqRYpooUcjClHRtEWLFmJtj6I2eOFBA1YoFhWqYqGSFjhWLgqlR20PoPArKZh4IwmSZLPZnd8fHOcY5LbJXvJNXs/HYx/97uzs7OfbuMN7Z77zHYdlWZYAAADixJnoAgAAQMdC+AAAAHFF+AAAAHFF+AAAAHFF+AAAAHFF+AAAAHFF+AAAAHFF+AAAAHHlTnQBhwuHw9q3b5+6dOkih8OR6HIAAMAJsCxLBw4cUFZWlpzOYx/baHPhY9++fcrOzk50GQAAoAUqKirUq1evY67T5sJHly5dJB0qPiUlJcHVAACAE1FbW6vs7Gz73/FjaXPh46tTLSkpKYQPALZgKKgX/vmCJOmm/jfJ4/IkuCIAR3IiQybaXPgAgCMJhoN6bMtjkqTrz7ie8AEYjPABwAhup1tX97nabgMwF99gAEZIciVp9rDZiS4DQBQwzwcAAIgrwgcAAIgrwgcAI9QF6zR0yVANXTJUdcG6RJcDoBUiCh+hUEgzZ85Ubm6u/H6/+vTpo4cffliWZdnrWJalBx54QD169JDf71d+fr527twZ9cIBdDwHggd0IHgg0WUAaKWIBpz+4he/0NNPP61FixbprLPO0rvvvqtJkyYpNTVVd955pyTp0Ucf1ZNPPqlFixYpNzdXM2fO1MiRI/WPf/xDPp8vJp0A0P753D69/t3X7TYAczmsrx+2OI4xY8YoIyNDCxYssJcVFBTI7/frhRdekGVZysrK0t1336177rlHklRTU6OMjAwtXLhQN9xww3E/o7a2VqmpqaqpqWGSMQAADBHJv98RnXYZOnSoysvLtWPHDknS3//+d7311lsaNWqUJGn37t2qrKxUfn6+/Z7U1FQNGTJEGzZsOOI2A4GAamtrmz0AAED7FdFpl/vvv1+1tbXq16+fXC6XQqGQZs+erfHjx0uSKisrJUkZGRnN3peRkWG/driSkhLNmjWrJbUD6ECC4aBe2vGSJOna06+Vx8kMp4CpIjry8fvf/16LFy/WkiVLtHXrVi1atEi/+tWvtGjRohYXUFxcrJqaGvtRUVHR4m0BaL+CoaDmbJqjOZvmKBgKJrocAK0Q0ZGPe++9V/fff789duOcc87Rv/71L5WUlGjixInKzMyUJFVVValHjx72+6qqqjRgwIAjbtPr9crr9bawfAAdhcvp0mW9L7PbAMwVUfioq6uT09n8YInL5VI4HJYk5ebmKjMzU+Xl5XbYqK2t1aZNm3TbbbdFp2IAHZLX5dVjA6ZLdZ9LVf88+orJ6VJadvwKAxCxiMLHVVddpdmzZysnJ0dnnXWW/va3v+mxxx7TLbfcIunQbXSnTZumn//85+rbt699qW1WVpbGjh0bi/oBdBTVFVLpYOl4E4x5kqXCzQQQoA2LKHw89dRTmjlzpm6//XZ98sknysrK0o9+9CM98MAD9jozZszQwYMHNWXKFFVXV2vYsGFatWoVc3wAaJ26zw8Fj3HzpZNOP/I6n+2QXpl8aF3CB9BmRTTPRzwwzweAI6mv2KQxf7pZ6nSyXr92tfxu/zdX2rdNem64NGWdlDUgzhUCHVsk/35HdOQDABLFkqVP3G4psF9t7DcTgAgRPgAYwetM0rK9H0sFC+R1cYUcYDLCBwAjuBxO9WsMSl16S1xqCxgtoknGAAAAWovwAcAIwXCTVnTupBX71isYZoZTwGScdgHQNlRXHLpE9iiCn/6PZp6cLn0wX5ef90Pu7QIYjPABIPFOYAIxl0P6VkamdMrFTK8OGI7wASDxTmACMa+keUydDrQLhA8AbcdJpzM5GNABMOAUAADEFeEDgBHqm+o1+pXRGv3KaNU31Se6HACtwGkXAEawLEt7Duyx2wDMRfgAYASvy6vnRz1vtwGYi/ABwAgup0sDuw9MdBkAooAxHwAAIK448gHACE3hJpXvKZckXZpzqdxOdl+Aqfj2AjBCY6hR96y7R5K06fubCB+Awfj2AjCC0+HUBRkX2G0A5iJ8ADCCz+1T2RVliS4DQBTw8wEAAMQV4QMAAMQV4QOAERqaGnTtq9fq2levVUNTQ6LLAdAKjPkAYISwFdb2/dvtNgBzET4AGMHr8urZy5612wDMRfgAYASX06WhWUMTXQaAKGDMBwAAiCuOfAAwQlO4SW/ve1uSNDRrKDOcAgbj2wvACI2hRhWWF0pienXAdHx7ARjB6XDqrPSz7DYAcxE+ABjB5/Zp6ZiliS4DQBTw8wEAAMQV4QMAAMQV4QOAERqaGjThjQma8MYEplcHDBdR+DjllFPkcDi+8SgsPDQCvaGhQYWFhUpPT1fnzp1VUFCgqqqqmBQOoGMJW2Ft+3Sbtn26jenVAcNFNOD0nXfeUSgUsp+///77uuyyy3TddddJkqZPn66VK1dq2bJlSk1N1R133KFx48bpr3/9a3SrBtDhJLmSNPeSuXYbgLkiCh8nn3xys+ePPPKI+vTpo+HDh6umpkYLFizQkiVLNGLECElSWVmZ+vfvr40bN+qiiy6KXtUAOhy3061Lcy5NdBkAoqDFYz4aGxv1wgsv6JZbbpHD4dCWLVsUDAaVn59vr9OvXz/l5ORow4YNUSkWAACYr8XzfKxYsULV1dW6+eabJUmVlZVKSkpSWlpas/UyMjJUWVl51O0EAgEFAgH7eW1tbUtLAtCOhcIhbf1kqyTp/O7ny+V0JbgiAC3V4iMfCxYs0KhRo5SVldWqAkpKSpSammo/srOzW7U9AO1TIBTQLatv0S2rb1EgFDj+GwC0WS0KH//617+0Zs0a/fCHP7SXZWZmqrGxUdXV1c3WraqqUmZm5lG3VVxcrJqaGvtRUVHRkpIAtHMOh0N9UvuoT2ofORyORJcDoBVadNqlrKxM3bt31+jRo+1lgwYNksfjUXl5uQoKCiRJ27dv1549e5SXl3fUbXm9Xnm93paUAaAD8bv9WjF2RaLLABAFEYePcDissrIyTZw4UW73/709NTVVt956q4qKitStWzelpKRo6tSpysvL40oXAABgizh8rFmzRnv27NEtt9zyjdcef/xxOZ1OFRQUKBAIaOTIkZo3b15UCgUAAO1DxOHj8ssvl2VZR3zN5/OptLRUpaWlrS4MAL6uoalBU/88VZL01Iin5HP7ElwRgJZq8aW2ABBPYSusjR9vtNsAzEX4AGCEJFeSSr5VYrcBmIvwAcAIbqdbY04dk+gyAERBiycZAwAAaAmOfAAwQigc0j+/+KckqX+3/kyvDhiM8AHACIFQQDeuvFGStOn7m5TsTE5wRQBaivABwAgOh0NZnbLsNgBzET4AGMHv9mv1tasTXQaAKGDAKQAAiCvCBwAAiCvCBwAjBEIB3fnnO3Xnn+9UIBRIdDkAWoExHwCMEAqH9JeKv9htcaUtYCzCBwAjeFwePZj3oN0GYC7CBwAjeJweXXv6tYkuA0AUMOYDAADEFUc+ABghbIW1q3qXJOnUtFPldPDbCTAV4QOAERqaGvTdV78r6X+nV/cwvTpgKsIHAGN09XZNdAkAooDwAcAIyZ5krb9hfaLLABAFnDQFAABxRfgAAABxRfgAYIRAKKD71t+n+9bfx/TqgOEIHwCMEAqH9MbuN/TG7jcOTa8OwFgMOAVgBI/LoxkXzrDbAMxF+ABgBI/TowlnTkh0GQCigNMuAAAgrjjyAcAIYSusjw9+LEnq0akH06sDBiN8ADBCQ1ODrnj5CklMrw6YjvABwBh+t//EVvxsx7FfT06X0rJbXxCAFiF8ADBCsidZm8dvPs5K6ZInWXpl8rHX8yRLhZsJIECCED4AtB9p2YdCRd3nR1/nsx2Hwknd54QPIEEIHwDal7RsQgXQxhE+ABihMdSoOZvmSJJ+POTHSnIlJbgiAC3FtWoAjNAUbtLLO1/WyztfVlO4KdHlAGiFiMPH3r17ddNNNyk9PV1+v1/nnHOO3n33Xft1y7L0wAMPqEePHvL7/crPz9fOnTujWjSAjsfj9GjqwKmaOnCqPE6mVwdMFtFpl/379+viiy/WJZdcoj/+8Y86+eSTtXPnTnXt2tVe59FHH9WTTz6pRYsWKTc3VzNnztTIkSP1j3/8Qz6fL+odANAxeFweTTl3SqLLABAFEYWPX/ziF8rOzlZZWZm9LDc3125blqW5c+fqpz/9qa655hpJ0vPPP6+MjAytWLFCN9xwQ5TKBgAAporotMurr76qCy64QNddd526d++ugQMHav78+fbru3fvVmVlpfLz8+1lqampGjJkiDZs2HDEbQYCAdXW1jZ7AMDhLMvSFw1f6IuGL2RZVqLLAdAKEYWPXbt26emnn1bfvn21evVq3Xbbbbrzzju1aNEiSVJlZaUkKSMjo9n7MjIy7NcOV1JSotTUVPuRnc0lcgC+qb6pXsNfHK7hLw5XfVN9ossB0AoRhY9wOKzzzz9fc+bM0cCBAzVlyhRNnjxZzzzzTIsLKC4uVk1Njf2oqKho8bYAAEDbF9GYjx49eujMM89stqx///56+eWXJUmZmZmSpKqqKvXo0cNep6qqSgMGDDjiNr1er7xebyRlAOiAkj3Jem/ie4kuA0AURHTk4+KLL9b27dubLduxY4d69+4t6dDg08zMTJWXl9uv19bWatOmTcrLy4tCuQAAwHQRHfmYPn26hg4dqjlz5uh73/ueNm/erOeee07PPfecJMnhcGjatGn6+c9/rr59+9qX2mZlZWns2LGxqB8AABgmovBx4YUXavny5SouLtZDDz2k3NxczZ07V+PHj7fXmTFjhg4ePKgpU6aourpaw4YN06pVq5jjA0CrNIYa9fiWxyVJ0wdNZ3p1wGAR39tlzJgxGjNmzFFfdzgceuihh/TQQw+1qjAA+LqmcJNe+OcLkqSpA6cSPgCDcWM5AEbwOD2afM5kuw3AXIQPAEbwuDy68/w7E10GgCjgrrYAACCuOPIBwAiWZdkzm/rdfjkcjgRXBKClOPIBwAj1TfUasmSIhiwZwvTqgOEIHwAAIK447QLACH63X5u+v8luAzAX4QOAERwOh5I9yYkuA0AUcNoFAADEFUc+ABghGArq6b8/LUm67bzb5HEx0RhgKsIHACMEw0HNf2++JOmH5/yQ8AEYjPABwAhup1s39b/JbgMwF99gAEZIciXpvsH3JboMAFHAgFMAABBXhA8AABBXhA8ARqgL1umcRefonEXnqC5Yl+hyALQC4QMAAMQVA04BGMHv9mvd9evsNgBzET4AGMHhcKibr1uiywAQBZx2AQAAccWRDwBGCIaCKvugTJI06axJzHAKGIzwAcAIwXBQT/3tKUnSTf1vInwABiN8ADCC2+lWQd8Cuw3AXHyDARghyZWknw39WaLLABAFDDgFAABxRfgAAABxRfgAYIS6YJ0GLx6swYsHM706YDjGfAAwRn1TfaJLABAFhA8ARvC5fVpVsMpuAzAX4QOAEZwOp3p27pnoMgBEAWM+AABAXHHkA4ARguGglv7PUknSDf1ukMfJDKeAqQgfAIwQDAX16DuPSpIK+hYQPgCDRXTa5Wc/+5kcDkezR79+/ezXGxoaVFhYqPT0dHXu3FkFBQWqqqqKetEAOh6X06Urc6/UlblXyuV0JbocAK0Q8ZGPs846S2vWrPm/Dbj/bxPTp0/XypUrtWzZMqWmpuqOO+7QuHHj9Ne//jU61QLosLwur37x7V8kugwAURBx+HC73crMzPzG8pqaGi1YsEBLlizRiBEjJEllZWXq37+/Nm7cqIsuuqj11QIAAONFfLXLzp07lZWVpVNPPVXjx4/Xnj17JElbtmxRMBhUfn6+vW6/fv2Uk5OjDRs2RK9iAABgtIiOfAwZMkQLFy7UGWecoY8//lizZs3St771Lb3//vuqrKxUUlKS0tLSmr0nIyNDlZWVR91mIBBQIBCwn9fW1kbWAwAdQl2wTle8fIUkaVXBKiV7khNcEYCWiih8jBo1ym6fe+65GjJkiHr37q3f//738vv9LSqgpKREs2bNatF7AXQs+wP7E10CgCho1SRjaWlpOv300/Xhhx8qMzNTjY2Nqq6ubrZOVVXVEceIfKW4uFg1NTX2o6KiojUlAWinfG6fll+9XMuvXs706oDhWhU+vvzyS3300Ufq0aOHBg0aJI/Ho/Lycvv17du3a8+ePcrLyzvqNrxer1JSUpo9AOBwTodTp3U9Tad1PU1OB5MzAyaL6LTLPffco6uuukq9e/fWvn379OCDD8rlcunGG29Uamqqbr31VhUVFalbt25KSUnR1KlTlZeXx5UuAADAFlH4+Pe//60bb7xRn3/+uU4++WQNGzZMGzdu1MknnyxJevzxx+V0OlVQUKBAIKCRI0dq3rx5MSkcQMcSDAf1hw//IEm65rRrmOEUMFhE4WPp0qXHfN3n86m0tFSlpaWtKgoADhcMBTVrw6HB6VfmXkn4AAzGvV0AGMHldOmS7EvsNgBzET4AGMHr8urJEU8mugwAUcCQcQAAEFcc+QAQW9UVUt3nx17nsx3xqQVAm0D4ABA71RVS6WApWHf8dT3JUnL6UV+ub6rX2BVjJUkrxq6Q392yWZUBJB7hA0Ds1H1+KHiMmy+ddPqx101Ol9Kyj/qyZVnad3Cf3QZgLsIHgNg76XQpa0CrNuF1efW70b+z2wDMRfgAYASX06WzTzo70WUAiAKudgEAAHHFkQ8ARmgKN2nV/1slSbrilCvkdrL7AkzFtxeAERpDjSr+r2JJ0ojsEYQPwGB8ewEYwelw6qIeF9ltAOYifAAwgs/t0/zL5ye6DABRwM8HAAAQV4QPAAAQV4QPAEb4anr1sSvGqr6pPtHlAGgFxnwAMIJlWfqo5iO7DcBchA8ARvC6vPrNyN/YbQDmInwAMILL6dKFmRcmugwAUcCYDwAAEFcc+QBghKZwk9b9e50kaXiv4cxwChiMby8AIzSGGjXtL9MkSZu+v4nwARiMby8AIzgdTg04eYDdBmAuwgcAI/jcPv32yt8mugwAUcDPBwAAEFeEDwAAEFeEDwBGaGhq0A2v36AbXr9BDU0NiS4HQCsw5gOAEcJWWB98/oHdBmAuwgcAIyS5klR6aandBmAuwgcAI7idbn2717cTXQaAKGDMBwAAiCuOfAAwQigc0qbKTZKkIZlD5HK6ElwRgJYifAAwQiAU0I/e/JGkQ9OrJzuTE1wRgJYifAAwgtPh1Bldz7DbAMzVqm/wI488IofDoWnTptnLGhoaVFhYqPT0dHXu3FkFBQWqqqpqbZ0AOjif26eXrn5JL139knxuX6LLAdAKLQ4f77zzjp599lmde+65zZZPnz5dr732mpYtW6Z169Zp3759GjduXKsLBQAA7UOLwseXX36p8ePHa/78+eratau9vKamRgsWLNBjjz2mESNGaNCgQSorK9Pbb7+tjRs3Rq1oAABgrhaFj8LCQo0ePVr5+fnNlm/ZskXBYLDZ8n79+iknJ0cbNmw44rYCgYBqa2ubPQDgcA1NDZq0apImrZrE9OqA4SIecLp06VJt3bpV77zzzjdeq6ysVFJSktLS0potz8jIUGVl5RG3V1JSolmzZkVaBoAOJmyF9W7Vu3YbgLkiCh8VFRW666679Oabb8rni86Ar+LiYhUVFdnPa2trlZ2dHZVtA2g/klxJ+tXwX9ltAOaKKHxs2bJFn3zyic4//3x7WSgU0vr16/XrX/9aq1evVmNjo6qrq5sd/aiqqlJmZuYRt+n1euX1eltWPYAOw+10a+QpIxNdBoAoiCh8XHrppXrvvfeaLZs0aZL69eun++67T9nZ2fJ4PCovL1dBQYEkafv27dqzZ4/y8vKiVzUAADBWROGjS5cuOvvss5st69Spk9LT0+3lt956q4qKitStWzelpKRo6tSpysvL00UXXRS9qgF0OKFwSP/92X9Lks496VymVwcMFvUZTh9//HE5nU4VFBQoEAho5MiRmjdvXrQ/BkAHEwgF9IM//kAS06sDpmt1+Fi7dm2z5z6fT6WlpSotLW3tpgHA5nA4lNMlx24DMBf3dgFgBL/br5XjVia6DABRwN2ZAABAXBE+AABAXBE+ABghEAro9jW36/Y1tysQCiS6HACtwJgPAEYIhUP6r73/ZbfFlbaAsQgfAIzgcXn08MUP220A5iJ8ADCCx+nR2NPGJroMAFHAmA8AABBXHPkAYIRQOKSd1TslSX3T+jK9OmAwwgcAIwRCAV332nWSmF4dMB3hA4ARHA6Huvu7220A5iJ8ADCC3+1X+ffKE10GgChgwCkAAIgrwgcAAIgrwgcAIwRCARWtLVLR2iKmVwcMR/gAYIRQOKQ3//Wm3vzXm4emVwdgLAacAjCCx+XRj4f82G4DMBfhA4ARPE6Pbux3Y6LLABAFnHYBAABxxZEPAEYIW2FVHKiQJGV3yZbTwW8nwFSEDwBGaGhq0JjlYyT97/TqHqZXB0xF+ABgjC6eLokuAUAUED4AGCHZk6y3v/92ossAEAWcNAUAAHFF+AAAAHFF+ABghMZQo37y1k/0k7d+osZQY6LLAdAKhA8ARmgKN+nVj17Vqx+9qqZwU6LLAdAKDDgFYASP06OiQUV2G4C5CB8AjOBxeTTp7EmJLgNAFHDaBQAAxBVHPgAYIWyF9Wndp5Kkk5NPZnp1wGCEDwBGaGhqUP5L+ZKYXh0wHeEDgDHcDnZZQHvANxmAEZI9yfrbD/6W6DIAREFEJ02ffvppnXvuuUpJSVFKSory8vL0xz/+0X69oaFBhYWFSk9PV+fOnVVQUKCqqqqoFw0AAMwVUfjo1auXHnnkEW3ZskXvvvuuRowYoWuuuUYffPCBJGn69Ol67bXXtGzZMq1bt0779u3TuHHjYlI4AAAwU0SnXa666qpmz2fPnq2nn35aGzduVK9evbRgwQItWbJEI0aMkCSVlZWpf//+2rhxoy666KLoVQ2gw2kMNerRdx6VJM24cIaSXEkJrghAS7X4WrVQKKSlS5fq4MGDysvL05YtWxQMBpWfn2+v069fP+Xk5GjDhg1H3U4gEFBtbW2zBwAcrincpBe3v6gXt7/I9OqA4SIecPree+8pLy9PDQ0N6ty5s5YvX64zzzxT27ZtU1JSktLS0pqtn5GRocrKyqNur6SkRLNmzYq4cAAdi8fp0W3n3Wa3AZgr4vBxxhlnaNu2baqpqdFLL72kiRMnat26dS0uoLi4WEVFRfbz2tpaZWdnt3h7ANonj8uj2wfcnugyAERBxOEjKSlJp512miRp0KBBeuedd/TEE0/o+uuvV2Njo6qrq5sd/aiqqlJmZuZRt+f1euX1eiOvHAAAGKnV8xOHw2EFAgENGjRIHo9H5eXl9mvbt2/Xnj17lJeX19qPAdDBWZal2sZa1TbWyrKsRJcDoBUiOvJRXFysUaNGKScnRwcOHNCSJUu0du1arV69Wqmpqbr11ltVVFSkbt26KSUlRVOnTlVeXh5XugBotfqmel38u4slMb06YLqIwscnn3yiH/zgB/r444+Vmpqqc889V6tXr9Zll10mSXr88cfldDpVUFCgQCCgkSNHat68eTEpHAAAmCmi8LFgwYJjvu7z+VRaWqrS0tJWFQUAh/O7/do6Yask7vECmI5vMAAjOBwOeRxcYgu0B60ecAoAABAJjnwAMEIwFNSTf3tSknTnwDvlcXEUBDAVRz4AGCEYDmrhBwu18IOFCoaDiS4HQCtw5AOAETxOj24+62a7DcBchA8ARvC4PLr7grsTXQaAKCB8AOiYPttx7NeT06U07jMFxALhA4ARLMtSk9Uk6dA8Hw6Ho2UbSk6XPMnSK5OPvZ4nWSrcTAABYoDwAcAI9U31GrJkiKRWTq+eln0oVNR9fvR1PttxKJzUfU74AGKA8AGg40nLJlQACUT4AGAEv9uvv974V7sNwFyEDwBGcDgcSklKSXQZAKKAScYAAEBcceQDgBGCoaDmvzdfkjT5nMlMrw4YjPABwAjBcFBP//1pSdLNZ91M+AAMRvgAYAS3063rz7jebgMwF99gAEZIciXppxf9NNFlAIgCBpwCAIC4InwAAIC4InwAMEJdsE4Dnx+ogc8PVF2wLtHlAGgFxnwAMMZXN5YDYDbCBwAj+Nw+rbl2jd0GYC7CBwAjOB1OZXTKSHQZAKKAMR8AACCuOPIBwAjBUFAv/PMFSdJN/W9ihlPAYIQPAEYIhoN6bMtjkqTrz7ie8AEYjPABwAhup1tX97nabgMwF99gAEZIciVp9rDZiS4DQBQw4BQAAMQV4QMAAMQV4QOAEeqCdRq6ZKiGLhnK9OqA4RjzAcAYB4IHEl0CgCggfAAwgs/t0+vffd1uAzBXRKddSkpKdOGFF6pLly7q3r27xo4dq+3btzdbp6GhQYWFhUpPT1fnzp1VUFCgqqqqqBYNoONxOpzqndJbvVN6y+ngjDFgsoi+wevWrVNhYaE2btyoN998U8FgUJdffrkOHjxorzN9+nS99tprWrZsmdatW6d9+/Zp3LhxUS8cAACYKaLTLqtWrWr2fOHCherevbu2bNmib3/726qpqdGCBQu0ZMkSjRgxQpJUVlam/v37a+PGjbrooouiVzmADiUYDuqlHS9Jkq49/Vp5nMxwCpiqVWM+ampqJEndunWTJG3ZskXBYFD5+fn2Ov369VNOTo42bNhwxPARCAQUCATs57W1ta0pCUA7FQwFNWfTHEnSNX2uIXwABmtx+AiHw5o2bZouvvhinX322ZKkyspKJSUlKS0trdm6GRkZqqysPOJ2SkpKNGvWrJaWAaCDcDlduqz3ZXYbgLlaHD4KCwv1/vvv66233mpVAcXFxSoqKrKf19bWKjs7u1XbBND+eF1ePfadxxJdBoAoaFH4uOOOO/T6669r/fr16tWrl708MzNTjY2Nqq6ubnb0o6qqSpmZmUfcltfrldfrbUkZABKtukKq+/zor3+2I361ADBGROHDsixNnTpVy5cv19q1a5Wbm9vs9UGDBsnj8ai8vFwFBQWSpO3bt2vPnj3Ky8uLXtUAEq+6QiodLB1vtlFPspScHp+aABghovBRWFioJUuW6A9/+IO6dOlij+NITU2V3+9Xamqqbr31VhUVFalbt25KSUnR1KlTlZeXx5UuQHtT9/mh4DFuvnTS6UdfLzldSmv9qdT6pnqNeWWMJOn1ca/L7/a3epsAEiOi8PH0009Lkr7zne80W15WVqabb75ZkvT444/L6XSqoKBAgUBAI0eO1Lx586JSLIA26KTTpawBMf8Yy7L0Sf0ndhuAuSI+7XI8Pp9PpaWlKi0tbXFRAHA4r8urZVcts9sAzMW9XQAYweV0qV+3fokuA0AUcIMEAAAQVxz5AGCEYDiolbtWSpJGnzqaGU4BgxE+ABghGApq5l9nSpIu73054QMwGOEDgBFcTpe+1fNbdhuAuQgfAIzgdXk1L5/L9oH2gAGnAAAgrggfAAAgrggfAIxQ31Sv0a+M1uhXRqu+qT7R5QBoBcZ8ADCCZVnac2CP3QZgLsIHACN4XV49P+p5uw3AXIQPAHG1t7pe+w82tui9Hp0qSfrnx19Gs6Rv8H32pU6T9OGnX6rBqlHXTknqmcZddIFoIXwAiJu91fXK/891qg+GEl3KMZ3l2K2VXumupdv0gVUjv8elNXcPJ4AAUUL4ABA3+w82qj4Y0tzrB+i07p0jem8o3KR3Pl0vSbrw5G/L5Yzd7sv3Waq0XHrihgF6P5yraS9u0/6DjYQPIEoIHwDi7rTunXV2z9SI3lMXrNNNf35AkrTp+5uU7EmORWmHOA4Fo9NO7qwGK7KQBOD4CB8AjOB0OHVBxgV2G4C5CB8AjOBz+1R2RVmiywAQBfx8AAAAcUX4AAAAccVpF8BgrZkzo7UOnwvjRHz4Scvn52hoatBNb9wkSXrhyhfkc/tavC0AiUX4AAyV6DkzDp8L40T5PS517ZQU8eeFrbC2799utwGYi/ABGKo1c2ZEw9fnwmg46ZwTfl9LZwv1urx69rJn7TYAcxE+AMO1ZM6MqPjaXBjKiv3nu5wuDc0aGvPPARB7hA8AOAGtGa+SKNyTBm0V4QOAEZrCTXp739uSpKFZQ+WO4fTqX9e1U5L8HpemvbgtLp8XTdyTBm0V4QOAERpDjSosL5R0aHr1eIWPnml+rbl7eMKuKmqpDz/5knvSoM0ifAAwgtPh1FnpZ9nteOqZ5ucfcCCKCB8AjOBz+7R0zNJElwEgCpjhFAAAxBXhAwAAxBXhA4ARGpoaNOGNCZrwxgQ1NDUkuhwArcCYDwBGCFthbft0m90GYC7CBwAjJLmSNPeSuXYbgLkIHwCM4Ha6dWnOpYkuA0AURDzmY/369brqqquUlZUlh8OhFStWNHvdsiw98MAD6tGjh/x+v/Lz87Vz585o1QsAAAwXcfg4ePCgzjvvPJWWlh7x9UcffVRPPvmknnnmGW3atEmdOnXSyJEj1dDAADEALRcKh/RO5Tt6p/IdhcKhRJcDoBUiPu0yatQojRo16oivWZaluXPn6qc//amuueYaSdLzzz+vjIwMrVixQjfccEPrqgXQYQVCAd2y+hZJh6ZXT3YmJ7giAC0V1Uttd+/ercrKSuXn59vLUlNTNWTIEG3YsOGI7wkEAqqtrW32AIDDORwO9Untoz6pfeRwOBJdDoBWiOqA08rKSklSRkZGs+UZGRn2a4crKSnRrFmzolkGgHbI7/ZrxdgViS4DQBQkfJKx4uJi1dTU2I+KiopElwQAAGIoquEjMzNTklRVVdVseVVVlf3a4bxer1JSUpo9AABA+xXV8JGbm6vMzEyVl5fby2pra7Vp0ybl5eVF86MAdDANTQ2a/KfJmvynyUyvDhgu4jEfX375pT788EP7+e7du7Vt2zZ169ZNOTk5mjZtmn7+85+rb9++ys3N1cyZM5WVlaWxY8dGs24AHUzYCmvjxxvtNgBzRRw+3n33XV1yySX286KiIknSxIkTtXDhQs2YMUMHDx7UlClTVF1drWHDhmnVqlXy+XzRqxpAh5PkSlLJt0rsNgBzRRw+vvOd78iyrKO+7nA49NBDD+mhhx5qVWEA8HVup1tjTh2T6DIAREHCr3YBAAAdCzeWA/7X3up67T/YmOgyTtiHn3yZ6BLiKhQO6Z9f/FOS1L9bf7mcrgRXBKClCB+ADgWP/P9cp/qgWfcM8Xtc6tqpY4x/CIQCunHljZKYXh0wHeEDkLT/YKPqgyHNvX6ATuveOdHlnLCunZLUM82f6DLiwuFwKKtTlt0GYC7CB/A1p3XvrLN7pia6DByB3+3X6mtXx/dDP9tx7NeT06W07PjUArQjhA8AOFxyuuRJll6ZfOz1PMlS4eY2HUBMGxvUkY7mdWSEDwA4XFr2oVBR9/nR1/lsx6FwUvd5mwwfXTslye9xadqL2xJdSkT8HpfW3D2cANLOET4AGCEQCujedfdKkn45/Jfyuryx/cC07DYZKk5UzzS/1tw93LgruKa9uE37DzYSPto5wgcAI4TCIf2l4i92W1xpe1w90/z8I442ifABwAgel0cP5j1otwGYi/ABwAgep0fXnn5tossAEAVMrw4AAOKKIx8AjBC2wtpVvUuSdGraqXI6+O0EmIrwAcAIDU0N+u6r35X0v9Ore5heHTAV4QOAMbp6uya6BABRQPgAYIRkT7LW37A+0WUAiAJOmgIAgLgifAAAgLgifAAwQiAU0H3r79N96+9TIBRIdDkAWoExH4iJvdX1xt1TAoeprjj+jdXiKBQO6Y3db0jSoZlOmV4dMBbhA1G3t7pe+f+5TvXBUKJLiYjf41LXTkmJLqNtqK6QSgdLwbpjr+dJPnT7+TjwuDyaceEMuw3AXIQPRN3+g42qD4Y09/oBOq1750SXc8K6dkriJlxfqfv8UPAYN1866fSjr5ecHrc7v3qcHk04c0JcPgtAbBE+EDOnde+ss3umJroMtMZJp0tZAxJdBYB2hvABwAhhK6yPD34sSerRqQfTqwMGI3wAMEJDU4OuePkKSUyv3t6ZNgCcU7aRI3wAMIbfzQ6+PevaKUl+j0vTXtyW6FIi4ve4tObu4QSQCBA+ABgh2ZOszeM3J7qMbzreJcdxHJRrup5pfq25e7hxl+lPe3Gb9h9sJHxEgPABAC2RnH7oUuNXJh97PU+yVLiZAHKCeqb5+Ue8AyB8AEBLpGUfChXHm4jtlcmH1iF8ADbCBwAjNIYaNWfTHEnSj4f8WEmuNjAhXFo2oQJoAcIH0BG1sanTT0RTuEkv73xZkjTjwhltI3wAaBHCB9DRtMGp00+Ex+nR1IFT7TYAcxE+DMBN2hBVbXDq9BPhcXk05dwpiS4DQBTELHyUlpbql7/8pSorK3Xeeefpqaee0uDBg2P1ce0WN2lDzDB1OoAEiUn4ePHFF1VUVKRnnnlGQ4YM0dy5czVy5Eht375d3bt3j8VHtlvcpA04xLIs7Q/slyR19XaVw+FIcEUAWiom4eOxxx7T5MmTNWnSJEnSM888o5UrV+o3v/mN7r///lh85Akz9RQGN2lDR1ffVK/hLw6XxPTqaHtMO92c6B+IUQ8fjY2N2rJli4qLi+1lTqdT+fn52rBhwzfWDwQCCgQC9vOamhpJUm1tbbRL077qel3967fUEAxHfdux5PM45Q41qLaWX3qIggNfSgHr0P/G4HsWK3XBOoXqD51+rK2tVZOnKcEVnQBD/7/GiXOHGpQUbtCdz7+d6FIi4vM49eodw5QVxQDy1b/blmUdf2Uryvbu3WtJst5+++1my++9915r8ODB31j/wQcftCTx4MGDBw8ePNrBo6Ki4rhZIeFXuxQXF6uoqMh+Hg6H9cUXXyg9Pd2oc7q1tbXKzs5WRUWFUlJSEl1OzHWk/nakvkr0tz3rSH2V6G+8WZalAwcOKCsr67jrRj18nHTSSXK5XKqqqmq2vKqqSpmZmd9Y3+v1yuv1NluWlpYW7bLiJiUlpUP8R/6VjtTfjtRXif62Zx2prxL9jafU1NQTWs8Z7Q9OSkrSoEGDVF5ebi8Lh8MqLy9XXl5etD8OAAAYJianXYqKijRx4kRdcMEFGjx4sObOnauDBw/aV78AAICOKybh4/rrr9enn36qBx54QJWVlRowYIBWrVqljIyMWHxcm+D1evXggw9+4xRSe9WR+tuR+irR3/asI/VVor9tmcOyTuSaGAAAgOiI+pgPAACAYyF8AACAuCJ8AACAuCJ8AACAuCJ8HENpaalOOeUU+Xw+DRkyRJs3bz7m+nPnztUZZ5whv9+v7OxsTZ8+XQ0NDfbrJSUluvDCC9WlSxd1795dY8eO1fbt22PdjRMS7b5+3SOPPCKHw6Fp06bFoPKWiUV/9+7dq5tuuknp6eny+/0655xz9O6778ayGycs2v0NhUKaOXOmcnNz5ff71adPHz388MMndk+HGIukr8FgUA899JD69Okjn8+n8847T6tWrWrVNuMt2v1ty/spKTZ/36+0tX1VLPraZvZTUbmhSzu0dOlSKykpyfrNb35jffDBB9bkyZOttLQ0q6qq6ojrL1682PJ6vdbixYut3bt3W6tXr7Z69OhhTZ8+3V5n5MiRVllZmfX+++9b27Zts6688korJyfH+vLLL+PVrSOKRV+/snnzZuuUU06xzj33XOuuu+6KcU9OTCz6+8UXX1i9e/e2br75ZmvTpk3Wrl27rNWrV1sffvhhvLp1VLHo7+zZs6309HTr9ddft3bv3m0tW7bM6ty5s/XEE0/Eq1tHFGlfZ8yYYWVlZVkrV660PvroI2vevHmWz+eztm7d2uJtxlMs+ttW91OWFZv+fqWt7ati0de2tJ8ifBzF4MGDrcLCQvt5KBSysrKyrJKSkiOuX1hYaI0YMaLZsqKiIuviiy8+6md88sknliRr3bp10Sm6hWLV1wMHDlh9+/a13nzzTWv48OFt4gttWbHp73333WcNGzYsNgW3Uiz6O3r0aOuWW25pts64ceOs8ePHR7HyyEXa1x49eli//vWvmy07vB+RbjOeYtHfw7WV/ZRlxa6/bXFfFYu+tqX9FKddjqCxsVFbtmxRfn6+vczpdCo/P18bNmw44nuGDh2qLVu22IfFdu3apTfeeENXXnnlUT+npqZGktStW7coVh+ZWPa1sLBQo0ePbrbtRItVf1999VVdcMEFuu6669S9e3cNHDhQ8+fPj21nTkCs+jt06FCVl5drx44dkqS///3veuuttzRq1KgY9ubYWtLXQCAgn8/XbJnf79dbb73V4m3GSyz6eyRtYT8lxba/bW1fFau+tqn9VKLTT1u0d+9eS5L19ttvN1t+7733WoMHDz7q+5544gnL4/FYbrfbkmT9x3/8x1HXDYVC1ujRo495ZCQeYtXX3/3ud9bZZ59t1dfXW5ZltZlfE7Hqr9frtbxer1VcXGxt3brVevbZZy2fz2ctXLgwJv04UbHqbygUsu677z7L4XBYbrfbcjgc1pw5c2LShxPVkr7eeOON1plnnmnt2LHDCoVC1p/+9CfL7/dbSUlJLd5mvMSiv4drK/spy4pdf9vivipWfW1L+ymOfETJ2rVrNWfOHM2bN09bt27VK6+8opUrV+rhhx8+4vqFhYV6//33tXTp0jhX2nrH62tFRYXuuusuLV68+BtJ3EQn8rcNh8M6//zzNWfOHA0cOFBTpkzR5MmT9cwzzySw8pY5kf7+/ve/1+LFi7VkyRJt3bpVixYt0q9+9SstWrQogZVH7oknnlDfvn3Vr18/JSUl6Y477tCkSZPkdLbPXWOk/TV5PyUdv7/taV91In/bNrWfinvcMUAgELBcLpe1fPnyZst/8IMfWFdfffUR3zNs2DDrnnvuabbst7/9reX3+61QKNRseWFhodWrVy9r165dUa27JWLR1+XLl1uSLJfLZT8kWQ6Hw3K5XFZTU1OsunNcsfrb5uTkWLfeemuzdebNm2dlZWVFr/gWiFV/e/Xq9Y3zyw8//LB1xhlnRK/4CLWkr1+pr6+3/v3vf1vhcNiaMWOGdeaZZ7Z6m7EWi/5+XVvaT1lWbPrbVvdVsfrbtqX9VPuM962UlJSkQYMGqby83F4WDodVXl6uvLy8I76nrq7uG78eXC6XJNmXH1qWpTvuuEPLly/Xn//8Z+Xm5saoBycuFn299NJL9d5772nbtm3244ILLtD48eO1bds2e91EiNXf9uKLL/7G5Yg7duxQ7969o1l+xGLV36OtEw6Ho1l+RFrS16/4fD717NlTTU1Nevnll3XNNde0epuxFov+Sm1zPyXFpr9tdV8Vq79tm9pPxT3uGGLp0qWW1+u1Fi5caP3jH/+wpkyZYqWlpVmVlZWWZVnWhAkTrPvvv99e/8EHH7S6dOli/e53v7N27dpl/elPf7L69Oljfe9737PXue2226zU1FRr7dq11scff2w/6urq4t6/r4tFXw/XFs6jfiUW/d28ebPldrut2bNnWzt37rQWL15sJScnWy+88ELc+3e4WPR34sSJVs+ePe1LbV955RXrpJNOsmbMmBH3/n1dpH3duHGj9fLLL1sfffSRtX79emvEiBFWbm6utX///hPeZiLFor9tdT9lWbHp7+Hayr4qFn1tS/spwscxPPXUU1ZOTo6VlJRkDR482Nq4caP92vDhw62JEyfaz4PBoPWzn/3M6tOnj+Xz+azs7Gzr9ttvb/aHl3TER1lZWfw6dRTR7uvh2soX+iux6O9rr71mnX322ZbX67X69etnPffcc3HqzfFFu7+1tbXWXXfdZeXk5Fg+n8869dRTrZ/85CdWIBCIY6+OLJK+rl271urfv7/l9Xqt9PR0a8KECdbevXsj2maiRbu/bXk/ZVmx+ft+XVvaV8Wir21lP+WwrDYwJSEAAOgwGPMBAADiivABAADiivABAADiivABAADiivABAADiivABAADiivABAADiivABAADiivABAADiivABAADiivABAADiivABAADi6v8DeyrgvNLFPp0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def f(x):\n",
    "    return np.exp(-np.abs(x**3))\n",
    "\n",
    "# Compute the integral using quad form comparison\n",
    "I0, err, info = scipy.integrate.quad(f, 0.0, np.inf, full_output = True)\n",
    "print(\"Integral from scipy.quad = %g with error %g (%d function evaluations)\" % \n",
    "      (I0, err, info['neval']))\n",
    "\n",
    "# Number of samples for the MC integration\n",
    "N = 10**4\n",
    "\n",
    "# Now do the integration 1000 times for both uniform and Gaussian\n",
    "n_trials = 1000\n",
    "I_arr = np.zeros(n_trials)\n",
    "I2_arr = np.zeros(n_trials)\n",
    "\n",
    "for i in range(n_trials):\n",
    "    I_arr[i], err = integrate_MC(f, N)\n",
    "    I2_arr[i], err2 = integrate_MC(f, N, sampling = 'gaussian')\n",
    "    \n",
    "I_mean = np.mean(I_arr)\n",
    "I_std = np.std(I_arr)\n",
    "print(\"Uniform: I = %g +- %g; frac err = %g; error estimate = %g\" % (I_mean, I_std, I_std/I_mean, err))\n",
    "\n",
    "I2_mean = np.mean(I2_arr)\n",
    "I2_std = np.std(I2_arr)\n",
    "print(\"Gaussian: I = %g +- %g; frac err = %g; error estimate = %g\" % (I2_mean, I2_std, I2_std/I2_mean, err2))\n",
    "\n",
    "counts, bins = np.histogram(I_arr, bins=10, density = True)\n",
    "plt.stairs(counts, bins)    \n",
    "counts, bins = np.histogram(I2_arr, bins=10, density = True)\n",
    "plt.stairs(counts, bins)\n",
    "plt.plot([I0, I0], (0,1.05*max(counts)), \":\")\n",
    "plt.ylim((0,1.05*max(counts)))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2344769",
   "metadata": {},
   "source": [
    "Notes:\n",
    "\n",
    "- if you try different $N$ values, you'll see that the fractional error scales $\\propto 1/\\sqrt{N}$.\n",
    "\n",
    "- since our estimate for the integral is the mean value of the $N$ samples of $f(x)/p(x)$, another way to estimate the error in the integration is to calculate the error in the mean (or ``standard error''),\n",
    "\n",
    "$$\\sigma_I^2 = {1\\over N} \\mathrm{Var}(f) = {1\\over N} \\left( \\langle f^2\\rangle - \\langle f\\rangle^2\\right).$$\n",
    "\n",
    "The function `integrate_MC` in the code above does this by returning $\\sigma_I$ as `np.std(f(x)/p(x))/np.sqrt(N)`\n",
    "\n",
    "The values agree well with the standard deviation calculated from the histograms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "661d095e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}