Create random number within an annulus

问题:

I am trying generate a random number that is within an annulus, i.e. we have a max and min radius. I tried doing:

while True:
    x=random.uniform(-maxR, maxR)
    y=random.uniform(-maxR, maxR)
    R=math.sqrt(x**2 + y**2)
    if R <= maxRadius and R >= minRadius:
        if x>= -maxRadius and x <= maxRadius and x<=-minRadius and x>= minRadius:
            print "passed x"
            if y>= -maxRadius and y <= maxRadius and y<=-minRadius and y>= minRadius: 
                break

But this is very slow. Is it possible to feed more contraints into random.uniform or is there another method?

回答1:

In general you can either draw the correct distribution directly or use rejection.

To draw directly use

• draw theta uniformly on [0,2pi): theta = random.uniform(0,2*pi)

• draw r from the power-law distribution r^1.

  The only complexity compared to doing this for a circle is that you PDF runs from [r_min,r_max] not [0,r_max]. This leads to

  CDF = A \int_{r_min}^{r} r' dr' = A (r^2 - r_min^2)/2

  for A the normalizing constant

  A = 2/(r_max*r_max - r_min*r_min)
  

  implying that

  r = sqrt(2*random.uniform(0,1)/A + r_min*r_min)
  

  and you can simplify slightly.

• then compute (x,y) by the usual transformation from radial coordinates
  x = r * cos(theta)
y = r * sin(theta)

This method of integrating the PDF, normalizing the CDF and inverting is general and is sometimes called the "Fundamental Theorem of Sampling".

Rejection

Draw (x,y) on a box big enough to contain the annulus, then reject all cases where `r = sqrt(xx + yy) exceeds r_max or is less than r_min.

This is reasonably efficient if the hole in the middle is small, and very inefficient if the hole is large.

回答2:

The method you're using should work pretty efficiently for a thick annulus (where r1 <<< r2). If you're dealing instead with a narrow annulus (r2 - r1 <<< r1), then you may want to use something like this instead:

r = random.uniform(r1, r2)
theta = random.uniform(0, 2 * PI)
x = r * math.sin(theta)
y = r * math.cos(theta)

Note that this gives mildly non-uniform results (there's a constant distribution of points per unit angle, not per unit area).