For the equivalence of 1 and 2, we can show that (left or right) translation of the unit quaternions by a unit quaternion is an isometry of $S^3$:
Given vectors $x,y\in S^3$, translating them both by $a\in S^3$ in the sense of quaternion multiplication leads to the dot product
$$
\begin{align}
(ax)\cdot(ay)
&=
\pmatrix{a_0x_0-\vec a\cdot\vec x\\a_0\vec x+x_0\vec a-\vec a\times\vec x}
\cdot
\pmatrix{a_0y_0-\vec a\cdot\vec y\\a_0\vec y+y_0\vec a-\vec a\times\vec y}
\\
&=a_0^2x_0y_0+\vec a^2x_0y_0+a_0^2\vec x\cdot\vec y+(\vec a\cdot\vec x)(\vec a\cdot\vec y)+
(\vec a\times\vec x)(\vec a\times\vec y)
\\
&=(a_0^2+\vec a^2)x_0y_0+(a_0^2+\vec a^2)\vec x\cdot\vec y
\\
&=x_0y_0+\vec x\cdot\vec y
\\
&=x\cdot y\;,
\end{align}
$$
so the translation is an isometry of $S^3$ and leaves the distribution 2 invariant. (I'm taking the more technical aspects of the Haar measure involving inner and outer regularity for granted in this context.)
If we split the rotations up as in 3, then for the distribution to be invariant under translation by arbitrary rotations, in particular the distribution of $\theta$ must be invariant under arbitrary rotations about $\hat w$, so $\theta$ must be uniformly distributed. Also the distribution of $\hat w$ clearly has to be rotationally symmetric about $\hat z$. Thus, it remains only to show that the uniform distribution of $\hat w\cdot\hat z$ in distribution 3 is correct.
I don't have an elegant argument without calculus for this, but here's a calculation that tries to do without too much brute force. For a rotation by $\alpha$ about $\hat u$, the rotated $z$ axis is
$$
\cos\alpha\,\hat z+\sin\alpha\,\hat u\times\hat z+(1-\cos\alpha)\hat u(\hat u\cdot\hat z)\;,
$$
and forming the dot product with $\hat z$ yields a relationship between the three cosines,
$$c=a+(1-a)b^2\;,$$
where $a=\cos\alpha$, $b=\hat u\cdot\hat z$ and $c$ is the cosine between the original and rotated $z$ axes.
To calculate the cumulative distribution function for $c$, note that $c\ge a$ and that $b$ is uniformly distributed according to distribution 2, so we can solve for the upper limit of $b$,
$$
b=\sqrt{\frac{c-a}{1-a}}\;,
$$
and find the density for $a$ according to distribution 2,
$$
\sin^2\frac\alpha2\mathrm d\alpha\propto\frac{1-\cos\alpha}{\sin\alpha}\mathrm d(\cos\alpha)=\frac{1-a}{\sqrt{1-a^2}}\mathrm da=\sqrt{\frac{1-a}{1+a}}\mathrm da\;,
$$
to find the cumulative distribution function for $c$ proportional to
$$
\int_{-1}^c\sqrt{\frac{c-a}{1-a}}\sqrt{\frac{1-a}{1+a}}\mathrm da=\int_{-1}^c\sqrt{\frac{c-a}{1+a}}\mathrm da\;.
$$
Now we can differentiate with respect to $c$ to get the density for $c$. The term from varying the integration limit vanishes since the integrand vanishes at $a=c$, so the result is proportional to
$$
\int_{-1}^c\frac1{\sqrt{(c-a)(1+a)}}\mathrm da=\left[2\arctan\sqrt{\frac{1+a}{c-a}}\right]_{-1}^c=\pi\;,
$$
a constant, which establishes the equivalence between distributions 2 and 3.