Hit-or-Miss Integration

Warning: All code was play-tested in the QB 4.5 IDE.

Introduction

This article discusses the simplest type of probabilistic numerical integration which is sometimes referred to as the hit-or-miss Monte Carlo method. This approach may not compete with deterministic methods in general, but it's something to think about when you wake up in a cold sweat at 4 a.m. I guess. <shrug>

Parabolas

Consider the definite integral of x 2 from x = 0 to x = 1. Anyone who's had an introduction to calculus will remember that the indefinite integral of x 2 is given by (1/3) x 3 + C, where C is a constant. It follows that the value of the definite integral is equal to [(1/3) 1 2 + C] - [(1/3) 0 2 + C] = 1/3. In other words, the magnitude of the area under the graph of y = x 2 from (0,0) to (1,1) is equal to 1/3.

That said, let's find the value of the definite integral in a very different fashion as outlined by the following algorithm. For the sake of discussion, random is an ideal random-number generator over [0,1] while num represents the number of randomly-generated points.

``` procedure estimate_area area = 0 loop from 1 to num {num > 0} x = random y = random if y <= x 2 then (area = area + 1) end loop return (area / num) end procedure ```

If all goes well, this algorithm should return an estimate to within 99% accuracy or better. This is not significant in terms of absolute error, but what do you expect? Fortunately, the rationale is quite simple. (If it wasn't, I'd have to find someone who could explain it to me.) For a sufficiently large number of randomly-generated points, the number of points in one unit-area should be nearly equal to the number of points in any other unit-area. This means that any area A n is proportional to P n the number of points it contains. So, if A 1 = k P 1 and A 2 = k P 2 then A 1 / A 2 = ( k P 1 ) / ( k P 2 ) or A 1 / A 2 = P 1 / P 2. See 11-x2.bas

The reader is invited to generalize the above algorithm for y = f(x) over the interval [a,b]. Is this practical? What do you need to know about f(x) over [a,b] before applying our hit-or-miss method?

Pi a la Rude

Consider a circle of diameter D inscribed within a square of area D 2. The ratio of the area of the circle to the area of the square is pi / 4. As with the previous example, we will use randomly-generated points to find an approximate value for pi. If D = 1 and the center-point of the circle is (0.5,0.5), then our test for points which lie on the disk is (x - 0.5) 2 + (y - 0.5) 2 <= 0.5 2. See 11-pi2d.bas.

Funny thing is, if we restrict our randomly-generated points to the line-segment from (0,0) to (1,1), the ratio of the total number of points on the line-segment to the number of points on the circle's diameter approaches 2 1/2. Why?

By now, you must be wondering if this experiment can be extended to three dimensions. Well, of course it can! <grin> This time, we'll need a sphere of diameter D within a cube of volume D 3. If Vs and Vc represent the volume of the sphere and the cube, respectively, then Vs / Vc = pi / 6. Our test for points that lie within the ball is (x - 0.5) 2 + (y - 0.5) 2 + (z - 0.5) 2 <= 0.5 2. See 11-pi3d.bas.

This time, if we restrict our randomly-generated points to the line-segment from (0,0,0) to (1,1,1), the ratio of the total number of points on the line-segment to the number of points on the sphere's diameter approaches 3 1/2. Why?

Hyperspace

The next logical(?) step is to generalize our experiment using n-dimensional hyperspatial objects, i.e., the hypersphere and hypercube.

Let Vs represent the "volume" of an n-dimensional hypersphere of diameter D.

Let Vc = D n represent the "volume" of an n-dimensional hypercube.

If n is even, Vs / Vc = pi m / ( m! 2 n) where m = n / 2.
If n is odd, Vs / Vc = ( m! pi m ) / n! where m = (n - 1)/ 2.

The reader is invited to write a program that will use our probabilistic calculation of Vs / Vc to estimate pi in the general case. Keep in mind that the formula for an n-dimensional hypersphere of radius R is given by

x12 + x22 + ... + xn2 = R 2.
Everyone who submits an entry to the Basix Fanzine will receive a lifetime subscription absolutely free! First prize is a lovely string of multi-colored paper-clips. Really.

Extra Credit

If we restrict randomly-generated points in the hypercube to the line-segment from (0,0,...,0) to (1,1,...,1), what can we say about the ratio of the number of points on the line-segment to the number of points lying on the hypersphere's diameter?

Deterministic algorithms can sometimes be hybridized to include probabilistic methodology. For example, the so-called "trapezoidal rule" may be improved(?) as follows. Rather than using the mean of f(x) taken at the endpoints of an interval [a,b], or [f(b) - f(a)]/2, we take the mean of f(x) for n randomly-generated values of x over the interval [a,b], i.e., [ f(x 1) + f(x 2) + ... + f(x n)] / n for a <= x i <= b. Once again, the reader is encouraged to submit either an algorithm or example program employing this technique to the Basix Fanzine for future publication. Let's face it, if you had anything better to do then you wouldn't be reading this!