While I prefer the discarding method for spheres, for completeness I offer the exact solution.
In spherical coordinates, taking advantage of the sampling rule:
phi = random(0,2pi)
costheta = random(-1,1)
u = random(0,1)
theta = arccos( costheta )
r = R * cuberoot( u )
now you have a (r, theta, phi)
group which can be transformed to (x, y, z)
in the usual way
x = r * sin( theta) * cos( phi )
y = r * sin( theta) * sin( phi )
z = r * cos( theta )