This is in relation to the following blog post by John Baez: Random Points on a Sphere (Part 1)
It is claimed that, for two unit quaternions $x,y\in S^3$, uniformly randomly selected, the probability distribution of $|xy-yx|$ (how "far apart" $x$ and $y$ are from commuting) is given by: $$ P(d) = d-\frac{2d\arcsin(d/2)}{\pi},\quad d = |xy-yx|,\quad x,y\in S^3. $$ More generally, the probability distribution for the distance between two uniformly randomly chosen points on the $n$-sphere is given by: $$P_n(d) = \frac{d^{n-2}(4-d^2)^{\frac{n-3}{2}}\Gamma(n/2)}{2^{n-3}\sqrt{\pi}\ \Gamma(\frac{n-1}{2})},\quad d=|x-y|,\quad x,y\in S^n.$$ My knowledge of probability theory is in all honesty, severely lacking. So my question is:
Q: How does one derive these probability distributions?
To be clear, the formulas are probability densities, not cumulative distributions. The hypersphere formula is not a generalization of the quaternion formula, though similar math is involved in both. Also note the hypersphere formula you gave is for $S^{n-1}$ not $S^n$.
Let $F(d)=\mathrm{Prob}(|x-y|\le d\mid x,y\in S^{n-1})$ be the cumulative distribution for the distance between two points chosen uniformly at random from the unit sphere in $\mathbb{R}^n$. By symmetry, we can let $y$ be any fixed unit vector. The set of $x$ for which $|x-y|\le d$ is a hyperspherical cap enclosed in a hyperspherical sphere of spherical radius $\theta$ around $y$ defined by
$$d^2=(1-\cos\theta)^2+\sin^2\theta=2(1-\cos\theta). $$
[Try drawing this for $n=3$ as a visual aid for where the Pythagorean theorem comes in.] Then $F(d)$ will be the area of this cap divided by the total "area" of the sphere.
Let $A_{n-1}=2\pi^{n/2}/\Gamma(n/2)$ be the total area of the $(n-1)$-dimensional sphere $S^{n-1}\subset\mathbb{R}^n$. [This formula is derived in many places using recursion, induction, substitution.] Since the cap is a hypersurface of revolution we may take advantage of symmetry to calculate its "area." It can be sliced into cross-sections which are $(n-2)$-spheres of Euclidean radius $\sin\phi$ as the angle $\phi$ from $y$ varies from $0$ to $\theta$. By integrating the content of these slices with respect to $\phi$ we get the CDF
$$ F(d)=\frac{1}{A_{n-1}}\int_0^\theta A_{n-2}(\sin\phi)^{n-2}\,\mathrm{d}\phi. $$
Differentiating WRT $\theta$ then yields the PDF
$$ f(d)\frac{\mathrm{d}d}{\mathrm{d}\theta} = \frac{A_{n-2}}{A_{n-1}}(\sin\theta)^{n-2}. $$
To start simplifying this we need $d$ and $\mathrm{d}d/\mathrm{d}\theta$ in terms of $d$. Differentiate $d^2=\cdots$ to get
$$ 2d\frac{\mathrm{d}d}{\mathrm{d}\theta}=2\sin\theta \quad\Rightarrow\quad \frac{\mathrm{d}d}{\mathrm{d}\theta}=\frac{\sin\theta}{d}. $$
Also $\cos\theta=1-d^2/2$ so $\sin\theta=\sqrt{1-\cos^2\theta}=\sqrt{1-(1-d^2/2)^2}=d\sqrt{1-(d/2)^2}$. Now
$$ f(d)\frac{\sin\theta}{d}=\frac{A_{n-2}}{A_{n-1}}(\sin\theta)^{n-2} \implies $$
$$ f(d)=\frac{A_{n-2}}{A_{n-1}}d(\sin\theta)^{n-3}$$
Substituting in for the $A$s and $\sin\theta$ yields
$$ f(d)=\frac{2\pi^{(n-2)/2}/\Gamma(\frac{n-1}{2})}{2\pi^{n/2}/\Gamma(\frac{n}{2})} d\left(d\sqrt{1-\Big(\frac{d}{2}\Big)^2}\,\right)^{n-3} $$
which simplifies to the given formula
$$ f_{|x-y|}(d)=\frac{\Gamma(\frac{n}{2})d^{n-2}(4-d^2)^{(n-3)/2}}{2^{n-3}\sqrt{\pi}\Gamma(\frac{n-1}{2})} \quad (0\le d\le2). \tag{$\bullet$}$$
To compute the density for $|xy-yx|$ with $x,y\in S^3$ (for versors, aka unit quaternions, drawn uniformly at random WRT the usual Lebesgue measure), write $x,y$ in polar form as $x=\exp(\alpha\mathbf{u})$, $y=\exp(\beta\mathbf{v})$ and let $\angle\mathbf{uv}=\gamma$ be the angle between $\mathbf{x}$ and $\mathbf{y}$. Using convex angles only, these are uniquely defined for all $x,y\ne\pm1$ (measure zero, so no worries). Note $\alpha,\beta,\mathbf{u},\mathbf{v}$ are all independent random variables in their own right, with $\mathbf{u},\mathbf{v}\in S^2$ drawn uniformly at random. Then the commutator magnitude is given by
$$|xy-yx|=2|\mathrm{Im}(x)\times\mathrm{Im}(y)|=2|\mathrm{Im}(x)||\mathrm{Im}(y)||\mathbf{u}\times\mathbf{v}|=2\sin\alpha\sin\beta\sin\gamma. $$
The CDFs of $\sin\alpha,\sin\beta,\sin\gamma$ will all be calculated using spherical caps again, the first two in $S^3$ and the third in $S^2$. The set of all unit vectors whose convex angle $\psi$ made with a given unit vector satisfies $\sin\psi\le d$ is a pair of antipodal spherical caps with spherical radius $\delta$ and Euclidean radius $\sin\delta=d$. For a general sphere $S^{n-1}$ we have thus have the CDF
$$ F_{\sin\delta}(d) = \frac{2}{A_{n-1}}\int_0^\delta A_{n-2}(\sin\psi)^{n-2}\,\mathrm{d}\psi. $$
and differentiating WRT $\delta$ yields
$$ f_{\sin\delta}(d)\cos\delta = \frac{2A_{n-2}}{A_{n-1}}(\sin\delta)^{n-2} $$
which means the PDF is given by
$$ f_{\sin\delta}(d)=\frac{2A_{n-2}}{A_{n-1}}\frac{d^{n-2}}{\sqrt{1-d^2}} $$
Setting $n=4,4,3$ (for $S^3,S^3,S^2$ resp.) we get PDFs
$$ f_{\sin\alpha}(a)=\frac{4}{\pi}\frac{a^2}{\sqrt{1-a^2}}, \quad f_{\sin\beta}(b)=\frac{4}{\pi}\frac{b^2}{\sqrt{1-b^2}}, \quad f_{\sin\gamma}(c)=\frac{c}{\sqrt{1-c^2}}.~ $$
Convolving these three PDFs gives us the original CDF
$$ F_{|xy-yx|}(d)=\Big(\frac{4}{\pi}\Big)^2\!\iiint\limits_{2abc\le d}\frac{a^2}{\sqrt{1-a^2}}\frac{b^2}{\sqrt{1-b^2}}\frac{c}{\sqrt{1-c^2}}\,\mathrm{d}c\,\mathrm{d}b\,\mathrm{d}a. $$
For better parametrization, consider the complement
$$ 1-\Big(\frac{4}{\pi}\Big)^2\!\!\iiint\limits_{abc\ge d/2}\frac{a^2}{\sqrt{1-a^2}}\frac{b^2}{\sqrt{1-b^2}}\frac{c}{\sqrt{1-c^2}}\,\mathrm{d}c\,\mathrm{d}b\,\mathrm{d}a $$
which can now be turned into an iterated integral
$$1-\Big(\frac{4}{\pi}\Big)^2\!\int_{d/2}^1 \frac{a^2}{\sqrt{1-a^2}} \left(\int_{d/2a}^1 \frac{b^2}{\sqrt{1-b^2}} \left(\int_{d/2ab}^1 \frac{c}{\sqrt{1-c^2}}\,\mathrm{d}c\right)\mathrm{d}b\right)\mathrm{d}a. $$
It is possible to evaluate this directly, but for the density we can now differentiate WRT $d$ using the Leibniz integral rule. Despite this being a triple integral, all resulting terms vanish except where we apply the FTC to the innermost integral. The PDF $f_{|xy-yx|}(d)$ is thus given by
$$ -\Big(\frac{4}{\pi}\Big)^2\!\int_{d/2}^1 \frac{a^2}{\sqrt{1-a^2}} \Bigg(\int_{d/2a}^1 \frac{b^2}{\sqrt{1-b^2}} \frac{\frac{d}{2ab}(-\frac{1}{2ab})}{\sqrt{1-(\frac{d}{2ab})^2}}\,\mathrm{d}b\Bigg)\mathrm{d}a. $$
We can simplify this to
$$ \frac{4d}{\pi^2}\int_{d/2}^1 \frac{1}{\sqrt{1-a^2}} \Bigg(\int_{d/2a}^1 \frac{b}{\sqrt{(1-b^2)(b^2-(\frac{d}{2a})^2)}}\,\mathrm{d}b\Bigg)\mathrm{d}a $$
To evaluate the inner integral, set $k=\frac{d}{2a}$ and substitute $\ell=\sqrt{b^2-k^2}$ then $m=\ell/\sqrt{1-k^2}$:
$$ \int \frac{1}{\sqrt{1-b^2}} \frac{b\,\mathrm{d}b}{\sqrt{b^2-k^2}}=\int \frac{\mathrm{d}\sqrt{b^2-k^2}}{\sqrt{1-b^2}}=\int \frac{\mathrm{d}\ell}{\sqrt{(1-k^2)-\ell^2}} = $$
$$ \int\frac{\mathrm{d}\Big(\frac{\ell}{\sqrt{1-k^2}}\Big)}{\sqrt{1-\big(\frac{\ell}{\sqrt{1-k^2}}\big)^2}}=\int\frac{\mathrm{d}m}{\sqrt{1-m^2}}=\sin^{-1}(m)=\sin^{-1}\!\Big(\sqrt{\frac{b^2-k^2}{1-k^2}}\,\Big). $$
Therefore the inner integral evaluates simply to
$$ \int_k^1 \frac{b\,\mathrm{d}b}{\sqrt{(1-b^2)(b^2-k^2)}}=\Bigg[\sin^{-1}\!\Big(\sqrt{\frac{b^2-k^2}{1-k^2}}\,\Big)\Bigg]_k^1=\frac{\pi}{2}. $$
So now what we have is
$$ f(d)=\frac{4d}{\pi^2}\int_{d/2}^1\frac{\mathrm{d}a}{\sqrt{1-a^2}}\frac{\pi}{2} $$
which evaluates simply to
$$ f_{|xy-yx|}(d)=\frac{2d}{\pi}\cos^{-1}\Big(\frac{d}{2}\Big) \quad(0\le d\le2), \tag{$\large\star$}$$
which matches the given formula.