How to find a random axis or unit vector in 3D?
You need to use an equal-area projection of the sphere onto a rectangle. Such projections are widely used in cartography to draw maps of the earth that represent areas accurately.
One of the simplest such projections is the axial projection of a sphere onto the lateral surface of a cylinder, as illustrated in the following figure:
This projection is area-preserving, and was used by Archimedes to compute the surface area of a sphere.
The result is that you can pick a random point on the surface of a unit sphere using the following algorithm:
Choose a random value of $\theta$ between $0$ and $2\pi$.
Choose a random value of $z$ between $-1$ and $1$.
Compute the resulting point: $$ (x,y,z) \;=\; \left(\sqrt{1-z^2}\cos \theta,\; \sqrt{1-z^2}\sin \theta,\; z\right) $$
Another commonly used convenient method of generating a uniform random point on the sphere in $\mathbb{R}^3$ is this: Generate a standard multivariate normal random vector $(X_1, X_2, X_3)$, and then normalize it to have length 1. That is, $X_1, X_2, X_3$ are three independent standard normal random numbers. There are many well-known ways to generate normal random numbers; one of the simplest is the Box-Muller algorithm which produces two at a time.
This works because the standard multivariate normal distribution is invariant under rotation (i.e. orthogonal transformations).
This has the nice property of generalizing immediately to any number of dimensions without requiring any more thought.
You can also do this. Generate three random numbers $(a,b,c)$ in $[-1,1]$; if $a^2 + b^2 + c^2\le 1$, then normalize them. Otherwise try again and pick triplets until you have a usable triplet. The volume of the cube we pick from is 8. The volume of the unit ball is $4/3\pi$, so typically you will choose roughly two triplets to get one good random vector on the sphere.