Finding the orientation of a noisy ellipse
Consider the coordinates as complex numbers; let u+iv be the sum of their squares; take ½ arctan(v/u).
Later: Or rather, atan2(v,u)/2; if u<0, the naïve form above will give you the short axis instead. In Python, I'd use cmath.phase(z)/2.
Much later: I should have said in the first place that this is how to find the line through the origin that minimizes the squared distance of the sample points to the line. If the samples lie on an ellipse centred at the origin, it stands to reason that this line is close to the ellipse's long axis. But if the axis is not one of the sample lines, there may be a bias; one of these years I'll look into that.
I would assume that your points $(x_i,y_i)$, $i=1,\dots, N$ come from some normal bivariate distribution with mean $\bar x= 0$ and $\bar y = 0$. In such case you can easily find the covariance matrix: $$ \Sigma = \begin{pmatrix}a & b \\ b & c\end{pmatrix} $$ with these formulas: $$ a = \frac 1 N \sum_{i=1}^N x_i^2, \qquad b = \frac 1 N \sum_{i=1}^N x_i y_i, \qquad c = \frac 1 N \sum_{i=1}^N y_i^2. $$ Then you should compute the two eigenvectors of $\Sigma$ (or $\Sigma^{-1}$? ...to be checked) and corresponding eigenvalues. The direction you look for should be the eigenvector corresponding to the largest eigenvalue.
(see https://stats.stackexchange.com/questions/24380/how-to-get-ellipse-region-from-bivariate-normal-distributed-data)
I would build an outer product tensor and find it's eigensystem. First to build the tensor: $${\bf T} = \sum_{samples}\left(\sum_{i=0}^{11} {\bf v}_i {\bf v}_i^T\right)$$
By the spectral theorem, since it will be symmetric (why?), $\bf T$ is then ensured to have ON system with eigenvalues $\lambda_k > 0$ and normalized eigenvectors ${\bf \hat e}_k$ $${\bf T} = \sum_{i=1}^{\dim({\bf T})} \lambda_i {\bf \hat e}_i{\bf \hat e}_i^T$$ Since they are real (non-negative even), we can sort them and any orientational bias would be if they differ from each other. You can measure this in many ways, maybe variance of eigenvalues, or just difference between smallest and largest and so on.
The benefit of using this method rather than many of the others proposed here is that it easily expands to any number of dimensions. Fourier transforms on spheres are tedious, complex numbers can treat two dimensions but already at three it become a nuisance and so on.