{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise N: Numerical Integration\n",
    "\n",
    "## Estimated time: 20 minutes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numerical integration is very useful. Although it won't replace the tedium of cranking out your IA integrals, it can come in handy for tackling integration over data you collect from experiment, or for computing the values of otherwise analytically intractable integrals.\n",
    "\n",
    "There are many ways of doing this. As an introduction to the area, let's look at the Trapezium rule - which you might be familiar with - and Simpson's rule - which you may not be."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Demonstration: Trapezium Rule"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The Trapezium rule is one of the simplest methods for numerical integration.\n",
    "\n",
    "It works by splitting the area underneath a curve into multiple trapezia.\n",
    "\n",
    "If we want $N$ trapezia between $a$ and $b$, we can split this domain into steps of,\n",
    "\n",
    "\\begin{equation}\n",
    "    \\Delta x = \\frac{b - a}{N}\n",
    "\\end{equation}\n",
    "\n",
    "By considering the area of a trapezium, we can approximate the integral in the following way,\n",
    "\n",
    "\\begin{equation}\n",
    "    \\int^b_a f(x)dx \\approx \\frac{1}{2}\\Delta x \\Big[f(x_0) + f(x_N) + 2(f(x_1) + f(x_{2}) + ... f(x_{N-1}) \\big]\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Now, let's do this in python."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "We'll need to generate a *linspace*, so let's import *numpy*. We can also use this to call functions like *sin*, or the value of $\\pi$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true,
    "hidden": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The function $f(x)$ is the one we'll integrate - feel free to play around with it.\n",
    "To start, here's a simple *sin* function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true,
    "hidden": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return np.sin(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Now, we need to set the values of $a$, $b$, and $N$. From these we will calculate the spacing between x points $\\Delta x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false,
    "hidden": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "a = 0\n",
    "b = np.pi\n",
    "N = 10 ** 6\n",
    "delta_x =  (b - a)/N"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Using *numpy* we generate a linspace. The keyword argument *endpoint* is set to false to let numpy know that we wish to terminate just before we reach the value of $b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false,
    "hidden": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "x = np.linspace(a, b, num=N, endpoint=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "We'll store the value of our integration in *result*. We'll do it in python in the following way:\n",
    "\n",
    "1. Store the initial and final values of our function at $f(x_0)$ and $f(x_N)$ in *result*.\n",
    "2. *Iterate* through the section of our $x$ array between those points. This is done in python by using the slice of the array [1:-1] (between the first value up to the last).\n",
    "3. During the iteration, we'll add two times the value of the function $f(x_i)$ to the result using the *+=* operator.\n",
    "4. Finally, we print the result, multiplied by $\\frac{\\Delta x}{2}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false,
    "hidden": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.99999999999\n"
     ]
    }
   ],
   "source": [
    "result = f(x[0]) + f(x[-1] + delta_x)\n",
    "for x_i in x[1:-1]:\n",
    "    result += 2*f(x_i)\n",
    "print(0.5 * delta_x * result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Does this result match what you expect? Why, why not?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Task: Simpson's Rule"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Simpson's rule is another way of approximating a definite integral.\n",
    "\n",
    "Instead of approximating the area under a curve by a trapezium, Simpson's rule approximates the area using a *quadratic interpolant*.\n",
    "\n",
    "![caption](https://upload.wikimedia.org/wikipedia/commons/thumb/c/ca/Simpsons_method_illustration.svg/266px-Simpsons_method_illustration.svg.png)\n",
    "\n",
    "Image from Wikipedia. The curve $P(x)$ is the quadratic interpolant.\n",
    "\n",
    "https://en.wikipedia.org/wiki/Simpson%27s_rule\n",
    "\n",
    "We're going to use it to calculate the confidence intervals from a Gaussian."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Define the function $g(x) = \\frac{e^{-\\frac{t^2}{2}}}{\\sqrt{\\pi}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def g(x):\n",
    "    ### BEGIN SOLUTION\n",
    "    sdev = 1\n",
    "    return np.exp(-x**2 / (2*sdev**2)) / np.sqrt(2*np.pi*sdev**2)\n",
    "    ### END SOLUTION"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Now, define the function $simpson$. Make sure it  can take the following arguments (in this order): $func$, $a$, $b$, $N$.\n",
    "\n",
    "Where $func$ is a function to be passed to the integrator (not necessarily $g(x)$). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def simpson(func, a, b, N):\n",
    "    ### BEGIN SOLUTION\n",
    "    # Hmm, seems easiest to implement this non-pythonically, will possibly refactor in future.\n",
    "    # Also, I appreciate that a Python solution is already on wikipedia and I don't want to plagairise.\n",
    "    # Need to check that N is even, or else the method is invalid.\n",
    "    if N % 2 is not 0:\n",
    "        raise ValueError(\"N must be an even number.\")\n",
    "    \n",
    "    delta_x = (b - a)/N\n",
    "    result = func(a) + func(b)\n",
    "    \n",
    "    # Sum 1\n",
    "    j = 1\n",
    "    while j <= int(N/2) - 1:\n",
    "        result += 2 * func(a + 2*j*delta_x) \n",
    "        j += 1\n",
    "        \n",
    "    # Sum 2\n",
    "    j = 1\n",
    "    while j <= int(N/2):\n",
    "        result += 4 * func(a + (2*j - 1)*delta_x)\n",
    "        j += 1\n",
    "    \n",
    "    return (delta_x/3) * result\n",
    "    \n",
    "    ### END SOLUTION"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Use the *simpson* function to find the $1\\sigma$, $2\\sigma$, $3\\sigma$, $5\\sigma$, and $6\\sigma$ probabilities. In the markdown cell below, write down the chances in terms of fractions (e.g. 1 in a million)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Print out the values for the confidence intervals here..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "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.5.2"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}