# Drawing uniform random numbers from within an ellipsoid

Here’s a Matlab function I wrote to draw a set of random numbers from a uniform distribution within an n-dimensional ellipsoid, where the size and orientation of the ellipsoid is defined by a covariance matrix.

function pnts = draw_from_ellipsoid( covmat, cent, npts )

% function pnts = draw_from_ellipsoid( covmat, cent, npts )
%
% This function draws points uniformly from an n-dimensional ellipsoid
% with edges and orientation defined by the the covariance matrix covmat.
% The number of points produced in the n-dimensional space is given by
% npts, with an output array of npts by ndims. The centre of each dimension
% is given in the vector cent.

% get number of dimensions of the covariance matrix
ndims = length(covmat);

% calculate eigenvalues and vectors of the covariance matrix
[v, e] = eig(covmat);

% check size of cent and transpose if necessary
sc = size(cent);
if sc(1) > 1
cent = cent';
end

rs = rand(npts,1);

% generate points
pt = randn(npts,ndims);

% get scalings for each point onto the surface of a unit hypersphere
fac = sum(pt(:,:)'.^2);

% calculate scaling for each point to be within the unit hypersphere
fac = (rs.^(1/ndims)) ./ sqrt(fac');

pnts = zeros(npts,ndims);

% scale points to the ellipsoid using the eigenvalues and rotate with
% the eigenvectors and add centroid
d = sqrt(diag(e));
for i=1:npts
% scale points to a uniform distribution within unit hypersphere
pnts(i,:) = fac(i)*pt(i,:);

% scale and rotate to ellipsoid
pnts(i,:) = (pnts(i,:) .* d' * v') + cent;
end

end


## 5 Replies to “Drawing uniform random numbers from within an ellipsoid”

1. Krishna says:

Thank you so much for the article Matthew. Really helped me a lot 🙂

2. John A. Major says:

Thanks – I needed that and translated it into R. Will send you code if you want.

3. Thanks very much for this. In case anyone wants the code for Python, here you go:

def draw_from_ellipsoid( covmat, cent, npts):
# random uniform points within ellipsoid as per: http://www.astro.gla.ac.uk/~matthew/blog/?p=368
ndims = covmat.shape[0]

# calculate eigenvalues (e) and eigenvectors (v)
eigenValues, eigenVectors = np.linalg.eig(covmat)
idx = (-eigenValues).argsort()[::-1][:ndims]
e = eigenValues[idx]
v = eigenVectors[:,idx]
e = np.diag(e)

rs = np.random.uniform(0,1,npts)

# generate points
pt = np.random.normal(0,1,[npts,ndims]);

# get scalings for each point onto the surface of a unit hypersphere
fac = np.sum(pt**2,axis=1)

# calculate scaling for each point to be within the unit hypersphere
fac = (rs**(1.0/ndims)) / np.sqrt(fac)

pnts = np.zeros((npts,ndims));

# scale points to the ellipsoid using the eigenvalues and rotate with
# the eigenvectors and add centroid
d = np.sqrt(np.diag(e))
d.shape = (ndims,1)

for i in range(0,npts):
# scale points to a uniform distribution within unit hypersphere
pnts[i,:] = fac[i]*pt[i,:]
pnts[i,:] = np.dot(np.multiply(pnts[i,:],np.transpose(d)),np.transpose(v)) + cent

return pnts

4. Felix says:

Hi, thanks a lot! Can you point me to the underlying theory?

I get the intuition:

You sample points from a normal distribution and normalize them to be on the unit sphere. That way, you get points with a uniform distribution in their orientation. Then you combine the uniform orientations with uniform radii (is that results still uniform?). Where is the exponent in r^*(1/ndim) coming from?
Finally, you rotate and rescale.

Looking forward to some hint,

Felix