Estimating kappa of von Mises distribution

Question

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.

Did you consult the wiki article [on estimation of parameters for the von Mises distribution](http://en.wikipedia.org/wiki/Von_Mises_distribution#Estimation_of_parameters)? — varty, Nov 21 '11 at 04:54
May also want to check out [Bayesian Estimation Of The von Mises Concentration Parameter](http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.48.5719) — emakalic, Nov 21 '11 at 06:03
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. — whuber, Nov 21 '11 at 13:48
Apologies, I edited the question to more specifically describe the problem. — Swiss Army Man, Nov 21 '11 at 16:13

score 6 · Answer 1 · edited Jun 12 '15 at 19:30

6

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.

edited Jun 12 '15 at 19:30

whuber

281,159
54
637
1,101

answered Jun 12 '15 at 16:16

mic

3,848
3
23
38

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

score 3 · Answer 2 · answered Nov 21 '11 at 12:06

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.

answered Nov 21 '11 at 12:06

Mike Lawrence

12,691
8
40
65

score 1 · Answer 3 · edited Jun 12 '15 at 23:13

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$.

Estimating kappa of von Mises distribution

3 Answers3