Estoy buscando para ser capaces de generar un aleatoria uniforme de la muestra de partículas lugares que caen dentro de una esfera de volumen.

La imagen a continuación (cortesía de http://nojhan.free.fr/metah/) muestra de qué estoy buscando. Este es un corte a través de la esfera, mostrando una distribución uniforme de los puntos:

Muestreo distribuidos de manera uniforme puntos al azar dentro de una esfera de volumen

Esto es lo que actualmente estoy recibiendo:

Muestreo distribuidos de manera uniforme puntos al azar dentro de una esfera de volumen

Usted puede ver que hay un conjunto de puntos en el centro debido a la conversión entre coordenadas Cartesianas y esféricas.

El código que estoy utilizando es:

def new_positions_spherical_coordinates(self):
   radius = numpy.random.uniform(0.0,1.0, (self.number_of_particles,1)) 
   theta = numpy.random.uniform(0.,1.,(self.number_of_particles,1))*pi
   phi = numpy.arccos(1-2*numpy.random.uniform(0.0,1.,(self.number_of_particles,1)))
   x = radius * numpy.sin( theta ) * numpy.cos( phi )
   y = radius * numpy.sin( theta ) * numpy.sin( phi )
   z = radius * numpy.cos( theta )
   return (x,y,z)

Por debajo de un cierto código de MATLAB que supuestamente crea un uniforme esférica muestra, la cual es similar a la ecuación dada por http://nojhan.free.fr/metah. Me parece que no puede descifrar o entender lo que hicieron.

function X = randsphere(m,n,r)

% This function returns an m by n array, X, in which 
% each of the m rows has the n Cartesian coordinates 
% of a random point uniformly-distributed over the 
% interior of an n-dimensional hypersphere with 
% radius r and center at the origin.  The function 
% 'randn' is initially used to generate m sets of n 
% random variables with independent multivariate 
% normal distribution, with mean 0 and variance 1.
% Then the incomplete gamma function, 'gammainc', 
% is used to map these points radially to fit in the 
% hypersphere of finite radius r with a uniform % spatial distribution.
% Roger Stafford - 12/23/05

X = randn(m,n);
s2 = sum(X.^2,2);
X = X.*repmat(r*(gammainc(s2/2,n/2).^(1/n))./sqrt(s2),1,n);

Les agradecería mucho cualquier sugerencia en la generación de un verdadero uniforme de la muestra a partir de una esfera de volumen en Python.

Parece ser que hay un montón de ejemplos que muestran cómo muestra de un uniforme de cáscara esférica, pero que parece ser más fácil un problema más fácil. El problema tiene que ver con la ampliación de la escala no debe ser menor número de partículas en un radio de 0,1 que en un radio de 1.0 a generar una muestra uniforme del volumen de la esfera.

Edición: Fijo y quita el hecho de que me pidió normalmente y que me refería uniforme.

  • dependiendo de sus propósitos, también puede ser útil para mirar cuasi-aleatoria de números en lugar de (software) de números aleatorios
InformationsquelleAutor Tim McJilton | 2011-03-23

8 Comentarios

  1. 36

    Mientras que yo prefiero el método de descartes para las esferas, para la integridad Ofrezco la solución exacta.

    Esféricas en coordenadas, tomando ventaja de la regla de muestreo:

    phi = random(0,2pi)
    costheta = random(-1,1)
    u = random(0,1)
    
    theta = arccos( costheta )
    r = R * cuberoot( u )

    ahora usted tiene un (r, theta, phi) grupo que puede ser transformado a (x, y, z) en la forma habitual

    x = r * sin( theta) * cos( phi )
    y = r * sin( theta) * sin( phi )
    z = r * cos( theta )
    • Me siento como un idiota ahora. Todo lo que tenía que hacer era tomar la raíz cúbica de la radio? Perfecto gracias por su paciencia y continua de la parte en la discusión 🙂
    • Como dije en los comentarios de Jim respuesta, la mayoría de los libros prefieren el método de descarte para las esferas. Incluso con soporte de hardware raíces cúbicas tomar un par de ciclos, y la trigonometría necesario para volver a coordenadas Cartesianas costo algo de tiempo también. También tenga en cuenta que he escondido una segunda aplicación de este método de dibujo costheta de manera uniforme.
    • Ah bueno voy a conseguir. El costo de funcionamiento de los cubos de la raíz de todos ellos es mayor que el costo de todo lo que tiraba. Creo que voy a hacerlo de las dos maneras y me permiten decidir sobre la utilización de la función. De todos modos gracias de nuevo por tu ayuda!
    • ¿qué es R, u, y costheta. Estoy asumiendo R es la magnitud, pero ¿por qué es u generado aleatoriamente?
    • La variable u es simplemente un muñeco para hacer los pasos clara. La razón para tomar la raíz cúbica de un uniforme de la muestra está enterrado en el teorema fundamental de muestreo que explico en el vinculado a las preguntas. Estás en lo correcto acerca R (para la radio, porque soy de los físicos) y costheta es el coseno de theta (lanzado este modo, debido al teorema de muestreo de nuevo).
  2. 14

    Generar un conjunto de puntos distribuidos uniformemente dentro de un cubo, luego descartar aquellos cuya distancia desde el centro supera el radio de la esfera.

    • Buena llamada. El descarte es bastante eficiente para las esferas, y se recomienda en muchos textos de ser más rápido que la transformación necesaria para muestra exactamente.
    • Pensé en esa idea, sino que perdería ~50% del total de los puntos creados. Así que para crear 16.000 partículas se crearía ~32,000. Desde el Área de un cubo es de r^3 y la esfera es de 4/3*pi*(r/2)^3 = por lo que el ratio = ~4/8 = .5
    • Esto no debería resultar una distribución normal, que es lo que @Tim estaba pidiendo. La distribución Normal no es la misma que la distribución uniforme.
    • Oh @Juanchopanza me equivocaba. Me refería uniforme. Yo voy a ir a arreglar esto ahora.
    • McJilton: me di cuenta de lo que significaba. El cubo de la solución es buena! Si usted está preocupado acerca de la eficiencia, usted puede tratar de perfiles de cubo-descartar vs transformar.
    • que requieren un conocimiento de la verdadera transformación de uno 😉
    • Tenga en cuenta que esto no es eficiente para altas dimensiones, ya que el volumen de la unidad de la bola llega a cero (=probabilidad para el muestreo de un punto dentro de la bola).

  3. 13

    Hay una manera brillante para generar de manera uniforme los puntos sobre la esfera en el espacio n-dimensional, y se han señalado esta en tu pregunta (me refiero código de MATLAB).

    Por qué funciona? La respuesta es: veamos la densidad de probabilidad de n-dimensiones de la distribución normal. Es igual a (a constante)

    exp(-x_1*x_1/2) *exp(-x_2*x_2/2)… = exp(-r*r/2),
    por lo que no depende de la dirección, sólo en la distancia! Esto significa que, después de normalizar el vector, el resultado de la distribución de la densidad será constante en toda la esfera.

    Este método debe ser sin duda el preferido debido a su simplicidad, generalidad y la eficiencia (y la belleza).
    El código que genera 1000 eventos en la esfera en tres dimensiones:

    size = 1000
    n = 3 # or any positive integer
    x = numpy.random.normal(size=(size, n)) 
    x /= numpy.linalg.norm(x, axis=1)[:, numpy.newaxis]

    Por CIERTO, el enlace bueno para tener en cuenta: http://www-alg.ist.hokudai.ac.jp/~ene/randsphere.pdf

    Como para tener una distribución uniforme dentro de una esfera, en lugar de la normalización de un vector, se debe multiplicar vercor por algunos f(r): f(r)*r se distribuye con la densidad proporcional a r^n [0,1], que se hizo en el código que has publicado

  4. 2

    Sería uniforme suficiente para sus propósitos?

    In []: p= 2* rand(3, 1e4)- 1
    In []: p= p[:, sum(p* p, 0)** .5<= 1]
    In []: p.shape
    Out[]: (3, 5216)

    Una rebanada

    In []: plot(p[0], p[2], '.')

    parece:
    Muestreo distribuidos de manera uniforme puntos al azar dentro de una esfera de volumen

    • Esto es lo mismo que Jim Lewis sugirió a la derecha? La creación de un uniforme cubo y desecha cualquier cosa que no sea en la esfera?
    • McJilton: Sí, yo estaba escribiendo mientras Jim respuesta venga en y puesto que es el código que voy a decidido publicar de todos modos. De todos modos nada sugiere en su pregunta que la generación de más puntos realmente necesaria sería de algún modo problemático. Atención a los detalles? Gracias
    • Estoy corriendo una simulación de más de 16.000 estrellas con la velocidad y la posición de los lugares a los que quiero ser uniforme. Tenía la esperanza de tener una manera en la que puede establecer la cantidad de puntos a tener, ya que necesito exactamente 16,000 o lo que sea. Su forma de trabajo y supongo que podría seguir sumando puntos 1 por 1 hasta que de 16.000 cerrado.
    • McJilton: FWIW, en mi muy modesta máquina de generación de 1e5 3d uniforme de puntos al azar y descartando los de afuera (lo que da a unos 5.3e4 puntos), lleva a unos 35 ms. Si este tipo de interpretación no es aplicable a usted, por favor, danos más detalles. Gracias
    • Tal vez yo estaba molesto. No tuve la oportunidad de probar los dos.
    • McJilton: ahora viendo dmckee‘s respuesta creo que es la correct respuesta para esta situación (el número exacto de puntos en la esfera). Sin embargo, en general no se escala muy fácilmente a las dimensiones superiores (, si es que alguna vez la materia). Gracias

  5. 2

    Normativa de gauss vector 3d es distribuido uniformemente en la esfera, ver http://mathworld.wolfram.com/SpherePointPicking.html

    Por ejemplo:

    N = 1000
    v = numpy.random.uniform(size=(3,N)) 
    vn = v / numpy.sqrt(numpy.sum(v**2, 0))
    • El uniform debe ser normal. uniform significa una distribución uniforme mientras normal medio de la distribución gaussiana
  6. 2

    Estoy de acuerdo con Alleo. He traducido su código de Matlab para Python y puede generar miles de puntos muy rápido (una fracción de segundo en mi ordenador en 2D y 3D). Incluso he corrió hasta 5D hyperspheres.
    He encontrado en su código, por lo útil que me estoy aplicando en un estudio. Tim McJilton, que debo agregar como referencia?

    import numpy as np
    from scipy.special import gammainc
    from matplotlib import pyplot as plt
    def sample(center,radius,n_per_sphere):
        r = radius
        ndim = center.size
        x = np.random.normal(size=(n_per_sphere, ndim))
        ssq = np.sum(x**2,axis=1)
        fr = r*gammainc(ndim/2,ssq/2)**(1/ndim)/np.sqrt(ssq)
        frtiled = np.tile(fr.reshape(n_per_sphere,1),(1,ndim))
        p = center + np.multiply(x,frtiled)
        return p
    
    fig1 = plt.figure(1)
    ax1 = fig1.gca()
    center = np.array([0,0])
    radius = 1
    p = sample(center,radius,10000)
    ax1.scatter(p[:,0],p[:,1],s=0.5)
    ax1.add_artist(plt.Circle(center,radius,fill=False,color='0.5'))
    ax1.set_xlim(-1.5,1.5)
    ax1.set_ylim(-1.5,1.5)
    ax1.set_aspect('equal')

    Muestreo distribuidos de manera uniforme puntos al azar dentro de una esfera de volumen

  7. 0

    HTML:

    import numpy as np
    import matplotlib.pyplot as plt
    
    
    
    
    
    r= 30.*np.sqrt(np.random.rand(1000))
    #r= 30.*np.random.rand(1000)
    phi = 2. * np.pi * np.random.rand(1000)
    
    
    
    x = r * np.cos(phi)
    y = r * np.sin(phi)
    
    
    plt.figure()
    plt.plot(x,y,'.')
    plt.show()

    esto es lo que quieres

    • Usted está respondiendo este años más tarde. Hay dos respuestas correctas ya se ha mencionado aquí. Su ejemplo es para un uniforme círculo, no un uniforme de la esfera. Usted no tiene ninguna Coordenada Z… ?
  8. -1

    Sólo puede generar puntos al azar en coordenadas esféricas (suponiendo que usted está trabajando en 3D): S(r, θ, φ ), donde r ∈ [0, R), θ ∈ [0, π ], φ ∈ [0, 2π ), donde R es el radio de la esfera. Esto también permite controlar directamente la cantidad de puntos se generan (es decir, no es necesario descartar cualquier puntos).

    Para compensar la pérdida de la densidad con la radio, que generaría la coordenada radial siguiendo una ley de potencia de distribución (ver dmckee la respuesta de para una explicación sobre cómo hacer esto).

    Si su código de necesidades (x,y,z) (es decir, cartesiano) de las coordenadas, entonces usted acaba de convertir al azar los puntos de esférica a coordenadas cartesianas como se explica aquí.

    • Que es lo que mi código está haciendo. El problema ocurre cuando se convierte en polar de coordenadas tiene una cantidad de partículas en radio R, como en la radio r < R que hace que el centro sea más denso que el de los bordes
    • Para realizar este trabajo, que debe dibujar la posición radial no uniformaly. Debido a que el elemento de volumen es r^2 dr d\phi d(cos\theta). Y esta consiste en extraer una raíz cúbica y un coseno inverso, por lo que tienden a ser más lento que el descarte procedimiento.
    • Interesante. Bien supongo que me va a chupar y hacer el montecarlo estilo de uno. Yo todavía soy curioso, aunque de un ejemplo de hacerlo de la otra manera. Veo el código de otras personas yo no lo entiendo -_-.
    • Lo consiguió. Luego de descartar los puntos es probablemente más fácil. Si usted todavía desea hacerlo sin descartar puntos, deberá generar el radial coordinar siguiente el poder de la ley de distribución de. Por desgracia, el muestreo de los no-trivial de las distribuciones no es fácil, pero si usted todavía está interesado, consulte el Metropolis-Hastings algoritmo para un método general (cualquier otro tipo de MCMC método también funcionaría).
    • No, repito no uso de la Metrópoli para esto! El muestreo de las distribuciones de ley de potencia es fácil, simplemente no es trivial.
    • Gracias @dmckee. Estoy interesado en esto mismo. ¿Sabes de algún fuentes de donde puedo aprender más acerca de muestreo, a partir de estas distribuciones?
    • A veces es llamada la Ley Fundamental de Muestreo. Me discutir en otros contextos en los dos enlaces en mi respuesta aquí. Normalizar el PDF, integrar, a continuación, invertir, y utilizar el resultado de la función de transformar los valores de dibujado de manera uniforme sobre [0,1). En este caso se trata de tomar una raíz cúbica.
    • Excelente. Gracias!

Dejar respuesta

Please enter your comment!
Please enter your name here