Generate random uniform points in an ellipse

I have previously written about how to efficiently generate points uniformly at random inside a sphere (often called a ball by mathematicians). The method uses a mathematical fact from multivariate statistics: If X is drawn from the uncorrelated multivariate normal distribution in dimensiond, then S = r*X / ||X|| has the uniform distribution in the unit ball, where r ~ (1/U)(1/d) and U ~ U(0,1). This is true in any dimension; see the previous article for pictures of a uniform sample inside a ball in 3-D.

It is worth mentioning that this procedure can also be used to generate random points uniformly in any ellipse. This article shows the mathematical algorithm and implements it in SAS. Examples are shown in 2-D (an ellipse in the plane) and in 3-D (an ellipsoid). The same process works for hyperellipsoids in high-dimensional Euclidean space.

From spheres to ellipsoids

Intuitively, you can move from spheres to ellipsoids by changing the covariance matrix that you use for the multivariate normal (MVN) distribution. When the covariance matrix is the identity matrix, the equi-probability contours of the density are spheres. However, for a general positive definite covariance matrix, the equi-probability contours of the density are (hyper)ellipsoids. A particularly simple choice for a covariance matrix is a diagonal matrix. For a diagonal matrix, the axes of the density ellipsoids are orthogonal and are aligned with the standard coordinate axes.

There is a simple geometric relationship between the equi-probability contours for an uncorrelated MVN distribution and for a correlated MVN distribution with covariance matrix Σ. Briefly, if you use the Cholesky factorization to decompose Σ = LLT, then the matrix L is a linear transformation that maps the spherical contours of the uncorrelated distribution to corresponding contours of the correlated distribution, as shown by the image to the right.

The geometry (and the computation) is simplest when you use a diagonal covariance matrix. If \(\Sigma = \mbox{diag}(\sigma_1^2, \sigma_2^2, \ldots, \sigma_d^2)\) is a diagonal covariance matrix, then L = diag(σ1, σ2, ..., σd) for σi > 0. The linear transformation, L, maps the unit sphere at the origin to an ellipsoid with axes that are aligned with the standard coordinate axes and whose radii are σ1, σ2, ..., σd.

The uniform distribution is preserved under a linear transformation

Of course, this problem requires more than mapping a ball to an ellipsoid. We want to generate uniform random variates in an ellipsoid. Fortunately, you can prove that the uniform distribution in a region is preserved under a linear transformation. More precisely, suppose that X is a continuous random variable with a uniform distribution on some subset S in d-dimensional Euclidean space, where S has positive measure. Let Y = μ + L*X be any affine transformation of X, where μ is a d-dimensional vector and L is an invertible d x d matrix. Then Y is uniformly distributed on the image of S, which is the set {y = μ + L*x | x ∈ S}.

Notice that the set S must be a set of positive measure. For example, it must be an area in 2-D, a volume in 3-D, and so forth.

This theorem tells you that you can generate random uniform variates in the unit ball (as shown in the previous article) and then linearly transform those variates to obtain uniform variates in an ellipsoid. The transformed variates will be uniformly distributed in the image of the ball, which is an ellipsoid.

A SAS program to generate random uniform variates in an ellipse

You can carry out this process in any matrix-vector language, such as R, MATLAB, or SAS IML. The following SAS IML program defines a function that takes two arguments: the number of points (n) to generate, and a vector of d positive values, which represent the elements of a diagonal covariance matrix. It returns n points in d-dimensional space that are in the ellipse {y = diag(σ1, σ2, ..., σd)*x | x is in the unit ball and σi > 0}.

