# 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#:

 1 2 3 4 5 6 7 8 9  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: