Throwing Buffon’s Needle with Mathematica

Enis Siniksaran

It has long been known that Buffon’s needle experiments can be used to estimate π. Three main factors influence these experiments: grid shape, grid density, and needle length. In statistical literature, several experiments depending on these factors have been designed to increase the efficiency of the estimators of π and to use all the information as fully as possible. We wrote the package BuffonNeedle to carry out the most common forms of Buffon’s needle experiments. In this article we review statistical aspects of the experiments, introduce the package BuffonNeedle, discuss the crossing probabilities and asymptotic variances of the estimators, and describe how to calculate them using Mathematica.

Introduction

Buffon’s needle problem is one of the oldest problems in the theory of geometric probability. It was first introduced and solved by Buffon [] in 1777. As is well known, it involves dropping a needle of length l at random on a plane grid of parallel lines of width d>l units apart and determining the probability of the needle crossing one of the lines. The desired probability is directly related to the value of π, which can then be estimated by Monte Carlo experiments. This point is one of the major aspects of its appeal. When π is treated as an unknown parameter, Buffon’s needle experiments can be seen as valuable tools in applying the concepts of statistical estimation theory, such as efficiency, completeness, and sufficiency. For instance, in order to obtain better estimators of π, Kendall and Moran [] and Diaconis [] examine several aspects of the problem with a long needle (l>d). Morton [] and Solomon [] provide the general extension of the problem. Perlman and Wishura [] investigate a number of statistical estimation procedures for π for the single, double, and triple grids. In their study, they show that moving from single to double to triple grid, the asymptotic variances of the estimators get smaller and hence more efficient estimators can be obtained. Wood and Robertson [] introduce the concept of grid density and provide an alternative idea. They show that Buffon’s original single grid is actually the most efficient if the needle length is held constant (at the distance between lines on the single grid) and the grids are chosen to have equal grid density (i.e., equal length of grid material per unit area). In [], Wood and Robertson investigate the ways of maximizing the information in Buffon’s experiments.

We organize this article as follows. In the first three sections, we review Buffon’s experiments on single, double, and triple grids and their statistical issues. In the next section, we introduce the features of the package BuffonNeedle. The functions in the package implement Monte Carlo experiments for the three types of grids. The results of each experiment are given in a table and in a picture. When the number of the needles thrown on each grid is large, very nice pictures exhibit the interface between chance and necessity. In the last two sections, we describe how to calculate the crossing probabilities in single- and double-grid experiments and the asymptotic variances of the estimators for each grid using Mathematica.

Single-Grid Experiment

The single-grid form is Buffon’s well-known original experiment. A plane (table or floor) has parallel lines on it at equal distances d from each other. A needle of length l (l<d) is thrown at random on the plane. Figure  shows a single grid with two needles of length l representing two possible outcomes. The probabilities of crossing zero lines and one line are

({{(Subscript[p, 0]), =, 1-2r θ,}, {(Subscript[p, 1]), =, 2r θ,}})

where r=l/d and θ=1/π. These results can be found in many probability and statistics textbooks, for instance [, , ]. They are also explained in the “Calculating the Crossing Probabilities in Single and Double Grids” section of this article.

