Let me explain it differently than the rest, with a simple, intuiitive derivation.

Again, the problem is perfectly symmetric in polar angle, so that part is trivial. All we need to understand is the distribution in the radius, r.

Assume you wish to see a uniform sampling on the annulus from a to b in radius. Now consider any fixed radius r, so now just a circle of radius r. We want to see proportionally 2*pi*r*d points around that infinitessimally narrow circular annulus. Essentially, this gives us the radial distribution. The PDF is simply:

Here d is the packing density, and r will vary from a to b. We need to normalize the pdf so it has unit area on that domain.

PDF = PDF/int(PDF,r,[a,b])

PDF =

I hope that makes sense, in a very intuitive way. Now that we have the PDF, we can compute the CDF of the radial distribution, as simply the integral over r.

CDF = int(PDF,r,[a,R])

CDF =

Again, r will live on the interval [a,b]. And also see that I used the definite integral from a to the radius R. This is important!

For example, with a=1, b=10, we would have:

fplot(subs(PDF,[a,b],[1,10]),[1,10])

fplot(subs(CDF,[a,b],[1,10]),[1,10])

As you can see, the PDF is linear, and the CDF is a simple quadratic polynomial.

How do we sample from a distribution when the CDF is known, and it has a simple inverse? This itself is simple enough. We use the inverse of the CDF. That is, we choose a UNIFORM random sampling on the interval [0,1]. rand does this perfectly. Then we solve for R.

CDFinv = solve(CDF == u,R)

CDFinv =

We clearly want the second solution there, since it is the positive one.

CDFinvfun = @(u,a,b) sqrt((b^2-a^2)*u + a^2)

CDFinvfun =

@(u,a,b)sqrt((b^2-a^2)*u+a^2)

I just wrote it as a function, since matlabFunction did something screwy to me, and it is too early in the AM right now for me to think.

Rsample = CDFinvfun(Usample,1,10);

histogram(Rsample,100,'norm','pdf')

And you see a perfectly linear sample PDF drop out ion the interval [a,b]. Now lets try it out.

Rsample = CDFinvfun(Usample,1,10);

And that is exactly what we want to see. A nice uniform sampling on the annulus.

If what I did does not make sense, then go through it step by step. I tried to be as complete as possible in my Answer, but I am sure I went too fast at some point.

Addendum: Some additional points of interest I may have glossed over above,

- The entire solution stems directly from the assumption of a uniform scattering on the annulus. That directly implies the PDF, and then the rest follows.
- The use of a definite integral to compute the CDF is important, since if you just compute an indefinite integral with no other flags to int, the constant of integration will be lost. The definite integral implicitly builds the (correct) constant of integration into the CDF. I think this is an easily missed point.
- The above analysis can also trivially be extended to 3 or more dimensions, so with annular spheres. There the transformation will be different, but as easily derived based on the resulting polynomial PDF on the finite annular domain.
- If you follow the solution I've posed, you can see that whan a==0, and we want to simply sample uniformly over the entire circle of radius b, then the intuitive square root transformation that Adam proposed now appears. But that clearly fails when a>0. When a==0, the inverse CDF reduces to simply b*sqrt(u),