10

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?

madtowneast
  • 2,350
  • 3
  • 22
  • 31
  • What are the last two if-statements for? Why is simply checking if `R` (the radius) is within range not sufficient? – Rob Wouters Jan 28 '12 at 20:04
  • You _could_ always generate the numbers beforehand and store them in a dict mapping. – Joel Cornett Jan 28 '12 at 20:25
  • There are already endless questions about drawing from non-unifrom PDFs, and from circles in particular. The better answers to those question contain the method for drawing from arbitrary invertable PDFS. – dmckee --- ex-moderator kitten Jan 28 '12 at 20:26
  • Obviously the angle has to be uniformly distributed between 0 and 2π. As for the radius distribution, you should use the [Inverse Transform Method](https://en.wikipedia.org/wiki/Inverse_transform_sampling). – jpmarinier Feb 20 '23 at 19:59

1 Answers1

24

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.

Community
  • 1
  • 1
dmckee --- ex-moderator kitten
  • 98,632
  • 24
  • 142
  • 234
  • Thanks for the info. Hadnt thought of that at all. – madtowneast Jan 28 '12 at 21:52
  • I have drawn this solution in Mathematica: http://i.imgur.com/7L04w.png It's looking nice and uniform! :D – Ricardo Sanchez-Saez Oct 25 '12 at 10:36
  • 1
    "and you can simplify slightly." Am I correct in thinking that in the general case (assuming your language of choice supports it) it can be simplified to this: `r = sqrt(random.uniform(r_min*r_min,r_max*r_max))` ? – dgnuff Jan 21 '16 at 23:52