Chicago Section AAPT Meeting 06 November 1999

Monte-Carlo Games on a Square Table

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 '58

Clarification: 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!]:

#1: (x1,y1)
#2: (x2,y2)
#3: (x3,y3)

When one picks a fourth point, (x4,y4), the probability that it lies inside the triangle shown is equal to the area of that triangle. How do we calculate the area of the triangle?

Put away those geometry books; you're physicists!

Area(x1, x2, x3, y1, y2, y3) = 1/2 | (r2 - r1) ( r3 - r1) |

= | y1(x2-x3) + y2(x3-x1) + y3(x1-x2) | /2

So, ... , we need only to evaluate the following six-dimensional integral [gulp!]:

ò 01 dx1 ò 01 dx2 ò 01 dx3 ò 01 dy1 ò 01 dy2 ò 01 dy3 Area(x1, ..., y3)

WAIT A MINUTE!

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?

The answer: "nicht einfach; aber falsch" -Wolfgang Pauli
[it's not simple, but it is wrong].

BACK TO THE SALT MINES!

Mathematica to the rescue?

P=0.5 Integrate[Abs[(x2-x1)(y3-y1) - (x3-x1)(y2-y1)], {x1,0,1},{x2,0,1},

{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.

SUBSTITUTES
ARE LIKE A GIRDLE
THEY FIND SOME JOBS
THEY JUST CAN'T HURDLE
-BURMASHAVE

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,000381,578 .0763156± .00045
10,000,000764,127 .0765347± .00031
50,000,0003,820,381 .07640762± .00014
100,000,0007,637,923 .07637924± .00010
200,000,00015,276,160 .07638080± .00007
Notice that the answers are close to, but never identical with, the fraction

11/144 = .0763888 ...

CAN WE PROVE THAT 11/144 IS/ISN'T CORRECT? NOT THIS WAY!

Another Attempt: worse than brute force:
ò 01 dy1 ò 01 dy2 ò 01 dy3 Area(x1, x2, x3, y1, y2, y3) =
[2(x1-x2)2 + 3(x1-x2)(x2-x3) + 2(x2-x3)2] / (12 (x1-x3))

subject only to the constraint: 1> x1 > x2 > x3 > 0.

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.

The result P = 11/144 follows:

      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