The area of the circle

Random points

Consider the following script pi.kts:

val random = java.util.Random()

fun pi0(n: Int) {
  for (i in 1..n) {
    val x = random.nextDouble()
    val y = random.nextDouble()
    println("Sample $i: ($x, $y)")
  }	  
}

val n = if (args.size == 1) args[0].toInt() else 100
pi0(n)

It generates \(n\) random points in the unit square \([0,1]^{2}\), for instance:

$ kts pi.kts 5
Sample 1: (0.5907325188943173, 0.09132164982022872)
Sample 2: (0.09249216713952713, 0.5679322425654839)
Sample 3: (0.62714728132915, 0.1969654689089445)
Sample 4: (0.49092573545422113, 0.23313234130565952)
Sample 5: (0.3615277810877552, 0.8278623847513092)

Counting points

Modify the script so that it counts how often the generated point has distance at most one from the origin. In other words, count the point if \(x^{2} + y^{2} \leq 1\).

After the for loop, print out the number of hits, and four times the ratio between the hits and the number of samples, like this:

$ kts pi.kts 200
Among 200 samples there were 162 hits
 4 * (162/200) = 3.24
$ kts pi.kts 2000
Among 2000 samples there were 1565 hits
 4 * (1565/2000) = 3.13

Can you explain this result? The probability of a point having distance at most one from the origin is exactly the area of a quarter circle of radius one.

How large do you need to make \(n\) to get more digits of \(\pi\) correct?

Buffon's needle

Let's try another experiment, the famous Buffon's needle experiment.

Suppose we have a floor made of parallel strips of wood, each of width one, and we drop a needle of length one onto the floor. What is the probability that the needle will lie across a line between two strips (instead of fitting completely inside a strip)?

Let's assume the strips are horizontal, so the lines are \(y = 0\), \(y = 1\), \(y = 2\), and so on. We will randomly generate a lower endpoint \((x_0, y_0)\) and an upper endpoint \((x_1, y_1)\) of the needle.

Since the strips are infinitely long, we actually don't need \(x_0\) and \(x_1\), and since the pattern repeats with period one, it suffices to generate \(y_0\) in the range \([0,1]\):

  val y0 = random.nextDouble()

To model the correct probability distribution of the needle, we need to generate a random angle of the needle, for instance in the range \([0, 90^{\circ}]\), in radians. We can do this as follows:

  val angle = random.nextDouble() * Math.PI / 2.0

The \(y\)-coordinate of the upper endpoint will then be

  val y1 = y0 + Math.sin(angle)

Finally, the needle crosses the line \(y = 1\) if and only if \(y_{1} > 1\) (and this is the only line that the needle can cross).

Count how often this happens. Print out \(2n/k\), where \(k\) is the number of line crossings observed:

$ kts pi.kts 2000000
Among 2000000 needle throws, 1273320 crossed a line
 2 * (2000000/1273320) = 3.1413941507240914

Buffon's needle, revisited

Actually, we cheated a bit in the previous experiment. When we created a random angle, we actually used the constant Math.PI. It's a bit unfair to use this constant when our goal is to estimate it, right?

So let's try to redo the experiment without using \(\pi\) and without trigonometry.

To obtain a random needle direction in the range \([0, 90^{\circ}]\), we generate a random point \((x,y)\) in the unit square \([0,1]^{2}\). As long as \(x^{2} + y^{2} > 1\), we discard the point and generate a new one. Then compute the norm \(m = \sqrt{x^{2} + y^{2}}\).

Now generate a random bottom end \(y_0\) for the needle, in the range \([0,1]\), and compute the upper end of the needle as \(y_1 = y_0 + y/m\).

Repeat Buffon's needle experiment with this method.