# Lomont.org

(Also lomonster.com and clomont.com)

Published

What's the average length of a segment chosen uniformly by area on a unit square?

# Average segment length

Here is an elementary derivation of the average length of a segment chosen randomly on the unit square. By random I mean $$x$$ and $$y$$ coordinates are chosen uniformly in $$[0,1]$$. Let the endpoints be $$(x_1,y_1)$$ and $$(x_2,y_2)$$. Then the average is given by

$\int_0^1 \int_0^1 \int_0^1 \int_0^1 \sqrt{(x_1-x_2)^2 + (y_1-y_2)^2} dx_1 dx_2 dy_1 dy_2$

which unfortunately is difficult to evaluate directly, so I'll use the deltas $$\Delta x = x_1-x_2$$ and $$\Delta y = y_1-y_2$$ instead. Now I need the probability distributions on the $$\Delta$$s to evaluate the resulting integral. The distribution of the difference of two uniform random variables is the triangle distribution. Here is a quick derivation:

Let $$X_1$$ and $$X_2$$ be uniform random variables on $$[0,1]$$. They have probability distribution functions (PDFs) of $$f_{X_i}(x)=1$$ and cumulative probability distribution functions (CDFs) of $$F_{X_i}(x)=x$$. The joint PDF of $$X_1$$ and $$X_2$$ is

$f_{X_1,X_2}(x_1,x_2)=1 \qquad 0\le X_1,X_2 \le 1$

Let $$Z=X_1-X_2$$ be the random variable of the difference. Using the method of cumulative probability functions,

\begin{aligned} F_Z(z) &= P(Z\le z) \\ &= P(X_1-X_2\le z) \\ &= \begin {cases} \int_0^{1+z} \int_{x_1-z}^1 f_{X_1,X_2}(x_1,x_2) dx_2 dx_1, & -1\le z < 0 \\ 1-\int_{z}^{1} \int_0^{x_1-z} f_{X_1,X_2}(x_1,x_2) dx_2 dx_1, & 0\le z \le 1 \end{cases} \\ &= \begin {cases} z^2/2 + z + 1/2, & -1\le z < 0 \\ -z^2/2 + z + 1/2, & 0\le z \le 1 \end{cases} \end{aligned}

Differentiating w.r.t $$z$$ then gives the PDF

$f_Z(z) = \begin {cases} 1+z, & -1\le z < 0 \\ 1-z, & 0\le z \le 1 \end{cases}$

which simplifies to $$f_Z(z)=1-|z|$$ for $$-1\le z \le 1$$.

To evaluate the average length, let $$x=x_1-x_2$$ and $$y=y_1-y_2$$, then the average length of a random segment is given by

$\int_{-1}^1\int_{-1}^1 \sqrt{x^2+y^2}(1-|x|)(1-|y|) dx dy$

Over the integration region of the square $$[-1,1]\times[-1,1]$$ this is symmetric over the four quadrants, so integrating over $$[0,1]\times[0,1]$$ obtains the value $$I$$ as

$I=4\int_0^1\int_0^1 \sqrt{x^2+y^2}(1-x)(1-y) dx dy$

Switch to polar coordinates, $$x=r \cos\theta$$, $$y=r\sin\theta$$, $$\sqrt{x^2+y^2}=r$$, using $$xy$$ symmetry, integrate over the region $$0\le\theta\le\pi/4$$ and $$0\le r\le \sec\theta$$. Note the right side of the square, $$x=1$$, gives the $$r$$ bound $$x=r\cos\theta=1$$, i.e., $$r=1/\cos\theta=\sec\theta$$. Remembering the Jacobian for the change of variables gives $$dx dy=r dr d\theta$$. The symmetry gives another factor of 2 in front, resulting in

\begin{aligned} I &= 2\times 4 \int_0^{\pi\over 4} \int_0^{\sec\theta} r (1-r\cos\theta)(1-r\sin\theta) r dr d\theta \\ &= 8 \int_0^{\pi\over 4} \int_0^{\sec\theta} r^2(1-r\cos\theta-r\sin\theta+r^2\sin\theta\cos\theta) dr d\theta \\ &= 8 \int_0^{\pi\over 4} \left. \frac{r^3}{3}-\frac{r^4}{4}\cos\theta-\frac{r^4}{4}\sin\theta+\frac{r^5}{5}\sin\theta\cos\theta \right|_0^{\frac{1}{\cos\theta}} d\theta \\ &= 8 \int \frac{1}{12\cos^3\theta} - \frac{\sin\theta}{20\cos^4\theta}d\theta \\ &= \frac{2}{3}\int \frac{1}{\cos^3\theta} d\theta + \frac{2}{5}\int \frac{-\sin\theta}{\cos^4\theta}d\theta \\ &= I_1 + I_2 \end{aligned}

