拒絶のない多次元リング内のサンプル均一

Question

このquestionのアルゴリズムは、多次元ボールから効率的にサンプルする方法を教えてくれます。同様に、多次元リングから同様に効率的にサンプルする方法はありますか？r1<r<r2

あまりにも複雑ではないスケーリング関数の変更が望まれます。 r*(gammainc(s2/2,n/2).^(1/n))./sqrt(s2)です。 （Mediocrity免責事項：元のスケーリング関数の代数/幾何学をまだ考えていない）。 copypasted

オリジナルMATLABコード：Daniel's answerからデモと

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

Answer 1

最後方法here^（1）は、任意の次元の球体のに適している：

球面上の任意の点を選択するには：
- Nガウスランダム変数をx1,x2..xN
を生成 - xのノルムを取得[I]

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

を次にUXベクトルの分布[i]がN-1 ^面S上

均一でありますリングに均一な分布を提供する：
は - 範囲内の一様乱数生成

R_NPow = RandomUniform(R_Inner^N, R_Outer^N）

と半径を取得する（等this 2D case）

R = R_NPow^1/N

次にポイントを得計算する座標：

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

（1）ミュラー、ME " 「次元球上で点を一様に生成する方法に関する注意」 Comm。 Assoc。計算。マッハ2、19-20、4月1959年

Answer 2

Iが実際に均一

よう球上に分布する点に適用される逆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 /=np.linalg.norm(x,axis=1)[:,np.newaxis] #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() 

これは@ Mboの回答と同じかもしれませんしかし、残念ながら私は本当にテストする時間がありません。誰かが答えをテストできるなら、私は喜んで受け入れるだろう。

Answer 3

いくつかの試行錯誤の後、私はgammaincのアプローチでそれを達成することができました。その背後にある数学は私の深みを超えていますが、私は基本的にgammaincの係数2を一様性を改善するための力zに二分しました。

また、3Dでテストしたところ正常に動作しているようです。

（これはアイデアのために、しばらくの間、私のToDoリストに感謝していた！）

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)