$\newcommand{\+}{^{\dagger}}
\newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle}
\newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace}
\newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack}
\newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,}
\newcommand{\dd}{{\rm d}}
\newcommand{\down}{\downarrow}
\newcommand{\ds}[1]{\displaystyle{#1}}
\newcommand{\expo}[1]{\,{\rm e}^{#1}\,}
\newcommand{\fermi}{\,{\rm f}}
\newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}
\newcommand{\half}{{1 \over 2}}
\newcommand{\ic}{{\rm i}}
\newcommand{\iff}{\Longleftrightarrow}
\newcommand{\imp}{\Longrightarrow}
\newcommand{\isdiv}{\,\left.\right\vert\,}
\newcommand{\ket}[1]{\left\vert #1\right\rangle}
\newcommand{\ol}[1]{\overline{#1}}
\newcommand{\pars}[1]{\left(\, #1 \,\right)}
\newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
\newcommand{\pp}{{\cal P}}
\newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,}
\newcommand{\sech}{\,{\rm sech}}
\newcommand{\sgn}{\,{\rm sgn}}
\newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}}
\newcommand{\ul}[1]{\underline{#1}}
\newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}
\newcommand{\wt}[1]{\widetilde{#1}}$
Lets ${\rm P}\pars{\vec{r}}$ the probability distribution of random points uniformly distributed in the surface of a sphere. $\ds{a}$ is the sphere radius. It satisfies ( note that $\ds{{\rm P}\pars{\vec{r}}}$ is, indeed, $\ds{\vec{r}}$-independent ):
$$
1 = \int_{S}{\rm P}\pars{\vec{r}}\,\dd S = {\rm P}\pars{\vec{r}}\int_{S}\dd S
={\rm P}\pars{\vec{r}}\pars{4\pi a^{2}}\quad\imp\quad{\rm P}\pars{\vec{r}}
={1 \over 4\pi a^{2}}
$$

Then,
$$
1 = \int_{S}{\dd S \over 4\pi a^{2}}
=\int_{0}^{\pi}\half\sin\pars{\theta}\,\dd\theta\int_{0}^{2\pi}{\dd\phi \over 2\pi}
$$
which we can visualize as two independent probability distributions $\ds{\half\,\sin\pars{\theta}}$ and $\ds{1 \over 2\pi}$.

Now, we generate random numbers $\ds{\xi_{\theta}}$ and $\ds{\xi_{\phi}}$
uniformly distributed in $\ds{\left[0,1\vphantom{\large A}\right)}$ such that

\begin{align}
\int_{0}^{\theta}\half\,\sin\pars{t}\,\dd t=\int_{0}^{\xi_{\theta}}\,\dd t\quad &\imp\quad\sin\pars{\theta \over 2} = \root{\xi_{\theta}}\quad\imp\quad\theta = 2\arcsin\pars{\root{\xi_{\theta}}}
\\[3mm]
\int_{0}^{\phi}{1 \over 2\pi}\,\dd t=\int_{0}^{\xi_{\phi}}\,\dd t\quad &\imp\quad
\phi = 2\pi\xi_{\phi}
\end{align}
Note that
$$
\sin\pars{\theta} = 2\sqrt{\xi_{\theta}\pars{1 - \xi_{\theta}}}
\qquad\mbox{and}\qquad\cos\pars{\theta} = 1 - 2\xi_{\theta}
$$

Point $\ds{\pars{x,y,z}}$ is given by
$$
\begin{array}{rcrcl}
x &= &a\sin\pars{\theta}\cos\pars{\phi}&=&2a\root{\xi_{\theta}\pars{1 - \xi_{\theta}}}
\cos\pars{2\pi\xi_{\phi}}
\\
y &=& a\sin\pars{\theta}\sin\pars{\phi}&=&2a\root{\xi_{\theta}\pars{1 - \xi_{\theta}}}
\sin\pars{2\pi\xi_{\phi}}
\\
z &=& a\cos\pars{\theta}&=&a\pars{1 - 2\xi_{\theta}}
\end{array}
$$
*Remember that $\ds{\xi_{\theta}\ \mbox{and}\ \xi_{\phi}}$ are random numbers uniformly distributed in $\ds{\left[0,1\right)}$.*

A typical 'run' of the $\tt C++$ script at the bottom yields:

( -0.08187 , -1.981 , -0.2642 )
( 1.562 , 0.7686 , -0.9843 )
( -0.2041 , 0.942 , 1.752 )
( -0.4777 , 1.891 , 0.4411 )
( 1.891 , 0.6153 , 0.2145 )

/* surfaceSphere_0.cc 18-jun-2014 Felix Marin
http://math.stackexchange.com/a/839189/85343
*/
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <iomanip>
#include <iostream>
#include <vector>
using namespace std;
const double RANDMAX1=double(RAND_MAX) + 1.0,TWOPI=2.0*M_PI;
void surfacePoint(double radius,vector<double> &p);
inline double rand01() // Return uniform values in [0,1)
{
return rand()/RANDMAX1;
}
int main()
{
vector<double> point;
point.resize(3U);
cout<<setprecision(4);
srand(size_t(time((time_t *)0)));
for ( unsigned char n=0 ; n<5 ; ++n ) {
surfacePoint(2.0,point);
cout<<"( "<<point[0]<<" , "<<point[1];
cout<<" , "<<point[2]<<" )\n";
}
return 0;
}
void surfacePoint(double radius,vector<double> &p)
{
static double temp,xiTheta,twoPiXiPhi;
xiTheta=rand01();
temp=2.0*radius*sqrt(xiTheta*(1.0 - xiTheta));
twoPiXiPhi=TWOPI*rand01();
p[0]=temp*cos(twoPiXiPhi);
p[1]=temp*sin(twoPiXiPhi);
p[2]=radius*(1.0 - 2.0*xiTheta);
}