phys512/montecarlo_solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6d46a72f",
   "metadata": {},
   "source": [
    "# Monte Carlo Integration"
   ]
  },
  {
   "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": 2,
   "id": "5582478b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrate_MC(f, N, sampling = 'uniform'):\n",
    "\n",
    "    if sampling == 'gaussian':\n",
    "        p = lambda x: np.exp(-x**2/2)/(np.sqrt(2*np.pi))\n",
    "        x = rng.normal(size = N)    # Generate x values between +-xmax\n",
    "    else:\n",
    "        xmax = 3\n",
    "        p = lambda x: np.ones_like(x) / (2*xmax)\n",
    "        x = xmax*(-1 + 2*rng.uniform(size = N))\n",
    "    \n",
    "    return np.sum(f(x)/p(x)) / N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5c39e032",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integral from scipy.quad = 1.78596 with error 5.49113e-09\n",
      "Uniform: I = 1.78656 +- 0.0220357; frac err = 0.0123341\n",
      "Gaussian: I = 1.78585 +- 0.00960208; frac err = 0.00537675\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAus0lEQVR4nO3de3hU1b3/8c/cMplwSSA0XIMHpIgVEQTBQJUqUapUpWCLVhCVB1sbrYitlFq11QpeWqXWgJZjUY9QECpWgeJBquCFmyj94aVQqpYUJCglBHKZzGX9/kid0wgoQ/Zkr2S/X88zjys7O2u+yzCTz6y9Zo3PGGMEAABgMb/bBQAAAHwRAgsAALAegQUAAFiPwAIAAKxHYAEAANYjsAAAAOsRWAAAgPUILAAAwHpBtwv4rGQyqd27d6tNmzby+XxulwMAAI6BMUYHDx5Uly5d5Pc7Px9iXWDZvXu3CgsL3S4DAAAch7KyMnXr1s3xfq0LLG3atJFUP+C2bdu6XA0AADgWlZWVKiwsTP0dd5p1geXTy0Bt27YlsACQJMUSMT313lOSpPEnj1coEHK5IgBHk6nlHNYFFgD4rFgypgc2PyBJGnfSOAIL4EEEFgDWC/qDuvjEi1NtAN7DIx+A9bICWbr7q3e7XQYAF7EPCwAAsB6BBQAAWI/AAsB61bFqDV0wVEMXDFV1rNrtcgC4gDUsAJqFg7GDbpcAwEUEFgDWyw5ma9k3l6XaALyHwALAen6fXye0PcHtMgC4iDUsAADAesywALBeLBnTku1LJEmX9r5UIT873QJeQ2ABYL1YIqYZG2ZIki458RICC+BBBBYA1gv4AzrvhPNSbQDeQ2ABYJeKMql6X4NDYUkP9L6y/ovy99LvMydfyitsfG0AXENgAWCPijKpdLDk9OZwoRypZCOhBWjGCCwA7FG9rz6sjJkrdejtTJ+fbJeemVzfN4EFaLYILADs06G31KV/6suaeI2+8cw3JEnLxixTJBhxqTAAbiGwALCeMUZ7a/am2gC8h8ACwHrhQFiLL1qcagPwHgILAOsF/AH1ad/H7TIAuIit+QEAgPWYYQFgvVgypuXvL5ckjeo5ip1uAQ8isACwXiwR022v3SZJOv+E8wksgAcRWABYL+AP6KyuZ6XaALyHwALAeuFAWLOLZ7tdBgAXsegWAABYj8ACAACsR2ABYL2aeI1GPTNKo54ZpZp4jdvlAHABa1gAWM8Yo50Hd6baALyHwALAeuFAWE9e8GSqDcB7CCwArBfwBzSgYIDbZQBwEWtYAACA9ZhhAWC9eDKu1TtXS5JGdB+hoJ+nLsBreNQDsF5dok4/XPNDSdKG72wgsAAe1KhLQvfcc498Pp+mTJmSOlZbW6uSkhLl5+erdevWGjt2rMrLyxtbJwAP8/v8GtRxkAZ1HCS/jyvZgBcd98uUTZs26dFHH1W/fv0aHL/pppu0fPlyLV68WLm5ubr++us1ZswYvfbaa40uFoA3ZQezNe/r89wuA4CLjuulyqFDh3TFFVdo7ty5ateuXer4gQMH9Nhjj+mBBx7Queeeq4EDB2revHl6/fXXtX79eseKBgAA3nJcgaWkpESjRo1ScXFxg+ObN29WLBZrcLxPnz7q3r271q1bd8S+otGoKisrG9wAAAD+U9qXhBYuXKg333xTmzZtOux7e/bsUVZWlvLy8hoc79ixo/bs2XPE/mbOnKmf//zn6ZYBwENq47Uav2K8JOmpC59SdjDb5YoANLW0ZljKysp04403av78+crOduYJY/r06Tpw4EDqVlZW5ki/AFqOpElq2/5t2rZ/m5Im6XY5AFyQ1gzL5s2btXfvXp1++umpY4lEQmvXrtXDDz+sF154QXV1daqoqGgwy1JeXq5OnTodsc9wOKxwmK22ARxdOBDWo+c9mmoD8J60AsuIESO0devWBseuvvpq9enTR9OmTVNhYaFCoZBWr16tsWPHSpK2bdumnTt3qqioyLmqAXhKwB/Q0C5D3S4DgIvSCixt2rRR3759Gxxr1aqV8vPzU8cnTZqkqVOnqn379mrbtq1uuOEGFRUV6cwzz3SuagAA4CmObxf54IMPyu/3a+zYsYpGoxo5cqRmz57t9N0A8JB4Mq7Xd78uSRraZSg73QIe1OhH/csvv9zg6+zsbJWWlqq0tLSxXQOApPqt+UtWl0hia37Aq3jUA7Ce3+fXKfmnpNoAvIfAAsB62cFsLfzGQrfLAOAiXqoAAADrEVgAAID1CCwArFcbr9WEFRM0YcUE1cZr3S4HgAtYwwLAekmT1JaPt6TaALyHwALAelmBLM06Z1aqDcB7CCwArBf0BzWi+wi3ywDgItawAAAA6zHDAsB6iWRCb+59U5J0esHpCvgDLlcEoKkRWABYL5qI6poXrpFUvzV/jj/H5YoANDUCCwDr+Xw+nZh7YqoNwHsILACsFwlG9OzoZ90uA4CLWHQLAACsR2ABAADWI7AAsF5tvFaT/3eyJv/vZLbmBzyKNSwArJc0Sa3/aH2qDcB7CCwArJcVyNLMs2am2gC8h8ACwHpBf1Df6PkNt8sA4CLWsAAAAOsxwwLAeolkQu/96z1J0sntT2ZrfsCDCCwArBdNRHX58sslsTU/4FUEFgDW8/l86tKqS6oNwHsILACsFwlG9MKlL7hdBgAXsegWAABYj8ACAACsR2ABYL1oIqof/PkH+sGff6BoIup2OQBcwBoWANZLJBN6qeylVFu8qxnwHAILAOuFAiHdUXRHqg3AewgsAKwX8od0ae9L3S4DgItYwwIAAKzHDAsA6yVNUu9XvC9J6pnXU34fr7UAryGwALBebbxW33zum5L+vTV/iK35Aa8hsABoFtqF27ldAgAXEVgAWC8nlKO1l611uwwALuJCMAAAsB6BBQAAWI/AAsB60URU09ZO07S109iaH/AoAgsA6yWSCa34YIVWfLCifmt+AJ7DolsA1gsFQrrljFtSbQDeQ2ABYL2QP6QJX5ngdhkAXMQlIQAAYD1mWABYL2mS+qjqI0lS51ad2Zof8CACCwDr1cZr9fU/fF0SW/MDXkVgAdAsRIIRt0sA4CICCwDr5YRytPGKjW6XAcBFXAgGAADWI7AAAADrcUkIwPGrKJOq9znX3yfbj3i4LlGnGRtmSJJ+MuQnygpkOXefAJoFAguA41NRJpUOlmLVzvYbypFy8hsciifj+sPf/iBJuuWMWwgsgAcRWAAcn+p99WFlzFypQ2/n+s3Jl/IKGxwK+UO6YcANqTYA7yGwAGicDr2lLv0zehehQEjX9rs2o/cBwG4sugUAANZjhgWA9Ywx2h/dL0lqF24nn8/nckUAmhqBBYD1auI1Gr5ouCS25ge8iktCAADAesywALBeTihHWydudbsMAC5ihgUAAFiPwAIAAKzHJSEA1qtL1OnBzQ9Kkm4aeBM73QIexAwLAOvFk3E99d5Teuq9pxRPxt0uB4ALmGEBYL2QP6TJp05OtQF4D4EFgPVCgZB+cPoP3C4DgIu4JAQAAKzHDAsA6xljVBOvkSRFghG25gc8iBkWANariddoyIIhGrJgSCq4APAWAgsAALAel4QAWC8SjGjDdzak2gC8h8ACwHo+n49PaAY8jktCAADAesywALBeLBHTnL/MkSRdd9p1CgXYPA7wmrRmWObMmaN+/fqpbdu2atu2rYqKivSnP/0p9f3a2lqVlJQoPz9frVu31tixY1VeXu540QC8JZaMae7WuZq7da5iyZjb5QBwQVqBpVu3brrnnnu0efNmvfHGGzr33HN1ySWX6J133pEk3XTTTXr++ee1ePFirVmzRrt379aYMWMyUjgA7wj6gxp/8niNP3m8gn4mhgEvSuuRf9FFFzX4+u6779acOXO0fv16devWTY899pgWLFigc889V5I0b948nXzyyVq/fr3OPPNM56oG4ClZgSxNGzzN7TIAuOi4F90mEgktXLhQVVVVKioq0ubNmxWLxVRcXJw6p0+fPurevbvWrVt31H6i0agqKysb3AAAAP5T2oFl69atat26tcLhsL73ve9p6dKl+spXvqI9e/YoKytLeXl5Dc7
      "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 = scipy.integrate.quad(f, -np.inf, np.inf)\n",
    "print(\"Integral from scipy.quad = %g with error %g\" % (I0, err))\n",
    "\n",
    "N = 10**4\n",
    "\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] = integrate_MC(f, N)\n",
    "    I2_arr[i] = 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\" % (I_mean, I_std, I_std/I_mean))\n",
    "\n",
    "I2_mean = np.mean(I2_arr)\n",
    "I2_std = np.std(I2_arr)\n",
    "print(\"Gaussian: I = %g +- %g; frac err = %g\" % (I2_mean, I2_std, I2_std/I2_mean))\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": "code",
   "execution_count": null,
   "id": "e3558622",
   "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
}