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.

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

3 Answers 3


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}    *)

enter image description here


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    *)
  • $\begingroup$ Thank you @Roman . Do you know how to extend this solution for a unit $n$-hemisphere (instead of a unit $n$-sphere)? $\endgroup$ Nov 21, 2020 at 16:24
  • 1
    $\begingroup$ I think that's much harder because hyperspherical symmetry is broken. $\endgroup$
    – Roman
    Nov 21, 2020 at 16:33
  • 1
    $\begingroup$ Yes I think that it should be a new question, it's too different from this one. $\endgroup$
    – Roman
    Nov 21, 2020 at 17:18
  • 1
    $\begingroup$ Well done getting what appears to be the exact form. How did you manage it? $\endgroup$
    – flinty
    Nov 21, 2020 at 18:01
  • 2
    $\begingroup$ @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$. $\endgroup$
    – Roman
    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}]]

sphere point distance

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.

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

(* result: 4/Pi *)
  • $\begingroup$ 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. $\endgroup$ Nov 21, 2020 at 16:17
  • 1
    $\begingroup$ 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. $\endgroup$
    – flinty
    Nov 21, 2020 at 16:17
  • 1
    $\begingroup$ @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; $\endgroup$
    – flinty
    Nov 21, 2020 at 16:19
  • $\begingroup$ Thank you @flinty . $\endgroup$ 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] := 
  {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....


Not the answer you're looking for? Browse other questions tagged or ask your own question.