# How to improve the accuracy of this Monte Carlo simulation

A friend of mine went on an interview and he was asked to calculate the angle between the hour hand and minute hand of an analog clock when the time was 3:15.

Arguably, this is a trivial problem that can be easily solved analytically. 15 minutes is the 1/4th of an hour, and an hour is the 1/12th of a full circle, therefore the angle is $(1/4)\times(1/12)\times 360 = 7.5 \deg$. I was wondering whether a Monte Carlo simulation could do well.

EDITED: As suggested by @kuba here is the image of a "real" clock hopefully at 3:15 precisely.

My code that follows operates on the mask of an ideal clock at 3:15 (angle between two hands at 7.5 degrees) and not on a mask generated by the above image. It takes more than $10^6$ points to derive an answer of $\sim 7 \deg$. Could the experts here improve the accuracy of the simulation?

ClearAll["Global*"];
(* Generate some random points. Use 0.98 instead of 1.00 or otherwise
some points will be half-inside the circle and half-outside.
This will hurt our results in left untreated. *)
points = Show[
Graphics[{Pink,
Point /@
Select[
Partition[RandomReal[{-1, 1}, 10000], 2],
({a, b} = #; a^2 + b^2 <= 0.98) &]}]]


(* In reality the mask below should have been generated by a real clock *)
disk = Graphics[{Pink, Disk[{0, 0}, 1, {0, 2π - 7.5π/180}]}]


finalImage = ImageSubtract[disk, Show[disk, points]]


(* Count the number of points inside the slice and outside the slice.
If I use "Count" as the property in ComponentMeasurements[], I get worse results.
EDITED: As noted by @MichaelE2 the following statement for 'outArea'
might calculate the whole rectangular area and not only the points! *)
inArea  = Last /@ ComponentMeasurements[finalImage, "Area"] // Total;
outArea = Last /@ ComponentMeasurements[points, "Area"] // Total;

(inArea/outArea) * 360 // N


EDITED: With the point set by @DanielLichtblau, I get $\sim 7.3$ degrees with $10^4$ points! finalImage looks like this:

-
The problem as I see it has two parts: 1) create the mask from a real clock and 2) generate the random points and calculate how many fall into the slice. I used an ideal mask for the 1st part of the solution. I did that because I wanted to check a la proof-of-concept if counting points in the 2nd part would work. If it didn't with the ideal slice, then it certainly wouldn't with the slice from a real clock. As it turns out even with an ideal slice it converges very slowly. (Any solution to the problem involving randomness is accepted.) –  Zet Oct 7 '13 at 10:40
Why generate two dimensional points when all you really require is one dimension along the circumference of the clock? You don't need to consider the whole clock but just a segment of the clock you know your small angle fits within, say 5 minutes which would improve efficiency by a factor of 12. –  Ymareth Oct 7 '13 at 11:11
Quick and dirty approach to a low discrepancy (uniformly distributed) point set: pts = Point[ Select[Flatten[ Table[2*Mod[{j*GoldenRatio, k*GoldenRatio}, 1.] - 1, {j, 2*100}, {k, 2*100}], 1], #[[1]]^2 + #[[2]]^2 < .99 &]]; –  Daniel Lichtblau Oct 7 '13 at 15:04
@DanielLichtblau that's awesome! With your point set I get ~7.3 degrees for 10^4 points! Congratulations! –  Zet Oct 7 '13 at 15:41

Could work with a low discrepancy rather than pseudorandom sequence. These tend to be better for avoiding approximation errors associated with "clumping" (and the corresponding "empty regions") that one gets with the latter. The code below will provide a sequence that is uniformly distributed modulo 1 in the unit square and seems to be low discrepancy, although there are, I believe, more refined ways of getting such a sequence. That is to say, one can almost certainly improve on the quality.

pts = Point[ Select[Flatten[Table[
2*Mod[{j*GoldenRatio, k*GoldenRatio}, 1.] - 1, {j, 2*100}, {k, 2*100}], 1],
#[[1]]^2 + #[[2]]^2 < .99 &]]


The process of selecting those that lie inside the unit circle could also introduce some smallish bias but I would still expect these to give better coverage of the disk than a set of pseudorandom points.

For reasons unclear I was not able to get quite the right pictures, or numerical values, using the code in the question, so I'll not try to show those computations.

-

The problem seems to be that

outArea = Last /@ ComponentMeasurements[points, "Area"] // Total;


estimates the area of the whole square image, not of the disk or the points.

For instance, with

SeedRandom[1];
points = Show[
Graphics[{Pink,
Point /@
Select[Partition[RandomReal[{-1, 1}, 100000],
2], ({a, b} = #; a^2 + b^2 <= 0.98) &]}]]


I get

inArea = Last /@ ComponentMeasurements[finalImage, "Area"] // Total;
outArea = Last /@ ComponentMeasurements[points, "Area"] // Total;

(inArea / outArea) * 360 // N
(* 5.93865 *)


If I adjust for the ratio of the areas of the square to the circle, I get

4 / Pi (inArea / outArea) * 360 // N
(* 7.56132 *)


which is much closer to the theoretical 7.5. The outArea is

outArea
(* 129419. *)


The image size is 360 and 360^2 is 129600.

For a Monte Carlo simulation, it seems to me to be appropriate to compute the proportion of points in the circle that land in the wedge.

trialpts[n_] :=
Pick[#, UnitStep[0.98 - Total[#^2, {2}]], 1] &@ RandomReal[{-1, 1}, {n, 2}];

est[pts_] :=
2. Pi Count[UnitStep[-#] UnitStep[# + 7.5 Degree] &[ArcTan @@ Transpose@pts],
1] / Length[pts] / Degree


Example

SeedRandom[1];
Table[est[trialpts[10^k]], {k, 3, 7}]
(* {9.43644, 6.6289, 7.20019, 7.60667, 7.53296} *)


The 95% confidence level margins of error for these sample sizes are respectively

p = (1/4)*(1/12);
Table[1.96 360 2 Sqrt[p (1. - p)/n] /. n -> 10^k, {k, 3, 7}]
(* {6.37377, 2.01556, 0.637377, 0.201556, 0.0637377} *)


These are proportional to 1/Sqrt[n]`, and that's as good as one should expect from a uniformly distibuted random sample.

-
Nice catch @MichaelE2! Thanks! This would also explain why @DanielLichtblau wasn't able to reproduce my numerical results. –  Zet Oct 9 '13 at 15:26