With[{points = 1000, samples = 500000, iterations = 20},
Nest[
With[{randoms = Join[#, RandomPoint[Sphere[], samples]]},
Normalize@Mean@randoms[[#]] & /@
Values@PositionIndex@Nearest[#, randoms]] &,
RandomPoint[Sphere[], points], iterations]] //
Graphics3D[{Sphere[{0, 0, 0}, 0.999], Red, Point@#}] &