Rozwiązanie ,,naiwne'' (i błędne): losujemy dwa kąty we współrzędnych sferycznych, a następnie obliczamy współrzędne wersora wskazującego kierunek ze wzorów na transformację wkładu sferycznego do kartezjańskiego:
x= sin(θ) sin(ϕ)
y=sin(θ) cos(ϕ)
z=cos(θ)

lista={}; (* tworzymy pustą listę do której będziemy dokładać wylosowane kierunki *)
For[k=0,k<20000,k++, (* główna pętla *)
θ=RandomReal[{0,Pi}]; (* losujemy kąt theta *)
ϕ=RandomReal[{0,2 Pi}]; (* losujemy kat phi *)
n={Sin[θ]Sin[ϕ],Sin[θ]Cos[ϕ],Cos[θ]}; (* obliczamy składowe wersora *)
AppendTo[lista,n]; (* dodajemy wersor na koniec listy *)
]

MaxMemoryUsed[]

Jak widać kierunki nie są całkiem losowe, ale koncentrują się na ,,biegunach'':

ListPointPlot3D[lista,PlotRange->{{-1,1},{-1,1},{-1,1}},BoxRatios->{1,1,1},PlotStyle->PointSize[Tiny],RotationAction->"Clip"]

Rozwiązanie poprawne, polegające na losowaniu punktów wewnątrz sześcianu, odrzuceniu tych poza kulą, i przeskalowaniu wektora do długości 1:

lista={};
For[k=0,k<40000,k++,
n=RandomReal[{-1,1},3];
If[n.n≤1,AppendTo[lista,n/Sqrt[n.n]]];
]

ListPointPlot3D[lista,PlotRange->{{-1,1},{-1,1},{-1,1}},BoxRatios->{1,1,1},PlotStyle->PointSize[Tiny],RotationAction->"Clip"]

Options[ListPointPlot3D]