How to average cyclic quantities? Looking on Internet, I mostly found this definition:
Given quantities on a cyclic domain D, first rescale the domain to [0;2$\pi$[, then, find $z_n$ the point on the unit circle corresponding to the $n$th value, and compute the average by:
$$z_m = \sum_{n=1}^N z_n$$
The average angle is then $\theta_m = \arg z_m$ and to obtain the average value you scale back to your original domain D.
I must say, I have problems with this definition. For simplicity, I will use oriented angles in degree for my examples (i.e. D = [0;360[). With this formula, having angles -90, 90 and 40 will give 40 as mean angle, when I would expect 13.33 as an answer (i.e. (90-90+40)/3).
For my own problems, I usually use:
$$v_m = \mathop{\rm arg\,min}_{v\in D} \sum_{n=1}^{N} d(v_n,v)^2$$
Where $d(x,y)$ is the distance in the cyclic domain, and $\{v_1, v_2, \ldots, v_n\}$ is the set of cyclic data I want to average of.
It has the advantage to work the same way whatever the domain (replace D by a non-cyclic domain and $d$ with the usual euclidean distance, and you find the usual definition of an average). However, it is expensive to compute and I don't know any exact method to do it in general.
So my question is: what is the appropriate way to deal with average of cyclic data? And do you have good pointers that explain the problem and its solutions?
 A: Like all averages, the answer depends upon the choice of metric. For a given metric $M$, the average of some angles $a_j \in [-\pi,\pi]$ for $j \in [1,N]$ is that angle $\bar{a}_M$ which minimizes the sum of squared distances $d^2_M(\bar{a}_M,a_j)$. For a weighted mean, one simply includes in the sum the weights $w_j$ (such that $\sum_j w_j = 1$). That is,
$$\bar{a}_M = \mathop{\rm arg\,min}_{x} \sum_{j=1}^{N}\, w_j\, d^2_M(x,a_j)$$
Two common choices of metric are the Frobenius and the Riemann metrics. For the Frobenius metric, a direct formula exists that corresponds to the usual notion of average bearing in circular statistics. See "Means and Averaging in the Group of Rotations", Maher Moakher, SIAM Journal on Matrix Analysis and Applications, Volume 24, Issue 1, 2002, for details.
http://lcvmwww.epfl.ch/new/publications/data/articles/63/simaxpaper.pdf
A: The problem with expecting the mean of 90°, −90°, and 40° to be (90°−90°+40°)/3 = 13.33° is that you would then expect the mean of 10° and 350° to be (10°+350°)/2 = 180°, and not 0° which is the more reasonable answer. It only gets worse when you have more than two angles (What is the mean of 340°, 350°, 360°, 10°, and 20°? What about 340°, 350°, 0°, 10°, and 20°?). Essentially, what you're doing there is equivalent to setting $z_n = e^{i\theta_n}$ and computing $$\bar z = (z_1 z_2 \cdots z_N)^{1/N},$$ and the problem is of course that it's not obvious a priori which of the $N$ possible roots of that equation is the right one, if any.
The "circular mean" definition is not so bad. In fact, it corresponds to the point which minimizes the sum of its squared distances to the points corresponding to the data, $$\bar z = \underset{\lvert z \rvert = 1}{\arg\min} \sum_{n=1}^N \lvert z - z_n \rvert^2.$$ So this is almost the same as the formula you like to use; you only have to define the "distance" between angles as the distance between the corresponding points on the unit circle. That is, $d(\theta, \phi) = \sqrt{2 - 2\cos(\theta-\phi)} = 2 \sin(\lvert\theta-\phi\rvert/2)$. This metric is close to $\lvert\theta-\phi\rvert$ when $\theta$ and $\phi$ are close, and has the advantage of being really easy to find the solution to.
A: The variables you mention are points belonging to the manifold, which is the circle. Therefore, they cannot be treated as if they belonged to Euclidean space.
I recommend the material that I have prepared on this subject and today I am sharing it on YouTube:
Circular means - Introduction to directional statistics.
There are two main types of circular mean: extrinsic and intrinsic.
Extrinsic mean is simply the mean calculated as the centroid of the points in the plane projected onto the circle.
$$
\bar{\vec{x}}=\frac{1}{N}\sum_{j=1}^N \vec{x}_j=\frac{1}{N}\sum_{j=1}^N [x_j,y_j]=\frac{1}{N}\sum_{j=1}^N [\cos{\phi_j},\sin{\phi_j}]
$$
$$
\hat{\bar{x}}=\frac{\bar{\vec{x}}}{|\bar{x}|}
$$
$$
\DeclareMathOperator{\atantwo}{atan2}
\bar{\phi}_{ex}=\atantwo(\hat{\bar{x}})
$$
It is NOT a mean calculated using the natural metric along the circle itself.
Intrinsic mean, on the other hand, does have this property. This mean can be obtained by minimizing the Fréchet function.
$$
\DeclareMathOperator*{\argmin}{argmin}
\bar{\phi}_{in}=\argmin_{\phi_0\in C} \sum_{j=1}^N (\phi_j-\phi_0)^2
$$
For discrete data, you can also analytically determine the $N$ points suspected of being the mean and then compare them using the Fréchet function.
$$
\bar{\phi}_k=\arg \sqrt[N]{\prod_{j=1}^N e^{i\phi_j} }=\bar{\phi}_0+k\frac{2\pi}{N}
$$
Where the N-th root is a N-valued function with outputs indexed with $ k\in\{1,\dots,N\} $. They are distributed evenly on the circle. And $ \bar{\phi}_0 $ is a usual mean calculated in an arbitrary range of angle values of length of $2\pi$.
If somebody dislikes the complex numbers
$$
\bar{\phi}_k=\frac{1}{N} \left(\sum_{j=1}^N \phi_j+k2\pi\right)
$$
The result is, of course, the same.
Then you have to compare the points suspected of being the mean using the Fréchet function.
$$ 
\DeclareMathOperator*{\argmin}{argmin}
\bar{\phi}_{in}=\argmin_{k\in\{1,\dots,N\}} \sum_{j=1}^N (\phi_j-\bar{\phi}_k)^2
$$
Where the search for minimum runs over $N$ discreet indices.
A: The angle is supposed to be the independent variable, not the dependent variable.  If your function is cyclic (using degrees), z(-90)=z(270) so it doesn't matter which you use.  Then the average value of the function is $$z_m=\sum_{n=1}^Nz(\theta_n)$$  You are right that if you average the angles, it matters which lap of the circle you use (-90 vs 270) because the difference of 360 gets divided by N.
A: For angles, one can adapt the iterative way of computing means to angles, that is:
given angles v[1] .. v[n]
m[1] = v[1]
m[i] = remainder( m[i-1] + remainder( v[i]-m[i-1], C)/i, C)   (i=2..n)
where remainder(x,y) gives the signed remainder on dividing x by y and C is the measure of a circle.
I suspect this gives your vm, but haven't been able to prove it.
