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

Question

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.

Answer 1

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

Answer 2

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

Answer 3

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....