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

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

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

Answer 3

在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]

BTW,看看的好链接:http://www-alg.ist.hokudai.ac.jp/~jan/randsphere.pdf

至于在球体内均匀分布,而不是对矢量进行归一化,你应该将vercor乘以一些f(r):f(r)*r的分布密度与[0,1]上的r ^ n成正比,这是在您发布的代码中完成