How to pick random uniformly distributed points in a sphere has been asked before. The difference is that I don't want uniform distribution, rather I would like the number density to scale by $r^{-1}$.
The only method I can think of is to generate random x, y and z coordinates, reject those that are outside the sphere and then do another test. For simplicity I chose to distribute points between radius 1 and the sphere's radius R, with these values I calculated a probability function $p(r)=\frac{1}{\rm{ln}(R)r}$. The test was that I generated another random number and rejected the point if that number was larger than p(r). However I don't know if this is correct, and I'm hoping for a better idea since it was very slow for large R.
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.