Monte Carlo Integration: Improper Integrals

Photo by Angely Acevedo on Unsplash

Improper integrals are a subcategory of definite integrals in which one or both of the limits goes to infinity or negative infinity.

There is a range of fast and accurate techniques to solve improper integrals. But I am going to focus on using change-of-variables to map an infinite range onto traditional Monte Carlo integration.

Unbound integration limits pose a difficulty to traditional Monte Carlo integration because the technique relies on representative random sampling of the function within the integration limits. Sampling from an infinite distribution is troublesome for several reasons. First, the “important” part of a function with the highest integrand may span a small range while the limits are unbound. Therefore to get an accurate integrand estimate we would need an incredibly large number of samples between arbitrarily large positive and negative numbers to appropriately sample the “important” range. To make matters worse, if we set arbitrary upper and lower limits to our sampling we may miss out on important data pieces that lie outside of the limits.

This doesn’t mean that brute force Monte Carlo integration can’t be used for estimating improper integrals — it just suggests that naively using Monte Carlo might give unreliably or unpredictable results.

For instance, take the function:

As shown, the analytical solution is π, so we have something to aim for. We’ll first compute it assuming we do not know the function’s shape, so we need to arbitrarily select our integration limits. So let’s assume we take 30,000 samples from the region x=[0,5000]. The code below performs the naive Monte Carlo integration and results are printed below.

Estimate for 0.0 -> 5000.0 is 5.90, (30000 iterations)

The brute force method simply isn’t very accurate with the arbitrarily imposed limits. One possibility would be to look at the function’s shape on a graph and select a cutoff point visually. In this case, the function has the shape shown below:

So maybe a lower arbitrary cutoff distance would be better — let’s try using an upper cutoff of 50 and rerunning the same code above.

Estimate for 0.0 -> 50.0 is 3.10, (30000 iterations)

The lower integration limit certainly improves the estimate in this case. But we can’t always just look at a plot of a function to determine which limits to impose, especially when the functions are complex or the minima and maxima spread across a wide range of values. Also, this approach would require visual inspection for each new integral, which could be extremely time-consuming.

Change of Variables

Rather than trying to sample an infinite distribution, we can use a change of variables to sample a finite distribution that maps to an infinite distribution. If we want to sample a [-∞, ∞] distribution, we can find a distribution with a finite input range and an output range that spans [-∞, ∞]. One possible distribution is tan(x), shown below between x=-π/2 and x = π/2 radians.

Because the tan(x) distribution goes to negative infinity as it approaches -π/2 radians and positive infinity as it approaches π/2, we can sample an infinite distribution by sampling tan(x) where x ∈ [π/2 rad, -π/2 rad]. We can use the traditional change-of-variables approach for integrals:

In the case of tan(u), we can first find the derivative of our change of variable function, dx/du:

From which we can convert any integral with two unbound limits as:

To use this approach on the function on which that we previously used the brute force Monte Carlo approach, the integral could be transformed as:

We can then use normal brute force Monte Carlo integration on this transformed function to see how it performs. The code below uses the above expression to execute the integration:

Estimate for 0.0 -> 90.0 is 3.18, (30000 iterations)

Much better. We also didn’t have to make any assumptions about the shape of the function.

Many other functions have similar properties as tan(x) in terms of tending towards +/- infinity within a finite range. So finding the best limit-transforming function may require some knowledge of the field of application. Additionally, if a steeper gradient towards -∞/∞ is needed the original transformation function can be put to an odd power to preserve the limits, such as tan(x)³ or tan(x)⁵.

Thanks for reading — any suggestions or corrections are welcome! Please reach out to me at mediumCameron@gmail.com

References

1. Lisovskaja, Vera. Mathematical Statistics, Chalmers — TMS150, Stochastic Data Processing and Simulation, www.math.chalmers.se/Stat/Grundutb/CTH/tms150/1112/.
2. Jarosz, Wojciech. “Efficient Monte Carlo Methods for Light Transport in Scattering Media.” UC San Diego, 2008, cs.dartmouth.edu/~wjarosz/publications/dissertation/.
3. “Monte Carlo Method.” Wikipedia, Wikipedia, https://en.wikipedia.org/wiki/Monte_Carlo_method