Standard Deviation of random unit vector

Question

I have a question i ve been thinking about for a few days:

I want to better understand the statistics of random 3 dimensional unit vectors. For example, i pick two angles $\theta$ and $\phi$ from $U(0,\pi)$ and $U(0,2\pi)$ respectively, i create the random vector $\mathbf {X} = \{\sin \theta \cos \phi, \sin \theta \sin \phi, \cos \theta\}$, repeat this $N$-times, and then i try to obtain the standard deviation over all $N$s of the components of $\mathbf{X}$. The mean is zero (the distribution is symmetric), so the standard deviation is:

$\sigma_i =\sqrt{\frac{1}{N-1} \sum^N X_i^2} \quad (1)$

for the component $i$. This is also how Matlab computes it. Now because of the spherical symmetry, the standard deviations ($N$ large) for all three components should be equal, i.e. $\sigma_1 = \sigma_2 = \sigma_3 = \sqrt{\frac{1}{3}}$, which i can prove with

$\sigma_i^2 = \frac{1}{4 \pi} \int_0^{\pi} \int_0^{2 \pi} X_i^2 \sin \theta d\theta d\phi \quad (2)$.

So equation (1) is wrong, because it misses the $\sin \theta$ weight. If this is implemented, then i obtain the correct result.

My question is now: If i don't know from which distribution i randomly pick the components, but still assume that it is symmetric around $0$ and that $\mathbf{X}$ is a unit vector, will $\sigma_i$ always be $\sqrt{\frac{1}{3}}$? And is there a way to obtain the weight just from equation (1) and the $X_i$'s?

Any help is appreciated !

EDIT: I will try to clarify: maybe i am mixing up the standard deviation (as used in the R program) with the second moment of the underlying distribution of the $X_i$'s. I want to find the latter. Here is my mathematica code:

n = 10000;
th = Table[Pi*RandomReal[], {n}]; ph = Table[2*Pi*RandomReal[], {n}];
x = Sin[th]*Cos[ph]; y = Sin[th]*Sin[ph]; z = Cos[th];
wx = WeightedData[x, Sin[ArcCos[z]]];
wy = WeightedData[y, Sin[ArcCos[z]]];
wz = WeightedData[z, Sin[ArcCos[z]]];
{StandardDeviation[x], StandardDeviation[wx]}
{StandardDeviation[y], StandardDeviation[wy]}
{StandardDeviation[z], StandardDeviation[wz]}

{0.499561, 0.578086} -> unweighted gives 1/2, weighted gives Sqrt[1/3]
{0.500725, 0.579126} -> unweighted gives 1/2, weighted gives Sqrt[1/3]
{0.706939, 0.574878} -> unweighted gives Sqrt[1/2], weighted gives Sqrt[1/3]

sqrt(1/3) is what i expect. My question now regards the case of arbitrary underlying distributions of the $X_i$'s

Answer 1

You are making everything more complicated than it needs to be, and not quite correct.

A standard formula for standard deviation will give the desired result if the distribution is uniform. There are various ways to get a uniform distribution over a sphere. Choosing $\theta$ uniformly between $0$ and $\pi$ is not part of any of the ways I've ever seen. Several correct methods are described in the answers to How to generate random points on a sphere.

It's not really meaningful to use $X_i$ in the integral, because there's nothing random about the integral. You simply integrate over all possible vectors: $$\frac{1}{4\pi} \int_0^\pi \!\!\int_0^{2\pi} x_i^2 \sin\theta\, d\theta\, d\phi.$$ It seems to me to be misguided to think of this integral as weighting $x_i$ by $\sin\theta$; I would rather say that $\iint \cdots\sin\theta\,d\theta\,d\phi$ represents the integration of something (whatever you put in place of "$\cdots$") over the surface of the sphere. In other words, $\sin\theta\,d\theta\,d\phi$ is an expression for $dS$ where $S$ measures surface area.

So what you have done is, first, to construct a non-uniform distribution of points on the surface of the sphere; then, because you have packed too many points near the poles, you have to apply weights to those points so that they have only the same relative influence that the correct (smaller) number of points in that region would have had. In doing this, you make it impossible to simply use formulas that were meant to apply directly to distributions of random data, because the thing you're trying to measure is not represented by the distribution you generated. Rather than fix the underlying distribution so that it works as it should, you try to work around it by hacking the formulas.

The hacking also obscures the question of what happens with "other" distributions that $X_i$ might have. If we don't really intend the distribution to be different, but instead mean to "hack" it until it gives results like a uniform distribution over the area of the sphere, then of course we will always get the answer $\frac1{\sqrt3}.$ On the other hand, if we actually mean to investigate some other distribution of $X_i$ according to how that distribution actually behaves, then we can easily get results that are not $\frac1{\sqrt3},$ just as you got results that were not $\frac1{\sqrt3}$ when you measured the actual standard deviations of the coordinates of your actual distribution of random vectors.

---

For the uniform distribution, the result $\frac1{\sqrt3}$ is a conseuqence of the well-known fact that when points are uniformly distributed over a three-dimensional sphere, each coordinate of the points is uniformly distributed over the interval $[-1,1].$ See Projection of a uniform distribution on a sphere.

The standard deviation of a uniform distribution on an interval $[a,b]$ is $\frac1{\sqrt{12}}(b - a),$ and taking $a = -1$ and $b = 1$ this evaluates to $\frac1{\sqrt3}.$

---

By the way, the formula $$\sigma_i =\sqrt{\frac{1}{N-1} \sum^N (X_i - \bar X_i)^2} $$ is an estimator for the standard deviation when the true mean of the distribution is unknown. When we do not know the true mean, we use the estimated mean instead, but this tends to give a smaller result than you would get by using the true mean. We compensate by using $\frac1{1-N}$ instead of $\frac1N.$ (This is not a perfect solution of that problem, but that's beyond the level of this discussion.)

In your case, since you know the mean is $0$ and do not need to use $\bar X_i$ to estimate the mean, the formula for standard deviation would be $$\sigma_i =\sqrt{\frac{1}{N} \sum^N X_i^2}. $$ Of course for large values of $N$ you may not notice the difference between the results of the formulas.