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

xyz = np.zeros((n_samples, 3))
xyz[:,0] = r * np.cos(phi) * np.sin(theta)
xyz[:,1] = r * np.sin(phi) * np.sin(theta)
xyz[:,2] = r * np.cos(theta)

plt.figure(figsize=(5,5))
plt.plot(xyz[:,0], xyz[:,1], linestyle='none',
         marker=',', alpha=0.25, color='k')
plt.xlim(-5, 5)
plt.ylim(-5, 5)
plt.xlabel('$x$')
plt.ylabel('$y$')