2

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.

Stephan Kolassa
  • 95,027
  • 13
  • 197
  • 357
Tyrion
  • 41
  • 3
  • Could you provide some details? We have no way to check whether your expectation or results are reasonable. The Bessel function issue can be handled by invoking a function to compute [exponentially scaled Bessel functions](https://math.stackexchange.com/questions/256985/scale-modified-bessel-functions-to-then-unscale-later). – whuber Nov 27 '18 at 16:13
  • Thank you for the tip with the scaled Bessel functions. Would you agree that in such a histogram the expected value of the von $\pi$ periodic Mises distribution should be close to $180^{\circ}$? This is my expectation and not fulfilled by $\mu$ the way I derived it. Additionally, I calculated a median as the angle which minimizes the sum of the distances to the other angles with the distance $\text{min}(d,\pi-d)$ with $d(\alpha,\beta)=|\alpha - \beta|$. This results in app. $179^{\circ}$ contrary to $\mu$. – Tyrion Nov 27 '18 at 17:27
  • Yes, the histogram is clear--but your method of estimating $\mu$ is not. Thus, about all we have to go on is that you have a problem, but we have almost no information to suggest how to address it. – whuber Nov 27 '18 at 17:52
  • I edited the first post, hope it can help. – Tyrion Nov 27 '18 at 18:49
  • Thank you. You should check your arctangent calculation, because the histogram suggests the numerator of its argument will be small while the denominator will be large, giving $\mu\approx 0.$ – whuber Nov 27 '18 at 19:30
  • Thank you, you are correct. I forgot to mention that I shift the mean by $\pi/2$ to get it into the range $[0,\pi]$. I edited the original question again with an example where this procedure works well. – Tyrion Nov 28 '18 at 07:38
  • Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/86318/discussion-between-tyrion-and-whuber). – Tyrion Nov 28 '18 at 08:46

1 Answers1

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.

Tyrion
  • 41
  • 3