/* Uniform random points in an ellipsoid in d-dimensional space.
The ellipsoid is the image of the unit ball under the transformation x -> L*x
where x is in the ball and L = diag(sigma_1, sigma_2, ..., sigma_d), sigma_i > 0.
*/
proc iml;
/* generate n random uniform points in the ellipse L*x, where x is in unit ball and L = diag(s) */
start RandEllipse(n, s);
d = ncol(s);
mu = j(1, d, 0);
X = randnormal(n, mu, I(d)); /* X ~ N(0, I) */
u = randfun(n, "Uniform"); /* U ~ U(0,1) */
r = u##(1/d); /* r proportional to d_th root of U */
Z = r # X / sqrt(X[,##]); /* uniform in disk or ball */
L = diag(s); /* Cholesky root of covariance matrix Sigma=diag(s##2) */
return Z * L`; /* uniform in ellipse/ellipsoid */
finish;

call randseed(12345);
/* 2-D example with n=200 points in an ellipse in the plane */
n = 1000;
sigmas = {4 3};
Y = RandEllipse(n, sigmas);
title "1,000 Random Points in Ellipse";
call scatter(Y[,1], Y[,2]) grid={x y}
procopt="aspect=0.75" option="transparency=0.5";

A scatter plot of the output shows that the points are distributed uniformly at random in the ellipse with one axis aligned with the vector (4,0) and the other axis aligned with the vector (0,3).

The same function supports generating points in higher-dimensional ellipsoids. For example, to get random points in an ellipsoid in 3-D, you can use the following statements.

/* 3-D example with points in an ellipsoid */
Y = RandEllipse(5000, {4 3 1}); /* axes aligned with (4,0,0), (0,3,0), and (0,0,1) */
/* write points to a data set and visualize */
create RandEllipsoid from Y[c={'x1' 'x2' 'x3'}];
append from Y;
close;
QUIT;

title "5,000 Random Points in an Ellipsoid";
proc sgplot data=RandEllipsoid;
scatter x=x1 y=x2 / colorresponse=x3 transparency=0.5 markerattrs=(symbol=CircleFilled);
xaxis grid;
yaxis grid;
run;

The graph projects the points onto the (x,y) plane. The markers are color-coded according to the z value, which ranges in [-1,1]. Because the points are in an ellipsoid, points that are close to z=0 are located near the planar ellipse that passes through (±4, 0) and (0, ±3). These points are gray in color. Points near the middle of this graph are colored blue, gray, or red, because the z value varies.

Ellipsoids not aligned with the coordinate axis

The same technique works even if you do not use a diagonal covariance matrix. Given any positive definite matrix, you can find the Cholesky decomposition Σ = LLT and use L to map random variate from the unit disk into an ellipse. The following SAS IML function implements this general technique:

proc iml;
/* generate n random uniform points in the ellipse L*x, where x is in the
unit ball and L satisfies Sigma = L*L` */
start RandEllipseCov(n, Sigma);
d = ncol(Sigma);
mu = j(1, d, 0);
X = randnormal(n, mu, I(d)); /* X ~ N(0, I) */
u = randfun(n, "Uniform"); /* U ~ U(0,1) */
r = u##(1/d); /* r proportional to d_th root of U */
Z = r # X / sqrt( X[,##] ); /* Z is uniform in disk or ball */
U = root(Sigma); /* Cholesky decomp Sigma=U`*U, U upper triangular */
L = U`; /* Sigma=L*L', L lower triangular */
return Z * L`; /* points are uniform in ellipse/ellipsoid */
finish;

call randseed(12345);
/* 2-D example with n=1,000 points in ellipse in plane */
n = 1000;
Sigma = {16 4,
4 9};
Y = RandEllipseCov(n, Sigma);
title "1,000 Random Points in a Tilted Ellipse";
/* slopes of reference lines are the directions of the eigenvectors of Sigma */
call scatter(Y[,1], Y[,2]) grid={x y}
procopt="aspect=0.75 noautolegend" option="transparency=0.5"
other="lineparm x=0 y=0 slope=0.4537682/ lineattrs=(color=gray); lineparm x=0 y=0 slope=-2.203768/lineattrs=(color=gray);";

The graph overlays two reference lines, which indicate the axes for the ellipse. They are computed by using the eigenvectors of the Σ matrix. The computation is not shown here, but you can read the details in a previous article about ellipses and covariance matrices.

Summary

This article shows how to generate points uniformly at random in an ellipse, in an ellipsoid, or in any hyperellipsoid in higher dimensions. It first generates points uniformly at random in the unit ball, then uses a linear transformation to map the points into an ellipsoid. If you want an ellipse whose axes are aligned with the coordinate axes, the linear transformation is given by a simple diagonal matrix with positive values on the diagonal. You can also generate points in ellipsoids whose axes are not aligned with the coordinate axes. If Σ is any positive definite matrix, you can decompose Σ = LLT. If you generate points in the unit ball, you can use L to map them into the ellipse whose axes are aligned with the eigendirections of Σ.