Mean distance between 2 points in a ball

(a) Two random points on the unit sphere $S^2$:

Assume the first point at the north pole ($\theta=0$) of $S^2$. Then the distance to a point at latitude $\theta$ is $2\sin{\theta\over2}$. Therefore the mean distance between the north pole and the second point is given by $${1\over 4\pi}\int_0^\pi 2\sin{\theta\over2}\cdot 2\pi\sin\theta\ d\theta={4\over3}\ .$$ $$ $$

(b) Two random points in the unit ball of ${\mathbb R}^3\ $:

Let ${\bf X}$ and ${\bf Y}$ be the two random points. Then $R:=|{\bf X}|$, $\ S:=|{\bf Y}|$, and $\Theta:=\angle({\bf X},{\bf Y})$ are independent random variables with densities $$f_R(r)=3r^2\quad (0\leq r\leq 1)\ ,\qquad f_S(s)=3s^2\quad(0\leq s\leq 1)\ ,$$ and $$f_\Theta(\theta)={1\over2}\sin\theta\quad(0\leq\theta\leq\pi)\ .$$ (Concerning $f_R$ and $f_S$ note that the volume included between $r$ and $r+dr$ is proportional to $r^2$. For $f_\Theta$ you may assume ${\bf X}$ pointing due north. The abstract surface area between $\theta$ and $\theta+d\theta$ is then proportional to $\sin\theta$, as in (a).)

It follows that the mean distance $\delta$ between ${\bf X}$ and ${\bf Y}$ is given by $$\delta=\int_0^1\int_0^1\int_0^\pi \sqrt{r^2+s^2-2rs\cos\theta}\ f_R(r) f_S(s)f_\Theta(\theta) d\theta\ ds\ dr\ .$$ The innermost integral computes to $$\eqalign{{1\over2}\int_0^\pi \sqrt{r^2+s^2-2rs\cos\theta}\ \sin\theta\ d\theta&={1\over 6rs}\bigl(r^2+s^2-2r s\cos\theta\bigr)^{3/2}\Biggr|_0^\pi \cr &={1\over 6rs}\bigl((r+s)^3-|r-s|^3\bigr)\ . \cr}$$ In the sequel we assume $s\leq r$ and compensate this by a factor of $2$. We are then left with $$\delta=\int_0^1\int_0^r 3r s(6r^2 s +2s^3)\ ds\ dr={36\over35}\ .$$ It should not be too difficult to set a similar computation up that is valid for a ball in any ${\mathbb R}^n$, $\ n\geq 2$.

The average distance can be calculated as a triple integral in polar coordinates:

