I have a one-dimensional function I(r)
defined as:
I(r) = A * exp( -2 * ( r / w ) ^ b )
where A
is arbitrary, w
and b
are variables that adjust the width and "rectangularity" of the function. I'm plotting this function with A = 1.0
, w = 1.0
, and different b
's for your reference:
Plots of I(r) with different b's
I want to be able to draw random points in polar coordinates in a circle (or a disk) of radius R
according to this function I(r)
, where r
would be the radius in the interval [0, R]
. The angle would be chosen uniformly in the interval [0, 2 * pi[
.
However, if I sample r
's straight away from I(r)
I get an accumulation of points close to the origin. How shall I properly sample my radii so that the radial distribution ends up being I(r)
?
At first, I didn't know how to draw samples from an arbitrary function such as I(r)
. Luckily, someone pointed me to this Python article, where Harry45 describes how draw samples from any distribution by subclassing the rv_continuous
class of scipy.stats
. Then, I generated a number of radii using the rvs(...)
method of rv_continuous
, and the same number of angles with numpy.random.rand
in the interval [0, 2 * pi[
. Finally, I transformed to Cartesian coordinates and plotted the results in a scatter plot. I had a feeling that too many points were clustered around the origin. I used a value b
of 1000, which is close to a uniform distribution (see my plot above). As it turns out, the circle wasn't uniformly filled, and more points were clustered around the origin. This is an issue that was addressed on Stackoverflow before at this link. Based on this link, I did the following to sample one point:
- Get a radius
r
fromI_norm(r) = A * exp ( -2 * r ^ b )
- Take the square root of
r
- Multiply
r
byw
- Get an angle from the uniform interval
[0, 2 * pi[
- Transform to Cartesian coordinates
At first, it seemed to work, but when I gave it more thought I realized. When I sample I_norm(r)
and b
is like 2, 4, or 8, there's a substantial chance to get a radius above 1. If the radius is above 1, taking the square root makes it even smaller and does the opposite of what we want (that is increase the number of points at larger radii to conserve the same denisty of points along the circle at that radius).
Moreover, with rv_continuous.rvs(...)
, I'm not sure how to limit myself to an interval [0, R]
. Naturally, the greater r
is the less likely it is to be sampled, but still, I think it remains an issue.
Long story short, I'm pretty sure what I'm doing currently is mathematically wrong, although I struggle to explain why. I'd appreciate if someone could shine some light on how to generate random points in a circle of radius R with a radial distribution that follows I(r)
.
Thank you in advance, and take care,
David
Aucun commentaire:
Enregistrer un commentaire