2

这个问题中的算法告诉我们如何有效地从多维球中采样。有没有办法从多维环中同样有效地采样，即有r1<r<r2

我希望可以对该缩放函数进行不太复杂的修改 r*(gammainc(s2/2,n/2).^(1/n))./sqrt(s2)。（平庸免责声明：甚至还没有计算出原始缩放函数的代数/几何）。

原始matlab代码复制粘贴：

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 代码：

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

3 回答 3

5

这里的最后一种方法（1） 适用于任何维度的球体：

要在球体上选择一个随机点：
- 生成 N 个高斯随机变量x1,x2..xN
- 获得 x[i] 的范数

 L = Sqrt(x1*x1 + x2*x2 + .. + xn*xn)
 ux1 = x1 / L
 ux2 = x2 / L
 ...

那么向量 ux[i] 的分布在表面 S N-1上是均匀的

在环中提供均匀分布：
- 在范围内生成均匀随机

R_NPow = RandomUniform(R_Inner无）, R_Outer_

并获得半径（像这个 2D 案例）

R = R_NPow1/N

然后计算结果点坐标：

res_x1 = R * ux1
res_x2 = R * ux2
...
res_xn = R * uxn

(1) Muller, ME“关于在 -Dimensional Spheres 上均匀生成点的方法的注释”。通讯。副教授。计算。马赫 1959 年 4 月 2 日，19-20 日。

于 2017-11-24T13:43:26.417 回答
2

我实际上最终使用了应用于球体上均匀分布的点的逆 cdf 方法

像这样

def random_uniform_ring(center=np.array([0,0]),R=1,r=0,nsamples=1):
    """
    generate point uniformly distributed in a ring
    """
    nd = len(center)
    x = np.random.normal(size = (nsamples,nd))
    x = x / np.linalg.norm(x,axis=-1,keepdims=True) #generate on unit sphere
    # using the inverse cdf method
    u = np.random.uniform(size=(nsamples))
    sc = (u*(R**nd-r**nd)+r**nd)**(1/nd) #this is inverse the cdf of ring volume as a function of radius
    return x*sc[:,None]+center

去测试

import numpy as np
from scipy.special import gammainc
from matplotlib import pyplot as plt
def test1():
    fig1 = plt.figure(1)
    ax1 = fig1.gca()
    # center = np.zeros((600))
    # center = np.array([0,0])
    center = np.array([2,1])
    r = 0.5
    R = 1.
    n = 1000
    p = random_uniform_ring(center,R,r,n)
    assert p.shape[0]==n
    ax1.scatter(p[:,0],p[:,1],s=0.5)
    ax1.add_artist(plt.Circle(center,R,fill=False,color='0.5'))
    ax1.add_artist(plt.Circle(center,r,fill=False,color='0.5'))
    ax1.set_xlim(-R-0.5+center[0],R+0.5+center[0])
    ax1.set_ylim(-R-0.5+center[1],R+0.5+center[1])
    ax1.set_aspect('equal')
    plt.show()


test1()

采样_uniformly_in_a_ring

这可能相当于@Mbo 的答案，但不幸的是我真的没有时间测试。如果有人可以测试他的答案，我很乐意接受。

于 2017-11-24T14:25:23.763 回答
1

经过一些试验和错误，我能够使用 gammainc 方法完成它。它背后的数学超出了我的深度，但我基本上将 gammainc 中的系数 2 调整为 z 次方以提高均匀性。

还在 3D 中对其进行了测试，它似乎工作正常。

（这在我的待办事项列表中已经有一段时间了，谢谢你的想法！）

import numpy as np
from scipy.special import gammainc
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def sample_ring(center,r1,r2,n_points):
    nd = center.size
    x = np.random.normal(size=(n_points, nd))
    sq = np.sum(x**2,axis=1)
    z = (r2-r1)/r2
    fr = (r2-r1)*gammainc(nd/2**z,sq/2**z)**(1/nd)/np.sqrt(sq) + r1/np.sqrt(sq)
    frtiled = np.tile(fr.reshape(n_points,1),(1,nd))
    p = center + np.multiply(x,frtiled)
    return p

fig1 = plt.figure(1)
ax1 = fig1.gca()
center = np.array([0,0])
r1 = 1.5
R2 = 3
p = sample_ring(center,r1,R2,5000)
ax1.scatter(p[:,0],p[:,1],s=0.5)
ax1.add_artist(plt.Circle(center,r1,fill=False,color='0.5'))
ax1.add_artist(plt.Circle(center,R2,fill=False,color='0.5'))
ax1.set_xlim(-4,4)
ax1.set_ylim(-4,4)
ax1.set_aspect('equal')

fig3 = plt.figure(3)
ax3 = plt.gca(projection='3d')
ax3.set_aspect("equal")
theta, phi = np.mgrid[0:2*np.pi:10j, 0:np.pi:10j]
c_3d = np.array([0,0,0])
r1_3d = 0.5
x1 = c_3d[0] + r1_3d*np.cos(theta)*np.sin(phi)
y1 = c_3d[1] + r1_3d*np.sin(theta)*np.sin(phi)
z1 = c_3d[2] + r1_3d*np.cos(phi)
r2_3d = 1.4
x2 = c_3d[0] + r2_3d*np.cos(theta)*np.sin(phi)
y2 = c_3d[1] + r2_3d*np.sin(theta)*np.sin(phi)
z2 = c_3d[2] + r2_3d*np.cos(phi)
ax3.plot_wireframe(x1, y1, z1, color="r")
ax3.plot_wireframe(x2, y2, z2, color="r")
p = sample_ring(c_3d,r1_3d,r2_3d,1000)
ax3.scatter(p[:,0],p[:,1],p[:,2], c='b', marker='o')
ax3.set_xlim(-1.5, 1.5)
ax3.set_ylim(-1.5, 1.5)
ax3.set_zlim(-1.5, 1.5)

环内均匀样品

统一样品 3d 环

于 2017-11-26T01:24:06.443 回答