¿Cómo generar puntos distribuidos uniformemente en la superficie de un elipsoide?

Question

Estoy tratando de encontrar una manera de generar puntos aleatorios distribuidos uniformemente en la superficie de un elipsoide.

Si se tratara de una esfera hay una forma muy clara de hacerlo: Generar tres $N(0,1)$ variables $\{x_1,x_2,x_3\}$ calcula la distancia al origen

$$d=\sqrt{x_1^2+x_2^2+x_3^2}$$

y calcular el punto

$$\mathbf{y}=(x_1,x_2,x_3)/d.$$

Puede demostrarse entonces que los puntos $\mathbf{y}$ se encuentran en la superficie de la esfera y se distribuyen uniformemente en la superficie de la esfera, y el argumento que lo demuestra es una sola palabra, "isotropía". No hay dirección preferida.

Supongamos ahora que tenemos un elipsoide

$$\frac{x_1^2}{a^2}+\frac{x_2^2}{b^2}+\frac{x_3^2}{c^2}=1$$

¿Qué le parece generar tres $N(0,1)$ como en el caso anterior, calcule

$$d=\sqrt{\frac{x_1^2}{a^2}+\frac{x_2^2}{b^2}+\frac{x_3^2}{c^2}}$$

y utilizando $\mathbf{y}=(x_1,x_2,x_3)/d$ como arriba. Así obtendremos puntos garantizados en la superficie del elipsoide, pero ¿estarán uniformemente distribuidos? ¿Cómo podemos comprobarlo?

Cualquier ayuda es muy apreciada, gracias.

PD Estoy buscando una solución sin aceptar/rechazar puntos, que es algo trivial.

EDITAR:

Pasando a coordenadas polares, el elemento de superficie es $dS=F(\theta,\phi)\ d\theta\ d\phi$ donde $F$ se expresa como $$\frac{1}{4} \sqrt{r^2 \left(16 \sin ^2(\theta ) \left(a^2 \sin ^2(\phi )+b^2 \cos ^2(\phi )+c^2\right)+16 \cos ^2(\theta ) \left(a^2 \cos ^2(\phi )+b^2 \sin ^2(\phi )\right)-r^2 \left(a^2-b^2\right)^2 \sin ^2(2 \theta ) \sin ^2(2 \phi )\right)}$$

Answer 1

