If integration is summing many small piles, we have to figure out how big the piles are. Their height is usually given by a mathematical function, and our first example will be the same as in the Feynman trick article.

\[f(x) = \frac{x - 1}{\ln{x}}\]

This is to be integrated from zero to one, i.e. we want to know the size of the shaded area in the plot below. You can think of each column of shaded pixels as one pile, and we sum the size of all of them to get the total area.1¹ Of course, this is an svg image so there are no columns of pixels. Alternatively, the more we zoom in, the thinner the columns become – but the more of them there are. This is why we need integration: it’s dealing with the limit case of infinitely many, infinitely thin columns.

We could imagine drawing six random numbers between zero and one, and plotting piles of the corresponding height at those locations. Since there are six piles, their width is one sixth of the width of the area we are integrating.

Even though some of these piles overlap by chance, and even though there are some random gaps between them, the sum of their areas (0.66) comes very close to the actual shaded area determined analytically (0.69). If we draw more piles, we have to make them correspondingly thinner, but the agreement between their sum and the total size of the area improves.

These are 100× as many piles, and they’re 1/100th as thick to compensate. Their total area is 0.70 – very close to 0.69. If we draw even more piles, we’ll get even closer.

This illustrates a neat correspondence between integrals and expected values. In the simple case, we can frame it mathematically as

\[\int_a^b f(x) \mathrm{d}x = E(f(x))\]

In words, this says that integrating the function \(f\) between \(a\) and \(b\) is the same as taking the expected value of \(f(x)\) at uniformly distributed random points between \(a\) and \(b\).

Here’s a JavaScript function that estimates the value of an integral in the most primitive way possible.

I = (B, lo, hi, f) => {
    // Generate B random values uniformly between lo and hi.
    let xs = Array.from({length: B}, _ => lo + (hi - lo) * Math.random());
    // Compute the value of f at each location.
    let ys = xs.map(f);
    // Return the total area of each corresponding pile.
    return (hi-lo)*ys.reduce((r, y) => r + y, 0)/ys.length;
}

To compute an approximation to the value of the integral we’ve seen, we run

I(10_000,
  0, 1,
  x => (x-1)/Math.log(x)
);

0.6916867623261724

This is fairly close to 0.69. And we got there in four lines of JavaScript, as promised.

We can try this on the next example too. Now we’re asking about the integral

\[\int_0^{\frac{\pi}{2}} \frac{\ln{(1 - \sin{x})}}{\sin{x}} \mathrm{d}x\]

which, translated to JavaScript, becomes

I(10_000,
  0, Math.PI,
  x => Math.log(1 - Math.sin(x))/Math.sin(x)
);

-3.67

This is again fairly close to the desired −3.7, but not quite there yet. The tricky shape of the function is the reason we aren’t getting as close as we want.

At the upper endpoint of the integration interval, this function goes to negative infinity. The random piles we draw come primarily from the well behaved region of the function, and thus don’t help the computer realise this behaviour.

There are clever ways to sample adaptively from the trickier parts of the function, but an easy solution is to just visually find a breakpoint, split the interval on that, and then estimate the sensible part separately from the crazy-looking part. Since the total area must be the sum of both areas, we can add their results together for a final estimation.

In this case, we might want to pick e.g. 1.5 as the breakpoint, so we combine the area estimations from 0–1.5 and then 1.5–\(\frac{\pi}{2}\). The result is

I(2_000, 0, 1.5, x => Math.log(1 - Math.sin(x))/Math.sin(x))
+ I(8_000, 1.5, Math.PI/2, x => Math.log(1 - Math.sin(x))/Math.sin(x));

-3.70

which is indeed much closer to the actual value of −3.7.

Note that we aren’t taking more samples, we’re just sprinkling them more wisely over the number line. We spend 2,000 samples in the relatively well-behaved region where the function takes values from −1 to −6, and then we spend the other 8,000 samples in the small region that goes from −6 to negative infinity. Here it is graphically:

The reason this helps us is that this latter region contributes a lot to the value of the integral, but it is so small on the number line that we benefit from oversampling it compared to the other region. This is a form of sample unit engineering, which we have seen before in different contexts.

We can continue with some more examples from the Feynman trick article. That gets us the following table.

Notes:

1. The integration is from zero to infinity, but the function practically only has a value between zero and three, so that’s the region we estimate over.
2. This is another case where the function goes to infinity near zero, so we split up the estimation into one for the range 0–0.1, and the other for 0.1–1.0. We have not increased the sample count, only reallocated the 10,000 samples.
3. Again, the integration is from zero to infinity, but the function practically only has a value between zero and three, so that’s the region we estimate over.

