Genera una muestra aleatoria de puntos distribuidos en la superficie de una esfera unitaria

Estoy tratando de generar puntos aleatorios en la superficie de la esfera usando numpy. He revisado el post que explica la distribución uniforme aquí . Sin embargo, necesita ideas sobre cómo generar los puntos solo en la superficie de la esfera. Tengo coordenadas (x, y, z) y el radio de cada una de estas esferas.

No estoy muy versado en matemáticas en este nivel y estoy tratando de darle sentido a la simulación de Monte Carlo.

Cualquier ayuda será muy apreciada.

Gracias parin

Basado en el último enfoque en esta página , simplemente puede generar un vector que consiste en muestras independientes de tres distribuciones normales estándar, luego normalizar el vector de tal manera que su magnitud sea 1:

import numpy as np def sample_spherical(npoints, ndim=3): vec = np.random.randn(ndim, npoints) vec /= np.linalg.norm(vec, axis=0) return vec 

Por ejemplo:

 from matplotlib import pyplot as plt from mpl_toolkits.mplot3d import axes3d phi = np.linspace(0, np.pi, 20) theta = np.linspace(0, 2 * np.pi, 40) x = np.outer(np.sin(theta), np.cos(phi)) y = np.outer(np.sin(theta), np.sin(phi)) z = np.outer(np.cos(theta), np.ones_like(phi)) xi, yi, zi = sample_spherical(100) fig, ax = plt.subplots(1, 1, subplot_kw={'projection':'3d', 'aspect':'equal'}) ax.plot_wireframe(x, y, z, color='k', rstride=1, cstride=1) ax.scatter(xi, yi, zi, s=100, c='r', zorder=10) 

introduzca la descripción de la imagen aquí

El mismo método también se generaliza para seleccionar puntos distribuidos uniformemente en el círculo unitario ( ndim=2 ) o en las superficies de las hiperesferas unitarias de mayor dimensión.

Los puntos en la superficie de una esfera se pueden express utilizando dos coordenadas esféricas, theta y phi , con 0 < theta < 2pi y 0 < phi < pi .

Fórmula de conversión en coordenadas cartesianas x, y, z :

 x = r * cos(theta) * sin(phi) y = r * sin(theta) * sin(phi) z = r * cos(phi) 

donde r es el radio de la esfera.

Así que el progtwig podría muestrear aleatoriamente theta y phi en sus rangos, en una distribución uniforme, y generar las coordenadas cartesianas a partir de él.

Pero luego los puntos se distribuyen más densamente en los polos de la esfera. Para que los puntos se distribuyan uniformemente en la superficie de la esfera, phi debe elegirse como phi = acos(a) donde -1 < a < 1 se elige en una distribución uniforme.

Para el código Numpy sería igual que en Muestreo de puntos aleatorios distribuidos uniformemente dentro de un volumen esférico , excepto que el radius variable tiene un valor fijo.

Después de una discusión con @Soonts, sentí curiosidad por el rendimiento de los tres enfoques utilizados en las respuestas: uno con la generación de angularjs aleatorios, uno que usa coordenadas distribuidas normalmente y uno que rechaza puntos distribuidos uniformemente.

Aquí está mi bash de comparación:

 import numpy as np def sample_trig(npoints): theta = 2*np.pi*np.random.rand(npoints) phi = np.arccos(2*np.random.rand(npoints)-1) x = np.cos(theta) * np.sin(phi) y = np.sin(theta) * np.sin(phi) z = np.cos(phi) return np.array([x,y,z]) def sample_normals(npoints): vec = np.random.randn(3, npoints) vec /= np.linalg.norm(vec, axis=0) return vec def sample_reject(npoints): vec = np.zeros((3,npoints)) abc = 2*np.random.rand(3,npoints)-1 norms = np.linalg.norm(abc,axis=0) mymask = norms<=1 abc = abc[:,mymask]/norms[mymask] k = abc.shape[1] vec[:,0:k] = abc while k 

Luego por 1000 puntos

 In [449]: timeit sample_trig(1000) 1000 loops, best of 3: 236 µs per loop In [450]: timeit sample_normals(1000) 10000 loops, best of 3: 172 µs per loop In [451]: timeit sample_reject(1000) 100 loops, best of 3: 13.7 ms per loop 

Tenga en cuenta que en la implementación basada en el rechazo primero npoints muestras y npoints los malos, y solo usé un bucle para generar el rest de los puntos. Parecía ser el caso que el rechazo directo paso a paso lleva más tiempo. También eliminé la verificación de división por cero para tener una comparación más sample_normals con el caso sample_normals .


Eliminar la vectorización de los dos métodos directos los coloca en el mismo campo de juego:

 def sample_trig_loop(npoints): x = np.zeros(npoints) y = np.zeros(npoints) z = np.zeros(npoints) for k in xrange(npoints): theta = 2*np.pi*np.random.rand() phi = np.arccos(2*np.random.rand()-1) x[k] = np.cos(theta) * np.sin(phi) y[k] = np.sin(theta) * np.sin(phi) z[k] = np.cos(phi) return np.array([x,y,z]) def sample_normals_loop(npoints): vec = np.zeros((3,npoints)) for k in xrange(npoints): tvec = np.random.randn(3) vec[:,k] = tvec/np.linalg.norm(tvec) return vec 
 In [464]: timeit sample_trig(1000) 1000 loops, best of 3: 236 µs per loop In [465]: timeit sample_normals(1000) 10000 loops, best of 3: 173 µs per loop In [466]: timeit sample_reject(1000) 100 loops, best of 3: 14 ms per loop In [467]: timeit sample_trig_loop(1000) 100 loops, best of 3: 7.92 ms per loop In [468]: timeit sample_normals_loop(1000) 100 loops, best of 3: 10.9 ms per loop 

Otra forma en que dependiendo del hardware podría ser mucho más rápido.

Elija a, b, c para que sean tres números aleatorios, cada uno entre -1 y 1

Calcula r2 = a^2 + b^2 + c^2

Si r2> 1.0 (= el punto no está en la esfera) o r2 <0.00001 (= el punto está demasiado cerca del centro, tendremos división por cero mientras proyectamos a la superficie de la esfera) descartamos los valores , y elige otro conjunto de al azar a, b, c

De lo contrario, tienes tu punto aleatorio (relativo al centro de la esfera):

 ir = R / sqrt(r2) x = a * ir y = b * ir z = c * ir 

(editado para reflejar las correcciones de los comentarios)

Investigué algunos enfoques de tiempo constante para este problema en 2004.

asumiendo que está trabajando en coordenadas esféricas donde theta es el ángulo alrededor del eje vertical (por ejemplo, longitud) y phi es el ángulo elevado desde el ecuador (por ejemplo, latitud), para obtener una distribución uniforme de puntos aleatorios en el hemisferio norte. el ecuador haces esto

  1. elige theta = rand (0, 360).
  2. elija phi = 90 * (1 – sqrt (rand (0, 1))).

para obtener puntos en una esfera en lugar de un hemisferio, simplemente niega phi 50% del tiempo.

para los curiosos, se aplica un enfoque similar para generar puntos distribuidos uniformemente en una unidad de disco:

  1. elige theta = rand (0, 360).
  2. elija el radius = sqrt (rand (0, 1)).

No tengo pruebas de la corrección de estos enfoques, pero los he utilizado con mucho éxito en la última década y estoy convencido de que son correctos.

Aquí se incluye una ilustración (a partir de 2004) de los diversos enfoques, que incluye una visualización del enfoque de elección de puntos en la superficie de un cubo y su normalización en la esfera.