We ran into it some years back when doing research on proteins (http://www.ncbi.nlm.nih.gov/pubmed/9514112).

(a) UNIT SPHERE

Let us recall the Laplace's multipole expansion formula into Legendre polynomials: $$\frac{1}{|\vec{r}'-\vec{r}|}=\frac{1}{r}\sum_{n=0}^\infty \left(\frac{r'}{r}\right)^n P_n(Y)\qquad ;r\geq r',$$

where I have denoted $Y:=\hat{r}\bullet\hat{r}'$. If we assume that $r=r'=1$, then $$\frac{1}{|\hat{r}'-\hat{r}|}=\sum_{n=0}^\infty P_n(Y)$$

The average distance between two uniformly distributed points on the unit sphere we denote as $\bar{L}$. Since $$|\vec{r}'-\vec{r}|=(\vec{r}'-\vec{r})^2/|\vec{r}'-\vec{r}|=(r'^2+r^2-2rr'Y)/|\vec{r}'-\vec{r}|,$$

we have for $r=r'=1$: $$|\hat{r}'-\hat{r}|=2(1-Y)/|\hat{r}'-\hat{r}|$$

and via Laplace's expansion therefore

$$\bar{L} = \mathbb{E}\left[|\hat{r}'-\hat{r}|\right]=\frac{1}{(4\pi)^2}\oint \oint |\hat{r}'-\hat{r}| d\Omega' d\Omega=\frac{1}{8\pi^2}\sum_{n=0}^\infty\oint \oint \left(1-Y\right) P_n(Y) d\Omega' d\Omega.$$

Due to orthogonality, we know as well that $$\oint P_n(Y)P_m(Y) d\Omega = \frac{4\pi}{2n+1} \delta_{nm}.$$

Since $1=P_0(Y)$ and $Y = P_1(Y)$, we immediately get

$$ \bar{L} = \frac{1}{8\pi^2}\sum_{n=0}^\infty\oint \oint \left(1-Y\right) P_n(Y) d\Omega' d\Omega = \frac{1}{2\pi} \oint \left(1-\frac{1}{3}\right) d\Omega = \frac{1}{2\pi}\frac{2}{3}4\pi = \frac{4}{3}$$

(b) UNIT BALL

The average distance between uniformly distrubuted two random points in the unit ball is given by$\bar{L} = \mathbb{E}\left[|\vec{r}'-\vec{r}|\right]$. Since the problem is symmetric with respect of labeling of the two points and changing their place, we can define a random variable $L' := |\vec{r}'-\vec{r}|\theta(r-r')$ which vanishes when $r'>r$. Clearly, we miss exactly half of the points, therefore

$$\bar{L} = 2\bar{L}' = 2\mathbb{E}\left[|\vec{r}'-\vec{r}|\theta(r-r')\right]=2\left(\frac{3}{4\pi}\right)^2\int \int |\vec{r}'-\vec{r}|\theta(r-r') dV' dV$$,

which can be expanded by the same trick as before. One gets $$ \bar{L} = \frac{9}{8\pi^2}\int_0^1 \int_0^r \oint \oint \left(r'^2+r^2-2rr'Y\right) \frac{1}{r} \sum_{n=0}^\infty \left(\frac{r'}{r}\right)^n P_n(Y) r^2r'^2 d\Omega' d\Omega dr' dr$$

Due to orthogonality, we get $$ \bar{L} = 18\int_0^1 \int_0^r\frac{1}{r}\left(r'^2+r^2-\frac{2}{3}rr'\frac{r'}{r}\right) r^2r'^2 dr'dr = 18\int_0^1 r^6\left(\frac{1}{5}+\frac{1}{3}-\frac{2}{3}\frac{1}{5}\right)dr = \frac{36}{35}$$

(c) BALL DISTRIBUTION FUNCTION

Denote the probablity density of $L$ as a function $g(\lambda)$, for this function one has $$g(\lambda) = \left(\frac{3}{4\pi}\right)^2\int\int \delta\left(|\vec{r}'-\vec{r}|-\lambda\right)dV'dV.$$

We expand the $\delta(L-\lambda)$ as a sum of Legendre polynomials. Note that $\delta(L-\lambda) \neq 0$ only when $|r'-r|<\lambda<r+r$ since $|\vec{r}'-\vec{r}|\in (|r'-r|,r'+r)$. Since $L = |\vec{r}'-\vec{r}|=\sqrt{r'^2 + r^2 - 2rr'Y}$ and $\delta(f(x))=\delta(x-x_0)/|f'(x_0)|$ with $f(x_0)=0$, we get

$$\delta(L-\lambda) = \frac{\sqrt{r'^2 + r^2 - 2rr' Y}}{rr'}\delta\left(Y-\frac{r'^2+r^2-\lambda^2}{2rr'}\right)=\frac{\lambda}{rr'}\delta\left(Y-\frac{r'^2+r^2-\lambda^2}{2rr'}\right).$$

Without loss of generality, we assume $r'<r$ (symmetry), $r-r'<\lambda<r+r'$ (vanishment), then

$$\delta(L-\lambda) = \sum_{n=0}^\infty A_n P_n(Y)$$

Integrating with $P_m(Y)$ over $d\Omega$ and using the orthogonality relation, we get

$$A_n = \frac{2n+1}{2} \int_{-1}^1 \delta(L - \lambda) P_n(Y) dY = \frac{2n+1}{2} \int_{-1}^1 \frac{\lambda}{rr'}\delta\left(Y-\frac{r'^2+r^2-\lambda^2}{2rr'}\right) P_n(Y) dY,$$

therefore for $|r-r'|<\lambda<r+r'$: $$\delta(L-\lambda)= \frac{\lambda}{2rr'} \sum_{n=0}^\infty (2n+1)P_n\left(\frac{r'^2+r^2-\lambda^2}{2rr'}\right)P_n(Y).$$

For the probability density therefore, with help of symmetry $$g(\lambda) = 2\lambda\left(\frac{3}{4\pi}\right)^2 \sum_{n=0}^\infty \underset{0<r-r'<\lambda<r+r'}{\int \int \oint \oint} \frac{2n+1}{2rr'} P_n\left(\frac{r'^2+r^2-\lambda^2}{2rr'}\right)P_n(Y) r'^2 r^2 d\Omega' d\Omega dr' dr$$

which is due to orthogonality (only $n=0$ survives) $$g(\lambda) = 9\lambda \int_{\lambda/2}^{1} \int_{\lambda-r}^{r}\!\! r'r\, dr' dr = \frac{9}{2}\lambda^2 \int_{\lambda/2}^{1} \!\! 2r^2\!-\!\lambda r\, dr = \frac{3}{16}\lambda^2 (2-\lambda)^2(4+\lambda)$$