8

Is there a way to calculate an estimate of the parameter $\kappa$ from data for the von Mises distribution?

It seems very easy to do in R, http://rgm2.lab.nig.ac.jp/RGM2/func.php?rd_id=CircStats:A1inv, but python doesn't have an A1inv function to calculate the ratio of the first and zeroth order Bessel functions of the first kind.

Stephan Kolassa
  • 123,354
Swiss Army Man
  • 771
  • 1
  • 8
  • 19

3 Answers3

10

According to Banerjee et al., Clustering on the Unit Hypersphere using von Mises-Fisher Distributions (J. Mach. Learning Res. 6 (2005)), you can estimate the von Mises-Fisher parameters $\mu$ and $\kappa$ with maximum likelihood.

Let $x_i$ be the $n$ points in dimension $d$ from your sample.

Let $r = \sum_i x_i$.

Let $\overline{r} = \frac{||r||_2}{n}$ (the Euclidean distance from the barycenter to the origin). Then

$$\hat{\mu} = \frac{r}{||r||_2}$$

(the unit vector in the direction of the barycenter) and

$$\hat{\kappa} \approx \frac{\overline{r}d - \overline{r}^3}{1 - \overline{r}^2}$$

approximates the Maximum Likelihood estimate, which to be found exactly is obtained (numerically) as the solution to

$$I_{d/2}(\kappa) / I_{d/2-1}(\kappa) = \overline{r}.$$

$I_m$ is the modified Bessel function of the first kind of order $m$. The approximation can be used as the starting point for Newton-Raphson iteration--but the authors point out that for "high-dimensional" data this can be "quite slow" compared to the cost of computing just the approximate value.

whuber
  • 322,774
mic
  • 4,318
  • 1
    The last equation is just a rough approximation, as the paper explains, so "precisely" is the wrong wording here. – Nick Cox Jun 12 '15 at 16:29
  • +1, these equations agree with my notes. However, won't this be pretty fast if you use the approximate value as a starting point for a numerical optimizer? What does the dimension of the data have to do with it? $\bar r$ is a scalar after all, and typically numerical optimization runs linearly with respect to the number of bits? – Neil G Jun 12 '15 at 19:55
  • @Neil Since I added the reference to that footnote, I'll try to respond. The authors do not explain. However, having had a run-in with numerical problems in calculating Bessel functions early in my career, I suspect that accurate computation of $I_m$ when $m$ is large might not be easy, so that even a few N-R iterations could be expensive compared to the obviously fast approximation (where almost all the work is done once $r$ has been found). Plotting the Bessel ratio in Mathematica 10 (in double precision) eventually scales quadratically and becomes prohibitively slow for $d\gg 10000$. – whuber Jun 12 '15 at 20:20
  • @whuber: Thanks! I incorrectly assumed that all of these functions were just Taylor series with a finite number of terms. – Neil G Jun 12 '15 at 20:24
  • 1
    @whuber: Just to note that the question refers to the "Von Mises" distribution, which is the two-dimensional special case of the "Von Mises-Fisher" distribution. – Neil G Jun 12 '15 at 20:26
3

Check out the est.kappa() function in the CircStats package for R:

Computes the maximum likelihood estimate of kappa, the concentration parameter of a von Mises distribution, given a set of angular measurements.

Mike Lawrence
  • 13,793
2

Yes, the Von-Mises distribution family is an exponential family, so you can find the maximum likelihood estimate of its parameters as you would for any exponential family: set the expectation parameters to the average of the sufficient statistics $T(x) = x$ whose magnitude we'll call $\bar r$. You'll have to convert these parameters to your parametrization after to get $\kappa$. See @mic's answer for the equation.

Just in case you're wondering how you implement @mic's solution in Python: I would use scipy.optimize to find the root of your function: the ratio of Bessel functions minus $\bar r$.


Edit: Since this answer was posted, I've written an exponential family library for Python:

import jax.numpy as jnp
from efax import VonMisesFisherEP
mean_r = jnp.asarray([0.2, 0.4, -0.3])  # The mean observation.
x = VonMisesFisherEP(mean_r)
print(x.to_nat().kappa())  # 2.007
Neil G
  • 15,219
  • I have a follow up question that I'm hoping you might know the answer to. What concentration maximizes the entropy of the von Mises-Fisher distribution? https://stats.stackexchange.com/questions/624855/what-concentration-kappa-maximizes-the-entropy-of-the-von-mises-fisher-distri – Rylan Schaeffer Aug 25 '23 at 06:15