我正在尝试使用 numpy 在球体表面上生成随机点。我在这里查看了解释均匀分布的帖子。但是，需要有关如何仅在球体表面上生成点的想法。我有坐标（x，y，z）和每个球体的半径。
我对这个级别的数学不是很精通，并试图理解蒙特卡洛模拟。
任何帮助都感激不尽。
谢谢， 帕林
问问题
33544 次
5 回答
46
基于此页面上的最后一种方法，您可以简单地生成一个由来自三个标准正态分布的独立样本组成的向量，然后对该向量进行归一化，使其大小为 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) 或更高维单位超球面上选取均匀分布的点。
于 2015-11-28T23:11:15.127 回答
8
球体表面上的点可以用两个球坐标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)where-1 < a < 1是在均匀分布上选择的。
对于 Numpy 代码，它与在球形体积内对均匀分布的随机点进行采样相同，只是变量radius具有固定值。
于 2015-11-28T22:18:22.173 回答
3
另一种取决于硬件的方式可能会更快。
选择a, b, c三个随机数，每个在 -1 和 1 之间
计算r2 = a^2 + b^2 + c^2
如果 r2 > 1.0（=该点不在球体中）或 r2 < 0.00001（=该点离中心太近，我们将在投影到球体表面时除以零）您丢弃这些值，并选择另一组随机a, b, c
否则，您将获得随机点（相对于球心）：
ir = R / sqrt(r2)
x = a * ir
y = b * ir
z = c * ir
于 2015-11-28T22:32:11.100 回答
3
在与@Sonts 进行一些讨论后，我对答案中使用的三种方法的性能感到好奇：一种生成随机角度，一种使用正态分布坐标，另一种拒绝均匀分布的点。
这是我尝试的比较：
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 range(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 range(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
于 2015-11-29T13:47:50.610 回答
1
（编辑以反映评论的更正）
我在 2004 年调查了一些解决这个问题的恒定时间方法。
假设您在球坐标中工作，其中theta是围绕垂直轴的角度（例如经度），并且phi是从赤道升起的角度（例如纬度），然后在赤道以北的半球上获得均匀分布的随机点你来做这件事：
- 选择
theta= rand(0, 360)。 - 选择
phi= 90 * (1 - sqrt(rand(0, 1)))。
要在球体而不是半球上获得点，然后简单地否定phi50% 的时间。
出于好奇，类似的方法适用于在单位磁盘上生成均匀分布的点：
- 选择
theta= rand(0, 360)。 - 选择
radius= sqrt(rand(0, 1))。
我没有证据证明这些方法的正确性，但在过去十年左右的时间里，我已经成功地使用了它们，并且确信它们的正确性。
各种方法的一些说明（从 2004 年开始）在这里，包括在立方体表面上选择点并将它们标准化到球体上的方法的可视化。
于 2015-12-01T07:10:59.193 回答