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

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

Я был бы очень признателен за любые предложения по созданию действительно однородного образца из сферического тома в 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

Существует блестящий способ генерировать равномерно точки на сфере в 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 3

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

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

Нормированный гауссовский трехмерный вектор равномерно распределен по сфере, см. 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
def sample(center, radius, n_samples, seed=None):

    # initial values
    d = center.shape[0]

    # sample n_samples points in d dimensions from a standard normal distribution
    rng = np.random.default_rng(seed)
    samples = rng.normal(size=(n_samples, d))

    # make the samples lie on the surface of the unit hypersphere
    normalize_radii = np.linalg.norm(samples, axis=1)[:, np.newaxis]
    samples /= normalize_radii

    # make the samples lie inside the hypersphere with the correct density
    uniform_points = rng.uniform(size=n_samples)[:, np.newaxis]
    new_radii = np.power(uniform_points, 1/d)
    samples *= new_radii

    # scale the points to have the correct radius and center
    samples = samples * radius + center
    return samples

Этот код работает следующим образом.

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

Во-вторых, он вычисляет радиус каждой точки данных с помощью . (Примечание: называет евклидовой нормой, которая равна радиусу, потому что r ^ 2 = x ^ 2 + y ^ 2 в 2D и, в более общем смысле, r^2 = \sum_{i=1}^d x_i^2, что является уравнением для евклидова расстояния.) Он делит каждую точку на ее радиус, чтобы нормализовать радиус до длины 1. Это эквивалентно тому, что все точки данных лежат на поверхности единичной гиперсферы с радиусом .

В-третьих, он создает новый радиус для каждой точки данных, генерируя точек из равномерного распределения, а затем из каждой из этих точек извлекается корень-й степени. Почему мы это делаем? Рассмотрим двумерный случай, когда мы хотим равномерно выбрать точки на окружности. Если мы разделим эту окружность на концентрические кольца, то на кольце справа от центра окружности будет мало точек, а на кольце по окружности будет гораздо больше точек. В общем, чем больше радиус, тем больше точек нам нужно создать. Точнее, для достижения однородной выборки по радиусам различной длины требуется распределение вероятностей, пропорциональное D-й степени радиуса r. Мы можем выразить это как куда это количество точек радиуса от начала координат, которые дают нам равномерное распределение. Начнем с равномерного распределения точек в и может изменить приведенное выше уравнение как , то есть мы должны взять корень из равномерно распределенных радиусов.

В-четвертых, этот код масштабирует все точки на поверхности единичной гиперсферы по этому новому радиусу, равномерно распределяя эти точки по объему гиперсферы. Затем он сдвигает эти точки, чтобы получить желаемое .

Функция Производительность

Мой код намного быстрее, чем код @Daniel, который опирается на х функция, которая довольно дорога для оценки. Я оцениваю производительность здесь:

      def time_it():

    # initial values
    d = 20
    center, radius = np.full(d, 2), 3
    n_samples_list = np.logspace(3, 6, num=10).astype(int)
    runtime_my_code, runtime_daniel_code = [], []

    for n_samples in n_samples_list:

        # time my code
        start = time.perf_counter()
        sample_hypersphere_uniformly(center, radius, n_samples)
        end = time.perf_counter()
        duration = end - start
        runtime_my_code.append(duration)

        # time Daniel's code
        start = time.perf_counter()
        daniels_code(center, radius, n_samples)
        end = time.perf_counter()
        duration = end - start
        runtime_daniel_code.append(duration)

    # plot the results
    fig, ax = plt.subplots()
    ax.scatter(n_samples_list, runtime_my_code)
    ax.scatter(n_samples_list, runtime_daniel_code)
    ax.plot(n_samples_list, runtime_my_code, label='my code')
    ax.plot(n_samples_list, runtime_daniel_code, label='daniel\'s code')
    ax.set(xlabel='n_samples',
        ylabel='runtime (s)', title='runtime VS n_samples')
    plt.legend()
    fig.savefig('runtime_vs_n_samples.png')

Эта диаграмма показывает, что мой код намного эффективнее, чем у Дэниела.

Другие источники

Ознакомьтесь с этим полезным руководством по выборке из гиперсферы здесь. Это ясно объясняет, почему эта математика работает.

Answer 8

import random
R = 2

def sample_circle(center):
    a = random.random() * 2 * np.pi
    r = R * np.sqrt(random.random())
    x = center[0]+ (r * np.cos(a))
    y = center[1] + (r * np.sin(a))
    return x,y

ps = np.array([sample_circle((0,0)) for i in range(100)])

plt.plot(ps[:,0],ps[:,1],'.')
plt.xlim(-3,3)
plt.ylim(-3,3)
plt.show()

Answer 9

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 10

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

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

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