{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Comparing the error rates for Monte Carlo and Riemann Integration\n",
    "A commonly mentioned advantage of Monte Carlo integration is that the error rate does not depend on the dimensionality of the input space; it only depends on the number of samples used to estimate the integral. Let $\\int_{\\mathcal{H}} f(x)$ be the integral we want to estimate. In Monte Carlo integration, we use samples drawn uniformly from $\\mathcal{H}$ to form our estimate as follows\n",
    "$$F_{mc} = \\frac{1}{N}(\\sum_i f(x_i))$$\n",
    "If we let $F$ to denote the true value of the integral, by the law of large numbers, $F_{mc}$ approaches $F$. Moreover, the variance of our estimate is $\\frac{\\sigma^2}{N}$ where $\\sigma^2$ is the variance of $f(x)$. Hence, the error of Monte Carlo estimate varies with $O(n^{-1/2})$.\n",
    "\n",
    "In Riemann integration, we simply partition $\\mathcal{H}$ into $N$ equally spaced intervals and estimate the integral simply by\n",
    "$$F_{r} = \\frac{1}{N} \\sum_i f(x_i) |H|$$\n",
    "where $|H|$ is the size of the domain of integration. The error of Riemann estimate varies with $O(n^{-1})$, i.e., a better rate compared to Monte Carlo integration.\n",
    "\n",
    "Let us first see that this is the case with a simple example. Let $f(x) = (1 - ||x||_2^2)$ where $x \\in \\mathbb{R}^D$. We use Monte Carlo and Riemann integration to estimate $\\int_0^1 f(x)$ and look at how the error changes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.01851852  0.01111111  0.00641026  0.0037037   0.00219298  0.00130208\n",
      "  0.00077882  0.00046555  0.00027871  0.00016683]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def grid(xl, xu, N, D):\n",
    "    \"\"\"\n",
    "    Create a grid of N evenly spaced points in D-dimensional space\n",
    "    xl: lower bound of x\n",
    "    xu: upper bound of x\n",
    "    N: number of points per dimension\n",
    "    D: number of dimensions\n",
    "    \"\"\"\n",
    "    xr = np.linspace(xl, xu, N)\n",
    "    g = np.zeros((N**D, D))\n",
    "    for n in range(N**D):\n",
    "        index = np.unravel_index(n, tuple([N]*D))\n",
    "        g[n] = [xr[i] for i in index]\n",
    "        \n",
    "    return g\n",
    "\n",
    "def f(x):\n",
    "    return (1 - (np.sum(np.square(x)) / x.size))\n",
    "\n",
    "def riemann(N, D):\n",
    "    # riemann integration\n",
    "    x = grid(0.0, 1.0, N, D)\n",
    "    dx = 1.0 / (N**D)\n",
    "    F_r = np.sum(np.apply_along_axis(f, 1, x) * dx)\n",
    "    return F_r\n",
    "\n",
    "def monte_carlo(N, D):\n",
    "    # monte carlo integration\n",
    "    x = np.random.rand(N**D, D)\n",
    "    F_mc = np.sum(np.apply_along_axis(f, 1, x)) / (N**D)\n",
    "    return F_mc\n",
    "\n",
    "D = 1\n",
    "N = np.logspace(1, 3, num=10, dtype=int)\n",
    "\n",
    "F_r = np.zeros(10)\n",
    "for i,n in enumerate(N):\n",
    "    F_r[i] = riemann(n, D)\n",
    "    \n",
    "# error in riemann estimate\n",
    "e_r = np.abs(F_r - 2.0/3.0)\n",
    "print(e_r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we plot the error in Riemann estimate with respect to number of sample points. As it is clear from the figure, this error varies with $O(N^{-1})$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(N, e_r)\n",
    "plt.plot(N, (e_r[0]*N[0])/N)\n",
    "plt.legend(['Error in Riemann estimate', '1/N'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.0740192   0.06172955  0.04757969  0.03471028  0.02688691  0.02152748\n",
      "  0.01563089  0.01264551  0.00963057  0.00727757]\n"
     ]
    }
   ],
   "source": [
    "repeats = 1000\n",
    "e_mc = np.zeros(10)\n",
    "for i,n in enumerate(N):\n",
    "    for r in range(repeats):\n",
    "        e_mc[i] += np.abs(monte_carlo(n, D) - 2.0/3.0)\n",
    "\n",
    "e_mc /= repeats\n",
    "print(e_mc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let us look at the error in the Monte Carlo estimate. Again, as we have seen above, the error in Monte Carlo estimate varies with $O(n^{-1/2})$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(N, e_mc)\n",
    "plt.plot(N, (e_mc[0]*np.sqrt(N[0]))/np.sqrt(N))\n",
    "plt.legend(['Error in Monte Carlo estimate', r\"$\\frac{1}{N^{-1/2}}$\"])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, why would anyone use Monte Carlo integration? Obviously, Riemann estimate decreases much more quickly as we increase the number of sample points. What is the deal with the error rate in Monte Carlo integration being independent of the number of dimensions? Do the Monte Carlo and Riemann errors change differently as we increase the number of dimensions? Let us look into that.\n",
    "\n",
    "First, we need to understand how the Riemann error changes if we go to a higher dimensional space. Let $N$ denote the number of points along each axis; hence, if we have $D$ dimensions, we have $N^D$ sample points. In a D-dimensional space, does the error change as $O(1/N)$ or $O(1/(N^D))$? Let us try $D=2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "D = 2\n",
    "N = np.array([5, 10, 20 , 50, 100])\n",
    "\n",
    "F_rd = np.zeros(5)\n",
    "for i,n in enumerate(N):\n",
    "    F_rd[i] = riemann(n, D)\n",
    "    \n",
    "e_rd = np.abs(F_rd - 2.0/3.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(N, e_rd)\n",
    "plt.plot(N, (e_rd[0]*N[0])/(N))\n",
    "plt.plot(N, (e_rd[0]*(N[0]**D))/(N**D))\n",
    "plt.legend(['Error in Riemann estimate', '1/N', '1/N^D'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks more like $O(1/N)$. Let us look at $D=3$ too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "D = 3\n",
    "N = np.array([3, 6, 12 , 25, 50])\n",
    "\n",
    "F_rd = np.zeros(5)\n",
    "for i,n in enumerate(N):\n",
    "    F_rd[i] = riemann(n, D)\n",
    "    \n",
    "e_rd = np.abs(F_rd - 2.0/3.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(N, e_rd)\n",
    "plt.plot(N, (e_rd[0]*N[0])/(N))\n",
    "plt.plot(N, (e_rd[0]*(N[0]**D))/(N**D))\n",
    "plt.legend(['Error in Riemann estimate', '1/N', '1/N^D'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again it looks more like $O(1/N)$. Another way to look at that is to see how the error changes as we change the number of dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "D = [1,2,3,4,5]\n",
    "N = 10\n",
    "\n",
    "F_rn = np.zeros(5)\n",
    "for i,d in enumerate(D):\n",
    "    F_rn[i] = riemann(N, d)\n",
    "    \n",
    "e_rn = np.abs(F_rn - 2.0/3.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(D, e_rn)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As it is clear, if we keep the number of points per dimension constant, the error stays constant. In other words, Riemann error behaves as $O(1/N)$. That also means that if we keep the number of sample (grid) points constant, Riemann error should increase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "D = np.array([1, 2, 4, 8, 16])\n",
    "# we need to keep the number of samples constant across dimensions\n",
    "N = np.array([65536, 256, 16, 4, 2])\n",
    "F_rc = np.zeros(5)\n",
    "for i,(n,d) in enumerate(zip(N,D)):\n",
    "    F_rc[i] = riemann(n, d)\n",
    "    \n",
    "# error in riemann estimate\n",
    "e_rc = np.abs(F_rc - (2.0 / 3.0))\n",
    "print(e_rc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(D, e_rc)\n",
    "plt.plot(D, (e_rc[0]*N[0])/N)\n",
    "plt.legend(['Error in Riemann estimate', '1/N'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Though the increase seems to be larger than $O(1/N)$.\n",
    "\n",
    "Nevertheless, for us the important point is to compare how the Monte Carlo error changes to Riemann error. Let us look at how the Monte Carlo error changes with respect to $N$ for $D=2$. Remember that $N$ is the number of points per dimension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.03473115  0.01715977  0.00924366  0.00337691  0.00177935]\n"
     ]
    }
   ],
   "source": [
    "repeats = 100\n",
    "D = 2\n",
    "N = np.array([5, 10, 20 , 50, 100])\n",
    "e_mcd2 = np.zeros(5)\n",
    "for i,n in enumerate(N):\n",
    "    for r in range(repeats):\n",
    "        e_mcd2[i] += np.abs(monte_carlo(n, D) - (2.0/3.0))\n",
    "\n",
    "e_mcd2 /= repeats\n",
    "print(e_mcd2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(N, e_mcd2)\n",
    "plt.plot(N, (e_mcd2[0]*(N[0]))/N)\n",
    "plt.legend(['Error in Monte Carlo estimate', r\"$\\frac{1}{\\sqrt{N^D}}$\"])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ahah! Note that the Monte Carlo error changes as $O(1/\\sqrt{N^D})$ which is $O(1/N)$ in this case. In other words, Monte Carlo error depends on the total number of sample points, not on the number of points per dimension (as Riemann error does). For example, if we keep $N$ constant as we increase the number of dimensions, the Monte Carlo error should decrease. Remember, Riemann error stays the same as long as $N$ is constant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.08284523  0.01660181  0.00478395  0.00116964  0.000303  ]\n"
     ]
    }
   ],
   "source": [
    "repeats = 100\n",
    "D = np.array([1, 2, 3, 4, 5])\n",
    "# we need to keep the number of samples constant across dimensions\n",
    "N = 10\n",
    "e_mcd = np.zeros(5)\n",
    "for i,d in enumerate(D):\n",
    "    for r in range(repeats):\n",
    "        e_mcd[i] += np.abs(monte_carlo(N, d) - (2.0/3.0))\n",
    "\n",
    "e_mcd /= repeats\n",
    "print(e_mcd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(D, e_mcd)\n",
    "D = np.array([1, 2, 3, 4, 5])\n",
    "plt.plot(D, (e_mcd[0]*np.sqrt(N))/np.sqrt(N**D))\n",
    "plt.legend(['Error in Monte Carlo estimate', r\"$\\frac{1}{\\sqrt{N^D}}$\"])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that how the Monte Carlo error changes as $\\frac{1}{\\sqrt{N^D}}$. Therefore, one only needs a constant number of  sample points to achieve a constant error rate. In contrast, for Riemann integration, we need to keep the number of points per dimension constant, which means an exponentially increasing number of sample points as $D$ increases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To recap, for Monte Carlo integration, the error varies with $O(1/\\sqrt{N^D})$ where $N$ is the number of points per dimension. In contrast, for Riemann integration, the error varies with $O(1/N)$. Therefore, if we keep the total number of sample points constant, Riemann estimate will get worse as $D$ increases, while Monte Carlo error should stay the same. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}