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 

% 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

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

  1. 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)
        
        # generate radii of hyperspheres
        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 
        # with radii rs
        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
    
  2. 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

Leave a Reply