Assume I need to generate (pseudo)random points uniformly distributed on the surface of a sphere with radius $1$ and center at coordinate's origin.
The obvious, but incorrect, approach would be to use spherical coordinate system: pick up values of $\phi\in]0, 2\pi]$ and $\theta\in]0,\pi]$ from uniform distribution, this would give us a random point $(\phi, \theta)$ on a sphere.
The not-so-obvious, but correct way of doing this (thanks to comments) as explained here: https://mathworld.wolfram.com/SpherePointPicking.html :
choose $u$ and $v$ from uniform distribution $(0, 1)$. Then
$$ \phi=2\pi u $$ $$ \theta=cos^{-1}(2v - 1) $$
gives me a correct distribution.
Now, what if instead I take points $\boldsymbol{p}$ with coordinates $x$,$y$ and $z$ from uniform distribution on $[-1, 1]$ interval (call it $U(-1, -1)$) and then normalize the results (neglect the case when all three coordinates are zero)?
That is:
$$ \boldsymbol{p}_s=\frac{\boldsymbol p}{\parallel \boldsymbol p \parallel} $$
Will it still give uniform distribution on the sphere?
What if $x$, $y$ and $z$ are taken from the Normal distribution with mean $0$ and standard deviation $1$ : $N(0, 1)$ ?
Is there any way to check the 'quality' of the distribution?
Suppose $a$ is a random point from the distribution with probability density function $p : S^2 \to \mathbb{R}$, and $X$ is any measurable subset of the sphere $S^2$. Then we have
$$\int_{b \in S^2} p(b)\, d\sigma = 1$$
$$ P(a \in X) = \int_{b \in X} p(b)\, d\sigma $$
"Uniform distribution" means that the probability the point is a member of any measurable subset of the domain is proportional to the measure of that set. For the sphere $S^2$, the standard measure is the area (the surface integral of the constant function $1$). So if $p$ has this property, then for every measurable $X \subseteq S^2$,
$$ \frac{P(a \in X)}{P(a \in S^2)} = \frac{\int_X d\sigma}{\int_{S^2} d\sigma} $$
$P(a \in S^2)=1$, and $\int_{S^2} d\sigma$ is the surface area of the whole sphere, which is $4 \pi$. With the formula for $P(a \in X)$ above, this gives
$$ P(a \in X) = \int_{b \in X} p(b)\, d\sigma = \frac{1}{4\pi} \int_X d\sigma $$
The only way this can be true for a continuous $p$ and for all subsets $X \subseteq S^2$ is if $p(b) = \frac{1}{4 \pi}$ everywhere.
In spherical coordinates with $0 \leq \phi < 2\pi$, $0 \leq \theta \leq \pi$, and $\rho=1$, we have
$$ P(\phi < \phi_0 \land \theta < \theta_0) = \frac{1}{4\pi} \int_0^{\phi_0} \int_0^{\theta_0} \sin \theta\, d\theta\, d\phi = \frac{1}{4\pi} (1-\cos \theta) \phi $$
We can satisfy this by
$$ \begin{align*} P(\theta < \theta_0) &= \frac{1-\cos \theta_0}{2} \\ P(\phi < \phi_0) &= \frac{\phi_0}{2\pi} \end{align*} $$
If we have independent uniformly random variables $u$ and $v$, we can transform them to these distributions by inverting those equations: $$ \begin{align*} u = \frac{\phi}{2\pi} &\implies \phi = 2\pi u \\ v = \frac{1-\cos \theta}{2} &\implies \theta = \cos^{-1}(1-2v) \end{align*} $$
The given correct transformation in the question is also fine, since it differs just by $\phi \mapsto 2\pi - \phi$.
The proposed distribution with spherical coordinates $\phi$ and $\theta$ independent and uniformly random gives
$$ P(\phi < \phi_0 \land \theta < \theta_0) = \frac{1}{2 \pi^2} \phi_0 \theta_0 $$
This does not match the $P(\phi < \phi_0 \land \theta < \theta_0) = \frac{1}{4\pi} (1-\cos \theta) \phi$ for the distribution uniform on the sphere.
Another proposed distribution takes $x,y,z$ from independent and uniformly random distributions, then takes the unit vector in the direction of $(x,y,z)$. If $\theta_0 < \frac{\pi}{4}$, then the probability $P(\theta < \theta_0)$ is the volume of the cone whose apex is the origin and whose base is the disk $x^2+y^2 \leq \sin^2 \theta, z=1$, divided by the total volume $8$ of the distribution's cube $(-1,1) \times (-1,1) \times (-1,1)$. This volume is $\frac{\pi}{3} \sin^2 \theta$. This does not match the $P(\theta < \theta_0) = (1-\cos \theta_0)/2$ for the distribution uniform on the sphere.
Finally, a proposed distribution takes $x,y,z$ from independent Gaussian distributions. So the probability density function for the original $(x,y,z)$ is
$$ p_G(x,y,z) = \left(\frac{1}{\sigma\sqrt{\pi}}\right)^3 e^{-x^2/\sigma^2} e^{-y^2/\sigma^2} e^{-z^2/\sigma^2} = \frac{1}{\sigma^3 \pi^{3/2}} e^{-\rho^2/\sigma^2} $$
Then given a measurable set $X \subseteq S^2$,
$$ \begin{align*} P(a \in X) &= \int_{\{(x,y,z) \mid (x/\rho, y/\rho, z/\rho) \in X\}} p_G(x,y,z)\, dx\, dy\, dz \\ P(a \in X) &= \int_{\{(\phi,\theta,\rho) \mid (\cos \phi \sin \theta, \sin \phi \sin \theta, \cos \theta) \in X\}} p_G(x,y,z)\, \rho^2 \sin \theta\, d\rho\, d\theta\, d\phi \\ P(a \in X) &= \frac{1}{\sigma^3 \pi^{3/2}} \int_{X} d\sigma \int_0^\infty e^{-\rho^2/\sigma^2} \rho^2 d\rho \\ P(a \in X) &= \frac{1}{4\pi} \int_{X} d\sigma \end{align*} $$
This method does give the correct distribution. (It might take a little more CPU time than the other correct methods if they're all converted to computer algorithms.)
A note: You used the $\theta$ and $\phi$ coordinate variables in the swapped meaning from the way I usually see and use them. I've followed your usage in this answer, but make sure you're using $x = \sin \theta \cos \phi, y = \rho \sin \theta \sin \phi, z = \cos \theta$.