def sphere(x, y, z):
    return x ** 2 + y ** 2 + z ** 2 <= 1
x, y, z = np.random.uniform(-2, 2, (3, 2000000))
sphere(x, y, z).mean() * 4 ** 3, 4 / 3.0 * np.pi
Python Numerical Integration for Volume of Region