Plot uniform points on sphere's surface

sphere.py

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

n = 5000
r = 1

pos_neg = [1, -1]

xs, ys, zs = np.zeros(n), np.zeros(n), np.zeros(n)
for i in range(n):
    ind = np.random.randint(0,2)
    theta = np.random.rand()*2*np.pi
    cos_phi = 2.0*(np.random.rand()-0.5)
    sin_phi = np.sqrt(1-cos_phi**2)*pos_neg[ind]
    x = r*sin_phi*np.cos(theta)
    y = r*sin_phi*np.sin(theta)
    z = r*cos_phi
    xs[i], ys[i], zs[i] = x, y, z

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(xs, ys, zs, c='m', s=1)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()