# Disturbing MATLAB Accuracy in Monte Carlo Simulation

I've been reading Nahin's book Inside Interesting Integrals and got up to the point where he is calculating the probability that a circle $C_2$ is entirely contained in circle $C_1$, given that three points are chosen at random (which are contained in $C_1$) that uniquely determine $C_2$. The analytic proof makes sense, but there is a contradiction between that number and the number achieved using his Matlab Monte Carlo Simulation Code.

Why the discrepancy considering that the code executes this a million times, and no matter how many times I do it, the numbers fluctuate somewhat between 3.9992 and 4.008, all underestimates of 0.418879...?

The math checks out for me, but not so much the code and the way a computer works with these simulations. If any can provide some insight to why this is, that would be very much appreciated.

• What if we don't have access to all the pages due to preview restrictions in Google Books? Please make an effort to make your question self-contained. – Rahul Jan 19 '18 at 9:52
• I don't think it's very likely they use pseudo-random numbers of such low quality. But the formulas for determining the center and radius of the circumcircle of three points aren't exactly trivial, so are you sure there's no mistake? – Professor Vector Jan 19 '18 at 9:56
• I am nearly certain that the correct probability is $2/5$ and not $2\pi/15$ as derived mathematically, as I have performed simulations using two independent methods and my results match your empirical results. – heropup Jan 20 '18 at 1:11
• I wrote a Monte Carlo simulation in Python with a result of 486,240 successes in $10^6$ trials, taking care that the three points are uniformly distributed inside the unit circle. I spot-checked my geometrical calculations (using Geogebra) and they appear correct. I suspect the discrepancy with your book is due to the problem of generating uniformly distributed points within a disk--it's easy to make a mistake, as noted by @heropup. – awkward Jan 21 '18 at 12:16

As far as I am able to tell, the author makes a critical error regarding the probability calculation: the distribution of the circumcenter's distance from the origin (a quantity he calls $r$) is not necessarily uniform with respect to the uniform selection of the vertices of the triangle in the disk, therefore the integration with respect to $r$ fails to account for this.

My simulation code is:

a[p_] := Det[Append[#, 1] & /@ p]
d[p_] := Join[{Total[#^2]}, #, {1}] & /@ p
bx[p_] := -Det[Delete[#, 2] & /@ d[p]]
by[p_] := Det[Delete[#, 3] & /@ d[p]]
c[p_] := -Det[Delete[#, 4] & /@ d[p]]
cc[p_] := -{bx[p], by[p]}/(2 a[p])
cr[p_] := Sqrt[bx[p]^2 + by[p]^2 - 4 a[p] c[p]]/(2 Abs[a[p]])

Manipulate[Graphics[{Line[{pt1, pt2, pt3, pt1}],
Circle[cc[{pt1, pt2, pt3}], cr[{pt1, pt2, pt3}]]}, Axes -> True],
{{pt1, {0, 0}}, Locator}, {{pt2, {1, 0}}, Locator}, {{pt3, {0, 1}}, Locator}]

L[p_] := And @@ (#.# <= 1 & /@ p)
R[n_] := RandomReal[{-1, 1}, {10^n, 3, 2}];
Parallelize[{Boole[Sqrt[cc[#].cc[#]] + cr[#] <= 1], Boole[L[#]]} & /@ R[6]] // Tally


Then count the number of results for which the output is {1,1}, divided by the number of outputs that are not {0,0}.

I also did the simulation in polar coordinates, but the above, although crude, is more transparent to the reader. The Manipulate[] shows that the formulas for the circumcircle are correct. The function L tests whether the vertices of the triangle all lie in the unit circle. The function R[n] generates $10^n$ realizations of three vertices selected uniformly at random in the square $[-1,1]^2$. The final function calculates whether each triangle satisfies the circumcircle criterion, and whether the triangle satisfies the unit circle criterion.