# Calculating the standard deviation of a circular quantity

This is partially an algorithm question, but I think it is best asked in this stackexchange.

I need to find the standard deviation of an angle bounded on the interval $(-\pi,\pi]$. Taking the mean if such a quantity is usually done by converting the angles to points on the unit circle, taking the mean of those the x and y values of those points and converting the mean of x and y back into an angle. I don't think it's as straight-forward to make such a conversion with the standard deviation. My idea was to get the standard deviation of the points on the circle and convert the standard deviations in x and y to an angle using this classic equation for error propagation where

$f(x, y) = \arctan{\left(\frac{y}{x}\right)}$

There seems to be a scipy function that is designed to do this, scipy.stats.circstd(). The problem is my implementation gives a very different result than the scipy function, and my result gives an answer closer to what I would expect for the physical situation I am applying it to. Additionally, I do not understand what the scipy function is doing. It's source code starts on line 2773 here.

What is happening in their algorithm? What is the proper way to do this?

EDIT: This is data generated from a simulation, and I don't know the underlying distribution of the angles.

• The scipy algorithm does the following: Feb 21, 2017 at 6:30
• (whoops, pressed enter too soon) 1. Scales your sample angles into the range [0, 2*pi]. 2. Converts each sample angle into a point on the unit circle. 3. Finds the mean of those points. 4. Computes the norm of that mean point, call it r. 5. Computes sqrt(log(1 / r^2)). 6. Scales this value back to the original range of samples, [low, high]. This doesn't completely answer your question, but hopefully it helps. Feb 21, 2017 at 6:44
• Why would sqrt(log(1/(r^2))) be what they are calculating?
– Alex
Feb 21, 2017 at 6:56
• Consider that the usual standard-deviation "formula" is really just an unbiased estimator of the parameter $\sigma$ for a normal distribution that you assume your data is sampled from. I imagine that the scipy formula comes from the same notion, applied to the wrapped normal distribution. But best would be to ask the guy who seems to have wrote it! haha Feb 21, 2017 at 8:33
• In the moments section of the wrapped normal distribution article they actually seem to state the exact equation. And in the estimation of parameters section, they essentially prove it. The difficult thing to wrap your head around (no pun intended), is the wrapped normal distribution itself. Perhaps one way to interpret it is simply as a distribution of phasors $e^{j\theta_i}$ on $\mathbb{SO2}$ instead of vectors on $\mathbb{R}^n$. Feb 21, 2017 at 8:45

I'm pretty sure Scipy's algorithm does the following:

1. Uniformly scales your sample angles into the range $[0, 2\pi]$
2. Converts each sample angle into a "point" on the unit circle
3. Finds the mean of those points
4. Computes the norm of that mean point, call it r
5. Computes $\sqrt{\log(1 / r^2)}$
6. Scales this value back to the original range of the samples

All the steps seem reasonable except maybe #5. Where does that formula come from? Well first we need to ask the question, where does this formula come from? $$\hat\sigma = \sqrt{\frac{1}{n-1}\sum_{i=1}^n(x-\hat\mu)^2},\ \ \ \ \ \hat\mu=\frac{1}{n}\sum_{i=1}^nx_i$$

I write $\hat\sigma$ and $\hat\mu$ with hats because they are really estimators of the true population values of standard deviation and mean. I'll let Wikipedia do the variance proof for you, i.e. prove that $\mathbf{E}(\hat\sigma^2)=\sigma^2$.

Their results make use of some operators that we take for granted when our samples are from $\mathbb{R}$, like addition, subtraction, and a commutative multiplication. You have random angles though, which live on $\mathbb{SO2}$ and so we cannot be so care-free.

Rather than trying to guess an estimator of variance for random variables on $\mathbb{SO2}$, we can take an alternative perspective. The normal distribution on $\mathbb{R}$ is special for many reasons, but one is that its variance is exactly one of the two parameters that defines its probability density function (pdf), $$\rho(x) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{\frac{-(x-\mu)^2}{2\sigma^2}}$$

We can "wrap" this pdf around the unit circle as described here by stacking up the density at $x$ with the density at $x+2\pi k,\ \ \forall k \in \mathbb{N}$. Once we do this, we must limit our sample-space to $[0, 2\pi]$ (or any interval of length $2\pi$) or else we won't be satisfying the definition of a pdf: that its integral over the whole sample-space is 1. Make sure you shift your samples into this range before you use a wrapped pdf!

Doing this results in the wrapped normal distribution; call its probability measure $P_o$ and its density $\rho_o$. Perhaps its variance can give us some insight to what a good "circular variance" estimator would look like. To compute the variance, we can't make use of the expectation over $\mathbb{R}$ like we usually do, $$\mathbf{E}(X) \neq \int_{0}^{\infty}x\rho_o(x)dx$$ for the reason I stated a moment ago. We'll take our samples as the unit phasors (unit-magnitude complex numbers), a valid representation of $\mathbb{SO2}$ complete over the real domain $[0, 2\pi]$. $$\mathbf{E}(X) = \int_{\mathbb{SO2}}x\ dP_o = \int_{0}^{2\pi}e^{j\theta}\rho_o(\theta)d\theta$$

This is what is done here where they are able to compute the nth moment and then yank out an expression for $\sigma$ from it, that looks much like the scipy equation. In my description of the scipy algorithm, I put the word "point" in quotes because they are really treating your samples as unit phasors, not vectors. In the next section of the Wikipedia article, they find an (albeit biased) estimator of this $\sigma$ which is exactly what scipy computes.

I would say that this method is "correct" in that it isn't ad-hoc. If you find that your estimator makes more sense for your application, perhaps it is really what you want, but I would consider the scipy function to be the "right" way to measure spread of angles.

A random angle is defined as $X_2-X_1$ with $X_1$ and $X_2$ independent, uniform random variables on $(-\pi,\pi]$. Therefore:

$$Var(X_2-X_1)=Var(X_1)+(-1)^2Var(X_2)$$

$$=2Var(X_1)=2 \dfrac{1}{12}(2 \pi)^2=\dfrac23 \pi^2$$

Taking the square root:

$$\sigma=\sqrt{\dfrac23} \pi$$

• This is informative, but I don't know the underlying distribution of the angles. The angles are data generated from a simulation, and I'm looking to characterize the spread of their values. I wasn't terribly clear about this though.
– Alex
Feb 21, 2017 at 6:55