Wednesday, January 5, 2011

Eclipse

Yesterday much of Europe, North Africa and central Asia had an opportunity to observe the first solar eclipse of 2011. NASA reports that in Kiev the obscuration (fraction of the Sun’s surface occulted by the Moon) reached the value of 0.731.

What we see during the eclipse is two overlapping circles of nearly equal radius.

Assuming the radius is the same, what is the observable distance between the centers of those circles?

We could compute this using pure geometry, but let’s make it a statistical exercise.

What we know by now is that the intersection of the circles makes a known fraction of each of them (0.73). Let’s drop points into the first circle at random. Then, on average, 73% of them will also fall into the second circle.

Let d be the distance between the centers. For simplicity let the both circle have the radius 1, the first circle have the center (0,0) and the second circle have center (d,0).

The probability of a random point from the first circle to be inside the second circle equals

P{distance((x,y),(d, 0)) ≤ 1} = P{(x-d)2+y2 ≤ 1} = P{d ≤ x + √(1-y2)} = 1-P{x + √(1-y2) ≤ d}.

The last probability is just the cumulative distribution function for the random variable x+√(1-y2). And the value of d we want to find is the 1–0.73=0.27-quantile of that CDF.

Now let’s make an estimate of that quantile using R. To generate uniformly distributed random points inside the circle we just generate random points inside the square [–1;1]2 and then use only those that are inside the circle.

> x  <-  runif(10000, min=-1, max=1)
> y  <-  runif(10000, min=-1, max=1)
> i  <-  x^2 + y^2 <= 1
> length(1[i])
[1] 7861

(Note that the fraction of the points inside the circle, 7861/10000, gives an approximation to the ratio of the areas of the circle and the square, π/4.)

For each of these points, compute the random value we’re interested in

> d  <-  x[i] + sqrt(1 - y[i]^2)

and estimate the quantile, which is basically the largest of the lowest 27% of our observations.

> quantile(d, 1 - 0.73)

  27% 
  0.4305514 

For comparison, the true value of d is 0.42737.