metodo - método montecarlo python

Muestreo de puntos aleatorios distribuidos uniformemente dentro de un volumen esférico

Question

Estoy buscando ser capaz de generar una muestra aleatoria uniforme de ubicaciones de partículas que caen dentro de un volumen esférico.

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

Esto es lo que estoy recibiendo actualmente:

Puede ver que hay un grupo de puntos en el centro debido a la conversión entre coordenadas esféricas y cartesianas.

El código que estoy usando 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)

A continuación se muestra un código de MATLAB que supuestamente crea una muestra esférica uniforme, que es similar a la ecuación dada por http://nojhan.free.fr/metah . Simplemente no puedo descifrarlo ni 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);

Agradecería cualquier sugerencia sobre la generación de una muestra verdaderamente uniforme a partir de un volumen esférico en Python.

Parece que hay muchos ejemplos que muestran cómo extraer muestras de un caparazón esférico uniforme, pero parece ser más fácil y un problema más fácil. El problema tiene que ver con la escala: debe haber menos partículas en un radio de 0.1 que en un radio de 1.0 para generar una muestra uniforme del volumen de la esfera.

Editar: corrigió y eliminó el hecho que pedí normalmente y quise decir uniforme.

Answer 1

¿Sería esto lo suficientemente uniforme 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 porción de ella

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

parece:

Answer 2

El vector 3d normed gaussiano se distribuye uniformemente en la esfera, consulte 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))

Answer 3

Estoy de acuerdo con Alleo. Traduje tu código Matlab a Python y puede generar miles de puntos muy rápido (una fracción de segundo en mi computadora para 2D y 3D). Incluso lo he ejecutado hasta 5D hyperspheres. Encontré tu código tan útil que lo estoy aplicando en un estudio. Tim McJilton, ¿a quién debería 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'')

Answer 4

Genere un conjunto de puntos distribuidos uniformemente dentro de un cubo, luego descarte aquellos cuya distancia desde el centro excede el radio de la esfera deseada.

Answer 5

Hay una manera brillante de generar puntos uniformes en la esfera en el espacio n-dimensional, y usted ha señalado esto en su pregunta (me refiero al código MATLAB).

¿Por qué funciona? La respuesta es: veamos la densidad de probabilidad de la distribución normal n-dimensional. Es igual (hasta 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, ¡solo de la distancia! Esto significa que, después de normalizar el vector, la densidad de distribución resultante será constante en toda la esfera.

Este método debe ser definitivamente preferido debido a su simplicidad, generalidad y eficiencia (y 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 buen enlace para mirar: http://www-alg.ist.hokudai.ac.jp/~jan/randsphere.pdf

En cuanto a tener una distribución uniforme dentro de una esfera, en lugar de normalizar un vector, debes multiplicar vercor por alguna f (r): f (r) * r se distribuye con densidad proporcional a r ^ n en [0,1], que era hecho en el código que publicaste

Answer 6

Puede generar puntos aleatorios en coordenadas esféricas (suponiendo que está trabajando en 3D): S (r, θ, φ), donde r ∈ [0, R), θ ∈ [0, π], φ ∈ [0, 2π), donde R es el radio de tu esfera. Esto también le permitiría controlar directamente cuántos puntos se generan (es decir, no necesita descartar ningún punto).

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

Si su código necesita (x, y, z) (es decir, cartesianas) coordenadas, entonces simplemente convierta los puntos generados aleatoriamente en coordenadas esféricas a cartesianas como se explica here .

Answer 7

Si bien prefiero el método de descarte de esferas, para completarlo ofrezco la solución exacta .

En coordenadas esféricas, aprovechando la regla de muestreo :

phi = random(0,2pi)
costheta = random(-1,1)
u = random(0,1)

theta = arccos( costheta )
r = R * cuberoot( u )

ahora tienes un grupo (r, theta, phi) que se puede transformar a (x, y, z) de la manera habitual

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

Answer 8

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