Disturbing MATLAB Accuracy in Monte Carlo Simulation

Question

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.

NOTE: Please check pages 25-30 in the Google Books version here

Answer 1

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.