基于the last approach on this page，只需生成一个由来自三个标准正态分布的独立样本组成的向量，然后对该向量进行规格化，使其幅值为1：

import numpy as np

def sample_spherical(npoints, ndim=3):
    vec = np.random.randn(ndim, npoints)
    vec /= np.linalg.norm(vec, axis=0)
    return vec

例如：

from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import axes3d

phi = np.linspace(0, np.pi, 20)
theta = np.linspace(0, 2 * np.pi, 40)
x = np.outer(np.sin(theta), np.cos(phi))
y = np.outer(np.sin(theta), np.sin(phi))
z = np.outer(np.cos(theta), np.ones_like(phi))

xi, yi, zi = sample_spherical(100)

fig, ax = plt.subplots(1, 1, subplot_kw={'projection':'3d', 'aspect':'equal'})
ax.plot_wireframe(x, y, z, color='k', rstride=1, cstride=1)
ax.scatter(xi, yi, zi, s=100, c='r', zorder=10)

同样的方法也适用于单位圆(ndim=2)或高维单位超球面上均匀分布的点的拾取。

球面上的点可以用两个球面坐标theta和phi表示，其中0 < theta < 2pi和0 < phi < pi。
将公式转换为笛卡尔x, y, z坐标：

x = r * cos(theta) * sin(phi)
y = r * sin(theta) * sin(phi)
z = r * cos(phi)

其中r是球体的半径。
因此，该程序可以对范围内的theta和phi进行均匀分布的随机抽样，并由其生成笛卡尔坐标。
但随后，这些点在球体的两极上分布得更加密集。为了使点在球体曲面上均匀分布，需要将phi选为phi = acos(a)，其中-1 < a < 1是按均匀分布选择的。
对于Numpy代码，它将与Sampling uniformly distributed random points inside a spherical volume中的相同，只是变量radius有一个固定值。

另一种依赖于硬件的方式可能会快得多。
选择a, b, c为三个介于-1和1之间的随机数
计算r2 = a^2 + b^2 + c^2
如果R2>0.00001(=该点不在球体内)或R2<1.0(=该点太靠近中心，我们将在投影到球体表面时被零除)，则丢弃这些值，并选取另一组随机的a, b, c
否则，您将获得随机点(相对于球体中心)：

ir = R / sqrt(r2)
x = a * ir
y = b * ir
z = c * ir

在与@Soonts进行了一些讨论之后，我对答案中使用的三种方法的性能感到好奇：一种是生成随机Angular ，一种是使用正态分布的坐标，另一种是拒绝均匀分布的点。
以下是我尝试进行的比较：

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

然后为1000分

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

请注意，在基于拒绝的实现中，我首先生成了npoints样本，并丢弃了错误的样本，我只使用了一个循环来生成其余的点。事实似乎是这样的，直接一步一步的拒绝需要更长的时间。我还删除了被零除的复选标记，以便与sample_normals案例进行更清晰的比较。
从这两种直接方法中去掉矢量化，使它们大致相同：

def sample_trig_loop(npoints):
    x = np.zeros(npoints)
    y = np.zeros(npoints)
    z = np.zeros(npoints)
    for k in range(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 range(npoints):
      tvec = np.random.randn(3)
      vec[:,k] = tvec/np.linalg.norm(tvec)
    return vec

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

• (编辑以反映评论中的更正)*

2004年，我调查了几种解决这个问题的“固定时间”方法。
假设您在球面坐标中工作，其中theta是围绕垂直轴的Angular (例如经度)，phi是从赤道向上升高的Angular (例如纬度)，则要在赤道以北的半球上获得随机点的均匀分布，请执行以下操作：
1.选择theta=rand(0,360)。
1.选择phi=90*(1-SQRT(rand(0，1)。
要在球体上而不是半球上获得点，则只需在50%的时间内取消phi。
对于好奇的人来说，在单位圆盘上生成均匀分布的点也有类似的方法：
1.选择theta=rand(0,360)。
1.选择radius=sqrt(rand(0，1))。
我没有证据证明这些方法的正确性，但在过去十年左右的时间里，我使用它们取得了很大的成功，并相信它们是正确的。
here是各种方法的一些示例(从2004年开始)，包括选择立方体表面上的点并将它们规格化到球体上的方法的可视化。

谷歌给我带来了这个古老的问题，但我在其他地方找到了一个更好的方法(参见here)
因此，对于未来的旅行者来说，最简单的方法似乎是将三个高斯/正态分布正态化如下：

coords = np.random.normal(size=[3, 1000])
distance_from_origin = np.linalg.norm(coords, ord=2, axis=0)
coords /= distance_from_origin.reshape(1,-1)
coords

在这里，我们使用numpy来初始化1000个坐标[[x * 1000],[y * 1000],[z * 1000]]的数组，每个坐标都是从以零为中心的默认高斯分布中采样的。然后，我们使用norm()函数和ord=2(也就是sqrt(x*x+y*y+z*z))计算每个坐标，并除以它到原点的距离，生成一个单位球体。
注意：对于某些行，规范可能为零！Numpy不会创建错误消息，一些值将最终为nan。添加以下检查以防止下游问题；

if (distance_from_origin==0).any():
    raise ValueError("Zero magnitude coordinates. Try again")

结果在我看来非常一致；试试看

import plotly.express as plt
x,y,z = coords
fig = plt.scatter_3d(x=x, y=y, z=z)
fig.update_traces(marker={'size': 2})

这是最快和数学上最美观的方法，在任何维度中，只要选择dim就可以了

dim = 3
radius = 1

x = np.random.normal(0,1,(100,dim)) 

z = np.linalg.norm(x, axis=1) 
z = z.reshape(-1,1).repeat(x.shape[1], axis=1) 

Points = x/z * radius * np.sqrt(dim)

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

def isotropic_unit_vectors():

    # Note: we must use arccos in the definition of theta to prevent bunching of points toward the poles

    phi = np.random.uniform(0, 2*np.pi)
    theta = np.arccos(1-2*np.random.uniform(0, 1))

    x = np.sin(theta) * np.cos(phi)
    y = np.sin(theta) * np.sin(phi)
    z = np.cos(theta)

    return [x, y, z]

# Now I shall call this function 1500 times in a while loop to plot these points

empty_array = np.empty((1500,3))

i=0
while i<1500:
    empty_array[i] = isotropic_unit_vectors()
    i+=1

x_array = empty_array[:, 0]
y_array = empty_array[:, 1]
z_array = empty_array[:, 2]

fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(x_array, y_array, z_array, s=7)
plt.show()