Выборка равномерно распределенных случайных точек внутри сферического объема

Question

Я надеюсь, что смогу сгенерировать случайную однородную выборку местоположений частиц, которые попадают в сферический объем.

Изображение ниже (любезно предоставлено http://nojhan.free.fr/metah/) показывает то, что я ищу.Это кусочек сферы, показывающий равномерное распределение точек:

Вот что я сейчас получаю:

Вы можете видеть, что в центре есть группа точек из-за преобразования между сферическими и декартовыми координатами.

Код, который я использую:

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)

Ниженекоторый код MATLAB, который предположительно создает однородную сферическую выборку, которая похожа на уравнение, заданное http://nojhan.free.fr/metah. Я просто не могу расшифровать его или понять, что они сделали.

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);

IБуду очень признателен за любые предложения по созданию действительно однородного образца из сферического тома в Python.

Кажется, есть множество примеров, показывающих, как производить выборку из однородной сферической оболочки,но это, кажется, легче, легче проблема.Проблема связана с масштабированием - в радиусе 0,1 должно быть меньше частиц, чем в радиусе 1,0, чтобы получить однородный образец из объема сферы.

Править: Исправлено и убрано то, что я обычно просил и имел в виду униформу.

Answer 1

Хотя я предпочитаю метод отбрасывания для сфер, для полноты я предлагаю точное решение .

В сферических координатах, используя преимущество правила выборки :

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

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

теперь у вас есть группа (r, theta, phi), которую можно преобразовать в (x, y, z) обычным способом

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

Answer 2

Создайте набор точек, равномерно распределенных внутри куба, затем отбросьте те, чье расстояние от центра превышает радиус желаемой сферы.

Answer 3

Существует блестящий способ создания равномерных точек на сфере в n-мерном пространстве, и вы указали это в своем вопросе (я имею в виду код MATLAB).

Почему это работает? Ответ таков: посмотрим на плотность вероятности n-мерного нормального распределения. Равно (до постоянного)

exp (-x_1 * x_1 / 2) * exp (-x_2 * x_2 / 2) ... = exp (-r * r / 2), так что это не зависит от направления, только от расстояния! Это означает, что после нормализации вектора полученная плотность распределения будет постоянной по всей сфере.

Этот метод определенно предпочтителен из-за его простоты, универсальности и эффективности (и красоты). Код, который генерирует 1000 событий на сфере в трех измерениях:

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

Кстати, хорошая ссылка для просмотра: http://www -alg.ist.hokudai.ac.jp / ~ jan / randsphere.pdf

Что касается равномерного распределения в пределах сферы, вместо нормализации вектора вы должны умножить vercor на некоторое f (r): f (r) * r распределено с плотностью, пропорциональной r ^ n на [0,1], что было сделано в коде, который вы разместили

Answer 4

Я согласен с Alleo. Я перевел ваш код Matlab на Python, и он может очень быстро генерировать тысячи точек (доли секунды в моем компьютере для 2D и 3D). Я даже запускал его для 5D гиперсфер. Я нашел ваш код настолько полезным, что я применяю его в исследовании. Тим МакДжилтон, кого я должен добавить в качестве ссылки?

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 5

Нормированный гауссовский 3d вектор равномерно распределен по сфере, см. http://mathworld.wolfram.com/SpherePointPicking.html

Например:

N = 1000
v = numpy.random.uniform(size=(3,N)) 
vn = v / numpy.sqrt(numpy.sum(v**2, 0))

Answer 6

Это будет достаточно для ваших целей?

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

Ломтик этого

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

выглядит так:

Answer 7

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()

это то, что вы хотите

Answer 8

Вы можете просто генерировать случайные точки в сферических координатах (при условии, что вы работаете в 3D): S (r, θ, φ), где r ∈ [0, R), θ ∈ [0, π], φ ∈ [0, 2π), где R - радиус вашей сферы.Это также позволит вам напрямую контролировать, сколько очков генерируется (т.е. вам не нужно сбрасывать какие-либо очки).

Чтобы компенсировать потерю плотности с радиусом, вы должны сгенерировать радиальную координату, следуя распределению по степенному закону (см. ответ dmckee для объяснения того, как это сделать).

Если вашему коду требуются (x, y, z) (т.е. декартовы) координаты, вы должны просто преобразовать случайно сгенерированные точки в сферические в декартовые координаты, как объяснено здесь .