在球形体积内采样均匀分布的随机点

Posted 2023-02-25

【中文标题】在球形体积内采样均匀分布的随机点

【英文标题】：Sampling uniformly distributed random points inside a spherical volume

【发布时间】：2011-07-21 11:10:15

【问题描述】：

我希望能够生成落在球形体积内的粒子位置的随机均匀样本。

下面的图片（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 的粒子少，才能从球体的体积中生成均匀的样本。

编辑：修复并删除了我通常要求的事实，我的意思是统一的。

【问题讨论】：

根据您的目的，查看准随机数而不是（软件）随机数也很有用

【参考方案1】：

虽然我更喜欢球体的丢弃方法，但为了完整性I offer the exact solution。

在球坐标中，利用sampling rule:

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 )

【讨论】：

我现在感觉自己像个白痴。我所要做的就是取半径的立方根吗？非常感谢您的耐心和持续参与讨论:)

@Tim：正如我在 cmets 中对 Jim 的回答所说，大多数书籍更喜欢球体的丢弃方法。即使有硬件支持，立方根也需要几个周期，并且返回笛卡尔坐标所需的三角函数也需要一些时间。另请注意，我通过统一绘制costheta 隐藏了此方法的第二个应用程序。

嗯，我明白了。为所有这些运行立方根的成本超过了所有东西的成本。我想我会两种方式都这样做，并允许我决定使用该功能。无论如何，再次感谢您的帮助！

R、u 和 costheta 是什么。我假设R 是幅度，但为什么u 是随机生成的？

@Archmede 变量u 只是一个虚拟变量，以使步骤清晰。取均匀样本的立方根的原因隐藏在我在链接问题中解释的抽样基本定理中。 R 是对的（对于半径，因为我是物理学家），costheta 是 theta 的余弦（再次因为采样定理而抛出）。

【参考方案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] 上，这是在您发布的代码中完成的

【讨论】：

如果您有空，可以看看这个问题吗？ ***.com/q/47472123/3791466

这将返回球体表面上的点，而不是体积内。

@JoanSolà 是的，它是用粗体写的

@Alleo 是的，我错过了！但它看起来像一个旁注。最好有球体积的代码，而不是表面的代码，然后是最后的注释。

【参考方案3】：

生成一组均匀分布在立方体内的点，然后丢弃与中心的距离超过所需球体半径的点。

【讨论】：

好电话。丢弃对于球体来说相当有效，并且在许多文本中被推荐为比精确采样所需的变换更快。

我想到了这个想法，但那样会损失大约 50% 的总积分。因此，要创建 16,000 个粒子，它将创建约 32,000 个。由于立方体的面积是 r^3 而球体是 4/3*pi*(r/2)^3 = 所以比率 = ~4/8 = .5

这不会导致正态分布，这是@Tim 所要求的。正态分布与均匀分布不同。

哦@Juanchopanza 我搞砸了。我的意思是制服。我现在要去解决这个问题。

请注意，这对于高维度而言效率不高，因为单位球的体积变为零（=在球内采样点的概率）。

【参考方案4】：

我同意 Alleo。我将您的 Matlab 代码翻译成 Python，它可以非常快地生成数千个点（在我的计算机中，2D 和 3D 只需要几分之一秒）。我什至已经将它用于高达 5D 的超球体。 我发现您的代码非常有用，因此我将其应用到研究中。 Tim McJilton，我应该添加谁作为参考？

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

【讨论】：

遗憾的是它不是 my matlab 代码。我想我是从mathworks.com/matlabcentral/fileexchange/… 那里得到的（2005 年写的）

@Daniel 如果您有空闲时间，请看一下这个问题吗？ ***.com/q/47472123/3791466

【参考方案5】：

这对于您的目的来说是否足够统一？

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], '.')

看起来像：

【讨论】：

这和 Jim Lewis 建议的一样，对吧？创建一个统一的立方体并扔掉任何不在球体中的东西？

@Tim McJilton：是的，我在 Jim 的回答进来时正在输入它，因为它是代码，我还是决定发布它。无论如何，在您的问题中没有任何建议表明产生更多实际需要的点会有些问题。需要详细说明吗？谢谢

我正在运行 16,000 多颗恒星的模拟，其中包括我想要统一的速度和位置位置。我希望有一种方法可以设置要保留的点数，因为我需要正好 16,000 或诸如此类。你的方法行得通，我想我可以继续逐个添加点，直到我们包含 16,000 个点。

@Tim McJilton：FWIW，在我的（非常）普通的机器中生成 1e5 3d 均匀随机点并丢弃外来者（产生大约 5.3e4 点），需要大约 35 毫秒。如果这种表现不适用于您，请告诉我们更多细节。谢谢

@Tim McJilton：现在看到dmckee 的答案我认为这是correct 这种情况的答案（球体中的确切点数）。然而，一般来说，它不会那么容易地扩展到更高的维度（如果这很重要的话）。谢谢

【参考方案6】：

范数高斯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))

【讨论】：

uniform 应该是 normal。 uniform 表示均匀分布，normal 表示高斯分布

【参考方案7】：

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

【参考方案8】：

您可以在spherical coordinates 中生成随机点（假设您在 3D 中工作）：S(r, θ, φ )，其中 r ∈ [0, R), θ ∈ [0, π ], φ ∈ [0, 2π )，其中 R 是球体的半径。这也将允许您直接控制生成的点数（即您不需要丢弃任何点）。

为了补偿半径造成的密度损失，您将根据幂律分布生成径向坐标（有关如何执行此操作的说明，请参阅 dmckee 的答案）。

如果您的代码需要 (x,y,z)（即笛卡尔）坐标，那么您只需按照 here 的说明将随机生成的球面坐标转换为笛卡尔坐标。

【讨论】：

这就是我的代码正在做的事情。当您将其转换为极坐标时，会出现问题，它在半径 R 处具有均匀数量的粒子，因为半径 r

要完成这项工作，您必须非均匀地绘制径向位置。因为体积元素是r^2 dr d\phi d(cos\theta)。这涉及提取立方根和一个反余弦，因此它往往比丢弃过程要慢。

知道了。那么丢弃点可能更容易。如果您仍想在不丢弃点的情况下这样做，则需要在power law distribution 之后生成径向坐标。不幸的是，从非平凡分布中采样并不容易，但如果您仍然感兴趣，请参阅Metropolis-Hastings 算法了解通用方法（任何其他MCMC 方法也可以）。

不要，我再说一遍不要为此使用 Metropolis！对幂律分布进行抽样很容易，但并不简单。

@AmV：它有时被称为抽样基本定律。我在此处回答的两个链接中的其他上下文中讨论它。您对 PDF 进行规范化，整合然后反转它，并使用生成的函数来转换在 [0,1) 上均匀绘制的值。在这种情况下，它归结为取立方根。

【参考方案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()

这就是你想要的

【讨论】：

你在几年后才回答。这里已经提到了两个正确的答案。你的例子是一个统一的圆圈，而不是一个统一的球体。你没有 Z 坐标...？