The world’s leading publication for data science, AI, and ML professionals.

A Simple Method for Numerical Integration in Python

In this article, we will introduce a simple method for computing integrals in python. We will first derive the integration formula and…

In this article, we will introduce a simple method for computing integrals in python. We will first derive the integration formula and then implement it on a few functions in python. This article assumes you have a basic understanding of probability and integral calculus, but if you don’t you can always skip ahead to the examples. Enjoy!

Derivation

We must first state the definition of the expected value of a continuous random variable.

Suppose X is a random variable with with probability density function f(x). The expected value of X is defined as follows:

Expectation of X
Expectation of X

Next, we use the expectation formula to derive a simple equation for computing an integral. We would like to estimate the following integral:

Integral from a to b of a function g(x)
Integral from a to b of a function g(x)

We first rewrite the integral as follows:

We can then define a function h(x) as:

This allows us to rewrite the integral in a familiar form:

All of the computation in the integral has been reduced down to an expectation, and we know how to find the expected value of a set of data. The final approximation becomes:

Python Example 1: Integrating an Elementary Function

We will start simple by integrating the quadratic function f(x) = x² from 0 to 1.

# Dependencies
import numpy as np
import scipy.integrate as integrate
import matplotlib.pyplot as plt
# Our integral approximation function
def integral_approximation(f, a, b):
    return (b-a)*np.mean(f)
# Integrate f(x) = x^2
def f1(x):
    return x**2
# Define bounds of integral
a = 0
b = 1
# Generate function values
x_range = np.arange(a,b+0.0001,.0001)
fx = f1(x_range)
# Approximate integral
approx = integral_approximation(fx,a,b)
approx

Our integral approximation comes out to be:

This is about what we would expect since the true value of the integral is 1/3. However, we can also compare our result to Scipy’s "quad" function.

# Scipy approximation
integrate.quad(f1,a,b)

The resulting value:

I’m sure this example has you on the edge of your seat, but let’s see if we can’t integrate a more complicated function.

Python Example 2: Integrating the Gaussian Function

The gaussian function is notorious for being extremely difficult to integrate. In this example, we will put our method to the test by integrating the standard normal distribution.

# Integrating a random function
np.random.seed(42)
def gaussian(x, mu, sigma):
    return np.exp((-1*(x - mu)**2) / (2 * sigma**2))
# Define mu and sigma
mu = 0
sigma = 1
# Define bounds of integral
a = -3
b = 3
# Generate function values
x_range = np.linspace(-3,3,200)
fx = gaussian(x_range, mu, sigma)

The resulting function looks like this:

The nice thing about our integral approximation is that the complexity of the function does affect the difficulty of the computation. In every case, all we need is the bounds of integration and the function values.

# Our approximation
approx = integral_approximation(fx,a,b)
approx

Comparing to the Scipy solution:

# Scipy approximation
integrate.quad(lambda x: np.exp((-1*(x - mu)**2) / (2 * sigma**2)),a,b)

I hope that you found this article easy to follow and interesting!

Thank you!

Like my articles? Buy me a coffee: https://www.buymeacoffee.com/HarrisonfhU


Related Articles