Generate random points within N-dimensional ellipsoid I'd like to generate uniformly distributed points within an N-dimensional ellipsoid, where the ellipsoid axes are randomly oriented with respect to the Cartesian axes, and the means along the ellipsoid axes do not necessarily lie on the Cartesian axes. For instance, in 2D, something like the following would suffice:
2D ellipsoid
(I don't have enough reputation to post images yet.) I would like to extend this to any number of dimensions, though, and I would like to be able to set the means and directions of the ellipsoid axes. There are many tutorials online for generating random points on a disk or within a sphere, but it's not obvious how to extend these to an ellipsoid.
Any ideas?
 A: Maybe this paper will help. 
I've not entirely understood the "how" but the paper's approach was useful in generating uniformly distributed samples in both elongated & high-dimensional ellipsoids. The matlab code for the sample generation process is also included in the paper. The code requires the centre vector & (inverse of) ellipsoid matrix form as inputs.
Uniform sample generation results in 3-ellipsoid interior
I also noticed that setting the term "Gamma_Threshold" as 1 ensured that all samples are generated within the ellipsoid interior.
A: Building on @optimum_jay, I wrote a python implementation of the code from that paper. The algorithm works fairly well for my purposes which has around 30 dimensions, but the points seem to start clustering around the midpoint as the number of dimensions increase. This can probably be 'quick-fixed' by increasing the Gamma_Threshold parameter, but it's something to keep in mind. Pics included below.
import numpy as np
from scipy.linalg import cholesky # computes upper triangle by default, matches paper


def sample(S, z_hat, m_FA, Gamma_Threshold=1.0):

    nz = S.shape[0]
    z_hat = z_hat.reshape(nz,1)

    X_Cnz = np.random.normal(size=(nz, m_FA))

    rss_array = np.sqrt(np.sum(np.square(X_Cnz),axis=0))
    kron_prod = np.kron( np.ones((nz,1)), rss_array)

    X_Cnz = X_Cnz / kron_prod       # Points uniformly distributed on hypersphere surface

    R = np.ones((nz,1))*( np.power( np.random.rand(1,m_FA), (1./nz)))

    unif_sph=R*X_Cnz;               # m_FA points within the hypersphere
    T = np.asmatrix(cholesky(S))    # Cholesky factorization of S => S=T’T


    unif_ell = T.H*unif_sph ; # Hypersphere to hyperellipsoid mapping

    # Translation and scaling about the center
    z_fa=(unif_ell * np.sqrt(Gamma_Threshold)+(z_hat * np.ones((1,m_FA))))

    return np.array(z_fa)




