40

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

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

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

4

9 回答 9

46

虽然我更喜欢球体的丢弃方法，但为了完整性，我提供了精确的解决方案。

在球坐标中，利用采样规则：

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 )
于 2011-03-23T16:55:02.577 回答
18

有一种绝妙的方法可以在 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 以与 [0,1] 上的 r^n 成比例的密度分布，即在您发布的代码中完成

于 2014-05-21T13:54:22.057 回答
15

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

于 2011-03-23T16:16:46.460 回答
3

我同意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')

均匀样品

于 2017-06-27T14:27:45.517 回答
2

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

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

好像： 在此处输入图像描述

于 2011-03-23T16:20:23.147 回答
2

范数高斯 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))
于 2014-05-24T12:13:11.387 回答
0 
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()

在此处输入图像描述

于 2019-11-04T12:59:03.977 回答
-1

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

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

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

于 2011-03-23T16:37:13.460 回答
-1 

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

这就是你想要的

于 2017-10-18T10:17:42.933 回答