Do the easier integral $$I_2$$ first via the substitution $$u=\cos\theta$$, $$du=-\sin\theta d\theta$$.

\begin{aligned} I_2 &= \frac{2}{5}\int_1^{\frac{1}{\sqrt{2}}} \frac{du}{u^4} \\ &= \frac{2}{5}\left.\frac{u^{-3}}{-3}\right|_1^{\frac{1}{\sqrt{2}}} \\ &= \frac{2}{15}(1-2\sqrt{2}) \end{aligned}

Evaluate $$I_1$$ using the substitution $$v=\sin\theta$$, $$dv=\cos\theta d\theta$$, and the identity $$\cos^2\theta=1-\sin^2\theta$$:

\begin{aligned} I_1 &= \frac{2}{3}\int_0^{\pi\over 4} \frac{d\theta}{\cos^3\theta} \\ &= \frac{2}{3}\int_0^{\pi\over 4} \frac{\cos\theta d\theta}{\cos^4\theta} \\ &= \frac{2}{3}\int_0^{1\over \sqrt{2}} \frac{dv}{(1-v^2)^2} \end{aligned}

Split the integrand using $$1-v^2=(1+v)(1-v)$$ and the method of partial fractions

\begin{aligned} \frac{1}{(1-v^2)^2} &= \frac{1}{(1-v)^2(1+v)^2} \\ &= \frac{A}{1-v}+\frac{B}{(1-v)^2}+\frac{C}{1+v}+\frac{D}{(1+v)^2} \end{aligned}

for some to be determined constants $$A$$, $$B$$, $$C$$, and $$D$$. Clearing fractions yields the identity

$1=A(1-v)(1+v)^2+B(1+v)^2+C(1-v)^2(1+v)+D(1-v)^2$

This can be expanded and powers of $$v$$ set equal on each side, and the resulting linear system solved. A useful trick here is to plug in values for $$v$$ to get enough relations.

\begin{aligned} v &=+1 &\implies& 1= 0A + 4B + 0C + 0D \\ v &=-1 &\implies& 1= 0A + 0B + 0C + 4D \\ v &=+2 &\implies& 1=-9A + 9B + 3C + D \\ v &=-2 &\implies& 1= 3A + B - 9C + 9D \end{aligned}

Lines 1 and 2 give $$B=D=1/4$$. Plugging into lines 3 and 4 gives

\begin{aligned} 0 &= \frac{3}{2} - 9A + 3C \\ 0 &= \frac{3}{2} + 3A - 9C \end{aligned}

Subtracting line 1 from line 2 gives $$0=12A-12C$$ so $$A=C$$. Plugging into line 1 and solving gives $$A=C=1/4$$. Plugging into $$I_1$$ and pulling the $$1/4$$ out front gives

\begin{aligned} I_1 &= \frac{1}{6}\int_0^{1\over\sqrt{2}} (1+v)^{-2}+(1-v)^{-2}+(1+v)^{-1}+(1-v)^{-1} dv \\ &= \frac{1}{6} \left(\left. \frac{(1+v)^{-1}}{-1} - \frac{(1-v)^{-1}}{-1} + \ln|1+v| - \ln|1-v| \right|_0^{1\over\sqrt{2}}\right) \\ &= \frac{1}{6}\left(\frac{1}{1-{1\over\sqrt{2}}}-\frac{1}{1+{1\over\sqrt{2}}} -1 + 1 + \ln\frac{1+{1\over\sqrt{2}}}{1-{1\over\sqrt{2}}}-\ln\frac{1}{1} \right) \\ &= \frac{\sqrt{2}}{3} + \frac{1}{6} \ln\left(3+2\sqrt{2}\right) \end{aligned}

Combining $$I_1$$ and $$I_2$$ gives the final answer

\begin{aligned} I &= I_1+I_2\\ &= \frac{\sqrt{2}}{3} + \frac{1}{6} \ln\left(3+2\sqrt{2}\right) + \frac{2}{15}(1-2\sqrt{2}) \\ &= \frac{2+\sqrt{2}}{15} +\frac{\ln(3+2\sqrt{2})}{6} \end{aligned}

As a sanity check simulate it. Note

$\frac{2+\sqrt{2}}{15} +\frac{\ln(3+2\sqrt{2})}{6} = 0.5214....$

Here is a simulation in F#:

open System
let r = System.Random()
let count = 100000
let mutable sum = 0.0
for j in 1..count do
let x1,x2,y1,y2 = r.NextDouble(),r.NextDouble(),r.NextDouble(),r.NextDouble()
let dx, dy = x1-x2, y1-y2
sum <- sum + sqrt (dx*dx + dy*dy)
printfn "%f" (sum / (float count))

returns 0.520409, agreeing with the symbolic result. (It can be tested online, for example, at http://tryfs.net/).

Categories:

Tags: