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
% generate radii of hyperspheres
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
% with radii rs
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

*Related*

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

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

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

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

It seems similar to what is done here, even though I was hoping for something even simpler:

https://ieeexplore.ieee.org/document/758215