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.
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.
The Trapezium rule is one of the simplest methods for numerical integration.
It works by splitting the area underneath a curve into multiple trapezia.
If we want $N$ trapezia between $a$ and $b$, we can split this domain into steps of,
\begin{equation} \Delta x = \frac{b - a}{N} \end{equation}By considering the area of a trapezium, we can approximate the integral in the following way,
\begin{equation} \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] \end{equation}Now, let's do this in python.
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$.
import numpy as np
The function $f(x)$ is the one we'll integrate - feel free to play around with it. To start, here's a simple sin function.
def f(x):
return np.sin(x)
Now, we need to set the values of $a$, $b$, and $N$. From these we will calculate the spacing between x points $\Delta x$.
a = 0
b = np.pi
N = 10 ** 6
delta_x = (b - a)/N
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$.
x = np.linspace(a, b, num=N, endpoint=False)
We'll store the value of our integration in result. We'll do it in python in the following way:
result = f(x[0]) + f(x[-1] + delta_x)
for x_i in x[1:-1]:
result += 2*f(x_i)
print(0.5 * delta_x * result)
1.99999999999
Does this result match what you expect? Why, why not?
Simpson's rule is another way of approximating a definite integral.
Instead of approximating the area under a curve by a trapezium, Simpson's rule approximates the area using a quadratic interpolant.
Image from Wikipedia. The curve $P(x)$ is the quadratic interpolant.
https://en.wikipedia.org/wiki/Simpson%27s_rule
We're going to use it to calculate the confidence intervals from a Gaussian.
Define the function $g(x) = \frac{e^{-\frac{t^2}{2}}}{\sqrt{\pi}}$
def g(x):
### BEGIN SOLUTION
sdev = 1
return np.exp(-x**2 / (2*sdev**2)) / np.sqrt(2*np.pi*sdev**2)
### END SOLUTION
Now, define the function $simpson$. Make sure it can take the following arguments (in this order): $func$, $a$, $b$, $N$.
Where $func$ is a function to be passed to the integrator (not necessarily $g(x)$).
def simpson(func, a, b, N):
### BEGIN SOLUTION
# Hmm, seems easiest to implement this non-pythonically, will possibly refactor in future.
# Also, I appreciate that a Python solution is already on wikipedia and I don't want to plagairise.
# Need to check that N is even, or else the method is invalid.
if N % 2 is not 0:
raise ValueError("N must be an even number.")
delta_x = (b - a)/N
result = func(a) + func(b)
# Sum 1
j = 1
while j <= int(N/2) - 1:
result += 2 * func(a + 2*j*delta_x)
j += 1
# Sum 2
j = 1
while j <= int(N/2):
result += 4 * func(a + (2*j - 1)*delta_x)
j += 1
return (delta_x/3) * result
### END SOLUTION
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).
# Print out the values for the confidence intervals here...
Answer.