motion on the surface of a 3-sphere I am trying to use motion on the surface of a 3-sphere to represent constrained, autocorrelated points in 4D.  The phenomena I'm trying to model is with the autocorrelation is fairly cyclic -- if a particular dimension is increasing, it tends to increase for awhile then decrease and vice versa.

Mechanically, this works fine by choosing an initial $(\phi_1,\phi_2,\phi_3)$, angular velocities and then doing a series of steps $\phi_1^{(n+1)}=\phi_1^{(n)}+v_1^{(n)}$, etc.  I also do multi-normal noising to the angular velocity on each step.  I simulate motion on the hypersphere surface for many steps, then translate the angles back into cartesian coordinates.

This generates values that comply with constraint to the surface, but when I check the distribution of points in the 4D space, I find that ${x_1}^2$ average is roughly 1/2, ${x_2}^2$ is roughly 1/4, and ${x_3}^2$ and ${x_4}^2$ are both roughly 1/8.

I need the averages to be roughly equal.  Looking at the equations for converting between coordinate systems, this outcome seems obvious in retrospect -- each cosine and sine of $\phi_i$ terms would have the same average value, and each subsequent $x_j$ multiplies by another cosine or sine term, and movement in the angular dimensions is independent.

Is there a straightforward way to modify this general approach to yield the correct average $x$ behavior?
Thanks to @IvanNeretin's comments, I got the right perspective.

My dynamic model results in uniform distributions for all of the $\phi_i$ terms, which is what I thought I wanted.

But when I thought in terms of randomly distributing points on the surface, I realized that uniformly distributed points only result in uniformly distributed $\phi_{n}$ for an $n$-sphere.  Thinking about the regular sphere (a 2-sphere), the azimuthal angle ($\phi_2$) is uniformly distributed, but the polar angle is not - this is easy to realize by recognizing the relative area strips associated with projecting $d\phi_1$ at different $\phi_1$.  Near the poles, the strip area is small, while near the equator, it is large.  The relative difference is a factor of $\sin\phi_1$.  So a uniformly distributed sample of $\phi_1$ on $(0,\pi)$ needs to be re-distributed to be a $\sin\phi_1$ weighted sample, which can be accomplished with the CDF method.  We know the associated CDF would be

$$
CDF(x) = \frac{\int_0^{x}\sin\phi d\phi}{\int_0^\pi \sin\phi d\phi}
$$

which can be solved, inverted, and applied to uniformly sampled $\phi_1/\pi$

For each additional sphere dimension, there is an extra power of $\sin\phi$ in the leading $\phi_1$.  So for the the 3-sphere, the distribution of $\phi_1$ needs to become $\sin^2\phi_1$, $\phi_2$ becomes $\sin\phi_2$ and $\phi_3$ remains uniform.

By squishing the $\phi_i$ dimensions to $\sin^{n-i}\phi_i$ distributions, I am able to keep the boundary constraint, get the autocorrelation, and meet the mean balance of the $x_i$s.
Thanks to @IvanNeretin's comments, I got the right perspective.

My dynamic model results in uniform distributions for all of the $\phi_i$ terms, which is what I thought I wanted.

But when I thought in terms of randomly distributing points on the surface, I realized that uniformly distributed points only result in uniformly distributed $\phi_{n}$ for an $n$-sphere.  Thinking about the regular sphere (a 2-sphere), the azimuthal angle ($\phi_2$) is uniformly distributed, but the polar angle is not - this is easy to realize by recognizing the relative area strips associated with projecting $d\phi_1$ at different $\phi_1$.  Near the poles, the strip area is small, while near the equator, it is large.  The relative difference is a factor of $\sin\phi_1$.  So a uniformly distributed sample of $\phi_1$ on $(0,\pi)$ needs to be re-distributed to be a $\sin\phi_1$ weighted sample, which can be accomplished with the CDF method.  We know the associated CDF would be

$$
CDF(x) = \frac{\int_0^{x}\sin\phi d\phi}{\int_0^\pi \sin\phi d\phi}
$$

which can be solved, inverted, and applied to uniformly sampled $\phi_1/\pi$

For each additional sphere dimension, there is an extra power of $\sin\phi$ in the leading $\phi_1$.  So for the the 3-sphere, the distribution of $\phi_1$ needs to become $\sin^2\phi_1$, $\phi_2$ becomes $\sin\phi_2$ and $\phi_3$ remains uniform.

By squishing the $\phi_i$ dimensions to $\sin^{n-i}\phi_i$ distributions, I am able to keep the boundary constraint, get the autocorrelation, and meet the mean balance of the $x_i$s.