Lecture12.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 12:  Numerical Integration and Monte Carlo\n",
    "\n",
    "**Overview:**\n",
    "* Numerical integration.\n",
    "* Simple Monte Carlo integration.\n",
    "* Importance sampling.\n",
    "\n",
    "**Next Lecture:**\n",
    "* Markov chain Monte Carlo and the Metropolis Algorithm.\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#matplotlib notebook\n",
    "import numpy as np\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import integrate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dart Board Estimate of $\\pi$\n",
    "\n",
    "The code in the cell below generates a set of random coordinates inside our unit square and calculates the magnitude of the vector defined by these coordinates.\n",
    "\n",
    "* Run the code in the cell and call out the number generated for Prof. Plumb to plot on the board."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = 2 * (np.random.random([1, 2])) - 1\n",
    "print(np.sqrt(np.sum(x**2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functions to integrate, and some exact results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gaussian(x):\n",
    "    return np.exp(-x**2)\n",
    "\n",
    "# A Function that is not well behaved \n",
    "def Fermi(x):\n",
    "    num = 1/np.sqrt(x)\n",
    "    den = np.exp(x)+1\n",
    "    return num/den\n",
    "\n",
    "def fp(x):\n",
    "    arg = 1/(x*(2-x))\n",
    "    return (np.sin(arg))**2\n",
    "\n",
    "# area of a unit circle\n",
    "def sphere(x): \n",
    "    \"\"\" \n",
    "    return 1 if point is inside radius, zero otherwise\n",
    "    x is a multidimensional vector, must have dimension greater than 1\n",
    "    \"\"\"\n",
    "    r = np.sum(x**2, axis = 1) \n",
    "    a = (r<=1).astype(int)\n",
    "    return 1.0*a\n",
    "\n",
    "# volume of a hypersphere in n dimensions\n",
    "hypersphere = lambda r, n: math.pi**(n / 2)/math.gamma(n / 2 + 1)*r**n\n",
    "\n",
    "# exact integral of a Gaussian\n",
    "analyticalIntegral = np.sqrt(np.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Numerically integrate a Gaussian \n",
    "\n",
    "def riemannSum(f):\n",
    "    width = 10.0\n",
    "    n = 100\n",
    "    dx = width/n\n",
    "    x = np.arange(-0.5*width, 0.5*width, width/n)\n",
    "\n",
    "    return np.sum(f(x)) * dx\n",
    "\n",
    "\n",
    "def Simpson(f):\n",
    "    width = 10.0\n",
    "    n = 100\n",
    "    dx = width/n\n",
    "    x = np.arange(-0.5*width, 0.5*width, width/n)\n",
    "    \n",
    "    s = (f(-0.5*width)+f(0.5*width)) \n",
    "    return (2*f(x[2:-2:2]).sum() + 4*f(x[1:-2:2]).sum()+s) * dx/3\n",
    "\n",
    "# using Scipy's built in integration schemes\n",
    "scipyIntegral = integrate.quad(gaussian, -100.0, 100.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#np.random.seed(256)\n",
    "def naiveMonteCarlo(f, limits = [-10,10],d = 1, n_points = 1000, NSamples =100):\n",
    "    \"\"\"\n",
    "    Implement a mean value Monte-Carlo Integration in d dimensions\n",
    "    \n",
    "    f is function to integrate, must take an input vector x of dimension d\n",
    "    \n",
    "    limits define the range of integration, this function only works for integration ranges that are the same\n",
    "    all dimensitons\n",
    "    \n",
    "    n_points are number of points to sample in domain\n",
    "    NSamples number of time to repeat integration, decrease statistical noise\n",
    "    \"\"\"\n",
    "    width = np.abs(limits[1] - limits[0])\n",
    "    samples = np.zeros(NSamples)\n",
    "    \n",
    "    for i in range(NSamples):\n",
    "        x = width * (np.random.random([n_points, d])) + limits[0]\n",
    "        samples[i] = width**d * np.sum(f(x))/n_points    \n",
    "    \n",
    "    return samples.mean(), samples.std()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Analytical (exact) integral = \", analyticalIntegral)\n",
    "print(\"Riemann Sum = \", riemannSum(gaussian))\n",
    "print(\"Simpsons Rule = \", Simpson(gaussian))\n",
    "print(\"SciPy Integral = \", scipyIntegral)\n",
    "print(\"Naive Mean Value Monte Carlo = \", naiveMonteCarlo(gaussian))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tasks\n",
    "* Use the mean value method Monte Carlo method to estimate the value of $\\pi$, (area of unit circle) to a higher accuracy than what was done in the demonstration.\n",
    "* In the cell below, we plot a Pathological function, first examine the plot. Now suppose we wished to determine the value of the integral of this function between x=0 and x=2. Can you evaluate this integral using Simpsons rule or Scipy? Do you expect these results to be accurate? What about using Monte Carlo integration? (the numerical value of the integral should be about 1.45."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "x = np.linspace(1e-8,2-1e-8,10000)\n",
    "ax.plot(x,fp(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### More Tasks\n",
    "* Use the  Monte Carlo integrator to n-dimensions to find the volume of a hypersphere in 5 dimensions. Compare this results with Simpsons rule and the exact value.\n",
    "* Can you confirm the error on the MC integration is independent of the number of dimensions?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importance sampling Monte Carlo\n",
    "* Review the importance sampling method below. Do you understand all of the steps?\n",
    "* Can you modify the method and integrate a 4 dimensional Gaussian function?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def p_normal(stdev, x): # normal distribution\n",
    "    s = 1.0 / stdev\n",
    "    s2 = s**2\n",
    "    return np.exp(-s2 * x**2) * s / np.sqrt(np.pi)\n",
    "\n",
    "def importanceSampledMonteCarlo(f, p, NSamples=10):\n",
    "    n = 10000\n",
    "    stdev = 1.0\n",
    "    samples = np.zeros(NSamples)\n",
    "\n",
    "    for i in range(NSamples):\n",
    "        # sample random values from a normal distribution\n",
    "        x = np.random.normal(loc = 0.0, scale = np.sqrt(0.5) * stdev, size = n)\n",
    "        samples[i] = (f(x) / p(stdev, x)).mean()       \n",
    "    \n",
    "    return samples.mean(), samples.std()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print(\"Naive Mean Value Monte Carlo = \", naiveMonteCarlo(gaussian))\n",
    "print(\"Importance Sampled Monte Carlo = \", importanceSampledMonteCarlo(gaussian, p_normal))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}