“Now,” the clever reader says, “this is all well and good when we have the actual value to compare to so we know the size of the error. What will we do if we’re evluating a brand new integral? What is the size of the error then, huh?”

This is why we sampled the function randomly. That means our approximation is a statistical average over samples, and for that we can compute the standard error of the mean. In the JavaScript implementation below, we use the quick variance computation, but we could perhaps more intuitively have used the spc inspired method.

Ic = (B, lo, hi, f) => {
    let xs = Array.from(
      {length: B}, _ =>
      lo + (hi - lo) * Math.random()
    );
    let ys = xs.map(f);
    // Compute the variance of the ys from the sum and
    // the sum of squared ys.
    let s = ys.reduce((r, y) => r + y, 0);
    let ssq = ys.reduce((r, y) => r + y**2, 0);
    let v = (ssq - s**2/B)/(B-1);
    // Compute the mean and the standard error of the mean.
    let m = (hi-lo)*s/B;
    let se = (hi-lo)*Math.sqrt(v/B);
    // Compute the 90 % confidence interval of the value of
    // the integral.
    return {
        p05: m - 1.645*se,
        p95: m + 1.645*se,
    }
}

If we run this with the first integral as an example, we’ll learn that

Ic(10_000,
  0, 1,
  x => (x-1)/Math.log(x)
)

Object {
  p05: 0.6896
  p95: 0.6963
}

Not only is this range an illustration of the approximation error (small!), it is also very likely to capture the actual value of the integral. Here are some more examples from the same integrals as above:

These are all built naïvely from 10,000 uniform samples. In other words, in none of the cases have the computation been split up to allocate samples more cleverly.

Again, we could spend a lot of time learning to integrate by hand … or we ask the computer for less than a second of its time first, and see if the accuracy it can do it with is appropriate for our use case. In my experience, it generally is.

What’s neat is we can still split up the computation like we did before, if we believe it will make the error smaller and the confidence interval narrower. Let’s use the following integral as an example.

\[\int_0^\infty \frac{\sin{x}}{x} \mathrm{d}x\]

This oscillates up and down quite a bit for small \(x\), and then decays but still provides significant contributions for larger \(x\). A naive evaluation would have a confidence interval of

Ic(10_000, 0, 100, x => Math.sin(x)/x)

Object {
  p05: 1.461
  p95: 1.884
}

and while this is certainly correct2² The actual value of the integral is half \(\pi\) or approximatey 1.571., we can do better. We’ll estimate the region of 0–6 separately from 6–100, using half the samples for each3³ Why put the break point at 6? The period of sin is a full turn, which is roughly 6 radians. This ensures we get roughly symmetric contributions from both integrals. That’s not necessary for the technique to work, but it makes the illustration a little cleaner.:

Ic(5_000, 0, 6, x => Math.sin(x)/x)

Object {
  p05: 1.236
  p95: 1.468
}

This contains the bulk of the value of the integral, it seems. Let’s see what remains in the rest of it.

Ic(5_000, 6, 100, x => Math.sin(x)/x)

Object {
  p05: 0.080
  p95: 0.198
}

We can work backwards to what the standard errors must have been to produce these confidence intervals.4⁴ The midpoint is the point estimation for each region, and the standard error is 1/1.645 times the distance between the 5 % point and the midpoint.

The estimation of the total area would be the values summed, i.e. 1.5457. The estimation of the standard error of this we get through Pythagorean addition and it is approximately 0.05143. We convert it back to a confidence interval and compare with when we did not break it up into multiple components.

Although in this case the two methods happen to share a lower bound, the upper bound has been dramatically reduced. The total range of the confidence interval is more than halved! This was because we allocated the samples more cleverly – concentrated them in the early parts of the function – rather than increased the number of samples.

That said, we’re at a computer, so we could try increasing the sample count. Or maybe both?

It seems like sampling more cleverly has almost the same effect as taking ten times as many samples.

We could play around with where to put the breakpoint, and how many samples to allocate to each side of it, and see which combination yields the lowest error. Then we can run that combination with a lot of samples to get the most accurate final result. That would take maybe 15 minutes of tooting about and exploring sensible-seeming alternatives, so it’s probably still quicker than integrating by hand.

It should be said that there are times when numeric solutions aren’t great. I hear that in electronics and quantum dynamics, there are sometimes integrals whose value is not a number, but a function, and knowing that function is important in order to know how the thing it’s modeling behaves in interactions with other things.

Those are not my domains, though. And when that’s not the case, the computer beats Feynman any day of the week.