Lomont.org

(Also lomonster.com and clomont.com)

Average segment length

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

Tags: Math Riddle Probability

Comments

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