(Also lomonster.com and clomont.com)

Average segment length


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


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: Math Riddle

Tags: Math Riddle Probability


Comment is disabled to avoid unwanted discussions from 'localhost:1313' on your Disqus account...