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$.
Here's an asymptotic bound - hopefully it's tight.
We throw $N$ little balls of diameter $D$ randomly (uniformly) inside the big sphere of radius $R$.
UPDATE: The new version (lower half) is better than what follows.
Neglecting border effects (reasonable if $N$ is large) the probability that the ball $i$ is "free" (no other overlaps with it) is
$$P(F_i)=\left(1-v(D)\right)^{N-1} \tag{1}$$
where $v(D) \triangleq D^3/R^3$
The probability that all balls are free can be bounded as
$$P(\cap F_i)=1- P(\cup F_i^c)\ge 1 - N(1-P(F_i)) \triangleq g(D,N) \tag{2} $$
For large $N$
$$g(D,N) \approx 1 - N^2 \, v(D) \tag{3} $$
in the range where this is positive, ie. $0\le D \le D_1 \triangleq R/N^{2/3} $
Now, let $t$ be the minimun distance between the sphere centers. Then
$$P(t \ge D) = P(\cap F_i) \ge g(D,N) \tag{4}$$
And then
$$E(t) = \int_0^{\infty} P(t \ge D) dD \ge\int_0^{D_1} g(D,N) \, dD \approx D_1- N^2 \frac{ D_1^4}{4 R^3} = \frac{3 }{ 4 } \frac{R}{N^{2/3}} $$
(Simulation data suggests that the order is right, and so is the bound, but the real coefficient is around $1.12$ - perhaps $9/8$)
Update: (Improved version)
A better approach can be obtained by considering instead of $F_i$ (free ball) the event $S_j\equiv$ "separated pair" (pair of balls are separated, they do not overlap) where $j$ indexes the $M=N(N-1)/2 \approx N^2/2$ pairs.
By the same reasoning:
$$P(S_j)=1-v(D) =1 - \frac{D^3}{R^3} \tag{5}$$
$$P(\cap S_i) \ge \max(1 - M(1-P(S_i)),0)= \max\left(1 - M \frac{D^3}{R^3},0\right) \triangleq h(D,M) \tag{6} $$
The range where $h(D,M)$ is positive, ie. $0\le D \le D_2 \triangleq R/M^{1/3} $
Now, let $t$ be the minimun distance between the sphere centers. Then
$$P(t \ge D) = P(\cap S_i) \ge h(D,M) \tag{7}$$
And then
$$E(t) = \int_0^{\infty} P(t \ge D) dD \ge\int_0^{D_2} h(D,M) \, dD =\\
= \frac{3}{4} \frac{R}{M^{1/3}}
\approx 0.945 \frac{R}{\sqrt[3]{N(N-1)}} \approx 0.945 \frac{R}{N^{2/3}} \tag{8}$$
Update 2 : A simple heuristic which seems to produce the correct coefficient:
Following the approach above, we could assume that $S_i$ are asympotically independent, and then:
$$P(\cap S_i) \approx \left(1-\frac{D^3}{R^3}\right)^M \tag{9}$$
Then
$$E(t) \approx \int_0^{R}\left(1-\frac{D^3}{R^3}\right)^M dD =\\= R \, \Gamma(4/3) \frac{\Gamma(M+1)}{\Gamma(M+4/3)} \approx R \, \Gamma(4/3) M^{-1/3} \approx 1.12508368 \frac{R}{N^{2/3}} \tag{10}$$
Update 3 : Regarding corrections for border effects.
(Lets assume $R=1$ to save notation, it's just a scale factor)
If we wished to include border effects we should replace $(5)$ (computing the balls intersection as here) by
$$1-D^3+\frac{9}{16}D^4 -\frac{1}{32}D^6 \hspace{1cm} 0\le D\le 2$$
The integral gets more complicated, but the (first order) asymptotic result is not altered:
Lemma: For any positive differentiable function $g(x)$ which , in $[0,+\infty)$, has global maximum at $g(0)=1$, and which has zero first and second derivates $g(x)=1-a x^3 + O(x^4)$ we have (variation of Laplace method, see eg here sec 2.1.3)
$$ \int_0^\infty g(x)^M dx = \frac{\Gamma(1/3)}{3 a^{1/3}} M^{-1/3}+ o(M^{-1/3})$$
which again leads as to $(10)$.
Best Answer
It is tempting to generate points (or equivalently, directions) uniformly distributed on the sphere in this way. Unfortunately it is not correct. For a short discussion of the issue, you may want to look at the following.
Geometric probability has many examples of this type, where different methods of picking a direction, or a point, "at random" yield quite different answers. And you should not "feel dumb." There is a fair history of different answers being obtained by using different, but not obviously different, notions of randomly chosen points.
The issue in this case is that we get a heavy biasing towards the two poles. A strip around the equator, of large area, has a much smaller probability of being selected than a region of similar area near the poles. This can be seen, among other ways, by looking at the Jacobian. But at the informal level, it can also be seen by examining the piece from $\phi=0$ to $\phi=0.1$, and then the piece from $\phi=\pi/2 -0.1$ to $\phi=\pi/2$. The disparity of areas is clear. Thus your integral does not give uniform weight to all directions, and therefore does not "average" $f$ in the intended way.