Monte Carlo Integration

Monte Carlo Integration is useful for calculating the area, volume, or n-dimensional content of a complicated space. With the facility of "good" random numbers, one can determine by experiment the probability of "landing" in the space. By multiplying the probability by the size of the encompassing domain, the size of the complicated domain is found.

My goal was to utilize Monte Carlo Integration to approximate the value of Pi. To do this, you must find a probability space that is dependent upon Pi, then perform and experiment to get the measured probability, then equate the expression for the probability with the measured probability and solve for Pi as if it were a variable. I found that the probability space of a pin falling on an infinite chessboard and not intersecting any square's border to be a probability space that is dependent on Pi and here is how.

Note: Although my decision was arbitrary in terms of effect, I chose to use a pin length equal to the length of a square for uniformity.

First I had to come up with a classification of the random variables. In this case there are a few different ways to look at the situation. One option is to take the landing point of one end of the pin and the angle of the stick as the random variables. This is helpful because the pin is asymmetrical with respect to the variables, so there won't be any problems with indistinguishable cases. However, I decided that the simplicity of taking the center of the pin as variable outweighed the asymmetrical factor.

Then I though about the position variable; where could the pin land and potentially not intersect? Here I simplified the situation to the case of one square, since all squares are the same on an infinite chessboard. If it was too close to the corner of this square it would be stuck, so the zero-probability regions would cluster there. I noticed that if the center of the pin was further than a half pinlength from a corner (but still in the same quadrant of the square) it would be able to put one end in the corner and the other end sticking out into the opposite quadrant. This does not work when closer than one half pinlength, so the quarter circles of radius 1/2 pinlength about each corner are the zero-probability regions. This means that the leftover regions must be counted- shown in black here:

Next I thought about the angle variable; for each point in this 4-pointed star, what angles would allow the pin to not intersect any border? After testing many points and visualizing the motion of the angles that do intersect, I noticed that these regions are always bordered by a line of half pinlength that connects the landing point to one of the square's borders. This is because if you rotate a pin about that location, it leaves the positive probability region when the end of the pin hits a border exactly. Now, the real trick is to notice that due to this uniform length, you can calculate the angles of the zero probability regions using a arccosine with the parameter of distance to the square's border. The distances from the borders are simply represented by the x and y coordinates of the landing position. There is one zero probability angle caused by each of the two closest borders - one for each of x and y. So, since arccos gives half of the angle cut off, and the opposite direction is also eliminated, the angle must be multiplied by 4 to get the full zero-probability angle sum. Therefore, the formula for the angle sum that is cut off by border intersection is (if pinlength is taken as 2): 4(arccos(x)+arccos(y)) To convert this into a probability, we must express it as a fraction of 2Pi, and the probability of not intersecting is 1 minus that value, so the probability of not intersecting at location (x,y) is: 1-(2/Pi)(arccos(x)+arccos(y))

With the probability expression, it is possible to calculate the probability volume. The previous expression gives the probability at each point. If the integral is of this expression over its two dimensional domain is taken, then the resulting volume when divided by the area represents the average probability. This averaging operation can be viewed as water taking the form of its container. Imagine the water is initially confined to the volume of the tent shape, and then let the water splash down to take the shape of a cube container. This is a very powerful natural averaging operation - the height of the top of the water is the average. The result of integration, the volume, is spread over a surface of 1 unit squared, so the height is equal to its volume (numerically). The integration can be done by integrating over the whole tent shaped region, but since it is symmetrical, it is possible, and much easier, to just integrate one quadrant. The process is slightly complicated by the fact that inside the zero probability region the expression still takes on values. During the generation of the expression these elements of the domain were not considered, and are thus garbage to this calculation. They can be eliminated by adding the integral of the absolute value and dividing the result by 2, or simply changing the limits of the integral to avoid this part. I chose the latter, and allowed my calculator to integrate: Int[Int[1-(2/Pi)(arccos(x)+arccos(y)),y,Sqrt[1-x^2],1],x,0,1] The output was (Pi-3)/Pi ; a pleasant result.

Finally it was time to measure the probability in the real world. I executed 1000 drops. Each digit below represents the count of non-intersections out of 10.

011011010010
110001021022
000001100100
000101110001
011101100101
101100100101
010011001010
000200321100
1000

There was a total of 52/1000 drops that did not intersect any borders -> Probability of 0.052

So then it was time to solve for Pi:

(Pi-3)/Pi = 0.052
1-(3/Pi) = 0.052
1-0.052 = 3/Pi
Pi = 3/0.948 ~= 3.16

So the result was not so accurate, but with a mere 1000 drops not much can be expected. Even if the probability was behaving perfectly (the number of non-intersecting was the precise number that has the highest probability of occurring) the approximation of Pi would only be accurate to 3 decimal places since any further accuracy would correspond to a non-integer number of non-intersections. Interestingly, if I count any set of ten with more than one non-intersect as only one, then I get that optimal number with 3 digit accuracy. My conclusion is that the Monte Carlo method can be very useful... as long as you don't have to perform it manually.

Note: I later found that this is called the Buffon-Laplace Needle Problem and a general formula has been found for the probability of intersection: P(l;a,b) = 2l(a+b)-l^2/(PIab) where l is the needle length, and a,b are the spacings between the horizontal and vertical lines (from Mathworld). With my specific case, this becomes P(1;1,1) = 2(1)(1+1)-1^2/(PI*1*1) = 3/PI, which is the same as my result. This implies something interesting. The original Buffon's Needle problem involved only horizontal lines, and the probability for that is 2/PI, so the probability for not intersecting one set AND not intersecting the other set would seem to be (1-2/PI)^2 = 1 - 4/PI + 4/PI^2 = (PI - 4)/PI + (4/PI)/PI = (PI-4+4/PI)/PI which is not equal to (PI-3)/PI. This implies that intersecting each set of lines is not independent. This is because the probabilities for the horizontal and vertical cases are related through the angle of the needle. Suppose we consider the location of the midpoint first and the angle second. After the location is set, there will be an arc of angles about 0 degrees for which the needle won't intersect the horizontal lines and an arc of angles about 90 degrees for which the needle won't intersect the vertical lines. The probability that the needle won't intersect either set of lines corresponds to the area of the intersection of the two arcs of angles. Since the probability is zero when the midpoint location is near the corners, the intersection is zero in these regions, so the probability in the case with both lines cannot be the product since both of the probabilities are nonzero, so they won't make a product of 0. Basically, when the needle is at an angle that makes it easier to not intersect one set of lines, it makes it more probable that it will intersect the other set of lines, so the probability of not intersecting both sets is much less.