{abysis1,abysis2,ordinates1,ordinates2}={{0.927`,3.102`},{1.639`,2.57`},{2.623`,2.237`},{2.159`,1.574`}};
{points1,points2}=Transpose[Mean/@{{abysis1,abysis2},{ordinates1,ordinates2}}];
fig1=Show[Graphics[{Table[Line[{{abysis1[[i]],ordinates1[[i]]},{abysis2[[i]],ordinates2[[i]]}}],{i,1,2}],{RGBColor[0,0,0],Thickness[0.001`],{Arrow[{{0.25`,1},{0.25`,2}}]},{Arrow[{{0.25`,2},{0.25`,1}}]},{Arrow[{{0.25`,2},{0.25`,3}}]},{Arrow[{{0.25`,3},{0.25`,2}}]}},Text[d,{0.05`,2.5`},{-1,0}],Text[d,{0.05`,1.5`},{-1,0}],{PointSize[(1/55)],RGBColor[0,1,0],Point[points1]},{PointSize[(1/55)],RGBColor[1,0,0],Point[points2]},{RGBColor[0.15`,0.35`,0.25`],Thickness[0.006`],Line[{{4.5`,1},{0.5`,1}}],Line[{{0.5`,1},{4.5`,1}}],Line[{{0.5`,2},{4.5`,2}}],Line[{{0.5`,3},{4.5`,3}}]}}],PlotRange->All,AxesOrigin->{0,0},AspectRatio->0.55`,Axes->None]

From equation (), we can write

θ=(1/π) =((Subscript[p, 1])/(2r)).

Let (Subscript[N, 1]) be the number of times in n independent tosses that the needle crosses any line. Then the proportion of crossings (Subscript[(OverHat[p]), 1]), a point estimator of (Subscript[p, 1]), becomes (Subscript[(OverHat[p]), 1])=(Subscript[N, 1])/n. Hence, we can write the point estimator of θ in equation () as

(OverHat[θ ])=((Subscript[(OverHat[p]), 1])/(2r))=((Subscript[N, 1])/(2 r n)).

The random variable (Subscript[N, 1]) is binomially distributed with the parameters n and (Subscript[p, 1]). (OverHat[θ ]) is the uniformly minimum variance unbiased estimator (UMVUE). Furthermore, it is the maximum likelihood estimator (MLE) of θ and hence has 100% asymptotic efficiency in this experiment []. The variance of (OverHat[θ ]) is then

Var((OverHat[θ ]))=Var(((Subscript[N, 1])/(2r n)))=(((Subscript[p, 1])(1-(Subscript[p, 1])))/(4(r^2)n))=((θ (1-2 rθ))/(2 n r))=(θ/(2n))((1/r)-2θ),

which is minimized by taking r as close as possible to 1. Choosing needle length l=d (r=1) ensures this purpose. In this case, the variance of (OverHat[θ ]) becomes

Var((OverHat[θ ]))=(((θ )^2)/(2n))((1/θ)-2).

In Buffon’s experiments, the parameter of main interest is π rather than θ. The estimator of this parameter, (OverHat[π])=1/(OverHat[θ]), is called Buffon’s estimator and can be obtained from equation () as

(OverHat[π])=((2r)/(Subscript[(OverHat[p]), 1])).

It can also be expressed in terms of (Subscript[(OverHat[p]), 0])=(Subscript[N, 0])/n

(OverHat[π])=((2r)/(1-(Subscript[(OverHat[p]), 0]))),

where (Subscript[N, 0]) is the number of times in n tosses that the needle does not cross any line. Using equation () or () and Monte Carlo methods, we can obtain empirical estimates of π for various values of r. The best estimate is expected at r=1 (l=d). Standard theory [] ensures that Buffon’s estimator is an asymptotically unbiased, 100% efficient estimator of 1/θ. Applying the delta method shows that its asymptotic variance is

AVar((OverHat[π]))=((π^( 2))/(2n))((π/r)-2),

which is, as expected, minimized at r=1. For this value of r, the asymptotic variance of (OverHat[π]) is

AVar((OverHat[π]))=(((π )^2)/(2n))(π-2).

If it is evaluated at π=3.1416, we have

AVar((OverHat[π]))=(5.63/n).

Double-Grid Experiment

In the double-grid experiment, also called the Laplace extension of Buffon’s problem, a plane is covered with two sets of parallel lines where one set is orthogonal to the other.

{abysis1,abysis2,abysis3,ordinates1,ordinates2,ordinates3}={{1.864`,1.215`},{2.701`,3.439`},{1.907`,2.4`},{2.943`,2.568`},{1.275`,1.143`},{1.373`,0.807`}};
{points1,points2,points3}={{Mean[abysis1],Mean[ordinates1]},{Mean[abysis2],Mean[ordinates2]},{Mean[abysis3],Mean[ordinates3]}};
fig2=Show[Graphics[{Line[Table[{abysis1[[i]],ordinates1[[i]]},{i,1,2}]],Line[Table[{abysis2[[i]],ordinates2[[i]]},{i,1,2}]],Line[Table[{abysis3[[i]],ordinates3[[i]]},{i,1,2}]],{PointSize[(1/55)],RGBColor[0,1,0],Point[points1]},{PointSize[(1/55)],RGBColor[0,0,1],Point[points2]},{PointSize[(1/55)],RGBColor[1,0,0],Point[points3]},Text[d,{0.05`,2.5`},{-1,0}],Text[d,{0.05`,1.5`},{-1,0}],Text[d,{1.5`,0.05`},{-1,0}],Text[d,{2.5`,0.05`},{-1,0}],Text[d,{3.5`,0.05`},{-1,0}],{Arrow[{{0.25`,1},{0.25`,2}}]},{Arrow[{{0.25`,2},{0.25`,1}}]},{Arrow[{{0.25`,2},{0.25`,3}}]},{Arrow[{{0.25`,3},{0.25`,2}}]},{Arrow[{{1,0.25`},{2,0.25`}}]},{Arrow[{{2,0.25`},{1,0.25`}}]},{Arrow[{{3,0.25`},{2,0.25`}}]},{Arrow[{{2,0.25`},{3,0.25`}}]},{Arrow[{{3,0.25`},{4,0.25`}}]},{Arrow[{{4,0.25`},{3,0.25`}}]},{RGBColor[0.15`,0.35`,0.25`],Thickness[0.006`],Line[{{1,1},{4,1}}],Line[{{1,2},{4,2}}],Line[{{1,3},{4,3}}],Line[{{1,1},{1,3}}],Line[{{2,1},{2,3}}],Line[{{3,1},{3,3}}],Line[{{4,1},{4,3}}]},{RGBColor[0.15`,0.35`,0.25`],Thickness[0.006`],Dashing[{0.01`}],Line[{{0.5`,1},{1,1}}],Line[{{0.5`,2},{1,2}}],Line[{{0.5`,3},{1,3}}],Line[{{1,3},{1,3.5`}}],Line[{{1,1},{1,0.5`}}],Line[{{2,3},{2,3.5`}}],Line[{{2,1},{2,0.5`}}],Line[{{3,1},{3,0.5`}}],Line[{{4,3},{4.5`,3}}],Line[{{4,2},{4.5`,2}}],Line[{{4,1},{4.5`,1}}],Line[{{4,3},{4,3.5`}}],Line[{{4,1},{4,0.5`}}],Line[{{3,3},{3,3.5`}}]}}],Axes->None]

In Figure , we see a double-grid plane and three needles of length l crossing zero, one, and two lines. These crossing probabilities are

({{(Subscript[p,  0]), =, 1-r(4-r)θ,}, {(Subscript[p,  1]), =, 2r(2-r)θ,}, {(Subscript[p,  2]), =, (r^( 2))θ.}})

Let (Subscript[N, i]) be the number of times in n tosses that the needle crosses exactly i lines (i=0,1,2). Perlman and Wichura [] showed that the random variable (Subscript[N, 1])+(Subscript[N, 2]) is distributed binomially with the parameters n and m θ and is completely sufficient for θ where m is 4 r-(r^2). They also showed that the random variable

(OverHat[θ ])=(((Subscript[N, 1])+(Subscript[N, 2]))/(m n))

is the UMVUE and has 100% asymptotic efficiency with the variance

Var((OverHat[θ ]))=(θ/n)((l/m)-θ).

As in the case of a single grid, the variance of (OverHat[θ ]) is minimized by r=1, or equivalently by l=d. Replacing (Subscript[N, 1])+(Subscript[N, 2]) with n-(Subscript[N, 0]) in the right-hand side of equation (), we have

(OverHat[θ ])=((n-(Subscript[N, 0]))/(m n))=(1/m)(1-((Subscript[N, 0])/n))=((1-(Subscript[(OverHat[p]), 0]))/(4r-(r^( 2)))).

Then Buffon’s estimator, (OverHat[π])=1/(OverHat[θ ]), can be expressed as

(OverHat[π])=((4r-((r )^2))/(1-(Subscript[(OverHat[p ]), 0]))),

which can be used to obtain empirical estimates of π. By the delta method, we can obtain the asymptotic variance of (OverHat[π]) as

AVar((OverHat[π]))=(( (π^( 2))(4r-(r^( 2))-π))/(n r(r-4))),

which is minimized at r=1. For this value of r, it becomes

AVar((OverHat[π]))=(( (π^( 2)))/(3 n))(π-3).

When evaluated at π=3.1416, it is

AVar((OverHat[π]))=(0.466/n).

Compare the last equation with equation (). Buffon’s estimator in the double-grid experiment is 5.63/0.466≃12 times as efficient as that in the single-grid experiment.

Triple-Grid Experiment

In the triple-grid experiment, a plane is covered with equilateral triangles of altitude d and hence of side 2d/(Sqrt[3]).

layn[{x1_,y1_},{x2_,y2_}]:=y1+(((y2-y1) (x-x1))/(x2-x1));
{point1,point2,point3,point4,point5,point6,point7}={{2.5`+(1/(Sqrt[3])),1},{2.5`-(1/(Sqrt[3])),1},{2.5`-(2/(Sqrt[3])),2},{2.5`-(1/(Sqrt[3])),3},{2.5`+(1/(Sqrt[3])),3},{2.5`+(2/(Sqrt[3])),2},{2.5`,2}};
points={point1,point2,point3,point4,point5,point6,point7};
frame=Table[{points[[i-1]],Insert[points,point1,7][[i]]},{i,2,7}];
mainframe=Table[{point7,points[[i]]},{i,1,6}]⋃Table[{points[[i-1]],Insert[points,point1,7][[i]]},{i,2,7}];
extend={{{2.5`-(2/(Sqrt[3])),2},{2.5`-((2.5`)/(Sqrt[3])),layn[point3,point4]/.x->2.5`-((2.5`)/(Sqrt[3]))}},{{2.5`-(2/(Sqrt[3])),2},{2.5`-((2.5`)/(Sqrt[3])),layn[point2,point3]/.x->2.5`-((2.5`)/(Sqrt[3]))}},{{2.5`-(2/(Sqrt[3])),2},{2.5`-(3/(Sqrt[3])),2}},{{2.5`-(1/(Sqrt[3])),3},{2.5`-((0.5`)/(Sqrt[3])),layn[point3,point4]/.x->2.5`-((0.5`)/(Sqrt[3]))}},{{2.5`-(1/(Sqrt[3])),3},{2.5`-((1.5`)/(Sqrt[3])),layn[point7,point4]/.x->2.5`-((1.5`)/(Sqrt[3]))}},{{2.5`-(1/(Sqrt[3])),3},{2.5`-(2/(Sqrt[3])),3}},{{2.5`+(1/(Sqrt[3])),3},{2.5`+((0.5`)/(Sqrt[3])),layn[point6,point5]/.x->2.5`+((0.5`)/(Sqrt[3]))}},{{2.5`+(1/(Sqrt[3])),3},{2.5`+(2/(Sqrt[3])),3}},{{2.5`+(1/(Sqrt[3])),3},{2.5`+((1.5`)/(Sqrt[3])),layn[point7,point5]/.x->2.5`+((1.5`)/(Sqrt[3]))}},{{2.5`+(2/(Sqrt[3])),2},{2.5`+((2.5`)/(Sqrt[3])),layn[point6,point1]/.x->2.5`+((2.5`)/(Sqrt[3]))}},{{2.5`+(2/(Sqrt[3])),2},{2.5`+(3/(Sqrt[3])),2}},{{2.5`+(2/(Sqrt[3])),2},{2.5`+((2.5`)/(Sqrt[3])),layn[point6,point5]/.x->2.5`+((2.5`)/(Sqrt[3]))}},{{2.5`+(1/(Sqrt[3])),1},{2.5`+((1.5`)/(Sqrt[3])),layn[point7,point1]/.x->2.5`+((1.5`)/(Sqrt[3]))}},{{2.5`+(1/(Sqrt[3])),1},{2.5`+(2/(Sqrt[3])),1}},{{2.5`+(1/(Sqrt[3])),1},{2.5`+((0.5`)/(Sqrt[3])),layn[point6,point1]/.x->2.5`+((0.5`)/(Sqrt[3]))}},{{2.5`-(1/(Sqrt[3])),1},{2.5`-((0.5`)/(Sqrt[3])),layn[point3,point2]/.x->2.5`-((0.5`)/(Sqrt[3]))}},{{2.5`-(1/(Sqrt[3])),1},{2.5`-((1.5`)/(Sqrt[3])),layn[point7,point2]/.x->2.5`-((1.5`)/(Sqrt[3]))}},{{2.5`-(1/(Sqrt[3])),1},{2.5`-(2/(Sqrt[3])),1}}};
lines=mainframe⋃extend;
{abysis1,abysis2,abysis3,abysis4,ordinates1,ordinates2,ordinates3,ordinates4}={{1.629`,2.329`},{2.418`,2.88`},{2.755`,2.126`},{2.854`,3.456`},{1.848`,1.845`},{2.839`,2.313`},{1.795`,1.489`},{3.03`,2.673`}};
{points1,points2,points3,points4}={{Mean[abysis1],Mean[ordinates1]},{Mean[abysis2],Mean[ordinates2]},{Mean[abysis3],Mean[ordinates3]},{Mean[abysis4],Mean[ordinates4]}};
framepic=Show[Graphics[{{RGBColor[0.15`,0.35`,0.25`],Thickness[0.006`],Table[Line[mainframe[[i]]],{i,1,Length[mainframe]}]},{RGBColor[0.15`,0.35`,0.25`],Thickness[0.006`],Dashing[{0.01`}],Table[Line[extend[[i]]],{i,1,Length[extend]}]}}],AspectRatio->1,DisplayFunction->Identity];
fig3=Show[Graphics[{Line[Table[{abysis1[[i]],ordinates1[[i]]},{i,1,2}]],Line[Table[{abysis2[[i]],ordinates2[[i]]},{i,1,2}]],Line[Table[{abysis3[[i]],ordinates3[[i]]},{i,1,2}]],Line[Table[{abysis4[[i]],ordinates4[[i]]},{i,1,2}]],{PointSize[(1/55)],RGBColor[0,1,0],Point[points1]},{PointSize[(1/55)],RGBColor[0,0,1],Point[points2]},{PointSize[(1/55)],RGBColor[1,1,0],Point[points3]},{PointSize[(1/55)],RGBColor[1,0,0],Point[points4]},{RGBColor[0,0,0],Thickness[0.001`],{Arrow[{{0.25`,1},{0.25`,2}}]},{Arrow[{{0.25`,2},{0.25`,1}}]},{Arrow[{{0.25`,2},{0.25`,3}}]},{Arrow[{{0.25`,3},{0.25`,2}}]},{Arrow[{{2.5`+(1/(Sqrt[3])),0.25`},{2.5`-(1/(Sqrt[3])),0.25`}}]},{Arrow[{{2.5`-(1/(Sqrt[3])),0.25`},{2.5`+(1/(Sqrt[3])),0.25`}}]}},Text[d,{0.05`,2.5`},{-1,0}],Text[d,{0.05`,1.5`},{-1,0}],Text["\!\(\*FractionBox[\(2  d\), SqrtBox[\(3\)]]\)",{2.5`,0.05`},{0,1}],{RGBColor[0.15`,0.35`,0.25`],Thickness[0.006`],Table[Line[mainframe[[i]]],{i,1,Length[mainframe]}]},{RGBColor[0.15`,0.35`,0.25`],Thickness[0.006`],Dashing[{0.01`}],Table[Line[extend[[i]]],{i,1,Length[extend]}]}}],AspectRatio->0.75`,Axes->None]

Figure  shows a triple-grid plane and four needles of length l crossing zero, one, two, and three lines. In [], the crossing probabilities are given as

({{(Subscript[p,  0]), =, 1+((r^( 2))/2)-(3/2)r(4-((Sqrt[3])/2)r)θ,}, {(Subscript[p,  1]), =, -(5/4)(r^( 2))+(3/2)r(4-((Sqrt[3])/2)r)θ,}, {(Subscript[p,  2]), =, (r^( 2))-((3(Sqrt[3]))/4)((r )^2)θ,}, {(Subscript[p,  3]), =, -(((r )^2)/4)+((3(Sqrt[3]))/4)((r )^2)θ.}})

Let (Subscript[N, i]) denote the number of times in n tosses that the needle crosses exactly i lines (i=0,1,2,3). For this experiment, Perlman and Wichura [] investigated the random variable (Subscript[N, 1])+(Subscript[N, 2])+(Subscript[N, 3]) which is distributed binomially with the parameters n and α θ-1/2 where α=3r/2(4-(Sqrt[3])r/2). They proposed, among others, the following unbiased estimator of θ as a function of (Subscript[N, 1])+(Subscript[N, 2])+(Subscript[N, 3])

(OverHat[θ ])=(1/α)((((Subscript[N, 1])+(Subscript[N, 2])+(Subscript[N, 3]))/n)+(1/2)).

By replacing (Subscript[N, 1])+(Subscript[N, 2])+(Subscript[N, 3]) with n-(Subscript[N, 0]), as in the other experiments, we obtain the same estimator as a function of (Subscript[N, 0])

(OverHat[θ ])=(1/α)(((n-(Subscript[N, 0]))/n)+(1/2))=(1/α)((3/2)-(Subscript[(OverHat[p]), 0])).

The variance of (OverHat[θ ]) is

Var((OverHat[θ ]))=(((2αθ-3) (2 αθ-1))/(4 (α^2)n)).

As in the cases of the single- and double-grid experiments, the variance of (OverHat[θ ]) is minimized by taking r=1 (l=d).

From equation (), Buffon’s estimator can be written as

(OverHat[π])=((2α)/(3-2(Subscript[(OverHat[p]), 0])))=((3r(8-(Sqrt[3])r))/(2(3-2(Subscript[(OverHat[p]), 0])))).

For this experiment, the asymptotic variance of Buffon’s estimator is

AVar((OverHat[π]))=-(2(π^2) (-3 (Sqrt[3])+π) (-3 (Sqrt[3])+4 π) (-24+3 (Sqrt[3]) r+5 π r) (3 r (-8+(Sqrt[3]) r)+2 π (2+(r^2))))/(9 n r (54 (64-16 (Sqrt[3]) r+3 (r^2))+(π^2) (512-128 (Sqrt[3]) r-270 (r^2)+96 (Sqrt[3]) (r^3)+27 (r^4))+6 π (-320 (Sqrt[3])+168 r+234 (Sqrt[3]) (r^2)-306 (r^3)+27 (Sqrt[3]) (r^4)))).

For r=1 and π=3.1416, it is

AVar((OverHat[π]))=(0.015781/n).

Comparing this with equations () and (), we can infer that Buffon’s estimator in the triple-grid experiment is 0.466/0.015781≃29 times as efficient as in the double-grid experiment and 5.63/0.015781≃356 times as efficient as in the single-grid experiment (see Table 1). Now, we can conclude that when we increase the complexity of the grid, we can obtain tighter estimators of π. Wood and Robertson [] investigated this conclusion. They introduced the notion of grid density, which is the average length of grid in a unit area and showed that when the experiments are standardized, Buffon’s estimator in a single grid is the most efficient. In their approach, when d=1, the single grid has unit grid density, the double grid has grid density of two, and the triple grid has grid density of three. Hence, the standardization of experiments corresponds to r=1 in the single grid, r=1/2 in the double grid and, finally, r=1/3 in the triple grid. Replacing these values of r in equations (), (), and () and evaluating them at π=3.1416 yields the values of AVar((OverHat[π])) given in Table 2. As Wood and Robertson claimed, the tightest estimator is obtained in the single-grid experiment.

The Package

The BuffonNeedle package is designed to throw needles on single, double, and triple grids. Copy the file BuffonNeedle.m (see Additional Material) into the Mathematica ⊳ AddOns ⊳ Applications folder. The following command loads the program.

<<BuffonNeedle.m

There are three functions in this package: SingleGrid[n,r], DoubleGrid[n,r], and TripleGrid[n,r]. Here, n is the number of needles and r is the ratio of needle length to grid height (i.e., r=l/d), where n can be any integer, while r is a real number on the interval (0,1].

SingleGrid[n,r] implements a single-grid Buffon’s experiment. It gives a table showing the number and frequency ratios of the two possible outcomes, together with the theoretical probabilities and the estimate of π defined in equation (). The function also gives a picture of the simulation results. In the picture, the midpoints of the needles crossing any line are colored red, while those of needles crossing no line are colored green. The functions DoubleGrid[n,r] and TripleGrid[n,r] carry out similar processes for double and triple grids, respectively. In the picture of a double-grid experiment, the midpoints of the needles are colored green, blue, and red to show the three possible outcomes of zero, one, and two crossings. The estimate of π in this experiment is defined in equation (). As in the other two cases, in a triple-grid experiment carried out by the function TripleGrid[n,r], the four possible outcomes of zero, one, two, and three crossings are represented by four different colors of the midpoints of the needles. The estimate of π given in the table is defined in equation ().

For each grid, as n gets larger, it could be expected that the difference between the estimated and actual values of π would get smaller. You can also check some statistical results discussed in the previous sections by trying different values of r. Additionally, for large values of n, very nice pictures that exhibit the interface between randomness and determinism can be obtained. Some examples for various values of n and r are given below.

SingleGrid[10,.87]

DoubleGrid[15,.76]

TripleGrid[12,.70]

SingleGrid[10000,.91]

DoubleGrid[38000,.85]

TripleGrid[28000,1]

Calculating the Crossing Probabilities in Single and Double Grids

In this section, we show how to calculate the crossing probabilities in single- and double-grid experiments using Mathematica.

Single-Grid Probabilities

In the single-grid experiment, two independent random variables with uniform distribution are defined to determine the relative position of the needle to the lines: the distance X of the needle’s midpoint to the closest line and the acute angle α formed by the needle (or its extension) and the line (Figure ). It is seen that X can take any value between 0 and d/2 and α can take any value between 0 and π/2. The density functions of X and α are then given by

gdist1=UniformDistribution[{0,(d/2)}];
(Subscript[f, X])=PDF[gdist1,x]

gdist2=UniformDistribution[{0,(Pi/2)}];
(Subscript[f, α])=PDF[gdist2,α]

Since X and α are independent, the joint density function is the product of the density function of X alone and the density function of α alone:

(Subscript[f, X,α])=(2/d)(2/π)

for 0≤x≤d/2, 0≤α≤π/2.

From Figure , it is clear that the needle crosses the line when X≤l/2 sin α. The probability of this event is then

(Subscript[p, 1])=(∫_0^(π/2))(∫_0^((l/2)*Sin[α]))(Subscript[f, X,α]) \[DifferentialD]x\[DifferentialD]α

As there are two possible outcomes in the single-grid experiment, the probability that the needle does not cross any line is given by

(Subscript[p, 0])=1-(Subscript[p, 1])

which can alternatively be calculated by

(Subscript[p, 0])=(∫_0^(π/2))(∫_((l/2)*Sin[α])^(d/2))(Subscript[f, X,α]) \[DifferentialD]x\[DifferentialD]α

The probabilities (Subscript[p, 0]) and (Subscript[p, 1]) can be written as a function of θ and r as in equation ().

(Subscript[p, 0])/.Pi->(1/θ)/.l->r*d

(Subscript[p, 1])/.Pi->(1/θ)/.l->r*d

{{abs1,ord1},{abs2,ord2}}={{4.084193`,2.952399`},{3.461541288`,2.1699`}};
poi={3.77286`,2.5611501`};
fig4=Show[Graphics[{{Thickness[0.01`],RGBColor[1,0,0],Line[{{abs1,ord1},{abs2,ord2}}]},{PointSize[(1/45)],RGBColor[0,0,1],Point[poi]},{Thickness[0.006`],RGBColor[0.15`,0.35`,0.25`],Line[{{4.5`,1},{0.5`,1}}],Line[{{0.5`,1},{4.5`,1}}],Line[{{0.5`,2},{4.5`,2}}],Line[{{0.5`,3},{4.5`,3}}]},{Dashing[{0.01`}],Thickness[0.001`],Line[{{3.81`,2.56`},{4.2`,poi[[2]]}}]},{Dashing[{0.01`}],Thickness[0.001`],Arrow[{{3.77`,2.6`},{3.77`,3}}]},{Dashing[{0.01`}],Thickness[0.001`],Arrow[{{3.77`,3},{3.77`,2.6`}}]},{RGBColor[0,0,0],Thickness[0.003`],{Arrow[{{1,1},{1,2}}]},Arrow[{{1,2},{1,1}}],Arrow[{{1,2},{1,3}}],Arrow[{{1,3},{1,2}}]},Text[d,{0.85`,1.5`},{-1,0}],Text[d,{0.85`,2.5`},{-1,0}],Text[X,{3.7`,2.8`},{1,0}],Text["l/2",{3.52`,2.4`},{-1,0},{1,1.25`}],Text["l/2",{3.88`,2.84`},{-1,0},{1,1.25`}],Text[α,{3.88`,2.62`},{-1,0}]}],PlotRange->All,Axes->None]

Double-Grid Probabilities

In the double-grid experiment, three independent random variables with uniform distribution can be defined to determine the relative position of the needle to the lines: the distance X of the needle’s midpoint to the closest horizontal line, the distance Y of the needle’s midpoint to the closest vertical line, and the acute angle α formed by the needle and the horizontal line, as in Figure . It is seen that X and Y can take any value between 0 and d/2 and α can take any value between 0 and π/2. The density functions of X, Y, and α are given by

gdist1=UniformDistribution[{0,(d/2)}];
(Subscript[f, X])=PDF[gdist1,x]

(Subscript[f, Y])=PDF[gdist1,y]

gdist2=UniformDistribution[{0,(Pi/2)}];
(Subscript[f, α])=PDF[gdist2,α]

As in the case of the single-grid experiment, the joint density function of X, Y, and α is the product of the density functions of X, Y, and α:

(Subscript[f, X,Y,α])=(2/d)(2/d)(2/π)

for 0≤x≤d/2, 0≤y≤d/2, 0≤α≤π/2.

In the double-grid experiment, there are four possible outcomes:

The needle crosses a horizontal line while not crossing a vertical line.

The needle crosses a vertical line while not crossing a horizontal line.

The needle crosses both a vertical line and a horizontal line or, equivalently, the needle crosses two lines.

The needle crosses neither a vertical line nor a horizontal line or, equivalently, the needle crosses no line.

The needle crosses a horizontal line but does not cross a vertical line when X≤l/2 sin α and Y>l/2 cos α. The probability of this event is given by

(Subscript[p, X(OverBar[Y])])=(∫_0^(π/2))(∫_((l/2)*Cos[α])^(d/2))(∫_0^((l/2)*Sin[α]))(Subscript[f, X,Y,α]) \[DifferentialD]x\[DifferentialD]y\[DifferentialD]α

The needle crosses a vertical line but does not cross a horizontal line when X>l/2 sin α and Y≤ l/2 cos α. The probability of this event is given by

(Subscript[p, (OverBar[X])Y])=(∫_0^(π/2))(∫_0^((l/2)*Cos[α]))(∫_((l/2)*Sin[α])^(d/2))(Subscript[f, X,Y,α])\[DifferentialD]x\[DifferentialD]y\[DifferentialD]α

Thus, the probability that the needle crosses exactly one line is

(Subscript[p, 1])=(Subscript[p, X(OverBar[Y])])+(Subscript[p, (OverBar[X])Y])

The needle crosses both the vertical line and the horizontal line when X≤l/2 sin α and Y≤ l/2 cos α. The probability of this event is

(Subscript[p, 2])=(∫_0^(π/2))(∫_0^((l/2)*Cos[α]))(∫_0^((l/2)*Sin[α]))(Subscript[f, X,Y,α]) \[DifferentialD]x\[DifferentialD]y\[DifferentialD]α

Finally, the needle crosses neither the vertical line nor the horizontal line when X>l/2 sin α and Y>l/2 cos α. The probability of this event is

(Subscript[p, 0])=(∫_0^(π/2))(∫_((l/2)*Cos[α])^(d/2))(∫_((l/2)*Sin[α])^(d/2))(Subscript[f, X,Y,α]) \[DifferentialD]x\[DifferentialD]y\[DifferentialD]α

The probabilities (Subscript[p, 0]), (Subscript[p, 1]), and (Subscript[p, 2]) can be written as functions of θ and r, as in equation ().

(Subscript[p, 0])/.Pi->(1/θ)/.l->r*d//Simplify

(Subscript[p, 1])/.Pi->(1/θ)/.l->r*d//Simplify

(Subscript[p, 2])/.Pi->(1/θ)/.l->r*d//Simplify

fig5=Show[Graphics[{{RGBColor[1,0,0],Thickness[0.01`],{Line[{{3.49451`,2.8633`},{3.04665`,2.2004`}}]}},{PointSize[(1/45)],RGBColor[0,0,1],Point[{3.270582`,2.53193`}]},{RGBColor[0.15`,0.35`,0.25`],Thickness[0.005`],Line[{{1,1},{4,1}}],Line[{{1,2},{4,2}}],Line[{{1,3},{4,3}}],Line[{{1,1},{1,3}}],Line[{{2,1},{2,3}}],Line[{{3,1},{3,3}}],Line[{{4,1},{4,3}}]},{RGBColor[0.15`,0.35`,0.25`],Thickness[0.005`],Dashing[{0.01`}],Line[{{0.5`,1},{1,1}}],Line[{{0.5`,2},{1,2}}],Line[{{0.5`,3},{1,3}}],Line[{{1,3},{1,3.5`}}],Line[{{1,1},{1,0.5`}}],Line[{{2,3},{2,3.5`}}],Line[{{2,1},{2,0.5`}}],Line[{{3,1},{3,0.5`}}],Line[{{4,3},{4.5`,3}}],Line[{{4,2},{4.5`,2}}],Line[{{4,1},{4.5`,1}}],Line[{{4,3},{4,3.5`}}],Line[{{4,1},{4,0.5`}}],Line[{{3,3},{3,3.5`}}]},{RGBColor[0,0,0],Dashing[{0.01`}],Thickness[0.001`],Line[{{3.230582`,2.53193`},{3.570582`,2.53193`}}]},{Arrow[{{0.25`,1},{0.25`,2}}]},{Arrow[{{0.25`,2},{0.25`,1}}]},{Arrow[{{0.25`,2},{0.25`,3}}]},{Arrow[{{0.25`,3},{0.25`,2}}]},{Arrow[{{1,0.25`},{2,0.25`}}]},{Arrow[{{2,0.25`},{1,0.25`}}]},{Arrow[{{3,0.25`},{2,0.25`}}]},{Arrow[{{2,0.25`},{3,0.25`}}]},{Arrow[{{3,0.25`},{4,0.25`}}]},{Arrow[{{4,0.25`},{3,0.25`}}]},{Arrowheads[Small],Arrow[{{3.270582`,2.57193`},{3.270582`,3}}]},{Arrowheads[Small],Arrow[{{3.270582`,3},{3.270582`,2.57193`}}]},{Arrowheads[Small],Arrow[{{3.230582`,2.53193`},{3,2.53193`}}]},{Arrowheads[Small],Arrow[{{3,2.53193`},{3.230582`,2.53193`}}]},{Text[d,{0.05`,2.5`},{-1,0}],Text[d,{0.05`,1.5`},{-1,0}],Text[d,{1.5`,0.05`},{-1,0}],Text[d,{2.5`,0.05`},{-1,0}],Text[d,{3.5`,0.05`},{-1,0}],Text[X,{3.15`,2.8`},{-1,0}],Text[Y,{3.08`,2.65`},{-1,0}],Text[α,{3.38`,2.6`},{-1,0}]}}],PlotRange->All,Axes->None,AspectRatio->0.75`]

Delta Method and Asymptotic Variance

Let a random variable Y be a function of the random variable X, that is, Y=H(X). When the function H(X) is nonlinear, it may not be possible to compute the true mean and the true variance of Y. One can, however, calculate estimates of the true mean and true variance. The delta method is a very useful way to find such estimates [, ] and is based on the Taylor expansion about the mean of X. Let the mean of X be μ and the variance (σ^2). Then the Taylor expansion of the function H(X) about μ to the third term is

H(X)=H(μ)+H'(X-μ)+(1/2)H''(((X-μ))^2).

Taking the expectation of both sides, we obtain the approximate mean of Y=H(X) as

AMean[H(X)]=E[H(X)]=H(μ)+(1/2)H''(μ)(σ^( 2)).

From the well-known identity of statistics

Var[H(X)]=E(({H(X)-E[H(X)]})^( 2)),

the approximate variance, also called asymptotic variance, of Y=H(X) is

AVar[H(X)]=(([H'(μ)])^( 2))(σ^( 2)).

Thus, we can say that the random variable Y is distributed with the approximate mean H(μ)+1/2H''(μ)(σ^( 2)) and the approximate variance (([H'(μ)])^( 2))(σ^( 2)).

Buffon’s estimator is a nonlinear function of the random variable (OverHat[θ ]) ((OverHat[π])=H((OverHat[θ ]))=1/(OverHat[θ])); hence, the delta method can be used to find its asymptotic variance. From the previous sections, for each grid, we know E((OverHat[θ ])) and Var((OverHat[θ ])). Then, in equation (), substituting H(X)=(OverHat[π]), μ=θ, and (σ^( 2))=Var((OverHat[θ ])) from equations (), (), and (), we can obtain the asymptotic variances of (OverHat[π]) in equations (), (), and () for each grid.

Alternatively, asymptotic variances of (OverHat[π]) can be computed as follows [, ]

AVar[(OverHat[π])]=((([H'(θ)])^( 2))/(n I(θ))).

Here I(θ) is the Fisher information number, given by

I(θ)=(∑_i)((([(Subscript[p, i]^″)(θ)])^( 2))/((Subscript[p, i])(θ))),

where (Subscript[p, i])(θ) is the probability that the needle crosses i lines. These probabilities given in equations (), (), and () actually define a list for each grid as follows:

ProbSingle={1-2 r θ,2 r θ};
ProbDouble={1-r(4-r)θ,2r(2-r)θ,(r^2)θ};
ProbTriple={1+((r^2)/2)-(3/2)r(4-((Sqrt[3])/2)r)θ,-(5/4)(r^2)+(3/2)r(4-((Sqrt[3])/2)r)θ,(r^2)-((3(Sqrt[3]))/4)(r^2)θ,-((r^2)/4)+((3(Sqrt[3]))/4)(r^2)θ};

From equations () and (), one can define the following functions to obtain the asymptotic variances:

deriv[expr_,var_]:=((((D[expr,var]))^2)/((expr)));
fisherInfo[list_,var_]:=Total[Table[deriv[list[[i]],var],{i,Length[list]}]];
aVar[list_,var_]:=((((D[1/var,var]))^2)/((n fisherInfo[list,var])));

For each grid, therefore, the asymptotic variances are

aVar[ProbSingle,θ]

aVar[ProbDouble,θ]

aVar[ProbTriple,θ]

Substituting θ=1/π and factoring the expressions, we have

aVar[ProbSingle,θ]/.θ->(1/π)//Factor
aVar[ProbDouble,θ]/.θ->(1/π)//Factor
aVar[ProbTriple,θ]/.θ->(1/π)//Factor

which were previously given in equations (), (), and (), respectively. For r=1 and π=3.1416, we obtain the same values given in Table 1 for each grid.

aVar[ProbSingle,θ]/.θ->(1/π)/.r->1/.π->3.1416

aVar[ProbDouble,θ]/.θ->(1/π)/.r->1/.π->3.1416

aVar[ProbTriple,θ]/.θ->(1/π)/.r->1/.π->3.1416

For the standardized experiments of Wood and Robertson [], we can obtain the asymptotic variances given in Table 2 by substituting r=1, 1/2, and 1/3 for single, double, and triple grids, respectively.

aVar[ProbSingle,θ]/.θ->(1/π)/.r->1/.π->3.1416

aVar[ProbDouble,θ]/.θ->(1/π)/.r->(1/2)/.π->3.1416

aVar[ProbTriple,θ]/.θ->(1/π)/.r->(1/3)/.π->3.1416

Acknowledgments

I would like to thank the reviewers whose comments led to a substantial improvement in this article. I would also like to thank my colleague, M. H. Satman, for creating some of the functions in the package BuffonNeedle. Thanks to Dr. Aylin Aktukun for reading and commenting on numerous versions of this manuscript. This work is supported by the Scientific Research Projects Unit of Istanbul University.

References

G. L. Buffon, “Essai d’arithmétique morale,” Histoire naturelle, générale, et particulière, Supplément 4, 1777 pp. 685–713.

M. G. Kendall and P. A. P. Moran, Geometrical Probability, New York: Hafner, 1963.

P. Diaconis, “Buffon’s Problem with a Long Needle,” Journal of Applied Probability, 13, 1976  pp. 614–618.

R. A. Morton, “The Expected Number and Angle of Intersections Between Random Curves in a Plane,” Journal of Applied Probability, 3, 1966 pp. 559–562.

H. Solomon, Geometric Probability, Philadelphia: Society for Industrial and Applied Mathematics, 1978.

M. D. Perlman and M. J. Wichura, “Sharpening Buffon’s Needle,” American Statistician, 29(4), 1975 pp. 157–163.

G. R. Wood and J. M. Robertson, “Buffon Got It Straight,” Statistics and Probability Letters, 37,1998 pp. 415–421.

G. R. Wood and J. M. Robertson, “Information in Buffon Experiments,” Journal of Statistical Planning and Inferences, 66,1998 pp. 21–37.

M. R. Spiegel, Probability and Statistics, Schaum’s Outline of Probability and Statistics, New York: McGraw-Hill, 1975 pp. 67–68.

J. V. Uspensky, Introduction to Mathematical Probability, New York: McGraw-Hill, 1937, p. 252.

C. R. Rao, Linear Statistical Inference and Its Applications, 2nd ed., New York: Wiley & Sons, 1973.

G. W. Oehlert, “A Note on the Delta Method,” American Statistician, 46, 1992 pp. 27–29.

J. A. Rice, Mathematical Statistics and Data Analysis, 2nd ed., Belmont, CA: Duxbury Press, 1994 p. 149.

Additional Material

BuffonNeedle.m

Available at www.mathematica-journal.com/data/uploads/2011/11/BuffonNeedle.m.