<p>在与@Soonts进行了一些讨论之后,我对答案中使用的三种方法的性能感到好奇:一种是生成随机角度,一种是使用正态分布坐标,还有一种是拒绝均匀分布的点。</p>
<p>以下是我试图进行的比较:</p>
<pre><code>import numpy as np
def sample_trig(npoints):
theta = 2*np.pi*np.random.rand(npoints)
phi = np.arccos(2*np.random.rand(npoints)-1)
x = np.cos(theta) * np.sin(phi)
y = np.sin(theta) * np.sin(phi)
z = np.cos(phi)
return np.array([x,y,z])
def sample_normals(npoints):
vec = np.random.randn(3, npoints)
vec /= np.linalg.norm(vec, axis=0)
return vec
def sample_reject(npoints):
vec = np.zeros((3,npoints))
abc = 2*np.random.rand(3,npoints)-1
norms = np.linalg.norm(abc,axis=0)
mymask = norms<=1
abc = abc[:,mymask]/norms[mymask]
k = abc.shape[1]
vec[:,0:k] = abc
while k<npoints:
abc = 2*np.random.rand(3)-1
norm = np.linalg.norm(abc)
if 1e-5 <= norm <= 1:
vec[:,k] = abc/norm
k = k+1
return vec
</code></pre>
<p>然后1000分</p>
<pre><code>In [449]: timeit sample_trig(1000)
1000 loops, best of 3: 236 µs per loop
In [450]: timeit sample_normals(1000)
10000 loops, best of 3: 172 µs per loop
In [451]: timeit sample_reject(1000)
100 loops, best of 3: 13.7 ms per loop
</code></pre>
<p>注意,在基于拒绝的实现中,我首先生成了<code>npoints</code>样本并丢弃了坏的样本,并且我只使用了一个循环来生成其余的点。似乎是这样的,直接一步一步的拒绝需要更长的时间。我还删除了除法为零的检查,以便与<code>sample_normals</code>情况进行更清晰的比较。</p>
<hr/>
<p>从两种直接方法中删除矢量化会将它们放在同一个问题上:</p>
<pre><code>def sample_trig_loop(npoints):
x = np.zeros(npoints)
y = np.zeros(npoints)
z = np.zeros(npoints)
for k in xrange(npoints):
theta = 2*np.pi*np.random.rand()
phi = np.arccos(2*np.random.rand()-1)
x[k] = np.cos(theta) * np.sin(phi)
y[k] = np.sin(theta) * np.sin(phi)
z[k] = np.cos(phi)
return np.array([x,y,z])
def sample_normals_loop(npoints):
vec = np.zeros((3,npoints))
for k in xrange(npoints):
tvec = np.random.randn(3)
vec[:,k] = tvec/np.linalg.norm(tvec)
return vec
</code></pre>
<pre><code>In [464]: timeit sample_trig(1000)
1000 loops, best of 3: 236 µs per loop
In [465]: timeit sample_normals(1000)
10000 loops, best of 3: 173 µs per loop
In [466]: timeit sample_reject(1000)
100 loops, best of 3: 14 ms per loop
In [467]: timeit sample_trig_loop(1000)
100 loops, best of 3: 7.92 ms per loop
In [468]: timeit sample_normals_loop(1000)
100 loops, best of 3: 10.9 ms per loop
</code></pre>