Semicircular mean

Question

I am running MCMC sampling to determine an angular parameter resulting in several thousand samples. The angle is restricted to $[0,\pi]$ as we cannot distinguish $\alpha$ and $\alpha+\pi$. For distributions close to 0° and 180° I am running in the problem of determining an appropriate location parameter of the distribution. I attached an example histogram: Example histogram. While the histogram is clipped to $[0,\pi]$, the data could actually be mirrored at $\pi$ so that there would be a smooth distribution located around $\pi$.

My first guess was to do MLE of the $\pi$-periodic Von Mises distribution: $p(x|\mu,k)=\frac{1}{\pi I_0(x)}\exp(k\cos(2(x-\mu))$

The log likelihood yields for a set of $N$ observations $\ln p=-N\ln \pi-N\ln(I_0(k))+k\sum_i\cos(2(x_i-\mu))$

Next, we try to find the extrema:

$\frac{\partial\ln p}{\partial \mu}=2k\sum_i\sin(2(x_i-\mu))=0$

With the identities $\sin(2x)=2\sin x\cos x$, usual addition theorems and $\cos^2(x)-\sin^2(x)=\cos(2x)$ I ended up with $\mu=\frac{1}{2}\arctan\left(\frac{\sum_i\sin(2x_i)}{\sum_i\cos(2x_i)}\right)$. To project $\mu$ into the range $[0,\pi]$, I shift it via $\mu=\mu+\frac{\pi}{2}$.

For a well behaved case with mean angle far from the boundaries, this works fine: Histogram, mean in red. For the example histogram with a mean close to the boundaires this results in $\mu=89^{\circ}$: Histogram with mean. I had expected something like $179^{\circ}$ resulting in a distribution similar to this plot of the $\pi$-periodic Von Mises dist.: Von Mises with $\mu\approx179^{\circ}$

Does anyone have an idea how to resolve this?

Another problem is the determination of $k$ for a very peaky distribution when a good location parameter is available. I determine $k$ via root-finding in Python, but for high $k$ the Bessel functions just become $\infty$ or too large to handle them in scipy's implementation.

Answer 1

If anyone else runs into this problem, there is actually a simple solution already implemented in scipy: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.circmean.html. The circular mean in a defined range $[min,max]$ (here $[0,\pi]$) is calculated by shifting the common circular mean:

$\mu = \arctan2\left(\sum_i\sin(x_i),\sum_i\cos(x_i)\right)\cdot\frac{max-min}{2\pi}+min$

To estimate the concentration parameter $k$ of the Von Mises distribution, one has to solve the equation

$N\frac{I_1(k)}{I_0(k)}-\sum_i\cos(2(x_i-\mu))=0$

with the number of samples $N$. As the Bessel functions are divided by each other, one can instead use the exponentially scaled Bessel functions $e^{-x}I_n(x)$ to avoid the problem of extremely high values for the Bessel functions. Then the equation can easily be solved using a root finding algorithm.