# Estimate the expected distance between two random points on the unit $n$-sphere [duplicate]

What is the best approach to estimate, with Wolfram Mathematica, the expected Euclidean distance in a $$(n+1)$$-dimensional space between two points selected uniformly at random on the unit $$n$$-sphere? The approach I have in mind uses an expression whose length is proportional to $$n$$, while I would like a simpler and more elegant approach.

• To clarify, do you mean distance in Euclidean n+1 space or distance on the sphere itself? Nov 21, 2020 at 16:18
• Are the points inside the sphere or on the surface? Nov 21, 2020 at 16:23
• @DanielHuber An $n$-sphere is the surface of an $(n + 1)$-dimensional ball. Nov 21, 2020 at 17:07
• @DanielLichtblau I am referring to the Euclidean distance in the $(n+1)$-dimensional space (not the great-circle distance). Nov 21, 2020 at 17:09
• Disregard my comment, I was not thinking it through clearly. Nov 22, 2020 at 15:01

Not an exact answer but a Monte-Carlo way of checking the exact answers.

Generate a random point on the unit $$n$$-sphere:

P[n_Integer?Positive] := Normalize[RandomVariate[NormalDistribution[], n]]


Measure the mean distance between a random point $$P_0$$ and another random point on the unit $$n$$-sphere, by averaging over $$m$$ random points:

M[n_Integer?Positive, m_Integer?Positive] := With[{P0 = P[n]},
Mean[Table[Norm[P[n] - P0], {m}]]]


Try for different values of $$n$$:

M[1, 10^6]
(*    0.998648    *)


The result is 1.

M[2, 10^6]
(*    1.27374    *)


This matches @flinty's result of $$4/\pi$$.

M[3, 10^6]
(*    1.33315    *)


This matches @flinty's result of $$4/3$$.

More values:

Table[M[n, 10^6], {n, 1, 10}]
(*    {0.998648, 1.27374, 1.33315, 1.35903, 1.37166,
1.37969, 1.38504, 1.38929, 1.39232, 1.39459}    *)


## Update

I think the exact answer is

d[n_] = 2^(n-1)*Gamma[n/2]^2/(Sqrt[π]*Gamma[n-1/2])


For large $$n$$ the mean distance is therefore approximately

Series[d[n], {n, ∞, 1}]

(*    Sqrt[2] - 1/(4 Sqrt[2] n) + O(1/n)^2    *)

• Thank you @Roman . Do you know how to extend this solution for a unit $n$-hemisphere (instead of a unit $n$-sphere)? Nov 21, 2020 at 16:24
• I think that's much harder because hyperspherical symmetry is broken. Nov 21, 2020 at 16:33
• Yes I think that it should be a new question, it's too different from this one. Nov 21, 2020 at 17:18
• Well done getting what appears to be the exact form. How did you manage it? Nov 21, 2020 at 18:01
• @flinty The distance between two points on the $n$-hypersphere is $2\sin(\theta/2)$ in terms of their angle (scalar product) $\vec{a}\cdot\vec{b}=\cos\theta$. Averaging over the Jacobian $J_n\propto\sin^{n-2}\theta$: Table[Integrate[2*Sin[θ/2]*Sin[θ]^(n-2), {θ, 0, π}]/Integrate[Sin[θ]^(n-2), {θ, 0, π}], {n, 2, 10}] and then use FindSequenceFunction[%, n-1] to discover the formula for arbitrary $n$. Nov 21, 2020 at 18:15

For $$n=3$$:

The PDF is $$f(d) = d/2$$ or in Mathematica TriangularDistribution[{0, 2}, 2] - a ramp shaped distribution. We can test this numerically and we get a high $$p$$-value of about 0.31 so it's a good fit:

points = RandomPoint[Sphere[], {1000000, 2}];
distances = EuclideanDistance @@@ points;
testdist = TriangularDistribution[{0, 2}, 2];
DistributionFitTest[distances, TriangularDistribution[{0, 2}, 2]]
Show[Histogram[distances, 1000, "PDF"],
Plot[PDF[testdist, x], {x, 0, 2}], Plot[d/2, {d, 0, 2}]]


The expected distance is Mean[testdist] which gives $$4/3$$. Or you can do this yourself as an integral

$$\int_{0}^{2} x\cdot\frac{x}{2} dx = \frac{4}{3}$$

For $$n=2$$:

You can integrate around the circle to find the average holding one point fixed.

Integrate[
EuclideanDistance[{0, 1}, {Cos[θ], Sin[θ]}],
{θ, 0, 2 π}]/(2 π)

(* result: 4/Pi *)

• Thank you @flinty ! My problem is to find an estimation when the number of dimensions $n$ is very large, without writing a list of constraints, i.e., one constraint per dimension when $n\gg1$. How can we solve this problem? Maybe it is worth to create a list of constraints first, and then use it with a wolfram command. Nov 21, 2020 at 16:17
• For other dimensions, the approximate distribution looks a bit like TransformedDistribution[2 x, x \[Distributed] BetaDistribution[a, b]] over domain $0<x<2$, where for $n=2, \alpha=1, \beta=1/2$, $n=3, \alpha=2, \beta=1$, $n=4, \alpha=3, \beta=3/2$, ... However, these are approximate - I'm not sure if those parameters are exact. Nov 21, 2020 at 16:17
• @PenelopeBenenati for $n>3$ to create a point on the $n$-sphere, just generate $n$ variates from a normal distribution, and Normalize the point: points = Partition[ Normalize /@ RandomVariate[NormalDistribution[0, 1], {1000000, n}], 2]; distances = EuclideanDistance @@@ points; Nov 21, 2020 at 16:19
• Thank you @flinty . Nov 21, 2020 at 16:19

RandomPoint can be used to sample from arbitrary Region definitions, and Sphere describes the unit sphere in $$R^n$$ (it can be both a geometric region and a graphics primitive). So @flinty's original solution can be generalized to arbitrary numbers of dimensions. Using this type of approach avoids having to know very much about the problem (as in @Roman's solution).

A Monte Carlo based way to estimate this can look like the following:

(*sample Euclidean distances of pairs of points*)
distanceDistributionOnSphere[dimensionality_, nSamples_ : 10^5] :=
With[
{randomPointsOnSurfaceOfNSphere = RandomPoint[Sphere[dimensionality], {nSamples, 2}]},
EuclideanDistance @@@ randomPointsOnSurfaceOfNSphere]

(*Evaluate mean of the sample*)
MeanAround /@ distanceDistributionOnSphere /@ Range[10]


This yields the same results as noted above.

It may also make it easier to generate a region that combines a Sphere with a HalfPlane in order to implement the subsequent question in the comment thread about hemisphere distances, but I don't know enough about $$n>3$$ dimensional geometry to implement that correctly....