I know how to pick a random point on a circle or on a sphere, but how could I pick a random point on a tilted unit circle?
By tilted, I mean that the circle is not in a xy, xz or yz plane. Instead, let $\vec u$ be the vector normal to the circle's plane.
One obvious solution would be to pick a random point on the xy unit circle, then apply a rotation matrix that would transform $\vec z$ to $\vec u$. However, this does not work because some pesky divisions by zero occur when $\vec u=\vec z$. When trying to achieve this numerically, machine precision errors also occurs when $\vec u\sim\vec z$. I could certainly separate the solution in two cases by rotating the system by $\pi/2$ when needed, but that is not practical.
Is there an elegant solution that works for any $\vec u$?
A simple solution would be to generate two unit vectors in the plane of the tilted circle. Since you already have the normal $\vec{u}=(u_x,u_y,u_z)$ you could use for instance $\vec{v} = (-u_y,u_x,0)$ and the other $\vec{w}=\vec{u} \times \vec{v} = (-u_x u_z,-u_y u_z,u_x^2+u_y^2)$ and normalise them.
Then use your favourite method to generate a point $(x,y)$ on the unit-circle and transform this to a point on the tilted circle $x \hat{u} + y \hat{v}$.
Adding the comment of Robert Israel. In the case that numerical issues would arise caused by $u_x,u_y$ being close to zero you can take a $\vec{v}=(u_z,0,-u_x)$ or $(0,u_z,-u_y)$ as an alternative. This might occur/become important if you have to generate many such points on different randomly tilted circles.
Outlined above is one of the standard ways of solving the problem, where in a computer code for particular problems that require consistent high numerical accuracy different cases will have to be considered.