python - two - Genera una muestra aleatoria de puntos distribuidos en la superficie de una unidad de esfera

random python

Question

Intento generar puntos aleatorios en la superficie de la esfera usando numpy. Revisé la publicación 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 bien versado en Matemáticas en este nivel y trato de darle sentido a la simulación de Monte Carlo.

Cualquier ayuda será muy apreciada.

Gracias, Parin

Answer 1

Basado en el último enfoque de esta página , puede simplemente generar un vector que consta de muestras independientes de tres distribuciones normales estándar, luego normalizar el vector de 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)

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

Answer 2

Después de un debate con @Soonts, me llamó la atención el rendimiento de los tres enfoques utilizados en las respuestas: uno con la generación de ángulos aleatorios, uno con coordenadas normalmente distribuidas y otro que rechaza los puntos distribuidos uniformemente.

Aquí está mi intento 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<npoints:
       abc = 2*np.random.rand(3)-1
       norm = np.linalg.norm(abc)
       if 1e-5 <= norm <= 1:  
           vec[:,k] = abc/norm
           k = k+1
    return vec

Entonces 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 rechazo primero npoints muestras de npoints y npoints las malas, y solo usé un bucle para generar el resto de los puntos. Parecía que el rechazo directo paso a paso requería una mayor cantidad de tiempo. También eliminé el cheque para la división por cero para tener una comparación más sample_normals con el caso de sample_normals .

Eliminar la vectorización de los dos métodos directos los coloca en el mismo estadio de béisbol:

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

Answer 3

Los puntos en la superficie de una esfera se pueden expresar usando 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.

Entonces, el programa 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 de Numpy sería el mismo que en Sampling puntos aleatorios distribuidos uniformemente dentro de un volumen esférico , excepto que el radius variable tiene un valor fijo.

Answer 4

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

Calcular 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 se proyecta a la superficie de la esfera) usted descarta los valores , y elige otro conjunto aleatorio 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

Answer 5

(editado para reflejar las correcciones de los comentarios)

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

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

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

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

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

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

No tengo pruebas de la exactitud de estos enfoques, pero los he utilizado con mucho éxito durante la última década y estoy convencido de su corrección.

alguna ilustración (desde 2004) de los diversos enfoques está aquí , incluida una visualización del enfoque de elegir puntos en la superficie de un cubo y normalizarlos en la esfera.