Solving an Integration Bee Problem – Jason Bernstein’s Website

I recently learned about integration bees, which are integration contests often held for undergraduates. I wondered if I could solve any of these integrals, so I searched for problems online and found this integral, \[ I=\int_{-\infty}^0x^{2022}\exp(4x)\;dx, \] from an integration bee held at the University of New South Wales (UNSW) in 2022 (problem 4 in knockout match 1). In this post, I’ll share how to compute this integral in closed-form and numerically in Python.

The first thought I had when I saw this integral is to make the change of variables from \(x\) to \(-x\) and to switch the integration limits, so that the integral becomes \[ I=\int_0^{\infty}x^{2022}\exp(-4x)\;dx. \] It seems more natural to me to have the integrand going from zero to positive infinity instead of from negative infinity to zero. Then I thought the exponential term reminded me of the unnormalized probability density function (PDF) of an exponential random variable with scale parameter \(\lambda=4\), or mean 1/4, so I rewrote the integral to include the normalizing constant as \[ I=\frac{1}{4}\int_0^{\infty}x^{2022}(4\exp(-4x))\;dx. \] Ah, so if \(X\) is an exponential random variable with mean \(1/4\), then \[ I=\frac{1}{4}E[X^{2022}]. \] I spent a few minutes computing the 2022nd moment of an exponential random variable from its moment generating function, but eventually looked on the Wikipedia page for exponential random variables to find that \[ E[X^{2022}]=\frac{2022!}{4^{2022}}. \] Plugging this expected value back into our last expression for \(I\) gives \[ \boxed{I=\frac{2022!}{4^{2023}}}. \] Woohoo, integral solved using higher-order moments!

After some more thought, I realized another way to compute the integral is to identify the integrand as the kernel of a gamma distribution with shape parameter 2023 and rate parameter 4. It follows that \[ I=\frac{\Gamma(2023)}{4^{2023}}\left(\frac{4^{2023}}{\Gamma(2023)}\int_0^{\infty}x^{2022}\exp(-4x)\;dx\right). \] The term in parentheses equals one since the integral is of a gamma distribution over its support, so again we get our previous value of \(I\) since \(\Gamma(2023)=2022!\). The normalizing constant follows from the definition of the gamma function, \[ \Gamma(z)=\int_0^{\infty}t^{z-1}\exp(-t)\;dt \] for complex \(z\) with \(\textrm{Re}(z)>0\), after making the change of variables from \(t\) to \(4t\). I’ll guess that most people would compute \(I\) by recognizing it as a special case of the gamma function with \(z=2023\) and a change of variables.

We can also compute the integral numerically in Python using the math module. It is safe to use the floor division operator // here since 2022! is evenly divisible by \(4^{2023}\).

import math
import sys

# Evaluate integral
I = math.factorial(2022) // 4**2023

# Format integral in scientific notation
sys.set_int_max_str_digits(6_000)
s = str(I)
pretty_I = f"{s[0]}.{s[1:3]}e{len(s)-1}"

The integral is 1.70e4590, which is 4,591 digits long. That’s a big number!

The integral can also be evaluated using the SymPy package.

import sympy as sp

# Evaluate integral symbolically
x = sp.symbols("x")
I2_sym = sp.integrate(x**2022 * sp.exp(4 * x), (x, -sp.oo, 0))

# Convert symbolic expression to numerical result
I2_number = I2_sym.evalf()

In this code chunk, I2_sym is a fraction, it contains a / and has type sympy.core.numbers.Rational. Evaluating the fraction using the evalf() method gives us the expected result, 1.70E+4590.

Now, time for a small side quest with SymPy. The numerator and denominator in I2_sym are coprime, meaning they have no common factors greater than one. Let’s check this claim by factoring the numerator and denominator returned by SymPy and obtaining their largest common factor. If there are no common factors, then the two numbers are coprime.

# Get prime factors of numerator and denominator
num_factors = sp.factorint(I2_sym.numerator).keys()
den_factors = sp.factorint(I2_sym.denominator).keys()

# Confirm numerator and denominator have no common factors
common_factors = set(num_factors) & set(den_factors)
if not common_factors:
    print("coprime!")

coprime!

Can we compute this integral using a method covered in a first course on integration like integration by parts? Yes! The Wikipedia page on integration by parts has this formula for repeated integration by parts, \[ \int_0^{\infty} u^{(0)}v^{(n)}\;dx=\sum_{k=0}^{n-1}\left.(-1)^nu^{(k)}v^{(n-1-k)}\right|_0^{\infty}+(-1)^n\int_0^{\infty} u^{(n)}v^{(0)}\;dx. \] If we take \(n=2022\), \(u^{(0)}=x^{2022}\), and \(v^{(n)}=\exp(-4x)\), then all the summands vanish, \(u^{(n)}=2022!\), and \(v^{(0)}=4^{-2022}\exp(-4x)\), so that \[ \int_0^{\infty} x^{2022}\exp(-4x)\;dx=\frac{2022!}{4^{2022}}\int_0^{\infty}\exp(-4x)\;dx. \] The integral on the right hand side equals 1/4, so this approach results in the previously obtained value of the integral of interest.

Overall, we computed the integral:

As an expectation with respect to an exponential distribution
As the normalizing constant of a gamma distribution
Using the math module in Python
Using the SymPy module in Python
Using repeated integration by parts

When I started writing this post I thought I would also estimate the integral using Monte Carlo by simulating \(X_1,...,X_n\stackrel{iid}{\sim}\textrm{Exp}(4)\) and then computing \(\hat{I}=n^{-1}\sum_i X_i^{2022}\). A quick simulation shows this leads to overflow errors though, which is not surprising considering we are raising numbers greater than one to the 2022nd power.

import numpy as np

np.random.seed(12)
x = np.random.exponential(size=3, scale=4)
x**2022

array([0.00000000e+000,             inf, 2.14161603e+176])

I would like to add this Monte Carlo approach to this blog post eventually.