Porter Johnson, IIT
DOUBLE BONUS
Tau Beta Pi Bent: Fall 1998
I hope that you do not misconstrue this as bragging, but my 3-year old twins are destined for greatness. Why just this morning they were investigating geometric probabilities. First, the older one drove three nails in random positions into the top of our rectangular [square] dining-room table. Then the younger one formed a triangle by wrapping the string, and then he drove a fourth random nail into the table top. What is the exact probability that the fourth nail was within the triangle? (My wife seems curiously unimpressed with their accomplishments.) -B R Adams, TX A '58Clarification: The word "random" is ill-defined, unless you specify the "probability measure" from which events are selected. Your chances of arriving to work after driving on the Dan Ryan Expressway are random [but pretty good], whereas your chances of winning the jackpot lottery are random [and lousy]. I interpret this problem to mean that the probability of choosing any point within some region of a unit square is equal to the area of that region.
Choose three points inside unit square
Pick three points at random [with probability measure induced by cartesian
metric!]:
Put away those geometry books; you're physicists!
= | y1(x2-x3) + y2(x3-x1) + y3(x1-x2) | /2
So, ... , we need only to evaluate the following six-dimensional integral
[gulp!]:
Isn't there an EASIER way to do this?
We know that it we choose two points at random inside the square, and draw a line passing through them, the third point will lie on either side of the line with equal probability. Therefore, if we consider all three lines on the triangle, isn't the probability of ending up inside just 1/8?
Isn't That Correct?
BACK TO THE SALT MINES!
Mathematica to the rescue?
{x3,0,1},{y1,0,1},{y2,0,1},{y3,0,1}]
Mathematica's answer: 0.070910993. [wrrrrrrrrong!]
Try it again: N[P,11] . Answer [unreliable].
The problem: mathematica's integration packages like smooth functions, which this is not. Mathematica lets you down if naively you expect too much of it. There is no substitute for understanding what is happening, as was noted a long time ago.
Another futile attempt to avoid thinking about this problem---Monte Carlo simulation.
It is convenient to use vectors to decide whether a point is inside a triangle. Let us choose one vertex [#1] as being special, and call the vectors from it to the other vertices as R2 and R3, respectively. Let the vector to the point in question be R4 . The criterion for being inside the triangle is R4 = a R2 + b R3 , where a and b are both positive, with a + b < 1.
Results: [FORTRAN program listing on next page]
| Number of Shots | Number of Hits | Probability | ± | Error |
|---|---|---|---|---|
| 5,000,000 | 381,578 | .0763156 | ± | .00045 |
| 10,000,000 | 764,127 | .0765347 | ± | .00031 |
| 50,000,000 | 3,820,381 | .07640762 | ± | .00014 |
| 100,000,000 | 7,637,923 | .07637924 | ± | .00010 |
| 200,000,000 | 15,276,160 | .07638080 | ± | .00007 |
A Challenge to MATHEMATICA aficionados: can you obtain this result from Mathematica?? If not, what good is it?
Actually, this is an unsubtle attempt on my part to learn to coax Mathematica to get useful, difficult results.
Program Square
c 4 points chosen at random with uniform weight inside unit square
c determine probability that fourth point lies inside triangle formed by other 3
c monte-carlo simulation
write (6,*) ' put in number of shots'
read (5,*) nshots
c seed program to set up ms fortran random number generator
nrand = 11
call seed(nrand)
c nhits is counter for number of 'hits' [4th point inside triangle]
nhits = 0
do 20 j = 1, nshots
c random number generator to get the four points: #1 (x1, y1) ...
call random(x1)
call random(y1)
call random(x2)
call random(y2)
call random(x3)
call random(y3)
call random(x4)
call random(y4)
c distances from point #1: #2: :r2 =(dx2, dy2) ...
dx2 = x2 - x1
dx3 = x3 - x1
dx4 = x4 - x1
dy2 = y2 - y1
dy3 = y3 - y1
dy4= y4 - y1
c scalar products of vectors
r22 = dx2**2 + dy2**2
r23 = dx2*dx3 + dy2*dy3
r33 = dx3**2 + dy3**2
r24 = dx2*dx4 + dy2*dy4
r34 = dx3*dx4 + dy3*dy4
c calculation of /alpha and /beta in the formula
c /vec{r4} = /alpha /vec{r2} + /beta /vec{r3}
det = r22*r33 - r23**2
c det = 0 if the three points are collinear, and we stop
if (det.eq.0) goto 20
alpha = (r24*r33 - r34*r23)/det
beta = (r22*r34 - r23*r24)/det
c tests for lying inside triangle: alpha>0 beta>0 alpha+beta<1
if (alpha.gt.0.) then
if (beta.gt.0) then
if ((alpha+beta) .lt. 1.0) nhits = nhits+1
end if
end if
20 continue
write(6,99) nshots, nhits
99 format (' shots: ',i10, ' hits:', i10)
stop
end