OK, look at it this way. Imagine that, after you generate the points randomly in the cube, and discard those outside the sphere, you then convert the remaining points into 3D polar coordinates (AKA spherical coordinates [0]). This doesn't change the distribution at all, just the numerical representation. So each point is described by three numbers, r, theta, and phi.

You're correctly pointing out that the values of r won't be uniformly distributed. There will be many more points where the value of r is close to 1 then there will be where the value of r is close to 0. This is a natural consequence of the fact that the points are uniformly distributed throughout the volume, but there's more volume near the surface than there is near the center. That's all true.

But now look at the final step. By projecting every point to the surface of the sphere, you've just overwritten every single point's r-coordinate with r=1. Any bias in the distribution of r has been discarded. This step is essentially saying "ignore r, all we care about are the values of theta and phi."

[0]https://en.wikipedia.org/wiki/Spherical_coordinate_system