Why do normalized gaussian points, with the same magnitudes, create a uniformly random spherical shell?

Question

I’m creating random x,y,z points on a sphere for a background I’m making. I first tried rectangular coordinates created from spherical coordinates by making angles phi (0,2pi) and theta (0,pi) uniformly random. I thought I noticed some clustering in my simulation on the z axis. I wrote some python to test it and my points are not uniformly distributed.

My next guess (with looking up anything) was to make every call to phi and theta a new random angle. This was a bad guess. It should have been obvious that the magnitudes would be different and points would not land on the sphere. It did create a cube like distribution with a neat plus sign when looking down on it.

Looking it up, I found that creating normalized rectangular coordinates from a gaussian distribution does in fact make a uniformly random sphere! I also read this is common knowledge for people who deal with probability. Reading “a multivariate standard normal distribution is spherically symmetric”. I have a very basic understanding of probability so that sentence doesn't mean much to me.

Is there an intuitive explanation for this?

import matplotlib.pyplot as plt
import numpy as np

pi = np.pi
tau = 2*np.pi
r = 1
n = 10

def randUniform(r,n):
    theta = np.random.uniform(0,1,size=n)*pi
    phi = np.random.uniform(0,1,size=n)*tau
    return r*np.array([np.sin(theta)*np.cos(phi),np.sin(theta)*np.sin(phi),np.cos(theta)]).T

def moreRandUniform(r,n):
    theta = np.random.uniform(0,1,size=(n,3))*pi
    phi = np.random.uniform(0,1,size=(n,3))*tau
    return r*np.array([np.sin(theta[:,0])*np.cos(phi[:,0]),np.sin(theta[:,1])*np.sin(phi[:,1]),np.cos(theta[:,2])]).T

def randGauss(r,n):
    xyz = np.random.normal(0, 1, size=(n,3))
    return r * xyz / np.outer(np.sqrt(np.sum(xyz**2,axis=1)),[1,1,1])


fig = plt.figure(0)
ax = fig.add_subplot(projection='3d')
points = randUniform(1,5000)
print(f'avg mag {np.mean(np.sqrt(np.sum(points**2,axis=1)))}')
ax.scatter(points[:,0], points[:,1], points[:,2], s=.1, label='random phi & theta')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.legend()


fig = plt.figure(1)
ax = fig.add_subplot(projection='3d')
points = moreRandUniform(1,10000)
print(f'avg mag {np.mean(np.sqrt(np.sum(points**2,axis=1)))}')
ax.scatter(points[:,0], points[:,1], points[:,2], s=.1, label='more random phi & theta')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.legend()

fig = plt.figure(3)
ax = fig.add_subplot(projection='3d')
points = randGauss(1,5000)
print(f'avg mag {np.mean(np.sqrt(np.sum(points**2,axis=1)))}')
ax.scatter(points[:,0], points[:,1], points[:,2], s=.1, label='gaussian xyz')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.legend()

first try

second try

third try