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+POTHSnIlJbg
      "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
}