对球形体积内均匀分布的随机点进行采样

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

下图（由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 处的粒子，以便从球体体积生成均匀样本。

Edit:修正并删除了我通常要求的事实，我的意思是制服。

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

在球坐标系中，利用抽样规则:

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 )