To generate a uniform distribution on a sphere $S^n$ in $\mathbb R^{n+1}$, we can normalize a vector whose entries are $n+1$ i.i.d normal random variables. If $\rho$ is a correlation, $|\rho|<1$, how can we generate a uniform distribution on the manifold
$$
M = \left \{ x, y \in S^n: x^Ty = \rho \right\}\ ?
$$



Uniform distribution on a manifold
It can be done in pretty much the same way as for a single vector by using the fact that if you fix $x$, then the conditional distribution of $y-x$ is uniform on the sphere of radius $\sqrt{1-\gamma^2}$ in the hyperplane perpendicular to $x$. Therefore, first you generate $x$ (as you say, by normalizing a vector with i.i.d. standard normal coordinates), then you complement $x$ to an orthonormal basis $(x,e\_2,\dots, e\_n)$, generate, as above, a vector $y'$ uniformly distributed on the unit sphere of the space spanned by $e\_2,\dots, e\_n$, and finally put $y=\gamma x + \sqrt{1-\gamma^2} y'$.


PS While I was typing you replaced $\gamma$ with $\rho$ and $n$ with $n+1$.


