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

Question

Я пытаюсь создать случайные точки на поверхности сферы с помощью numpy. Я рассмотрел сообщение, объясняющее единообразное распределение here. Однако нужны идеи о том, как создавать точки только на поверхности сферы. У меня есть координаты (x, y, z) и радиус каждой из этих сфер.

Я не очень хорошо разбираюсь в математике на этом уровне и стараюсь понять симуляцию Монте-Карло.

Любая помощь будет очень признательна.

Благодаря, Parin

Answer 1

точки на поверхности сферы могут быть выражены с помощью двух сферических координат, theta и phi, с 0 < theta < 2pi и 0 < phi < pi.

формула преобразования в декартовых координатах x, y, z:

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

, где r радиус сферы.

Таким образом, программа может произвольно выбирать theta и phi в своих диапазонах при равномерном распределении и генерировать декартовы координаты от него.

Но тогда очки распределяются по более длинным на полюсах сферы. Для того чтобы точки могли равномерно распределяться по поверхности сферы, phi должен быть выбран как phi = acos(a), где -1 < a < 1 выбран по равномерному распределению.

Для кода Numpy он будет таким же, как в Sampling uniformly distributed random points inside a spherical volume, за исключением того, что переменная radius имеет фиксированное значение.

Answer 2

Другой способ, который в зависимости от аппаратного обеспечения может быть намного быстрее.

Выберите a, b, c быть три случайные числа друг от -1 до 1

Вычислить r2 = a^2 + b^2 + c^2

Если г2> 1,0 (= точка не находится в области) или r2 < 0.00001 (= точка слишком близко к центру, у нас будет деление на ноль при проецировании на поверхность сферы), вы отбрасываете значения и выбираете другой набор случайных a, b, c

В противном случае у вас есть ваша случайная точка (по отношению к центру сферы):

ir = R/sqrt(r2) 
x = a * ir 
y = b * ir 
z = c * ir 

Answer 3

на основе the last approach on this page, вы можете просто сгенерировать вектор, состоящей из независимых выборок из трех стандартных нормальных распределений, то нормализовать вектор таким образом, что его величина 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 

Например:

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) 

Тот же метод обобщается выбора равномерно распределенных точек на единичной окружности (ndim=2) или на поверхностях многомерных элементарных гиперсфер.

Answer 4

После некоторого обсуждения с @Soonts мне стало интересно узнать о трех подходах, используемых в ответах: один с генерированием случайных углов, один с использованием нормально распределенных координат и один отказ от равномерно распределенных точек.

Вот моя попытка сравнения:

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 

Тогда для 1000 точек

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 

Обратите внимание, что в осуществлении отказа на основе я первым сгенерированных npoints образцы и выбросил плохие, и я только использовал цикл для генерации остальных точек. Казалось, что прямое поэтапное отторжение занимает больше времени. Я также удалил чек для деления на ноль, чтобы получить более чистое сравнение с футляром sample_normals.

---

Удаление векторизации из двух прямых методов ставит их в тот же стадион:

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 5

(отредактированный с учетом поправок от комментариев)

я исследовал несколько постоянное время подходы к t его проблемы в 2004 году

при условии, что вы работаете в сферических координатах, где theta угол вокруг вертикальной оси (например, долгота) и phi угла поднимается от экватора (например, широта), затем для получения равномерное распределение случайных точек на полусфере к северу от экватора, вы делаете это:

1. выбрать theta = рандов (0, 360).
2. выбрать phi = 90 * (1 - sqrt (rand (0, 1))).

, чтобы получить точки на сфере вместо полусферы, а затем просто отрицать phi 50% времени.

для любопытных, подобный подход имеет место для генерации равномерно распределенных точек на единицу диске:

1. выбирают theta = RAND (0, 360).
2. выбрать radius = sqrt (rand (0, 1)).

я нет доказательств для правильности этих подходов, , но я использовал их с большим количеством успеха за последнее десятилетие или около того, и я убежден в их правильности.

Некоторые иллюстрации (из 2004 г.) различных подходов составляют here, включая визуализацию подхода выбора точек на поверхности куба и нормализации их на сферу.