Taking points uniformly inside a general finite geometric domain

Question

It is well known that if we want to take $n$ uniformly and randomly points inside a circle of radius $r$ and centered at the origin the following apparently correct approach for generating $x$ and $y$ $$ x= U \cos(\theta), \;y= U \sin(\theta)$$ where $U$ is a uniform variate in $(0,1)$ and $\theta$ is uniform variate in $(0, 2 \pi)$, does not work. Rather the correct way to generate $x$ and $y$ coordinates is to use: $$ x= \sqrt{U}\cos(\theta),\; y= \sqrt{U} \sin(\theta) $$ with $U$ and $\theta $ as mentioned above. I wonder is there a some way of generating a of $n$ points inside the bounded region defined by $$\mathbf{D}=\{(x,y):f(x,y)=0\},$$ given that $ \mathbf{D}$ forms a closed region with smooth boundary. For instance let's say $$\mathbf{D}=\{ (x,y):\frac{x^2}{a^2}+\frac{y^2}{b^2}=1\},$$ how can we create $n$ points that are uniformly distributed inside the ellipse. Thanks for any hints/responses in advance!

Answer 1

For an ellipse, one can rescale the coordinates so that the region becomes a disk and then sample in the way you mentioned.

However, in general sampling efficiently from irregular regions (or distributions) is a really hard problem. If the dimension is low and the region is not too crazy, you can find some hypercube $C$ which contains $D$ and sample uniformly from $C$ while rejecting draws that don't belong to $D$. However, if the dimension is high or the shape is too irregular, you end up throwing away the vast majority of your draws, which really limits the effectiveness of this brute force approach.

Instead, the common thing to use in practice is some variant of Markov chain Monte-Carlo, which attempts to walk around the space in a random way. This is much more efficient to implement and by various ergodic theorems should converge to uniform sampling in the limit. Unfortunately, in practice it's very difficult to determine whether you've let it run long enough for the convergence to actually occur. Statisticians have developed a whole slew of heuristics for this problem, but it's very hard to say things rigorously.

Edit: in fact, it is possible to sample uniformly from an ellipsoid in any dimension by taking an affine change of coordinates so that the region is a ball. You then sample the radius $\sim r^{n-1}$. To determine the angle, you take $n$ independent draws from the standard univariate normal distribution. It turns out that the angle of the resulting vector will be equidistributed in $\mathbb{S}^n$, which you can use to choose an angle quickly. This example is a bit magical, and not at all what you should expect for more general regions.

Answer 2

The approach (already alluded to in one of the answers) of sampling uniformly from a larger set and then throwing out the samples that you don't want is known as rejection sampling. You'll find all sorts of useful ideas in the linked Wikipedia article and references thereof.

Answer 3

A quick way to generate a random point uniformly distributed in a bounded region $D$ is to generate a random point $P$ uniformly distributed in a rectangle $R$ containing $D$ and, if $P\notin R$, then discard $P$ and continue until you have as many random points as you want.

For instance, here is the generation (in Mathematica) of $3000$ random points uniformly distributed inside the ellipse centered at the origin with half-axes $a=2$ and $b=1$:

Here the "waste" fraction is $1-\frac{\pi ab}{4ab}\approx0.21$, about $21\%$, no problem at all.

---

Alternatively, one may generate a random point uniformly distributed in an arbitrary measurable plane region $D$ of positive area without any waste, as follows. For real $x$, let $$F(x):=F_D(x):=\frac{A(x)}{A(\infty)},$$ where $A(x)$ is the area of the region $\{(s,t)\in D\colon s\le x\}$. So, $F$ is a cumulative probability distribution function (cdf), which is actually the cdf of the abscissa of a random point uniformly distributed in $D$. For any $u\in(0,1)$, let $$F^{-1}(u):=\min\{x\in\mathbb R\colon F(x)\ge u\},$$ the quantile function corresponding to $F$. So, if $U$ is a random variable (r.v.) uniformly distributed on the interval $(0,1)$, then the distribution of the r.v. $$X:=F^{-1}(U)$$ will coincide with the distribution of the abscissa of a random point uniformly distributed in $D$.

If now the conditional distribution of a r.v. $Y$ given $X=x$ is the uniform distribution on the one-dimensional set $$D_x:=\{y\in\mathbb R\colon(x,y)\in D\},$$ then the random point $(X,Y)$ will be uniformly distributed in $D$.

So, the generation of a $(X,Y)$ uniformly distributed in $D$ is reduced to the generation of two random points on the real line.

For instance, here is the generation (in Mathematica) of $2000$ random points uniformly distributed in the ellipse $E:=\{(x,y)\in\mathbb R^2\colon100 (x - y)^2 + (x + y)^2 \le4\}$:

Here there is no "waste" at all, but the volume of calculations is much greater than in the previous example. In this particular case, it would be more economical to rotate the ellipse appropriately to make its axes horizontal and vertical and then use the approach of the previous example (without rotation, there could be too much waste).

---

In response to the comment by Timothy Budd, who wrote: "RandomPoint[Disk[{0, 0}, {a, b}], n] achieves the same result but is over a hundred times faster (for $n=3000$). Of course, this may be just due to low-level optimization."

I think the advantage of Mathematica's command RandomPoint[] over the function QQ[] defined above is mainly due to two things: (i) QQ[] produces (pseudo-)random points one-by-one, whereas RandomPoint[] apparently works with entire lists/arrays and (ii) RandomPoint[Disk[{0, 0}, {a, b}], n] takes into account the knowledge that the region is a (stretched) disk.

The image below of a Mathematica notebook shows that, when QQ[] is modified to a command QQQ[] operating on entire lists, RandomPoint[Disk[{0, 0}, {2, 1}], 3000] is only $0.0051967/0.0028878<2$ times faster than QQQ[2, 1, 3000]. Moreover, the command RandomPoint[ImplicitRegion[x^2/4 + y^2 <= 1, {x, y}], 3000], which does not let Mathematica know that the region is a (stretched) disk, is $0.157363/0.0051967>30$ times slower than QQQ[2, 1, 3000]: