6

I have a sphere of radius $R_1$, and a smaller, concentric sphere of radius $R_2$. Let them be centered at the origin $(0,0,0)$. I need to generate random points with uniform density in the volume between the two sphere surfaces.

Here I found a solution for picking random points in a sphere volume, using

$$\frac{R_s U^{1/3}}{\sqrt{X_1^2 + X_2^2 + X_3^2}} (X_1, X_2, X_3)$$

where $X_1,X_2,X_3$ are independent normal random variables with mean 0 and variance 1, and $U$ is uniformly distributed between 0 and 1.

That strategy could be used for my problem as well, combining it with a rejection method:

  1. Generate a random point $(x,y,z)$ within $R_1$, using the abovementioned solution
  2. If $x^2 + y^2 + z^2 \leq R_2^2$, go back to step 1

But this is quite inefficient when the spherical shell is small. I guess it is possible to use a transformation method to generate points between $R_1$ and $R_2$. Can this be done?

Community
  • 1
Geert Eikelboom
  • 65
  • You can use generate a random spherical coordinate $(\rho, \theta, \phi$, $R_2 \leq \rho \leq R_1$, where each of the variables is weighted according to their appearance in the volume form $\rho^2 \sin \phi , d\rho , d\theta , d\phi$. – Travis Willse Jan 20 '15 at 12:29
  • @Travis That is the correct idea. – Felix Marin Jan 21 '15 at 11:05

2 Answers2

8

No need to use accept-reject methods here.

First generate a point on the unit sphere ($r=1$). There are several ways, for example:

$$\theta \sim U[0,2\pi]$$ $$Z_0 \sim U[-1,1]$$ $$X_0=\sqrt{1-Z_0^2} \cos(\theta)$$ $$Y_0=\sqrt{1-Z_0^2} \sin(\theta)$$

Then generate a random radius $R$ in $[r_1,r_2]$ with density proportional to $r^2$ : first generate $T \sim U[r_1^3,r_2^3]$ and then

$$R= T^{1/3}$$

(so that $f_R(r)=\frac{3}{r_2-r_1}{r^2}$ as desired)

Finally scale: $X=X_0 R$, $Y=Y_0 R$, $Z=Z_0 R$

It's also possible to generate the random point in the unit sphere volume, like in the question body, and then traslate+rescale radially.

Community
  • 1
leonbloy
  • 63,430
  • @Rahul Of course, fixed - thanks – leonbloy Jan 21 '15 at 13:29
  • 1
    In practice, using current technology, I think it's faster to compute an extra value or two via a Mersenne Twister (due to rejection) than to compute a sine, cosine, or square root. That's one reason to use rejection to generate a normal pseudorandom variable although a non-rejection method is available. Another reason is that on a computer, the alleged "real" uniform variable is really only a discrete variable with a large number of possible outcomes. In practice I would probably still use my answer to this question. Mathematically, however, I think I prefer this answer. – David K Jan 21 '15 at 13:31
  • @DavidK I would be somewhat surprised by that - the trig functions and square root have all been highly optimized too, and the compares and branching necessary to handle the rejection algorithm could easily swamp the trig computations. – Steven Stadnicki Dec 22 '15 at 20:52
  • @StevenStadnicki I periodically check the performance of functions in my code and find that trig functions and square root are suprisingly expensive; on the other hand, modern computer architectures do suprisingly well with branches. Looking at both algorithms again, I'm making somewhere between two and three expensive function calls per successful outcome, the solution above makes four; so it really comes down to how badly the "failed" branches affect processing. – David K Dec 22 '15 at 21:18
  • @StevenStadnicki FWIW, I did vote for this answer, because I really do like it. – David K Dec 22 '15 at 21:33
1

A very good way to generate normal variables is to start with pairs of uniform variables and apply a rejection method. If you're considering these normal variables as input to a rejection method, then, why not start with uniform variables so you only have to do one level of rejection?

You could generate a random point uniformly distributed within the volume of the sphere of radius $R_1$ (using three uniform variables and a rejection method) and then move that point further from the origin by applying a transformation of the form $r' = f(r)$ ($r$ the original distance from the origin, $r'$ the distance after moving the point) so that the result is uniformly distributed in the spherical shell.

The transformation function $f$ should satisfy $$\int_r^{R_1} t^2 dt = k \int_{f(r)}^{R_1} t^2 dt$$ where $k = \dfrac{R_1^3}{R_1^3-R_2^3}$ is the ratio of the volume of the entire sphere to the volume of the spherical shell. I get the relationship $$(f(r))^3 = \frac 1k r^3 + \frac{k-1}{k} R_1^3 = \frac{R_1^3-R_2^3}{R_1^3} r^3 + R_2^3 $$ by computing both integrals and rearranging terms. Note that $f(0) = R_2$ and $f(R_1) = R_1$, as desired. Another way of writing this is $$\frac{f(r)}{r} = \left(\frac{R_1^3-R_2^3}{R_1^3} + \frac{R_2^3}{r^3} \right)^{\frac13}.$$

So the procedure is to generate three uniform variables, $U_1$, $U_2$, and $U_3$ (between $-R_1$ and $R_1$), and compute $r = \sqrt{U_1^2 + U_2^2 + U_3^2}.$ If $r \leq R_1$, the desired point within the spherical shell is $$(x,y,z) = \left(\frac{R_1^3-R_2^3}{R_1^3} + \frac{R_2^3}{r^3} \right)^{\frac13} (U_1, U_2, U_3).$$ If $r > R_1$ you reject the three variables and try again. You will be able to use $\frac\pi6$ of the triples you generate, which is worse than the $\frac\pi4$ of pairs of uniform variables that you would be able to use when generating pairs of normal variables by the usual rejection method, but a lot better than rejecting everything inside the shell if the shell is thin.

Geert Eikelboom
  • 65
David K
  • 98,388