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.