単位球内に一様に分布する点の座標(x,y,z)は、次の３つの一様乱数

\begin{equation}
\begin{cases}
\theta & (-1\leq \theta\leq 1) \\
\phi & (0\leq \phi\leq 2\pi) \\
r & (0 \leq r \leq 1)

\end{cases}
\end{equation}

を用いて以下のように求められます。

\begin{equation}
\begin{cases}
x = r^{\frac{1}{3}}\sqrt{1-\theta^{2}}\cos\phi \\
y = r^{\frac{1}{3}}\sqrt{1-\theta^{2}}\sin\phi \\
z = r^{\frac{1}{3}}\theta

\end{cases}
\end{equation}

これをpythonで実装すると以下のようになります。

import numpy as np

xl=[]
yl=[]
zl=[]

for i in range(5000):
    theta = np.random.uniform(-1, 1)
    phi = np.random.uniform(0, 2*np.pi)
    r = np.random.uniform(0, 1)
    x=r**(1/3)*(1-theta**2)*np.cos(phi)
    y=r**(1/3)*(1-theta**2)*np.sin(phi)
    z=r**(1/3)*theta
    xl.append(x)
    yl.append(y)
    zl.append(z)

３次元の散布図は以下のコードで描くことができます。

import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
ax.scatter(xl, yl, zl, s=1)
plt.show()

次のようにして球内に分布する点の距離のリストを得ることができます。

distance=[]
for i in range(len(xl)):
    dis = np.sqrt(xl[i]**2 + yl[i]**2 + zl[i]**2)
    distance.append(dis) 

dist = 1000
distlist = list(map(lambda x: x *dist , distance))

参考