Una forma de proceder es generar un punto uniformemente en la esfera, aplicar el mapeado $f : (x,y,z) \mapsto (x'=ax,y'=by,z'=cz)$ y luego corregir la distorsión creada por el mapa descartando el punto aleatoriamente con cierta probabilidad $p(x,y,z)$ (después de descartar se reinicia todo).

Cuando aplicamos $f$ un área pequeña $dS$ alrededor de algún punto $P(x,y,z)$ se convertirá en una pequeña zona $dS'$ en torno a $P'(x',y',z')$ y necesitamos calcular el factor multiplicativo $\mu_P = dS'/dS$ .

Necesito dos vectores tangentes alrededor de $P(x,y,z)$ así que elegiré $v_1 = (dx = y, dy = -x, dz = 0)$ y $v_2 = (dx = z,dy = 0, dz=-x)$

Tenemos $dx' = adx, dy'=bdy, dz'=cdz$ ; $Tf(v_1) = (dx' = adx = ay = ay'/b, dy' = bdy = -bx = -bx'/a,dz' = 0)$ y del mismo modo $Tf(v_2) = (dx' = az'/c,dy' = 0,dz' = -cx'/a)$

(podemos hacer una comprobación de cordura y calcular $x'dx'/a^2+ y'dy'/b^2+z'dz'/c^2 = 0$ en ambos casos)

Ahora, $dS = v_1 \wedge v_2 = (y e_x - xe_y) \wedge (ze_x-xe_z) = x(y e_z \wedge e_x + ze_x \wedge e_y + x e_y \wedge e_z)$ así que $|| dS || = |x|\sqrt{x^2+y^2+z^2} = |x|$

Y $dS' = (Tf \wedge Tf)(dS) = ((ay'/b) e_x - (bx'/a) e_y) \wedge ((az'/c) e_x-(cx'/a) e_z) = (x'/a)((acy'/b) e_z \wedge e_x + (abz'/c) e_x \wedge e_y + (bcx'/a) e_y \wedge e_z)$

Y finalmente $\mu_{(x,y,z)} = ||dS'||/||dS|| = \sqrt{(acy)^2 + (abz)^2 + (bcx)^2}$ .

Es rápido comprobar que cuando $(x,y,z)$ está en la esfera el extremo de esta expresión sólo puede darse en uno de los seis "polos" ( $(0,0,\pm 1), \ldots$ ). Si suponemos $0 < a < b < c$ su mínimo está en $(0,0,\pm 1)$ (donde la superficie se multiplica por $ab$ ) y el máximo está en $(\pm 1,0,0)$ (donde la superficie se multiplica por $\mu_{\max} = bc$ )

Cuanto menor sea el factor de multiplicación, más puntos tendremos que eliminar, por lo que después de elegir un punto $(x,y,z)$ uniformemente sobre la esfera y aplicando $f$ tenemos que mantener el punto $(x',y',z')$ con probabilidad $\mu_{(x,y,z)}/\mu_{\max}$ .

De este modo, obtendrá puntos distribuidos uniformemente en el elipsoide.

Answer 2

¿Y si utilizas el método de la distribución normal, como propones, pero utilizas los radios respectivos como varianzas? Entonces escala con d, como se propone. Empíricamente, esto no parece más agrupado (al menos para mí) en los "polos de cigarros" que una selección de fuerza bruta a partir de puntos distribuidos uniformemente en $\mathbb{R}^3$ .

Por lo tanto, el método es generar variables aleatorias:

$$x_1 \sim N(0,a^2)$$ $$x_2 \sim N(0,b^2)$$ $$x_3 \sim N(0,c^2)$$

El resto es como se sugiere:

$$d=\sqrt{\frac{x_1^2}{a^2}+\frac{x_2^2}{b^2}+\frac{x_3^2}{c^2}}$$

y puntos de uso $\mathbf{y}=(x_1,x_2,x_3)/d$

¿La prueba de la selección esférica original es simplemente que las distribuciones normales son radialmente simétricas? ¿No sigue siendo así con varianzas desiguales?

NOTA: Esto es diferente de la respuesta en el comentario de @jibounet que propone escalar los puntos que se generaron en la superficie de una esfera, porque los términos en $d$ son diferentes, es decir, para

$$x \sim N(0,n^2)$$ $$ \frac{x^2}{n^2} \neq N(0,1) $$

Answer 3

Esta semana he estado trabajando en este problema y me ha resultado eficaz el siguiente planteamiento.

El elipsoide está parametrizado por

$$\vec f(t, u) = \begin{bmatrix}a \cos u \cos t \\ b \cos u \sin t \\ c \sin u\end{bmatrix}, \quad \begin{align}t &\in [0, 2 \pi] \\ u &\in [0, \pi]\end{align}$$

Ingenuamente, podríamos elegir $T \sim U(0, 2 \pi)$ y $U \sim U(0, \pi)$ y elija $\vec f(T, U)$ como nuestro punto aleatorio, pero esto equivale a estirar la esfera, lo que crea distorsión como otros han señalado.

En su lugar, defina las regiones

$$\begin{align} D_{t'} :\quad & t \in [0, t'] , && u \in [0, \pi] \\ D_{u'} :\quad & t \in [0, 2 \pi] , && u \in [0, u'] \end{align}$$

Y las funciones

$$\begin{align} S_t(t') \equiv & \iint_{D_{t'}} dS \\ S_u(u') \equiv & \iint_{D_{u'}} dS \end{align}$$

Eso es, $S_t$ se obtiene la superficie acumulada en $t\in [0,t']$ y para todos $u$ mientras que $S_u$ hace lo contrario.

A continuación, elija dos números $X_1, X_2 \sim U(0,A)$ donde $A$ es la superficie total $$A = S_t(2\pi) = S_u(\pi)$$

Desde $S_t$ y $S_u$ son monotónicamente crecientes en el espacio de parámetros, pueden evaluarse numéricamente e invertirse mediante interpolación.

Nuestros parámetros aleatorios son entonces

$$\begin{align} T = S_t^{-1}(X_1) \\ U = S_u^{-1}(X_2) \\ \end{align}$$

Y un punto cualquiera de la superficie es $$\vec f(T, U)$$

La ventaja de este método es que, al no tener que descartar ningún punto, no necesitamos calcular un factor multiplicativo particular para esta elipse (como hace otro contestador). De hecho, este método elige puntos distribuidos uniformemente en cualquier superficie paramétrica cuya región de parámetros es cuadrada.

También es eficiente desde el punto de vista computacional, porque podemos utilizar los mismos datos (a saber, una cuadrícula de parámetros y los valores de la función en esos parámetros) para determinar la superficie global y, por interpolación, las funciones de superficie acumulativa y sus inversas.

He desarrollado un pequeño módulo de Python que aplica este enfoque como r_surface() . El ejemplo 3 del módulo utiliza un elipsoide:

Esta aplicación está en deuda con las respuestas a estas preguntas de la SE:

• Generar puntos aleatorios en el perímetro de la elipse (Matemáticas SE)
• Generar puntos aleatorios en el perímetro de la elipse (Code Review SE)

Answer 4

Aquí está el código, Esto funciona en CUALQUIER dimensión :

import numpy as np

dim = 2
r=1
a = 10
b = 4
A = np.array([[1/a**2, 0],
               [0, 1/b**2]])
L = np.linalg.cholesky(A).T 

x = np.random.normal(0,1,(200,dim)) 

z = np.linalg.norm(x, axis=1) 
z = z.reshape(-1,1).repeat(x.shape[1], axis=1) 
y = x/z * r
y_new = np.linalg.inv(L) @ y.T

plt.figure(figsize=(15,15))
plt.plot(y[:,0],y[:,1], linestyle="", marker='o', markersize=2)
plt.plot(y_new.T[:,0],y_new.T[:,1], linestyle="", marker='o', markersize=5)
plt.gca().set_aspect(1)
plt.grid()

Y aquí está la teoría : https://freakonometrics.hypotheses.org/files/2015/11/distribution_workshop.pdf

Answer 5

Idea para una solución aproximada: dividir el elipsoide en cuasi-rectángulos suficientemente pequeños y planos. $P_i$ elija una parametrización casi preservadora del área de cada pieza: $$f_i:R_i\longrightarrow P_i$$ con $R_i\subset\Bbb R^2$ un rectángulo. Elige al azar un índice $j$ con probabilidad $\text{area}(R_j)/\sum_i\text{area}(R_i)$ elija un punto en $R_j$ y finalmente aplicar $f_j$ .