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 ² from x = 0 to x = 1. Anyone who's had an introduction to calculus will remember that the indefinite integral of x ² is given by (1/3) x ³ + C, where C is a constant. It follows that the value of the definite integral is equal to [(1/3) 1 ² + C] - [(1/3) 0 ² + C] = 1/3. In other words, the magnitude of the area under the graph of y = x ² 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 ₁ = k P ₁ and A ₂ = k P ₂ then A ₁ / A ₂ = ( k P ₁ ) / ( k P ₂ ) or A ₁ / A ₂ = P ₁ / P ₂. 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 ². 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) ² + (y - 0.5) ² <= 0.5 ². 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 ³. If V_s and V_c represent the volume of the sphere and the cube, respectively, then V_s / V_c = pi / 6. Our test for points that lie within the ball is (x - 0.5) ² + (y - 0.5) ² + (z - 0.5) ² <= 0.5 ². 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 V_s represent the "volume" of an n-dimensional hypersphere of diameter D.

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

If n is even, V_s / V_c = pi ^m / ( m! 2 ⁿ) where m = n / 2.

If n is odd, V_s / V_c = ( m! pi ^m ) / n! where m = (n - 1)/ 2.

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

x₁² + x₂² + ... + x_n² = R ².

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?

Addendum

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 ₁) + f(x ₂) + ... + 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!

Addendum ²

Okay, an Addendum to the Addendum may be a bit much, but I just remembered something. The mathematical mean of a series of values can be extremely misleading if a few of those values are unusually high or low in comparison to the rest of the series. In response to this problem, somebody came up with the idea of the sample median.

To find the sample median, first sort the series of values in question. If n is odd, the central element of the series is our sample median. If n is even, our sample median is the mathematical mean of the two values at the center of the series. I would very much like to see someone incorporate the use of the sample median into a hit-or-miss Monte Carlo program. As usual, it would make a smashing article for a future issue of the Basix Fanzine!

C'ya,
RudeJohn
"I'm rude. It's a job."

PS: If it looks as though I've been a bit heavy-handed in my constant solicitation of articles for the 'zine all I can say is "Glad you noticed!" Now why don't you write something up and send it in? Got an idea, an opinion, a story, or even a good joke? Send it in! Heck, think of the 'zine as a multi-newsgroup bulletin board. Advertise your website, your proggies, and even your latest game! Please, talk to us. We're so lonely. <sniff>