Generar puntos aleatorios dentro de un elipsoide de N dimensiones

Question

Me gustaría generar puntos distribuidos uniformemente dentro de un elipsoide de N dimensiones, donde los ejes del elipsoide están orientados aleatoriamente con respecto a los ejes cartesianos, y los medios a lo largo de los ejes del elipsoide no se encuentran necesariamente en los ejes cartesianos. Por ejemplo, en 2D, bastaría con algo como lo siguiente:

Elipsoide 2D

Me gustaría ampliarlo a cualquier número de dimensiones, sin embargo, y me gustaría poder establecer los medios y las direcciones de los ejes del elipsoide. Hay muchos tutoriales en línea para generar puntos aleatorios en un disco o dentro de una esfera, pero no es obvio cómo extenderlos a un elipsoide.

¿Alguna idea?

Answer 1

Tal vez esto papel ayudará.

No he entendido del todo el "cómo", pero el enfoque del artículo fue útil para generar muestras distribuidas uniformemente en elipsoides alargados y de alta dimensión. El código matlab para el proceso de generación de muestras también se incluye en el artículo. El código requiere el vector central y la forma matricial (inversa) del elipsoide como entradas.

La generación uniforme de muestras da lugar a un interior de 3 elipsoides

También me di cuenta de que establecer el término "Gamma_Threshold" como 1 aseguraba que todas las muestras se generaran dentro del interior del elipsoide.

Answer 2

Basándome en @optimum_jay, escribí una implementación en python del código de ese papel . El algoritmo funciona bastante bien para mis propósitos, que tienen alrededor de 30 dimensiones, pero los puntos parecen empezar a agruparse alrededor del punto medio a medida que el número de dimensiones aumenta. Esto probablemente se puede "arreglar rápidamente" aumentando el parámetro Gamma_Threshold, pero es algo a tener en cuenta. Las fotos se incluyen a continuación.

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)