# Monte-Carlo estimation of the mean chord length in a polygon

This is a more "statistical" version of my Math.SE problem .

Question Let $\mathbf{P}$ be a convex polygon of $m$ sides. Pick two of its edges at random, and further pick a random point on each of these edges, and compute the length of the chord we get. If we repeat this process $n$ times,and average the chord lengths we thus obtained, we get an estimate $\hat{\ell}_n$ of the expected chord length $\ell$ of $\mathbf{P}$. Then,how large should $n$ be to ensure that $\hat{\ell}_n$ differs from $\ell$ by (say) at most 5%,with a confidence level of (say) 95%, and precisely how should we collect our samples?

Note that determining the expected chord length analytically leads to intractable integrals as shown in my SE post. Also,note that there are multiple inequivalent ways to define the expected chord length of $\mathbf{P}$; the one we need here is $\ell = \lim_{n\rightarrow \infty} \hat{\ell}_n$.

-
I've tried to edit your post to (a) make the math slightly more readable and (b) change what I believe were some typos in the math. Please make sure I have not messed things up in the process. –  cardinal Apr 7 '11 at 2:16
Also, what do you mean by "ensure"? Strictly speaking, Monte Carlo can't (quite) provide such guarantees. –  cardinal Apr 7 '11 at 2:16
@cardinal: Thanks for your editing. I think I have fixed the guarantees part. –  PKG Apr 7 '11 at 2:24
The integrals do not appear to be "intractable"; indeed, they appear to be solvable with elementary functions. The mean distance from a point to a line segment involves linear expressions times quadratic expressions and quadratic expressions times logarithms of linear plus quadratic expressions. Integrating the latter by parts gives rational functions of quadratic algebraic expressions. Those all can be done by completing the square and trig substitutions. –  whuber Apr 7 '11 at 5:18
E.g., 48 times the mean distance between two non-parallel, non-adjacent sides of a regular hexagon of side 2 equals $32-48 \sqrt{3}-22 \sqrt{7}+46 \sqrt{13}+27 \text{ArcCsch}(3 \sqrt{3})-108 \text{ArcSinh}(\frac{1}{3 \sqrt{3}})+81 \text{ArcSinh}(\frac{5}{3 \sqrt{3}})$ + $66 \text{ArcSinh}(\sqrt{3})-27 \text{Log}(3)-18 \text{Log}(27)+36 \text{Log}(432)+36 \text{Log}(7-4 \sqrt{3})+6 \text{Log}(-135+78 \sqrt{3})$ + $48 \text{Log}(2+\sqrt{7})-72 \text{Log}(-\frac{2}{3} (2+\sqrt{7}) (-7+\sqrt{13}))-24 \text{Log}(-1+\sqrt{13})+24 \text{Log}(1+\sqrt{13})$ (yielding a mean of about 2.8331). –  whuber Apr 7 '11 at 5:25

This is really a question about numerical integration because it reduces immediately to the problem of finding the mean distance between two line segments in the plane.

Monte-Carlo integration would be wasteful and inefficient compared to other methods that are readily available, including

1. Direct analytic integration. This is messy but can be done.

3. Analytic integration of the mean distance between a point and a line segment and numeric quadrature of that mean distance.

4. Discrete approximations, such as Riemann sums, the Trapezoidal Rule, etc.

Because the mean distance function in (3) is differentiable and varies relatively slowly, simple numerical methods such as Simpson's Rule work well. It's best to do the numeric integration over the shorter of the two line segments. In experiments I find that using the endpoints and midpoint in Simpson's Rule typically gives better than 1% accuracy and using five equally spaced points is far better than that. (The toughest cases occur where two sides are equally long, parallel, and close together.) If you need to, you can do the analysis to estimate the Simpson's Rule error and adaptively decrease the point spacing.

For the original problem, the approximations sometimes overestimate and sometimes underestimate the mean distance, so we could expect some cancellation of errors.

At any rate, you will definitely achieve 5% accuracy using method (3). For an $n$ sided polygon it requires $n$ computations (each of constant cost) per vertex and $n-1$ computations per midpoint for a cost of $O(n^2)$. The individual computations of mean point-segment distances require some basic arithmetic, four square roots, and two logarithms, so they're quick and easy. (Mathematica computes about 7,500 mean segment-segment distances per core per second.)

If you want to avoid all but the simplest arithmetic, you can approximate each segment-segment integral with a double Simpson's Rule pass, one for each segment (method 4). With three division points (requiring the computation of nine distances) you will achieve 5% accuracy. With five division points (25 distances) you will be better than 1% accurate (usually better than 0.05%) except for extremely close parallel segments. In Mathematica, with its interpreted overhead, this technique actually is only half as fast as method (3).

This plot shows (negative) log (base 10) relative errors for 1000 segments randomly distributed near the endpoint of a unit segment (which is the larger of the two). E.g., a value of 2 is a 1% error; a value of 6 is a 0.0001% error. This is a severe test, because the worst performances occur when the segments are crossing or close to doing so, situations that are impossible or unlikely in convex polygons. (Note that the smallest possible mean distance is 1/4 (for an infinitesimally short segment near the middle of the unit segment).)

Where the approximation is high, the error is plotted on the right, and where the approximation is low, the error is plotted on the left. Blue dots are the three-point by three-point Simpson's Rule calculation (method 4) and red dots are the three-point Simpson's Rule calculation (method 3).

Typically method 3 (red) has much better than 2% accuracy and is better than method 4 (blue), which in a few cases has almost 10% error. Method 3 has a tendency to overestimate slightly (more of its values are at the right). Both methods almost always obtain better than 0.1% accuracy (a height of 3 or better on the plot), especially at mean distances greater than the longer side. This implies that if you want to improve the accuracy, focus first on nearby edges of a convex polygon.

In this plot of the 1000 segments, the brighter, thicker, redder ones have the greatest relative error for Method 3. The black segment is the reference (unit) segment.

-
Many thanks for the extremely detailed response. –  PKG Apr 7 '11 at 17:40