Let's say your sphere is centered at the origin $(0,0,0)$.
For the distance $D$ from the origin of your random pointpoint, note that you want $P(D \le r) = \left(\frac{r}{R_s}\right)^3$. Thus if $U$ is uniformly distributed between 0 and 1, taking $D = R_s U^{1/3}$ will do the trick.
For the direction, a useful fact is that if $X_1, X_2, X_3$ are independent normal random variables with mean 0 and variance 1, then
$$\frac{1}{\sqrt{X_1^2 + X_2^2 + X_3^2}} (X_1, X_2, X_3)$$
is uniformly distributed on (the surface of) the unit sphere. You can generate normal random variables from uniform ones in various ways; the Box-Muller algorithm is a nice simple approach.
So if you choose $U$ uniformly distributed between 0 and 1, and $X_1, X_2, X_3$ iid standard normal and independent of $U$, then
$$\frac{R_s U^{1/3}}{\sqrt{X_1^2 + X_2^2 + X_3^2}} (X_1, X_2, X_3)$$
would produce a uniformly distributed point inside the ball of radius $R_s$.
$\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{\ds}[1]{\displaystyle{#1}}%
\newcommand{\equalby}[1]{{#1 \atop {= \atop \vphantom{\huge A}}}}%
\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]{\,#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}$
Let's assume we have a mechanism to generate random numbers in $\left[0,1\right)$ and ${\rm P}\pars{\Omega_{\vec{r}}}$ is a distribution function for random points in a sphere of radius $a > 0$. $\Omega_{\vec{r}}$ is the solid angle. In this case, ${\rm P}\pars{\Omega_{\vec{r}}}$ is, indeed, $\Omega_{\vec{r}}\,\,$-independent:
$$
1 = \int{\rm P}\pars{\Omega_{\vec{r}}}\,\dd\Omega_{\vec{r}}
= {\rm P}\pars{\Omega_{\vec{r}}}\int\dd\Omega_{\vec{r}}
= {\rm P}\pars{\vec{r}}\pars{4\pi}
\quad\imp\quad{\rm P}\pars{\vec{r}} = {1 \over 4\pi}
$$
Then,
$$
1 = \int_{0}^{\pi}\half\,\sin\pars{\theta}\,\dd\theta\int_{0}^{2\pi}
\,{\dd\phi \over 2\pi}
$$
We can generate random numbers $\xi_{\theta}$ and $\xi_{\phi}$ such that:
$$
\bracks{~\half\,\sin\pars{\theta}\,\dd\theta = \dd\xi_{\theta}\,,
\quad\xi_{0} = 0 \imp \theta = 0~}\
\mbox{and}\
\bracks{~{\dd\phi \over 2\pi} = \dd\xi_{\phi}\,,\quad\xi_{0} = 0 \imp \phi = 0~}
$$
Those relations yield: $\ds{\half\bracks{-\cos\pars{\theta} + 1} = \xi_{\theta}}$
$\ds{\pars{~\mbox{or/and}\ \sin\pars{\theta/2} = \root{\xi_{\theta}}~}}$ and $\ds{\phi = 2\pi\,\xi_{\phi}}$:
$$\color{#0000ff}{\large%
\theta = 2\arcsin\pars{\root{\xi_{\theta}}}\,,\qquad \phi = 2\pi\xi_{\phi}}
$$
As we mentioned above, $\xi_{\theta}$ and $\xi_{\phi}$ are uniformly distributed in $\left[0, 1\right)$.
For a sphere of radio $a$ the random points are given by:
$$
\left\lbrace%
\begin{array}{rcl}
\color{#0000ff}{\large x} & = & a\sin\pars{\theta}\cos\pars{\phi}
= \color{#0000ff}{\large 2a\root{\xi_{\theta}\pars{1 - \xi_{\theta}}}\cos\pars{2\pi\xi_{\phi}}}
\\
\color{#0000ff}{\large y} & = & a\sin\pars{\theta}\sin\pars{\phi}
= \color{#0000ff}{\large 2a\root{\xi_{\theta}\pars{1 - \xi_{\theta}}}\sin\pars{2\pi\xi_{\phi}}}
\\
\color{#0000ff}{\large z} & = & a\cos\pars{\theta} = \color{#0000ff}{\large a\pars{1 - 2\xi_{\theta}}}
\end{array}\right.
$$
Best Answer
If you look at the accepted answer for the post you linked, you see that the direction is generated by sampling from a multivariate normal distribution, and then the radius is chosen according to the desired distribution. In your case, you don't need to define a minimal radius $R_0 > 0$ because when you multiply the scaling factor $1/r$ for probability of the radius by the surface area of the sphere $\Theta(r^2)$ you get $\Theta(r)$, and so you have $P(r \leq R) \propto \int_{0}^R (1/r)r^2 dr = \Theta(R^2)$ where you have a maximum radius $R_{max}$. This gives you your cdf (cumulative distribution function ) of $f(R) = P(r \leq R) = R^2/R_{max}^2$ for sampling the radius. The inverse of the cumulative distribution function is $g(X) = R_{max} \sqrt{X}$. If you sample a uniform random variable $X$ from $[0,1]$ and compute $r = g(X)$, this gives you the desired distribution on the radius $r$ that has cumulative distribution function $f(R)$. See http://en.wikipedia.org/wiki/Inverse_transform_sampling for an explanation on inverse transform sampling.