8
$\begingroup$

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.

$\endgroup$
4
  • 1
    $\begingroup$ Did you consult the wiki article on estimation of parameters for the von Mises distribution? $\endgroup$
    – varty
    Commented Nov 21, 2011 at 4:54
  • 1
    $\begingroup$ May also want to check out Bayesian Estimation Of The von Mises Concentration Parameter $\endgroup$
    – emakalic
    Commented Nov 21, 2011 at 6:03
  • $\begingroup$ We do appreciate questions that reflect some research (which is why the downvote hover text refers to "does not show any research effort"). Questions that have good answers accessible with an obvious search or on the Wikipedia site tend to indicate no research was done or that perhaps a slightly different question was intended. $\endgroup$
    – whuber
    Commented Nov 21, 2011 at 13:48
  • $\begingroup$ Apologies, I edited the question to more specifically describe the problem. $\endgroup$ Commented Nov 21, 2011 at 16:13

3 Answers 3

10
$\begingroup$

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.

$\endgroup$
5
  • 1
    $\begingroup$ The last equation is just a rough approximation, as the paper explains, so "precisely" is the wrong wording here. $\endgroup$
    – Nick Cox
    Commented Jun 12, 2015 at 16:29
  • $\begingroup$ +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? $\endgroup$
    – Neil G
    Commented Jun 12, 2015 at 19:55
  • $\begingroup$ @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$. $\endgroup$
    – whuber
    Commented Jun 12, 2015 at 20:20
  • $\begingroup$ @whuber: Thanks! I incorrectly assumed that all of these functions were just Taylor series with a finite number of terms. $\endgroup$
    – Neil G
    Commented Jun 12, 2015 at 20:24
  • 1
    $\begingroup$ @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. $\endgroup$
    – Neil G
    Commented Jun 12, 2015 at 20:26
3
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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
$\endgroup$
1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.